Improved GNSS-R Altimetry Methods: Theory and Experimental Demonstration Using Airborne Dual Frequency Data from the Microwave Interferometric Reflectometer (MIR)

Abstract

Altimetric performance of Global Navigation Satellite System - Reflectometry (GNSS-R) instruments depends on receiver’s bandwidth and signal-to-noise ratio (SNR). The altimetric delay is usually computed from the time difference between the peak of the direct signal waveform and the maximum of the derivative of the reflected signal waveform. Dual-frequency data gathered by the airborne Microwave Interferometric Reflectometer (MIR) in the Bass Strait, between Australia and Tasmania, suggest that this approach is only valid for flat surfaces and large bandwidth receivers. This work analyses different methods to compute the altimetric observables using GNSS-R. A proposed novel method, the Peak-to-Minimum of the 3rd Derivative (P-Min3D) for narrow-band codes (e.g., L1 C/A), and the Peak-to-Half Power (P-HP) for large bandwidth codes (e.g., L5 or E5a codes) show improved performance when using real data. Both methods are also compared to the Peak-to-Peak (P-P) and Peak-to-Maximum of the 1st Derivative (P-Max1D) methods. The key difference between these methods is the determination of the delay position in the reflected signal waveform in order to compute the altimetric observable. Airborne experimental results comparing the different methods, bands and GNSS-R processing techniques show that centimeter level accuracy can be achieved.

1. Introduction

The concept of what is called today interferometric GNSS-R (iGNSS-R) was proposed in 1993 for mesoscale altimetry [1]. It consist of the cross-correlation of the collected direct and reflected signals. After some encouraging results [2,3,4,5,6], the concept was refined in [7], and the technique was first demonstrated at the Zeeland bridge experiment in 2011 [8]. In the mid 90’s another GNSS-R technique was devised, which is called today conventional GNSS-R (cGNSS-R) [9] in which the reflected signals are cross-correlated with a locally generate replica of the transmitted signal.

In GNSS-R there are three main types of altimetry techniques: delay altimetry, phase altimetry, and the interference pattern technique. In delay altimetry, precision (${\sigma }_h$) is related to the delay precision ($sigm{a}_{tau}$), and in the presence of white thermal noise, the achievable delay precision is given by the Cramér Rao bound, which mainly depends on the signal-to-noise ratio (SNR), the instrument’s bandwidth (B) and the Gabor bandwidth ($\beta$) [10]. In phase altimetry, the phase of the complex waveform of the received signal is estimated, and if the GNSS-R receiver includes a phase-locked loop, then phase variations are dynamically compensated, and the delay can be readily estimated from the pseudo-ranges and a model of the scattering geometry. Phase altimetry is feasible if the scattering takes place over a specular surface, or at near grazing angles, which ultimately limits its range of applicability. Finally, the interference pattern technique can be used in ground-based receivers, and it is based on the observation of the fringes induced by the coherent interference of the direct and reflected signals.

Phase [11,12,13,14,15,16,17,18,19] and delay [20,21,22,23,24,25,26,27] altimetry performance studies for GNSS-R instruments have been published using spaceborne and airborne [28,29,30,31,32,33] data, or using the interference pattern between the direct and reflected signals [34,35,36,37]. Phase altimetry suffers from ambiguities difficult to correct for. So far, the most widely used technique to compute the altimetric delay the altimetric delay is typically based on the position of the maximum of the first derivative of the reflected signal waveform. The waveform is the cross-correlation of the direct or reflected signal with either a locally generated replica of the transmitted signal if using conventional GNSS-R (cGNSS-R) technique, or the cross-correlation of both the direct and the reflected signals with each other if using interferometric GNSS-R (iGNSS-R) technique [38]. The achievable performance has been studied in [10,28,39,40,41,42], including the retracking considerations [43,44,45,46], and the electromagnetic bias [47,48]. The relative delay estimation technique can be done using the maximum of the first derivative (P-Max1D), is well explained in [49], and a known disadvantage is that this position shifts depending on the receiver’s bandwidth [50,51].

In order to assess the performance of different GNSS-R altimetry methods, data from a MIR flight over the Bass Strait, which separates Australia and Tasmania, are used. MIR provides excellent data to fairly assess different GNSS-R methods, as the dual-frequency data are from different constellations (GPS and Galileo) sharing the same frequency bands, and can be processed with different parameters (the coherent and incoherent integration times), and processing techniques (cGNSS-R or iGNSS-R). In order to make a fair comparison, all the methods (P-Max1D, P-P, P-Min3D and P-HP) are evaluated using exactly the same MIR data-set. As shown in this paper, the delay of the maximum of the first derivative presents a number of limitations. To overcome them, a novel method for GNSS-R altimetry is introduced in this work: the Peak-to-Minimum of the third derivative (P-Min3D), specially useful for narrow-band signals as it is more resilient to the receiver’s finite bandwidth. For large bandwidth signals such as L5 or E5a, it is not clear which method would perform best as there is fewer literature. Therefore, all the methods are tested again for these signals, and it is found that the best performance is achieved by the Peak-to-Half power (P-HP).

