A Novel Analytical Formulation of the Correlation of GNSS-R Signals Scattered by a Natural Fractal Surface

Natural surfaces exhibit scale-invariance properties and power-law spectra over a wide range of scales. For this reason, they are well described in the frame of fractal geometry and, in particular, using fractional Brownian motion (fBm) 2-D statistical processes. In this letter, we present the derivation of the correlation of Global Navigation Satellite System reflectometry (GNSS-R) signals scattered by a fBm surface. We here show that this correlation depends on a parameter related to the root-mean-square (rms) surface slope as measured at the electromagnetic wavelength scale. When this parameter increases, the correlation smoothly decreases from a value close to unity reaching the value of the roughness-independent expression already available in the literature. In our experiments, based on roughness measurements available in the literature, we first show how the description of the roughness of natural surfaces can be conveniently obtained via fBm parameters. Then, we illustrate the behavior of the correlation time for different roughness regimes. For a high-altitude airborne receiver, values ranging from about 6 ms (as predicted by the expression already available in the literature) to several tens of milliseconds are obtained.


I. INTRODUCTION
I N THIS letter, we focus on the evaluation of the corre- lation of Global Navigation Satellite System reflectometry (GNSS-R) signals scattered by natural rough surfaces.This topic is of great interest not only in GNSS-R [1] but also in monostatic and bistatic synthetic aperture radar (SAR) interferometric applications [2].Therefore, this problem has been long since tackled for very rough surfaces, for which a classical roughness-independent expression of the correlation is available [1], [2].For the ocean surface, a dependence of signal correlation on wind speed is observed only in case of low-altitude low-velocity GNSS-R receivers [3], [4], [5].
Stimulated by the interest in novel GNSS-R land applications and by the finding that a dominant near specular incoherent component can be appreciated even in case of very flat gently undulating land surfaces [6], [7], [8], this topic has been recently reconsidered, especially for near-specular scattering (which is of interest for GNSS-R applications).In this case, a roughness-dependent expression of the correlation is expected, as first confirmed through simulations in [7] and [9] and then by the analytical formulation of the correlation obtained in [10] by the authors of this letter.In particular, the obtained expression of the correlation coefficient decreases for increasing roughness from a value close to unity (for gently undulating surfaces) to the value predicted by the classical roughness-independent expression [1], [2] (for very rough surfaces).
In [10], as well as in most of the relevant literature, the roughness is assumed to be well described through a stationary 2-D random process, so that the process can be characterized using two statistical parameters, namely, the root-mean-square (rms) height σ and the correlation length L, upon choosing an appropriate model for the autocorrelation function (usually, Gaussian).However, when measuring the roughness of natural surfaces, one finds out that the obtained values of σ and L strongly depend on the size of the measured profile and, in particular, their values increase with increasing profile lengths [8], [11], [12].This makes their use in the modeling of scattering from natural surfaces problematic [11].It is widely recognized that natural surfaces exhibit power-law spectra and scale-invariance properties over a wide range of scales and that this behavior should be described through fractal geometry, e.g., using fractional Brownian motion (fBm) surfaces [12].Hence, the availability of an expression of the scattered field correlation coefficient in terms of fBm parameters is of great practical relevance.
In this letter, we extend the approach of [10] to the case of fBm surfaces.Preliminary results have been discussed in [13].In the experimental section, based on recent roughness measurement data available in [8], we first show that flat (at the topographic scale) natural surfaces are well modeled by fBm surfaces at the scale of interest for the considered sensors and applications.Then, the behavior of the correlation coefficient and of the correlation time as a function of surface roughness is investigated.
The fBm process is statistically nonstationary, with stationary increments, and has infinite variance [12].However, actual natural surfaces exhibit the fBm behavior only up to an outer scale l that may be their linear size or, for a sea surface, the dominant wavelength.These are known as physical (or bandlimited) fBm random surfaces and they are statistically stationary [12], [14], [15], with finite variance where C( ) is the surface normalized autocorrelation function.It is clear from (1) that C( ) is not Gaussian and that its second derivative C ′′ ( ) diverges in the origin.Since the results of [10] can be applied to random rough surfaces satisfying the conditions discussed in [10, Appendix B] (among them, regularity of the derivatives of the correlation function), in our case of interest, they cannot be used.Therefore, it is necessary to devise a new, original analytical formulation of the correlation coefficient.Its rationale is based on proper modifications of the approach of [10], and it is detailed in Section II-B.

