Evaluation of Three Techniques for Correcting the Spatial Scaling Bias of Leaf Area Index

The correction of spatial scaling bias on the estimate of leaf area index (LAI) retrieved from remotely sensed data is an essential issue in quantitative remote sensing for vegetation monitoring. We analyzed three techniques, including Taylor’s theorem (TT), Wavelet-Fractal technique (WF), and Fractal theory (FT), for correcting the scaling bias of LAI with empirical models in different functions (i.e., power, exponential, logarithmic and polynomial) on both simulated data and a real dataset over a cropland site. The results demonstrated that the scaling bias became greater when the model non-linearity increased. The spatial heterogeneity, which was characterized by the class-specific proportion, the between-class spectral difference and the number of classes within each coarse pixel, was found to be the primary factor in the scaling effect. These factors influenced the scaling effect collectively and existed dependently. With the RMSE less than 0.3 × 10−6 m2/m2, TT was suggested for the correction with a polynomial LAI-NDVI functions. WF was preferred for neighboring scales rather than continuous scales. FT was not recommended for correcting the scaling bias caused by the significant non-linearity in LAI estimation models. This study illustrates the main causes of the scaling effect and provides a reference of technique selection for scaling bias correction to improve the application of remotely sensed estimates.


Introduction
Leaf area index (LAI), defined as half the total foliage area per unit ground horizontal surface area [1], is an essential biophysical property of vegetation for environmental assessment [2,3], crop growth monitoring [4,5] and yield prediction [6][7][8].The accurate monitoring of LAI on multiple scales has become a key technique in those applications.However, the description of surface features based on quantitative remote sensing must have a multi-scale problem, scaling effect [9].Due to the difference of spatial scale (the resolution of an image), the properties or laws of LAI obtained from one scale would change on another scale even with the same observation conditions [10][11][12][13].Without correcting the scaling bias, the relative errors of LAI estimates from SPOT5 at a large scale were 2% in woodland and 7% in a crop-water area [14].A case study of the middle reach of the Heihe River Basin showed that the scaling bias in LAI retrieval was 26% if the scaling effect was ignored [15].Wu et al. [16] found that the relative scaling bias of LAI could be 55% in an uncorrected, large-scale mixed pixel.LAI of mixed forests in Canada were overestimated by approximately 200% based on the MODIS LAI product [17].The maximum relative scaling error of 400% would be introduced in the mixing pixel composed of vegetation and water [18].Huang et al. [19] assimilated the LAI estimates at various resolutions into the WOFOST model to monitor heavy metal stress in rice.They found the uncertainties from 16 m to 256 m changed from 6% to 90% compared to the original GF-1 image with 8-m spatial resolution for the area with severe metal stress.Therefore, it is crucial to correct the scaling bias of LAI for improving the authenticity and availability of quantitative remote sensing monitoring.
Before correcting the scaling bias of remotely sensed LAI, the physical mechanism of the scaling effect should be investigated.Many studies have focused on the spatial scaling effect of LAI [10,11,[20][21][22][23][24][25][26][27][28][29].For continuous canopies, Milne and Cohen [27] analyzed the variables including LAI at multiple scales, and Xu et al. [26] found that the nonlinear relationship of the averaged vegetation area fraction within mixed pixel was the main cause of scaling effect.For discrete canopies, Fan et al. [10] explored the reason for spatial scaling effect of LAI based on the directional second-derivative method of effective LAI retrieval.For some typical crops, such as wheat and corn, Zhu et al. [29] analyzed the scaling effect of LAI at maturity stage using a two-layer canopy reflectance model.For reeds, Chen et al. [28] discussed the main cause of LAI scaling effect based on Thematic Mapper and Moderate-resolution Imaging Spectroradiometer data.It has been concluded that the primary causes of scaling effect are the non-linearity of LAI estimation model and the spatial heterogeneity of land surface.However, previous studies only focused on a specific condition (i.e., specific vegetation, particular remotely sensed data).How the model non-linearity and spatial heterogeneity affect the scaling effect remains unclear.
To correct the scaling bias, some statistical regression relationships among LAIs on different scales have been established [30][31][32][33][34]. Since the empirical regression method is generally susceptible to external factors, it has great limitations in applicability.Therefore, several generic algorithms were established based on mathematical and geostatistical theory [11,24,35,36].In 1992, a spatialization model named CGM (computational geometry method) was proposed to reduce the scaling bias [35].For the application of CGM, the retrieval model should follow a uniform distribution of the radiometric variables between the convex and concave function.However, the assumption was argued to be inappropriate within a moderate resolution pixel [11], which led to significant overestimations of scaling bias [37].In consideration of this uncertainty, Hu and Islam [24] proposed the Taylor series expansion method.To improve the practical utility of this method, Garrigues et al. [11] extended the Taylor expansion to correct the scaling bias by a multivariate transfer function based on the Taylor's theorem (TT) of the relationship between variance, covariance and the mean values of parameters on the pixel scale.However, the multivariate model strongly relies on the second-order stationarity hypothesis.TT might fail to correct the scaling bias when the hypothesis does not hold.
To cope with the problem, Jiang et al. [38] proposed the Wavelet-Fractal technique (WF).The fractal theory was first introduced by Mandelbrot [39] to decompose a function or signal to infinite dimensional space.The exponent of each item x indicates the dimension of the space (such as x 2 representing a two-dimensional space).The coefficient of x is the energy of the function or signal in the corresponding dimensional space.Based on the wavelet transform, each component is orthogonal.That means the energy in different dimensional space has no intersection and the components from wavelet transform are more representative than those from TT.By combining these two methods, it is found that the scaling bias and the high frequency components decomposed by Haar wavelet were fractal relationship.The scaling bias of LAI for both univariate and bivariate models were effectively corrected on multi-scales [38].
Given the fractal characteristics (self-similarity and fractional dimensionality) of remotely sensed images, the fractal relationship is one of the transformation tendencies that can be used to describe the scaling effect.Zhang et al. [12] defined the information fractal of quantitative remote sensing products to describe the scaling transformation law of LAIs on different scales, but it could not be applied to scaling effect correction.Subsequently, an image-based LAI scaling transfer model was established [40], which could only correct the scaling bias of the entire image.With extension from the image-based method, Wu et al. [41] developed a pixel-based model based on the Fractal Theory (FT) to correct the LAI scaling effect.
Previous studies have been dedicated to answer the key questions on the scaling effect of LAI, which include: What are the fundamental causes of the scaling effect?What kind of methods could be effectively used to correct the scaling bias?Many of them have focused on the scaling effect from different perspectives.However, the analysis of LAI scaling effect has been mostly confined to a specific crop type or particular kinds of remotely sensed data.It has been difficult to determine the underlying mechanism and interactions among the process of the scaling effect.Moreover, to the best of our knowledge, no research in the literature has provided a comprehensive comparison of these methods for correcting LAI scaling bias.To address these research gaps, this study uses numerical experiments with a large amount of randomly simulated data to explore the fundamental causes of scaling effect and their influences.It presents a comparative assessment of three representative techniques (TT, WF, and FT) applied to simulated datasets and remotely sensed observations over a cropland site.The main objectives were: (1) to investigate the influences of landscape spatial heterogeneity and the non-linear estimation of LAI on the scaling effect; and (2) to evaluate the performance and robustness of the three scaling effect correction techniques in various conditions.

Principle of Scaling Effect
Figure 1 indicates two different processes of retrieving a biophysical parameter through remote sensing for the same ground object in the same area using the same model f that is a function relating x to y. Y exa is the result obtained by retrieving at the high spatial resolution using f and then aggregating the values from four pixels.Y app is computed from the direct application of f to the average value of four pixels (X).In the application of LAI estimation, x represents the vegetation index (VI) of a high resolution pixel with spatial heterogeneity.If x 1 = x 2 = VI a and x 3 = x 4 = VI b , then VI m is the average representing the VI value of the coarse resolution pixel.f is the non-linear estimation model relating LAI to VI. LAI exa and LAI app are the exact and approximated values on the coarse scale, respectively.It is obvious that the estimation process would be scale invariant if the land surface is homogeneous or the model is linear.With the effects of both spatial heterogeneity and model non-linearity, LAI app (from Path_2, "aggregation-estimation") overestimates the exact value of LAI (LAI exa from Path_1, "estimation-aggregation"), which results in the scaling effect and introduces the scaling bias.Therefore, the non-linearity of LAI estimation model and the landscape spatial heterogeneity are the sufficient and necessary conditions of scaling effect.not be applied to scaling effect correction.Subsequently, an image-based LAI scaling transfer model was established [40], which could only correct the scaling bias of the entire image.With extension from the image-based method, Wu et al. [41] developed a pixel-based model based on the Fractal Theory (FT) to correct the LAI scaling effect.
Previous studies have been dedicated to answer the key questions on the scaling effect of LAI, which include: What are the fundamental causes of the scaling effect?What kind of methods could be effectively used to correct the scaling bias?Many of them have focused on the scaling effect from different perspectives.However, the analysis of LAI scaling effect has been mostly confined to a specific crop type or particular kinds of remotely sensed data.It has been difficult to determine the underlying mechanism and interactions among the process of the scaling effect.Moreover, to the best of our knowledge, no research in the literature has provided a comprehensive comparison of these methods for correcting LAI scaling bias.To address these research gaps, this study uses numerical experiments with a large amount of randomly simulated data to explore the fundamental causes of scaling effect and their influences.It presents a comparative assessment of three representative techniques (TT, WF, and FT) applied to simulated datasets and remotely sensed observations over a cropland site.The main objectives were: (1) to investigate the influences of landscape spatial heterogeneity and the non-linear estimation of LAI on the scaling effect; and (2) to evaluate the performance and robustness of the three scaling effect correction techniques in various conditions.

Principle of Scaling Effect
Figure 1 indicates two different processes of retrieving a biophysical parameter through remote sensing for the same ground object in the same area using the same model f that is a function relating x to y. Yexa is the result obtained by retrieving at the high spatial resolution using f and then aggregating the values from four pixels.Yapp is computed from the direct application of f to the average value of four pixels (X).In the application of LAI estimation, x represents the vegetation index (VI) of a high resolution pixel with spatial heterogeneity.If x1 = x2 = VIa and x3 = x4 = VIb, then VIm is the average representing the VI value of the coarse resolution pixel.f is the non-linear estimation model relating LAI to VI. LAIexa and LAIapp are the exact and approximated values on the coarse scale, respectively.It is obvious that the estimation process would be scale invariant if the land surface is homogeneous or the model is linear.With the effects of both spatial heterogeneity and model non-linearity, LAIapp (from Path_2, "aggregation-estimation") overestimates the exact value of LAI (LAIexa from Path_1, "estimation-aggregation"), which results in the scaling effect and introduces the scaling bias.Therefore, the non-linearity of LAI estimation model and the landscape spatial heterogeneity are the sufficient and necessary conditions of scaling effect.is the value of each pixel at the original high resolution and X is their average.Suppose x 1 = x 2 = VI a , x 3 = x 4 = VI b , then X is VI m representing the vegetation index value of the coarse spatial resolution pixel, and y i is the corresponding LAI value.The subscript exa means the exact value and app represents the approximated value of the coarse resolution pixel.

Real Data
The study area for the real dataset was located within the city of Huai'an, Jiangsu Province and in the Yangtze-Huai Plain, one of the major agricultural regions in China (119 • 6 16 E, 33 • 23 45 N) (Figure 2).The crops in this area are often grown in a rice-wheat double cropping system.The dataset was acquired in the winter wheat season of 2014-2015 when the cultivated winter wheat cultivars were Huai-35 and Yangfu-4.To acquire field measurement data over the study area, five times of sampling were performed during the growth period of winter wheat (early jointing, jointing, booting, heading and filling stages) from March to May 2015.Throughout the study area, 24 samples were collected within an area of 2 m × 2 m each.Three subsample plots (0.5 m × 0.5 m) were enclosed in each sample.The plots were relatively homogeneous and distributed as shown in Figure 2c.The spectral measurements expressed as NDVI (Normalized Difference Vegetation Index) were measured with GreenSeeker RT100 (NTech Industries, Ukiah, CA, USA).On the same days, LI-3000C Portable Area Meter (LI-COR Biosciences, Lincoln, NE, USA) was used to measure LAI with 1 mm 2 resolution by destructive sampling of 30 plant individuals.NDVI and LAI of all subsamples from each sample were collected and then averaged to represent the sample means.and yi is the corresponding LAI value.The subscript exa means the exact value and app represents the approximated value of the coarse resolution pixel.

Real Data
The study area for the real dataset was located within the city of Huai'an, Jiangsu Province and in the Yangtze-Huai Plain, one of the major agricultural regions in China (119°6′16″E, 33°23′45″N) (Figure 2).The crops in this area are often grown in a rice-wheat double cropping system.The dataset was acquired in the winter wheat season of 2014-2015 when the cultivated winter wheat cultivars were Huai-35 and Yangfu-4.To acquire field measurement data over the study area, five times of sampling were performed during the growth period of winter wheat (early jointing, jointing, booting, heading and filling stages) from March to May 2015.Throughout the study area, 24 samples were collected within an area of 2 m × 2 m each.Three subsample plots (0.5 m × 0.5 m) were enclosed in each sample.The plots were relatively homogeneous and distributed as shown in Figure 2c.The spectral measurements expressed as NDVI (Normalized Difference Vegetation Index) were measured with GreenSeeker RT100 (NTech Industries, Ukiah, CA, USA).On the same days, LI-3000C Portable Area Meter (LI-COR Biosciences, Lincoln, NE, USA) was used to measure LAI with 1 mm 2 resolution by destructive sampling of 30 plant individuals.NDVI and LAI of all subsamples from each sample were collected and then averaged to represent the sample means.The remotely sensed image used in this study was a scene of the GF-1, the first series of Chinese High Resolution Earth Observation System launched on 26 April 2013.The sensor is equipped with two 2-m resolution panchromatic/8-m multispectral cameras (PMS) and four 16-m wide-field-of-view cameras (WFV).The GF-1 image was acquired on 22 April 2015 quasi simultaneously with the collection of field measurements on 27 April at heading stage.To reduce the estimation bias due to the time gap between satellite data and field measurements acquisitions, we used the field-measured LAI and GreenSeeker NDVI over all of the major growth stages of winter wheat to derive four basic statistical functions (Figure 3, Table 1).This ensures the universality of the LAI-NDVI models for the simulation experiments.The radiometric correction was performed for all the panchromatic and multispectral data.The atmospheric correction was taken by using the Fast Line-of-Sight Atmospheric Analysis of Spectral Hyper-cubes (FLAASH) model to obtain the reflectance in each multispectral band.The geometric correction was performed for all the images through the RPC The remotely sensed image used in this study was a scene of the GF-1, the first series of Chinese High Resolution Earth Observation System launched on 26 April 2013.The sensor is equipped with two 2-m resolution panchromatic/8-m multispectral cameras (PMS) and four 16-m wide-field-of-view cameras (WFV).The GF-1 image was acquired on 22 April 2015 quasi simultaneously with the collection of field measurements on 27 April at heading stage.To reduce the estimation bias due to the time gap between satellite data and field measurements acquisitions, we used the field-measured LAI and GreenSeeker NDVI over all of the major growth stages of winter wheat to derive four basic statistical functions (Figure 3, Table 1).This ensures the universality of the LAI-NDVI models for the simulation experiments.The radiometric correction was performed for all the panchromatic and multispectral data.The atmospheric correction was taken by using the Fast Line-of-Sight Atmospheric Analysis of Spectral Hyper-cubes (FLAASH) model to obtain the reflectance in each multispectral band.The geometric correction was performed for all the images through the RPC (rational polynomial coefficient) ortho-rectification with a DEM (digital elevation model) file.Finally, we used the Nearest Neighbor Diffusion (NNDiffuse) pan sharpening algorithm to merge the panchromatic and multispectral images into a 2-m resolution fused image set.All of these operations were performed in ENVI 5.2.According to Figure 1, the high-spatial resolution data GF-1 with 2 m (1024 × 1024 pixels) was up-scaled to six large scales: 4 m (512 × 512 pixels), 8 m (256 × 256 pixels), 16 m (128 × 128 pixels), 32 m (64 × 64 pixels), 64 m (32 × 32 pixels), and 128 m (16 × 16 pixels).A land-use map in the study area was obtained using the Maximum Likelihood Classification with 24,517 pixels representing the five classes of land use classes built-up, water, bare soil, winter wheat and other vegetation (Figure 4a).The overall accuracy of the classification was 90.1% and the Kappa coefficient was 0.87.The land use map was used to mask all the pixels that were not classified as winter wheat (Figure 4c).A land-use map in the study area was obtained using the Maximum Likelihood Classification with 24,517 pixels representing the five classes of land use classes built-up, water, bare soil, winter wheat and other vegetation (Figure 4a).The overall accuracy of the classification was 90.1% and the Kappa coefficient was 0.87.The land use map was used to mask all the pixels that were not classified as winter wheat (Figure 4c).A land-use map in the study area was obtained using the Maximum Likelihood Classification with 24,517 pixels representing the five classes of land use classes built-up, water, bare soil, winter wheat and other vegetation (Figure 4a).The overall accuracy of the classification was 90.1% and the Kappa coefficient was 0.87.The land use map was used to mask all the pixels that were not classified as winter wheat (Figure 4c).

Simulated Data
The model non-linearity and spatial heterogeneity are the main factors of scaling effect.The spatial heterogeneity is manifested as the mixing or purity of the interior, including the class-specific proportion, between-class spectral difference, and the number of classes within coarse pixels.To explore the influences and interactions of these factors on scaling effect, numerical simulation experiments were conducted based on the variable control principle (Table 2).

•
Experiment 1, model non-linearity.We designed a simulation experiment encompassing 10,000 groups of matrices in the dimension of 16 × 16.These matrices comprised randomly generated numbers from 0 to 1 as NDVI values and represented subpixel images of the coarse pixels.We used the exponential function model (f : LAI = a NDVI ) and the base number a to characterize the variation in linearity of the LAI estimation model.For random pixels of each group, the results were calculated by both "estimation-aggregation" (Path_1) and "aggregation-estimation" (Path_2) as shown in Figure 1.The difference was the scaling bias.

•
Experiment 4, the number of classes.Three classes were defined as NDVI 1 , NDVI 2 and NDVI 3 , with the value of 0.01, 0.5 and 0.9 respectively.Firstly, calculating the scaling bias of pixels with two classes (NDVI 1 &NDVI 2 , NDVI 1 &NDVI 3 , and NDVI 2 &NDVI 3 ).Then, we added the third class to this coarse pixel and maintained the equal proportion of each class (p i = 33.33%).Such that, the number of classes within the coarse pixel increased to three (NDVI 1 &NDVI 2 &NDVI 3 ).
In Experiment 2-4, four types of estimation models (power function, exponential function, logarithmic function, and polynomial function) were applied to calculate LAI (Table 2).To examine the effect of the three techniques (TT, WF, and FT) on spatial heterogeneity (Experiment 2-4), we kept the linearity of each model approximately the same by uniformly distributing classes within each coarse pixel randomly from 0 to 1.After the elimination of model non-linearity effect, the scaling effect could be considered to be mainly caused by spatial heterogeneity represented by standard deviation of NDVI (σ vi ) within coarse pixels.According to Figure 1, the NDVI m at the coarse resolution is considered to be a simple average of the NDVIs within the coarse pixel at the high resolution (NDVI i ), where N is the number of high resolution pixels within the corresponding coarse pixel.Then, the application of f at the coarse resolution leads to LAI app (Path_2 in Figure 1): LAI exa is calculated by first applying f at the high resolution and then by the aggregating at the coarse resolution (Path_1 in Figure 1): The difference between LAI app and LAI exa is the scaling bias: Based on the Taylor's theorem, the bias could be approximated by a second order Taylor development of f around VI m : following Equation ( 5), where σ vi is the standard deviation of NDVI i within the coarse resolution pixel.Therefore, the scaling bias of LAI could be corrected based on the variance of the original image and the second order differential equation of the estimation model.
where LAI cor is the corrected value of LAI on the large scale.

The Wavelet-Fractal Technique (WF)
In the process of aggregation, some details will inevitably be lost.According to the wavelet transform theory, the information at the high resolution contains the information at the coarse resolution and some extra data.It is demonstrated that the additional information is the main reason for the scaling bias, which could be extracted by wavelet decomposition [31].As a quantitative representation of high-frequency information, wavelet coefficients from the detail components are suitable for analyzing the changes of remotely sensed data on different scales.Two-dimensional discrete wavelet transform was used to decompose a remotely sensed image on the small scale.Therefore, the high-frequency coefficients (high) can be defined as follows: where cD h j+1 , cD v j+1 , cD d j+1 and cD1 are the component in the horizontal, vertical and diagonal direction, respectively.With self-similarity and fractional dimensionality, remotely sensed images are proven to have fractal characteristics.Based on the advantages of the two methods, the wavelet-fractal technique was first developed by [38] to analyze and correct the spatial scaling bias of LAI.It is found that the relationship between high-frequency coefficient and scaling bias can be built in log-log coordinate, and the slope of the fit line is the fractal dimension.bias = a × high b (10) where a and b are the parameters of error correction equation.Therefore, the corrected LAI on the large scale could be obtained according to:

The Fractal Theory Technique (FT)
The concept of fractal dimension (D) is a statistical index of complexity comparing how detail in a pattern changes with the scale at which it is measured and it rests in unconventional views of scaling and dimension.Based on the fractal theory, the measurements perform power function correlation with the scales [42]. ) where LAI m,j is the LAI value of the jth (j = 1, 2, ..., (S/n) 2 ) coarse pixel on the scale of m × 2 m (m = 1, 2, ..., n) and S × S is the size of the original remotely sensed image with the high spatial resolution.VI i,g is the ith (i = 1, 2, ..., m 2 ) high spatial resolution NDVI value within the gth (g = 1, 2, ..., (n/m) 2 ) pixels with a size of m × m original pixels.d n,j is the slope of the linear fitting curve of the scale m versus the LAI value LAI m,j in log-log coordinate, which relates to the information fractal dimension D of the jth coarse pixel [43].When scale m equals 1, LAI 1,j is the exact LAI value of the jth coarse resolution pixel.
In addition, LAI 1,j also equals K in Equation (13).Hence, the transformation of Equations (2), (3) and ( 13) are as follows: Because both the information fractal dimension (D) and the standard deviation (σ vi ) can reflect the heterogeneity characteristic of coarse mixed pixel, a correlation should exist between them [41].It is found that D-2 has power law dependence at σ vi , so a linear empirically determined function is used to estimate D: where a and b are the fitting parameters of D-2 and σ vi in log-log coordinate.Note that the Fractal theory (FT) used in this study refers in particular to the pixel-based method by [41].

