Normalized GNSS Interference Pattern Technique for Altimetry

It is well known that reflected signals from Global Navigation Satellite Systems (GNSS) can be used for altimetry applications, such as monitoring of water levels and determining snow height. Due to the interference of these reflected signals and the motion of satellites in space, the signal-to-noise ratio (SNR) measured at the receiver slowly oscillates. The oscillation rate is proportional to the change in the propagation path difference between the direct and reflected signals, which depends on the satellite elevation angle. Assuming a known receiver position, it is possible to compute the distance between the antenna and the surface of reflection from the measured oscillation rate. This technique is usually known as the interference pattern technique (IPT). In this paper, we propose to normalize the measurements in order to derive an alternative model of the SNR variations. From this model, we define a maximum likelihood estimate of the antenna height that reduces the estimation time to a fraction of one period of the SNR variation. We also derive the Cramér–Rao lower bound for the IPT and use it to assess the sensitivity of different parameters to the estimation of the antenna height. Finally, we propose an experimental framework, and we use it to assess our approach with real GPS L1 C/A signals.


Introduction
Global Navigation Satellite Systems reflectometry (GNSS-R) is a well-established method for remotely sensing many relevant geophysical properties of the reflection surfaces. GNSS-R was first proposed within the frame of the PAssive Reflectometry Interferometric System (PARIS) project as a bistatic radar remote sensing technique for ocean altimetry using the L-band GPS signal [1]. Since then, GNSS-R has been demonstrated to be useful in other applications, such as monitoring water levels and snow height with a ground approach [2][3][4]. In this approach, the antenna, situated on a mast, receives a direct GNSS signal coming from the satellite and a nadir signal reflected by the observed surface. Assuming that the antenna position is known, we can compute the position of the surface of reflection.
This approach provides precise localization and dating of the measures that allows for spatio-temporal comparison of water levels and snow cover, respectively [5][6][7][8]. These parameters are very important for flood monitoring and avalanche prevention, as well as for hydroelectric companies. Furthermore, the approach is noninvasive and can be easily implemented on a portable instrument and embedded in a vehicle with a mast [9].
GNSS-R altimetry can be carried out in two different ways, depending on the ranging principle: code altimetry and phase altimetry [10]. With code altimetry, only the GNSS code delay difference between the direct and reflected signals is used. With phase altimetry, the phase of the signal is also used for computing this delay difference [9]. The interference pattern technique (IPT) considers the behavior of the signal-to-noise ratio (SNR) of the received GNSS signal as a function of the satellite elevation [2,11]. The direct and reflected GNSS signals are combined at the antenna. Due to their different phase variations, the SNR oscillates at a rate proportional to the distance between the antenna and the surface of specular reflection. Unlike satellite or airborne reflection scenarios, ground GNSS receivers observe a coherent interference pattern if we consider a flat reflecting surface (compared to the carrier wavelength) on an area corresponding to the first Fresnel zone. A few previous works can be found analyzing the accuracy of these GNSS-R altimetry techniques [12][13][14][15]. Initial works proposed simple analytical models to describe the altimetry precision as a function of system/instrument parameters [12]. These methods rely on a considerable number of assumptions that might hold only for some specific scenarios and provide a pessimistic bound on the achievable precision. In [13], the authors proposed a Cramér-Rao lower bound (CRLB) closed expression for code altimetry. The CRLB is a statistical tool that provides a lower bound on the achievable estimation error for any unbiased estimator. A new derivation was proposed in [15] to compute the CRLB for code altimetry and a specific set of measurement data under multiple SNR scenarios.
Unfortunately, one of the main drawbacks of the IPT is that very long measurement times are usually needed. The observed SNR oscillates with the variation of the satellite elevation, but satellite elevation varies slowly, thus long measurement times are required to estimate the SNR frequency of oscillation. To reduce the estimation time to a fraction of one period of the SNR variations, we propose an alternative model for the measured SNR. This normalized model is based on the normalization of the measured signal amplitudes and is possible only after an initial calibration step. This calibration step consists of varying the antenna height of the receiver a value dh in order to obtain the minimum and maximum value of SNR for a given satellite elevation. Using the normalized model, we define a maximum likelihood estimate of the antenna height that allows for the reduction of the required estimation time to a fraction of one period of the SNR variation. We also derive the minimum antenna variation range dh as a function of the satellite elevation and deduce from this function the minimum observation time as a function of the satellite elevation rate. In addition, we derive the CRLB for the IPT and use it to assess the sensitivity of different parameters to the estimation of the antenna height. Finally, we propose a novel experimental framework, which we use to assess our approach with real signals.
This paper is organized as follows: Section 2 describes the interference pattern problem. The considered signal model is introduced as a function of the receiver height. The proposed estimator is presented in Section 3. Section 4 includes the derivation of the CRLB for the IPT and the proposed estimator performance assessment using synthetic signals. In Section 5, the proposed experimental framework is described, and the results obtained with real GPS L1 C/A signals within this framework are presented. Finally, Section 6 summarizes the paper, highlighting its main conclusions.

