The InVEST carbon storage model was applied to calculate the carbon stocks of rubber plantations under actual planted and hypothetical non-rubber-planted conditions in this study. However, this model is based on a simplified carbon cycle, where the carbon storage is a static inventory, assuming that each hectare of land is identical and constant [52]. This may bias the carbon storage estimates for rubber and natural forests, whose carbon storage gradually accumulates and increases with time [80]. To this end, based on obtaining the rubber forest stand age and initial land-cover types before rubber planting and referring to the IPCC calculation methods and result criteria for carbon storage of agricultural and forestry land [67], each carbon storage component of rubber and natural forests in different years over time was given a new definition and assignment, and then was added to the InVEST carbon storage model to improve the accuracy of carbon storage simulation to some extent. This theoretical approach has not been reported in previous related InVEST carbon storage simulation studies [9,11,51].

According to the carbon storage results under the actual rubber forest plantation scenario, we found that rubber forest plantation can increase regional carbon density when occupying cropland, but through deforestation and clearing it would cause a rapid decrease in carbon density, and then gradually increase and recover to the storage stock level before deforestation. This result was consistent with the findings of some existing studies [9,69,81]. By comparing the hypothetical scenario of no rubber planting, we found that rubber planting reduced the regional carbon storage and caused regional carbon emissions. Moreover, with the expansion of rubber forests, the difference between the actual carbon storage and the hypothetical scenario became larger and larger, causing a large carbon storage gap of −0.15 million tons in 2018. This result was also consistent with the findings of some studies on carbon emissions caused by rubber plantations in tropical regions [7,80,82,83]. In order to achieve sustainable rubber forest plantation and management, the carbon storage gap caused by deforestation and clearing must be bridged by partial replacement of cropland with rubber in the future to achieve the carbon balance of rubber forest plantation and environmental sustainability.

#### *4.4. Limitations*

In this study, the impact of rubber forest expansion on regional carbon storage was investigated through rubber forest extraction, forest age estimation, pre-planting land cover type identification, carbon storage estimation, and different scenario simulations. However, it should be noted that eight spatial texture features and four structural features of the optical and SAR images were combined in the rubber forest mapping. How to further combine the diverse textural and structural features of both to improve the accuracy of fine rubber forest mapping will be the focus of later research.

Second, the tree age estimation of a rubber forest relies only on the remote sensing NDVI to construct time-series curves. During the planting and growing process of rubber forests, the land-cover state undergoes a transformation from natural forest/cropland to bare land, rubber with a low canopy cover, and rubber with a high canopy cover. Different land-cover states also have obvious differences in texture characteristics, and the application of texture characteristic time series to the tree age estimation of rubber forests will be explored in the future.

Finally, the application of the InVEST model for carbon storage estimation is more sensitive to the input data. Researchers generally collect carbon density data from the literature or field experimental observations in their study area [9,11,52]; because of the lack of observation data, only the relevant research algorithms and results of the IPCC [66,67] were used as input for carbon storage simulation in this study. This may cause uncertainty in the simulation results, and it poses difficulties in the validation of the model results owing to the lack of measured data [11,19,50]. Although the results simulated by the carbon storage module of the InVEST model have been shown to be reasonably accurate [51] and representative [9,11], validation of the measured data can increase the effect of the carbon density simulation to some extent [52]. In the future, when conditions permit, partial field observations will be carried out to compare the field data with IPCC data for testing whether the data consistent, to then further correct the results of the carbon densityrelated components and improve the simulation accuracy. On the other hand, the carbon accumulation process of cropland was not considered in the carbon storage simulation. Later, the law and process of cropland soil carbon sequestration in tropical regions will be further explored to improve the simulation accuracy and eventually improve the simulation accuracy of regional carbon storage caused by rubber expansion.

#### **5. Conclusions**

Based on the multi-source satellite time series data of Sentinel-1, Sentinel-2, and Landsat, coupled with random forest, shapelet, and InVEST carbon storage model, a systematic approach was proposed for estimating the carbon storage of regional rubber forests, and then explored the impact of rubber forest expansion on regional carbon storage. The conclusions of the study are as follows.

old). Before 2004, rubber was mainly grown through encroachment on cropland. After that, rubber conversion from natural forests started to increase.

**Supplementary Materials:** The following supporting information can be downloaded at: https://www. mdpi.com/article/10.3390/rs14246234/s1. Table S1. List of spectral indices generated for Sentinel-2 imagery; Table S2: Carbon density of each component in different land-cover types in different years (units: t/ha); Table S3: Confusion matrix of the forest and rubber plantation mapping; Table S4: Confusion matrix of identification results of rubber planting years.; Figure S1. Importance statistics of spectral features. References [84–94] are citied in the Supplementat Materials.

**Author Contributions:** Conceptualization, C.H. and H.L.; methodology, C.Z. and H.L.; data curation, C.Z. and H.L.; writing—original draft preparation, H.L.; writing—review and editing, C.H.; supervision, C.H.; project administration, C.H.; funding acquisition, C.H. All authors have read and agreed to the published version of the manuscript.

**Funding:** This study was supported by the National Natural Science Foundation of China (42130508) and the CAS Earth Big Data Science Project (XDA19060302).

**Acknowledgments:** We acknowledge all who contributed to the data collection and processing, as well as the constructive and insightful comments by the editor and anonymous reviewers.

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

#### **References**

## *Article* **Modelling the Dynamics of Carbon Storages for** *Pinus densata* **Using Landsat Images in Shangri-La Considering Topographic Factors**

**Yi Liao 1, Jialong Zhang 1,\*, Rui Bao 2, Dongfan Xu 3 and Dongyang Han 4**

**Abstract:** Accurate estimation of forest carbon storage is essential for understanding the dynamics of forest resources and optimizing decisions for forest resource management. In order to explore the changes in the carbon storage of *Pinus densata* in Shangri-La and the influence of topography on carbon storage, two dynamic models were developed based on the National Forest Inventory (NFI) and Landsat TM/OLI images with a 5-year interval change and annual average change. The three modelling methods used were partial least squares (PLSR), random forest (RF) and gradient boosting regression tree (GBRT). Various spectral and texture features of the images were calculated and filtered before modelling. The terrain niche index (TNI), which is able to reflect the combined effect of elevation and slope, was added to the dynamic model, the optimal model was selected to estimate the carbon storage, and the topographic conditions in areas of change in carbon storage were analyzed. The results showed that: (1) The dynamic model based on 5-year interval change data performs better than the dynamic model with annual average change data, and the RF model has a higher accuracy compared to the PLSR and GBRT models. (2) The addition of TNI improved the accuracy, in which R2 is improved by up to 10.48% at most, RMSE is reduced by up to 7.32% at most, and MAE is reduced by up to 8.89% at most, and the RF model based on the 5-year interval change data has the highest accuracy after adding TNI, with an R2 of 0.87, an RMSE of 3.82 t-C·ha−1, and a MAE of 1.78 t-C·ha−1. (3) The direct estimation results of the dynamic model showed that the carbon storage of *Pinus densata* in Shangri-La decreased in 1987–1992 and 1997–2002, and increased in 1992–1997, 2002–2007, 2007–2012, and 2012–2017. (4) The trend of increasing or decreasing carbon storage in each period is not exactly the same on the TNI gradient, according to the dominant distribution, as topographic conditions with lower elevations or gentler slopes are favorable for the accumulation of carbon storage, while the decreasing area of carbon storage is more randomly distributed topographically. This study develops a dynamic estimation model of carbon storage considering topographic factors, which provides a solution for the accurate estimation of forest carbon storage in regions with a complex topography.

**Keywords:** Landsat; *Pinus densata*; terrain niche index; dynamic model; carbon storage

#### **1. Introduction**

Forests are the mainstay of terrestrial ecosystems, which store 60% of the carbon in the terrestrial ecosystem [1]. Forests play an indispensable role in balancing and regulating CO2 in the atmosphere with their powerful carbon sink function. Global warming is a great threat to humanity today, and the excessive emissions of greenhouse gases, of which CO2 is an important member, is one of the main reasons for global warming. Atmospheric CO2

**Citation:** Liao, Y.; Zhang, J.; Bao, R.; Xu, D.; Han, D. Modelling the Dynamics of Carbon Storages for *Pinus densata* Using Landsat Images in Shangri-La Considering Topographic Factors. *Remote Sens.* **2022**, *14*, 6244. https://doi.org/ 10.3390/rs14246244

Academic Editor: Wataru Takeuchi

Received: 12 November 2022 Accepted: 6 December 2022 Published: 9 December 2022

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

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

concentrations have risen dramatically over the past 100 years (they have risen by approximately 110 ppm) [2], a phenomenon that has accelerated the rate of global warming. The accurate estimation of forest carbon storage is related to the ability to reduce atmospheric CO2 concentration. Against this background, the necessity of the study of forest carbon storage is revealed and its accurate estimation is an important basis for regional and even global carbon change studies.

The methods for estimating forest carbon storage includes both traditional direct measurement methods and indirect estimation methods based on remote sensing [3,4]. The traditional method is based on sample site survey and data statistics, which require a large investment of human and financial resources [5], and due to cost and time constraints, there are significant limitations regarding the large range of study objectives and long study periods. In contrast, remote sensing-based methods have the characteristics of being fast, low-cost, large-scale, and less destructive [1,6]. With the rapid development of quantitative remote sensing technology, the results of carbon storage estimation based on remote sensing have become more and more stable and reliable. Therefore, the use of remote sensing methods for the quantitative estimation of forest carbon storage on a larger scale has become a major current trend [7,8].

At present, there is a significant amount of research in remote sensing estimation, and while scholars in various fields have constructed different remote sensing models [9], most of these are static models. The static model is a model developed directly from the attribute values of the sample sites and the remotely sensed feature values of the corresponding points, which are uncalculated state data. If these values are further processed and calculated to produce different types of change values (periodic change, annual average change, rate of change, etc.), then the model developed using these change data is the dynamic model. Change data can describe the different processes of change in forest ecosystems at a given time and have the ability to monitor forest dynamics [10]. Gómez et al. [11] and Zhang et al. [12] compared static models with dynamic models and showed that dynamic models have higher accuracy and better predictive power than static models, with the highest R2 for dynamic models in their study being 0.90 and 0.94, respectively. In the study of forest dynamics, the dynamic model is able to directly estimate the corresponding level of forest change and directly respond to forest dynamics, which improves the efficiency of the long-term forest dynamics study. The calculation of change data needs to be supported by multi-period data. The current National Forest Inventory (NFI) projects in many countries provide the possibility for long-term forest dynamics modelling studies. The ongoing NFI projects in many countries and open-access Landsat time series data offer the potential for long-term forest dynamics modelling studies.

Trees accumulate carbon mainly through photosynthesis during growth, and topographical factors influence the photosynthetic effect by affecting light and water conditions in the environment, which eventually affects the carbon storage of trees [13]. Therefore, topographic factors are an important concern in forest carbon storage studies. The digital elevation model (DEM) provides valuable topographic information including elevation, slope, and aspect [14]. These topographic factors are often used as independent factors in various studies because they do not require complex calculations and processing, are relatively easy to obtain, and reflect the effects of topography on different aspects of vegetation to some extent [15]. However, the influence of topography on forest carbon storage is the result of the combined effect of various topographic factors, and it is difficult to reflect this combined effect by a single elevation or slope. The terrain niche index (TNI) combines elevation and slope information, which is able to reflect the spatial variability of regional elevation and slope [16,17]. Therefore, in past studies, it has been frequently applied to land use change, landscape patterns and ecological effects [18,19]. The topography of forest ecosystems is complex and diverse, so the comprehensive influence of topographic factors should also be considered when modelling carbon storage. The change data can reflect the degree of change in forest carbon storage, and TNI can reflect the comprehensive influence of topographic factors on forest carbon storage, so the combination of the two

is important for the accurate estimation of change in forest carbon storage in areas with complex topography.

Shangri-La, Yunnan Province, is rich in forest resources, and *Pinus densata* is one of the major tree species, which is widely distributed and has an essential impact on the carbon balance of the region. This region is located at the edge of the Tibetan Plateau, with an overall high altitude, and the terrain comprises high mountains and deep valleys with large topographic undulations, making it more difficult to estimate the forest carbon storage in this region accurately. Therefore, in this study, we used the remote sensing dynamic model combined with TNI to estimate the carbon storage of *Pinus densata* in Shangri-La. The main objectives are the following: (1) developing the carbon storage dynamics model based on two different kinds of change data; (2) exploring the impact of TNI on carbon storage models; and (3) estimating the historical carbon storage changes of *Pinus densata* in Shangri-La and analyzing the topographic effects of carbon storage changes.

#### **2. Materials and Methods**

#### *2.1. Study Area and Study Process*

The study area, Shangri-La, is located in Yunnan Province in southwestern China (Figure 1), at the border of the Yunnan Province and the Sichuan Province, and the geographical coordinates range from 99◦20 to 100◦19 eastern longitude and 26◦52 to 28◦52 northern latitude, the total area is 1.16 × 106 ha. The area is located in the southeastern part of the Tibetan Plateau (the average elevation is 3459 m), with the Hengduan Mountains running north and south through the whole territory, the regional topography is undulating [20], and the whole terrain is high in the northwest and low in the southeast. Shangri-La has a remarkable monsoon climate, with rainfall concentrated in the months of June to October each year [21], the average annual precipitation is 268~945 mm, and the annual sunshine is 1742.9~2186.6 h. This region is very rich in forest resources, with 76% of forest cover [22]. *Pinus densata*, *Pinus yunnanensis,* and *Picea asperata* are the main tree species in the region [23]. The target species for this study is the *Pinus densata*, which includes both artificial and natural forests distributed within the study area.

