*3.2. SMC Retrieval*

#### 3.2.1. Choice of Elevation Angle

Given the lack of measured soil moisture data, we selected the soil moisture value from the P041 station estimated with the new GPS L2C carrier SNR as the reference value. The soil moisture value of the PBO site is based on a value per day and is the median value

estimated by all satellites [22]. As soil moisture changes more during the day than that at night and more GPS satellites are visible during the day than at night, to simplify the calculation, generally, only four and eight observation periods are selected at night and day, respectively. As such, we used twelve observation periods every day. To ensure the solution's reliability and avoid excessive differences in multipath errors from affecting the parameter solution, we assumed that the soil moisture remained constant for a short time (e.g., within 5 min). We used the multipath error of the element used as the observation value for 11 multipath errors in total.

To provide a more intuitive view of the number of GPS satellites visible during the day compared to at night, Figure 6 shows the change in elevation of almost all GPS satellites at site P041 on 19 February 2014 (DOY: 2014–050). Figure 6 shows that more GPS satellites are visible during the day than at night, and the duration of a single GPS satellite being visible during l day is generally about 2.5–8 h, and the period at a low satellite elevation angle (<30◦) is shorter, being maintained at about 0.5–2 h.

**Figure 6.** Changes in the elevation angles of almost all GPS satellites at P041.

Figure 7 shows the linear combination of GPS satellite G22 (DOY: 2014-050) observations at the test site. The multipath error varied with the elevation angle of the satellite. The black curve represents the change trend in MP2 multipath error and L4\_IF multipath error over time, and the red parabola represents the change trend in satellite elevation angle with time. Figure 8 shows that the satellite elevation angle occurs between 13:30 and 21:30. The time gradually increases to the maximum value and then slowly decreases until the satellite disappears. The changes in the MP2 multipath error and L4\_IF multipath error are closely related to the satellite elevation angle. When the satellite elevation angle is low, the MP2 value oscillation amplitude can reach several meters, and the L4\_IF value oscillation amplitude is relatively large. When it is high, the oscillation amplitudes of the MP2 value and L4\_IF value decrease accordingly. For increased accuracy, we selected the data at the lower elevation angle as the experimental data, as shown in Figure 8.

Figure 8 shows the multipath error of the linear combination observations of the GPS satellite G22 (DOY: 2014-050) at a satellite elevation of 5–25◦ for the test site GPS satellite. Figure 8 shows that as the satellite elevation angle increased, the multipath error showed a downward trend, which indirectly indicated that the carrier pseudorange and carrier phase observations of low-elevation satellites are more susceptible to the influence of multipath

errors; it provides a theoretical basis for the selection of GNSS-IR satellite altitude angle. The satellite elevation angles in this experiment were all limited to 5–25◦.

**Figure 7.** The multipath error varies with satellite elevation angle: (**a**) dual-frequency pseudorange multipath error; (**b**) dual-frequency carrier phase linear combination of multipath errors.

**Figure 8.** The multipath error of low elevation angle varies with the sine value of satellite elevation angle: (**a**) the dual-frequency pseudo multipath error MP2; (**b**) the dual-frequency carrier phase linear combination multipath error L4\_IF for removing ionospheric delay.

#### 3.2.2. Choice of Azimuth

Satellite signal reflection point trajectory is a function of satellite elevation, azimuth, and antenna height. The reflection point trajectory reflects the position of the satellite relative to the receiver at a certain moment, from which the spatial information of the reflected soil can be obtained. Therefore, the advantages of the satellite reflection point trajectory map can be used, and the satellite azimuth angle can be considered to ensure as much consistency of the multiperiod multipath environment as possible. We selected the area with more visible satellites as the azimuth angle range at the experimental site.

Given the spatial differences in soil moisture and surface environment in the study area, the soil moisture at different locations was not entirely consistent, which led to the difference between the retrieval results based on the observation values at different reflection locations. In addition, not all DFP observations and L4\_IF observations of the satellite elevation state contain the physical information of the surface reflector. Therefore, we needed to select observation values based on the study area environment, satellite elevation angle, and satellite azimuth angle to improve the accuracy of GNSS-IR soil moisture retrieval. As shown in Figure 9, as the reflection area of the soil was roughly the same, and on the premise of ensuring the number of satellites in each period, we selected an azimuth angle range of 30–330◦. The satellites all had continuous observation values when their elevation angle was 5–25◦, providing relatively sufficient observation data.

