Communication

# Using the Negative Soil Adjustment Factor of Soil Adjusted Vegetation Index (SAVI) to Resist Saturation Effects and Estimate Leaf Area Index (LAI) in Dense Vegetation Areas

Zhijun Zhen 1,2 ID, Shengbo Chen 1,*, Tiangang Yin 3 ID, Eric Chavanon 2, Nicolas Lauret 2, Jordan Guilleux 2, Michael Henke 4 ID, Wenhan Qin 1, Lisai Cao 1, Jian Li 1, Peng Lu 1 and Jean-Philippe Gastellu-Etchegorry 2,∗

1 College of Geoexploration Science and Technology, Jilin University, Changchun 130026, China; zhijun.zhen@univ-tlse3.fr (Z.Z.); qinwenhan@jlu.edu.cn (W.Q.); cls@jlu.edu.cn (L.C.); orinatureli@jlu.edu.cn (J.L.); lupeng@jlu.edu.cn (P.L.)

2 CESBIO—UPS, CNES, CNRS, IRD, Université de Toulouse, CEDEX 9, 31401 Toulouse, France; eric.chavanon@cesbio.cnes.fr (E.C.); nicolas.lauret@univ-tlse3.fr (N.L.); guilleuxj@cesbio.cnes.fr (J.G.)

3 Earth System Science Interdisciplinary Center, University of Maryland, College Park, MD 20740-3823, USA; tiangang.yin.85@gmail.com

4 Leibniz Institute of Plant Genetics and Crop Plant Research (IPK Gatersleben), OT Gatersleben, Corrensstr 3, Stadt Seeland, D-06466 Gatersleben, Germany; mhenke@uni-goettingen.de

* Correspondence: chensb@jlu.edu.cn (S.C.); jean-philippe.gastellu@iut-tlse3.fr (J.-P.G.-E.)

<figure><img src="image_1.png" /></figure>

check for
updates

**Citation:** Zhen, Z.; Chen, S.; Yin, T.; Chavanon, E.; Lauret, N.; Guilleux, J.; Henke, M.; Qin, W.; Cao, L.; Li, J.; et al. Using the Negative Soil Adjustment Factor of Soil Adjusted Vegetation Index (SAVI) to Resist Saturation Effects and Estimate Leaf Area Index (LAI) in Dense Vegetation Areas. *Sensors* **2021**, *21*, 2115. https://doi.org/10.3390/s21062115

Academic Editor:
Barmpoutis Panagiotis

Received: 18 February 2021

Accepted: 15 March 2021

Published: 17 March 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_2.png" /></figure>

**Copyright:** © 2021 by the authors.
Licensee MDPI, Basel, Switzerland.
This article is an open access article
distributed under the terms and
conditions of the Creative Commons
Attribution (CC BY) license (https://
creativecommons.org/licenses/by/
4.0/).

**Abstract:** Saturation effects limit the application of vegetation indices (VIs) in dense vegetation areas. The possibility to mitigate them by adopting a negative soil adjustment factor X is addressed. Two leaf area index (LAI) data sets are analyzed using the Google Earth Engine (GEE) for validation. The first one is derived from observations of MODerate resolution Imaging Spectroradiometer (MODIS) from 16 April 2013, to 21 October 2020, in the Apiacás area. Its corresponding VIs are calculated from a combination of Sentinel-2 and Landsat-8 surface reflectance products. The second one is a global LAI dataset with VIs calculated from Landsat-5 surface reflectance products. A linear regression model is applied to both datasets to evaluate four VIs that are commonly used to estimate LAI: normalized difference vegetation index (NDVI), soil adjusted vegetation index (SAVI), transformed SAVI (TSAVI), and enhanced vegetation index (EVI). The optimal soil adjustment factor of SAVI for LAI estimation is determined using an exhaustive search. The Dickey-Fuller test indicates that the time series of LAI data are stable with a confidence level of 99%. The linear regression results stress significant saturation effects in all VIs. Finally, the exhaustive searching results show that a negative soil adjustment factor of SAVI can mitigate the SAVIs' saturation in the Apiacás area (i.e., X = -0.148 for mean LAI = 5.35), and more generally in areas with large LAI values (e.g., X = -0.183 for mean LAI = 6.72). Our study further confirms that the lower boundary of the soil adjustment factor can be negative and that using a negative soil adjustment factor improves the computation of time series of LAI.

**Keywords:** dense forest; google earth engine (GEE); leaf area index (LAI); remote sensing (RS); soil adjustment factor; soil adjusted vegetation index (SAVI)

## 1. Introduction

As major land cover of the planet, forests have become a key priority in studies of the biodiversity and the carbon cycle of terrestrial ecosystems [1,2]. It is essential to record forest dynamics to understand the terrestrial carbon cycle better and improve forest management practices [3,4]. Leaf area index (LAI), defined as the one-sided green leaf area per unit ground area in broadleaf canopies and one-half the total needle surface area per unit ground area in coniferous canopies, is an essential indicator for describing the canopy structure of forest ecosystems, estimating the primary productivity of the stand, and evaluating forest condition over large areas [5,6].
A number of studies link the LAI and radiometric measurements of vegetation. In particular, Vegetation indices (VIs) are often used to estimate LAI from broad spectral bands [7,8]. Although their analytical expressions differ significantly, the implementations of these indices can be divided roughly into three categories: (1) Intrinsic VIs such as simple ratio (SR) [9] and normalized difference vegetation index (NDVI) [10]. (2) Soil adjusted VIs such as soil adjusted vegetation index (SAVI) [11], transformed SAVI (TSAVI) [12], modified SAVI (MSAVI) [13], modified transformed SAVI (MTSAVI) [14], optimized SAVI (OSAVI) [15], and generalized SAVI (GESAVI) [16]. (3) Atmospheric corrected VIs such as atmospherically resistant vegetation index (ARVI) [17], global environment monitoring index (GEMI) [18], enhanced vegetation index (EVI) [19]. Intrinsic VIs [9,10] are widely used because of their simplicity [8]. They are usually the linear combination of reflectance values in the red and near-infrared (NIR) bands. Based on the intrinsic VIs, soil-adjusted VIs [11–16] add a soil adjustment factor and/or soil line parameter (slope and intersection) to mitigate soil noise’s impact on VIs. Atmospheric corrected VIs [17–19] add another atmosphere adjustment factor based on soil adjusted indices to mitigate the impact of atmosphere on VIs.

