## Earth

Fig. 1 Schematic drawing of the occultation geometry. The various quantities are discussed in the main text nr sin ft = p = const

where f is the angle between the ray direction and the radius vector r (Fig. 1) and n is the refractive index. Since the GPS and LEO satellites are located outside of the atmosphere where n = 1 (here we neglect the ionosphere), it follows that p = rG sin fG = rL sin fL, which equals the leveling distance between the ray and the Earth's curvature center.

When considering an atmosphere with horizontal gradients, the situation changes. The complex wave field in the inhomogeneous medium can be written in the following form: u(x) = A(x)exp(ik&(x)), where x = (x') is the coordinate vector, k = 2n/X is the wave number, X is the wavelength, A(x) is the amplitude, and &(x) is the eikonal. The geometrical optics of inhomogeneous media is based on the eikonal equation (Kravtsov and Orlov 1990):

Rays are described by the Hamilton system, where the Hamilton function follows from Eq. (2): H(p, x) = 1/2 ■ (p2 — n2(x)), where p = V& is the momentum. The associated Hamilton equations for a ray take the following form:

dx dp

where t is the trajectory parameter. Along the ray trajectories, H(p, x) = 0, therefore |dx/dr | = |p| = n(x), which allows for the conclusion that dT = di/n (since di = |dx|). Consider now polar coordinates (r, 0) and the corresponding metrics diag(1, r2). The equations for the angular component of the momentum p0 take the following form:

From the first equation we see that in the polar coordinates the angular component of momentum equals the ray impact parameter as it is defined for a spherically symmetrical medium, i.e., p = pe. The second equation generalizes Snell's law (Eq. 1) and indicates that p is only invariant in a spherically symmetrical medium (dn/dO = 0). When discussing the RO sounding of the atmosphere with horizontal gradients we must consider two impact parameters: pG = rG sin fG and pL = rL sin fL. The relation between them follows from Eq. (4):

J GPS de

In processing radio occultation data it is not possible to directly determine the two impact parameters, at the receiver and at the transmitter. Instead, an effective impact parameter can be found. It is computed from a, the time derivative of the optical path &(t) measured along the LEO observation trajectory as a function of time t, using the same formulas as for the unique impact parameter in a spherically symmetrical atmosphere. The general expression for a is as follows:

= rG cos f G + rGOG sin f g + rL cos f L - rLOL sin f l (6)

rcy rL v

Here, in the final expression for a we neglect horizontal gradients perpendicular to the occultation plane. For the effective impact parameter P we have the following implicit definition:

Using the following formula for the optical path:

where AS is the phase excess, we can obtain the effective impact parameter from orbit data and phase excess by numerically solving Eq. (7) in single ray areas, or by applying the Canonical Transform (CT) (Gorbunov 2002; Gorbunov and Lauritsen 2004), Full-Spectrum Inversion (FSI) (Jensen et al. 2003), or Phase Matching (Jensen et al. 2004) methods for unfolding multipath propagation. These techniques are based on the assumption that the effective impact parameter as defined by Eq. (7) is a unique coordinate of the ray manifold, in other words, that different rays have different p's. The (implicit) relation between pG, pL, and p is given by Eqs. (5), (6), and (7). It includes the horizontal gradient of refractivity dn/dO.

In the presence of strong horizontal gradients varying with height the effective impact parameter may no longer represent a unique coordinate of the ray manifold. One possible way of dealing with this situation would be to modify the CT method in such a way that the ray manifold is projected to a coordinate different than p. By itself, this might be straightforward except for the fact that this new unique coordinate is unknown, because it depends on the unknown horizontal gradients of refractivity. Here we will follow another prospect: instead of modifying the CT method, we will estimate bending angle errors due to strong horizontal gradients.

### 3 Radio Holographic Error Estimation

Radio holographic analysis of the wave signal can be performed both in t- and p-domains. The radio holographic (or sliding spectrum) analysis in the t-domain is widely used for the visual analysis of RO data in order to identify reflected rays, different sorts of problems etc. (Igarashi et al. 2000; Gorbunov et al. 2005). The radio holographic analysis in p-domain allows for the error estimate of retrieved bending angles (Gorbunov et al. 2005).

Here we will apply the method for the estimate of the bending angle error due to non-fully unfolded multipath structure (Fig. 2).

In the CT method, the wave field u (t) is mapped to the impact parameter representation by the Fourier Integral Operator (FIO) &. The field in the transformed space, &u(p) = A'(p) exp(ik&'(p)) can be subjected to the sliding spectral analysis (Gorbunov et al. 2005):

where Ap is the window width, & (p) is the model of the phase variation in the p-domain, and % is the bending angle variation. m(p, %) has a maximum near the true bending angle % = e(p) - e(p), where e(p) is the smooth model of the bending angle corresponding to & (p).

We used Ap = 250 m and model & (p) was computed as &'(p) smoothed with window Ap. The bending angle error is then estimated as the spectral width (cf. the schematic widths indicated in the two left panels of Fig. 2) (Gorbunov et al. 2005): Fig. 2 The principle of the radio holographic (sliding spectrum) estimate of bending angle errors. Top: the situation where the impact parameter is a unique coordinate and where the spectrum is narrow; bottom: the situation where horizontal gradients result in p not being the unique coordinate of the ray manifold. The spectral width allows for the estimate of the bending angle errors Fig. 2 The principle of the radio holographic (sliding spectrum) estimate of bending angle errors. Top: the situation where the impact parameter is a unique coordinate and where the spectrum is narrow; bottom: the situation where horizontal gradients result in p not being the unique coordinate of the ray manifold. The spectral width allows for the estimate of the bending angle errors

In non-unique ray regions the sliding spectrum will exhibit a broad structure, which eventually is mapped into a larger bending angle error than in the case where there is only one ray present for a given value of impact parameter.

### 4 Numerical Simulations

