Derivation of the Cramér-Rao Bound in the GNSS-Reflectometry Context for Static, Ground-Based Receivers in Scenarios with Coherent Reflection

The use of the reflected Global Navigation Satellite Systems’ (GNSS) signals in Earth observation applications, referred to as GNSS reflectometry (GNSS-R), has been already studied for more than two decades. However, the estimation precision that can be achieved by GNSS-R sensors in some particular scenarios is still not fully understood yet. In an effort to partially fill this gap, in this paper, we compute the Cramér–Rao bound (CRB) for the specific case of static ground-based GNSS-R receivers and scenarios where the coherent component of the reflected signal is dominant. We compute the CRB for GNSS signals with different modulations, GPS L1 C/A and GPS L5 I/Q, which use binary phase-shift keying, and Galileo E1 B/C and E5, using the binary offset carrier. The CRB for these signals is evaluated as a function of the receiver bandwidth and different scenario parameters, such as the height of the receiver or the properties of the reflection surface. The CRB computation presented considers observation times of up to several tens of seconds, in which the satellite elevation angle observed changes significantly. Finally, the results obtained show the theoretical benefit of using modern GNSS signals with GNSS-R techniques using long observation times, such as the interference pattern technique.


Introduction
During the past two decades, researchers have been studying the use of the reflected Global Positioning System (GPS) and other Global Navigation Satellite Systems' (GNSS) radio signals for remote sensing purposes [1]. This is commonly referred to as GNSS reflectometry (GNSS-R). The GNSS signals are transmitted at the L-band frequencies, which can penetrate cloud cover and are particularly sensitive to soil moisture, snow water content and sea-ice salinity. The GNSS-R concept was first proposed by Martin-Neira in [2], in 1993, for its application on ocean altimetry from airborne and spaceborne platforms. Since then, the use of GNSS-R for ocean altimetry, wind speed and sea state determination has gained considerable traction. Several works have been published describing specific GNSS-R instruments and architectures [3][4][5], retrieval techniques and algorithms [1,6], theoretical signal models [7][8][9][10], simulators [11][12][13] and several measurement campaigns carried out from different platforms [1,14,15]. Even specific GNSS-R space-based missions have been planned by different space agencies [3,16,17]. Furthermore, more remote sensing applications have been also proposed and validated through experimental campaigns, such as: soil moisture monitoring [7,[18][19][20][21][22][23][24][25], snow others such as [59,60], have focused on GNSS (mostly GPS) ocean scattered signals, which are mostly incoherent, i.e., with significant speckle noise.
In an attempt to fill this existing gap, in this paper, we make use of the CRB to provide the theoretical estimation bound for scenarios dominated by coherent reflection when the LOS and the reflected signals are received by the same antenna, both over short and long observation periods and for different GNSS signals: GPS L1 C/A and L5 and Galileo E1 and E5. In our previous work [57,58], the CRB was computed considering the samples at the output of the receiver's prompt correlator during the tracking of the composite signal (i.e., LOS and the reflected component) as our measurements. In this paper, we assume no specific architecture for the GNSS-R receiver. The CRB is derived for the signal model corresponding to the samples obtained at the output of the receiver front-end. This implies that we are now considering all of the information contained in the received signal when computing the CRB, while in previous works, some information could be lost due to the processing of the received signal. Thus, the CRB expressions derived in this paper represent a more fundamental performance bound, independent of the GNSS receiver architecture considered.
This paper is organized as follows. Section 2 describes the GNSS-R scenario under consideration, summarizing all of the relevant assumptions made and later introducing the corresponding signal model for this scenario. Section 3 contains the derivation of the CRB for the LOS and an arbitrary number of multipath components over short measurement intervals, which is further adapted to our specific scenario using the Fisher information matrix (FIM) transformation approach. In Section 4, we analyze the impact of different parameters, e.g., the propagation path difference, over the CRB obtained. Section 5 extends the derivation of the CRB for long observation times, while in Section 6, we discuss the results obtained when computing this extended CRB as a function of the scenario parameters. Finally, Section 7 summarizes the results and highlights the paper's main conclusions.

Scenario Definition
In this section, we describe the scenario for GNSS-R coherent reflections that we have considered for this work. By doing so, we summarize all of the relevant assumptions that are needed to define the received GNSS signal model at the output of the receiver front-end.
As expected, the properties of the reflected signal when impinging the receiver antenna will strongly depend on the scenario geometry and on the reflecting surface itself [61]. It is precisely that dependency that allows us to retrieve useful information from the signal. In this paper, we study a scenario similar to the one defined in [7].
In such a scenario, the GNSS-R receiver is ground-based and stationary during the data acquisition. This receiver has a single antenna simultaneously receiving the line-of-sight (LOS) signal and a single reflected signal, where the latter has been forwardly scattered from the specular reflection point on the surface. The geometry of our scenario is shown in Figure 1. The reflection surface is assumed to be a horizontal planar surface with a homogeneous layered structure, i.e., the dielectric permittivity of the medium will only depend on the depth, i.e., the z-axis. These conditions could be fulfilled by surfaces, such as bare soil, snow-or vegetation-covered soil, or a calm water surface. Because of the large distance separating the satellite and the receiver, the elevation angle of the incident signal at the specular reflection point and at the receiver antenna can be considered both equal to the angle θ el . In addition, we consider the reflection surface in our scenario to be smooth according to the Rayleigh criterion [61], and thus, the standard deviation of the surface height σ s h , i.e., the surface roughness, satisfies: where λ is the carrier wavelength. By default, we will assume a σ sh = 0.01 m, low enough to fulfill the Inequality (1) for the GNSS L1 and L5 carrier frequencies and θ el ≤ 35 • . The reflected signals are usually characterized as the sum of two contributions or components: one specular and another one diffuse [61]. Under the smooth surface assumption, the coherent component reflected from the specular reflection point will clearly dominate, and as a consequence, the contribution of the diffuse component will be disregarded in our analysis. Under such an assumption, we can consider the energy reflected mainly coming from the first Fresnel zone [7,61], and the reflected signal can be well approximated by a single delayed and complex-attenuated replica of the LOS signal. The first Fresnel zone describes an ellipse on the reflection surface (denoted as S in Figure 1) with: where a represents the semi-major axis, oriented towards the satellite's azimuth; and b is the ellipse's semi-minor axis. As was also pointed out in [7,61], the waves carrying the GNSS signals will acquire some degree of depolarization upon their reflection on dielectric mediums. While the wave corresponding to the LOS signal will be considered as purely right-hand circularly polarized (RHCP), the reflected wave can be described as a mix of both RHCP and left-hand circular polarization (LHCP) components. The receiver's antenna radiation pattern is assumed to be known, and its effect under different polarizations of the incoming signals is included within the modeling of the amplitude of the received signals. Furthermore, we make no assumption about the maximum height for the receiver's antenna phase center position in our scenario, as long as the smooth criterion defined by Inequality (1) is still respected. Without loss of generality, we constrain our analysis to the single-satellite and single-signal scenario in each case. This assumption is reasonable given the low crosstalk among different GNSS signals, as well as the low cross-correlation among the different satellite spreading codes [62,63]. In this paper, we have focused only on the GPS L1 C/A and GPS L5 (both L5I and L5Q) signals, as well as on the Galileo E1 (E1B: data; and E1C: pilot) and the full Galileo E5 signals. Table 1 summarizes the characteristics of the signals considered, and Figure 2 shows their normalized autocorrelation (real part).   Table 1 for more details.
Although the GNSS signals in Table 1 have different typical received powers [64][65][66], for comparison purposes, we assume a default carrier-to-noise ratio (C/N 0 ) of 45 dB-Hz for the LOS signal in all of the cases, before considering the antenna gain.