Despite their wide variety and that they are usually designed to correlate with LAI [12], all VIs inevitably suffer from saturation effects [19]. However, for LAI values more than a certain threshold, the derivative of VIs is relative to LAI decreases. This is the so-called saturation effect [12,20]. Saturation is primarily due to the low sensitivity of reflectance in the red band [8,20]. It directly limits the application of VIs in areas of dense vegetation areas. A number of methods have been developed to address saturation effects on NDVI. Ünsalan and Boyer [21] suggested transforming the NDVI by using an inverse tangent function. However, the transformed NDVI does not boost VIs’ sensitivity for vegetation fractions greater than 0.6 [22]. Gitelson [23] and Vaiopoulos et al. [24] further proposed to adjust the NIR and red reflectances relative contributions to NDVI by adding weighting factors to the NIR reflectance term. However, these weighting factors do not account for the influence of the soil noise and alter the dynamic range of the NDVI, resulting in a range of −0.6 to 0.6 [23].

Therefore, in this paper, we addressed the possibility of adjusting the soil adjustment factor to mitigate the saturation effects. Time series of NDVI, SAVI, TSAVI, and EVI calculated from Sentinel-2, Landsat-5, and Landsat-8 surface reflectance products were used to estimate the time series LAI. The latter ones were validated with two LAI data sets, including a time series of MODIS LAI data over the Apiacás area in Brazil and a global LAI dataset. The positive and negative soil adjustment factor performances were then compared based on the correlation coefficients of the linear regression model.

## 2. Data and Methods

### 2.1. Data

Our study uses two data sets: MODIS 4-Day global 500 m LAI product (MCD15A3H V6 level 4, provided by the National Aeronautics and Space Administration Land Processes Distributed Active Archive Center at the U.S. Geological Survey Earth Resources Observation and Science Center which is located outside Sioux Falls, South Dakota, USA) from time-series observation data and the global LAI product global leaf area index from field measurements (GLAIFM) from field measurements. MODIS LAI data were used to study a region, Apiacá, for which VIs were calculated from Sentinel-2 and Landsat-8 surface reflectance. For the global field-measured LAI dataset, the VIs were calculated with Landsat-5 surface reflectance products.

#### 2.1.1. Local MODIS Time Series LAI Data

Our study includes the entire Apiacás (Latitude: -7.35°--9.82° N, Longitude: 57.04°-58.57° W) in the northernmost of Mato Grosso, Brazil. It was chosen because of its dense vegetation (Figure 1). It contains part of the 19,582 square kilometers of Juruena National Park, one of the largest conservation units in Brazil. Its elevation is around 200 m [25].
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Base map of the study area from the MODIS classification product (MCD12Q1 V6). It has five types of land cover: water body (blue), evergreen broadleaf vegetation (dark green), deciduous broadleaf vegetation (yellow), annual grass vegetation (yellow-green), and urban and built-up lands (dark gray).</figcaption></figure>