2. Materials

2.1. Receiver Specifications

The Microwave Interferometric Reflectometer (MIR) instrument [52] (Figure 1) was designed, built and operated by the Universitat Politècnica de Catalunya. MIR is an airborne multi-constellation multi-beam dual-band GNSS-Reflectometer, with beam-steering capabilities to automatically track the direct and reflected GNSS signals from GPS and Galileo. MIR can acquire the direct and reflected signals at L1/E1 (1575.42 MHz), and at L5/E5a (1176.45 MHz) simultaneously using 2 + 2 up-looking beams, and 2 + 2 down-looking beams with a measured average bandwidth at −10 dB of ≈20 MHz at L1/E1, and ≈ 34 MHz at L5/E5a. The beams are automatically steered to track the direct and reflected GNSS signals, and compensate for aircraft attitude changes, and GNSS satellites movement. The signals are sampled at 32.736 MSps, with 1 bit for both the in-phase and quadrature (I and Q) components, and stored for offline processing [52], using either cGNSS-R or iGNSS-R techniques and configuring other parameters as the coherent and incoherent integration times.

2.2. Data

The data used in this analysis correspond to one of the transects of the MIR flight on 6 June 2018 over the Bass Strait. This transect has been selected because it tracks exactly the same GPS satellite simultaneously at L1 and L5 bands. This allows for a consistent comparison between bands. In addition, simultaneously to the tracking of the L1 and L5 bands, another beam was tracking a Galileo satellite at the E5a band. The aircraft flew at an average stable altitude of 1535 m with respect to the WGS84 ellipsoid. The ground track is shown in Figure 2. The transect started at 4:14 UTC at (38.9176° S, 149.0996° E), and ended at 4:36 UTC in (38.1843° S, 149.1952° E). The GPS satellite was PRN 3. It changed its azimuth angle from −44.5° to −46.3° with respect to the aircraft, and the incidence angle changed from 31.1° to 21.0°. The Galileo satellite was the PRN E08. In this case, the satellite changed its azimuth angle from −45.5° to −34.1° with respect to the plane, and the incidence angle changed from 33.7° to 27.8°.

The flight track and time were chosen so that it overlapped precisely with the track of CryoSat-2 [53] for comparison, starting at 4:11:29 UTC in (38.1846° S, 149.1982° E) and ending at 4:11:51 UTC in (38.9182° S, 149.1052° E) (Figure 2). CryoSat-2 altimeter offers centimeter level resolution, its data are freely available, so they can be used as ground truth of the closest MIR GNSS-R specular reflection point.

For this scenario the size of the first Fresnel zone at L5 is between 21 m and 33 m [54]. The plane speed was ∼67 m/s, which translates into an average distance between samples of ∼6.7 m for an integration time of 100 ms. For the L1 C/A signal, the resolution is limited by the width of its auto-correlation function, and the footprint is ∼300 m. In comparison, CryoSat-2 has a cross-track resolution of 15 km [53]. The along-track base resolution is also 15 km, although CryoSat-2 performs a SAR-processing technique, which significantly improves the resolution to ∼250 m in the along-track direction. CryoSat-2 high quality data clearly allows the observation of the shape of the geoid (Figure 3).

Equation (1) shows the cross-correlation to process the data in cGNSS-R; $y_s$ corresponds to either the direct or reflected signal in the time domain, and $y_r$ is the locally-generated PRN code replica of the transmitted signal. The term $e^{j2\pi f_{d}t}$ is a complex phasor at the Doppler frequency $f_d$; $T_{coh}$ is the coherent integration time; $Y(\tau ,f_d)$ is the so-called complex Delay-Doppler Map (DDM), and $|Y(\tau ,f_d){|}^2$ is the power DDM, which can be incoherently integrated to reduce speckle noise (Equation (2)). In Equation (2), $N_{incoh}$ is the number of correlations incoherently averaged corresponding to an incoherent integration time ($T_{incoh}$, Equation (3)). The cross-correlations are computed using Fast Fourier Transforms (FFT) as described in [55]. In this work, the GNSS-R observable is the waveform or delay profile at the particular Doppler frequency that passes through the maximum of the amplitude of the DDM, $|Y(\tau ,f_d=f_d{|}_{{\left|Y\right|}_{max}^2}){|}^2$. The computation time is limited, so in this study, L1 and L5 waveforms which come from the same satellite are prioritized and computed using the cGNSS-R technique with coherent integration times of 1, 4 and 10 ms for GPS data. While Galileo E5a data is only computed with 1 ms. Using the iGNSS-R technique the coherent integration time is also set to 1 ms. In all cases the data are incoherently averaged up to 100 ms, as shown in

