Intercomparison of Electromagnetic Scattering Models for Delay-Doppler Maps Along a CYGNSS Land Track With Topography

A comparison of three different electromagnetic scattering models for land surface delay-Doppler maps (DDMs) obtained from global navigation satellite system reflectometry (GNSS-R) along a Cyclone Global Navigation Satellite System (CYGNSS) track in the San Luis Valley, Colorado, USA, is presented. The three models are the analytical Kirchhoff solutions (AKS), the Soil And VEgetation Reflection Simulator (SAVERS), and the improved geometrical optics with topography (IGOT). Common inputs to the three models were defined by using field samples of soil moisture and texture, soil surface roughness measurements, and a digital elevation model (DEM). The resulting peak reflectivity profiles of the models and the CYGNSS data all had a range of 10 dB along the selected track, mainly due to the influence of topography. The reflectivities obtained from all three models agreed with one another within 2.4 dB along the full length of the track. The models also showed general agreement with the corresponding CYGNSS data, although the modeled profiles were higher than CYGNSS Science Data Record Version 3.1 by an average of 5 dB and also smoother. Additional characterization of fine-scale surface roughness is identified as an area for future work to improve model fidelity. An intercomparison of DDM structure for three selected acquisitions is also provided.


I. INTRODUCTION
R EFLECTIONS of satellite signals from land surfaces are sensitive to a variety of biogeophysical parameters of interest for environmental monitoring from space [1], [2], [3]. With the proliferation of spaceborne global navigation satellite system reflectometry (GNSSR) missions and experiments in recent years, there is a growing need for the development and validation of electromagnetic scattering models to describe the delay-Doppler map (DDM) data generated by these sensors.
This work provides an intercomparison of global navigation satellite system reflectometry (GNSS-R) delay-Doppler map (DDM) models for a cyclone global navigation satellite system (CYGNSS) [4] track over a validation site in the San Luis Valley (SLV), Colorado, USA. As the validation site has little vegetation, the intercomparison focuses on modeling 1558-0644 © 2022 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission.
the effects of topography, microwave-scale surface roughness, and soil dielectric constant. The intercomparison includes the following three models: • an implementation of the analytical Kirchhoff solutions (AKSs) [5]; • the Soil And VEgetation Reflection Simulator (SAVERS) [6]; • an implementation of the improved geometrical optics with topography (IGOT) method [7]. Each of these models is based on the Kirchhoff integral, and each uses a different approach for its evaluation.
In addition to the three models and simulators identified earlier and studied in this work, another end-to-end simulator for spaceborne GNSS-R land applications that includes the effects of topography is SIM4Land, which grew out of the European Space Agency (ESA)-funded GNSS-R Assessment of Requirements and Consolidation of Retrieval Algorithms (GARCA) and GNSS REflectometry, Radio Occultation, and Scatterometry (GEROS) projects. Matchups between GARCA/ GEROS-SIM4land and TechDemoSat-1 (TDS-1) GNSS-R data, including a track with high topographic variation, have been reported in [8].
The methodology of the model intercomparison and the corresponding validation site are described in Section II. Results are presented in Section III and further discussed in Section IV. A preliminary version of this work was presented in July 2021 [9].

II. METHODOLOGY A. Summary of Models
Given the near-specular observations of GNSS-R systems, the physical optics approximation (or Kirchhoff approach) is expected to provide reasonable accuracy for predicting the scattering of a rough land surface under low vegetation conditions. All the three models considered in this work are fundamentally based on the Kirchhoff approximation for surface scattering. The models differ in their description of land surface roughness properties.
Land surface roughness can occur over a variety of length scales. Roughness components whose root mean square (rms) heights are very small compared with the electromagnetic wavelength are referred to as "microwave scale" roughness in what follows; such roughness will occur only up to a related horizontal length scale for a given surface. The scattering caused by roughness on greater horizontal length scales is characterized mainly by its surface slope statistics. For the electromagnetic wavelengths of GNSS-R in the L band, the horizontal scale of this transition is on the order of 1 m, for which typically millimeter-to centimeter-scale heights are used to represent the microwave scale roughness rms height. At the other end of the spectrum, the horizontal resolution of the available digital elevation model (DEM) represents the scale above which surface heights and slopes can be represented deterministically rather than statistically. This deterministic regime will be called topographic roughness. In between the microwave and topographic scales lies an intermediate scale where scattering depends on slopes that cannot be resolved by the DEM. Roughness on all of these scales can impact specular scattering and should therefore be characterized in any model to be applied. Note that information on intermediate-scale roughness is frequently unavailable, so that its effects can often be considered as a tuning parameter in the modeling process.
The wide range of length scales involved motivates approaches that decouple prediction of the surface specular scattering with respect to the roughness length scale. The three models compared in this article approach this decomposition in distinct ways. It is also noted that specular reflections from Earth's land surface can include coherent contributions, where reflections from many surface points add constructively at the receiver, and also incoherent contributions, where reflections have random phase, with coherent contributions being more likely in cases with extremely low surface roughness. The three models considered also have distinct means for treating coherent contributions. In the special case of a smooth planar surface, each model reduces to the Friis equation for reflection, as given, for example, in [10, (5.123)] for normal incidence.