For the numerical simulations we choose an artificial occultation example, where the situation with a non-unique projection of the ray manifold to the impact parameter axis occurred. The simulated occultation was based on an ECMWF re-analysis field from February 5,1997, UTC 00:00. We simulated soundings of the same region from different azimuths. Azimuths are characterized by the angle between the occultation plane and local north direction. We changed the azimuth from 0° to 330° with a step of 30°. The resulting 12 geometric optical profiles (the simulated "truth") are shown in Fig. 3. These profiles were obtained by the geometric optical ray-tracing. Here bending angles e are shown as functions of the ray impact heights defined as p - rE, where rE is the Earth's curvature radius. The strongest effect of a  Fig. 3 True simulated bending angle profiles for sounding of the same location 12.4°N, 170.9° W from different azimuths (from 0° to 330°, with a step of 30°). Azimuths are characterized by the angle between the occultation plane and local north direction. The modeling was based on ECMWF fields February 5, 1997, UTC 00:00

Fig. 3 True simulated bending angle profiles for sounding of the same location 12.4°N, 170.9° W from different azimuths (from 0° to 330°, with a step of 30°). Azimuths are characterized by the angle between the occultation plane and local north direction. The modeling was based on ECMWF fields February 5, 1997, UTC 00:00

multi-valued bending angle profile is observed for the azimuths of 120°, 150°, and 180°. A small area of the non-unique ray manifold projection can also be noticed for the azimuth of 0°.

Figure 4 presents the CT-retrieved bending angle profiles. These profiles were retrieved from the artificial RO data simulated by wave optics (Gorbunov 2002). Because the standard CT technique relies upon the impact parameter p being a unique coordinate for the ray manifold, it is incapable of reproducing the true multivalued bending angle profiles for the azimuths of 120°, 150°, and 180°.

A small area of multi-valued bending angle profile can also be noticed for the azimuth of 0°. Single-valued profiles are retrieved instead. The azimuths, where the assumption of the uniqueness of the ray manifold parameterization with the impact parameter p is broken, are marked by the increased level of the estimated errors based on Eq. (10). In Fig. 5 we show the radio holographic running spectra in the p-domain. It is observed that for the three cases, the spectral width increases in areas of the multi-valued projection of the ray manifold, in accordance with the theoretical expectations. The error bars are, however, smaller than the difference between the  0.03 0.05 0.03 0.05 0.03 0.05 0.03 0.05 0.03 0.05 0.03 0.05

Fig. 4 True simulated bending angle profiles (GO, dashed line), retrieved simulated bending angle profiles (solid lines) and their error estimates for sounding at the same location 12.4°N, 170.9°W from different azimuths (from 0° to 330°, with a step of 30°). Azimuths are characterized by the angle between the occultation plane and local north direction. The modeling was based on ECMWF fields February 5, 1997, UTC 00:00

0.03 0.05 0.03 0.05 0.03 0.05 0.03 0.05 0.03 0.05 0.03 0.05

Fig. 4 True simulated bending angle profiles (GO, dashed line), retrieved simulated bending angle profiles (solid lines) and their error estimates for sounding at the same location 12.4°N, 170.9°W from different azimuths (from 0° to 330°, with a step of 30°). Azimuths are characterized by the angle between the occultation plane and local north direction. The modeling was based on ECMWF fields February 5, 1997, UTC 00:00

120 1 50 180 Fig. 5 Sliding spectra in the p-domain for three occultation plane directions, where the bending angle is a multi-valued function of the impact parameter (the text above the figures corresponds to the azimuth value). One observes the broadening of the spectra in the multi-valued multipath regions around 2.75-3.25 km (cf. Fig. 3)

Fig. 5 Sliding spectra in the p-domain for three occultation plane directions, where the bending angle is a multi-valued function of the impact parameter (the text above the figures corresponds to the azimuth value). One observes the broadening of the spectra in the multi-valued multipath regions around 2.75-3.25 km (cf. Fig. 3)

true multi-valued and the CT-retrieved bending angle. This is because the spike of the bending angle contains a small portion of energy, which is distributed over a stretched fragment of the ray manifold.

### 5 Conclusions

The CT/FSI methods use the fact that the impact parameter p is the unique coordinate in the ray space. This condition is universally valid for a spherically symmetrical atmosphere. In an atmosphere with horizontal gradients of refractivity this is not always the case. Strong horizontal gradients result in the variation of the impact parameters along rays rather than being invariant. This is the reason why it is generally impossible to introduce the concept of impact parameter for the atmosphere with horizontal gradients. Instead, we introduce the so-called effective impact parameter using the same definition from the Doppler frequency as for the unique impact parameter for a spherically symmetric atmosphere. In the presence of strong horizontal gradients varying with height, perturbations of the effective impact parameter may result in the non-unique projection of the ray manifold to the impact parameter axis. In this case the bending angle can be a multi-valued function of the effective impact parameter. Because the CT/FSI methods rely upon the impact parameter being the unique coordinate along the ray manifold, it follows that in this situation this method will be incapable of correctly retrieving the bending angle profile.

In the general CT formalism it may be possible to introduce a new coordinate, other than impact parameter, such that the ray manifold should have a unique projection to the corresponding axis. Still, this definition of the new coordinate depends on the specific atmospheric fields and, therefore, it cannot be specified a priori. Another opportunity is the estimation of the errors of the retrieved bending angles by means of the radio holographic (sliding spectrum) analysis in the p-domain. Numerical simulations with realistic atmospheric fields from re-analyses of ECMWF show that the error estimate increases in the presence of strong horizontal gradients resulting in multi-valued geometric optical bending angle profiles. This makes the standard CT/FSI techniques complemented with radio holographic error estimates a valuable tool for the analysis of RO data in all situations.

### References

