Analysis of reflections in GNSS radio occultation measurements using the phase matching amplitude

Abstract. It is well-known that in the presence of super-refractive (SR) layers in the lower troposphere inversion of GNSS radio occultation (RO) measurements using the Abel transform yields biased refractivity profiles. As such it is problematic to reconstruct the true refractivity from the RO signal. Additional information about this lower region of the atmosphere might be embedded in reflected parts of the signal. To retrieve the bending angle, the phase matching operator can be used. This operator produces a complex function of the impact parameter, and from its phase we can calculate the bending angle. Instead of looking 5 at the phase, in this paper we focus on the function’s amplitude. The results in this paper show that the signatures of surface reflections in GNSS RO measurements can be significantly enhanced when using the phase matching method by processing only an appropriately selected segment of the received signal. This signature enhancement is demonstrated by simulations and confirmed with ten hand-picked MetOp-A occultations with reflected components. To validate that these events show signs of reflections, radio holographic images are generated. Our results suggest that the phase matching amplitude carries information 10 that can improve the interpretation of radio occultation measurements in the lower troposphere.

1 Introduction GNSS radio occultation (RO) is a technique used for sounding the Earth's atmosphere.Assuming spherical symmetry of the atmosphere, the bending angles of GNSS signals passing through the atmosphere can be found and assimilated into nu-merical weather prediction systems.The bending angle measurements contain valuable information due to their relation to the atmosphere's refractivity, which can yield information about humidity, temperature and pressure (e.g., Kursinski et al., 2000;Yunck et al., 2000).The geometry of the transmitter, Earth and receiver, as well as the short wavelength of the signals, results in a measurement with high vertical resolution.In many RO events, the instrument in orbit receives reflected components of the signal as well as direct ones.Boniface et al. (2011) have shown that reflected signals contain meteorological information.A method to detect these reflected components was suggested by Hocke et al. (1999) and has later been used on real data and shown to work (e.g., Beyerle et al., 2002;Pavelyev et al., 2002).This method uses a radio hologram generated by subtracting a ray-traced reference field from the received signal.An effort to flag occultation events where reflections are present is described by Cardellach and Oliveras (2016), based on a supervised learning approach classifying such radio holographic images.It has since then been employed by the Radio Occultation Meteorology Satellite Application Facility (ROM SAF) to flag millions of occultation events.Cardellach and Oliveras also investigate whether knowledge of these reflections can improve the quality of RO data, but Healy (2015) concludes that a binary reflection flag is probably inappropriate for assimilation purposes.Gorbunov (2016) proposes a technique based on the canonical transform to retrieve bending angle profiles of reflected rays and achieves a good agreement with the ROM SAF database.
When processing a RO signal, we use the phase matching (PM) operator (Jensen et al., 2004), which outputs a complex |U | for a simulated signal that is propagated to a shallow orbit and one that is propagated to a deeper orbit (blue).For illustration the truncated signal is not tapered, producing dramatic oscillations around 4 km.
function of impact parameter whose phase is proportional to the signal's bending angle.Although the amplitude of this function may contain valuable information as well, it has not been appropriately investigated.
In this paper, we compare the PM amplitudes of simulated measurements to real measurements made by the GRAS instrument aboard MetOp-A.We demonstrate signatures corresponding to surface reflections and truncate the received signal to distinguish them from the much stronger signatures of the direct components.We justify this truncation using simulated data and a simple geometric model for which reflected components are expected to appear in the signal, and we provide radio holographic images as a means of validation.Finally we discuss the difference in structure between simulated and real signals, as well as the potential future uses of the PM amplitude.In the Appendix we show that the PM operator admits reflected signal components.

Phase matching
Jensen et al. define PM as an operator that transforms a complex signal of time u(t) to impact parameter space: (1) Here, k 0 is the wavenumber in vacuum associated with the carrier frequency of the GNSS transmitter, S(t, a) is the optical path length of a model ray path, and a is impact parameter.The derivative of arg(U (a)) with respect to a is proportional to the bending angle.As this integral is defined for any impact parameter a, it is important to determine at which a the function does not contain relevant information anymore.
The lowest a (which we call a min in this paper) corresponds to a ray that is tangential to the Earth's surface.To make sure that all information in the signal is mapped to impact parameter space, we compute the U for impact parameter values going all the way to the Earth's surface.This also ensures that we include reflected rays in the U function.It is not obvious that the PM method should work for reflected rays, as the geometrical model ray path is constructed for a direct ray, but in the Appendix it is shown that the standard PM technique admits reflected rays as well, provided the Earth surface is smooth.

