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

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:

$$\begin{array}{}\text{(1)}& {\displaystyle }{\displaystyle }u\left(t\right)=A\left(t\right)\exp \left(ik\left(S_{\mathrm{0}}\left(t\right)+\mathrm{\Delta }S\left(t\right)\right)\right),\end{array}$$

where A(t) is the observed amplitude, S₀(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{‾}}{\mathrm{\Delta }S}\left(t\right)$ and its derivative $\stackrel{\mathrm{‾}}{\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:

$$\begin{array}{}\text{(2)}& {\displaystyle }{\displaystyle }\overline{d}\left(t\right)=d_{\mathrm{0}}\left(t\right)-{\displaystyle \frac{\stackrel{\mathrm{‾}}{\mathrm{\Delta }{S}^{\prime }}\left(t\right)}{c}},\end{array}$$

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}\overline{p}\left(t\right)/\mathrm{d}\overline{d}$ (Vorob'ev and Krasil'nikova, 1994). The ancillary function f(t) is evaluated as follows (Gorbunov and Lauritsen, 2004):

$$\begin{array}{}\text{(3)}& {\displaystyle }{\displaystyle }f\left(t\right)=\overline{p}\left(t\right)-\overline{d}\left(t\right){\displaystyle \frac{\mathrm{d}\overline{p}\left(t\right)}{\mathrm{d}\overline{d}}}.\end{array}$$

The new coordinate is determined as follows (Gorbunov and Lauritsen, 2004):

$$\begin{array}{}\text{(4)}& {\displaystyle }{\displaystyle }\mathrm{{\rm Y}}\left(t\right)={\mathrm{{\rm Y}}}_{\mathrm{0}}-c\underset{t_{\mathrm{0}}}{\overset{t}{\int }}{\left({\displaystyle \frac{\mathrm{d}\overline{p}\left(t\right)}{\mathrm{d}\overline{d}}}\right)}^{-\mathrm{1}}\mathrm{d}{t}^{\prime },\end{array}$$

where Υ₀ is a constant determined in such a way that Υ(t)≥0 for the observation time interval. We evaluate the integral of f(t):

$$\begin{array}{}\text{(5)}& {\displaystyle }{\displaystyle }{f}_{\mathrm{I}}\left(t\right)=\underset{\mathrm{{\rm Y}}\left(t_{\mathrm{0}}\right)}{\overset{\mathrm{{\rm Y}}\left(t\right)}{\int }}f\left(t^{\prime }\right)\mathrm{d}\mathrm{{\rm Y}}\left(t^{\prime }\right).\end{array}$$

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

$$\begin{array}{}\text{(6)}& {\displaystyle }& {\displaystyle }{S}_{\mathrm{0}}^{\left(M\right)}\left(t\right)=S_{\mathrm{0}}\left(t\right)-R_{\mathrm{E}}\mathrm{{\rm Y}}\left(t\right)+f_{\mathrm{I}}\left(t\right),\text{(7)}& {\displaystyle }& {\displaystyle }{S}^{\left(M\right)}\left(t\right)=S_{\mathrm{0}}^{\left(M\right)}\left(t\right)+\mathrm{\Delta }S\left(t\right),\end{array}$$

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

Figure 1 Reflection and multipath propagation geometry. (a) Upper panel shows the radio occultation geometry. At moment of time t₁, one ray is observed, at moment of time t₂ 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. p_E is the impact parameter of the ray that touches the Earth's surface. Reflected rays have impact parameters p<p_E. 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 Lauritsen, 2004):

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\widehat{\mathrm{\Phi }}}_{\mathrm{2}}u\left(\stackrel{\mathrm{̃}}{p}\right)=\\ \text{(8)}& {\displaystyle }& {\displaystyle }\sqrt{{\displaystyle \frac{k}{\mathrm{2}\mathit{\pi }}}}a\left(\stackrel{\mathrm{̃}}{p}\right)\int A\left(\mathrm{{\rm Y}}\right)\exp \left(ik{S}^{\left(M\right)}\left(\mathrm{{\rm Y}}\right)-ik\stackrel{\mathrm{̃}}{p}\mathrm{{\rm Y}}\right)\mathrm{d}\mathrm{{\rm Y}},\end{array}$$

where a(p) is the amplitude function, whose definition can be found in Gorbunov and Lauritsen (2004). The variable $\stackrel{\mathrm{̃}}{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{̃}}{p}\right)$ is represented as follows:

$$\begin{array}{}\text{(9)}& {\displaystyle }{\displaystyle }{\widehat{\mathrm{\Phi }}}_{\mathrm{2}}u\left(\stackrel{\mathrm{̃}}{p}\right)=A^{\prime }\left(\stackrel{\mathrm{̃}}{p}\right)\exp \left(i{\mathit{\phi }}^{\prime }\left(\stackrel{\mathrm{̃}}{p}\right)\right),\end{array}$$

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