Gorbunov ME (2002) Canonical transform method for processing radio occultation data in lower troposphere. Radio Sci 37(5), doi:10.1029/2000RS002592 Gorbunov ME, Lauritsen KB (2004) Analysis of wave fields by Fourier Integral Operators and its application for radio occultations. Radio Sci 39(4), doi:10.1029/2003RS002971 Gorbunov ME, Lauritsen KB, Rhodin A, Tomassini M, Kornblueh L (2005) Analysis of the CHAMP experimental data on radio-occultation sounding of the earth's atmosphere. Izv Atmos Ocean Phys 41(6):726-740

Gorbunov ME, Lauritsen KB, Rhodin A, Tomassini M, Kornblueh L (2006) Radio holographic filtering, error estimation, and quality control of radio occultation data. J Geophys Res 111(D10105), doi:10.1029/2005JD006427 Igarashi K, Pavelyev A, Hocke K, Pavelyev D, Kucherjavenkov IA, Matyugov S, Zakharov A, Yakovlev O (2000) Radio holographic principle for observing natural processes in the atmosphere and retrieving meteorological parameters from radio occultation data. Earth Planets Space 52(11):893-899

Jensen AS, Lohmann MS, Benzon HH, Nielsen AS (2003) Full spectrum inversion of radio occultation signals. Radio Sci 38, doi:10.1029/2002RS002763 Jensen AS, Lohmann MS, Nielsen AS, Benzon HH (2004) Geometrical optics phase matching of radio occultation signals. Radio Sci 39(3), doi:10.1029/2003RS002899 Kravtsov YA, Orlov YI (1990) Geometrical Optics of Inhomogeneous Media. Springer-Verlag,

Berlin, Heidelberg, New York Lohmann MS (2006) Dynamic error estimation for radio occultation bending angles retrieved by the full spectrum inversion technique. Radio Sci 41(RS5005), doi:10.1029/2005RS003396 Sokolovskiy SV, Rocken C, Lenschow DH, Kuo YH, Anthes RA, Schreiner WS, Hunt DC (2007) Observing the moist troposphere with radio occultation signals from COSMIC. Geophys Res Lett 34(L18802), doi:10.1029/2007GL030458

Phase Transform Algorithm for Radio Occultation Data Processing

Abstract The phase transform algorithm is described and derived in the context of GRAS instrument data processing. This algorithm allows radio occultation neutral bending angles to be recovered in a locally spherically symmetrical atmosphere in the presence of multipath. The instrument operation and carrier phase reconstruction is described for both open loop mode and closed loop mode.

### 1 Introduction

The Global Navigation Satellite Systems (GNSS) Radio Occultation Receiver for Atmospheric Sounding (GRAS) instrument (Luntama 2006; Wilson and Luntama 2007) is one of the instruments carried by the ESA/EUMETSAT METOP satellites (METOP-A, METOP-B, and METOP-C). Data from GRAS, in conjunction with other data, is used to estimate atmospheric temperature, pressure, and humidity profiles (for numerical weather prediction) and also to investigate ionospheric total electron content (TEC). An advanced processing algorithm for estimating neutral bending angles from the GRAS data is described and discussed; this algorithm is the phase transform algorithm.

### 2 The GRAS Instrument and Signal Regeneration

The GRAS instrument consists of three antennas (the GVA, the GAVA, and the GZA), three RF Conditioning Units (RFCUs), one GRAS Electronics Unit (GEU), and harness items (Loiselet et al. 2000). Each antenna is connected to its RFCU and each RFCU is connected to the GEU. The GVA and GAVA view the Earth's limb in the satellite velocity direction and satellite anti-velocity direction, respectively.

EUMETSAT, Postfach 10 05 55, 64205 Darmstadt, Germany e-mail: [email protected]

A. Steiner et al. (eds.), New Horizons in Occultation Research,

DOI 10.1007/978-3-642-00321-9.3, © Springer-Verlag Berlin Heidelberg 2009

These two antennas have 18 dual frequency receiving elements in a 3 by 6 array. The GZA is a wide beam, low gain antenna viewing the region of space above the satellite; this antenna has two receiving elements-one for each frequency. The GRAS receiver has 16 complex channels assigned to the GZA antenna, 4 complex channels assigned to the GVA antenna, and 4 complex channels assigned to the GAVA antenna. Thus the receiver can track up to 8 GPS satellites (on both L1 and L2) for navigation purposes with the GZA antenna and simultaneously observe 2 occultations with the GVA antenna (on both L1 and L2) and 2 occultations with GAVA antenna (on both L1 and L2). Each channel can be assigned to a frequency and a satellite.

The signal received at an antenna is subject to RF filtering, low noise amplification, and down conversion from RF to IF in the associated RFCU. It is then transmitted to the GEU where it is digitized, filtered, subjected to further digital down conversion, de-spread, and finally subjected to code and carrier phase tracking. The most important elements within the GEU are the A/D convertors, the DISC ASICs, the AGGA-2A ASICs, the Digital Signal Processor (DSP), the Frequency Generator (FG), and the Ultra Stable Oscillator (USO). The IF signals from the RFCUs are digitized by the A/D convertors and digitally filtered and down converted by the DISC ASICs. The AGGA-2A ASICs perform the required correlations and accumulations operations under the control of the DSP. The receiver measures, under suitable conditions, the following observable parameters for each GPS satellite tracked: (1) L1 C/A-code phase, (2) L1 P-code phase, (3) L1 carrier phase, (4) L1 signal amplitude, (5) L2 P-code phase, (6) L2 carrier phase, and (7) L2 signal amplitude. These measurements are time stamped with the GRAS Instrument Measurement Time (IMT). The carrier phases are tracked using a tracking loop to give a very narrow band output (50 Hz) except under poor phase noise conditions where the carrier phase is brought to base-band using a tracking law, based on expected behavior, which results in a wider bandwidth output (1000 Hz).