1) Analytical Kirchhoff Solutions:
In the AKS approach [5], the land surface terrain is represented by a three-scale surface model where f 1 (x, y) and f 3 (x, y) represent the microwave and topographic scales, respectively, with the former described statistically using parameters from in situ profile measurements, and the latter obtained from a 30 m digital elevation model (DEM). The DEM elevations are used to construct 30-m tilted planar patches, whose slopes are determined from the derivatives of f 3 (x, y). The f 2 (x, y) profile represents an intermediate scale of roughness, which is alternatively called fine-scale topography. Recently, light detection and ranging (lidar) measurements have been taken from which f 2 (x, y) can be reconstructed deterministically in the future. In this work, f 3 (x, y) is represented by deterministic planar patches, whereas both f 1 (x, y) and f 2 (x, y) are treated statistically. Let f 12 (x, y) = f 1 (x, y) + f 2 (x, y), which is the combination of microwave roughness and fine-scale topography. Stochastic descriptions for f 1 (x, y), f 2 (x, y), and f 12 (x, y) are given in [5]. Salient features of the analytical Kirchhoff solutions (AKS) approach include the following.
• Analytical expressions are derived based on the Kirchhoff integral for both coherent and incoherent waves. • Monte Carlo simulations are not required, making it computationally efficient. • Analytical solutions are given in terms of the spectrum of f 12 (x, y). In this approach, there is no need to divide microwave roughness and fine-scale topography, and the surface spectrum derived from lidar measurements can be incorporated directly. For both coherent and incoherent waves, the AKS approach gives results that are indistinguishable from the numerical Kirchhoff approach (NKA), which is a brute force accurate method that carries out the Kirchhoff integral directly using 2 cm discretization. Results show that f 2 (x, y) has significant effects [5].
To calculate the coherent waves for a certain area, the coherent field is first obtained from each 30 m DEM patch while accounting for the impact of the f 12 roughness on the patch. Then, the total coherent field is obtained by the complex summation of the coherent field over patches, and the absolute value squared is finally taken to find the coherent intensity of the area. For the incoherent waves, the incoherent intensity from each 30 m DEM patch is computed. Then, the total intensity is obtained by summation of all the contributions from patches. With calculated coherent and incoherent intensity, the total scattering from the area can be obtained [5].
2) SAVERS: SAVERS is based on the original formulation in [11], which was designed to simulate low-altitude receivers. It was updated as described in [6] to include the topography of the illuminated area, which cannot be neglected at satellite altitude. To account for topography, a DEM is considered, each of whose elements is represented by a facet with its individual orientation, above which the roughness at wavelength scale is superimposed.
SAVERS implements the integral bistatic radar equation [12], independently evaluating the incoherent diffuse and coherent near-specular scattering components. The former is computed through the advanced integral equation method (AIEM) [13], whereas the latter is simulated using the approach described in [14], which relies on the definition of a coherent normalized radar cross section (NRCS), associating a prescribed scattering pattern with the DEM facets. The beamwidth of the pattern of the coherent component is a parameter, denoted by β, that depends on the system geometry and configuration, such as frequency and distance. Since it drives the directivity of the quasi-specular scattering component [14], it can be connected to the effect of the intermediate-scale roughness on the patch NRCS angular pattern. In the SAVERS simulations performed in this study, the parameter β has been set to 0.03 • , which is the same order of magnitude of the one retrieved in [14] for a satellite geometry.
The contributions of each DEM facet are combined through the radar equation considering the Cyclone Global Navigation Satellite System (CYGNSS) antenna pattern and the transformation between the local incidence and observation angles and polarization state going from the global reference frame linked to the instrument to the local facet reference frame and vice versa.
SAVERS also includes a module for simulating vegetated areas, evaluating the scattering from forest and agricultural vegetation and setting the relevant geometric features of the vegetation elements through growth models. The modeling implements the radiative transfer equation, considering media constituted of randomly distributed scatterers representing different vegetation elements (namely, leaves, branches, and trunks) [6]. In the work reported in this article, the effects of vegetation are not included, since the surfaces near the validation site were nearly bare at the time of data acquisition.
3) Improved Geometrical Optics With Topography: The IGOT method, first described in [7], follows closely the approach of [12] with the assumption that the Rayleigh roughness parameter is large on the horizontal scale of the footprint of each delay-Doppler bin [15] so that coherent contributions vanish. Unlike [12], however, the improved geometrical optics with topography (IGOT) surface height is considered to be not a purely random field but rather decomposed into a deterministic part obtained from a DEM and a random part representing the residual height between the DEM and the surface. The random part of surface height is further decomposed into a long-wave process and a shortwave process, following the improved geometric optics model [16]. A derivation from first principles that includes the shortwave diffraction process is provided in the video presentation of [17]. Thus, the IGOT model is parameterized by the following: 1) DEM heights and gradients representing large-scale topography; 2) an rms slope characterizing the roughness scale between the DEM resolution and the geometric optics cutoff; 3) an rms height accounting for attenuation of the geometrical optics scattering due to shortwave diffraction. 4) Discussion: The three models described have similar approaches for describing the microwave-and topographicscale roughness, but all represent the intermediate-scale land roughness and treat coherent scattering contributions in distinct ways. The AKS model represents this roughness statistically in terms of an associated surface covariance function and computes the Kirchhoff expression for the surface NRCS without resorting to the geometrical optics approximation. The SAVERS model similarly requires knowledge of the surface covariance function for the diffuse incoherent component and relies on the beamwidth parameter introduced in Section II-A2 for the evaluation of the near-specular contribution. The IGOT model in contrast assumes that geometrical optics holds for all surface DEM patches, which implies an assumption regarding the amplitude of the intermediate-scale roughness rms heights, so that description only of the rms slopes of the patch roughness is required. The intercomparisons to be shown will provide insight into the applicability of these assumptions.

