GNSS Ocean Bistatic Statistical Scattering in the Time-Varying Regime: Simulation and Model Validation

This work develops around the problem of modeling statistics and correlation properties of the Global Navigation Satellite System (GNSS) L-band signal scattered by the moving ocean surface, presented in the companion article by Principe et al. (2021). The modeling work is here completed by deriving an effective method for simulation of the received signal as well as by making a validation study based on simulated and real data. The simulation of the stochastic process embedding the design statistical properties, the temporal correlation model and the system and surface parameters is intricate, in general, and provides a useful tool for researchers interested in modeling and simulation of GNSS reflected signals from a time-evolving sea surface. The GNSS received waveforms after 1 ms coherent integration are simulated and the correlation properties of the specular and near-specular points are compared to the theoretical model and to results from real data. Real raw data were collected during an aircraft mission in the Gulf of Finland and in a spaceborne scenario with Cyclone Global Navigation Satellite System (CYGNSS) satellites. Comparisons show good agreement among simulations, model, and real data, and demonstrate that the estimated correlation time of the GNSS-R signal fits real data with higher accuracy with respect to the case in which sea correlation is not accounted, especially when the receiving platform velocity is moderate.


I. INTRODUCTION
T HE design and assessment of airborne and spaceborne Global Navigation Satellite System GNSS reflectometers and the development of retrieval algorithms requires testing stages based on simulated observables under the possible configurations of the satellites' geometry and of the sensed environment. Since the studies in [1], investigations have been carried out on how to incorporate in the scattering model all possible relevant sea surface features and their interaction with the wind fields. In this context, the transition from strong to weak diffuse scattering is also relevant when the wind speed is lower than 4 m/s, as reported in [2] and [3].
Direct electromagnetic simulation of the scattering from a discretized realistic sea surface with the appropriate facet Manuscript  size provides accurate results but the computational time required by this approach is not compatible with simulation of statistically significant ensembles of delay-Doppler maps (DDMs) [4]. Other simulation techniques have been developed to overcome the computational complexity of the full electromagnetic approach. In the Cyclone Global Navigation Satellite System (CYGNSS) end-to-end simulator, the surface is sampled with a grid size that is much larger than the wavelength and the bistatic radar cross section (BRCS) is calculated by integration over the interested surface area [5].
The final representation of the received signal is obtained after multiplication by a complex exponential function with uniformly distributed random phase and 2-D convolution with the ambiguity function (AF) of the GNSS signal. The Doppler variation due to transmitter and receiver motion is introduced for each DDM resolution cell. A similar approach is considered in the PAU/PARIS end-to-end performance simulator (P 2 EPS) [6] but the surface grid is defined in the delay-Doppler domain. Using inversion equations and consistency conditions the points are mapped into the coordinates defined by the specular point and the incidence plane.
The simulation framework developed in [7] performs an elegant and computationally efficient way to overcome the problem of coordinates transformation and to incorporate the Doppler changes due to satellites orbital motion into the simulation model. The simulator in [8] is based on the statistical characterization of the received signal, developed in [9]. Here, indeed, the scattered signal is modeled as the random sum of contributions from specular facets inside the delay-Doppler cell. By merging together the statistical and electromagnetic approach, the model provides a way to define the overall contribution from the ensemble of scattering facets without calculation of the single contributions at a microscopic level. The DDM simulation approach is similar to that in [6] and [10] with an improved set of equations for coordinates transformation and computation of the surface area of the resolution cell.
The common drawback of all these approaches is that they do not incorporate the sea surface dynamics into the model. However, it has been demonstrated in [11] and [12] that this component cannot be neglected when the receiver is placed on a low speed or fixed platform or when the platform motion is compensated in some way.
In the context of GNSS land scattering, the modeling and simulation of temporally correlated GNSS signals has been This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ recently introduced, where the microscopic roughness as well as the large-scale roughness have a large impact on assessing the coherency of the reflected signal [13], [14], [15].
Starting from the work in [11] this article has two main purposes: a) to present a description of a simulator that includes the sea surface dynamics and b) to validate the statistical scattering model with real airborne and spaceborne data. The sea surface motion is implemented by introducing a timevarying non-Gaussian process. The model is derived in [11] by approximating the scattered field as the sum of discrete elementary contributions by specular facets in the geometric optic limit of Kirchhoff approximation. The corresponding statistical model is a random walk in two dimensions with stochastic number of steps that, asymptotically, converges to a compound Gaussian process having the form Z (t) = (S(t)) 1/2 G(t) with S(t) a Gamma-distributed process and G(t) a complex Gaussian process. The variance of the Gaussian component depends on the curvature of the scattering facets and on the mean number of scatterers in a cell; the shape parameter of the Gamma process depends on the temporal fluctuations of the number of facets in a cell [9], [11]. The proposed simulation approach is based on linear filtering with autoregressive (AR) systems and is able to exactly reproduce the 2-D probability density function (pdf) and the correlation properties of the compound Gaussian process.
The article is organized as follows. In Section II the details for the simulation of time-correlated scattering are presented, using theorems and properties summarized in Appendix A. The problem of generating a compound Gaussian sequence with exponential wind-dependent autocorrelation function (ACF) is considered and the simulation procedure described in [8] is modified to include platforms motion in the generation of the instantaneous DDM. Section III shows simulated results and their agreement with the model defined in [11]. The validation through comparisons with real airborne data from a 2011 European Space Agency (ESA) campaign and with real spaceborne CYGNSS data is presented in Section IV. Conclusion and future works are shown in Section V.