Table 1. Summary of the MIR data processed in this work. Available bandwidth (BW) refers to the bandwidth limited by the public codes in cGNSS-R cases or by the total transmitted bandwidth in the iGNSS-R cases [10]. MIR has a measured average bandwidth at −10 dB of ≈20 MHz at L1/E1, and ≈34 MHz at L5/E5a.

Sat.	Band	Available	MIR	Process.	${\mathit{T}}_{\mathbf{coh}}$	${\mathit{N}}_{\mathbf{incoh}}$	${\mathit{T}}_{\mathbf{incoh}}$	Average SNR
PRN		BW [MHz]	beam	Tech.	[ms]		[ms]	[dB]
3	L1	2.046	2	cGNSS-R	1	100	100	20.5
3	L1	2.046	2	cGNSS-R	4	25	100	26.5
3	L1	2.046	2	cGNSS-R	10	10	100	28.7
3	L1	30.69	2	iGNSS-R	1	100	100	10.9
3	L5	20.46	3	cGNSS-R	1	100	100	22.6
3	L5	20.46	3	cGNSS-R	4	25	100	25.3
3	L5	20.46	3	cGNSS-R	10	10	100	27.5
3	L5	24.00	3	iGNSS-R	1	100	100	18.8
E08	E5a	20.46	4	cGNSS-R	1	100	100	20.8

.

$$|Y(\tau ,f_d){|}^2={|{\int }_0^{T_{coh}}{y}_r\left(t\right)y_s^{*}(t-\tau )e^{j2\pi f_{d}t}dt|}^2$$

(1)

$$\overline{|Y(\tau ,f_d){|}^2}=\frac{1}{N_{incoh}}\sum _{n=1}^{N_{incoh}}{\left|Y(\tau ,f_d)\right|}_n^2$$

(2)

$$T_{incoh}=N_{incoh}·T_{coh}$$

(3)

The SNRs are computed as the ratio of the waveform peak and the noise floor level, estimated using 30 delay bins prior to the start of the waveform, where the signal power is not present. cGNSS-R uses a locally-generated replica of the code and the SNR is higher than in iGNSS-R in which the correlation is computed with the direct noisy signal. However, it has the advantage that it is not limited to public codes and their limited bandwidths (Table 1). Increasing the coherent integration time increases the SNR [56]. Comparing GPS L1 and L5, which are transmitted from the same satellite, L1 signals with $T_{coh}=1$ ms have lower SNR than L5, as the transmitted power at L1 C/A is lower than at L5 [10].

3. Altimetric Methods

The waveforms contain information on the time of arrival of the signal. For airborne heights, the flat Earth approximation holds, and the height measurement ($h_{measured}$) can be computed from the delay difference between the direct ($T_{direct}$) and the reflected ($T_{reflected}$) waveforms as in Equation (4) (from p.171 of [57]), where c is the speed of light, and ${\theta }_{elev}$ is the elevation angle, which is the complementary of the incidence angle (${\theta }_{inc}$) (Figure 4). In the interferometric case, the delay can be derived from the interferometric waveform [58]. Equation (4) can be transformed into Equation (5) to use the difference in samples between waveforms ($dif{f}_{samp}$), where $Fs$ is the sampling frequency. Finally, subtracting the height of the platform ($h_{planeWGS84}$) referenced to the WGS84 ellipsoid, the sea surface height ($SS{H}_{measured}$) referenced to the WGS84 is obtained (Equation (6)). Note that, because of the low flight height, atmospheric and ionospheric delay corrections can be neglected. Variations in the pitch (${\sigma }_{pitch}=0.46$° ), roll (${\sigma }_{roll}=0.65$° ) and yaw (${\sigma }_{yaw}=0.69$° ) angles from the plane are also negligible during the integration time.

$$h_{measured}=\frac{T_{reflected}-T_{direct}}{2sin\left({\theta }_{elev}\right)}·c$$

(4)

$$h_{measured}=\frac{dif{f}_{samp}}{2cos\left({\theta }_{inc}\right)}·\frac{c}{Fs}$$

(5)

$$SS{H}_{measured}=h_{planeWGS84}-h_{measured}$$

(6)

One of the key points in GNSS-R altimetry is the determination of the tracking point of the reflected waveform to compute the delay difference (Figure 5). In this work different approaches are analyzed:

• In Peak-to-Peak (P-P), the peak (maximum) of the waveform is optimal for the direct signal. However, for the reflected case, it will only be optimal if the reflection is specular. As seen in [54], during the analyzed flight the wind-driven waves were ∼1.7 m high with a period of 5 s, and the swell waves were ∼1.2 m high with a period of 10 s. Therefore, the flat surface condition (Rayleigh criterion, Equation (7), [35,59]) is not satisfied for the frequencies and incidence angles considered, and the reflection was not specular.

  $${\sigma }_h<\frac{\lambda }{8cos\left({\theta }_{inc}\right)}$$

  (7)
• The second method considered uses the peak of the direct signal waveform, and the maximum of the first derivative of the reflected signal waveform (P-Max1D) [49].
• The third method tracks the Half-power position of the waveform (P-HP), as in radar altimetry [60]. It is adapted to GNSS-R using the peak of the direct signal waveform and the half-power point of the leading-edge of the reflected signal waveform. The main difference with respect to radar altimetry is that the bandwidths used in GNSS-R are much smaller, at least by an order of magnitude. Based on the model of the position of the maximum of the first derivative [49], some GNSS-R experiments [21,41] have also used this technique, but using the point at 70% instead.
• Finally, to overcome the limitations of the previous methods, a new method is proposed: the use of the peak of the direct signal waveforms and the minimum of the 3rd derivative of the reflected signal waveform (P-Min3D).

3.1. GNSS-R Altimetry Using the Peak-to-Minimum of the 3rd Derivative

This novel method can be understood as a refinement of the analysis made in [49]. In Case 1 the simulation in [49] is repeated to study the minimum of the 3rd derivative, and compared it with the maximum of the first derivative. Then, in Case 2, the consideration of having a rougher sea altered by waves is included in the simulation to make it more realistic.

3.1.1. Case 1

Following the procedure in [49], the ideal L1 C/A waveform (Figure 6a), and an off-specular (or scattered) reflection impulse response (Figure 6b) are generated first.

The reflected waveform (Figure 6c) is generated as the convolution of the previous two functions (Figure 6a,b). As expected, the maximum of the first derivative corresponds to the inflexion point (zero of the second derivative) at delay 0, and so it is the minimum of the third derivative.

To include the effect of the receiver’s finite bandwidth, the waveform has to be smoothed (low-pass filtered). In doing so, the maximum of the first derivative and the minimum of the third derivative shifts towards negative delays in the non-reflected case (Figure 7a, bandwidth ≈ 5.115 MHz), as compared to the ideal case (Figure 6a, infinite bandwidth).

When convolved with the reflected power function (Figure 7b) the first derivative shifts towards a negative delay [49,50]. The zero in the second derivative shifts with the maximum in the first derivative, while the minimum shifts to the positive side, closer to the peak of the waveform, at the point of maximum curvature. Between these two points, the second derivative describes an abrupt slope change because the waveform’s curvature changes drastically. The third derivative tracks the rate of change of the curvature, which has a minimum when the slope of the second derivative is minimum (maximum negative change). It has to be noted that the maximum negative rate of change in the waveform curvature (minimum of the 3rd derivative) occurs at delay 0, in both the ideal case and in the smoothed one (finite bandwidth case).

The position change of the maximum of the first derivative, and the minimum of the third derivative, is now computed for several bandwidths (i.e., low-pass filtered versions). Figure 8 shows that the position of the first derivative maximum is only optimum for very large bandwidths (near infinite), otherwise a non-negligible shift in the tracking position of the waveform is introduced. Conversely, the third derivative minimum correctly tracks the delay for a wider range of bandwidths (larger than ∼5 MHz), and then the shift is less abrupt.

3.1.2. Case 2

In the previous case the leading edge of the reflected power function had an infinite slope. However, when the reflection takes place over a rough sea (i.e., waves) as in radar altimetry the leading edge of the returned power has a finite slope [60].

To analyze this effect the increase of the rise time of the leading edge with the significant wave height (SWH) is modelled as $2·SWH/c$ [61]. Considering the swell waves of 1.7 m, the rise time is not a step function anymore at delay 0, but takes about 11.3 ns, or 0.0116 C/A chips. Since the power function is normalized, the slope can be estimated from the rise time delay as $1/0.0116=86.2$ u/(C/A chips). The maximum point of the triangle (rise time, as it starts at 0) is 0.0116 C/A chips, and the mid-rise position of the leading-edge is at 0.0058 C/A chips (Figure 9a). Following the same procedure, Figure 9a is convoluted with the ideal C/A waveform (Figure 6a). Then, as in the previous case, low-pass filtered cases are computed to simulate a bandwidth reduction. Figure 9c shows the simulated case for a bandwidth ≈ 5.115 MHz.

