Estimating Soil Moisture over Winter Wheat Fields during Growing Season Using RADARSAT-2 Data

Abstract

Soil moisture content (SMC) is a significant factor affecting crop growth and development. However, SMC estimation, based on synthetic aperture radar (SAR), is influenced by a variety of surface parameters, such as vegetation cover and surface roughness. As a result, determining the SMC across agricultural areas (e.g., wheat fields) remotely (i.e., without ground measurement) is difficult to achieve. In this study, a model-based polarization decomposition method was used to decompose the original SAR signal into different scattering components that represented different scattering mechanisms. The different volume scattering models were applied, and then the results were compared in order to remove the scattering contribution from vegetation canopy, and extract the surface scattering components related to the soil moisture. Finally, by combining extensively used surface scattering models (e.g., CIEM and Dubois), and a method of roughness parameters optimization, a lookup table was developed to estimate the soil moisture during the wheat growth period. When CIEM is applied, the R² and RMSE of the SMC are 0.534, 5.62 vol.%, and for the Dubois model, 0.634, 5.16 vol.%, respectively, which indicates that this approach provides good estimation performance for measuring soil moisture during the wheat growing season.

1. Introduction

Soil moisture content (SMC) is a crucial parameter of the land surface ecological water cycle, which controls surface runoff and the evaporation of surface water [1,2,3,4]. In agricultural regions, SMC is an important water source for crops, is associated with their survival, and is also a key factor affecting soil fertility [5,6,7,8]. Therefore, farmland SMC is used as an effective index of farmland health. The timely acquisition of soil moisture across a wide range of farmlands is of great significance to guide agricultural production. There are two traditional soil moisture monitoring methods. One method uses a soil moisture meter (e.g., a TDR probe) to measure the volumetric SMC, by inserting it into the soil at a certain depth. Another method is to submit the soil sample for laboratory analysis, through the drying processing, to measure the gravimetric SMC [9]. However, since the traditional methods are time consuming and laborious, they are not practical for large-scale SMC monitoring [10,11].

The emergence of remote sensing technology provided an effective way to estimate soil moisture from a limited representative “point” to a regional “area” [12]. However, since optical images are often affected by weather conditions, such as fog and clouds, in contrast, the usage of microwave remote sensing is preferred, since it has stronger penetration abilities, which are less affected by light, clouds, and fog, and it registers all weather conditions. More importantly, surface microwave radiation is sensitive to moisture changes in the soil, so microwave remote sensing has gradually become the main means of monitoring soil moisture [5,13,14]. The methods of microwave remote sensing for monitoring SMC are divided into active and passive microwave remote sensing technologies. Passive microwave remote sensing is based on the radiometer to construct an SMC inversion method. At present, the existing spaceborne microwave radiometers include, among others, AMSR-E and SMOS [15,16]. Passive microwave remote sensing has a high temporal resolution, but its spatial resolution is low, which limits its ability to accurately monitor soil moisture in agricultural areas [17]. Active microwave remote sensing is based on radar, or a scatterometer, for SMC estimation [18,19]. Its temporal resolution is lower than that of passive remote sensing, but its spatial resolution is higher. In addition, SAR benefits from its higher spatial resolution, and provides multi-angle and multi-polarization data, which are the main means of monitoring SMC in active microwave remote sensing fields [14,20,21,22,23,24].

Many SAR-based scattering models were developed for SMC estimation of bare surfaces. This study focused on constructing the relationship between soil moisture and radar backscattering coefficients [2]. Representative models were empirical/semi-empirical models, such as the Oh and Dubois models [25,26,27], and physical models, such as IEM [28]. Machine learning models, such as the random forest (RF) model and the support vector machine (SVM), were used for soil moisture estimation based on SAR parameters [29,30,31]. For most physical models and empirical/semi-empirical models, the construction of these models often required many parameters, such as radar incidence angle, dielectric constant, and roughness parameters (e.g., correlation length (L) and root mean square (RMS) height). Roughness parameters were not easy to obtain accurately, especially in dense vegetation areas, resulting in high levels of uncertainty [32]. To solve the problem, some studies replace certain roughness parameters that are difficult to measure with fitted roughness parameters, or by using a single roughness parameter value as the optimal roughness parameter for the whole study area [5,33,34]. Additionally, in highly vegetated areas, the vegetation canopy not only generates volume scattering, but also attenuates scattering from the surface, resulting in a mixed backscattering consisting of components from the vegetation, the surface, and the interaction between the vegetation and the ground surface. Moreover, when the vegetation cover exceeded a threshold, it highly attenuated the SAR signal, so that the scattering contribution from the soil received by the sensor was negligible [24,35,36,37]. This made the soil moisture estimation under vegetation difficult to measure [38]. Generally, in order to estimate the SMC under vegetation cover, the volume scattering from the vegetation canopy should be considered; therefore, the common solution is to couple vegetation models to differentiate the contributions of bare soil and vegetation [39,40], but such models also introduce more unknown parameters.

SAR polarization decomposition technology decomposes the full polarization signal into different scattering components and, thus, eliminates the scattering from vegetation, which is gradually applied to the SMC retrieval [41]. Different scholars tried to use different polarization decomposition methods to estimate soil moisture, such as H-A-α decomposition [42,43,44] and Freeman–Durden decomposition [45]. Facts prove that Freeman–Durden decomposition shows potential in soil moisture retrieval [45]. Hajnsek et al., first establish the relationship between the surface and dihedral scattering components and soil moisture through the Freeman–Durden decomposition [46], and then estimate the soil moisture of three crops during the full vegetation growth period. The achieved average root mean square error (RMSE) is about 10 vol.%. Jagdhuber et al., develop a multi-angular physical-based decomposition method to retrieve soil moisture across a bare surface with sparse vegetation, using the X-Bragg model and Freeman–Durden decomposition with a high inversion rate and low RMSE [47]. Wang et al., estimate the soil moisture under different crop covers, based on the surface and volume scattering of Freeman–Durden decomposition. In their research, the dihedral scattering is ignored, so the estimation of the soil moisture is simplified. Finally, the correlation coefficient and the estimated RMSE during the entire phenological period are 0.63–0.76 and 5.8–7.4 vol.%, respectively [48].