$$\begin{array}{}\text{(10)}& {\displaystyle }{\displaystyle \frac{\mathrm{d}{\overline{{\mathit{\phi }}^{\prime }}}_{\mathrm{D}}\left(\stackrel{\mathrm{̃}}{p}\right)}{\mathrm{d}\stackrel{\mathrm{̃}}{p}}},{\displaystyle \frac{\mathrm{d}{\overline{{\mathit{\phi }}^{\prime }}}_{\mathrm{R}}\left(\stackrel{\mathrm{̃}}{p}\right)}{\mathrm{d}\stackrel{\mathrm{̃}}{p}}},\end{array}$$

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

$$\begin{array}{}\text{(11)}& {\displaystyle }& {\displaystyle }\stackrel{\mathrm{̃}}{p}=f\left(\mathrm{{\rm Y}}\right)+\mathit{\eta },\text{(12)}& {\displaystyle }& {\displaystyle }\mathit{\xi }=-\mathrm{{\rm Y}},\end{array}$$

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:

$$\begin{array}{}\text{(13)}& {\displaystyle }{\displaystyle }{\stackrel{\mathrm{̃}}{p}}^{\prime }=\stackrel{\mathrm{̃}}{p}+\mathit{\beta }\mathrm{{\rm Y}},\end{array}$$

