Spatiotemporal Changes of Ecosystem Service Values in Response to Land Cover Dynamics in China from 1992 to 2020

Abstract

Global land cover changed significantly in the last several decades due to strong climate warming and intensive human activities, and those changes greatly affected ecosystem services all over the world. Using CCI-LC land cover data from 1992 to 2020, the spatiotemporal characteristics of land cover change in China were investigated, and the annual ecosystem service values (ESVs) were estimated with the equivalent factor method. The results showed that: (1) The overall accuracy and Kappa coefficient of CCI-LC products in China were 71.1% and 0.65, respectively. (2) From 1992 to 2020, the area of cropland in China increased generally first before 2004 then decreased after 2008; the area of forest land decreased before 2003 then increased after 2015; the area of grassland and bare land consistently decreased; and the area of built-up land continuously increased, with a total increase of 113,000 km². The primary characteristics of land cover transitions in China were the mutual conversion of cropland, forestland, and grassland as well as the continuous increase of built-up land. (3) Forest land was the most significant contributor of ESV in China, making 62.9% of the total ESV by multi-year average, followed by grassland (18.5%) and water (10.3%); the ESV was roughly high in the southeast China and low in the northwest. (4) The total ESV in China decreased generally before 2015 and got stable in the last five years. The hot spots with rising ESV were mainly concentrated in the western, northern and southwestern parts of China, while the cold spots with declining ESV were mainly concentrated in the economically developed eastern and southern China. (5) Cropland, forest land, grassland, and water were the positive contributors to ESV change in China, while built-up land and bare land were the negative contributors. The findings provide a theoretical foundation for China’s harmonized socioeconomic and environmental development.

1. Introduction

Ecosystem services are the well-being and benefits that humans derive from ecosystems [1], which are not only the conditions and processes provided by natural ecosystems and their species to sustain human survival and meet human development needs, but also an important indicator of the quality of the regional ecological environment [2]. In 1997, Costanza et al. classified the ecosystem services into 17 categories and estimated the global ecosystem service value (ESV) for the first time [3]. In 2005, the United Nations released the Millennium Ecosystem Assessment (MEA) report [4], which classified ecosystem service functions into four categories: provisioning, supporting, regulating, and cultural services [5,6]. In 2012, the United Nations Statistical Commission established the “System of Environmental-Economic Accounting (SEEA) Core Framework” [7], which set up the principles and methods for ecosystem services accounting. The framework had been applied in ecosystem services in various regions over the world [8,9,10], and it was supplemented in 2021 [11]. To provide robust and easily accessible information on the economic benefits of ecosystems and biodiversity globally, the Foundation for Sustainable Development [12] developed the Ecosystem Services Valuation Database (ESVD), which currently contains over 6700 value records from over 950 studies distributed across all biomes, ecosystem services, and geographic regions.

Ecosystem service values can be estimated through two types of methods, i.e., (1) market valuation, which measures ESVs using market prices, such as the market price method, hedonic pricing method [13], and travel cost method [14]; and (2) non-market valuation, which uses survey-based or cost-based methods to elicit people’s willingness to pay (WTP) for an ecosystem service, such as the contingent valuation method [15], contingent choice method [16], and the replacement cost method [17]. Market valuation methods are applicable to the evaluation of some ecosystem services that can be bought and sold in commercial markets, but they are usually applicable only at local scales and neglect critical externalities and “non-use” benefits, such as human health impacts associated with common agricultural practices [18], changes in greenhouse gases storage and wild-species diversity [19], or impacts on cultural services [20]. By contrast, non-market valuation is capable of quantifying both use and non-use benefits, therefore they are more widely used for ESV estimation for regional and global studies. When it is too expensive and/or there is too little time available to conduct an original valuation study, the benefit transfer method (such as unit value transfer and benefit function transfer) is often used to estimate ESVs [21], which transfers available valuation results from pre-existing primary studies completed in another location and/or context with either market or non-market valuation methods. Costanza et al. [3] established a valuation scheme based on the benefit transfer method to estimate the global average unit values of 17 ecosystem services for 16 biomes. To make Costanza’s valuation scheme applicable to regional applications, many researchers modified its value coefficients considering local contexts, usually based on the environmental valuation database (e.g., Kindu et al. [22]). Following Costanza’s valuation scheme [3], Xie et al. [23] developed a table of the equivalent factor for per unit area of the ecosystem based on literature research, expert knowledge, statistical information, remote sensing monitoring, and other data sources through model calculations and geographic information spatial analysis. The equivalent factor method was further refined by Xie et al. [6] and widely used in Chinese applications (e.g., Gu et al. [24] and Yang et al. [25]).

Due to intensive human activities, global land cover (LC) underwent some major changes in the last several decades, including urbanization, agricultural land loss, agricultural land expansion, deforestation, afforestation, and desertification [26]. LC changes were direct drivers of ESV changes [27] and have contributed to comparable changes in the ecological environment, along with population growth and rapid economic development [28]. Many studies showed that ESV in China declined in the last several decades in response to LC changes. For instance, Wei et al. [29] studied the response of ESV to land cover change in China and concluded that the ESV in China was declining, with a total decline of 6.67% from 1990 to 2018, and that the expansion of built-up land under economic development was the most important factor contributing to the decline of ESV in China; Li et al. [30] estimated that from 2000 to 2020, the total ESV in China decreased by 0.98% (from 28.05 to 27.77 trillion yuan); Pan et al. [31] showed that from 2000 to 2018, the ESVs in China decreased by 0.71% (from 112.39 trillion yuan to 111.60 trillion yuan). Some other results showed less decreases of ESV. For instance, Song and Deng [32] investigated the ESV in China in 1988, 2000, and 2008 and found that ESVs decreased by 0.45% and 0.10% during the periods 1988–2000 and 2000–2008, respectively; Yang et al. [33] estimated that the total ESV of China fluctuated and decreased from USD 3265.3 billion in 1992 to USD 3253.3 billion in 2018, a decrease of about 0.4%.

The objective of the present study is to investigate the LC changes in China on a yearly basis during 1992 to 2020 using the latest ESA CCI-LC land cover product, estimate the spatiotemporal characteristics of the value of ecosystem as the consequence of LC changes in space and time using the equivalent factor method, and explore the relationship between ESV change and economy development.

2. Data Used and Methodology

2.1. Data Used and Pre-Processing

China’s land area is about 9.6 million km², with high mountains and plateaus in the west and low lands in the east. China has three distinguishable geographical regions, i.e., the eastern China dominated by humid/semi-humid monsoon climate, the northwestern China dominated by arid/semi-arid temperate continental climate, and southwestern China dominated by alpine climate over Tibetan Plateau. As the consequence of the complex climate, eastern China has much better vegetation coverage with vast areas of forestlands and croplands than the western parts of China, and the most human population concentrated in the eastern China. The arid/semi-arid northwestern China is dominated by grasslands and deserts.

The data used included land cover products derived from satellite images, sample data for verifying the accuracy of land cover products, geographic data for making spatial corrections to the equivalent factors, and socioeconomic data for economic value calculation of the equivalent factors.

2.1.1. Land Cover Data

