Spatial variability of biophysical multispectral indexes under heterogeneity and anisotropy for precision monitoring

HIGHLIGHTS: The multispectral indexes successfully represented the spatial distribution of land use in a highly heterogeneous field. Better results were obtained when anisotropy was considered. Soil-Adjusted Vegetation Index represented better the spatial structure of the anisotropic field. ABSTRACT The study aimed to characterize the spatial structure of variability of biophysical indexes of vegetation through images obtained by Unmanned Aerial Vehicles under strong heterogeneity and anisotropy, using geostatistical procedures. Plots with different types and densities of culture were evaluated in a didactic vegetable garden. Five vegetation indexes obtained from aerial multispectral camera images were evaluated parallel with geostatistical analysis and anisotropy investigation for multiscale spatial modeling. For the studied domain, geometric anisotropy was identified for the biometric indexes. The spherical model presented a better fit when anisotropy was not considered, whereas the exponential model had the best performance in the anisotropic analysis. Contrasting targets were better identified in multispectral images and considering anisotropy. The Soil-Adjusted Vegetation Index is recommended for similar applications.


Introduction
Spatial variability analysis of pedological, climatic, and biological attributes is essential for the development of precision agriculture (Pradipta et al., 2022).Moreover, the application of remote sensing techniques allows high-resolution data to be obtained, thus, enhancing the understanding of spatial variability and heterogeneities of agricultural fields.
Spectral indexes based on remote sensing applying UAVs (Unmanned Aerial Vehicles) have been widely applied to the assessment of vegetation and soils (Sishodia et al., 2020), crop classification (Chen et al., 2020), long-term changes in vegetation growth (Lu et al., 2020), surface and subsurface hydrology (Vélez-Nicolás et al., 2021), predicting yield (Andrade Júnior et al., 2022), and crop diseases (Belmonte et al., 2023).Data from UAV images usually encompasses multiple spatial scales of variability.Furthermore, pixels of a high-resolution UAV image cannot be considered independent of each other.Geostatistics is largely used to account for such spatial correlations and allows the understanding of structural and vector characteristics.The incorporation of anisotropy is necessary in several studies of surface reflectance imaging in cases where the vegetative indexes are conditioned to the adopted cartesian direction, interfering with the optical and structural properties of plants and soil (Lin et al., 2021).Despite many studies applying vegetation indexes from UAV images for mapping, investigations addressing small-scale autocorrelation structures in anisotropic fields are still incipient, particularly in northeast of Brazil.
Given the above, the present study aimed to analyze the spatial structure of variability and dependence of biophysical vegetation indexes obtained using UAV, investigating the small-scale variability in the presence of strong anisotropy and the efficiency of different indexes in identifying zones of high heterogeneity.

