*2.1. GNSS-IR SMC Retrieval Principle*

2.1.1. GNSS Multipath Error Principle

In the actual measurement, the signal received by the GNSS receiver antenna is not only the signal directly from the satellite but also the signal reflected by the surface. The direct and the reflected signals entering the receiver antenna interfere with each other, which causes the observation value to deviate from the actual value and produce the multipath error. According to the roughness of the reflective surface of the ground, the two main types of reflection are diffuse reflection and specular reflection [22,24]. For convenience, we assume that multipath error is only caused by specular reflections. Figure 1 shows that

direct signal and reflected signal generate corresponding interference effects at the receiver to form a composite interference signal.

**Figure 1.** Schematic diagram of GNSS-IR SMC retrieval. After the satellite sends the signal, the right-handed circular polarized (RHCP) antenna receives the direct signal and the surface reflected signal, producing an interference effect at the receiver. A is the GNSS receiver antenna; B is the reflection point position of GNSS satellite signals passing through the ground; O is the footing point of the vertical line between the over-reflection point and the GNSS satellite direct signal; *θ* is the elevation angle of the satellite; *S*<sup>1</sup> is the direct satellite signal received by the receiver antenna; *S*<sup>2</sup> is the reflected signal reflected by the object's surface around the receiver antenna and enters the antenna; *H* is the vertical height from the antenna phase center to the ground.

The elevation of the GPS satellite (*θ*) determines the path delay Δ*S*, which is the extra distance travelled by the reflected signal compared to that of the direct signal.

$$
\Delta S = 2H \sin \theta\_\prime \tag{1}
$$

By substituting carrier wavelength *λ* into Equation (1), the phase delay *δϕ*(*t*) corresponding to the path delay can be expressed as:

$$
\delta\varphi(t) = 2\pi \frac{\Delta S}{\lambda} = \frac{4\pi H}{\lambda} \sin \theta(t),
\tag{2}
$$

where *λ* is the carrier wavelength; *t* is the observation epoch. The phase delay is related to the antenna height, carrier wavelength and satellite elevation angle; that is, for a given antenna height and carrier signal, the phase delay is a function of the elevation angle. Equations (1) and (2) only consider the geometric delay and ignore the phase contribution of the Fresnel reflection and antenna radiation.

The amplitude *AC* and phase difference *β*(*t*) of the reflected signal compared to the direct signal can be expressed as:

$$\begin{cases} A\_C = \sqrt{A\_d^2 + A\_m^2 + 2A\_d \cdot A\_m \cdot \cos \delta \rho(t)} \\ \qquad \beta(t) = \tan^{-1} \left[ \frac{\kappa \cdot \sin \delta \rho(t)}{1 + \kappa \cdot \cos \delta \rho(t)} \right] \end{cases} \tag{3}$$

where *Am* = *κ* · *Ad*, *κ* is the amplitude attenuation factor (AAF); *SC* is the composite signal; *AC* is the amplitude of *SC*. The factor *β*(*t*) is the composite excess phase with respect to the direct phase, which can be approximately expressed as [21]:

$$
\beta(t) = \left(\frac{A\_m}{A\_d}\right) \cdot \sin \delta \varphi(t) = \kappa \cdot \sin \delta \varphi(t),
\tag{4}
$$

When only specular reflection is considered, the composite signal formed by the superposition of the direct signal and the reflected signal *SC* can be expressed as [21]:

$$S\_{\mathbb{C}} = S\_L + S\_M = A\_d \cos(\omega\_0 t) + A\_m \cos[\omega\_0 t + \delta \varphi(t)] = A\_{\mathbb{C}} \cos[\omega\_0 t + \beta t],\tag{5}$$

where *SL* and *SM* are the direct and reflected signals, respectively; *Ad* and *Am* are the amplitudes of the direct and reflected signals, reflectively; *ω*<sup>0</sup> is the angular frequency of the signal.

#### 2.1.2. Calculation of Multipath Error of Linear Combination of Observations

If only the influence of atmospheric delay error is considered, the GNSS code measurement pseudorange and carrier phase observation equation can be approximately described as [25]:

$$P\_i = \rho + \varepsilon (\delta t\_R - \delta t\_S) + T + I\_i + E + MP\_i + \varepsilon\_{pi} \tag{6}$$

$$
\lambda\_i \circ \rho\_i = \rho + \varepsilon (\delta t\_R - \delta t\_S) + T - I\_i + E + M\rho\_i + \varepsilon\_{\varphi\_i} - \lambda\_i N\_{i\prime} \tag{7}
$$

