Open AccessTechnical Note

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

Thomas Sievert

^1,*

Joel Rasch

²,

Anders Carlström

Vinícius Ludwig Barbosa

Mats I. Pettersson

and

Viet Vu

Department of Mathematics and Natural Sciences, Blekinge Institute of Technology, 371 79 Karlskrona, Sweden

Molflow, 411 01 Gothenburg, Sweden

RUAG Space AB, 405 15 Gothenburg, Sweden

Author to whom correspondence should be addressed.

Remote Sens. 2021, 13(5), 970; https://doi.org/10.3390/rs13050970

Submission received: 18 February 2021 / Revised: 26 February 2021 / Accepted: 1 March 2021 / Published: 4 March 2021

(This article belongs to the Section Atmospheric Remote Sensing)

Download

Browse Figures

Versions Notes

Abstract

Global Navigation Satellite System Radio Occultation (GNSS-RO) is a technique used to sound the atmosphere and derive vertical profiles of refractivity. Signals from GNSS satellites are received in a low-Earth orbit, and they are then processed to produce bending angle profiles, 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 reflections 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.

Keywords:

radio occultation; GNSS; radio-holographic analysis; phase matching

Graphical Abstract

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

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

u (t) = v (t) exp (i ϕ (t)),

(1)

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

ϕ (t) = k S (t) + η (t),

(2)

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

S (t) = \sqrt{r_{L} {(t)}^{2} - a {(t)}^{2}} + \sqrt{r_{G} {(t)}^{2} - a {(t)}^{2}} + α (t) a (t) + \int_{a (t)}^{\infty} α (x) d x,

(3)

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:

α (t) = θ (t) + \arcsin (\frac{a (t)}{r_{L} (t)}) + \arcsin (\frac{a (t)}{r_{G} (t)}) - π,

(4)

where

θ (t)

is referred to as the separation angle between LEO and GNSS. The time derivative of the OPL is

\dot{S} (t) = \frac{{\dot{r}}_{L} (t)}{r_{L} (t)} \sqrt{r_{L} {(t)}^{2} - a {(t)}^{2}} + \frac{{\dot{r}}_{G} (t)}{r_{G} (t)} \sqrt{r_{G} {(t)}^{2} - a {(t)}^{2}} + a (t) \dot{θ} (t) .

(5)

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,

S (t, t_{0}) \approx S (t_{0}) + (t - t_{0}) \dot{S} (t_{0}),

(6)

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

T ≪ \frac{| \dot{S} (t_{0}) |}{| \ddot{S} (t_{0}) |} \sim 1000,

(7)

making the approximation very good for window lengths of a few seconds. In this time window, the signal is thus approximated as

u (t; t_{0}, T) \approx v (t) exp (i k S (t_{0}) + i k \dot{S} (t_{0}) (t - t_{0}) + i η (t)) .

(8)

It is clear that the term

