Assessment of Forest Ecological Function Levels Based on Multi-Source Data and Machine Learning

Abstract

Forest ecological function is one of the key indicators reflecting the quality of forest resources. The traditional weighting method to assess forest ecological function is based on a large amount of ground survey data; it is accurate but costly and time-consuming. This study utilized three machine learning algorithms to estimate forest ecological function levels based on multi-source data, including Sentinel-2 optical remote sensing images and digital elevation model (DEM) and forest resource planning and design survey data. The experimental results showed that Random Forest (RF) was the optimal model, with overall accuracy of 0.82, recall of 0.66, and F1 of 0.62, followed by CatBoost (overall accuracy = 0.82, recall = 0.62, F1 = 0.58) and LightGBM (overall accuracy = 0.76, recall = 0.61, F1 = 0.58). Except for the indicators from remote sensing images and DEM data, the five ground survey indicators of forest origin (QI_YUAN), tree age group (LING_ZU), forest category (LIN_ZHONG), dominant species (YOU_SHI_SZ), and tree age (NL) were used in the modeling and prediction. Compared to the traditional methods, the proposed algorithm has lower cost and stronger timeliness.

1. Introduction

As one of the most important components of the ecosystem, forests provide the basic natural resource foundation for the sustainable development of human beings [1]. The forest ecological function index can comprehensively reflect the structure and ecological benefits of forests. Therefore, establishing a scientific and dynamic comprehensive system to evaluate forest ecological functions plays an important role in accurately addressing ecological and economic development [2]. As research progresses, a single analysis of a particular forest characteristic no longer meets the current requirements. Therefore, it is important to assess forests’ ecological functions by integrating the synergistic effects of multiple factors [3,4]. Many scholars have attempted to present concepts or tools for forest ecological function analysis [5,6,7], greatly helping to measure and quantify species’ actions.

Currently, some progress in the quantitative estimation of forest ecological functions has been made in China. In terms of research methods, Li Ma et al. established the Beijing Forest Ecosystem Health Evaluation Index System (EIS-BFEH) to evaluate the health function of forest ecosystems and used the hierarchical analysis process (AHP) to obtain a comprehensive index (CI) representing the health status of forest ecosystems [8]. Fang Xiaomin et al. used a comprehensive index method and a statistical grouping method to evaluate forest functions and calculated the forest ecological function index to compare forest ecological functions for different age groups, species structure, and origin [9]. Du Qun et al. proposed that the three factors of forest quantity, quality, and spatial distribution should be used for forest ecological function evaluation [10]. Although there has been some progress in research on computational methods and evaluation factors, the current studies are mainly focused on simple statistical methods, and there is a lack of attempts of using machine learning methods; this procedure cannot meet the requirements of accurately evaluating and monitoring forest ecological functions based on multi-data [11].

Among machine learning algorithms, ensemble learning algorithms can combine multiple classifiers together to improve the accuracy and generalization ability of classifiers. Overall, based on the presence or absence of dependencies between base classifiers, ensemble learning algorithms are divided into two types: boosting and bagging. Boosting algorithms have a strong dependency between base classifiers, and a series of base classifiers needs to be generated serially, represented by AdaBoost and GBDT. In fact, GBDT is more suitable for multi-category classification, and LightGBM and CatBoost are two important improved algorithms based on GBDT. In contrast, bagging algorithms do not have a strong dependency between base classifiers, and a series of base classifiers can be generated in parallel, represented by Random Forest. Regarding the research data, Wang Daling et al. evaluated the forest ecological function of arboreal forests based on the data of subcompartments from the 2018 forest resources planning and design survey of the Sanchazi Forestry Bureau [12]. Jun Yang et al. conducted a qualitative analysis of the health and ecological status of Chinese forests in 2009 based on the data of the seventh national forest resource inventory [13]. Liu Lixia et al. evaluated the ecological functions of forests based on the data of the forest resources planning and design survey [14]. Zhang Xianwu et al. used the results of the three continuous forest inventories from 2004, 2009, and 2014 to analyze the current situation, changes, and reasons of changes in the comprehensive index of forest ecological functions in Shanghai using nine indicators including forest stock, forest naturalness, and proportion of forest land area to national land area [15]. Even though forest resource inventory data has high accuracy, there are still issues of high cost and difficulty in data acquisition. Forest managers, decision makers, and politicians need to be able to make data-driven rapid decisions based on short-term and long-term monitoring information, complex modeling, and analysis approaches [11]. Thus, researchers are increasingly considering incorporating lower-cost data such as remote-sensing images and lower-cost strategies into our study.