II. SIMULATION OF THE CORRELATED SEA SCATTERING
The received signal, after conventional delay-Doppler processing, at the time t = nT , can be represented as [11] where T is the coherent processing interval and the sum is over all delay-Doppler resolution cells in the glistening zone. In the expression (1), all quantities with subscript m are evaluated at the center of the resolution cell, τ m = τ − τ m , f m = f − f m , k is the GNSS signal wavenumber, D m is the footprint function that depends on the transmitting and receiving antenna radiation pattern projected onto ground, R m is the reflection coefficient, R 0m and R m are the distances from transmitter and receiver to the mth resolution cell, respectively. The AF of the transmitted waveform is denoted as χ (·, ·) and S(·) and G(·) are the texture and speckle components of a compound Gaussian process, modeling the sea surface scattering. As discussed in [11], the variables Z m = (S m ) 1/2 G m can be reasonably assumed independent of cell to cell and the temporal ACF of the process has been evaluated as where P(·) is the pdf of the sea surface slopes, q m = [q ⊥m , q z m ] is the scattering vector, and A m is the resolution cell area. Therefore, the ACF of the process is composed by the product of two factors: one is a weighted sum of Dopplershifted terms, due to the platforms' velocity, and the other ρ G (nT )e −μ|n|T , is due to sea surface motion and describes the correlation contribution given by the speckle and the texture components of the compound Gaussian process, respectively. As shown in [11] and [12], it is reasonable to assume that the correlation term ρ G (nT ) is a negative exponential function; therefore, it can be written as ρ G (nT )e −μ|n|T = ρ |n| c in the following. The analysis of the ACF (2) is reported in [11] and shows that the decorrelation time is dominated by the platform velocity, but some effect of the sea surface motion is appreciable in the airborne case.
The main difficulty in simulating a temporal sequence of DDMs, as defined in (1), is the generation of the compound-Gaussian process with assigned pdf and ACF. The problem cannot be solved in general and ad hoc techniques can be found in [16] with reference to spherically invariant random process simulation. For an effective design of the simulation scheme, the normalization is helpful, where the sampling period T is omitted for shortness and σ c is the standard deviation of the complex compound Gaussian process. The joint pdf of the variables Z (n − 1) and Z (n) is conditionally Gaussian, given S(n − 1) and S(n) [1,Appendix B], and the covariance matrix of the Gaussian Indeed, the simulation problem is reduced to generate the random process S(n) with a second-order characterization given by the joint Gamma pdf reported in [1, eq. (24)] and to use such variables in the simulation of the Gaussian vector with covariance matrix (4). The details for generation of a correlated compound-Gaussian process are presented in Sections II-A and II-B.

