Articles | Volume 11, issue 2
Research article
01 Mar 2018
Research article |  | 01 Mar 2018

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

Michael E. Gorbunov, Estel Cardellach, and Kent B. Lauritsen

Linear and non-linear representations of wave fields constitute the basis of modern algorithms for analysis of radio occultation (RO) data. Linear representations are implemented by Fourier Integral Operators, which allow for high-resolution retrieval of bending angles. Non-linear representations include Wigner Distribution Function (WDF), which equals the pseudo-density of energy in the ray space. Representations allow for filtering wave fields by suppressing some areas of the ray space and mapping the field back from the transformed space to the initial one. We apply this technique to the retrieval of reflected rays from RO observations. The use of reflected rays may increase the accuracy of the retrieval of the atmospheric refractivity. 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 reflection index that characterizes the relative intensity of the reflected ray with respect to the direct ray. The index allows for the automatic identification of events with reflections. We use the radio holographic estimate of the errors of the retrieved bending angle profiles of reflected rays. A comparison of indices evaluated for a large base of events including the visual identification of reflections indicated a good agreement with our definition of reflection index.

1 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 Hocke2001; Beyerle et al.2002). Similar patterns were also found in Microlab-1 GPS/MET data (Gorbunov2002b, 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 (Gorbunov2002c; Cardellach et al.2009, 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, 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) (Gorbunov2002a; Gorbunov and Lauritsen2004), 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, 2010). For this projection we use the inverse FIO. Generally, our approach is similar to that developed in Cardellach et al. (2009, 2010), where the spectrograms were filtered in order to remove direct rays and inverted. However, the CT-based 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 Hocke2001; Beyerle et al.2002), as it is more accurate and computationally more efficient (Gorbunov2002c). This work was published as a technical report (Gorbunov2016).

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.

2 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.

2.1 Reflected rays retrieval in impact parameter domain

The measured wave field has the following form:

(1) u ( t ) = A ( t ) exp i k S 0 ( t ) + Δ S ( t ) ,

where A(t) is the observed amplitude, S0(t) is the satellite-to-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:

(2) d ¯ ( t ) = d 0 ( t ) - Δ S ( t ) c ,

where d0(t)=S0(t)/c is the vacuum relative Doppler frequency shift. From d¯(t) and satellite orbit data we evaluate the smooth impact parameter model p¯(t), bending angle model ε¯(t), and the derivative of the impact parameters over Doppler shift dp¯(t)/dd¯ (Vorob'ev and Krasil'nikova1994). The ancillary function f(t) is evaluated as follows (Gorbunov and Lauritsen2004):

(3) f ( t ) = p ¯ ( t ) - d ¯ ( t ) d p ¯ ( t ) d d ¯ .

The new coordinate is determined as follows (Gorbunov and Lauritsen2004):

(4) Υ ( t ) = Υ 0 - c t 0 t d p ¯ ( t ) d d ¯ - 1 d t ,

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):

(5) f I ( t ) = Υ t 0 Υ ( t ) f t d Υ t .

Using this function, we evaluate the vacuum and observed phase path with subtracted model as follows:


where RE is the Earth's local curvature radius. Subtraction of REΥ(t), a linear function of Υ, from the phase corresponds to the reduction of the frequency, which equals the impact parameter, by a constant of RE. All the functions of time t can also be looked at as functions of the new coordinate Υ.

Figure 1Reflection and multipath propagation geometry. (a) Upper panel shows the radio occultation geometry. At moment of time t1, one ray is observed, at moment of time t2 the receiver is in a multipath area and observes 3 direct rays and one reflected ray. (b) Each moment of time corresponds to a line in the impact parameter – bending angle coordinate plane; the line indicating possible values of impact parameters and bending angles that can be observed by the receiver at this moment of time. pE is the impact parameter of the ray that touches the Earth's surface. Reflected rays have impact parameters p<pE. The bending angle profile of reflected rays is shown as a dashed line. The lines shown in red represent the constant modified impact parameter.


The FIO is defined as follows (Gorbunov and Lauritsen2004):


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 RE due to the definition of S0M(t)). The transformed field Φ^2up̃ is represented as follows:

(9) Φ ^ 2 u p ̃ = A p ̃ exp i φ p ̃ ,

