Reflected ray retrieval from radio occultation data using radio holographic filtering of wave fields in ray space

Linear and non-linear representations of wave fields consti tute the basis of modern algorithms for analysis of radio occultation (RO) data. Linear representations are impleme nted by Fourier Integral Operators, which allow for high-re solution retrieval of bending angles. Non-linear representations i nclude Wigner Distribution Function (WDF), which equal the p seudodensity of energy in the ray space. Representations allow fo r filtering wave fields by suppressing some areas of the ray spa ce and mapping the field back from the transformed space to the in itial one. We apply this technique to the retrieval of reflect ed 5 rays from RO observations. The use of reflected rays may incre ase the accuracy of the retrieval of the atmospheric refract ivity. Reflected rays can be identified by the visual inspection of WDF or spectrogram plots. Numerous examples from COSMIC data indicate that reflections are mostly observed over oceans or snow, in particular, over Antarctica. We introduce the refle ction index that characterizes the relative intensity of the refle cted ray with respect to the direct ray. The index allows for t he automatic identification of events with reflections. We use t he radio holographic estimate of the errors of the retrieved bending 10 angle profiles of reflected rays. A comparison of indices eval uated for a large base of events including the visual identifi ca on of reflections indicated a good agreement.


Introduction
A clear signature of signals reflected by the Earth's surface was revealed as early as the beginning of 21st century, by means of the radio holographic analysis of CHAMP radio occultation (RO) data (Beyerle and Hocke, 2001;Beyerle et al., 2002).Similar patterns were also found in Microlab-1 GPS/MET data (Gorbunov, 2002b, c).It was pointed out that the utilization of reflected signals can be useful for the enhancement of the retrievals.Reflections are mostly observed above water (in this case, ocean) or snow (Antarctica).Another application of reflected signals is linked to the altimetry (Cardellach et al., 2004).
Currently, the main means of identification of reflections remains the radio holographic analysis in sliding apertures (Gorbunov, 2002c;Cardellach et al., 2009Cardellach et al., , 2010)), which has also been used to extract the reflected fields (Aparicio et al., 2017).Alternatively, the reflected signals can also be identified with techniques based on the Wigner Distribution Function (WDF) (Gorbunov et al., 2010(Gorbunov et al., , 2012)).However, it is known that the techniques based on different approximations for the Fourier Integral Operator, developed for retrieval of RO in multipath areas, are also capable of retrieving the reflected part of the bending angle (BA) profile.These techniques include: canonical transform (CT) methods based on the concept of Fourier integral operators (FIO) (Gorbunov, 2002a;Gorbunov and Lauritsen, 2004), full spectrum inversion (FSI) (Jensen et al., 2003), and phase matching (PM) (Jensen et al., 2004).
In this paper, we discuss the algorithm of reflected ray retrieval based on the modification of the CT method.We use the WDF technique for the visualization purposes only, because it does not provide accuracy of the bending angle retrieval comparable with the CT technique.The algorithm is based on the filtering in ray space.In the CT method, the original wave field observed as a function of observation time t is projected into the impact parameter domain.The field in the transformed space is a function of impact parameter p, and the direct and reflected rays can be classified according to their value of p.The reflected rays have an impact parameter value below, corresponding to the geometric optical shadow.Therefore, the reflected field component can be separated in the impact parameter domain.For further processing, it is more convenient to project the reflected field component back into the time domain (Cardellach et al., 2009(Cardellach et al., , 2010)).For this projection we use the inverse FIO.Generally, our approach is similar to that developed in Cardellach et al. (2009Cardellach et al. ( , 2010)), where the spectrograms were filtered in order to remove direct rays and inverted.However, the CTbased approach allows for constructing efficient numerical algorithms.In relation to the problem of the reflected ray retrieval, this approach has the same advantages of the radio holographic (sliding spectral) method (Beyerle and Hocke, 2001;Beyerle et al., 2002), as it is more accurate and computationally more efficient (Gorbunov, 2002c).This work was published as a technical report (Gorbunov, 2016).
The paper is organized as follows.In Sect.2, we discuss possible approaches to the reflection retrieval.We present typical examples of reflections observed in COSMIC data, which allow us to choose the most adequate method of reflection retrieval.In Sect. 3 we discuss the phase model for reflected rays, and its use for the definition of the reflection index that is a measure of the reflection strength.Next we discuss the filtering algorithm based on FIO that allow separating reflection from RO observations.Section 4 shows a few examples of COSMIC events with different reflection indices.In Sect.5, we make our conclusions.