A. Generation of the Sequence of Correlated Gamma Variables
In this subsection, we consider the problem of simulation of a correlated sequence S(n) with first-order Gamma pdf and second-order pdf given through the generating function in [1, eq. (23)]. The procedure that we develop uses some known properties of Gamma random variables, that for convenience are reported in Appendix A.
Therefore, a sequence of correlated Gamma random variables is generated as follows.
1) A vector of L independent Gaussian random seque- . . , L, with exponentialtype temporal correlation and 1-lag correlation coefficient ρ, is generated through AR filtering (Appendix A, Property 3). The procedure yields correlated Gaussian sequences with ACF 2) The sum is a correlated process with first-order Gamma pdf, S(n) ∼ (L/2, 2σ 2 ) (see Appendix A, Property 4). The scheme for generation is sketched in Fig. 1. By choosing σ 2 = 1/L we obtain a unit-mean first-order Gamma distribution. The joint MGF of S(n) and S(n − 1) has the form derived in Appendix B that matches [1, eq. (23)] for L = 2α and ρ 2 = θ(T ) = exp(λ − μ)T where λ and μ are the birth and death rates and α is the shape parameter of the Gamma pdf, as defined in [1, eq. (17) and (22)].

B. Generation of the Compound Gaussian Sequence
The compound Gaussian sequence with exponential correlation is generated according to (3). The entries of matrix (4) show that all covariance terms between the real and imaginary parts of the complex Gaussian variables G 0 (n) and G 0 (n − 1) are zero, whereas the real part and the imaginary part have the same covariance depending on ρ c and on the stochastic ratio between the correlated Gamma random variables S(n−1) and S(n). In our setup, the desired exponential correlations are implemented through AR filtering of Gaussian variables. The simulation is carried out as follows.
1) A sequence of Gamma distributed random variables S(n) is generated according to the algorithm defined in Section II-A. 2) Two sequences of complex random variables are generated, with W 1 (n) and is produced, with the ACF The generation scheme is sketched in Fig. 2.

C. Parameters Setting and Discussion
The temporal evolution and statistics of the simulated process are determined by the design parameters α, the shape parameter of the Gamma distribution, ρ 2 , the correlation coefficient of S(n), and ρ c , the correlation coefficient of the compound process. Guidelines for choosing such values are given in the following.
1) Shape Parameter of the Gamma Distribution: This parameter rules the decay of the upper tail of the Gamma pdf and, therefore, the non-Gaussianity of the scattered signal. Small values of α denote a marked non-Gaussian statistic of the process, whereas for large values of α the Gamma distribution approaches a Dirac delta function centered on the mean value and the process becomes Gaussian. As explained in [9], the shape parameter does not play a role in the 1 ms intensities (DDMs), but it modifies the variance of the DDMs. In Section II-A, the value of α determines the two simulation parameters L = 2α and σ 2 = 1/2α.
2) Correlation Coefficient of the Gamma Process: This parameter rules the temporal decay of the exponential correlation of the Gamma process. In radar backscattering models, this is usually referred to as the correlation of the sea surface texture component and exhibits much slower decay than the correlation of the Gaussian component. In the GNSS reflectometry forward-scattering model, the correlation of the Gamma process is related to the lifetime of the specular facets whereas the correlation of the Gaussian component is governed by the temporal amplitude and phase coherence of the specular facets. Theoretical and experimental determination of specular facet lifetime is complex and beyond the scope of this work. The argument was considered in the articles [17] and [18] in the context of radar scattering from the sea surface. On the other hand, a precise characterization of the correlation coefficient ρ 2 is not crucial because of the product form of the overall correlation of the compound process that will be discussed in the next paragraph. 3

) Correlation Function of the Compound Process:
The correlation function of the compound Gaussian process ρ |n| c = ρ G (nT )e −μ|n|T includes the statistics of the correlated Gamma process and the statistics of the Gaussian process. As discussed in [11], an empirical law for the correlation function of the compound process was considered, following the study proposed in [12]. The correlation model is exponential with parameter 3λ c /U , where λ c is the carrier wavelength and U the wind speed. The validity of the approximation depends on the ratio between the spatial resolution and wind speed. For case of airborne GNSS-R resolution (which is footprint limited), the approximation is valid starting from 4 m/s wind speed and remains applicable for very high winds.