The GRAS instrument performs navigation, satellite prediction and selection, acquisition and tracking, and data handling. Navigation entails calculation of the METOP position and velocity from the data acquired by the GZA antenna. Satellite prediction and selection entails determination of which GPS satellites are to be tracked and when this is to be done for both navigation and for occultation. Satellite acquisition involves synchronizing the carrier frequency and code generators of the receiver to the received transmitted signal. C/A-code acquisition is based on a two dimensional search over defined ranges of code phase and carrier frequency using a maximum likelihood decision criteria; the search is simplified if a good navigation solution is available. When the C/A-code phase has been determined the C/A-code tracking and the L1 frequency tracking loops are started. Code tracking utilizes a normalized dot product discriminator operating on the in-phase and quadrature components of the early, late, and punctual outputs of correlators within the AGGA-2A ASIC to estimate code offset. Frequency tracking is done with a second order loop filter using the code offsets. Before starting the L1 carrier phase tracking loop a FFT based frequency estimator is employed to reduce the frequency error to an acceptable level. The carrier phase tracking is performed with an arctangent phase discriminator and a third order Costas loop. When L1 carrier tracking has been achieved the P-code on the L1 carrier and the L2 carrier can be searched; the search makes use of the hand over word from the navigation data from the tracked C/A-code. If the P-code is encrypted a code-less acquisition and tracking scheme must be used employing cross-correlation of the P-codes on the L1 and L2 channels. When the P-code on the L2 carrier is tracked, L2 carrier phase tracking can be started.

The first part of the GRAS data processing on-ground involves the reconstruction of the carrier phase received on a given channel at the antenna phase center (for that channel) from the data produced by the instrument. This involves adding back the signals removed by fixed down conversions and by the tracking loop or law-based final down conversion. This reconstruction is quite a complex and instrument specific process. The resulting signal is called the regenerated signal.

Let the set of N regenerated phase samples and their corresponding GRAS IMT for a given channel and a particular occultation be:

These IMT times are times stamped at the ADC. The regenerated phase is given by the following equation:

2n fGps(n)

+ ^constant

In Eq. (2) for the regenerated phase the terms have the following significance. The first and second terms are associated with geometric delay between the transmitter and the receiver; the first term is the normal geometrical distance and the second term is a correction for general relativistic effects due to the Earth's gravitational field. The third term is associated with the clock offset of the transmitter and receiver with respect to the inertial reference time, iREF. The fourth and fifth terms are the phases associated with the troposphere and the ionosphere, that it is desired to recover. The sixth term is an instrument characterization term. The seventh term is the cycle slip term; this term is ideally zero however cycle slips will occur if phase changes occur at a faster rate than the receiver can track during the occultation. The eighth and last term is a constant phase offset.

The measurement of the lowest part of the troposphere is difficult due to the low signal-to-noise ratio (SNR) of the GPS signal and strong atmospheric modulation. The receiver has problems in maintaining phase lock on the GPS carrier phase in a setting occultation and in acquiring the phase lock in a rising occultation. The GRAS receiver enters a "Raw Sampling" (RS) measurement mode (also called open loop mode) when the tracking of the L1 C/A-code is possible, but the phase lock on the carrier phase is not available. In the raw sampling mode the GRAS measurement data contain the output of the correlator unit of the receiver sampled at a 1000 Hz rate.

The GRAS receiver makes the switch between a phase locked tracking mode and an RS mode automatically. In a setting occultation the GRAS receiver switches normally into a single frequency tracking mode tracking the L1 code phase, carrier phase, and amplitude when the phase lock on the L2 signal is lost. At some point when the ray path tangent height is in the lower troposphere, the receiver may also lose the phase lock on the L1 carrier. The receiver software predicts this by monitoring the L1 signal and enters the RS mode if the instability of the L1 phase lock exceeds a predetermined level. This ensures the continuity of the observation data. The measurement of this particular setting occultation will end when the lock on the L1 code phase is lost.

A measurement of a rising occultation is slightly more complex as there is nothing for the receiver to track until the GPS signal becomes available from behind the limb of the Earth. The start time of a rising occultation is predicted by the receiver on-board software based on the propagated GPS and METOP orbits. The GPS orbits are calculated by using the GPS navigation message and the METOP orbit is determined by the on-board navigation solution. Normally the GRAS receiver should be able to acquire the L1 C/A-code phase very quickly after the signal becomes available. At this point phase lock on the L1 carrier is not possible and measurement of the rising occultation starts with the GRAS receiver in RS mode. The RS mode for a rising occultation continues until the receiver has acquired a stable phase lock on the L1 carrier phase. When the L1 carrier phase lock has been achieved, the receiver enters a single-frequency tracking mode. Eventually the L2 code phase signal will also become strong enough to be tracked by the receiver. After a stable lock on the L2 carrier phase has been acquired, the receiver will enter a two frequency tracking mode.

In RS mode the final down conversion is done by using a law whose parameters are given as part of the output data. This allows the carrier phase to be reconstructed in the ground processing.

There are two main differences between the carrier phase measurements provided in the RS mode and those provided when the tracking loop is locked. In raw sampling mode the sampling rate is 1000 Hz and the Navigation Signal Message (NSM) is not demodulated from the sampled data. This NSM signal is a binary signal modulated on the raw sampling mode carrier phase at 50 Hz. Where a 0 to 1 or a 1 to 0 transition occurs the carrier phase is shifted by n radians. In the ground processing the raw sampling mode data is combined with the closed loop data. The raw sampling mode data must agree with the closed loop data at the beginning of the overlapping region. To make this occur all the raw sampling mode data is rotated by 0 or n radians. The last part (or first part in the case of a rising occultation) of the single frequency tracking data is unreliable because the receiver is losing the phase lock. Thus, the overlapping closed loop data is discarded.

