Mapping Soil Surface Moisture of an Agrophytocenosis via a Neural Network Based on Synchronized Radar and Multispectral Optoelectronic Data of SENTINEL-1,2—Case Study on Test Sites in the Lower Volga Region

Abstract

In this article, the authors developed a novel method for the moisture mapping of the soil surface of agrophytocenosis using a neural network based on synchronized radar and multispectral optoelectronic data from Sentinel-1,2. The significance of this research lies in its potential to enhance precision farming practices, which are increasingly vital in addressing global agricultural challenges such as water scarcity and the need for sustainable resource management. To verify the developed method, data from two experimental plots were utilized. These plots were located on irrigated soybean crops, with the first plot situated on the right bank (plot No. 1) and the second on the left bank (plot No. 2) of the lower Volga River. Two experimental soil moisture geodatasets were created through measurements and geo-referencing points using the gravimetric method (for plot No. 1) and the proximal sensing method (for plot No. 2) employing the Soil Moisture Sensor ML3-KIT (THETAKIT, Delta). The soil moisture retrieval algorithm was based on the use of a neural network to predict the reflection coefficient of an electro-magnetic wave from the soil surface, followed by inversion into soil moisture using a dielectric model that takes into account the soil texture. The input parameter of the neural network was the ratio of the microwave radar vegetation index (calculated based on Sentinel-1 data) to the index (calculated based on the data of multispectral optoelectronic channels 8 and 11 of Sentinel-2). The retrieved soil moisture values were compared with in situ measurements, showing a determination coefficient of 0.44–0.65 and a standard deviation of 2.4–4.2% for plot No. 1 and similar metrics for plot No. 2. The conducted research laid the groundwork for developing a new technology for remote sensing of soil moisture content in agrophytocenosis, serving as a crucial component of precision farming systems and agroecology. The integration of this technology promotes sustainable agricultural practices by minimizing water consumption while maximizing crop productivity. This aligns with broader environmental goals of conserving natural resources and reducing agricultural runoff. On a larger scale, data derived from such studies can inform policy decisions related to water resource management, guiding regulations that promote efficient water use in agriculture.

1. Introduction

The soil moisture content and the characteristics of crop cover are the main parameters of irrigated agrophytocenoses, reflecting their state in terms of both the soil water regime and physiological development [1,2]. Assessing these parameters at the level of an agrophytocenosis and its components is essential for effective management and agroecological control of environmental impacts [3,4]. Simultaneously, the moisture content in the surface layer of the soil cover and the vegetation indices of the vegetation cover, measured remotely, along with the surface roughness of the soil cover of agrophytocenoses and the scattering elements of the vegetation cover, are key factors influencing the radar backscatter coefficient (RBC), measured via the radar of the Sentinel-1 satellite at a frequency of 5.4 GHz. The moisture content of the soil surface, reconstructed based on existing well-known scattering models such as Oh [5], Duboa [6], and the Fung integral equation method (IEM) [7], often leads to significant discrepancies compared to moisture values measured via proximal methods on subsatellite test plots [8,9]. This discrepancy is largely due to the variability of soil moisture and roughness, which are dynamic over time and across different scales within an agrophytocenosis. Such variability poses significant challenges in organizing periodic and highly detailed subsatellite monitoring of these parameters over large areas. Indeed, regarding the issue of restoring soil moisture according to [10], the standard deviation, correlation length, and form of the autocorrelation function of irregularity heights depend not only on the degree of surface roughness but also on the length of the measured profile (0.5–25 m). Consequently, difficulties arise when using statistical characteristics from ground-based roughness measurements as input parameters for existing radar scatter models and applying them in global algorithms for multi-scale satellite radar moisture sensing.

Alternative approaches to solving the problem of restoring soil moisture are based on developed semi-empirical methods using integral equations. These methods allow for consideration of local features in the statistical characteristics of roughness. However, this requires calibrating the input parameter (correlation length) as a function of the root-mean-square deviations of irregularity heights and the angle of incidence of the wave obtained in the study area [11,12]. For this purpose, approaches based on neural network (NN) training using scattering models are widely employed [13,14]. In these NNs, soil moisture serves as an output parameter, while combinations of backscattering cross-sections measured at different polarizations, along with soil surface roughness and sensing angle, are treated as input parameters. Despite the complexity involved in solving inverse problems with large spatial arrays of radar data, the results obtained from combining semi-empirical methods and NNs are highly promising for addressing issues related to moisture restoration [15,16].

