Using A Sliding Window Phase Matching Method for Imaging of GNSS Radio Occultation Signals

: Global Navigation Satellite System Radio Occultation (GNSS-RO) is a technique used to sound the atmosphere and derive vertical proﬁles of refractivity. Signals from GNSS satellites are received in a low-Earth orbit, and they are then processed to produce bending angle proﬁles, from which meteorological parameters can be retrieved. Generating two-dimensional images in the form of spectrograms from GNSS-RO signals is commonly done to, for instance, investigate reﬂections or estimate signal quality in the lower troposphere. This is typically implemented using, e.g., the Short-Time Fourier Transform (STFT) to produce a time-frequency representation that is subsequently transformed to bending angle (BA) and impact height (IH) coordinates by non-linear mapping. In this paper, we propose an alternative method based on a straightforward extension of the Phase Matching (PM) operator to produce two-dimensional spectral images in the BA-IH domain by applying a sliding window. This Sliding Window Phase Matching (SWPM) method generates the spectral amplitude on an arbitrary grid in BA and IH, e.g., along the coordinate axes. To illustrate, we show both SWPM and STFT methods applied to operational MetOp-A data. For SWPM we use a constant window in the BA-dimension, whereas for STFT we use a conventional constant time window. We show that the SWPM method produces the same result as STFT when the same window length is used for both methods. The sample points in impact parameter and bending angle are those generated by and the main advantage is that SWPM offers the user a convenient way to freely sample the BA-IH space. The cost for this is processing time that is somewhat longer than implementations based on the Fast Fourier Transform, such as the STFT method.


Introduction
Since their inception, Global Navigation Satellite System (GNSS) constellations have come to play an important part not only for geographic positioning, but also for Earth observation purposes [1]. GNSS Radio Occultation (GNSS-RO) is an atmospheric sounding technique that evaluates the refraction of GNSS signals as they propagate through the Earth's atmosphere [2]. Measurements are retrieved in the form of one-dimensional (vertical) bending angle (BA) profiles as a function of impact height (IH). These can later be inverted to profiles describing the refractive index of the atmosphere by applying the inverse Abel transform [3].
Time-frequency representations of GNSS-RO measurements have been widely used for a long time, as they visualize BA profiles as two-dimensional spectral images rather than one-dimensional profiles (see [4][5][6][7][8]). They have since found a wider use. For instance, Beyerle et al. [9] applied the multiple signal classification algorithm (MUSIC) to identify reflected components in CHAMP [10] data. They described how surface reflections result in distinct deviations in Doppler frequency, visualized in spectrograms generated from RO signals. This idea has been developed further, leading to automatic classification of reflections [11][12][13] using supervised machine learning techniques. These classification methods rely on formulating a computer vision problem, where machine learning algorithms are trained to recognize the features of reflections in spectrogram images. Aparicio et al. [14] used a time-frequency representation of GNSS-RO signals to filter out direct rays from RO signals. After inverting the spectrogram back to RO signals, they could retrieve the bending angles of reflected rays. The use of one-sided windows to remove the deeper direct rays while keeping reflected rays in the signal was proposed by Sievert et al. [15]. Since reflected rays are received at a higher point in orbit than some direct rays, they were able to show that reflected rays could potentially be separated from the direct rays by truncating the RO signal. Interference from other GNSS transmitters has been visualized in spectrogram images [16,17]. One method of error estimation and quality assessment based on the spectral analysis of BA profiles was suggested by Gorbunov et al. [18], and it has recently been revisited by Liu et al. [19]. The core idea is to use the local spectral width as a metric of quality, where a narrow local spectrum implies high quality.
Although it is very common to make time-frequency representations using the Short-Time Fourier Transform (STFT) or related methods, Gorbunov et al. [20] instead applied the Wigner Distribution Function (WDF) [21]. Although WDF offers a sharper resolution than STFT, it also introduces cross terms that could interfere with the main signal. These cross terms are minimized by applying a sliding window, which in turn degrades the resolution. WDF is computationally complex, which was recently addressed by applying the Kirkwood Distribution Function (KDF) [22], which can be reduced to a single Fourier transform. However, KDF differs from aforementioned methods in that rays are not straightforwardly observed in the amplitude, but instead by finding the stationary phase of the output function.
Given a time-frequency representation, non-linear mapping is needed to transform the image to the BA and IH domain [23]. In this paper, we propose a method of extending the Phase Matching (PM) [24] operator to visualize RO profiles as two-dimensional images. The proposed extension applies a sliding window to PM, and it does not need a separate mapping step. This method allows us to define the length of the sliding window in terms of BA instead of time. The rest of this paper is organized as follows: In Section 2, we describe the application of STFT to GNSS-RO measurements. In Section 3, we describe PM and the modification that enables the application of a sliding window, and we show its relation to STFT. In Section 4, we show the application of this method on MetOp data, and in Section 5, we provide further discussion and a conclusion.