Signal Model at the Output of the Receiver Front-End
For the single-satellite scenario described, we collect x ∈ C N×1 , a vector representing the complex baseband signal received, sampled at the output of the receiver's front-end over a total observation time interval of T obs seconds. We assume the sampling rate f s , such that T obs = N f −1 s . The n-th sampling instant is defined as t n = t 0 + n f −1 s , where t 0 is the time where the first measurement was taken. The vector x can then be expressed as: where the GNSS signal, described by the vector r, has been expressed as M superimposed independent signal components r m , each one corresponding to a different propagation path. In Equation (3), n is circularly-symmetric complex Gaussian noise with covariance C. Thus, x is a random vector distributed as CN (r, C). Each r m component can be expressed as: where u m and φ m are, respectively, the amplitude and the phase at the instant t 0 (referenced to the receiver's local oscillator) of the m-th individual component. u m and φ m are assumed to remain constant during observation intervals of up to a 1-s duration. These parameters already include the effects of the receiver antenna's gain for each signal component. In Equation (4), the vector has been defined as: with each of its component being equal to: s (t) ∈ C represents the received filtered baseband navigation GNSS signal spread by the pseudorandom code of the satellite considered. In addition, for each component, we have that: τ m is the time-delay at the instant t 0 ; α m = (1 + v m /c) is the coefficient capturing the effect of the Doppler shift, where c is the signal propagation speed and v m the equivalent relative velocity between the transmitter and the receiver; and ω c = 2π f c , where f c is the signal carrier frequency. τ m and α m will also be assumed constant over time intervals of up to 1 s. This assumption for α m is justified by the very low rate of change of the Doppler shift observed by a static receiver [67]. Given the maximum rate of change of the relative speed possible in this case, i.e., dv m /dt| Max ≈ 0.178 m/s 2 [67], α m can vary as little as ≈5.94×10 −10 within a 1-s interval. In this paper, the front-end's low-pass filter is implemented as a fifth-order Butterworth filter. The distortion introduced by this filter is already taken into account in the signal model s(t). The sampling frequency of the front-end is considered equal to the Nyquist rate with complex sampling, i.e., f s = 2B f e , where B f e is the one-sided bandwidth (3 dB) of the low-pass filter. Finally, no quantization loss is considered in our analysis.

Derivation of the CRB for Short Observation Intervals
In this section, we derive the CRB given the output samples of the receiver's front-end, whose signal model was defined in the previous section, over a total observation time T obs of up to 1 s. After introducing the CRB, we derive the general solution for M propagation paths. Then, we particularize the solution obtained for our scenario, described in Section 2, i.e., with M = 2, corresponding to the LOS signal plus a single coherent specular reflection.

Definition and General Case for M Propagation Paths
Estimation of the error bounds plays a major role in parameter estimation in general, since they provide benchmarks to assess the performance of practical estimators. The workhorse in this respect is the well-known Cramér-Rao bound. The CRB defines a lower bound on the variance of any unbiased estimator for a set of deterministic parameters when the probability density function (pdf) of the measurements is known ( [68], Chapter 3). The CRB quantifies the dependency of the measurements' pdf on the unknown parameters to estimate. The stronger the pdf's dependency on these parameters, the lower the CRB. The CRB for a vector of L unknown real-valued parameters ξ ∈ R L×1 states that the covariance matrix of the estimates, Cξ, must satisfy (under standard regularity conditions): where J (ξ) is the Fisher information matrix (FIM), whose inverse is the CRB matrix. The elements in the FIM are defined as: where ln p(x; ξ) is the log-likelihood function of ξ, given the vector of random measurements x.
Taking the inverse from Equation (8), we can express the CRB of the i-th parameter in ξ as: In Section 2, we have introduced the vector x, modeling the discrete output of the receiver front-end during the time interval T obs . As described there previously, the vector x is assumed to obey a circularly-symmetric complex Gaussian distribution. Its expected value, defined as µ, will then be r = ∑ r m , as defined in Equation (3). We can express each signal component r m as a function of the following unknown parameters, grouped into the vector ξ m as: where, for the m-th signal component, u m represents its amplitude, φ m its phase shift, τ m its time-delay or code phase and α m its Doppler shift coefficient. In this case, we want to jointly estimate the vector ξ ∈ R 4M×1 , which has been defined as the vector concatenating all of the ξ m vectors. Now, given the Gaussian distribution of the vector x assumed in our case, it is possible to make use of the Slepian-Bangs formula ( [68], Chapter 3), which provides a convenient way to compute Equation (8) as: where {·} denotes the real part of a complex vector or scalar, tr{·} denotes the trace operator, the superscript H is used to denote the conjugate transpose operator and the C represents the covariance matrix of x. Given that the noise component is independent from any of the parameters in ξ, then the first term of the sum in Equation (11) will be equal to zero, due to ∂C/∂ξ i = 0. In order to compute the remaining term in Equation (11), we can express the elements of the gradient of µ as: where the vector d m is defined as in Equation (5). It is worth noticing that J(ξ) can be conveniently expressed as the following symmetric block matrix: where every element of each sub-matrix J ξ m ξ l is computed as: with ξ i ∈ {ξ m } and ξ j ∈ {ξ l }. Equation (17) can be used to compute the elements of J(ξ), and the CRB of ξ is straightforwardly obtained as in Equation (9). This case corresponds to the CRB under the assumption of Gaussian noise, independent of the received GNSS signal.
Let us now examine the case where only thermal noise is present. Since we have considered the sampling rate to be a multiple of the Nyquist sampling frequency, the samples in x can be considered uncorrelated, which implies a diagonal covariance matrix, i.e., C = σ 2 w I. The noise variance σ w can be computed as: where N 0 /2 is the assumed white noise spectral density in W/Hz, H Rx ( f ) is the front-end's normalized baseband filter response in the frequency domain and B n is defined as the equivalent noise bandwidth (one-sided). Thus, for complex white Gaussian noise (CWGN), we can express Equation (17) as: If the narrowband approximation for the received signal is taken, i.e., α {m,l} ≈ 1, then the CRB obtained as a result of inverting the matrix defined in Equation (16) is analogous to the results described in [69] (Chapter 4), where the CRB was derived for two propagation paths using a similar approach. In addition, Skournetou et al. used in [70] an approach similar to the one described in this paper, but for the case of up to four signal propagation paths.