Material and Methods
Images of the flight were obtained from a dedicated Unmanned Aircraft System (UAS) flight conducted on September 20, 2020, over an experimental vegetable garden at the Department of Agronomy of the Universidade Federal Rural de Pernambuco (UFRPE) (8° 1' 3.06" S and 34° 56' 44.69" W, and altitude of 38 m) in Recife, Pernambuco state, Brazil.The vegetable garden is adopted for university teaching, research, and extension services.In the area, there are rectangular plots (11 × 1 m²) in which plants of different vegetable crops are grown and on which imaging was carried out.
A DJI Phantom 4 Multispectral equipment was used to perform the flight at the height of 30 m, over an area of approximately 2,582 m² (Figure 1A), adopting a spatial image resolution equal to 0.5 m.The drone is equipped with an RGB camera and a multispectral camera array with five cameras covering the Blue, Green, Red, Red Edge, and Near Infrared ranges -all at 2 MP with a global shutter on a 3-axis stabilized gimbal.
A 200 m² polygon was selected in the image for the spectral analysis.In particular, the polygon presents plots with different species and plant architectures, which allows working with the analysis of different vegetation cover and behavior of the targets.In the area of the selected polygon, a sampling grid was then established with 0.5 spacing between points, preserving the spatial image resolution, which totaled 794 observations, in which the index values were extracted.They can be observed in the grid of the sample mesh image (Figure 1B).
The extraction of images was performed with the DroneDeploy application software, and the processing was performed in the cloud through that platform, generating an orthomosaic by merging the aerial images.The QGIS software (QGIS Geographic Information System Open-Source Geospatial Foundation Project) was adopted to determine the biophysical indexes and extract the values of the variables at each sampled point.
The irrigation system consists of micro-sprinklers, and there was no irrigation in the month before the experiment, with a recorded rainfall of 69.96 mm, average temperature of 25.63 °C, maximum temperature of 29.23 °C, minimum temperature of 22.02 °C, relative air humidity of 78.53%, and solar radiation of 19.75 MJ m -2 per day (Figure 1C).The substrate used at the experimental site is a clayey organic garden substrate composed of 42% clay, 28% silt, and 18% sand, with an electrical conductivity of 0.3 dS m -1 and no salinity issues.
Predominantly the data acquired by drones are in the visible spectral range (blue, green, and red -400 to 700 nm) (Fernandez-Gallego et al., 2019) and near-infrared (Marty et al., 2022).As a result, spectral indexes were selected to evaluate the vegetation cover for these two spectral bands.
The spectral indexes evaluated that consider the nearinfrared were the Soil-Adjusted Vegetation Index (SAVI), developed by Huete (1988), which consists of an improvement of the Normalized Difference Vegetation Index (NDVI) according to the application of the L constant to minimize the effects of soil color on the results.The L constant ranges between 0 and 1, according to the vegetation density, with the value 1 being adopted for areas with little vegetation cover, whereas in areas with intermediate cover, it is used L = 0.5; when L is equal to 0, SAVI is identical to NDVI (Huete, 1988).For such a case, the L value adopted was 0.5.The index was calculated for the area of interest using Eq. 1.

( )( )
1 L nir red SAVI nir red L where: nir -refers to the near-infrared spectral band; red -red band; and L -the parameter associated with the vegetation cover.
The NDVI was also calculated (Rouse et al., 1974), considering that it is the most usual index with significant robustness in the literature.Values were obtained by applying Eq. 2.

( ) ( )
nir red NDVI nir red Regarding the visible spectrum indexes applied in the study of the area, the Normalized Red Green Difference Index (NGRDI) (Tucker, 1979) and the Modified Green Red Vegetation Index (MGRVI) (Bendig et al., 2015) were calculated.The values were obtained by Eqs. 3 and 4, respectively.
where: Fa is the anisotropy factor and consists of the direction of greater continuity, obtained by dividing the range (A) of the models in the direction of greater and lesser spatial continuity, and θ is the angle of the direction of the greatest continuity.
The theoretical model was then adjusted to the experimental omnidirectional semivariogram and the corrected directional semivariograms for each index, thus determining the nugget effect (C 0 ), threshold (Cs), and range (a) values.The goodness of fit of the models was evaluated using the Sum of the Q-squares of the Errors or Residuals (SSErr) determined from Eq. 9.

( ) ( )
green red NGRDI green red green red MGRVI green red where: green -green band; and red -red band.
The data were submitted to descriptive statistical analysis to determine the mean, variance, standard deviation, maximum and minimum value, coefficient of variation, symmetry, and kurtosis.
The methodology presented by Li ( 2017) was followed for the spatial characterization of the data and subsequent geostatistical analysis of the spectral indexes.Initially, an experimental omnidirectional semivariogram was determined with adjustment for the distance -Lag (h), the tolerance angle (∆θ), and the maximum distance (Max dist.) through Eq. 5.
where: γ(h) is the estimated value of the semivariance of the experimental data; Z(X i +h) and Z( i ) are observed values of the regionalized variable; and N(h) denotes the number of pairs of measured values separated by distance h.Then, the experimental semivariogram for each index was plotted as a preliminary study of the anisotropy of the experimental semivariogram values of the studied variable.Afterward, the semivariogram values were calculated firstly in the four main directions with azimuth angle (θ) equal to 0°, 45°, 90°, and 135° separately to identify the existence and type of anisotropy (geometric, zonal, or both) and investigate the major and minor anisotropy in each direction.Afterward, analysis was carried out for other complementary directions.
After fitting the models for each direction on the same graph, the directional models were rotated into a single consistent model for all directions, according to the methodology described by Almeida & Bettini (1994).The correction for the case of geometric anisotropy is performed by linear transformations, in which procedures of dilation and rotation of spatial coordinates are used through Eq. 6.

( ) r d h ' x, y M M =
where: h' is the transformation of h, and the rotation (Mr) and expansion (Md) matrices are given by: ∑ where: γ(h i ) is the predicted semivariogram value using the fitted theoretical model, and γ (h i ) is the mean value of the experimental semivariogram at a h i -distance lag.
The spatial dependence evaluated according to C 0 is considered strong when the nugget effect is lower than or equal to 25% of the sill; moderate when the nugget effect is between 26-75% of the sill; and weak when it is above 75% (Cambardella et al., 1994).
Cross-validation was also performed with the adjusted semivariogram model, and the Mean Reduced Error (Rε), and the Reduced Variance (S 2 Rε ) were calculated as evaluation metrics based on Eqs. 10 and 11 (Vauclin et al., 1983).Crossvalidation is achieved by removing the first Z data value at the S i location (i = 1 to N) and using the remaining values (N-1) of data sampled over the area of study to fit the theoretical and predicted semivariogram model Z*(S i ), and σ k Kriging Standard Deviation.The cross-validation process is repeated for all data values sampled at all N S i sites.

( ) (
) Finally, ordinary kriging was performed based on the adjusted models to determine the values at unsampled locations among the points and contour plotting to assess the results of each index.The procedures were performed in the R Core Team software.

Results and Discussion
The results of the descriptive statistical analysis of the vegetation spectral indexes evaluated in the study are presented in Table 1.The mean values of SAVI, NDVI, MGRVI, and NGRDI were 0.34, 0.23, 0.13, and 0.07, respectively.
From the descriptive data characterization, it is verified that the distribution of indexes does not follow a normal function based on the obtained asymmetry coefficient, with absolute values higher than 0.5.According to the coefficient of variation (CV) classification proposed by Warrick & Nielsen (1980), the data are characterized as presenting high variation (CV > 60%).Furthermore, the evaluated kurtosis coefficient, the data distribution is Platykurtic with short tails (Ghasami et al., 2020), consistent with the data graphic characterization behavior (Figure 2).
The infrared spectrum indexes are close to those obtained by Zanzarini et al. (2013), characterizing spectral indexes from orbital sensors for assessing land use and occupation classes.However, due to scale effects, the standard deviation values differed from those obtained by the authors.
It can be observed that the amplitude of the values of infrared indexes obtained differs from those of the visible Table 1.Descriptive statistical analysis of the spectral indexes SAVI -Soil Adjusted Vegetation Index; NDVI -Normalized Difference Vegetation Index; MGRVI -Modified Green Red Vegetation Index; NGRDI -Normalized Red Green Difference Index; CV -Coefficient of Variation; SD -Standard Deviation Although they all present values ranging from -1 to 1 depending on the different spectral targets, a distribution of such values is observed along the range of variation, mainly of the infrared indexes, being possible to verify a persistence of values between -0.4 and 0.4 for the NGRDI.This index consists of an adaptation of the NDVI for the application with sensors with only the visible range, replacing the infrared with green.
However, the NGRDI index is not linearly correlated with the NDVI, and the index presents a satisfactory performance when evaluating areas with exposed soil or low vegetation density, corresponding to negative values close to zero (Figure 2).Such values represent these targets, similar to the others.However, in areas with higher vegetation density (NDVI > 0.7), the use of the index is hampered by the saturation of its capacity (Hunt et al., 2005).
Despite the narrower amplitude observed for the visible spectral indexes, all values were able to highlight areas with vegetation cover, contrasting with those with exposed soil.However, as observed by Agapiou (2020), the performance of visible indexes is influenced by the area size and heterogeneity.Indeed, the author observed better performances for homogeneous targets, which differs from the condition of the imaged area in this study.
For the analysis of spatial variability, the frequency distribution pattern of the spectral index was investigated.Although high variability is verified, a stationary variance distribution in the area is observed for all indexes (Figure 2).The experimental omnidirectional semivariogram was then evaluated, initially disregarding the influence of anisotropy on the data.The fitted theoretical models and the goodness of fit can be seen in Figure 3.The model that presented the best fit for all indexes, with the lowest SSErr value, was the spherical one, with errors between 0.0008 and 0.097, as also obtained by Teixeira et al. (2018), applying geostatistics for soil mapping.However, the model adjustment error obtained in this study was much lower than that observed by these authors.The spherical model, adjusted to the experimental variograms, has the characteristic of adapting to attributes with abrupt changes in space (Isaaks & Srivastava, 1989), as is the case of the data studied, due to the presence of beds, which produce sharp spectral targets.
According to Cambardella et al. (1994), the spherical model best describes the spatial distribution of soil attributes showing high variability.On the other hand, the model that presented the highest adjustment error was the Gaussian, with 0.471 for the SAVI index.
None of the models presented the nugget effect (C 0 ) equal to zero, accounting for the sum of variabilities due to measurement errors and those occurring on a scale smaller than the adopted one (Isaaks & Srivastava, 1989).
In this case, the exponential model showed a strong to moderate spatial dependence for all indexes, while the spherical model was on the threshold for moderate dependence for the SAVI index and showed low dependence for the visible spectrum indexes.For the NDVI and NGRDI, the spherical model presented a strong spatial dependence.As for the range parameter (A), the highest values were obtained for the spherical model, which ranged from 1.95 to 2.02 m.This parameter refers to the maximum distance at which spatial dependence occurs (Isaaks & Srivastava, 1989).
Evaluating the spectral indexes of satellite images, Zanzarini et al. (2013) obtained moderate spatial dependence for the spherical model.According to Montenegro & Montenegro (2006), better estimates are obtained when the models are based on experimental semivariograms that present the lowest "nugget/sill effect" ratio and the greatest range.The spherical model presents the lowest SSErr and the highest Range (A) for all indexes.However, it is not the model that presents the highest value of spatial dependence and the lowest ratio (C 0 / (C 0 + C 1 ) for all cases, which indicates that not necessarily only one model should be considered as the best fit for the experimental semivariograms of spectral indexes (Li, 2017).
In this study, due to the greater combination of favorable fit evaluation attributes, it was decided to proceed with the kriging simulations using the spherical model for all indexes.It should be noted that none of the spherical omnidirectional models presented range values consistent with the scales of variability inherent to the crop beds, as expected.Such a result is related to the mean directional behavior of those semivariogram models, which smooth variability along different directions in anisotropic fields, as Friedland et al. (2017) pointed out.
In order to investigate the structural spatial distribution behavior of the indexes across directions, as a way of accounting for the anisotropy, the experimental directional semivariograms were calculated for several directions, plotting them in a single graph aiming to visually identify the degree of directional variation (Figure 4).
Anisotropic behavior of the geometric type is observed for the evaluated data of the spectral indexes.In geometric anisotropy (Figure 4), the range varies according to the directions but under the same sill (Almeida & Bettini, 1994).Based on the generated directional semivariogram map, directions of minimum and maximum spatial continuity of the indexes under study were identified, with the maximum between 103 and 104º and the minimum at 15º.The anisotropic model was then corrected according to the methodology described in Li (2017) through linear transformations, in which procedures of dilation and rotation of spatial coordinates were used based on angular properties of tangent lines for the anisotropic semivariograms (Barbosa et al., 2019).
After the anisotropy correction, the theoretical model that presented the best fit for the indexes was then the exponential, showing SSErr values of 0.04366 (SAVI), 0.05443 (MGRVI), The range parameters varied from 1.4 to 5.4 m, consistent with the width and length of the cultivation beds.Moreover, it was verified that the maximum spatial correlation scale along the beds is lower than the bed length.All the obtained directional semivariograms reached a sill.Observations cannot be assumed spatially independent for separation distances smaller than the range value.Classical statistical methods must be modified accordingly for such scales, as Belmonte et al. (2023) pointed out.
The corrected anisotropic model for each index was then applied to generate their spatial distribution using ordinary kriging.The generated anisotropic maps can be seen in Figure 5, compared to the omnidirectional maps.
It can be observed that the SAVI, NDVI, MGRVI, and NGRDI indexes show a remarkable difference in the spatial distribution, with a predominance of the distribution of values in the main direction (103.97º) in the anisotropic model, corresponding to the direction of the beds.This behavior would be ignored if only the isotropic models were considered, which does not represent this predominance of distribution in the direction of the angle of highest variation of the variable, even under an adequate adjustment to the experimental data.
For the visible indexes, it is possible to observe an anisotropic distribution of values in the greater variance direction, even at a lower intensity.This behavior may be associated with a lower ability to characterize the spectral response of visible indexes (Breunig et al., 2019).It can be verified as a spectral mixture in the pixels hampering vegetation delimitation in the beds.
The spectral response of the spontaneous vegetation that borders the beds may have been mixed with that of cultivated vegetable crops.Hence, the variance and anisotropic behavior of the indexes were partially altered.Another possibility is the intensification of anisotropy due to the distortion caused by the reflectance of the exposed soil, better characterized by the infrared indexes.According to Lin et al. (2021), only vegetation coverings higher than 50% of the area would suffer distortions close to zero from the influence of the soil, reducing the model uncertainties.
The cross-validation parameters for the estimates produced by the isotropic and anisotropic models can be seen in Table 2.For all indexes, both omnidirectional and anisotropic models presented a reduced error close to zero and the Reduced Variance close to 1, which are required for proper experimental data modeling.
For the mean reduced error (Rε) of the anisotropic models, clear increments in the quality obtained with the correction are observed.The anisotropic model adjusted for the SAVI index presented an error of -0.0003495, lower in absolute value than that associated with the omnidirectional model (-0.0005322).It means that particular aspects of the dataset are better represented when adopting the anisotropic modeling, hampering smoothing the spatial distribution of spectral index data.
Improvements in the model performance are more evident for the SAVI index, which better represented the

Conclusions
1.The biophysical indexes obtained using Unmanned Aerial Vehicles present high spatial autocorrelation in the study area across different directions.The data-adjusted directional semivariogram indicated the presence of geometric anisotropy.
2. The adopted experimental resolution for the Unmanned Aerial Vehicles flight allowed a sill value to be achieved for all directional semivariogram models for the biophysical indexes.
3. Multispectral images enabled the identification of zones of high heterogeneity, allowing a more precise classification of areas with exposed soil.
4. The Soil-Adjusted Vegetation Index provided the best representation of the heterogeneous targets across directions of the anisotropic field.
5. The obtained results provide guidance for larger-scale fields, particularly in areas with exposed soils (bare soil).
Figure 1.Study area in the experimental vegetable garden of UFRPE (A), sampling grid and distribution of points (B), and rainfall and temperature (C)

Figure 2 .
Figure 2. Analysis of the spatial distribution of data variance.Histogram of spectral indexes values and the values obtained for the study area X Coord and Y Coord -X and Y coordinates, respectively; SAVI -Soil Adjusted Vegetation Index; NDVI -Normalized Difference Vegetation Index; MGRVI -Modified Green Red Vegetation Index; NGRDI -Normalized Red Green Difference Index

Figure 3 .
Figure 3. Theoretical semivariograms (spherical, exponential, and Gaussian models) fitted to the data for the SAVI, MGRVI, NDVI, and NGRDI spectral indexes, and their adjustment parameters SAVI -Soil Adjusted Vegetation Index; NDVI -Normalized Difference Vegetation Index; MGRVI -Modified Green Red Vegetation Index; NGRDI -Normalized Red Green Difference Index; C 0 -Nugget effect; C 0 + C -Sill; A -Range and model adjustment evaluated by the Sum of Squares of Error or Residuals (SSErr)