Interference Pattern Problem
We show in Figure 1 the reflectometry principle for an antenna situated on a mast of height h. In our approach, we take into account the subset of satellites M that have a specular reflection on the surface to analyze. The antenna receives M scaled, time-delayed and Doppler-shifted signals with known signal structures. Each signal corresponds to the line-of-sight or direct signal (s D i ) plus its corresponding specular reflection (s R i ). The overall received signal can be modeled as: where f RF is the carrier frequency, f d i the Doppler frequency shift of the i-th satellite, ϕ D i the receiver clock phase offset, ϕ R i (t) the phase delay between the direct and reflected signals as a function of time, τ D i , τ R i the time-delays, c i (t) the pseudorandom code sequence, A D i and A R i the amplitudes of the direct and reflected received signals and n(t) zero-mean additive Gaussian noise with variance σ 2 n . In this paper, we will assume that the variations of A D i and A R i are negligible during the short periods of observation considered, e.g., a few minutes. For long observation times, A D i and A R i will change with time as a function of the satellite elevation due to the variation of the received power and the antenna footprint. The time dependence of τ D i , τ R i , f d i and ϕ D i has been neglected for simplicity in Equation (1), since their variation over time will be compensated for by the receiver's tracking stage with little impact on the proposed analysis.

GNSS th satellite
In our current study, we will consider only the processing of the GPS L1 C/A signals. In this case, for an antenna mounted on a mast a few meters above the reflecting surface, the difference in GNSS signal path between the direct and the reflected signals will be small compared to the chip duration of the code. Thus, we can assume that According to the geometry depicted in Figure 1, it is easy to show that: where θ el i (t) is the elevation of the i-th satellite at instant t, λ L1 = 19.042 cm is the GPS L1 C/A carrier wavelength, h is the height of the antenna and c represents the speed of light. In the following, for notation simplicity, we will drop the satellite index i, since we will focus on the processing of the direct and reflected signals coming from a single satellite. In this case, by using some trigonometry, we can express x(t) as: where: and A G (t) represents the magnitude of the composite signal, while ϕ G (t) represents its phase. In practice, in a GNSS receiver, the SNR is estimated after the signal down-conversion and correlation with a local code replica. In this case, the SNR is proportional to A 2 G (t), the squared amplitude of the received signal. In this context, we see from Equations (2) and (6) how A 2 G (t) evolves as a cosine of the sine of the satellite elevation. The frequency of this cosine, 2h λ L1 , is proportional to the antenna height. This means that by estimating the frequency of the observed SNR, we can obtain the height of the receiver.
In order to get an accurate estimate of the frequency of cos(ϕ R (t)) with classic approaches, one must observe at least one period of the signal. For a given initial elevation θ el (t 0 ), we define ∆θ el , the satellite elevation variation required to observe one period of the signal. Based on Equation (2) and using trigonometric identities, we can thus write: = 0 (7) Figure 2 shows the corresponding elevation variation required according to antenna height, for different θ el (t 0 ) values. In particular, we can see that for a height of 3 m, one period of the cosine can be observed for a satellite elevation variation of at least 2 • when θ el (t 0 ) is close to 0 • . If we consider, for example, a mean satellite elevation speed ω el = 10 −3 • /s, we must thus wait at least 33 min to observe one period of the signal. In the next section, we propose a normalization procedure to decrease the required observation period.