Two Propagation Path Case for a Ground-Based Static Receiver
For the scenario described in Section 2, we assume M = 2, with the r 0 component used to represent the line-of-sight (LOS) signal and r 1 used for the specular reflection. In this scenario, we consider a stationary receiver, where subsequently, the LOS and the specular reflection will experience the same Doppler shift; thus, α m = α l := α. This simplifies the FIM J(ξ) defined in Equation (16). Its elements related to α [m,l] and a different parameter in ξ, as defined in Equation (19), will now become zero. This implies that the estimates of α will be independent from the estimates of the rest of parameters, i.e., they will decouple. The same conclusion was reached in [69] (Chapter 4) when computing the CRB for two propagation paths. Having decoupled parameters allows us to disregard the estimation of the α for the rest of our analysis. The reason is that we are only interested in estimating the geometry of the scenario and the characteristics of the reflecting surface. In a stationary environment, α cannot be related with any useful information about the scenario that we want to retrieve. Thus, we are now redefining the vector of unknown parameters ξ, by redefining each ξ m to just [u m , φ m , τ m ] T .
For the scenario described in Section 2, we have only two propagation paths; using Equation (16), the FIM can be expressed as the 6 × 6 block matrix: where each submatrix, computed using Expression (19), can be expressed as: with m, l ∈ {0, 1} and ∆φ φ l − φ m . Let us now introduce the discrete autocorrelation function (ACF) of s(t), defined as: The superscript * is used to denote the complex conjugate. Now, we introduce the following identities: Notice that in order to keep our notation as simple as possible, we use R (n) (∆τ) to indicate the n-th order partial derivative of R with respect to τ evaluated at the delay difference ∆τ τ l − τ m . The Identities (23)- (25) are obtained straightforwardly by applying the chain rule, and they allow us to express Equation (21) as: Now, since we can express J ξ m ξ l as in Equation (26), therefore J(ξ) will be a function of the ACF of s(t) and first and second order derivatives with τ evaluated at ∆τ. It is helpful to notice that ACF and its derivatives can be expressed using their relationship with the signal's power spectral density (PSD). By making use of the well-known Wiener-Khinchin theorem [68] (pp. 576-577) and the properties of the Fourier transform, we have that: where S s ( f ) represents the signal s(αt)'s PSD. Notice that s(αt) represents the filtered signal at the output of the receiver's front-end, scaled in the time domain due to the Doppler effect. It already takes into account any distortion effect on the signal due to the front-end's filtering. Thus, in our derivation, these effects are taken into account within the R() and its derivatives. Nonetheless, if we are interested in evaluating the influence of the front-end's filtering on the resulting CRB, it might be more convenient to express S s ( f ) as: where S( f ) is the PSD of the signal received at the antenna. In the context of GNSS-R, in [54], the authors followed this approach to study the effects of using different front-end bandwidths and filter types in the estimation performance of the interferometric GNSS-R technique. Moreover, they considered also the band limitations of the signal transmitted in practice. In contrast, in our approach, we generate a model of s(t) simulating the samples at the output of the front-end, and then, we compute its ACF and its derivatives. The selected front-end bandwidth, B f e , will define the sharpness of the normalized autocorrelation peaks, which will have a great impact on its derivatives, R (n) . In addition, the total number of samples considered within T obs will also depend on B f e , as N = T obs /(2B f e ). If we consider the support of R(∆τ) to be approximately limited to the interval [−T c , T c ], where T c is the chip duration, then, according to Equation (26), we have J ξ m ξ l ≈ 0 for propagation path differences greater than ρ chip = T c c, where c is the signal propagation speed. This is a consequence of the autocorrelation properties of PRNcodes [62]. Moreover, J (ξ) becomes a diagonal matrix, which implies that all of the parameters in ξ are now uncoupled. From the FIM calculation point of view, this is equivalent to the scenario of having a receiver with two synchronized RF front-end channels, with uncorrelated noise components, and two different antennas: an up-looking one, exclusively receiving the LOS signal, and a second down-looking antenna, exclusively receiving the reflected component.