It is necessary to attempt to integrate multi-source data to evaluate the forest ecological function. In recent years, the application of multi-source data fusion in the forestry field has become increasingly common, indicating that multi-source data have broad application prospects. For instance, Wang et al. used multi-source remote sensing data (Gaofen 1, Sentinel-2, Landsat 9, and Gaofen 3) to classify mangrove species in urban areas of Leizhou City, Guangdong Province [16]; Abd Rahman Kassim et al. used hyperspectral images and airborne LiDAR data to evaluate the ecological status of the FRIM campus forest ecosystem [17]. Some of the studies on multi-source data fusion in the field of forestry have achieved good results, but there is still more room to explore its application in forest ecological function level assessment.

This paper aimed to combine multi-source data with machine learning algorithms to assess the forest ecological function levels based on the unit of forest subcompartment. Spectral features and topographic features of the study area were collected from remotely sensed images and DEM. Some ground data with low acquisition cost, such as forest origin, tree age group, forest category, dominant species, and tree age, were provided by the forest resource planning and design survey. Based on the multi-source data scheme, three classic machine learning algorithms, i.e., Random Forest, LightGBM, and CatBoost, were involved in the study.

2. Materials and Methods

2.1. Schematic Framework of Materials and Methods

The schematic framework of the Materials and Methods of this study is shown in Figure 1, including multi-source data, preprocessing methods, models to predict forest ecological function levels.

2.2. Overview of the Study Area

Lin’an District (118°51′~119°52′ E, 29°56′~30°23′ N) is located in the northwest of Zhejiang Province, with a total area of 3126.8 square kilometers, shown in Figure 2. It is in the central subtropical monsoon climate zone, rich in plant resources. As a typical southern forest city, the forestry land area of Lin’an is up to 263,868.79 hm², with 1603.88 hm³ of forest standing stock and 81.99% of forest canopy density.

2.3. Processing of Label Dataset

According to the Technical Operation Rules for forest resources planning and design survey of Zhejiang Province of 2014 [18] and the Technical Regulations for the Continuous Inventory of Forest Resources of 2020 (GB/T 38590-2020) [19], there are eight investigation factors (

Table 1. Factors to evaluate forest ecological function levels.

Code	Factors	Classification Standards			Weight	References
		Ⅰ	Ⅱ	Ⅲ
1	Forest biomass (t/hm2)	≥150	50~149	<50	0.20	[19]
2	Forest naturalness	1, 2	3, 4	5	0.15
3	Forest community structure	1	2	3	0.15
4	Tree species structure	6, 7	3, 4, 5	1, 2	0.15
5	Vegetation coverage (%)	≥70	50~69	<50	0.10
6	Canopy density	≥0.70	0.40~0.69	0.20~0.39	0.10
7	Mean tree height/m	≥15.0	5.0~14.9	<5.0	0.10
8	Thickness of dead leaves	1	2	3	0.05

), i.e., forest biomass, forest naturalness, forest community structure, tree species structure, vegetation coverage, forest canopy density, mean tree height, and thickness of dead leaves, involved into the traditional algorithm to evaluate a forest ecological function levels. Their weights were determined according to their importance to the evaluation of forest ecological functions (Table 1).

To facilitate the data standardization, the values of the above eight evaluation factors were uniformly classified into three classes, i.e., Ⅰ, Ⅱ, and Ⅲ. According to the data type and the distribution of the values, the above eight evaluation factors were divided into three categories. The first category, including vegetation coverage, canopy density, mean tree height, forest community structure, and thickness of dead leaves, was directly classified into types Ⅰ, Ⅱ, and Ⅲ (as shown in Table 1). The second category. consisting of forest naturalness and tree species structure, was classified first by the division standards of

Table 2. Criteria and codes for the classification of naturalness in the continuous inventory of forest resources.

Naturalness	Division Standard	Code	References
Ⅰ	Forest types are pristine or in a largely untouched state, with little human influence.	1	[19]
Ⅱ	Natural forest types with obvious human interference or secondary forest types in the later stage of succession, mainly consisting of tree species with high adaptability at the top level of zonality.	2
Ⅲ	A secondary forest type with great human disturbance, in the late stage of secondary succession. In addition to pioneer species, top-level species can also be seen.	3
Ⅳ	Highly disturbed by humans, succession retrograde, is in an extremely fragile secondary forest stage.	4
Ⅴ	Highly and continuously disturbed by humans, with the destruction of almost all zonal forest types, in the late stage of difficult-to-recover retrograde succession.	5