Short-Time Fourier Transform
As mentioned in the introduction, STFT has been used in the context of GNSS-RO for a long time. In this section, we provide a description of the procedure along with derivations that explain how it maps the time signal to BA and IH. This procedure is well known in the literature, see e.g., [20,23].
In processing an RO event, ephemeris data describing the geometry of transmitter and receiver relative to the Earth's local curvature is required in conjunction with the complex phase and amplitude of the received signal. The geometry of an RO event is shown in Figure 1, along with the nomenclature used for the ephemeris data in this paper. To improve the readability in figures, we define the term IH as a − r c , where r c is the Earth's local radius of curvature.
We define a GNSS-RO signal u(t) as where v(t) is real, and φ(t) is the signal phase in radians. Under the assumption that the signal is a ray that passed through a spherically symmetric atmosphere, we can say that where k is the wave number corresponding to the carrier frequency, S(t) is the optical path length (OPL) of a ray, and η(t) is a general noise term. The OPL (described in e.g., [24]) is given by where r L (t) is the LEO radius, r G (t) is the GNSS radius, a(t) is the impact parameter, and the bending angle α(t) is derived from the angles of the quadrilateral in Figure 1, namely: where θ(t) is referred to as the separation angle between LEO and GNSS. The time derivative of the OPL iṡ Considering the signal in a small time window, [t 0 − T, t 0 + T], we can assume that the OPL can be approximated by the first two Taylor terms, where S(t 0 ) is shorthand for the OPL evaluated at time t 0 , i.e., S(t)| t 0 . This approximation is valid when making the approximation very good for window lengths of a few seconds. In this time window, the signal is thus approximated as It is clear that the term kṠ(t 0 ) represents a frequency that is constant over the time window. If the signal at a given time comprises multiple rays, so-called multipath, the exponent would involve several OPL terms.
When applying STFT, we look at small time windows, [t 0 − T, t 0 + T], and we perform the Fourier transform on the down-converted signal to identify rays that correspond to different values of impact parameter. The down-conversion is performed by subtracting some range model, R M (t), namely The range model must have the same dynamics as the OPL of a ray, such that in a small time window it can also be approximated using the first two Taylor terms, i.e., This model is commonly selected from a climatological model [21,25], or a smooth BA profile retrieved with geometrical optics. In the time window, we have The Fourier transform is performed, generating a complex function The majority contribution to the amplitude at a specific value for the frequency ω comes from the signal phase components that match the applied down-conversion frequency and the frequency in the Fourier transform. Recalling the Taylor approximation in Equations (6) and (8), the specific frequency ω can then be connected to the presence of a ray through Substituting Equation (13) into Equation (5) allows us to solve for the impact parameter a(t 0 , ω) corresponding to time t 0 and frequency ω, througḣ For this impact parameter, we can then use Equation (4) to find the corresponding BA. Thus, for every time window, we calculate pairs of BAs and impact parameters corresponding to the frequencies used in the Fourier transform. To implement this numerically entails a set of discrete frequencies, forming a band of BAs and impact parameters around the range model, where the band width is given by the sampling frequency of the signal. The smeared features in Figure 2 demonstrate how these bands are arranged along the range model. Figure 3 shows the spectrogram that results from applying STFT to a GNSS-RO measurement in time-frequency coordinates (left) and BA-IH coordinates (right) respectively.
When using FFT for performing STFT, the resolution is limited by the number of sample points in the frequency domain. The sample spacing in time is where N is the number of sample points. The maximum frequency represented is ω max = 2π/2∆t, and the frequency range [−ω max , ω max ] is divided into N samples, meaning that the frequency resolution is This means that the height of the pixels in the mapped spectrogram is constrained by the length of the signal segments used in FFT. This issue can be solved by implementing the discrete Fourier transform on a customized vector of ω-values. The BA, on the other hand, may be sampled as close as desired by using overlapping time segments. One further, practical drawback of the STFT technique, regardless of whether FFT is used or not, is that the procedure results in three matrices: one for amplitude, one for BA, and one for IH, where the ordering of the values for BA and IH is not trivial. This is not a problem for plotting the amplitude, but to further process the information in the image it is likely that interpolation to a regular grid in BA and IH is needed. We show in Section 3 that the SWPM method does not have these problems.

