Generalized canonical transform method for radio occultation sounding with improved retrieval in the presence of horizontal gradients

By now, a series of advanced wave optical approaches to the processing of radio occultation (RO) observations are widely used. In particular, the canonical transform (CT) method and its further developments need to be mentioned. The latter include the full spectrum inversion (FSI) method, the geometric optical phase matching (PM) method, and the general approach based on the Fourier integral operators (FIOs), also referred to as the CT type 2 (CT2) method. The general idea of these methods is the application of a canonical transform that changes the coordinates in the phase space from time and Doppler frequency to impact parameter and bending angle. For the spherically symmetric atmosphere, the impact parameter, being invariant for each ray, is a unique coordinate of the ray manifold. Therefore, the derivative of the phase of the wave field in the transformed space is directly linked to the bending angle as a single-valued function of the impact parameter. However, in the presence of horizontal gradients, this approach may not work. Here we introduce a further generalization of the CT methods in order to reduce the errors due to horizontal gradients. We describe, in particular, the modified CT2 method, denoted CT2A, which complements the former with one more affine transform: a new coordinate that is a linear combination of the impact parameter and bending angle. The linear combination coefficient is a tunable parameter. We derive the explicit formulas for the CT2A and develop the updated numerical algorithm. For testing the method, we performed statistical analyses based on RO retrievals from data acquired by the Constellation Observing System for Meteorology, Ionosphere, and Climate (COSMIC) and collocated analysis profiles of the European Centre for Medium-Range Weather Forecasts (ECMWF). We demonstrate that it is possible to find a reasonably optimal value of the new tunable CT2A parameter that minimizes the root mean square difference between the RO retrieved and the ECMWF refractivity in the lower troposphere and allows the practical realization of the improved capability to cope with horizontal gradients and serve as the basis of a new quality control procedure.

1 Introduction

The first step in the development of the wave optical (WO) approach to the processing of radio occultation (RO) observations was made by Melbourne et al. (1994) who used the thin screen approximation for the atmosphere combined with the back propagation (BP) technique. This approach was further developed under the name of Fresnel inversion by Mortensen and Høeg (1998). Although the accuracy of this approximation in the lower troposphere was insufficient for practical application, its basic idea was correct. It consisted in the reduction of the influence of the diffraction by using the BP, which made the inversion results independent from the observation distance and canceled the resolution restriction due to the Fresnel zone size.

Later works (Gorbunov et al., 1996a; Karayel and Hinson, 1997; Gorbunov and Gurvich, 1998a, b) developed a different understanding of the BP technique. The BP wave field evaluated in some plane was not considered as the actual wave field but as a representation of the original field observed at low Earth orbit (LEO); in this representation, the effects of diffraction and multipath propagation were significantly reduced. This, in a straightforward way, allowed the evaluation of the geometric optical (GO) bending angle profile, which was inverted in the framework of the standard GO scheme (Ware et al., 1996; Kursinski et al., 1997).

The further development of the WO approach based on the representation view relied upon the concept of the canonical transform (CT) originating from classical mechanics (Arnold, 1978; Goldstein et al., 2014), generalized for quantum mechanics by Fock (1978) and mathematically substantiated by Egorov (1985) and Egorov and Shubin (1993). Later on, this concept obtained an extensive mathematical development (Treves, 1982a, b; Hörmander, 1985a, b). The correspondence between quantum and classical mechanics is the same as the link between wave optics and geometrical optics.

In both cases, there is a strict mathematical representation (quantum mechanics or wave optics) and its asymptotic solution (classical mechanics or geometrical optics). While the evolution of de Broglie waves of probability or electromagnetic waves is described by the Hamilton operator, the evolution of rays or classical trajectories of particles is described by the Hamilton system. The Hamilton operator is obtained from the Hamilton function by the substitution of the momentum operator instead of classical momentum. Accordingly, for the classical problem, the phase space is introduced, the dimension of which equals double the geometric dimension because to each geometrical coordinate, we can conjugate the corresponding momentum. For the wave problems, momentum is understood as the ray direction vector.

The canonical transforms arise when we consider the class of the transforms of the phase space that conserve the canonical form of the Hamilton dynamical system. It was first demonstrated by Fock (1978) that these transforms have a very simple implementation in quantum mechanics: they correspond to linear transforms of the wave function. The kernel of this transform is derived in classical terms, but still it describes a short-wave asymptotic solution of the wave problem. This idea was later mathematically developed first by Egorov (1985) and Egorov and Shubin (1993) and then by Treves (1982a, b) and Hörmander (1985a, b).

The application of the CT approach for the RO observation processing was pioneered by Gorbunov (2002), who combined it with the BP. The idea of the CT without BP was first developed by Jensen et al. (2003, 2004), and later the general view of these results in the framework of the CT approach was developed by Gorbunov and Lauritsen (2004a, b). Finally, it was recognized that the different methods – CT (Gorbunov, 2002), full spectrum inversion (FSI) (Jensen et al., 2003), phase matching (PM) (Jensen et al., 2004), and CT of the second type (CT2) (Gorbunov and Lauritsen, 2004a) – were, in fact, different approximations of the same solution, for which Fourier integral operators (FIOs) provided the general transform approach (Gorbunov and Lauritsen, 2004a).

The idea of the CT approach is as follows. Given the observations or RO complex signal u(t) as a function of time t, which can be represented through its amplitude A(t) and phase ϕ(t), u(t)=A(t)exp (iϕ(t)). It is convenient to use eikonal or phase path $\mathrm{\Psi }\left(t\right)=\mathit{\varphi }\left(t\right)/k$, where $k=\mathrm{2}\mathit{\pi }/\mathit{\lambda }$ is the wavenumber, and λ is the wavelength. Thus, u(t)=A(t)exp (ikΨ(t)), and k is the large parameter. The signal is composed of multiple sub-signals u_i(t)=A_i(t)exp (ikΨ_i(t)) corresponding to interfering rays. For each sub-signal, it is possible to introduce the instantaneous frequency $k\dot{{\mathrm{\Psi }}_i}=k{\mathit{\sigma }}_i$. However, the instantaneous frequency cannot be introduced for their composition.

The multipath propagation problem consists in the decomposition of the signal equal to the sum or different sub-signals to retrieve the ray structure of the observed field. The solution of this problem discussed in the aforementioned papers consisted in the transform of the observed wave field u(t) into a different representation. The new coordinates in the transformed space were the ray impact parameter p and bending angle ϵ. The transform $\left(t,\mathit{\sigma }\right)\to \left(p,\mathit{ϵ}\right)$ is canonical (Gorbunov and Lauritsen, 2004a), which allows for writing the corresponding linear transform ${\widehat{\mathrm{\Phi }}}_{\mathrm{2}}$, where the subscript 2 indicates that it is a CT of the second type (Arnold, 1978; Goldstein et al., 2014) that maps the original field u(t) to a field in the impact parameter representation $\widehat{u}\left(p\right)={\widehat{\mathrm{\Phi }}}_{\mathrm{2}}\left[u\left(t\right)\right]\left(p\right)$. The idea of the choice of the ray impact parameter as the new coordinate is based on the fact that in a spherically symmetric medium, the ray impact parameter is the ray invariant, which is known as Bouguer's law. The locally spherically symmetric medium is the basic approximation used in the inversion of RO data. For the real atmosphere with horizontal gradients, the dynamic equation for p was derived by Gorbunov and Kornblueh (2001), who demonstrated that the derivative of p with respect to the ray arc length is equal to the horizontal component of the refractivity gradient in the occultation plane. Strong horizontal gradients may result in the situation when dependence ϵ(p) becomes multivalued (Healy, 2001; Gorbunov and Lauritsen, 2009), which was referred to as the impact parameter multipath (Zou et al., 2019).