for naturalness and

Table 3. Criteria and codes for the classification of forest species structure.

Tree Species Structure Type	Division Standard	Code	References
Ⅰ	Pure coniferous forests, where the volume of individual coniferous species is greater than or equal to 90% of the total volume.	1	[19]
Ⅱ	Pure broadleaved forests, where the volume of individual broadleaved species is greater than or equal to 90% of the total volume.	2
Ⅲ	Relatively pure coniferous forest, where the volume of individual coniferous species is greater than or equal to 65% and less than 90% of the total volume.	3
Ⅳ	Relatively pure broad-leaved forests, where the volume of individual broad-leaved species is greater than or equal to 65% and less than 90% of the total volume.	4
Ⅴ	Mixed coniferous forests, where the volume of total coniferous species is greater than or equal to 65% of the total volume.	5
Ⅵ	Mixed coniferous and broad-leaved forests, where the volume of total coniferous species or total broad-leaved species is greater than or equal to 35% and less than 65% of the total volume.	6
Ⅶ	Broad-leaved mixed forests, where the volume of total broad-leaved species is greater than or equal to 65% of the total volume.	7

for tree species structure, and then into types of Ⅰ, Ⅱ, and Ⅲ. The third category merely included forest biomass, which was firstly calculated according to the amount of forest volume for various dominant species by equations (shown in

Table 4. Biomass models of major tree species and vegetation types.

Code	Tree Species/Vegetation Type	Biomass Model	References
1	Cunninghamia lanceolata	W = 0.3999 V + 22.5410	[20]
2	P. massoniana	W = 0.5101 V + 1.0451
3	Other pine and conifer tree species (besides P. massoniana, Tsuga, Cryptomeria, and Keteleeria), coniferous mixed forest	W = 0.5168 V + 33.2378
4	Cypress	W = 0.6129 V + 46.1451
5	Mixed conifer and deciduous forests	W = 0.8019 V + 12.2799
6	Betula	W = 0.9644 V + 0.8485
7	Deciduous oaks	W = 1.3288 V – 3.8999
8	Eucalyptus	W = 1.0357 V + 8.0591
9	Mixed deciduous and Sassafras	W = 0.6255 V + 91.0013
10	Tsuga, Cryptomeria, Keteleeria	W = 0.4158 V + 41.3318

Note: W is the biomass of forest stand measured in t/hm², V is the forest volume per hectare measured in m³/hm².

) and then classified into types Ⅰ, Ⅱ, and Ⅲ.

2.4. Data Sources and Pre-Processing

2.4.1. Data Sources

The data sources mainly included remote sensing images and Digital Elevation Model and ground survey data. The remote sensing images from the satellite Sentinel-2 (with 13 Bands, a spatial resolution of 10 m, 20 m, and 60 m) and DEM (with a spatial resolution of 30 m) with a format of ASTER GDEM, were all downloaded from the International Science and Technology Data Mirror of the Computer Network Information Centre of the Chinese Academy of Sciences (www.gscloud.cn) on 27 September 2021. The ground data were obtained from the forest resources planning and design survey provided by the Lin’an District Forestry Bureau in 2019.

2.4.2. Data Pre-Processing

(1)

Forest Resources Planning and Design Survey Data

The original data of forest resource planning and design survey consist of 119,792 sub-compartments. To eliminate erroneous data, two steps were taken. Firstly, the data with non-forest, zero volume, and null value of forest ecological function level were removed. Secondly, according to the Pauta criterion [21], the abnormal data exceeding the mean (u) ± three times the standard deviation (3σ) were also removed. Consequently, 47,596 valid samples were retained (

Table 5. Number of experimental samples.

Number of Original Samples	Number of Valid Samples	Number of Tree Species
119,792	47,596	26

), consisting of 26 dominant tree species, that is, broad-leaved mixed forest, horsetail pine, fir, coniferous mixed forest, coniferous mixed forest, oak, other hard broad-leaved forest, maple, camphor, yellow pine, etc.

(2)

Extraction and processing of characteristic factors based on images from remote sensing and DEM

For the Sentinel-2 remote sensing images, Sen2Cor^® (v2.8, European Space Agency, Paris, France) was used for atmospheric correction, and SNAP^® (v6.0, European Space Agency, Paris, France) was exploited to resample the bands (band1, band5, band6, band7, band8A, band9, band10, band11, and band12) to fuse the lower-resolution images of 20 m and 60 m with higher-resolution images of 10 m by the nearest neighbor method. After that, by the operations of image mosaic and clipping in ArcGIS^® (v10.8, Environmental Systems Research Institute, Inc., Redlands, CA, USA), a valid and complete remote sensing image of the Lin’an District (Figure 3) was produced.

