Discover Ecosystem Service Value Growth Characteristics of a Subtropical Soil Erosion Area Using a Remote Sensing Driven Mountainous Equivalent Factor Method

1. Introduction

Soil erosion is a naturally occurring process that affects all landforms and threats ecosystem viability [1]. It results in soil property changing, land degradation, vegetation destruction, agricultural productivity drop and ecosystem service function declination [2,3]. Correspondingly, a comprehensive ecological restoration is practiced to control soil erosion, and brings vegetation recovery, reforestation, land upgradation, and soil reformation, which underpin ecosystem services [4]. Ecosystem services are the benefits human populations derive, directly or indirectly, from ecosystem functions [5]. The monetary value of ecosystem services has been explored in the past decades. Following the academic studying [6,7], international forums have conducted it, such as Millennium Ecosystem Assessment (MEA), Intergovernmental Science-Policy Platform on Biodiversity and Ecosystem Services (IPBES) and The Economics of Ecosystems and Biodiversity (TEEB) [8].

Among various monetization methods for estimation ecosystem service value (ESV) [9,10], two kinds of approaches as the primary data based and the unit-value based are utilized widely [5,11,12]. The primary data based approaches combine physical quantities of ecosystem services with market price, travel cost, opportunity cost, shadow price, and ecosystem services replacement cost, so they are complicated and require a number of parameters. Such methods are often performed on one or a few kinds of service, and usually are applied on small area or single ecosystem [13,14]. The unit-value based approaches use the monetized value per unit area of ecosystems [5]. They are suitable for the ESV evaluation under land use changes from urbanization, ecological restoration and so on [15,16,17]. One of them is the equivalent factor method, in which the framework modified by Xie et al uses equivalence coefficients to reflect the relative weights of the provisioning, regulating, habitat, and cultural service value for a certain ecosystem [18]. Such unit-value based method is widely used in different scales of the local, municipal, regional, provincial and national [19,20,21,22,23].

Besides the classic indicator of soil and water conservation rate, ESV can be used to evaluate ecological restoration for soil erosion regions under the background of “lucid waters and lush mountains are invaluable assets”. For example, Changting county, China, is a typical ecological restoration case selected at the 15th Conference of the Parties to the United Nations Convention on Biological Diversity (COP15). After several decades’ efforts in soil and water conservation driven by government [24], especially, the intensified efforts since 12^th Five-Year Plan, the county has made significant improvements in stemming soil erosion and ecological restoration. The soil and water conservation rate was improved, specifically, it was increased from 89.74% to 93.43% in the period from 2010 to 2022. However, it is unknown ESV variation from 2010 to 2022, ESV spatial distribution pattern, and what contribution of land use changing and vegetation recovery in ESV variation.

In addition, the existing equivalent factor method can be improved for accurate estimation of ESV in subtropical mountainous soil erosion area [25,26]. One important factor is spatial adjustment coefficient that illuminates the spatial variation of the ESV from any kind of land use. In most cases, net primary productivity (NPP) is used to represent spatial adjustment coefficient [27]. However, NPP mainly stands for the carbon deposited in the stem, branches, and roots of forest and underlying soil, but undervalues the canopy of vegetation, e.g., green leaves for air quality and climate regulation, and aesthetic landscape. Meanwhile, the spatial resolution of the NPP products is usually 500 m or 1000 m, which is difficult to assist the ESV evaluation for small patches of soil erosion controlling. Therefore, a spatial adjustment coefficient with more physical based and higher resolution will be preferable for more accurate monetization of regional ESV [11,28]. The other one is that topographic effect hinders accurate classification of land uses and estimation of vegetation canopy in mountains using remote sensing, due to the phenomena of “same body with different spectra” and “different bodies with same spectrum” caused by anisotropic solar illumination [29,30,31]. Surface reflectance and vegetation indices are always under estimated in topographic shadow, including self-shadow and cast shadow, where the pixels are obstructed from solar direct irradiance [32,33]. Therefore, accurate remote sensing parameters removed topographic effect are crucial for land uses classification, vegetation canopy monitoring, and ESV estimation in mountainous soil erosion area.