B. Calibration/Validation Site
The SLV of South Central Colorado, USA, was selected as the location of the first calibration and validation (cal/val) site for land applications of the CYGNSS mission, in part due to the elevation being sufficiently high to experience freeze-thaw cycles but not so high as to exceed the operating envelope of the instrument. Located at the headwaters of the Rio Grande near the northern limit of the CYGNSS coverage zone and surrounded by mountainous terrain, the SLV is generally flat and sandy with croplands as the dominant land cover.
Two Soil moisture Sensing Controller and oPtimal Estimator (SoilSCAPE) in situ wireless sensor networks (WSNs) named Z1 and Z4, whose locations are shown in Fig. 1, were installed in the SLV in late October 2019 [18]. The SoilSCAPE project aims to provide surface-to-depth observations of soil moisture and temperature on a local scale with optimal sampling [19], [20]. The locations of the two WSNs were selected such that both represent similar weather and climate conditions. The site Z1 is a flat open pasture, whereas Z4 samples more hilly terrain. Additional details are available at the project data server https://soilscape.usc.edu.

C. Track Selection
Criteria for selection of a CYGNSS track for the model intercomparison included the following: 1) time of acquisition around October 25-27, 2019, when field samples were collected; 2) readings from all SLV SoilSCAPE temperature sensors above freezing; 3) at least one acquisition located within 5 km of Z1 or Z4. Item 1) was included so that soil moisture from oven-drying of the field samples could be used rather than soil moisture from the SoilSCAPE sensors, whose calibration was delayed due to restrictions related to COVID-19. Additionally, the water content of the vegetation at the time of the field samples was observed to be low, e.g., as seen in Fig. 2. The monthly average of leaf area index (LAI) from the moderate resolution imaging spectroradiometer (MODIS) mission along CYGNSS tracks over the SLV for October 2019 ranges between 0.2 and 0.6. This is an indication of very low vegetation, which is aligned with the focus of this study on sensitivity to topography rather than to land cover.
After analysis of CYGNSS data near these field samples, the track around 2019-10-28 14:04:58.5 UTC from channel 3 of spacecraft 2 was selected for this study, as this track satisfied the three aforementioned selection criteria. The track, whose reported specular points are plotted over an elevation map in Fig. 1, passes over Z1, which has coordinates of (37.190 • , −105.992 • ). The start and end of the track were determined by requiring CYGNSS signal-to-noise ratio (SNR) to be greater than 2 dB.

