The research area is Shangri-La City ( Figure 1 ), Diqing Tibetan Autonomous Prefecture, Yunnan Province, China (Latitude: 26°52′~28°52′N, Longitude: 99°20′~100°19′E). The area has significant changes in altitude and is a key forest area, ecological protection area, and tourist area. The dominant tree species include Spruce (Picea asperata), Fir (Abies fabri), Oak (Quercus semecarpifolia), Pinus (Pinus yunnanensis), etc (Xu et al., 2021).

In this study, four ML algorithms were used to build the estimation model of NFCC based on light spot footprints, and then the NFCC of all light spot footprints was predicted and used as a spatial attribute of spatial heterogeneity analysis. Our framework approach comprises three main components ( Figure 2 ): (1) the process of preparing and preprocessing data, including data preprocessing for ATL08, resampling of terrain factors and extraction of slope and aspect; (2) NFCC model construction and evaluation based on four ML models (k-NN, SVM, RF, GBRT), and (3) NFCC spatial effects based on semi-variogram function and terrain impact analysis based on Pearson correlation.

In statistics, the Pearson correlation coefficient measures the linear correlation between two variables, with its value ranging from -1 to 1. The closer the coefficient’s absolute value is to 0, the weaker the linear correlation between the two variables. Conversely, an absolute value closer to 1 indicates a stronger linear correlation. The basic principle of Pearson correlation can be seen in the article of Yang et al (Yang et al., 2021). In this study, the correlation analysis was used to screen model independent variables and explain the effect of terrain factors on NFCC.

k-nearest neighbor (k-NN), a simple and efficient non-parametric method, effectively circumvents the issue of collinearity among independent variables. This algorithm is applicable to the parameter estimation of remote sensing data characterized by non-normal distributions and unknown density functions, and is extensively utilized in forestry investigations globally (Chirici et al., 2008, 2016). The fundamental concept of this algorithm involves using a point in the feature space as the reference object, capturing the attribute values of the k nearest sample points relative to this point, and determining the predicted value of this object by calculating the average of its inverse distance weights.

Support vector machine (SVM) algorithm originates from the VC dimension theory and the structural risk minimization principle (Chirici et al., 2008, 2016). The fundamental principle of SVM involves mapping training data features to a high-dimensional space through a defined kernel function, and identifying an optimal linear regression hyperplane in this space that best fits the eigenvalues.

RF proposed by Breiman (Breiman, 2001), is a method that combines weak classifiers to create strong classifiers. Its fundamental concept lies in the ability of the ensemble to compensate for incorrect predictions made by individual weak classifiers. Originally developed as an extension of classification and regression trees (CART), RF enhances predictive models by generating aggregate predictors (Breiman, 1996).

Gradient Boost Regression Tree (GBRT), as an ensemble learning method, builds a strong learning model by sequentially aggregating a set of weak CART regression tree submodels (Opitz and Maclin, 1999; Friedman, 2001). The key concept of GBRT is that each new regression tree submodel is built in the gradient direction of residuals reduction to reduce residuals from previous models (Liu et al., 2020).