As in Figure 8, the tracking point is computed as a function of the equivalent receiver’s bandwidth for the maximum of the first derivative, and the minimum of the third derivative (Figure 10a). The position of the delay of the minimum of the third derivative is coincident with the mid-rise position of the leading-edge of the reflected power function (marked with a dashed green line) for a wide range of bandwidths (larger than ∼5 MHz), which is optimal to compute the altimetry delay. For the first derivative, the tracking position is close to the zero delay only for very large bandwidths (near infinite), but starts to drift rapidly as in Case 1 (Figure 8).

In order to generalize the previous results, and include other factors like the incidence angle, a smaller leading edge slope is simulated for the power reflection function. The slope considered is 60 u/(C/A chips) which models the reflection over a rough sea surface with a SWH ≈2.4 m. The maximum for this is at about 0.0167 C/A chips, and the mid-rise position of the leading-edge is at 0.0083 C/A chips.

For this slope (Figure 10b), the behavior is similar to the previous case. The minimum of the third derivative still tracks the mid-rise position of the leading-edge of the reflected power function. Other lower slopes have been tested (e.g., 20 u/(C/A chips), SWH ≈ 7.3 m) and the same behavior has been observed.

4. Results

All the described altimetry methods are computed for the time series of all the waveforms computed with $T_{coh}=1$ ms and $T_{incoh}=100$ ms. For each different method, a linear regression with the CryoSat-2 data is computed. Once the best method is found (slope of the linear fit closer to 1, as offsets can be calibrated) processing for other $T_{coh}$ at the same frequency band and processing type are conducted (Table 1).

For each configuration (Table 1), the Allan Variance [62] and the unbiased root mean squared altimetry error (ubRMSE) with respect to the CryoSat-2 data are computed for different bands and processing methods.

The computation of the Allan Variance requires that the altimetric measurement series are detrended so that they have a zero mean. For the computation of the ubRMSE, the offset of the linear fit regression is subtracted from all the values. Missing data samples are omitted in the computations, and values of the track edges closer to the moving window length discarded.

4.1. Altimetric Methods

Figure 11 shows the SSH obtained from L1 computed using the cGNSS-R technique, and 3 different methods for the height estimation. The linear fit of the estimates is stable with increasing averaging. Therefore, to have a better visualisation, the figure shows the incoherently averaged data at 100 s. In Figure 11a, SSH estimates are plotted against time. In Figure 11b the unbiased SSH estimates from MIR are plotted with respect to the CryoSat-2 SSH data. Choosing a sub-optimal method not only introduces an offset, but it also affects the slope of the regression, which ideally should be equal to one.

Table 2. Slopes and offsets of the regressions for the different GNSS-R altimetry methods, bands, and processing techniques with respect to the CryoSat-2 data. Using $T_{coh}=1$ ms and $T_{incoh}=100$ ms.

Process. Tech.	Band	Altimetry Method	Slope [m/m]	Offset [m]
cGNSS-R	L1	P-P	1.36	−9.42
	L1	P-Max1D	0.71	12.25
	L1	P-HP	0.57	46.52
	L1	P-Min3D	1.03	2.78
cGNSS-R	L5	P-P	1.35	−8.63
	L5	P-Max1D	1.12	1.97
	L5	P-HP	1.03	3.87
	L5	P-Min3D	1.17	−0.97
cGNSS-R	E5a	P-P	1.20	−8.39
	E5a	P-Max1D	1.05	1.52
	E5a	P-HP	1.02	3.04
	E5a	P-Min3D	1.07	−0.57
iGNSS-R	L1	P-P	1.14	−4.06
	L1	P-Max1D	1.17	1.48
	L1	P-HP	0.50	11.17
	L1	P-Min3D	1.16	1.01
iGNSS-R	L5	P-P	1.20	−6.84
	L5	P-Max1D	1.19	1.73
	L5	P-HP	1.05	4.09
	L5	P-Min3D	1.16	0.63

shows the results of the fit of the GNSS-R altimetry products with CryosSat-2 for the different altimetry methods, using L1, L5 and E5a, cGNSS-R and iGNSS-R techniques, for $T_{coh}=1$ ms and $T_{incoh}=100$ ms. The following conclusions can be extracted:

• In the first part of the table L1 cGNSS-R cases are shown.

  -

  The P-P method is not optimal as the slope is significantly larger than 1.

  -

  The P-Max1D is not optimal either. As explained in Section 3.1, results match with the simulation of P-Min3D, where, due to the finite bandwidth, the tracking position shifts.

  -

  The P-HP shows the worst results. L1 C/A signal limited bandwidth (Table 1) makes a wide waveform with a gradual leading edge. Figure 9c) shows that the 50% position is shifted almost 0.2 C/A chips (∼58 m) from the optimal zero delay, which explains the large offset and that the behaviorfrom the estimates is so far away from the optimal (slope 1). As it will be seen, this method extrapolated from radar altimetry works better with larger bandwidth signals with a steeper leading edge.

  -

  As in the simulations, the P-Min3D tracks well the desired mid-rise leading-edge position of the reflected power function, and has a slope very close to 1. L1 MIR waveforms (Figure 12a) have similar shape to the simulated cases (Figure 9c).

• The next section of the Table 2 shows the main results for the cGNSS-R at L5:

  -

  As in L1, the P-P method does not offer the best performance. As it has been stated, this method works best with specular reflections and the sea waves present on that day induce a non-specular reflection.

  -

  The P-Max1D performs slightly better than at L1. In the L5 case, the first derivative maximum is close to the half power position (Figure 12b), because of the steepness of the waveform.

  -

  The P-HP achieves the best slope (i.e., closest to 1). L5 codes have a large bandwidth, 10 times larger than L1 C/A codes. which leads to a sharper waveform (Figure 13).

  -

  Finally, the P-Min3D presents a similar behaviorto the P-Max1D. Its position is closer to the half power position as now the leading edge is much steeper, due to the larger bandwidth.

• The slopes of the E5a waveforms have a similar behaviorto those of at L5. The waveform shapes are very similar (Figure 13b,c), and the code bandwidth is the same. The P-HP method is the one with the closest slope to 1. Galileo E5a transmitted power is also 3 dB higher than GPS L5 [10], this is likely to be the reason why, in general, all the slopes are better than at L5.
• As shown in Table 1, iGNSS-R cases exhibit a lower SNR than cGNSS-R, and lower at L1 than at L5. L1 waveforms also present multiple correlation peaks overlapping in the same composite waveform, formed by the correlation of the different signals present at L1 (e.g., C/A, P(Y) or M-codes). The low SNR, and the different shapes present in the waveforms are actually a challenge for these methods. The P-P method is the one performing the best for this case. It usually tracked the peak of the correlation of the encrypted codes, although not in the optimal position for non-specular reflections. For these reasons, the slope of the fit is worse than for the other cases.
• The iGNSS-R L5 case is quite similar to the conventional case, as it is already using large bandwidth codes. The secondary peaks from secondary reflections which sometimes appear in the L5 waveforms [54] are also present. As in the L5 conventional case the P-HP achieves the slope closest to 1.

All things considered, to get the best of iGNSS-R processing, different methods than those used for cGNSS-R are needed to properly handle multiple reflections, correlation peaks, or possibly different signals transmitted by the GNSS satellites.

Offsets theoretically should tend towards 0, but they will normally differ a bit for different factors (e.g., antenna phase centers are not in the same exact point). However, it is important to notice that a wrong method with a slope different than 1 also introduces error in the offset (e.g., cGNSS-R L1 P-HP). Once the slope is correct, this can be easily calibrated as it is constant.

4.2. Allan Variance

The computed Allan variance is shown in Figure 14. In all cases, the general decreasing trend of the variance is proportional to $\sqrt{T_{incoh}}$. The following conclusions can be extracted:

• L1 cGNSS-R with $T_{coh}=1$ ms measurements with few averaged samples have less variance than those from the L5. This might be due to the oversampling at L1. MIR uses 32 samples per C/A chip, at L5 due to the narrower waveforms the number of samples defining the peak is smaller. As the averaging increases, L5 starts to show a smaller variance than the L1 ones, although both results are quite close.
• L1 and L5 cGNSS-R at long coherent integration times, $T_{coh}=4$ ms and $T_{coh}=10$ ms exhibit better performance, at L5 better than L1. At L1 cGNSS-R the SNR improved with longer coherent integration times, and the variance decreased as the SNR increased. At L5, the SNR also increased with longer coherent integration times, but for $T_{coh}=4$ ms the variance is smaller than for $T_{coh}=10$ ms. This behavioris attributed to the coherence time of a surface (${\tau }_s$) which can be estimated as Equation (8) [56].

  $${\tau }_s=\frac{\lambda }{2·v}\sqrt{\frac{h}{2·c·{\tau }_c·cos{\theta }_{inc}}}$$

  (8)

  where $\lambda$ is the wavelength; v the platform speed; h is the height; c the speed of light; ${\tau }_c$ is the duration of a chip; and ${\theta }_{inc}$ the incidence angle. For the conditions of this flight, at L1, the coherency would be lost after 2 ms for ${\theta }_{inc}=0$° or 2.8 ms for ${\theta }_{inc}=45$° [63], so the improvement in 4 and 10 ms is purely linked to the SNR improvement. However, at L5 coherency would be lost after 7.7 ms for ${\theta }_{inc}=0$° or 10 ms for ${\theta }_{inc}=45$° [63]. Since ${\theta }_{inc}$ varied from 31° to 21°, it never reached 10 ms, but the coherence of the surface was always longer than for 4 ms. This is clearly shown by the L5 variance which reached the better value for $T_{coh}=4$ ms.