D. Generation of the Instantaneous Correlated DDM
Single 1-ms DDMs are generated following the same procedure as in [9]. The delay-Doppler domain is uniformly discretized into a matrix whose values correspond to the cells centers (τ i , f j ), with i = 1, . . . , I and j = 1, . . . , J , and for each cell the corresponding spatial contribution is generated, having in mind that, for the ambiguity problem, the signal relative to the point (τ i , f j ) is in many cases the sum of signals coming from two different spatial cells. The term in (1) that depends only on geometry is generated as The sum at the right-hand side of (1) can be rewritten in form of a convolution between the scattering function and the AF [8], [12]. Then, the correlated compound Gaussian process Z (t) is generated for each resolution cell with index m as previously explained, and the time-varying scattering function The temporal variation of the scattering function depends solely on sea surface motion through the compound process Z n (nT, τ i , f j ). Finally, the instantaneous DDM can be computed as the 2-D discrete-time convolution between the scattering function and the product of the AF with the complex exponential exp(i 2π f j nT ), that is, where τ i = τ − τ i and f j = f − f j . One can clearly see here that the simulator accounts for dynamic range variations through the Doppler term exp(i 2π f j t) as well as the sea surface motion, through the correlation properties of the terms Z h (nT, τ i , f j ), with a negligible increment of the computational load.

III. NUMERICAL SIMULATION
The main purpose of this section is to present significant examples of simulated zero-Doppler waveforms and to estimate the correlation function in the specular point bin and in other delay bins. Results are compared with theoretical correlation functions obtained under various conditions of platforms setup and wind speed. Simulation parameters have been chosen using geometry and system setup matching two airborne tracks (POPI0038 denoted as Airborne 1 and POPI0018 denoted as Airborne 2) and one CYGNSS satellite track (Florida overpass on day 190, year 2019, space vector (SV) 04, pseudo random noise (PRN) code 02), the same used for analysis on real data which will be discussed in Section IV.
The simulation analysis presented in this section is focused on the airborne experiments. As regards the satellite configuration, the overall signal decorrelation time is dominated by the term e i2π f m t which is of the order of 1 ms [19] and an experiment is shown only in comparison with real data in Section IV. Values of the main parameters used in simulations are reported in Table I where the Earth-Centered Earth-Fixed (ECEF) coordinates of satellite positions and velocities correspond to the platform's initial positions. The ACF is estimated using a set of temporally correlated zero-Doppler waveforms obtained after T = 1 ms coherent integration time. Indeed, the ACF is computed by averaging 30 realizations of the estimator with Y (n) = Y (nT, τ, f = 0), and is compared to the theoretical expression in (2) computed over the same geometrical scenario with N = 1000. The following analysis is carried out considering the nominal specular point bin and two other delay bins, one located on the leading edge and the other on the trailing edge of the delay waveform. Fig. 3 shows a set of 30 realizations of the normalized ACF magnitude, their average, and the model from (2). The simulations are performed for a wind speed U = 9 m/s, in Airborne 1 configuration, from: bin 20, located on the leading edge of the delay waveform; bin 33, corresponding to the nominal specular point bin; and bin 40, located on the trailing edge of the zero-Doppler waveform. Plots in panel 3(a) clearly show the dominant impact of the additive white Gaussian noise over the signal. As reported in [1, eq. (21)], the noise ACF has a triangular shape with peak value N 0 T , which is clearly emerging in the ACF for all delay bins lying more than 1 chip before the specular point. Plots in panels 3(b) and (c) show that the ACF of the signal is now dominant with respect to noise. The correlation time of the process is investigated as a global indicator of the signal time coherence. It is computed through a nonlinear least squares Gaussian fitting, as in [11], [19] and [20] and results are reported in Table II. This confirms a good agreement between model and simulations, resulting in very similar values of the correlation timet cor computed on the simulated ACF and correlation time t cor computed on the model ACF.
A plot of the correlation time versus the delay bin in chips is reported in Fig. 4, for wind speed U = 9 m/s. For delay values less than −1 chip, the correlation time is ruled by the thermal noise and, therefore is less than 1 ms. The maximum value of correlation time is achieved before the specular point and remains nearly constant until a breakpoint value located in the specular point, after which a fast decreas- ing trend is evident. These results are also in agreement with the analysis reported in [20]. The sensitivity of the correlation time to wind speed is analyzed in Fig. 5. Plots show that, as expected, the correlation time decreases with wind speed. Results are also in agreement with [20] and [7]. Looking at the trailing edge of the plots in the figure, it can be noted that the slope decreases with the wind speed. This is possibly caused by the impact of the wind speed that extends over a wider range of delay bins.
An additional analysis has been carried out to investigate the sensitivity of the correlation time to receiver velocity. Results are reported in Table III for several values of the wind speed. As expected, a reduction of the correlation time for higher values of wind speed and platform velocity is observed. Looking into (2), it is found that the ACF depends on the wind speed (through the sea surface slopes pdf as well as on  the statistical correlation of the compound Gaussian process) and on temporal variations of platforms ranges R m and R 0m . Therefore, the decorrelation time has two components: 1) a statistical component, due to sea surface motion, and 2) a geometric component. The statistical decorrelation is weaker, as confirmed by results in Table III, and for platform velocity greater than 100 m/s the decorrelation time is almost entirely due to range variations.

