A full-polarization GNSS-R Delay-Doppler-Map (DDM) simulator for bare soil freeze/thaw process detection

The soil freeze/thaw processes have a significant impact on the surface energy and moisture balance, which play a key role in ecosystem diversity and productivity. The air-/space-borne Global Navigation Satellite System (GNSS)-Reflectometry (GNSS-R) is a bistatic radar with receiving power as the Delay Doppler Map (DDM), which may monitor soil freeze/thaw processes. However, the scattering mechanism of monitoring soil freeze/thaw processes by GNSS-R DDM is not clear. In this paper, it is the first time to simulate full-polarization GNSS-R Delay-Doppler-Map (DDM) for bare soil freeze/thaw process, not only the linear polarization but also the circular polarizations. As the bare soil freeze/thaw process occurs, the corresponding DDM variations are able to present by this simulator. Other geophysical parameters, such as soil moisture and surface roughness, are also two important parameters affecting the final GNSS-R receiver power and their effects on DDM are also presented by the simulator. This simulator can be a potentially efficient tool for data analysis and interoperation of GNSS-R received power as well as the in situ experimental design.


Introduction
In the past two decades, Global Navigation Satellite System (GNSS)-Reflectometry (GNSS-R) has emerged as a new promising remote sensing technique. The multipath observables of geodetic or geophysical GNSS receivers can be used for geophysical parameters' retrieval (Larson 2016), such as near surface soil moisture (Chew 2015;Vey et al. 2016), vegetation water content or growth Wan et al. 2015), snow depth (McCreight et al. 2015) and sea level estimations (Löfgren et al. 2014;Jin et al. 2017). Due to the unique advantages of low cost, small power, wide coverage and high spatial/temporal resolutions, specially designed GNSS-R receivers on airborne or space-borne platforms are able to collect the ocean or land surface reflected signals. They have been widely used in various terrestrial remote sensing, which can be used to retrieve ocean or land surface characteristics (Cardellach et al. 2011;Jin et al. 2011;Zavorotny et al. 2015). Since the first UK-DMC mission, the successively PARIS IoD and the GEROS-ISS have been put forward by ESA. TechDemoSat-1 has launched in July of 2014 and its data are available for community (Unwin et al. 2017). Researchers have used the TDS-1 data for the remote sensing detections of sea ice (Yan and Huang 2016), soil moisture (Camps et al. 2016;Chew et al. 2016), vegetation (Camps et al. 2016), and wetland , while NASA's Cyclone GNSS (CYGNSS) mission has been launched on 15 Dec 2016 (Ruf et al. 2012), it will provide more chances. The transmitters of the GNSS constellations and the corresponding GNSS-R receivers form the bistatic radar mode and the received signals are the Delay Doppler Maps.
Soil freeze-thaw process detection is a very important part of earth system monitoring. About 35% of the total Open Access *Correspondence: xrwu@shao.ac.cn 1 Shanghai Astronomical Observatory, Chinese Academy of Sciences, Shanghai 200030, China Full list of author information is available at the end of the article land area is covered by the frozen soil and permafrost in high latitude and altitude. The understanding of the soil freeze/thaw process has a very significant impact on the surface energy and moisture balance, and it plays a very important role in ecosystem diversity and productivity. Traditionally, investigation of the process is conducted by in situ point measurement, which is labor cost and limited representative. Microwave remote sensing (either radar or radiometer) has been used for the detection (Judge et al. 1997). However, its time and spatial resolutions are out of the scientific requirements. The newly developed GNSS-R technique has the potential to be an efficient remote sensing tool. At present, few of the GNSS-R applications no matter ground based GNSS-IR or the airborne/space-borne are aimed for soil freezethaw process detection, especially the theoretical studies (Wu and Jin 2014). The first GPS-reflected power model is the Z-V model, which takes the GNSS signal structure and the geometry into considerations (Zavorotny and Voronovich 2000). However, it originally aimed at the ocean surface applications. The authors also have used this theoretical model for bare rough surface study, but it is not sufficient to satisfy the real experimental study (Masters 2004). Within the framework of the Land Monitoring with Navigation signal (LEiMON), a simulator of GNSS reflections aimed for bare and vegetated soils (SAVERS) is developed and presented (Pierdicca et al. 2014). The basic framework of SAVERS is Z-V model (Zavorotny and Voronovich 2000); the formulation proposed by Fung and Eom (1983) is employed for the calculations of bare soil coherent scattering properties, while the advanced integral equation model (AIEM) is used for its incoherent scattering calculations (Wu and Chen 2004). A discrete microwave scattering model, i.e., the Tor Vergata model, is employed to characterize interactions between GNSS electromagnetic waves and vegetation cover (Bracaglia et al. 1995).
In this paper, we put forward a fully polarization GNSS reflection model aimed for the bare soil freeze/thaw process detection. Theoretical fundamentals are presented in "Theories and models" section, while the simulation and results are presented in "Simulations and analysis" section as well as discussions in "Discussion" section and finally conclusions are given in "Conclusions" section.