It is noteworthy that a generalized model has yet to be developed to address radar scattering of electromagnetic waves on elements of vegetation cover. However, extinction coefficient (EC) models have already been established to describe the attenuation of electromagnetic waves in various types of vegetation [17,18]. It has been shown that EC is proportional not only to the volumetric water content in plants but also to some empirical variables, for which only approximate relationships with electromagnetic wave frequency and type of vegetation cover have been established thus far [18,19]. Recently, an empirical model has been widely used to describe the radar backscatter coefficient (RBC) of a soil surface covered with vegetation. This model accounts for the attenuation and scattering of a wave in a layer represented by uniformly dispersed particles (the “cloud” model) [20,21,22,23,24]. The parameters of this model—such as the effective scattering amplitude on particles, the effective value of the layer extinction coefficient, and the proportionality coefficients—are specific to a given plant type. These parameters are calibrated either using analytical models of scattering on vegetation cover elements [22], or, more commonly, using satellite parameters (vegetation indices; leaf area indices) alongside corresponding subsatellite measurements (biomass and vegetation height) [20,21,22,23,24]. However, to date, a generalized relationship between the parameters of the “cloud” model and various types of vegetation, depending on the sensing wave frequency, altitude, and vegetation biomass, has not been established.

Due to significant difficulties in calibrating existing scattering models for a wide variety of combinations of soil and vegetation covers, neural network (NN) methods have recently been developed to predict the moisture content of soils covered with vegetation. In this case, either radar data (RBC at different polarizations) [15,25] or combinations with multispectral optical data (vegetation indices; leaf surface indices) measured via various survey systems on space platforms are used as input parameters for the NN [26,27,28]. Soil moisture is now taken as the direct output parameter of the NN data. At the same time, no additional connections are established between the combination of multispectral and radar indices and the biophysical parameters of vegetation [28,29,30], which would increase the unambiguity of the problem and possibly simplify the structure and number of layers (neurons) in an NN.

In contrast to these approaches, this work utilizes the NN in two stages. In the first stage, the NN predicts the reflection coefficient of the soil cover. For this purpose, the ratio of the multispectral index of vegetation cover in the optical range to the microwave index of the same cover is used as an input parameter. In the second stage, based on the dielectric soil model [31], soil cover moisture is restored by solving the inverse problem using the reflection coefficient estimated via the NN. The proposed algorithm helps minimize the influence of vegetation cover on the restored value of soil cover moisture.

2. Materials and Methods

To test the proposed method, data sets from two test plots located in the southern part of the European region of Russia were utilized. Test plot No. 1 (44.1568 N, 48.6026 E) was situated in one of the fields of the experimental production facility of the All-Russian Research Institute (VNIIOZ) (see Figure 1a). During field measurements in July 2020, the test plot was sown with soybeans in the heading phase and irrigated with a Bayer frontal sprinkler machine. Test plot No. 2 (51.1258 N, 46.0001 E) was located on one of the agricultural fields of the Educational Research and Production Association “Povolzhie” at Saratov State Vavilov Agrarian University (SSAU) (see Figure 1b). During the test period in the second half of August 2022, the corresponding field containing this test plot was sown with soybeans in the ripening phase, which was irrigated with a “Kascad” circular sprinkler.

The soil cover of test plot No. 1 consisted of light chestnut irrigated alkaline calcareous soils on yellow-brown loams (according to the soil classification of the USSR). Soils are classified according to the WRB classification as Luvic Kastanozems (Loamic, Aric, Protosodic, and Bathygypsic). The texture of the soil cover was determined using a combination of sieve and pipette (Kaczynski version) analyses. The obtained results indicated that the soil particle distribution was classified as silty loam (clay loam) according to FAO soil classification.

