*Article* **Ecosystem and Driving Force Evaluation of Northeast Forest Belt**

**Zhihong Liao 1,†, Kai Su 1,\*,†, Xuebing Jiang 2, Xiangbei Zhou 1, Zhu Yu 3, Zhongchao Chen 4, Changwen Wei 1, Yiming Zhang <sup>1</sup> and Luying Wang <sup>1</sup>**


**Abstract:** The ecosystem in the Northeast Forest Belt (NFB) can provide various ecosystem services, such as soil conservation, habitat provision, water conservation, and so on. It is essential for maintaining the ecological environment in Northeast China and the entire country. In the face of increasingly severe environmental problems, the comprehensive and accurate evaluation of ecosystem conditions and their changes is significant for scientific and reasonable recovery and protection measures. In this study, the NFB was taken as the research area. The spatio-temporal changes in ecological quality from 2005 to 2015 and the main driving factors behind them were analyzed by constructing the comprehensive ecosystem evaluation index. The results showed that: The landscape types of the NFB were mainly forest, cropland, and grassland. And the better ecological environment of the NFB was mainly distributed in the south of Changbai Mountains (CBM), the middle of Lesser Khingan Mountains (LKM), and the northwest of Greater Khingan Mountains (GKM). In contrast, the northeast of CBM, the southwest of LKM, and the edge of southern GKM were relatively poor. During 2005–2015, the ecosystem in the NFB was in a relatively good state as a whole, showing a steady-to-good development trend. However, more attention needed to be paid to some areas where degradation still existed. Land use/cover, climate (annual average rainfall, etc.), and human disturbance were potential factors affecting ecosystem evolution in the NFB. This study aims to provide an effective scientific basis and policy reference for the environmental protection and construction of the NFB.

**Keywords:** land use; land cover; temporal and spatial change; ecosystem assessment; GIS; northeast forest belt

#### **1. Introduction**

As the largest ecosystem on land, the forest ecosystem can provide a variety of ecosystem services (ES) such as soil conservation service (SCS), carbon sequestration service (C), habitat provision service (HP), sand-stabilization service (SSS), water conservation service (WCS) and so on. The Northeast Forest Belt (NFB) contains China's typical and essential forest ecosystem [1]. Due to the rich forest and biological resources, it is an important forest resource and biodiversity protection base in China. At the same time, because of the fertile soil and rich mineral resources, the NFB is an important food base and an old industrial base in China. To further protect the environment and promote regional sustainable development, China proposed the strategic pattern of building "two ecological barriers and three shelters" ecological stability in 2006 [2].

However, under the dual influence of global climate change and human disturbance, a series of ecological environment problems have been triggered, affecting the ecological environment stability in Northeast China. Understanding ecosystem evolution and driving factors is significant for ecosystem protection and sustainable development [3–10]. In the

**Citation:** Liao, Z.; Su, K.; Jiang, X.; Zhou, X.; Yu, Z.; Chen, Z.; Wei, C.; Zhang, Y.; Wang, L. Ecosystem and Driving Force Evaluation of Northeast Forest Belt. *Land* **2022**, *11*, 1306. https://doi.org/10.3390/ land11081306

Academic Editors: Jay Mar D. Quevedo, Norie Tamura, Yuta Uchiyama and Ryo Kohsaka

Received: 18 July 2022 Accepted: 11 August 2022 Published: 13 August 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/).

face of increasingly serious ecological problems, it is urgent to scientifically evaluate the ecosystem status, changing trends, and driving factors of the NFB. With the construction of the pattern of "two ecological barriers and three shelters" in China, many studies of the ecosystem of the NFB have also been carried out. Sun [1] assessed the ecosystem pattern, quality, and service function of the NFB by remote sensing data combined with monitoring statistical data through model simulation and sample surveys, and analyzed the changing trend between 2000 and 2010. It was found that there was no significant change in ecosystem quality (EQ) and service function of the NFB during the period. Still, the interference of social activity factors was increasing continuously. Based on this study, the calculation model of net ecosystem production was constructed to evaluate the carbon sequestration in the forest ecosystem of NFB [2]. It was concluded that the overall performance of carbon sequestration was carbon sinks, and human disturbances such as urban expansion, poor vegetation growth, and high temperature were significant factors affecting carbon sequestration. It provided an effective scientific basis for formulating ecological and environmental protection measures and policies in the NFB. After that, Su et al. [3] analyzed the changes in landscape pattern in the NFB from 2000 to 2015, and pointed out that the ecosystem in the NFB remained stable as a whole for 15 years, and simulated the changing trend in 2020. Qi et al. [11] explored the spatial differentiation characteristics and mechanism of trade-offs and synergies among six different ES in the NFB, indicating that the six ES showed a significant aggregation distribution, and the trade-offs and synergies had apparent spatial differentiation. Zhu et al. [12] evaluated the spatial-temporal distribution and changes in ecological vulnerability in the NFB in two different periods. They concluded that the overall ecological vulnerability was at a good level. Net primary productivity and land use types were the main driving factors affecting the spatial differentiation pattern of ecological vulnerability. This series of studies have provided a certain theoretical and scientific basis for the ecological system protection and management and its sustainable development in the NFB. However, previous studies have not constructed a comprehensive evaluation index for the NFB, so as to conduct a more integrated and comprehensive evaluation of its ecosystem.

Due to the limitations of single-factor evaluation results, more and more studies have begun to construct remote sensing ecological index (RSEI) to comprehensively evaluate and analyze the ecosystem changes in different regions [13–16]. In recent years, analyzing ecosystem evolution based on RSEI is also a hot topic worldwide. Mohammad et al. [17] used RSEI and impervious surface percentage feature space to map the urban surface ecological poorness zone (USEPZ) by combining land surface temperature (LST), humidity, normalized difference vegetation index (NDVI), and normalized difference soil index (NDSI), and thus proposed a new method for quantifying USEPZ. The results showed that the significant differences between surface ecological status (SES) and USEPZI in different cities were mainly caused by the physical characteristics of surface organisms. Based on this study, Mohammad et al. [18] continued to develop a new surface ecological condition index– land surface ecological status composition index (LSESCI) and compared it with RSEI. It was concluded that LSESCI was superior to RSEI in simulating the spatial and temporal changes of SES. At the same time, Karbalaei et al. [19] also synthesized RSEI under the PSR (pressure-state-response) framework by using principal component analysis combined with NDVI, LST, land surface moisture (LST) and normalized differential built-up, and bare soil index (NDBSI). The results showed that the disturbance of human activities and climate change disturbed the ecological balance, resulting in the reduction of urban EQ. Then based on RSEI, Maity et al. [20] used Landsat image data to evaluate the temporal and spatial variation of ecological environment quality in Kolkata urban agglomeration. The results accurately described the ecological environment quality of the study area, which could help the policy makers to scientifically guide ecological management decisions. In China, compared with other relevant studies, Xu et al. [21] early proposed a new remote sensing ecological index by coupling four evaluation indexes: vegetation index, surface temperature, humidity component, and soil index, and integrating each index with