where Ap̃ 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 up̃ 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 (Gorbunov2002a; Gorbunov and Lauritsen2004; 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:

(10) d φ ¯ D p ̃ d p ̃ , d φ ¯ R p ̃ d p ̃ ,

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.

2.2 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 (6) corresponds to the following linear canonical transform (Gorbunov and Lauritsen2004):


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:

(13) p ̃ = p ̃ + β Υ ,

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:

(14) d ε d p ̃ = d ε d p ̃ 1 + β d ε d p ̃ .

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.

Figure 2Reflection 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.


2.3 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 non-monotonic. 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 uR(t) in the time domain that only contains the reflections:

(17) u R ( t ) = Φ ^ 2 - 1 Φ ^ 2 u ( t ) θ p ̃ E - p ̃ ,

where Φ^2 is the FIO, Φ^2-1 is its inverse, and θ is the theta-function, 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.

Figure 3Reflection over ice and snow in Antarctica.


Figure 4Reflection over ocean.


Figure 5Reflection over ocean.


Figure 6Reflected bending angle model for occultation event 1 January 2008, 01:02:23 UTC, 70.28 N 121.87 W (the same event as in Fig. 2).


2.4 Examples of reflections in COSMIC observations and discussion

Figure 2 through Fig. 5 show 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. See Aparicio et al. (2017) for a statistical analysis of the reflection events.

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 Kornblueh2001; Gorbunov and Lauritsen2009). 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.

3 Reflection retrieval implementation

3.1 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 nM(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):

(18) ε MR ( p ) = - 2 p p E d ln n M d x d x x 2 - p 2 - 2 arccos p p E ,