In summary, if the penetration of SAR signal is sufficient, two issues need to be resolved to estimate soil moisture based on SAR data: removing the contribution of vegetation scattering, and obtaining the required surface parameters. In view of these challenges, this study proposed a soil moisture retrieval method combining Freeman–Durden decomposition with different surface scattering models (e.g., CIEM and Dubois), and an optimal roughness parameter, to estimate soil moisture at different stages of wheat growth. The polarization decomposition method was used to remove volume scattering from the vegetation canopy. In the absence of the measured roughness parameters, the optimal roughness parameter was used to reduce the dependence of the model on the measured roughness parameters. Therefore, this study had the following two objectives: (1) investigating the feasibility of the proposed method for estimating soil moisture in wheat during different growth stages, and provide an effective means for monitoring crop soil moisture; and (2) investigating and comparing the performances of different surface scattering models (e.g., CIEM and Dubois) in soil moisture estimation.

2. Materials

2.1. Study Area

The study area is located in southwestern Ontario, Canada (Figure 1). The selected area is a rain-fed agricultural area. The main crops include winter wheat, corn, and soybeans. In this study, the estimation method was verified only for the winter wheat-growing area. An L-shaped wheat growing area of approximately 27 hectares was selected as the study region (Figure 1). The sowing date of winter wheat in the study area is generally from October to November of the previous year, and enters the germination period approximately a week after sowing. From the end of December to January of the next year, due to the decrease in temperature, the wheat ceases to grow and enters the overwintering period. When temperatures begin to rise in March of the following spring, wheat enters the returning green period and starts to grow. Tillering and stem elongation stages typically commence in May when the wheat grows rapidly. The booting and heading stages occur in June, and the ripening and harvest stages are in July.

2.2. Data Used

2.2.1. Ground Truth Data Collection

During the growth stage of winter wheat in the study area, we conducted a total of 8 field surveys from May to July in 2019 [49]. The specific sampling date and the measured SMC range on each sampling date are listed in

Table 1. Collected SMC details.

Measurement Date	Growth Stage	Soil Moisture Range (Vol.%)	Mean Values (Vol.%)	Std. Dev. (Vol.%)
9 May 2019	Tillering	20.88–40.20	33.20	5.40
16 May 2019	Tillering	14.20–37.20	26.76	6.24
20 May 2019	Stem elongation	14.43–36.02	26.42	5.88
29 May 2019	Stem elongation	17.47–35.45	27.42	5.08
2 June 2019	Booting	17.65–34.25	27.40	4.96
9 June 2019	Heading	11.35–31.52	21.36	5.49
16 June 2019	Grouting	22.02–37.87	28.87	3.98
10 July 2019	Ripening	4.17–20.97	8.98	4.04

. Except for the number of sampling points (16) on the last sampling date (July 10), the number of sampling points on the remaining dates was 32, and a total of 240 sampling points were obtained. For the measurement of SMC at each sampling site, the spacing between any adjacent sampling points was at least 30 m. The surface soil moisture was measured using the Theta Probe soil moisture sensor, at a depth of 5 cm. In order to avoid accidental errors, we measured, on average, six times at each sampling point, and took their mean value as the final measured SMC of the sampling point. Meanwhile, real-time kinematic (RTK) was used to record the longitude and latitude information, so that the parameter information corresponding to the sampling point position could be extracted from the geocoded SAR images. The location of the RTK device was 1.5 cm.

2.2.2. RADARSAT-2 Data Collection and Preprocessing

RADARSAT-2 was equipped with a C-band sensor, with a center frequency of approximately 5.4 GHz. The nominal spatial resolution of the RADARSAT-2 images was approximately 8 m. The satellite revisit period was 24 days, and the orbit period was 100.7 min. In this study, we obtained 8 fine quad-polarized RADARSAT-2 images covering the study area, and the image acquisition date was consistent with the sampling date.

Table 2. Characteristics of collected RADARSAT-2 data.

Acquisition Date	Orbit	Acquisition Time (UTC)	Central Incidence Angle (Study Area)
9 May 2019	Ascending	23:14	30.0°
16 May 2019	Ascending	23:10	25.4°
20 May 2019	Descending	11:34	34.2°
29 May 2019	Ascending	23:30	48.2°
2 June 2019	Ascending	23:14	30.0°
9 June 2019	Ascending	23:10	24.2°
16 June 2019	Ascending	23:05	19.3°
10 July 2019	Ascending	23:05	19.3°

shows the acquisition date and time, as well as the orbit and incidence angles of these 8 SAR images.

The preprocessing of RADARSAT-2 data was completed using SNAP 7.0 software, provided by ESA, and included the following steps:

(1)

Radiometric calibration: the quad-polarization complex images of backscattering coefficient were generated;

(2)

Based on the complex images obtained in step 1, we generated the coherency matrix (T3 matrix);

(3)

Polarization filtering: the refined Lee filter was used to reduce the noise generated in step 2, and its window size was set to 7 × 7;

(4)

Terrain correction and projection transformation: the filtered data was geocoded with an SRTM-3 digital elevation model, and the output coordinate system was projected to WGS84.

3. Methods