**Figure 9.** The projected trajectory of the P041 (DOY: 2014-065) satellite at an elevation of 5–25◦.

3.2.3. Selection of Effective Satellites

The Fresnel reflection region of a GNSS signal is a set of ellipses related to the elevation, azimuth, and antenna height of a satellite. The first Fresnel reflection region contributes the most to the reflected signal, and the reflection medium change in the reflection zone strongly affects the relevant physical characteristics of the reflected signal.

The semimajor and semiminor axes of the ellipse of the first Fresnel reflection region of the ground-based GNSS-R can be determined by [28]:

$$S\_{\mathbf{x}} = h \cdot \cot \theta\_{\mathbf{\cdot}} \tag{30}$$

$$S\_y = 0,\tag{31}$$

$$a = \frac{\sqrt{\lambda h \sin \theta}}{\sin^2 \theta},\tag{32}$$

$$b = \frac{\sqrt{\lambda h \sin \theta}}{\sin \theta},\tag{33}$$

where *Sx* is the projection position of the reflection point on the ground connecting the satellite and the receiver; *Sy* is the position of the receiver projection point; *a* and *b* are semimajor and semiminor axes of the first Fresnel reflection region; *h* is the distance from the antenna phase center to the reflection surface; *λ* is the carrier wavelength; *θ* is the elevation angle of the GNSS satellite.

Figure 10 shows the first Fresnel reflection region ellipse group of the GPS *L*<sup>2</sup> carrier when the satellite azimuth was 0◦, the satellite elevation was 5–25◦, and the receiver antenna height was 1.9 m. The horizontal axis represents the distance from the receiver antenna to the direction of the satellite, and the vertical axis represents the distance perpendicular to the direction of the receiver antenna and satellite. When the satellite elevation angle increases, the Fresnel reflection region shrinks, and the reflection center is closer to the GNSS receiver. The area of the ellipse can be obtained from the semimajor and semiminor axes of the ellipse, that is, the area of the first Fresnel reflection region. Using the first Fresnel reflection region can avoid the influence of other multipath sources, such as tall buildings, vegetation, and water; maximize the unity of the reflection medium in the reflection zone; and be used for the design of the experimental plan. Therefore, the Fresnel reflection region provides a certain basis for determining the location of soil moisture sample collection points and selecting effective satellites. The change in the reflection medium in the Fresnel reflection region considerably affects the signal received by the receiver.

**Figure 10.** First Fresnel reflection region for GPS *L*<sup>2</sup> frequency at different elevation angles.

Figure 11 shows the first Fresnel reflection region map of station P041 (DOY: 2014-065). When the satellite elevation angle was 10◦, the reflection point was up to 21 m away from the receiver. If slopes or other interference sources are near the measurement area, the first Fresnel reflection region should be introduced to reduce the excessive soil moisture retrieval error caused by different soil reflection media. The lower the satellite elevation angle, the larger the region of the first Fresnel reflection. When the elevation angle range is limited, the azimuth angle of each GPS satellite does not change much. Therefore, the calculation of the first Fresnel reflection region can ensure the unity of the medium in the reflection region, meaning that the homogeneity of the soil can be more accurately assessed and can provide a certain basis for other factors such as experimental planning and satellite selection.

The maximum power spectral density can be used to characterize the quality of a multipath error signal to a certain extent. Spectrum analysis using the Lomb–Scargle method can be used to judge satellite signal quality [22,29]. A satellite track should have a relatively stable singular dominant frequency for the same reflecting surface. In general, the dominant frequency's power spectral density (PSD) should be at least twice as high as the power of the noise or the second most powerful frequency in the periodogram [22,30].

**Figure 11.** First Fresnel reflection region map of station P041 (DOY: 2014-065): GPS satellite distribution in the first Fresnel reflection region of station P041 when satellite elevation angle was (**a**) 10◦; (**b**) 25◦.

Figure 12a,b show that the power that meets the dominant frequency should be at least twice the power of the secondary main peak frequency, and that the quality of satellite data is higher. We used the PSD of dual-frequency pseudorange and dual-frequency carrier phase combined multipath error as the main auxiliary tool for satellite selection. Figure 12c,d lack any dominant frequency, which could have introduced substantial errors in the subsequent parameter estimation. Therefore, we also discarded such satellites when selecting satellites.

