the Creative Commons Attribution 4.0 License.

Special issue: Observing Atmosphere and Climate with Occultation Techniques...

**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

^{1,2},

^{3,4},

^{5}

**Michael E. Gorbunov et al.**Michael E. Gorbunov Estel Cardellach and Kent B. Lauritsen

^{1,2},

^{3,4},

^{5}

^{1}A. M. Obukhov Institute of Atmospheric Physics RAS, Moscow, Russia^{2}Hydrometeorological Research Centre of Russian Federation, Moscow, Russia^{3}Institut d'Estudis Espacials de Catalunya, Barcelona, Spain^{4}Instituto de Ciencias del Espacio (ICE, CSIC), Barcelona, Spain^{5}Danish Meteorological Institute, Copenhagen, Denmark

^{1}A. M. Obukhov Institute of Atmospheric Physics RAS, Moscow, Russia^{2}Hydrometeorological Research Centre of Russian Federation, Moscow, Russia^{3}Institut d'Estudis Espacials de Catalunya, Barcelona, Spain^{4}Instituto de Ciencias del Espacio (ICE, CSIC), Barcelona, Spain^{5}Danish Meteorological Institute, Copenhagen, Denmark

**Correspondence**: Michael E. Gorbunov (gorbunov@ifaran.ru)

**Correspondence**: Michael E. Gorbunov (gorbunov@ifaran.ru)