Performance Assessment
To evaluate the performance of the three techniques (TT, WF, FT) for correcting the scaling bias, two indicators root-mean-square error (RMSE) and relative error (RE) were defined as follows: where LAI est is the estimated LAI (LAI app or LAI cor ).In this study, RMSE was used to characterize the average error of the whole image and RE was used to quantify the error for a single pixel.Higher values of RMSE and RE indicated lower estimation of LAI or the poorer performance of techniques on the correction of scaling bias.The numerical simulation experiments and the evaluation of the three techniques were implemented using custom-written MATLAB (2016a) scripts.

Model Non-Linearity
The characteristics of model non-linearity and the relationship of model non-linearity with scaling bias are shown in Figure 5.The non-linearity of the power function relating LAI to NDVI (LAI = a NDVI ) became more significant when the base number a increased from 1 to 10 (Figure 5a).Therefore, we used the base number a to represent the non-linearity of the power model.It is obvious that the scaling bias increases with the base number a (Figure 5b).That means the more significant non-linearity of LAI estimation model would lead to greater scaling effect.Because both the information fractal dimension (D) and the standard deviation ( ) can reflect the heterogeneity characteristic of coarse mixed pixel, a correlation should exist between them [41].It is found that D-2 has power law dependence at , so a linear empirically determined function is used to estimate D: where a and b are the fitting parameters of D-2 and in log-log coordinate.Note that the Fractal theory (FT) used in this study refers in particular to the pixel-based method by [41].