Based on the above analysis, after limiting the satellite elevation angle to about 5◦–25◦ and after further considering the dual-frequency multipath error for satellite selection, we used twelve time periods every day: eight evenly distributed in the daytime and four and evenly distributed in the nighttime. Although the multipath error duration of a single GPS satellite in a day is generally about 2–8 h, the corresponding duration of the low satellite elevation angle multipath error is shorter, basically maintained at about 1 h. As a result, even fewer GPS satellites are available at low elevations. Therefore, we selected only one GPS satellite for each observation period. Combining satellite elevation angle, satellite reflection trajectory, satellite azimuth, first Fresnel reflection area, station environment, spectrum analysis, etc., Table 2 shows the results of the twelve time periods on 6 March 2014.

**Figure 12.** Multisatellite dual-frequency combined multipath error Lomb–Scargle periodogram. Spectrum analysis diagrams of the (**a**,**c**) dual-frequency carrier phase combined multipath error (L4\_IF;) (**b**,**d**) dual-frequency pseudorange multipath error (MP2).



#### *3.3. GNSS-IR SMC Retrieval Results*

To obtain the parameters that characterize the change trend in soil moisture, we used the multipath error as the system input, and the amplitude attenuation factor and phase delay were the parameters to be estimated. We ignored the impact of the satellite elevation, soil moisture changes in the short term, and the amplitude attenuation factor and phase. The delay was a constant in the short term, and we obtained the error equation by linearizing Equation (5). The initial value of the phase delay was determined using Equation (2), where the satellite elevation angle is the satellite elevation angle corresponding to the first soil moisture sample collection time, considering that the soil surface reflectivity is between 0.3 and 0.8 [21]. In this study, we focused on the amplitude and phase of multipath error, so the actual value of *κ* is not required. Taking an *κ* initial value *κ*<sup>0</sup> = 0.3, we ignored the complex effect of antenna gain mode.

Considering the small changes in the sine values of antenna height and satellite height and our focus on the changing trend in the phase delay, we used the same initial phase delay value for all periods of data processing. For the solution of the parameters, we adopted the indirect adjustment method based on least squares, and the data processing followed these principles:


Combining the results of daily satellite selection, for the selected single satellite, we calculated the multipath error according to the dual-frequency signal of 11 epochs before and after the corresponding time. Combining the above principles and methods, we obtained the delay phase corresponding to the time of soil moisture acquisition through adjustment. Therefore, one phase delay could be calculated for each observation period, and 12 phase delays could be calculated for 12 observation periods in a day. By averaging the 12 results, we obtained the average phase delay, which is the daily phase delay. On a daily time scale, the L4\_IF method and the DFP method were calculated separately to obtain the corresponding daily phase delay, and Time trends of soil moisture and delayed phase at P041 and MFLE stations were plotted (Figures 13a,b and 14a,b).

(**b**) Correlation between DFP multipath model and SMC.

**Figure 13.** P041 station; correlation between Delay Phase and SMC.

To further test the L4\_IF method and the DFP method, 50 d of the P041 and MFLE SMC, respectively, and the estimated phase delay were employed to build an empirical model, and the other 28 d of data were used to verify the model. The modelling method employed was unary linear. In the relationship between soil moisture and phase delay, in this study, we used the phase delay as the independent variable (x) and the soil moisture as the dependent variable (y) to perform univariate linear regression modeling. We obtained the scatter plot and regression shown in Figure 15.

**Figure 14.** MFLE station; correlation between Delay Phase and SMC.

To better reflect the soil moisture retrieval performance of the two methods of P041 and MFLE stations, we performed precision statistical analysis on each method separately, and the statistical results are shown in Table 3.

Soil moisture is often affected by vegetation cover, soil temperature, air humidity and other factors. In addition, due to the difference between the actual reflection conditions of GNSS satellite signals and the hypothetical ideal reflection conditions, the relationship between multi-path induced phase delay and soil moisture becomes complicated, which makes the ordinary linear model not necessarily able to describe the trend of soil moisture. Compared with other analytical models, the neural network model is less sensitive to the influence of soil flatness, surface vegetation and soil surface roughness. Therefore, in order to improve the inversion accuracy of GNSS-IR soil moisture, BPNN and RBFNN were used to construct soil moisture prediction models, and the prediction results were compared with those of ULR model, and the accuracy of each model was evaluated.