The characteristic factors extracted from the Sentinel-2 optical remote sensing images consisted of two main types: original factors and derived factors. The former consisted of 13 bands [22,23,24]. There were three bands, namely, Band1 for the coastal band, Band9 for the water vapor band, and Band10 for the cirrus band, which were not relevant to the experiment and were removed. Therefore, the remaining ten original bands (as shown in

Table 6. Vegetation index formula.

Code	Vegetation Index	Formula
1	Atmospherically resistant vegetation index (ARVI)	ARVI = (NIR – (2 * R) + B)/(NIR + (2 * R) + B)
2	Enhanced vegetation index (EVI)	EVI = 2.5 × (NIR − R)/(NIR + 6 × R − 7.5 × B + 1)
3	Differential environmental vegetation index (DVI)	DVI = NIR − R
4	Normalized vegetation index (NDVI)	NDVI = (NIR − R)/(NIR + R)
5	Ratio red-edge vegetation index (RVIre)	RVIre = NIR/Re
6	Inverted red-edge chlorophyll index (IRECI)	IRECI = (Re3 − R)/(Re1 − Re2)
7	Normalized red-edge vegetation index1 (NDVIre1)	NDVIre1 = (NIR − Re1)/(NIR + Re1)
8	Normalized red-edge vegetation index2 (NDVIre2)	NDVIre2 = (NIR − Re2)/(NIR + Re2)
9	Non-linear red-edge index (NLIre)	NLIre = ((NIR * NIR) − Re1)/((NIR * NIR) + Re1)
10	Improved normalized red-edge vegetation index (mNDVIre)	mNDVIre = (NIR − Re1)/(NIR + Re1 − 2 * B)
11	Red-edge chlorophyll index (CIre)	CIre = (NIR/Re1) − 1

Note: R represents the red band, B represents the blue band, NIR represents the near-infrared band, and Re represents the red-edge band.

) and eleven vegetation indices-derived factors (as shown in

Table 7. Evaluation factors.

No.	Factor Name	Explanation	Source of Data
1	Band 2	Bule	Sentinel-2
2	Band 3	Green
3	Band 4	Red
4	Band 5	VNIR1
5	Band 6	VNIR2
6	Band 7	VNIR3
7	Band 8	NIR
8	Band 8A	Narrow NIR
9	Band 11	SWIR 1
10	Band 12	SWIR 2
11	HAI_BA	Elevation	DEM
12	PO_DU	Slope
13	PO_XIANG	Aspect
14	LIN_ZHONG	Forest category	Forest Resources Planning and Design Survey Data
15	QI_YUAN	Forest origin
16	YOU_SHI_SZ	Dominant species
17	NL	Tree age
18	LING_ZU	Tree age group
19–29	Refer to Table 6		Vegetation indices generated from optical remote sensing images

), that is, a total of 21 spectral feature factors, were used as independent variables.

In this study, the remote sensing images and DEM were preprocessed by converting, stitching, and cropping in ArcGIS. The three topographic factors—elevation, slope, and aspect—were obtained from DEM images (Figure 4), with a spatial resolution of 30 m.

2.5. Extraction of Feature Factors from Ground Survey Data

Even if the forest Resources Planning and Design Survey data have authenticity and reliability, many factors, such as forest biomass, vegetation coverage, forest canopy density, and mean tree height, with very high investigation costs, are not suitable for evaluating forest ecological function levels. Therefore, the eight factors (shown in Table 1) which were traditionally used to compute forest ecological function levels were removed, and unrelated factors such as county code and county name were also removed. After that, 36 variables remained. To further reduce the dimensionality of the features and improve the efficiency of the model, the feature importance was ranked (as shown in Figure 5). Furthermore, five feature factors, that is, QI_YUAN, LING_ZU, NL, YOU_SHI_SZ, and LIN_ZHON, were selected as independent variables of the model from top to bottom.

2.6. Multi-Source Data Integration

A total of 29 evaluation factors (as shown in Table 7) were involved in the study, including 10 spectral bands, 11 vegetation indices, 5 ground survey factors, and 3 topographic factors. The integration of the multi-source data was implemented according to the FID.