where β is a tunable parameter. Lines of constant coordinate ${\stackrel{\mathrm{̃}}{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:

$$\begin{array}{}\text{(14)}& {\displaystyle }{\displaystyle \frac{\mathrm{d}\mathit{\epsilon }}{\mathrm{d}{\stackrel{\mathrm{̃}}{p}}^{\prime }}}={\displaystyle \frac{\frac{\mathrm{d}\mathit{\epsilon }}{\mathrm{d}\stackrel{\mathrm{̃}}{p}}}{\mathrm{1}+\mathit{\beta }\frac{\mathrm{d}\mathit{\epsilon }}{\mathrm{d}{\stackrel{\mathrm{̃}}{p}}^{\prime }}}}.\end{array}$$

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:

$$\begin{array}{}\text{(15)}& {\displaystyle }& {\displaystyle }{\stackrel{\mathrm{̃}}{p}}^{\prime }=f\left(\mathrm{{\rm Y}}\right)+\mathit{\beta }\mathrm{{\rm Y}}+\mathit{\eta }\equiv f^{\prime }\left(\mathrm{{\rm Y}}\right)+\mathit{\eta },\text{(16)}& {\displaystyle }& {\displaystyle }\mathit{\xi }=-\mathrm{{\rm Y}}.\end{array}$$

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.

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.

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

$$\begin{array}{}\text{(17)}& {\displaystyle }{\displaystyle }{u}_{\mathrm{R}}\left(t\right)={\widehat{\mathrm{\Phi }}}_{\mathrm{2}}^{-\mathrm{1}}\left[{\widehat{\mathrm{\Phi }}}_{\mathrm{2}}\left[u\left(t\right)\right]\mathit{\theta }\left({\stackrel{\mathrm{̃}}{p}}_{\mathrm{E}}-\stackrel{\mathrm{̃}}{p}\right)\right],\end{array}$$

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.

Figure 3 Reflection over ice and snow in Antarctica.

Figure 4 Reflection over ocean.

Figure 5 Reflection over ocean.

Figure 6 Reflected 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 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 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 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):

$$\begin{array}{}\text{(18)}& {\displaystyle }{\displaystyle }{\mathit{\epsilon }}_{\text{MR}}\left(p\right)=-\mathrm{2}p\underset{p_{\mathrm{E}}}{\overset{\mathrm{\infty }}{\int }}{\displaystyle \frac{\mathrm{d}\ln n_{\mathrm{M}}}{\mathrm{d}x}}{\displaystyle \frac{\mathrm{d}x}{\sqrt{x^{\mathrm{2}}-p^{\mathrm{2}}}}}-\mathrm{2}\arccos \left({\displaystyle \frac{p}{p_{\mathrm{E}}}}\right),\end{array}$$

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

$$\begin{array}{}\text{(19)}& {\displaystyle }{\displaystyle }\mathit{\theta }\left(t\right)={\mathit{\epsilon }}_{\text{MR}}\left(p\right)+\arccos {\displaystyle \frac{p}{r_{\text{Tx}}\left(t\right)}}+\arccos {\displaystyle \frac{p}{r_{\text{Rx}}\left(t\right)}},\end{array}$$

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,

$$\begin{array}{ll}\text{(20)}& {\displaystyle }{\mathit{\epsilon }}_{\mathrm{R}}\left(p\right)& {\displaystyle }=-\mathrm{2}p\underset{p_{\mathrm{E}}}{\overset{\mathrm{\infty }}{\int }}{\displaystyle \frac{\mathrm{d}\ln n}{\mathrm{d}x}}{\displaystyle \frac{\mathrm{d}x}{\sqrt{x^{\mathrm{2}}-p^{\mathrm{2}}}}}-\mathrm{2}\arccos \left({\displaystyle \frac{p}{p_{\mathrm{E}}}}\right){\displaystyle }& {\displaystyle }\approx -\sqrt{\mathrm{2}{r}_{\mathrm{E}}}\underset{r_{\mathrm{E}}}{\overset{\mathrm{\infty }}{\int }}{\displaystyle \frac{\mathrm{d}\ln n\left(r\right)}{\mathrm{d}r}}{\displaystyle \frac{\mathrm{d}r}{\sqrt{rn\left(r\right)-p}}}-\mathrm{2}\arccos \left({\displaystyle \frac{p}{p_{\mathrm{E}}}}\right),\end{array}$$

where r_E is the Earth's radius. Then we can derive the expression for the derivative of the bending angle as:

$$\begin{array}{ll}{\displaystyle \frac{\mathrm{d}{\mathit{\epsilon }}_{\mathrm{R}}\left(p\right)}{\mathrm{d}p}}\approx & {\displaystyle }\sqrt{{\displaystyle \frac{r_{\mathrm{E}}}{\mathrm{2}}}}\underset{r_{\mathrm{E}}}{\overset{\mathrm{\infty }}{\int }}{\displaystyle \frac{\mathrm{d}\ln n\left(r\right)}{\mathrm{d}r}}{\displaystyle \frac{\mathrm{d}r}{{\left(rn\left(r\right)-p\right)}^{\mathrm{3}/\mathrm{2}}}}\\ \text{(21)}& {\displaystyle }& {\displaystyle }+\sqrt{{\displaystyle \frac{\mathrm{2}}{r_{\mathrm{E}}}}}{\displaystyle \frac{\mathrm{1}}{\sqrt{p_{\mathrm{E}}-p}}}.\end{array}$$

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,

$$\begin{array}{ll}{\displaystyle \frac{\mathrm{d}{\mathit{\epsilon }}_{\mathrm{R}}\left(p\right)}{\mathrm{d}p}}& {\displaystyle }\approx -\sqrt{{\displaystyle \frac{\mathrm{1}}{\mathrm{2}{r}_{\mathrm{E}}}}}{\displaystyle \frac{\mathrm{\Delta }r}{h_{\text{SR}}^{\mathrm{3}/\mathrm{2}}}}+\sqrt{{\displaystyle \frac{\mathrm{2}}{r_{\mathrm{E}}}}}{\displaystyle \frac{\mathrm{1}}{\mathrm{\Delta }{p}^{\mathrm{1}/\mathrm{2}}}}\\ \text{(22)}& {\displaystyle }& {\displaystyle }=\sqrt{{\displaystyle \frac{\mathrm{1}}{\mathrm{2}{r}_{\mathrm{E}}{h}_{\text{SR}}}}}\left(-{\displaystyle \frac{\mathrm{\Delta }r}{h_{\text{SR}}}}+\mathrm{2}\sqrt{{\displaystyle \frac{h_{\text{SR}}}{\mathrm{\Delta }p}}}\right),\end{array}$$

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:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\mathit{x}}_{\text{Rx}}\left(t\right)\times {\mathit{u}}_{\text{Rx}}\left(t\right)={\mathit{x}}_{\text{Tx}}\left(t\right)\times {\mathit{u}}_{\text{Tx}}\left(t\right),\\ \text{(23)}& {\displaystyle }& {\displaystyle }\left|{\mathit{x}}_{\text{Rx}}\left(t\right)\times {\mathit{u}}_{\text{Rx}}\left(t\right)\right|=\left|{\mathit{x}}_{\text{Tx}}\left(t\right)\times {\mathit{u}}_{\text{Tx}}\left(t\right)\right|=p,\end{array}$$

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

$$\begin{array}{}\text{(24)}& {\displaystyle }{\displaystyle }{\mathit{V}}_{\text{Tx}}\left(t\right)\cdot {\mathit{u}}_{\text{Tx}}\left(t\right)-{\mathit{V}}_{\text{Rx}}\left(t\right)\cdot {\mathit{u}}_{\text{Rx}}\left(t\right)=c{d}_{\text{MR}}\left(t\right).\end{array}$$