Theories and models
Most of the transmitted signals of GNSS constellation work at microwave L bands, and GNSS-R is essentially a kind of bistatic radar. Therefore, the fundamental of microwave remote sensing is suitable for the development of GNSS-R geophysical parameters' detection. As the geophysical parameters change, it will result in the variation of the permittivity, which will influence the scattering properties of the target objects. In this way, the received GNSS scattering power changes, from which we can used for the geophysical parameters' detection.

Permittivity
Permittivity is the unique character of the target objects. As for the calculations of the bare (either frozen or thawn) soil, dielectric mixing models are employed for the calculations. As for the thawn soil, a four-component (dry soil solid, bound water, bulk water and air) dielectric mixing model is used Hallikainen et al. 1985). While for the frozen soil, one part of ice is added in the simulation model (Zhang et al. 2003). The dielectric properties of soil are influenced by the electric magnetic wave frequency, soil temperature, soil salinity, volumetric content of water in the soil, percentage of the bound water and free water, soil volumetric density, shape of the soil particles and the impurity. The final permittivity of the mixture ε m is the sum of all component parts with their own weights.
where α (a superscript index, not a power exponent) is the shape factor, i implies the mixed component number, and the ith component permittivity is ε i and its volume is V i . As for frozen soil, an empirical function W u = A · |T − 273.2| −B is used for the estimations of the fractions of liquid water and ice in frozen soil, where W u is the content of unfrozen water, A and B are the soil texture parameters, and T is soil temperature.

Bistatic scattering cross-section
The GNSS constellations are the transmitters and the specially designed GNSS-R receivers can form the bistatic radar mode. Different from the emissivity models used in the passive microwave radiometer or the backscattering models used in the active microwave radar, a bistatic scattering model is needed for GNSS-R remote sensing technique. Here, the bisatic advanced integral equation model is employed.

AIEM model
The bare soil surface scattering problem is related to the random rough surface scattering, whose traditional models are the small perturbation method (SPM) and the Kirchhoff approximation (KA). For different roughness scales, Geometrical Optics model (GO) or Physical Optics (PO) models are used. To solve the discontinuous roughness problem of the traditional surface scattering models, Integral Equation Model (IEM) is developed. IEM model can cover very large range of natural surface roughness and it has been widely used in practical applications. However, due to several assumptions and the imperfect expressions of the surface roughness, it has been improved and the advanced Integral Equation model (AIEM) is proposed (Wu and Chen 2004). For simplify, AIEM model can be expressed as the sum of three components: the Kirchhoff term ( σ k qp ), the complementary term ( σ c qp ) and the cross-term of them ( σ kc qp ): where σ • qp is the bistatic scattering cross-section, the subscript p indicates the transmitted polarization and q is the received polarization.

Wave synthesis technique
The traditional microwave scattering models, either the emissivity models or the backscattering models, are oriented to the linear polarizations. To mitigate the influence of ionospheric effects, the transmitted polarizations of the GNSS constellations are circular polarizations, i.e., right-hand circular polarization (RHCP). After reflected from the ground surface, the polarization states of the electromagnetic wave are conversed. Several works have been focused on the study of the polarizations, both the experiment campaign and the theoretical analysis. To get the polarizations at any combinations, wave synthesis technique is employed (Ulaby and Elachi 1990).
where the subscripts and superscript of r and t are the polarization states of received power and transmitted power, respectively, ψ is the orientation angle, χ is the ellipticity angle, and M m is the 4 × 4 real modified Mueller matrix. Using the different forms of normalized stokes vectors of Y m , the corresponding polarization scattering cross-sections at any polarization combinations can be calculated.