Normalization of the GNSS Signal Amplitudes
As described in Equation (2), the phase ϕ R (t) is a function of the satellite elevation and of the antenna height. Since the satellite elevation evolves slowly, we propose a calibration procedure that uses instead a variation of the antenna height. From Equation (6), the minimum and maximum values of A G (t) can be obtained when cos(ϕ R (t)) is equal to −1 or one, respectively. They are defined by the following expressions: From these two equations, we can deduce the sum of the square of the amplitudes and their product as: Therefore, upon substituting Equations (11) and (10) into Equation (6), the single unknown parameter to estimate will be the phase delay ϕ R (t), which is proportional to the height of the antenna and the sine of the known satellite elevation angle.
In order to always get the maximum and minimum value of A G , the variation of ϕ R should be greater than or equal to 2π. According to Equation (2), the minimum variation of the antenna height dh should thus be equal to: In Figure 3, we show the value of dh min as a function of the satellite elevation. From this figure, we can see that a variation in the antenna height of 0.5 m is sufficient to observe a maximum and a minimum value for A G when the satellite elevation is higher than 12 • .

Frequency Estimation
After down-conversion, the received signal is correlated with a local code replica. For the following derivation, we define the samples y[n] as the noisy post-correlation measurements obtained every t = nT int , where T int is the coherent integration time. From Equation (6), we can define y[n] as: where w[n] is the resulting zero-mean additive white Gaussian noise (AWGN) with power σ 2 n . Note that if the pre-correlation noise samples are colored (e.g., due to front-end imperfections or interfering signals), we assume that spectral whitening has been used to optimally process the RF front-end output samples (see, e.g., [16]). From Equation (2), we note that ϕ R [n] evolves linearly as a function of Therefore, the factor β defines the frequency of cos(ϕ R [n]), and its evolution is defined as a function of the sine of the elevation. The satellite elevation θ el [n] is obtained from the current GPS ephemeris data and the estimated position of the GNSS receiver. In order to estimate β, after calibration, we can define the following model for A G [n]: and derive the maximum likelihood estimate of β for N measurements as: Finallyĥ is a function ofβ defined byĥ In the next section, we will derive the CRLB for Equation (14) in order to make a feasibility study and assess the expected performance of the proposed approach.

Cramér-Rao Lower Bound
We are interested in assessing the maximum theoretical accuracy that can be obtained when estimating the receiver height h. Unfortunately, Equation (14) is highly nonlinear, which makes it difficult to directly assess the impact of its different parameters over the estimation error. This nonlinearity is due mainly to the cosine function in the expression and is exacerbated by the presence of the root mean square. Instead, we propose to compute the CRLB for the signal model under consideration. The CRLB provides a lower bound on the variance of any unbiased estimator and, thus, will allow us to assess the performance of our estimator [17].
The signal model considered for A G [n] is provided in Equation (14). In order to provide more insightful results, we will express the reflected signal amplitude as A R = αA D , where α represents the attenuation coefficient due to reflection, assumed to be real and less than or equal to one. In addition, we define γ [n] 4π λ sin (θ el [n]). Thus, we obtain: where ξ = [A D , α, h] T is our unknown deterministic parameter vector and w[n] are zero-mean AWGN samples with variance σ 2 n . The CRLB for ξ can be expressed as [17]: ii (20) where I (θ) is the Fisher information matrix (FIM). A full account of the derivation of the FIM for the considered signal model of y[n; ξ] can be found in Appendix 1. The obtained FIM is: The CRLB for h can be obtained by computing [I −1 (θ)] 33 , resulting in: where is the post-correlator SNR when only the direct signal is received. For simplicity, the g() function is used to represent the multiple terms of The purpose of the following discussion is to identify the effects of {h, α, λ, ∆θ el } parameters through simulation. In order to highlight the effect of a specific parameter, the CRLB is computed for different values of the parameter of interest, while the rest are kept fixed. We have set λ = λ L1 , corresponding to the GPS C/A L1 wavelength, a sampling period T s = 1 sample/s and α = √ 0.7 (water surface scenario for a typical smooth sea [18]) in all of the considered cases, unless otherwise specified. The CRLB is computed from the probability density function (pdf) of the data observations. The estimation accuracy, lower-bounded by the CRLB, is related to the dependency of the data pdf on the unknown parameters (ξ). The CRLB is higher when the dependency between the observations and the parameter to estimate is weak.  Initial elevation:

Assessing the Initial Elevation on the CRLB
Remarks: • When θ el [0] is close to zero, the CRLB is several times higher than for higher θ el [0] values and almost independent of h. In order to explain this behavior, we have shown in Figure 5Left  • When A G [n; ξ] (see Equation (19)) is at its maximum or minimum value, for a small period of observation, the signal tends to a constant equal to A D √ 1 + α 2 + 2α or A D √ 1 + α 2 − 2α, respectively. In this region, the signal evolution can be assumed to be monotone with almost null variation. The dependency between the observations and h decreases with the absolute value of the slope of A G [n; ξ] during the period of observation. In this case, the CRLB is as high as the slope of the signal is low.
• We show in Figure 5Right that the frequency of the cosine evolution of A G [n; θ] decreases when θ el tends to 90 • . For h = 2 m and an observation interval of 600 s, we observe periodic CRLB increases with the satellite elevation when Equation (19) includes a minimum or a maximum. This situation corresponds to the peaks that appear in Figure 4 for From the observed behavior of the CRLB, we prefer using the proposed estimator with satellite elevations between 10 • and 70 • and antenna heights greater than 2 m to obtain better estimation accuracy. For low elevations, the CRLB is indeed high, independent of the antenna height. For elevations close to 90 • , the CRLB increases for low antenna heights, because just a small portion of the signal A G [n; ξ] oscillation, considerably less than a period, is observed. As expected, this effect is exacerbated when the ω el is lower, and a shorter evolution of A G [n; ξ] is observed for the same measurement period.  since it approximately corresponds to a minimum of the A G [n; ξ] evolution, which corresponds to the worst scenario for the estimation of its frequency of oscillation. These ∆θ el values were generated by assuming constant ω el up to 10 × 10 −3 • /s for a fixed measurement period of 600 s and a sampling rate of one sample/s. From the figure, we observe that by increasing ∆θ el , the variance decreases quickly at first, until an entire period of oscillation of the signal model is covered. This relation can be used to set the duration of the measurement time for our estimator, since we are interested in achieving the maximum possible accuracy with the minimum number of data observations. Unfortunately, the period of oscillation of our model depends on the true value of h and ω el . The measurement time required to achieve a certain target accuracy will depend on ω el , and in general, for faster ω el , shorter measurement periods will be required to achieve a similar estimation accuracy.