The excess phase is obtained by integrating the Doppler shift:

$$\begin{array}{}\text{(25)}& {\displaystyle }{\displaystyle }{S}_{\text{MR}}\left(t\right)=c\int \left(d^{\left(\mathrm{0}\right)}\left(t\right)-d_{\text{MR}}\left(t\right)\right)\mathrm{d}t,\end{array}$$

where d⁽⁰⁾(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.

Figure 7 Reflected 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 8 Radio holographic spectrum amplitude $\left|{\stackrel{\mathrm{̃}}{u}}_{\mathrm{R}}\left(\mathit{\omega }\right)\right|$ for occultation event 1 January 2008, 01:02:23 UTC, 70.28^∘ N 121.87^∘ W.

Figure 9 Occultation 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 10 Occultation 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 11 Occultation 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 12 Occultation 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 13 Occultation 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 14 Probability distribution function of I_R 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 (ikS_MR(t)) as the reference signal and evaluate the radio holographic spectrum as follows:

$$\begin{array}{}\text{(26)}& {\displaystyle }{\displaystyle }{\stackrel{\mathrm{̃}}{u}}_{\mathrm{R}}\left(\mathit{\omega }\right)=\int A\left(t\right)\exp \left(ik\left[S\left(t\right)-S_{\text{MR}}\left(t\right)\right]-i\mathit{\omega }t\right)\mathrm{d}t,\end{array}$$

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{‾}}{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₀ of the time interval. Moreover, it is convenient to introduce the reference value p₀ 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{̃}}{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:

$$\begin{array}{}\text{(27)}& {\displaystyle }& {\displaystyle }{\stackrel{\mathrm{̃}}{u}}_{max}=\underset{\left[-\mathrm{0.1},\mathrm{0.1}\right]}{max}{\left|\stackrel{\mathrm{̃}}{u}\left(\mathrm{\Delta }p\right)\right|}^{\mathrm{2}},\text{(28)}& {\displaystyle }& {\displaystyle }{\stackrel{\mathrm{̃}}{u}}_{\text{ave}}={〈{\left|\stackrel{\mathrm{̃}}{u}\left(\mathrm{\Delta }p\right)\right|}^{\mathrm{2}}〉}_{\left[p_{max}-\mathrm{0.3},p_{max}+\mathrm{0.3}\right]},\text{(29)}& {\displaystyle }& {\displaystyle }{\stackrel{\mathrm{̃}}{u}}_{\text{bkg}}={〈{\left|\stackrel{\mathrm{̃}}{u}\left(\mathrm{\Delta }p\right)\right|}^{\mathrm{2}}〉}_{\left[\mathrm{1.0},\mathrm{2.0}\right]},\text{(30)}& {\displaystyle }& {\displaystyle }{I}_{\mathrm{R}}={\displaystyle \frac{{\stackrel{\mathrm{̃}}{u}}_{max}^{\mathrm{2}}}{{\stackrel{\mathrm{̃}}{u}}_{\text{ave}}\left({\stackrel{\mathrm{̃}}{u}}_{max}+\mathit{\alpha }{\stackrel{\mathrm{̃}}{u}}_{\text{bkg}}\right)}}〈\exp \left(-{\left[{\displaystyle \frac{p\left(t\right)-p_{\mathrm{M}}\left(t\right)}{\mathrm{2}\mathit{\delta }p\left(t\right)}}\right]}^{\mathrm{2}}\right)〉,\end{array}$$

where ${\stackrel{\mathrm{̃}}{u}}_{max}$ is the maximum of the spectral density taken within the interval of $\mathrm{\Delta }p\in \left[-\mathrm{0.1}\text{km},\mathrm{0.1}\text{km}\right]$, p_max is the location of the spectral maximum of the reflection, ${\stackrel{\mathrm{̃}}{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{̃}}{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}\text{km},\mathrm{2.0}\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{̃}}{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{̃}}{u}}_{max}$ significantly exceeds both ${\stackrel{\mathrm{̃}}{u}}_{\text{ave}}$ and $\mathit{\alpha }{\stackrel{\mathrm{̃}}{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:

$$\begin{array}{}\text{(31)}& {\displaystyle }{\displaystyle }{u}_{\mathrm{R}}\left(t\right)={\widehat{\mathrm{\Phi }}}_{\mathrm{2}}^{-\mathrm{1}}\left[{\widehat{\mathrm{\Phi }}}_{\mathrm{2}}\left[u\left(t\right)\right]\mathit{\chi }\left(\stackrel{\mathrm{̃}}{p}\right)\right],\end{array}$$

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

4 Examples of processing COSMIC data

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.

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.