Throughout this study, the Freeman–Durden decomposition technique, optimal roughness parameters, and surface scattering models (e.g., CIEM and Dubois) were used to estimate the SMC across the wheat field in our study area. Figure 2 shows the flowchart of the SMC estimation.

(1)

RADARSAT-2 data were preprocessed, and T3 matrix was obtained;

(2)

Scattering mechanism analysis was used to analyze the main scattering mechanisms of winter wheat in different growth stages in the study area;

(3)

The backscattering coefficient of the bare surface was calculated based on the surface scattering matrix (TG). The surface backscattering coefficient and the measured SMC were used to construct the SMC estimation dataset. For each image, we randomly selected 70% as the training set, and the remaining 30% as the validation set. In order to avoid the statistical differences of soil moisture between the training set and the validation set, which would affect the accuracy from the validation set, a selection of statistical characteristics of the training and the validation sets were calculated. Figure 3 shows that the mean, median, and standard deviation of the training and validation sets are significantly close to the statistical characteristics of the total data set, which meets the requirements for data splitting;

(4)

On the training set, based on the surface backscattering coefficient simulated by the CIEM or the Dubois model under different roughness conditions in a given range, and the surface backscattering coefficient obtained by polarization decomposition, the cost function between them was constructed to estimate the SMC of each sampling point;

(5)

The roughness parameters were optimized when the RMSE between the estimated SMC and the measured SMC on the training set was minimized;

(6)

The simulated backscattering coefficient from the CIEM or the Dubois model with optimal roughness parameters was obtained, and then a same cost function was constructed between the simulated backscattering coefficient and the backscattering coefficient obtained by polarization decomposition on the validation set. Then, the estimated SMC was retrieved on the validation set, and the accuracy was also verified on the validation set;

(7)

The model with the highest accuracy on the validation set was selected to draw the regional SMC map of the winter wheat area.

3.1. Polarimetric Decomposition Methods

3.1.1. Cloude–Pottier Decomposition

Based on the eigenvector analysis of the polarization coherence matrix, Cloude and Pottier propose a decomposition theorem that includes all scattering mechanisms [50]. The coherence matrix is decomposed into the weighted sum of three components, using eigenvalue decomposition, each of which corresponds to a scattering mechanism. By calculating the eigenvalues and eigenvectors from the decomposition, three parameters are obtained: polarization entropy “H”, average scattering angle “α”, and anti-entropy “A”; the values of H and α provide a way to understand the mechanism of scattering [51]. The calculation formulas are as follows:

$$\begin{array}{l}H={\displaystyle \sum _{i=1}^3-p_i{\log }_3}{p}_i\\ \alpha ={\displaystyle \sum _{i=1}^3{p}_i\arccos (|n_{i1}|)}\\ p_i=\frac{{\lambda }_i}{{\displaystyle \sum _{i=1}^3{\lambda }_i}}\\ ni=[n_{i1} n_{i2} n_{i3}]{ }^T\end{array}$$

(1)

where ${\lambda }_i$ is the eigenvalue of coherence matrix T3, and $n_i$ is the eigenvector.

3.1.2. Freeman–Durden Decomposition

Freeman and Durden propose a three-component scattering model based on VanZyl’s research [52]. The main idea of the model is to decompose the polarization coherence matrix (T3) into three main scattering mechanisms, namely, the scattering from the surface (TG), volume scattering from the vegetation canopy (TV), and dihedral scattering (TD) generated by the interaction between the vegetation and the surface. T3 can be decomposed into the following form:

$$\begin{array}{ll}T3& =\left[\begin{array}{ccc}{T}_{11}& T_{12}& 0\\ T_{12}^{*}& T_{22}& 0\\ 0& 0& T_{33}\end{array}\right]=TG+TD+TV\\ & =f_s\left[\begin{array}{ccc}1& {\beta }^{*}& 0\\ \beta & {\left|\beta \right|}^2& 0\\ 0& 0& 0\end{array}\right]+f_d\left[\begin{array}{ccc}{\left|\alpha \right|}^2& \alpha & 0\\ {\alpha }^{*}& 1& 0\\ 0& 0& 0\end{array}\right]+f_{v}V\end{array}$$

(2)

where the “*” represents the conjugate transpose; TG is simulated by the Bragg model: β denotes the normalized difference of the Bragg scattering in H and V polarization, and $f_s$ is the intensity coefficient of the surfaces scatter component matrix; TD is simulated by the Fresnel model, where α is the normalized difference of the combined ground-stalk Fresnel reflection in H and V polarizations, and $f_d$ is the intensity coefficient of the dihedral scatter component matrix; and TV depends on the selected form of the volume scattering matrix ([V]), and the intensity coefficient of [V] ($f_v$). Yamaguchi proposes three kinds of volume scattering matrices for different vegetation orientations (i.e., vertical, random, and horizontal) according to the P_r value [46,53].

$$P_r=10\cdot \log \frac{<|S_{VV}{|}^2>}{<|S_{HH}{|}^2>}$$

(3)

$$\begin{array}{l}\mathrm{Tv}=f_{v}V=f_v\left[\begin{array}{ccc}{V}_{11}& V_{12}& 0\\ V_{12}& V_{22}& 0\\ 0& 0& V_{33}\end{array}\right]\to { \mathrm{Tv}}_v=\frac{f_v}{30}\left[\begin{array}{ccc}15& 10& 0\\ 10& 8& 0\\ 0& 0& 7\end{array}\right],\\ {\mathrm{Tv}}_r=\frac{f_v}{4}\left[\begin{array}{ccc}2& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right],{ \mathrm{Tv}}_h=\frac{f_v}{30}\left[\begin{array}{ccc}15& -10& 0\\ -10& 8& 0\\ 0& 0& 7\end{array}\right]\end{array}$$