Sliding Window Phase Matching
Phase matching [24] belongs to a class of operators referred to as Fourier integral operators, along with Full Spectrum Inversion (FSI) [26], Canonical Transform and the second Canonical Transform (CT and CT2, respectively) [27][28][29]. Along with a higher resolution than geometrical optics retrieval, this class of operators addresses the shortcoming of the standard geometrical optics approach [2]; that is, they resolve multipath phenomena, which geometrical optics does not. Although FSI and CT2 rely on a Fourier transform and a subsequent mapping from the frequency domain to the impact parameter domain, PM instead maps the RO signal directly to impact parameter. Practically, FSI and CT2 make use of approximations of impact parameter to enable the use of the Fourier transform, and thus the efficient FFT. On the other hand, PM cannot be reduced to a Fourier transform, resulting in higher computational complexity.
In extending PM to SWPM, we consider a small window in BA, [α 0 − ∆α, α 0 + ∆α], and perform an integral over a down-converted signal. The down-conversion is performed by subtracting a specific range model R N (t, a x ), which depends both on time and on a fixed value for the impact parameter a x , namely where The BA window can be mapped to a time window using Equation (18), meaning that The time derivative of the phase model iṡ It has the same property as the OPL in Equation (5), in that over a small window of time it can be approximated by the first two Taylor terms, The down-converted signal here becomes and in the time window we have The SWPM integral is performed: Similar to STFT, the amplitude value of the integral is determined mainly by the component of the signal that has a frequency of the OPL that corresponds to the frequency of the range model applied. This occurs when orṙ which in turn means that a(t 0 ) = a x . The amplitude of this function is thus directly identified to belong to a ray with BA α 0 and impact parameter a x . Essentially, the most intuitive way to use the SWPM technique is to calculate the spectral amplitude values in a coordinate system where IH is on the y-axis, and BA is on the x-axis.

Connection with STFT
Suppose we apply the range model used in SWPM (R N (t, a x )) as we perform STFT. We then havê When ω = 0, this is exactly the same integral as Equation (24), i.e., which means that a(t 0 , ω = 0) = a x .
When ω = 0, the equality will be different. The STFT integral for ω = 0 will correspond to the SWPM integral for a slightly different value of impact parameter, i.e., u(t 0 , T, ω, a x ) =ũ(t 0 , T, a x + δa x ), (30) where This equation can be solved for δa x (ω). Therefore, to compare the two methods one simply needs to perform SWPM for all the values of δa x (ω) that correspond to the frequency range used in the STFT integral. An even easier way to perform this comparison is to use STFT for a certain time window and take note of the array of impact parameter values generated. This array is then used as values for a x in the SWPM integral. The comparison can be performed without even calculating the BA values corresponding to the impact parameter array and time t 0 . To illustrate, we present an example of such a comparison in Figures 4-6. Figure 4 shows a mapped spectrogram with various cross sections highlighted. Figures 5 and 6 presents comparisons of the horizontal and vertical cross sections, respectively. The amplitudes from STFT and SWPM are overlaid, and they highlight the above derivation. The minute differences amount to numerical noise.

Resolution in Impact Parameter
We have already concluded that we find the impact parameter in a signal segment by matching the frequency in the segment, using a range model with a fixed value for the impact parameter. Now, this matching is not exact, and there is an interval around this impact parameter that contributes significantly to the amplitude. We need to estimate the width of this interval. Consider a signal with unity amplitude in a short time segment [t 0 − T, t 0 + T]. It is given by We assume that the value of the impact parameter that characterizes this segment is a x . If we apply the range model with a slightly different value for the impact parameter a x = a x + δa x , the SWPM integral (Equation (24)) becomes where δω = k[Ṡ(t 0 ) −Ṙ N (t 0 , a x )], and it is understood thatṠ(t 0 ) −Ṙ N (t 0 , a x ) = 0. Since the phase derivative is given by Equations (5) and (20), we have The last two terms may be neglected, and we have δω ≈ kθδa x . The integral evaluates to |ũ(t 0 , T, a x )| = 2T sin(δωT) This function is well known, and it becomes zero when δωT = π. This means that we cannot distinguish a frequency in the signal more precisely than ∆ω U = 2δω = 2π/T, giving us the uncertainty in frequency, ∆ω U ≥ 2π/T. Translating this into impact parameter we get

Resolution in Bending Angle
Using the SWPM method, we specify a BA value and integrate over the signal in a region around this value, using the mapping given by Equation (4). We have seen that the integral will give us an amplitude value provided there is any signal segment inside this time domain that matches the frequency of the range model. The uncertainty of the BA is thus directly given by the BA window used for the integration, [α 0 − ∆α, α 0 + ∆α]. This means that the uncertainty in BA is ∆α U = 2∆α. This BA window is connected to the time window through Equation (4), so The last two terms can again be safely neglected, which means that the uncertainty in BA is equal to an uncertainty in time through

