In this section, we mainly describe the proposed method. The selection of the optimal arc is transformed into the solution of the optimal value in two-dimensional coordinates, as shown in Fig. 1a. First, the starting elevation angle and the ending elevation angle of the arcs were determined. Three indicators were selected as important evaluation criteria: the correlation coefficient between SRC and the fitted cosine function curve, the RMSE between SRC and the fitted cosine function curve, and the amplitude of the fitted cosine function curve. Then, the entropy method is used to evaluate the fitting accuracy of the curves for different elevation angle ranges, and the PSO algorithm is used to find the arcs with typical waveforms. These two methods are the focus of this work. Figure 1b shows the process of finding the optimal elevation angle range from 5° to 30°. The proposed method is explained step by step below.
b Process of finding the optimal arc from raw SRC. The cyan dotted lines are the cosine function curves obtained by fitting different arcs. The pink line is the cosine function curve obtained by fitting the arc of the optimal elevation angle range
Step 1-Selecting anomalous tracks based on the H time series
As shown in Fig. 2, when the terrain at the site is flat, the reflector height is close to the antenna height. In contrast, in undulating terrain there is a significant difference between the reflector height and the antenna height. The calculated phases are not suitable for soil moisture inversion when the H time series of a track fluctuates greatly because soil moisture inversion relies on cosine functions with the same angular frequency (Chew et al. 2015). However, within the specified elevation angle range of this track, the SRC is likely to have a typical waveform and we can still calculate the valid SNR metrics from this track. Therefore, we select tracks with large fluctuations in the H time series as research objects.
Step 2-Setting the initial elevation angle values
The starting elevation angle and the ending elevation angle of the SRC are set as the two dimensions of the particles in the PSO algorithm. The ending elevation angle must be T periods greater than the starting elevation angle (T can vary between 1 and 3 depending on terrain conditions). Based on this special relationship, the initial values of the particles are set to (5°, 15°), (7.5°, 17.5°), (10°, 20°), (12.5°, 22.5°), (15°, 25°), (17.5°, 27.5°), (20°, 30°).
Step 3-Fitting the curve using trigonometric least squares method
Classical methods first use the LSP method to extract the reflector height of 5°-30° SRC or directly use the antenna height as the known parameter of the cosine function equation, and then use the nonlinear least squares method to calculate SNR metrics (A and φ). However, for SRC with poor data quality, the reflector height extracted by the LSP method may be abnormal, resulting in poor fitting quality of the SRC and the cosine function curve. Therefore, the LSP method to extract the reflector height is not adopted, and the reflector height is considered as a parameter to be calculated. It should be noted that the SNR metrics calculated here are only used to assess how much the arcs in different elevation angle ranges are affected by noise and are not used for soil moisture inversion. The nonlinear least squares estimation method of trigonometric functions is as follows.
Assuming the closest trigonometric function to the SRC (with discrete coordinates (xi, yi)) is expressed as:
$$F(x)=B \cdot \cos (\omega x)+C \cdot \sin (\omega x)+D$$
7
where ω, B, C, D are the parameters of the trigonometric function. Then for a fixed ω0, using the sum of the squared distances can be expressed as an objective function as follows (Jiao et al. 2023):
$$W(B,C,D)={\sum\nolimits_{{i=1}}^{N} {\left[ {{y_i} - B\cos ({\omega _0}{x_i}) - C\sin ({\omega _0}{x_i}) - D} \right]} ^2}$$
8
where N is the number of discrete points in the SRC. Using the least squares method to find the optimal parameters B, C, D corresponding to ω = ω0, and varying the ω values, the best ω, B, C, D are obtained through iterative solving. That is
$$W({\omega _{{\text{best}}}},{B_{{\text{best}}}},{C_{{\text{best}}}},{D_{{\text{best}}}})=\hbox{min} \left( {W{{({\omega _i})}_{\hbox{min} }}} \right),i=1,2,3 \ldots$$
9
where\({\omega _{{\text{best}}}}\), \({B_{{\text{best}}}}\), \({C_{{\text{best}}}}\) and \({D_{{\text{best}}}}\) are the parameters corresponding to the optimal trigonometric functions. Therefore, the amplitude and phase of the cosine function curve closest to the SRC can be expressed as follows:
$$A=\sqrt {B_{{{\text{best}}}}^{2}+C_{{{\text{best}}}}^{2}} ,\varphi =\arctan \left( {\frac{{{B_{{\text{best}}}}}}{{{C_{{\text{best}}}}}}} \right)$$
10
Step 4-Entropy method for evaluating curve fitting accuracy
Although arcs with different elevation angle ranges can obtain a best-fitting cosine function curve, their fitting accuracy is different. It is necessary to establish objective indicators to evaluate the effectiveness of curve fitting. Here we count the RMSE of the SRC versus the fitted cosine function curve, the correlation value of the SRC versus the fitted cosine function curve, and the amplitude of the fitted cosine function curve for different elevation angle ranges. These three metrics are used as an important basis for evaluating the curve fitting. The following explains the process of evaluating curve fitting accuracy based on the entropy method.
A smaller RMSE means that the fitted cosine function curve is closer to the SRC, while a larger R means that the cosine function curve is more closely correlated with the SRC. In addition, smaller amplitudes affect the evaluation of the curve fitting. Therefore, RMSE is a negative indicator, correlation value and amplitude are positive indicators. They are standardized as follows:
$${y_{ij}}={{(\hbox{max} ({X_j}) - {x_{ij}})} \mathord{\left/ {\vphantom {{(\hbox{max} ({X_j}) - {x_{ij}})} {(\hbox{max} ({X_j}) - \hbox{min} ({X_j}))}}} \right. \kern-0pt} {(\hbox{max} ({X_j}) - \hbox{min} ({X_j}))}}$$
11
$${y_{ij}}={{({x_{ij}} - \hbox{min} ({X_j}))} \mathord{\left/ {\vphantom {{({x_{ij}} - \hbox{min} ({X_j}))} {(\hbox{max} ({X_j}) - \hbox{min} ({X_j}))}}} \right. \kern-0pt} {(\hbox{max} ({X_j}) - \hbox{min} ({X_j}))}}$$
12
where \({x_{ij}}\) is the jth index of the ith particle, \({X_j}\) is the array composed of the jth index of all particles, and \({y_{ij}}\) is the standardized value of the jth indicator of the ith particle. The proportion of the jth indicator of the ith particle is calculated by the formula
$${f_{ij}}=({{{y_{ij}}} \mathord{\left/ {\vphantom {{{y_{ij}}} {\sum\limits_{{i=1}}^{{n+1}} {{y_{ij}}} }}} \right. \kern-0pt} {\sum\limits_{{i=1}}^{{n+1}} {{y_{ij}}} }})$$
13
where n is the number of particles, the n + 1th particle is the global optimal particle. Then the index entropy calculation formula of the jth index is as follows (Zhou et al. 2019):
$${h_j}={{ - \sum\limits_{{i=1}}^{{n+1}} {\left( {{f_{ij}} \cdot \ln {f_{ij}}} \right)} } \mathord{\left/ {\vphantom {{ - \sum\limits_{{i=1}}^{{n+1}} {\left( {{f_{ij}} \cdot \ln {f_{ij}}} \right)} } {\ln \left( {n+1} \right)}}} \right. \kern-0pt} {\ln \left( {n+1} \right)}}$$
14
Then, the weight \({w_j}\) of each indicator is calculated as follows:
$${w_j}={{(1 - {h_i})} \mathord{\left/ {\vphantom {{(1 - {h_i})} {(N - \sum\limits_{{i=1}}^{N} {{h_i}} )}}} \right. \kern-0pt} {(N - \sum\limits_{{i=1}}^{N} {{h_i}} )}}$$
15
where N is the number of evaluation indicators. Finally, the scores of all particles are calculated using the following formula
$${S_i}=\sum\limits_{{j=1}}^{N} {{w_j} \cdot } {y_{ij}}$$
16
We use the particle score to represent the accuracy of the curve fitting. The higher the particle score, the better the curve fitting.
Step 5-Updating the local best and global best of particles
During the iterative process, each particle updates itself by tracking both local and global optima. The local optimum is the best position that a particle searches during the search process, and the global optimum is the best position sought by all particles during the search process. Their updating formulas are as follows (Clerc et al. 2002):
$$lbest_{i}^{{k+1}}=\left\{ \begin{gathered} x_{i}^{{k+1}},S(x_{i}^{{k+1}})>S(lbest_{i}^{k}) \hfill \\ lbest_{i}^{k},S(x_{i}^{{k+1}}) \leqslant S(lbest_{i}^{k}) \hfill \\ \end{gathered} \right.$$
17
$$gbes{t^k}=\arg \mathop {\hbox{min} }\limits_{{1 \leqslant i \leqslant n}} (S(lbest_{i}^{k}))$$
18
where \(lbest_{i}^{k}\) is the local best value of the ith particle during the kth iteration, \(gbes{t^k}\) is the global best value for all particles during the kth iteration, and \(S\left( \cdot \right)\) is the scoring function obtained through the entropy method in formula (16).
Step 6-Updating the velocity and position of particles
Particle velocity and position are updated based on the two optimal values mentioned above. The velocity and position update formulas are as follows (Coello et al. 2004):
$$v_{i}^{{k+1}}=\omega v_{i}^{k}+{c_1}{r_1}(pbest_{i}^{k} - x_{i}^{k})+{c_2}{r_2}(gbes{t^k} - x_{i}^{k})$$
19
$$x_{i}^{{k+1}}=x_{i}^{k}+v_{i}^{{k+1}}$$
20
where \(v_{i}^{k}\) is the updated velocity of the ith particle in the kth iteration, \(x_{i}^{k}\) is the position of the ith particle in the kth iteration, \({c_1}\) and \({c_2}\) are learning factors, and \({r_1}\) and \({r_2}\) are random numbers between 0 to 1. If the maximum iteration count is reached, the optimal elevation angle range is output; otherwise, the process returns to step 3.
The detailed steps for estimating soil moisture using the OAS method are summarized in Fig. 3.
Study areas and data
The International GNSS Service (IGS) and PBO provide high-quality GNSS data for various scientific research projects. Meanwhile, the International Soil Moisture Network (ISMN) maintains a global in-situ soil moisture database. After comparing the locations of soil moisture monitoring stations and GNSS sites, we selected site SGPO (located in Billings, Oklahoma) provided by IGS, as well as the sites KSU1 (located in Manhattan, Kansas) and P183 (located in Manhattan, Kansas) provided by PBO. Information about these 3 GNSS stations is listed in Table 1.
The adjacent soil moisture monitoring stations for sites SGPO, KSU1, and P183 are Lamont_CF1, Manhattan_6_SSW, and Bodega_6_WSW, respectively. The former is operated by the Atmospheric Radiation Measurement, while the latter two are operated by the US Climate Reference Network. Considering the integrity of the satellite data collected by the GNSS receivers and the soil moisture data collected by the probes, 2022 was selected as the analysis period for this study. All three soil moisture monitoring stations provide in-situ measurements from probes at a depth of 5 cm, which are used for validation in soil moisture estimation.
The left column of Fig. 4 shows the surroundings of the sites and the right column shows the slope maps within 100 m of the sites. The slope maps were created from lidar point cloud data provided by the US Geological Survey (USGS). To understand the detailed condition of the terrain around the GNSS stations, the average slope, the maximum height difference, and the height standard deviation were counted for the unit area with the receiver as the center of the circle, a radius of 25m, and a circle angle of 30°, as shown in Table 2. It can be seen from Fig. 4 and Table 2 that the three sites have different levels of terrain complexity, ranked from simple to complex as follows: SGPO, KSU1, P183.
Table 1
Information for sites SGPO, KSU1, and P183
Site name | Receiver type | Antenna type | Utilizing data | Sampling interval | Analysis period | soil texture |
SGPO | JAVAD TRE_3 | JAVRINGANT_G5T | G + C + E + R | 30s | 2022.1.1-12.31 | Silty clay |
KSU1 | SEPT POLARX5 | TRM59800.99 | G + C + E + R | 15s | 2022.1.1-12.22 | Silty clay |
P183 | TRIMBLE NETR9 | TRM59800.80 | G + E + R | 15s | 2022.1.1-12.29 | Coarse sandy loam |
Table 2
Terrain overview of sites SGPO, KSU1, P183
Sites | SGPO | KSU1 | P183 |
Azimuth range (°) | Average slop (°) | Maximum HD (m) | Height STD (m) | Average slope (°) | Maximum HD (m) | Height STD (m) | Average slope (°) | Maximum HD (m) | Height STD (m) |
30–60 | 3.81 | 1.31 | 0.32 | 16.27 | 4.84 | 1.39 | 7.18 | 2.36 | 0.55 |
60–90 | 2.78 | 1.01 | 0.20 | 16.02 | 5.05 | 1.66 | 5.45 | 1.70 | 0.40 |
90–120 | 7.01 | 2.95 | 0.49 | 17.19 | 4.98 | 1.56 | 5.56 | 1.43 | 0.29 |
120–150 | 7.69 | 2.65 | 0.67 | 11.56 | 3.33 | 0.86 | 7.94 | 2.63 | 0.46 |
150–180 | 5.85 | 2.18 | 0.47 | 2.65 | 0.85 | 0.14 | 10.33 | 2.78 | 0.52 |
180–210 | 2.95 | 0.94 | 0.21 | 1.07 | 0.34 | 0.08 | 13.34 | 3.10 | 0.77 |
210–240 | 2.09 | 0.86 | 0.12 | 1.34 | 0.36 | 0.08 | 10.71 | 3.21 | 0.76 |
240–270 | 1.02 | 0.30 | 0.06 | 1.25 | 0.23 | 0.06 | 8.06 | 3.18 | 0.66 |
270–300 | 2.11 | 0.88 | 0.13 | 1.06 | 0.20 | 0.05 | 7.44 | 3.22 | 0.73 |
300–330 | 1.91 | 0.47 | 0.10 | 1.34 | 0.39 | 0.09 | 10.79 | 3.16 | 0.72 |
HD Height difference STD Standard deviation |