In our research, we randomly divided the plot data into two sets: training dataset (70%) and validation dataset (30%). The training set served to train and develop the models, whereas the validation set, not participating in the model-building process, was used to evaluate model performance. Root mean square error [RMSE; Equation (2)] and coefficient of determination [(R²; Equation (3)], as two commonly used evaluation indexes, were used to evaluate the prediction performance of regression models. A higher R² value indicates greater model accuracy, while a lower RMSE value signifies enhanced accuracy of the regression model. The calculation formulas for each indicator are as follows:

where n is the number of samples, y ^ i is the predicted by the ML models, y i is the observed FCC from the ground measurements, y ¯ is the arithmetic mean of observed values.

Moran’s I can effectively capture differences and correlations in the spatial distribution of observations, as well as reflecting the overall clustering pattern of objects in the study area (Zhang et al., 2023). The value interval of Moran’s I is [-1, 1]. When the value is less than 0, spatial objects have a negative correlation; when the value is equal or close to 0, spatial autocorrelation does not exist. When the value is greater than 0, it indicates that there is spatial autocorrelation, and spatial objects show a clustered distribution (Moran, 1950). Moran’s I formula [Equation (4)] is as follows:

where n is the number of observed values; X ¯ is the average of the variable X ; x p and x q refer to the observation values at plot p and plot q , respectively; w p q ( d ) is the spatial weight matrix value between plot p and plot q .

After calculating the Moran’s I value, the significance of Moran’s I can be tested by a Z-score. A positive Z-score points to a cluster of high values, whereas a negative Z-score suggests clusters of low values. The degree of clustering is greater (or lesser) with a higher Z-score. Conversely, if the Z-score is close to zero, it indicates the absence of significant clustering in the area. Equation (5) was used to calculate the Z-score.

where Z ( S ) represents an index that measures the intensity of a spatial agglomeration pattern; E ( S ) denotes the expected value of the index value I, while V a r ( S ) represents its variance.

Furthermore, local spatial autocorrelation can elucidate the level of spatial correlation existing between a research object and its neighboring units. Equation (6) was used to calculate the local Moran"s I.

where m is the number of plots; x p and x q are the observation values at plot p and plot q , x ¯ is the average of all NFCC values; w p q   ( d ) is the weight matrix value. The local Moran’s I differs from the global Moran’s I in terms of value range. Unlike the global Moran’s I, the local Moran’s I is not limited to the range of (-1, 1). If the I l value is positive, it indicates a positive correlation in the location of the plot and reflects the aggregation of similar values. Conversely, if the I l is negative, it signifies a negative correlation in the location of the plot and reflects the aggregation of different values.

Semi-variogram is often used to describe spatial heterogeneity (Wang et al., 2000). In the semi-variogram function parameter, nugget (C₀) reflects the possible degree of randomness within the regionalized variable, and explains the discontinuous variation of the regional variable at a small scale, which is due to the measurement error and random variation of the sampling scale. Sill (C₀+C) was used to measure the degree of spatial heterogeneity and reflected the maximum degree of variation of the variable. Range (a) refers to the average maximum distance of spatial autocorrelation between variables (Chiles and Delfiner, 2012). The semi-variogram formula [Equation (7)] (Chiles and Delfiner, 2012) is as follows:

where r ( h ) is the semi-variogram of NFCC; N ( h ) is the total logarithm of sample points spaced h in one direction; Z ( x i ) is the measured NFCC at spatial position x i ; Z ( x i + h ) is the NFCC value at h distance from point x i .

The relationship between the semi-variogram value and the distance usually requires a mathematical model to fit. The common mathematical models include spherical model, exponential model, Gaussian model, and linear model. According to the principle that the R² is large and the RSS is tiny, it is found that the exponential model is more suitable for revealing the spatial heterogeneity of the NFCC. The expression of the exponential model [Equation (8)] is as follows:

where γ ( h ) is the semi-variogram of NFCC; a is the range; C is the partial sill value; C₀ is nugget value; h is distance.

Spatial heterogeneity is not only related to scale, but also to direction (Li and Reynolds, 1995), and the spatial distribution of NFCC is also different depending on the direction. Anisotropic semi-variogram was used to analyze the direct change of spatial heterogeneity of NFCC (Habin et al., 1998). In general, the anisotropy ratio [ K ( h ) ] between the semi-variogram functions in different directions is used to describe the anisotropic structural characteristics, and the formula [Equation (9)] is as follows:

where λ ( h , θ 1 ) and λ ( h ,   θ 2 ) are the semi-variogram functions in the directions θ 1 and θ 2 , respectively. If K ( h ) is equal to or close to 1, the spatial heterogeneity is isotropic, otherwise it is anisotropic.

Ordinary Kriging (OK), its essence is to infer the regionalized spatial distribution of variables by the variable in the spatial regionalization of a finite number of sample points. Based on the information of several measured points in the search field where the point to be estimated is the center of the circle, it uses the semi-variogram function as a tool to calculate the weighted value of the measured points around the point to be estimated, and finally makes the optimal and unbiased estimation of the estimated points (Christakos, 2000). The formula [Equation (10)] is as follows:

where Z e # ( x 0 ) represents the predictive NFCC at the point to be predicted; Z ( x i ) stands for the observed NFCC at the point to be predicted; λ i represents the weight of each known parameter value, and m represents the number of spot footprints.

Sequential Gaussian Conditional Simulation (SGCS) is a spatial stochastic simulation method that constructs a Gaussian function based on known data and treats each value of the regionalized random variable Z(x) as a random realization of the Gaussian function (i.e. the normal distribution function) F(x). It is mainly used to generate spatial explicit estimates of interest variables based on regionalized variable theory and spatial autocorrelation theory measured by the semi-variogram functions. More information and detailed processes about the SGCS can be found in the articles by Luo et al (Luo et al., 2023), and Zhao et al (Zhao et al., 2010).

Zhao et al (Zhao et al., 2010). derived the connection between the OK and the SGCS, and compared the statistical parameters of both computational results and the original data. The results showed that the values from SGCS actually consist of two parts: one is the result of Kriging interpolation, and the other is a random deviation with a mathematical expectation of 0 and a variance equal to the Kriging error variance S (X_m). The difference lies in that Kriging interpolation solely uses the original known point data as a basis for estimating unknown points, whereas in SGCS, each simulated value for a location not only applies known point data but also previous simulation data. In this section, two interpolation methods, Kriging and SGCS, are applied. The motive is to find a suitable interpolation method for the footprint NFCC obtained by inversion.

In the evaluation of interpolation results, we randomly divided the footprint NFCC predicted by the optimal model into two sets: interpolation dataset (80%) and validation dataset (20%). The interpolation dataset is used for spatial interpolation, and the validation dataset is used to evaluate the interpolation results.

With the help of SPSS 27.0, Pearson correlation analysis was used to evaluate the correlation between ATL08 parameters and NFCC, and the ATL08 parameters were selected as independent variables of the model according to the value and significance of the correlation coefficients.

Based on Python 3.7 environment, the four machine learning algorithms (k-NN, SVM, RF, GBRT) in this study were implemented using the scikit-learn package in the Python library.

Global and local Moran’s I were calculated in the spatial analysis toolbox of ArcGIS 10.8.

GS⁺ 9.0 (GeoStatistics for the Environmental Sciences, version 9.0), a comprehensive geostatistics program, provides all geostatistics components, from variogram analysis through Kriging and mapping, in a single integrated program that is widely praised for its flexibility and user-friendly interface. The semi-variogram analysis and the corresponding parameters (C₀, C₀+C, a) were obtained in this software. To further know the spatial distribution of NFCC within the study area, based on the NFCC of ATLAS footprints and the fitted semi-variogram function, GS⁺ 9.0 and ArcGIS 10.8 were used to achieve the OK interpolation and SGCS for NFCC ( Figures 5A, B ). The number of SGCS simulations was set to 50 times (Luo et al., 2023).

The result of Pearson correlation analysis showed seven parameters from the ATL08 product (asr, landsat_perc, photon_rate_can, toc_roughness, n_toc_photons, h_canopy, h_dif_canopy) were significantly correlated with NFCC at the 0.05 confidence level ( Figure 6 ). Then the seven parameters were selected as the independent variables (the description is shown in Table 2 ), and the NFCC measurement serves as the dependent variable for constructing the ML models.

Two statistical metrics (R², RMSE) were applied to evaluate the models constructed, utilizing the reserved 30% of field plot data ( Table 3 ). The predicted performance of the ML models was ranked as follows (descending order): RF (R² = 0.75, RMSE = 0.09), GBRT (R² = 0.60, RMSE = 0.12), SVM (R² = 0.45, RMSE = 0.14), k-NN (R² = 0.43, RMSE = 0.15). Figure 7 showed the comparison of the NFCC predicted values with the measured values in test set.

Due to its greater predicted performance, the RF model was used to predict NFCC within the natural forest footprint, and then the footprint NFCC was visualized in Figure 8 . Most of the natural forest footprint CC was above 0.5. The areas with high-values were mainly distributed in the northwest, middle and south of Shangri-La ( Figure 8 ).

Table 4 showed the footprint’s descriptive statistics of NFCC and topographic factors within the footprint. The P-P Plot of NFCC ( Figure 9 ) showed a normal distribution, which meets the requirements of structural analysis of semi-variogram.

Table 5 showed that Z-score = 6.47 and P value< 0.01, indicating that Moran’s I passed the test with 99% confidence level. Moran’s I of NFCC in the study area is positive (Moran’s I = 0.36), indicating that the NFCC has a positive spatial correlation and belongs to spatial agglomeration distribution.

Local spatial autocorrelation analysis can capture local spatial elements’ clustering and difference characteristics. As a common index of local spatial autocorrelation, the local Moran index was used to continue exploring the NFCCs’ spatial relationships of each footprint. As shown in Figure 10 , NFCC located in the central and northern parts of the study area showed significant HH clustering, while the spatial clustering pattern of the NFCC situated on the west and east sides of the study area showed HL outliers. LH outliers were mainly concentrated in the middle of the study area. Besides, LL clusters were primarily concentrated in the eastern part of the study area.

The fitting of the semi-variogram function in this study was implemented in GS⁺ 9.0 software. The semi-variogram function revealed the regional differences in the spatial structure of NFCC, and the fitting results were shown in Figure 11 . According to the principle that the Residual Sum of Squares (RSS) is minimum and the coefficient of determination (R²) is maximum, the exponential model (R² = 0.61, RSS = 1.96×10⁶) is best fit to describe the relationship between values and distances. The abutment value (C₀ + C) of the exponential model is 0.89×10^-2, the partial sill value (C) is 0.77×10^-2, the variable range (A₀) is 10200 m, the nugget value (C₀) is 0.12×10^-2, and. The NSR of NFCC is 13.40%, indicating that the variables have strong spatial autocorrelation within the range. The expression of the exponential model [Equation (11)] is as follows:

The results of the anisotropic semi-variogram function showed that the NFCC changes in all directions on the scale of 106 km ( Figure 12 ). The anisotropy of NFCC in the northwest-southeast direction (θ = 135°) was the most obvious, followed by the north-south direction (90°). However, the anisotropy of NFCC in the east-west (0°) direction was relatively low.

To realize the influence degree of topography on NFCC, Pearson correlation analysis was conducted between the NFCC and the topographic factors ( Table 6 ). Table 6 showed that the NFCC in the study area is significantly correlated with elevation, slope, and aspect at 0.01 level. The order according to the correlation coefficient’s absolute value absolute value of the correlation coefficient was as follows: elevation > slope > aspect.

The interpolation results showed that the spatial distribution map of NFCC obtained by the OK method ( Figure 5A ) was roughly similar to the SGCS interpolation map ( Figure 5B ). High values were concentrated mainly in the study area’s northern, central, and southern regions. As can be seen from Figure 5A , the spatial distribution of NFCC obtained by the OK method is relatively continuous and has obvious smoothing effect. In Figure 5B , the overall distribution using the SGCS method is relatively discrete and less affected by smoothing effect. In addition, the spatial predicted values obtained by SGCS were in good agreement with NFCC footprint values ( Figure 5D , R² = 0.59). In contrast, the spatial predicted values obtained by OK were less consistent with the NFCC footprint values ( Figure 5C , R² = 0.43).

The results of this study confirm that the combination of ATLAS data, plot survey, ML algorithm, and geostatistical method can provide a valuable framework for county-scale NFCC spatial effects analysis. Before the spatial effects analysis, the NFCC of 1106 footprints was estimate by the RF model. As the input of the model, the measured sample plots play an important role in the modeling. In general, the more samples used, the more reliable the model. However, Shangri-la has many high-altitude mountains and complex terrain, which makes it not easy to collect samples based on LiDAR footprints. In addition, the existing remote sensing estimation of forest parameters is based on the traditional empirical sample size, that is, below 30 is a small sample size, and above 50 is a large sample size (Shu et al., 2022). To minimize the labor and time required, exploring the optimal sample data is needed. Shu et al (Shu et al., 2022). solved the optimal sample size by integrating the statistical variance function and value coefficient, which was reconstructed using the model accuracy evaluation index RMSE and the model sample cost. Therefore, the optimization of the sample size can be further performed in the future to minimize costs.

Although the number of sample plots is small, the prediction accuracy of the estimation model based on the measured samples was great. However, the phenomenon of high underestimation in all models is relatively easy to find. From the scatter plot ( Figures 7A–D ), when the NFCC is above 0.7, the predicted value below the 1:1 line can be visually seen, which means that the model is still underestimating at higher NFCC. Previous studies (Xing et al., 2010; Xi et al., 2022) have shown that incorporating distinct forest types into the modeling process can enhance performance and decrease the model’s reliance on training samples. However, because of the absence of sample plot data, it was impossible to distinguish between forest types or NFCC levels for modeling. In order to reduce uncertainties in the modeling process, it is recommended that sufficient sample plot data be collected in the future. Furthermore, the importance of physical geography, bioclimate, and biology in estimating forest parameters has been demonstrated (Su et al., 2016; Fayad et al., 2021). Therefore, in future studies, it is suggested that remote sensing data should be combined with a forest physiological process model to enhance the generalization and accuracy of the predictive model.

The spatial effect analysis of NFCC obtained by combining machine learning algorithms, relatively new remote sensing data sources, and measured samples is one of the innovations of this study. Many spatial heterogeneity studies (Yao et al., 2015; Li et al., 2017; Liu et al., 2018) can only be carried out on a small scale due to labor costs, resulting in little difference in environmental factors. But the spatial heterogeneity of NFCC is often the result of the interaction of topography, climate, soil, stand, external disturbance, and other random factors. The distribution of light, temperature, water, and other climatic factors is determined by topographic differences. Analyzing the influences of topographic factors on the spatial heterogeneity of forest parameters can provide a better understanding of the mechanism of climate-forest interaction, which is often overlooked in current studies.

The results of this study showed that the canopy cover of Shangri-la natural forest is moderately variable, indicating that it is susceptible to structural and random factors. Spatial variation is mainly divided into structural factor variation and random factor variation (Zou et al., 2021). The NSR of NFCC is 0.13 in section 3.4, showing strong spatial autocorrelation, indicating that the influence of natural factors was dominant. Since 1999, the China National Forestry and Grassland Administration has implemented many large-scale forest conservation projects, such as the Natural Forest Protection Project and Grain for Green Project. Shangri-La is the key area for the implementation of these projects. The natural forests are less disturbed by human factors.

With the support of large-scale spaceborne LiDAR data, this study found that the elevation itself has the most significant influence on the spatial distribution of NFCC, followed by the slope and aspect. However, this study lacks the relationship between topography and climate factors and their joint influence on the spatial distribution of NFCC, which needs further analysis in the future.

As shown in Figure 8 , with the assistance of ATLAS, the ability to estimate large-scale NFCC is obtained, which is limited to the predetermined ground track footprint range. In view of the feature of ATLAS discontinuous sampling, the discontinuous spatial attributes (NFCC) were used to analyze spatial effects, and extend to the extent of natural forest land throughout the study area by the spatial statistical method. Because spatial interpolation uses known spatial attributes for prediction, the limitations (e.g., climate effect, saturation effect) associated with the using optical images are significantly eliminated (Liu et al., 2022; Yu et al., 2023).

However, the spatial interpolation based on LiDAR spot footprint still faces many problems, such as banding effect, and smoothing effect. Compared with OK, the SGCS method overcomes the shortcomings of Kriging’s smoothing effect (Luo et al., 2023). However, since the spot footprints are distributed along the track, the spatial interpolation results located around the track may have a strong banding effect. Increasing the randomness of spot footprint distribution often has a great effect on avoiding banding effect. In this study, spatial interpolation based on the 1106 footprints obtained through systematic sampling from 11,060 footprints located within natural forests did not show a significant banding effect. Therefore, the sampling of the footprint can be used as an alternative scheme to increase the randomness of the footprint. In addition, Liu et al (Liu et al., 2022). integrated ICESat-2 and GEDI data to carry out spatial interpolation of forest canopy height, and the interpolation results did not show obvious banding effect. Therefore, adding other spaceborne LiDAR data sources can also avoid banding effects.

In this study, the research object only focuses on the NFCC in Shangri-La. Still, the proposed method can be extended to other areas or forest parameters after the same treatment. The Earth will be observed further by ICESat-2/ATLAS, yielding additional high-precision orbital observation data. In order to obtain more and denser space observation footprints, another spaceborne LiDAR named GEDI (Global Ecosystem Dynamics Investigation) can be introduced in future research.