B. Evaluation of the Correlation Coefficient
In the following, we refer to the geometry illustrated in Fig. 1 (3-D view) and Fig. 2 (projections on the x −z and x − y planes), where the two receiver positions are those assumed by a single moving GNSS-R receiver at two slightly different acquisition times.At time t, the transmitter is placed at T ≡ (x T ,0, z T ), with x T = −r T sinϑ 0 and z T = r T cosϑ 0 ; and the receiver is placed at R 1 ≡ (x R1 ,0, z R ), with x R1 = r R1 sinϑ 0 and z R = r R1 cosϑ 0 , so that the origin O is the specular point at time t for the mean plane and the plane x-z is the vertical plane containing transmitter and receiver at time t.Both transmitter and receiver move with velocities v T and v, respectively; but since for all spaceborne and airborne GNSS-R systems r R1 ≪ r T , so that v T /r T ≪ v/r R1 , we can neglect the movement of the transmitter, assuming that at time t + t, the transmitter is still in T , whereas the receiver has moved to the point R 2 ≡ (x R2 , y R2 , z R ), with x R2 = r R2 sinϑ R2 cosφ R2 = x R1 + v x t and y R2 = r R2 sinϑ R2 sinφ R2 = v y t, where v x and v y are the xand y-components of the receiver velocity, so that v = (v 2 x + v 2 y ) 1/2 .We further assume that |v x | t≪ r R1 , |v y | t≪ r R1 , so that, see Fig. ( Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
We consider the resolution cell that includes the specular point, so that the observed area is centered in the origin and its size is related to the sensor spatial resolution [10].Finally with and with the subscript X that must be replaced by T , R1, or R2 as needed.Therefore, while r X are the distances of sensors from the origin, R X are their distances from the generic point (x, y,0) of the mean plane, and R X are their distances from the generic point (x, y, z(x, y)) of the rough surface.
The correlation coefficient of the fields is defined as where E(R 1 ) and E(R 2 ) are the generic components of the electric fields at R 1 and R 2 , evaluated according to the Kirchhoff approximation (KA).By following the same procedure as in [10], we get, see [10, eq. ( 10)]: where k = 2π/λ is the wavenumber: is the sensor illumination function, peaked around the origin, so that A x and A y are the x and y sensor resolutions that are assumed much smaller than r R1 ; F(x, y) is a slowly varying function whose expression, as shown in [7], is not of interest for evaluating the correlation coefficient, and, see [10, eq. ( 12)] where u z1,2 (x, y) = (z T /R T (x, y)) + (z R /R R1,2 (x, y)).
Here, at variance with [10], we use (1) in ( 8) and assume ≫ 1 (i.e., σ larger than the electromagnetic wavelength λ , that is about 20 cm for GNSS-R), thus obtaining Proceeding again as in [10], we get cov where and J 0 is the zero-order Bessel function.The integral over in (10) also appears in [14], where it is shown that it can be expressed via a series expansion around u x y = 0 (which is the value of u x y at the specular point, i.e., in the origin).By arresting the expansion at the second order, we get where is the gamma function.Similar to w 2 (x, y), the exponential function in (11) is peaked around the origin and its width, which must be compared with the sensor resolution, can be evaluated by expanding the exponent around the origin: by using the same approach employed in [10, Appendix C] to obtain [10, eqs. ( 22) and ( 23)] we get exp x y (x, y) where It can be noted that k 1−H s ∼ sλ H λ , which is the rms slope at the wavelength scale.
From this point on, the results of [10] apply, in which ( 13) replaces [10, eq. ( 23)].Accordingly, we find with W x,y = A x,y G x,y /(A 2 x,y + G 2 x,y ) 1/2 .The correlation time τ is conventionally defined as the value of t such that ρ has decreased to the value 1/e and, in this case, it is equal to For A x,y ≪G x,y , we get W x,y ∼ =Ax,y and which is the classical solution of [1].This leads to a correlation time τ Conversely, for G x,y ≪ A x,y , we have W x,y ∼ = G x,y and, using ( 13) in ( 14) and Equations ( 18) and ( 19) reduce to the results of [10] for H →1, with s assuming the meaning of rms slope.

III. EXPERIMENTAL RESULTS
As a first step, we here show that the fBm model is well-suited for the description of natural surfaces, even in the case of gently undulating terrains.To this aim, we consider the results obtained in [8] from airborne Lidar data acquired in May 2020 over the CYGNSS cal/val sites in San Luis Valley, Colorado.Two areas are analyzed there: the Z1 site, which is a very flat agricultural area, and the Z4 site, which is a flat area containing some terrain undulations.The estimation of the rms height using Lidar patches of different linear sizes, ranging from about 1 m to about 140 m, was performed in [8], thus producing plots of the rms height σ as a function of the linear size l for the two sites.They are reported in [8,Fig. 8].We extrapolated these data plotting them on a logarithmic plane (see Fig. 3).These graphs show an almost perfect fit with a linear behavior, as predicted by modeling the surfaces with a physical fBm process.In fact, for an fBm, as reported in Section II-A i.e., moving to logarithmic quantities

TABLE I HIGH-ALTITUDE AIRBORNE GNSS-R SYSTEM PARAMETERS
Therefore, we can conclude that fBm is a good model for the considered sites.In addition, according to (21), it is possible to obtain an estimate of the fBm parameters s and H by applying a simple linear regression to the data shown in Fig. 3.The obtained estimates for the two sites are H = 0.93 and s = 0.06 m 0.07 for Z4 and H = 0.83 and s = 0.01 m 0.17 for Z1.
Let us now move to analyze the behavior of the correlation coefficient and of the correlation time as a function of surface roughness.We consider the high-altitude airborne GNSS-R system, whose parameters are reported in Table I [10].Note that the sensor is moving along the x-axis, so that v y = 0.In Fig. 4, we show the graphs of the correlation coefficient (setting t = 7 ms) and of the correlation time as a function of the parameter s of the observed surface, obtained using (14) and (15) with the two values of H estimated for sites Z4 and Z1.It can be noted that in both cases, as s increases, the correlation coefficient decreases from a value close to unity to the value obtained via the classical, roughness-independent expression (16), i.e., about 0.336.Also, the correlation time decreases from several milliseconds to the value obtained by (17).This value is achieved for s larger than about 0.04 m 1−H , and for the considered system, it is equal to 6.7 ms.Accordingly, it is expected that for the Z4 site (s = 0.06 m 0.07 ), a value of the correlation time very close to the "classical" one should be obtained.This is confirmed by the fact that using (13), we obtain G x = 1975 m and G y = 1711 m, which are larger than resolution, so that ( 17) can be applied.In fact, the value obtained by ( 17) is in rather good agreement with that obtained for the Z4 site using the exact expression (15), i.e., 7.0 ms.
Conversely, for the Z1 site (s = 0.01 m 0.17 ), we expect to obtain a value of the correlation time significantly larger than 6.7 ms.Indeed, the value of correlation time obtained using the exact expression (15) is 15.2 ms.Note that for the Z1 site, the use of ( 13) provides G x = 230 m and G y = 199 m, which are smaller than resolution, so that in this case, (19) can be applied, providing a correlation time of 14.8 ms, close to the exact value.
These results confirm that for very flat, gently undulating land surfaces, the value of the correlation time might be significantly higher than the one foreseen by the classical roughness-independent expression available in the literature, at least for airborne GNSS-R systems.
IV. CONCLUSION A novel analytical formulation for the expression of the correlation coefficient and correlation time of GNSS-R signals scattered by natural land surfaces modeled as physical fractal fBm surfaces has been obtained.The obtained expressions depend on the parameter k 1−H s, which is the rms slope at wavelength scale.For sufficiently low values of the parameter, ρ and τ depart from the classical roughness-independent values reported in the literature, paving the way to their use in support of coherent GNSS-R processing or roughness retrieval.When this parameter decreases, the correlation coefficient smoothly increases from the value obtained by the roughness-independent expression to a value close to unity.
As a meaningful side result, our analysis also provides a first preliminary evidence of the appropriateness of a fractal description of natural land surfaces also in flat, gently undulating areas, by exploiting roughness measurements independently reported in the scientific literature.This is a new result, since the use of fBm for natural surfaces was previously demonstrated for sea surfaces [16], rocky terrains [17], areas with significant topography [18], and lava tubes [19].

Fig. 2 .
Fig. 2. Geometry of the problem.(a) Projection onto the x-z plane.(b) Projection onto the x-y plane.

Fig. 3 .
Fig. 3. Logarithm of the height standard deviation σ [m] versus logarithm of the profile length l [m] for CYGNSS cal/val site (a) Z1 and (b) Z4.Solid line is the obtained regression line, whereas the dots are the measured values obtained from [8].