The soil cover of test plot No. 2 was characterized by a complex of medium and thin dark chestnut soils with medium loamy and light loamy granulometric composition. The content of physical clay in the plow horizon ranged from 36% to 38%. Additionally, the volumetric mass of the topsoil was measured at 1.34 g/cm³, while the density of the solid phase ranged between 2.62 and 2.65 g/cm³. The corresponding value of porosity in the arable layer lies between 0.49 and 0.53 cm³/cm³, and the value of field capacity (FC) is between 0.25 and 0.27 cm³/cm³.

At test plot No. 1, satellite monitoring included an aerial survey of the moisture in a 0–5 cm layer of soil, as well as measurements of the height of soybean plants. Both types of measurements were conducted synchronously and corresponded to the times of the Sentinel-1 radar survey on two dates, 9 July and 21 July 2020. Simultaneously, subsatellite moisture monitoring involved collecting undisturbed soil samples using a special sampler, followed by determining their volumetric moisture content in laboratory conditions [32]. On both specified dates, samples were collected at 45 points in the nodes of a uniform rectangular grid, with a distance of about 10 m between nodes and covering an area of approximately 0.6 ha (see Figure 2a).

At test plot No. 2, subsatellite monitoring included a one-time measurement of volumetric soil moisture and soybean plant height at the time of the Sentinel-1 radar survey on 22 August 2022. These measurements were carried out at 201 points within one of the irrigated soybean fields, where moisture content had been established through irrigation with the “Kaskad” system on 21 August 2022. The geometric boundaries of this section began at the center of rotation of the “Kaskad” and extended northwestward, protruding 25 m beyond the outer boundary. Within this section, two parallel routes were laid out, with transverse distances between them ranging from 10 to 15 m. To measure the volumetric moisture content of the soil, an ML3-KIT (THETAKIT) device manufactured by Delta was used. The device was tested prior to use [33]. The average deviations of the moisture content obtained with this device from those derived from soil samples conformed to its technical specifications and were within 2%. Moisture monitoring in the second area was conducted approximately 24 h after it was irrigated with the “Kaskad” at a rate of 17 mm.

To create a geodatabase of the measured values for soil moisture and plant heights corresponding to both test plots, the freely distributed mobile software application GPS MapCamera 1.4.17 was utilized. This application is available for smartphones running iOS and Android platforms. It enabled the collection of necessary video images, which were then used to form the corresponding layers of the geodatabase, including date and time of measurements, coordinates of sampling/measurement locations, plant height at those points, and either the sampling box number or the measured moisture value [30]. Volumetric soil moisture values measured in test plot No. 1 ranged from 6% to 26% (9 July 2020) and from 11% to 23% (21 July 2020), while plant heights varied from 55 cm to 80 cm (9 July) and from 70 cm to 110 cm (21 July). In test plot No. 2, volumetric soil moisture measured on 23 August 2022, varied from 5.2% to 36.1%, with plant heights ranging from 85 cm to 100 cm.

Maps of the interpolation results of the measured values of soil moisture, with the places of the measurements plotted on them, as well as maps of NDVI = (К₈ − К₄)/(К₈ + К₄), where К_4,8 is the corresponding Sentinel-2 channels, calculated based on the results of the Sentinel-2 survey, of both test plots are shown in Figure 2. The NDVI (Sentinel-2) maps for 11 July 2020 and 22 August 2022 are used as a background in Figure 2c (test plot No. 1) and Figure 2d (test plot No. 2), respectively. The local locations of soil sampling for measuring soil density and soil moisture are shown as red and black dots in Figure 2c (test plot No. 1) and white dots in Figure 2d (test plot No. 2). In addition, circles (see dashed lines) are added to Figure 2c, which correspond to the support points of the irrigation machine.