A model for reflected rays
In Fig. 1 we use data from two real MetOp-A measurements to demonstrate that some of the features we see in the amplitude for the complex function U are caused by reflections.By overlaying the predicted straight line tangent altitude (SLTA) for the direct rays (blue line) and reflected rays (red line) as a function of impact height the reflection is illustrated.To generate these SLTA plots we use the co-located ECMWF refractivity profiles.The black lines show the amplitude of the U function when passing segments of the signal to the PM operator using a 10 km sliding window.The relationship between reflected bending angle and impact parameter is wellknown (see, e.g., Pavelyev et al., 2011;Gorbunov, 2016) and can be described as where α is the bending angle, R E the Earth radius of curvature, and n the refractive index as a function of radius.The bending angle of a direct ray is given by the same expression without the last two terms.The integral can be evaluated numerically using a number of techniques, and we employ the method described in Rasch (2014).The SLTA for fixed values of the LEO and GNSS orbital radii is given by  where r L and r G are the LEO and GNSS orbit radii, and θ is the separation angle between the satellites, given by 4 MetOp-A data The data from occultation events are collected from the COS-MIC Data Analysis and Archive Center (CDAAC) web interface, specifically day 2007.274with the metopa2016 designation, indicating reprocessed measurements from MetOp-A.The signal amplitude, excess phase and orbit data needed for PM are all found in the atmPhs files.In these files, the orbit coordinates are given with the Earth's center of mass as the point of origin.To fulfil the assumption of local spherical symmetry, we translate the coordinates so that they instead consider the center of curvature of the Earth at the occultation point.This is done by collecting center of curvature data from the corresponding atmPrf files.As the atmPrf files contain bending angle and impact height values, these were used to control the accuracy of the PM implementation.
For simulating a GNSS signal as it passes through the atmosphere, a wave optics propagator (WOP) is used with the multiple phase screen technique (see, e.g., Benzon et al., 2003;Benzon and Gorbunov, 2012;Rasch, 2014).The WOP uses a simpler, two-dimensional geometry where the GNSS transmitter is stationary.The three-dimensional geometry data in the atmPhs files are thus projected onto a plane where the LEO orbit is defined by the separation angle.Both transmitter and receiver are assigned a constant radius, and the Earth is assigned the radius of curvature for the particular event.The separation angles can be modified to simulate an occultation event tracked to a lower SLTA than its  corresponding measurement.For an atmosphere, the highresolution, co-located refractivity profile from ECMWF also provided by CDAAC, is used (echPrf files).These profiles do not go all the way to the ground -the last bit is extrapolated linearly.To simulate surface reflections, we set the electromagnetic field to zero on all parts of the phase screens that lie inside the Earth.Although it is not clear whether this method has a solid physical basis it appears to give quite accurate results and is routinely used in WOP simulations (Gorbunov, 2016;Levy, 2000).
As a frame of reference, radio holographic images are produced by constructing a "beating function" from the excess phase, signal amplitude and a reference field and using the short-time fast Fourier transform with a sliding window, similar to Boniface et al. (2011).

Surface reflections
In the lowest part of an occultation (SLTA around −80 km, impact height around 2.1 km), where the signal becomes shadowed by the Earth, the magnitude of U also decreases.The lowest direct ray becomes diffracted by the Earth's surface and gradually decreases in magnitude over a region corresponding to the Fresnel zone size, which is seen clearly in simulated data and frequently less clearly in measured data.This transition occurs over a few hundred meters (Kursinski et al., 2000), around SLTA min , and a min , whose values are determined by the Earth radius and the refractive index at ground.Quite often when the signal is lost at an SLTA value that is substantially higher than SLTA min we see a spike in |U | around a min .This spike is most easily explained as a reflection.If tracking of the signal goes all the way to the surface this spike overlaps with the direct signal and is obfuscated.In Fig. 1 we see that the reflected signal and the di-  The nature of the PM operator is such that sharp discontinuities in the signal introduce significant noise in the |U | function.To avoid these, we taper the signals using a onesided Tukey window (Bloomfield, 2004).In the case of the sliding window used in Fig. 1, we use a two-sided window.The noise generated by omitting any sort of tapering is demonstrated in Fig. 2, where the simulated signal ends abruptly.