**Figure 1.** Study area and sample plots.

The main content of the study is shown in Figure 2. The research process began with the pre-processing of all the data; then, two types of change data (annual average change and 5-year interval change) were calculated, three methods were used to construct the model, and TNI was added to the model. Finally, in order to compare the effect of the dynamic model with different change data and to compare the effect of the model before and after the addition of TNI, different indicators were used to evaluate the accuracy.

**Figure 2.** Flow chart of the main content in this study.

#### *2.2. Ground Survey Data and Carbon Storage Calculation*

This study focused on *Pinus densata*, and the ground survey data were obtained from the National Forestry Inventory, which contains a total of 136 *Pinus densata* sample plots (comprising pure forests of *Pinus densata* or forests with *Pinus densata* as the main species). The years of the collection were 1987, 1992, 1997, 2002, 2007, 2012 and 2017, where the number of sample plots was 19, 22, 23, 16, 16, 17 and 23 for each year, respectively (Figure 3). The dataset records information on tree height, diameter at breast height (DBH), number of trees, coordinate location, and major tree species. The distance between the sample plots consisted of regular distributions of 6 km × 8 km, and each sample plot was a rectangle of 28.28 m × 28.28 m (0.08 ha).

**Figure 3.** Sample plot data records by year: different colors represent sample plots with different locations (different numbered sample plots).

The aboveground biomass (AGB) of the sample plots was calculated based on the allometric growth equation of *Pinus densata* [24]. The average AGB was first calculated using the average tree height and average DBH of the sample plots. Then, the total AGB of the sample plots was calculated based on the average AGB and the number of *Pinus densata*. The allometric growth equation is as follows:

$$\text{AGB}\_1 = 0.073 \times \text{DBH}^{1.79} \times \text{H}^{0.880} \tag{1}$$

where AGB1 is single wood aboveground biomass (t), DBH is the diameter at breast height (cm), and H is tree height (m). We filtered the sample plot data. In this process, five sample plots with an AGB (AGB<1t·ha−1) that was too small were removed, followed by six outliers that were screened out and removed using the Pauta criterion [25]. According to the Pauta criterion, if a value is outside the range of three times the standard deviation of the mean (outside the range of x ± 3σ, where x is the mean and σ is the standard deviation), it is considered an outlier. This is a common method for outlier screening [26]. In the end, 125 sample plots of *Pinus densata* remained.

The carbon storage of the sample sites was calculated using the AGB multiplied by the carbon content coefficient. According to the guidelines for measuring carbon storage in forest ecosystems issued by the State Forestry and Grassland Administration of China [27], the average carbon content in the dry matter of *Pinus densata* is 0.501. The equation is as follows:

$$\text{CS} = 0.501 \times \text{AGB}\_2 \tag{2}$$

where CS is the carbon storage (t-C·ha−1), and AGB2 is the aboveground biomass of the sample plots (t·ha−1). The basic sample data is shown in Figure 4. The box plots (Figure 4) depict information on the diameter at breast height, tree height, and carbon storage for each year of the sample plots. The overall range of DBH is 6.0–91.7 cm, with a maximum annual mean of 32.0 cm and a maximum annual standard deviation of 27.8 cm; the overall range of tree height is 1.5–24.8 m, with a maximum annual mean of 11.4 m and a maximum annual standard deviation of 5.5 m; and the overall range of carbon storage is 0.7–84.2 t-C·ha−1, with a maximum annual mean of 38.6 t-C·ha−1 and a maximum annual standard deviation of 41.9 t-C·ha−1.

**Figure 4.** Boxplots of basic information on sample plots by year.

#### *2.3. Remote Sensing Images and Obtained Features*

#### 2.3.1. Remote Sensing Data