Finally, the dataset was randomly divided into a training set with 80% of the samples for modeling and a test set with 20% of the samples for testing.

2.7. Methods

2.7.1. Grid SearchCV

To prevent overfitting and underfitting, a hyper-parameter optimization method—grid search (Grid SearchCV [25])—was used to select the optimal hyper-parameter values for the three models. Grid SearchCV allows performing hyperparameter tuning in order to determine the optimal values for a given model. Specifically, based on a specified parameter range and a validation dataset, the parameters are gradually adjusted based on a pre-set step size, and finally the optimal parameter value was selected with the highest accuracy.

2.7.2. Random Forest (RF)

Random forest [26,27] is a collection of multiple decision tree algorithms with random sampling, which is a combination of Breiman’s “bagging” idea and a random selection of features. The procedure consists in making a precise prediction by taking the average or mode of the output of multiple decision trees (shown in Figure 6). Typically, the greater the decision trees’ number, the more precise the output and the greater the overhead. In our study, to balance accuracy and overhead, the default number of decision trees was set to 100, and the parameter range of the grid search was set to 100–500 with a step size of 5. The final experimental results showed that when the number of subtrees reached 200, the increase in the number of subtrees had a minimal effect on the model enhancement. Therefore, the optimal parameter value was finally set to 200. Also, the maximum number of features was set to 195, which was the square root of the number of training samples (47,596 * 0.8).

2.7.3. Light Gradient Boosting Machine (LightGBM)

The Light Gradient Boosting Machine, XGBoost, and Catboost are lifting algorithms [28]. For the LightGBM, during the training process, the decision tree algorithm of Histogram was adopted, which greatly reduced the calculation amount of the model. Meanwhile, the leaf-wise growth strategy (shown in Figure 7) was introduced into the growth process of the subtree, which reduced the splitting of invalid nodes.

There are three hyperparameters in LightGBM that need to be determined through the grid search method, namely, n_estimators, max_depth, and learning_rate. Here, n_estimators is the maximum number of base learners, max_depth is the maximum depth of the tree, and learning_rate indicates the magnitude of each parameter update. Correspondingly, the parameter ranges were set to [100, 300], [2, 10], and [0.05, 0.2], with step sizes of 5, 1, and 0.01, respectively. Finally, the optimal values of the three parameters were 200 for n_estimators, 3 for max_depth, and 0.1 for learning_rate, respectively.

2.7.4. CatBoost

CatBoost [29,30] is a GBDT algorithm based on a symmetric binary tree, which can automatically process category-based features and effectively solve gradient bias and prediction shift problems and has excellent accuracy and generalization capabilities. In addition, CatBoost combines multiple categorical features by adding a priori distributed modified Greedy TS approach to reduce the effect of noise and low-frequency categorical data on the data distribution (as calculated in Equation (1)).

$$X_{{\sigma }_p,k}=\frac{{\sum }_{p-1}^{j=1}\left[X_{{\sigma }_j,k}=X_{{\sigma }_p,k}\right]Y_{{\sigma }_j}+a\times p}{{\sum }_{p-1}^{j=1}\left[X_{{\sigma }_j,k}=X_{{\sigma }_p,k}\right]+a}$$

(1)

where $a$ is a weighting factor greater than 0; $p$ is the priori.

The default value of n_estimators in CatBoost package we used was 500, which was already large enough. For the other two parameters, tree depth and learning rate, their ranges were set to [5, 12] and [0.05, 0.2], with step sizes of 1 and 0.01, respectively. Consequently, the optimal tree depth was 11, and the learning rate was 0.05.

2.7.5. Performance Metrics

The performance metrics of the model are generally calculated based on a confusion matrix [31], as shown in

Table 8. Confusion matrix for multi-classification models.

Confusion Matrix		Predicted Value
		Category 1	Category 2	Category k	Total
Measured value	Category 1	$a_{11}$	$a_{21}$	$a_{1k}$	$a_{1+}$
	Category 2	$a_{21}$	$a_{22}$	$a_{1k}$	$a_{2+}$
	Category k	$a_{k1}$	$a_{k2}$	$a_{kk}$	$a_{3+}$
	Total	$a_{+1}$	$a_{+2}$	$a_{+3}$	N

. Here, $a_{ij}$ denotes the number of samples, the measured value is denoted by $i$, and the predicted value is denoted by $j$, $N$ is the total number of samples, $k$ is the number of target categories, and $a_{i+}=\sum _j{a}_{ij}$, $a_{+j}=\sum _i{a}_{ij}$.