D. Soil Dielectric Constant
The soil moisture content measured from oven-drying of the Z1 field samples was 0.0259 m 3 m −3 . The clay fraction measured from the corresponding analysis of soil texture was 18%. With these values of soil moisture and clay fraction,  Special test equipment for measuring microwave-scale surface roughness: laser range finder that slides along a spirit level mounted on a pair of tripods. Photo by Amer Melabari. a soil dielectric constant of 2.987 + 0.173 i was calculated by the Mironov model [21], [22]. The same value was used by all three GNSS-R DDM models over the entire footprint of the selected track.

E. Soil Surface Roughness
A laser range finder mounted on a spirit level supported by a tripod at each end, as shown in Fig. 3, was used to measure microwave-scale soil surface roughness along the baseline of the level at a total of 13 locations and orientations around Z1, excluding a 14th measurement over a trench. Of the 13 characterizations, seven used a baseline of approximately 0.5 m, and the other 6 had a baseline of approximately 1 m. Since the laser signal could reflect from vegetation, care was taken to avoid vegetated areas. The rms of each characterization was computed relative to the mean for flat areas and relative to a linear regression for hilly locations. As shown in Fig. 4, a wide range of rms values was observed. In particular, the values of rms surface roughness measured during fieldwork at Z1 in the ascending order of magnitude. To further assess microwave-scale surface roughness, the field bare_soil_roughness_retrieved was extracted from the soil moisture active passive (SMAP) Level-3 product Radar Global Daily 3 km equal area scalable Earth (EASE)-Grid Soil Moisture, Version 3 [23] for the week beginning July 1, 2015. This field is plotted against distance from Z1 in Fig. 5, which shows that points farther away from Z1 have a roughness value around 2 cm and sometimes higher, whereas the closest points have a roughness value around 0.5 cm. Using a radius of 3 km, the average roughness is 0.45 cm. A 5 km radius gives an average roughness of 1.03 cm, and an 8 km radius yields an average roughness of 1.48 cm. The increase in roughness with distance from Z1 is expected, since Z1 is located in a relatively flat area that is surrounded by mountains. Note that the SMAP-derived roughness parameter is associated with a backscatter radar measurement at 40 • incidence angle, and therefore may not be completely applicable to describing near-specular surface scattering. The range of values nevertheless provides some level of insight into roughness properties in this region.
In addition to rms surface roughness retrieved by the SMAP radar, rms surface roughness from SMAP ancillary data [24] was examined, and the values over SLV were found to be about 1 cm. These values are based on SMAP passive data and the tau-omega model, which may not accurately reflect actual surface conditions.
In light of the data presented earlier on the microwave-scale soil surface roughness at Z1, the rms height parameter was set to 1 cm for all three models. Since the AKS method represents f 1 not only by an rms height parameter but also with a correlation function, an exponential correlation function with a correlation length of 10 cm was applied. The same correlation length was used in SAVERS for the incoherent component. These parameters were held constant over the entire footprint of the track selected for the reflectometry model intercomparison study of Sections III-A and III-B.

