Application of the full spectrum inversion algorithm to simulated airborne GPS radio occultation signals

With a GPS receiver on board an airplane, the airborne radio occultation (ARO) technique provides dense lower-tropospheric soundings over target regions. Large variations in water vapor in the troposphere cause strong signal multipath, which could lead to systematic errors in RO retrievals with the geometric optics (GO) method. The spaceborne GPS RO community has successfully developed the full-spectrum inversion (FSI) technique to solve the multipath problem. This paper is the first to adapt the FSI technique to retrieve atmospheric properties (bending and refractivity) from ARO signals, where it is necessary to compensate for the receiver traveling on a non-circular trajectory inside the atmosphere, and its use is demonstrated using an end-to-end simulation system. The forward-simulated GPS L1 (1575.42 MHz) signal amplitude and phase are used to test the modified FSI algorithm. The ARO FSI method is capable of reconstructing the fine vertical structure of the moist lower troposphere in the presence of severe multipath, which otherwise leads to large retrieval errors in the GO retrieval. The sensitivity of the modified FSI-retrieved bending angle and refractivity to errors in signal amplitude and errors in the measured refractivity at the receiver is presented. Accurate bending angle retrievals can be obtained from the surface up to ∼ 250 m below the receiver at typical flight altitudes above the tropopause, above which the retrieved bending angle becomes highly sensitive to the phase measurement noise. Abrupt changes in the signal amplitude that are a challenge for receiver tracking and geometric optics bending angle retrieval techniques do not produce any systematic bias in the FSI retrievals when the SNR is high. For very low SNR, the FSI performs as expected from theoretical considerations. The 1 % in situ refractivity measurement errors at the receiver height can introduce a maximum refractivity retrieval error of 0.5 % (1 K) near the receiver, but the error decreases gradually to ∼ 0.05 % (0.1 K) near the surface. In summary, the ARO FSI successfully retrieves the fine vertical structure of the atmosphere in the presence of multipath in the lower troposphere.