Furthermore, the performance of the predicted results was evaluated by accuracy (Formulas (3)), recall (Formulas (5)), and F1_score (Formulas (8)).

$$accurac{y}_i=\frac{a_{ii}}{N}$$

(2)

$$accuracy={\displaystyle \sum }_{i=1}^{k}accurac{y}_i$$

(3)

$$recal{l}_i=\frac{a_{ii}}{a_{+i}}$$

(4)

$$recall={\displaystyle \sum }_{i=1}^{k}recal{l}_i$$

(5)

$$precisio{n}_i=\frac{a_{ii}}{a_{i+}}$$

(6)

$$F1\_scor{e}_i=2\ast \frac{precisio{n}_i\ast recal{l}_i}{precisio{n}_i+recal{l}_i}$$

(7)

$$F1\_score={\displaystyle \sum }_{i=1}^{k}F1\_scor{e}_i$$

(8)

3. Results

3.1. Labeling of the Data

The quantification scores for Type I, II, and III in Table 1 were assigned to 1, 2, and 3, respectively [32], and the composite score of forest ecological functions were calculated according to Formula (9) [19]

$$Y={\displaystyle \sum }_{i=1}^8{W}_i{X}_i$$

(9)

where $Y$ is the composite score of the forest, $X_i$ is the standardized score of the i-th evaluation factor, and $W_i$ is the weight of the i-th evaluation factor.

Hence, the ecological function index was calculated according to the composite score and is represented by Equation (10)

$$K=\frac{1}{Y}$$

(10)

where $K$ is the the ecological function index with a value less than or equal to 1. The larger the value of $K$, the better the forest’s ecological function.

Further, according to the values of $K$, the ecological function levels were divided into three groups, i.e., Good, Medium, and Poor (as shown in

Table 9. Criteria and codes for rating forest ecological functions.

Ecological Function Level	Comprehensive Score Value (Y)	Forest Ecological Function Index (K)	Code	References
Good	<1.5	>0.6667	1	[19]
Medium	1.5~2.4	0.6667~0.4167	2
Poor	≥2.5	≤0.4	3

).

Consequently, the classification results based on Equation (10) acted as the labeled data and were plotted in Figure 8. As shown, 10,844 forest subcompartments were graded as “good”, 36,365 forest subcompartments were graded as “medium”, and 386 forest subcompartments were graded as “poor”. In terms of area, the ecological function levels of “good”, “medium”, and “poor” were 59,094.3 hectares, 92,011.9 hectares, and 487.6 hectares, respectively, accounting for 38.98%, 60.70%, and 0.32% of the total area of forest land, respectively. The average forest ecological function index of forested land in Lin’an was 0.63, which was only 0.04 lower than the good ecological function index of 0.67, and the overall ecological function was at a medium, tending to good, level.

There was an upward tendency for forest ecological function levels from east to west, which was consistent with the zoning of construction land in Lin’an District, with more urban construction land to the east and more ecological forest areas to the west. The extent to which the experimental results were consistent with the actual situation from the cross-reference map of the focus areas (Figure 9) was determined. For example, the two national nature reserves in Figure 9a,b corresponded to areas where the ecological function of the forest was rated as “good” on a large scale and, partly, as “medium”, with no areas rated as “poor”. Closer to urban areas, there was a greater chance of a “poor” forest ecological function level (as shown in Figure 9c). To some extent, this reflected the impact of human activities on the levels of forest ecological function. The regions with lower human activity often had higher levels of forest ecological function, while the regions with higher human activity had lower levels of forest ecological function.

3.2. Design of the Data Scheme

Four data combination schemes (as shown in

Table 10. Data combination schemes.

Data Combination Scheme	Data Source
A	Sentinel-2
B	Sentinel-2, DEM
C	Sentinel-2, forest resource planning and design survey data
D	Sentinel-2, DEM, forest resource planning and design survey data

)—A, B, C, and D—were designed according to the three data sources.

3.3. Testing Results

The results of the optimal hyperparameters for three models were obtained by grid search (Section 2.7 for details), as shown in

Table 11. Combinations of optimal hyperparameters for three models.

Model	Optimal Values of Hyperparameters
RF	n_estimators = 200, max_features = 195
LightGBM	n_estimators = 200, max_depth = 3, learning_rate = 0.1
CatBoost	n_estimators = 500, depth = 11, learning_rate = 0.05

. Furthermore, the four different data combination schemes (shown in Table 10) were modeled and analyzed using the random forest, LightGBM, and CatBoost algorithms, and the experimental results are shown in