The variation in NDVI indices for both plots fell within close limits of 0.4–0.8 (see Figure 2). The Sentinel-1 satellite measured radar backscatter coefficients (RBCs) in interferometric wideband mode (IW) at a frequency of 5.4 GHz with VH and VV polarizations over the territories of both test plots (9 July and 21 July 2020 and 22 August 2022, at 7:06 local time, UTC+4). Using ESA SNAP software 10.0.0, standard processing of Sentinel-1 data was performed: this included utilizing precision orbits, calibration, and speckle noise filtering through successive application of two Gamma map filters sized at 3 × 3 pixels. After filtering of Sentinel-1 data, the effective spatial resolution became about ~10 × 10 m. Radiometric accuracy is about 1 dB (3σ).

The RBC was normalized to a reference probing angle of 30° for both test plots according to the method outlined in [25,34]. Sentinel-2 (MSIL2A) multispectral survey data were collected on 11 July and 21 July 2020, for test plot No. 1, and on 22 August 2022 (11:56 local time), for test plot No. 2. For consistent processing, the Sentinel-2 multispectral survey data were recalculated using inverse distance-weighted interpolation on the Sentinel-1 radar data grid. (Note that the resolution of Sentinel-2 images was reduced to the resolution of channel 11–20 m).

In the study, the created NN model was calibrated and verified according to the results obtained in the first test section, and its additional verification prior to training of the NN was carried out according to the results obtained in the second test section. The RBC, as well as the NDVI (normalized difference vegetation index) values were calculated, respectively, based on data from Sentinel-1,2 satellites covering test plot No. 1 (see Figure 3). The Pearson correlation coefficient between the RBC at vertical–vertical (σ_VV) and vertical–horizontal (σ_VH) polarizations and soil volumetric moisture is no more than 0.227 and 0.084, respectively (see Figure 3a). The Pearson correlation coefficient between the NDVI (Sentinel-2) and vegetation height was 0.297 (see Figure 3b). Due to the fact that the RBC is more susceptible to volume scattering by vegetation cover elements at cross polarization (σ_VH), the correlation between the RBC at matched polarization σ_VV and soil volumetric moisture is stronger for cross polarization (see Figure 3a).

The weak correlation between the NDVI and vegetation height, h_l, (see Figure 3b) is apparently due to the fact that the NDVI is more related to the reflective characteristics of the vegetation cover, which depends on its photosynthetic activity, rather than on the general volume of biomass associated with the height of plants. At the same time, we found a significantly greater correlation between the vegetation index (Sentinel-2) I₀ = (К₈ − К₁₁)/(К₈ + К₁₁) as well as the radar vegetation index RVI = 4σ_VH/(σ_VH + σ_VV) and plant height in test plot No. 1 (Sentinel-1) (see Figure 4).

The Pearson correlation coefficient (0.798) between the multispectral vegetation index I0 and plant height hv is higher than that between the microwave index RVI and hv (0.334). This is due to the fact that the multispectral index I0 contains information about the interaction of reflected solar radiation with the surface of vegetation cover, while the microwave index RVI provides information about the interaction of an electromagnetic wave with the reflectivity characteristics of both the vegetation cover and the soil surface (see Figure 5).

Figure 6 shows a simple feed-forward NN with one hidden layer containing N neurons.

The input parameter of the used NN (see Figure 6) is the ratio ξⁱⁿ = RVI/I₀. In contrast to the existing approaches, instead of soil moisture, the output parameter was represented by the modulus of the Fresnel reflection coefficient of an electromagnetic wave with a flat front from the soil surface with a smooth boundary ξ^out = |R₀(ε_s)|, where ε_s is the complex dielectric permittivity (CDP) of soil. Using |R₀(ε_s)| as the output value of NN, we do not have to train it every time for a new type of soil cover, but to use a dielectric model that takes into account the dependence of the CDP on the type of soil cover ε_s = ε_s(W, m_c) [5,24,31]; here, W and m_c are the volumetric soil moisture and the clay fraction share in the soil cover.

NN was modeled by means of Matlab (see function «feedforwardnet»). A feed forward NN was used, consisting of one hidden layer, in which from 1 to 65 neurons were specified. A hyperbolic tangent was used as neural transfer function for hidden layers. A linear transfer function was used for output layer. The root-mean-square deviations between the output true values and the output values predicted via the NN during training were minimized based on the Levenberg–Maquard algorithm. The result of NN training depending on the number of neurons is shown in Figure 7. When calculating the true values of the reflection coefficients, we used the dielectric model [24] and data from ground-based measurements of the volumetric soil moisture at the soil sampling points in test plot No. 1 (see Figure 2a).