Introducing Phase Coherence Using FIM's Parameter Transformation
So far, we have derived an expression for the FIM of ξ given the scenario defined in Section 2. However, in the GNSS-R context, we are actually more interested in a different set of unknown parameters that can be related to the geometry of the scenario, as well as to the properties of the reflecting surface. Moreover, up to this point, we have not taken into account in our signal model the phase coherence presented by specular reflection. This implies that the LOS and the reflected signals maintain a fixed phase relationship with each other and with respect to the local oscillator clock, used as the reference, i.e., φ 1 = g(φ 0 ).
In our approach, we have decided to compute the CRB for ξ and then transform it for a new set of parameters Ψ. This transformation of J(ξ) allows us to conveniently introduce the phase coherence assumption. We define the new vector of unknown parameters as: where u 0 , φ 0 , ρ 0 are the amplitude, phase and propagation distance (ρ = τ c) for the LOS signal component. Γ = |Γ|e jφ Γ is the reflected signal coefficient, i.e., the complex attenuation experienced by the reflected signal when compared to the LOS. Its module is defined as: Γ takes into account the surface reflectivity coefficient, the antenna complex gain and the polarization mismatch. Thus, it will depend on the satellite elevation angle (θ el ) and on the azimuth angle of arrival of the reflected signal to the antenna. Since θ el varies very slowly [67], we will assume it to remain constant for observation periods up to 1 s. h is the height of the receiver, defined as the orthogonal distance between the antenna phase center and the reflection surface. In our flat surface scenario, this height h relates to the propagation path difference between the LOS and the reflected component ∆ρ as: where ∆ρ is defined as ∆ρ = (τ 1 − τ 0 ) c. In addition, because of the coherence assumption on the reflected component, the phase difference between the two components, i.e., ∆φ = φ 1 − φ 0 , can be expressed as: where λ is the carrier wavelength.
Following the approach described in [71], we can express the FIM for Ψ as: where G is the Jacobian transformation matrix defined as: The expressions for each of the elements of G are provided in Appendix A. The CRB for the parameters in Ψ can be directly obtained as the diagonal elements of J −1 (Ψ). This two-step approach, i.e., first computing the FIM for ξ and then transforming it, has been used for its simplicity. It allows us to redefine the bound for a different set of parameters and to introduce more assumptions into our signal model without having to derive the FIM expression from scratch.

Analysis of the CRB for Short Observation Intervals
The matrix J −1 (Ψ) derived in the previous section shows that CRB(Ψ) in the specular reflection scenario considered will depend on: i The scenario geometry: through the propagation path difference between the two signal components ∆ρ, which will be a function of h and the satellite elevation angle θ el (and its variation) during the observation time. ii The properties of the reflecting surface: through its reflection coefficient, Γ, that can be modeled as a function of the signal's incident angle and the electrical properties of the surface and its roughness. iii The GNSS signal considered and the receiver features: the signal modulation and its transmission bandwidth, the C/N 0 of the LOS signal, the antenna's radiation pattern and the receiver front-end's bandwidth (and the filtering scheme considered) will all have an impact on the resulting CRB.
However, understanding these dependencies only by examining J −1 (Ψ) is not straightforward. The expressions for their diagonal elements can be, in the general case, quite cumbersome with several terms. Nonetheless, we examine the specific case when ∆ρ > ρ chip , where the initial FIM J(ξ) is a diagonal matrix, which greatly simplifies the terms in the final J −1 (Ψ).

CRB when ∆ρ > ρ chip
In this particular case, not only the diagonal terms for the resulting J −1 (Ψ) matrix, i.e., [J −1 (Ψ)] ii , are much simpler, but this assumption is particularly relevant because it represents the best estimation case, i.e., when its CRB lower bounds all other cases with shorter ∆ρ. As shown in Appendix A, for h and |Γ| in Ψ, after some simple algebra, we obtain the following CRB expressions: where the coefficient γ introduced has been defined as: and SNR 0 represents the equivalent SNR as if only the LOS signal was present and, as such, can be expressed as: Notice that the amplitude u 0 takes into account the effects of the receiver's antenna gain pattern for the received LOS signal, which is dependent on the signal incident angle, i.e., a function of θ el . From examining Equations (35) and (36), we can directly deduce the following: 1.
With little surprise, both CRB expressions are inversely proportional to SNR 0 . However, CRB(h) is also inversely proportional to sin 2 (θ el ). This poses a trade-off: we require low elevations for the smooth surface assumption to hold, but at the same time, the lower the elevation, the higher the CRB(h). Moreover, as shown by Equation (2), low elevation angles imply larger first Fresnel zone areas, which decrease the spatial resolution of our estimates.

3.
In both CRB expressions, the effects of the assumed received signal bandwidth and the front-end filtering (possible losses and distortions) are modeled within the R (0) and the ACF's peak, i.e., R(0). R (0) can be understood as the curvature or the sharpness of the R(0). The higher the front-end bandwidth, the sharper the resulting ACF peak and, thus, the higher R (0). This holds as long as the front-end bandwidth is narrower than the received signal bandwidth, otherwise there will be no sharpening on the resulting R(0). Additionally, for a higher front-end bandwidth, the SNR 0 will also increase if the Nyquist sampling rate assumption is maintained, as a consequence of having more samples for the same T obs .
In Figure 3, we show the results obtained when computing the CRB for h using Equation (35), as a function of the front-end's bandwidth B f e for the different GNSS signals considered. Note that in Figure 3, the y-axis represents (the convention of representing the √ CRB in the figures is followed consistently along the entire paper) σ Ψ i = CRB(Ψ i ). An isotropic antenna pattern was assumed, and |Γ| was set to √ 0.1, which corresponds to an attenuation of 10 dB in the reflected signal power with respect to the LOS. The T obs considered was equal to 1 s. The ∆ρ was fixed to 600 m to ensure isolation between the LOS and the reflected signal components provided by the code correlation properties. We observe how the CRB of h consistently decreases with the increase of B f e for all GNSS signals studied, as long as B f e remains lower than the (one-sided) signal bandwidth B Tx . The slight decrease of the CRB of h for all GNSS signals observed in Figure 3 for B f e > B Tx /2 is explained due to the use of a non-ideal low-pass filter, like the one explained in Section 2, to model the limited transmission bandwidth. Also from Figure 3, we see how the CRB(h) stays in the meter order, except for the Galileo E5 signal, which reaches decimeter order precision when B f e > 40 MHz, due to its highest bandwidth.  The CRB for the rest of parameters in Ψ is jointly obtained when computing J −1 (Ψ). In Equation (32), we have expressed the phase difference between the two signal components as the sum of two contributions: φ Γ and (2π/λ)∆ρ, with the latter capturing the shift due to the extra propagation path corresponding to the specular reflection. For short observation intervals, where both of these contributions to ∆φ remain constant, it is not possible to individually estimate them precisely. This issue is easily spotted by examining the CRB (φ Γ ) expression, i.e., the fifth element in [J −1 (Ψ)] ii , derived in Appendix A, that can be approximated as: where k = 2π/λ is the wavenumber for the received signal carrier. Since the numerator in Equation (39) will usually be many orders of magnitude higher than its denominator, resulting in CRB (φ Γ ) 2π rad, this implies that we cannot estimate φ Γ properly. Given this ambiguity in individually estimating the two contributions to ∆φ, performing precise phase altimetry is only possible if φ Γ is assumed to be known or small enough that it can be neglected. Otherwise, a bias in the altimetry estimates is introduced. Thus, even if we considered the phase coherence in our derivation of CRB(h), Expression (35) is equivalent to the CRB for pure GNSS-R code altimetry only, i.e., as if no phase information were used. The CRB for code altimetry can be also derived following the same proposed approach, described in the previous section, but just assuming no phase coherence between the LOS and the reflected signal, i.e., considering ∆φ as the difference between two independent phases when computing the G transformation matrix defined in Equation (34), which leads to the same results as in Equation (35). In addition, in [72,73], in an effort made by the authors to characterize the code-tracking accuracy of GPS receivers, they computed the minimum code thermal noise jitter when using noncoherent delay-lock loop (DLL) discriminators, i.e., σ 2 DLL . The expressions provided in their results are very frequently used to characterize the tracking performance of different receiver architectures [62] (Chapter 5). Additionally, we can use these same expressions to compute the minimum variance on the h estimate, i.e., σ 2 h , in a straightforward way, as: where σ 2 DLL−0 and σ 2 DLL−1 are the minimum variances of the LOS and the reflected signal code delays, respectively. In Equation (40), two independent tracking channels are assumed: one tracking the LOS signal and the other tracking the reflected signal. The results when computing Equation (40) match the CRB(h) values obtained using Equation (35) when ∆ρ > ρ chip , which correspond to the flat part of the curves later shown in Figure 4. As a consequence, it follows that a two-tracking channel architecture for code altimetry is efficient in terms of the bound defined in Equation (35) given a sufficiently large propagation path difference. In practice, such a scenario might appear to be unlikely for ground-based receivers. Nonetheless, many of the new GNSS signals have higher bandwidths, which implies lower ρ chip values (e.g., 293 m for GPS L1 C/A signal vs. 29.3 m for GPS L5) [62].   Table 1). The two vertical lines represent the chip lengths: 293.25 m, of the GPS L1 CA and Galileo E1 BC signals; and 29.32 m, of the GPS L5 and Galileo E5 signals.
In the previous discussion, we have highlighted the problem arising when trying to separately estimate the two contributions to ∆φ and how the CRB for the joint estimation of Ψ captures this effect. However, what if we assume some prior knowledge of one of the two terms in ∆φ? We examine now the impact on the CRB if we assume that one of these terms is known, which represents the best case scenario given such an assumption. We study two cases: (1) when φ Γ is assumed to be known, referred to as "phase altimetry" case; and (2) when h is known, referred to as "reflection coefficient estimation" case.