Results
We present 10 cases where PM is performed on real signals alongside their simulated counterparts based on ECMWF's co-located refractivity profiles.These cases are classified in  Simulations on cases classified as category 3 and 4shown in Figs. 8 through 12 -give rise to a sharply negative spike at an impact height corresponding to the sharp gradient in the refractivity profile.This structure cannot be found in the real data.Moreover, the real data show a high level of noise that is not found in simulated data.
We note that there are peculiar oscillating structures in the real data, particularly in Fig. 4  Furthermore we note that in all cases but Fig. 11, the reflection spike appears some distance below the a min calculated from the ECMWF profiles.
While these 10 hand-picked cases all clearly contain reflected components, there are several measurements in which the reflected parts cannot be distinguished.This is typically either because the measurement was not deep enough or be- .This is corroborated by the radio holographic images, which also show very clear reflection patterns in those cases.The bright yellow trail around 0 Hz is the direct signal, whereas the more faint trail going off to negative frequencies (and appearing again at positive frequencies due to aliasing) is the reflected signal.The signal has been divided everywhere by the average signal amplitude in the window.This is done to make the reflection and the direct signal for low SLTA stand out against the strong signature of the direct signal for high SLTA.This also often causes the hologram to give the appearance of containing strong broadband noise for low SLTA, but it is actually the amplification of weak noise.

Conclusions
The results in this paper show that the signatures of surface reflections in GNSS RO measurements can be significantly enhanced when using the PM method by processing only an appropriately selected segment of the received signal.We can then identify reflection signatures even in cases where they are normally obscured by the direct signal's influence on the PM amplitude.This signature enhancement is demonstrated by simulations and confirmed with MetOp-A measurements.
For further validation of the reflections, radio holographic images are provided for comparison.For the cases presented, clear and sharp reflection spikes in the PM amplitude are corroborated with clear reflection patterns in their corresponding radio holograms.Weaker or less distinct spikes have more noise in their corresponding radio holograms.
The events containing reflected signals presented here are hand-picked to illustrate that the approach can be useful for identifying reflections in real measurements.There are still questions that can only be answered by larger volumes of data, e.g., if the PM amplitude can possibly be used to retrieve information related to the reflected signal.Furthermore, we observe that the reflection spikes vary greatly in shape.Additionally, simulated reflection spikes always occur very close to the surface impact parameter, but in real measurements this is not always the case.The reason for these variabilities needs to be investigated.At this point we have only analyzed events over water, since we expect those reflections to be much clearer and numerous compared to events over land.
When comparing the PM amplitudes of real and simulated signals, some interesting discrepancies were found.First, the characteristic dip in |U | associated with a region of sharp gradient in the refractivity could not be identified in the real data, being completely drowned in noise or simply not present.Second, the levels of noise in the |U | function tend to increase close to the Earth surface.As no noise was added to the simulations, it is not clear whether this is due to instrument noise or atmospheric disturbances.Third, peculiar oscillatory features at quite distinct heights were seen in the real data but not at all in the simulated.These features are not understood at present, but it is likely that they have to do with horizontal gradients in the refractive index and with aliasing of the reflected signal onto the direct signal due to low sampling rate.What horizontal gradients do to the PM amplitude has not been investigated, but it is likely that it will decrease the amplitude and cause additional noise, perhaps resulting in these peculiar oscillations.Regarding the aliasing of reflected signals, it is quite clearly seen in some of the radio holograms that the reflected signal is still strong when it has been aliased and again approaches the direct signal in frequency.This may cause a degradation of the PM method (the integral over the stationary phase point) at the specific impact height that belongs to the direct signal at the point where the direct and reflected aliased signal match in frequency.
Data availability.The data from occultation events are collected from the COSMIC Data Analysis and Archive Center (CDAAC) web interface, found at the URL http://cdaac-www.cosmic.ucar.edu/cdaac/index.html.The time derivative of S is given by Likewise, the time derivative of S g is Hence, the stationary phase point occurs where a(t g ) = a g .At that point the difference in optical path length becomes and the phase matching integral is given by where C(t g ) is an amplitude factor depending on the signal amplitude and phase in the region around the stationary phase point.The bending angle as a function of impact parameter is thus found by taking the derivative of the phase of the function U with respect to a g , i.e., A2 Reflected rays For rays suffering reflection the Bouger's rule still applies, but the ray never reaches the point where φ = π/2.Instead the ray is reflected at the point where r = R E , where R E is the Earth radius of curvature.Using the definition R E n(r E ) = a E we find the angle the ray makes with the radial vector at reflection to be Here we naturally assume that a < a E ; otherwise the ray would never reach the surface and be reflected.Since we assume the surface to be completely smooth, the radial vector is parallel to the surface normal, and since the incidence angle with respect to the surface normal is equal to the reflected ray angle with respect to the surface normal, we find that the ray angle after reflection is We conclude that the ray suffers a negative bending of π − 2φ E radians due to the reflection.The total bending angle for a reflected ray is therefore given by α(a) = −2a The integral for the optical path length becomes more complicated (although the derivation is straightforward): where dot signifies derivative with respect to time.The last term can be rewritten using Eq.(A21): Figure 1.Using a sliding window for the signal yields PM amplitudes that correspond to the reflection model.Panel (a) shows a case with weak refractive gradients; panel (b) shows a case with strong refractive gradients.