The coefficient of determination, R², and the standard deviation (RMS) between the predicted NN model ${|R}_0^N|$ and calculated |R₀(ε_s)| values of the modulus of the reflection coefficient varies from R² = 0.31 and RMS = 0.039 to R² = 0.63 and RMS = 0.05 as the number of neurons increase from N = 1 to N = 65. Due to the fact that with an increase in the number of neurons the values of R² and RMS are more and more random, for further calculations, the number of neurons in the hidden layer was set equal to N = 20 Subsequently, volumetric soil moisture, $W_{retr}$, can be determined by minimizing the norm of the discrepancy between the informative features of the estimated reflection coefficient $\left|R\right( {\mathsf{\epsilon }}_s(W_{.retr}, m_{с})\left)\right|$ and the value of ${|R}_0^N|$, predicted via the NN model based on the observational data of the Sentinel-1,2 satellites.

$$W_{retr}=min F\left(W_{retr}\right), F\left(W_{retr}\right)={\sum }_{n=1}^{n=N}\left|{\displaystyle \frac{{|R}_0^N\left|-\left|R\left({\mathsf{\epsilon }}_s\left(W_{retr},m_{с}\right)\right)\right|\right)}{{|R}_0^N|}}\right|,$$

(1)

The minimization task in (1) was solved via a direct method by selecting W_retr from the range of W_retr ∈ [0%, 50%] with a step of 1% for the central coordinate of each pixel.

3. Results and Discussion

3.1. Results of Test Plot No. 1

Soil moisture values restored from the combined radar and optical data of the Sentinel-1 and Sentinel-2 satellites, relative to the moisture values measured on 9 July and 21 July 2020 at the locations where soil samples were taken on test plot No. 1 are shown in Figure 8.

With a coefficient of determination and RMS equal to 0.435 and 2.4%, respectively, the reconstructed values of soil moisture from the remote sensing results align with the soil moisture measured using the sampling method in the 0 to 5 cm layer under the vegetation cover on test plot No. 1. Figure 9a,b show, as an example, maps of soil surface moisture built based on the proposed method using Sentinel-1,2 satellite data, as well as similar maps built based on the results of measurements on 9 July 2020 at test plot No. 1. The maps in Figure 9a,b were constructed for each pixel of the Sentinel-1 observation (size ~10 × 10 m) in test plot No. 1. Figure 9c,d show a map of the absolute error of the retrieved soil moisture values from Sentinel-1 data using a neural network relative to the measured values in situ (using the thermostat-weighting method).

The comparison of the maps presented in Figure 9 demonstrates a fairly good agreement between the spatial variations in the measured and predicted soil moisture values in different local areas of the field. In this case, the standard deviation between the reconstructed and determined soil moisture values is approximately 2.4%, with maximum and minimum absolute errors of +5.5% and −3.1%, respectively (see Figure 9c,d).

The maximum and minimum absolute errors of the reconstructed soil moisture values relative to the measurements taken on 21 July 2022 across the entire area of test plot No. 1 were +1.9% and −2.7%, respectively. The maps presented in Figure 9 were created using inverse-distance-weighted interpolation of various-scale radar and ground data over a 9 × 17 section within a rectangular area measuring 116 m × 58 m.

3.2. Results of Test Plot No. 2