F. Surface Topography
All the three models used the 1-arcsecond Shuttle Radar Topography Mission (SRTM) DEM having approximately 30-m horizontal sampling for elevation. The models also used a common set of gradients estimated from the 1-arcsecond SRTM DEM using linear least squares with a Hann window of size 15 by 15 samples. At the latitude of the SLV, this window corresponds to a region on the ground of approximately 370 by 460 m. The window size needed to be large enough to reduce the effects of speckle noise from the SRTM product on the estimate of the gradients but not so large as to degrade their resolution. The speckle noise, which is the dominant SRTM error source, is significantly correlated for separations of 50 m, with correlation dropping off to a plateau for separations greater than approximately 100 to 200 m [25]. It is noted that, since the scattering coefficient in the three models being studied is sensitive to surface gradients and since the 15 by 15 sample Hann window of the filter used to estimate the surface gradients has a resolution of roughly 300 m, the upper limit of the intermediate scale of surface roughness in this study is roughly 300 m. A comparison of results using two different window sizes for gradient estimation is provided in Section III-C.
The AKS and IGOT models both ran at the 1-arcsecond resolution of the SRTM DEM, while SAVERS resampled to a spacing of 9 arcsecond, or approximately 300 m.
As described in Section II-A, all three models include parameters for characterizing surface roughness in the intermediate scale between the resolution of the DEM and the microwave scale. For this study, the AKS model represented f 2 by a Gaussian correlation function with an rms height of 5 cm and a correlation length of 125 times the rms height. The SAVERS simulation sets β to 0.03 • . The IGOT model was run with a relative rms slope of 0.4 • . These parameters were held fixed over the entire footprint of the selected track.

G. Antenna Pattern
In the interest of reproducibility from publicly available data, all three models used an isotropic antenna pattern. The convention for antenna gain is described in Section II-I.
A comparison of results with isotropic and anisotropic (i.e., directional) antenna patterns is provided in Section III-C.

H. CYGNSS Data
Version 3.1 of the CYGNSS level 1 science data record was selected for comparison with the model results [26], [27]. This was the most recent version at the time of writing. In particular, calibrated CYGNSS DDMs of bistatic radar cross section (BRCS) with units of square meters were obtained from the Network Common Data Form (NetCDF) variable brcs and converted to reflectivities using the following procedure.