Performance Assessment
To evaluate the performance of the three techniques (TT, WF, FT) for correcting the scaling bias, two indicators root-mean-square error (RMSE) and relative error (RE) were defined as follows: where LAIest is the estimated LAI (LAIapp or LAIcor).In this study, RMSE was used to characterize the average error of the whole image and RE was used to quantify the error for a single pixel.Higher values of RMSE and RE indicated lower estimation of LAI or the poorer performance of techniques on the correction of scaling bias.The numerical simulation experiments and the evaluation of the three techniques were implemented using custom-written MATLAB (2016a) scripts.

Model Non-Linearity
The characteristics of model non-linearity and the relationship of model non-linearity with scaling bias are shown in Figure 5.The non-linearity of the power function relating LAI to NDVI (LAI = a NDVI ) became more significant when the base number a increased from 1 to 10 (Figure 5a).Therefore, we used the base number a to represent the non-linearity of the power model.It is obvious that the scaling bias increases with the base number a (Figure 5b).That means the more significant nonlinearity of LAI estimation model would lead to greater scaling effect.To correct the scaling bias caused by model non-linearity, the three techniques (TT, WF, and FT) were applied to random datasets with the same standard deviation.Figure 6 shows that TT and WF had better performance of correcting the scaling bias than FT.When the model became more nonlinear (represented by the base number a), the scaling bias became much larger after the correction by FT.To correct the scaling bias caused by model non-linearity, the three techniques (TT, WF, and FT) were applied to random datasets with the same standard deviation.Figure 6 shows that TT and WF had better performance of correcting the scaling bias than FT.When the model became more nonlinear (represented by the base number a), the scaling bias became much larger after the correction by FT.