where *i* is the carrier number; *ϕ* is the carrier phase observation value; *ρ* is the geometric distance between the receiver and the satellite; *c* is the propagation velocity of electromagnetic waves in a vacuum; *δtR* and *δtS* are the clock difference between the receivers and satellite clock, respectively; *N* is the ambiguity of the whole week; *λ* is the carrier wavelength; *MP* and *M<sup>ϕ</sup>* are the pseudorange and carrier multipath error, respectively; *T* and *I* are the tropospheric and ionospheric delay errors, respectively; *ε* represents other unmodeled errors.

#### 2.1.3. Error Calculation of Dual-Frequency Pseudorange Multipath

As shown in Equations (6) and (7), when *i* is 1 and 2, Equations (6) and (7) can be differentiated, respectively. Considering |*MPi*||*Mϕi*|, *<sup>ε</sup> pi*||*εϕ<sup>i</sup>* is as shown in Equations (8) and (9).

$$P\_1 - \lambda\_1 \varphi\_1 = 2I\_1 + MP\_1 + \varepsilon\_{\varphi1} + \lambda\_1 N\_1. \tag{8}$$

$$P\_2 - \lambda\_2 \varphi\_2 = 2I\_2 + MP\_2 + \varepsilon\_{\varphi2} + \lambda\_2 N\_{2\prime} \tag{9}$$

The relationship between carrier frequency *f* and ionospheric delay *I* is:

$$I\_2 = \frac{f\_1^2}{f\_2^2} I\_{1\prime} \tag{10}$$

By subtracting the *L*<sup>1</sup> and *L*<sup>2</sup> carrier phase observation equations and taking into account Equation (10), we obtain:

$$
\lambda\_1 \varphi\_1 - \lambda\_2 \varphi\_2 = I\_2 - I\_1 + \lambda\_1 N\_1 - \lambda\_2 N\_2 = \left(\frac{f\_1^2}{f\_2^2} - 1\right) I\_1 + \lambda\_2 N\_2 - \lambda\_1 N\_1,\tag{11}
$$

Therefore, as shown in Equation (12):

$$I\_1 = \frac{(\lambda\_1 \varrho\_1 - \lambda\_2 \varrho\_2 + \lambda\_1 N\_1 - \lambda\_2 N\_2) \times f\_2^2}{f\_1^2 - f\_2^2},\tag{12}$$

Therefore, the pseudorange multipath error equation at *L*<sup>1</sup> band is:

$$\begin{split}MP\_1 &= P\_1 - \lambda\_1 \varrho\_1 - 2I\_1 - \varepsilon\_{P1} - \lambda\_1 N\_1\\ &= P\_1 - \frac{f\_1^2 + f\_2^2}{f\_1^2 - f\_2^2} \lambda\_1 \varrho\_1 + \frac{2f\_2^2}{f\_1^2 - f\_2^2} \lambda\_2 \varrho\_2 + K(N\_{1\prime}, N\_{2\prime}, \varepsilon\_{P1}),\end{split} \tag{13}$$

Similarly, the pseudorange multipath error equation in *L*<sup>2</sup> band is:

$$\text{MP2} = \text{P2} - \frac{2f\_1^2}{f\_1^2 - f\_2^2} \lambda\_1 \varphi\_1 + \frac{f\_1^2 + f\_2^2}{f\_1^2 - f\_2^2} \lambda\_2 \varphi\_2 + K(\mathcal{N}\_1, \mathcal{N}\_2, \varepsilon\_{P2}), \tag{14}$$

where *MP*<sup>1</sup> and *MP*<sup>2</sup> are the multipath errors of the *L*<sup>1</sup> and *L*<sup>2</sup> carriers, respectively; *P*<sup>1</sup> and *P*<sup>2</sup> are the pseudorange observations of the *L*<sup>1</sup> and *L*<sup>2</sup> carriers, respectively; *f*<sup>1</sup> and *f*<sup>2</sup> are the frequencies of the *L*<sup>1</sup> and *L*<sup>2</sup> carriers, respectively; *λ*<sup>1</sup> and *λ*<sup>2</sup> are the *L*<sup>1</sup> and *L*<sup>2</sup> carrier wavelengths, respectively; *ϕ*<sup>1</sup> and *ϕ*<sup>2</sup> are the observed values of *L*<sup>1</sup> and *L*<sup>2</sup> carrier phases, respectively; *N*<sup>1</sup> and *N*<sup>2</sup> are the ambiguities of *L*<sup>1</sup> and *L*<sup>2</sup> carriers, respectively; *εP*<sup>1</sup> and *εP*<sup>2</sup> are other unmodeled errors.