Possible approaches to reflection retrieval
In this Section, we will discuss possible retrieval algorithms for the retrieval of reflected rays.As compared to direct rays, reflected rays are characterized by a lower amplitude and by a bending angle profile rapidly increasing with impact height.This makes it difficult to apply the CT algorithm directly.Below we discuss three algorithms: (1) direct application of CT technique, (2) the modified CT technique where the impact parameter is replaced with a linear combination of impact parameter and bending angle, and (3) composition of a filter in the impact parameter space that suppresses the direct rays with subsequent geometric optical inversion in the time domain.

Reflected rays retrieval in impact parameter domain
The measured wave field has the following form: where A(t) is the observed amplitude, S 0 (t) is the satelliteto-satellite distance evaluated from the orbit data, S(t) is the observed atmospheric excess phase.The smoothed excess phase S(t) and its derivative S (t) are obtained by the filtering of the measured excess phase S(t).The smoothed relative Doppler shift d(t) is obtained as follows: where  1994).The ancillary function f (t) is evaluated as follows (Gorbunov and Lauritsen, 2004): The new coordinate is determined as follows (Gorbunov and Lauritsen, 2004): where ϒ 0 is a constant determined in such a way that ϒ(t) ≥ 0 for the observation time interval.We evaluate the integral of f (t): Using this function, we evaluate the vacuum and observed phase path with subtracted model as follows: where R E is the Earth's local curvature radius.Subtraction of R E ϒ(t), a linear function of ϒ, from the phase corresponds to the reduction of the frequency, which equals the impact parameter, by a constant of R E .All the functions of time t can also be looked at as functions of the new coordinate ϒ.The FIO is defined as follows (Gorbunov and Lauritsen, 2004):

GPS
where a (p) is the amplitude function, whose definition can be found in Gorbunov and Lauritsen (2004).The variable p is the approximate impact height (impact parameter with subtracted R E due to the definition of S (M) 0 (t)).The transformed field ˆ 2 u ( p) is represented as follows: where A ( p) is the amplitude of the transformed field and ϕ ( p) is its accumulated phase.The frequency variable ϒ is defined in such a way that it is always positive, where any rays may be expected.This simplifies the evaluation of the accumulated phase ϕ ( p).The amplitude function is evaluated using ϕ ( p), which is the reason why field u ( p) is first evaluated up to this factor.The analysis of the amplitude of the field in the transformed space allows for the determination of the shadow border impact height p E (Gorbunov, 2002a;Gorbunov and Lauritsen, 2004;Jensen et al., 2004), as shown in Fig. 1.Practically, because the energy of reflected rays is much smaller than that of direct rays, this will also be the border between direct and reflected rays.The CT2 algorithm evaluates the filtered phase derivative, separately for p < p E and for p > p E , which we denote as follows: where subscript D stays for direct rays, and subscript R stays for reflected rays.However, it will be necessary to implement an additional option that specifies the filter width for reflected rays.This is explained by the fact that the impact parameter interval for reflected rays is usually as narrow as 100-200 m.This requires a narrow filter window of about 20 m, while the typical setting for processing direct rays in the lowest troposphere is 250 m.Such a filter will not be able to effectively suppress random noise.

Reflected rays retrieval in modified impact parameter domain
In order to circumvent the problem of the reflected ray retrieval in the impact parameter space, caused by steep increase of bending angle profile for reflected rays, we consider the following modification of the CT method.
The FIO (8) corresponds to the following linear canonical transform (Gorbunov and Lauritsen, 2004): where η is the eikonal derivative (momentum) of the original observed field u(t), and ξ is the momentum of the transformed field.This transform can be modified in order to use another coordinate: where β is a tunable parameter.Lines of constant coordinate p are shown in Fig. 1 in red.This indicates that the bending angle profile of reflected ray is less steep in this coordinate space for β > 0, because: www.atmos-meas-tech.net/11/1181/2018/Conversely, β cannot be made too large, because in this case, the direct rays may overlap with the reflected rays in the modified space.
The modification of the integral transform is straightforward.The modified canonical transform is written as follows: Using the modified function f (ϒ) instead of the original one defined in Eq. ( 3), we obtain the expression for the modified FIO ˆ 2 .The advantage of this approach is that it can be implemented by a relatively small modification of the existing CT2 algorithm.Its disadvantage is the presence of a tunable parameter β, whose optimal value is unknown in advance and may vary from event to event.