Consequences for SWPM
We have seen that using FFT binds the height of the pixels in the spectrogram to the integration time interval. To sample the IH domain more freely, one should instead use STFT with the discrete Fourier transform or the SWPM method. On the other hand, the resolution in IH and BA are not independent; the uncertainty in IH is decreased by increasing the integration time interval, while the opposite is true for the uncertainty in BA. Combining the expressions for uncertainty in BA and IH, we get where λ is the wavelength. Incidentally, this is basically the same relation as Equation (11) in [30] in the case of a slowly varying atmosphere. Throughout the rest of this paper, we use a BA window of 2 mrad, which gives us a resolution in IH of around 200 m. In Figure 7, we see a comparison with the window lengths 0.5 mrad and 10 mrad, which correspond to IH uncertainties of approximately 800 m and 40 m, respectively. Clearly, the above derivations fit well with what we see in the figure.

Application of SWPM on Real RO Data
In this section, we show the results of applying SWPM to MetOp-A measurements. For reference, we also present the same measurements processed by STFT. The window applied for SWPM is a Hanning window with a length of 2 mrad. Although it is not within the scope of the paper to find an optimal window length for either method, window lengths for STFT images are chosen manually to achieve a resolution similar to those of the SWPM images. The amplitudes are normalized and presented in a dB scale, and a threshold is applied at −50 dB, i.e., all values below the threshold are set to −50 dB. In addition, in the figures, we supply the BA profiles of co-located refractivity profiles from the European Centre for Medium-range Weather Forecasts (ECMWF) as well as the ones retrieved when we applied the one-dimensional PM operator. The selection of profiles is arbitrary to highlight certain common features in measurements. The data are collected from the COSMIC Data Analysis and Archive Center. Figures 8 and 9 show clear cases of reflection around IHs of 2 km. Although PM does resolve reflected rays to an extent, the smoothing procedure applied in the impact parameter domain causes the direct and reflected parts of the BA profile to bleed into each other. Since neither SWPM nor STFT relies on any differentiation, no smoothing is performed in the BA-IH domain, and the reflections can be resolved more accurately. The PM retrievals' good agreement with the ECMWF profiles implies an event where the conditions of spherical symmetry are not violated, i.e., there is little to no impact multipath. In Figure 10 we see impact multipath, potentially implying a large-scale horizontal gradient, breaking the condition of symmetry. This is demonstrated by the fact that the BA amplitude features a spike that bends upwards in IH, which cannot be represented in a one-dimensional profile. The ECMWF profile instead features a spike in BA around 3 km that tends to infinity. Figures 11 and 12 show smeared out BA amplitudes over large spans of IH, implying significant small-scale horizontal gradients.     The figures illustrate that SWPM produces images equivalent to those from STFT. In Section 3, we showed that the methods can produce equivalent results. Despite this, the SWPM and STFT images in this section appear to differ slightly. The reason for this is twofold: (1) the SWPM images are sampled on a uniform BA-IH grid rather than on the exact same coordinates as STFT and (2) the SWPM images are produced with a window of fixed length in BA, resulting in time windows that vary slightly in length. Although arbitrary sampling and variable window length is possible using STFT, it requires abandoning the efficient FFT and mapping BA and IH to time and frequency first. By contrast, SWPM enables this simply by extending the PM operator.

Conclusions
In this paper, we present a straightforward extension to the Phase Matching (PM) operator, producing spectral images of bending angle (BA) and impact height (IH) from Global Navigation Satellite System Radio Occultation (GNSS-RO) signals. The extension involves applying a sliding window defined in terms of BA, which is then mapped to the time axis using the optical path length function for a ray. We show that the Sliding Window Phase Matching (SWPM) operator is equivalent to the Short-Time Fourier Transform (STFT), by using the kernel of PM as range model. We show the results when applying SWPM on MetOp-A data, and supply images where we applied STFT for reference. Unlike methods such as STFT, where sub-apertures are defined in the time domain, SWPM defines subapertures that slide over the BA axis instead. The resulting operator lets the user sample the BA-IH spectrum freely with a very simple addition to the PM algorithm. The most likely scenario would be to use the method in such a way as to get the spectral amplitude values onto a grid in the BA-IH domain that is both aligned with the coordinate axes and uniformly spaced. By contrast, STFT either confines the user to the slanted grid provided by the output frequencies from the Fast Fourier Transform (FFT) algorithm, forcing a compromise between resolution in BA or IH, or an implementation of the discrete Fourier transform combined with mapping every BA-IH pair to the time-frequency domain, as well as mapping the BA window. This makes SWPM a convenient tool for exploring specific regions of an RO signal, e.g., reflections or regions with impact multipath.
Since the computational complexity of FFT is lower than quadratic, execution time is not an issue for STFT. SWPM, on the other hand, does have quadratic computational complexity. However, the method is trivial to run on parallel cores, and optimized implementation could bring SWPM to near real-time execution if desired. Furthermore, the problem is bounded not only by the length of an RO signal, but also by the fact that timefrequency analysis such as STFT or SWPM is useful only in the lower parts of the atmosphere.