principal component transformation. Based on previous studies, the proposed remote sensing ecological index was improved, and the time series of ecological status images were generated by using the improved index and surface temperature images to realize the detection of ecological changes at different scales [22]. Then Yuan et al. [23] constructed remote sensing ecological index by integrating NDVI, humidity, LST, and NDBSI, and analyzed the temporal and spatial changes of EQ in Dongting Lake Basin from 2001 to 2019. Wang et al. [24] selected four indicators of humidity, greenness, dryness, and heat, and used the principal component analysis method to construct the ecological environment quality evaluation index. Combined with the land use data of Xinjiang and the central urban area of Urumqi, they explained the ecological environment quality and land use status in each period, thus revealing its dynamic development trend. Although the RSEI has high effectiveness and reliability and is widely used in different regions and ecosystem assessments, it ignores the impact of land use/ cover change on the structure and function of the ecosystem [23,25,26]. Single-factor (or specific-factor) assessments are oriented and can be used to assess ecosystem conditions from a specific perspective. Although multi-factor (pattern-quality-function) multi-level assessments can systematically and comprehensively assess ecosystems, the results of different levels of assessments have opposite conclusions among certain regions [27–29].

To sum up, using a comprehensive evaluation index to evaluate the changes in the ecosystem in each region effectively reduces the uncertainty based on single-factor evaluation results. In addition, combined with the analysis of various influencing factors, we can more comprehensively and accurately understand the impact of natural and human disturbance on the ecosystem. Therefore, this study used remote sensing (RS) and GIS technology, based on the "pattern-quality-service" evaluation framework system, using the analytic hierarchy process (AHP) to determine the weight of each index, and constructed a comprehensive index to evaluate the ecosystem of the NFB. The purpose of this study was: (1) To assess the long-term dynamic changes of EQ in the NFB from 2005 to 2015; (2) To explore the driving forces of these changes. This study can provide decision support and a scientific basis for the ecosystem stability of the NFB.

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

#### *2.1. Study Area*

The study area is located in the NFB. According to its geographical characteristics, it can be divided into three regions: Greater Khingan Mountains (GKM) area, Lesser Khingan Mountains (LKM) area, and Changbai Mountains (CBM) area. The NFB is located in 118◦48- ~134◦22- East longitude and 40◦52- ~53◦34- North latitude, with a total area of about 670,000 km2, of which the forest area accounts for about 400,000 km2. The terrain is diverse, mainly mountainous and hilly, the landscape is undulating, and the soil is fertile. The climate belongs to the temperate monsoon climate, hot and rainy in summer, dry and less rain in winter, and relatively cold; four seasons are distinct; annual precipitation is about 300~1000 mm. As the ecological stability barrier in Northeast China, the NFB is a forest region with a high conservation value in Northeast China [3,12]. The location of the study area is shown in Figure 1.

**Figure 1.** Location of the study area.

#### *2.2. Data Source and Processing*