Reflected rays retrieval in time domain
The CT2 algorithm is designed for the retrieval of bending angle profiles in multipath areas, where the profiles are nonmonotonic.This is not the case for bending angle profiles of reflected rays, which for spherically symmetric media, must monotonically increase, as will be discussed in Sect.3.1.This makes it convenient to retrieve the dependence p (ε) rather than ε(p).On the other hand, it is more straightforward to formulate the retrieval algorithm in the time domain, as illustrated by Fig. 1.In the time domain, in presence of reflection, there is always multipath propagation due to the interference of direct and reflected rays.But the field component related to the direct rays can be effectively removed.To this end, we can use the impact parameter domain, where the direct and reflected rays are clearly separated by the border impact height of p E .Therefore, we can form the following field u R (t) in the time domain that only contains the reflections: where ˆ 2 is the FIO, ˆ −1 2 is its inverse, and θ is the thetafunction, which takes the value of 1 for positive arguments and 0 for negative arguments.This function can then be processed using the standard geometric optical (GO) technique.The advantage of this approach is that it is free of tunable parameters, and the final retrieval of reflected rays is performed in the time domain, which is the optimal coordinate for the manifold of reflected rays.This is illustrated by Fig. 1, which shows a very typical bending angle profile of reflected rays, where at each moment of time only one reflected ray is observed.This will also be illustrated by examples of experimental data.The examples of reflections allow for the following conclusion.Observed reflections indicate a very rapid increase of the bending angle of reflected rays, ε R as a function of impact parameter p. Dependence ε R (p) is mostly confined in a narrow impact parameter interval of about 100 m.Often ε R (p) is a multi-valued function, which is due to the impact parameter perturbations in a horizontally non-uniform atmosphere (Gorbunov and Kornblueh, 2001;Gorbunov and Lauritsen, 2009).This indicates the method of choice should be the time domain retrieval, preceded by the extraction of the reflected signal, by using the filtering in the impact parameter space.

Examples of reflections in COSMIC observations and discussion
3 Reflection retrieval implementation

Phase model for reflected RO signals
The phase model will play an important role in our further discussion.Here we describe the algorithm for the evaluation of the phase model for reflected RO signals.Given a spherically symmetric model of the neutral atmosphere n M (r), where r is the distance from the Earth's curvature center, the corresponding bending angle profile for reflected rays is expressed as follows (Aparicio et al., 2017): where θ (t) is the satellite-to-satellite angle with respect to the local curvature center, r Tx,Rx (t) are the radial coordinates of the satellites, hereinafter index Tx indicating the transmitter and index Rx indicating the receiver.The equation is solved for time t for each prescribed impact parameter.This allows for the determination of impact parameters as function of time, p MR (t).Dependence p MR (t) is always single-valued for reflected rays, because reflected bending angle profiles are monotonic and do not result in multipath propagation.This is illustrated by Fig. 1 and explained by Eq. ( 18), where the derivative of the second, reflective term proves to be much stronger than that of the first, refractive term, for realistic atmospheric conditions.In order to show this, we recall that strong multipath effects are caused by superrefraction layers, and they are the strongest for spherically layered medium.Under these assumptions, we can write, where r E is the Earth's radius.Then we can derive the expression for the derivative of the bending angle as: Now, assuming that the strongest perturbation comes from a superrefraction layer with a thickness of r, critical refractivity gradient of −r −1 E and located at an altitude of h SR , we arrive at the expression, where p = p E − p.Assuming that h SR is about the PBL height (i.e.1.5 km), and the superrefraction layer thickness is about 0.2 km, and p < 0.2 km, we see that r/h SR 2 √ h SR / p and therefore, dε R (p)/dp > 0. Because the bending angle profile of reflected rays is monotonic, there is only one ray at each moment of time, as illustrated by Fig. 1.Given satellite coordinates x Tx,Rx (t), the ray directions at the satellites, unit vectors u Tx,Rx (t) are inferred from p(t) using the geometrical relationships: which express the fact that rays lie in the vertical occultation plane, and the impact parameter has the same value at the transmitter and at the receiver.Using the satellite velocities V Rx,Tx (t), we find the relative Doppler frequency shift d MR (t): The excess phase is obtained by integrating the Doppler shift: where d (0) (t) is the vacuum Doppler shift for the direct rays, evaluated from Eq. ( 24), by inserting unit vector u Tx,Rx (t) corresponding to satellite-to-satellite straight-line direction.An example of reflected excess model phase is shown in Fig. 7.