(4)

where Tv_v, Tv_r, and Tv_h represent the volume coherence matrix components of vertical, random, and horizontal vegetation orientations, respectively. In this study, the estimation performances of SMC under these three volume scattering matrices were compared.

The power of each scattering component was obtained as follows, based on the decomposition result of T3 [54]; P_s, P_d, and P_v correspond to the power of surface scattering, dihedral scattering, and volume scattering, respectively.

$$\begin{array}{c}{P}_s=f_s(1+|\beta {|}^2)\\ P_d=f_d(1+|\alpha {|}^2)\\ P_v=f_v\end{array}$$

(5)

When the dihedral scattering component is ignored, TG is expressed in the following form:

$$TG=T3-f_{v}V$$

(6)

When the volume scattering matrix was selected, the nonnegative eigenvalue decomposition (NNED) method was used to determine the value of $f_v$ [55], which met the condition that the minimum eigenvalue of TG was zero after removing the dihedral scattering component. Based on the coherency matrix corresponding to the surface scattering component, the surface backscattering coefficient was obtained [51,56].

The local incidence of SAR usually has a significant effect on the surface backscattering coefficient, which cannot be ignored in soil moisture estimation. In order to limit the influence of different local incidence angles, the theoretical method proposed by Ulaby et al. [57] was used to reduce the influence of local incidence angles after decomposition:

$${\sigma }_{{\theta }_{ref}}^{\circ }={\sigma }_{\theta }^{\circ }\times \frac{{\cos }^2{\theta }_{ref}}{{\cos }^2\theta }$$

(7)

where θ is the incidence angle, ${\theta }_{ref}$ is the reference incidence angle (30°), ${\sigma }_{\theta }^{°}$ is the surface backscattering coefficient obtained by polarization decomposition, and ${\sigma }_{\theta ref}^{°}$ is the normalized backscattering coefficient according to the reference incidence angle. In this study, the reference incidence angle was set to 30°.

3.2. Soil Surface Backscattering Modeling

3.2.1. Dubois Model

Dubois et al. [26] establish an empirical scattering model based on multi-frequency, multi-polarization, and multi-angled scatterometer data measured on the ground with different surface roughness [26], and its specific expression is as follows:

$$\begin{array}{l}{\sigma }_{hh}^{\circ }={10}^{-2.75}\frac{{\cos }^{1.5}\theta }{{\sin }^5\theta }{10}^{0.028\epsilon \tan \theta }(ks\times \sin \theta ){ }^{1.4}{\lambda }^{0.7}\\ {\sigma }_{vv}^{\circ }={10}^{-2.35}\frac{{\cos }^3\theta }{{\sin }^3\theta }{10}^{0.046\epsilon \tan \theta }(ks\times \sin \theta ){ }^{1.1}{\lambda }^{0.7}\end{array}$$

(8)

The equation above establishes the functional relationship between the backscattering coefficient and surface parameters (e.g., RMS height (s) and soil dielectric constant (ε)), and the sensor configuration parameters (e.g., incident angle (θ) and wavelength (λ)). Where ${\sigma }_{hh}^{°}$ and ${\sigma }_{vv}^{°}$ are the backscattering coefficients of the bare land surface under HH and VV channels simulated by the Dubois model, respectively, k is the free space wavenumber (k = 2π/λ).

The Dubois model is not suitable for the rough surface conditions (ks > 2.5) [58]. Therefore, the RMS height of the lookup table constructed by the Dubois model in this study is between 0–2.2 cm. The SMC was calculated from the ε using the Topp model, shown below [59]:

$$m_v=-5.3\times {10}^{-2}+2.92\times {10}^{-2}\epsilon -5.5\times {10}^{-4}{\epsilon }^2+4.3\times {10}^{-6}{\epsilon }^3$$

(9)

3.2.2. CIEM Model

The integral equation model (IEM) is a theoretical radar backscatter model with a high validity range [28]. The IEM model establishes the functional relationship between the co-polarized backscattering coefficient and the roughness parameters (correlation length (L), RMS height (s)), dielectric constant, incident angle, and wavelength, etc.). Given the uncertainty of the field measurement of the surface correlation length (L), Baghdadi et al., (2006) establish the semi-empirical relationship between the corrected L and the s, based on abundant SAR images and measured data based on an IEM model, reduce the parameters describing the surface roughness of the soil, and develop a calibrated IEM model (CIEM). The empirical relationship is expressed as [60]:

$$L(s,\theta ,pp)=\delta {(\sin \theta )}^{\mu }{s}^{(\eta \theta +\xi )}$$

(10)

where L is the corrected correlation length, s is the RMS height, θ represents the incident angle, and pp represents the polarization mode (HH, VV). The coefficients “$\delta$” and “$\xi$” are related to the polarization mode, and the coefficients “$\mu$” and “$\eta$” are independent of the polarization mode. The specific coefficients of the semi-empirical correction model for C-band SAR data are as follows [60]:

$$\begin{array}{l}{\delta }_{hh}=4.026,{\xi }_{hh}=1.551,{\delta }_{vv}=3.289,{\xi }_{vv}=1.222\\ {\mu }_{hh}={\mu }_{vv}=-1.774,{\eta }_{hh}={\eta }_{vv}=-0.0025\end{array}$$

(11)

In this study, the RMS height range of the lookup table constructed by CIEM model is consistent with the Dubois model, i.e., 0–2.2 cm.

3.2.3. Optimum Surface Roughness Parameter