Phase Altimetry
In this case, we assume the angle φ Γ as known. Now, in order to obtain the new expression CRB alt (h), we remove the row and the column of J(Ψ) associated with φ Γ and then compute the matrix inversion J −1 (Ψ), as is described in Appendix A. We obtain the following expression: by making the approximation of R (0)/c 2 very small (e.g., ≈0.1 for the worst case, i.e., the Galileo E5 signal with B f e = B Tx ) compared to k 2 R(0) (e.g., ≈623), given the signal bandwidths considered. The phase information is now effectively used and the bound multiple orders of magnitude lower. As expected, it is directly proportional to λ 2 , i.e., the lower the carrier wavelength, the higher the possible precision of the h estimate. In addition, we observe that when comparing Equation (41) with Equation (35), the dependence on R (0) can be neglected now. This implies that in this phase altimetry case, the estimation performance of the phase altimetry techniques should be negligibly impacted by the bandwidth of the signal used. As a consequence, the use of the new GNSS with higher bandwidths, but with the same carrier frequency, should bring little improvement to the achievable precision of the phase altimetry techniques for short T obs when ∆ρ > ρ chip .

Reflection Coefficient Estimation
Analogously, we compute now CRB r for |Γ| and φ Γ (subindex r), but in this case, we proceed by removing the row and the column associated with h, as shown, once again, in Appendix A. We obtain the following expressions: As expected, the CRB r (|Γ|) is the same as in the joint estimation of Ψ, while CRB r (φ Γ ) is now orders of magnitude smaller when compared to Equation (39) and well bellow 2π, thus implying that in this case, it is possible to obtain useful phase estimates.

CRB when ∆ρ < ρ chip : Interference Case
When ∆ρ < ρ chip , we cannot process anymore the two signal separately, which was otherwise possible thanks to the correlation properties of the PRN codes. Now, the initial FIM J(ξ) is no longer a diagonal matrix, and the analytic expressions for [J −1 (Ψ)] ii become quite cumbersome. Thus, from now on, we will rely on numerical evaluation of [J −1 (Ψ)] ii to compute the CRB.
In Figure 4, we show the results of computing CRB(h) as a function of the ∆ρ for the different GNSS signals under study. Once again, we have considered a |Γ| 2 = 0.1, and T obs = 1 s. As expected, we have observed that when ∆ρ increases, the bound reaches the values obtained using Equation (35). Moreover, we also observe how that the CRB reaches these values even for ∆ρ ρ chip . In summary, we can identify two consequences of having higher bandwidth signals: (a) the asymptotic value for the bound is reached faster; and (b) the value of the bound is lower for any ∆ρ. In the case of (a), this can be clearly related with having a sharpened peak of the ACF for higher front-end bandwidths.
In Figure 5, we show the results of computing CRB(h) for different ∆ρ and |Γ| values. The plots presented allow us to visually determine, approximately, the order of the best achievable precision given a ∆ρ and |Γ| combination. The y-axis represents |Γ| 2 , which in practice is equivalent to the gain experienced by the reflected signal compared to the LOS. As expected, for ∆ρ values higher than ρ chip , the CRB agrees with the one obtained with Equation (35). Figure 4 confirms that the CRB when ∆ρ > ρ chip lower bounds all of the other cases, as we hypothesized at the beginning of this section.