There are two ways to remove the NSM from the raw sampling mode data: the NSM is either removed automatically or the NSM is removed using knowledge of what it is from an alternative source. The automatic NSM removal is based on the possibility to locate one of the NSM n jumps in the region of the raw sampling mode data near to the closed loop data. Once the location of one carrier phase jump by n is found the positions of all other potential transitions are determined. The raw sampling mode carrier phase may therefore be averaged over regions (e.g., 20 samples) where the NSM is constant to reduce noise. In this way the NSM can be estimated and removed from the raw sampling mode carrier phase by rotating the phase by 0 or n according to the value of the NSM at that sample. Automatic removal will become more difficult the further one goes away from the closed loop data and will fail when the phase noise associated with the atmosphere approaches the n jumps associated with NSM bit transitions. Raw sampling mode data from which the NSM cannot be removed must be discarded. The knowledge based removal uses an NSM bit stream recorded by an independent GPS receiver (e.g., network of ground based receivers). The complete NSM sequence from the occulting satellite is then correlated in time with the received RS data and phase jumps matching the NSM bit transitions are removed. Knowledge based removal is more reliable than the automatic removal scheme, but requires continuous recording of the NSM sequences of all GPS satellites.

### 3 The Phase Transform Algorithm

The classical geometric optics algorithm is applicable provided no multi-path propagation occurs in the atmosphere. If multi-path propagation occurs, a more sophisticated algorithm is required to recover the neutral bending angle as a function of impact parameter, namely the phase transform algorithm. The phase transform algorithm was originally derived by A. S. Jensen (Jensen et al. 2003, 2004, 2006), in a slightly different form. The derivation given here is explicit.

To perform the phase transform algorithm the regenerated phase must be first corrected for (a) instrument characteristics, (b) general relativity, and (c) clock drifts. After these corrections have been performed a complex signal with unit amplitude is constructed from the resulting regenerated phase ^REG-B using the following equation:

Next the phase transform is performed:

u(p) = « (rRMT) ^REG-B (4MT) exp (-i ß (p, tRMT)) drRMT, (4) Jo where m is a weighting function introduced to prevent ringing and ¡3(p, t) is the reference phase function for impact parameter p. The bending angle is given by the following equation:

The bending angles from the L1 and L2 frequencies can be smoothed and then combined to give the neutral bending angle.

fL1 fL2

The reference phase function is given by the following equations where the geometry is referred to the time in the phase transform integral (IMT time or reference time).

ß(P, T) = +T I V rLEc(T) - P2 + V rGps(T) - p2 + pr(T)

- p arctan

- p arctan

The phase transform algorithm is essentially a matched filter that selects from the signal received by the GRAS receiver that portion of the signal associated with a specific impact parameter. The phase transform algorithm thus allows the bending angle to be recovered in the presence of multi-path provided the atmosphere is (locally) spherically symmetric.

4 Derivation of the Phase Transform Algorithm

The phase transform algorithm is best understood by reviewing its derivation; this is given below. Consider the Fig. 1 (Born and Wolf 1983). From Fig. 1 it is seen that:

For a spherically symmetric atmosphere the formula of Bouguer, Eq. (10), is valid p = n(r )r sin(^). (10)

Therefore

The phase along the ray path from a GPS satellite to the GRAS receiver on a METOP satellite is given by:

p = —J n(s )ds = —J n(r )J 1 + r 2 ( —) dr. (12)

And thus

2n fr

2n f2 n2(r)rdr

2n fri n2(r)rdr 2n fr2 n2(r)rdr

Above r0 is the point of closest approach. Now consider one of these integrals:

J rQ

Using the identity:

1 dn

2 dr n dr n dr and the equation:

The integral becomes:

2np rrA

X Jrn dr

X Jro jn2(r)r2 - p2 dr Now from Fig. 1 and from the formula of Bouguer it is seen that:

Therefore

X Jro V

And therefore the phase is given by

-2 r ynw-pdilnnu

Jro dr J

Where, with reference to Fig. 2:

^gps p2

Now consider the integral:

rTiMr u(p) = „ (zRXT) ^reg-b (r R|T) exp (-ifi (p, tRMT)) d tRMt Jo

Thus

/»Timt u(p) — « (tRXt) Areg-b (tRXT) exp (i <p (p, tRMt) - i ß (p, tRXt)) d tRMt Jo

The integral can be evaluated (approximately) by the method of stationary phase to be:

^« (rRXT0 Areg-b (tRXt0 exp (iv (p, - i ß (p, tRXTi)) Fig. 2 Illustrating the occultation plane geometry

where

Now P(p, t) can be chosen arbitrarily; if P(p, t) is chosen to be:

P(p, t) = t(^eo(t) - p2 + 7rGps(T) - p2 + p a(p,T^j , (31)

where

Then with reference to Eq. (24):

Therefore we have:

— (p(p, T1) - P(p, T1)) = —- — Vn2(r)r2 - p2---dr

A 2n

Equation (34) establishes the desired result. The main steps in the derivation are: (i) expressing the corrected phase between the transmitter and the receiver in the form given in Eq. (24), (ii) solving the Eq. (28) for u(p) by the method of stationary phase, (iii) noticing that if p(p, t) is chosen to be given by Eq. (31) then Eq. (33) will give arg(u(p)) and (iv) observing that the derivative of arg(u(p)) with respect to p is proportional to the bending angle as seen from Eq. (34). It is seen from the derivation that all the information about the bending angle is contained in the (corrected) regenerated phase measured by the GRAS receiver. The amplitude of the signal is not used.

The phase transform algorithm is valid for any spacecraft orbits provided an occultation plane exists (at least approximately). It is an exact algorithm and not an approximation. The Full Spectrum Inversion (FSI) algorithm (Jensen et al. 2003) is an approximation of the phase transform algorithm which is applicable assuming a spherical Earth and circular spacecraft orbits in the same plane. In this situation Eq. (4) can be approximated by a Fourier transform.

### References

Born M, Wolf E (1983) Principles of Optics, 6th edn. Cambridge University Press, England Jensen AS, Lohmann MS, Benzon HH, Nielsen AS (2003) Full spectrum inversion of radio occultation signals. Radio Science 38(3), doi:10.1029/2002RS002763 Jensen AS, Lohmann MS, Nielsen AS, Benzon HH (2004) Geometrical optics phase matching of radio occultation signals. Radio Science 39(RS3009), doi:10.1029/2003RS002899 Jensen AS, Benzon HH, Lohmann MS, Nielsen AS (2006) Processing radio occultation data by full spectrum inversion techniques: An overview and recent developments. In: Foelsche U, Kirchengast G, Steiner AK (eds) Atmosphere and Climate: Studies by Occultation Methods, Springer-Verlag, Berlin Heidelberg New York, pp 95-112 Loiselet M, Stricker N, Menard Y, Luntama LP (2000) GRAS - Metop's GPS-Based Atmospheric

Sounder. ESA Bulletin 102, pp 38-44 Luntama JP (2006) The operational EPS GRAS measurement system. In: Foelsche U, Kirchengast G, Steiner AK (eds) Atmosphere and Climate: Studies by Occultation Methods, Springer, Berlin Heidelberg, pp 147-156, doi:10.1007/3-540-34121-8 Wilson JJW, Luntama JP (2007) Advanced processing algorithms for GRAS instrument data. In: Geoscience and Remote Sensing Symposium, IGARSS 2007, IEEE International, Barcelona, Spain, pp 1049-1054, doi:10.1109/IGARSS.2007.4422981

Using Airborne GNSS Receivers to Detect Atmospheric Turbulence

L.B. Cornman, A. Weekley, R.K. Goodrich, and R. Frehlich

Abstract A methodology for estimating turbulence intensity from GNSS-aircraft occultations is presented. The theoretical underpinnings are from standard weak-scattering theory for electromagnetic wave propagation in random media. These techniques are modified to deal with a transmitter and receiver moving relative to each other. A simulation method is then used to evaluate the sensitivity of the intensity estimates to two other unknown parameters, the turbulence length scale and the distance of the turbulence from the receiver. It is shown that the estimation is highly sensitive to the latter and relatively insensitive to the former. An iterative technique is presented that uses estimates of the distance parameters to improve the intensity estimation. It is shown that given the assumptions in the problem, the iterative technique provides relatively accurate estimates of the turbulence intensity parameter.

### 1 Introduction

There is a long and distinguished history in the study of electromagnetic (EM) wave propagation through random media (cf. Tatarskii 1971; Yeh and Liu 1982; Ishimaru 1997). These works have provided a firm theoretical foundation for estimating statistical properties of the neutral atmosphere and ionosphere via the statistical properties of the received EM signals. That is, the characteristics of the turbulent atmosphere can be deduced from, for example, correlation and/or spectral analysis of the phase and/or amplitude of the received signal. In the past, the bulk of the experimental analysis in this area has been performed with ground-based transmitters and receivers (e.g., radars and lidars), as well as with ground-based receivers and space-based transmitters. With the advent of Global Navigation Satellite System (GNSS) constellations, a new avenue has become available to investigate the turbulence properties of the Earth's atmosphere. The deployment of an ever-increasing number of Low Earth Orbiting (LEO) satellites—with high-quality, high-sampling rate receivers—provides a very valuable new source of turbulence measurements:

National Center for Atmospheric Research, Boulder, USA e-mail: [email protected]

A. Steiner et al. (eds.), New Horizons in Occultation Research,

DOI 10.1007/978-3-642-00321-9.4, © Springer-Verlag Berlin Heidelberg 2009

GNSS-LEO occultations. There are numerous other practical applications which can benefit by having turbulence measurements (and resultant climatologies) from GNSS-aircraft occultations. For example: space weather (Hajj et al. 1994), trans-ionospheric communication links (Secan et al. 1997), aviation safety and navigation systems (Conker et al. 2003; Cornman et al. 2004), the accuracy of ground-based GNSS receivers (Ganguly et al. 2004), as well as the effect of turbulence on measuring atmospheric state variables from GNSS-LEO occultations (Kursinski et al. 1997; Gorbunov and Kirchengast 2005).

The methods and results presented below are extracted from a much more thorough report prepared for the Boeing Corporation.

2 Theory—Weak Scattering

The essential theory for the problem of wave propagation through a random media is described in a variety of references, e.g., Tatarskii (1971) or Ishimaru (1997). Specific application to ionospheric measurements which follows directly from the basic theory has also been developed (cf. Yeh and Liu 1982; Bhattacharyya et al. 1992). Further application of this theory to the GNSS-LEO problem is discussed in Vorob'ev and Kan (1999) and Cornman et al. (2004). In the following, this latter reference is used as a starting place.

Since the transmitter and receiver are moving relative to each other, the origin of the coordinate system for the problem is taken to be at the receiver. In this case, a displacement vector p in the plane perpendicular to the line-of-sight (LOS) between the transmitter and receiver is given to the first order by p = Veffr, where t is the change in time over which the displacement occurs. Assuming that the turbulence is isotropic in the perpendicular plane, it can only be a function of the length of the displacement vector, ||p|| = ||Veff|| t, where ||Veff|| for a non-moving atmosphere and spherical wave propagation is given by

HVeffWII =

where n is the distance along the LOS from the transmitter (n = 0) to the receiver (n = R). VT is the velocity of the transmitter and VR is the velocity of the receiver, both computed as the projection of the full 3-D vectors into the plane transverse to the LOS. It can be shown that the frequency spectrum of a spherical wave measured in the plane perpendicular to the LOS, for the log-amplitude fluctuations of the index of refraction field, x, is given by ni+An/2

l|Veff(n)ll

dn, where it is assumed that the turbulence resides in a patch centered at n1 and of length An along the LOS; k is the wavenumber for the transmitted signal; Re refers to the real part; C2 is the turbulence intensity (also referred to as the structure constant); U is the confluent hypergeometric function of the second kind, and r is the gamma function (Abramowitz and Stegum 1972); x = 2nf/||Veff||; y = x2 + L—2; L0 is the integral length scale of the turbulence; a = n(R - n)/kR; and A = 0.033. A von Karman model for the isotropic index of refraction spectrum is used in Eq. (2) (Ishimaru 1997). Note that the Kolmogorov case is obtained by setting L0 ^ to. The Kolmogorov model is a power law at all wavenumbers. It is valid in the so-called inertial subrange of wavenumbers—far away from both the energy-producing scales (small wavenumbers) and the dissipation scales (large wavenum-bers, Tatarskii (1971)). So, theoretically it cannot be a valid model at the smallest wavenumbers (it goes to infinity). In this subrange, the index of refraction spectrum is only a function of the single intensity parameter, C2. The von Karman model has the same large wavenumber limit as the Kolmogorov model (power-law), but at small wavenumbers, it asymptotes to a constant (a function of both C;2 and L0), which ensures finite energy. The plane wave propagation case is given by setting a = (R — n)/k. For turbulence in the ionosphere, a different set of coefficients are required as well as changing the k2 factor to k—2, but the rest of the formulation is the same. If it is further assumed that the length of the turbulence patch An is small enough such that a mid-point approximation to the integral in Eq. (2) is an accurate one (given some criteria), then

( = An^k2ACfrùAn 4/3. j 2TO/3) _ ^ x VefKm) \sr (5/6)

It should be noted that typically An will not be small enough such that a midpoint approximation is accurate; hence, the patch must be broken into a number of smaller segments such that the approximation is acceptable within each one. A detailed discussion on this topic is beyond the scope of this paper (a comment is made in Sect. 4). However, to simplify the problem, for the results presented below it is assumed that the mid-point approximation is accurate enough.

3 Results

3.1 Log-Amplitude Spectra: Analytic Forms

Next, consider the sensitivity of Eq. (3) to the unknown parameters, L0 and n1. The other unknown, the combined parameter C^(n1)An will merely change the overall level of the curves. Figure 1a shows the log-amplitude spectra when L0 is held fixed and R — n1 (the distance along the LOS from the aircraft to the turbulence patch), is varied. The parameter values are: L0 = 3.0 km and R — n1 = {10, 50, 100, 200} km. Note that in the following analyses, it has been assumed that the angle between Fig. 1 (a) Log-amplitude spectra—showing the variation over different R — n values, for fixed L0 = 3.0 km. (b) Log-amplitude spectra—showing the variation over different Lo, for fixed R — m = 100 km

the projections of the transmitter and receiver velocity vectors onto the plane perpendicular to the LOS is zero. That is, the projections are parallel. Clearly, this is not the general case, but for this heuristic analysis it will suffice. Next, consider holding n1 fixed and varying L0. This can be seen in Fig. 1b, where the parameter values are R — n1 = 100 km and L0 = {0.1, 1.0, 3.0, 7.0} km. From Fig. 1a, it is clear that at low-frequencies the spectrum is highly sensitive to n1, with lessening sensitivity for increasing frequencies. From Fig. 1b, it can be seen that excepting for L0 = 0.1 km—a very small value for the upper atmosphere—the log-amplitude spectrum is fairly insensitive to the turbulence length scale. These considerations will be important when it comes to parameter estimation.

### 3.2 Simulation Studies and Parameter Estimation

An investigation into parameter estimation was performed, focusing on the logamplitude spectra. Assuming that the received signal is a sample from a Gaussian process, its spectrum will have an exponential distribution. Using this fact, numerous realizations of the spectra can be generated, i.e., by replacing each spectral point with a sample from an exponentially-distributed random number sequence. An example of this is shown in Fig. 2. The simulation is intended to represent real-world situations; hence, it was assumed that the received signal would be available at 50 Hz, and approximately 10 s worth of data would go into a single spectrum (512 samples were used).

For this exercise, the parameter to be estimated is C^nO^n. Two cases were examined, (i) holding n1 fixed and varying L0, and (ii) holding L0 fixed and varying n1. In each case, a single-parameter maximum-likelihood (ML) method was used to estimate C2(n1)^n. Let @g(f) = p\$g(f); where the subscript g refers to the function (here, the log-amplitude spectrum), and p is the parameter to be estimated. In this case, (f) is Eq. (3) with Cl(n1)^n set to one. Then given a measured—or

Fig. 2 Theoretical log-amplitude spectrum (black curve), and with values replaced by exponentially-distributed noise (gray dots)

10-2

10-2 10-1 100 101 frequency (Hz)

in this case simulated—set of spectral points, &g (fi), i = 1,..., N the maximum-likelihood estimate of p is given by

where the calculation is done over a set of N frequency values (Smalikho 1997).

Figure 3a illustrates results of the simulation exercise for the estimation of C2(n1)An values, given fixed (and known) n1, and varying L0 around the true value. For the data in Fig. 3b, L0 is fixed and n1 is varied. The "true" values for L0 and R - n1 are 3 km and 100 km, respectively. This corresponds to the curve with diamond symbols in Fig. 1a. The concept is to determine the sensitivity in estimating C2(n1)An given an uncertainty in the knowledge of one of the two other unknown parameters. That is, a realization of the log-amplitude spectrum is generated for a given set of input parameters {C2(n1)An, n1, L0}. This is <t'x(f) in Eq. (4). The model function, (f), is calculated from Eq. (3), (with C^(n1) An = 1), using the same value of n1 that was employed in generating <PX(f); but with values of L0 that range around the fixed value that was used in generating &x (f). These L0 values constitute the ordinate in Fig. 3a. The black horizontal and vertical lines are the true values of L0 and C^nOAn, respectively. Figure 3a is a two-dimensional histogram for the {C;2(n1)An, L0} parameter space, with the color representing the number of counts in that bin. One thousand realizations, such as the one shown in Fig. 2, were generated for each L0 value. It can be seen that over a large range of length scales, the errors in estimating C2(n1)An are mostly within ±5% of the true value—the spread being dominated by the statistics of the random realizations. These results confirm the qualitative discussion above regarding the relatively insensitivity of the log-amplitude spectrum with respect to the length scale. That is, given a two-parameter model where one of the parameters is known, and it is insensitive to the other one, it is not surprising to see these results. Fig. 3 (a) Distribution of C^CniMn estimates as a function of varying L0, while holding R — n1 = 100 km fixed. (b) Distribution of Cj2(n1)^n estimates as a function of varying R — n1, while holding L0 = 3 km fixed. The black horizontal and vertical lines are the true values of L0 and C^(n1)An, respectively. Note that the scale for the abscissa on the two plots is different

Fig. 3 (a) Distribution of C^CniMn estimates as a function of varying L0, while holding R — n1 = 100 km fixed. (b) Distribution of Cj2(n1)^n estimates as a function of varying R — n1, while holding L0 = 3 km fixed. The black horizontal and vertical lines are the true values of L0 and C^(n1)An, respectively. Note that the scale for the abscissa on the two plots is different

A positive bias is visible at the smallest L0 value (L0 = 0.1 km). While this is a somewhat unrealistic value for the upper troposphere, the analysis of this case is instructive. From Fig. 1b, it can be seen that at low frequencies, the spectrum for this length scale value is about a factor of two smaller than that for the spectrum with L0 = 3 km. Recall that the simulated spectrum is generated using the true value (3 km) and it is the model function that is given the erroneous value (0.1 km). Referring to Eq. (4), it can be seen why the high bias occurs. The model spectrum with the erroneous length scale is divided into the simulated spectrum, and since the former is smaller over a handful of low-frequency points, the ML estimator will produce a small positive bias.

Figure 3b illustrates the same procedure as above, but keeping L0 fixed and varying n1. For convenience, the ordinate of this figure is R — n1, the distance to the turbulent patch from the receiver. The black horizontal and vertical lines are the true values of R — n1 and C^nO^n, respectively. Due to the larger sensitivities of the spectrum to n1, the biases are even more apparent in this scenario. From Fig. 1a, it is obvious that the spectral points at low frequencies for the model spectra with the smaller R — n1 value (10 km) are significantly lower than those from the simulated spectrum (R — n1 = 100 km). In going from the low-frequency values on the 10 km curve to those on the 100 km curve, there is an increase by a factor of around 20. From the 50-100 km curve would result in an increase by a factor of around 2.5. Hence, the increasingly large positive biases for decreasing R — n1 values. This is seen in the R — n1 < 100 km region (above the black horizontal line) of Fig. 3b. Of course, the opposite is true: as the R — n1 values used in the model spectrum increase over those in the true (e.g., simulated) spectrum, the bias will be increasingly negative. For example, in going from the low-frequency values on the 200 km curve to those on the 100 km curve, there is a decrease by a factor of around 2.5. This is seen in the R — n1 > 100 km portion of the figure.

Estimate over all Freq. Estimate over High Freq. Data True Value C2n delta eta

Estimate over all Freq. Estimate over High Freq. Data True Value C2n delta eta :

Estimate over all Freq. Estimate over High Freq. Data -True Value C delta eta

:

Fig. 4 Histograms of the distribution of C;2(n1)An values from the ML estimator, using logamplitude spectrum. For the histograms in dark gray, all of the spectral points are used in the estimation; for the ones in light gray only the high frequency spectral points are used. In this case the true value of R — n1 was 100 km and the model value was 10 km. (a) Prior to iteration using n1 estimation, (b) After iteration using n1 estimation

A partial solution to the bias problem is to restrict the spectral points used in the ML estimator to the high-frequency regime—where the effect of n1 is lessened. The result of this exercise is shown in Fig. 4a. The histogram in dark gray is the distribution of C^(n1)An estimates, over 1000 realizations, using all of the spectral points. In this case, a value of R — n1 = 10 km is used for the model spectrum whereas the value for the simulated spectra is 100 km. This corresponds to the top line in Fig. 3b, where it can be seen that there is a significant positive bias, as well as a broader distribution. When the set of points used in the ML estimator is restricted to higher-frequency points, the bias and the spread are significantly reduced. This can be seen in the distribution in light gray. The input C;2(n1)An value is indicated by the vertical black line.

### 3.3 Iterative Parameter Estimation

As can be seen from Fig. 4a, there is still a non-trivial bias in the estimation of C2(n1)An even when restricting the spectral points used in the ML estimator to those at higher frequency values. It can be shown that the bias in the ML estimator is directly proportional to the right-hand side of Eq. (4), minus one. At the higher frequencies, the ratios on the right-hand side of Eq. (4) are smaller than those at lower frequencies, but there are many more points at the higher frequencies so that the sum of these ratios can be significant. The reason that these ratios are not one in this case is that n1 does not equal n. One approach to ameliorating this influence is to use a two-step method for calculating the estimates of ^(nOAn. The first step is to estimate C;2(n1)An from the method described in the previous section. That is, a value of n1 is guessed, and used in the ML algorithm. Since the length scale value is not very pertinent, choose a reasonable physical value. Given these two parameters, the log-amplitude spectrum is now a function of the single unknown parameter n1. A least-squares minimization method is then used to estimate n1, which in turn is utilized in place of the "guessed" value in a second iteration of the ML estimation for C2(n1)^n. In this case, the true value of R — n1 is 100 km and the model value is 10 km. The n1 estimates computed from the minimization method have a high bias. This is due to the positive bias in the first estimates of C^nOAn. Referring to Eq. (4), if the C^(n1)An estimates are too large, then the model function, (f, i^), must be smaller in order to minimize the difference between the data, &g (fi), and the product of the C2(n1)An estimates and the model. Making the model function smaller implies that larger n1, or equivalently, smaller R — n1 values, must be used (cf. Fig. 1a). Note that since the true value of R — n1 is larger than the original model one, (100 km vs. 10 km); estimating larger n1 values should redu