The multispectral satellite datasets and MODIS LAI datasets used in this study were provided and processed by the Google Earth Engine [26–28] (GEE, https://earthengine.google.com/, accessed on 11 March 2021). LAI data were from the MODIS 4-Day global 500 m LAI products (MCD15A3H V6 level 4). They were used to validate LAI estimation from VIs calculated from surface multispectral reflectance datasets computed with Sentinel-2 images from 4 February 2019, to 10 November 2020, and Landsat-8 Tier 1 images from 16 April 2013, to 21 October 2020. The use of multiple sources of data helped us to increase the sampling frequency of the time series. Because of the similarity of satellite resolution of Sentinel-2 (10 m) and Landsat-8 (30 m), the VIs differences caused by the two sensors were neglected.

First, the multispectral images were processed. All available Sentinel-2 and Landsat-8 Tier 1 surface reflectance data in GEE before 10 November 2020, were selected and then filtered by two rules: selecting images that cover the study area and an average cloud cover lower than 20%. Then, each remaining image was clipped to ensure that only the study area's pixels were retained. After that, two boundary masks were used in a preprocessing stage. The first was the cloud mask, and the second was the vegetation mask. GEE provides cloud masks of Sentinel-2 and Landsat-8. Pixels in the cloud mask were removed in the subsequent data processing. The vegetation mask of Sentinel-2 was created through the classification of the NDVI and simple ratio NIR/Green Ratio Vegetation Index (GRVI) [29] images using thresholds indicated in a technical report of Sentinel-2 [30]: pixels with (NDVI > 0.40) or (0.36 ≤ NDVI ≤ 0.40 and GRVI > 2.50) were marked as vegetation pixels. The same rule was implemented to create a Landsat-8 vegetation mask. VIs and reflectance values were calculated pixelwise, and the final VIs and reflectance of each satellite image were derived by aggregating VIs and reflectance calculated per pixel. In the next step, LAI was extracted for each corresponding satellite image.

In this work, we used MODIS leaf area index 4-day global 500 m (MCD15A3H V6, https://doi.org/10.5067/MODIS/MCD12Q1.006, accessed on 11 March 2021) [31]. LAI data were searched for each surface reflectance image within eight days of the acquisition time of the surface reflectance image. If no LAI product was available, LAI data within
16-day intervals was searched instead. If no LAI data was available within 16 days, this particular surface reflectance image was ignored. After getting the LAI data in the time range of satellite images, we extracted pixels within both the study area boundary and reflectance image boundary, and a mean LAI value was derived by aggregating all pixels' LAI inside these two boundaries. Finally, for each multispectral image, we had a mean value of different VIs, a mean value of reflectance in various bands, and a mean value of LAI.

Depending on the type and location of sensors, the study area could correspond to several remote sensing images on the same day. In that case, mean values of VIs, reflectance, and LAI were calculated using the images acquired on that day. The reason for considering the entire study area rather than pixels as one study object was to smooth out inconsistencies caused by differences in the resolution of the remotely sensed data.

### 2.1.2. Global Field Measured LAI Data

We used global LAI data from the GLAIFM dataset [32]. This dataset is compiled from around 1000 published estimates of LAI covering the period from 1932 to 2000. These historical LAI data include natural and seminatural (managed) ecosystems and some cultivated vegetation. It contains a wide range of LAI values of 15 biome/land cover classes, from 0.46–2.16 for deserts to 4.4–13.04 for tree plantations. We downloaded it from the website (https://daac.ornl.gov/cgi-bin/dsviewer.pl?ds_id=584, accessed on 11 March 2021).

We extracted the Landsat-5 Tier 1 surface reflectance product for each LAI pixel using the LAI measurement date and geographical information (latitude and longitude). Then, we calculated VIs from these reflectance values. The specific steps are as follows:

(1) Data cleaning: data with null values for latitude, longitude, or date of measured LAI were removed. Also, data with LAI as a range value instead of a specific value were removed; data with LAI measurements spanning over one month were removed.

(2) Screening of data corresponding to dense vegetation: data with LAI values less than 4.0 were removed.

(3) Screening of time range: data with a LAI measurement date range outside the Landsat-5 coverage time range were removed.

(4) Removing duplicate data: some data had the same measurement time and geographical location with different LAI values. On closer examination, we found that some of these data were from measurements of different biomes. Therefore, these data were also removed as we could not establish a one-to-one relationship between LAI and VIs.

(5) Matching LAI with satellite reflectance image: because the temporal accuracy of all measured LAI data is only one month, we searched all surface reflectance data for the month of Landsat-5 Tier 1 data and averaged them to obtain the final surface reflectance data. Some reflectance data could not be retrieved due to cloud cover, so their corresponding LAI data were removed.

Table 1 shows the LAI data and corresponding station information that remained after the above steps. These data were used for evaluating the LAI estimation accuracy of VIs. The LAI value ranges from 4.06 to 10.59, and the sites are mainly located in Canada, Japan, and the USA.
<table><caption>**Table 1.** Field measured LAI data from GLAIFM for evaluating the LAI estimation accuracy of VIs.</caption><thead><tr><th>Site Name</th><th>Latitude (°)</th><th>Longitude (°)</th><th>LAI</th><th>Date</th></tr></thead><tbody><tr><td>BOREAS NSA/OJP, Thompson</td><td>55.92</td><td>-98.62</td><td>4.38</td><td>July 1994</td></tr><tr><td>BOREAS NSA/OBS, Thompson</td><td>55.91</td><td>-98.45</td><td>4.06</td><td>July 1994</td></tr><tr><td>BOREAS NSA, Thompson</td><td>55.91</td><td>-98.52</td><td>8.41</td><td>July 1994</td></tr><tr><td>BOREAS NSA, Thompson</td><td>55.80</td><td>-98.00</td><td>6.21</td><td>July 1994</td></tr><tr><td>BOREAS NSA, Thompson</td><td>55.75</td><td>-97.80</td><td>5.44</td><td>July 1994</td></tr><tr><td>BOREAS SSA, Prince Albert</td><td>54.06</td><td>-105.93</td><td>10.59</td><td>August 1994</td></tr><tr><td>Arakawa River, Urawa</td><td>35.83</td><td>139.62</td><td>4.24</td><td>September 1985</td></tr><tr><td>Westvaco, Summerville, SC</td><td>33.20</td><td>-80.25</td><td>10.4</td><td>February 1991</td></tr></tbody></table>

## 2.2. Methods

### 2.2.1. Vegetation Indices

Four VIs, including NDVI, SAVI, TSAVI, and EVI, were studied. The NDVI is the most commonly used VI. It is defined as:

$$NDVI = \frac{N - R}{N + R} \qquad (1)$$

where *R* and *N* are the atmospherically corrected surface reflectance in the red and NIR bands.

Conversely to NDVI, SAVI and TSAVI have a soil adjustment factor designed to mitigate the soil noise by considering multiple soil scattering. Here, this factor is referred to as “X” with the VIs name as a suffix. SAVI is defined as:

$$SAVI = \frac{(N - R)(1 + X_{SAVI})}{N + R + X_{SAVI}} \qquad (2)$$

where $X_{SAVI}$ is the soil adjustment factor of SAVI. A value equal to 0.5 is recommended in its original paper [11].

TSAVI is defined as:

$$TSAVI = \frac{a \cdot (N - a \cdot R - b)}{a \cdot N + R - a \cdot b + X_{TSAVI} \cdot (1 + a^2)} \qquad (3)$$

where *a* and *b* are the slope and interception of the soil line, respectively. Here these two parameters were set as the constant value of 1.2 and 0.04, respectively, which are considered global soil line parameters [12]. $X_{TSAVI}$ is the soil adjustment factor of TSAVI, recommended to be equal to 0.08 in the original paper [12].

EVI, in addition to a soil adjustment factor, uses atmosphere resistance factors:

$$EVI = G \frac{N - R}{N + C_1 R - C_2 B + X_{EVI}} \qquad (4)$$

where *B* is the atmospherically corrected surface reflectance in the blue band. *G* is a factor gain equal to 2.5. $C_1$, $C_2$ are the aerosol resistance coefficient ($C_1$ = 6, $C_2$ = 7.5), using the blue band reflectance to correct aerosol influences in the red band. $X_{EVI}$ is a soil adjustment factor initially equal to 1 [19]. Since the surface reflectance products were used as multispectral images, we assumed that atmospheric influences had been almost completely removed. Therefore, atmosphere resistance was not considered in this paper.

### 2.2.2. Time Series Analysis

Time series analysis is a statistical method that processes time series data to forecast, control, and understand features of the data. [33]. Time series of LAI can be separated into the trend, seasonality, and residuals. Trend is the increasing (or decreasing) value in the series data. Seasonality is the repeating short-term cycle in the data series. Residuals are the
time series data after the trend, and seasonal components are removed. The decomposition was achieved using classical seasonal decomposition by moving averages [34]. The seasonal multiplication component was selected, and the base-level static analysis was conducted with the rolling mean and Dickey-Fuller test. As the threshold for this analysis, a 99% significance level was employed for the null hypothesis (i.e., no trend). This test explored the trend from a quarter of a year to the next. In the time series analysis, the independent variable was date, and the dependent variable was VIs or LAI. Considering that the field measured LAI data were not continuously observed, only MODIS LAI data were used for time series analysis.

### 2.2.3. Linear Regression

Linear regression is a method to build the linear relationship between a dependent variable and one or more independent variables [35]. We used a unitary linear regression with VI as the independent variable and LAI as the dependent variable. The linear regression equation is defined as:

$$LAI = c \cdot VIs + d \quad (5)$$

where c is the slope and d is the interception of the regression line. The linear regression is widely used to estimate LAI from VIs for its simplicity [8,36]. When LAI increases to a certain threshold value and no longer maintains a linear relationship with VIs, the saturation effect is observed.

### 2.2.4. Exhaustive Search

For discrete problems where no effective solution method is known, it may be necessary to test every possibility sequentially to determine whether it is the solution. Such exhaustive examination of all possibilities is called exhaustive search, direct search, or the "brute force" method [37]. The exhaustive search was used to evaluate the performance of SAVI calculated with exhaustive optimal soil adjustment factor $X_{SAVI}$. Here the independent variable was $X_{SAVI}$, and the dependent variable was $R^2$, slope, and P of the linear regression equation, respectively. The work of Ren et al. [38] demonstrated that a negative $X_{SAVI}$ value (e.g., -0.2) is acceptable, which allowed us to set the exhaustive interval of $X_{SAVI}$ as [-0.3, 1] and to divide the boundary into 2000 parts at intervals of 0.001. The ability to estimate LAI using SAVI was investigated using a $X_{SAVI}$ sequence over the whole study period. The estimation was achieved using a linear regression model between LAI and SAVI calculated with the $X_{SAVI}$ sequence. Both MODIS and field measured LAI data were used here. Finally, the $R^2$ and P were used to evaluate the linear regression performance of various SAVI values calculated with positive or negative $X_{SAVI}$.

## 3. Time Series Analysis

### 3.1. Seasonal Decomposition of LAI

Time series analysis of LAI was conducted to ensure the stability of vegetation conditions in the study area. The result of the decomposition analysis of LAI time series can be seen in Figure 2. All LAI data (labeled as observation in the figure) were averaged using a quarter (i.e., winter, spring, summer, autumn). However, due to cloud cover, few quarters had no corresponding LAI. In that case, a mean value between the preceding and the following values of this null value was assigned. All quarterly LAI data were then decomposed into three parts: trend, seasonal, and residual with a quarterly frequency.

The residual LAI stationarity was validated using the rolling mean (Figure 3) and the Dickey-Fuller test (Table 2). As shown in Figure 3, the means and variances of LAI tended to be constant with little volatility. According to the Dickey-Fuller Test results (Table 2), the statistic value was much smaller than the critical value at 1%, meaning that the data can be considered stable at the 99% confidence level.
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Decomposition analysis for the long-term LAI in the Apiacás (2013–2020). (a) LAI is decomposed into (b) trend, (c) seasonal, and (d) residual. The full-year LAI is separated into four quarters. Then, the mean value of each quarter is calculated and is used for time series analysis.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Stationarity analysis of MODIS LAI residual by rolling mean and standard deviation over the whole observation period.</figcaption></figure>

<table><caption>**Table 2.** Results of Dickey-Fuller Test.</caption><thead><tr><th>Parameters</th><th>Value</th></tr></thead><tbody><tr><td>Test Statistic</td><td>-6.979</td></tr><tr><td>MacKinnon's approximate p-value</td><td>8.281 × 10<sup>-10</sup></td></tr><tr><td>Lags Used</td><td>1</td></tr><tr><td>Number of Observations Used</td><td>26</td></tr><tr><td>Critical value (1%)</td><td>-3.711</td></tr><tr><td>Critical value (5%)</td><td>-2.981</td></tr><tr><td>Critical value (10%)</td><td>-2.630</td></tr></tbody></table>

### 3.2. Time Series Analysis of LAI and VIs

The variations of VIs and LAI time series are plotted in Figure 4. It allows one to compare the performance of VIs on LAI estimation. The VIs were calculated with the default soil adjustment factor ($X_{SAVI} = 0.5$, $X_{TSAVI} = 0.08$, and $X_{EVI} = 1$). LAI values ranged
from 2 to 6.5. NDVI had the highest value over the entire period (i.e., values from 0.7 to 0.9), followed by EVI (i.e., values from 0.4 to 0.7). SAVI and TSAVI shared almost the same range, from 0.2 to 0.6. These high values suggest that the area has high foliage coverage. Besides, Figure 4 shows that VIs can be indicators or not of LAI's trends depending on time. The solid circles highlight a period when LAI has a trend opposite to that of all VIs. The slash circles highlight a period when the trend of LAI is consistent with NDVI while being opposite to other VIs. The dashed circles highlight a period when the trend of LAI is consistent with all VIs.

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Time series of VIs (left axis) and LAI (right axis). The solid circles highlight a period when the LAI has a trend opposite to that of all VIs. The slash circles highlight a period when the trend of LAI is consistent with NDVI while being opposite to other VIs. The dashed circles highlight a period when the trend of LAI is consistent with all VIs.</figcaption></figure>

# 4. Optimal Soil Adjustment Factor X

## 4.1. Linear Regression between VIs and LAI

We tested the ability of VIs for LAI estimation. Figure 5 shows the scatter plot of LAI (y-axis) and VIs (x-axis). The linear regression model was used to compare the VIs. Table 3 gives the coefficients of the linear regression models. Figure 5 and Table 3 clearly show strong saturation effects for all VIs. The VIs hardly showed an upward trend with the increase of LAI. Besides, all VIs showed a weak correlation relationship with LAI. In general, NDVI showed relatively good performance in both cases, with the highest R² (0.1632 in MODIS LAI and 0.4313 in field measured LAI) and the lowest P (0.0173 in MODIS
LAI and 0.2860 in field measured LAI), followed by TSAVI. EVI and SAVI performances were unexpected, especially with the MODIS LAI data, with which they were negatively correlated. In all, the use of VIs to estimate LAI in very densely vegetated areas could lead to large errors.

<figure><img src="image_7.png" /><figcaption>**Figure 5.** Scatter plot between VIs and (a) MODIS and (b) field measured LAI. The solid lines represent the trendlines of the linear regression model.</figcaption></figure>

<table><caption>**Table 3.** Linear regression (slope, interception, R², P) of the four VIs and LAI.</caption><thead><tr><th>VIs</th><th>LAI Type</th><th>Slope</th><th>Interception</th><th>R<sup>2</sup></th><th>P*</th></tr></thead><tbody><tr><td rowspan="2">NDVI</td><td>MODIS LAI</td><td>2.4769</td><td>3.1819</td><td>0.1632</td><td>0.0173</td></tr><tr><td>Field measured LAI</td><td>6.2097</td><td>3.4104</td><td>0.4313</td><td>0.2860</td></tr><tr><td rowspan="2">SAVI</td><td>MODIS LAI</td><td>-2.4165</td><td>6.4499</td><td>-0.0904</td><td>0.1894</td></tr><tr><td>Field measured LAI</td><td>10.1760</td><td>3.9748</td><td>0.2915</td><td>0.4836</td></tr><tr><td rowspan="2">TSAVI</td><td>MODIS LAI</td><td>1.0364</td><td>4.7082</td><td>0.0454</td><td>0.5103</td></tr><tr><td>Field measured LAI</td><td>7.6553</td><td>4.9714</td><td>0.3386</td><td>0.4120</td></tr><tr><td rowspan="2">EVI</td><td>MODIS LAI</td><td>-3.0249</td><td>6.8936</td><td>-0.1504</td><td>0.0285</td></tr><tr><td>Field measured LAI</td><td>4.4827</td><td>5.3874</td><td>0.1017</td><td>0.8106</td></tr></tbody></table>

<footnote>* Two-sided *P* for a hypothesis test whose null hypothesis is that the slope is zero, using Wald Test with t-distribution of the test statistic.</footnote>
## 4.2. LAI Estimation Using SAVI with Negative Soil Adjustment Factor $X_{SAVI}$

We also conducted a linear regression analysis of SAVI and LAI by varying $X_{SAVI}$ from -0.3 to 1 with a step of 0.001 to validate SAVI's robustness using a negative soil adjustment factor. Both MODIS and field measured LAI data were used for the validation. The results of the linear regression are shown in Figure 6.

From both LAI data, we can observe that the best results mainly appeared for $X_{SAVI}$ between -0.18 and 0, indicating that negative $X_{SAVI}$ values were very satisfactory and behaved well compared to positive values. In the MODIS LAI data (Figure 6a), the negative $R^2$ and slope appeared in the interval from 0.225 to 1 for $X_{SAVI}$, indicating that the correlation between SAVI and LAI was extremely weak. Besides, the negative slope meant a negative relationship between LAI and SAVI. The maximin $R^2$ (0.2472) and minimum $P$ (0.0003) appeared at $X = -0.148$. The negative $R^2$ and slope in the field measured data (Figure 6b) appeared for $X_{SAVI}$ between -0.3 and -0.2. Besides, $R^2$ and $P$ were found to be highly volatile in this region. The maximin $R^2$ (0.6417) and minimum $P$ (0.0863) appeared at $X_{SAVI} = -0.183$. With $X_{SAVI}$ from -0.184 to 1, the performance of SAVI gradually worsened, with a decreasing $R^2$ and an increasing $P$.

<figure><img src="image_8.png" /><figcaption>**Figure 6.** Linear regression between SAVI and (a) MODIS and (b) field-measured LAI with varying $X_{SAVI}$. The optimal results are observed in the negative $X_{SAVI}$ region.</figcaption></figure>

# 5. Discussion and Conclusions

## 5.1. Discussion

### 5.1.1. Application Condition of the Negative Soil Adjustment Factor

According to a earlier study [11], the optimal value of the soil adjustment factor is related to the general condition of the vegetation, with $X_{SAVI}$ close but smaller than one for sparse vegetation and close to but not smaller than 0 for dense vegetation. However, the lower boundary zero for dense vegetation was questioned by Ren et al. [38] based on field-measured data in arid grasslands and by Zhen et al. [8,14] based on simulated
data by the discrete anisotropic radiative transfer [39–48] (DART, https://dart.omp.eu/#/, accessed on 11 March 2021) model. Remote sensing models such as DART that simulate the reflectance of land surfaces are very useful tools to investigate the domain of validity of VIs. More generally, they have a great potential to invert satellite observations in terms of LAI maps.

Similar to Ren et al. [38], we found an optimal performance of the negative soil adjustment factor. However, the vegetation conditions differed from the previous study [38] where arid grasslands gave had a NDVI of around 0.2, indicating very sparse vegetation. In our study area, the NDVI values were mainly distributed around 0.8, indicating very dense vegetation. We assumed that the better performance of negative soil adjustment factors in the two situations is because the intersections between the soil line and vegetation isolines were located in the first quadrant, even though their vegetation density varied dramatically.

This paper further confirmed that the lower boundary of the soil adjustment factor could be negative, based on field-measured and remote sensing data (Figure 6). Also, as was already pointed out by Zhen et al. [8] based on simulated data, the determining factors of the location of intersections between the soil line and vegetation isolines are determined by the interception of the vegetation isolines and the soil lines. When the vegetation isoline intercept is larger than the soil line intercept, the intersection point is located in the second or third quadrant of the red-NIR plane. When the vegetation isoline intercept is smaller than the soil line intercept, the intersection point is located in the first quadrant. Moreover, the intercept depends also largely on the spectral properties of the leaves and soil apart from canopy structure parameters such as LAI. Our study further confirmed this conclusion using both field-measured and remote sensing data.

Besides, we noticed that NDVI had a better performance than SAVI, which is considered an improved version of NDVI. Similar undesired SAVI performances have also been reported by [38,49]. Reasons for this poor performance may contribute to the fact that NDVI assumes vegetation isolines converge to the original point, and SAVI assumes vegetation isolines converge to one point (−0.5, −0.5) located in the third quadrant. In our study (very dense vegetation areas), the vegetation isolines converged to one point in the first quadrant. As a result, the assumption based on the original point is much better than the point (−0.5, −0.5). This reason can also be seen in the TSAVI slightly better performance than SAVI because TSAVI assumes that the vegetation isolines' intersection lies between the origin and the point (−0.5, −0.5). Another possibility is that NDVI was used in the backup algorithm for inversion in MODIS LAI data when the radiative transfer inversion method fails, which might also influence the results.

### 5.1.2. Saturation Effect

Numerous studies [8,11–13,15,16] have indicated that VIs show significant saturation effects with LAI growth. Several methods have been implemented to enhance the saturation resistance by adding spectral information or geometry information [19,36,50]. For example, EVI added a blue band which is more sensitive than the red band for high LAI [19]. The absorption intensity of the leaf in the blue band is less than that in the red band. When the LAI exceeded the threshold value, the canopy reflectance had almost no change in the red band due to the strong absorption intensity, and the canopy reflectance still decreased with the LAI because of the weaker absorption intensity in the blue band. Also, modified normalized difference vegetation index (MNDVI) [50] involved a combination of narrow bands in the shorter wavelengths of the red edge (700–750 nm) and longer wavelengths of the red edge (750–780 nm). Besides, NHVI [36] added an angular index representing the distribution of leaves inside the canopy to enhance the saturation resistance ability.

However, the increased information required leads to higher demands on the sensors. For example, EVI products cannot be produced from a sensor without the blue band (i.e., the advanced very high resolution radiometer, AVHRR). The narrowband requirement makes MNDVI hard to be calculated from some common sensors (i.e., Landsat 8). NHVI
requires the sensor can capture multiple-angler observation. Therefore, it is a challenge
for these VIs to produce long-term products across sensor systems with variable spectral
response functions, swath width, and orbiting geometry. Compared with these methods,
adjusting the soil adjustment factor may be effortless because it requires less information
from the sensors. However, although the method using negative soil adjustment factor
achieves better advantages in LAI estimation than other VIs, the accuracy is still lacking
compared to these LAI estimation methods. Our method is, therefore, more suitable for
use in situations where sensor conditions are limited. For example, when only an optic
sensor with a fixed viewing direction is available.

In addition, we found that EVI and SAVI were negatively correlated with LAI in the
MODIS LAI dataset. This negative correlation was not observed in the field measured
data. It leads to inaccurate LAI estimation with errors that cannot be easily estimated.
The negative correlation between VIs and the LAI is few reported in other VIs saturation
studies, and we did not find a solid explanation for it. We simply assume that it is related
to the insensitivity of reflectance values in the red band in the high LAI region. It will need
to be further investigated in the future.

## 5.2. *Conclusions*

The saturation effects limit the use of VIs. In this study, we explored the possibility of
using the negative soil adjustment factor to mitigate the saturation effects. Two data sets
were used for the validation, including a long time series observation of MODIS LAI data
in the Apiacás area and a global LAI dataset. The surface reflectance data derived from
Sentinel-2, Landsat-8, and Landsat-5 were used for VIs calculation. An exhaustive search
was employed to estimate the optimal soil adjustment factor for LAI estimation.

Results show that LAI was highly cyclical with seasonal changes. Besides, VIs calculated using the positive soil adjustment factor had serious saturation effects, and the correlation between VIs and LAI was very low. In the MODIS LAI data, SAVI and EVI even showed a negative growth with LAI. The exhaustive search suggested that a negative $X_{SAVI}$ (e.g., -0.148) in MODIS LAI and $X_{SAVI}$ (e.g., -0.183) in the field measured dataset might be better than the commonly used $X_{SAVI}$ (e.g., 0.5) in a high LAI area.

So far, apart from the arid grasslands and simulated data that have been mentioned,
the feasibility of using negative soil adjustment factors, our study indicates the feasibility of
using a negative soil adjustment factor for saturation resistance in dense vegetation areas.
Our study further confirms that the lower bound of the soil adjustment factor can be less
than zero and reach negative regions. However, although the accuracy of using negative
soil conditioning factors is much higher than positive values, the accuracy of our estimates
still falls short compared to some other methods that are more demanding on the sensor.
Therefore, the method proposed in this paper is more suitable for use with common optical
satellites with single fixed observation direction.

**Author Contributions:** Conceptualization and Writing—original draft, Z.Z.; software, E.C., N.L. and J.G.; review and editing, T.Y., W.Q. and M.H.; supervision, S.C. and J.-P.G.-E.; data, L.C., J.L. and P.L. All authors have read and agreed to the published version of the manuscript.

**Funding:** This paper was supported partly by the three-year action plan for nurturing and developing new industries in the northeastern region of the National Development and Reform Commission ([2016]512), partly by the central government budget, special construction plan for provincial and university, the program for JLU Science and Technology Innovative Research Team (JLUSTIRT, 2017TD-26) which is funded by the Fundamental Research Funds for the Central Universities, China, and partly by China Scholarship Council (CSC, 201906170236), and partly by Graduate Innovation Fund of Jilin University (101832020CX231).

**Data Availability Statement:** Landsat data are available at https://developers.google.com/earth-engine/datasets/catalog/landsat, accessed on 11 March 2021. Sentinel-2 data are available at https://developers.google.com/earth-engine/datasets/catalog/COPERNICUS_S2_SR#description, accessed on 11 March 2021. MODIS leaf area index 4-day global 500 m (MCD15A3H V6) are available at https://doi.org/10.5067/MODIS/MCD12Q1.006, accessed on 11 March 2021. The global leaf area
index from field measurements, 1932–2000 data are available at https://daac.ornl.gov/cgi-bin/dsviewer.pl?ds_id=584, accessed on 11 March 2021.

**Acknowledgments:** The authors would like to thank the editor and reviewers for their valuable help in improving this paper.

**Conflicts of Interest:** The authors declare no conflict of interest.

## References

1. Henke, M.; Huckemann, S.; Kurth, W.; Sloboda, B. Reconstructing leaf growth based on non-destructive digitizing and low-parametric shape evolution for plant modelling over a growth cycle. *Silva. Fenn.* **2014**, *48*, 1019. [CrossRef]

2. Gastellu-Etchegorry, J.; Wang, Y.; Regaieg, O.; Yin, T.; Malenovsky, Z.; Zhen, Z.; Yang, X.; Tao, Z.; Landier, L.; Al Bitar, A. Why to Model Remote Sensing Measurements in 3d? Recent Advances in dart: Atmosphere, Topography, Large Landscape, Chlorophyll Fluorescence and Satellite Image Inversion. In Proceedings of the 2020 5th International Conference on Advanced Technologies for Signal and Image Processing (ATSIP), Sousse, Tunisia, 2–5 September 2020; pp. 1–6.

3. Yin, T.; Gastellu-Etchegorry, J.-P.; Lauret, N.; Grau, E.; Rubio, J. A new approach of direction discretization and oversampling for 3d anisotropic radiative transfer modeling. *Remote Sens. Environ.* **2013**, *135*, 213–223. [CrossRef]

4. Li, J.; Chen, S.; Qin, W.; Murefu, M.; Wang, Y.; Yu, Y.; Zhen, Z. Analysis of accuracy of modis brdf product (mcd43 c6) based on misr land surface brf product-a case study of the central part of northeast asia. *Int. Arch. Photogramm. Remote Sens. Spat. Inf. Sci.* **2018**, *42*, 819–823. [CrossRef]

5. Yan, G.; Hu, R.; Wang, Y.; Ren, H.; Song, W.; Qi, J.; Ling, C. Scale effect in indirect measurement of leaf area index. *IEEE Trans. Geosci. Remote Sens.* **2016**, *54*, 3475–3484. [CrossRef]

6. Gastellu-Etchegorry, J.P.; Yin, T.; Lauret, N.; Cajgfinger, T.; Gregoire, T.; Grau, E.; Feret, J.B.; Lopes, M.; Guilleux, J.; Dedieu, G. Discrete anisotropic radiative transfer ( dart 5) for modeling airborne and satellite spectroradiometer and lidar acquisitions of natural and urban landscapes. *Remote Sens.* **2015**, *7*, 1667–1701. [CrossRef]

7. Yan, G.; Ronghai, H.U.; Luo, J.; Xihan, M.U.; Xie, D.; Zhang, W. Review of indirect methods for leaf area index measurement (in chinese with english abstract). *J. Remote Sens.* **2016**, *20*, 958–978.

8. Zhen, Z.; Chen, S.; Qin, W.; Yan, G.; Gastellu-Etchegorry, J.; Cao, L.; Murefu, M.; Li, J.; Han, B. Potentials and limits of vegetation indices with brdf signatures for soil-noise resistance and estimation of leaf area index. *IEEE Trans. Geosci. Remote Sens.* **2020**, *58*, 5092–5108. [CrossRef]

9. Jordan, C.F. Derivation of leaf-area index from quality of light on the forest floor. *Ecology* **1969**, *50*, 663–666. [CrossRef]

10. Tucker, C.J.; Vanpraet, C.L.; Sharman, M.J.; Van Ittersum, G. Satellite remote sensing of total herbaceous biomass production in the senegalese sahel: 1980–1984. *Remote Sens. Environ.* **1985**, *17*, 233–249. [CrossRef]

11. Huete, A.R. A soil-adjusted vegetation index (savi). *Remote Sens. Environ.* **1988**, *25*, 295–309. [CrossRef]

12. Baret, F.; Guyot, G. Potentials and limits of vegetation indices for lai and apar assessment. *Remote Sens. Environ.* **1991**, *35*, 161–173. [CrossRef]

13. Qi, J.; Chehbouni, A.; Huete, A.R.; Kerr, Y.H.; Sorooshian, S. A modified soil adjusted vegetation index. *Remote Sens. Environ.* **1994**, *48*, 119–126. [CrossRef]

14. Zhen, Z.; Chen, S.; Qin, W.; Li, J.; Mike, M.; Yang, B. A modified transformed soil adjusted vegetation index for cropland in jilin province, china. *Acta Geol. Sin. Engl. Ed.* **2019**, *93*, 173–176. [CrossRef]

15. Rondeaux, G.; Steven, M.; Baret, F. Optimization of soil-adjusted vegetation indices. *Remote Sens. Environ.* **1996**, *55*, 95–107. [CrossRef]

16. Gilabert, M.A.; González-Piqueras, J.; García-Haro, F.J.; Meliá, J. A generalized soil-adjusted vegetation index. *Remote Sens. Environ.* **2002**, *82*, 303–310. [CrossRef]

17. Kaufman, Y.J.; Tanre, D. Atmospherically resistant vegetation index (arvi) for eos-modis. *IEEE Trans. Geosci. Remote Sens.* **1992**, *30*, 261–270. [CrossRef]

18. Pinty, B.; Verstraete, M.M. Gemi: A non-linear index to monitor global vegetation from satellites. *Vegetatio* **1992**, *101*, 15–20. [CrossRef]

19. Huete, A.R.; Liu, H.Q.; Batchily, K.; Leeuwen, W.V. A comparison of vegetation indices over a global set of tm images for eos-modis. *Remote Sens. Environ.* **1997**, *59*, 440–451. [CrossRef]

20. Huete, A.R.; Liu, H.; van Leeuwen, W.J. The Use of Vegetation Indices in Forested Regions: Issues of Linearity and Saturation. In Proceedings of the IGARSS’97. 1997 IEEE International Geoscience and Remote Sensing Symposium Proceedings. Remote Sensing-A Scientific Vision for Sustainable Development, Singapore, 3–8 August 1997; pp. 1966–1968.

21. Unsalan, C.; Boyer, K.L. Linearized vegetation indices based on a formal statistical framework. *IEEE Trans. Geosci. Remote Sens.* **2004**, *42*, 1575–1585. [CrossRef]

22. Jiang, Z.; Huete, A.R.; Li, J.; Chen, Y. An analysis of angle-based with ratio-based vegetation indices. *IEEE Trans. Geosci. Remote Sens.* **2006**, *44*, 2506–2513. [CrossRef]

23. Gitelson, A.A. Wide dynamic range vegetation index for remote quantification of biophysical characteristics of vegetation. *J. Plant Physiol.* **2004**, *161*, 165–173. [CrossRef]
24. Vaiopoulos, D.; Skianis, G.; Nikolakopoulos, K. The contribution of probability theory in assessing the efficiency of two frequently used vegetation indices. *Int. J. Remote Sens.* **2004**, *25*, 4219–4236. [CrossRef]

25. Unidade de Conservação: Parque Nacional do Juruena. (In Portuguese). Available online: http://sistemas.mma.gov.br/cnuc/index.php?ido=relatorioparametrizado.exibeRelatorio&relatorioPadrao=true&idUc=281 (accessed on 11 March 2021).

26. Mutanga, O.; Kumar, L. Google earth engine applications. *Remote Sens.* **2019**, *11*, 591. [CrossRef]

27. Kumar, L.; Mutanga, O. Google earth engine applications since inception: Usage, trends, and potential. *Remote Sens.* **2018**, *10*, 1509. [CrossRef]

28. Gorelick, N.; Hancher, M.; Dixon, M.; Ilyushchenko, S.; Thau, D.; Moore, R. Google earth engine: Planetary-scale geospatial analysis for everyone. *Remote Sens. Environ.* **2017**, *202*, 18–27. [CrossRef]

29. Henrich, V.; Krauss, G.; Götze, C.; Sandow, C. The Indexdatabase. Available online: https://www.indexdatabase.de/ (accessed on 12 March 2021).

30. Richter, R.; Louis, J.; Müller-Wilm, U. Sentinel-2 msi—level 2a products algorithm theoretical basis document. *ESASP* **2012**, *49*, 1–72.

31. Li, J.; Chen, S.; Qin, W.; Useya, J.; Zhen, Z.; Wang, Y. A fast reprojection method for modis products with sinusoidal projection. *J. Indian Soc. Remote Sens.* **2018**, *46*, 1563–1567. [CrossRef]

32. Scurlock, J.; Asner, G.; Gower, S. Global Leaf Area Index from Field Measurements, 1932–2000; Oak Ridge National Laboratory Distributed Active Archive Center: Oak Ridge, TN, USA, 2001.

33. Hamilton, J.D. *Time Series Analysis*; Princeton University Press: Princeton, NJ, USA, 1994.

34. Hirsch, R.M.; Slack, J.R. A nonparametric trend test for seasonal data with serial dependence. *Water Resour. Res.* **1984**, *20*, 727–732. [CrossRef]

35. Montgomery, D.C.; Peck, E.A.; Vining, G.G. *Introduction to Linear Regression Analysis*; Wiley-Interscience: Hoboken, NJ, USA, 2016.

36. Hasegawa, K.; Matsuyama, H.; Tsuzuki, H.; Sweda, T. Improving the estimation of leaf area index by using remotely sensed ndvi with brdf signatures. *Remote Sens. Environ.* **2010**, *114*, 514–519. [CrossRef]

37. Weisstein, E.W. Exhaustive Search. Available online: https://mathworld.wolfram.com/ExhaustiveSearch.html (accessed on 11 March 2021).

38. Ren, H.; Zhou, G.; Zhang, F. Using negative soil adjustment factor in soil-adjusted vegetation index (savi) for aboveground living biomass estimation in arid grasslands. *Remote Sens. Environ.* **2018**, *209*, 439–445. [CrossRef]

39. Gastellu-Etchegorry, J.-P.; Lauret, N.; Yin, T.; Landier, L.; Kallel, A.; Malenovský, Z.; Al Bitar, A.; Aval, J.; Benhmida, S.; Qi, J. Dart: Recent advances in remote sensing data modeling with atmosphere, polarization, and chlorophyll fluorescence. *IEEE J. Sel. Top. Appl. Earth Observ.* **2017**, *10*, 2640–2649. [CrossRef]

40. Gastellu-Etchegorry, J.-P.; Yin, T.; Lauret, N.; Grau, E.; Rubio, J.; Cook, B.D.; Morton, D.C.; Sun, G. Simulation of satellite, airborne and terrestrial lidar with dart (i): Waveform simulation with quasi-monte carlo ray tracing. *Remote Sens. Environ.* **2016**, *184*, 418–435. [CrossRef]

41. Gastellu-Etchegorry, J.P.; Lauret, N.; Yin, T.; Landier, L.; Al, B.A.; Guilleux, J.; Jan, C.; Chavanon, E. Dart: Radiative transfer modeling for simulating terrain, airborne and satellite spectroradiometer and lidar acquisitions and 3d radiative budget of natural and urban landscapes. In Proceedings of the 2016 IEEE International Geoscience and Remote Sensing Symposium (IGARSS), Beijing, China, 10–15 July 2016; pp. 3632–3635.

42. Gastellu-Etchegorry, J.P.; Demarez, V.; Pinel, V.; Zagolski, F. Modeling radiative transfer in heterogeneous 3-d vegetation canopies. *Remote Sens. Environ.* **1996**, *58*, 131–156. [CrossRef]

43. Gastellu-Etchegorry, J.P.; Bruniquel-Pinel, V. A modeling approach to assess the robustness of spectrometric predictive equations for canopy chemistry. *Remote Sens. Environ.* **2001**, *76*, 1–15. [CrossRef]

44. Yin, T.; Lauret, N.; Gastellu-Etchegorry, J.P. Simulation of satellite, airborne and terrestrial lidar with dart (ii): Als and tls multi-pulse acquisitions, photon counting, and solar noise. *Remote Sens. Environ.* **2016**, *184*, 454–468. [CrossRef]

45. Yin, T.; Feret, J.B.; Gastellu-Etchegorry, J.P.; Lauret, N. Data simulation and fusion of imaging spectrometer and lidar multi-sensor system through dart model. In Proceedings of the 2016 IEEE International Geoscience and Remote Sensing Symposium (IGARSS), Beijing, China, 10–15 July 2016; pp. 3640–3643.

46. Yin, T.; Kotthaus, S.; Gastellu-Etchegorry, J.-P.; Morrison, W.; Norford, L.; Grimmond, C.; Lauret, N.; Chrysoulakis, N.; Al Bitar, A.; Landier, L. Atmospheric correction of ground-based thermal infrared camera through dart model. In Proceedings of the IGARSS 2017—2017 IEEE International Geoscience and Remote Sensing Symposium, Fort Worth, TX, USA, 23–28 July 2017; pp. 5885–5888.

47. Gascon, F.; Gastellu-Etchegorry, J.P.; Leroy, M. Using multi-directional high-resolution imagery from polder sensor to retrieve leaf area index. *Int. J. Remote Sens.* **2007**, *28*, 167–181. [CrossRef]

48. Gascon, F.; Gastellu-Etchegorry, J.P.; Lefevre-Fonollosa, M.J.; Dufrene, E. Retrieval of forest biophysical variables by inverting a 3-d radiative transfer model and using high and very high resolution imagery. *Int. J. Remote Sens.* **2004**, *25*, 5601–5616. [CrossRef]

49. Ren, H.; Zhou, G.; Zhang, X. Estimation of green aboveground biomass of desert steppe in inner mongolia based on red-edge reflectance curve area method. *Biosys. Eng.* **2011**, *109*, 385–395. [CrossRef]

50. Mutanga, O.; Skidmore, A.K. Narrow band vegetation indices overcome the saturation problem in biomass estimation. *Int. J. Remote Sens.* **2004**, *25*, 3999–4014. [CrossRef]