Roughness is an important parameter in the surface scattering model, but it is difficult to measure, especially in areas with dense vegetation. Due to the lack of measured roughness parameters in the study area, the surface scattering models cannot be used to simulate the backscattering coefficients of bare soil surface with the actual roughness of sampling points. Therefore, the optimum surface roughness parameter proposed by Bai et al., (2016) [5] was used to parameterize the surface scattering model. This method assumes that the roughness parameter in the study area is a fixed value. SMC is estimated within a certain roughness range. The roughness parameter is taken as the optimum surface roughness parameter of the study area when the RMSE between the estimated SMC values and the measured SMC values is minimized. In this study, only the RMS height is used as the roughness parameter in the surface scattering models.

3.3. Soil Moisture Estimation and Performance Assessment

From the training set, the ground backscattering coefficient (${\sigma }_{ground}^{°}$) was obtained, using the polarization decomposition. Under different roughness conditions, the lookup table constructed by the Dubois and CIEM models is obtained: the simulated ground backscattering coefficient at different SMC conditions. The following cost function was constructed to estimate SMC:

$$F=\left\{\begin{array}{l}\sqrt{{({\sigma }_{vv}^{\circ }(ground)-{\sigma }_{vv}^{\circ }(simu))}^2}, \left[\mathrm{F}\right(1\left)\right]\\ \sqrt{{({\sigma }_{hh}^{\circ }(ground)-{\sigma }_{hh}^{\circ }(simu))}^2}, \left[\mathrm{F}\right(2\left)\right]\\ \sqrt{{({\sigma }_{vv}^{\circ }(ground)-{\sigma }_{vv}^{\circ }(simu))}^2+{({\sigma }_{hh}^{\circ }(ground)-{\sigma }_{hh}^{\circ }(simu))}^2}, \left[\mathrm{F}\right(3\left)\right]\end{array}\right.$$

(12)

F (1), F (2), and F (3) represent the cost functions constructed using the VV polarization backscattering coefficient, HH polarization backscattering coefficients, and the combination of VV polarization and HH polarization backscattering coefficients, respectively. When the cost function is minimized, the SMC value corresponding to the backscattering coefficient simulated by the scattering model is taken as the SMC value of the sample point.

The performance of different models was compared and evaluated using statistical indicators, including the coefficient of determination (R², Equation (13)) and root mean square error (RMSE, Equation (14)). The higher the R² and the lower the RMSE mean, the better the fitting effect is between the estimated and the measured values, with better model estimation performance and higher accuracy.

$${\mathrm{R}}^2=1-\frac{{\displaystyle \sum _{i=1}^n{(y_i-{\widehat{y}}_i)}^2}}{{\displaystyle \sum _{i=1}^n{(y_i-\overline{y})}^2}}$$

(13)

$$\mathrm{RMSE}=\sqrt{{\displaystyle \sum _{i=1}^n\frac{{(y_i-{\widehat{y}}_i)}^2}{n}}}$$

(14)

where n is the number of sampling points; $y_i$ and ${\widehat{y}}_i$ represent the measured SMC value and estimated SMC value of the i-th sample site, respectively; and $\overline{y}$ is the mean values of the measured SMC values.

4. Results and Discussion

4.1. Scattering Mechanisms Analysis

The H/α scatter diagram obtained using the Cloude–Pottier decomposition of the pixels of the sampling points in the study area is shown in Figure 4.

As shown in Figure 4, the scatter points of wheat sampling points during the growth period are relatively concentrated in the H/α scatter diagram. The scattered points are distributed primarily among the Z2, Z5, and Z6 areas of the graphs, which indicates that the backscattering composition of the study area is mainly composed of vegetation body scattering, and surface scattering attenuated by vegetation. Therefore, the bare surface scattering model is not directly applied to the soil moisture inversion of wheat-covered farmland. As the wheat grows and develops, its canopy coverage increases gradually, so its volume scattering grows stronger and stronger, which is consistent with the trend that the scattered points show, gradually moving from the Z6 area to the Z5 and Z2 areas. On 16 May, as compared to 9 May, the scattered points are concentrated in the Z6 area, which is due to the incident angle of the SAR image on 16 May being lower than it is on 9 May, and the surface scattering is dominant at a low incident angle [61].

In order to further analyze the scattering mechanism, the percentages of the surface scattering components, the dihedral scattering, and the volume scattering during different phases are obtained based on Freeman–Durden decomposition, as shown in Figure 5. Each point corresponding to the left Y-axis in the figure is obtained by calculating the mean value of different scattering component proportions of all the sampling points in the study area on a certain date. For each sampling date, we calculate the mean value of the volume, surface, and dihedral scattering proportions of all the sampling points. The NDVI in the figure is calculated based on six cloud-free Sentinel-2 images, with acquisition dates near the sampling time. As shown in Figure 5, NDVI first indicate a rapid growth, followed by a slower increase, and then a rapid decrease on the last date. This is consistent with the growth process of winter wheat. The rapid decline in the NDVI on the last sampling date is due to the wheat being fully mature and yellow, resulting in a weakening of the absorption of the red light band, but the coverage of the wheat is still high. On the first two sampling dates, when the NDVI is low, the surface scattering ratio is high, while with the increase in the NDVI, or vegetation cover, the surface scattering ratio decreases significantly. This is consistent with previous studies; the increase in wheat coverage results in more backscattering attenuation [36]. The proportion of scattering components obtained by polarization decomposition is consistent with the growth of the vegetation. The proportion of the dihedral scattering energy is deficient during the wheat growth period, and the backscattering is from both volume scattering and surface scattering. Therefore, ignoring the dihedral scattering component simplifies the calculation of the backscattering coefficient by polarization decomposition, as well as the subsequent estimation of soil moisture. As shown in Figure 5, on 16 May, the surface scattering ratio increases, and the volume scattering ratio decreases, which is consistent with the results shown in Figure 4, due to the influence of the low incidence angle. In general, the proportion of the volume scattering components of wheat in the study area is significantly larger than that of the surface scattering components. Therefore, wheat in growing season has a strong attenuation effect on the total SAR signal, resulting in a relatively low proportion of scattering components reaching the surface.

4.2. Simulated Backscattering Coefficient

Figure 6 illustrates the backscattering coefficients (VV, HH) simulated using the CIEM and Dubois models at a 30° incident angle, with different RMS heights and SMCs. It is evident that the backscattering coefficients simulated by the two surface scattering models increase as the RMS height or SMC content increase. The simulated values of the VV polarization are higher than that of the HH polarization when the roughness and soil moisture content are the same. When the roughness is very low, the two models are not sensitive to the changes in soil moisture. As compared with the CIEM model, the simulated backscattering coefficient of the Dubois model is not sensitive to the change in soil moisture when it is low (SMC < 10 vol.%). However, when the roughness is higher than a certain range, it is very sensitive to the change in soil moisture at a high SMC (SMC greater than 30 vol.% or 40 vol.%). Therefore, when the soil is too dry or too wet, there is uncertainty.

4.3. CIEM and Dubois Estimation Results

In this study, the dihedral scattering component is ignored due to its low contribution. Different volume scattering matrices (e.g., vertical, random, horizontal, and volume matrix selected according to the P_r value) are selected to compute the surface backscattering coefficient based on polarization decomposition. After solving the surface backscattering coefficients, and combining the scattering models (i.e., CIEM and Dubois), the optimal roughness parameter, and the corresponding cost function, the soil moisture of the validation set is estimated according to the workflow shown in Figure 2. The results are listed in

Table 3. Estimation results using CIEM model.

CIEM	F (1)		F (2)
	R2	RMSE (Vol.%)	R2	RMSE (Vol.%)
vertical	0.400	7.02	0.156	18.24
horizontal	0.166	19.09	0.360	8.71
random	0.142	16.46	0.183	12.94
Pr-based	0.029	16.09	0.046	19.56

and

Table 4. Estimation results using Dubois model.

Dubois	F (1)		F (2)
	R2	RMSE (Vol.%)	R2	RMSE (Vol.%)
vertical	0.510	5.80	0.151	18.65
horizontal	0.181	20.82	0.480	8.00
random	0.213	15.96	0.211	14.14
Pr-based	0.154	17.11	0.052	19.82

. The horizontal label in the table represents which cost function in Equation (12) is used, and the vertical label represents which volume scattering matrix selection is used. Furthermore, Table 3 reports the estimated results of the validation set using the CIEM model. When using the fixed volume scattering matrix, including either the vertical, random, or horizontal volume scattering matrix, a relatively high estimation accuracy is obtained on the test set only in the following two cases, as shown in bold in Table 3: when the vertical volume scattering matrix is combined with the cost function F (1), and when horizontal volume scattering matrix is combined with F (2), and the estimation accuracy is R² = 0.400 and RMSE = 7.02 vol.%; and R² = 0.360 and RMSE = 8.71 vol.%, respectively. In other cases, including the strategy of dynamically selecting the volume scattering matrix according to the P_r value, the estimation effect on the test set is poor, which suggests a lower R value and a higher RMSE value. Table 4 reports the results of the Dubois model on the validation set, and the same results as the CIEM model are obtained. When the vertical volume scattering matrix is combined with F (1), and the horizontal volume scattering matrix is combined with F (2), a relatively high estimation accuracy is obtained, which are R² = 0.510 and RMSE = 5.80 vol.%; and R² = 0.480 and RMSE = 8.00 vol.%, respectively.

Figure 7a,b illustrate the field mean of the surface backscattering coefficients of the VV polarization and HH polarization obtained by the polarization decomposition method using the vertical volume scattering matrix and the horizontal volume scattering matrix, respectively. As shown in Figure 7a, the surface backscattering coefficient of the HH polarization is very low, as compared to the VV polarization, which is basically less than −16 dB. Therefore, there is a large estimation error when using the lookup table and cost function constructed under different roughness conditions to estimate soil moisture. For example, when the RMS height is low, the simulated backscattering coefficient is not sensitive to the change in soil moisture, and the estimated SMC is more likely to be larger or smaller than the actual value. When the RMS height exceeds the threshold, the minimum value of the simulated backscattering coefficient exceeds −16 dB, as shown in Figure 6, and, therefore, the soil moisture is significantly underestimated, due to the low backscattering coefficient. Moreover, Figure 7b reflects the same problem; after using the horizontal volume scattering matrix, the value of the surface backscattering coefficient of VV polarization is very low, resulting in a large estimation error. The estimation accuracies of the other two volume scattering matrix selection strategies (i.e., random and P_r-based) are relatively low. However, the calculated surface backscattering coefficient is low as a result of not being within the effective estimation range of the lookup table. Therefore, those two matrix selection strategies are not be suitable for extracting the surface backscattering coefficient of the study area in order to estimate the SMC.

According to the estimation results of the validation set, the two scattering models obtain better estimation performances for the VV surface backscattering coefficient calculated based on the vertical volume scattering matrix, and the HH surface backscattering coefficient calculated based on the horizontal volume scattering matrix. Therefore, these two scattering components are regarded as effective surface backscattering coefficients. Therefore, we combine the two effective surface scattering components, and then estimate the soil moisture of the validation set based on the optimal roughness parameter and the cost function F (3) in Equation (12); the results are shown in Figure 8.

As shown in Figure 8, the accuracy of combining the VV and HH effective backscattering components is higher than using only one of them on the validation set: the accuracy of the CIEM model is R² = 0.534 and RMSE = 5.62 vol.%, and the accuracy of the Dubois model is R² = 0.634 and RMSE = 5.16 vol.%. This is likely due to the combination of the HH and VV effective surface scattering that contains more backscattering information than a single polarization, which reduces the estimation error of the scattering model under the optimal roughness parameters. Overall, the estimation performance of the Dubois model is better than that of the CIEM model. To show a further comparison of the differences between the two models, Figure 9 displays the estimation accuracy (RMSE) of the CIEM and Dubois models on different sampling dates. As shown in Figure 9, the accuracy of the Dubois model estimation is generally better than that of the CIEM model on seven sampling dates, from 9 May to 16 June. Specifically, the RMSE is close to, or lower than, that of the CIEM. However, on the last sampling date, the accuracy of the Dubois model is significantly lower than that of the CIEM model. This is due to the SMC at the sampling point on 10 July being low (the average is less than 10 vol.%). As shown in Figure 6, when the SMC is low, the backscattering coefficient simulated by the Dubois model is not sensitive to the change in soil moisture, so there is a large estimation error.

4.4. SMC Dynamics during Wheat Growing Seasons

The dynamic changes in the measured and estimated values on different dates in the validation set when the CIEM or Dubois model is applied to achieve the best SMC estimation performance are shown in Figure 10. The estimated SMC can better restore the change trend of measured soil moisture. There is no artificial irrigation for wheat in the study area, so the SMC is only affected by precipitation, vegetation absorption, and soil evapotranspiration. In order to better explain the SMC temporal variability in Figure 10 and the soil moisture maps, the rainfall data in the study area is obtained from the ERA5-Land dataset (https://www.ecmwf.int/ (accessed on 20 December 2021)). ERA5-Land provides a consistent view of the water and energy cycles at the surface level for several decades. It contains a detailed record from the year 1950 onwards, with a temporal resolution of 1 h. The native spatial resolution of the ERA5-Land reanalysis dataset is 9 km on a reduced Gaussian grid. The daily precipitation from 1 May to 10 July in the study area is shown in Figure 10.

As shown in Figure 10, there is obvious precipitation in the two days prior to 9 May. In addition, the wheat plants are sparse, and water absorption is low, so the SMC on 9 May is high. From 16 May to 2 June, the mean values of SMC are roughly the same. This is most likely due to less rainfall. The soil moisture supplemented by rainfall reaches a dynamic balance with the soil water absorbed by wheat and consumed by soil evapotranspiration. There is also obvious heavy precipitation on 2 June, but the soil moisture does not increase significantly. This is due to the precipitation time occurring after the sampling was completed. From 2 June to 16 June, the mean SMC first decreases and then increases, probably due to the difference in precipitation near the sampling date. On the last sampling date, the soil moisture is the lowest, which could be due to the high temperatures in July, the strong soil evapotranspiration, and the lack of precipitation a few days before sampling.

4.5. Variation and Analysis of Optimal Roughness Parameters

In order to verify whether the optimal roughness parameter adjusts the estimation accuracy of soil moisture, and the feasibility of parameterized surface scattering models, a time-series dynamic variation of the optimal roughness parameter is obtained. The roughness parameters of the scattering model used in this study are the RMS height of the soil surface. According to the above estimation results, it is found that better estimation performance is obtained under the following scattering component conditions: (1) an effective VV scattering component, (2) an effective HH scattering component, and (3) a combination of the VV and HH scattering components. Therefore, Figure 11 shows the variation of the optimal roughness parameter (RMS height) when estimating soil moisture under the above three conditions. The estimation results of the two scattering models indicate that the RMS height first presents a decreasing trend, and then increases. The RMS heights on the first and last days are higher than those on the other days. For each sampling date, the RMS height based on the VV effective scattering component is the largest, the RMS height based on the HH effective scattering component is the smallest, and the RMS height based on the combination of VV and HH is somewhere between the first two.

Figure 12 shows the mean values of effective surface scattering and measured soil moisture on different sampling dates. We find that the HH effective surface scattering component based on polarization decomposition is higher than the VV effective surface scattering component. As shown in Figure 6, when the RMS height is constant, the VV polarization backscattering coefficient simulated by the models is greater than the HH polarization backscattering coefficient simulated by the CIEM or the Dubois model. Due to the simulated backscattering coefficient having a positive correlation with the roughness parameters, when the simulated value is close to the effective scattering component value, higher effective roughness parameters are obtained using HH polarization. Among all the sampling dates, the SMC mean value on 9 May is the highest. As shown in Figure 6, when the SMC is constant, the simulated backscattering coefficient is positively correlated with the RMS height. When the RMSE between the estimated SMC and the measured SMC is the smallest, and the backscattering coefficient difference between the two dates is small, a higher RMS height is obtained when the SMC is large. Therefore, the optimal roughness parameters obtained on 9 May are higher. The average SMC on 10 July is the lowest, which is less than 10 vol.%. When the effective surface scattering component hardly changes, in order to obtain a lower SMC value, the corresponding RMS height value increases, which also indicates why the optimal roughness parameter is high on 10 July.

Taking the CIEM model as an example, Figure 13 illustrates the cost function values obtained based on different effective surface scattering components at one sampling point on 9 May, when the RMS height is 0.29 cm and 0.89 cm. F1, F2, and F3, respectively, correspond to the cost function in Equation (12), and use the corresponding effective scattering component (VV, HH, VV + HH). The horizontal axis represents the SMC. When different cost functions are used, the SMC corresponding to the minimum value is taken as the estimated soil moisture at this sampling point. As shown in Figure 13, when F1, F2, and F3 obtain minimum values, the SMC value corresponding to F3 is between the SMC corresponding to F1 and F2. Moreover, with the increase in roughness, the SMC corresponding to the minimum of the three cost functions decreases. Therefore, when the estimated soil moisture value using different cost functions is close to the actual value (F1_min ≈ F2_min ≈ F3_min), the RMS height obtained by combining the VV and HH effective scattering components is between the RMS heights obtained using a single effective scattering component. In conclusion, it proves that the optimal roughness parameter can be used as an empirical coefficient of the surface scattering models, in order to adjust the estimation accuracy of soil moisture. Therefore, in the absence of measured roughness parameters, it is feasible to use the optimal roughness parameter to parameterize the surface scattering model for the soil moisture estimation method in this study.

4.6. Soil Moisture Map

Based on these results, an optimal soil water estimation model is employed, combining the effective surface scattering components VV and HH to complete the mapping of the SMC in the study area, as shown in Figure 14. Figure 14a,b are maps using the CIEM and Dubois models as surface scattering models, respectively. Compared to the measured mean value of soil moisture, the soil moisture map agrees with the measured values. For example, Figure 12 shows that the soil moisture of the sampling points on 16 June is lower than that of the two adjacent sampling dates, seen in Figure 14. When comparing the soil moisture maps under two scattering models, the CIEM soil moisture map and the Dubois soil moisture map are generally consistent in the high- and low-value soil moisture areas. The soil moisture in the study area shows significantly high and low value distribution rather than uniform distribution, which may be due to two factors. First, the terrain of the study area is not flat and has some slope, resulting in an uneven soil moisture distribution. Second, the growth of wheat and the soil moisture absorption in different regions are not the same. Taking the Dubois-SMC map as a reference, the low-value area on the Dubois-SMC map is orange or red, and its value is lower than the same area on the CIEM-SMC map. The high-value area on the CIEM-SMC map is green in color, and its value is higher than the same area on the Dubois-SMC map. Except for the last sampling date, the patch area of the low-value area on the Dubois-SMC map is also smaller than the low-value area on the CIEM-SMC map. When combined with Figure 6 and Figure 11, we find that, except for the last sampling date, when constructing the lookup table of the surface scattering model based on the optimal roughness, the backscattering coefficient simulated by the CIEM model is not sensitive to the change in soil moisture, as compared to the Dubois model. Therefore, the comparison between the dry area and the wet area in the CIEM-SMC diagram is not apparent, as compared to the Dubois-SMC diagram. On the last sampling date, the Dubois-SMC map shows a large area of zero value. The main reason is that when the soil moisture is low, the Dubois model is not sensitive to the change in soil moisture, resulting in a significant underestimation of the SMC. Therefore, the Dubois model is not suitable for the estimation of the soil moisture over dry areas.

The range of the soil moisture is large (0–50 vol.%), which is an obvious deviation from the actual situation. This is attributed to the use of optimal roughness parameterized scattering models in this study. For each sampling date, using the optimal roughness parameter to characterize the roughness conditions of all sampling points adjusts the estimation accuracy of soil moisture as a whole, but it has limitations. For example, when the optimal roughness parameter is far less than the actual roughness parameter, the estimated value based on the optimal roughness parameter is significantly overestimated when using the same surface backscattering coefficient. An estimate of 50 vol.% above the critical value occurs. According to Table 1, the range of measured soil moisture is relatively large. The soil moisture of most pixels in the study area is accurately estimated, but the soil moisture estimated based on the optimal roughness parameters is overestimated and underestimated in some pixels. Therefore, the SMC in the figures ranges from 0 to 50 vol.%.

4.7. Limitation and Potential Improvements

Although this study provides a good estimation performance for the SMC in the study area, it has some limitations that should be improved in future research. Firstly, the measured data in the field only includes part of the wheat growth period. In future, it is necessary to collect and analyze the measured soil moisture for the entire growing season to evaluate the effective range of the proposed method. Secondly, the proposed method captures the dynamic changes of soil moisture in the study area. However, this is the result of constructing the model based on the sampling points in the study area. For the applicability of the final model, constructed based on a specific area, to estimate the soil moisture content of other wheat fields, or the same field in the same growing stage but in different years, it needs further verification in future research. Thirdly, volume scattering is the main scattering component for the wheat field in this study, but only a relatively simple volume scattering matrix is used in our study. Therefore, it is necessary to explore a complicated volume scattering matrix, or a surface scattering component extraction method, in future. In addition, this study focuses on wheat crops. Consequently, this method should also be tested on other crops, to validate its generalizability.

5. Conclusions

In this study, based on RADARSAT-2 fully polarized data and polarization decomposition, the surface backscattering coefficients with different volume scattering matrices were extracted. The SMC of the wheat fields was estimated by combining surface scattering models, optimal roughness parameters, and three cost functions. From our experiment, we drew the following conclusions:

(1)

Based on the H/α decomposition and the Freeman–Durden decomposition, the scattering components in the study area are mainly volume scattering component and surface scattering component.

(2)

The surface backscattering coefficients were extracted using four volume scattering matrix selection strategies (i.e., horizontal, vertical, random, and based on the P_r value), and combined with the optimal roughness parameters; then the CIEM and Dubois models, as well as different cost functions, were used to retrieve the SMC. Finally, we found that the VV scattering component, based on the vertical volume scattering matrix, and the HH volume scattering component, based on the horizontal volume scattering matrix, achieve the best performance. Using these two scattering components as the effective surface scattering components of VV and HH, combined with the cost function F (3), the highest estimation accuracy is obtained.

(3)

From the optimal soil moisture estimation model, the spatial distribution map of SMC in the winter wheat region was completed.

This study achieves encouraging results. Once the limitations are addressed, it has the potential to be an effective means for the remote monitoring of soil moisture in winter wheat farmland, and may be used as a basis to diagnose the health status of farmland.