Using Sentinel-1 Time Series Data for the Delineation of Management Zones

1. Introduction

Information regarding the spatial and temporal variability of soil attributes plays a crucial role in the development of effective soil management strategies. By examining these data, farmers can adopt the most suitable cultivars and plant population densities for each specific point within the production area. This, in turn, facilitates the precise determination of the required amounts of fertilizers and soil acidity correctives, not only to maximize financial returns but also to promote more sustainable production [1].

However, developing an effective strategy for collecting data to characterize the spatial and temporal variability of soil attributes is a complex and challenging task. Research has highlighted the importance of establishing dense sampling grids, with a minimum density of one sample per hectare, to adequately capture the variability of soil attributes [2,3,4]. [5] demonstrated that variograms used to infer soil attributes at unsampled points are unreliable when based on fewer than 100 data points, potentially leading to inaccurate estimates with significant margins of error. Therefore, grid sampling can provide a precise basis for variable rate application, but the costs and labour requirements, especially in extensive areas with high variability, suggest that other approaches may be more economical [6].

To understand the spatial and temporal variability of soil attributes without the need to establish dense sampling grids, studies have demonstrated the potential of using soil sensors or the crop itself as a 'soil sensor' [1]. Apparent soil electrical conductivity sensors, yield maps, and canopy reflectance indices can provide maps with different spatial and temporal variability patterns and be used to delineate homogeneous areas known as Management Zones (MZs) [7,8,9,10,11,12]. Within each MZs, low variability of soil attributes is assumed, recommending the collection of a single composite sample. Based on specific levels of these attributes, targeted management practices for each MZs are established. This strategy reduces soil sampling costs compared to dense sampling grids, while simultaneously providing a better distribution of management practices (cultivars, plant density, fertilizers) compared to the conventional soil sampling method, in which only a single attribute level and consequently a single management strategy is determined for the entire area.

Although the development of MZs through these methods represents an advancement in precision agriculture, their adoption among farmers remains limited. This limitation is largely due to difficulties in accessing reliable historical yield maps, electrical conductivity data, and multispectral satellite image time series with high temporal resolution. For example, yield maps have been available since the early 1990s, yet their adoption is still limited to only 5% to 25% of the total cultivated area in the United States for crops such as winter wheat, cotton, sorghum, and rice, and 45% for corn and soybean crops [13]. Apparent soil electrical conductivity presents itself as an attractive alternative because it can be quickly and easily measured for fields using electromagnetic induction instruments. However, this type of data collection strongly depends on specialized service providers for data acquisition and interpretation, whose availability varies across agricultural regions, complicating the implementation of this technology.

The use of multispectral optical images, freely available from orbital platforms such as Landsat-8 and Sentinel-2, enables remote service delivery and extensive spatial coverage. However, its application faces significant challenges, such as cloud cover, which compromises consistent data acquisition. This issue is particularly critical in tropical regions, where average annual cloud cover can reach approximately 66%, hindering the construction of representative historical time series [14,15]. Therefore, to expand farmers' adoption of MZs, it is essential to develop alternative methods capable of efficiently characterizing the spatial and temporal variability of soil attributes, with lower cost per unit area and broader spatial coverage.

A promising line of research for characterizing the spatial and temporal variability of soil attributes through MZs is the use of Synthetic Aperture Radar (SAR) data. The Sentinel-1 mission, part of the European Union's Copernicus program, currently consisting of the Sentinel-1A sensor, freely provides SAR imagery with a spatial resolution of 20 × 22 meters and a temporal resolution of 12 days [16]. Equipped with an active C-band SAR sensor operating at a central frequency of 5.405 GHz with dual polarization (Vertical-Vertical and Vertical-Horizontal), this satellite can penetrate cloud cover and acquire imagery both day and night [17,18,19,20]. Moreover, its electromagnetic waves, characterized by a longer wavelength, can penetrate the superficial vegetation layers and, in some cases, reach deeper soil layers. In agricultural contexts, SAR backscatter data have been used, either alone or in combination with multispectral data, for various applications, including soil moisture estimation [21,22,23], assessment of soil physical properties [24,25,26,27], and estimation of multispectral indices such as the Normalized Difference Vegetation Index (NDVI) [28,29,30], among other applications.