Spatial Heterogeneity
The Class-Specific Proportion within A Coarse Pixel Assuming there were only two classes (NDVI1 and NDVI2) within each coarse pixel, the relationships between scaling bias and proportion (p1) of the class NDVI1 within a coarse pixel are shown in Figure 7.By changing p1 from 0% to 100%, scaling bias firstly increased and then decreased for all types of estimation models (power, exponential, logarithmic, and polynomial functions).Around 50% of p1, the scaling bias reached the peak.However, when using different models, the class-specific proportions corresponding to the maximum scaling bias were not the same.

Spatial Heterogeneity
The Class-Specific Proportion within A Coarse Pixel Assuming there were only two classes (NDVI 1 and NDVI 2 ) within each coarse pixel, the relationships between scaling bias and proportion (p 1 ) of the class NDVI 1 within a coarse pixel are shown in Figure 7.By changing p 1 from 0% to 100%, scaling bias firstly increased and then decreased for all types of estimation models (power, exponential, logarithmic, and polynomial functions).Around 50% of p 1 , the scaling bias reached the peak.However, when using different models, the class-specific proportions corresponding to the maximum scaling bias were not the same.To correct the scaling bias caused by model non-linearity, the three techniques (TT, WF, and FT) were applied to random datasets with the same standard deviation.Figure 6 shows that TT and WF had better performance of correcting the scaling bias than FT.When the model became more nonlinear (represented by the base number a), the scaling bias became much larger after the correction by FT.