Derivation of the CRB for Long Observation Intervals: CRB(Ψ long )
So far, we have assumed total observation times of up to one second. However, the IPT uses much longer observation times, usually on the order of several minutes. Can we use the CRB expressions obtained in Sections 3.3 and 4 for such long observation times? Unfortunately, we cannot use them for observation times T obs longer than one second. This is due to the fact that, for T obs > 1 s, the approximation described in Section 2, of considering the parameters in Ψ to be constant over T obs , is no longer valid. In this case, the change in the delay of the received LOS signal cannot be well approximated anymore as a linear expression with an initial value of τ 0 s. In addition, neither the Doppler shift observed, nor the phase shift φ 0 with respect to the receiver oscillator can be considered as constant. That is due to the non-linear relative movement of the satellite and the receiver over long observation times and due to the signal propagating through the atmosphere, which effects the signal as it varies with time [62]. In addition, the direction in which the signals impinge on the antenna will vary over time as a consequence of the variation of the satellite's elevation, i.e., θ el , and the azimuth angles observed. Subsequently, this variation can result in a change in the LOS signal amplitude u 0 and an extra phase shift, i.e., a variation of φ 0 , since these two parameters also model the effect the receiver's antenna radiation pattern.
In realistic scenarios, modeling the reflection coefficient Γ, even for homogeneous and horizontal flat surfaces, can be a complex matter. Nonetheless, we still can characterize Γ as a function of the vector s Γ , containing the set of parameters related to the reflection surface physical properties; the signal's impinging angle on the surface; and the antenna gain towards the specular reflection point [61]. The two last terms can be expressed as a function of θ el . Thus, we have that: If we consider h and the parameters in s Γ constant over T obs 1 s, then we can compute the CRB for a new vector of unknown parameters Ψ long , defined as: where u 0 , φ 0 and τ 0 are column vectors of L elements, which we refer to as "fast-varying" parameters. In our analysis, these parameters are associated with the values of u 0 , φ 0 and τ 0 at the l-th time interval of T coh =1 s duration, with L = T obs /T coh being the total number of these intervals for T obs . We assume that θ el remains constant over each T coh interval. The previous assumption is based on the fact that the rate of change of θ el over time observed by ground-based static receivers is very low, i.e. on the order of 10 −3 • /s [67].
We can now use the same FIM parameter transformation approach to compute the CRB(Ψ long ), which was introduced in Section 3.3. This approach allows us to obtain J(Ψ long ) as a function of the J(Ψ) for each of the different L short time intervals T coh . In order to do so, we first define the vector Ψ ext of unknown parameters as: where each subvector Ψ l corresponds to the l-th T coh observation interval. We build now the matrix J (Ψ ext ) as: where each J(Ψ l ) is the FIM computed using Equations (20)-(34) for the l-th T coh observation interval.
The J (Ψ ext ) generated is equivalent to assuming all of the parameters in any Ψ l are independent from any of the parameters in a different subvector Ψ k . This assumption is only used as a means to introduce the parameter dependencies between observation intervals by transforming J (Ψ ext ), once again, using a new Jacobian transformation matrix A ∈ R 6L×P , defined as: where i and j are the indexes used for the rows and columns, respectively. The number of columns in the matrix A, P, is equal to the number of parameters in Ψ long , i.e., 3L + 1+ the number of elements in s Γ . Since the parameter h was already in Ψ and is considered constant for all of the L intervals, ∂h l /∂h = 1, for any l. Finally, J(Ψ long ) is obtained as: A more detailed description of how the matrix A is constructed is provided in Appendix A. The CRB can be obtained directly as [J −1 (Ψ long )] ii . Furthermore, we shall point out that the previous derivation is independent of the model of Γ considered, as long as the gradients of |Γ| and φ Γ with respect to s Γ exist. No assumption was made regarding the number of parameters in s Γ . Nonetheless, it is well known that the CRB can only increase with the number of unknown parameters considered, which intuitively makes sense, given the increased uncertainty.

Analysis of the CRB(Ψ long )
In this section, we compute the CRB, using Equation (50), for different simulation cases. Our goal is to study the effect of the received signal's satellite elevation span covered during the observation time (∆θ el ), as well as the impact of the GNSS signal modulation, over the theoretical achievable precision in the estimation of the parameters in Ψ long .
The variation of θ el during the observation time has the following implications: 1. The propagation path difference ∆ρ will change due to the displacement of the specular reflection point. In addition, the angle of arrival of both the LOS and the reflected component to the antenna will also vary.

2.
If the reflection coefficient Γ is assumed to depend on θ el , this will also change during the observation time.
Up to this point, we have not selected any specific model for the reflection coefficient Γ. From now on, we will use the model for a single-layer bare soil introduced in [7], which we will call the Zavorotny-Larson model (Z-L). This model is described in more detail in Appendix B, and it allows us to specify Γ as a function of θ el and the following parameters: where ε r and σ are the real part of the surface's relative permittivity and the surface conductivity, respectively, which are assumed constant over the entire observation interval. Notice that this model for Γ is used as an example in our analysis and that the approach proposed to obtain J(Ψ long ) can accommodate other models for Γ(θ el ; s Γ ).
In addition, the gain of the receiving antenna will have an important impact over the CRB. To assess this impact, we consider two possible antenna cases in our analysis: 1.
An ideal case with an isotropic antenna model.

2.
A Leica GNSS AR25 antenna [74] model. Only the antenna gain's amplitude has been modeled. The antenna phase center is assumed to be located at the distance h perpendicular to the reflection surface. The antenna boresight is pointing upwards, parallel to the surface normal vector (i.e., an elevation of 90 • ). We have used the pattern model provided as part of the open source GPS multipath simulator [75].
In both cases, we make the assumption that the antenna complex gain is accurately known and that all of the available signal's energy is considered by defining the receiver's front-end bandwidth as B f e = B Tx /2.