Figure 3 .
Figure 3.An event classified as category 1. SLTA min is at approximately −73 km, and the signal is truncated at SLTA = −50 km.The reflection spike is located at an impact height of 1919 m, and a min is located at an impact height of 2017 m.

Figure 4 .
Figure 4.An event classified as category 1. SLTA min is at approximately −73 km, and the signal is truncated at SLTA = −60 km.The reflection spike is 1881 m, and a min is 1992 m.

Figure 5 .
Figure 5.An event classified as category 1. SLTA min is at approximately −86 km, and the signal is truncated at SLTA = −60 km.The reflection spike is located at an impact height of 1948 m, and a min is 2008 m.

Figure 6 .
Figure 6.An event classified as category 1. SLTA min is at approximately −65 km, and the signal is truncated at SLTA = −50 km.The reflection spike is located at an impact height of 1849 m, and a min is 1888 m.

Figure 7 .
Figure 7.An event classified as category 1. SLTA min is at approximately −72 km, and the signal is truncated at SLTA = −60 km.The reflection spike is located at an impact height of 1940 m, and a min is 1993 m.

Figure 8 .
Figure 8.An event classified as category 3. SLTA min is at approximately −87 km, and the signal is truncated at SLTA = −70 km.The reflection spike is located at an impact height of 2075 m, and a min is 2204 m.
Figure 9.An event classified as category 3. SLTA min is at approximately −84 km, and the signal is truncated at SLTA = −50 km.The reflection spike is located at an impact height of 2057 m, and a min is 2182 m.

Figure 10 .
Figure 10.An event classified as category 3. SLTA min is at approximately −85 km, and the signal is truncated at SLTA = −50 km.The reflection spike is located at an impact height of 1928 m, and a min is 2200 m.

Figure 11 .
Figure 11.An event classified as category 3. SLTA min is at approximately −99 km, and the signal is truncated at SLTA = −50 km.The reflection spike is located at an impact height of 2363 m, and a min is 2301 m.

Figure 12 .
Figure 12.An event classified as category 4. SLTA min is at approximately −103 km, and the signal is truncated at SLTA = −50 km.The reflection spike is located at an impact height of 1994 m, and a min is 2232 m.
at 6 km, Fig. 7 at 4 and 6 km, Fig. 8 at 7.5 km, Fig. 11 at 5 and 9 km and Fig. 12 at 5 and 9 km.These oscillating structures are not found in the simulated data.