Radio holographic index of reflections
The idea of flagging radio occultation with an index of the strength of the reflected ray can be expressed as follows.Although the amplitude of the reflected signal is weak as compared to the direct signal, the instant frequencies of the reflected signal concentrate around the instant frequencies of the model reflected signal.Therefore, we can use the model reflected signal exp (ikS MR (t)) as the reference signal and evaluate the radio holographic spectrum as follows: where A (t) and S (t) are the amplitude and the excess phase of the observed signal.As the reference signal, we use the smoothed reflected signal excess phase S R (t) rather than the model S MR (t).We apply the sliding polynomial smoothing with a window of about 1 s.This modification makes the radio holographic spectrum sharper, while its maximum reflection is located closer to the zero frequency.The integration here covers the time interval, for which we can evaluate the reflected excess phase model.Each frequency ω can be transformed to the equivalent impact parameter.This allows for considering the spectrum as a function of impact parameter related to middle point t 0 of the time interval.Moreover, it is convenient to introduce the reference value p 0 of the impact parameter, corresponding to frequency ω 0 = ṠMR (t 0 ).The spectrum can be considered as function u R ( p) of relative impact height p = p − p 0 .An example of radio holographic spectrum is shown in Fig. 8.The spectrum indicates a distinct spike near p = 0, corresponding to reflection.The presence of reflection is also confirmed by the Wigner function plot in Fig. 2.
We define the reflection index I R as a functional of the radio holographic spectrum, depending on a series of empirical parameters: where u max is the maximum of the spectral density taken within the interval of p ∈ [−0.1 km, 0.1 km], p max is the location of the spectral maximum of the reflection, u ave is the spectral density averaged over the interval of p max − 0.3, p max + 0.3 , u bkg is the background (noise level) spectral density estimated by averaging over the interval of p ∈ [1.0 km, 2.0 km], and α is the regularization parameter, p M (t) is the dependence of the impact parameter on the model reflected signal vs. time, and δp(t) is the radio holographic error estimate of the impact parameter.The regularization allows for suppressing random maxima at the noise level, if the reflection is weak or absent.We estimate the background spectrum density u bkg from the impact parameter interval of [1.0, 2.0] km, where a signal from direct ray is present.The regularization strength is controlled by the parameter α reflection is only identified if u max significantly exceeds both u ave and α u bkg .The optimal value of α was empirically estimated to be about 0.2.The additional exponential factor in the definition of I R penalizes profiles deviating too much from the model.The averaging in this factor is spread over the whole domain, where the reflected bending angle profile is evaluated.This reflection index definition is easy to implement and computationally inexpensive.
The index characterizes the strength of the spectral spike and suppresses random spikes at noise level.The value of I R = 0.25 corresponds to a flat radio holographic spectrum, i.e. a definite absence of reflection.For the illustrative event considered above, the index has a value of 17.9.This index can therefore be used to identify presence of reflected signals in RO data.

Filtering in impact parameter space
In order to extract the reflected field component u R (t), we implemented the filtering in the impact parameter space.In Fig. 2, we see that the reflected ray is observed both around impact height of 2 and 10.5 km.The latter originates from aliasing, where the Doppler frequency shift of the reflected ray deviates from the direct ray excess phase model by more than a half of the receiver band width, which equals 50 Hz.The impact parameters difference p alias between non-aliased and aliased components for typical observation geometry is about 8-10 km.The exact value p alias for a specific event is evaluated by using the GO relationship ( 23) and (24).To this end, we evaluate impact height from the where χ ( p) is equal to unity inside the impact height interval of p E − p R , p E and the corresponding aliased interval of p E + p alias − p R , p E + p alias .The width p R of these intervals is set to 1 km.Outside these intervals, we employ the Gaussian apodization with a characteristic width of δp R = 0.2 km.Apodization allows for avoiding sharp boundaries of the filtering function, improving filter quality.For the estimate of impact parameter uncertainty, we employ radio holographic analysis.Gorbunov et al. (2006): the uncertainties are estimated as the widths of sliding radio holographic spectra, in terms of the impact parameter.In order to validate our retrieval algorithm and reflection index definition statistically, we performed a comparison of our retrievals with the ROM SAF reflection flag database (ROM SAF, 2016;Cardellach and Oliveras, 2016).The database contains occultation events classified into three categories: (1) no reflection, (2) reflection, and (3) unclear.The events are accompanied by the Supporting Vector Machine (SVM) index (Cardellach et al., 2009(Cardellach et al., , 2010) ) based on the radio holographic analysis and supervised learning method.The reflection index is a functional of a process containing a random component (noise, turbulence effects etc.) and deterministic regular structures (direct and reflected ray).The index characterizes the intensity and the sharpness of the reflected ray.Being a function of a random process, the index itself is a random quantity with its own distribution.No index under these conditions, can exactly characterize the regular structure in 100 % of cases.Instead, it characterizes the probability of reflection occurrence.Practically, the use of the index is accompanied by setting a threshold.The events with the index below the threshold are rejected, the remaining events are treated as those containing a reflection.The higher the threshold chosen, the higher the probability, and thus fewer events will pass the threshold.Practically, the threshold is chosen from the comparison of the index with the visual investigation of an ensemble of events that is large enough to provide statistically significant results.
Figure 14 shows the probability distribution function (PDF) of the reflection index I R for the three categories of events.For the no-reflection category, the PDF has a strong maximum for events with the index below 1.For the reflection cases, the PDF has a tail for indexes below 3.At I R = 3, the PDFs for no-reflection and reflection cases have an equal magnitude.This allows for taking the value of 3 as a lowest threshold.About 5 % of events classified as clear reflection will be rejected by this threshold.At I R = 5, the PDF of noreflection cases reaches 0. This allows for taking the value of 5 as the highest (safe) threshold.About 10 % of events classified as clear reflection will be rejected by this threshold.