where x(r)=rnM(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, pE=rEnM(rE), and rE is the Earth's curvature radius with the account of the surface height above the geoid. Our neutral atmospheric model nM(r) is based on the MSIS-90 model complemented with 80 % relative humidity below 15 km as described in Gorbunov 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:

(19) θ ( t ) = ε MR ( p ) + arccos p r Tx ( t ) + arccos p r Rx ( t ) ,

where θ(t) is the satellite-to-satellite angle with respect to the local curvature center, rTx,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, pMR(t). Dependence pMR(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. (11), 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 rE 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 -rE-1 and located at an altitude of hSR, we arrive at the expression,


where Δp=pE-p. Assuming that hSR is about the PBL height (i.e. 1.5 km), and the superrefraction layer thickness is about 0.2 km, and Δp<0.2km, we see that Δr/hSR2hSR/Δ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 xTx,Rx(t), the ray directions at the satellites, unit vectors uTx,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 VRx,Tx(t), we find the relative Doppler frequency shift dMR(t):

(24) V Tx ( t ) u Tx ( t ) - V Rx ( t ) u Rx ( t ) = c d MR ( t ) .

The excess phase is obtained by integrating the Doppler shift:

(25) S MR ( t ) = c d ( 0 ) ( t ) - d MR ( t ) d t ,

where d(0)(t) is the vacuum Doppler shift for the direct rays, evaluated from Eq. (24), by inserting unit vector uTx,Rx(0)(t) corresponding to satellite-to-satellite straight-line direction. An example of reflected excess model phase is shown in Fig. 7.

Figure 7Reflected excess phase model for occultation event 1 January 2008, 01:02:23 UTC, 70.28 N 121.87 W (the same event as in Fig. 2).


Figure 8Radio holographic spectrum amplitude ũRω for occultation event 1 January 2008, 01:02:23 UTC, 70.28 N 121.87 W.


Figure 9Occultation event 1 January 2008, 14:59:50 UTC, 72.34 N 164.30 W. Expressed as a: (a) Wigner distribution function; (b) radio holographic spectrum; (c) model and retrieved reflection bending angles, the latter with error bars indicating the uncertainty estimate. Reflection index 27.2.


Figure 10Occultation event 1 January 2008, 08:46:37 UTC, 61.38 S 69.50 E. Expressed as a: (a) Wigner distribution function; (b) radio holographic spectrum; (c) model and retrieved reflection bending angles, the latter with error bars indicating the uncertainty estimate. Reflection index 20.8.


Figure 11Occultation event 1 January 2008, 10:07:11 UTC, 48.68 S 134.46 W. Expressed as a: (a) Wigner distribution function; (b) radio holographic spectrum; (c) model and retrieved reflection bending angles, the latter with error bars indicating the uncertainty estimate. Reflection index 13.6.


Figure 12Occultation event 1 January 2008, 04:39:47 UTC, 44.78 N 44.75 W. Expressed as a: (a) Wigner distribution function; (b) radio holographic spectrum; (c) model and retrieved reflection bending angles, the latter with error bars indicating the uncertainty estimate. Reflection index 8.1.


Figure 13Occultation event 1 January 2008, 02:51:36 UTC, 49.16 N 76.55 W. Expressed as a: (a) Wigner distribution function; (b) radio holographic spectrum; (c) model and retrieved reflection bending angles, the latter with error bars indicating the uncertainty estimate. Reflection index 0.07.


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


3.2 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 (ikSMR(t)) as the reference signal and evaluate the radio holographic spectrum as follows:

(26) u ̃ R ω = A ( t ) exp i k S ( t ) - S MR ( t ) - i ω t d t ,

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 SR(t) rather than the model SMR(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 t0 of the time interval. Moreover, it is convenient to introduce the reference value p0 of the impact parameter, corresponding to frequency ω0=S˙MRt0. The spectrum can be considered as function ũRΔp of relative impact height Δp=p-p0.

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 IR as a functional of the radio holographic spectrum, depending on a series of empirical parameters:


where ũmax is the maximum of the spectral density taken within the interval of Δp-0.1km,0.1km, pmax is the location of the spectral maximum of the reflection, ũave is the spectral density averaged over the interval of pmax-0.3,pmax+0.3, ũbkg is the background (noise level) spectral density estimated by averaging over the interval of Δp1.0km,2.0km, and α is the regularization parameter, pM(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 ũ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 ũmax significantly exceeds both ũave and αũbkg. The optimal value of α was empirically estimated to be about 0.2. The additional exponential factor in the definition of IR 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 IR=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.

3.3 Filtering in impact parameter space

In order to extract the reflected field component uR(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 Δpalias between non-aliased and aliased components for typical observation geometry is about 8–10 km. The exact value Δpalias for a specific event is evaluated by using the GO relationship (13) and (13). To this end, we evaluate impact height from the original relative Doppler shift, and from the relative Doppler shift corresponding to the aliased frequency shifted by the sampling rate. In order to retain the aliased component, we modify filter (10) as follows:

(31) u R ( t ) = Φ ^ 2 - 1 Φ ^ 2 u ( t ) χ p ̃ ,

where χp̃ is equal to unity inside the impact height interval of p̃E-ΔpR,p̃E and the corresponding aliased interval of p̃E+Δpalias-ΔpR,p̃E+Δpalias. The width ΔpR of these intervals is set to 1 km. Outside these intervals, we employ the Gaussian apodization with a characteristic width of δpR=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.

4 Examples of processing COSMIC data

Figures 912 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 SAF2016; Cardellach and Oliveras2016). 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, 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 IR 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 IR=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 IR=5, the PDF of no-reflection 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.

5 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 index 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 EUMETSAT ROM SAF to assess their added value and impact in NWP assimilation scenarios.

Code availability

A public version of the code used for this study is being prepared for the ROM SAF Radio Occultation Processing Package (ROPP);

Data availability

The data used can be sent by request.

Competing interests

The authors declare that they have no conflict of interest.

Special issue statement

This article is part of the special issue “Observing Atmosphere and Climate with Occultation Techniques – Results from the OPAC-IROWG 2016 Workshop”. It is a result of the International Workshop on Occultations for Probing Atmosphere and Climate, Leibnitz, Austria, 8–14 September 2016.


Part of this work (Sect. 2 and part of Sect. 3) was conducted as part of the Visiting Scientist program of the Radio Occultation Meteorology Satellite Applications Facility (ROM SAF), which is a decentralized operational radio occultation processing center under EUMETSAT. Michael E. Gorbunov was a ROM SAF Visiting Scientist for the CDOP-2 VS27 project and Estel Cardellach and Kent B. Lauritsen are members of the ROM SAF. The work on the rest of Sect. 3 and Sect. 4 was co-supported by Russian Foundation for Basic Research (grant RFBR No. 16-05-00358 A).

Edited by: Jens Wickert
Reviewed by: two anonymous referees


Aparicio, J. M., Cardellach, E., and Rodriguez, H.: Information content in reflected signals during GPS radio occultation observations, Atmos. Meas. Tech. Discuss.,, in review, 2017. a, b, c

Beyerle, G. and Hocke, K.: Observation and simulation of direct and reflected GPS signals in radio occultation experiments, Geophys. Res. Lett., 28, 1895–1898, 2001. a, b

Beyerle, G., Hocke, K., Wickert, J., Schmidt, T., Marquardt, C., and Reigber, C.: GPS radio occultations with CHAMP: a radio holographic analysis of GPS signal propagation in the troposphere and surface reflections, J. Geophys. Res., 107, 27–1–27–14,, 2002. a, b

Cardellach, E., Ao, C. O., de la Torre-Juárez, M., and Hajj, G. A.: Carrier phase delay altimetry with GPS-reflection/occultation interferometry from low Earth orbiters, Geophys. Res. Lett., 31, L10402,, 2004. a

Cardellach, E., Oliveras, S., and Rius, A.: GNSS signal interference classified by means of a supervised learning method applied in the time-frequency domain, in: IEEE Proceedings of 2009 2nd International Congress on Image and Signal Processing, ISBN 13: 978-1-4244-4130-3, Institute of Electrical and Electronic Engineers, Tijuan, China, 17–21 October 2009, 2009. a, b, c, d

Cardellach, E., Oliveras, S., and Rius, A.: Tropospheric information content embedded in GNSS RO reflected signals, in: The International Beacon Satellite Symposium, 7–11 June 2010, CIMNE, Barcelona, BSS2010, 2010. a, b, c, d

Cardellach, E. and Oliveras, S.: Assessment of a potential reflection flag product, SAF/ROM/METO/REP/RSR/023, IEEC, Barcelona, Spain, available at:, last access: 23 January 2016. a

Gorbunov, M.: Development of wave optics code for the retrieval of bending angle profiles for reflected rays, ROM SAF CDOP-2, Visiting Scientist Report 27, Danish Meteorological Institute, European Centre for Medium-Range Weather Forecasts, Institut d'Estudis Espacials de Catalunya, Met Office, available at: (last access: 20 February 2018), (SAF/ROM/DMI/REP/VS27/001), 2016.  a

Gorbunov, M. E.: Canonical transform method for processing radio occultation data in the lower troposphere, Radio Sci., 37, 9–1–9–10,, 2002a. a, b

Gorbunov, M. E.: Radio-holographic analysis of Microlab-1 radio occultation data in the lower troposphere, J. Geophys. Res.-Atmos., 107, 7–1–7–10,, 2002b. a

Gorbunov, M. E.: Radioholographic analysis of radio occultation data in multipath zones, Radio Sci., 37, 14–1–14–9,, 2002c. a, b, c

Gorbunov, M. E. and Kornblueh, L.: Analysis and validation of GPS/MET radio occultation data, J. Geophys. Res., 106, 17161–17169, 2001. a

Gorbunov, M. E. and Lauritsen, K. B.: Analysis of wave fields by Fourier integral operators and its application for radio occultations, Radio Sci., 39, RS4010,, 2004. a, b, c, d, e, f, g

Gorbunov, M. E. and Lauritsen, K. B.: Error estimate of bending angles in the presence of strong horizontal gradients, in: New Horizons in Occultation Research, edited by: Steiner, A., Pirscher, B., Foelsche, U., and Kirchengast, G., 17–26, Springer, Berlin, Heidelberg, 2009. a

Gorbunov, M. E., Lauritsen, K. B., Rhodin, A., Tomassini, M., and Kornblueh, L.: Radio holographic filtering, error estimation, and quality control of radio occultation data, J. Geophys. Res., 111, D10105,, 2006. a

Gorbunov, M. E., Lauritsen, K. B., and Leroy, S. S.: Application of Wigner distribution function for analysis of radio occultations, Radio Sci., 45, RS6011,, 2010. a, b

Gorbunov, M. E., Shmakov, A. V., Leroy, S. S., and Lauritsen, K. B.: COSMIC radio occultation processing: Cross-center comparison and validation, J. Atmos. Ocean. Tech., 28, 737–751,, 2011. a

Gorbunov, M. E., Lauritsen, K. B., and Leroy, S. S.: Analysis of RO data retrieved from the Wigner distribution function, in: International Radio Occultation Working Group, 2nd Workshop, 28 March–3 April 2012, Stanley Hotel in Estes Park, CO, USA, available at: (20 February 2018), 2012. a

Jensen, A. S., Lohmann, M. S., Benzon, H.-H., and Nielsen, A. S.: Full spectrum inversion of radio occultation signals, Radio Sci., 38, 6–1–6–15,, 2003. a

Jensen, A. S., Lohmann, M. S., Nielsen, A. S., and Benzon, H.-H.: Geometrical optics phase matching of radio occultation signals, Radio Sci., 39, RS3009,, 2004. a, b

ROM SAF: ROM SAF Reflection Flag database, available at: (20 February 2018), 2016. a

Vorob'ev, V. V. and Krasil'nikova, T. G.: Estimation of the accuracy of the atmospheric refractive index recovery from Doppler shift measurements at frequencies used in the NAVSTAR system, Izv. Atmos. Ocean. Phy., 29, 602–609, 1994. a

Short summary
We apply linear and non-linear representations of wave fields, based on Fourier integral operators and Wigner distribution function, to the retrieval of reflected rays from radio occultation observations. We introduce a reflection index that characterizes the relative intensity of the reflected ray. A comparison of indices evaluated for a large base of events including the visual identification of reflections indicated a good agreement with our definition of reflection index.