AM–GM Algorithm for Evaluating, Analyzing, and Correcting the Spatial Scaling Bias of the Leaf Area Index

Abstract

The leaf area index (LAI) is a crucial variable in climate, ecological, and land surface modeling. However, the estimation of the LAI from coarse-resolution remote sensing data can be affected by the spatial scaling bias, which arises from the nonlinearity of retrieval models and the heterogeneity of the land surface. This study provides an algorithm named Arithmetic Mean and Geometric Mean (AM–GM) to correct the spatial scaling bias. It is established based on negative logarithmic functions and avoids second-order stationarity. In this algorithm, relationships are derived between the scaling bias of LAI and the arithmetic and geometric means of directional gap probability for two commonly used remote sensing models, the Beer–Lambert law and a semi-empirical transfer function, respectively. According to the AM–GM algorithm, the expression representing the model nonlinearity is derived and utilized for the analysis of LAI scaling bias. Furthermore, the AM–GM algorithm is simplified by a linear relationship, which is constructed between two quantities related to the directional gap probability between two specific resolutions. Two scenes simulated by the LargE-Scale remote sensing data and image Simulation framework (LESS) model and three sites are used to evaluate the proposed algorithm and analyze the scaling bias of LAI. The validation results show that the AM–GM algorithm provides accurate correction of LAI scaling bias. The analyses based on the AM–GM algorithm demonstrate that the scaling bias of LAI increases with the increase in the LAI value, with stronger surface heterogeneity and coarser spatial resolution. The validation results of the simplified AM–GM algorithm demonstrate that at the Sud-Ouest site, the absolute value of the bias for the estimated LAI decreases from 0.10, 0.22, 0.29, and 0.31 to 0.04, 0.01, 0.04, and 0.05 at 200 m, 500 m, 1000 m, and 1500 m resolutions, respectively. In conclusion, the proposed algorithm is effective in the analysis and correction of the scaling bias for coarse-resolution LAI.

1. Introduction

The leaf area index (LAI) is an essential variable in ecological, land surface, and climate modeling, as it explains the transfer of water, carbon, and energy between vegetation and the atmosphere [1,2]. The LAI is defined as one-half of the total green leaf area per unit of horizontal ground surface area [3]. In recent years, much research has applied coarse-resolution LAI (250 m–1000 m) to comprehend global vegetation change [4,5] and ecosystem productivity [6,7]. However, the accurate estimation of LAI at coarse resolution using remote sensing methods has been challenging due to the scaling effect [8].

The scaling effect is caused by model nonlinearity and surface heterogeneity. Most retrieval functions are nonlinear and built at a local scale (fine resolution) [9,10], which only interprets the surface heterogeneity at the corresponding scale [11,12]. When these functions are applied to a coarse spatial resolution, the scaling effect causes spatial scaling bias of LAI estimation [13,14]. Yin et al. reported that the scaling bias in LAI inversion can reach up to 26% [15]. Hence, scaling bias is an important consideration for the LAI estimated from coarse-resolution remote sensing data, especially on a heterogeneous surface [16,17,18].

To deal with this task, previous studies have proposed several methods. On the one hand, statistical regression relationships have been developed to calculate the spatial scaling bias [19,20]. For example, Zhao and Fan established a statistical regression relationship between a spatial heterogeneity coefficient and the scaling bias of LAI. However, empirical methods require external factors that limit their applicability. On the other hand, a series of universal methods have been provided [8,21,22,23,24,25,26]. For instance, the computational geometry method (CGM) sets the upper and lower envelopes of the retrieval model to achieve accurate parameter estimation at coarse resolutions [21]. The advantage of this method is that it is not affected by the continuity of the remote sensing model and does not require high-resolution remote sensing data. However, the true distribution of parameters is not considered. Chen et al. addressed the limitations of the CGM algorithm by introducing variable weights instead of fixed weights [26]. Nevertheless, the inherent flaw of the CGM method—its inability to express the scaling bias using a formula—has yet to be resolved. Similarly, machine learning models estimated parameters at coarse resolution based on the relationships built at a fine resolution, and the accuracy of the inversion parameter was controlled by the reference parameters at a coarse resolution [27]. Such methods are not suitable for analyzing scaling bias because a relevant formula cannot be established. The Taylor series expansion method (TSEM) is another universal method in which each term has a physical meaning [22]. This method is helpful for analyzing scaling bias. However, the arithmetic mean of the input variable for the LAI retrieval models is used to derive the second derivative term in the TSEM, which also ignores the true distribution of parameters [10,28]. This method also requires synchronous fine-resolution remote sensing data. Later, Garrigues et al. extended the TSEM algorithm to a multivariate transfer function but still left the problem unresolved [8]. Recently, several upscaling models are proposed based on fractal theory [23,24,25]. A number of studies have attempted to explain the fractal dimension [29], but no consensus has been reached. In summation, the TSEM algorithm is more suitable for analyzing and correcting the scaling bias, but the algorithm itself is limited by second-order stationarity and fine-resolution data.