**Figure 15.** For the P041 station; scatter diagram of the SMC and regression equation based on (**a**) the L4\_IF method; (**b**) the DFP multipath model. For the MFLE station; scatter diagram of the SMC and regression equation based on (**c**) the L4\_IF method; (**d**) the DFP multipath model.


**Table 3.** Statistical table for accuracy comparison between L4\_IF and DFP methods.

Note: The STD in Table 3 is the error standard deviation; the MAE in Table 3 is mean absolute error; the RMSE in Table 3 is root-mean squared error.

The phase delay and measured soil moisture of P041 and MFLE stations calculated based on L4\_IF method and DFP method were used as data sources, and the training set and test set were divided according to the ratio of 7:3. That is, 50 d data were used as training set data, and 28 d data were used as test set data. BPNN and RBFNN are used to model. Figures 16 and 17 show the comparison between the predicted values of the three models of the two methods corresponding to the P041 station and the MFLE station and the measured values of soil moisture.

**Figure 16.** For P041 station; comparison between the predicted values of three models and measured values of soil moisture (**a**) the L4\_IF method; (**b**) the DFP multipath model.

In order to better reflect the prediction results of three prediction models for soil moisture and the comparison results between different models, the accuracy of L4\_IF and DFP methods corresponding to the three models of P041 and MFLE stations were statistically analyzed. The statistical results are shown in Table 4.

**Figure 17.** For the MFLE station, comparison between predicted values of three models and measured values of soil moisture (**a**) the L4\_IF method; (**b**) the DFP multipath model. Note: The P041\_SMC in Figure 16a,b are SMC value of the P041 station; the MFLE\_SMC in Figure 17a,b are SMC values of the MFLE station; the ULR\_SMC in Figures 16 and 17 are unitary linearity regression model predictive value; the BPNN\_SMC in Figures 16 and 17 are back propagation neural network model predictive value; the RBFNN\_SMC in Figures 16 and 17 are radial basis function neural network model predictive value.



Note: The STD in Table 4 is error standard deviation; the MAE in Table 4 is mean absolute error; the ULR in Table 4 is unitary linearity regression; the BPNN in Table 4 is back propagation neural network; the RBFNN in Table 4 is radial basis function neural network.

#### **4. Discussion**

From the above experimental results can be drawn, Figures 13a,b and 14a,b show multiple apparent peaks in the phase delay and soil moisture, and the overall change trend has a strong consistency. For the P041 station, the correlation coefficients between the soil moisture and phase delay calculated based on the L4\_IF multipath error and the DFP multipath error are 0.97 and 0.91, respectively. For the MFLE station, the correlation coefficients between the soil moisture and phase delay calculated based on the L4\_IF multipath error and the DFP multipath error are 0.93 and 0.86, respectively, which are statistically significant. Figure 15a–d shows that the soil moisture values obtained using the L4\_IF and DFP methods both fluctuate around the linear regression equation built by each, and the deviation is slight. The L4\_IF and DFP methods have a statistical significance level of 0.01 and 0.05, respectively, regarding soil moisture retrieval. The probability values of the two methods are close to zero (*P* ≈ 0), indicating that the correlation between the estimated phase delay and soil moisture at the P041 and MFLE stations is significant. Figures 16 and 17 show that the fitting degree between the predicted and measured values of the BPNN model and RBFNN model is better than that of the ULR model, and the error is smaller.

Table 3 shows that the standard deviation and absolute average error of the L4\_IF and DFP methods are the same, but the root mean square error of the L4\_IF method is smaller than that of the DFP method. Therefore, the L4\_IF method is generally more accurate than the DFP method. The observation precision of the carrier phase is much higher than that of the pseudorange. To better compare the L4\_IF and DFP methods, we used the error propagation law to compare the error level of L4\_IF and DFP. Assuming that the observation error of the dual-frequency carrier phase has the same standard deviation of 1 mm (*σ*<sup>0</sup> = *σ*<sup>1</sup> = *σ*<sup>2</sup> = 1mm), the standard deviation of the GPS C/A code measurement pseudorange observation value is 2.93 m. Therefore, the error standard deviation of the L4\_IF method is *<sup>σ</sup>L*<sup>4</sup> <sup>=</sup> <sup>√</sup>2*σ*<sup>0</sup> <sup>=</sup> 1.41 mm, and that of the DFP method is *σMP*<sup>2</sup> = 2.93 m *σL*<sup>4</sup> < *σMP*<sup>2</sup> , also means that the quality of the multipath error of the L4\_IF method is higher than that of the DFP method, so the L4\_IF method obtained a higher correlation coefficient and root mean square error.