European Space Agency Climate Change Initiative Land Cover (ESA CCI-LC) data are available for the period from 1992 to 2020 with a spatial resolution of 300 m (https://cds.climate.copernicus.eu/cdsapp#!/dataset/satellite-land-cover, accessed on 15 February 2022). The dataset is the longest annual LC data available to-date and is considered the first consistent time-series of global LC coverage at relatively fine spatial resolution [26]. CCI-LC data classified land cover types into 22 categories according to the land cover classification system developed by the FAO. Due to the great difference between the CCI-LC classification and the classification scheme of secondary ecosystems in the Xie’s equivalent factor table [6], for the convenience of ESV calculation, the original 22 land cover types were aggregated to six types of primary ecosystems, namely, cropland, forest land, grassland, water, built-up land, and bare land. For the purpose of analyzing the spatiotemporal variation of land cover locally, a fishnet consisting of 4258 rectangular cells of 50 km × 50 km size was created using the Fishnet tool of ArcGIS 10.2 software, then the fraction of every land cover type within each cell was calculated using the Tabulate Area tool for each year.

The FROM-GLC2010 global land cover validation point dataset, which is available at the National Earth System Science Data Center (http://www.geodata.cn/data/datadetails.html?dataguid=210340446085114&docid=32570, accessed on 14 March 2022) [34], was used to validate the accuracy of CCI-LC data, taking the year 2010 as a representative year. By visually comparing the dataset with the 2010 images in Google Earth, we removed some sample points and added some new points, so that in total, 3129 validation points of six land cover types covering China in 2010 were obtained (Figure 1).

2.1.2. Data for Ecosystem Service Valuation

Net primary productivity (NPP), precipitation, and soil erosion data were used for adjusting the equivalent factors. 2000–2010 NPP data and 1995 soil erosion data were obtained from the Resource and Environmental Science Data Center of CAS (http://www.resdc.cn, accessed on 14 March 2022) with a spatial resolution of 1 km, which were produced by the Institute of Geographical Sciences and Resources of the CAS. NPP data were calculated based on the light energy utilization model GLO_PEM; soil erosion data were compiled according to the general requirements of the industry standard of the People’s Republic of China (SL190-96). The 1992–2020 monthly precipitation data were obtained from the National Earth System Science data center (http://www.geodata.cn, accessed on 14 March 2022) with a spatial resolution of 1 km, which was produced by Peng et al. [35] from Northwest Agriculture and Forestry University.

Socio-economic data: including year-by-year data on plant area and net profit per unit area of rice, wheat, and corn in mainland China from 1992 to 2020. The plant area data are from the “China Statistical Yearbook” (1993–2021) [36] and the net profit per unit area data are from the “National Agricultural Products Cost Return Assembly” (1993–2021) [37].

2.2. Methodology

2.2.1. Accuracy Assessment of CCI-LC Product

The confusion matrix was used for land cover product accuracy assessment. The accuracy evaluation metrics derived from the confusion matrix are Overall Accuracy (OA), Producer Accuracy (PA), User Accuracy (UA), and Kappa coefficient (K). Among them, OA indicates the proportion of correctly classified area in all types; PA indicates the proportion of correctly classified area in a type to the area of that type in the data to be evaluated; UA indicates the proportion of correctly classified area in a type to the area of that type in the reference data; K is used to evaluate the accuracy and consistency of classification results. The equations for calculating these indicators are as follows:

$$\begin{array}{c}OA=\frac{{{\displaystyle \sum }}_{i=1}^r{n}_{ii}}{N}\end{array}$$

(1)

$$\begin{array}{c}P{A}_i=\frac{n_{ii}}{n_{+i}}\end{array}$$

(2)

$$\begin{array}{c}U{A}_i=\frac{n_{ii}}{n_{i+}}\end{array}$$

(3)

$$\begin{array}{c}K=\frac{N{{\displaystyle \sum }}_{i=1}^r{n}_{ii}-{{\displaystyle \sum }}_{i=1}^r\left(n_{i+}{n}_{+i}\right)}{N^2-{{\displaystyle \sum }}_{i=1}^r\left(n_{i+}{n}_{+i}\right)}\end{array}$$

(4)

where $N$ is the total number of pixels; $n_{ii}$ is the number of correctly classified pixels; $n_{i+}$ is the number of pixels of a LC type in the data to be evaluated; $n_{+i}$ is the number of pixels of the LC type in the reference data; $r$ is the number of LC classes.

2.2.2. Estimation of ESV

In this study, the equivalent factor method by Xie et al. [6] was used to estimate the ESV. The method estimating the economic value of each ecosystem service from an ecosystem as the product of an equivalence coefficient (dimensionless) and the economic value represented by its equivalence factor, which is defined as the economic value of the natural grain output per unit area of farmland based on the national average yield. Then, the total ESV is the sum of the values of all ecosystem services in the region.

(1)

Determination of equivalent factors and their economic value

Xie et al. established an ESV equivalent factor table for 9 ecosystem service types of 6 primary terrestrial ecosystems (i.e., forestland, grassland, cropland, wetland, water and desert) in 2003 [23], and proposed an improved equivalent factor table of 11 ecosystem service types for 14 secondary ecosystems of 6 primary ecosystems in 2015 [6]. As the classification of the 14 secondary ecosystems is not comparable to the 22 categories of CCI-LC data, we aggregate the 14 secondary ecosystems to 5 primary ecosystems (i.e., forestland, grassland, cropland, water, and bare land) according to the portion of each secondary ecosystem area in its primary ecosystem. The ESV equivalent factors per unit area of the five ecosystems are presented in

Table 1. ESV equivalent factors per unit area of the five ecosystems for 11 ecosystem services according to Costanza’s scheme [3] and Xie’s scheme [6].

Service Type	Ecosystem Service	Cropland	Forest Land	Grassland	Water	Bare Land
Provisioning	Food production	1/1	0.8/0.23	1.24/0.18	0.83/0.62	0/0
	Raw material production	0/0.31	2.56/0.54	0/0.26	0.31/0.32	0/0.01
	Water supply	0/−0.78	0.06/0.28	0/0.15	14.79/5.58	0/0.01
Regulating	Gas regulation	0/0.8	0/1.77	0.13/0.93	0.39/1.16	0/0.05
	Climate regulation	0/0.42	2.65/5.3	0/2.46	0/2.66	0/0.04
	Hydrological regulation	0/1.01	0.04/3.81	0.06/1.8	9.12/64.21	0/0.09
	Purify environment	0/0.12	1.61/1.57	1.61/0.81	13.48/4.36	0/0.18
Supporting	Soil conservation	0/0.72	1.90/2.16	0.56/1.13	0/1.39	0/0.06
	Maintain nutrient cycle	0/0.14	6.69/0.16	0/0.09	0/0.11	0/0
	Biodiversity conservation	0.7/0.15	0.33/1.96	0.89/1.03	0/4.42	0/0.06
Cultural	Aesthetic landscape	0/0.07	1.26/0.86	0.04/0.46	2.99/2.85	0/0.02
Total		1.7/3.97	17.9/18.64	4.53/9.3	41.92/87.68	0/0.52

Note: The value after “/” in the table was proposed by Xie et al. [6] for China in 2015. The equivalent factor before “/” was calculated based on the results of Costanza et al. [3] at a global scale, by taking the value of food production on cropland of 54 $/hm²/year in 1994 as the economic value of one equivalent factor and dividing the values per unit area of different service functions of all different ecosystems by 54. The values of wetlands, lakes/rivers, and ice are merged into water according to the area proportion of the three types in water.

. The 8 terrestrial ecosystems (i.e., forest, grassland, cropland, wetlands, lakes/rivers, desert, tundra, ice/rock) are also regrouped into 5 primary ecosystems and as the ESV of built-up land vary greatly from negative values to positive values depending on the specific land use and management practices, which makes it very hard to give reliable estimates; it is common to ignore the ESV of built-up land [3,6]. Therefore, we set the ESV of built-up land to be 0 in the present study.

When Costanza et al. (1997) [3] estimated global total ESV, they estimated unit area values for the ecosystem services and then multiplied by the total area of each biome. This can only be considered a crude first approximation and can introduce errors depending on the type of ecosystem service and its spatial heterogeneity. To reflect the spatial differences in ecosystem service functions, based on related studies, three regulating factors—NPP, precipitation, and soil erosion—are often used to spatially correct the original equivalence factors [38,39], which were also adopted in this study. The correction equations for the equivalent factor within each grid are as follows:

$$\begin{array}{c}{V}_{kf}^i=\left\{\begin{array}{c}{P}_k^i\times V_{kf1} or\\ R_k^i\times V_{kf2} or\\ S_k^i\times V_{kf3} \end{array}\right.\end{array}$$

(5)

where $V_{kf}^i$ refers to the equivalent factor of class f service function (f = 1, 2,..., 11) of type k ecosystem (k = 1,2,...,5) in grid i; $V_{kf}$ refers to the national average equivalent factor of class f service function of the type k ecosystem (f₁ represents the service functions associated with NPP, namely food production, raw material production, gas regulation, climate regulation, maintain nutrient cycle, purify environment, biodiversity conservation, and aesthetic landscape; f₂ represents service functions related to precipitation, namely water supply and hydrological regulation; f₃ represents a service function related to soil erosion, namely soil conservation); $P_k^i$, $R_k^i$, and $S_k^i$, respectively, refer to the NPP regulatory factors, precipitation regulatory factors, and soil erosion regulatory factors of type k ecosystems in grid I. These three regulatory factors correspond to the correction of service functions of f₁, f₂, and f₃, respectively.

The NPP regulatory factor was calculated as follows:

$$\begin{array}{c}{P}_k^i=\frac{B_i}{B}\end{array}$$

(6)

where $B_i$ refers to the annual average NPP of grid i; $B$ refers to the national annual average NPP.

The rainfall regulatory factor was calculated as follows:

$$\begin{array}{c}{R}_k^i=\frac{C_i}{C}\end{array}$$

(7)

where $C_i$ refers to the annual average rainfall of grid i; $C$ refers to the national annual average rainfall.

The soil erosion regulatory factor was calculated as follows:

$$\begin{array}{c}{S}_k^i=\frac{E}{E_i}\end{array}$$

(8)

where $E_i$ refers to the average soil erosion intensity of grid i; $E$ refers to the national average soil erosion intensity.

For the measurement of the economic value of the equivalent factor, this study used the net profit of food production per unit area of the cropland ecosystem (considering rice, wheat, and maize) as the economic value of one standard equivalent factor, calculated as follows:

$$\begin{array}{c}D=S_r\times F_r+S_w\times F_w+S_c\times F_c\end{array}$$

(9)

where D is the economic value of a standard unit of ESV equivalent factor; $S_r$, $S_w$, and $S_c$ represent the percentage (%) of the planting area of rice, wheat, and maize to the total planting area of the three crops in that year, respectively; $F_r$, $F_w$, and $F_c$ represent the average net profit per unit area (CNY/hm²) of rice, wheat, and maize nationwide in that year, respectively.

(2)

Calculation of regional total ESV

After specifying the equivalent factors and their economic values, the total ESV can be calculated based on the area of each land cover type as follows:

$$ESV={\displaystyle \sum }_{i=1}^l{\displaystyle \sum }_{k=1}^m{\displaystyle \sum }_{f=1}^n(V{C}_{kf}^i\times A_k^i)$$

(10)

$$\begin{array}{c}V{C}_{kf}^i=D\times V_{kf}^i\end{array}$$

(11)

where ESV is the total value of ecosystem services in China (CNY); $A_k^i$ is the area of the kth type of the land ecosystem in grid i; $V{C}_{kf}^i$ is the value coefficient of the type f function of the type k land ecosystem in grid i; l is the number of grids, which refers to 4258 grids in the study; m is the number of ecosystem types, with 5 categories in the study; n is the number of ecological service function categories, which is 11 categories in the study; D is the economic value of a standard equivalent factor, calculated using Equation (9).

2.2.3. Analysis of Spatiotemporal Changes

(1)

Transition analysis of land cover types

The transition matrix is used to reflect the spatial changes of LC types in China from 1992 to 2020. A transition probability matrix (A) describes the probability of change from a land cover type to another between two LC maps of the different time (year t and the subsequent year t + c) as follows [40]:

$$X_{t+c}=X_{t}A$$

(12)

$$A=\left[a_{ij}\right]=\left[\begin{array}{ccc}{a}_{11}& \cdots & a_{1n}\\ ⋮& \ddots & ⋮\\ a_{n1}& \cdots & a_{\mathrm{nn}}\end{array}\right]$$

(13)

where $X_t$ is an m × n matrix representing m pixels with n land cover types at time t. A is an n × n matrix with each element $a_{ij}$ representing the conditional transition probability of a pixel from LC category i at time t to change into LC category j by time t + c.

(2)

Trend change and hotspot analysis

The Mann–Kendall (MK) nonparametric statistical method is used for testing the trend of time series. The method is not affected by outliers and does not require the data to follow a specific distribution. The null hypothesis for the MK test is that there is no monotonic trend in the series. We can judge the strength of the trend based on the correlative change rate τ and its p-value of the MK test results, which are presented in the form shown in

Table 2. Trend test p-value and its corresponding trend category.

$\mathsf{\tau }$	p-Value	Trend Category	Trend Characteristics
$\mathsf{\tau }$ > 0	p ≤ 0.01	3	Extremely significant increase
	0.01 < p ≤ 0.05	2	Significant increase
	0.05 < p ≤ 0.1	1	Slight increase
$\mathsf{\tau }\approx 0$	p > 0.1	0	No change
$\mathsf{\tau }$ < 0	0.05 < p ≤ 0.1	−1	Slight decrease
	0.01 < p ≤ 0.05	−2	Significant decrease
	p ≤ 0.01	−3	Extremely significant decrease

[41].

The trend of ESV changes is tested over China at the grid scale, and the hotspots and cold spots analyses were further used to identify the spatial clustering of high values (hot spots) and low values (colds pots) to reflect the clustering effect of ESV change trends in the grid cells in China. As a tool integrated in ArcGIS 10.2, this approach takes each raster pixel within the context of neighboring features into the calculation and outputs a new feature class with z-score, p-value, and confidence level. Features with high z-score and small p-value indicate statistically significant hotspots, and features with low negative z-score and small p-value demonstrate statistically significant cold spots. The magnitude of the absolute value of the z-score explains the intensity of the clustering [42]. The principle of this method is shown as follows [43]:

$$\begin{array}{c}{G}_i^{*}=\frac{{{\displaystyle \sum }}_{j=1}^n{w}_{i,j}{x}_j-X{{\displaystyle \sum }}_{j=1}^n{w}_{i,j}}{\sqrt[s]{\frac{[n{{\displaystyle \sum }}_{j=1}^n{w}_{i,j}^2-{\left({{\displaystyle \sum }}_{j=1}^n{w}_{i,j}\right)}^2]}{\left(n-1\right)}}}\end{array}$$

(14)

$$X=\frac{1}{n}{\displaystyle \sum }_{j=1}^n{x}_i$$

(15)

$$S=\sqrt{\frac{1}{n}{\displaystyle \sum }_{j=1}^n{x}_j^2-{\left(X\right)}^2}$$

(16)

where $G_i^{*}$ is a z-score of patch i; X is the mean of the element values; S is the standard deviation of the element values; n is the total number of grid cells; $x_i$ and $x_j$ are the attribute values for patch i and j, respectively; $w_{i,j}$ is the spatial weight between patch i and patch j, if the distance from a neighbor j to the feature i is within the distance, $w_{i,j}$ = 1; otherwise, $w_{i,j}$ = 0. Identifying and mapping the hotspots and cold spots can visualize priority areas spatial-explicitly, which is helpful for targeted policy making.

(3)

Bivariate spatial autocorrelation analysis

The bivariate Moran’s I, including the global bivariate spatial autocorrelation and local bivariate spatial autocorrelation, is used to analyze the spatial agglomeration and spatial dispersion pattern between the area change of LC types and ESV change across China in this study. We used the $\mathsf{\tau }$ (Table 2) of LC area change and the $\mathsf{\tau }$ of ESV change as two variables. The global bivariate Moran’s I was used to detect the spatial autocorrelation between the two variables in the whole region and the calculation Equations were as follows [44]:

$$\begin{array}{c}I=\frac{{{\displaystyle \sum }}_{i=1}^n{{\displaystyle \sum }}_{j=1}^n{w}_{ij}\left(x_i-\overline{x}\right)\left(y_j-\overline{y}\right)}{{{\displaystyle \sum }}_{i=1}^n{{\displaystyle \sum }}_{j=1}^n{w}_{ij}\sqrt{S_x^2{S}_y^2}}\end{array}$$

(17)

where $I$ is the global bivariate Moran’s I; $n$ is the number of spatial units (grids); $w_{ij}$ is the element of the normalized neighborhood matrix; $x_i$ and $y_j$ are the attributed values of independent variable x and dependent variable y in space units i and j, respectively; $\overline{x}$ and $\overline{y}$ are the mean values of independent variable x and dependent variable y, respectively; $S_x^2$ and $S_y^2$ are the variances of independent variable x and dependent variable y, respectively.

The local bivariate Moran’s I was used to detect the aggregation and differentiation of local spatial elements given by [44]:

$$\begin{array}{c}{I}_i=\frac{x_i-\overline{x}}{S_x^2}\end{array}{{\displaystyle \sum }}_{j=1}^n{w}_{ij}\frac{y_j-\overline{y}}{S_y^2}$$

(18)

where $I_i$ is the local bivariate Moran’s I. When the value of $I_i$ is greater than 0, it indicates that there is a positive spatial correlation between these two variables; when $I_i$ is less than 0, it indicates that there is a negative spatial correlation between these two variables. The $I_i$ can be used in conjunction with the Moran scatter plot to screen out the spatial units that show significant spatial dependence through the z test (p < 0.01). Moran scatter plot refers to the scatter plot which takes the observed value ($\mathsf{\tau }$ of land cover types) of each space unit as the abscissa and its spatial lagged term ($\mathsf{\tau }$ of ESV) as the ordinate. The spatial lagged term refers to the weighted mean of another variable around the observed value, which in this study reflects to the degree of ESV change in the surrounding area, with the change of a certain land cover type at a certain location. The four quadrants of the scatter plot express four specific patterns of spatial autocorrelation between a certain point (region) and surrounding points (region), namely, “high-high” (first quadrant), “low-high” (second quadrant), “low-low” (third quadrant), and “low-high” (fourth quadrant).The four patterns are represented on the map to make the LISA cluster map [45].

3. Results

3.1. Assessment of CCI-LC Land Cover Product Accuracy

The accuracy of CCI-LC data in 2010 was assessed using the ground-truth samples in 2010. The Chinese land cover map based on CCI-LC data in 2010 is shown in Figure 2. The confusion matrix for CCI-LC is presented in

Table 3. CCI-LC classification accuracy in China in 2010.

	Validation Sample Type
LC Type	Cropland	Forest Land	Grassland	Water	Built-Up Land	Bare Land	Total	User Accuracy (%)
Cropland	478	159	48	13	36	66	800	59.8
Forest land	25	471	21	1	2	10	530	88.9
Grassland	37	67	295	40	4	321	764	38.6
Water	1	0	1	142	0	5	149	95.3
Built-up land	5	1	1	7	411	0	425	96.7
Bare land	0	5	23	4	1	428	461	92.8
Total	546	703	389	207	454	830	3129
Producer accuracy (%)	87.6	67.0	75.8	68.6	90.5	51.6

, which shows that the overall accuracy of the CCI-LC land cover product in China was 71.1% with a kappa coefficient of 0.65. The best accuracy was achieved for built-up land, with producer accuracy of 90.5% and user accuracy of 96.7%, indicating that most of the built-up land was correctly identified with rare misclassification and omission. The producer accuracy and user accuracy of the cropland were 87.6% and 59.8%, respectively, indicating that there was misclassification of this type, and omission was relatively rare, which was due to the similar spectral characteristics of cropland and forest land, and it was easy to misclassify forest land as cropland. The producer accuracy and user accuracy of the forest land type were 67% and 88.9%, respectively, indicating that there was a phenomenon of omission for this type and the phenomenon of misjudgment was relatively rare, which was due to the fact that forest land is mostly distributed in undulating terrain and highly mixed with various vegetation types, such as grassland and cropland, so the automatic identification method was more difficult to discern the woodland type and there will be an omission. The producer accuracy and user accuracy of the grassland type were 75.8% and 38.6%, respectively, which indicated that the misclassification of this type was serious and the omission was rare. The producer accuracy and user accuracy of the bare land type were 51.6% and 92.8%, respectively, which indicated that the misclassification of this type exists and the omission was rare, which was mainly caused by the misclassification of bare land as grassland. The producer accuracy and user accuracy of the water type were 68.6% and 95.3%, respectively, which indicated that the omission of this type existed and the misclassification was rare, which was due to the error caused by the inconsistency between the temporal phase of the classification data and the manual identification and verification data for water bodies, such as rivers and lakes, which were influenced by the season. The accuracy of CCI-LC products has also been evaluated by scholars, and the results of Yang et al. [46] showed that the overall accuracy of CCI-LC products was 72.0% with a Kappa coefficient of 0.65. In general, the assessment results were more consistent with those of previous authors.

3.2. Spatiotemporal Characteristics of Land Cover Change

3.2.1. Temporal Change of LC Areas

Figure 3 shows the area changes of six LC types in China from 1992 to 2020. The area of cropland increased significantly before 2000, became stable in the next 10 years, then decreased after 2010, with an overall increase of 3484 km². The area of forest land exhibited a trend opposite to that of cropland, decreased before 2003, then increased after 2015, with a total decrease of 8088 km². The area of grassland and bare land both decreased continuously, with a total decrease of 38,000 km² and 73,000 km², respectively. The area of water changed very little before 2012, but increased significantly afterwards, with a total increase of 2488 km². The area of built-up land kept increasing steadily, with the highest changing rate during the study period among all different types, with a total increase of 113,000 km² and an average annual growth rate of 11.5%. It can be seen that although the total amount of cropland and forest land did not change much, both experienced large ups and downs during the study period, with the area of arable land increasing sharply from 1994 to 2000 and the area of forest land and grassland decreasing sharply during this period, due to the fact that the continuous increase of the Chinese population led to an increased demand for grain consumption [47] before 2000, and a large area of forest land and grassland was reclaimed. In 2000, China’s Grain for Green Program was carried out nation-wide after it being tested in three provinces in northern China [48]. The area of cropland remained stable and started to decline after 2000, while the area of forest land started to rebound after 2003, which is strongly related to the Grain for Green Program. With the fast economic boom in the last four decades, China experienced rapid urbanization, leading to a continuous expansion of built-up land. In general, the change in areas of various LC types in China is characterized by an early increase and late decrease in cropland, an early decrease and a later increase in forest land, significant decreases in bare land and grassland, a fluctuating increase in water, and a continuous expansion of built-up land. The results are consistent with those of Liu et al. [47].

3.2.2. Characteristics of LC transitions

(1)

Temporal variation of LC transitions

The Sankey diagram was produced based on the LC type transition matrix for the three time periods 1992–2000, 2000–2010, and 2010–2020 (Figure 4). Only the converted proportions are shown because the unconverted area of each LC type accounted for a large proportion.

Combining the annual loss and gain of each category (Figure 5) from 1992 to 2000, the gain of croplands was more than its loss, and the main sources of the gain of croplands were forestland and grassland. Forest lands maintained a high amount of losses, which was mainly converted to croplands and grasslands, and the main sources of gain were grassland. The amount of losses and gains of grasslands are highly positively correlated, with most losses being cropland, forest land, and bare land, and the main sources of gains were forest land and bare land.

From 2000 to 2010, the loss and gain of cropland was relatively similar, which was mainly converted to built-up land, and the main sources of gain were forest land, grassland, and bare land; forest land was mainly converted to cropland and grassland, and the main source of gain was grassland; grassland was mainly converted to cropland, forest land, built-up land, and bare land, and the main sources of gain were forest land and bare land; bare land was mainly converted to grassland.

From 2010 to 2020, cropland maintained a high volume of loss, which was mainly converted to forest land, grassland, and built-up land, and the main sources of gain were forest land and grassland; forest land maintained a high volume of gain, and the main sources of loss and gain were both cropland and grassland; grassland maintained a high volume of loss and gain, which was mainly converted to cropland, forestland, built-up land, and bare land, and the main sources of gain were cropland, forestland, and bare land; the vast majority of loss of bare was converted to grassland.

In general, the mutual conversion of cropland, forest land, grassland, and the continuous gain of built-up land are the main features of transition of LC types in China. The total area of LC change between 1992 and 2020 in China was 824,013 km², equal to 8.7% of the total area, and the largest transition was bare land converting to grassland, accounting for 15.7% of LC transitions in China, followed by the conversion from grassland to forest land, accounting for 12.1% of LC transitions.

(2)

Spatial characteristics of LC transitions

In this study, the MK trend test was used to judge the variation trend of variables in each grid in time series. Figure 6 shows the results of the MK trend test about the proportion of every LC types in each 50 km × 50 km grid cells over the period from 1992 to 2020. According to Figure 6, it can be seen that the regions with significant increase in cropland area were mainly concentrated in the southeast (Hunan, Jiangxi, Fujian, and northern Guangdong provinces), northeast (northeastern part of Inner Mongolia), and northwest (western and northern part of Xinjiang) regions, and the regions with significant decrease were mainly concentrated in the central, eastern coastal provinces (Jiangsu, Shandong, Hebei, and Liaoning), and some southern regions (southern Tibet, Yunnan, Guizhou, eastern Sichuan, Guangxi, southern Guangdong, and Hainan). The majority of the regions with decreasing (increasing) forest land area correspond to those with increasing (decreasing) cropland area in eastern China, while in Xinjiang, forestlands increased together with cropland. The grassland area exhibited a large-scale decrease in the southeastern half of China due to the significant increase of either croplands or forelands, while the grassland area in the arid/semi-arid regions in northwestern China (such as Gansu, central Inner Mongolia, southwestern Tibet, northern Qinghai) increased significantly due to joint effects of the Grain to Green Program and the climate wetting and warming in those areas [49]. While the overall change in the area of water was not significant, there was a significant change in some local areas, such as the Qinghai-Tibet Plateau due to the increase in glacial meltwater and precipitation caused by global warming [50], while the area of water in the southeast has decreased significantly, mainly due to human activities, such as the large amount of irrigation diversions and the construction of upstream reservoirs, which have led to a decrease in the downstream water supply. There was a large and highly significant increase in the area of built-up land in the southeastern half of China (from Heilongjiang to Yunnan) where the economy developed very fast and the population density is very high, whereas the northwestern half experienced little increase in built-up lands except for scattered locations. The bare lands are mostly distributed in the northwestern half of China, and their area significantly decreased (or increased), corresponding to the increase (or decrease) of grasslands increasing (or decreased), possibly due to the very significant climate change in this part of China. As shown in Figure 4, there were large amounts of a mutual conversion between the bare land and grassland largest transition, and by visually inspecting Figure 6, we know such kind of conversions occurred mostly in the northwestern half of China.

3.3. Spatiotemporal Change Characteristics of ESV

3.3.1. Temporal Change Characteristics of ESV

The first thing to do for calculating ESV is to estimate the economic value of a standard equivalent factor, D, using Equation (9). Xie et al. [6] proposed to set D equal to 1/7 of the grain yield per unit area multiplied by the price per kilogram in that year. However, the grain price may vary greatly due to market fluctuations. Therefore, we took the average value of grain net profit from 1992 to 2020 as D in our study. Using the data of planting area of rice, wheat, and maize retrieved from the China Statistical Yearbook issued by the National Bureau of Statistics of the People’s Republic of China, and the average net profit per unit area (CNY/hm²) of rice, wheat, and maize retrieved from the National Agricultural Products Cost Return Assembly [37] issued by the Department of Price of China National Development and Reform Commission (NDRC) during 1993–2021, the value of D was estimated to be CNY 1573.6/(hm²·year) with Equation (9).

The changes in the ESVs for different LC types in China were obtained based on the estimation method of ESV (Figure 7). The total ESV in China between 1992 and 2020 showed a variation of first decreasing, then rebounding, and getting stable in the last five years. The significant decrease in the total ESV in 1994–1995 and 1998–1999 was mainly related to the sharp decrease in forest land area, while the significant rebound in 2015–2016 was also related to the significant increase in forest land area, which indicates that the change in forest land area is the dominant factor in the change of the total ESV in China. In general, the ecological environment in China was in a state of continuous deterioration before 2015, which is consistent with the findings of Song et al. [32] and Li et al. [28]. The total ESV decreased by CNY 307.93 billion during the whole study period. As expected, the variations of ESVs of cropland, forest land, grassland, water, and bare land have almost the same temporal pattern as the variations of those LC types. The ESVs of cropland increased sharply from 1994 to 2000 and decreased after 2005, with an overall decrease of CNY 4.61 billion. Forest land is the most important part of China’s ESV system, making more than 62% of the total ESV in China. Over the study period, the ESV of forest land largely decreased by CNY 69.07 billion before 2015 but increased in the last five years. Grassland is also an important part of China’s ESV system. The value of grassland kept decreasing during the study period, with a total reduction of CNY 208.52 billion. During the study period, the value of water exhibited a decreasing trend before 2012, then increased afterwards, with an overall decrease of CNY 23.98 billion. The contribution proportion of bare land is the smallest and the ESV of which showed a continuous decrease trend. Overall, the decrease in ESVs indicates ecological degradation due to LC dynamics, but ecological degradation stopped after 2013, and improvements were clearly seen during 2015~2020, mainly due to the increase of forest land area.

3.3.2. Spatial Change Characteristics of ESV

(1)

Spatial Distribution of ESVs in China

To reflect the spatial differences in ecosystem service functions, we use three regulating factors, i.e., NPP, precipitation, and soil erosion, to spatially correct the equivalence factors for ESV calculation. Figure 8 shows the spatial distribution of the regulators at the grid scale in China. It can be seen that the spatial distribution of these three regulators was highly similar, all decreasing from the southeastern coastal area to the northwestern inland areas, indicating that these three regulators had a strong relationship and interaction.

Considering the area of China and the scale effect of the minimum plasticity unit on the measured results. After many tests, we chose the 50 km × 50 km grid as the basic unit and calculated the mean ESV of each 50 km × 50 km grid over China in the period from 1992 to 2020. The 50 km × 50 km grid-based average ESV over China is presented in Figure 9. In general, the spatial distribution of ESV in China varied significantly, basically consistent with the spatial distribution of NPP and rainfall regulating factors, roughly high in the southeast and low in the northwest, which was consistent with the findings of Yang et al. [33] and Xue et al. [39]. The high-value areas were mainly distributed in the southern region, especially in Fujian, Taiwan and the southwestern part of Yunnan provinces, where rainfall is abundant and solar radiation is sufficient for vegetation growth. The low-value areas were mainly distributed in the northwest, including Tibet, Xinjiang, Qinghai, Gansu, Ningxia, northern part of Shaanxi, and most of Inner Mongolia provinces, where precipitation is low but evaporation is high, and forest and grass vegetation coverage is low.

(2)

Trends and hot spots of ESV changes

To analyze the temporal trend of ESV at the grid scale, and show the clustering pattern of the trend, we conducted the MK trend analysis for each 50 km × 50 km grid and hotspot analysis for the trend analysis results. Figure 10 shows the MK trend analysis results and the distribution of hot and cold spots of the ESV change trend from 1992 to 2020. It can be seen that there are large ESV increasing hotspots in northwest China, primarily in Xinjiang, Gansu, northern Qinghai, and the western part of Inner Mongolia provinces. The increase of ESV in these regions was mostly strongly tied to the increase of cropland, forest land, and grassland (Figure 6). There were significant hot spot gathering areas of ESV in Yunnan, northern Heilongjiang, and some central areas (Shaanxi, Hubei, Chongqing, and Sichuan borders), where a significant increase in forest land area (Figure 6) had led to an increase in ESV, while an increase in water area (Figure 6) was the main reason for the increase in ESV in central western Tibet. Cold spots prevail in the plains (Northeast Plain, North China Plain, Middle and Lower Yangtze River Plain) and hilly areas (Liaodong Hills and Southeast Hills) in the eastern half of China, where rapid economic development and rapid expansion of built-up land had occupied a large amount of cropland and forest land, leading to a large decrease in ESV. In addition, there were also obvious cold spot gathering areas at the central eastern Tibetan Plateau, mainly due to the LC transition from grassland to cropland.

3.4. Spatial Autocorrelation between Land Cover Change and ESV Change

The spatial correlation characteristics of the trends of LC area change and ESV change across China were analyzed using the bivariate spatial autocorrelation model. $\mathsf{\tau }$ value (Table 2) in trend analysis was used to represent the variation trends of LC area change and ESV change. The bivariate global Moran’s I values were counted based on the obtained Moran scatter plots (Figure 11) and significance tests were performed to obtain z scores and p values. After calculation, all p values were 0.001 and the absolute value of z was larger, indicating that the null hypothesis can be rejected when the confidence level is 99%, and the results of the confidence test are highly significant. Global Moran’s I results show that the temporal variations of cropland, forest land, grassland, and water have significant positive spatial correlations with ESV change trends, among which the positive correlation between forest land and ESV is the strongest. Significant negative spatial correlations present between the trends in ESV and the trends in the area of built-up land and bare land.

Bivariate local spatial autocorrelation LISA clustering plots for ESV and the area variation of each LC type are presented in Figure 12. It can be seen that the local spatial association between the trends of ESV and the trends of cropland, forest land, grassland, and water, which have positive global spatial correlation with ESV, exhibit similar spatial pattern. The high-high aggregation areas of these four LC types changes and ESV changes were mainly distributed in the northwestern China and central Inner Mongolia. Some large areas of high-high aggregation of forest land and water changes and ESV changes are also found in Yunnan Province. Large areas of low-low aggregation of cropland, forest, grassland changes, and ESV changes present in the middle and lower reaches of the Yangtze River Plain, North China Plain, and Northeast China Plain, whereas large areas of low-low aggregation of forest land changes and ESV changes present in the hilly areas of Southeast China, and low-low aggregation of grassland changes and ESV changes in the eastern part of the Qinghai-Tibet Plateau (at the junction of Tibet, Qinghai, and Sichuan provinces). The spatial association between the area changes of built-up land, bare land, and ESV changes was highly similar, with the low-high aggregation areas of both two LC types changes and ESV changes distributed in northwest China and central Inner Mongolia, also in central and southwest Tibet. The high-low aggregation areas of land changes and ESV changes mainly occurred in the southeastern hilly areas, the middle and lower reaches of the Yangtze River Plain, and the North China Plain, but there were also large areas of high-low aggregation areas of built-up land change and ESV change in the Northeast Plain.

4. Discussion

4.1. Uncertainty Analysis of ESV Estimation Results

Since the equivalent factor method has great advantages in terms of convenience in data acquisition and calculation [3,51,52], it is used in the present study to estimate ESV. When using the equivalent factor method, the reliability of ESV estimation is controlled by four factors: the set-up of the equivalent factor table, the estimate of the economic value (D) of equivalent factors, the spatial correction method, and the accuracy of LC classification.

(1)

Estimate of equivalent factors

Although the valuation of Xie’s equivalent factor [6] is still controversial [53], it has been widely used at national [33], provincial [54], municipal [55], and topographic area [56] scales in the Chinese context due to its advantages in terms of estimation cost, estimation time, difficulty of data acquisition, and operability [52]. Therefore, in this study, Xie’s equivalent factor table was adopted. We calculated the equivalent factor table based on Costanza’s ESV valuation scheme as well, which is also presented in Table 1 for comparison. We can find large difference between the two equivalent factor tables due to the fact that both equivalent factor tables are highly subjective, based on the scores given by respondents to the assessment of willingness to pay for various types of ecosystem services, and experts in different socio-geographic contexts and at different times had different definitions and valuations of ecosystem services. On the one hand, values of some ecosystem services of a few LC types were ignored in Costanza’s equivalent factor table, such as most ESVs of cropland, the gas regulation service values of forest land, and all ESVs of bare land; on the other hand, although the values of some ecosystem services of some LC types were estimated, the differences were very large. For example, the ESV of water’s hydrological regulation function by Xie’s table is about seven times that of Costanza’s estimate. In terms of the aggregated equivalent values, Xie’s equivalent factors tend to overestimate total ESV compared to Costanza’s equivalent factors.

One special item for ESV estimation is the value of built-up lands for ecosystem services, which is complex and may vary in different regions. Some quantitative studies were conducted on the impacts of built-up land on the ecosystem. Duan et al. [57] concluded that the ESV of built-up land in Beijing was negative by the willingness-to-pay method; Wratten et al. [58] concluded that there was a degree of a positive role of built-up land in climate regulation, water regulation, and recreation functions; Li et al. [59] found that the negative ESV by built-up land landscape in Xi’an was much higher than the positive ESV. In general, built-up land had a negative effect on the ecosystem, but there is a lack of consensus about the exact ESV of build-up land. Taking the study of Duan et al. [60] as an example, the equivalent factor of built-up land is approximately −6.07, then roughly the ESV was overestimated by about CNY 33.48 billion in 1992 and by about CNY 141.39 billion in 2020. In this study, the equivalent factor of built-up land was treated as 0 and the negative impact of built-up land on the ecological environment was ignored, which may result in an overestimate of total ESV.

(2)

Estimate of the economic value D of the equivalent factor

In this study, we took the net profit of grain production per unit area of the cropland ecosystem as the D value (Equation (9)) following the method of Xie et al. [6], while Costanza et al. [3] derived the value of food production per unit area of cropland as 54 CNY/hm²/year using benefit transfer method based on a large number of case studies worldwide. If we adopted Costanza’s estimate (CNY 54 in 1994) instead to set the D value, D would equal 551.6 CNY/hm²/year, by adjusting $54 according to the CPI of USA in 1994 and 2010 (https://stats.areppim.com/calc/calc_usdlrxdeflxcpi.php, accessed on 2 November 2022), and converting USD to CNY using the average exchange rate in 2010, which was much lower than 1573.6 CNY/hm²/year in this study. Compared the D values with others, our D value is similar to those of Yang et al. [33] (1847.46 CNY/hm²/year) and Chen et al. [61] (1231.9 CNY/hm²/year), but less than Xie et al. [6] (3406.5 CNY/hm²/year). The different setting of D value greatly impacts the estimation result of ESV (see

Table 4. Comparison with related studies.

Area	Year	ESV (Trillion CNY)	ESV/GDP	D Value (CNY/(hm2·Year))
Global	1997	120.84	0.37	551.60	Costanza et al. [3]
	2011	533.81	1.00	16,511.83	Costanza et al. [51]
China	2010	38.11	0.87	3406.5	Xie et al. [6]
		22.02	0.5	1847.46	Yang et al. [33]
		110.57	2.52	/	Pan et al. [31]
		28.05	0.64	/	Li et al. [30]
		16.13	0.37	1231.9	Chen et al. [61]
		18.76	0.43	1573.6	This study

Note: The values in Costanza et al.’s study were estimated for 1997 [3] and 2011 [51], respectively, while the ESVs for 2010 in other studies were adopted here for comparison. Only ESVs of terrestrial ecosystems in their studies were considered here for comparing with studies in China, which were focused on terrestrial ecosystems. According to the exchange rate of USD to CNY in 2010 (6.77), the USD in Costanza et al. [51] and Yang et al. [33]’s paper was converted into CNY.

for details).

(3)

Spatial correction of equivalent factors

To represent the spatial heterogeneity of ESV, the equivalent factors need to be spatially corrected. Xie et al. [6] showed that ecosystem food production, raw material production, gas regulation, climate regulation, environmental purification, nutrient cycle maintenance, biodiversity conservation, and aesthetic landscape functions were positively correlated with the biomass in general, and water supply and hydrological regulation were positively correlated with precipitation variation. Xue et al. [39] and Sun et al. [62] showed that the soil conservation function was negatively correlated with soil erosion degree. According to their results and considering the availability of data and the feasibility of calculation, three factors, i.e., NPP, precipitation, and soil erosion, were selected for conducting spatial correction at the grid scale in this study. However, the factors influencing the ESVs are multiple and complex. In addition to the above three factors, many other factors affect the spatial pattern of ESV. For instance, habitat quality [63] and accessibility [64] affect the biodiversity and aesthetic landscape functions of ecosystems respectively, and specific ecosystem service functions are not influenced by only a single factor; Zhou et al. [65] considered the normalized vegetation index (NDVI), leaf area index (LAI), and net primary productivity (NPP) to regulate the original equivalent factors; Dong et al. [66] calculated that the soil conservation function of the forest was affected by multiple factors, such as rainfall, NDVI, and soil texture. Moreover, people have different preferences in different countries and regions with different levels of economic development and residential income, which will lead to great differences in people’s willingness to pay for ecosystem services [67], so the socio-economic factors such as population and GDP [68,69] will influence the estimated results from the aspect of economic value. Therefore, appropriate selection of correction factors needs to be further explored. At the grid scale, multiple regulatory factors can be selected to focus on the ecosystem service function, and not limited to a single factor to modify a service function, but can be modified according to the influence degree of different factors in proportion; while under the administrative unit, correction factors can be selected focusing on socio-economic aspects to improve the accuracy of the estimation results.

(4)

Accuracy of land cover classification

The regional ESV is a summation of the ESV of each individual ecosystem, which depends on the classification of LC types, so the scheme of LC classification and the accuracy of the LC data products used also have a large impact on ESV estimation. Errors in LC products will be transferred to the accounting of ESV. According to the accuracy assessment results of LC products in this study, the user accuracy of cropland and grassland and the producer accuracy of forest land and water are not high; and there was also the phenomenon of omission of water, thus making the area of forest land and water with high ecological service value types small compared to the actual, which led to the underestimation of ESV. In addition, the ecosystem classification system in Xie’s equivalent factor table is not completely consistent with the classification system of CCI-LC land cover data. Fortunately, LC data are constantly updated and improved. In this study, long-term LC data with the highest temporal resolution were adopted, with a spatial resolution of 300 m, which was superior to the existing research in terms of temporal resolution.

(5)

Temporal variations of equivalent factors and D values

The temporal variations of equivalent factors and D values make the valuation of ecosystem services more complicated. ESV may vary over time due to changes in the ecological system, societal needs, and economic conditions. For example, as the aquatic ecological environment has improved in China in recent years [70], the service function of water should be improved; the value of a wetland ecosystem may increase over time if it is used for recreational purposes as the local population grows; with the increase of incomes and environmental awareness, people would have a higher willingness to pay for ESs (e.g., Kalfas et al. [71]). In those cases, we need to adjust equivalent factors of those ecosystem services as the aquatic environment and social need change. Moreover, the D value is subject to changes in the grain market. According to the line fitted in Figure 13, the average price of China’s three major agricultural products increases by about CNY 0.06 per year. As a result, even if the condition of ecosystems did not change, the total ESV would vary with the changes in the grain value. Shi et al. [72] introduced the concept of comparable prices and comparable growth rates from economic sciences and used the CPI index to uniformly adjust the prices in different years to the constant prices in 1999. In this study, the multi-year average was adopted to eliminate the influence of price variation, without taking into account the dynamics of ESV over time. Such a practice can eliminate market influences but ignore the temporal variation of economic values itself. One possible solution is to take the moving average over several years (e.g., 5 years or 9 years), which is worth further investigation in the future.

4.2. The Relationship between ESV and Economic Development

With the continuous growth of the population, the continuous expansion of built-up land, which has little or even negative ecological service functions, continuously consumes high-valued land covers, such as grassland and forest land, resulting in the decline of ESV. According to the results in Section 3.2, the rapid expansion of built-up land was one of the most obvious characteristics of LC change in China. Based on the transition matrix from 1992 to 2020, the ESV loss caused by built-up land expansion was estimated to be CNY 188.55 billion, making 61.2% of the total loss. China’s GDP data from 1992 to 2020, obtained from the data of the National Bureau of Statistics (http://www.stats.gov.cn/, accessed on 16 November 2022), were compared with China’s total ESV, as shown in Figure 14a. It can be seen that GDP and ESV basically exhibited opposite trends, indicating that economic development led to negative changes in the ecosystem to a certain extent. Before 2000, the surge of population led to increased demand for grain consumption, and the large expansion of built-up land and cropland occupied a large amount of forest land and grassland, which resulted in the decrease of ESV. The decrease of ESV from 2000 to 2015 was mainly due to the rapid expansion of built-up land mostly converted from cropland, forest land, and grassland. It is also worth noting that ESV did not continue to decline with the increase of GDP, and has rebounded since 2015, mainly due to the implementation of a series of ecological protection policies in China (grain for green policy, ecological red line policy, etc.), which led to a significant increase in the area of forest land. In general, the expansion of built-up land under the current stage of economic development was the most important factor leading to the decline of ESV in China.

A binomial equation was fitted between the ESV and GDP of China from 1992 to 2020, as shown in Figure 14b. It can be seen that with the growth of GDP, the ESV presented a sharp decline at first and then recovered somewhat, showing a positive U-shape, which may be similar to the Kuznets curve hypothesis proposed by environment economists [73]. In the process of transition from a low economic level to a high economic level, the ecological environment first declines and then rises with economic development after a certain point. This phenomenon was also analyzed in the study of Pan et al. [31]. GDP is bound to rise continuously in the future, but built-up land cannot keep expanding due to the limited supply of land resources. International experience showed that, with the development of the economy, the proportion of built-up land increased gradually, while the cropland and bare land decreased gradually. When the proportion of the tertiary industry rose rapidly and took the main position, the structural ratio of cropland and built-up land also tended to get stable, and the growth rate of built-up land was slower at this time [74]. Therefore, in the process of economic development, the strong dependence on built-up land expansion turns weak, that is, the intensity of demand for built-up land gradually weakens [75,76], which has been shown in Japan, South Korea, and other developed countries. As the largest developing country in the world, China is still at the stage where economic growth is mainly driven by enhancing the productive factors, and the trend of built-up land expansion changes almost simultaneously with GDP. Referring to international experience, after China’s economic development reaches a certain stage, the dependence of economic growth on the expansion of built-up land will gradually get weak, and the built-up land will gradually approach the limit and stop growing [77]. Then, the ESV will no longer decline or even start to rise.

Meanwhile, the imbalance of regional development makes the ESV of different regions in China vary with the change trend of GDP. Figure 14b also shows the relationship between ESV and GDP in China’s two provinces with the highest total GDP (Guangdong and Jiangsu) and two provinces with the lowest total GDP (Ningxia and Tibet). It can be seen that the relationship between ESV and GDP in developed areas (Guangdong and Jiangsu) also presents a positive U-shape, but is different in underdeveloped areas (Ningxia and Tibet). The reason is that when the economic level of underdeveloped areas is very low, ESV increases with the increase of GDP, which may due to the low urbanization level in this period, the slow expansion of built-up land (Figure 3), the small impact of human activities on the natural ecosystem, and climate change also leads to the warming and wetting of western China and the increase of vegetation coverage [78]. With the improvement of the development level, the underdeveloped areas will also show a similar change trend to the developed areas. Therefore, the relationship between ESV and GDP in China’s underdeveloped provinces may show an N-shape.

We compared reports regarding ESV in China in the literature and also the estimates of global ESVs by Costanza et al. [3] and Costanza et al. [51], presented in Table 4. We find that the ratio of ESV to GDP in China varies greatly among different studies, but mostly less than 1, except Pan et al. [31].

We further investigated the spatial difference and temporal variation of the ratio of ESV to GDP in China. As shown in Figure 15a, provinces in the eastern parts of China, which are economically more developed than the western provinces, had a lower ratio of ESV to GDP than less developed western provinces. From the variation of ESV/GDP in different areas of China illustrated in Figure 15b, the variation of ESV/GDP in different areas had a similar temporal pattern, that is, ESV/GDP declined continuously from high values at the beginning, but the decline rate became smaller and smaller, and finally approached a low value and got stable. The rate of decline in ESV/GDP was higher in underdeveloped areas than that in developed areas, which indicated that natural resources were more consumed for economic development in underdeveloped areas.

5. Conclusions

In this study, we assessed the spatiotemporal evolution of LC based on annual CCI-LC land cover data from 1992 to 2020, estimated the ecosystem service value (ESV) in China, and analyzed the spatiotemporal dynamics and its relationship with LC changes and GDP changes.

The overall accuracy of CCI-LC land cover products in the China region was high at 71.1% with a Kappa coefficient of 0.65. However, the user accuracy was low for cropland (59.8%) and grassland (38.6%), and the producer accuracy was low for forest land (67%), water (68.6%), and bare land (51.6%).

In China, grassland had the largest proportion among all LC types, with a multi-year average of 37.2% of the total area, followed by forest land (23.9%) and bare land (17.9%). The area of cropland in China increased generally first before 2004, then decreased after 2008; the area of forest land decreased before 2003, then increased after 2015; the area of grassland and bare land consistently decreased; and the area of built-up land continued to increase by 113,000 km², with an average annual growth rate of 5.3%. The primary characteristics of land cover transitions in China were the mutual conversion of cropland, forestland, and grassland as well as the ongoing increase of built-up land.

Forest land was the most important contributor of ESV in China, which makes 62.9% of the total ESV, followed by grassland (18.5%) and water (10.3%). The total ESV showed a sharp decline from CNY 19.07 trillion in 1994 to CNY 18.72 trillion in 2013, followed by a rebound and finally became stabilized at around CNY 18.77 trillion, with a total decrease of CNY 307.93 billion. The ESV was typically high in the southeast and low in the northwestern China, with the ESV in the three coastal regions (southeast, south-central, and northeast) being higher than that in other regions due to their water and heat conditions; most of the northern and northwestern areas in arid and semi-arid climates were constrained by natural conditions and had a high percentage of bare land, so there were large areas of low ESV value; forest land and water occupied a high proportion in southwest China, and there were also a small number of high ESV value areas.

Cropland, forest land, grassland, and water were the positive contributors to ESV change in China, while built-up land and bare land were the negative contributors to ESV change in China. The cold spots of ESV decrease were concentrated in the more economically developed areas in the east and south, which were mainly due to the transition of LC types from high ESV types (cropland, forest land) to low ESV types (built-up land) in these areas; the hotspots of rising ESV were mainly concentrated in the west, north, and some parts of southwest China. The conversion of LC types from low ESV types (bare land) to high ESV types (cropland, forest land, and grassland) was the main reason for the rising ESV in northwest China, while the increase of water and grassland was the main reason for the rising ESV in the central Tibetan Plateau.

The expansion of built-up land under the current stage of economic development was the most important factor leading to the decline of ESV in China. With the improvement of development level, the negative effect of economic development on ESV will gradually weaken. The relationship between ESV and GDP in eastern developed areas may present a positive U-shape, while that between ESV and GDP in western underdeveloped areas may be an N-shaped. The ratio of ESV to GDP and the rate of decline in ESV/GDP were higher in the western underdeveloped areas than in the eastern developed areas, which indicated that natural resources were more consumed for economic development in underdeveloped areas.