The idea explored in the present manuscript consists in the further development of the CT approach by using a generalized transform with the coordinate $p^{\prime }=p+\mathit{\beta }\mathit{ϵ}$. Unlike the standard CT approach, for which the form of the new coordinates in the phase is known in advance, this transform has the tunable parameter β that can take into account the statistical impact parameter multipath effect.

The paper is organized as follows. In Sect. 2, we discuss the canonical transform in wave optics and quantum mechanics in general terms, including a brief review of FIOs. Based on this context, we discuss in Sect. 3 the application of the CT method for RO and introduce the particular phase space and the specific choice of coordinates, as well as the new generalization, adding an affine transform with a tunable parameter for improving the coping capability with horizontal gradients. In Sect. 4, we discuss the practical modifications needed to readily advance existing numerical implementations of the CT algorithm and present results of our performance evaluation from processing actual observed Constellation Observing System for Meteorology, Ionosphere, and Climate (COSMIC) RO data, including how to find an optimal value of the tunable parameter to minimize the retrieval errors in the lower troposphere. Section 5 finally provides the summary and main conclusions of the paper.

2 General concept of canonical transform in wave optics

The canonical transforms (CTs) in classical mechanics are a class of transforms of the coordinates and momenta conserving the Hamiltonian form of the dynamical equation (Arnold, 1978; Goldstein et al., 2014). Fock (1978) introduced the CT into quantum mechanics. Note that the first Russian edition of the monograph by Fock appeared as early as 1929. Because the relation between classical and quantum mechanics, on one side, and the relation between geometrical and wave optics, on the other side, are the same, we can immediately apply the approach introduced by Fock (1978).

We assume that the wave field can be represented in the standard form:

$$\begin{array}{}\text{(1)}& u\left(t\right)=A\left(t\right)\exp \left(ik\mathrm{\Psi }\left(t\right)\right),\end{array}$$

where t is the observation time, Ψ(t) is the eikonal, $k=\mathrm{2}\mathit{\pi }/\mathit{\lambda }$ is the wavenumber, λ is the wavelength, and A(t) is the amplitude. The time t can be associated with a specific spatial location of the observation, as is the case in RO, but u(t) can also be looked at as a generic signal.

The amplitude A(t) and the derivative of Ψ(t) are assumed to be slowly changing within an oscillation period. In this case, the wave field is termed quasi-monochromatic with an instant amplitude A(t) and frequency $\mathit{\omega }\left(t\right)=k\dot{\mathrm{\Psi }}\left(t\right)$. Otherwise, more generally, the field should be equal to a super-position of quasi-monochromatic components:

$$\begin{array}{}\text{(2)}& u\left(t\right)=\sum _j{A}^{\left(j\right)}\left(t\right)\exp \left(ik{\mathrm{\Psi }}^{\left(j\right)}\left(t\right)\right),\end{array}$$

where the upper index j enumerates the components, A^(j)(t) are their amplitudes, and Ψ^(j)(t) are their eikonals. Each component has its own instant amplitude and frequency.

When discussing the CTs, it is necessary to bear in mind that most of the relations have an asymptotic nature, in which k is the large parameter (or λ is the small parameter). The reason is as follows. Given measurements of wave field, each monochromatic component can be interpreted in terms of wave fronts and rays defined in terms of instant tones of the signal. At the observation point at time moment t, each component has a single ray, and its direction is linked to the normalized frequency $\mathit{\sigma }\left(t\right)=\dot{\mathrm{\Psi }}\left(t\right)$ through the geometry of the observation trajectory.

Therefore, for a specific class of signals, including quasi-monochromatic ones and their superpositions, it is possible to introduce a phase space (t,σ). Although the original signal is 1-D, this space is 2-D, and the structure of the signal can be described in terms of the function σ(t) which can be both single valued for quasi-monochromatic signals or multivalued for their superpositions.

Consider RO observations. The original signal corresponds to a range of rays starting at the transmitter and the phase space σ(t) is a very smooth continuous line. As the signal propagates through the atmosphere, its structure gets more and more complicated. Still, in the phase space, its topological structure remains the same; it is always a single continuous line, although it may not be single-valued with respect to time t, which corresponds to multipath propagation (Gorbunov, 2002; Gorbunov and Lauritsen, 2004a). Such a line representing the signal structure is referred to as the ray manifold (Mishchenko et al., 1990).

The outstanding and still simple idea of Fock (1978) was that the classical CTs correspond to linear integral transforms of the wave field with oscillating kernels. This class of transform was later named Fourier integral operators (FIOs) (Egorov, 1985; Egorov and Shubin, 1993; Treves, 1982a, b; Hörmander, 1985a, b). The general form of such an operator first discussed by Fock (1978) has the following form:

$$\begin{array}{}\text{(3)}& \begin{array}{rl}\widehat{u}\left(p\right)& =\sqrt{-{\displaystyle \frac{ik}{\mathrm{2}\mathit{\pi }}}}\int a_{\mathrm{2}}\left(p,t\right)\exp \left(ik{S}_{\mathrm{2}}\left(p,t\right)\right)u\left(t\right)\mathrm{d}t\\ & \equiv {\mathrm{\Phi }}_{\mathrm{2}}\left[u\left(t\right)\right]\left(p\right),\end{array}\end{array}$$

where p is a new coordinate in the mapped space. We use notation Φ₂ and, accordingly, a₂ and S₂ because this type of operator was referred to as the FIO of the second type (Gorbunov and Lauritsen, 2004a), while the FIO of the first type is the composition of a Fourier transform and a second type FIO (Egorov, 1985; Egorov and Shubin, 1993). This type of operator is linked to the corresponding type of the CT generating function (Arnold, 1978; Goldstein et al., 2014). Note that historically FIOs of the second type appeared first, but in mathematical works, it was FIOs of the first type that were discussed first.

Considering now u(t) as a quasi-monochromatic signal, we can derive the asymptotic form of transform (Eq. 3) using the stationary phase principle:

$$\begin{array}{}\text{(4)}& \begin{array}{rl}\widehat{u}\left(p\right)& =\sqrt{-{\displaystyle \frac{ik}{\mathrm{2}\mathit{\pi }}}}\int a_{\mathrm{2}}\left(p,t\right)A\left(t\right)\exp \left(ik(S_{\mathrm{2}}\left(p,t\right)+\mathrm{\Psi }\left(t\right))\right)\mathrm{d}t\\ & \equiv {\mathrm{\Phi }}_{\mathrm{2}}\left[u\left(t\right)\right]\left(p\right).\end{array}\end{array}$$

The stationary phase point t_s(p) of this integral satisfies the equation

$$\begin{array}{}\text{(5)}& {\displaystyle \frac{\partial }{\partial t}}{S}_{\mathrm{2}}\left(p,t\right)+\dot{\mathrm{\Psi }}\left(t\right)=\mathrm{0}.\end{array}$$

Accordingly, the transformed field, under the assumption that Eq. (5) has a single solution t_s(p), is also quasi-monochromatic and can be written as follows:

$$\begin{array}{}\text{(6)}& \begin{array}{rl}\widehat{u}\left(p\right)& =A^{\prime }\left(p\right)\exp \left(ik{\mathrm{\Psi }}^{\prime }\left(p\right)\right)\\ & =A^{\prime }\left(p\right)\exp \left(ik\left(S_{\mathrm{2}}\left(p,t_{\mathrm{s}}\left(p\right)\right)+\mathrm{\Psi }\left(t_{\mathrm{s}}\left(p\right)\right)\right)\right).\end{array}\end{array}$$

Its instantaneous frequency equals the following:

$$\begin{array}{}\text{(7)}& \begin{array}{rl}\mathit{\xi }\left(p\right)& =\dot{{\mathrm{\Psi }}^{\prime }}\left(p\right)={\displaystyle \frac{d}{\mathrm{d}p}}\left(S_{\mathrm{2}}\left(p,t_{\mathrm{s}}\left(p\right)\right)+\mathrm{\Psi }\left(t_{\mathrm{s}}\left(p\right)\right)\right)=\\ & ={\displaystyle \frac{\partial }{\partial p}}{S}_{\mathrm{2}}\left(p,t_{\mathrm{s}}\left(p\right)\right)+{\left.\left({\displaystyle \frac{\partial }{\partial t}}{S}_{\mathrm{2}}\left(p,t\right)+\dot{\mathrm{\Psi }}\left(t\right)\right)\right|}_{t=t_{\mathrm{s}}\left(p\right)}{\displaystyle \frac{\mathrm{d}{t}_s}{\mathrm{d}p}}\\ & ={\displaystyle \frac{\partial }{\partial p}}{S}_{\mathrm{2}}\left(p,t_{\mathrm{s}}\left(p\right)\right),\end{array}\end{array}$$

by virtue of Eq. (5). Recalling that $\dot{\mathrm{\Psi }}\left(t\right)=\mathit{\sigma }$, which is the original momentum, we have the following relation between the canonical coordinates (t,σ) and (p,ξ) in the original and mapped spaces:

$$\begin{array}{}\text{(8)}& {\displaystyle \frac{\partial S_{\mathrm{2}}}{\partial t}}=-\mathit{\sigma },{\displaystyle \frac{\partial S_{\mathrm{2}}}{\partial p}}=\mathit{\xi },\end{array}$$

which can be expressed in terms of the differential dS₂:

$$\begin{array}{}\text{(9)}& \mathrm{d}{S}_{\mathrm{2}}=\mathit{\xi }\mathrm{d}p-\mathit{\sigma }\mathrm{d}t,\end{array}$$

and, vice versa, the requirement that the right-hand part in Eq. (9) should be equal to a full differential dS₂ of a function S₂(p,t) is a necessary and sufficient condition for the transform $(t,\mathit{\sigma })\to (p,\mathit{\sigma })$ to be canonical (Arnold, 1978; Goldstein et al., 2014). The function S₂(p,t) is then termed the generating function of the canonical transform.

In terms of FIO, S₂(p,t) is referred to as its phase function, and a₂(p,t) is its amplitude function. The phase function, which specifies the canonical transform, is of primary importance, while the amplitude function is derived using energy conservation (Gorbunov and Lauritsen, 2004a). We see, therefore, that by using classical or geometric optical concepts, it is possible to write down the asymptotic form of the quantum or wave optical operator, implementing the transformation of the original signal into a different representation. If the structure of the original signal is represented as a ray manifold in the phase plane, such a transform is applied to the coordinates in this space. In particular, it may be possible to find such a coordinate system in which the ray manifold geometry will be exceptionally simple.

Figure 1 Radio occultation observation geometry with relevant geometrical variables indicated (for description see Sect. 3.1)

3 The canonical transform method for RO and its generalization

Here we discuss the application of the CT technique for the analysis of RO observations (Fig. 1) by first reviewing the different existing variants (Sect. 3.1) and then introducing the new generalized CT method (Sect. 3.2) and an application-relevant formulation for readily updating existing algorithms (Sect. 3.3).

3.1 Canonical transform method in different existing variants

The RO observation geometry is schematically represented in Fig. 1. The wave emitted by a transmitter Tx is received by a receiver Rx in low Earth orbit. The transmitter is borne by a satellite belonging to one of the modern Global Navigation Satellites Systems (GNSS), including GPS, GLONASS, Galileo, etc. Due to the movement of the transmitter and receiver, the ray descends or ascends in the atmosphere, which allows for the derivation of the atmospheric profiles from the bending angles ϵ(p) (Ware et al., 1996; Kursinski et al., 1997). The CT technique is used for the retrieval of the bending angle profile from the wave field measurements.

The first approach of processing RO data, belonging to the class of CT, was back propagation (BP) (Gorbunov et al., 1996a; Karayel and Hinson, 1997; Gorbunov and Gurvich, 1998a, b). In this technique, the field was linearly transformed to be recalculated to the BP plane locate at coordinate x_B:

$$\begin{array}{}\text{(10)}& \begin{array}{rl}{u}_B\left(y\right)& =\sqrt{{\displaystyle \frac{ik}{\mathrm{2}\mathit{\pi }}}}\int {\displaystyle \frac{u\left(t\right)\exp \left(-ik\left|{\mathbf{r}}_B\left(y\right)-{\mathbf{r}}_R\left(t\right)\right|\right)}{{\left|{\mathbf{r}}_B\left(y\right)-{\mathbf{r}}_R\left(t\right)\right|}^{\mathrm{1}/\mathrm{2}}}}\times \\ & \times \left|\sin \mathit{\varphi }\left({\dot{\mathbf{r}}}_R\left(t\right),{\mathbf{r}}_B\left(y\right)-{\mathbf{r}}_R\left(t\right)\right){\dot{\mathbf{r}}}_R\left(t\right)\right|\mathrm{d}t,\end{array}\end{array}$$

where 2-D vector r_B(y) equals (x_B,y), and ϕ(a,b) is the angle between vectors a and b. This transform is preceded by the stationarization of the transmitting satellite and projection of the satellite movement to the vertical plane. It is important that the BP field is not the real field in the BP plane because the BP procedure assumes vacuum propagation. This procedure results in some representation of the original wave field with reduced diffraction effects due to the reduction of the propagation distance. The new coordinate y is more favorable for finding a unique projection of the ray manifold that disentangles the multipath propagation. Still, this coordinate is not the best choice.

A much better coordinate for the new representation should be the impact parameter p because in a spherically symmetric medium, it is an invariant for each ray due to Bouguer's law, and thus it is unique for each ray. A dynamic equation for the variation of p along the ray as a function of the horizontal gradient of refractivity was obtained by Gorbunov and Kornblueh (2001). The idea of complementing the BP technique with one more transform that maps the field to the impact parameter representation was pioneered by Gorbunov (2002). It was the first application of the FIO of the first type, which has the form

$$\begin{array}{}\text{(11)}& \begin{array}{rl}\widehat{u}\left(p\right)& =\sqrt{-{\displaystyle \frac{ik}{\mathrm{2}\mathit{\pi }}}}\int a_{\mathrm{1}}\left(p,\mathit{\sigma }\right)\exp \left(ik{S}_{\mathrm{1}}\left(p,\mathit{\sigma }\right)\right)\stackrel{\mathrm{̃}}{u}\left(\mathit{\sigma }\right)\mathrm{d}\mathit{\sigma }\\ & \equiv {\mathrm{\Phi }}_{\mathrm{1}}\left[u\left(t\right)\right]\left(p\right),\end{array}\end{array}$$

where the only difference from the second type operator is that it acts upon the Fourier-transformed field $\stackrel{\mathrm{̃}}{u}\left(\mathit{\sigma }\right)$. It can be looked at as the composition of the Fourier transform, which itself is a second type FIO and another second type FIO. Because the Fourier transform is a simple rotation of the phase space by π∕2: $(t,\mathit{\sigma })\to \left(\mathit{\sigma },-t\right)$, the equation for the phase function of the first type takes the following form (Arnold, 1978; Goldstein et al., 2014):

$$\begin{array}{}\text{(12)}& \mathrm{d}{S}_{\mathrm{1}}=\mathit{\xi }\mathrm{d}p+t\mathrm{d}\mathit{\sigma }.\end{array}$$

Gorbunov (2002) applied this operator to the back-propagated field. To this end, using the normal vector $\mathit{\nu }=\left(\mathit{\eta },\sqrt{\mathrm{1}-{\mathit{\eta }}^{\mathrm{2}}}\right)$ to the straight ray, we express the impact parameter:

$$\begin{array}{}\text{(13)}& p\left(y,\mathit{\eta }\right)=-x\mathit{\eta }+y\sqrt{\mathrm{1}-{\mathit{\eta }}^{\mathrm{2}}}.\end{array}$$

Now it is necessary to find the canonical transform $(y,\mathit{\eta })\to \left(p,\mathit{\xi }\right)$ whose characteristic property in the 2-D case is the conservation of the volume element, as follows from Eq. (12):

$$\begin{array}{}\text{(14)}& {\displaystyle \frac{\partial \mathit{\xi }}{\partial \mathit{\eta }}}{\displaystyle \frac{\partial p}{\partial y}}-{\displaystyle \frac{\partial \mathit{\xi }}{\partial y}}{\displaystyle \frac{\partial p}{\partial \mathit{\eta }}}=\mathrm{1}.\end{array}$$

It is enough to consider solutions ξ=ξ(η). Then, from Eq. (14), we readily derive

$$\begin{array}{}\text{(15)}& \begin{array}{rl}{\displaystyle \frac{\partial \mathit{\xi }}{\partial \mathit{\eta }}}& ={\left({\displaystyle \frac{\partial p}{\partial y}}\right)}^{-\mathrm{1}}={\displaystyle \frac{\mathrm{1}}{\sqrt{\mathrm{1}-{\mathit{\eta }}^{\mathrm{2}}}}},\\ & \mathit{\xi }=\arcsin \mathit{\eta }.\end{array}\end{array}$$

This results in the solution for the phase and amplitude functions:

$$\begin{array}{}\text{(16)}& \begin{array}{rl}& S_{\mathrm{2}}\left(p,\mathit{\eta }\right)=p\arcsin \mathit{\eta }-x\sqrt{\mathrm{1}-{\mathit{\eta }}^{\mathrm{2}}},\\ & a_{\mathrm{2}}\left(p,\mathit{\eta }\right)=\sqrt{{\displaystyle \frac{{\partial }^{\mathrm{2}}S}{\partial p\partial \mathit{\eta }}}}={\left(\mathrm{1}-{\mathit{\eta }}^{\mathrm{2}}\right)}^{-\mathrm{1}/\mathrm{4}}.\end{array}\end{array}$$

This defines the FIO, which is applied to the backpropagated wave field u_B(y) and produces the mapped field

$$\begin{array}{}\text{(17)}& \widehat{u}\left(p\right)=A^{\prime }\left(p\right)\exp \left(ik\int \mathit{\xi }\left(p\right)\mathrm{d}p\right).\end{array}$$

The derivative ξ(p) of its eikonal is algebraically linked to the bending angle:

$$\begin{array}{}\text{(18)}& \mathit{ϵ}\left(p\right)=-\mathit{\xi }\left(p\right)-\arcsin \left({\displaystyle \frac{x_{T}p+y_T\sqrt{r_T^{\mathrm{2}}-p^{\mathrm{2}}}}{r_T^{\mathrm{2}}}}\right),\end{array}$$

where $\left(x_T,y_T\right)={\mathbf{r}}_T$ is the transmitter position in the occultation plane. Because the cross term in S₂, which depends both on p and η, is linear with respect to p, the integration over new coordinate ξ=arcsin η turns it to pξ, and, therefore, the operator is reduced to the Fourier transform in combination with a nonlinear change in coordinates. This indicates that this operator allows for a fast implementation. A similar idea will be applied below.

The complicated nature of the BP+CT algorithm stimulated further studies (Gorbunov and Lauritsen, 2002, 2004b) where the idea was expressed of applying the FIO directly to the observed wave field u(t), without intermediate and numerically expensive steps like BP. Full spectrum inversion (FSI), developed by Jensen et al. (2003), was the first solution of this type, although with some restrictive assumptions. However, the general solution was just 1 year away; the phase matching (PM) method was developed by Jensen et al. (2004) and then put into the context of the CT approach by Gorbunov and Lauritsen (2004a), who introduced an approach based on the linearized canonical transform that reduced the FIO to the composition of nonlinear coordinate changes and the Fourier transform. This algorithm was termed the second type CT or CT2. An important advantage of the PM and CT2 methods consists in the fact that they operate with the real transmitter and receiver orbits without the stationarization.

In order to arrive at the phase function of the FIO of the second type, consider the expression for the derivative of the phase of the observed wave field:

$$\begin{array}{}\text{(19)}& \dot{\mathrm{\Psi }}=\mathit{\sigma }(p,y)=p\dot{\mathit{\theta }}+{\displaystyle \frac{{\dot{r}}_T}{r_T}}\sqrt{r_T^{\mathrm{2}}-p^{\mathrm{2}}}+{\displaystyle \frac{{\dot{r}}_R}{r_R}}\sqrt{r_R^{\mathrm{2}}-p^{\mathrm{2}}}.\end{array}$$

Using Eq. (8), we derive the phase function:

$$\begin{array}{}\text{(20)}& \begin{array}{rl}{S}_{\mathrm{2}}& \left(p,t\right)=-\int \left(p\dot{\mathit{\theta }}+{\displaystyle \frac{{\dot{r}}_T}{r_T}}\sqrt{r_T^{\mathrm{2}}-p^{\mathrm{2}}}+{\displaystyle \frac{{\dot{r}}_R}{r_R}}\sqrt{r_R^{\mathrm{2}}-p^{\mathrm{2}}}\right)\mathrm{d}y=\\ & =-\int \left(pd\mathit{\theta }+{\displaystyle \frac{\mathrm{d}{r}_T}{r_T}}\sqrt{r_T^{\mathrm{2}}-p^{\mathrm{2}}}+{\displaystyle \frac{\mathrm{d}{r}_R}{r_R}}\sqrt{r_R^{\mathrm{2}}-p^{\mathrm{2}}}\right)=\\ & =-p\mathit{\theta }-\sqrt{r_T^{\mathrm{2}}-p^{\mathrm{2}}}+p\arccos {\displaystyle \frac{p}{r_T}}-\sqrt{r_R^{\mathrm{2}}-p^{\mathrm{2}}}\\ & +p\arccos {\displaystyle \frac{p}{r_R}},\end{array}\end{array}$$

where θ, r_T, and r_R are functions of time t. We do not reproduce here the derivation of the amplitude function a₂(p,t), which uses simple geometrical considerations (Gorbunov and Lauritsen, 2004a). This phase function, although providing the accurate solution, has a disadvantage: its cross term depending on both p and t, generally speaking, cannot be decomposed as g₁(p)g₂(t), and the FIO cannot be reduced to a Fourier transform in compositions with nonlinear coordinate changes. This is only possible in some particular cases, e.g., for circular orbits, when the phase function equals pθ, and using θ as a new coordinate instead of time reduces the operator to the Fourier transform. This method was referred to as FSI.

To find an approximate solution that significantly reduces the computational costs at the expense of an insignificant reduction of accuracy, the representation of the approximate impact parameter was introduced. The impact parameter p is a function of t,σ: $p=p(t,\mathit{\sigma })$. We introduce its approximation $\stackrel{\mathrm{̃}}{p}$:

$$\begin{array}{}\text{(21)}& \begin{array}{rl}\stackrel{\mathrm{̃}}{p}& \left(t,\mathit{\sigma }\right)=p_{\mathrm{0}}\left(t\right)+{\displaystyle \frac{\partial p_{\mathrm{0}}}{\partial \mathit{\sigma }}}\left(\mathit{\sigma }-{\mathit{\sigma }}_{\mathrm{0}}\left(t\right)\right)=f\left(t\right)+{\displaystyle \frac{\partial p_{\mathrm{0}}}{\partial \mathit{\sigma }}}\mathit{\sigma },\\ & f\left(t\right)=p_{\mathrm{0}}\left(t\right)-{\displaystyle \frac{\partial p_{\mathrm{0}}}{\partial \mathit{\sigma }}}{\mathit{\sigma }}_{\mathrm{0}}\left(t\right)=\\ & =p_{\mathrm{0}}-{\left({\displaystyle \frac{\mathrm{d}\mathit{\theta }}{\mathrm{d}t}}-{\displaystyle \frac{\mathrm{d}{r}_G}{\mathrm{d}t}}{\displaystyle \frac{p_{\mathrm{0}}}{r_G\sqrt{r_G^{\mathrm{2}}-p_{\mathrm{0}}^{\mathrm{2}}}}}-{\displaystyle \frac{\mathrm{d}{r}_L}{\mathrm{d}t}}{\displaystyle \frac{p_{\mathrm{0}}}{r_L\sqrt{r_L^{\mathrm{2}}-p_{\mathrm{0}}^{\mathrm{2}}}}}\right)}^{-\mathrm{1}}{\mathit{\sigma }}_{\mathrm{0}},\end{array}\end{array}$$

where σ₀(t) is a smooth model of normalized Doppler frequency, $p_{\mathrm{0}}\left(t\right)=p(t,{\mathit{\sigma }}_{\mathrm{0}}(t\left)\right)$, and $\partial p_{\mathrm{0}}/\partial \mathit{\sigma }=\partial p/\partial \mathit{\sigma }{|}_{\mathit{\sigma }={\mathit{\sigma }}_{\mathrm{0}}\left(t\right)}$. We now parameterize the trajectory with the coordinate Υ=Υ(t). For brevity, we use the notation u(Υ) instead of u(t(Υ)). For the coordinate Υ and the corresponding momentum η, we use the following definitions:

$$\begin{array}{}\text{(22)}& \begin{array}{rl}\mathrm{d}\mathrm{{\rm Y}}& ={\left({\displaystyle \frac{\partial p_{\mathrm{0}}}{\partial \mathit{\sigma }}}\right)}^{-\mathrm{1}}\mathrm{d}t={\displaystyle \frac{\partial \mathit{\sigma }}{\partial p_{\mathrm{0}}}}\mathrm{d}t,\\ & \mathit{\eta }={\displaystyle \frac{\partial p_{\mathrm{0}}}{\partial \mathit{\sigma }}}\mathit{\sigma }.\end{array}\end{array}$$

Finally, we arrive at the following linear canonical transform $(\mathrm{{\rm Y}},\mathit{\eta })\to (p,\mathit{\xi })$:

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

The generating function of this canonical transform is easily computed from the differential equation:

$$\begin{array}{}\text{(24)}& \begin{array}{rl}\mathrm{d}{S}_{\mathrm{2}}& =\mathit{\xi }\mathrm{d}\stackrel{\mathrm{̃}}{p}-\mathit{\eta }\mathrm{d}\mathrm{{\rm Y}}=-\mathrm{{\rm Y}}\mathrm{d}\stackrel{\mathrm{̃}}{p}-(\stackrel{\mathrm{̃}}{p}-f(\mathrm{{\rm Y}}\left)\right)\mathrm{d}\mathrm{{\rm Y}}\\ & S_{\mathrm{2}}(\stackrel{\mathrm{̃}}{p},\mathrm{{\rm Y}})=-\stackrel{\mathrm{̃}}{p}Y+\underset{\mathrm{0}}{\overset{\mathrm{{\rm Y}}}{\int }}f\left({\mathrm{{\rm Y}}}^{\prime }\right)\mathrm{d}{\mathrm{{\rm Y}}}^{\prime }.\end{array}\end{array}$$

For the new coordinate Υ, we have the following relation:

$$\begin{array}{}\text{(25)}& \mathrm{d}\mathrm{{\rm Y}}=\mathrm{d}\mathit{\theta }-{\displaystyle \frac{\mathrm{d}{r}_G}{r_G}}{\displaystyle \frac{p_{\mathrm{0}}}{\sqrt{r_G^{\mathrm{2}}-p_{\mathrm{0}}^{\mathrm{2}}}}}-{\displaystyle \frac{\mathrm{d}{r}_L}{r_L}}{\displaystyle \frac{p_{\mathrm{0}}}{\sqrt{r_L^{\mathrm{2}}-p_{\mathrm{0}}^{\mathrm{2}}}}}.\end{array}$$

For circular orbits, this approximation, once again, reduces to FSI. To evaluate the bending angle, we use the fact that the momentum of the field in the mapped space equals −Υ. We also evaluate the accurate impact parameter p as follows. Given the dependence $\mathrm{{\rm Y}}\left(\stackrel{\mathrm{̃}}{p}\right)$, it is possible to find the corresponding time $t\left(\stackrel{\mathrm{̃}}{p}\right)$. Using Eq. (21), we infer

$$\begin{array}{ll}{\displaystyle }\mathit{\sigma }\left(\stackrel{\mathrm{̃}}{p}\right)& {\displaystyle }=\left(\stackrel{\mathrm{̃}}{p}-p_{\mathrm{0}}\left(t\left(\stackrel{\mathrm{̃}}{p}\right)\right)\right){\left({\displaystyle \frac{\partial p_{\mathrm{0}}}{\partial \mathit{\sigma }}}\right)}^{-\mathrm{1}}+{\mathit{\sigma }}_{\mathrm{0}}\left(t\left(\stackrel{\mathrm{̃}}{p}\right)\right),\\ \text{(26)}& {\displaystyle }p\left(\stackrel{\mathrm{̃}}{p}\right)& {\displaystyle }=p(t\left(\stackrel{\mathrm{̃}}{p}\right),\mathit{\sigma }\left(\stackrel{\mathrm{̃}}{p}\right)).\end{array}$$

Finally, for each impact parameter p, we determine the coordinate $\mathrm{{\rm Y}}\left(p\right)=-\mathit{\xi }\left(p\right)$ and, therefore, the corresponding moment of time t=t(Υ(p)) when this ray was observed; the bending angle is then evaluated from the geometrical relation:

$$\begin{array}{}\text{(27)}& \begin{array}{rl}\mathit{ϵ}\left(p\right)& =\mathit{\theta }\left(t\left(\mathrm{{\rm Y}}\left(p\right)\right)\right)-\arccos {\displaystyle \frac{p}{r_T\left(t\left(\mathrm{{\rm Y}}\left(p\right)\right)\right)}}\\ & -\arccos {\displaystyle \frac{p}{r_R\left(t\left(\mathrm{{\rm Y}}\left(p\right)\right)\right)}}.\end{array}\end{array}$$

This method, termed CT2, indicates both a high accuracy and numerical performance. This discourse leads us to the conclusion that there is a family of closely related WO methods that are based on the same principle. The observed wave field is subjected to a linear integral operator with an oscillating kernel that transforms the field into a different representation. The representation is chosen in such a way that the projection of the ray manifold to the new coordinate axis is unique. The operation is also referred to as the unfolding of multipath. Finally, such methods as CT, FSI, PM, and CT2 involve the evaluation of the same integral transform under different assumptions and approximations. The difference in the results of the application of these WO methods is less significant than the difference coming from other parts of RO data processing systems, including cut-off, filtering, and quality control (QC) procedures (Gorbunov et al., 2004, 2011).

3.2 Generalized canonical transform method

All the modifications of the CT approach discussed above relied upon impact parameter p as the unique coordinate of the ray manifold. However, the impact parameter is, generally speaking, not invariant for each ray, and its perturbations due to horizontal gradients may result in breaking the above condition. To see this, consider the ray equations in the Hamilton form. They are derived from the Hamilton function