Table 12. Performance of the three models based on the four data combination schemes.

Program	Overall Accuracy Rate	Category Accuracy Rate			Recall	F1 Score
		Good	Medium	Poor
RF-A	0.46	0.57	0.80	0.35	0.39	0.39
RF-B	0.47	0.62	0.80	0.40	0.40	0.41
RF-C	0.82	0.73	0.89	0.80	0.54	0.57
RF-D	0.82	0.76	0.89	0.83	0.66	0.62
LightGBM-A	0.47	0.61	0.80	0.32	0.41	0.40
LightGBM-B	0.47	0.62	0.80	0.33	0.40	0.41
LightGBM-C	0.73	0.71	0.90	0.58	0.52	0.55
LightGBM-D	0.76	0.73	0.90	0.64	0.61	0.58
CatBoost-A	0.46	0.59	0.80	0.35	0.42	0.41
CatBoost-B	0.48	0.62	0.81	0.42	0.42	0.43
CatBoost-C	0.73	0.73	0.90	0.57	0.55	0.56
CatBoost-D	0.82	0.75	0.90	0.80	0.63	0.58

.

Comparing the performance metrics of the four data schemes in the three models in Table 12, the data scheme A (single-data combination scheme) had the worst performance, with an overall accuracy rate of only 0.46~0.47 and a classification accuracy of 0.32~0.35 for the “poor” category samples. However, the data scheme D (multi-source-data combination scheme) performed the best, with an overall accuracy rate of 0.76~0.82. The accuracy of the “good”, “medium”, and “poor” categories reached 0.73~0.76, 0.89~0.90, and 0.64~0.80, respectively, and the F1 score was 0.58~0.62. When comparing the data scheme B (after adding the DEM data to scheme A) with A, it was found that the DEM had an insignificant contribution to the model, with an overall accuracy improvement of only 0.01~0.02, and the accuracy of the “good”, “medium”, and “poor” categories, respectively, improved by 0.01~0.16, 0~0.01, and 0.01~0.07. However, when comparing the data scheme C (after adding the ground survey data to scheme A) with A, it was found that the addition of the forest resource planning and design survey data made a significant positive contribution to the model, with an overall accuracy improvement of 0.26~0.36, significantly improving the performance metrics of “good”, “medium”, and “poor” categories, respectively.

Ultimately, the RF-D executed the optimal program with overall accuracy of 0.82, recall of 0.66, and F1 score of 0.62, and the classification accuracy was significantly improved, especially for the small sample category of “poor”.

3.4. Ranking of Features’ Importance

Based on the optimal data scheme D, the performance metrics (as shown in Table 12) were calculated by the three machine learning algorithms of RF, LightGBM, and CatBoost, and the ranking of the feature importance was obtained and shown in Figure 10.

As shown in Figure 10, there were five factors from the ground survey data that led to a higher ranking of feature importance and played an important role in the model. Nevertheless, there were three factors from the DEM data with a lower ranking of feature importance that played an unimportant role in the model, which is consistent with the results in Table 11. In the optical remote sensing data, the factors of b12, NDVIre2, EVI, IRECI, b2, b11, and b4 ranked relatively high in the model.

4. Discussion

4.1. Performance Metrics

Even if many studies have been aimed to the forest ecological function rating methodology, some deficiencies still exist, as reported below.

(1)

The evaluation indicators are heavily influenced by foresters’ experience

In the evaluation of forest ecological functions, due to many qualitative indicators such as forest naturalness, tree species structure, and thickness of dead leaves, the original way of obtaining data is heavily influenced by foresters’ experience, which makes it difficult to establish an objective indicator system and may lead to inconsistent forest ecological function levels evaluated by different foresters for the same forest stand [33,34].

(2)

High cost of data acquisition for evaluation indicators

The data of the eight evaluation factors (forest naturalness, forest community structure, tree species structure, vegetation coverage, forest canopy density, mean tree height, and thickness of dead leaves) used to calculate the forest ecological function levels are obtained from ground surveys, which leads to high acquisition costs and time consumption [35].

The proposed evaluation method combined the advantages of multiple sources of data and machine learning algorithms to effectively reduce the human influence on the evaluation system and save data acquisition expenses and time. Specifically, the addition of remote sensing data effectively reduced the influence of human subjective experiences and increased the frequency of data acquisition. The costs of acquiring remote sensing data are also much lower than that of ground surveys [36], especially in areas that are difficult to access by humans, such as deep forests and cliffs.