CRB(Ψ long ) vs. ∆θ el
In this first case, we compute the CRB for different h values and elevation spans ∆θ el , covering from θ el = 10 • -35 • . In the results obtained, the first individual ∆θ el value represents a special case, since it corresponds to a single short observation interval of T obs = 1 s, with no change in θ el . In this case, the bound matches the results obtained using the short-period CRB derived in Section 3. Figure 6 shows the CRB of h results for the isotropic antenna and the AR25 antenna cases, with the GPS L1 C/A signal. In the figure, the highest CRB values are located at the left-upper corner of the plot, which correspond to the lower values of antenna heights (below 5 m) and elevation spans (below 15 • ). The CRB values decrease almost monotonically with both the increase of h and the elevation span covered. As expected, it becomes clear from the plots that the antenna pattern has a very significant impact on the CRB of h given the same h and ∆θ el . Figure 6a shows that sub-millimeter precision could be achieved in theory in the ideal isotropic antenna case. On the other hand, in the AR25 antenna case, the plot in Figure 6b shows that now, only centimeter precision can be achieved. These CRB values now obtained are significantly lower than the CRB values computed previously in Section 3. Firstly, given the satellite elevation rate of 5 × 10 −4 • /s considered, the observation time is above the hour duration for ∆θ el > 18 • . The longer the observation time considered, the more independent samples we collect, thus the lower CRB values obtained. Secondly, the CRB computed now takes into account that the phase difference information can be used for the estimation of h in this case. Intuitively, this can be explained by the different rate of change of the two contributions to the phase of the reflected component: the phase shift due to the scenario geometry (related to h) and the phase shift due to the polarization change due to the reflecting surface properties, i.e., to Γ.  Having a higher true h has two main implications. First, the higher the h, the faster the oscillation rate of the interference pattern observed in the composite amplitude of the received signal with the change of satellite elevation. Second, the ∆ρ will be higher, which makes it easier to differentiate between the two signal components, given the code autocorrelation properties, and the resulting CRB can be lower, as we have shown in Section 3. However, techniques like the IPT only use the output of the estimated SNR or the prompt correlator (single output), and they can only work if there is interference between the LOS and the reflected signal, i.e., if ∆ρ < ρ chip [21]. Thus, with IPT, we face a trade-off: we require low ∆ρ values in order to observe the interference pattern, but according to the CRB expressions obtained, low ∆ρ values degrade the achievable estimation precision.

Effects of the Signal Modulation on the CRB(Ψ long )
We now present the results of computing the CRB of s Γ and h, given ∆θ el spans of different lengths and the GNSS signals considered. For each signal, we have considered a particular B f e = B Tx /2. We have considered the Leica AR25 antenna model and a true receiver height h = 3 m. We also have used the Z-L model for Γ(θ el ) given a dry ground surface with a permittivity (real part) ε r = 30 and a conductivity σ = 0.2 S/m [76]. The initial θ el is 10 • for all of the elevation spans considered. The results obtained are summarized in Figure 7. As expected, the longer the ∆θ el covered during the observation time, the lower the CRB values obtained, for all three parameters. However, this improvement is not linear. We observe how the CRB decreases significantly for ∆θ el covering up to ∼2 • , while it decreases at a much slower rate when extending ∆θ el beyond that value. This observed behavior can be explained as follows. For our scenario and given the GNSS signals' carrier wavelength, i.e., from 0.19-0.25 m, a ∆θ el of ∼2 • is required to observe a 2π phase shift in the reflected signal phase, caused exclusively due to the change in ∆ρ. In the IPT context, this corresponds to observing a single complete oscillation of the SNR interference pattern. Longer ∆θ el , i.e., longer observation times still provide some improvement. It is important to remark about the existing trade-off between extending ∆θ el , i.e., extending the observation time, and the displacement of the specular reflection point over the reflection surface. The change in θ el also changes the size of the first Fresnel zone. For low elevation angles, such as the ones that we need to consider to ensure the smooth surface assumption, the change in the area covered by the first Fresnel zone and the displacement of the specular reflection point is significant, as shown in Figure 8. In practice, it is unlikely to find real surfaces that are that flat and homogeneous (and with no obstacles), except maybe for the case of calm water surfaces. On the other hand, the antennas used in practice (e.g., in geodetic sites) have lower cross-coupling at lower elevation angles, i.e., less attenuation on the LHCP signal for these elevations. This is why higher oscillation amplitudes of the SNR interference pattern are observed for θ el angles in standard GNSS geodetic receivers [7,55]. We shall also notice that even for long observation times, where the signal's phase information can be used in the parameter estimation, the greater the signal's bandwidth, the lower the CRB for all of the parameters. This is due to the fact that for a h = 3 m, ∆ρ is small enough for the two signal components to interfere, and the separability provided by the PRN code's correlation properties acts only partially. This separability is needed to measure the code delay difference between the two components, which is required to unambiguously estimate h. This cannot be done by only using the phase difference information.
We would like to highlight that, if during the ∆θ el observed, the Γ coefficient changes significantly, e.g., when the Brewster angle is covered by ∆θ el [21,61], lower values of the CRB for ε r and/or σ can be achieved with shorter ∆θ el . From the CRB point of view, if the pdf of the observed signal changes significantly with the variation of the unknown parameters, such as ε r and/or σ, then it is easier to estimate these parameters from the observed samples, resulting in lower CRB values. Likewise, when the Γ coefficient barely changes during the observation time considered, it becomes harder to unambiguously retrieve the ε r and/or σ parameters from the observed samples, thus resulting in higher CRB values. Figure 7 shows that centimeter precision can be theoretically achieved when estimating h over elevation spans larger than 20 • for all GNSS signals. The CRB also shows that millimeter precision can be theoretically achieved for the ∆θ el = [10 • , 25 • ] in the case of the Galileo E5 signal, i.e., using its full bandwidth. In practice, however, it is unlikely to find real surfaces that are that flat and homogeneous (and with no obstacles), maybe except for very calm water surfaces. Thus, the millimeter precision is unlikely to be achieved by any estimator in a practical scenario, even with the assumptions made in this paper (e.g., perfect knowledge of the antenna radiation pattern, the satellite elevation angle, etc.). The asymptotic behavior observed over long observation times of the CRB of the conductivity σ can be related to the small dependency and change of Γ during the considered elevation spans for the dry ground surface.