For soil moisture predictions in test plot No. 2, the pre-trained neural network (NN) model obtained for test plot No. 1 was utilized (see Figure 6). A sample containing 50% of the subsatellite soil moisture content for test plot No. 2, obtained on 22 August 2022, was used for this purpose. In accordance with the proposed method (see Paragraph 3), the modulus of the reflection coefficient was predicted using the RVI and I0 values measured via the Sentinel-1 and Sentinel-2 satellites through the pre-trained NN. Subsequently, using the dielectric model [24], we solved the inverse problem (1) to restore soil moisture. Additionally, to assess the significance of each of the I0 and RVI parameters in the overall soil moisture prediction via the NN, it was also pre-trained based on inputs from either 50% of just the I0 or RVI. The results of soil moisture restoration in test plot No. 2 are presented in the cartograms (see Figure 10). Maps (see Figure 10) of the retrieved soil moisture values from Sentinel-1 data using a neural network were constructed for each pixel of the Sentinel-1 observation (size ~10 × 10 m) in test plot No. 2. The dotted circular lines in Figure 10 indicate the lines of movement of the wheels of the watering machine. An analysis of these maps shows that both maps—one constructed using the I0 index (see Figure 10a) and the other using the RVI (see Figure 10b)—are similar. However, the soil moisture map created using the I0 index (see Figure 10a) appears to be more detailed, with clear contours and well-defined sprinkler stops.

In contrast, the RVI radar image obtained through antenna aperture synthesis exhibits interference speckle noise, which accounts for the noisiness of the soil moisture map (see Figure 10b). Furthermore, after filtering the speckle noise, the effective resolution of the radar image decreased by a factor of three to 30 m × 30 m.

From the comparison of soil moisture maps (see Figure 10), it is evident that the I0 index (see Figure 10a) has a greater corrective effect on the soil moisture map (see Figure 10c) constructed using the RVI/I0 index ratio. At the same time, it can be observed that the soil moisture values restored based on the I0 index are somewhat “shielded” and fairly averaged (displaying similar colors). Conversely, the map constructed based on the RVI/I0 index ratio is much more contrasting and reflects a broader range of soil moisture variations.

(Note that the Sentinel-1 survey at 07:06 on 23 August 2022 was conducted 19 h after the Sentinel-2 survey at 11:56 on 22 August 2022 and that the Sentinel-2 survey occurred four hours after an irrigation event at test plot No. 2.) In this case, concentric ring patterns with higher moisture values were at the center and along the outer boundary, while lower values were recorded in the middle. This phenomenon was attributed to the varying irrigation rates applied via the “Kaskad” irrigation machine along the pressure front [35]. Indeed, when comparing the retrieved soil moisture values through three methods with in situ measurements, it is evident that the highest determination coefficient is observed when using the RVI/I0 index as an input parameter for the NN (see Figure 11).

4. Conclusions

This study demonstrates the promising potential of using neural networks (NNs) to establish adaptive relationships between multispectral and microwave indices, as well as the reflective properties of vegetated soil cover, in order to restore soil moisture. In this context, there is no need to calibrate the parameters of scattering models based on the height (biomass) of the vegetation cover obtained through ground-based measurements. By applying this combined approach for training the NN, only ground-based soil moisture measurements are required, in addition to Sentinel-1 satellite radar polarimetric observations (radar backscatter cross-sections on VV and VH polarizations) and multispectral measurements from the Sentinel-2 satellite (channels 8–11).

One drawback of the proposed approach is that it does not account for the probing angle, an effect that needs to be studied in detail in future research on test plots located at significant distances from one another. Additionally, it is essential to evaluate the efficiency of the proposed method in minimizing the impact of vegetation cover by utilizing the ratio of multispectral optical and microwave plant indices for various types of irrigated crops growing on soils with different granulometric compositions and organic content. It is also worth noting that Sentinel-2 optical data are not always available for the area of interest due to dense cloud cover.

The results indicated that when restoring soil moisture in two test plots, the combined use of RVI and I0 proved to be more informative, as their ratio resulted in a more linear relationship with the restored soil moisture values compared to the measured in situ values. The developed method allowed for high reliability in identifying patterns of high and low soil moisture values within irrigation boundaries created with sprinkler systems, as well as patterns formed due to irrigation runoff beyond these boundaries, indicating potential risks of negative environmental consequences.

Of particular interest is also the investigation of other types and architectures of neural networks to improve the proposed method. This proposed method can be utilized to monitor the progress and outcomes of irrigation practices implemented using various technologies, including spatially differentiated irrigation techniques. However, the proposed methodology in this work requires further validation on various agricultural crops at different stages of vegetation growing in different regions.