4.2. Complementarity of Multi-Source Data

As machine learning algorithms are data-driven, variations in the data could greatly affect the accuracy of a classification [37]. When providing insufficient data, such as for data scheme A (merely optical remote sensing data from Sentinel-2), poor results were obtained (overall accuracy of 0.46~0.47 for scheme A in Table 11). Due to influences by environmental conditions such as different angles and intensities of sunlight, topography, water content, and other factors, it is possible that the same object may have different spectrums, and different objects may have the same spectrum [38,39]. For instance, the spectral information for the same vegetation on sunny and shady slopes could be different, and the height of lower vegetation is easily obscured by shadows. This tends to increase the errors during training. Complementarity of multiple sources of data is generally used to address this problem.

Compared to the data scheme A, the data scheme D (addition of DEM data and some ground survey data) significantly increased the overall accuracy by 0.29~0.36, with the accuracy rates for the “good”, “medium”, and “poor” categories increasing by 0.12~0.19, 0.09~0.1, and 0.32~0.48, respectively. The addition of DEM data complements vertical structure parameters which are lacking in optical remote sensing data [40], allowing areas of deciduous trees to be distinguished from areas of vegetation with similar spectral characteristics (e.g., high-density grassland) [41], which in turn has an impact on the accuracy of the results. The addition of ground survey data further complements the growth status information for the vegetation—such as the age of trees that can lead to changes in their growth rate—thereby improving the accuracy of the model.

4.3. The Feasibility of Machine Learning

Different from the traditional statistical methods used by Huafu Liu et al. [42], Hailong Yin et al. [43], and Kassim et al. [17], this paper exploited a machine learning algorithm to develop a comprehensive model to evaluate forest ecological function levels, which has a higher flexibility and faster processing speed for high-dimensional data with more complex relationships among feature factors. Most indicators are non-linearly related to forest ecological function levels, including ground survey factors and spectral characteristics. For instance, with the increasing NL, the forest ecological function levels increase first and then decrease. Machine learning algorithms are non-linear approximations to an objective function, different from than the traditional comprehensive evaluation methods bound to a linear function, such as the scoring method, principal component analysis, etc. Therefore, their powerful fitting ability could make the predicted results closer to the reality and improve the evaluation accuracy of the model.

The Random Forest showed the best performance among the above three models, with an overall accuracy of 0.82 and an F1 score of 0.62. As a non-linear, parametric classifier, Random Forest is robust with non-equilibrium data and can randomly generate multiple decision trees to form a forest, effectively avoiding overfitting [44,45]. It allows the fusion of high-dimensional data from multiple sources [46] and has a high tolerance for missing values and outliers, so that it can effectively reduce the interference of noise in the data. In addition, it can automatically determine the importance of variables, which in turn improves its accuracy and usability.

4.4. Limitations of this Study

This study, based on the multi-source data of Sentinel-2 remote sensing images and DEM and partial data from the forest resource planning and design survey, as well as the three machine learning algorithms RF, LightGBM, and CatBoost, evaluated the forest ecological function levels. Overall, our results may promote research on the evaluation of forest ecological function levels. However, the 10 m resolution of the Sentinel-2 images may limit a further improvement of the performance metrics. If higher-resolution images can be acquired in the future, such as remote sensing images from Gaofen series satellites or UAV images, it will be possible to further increase the performance of our model and even further reduce the participation of ground survey indicators. On the other hand, from the perspective of research methods, deep learning algorithms such as YOLO are also worth a try in the future.

5. Conclusions

Optical remote sensing data, DEM data, and forest resource planning and design survey data were used in this study to evaluate the forest ecological function levels of Lin’an District using three machine learning algorithms, i.e., RF, LightGBM, and CatBoost.

In the three models, Random Forest was the best-performing model, with an overall accuracy rate of 0.82 (the accuracy rates for the “good”, “medium”, and “poor” categories being 0.76, 0.89, and 0.83, respectively) and with an F1 score of 0.62.

The multi-source data significantly improved the performance metrics. Furthermore, the acquisition of ground survey data such as QI_YUAN, LING_ZU, LIN_ZHONG, YOU_SHI_SZ, and NL, was achieved at lower costs than those required for the traditional eight indicators of forest biomass, forest naturalness, forest community structure, tree species structure, vegetation coverage, forest canopy density, mean tree-height, and thickness of dead leaves.

If more data sources are used, such as higher-resolution remote sensing images, LiDAR remote sensing images, etc., the estimation performance might further improve in the future.