• The iGNSS-R cases achieved the smallest variances since they have largest bandwidths, despite the lower SNRs. In this particular case L1 is a bit better than L5. iGNSS-R cases are limited by the receiver’s bandwidth and the maximum bandwidths transmitted at the band. In L5 MIR bandwidth ∼34 MHz is larger than the total transmitted bandwidth at L5, 24 MHz, so noise is being added to the signal. At L1 MIR bandwidth ∼20 MHz is lower than the 30.69 MHz transmitted, in this case the maximum resolution [10] is limited by the bandwidth, but there is no addition of noise due to a larger bandwidth.
• E5a measurements (cGNSS-R with $T_{coh}=1$ ms), start with a variance similar to L1, and as the averaging increases it shows similar variances to the L1 and L5 cGNSS-R ones with $T_{coh}=1$ ms.

After ∼40 to 90 samples of averaging the variance stops decreasing. This corresponds to remaining differences between the aircraft height variations and the measured height. Especially in the iGNSS-R cases, that the variance also stopped decreasing at ∼5 to 8 samples. The equivalent distance of these amount of samples corresponds to ∼30 m to ∼50 m, which corresponds to the averaging distance of swell waves.

4.3. ubRMSE

Figure 15 shows the unbiased root mean squared error. At $T_{incoh}=100$ ms (1 sample, 6.7 m) the worst ubRMSE is ∼6.6 m for L5 cGNSS-R with $T_{coh}=1$ ms case, and the best one is ∼2.7 m for the iGNSS-R case. At about 10 s integration (100 samples, 670 m) for all cases it is about 80 cm–1.3 m, and at 100 s, between ∼30–40 cm, except for L1 iGNSS-R case which is ∼60 cm. L1 iGNSS-R has the slope fit with CryoSat-2 data most different from 1 among the considered cases, although having the lowest Allan variance for $\tau$[average samples] > 300 s.

Coincident with the minimum of Allan variance, most cGNSS-R cases reached the minima from ∼7 cm to ∼17 cm at around 500–550 s. Except, L1 iGNSS-R which only reached ∼30 cm. The drawback of this long averaging times is the coarse spatial resolution (∼37 km).

5. Discussion

The analysis of the slope of the regression for each method (P-Max1D, P-P, P-Min3D and P-HP) has shown that for L1 cGNSS-R, the best method is our proposed P-Min3D (Table 2). Unlike P-P and P-HP, this method computes the delay estimate in the waveform near the optimal place as P-Max1D, but it has shown to be more resilient to the receiver’s finite bandwidth. For waveforms derived from larger bandwidth signals, such as the codes at L5 or E5a, the methods that had performed best at L1 (P-Max1D and P-Min3D) may work, but are not optimal as the waveform is much narrower. The P-P does not perform well as the reflections is not specular. P-HP method performs the best for these cases, thanks to the larger bandwidth of the signal.

Moreover, analyzing the magnitude and behavior from the precision results (

Table 3. Unbiased root mean squared error (ubRMSE, ${\sigma }_h$), with respect CryoSat-2, for the computed bands, altimetric methods and coherent times. For $T_{incoh}=$ 0.1, 1 and 10 s. The slopes are computed from the fit with respect to the CryoSat-2 data using $T_{incoh}=100$ ms, as in Table 2.

Process.	Band	Altimetry	${\mathit{T}}_{\mathbf{coh}}$	Slope	${\mathit{\sigma }}_{\mathit{h}}$ [m]	${\mathit{\sigma }}_{\mathit{h}}$ [m]	${\mathit{\sigma }}_{\mathit{h}}$ [m]
Tech.		Method	[ms]		at 0.1 s	at 1 s	at 10 s
cGNSS-R	L1	P-Min3D	1	1.03	6.11	2.92	1.35
cGNSS-R	L1	P-Min3D	4	1.06	4.78	2.37	1.22
cGNSS-R	L1	P-Min3D	10	0.99	4.74	2.34	1.14
iGNSS-R	L1	P-P	1	1.14	2.87	1.89	1.10
cGNSS-R	L5	P-HP	1	1.03	6.53	2.78	1.14
cGNSS-R	L5	P-HP	4	1.05	4.03	2.07	0.97
cGNSS-R	L5	P-HP	10	1.04	4.20	2.09	0.92
iGNSS-R	L5	P-HP	1	1.05	2.75	1.67	0.85
cGNSS-R	E5a	P-HP	1	1.02	5.99	2.79	1.25

) match with those found in other research of airborne experiments [31]. For 1 s (67 m) incoherent averaging, the obtained ubRMSE is of ∼2–3 m, ∼1 m for 10 s (670 m) incoherent averaging, and up to centimeter level, 7–17 cm, with longer averages.

