Detection of reflections in GNSS radio occultation measurements using phase matching

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. We can then identify reflection signatures even in cases where they are normally obscured by the direct signal’s influence on the phase matching amplitude. This signature enhancement is demonstrated by simulations and confirmed with real MetOp-A data. 10

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 numerical weather prediction (NWP) 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 transmitter, Earth, and receiver, as well as the short wavelength of the signals, result in a measurement with high vertical resolution.In many RO events, the atmosphere is sounded above a region of the Earth where the surface reflects the signal.This results in reflected parts of the signal also being received by the RO instrument in orbit.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 analysis of reflection signatures in GNSS RO data was performed by Gorbunov (Gorbunov, 2016), where the canonical transform method was used for processing the data.Some investigations have been made into whether these reflections can improve the quality of RO data (Cardellach and Oliveras, 2016).While occultations with recorded reflections that the orbital radii for the satellites are fixed, and that the Earth surface is completely smooth.In the appendix it is shown that the bending angle for a reflected ray (provided the Earth surface is smooth) is given by 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 The data from occultation events is collected from the COSMIC 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 get more accurate values of impact height, 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 translation 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).To make simulations similar to the simpler geometry of the WOP, the GNSS transmitter position is set at a distance corresponding to the mean value of the distance between the center of curvature and the actual GNSS position in the CDAAC data file.Similarly, the receiver satellite orbit is assigned a constant radius.To determine the length of the orbit (how deep the occultation goes), the minimum and maximum separation angles between transmitter and receiver are computed and used.These values can be modified to simulate a deeper occultation event.For an atmosphere, the high-resolution, 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 0 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 routinley used in WOP simulations (Gorbunov, 2016;Levy, 2000).

Surface reflections
In the low part of the occultation (SLTA around −80 km, impact height around 2.1 km), where the signal becomes blocked by the Earth, the magnitude of U also decreases.This decrease is very sharp in simulated data, and more varied in real signals.
Around the a value that corresponds to the surface refractivity, there is sometimes a spike in |U |, even if the tracking was lost higher up.This spike is an indicator of the reflected signal.If tracking of the signal goes all the way to the surface, the spike overlaps with the direct signal and is obfuscated, as seen in Fig.
(2). Figure ( 1) illustrates the relationship between the SLTA of the received signal and the specific rays for an exponential profile, represented by their impact heights.This shows that reflected rays can be embedded at higher SLTAs than the lowest direct rays.It also implies that reflections that are otherwise hidden below a strong direct signal can be revealed by only performing PM on an orbital segment above a certain SLTA.To clearly identify the reflection spike and suppress the influence of the lowest direct rays, we pass a segment of the received signal to the PM operator.The reflection spike then arises on its own.The reflection spike then becomes clearly separated from the other parts of the signal.The sliding window used in Fig.
(1) has a width of 10 km SLTA, and to avoid artifacts due to the abrupt cutting of the signal the outer 5% of the window edges were tapered using a Tukey window (Tukey, 1967).Not

Results
We present ten cases where PM is performed on real signals alongside their simulated counterparts based on ECMWF's colocated refractivity profiles.These cases are classified based on the sharpest gradient above 100 m, in the same manner as the reference dataset from ECMWF (Healy, 2012).Overall, the structure of |U | is similar to that of a step function, both  8) through ( 12) -give rise to a sharp, 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 shows 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.

Conclusions
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.
We can then identify reflection signatures even in cases where they are normally obscured by the direct signal's influence on the phase matching amplitude.This signature enhancement is demonstrated by simulations and confirmed with real MetOp-A data.
The events containing reflected signals presented here are hand-picked to illustrate that the method works on real data.There are still questions that can only be answered by larger volumes of data, e.g.how this method compares with other methods used to identify reflections.Furthermore we observe that the reflection spikes vary greatly in shape.Additionally, simulated reflection spikes always occur very close to the surface impact parameter, however 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 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 tends to increase close to the Earth surface.As no noise was added to the simulations, it is not clear if 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.
All these observations need further investigation before any conclusions can be drawn.

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 π radians due to the reflection.The total bending angle for a reflected ray is therefore given by α The integral for the optical path length becomes more complicated (although the derivation is straightforward) Taking the time derivative of this expression leads to the very same expression as Eq.A14.Hence, the stationary phase point again occurs where a(t g ) = a g .At this point we have This is the term that appears in the phase of the phase matching function U (a g ).Taking the derivative with respect to a g leads to d da g (S(t g ) − S g (t g , a g )) = −2a g Which is the bending angle for a reflected ray as given in Eq. (A21) Consequently the phase matching method works in the exact same way for direct and reflected rays.It should be stressed that these derivations are only valid when the Earth surface can be considered smooth.When the surface is not smooth the incoming ray will change impact parameter upon reflection.
Due to this the expression for the optical path length becomes a function of the old and new impact parameter, and the simple geometrical ray model used in the phase matching method cannot lead to a stationary phase point.This is basically a case of multipath in the impact parameter domain.It may be argued though that this is of little consequence for real measurements since the occultation measuring instrument will not record signals that deviate too strongly from direct rays, as they quite rapidly become heavily Doppler shifted with increasing reflection angle.For this reason reflected signals will only be seen at impact parameters that are very close to the value at the Earth surface.These rays are of grazing incidence, and under such circumstances the surface may always be considered as flat (Beekmann and A, 1963).

Figure 1 .
Figure 1.Using a sliding window for the signal yields PM amplitudes that correspond to the reflection model.Left shows a case with weak refractive gradients, right shows a case with strong refractive gradients.
performing such tapering would result in unwanted noise, seen in the truncated |U | in Fig. (2).
Figure2.|U | for a simulated signal that is propagated to a similar orbit as the real measurement (red), and one that is propagated to a deep orbit (blue).For illustration the truncated signal is not tapered, the dramatic oscillations seen slightly around 4 km is produced.