Conclusions and Future Work
In this work, we have derived a theoretical precision bound for the GNSS-R coherent reflection scenario, for ground-based static receivers. The CRB has been derived first for observation times of up to 1 s (short observation intervals), and later on, we have extended this bound for long observation intervals of up to hundreds of seconds, where the satellite elevation angle was changing during the observation time. The CRB was computed for GPS and Galileo signals with different bandwidths and different modulations, i.e., BPSK, CBOC (6,1,1/11) and AltBOC (15,10). For short observation times and ∆ρ > ρ chip , we obtained simple analytic expressions of the CRB for the parameters of interest. Later on, the CRB was extended for long observation times by using the FIM transformation approach. The approach followed can easily accommodate different models for the surface's reflection coefficient Γ. In this paper, we selected the Z-L model to describe Γ, and then, we have computed the CRB as a function of satellite elevation spans of different sizes covered during the observation time. The impact of the antenna radiation pattern was considered within the modeling of Γ. The specific case that we have studied (receiver height, antenna pattern and surface properties) is an example selected to understand the dependency of the CRB on the signal properties and the scenario parameters.
For short observation times of up to 1 s, the results obtained have shown that even for antennas designed to mitigate multipath and these antennas being installed with their boresight perpendicular to the ground, it is still theoretically possible to obtain meter precision (or even decimeter in the case of using the full Galileo E5) given a front-end with a large enough bandwidth. In the case of long observation intervals, the CRB computed confirms the link between the variation of the satellite elevation and the achievable precision. Moreover, the CRB evaluation presented confirms the cm level precision that has been reported by many works, e.g., [29,32,38,41,[77][78][79], among others, some of them using standard GNSS geodetic stations [80]. Reasonable accuracy was obtained in the estimation (joint) of ε r and σ over observation times covering satellite elevation variations above 10 • . According to the CRB results computed, the modernized GNSS signals considered in this work, i.e., GPS L5, Galileo E1 and Galileo E5, can potentially provide better estimation precision for all of the parameters when compared to the GPS L1.
In this work, the CRB results presented have been mainly focused on the h estimation. Using the approach described in this paper, further study of the CRB for the parameters related to the reflection surface's properties is planned, in order to characterize the precision achievable in other GNSS-R applications, such as soil moisture determination. In addition, two main future activities are planned to continue with the work described in this paper: first, to further explore the CRB for similar scenarios with different models for the Γ and, second, to compute the CRB for the scenarios where experimental campaigns were carried out, as a way to fully validate the bound obtained and assess the performance of the techniques used.
As a final remark, we want to point out that a practical limitation of the precision bounds provided in this paper is the fact that our signal model (for coherent reflection) might not accurately match the true reflected signal. In many practical scenarios, the reflected signal will be better described as a combination of the coherent and incoherent components. The importance of each component will depend on the considered scenario. In any case, the CRB results provided in this paper will lower bound those obtained with a signal model considering also the presence of incoherent reflections. We are considering as future work the derivation of CRB expressions for different kinds of ground-based receiver scenarios, which is not straightforward from the analytic point of view.
Author Contributions: Miguel Angel Ribot and Cyril Botteron conceived of the idea of deriving CRB within this context. Miguel Angel Ribot carried out the theoretical derivations, designed and implemented the MATLAB software, performed the simulations and wrote the paper. Cyril Botteron contributed the analysis and the discussion of the results. All of the authors contributed to the revision of the manuscript.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Expressions for the Transformation Matrices Used and Analytic Expressions for the CRB in Short Observation Times
Here, we provide more detail about how the Jacobian transformation matrices, G and A, which were described in Sections 3.3 and 5, are computed. In addition, we show how the CRB expressions for short observation times discussed in Section 4 were obtained.
As discussed in Section 4.1, the FIM J(ξ), which is expressed using Equations (20) and (26), becomes a diagonal matrix when ∆ρ > ρ chip , i.e.: However, we are interested in computing the CRB for the set unknown of parameters Ψ, already defined in Equation (29) as: by transforming J(ξ) into J(Ψ) (which is not diagonal anymore), as described by Equation (33), i.e.: As discussed in Section 3.3, the Jacobian transformation matrix G for short observation intervals and this specific case with phase coherence between both signal components is expressed as: With some algebra, the CRB for Ψ is then obtained as: CRB(h), CRB(|Γ|), CRB(φ Γ ), defined in Expressions (35), (36) and (39), are obtained in a straightforward manner from the 3rd, 4th and 5th elements of Expression (A5), respectively. The CRB expression presented in Equation (41) for the phase altimetry case, described in Section 4.1.1, is obtained by removing the fifth row and the fifth column (both related to Φ Γ ) from J(Ψ) and then computing the inverse of this matrix. Now, the CRB alt (h) is the element related to h in the main diagonal, i.e., the fifth element. In the same way, the CRB for the reflection coefficient estimation case, i.e., CRB r , described in Section 4.1.2, with Equations (42) and (43), is obtained following the same procedure, but now removing the last row and the last column (both related to h) from J(Ψ).
In Section 5, the Jacobian transformation matrix A ∈ R 6L×P was introduced, with L being the number of T coh intervals considered and P the number of parameters in Ψ long . For a s Γ = [ε r , σ] T , as in Section 6, we have the following: where the vector δ k ∈ R 1×6 has been defined as a row-vector with all elements equal to zero, except for the position indicated by the subindex k ∈ {0, . . . , 5}, which has a value of one. |Γ(θ el [l])| and φ Γ (θ el [l]) will change together with θ el in each of the L intervals of T coh duration, over the total observation time.

Appendix B. Zavorotny-Larson Model for the Surface Reflection Coefficient Γ
In [7], the authors proposed a model for the complex reflection coefficient Γ = |Γ| exp {φ Γ } of the GPS L1 signal for the case of bare soil with a plane surface. Γ is defined like in Expression (30) as the amplitude ratio between the reflected and the LOS signals. The model treats the soil as a single-layer dielectric medium with losses. It also takes into account the polarization dependent antenna gain effect on the reflected signal. The main geometric optic contribution to the amplitude of the reflected signal comes from the first Fresnel zone around the specular reflection point. According to the model, when the antenna boresight is perpendicular to the ground, Γ can be expressed as a function of the satellite elevation θ el as: where the effect of the surface roughness is incorporated by the exponential term, the so-called coherence loss factor [61]. k = 2π/λ is the carrier wavenumber, and σ s h is the standard deviation of the surface height. F R (θ el ) and F L (θ el ) represent the antenna complex amplitude gains for right-hand and left-hand circular polarizations (RHCP and LHCP) at θ el , respectively; and are the polarization-dependent reflection coefficients, where the first subindex stands for the polarization state (R for RHCP and L for LHCP) of the incident wave, while the second subindex stands for the polarization state of the reflected wave. V v and V h are the coefficient for the linear vertical and horizontal polarizations, which can be expressed as: with: the complex dielectric permittivity of the material [81,82]. In Equation (B4), ε r represents the real-part of the permittivity, and σ is the conductivity of the reflection surface (in S/m units). The model for Γ can be extended for the multiple layer ground case, where Γ will also depend on the thickness of each layer and their relative permittivity [81,83]. For example in [84,85], the authors use a three-layer (air + two soil layers) model for GNSS-R applications working with the GPS L1 band signal. For more details on all of the assumptions made to obtain (B1), we refer the interested reader to [7].