This means that when the same data is used to compare between the different altimetry methods the new ones considered (P-Min3D and P-HP) perform the best and, in addition, when these are compared with other airborne experiments they show consistency in the measurements.

Further analyzing the ubRMSE of MIR SSH data when compared to CryoSat-2 SSH data:

• iGNSS-R cases achieved the best results showing an ubRMSE improvement of a factor of ∼2.2 with respect to the same band and $T_{coh}$ using the cGNSS-R technique. From the simulated precisions for GEROS-ISS [42], more improvement was predicted when using the iGNSS-R technique as compared to the cGNSS-R technique. However, this was a space-borne prediction, and considering dual-frequency observables, among many other corrections, which makes the results hardly comparable. Differently, it was experimentally shown [32] that the improvement factor when using P(Y) with the reconstructed GNSS-R (rGNSS-R) technique compared to C/A with the cGNSS-R technique was ∼2.4. Considering that rGNSS-R will have a higher SNR than iGNSS-R it matches that the improvement factor is higher for rGNSS-R, but with very similar magnitudes.
• Comparing L1 and L5 iGNSS-R cases, they performed similar with low averages, but L5 better overall. L5 iGNSS-R signal had almost 8 dB more SNR than L1 iGNSS-R, and MIR bandwidth is more limited at L1 than at L5. Therefore, the expected achievable resolution was lower at L1 than at L5, and the results show this behavior. Theoretical optima have been analyzed in [10], showing that L1 should perform better than L5. However, the conditions considered for the theoretical analysis are distinct than the ones provided by MIR, as in that analysis the signals are considered to have the optimal bandwidth and same 20 dB SNR, while MIR signals have different bandwidths than the theoretical optimal and SNR varies with the signal. For these reasons both iGNSS-R cases have these performances.
• A similar thing happens with the E5a signal. It does perform the best for cGNSS-R with $T_{coh}=1$ ms and $T_{incoh}=100$ ms, but when increasing the incoherent averaging it does not have a clear advantage compared with L5. According to the GEROS-ISS simulations [42] Galileo signal should perform the best, but again it is considering dual-frequency observables. And looking at the optimal parameters again [10], the most optimal parameters are not met. Also, the signals come form different satellites which may generate significant differences for its fair comparison. Between the L1 and L5 cGNSS-R cases, L5 performs better than L1, because L1 only uses the C/A code, which has a narrower bandwidth compared to L5 code.
• Increasing the coherent integration time in cGNSS-R cases increases the SNRs, which improves the ubRMSE. An improvement of a factor of ∼1.3–1.6 is observed with respect to the $T_{coh}=1$ ms. The best performance with coherent integration time is observed when the coherency of the surface is considered, as at L5.

6. Conclusions

This work has analyzed different methods for GNSS-R altimetry, and compared them using dual frequency data from a flight of the MIR instrument over the Bass Strait. The proposed P-Min3D for L1 C/A signals, and P-HP for larger bandwidth signals (e.i. L5/E5a) showed the best altimetry performance.

For 10 s of incoherent integration time, an ubRMSE with respect to CryoSat-2 data of ∼0.85 m to ∼1.35 m was achieved. The minimum ubRMSE reached ranged from ∼7 cm to ∼17 cm for all cases considering averages ∼550 s which is consistent with other airborne experiments [31], but spatial resolution is very coarse for these long averages. Only L1 iGNSS-R which have slightly more error. L5 and E5a signals showed better performance than L1 C/A signals, which have narrower bandwidth. The iGNSS-R performance was ∼2.2x times better than cGNSS-R showing the potential of this technique. This improvement coincides with the one found in [32]. Increasing the coherent integration time in cGNSS-R cases achieved an improvement of factor up to ∼1.6x times. Future work will explore the feasibility of using higher waveform derivatives to better track the altimetric delay, the impact of noise amplification and the reliability of the methods for space-borne missions.

Despite, the methods presented have been derived from airborne data, they can also be applied for Low-Earth Orbit (LEO) data as well with the proper geometry modifications to account for Earth’s curvature, and in single frequency GNSS-R instruments using ionospheric models (e.g., Klobuchar or NeQuick) to compensate for the ionospheric delay.