GPS scattering power model
The frame of the GPS scattering model is the first GPS scattering model developed by Zavorotny and Vorinovich (2000), although the original model is oriented for ocean wind filed study, but it has been extended to some other applications, such as oil slick study. The output power, i.e., DDM, can be represented as a function of the time delay and of the Doppler shift, is where P T is the transmitted power of the GNSS satellite and its wavelength is represented by , G is the gain and R is the range between scattering point to the transmitter or receiver, while the subscript T or R represents the transmitter or the receiver, and σ • is the bistatic scattering cross-section of the objective targets. The triangular function and sinc function are presented by Λ and S. The bistatic scattering cross-section σ 0 rt in the integration area A is calculated using the above-mentioned models.
It can be seen that Z-V model is the integral form of bistatic radar equation. The flowchart of the fully polarization DDM simulator for soil freeze/thaw process detection is presented in Fig. 1.
There are two modes for the calculations of the model. One is for the calculations of satellite geometry, the position and velocity of GPS constellation and (4) the corresponding receiver are used to get the satellite geometry, while the second mode is for the calculations for bare soil (either frozen or thawn) bistatic scattering cross-sections. By combing the two modes, a fully polarization DDM simulator is got and one of the most important final outputs is the Delay Doppler Map (either frozen or thawn) power. Figure 1 presents the permittivity variations versus soil temperature. As can be seen from Fig. 1, when the soil temperature is below 0 °C, both the real part and the imaginary part increase as the soil temperature increases, but the trends are just opposite as the soil temperature is above 0 °C. One important phenomenon presented in Fig. 1 should be noticed: as the soil temperature changes from below 0 °C to above 0 °C, there is an abrupt change for the soil permittivity (both the real part and the imaginary part). These variations are the fundamentals for soil freeze/thaw process detection using GNSS-R technique. As in our GPS scattering power (DDM) simulations, the carrier frequency is GPS L1 band, while the wavelength is 0.19 cm, and integration time is set as 0.001 s. For the sake of simplification, the antenna gain of the transmitter and receiver is all set 1, while the incident angel is 13°.