Spatial Heterogeneity
The Class-Specific Proportion within A Coarse Pixel Assuming there were only two classes (NDVI1 and NDVI2) within each coarse pixel, the relationships between scaling bias and proportion (p1) of the class NDVI1 within a coarse pixel are shown in Figure 7.By changing p1 from 0% to 100%, scaling bias firstly increased and then decreased for all types of estimation models (power, exponential, logarithmic, and polynomial functions).Around 50% of p1, the scaling bias reached the peak.However, when using different models, the class-specific proportions corresponding to the maximum scaling bias were not the same.

The Between-Class Spectral Difference within a Coarse Pixel
The scaling bias became larger when the difference between classes (NDVI 1 and NDVI 2 ) within each coarse pixel increased (Figure 8).The scaling bias for the exponential function was the most sensitive to the between-class spectral difference and that for polynomial function was relatively stable.Regardless of function types, the greater difference of classes within a coarse pixel would introduce higher scaling bias and lead to more significant scaling effect.The Between-Class Spectral Difference within a Coarse Pixel The scaling bias became larger when the difference between classes (NDVI1 and NDVI2) within each coarse pixel increased (Figure 8).The scaling bias for the exponential function was the most sensitive to the between-class spectral difference and that for polynomial function was relatively stable.Regardless of function types, the greater difference of classes within a coarse pixel would introduce higher scaling bias and lead to more significant scaling effect.

The Number of Classes within a Coarse Pixel
The scaling bias varied by the number of classes within a coarse pixel (Table 3).When a coarse pixel contained two classes, the greatest difference between NDVI1 (0.01) and NDVI3 (0.9) resulted in the largest scaling bias, which was consistent with the finding on the experiment of "the betweenclass spectral difference".When the initial classes were NDVI1&NDVI2 or NDVI2&NDVI3, the scaling bias became larger after the third class (NDVI3 or NDVI1) was added.When NDVI2 as the third class was added to the coarse pixel which initially contained NDVI1 and NDVI3, the scaling bias decreased.The scaling bias would increase when the value of newly added class (NDVInew) was more than or less than the two existing classes (NDVInew < NDVI2 < NDVI3 or NDVI1 < NDVI2 < NDVInew).It would decrease when the value of the third class was between the initial classes (NDVI1 ≤ NDVInew < NDVI3).

Correction of Scaling Bias Caused by Spatial Heterogeneity
The scaling bias values for different estimation models were all reduced by using the three techniques (Figure 9).As for WF, when the spatial heterogeneity ( ) became larger, the relative error (RE) firstly decreased to 0 and then rebounded to the peak.The valley values of RE were distributed approximately in the range of standard deviation from 0.27 to 0.29, especially for the power and polynomial models.The greatest difference in correction performance occurred to TT, which strongly depended on the LAI estimation model.TT almost eliminated the values of RE for the polynomial model, but had lowest effect for the power model among the three techniques.The performance of FT was more stable than that of others for the spatial heterogeneity of different levels.

The Number of Classes within a Coarse Pixel
The scaling bias varied by the number of classes within a coarse pixel (Table 3).When a coarse pixel contained two classes, the greatest difference between NDVI 1 (0.01) and NDVI 3 (0.9) resulted in the largest scaling bias, which was consistent with the finding on the experiment of "the between-class spectral difference".When the initial classes were NDVI 1 &NDVI 2 or NDVI 2 &NDVI 3 , the scaling bias became larger after the third class (NDVI 3 or NDVI 1 ) was added.When NDVI 2 as the third class was added to the coarse pixel which initially contained NDVI 1 and NDVI 3 , the scaling bias decreased.The scaling bias would increase when the value of newly added class (NDVI new ) was more than or less than the two existing classes (NDVI new < NDVI 2 < NDVI 3 or NDVI 1 < NDVI 2 < NDVI new ).It would decrease when the value of the third class was between the initial classes (NDVI 1 ≤ NDVI new < NDVI 3 ).

Correction of Scaling Bias Caused by Spatial Heterogeneity
The scaling bias values for different estimation models were all reduced by using the three techniques (Figure 9).As for WF, when the spatial heterogeneity (σ vi ) became larger, the relative error (RE) firstly decreased to 0 and then rebounded to the peak.The valley values of RE were distributed approximately in the range of standard deviation from 0.27 to 0.29, especially for the power and polynomial models.The greatest difference in correction performance occurred to TT, which strongly depended on the LAI estimation model.TT almost eliminated the values of RE for the polynomial model, but had lowest effect for the power model among the three techniques.The performance of FT was more stable than that of others for the spatial heterogeneity of different levels.

Image Level
Based on the numerical simulation experiments, the influence of model non-linearity and spatial heterogeneity were analyzed separately with different LAI estimation models.In applications to real data, both factors needed to be considered.The RMSE values after the correction of scaling bias were illustrated on different scales (Figure 10).At 4-m, 8-m, and 16-m resolution, the RMSE values from the three techniques for different models were similarly approaching to zero.After 16 m, the RMSE of WF increased to the peak at 64 m, and then decreased as the spatial scale became larger.This trend was the same for the different LAI estimation models and was similar to the results of correcting scaling bias mainly caused by spatial heterogeneity on the numerical simulation experiment (Figure 9).In contrast to WF, the performances of TT and FT were highly dependent on the type of LAI estimation model.For the exponential model, the RMSE was significantly reduced (<0.0052m 2 /m 2 ) after correction by FT.For the polynomial model, the LAI scaling bias was almost completely eliminated by TT at all resolutions.

Image Level
Based on the numerical simulation experiments, the influence of model non-linearity and spatial heterogeneity were analyzed separately with different LAI estimation models.In applications to real data, both factors needed to be considered.The RMSE values after the correction of scaling bias were illustrated on different scales (Figure 10).At 4-m, 8-m, and 16-m resolution, the RMSE values from the three techniques for different models were similarly approaching to zero.After 16 m, the RMSE of WF increased to the peak at 64 m, and then decreased as the spatial scale became larger.This trend was the same for the different LAI estimation models and was similar to the results of correcting scaling bias mainly caused by spatial heterogeneity on the numerical simulation experiment (Figure 9).In contrast to WF, the performances of TT and FT were highly dependent on the type of LAI estimation model.For the exponential model, the RMSE was significantly reduced (<0.0052m 2 /m 2 ) after correction by FT.For the polynomial model, the LAI scaling bias was almost completely eliminated by TT at all resolutions.