The remote sensing images used in this study were obtained from the USGS website (http://glovis.usgs.gov/ (accessed on 28 October 2021)). Landsat 5 TM images from 1987, 1992, 1997, 2002, 2007, and 2012 and Landsat 8 OLI images from 2017 of the Shangri-La region were downloaded from the website (Table 1), comprising 21 views in total and a spatial resolution of 30 m. When an image contains too much cloud or is of poor quality, we chose an image from a neighboring timepoint to replace it. Most of the selected images contained less than 5% cloud, and the cloud did not cover the study area when the cloud content was greater than 5%.

**Table 1.** Landsat time-series of imagery used in the study.

In order to improve the image quality, all the images were preprocessed: firstly, radiometric calibration was performed, which converts the original DN (digital number) value of the image to a consistent radiometric brightness in order to eliminate the influence of the sensor [28]; secondly, the Fast Line-of-sight Atmospheric Analysis of Spectral Hypercubes (FLAASH) method was used to eliminate the influence of the atmosphere [29]; thirdly, geometric correction was performed on each image with reference to the existing standard SPOT-5 images to eliminate geometric errors, where the error in geometric correction was

less than one pixel; and finally, in order to eliminate the change of image grey value caused by the topographic undulation, topographic correction was performed using the slope matching model [30]. After topographic correction, the brightness of the illuminated side was suppressed and the brightness of the shadow side was enhanced, which better restored the reflectance of the hidden features in the shadow of the original image [31].

#### 2.3.2. Remote Sensing Features

The remote sensing pixel values were evaluated by considering bands with the same name and close wavelengths in Landsat 8 and Landsat 5. In this study, two types of remote sensing features, spectral features and texture features, are calculated. Spectral features are the most basic and direct features of remote sensing images, but the phenomena of "the same objects with different spectra" and "different objects with the same spectra" commonly exist in remote sensing images, so it may be difficult to obtain good research results when only relying on spectral features [32,33]. The texture features can describe the grey scale information of an image, express the spatial distribution of the grey scale of pixels in the image, and reflect the structural information of the image. Compared with spectral features, texture features are less affected by the environment and more stable in information representation [34]. Therefore, the combination of spectral and textural features can effectively reflect the characteristics of the ground objects and their changes [35,36]. A total of 35 spectral features and 540 texture features were obtained in this study (the calculations of various types of textures from all single bands of each image in odd windows from 1 to 19) (Table 2).

**Table 2.** Extracted remote sensing feature information.

Where B1 is the blue band, B2 is the green band, B3 is the red band, B4 is the near-infrared band, B5 is the shortwave infrared-1 band, and B7 is the shortwave infrared-2 band.

#### *2.4. Terrain Niche Index*

To explore the influence of elevation and slope on the carbon storage dynamics model, the terrain niche index (TNI), a composite factor combining elevation and slope, was established in this study [41]. Elevation and slope were extracted from the digital elevation

model (DEM). The DEM used in this study was an ASTER GDEM data product with a spatial resolution of 30 m. The DEM was downloaded from the Geospatial Data Cloud website (http://www.gscloud.cn/ (accessed on 18 August 2022)). The equation for TNI is as follows:

$$\text{TNI} = \log\left[ \left( \frac{\text{E}}{\text{E}} + 1 \right) \times \left( \frac{\text{S}}{\text{S}} + 1 \right) \right] \tag{3}$$

where E and S are the elevation and slope of any point, respectively, and E and S are the average elevation and average slope of the study area, respectively. A larger TNI value indicates a higher elevation and slope; a smaller TNI value indicates a lower elevation and gentler slope; and a medium TNI value indicates a high elevation gentle slope, or a low elevation steep slope, or a medium elevation and slope [42].

#### *2.5. Distribution Index*

The distribution index (DI) is often used to show the distribution of various land use types. In this study, the areas of the increasing and decreasing carbon storage of *Pinus densata* were considered as two different land classes. Then, the DI was used to explore the distribution of carbon storage changes in different TNI gradients. The DI is a standardized, dimensionless metric that eliminates the effects of area differences and allows the comparison of distribution characteristics between carbon storage changes of different area proportions [18,43]. The equation is as follows:

$$\mathbf{P} = (\mathbf{S}\_{\text{ie}} / \mathbf{S}\_{\text{i}}) / (\mathbf{S}\_{\text{e}} / \mathbf{S}) \tag{4}$$

where P is the DI, e is the gradient of TNI, Sie is the area of i changes on the e TNI gradient, Si is the area of the i change, Se is the area of the TNI for gradient e, and S is the total area of the study area. A larger P indicates a higher frequency of this type of change, and P > 1 indicates that the type of change under this TNI gradient belongs to the dominant distribution. Therefore, the analysis of the distribution of each carbon storage change type on different TNI levels is able to reveal the influence of different topographic conditions on carbon storage change.

#### *2.6. Modelling Process*

#### 2.6.1. Calculation of Change Data

In this study, two kinds of change data were used to develop the forest carbon storage model, one is the 5-year interval change, and the other is the annual average change. The change data were calculated based on a common sample plot for both years: we calculate the change value of a sample plot when data were recorded in both adjacent years, otherwise the change value for the sample plot during that specific time period is not calculated. The equations are as follows:

$$
\Delta \mathbf{I} = \mathbf{I}\_\mathbf{n} - \mathbf{I}\_\mathbf{m} \tag{5}
$$

$$
\Delta \mathbf{R} = \frac{\Delta \mathbf{I}}{\mathbf{n} - \mathbf{m}} \tag{6}
$$

where ΔI is the interval change value, n and m are two different years in which the change is calculated (n > m), In and Im are the data values for years n and m, respectively, and ΔR is the average annual change value. The above two equations were used to calculate the change data of carbon storage and the corresponding remote sensing features in the sample sites. The calculation of change values is supported by the NFI data records every 5 years from 1987–2017. The shortest time interval for which change data can be calculated is 5 years and the longest time interval is 30 years. In generally, the shorter the time interval, the stronger the continuity of the data and the more accurately the dynamic process of forest change can be expressed [10,44]. Zhang et al. [12] tested the change data for 5, 10, and 15 year intervals and found that the 5 year interval resulted in the smallest error. Lunetta et al. [45] compared the effect of 3, 5, and 7 year intervals on monitoring change and showed that the 3 year interval had the highest accuracy. Therefore, only the change data for the shortest interval (5 years) were calculated in this study: the 5-year interval change and the annual average change during the 5-year period.

#### 2.6.2. Remote Sensing Features Selection

The efficiency of the model can be affected by the number of features, because too many features reduce the speed of model fitting [46]. In addition, the original features contain some redundant information, which has a negative impact on the model accuracy [47]. Therefore, the feature screening was performed. In this study, 575 extracted remote sensing features were screened using the feature importance assessment method of the random forest [48]. In order to show the combined effect of spectral features and texture features in modelling, the spectral feature and texture feature were screened independently in this study. These two features were then combined and used for modelling.

Since most of the texture values of the extracted 1 × 1 windows were 0 or 1, they were not usable and thus removed, and the remaining 486 texture features were used for the subsequent screening. The 5-year interval change data and the annual average change data were calculated using Equations (5) and (6), respectively. Then, the feature importance assessment method of the random forest was used to filter the feature variables in the two different change types, and the top five features with contributions greater than 5% were selected from the spectral features and the texture features, respectively. Considering the problem of collinearity between features, the variance inflation factor (VIF) was used to check all the features selected in this study. When the VIF value is less than 3, it means that the collinearity between features is weak [49]. According to this principle, two features with VIF values greater than 3 in the 5-year interval change were removed, and the 5-year interval change and annual average change resulted in 8 and 10 features for the modelling study, respectively (Table 3).

**Table 3.** Results of remote sensing feature screening.

The expression of the texture features is "RXBYTT", RX is the size of the texture window, BY is a certain single band of the image, and TT is the abbreviation of a certain texture feature.

#### 2.6.3. Modelling Methods

In this study, carbon storage dynamics models were developed using change data, and the modelling methods included partial least squares regression (PLSR), gradient boosting regression tree (GBRT), and random forest (RF).

PLSR is a method used to study the correlation and quantitative relationship between the response variable *Y* and a set of explanatory variables, *X* = *x*1, *x*2, ··· , *xn*, and the set of independent variables *X* can also be used to predict *Y* [50]. The modelling principle of PLSR is the development of multivariate linear regression and principal component analysis. Compared with these two methods, PLSR has better stability and can solve the problem of collinearity between multiple explanatory variables *X* [51–53]. The implementation of the PLSR model in this study was based on the Minitab 20 software [54].

GBRT is an algorithm for iterative regression trees proposed by Friedman [55,56]. All the regression trees in this model are interconnected, and at its core, each tree is fitted based on the residuals and conclusions of the previous tree [57]. GBRT uses a forward-distributed algorithm to minimize the loss function by selecting the appropriate decision tree function from the current model and the fitted function. Based on the characteristics described above, the GBRT model has low computational complexity, is able to reduce errors, and has the ability to handle unevenly distributed data [58]. The implementation of the GBRT model in this study is based on the "Gradient Boosting Regressor" algorithm provided in the "Scikit-learn" package for the Python language.

RF is an integrated learning model, first proposed by Breiman, which consists of many aggregated randomly generated trees [59,60]. In the classification problem, RF outputs the type with the most votes, and in the regression problem, RF outputs the mean of all decision trees [61]. RF can reduce variance and effectively reduce overfitting by assembling different trees, and it has excellent classification and regression performance; thus, it is currently used in many fields of research [58,61,62]. The implementation of the RF model in this study is based on the "Random Forest Regressor" algorithm provided in the "Scikit-learn" package for the Python language.

In this study, the above three methods were used to develop 5-year interval change and annual average change models, respectively. In order to explore the influence of TNI on the dynamic model of forest carbon storage, the TNI was added to each model in this study, then the accuracy of different dynamic models before and after adding the TNI was compared, and the optimal model was selected to estimate the carbon storage of *Pinus densata* in Shangri-La.

#### *2.7. Accuracy Evaluation*

Ultimately, this study obtained 92 groups of change values of each type of change; 70% (64 groups) of the data were randomly selected for model fitting, and the remaining 30% (28 groups) were used for validation, and cross-validation was performed during model fitting. The indicators used to evaluate the accuracy of the model were the coefficient of determination (R2), the root mean square error (RMSE), and the mean absolute error (MAE). In order to ensure that the model results were as objective as possible, each model was fitted 20 times in this study to allow take the mean values of evaluation indicators to be used for comparison. The equations are as follows:

$$\mathbf{R}^2 = 1 - \frac{\sum\_{i=1}^n (\mathbf{y}\_i - \hat{\mathbf{y}}\_i)^2}{\sum\_{i=1}^n (\mathbf{y}\_i - \mathbf{y})^2} \tag{7}$$

$$\text{RMSE} = \sqrt{\frac{\sum\_{i=1}^{n} (\mathbf{y}\_i - \mathbf{y}\_i)^2}{n}} \tag{8}$$

$$\text{MAE} = \frac{1}{n} \sum\_{i=1}^{n} |\mathbf{y}\_i - \mathbf{\hat{y}}\_i| \tag{9}$$

where yi is the observed value of the sample site, yˆi is the predicted value of the model, y is the observed mean value of the sample site, and n is the number of samples.

#### **3. Results**

#### *3.1. Analysis of Modelling Results*

The PLSR, GBRT, and RF methods were used to develop the 5-year interval change and annual average change models of carbon storage, respectively, TNI was added to each model subsequently, and the results of accuracy evaluation indexes (20 times fitted mean values) of each type of model are shown in Table 4. As can be seen from Table 4, among the models developed with the same type of change data, the PLSR model performed the worst and the RF model performed the best. The change of model indicators (percentage) after adding TNI is shown in Figure 5, from which it can be seen that R2 increases and RMSE decreases for all models, and the MAE values decreased for all models except for the PLSR model with annual average change. The maximum increase in R2 is 10.48%, the maximum decrease in RMSE is 7.32%, and the maximum decrease in MAE is 8.89%. Considering the changes of all accuracy indicators together, the addition of TNI improves the model accuracy. The mean values of carbon storage calculated for the 5-year interval change and the mean annual change of the sample plot data were 5.14 t-C·ha−1 and 1.03 t-C·ha−1, respectively. After comparison, it can be seen that the RMSE and MAE values of both the RF and GBRT models are smaller than the mean values, except for the PLSR model, indicating that the errors of the RF and GBRT models are smaller and the model results are more reasonable.

**Table 4.** Comparison of modelling results.

PLSR TNI, GBRTTNI, and RFTNI in the table indicate the PLSR model, GBRT model, and RF model after adding TNI, respectively.

**Figure 5.** Changes in model indicators after adding TNI: positive values on the right side of the *y*-axis represent an increase, and negative values on the left side of the *y*-axis represent a decrease.

The 5-year interval change represents the change in carbon storage over five years, and the annual average change represents the change in carbon storage in one year. In order to compare the effects of dynamic models with 5-year interval changes and annual average changes, the RMSE and MAE values of the annual average change model were adjusted to five times their original values in this study, so that the error indicators of both models are on the same time scale. The adjusted results are shown in Table 5. Comparing the accuracy indicators of the two models in Tables 4 and 5, it can be seen that the R2 of all 5-year interval change models is higher than that of the annual average change model, the RMSE values are lower than those of the 5-year interval change model except for the PLSR model, and the difference in MAE values is not significant. Therefore, the 5-year interval change model works better overall.

According to the confidence interval plot (Figure 6) of the 20 fittings for each model, the confidence range of the RF model among the types is smaller than that of the GBRT and PLSR models, indicating that the RF model in this study has the least uncertainty and the results obtained are more stable and reliable. In summary, the accuracy of 5-year interval change RF model with the addition of TNI is the highest.

**Table 5.** Error indicator for the adjusted annual average change.

**Figure 6.** Confidence intervals (95%) of the accuracy evaluation indicators for different models: (**a**) the model of 5-year interval change before adding TNI; (**b**) the model of 5-year interval change after adding TNI; (**c**) the model of annual average change model before adding TNI; and (**d**) the model of annual average change after adding TNI.

#### *3.2. Mapping Carbon Storage*

According to the model evaluation results, the 5-year interval change RF model with the highest accuracy was used to estimate the carbon storage change. The scatter plot of this estimation model is shown in Figure 7. The direct estimation result of the model is the 5-year interval change of carbon storage of *Pinus densata* in Shangri-La. The values of the carbon storage change for the six change periods are ΔCS1 = CS1992 − CS1987, ΔCS2 = CS1997 −CS1992,ΔCS3 = CS2002 −CS1997, ΔCS4 = CS2007 −CS2002,ΔCS5 = CS2012 −CS2007 and ΔCS6 = CS2017 − CS2012.

**Figure 7.** Scatter plot of the estimation model.

The estimated carbon storage changes in the six periods are as follows: ΔCS1 = −53.327 × 104 t, ΔCS2 = 13.887 × 104 t, ΔCS3 = −94.041 × 104 t, ΔCS4 = 14.602 × 104 t, ΔCS5 = 20.869 × 104 t, and ΔCS6 = 5.602 × 104 t.  The estimation results indicate that the carbon storage of *Pinus densata* in Shangri-La decreased in the periods 1987–1992 and 1997–2002 and increased in the remaining four periods.

In order to obtain carbon storage maps for each year, the published AGB estimates for each year from our group were used [63], which contain the aboveground biomass estimates of *Pinus densata* in Shangri-La for 1987, 1992, 1997, 2002, 2007, 2012, and 2017. These data were combined with the results of carbon storage change values estimated in this paper to calculate the total carbon storage values for each year. Firstly, the AGB of each year in the published study was converted to carbon storage by multiplying the carbon content coefficient (0.501) and the converted values were expressed as T1987, T1992, T1997, T2002, T2007, T2012, and T2017; then the carbon storage change values ΔCS1, ΔCS2, ΔCS3, ΔCS4, ΔCS5 and ΔCS6 were added to the conversion values of the smaller years in each change period in order to obtain the carbon storage values for the six years 1992, 1997, 2002, 2007, 2012 and 2017, respectively (e.g., CS1992 = ΔCS1 + T1987), and the carbon storage value for 1987 was obtained by subtracting ΔCS1 from T1992 (CS1992 = T1987 − ΔCS1). The statistical results for carbon storage are shown in Table 6, and Figures 8 and 9 (the spatial resolution is 30 m) were obtained by mapping the values of the carbon storage change for each period and each year according to the *Pinus densata* distribution range. The distribution range and area of *Pinus densata* are derived from Forest Manager Inventory (FMI) data.

**Table 6.** Statistics of *Pinus densata* carbon storage.

**Figure 8.** Change values for carbon storage for six time periods.

**Figure 9.** Carbon storage values for seven years.

#### *3.3. Spatial Distribution Characteristics for the Changes in Carbon Storage on Different TNI Gradients*

In this study, the estimated changes in carbon storage were classified as increases (positive change values) and decreases (negative change values), and the TNI (0.30 to 0.93) was classified into ten levels from low to high using the natural breaks (Jenks) method: low (1–3), medium–low (4–5), medium–high (6–7), high (8–10). The DI was then calculated according to Equation (4) and used to explore the distribution of carbon storage changes along different TNI gradients. The distribution results are shown in Figure 10.

**Figure 10.** Distribution of increased carbon storage (**a**) and decreased carbon storage (**b**) on the TNI gradient.

As can be seen from Figure 10, the increase and decrease in the carbon storage of *Pinus densata* in Shangri-La at each period have a certain regularity on the TNI gradient. The DI curve for ΔCS1(1987–1992) decreases and then rises when carbon storage increases (Figure 10a) with a dominant distribution in the low and high TNI gradients; the DI curve for ΔCS2 (1992–1997), ΔCS4 (2002–2007), ΔCS5 (2007–2012), and ΔCS6 (2012–2017) rises and then reduces with the dominant distribution areas of ΔCS2 being in the low and medium TNI gradients, and the dominant distribution areas of ΔCS4, ΔCS5, and ΔCS6 being concentrated in the medium–high TNI gradients; and ΔCS3 (1997–2002) reduces with the TNI gradient, and its dominant distribution area is in the low and medium–low gradients. The corresponding curve change period displays a decreasing trend when the carbon storage decreases (Figure 10b). Overall, the dominant distribution areas for each period when carbon storage increases are mainly found in the lower or middle TNI gradients, while the dominant distribution areas for each period when carbon storage decreases are found in all TNI gradients.

#### **4. Discussion**

In long time series data, the uncertainty of the historical data, due to the events that took place at the time being taken into account, cause some of the data to display large deviations in values, which are known as anomalies or outliers [64]. The stability of the model can be affected by outliers, and at the same time, the prediction accuracy of the model may be reduced, so the detection and handling of outliers is very important in the modelling process [65]. The sample plot data were screened twice in this study, the first time removing plots with an AGB of less than 1 t·ha−1 because they had an average diameter at breast height of less than 5 cm and an average tree height of less than 1.5 m. They are young forests, too low in canopy density, and poorly characterized by the appearance of the forest floor (similar to bare ground) on 30 m resolution imagery. Values exceeding the range of *x* ± 3*σ* were selected for the second time using the Pauta criterion. By removing values that are too small and too large from the data in these two steps, respectively, the overall uncertainty of the data is reduced and the reliability of the data is improved.

The remote sensing-based estimation model is one of the main methods used in current carbon storage studies [7]. Based on the accumulation of ground survey data, the results of remote sensing estimation models are considered to have certain advantages and reliability in forest carbon storage studies. Time series images also provide an important basis for describing change [66]. The NFI and Landsat time series data over a 30-year period were combined to develop a dynamic model to study forest carbon storage in this paper. Although dynamic models based on change data are currently less commonly used in forestry, their advantages in terms of accuracy have been demonstrated [11,12]. In this study, although the effects of dynamic models and static models are not directly compared through experiments, we can refer to some of the existing studies related to forest carbon storage modelling in recent years [8,67,68]. The R2 of the remote sensing static models developed in these studies ranges from 0.64 to 0.73 as the highest value, while the R2 of the optimal model in this study is 0.87. Thus, the advantages of the dynamic model in terms of accuracy can be seen.

The DI is an area evaluation indicator. When the distribution curve is flatter, it indicates that the distribution of such changes deviates less from the standard distribution, and its adaptability to the terrain is wider [16]. In this study, the areas of increasing or decreasing *Pinus densata* carbon storages were calculated separately, and the main topography of these two areas was evaluated according to the DI. The results of this study demonstrate that the increase in carbon storage in the periods ΔCS4, ΔCS5, and ΔCS6 is much greater than the decrease, leading to an increase in total carbon storage from 2002 to 2017, which is related to the policy of "returning farmland to the forest" that has been implemented in Shangri-La since 2000. In the area of increasing carbon storage, the distribution curves of these three periods are relatively flat, showing a general adaptation to different levels of TNI and suggesting that the implementation of this government policy has led to the widespread planting of *Pinus densata* in Shangri-La, thus resulting in a more balanced topographic distribution of increasing carbon storage across the region.

Different TNI levels reflect different elevation and slope conditions, and the dominant distribution reflects the main topographic range of carbon storage changes. Although the dominant distribution of carbon storage changes over the TNI gradient can be clearly seen throughout a particular period of change, the six periods of change in carbon storage are more complex, and the distribution trends across the periods of change are not entirely consistent along the TNI gradient. What can be seen from the dominant distribution is that the regions where carbon storage increased are mainly located at low altitudes or on gently sloping terrain or both, as such terrain is conducive to artificial tending. In contrast, in regions with decreased carbon storage, the dominant distribution occurs across the TNI gradient, and the overall dominance of the terrain is not obvious, which is due to the randomicity of deforestation or forest destruction. In the past, TNI has not been applied to separate forest land class studies, but this study shows that its combination with the DI can be used to evaluate the increase and decrease effects of a single land type.

Some studies have suggested that human activities are the main driver of forest carbon [69], but as Shangri-La is a highland region with a low population density, the impact of human life on the surrounding forests is also relatively low, and it is environmental factors, including topography, that mainly impact forest carbon storage. Therefore, it is essential to consider topographic factors when studying forest carbon storage in this region. The results of this paper show that adding TNI to the dynamic model can reflect the topographic effect of the model and enable the improvement of the dynamic model accuracy, and by analyzing the distribution of forest carbon storage changes along the TNI gradient, it is possible to reveal the changes in forest carbon storage under different topographic conditions.

At present, the main approach to the study of historical forest dynamics is still based on the estimation of forest biomass or carbon storage over several single years, and the values of changes in different periods are calculated in order to derive forest dynamics results [70]. The level of periodic change in the forest can be estimated directly from models developed from interval change data, and due to the advantages of dynamic models in terms of accuracy, the resulting change values are estimated more precisely. Different dynamic models can be obtained for different types of change data. In this paper, we only compared two dynamic models for the 5-year interval change data and annual average change data, and more types of dynamic models are yet to be explored.

Although it was difficult to derive carbon storage results for a single year in this study without the support of carbon storage inventory data or graphs related to the study years, the results of this study demonstrate that the dynamic model is able to estimate the change in forest carbon storage and is well suited for the study of forest carbon storage change. Dynamic models are capable of quantitatively describing changes in specific properties, so they need not be limited to the field of forestry research but have the potential for application in other areas of natural or manufactured landscapes where regular variation exists. This study has so far only considered the influence of topographic factors, and there remains the possibility of other environmental factors, such as climate and soil, also influencing the model, which could be the next direction for future research.

#### **5. Conclusions**

In this study, the carbon storage dynamic model was developed based on NFI data and Landsat time series images over a 30-year period. PLSR, RF, and GBRT were used for modelling, and the accuracy of the two dynamic models was compared; then TNI, which represents topographic factors, was introduced into the model, and the distribution of carbon storage changes in Shangri-La on different TNI gradients was analyzed. The main conclusions are as follows: (1) the model effects of the non-parametric methods, RF and GBRT, are much better than those of the parametric method, PLSR, and the accuracy of the 5-year interval change model is better than that of the annual average change model; (2) TNI can improve the accuracy of the dynamic carbon storage estimation model, and the accuracy of the dynamic model with RF is the highest after adding TNI; (3) the dynamic model displays good performance regarding the estimation of carbon storage changes, and the results of the interval change model show that the total carbon storage of *Pinus densata* in Shangri-La decreased in 1987–1992 and 1997–2002, and increased in 1992–1997, 2002–2007, 2007–2012, and 2012–2017; and (4) the DI can be used to evaluate the main topography for the regions that display a change in *Pinus densata* carbon storage, the predominant topography in the regions of increasing carbon storage is that of low elevations or low slopes, or a combination of both conditions, while the topography in regions of decreasing carbon storage is more random. The results of this study can be used as a reference for forest carbon storage estimation using Landsat images in areas with complex topography.

**Author Contributions:** Y.L. participated in the data analysis and wrote the draft article, J.Z. helped with the idea and draft revision, R.B. helped with the data analysis, and D.X. and D.H. participated in the draft revision. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the National Natural Science Foundation of China (Nos. 31860207 and 32260390), and the "Young Top Talents" special project of the high-level talent training support program of Yunnan province, China, in 2020 (No. YNWR-QNBJ-2020-164).

**Data Availability Statement:** The data presented in this study are available on request from the corresponding author. The data are not publicly available due to the confidentiality of the NFI dataset.

**Acknowledgments:** We would like to thank the editor and anonymous reviewers for their comments, which helped improve the manuscript. We also would like to acknowledge all the others individuals who contributed to this paper.

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

#### **References**

## *Article* **Estimation of Urban Forest Characteristic Parameters Using UAV-Lidar Coupled with Canopy Volume**

**Bo Zhang 1,2,3, Xuejian Li 1,2,3, Huaqiang Du 1,2,3,\*, Guomo Zhou 1,2,3, Fangjie Mao 1,2,3, Zihao Huang 1,2,3, Lv Zhou 1,2,3,4, Jie Xuan 1,2,3, Yulin Gong 1,2,3 and Chao Chen 1,2,3**

**Abstract:** The estimation of characteristic parameters such as diameter at breast height (DBH), aboveground biomass (AGB) and stem volume (V) is an important part of urban forest resource monitoring and the most direct manifestation of the ecosystem functions of forests; therefore, the accurate estimation of urban forest characteristic parameters is valuable for evaluating urban ecological functions. In this study, the height and density characteristic variables of canopy point clouds were extracted as Scheme 1 and combined with the canopy structure variables as Scheme 2 based on unmanned aerial vehicle lidar (UAV-Lidar). We analyzed the spatial distribution characteristics of the canopies of different tree species, and multiple linear regression (MLR), support vector regression (SVR), and random forest (RF) models were used to estimate the DBH, AGB, and V of urban single trees. The estimation accuracy of different models was evaluated based on the field-measured data. The results indicated that the model accuracy of coupling canopy structure variables (R2 = 0.69–0.85, rRMSE = 9.87–24.67%) was higher than that of using only point-cloud-based height and density characteristic variables. The comparison of the results of different models shows that the RF model had the highest estimation accuracy (R2 = 0.76–0.85, rRMSE = 9.87–22.51%), which was better than that of the SVR and MLR models. In the RF model, the estimation accuracy of AGB was the highest (R2 = 0.85, rRMSE = 22.51%), followed by V, with an accuracy of R2 = 0.83, rRMSE = 18.51%, and the accuracy of DBH was the lowest (R2 = 0.76, rRMSE = 9.87%). The results of the study provide an important reference for the estimation of single-tree characteristic parameters in urban forests based on UAV-Lidar.

**Keywords:** urban forest; UAV-Lidar; canopy volume; diameter at breast height (DBH); aboveground biomass (AGB); stem volume (V)

#### **1. Introduction**

Urban forests are an important part of urban ecosystems, and they are the foundation and guarantee of urban sustainable development [1,2]. They can effectively reduce the urban heat island effect and improve air quality and other environmental conditions as well as ecosystem services [3,4]. The single tree is the basic unit of the forests [5]; their characteristic parameters such as tree height (H), diameter at breast height (DBH), crown width, aboveground biomass (AGB), and tree volume (V) can effectively reflect the growth status, spatial distribution, and structural characteristics of forest resources, which are important elements of forest resource investigation and reliable diversity indicators of forest succession stages [6], as well as being the focus of research on urban ecosystems and their functions [7].

**Citation:** Zhang, B.; Li, X.; Du, H.; Zhou, G.; Mao, F.; Huang, Z.; Zhou, L.; Xuan, J.; Gong, Y.; Chen, C. Estimation of Urban Forest Characteristic Parameters Using UAV-Lidar Coupled with Canopy Volume. *Remote Sens.* **2022**, *14*, 6375. https://doi.org/10.3390/rs14246375

Academic Editor: Klaus Scipal

Received: 6 November 2022 Accepted: 13 December 2022 Published: 16 December 2022

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

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

Single trees in cities are highly fragmented and unevenly distributed in urbanized areas with high population densities and high concentrations of artificial landscapes [8], so it is time-consuming and labor-intensive to obtain information on the characteristic parameters of urban single trees through traditional forest resource investigation methods [9]. Remote sensing technology has rapidly developed and can be used to quickly and accurately obtain multiscale and multitemporal information on forest structure characteristics, effectively making up for the shortcomings of traditional forest resource monitoring methods and greatly improving the work efficiency [10]. However, passive remote sensing technology only provides spectral and textural information of the forest canopy surface, which is susceptible to atmospheric conditions and other factors, and it is difficult to obtain the three-dimensional structure of the vegetation canopy [11,12].

Light detection and ranging (LiDAR) is an active remote sensing technology that obtains the distance between a sensor and a target by calculating the time difference between the laser pulse emitted by the sensor and the received echo pulse. Because LiDAR has strong penetration into the forest, it can accurately obtain the three-dimensional structure information of forest tree height and canopy [13–15], thus realizing the leap from two-dimensional to three-dimensional forest canopy structure information, resulting in the emergence of canopy parameters extracted from lidar data as a hot research topic globally. Previous studies have shown that LiDAR-based estimation of diameter at breast height, biomass, stock volume, and forest distribution mapping can be effective [16–19]. For example, Cao et al. [20] used full-waveform unmanned aerial vehicle (UAV) lidar data to extract point cloud metrics and waveform metrics calculated based on voxel-based methods to estimate the single tree AGB of plantation forests in the coastal region of east China. The results indicated that full-waveform lidar data can effectively estimate the AGB of single trees. Liu et al. [21] used the constant allometric ratio model to estimate the forest single tree biomass based on UAV lidar data to obtain single tree canopy characteristic parameters (tree height, crown width, canopy projection area, and canopy volume) and achieved a good fit and high prediction accuracy. Qin et al. [22] used UAV lidar to estimate the subtropical single tree carbon stock in Shenzhen, southern China, and the results indicated that the height variable can explain the variation of tree carbon stock and estimate the single tree carbon stock well. Therefore, it is important to examine the application of UAV lidar for estimating single tree characteristic parameters in urban areas.

Forest canopy structure includes the horizontal and vertical directions of branches and leaves, canopy width and height, and canopy light transmission [23,24]. The canopy width and cross-sectional area can be used to measure the horizontal extension size of the canopy. The vertical structure of the canopy is mainly the spatial distribution and hierarchical characteristics of the forest vegetation, and the performance is comparatively complex [25–27]. To quantify the three-dimensional structure of the forest canopy, Lefsky et al. proposed a voxel-based canopy volume model (CVM) to characterize the differences in the volume and vertical spatial distribution of the canopy [28]. The basic principle of the CVM model is reflecting the spatial heterogeneity of the forest structure arising from the difference in the light environment within the canopy by dividing the canopy into two parts, the photosynthetically active zone and the inactive zone, to realize the spatial arrangement of elements within the canopy structure and the distinction of volume structure [29–32]. Therefore, the CVM model is an important method for obtaining the parameters of tree canopy structure.

In summary, the estimation of DBH, AGB, and stem volume of trees is an important element of urban forest resource monitoring and the most direct manifestation of the ecosystem functions of forests, and lidar is an advanced technical tool for detecting the three-dimensional structure of forests. Therefore, in this study, three urban tree species, Ginkgo (*Ginkgo biloba* L.), Cinnamomum camphora (*Cinnamomum camphora (Linn.) Presl*), and Metasequoia glyptostroboides (*Metasequoia glyptostroboides Hu et Cheng*) were used as examples. UAV-Lidar data were used to obtain the canopy point clouds of the three kinds of single trees, coupled with point clouds and canopy structure variables, and three

methods, multiple linear regression (MLR), support vector regression (SVR) and random forest (RF), were used to establish models for estimating the DBH, AGB, and V of single trees in urban forests based on UAV-Lidar data. The model and estimation results were validated using ground-measured data. This provides an important technical tool for rapid and accurate monitoring of single tree parameters in urban forests.

#### **2. Materials and Methods**

*2.1. Study Area*

The study area is in Lin'an (29◦56 to 30◦23 N latitude and 118◦51 to 119◦52 E longitude) (Figure 1), Hangzhou City, Zhejiang Province. Lin'an belongs to the subtropical monsoon climate with warm and humid conditions, abundant light and rainfall, and four distinct seasons. The average annual temperature is 16.4 ◦C, the frost-free period is 237 days, the sunshine time is 1847.3 h, and the annual precipitation is 1613.9 mm. The area is dominated by hills and mountains, the terrain inclines from northwest to southeast, and the three-dimensional climate is obvious. The climax vegetation is subtropical evergreen broad-leaved forest, and the main tree species planted in the urban area of Lin'an City include *Metasequoia glyptostroboides* (*M. glyptostroboides*), *Ginkgo biloba* (*G. biloba*), *Cinnamomum camphora* (*C. camphora*), etc.

**Figure 1.** Overview of the study area: (**a**) location of Lin'an, (**b**) location of the study area, (**c**) UAV-Lidar point clouds of the study area, (**d**) Metasequoia glyptostroboides point cloud profile, (**e**) Ginkgo biloba point cloud profile, (**f**) Cinnamomum camphora point cloud profile.

#### *2.2. Field Measurements*

In July 2021, the DBH, tree height, crown height, and crown width of 64 stems of *G. biloba*, 74 stems of *C. camphora*, and 55 stems of *M. glyptostroboides* in the study area were measured in detail, as shown in Figure 1. The single trees were positioned with Huace Smart Real Time Kinematic. The AGB consisted of stem biomass (WS), branch biomass (WB), and foliage biomass (WF). In this study, the biomass of each component was calculated according to the biomass allometric equations of different tree species and summed to obtain the AGB of a single tree [33–35], as shown in Table 1. The stem volume was calculated according to the single-entry stem volume table of Zhejiang Province and the measured single tree DBH data [36]. Table 2 shows the statistical characteristics of the parameters of the three tree species.

**Table 1.** Biomass allometric equations for each biomass component of the three tree species.

Note: H is tree height (m), DBH is diameter at breast height (cm).

**Table 2.** Summary of information on measured characteristics parameters of the three tree species.

#### *2.3. Lidar Data*

The DJI Matrice 600 Pro six-rotor UAV with a lightweight Velodyne Puck LITE™ laser scanner was used to acquire the original lidar point clouds in the study area (Figure 2). The flight height of the UAV is 60 m above ground level, with a flight speed of 8 m/s, a route spacing of 25 m, and a lateral overlap rate of data sampling of 50%. The sensor records the first echo information of the laser pulse with a wavelength of 903 nm, a maximum scanning angle of ±15◦, a scanning frequency of 20 Hz, and a scanning speed of 300,000 points/s. The final average point cloud density obtained is approximately 230 points/m2.

**Figure 2.** UAV and LiDAR system.

#### 2.3.1. Lidar Data Preprocessing

The original lidar point cloud data were denoised using the height thresholding method, and the point cloud data after noise removal were filtered and separated into ground points and nonground points. First, the ground points were extracted by filtering with the improved progressive TIN densification algorithm [37]. Then, the height average of the laser points within a cell was calculated via the inverse distance weighting method to obtain a digital elevation model (DEM) with a spatial resolution of 0.5 m. Finally, the DEM was used to normalize the point cloud data to obtain the normalized point cloud data. In this study, the point cloud segmentation algorithm was used to segment individual trees based on normalized point cloud data [38,39]. This algorithm identified single trees via region growing combined with thresholding, and then identified the top of the tree to determine the distance between the surrounding points and the vertex, and expanded the region to segment the first tree. Successive iterations were made on this basis until all trees were segmented. The characteristic variables of a single tree were extracted based on the segmented single tree point cloud.

#### 2.3.2. Lidar Metrics

The characteristic variables based on the lidar data can be used to estimate the forest characteristic parameters, and the point cloud characteristics extracted from the first returns have a remarkable correlation with the height, which is more suitable for estimating the forest characteristic parameters [13,40]. Of course, to reduce the influence of low ground vegetation on the data, the data after filtering the point clouds below 2 m were used as the crown point clouds, and characteristic variables were extracted from the first returns of the lidar point cloud [41,42]. In this study, the lidar data characteristic variables included: height-based metrics (HB) describing the parameters related to the lidar point cloud height; density-based metrics (DB) describing the canopy return density variable, which is the ratio of the number of height percentile point clouds to the total number of point clouds; the canopy area (S), which is the projected area of the canopy point cloud calculated based on the two-dimensional convex packet algorithm; the crown diameter (CD), which is the average of the east–west and north–south crown diameters of the point clouds. The canopy volume variables include OG, CG, EV, and OV. The metrics and descriptions are shown in Table 3.

#### 2.3.3. Calculation of the Volume of Single Tree Canopy

Urban trees are frequently pruned and truncated, [43], making the crown of the pruned single tree change, often with a special crown shape. Therefore, canopy volume was extracted as a metric for single tree characteristic parameter estimation, as shown in Table 3. As shown in Figure 3, in this study, the voxel-based canopy volume method was used to calculate the canopy volume metrics for lidar point clouds of tree crowns:

**Figure 3.** Illustration of the voxel-based canopy volume model.

The space where the canopy point clouds was located was divided into 0.5 × 0.5 × 0.5 m voxels [44], which were divided into vertical columns, and each column was further layered into four canopy structures. First, each voxel was classified as "filled" or "empty" according to whether there was a point cloud in the voxel, that is, the volume in unit area (m3/m2).

**Table 3.** Description of metrics derived from lidar data.

Then, according to the distribution of the filled position, the upper 65% of the filled zone was defined as "euphotic", in every column of voxels, and the remaining 35% was defined as "oligophotic". According to the spatial distribution and location of the empty voxel, in each voxel column, the empty voxels between the top of the canopy and the first filled voxels were defined as the "open gap", and the empty voxels between the filled voxels and the ground were called the "closed gap". The three-dimensional canopy volume distribution was converted into a two-dimensional canopy volume profile (CVP) according to the percentage of the volume of the four classified canopy volume characteristics in each height interval. The canopy volume distribution indicated the distribution of elements arranged in the vertical spatial extent of the canopy [30].

#### *2.4. Model Construction Methods and Scheme*

In this study, three modeling methods, MLR, SVR and RF, were used to construct the estimation models of urban single tree characteristic parameters based on the obtained lidar data characteristic variables. To study the influence of canopy structure on the accuracy of single tree characteristic parameters, the models were constructed in two schemes. The model excluding canopy structure variables is referred to as "Scheme 1", and the model including canopy structure variables is referred to "Scheme 2".

#### 2.4.1. MLR Model

MLR is the most commonly used parameterization method for estimating forest characteristic parameters from remote sensing information, and can quickly establish a linear relationship between two or more independent variables and dependent variables to achieve parameter estimation. The MLR is generally expressed as follows:

$$Y = a\_0 + a\_1 \mathbf{x}\_1 + a\_2 \mathbf{x}\_2 + \dots + a\_n \mathbf{x}\_n \tag{1}$$

where *a*0 is the constant term; *a*1, *a*2, ... , *an* are regression coefficients representing the degree of contribution of the respective variables to the dependent variable; and *x*1, *x*2, ... , *xn* are the independent variables, which are the characteristic variables shown in Table 3. *Y* is the dependent variable, which is the estimated characteristic parameter of this study. In this study, all possible combinations of variables were evaluated using "all subsets" regression, and the best combination of variables was selected to build the MLR model to estimate the three characteristic parameters [48].

#### 2.4.2. SVR Model

SVR is a machine learning model that uses support vector machines to perform regression analysis [49,50]. It applies classification methods to solve regression problems with finite samples, mainly based on a given sample data set, by seeking a function to fit all sample points so that the total variance of sample points from the hyperplane is minimized [51]. SVR transforms the nonlinear problem into a linear problem in high-dimensional space via kernel functions for nonlinear separable samples in low-dimensional input space, replacing the inner product operation in high-dimensional space, and ensuring good generalization ability [52]. Therefore, the SVR model has high accuracy, good ability to handle high-dimensional and small sample data, good generalization ability, and robustness. In this study, four kernel function models, including linear, polynomial, radial basis function (RBF), and multilayer perceptron (Sigmoid), were used, and the best penalty coefficient (C) with gamma value (g) was selected via grid search cross-validation to construct the SVR model to estimate the three characteristic parameters.

#### 2.4.3. RF Model

The random forest algorithm is another commonly used machine learning method. The algorithm is based on modified nonparametric modeling of decision trees [53], and constructs a decision tree by bootstrapping from the original sample set with put-back randomly selected N samples to predict the results. The RF algorithm has good noise resistance and can handle high-dimensional data with relatively high prediction accuracy. An unbiased estimate of the error can be generated during the RF calculation, and the importance of each variable involved in the model can be evaluated. There are three important parameters in the estimation of single tree characteristic parameters using the RF algorithm: Mtry is the number of variables used randomly at the nodes of each tree, and Ntree is the number of regression trees in the RF. Nodesize is the minimum number of terminal nodes in the regression analysis, and the default value is 5 [54,55]. In this study, three characteristic parameters estimation models are developed based on the optimization of RF parameters.

#### *2.5. The Flow Chart and Accuracy Validation*

The flow chart of this study is shown in Figure 4. First, field measurements, lidar data processing, and characteristic variable extraction were conducted. Second, two-thirds of all measured samples were selected into training samples, and one-third were divided into test samples. Finally, three modeling methods, MLR, SVR and RF, were used to construct single tree characteristic estimation models according to Scheme 1 and Scheme 2, and the accuracy of the models was evaluated.

**Figure 4.** The flow chart of this study.

The model accuracy evaluation metrics include the determination coefficient (R2), root mean square error (RMSE), and relative root mean square error (rRMSE). Generally, higher values of R2 and lower values of RMSE and rRMSE indicate better performance of the model. R2, RMSE, and rRMSE are calculated as follows:

$$\mathbf{R}^2 = 1 - \frac{\sum\_{i=1}^n \left(\mathbf{x}\_i - \mathbf{x}\_i\right)^2}{\sum\_{i=1}^n \left(\mathbf{x}\_i - \overline{\mathbf{x}}\right)^2} \tag{2}$$

$$\text{RMSE} = \sqrt{\frac{1}{n} \sum\_{i=1}^{n} \left(\mathbf{x}\_i - \hat{\mathbf{x}}\_i\right)^2} \tag{3}$$

$$\text{rRMSE} = \frac{\text{RMSE}}{\overline{\text{x}}} \times 100\% \tag{4}$$

where *n* is the number of samples, *x* is the measured value of the sample canopy characteristic parameters, *x* is the mean value of the sample canopy characteristic parameters, and *x*ˆ*i* is the predicted value of the sample canopy characteristic parameters.

#### **3. Results**

#### *3.1. Canopy Volume and Profile Analysis*

The canopy volume profile can directly show the spatial heterogeneity of the forest canopy structure and the distribution and change in elements in the vertical direction. Comparing the canopy volumes and profiles of the three species, the euphotic volume was significantly larger than the oligophotic volume in the filled volume of the *G. biloba* canopy (Figure 5); the closed gap volume was larger than the open gap volume in the empty volume, and the closed gap volume occupied the largest volume. The canopy distribution of *C. camphora* was consistent with that of *G. biloba* (Figure 6), and the euphotic volume was greater than the oligophotic volume; the closed gap volume was larger than the open gap volume in the empty volume, and the closed gap volume occupied the largest volume. In the canopy distribution of *M. glyptostroboides* (Figure 7), the euphotic volume was close to the oligophotic volume; the closed gap volume was slightly larger than the open gap volume in the empty volume.

**Figure 5.** (**a**) *G. biloba* point cloud, (**b**) canopy volume distribution, which shows the distribution of canopy structure classes after the expansion of all columns in the panel; and (**c**) canopy volume profile, which shows the volume percentage of each class of total volume in each height interval.

**Figure 6.** (**a**) *C. camphora* point cloud, (**b**) canopy volume distribution, which shows the distribution of canopy structure classes after the expansion of all columns in the panel, and (**c**) canopy volume profile, which shows the volume percentage of each class of total volume in each height interval.

**Figure 7.** (**a**) *M. glyptostroboides* point cloud, (**b**) canopy volume distribution, which shows the distribution of canopy structure classes after the expansion of all columns in the panel; and (**c**) canopy volume profile, which shows the volume percentage of each class of total volume in each height interval.

#### *3.2. Variable Importance Analysis*

In this study, the point cloud characteristic variables extracted from lidar data were run 100 times using the RF model according to two modeling schemes. Figures 8 and 9 show the importance scores of Scheme 1 and Scheme 2 input variables, respectively. As shown in Figure 8, the importance score of the height characteristic variable in Scheme 1 is significantly higher than that of the density characteristic variable. Among the three single tree characteristic parameter models, Hcv was the largest in DBH and AGB, with values of 27.89 and 44.05, respectively; Hsq was the largest in V, with a value of 21.12. As shown in Figure 9, the three parameters with the highest importance in Scheme 2 were all canopy structure variables. In the DBH and V models, CG had the highest influence, with 31.96 and 28.97, respectively, and in the AGB model, CD had the highest influence, with 40.57. The results show that the canopy structure variable was significantly more important than other variables and was the key variable in the estimation of single tree structure parameters.

**Figure 8.** Importance of input variables based on Scheme 1 (%InMSE: the percentage increase in the mean square error, ((**a**): DBH, (**b**): AGB, (**c**): V)).

**Figure 9.** Importance of input variables based on Scheme 2 (%InMSE: the percentage increase in the mean square error, ((**a**): DBH, (**b**): AGB, (**c**): V)).

Figure 10 shows the path coefficients (absolute values) calculated using the structural equations for the direct effect of the three characteristic variables of height characteristics, density characteristics, and canopy structure on the characteristic parameters. As shown in Figure 10, CS had the largest direct effect on DBH and AGB, while HD had the largest direct effect on V, and DB had the smallest effect on characteristic parameters. For DBH, the path coefficient values, in order, were as follows: CS (0.647) > HD (0.318) > DB (−0.095). For AGB, the path coefficient values were as follows: CS (0.871) > HD (0.116) > DB (−0.101). For V, the path coefficient values were as follows: HD (0.593) > CS (0.389) > DB (0.038). The canopy structure variables had important influence on the estimation of the three single tree characteristics parameters.

**Figure 10.** Relative importance of input variables to characteristic parameters based on Scheme 2 ((**a**): DBH, (**b**): AGB, (**c**): V).

#### *3.3. Model Construction and Evaluation*

#### 3.3.1. MLR Estimation Results

To evaluate the accuracy of the established prediction models of single tree characteristic parameters, three independent variables were selected to construct the MLR models in this study (Table 4). Among all the metrics selected in the MLR model of Scheme 1, the height percentiles (Hcv, Hmean, Hsq, Hstd, H99, and Hvar) were frequently selected by the models. Among all the metrics selected by the MLR model in Scheme 2, height percentile (Hcm, H95, Hcv and Hvar), canopy structure variables (CD, CG, EV) were frequently selected by the models. The accuracy value R2 of the models was improved after coupling the canopy

structure variables. The canopy structure variables played an important role in the construction of the model, indicating that these variables are sensitive in estimating forest structure.

**Table 4.** The MLR prediction models and their accuracy assessment under different schemes.

Figures 11 and 12 show the correlations between the predicted and field-measured values of characteristic parameters estimated using the MLR models of Scheme 1 and Scheme 2, respectively. As shown in Figure 11, the R2 values of the model training accuracy for single tree DBH, AGB, and V in Scheme 1 were 0.76, 0.78, and 0.83, respectively, with rRMSE = 9.67 to 27.23%; the R2 values of the model testing accuracy were 0.64, 0.69, and 0.62, respectively, with rRMSE = 11.02 to 39.01%. As shown in Figure 12, the training and testing accuracy of Scheme 2 single tree characteristic parameters were improved. The R2 values of the model training accuracy of DBH, AGB, and V were 0.78, 0.85, and 0.84, respectively, with rRMSE = 9.33 to 30.82%. The R2 values of the model testing accuracy were 0.69, 0.82, and 0.80, respectively, with rRMSE = 10.41–24.67%. The comparison of the results shows that the R2 values of the estimation accuracy of the characteristic parameters of Scheme 2 were all improved after adding the canopy structure variable. Among them, the estimation accuracy of V improved by 29.0%, which was the greatest improvement, followed by AGB with an accuracy improvement of 18.8% and DBH with an accuracy improvement of 7.8%. In addition, the rRMSE values of the AGB and V estimation results also decreased by 36.8% and 27.7%, respectively, which were large decreases.

**Figure 11.** Correlation between measured and estimated values of single tree DBH (**a**), AGB (**b**), and V (**c**) of the MLR model based on Scheme 1.

**Figure 12.** Correlation between measured and estimated values of single tree DBH (**a**), AGB (**b**), and V (**c**) of the MLR model based on Scheme 2.

#### 3.3.2. SVR Estimation Results

Figure 13 shows the comparison of the R2 values of the estimated characteristic parameters for the four kernel function SVR models under the two modeling schemes. As shown in Figure 13, the R2 values of the single tree characteristic parameters in both Scheme 1 and Scheme 2 were the highest for the SVR model with RBF as the kernel function. Therefore, the RBF kernel function was chosen to construct the SVR model in this study.

**Figure 13.** Results of single tree characteristic parameters using different kernel functions ((**a**): Scheme 1, (**b**): Scheme 2).

Figures 14 and 15 show the correlations between the predicted and field-measured values of characteristic parameters estimated using the SVR models of two modeling schemes, respectively. As shown in Figure 14, the R2 values of the model training accuracy for single tree DBH, AGB, and V in Scheme 1 were 0.83, 0.84, and 0.85, respectively, with rRMSE = 8.14 to 24.10%; the R2 values of the model testing accuracy were 0.67, 0.77, and 0.75, respectively, with rRMSE = 10.76 to 26.78%. As shown in Figure 15, the R2 values of the model testing accuracy of DBH, AGB, and V were 0.72, 0.82, and 0.81, respectively, with rRMSE = 10.24 to 23.37%. The comparison of the results shows that the R2 values of the training and testing accuracy of single tree characteristic parameters of the SVR model were all improved by adding canopy structure variables in Scheme 2. Among them, the estimation accuracy of V improved by 8%, which was the greatest improvement, followed by DBH with an accuracy improvement of 7.4% and AGB with an accuracy improvement of 6.5%. In addition, the rRMSE values of the AGB and V estimation results also decreased by 12.7% and 11.6%, respectively, which were large decreases, indicating that coupling canopy parameters can improve the estimation accuracy of urban single tree parameters.

**Figure 14.** Correlation between measured and estimated values of single tree DBH (**a**), AGB (**b**), and V (**c**) of the SVR model based on Scheme 1.

**Figure 15.** Correlation between measured and estimated values of single tree DBH (**a**), AGB (**b**), and V (**c**) of the SVR model based on Scheme 2.

#### 3.3.3. RF Estimation Results

The training data were input into the random forest model to traverse all variable values and eventually obtain the optimal parameters. Figure 16 shows Mtry, which was used to determine the minimum variable for each tree in the RF model, and the minimum Mtry value was required when the model error was minimal. As shown in Figure 17, the RMSE of the model error tended to be stable after Ntree reached 1500. Therefore, the values of Ntree in the optimized random forest model were set to 1500 in this study. Table 5 lists the specific settings for different parameter values for the two schemes.

**Figure 16.** Influence of Mtry on model error ((**a**): DBH, (**b**): AGB, (**c**): V).

**Figure 17.** Influence of Ntree on model error ((**a**): DBH, (**b**): AGB, (**c**): V).

**Table 5.** Results of the optimization of model parameters for different schemes.

Figures 18 and 19 show the correlations between the predicted and field-measured values of characteristic parameters estimated using the RF models of Scheme 1 and Scheme 2, respectively. As shown in Figure 18, the R2 values of the model training accuracy for single tree DBH, AGB, and V in Scheme 1 were 0.74, 0.80, and 0.81, respectively, with rRMSE = 10.46 to 26.36%; the R2 values of the model testing accuracy were 0.67, 0.74, and 0.76, respectively, with rRMSE = 10.95 to 29.69%. As shown in Figure 19, the training and testing accuracy of Scheme 2 single tree characteristic parameters were improved. The R2 values of the model training accuracy of DBH, AGB, and V were 0.80, 0.89, and 0.86, respectively, with rRMSE = 8.9 to 20.77%. The R2 values of the model testing accuracy were 0.76, 0.85, and 0.83, respectively, with rRMSE = 9.87–22.51%. Comparing the RF model results, the R2 values of the estimation accuracy of the characteristic parameters of Scheme 2 were all improved after adding the canopy structure variable. Among them, the estimation accuracy of AGB improved by 14.9%, which was the greatest improvement, followed by DBH with an accuracy improvement of 13.4% and V with an accuracy improvement of 9.2%. In addition, the rRMSE values of the AGB and V estimation results also decreased by 24.2% and 15.2%, respectively, which were large decreases.

**Figure 18.** Correlation between measured and estimated values of single tree DBH (**a**), AGB (**b**), and V (**c**) of the RF model based on Scheme 1.

**Figure 19.** Correlation between measured and estimated values of single tree DBH (**a**), AGB (**b**), and V (**c**) of the RF model based on Scheme 2.

#### *3.4. Comparison of Model Results*

Table 6 shows the summary of the field measured data and the estimation results of different model characteristic parameters. Compared with the field-measured data, the estimated CV values of MLR, SVR, and RF models range from 16.6 to 18.11%, the estimated CVs of AGB range from 49.15 to 54.22%, and the estimated CVs of V range from 44.67 to 46.55%, with smaller variations and mean values closer to the measured values. Appendix A (Table A1) shows the comparison of the model training accuracy and testing accuracy of three models for estimating three single tree parameters in two schemes. An analysis of the results in Appendix A (Table A1) shows that all models of the two schemes achieved higher accuracy estimation of single tree parameters in urban forests. Figure 20 shows the distribution of the normalized residuals for the testing phase of the characteristic parameters for different models. The normalized residuals of the characteristic parameters of the testing samples of the three models were in the range of −2 to 2, indicating that all models had good stability and reliability in predicting the characteristic parameters of a single tree.

**Table 6.** Summary of information on measured and predicted characteristic parameters of the samples.

However, for the two schemes, the model accuracy and testing accuracy of the single tree parameters for all models of Scheme 2 were improved, and the error was decreased compared with Scheme 1, indicating that the coupled canopy parameters could improve the estimation accuracy of the urban single tree parameter. In addition, from the three models, the performance of the two machine learning models was better than that of the MLR model, where the training accuracy of single tree characteristic parameters of the SVR model was slightly higher than that of the RF model, but the testing accuracy of single tree characteristic parameters of the RF model was better than that of the SVR model.

**Figure 20.** Distribution of normalized residuals for the testing phase of characteristic parameters for different models.

#### **4. Discussion**

The forest canopy is a key component of forests that affects ecosystem processes and functions [56], and the quantification and analysis of canopy distribution is one of the methods used to characterize the spatial structure of forests. In this study, voxel-based canopy volume was used to characterize the canopy spatial structure of different urban trees, and the canopy volume profile was derived to more intuitively reflect the spatial heterogeneity of the canopy structure and the variation in the arrangement of elements within the canopy. As shown in Figures 5–7, the open gap volume of *G. biloba* was the largest, while the closed gap volume of *C. camphora* was larger than that of both *G. biloba* and *M. glyptostroboides*. First, this had a direct influence on the structure and shape of the canopy, which was relatively regular for *C. camphora* and *M. glyptostroboides*, while the canopy shape of *G. biloba* was more complex and the distribution of branches and leaves was more dispersed, resulting in a larger volume share in the open gap of *G. biloba*. In addition, closed gap volume was also strongly related to proper pruning in cities, which causes higher crown base height [43], thus leading to a larger percentage of closed gap volume. On the other hand, coniferous trees allow more light penetration into the lower canopy compared to broadleaved trees [57], and *M. glyptostroboides* belongs to the coniferous species, thus allowing *M. glyptostroboides* to form more areas of oligophotic zone.

Lidar data can provide parameter information directly connected to forest canopy structure, and the estimation of forest characteristic parameters using different regression methods can produce satisfactory prediction results [58]. The results of this study show that the estimation accuracy of single tree characteristic parameters with coupled canopy structural variables was improved compared to using only height and density characteristic variables. The model accuracy was higher than the prediction accuracy of stand stock and aboveground biomass in urban broadleaf forest areas estimated using ALS data [59]. This is due to the fact that high-density ULS data had richer canopy structure information. Therefore, high-density lidar data were more advantageous in estimating stand volume and aboveground biomass.

This research shows that the performance of the SVR and RF machine learning models was better than that of the MLR model. As a statistical regression model, MLR is not suitable for representing the complexity of high data and is sensitive to noise, and the MLR model is often prone to underfitting, making the model performance poor. SVR seeks linear regression hyperplanes and solves nonlinear problems in low-dimensional spaces by mapping kernel functions from low to high dimensions. Although the SVR model fits the training data well, there may be overfitting of the model. In addition, the SVR model also needs to find the optimal penalty coefficient and gamma value to obtain the optimal model [60,61]. The RF was able to handle high-dimensional data and had good noise resistance. During the model operations, unbiased estimates of errors can be obtained, and the importance of each variable can be evaluated in the RF model [62]. The accuracy of training the single tree characteristic parameters of the SVR model in this study was slightly higher than that of the RF model, but the validation accuracy was not as good as that of the RF model in the study, which may be caused by this reason.

Related studies around the world also indicate that RF has good predictive ability in forest parameter estimation [58,60,63]. For example, Zhou et al. [64] used a RF model to estimate the AGB of urban single trees based on UAV lidar data and achieved high estimation accuracy. Zhang et al. [65] combined lidar and high-resolution remote sensing images by comparing different models (SLR, LNN, BPNN, SVR, RF) for the quantitative estimation and inversion of biomass, and the results indicated that the RF model had the highest fitting accuracy. Cao et al. [63] indicated that the accuracy of the RF model was higher than that with SVR, backpropagation neural networks, k-nearest neighbor, and the generalized linear mixed model in the remote sensing estimation of forest biomass based on satellite remote sensing. Peng Xi et al. [66] established different models for estimating the s characteristic parameters of tropical forests in China based on UAV lidar data, and the study indicated that the RF model had good accuracy in estimating forest characteristic parameters. Although this study provides a reference for the application of UAV lidar in urban forest characteristic parameter estimation, there are still some limitations. The urban forest is unevenly distributed, and the estimation of large-scale urban forest characteristic parameters using UAV lidar is still a challenge.

#### **5. Conclusions**

In this study, we used UAV lidar to obtain three kinds of single tree canopy point clouds coupled with point cloud and canopy structure variables. MLR, SVR, and RF models were used to estimate the characteristic parameters of DBH, AGB, and V of single trees in urban forests based on UAV-Lidar data. The results indicate that canopy volume profiles can visualize the spatial heterogeneity of forest canopy structure and the distribution and variation of elements in the vertical direction, and canopy structure variables such as CG, OV, EV, S, and CD had important effects on single tree characteristic parameters. The model training accuracy and testing accuracy of the single tree parameters of the MLR, SVR, and RF models were improved by incorporating canopy structure variables. In comparison, the two machine learning models, SVR and RF, outperformed MLR, but the testing accuracy of single tree characteristic parameters of the RF model was better than that of the SVR model. The results of the study provide an important reference for the estimation of single tree characteristic parameters in urban forests based on UAV-Lidar data, which is necessary and useful for urban managers to understand the functions and values of urban forests and to maximize the environmental benefits of urban forests.

**Author Contributions:** Conceptualization, H.D.; methodology, B.Z.; validation, B.Z.; formal analysis, B.Z. and Z.H.; investigation, B.Z., L.Z., J.X., Y.G. and C.C.; data curation, B.Z.; writing—original draft preparation, B.Z.; writing—review and editing, X.L., F.M., G.Z. and H.D.; visualization, B.Z.; supervision, H.D. All authors have read and agreed to the published version of the manuscript.

**Funding:** The authors gratefully acknowledge the support of the National Natural Science Foundation of China (U1809208, 32171785), the Key Research and Development Program of Zhejiang Province (2021C02005), and the State Key Laboratory of Subtropical Silviculture (No. ZY20180201).

**Data Availability Statement:** Not applicable.

**Acknowledgments:** The authors gratefully acknowledge the support of various foundations. The authors are grateful to the editor and anonymous reviewers whose comments have contributed to improving the quality of this manuscript.

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

#### **Appendix A**

**Table A1.** Accuracy evaluation of the parametric model of a single tree characteristic for different schemes.

#### **References**

## *Article* **High-Resolution Remote Sensing Images Can Better Estimate Changes in Carbon Assimilation of an Urban Forest**

**Qing Huang 1, Xuehe Lu 2,\*, Fanxingyu Chen 2, Qian Zhang 3,4 and Haidong Zhang 5**

**Abstract:** Urban forests have the potential to sink atmospheric CO2. With the improvement of coverage of vegetation in urban environments, more attention has been paid to the carbon sequestration potential of the urban forest. However, the high fragmentation of urban forests makes it difficult to evaluate their carbon budget on a regional scale. In this study, the GPP-NIRv relationship model was employed to estimate GPP in Suzhou by MODIS, Landsat-8 and Sentinel-2 remote sensing data, and to further explore what kind of remote images can figure out the spatial-temporal pattern of GPP in urban forests. We found that the total GPP of the terrestrial ecosystem in Suzhou reached 8.43, 8.48, and 9.30 Tg C yr-1 for MODIS, Landsat-8, and Sentinel-2, respectively. Monthly changes of GPP were able to be derived by MODIS and Sentinel-2, with two peaks in April and July. According to Sentinel-2, urban forests accounted for the majority of total GPP, with an average of about 44.63%, which was larger than the results from GPP products with coarser resolutions. Additionally, it is clear from the high-resolution images that the decline of GPP in May was due to human activities such as the rotation of wheat and rice crops and the pruning of urban forests. Our results improve the understanding of the contribution of the urban forest to the carbon budget and highlight the importance of high-resolution remote sensing images for estimating urban carbon assimilation.

**Keywords:** gross primary productivity; near-infrared reflectance of vegetation; urban forest; carbon budget

#### **1. Introduction**

Gross Primary Productivity (GPP) quantifies the total amount of carbon assimilated by plants through photosynthesis per unit time. As the critical variable of carbon cycling of terrestrial ecosystems, GPP is the initial amount of energy and material entering the terrestrial ecosystem and plays an essential role in regulating the global carbon cycle [1–3]. GPP can be observed at the leaf level or ecosystem level [4–6], and also can be simulated by process-based terrestrial ecosystem models [7–9] or estimated through remote sensing [10–13]. Benefiting from the eddy-covariance techniques and the multisource data of remote sensing, satellite-based GPP models were developed to estimate regional and global GPP by establishing the empirical relationship between the vegetation index (VI) and surface observation data [14–18]. Previous studies have shown that the correlation coefficient between Enhanced Vegetation Index (EVI) and surface observation data in global GPP estimation can reach 0.52–0.92 [19], and the correlation coefficient between Solar-induced Chlorophyll Fluorescence (SIF) and GPP for crops and grass can achieve 0.87 [20]. However, there are several limitations in these methods; for example the adaptability of the empirical method is affected by different ecosystem structures and climate conditions [21].

**Citation:** Huang, Q.; Lu, X.; Chen, F.; Zhang, Q.; Zhang, H. High-Resolution Remote Sensing Images Can Better Estimate Changes in Carbon Assimilation of an Urban Forest. *Remote Sens.* **2023**, *15*, 71. https://doi.org/10.3390/rs15010071

Academic Editor: Arturo Sanchez-Azofeifa

Received: 13 November 2022 Revised: 19 December 2022 Accepted: 20 December 2022 Published: 23 December 2022

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

Additionally, the coarse spatiotemporal resolution of satellite SIF is generally a limiting factor in the regional estimation of GPP [22,23].

In recent years, Near-infrared Reflectance of Vegetation (NIRv), which integrates the advantages of NDVI and near-infrared reflectance, has provided a new avenue for GPP estimation [24–27]. NDVI can intuitively reflect vegetation coverage, which is tightly correlated to vegetation productivity, while near-infrared reflectance can embody the information on vegetative canopy structure and leaf area [24,28]. Therefore, NIRv takes into consideration both the pigment and structure of vegetation in remote sensing images to avoid the saturation effects of NDVI [29–31]. Meanwhile, a strong and stable correlation was found between NIRv and GPP in the same ecosystem [21,23,32]. In the study of regional and global carbon budgets of terrestrial ecosystems, NIRv has been suggested as the effective substitution for satellite SIF, based on theoretical derivations and radiative transfer simulations, and has gained a great deal of attention in recent years [24,33].

Urban forest is considered as a large green infrastructure in the urban area, which consists of trees, shrubs and grasses in streets, parks, gardens, alongside rivers and so on [34–36]. Additionally, urban forest in this study is also defined as the above description. Urban forest usually plays an indispensable role for people living in urban environments, such as reducing air pollution and heat islands effect [37], and increasing biodiversity [38]. As global CO2 concentrations continue to increase, there is also increasing interest in the carbon storage and sequestration of the urban forest [39–41]. Vegetation in urban areas is offered for environments with higher temperatures and CO2 concentrations, highly intensive use of water and pesticides [42,43], and its photosynthetic activity would play an important role in the carbon cycle of terrestrial ecosystems.

However, in some studies, the vegetation productivity in urban areas was set as zero [44]; the role of urban vegetation might be underestimated in the global carbon budget. With the improvement of urban green coverage (e.g., the green coverage in the urban built areas in Jiangsu province increased from 19.5% in 1990 to 42.1% in 2020), more attention should be paid to the effect of the urban forest in the global and regional carbon assimilation [45–47]. Although the ability of carbon assimilation can be estimated through GPP, the contribution of the urban forest to regional carbon assimilation is still unclear. A number of difficulties in the current estimations of urban GPP need to be overcome; for example the lack of surface observation data in urban areas makes it difficult to conduct empirical estimation on the urban scale by using traditional VI. At the same time, the lack of driving and parameterized data with high spatiotemporal resolution significantly reduced the performance of process-based ecological models in the simulation of GPP in urban areas.

As NIRv has successfully estimated GPP in other ecosystems, it may be able to provide a new approach to estimating carbon assimilation in urban forests. To the best of our knowledge, the contribution of urban forests to regional carbon assimilation has not yet been studied according to this method. Meanwhile, it is not yet clear which resolution of satellite data is appropriate for determining urban forest carbon assimilation capacity. Therefore, in this study, the GPP-NIRv relationship model is employed to estimate the GPP of terrestrial ecosystems in Suzhou by MODIS, Landsat-8, and Sentinel-2 remote sensing data, and to further explore what kind of remote images can figure out the spatiotemporal pattern of GPP of urban forests.  The expected results will improve our understanding of the ability of urban forests to affect carbon assimilation and will provide an appropriate reference for selecting the right resolution of remote sensing images for carbon budget research in urban areas.

#### **2. Materials and Methods**

#### *2.1. Study Area*

Suzhou is one of the most economically developed cities in eastern China, and is located in the southeast of Jiangsu and the middle of the Yangtze River Delta, with the region ranging from eastern longitude 119◦55 to 121◦20 and northern latitude from 30◦47 to 32◦02 . Suzhou connects to Shanghai to the east, Jiaxing and Huzhou of Zhejiang province to the north, Wuxi to the west, and the Yangtze River to the north, with a total area of 8657.32 km2. As one of the largest industrial cities in China, Suzhou has an urban population of 12.84 million, and the gross output value of all the above designated-size industrial enterprises in Suzhou exceeded CNY 4 trillion in 2021 (http://tjj.suzhou.gov.cn/ sztjj/tjnj/2021/zk/indexce.htm (accessed on 1 May 2022)).

Suzhou lies in the subtropical monsoon climate zone, where abundant precipitation and warm temperatures are suitable for vegetation growth. The annual precipitation and mean temperature in the year 2021 were 1318.6 mm and 18.3 ◦C, respectively. Suzhou is low and even, with a general elevation range from 3.5 to 5 m above sea level. The southeast of Suzhou is lower, with the lowest elevation below 2 m, and the southwest is a hilly area where vegetation grows well. The terrain of this area is shown in Figure 1. Suzhou is a famous water country region with a dense river network and numerous lakes, and rivers, lakes, tidal flats and wetlands account for 34.6% of the total area. Cropland is the dominant land use type in Suzhou, with an area of about 2871.6 km2, accounting for 33.17% of the total area.

**Figure 1.** The geographic location (**a**) and land cover (**b**) derived from Sentinel-2 by ESRI of Suzhou in Jiangsu province in China.

#### *2.2. Data Sources and Processing*

Remote sensing data, climate data and land cover data were used in this research to estimate GPP in Suzhou and analyze the importance of urban forests for the regional carbon budget. Three kinds of resolution remote sensing images are used to estimate GPP. Of these, 30 m-Landsat 8 and 10 m-Sentinel 2 surface reflectance datasets were downloaded from the United States Geological Survey (https://earthexplorer.usgs.gov/ (accessed on 10 April 2022)) and the Copernicus Open Access Hub (https://scihub.copernicus.eu/dhus/#/home (accessed on 15 April 2022)) in 2021 to calculate NIRv. In order to avoid the impact of cloud on data quality, the maximum value composite (MVC) method was set to generate the monthly NIRv.

MODIS 500 m surface reflectance data (MOD09A1) and NDVI data (MOD13A1) downloaded from NASA's Distributed Active Archive Center (DAAC) have eliminated the

impact of cloud. In order to maintain the consistency of the production of NIRv, the MVC method was also used to generate the MODIS monthly NIRv.

Climate data, including monthly temperature and precipitation, were used to investigate GPP changes in response to the local climate; climate data were obtained from the National Meteorological Science Data Center (https://data.cma.cn/ (accessed on 2 May 2022)).

MODIS, Landsat and Sentinel-2 land cover datasets were downloaded to summarize the GPP of different land cover types according to the corresponding estimation of GPP with the same spatial resolution. A MODIS land cover data (MCD12Q1) product was downloaded from NASA DAAC. The MCD12Q1 product supplies maps of land cover at annual steps, and the land cover in Suzhou was classified according to the International Geosphere Biosphere Programme (IGBP) land cover classification scheme as urban and built-up lands, croplands, grasslands, water bodies, evergreen broadleaf forests (EBF), deciduous broadleaf forest (DBF), mixed forest (MF) and so on.

Global land cover data producing 30 m resolution were downloaded from http://data. ess.tsinghua.edu.cn/ (accessed on 17 April 2022). This product used the amount of training samples across the world to optimize many kinds of classifiers, eg. maximum likelihood, decision tree, random forest, etc. A unique land-cover classification system was used in this product. In Suzhou, the typical land cover types were crop, forest, grass, shrub, water, and impervious. Additionally, impervious was considered as urban built-up areas in this study.

ESRI generated a global map of land use and land cover (https://livingatlas.arcgis. com/landcover/ (accessed on 17 April 2022)) derived from Sentinel-2 imagery at 10 m resolution by using a deep learning AI land classification model trained by billions of human-labeled image pixels [48]. This product has 9 classes, of which water, tree, crop, shrub land, built area and grass are the dominant ones in Suzhou.

In order to better calculate the changes in carbon assimilation of urban forest, the urban area of Suzhou was divided into built area (black areas in Figure 1b) and non-built area (remaining areas besides built area).

#### *2.3. Estimation of GPP*

NIRv has been found to accurately capture both the seasonal and annual variation in GPP at flux sites [25,26]. GPP correlates linearly with NIRv among different vegetation types (Table S1), and global GPP was estimated with high accuracy on a monthly basis by upscaling the relationships between NIRv and GPP [23]. In this study, we also use these correlations to determine the Suzhou monthly GPP based on NIRv, and the calculation of NIRv and GPP are shown in Equations (1) and (2):

$$NIRv = \frac{NIR - R}{NIR + R} \times NIR \tag{1}$$

$$GPP = a \times NIRv + b \tag{2}$$

*R* and *NIR* are the red and the near-infrared bands, *a* and *b* are derived from linear regression, and the values can be found in Table S1. We used the MVC method to generate the monthly NIRv. Landsat-8 has a limited frequency of revisits, so some images were lost during the rainy season, i.e., in July and August. MODIS and Sentinel-2 had sufficient data to produce a full image each month.

For different kinds of forests, such as deciduous broadleaf forests (DBF), evergreen broadleaf forests (EBF), evergreen needle forests (ENF), and mixed forests (MF), Wang et al. [23] used different coefficients for calculating GPP. However, the land cover data of Sentinel-2 used in this study do not distinguish forest type. As a result, based on the studies of [49], changes in tree NDVI between summer and winter were examined to determine the forest types, including deciduous broadleaf forest, mixed forest and evergreen broadleaf forest, according to the criteria listed in Criteria 2 of Table 1. Additionally, we did not include evergreen needle forest in calculating GPP, for the reason that its distribution in

Suzhou is limited and hard to separate from evergreen broadleaf forests by only the change in NDVI. Considering that the grassland and different kinds of forests in urban area are always integrated but still able to be identified according to dominated coverage of plants, vegetation types of urban forest in built area were also identified as grass (GRA), DBF, MF and EBF.

**Table 1.** NDVI for identifying vegetation type in land cover of tree and built area.

Note: The subscript is month, NDVI8 and NDVI12 are monthly NDVI in August and December. C1 and C2 are criteria 1 and criteria 2, respectively. Veg is the vegetation type. TR and BA are the land cover of tree and built area. BRE and GRA are bare land and grass land, respectively.

In the built area, the contribution of urban forests to total regional carbon assimilation cannot be ignored as they are widely distributed in Suzhou. However, they are usually buried in lower-resolution imagery. Considering the fragmentation of urban forests, three different resolutions of remote images, i.e., MODIS, Landsat and Sentinel-2, were used to evaluate the ability to detect urban forests. For calculating GPP in urban forest, NDVI changes in the built area were used to determine vegetation types as shown in Criterions 1 and 2 of Table 1.

#### *2.4. Data Analysis*

Based on our estimation, the GPPs of different land cover types were summarized according to land cover data. Then, we identified the contribution of urban forests to the regional total GPP. Vegetation in urban areas is inevitably affected by human activities. Therefore, the monthly changes of GPP were examined to find out the abnormal fluctuation of monthly GPP and identify what kind of human activities disturbed the change in GPP.

#### **3. Results**

#### *3.1. Variance of NIRv in Different Remote Sensing Data with Distinct Spatial Resolutions*

The NIRv can be calculated based on MODIS, Landsat-8 and Sentinel-2 by using Equation (1). Figure 2 shows the spatial distribution of Sentinel-2 NIRv monthly.

The area along the Yangtze River in the northern part of Suzhou was a high-value area of NIRV. Additionally, the changes in NIRv in this area were significantly greater than in other parts of the city. During the growing season, in April and from July to September, NIRv reached more than 0.45. The high NIRv value in this area was a result of the high photosynthetic ability of crops planted in this region, for example, wheat and rice. Furthermore, around Tai Lake, the large water body in the southwest of the city, there was another area with a high NIRv. Forests and crops were spread throughout this area (Figure 1b), and the peak NIRv was around 0.4 from June to August.

The NIRv in Suzhou's built area ranged from 0 to 0.1, generally, without an evidently monthly dynamic. Nevertheless, there were many points scattered in the built area with high NIRvs (ranging from 0.2 to 0.35) and also seasonal changes. The urban forest in built areas contributed to the scattered high-value points. This kind of vegetation is usually found in small urban green spaces, and was represented as discrete points by remote sensing data with high resolution.

**Figure 2.** Spatial distribution of monthly Sentinel-2 NIRv in 2021, Suzhou.

For comparing the differences of NIRv derived from MODIS, Landsat-8 and Sentinel-2, the monthly averages and standard deviations of NIRv were listed in Figure 3. Sentinel-2 NIRv showed two peaks in April and July, with values of 0.13 and 0.18, respectively. The values of MODIS NIRv were generally lower than Sentinel-2 and had similar trends to

Sentinel-2. Due to the impact of the cloud, the continuous change of Landsat-8 NIRv was hard to obtain. Additionally, current valid data of Landsat-8 were similar to MODIS NRIv and also showed a peak in September.

**Figure 3.** Monthly averages (red, blue and green cycle) and standard deviations (red, blue and green bar) of NIRv in Suzhou derived from MODIS, Landsat-8 and Sentinel-2.

#### *3.2. Comparing the GPPs Estimated by Different Resolutions of Remote Sensing Data*

The GPPs estimated by NIRv derived from MODIS, Landsat-8 and Sentinel-2 were shown in Table 2. The total GPP in Suzhou in 2021 estimated by MODIS, Landsat-8 and Sentinel-2 was 8.43, 8.48, and 9.30 Tg C yr−1 (Tg = 1012 g), respectively. Our estimations of GPP were higher than the results of MOD17A2 (4.37 Tg C yr−1), which ignores the contribution of the urban forest in built areas (black areas in Figure 4a). In contrast, our results included this ignoring component and fell within the range of previous studies which considered the contribution of urban forest (Table 2).

**Table 2.** Comparison of urban contribution to GPP in different GPP products.

**Figure 4.** The spatial distribution of GPP from MOD17A2 product (**a**) and estimated by MODIS (**b**), Landsat-8 (**c**), Sentinel-2 (**d**). Areas in panel (**a**) with no data are due to the MODIS product not calculating GPP if the land cover type is built area. Additionally, the no vegetation (veg.) areas in panels (**b**–**d**) are the regions where the NDVI is less than 0.2 and, as a result, is recognized as bare land with no GPP.

The spatial patterns of GPP in Suzhou estimated by MODIS, Landsat-8 and Sentinel-2 in 2021 were shown in Figure 4b–d. All three kinds of GPP showed the same high GPP region located around Tai Lake (grey area in the southwest of Figure 4). In many areas of this region, the GPP ranged from 2500 to 3000 g C m−2 yr−1. The estimated GPPs had some differences in the northeast, along the Yangtze River (grey area in the northeast of Figure 4). In this region, the GPP of MODIS was lower than the GPP of Landsat-8 and Sentinel-2. According to our MODIS GPP results, most built areas were capable of photosynthesis from the vegetation within them. Additionally, GPP in built area of Suzhou ranged from 250 to 500 g C m−2 yr−1. Ignoring the contribution of GPP from the built area (black area in Figure 4a) induced the lower GPP in MOD17A2. Furthermore, high-resolution results from Landsat-8 and Sentinel-2 also indicated the contribution of urban forests, which were scattered over the built area due to their spatial distribution being highly discrete in Suzhou (Figure 4c,d).

#### *3.3. Monthly Change of GPP in the Year 2021*

The Landsat-8 GPP failed to overcome the effects of clouds every month since its long revisit. Additionally, with the help of the MVC method, monthly NIRvs derived by MODIS and Sentinel-2 were generated. As a result, Figure 5 shows only the MODIS and Sentinel-2 GPPs, as well as the changes in monthly temperature and precipitation for the year 2021.

**Figure 5.** The monthly change of Suzhou's GPP in the year 2021.

The monthly GPP of MODIS and Sentinel-2 showed two peaks in 2021. One peak was in April, with GPPs 0.85 and 0.95 Tg C m−1 for MODIS and Sentinel-2, respectively. The other peak was in July, when the GPPs of MODIS and Sentinel-2 reached their summit for the whole year: 1.43 and 1.12 Tg C m−1, respectively. Between the two peaks was an evident decrement in GPP in May. GPP decreased by 16.4% and 25.9% according to the estimations of MODIS and Sentinel-2, respectively.

In the growing season, the MODIS GPP was generally lower than Sentinel-2 GPP. For example, from July to September, the average of Sentinel-2 GPP was 1.31 Tg C m−1, which was 21.3% higher than MODIS GPP in the same period. Additionally, a similar difference was also found in the growing peak in April. Moreover, in winter and early spring MODIS GPP and Sentinel-2 GPP were similar.

As Figure 6 showed, a significant (*p* < 0.001) linear relationship was found between the air temperature and GPP. R2 for the temperature–GPP relationships were 0.81 and 0.89 for Sentinel-2 and MODIS GPP, respectively (Figure 6a,b), which was higher than the precipitation–GPP relationship (Figure 6c,d). Therefore, the change in GPP was mainly determined by the air temperature in the year 2021. Additionally, the positive coefficient of slopes for the linear relationships indicated that temperature and precipitation had positive effects on GPP. However, the GPP declined when temperature and precipitation increased in May. This decrement in GPP implies that, besides climatic factors, some other factors determined the change in GPP in May. In urban areas, the impact of anthropogenic activities on vegetation cannot be ignored.

**Figure 6.** The relationship between the estimated GPP and climate factors. (**a**) and (**b**) are the correlations between temperature and monthly GPP for Sentinel2 and MODIS, respectively. (**c**) and (**d**) are the correlations between precipitation and monthly GPP for Sentinel2 and MODIS, respectively.

#### *3.4. Changes in GPP of Different Land Cover Types*

To better understand the importance of urban forests to the regional GPP, the Sentinel-2 and MODIS GPP changes in different land cover types were shown in Figure 7. According to the result of Sentinel-2, the average monthly GPP of vegetation in the built area (e.g., urban forest) and non-built area (e.g., crop, tree, grassland) was 0.35 Tg C and 0.43 Tg C (0.32 for crops, 0.08 for trees, and 0.03 Tg C for grassland), respectively. Additionally, for MODIS, GPP of the built area and non-built area was 0.28 Tg C and 0.42 Tg C (0.18, 0.01, and 0.23 Tg C for crop, tree, and grassland, respectively). The MODIS GPP of grassland was 0.20 Tg C higher than Sentinel-2 GPP, and the MODIS GPP of cropland was 0.14 Tg C lower than Sentinel-2 GPP. This may be caused by the use of different land cover datasets in statistics. MODIS land cover products might be misplaced between farmland and grassland due to the lower spatial resolution. In fact, as an industrialized city, Suzhou does not have a lot of grasslands.

The GPP of each land cover type varied throughout the year (Figure 7a,c). Similar to the change in total GPP, each land cover type had two peaks in spring and summer. Additionally, they also showed decrements in GPP in May. These decrements in GPP for crops and urban forests were obviously about 32.33% and 23.57%, respectively, according to Sentinel-2 GPP.

Sentinel-2 results indicated that crops and urban forests accounted for a majority of the total GPP in Suzhou in 2021, with an average of about 38.15% and 44.43%, respectively. Trees and grasses in non-built area made relatively small contributions to total GPP, accounting for approximately 11.87% and 4.17%, respectively. According to MODIS GPP, urban forests still accounted for the majority (39.66%) of total GPP, followed by grassland and crops at 33.66% and 25.60%, respectively.

**Figure 7.** The Sentinel-2 and MODIS GPP changes in different land cover types in 2021 (**a**,**c**) and their contribution to total GPP (**b**,**d**).

The contribution of different land cover types also varied throughout the year (Figure 7b,d). Taking the most contributed urban forests as an example, their Sentinel-2 GPP exceeded 50% of the total GPP for five months of the year, with the highest percentage at 55.39% in June. Additionally, the contribution made by crops was also large, even more than the contribution of urban forests in growing peak seasons, such as April and August.

#### *3.5. Changes in GPP by Anthropogenic Factors*

The spatial distribution of this decrement in GPP in May is shown in Figure 8. The evident decrement of GPP for the crop was mainly in the rice and wheat rotation area along the Yangtze River. In May, as wheat matures, its photosynthetic capacity also decreases significantly. Additionally, in June, after the wheat harvest, rice has just been planted. In consequence, the photosynthetic capacity had not yet recovered.

**Figure 8.** The difference in GPP between April and May in Suzhou.

Furthermore, urban forests' GPP also decreased significantly in May, but recovered quickly in June (Figure 7a,c). For making urban vegetation more attractive, urban forests in Suzhou undergo a pruning process in May. This process removes a considerable amount of leaves from the canopy. Usually, these leaves are current-year leaves and are located at the top of the canopy. Additionally, the photosynthetic capacity of these kinds of leaves is greater than that of old leaves and lower canopy leaves [52,53]. The pruning process is usually conducted on the grassland and shrub land, which were the main parts of urban green spaces in Suzhou, approximately 55.67% (Figure 9) according to our criteria in Table 1. Haberl [54] proposed that the averaged biomass loss during gardening (such as pruning) or park and infrastructure maintenance amounted to 50% of the aboveground Net Primary Productivity (the remainder of GPP deducts autotrophic respiration). Therefore, the pruning process removes the most photosynthetic part of the canopy, which results in the decline of GPP of urban forests.

#### **4. Discussion**

#### *4.1. Uncertainty of Estimating GPP by NIRv*

Although, as we know, few studies have examined vegetation productivity in Suzhou, studies in similar cities may provide useful information. From 2000 to 2014, the average rate of GPP in Shanghai, the city next to Suzhou, was between 800 and 1050 g C m−2 yr−1 based on the simulations by Vegetation Photosynthesis Model (VPM) [55]. This result is close to our 973.74, 979.51, and 1074.21 g C m−2 yr−1 for MODIS, Landsat-8, and Sentinel-2, respectively.

According to the MOD17A2 product, the previous study also indicated the GPP of crops along the Yangtze River in the northeast of Suzhou as around 1200–1400 g C m−2 yr−1 [56], which is lower than our Sentinel-2 GPP, 1250–2000gCm−2 yr−1, in the same region (Figure 4d). Limited by mixed pixel effects, the MOD17A2 product usually overestimates in lower value and underestimates in higher value [57,58]. The high resolution of Sentinel-2 eliminated the many effects of the mixed pixel. Therefore, the Sentinel-2 GPP was obviously higher than the GPP calculated by MODIS reflectance. Furthermore, Landsat-8 GPP would have to be greater than MODIS GPP if it did not miss images during July and August, when carbon assimilation rates are at their peak.

Compared to deciduous forests, evergreen forests are inflexible to short-term changes in environmental conditions. As a result, their NIRv is suitable for predicting GPP over a longer period, such as 90 days [25]. Consequently, the NIRv–GPP relationship is a challenge in estimating monthly GPP for evergreen forests. Additionally, for this reason, our GPP of evergreen forests in the southwest of Suzhou was more than 2500 g C m−2 yr−1, which is higher than the previously observed and simulated results [59,60]. Additionally, the original NIRv-GPP relationships across the different land cover types are generated by 0.05 degree AVHRR reflectance [23]. The uncertainty induced by the different remote sensing data in terms of field of view angle, signal-to-noise ratio, and spectral width needs to be further evaluated.

#### *4.2. Evaluation of High-Resolution Remote Sensing Images in Urban Carbon Research*