Simulations and analysis
The real part of the permittivity manifests the reflection and refraction of the electromagnetic waves in the medium, while its imaginary part represents its attenuation (absorption and conversion). Figure 2 shows that as the soil temperature changes, which will result in the variations of the soil permittivity, not only the real part, but also the imaginary part. When the soil changes from below 0° to above 0°, the soil permittivity suddenly changes, both the real part and the imaginary part increase suddenly. When soil is very dry (vms = 0.1), the variation is not very large, while as the soil moisture increases, the differences increase, i.e., larger soil moisture corresponds to larger variations of permittivity. This phenomenon is just because of the larger permittivity of liquid water than solid ice. Figure 2 also demonstrates that as the soil changes from below 0 °C to above 0 °C, soil moisture will affect the permittivity, since it is a function of soil moisture. Therefore, when we take the soil freeze/thaw process detection into consideration, we need to get rid of the soil moisture influence.
The soil permittivity is a function of frequency. The frequency bands for the Radio Navigation Satellite System (RNSS) cover from lower L band to upper L band (http:// www.navip edia.net/index .php/GNSS_signa l). However, they all work at L band, and frequency effects on the permittivity can be ignored as for different L band GNSS systems. Model inputs of the soil permittivity for the simulation parts are presented in Table 1.
Two representative soil temperatures are selected for comparisons, one is − 1 °C, another one is 1 °C. While the real part of soil permittivity is from 8.7 to 21.1 and the imaginary part has changed from 1.03 to 3.54. These variations of soil permittivity correspond to the changes of bistatic scattering cross-sections. When the soil temperature changes from below 1° to above 1°, Fig. 2 shows the variations of bistatic radar cross-section. L1 carrier (f1 = 1.57542 GHz) is employed for our simulation. The incidence angle is 40°, while the scattering azimuth angle varies from 0.1° to 180° and the zenith angle varies from 0.1° to 360°. For different polarizations at different scattering angles (azimuth angle and zenith angle), the scattering properties vary a lot, while the variations of VV, HH and RR pol are more obvious. There is no scattering as the incidence angle is 0°; therefore, there are blank areas in all the subplots of Fig. 2. The minimum and maximum BRCS differences are also presented in Fig. 3; the positions of these values are related to the scattering azimuth and zenith angle.
According to the plot, Table 2 shows the minimum and maximum BRCS differences at different polarizations. The largest variations are occurred at VH pol; the BRCS difference variation for VH pol is 31.23. As for circular polarization, BRCS difference at RR pol is larger than the one of RL pol and the max BRCS at RR pol is 1.84, while the min BRCS difference at RR pol is − 12.28.
The DDM differences, which are simulated using the above-mentioned models, are presented in Fig. 4. The integral time is 0.001 s. The positions (Unit km) and velocities (Unit m/s) of the transmitter (ts) and receiver (rx) are like the followings: When the soil changes from frozen state to thawn state, i.e., the soil temperature is changed from − 1 to 1 °C. Figure 2 shows the DDM differences at linear polarization and circular polarization under the geophysical reference frame. The DDM difference is like horseshoe, the apparent differences can be seen in Fig. 4.
From the simulations, we can see that as the soil change from frozen state to thawn state, there are apparent changes for the DDMs. The maximum and minimum DDM differences at linear and circular pol are presented in Table 3.
As the soil changes from frozen state to thawn state or vice versa, surface roughness and soil moisture are also two important parameters that affect the final GPS DDM powers. We have employed three kinds of soil for comparisons. As for soil 1 and soil 2, the surface roughness (rms height and correlation length) is the same. While the volumetric soil moisture (vms) is different and they   are used for the evaluations of soil moisture effects; while the soil moisture contents of soil 2 and soil 3 are the same, the surface roughness are different and they are used for the evaluations of surface roughness effects. The final permittivities of different soils are also presented in Table 4. Table 4 presents the soil permittivity differences; soils 1 and 2 are used for the comparisons of soil moisture effects, while soil 2 and 3 are provided for the evaluations of soil surface roughness effects. Different soil moistures result in the variations of final power, which will reflect in the variations of DDM. In addition, they are presented in Fig. 3. The variations of final DDMs at linear polarizations and circular polarizations are different.
The maximum and minimum DDMs differences are presented in Table 4. The maximum difference of DDM is at VV pol. While the minimum DDM difference is at RL pol. The largest variations of different polarizations occur at RR pol. Figure 5 presents the soil moisture effects on the final DDMs at linear and circular polarizations. It can be seen the effects of soil moisture on DDMs are different for different polarizations. Table 5 gives the maximum and minimum DDM differences at linear and circular polarizations for different soil moisture. The DDM difference at RL polarization is maximum, while the difference at VH polarization is minimum. As for the exploration of the final soil freeze/thaw process detection, soil moisture effects on DDMs should be considered.
As for the evaluation of bare soil freeze/thaw process detection using GNSS-R, soil moisture is not the only factor that affects received power. The natural surface is not flat and it is random rough. Roughness is one of the important effects that affect the final scattering properties and DDMs. Surface roughness effects are presented in Fig. 6. Bare soil model inputs for different roughness are given in Table 3. As can be seen from the simulations, soil surface roughness will affect DDMs. As the soil Fig. 4 The DDM differences at linear (HH, VV and HV) and circular (RR and RL) polarization, when the bare soil changes from frozen sate and thawn state

Discussion
This paper presents the bare soil freeze/thaw process detection based on the theoretical simulations and analysis using the GNSS-R. One work that needs further test is the model validation using the experimental data. Due the lack of exactly pertinence experiments, we have employed the method of model validation by model. Our developed model is based on the Z-V model and it is first developed to use for the ocean applications. We have modeled the DDM power of the ocean surface, which is provided in the supplementary file. The result is the same with the original model. The developed model is thought to be composed by the GPS scattering power part and the bare soil surface (frozen and thawn) bistatic scattering part. As for the validation, one is to test the conjunctions between the two parts. And as have mentioned, we have used our methods to get the ocean surface scattering and merges into the Z-V model, the final results are the same with the original model. The permittivity models for frozen and thawn soil have been validated and used as common models, and have been used and tested in the past years. Now the bistatic model is also validated using the Monto Carlo Method, while as for the scattering of different polarizations, we have switched the modified Stokes vector to linear polarizations and they have the same   results with the simulations and therefore, the polarizations at any other modified stokes vectors are thought to be the correct. Since the environment effects of the soil freeze/thaw are complex, this model only takes the bare soil into consideration, while works of snow and vegetation cover should be taken into account in the future.

Conclusions
This paper has developed a fully polarizations GNSS-R Delay-Doppler-Map (DDM) simulation for bare soil freeze/thaw process detection. Theoretical models based on the original Z-V model are developed with taking the GPS geometry and the bistatic scattering properties of the object targets into account. Soil freeze/thaw process, soil moisture and surface roughness are all evaluated in the model. The presented theoretical model indicates the possibility of soil freeze/thaw process detection using airborne/spaceborne GNSS-R and the detections will benefit from the combinations of different polarizations. In the future, some experiments and observations are further tested and analyzed.