The DEM data used in this study is original elevation products with STRMDEM 90 m resolution, downloaded from the Geospatial Data Cloud Platform of the Computer Network Information Center of the Chinese Academy of Sciences (http://www.gscloud.cn) and accessed on 3 March 2022. The land use/cover data used are MODIS land use products from the NASA official website (https://ladsweb.nascom.nasa.gov/search), accessed on 3 March 2022. The soil data is derived from the harmonious world soil database China subset of the National Qinghai-Tibet Plateau/Third Pole Environmental Data Center, which includes soil texture, sandy soil, silt, clay, and organic carbon. Climate data comes from the National Meteorological Information Center (http://data.cma.cn) (accessed on 7 March 2022), which includes annual precipitation and temperature, relative humidity, sunshine hours, etc. The Data Centre for Resources and Environmental Sciences (RESDC) of the Chinese Academy of Sciences (http://www.resdc.cn) (accessed on 7 March 2022) provides population density and GDP data (spatial resolution 1 km). The road data used are downloaded on OpenStreetMap. The road type is national highway and expressway. To facilitate spatial analysis and comparison, all data used in this study is resampled to grid data with a 1 km resolution.

In this study, after completing the collection and pretreatment of various types of data, based on the "pattern-quality-service" evaluation framework system, the AHP was used to determine the weight coefficient of each index. Then ArcGIS was used to calculate the ecosystem comprehensive evaluation index, which was used to comprehensively evaluate the ecosystem evolution in the NFB. At the same time, the forward and reverse conversion index (PNTIL) model was used to analyze the evolution of its ecosystem pattern. The transfer matrix was used for quantitative spatio-temporal analysis of EQ, and the geographical detector model was used to analyze its driving force. The technical roadmap is shown in Figure 2.

**Figure 2.** Technology roadmap.

#### *2.3. Methods*

#### 2.3.1. Construction of Ecosystem Comprehensive Evaluation Index

Changes in land use patterns will directly affect the type and intensity of services provided by ecosystems. As the supply area of ES such as WCS, C, and HP, the introduction of ES accounting into land use decision-making in the NFB can better promote the rational development of natural resources and achieve sustainable land use. In addition, as an important part of the ecosystem, the in-depth assessment of the growth status and vitality of vegetation should also be an essential parameter for ecosystem assessment. Through the analysis of various influencing factors, the comprehensive evaluation index is used to evaluate the changes in the ecosystems in various regions, which can more comprehensively and accurately understand the influence of natural and human disturbance on the ecosystem. Therefore, this study used a "quality-service-pattern" assessment framework to evaluate the evolution and driving forces of the NFB ecosystems. Based on the ecosystem's "quality-service-pattern" evaluation framework, a comprehensive index evaluation model (Table 1) was constructed. The indicators were calculated using remote sensing, land use, other data, and related ecological process models.


The weight coefficient was determined by the AHP [30–32]. The main steps of AHP were as follows: Firstly, the hierarchical model of the NFB ecosystem was established, which was divided into three levels, including goal layer A, labeling layer B, and indicator layer C. Secondly, seven experts in this study scale the factors at each level and constructed judgment matrix, that is, the number 1~5 and its reciprocal were used as the scale to objectively judge the factors of each level. Subsequently, the eigenvector and the maximum eigenvalue of the matrix were obtained, and the consistency test of the matrix was carried out to obtain the weight coefficients. In this study, *CR* < 0.1 passed the consistency test. The judgment matrix of different levels is shown in Table 2. The consistency test index formulas are [33,34]:

$$CR = CI / RI \tag{1}$$

$$CI = \frac{\lambda\_{\text{max}} - n}{n - 1} \tag{2}$$

where *C I* is the deviation consistency index; *RI* is the average consistency index; *λ*max is the maximum eigenvalue; *n* is the order of the judgment matrix; ensuring that the ratio of *CI* to *RI* is less than 0.1 [34].


**Table 2.** Judgment matrix of different levels.

Note: A is the goal layer; B is the labeling layer, and B1–B3 are ecosystem quality, ecosystem service, and ecosystem pattern; C is the indicator layer, C1–C8 are FVC, NPP, LAI, WCS, SCS, SSS, C, HP; Wi is the weight value.

Then the obtained weight coefficients are used to calculate the comprehensive index. The calculation method is [35]:

$$V\_{\rm cei} = \sum\_{i=1}^{n} V\_i \times \varepsilon\_i \tag{3}$$

where, *Vcei* is the composite index; *Vi* is the evaluation result of indicator *i*; *ε<sup>i</sup>* is the weight of the first index.

#### 2.3.2. Ecosystem Service

One of the important purposes of the construction of the NFB is to improve regional ES and promote regional ecological stability. The rich ES assessment in the NFB is added to the assessment framework system, which can better realize the protection, improvement, and sustainable development of the ecosystem. Therefore, this study selected SSS, SCS, WCS, C, and HP five ES for analysis and evaluation.

SSS is usually estimated by the amount of sand-stabilization of vegetation. By calculating the difference between the potential and the actual wind erosion, the sand stabilization of vegetation can be estimated. The potential wind erosion refers to the wind erosion when there is no vegetation, and the actual wind erosion refers to the wind erosion when there is vegetation. The calculation formulas are [36]:

$$SR = S\_{L1} - S\_{L2} \tag{4}$$

$$S\_{L1} = \frac{2 \cdot x}{S\_1^2} Q\_{\text{max}1} \cdot e^{-\left(\frac{x}{S\_1}\right)^2} \tag{5}$$

$$S\_{L2} = \frac{2 \cdot x}{S\_2^{\ast 2}} Q\_{\text{max2}} \cdot e^{-\left(\frac{x}{S\_2}\right)^2} \tag{6}$$

$$Q\_{\text{max1}} = 109.8 \left[ \text{WF} \times EF \times \text{SCF} \times K' \right] \tag{7}$$

$$S\_1 = 150.71 \left( WF \times EF \times SCF \times K' \right)^{-0.5711} \tag{8}$$

$$Q\_{\text{max2}} = 109.8 \left[ WF \times EF \times SCF \times K' \times COG \right] \tag{9}$$

$$S\_2 = 150.71 \left(\text{W}F \times EF \times \text{SCF} \times \text{K'} \times \text{COG}\right)^{-0.3711} \tag{10}$$

where *SR* is the sand-stabilization, kg/m2; *SL*<sup>1</sup> and *SL*<sup>2</sup> are potential wind erosion and actual wind erosion, kg/m2, respectively; *Q*max1 and *Q*max2 are the maximum migration capacities of potential and actual soil erosion, kg/m; *S*<sup>1</sup> and *S*<sup>2</sup> are the critical field lengths of potential and actual soil erosion, respectively, m; *x* is the maximum wind erosion distance, m; *WF* is the meteorological factor; *EF* is the soil erodibility factor; *SCF* is a soil crust factor; *K*is the surface roughness factor; *COG* is the vegetation cover factor.

SCS is usually estimated by the soil conservation of vegetation. The soil conservation of vegetation can be estimated by calculating the difference between the potential and the actual soil conservation. The potential soil conservation refers to the amount of soil erosion generated without vegetation, and the actual soil conservation refers to the amount of soil erosion generated with vegetation. This study selected USLE to evaluate SCS in the NFB ecosystems. The calculation formula is as follows [37]:

$$\text{SC} = \text{SE}\_p - \text{SE}\_a = \text{R} \cdot \text{K} \cdot \text{LS} \cdot (1 - \text{COG}) \tag{11}$$

where, *SC* is soil conservation, [t/(hm2 a)]; *SEp* and *SEa* are potential and actual soil conservation, respectively, [t/(hm2 a)]; *<sup>R</sup>* is rainfall erosivity factor, MJ·mm/(hm2·ha); *<sup>K</sup>* is soil erodibility factor, t·hm2·h/(hm2·MJ·mm); *LS* and *COG* are terrain factor, vegetation cover factor, no dimension.

WCS is usually estimated by the water balance equation. WCS can be estimated using precipitation minus evaporation, evapotranspiration, and rainstorm runoff. The calculation formula is as follows [38]:

$$WR = PRE - ET - QF \tag{12}$$

where *WR* is water conservation, mm; *PRE* is annual precipitation, mm; *QF* is rainstorm runoff, mm. *ET* is the actual evaporation and emission, mm.

*C* is usually calculated using above-ground biomass multiplied by the biomass-carbon conversion coefficient. The main calculation formula is as follows [39]:

$$\text{COS} = \sum\_{i=1}^{j} AGB\_i \times \text{C}\_i \tag{13}$$

where *COS* is the ecosystem carbon storage; *i* is the type *i* ecosystem; *j* is the total number of ecosystem types; *AGBi* is the above-ground biomass of ecosystem type *i*; *Ci* is the biomass-carbon conversion coefficient of *i* ecosystem types.

HP is usually estimated by the distribution of indicator species in county units. The calculation formula is [11]:

$$E\_{\mathbf{a}} = \sum\_{i=1}^{n} E\_i \tag{14}$$

where *E*<sup>a</sup> indicates the total number of species (number); *Ei* is the number of species (number).

#### 2.3.3. Ecosystem Quality

Vegetation is an important part of the ecosystem, which connects the natural processes such as the atmosphere, water, and soil. Its change will directly affect the regional climate, hydrology, and soil conditions, and has an important impact on the regional energy cycle

and the biochemical cycle of substances. It is an important indicator of regional EQ change. It is also significant to add the in-depth evaluation of vegetation growth status and growth vigor into the evaluation framework system. Therefore, this study selected fractional vegetation coverage (FVC), net primary productivity (NPP), leaf area index (LAI), three EQ analyses, and evaluations.

FVC refers to the ratio of the vertical projection area of plants in a region to that of the region, which is one of the most essential indicators to measure EQ and surface vegetation status. Normally, the FVC is usually calculated by NDVI. The calculation formulas are [40]:

$$\text{NDVI} = (\rho\_{\text{NIP}} - \rho\_{\text{R}}) / (\rho\_{\text{NIP}} + \rho\_{\text{R}}) \tag{15}$$

where ρNIP and ρ<sup>R</sup> are the reflectivity of near-infrared band and red band, respectively.

$$\text{FVC} = (\text{NDVI} - \text{NDVI}\_{\text{soil}}) / \left( \text{NDVI}\_{\text{veg}} - \text{NDVI}\_{\text{soil}} \right) \tag{16}$$

where FVC is fractional vegetation coverage, NDVIveg is the NDVI value of pure vegetation pixel, and NDVIsoil is the NDVI value of no vegetation pixel.

NPP reflects the efficiency of fixing and converting light energy into compounds in plants. NPP can well reflect the vegetation productivity and EQ. At present, there are many models to calculate NPP. In this study, the CASA model estimates the NPP of the NFB. In the CASA model, NPP is mainly determined by the two variables of APAR and light energy conversion rate (*ε*) absorbed by vegetation. The calculation formulas are [41]:

$$\text{NPP} = \text{APAR}(t) \times \varepsilon(t) \tag{17}$$

$$\text{APAR} = \text{FPAR} \times \text{PAR} \tag{18}$$

$$
\varepsilon(t) = \varepsilon \times T1(t) \times T2(t) \times \mathcal{W}(t) \tag{19}
$$

where APAR(*t*) is the photosynthetically active radiation absorbed by plants; *ε*(*t*) is the conversion efficiency of APAR into organic carbon; FPAR is the effective absorptivity of APAR by plants; PAR is the driving energy of plant photosynthesis. *ε* is the maximum light energy conversion rate under ideal conditions; *T*1(*t*) and *T*2(*t*) are the effects of temperature and photosynthesis in varying degrees, respectively; *W*(*t*) is the influence coefficient of water stress.

LAI refers to the ratio of the total leaf area of plants to the land area per unit of land area. LAI can reflect the vegetation growth in vertical structure, which is also one of the important indicators for evaluating EQ. There are many methods to calculate LAI, among which remote sensing technology can realize dynamic real-time and large-scale monitoring of LAI. The vegetation canopy radiative transfer model is one of the common models to calculate LAI by using remote sensing technology. The corresponding radiative transfer process is described as follows [42]:

$$-\mu \frac{\partial L(z,\Omega)}{\partial \tau} + G(\tau,\Omega)L(z,\Omega) = \frac{\omega}{4\pi} \int\_{4\pi} p\left(\Omega' \to \Omega\right) G\left(\Omega'\right) L(z,\Omega') d\Omega' \tag{20}$$

where *L* is brightness; *μ* = *cosθ* is the chord value of transmission direction zenith angle; *G* is the leaf angle distribution function; *ω* is the leaf albedo; *p*(Ω- → Ω) is the vegetation canopy phase function.

#### 2.3.4. Evolution of Ecosystem Pattern

The land use/cover change has an important function in maintaining ecosystem service function. Land use/cover change will affect the function and structure of the ecosystem, and the evolution of land use/cover mode will directly affect the types and intensity of ES. It is significant to study the evolution of ecosystem service value driven by land use/cover change, which is also an essential quantitative index of environmental effects of land use/cover change. Therefore, this study selected the ecosystem transfer matrix and PNTIL model to analyze and evaluate the land use/cover change in the NFB.

The ecosystem transfer matrix can effectively list the conversion relationship between different ecosystems in different periods in the form of a matrix to reflect the change characteristics of ecosystems and the flow direction between systems in detail, and quantify the conversion status between systems [43–46]. Using the following formula to calculate the proportion of different types of landscape areas [47]:

$$P\_{\vec{\eta}} = S\_{\vec{\eta}} / TS(\vec{\iota} = 1, \, 2, \, \dots, \, n; \, \vec{\eta} = 1, \, 2, \, \dots, \, n) \tag{21}$$

where *Pij* is the area ratio of type *i* ecosystem based on classification at all levels in the first year of the ecosystem classification system; *Sij* is the area of category *i* ecosystems in the ecosystem classification system at all levels in the year *j*; *TS* is the total area of the evaluation area; *n* is the number of ecosystem types.

PNTIL model can analyze the evolution trend of ecosystem patterns [47–49]. Based on the stability and biodiversity of each ecosystem, the ecosystem transformation was divided into positive change and reverse change, and the positive/reverse transformation rules of the ecosystem were established (Table 3). The ecosystem transfer matrix of the NFB from 2005 to 2015 was used to quantitatively evaluate the ecosystem PNTIL of the NFB. The calculation formula of the ecosystem conversion index is [47]:

$$LCTR\_{i,k} = (\Delta S\_{i,k} / S\_i) \times 1/t \times 100\tag{22}$$

where *LCTRi*,*<sup>k</sup>* is the conversion rate of ecosystem type to *k* in the *i* region; *Si* is the total area of all ecosystems in region *i*; Δ*Si*,*<sup>k</sup>* is the total area converted into *k* type; *k* = 1 represents the positive conversion rate of ecosystems in the region; *k* = 2 represents the inverse conversion rate of the ecosystem, *t* being a time variable.


**Table 3.** Rules of forward/reverse ecosystem transformation.

2.3.5. Quantitative Spatio-Temporal Analysis of EQ

A quantitative spatio-temporal analysis of EQ in the NFB was carried out. The results of EQ were divided into five grades, and the transfer matrix of five grades of EQ in three periods was counted to reflect the changes in EQ in time and space. Using the following formula to calculate the proportion of different levels of EQ area [47]:

$$P\_{\vec{\eta}} = \mathcal{S}\_{\vec{\eta}} / TS(i = 1, \ 2, \ \ldots, n; \ j = 1, \ 2, \ \ldots, n) \tag{23}$$

where *Pij* is the area ratio of type *i* EQ in year *j* based on classification at all levels in the EQ classification system; *Sij* is the area of category *i* EQ in year *j* based on classification at all levels in the EQ classification system; *TS* is the total area of the evaluation area; *n* is the number of types of EQ at different levels.

#### 2.3.6. Ecosystem Driving Force Analysis

The geo-detector model is mainly a new statistical method to reveal the driving factors behind spatial differentiation, which can detect whether a certain factor is a reason for the evolution of the ecosystem. This model mainly uses the q-value measurement factor to detect the results. The larger the q-value, the stronger the driving effect of independent variables for the spatial differentiation of dependent variables and vice versa [50–55]. Based on the availability of relevant research [12,35] and regional data, 18 potential socio-ecological driving factors were selected as independent variables X. The comprehensive evaluation indexes of the ecosystem in the three periods of 2005, 2010, and 2015 were selected as the dependent variable Y, including terrain-related factors (elevation, slope, aspect), climate (annual average precipitation, annual average temperature), social economics (population density, GDP, road), and land cover (cropland, forest, grassland, wetland, built-up land, desert, other types), soil texture (silt, clay, sand). Next, the natural breakpoint method was used to change the independent variable X from numerical value to category quantity. The independent variable X was matched with the dependent variable Y. Then, the geographical detector model was used for factor detection analysis. Finally, the influence value and driving effect of each variable X for the dependent variable Y were obtained.

#### **3. Results and Analysis**

#### *3.1. Evolution of Ecosystem Pattern*

From 2005 to 2015, the classification result of ecosystem patterns in the NFB is shown in Figure 3 and Table 4. It can be seen that forest, grassland, and cropland were the main ecosystems. Statistics showed that the area of these three ecosystems was about 580,000 km2, accounting for more than 94%. Wetland, built-up land, desert, and other ecosystems were scattered in the NFB. During 2005–2015, the structure of ecosystems in the NFB was relatively stable, but conversions between different types were frequent due to climate change, human activities, and other interference factors. Ecosystem transfer in the NFB during 2005–2015 is shown in Figure 4.

**Figure 3.** Spatial distribution of different land use types in the Northeast Forest Belt (NFB) from 2005–2015.


**Table 4.** Statistics of different land use types of ecosystem patterns of the NFB from 2005–2015.

**Figure 4.** Area transfer matrix chord of landscape pattern types in the NFB from 2005–2015.

From 2005 to 2010, the areas of cropland and forest decreased, and the areas of grassland, wetland, and built-up land increased. Cropland was mainly converted into forest, about 378 km2, followed by grassland, about 219 km2. However, some grassland and forest were still transformed into cropland, with a conversion area of about 537 km2. This showed that the project of returning cropland to forest and grassland in the important ecological function areas of the NFB had achieved good results. However, there were still some cases of deforestation and increasing human disturbance.

From 2010 to 2015, large areas of grassland, forest, and wetland were transformed into cropland and built-up land. During this period, cropland conversion was the main driver of the increase in built-up land areas, with about 153 km2 of cropland converted into built-up land. However, there were still 951 km2 of grassland, forest, and wetland transformed into cropland. During 2010–2015, the area of cropland and built-up land increased steadily, and the interference of social factors such as human activities was significantly higher than that between 2005–2010.

Between 2005 and 2015, there were more mutual transformations between cropland and forest, forest and grassland, cropland and grassland, built-up land and cropland, cropland and wetland. Among them, the forest type was transferred out of 1214 km2, mainly to cropland and grassland. The area transferred to other land was relatively small, indicating that the protection of forests had a certain effect in 10 years. A total of 854 km<sup>2</sup> area of grassland changed, mainly converted to forest and cropland. The outflow of cropland mainly flowed to forest, grassland, built-up land, and wetland, with areas of 449 km2, 260 km2, 230 km2, and 112 km2, respectively. The built-up land area increased yearly, and the conversion speed was gradually accelerating. The increasing area mainly came from cropland, forest, and grassland. The wetland area mainly flowed to cropland, and other landscape types had almost no change. Overall, the project of returning cropland to forest and grassland had achieved a good result. Still, the areas of cropland and the built-up land areas had gradually increased in the past 10 years, indicating that the impact of human interference, such as the blind expansion of cultivated land and construction land on the NFB, was gradually increasing, which required us to pay more attention to and the protection of its ecosystem needed to be continued.

Forward/reverse conversion rates (TFR and RFR) of ecosystems in the NFB for 2005–2010 and 2010–2015 are shown in Table 5. During these two periods, the RFR of the forest and wetland was greater than TFR, indicating that these two ecosystem types were more transformed to lower ecosystem types. The TFR of built-up land and other ecosystem types was higher than that of RFR, indicating that ecological control projects had also achieved good results and progress in the NFB. During the two periods, RFRs of the forest, grassland, and wetland were high and showed an upward trend between 2010 and 2015, indicating that there were still some problems to be solved in the NFB. At the same time, among all ecosystems, the TFR of built-up land, cropland, and other ecosystem types was significantly higher than that of the other four ecosystems, indicating that the implementation of policies such as returning cropland to forest and grassland and limiting the scope and extent of economic activities had a great impact on the protection and improvement of ecosystems.


**Table 5.** Forward/reverse conversion rate in the NFB.

Note: TFR is the forward transformation rate, 100%; and RTR is the reverse transformation rate, 100%.

#### *3.2. Spatial Patterns and Variation of EQ*

The results of NPP, LAI, and FVC were divided into five grades (excellent, good, middle, general, and poor) according to the natural breakpoint method. The spatial distribution map of EQ and the area transfer matrix of each EQ grade were obtained by spatial statistical analysis, and the quantitative spatial and temporal analysis of EQ was carried out. The spatial distribution of NPP, LAI, and FVC at different levels in the NFB from 2005 to 2015 is shown in Figure 5.

**Figure 5.** Spatial distribution of three types of EQ at different levels in the NFB from 2005–2015. Note: FVC is fractional vegetation coverage, NPP is net primary productivity, and LAI is leaf area index.

#### 3.2.1. Dynamic Characteristics of NPP

It can be seen from Figure 5 that in 2005, the areas with high NPP were mainly distributed in the south of GKM, the north of GKM, the north of LKM, and the east of CBM, while other areas were mainly dominated by middle NPP. In addition to the decline in northern GKM in 2010, other regions had increased significantly; however, by 2015, the GKM, the LKM, and the northeast of CBM showed a sharp downward trend, and only the central part of CBM remained stable, and the southern part improved. Figure 6 and Table 6 showed that NPP in the NFB was mainly middle-quality and good-quality areas. During 2005–2010, the area of general and middle grades showed a significant downward trend, while that of poor, good, and excellent grades showed an increasing trend. The area ratio of good and excellent NPP increased rapidly. The largest transformation area was the middle grade, mainly to good grade, and the transfer area was 127,125 km2. In general, the period showed a trend of transformation in a good direction. During 2010–2015, the overall trend was opposite to that before, and the areas of poor, good, and high-quality grades mainly showed a downward trend, while those of middle and general grades showed an increasing trend. The largest transformation area was the good grade, mainly middle grade, and the transfer area was 135,889 km2. The overall period showed a certain degree of quality decline. This indicated that the NPP in the NFB increased from 2005 to 2010, and it may decrease from 2010 to 2015 due to ecological damage and environmental degradation.

**Figure 6.** Sankey diagram of the area transfer matrix of different NPP levels in the NFB from 2005–2015.


**Table 6.** Area transfer matrix of different NPP levels in the NFB from 2005–2015.

3.2.2. Dynamic Characteristics of LAI

It can be seen from Figure 5 that the distribution of high LAI areas was similar to that of FVC, which was mainly distributed in the middle of the three mountains in the NFB, and the edge areas of northeast CBM, southwest LKM and southern GKM were relatively low. In 2015, the area with an LAI of greater than 48 was mainly distributed in the central part of LKM and the central and southern part of CBM, with an area of 206,334.50 km2, accounting for 33.50%. The area with an LAI of less than 32 was mainly located northeast of CBM and on the edge of the southern GKM, with an area of 146,782.5 km2, accounting for 23.82%.

It can be seen from Figure 7 and Table 7 that in addition to the small area of poor quality, the other four grades were almost similar, and the area of excellent grade had always been the highest. During 2005–2010, areas at the general and middle levels showed a decreasing trend, while areas at the poor, good and excellent levels showed an increasing trend. The largest transformation area was from fair grade to middle grade, and the transfer area was 44,022.50 km2. The period generally showed a trend of transformation in a good direction. During 2010–2015, the area of some middle and general grades increased and decreased, but the area of excellent and good grades continued to increase. The largest transformation area was from middle to fair grade, and the transfer area was 127,125 km2. The overall trend was relatively stable. It showed that the LAI of the NFB had been continuously improved and enhanced during 2005–2015, and the policy of returning cropland to forest and grassland had protected and improved its ecosystem.

**Figure 7.** Sankey diagram of the area transfer matrix of different LAI levels in the NFB from 2005–2015.



3.2.3. Dynamic Characteristics of FVC

It can be seen from Figure 5 that the high-level regions were mainly distributed in most parts of the NFB, and most of them were concentrated in the middle of the three mountains of the NFB. The areas with low FVC were mainly distributed in the edge area of the southern GKM and the northeast of CBM. The cropland area in this area was large, and the FVC value was low. During 2005–2015, except for the southern edge of the GKM, FVC decreased, and the other regions showed an upward trend.

It can be seen from Figure 8 and Table 8 that the area with good and excellent grades quality FVC accounted for a large proportion, accounting for about 81% of the total area. During 2005–2010, the area of excellence level had increased significantly, and its proportion exceeded half of the total area. The poor grade mainly turned to the general grade; the general grade was mainly transferred to the middle grade; the middle grade mainly turned to the good grade; the good grade mainly turned to the excellent grade, and the excellent grade mainly turned to a good grade. During 2010–2015, the area at the excellent and middle levels continued to increase, while the area at the poor and general levels decreased; the poor grade mainly turned to the general grade, the general grade primarily turned to the middle grade, the middle grade mainly turned to the good grade, the good grade mainly turned to the excellent grade, the excellent grade mainly turned to good grade. During 2005–2015, it generally showed a trend of good transformation, indicating that FVC in the NFB had been in a better situation during this period and had continued to improve and enhance.

**Figure 8.** Sankey diagram of area transfer matrix of different FVC levels in the NFB from 2005–2015.



#### *3.3. Spatial Patterns and Variation of ES*

#### 3.3.1. Habitat Provision Service

The high value of HP was mainly distributed in the middle of GKM, LKM, and CBM, while the low value was mainly distributed around them (Figure 9). Overall, HP in forest areas was higher than other spatial pattern types such as cropland and built-up land and showed a decreasing trend from inside to outside. From the point of view of the time change, it was mainly to maintain a relatively stable state, the overall change was small, from 2005 to 2015 showed a slight upward trend, and 2010–2015 compared to 2005–2010 upward trend was more obvious. During 2005–2010, the overall decreasing and increasing areas accounted for about 37.60 % and 45.05%, respectively. The decreasing areas were mainly located where cropland and built-up land were located, while the increasing areas were mainly located where LKM and CBM forests were distributed. From 2010 to 2015, the areas where HP increased and decreased accounted for about 22.54% and 55.02%, respectively. The increased areas were mainly located in the eastern and southern parts of the GKM and the areas where cropland was distributed. The main reason may be that forest protection measures such as returning cropland to forest and grassland were implemented there, while the reduced areas were mainly distributed in some areas of the western GKM.

**Figure 9.** Spatial distribution of five types of ES at different levels in the NFB from 2005–2015. Note: HP is habitat provision service, SCS is soil conservation service, C is carbon sequestration service, SSS is sand-stabilization service, and WCS is water conservation service.

#### 3.3.2. Soil Conservation Service

The SCS of the NFB was generally low. The high value of SCS was mainly distributed in the south of CBM, and the middle of the three mountains was also mainly about the middle value. Most areas of SCS were low value (Figure 9). It remained stable over time, with very small overall changes, with a slight upward trend in 2005–2010 and a slight downward trend in 2010–2015.

#### 3.3.3. Carbon Sequestration Service

The spatial distribution of C in the southeast was higher than that in the northwest (Figure 9). The highest value appeared in the south of CBM, and the lowest value appeared in the southwest of GKM. For the time change, the whole area of GKM had a low value in 2005, and then the value in the eastern and northern regions improved. From 2005 to 2010, the increased area of C (71.68%) was significantly larger than the decreased area (27.42%), showing an upward trend. The increase of C was mainly concentrated in the northern GKM and the whole area of LKM and CBM, and only the southwestern GKM showed a significant downward trend. The southwestern GKM had more cropland area, so this sharp decline may be related to increased disturbance of human activities during the period. From 2010 to 2015, the decrease in the C area (54.75%) was greater than the increase (44.36%), showing a downward trend in general. The C value in the southwest GKM increased, and the decrease was mainly distributed in the northern GKM and the southern CBM.

#### 3.3.4. Sand-Stabilization Service

The high SSS value was mainly distributed in the northeast of CBM and the southwest of GKM, while the low SSS value was mainly distributed in the north of GKM, the middle of LKM, and the south of CBM (Figure 9). During 2005–2010, SSS showed a downward trend, but from 2010–2015 showed an upward trend. From 2005 to 2010, the proportion of SSS increased (11.03%) and decreased (15.99%) relatively small, most of which remained stable (72.98%). From 2010 to 2015, the increased area (19.36%) was improved compared with the previous, but most areas remained relatively stable (73.02%). During the two periods, the increased or decreased areas were mainly distributed in the areas where other spatial pattern types except forests were located, while the areas where forests were located remained relatively stable.

#### 3.3.5. Water Conservation Service

The high value of WCS was mainly distributed in the southern region of the CBM. The middle of the three mountains of the GKM, the LKM, and the CBM also maintained a middle value, while the low value was mainly distributed in the northeast of the CBM, the southwest of the LKM, and the southern edge of the GKM. The WCS of the whole area of the GKM was lower (Figure 9). From the perspective of the time change, from 2005 to 2010, WCS showed a significant upward trend, and 59.75% of the regions increased. The regions with obvious increase were the central GKM, the central LKM, and the southern CBM. From 2010 to 2015, the WCS showed a significant downward trend, and the decreased area (54.13%) was greater than the increased area (23.90%). The decreased area was mainly distributed in the central part of the GKMs and the southern part of the CBM, especially on the southern edge of the CBM.

#### *3.4. Comprehensive Ecosystem Index Evaluation*

Firstly, the weight coefficients of the relevant indicators of the ecosystem in the NFB were determined by the AHP, which passed the consistency test, and the weight setting was reasonable. Then, ArcGIS was used to obtain the spatial distribution map of the comprehensive evaluation index of the ecosystem in different years. Finally, it was divided into five grades according to the natural breakpoint method and manual adjustment, and the data were statistically analyzed (as shown in Figure 10 and Table 9).

**Figure 10.** Spatial distribution of different levels of the comprehensive index in the NFB from 2005–2015.


**Table 9.** Statistics of different levels of comprehensive index data of the NFB from 2005–2015.

The spatial and temporal distributions of the ecosystem comprehensive evaluation index in the NFB in 2005, 2010, and 2015 are shown in Figure 10. The areas with a high spatial distribution value of the ecosystem in the NFB were mainly distributed in the southern part of CBM, the central part of LKM, and the northwestern part of GKM, while the northeastern part of CBM, the southwestern part of LKM and the marginal area of southern GKM were relatively low. During 2005–2010, the value of the southern edge of GKM showed a significant decline, and the eastern and central CBM and the central LKM also showed a slight decline. During 2010–2015, the change in the ecosystem comprehensive evaluation index showed an almost opposite trend; that is, the values of southern GKM and central LKM increased significantly, and other regions also increased to varying degrees. During the whole study period from 2005 to 2015, except for the decrease in the value around the marginal area in the southern GKM, the value in other regions showed an upward trend, which was closely related to a series of protection and improvement measures such as the "two ecological barriers and three shelters" national ecological stability pattern strategy, grain for green project, and comprehensive management of soil and water loss carried out in the NFB.

According to the change level of the ecosystem comprehensive evaluation index (0.013–0.746, Table 9), the change of the ecosystem comprehensive evaluation index was divided into five grades: excellent, good, middle, general, and poor. The area of the five grades of the ecosystem comprehensive evaluation index in 2005, 2010, and 2015 is shown in Table 4. Specifically, in 2005, the area with good and middle grades accounted for a large proportion (total 63.24%), and the proportion of poor grade was very small (5.21%). This indicated that the overall quality of ecosystems in 2005 was good. The grade area ratio in 2010 was consistent with 2005, and the overall area changed little. However, on the whole, some areas fell from the middle, good and excellent grades to the poor and general grades, indicating that during this period, the NFB faced some problems such as ecological destruction and environmental degradation. In 2015, the good grade area was the largest (230,204 km2), and the poor and middle-grade areas significantly decreased. Good and

excellent areas had been significantly increased, mainly from poor, general, and middle regions, indicating the effectiveness of the "two ecological barriers and three shelters" national ecological stability pattern strategy implemented in the NFB.

#### *3.5. Driving Force Analysis*

In this study, a 10 km × 10 km grid was used to generate 7119 points in the study area in ArcGIS. After analysis, the influence value and driving effect of each index factor on the ecosystem evolution in the NFB were obtained (as shown in Figure 11).

**Figure 11.** Geo-detector model results of the original index system in the NFB. Note: 1 is elevation, 2 is slope, 3 is slope direction, 4 is temperature, 5 is precipitation, 6 is NPP, 7 is cropland, 8 is forest, 9 is wetland, 10 is grassland, 11 is built-up land, 12 is clay, 13 is silty sand, 14 is sand, 15 is soil type, 16 is population density, 17 is GDP, and 18 is road. \*\*\* indicates that the correlation is significant at the 0.001 level; \*\* expresses relevance at the 0.01 level; \* indicates that the correlation is significant at the 0.05 level.

It can be seen from Figure 11 that land use/cover, climate (annual average rainfall, etc.), and human disturbance (GDP, etc.) were related potential factors affecting the ecosystem evolution of the NFB. Most of the driving factors in the three periods of 2005, 2010, and 2015 were sufficient for the evolution of the NFB ecosystem (*p* < 0.05), indicating that the evaluation system constructed in this study was reasonable. In 2005, the three most influential factors in the evolution of the NFB ecosystem were forest, cropland, and elevation, and their influence values were 0.299, 0.287, and 0.223, respectively. The influence of driving factors from high to low were terrain factors, soil factors, meteorological factors, and human social factors. The influence of some driving factors in 2010 and 2015 differed from 2005, but the impact of the overall driving factors showed little change. Land use types had

the greatest influence on ecosystem evolution in the NFB, and cropland and forest had a more significant influence than wetland, grassland, built-up land, and other types. The influence value of NPP remained at a middle-to-high level during the study period, and its influence value was more obvious than other factors. Regarding topographic factors, the influence value of elevation and slope was relatively significant, while the influence of slope direction was not big. In terms of soil factors, the influence of sand was larger than the other three factors. The influence of temperature on climate factors had been small, while the influence of precipitation is more obvious; even in 2010, it was only lower than that of the forest. In terms of social factors, the influences of human activity factors such as GDP and population density were small. Still, their visibility was extremely significant (*p* < 0.001, except for the GDP factor in 2005), indicating that human disturbance factors such as GDP and population density were also potential factors affecting the evolution of the NFB ecosystem.

#### **4. Discussion**

#### *4.1. Advantages of the Integrated Ecosystem Assessment Index*

In this study, the AHP was used to obtain the comprehensive evaluation index of the ecosystem, and a variety of methods were used to evaluate the ecosystem of the NFB comprehensively. The results showed that the landscape types of the NFB were mainly forest, cropland, and grassland, and the study showed a stable state as a whole. However, the areas of cropland and built-up land were still expanding, which needed more attention. This result was consistent with the research of others [3]. Compared with the principal component analysis method [56–60] that has been studied more before, the results showed that the ecosystem comprehensive evaluation index created by AHP in this study can also carry out a comprehensive ecosystem assessment. The results showed that the areas with a good ecosystem environment in the NFB were mainly distributed in the south of CBM, the middle of LKM, and the northwest of GKM and had a significant improvement and improvement trend with time. In contrast, the areas around the southeast and southwest of GKM were relatively poor, similar to some previous research results [2,12]. The reason may be that the good region was mainly located in the main part of the three mountains and was dominated by forests, with lush vegetation and less human disturbance, so the comprehensive index was high. The poor area may be due to the northeast plain, arable land to the NFB continued to expand, and more serious interference from human activities, resulting in the poor composite index. Throughout 2005–2015, the ecosystem in the NFB was in a relatively good state as a whole, and maintained a trend of improvement, which was closely related to a series of protection and improvement measures such as the "two ecological barriers and three shelters" national ecological stability pattern strategy [61], the project of returning cropland to forest and grassland [62–64], and the comprehensive control of soil and water loss. By comparing the results of landscape pattern, ES, and EQ in the NFB in this study, it could be seen that the evaluation results based on a single factor were uncertain and limited. And there were conflicts and contradictions between the evaluation results of different factors. Thus the ecosystem cannot be accurately evaluated. Based on these shortcomings and deficiencies, some researchers chose the method of constructing RSEI to analyze the ecosystem evolution [23,65,66] to obtain more accurate and comprehensive results. However, the construction of the remote sensing ecological index will ignore the impact of the land use/cover change process on the structure and function of the ecosystem. Therefore, in this study, combined with multiple indicators such as land use/cover, ES, and EQ, it was more accurate, systematic and comprehensive to evaluate the evolution of the NFB ecosystem by constructing the ecosystem comprehensive evaluation index.

#### *4.2. Uncertainty of Driving Force Analysis*

The driving force analysis showed that the natural factors had a great influence on the ecosystem of the NFB. Among them, land use type was an important factor affecting the ecosystem evolution of the NFB, which was consistent with the research results of Su and Zhu et al. [3,12]. However, the influence value of NPP in the study of Zhu et al. was more significant than in this study, which may be related to the difference in data sources and the number of grid points. The influence of elevation and slope was larger, probably because the higher and steep areas were covered by vegetation or less affected by human activities; the ranking of precipitation had been ahead, and its influence value rose to second in 2010, which was consistent with some studies [67–69]. Precipitation was one of the main factors affecting vegetation growth, so it also affected the evolution of ecosystems. The impact of human activity factors such as population density and GDP on the evolution of the ecosystem in the NFB was small, consistent with the results of Zhu et al. [12]. This may be because the overall social and economic activities in the NFB were less, and the population density and GDP were relatively sparse. In addition, this study only selected population density, GDP, and road as the three human activity factors for analysis. It did not comprehensively and specifically analyze what human social factors will affect the evolution of the NFB ecosystem, resulting in uncertainty in the final analysis of the results. Currently, there is no consistent method to evaluate the driving force that affects the evolution of the ecosystem. The methods for analyzing the driving force of the ecosystem mainly included correlation analysis, residual or redundancy analysis, and different methods may produce different results. Zhu et al. [12] found that the influence of social factors such as GDP and population density on the ecosystem of the NFB was not obvious; Sun et al. [3] indicated that human disturbance factors such as urban expansion were the main reason for the decrease of carbon sequestration in the NFB. As for the drought index, Wang et al. [70] considered that it was the main factor driving the growth of forest vegetation and the change of carbon sink function in Northeast China. The results of driving force analysis in this study were similar to those in Zhu et al. [12], although there were still differences in some factors and effects. The selection of driving force indicators and methods may have an important impact on the final results and the subsequent management of the ecosystem. Therefore, the accurate analysis of driving force still needs further continuous research.

#### *4.3. Limitations of Current Research*

This study still had some deficiencies and needed to be further studied. Although the landscape pattern of the NFB ecosystem analyzed in this study could roughly reflect the changes during the entire study period. It cannot accurately reflect the evolution of the NFB ecosystem pattern and the relationship between horizontal and vertical pattern structure and quality function. This still needs more detailed research. And the remote sensing data used in this study, due to sources and other limitations, only selected 2005, 2010, and 2015 data, and remote sensing data resolution was relatively low, which may affect the accuracy of the results. In addition, this study only conducted quantitative spatio-temporal assessments of ES and EQ in the NFB. In the future, the collaboration and trade-off between ES and quality will be explored. In addition, in the driving force analysis, due to the limitation of data sources, only three human activity factors, namely GDP, population density, and road, were selected for analysis. However, these three factors were relatively macroscopic, which cannot reflect the specific human disturbance activities that affect the evolution of the ecosystem in the NFB, which was not conducive to the application in a specific practice. In future research, more comprehensive driving factors can be added to analyze, and the evaluation system framework can be further optimized to more comprehensively and accurately evaluate the changes in the ecosystem environment and the impact of driving factors.

#### **5. Conclusions**

In this study, the multi-source remote sensing data and multiple indicators were used to construct the comprehensive evaluation ecosystem index of the NFB. The changes and driving forces of the NFB ecosystem from 2005 to 2015 were evaluated based on the PNTIL model, transfer matrix, geographical detector, and other methods. Using an ecosystem comprehensive evaluation index can avoid the shortcomings of uncertainty and limitation based on single factor evaluation results to evaluate the evolution of the NFB ecosystem more systematically, comprehensively, and accurately. The main conclusions are as follows:

The landscape pattern types of the NFB were mainly forest, cropland, and grassland, and the overall landscape pattern of the NFB was stable during 2005–2015. The comprehensive assessment of the ecosystem index showed that the areas with a good ecosystem environment in the NFB were mainly distributed in the southern CBM, the central LKM, and the northwestern GKM, while the northeastern CBM, the southwestern LKM, and the marginal areas of southern GKM were relatively poor. During 2005–2015, the overall ecosystem of the NFB was in a relatively good state and maintained a trend of improvement in a good direction. However, during 2005–2010, there were still some anthropogenic disturbances that led to a decline in the comprehensive index of some regions. The driving force analysis showed that natural factors had a great influence on the ecosystem of the NFB. Land use/cover was an important driving factor affecting the ecosystem evolution of the NFB. Climate (annual average rainfall, etc.) and human disturbance (GDP, etc.) were potential factors affecting the ecosystem evolution of the NFB. The results showed that in the past 10 years, the EQ of the NFB was generally in a good state and continued to improve, but the expansion of cropland and built-up land in some areas continued to increase, which needed to be paid more attention.

In the future, relevant research will continue to analyze the synergy and trade-off relationship of various ES and their spatial differentiation and further regionalize ES to explore the formation mechanism of the synergy and trade-off relationship of different ES. The correlation between synergy and trade-off relationship with the main driving forces obtained in the driving force analysis of this study can also be further studied. At the same time, the optimization of the evaluation system framework and the reasonable and comprehensive selection of impact factors still need to be further improved, including the prospect of the change characteristics of the ecosystem comprehensive index with the above factors. It is hoped that a series of studies will provide a more scientific and comprehensive reference for the ecosystem stability and improvement of the NFB.

**Author Contributions:** Conceptualization, Z.L. and K.S.; methodology, Z.L. and K.S.; software, Z.L.; validation, X.Z., Z.Y. and Z.C.; formal analysis, Z.L. and K.S.; resources, K.S.; data curation, L.W., C.W., Z.L. and Y.Z.; writing-original draft preparation, Z.L. and K.S.; writing-review and editing, K.S. and X.J.; visualization, L.W., C.W., Z.L. and Y.Z.; supervision; project administration, K.S.; funding acquisition, K.S. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by [Youth Science Foundation of Natural Science Foundation of Guangxi] grant number [2022GXNSFBA035570] and [Talent Introduction Program of Guangxi University] grant number [A3360051018].

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The DEM data, meteorological data, soil data, land cover data are available in this research, upon any reasonable request, by emailing the authors.

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

#### **References**