The main objectives of this study are to (1) evaluate ecological restoration in a subtropical mountainous soil erosion area from an ESV perspective, and obtain ESV variation in spatial-temporal pattern, and contributions from land uses changing and vegetation canopy recovery; and (2) develop a remote sensing driven mountainous equivalent factor (RS-MEF) method for comprehensive assessment ESV in rugged mountains, which coupled a modified equivalent factor framework, the remote sensing techniques for mountains, and carbon sink and vegetation canopy.

2. Study Area and Data

2.1. Study Area

The study area is Changting county (116°00′45″-116°39′20″E, 25°18′40″-26°02′05″N), Longyan city, east China’s Fujian Province (Figure 1). It covers land elevation ranging from 156 m to 1443 m with a mean of approximately 513 m and a slope ranging from ~1° to 88° with a mean of approximately 18°. The major land cover in this area is forest, which has prominent topographic effect in remote sensing images. Due to ecological vulnerability [34,35], Changting county had suffered intensive soil and water erosion in last century, which restricted local sustainable development. According to local survey, most of soil and water erosion locates in seven towns in the middle-east of the county, i.e., major soil erosion area (MSEA). After several decades’ efforts in soil and water conservation, especially, the intensified efforts since 12^th Five-Year Plan, it has made significant and continuous improvements in stemming soil erosion and ecological restoration. The soil and water conservation rate was improved from 68.55% in 1985 to 93.43% in 2022. At the 15th Conference of the Parties to the United Nations Convention on Biological Diversity (COP15), the soil erosion control practices in Changting were selected as a typical case of ecological restoration.

2.2. Data Source