I. Conversion to Reflectivity
To facilitate meaningful comparisons among the model results and the CYGNSS data, a common convention for antenna gains and DDM units needed to be defined. Typically, DDMs are expressed in one of the following ways.
• Power in units of watts. Whereas this convention is valid for both coherent and incoherent reflections, it has no interpretation as a physical property of the reflecting surface. • BRCS σ in units of square meters. This convention has physical significance for incoherent reflections, such as those from ocean surfaces, and it is used for the calibrated CYGNSS Level 1b (L1b) product [28]. • Normalized bistatic radar cross section (NBRCS) σ 0 in units of square meters per square meter. Having physical significance for incoherent reflections, this convention is used to average multiple BRCS bins together over their corresponding effective areas, as in [28]. • Surface reflectivity [10, Secs. [2][3][4][5][6][7][8]. This convention has physical meaning for coherent reflections, and it has traditionally been used for land applications, since the reflections in early tower-based and airborne GNSS-R experiments were primarily coherent. It is also being included in the CYGNSS land product [29]. For this study, the last alternative surface reflectivity was adopted for consistency with previous work. A convention for surface reflectivity is developed in the following to permit convenient comparison among the three models and CYGNSS data products. The reflectivity parameter defined in this article can be applied for both coherent and incoherent scattering and should be intended as an equivalent reflectivity rather than the actual reflectivity of a specular surface.
Each of the three forward models considered in this study can be expressed in the form of the bistatic radar equation given as [28, eq. (1)] x,y 2 2 τ ;x,y S 2f ;x,y dx dy (2) where superscript FM distinguishes forward model parameters from the corresponding CYGNSS L1b calibration parameters. The forward model antenna gains will be set equal to the CYGNSS gains in the following development. Here, P ĝ τ ,f is modeled power received after coherent processing at delayτ and Dopplerf by way of scattering of the global navigation satellite system (GNSS) source from the rough surface, P T is transmit power, λ is wavelength, (x, y) ∈ A are the variables of surface integration, A is the region of diffuse scattering, G T x,y and G R x,y are transmit and receive antenna gain patterns, respectively, as a function of (x, y), R T x,y is the distance from the transmitter to the surface at (x, y), and R R x,y is the distance from the surface at (x, y) to the receiver, σ 0 x,y is NBRCS at the bistatic scattering geometry defined by the transmitter position, the surface at (x, y), and the receiver position, and τ ;x,y and Sf ;x,y are the Woodward ambiguity functions (WAFs) in delay and Doppler, respectively.
For the AKS, (2) represents the incoherent component only. Since the coherent component of the AKS was found to be negligible for the track and the input parameters selected in this study, the coherent component is not included here. (Nevertheless, we note that the convention developed in the following can be extended in a straightforward manner to include a coherent component if necessary by replacing (2) with the corresponding equation for coherent power in terms of reflectivity.) In SAVERS, (2) is always applicable, since it includes both the diffuse and near-specular component. Likewise, under the assumptions of IGOT, (2) is obtained from first principles [7, eq. (42)].
Calibrated CYGNSS DDMs are provided in the NetCDF variable brcs with units of square meters. From [28, eq. (4)], these L1b BRCSs are computed as where superscripts L1a and L1b have been added to emphasize certain factors that are specific to the calibration, as distinguished from the corresponding factors in the forward models, and where subscript SP denotes evaluation at the specular point reported by the L1b calibration. Here, P g,L1â τ ,f is the measured power DDM in units of watts from the Level 1a (L1a) calibration.
Since (3) represents a CYGNSS version-specific convention for conversion from units of watts to units of square meters, the same convention can also be used to convert a forward model result from units of watts (whether coherent or incoherent) to units of square meters for comparison with a particular version of L1b CYGNSS data Substituting (2) into (4), we obtain x,y 2 2 τ ;x,y S 2f ;x,y dx dy. (5) Approximating the range losses and the antenna gains inside the integral by their values at the specular point (i.e., assuming the antenna pattern is isotropic over A), we find x,y 2 τ ;x,y S 2f ;x,y dx dy. (6) In the following, we adopt the conventions: so that (6) becomes Since (8) is independent of calibration-specific factors, this convention allows a single run of a forward model to be compared with multiple versions of CYGNSS L1b data. The convention also eliminates the need to favor any particular version of effective isotropic radiated power (EIRP) and receive antenna gain in the forward model run when comparing with multiple versions of CYGNSS L1b data. Finally, to convert BRCS to reflectivity, we use the relationship both for forward model results and for CYGNSS data. This relationship follows directly from (4) and the Friis formula for reflections [10,]. Here, the relationship is effectively independent of CYGNSS version, since the relative difference in the reported ranges of the specular point between two versions is small.

A. Along-Track Analysis
A comparison of peak reflectivity along the selected track is shown in Fig. 6. Each plotted point represents the maximum value of the DDM generated by the corresponding model or CYGNSS product, converted to units of reflectivity. The horizontal axis is expressed in units of longitude so that each point can be associated with the location of its corresponding reported specular point in Fig. 1. The plot also includes surface elevation at the specular points.
The reflectivity profile of all three models and of the CYGNSS data is highest over the valley floor and decreases in the mountainous terrain on either side, and to a lesser degree over some hilly terrain in the middle, with a range of 10 dB.
All three models agree with one another to within 2.4 dB over the entire track. However, all three models overestimate version 3.1 of the CYGNSS data. In particular, the average of the three models is 5.0 dB higher than version 3.1 when averaged along the track. Furthermore, all three models appear to be generally smoother than the CYGNSS data. Comparison of peak reflectivity from the three models and the corresponding CYGNSS level 1 science data along the selected track.

B. Selected DDM Matchups
A comparison of DDMs for the three acquisitions identified in Table I from the selected track is shown in Fig. 7. The middle acquisition was chosen to be located as close as possible to Z1, where the soil samples were collected and microwave-scale surface roughness measurements were made. The other two acquisitions were selected for their locations on opposite sides of the valley. As shown in Fig. 7, the west acquisition has a positive Doppler tail. The middle acquisition is relatively compact in delay with a weak positive Doppler tail. The east acquisition has a negative Doppler tail. These structures are represented in both the CYGNSS data and the three models. Interpretation of the observed structures is provided in Section IV. Additionally, with the sole exception of Fig. 7(a), the location of the maximum value of each model DDMs is within one bin of the location of the peak of the corresponding CYGNSS DDM for each of the three selected acquisitions.

C. Sensitivity Analyses
Results of the sensitivity studies referenced in Sections II-F and II-G are provided in the following. The purpose of these studies is not to compare sensitivities among the three models but rather to assess their general sensitivity to specific assumptions underlying the model intercomparison in Fig. 6. Only results from IGOT are shown due to space limitations, although similar results were also obtained with SAVERS. (No corresponding sensitivity results were generated with AKS.) The results are discussed in Section IV.
A comparison of peak reflectivity from the IGOT model with two different window sizes for gradient estimation and Comparison of DDMs from the three models and from the corresponding CYGNSS level 1 science data for the selected acquisition. Each DDM is normalized to its peak value to facilitate intercomparison of delay-Doppler structure. The units of the color map are decibels relative to the peak.
the corresponding CYGNSS level 1 science data along the selected track is shown in Fig. 8. As seen in the figure, the larger window provides a range that is more than 1 dB greater than that of the smaller window.
A comparison of peak reflectivity from the IGOT model with anisotropic and isotropic antenna patterns and from Comparison of peak reflectivity from the IGOT model with two different window sizes for gradient estimation and the corresponding CYGNSS level 1 science data along the selected track. the corresponding CYGNSS level 1 science data along the selected track is shown in Fig. 9. Here, the anisotropic antenna pattern was taken to be the average over the eight CYGNSS observatories of the receive antenna patterns that were used to generate CYGNSS Version 2.1 data products [30]. The isotropic and anisotropic results are within 0.64 dB of each other.

IV. DISCUSSION AND CONCLUSION
Three different electromagnetic scattering models for DDM of GNSS-R from land have been compared along a CYGNSS track with varying topography crossing a validation site. Insofar as possible, a common set of input parameters, input datasets, and output conventions was used to facilitate the comparison, but the three approaches used distinct methods for describing the intermediate-scale roughness for which ancillary information was not available. The models are particularly sensitive to soil surface roughness over a wide range of scales, which is also consistent with results in [5,Sec. 8.5].
At the topographic scale, variation in surface heights and surface gradients computed from SRTM alone is able to account for the most relevant part of the range of reflectivity observed along the transect, since all other input parameters were held fixed. Although the models were unable to resolve the finer-scale features of the CYGNSS reflectivity profile, particularly over the hilly region in the middle of the valley, preliminary investigation with a high-resolution airborne lidar survey from 2020 suggests that the problem may lie not with the models themselves but rather with the low resolution of the surface gradients estimated from SRTM.
At the microwave scale, attenuation due to shortwave diffraction can introduce a bias into the results. A comparison of the sensitivity of the bistatic scattering coefficient of AKS and IGOT to microwave-scale surface roughness was performed in [5,Sec. 8.5], in which IGOT was represented by the model referred to as geometric optics with attenuation (GO-Att). As rms height was increased from 1.0 cm to 3.5 cm, the response of both AKS and IGOT dropped by roughly 14 dB. The sensitivity analysis of SAVERS to surface roughness was reported in [6], where the results showed a reflectivity decrease of approximately 12 dB when the roughness increased from 0.5 to 3.0 dB at an incidence angle of 30 • . The analysis was repeated for the selected SLV track by SAVERS and IGOT, and it was found that increasing the rms height from 1.0 to 2.0 cm decreased the reflectivity of both models by approximately 4 dB, thereby approaching the CYGNSS observations. Thus, uncertainty in microwave scale surface roughness may account for some of the offset observed between the models and the CYGNSS data, and further investigation is needed. Characterization of microwave-scale surface roughness to sufficient accuracy over sufficiently large areas required for model validation against real data is a challenge.
Additional effort is also needed to characterize soil surface roughness in the intermediate regime between the topographic scale and the microwave scale. Some characterization may be possible by comparing the peakedness of the synthesized DDM with that of the satellite data. Additionally, work is in progress to make use of the 2020 lidar survey, since it has provided a DEM resolved at approximately 30 cm sampling and approximately 5 cm vertical uncertainty. Thus, the resolution of the lidar product is about two orders of magnitude better than that of the SRTM in each horizontal dimension.
Results are also sensitive to gradients at the topographic scale. As seen in Fig. 8, a change in the window size for gradient estimation alone can have an effect of more than 1 dB.
Error due to approximation of the antenna pattern as isotropic (i.e., nondirective) over the glistening zone within each acquisition, which was made to expedite the intercomparison runs and which is described in detail in Section II-I, appears to be low relative to errors from the other sources discussed earlier. The differences observed in Fig. 9 are well within 1 dB.
Another source of uncertainty in the comparison is inhomogeneity of soil moisture and soil texture across the SLV. The comparison study assumed a fixed value for soil dielectric constant. However, a cursory look at satellite imagery of the valley reveals a complex patchwork of croplands in various states of irrigation from center pivots, surrounded by the nonagricultural lands and mountains that form the drainage system at the headwaters of the Rio Grande. In situ soil moisture sensor networks are able to cover only a tiny fraction of the valley, and they generally cannot be installed on private cropland with irrigation pivots.
Uncertainty in the various calibration parameters for Version 3.1 CYGNSS Level 1 (L1) data is currently being evaluated by the science team [31]. Version 3.1 includes improvements to the CYGNSS science antenna gain patterns and a correction for coarse quantization effects from the onboard digital processor [26]. By comparison, peak reflectivity of version 3.0 [32] is an average of 4.2 dB higher across the selected track than peak reflectivity of Version 3.1. Peak reflectivity of version 2.1 is an average of 1.8 dB higher across the selected track than that of version 3.1.
Future comparison efforts could include forthcoming CYGNSS land products [29] and calibrated raw intermediate frequency (IF) products, with the latter allowing for shorter noncoherent integration times, finer spatial sampling, and direct measurement of signal coherence.
As noted in Section III-B, both positive and negative Doppler tails are observed in the DDMs of the three models and the CYGNSS data, similar to the Doppler tails previously observed with TechDemoSat-1 [6]. Analysis indicates that the observed Doppler asymmetries arise from structures in the underlying topography. In particular, since the specular point moves from west to east for the selected track, scatterers lying to the west of the specular point have Doppler frequencies less than the Doppler of the specular point, whereas scatterers lying to the east of the specular point experience positive relative Doppler. For the west acquisition, scatters with positive Doppler lie over the valley where scattering cross sections are relatively high, whereas scatterers with negative Doppler lie over the mountains to the west of the valley where scattering cross sections are relatively low. This accounts for the positive Doppler tail seen in the west acquisition. For the east acquisition, the situation is reversed with scatterers of positive Doppler located over attenuating mountains to the east and those of negative Doppler located over the more strongly reflecting valley floor, thereby accounting for the negative Doppler tail.
As observed in Section III-B, the location of the peak of the modeled DDMs is generally within one bin of the highest pixel of the corresponding CYGNSS DDM. Exact agreement is not necessarily expected, since an accurate registration of the delay coordinate of the CYGNSS DDMs requires the estimation of various clock biases, which is outside the scope of the present work. A more detailed assessment of delay and Doppler location could be a topic for future work.
The general agreement among the models used in the comparison shows that each of the approaches can be considered applicable to the extent that the intermediate-scale roughness values implicitly assumed can be related to other sources of information, thus avoiding the need to tune the related model parameter. The success of the IGOT model, for example, could be taken as implying that the surface is sufficiently rough in this region as to result in the geometrical optics approximation holding, so that only knowledge of rms slopes is required and so that coherent contributions can be neglected. These assumptions, however, may not hold in other regions where intermediate-scale roughness and topographic roughness are smaller. Further assessment of these questions must await more detailed assessments using the airborne lidar data and studies at sites having a variety of roughness characteristics.
Leila Guerriero (Member, IEEE) received the Laurea degree in physics from Sapienza University, Rome, Italy, in 1986, and the Ph.D. degree in electromagnetism from Tor Vergata University, Rome, in 1991.
Since 1994, she has been a Permanent Researcher with Tor Vergata University, where she is currently an Associate Professor holding the courses on Earth satellite observation and on geoinformation. Her research interests include modeling microwave scattering and emissivity from agricultural and forested areas. She participated in several international projects, among them the ESA projects Soil Moisture and Ocean Salinity Satellite, Development of SAR Inversion Algorithms for Land Applications, Use of Bistatic Microwave Measurements for Earth Observation, and Microwaves Observation Satellite-Companion Satellite (SAOCOM-CS) Bistatic Imaging, Radiometry and Interferometry over Land. Later, she has been involved in the modeling of GNSS-R signals for ESA projects and in the European FP7 and H2020 Programs.
Dr. Guerriero is a member of the Permanent Steering Scientific Committee of MicroRad. She is the Secretary and the Treasurer of the GRSS North-Central Italy Chapter.
Erik Hodges (Member, IEEE) received the B.S. degree in electrical engineering from the University of California Los Angeles, Los Angeles, CA, USA, in 2019. He is currently pursuing the Ph.D. degree in electrical engineering with the University of Southern California, Los Angeles, CA, USA.
His research interests include GNSS-R and electromagnetic models for land surfaces.
Dr. Hodges is a member of the CYGNSS Science Team and the IEEE GRSS.