Conclusions
In this paper, we described our modification of the CT technique for the retrieval of bending angle profiles of reflected  rays.Our approach uses a combination of filtering in the impact parameters space with the standard GO retrieval.The filtering uses FIO in order to map the observed wave field to the impact parameter space.The field in the transformed space is multiplied with the filter function which suppresses the direct ray and let only pass both the not aliased and aliased components of the reflected ray.The filtered field is mapped back to the time domain.The phase of the resulting field is re-accumulated in the vicinity of the phase model of the reflected ray.We use the radio holographic spectra in order to estimate the reflection index and the expected error of the impact parameter.The reflection index characterizes the strength of the reflection.We validated our reflection in-dex definition by a comparison with the ROM SAF reflection flag database.In general, our reflection index indicates a good agreement with the database.Some discrepancies are partly explained by the misclassification of some events in the database, and partly by the random nature of RO signals, resulting in an overrated reflection index for some tropical events.These events are located on the distribution function.Based on this comparison, it is possible to estimate the threshold values of the reflection index.Its values exceeding 5 allow for a conclusion about a definite presence of reflection.Its values below 3 are typical for the absence of reflection.Values between 3 and 5 may correspond to different cases with a likely or unlikely reflection.The extracted profiles of reflected bending angles and impact parameters have the potential to be assimilated into NWP models through the forward operator in Eq. ( 17) (Aparicio et al., 2017).They might contribute anchoring the atmospheric conditions at the surface level.Studies are being conducted within the EU-METSAT ROM SAF to assess their added value and impact in NWP assimilation scenarios.
of the International Workshop on Occultations for Probing Atmosphere and Climate, Leibnitz, Austria, 8-14 September 2016.

Figure 2 .
Figure 2. Reflection over ice and snow.The branch of bending angle profile corresponding to the reflection looks like a nearly horizontal line at the impact height of about 2 km.Cf.Fig. 1.

Figure 3 .
Figure 3. Reflection over ice and snow in Antarctica.

Figure 2
Figure2through Fig.5show examples of reflections detected in COSMIC observations.Each Figure includes the event time, location, and the 2-D plot of the WDF for the observed wave field(Gorbunov et al., 2010) (cf.Fig.1).As illustrated by the figure, reflections can be observed over ocean or snow and ice.Many interesting examples are observed over Antarctica.SeeAparicio et al. (2017) for a statistical analysis of the reflection events.
18)where x (r) = rn M (r) is the refractive radius, the first term describing the refraction due to the refractivity gradient, and the second term describing the ray bending due to the reflection at the surface, p E = r E n M (r E ), and r E is the Earth's curvature radius with the account of the surface height above the geoid.Our neutral atmospheric model n M (r) is based on the MSIS-90 model complemented with 80 % relative humidity below 15 km as described inGorbunov et al. (2011).An example of reflected bending angle model ε MR (p) is shown in Fig.6.Together with the satellite orbit data, the model bending angle profile ε MR (p) allows for the determination of the excess phase for the reflected rays.To this end, we have to first numerically solve the following equation:

Figures 9 -
Figures 9-12 show examples of processing COSMIC data.The examples start with strong reflections, demonstrate further events with a decreasing reflection index, and lastly show an event with no reflection.These examples demonstrate that the reflection index can serve as a measure of reflected ray strength.In order to validate our retrieval algorithm and reflection index definition statistically, we performed a comparison of our retrievals with the ROM SAF reflection flag database(ROM SAF, 2016;Cardellach and Oliveras, 2016).The database contains occultation events classified into three categories: (1) no reflection, (2) reflection, and (3) unclear.The events are accompanied by the Supporting Vector Machine (SVM) index(Cardellach et al., 2009(Cardellach et al., , 2010) ) based on the radio holographic analysis and supervised learning method.

Figure 14 .
Figure 14.Probability distribution function of I R for the three categories of events.