Pixel Level
For clarity, we selected the Landsat-like upscaled image with 32-m spatial resolution to compare the performance of the three techniques (TT, WF, FT) at pixel level (Figure 11).The results of other scales were similar.Except for logarithmic model (Figure 11c), the approximated LAI values (LAI app ) of the other three models were all distributed above the 1:1 line before correction.This means the scaling effect would lead to an underestimation (for power, exponential, and logarithmic models) or an overestimation (for logarithmic model) of LAI.In general, the values of RMSE were reduced by an average of 90% after the correction of LAI scaling bias and the corrected LAI (LAI cor ) were closer to the 1:1 line compared with LAI app .Particularly, the correction effect of LAI estimated the polynomial model (in Figure 11(d-1)) was remarkable with RMSE of 0.2 × 10 −6 m 2 /m 2 based on TT.

Pixel Level
For clarity, we selected the Landsat-like upscaled image with 32-m spatial resolution to compare the performance of the three techniques (TT, WF, FT) at pixel level (Figure 11).The results of other scales were similar.Except for logarithmic model (Figure 11c), the approximated LAI values (LAIapp) of the other three models were all distributed above the 1:1 line before correction.This means the scaling effect would lead to an underestimation (for power, exponential, and logarithmic models) or an overestimation (for logarithmic model) of LAI.In general, the values of RMSE were reduced by an average of 90% after the correction of LAI scaling bias and the corrected LAI (LAIcor) were closer to the 1:1 line compared with LAIapp.Particularly, the correction effect of LAI estimated from the polynomial model (in Figure 11(d-1)) was remarkable with RMSE of 0.2 × 10 −6 m 2 /m 2 based on TT.

Subpixel Level
In Figure 12, the spatial heterogeneity of study area on the large scale of 32 m was represented by the standard deviation of NDVI (σ vi ).There were 16 × 16 sub pixels of 2-m resolution within each coarse pixel.In the extreme case (the minimum and the maximum) of spatial heterogeneity, the factors (the class-specific proportion, the between-class spectral difference, and the number of classes within coarse pixel) were extracted (Figure 12b-g).For the pixel with minimum σ vi , it was relatively pure with 100% of winter wheat (Figure 12d) and the range of NDVI was only 0.017 (Figure 12b).For the four LAI estimation models, the RE values from the three techniques were all less than 0.00012 m 2 /m 2 .However, the performance for the maximum heterogeneity was more complex.In the pixel of maximum σ vi , the values of NDVI ranged from 0.251 to 0.753 (Figure 12c).Although the winter wheat occupied the largest proportion (64.5%), there were two other classes (Figure 12e) within the coarse pixel.The built-up resulted in the low values of NDVI (Figure 12c).If we used the uncorrected value of LAI on large scale to replace the exact LAI of winter wheat in this pixel, it would cause a RE up to 0.28 m 2 /m 2 .After correction, the RE values were reduced to less than 0.15 m 2 /m 2 .Particularly, despite the different heterogeneity, the lowest RE (~0 m 2 /m 2 ) occurred in the polynomial model by using TT, which was the same as the previous results.Moreover, it is worth noting that more classes within one pixel did not correspond to greater heterogeneity.The scaling effect produced by the three classes (winter wheat, other vegetation, and built-up) was even greater than that caused by complex multiple classes.

Subpixel Level
In Figure 12, the spatial heterogeneity of study area on the large scale of 32 m was represented by the standard deviation of NDVI ( ).There were 16 × 16 sub pixels of 2-m resolution within each coarse pixel.In the extreme case (the minimum and the maximum) of spatial heterogeneity, the factors (the class-specific proportion, the between-class spectral difference, and the number of classes within coarse pixel) were extracted (Figure 12b-g).For the pixel with minimum , it was relatively pure with 100% of winter wheat (Figure 12d) and the range of NDVI was only 0.017 (Figure 12b).For the four LAI estimation models, the RE values from the three techniques were all less than 0.00012 m 2 /m 2 .However, the performance for the maximum heterogeneity was more complex.In the pixel of maximum , the values of NDVI ranged from 0.251 to 0.753 (Figure 12c).Although the winter wheat occupied the largest proportion (64.5%), there were two other classes (Figure 12e) within the coarse pixel.The built-up resulted in the low values of NDVI (Figure 12c).If we used the uncorrected value of LAI on large scale to replace the exact LAI of winter wheat in this pixel, it would cause a RE up to 0.28 m 2 /m 2 .After correction, the RE values were reduced to less than 0.15 m 2 /m 2 .Particularly, despite the different heterogeneity, the lowest RE (~0 m 2 /m 2 ) occurred in the polynomial model by using TT, which was the same as the previous results.Moreover, it is worth noting that more classes within one pixel did not correspond to greater heterogeneity.The scaling effect produced by the three classes (winter wheat, other vegetation, and built-up) was even greater than that caused by complex multiple classes.and (e) maximum σ vi .By using TT (Taylor's theorem), WF (Wavelet-Fractal technique), and FT (Fractal theory), the correction results of pixels with: (f) minimum σ vi ; and (g) maximum σ vi for the four LAI estimation models (power, exponential, logarithmic model and polynomial models) were obtained.

The Influences of Model Non-Linearity and Spatial Heterogeneity on Scaling Effect
The main causes of scaling effect are the model non-linearity and the spatial heterogeneity.Before correcting the scaling bias, their influences on the scaling effect should be investigated.This study used numerical experiments to simulate different conditions.
The model non-linearity has not been quantitatively analyzed in previous studies.Because most studies focused on correcting the scaling bias on real data, the influence of model non-linearity was difficult to be extracted.Our results proved the common view that the scaling bias became larger when the non-linearity of LAI estimation model increased.Because the NDVIs in both Path_1 and Path_2 (Figure 1) were calculated at the fine resolution, the scaling bias caused by the non-linearity of NDVI as a function of red and near infrared reflectances was not considered.However, another way to upscale LAI is firstly aggregating the reflectances and then calculating NDVI and LAI at the coarse resolution.In this case, the scaling bias would be caused not only by the nonlinear relationship between LAI and NDVI, but also by the non-linearity of NDVI equation relating to the reflectances.Therefore, how the non-linearity of the two equations concerning NDVI-reflectance and LAI-NDVI relationships influences the scaling effect needs to be investigated in the future to understand the overall non-linearity of the LAI-reflectance model.
For the spatial heterogeneity, different heterogeneous structures within a coarse pixel would result in different scaling effect.The scaling bias firstly increased and then decreased when the proportion of one class within a coarse pixel became larger.Usually the class that occupies a larger area is called "the dominant class" and the low-proportion class is considered as "the heterogeneous class".In the numerical experiment, we simulated two contrasting classes (NDVI 1 and NDVI 2 ) within a coarse pixel.The proportion of the dominant class (NDVI 2 ) was relatively reduced when the proportion of the heterogeneous class (NDVI 1 ) increased from 0% to ~50%.The scaling bias became larger during this process.When the proportion of NDVI 1 continued to increase to 100%, it became the dominant class instead of NDVI 2 and the scaling bias dropped to zero (Figure 7).This is consistent with the finding in [36].He selected nine areas with different proportions of water from 0% to 93% in Canada from Landsat TM and found that the scaling bias was smaller than 2% for a pure forest pixel and would be up to 45% for a pixel containing around 50% open water [36].Therefore, we suggested that the scaling problem for pure pixels be ignored for LAI estimations.The larger difference of classes would lead to higher heterogeneity and more significant scaling effect.This is easy to understand and is consistent with various studies [15][16][17][18]36].The scaling bias would not monotonically increase or decrease with the number of classes within a coarse pixel.In the numerical simulation experiment (Table 3), the scaling bias increased when the newly added class with the highest or the lowest value.On the contrary, the scaling bias would decrease when the third class was added as "transitional class" with the value between the initial two classes.It is possible because a "compensation effect" weakens the heterogeneity.For the real data (Figure 12), the greatest scaling effect appeared in the pixel with three classes rather than the pixel with the maximum number of classes, as found by previous studies [11,44].
None of these factors existed independently in the process of scale transformation and they collectively influenced on the scaling effect.The non-linearity is related not only to the model form but also to the between-class spectral difference.For example, there are four subpixels with values of 0.11, 0.12, 0.13 and 0.14 in a coarse pixel and the other contains four subpixels with values of 0.1, 0.3, 0.5 and 0.7.We applied a simple non-linearity function with form of y = x 2 to the two coarse pixels.It is obvious that the model non-linearity for the two pixels is markedly different.In addition, internal composition of subpixels within coarse pixel can also further influence the model non-linearity such as the proportions and the number of classes.That means the model non-linearity is affected by the spatial heterogeneity.This implies that the influence of spatial heterogeneity on scaling effect is greater than the non-linearity of LAI estimation model [28,29,45].