IV. EXPERIMENTAL ANALYSIS
The analysis consists in comparing the ACFs and correlation times for real data, simulated data, and for the theoretical model defined in (2). Simulations are designed with the aim of reproducing geometry and wind parameters of the real acquisitions. The satellite acquisitions consist of a 1000 ms segment of the CYGNSS RAW-IF track acquired on July 09, 2019, processed to obtain complex delay waveforms with 16.0362 MHz sampling rate. Real data for the airborne experiment were collected during an ESA aircraft mission in the Gulf of Finland on November 11, 2011 [21]. The airborne platform was equipped with the GOLD-RTR receiver installed on board the Aalto University's research aircraft SKYVAN Short SC-7. The acquisitions consist of complex delay waveforms recorded at a sample rate of 20 MHz with 1 ms coherent integration time.

A. Comparison Using CYGNSS Data
After simulation of the set of complex delay waveforms, using wind and mission parameters as recorded in the Level 1 Physical Oceanography Distributed Active Archive Fig. 6. Comparison between the magnitude of the ACFs in the nominal specular point, from real CYGNSS ocean data, from simulated data and from the model. Center (PO.DAAC) science data record [22], the correlation function from the nominal specular point was estimated according to (12) for the simulated waveforms and for the CYGNSS complex delay waveforms. The normalized amplitude of the ACF from simulated and real data is shown in Fig. 6. In the same figure, it is also reported the ACF obtained from the model (2). All plots look quite similar and the estimated correlation time is of the order of 1 ms, in agreement with comments in [23]. As expected, for satellite data, it is not strictly necessary to introduce the temporal evolution into the model because the time-independent formulation provides a sufficiently accurate description of the scattered signal.