Received: 13 Jun 2017 – Discussion started: 23 Aug 2017 – Revised: 18 Dec 2017 – Accepted: 09 Jan 2018 – Published: 01 Mar 2018

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.

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., 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) (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, 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 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.

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:

where *A*(*t*) is the observed amplitude, *S*_{0}(*t*) is the
satellite-to-satellite distance evaluated from the orbit data, Δ*S*(*t*)
is the observed atmospheric excess phase. The smoothed excess phase
$\stackrel{\mathrm{\u203e}}{\mathrm{\Delta}S}\left(t\right)$ and its derivative $\stackrel{\mathrm{\u203e}}{\mathrm{\Delta}{S}^{\prime}}\left(t\right)$ are
obtained by the filtering of the measured excess phase Δ*S*(*t*). The
smoothed relative Doppler shift $\overline{d}\left(t\right)$ is obtained as follows:

where ${d}_{\mathrm{0}}\left(t\right)={{S}^{\prime}}_{\mathrm{0}}\left(t\right)/c$ is the vacuum relative Doppler
frequency shift. From $\overline{d}\left(t\right)$ and satellite orbit data we evaluate the
smooth impact parameter model $\overline{p}\left(t\right)$, bending angle model
$\overline{\mathit{\epsilon}}\left(t\right)$, and the derivative of the impact parameters over Doppler
shift $\mathrm{d}\phantom{\rule{0.125em}{0ex}}\overline{p}\left(t\right)/\mathrm{d}\phantom{\rule{0.125em}{0ex}}\overline{d}$ (Vorob'ev and Krasil'nikova, 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):

where *a*(p) is the amplitude function, whose definition can be
found in Gorbunov and Lauritsen (2004). The variable $\stackrel{\mathrm{\u0303}}{p}$ is the approximate
impact height (impact parameter with subtracted *R*_{E} due to
the definition of ${S}_{\mathrm{0}}^{\left(M\right)}\left(t\right)$). The transformed field
${\widehat{\mathrm{\Phi}}}_{\mathrm{2}}u\left(\stackrel{\mathrm{\u0303}}{p}\right)$ is represented as follows:

where ${A}^{\prime}\left(\stackrel{\mathrm{\u0303}}{p}\right)$ is the amplitude of the transformed field and ${\mathit{\phi}}^{\prime}\left(\stackrel{\mathrm{\u0303}}{p}\right)$ 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 ${\mathit{\phi}}^{\prime}\left(\stackrel{\mathrm{\u0303}}{p}\right)$. The amplitude function is evaluated using ${\mathit{\phi}}^{\prime}\left(\stackrel{\mathrm{\u0303}}{p}\right)$, which is the reason why field $u\left(\stackrel{\mathrm{\u0303}}{p}\right)$ 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 ${\stackrel{\mathrm{\u0303}}{p}}_{\mathrm{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 $\stackrel{\mathrm{\u0303}}{p}<{\stackrel{\mathrm{\u0303}}{p}}_{\mathrm{E}}$ and for $\stackrel{\mathrm{\u0303}}{p}>{\stackrel{\mathrm{\u0303}}{p}}_{\mathrm{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.

## 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 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
${\stackrel{\mathrm{\u0303}}{p}}^{\prime}$ 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:

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 ${{\widehat{\mathrm{\Phi}}}^{\prime}}_{\mathrm{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.

## 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
${\stackrel{\mathrm{\u0303}}{p}}_{\mathrm{E}}$. Therefore, we can form the following field
*u*_{R}(*t*) in the time domain that only contains the
reflections:

where ${\widehat{\mathrm{\Phi}}}_{\mathrm{2}}$ is the FIO, ${\widehat{\mathrm{\Phi}}}_{\mathrm{2}}^{-\mathrm{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.

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

## 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
*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 *x*(r)=*r**n*_{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 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:

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. (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 *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}_{\mathrm{E}}^{-\mathrm{1}}$ and located at an altitude of *h*_{SR}, we
arrive at the expression,

where $\mathrm{\Delta}p={p}_{\mathrm{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 $\mathrm{\Delta}r/{h}_{\text{SR}}\ll \mathrm{2}\sqrt{{h}_{\text{SR}}/\mathrm{\Delta}p}$ and therefore,
$\mathrm{d}{\mathit{\epsilon}}_{\mathrm{R}}\left(p\right)/\mathrm{d}p>\mathrm{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
${\mathit{u}}_{\text{Tx},\text{Rx}}^{\left(\mathrm{0}\right)}\left(t\right)$ corresponding to
satellite-to-satellite straight-line direction. An example of reflected
excess model phase is shown in Fig. 7.

## 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 (*i**k**S*_{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 $\stackrel{\mathrm{\u203e}}{{S}_{\mathrm{R}}}\left(t\right)$
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 ${\mathit{\omega}}_{\mathrm{0}}={\dot{S}}_{\text{MR}}\left({t}_{\mathrm{0}}\right)$. The spectrum can be
considered as function ${\stackrel{\mathrm{\u0303}}{u}}_{\mathrm{R}}\left(\mathrm{\Delta}p\right)$ of
relative impact height $\mathrm{\Delta}p=p-{p}_{\mathrm{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 ${\stackrel{\mathrm{\u0303}}{u}}_{max}$ is the maximum of the spectral density taken
within the interval of $\mathrm{\Delta}p\in \left[-\mathrm{0.1}\phantom{\rule{0.125em}{0ex}}\text{km},\phantom{\rule{0.125em}{0ex}}\mathrm{0.1}\phantom{\rule{0.125em}{0ex}}\text{km}\right]$, *p*_{max} is the location of the spectral maximum of the
reflection, ${\stackrel{\mathrm{\u0303}}{u}}_{\text{ave}}$ is the spectral density averaged
over the interval of $\left[{p}_{max}-\mathrm{0.3},{p}_{max}+\mathrm{0.3}\right]$,
${\stackrel{\mathrm{\u0303}}{u}}_{\text{bkg}}$ is the background (noise level) spectral density
estimated by averaging over the interval of $\mathrm{\Delta}p\in \left[\mathrm{1.0}\phantom{\rule{0.125em}{0ex}}\text{km},\phantom{\rule{0.125em}{0ex}}\mathrm{2.0}\phantom{\rule{0.125em}{0ex}}\text{km}\right]$, 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
${\stackrel{\mathrm{\u0303}}{u}}_{\text{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 ${\stackrel{\mathrm{\u0303}}{u}}_{max}$ significantly exceeds both
${\stackrel{\mathrm{\u0303}}{u}}_{\text{ave}}$ and $\mathit{\alpha}\phantom{\rule{0.125em}{0ex}}{\stackrel{\mathrm{\u0303}}{u}}_{\text{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.

## 3.3 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 (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:

where $\mathit{\chi}\left(\stackrel{\mathrm{\u0303}}{p}\right)$ is equal to unity inside the impact
height interval of $\left[{\stackrel{\mathrm{\u0303}}{p}}_{\mathrm{E}}-\mathrm{\Delta}{p}_{\mathrm{R}},{\stackrel{\mathrm{\u0303}}{p}}_{\mathrm{E}}\right]$ and the
corresponding aliased interval of $\left[{\stackrel{\mathrm{\u0303}}{p}}_{\mathrm{E}}+\mathrm{\Delta}{p}_{\text{alias}}-\mathrm{\Delta}{p}_{\mathrm{R}},{\stackrel{\mathrm{\u0303}}{p}}_{\mathrm{E}}+\mathrm{\Delta}{p}_{\text{alias}}\right]$. 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.

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

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.

A public version of the code used for this study is being prepared for the ROM SAF Radio Occultation Processing Package (ROPP); http://www.romsaf.org/ropp/.

The data used can be sent by request.

The authors declare that they have no conflict of interest.

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., https://doi.org/10.5194/amt-2017-241, 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, https://doi.org/10.1029/2001JD001402, 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, https://doi.org/10.1029/2004GL019775, 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: http://www.romsaf.org/general-documents/rsr/rsr_23.pdf, 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: http://www.romsaf.org/Publications/reports/romsaf_vs27_rep_v10.pdf (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, https://doi.org/10.1029/2000RS002592, 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, https://doi.org/10.1029/2001JD000889, 2002b. a

Gorbunov, M. E.: Radioholographic analysis of radio occultation data in multipath zones, Radio Sci., 37, 14–1–14–9, https://doi.org/10.1029/2000RS002577, 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, https://doi.org/10.1029/2003RS002971, 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, https://doi.org/10.1029/2005JD006427, 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, https://doi.org/10.1029/2010RS004388, 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, https://doi.org/10.1175/2011JTECHA1489.1, 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: http://irowg.org/wpcms/wp-content/uploads/2013/12/lauritsen.pdf (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, https://doi.org/10.1029/2002RS002763, 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, https://doi.org/10.1029/2003RS002899, 2004. a, b

ROM SAF: ROM SAF Reflection Flag database, available at: http://www.romsaf.org/pub/demo/reflection_flag/ (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