TSEM is closely related to the model type [30]. In mathematics, scaling bias is referred to as negative Holder’s defect. When the model form is a negative logarithmic function, Holder’s defect can be accurately expressed by utilizing the arithmetic and geometric means of the input variables [31]. Hence, the negative logarithmic functions can be employed to develop a method for calculating and correcting the scaling bias. Two commonly used LAI inversion models with negative logarithmic form, the Beer–Lambert law [32] and a semi-empirical transfer function that establishes the relationship between the normalized difference vegetation index (NDVI) to LAI [33], are employed in this research. Notably, the equation for NDVI related to reflectance is nonlinear, leading to the scaling bias of NDVI on heterogeneous surfaces [8,34]. This bias not only affects the scaling bias calculation but also hinders the analysis of the scaling bias and the comprehensive understanding of the scaling effect.

In this study, an Arithmetic Mean and Geometric Mean (AM–GM) algorithm is proposed to deeply analyze and correct the LAI scaling bias. The algorithm aims to solve the inherent flaws of the TSEM. The practicality of this algorithm is enhanced by establishing a linear relationship between two specific resolutions, eliminating the need for fine-resolution remote sensing data.

2. Materials

In this study, data are obtained from both simulated scenes and in situ sites to validate the algorithm and analyze the scaling bias of LAI. The LargE-Scale remote sensing data and image Simulation framework (LESS) model were employed to simulate heterogeneous 3D scenes [35]. The in situ data were downloaded from the VAlidation of Land European Remote-sensing Instruments (VALERI) project (http://w3.avignon.inra.fr/valeri/, accessed on 15 March 2023). Two simulated scenes and three in situ sites are used to analyze and validate the proposed algorithm, and the land cover type includes cropland and forest. Four additional cropland sites are selected to simplify the AM–GM algorithm. At simulated scenes, the data are simulated at 1 m and then aggregated to 5 m, 9 m, 15 m, and 45 m. At the VALERI sites, the data are obtained at 20 m spatial resolution and aggregated to 200 m, 500 m, 1000 m, 1500 m, and 3000 m.

2.1. Simulated Data

Two 3D forest scenes are constructed using the LESS model in this paper, as shown in Figure 1. Both scenes have an area of 45 m × 45 m. The first scene consists of 581 trees. These trees are derived from 42 single-tree models, simulated by 6 different canopy shapes and 7 levels of leaf area index. Each tree’s position is determined based on field measurements. The LAI of the whole scene is 2.67. For more detailed information on the scene construction, please refer to Scene S4, which was constructed by Xu et al. [36]. Additionally, a virtual forest scene with stronger surface heterogeneity is constructed from the same single-tree models and different numbers of trees. The LAI of this virtual scene is 2.14. The spectral information of leaves and soil has been measured on the ground by Xu et al. [36]. The data simulations are performed with an observation zenith angle of 0°. Detailed parameter settings can be found in

Table 1. Parameter settings of the LESS model.

Sensor	Parameter Value	Optical Database	Parameter Value
Type	Orthographic	Larch Branch	Ground Measurements
Width (pixels)	45	Brown Loam	Ground Measurements
Height (pixels)	45	Larch Leaf	Ground Measurements
Samples (/pixel)	64	Terrain	Parameter Value
Spectral Bands	482:60, 561.5:57, 654.5:37, 865:28	Type	Plane
Image Format	Spectrum	XSize (m)	45
NoData Value	−1	YSize (m)	45
Width Extent (m)	45	BRDF Type	Lambertian
Height Extent (m)	45	Optical Property	Brown Loam
Four Components Product	Tick
Observation	Parameter Value	Objects	Parameter Value
View Zenith (°)	0	Single-Tree Models of Larch	Constructed by Xu et al. [36]
View Azimuth (°)	180
Sensor Height (m)	30
Illumination and Atmosphere	Parameter Value	Advanced	Parameter Value
Sun Zenith (°)	30	Minimum Iterations	5
Sun Azimuth (°)	90	Number of Cores	20
Sky-Type	SKY_TO_TOTAL
Sky-Percentage	0, 0

.

Over the simulated scenes, the directional gap probabilities are calculated from four component products, which are simulated using the LESS model:

$$p\left(\theta \right)=p_{is}+p_{ss}$$

(1)

where $p_{is}$ and $p_{ss}$ represent the proportions of illuminated soil and shaded soil within each pixel, respectively. At a spatial resolution of 1 m, the directional gap probabilities are simulated. The results are then aggregated to coarser spatial resolutions of 5 m, 9 m, 15 m, and 45 m.

2.2. Site Data

Seven VALERI sites are used and the data of these sites are downloaded from the VALERI website (Figure 2). To validate the algorithm and analyze the scaling bias, three sites are used. In particular, the vineyard dominates the Plan-de-Dieu site; Quercus ilex dominates the Puéchabon site; and the Sud-Ouest site includes nine crops, including corn, soya, sunflower, harvest wheat, etc. To simplify the AM–GM algorithm, four additional cropland sites covered by a variety of crops are chosen, including Les Alpilles, Barrax, Demmin, and Haouz.

Table 2. Detailed information of the VALERI sites.

Site Name	Land Cover	Day of Year	Image Year	Spatial Resolution	Site Size	Location
Plan-de-Dieu	Crops	181	2004	20 m	3 km × 3 km	44°11′N, 4°56′E
Puéchabon	Mediterranean Forests	163	2001	20 m	3 km × 3 km	43°43′N, 3°38′E
Sud-Ouest	Nine Crops	201	2002	20 m	3 km × 3 km	43°30′N, 1°14′E
Les Alpilles	Crops	204	2002	20 m	3 km × 3 km	43°48′N, 4°42′W
Barrax	Crops	195	2003	20 m	5 km × 3 km	39°40′N, 2°60′W
Demmin	Crops	164	2004	20 m	5 km × 3 km	53°53′N, 13°12′E
Haouz	Crops	73	2003	20 m	3 km × 3 km	31°39′N, 7°36′W

provides detailed information about these sites, including land cover, image acquisition time, site size, and location.

Near-infrared (${\rho }_{NIR}$) and red (${\rho }_R$) reflectances are obtained from SPOT images, which have a spatial resolution of 20 m. These data are aggregated to coarser resolutions (200 m, 500 m, 1000 m, 1500 m, and 3000 m), then NDVI and directional gap probability are computed at those spatial resolutions. The in situ LAI data are also downloaded to build the retrieval model at 20 m spatial resolution.

3. Methods

In this study, a novel algorithm called AM–GM is proposed to more accurately evaluate, analyze, and correct the scaling bias of LAI. This algorithm utilizes the advantage of negative logarithmic functions to avoid the second-order stationarity hypothesis, and this algorithm considers the scaling bias of NDVI. Moreover, the practicality of this algorithm is improved by eliminating the requirement for fine-resolution remote sensing data. In particular, the relationship between the LAI scaling bias and the arithmetic and geometric means of directional gap probability is investigated based on negative logarithmic functions. This relationship is exploited to derive the LAI scaling bias. In addition, the AM–GM algorithm uses Taylor’s theorem [22,31] to indirectly calculate the accurate values of the second derivative term, obtaining conditions for in-depth analysis of the model nonlinearity. Lastly, statistical relationships are established between two quantities related to directional gap probability between two specific scales for cropland. The flowchart depicting this process is shown in Figure 3.

3.1. Mathematical Theory of the AM–GM Algorithm

3.1.1. Holder’s Defect and AM–GM Inequality

In the field of mathematics, suppose that $X=\left\{x_1,x_2,\dots ,x_N\right\}$ is a random variable and model $f$ is twice continuously differentiable. The value estimated by the expected value of X based on the function $f$ is denoted by $f\left(E\left(X\right)\right)$, and the expected value of the function is denoted by $E\left(f\left(X\right)\right)$. Holder’s defect ($\Delta$) is used to express the discrepancy between the two terms [31]:

$$\Delta =E\left(f\left(X\right)\right)-f\left(E\left(X\right)\right)=\frac{1}{2}\mu V\mathrm{ar}\left(X\right)$$

(2)

where $Var\left(X\right)$ is the variance of X, and there exists a real value $\mu$ for the above formula. When $f$ is a negative logarithmic function, the AM–GM inequality associated with random variable X can be given an explicit formula following Holder’s defect [31].

$$\begin{array}{l}-\ln (\frac{x_A}{x_G})=-\ln \left(x_A\right)-\left(-\ln \left(x_G\right)\right)=-\ln \left(\frac{1}{N}{\displaystyle \sum _{k=1}^N{x}_k}\right)-\left(-\ln \left(\sqrt[N]{{\displaystyle \prod _{k=1}^N{x}_k}}\right)\right)\\ =-\ln \left(\frac{1}{N}{\displaystyle \sum _{k=1}^N{x}_k}\right)-\frac{1}{N}{\displaystyle \sum _{k=1}^N\left(-\ln \left(x_k\right)\right)}=f\left(E\left(X\right)\right)-E\left(f\left(X\right)\right)=-\Delta \end{array}$$

(3)

where $x_A$ is the arithmetic mean of $X$; $x_G$ is the geometric mean of $X$; $x_k\in X$; N is a natural number; $k$ is equal to 1 to N. Accordingly, Holder’s defect can be expressed by the AM–GM inequality, as follows:

$$\Delta =\ln \left(\frac{x_A}{x_G}\right)$$

(4)

3.1.2. Calculation of LAI Scaling Bias

The LAI retrieval model f is constructed at fine spatial resolution. Applying $f$ to fine-resolution data, then aggregating the resulting data to coarse resolution, yields the exact LAI ($LA{I}_{exa}$). The $LA{I}_{exa}$ is the reference LAI at coarse resolution because leaf area follows the law of conservation of matter [10]. The approximate LAI ($LA{I}_{app}$) is obtained by applying model f to coarse-resolution data, which are aggregated from fine-resolution data. The difference between $LA{I}_{app}$ and $LA{I}_{exa}$ is the spatial scaling bias [23], as shown in Figure 4.

In the field of remote sensing, according to Equations (2) and (4), the real scaling bias of the inverted parameter can be expressed as follows:

$$bia{s}_{Real}=LA{I}_{app}-LA{I}_{exa}=f\left(E\left(X\right)\right)-E\left(f\left(X\right)\right)=-\Delta =-\ln \left(\frac{x_A}{x_G}\right)$$

(5)

Next, two LAI retrieval models with negative logarithmic form are used to derive the scaling bias of LAI.

The Beer–Lambert law is a negative logarithmic function [37,38]:

$$LAI=-\frac{\cos \theta }{\mathsf{\Omega }G\left(\theta \right)}\ln \left(p\left(\theta \right)\right)$$

(6)

where $G\left(\theta \right)$ represents the leaf projection coefficient on a plane perpendicular to the observation direction $\theta$. $\theta$ is set to 30 degrees in Table 1. Assuming that foliage angles are distributed spherically, $G\left(\theta \right)$ is set to 0.5 [39]. $\mathsf{\Omega }$ is the foliage clumping index, which is set as a constant according to the vegetation type [40,41]. Here, fine-resolution directional gap probability $p\left(\theta \right)$ is measured based on geometric features, and there is no scaling effect [42]. According to Equations (5) and (6), the scaling bias of LAI is:

$$bia{s}_{AM-GM}=-\frac{\cos \theta }{\mathsf{\Omega }G\left(\theta \right)}\ln \left(\frac{A\left[p\left(\theta \right)\right]}{G\left[p\left(\theta \right)\right]}\right)$$

(7)

where $G\left[p\left(\theta \right)\right]$ and $A\left[p\left(\theta \right)\right]$ are the geometric and arithmetic means of fine resolution $p\left(\theta \right)$ within a coarse-resolution pixel, respectively.

The NDVI-LAI semi-empirical transfer function also has a negative logarithmic form [43],

$$LAI=-\frac{1}{K_{LAI}}\ln \left(\frac{NDVI-NDV{I}_{\max }}{NDV{I}_{\min }-NDV{\mathrm{I}}_{\max }}\right)$$

(8)

where $NDV{I}_{\max }$ is the asymptotic value of NDVI, which is obtained when LAI is greater than 8; $NDV{I}_{\min }$ is the NDVI value for bare soil; $K_{LAI}$ represents the extinction coefficient, which relates to vegetation spatial structure and vegetation canopy structure [33,44]. Here, the directional gap probability is obtained using Equation (9) [41]:

$$p\left({\rho }_{NIR},{\rho }_R\right)=\frac{\frac{{\rho }_{NIR}-{\rho }_R}{{\rho }_{NIR}+{\rho }_R}-NDV{I}_{\max }}{NDV{I}_{\min }-NDV{\mathrm{I}}_{\max }}$$

(9)

where ${\rho }_{NIR}$ and ${\rho }_R$ represent near-infrared and red reflectances, respectively. To avoid introducing the scaling bias of NDVI into the directional gap probability calculation at coarse resolution, we use the aggregated reflectance rather than the aggregated NDVI. According to Equations (5) and (8), the scaling bias of LAI is:

$$bia{s}_{\rho \_AM-GM}=-\frac{1}{K_{LAI}^{}}\ln \left(\frac{p\left(A\left[{\rho }_{NIR}^{}\right],A\left[{\rho }_R^{}\right]\right)}{G\left[p\left({\rho }_{NIR}^{},{\rho }_R^{}\right)\right]}\right)$$

(10)

where $G\left[p\left({\rho }_{NIR}^{},{\rho }_R^{}\right)\right]$ is the geometric mean of fine resolution $p\left({\rho }_{NIR}^{},{\rho }_R^{}\right)$ within a coarse-resolution pixel, and $p\left(A\left[{\rho }_{NIR}^{}\right],A\left[{\rho }_R^{}\right]\right)$ is the directional gap probability at coarse resolution estimated based on the arithmetic mean of the fine-resolution reflectance.

3.2. Calculate Factor μ

The second derivative term $\mu$ absorbs the impact of all higher-order moments associated with $X$ in Equation (2). $p\left(\theta \right)$ is measured according to the canopy geometry [45], the value at coarse resolution is equal to the arithmetic mean of fine-resolution directional gap probability $p\left(\theta \right)$. The value of $\mu$ can be derived directly from Equation (2). It can also be obtained using the second derivative. Therefore, the $\mu$ value based on Equation (6) is shown as follows:

$${\mu }_{AM-GM}=-\frac{2bia{s}_{AM-GM}}{Var\left(X\right)}=\frac{1}{\frac{\mathsf{\Omega }G\left(\theta \right)}{\cos \theta }\times c_k{}^2}$$

(11)

where $c_k$ is related to $p\left(\theta \right)$.

When the directional gap probability calculation involves a nonlinear model, such as Equation (8), the value of $p\left(A\left[{\rho }_{NIR}^{}\right],A\left[{\rho }_R^{}\right]\right)$ is not equal to the value of $A\left[p\left({\rho }_{NIR}^{},{\rho }_R^{}\right)\right]$ (the arithmetic mean of fine-resolution $p\left({\rho }_{NIR}^{},{\rho }_R^{}\right)$). The first-order derivative term cannot be eliminated. Thus, Equation (2) is modified as follows:

$$\Delta =\frac{1}{2}{\mu }_{\rho \_AM-GM}Var(X)+f^{\prime }\left(x^M\right){\displaystyle \sum _{k=1}^N\frac{\left(x_k-x^M\right)}{N}}$$

(12)

where $x_k$ is the kth value of $p\left({\rho }_{NIR}^{},{\rho }_R^{}\right)$ within a coarse-resolution pixel; $x^M$ represents the directional gap probability at coarse resolution, i.e., $p\left(A\left[{\rho }_{NIR}^{}\right],A\left[{\rho }_R^{}\right]\right)$; ${\mu }_{\rho \_AM-GM}$ is the second derivative term and $Var\left(X\right)$ is the mean square error. According to Equations (10) and (12), the value of $\mu$ based on Equation (8) can be expressed as follows:

$${\mu }_{\rho \_AM-GM}=-\frac{2\left(bia{s}_{\rho \_AM-GM}+f^{\prime }\left(x^M\right)\frac{{\displaystyle \sum _{K=1}^N\left(x_k-x^M\right)}}{N}\right)}{ Var(X)}=\frac{1}{k_{LAI}\times c_k{}^2}$$

(13)

where $c_k$ is related to $p\left({\rho }_{NIR}^{},{\rho }_R^{}\right)$.

3.3. Simplify the AM–GM Algorithm

The application of the AM–GM algorithm is limited by its use of fine-resolution remote sensing data. In this section, this algorithm was simplified. According to Equation (A4) in Appendix A and Equation (10), the following formula can be derived:

$$bia{s}_{\rho \_AM-GM}=LA{I}_{app}\times \left(1-\frac{\ln \left(G\left[p\left({\rho }_{NIR}^{},{\rho }_R^{}\right)\right]\right)}{\ln \left(p\left(A\left[{\rho }_{NIR}^{}\right],A\left[{\rho }_R^{}\right]\right)\right)}\right)$$

(14)

where $\ln \left(p\left(A\left[{\rho }_{NIR}^{}\right],A\left[{\rho }_R^{}\right]\right)\right)$ represents the logarithm of the directional gap probabilities at coarse resolution and $\ln \left(G\left[p\left({\rho }_{NIR}^{},{\rho }_R^{}\right)\right]\right)$ represents the logarithm of the geometric mean of fine resolution $p\left({\rho }_{NIR}^{},{\rho }_R^{}\right)$ within a coarse-resolution pixel.

A linear relationship between $\ln \left(p\left(A\left[{\rho }_{NIR}^{}\right],A\left[{\rho }_R^{}\right]\right)\right)$ and $\ln \left(G\left[p\left({\rho }_{NIR}^{},{\rho }_R^{}\right)\right]\right)$ was found between two specific resolutions; the relationship will be described in further detail in the Results section. Equation (14) can be simplified according to the linear relationship, as follows:

$$bia{s}_{\rho \_AM-GM}=LA{I}_{app}\times \left(\frac{b}{\ln \left(p\left(A\left[{\rho }_{NIR}^{}\right],A\left[{\rho }_R^{}\right]\right)\right)}-a\right)$$

(15)

where parameters $a$ and $b$ are constants between two specific resolutions. According to Equation (15), the scaling bias of LAI can be calculated based on the coarse resolution data.

4. Results

4.1. Validation Results of the AM–GM Algorithm

Firstly, this study validates the correction effectiveness of the AM–GM algorithm when the Beer–Lambert law is the LAI retrieval model. The computation of LAI scaling bias based on the Beer–Lambert law is derived in Equation (7). Figure 5 and

Table 3. LAI scaling bias before and after correction at 500 m spatial resolution.

Site Name	Before Correction		After Correction
	RMSE	Bias	RMSE	Bias
Genhe	0.71	−0.67	0.00	0.00
Virtual	0.74	−0.62	0.00	0.00
Plan-de-Dieu	0.05	−0.03	0.00	0.00
Puéchabon	0.19	−0.14	0.00	0.00
Sud-Ouest	0.25	−0.22	0.00	0.00

show the magnitude of the scaling bias before and after correction. Before the correction, the retrieval $LA{I}_{app}$ at coarse resolution exhibits significant negative scaling biases (bias = −0.67 and bias = −0.62). After correction, the bias and RMSE values between $LA{I}_{cor}$ and $LA{I}_{exa}$ are both reduced to almost 0. Furthermore, the spatial distribution of the corrected values $LA{I}_{cor}$ coincides with $LA{I}_{exa}$. These results indicate that the AM–GM algorithm is effective when the Beer–Lambert law is the retrieval model.

We also assess the correction effectiveness when the NDVI-LAI semi-empirical transfer function is the LAI retrieval model. The equation for the calculation of LAI scaling bias is shown in Equation (10). Figure 6 and Table 3 show that the scaling biases are small at the Plan-de-Dieu site, with RMSE = 0.05 and bias = −0.03; moderate at the Puéchabon site, with RMSE= 0.19 and bias = −0.14; and high at the Sud-Ouest site, with RMSE = 0.25 and bias = −0.22. After correction, the values of RMSE and bias between $LA{I}_{cor}$ and $LA{I}_{exa}$ are also deduced to almost 0, and all spots are clustered on the 1:1 line. Meanwhile, the spatial distributions of $LA{I}_{cor}$ coincide with $LA{I}_{exa}$. These results indicate that the AM–GM algorithm is effective when the NDVI-LAI semi-empirical transfer function is the retrieval model.

4.2. Analysis of LAI Scaling Bias

Two figures are used to analyze the LAI scaling bias based on three aspects: spatial resolution, model nonlinearity, and surface heterogeneity. In Figure 7, we calculate the average of LAI scaling bias for each scene and site using the AM–GM algorithm at multiple scales. In Figure 8, the values of $\mu$ and $Var\left(X\right)$ are calculated for each coarse pixel based on Equations (11) and (13).

• Spatial resolution.

When the spatial resolution for LAI estimation becomes coarser, the scaling bias of LAI increases. The scaling bias of LAI between two close resolutions may also be large. In Figure 7, as the spatial resolution moves from 1 to 45, the value of LAI scaling bias decreases and the LAI scaling bias increases. For the Beer–Lambert law, the scaling bias of $LA{I}_{app}$ at 5 m spatial resolution relative to $LA{I}_{exa}$ estimated at 1 m spatial resolution is large, as shown in Figure 7a.

• Model nonlinearity.

When the LAI value increases, the degree of model nonlinearity also rises. This is because the model nonlinearity is related to $1/c_k{}^2$ in Equations (11) and (13), and the value of $c_k$ is related to directional gap probability. When the value of $c_k$ for a coarse pixel is small, the directional gap probability of this pixel is also small, and the LAI value is large [46]. For example, in Figure 6, the exact LAI values for row 2 to 4 in the fourth column of the Sud-Ouest site are relatively large. Correspondingly, in Figure 8, the values of ${\mu }_{\rho \_AM-GM}$ calculated by the AM–GM algorithm for these pixels are also relatively high.

• Surface heterogeneity.

When vegetation in a coarse pixel has different growing seasons, the degree of surface heterogeneity is high. The Sud-Ouest site includes nine crops with various growing seasons, while both the Plan-de-Dieu and Puéchabon sites have a relatively singular growing season. Therefore, Figure 8 shows that the $Var\left(X\right)$ values of the Sud-Ouest site are larger than those of the other two sites.

4.3. Scaling Bias Calculated by the Simplified Algorithm

Although the AM–GM algorithm can correct the LAI scaling bias more accurately, it requires fine-resolution remote sensing data, which limits the applicability of this algorithm. The simplified AM–GM algorithm does not use fine-resolution data.

Figure 9 shows that linear relationships between the two quantities $\ln \left(p\left(A\left[{\rho }_{NIR}^{}\right],A\left[{\rho }_R^{}\right]\right)\right)$ and $\ln \left(G\left[p\left({\rho }_{NIR}^{},{\rho }_R^{}\right)\right]\right)$ are observed at four croplands between the fine and multiple resolutions (200 m, 500 m, 1000 m, and 1500 m), respectively. Due to the limited sample size, this relationship is not built at the spatial resolution of 3 km. The parameters $a$ and $b$ for simplified AM–GM algorithm in Equation (15) are shown in

Table 4. The parameters $a$ and $b$ at multiple scales for the simplified AM–GM algorithm.

	Spatial Resolution	200 m	500 m	1000 m	1500 m
Parameters
$a$		0.052	0.089	0.056	0.043
$b$		0.011	0.022	0.063	0.081

.

The simplified AM–GM algorithm corrected the scaling bias of LAI at multiple scales over the Sud-Ouest site. Figure 10 shows that after correction, the values of $LA{I}_{cor}$ are close to the values of $LA{I}_{exa}$. The absolute values of bias at the Sud-Ouest site decrease from 0.10, 0.22, 0.29, and 0.31 to 0.04, 0.01, 0.04, and 0.05 at 200 m, 500 m, 1000 m, and 1500 m, respectively. Moreover, Figure 11 shows that the spatial distributions of $LA{I}_{cor}$ are well corrected, particularly at spatial resolutions of 1000 m and 1500 m. Therefore, the simplified AM–GM algorithm, which does not use fine-resolution images, is feasible.

5. Discussion

5.1. Algorithm Comparison

Prior work has documented the flaws in the TSEM algorithm resulting from the second-order stationary hypothesis, which assumes that the mean of a random function is a constant [30]. The AM–GM algorithm considered the higher-order moments, i.e., the skewness that measures the asymmetry, and the kurtosis that measures the heaviness of the tails of distribution without discarding [31]. We compared the values of $\mu$ calculated using the AM–GM and TSEM algorithms in

Table 5. Model nonlinearity factor ($\mu$) and input variable for the AM–GM and TSEM algorithms at 500 m resolution.

Site Name	$\mathit{\mu }$		Input Variable
	RMSE	Bias	RMSE	Bias
Plan-de-Dieu	5.56	−3.99	0.01	−0.00
Puéchabon	5.48	−4.58	0.02	−0.01
Sud-Ouest	13.43	7.16	0.02	0.01

. The results show bias in $\mu$ up to 7.16. However, the numerical discrepancies of the input variable used for $\mu$ calculation are less than 0.02. It is inferred that the $\mu$ value is sensitive to the input variable. On the heterogeneous surfaces, the skewness and kurtosis of the random variable may be large, and the AM–GM algorithm may exhibit a better performance for the computation of $\mu$ and LAI scaling bias.

In the improved CGM, the relative errors decreased from 3.35%, 11.01%, and 19.62% to an average of 0.28%, 1.48%, and 5.16%, respectively, at kilometer scales [26]. Jiang et al. indicated that the LAI scaling bias decreased to around 0.27 by the TSEM, Wavelet-Fractal (WF), and Fractal theory (TT) [30]. We determine that the AM–GM algorithm decreases the LAI scaling bias to almost zero. The simplified AM–GM algorithm decreases relative error from 16.07% to 4.3% at 1 km resolution. The LAI scaling bias decreases from 0.10, 0.22, 0.29, and 0.31 to 0.04, 0.01, 0.04, and 0.05 at 200 m, 500 m, 1000 m, and 1500 m resolutions, respectively. On the heterogeneous surfaces, the correction effect is comparable or even better. The simplified algorithm, TSEM, WF, or TT, does not perform well at homogeneous sites.

5.2. The Influence of NDVI Aggregation on the AM–GM Algorithm

To minimize errors of LAI estimation, care should be taken in relation to the scale issue of NDVI. In the past, studies usually used the vegetation index (VI) to calculate the scaling bias of LAI [30]. However, most of the equations for VIs are nonlinear, including the NDVI equation related to reflectance. The scaling bias of NDVI complicates the evaluation of the LAI scaling bias on the heterogeneous surface.

Table 6. Influence of NDVI aggregation on the AM–GM algorithm at 500 m resolution.

Site Name	Considering the Scaling Bias of NDVI		Ignoring the Scaling Bias of NDVI
	RMSE	Bias	RMSE	Bias
Plan-de-Dieu	0.00	0.00	0.18	0.01
Puéchabon	0.00	0.00	0.49	0.04
Sud-Ouest	0.00	0.00	0.61	−0.08

shows that the aggregated NDVI results in inaccurate estimates of LAI scaling bias based on the AM–GM algorithm. Additionally, the bias may be a positive or a negative value. In particular, several approaches have been provided to reconstruct time-series data for NDVI using both fine- and coarse-resolution NDVI data. The discussion in 5.2 may indicate that if the time-series data are used for LAI estimation, the scaling bias of NDVI should be considered during data reconstruction.

5.3. Limitation of the AM–GM Algorithm

This algorithm is only suitable for negative logarithmic functions because Holder’s defect provides an accurate expression of the AM–GM inequality for negative logarithmic functions [31]. This limitation restricts the application of the AM–GM algorithm. In addition, this study corrected LAI scaling bias due to the spatial heterogeneity of vegetation growth but did not consider the model difference caused by vegetation types. In the future, we plan to further improve this algorithm for use with mixed land surfaces.

6. Conclusions

We found that scaling bias in remote sensing can be accurately expressed by the AM–GM inequality for negative logarithmic function. Based on this relationship, the AM–GM algorithm was proposed. This algorithm avoids the influence of the second-order stationary hypothesis. In addition, the aggregation of NDVI, which was ignored in previous research about scaling bias correction, is considered in the AM–GM algorithm. In this study, we validated the AM–GM algorithm at two simulated scenes and three sites, proving its effectiveness in the correction and analysis of the LAI scaling bias at coarse resolution.

In addition, upon analyzing the second derivative term, we identified a positive correlation between the LAI value and the degree of model nonlinearity. This means that as the LAI value increases, the nonlinearity in the model also increases. It was also observed that when the spatial resolution is coarse and the surface heterogeneity increases, there is a greater scaling bias in the LAI estimation. Adjusting for the LAI scaling bias was critical in the above three cases.

Furthermore, we identified a linear statistical relationship between two quantities related to directional gap probability between two specific scales. The AM–GM algorithm is simplified by the relationship. The results indicate that the simplified algorithm can practically correct the scaling bias without the requirement of fine-resolution remote sensing data. In conclusion, the correction effect of this algorithm is reliable. Future work includes the improvement of this algorithm on mixed land surfaces.