Estimator Performance Evaluation with Synthetic Data
In order to complement the CRLB analysis, we want to assess the performance of the estimator proposed in Section 3 with synthetic data generated using real measurements of the satellites elevation (θ el ). These data were generated using the signal model defined in Equation (19) with a constant SNR D =18 dB for the direct signal and a sampling period of 1 s. The SNR D value selected is typically achieved at the post-correlation stage and matches the SNR corresponding to a carrier-to-receiver noise density C/N 0 = 45 dB-Hz with a front-end bandwidth of 3 MHz, sampling at the Nyquist frequency, and a coherent integration time of 1 ms. The power ratio between the reflected and direct signals was set to α 2 W = 0.7. The satellite elevation measurements correspond to the satellites in view at Calais, France (50 • 55 ′ 14.093" N, 1 • 56 ′ 59.44" E), on 25 September 2013. The root-mean-square error (RMSE) was computed for 120 observation periods of 10 min, each period starting one minute later than the previous one, with N = 600 observations taken at one sample/s. We assume here that we are in a classic bistatic configuration depicted in Figure 1, where the receiver antenna is located at h = 2 m above the ground. The parameter ϕ R is defined by Equation (2), and the height h is sought with a resolution step of 1 × 10 −3 m in a bounded search space h space = [0, 5] m.
In Figure 8a,b, we show the satellites' elevations (θ el ) and their elevation rates (ω el ), respectively, as a function of time. We present in Figure 8c,d the RMSE ofĥ for the proposed estimator as a function of time, computed using the Monte Carlo method with 100 realizations. Figure 8c shows the RMSE of the satellites reaching an elevation rate close to zero at some instant during the full observation interval. Figure 8d shows the RMSE of the satellites with a high elevation rate, except for the end of the observation interval. We observe in Figure 8c that the RMSE increases when the satellite elevation rate decreases to less than 2 × 10 −3 • /s. In the case of Satellite 6, the RMSE is higher, because the satellite elevation is close to 90 • . We observe in Figure 8d that the RMSE of Satellites 3 and 22 increases for high elevations. This is due to a low elevation rate and, so, a low frequency of the SNR evolution. These results are in accordance with the CRLB study in Section 4.1.
Let us now compute the RMSE ofĥ for different SNR values and observation interval durations. As before, the SNR refers to the signal-to-noise ratio for the direct signal, and the power ratio between the direct and reflected signals is kept at α 2 W = 0.7. In Tables 1-3, we report the RMSE obtained with 1000 realizations of the noisy A G signal, for observation intervals of 600 s, 300 s and 150 s, respectively. The duration of these observation intervals was selected considering that for a satellite elevation of 35 • and an elevation rate ω el = 6.8 × 10 −3 • /s (e.g., Satellite 3 at 12h02 UTC, in Figure 8), a complete period of signal is observed after 600 s, a half a period after 300 s and a quarter of a period after 150 s.
From the results in Tables 1-3, we observe that the proposed estimator is consistent and that the RMSE increases when the SNR decreases, which was expected. In these tables, we have defined two different sets of satellites. The first set includes Satellites 3, 6, 21 and 22 (in bold in the Tables) with a high absolute elevation rate ω el ≥ 6 × 10 −3 • /s. Satellites 3, 21 and 22 have low initial elevations between 20 • and 35 • , and Satellite 6 has a high elevation, superior to 70 • . In the second set, we consider the Satellites 7, 16 and 18, with a low ω el . Satellite 7 has a low initial elevation of 18 • ; Satellite 18 has an initial elevation of 44 • ; and Satellite 16 has a high initial elevation of 75 • .   Tables 1 and 2 approach the CRLB values computed in Section 4.1. In these cases, we reach centimeter accuracy due to the high elevation rate of these satellites. However, for a period of observation of 300 s and a C/N 0 = 35 dB-Hz, we reach just decimeter accuracy, showing the limits of this technique. In Table 3, we see that centimeter accuracy is not reached for an observation interval of 150 s, even with C/N 0 = 50 dB-Hz. For 150 s, the SNR evolution covers only a short part of the oscillation period during the interval of observation and can be considered monotonic. In this context, we observe in Table 3 that the RMSE strongly depends on the considered part of the SNR cosine evolution rather than on the satellite's ω el (e.g., Satellites 3, 6, 21 and 22).
For Satellite 6, the results shown in Tables 1 and 2 are less accurate despite the fact that its ω el is similar to those of Satellites 3, 21 and 22. This can be explained by the higher elevation of Satellite 6. This result is indeed in accordance with the CRLB study (Section 4.1), because in this case, the SNR evolves with a lower frequency, so we do not cover a full period of the SNR variation over the observation interval.
For Satellites 7, 16 and 18, the results shown in Tables 1-3 are also in accordance with the expected accuracy from the CRLB study. In this case, we reach just meter accuracy, because we observe only a small fraction of the SNR variation period due to the the low θ el of the satellites.