$$\begin{array}{}\text{(28)}& H\left(\mathbf{r},\mathbf{p}\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left({\mathbf{p}}^{\mathrm{2}}-n^{\mathrm{2}}\left(\mathbf{r}\right)\right),\end{array}$$

where p is the momentum, and n(r) is the refractivity field. The Hamilton system has the following form:

$$\begin{array}{l}{\displaystyle }\dot{\mathbf{r}}={\displaystyle \frac{\partial H}{\partial \mathbf{p}}},\dot{\mathbf{p}}=-{\displaystyle \frac{\partial H}{\partial \mathbf{r}}},\dot{\mathrm{\Psi }}=\mathbf{p}\dot{\mathbf{r}},\\ \text{(29)}& {\displaystyle }\dot{\mathbf{r}}=\mathbf{p},\dot{\mathbf{p}}=n\mathrm{\nabla }n,\dot{\mathrm{\Psi }}=n^{\mathrm{2}},\end{array}$$

where p is the classical momentum. Because $\left|\mathbf{p}\right|=\left|\mathrm{\nabla }\mathrm{\Psi }\right|=n$, we arrive at the following differential relation between the parameter τ of this system, the ray arc length s, and the eikonal

$$\begin{array}{}\text{(30)}& \mathrm{d}\mathit{\tau }={\displaystyle \frac{\mathrm{d}s}{n}},\mathrm{d}\mathrm{\Psi }=n\mathrm{d}s.\end{array}$$

Equation (29) has a form that is specific for the Cartesian coordinates. Consider an arbitrary coordinate system with the metric tensor g_ij: ds²=dxⁱg_ijdx^j, where xⁱ are the components of vector r, and we follow the Einstein tensor notation implying the summation over each pair of upper and lower indexes of the same name. If we define the momentum by the relation $p_i=g_{ij}{\dot{x}}^j$, the form pdr is invariant, the transform to the new coordinates (pⁱ, xⁱ) is canonical, and the canonical form of the Hamilton system also remains invariant (Arnold, 1978), provided that the Hamilton function is defined as follows:

$$\begin{array}{}\text{(31)}& H\left(\mathbf{r},\mathbf{p}\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left(p_i{g}^{ij}{p}_j-n^{\mathrm{2}}\left(\mathbf{r}\right)\right),\end{array}$$

where g^ij is the matrix inverse to g_ij. This results in the following form of the ray equations:

$${\dot{x}}^i={\displaystyle \frac{\partial H}{\partial p_i}}=g^{ij}{p}_j,{\dot{p}}_i=-{\displaystyle \frac{\partial H}{\partial x^i}}=n{\displaystyle \frac{\partial n}{\partial x^i}}-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{p}_k{\displaystyle \frac{\partial g^{kj}}{\partial x^i}}{p}_j.$$

The 2-D approximation (Zou et al., 2002) allows us to treat rays as plane curves. Consider polar coordinates (r,θ) with the metric tensor:

$$\begin{array}{}\text{(32)}& g_{ij}=\left(\begin{array}{cc}\mathrm{1}& \mathrm{0}\\ \mathrm{0}& r^{\mathrm{2}}\end{array}\right),g^{ij}=\left(\begin{array}{cc}\mathrm{1}& \mathrm{0}\\ \mathrm{0}& r^{-\mathrm{2}}\end{array}\right).\end{array}$$

Then we have the following equations:

$$\begin{array}{ll}{\displaystyle }{p}_{\mathit{\theta }}& {\displaystyle }=r^{\mathrm{2}}\dot{\mathit{\theta }}=nr{\displaystyle \frac{r\mathrm{d}\mathit{\theta }}{\mathrm{d}s}}=nr\sin \mathit{\psi },\\ {\displaystyle }{\dot{p}}_{\mathit{\theta }}& {\displaystyle }=n{\displaystyle \frac{\partial n}{\partial \mathit{\theta }}},\\ \text{(33)}& {\displaystyle }{\dot{p}}_r& {\displaystyle }=\ddot{r}=n{\displaystyle \frac{\partial n}{\partial r}}+{\displaystyle \frac{p^{\mathrm{2}}}{r^{\mathrm{3}}}},\end{array}$$

where ψ is the angle between vectors $\dot{\mathbf{r}}$ and r. The angular component of the momentum p_θ coincides with the ray impact parameter p, which is invariant in a spherically layered medium but is perturbed by the horizontal gradients (Gorbunov et al., 1996b; Gorbunov and Kornblueh, 2001; Healy, 2001; Gorbunov and Lauritsen, 2009).

Figure 2 Impact parameter multipath, old coordinate (impact parameter) lines, and modified coordinate lines.

The variations in the ray impact parameter, which is no longer an invariant coordinate in the ray space, seem to undermine the elegant idea of the CT approach. Still, the CT method can be applied using the same formulas, but the coordinate p will now acquire a different meaning: it will be understood as the “effective impact parameter”, i.e., the impact parameter which would result in the observed Doppler frequency shift if the atmosphere were spherically layered (Gorbunov et al., 2019; Gorbunov and Lauritsen, 2009). Accordingly, the evaluated bending angle will also be the “effective” bending angle. The reason is that for the evaluation of the real bending angle, understood as the angle between the ray directions at the transmitter and receiver, two corresponding values of the impact parameter are required, which cannot be derived from the single variable, the Doppler frequency. This, by itself, is not a significant problem because the assimilation of bending angle profiles can be based on the effective values (Gorbunov et al., 2019) provided that the observation operator correctly implements their evaluation.

More importantly, horizontal gradients may result in multivalued ray manifold projections when using the effective impact parameter p as the coordinate in the mapped space. This situation is termed “impact parameter multipath” (Zou et al., 2019). Theoretically, for any ray manifold perturbation, there always exists an unfolding coordinate transform. This follows from the fact that topologically the ray manifold is always a continuous line without self-crossing. However, this coordinate transform depends on the a priori unknown horizontal gradients of refractivity.

A typical multivalued bending angle profile (Gorbunov and Lauritsen, 2009; Zou et al., 2019) is shown in Fig. 2. From numerical simulations, it can be inferred that there is a kind of asymmetry; the impact parameter multipath manifests itself mostly in ascending spikes but hardly in descending spikes. Accordingly, in order to better unfold the multipath, it must be possible to use another coordinate in such a way that the modified coordinate lines are sloped. Therefore, we modify the transform Eq. (23) in order to use another coordinate,

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

where β is a tunable parameter and has dimensions of length/angle (km/rad). Although the optimal value of this parameter should be different for individual events, the aforementioned asymmetry results in the conclusion that the preferred value of β is expected to be negative. Therefore, it may be possible to find its optimal value that, in the statistical sense, will minimize errors due to the impact parameter multipath.

The modified canonical transform Eq. (23) is written as follows:

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

Using the modified function f^′(Υ) instead of the original one, we will 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 very simple modification of the existing CT2 algorithm. Using the numerical implementation of the modified CT will allow us to study its influence upon the RO inversion statistics in the lower troposphere.

We denote the generalized FIO ${\widehat{\mathrm{\Phi }}}_{\mathrm{2}}^{\left(\mathit{\beta }\right)}u\left(\stackrel{\mathrm{̃}}{p}\right)$. Now we can consider the wave field in the impact parameter representation, $\widehat{u}\left(\mathit{\beta };{\stackrel{\mathrm{̃}}{p}}^{\prime }\right)={\widehat{\mathrm{\Phi }}}_{\mathrm{2}}^{\left(\mathit{\beta }\right)}u\left({\stackrel{\mathrm{̃}}{p}}^{\prime }\right)$. The standard CT algorithm corresponds to the evaluation of $\widehat{u}\left(\mathrm{0};\stackrel{\mathrm{̃}}{p}\right)={\widehat{\mathrm{\Phi }}}_{\mathrm{2}}^{\left(\mathrm{0}\right)}u\left(\stackrel{\mathrm{̃}}{p}\right)$ with β=0.

It is possible to arrive at a quantitative estimate of β based on Gorbunov and Kornblueh (2001), Gorbunov and Lauritsen (2009), and Zou et al. (2019). We expect that $\left|\mathit{\beta }\right|\mathit{\lesssim }\mathit{\delta }p/\mathit{\delta }\mathit{ϵ}$, where δp is the typical variation in impact parameter due to the horizontal gradients, and δϵ is the corresponding bending angle variation. Assuming that δp≈0.1 km, and δϵ≈0.01 rad, we arrive at a first quantitative estimate of $\mathit{\beta }\approx -\mathrm{10}$ km/rad.

3.3 Affine transform for updating existing CT algorithms

Modification of existing numerical algorithms may not be so straightforward as it follows from the above mathematical considerations. In order to avoid this, it is possible to complement an existing implementation of any WO-based numerical algorithm by an additional affine transform.

We will now derive the transform between $\widehat{u}\left(\mathrm{0};\stackrel{\mathrm{̃}}{p}\right)$ and $\widehat{u}\left(\mathit{\beta };{\stackrel{\mathrm{̃}}{p}}^{\prime }\right)$. We can write the following transform between these representations:

$$\begin{array}{ll}{\displaystyle }{\stackrel{\mathrm{̃}}{p}}^{\prime }& {\displaystyle }=\stackrel{\mathrm{̃}}{p}-\mathit{\beta }\left(\mathit{\xi }-{\mathit{\xi }}_{\mathrm{0}}\right),\\ \text{(36)}& {\displaystyle }{\mathit{\xi }}^{\prime }& {\displaystyle }=\mathit{\xi },\end{array}$$

where ξ₀ is the reference point. This is an affine transform in the $(\stackrel{\mathrm{̃}}{p},\mathit{\xi })$ plane. This suggests the abbreviation CT2A for the new generalized form, which stands for the CT2 complemented by the affine transform.

The generating function of transform (36) $S\left({\stackrel{\mathrm{̃}}{p}}^{\prime },\mathit{\xi }\right)$ is defined by

$$\begin{array}{}\text{(37)}& \mathrm{d}{S}^{\left(\mathit{\beta }\right)}=\mathit{\xi }\mathrm{d}{\stackrel{\mathrm{̃}}{p}}^{\prime }+\stackrel{\mathrm{̃}}{p}\mathrm{d}\mathit{\xi },\end{array}$$

which is equivalent to the following system:

$$\begin{array}{ll}{\displaystyle \frac{\partial S^{\left(\mathit{\beta }\right)}}{\partial {\stackrel{\mathrm{̃}}{p}}^{\prime }}}& {\displaystyle }=\mathit{\xi },\\ \text{(38)}& {\displaystyle \frac{\partial S^{\left(\mathit{\beta }\right)}}{\partial \mathit{\xi }}}& {\displaystyle }=\stackrel{\mathrm{̃}}{p}={\stackrel{\mathrm{̃}}{p}}^{\prime }+\mathit{\beta }\left(\mathit{\xi }-{\mathit{\xi }}_{\mathrm{0}}\right).\end{array}$$

From this, we can conclude that

$$\begin{array}{}\text{(39)}& S^{\left(\mathit{\beta }\right)}\left({\stackrel{\mathrm{̃}}{p}}^{\prime },\mathit{\xi }\right)={\stackrel{\mathrm{̃}}{p}}^{\prime }\mathit{\xi }+\mathit{\beta }{\displaystyle \frac{{\left(\mathit{\xi }-{\mathit{\xi }}_{\mathrm{0}}\right)}^{\mathrm{2}}}{\mathrm{2}}}.\end{array}$$

This phase function defines the FIO of the first type:

$$\begin{array}{}\text{(40)}& \begin{array}{rl}\widehat{u}\left(\mathit{\beta };{\stackrel{\mathrm{̃}}{p}}^{\prime }\right)& =\sqrt{-{\displaystyle \frac{ik}{\mathrm{2}\mathit{\pi }}}}\int \exp \left(ik{S}^{\left(\mathit{\beta }\right)}\left({\stackrel{\mathrm{̃}}{p}}^{\prime },\mathit{\xi }\right)\right)\stackrel{\mathrm{̃}}{u}\left(\mathit{\xi }\right)\mathrm{d}\mathit{\xi }\\ & \equiv {\widehat{{\mathrm{\Phi }}^{\prime }}}_{\mathrm{1}}^{\left(\mathit{\beta }\right)}\left[u\left(t\right)\right]\left(\stackrel{\mathrm{̃}}{p}\right).\end{array}\end{array}$$

Finally, we can write the operator relation

$$\begin{array}{}\text{(41)}& {\widehat{\mathrm{\Phi }}}_{\mathrm{2}}^{\left(\mathit{\beta }\right)}={\widehat{{\mathrm{\Phi }}^{\prime }}}_{\mathrm{1}}^{\left(\mathit{\beta }\right)}{\widehat{\mathrm{\Phi }}}_{\mathrm{2}}^{\left(\mathrm{0}\right)},\end{array}$$

which can be used for the modification of the existing version of operator ${\widehat{\mathrm{\Phi }}}_{\mathrm{2}}^{\left(\mathrm{0}\right)}$.

The above derivation allows for one more generalization. We can consider β=β(ξ). In this case, the phase function is derived in a straightforward way:

$$\begin{array}{}\text{(42)}& S^{\left(\mathit{\beta }\right)}\left({\stackrel{\mathrm{̃}}{p}}^{\prime },\mathit{\xi }\right)={\stackrel{\mathrm{̃}}{p}}^{\prime }\mathit{\xi }+\int \mathit{\beta }\left(\mathit{\xi }\right)\left(\mathit{\xi }-{\mathit{\xi }}_{\mathrm{0}}\right)\mathrm{d}\mathit{\xi }.\end{array}$$

Using $\mathit{\beta }\left(\mathit{\xi }\right)=\sum {\mathit{\beta }}_j{\mathit{\xi }}^j$ results in a simple analytical expression for S^(β) with a set of tuning parameters β_j. In this work, we, however, use a constant β.

Figure 3 Statistics for latitude band 0–10^∘. (a) Amount of data. (b) RMS relative difference in refractivity COSMIC–ECMWF ΔN_CE. (c) RMS relative difference in refractivity CT2A–CT2. All are functions of the parameter β.

Figure 4 Statistics for latitude band 10–20^∘. (a) Amount of data. (b) RMS relative difference in refractivity COSMIC–ECMWF ΔN_CE. (c) RMS relative difference in refractivity CT2A–CT2. All are functions of the parameter β.

Figure 5 Statistics for latitude band 20–30^∘. (a) Amount of data. (b) RMS relative difference in refractivity COSMIC–ECMWF ΔN_CE. (c) RMS relative difference in refractivity CT2A–CT2. All are functions of the parameter β.

4 Implementation and numerical performance evaluation

Our implementation of the CT2A algorithm was based on the existing program code with the addition of the parameter β and using the modified function f^′(Υ) as defined by Eq. (35). Practically, this only required the modification of a few lines in the program code that implements the CT2 method, as well as the implementation of one more command line parameter.

In our numerical validation, we retrieved COSMIC refractivity profiles N_C using COSMIC data from the year 2008 on the 1st and 15th of every month, leading to a total of 24 d and altogether around 60 000 RO events. We used collocated ECMWF refractivity profiles N_E, i.e., interpolated to the corresponding COSMIC RO event location, as the reference. To this end, we employed ECMWF analyses at 1^∘ latitudinal and longitudinal resolution with 91 vertical levels covering the altitude range up to about 80 km. The refractivity was evaluated from pressure, temperature, and humidity fields. The tangent point drift was taken into account. We used the root mean square (RMS) relative difference in COSMIC from ECMWF ΔN_CE, defined as $\sqrt{〈{\left(\left(N_{\mathrm{C}}-N_{\mathrm{E}}\right)/N_{\mathrm{E}}\right)}^{\mathrm{2}}〉}$, which includes both systematic and random deviations.

Figure 6 Statistics for latitude band 30–40^∘. (a) Amount of data. (b) RMS relative difference in refractivity COSMIC–ECMWF ΔN_CE. (c) RMS relative difference in refractivity CT2A–CT2. All are functions of the parameter β.

Figure 7 Statistics for latitude band 40–50^∘. (a) Amount of data. (b) RMS relative difference in refractivity COSMIC–ECMWF ΔN_CE. (c) RMS relative difference in refractivity CT2A–CT2. All are functions of the parameter β.

Figure 8 Statistics for latitude band 50–60^∘. (a) Amount of data. (b) RMS relative difference in refractivity COSMIC–ECMWF ΔN_CE. (c) RMS relative difference in refractivity CT2A–CT2. All are functions of the parameter β.

Figures 3 through 11 show the statistical values of ΔN_CE as a function of latitude and parameter β. We averaged over 10^∘ wide latitude bands including both the Southern Hemisphere and Northern Hemisphere. The parameter β changed in the interval from −4 to −12 km/rad with a step of −1.

These results indicate that for the latitudes of 0–30^∘ and for the altitudes of 2.5 km, the application of the CT2A algorithm, together with our quality control (QC) procedure, results in the reduction of the RMS relative difference in refractivity profiles COSMIC–ECMWF ΔN_CE. The amount of data in the altitude range below 3 km that passes the QC slightly decreases with increasing β. Above the height of about 0.4 km, the increasing value of β reduces ΔN_CE. Below 0.4 km, there is an optimal value of β in the interval from −4 to −9 km/rad, depending on the altitude and latitude. For the latitudes higher than 30^∘, the application of CT2A is not expedient.

Figure 9 Statistics for latitude band 60–70^∘. (a) Amount of data. (b) RMS relative difference in refractivity COSMIC–ECMWF ΔN_CE. (c) RMS relative difference in refractivity CT2A–CT2. All are functions of the parameter β.

Figure 10 Statistics for latitude band 70–80^∘. (a) Amount of data. (b) RMS relative difference in refractivity COSMIC–ECMWF ΔN_CE. (c) RMS relative difference in refractivity CT2A–CT2. All are functions of the parameter β.

Figure 11 Statistics for latitude band 80–90^∘. (a) Amount of data. (b) RMS relative difference in refractivity COSMIC–ECMWF ΔN_CE. (c) RMS relative difference in refractivity CT2A–CT2. All are functions of the parameter β.

The statistics for different β, presented in the above figures, was evaluated independently, i.e., the statistical ensembles were different. Figures 12 and 13 show the statistics for the common datasets including events passing QC for both β=0 and current values of β. The statistical differences between refractivity retrieved with β=0 and other values of β are indistinguishably small (never exceeding a level of 0.0005 %). Here the reduction of ΔN_CE is slightly more than in Figs. 3 and 4

Figure 12 Statistics for latitude band 0–10^∘ evaluated for subsets common for β=0 and each other value of β. (a) Amount of data. (b) RMS relative difference in refractivity COSMIC–ECMWF ΔN_CE. (c) RMS relative difference in refractivity CT2A–CT2. All are functions of the parameter β.

Figure 13 Statistics for latitude band 10–20^∘ evaluated for subsets common for β=0 and each other value of β. (a) Amount of data. (b) RMS relative difference in refractivity COSMIC–ECMWF ΔN_CE. (c) RMS relative difference in refractivity CT2A–CT2. All are functions of the parameter β.

This indicates that CT2A acts as a QC procedure not involving any external data and only based on the internal properties of observed signals. On average, CT2A provides a higher cut-off height which is estimated from the CT amplitude by correlating it with the θ function (Gorbunov et al., 2006). By looking at the ray manifold in the phase space from different directions, it is possible to choose ray manifold pieces where its structure is most stable.

5 Summary and conclusions

In this study, we discussed the general idea of the canonical transform (CT) method and provided a new generalization adding more flexibility for application in RO processing. CTs in classical mechanics (geometrical optics) are implemented in quantum mechanics (wave optics) by linear operators with oscillating kernels. Such operators are referred to as Fourier integral operators (FIOs). During the past century, this approach acquired a solid theoretical basis. In numerous mathematical monographs, one finds the advanced theory of FIOs. The central role in this theory is played by the concept of the ray manifold and its projections.

In quantum mechanics and wave optics, FIOs were employed for the quantization procedure, i.e., the construction of the asymptotic quantum (quasi- or semiclassical) solutions on the basis of the classical (geometric optical) ones. The idea of the CT method for processing RO observations is the reverse: the reconstruction of the geometric optical solution from the wave optical one, which can be referred to as the dequantization.

Although there have been many modifications, like original CT combined with back propagation (BP), full spectrum inversion (FSI), phase matching (PM), and CT of type 2 (CT2), there is no essential difference between these FIO-based methods. The difference consists of the approximation of the phase function of the FIO leading to the corresponding approximate representation of the impact parameter and bending angle and in the specific implementation (such as cut-off, filtering, and quality control procedures). All these methods map the wave field into the representation of the impact parameter p. This choice of the coordinate in the mapped space has its reasons: in the case of a spherically symmetric medium, the impact parameter is always a unique coordinate of the ray manifold.

The implementation of this idea in the real, non-spherically symmetric atmosphere encounters some difficulties. First, in the strict sense, there is no such quantity as the impact parameter as a unique variable any more, but it is still possible to operate with the effective impact parameter derived from Doppler frequency shift using the same relations as for a spherically symmetric medium. This quantity can be implemented in the observation operator for the variational assimilation of RO observations, canceling errors due to horizontal gradients. However, the above property of the impact parameter, which is supposed to be a unique coordinate of the ray manifold, does not always hold for the effective value. In some cases, the situation referred to as the impact parameter multipath may occur, resulting in retrieval errors in atmospheric profiles derived from RO data.

In order to mitigate this fundamental shortcoming, we introduced a generalization of the CT approach. We used a generalized definition of the coordinate in phase space, defined as a linear combination of the impact parameter and bending angle. Because this can be understood as an affine transform of the phase space, we coined the abbreviation CT2A for the new method. This transform has a parameter β which can be tuned to optimize the algorithm's performance.

We implemented the CT2A algorithm by modifying our existing program code for the CT2 method. In order to evaluate its statistical performance under real RO observation conditions, including challenging horizontal gradients in the lower troposphere, we processed a large ensemble of COSMIC RO data for the year 2008 on the 1st and 15th day of every month, adding up to a total of about 60 000 RO events. We used the total relative difference in COSMIC from collocated ECMWF analysis profiles over the lower troposphere as the metric for this evaluation and the tuning parameter estimation.

For the latitudes of 0–30^∘ and for the altitudes between 0.4 and 2.5 km, the application of the CT2A algorithm decreased the COSMIC–ECMWF difference metric with increasing parameter β. For the altitudes below 0.4 km, the optimal value of parameter β is found to be −4 to −9 km/rad. This was achieved on account of a slight decrease in the amount of data passing through the whole retrieval chain including the QC. This indicates that the CT2A itself implements a QC procedure that does not involve any external information about the atmospheric refractivity but is only based on the analysis of the structure of the observed signals.

Overall, these results suggest that the CT2A method is not only theoretically an innovative generalization of the CT or FIO class of methods but also practically a valuable advancement for RO processing in that it can improve the capability to cope with challenging horizontal gradient conditions in the lower troposphere and serve as the basis of a new QC procedure.