k \dot{S} (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

u_{d o w n} (t) = u (t) exp (- i k R_{M} (t)) .

(9)

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

R_{M} (t; t_{0}, T) \approx R_{M} (t_{0}) + (t - t_{0}) {\dot{R}}_{M} (t_{0}) .

(10)

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

u_{d o w n} (t; t_{0}, T) \approx u (t) exp (- i k R_{M} (t_{0}) - i k {\dot{R}}_{M} (t_{0}) (t - t_{0})) .

(11)

The Fourier transform is performed, generating a complex function

\begin{matrix} \hat{u} (t_{0}, T, ω) = \int_{t_{0} - T}^{t_{0} + T} u_{d o w n} (t; t_{0}, T) exp (- i ω t) d t \\ \approx exp (- i k R_{M} (t_{0})) \int_{t_{0} - T}^{t_{0} + T} u (t) exp (- i k {\dot{R}}_{M} (t_{0}) (t - t_{0})) exp (- i ω t) d t . \end{matrix}

(12)

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

k \dot{S} (t_{0}) - k {\dot{R}}_{M} (t_{0}) - ω = 0 .

(13)

Substituting Equation (13) into Equation (5) allows us to solve for the impact parameter

a (t_{0}, ω)

corresponding to time

t_{0}

and frequency

ω

, through

\frac{{\dot{r}}_{L} (t)}{r_{L} (t)} \sqrt{r_{L} {(t)}^{2} - a {(t_{0}, ω)}^{2}} + \frac{{\dot{r}}_{G} (t)}{r_{G} (t)} \sqrt{r_{G} {(t)}^{2} - a {(t_{0}, ω)}^{2}} + a (t_{0}, ω) \dot{θ} (t_{0}) = {\dot{R}}_{M} (t_{0}) + \frac{ω}{k} .

(14)

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

Δ t = \frac{2 T}{N},

(15)

where N is the number of sample points. The maximum frequency represented is

ω_{m a x} = 2 π / 2 Δ t

, and the frequency range

[- ω_{m a x}, ω_{m a x}]

is divided into N samples, meaning that the frequency resolution is

Δ ω = \frac{π}{T} .,

(16)

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.

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

R_{N} (t, a_{x}) = \sqrt{r_{L} {(t)}^{2} - a_{x}^{2}} + \sqrt{r_{G} {(t)}^{2} - a_{x}^{2}} + α (t, a_{x}) a_{x},

(17)

where

α (t, a_{x}) = θ (t) + \arcsin (\frac{a_{x}}{r_{L} (t)}) + \arcsin (\frac{a_{x}}{r_{G} (t)}) - π .

(18)

The BA window can be mapped to a time window using Equation (18), meaning that

[α_{0} - Δ α, α_{0} + Δ α] \leftrightarrow [t_{0} - T, t_{0} + T .]

(19)

The time derivative of the phase model is

{\dot{R}}_{N} (t, a_{x}) = \frac{{\dot{r}}_{L} (t)}{r_{L} (t)} \sqrt{r_{L} {(t)}^{2} - a_{x}^{2}} + \frac{{\dot{r}}_{G} (t)}{r_{G} (t)} \sqrt{r_{G} {(t)}^{2} - a_{x}^{2}} + a_{x} \dot{θ} (t) .

(20)

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,

R_{N} (t; t_{0}, T, a_{x}) \approx R_{N} (t_{0}, a_{x}) + (t - t_{0}) {\dot{R}}_{N} (t_{0}, a_{x}) .

(21)

The down-converted signal here becomes

u_{d o w n} (t) = u (t) exp (- i k R_{N} (t, a_{x})),

(22)

and in the time window we have

u_{d o w n} (t; t_{0}, T, a_{x}) \approx u (t) exp (- i k R_{N} (t_{0}, a_{x}) - i k {\dot{R}}_{N} (t_{0}, a_{x}) (t - t_{0})) .

(23)

The SWPM integral is performed:

\tilde{u} (t_{0}, T, a_{x}) = \int_{t_{0} - T}^{t_{0} + T} u_{d o w n} (t; t_{0}, T, a_{x}) d t \approx exp (- i k R_{M} (t_{0}, a_{x})) \int_{t_{0} - T}^{t_{0} + T} u (t) exp (- i k {\dot{R}}_{N} (t_{0}, a_{x}) (t - t_{0})) .

(24)

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

k \dot{S} (t_{0}) - k {\dot{R}}_{N} (t_{0}, a_{x}) = 0,

(25)

\begin{matrix} \frac{{\dot{r}}_{L} (t_{0})}{r_{L} (t_{0})} \sqrt{r_{L} {(t_{0})}^{2} - a {(t_{0})}^{2}} + \frac{{\dot{r}}_{G} (t_{0})}{r_{G} (t_{0})} \sqrt{r_{G} {(t_{0})}^{2} - a {(t_{0})}^{2}} + a (t_{0}) \dot{θ} (t_{0}) \\ = \frac{{\dot{r}}_{L} (t_{0})}{r_{L} (t_{0})} \sqrt{r_{L} {(t_{0})}^{2} - a_{x}^{2}} + \frac{{\dot{r}}_{G} (t_{0})}{r_{G} (t_{0})} \sqrt{r_{G} {(t_{0})}^{2} - a_{x}^{2}} + a_{x} \dot{θ} (t_{0}), \end{matrix}

(26)

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.

3.1. Connection with STFT

Suppose we apply the range model used in SWPM (

R_{N} (t, a_{x})

) as we perform STFT. We then have

\hat{u} (t_{0}, T, ω, a_{x}) \approx exp (- i k R_{N} (t_{0}, a_{x})) \int_{t_{0} - T}^{t_{0} + T} u (t) exp (- i k {\dot{R}}_{N} (t_{0}, a_{x}) (t - t_{0})) exp (- i ω t) d t .

(27)

When

ω = 0

, this is exactly the same integral as Equation (24), i.e.,

\hat{u} (t_{0}, T, ω = 0, a_{x}) = \tilde{u} (t_{0}, T, a_{x}),

(28)

which means that

a (t_{0}, ω = 0) = a_{x} .

(29)

When

ω \neq 0

, the equality will be different. The STFT integral for

ω \neq 0

will correspond to the SWPM integral for a slightly different value of impact parameter, i.e.,

\hat{u} (t_{0}, T, ω, a_{x}) = \tilde{u} (t_{0}, T, a_{x} + δ a_{x}),

(30)

where

k {\dot{R}}_{N} (t, a_{x} + δ a_{x}) - k {\dot{R}}_{N} (t_{0}, a_{x}) = ω .

(31)

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 Figure 4, Figure 5 and Figure 6. Figure 4 shows a mapped spectrogram with various cross sections highlighted. Figure 5 and Figure 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.

3.2. Resolution

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

u (t) \approx exp (i k S (t_{0})) exp (i k \dot{S} (t_{0}) (t - t_{0})) .

(32)

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

\tilde{u} (t_{0}, T, a_{x}^{'}) = exp (i k [S (t_{0}) - R_{N} (t_{0}, a_{x}^{'})]) \int_{t_{0} - T}^{t_{0} + T} u (t) exp (i δ ω (t - t_{0})) d t,