Multiple source data for research area were collected, including satellite images, unmanned aerial vehicle (UAV) images, digital elevation model (DEM), field survey data, NPP open product, Changting Statistical Yearbook 2023, Compilation of National Agricultural Cost-benefit Data 2023 (Table 1). We acquired one scene of Landsat 5 TM in 29 October, 2010 and the other one of Landsat 9 OLI in 22 October, 2022, which are multi-spectral images with 30 m spatial resolution. The corresponding DEM data with 30 m resolution of ASTER GDEM V2 were also obtained. These data were all downloaded from the Computer Network Information Center, Chinese Academy of Sciences (Geospatial Data Cloud. Available online: http://www.gscloud.cn (accessed on 26 October, 2023)). The NPP product was downloaded from the Earth Science Data Systems (ESDS, Available online: https://appeears.earthdatacloud.nasa.gov/ (accessed on 6 May, 2024)). In addition, we made a field survey in 18 to 21 September, 2022 to check soil and water conservation, land use and vegetation canopy. Especially, we surveyed three sub-areas of Huangwuqian (HWQ), Laiyoukeng (LYK), and Xianggongting (WGT), using the DJI Phantom 4 pro multi-spectral UAV (Figure 1).

3. Methods

The RS-MEF method integrates a modified equivalent factor framework [18] and the remote sensing techniques for mountains, such as the integrated topographic correction (ITC) [30] and the shadow-eliminated vegetation index (SEVI) [32]. Its main procedures are image processing including topographic correction for spectral reflectance and vegetation index, and ESV calculation including a standard equivalent factor calculation, land uses classification, equivalent coefficients set for diverse ecosystem services from different land uses, and adjustment coefficients specification for dynamic spatial variation of ESV. Two major improvements of the RS-MEF method are spatial adjustment coefficient using NPP coupled with the SEVI and land uses classification using the satellite images removed topographic effect (Figure 2). Finally, a spatial-temporal analysis of ESV is conducted, specifically, evaluation the ESV variation since 12^th Five-Year Plan, i.e., from 2010 to 2022.

3.1. Image Processing

One feature of the RS-MEF method is utilization of the remote sensing techniques for mountains, including the SEVI coupled with NPP for spatial adjustment coefficient, and the ITC for land uses classification. The SEVI is calculated using the surface reflectance of the red band and the near-infrared (NIR) band (Eq. 1), in which the adjustment factor is calculated using the block information entropy algorithm (BIE-algorithm) [36]. First, a DEM of the study area is divided into blocks of 6 km × 6 km, and the steep blocks are selected from the 5% highest slope mean of the blocks. Then, the SEVI information entropy of any steep block is calculated (Eq. 2–3), and the target block that achieves the highest entropy is used to determine the final adjustment factor (Eq. 4–5). Finally, the SEVI of 2010 and 2022 is normalized using the fixed land objects and statistics.

$$\mathrm{S}\mathrm{E}\mathrm{V}\mathrm{I}={\displaystyle \frac{{\rho }_{nir}}{{\rho }_r}}+f\left(\mathsf{\Delta }\right)·{\displaystyle \frac{1}{{\rho }_r}}$$

(1)

$$H=-{\displaystyle \frac{{\sum }_{i=1}^n{p}_i·\mathrm{l}\mathrm{n}\left(p_i\right)}{\ln \left(n\right)}}$$

(2)

$$p_i={\displaystyle \frac{{SEVI}_i}{{\sum }_{i=1}^n{SEVI}_i}},{\sum }_{i=1}^n{p}_i=1,0<p_i<1$$

(3)

$$f_b=\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{a}\mathrm{x}\left(H\right),f(∆)\in (0.001,1.000)$$

(4)

$$f_s=\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{a}\mathrm{x}\left(H_b\right),f(∆)\in (f_{bj},j=1,\dots ,m)$$

(5)

where SEVI is the shadow-eliminated vegetation index, ${\rho }_{nir}$ is the near-infrared band reflectance, ${\rho }_r$ is the red band reflectance, $f\left(∆\right)$ is the adjustment factor, $H$ is the information entropy of SEVI in a block, $p_i$ is the percentage of a SEVI pixel value in a selected block, $n$ is the number of pixels in a selected block, ${SEVI}_i$ is a pixel value of SEVI, $f_b$ is an optimized adjustment factor for a block, $f_s$ is the final adjustment factor for an entire scene, $H_b$ is the maximum information entropy of SEVI in a block, and $m$ is the number of selected steep blocks in an entire scene.

The ITC is a data driven topographic correction method, including the SEVI information and the reflectance simulated from a physical based model (Figure 3), e.g., the Sun-canopy-sensor (SCS) + C correction (Eq. 6–7). Firstly, the total topographic shadow in the SCS+C corrected image is extracted. Then, randomly generated pixels in sunny areas are used as training samples, including dependent variable from the SCS + C corrected reflectance, and independent variable of the SEVI. Finally, the regressions generated from the Random Forest (RF) regressor between the SEVI and the waveband reflectance of sunny areas are used to correct the waveband reflectance of topographic shadow of the RGB wavebands. As for the NIR waveband, the inverse calculation method is used (Eq. 8), since the SEVI degrades into the ratio vegetation index (RVI) when topographic effect is eliminated in an image.

$${\rho }_{\mathrm{S}\mathrm{C}\mathrm{S}+\mathrm{C}}={\rho }_{\mathrm{T}}·\left({\displaystyle \frac{\cos \theta ·\cos \sigma +c}{\cos i+c}}\right)$$

(6)

$$\cos i=\mathrm{c}\mathrm{o}\mathrm{s}\sigma ·\mathrm{c}\mathrm{o}\mathrm{s}\theta +\mathrm{s}\mathrm{i}\mathrm{n}\sigma ·\mathrm{s}\mathrm{i}\mathrm{n}\theta ·\mathrm{c}\mathrm{o}\mathrm{s}\left(\beta -\omega \right)$$

(7)

where ρ_SCS+C is the reflectance observed for a horizontal surface, ρ_T is the reflectance observed over an inclined surface, and c, which equals the quotient of intercept b and inclination a of an observed empirical linear correlation between ρ_T and the cos i, is assumed to be constant for a given wavelength; i is solar incidence angle, defined as the angle between the normal to the ground and sun’s rays; and θ is solar zenith angle, σ is slope angle, β is aspect angle, and ω is solar azimuth angle.

$${\rho }_{nir-tc}=\mathrm{R}\mathrm{V}\mathrm{I}·{\rho }_{r-tc}=\mathrm{S}\mathrm{E}\mathrm{V}\mathrm{I}·{\rho }_{r-tc}$$

(8)

where ρ_nir-tc is the topographic corrected reflectance of the NIR waveband, RVI is the ratio vegetation index, and ρ_r-tc is the topographic corrected reflectance of the red waveband.

3.2. ESV Calculation

(1)

Standard equivalent factor

An empirical formula for standard equivalent factor was calculated (Eq. 9), using the data from the Changting Statistical Yearbook 2023 and the Compilation of National Agricultural Cost-benefit Data 2023. In this study, the standard equivalent factor of 2010 was modified using the consumer price index (CPI) to match the benchmark of that of 2022. So, these prices come more in line with the present currency than the standard equivalent factors modified to match that of one earlier year, e.g., 2000 or 2010.

$$E={\displaystyle \frac{1}{7}}\times {\sum }_{i=1}^{\mathrm{n}}{\displaystyle \frac{a_i{p}_i{q}_i}{A}}\times C$$

(9)

where E is a standard equivalent factor (CNY/hm²), a_i is the area of i-th crop in a study year (hm²), p_i is the average price of i-th crop in a study year (CNY/kg), and q_i is the crop yield of i-th crop in a study year (kg/hm²), A is the total area of all crops in a study year (hm²), C is the adjustment coefficient of consumer price index (CPI).

(2)

Equivalent coefficients

Equivalent coefficients provide a weight coefficient of valuation for each ecosystem service. An expert-based equivalent coefficients table was referenced (Table 2) [18].

(3)

Land uses

Based on the land use types in Table 2 and the actual situation of study area, the Landsat images were classified as coniferous forest, broad-leaved forest, mixed forest, grassland, farmland, construction land, unused land, and water. Considering the prominent topographic effect in Landsat images, we used the ITC rectified images to classify land uses. Since the ITC is better to remove topographic shadow, especially cast shadow in rugged mountains [30]. Firstly, the training samples are selected from the images, with the aids from field survey and the Google Earth (Table 3). Next, the features of spectral bands, vegetation indices, textures, and primary component calculated from the images are set and trained using a machine learning of the RF technique. Then, the classification accuracy is assessed using confusion matrix with overall accuracy and Kappa coefficient. Finally, eight land uses are classified to match the equivalent coefficients.

(4)

Spatial adjustment coefficient

Considering vegetation indices mainly reflect vegetation canopy of forest and grass [37], we improved spatial adjustment coefficient using NPP coupled with vegetation index to incorporate the stem, branches, roots and canopy of forest in mountains (Eq. 10). Specifically, we used the SEVI to improve the spatial adjustment coefficient for rugged mountains in this study (Eq. 11), since it has additional benefits of topographic effect elimination, and weak target recognition in green mountains [38]. In addition, the spatial resolution of the revised spatial adjustment coefficient is improved to 30 m from 500 m of the open NPP product.

$$S_i=\sqrt{{NPP}_i\times {VI}_i}$$

(10)

$$S_{i,30m}=\sqrt{{{NPP}_{i,500m}\times SEVI}_{i,30m}}$$

(11)

where $S_i$ is the spatial adjustment coefficient of i-th land use, ${NPP}_i$ is the NPP of i-th land use, and ${VI}_i$ is the vegetation index of i-th land use; $S_{i,30m}$ is the spatial adjustment coefficient of i-th land use with 30 m, ${NPP}_{i,500m}$ is the downloaded NPP product of i-th land use with 500 m, and ${SEVI}_{i,30m}$ is the SEVI of i-th land use with 30 m.

(5)

ESV calculation

Finally, an improved equivalent factor framework is developed (Eq. 12). Then, the RS-MEF is suggested to calculate the ESV of mountains (Eq. 13), which is a hybrid of the improved equivalent factor framework, the SEVI, and the classified land uses from the Landsat images removed topographic effect. Considering the “general labor” from ecosystems, i.e., general ecosystem service (resource), we used the spatial adjustment coefficient mean of non-major soil erosion area (n-MSEA) to match the standard equivalent factor. The n-MSEA is the county area subtracted the MSEA.

$$ESV=\sum \left(E{\times T}_i\times A_i\times S_i\right)=\left(E{\times T}_i\times A_i\times \sqrt{{NPP}_i\times {VI}_i}\right)$$

(12)

$$ESV=\sum \left(E{\times T}_i\times A_i^{’}\times S_{i,30m}^{’}\right)=\sum \left(E{\times T}_i\times A_i^{’}\times {\displaystyle \frac{S_{i,30m}}{\overline{\sqrt{{{(NPP}_{n,500m}\times SEVI}_{n,30m}}}}}\right)=\sum \left(E{\times T}_i\times A_i^{’}\times {\displaystyle \frac{\sqrt{{{(NPP}_{i,500m}\times SEVI}_{i,30m}}}{\overline{\sqrt{{{(NPP}_{n,500m}\times SEVI}_{n,30m}}}}}\right)$$

(13)

where ESV is ecosystem services value (CNY), E is standard equivalent factor (CNY/hm²), T_i is the conventional equivalent coefficients of i-th land use, $A_i$ is the area of i-th land use (hm²), $A_i^{’}$ is the area of i-th land use classified from the image removed topographic effect (hm²), $S_{i,30m}^{’}$ is a relative spatial adjustment coefficient of i-th land use, ${NPP}_{n,500m}$ is the NPP of n-MSEA, and ${SEVI}_{n,30m}$ is the SEVI of n-MSEA.

3.3. ESV Analysis

Using the calculated ESV result, we make a spatial-temporal analysis to discover the spatial distribution pattern of ESV and the ESV changing characteristics from 2010 to 2022. In spatial pattern analysis, we divide the county as MSEA and n-MSEA for comparison ESV. In addition, topographic gradient effect is analyzed using the topographic indicators of elevation, slope, relief amplitude, and terrain niche, with five grades of them set using quantile approach (Table 4). The ESV of land use is also analyzed to recognize the major contributor. In temporal analysis, the ESV increment and growth rate from 2010 to 2022 are computed. In order to achieve the contribution from land use changing in ESV variation, equation 13 is inversed to get a growth rate of land use changing contribution from 2010 to 2022 (Eq. 14).

$${\mathrm{C}}_{\mathrm{l}\mathrm{u}}={\displaystyle \frac{{T_{2022}\times A}_{2022}^{’}}{{T_{2010}\times A}_{2010}^{’}}}-1={\displaystyle \frac{{ESV}_{2022}/\left(E_{2022}\times S_{2022}^{’}\right)}{{ESV}_{2010}/\left(E_{2010}\times S_{2010}^{’}\right)}}-1$$

(14)

where C_lu is a growth rate of contribution from land use changing in ESV variation from 2010 to 2022.

4. Results

4.1. Standard Equivalent Factor

The standard equivalent factor for Changting county, China was calculated as 2858.42 CNY/hm² in 2010 and 5535.70 CNY/hm² in 2022, respectively.

4.2. Land Uses

The overall accuracy and the Kappa coefficient of land uses classification are 0.812 and 0.780 in 2010, and 0.819 and 0.782 in 2022, respectively. The classified result shows that major land uses are mixed forest, coniferous forest, grassland, and farmland (Figure 4 and Figure 5). From 2010 to 2022, the areas of mixed forest and coniferous forest in MSEA have enlarged from 393.97 km² to 433.74 km² and from 283.27 km² to 380.56 km², with area increments of 39.77 km² (growth rate: 10.09%) and 97.29 km² (growth rate: 34.35%), respectively. Meanwhile, those in n-MSEA have enlarged from 744.44 km² to 805.81 km² and from 436.73 km² to 523.31 km², with area increments of 61.37 km² (growth rate: 8.24%) and 86.58 km² (growth rate: 19.82%), respectively. In general, the largest increased area (and growth rate) in the county is coniferous forest with an increment of 183.87 km² (from 720.00 km² to 903.87 km²) and growth rate of 25.54%. The second is mixed forest with an increment of 101.14 km² (growth rate: 8.88%), enlarged from 1138.41 km² to 1239.55 km². On the contrary, the areas of grassland and farmland have shrunk in the same period from 475.81 km² to 268.07 km², and from 346.09 km² to 304.92 km², with area decrements of 207.75 km² (growth rate: −43.66%) and 41.17 km² (growth rate: −11.90%), respectively. The matrix of land use change shows that mixed forest was majorly changed from coniferous forest. Meanwhile, the enlargement of coniferous forest area comes from multi-lands, e.g., unused land, grassland, and mixed forest. According to the expert-based equivalent coefficients (Table 2), the land use changing from coniferous forest to mixed forest can make a great increasement of ESV.

4.3. Spatial Adjustment Coefficients

The adjustment coefficients of NPP, SEVI and $\sqrt{NPP\times SEVI}$ are displayed in Figure 6. The adjustment coefficient of NPP illuminates rough pixels (500 m) with a mean of 2010 as 789.09 gC/m²a and that of 2022 as 827.13 gC/m²a (growth rate: 4.82%), and a general low-medium value in MSEA. Meanwhile, the adjustment coefficient of SEVI shows smooth pixels (30 m) with a mean of 2010 as 0.56 and that of 2022 as 0.64 (growth rate: 12.95%), and a low-medium value in MSEA. Apparently, the value distribution of SEVI is more continuous than that of NPP. Finally, the adjustment coefficient of $\sqrt{NPP\times SEVI}$ illustrates the merged features in spatial resolution and value distribution, with a mean of 2010 as 20.16 and that of 2022 as 21.82 (growth rate: 8.23%). Specifically, the$\sqrt{NPP\times SEVI}$ mean of MSEA was 18.51 in 2010 and 21.08 in 2022, with a growth rate of 13.88%; meanwhile, that of n-MSEA was 21.51 in 2010 and 22.43 in 2022, with a growth rate of 4.25%. From 2010 to 2022, the growth rate of$\sqrt{NPP\times SEVI}$ of MSEA was 3 times higher than that of n-MSEA.

4.4. ESV and Distribution Characteristics

Using the equation 13, the ESV of Changting county was estimated (Figure 7). The total ESV and the ESV per unit area in 2010 were 15.80 billion CNY and 50913 CNY/hm², and those in 2022 were 34.83 billion CNY and 112244 CNY/hm², respectively. Apparently, low ESV locates mainly in downtown and MSEA. Specifically, the ESV per unit area of MSEA was about 43885 CNY/hm² in 2010 and 103186 CNY/hm² in 2022, less than that of n-MSEA with about 56606 CNY/hm² in 2010 and 119523 CNY/hm² in 2022.

From 2010 to 2022, the total ESV of the county increased 120.47%. Specifically, the total ESV of MSEA in 2010 and 2022 was 6.14 billion CNY and 14.44 billion CNY, with growth rate of 135.13%. While the total ESV of n-MSEA in 2010 and 2022 was 9.66 billion CNY and 20.39 billion CNY, with growth rate of 111.10%. Evidently, the ESV growth rate of MSEA was faster than that of n-MSEA. The difference rate between ESV per unit area of MSEA and that of n-MSEA in 2010 was 28.99%, while that in 2022 was declined to 15.83%. The gap between them was shortened, with low ESV area shrank obviously in MSEA (dark red in Figure 7).

Topographic gradient effect analysis shows that the ESV of the county displays an inverted U curve along elevation, while the incremental curves along slope, relief amplitude and terrain niche (Figure 8a,b). Specifically, the ESV of MSEA displays an inverted V curve along elevation with the summit in gradient II, while the incremental curves along slope, relief amplitude and terrain niche (Figure 8c,d). The ESV of n-MSEA displays an inverted V curve along elevation with the summit in gradient III, the inverted U curves along relief amplitude and terrain niche, and an incremental curve along slope (Figure 8e,f). From 2010 to 2022, the ESV of the county along elevation shows a change with the summit altered from gradient III to II, which indicates soil erosion control and ecological restoration were improved in the relative low elevation area in gradient II (385 ~ 513 m). These analyses indicate that the area with high elevation and small slope is one control focus in the next step.

Further analysis shows that major ESV was deposited in mixed forest and coniferous forest, with about 8.57 billion CNY (54.28%) and 3.50 billion CNY (22.21%) in 2010, and 19.85 billion CNY (57.13%) and 9.38 billion CNY (26.99%) in 2022 (Figure 9). It is in accordance with the dominant land uses and their changes (Figure 4 and Figure 5). Specifically, the ESV of mixed forest in MSEA was 3.04 billion CNY in 2010 and 7.14 billion CNY in 2022, with a growth rate of 134.72%; meanwhile, that in n-MSEA was 5.40 billion CNY in 2010 and 12.30 billion CNY in 2022, with a growth rate of 127.75%. The ESV of coniferous forest in MSEA was 1.38 billion CNY in 2010 and 4.05 billion CNY in 2022, with a growth rate of 194.67%; meanwhile, that in n-MSEA was 2.13 billion CNY in 2010 and 5.16 billion CNY in 2022, with a growth rate of 142.32%. Apparently, the coniferous forest in MSEA attains the fastest growth rate of ESV.

In addition, the result calculated using equation 14 indicates that the growth rate of contribution from land use changing in the ESV variation from 2010 to 2022 was 5.18% for the county, 6.62% for MSEA and 4.59% for n-MSEA, respectively. Combined with the contribution from spatial adjustment coefficient of $\sqrt{NPP\times SEVI}$, the growth rate of contribution from ecosystem resource (land uses changing coupled with NPP and SEVI changing) was 13.84% for the county, 21.41% for MSEA and 9.03% for n-MSEA, respectively.

5. Discussion

5.1. ESV of Soil Erosion Area

Under the background of “lucid waters and lush mountains are invaluable assets”, it is valuable to assess ecological restoration in the perspective of ESV for subtropical mountainous soil erosion area. In this study, the ESV of Changting county, China was estimated and analyzed. In spatial-temporal heterogeneity, the ESV per unit area of MSEA was less than that of n-MSEA, however, the ESV growth rate of MSEA was faster than that of n-MSEA from 2010 to 2022. Therefore, the difference rate between ESV per unit area of MSEA and that of n-MSEA was decreased nearly by half (from 28.99% to 15.83%). The main reasons are that the growth rate of $\sqrt{NPP\times SEVI}$ as spatial adjustment coefficient of MSEA was three times higher than that of n-MSEA (13.88% vs. 4.25%), and the growth rate of contribution from land use changing in the ESV variation of MSEA was also higher than that of n-MSEA (6.62% vs. 4.59%). These changings illustrate that vegetation coverage, land use structure, carbon sink and canopy density in MSEA were inferior to those in n-MSEA, however, the quality improvement of vegetation coverage, land use structure, carbon deposit and canopy density in MSEA was higher than those in n-MSEA. In n-MSEA, the growth rate of $\sqrt{NPP\times SEVI}$ (4.25%) and that of contribution from land use changing (4.59%) was similar; while in MSEA, that of $\sqrt{NPP\times SEVI}$ (13.88%) was more than twice as that of contribution from land use changing (6.62%). It indicates that improvement of carbon sink and vegetation canopy density are useful and still focus for soil erosion and water conservation at present and in future. In addition, topographic gradient effect analysis also illustrates that the area with high elevation and small slope is one control focus in the next step.

As a secondary succession of the ecological restoration in subtropical soil erosion area, the process is still poorly understood, since the intermingled biophysical and societal factors and drivers are often complex [34,39,40]. Theoretically, the successional pathway is from naked soil to grassland, coniferous forest, mixed forest, and evergreen broad-leaved forest in subtropical mountainous area. In this study, we observed that the dominant ESV in Changting county is from mixed forest, while the most increasing ESV is from coniferous forest, specifically, the vast expanding coniferous forest in MSEA. It is in accordance with that vegetation plays an important role in ESV supply, especially, forest provides the highest contribution [16]. These characteristics in spatial variation and ESV changing indicate that the present vegetation successional phase in Changting county is in a critical transition stage in a theoretical pathway, since both positive and negative succession [41] is still coexisted. ESV will be improved if positive succession of land uses is continued, or else, ESV will fluctuate if positive succession of land uses is disturbed. Therefore, the restoration measures combined ecological benefit with economic costs are priority in next step, e.g., closing hillsides for afforestation [24].

Finally, this study indicates that the achieved knowledge of the regional ESV growth characteristics and inherent factors is valuable and instructive for soil erosion control and ecological restoration, despite the difficulty to validate ESV in current socio-economic conditions.

5.2. Pros and Cons of the RS-MEF Method

Aiming to accurate evaluation of ESV in mountainous soil erosion area, we developed the RS-MEF method, which couples a general equivalent factor framework with the remote sensing techniques for mountains. It achieved several highlights of spatial adjustment using carbon sink and vegetation canopy, improvement spatial resolution, and removal topographic effect.

Firstly, an improved spatial adjustment coefficient using NPP coupled with vegetation index is developed. Specifically, the spatial adjustment coefficient coupled NPP with the SEVI is suitable for quantifying spatial heterogeneity of the ESV in mountains. NPP mainly represents the dry matter increment of the stem, branches, and roots of forest, while the SEVI reflects forest green leaves and grass and implies some ecological aesthetics and cultural service value. Generally, lucid waters and lush mountains are our favorite, so the greenness reflected by vegetation index complements to NPP. Using the improved adjustment coefficient, the spatial variation of the ESV in mountains is more reasonable and accurate. As for adjustment coefficient of socio-economic factors, it loses efficacy within a county, due to only a unique statistic for the Engel coefficient, percentage of urban population, or per-capita GDP in a county, respectively. Therefore, we do not consider the adjustment coefficient of socio-economic factors in the ESV assessment in this study.

Secondly, the spatial resolution of the critical parameters calculated from remote sensing images is improved. These parameters include land uses classified and the spatial adjustment coefficient of the SEVI, which derived from the used images of Landsat with 30 m resolution. Therefore, we do not depend on the external products of land use classification or vegetation indices, which resolution is 500 m or lower in general [18,42,43,44]. The RS-MEF method reduces major uncertainties from heterogeneous data, and provides better results for the regional ESV evaluation. It is especially benefit to the ESV evaluation from land use changes. For example, we detected the mixed forest was majorly changed from coniferous forest, which brings the ESV increment and is in accordance with ecological succession rule. Evidently, remote sensing can play an important role in ESV evaluation, which spatial resolution depends on the image used.

Thirdly, effective topographic correction improves the ESV evaluation accuracy in rugged mountains. Spectral characteristics variation of a land use between sunny area and topographic shadow decreases after the integrated topographic correction, which improves spectral homogeneity of the same land use in rugged mountains [30]. For example, the spectral reflectance in topographic shadow and that in sunny area are corrected to that in the adjacent flat area. Therefore, land use classification accuracy was improved, which better quantifies the ESV in mountains. In addition, a spatial adjustment coefficient of SEVI ($\sqrt{NPP\times SEVI}$) improves vegetation canopy quantification in mountains, since it has double benefits with topographic effect elimination [32] and weak target recognition in rugged green mountains [38]. They are benefit to the ESV evaluation for rugged mountains.

These improvements are effective to estimate quantitatively the ecosystem resource and the ESV in mountains. However, the disadvantages of the RS-MEF method still exist. Firstly, the standard equivalent factor calculation and the expert-based equivalent coefficients are empirical. These factors can be improved or modified locally, e.g., the equivalent coefficient of regulation of water flows for water area (102.24) can be adjusted for subtropical mountains. Secondly, spatial adjustment coefficients and land uses classification can be improved using higher resolution images, e.g., developing the NPP with 30 m resolution. Finally, the ESV validation is a great challenge, so the value transformation is urgent to study and attempt.

6. Conclusions

ESV is an essential indicator for evaluation the ecological restoration in mountainous soil erosion area. The estimated ESV of Changting county, China was increased drastically from 2010 to 2022, and displays the distinguished spatial-temporal characteristics. The ESV per unit area of MSEA in the county was less than that of n-MSEA, however, the ESV growth rate of MSEA was faster than that of n-MSEA. Especially, the growth rate of spatial adjustment coefficient ($\sqrt{\mathrm{N}\mathrm{P}\mathrm{P}\times \mathrm{S}\mathrm{E}\mathrm{V}\mathrm{I}}$) of MSEA was three times higher than that of n-MSEA. The gap of ESV per unit area between two areas was shortened about 45%. In all, the achieved knowledge of the ESV growth and inherent factors is valuable and instructive for the next step of restoration, e.g., improvement of vegetation canopy and carbon sink, and focus the area with high elevation and small slope. In addition, we proposed a RS-MEF method, which obtains several highlights in estimation of the ESV in mountainous soil erosion area, such as spatial adjustment using carbon sink and vegetation canopy, improvement spatial resolution, and removal topographic effect.