Therefore, the previously mentioned properties highlight the potential of Sentinel-1 SAR imagery as a rich source of spatiotemporal information, making it promising for mapping soil attributes and delineating MZs. A methodology can be applied to create temporal profiles of backscatter with dual polarization—VV (vertical-vertical) and VH (vertical-horizontal)—from SAR data, complemented by the calculation of specific SAR indices. These temporal profiles can be analyzed using unsupervised classification techniques to identify regions with similar backscatter responses, potentially associated with variations in soil attributes. Thus, the objective of this study was to develop a method for mapping soil attributes through the delineation of MZs using SAR data provided by Sentinel-1.

2. Materials and Methods

This study was conducted in two commercial grain production fields (Field A and Field B) that exhibit different soil texture characteristics (Figure 1). Field A covers an area of 117 hectares and is situated in the municipality of Sinop, Mato Grosso, Brazil (11°8'20" S and 56°19'18" W). Field B spans an area of 106 hectares and is situated in the municipality of Chapadão do Céu, Goiás, Brazil (18°20'10" S and 52°37'12" W). According to the Brazilian Soil Classification System, Field A is identified as a Dystrophic Red-Yellow Latosol, while Field B is classified as a Dystrophic Red Latosol [31]. For soil sampling, a grid-point sampling method was adopted with an approximate spacing of 100 meters. Field A was represented by 113 samples, and Field B by 104 samples. The values of 9 soil attributes were determined in the laboratory from the grid of points established in each field. Descriptive statistics of these attributes are summarized in Table 1.

The SAR data used in this study were freely obtained from the Sentinel-1A sensor of the European Union's Copernicus program. The Sentinel-1 mission provides global SAR data in the C-band (central frequency of 5.405 GHz) with dual polarization (VV and VH). The temporal resolution is 12 days, although it can be higher in some cases due to overlapping sensor passes. In this study, the Sentinel-1 collection available on Google Earth Engine was used, comprising Ground Range Detected (GRD) format images processed with the Sentinel-1 toolbox to produce calibrated and orthorectified products. All images were acquired in descending orbits using the Interferometric Wide (IW) swath mode and dual polarization (VV and VH). They have a pixel spacing of 10 meters but a spatial resolution of 20 × 22 meters [16].

The preprocessing steps for the SAR data included border noise removal, speckle filtering, terrain radiometric normalization, and conversion of the backscatter coefficient to decibels. Image border noise results from the process of converting acquisitions from GRD format to IW, and its presence is an undesired processing artifact that limits its full exploitation in various applications [32]. The speckle phenomenon, common in SAR images due to the interference of radar waves reflected by surfaces smaller than the radar resolution, was addressed through multitemporal filtering [33]. The Refined Lee filter, 3x3 [34], was used with a multitemporal filtering structure of 10 images. Terrain radiometric normalization corrects variations in the received signal due to terrain slope. For this, the Shuttle Radar Topography Mission digital elevation model with a 1-arc-second resolution (~30 m) [35] was employed, deriving elevation, slope, and aspect values for normalization. Finally, as the last preprocessing step, the terrain-corrected radiometric backscatter coefficient is converted to decibels through a logarithmic transformation. Table 2 summarizes the parameters and specifications for image acquisition and preprocessing. The data were preprocessed using Google Earth Engine (GEE) [36]. All GEE codes for Sentinel-1 data preprocessing were provided by [37] and are available at https://github.com/adugnag/gee_s1_ard.

After preprocessing, each SAR image was converted into a dataframe and stored to generate temporal backscatter profiles. For each analyzed field, a grid of 40m x 40m quadrilaterals was established. In each image, the average pixel values within these quadrilaterals were calculated, thereby constructing the temporal backscatter profiles for the VV and VH bands. This method helps minimize potential residual noise in the images and reduces the computational load required for subsequent analyses. Based on these profiles, six indices were calculated using Sentinel-1 data, as detailed in Table 3. Each VV and VH temporal backscatter series, as well as the indices for each quadrilateral, were standardized using z-score normalization. This process adjusts the data so that each set has a mean of 0 and a standard deviation of 1. To avoid issues of collinearity among the indices, Pearson correlation (r) was calculated between them. Based on this analysis, only those indices that showed lower collinearity (r < 0.95) were included, thus ensuring the independence and relevance of each chosen index for the definition of MZs.

To understand the seasonal variability of the SAR indices, rainfall data from the NASA-POWER system (https://power.larc.nasa.gov) were used, considering that radar data are sensitive to soil moisture [21,22,23].This system was developed to provide meteorological information directly applicable to fields such as architecture, energy generation, and agrometeorology. It compiles information from various data sources, including grid-derived data, to offer a comprehensive view of climate and weather conditions (Maldonado et al. 2019).

To delineate the MZs, all SAR indices were analyzed individually, and two approaches were proposed. The first approach (Approach 1) involved the direct application of the fuzzy k-means clustering algorithm (Bezdek et al. 1984) on the temporal series of the SAR indices. In the second approach (Approach 2), a machine learning method known as autoencoders was implemented to reduce the dimensionality of the temporal series. Autoencoders are a type of neural network often used in unsupervised machine learning tasks, such as feature extraction (Hoang and Kang 2019). The basic architecture of an autoencoder is divided into three parts: encoder, bottleneck layer, and decoder. The encoder receives the input data (in this case, SAR indices) and compresses it into a lower-dimensional representation (bottleneck layer). The decoder then takes this compressed representation and attempts to reconstruct the original data from it. This process is carried out during the training of the autoencoder network. Thus, after training, the bottleneck layer was used as input for the fuzzy k-means algorithm to cluster the SAR data time series and, consequently, define the MZs. To define the architecture and training parameters, a k-fold cross-validation procedure, with k=5, was implemented aiming to identify the most suitable hyperparameters. The selection of the final model was based on the lowest mean squared error obtained during the validation process. The selected architecture and training parameters are detailed in Table 4 and Table 5.

The simulations were conducted considering the number of clusters, in this case MZs, equal to three for both approaches and fields. After clustering, QGIS geoprocessing tools were used for refinements. Clusters with an area smaller than 3 hectares were integrated into the larger contiguous clusters, while clusters larger than 3 hectares but geographically disconnected and sharing the same label were considered as distinct clusters.

The most common soil sampling methods include the conventional soil sampling method (CONV), the cell-based soil sampling method (CEL), the uniform grid soil sampling method (GRID-1), and the MZs-based soil sampling method. The CONV method involves collecting several samples to form a single composite sample, which is considered representative of the entire field. In the present study, the CONV method was considered as the average of all collected samples. The CEL method, on the other hand, divided the area into cells (polygons of 5 to 10 hectares), ensuring the presence of at least four GRID-1 samples in each cell. The attribute estimation in each cell was the average of the samples collected within each cell. A fifth comparison method is proposed where cells are created randomly, called random cell-based soil sampling method (CEL-RND) (Figure 2). The number of random cells was defined to match the number of MZs established by the SAR data-based method. In total, 1000 random cell scenarios were generated for each study area. The minimum size of each random cell was 4 hectares, to ensure the presence of at least four soil samples per cell. The random cells were generated using a Python script based on Voronoi Diagrams, proposed by Georgy Voronoi (Voronoi 1908), and the soil attributes within each cell were estimated by the average of the samples collected within each cell.

To evaluate the performance of each sampling method, the 'leave-one-out' cross-validation (LOOCV) method was used, always based on the GRID-1 data. In this way, we ensured the participation of all points in the error calculation for all evaluated methods. In LOOCV, each point from the GRID-1 dataset was successively removed, and its estimate was made using the sampling methods: CONV, CEL, MZs, GRID-1, and CEL-RND. After the estimation, the point was reintegrated into the dataset. This process continued until all GRID-1 points were evaluated. For the CONV method, the estimate of the removed point was calculated from the average of the remaining points. For the CEL, MZs, and CEL-RND methods, the estimate was based on the average of the remaining samples within their respective areas. Finally, to evaluate the GRID-1 method, the estimate of the removed point was a value interpolated using ordinary kriging based on the fitting of semivariograms. Semivariograms are tools that allow for the characterization and determination of distribution patterns, such as randomness, uniformity, and spatial trends [41]. For this, equation (1) was used to calculate the semi-variance:

$${\gamma }^{*}\left(h\right)={\displaystyle \frac{1}{2N\left(h\right)}}{\sum }_{i=1}^{N\left(h\right)}{(z\left(x_i\right)-z(x_i+h\left)\right)}^2$$

(1)

where ${\gamma }^{*}\left(h\right)$ is the value of the experimental semivariance at the distance interval h; $z\left(x_i\right)$ is the sample value measured at the sampling points $x_i$, where data exist at $x_i$ and $x_i+h$; and $N\left(h\right)$is the total number of sample pairs within the distance interval h.

During the adaptation of theoretical models to the experimental semivariograms, coefficients were determined that describe the nugget effect (C0), sill (C0 + C), partial sill (C), and range (A). The models tested for adaptation included the spherical, exponential, Gaussian, and linear models, and they were selected based on maximizing coefficient of determination (R²), minimizing sum of squared residuals, and maximizing the correlation coefficient obtained through cross-validation. These metrics are used to evaluate how well the fitted model matches the experimental data. Spatial Dependence Index (SDI) was analyzed using the ratio C0/(C0 + C), and the intervals proposed by [42] were employed to classify spatial dependence into three categories: strong dependence (SDI < 25%), moderate dependence (25% ≤ SDI < 75%), and weak dependence (SDI ≥ 75%). Semivariograms of soil chemical and physical attributes were modelled using SmartMap [43] Version 1.4, an open-source plugin developed for QGIS.

By comparing the estimated values of soil attributes obtained by the CONV, CEL, MZs, GRID-1, and CEL-RND methods with the sampled values at each corresponding point, the Root Mean Square Error (RMSE) was calculated. This analysis was conducted individually for each study field, following the methodology established in equation (2):

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

(2)

where: $\widehat{x_i}$ represents the estimated value of the soil attribute at point i; $x_i$ is the observed value of the soil attribute at point i and n is the number of sampled points.

To compare the MZs defined from SAR data with the other soil sampling methods used in this study, the combination of approach (1 or 2) and SAR index that resulted in the lowest RMSE was selected. After this selection, the average of each attribute was calculated for each zone, and then the attributes were compared using the F-test (with a significance level of p-value < 0.05). The RMSE values of the MZs sampling method for each field were compared with the other values resulting from the sampling methods evaluated in this study. To evaluate the CEL-RND method, we quantified, across 1000 generated scenarios, the frequency with which the MZs method showed a higher RMSE than CEL-RND. Then, the percentage of these scenarios in which the MZs method performed worse, in terms of RMSE, compared to the CEL-RND method was calculated.

3. Results

3.1. Exploratory Analysis of the SAR Dataset

Between January 1, 2018, and March 31, 2023, 319 images were obtained for Field A and 266 images for Field B, resulting in average temporal resolutions of 6.0 and 7.4 days, respectively. Figure 3 shows the Peason correlation estimated between the VV and VH backscatter values and the SAR indices. The correlations obtained were found to be significant, with p-values equal to zero. The indices RVI, NRPB, VH/VV, and VV/VH show high r with each other, as well as PRVI and VH, exhibiting values above 0.95. On the other hand, the RVI4SI index showed the lowest level of r with the other indices. Based on these results, four indices (VV, VH, VV/VH, RVI4SI) were selected for the continuation of the study.

The monthly averages of the time series for the VV and VH backscatter coefficients, as well as the VV/VH and RVI4SI indices, were plotted for both fields in each year (Figure 4). Considering that the data availability extended only until March 2023, the graphical analysis was restricted to the period from 2018 to 2022. The time series of the calculated indices exhibited seasonal trends. Except for the VV/VH index, an increase was observed between September and December, followed by a decline in February. This behavior coincides with the period of increased monthly accumulated precipitation (Figure 5), suggesting a possible relationship between higher soil moisture or vegetation growth and elevated backscatter levels. From February to March, although precipitation tends to remain relatively constant, there is a new increase in backscatter, likely due to the growth of second-crop vegetation. Subsequently, from May to September, there was a decrease in the SAR indices, in line with the reduction in monthly accumulated precipitation and crop biomass. For the VV/VH index, there was an increase between January and February, which then gave way to a decrease from February to April, except for the years 2019 and 2021. Additionally, only for the VV/VH index, an increase in values was recorded from April to August.

3.2. Analysis of Spatial Variability of Soil

The analysis of the experimental semivariograms in GRID-1 confirmed the spatial variability of soil attributes for both fields (Table 6). When analyzing the SDI, it was found that clay (CLA) content and soil organic carbon (C) in Field A and CLA and potential acidity (H+ + Al3+) in Field B had SDI values lower than 25%, indicating a high spatial dependence [42]. All other elements, except for phosphorus (P) in both fields, exhibited SDI values between 25% and 75%, which indicates moderate spatial dependence. Figure 6 and Figure 7 show the maps constructed using ordinary kriging after semivariogram fitting. In Fields A and B, the maps of CLA and C, as well as those of Ca²⁺ and magnesium (Mg²⁺), respectively, display visual similarities that indicate high correlations between these soil attributes.

3.3. Delineation of Management Zones with SAR

Regarding the design and quantity of MZs, visual variations were observed depending on the SAR index used and the methodology adopted (Figure 8). These variations in size, shape, and number of MZs become even more noticeable when contrasting the proposed methodologies. The total number of MZs in Approach 2 exceeded that of Approach 1. In the case of Field B, the number of MZs according to Approach 1 was 2, lower than the initially stipulated value of 3. This result is because some clusters created by the Fuzzy C-means algorithm did not have significant associations with data points or exhibited an extremely low degree of membership in relation to all points for a specific cluster. Additionally, the RVI4SI index resulted in excessively fragmented clusters compared to the other algorithms evaluated. Due to this fragmentation, it was considered inappropriate to use this index for the creation of MZs in Field B. Therefore, its results were not considered in this study.

Table 7 presents the RMSE obtained through LOOCV, representing the accuracy of soil attribute estimates for Fields A and B using Approaches 1 and 2 with the MZs method. In the evaluation of errors associated with soil attribute estimation, the VV/VH index consistently stood out, achieving the lowest RMSE values compared to the other indices. For example, in Field A, when estimating the clay content attribute using Approach 2, the VV backscatter band recorded an RMSE of 8.33, while the VV/VH index showed a significantly lower RMSE of 5.74. When analyzing the two approaches, it is observed that Approach 2 has an advantage in terms of accuracy over Approach 1, especially when adopting the VV/VH ratio for soil attribute estimation. In all fields studied, the attributes CLA, P, Ca^2+, and C recorded lower RMSE with Approach 2 when using the VV/VH index.

The results indicated significant differences in soil attribute values between the MZs delineated by Approach 2 and the SAR index VV/VH, which had the lowest RMSE (Table 8). At least one mean of each soil attribute from the MZs differed statistically at a p-value < 0.05, except for the V attribute in both fields. For Field A, MZ₃ was characterized as the zone with the highest average values of CLA, potassium (K⁺), Ca²⁺, Mg²⁺, C, and potential acidity (H⁺ + Al³⁺), while MZ₂ had the lowest values. For Field B, MZ₂ was characterized as the zone with the lowest average values of CLA, P, Ca²⁺, Mg²⁺, and C.

Figure 9 presents a comparison between the error (RMSE) obtained by the MZs method using Approach 2 and the SAR index VV/VH and the error generated from 1000 different CEL-RND scenarios. The value indicated on each bar represents the percentage of times the MZs method, based on SAR data, had a higher error than the random cell method. Therefore, the lower the percentage value displayed on the bar, the better the performance of the MZs method compared to CEL-RND. Values equal to or greater than 50% indicate that MZs do not contribute to the representativeness of the soil attribute variability compared to a random process. As can be observed, for most soil attributes in both fields, the RMSE resulting from the CEL-RND sampling method was higher than the RMSE obtained by the MZs sampling method, resulting in values below 50%. In Field A, the attributes with the lowest percentage of errors (less than 5%) are CLA, K⁺, C, H⁺ + Al³⁺, and pH. This indicates that the MZs method performed well for these attributes in this field. The V attribute, on the other hand, showed a significantly higher error percentage of 43%, indicating that the MZs method may not be as efficient for this specific attribute. In Field B, the CLA attribute had a relatively low error percentage, around 3.9%. On the other hand, pH, K⁺, and V showed considerably higher errors, around 48%, 72%, and 77%, respectively.

3.4. Comparing Soil Sampling Methods

Table 9 highlights the RMSE obtained through LOOCV cross-validation, demonstrating the accuracy of soil attribute estimates in Fields A and B using different sampling strategies: GRID-1, CONV, CEL, and MZs using Approach 2 and the SAR index VV/VH. The GRID-1 method recorded the lowest RMSE values in both fields, indicating superiority in the precision of soil attribute estimates. In both fields, the CONV method stood out with the highest RMSE values, indicating lower precision in soil attribute estimates. The MZs delineated from SAR data showed superiority compared to the CONV and CEL methods in Field A. In Field B, MZs outperformed the conventional method.

4. Discussion

In the context of remote sensing applied to agricultural fields dedicated to grain production, SAR backscatter time series were used in this study to delineate MZs. Although the Sentinel-1 satellite has a nominal temporal resolution of 12 days for South America, a shorter revisit interval was observed, possibly due to the overlapping of imaging swaths during consecutive satellite passes. Previous studies have demonstrated that higher temporal resolution enhances trend detection, and the identification of spatiotemporal patterns related to crop phenology [44]. [45] pointed out that, in the case of C-band SAR data, temporal resolutions between 3 and 6 days are more suitable for distinguishing crop types and monitoring their phenology, while daily monitoring is necessary to capture rapid changes in soil moisture conditions. Therefore, the higher temporal resolution observed provides improved conditions for understanding and interpreting variations in backscatter over time, potentially contributing to a more accurate delineation of MZs and greater accuracy in estimating soil attributes.

The analysis of the relationship between VV and VH backscatter values and SAR indices revealed that certain indices exhibit strong correlations, indicating potential redundancies. This finding aligns with the study by [46], which identified that the RVI, NRPB, VH/VV, and VV/VH indices show high mutual correlation, with values greater than 0.95 or less than -0.95. Additionally, in both our analysis and the cited study, the RVI4SI index displayed the lowest correlation compared to other indices. Therefore, these findings suggest that, regardless of the agricultural fields studied, the relationship between VV and VH backscatter values and SAR indices tends to follow similar patterns.

From September to December, an increase in VV and VH indices was observed. This increase may be related to the return of the rainy season, which raises soil moisture content. Indeed, during this same period, there is an increase in the monthly accumulated precipitation, elevating soil moisture levels. Additionally, the planting period for agricultural crops, which occurs between September and October, also influences this phenomenon, as increased biomass intensifies signal backscatter [47]. However, during the same period, the VV/VH index showed a decline. According to studies, the VV polarization band is particularly more sensitive to soil moisture compared to the VH band, leading to a reduction in the VV/VH index during this period [48,49].

Between December and January, a stabilization of VV and VH backscatter values is observed. This phenomenon occurs because, with the crop biomass fully developed, there is an attenuation effect from the canopy on the bands, reducing their sensitivity to soil moisture variation. [49] showed that the sensitivity of VV and VH bands to soil moisture variation decreases with the increase in vegetation cover growth (NDVI) and is stronger in the VV polarisation than in the cross-polarization VH. [47] demonstrated that the VV polarization C-band penetrates the maize canopy even when the crop is at its biomass peak (NDVI > 0.7). However, penetration was limited in wheat and pastures. Therefore, during the crop canopy development, vegetation may become the primary component contributing to the volume scattering of the backscattered signal, while the influence of soil may become secondary. Finally, between April and August, there is a strong downward trend in backscatter values for both VV and VH polarizations. This behavior may be associated with the decrease in precipitation during this period, resulting in lower soil moisture content. Since the decrease in backscatter is more pronounced in VV polarization compared to VH, an increase in the VV/VH ratio is observed.

The evaluation of experimental semivariograms in GRID-1, for both fields, highlighted the spatial dependence of soil attributes. The SDI, which relates the nugget effect to the sill to quantify the spatial dependence of these attributes, was found to be less than 75% for most attributes. This indicates strong spatial dependence (less than 25%) and moderate spatial dependence (between 25% and 75%), as suggested by [42]. In this context, kriging emerges as an excellent method for interpolation and estimation of soil attributes in unsampled locations.

The VV and VH backscatter bands, along with the VV/VH and RVI4SI indices, showed variations in the size, shape, and number of MZs when subjected to Approaches 1 and 2. Approach 2, which applies clustering on features extracted from SAR time series via autoencoders, tended to generate more MZs in both fields compared to Approach 1, which performs clustering directly on the time series. Autoencoders belong to a specific class of deep artificial neural networks. They are designed to compress an input into a more compact representation and then reverse that compression, aiming for the reconstructed input to resemble the original as closely as possible [50]. The features extracted by the autoencoder, represented by the compact part, can capture nuances and patterns in the data that raw representation cannot. This leads to a more detailed segmentation of the fields, resulting in a higher number of MZs. Another point to consider is that SAR images are characterized by high levels of noise [51]. Therefore, the use of features extracted by autoencoders represents a less noisy version of the original data, as the learning process of architecture extracts patterns that explain the temporal behavior of the backscatter. This factor may result in more accurate clustering and an increase in the number of MZs.

The VV/VH index, combined with Approach 2 based on autoencoders, tended to exhibit lower RMSE values for soil attribute estimation using the LOOCV strategy. Thus, it was able to produce MZs with greater precision compared to other SAR indices. The integration of VV and VH backscatter band information has shown superior performance compared to the isolated use of each band in various applications [52,53]. This phenomenon is justified by the fact that the VV/VH ratio minimizes acquisition system errors and provides more consistent indications over time than the isolated VH or VV backscatter, as pointed out by [30]. Additionally, certain studies indicate that the VV/VH index correlates more closely with NDVI in specific phenological stages of the crop. This suggests that this index helps in understanding not only the spatial variability of soil moisture but also the canopy structure and crop biomass—crucial aspects for defining MZs [54].

When analyzing the MZs derived from the VV/VH index using autoencoders in both fields, a statistical distinction was observed in at least one mean of each soil attribute originating from the MZs, except for V. Despite the high variability of clay in the region, statistical differences were also detected in temporally unstable attributes, such as Ca²⁺, Mg²⁺, K⁺, etc. The ability of plants to access these attributes is strongly influenced by the soil's ability to retain water in its macro- and micropores. Thus, the sensitivity of SAR data to soil moisture, as evidenced in several studies, is crucial in identifying the variability of these macro- and micronutrients present in the soil [21,49,55].

When analyzing the MZs generated by the VV/VH index using autoencoders in comparison to the randomly created cells (CEL - RND), the potential of the SAR index to delineate MZS was highlighted. Only for the soil attributes K⁺ and V in Field A was observed that in more than 50% of the scenarios, the RMSE of CEL - RND was lower than the RMSE estimated by the MZs. The best performance scenario of the VV/VH index against CEL - RND was observed for clay in both fields. The clay fraction of the soil is intrinsically linked to water retention [56]. Therefore, the sensitivity of SAR data to soil moisture may be one of the explanations for the high correlation observed between the MZs and clay variability in the fields.

When evaluating the various sampling methods, it was observed that the method based on GRID-1 stood out, recording the lowest errors (lower RMSE) for all soil attributes. This result is justified by the fact that the fields investigated in this study exhibited high and moderate spatial dependencies for soil attributes, as indicated by the SDI. In such contexts, the fitting of semivariograms combined with kriging interpolation, the approach adopted in our study, tends to provide good estimates. In contrast, the significant spatial variability suggests that the CONV method, which attempts to represent the field through a single soil sample, may not be efficient. This observation is reinforced by noting that in the fields analyzed in this study, the CONV method had the highest RMSE values, indicating lower accuracy in the estimates of soil attributes.

The MZs delineated from SAR data showed superiority compared to the CONV, CEL, and CEL-RND methods, being occasionally surpassed only by the CEL and CEL-RND methods. Therefore, in scenarios with limited financial resources where conventional sampling is chosen, SAR data can be used to guide sampling through MZs. This methodology presented in this study offers specialists the opportunity to provide services remotely, eliminating the need for field trips. This results in cost savings facilitates the implementation of precision agriculture, even for small farmers. However, future research should be conducted to investigate the impact of reducing the time series in areas without overlapping satellite passes, where temporal resolution consequently decreases. Additionally, as evidenced, there is significant variation in backscatter intensity throughout the year, primarily influenced by fluctuations in precipitation. Thus, it is also suggested that future studies assess the possibility of using images acquired during specific periods of the year.

5. Conclusions

The strategy combining autoencoders with the VV/VH index resulted in more accurate estimates of soil attributes compared to other Synthetic Aperture Radar (SAR) indices. The GRID-1 method, which uses a high-density point grid followed by kriging interpolation, stood out as the most effective technique for mapping soil attributes, while the conventional soil sampling method (CONV) performed the least satisfactorily. The Management Zones (MZs) delineated using the VV/VH index based on autoencoders outperformed the CONV method, the random cell-based soil sampling method (CEL-RND), and, in many cases, the rectangular cell-based soil sampling method (CEL). These findings are encouraging and indicate the potential of SAR data in analyzing soil variability and defining MZs.