Experimental Results
In order to assess the proposed method, we constructed the following experimental setup to measure the height difference between two antennas, as depicted in Figure 9. The height difference is estimated with the interferometric approach described in Section 3. The advantage of the proposed assessment method is that we can have centimeter knowledge of the system geometry. Next, we derive the link between the height difference of the two antennas and the frequency of variation of the GNSS signal power. For a height difference h between the two antennas, we have (see Appendix 2): where a is the path difference between the GNSS signals for the two antennas, ∆x the distance between the two antennas in the ground plane, θ el the elevation for the satellite in view and φ the angle between the two antennas and the ground plane. The paraxial approximation is assumed, so that the satellite signal arrives at both antennas with the same elevation angle, θ el . ∆θ Az is the angle between the vertical plane containing the two antennas and the vertical plane containing the satellite and any of the antennas. Then, it follows that: with: where arctan * () is the quadrant-specific inverse of the tangent. Finally, we can write: The phase delay between the direct signal received by the two antennas, ϕ exp R , can be expressed as: ϕ exp R depends on θ el and on the satellite azimuth, with ∆θ Az assumed to be constant over the entire observation time. We then can conclude that ϕ exp R evolves linearly, with a slope β exp = 2π c f L1 h, as a function of K sin(θ el + K 0 ). Finally, β exp = 1 2 β, where β is the frequency of the variation of the SNR as a function of sin(θ el ) in the reflection scenario described in Section 2. Figure 10 shows the experimental vehicle and the telescopic mast used. The figure also shows the system that gives us the ability to precisely change the height of the antenna installed on the top of the mast. The second direct antenna is situated on the roof of the car at a horizontal distance of ∆x = 1.92 m. The height h between the two antennas is known, so we can derive the value of φ and define the complete geometry of the experimental system. We used the Novatel OEM4-G2 ProPak RT2W (GPS + WAAS/EGNOS) commercial receiver [21]. The two antennas were connected using a passive RF-combiner. The C/N 0 measurements, as well as the satellite elevation values were provided by the receiver at a rate of 1 Hz. Figure 10. The vehicle and its telescopic mast.

Assessment with Real Data
We show in Figure 11 an example of C/N 0 evolution as a function of K sin(θ el + K 0 ). In this figure, we differentiate two periods of time in the signal. These two periods correspond to two different processing steps: the calibration step and the estimation step.
We report in Tables 4 and 5       We can conclude that, after the calibration step, the proposed estimator can achieve centimeter accuracy under the experimental setup for a period of observation of 600 s.

Conclusions
In this article, we used an IPT to estimate the height between an antenna and a ground surface, where a GNSS signal has been reflected. We proposed to normalize the SNR measurements in order to construct a model of its evolution over time. The proposed estimator is based on two steps: A calibration step and an estimation step. The aim of the calibration step is to measure the maximum and minimum values of the SNR (or, equivalently, the C/N 0 ) amplitude, in order to model the SNR variations for a bounded interval of possible surface heights.
The maximum likelihood estimate of the antenna height constructed with this nonlinear model was assessed in a study of the CRLB of the model. In this study, we showed that the accuracy of the estimation can be defined as a function of the satellite elevation, the elevation rate, the C/N 0 and the power ratio between the direct and reflected signal.
In order to assess the method, we used synthetic data and verified that the proposed estimator is consistent with the number of observation samples and the C/N 0 of the GNSS signal. We also showed that, using the proposed calibration step, we can expect centimeter accuracy with half a period of the SNR oscillation when we are in the best scenario. These conditions were identified using the CRLB.
Finally, we proposed an experimental framework that uses two direct signals received in different locations. The results using real data from field measurements showed that the proposed estimator can provide centimeter accuracy for a period of observation of 10 minutes within this framework.
In this way, the FIM is defined as a symmetric 3 × 3 matrix as follows: The path difference a is defined with point H, which corresponds to the orthogonal projection of point A over BD, where BD is the line from the mast antenna to the satellite in view. We thus have Y H = 0, so y H = x H tan θ el and X H = x H cos ∆θ Az .