The Reliability and Practicability of the Three Techniques
The correction results of TT strongly depended on the model type and it was perfect for the polynomial model.Perhaps because the Taylor's theorem was proposed by decomposing a function, the function form directly affected correction results.The performance of WF exhibited lower stability.WF could reduce LAI scaling bias ideally on relatively small scale but not for large scale compared with TT and FT.The results of simulations (Figure 9) showed that WF could achieve great correction results in a certain range of heterogeneity (standard deviation of VI, σ vi ).It could be explained that the standard deviation of remotely sensed VI image with 64-m resolution was out of the optimal range for WF.This is probably because WF is based on the adjacent scales.When the gap between the original scale and the coarse scale becomes bigger, the bias of each scale on the gap would be accumulated.Therefore, WF is more suitable for neighboring scales.It is found that TT and WF had a better performance of correcting the scaling bias caused by model non-linearity than FT.This is because FT was proposed by considering the spatial heterogeneity but not the non-linearity of estimation model [41].According to the results in Figure 6, this property of FT was prominently proven.However, FT had good performance of correcting the scaling bias for the real data.
The advantages and disadvantages of the three techniques (TT, WF, FT) were summarized in Table 4.With considering both model non-linearity and spatial heterogeneity, TT and WF can be used to correct the scaling bias of LAI.By extending the application of TT, WF has been successfully used for the model without a clear function expression such as the radiation transfer model [38].Therefore, WF is of high universality for different models.However, WF is not suggested for continuous scales possibly due to the dyadic property of the Haar wavelet.Generally, TT is a better selection than WF if the second-order stationarity hypothesis holds.Although FT showed a potential capability of correcting the scaling bias, it was not recommended to correct the scaling bias caused by great model non-linearity.It is necessary to assess its performance for physical models in future work.Moreover, the uncertainties of the three techniques should be analyzed in different landscape patterns.Additionally, it is necessary to further investigate the application of the optimal technique selected in this study to promote the accuracy and integrity of the precision agriculture, such as multi-scale monitoring and prediction, data fusion, and data assimilation.

Conclusions
The scaling bias became greater when the non-linearity of LAI estimation model increased.For the spatial heterogeneity, different heterogeneous structures within a coarse pixel would lead to different scaling effect.None of these factors existed independently and they collectively influenced on the scaling effect.Compared with spatial heterogeneity, the model non-linearity was not the key factor of scaling effect and could be ignored in some applications especially for pure pixel with only one class.
The Taylor's theorem (TT) can be used when the second-order stationarity hypothesis is satisfied.With low RE and RMSE (~0 m 2 /m 2 ) values for both simulated and real data, TT was suggested for correcting the scaling bias with a polynomial function relating LAI to NDVI.The Wavelet-Fractal technique (WF) was preferable for neighboring scales rather than continuous scales.The Fractal theory technique (FT) only focused on the description of spatial heterogeneity, and it was not recommended for correcting the scaling bias caused by significant model non-linearity.

Figure 1 .
Figure 1.Schematic flowchart for the up-scaling process: f is a function relating y to x; xi (i = 1, 2, 3, 4) is the value of each pixel at the original high resolution and X is their average.Suppose x1 = x2 = VIa, x3 = x4 = VIb, then X is VIm representing the vegetation index value of the coarse spatial resolution pixel,

Figure 1 .
Figure 1.Schematic flowchart for the up-scaling process: f is a function relating y to x; x i (i = 1, 2, 3, 4) is the value of each pixel at the original high resolution and X is their average.Suppose x 1 = x 2 = VI a , x 3 = x 4 = VI b , then X is VI m representing the vegetation index value of the coarse spatial resolution pixel, and y i is the corresponding LAI value.The subscript exa means the exact value and app represents the approximated value of the coarse resolution pixel.

Figure 2 .
Figure 2. Location of the study area in Huai'an, Jiangsu Province (b), China (a) and 2-m resolution GF-1 satellite image of the study area (c).The band combination of this false color composite image is NIR (0.77 μm to 0.89 μm), R (0.63 μm to 0.69 μm), and G (0.52 μm to 0.59 μm).Shown on the right of the satellite image (c) is the layout of subsamples in each sample.

Figure 2 .
Figure 2. Location of the study area in Huai'an, Jiangsu Province (b), China (a) and 2-m resolution GF-1 satellite image of the study area (c).The band combination of this false color composite image is NIR (0.77 µm to 0.89 µm), R (0.63 µm to 0.69 µm), and G (0.52 µm to 0.59 µm).Shown on the right of the satellite image (c) is the layout of subsamples in each sample.