B. Comparison Using GOLD-RTR Airborne Data
Driven by the previous simulations and by the model investigation, the comparison among simulation, model, and real data for the airborne experiment requires some more discussion. Two different studies were considered: 1) the first is relevant to the autocorrelation of the complex time series Y (t, τ ) and on the corresponding estimated correlation time; and 2) the second consists in the estimation of the correlation time of the specular point during a long data time series.
In the first case, the POPI0038 dataset was considered. It includes 1 s of acquisitions with geometric conditions reported in Table I (Airborne 1 configuration), with signal received by the Nadir pointing antenna observing reflected global positioning system (GPS) waveforms from the space vehicle with PRN code number 12 and elevation higher than 80 • . For the acquired complex waveforms, the autocorrelation of the time series was calculated for delay bins corresponding to the specular point (bin 33), to the leading edge (bin 20) and to the trailing edge (bin 40). All delay bins are in the range between −1.5 and 1.5 chips, with the delay of the nominal specular point corresponding to τ = 0. Fig. 7 shows plots of 30 simulation trials for the normalized ACF magnitude, for the model ACF and for the ACF estimated from real data. Panels 7(a)-(c), show results for the leading edge, specular point, and trailing edge, respectively. The only additional comment with respect to considerations in Section III is that the real data ACF is very well included in the family of simulations,  close to the mean ACF and to the model ACF. A comparison in terms of the correlation time was also considered. Values estimated from mean ACF (denoted ast cor ), from model equation (denoted as t cor ) and from real data (denoted ast cor ), are reported in Table IV for the considered delay bins. It is worth to observe a good agreement among all values.
In the second case, the analysis of the correlation time was carried out for the bin corresponding to the nominal specular point, using two different datasets, POPI0038 and POPI0018, over 100 s of acquisition, in comparison with the correlation time estimated from simulation and the correlation time evaluated from the model. During acquisitions, small variations of the measured wind speed were registered, ranging between 8 and 10 m/s. Fig. 8(a) and (b) show results for the two datasets, with model and simulated wind speed set to a constant value of 9 m/s. It is interesting to note that data from real acquisitions and from simulations exhibit close patterns with comparable variance and mean value corresponding to the model data. Dataset POPI0038 shows a nearly constant trend for the correlation time whereas dataset POPI0018 shows some variability, in the range 6-4 ms. Such variations cannot be caused by changes in wind speed and direction because values are kept constant in the simulation and model. Reasonably they can be associated with geometry variations in the forward reflection path, that changes along the trajectory. The analysis was extended to other datasets to explore the consistency between simulations and real data. The model data were replaced by the mean over 30 simulated realizations and the correlation time was examined in Fig. 9(a)-(f) for datasets reported in the caption.
To understand the effect of the wind speed on the correlation time, a simplified model was considered by neglecting the ACF of the compound Gaussian process ρ G (nT )e −μ|n|T in (2). This is equivalent to the discrete version of the correlation function derived in [19]. A comparison among correlation time estimated from real data and from model, with and without the term corresponding to sea surface correlation, is reported in Fig. 10. The evidence is that the effect of sea surface  decorrelation is a downward shift of the correlation time, providing a better match with the real data. To assess the quality of the model, the root mean square error (RMSE) between data and model, with and without the sea surface correlation term, is calculated. In (13)t cor i is the correlation time evaluated on the experimental data for the i th second of acquisition, t cor i is the corresponding correlation time estimated on model and N is the number of seconds of acquisition. Results are reported in Table V. As expected, including the sea surface motion in the correlation model provides a better description of the correlation time. Starting from the GNSS reflected signal model presented in the companion work [11], this article provides a procedure for simulation of the GNSS-R signal scattered off the ocean surface, in which temporal correlation properties due to geometry variations as well as the sea surface motion are both incorporated. Stochastic sea surface dynamics are governed by a compound-Gaussian process which describes the time evolution through a 2-D probability distribution arising from a population model for the number of active scattering facets illuminated by the transmitted signal and reflecting in the receiver direction. The simulation procedure, which is based on AR filtering of Gaussian processes, generates a sequence of temporally correlated random variables with the exact dynamics of the scattered signal, with the only limitation that the shape parameter of the Gamma distribution must be semi-integer. The additional computational cost of the procedure with respect to generation of independent samples is negligible. Comparisons among model prediction, simulated time series and real data provide an indirect demonstration of the model validity, in addition to act as a validating procedure for the simulation method. Results confirm that temporal decorrelation in satellite platforms is dominated by geometry variations with correlation function that fades out on a time scale of the order of some milliseconds. For airborne-based receivers, the observation of the sea surface decorrelation discloses some dependence of the correlation time length on wind speed as well as additional features due to variations of the observational geometry that were never observed before.

APPENDIX A
In this Appendix, we recall some properties of Gaussian and Gamma random variables and introduce notation that will be used in the article.
The sequence G(n) is exponentially correlated with autocorrelation Property 4: Let G 1 (n), . . . , G L (n) a vector of independent sequences with G i (n) generated by AR filtering as in (17). The pdf of the sequence is (L/2, 2σ 2 ) with correlation coefficient ρ Z (m) = ρ 2|m| .