(33)

where

δ ω = k [\dot{S} (t_{0}) - {\dot{R}}_{N} (t_{0}, a_{x}^{'})]

, and it is understood that

\dot{S} (t_{0}) - {\dot{R}}_{N} (t_{0}, a_{x}) = 0

. Since the phase derivative is given by Equations (5) and (20), we have

δ ω = \frac{\partial {\dot{R}}_{N}}{\partial a_{x}} δ a_{x} = k (\dot{θ} - \frac{a_{x} {\dot{r}}_{L}}{r_{L} \sqrt{r_{L}^{2} - a_{x}^{2}}} - \frac{a_{x} {\dot{r}}_{G}}{r_{G} \sqrt{r_{G}^{2} - a_{x}^{2}}}) δ a_{x} .

(34)

The last two terms may be neglected, and we have

δ ω \approx k \dot{θ} δ a_{x}

. The integral evaluates to

| \tilde{u} (t_{0}, T, a_{x}^{'}) | = 2 T |\frac{sin (δ ω T)}{δ ω T}| .

(35)

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} \geq 2 π / T

. Translating this into impact parameter we get

Δ a_{u} \geq \frac{2 π}{k \dot{θ} T} .

(36)

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

Δ α = \frac{\partial α}{\partial t} Δ t = T (\dot{θ} - \frac{a_{x} {\dot{r}}_{L}}{r_{L} \sqrt{r_{L}^{2} - a_{x}^{2}}} - \frac{a_{x} {\dot{r}}_{G}}{r_{G} \sqrt{r_{G}^{2} - a_{x}^{2}}}) .

(37)

The last two terms can again be safely neglected, which means that the uncertainty in BA is equal to an uncertainty in time through

Δ α_{U} \approx 2 T \dot{θ} .

(38)

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

Δ α_{U} Δ a_{U} \geq \frac{4 π}{k} = 2 λ,

(39)

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.

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

Figure 8 and Figure 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. Figure 11 and Figure 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.

5. 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 sub-apertures 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 time-frequency analysis such as STFT or SWPM is useful only in the lower parts of the atmosphere.

Author Contributions

Conceptualization, T.S. and J.R.; methodology, T.S. and J.R.; software, T.S. and J.R.; validation, T.S. and V.L.B.; formal analysis, T.S.; investigation, T.S.; resources, M.I.P.; data curation, T.S.; writing–original draft preparation, T.S. and J.R.; writing–review and editing, T.S., J.R., A.C., V.L.B., M.I.P. and V.V.; visualization, T.S. and V.L.B.; supervision, M.I.P., J.R. and V.V.; project administration, M.I.P. and A.C. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Data Availability Statement

The data used in this paper are provided by the COSMIC Data Analysis and Archive Center (CDAAC), available at https://cdaac-www.cosmic.ucar.edu/cdaac/index.html (accessed on 26 February 2021).

Conflicts of Interest

The authors declare no conflict of interest.

References

Hein, G.W. Status, perspectives and trends of satellite navigation. Sat. Nav. 2020, 1, 22. [Google Scholar] [CrossRef]
Kursinski, E.; Hajj, G.; Schofield, J.; Linfield, R.; Hardy, K.R. Observing Earth’s atmosphere with radio occultation measurements using the Global Positioning System. J. Geophys. Res. Atmos. 1997, 102, 23429–23465. [Google Scholar] [CrossRef]
Fjeldbo, G.; Kliore, A.J.; Eshleman, V.R. The neutral atmosphere of Venus as studied with the Mariner V radio occultation experiments. Astron. J. 1971, 76, 123. [Google Scholar] [CrossRef]
Pavelyev, A. On the possibility of radio holographic investigation on communication link satellite-to-satellite. Radiotekhn. Elektron. 1998, 43, 939–944. [Google Scholar]
Igarashi, K.; Pavelyev, A.; Hocke, K.; Pavelyev, D.; Kucherjavenkov, I.; Matyugov, S.; Zakharov, A.; Yakovlev, O. Radio holographic principle for observing natural processes in the atmosphere and retrieving meteorological parameters from radio occultation data. Earth Planets Space 2000, 52, 893–899. [Google Scholar] [CrossRef] [Green Version]
Sokolovskiy, S.V. Modeling and inverting radio occultation signals in the moist troposphere. Radio Sci. 2001, 36, 441–458. [Google Scholar] [CrossRef] [Green Version]
Gorbunov, M. Radio-holographic analysis of Microlab-1 radio occultation data in the lower troposphere. J. Geophys. Res. Atmos. 2002, 107. [Google Scholar] [CrossRef]
Gorbunov, M.E. Radioholographic analysis of radio occultation data in multipath zones. Radio Sci. 2002, 37. [Google Scholar] [CrossRef] [Green Version]
Beyerle, G.; Hocke, K.; Wickert, J.; Schmidt, T.; Marquardt, C.; Reigber, C. GPS radio occultations with CHAMP: A radio holographic analysis of GPS signal propagation in the troposphere and surface reflections. J. Geophys. Res. Atmos. 2002, 107. [Google Scholar] [CrossRef] [Green Version]
Wickert, J.; Reigber, C.; Beyerle, G.; König, R.; Marquardt, C.; Schmidt, T.; Grunwaldt, L.; Galas, R.; Meehan, T.K.; Melbourne, W.G.; et al. Atmosphere sounding by GPS radio occultation: First results from CHAMP. Geophys. Res. Lett. 2001, 28, 3263–3266. [Google Scholar] [CrossRef] [Green Version]
Cardellach, E.; Oliveras, S.; Rius, A. GNSS signal interference classified by means of a supervised learning method applied in the time-frequency domain. In Proceedings of the 2009 2nd International Congress on Image and Signal Processing, Tianjin, China, 17–19 October 2009; pp. 1–5. [Google Scholar]
Cardellach, E.; Oliveras, S. Assessment of a Potential Reflection Flag Product. ROM SAF Report 23, Radio Occultation Meteorology Satellite Application Facility. 2016. Available online: https://www.romsaf.org/general-documents/rsr/rsr_23.pdf (accessed on 24 February 2021).
Hu, A.; Wu, S.; Wang, X.; Wang, Y.; Norman, R.; He, C.; Cai, H.; Zhang, K. Improvement of reflection detection success rate of GNSS RO measurements using artificial neural network. IEEE Trans. Geosci. Remote Sens. 2017, 56, 760–769. [Google Scholar] [CrossRef]
Aparicio, J.M.; Cardellach, E.; Rodríguez, H. Information content in reflected signals during GPS Radio Occultation observations. Atmos. Meas. Tech. 2018, 11, 1883–1900. [Google Scholar] [CrossRef] [Green Version]
Sievert, T.; Rasch, J.; Carlström, A.; Pettersson, M.I. Analysis of reflections in GNSS radio occultation measurements using the phase matching amplitude. Atmos. Meas. Tech. 2018, 11, 569–580. [Google Scholar] [CrossRef] [Green Version]
Bonnedal, M.; Christensen, J.; Carlström, A.; Berg, A. Metop-GRAS in-orbit instrument performance. GPS Solut. 2010, 14, 109. [Google Scholar] [CrossRef]
Sokolovskiy, S.; Schreiner, W.; Zeng, Z.; Hunt, D.; Lin, Y.C.; Kuo, Y.H. Observation, analysis, and modeling of deep radio occultation signals: Effects of tropospheric ducts and interfering signals. Radio Sci. 2014, 49, 954–970. [Google Scholar] [CrossRef]
Gorbunov, M.; Lauritsen, K.; Rhodin, A.; Tomassini, M.; Kornblueh, L. Radio holographic filtering, error estimation, and quality control of radio occultation data. J. Geophys. Res. Atmos. 2006, 111. [Google Scholar] [CrossRef] [Green Version]
Liu, H.; Kuo, Y.H.; Sokolovskiy, S.; Zou, X.; Zeng, Z.; Hsiao, L.F.; Ruston, B.C. A quality control procedure based on bending angle measurement uncertainty for radio occultation data assimilation in the tropical lower troposphere. J. Atmos. Ocean. Technol. 2018, 35, 2117–2131. [Google Scholar] [CrossRef]
Gorbunov, M.; Lauritsen, K.; Leroy, S. Application of Wigner distribution function for analysis of radio occultations. Radio Sci. 2010, 45. [Google Scholar] [CrossRef]
Gorbunov, M.; Lauritsen, K.; Leroy, S. Analysis of RO data retrieved from the Wigner distribution function. Presented at the 2nd Workshop, International Radio Occultation Working Group (IROWG), Estes Park, CO, USA, 28 March–3 April 2012. [Google Scholar]
Gorbunov, M.; Obukhov, O.; Kirchengast, G. Kirkwood Distribution Function and its Application for the Analysis of Radio Occultation Observations. Presented at the 7th Workshop, International Radio Occultation Working Group (IROWG), Elsinore, Denmark, 19–25 September 2019. [Google Scholar]
Gorbunov, M.E.; Kornblueh, L. Analysis and validation of GPS/MET radio occultation data. J. Geophys. Res. Atmos. 2001, 106, 17161–17169. [Google Scholar] [CrossRef]
Jensen, A.S.; Lohmann, M.S.; Nielsen, A.S.; Benzon, H.H. Geometrical optics phase matching of radio occultation signals. Radio Sci. 2004, 39. [Google Scholar] [CrossRef] [Green Version]
Sokolovskiy, S.; Rocken, C.; Hunt, D.; Schreiner, W.; Johnson, J.; Masters, D.; Esterhuizen, S. GPS profiling of the lower troposphere from space: Inversion and demodulation of the open-loop radio occultation signals. Geophys. Res. Lett. 2006, 33, 14. [Google Scholar] [CrossRef] [Green Version]
Jensen, A.S.; Lohmann, M.S.; Benzon, H.H.; Nielsen, A.S. Full spectrum inversion of radio occultation signals. Radio Sci. 2003, 38. [Google Scholar] [CrossRef]
Gorbunov, M. Canonical transform method for processing radio occultation data in the lower troposphere. Radio Sci. 2002, 37. [Google Scholar] [CrossRef]
Gorbunov, M.E.; Lauritsen, K.B. Canonical Transform Methods for Radio Occultation Data; Scientific Report 02-10; Danish Meteorological Institute: Copenhagen, Denmark, 2002; Available online: https://www.romsaf.org/Publications/reports/sr02-10.pdf (accessed on 24 February 2021).
Gorbunov, M.; Lauritsen, K. Analysis of wave fields by Fourier integral operators and their application for radio occultations. Radio Sci. 2004, 39. [Google Scholar] [CrossRef]
Gorbunov, M.E.; Gurvich, A.S.; Kornblueh, L. Comparative analysis of radioholographic methods of processing radio occultation data. Radio Sci. 2000, 35, 1025–1034. [Google Scholar] [CrossRef]

Figure 1. Overview of the geometry of an RO event.

Figure 2. Occultation event from UTC 0003 2012-09-11,

20.16^{\circ}

114.93^{\circ}

E. By applying e.g., Hanning windows of various lengths we can observe how the mapped coordinates are arranged. Smeared features along the frequency axis (left) or along the time axis (right) are instead smeared along the band around the range model, or along the range model itself. The center panel shows an appropriate window length for this occultation.

Figure 2. Occultation event from UTC 0003 2012-09-11,

20.16^{\circ}

114.93^{\circ}

Figure 3. Occultation event from UTC 0003 2012-09-11,

20.16^{\circ}

114.93^{\circ}

E. Spectrogram produced by STFT, where values ≤−50 dB are set to

- 50

dB. Left shows the spectrogram in time-frequency coordinates, right shows the spectrogram after mapping to BA and IH.

Figure 3. Occultation event from UTC 0003 2012-09-11,

20.16^{\circ}

114.93^{\circ}

E. Spectrogram produced by STFT, where values ≤−50 dB are set to

- 50

dB. Left shows the spectrogram in time-frequency coordinates, right shows the spectrogram after mapping to BA and IH.

Figure 4. Occultation event from UTC 0003 2012-09-11,

20.16^{\circ}

114.93^{\circ}

E. The dashed black lines show the various cross sections where we compare the amplitudes between STFT and SWPM.

Figure 4. Occultation event from UTC 0003 2012-09-11,

20.16^{\circ}

114.93^{\circ}

E. The dashed black lines show the various cross sections where we compare the amplitudes between STFT and SWPM.

Figure 5. Five horizontal cross sections of STFT (solid blue) and SWPM (dashed red) applied to the same measurement.

Figure 6. Five vertical cross sections of STFT (solid blue) and SWPM (dashed red) applied to the same measurement.

Figure 7. Occultation event from UTC 0003 2012-09-11,

20.16^{\circ}

114.93^{\circ}

E. An illustration of the effects of different window lengths. With shorter windows, uncertainty in IH grows (left). Longer windows cause more uncertainty in BA (right).

Figure 7. Occultation event from UTC 0003 2012-09-11,

20.16^{\circ}

114.93^{\circ}

E. An illustration of the effects of different window lengths. With shorter windows, uncertainty in IH grows (left). Longer windows cause more uncertainty in BA (right).

Figure 8. A MetOp-A occultation from UTC 0020 2015-02-10,

53.90^{\circ}

3.86^{\circ}

E. Left shows SWPM, right shows STFT. ECMWF BA is shown in solid black, standard PM is shown in dashed black. The window length for STFT is

1.5

Figure 8. A MetOp-A occultation from UTC 0020 2015-02-10,

53.90^{\circ}

3.86^{\circ}

E. Left shows SWPM, right shows STFT. ECMWF BA is shown in solid black, standard PM is shown in dashed black. The window length for STFT is

1.5

Figure 9. A MetOp-A occultation from UTC 0104 2015-02-10,

73.97^{\circ}

143.19^{\circ}

W. Left shows SWPM, right shows STFT. ECMWF BA is shown in solid black, standard PM is shown in dashed black. The window length for STFT is

1.5

Figure 9. A MetOp-A occultation from UTC 0104 2015-02-10,

73.97^{\circ}

143.19^{\circ}

W. Left shows SWPM, right shows STFT. ECMWF BA is shown in solid black, standard PM is shown in dashed black. The window length for STFT is

1.5

Figure 10. A MetOp-A occultation from UTC 0319 2015-02-10,

28.61^{\circ}

71.14^{\circ}

E. Left shows SWPM, right shows STFT. ECMWF BA is shown in solid black, standard PM is shown in dashed black. The window length for STFT is

1.9

Figure 10. A MetOp-A occultation from UTC 0319 2015-02-10,

28.61^{\circ}

71.14^{\circ}

E. Left shows SWPM, right shows STFT. ECMWF BA is shown in solid black, standard PM is shown in dashed black. The window length for STFT is

1.9

Figure 11. A MetOp-A occultation from UTC 0226 2015-02-10,

21.98^{\circ}

65.44^{\circ}

W. Left shows SWPM, right shows STFT. ECMWF BA is shown in solid black, standard PM is shown in dashed black. The window length for STFT is

2.0

Figure 11. A MetOp-A occultation from UTC 0226 2015-02-10,

21.98^{\circ}

65.44^{\circ}

W. Left shows SWPM, right shows STFT. ECMWF BA is shown in solid black, standard PM is shown in dashed black. The window length for STFT is

2.0

Figure 12. A MetOp-A occultation from UTC 0623 2015-02-10,

22.28^{\circ}

73.64^{\circ}

E. Left shows SWPM, right shows STFT. ECMWF BA is shown in solid black, standard PM is shown in dashed black. The window length for STFT is

2.0

Figure 12. A MetOp-A occultation from UTC 0623 2015-02-10,

22.28^{\circ}

73.64^{\circ}

E. Left shows SWPM, right shows STFT. ECMWF BA is shown in solid black, standard PM is shown in dashed black. The window length for STFT is

2.0

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

Share and Cite

MDPI and ACS Style

Sievert, T.; Rasch, J.; Carlström, A.; Ludwig Barbosa, V.; Pettersson, M.I.; Vu, V. Using A Sliding Window Phase Matching Method for Imaging of GNSS Radio Occultation Signals. Remote Sens. 2021, 13, 970. https://doi.org/10.3390/rs13050970

AMA Style

Sievert T, Rasch J, Carlström A, Ludwig Barbosa V, Pettersson MI, Vu V. Using A Sliding Window Phase Matching Method for Imaging of GNSS Radio Occultation Signals. Remote Sensing. 2021; 13(5):970. https://doi.org/10.3390/rs13050970

Chicago/Turabian Style

Sievert, Thomas, Joel Rasch, Anders Carlström, Vinícius Ludwig Barbosa, Mats I. Pettersson, and Viet Vu. 2021. "Using A Sliding Window Phase Matching Method for Imaging of GNSS Radio Occultation Signals" Remote Sensing 13, no. 5: 970. https://doi.org/10.3390/rs13050970

APA Style

Sievert, T., Rasch, J., Carlström, A., Ludwig Barbosa, V., Pettersson, M. I., & Vu, V. (2021). Using A Sliding Window Phase Matching Method for Imaging of GNSS Radio Occultation Signals. Remote Sensing, 13(5), 970. https://doi.org/10.3390/rs13050970

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Menu

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

Abstract

1. Introduction

2. Short-Time Fourier Transform

3. Sliding Window Phase Matching

3.1. Connection with STFT

3.2. Resolution

3.2.1. Resolution in Impact Parameter

3.2.2. Resolution in Bending Angle

3.2.3. Consequences for SWPM

4. Application of SWPM on Real RO Data

5. Conclusions

Author Contributions

Funding

Data Availability Statement

Conflicts of Interest

References

Share and Cite

Article Metrics

Article Access Statistics

Further Information

Guidelines

MDPI Initiatives

Follow MDPI