Figure 3 .
Figure 3. Relationships between field-measured NDVI and LAI as represented by four fitting functions.

Figure 3 .
Figure 3. Relationships between field-measured NDVI and LAI as represented by four fitting functions.

Figure 3 .
Figure 3. Relationships between field-measured NDVI and LAI as represented by four fitting functions.

Figure 4 .
Figure 4. Mapping of leaf area index (LAI) from the GF-1 satellite image: (a) the land use map of the study area; (b) the NDVI image derived from the GF-1 image; and (c) the LAI map generated from the power function relating LAI to NDVI (Table2).

Figure 5 .
Figure 5.The relationship of model non-linearity with scaling bias.(a) The non-linearity of LAI estimation model (LAI = a NDVI ) represented by the base number a.(b) The relationship between base number a and the scaling bias.The grey lines represent the simulations of 10,000 groups and the black line represent the average scaling bias.

Figure 5 .
Figure 5.The relationship of model non-linearity with scaling bias.(a) The non-linearity of LAI estimation model (LAI = a NDVI ) represented by the base number a.(b) The relationship between base number a and the scaling bias.The grey lines represent the simulations of 10,000 groups and the black line represent the average scaling bias.

Figure 6 .
Figure 6.Correction of the scaling bias caused by model non-linearity using the three techniques TT, WF, and FT.The non-linearity of LAI estimation model (LAI = a NDVI ) was represented by the base number a.

Figure 7 .
Figure 7.The relationship between scaling bias and the proportion of class (NDVI1) within coarse pixels for four types of LAI estimation models.

Figure 6 .
Figure 6.Correction of the scaling bias caused by model non-linearity using the three techniques TT, WF, and FT.The non-linearity of LAI estimation model (LAI = a NDVI ) was represented by the base number a.

Figure 6 .
Figure 6.Correction of the scaling bias caused by model non-linearity using the three techniques TT, WF, and FT.The non-linearity of LAI estimation model (LAI = a NDVI ) was represented by the base number a.

Figure 7 .
Figure 7.The relationship between scaling bias and the proportion of class (NDVI1) within coarse pixels for four types of LAI estimation models.

Figure 7 .
Figure 7.The relationship between scaling bias and the proportion of class (NDVI 1 ) within coarse pixels for four types of LAI estimation models.

Figure 8 .
Figure 8.The relationship between scaling bias and differences of classes (NDVI1 and NDVI2) within a coarse pixel.

Figure 8 .
Figure 8.The relationship between scaling bias and differences of classes (NDVI 1 and NDVI 2 ) within a coarse pixel.

Figure 10 .
Figure 10.Correction of scaling bias with the four LAI estimation models at different spatial resolutions using the three techniques.The LAI estimation models are: (a) power model; (b) exponential model; (c) logarithmic model; and (d) polynomial model.The three techniques used to correct scaling bias are: TT (Taylor's theorem), WF (Wavelet-Fractal technique), and FT (Fractal theory).The spatial scales are: 4 m, 8 m, 16 m, 32 m, 64 m, and 128 m.

Figure 10 .
Figure 10.Correction of scaling bias with the four LAI estimation models at different spatial resolutions using the three techniques.The LAI estimation models are: (a) power model; (b) exponential model; (c) logarithmic model; and (d) polynomial model.The three techniques used to correct scaling bias are: TT (Taylor's theorem), WF (Wavelet-Fractal technique), and FT (Fractal theory).The spatial scales are: 4 m, 8 m, 16 m, 32 m, 64 m, and 128 m.

Figure 11 .
Figure 11.Correction of scaling bias with the four LAI estimation models at 32-m spatial resolution using the three techniques.The LAI estimation models are: (a) power model; (b) exponential model; (c) logarithmic model; and (d) polynomial model.The three techniques used to correct scaling bias are: (1) Taylor's theorem (TT); (2) Wavelet-Fractal technique (WF); and (3) Fractal theory (FT).LAIapp is the approximated LAI before correction, and LAIcor is the corrected LAI.The unit for RMSE is m 2 /m 2 .

Figure 11 .
Figure 11.Correction of scaling bias with the four LAI estimation models at 32-m spatial resolution using the three techniques.The LAI estimation models are: (a) power model; (b) exponential model; (c) logarithmic model; and (d) polynomial model.The three techniques used to correct bias are: (d-1) Taylor's theorem (TT); (d-2) Wavelet-Fractal technique (WF); and (d-3) Fractal theory (FT).LAI app is the approximated LAI before correction, and LAI cor is the corrected LAI.The unit for RMSE is m 2 /m 2 .

Figure 12 .
Figure 12.Comparison of LAI scaling bias correction in the extreme case of spatial heterogeneity.(a) The spatial heterogeneity of study area at 32-m resolution was obtained by the standard deviation of NDVI ( ).Within each pixel, there are 16 × 16 subpixels with spatial resolution of 2 m.NDVI distribution in the pixel of: (b) minimum ; and (c) maximum .The classes in the pixel of: (d) minimum; and (e) maximum .By using TT (Taylor's theorem), WF (Wavelet-Fractal technique), and FT (Fractal theory), the correction results of pixels with: (f) minimum ; and (g) maximum for the four LAI estimation models (power, exponential, logarithmic model and polynomial models) were obtained.

Figure 12 .
Figure 12.Comparison of LAI scaling bias correction in the extreme case of spatial heterogeneity.(a) The spatial heterogeneity of study area at 32-m resolution was obtained by the standard deviation of NDVI (σ vi ).Within each pixel, there are 16 × 16 subpixels with spatial resolution of 2 m.NDVI distribution in the pixel of: (b) minimum σ vi ; and (c) maximum σ vi .The classes in the pixel of: (d) minimum σ vi ;and (e) maximum σ vi .By using TT (Taylor's theorem), WF (Wavelet-Fractal technique), and FT (Fractal theory), the correction results of pixels with: (f) minimum σ vi ; and (g) maximum σ vi for the four LAI estimation models (power, exponential, logarithmic model and polynomial models) were obtained.

Table 1 .
Statistical estimation models of LAI.

Table 1 .
Statistical estimation models of LAI.

Table 1 .
Statistical estimation models of LAI.
• Experiment 3, the between-class spectral difference.A coarse pixel contained two classes (NDVI 1 , NDVI 2 ) with the same proportions (p 1 , p 2 ) of 50%.The spectral values of the two classes were different, but the sum of spectral values remained the same (NDVI 1 + NDVI 2 = 1).It meant the value of the coarse pixel (average values of fine subpixels) was unchanged.

Table 2 .
Design of numerical simulation experiments.

Table 4 .
Summary of characteristics of the three techniques.