Introduction
Global Positioning System (GPS) satellites transmit radio signals that undergo refractive bending and a Doppler shift due to variations in the refractive index of the Earth's atmosphere.With a GPS receiver on board an aircraft, the airborne radio occultation (ARO) receiver tracks the occulting GPS signals traversing progressively lower (or higher) atmospheric layers when the GPS satellite sets below (or rises above) the local horizon of the receiver (Healy et al., 2002;Xie et al., 2008;Haase et al., 2014).Contrary to the case for spaceborne RO receivers, the ARO receiver is located inside the atmosphere, where the non-negligible atmospheric refractive index near the receiver must be taken into account.Moreover, the RO signals from above the local horizon also need to be recorded, or otherwise accounted for, to allow the retrieval of atmospheric properties below the ARO receiver.Similar to spaceborne RO, the measurements of the ARO carrier wave phase and amplitude can be inverted to retrieve the bending angle (the cumulative atmospheric refractive bending along each ray path) as a function of the impact parameter.In a spherically symmetric atmosphere, the impact parameter, which is the product of the radius and the refractive index at the tangent point (Kursinski et al., 2000), is a con-L.Adhikari et al.: Application of the full-spectrum inversion algorithm servative quantity for each signal ray.The tangent point is the point along a ray path where the radius vector from the center of curvature is normal to the ray, and it is the closest point on the ray path to the Earth surface.The bending angle can then be converted to refractivity through the inverse Abel transformation (Fjeldbo et al., 1971).The refractivity (N ) or the refractive index (n) of the neutral atmosphere depends on the atmospheric temperature (T in Kelvin), total pressure (P in hPa), and water vapor pressure (e in hPa) (Kursinski et al., 1997(Kursinski et al., , 2000)), as follows: N = (n − 1) × 10 6 = 77.6P T + 3.73 × 10 5 e T 2 . (1) The fundamental observations made during an ARO event are the time series of phase and signal-to-noise ratio (SNR) of the RO signals.After the precise positions of the GPS and the receiver are known (e.g., Muradyan et al., 2010), the excess phase delay due to atmospheric refraction can be derived by differencing the measured signal total phase (with some initial ambiguity) with the GPS-receiver line-ofsight (LOS) distance.In this study we simulate GPS L1 signals (1575.42 MHz) for an airborne receiver at 14 km.We neglect ionospheric effects, which can be removed through the linear combination of dual frequency measurements (e.g., Vorobev and Krasil'nikova, 1994;Hajj et al., 2002).Ionospheric errors dominate spaceborne RO retrievals above 30 km (Kursinski et al., 1997); however, ARO retrievals are not possible above the height of the aircraft.In theory, the ionospheric effects are negligible for ARO retrievals in a spherically symmetric atmosphere, because the ARO retrieval requires the differencing between the RO signals originating from below (negative elevation angle) and above (positive elevation angle) the local horizon, which will cancel out most of the ionospheric errors (Xie et al., 2008) (hereafter referred to as X08).In practice, the largest part of the ionospheric error is compensated by an initial code delay in the closed-loop tracking before transitioning to open-loop tracking, and the remaining decrease in ionospheric delay (advance) over the course of the occultation is much smaller than the neutral delay (Wang et al., 2016) and other error sources.This will minimize the need for dual GPS frequency ARO measurements.The partial bending angle, defined as the difference between the bending angle measurements at the negative elevation and the bending at the positive elevation, corresponding to the same impact parameter, is then used to derive the refractivity through the inverse Abel transform.The derivative of the excess phase represents the Doppler shift of the carrier signal.The commonly used geometric optics (GO) method uses the measured Doppler and the GPSreceiver positions and velocities to retrieve the bending angle (Vorobev and Krasil'nikova, 1994).One major limitation of the GO method is its inability to account for signal interference, known as multipath, that frequently occurs in the moist lower troposphere due to sharp water vapor gradients.When multipath occurs, the signal at the receiver at a given time consists of the superposition of multiple rays, each having a unique impact parameter, and the Doppler shift derived from the signal phase no longer corresponds to a unique ray path with one impact parameter.As a result, the GO method can lead to large retrieval errors.
Various radio holographic methods have been proposed to overcome the limitations of the GO method in the spaceborne RO retrievals (e.g., Gorbunov et al., 1996;Gorbunov and Gurvich, 1998;Sokolovskiy, 2001;Gorbunov, 2002;Jensen et al., 2003Jensen et al., , 2004;;Gorbunov and Lauritsen, 2004).The fullspectrum inversion (FSI) proposed by Jensen et al. (2003) (hereafter referred to as J03) has been applied to invert spaceborne RO signals and outperforms the GO method in the presence of multipath.However, prior to this work it had not been implemented for airborne RO because of the need to address the unique characteristics of ARO occultation measurements.In particular, due to the the asymmetry of the ray path from GPS to a receiver inside the atmosphere, ARO retrieval requires additional measurement of RO signals above the local horizon.Moreover, the irregular (non-circular) flight path of the airborne platform must also be taken into account.
In this paper, the development and implementation of the FSI retrieval method for ARO measurements are presented.An end-to-end simulation system is developed to carry out a sensitivity analysis of the modified FSI algorithm, and it is compared with the GO method.The occultation geometry of ARO events used in this simulation study is based on characteristics of the flights in the PRE-Depression Investigation of Cloud-systems in the Tropics (PREDICT) field campaign over the tropical Atlantic during August-September of 2010 (Montgomery et al., 2012;Haase et al., 2014;Murphy et al., 2015).During the PREDICT campaign, the Global Navigation Satellite System (GNSS) Instrument System for Multistatic and Occultation Sensing (GISMOS) was deployed on the Gulfstream-V High-performance Instrumented Airborne Platform for Environmental Research (GV HIAPER) aircraft (Garrison et al., 2007).Haase et al. (2014) and Murphy et al. (2015) presented very promising initial results of the airborne RO observations for the upper troposphere above 7 km based on the GO retrieval.However, significant refractivity retrieval errors in the middle and lower troposphere were detected in the GO retrievals for open-loop tracking data.These retrieval errors in the lower troposphere are likely caused by the multipath problem that plagues the GO method.The FSI retrieval method presented in this paper is expected to solve the multipath problem and offer high-quality lowertropospheric ARO soundings with high vertical resolution to enable more comprehensive studies near tropical storms of the hurricane genesis process by including near-surface moisture observations.In this study, high-vertical-resolution ERA-Interim reanalysis profiles (60 vertical levels) from the European Centre for Medium Range Weather Forecasts (ECMWF) over the campaign area are used to represent the atmospheric conditions.This paper is organized as follows: Sect. 2 describes the key implementation steps for the FSI method for ARO.An end-to-end simulation system is presented in Sect.3. Section 4 presents simulated ARO observations under severe multipath conditions, where the FSI is expected to provide major improvement.The sensitivity of the FSI retrievals to the ARO measurement errors in signal amplitude and the refractivity at the receiver are explored in Sect. 5.The conclusions are summarized in Sect.6.

Theoretical derivation of FSI for airborne RO measurements
The FSI method operates directly on the measured signal along the receiver trajectory, recognizes the recorded RO signal as radio waves of different frequencies determined by the refractive index of the media through which they pass, and accounts for interference of waves with different frequencies.In contrast, the GO retrieval does not account for the possible superposition in time of multiple waves.Jensen et al. (2003) demonstrated that the derivative of the phase of the Fourier transform of the measured signal as a function of open angle distinguishes the different frequency components by using the method of stationary phase (MSP).The open angle (θ) is defined as the angle between the two radius vectors from the center of curvature to the GPS transmitter and the airborne receiver.Therefore, the time series of phase and amplitude must be resampled with respect to equally spaced open angle rather than time.For a realistic occultation with an oblate Earth and non-coplanar near-circular trajectories for the GPS and the receiver, a correction must also be applied to the observed phase to project the observations onto circular trajectories.In the case of ARO, the GPS time series is split into two parts, corresponding to the positive and negative elevation angle measurements.Each time segment can then be separately processed with the FSI.
The GPS signal composed of several narrowband subsignals resampled with respect to θ can be expressed as where A p and ϕ p are the amplitude and phase of the pth subsignal, whose Fourier transform can be expressed as where θ 1 and θ 2 are the open angles at the beginning and the end of occultation, respectively.The MSP assumes that at an instantaneous pseudo frequency in the Fourier integral in Eq. ( 3) is dominated by the contribution from a single subsignal near the stationary point (Born and Wolf, 1999).Equation (3) can then be simplified to Eq. ( 4): where B is amplitude that is approximately constant and ωq is the pseudo frequency corresponding to the subsignal at the stationary phase that satisfies ωq = dϕ q dθ | θ =θ s . (5) The derivative of the phase (ψ = ϕ q − ωq θ s ) of the Fourier transform with respect to ω gives the open angle (θ s ) that corresponds to the pseudo frequency, i.e., Using the propagation path of a ray from the transmitter to the receiver, which corresponds to the individual subsignal q mentioned above, it can be shown following Appendix A of J03 that the derivative of the phase (φ q ) of the subsignal with respect to θ (the pseudo frequency ωq ) for the ARO signal can be expressed as where k, a, n rec , R rec , and R GPS are respectively the wave number of the GPS L1 carrier signal, impact parameter, refractive index at the receiver, radius of receiver, and radius of GPS from the center of the Earth.With the ARO receiver located within the atmosphere, the difference between Eq. ( 7) above and Eq. ( 13) in J03 arises from the inclusion of the refractive index at the ARO receiver (n rec ), which is greater than 1 and is a function of receiver position.Similar to the the spaceborne RO case with receivers outside the Earth's neutral atmosphere and only measuring at negative elevation angle, the sign of the second term in Eq. ( 7) is positive at negative elevation angle, but it becomes negative at positive elevation angle when the ARO receiver and the transmitter are located on the same side of the tangent point.Note that Eq. ( 7) is derived from the calculation of the path length of a signal ray.Therefore, at positive elevation, the tangent point to the receiver path length should be subtracted from the tangent point to GPS path length to get the total path length.
The last two terms in Eq. ( 7) are frequency changes caused by the radial variations due to the non-circular trajectories of the GPS and the receiver, which will be removed through a correction described in the following section.When both trajectories are circular, the equation can be simplified to www.atmos-meas-tech.net/9/5077/2016/Atmos.Meas.Tech., 9, 5077-5087, 2016 L. Adhikari et al.: Application of the full-spectrum inversion algorithm which shows that the pseudo frequency is proportional to the impact parameter of the subsignal when the GPS and the receiver travel in a circular orbit.
After identifying the contribution from individual subsignals, these can be attributed to distinct ray paths, and subsequently the impact parameter (a) can be derived based on Eq. ( 8).Note that the impact parameter (a) of a ray is defined as (Kursinski et al., 1997) where φ GPS and φ rec are the angle between the ray tangential direction and the radius vector at the GPS and the receiver, respectively (see Fig. 1).Using Eq. ( 9) and taking n GPS = 1 at the GPS position, the angle φ GPS can be calculated as In the case of ARO, φ rec = π/2 refers to the local horizon or zero elevation angle, whereas φ rec > π/2 refers to the positive elevation angle and φ rec < π/2 refers to the negative elevation angle (see X08 for detailed description of positive and negative elevation angles).Therefore, φ rec for positive and negative elevation angles are given by Eqs. ( 10b) and ( 10c), respectively as The derivative of Eq. ( 8) gives The value of d ω is the spectrum resolution (Fowles, 1989) of the phase of the Fourier transform and is given by where θ = θ 2 − θ 1 is the difference in the open angle from the start to the end of the occultation.
The open angle can be computed using Eq. ( 6).The bending angle (α) can then be calculated using Subsequently, the refractivity profile as a function of geometric height can be derived using the inverse Abel transformation (see X08 and Healy et al., 2002, for details).

Correction for non-spherical trajectory
The FSI retrieval equations in the preceding section are derived assuming a circular orbit of the GPS and the receiver relative to the center.When the GPS and the receiver trajectories become non-circular, the radial velocity terms in Eq. ( 7) are nonzero.
In real airborne occultation measurements, the perfectly circular trajectory assumption is not valid in part because of the oblateness of the Earth, as well as any height variations of the receiver.In addition, the asymmetry of the ray path from the source to a receiver inside the atmosphere requires additional measurement of RO signals above the local horizon and correction for the irregular (non-circular) flight path of the airborne platform.To take into account the oblateness of the Earth, Syndergaard (1998) showed that the inversion of the RO data should be performed assuming local spherical symmetry tangential to the Earth's ellipsoid.
In our current approach, we account for the oblateness of the Earth by calculating the local center of curvature for each occultation event and transforming the coordinates of receiver and transmitter to that center of curvature.After the oblateness correction, a correction is applied to account for the non-spherical trajectories of the receiver and the transmitter.In the current algorithm, the correction for non-spherical trajectories has been performed by projecting the position of both the receiver and the transmitter at each epoch to a circular trajectory.Figure 1 shows the schematic diagram of the projection of the GPS signal from a non-circular receiver trajectory onto a fixed radius circular trajectory inside the atmosphere.Similarly, the method is applied to the GPS orbit to project its position (in a vacuum) onto a circular orbit.The figure shows the receiver at position P with radius R rec relative to the local center of curvature O.The measured/simulated phase at P before the projection is denoted by s 0 .The projection is done along the direction vector of the ray at P , which is determined from the occultation geometry and the ray tangent angle at or near the receiver height obtained from the CIRA86aQ_UoG model (hereafter referred as CIRA+Q) (Kirchengast et al., 1999), as follows.The first estimate of the projected position is determined using the tri-angle formed by joining the origin, O, with P and the point where the direction vector at P intersects the reference circular trajectory with radius R rec,0 .The value of R rec,0 is taken as the highest receiver altitude for the duration of the occultation and is a constant.However, the atmospheric medium between P and the projected position P 0 causes additional bending of the ray path, which can be estimated from R rec n rec sin φ rec = R rec,0 n P 0 sin φ rec,0 , ( where n P 0 is the refractivity at the projected point, P 0 . In practice, the additional phase resulting from the bending from P to P 0 is estimated as the straight line projection multiplied by the mean refractive index (n rec ), obtained from the CIRA+Q refractivity climatological model, between the original position and the projected position.The correction using a straight line approximation is sufficient at the height of the aircraft (∼ 14 km) because the vertical refractivity gradients are small at these heights and do not induce significant bending until approximately 6-7 km in the troposphere (Murphy et al., 2015).As shown in the figure, this projection from position P to P 0 leads to a change in the total phase (s), zenith angle (φ), and open angle (θ ).The changes in phase (ds rec ) and the open angel (dθ rec ) at the receiver can be calculated as Since the GPS is located outside the Earth's atmosphere, the refractive index at the GPS altitude is 1.Therefore, the corresponding change in phase at the GPS (ds GPS ) is given by The total phase after projection then becomes where s 0 is the simulated phase at P 0 before the projection.
The excess phase (s e ) can then be defined as where d is the geometric phase, i.e., the GPS-receiver LOS distance.Note that the excess phase is used to calculate the excess Doppler used in the GO retrievals.
After the projection at each sample time, the new trajectories of both GPS and receiver are circular relative to the local center of curvature, and both R rec,0 and R GPS,0 are constants, so the two radial terms in Eq. ( 7) become zero.No adjustment is made to the signal amplitude because the atmospheric attenuation at these heights over short distances is insignificant.After the adjustment, FSI is applied to the modified signal phase and the original signal amplitude with both GPS and receiver on circular trajectories.An alternative method for correcting to a spherical orbit is to apply a Doppler shift (as in J03) or a diffraction correction (Gorbunov and Lauritsen, 2004); however the phase correction in the time domain is the most straightforward to apply to the irregular flight trajectories of the research aircraft.

Estimation of the bending angle at the local horizon
In the case of airborne RO measurements, the bending angle is not a unique function of the impact parameter (a).The same a occurs twice, once each at negative and positive elevation angle relative to the local horizon where φ rec = π/2 as seen in Fig. 1 (Zuffada et al., 1999;Healy et al., 2002).To avoid the non-unique relation between a and bending angle, the GPS total phase time series is split into two parts, i.e., the positive and negative elevation angle measurements.The FSI retrieval is then applied to each part separately.The time epoch when the occulting GPS is at the local horizon relative to the ARO receiver (i.e., the separation point between positive and negative elevation angles) is estimated by raytracing through the CIRA+Q refractivity model to get the total bending angle (α) and a at each epoch.The time epoch when a reaches a maximum was selected as the separation between positive and negative elevation angle.The effect of any small time shift due to differences between the real atmosphere and the CIRA+Q climatological model will be minimized when the time series is tapered prior to the Fourier transform.
3 End-to-end simulation system for airborne RO soundings An end-to-end simulation system (Fig. 2) was developed to investigate the performance of the modified FSI algorithm for airborne RO retrievals.The simulation system consists of two major components: (i) a forward simulator and (ii) an inverse simulator, i.e., FSI retrieval.The forward simulator is used to simulate the phase and the amplitude for an ARO signal www.atmos-meas-tech.net/9/5077/2016/Atmos.Meas.Tech., 9, 5077-5087, 2016 given an atmospheric refractivity model and occultation geometry.The inverse simulator (or retrieval component) processes the simulated ARO signal phase and amplitude to retrieve the atmospheric bending angle and refractivity profiles.
3.1 Forward ARO simulator (full-spectrum forward, FSF) Two different types of forward simulators were used in the study.The first option is a ray tracer that simulates the ARO signals as geometric optics rays (e.g., Xie et al., 2008).The initial ray direction from the transmitter is iteratively perturbed until the ray is found that reaches the receiver.Forward simulation using the ray-tracing technique becomes problematic when the atmosphere has sharp refractivity gradients.If the phase fluctuates rapidly for small perturbations in initial ray direction, the technique will not converge and cannot find a ray path connecting the GPS and receiver.Often this occurs when there are multiple solutions (rays) for a given GPS-to-receiver geometry, i.e., atmospheric multipath.
The second type of forward simulator that is used to simulate the GPS signal in the presence of sharp refractivity gradients is a two-step procedure that combines an Abel integral forward model followed by a FSF simulator based on J03.The refractivity profile is input to the forward Abel integral (e.g., Fjeldbo et al., 1971;Xie et al., 2008), which calculates the bending as an integral function of radius from the surface to the aircraft height and GPS transmitter height.The resulting bending angle profile, α(a), which is geometry independent, is then related to the given source-receiver geometry θ (a) using Eqs.( 10) and ( 13).The bending angle profile as a function of the impact parameter and the corresponding θ (a) are the input atmospheric conditions to this FSF forward model to produce simulated phase and amplitude as a function of time.
The complex signal after Fourier transform in Eq. ( 4) can be expressed as a function of a as follows: where B (a) can be approximated as a constant.Each pseudo frequency corresponds to a single ray path having a unique impact parameter (a).As the pseudo frequency is proportional to the impact parameter as shown in Eq. ( 8), the derivative of the phase (ψ = ϕ q − ωq θ s ) with respect to the impact parameter yields Therefore, the phase function ψ(a) of the Fouriertransformed occultation signals can be obtained by integrating Eq. ( 7): where a 1 and a 2 are the impact parameters at the beginning and start of the occultation, and C is an unknown constant.
In this FSF forward model, the input atmospheric condition is represented by a bending angle profile.Given α(a), the θ (a) can be calculated using Eqs.( 10) and ( 13).The phase function ψ(a) and thus F (a) can be be calculated through Eqs. ( 22) and (20), respectively.
Therefore, the complex occultation signal phase and amplitude can be constructed through the inverse Fourier transform of F (a), which can then be expressed as a) da. ( The complex signal, u(θ ), after the inverse Fourier transform, has a spectral resolution in θ given by 2π k(a 2 −a 1 ) .The phase, φ(θ ), of the complex signal represents the total phase of the carrier signal with wave number k plus an unknown constant C from Eq. ( 22).The signal amplitude, A(θ ), is calculated as the sum of the conjugate components of the complex signal.The pair (A(θ ), φ(θ )) and the corresponding θ are thus evaluated from the inverse Fourier transform of the complex GPS signal in the impact parameter domain.
Subsequently, the excess phase as a function of the open angle θ, with an unknown constant, can be calculated by subtracting the LOS distance between GPS and receiver.Given the known occultation geometry for the GPS and the receiver, the relationship between the time and θ can then be applied to convert the excess phase into time space.The complex signal can therefore be generated through this full-spectrum forward operator, which eliminates the limitation of the geometric optics ray-tracing technique in the presence of multipath.

Inverse ARO simulator (FSI)
The inverse simulators in the end-to-end simulation comprise both the GO retrieval (e.g., Xie et al., 2008) and the newly developed FSI.Both the GO and FSI retrievals derive the bending angle profiles as a function of the impact parameter from the input excess Doppler as a function of time (for GO) or the combination of both the total phase and amplitude as a function of open angle (for FSI).The inverse Abel transform is then applied to retrieve refractivity from the bending angle.
In the following section, the input atmospheric refractivity and/or bending angle profiles to the forward simulator are directly compared to the output from the inverse simulators to assess the performance of the inversion technique, and to quantify the sensitivities of the ARO FSI to the potential errors presented in several key input variables such as the SNR and refractivity at the receiver.

Application of the FSI retrieval for simulated ARO measurements
To assess the performance of the ARO FSI retrieval algorithm, we used the occultation characteristics from an actual aircraft flight and the atmospheric profile of temperature and water vapor from the ERA-Interim reanalysis at the flight location and time.One specific occultation involves the GPS satellite PRN24 (the pseudo-random number identifies the satellite) and the airborne receiver during the PREDICT flight from 18:20 to 19:00 UTC on 14 September 2010 (research flight no.19).In the simulation, the radius of the curvature of the Earth was found at the occultation location; then the aircraft height was set to a constant 14 km to produce an occultation geometry with a circular orbit for the receiver at the time when PRN24 was setting.The radius of the GPS orbit was set to a constant 26 000 km above the center of curvature.The grid profiles of ERA-I temperature and water vapor mixing ratio and the calculated refractivity profiles from the ARO sounding region are shown in Fig. 3a and e, respectively.Very moist atmospheric conditions with high mixing ratio (∼ 20 g kg −1 ) are seen near the surface, and moisture decreases rapidly at higher altitude; e.g., at 10 km, the temperature reduces to around 250 K (−23.15 • C).Observations from these flights showed that the contribution of water vapor to atmospheric refractivity becomes negligible in comparison to that of temperature above a height of about 9 km (Murphy et al., 2015).Figure 3b shows the excess phase and excess Doppler obtained from the FSF forward simulator, where the initial time is referenced to zero when the open angle is 66.7 • , at ∼ 5 • elevation angle.The excess phase increases monotonically as a function of time, whereas its derivative, the excess Doppler, becomes a non-monotonic function of the time starting at ∼ 1500 s (see inset figure).Such behavior in Doppler is a strong indication of signal interference due to multipath.The multipath is further illustrated by the time series of the signal amplitude in Fig. 3c, which shows large variations around 1500 s.This signal amplitude variation is caused by superposition of multiple signals with varying frequencies.The occultation phase and amplitude time series are divided into positive and negative elevation angle parts for the FSI retrieval.The time epoch of the local horizon (zero elevation angle) is estimated by a ray-tracing simulation with the CIRA+Q refractivity model and the given occultation geometry.The phase and amplitude time series are mapped to phase and amplitude as a function of open angle (θ ), which we refer to as "sample series".The sample series require zero-padding to a power of 2 samples before the FFT can be applied.To avoid Gibbs phenomenon ringing in the transformed signal due to the finite-length sample window, a cosine taper is applied to the sample series.The zero-padding and taper reduce the amplitude at the separation of the two tapered segments (Fig. 3c); therefore small errors in determining the exact time (θ ) separating positive and negative elevation angle due to the assumed CIRA+Q refractivity model have little effect on the FFT.In real observations this tapering will reduce the influence of observations at zero elevation angle, which in any case is constrained by the in situ refractivity observations.
The bending angle retrievals from GO (blue) and FSI (red) are plotted in Fig. 3d along with the "true" bending angle profile (black) from the Abel integral forward simulation.Two distinct and important features of the ARO retrievals are shown.The first is the large error in the retrieved bending angle near zero elevation when the tangent point is near the receiver.This feature is present for both the GO and FSI methods.There is a singularity in both the GO and FSI retrievals near zero elevation angle where small errors in ray tangent angle (near π/2) can result in large bending angle errors (e.g., Eqs. 10 and 13).In the case of FSI retrievals, there can be an additional error near the receiver height due to any uncertainties in the phase correction during the projection from the non-spherical trajectory to the spherical trajectory in Eq. ( 15).However, this is resolved by knowing the refractivity at the receiver height, and hence impact parameter, from in situ measurements from the aircraft.The second feature is the large error in the GO bending angle retrieval associated with multipath in the lower troposphere.The inset in Fig. 3d shows that the GO-retrieved bending angle below impact height of 3 km (corresponding to geometric tangent point height ∼ 1 km) deviates significantly from the known (forward Abel integral) bending angle, whereas the FSI retrieval closely follows the known bending angle.The FSI is capable of resolving the sharp bending angle structure in the presence of multipath that may be caused by significant changes in moisture and/or temperature gradients near the surface.
For both the GO-and FSI-retrieved bending angle, the refractivity below the aircraft is obtained through the inverse Abel transform, by integrating the partial bending angle (defined as the difference in bending angle between the negative and positive elevation at each impact parameter) from the tangent point height up to the receiver height (Fig. 3e) (X08).The retrieved refractivity has high errors for both GO and FSI immediately below the aircraft due to higher bending angle errors discussed above.In the inverse Abel integral, the effect of bending angle errors at the receiver height propagates downward to lower levels; however the bending angle increases exponentially downwards, so the refractivity errors also decrease exponentially downward (Fig. 3f, solid lines).This is consistent with the GO simulation study in X08.
In the lower troposphere, the large refractivity errors in the GO retrieval in the lowest 1 km are due to the bending angle retrieval error in the presence of multipath.The FSI retrieval, on the other hand, successfully resolves the fine vertical structure of both bending angle and refractivity in the presence of the multipath in the moist atmosphere near the surface without introducing retrieval biases.

Sensitivity to signal amplitude and refractivity at the receiver
The accuracy of the FSI retrieval depends on the accuracy of the signal phase and amplitude measurements, the occultation geometry, and the refractivity observation at the receiver.The sensitivity of the ARO retrieval to the excess phase or Doppler error in the geometrical determination of the aircraft position and velocity has been explored in the GO retrieval system in Xie et al. (2008) and Muradyan et al. (2010), and we expect the same sensitivity for the FSI retrievals.In this section, we will quantify the sensitivity of FSI retrievals to the errors in signal amplitude (not used in the GO retrieval), amplitude-dependent phase error, and the refractivity error at the receiver.
In the ARO measurements, the amplitude is affected by the aircraft heading and attitude relative to the line of sight when using a directionally focused antenna gain pattern.Sharp amplitude variations are introduced by aircraft turns during the ARO measurements (Murphy et al., 2015).Therefore, it is important to know how sensitive the FSI retrieval is to the uncharacterized changes in signal amplitude.In addition, under low-SNR conditions, phase measurements have greater uncertainties.To test the sensitivity of the FSI retrieval algorithm to the variations in signal amplitude, we need to account for (1) the effect of sharp fluctuations in the SNR and (2) the possible phase errors that may arise under low-SNR conditions.To accomplish this task, a baseline signal phase and amplitude were first simulated using the raytracing method in a smooth refractivity model that accounts for atmospheric losses due to attenuation.Then a sinusoidal amplitude function (Eq.24) was added to simulate the sharp amplitude jumps produced by the changing aircraft direction, and finally Gaussian noise was added to the amplitude to represent the variations in the ARO amplitude measurements.
where SNR 0 is the amplitude simulated using the ray-tracing method; and k 1 , k 2 , a, and b are constants that determine the shape of the resulting amplitude; and ε is random noise.Figure 4a shows the amplitude simulated using the ray-tracing method (blue) and the modified noise-added amplitude (red).
In the simulation, the Gaussian noise power is assumed to be a constant percentage of the peak signal power during the occultation.The observed peak SNR prior to occultation is typically 200 V/V, and the observed noise floor is typically ∼ 15 V/V (Wang et al., 2016).This implies noise power to be ∼ 0.56 % of the peak power.For simplicity we used an upper bound of 1 % of the peak signal power to represent a reasonable noise variance of the simulations.Note that the residual phases with and without the amplitude noise are close to each other and appear to be on top of each other in the figure .(c) FSI-retrieved bending angle error, (d) fractional refractivity error.The radius of the Earth has been subtracted from the impact parameter in (c), where the Earth's surface is at ∼ 2.5 km.Wang et al. (2016) have shown that, at low SNR, increased phase variance results in larger errors in the unwrapped phase of the signal.Therefore, to test the impact of signal amplitude errors on the FSI retrievals, it is important to assess its impact on the measured signal phase.Wang et al. (2016) developed and tested a realistic model that relates the phase error from ARO open-loop signal processing to SNR, which we use here to estimate phase error and add to the simulated excess phase.Two different model atmospheric profiles were used in the simulations.One ERA-Interim profile (12:00 UTC, 13 September 2010 at 15 • N, 77 • W) is used to represent the true atmospheric state, and a CIRA+Q climatological model profile is used to provide the initial prediction of the excess phase and Doppler of the expected ARO signals that would be used, for example, in the open-loop signal processing.Given the ARO geometry, the excess phases are simulated based on the two model profiles using ray tracing.The phase difference between the two, i.e., the open-loop residual phase, is then generated.In the presence of measurement noise, the residual phase (ϕ) and amplitude A(t) of the received signal can be expressed as the in-phase (I ) and quadrature (Q) components as follows: where the amplitude of the noise in the two components I n and Q n is given by Eq. ( 24), and the random noise component is assumed to be independent and normally distributed with zero mean and 1 % variance.The two modified signal components in Eqs. ( 25) and ( 26) were then used to reconstruct the noise-loaded residual phase (ϕ n ) and amplitudes (A n ) as The new residual phase ϕ n (Fig. 4b in red) was then added to the CIRA+Q climatological model phase, to represent the phase of the noise-added signal (Fig. 4b, in red).
Figure 4c shows the difference between the FSI-retrieved bending angle of the noisy signal and the true bending angle, calculated by forward Abel integration of the ERA-I refractivity profile.Similarly, Fig. 4d shows the percentage error of the FSI-retrieved refractivity compared to the input refractivity profile.Both bending and refractivity errors show nearzero mean with small variations, which indicate that large variations in the amplitude measurement do not introduce systematic bias in the FSI bending and refractivity retrievals when the SNR is high.It is worth noting that in very low SNR conditions the amplitude error could potentially lead to integer cycle unwrapping errors manifesting as cycle slips, which could lead to a systematic bias in the signal phase or Doppler observation if a climatological profile is not used (e.g., Wang et al., 2016).Such biased phase or Doppler will lead to biases in both the bending and refractivity retrievals.In a simulation using noise variance of 10 % of the peak signal power (not shown), the reconstructed residual phase deviates from the original residual phase below the 3.5-4 km height range due to the large unwrapping errors.The retrieved bending angle errors increase below this height, and the refractivity errors exceed ±2 %, causing large uncertainties of the retrieved quantities below 4 km.However, such errors are a result of the inability of the open-loop tracking software receiver to retrieve a phase observation at very low signal strength and are not introduced by the retrieval process (e.g., the FSI retrieval).When using a climatological model in the open-loop tracking, the refractivity biases are reduced, resulting in biases that exceed ±2 % only below approximately 2 km (Wang, 2015).In practice, SNR is adopted as the quality control parameter determining the lower limit of reliable ARO observations.
We now assess the sensitivity to uncertainties in the in situ refractivity at the receiver.The refractivity at the receiver can be obtained from the in situ temperature, pressure, and water vapor mixing ratio measurements at the aircraft flight level.It is also one of the key parameters required in both the GO and FSI retrievals.With an ARO receiver flying at ∼ 14 km, a typical flight level during the PREDICT campaign, the water vapor contribution relative to other errors is negligible, so the refractivity at the receiver can be assumed to be only a function of temperature and pressure.To quantify the sensitivity of the FSI retrieval to the refractivity measurement error at the receiver, in addition to the phase noise, Gaussian noise of 1 % in the refractivity (∼ 2 K error in temperature) at the receiver is added to the base refractivity at the receiver.Then the noisy FSF signal phase and amplitude from the ERA-I model profile in Fig. 3 were inverted using FSI assuming the perturbed value of refractivity at the receiver.The bending angle and the retrieved refractivity were then compared to the input profiles.Fifty realizations of random Gaussian noise on the refractivity at the receiver were carried out, and the statistics of the FSI retrieval errors were then compiled.
Figure 5a and b show the absolute bending angle error and fractional refractivity error of the FSI retrieval considering the phase noise combined with the in situ refractivity error.The bending angle error is dominated by the phase error of about 0.02 • across all altitudes with a standard deviation maximum of ∼ 0.007 • near the receiver height due to the in situ refractivity error.The mean refractivity error due to phase noise in the simulation is less than 0.2 % (Fig. 5b).However, the standard deviation of the refractivity error due to the in situ refractivity error reaches a maximum of ∼ 0.5 % and decreases to ∼ 0.05 % near the surface.However, no systematic bias is introduced by such random in situ measurement errors.
The refractivity at the receiver is also used to reduce the refractivity errors in the Abel inverse.Figure 3d shows that only the top ∼ 250 m of the bending angle retrieval are noisy, due to the high sensitivity of both the GO and FSI retrieval to the measurement noise in Doppler or phase near the zero elevation.To reduce the propagation of larger refractivity errors at the receiver height to lower levels in the refractivity retrieval, the noisy bending angle observations (e.g., top 250 m below the receiver height) are replaced with a model estimate of bending angle, constrained by the in situ refractivity at the receiver extrapolated exponentially with a scale height of 7 km (Murphy et al., 2015).At high flight altitudes where water vapor is very low (e.g., above 9 km), this assumption is justified and the downward propagation of flight level refractivity errors is almost completely removed for both GO and FSI retrievals.Ultimately the retrieval accuracy of the airborne technique at flight level is limited by the Doppler errors in the aircraft position (Muradyan et al., 2010) regardless of retrieval method.

Conclusions and discussions
In this study, a FSI algorithm is developed and successfully applied to simulated airborne GNSS RO (ARO) measurements for the first time.The simulation study demonstrates the capability of the FSI method to retrieve the atmospheric vertical structure in the lower moist troposphere where frequent multipath occurs.
In the FSI retrieval process, the oblateness correction is applied to transform the original occultation geometry to a local spherical radius of curvature to fulfill the local spherical symmetric assumption.Then the non-spherical trajectories of both the ARO receiver and the GPS satellite are projected onto circular trajectories relative to local center of curvature.Additional phase correction terms as a result of projection are then added to the measured phase.Afterward, the occultation phase and amplitude time series are divided into positive and negative elevation angle sections.The separation point at the local horizon can be estimated by ray tracing through the CIRA+Q climatological model and selecting the epoch when the impact parameter is maximum at the receiver height.The FSI algorithm is then applied to the amplitude and the modified signal phase as a function of open angle for the new circular GPS-receiver trajectories for each time segment.An end-to-end simulation system is used to test the FSI retrieval using the realistic airborne occultation geometry obtained from the PREDICT field campaign.
The end-to-end simulation system was used to quantify the sensitivity of the FSI bending and refractivity retrievals to the noise in two key parameters: the signal amplitude (which induces phase errors as well as amplitude errors) and the refractivity at the receiver.The FSI retrieval showed a weak sensitivity to signal amplitude errors.Even the abrupt changes in signal amplitude due, for example, to aircraft turns do not introduce any systematic bias to the retrieval, as long as the SNR is high.The sensitivity to refractivity errors at the receiver is greater.The 1 % in situ refractivity errors at the receiver height could introduce a maximum refractivity retrieval error of 0.5 % (1 K) near the receiver, but that error decreases gradually to ∼ 0.05 % (0.1 K) near the surface.In very low SNR conditions, large amplitude variations could introduce large errors in the phase measurement and result in large uncertainties in the ARO-retrieved bending angle and refractivity profiles below ∼ 2 km.These low-SNR conditions can occur close to the surface just before the GPS is oc-culted by the Earth.However, such errors resulting from the open-loop tracking at low SNR are not introduced by the retrieval process (e.g., the FSI retrieval) and are best addressed by improvements in antenna design.

Figure 1 .
Figure 1.Schematic plot of airborne radio occultation geometry.Projection of the receiver from its original position P with radius R rec onto a new position P 0 on a circular trajectory with radius R rec,0 relative to local center of curvature O.The radius of curvature, R e , is shown as a solid black line.

Figure 2 .
Figure 2. End-to-end simulation system for airborne RO data processing.(Note: derivation of bending from refractivity in forward Abel does not need the occultation geometry information but only the receiver height and refractivity at the receiver.)

Figure 3 .
Figure 3. (a) ERA-I temperature and water vapor mixing ratio at 18:00 UTC on 13 September 2010 at 16.5 • N, 76.5 • W; (b) simulated excess phase and derived Doppler velocity (inset shows the Doppler velocity after 1500 s); (c) simulated signal amplitude; (d) bending angles retrieved using GO and FSI, and simulated bending observation with ERA-I refractivity profile (inset shows the close-up of the profile near the surface); (e) refractivity profiles from GO, FSI retrievals and the ERA-I; (f) fractional refractivity error of GO (blue), and FSI (red) compared to ERA-I.Dashed and solid lines are respectively for refractivity retrievals before and after the simple exponential atmospheric model near the receiver is applied.

Figure 4 .
Figure 4. (a) The amplitude of the received signal simulated by the ray-tracing model (blue) and the noise-added amplitude (red), (b) residual phase with (red) and without (blue) the amplitude noise.Note that the residual phases with and without the amplitude noise are close to each other and appear to be on top of each other in the figure.(c)FSI-retrieved bending angle error, (d) fractional refractivity error.The radius of the Earth has been subtracted from the impact parameter in (c), where the Earth's surface is at ∼ 2.5 km.

Figure 5 .
Figure 5. (a) Mean (black) and standard deviation of FSI-retrieved bending angle error given 1 % Gaussian refractivity error at the receiver.(b) Same as (a) but for fractional refractivity error.