Table 4 shows that for the L4\_IF and DFP methods of P041 and MFLE stations, the prediction results of the BPNN model are better than the RBFNN model and ULR model, and the RBFNN model is slightly better than the ULR model. The reasons may be as follows. Firstly, soil moisture is affected by vegetation cover, soil temperature, air humidity and other factors, and multipath reflection is not an imaginary mirror reflection, which makes it difficult to characterize other nonlinear characteristics between delay phase and soil moisture by using simple linear or exponential models. Secondly, the BPNN algorithm can correct the dielectric constant difference caused by weakening soil roughness, and has a certain inhibitory effect on surface fluctuation, soil roughness and vegetation.

#### **5. Conclusions**

Based on the analysis of the multipath error generation mechanism and the calculation model of GNSS measurement, we constructed a soil moisture retrieval method based on multisatellite dual-frequency combined multipath error. We proposed a method of estimating soil moisture using two combined dual-frequency signals. We used the L4\_IF and DFP methods, which are independent of geometric distance, to estimate soil moisture. We verified the proposed method using the archived data set and SMC data of the P041 and MFLE stations. We drew the following conclusions through experimental analysis:

(1) The delay phases obtained by the multipath error solution and the soil moisture are strongly correlated. For the P041 station, the R values of the L4\_IF and DFP methods are 0.97 and 0.91, respectively. For the MFLE station, the R values of the L4\_IF and DFP methods are 0.93 and 0.86, respectively. Because the observation error of the L4 linear combination is low and the change in the ionospheric delay in the short term is small, we used the high-order fitting to further weaken the influence of the ionospheric

delay. The L4\_IF method has higher R and soil moisture estimation accuracy than the DFP method. When the BPNN model, RBFNN model and ULR model are used to predict soil moisture, the results show that the prediction results of the BPNN model are better than the RBFNN model and ULR model, and the RBFNN model is slightly better than the ULR model. The results show that BPNN can improve the inversion accuracy of GNSS-IR soil moisture.


If SNR does not exist in the original file of GNSS and the number of GNSS satellites with triple-frequency signals is still limited, the L4\_IF method and DFP method proposed in this paper can be used as alternative methods for monitoring through GNSS-IR so as to enrich the data source of GNSS-IR soil moisture inversion and improve the ability of GNSS to serve environmental monitoring. However, it should be noted that not all satellite observations are suitable for estimating soil moisture, and different satellite selections will lead to different results. Therefore, the effective selection of satellites is still a challenge. Because this method does not need to diagnose the signal frequency and only needs less epoch multipath error to calculate the delay phase, it is easier to achieve higher time resolution of GNSS-IR soil moisture inversion. Soil moisture inversion based on a multipath error enriches GNSS-IR data sources and enhances the reliability of GNSS-IR.

**Author Contributions:** Conceptualization, Y.W., J.T. and P.L.; methodology, S.N.; formal analysis, J.X., N.L., M.W. and P.L.; data curation, J.X., D.H. and J.S.; writing—original draft preparation, S.N.; writing—review and editing, Y.W. and J.T.; supervision, Y.W. and J.T. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Anhui Educational Commission Key Project, grant numbers: KJ2020A00706, KJ2021A1077, KJ2021A1078, KJ2021A1084; the Chuzhou Science and Technology Guidance Program Project, grant number: 2021ZD009; the Anhui Educational Commission General Project, grant number: KJ2021B08, KJ2021B07.

**Data Availability Statement:** The GNSS site data were provided under the PBO Observation Program of the United States, and the measured soil moisture data were obtained from https: //cires1.colorado.edu/portal (accessed on 10 September 2021) and https://data.unavco.org/archive/ gnss/products (accessed on 10 September 2021).

**Acknowledgments:** We thank the PBO Observation Program of the United States for providing the GNSS data, and the University of Colorado for providing data for soil moisture comparison analysis.

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

#### **References**