In Equations (13) and (14), *K*(*N*1, *N*2) is the integer ambiguity combination, which is usually a constant and does not affect the overall trend in the pseudorange multipath error. If no cycle slip occurs, it can be omitted. Therefore, the main factors affecting the quality of pseudorange multipath error are the noise in pseudorange observations and the accuracy of carrier phase observations and code pseudorange observations. If cycle slip occurs in carrier phase observations, it affects the value of pseudorange multipath error. Therefore, for the pseudo-range multipath error in each period, cycle slip detection and repair are required. We selected the multipath error *MP*<sup>2</sup> of the *L*<sup>2</sup> carrier as the research object.

The dual-frequency pseudorange multipath can be calculated by the linear combination of the pseudorange observations and the carrier phase observations, which can be expressed as [26]:

$$MP\_1 = P\_1 - \frac{f\_1^2 + f\_2^2}{f\_1^2 - f\_2^2} \lambda\_1 \varphi\_1 + \frac{2f\_2^2}{f\_1^2 - f\_2^2} \lambda\_2 \varphi\_2,\tag{15}$$

$$MP\_2 = P\_2 - \frac{2f\_1^2}{f\_1^2 - f\_2^2} \lambda\_1 \varphi\_1 + \frac{f\_1^2 + f\_2^2}{f\_1^2 - f\_2^2} \lambda\_2 \varphi\_2. \tag{16}$$

2.1.4. Multipath Error Calculation of Dual-Frequency Carrier Phase Linear Combination

For carrier phase observables, isolating the multipath from several other effects, such as the satellite antenna pseudorange and ionospheric and tropospheric delay, is difficult. However, when considering the multipath, the carrier phase contains two quantities: the phase of the direct signal and the composite excess phase with respect to the direct phase. Thus, the GNSS phase observation for carriers *L*<sup>1</sup> and *L*<sup>2</sup> can be expressed as [27]:

$$\begin{cases} L\_1 = \lambda\_1 \psi\_1 + \lambda\_1 \beta\_1(t) \\ L\_2 = \lambda\_2 \psi\_2 + \lambda\_2 \beta\_2(t) \end{cases} \tag{17}$$

where *ψ*<sup>1</sup> and *ψ*<sup>2</sup> represent the phase of the direct signal; *β*1(*t*) and *β*2(*t*) are the composite excess phase with respect to the direct phase.

In GNSS positioning, the carrier phase of the direct signal can be expressed as:

$$
\lambda\_i \psi\_i = \rho + I\_i + T + \Delta\_\prime \tag{18}
$$

where *i* = 1 or 2, which is the two signals with different frequencies. The geometric distance between the GNSS satellite and antenna is represented by *ρ*, *Ii* represents ionospheric delays, *T* represents tropospheric delays, and Δ accounts for all other carrier-independent effects. To isolate the multipath, the dual-frequency carrier-phase combination *L*<sup>4</sup> is the *L*<sup>1</sup> − *L*<sup>2</sup> combination of carrier phase measurement, and its mathematical expression is [9,27]:

$$L\_4 = L\_1 - L\_2 = I\_1 - I\_2 + \lambda\_1 \beta\_1(t) - \lambda\_2 \beta\_2(t),\tag{19}$$

Equation (19) eliminates the satellite clock error, receiver clock error, and tropospheric delay, as well as the geometric distance between the satellite and the receiver. That is the observed value of *L*<sup>4</sup> is affected by both ionospheric delay and phase error. Therefore, when analyzing the multipath error, a suitable method must be adopted to weaken the influence of the ionospheric delay on the multipath error. The high-order polynomial is used to fit the *L*<sup>4</sup> multipath error, and the dual-frequency linear combination multipath error that is not affected by ionospheric delay is obtained, which is the L4-free (L4\_IF) observation value.

#### *2.2. Establishment of Model Error Equation*

According to Equations (1) and (2), the phase delay and path delay are functions of the elevation angle and reflector height. The reflector height changes with the SMC, which contributes to the change in the phase delay and path delay [12]. demonstrated that the agreement between *H* and the SMC was not as strong as that between *δϕ*(*t*) and the SMC. Therefore, the phase delay can be used to characterize the change in the SMC. To accomplish this, Equation (5) must be linearized first. Thus, the error equations can be expressed as:

$$\beta(t) + V\_{\beta(t)} = \kappa \rho \sin \delta \rho \wp + \sin \delta \wp \wp V\_{\mathbb{K}} + \kappa \wp \cos \delta \wp \wp V\_{\delta \wp \var{ \searrow}} \tag{20}$$

where *β*(*t*) can be easily calculated through Equations (16), (17) and (19). *Vβ*(*t*) is the residuals corresponding to the DFP multipath error and L4\_IF multipath error. *κ*<sup>0</sup> denotes the initial value of the AAF. The initial value of the phase delay, i.e., *δϕ*0, can be calculated by incorporating the wavelength (*λ*), antenna height (*H*), and elevation angle (*θ*) into Equations (1) and (2), respectively. *V<sup>κ</sup>* and *Vδϕ* are two unknown correction parameters corresponding to *κ*<sup>0</sup> and *δϕ*<sup>0</sup> that need to be solved. To ensure the reliability of the solution and to avoid the influence of the excessive difference in multipath errors on the correction parameter solving, we assumed that the SMC remained the same in the short term (e.g., within 5 min) and employed a total of 11 multipath errors to solve the correction parameters through unweighted least squares adjustment. It is worth mentioning that the number of multipath errors is not limited to 11, but should be greater than or equal to the number of the correction parameters to be solved.

#### *2.3. Solving the Phase Delay*

The least-squares adjustment method is used to solve the delay phase. The L4\_IF method is taken as an example to illustrate the solving process of the phase delay. For the 11 multipath errors of one observation session for a single GPS satellite, 11 error equations, such as Equation (20), can be obtained. Then, these error equations can be expressed by the following matrix:

$$\begin{array}{c} V\_{\beta(t)} = \frac{A}{11 \times 22 \times 1} - \frac{l}{11 \times 1'}\\\end{array} \tag{21}$$

where

$$\begin{array}{c} V\_{\not\in\{t\}} = \left[ V\_{\not\in1} \; V\_{\not\in2} \; V\_{\not\in3} \; \cdots \; V\_{\not\in11} \right]^T \; \end{array} \tag{22}$$

$$\begin{array}{c} \stackrel{X}{2} \times 1 = \left[ V\_k \, V\_{\delta \varphi} \right]^T \end{array} \tag{23}$$

$$\begin{array}{llll} A & = \begin{bmatrix} \sin\delta\varrho\_0 & \sin\delta\varrho\_0 & \sin\delta\varrho\_0 & \cdots & \sin\delta\varrho\_0\\ \kappa\_0\cos\delta\varrho\_0 & \kappa\_0\cos\delta\varrho\_0 & \kappa\_0\cos\delta\varrho\_0 & \cdots & \kappa\_0\cos\delta\varrho\_0 \end{bmatrix}^T \\ & & \tag{24}$$

$$\frac{1}{11 \times 1} = \begin{bmatrix} \beta\_1 - \kappa\_0 \sin \delta \varrho\_0 \ \beta\_2 - \kappa\_0 \sin \delta \varrho\_0 \ \beta\_3 - \kappa\_0 \sin \delta \varrho\_0 \ \cdots \ \beta\_{11} - \kappa\_0 \sin \delta \varrho\_0 \ \end{bmatrix}^T,\tag{25}$$

where *β*(*t*) represents the L4\_IF multipath error, *Vβ*(*t*) is the residual for *β*(*t*). The subscripts 1, 2, 3, ... , 11 denote the serial number of the selected 11 multipath errors. Thus, based on unweighted least-squares adjustment (*VTV* = min). *X* can be obtained by the following formula.

$$X = \begin{bmatrix} V\_{\mathbf{x}} \\ V\_{\delta \rho} \end{bmatrix} = \left( A^T A \right)^{-1} \times \left( A^T l \right) = \begin{array}{c} N^{-1} & W \\ \mathbf{2} \times \mathbf{2} \, \mathbf{2} \times \mathbf{1}' \end{array} \tag{26}$$

where *N* = *ATA*;*W* = *ATl*. Because its rank equals one, *N* is rank deficient, which leads to a non-unique solution of *X*. Therefore, to obtain a unique solution, pseudo-inverse adjustment was adopted.

$$X = \begin{bmatrix} V\_{\mathbf{x}} \\ V\_{\delta \varphi} \end{bmatrix} = \begin{array}{c} N^- & W \\ \mathbf{2} \times \mathbf{2} \ \mathbf{2} \times \mathbf{1}' \end{array} \tag{27}$$

where *N*−<sup>1</sup> and *N*<sup>−</sup> is the Kelley inverse and the pseudo-inverse of *N*, respectively. Thus, the adjusted phase delay and AAF can be obtained:

$$
\begin{bmatrix}
\hat{\mathfrak{k}} \\
\hat{\delta\varphi}
\end{bmatrix} = \begin{bmatrix}
\kappa\_0 + V\_\mathbf{x} \\
\delta\varphi\_0 + V\_{\delta\varphi}
\end{bmatrix}'\tag{28}
$$

The diurnal phase delay representing the trend of the soil moisture change, namely the mean of the adjusted phase delay of twelve observation sessions, can be obtained by the following formula:

$$
\overline{\delta\varphi}\_i = \frac{\overline{\delta\hat{\varphi}\_{i,1}} + \overline{\delta\hat{\varphi}\_{i,2}} + \overline{\delta\hat{\varphi}\_{i,3}} + \dots + \overline{\delta\hat{\varphi}\_{i,12}}}{12},
\tag{29}
$$

where the subscript *i* equals 1, 2, 3, ... , 78, and represents the *i*th day; the numbers 1, 2, 3, . . . , 12 represent the serial number of observation sessions on a single day.

#### *2.4. Data Sources*

The experimental data are from the USA. Plate Boundary Observatory (PBO). We selected the GPS observation data of the P041 station (Figure 2) and the SMC data collected by the P041 and MFLE stations (https://cires1.colorado.edu/portal/) (accessed on 10 September 2021). The days of the year (DOYs) of selected observation data are between 45 and 131 of 2014 (Note: 57, 67–69, 71, 82, 94, 104, 105 days of soil moisture data are missing and have been eliminated).

The PBO SMC consisted of the site averaged SMC on a daily timescale, and the median SMC value of all satellite tracks for each day was used [22]. (https://data.unavco.org/ archive/gnss/products/) (accessed on 15 September 2021). Located in Colorado, USA, and positioned at an altitude of 1728.8 m, the P041 station (39.9495◦N, 105.1943◦W). The MFLE station is about 210 m away from the P041 station, its altitude is 1727.3 m, and its position is longitude and latitude (39.9476◦N, 105.1944◦W). The soil moisture data of the MFLE station is used to further verify this method's reliability. In view of this, this paper takes

P041 station as the research station, MFLE station as the verification station, and MFLE station is very close to P041 station, as shown in Figure 2. The surrounding environment of P041 and MFLE stations are shown in Figure 3 (https://www.unavco.org) (accessed on 25 September 2021). The GNSS receiver and related parameters of the P041 station are shown in Table 1.

**Figure 2.** GNSS test site locations. P041 and MFLE stations, Boulder, CO, USA.

**Figure 3.** Station conditions: (**a**) Environment around P041 station; (**b**) environment around MFLE station. Note: GNSS test site locations: P041 and MFLE, Boulder, CO, USA (https://www.unavco.org) (accessed on 25 September 2021).


**Table 1.** Receiver and its related parameters at the P041 site.

The P041 and MFLE stations can be seen in Figure 3a,b, respectively. They are located in a flat and open area, and the surrounding vegetation is scarce, and the land cover consists of exposed soils and short grasses. This paper selects the experimental data for the late winter and early spring. Therefore, the surface reflection signal is less affected by vegetation attenuation.

Figure 4 shows the observed SMC and precipitation data series corresponding to the time (DOYs: 45-131,2014) of selecting GNSS observation data from P041 and MFLE stations, which are presented in a line graph and a histogram, respectively. As demonstrated in Figure 4, six significant indigenous precipitation events occurred during the experiment, with a maximum precipitation of 24.8 mm, mainly during DOYs 63–66, 92–93, 102–103, 109–110, 127–128 and 130–131. Continuous precipitation led to a significant nonlinear increase in the SMC. As precipitation decreased or stopped, there was a decrease in the SMC. Evidently, precipitation was the primary factor that caused sudden changes in the SMC. The precipitation at the P041 and MFLE stations during the experimental period was appropriate and suitable for SMC retrieval.

**Figure 4.** The soil moisture rainfall diagram during the experimental period.
