The Spatiotemporal Correlation between Human Activity Intensity and the Evolution of Ecosystem Service Value in the Songnen Plain, China

Abstract

For the main grain-producing areas worldwide that balance multi-tasks of grain production, ecological protection, and economic development, quantitatively revealing the correlation between human activity intensity (HAI) and ecosystem service value (ESV) is conducive to formulating adapted ecological protection policies and promoting the coordinated development of the regional economy, society, and ecosystem. In this paper, we took the Songnen Plain of China as an example, employed a modified Equivalent Factor Method (integrating socio-economic data, the normalized difference vegetation index (NDVI), and land use data), and the HAI Assessment Model (based on the data of land use, night-time light, and population spatial distribution) to measure the ESV and HAI in the Songnen Plain of China for the years 1990, 1995, 2000, 2005, 2010, 2015, and 2020. We further applied the standard deviational ellipse method, the coupled coordination degree model, and the bivariate spatial autocorrelation models to reveal the spatiotemporal dynamics and correlation characteristics of ESV and HAI. The results showed the following: (1) Temporally, the ESV declined from 950.96 billion yuan in 1990 to 836.31 billion yuan in 2015, and then increased to 864.60 billion yuan in 2020, with the total loss attributed to the significant decline in the ESV of the natural ecosystem. Spatially, the ESV in the western and northeastern regions was relatively high, with a significant increase in the northeast. (2) HAI showed an upward trend from 1990 to 2020, while the high HAI levels gradually shrank after reaching the peak in 2000. Low HAI levels were mainly distributed in the northeast and southwest, aligning with the ecological space, while high HAI levels were distributed in the middle and southeast. (3) The interaction between ESV and HAI was marked by a negative correlation, transitioning from conflict to coordination. The spatial pattern of HAI and ESV showed H (HAI)-L (ESV) and L-H clustering, with H-H and L-L scattered distributions. This study contributes to providing a framework, methods, and suggestions for the sustainable planning and utilization of land and ecological protection in order to offer scientific references for the Songnen Plain, other major grain-producing areas, and related studies.

1. Introduction

Ecosystem services, including provision services, regulation services, support services, and cultural services, refer to the life-supporting goods and services that humans obtain from the properties, processes, and functions of natural ecosystems directly or indirectly [1,2,3,4]; meanwhile, they are also the natural environmental conditions and benefits that human beings depend on for survival and development [5,6]. However, with urbanization, industrialization, and agricultural modernization as well as the rapid development of economic and population aggregation, human activities have altered the direction, structure, and rate of land use, which has damaged the supply and health of natural ecosystem services, and in turn, greatly weakened the value of ecosystem services (ESV) [7,8,9,10]. The global assessment report from the Intergovernmental Science-Policy Platform on Biodiversity and Ecosystem Services (IPBES) noted that three-quarters of the terrestrial environment and about 66% of the marine environment have been affected by human activities, which has and will continue to pose a significant threat to ecosystem services [11,12,13]. Meanwhile, the deterioration of ecosystems could, in turn, restrict the sustainable development of human societies [13]. Therefore, the 15th goal of the United Nations 2030 Agenda for Sustainable Development has been proposed to control the impact of human activities on terrestrial ecosystems and promote the harmonious coexistence between humanity and nature [14,15]. Accordingly, the study of the correlation between human activities and ecosystem services has become an important issue in the field of human–land interactions [16]. Based on the background of human–land conflict, quantitatively evaluating the changing trend of ESV and HAI, and portraying the correlation between human activity and ecosystem services are of great significance for revealing the degradation mechanisms and improvement methods of regional ecosystems and formulating policies that promote harmony between ecosystems and humans.

This study on the relationship between human activity intensity (HAI) and ESV depends on the accurate evaluation of ESV, the quantification of human activities, and the revelation of the relationship between them.

The calculation of ESV is to quantify the products and services provided by the ecosystem, which is crucial to improving the public awareness of ecosystem preservation and promoting regional sustainable development [17]. At present, many scholars have evaluated ESV by using the physical assessment method [18], energy assessment method [19], and monetary valuation method [1,17], and each method has its own advantages and disadvantages [6]. The monetary valuation method is widely used because its research results can be adopted by the socio-economic management system. Meanwhile, it can be used to conduct large-scale and long time series comparative studies, and is more friendly to regions that lack relevant basic data [17,20]. The monetary valuation method was first proposed by Costanza et al. [1], which is based on land use data and the classification of ecosystem service functions and evaluates the value of ecosystem services at the global scale by establishing the table of equivalence values per unit area of ecosystem services. Subsequently, in light of the findings of Costanza et al. [1], Chinese scholars Xie et al. [21] proposed the table of equivalent values per unit area of ecosystem services in China depending on the characteristics of China’s terrestrial ecosystems and the magnitude of biomass and is continuously improved to reflect the health status of Chinese ecosystems [17,20,22,23]. This kind of method uses land use data to calculate ESV, but the reality remains that the ESV per unit area varies even with the same land use types due to different ecosystem health and vegetation growth conditions. Therefore, the NDVI (normalized difference vegetation index) and the NPP (net primary productivity) have been applied by scholars to further modify the monetary valuation method in line with the table of equivalence values per unit area of ecosystem services [14,24], aiming at reflecting the change in ESV caused by the ecosystem’s structure and function. In conclusion, the existing and improved ecosystem service value assessment methods can provide useful references for this study.

Human activity, which is commonly defined by human activity intensity, is a general term for a range of behaviors or activities that humans adopt to meet their needs for survival and development [25,26]. Among the existing measures of HAI, some studies have adopted a single indicator to analyze human activities, such as land use data [27], night-time lighting data [28], population density [29], agriculture [30], and urbanization level [31]. Since the methods relying on one indicator can only reflect HAI from a certain aspect, which is a simplification of human activities, it is difficult to present HAI in different regions and stages in the process of socio-economic development, and the accuracy of this assessment is still to be discussed [8,32]. Conversely, other scholars have developed an integrated indicator system to evaluate HAI, including a variety of indicators. For example, Correa Ayram et al. [33] assessed the intensity of human impact by selecting seven variables, such as land use type, rural population density, distance to roads, distance to settlements, the fragmentation index of natural vegetation, the biomass index relative to natural potential, and the time of intervention on ecosystems in years, and the results demonstrated that HAI increased during the study period. Zhang et al. [13] selected 27 indicators involving five aspects, such as agricultural activities, industrial activities, tertiary industry activities, social development, and ecological engineering, to construct human activity impact indicators in order to evaluate the interaction between human activities and the ecological environment. However, the method of constructing comprehensive indicators encounters several challenges: variations in indicator selection across different locations and individuals, difficulties in obtaining extensive and long-term data series, discrepancies in weighting methods, and obstacles to regional comparisons [25,32]. Recently, considering the limitations of the two methods as well as the complexity and diversity of human activities, some scholars have begun to calculate HAI from the perspective of the power sources and impact of human activities by using three types of data: population density, night-time light, and land use [34,35]. This method can provide reliable and comparable data support for analyzing ecological problems and formulating policies in different regions, owing to the availability of data (night-time light data and land use data are derived from remote sensing data, population density data as regional basic data have data accessibility), limited indicators, comparability, and revealing characteristics of HAI.

In the study of the interaction between human activities and ESV, previous studies have employed the spatial autocorrelation analysis [36], the grey correlation analysis [13], and the coupling coordination degree analysis [37] to reveal the spatial–temporal correlation and coupling coordination relationship of the two systems at different scales (including the global scale [38], national scale [10], and regional scale [39]) based on Geoda [40] and ArcGIS software [33]. The results suggest that the relationship between human activities and ecosystem services is diverse, non-linear, or complex in different regions and stages of development [16,41,42]. For example, Huang et al. [37] found that there was an inverted U-shaped curve between human footprint and ESV in the Pearl River Delta Urban Agglomeration (PRDUA) from 2000 to 2020, indicating a non-linear coupling relationship between the two indicators. Ai et al. [40] elucidated the relationship between human disturbance intensity (HDI) and ESV using the Bivariate Local Moran Index and spatial autocorrelation analysis, and the results showed that areas with high ecosystem service values were relevant to low levels of HDI, while regions with low ecosystem service values were impacted by high levels of HDI. Quintas-Soriano et al. [43] discussed the relationship between land use change and ecosystem services in Spain and pointed out that different human activities have diverse impacts on ecosystem services, such as the fact that greenhouse horticulture negatively affects all ecosystem services, while intensive agricultural services have a positive effect on it. The variations in results are associated with the scale of the study and the regional natural, social, and economic characteristics. For China, a country undergoing rapid socio-economic development and transformation, and at the same time affected by global climate change, scholars have been actively exploring the relationship between human activities and ecosystem services in ecologically fragile areas such as arid and semi-arid regions and the Qinghai-Tibet Plateau [16,24], and the economically developed regions such as the Pearl River Delta Urban Agglomeration and the Chengdu-Chongqing Urban Agglomeration [37,40], and attempting to reveal the spatial and temporal correlation between them. However, there is a lack of theoretical and empirical studies on the major grain-producing areas, particularly in Northeast China, which plays an important role in guaranteeing food security and ecological security in the north. This study may provide a useful reference for other similar areas to promote socio-economic development, improve human well-being, and realize the virtuous cycle of ecosystems. In addition, the results of previous studies have focused more on the quantitative analysis of results [44], and the application of the results concentrated on the extraction of the main controlling factors [13] and the measures to promote the improvement of ESV [40]. Currently, there are few thorough studies on ecological policies that aim to promote socio-economic and ecological development in the region from the perspective of serving government management, especially in the context of climate change and the spatial variations of ecosystems.

The Songnen Plain, located in the world’s four major black soil belts, is also an important grain-producing area and commodity grain base in China [45], yet the ecological environment is relatively fragile [46]. With the rapid development of the economy and global warming in the past 30 years, the region is experiencing a series of problems such as the expansion of build-up land, the sharp decline of ecological land, serious soil erosion, the degradation of black soil, salinization, and the deterioration of the ecological environment [45,47]. Among them, the irrational land use (land use transformation) triggered by the high intensity of human disturbances has been proven to be the main reason for the decline in ecological service functions and quality in the study area [48]. As the cornerstone for ensuring national food security and an important ecological function area in northern China, it is necessary to balance the relationship between economic development and ecological protection and then promote the improvement of the ecological environment. However, existing studies on the ecosystem services of Songnen Plain have focused on wetland ecosystems [49,50], the impact of land use change on ESV [48], and the ecological environment of a certain part of the study area [36]. There are research gaps in this area. Firstly, a long time series study on the spatial and temporal patterns of ESV evolution for the whole region is insufficient. Secondly, there are few studies on the mechanisms of ESV change from the perspective of the interaction between human activities and ecosystems. Thirdly, there is a lack of policy recommendations to serve the high-quality development of the main grain-producing areas, considering grain production, ecological protection, and economic revitalization in the context of climate change. The objectives of this research are to (1) calculate the spatiotemporal changes of ESV and HAI in the Songnen Plain; (2) quantitatively reveal the spatial and temporal relationship between ESV and HAI over the past 30 years in the Songnen Plain; and (3) put forward the policy recommendations on ecological protection to realize the sustainable development of the region’s society, economy, and ecology. The research framework and methods employed in this study are suitable for other regions worldwide that are investigating ESV, HAI, and their interactive feedback relationships. The policy implications provided in this study for balancing food security, protecting the ecological environment, and promoting economic development, as well as addressing climate warming, are applicable to other major grain-producing areas facing similar challenges and threats.

2. Materials and Methods

2.1. Study Area

The Songnen Plain, one of the three major plains in Northeast China, located in the central part of Northeast China (42°30′–51°20′ N,121°40′–128°30′ E), has a total area of 230,300 km². It includes the provinces of Heilongjiang and Jilin and looks like a rhombus shape formed by the Songhua and Nenjiang rivers’ alluvial deposits (Figure 1), with an elevation ranging from 92 m to 1696 m. The terrain is characterized by a high value in the surroundings and a low value in the middle. It has a temperate semi-humid and semi-arid continental monsoon climate with four distinct seasons, and rain and heat simultaneously. The mean annual temperature ranges from −2 °C to 7 °C, and the mean annual precipitation occurs from 350 mm to 600 mm [45]. The main soil types are black soil, black calcium soil, and meadow soil [45]. The land use types include cropland, forestland, wetland, grassland, water area, built-up land, and unused land, among which wetlands are widespread, and most wetlands are inland saline wetlands concentrated in ecologically fragile areas transitioning from semi-humid to semi-arid areas. This area, with adequate conditions of light, heat, and water resources, is one of China’s five major grain-producing areas and an important commercial grain production base, which plays an important role in guaranteeing national food security [46]. Moreover, it is also an important energy and heavy industry base in China. However, during the last 30 years, with the rapid development of the economy and global warming, the region has experienced deforestation, the return of farmland to forests, the conversion of dry land to paddy land, overgrazing, urbanization, and the outflow of population, etc. The ecological environment in this region is highly fragile due to the continuing disturbances from human activities [46].

2.2. Data Sources and Processing

The data in this study mainly included land use data, socio-economic statistical yearbook data, night-time light data, and population spatial distribution data, etc.

2.2.1. Land Use Data

The land use remote sensing monitoring data for 1990, 1995, 2000, 2005, 2010, 2015, and 2020, with a resolution of 30 m × 30 m, were obtained from the Data Center for Resources and Environmental Sciences, Chinese Academy of Sciences (RESDC) [51], and projection transformation, mask extraction, and reclassification were performed in the ArcGIS platform. The land use types included cropland (dry land and paddy land), forestland (broad-leaved and shrub), grassland, water area, wetland, built-up land, and unused land, and the overall classification accuracy was above 90% [52].

2.2.2. Socio-Economic Data

The socio-economic data such as sown area and the production data of grain crops from 1990 to 2020 were acquired from the Statistical Yearbook of Heilongjiang Province and Statistical Yearbook of Jilin Province, the China Economic and Social Big Data Research Platform [53]. The data were used to calculate the value of ecosystem services per unit area of cropland.

2.2.3. The Normalized Difference Vegetation Index (NDVI)

The NDVI data from 1990 to 2020 (five-year interval) came from RESDC [51], which is a year-by-year maximum NDVI dataset calculated based on Landsat 5/8 remote sensing imagery from the United States of America, reflecting the distribution and change of vegetation cover. The data, at a 1 km spatial resolution, was used to adjust ESV in the study area.

2.2.4. Night-Time Light Data

The night-time light data for 1992, 1995, 2000, 2005, 2010, 2015, and 2020, reflecting the intensity of human social and economic activities on the land surface, were provided by the Harvard Dataset Platform [24,54,55]. Night-time light data in 1990 was replaced by 1992 data because of data unavailability. The spatial resolution of the night-time light data is 1 km.

2.2.5. Population Density Data

The population density data from 1990 to 2019 (five-year interval), revealing the spatial distribution pattern of population in a certain area, were derived from the China Population Spatial Distribution Kilometer Grid Dataset of RESDC [51]. The data have a spatial resolution of 1 km. Due to the lack of data in 2020, the data from 2019 were used as a substitute.

To ensure the accuracy of the data, we preprocessed the spatial data in ArcGIS. The spatial coordinate systems were uniformed to the Krasovsky_1940_Albers and the resolution of all data was unified to 30 m. The detailed information of the data is shown in

Table 1. Details about the data in this study.

Data	Description	Resolution	Time Period
Land use data	Land use data types included cropland (dry land and paddy field), forestland (broad-leaved and shrub), grassland, water area, wetland, build-up land, and unused land	30 m	Every five years from 1990 to 2020
NDVI	Reflecting changes in vegetation coverage and biomass	1 km	Every five years from 1990 to 2020
Night-time light data	Revealing the intensity of human social and economic activities on the land surface	1 km	Every five years from 1990 to 2020
Population density data	Spatial distribution of population	1 km	Every five years from 1990 to 2020
Crops	Calculation of the value of a standardized equivalence factor based on sown area and yield data	City scale	1990–2020

.

2.3. Methods

2.3.1. Evaluation of Ecosystem Service Value

This study evaluated ESV based on the table of equivalent values per unit area of ecosystem services proposed by Xie et al. [20] and coupled with the conditions of the Songnen Plain. We took the cell as the study unit to measure ESV and HAI. After many repeated experiments, the study area generated a total of 26,339 cells, and each cell had a size of 3 km × 3 km by using ArcGIS. In the calculation steps, the key of the method was to determine the value of a standard equivalent factor, the equivalent factor coefficient, and the correction coefficient of ESV, on the basis of which ESV was estimated.

(1) Value of the standard equivalent factor

The standard equivalent factor in the study area was determined by 1/7 of the economic value of the grain yield per unit area [21]; the expression is as follows:

$$\mathrm{D}=1/7 \times \mathrm{P} \times \mathrm{Q}$$

(1)

where $\mathrm{D}$ stands for the standard equivalent factor, P denotes the price per unit of grain yield (yuan/kg), and Q represents the grain yield per unit area (kg/hm²).

In this study, Q was the grain yield per unit area of the Songnen Plain from 1990 to 2020, which is 5479.02 kg/hm². Meanwhile, to eliminate the influence of grain price changes on ESV each year and to reflect the real ecosystem service function and ecological environment quality in the study area, P was determined by the price per unit of grain yield of the region in 2020, which is 2.73 yuan/kg. Thus, the standard equivalent factor for the area was calculated to be CNY 2136.82/hm².

(2) Equivalent factor coefficients and the correction coefficient of ESV

The equivalent factor coefficient of Xie et al. [20] was determined according to the vegetation biomass of China’s terrestrial ecosystem. For the Songnen Plain, the equivalent factor coefficient would be adjusted in accordance with the biomass of the region. The formula of the correction coefficient is as follows [16]:

$$\mathrm{F}={\displaystyle \frac{\overline{{\mathrm{NDVI}}_{\mathrm{s}}}}{\overline{\mathrm{ND}{\mathrm{VI}}_{\mathrm{c}}}}}$$

(2)

where $\overline{{\mathrm{NDVI}}_{\mathrm{s}}}$ indicates the average value of vegetation cover in the study area, $\overline{{\mathrm{NDVI}}_{\mathrm{c}}}$ designates the average value of vegetation cover in China, and $\mathrm{F}$ refers to the correction coefficient of ESV in the study area. The average correction coefficient of the plain was finally calculated to be 1.33.

In addition, since broad-leaved forests accounted for more than 85% of the arbor (coniferous, mixed coniferous, and broad-leaved) [56,57], the equivalent factor coefficient of forests other than shrubs adopted the value of broad-leaved forests. The grassland in the region is dominated by meadow [58], so the equivalent factor coefficient of grassland was equal to the meadow. Referencing the research by Chu et al. [59], we set the built-up land equivalent factor coefficient as the value of desert due to the lack of the equivalent factor coefficient for built-up land in the table of equivalent values per unit area of the ecosystem services of Xie et al. [20].

Finally, the ESV coefficients of different ecosystem service types in the study area were obtained by multiplying the standard equivalent factor of the region with the coefficients of the equivalent factor (

Table 2. The ESV coefficients of different ecosystem service types in the Songnen Plain (CNY/hm²).

Primary Types	Secondary Types	Pf	Dl	Bl	Su	Gl	Wa	Wl	Bl	Ul
Provision services	FP	3777.90	2361.19	805.58	527.79	611.13	2222.29	1416.71	27.78	0.00
	MP	250.01	1111.15	1833.39	1194.48	916.70	638.91	1388.93	83.34	0.00
	WS	−7305.79	55.56	944.47	611.13	500.02	23,028.51	7194.67	55.56	0.00
Regulation services	GR	3083.43	1861.17	6027.97	3916.79	3166.77	2138.96	5277.95	305.57	55.56
	CR	1583.38	1000.03	18,056.13	11,750.37	8389.16	6361.31	10,000.32	277.79	0.00
	EP	472.24	277.79	5361.28	3555.67	2777.87	15,417.16	10,000.32	861.14	277.79
	HR	7555.80	750.02	13,167.08	9305.85	6139.08	284,009.02	67,307.69	583.35	83.34
Support services	SR	27.78	2861.20	7361.34	4777.93	3861.23	2583.42	6416.87	361.12	55.56
	NC	527.79	333.34	555.57	361.12	305.57	194.45	500.02	27.78	0.00
	BP	583.35	361.12	6694.66	4361.25	3527.89	7083.56	21,861.81	333.34	55.56
Cultural services	AL	250.01	166.67	2944.54	1916.73	1555.60	5250.17	13,139.31	138.89	27.78

Abbreviations: FP, MP, WS, GR, CR, EP, HR, SR, NC, BP, and AL refer to food production, material production, water supply, gas regulation, climate regulation, environmental purification, hydrological regulation, soil retention, nutrient cycling, biodiversity protection, and esthetic landscape, respectively. Pf, Dl, Bl, Su, Gl, Wa, Wl, Bl, and Ul refer to paddy field, dry land, broad-leaved, shrub, grassland, water area, wetland, built-up land, and unused land, respectively.

).

(3) Calculation of ESV

ESV was measured using the following equation:

$$\mathrm{ESV}=\underset{\mathrm{k}=1}{\sum ^{\mathrm{n}}}\left({\mathrm{VC}}_{\mathrm{k}} \times {\mathrm{A}}_{\mathrm{k}}\right)$$

(3)

$${\mathrm{ESV}}_{\mathrm{f}}=\underset{\mathrm{k}=1}{\sum ^{\mathrm{n}}}\left({\mathrm{VC}}_{\mathrm{kf}} \times {\mathrm{A}}_{\mathrm{k}}\right)$$

(4)

In the formula, $\mathrm{ESV}$ indicates the estimated ecosystem service value; ${\mathrm{VC}}_{\mathrm{k}}$ refers to the ESV coefficient of land use type $\mathrm{k}$ (CNY/hm²); ${\mathrm{A}}_{\mathrm{k}}$ stands for the area of land use type $\mathrm{k}$ (hm²); ${\mathrm{ESV}}_{\mathrm{f}}$ is the value of the ecosystem service function $\mathrm{f}$; and ${\mathrm{VC}}_{\mathrm{kf}}$ denotes the value coefficient of ecosystem service function $\mathrm{f}$ for the land use type $\mathrm{k}$.

2.3.2. Standard Deviational Ellipse Method

The standard deviation ellipse was used to analyze the spatial distribution trend of ESV in different periods, in which the mean center, long, and short axes of the ellipse can reveal the dynamics of its center of gravity, extension direction, and range [60]. The formula is expressed as follows:

$$\stackrel{\mathrm{-}}{\mathrm{X}}={\displaystyle \frac{{\sum }_{\mathrm{i}=1}^{\mathrm{n}}{\mathrm{X}}_{\mathrm{i}}{\mathrm{w}}_{\mathrm{i}}}{{\sum }_{\mathrm{i}=1}^{\mathrm{n}}{\mathrm{w}}_{\mathrm{i}}}}, \stackrel{\mathrm{-}}{\mathrm{Y}}={\displaystyle \frac{{\sum }_{\mathrm{i}=1}^{\mathrm{n}}{\mathrm{y}}_{\mathrm{i}}{\mathrm{w}}_{\mathrm{i}}}{{\sum }_{\mathrm{i}=1}^{\mathrm{n}}{\mathrm{w}}_{\mathrm{i}}}}$$

(5)

$$\mathsf{\theta }={\displaystyle \frac{\arctan /({\sum }_{\mathrm{i}=1}^{\mathrm{n}}{\mathrm{x}}_{\mathrm{i}}^{‘2}-{\sum }_{\mathrm{i}=1}^{\mathrm{n}}{\mathrm{y}}_{\mathrm{i}}^{‘2})+\sqrt{({\sum }_{\mathrm{i}=1}^{\mathrm{n}}{\mathrm{x}}_{\mathrm{i}}^{‘2}-{\sum }_{\mathrm{i}=1}^{\mathrm{n}}{\mathrm{y}}_{\mathrm{i}}^{‘2})+4{({\sum }_{\mathrm{i}=1}^{\mathrm{n}}{\mathrm{x}}_{\mathrm{i}}^{’}{\mathrm{y}}_{\mathrm{i}}^{’})}^2}}{2\mathrm{x}{\sum }_{\mathrm{i}=1}^{\mathrm{n}}{\mathrm{x}}_{\mathrm{i}}^{’}{\mathrm{y}}_{\mathrm{i}}^{’}}}$$

(6)

$${\mathsf{\delta }}_{\mathrm{x}}=\sqrt{{\sum }_{\mathrm{i}=1}^{\mathrm{n}}{({\mathrm{x}}_{\mathrm{i}}^{’}\cos \mathsf{\theta }-{\mathrm{y}}_{\mathrm{i}}^{’}\mathrm{s}\mathrm{in}\mathsf{\theta })}^2/\mathrm{n}} {\mathsf{\delta }}_{\mathrm{y}}=\sqrt{{\sum }_{\mathrm{i}=1}^{\mathrm{n}}{({\mathrm{x}}_{\mathrm{i}}^{’}\sin \mathsf{\theta }-{\mathrm{y}}_{\mathrm{i}}^{’}\cos \mathsf{\theta })}^2/\mathrm{n}}$$

(7)

where $(\stackrel{\mathrm{-}}{\mathrm{X}}, \stackrel{\mathrm{-}}{\mathrm{Y}})$ signifies the weighted mean center; $({\mathrm{x}}_{\mathrm{i}},{\mathrm{y}}_{\mathrm{i}}$) represents the geometric center coordinates; ${\mathrm{w}}_{\mathrm{i}}$ represents the ESV of the i-th grid; $\mathsf{\theta }$ denotes the azimuth angle of the standard deviation ellipse; $\mathrm{x}’$ and $\mathrm{y}’$ refer to the deviation of the coordinates from the center of the study area to the mean center, respectively; and ${\mathsf{\delta }}_{\mathrm{x}}$ and ${\mathsf{\delta }}_{\mathrm{y}}$ stand for the $\mathrm{x}$ and $\mathrm{y}$ axes of the standard deviation ellipse’s spatial distribution for ESV.

2.3.3. Assessment of Human Activity Intensity

Considering the power sources of human disturbance to the ecosystem and its consequences [34], the HAI comprehensive evaluation model was constructed using the night-time light index, population density, and conversion coefficients of different land use types for the built-up land equivalent (CI) proposed by Xu et al. [25] to depict the intensity of human activity in the study area. The model is more applicable for calculating the ESV of the Songnen Plain, which has undergone urbanization, industrialization, LUCC, and population loss.

The CI in the study area was determined according to the land use types of the Songnen Plain (shown in

Table 3. Conversion coefficients of different land use types for the built-up land equivalent of the Songnen Plain [25].

Land Use Types		Explanations of Characteristic Signs	$\mathbf{CI}$
Cropland	Paddy field Dry land	Natural cover of land surface changes and annual crops are planted.	0.2
Forestland	/	Natural cover of land surface does not change and is not used.	0
Grassland	/	Natural cover of land surface does not change but is used.	0.067
Water area	River/lake	Natural cover of land surface does not change and is not used.	0
	Reservoir	Natural cover of land surface changes. The exchanges of air and heat are blocked.	0.6
Wetland	/	Natural cover of land surface does not change and is not used.	0
Built-up land	/	There are artificial insulation layers on the surface. The exchanges of water, nutrients, air, and heat are blocked.	1
Unused land	/	Natural cover of land surface does not change and is not used.	0

) [25]. The HAI could be calculated with Equation (8):

$$\mathrm{HAI}=\mathrm{nNTL}+\mathrm{pPD}+\mathrm{lCI}$$

(8)

where NTL represents the normalized night-time lighting index; PD denotes population density; CI refers to the conversion coefficients of different land use types for built-up land equivalents; and $\mathrm{n}$, $\mathrm{p}$, and $\mathrm{l}$ stand for the indicator weights of NTL, PD, and CI. Referring to the research results of Chen et al. [34], $\mathrm{n}$, $\mathrm{p}$, and $\mathrm{l}$ are set as 0.3, 0.3, and 0.4, respectively, and the HAI is standardized to [0, 100].

2.3.4. Coupled Coordination Degree Model

The degree of coupled coordination is widely used to evaluate the degree of coordinated dependence between different systems in an interaction relationship [61]. In this paper, the coupled coordination degree model was applied to quantify the degree of interaction between ESV and HAI in the study area. The specific formula is as follows:

$$\mathrm{C}=\sqrt{{\displaystyle \frac{\mathrm{ESV}\cdot \mathrm{HAI}}{{\left[\frac{\mathrm{ESV}+\mathrm{HAI}}{2}\right]}^2}}}$$

(9)

$$\mathrm{T}=\mathsf{\alpha }\mathrm{ESV}+\mathsf{\beta }\mathrm{HAI}$$

(10)

$$\mathrm{D}=\sqrt{\mathrm{C}\cdot \mathrm{T}}$$

(11)

where C refers to the degree of coupling between ESV and HAI, 0 ≤ C ≤ 1; a higher value indicates a closer relationship between the two systems. ESV and HAI are standardized values, T signifies the integrated development coefficient, α and β stand for the degree of importance of ESV and HAI; both α and β are assigned a value of 0.5 due to their equal importance [62], D designates the degree of coordination between ESV and HAI, 0 ≤ D ≤ 1; a larger value indicates that the coordination of the two systems is higher, and vice versa.

2.3.5. Bivariate Spatial Autocorrelation Models

Spatial autocorrelation is an important indicator reflecting the correlation between the attribute values or phenomena of adjacent spaces, and the characteristics of local spatial agglomeration, which is specifically divided into global spatial autocorrelation and local spatial autocorrelation (denoted by global Moran’s I and local Moran’s I, respectively) [62]. Given the great spatial correlation between HAI and ESV, this paper applied the bivariate spatial autocorrelation model to analyze the spatial correlation between the two indicators. The global Moran’s I index was calculated as follows:

$$\mathrm{I}={\displaystyle \frac{\sum _{\mathrm{i}=1}^{\mathrm{n}}\sum _{\mathrm{j}=1}^{\mathrm{n}}{\mathrm{W}}_{\mathrm{ij}}({\mathrm{Y}}_{\mathrm{i}}-\stackrel{\mathrm{-}}{\mathrm{Y}})({\mathrm{Y}}_{\mathrm{j}}-\stackrel{\mathrm{-}}{\mathrm{Y}})}{{\mathrm{S}}^2\sum _{\mathrm{i}=1}^{\mathrm{n}}\sum _{\mathrm{j}=1}^{\mathrm{n}}{\mathrm{W}}_{\mathrm{ij}}}}$$

(12)

$${\mathrm{S}}^2={\displaystyle \frac{1}{\mathrm{n}}}\sum _{\mathrm{i}=1}^{\mathrm{n}}({\mathrm{Y}}_{\mathrm{i}}-\stackrel{\mathrm{-}}{\mathrm{Y}}{)}^2$$

(13)

where I refers to the bivariate global spatial autocorrelation index; n is the number of spatial units; ${\mathrm{Y}}_{\mathrm{i}}$ and ${\mathrm{Y}}_{\mathrm{j}}$ denote the attribute values of spatial units i and j; $\stackrel{\mathrm{-}}{\mathrm{Y}}$ represents the average value of the observed values of all grids; and ${\mathrm{W}}_{\mathrm{ij}}$ stands for the spatial weight matrix based on the proximity relationship of spatial units i and j.

The local Moran’s I index is expressed as follows:

$${\mathrm{I}}_{\mathrm{lm}}^{\mathrm{i}}={\mathrm{z}}_{\mathrm{l}}^{\mathrm{i}}{\sum }_{\mathrm{j}=1}^{\mathrm{n}}{\mathrm{W}}_{\mathrm{ij}}{\mathrm{z}}_{\mathrm{m}}^{\mathrm{j}}={\displaystyle \frac{{\mathrm{X}}_{\mathrm{l}}^{\mathrm{i}}-{\stackrel{\mathrm{-}}{\mathrm{X}}}_{\mathrm{l}}}{{\mathsf{\sigma }}_{\mathrm{l}}}}{\sum }_{\mathrm{j}=1}^{\mathrm{n}}{\mathrm{W}}_{\mathrm{ij}}{\displaystyle \frac{{\mathrm{X}}_{\mathrm{m}}^{\mathrm{j}}-{\stackrel{\mathrm{-}}{\mathrm{X}}}_{\mathrm{m}}}{{\mathsf{\sigma }}_{\mathrm{m}}}}$$

(14)

where ${\mathrm{X}}_{\mathrm{l}}^{\mathrm{i}}$ designates the observed value of the attribute $\mathrm{l}$ of spatial unit i; ${\mathrm{X}}_{\mathrm{m}}^{\mathrm{j}}$ denotes the observed value of the attribute $\mathrm{m}$ of spatial unit $\mathrm{j}$; ${\mathsf{\sigma }}_{\mathrm{l}}$ and ${\mathsf{\sigma }}_{\mathrm{m}}$ stand for the variances of the observed values of attributes $\mathrm{l}$ and $\mathrm{m}$, respectively; ${\mathrm{I}}_{\mathrm{lm}}^{\mathrm{i}}$ represents the bivariate local spatial autocorrelation index; and ${\mathrm{z}}_{\mathrm{l}}^{\mathrm{i}}$ and ${\mathrm{z}}_{\mathrm{m}}^{\mathrm{j}}$ are the normalized values of the standard deviation of spatial cell $\mathrm{i}$ and cell $\mathrm{j}$ for attribute value $\mathrm{l}$ and attribute value $\mathrm{m}$.

This study’s technical framework is depicted in Figure 2.

3. Results

3.1. Land Use Changes in 1990–2020

The dominant land types in the Songnen Plain from 1990 to 2020 were dry land (about 49.2–52.4% of the total area), forestland (about 14.0–14.6% of the total area), grassland (about 7.5–10.7% of the total area), and wetland (about 7.0–8.9% of the total area), among which dry land made up the largest proportion, accounting for approximately 50% of the total area (

Table 4. Land use type changes in the Songnen Plain in 1990–2020.

Year	1990	1995	2000	2005	2010	2015	2020	1990–2020
Land Use Types	Area (km2)	Area (km2)	Area (km2)	Area (km2)	Area (km2)	Area (km2)	Area (km2)	Change (km2)	Change Rate (%)
Pf	9932.76	15,428.62	14,394.68	13,944.07	13,380.99	14,249.07	15,537.7	5604.94	56.43
DL	113,205.15	115,627.89	117,921.7	119,246.69	120,671	120,111.22	117,663.46	4458.30	3.94
Fl	33,699.9	32,750.95	32,795.18	32,737.8	32,331.93	32,291.23	33,093.11	−606.79	−1.80
Gl	24,692.52	23,727.24	18,589	18,565.31	19,106.9	18,940.74	17,373.05	−7319.47	−29.64
Wa	6600.36	6090.28	5662.97	5347.84	5430.60	5424.73	5471.97	−1128.39	−17.10
Wl	20,560.04	18,386.98	19,090.21	18,486.14	16,389.79	16,224.22	18,392.06	−2167.97	−10.54
Bl	9511.35	9614.69	9707.92	9807.06	10,482.53	10,688.42	10,988.51	1477.16	15.53
Ul	12,120.48	8696.16	12,161.02	12,188.61	12,530.07	12,393.71	11,781.35	−339.13	−2.80

Abbreviations: Pf, DL, Fl, Gl, Wa, Wl, Bl, and Ul refer to paddy field, dry land, forestland, grassland, water area, wetland, built-up land, and unused land, respectively.

). Regarding to the degree of change, the paddy fields, dry land, and built-up land showed an overall upward trend during 1990–2020, with growth rates of 56.43%, 3.94%, and 15.53% respectively, while grassland area had the highest loss of 29.64% with a decrease in area of 7319.47 km², followed by water area and wetland, which decreased by 17.10% and 10.54%, respectively.

In terms of the direction and intensity of land conversion during the period of 1990 to 2020 (

Table 5. Conversion matrix of land use change in the Songnen Plain during 1990–2020 (km²).

Year						2020
	Land Use Types	Pf	DL	Fl	Gl	Wa	Wl	Bl	Ul	Decreased
	Paddy field	-	1694.65	28.35	53.00	22.64	349.63	176.28	3.82	2328.37
	Dry land	4375.20	-	2066.61	861.61	241.58	524.81	2328.59	314.37	10,712.77
	Forestland	93.25	3140.08	-	402.53	96.49	391.30	78.49	30.72	4232.86
1990	Grassland	1190.63	6415.21	951.48	-	101.83	659.33	168.61	1296.17	10,783.25
	Water area	119.34	263.73	58.46	238.80	-	902.77	16.11	655.60	2254.80
	Wetland	1861.87	1286.72	214.09	867.03	514.73	-	67.25	401.14	5212.83
	Built-up land	152.88	1199.48	60.47	38.29	13.64	14.48	-	31.50	1510.75
	Unused land	140.15	1172.22	250.79	1003.44	148.74	205.00	152.50	-	3072.82
	Increased	7933.31	15,172.08	3630.25	3464.69	1139.66	3047.33	2987.81	2733.33	40,108.46

Abbreviations: Pf, DL, Fl, Gl, Wa, Wl, Bl, and Ul refer to paddy field, dry land, forestland, grassland, water area, wetland, built-up land, and unused land, respectively.

), the conversion between paddy fields and dry land was one of the significant land use change characteristics. Specifically, the area of paddy fields that was converted into dry land is 1694.65 km², while 4375.20 km² of dry land was converted into paddy fields, meaning that the net inflow of paddy fields was massive. Moreover, most of the increased cropland (paddy field and dry land) was converted from ecological land (forestland, grassland, and wetland). For instance, 3140.08 km² of forestland, 6415.21 km² of grassland, and 1286.72 km² of wetland were transferred into dry land; meanwhile, 1190.63 km² of grassland and 1861.87 km² of wetland were converted into paddy fields. Thirdly, the expansion of built-up land came primarily from the transfer of dry land, with an area of 2328.59 km². Moreover, the conversion of dry land to forestland, water area, and wetland (Ecological De-farming) could not be ignored (Table 2). Overall, the region represented the characteristics of cropland converted into build-up land, ecological land converted into cropland, and dry land transformed into paddy fields.

3.2. Temporal and Spatial Evolution of ESV

3.2.1. Temporal and Spatial Patterns of Change in ESV

The ESV of the Songnen Plain from 1990 to 2020 was calculated (

Table 6. Ecosystem service value and its change in the Songnen Plain from 1990 to 2020.

Land Use Types	ESV/(1010 CNY·hm−2)							ESV Change 1990–2020 (%)						ESV (%)
	1990	1995	2000	2005	2010	2015	2020	1990–1995	1995–2000	2000–2005	2005–2010	2010–2015	2015–2020	2020
Pf	1.07	1.67	1.56	1.51	1.45	1.54	1.68	55.33	−6.70	−3.13	−4.04	6.49	9.04	1.94
DL	12.61	12.88	13.14	13.28	13.44	13.38	13.11	2.14	1.98	1.12	1.19	−0.46	−2.04	15.16
Fl	20.47	19.33	19.80	19.76	19.95	19.93	20.09	−5.60	2.45	−0.21	0.98	−0.12	0.79	23.23
Gl	7.84	7.53	5.90	5.89	6.07	6.01	5.52	−3.91	−21.66	−0.13	2.92	−0.87	−8.28	6.38
Wa	23.03	21.25	19.76	18.66	18.95	18.93	19.09	−7.73	−7.02	−5.56	1.55	−0.11	0.87	22.08
Wl	29.71	26.57	27.59	26.71	23.68	23.44	26.58	−10.57	3.82	−3.16	−11.34	−1.01	13.36	30.74
Bl	0.29	0.29	0.30	0.30	0.32	0.33	0.34	1.09	0.97	1.02	6.89	1.96	2.81	0.39
Ul	0.07	0.05	0.07	0.07	0.07	0.07	0.07	−28.25	39.84	0.23	2.80	−1.09	−4.94	0.08
Total	95.10	89.57	88.10	86.18	83.93	83.63	86.46	−5.81	−1.64	−2.18	−2.62	−0.36	3.38	100

Abbreviations: Pf, DL, Fl, Gl, Wa, Wl, Bl, and Ul refer to paddy field, dry land, forestland, grassland, water area, wetland, built-up land, and unused land, respectively.

) using Formulas (3) and (4). The results indicated that ESV fell continuously and then rose in 1990–2020. The total value decreased from 950.96 billion yuan in 1990 to 836.31 billion yuan in 2015 and then increased to 864.60 billion yuan in 2020. In terms of the contribution of various land use types to ESV, taking 2020 as an example (Table 6, Figure 3), wetlands had the highest contribution of 30.74% (although only accounting for 7.99% of the total area), followed by forestland (23.23%) and water area (22.08%). It is noticeable that wetlands, forestland, and water areas were the primary contributors to ESV in the region. During the period of 1990–2020, the natural ecosystem of the study area (i.e., forestland, grassland, water area, and wetland) declined from 810.55 billion yuan to 712.73 billion yuan, which dropped by 12.07% (Table 6). Among them, the reduction in ESV for water area, wetland, and grassland made up 40.25%, 32.03%, and 23.76% of the total ESV loss, respectively. The artificial (semi-artificial, semi-ecological) ecosystem (i.e., paddy field, dry land, and built-up land) increased by 8.21%, from 139.74 billion yuan to 151.22 billion yuan, with paddy field and dry land contributing 52.78% and 43.28%, respectively.

The values of ecosystem service functions in the Songnen Plain are displayed in Figure 4. Taking the value in 2020 as an example, among the value of primary services, regulation services contributed the highest (66.73%), followed by support services (18.76%), provision services (9.71%), and cultural services (4.80%). Of the secondary service values, hydrological regulation service had the highest value, with a proportion of 40.82%, making it the primary ecosystem service function in the region. Moreover, ecosystem service functions such as biodiversity, soil protection and retention, gas regulation, environmental purification, and food production also contributed markedly, accounting for 8.89%, 8.93%, 7.22%, 6.23%, and 4.75%, while the contribution of nutrient cycling service was lowest at only 0.93%. According to the rate of change from 1990 to 2020 (Figure 4), the overall ecosystem service function showed a downward trend except for food production (which increased by 5.44%). The ESV of water supply decreased significantly by 31.23%, followed by hydrological regulation (−11.81%), esthetics landscape (−9.82%), and biodiversity (−9.42%). The declining rate of material production was the lowest (−2.09%).

To exhibit the spatial and temporal distribution characteristics of ESV, we used the natural breaks method and classified ESV into five levels, namely low value (ESV ≤ 2373), medium-low value (2373 < ESV ≤ 5023), medium value (5023 < ESV ≤ 9708), medium-high value (9708 < ESV ≤ 19,078), and high value (ESV > 19,078) (Figure 5). Overall, the ESV in the western and northeastern regions was relatively high, and was the opposite in other regions. The high-value areas of ESV mainly occurred in the ecological areas, including water areas, wetlands, and forested areas, while the low value of ESV overlapped with areas of intensive human activity and ecological vulnerability, especially located in the dry land, built-up land, and unused land. Moreover, there were no significant changes in most areas during 1990–2020 (Figure 5), and the areas with a significant increment in ESV were mostly clustered in water areas and wetlands, while the areas with a significant decline in ESV were centered on the northeast of Baicheng, the west of Daqing, and certain areas in Qiqihar.

3.2.2. The Spatial Pattern Evolution Characteristics of ESV

The spatial change characteristics of ESV in the Songnen Plain from 1990 to 2020 were analyzed based on Formulas (5)–(7), and the results are shown in

Table 7. Standard deviation ellipse parameters of ESV.

Year	Center of Gravity Longitude	Center of Gravity Latitude	Long Semi-Axis (km)	Short Semi-Axis (km)	Rotation (°)	Area (104 km2)
1990	125.42	46.43	157.80	241.09	14.87	11.95
1995	125.43	46.41	156.28	246.25	14.54	12.08
2000	125.49	46.45	158.70	242.64	14.56	12.09
2005	125.51	46.46	158.63	243.76	14.41	12.14
2010	125.50	46.45	161.15	244.10	14.45	12.35
2015	125.50	46.45	161.28	244.33	14.46	12.37
2020	125.52	46.47	160.06	244.39	13.51	12.28

, Figure 6. The area of the standard deviation ellipse showed an overall expanding trend, meaning that the difference in ESV between the inside and outside of the ellipse was reduced, and areas of poor ecological quality were improved.

The centroid of ESV was located in the central part of the study area, which is the southwest of Suihua City. Temporally, the centroid of ESV first migrated to the southeast by 1.98 km and then moved to the northeast by 12.56 km as a whole, which confirmed that the ESV in the northeast has increased significantly and the ecological structure and function have improved. Meanwhile, the rotation angle of the standard deviation ellipse fluctuated around 14°, indicating that the spatial distribution of ESV followed a northeast–southwest direction, which is basically consistent with the distribution pattern of the study area.

3.3. Temporal and Spatial Distribution Patterns of HAI

In this paper, we measured HAI in the Songnen Plain from 1990 to 2020 using Formula (8), and analyzed the spatiotemporal pattern characteristics at the grid scale (Figure 7 and Figure 8). According to the data distribution, HAI was classified into five levels: low impact (HAI ≤ 7), medium-low impact (7 < HAI ≤ 15), medium impact (15 < HAI ≤ 20), medium-high impact (20 < HAI ≤ 30), and high-impact (HAI > 30).

For the temporal dimension, the highest proportion of the low HAI level showed an overall slowly declining tendency (decreased by 10.14%), and the medium-low HAI level exhibited a slight decrease (reduced by 4.27%). However, the medium-high and medium HAI levels showed a fluctuating upward trend generally (increased by 25.78%), confirming that the study area’s HAI gradually increased. Moreover, the high HAI level gradually shrank after reaching the peak in 2000, occupying 4.16% (Figure 7).

For the spatial dimension, HAI was generally characterized by a low level in the northeast and southwest, and a high level in the middle and southeast during 1990–2020 (Figure 8). The high HAI level mostly occurred in Qiqihar, Suihua, Harbin, and Changchun and was dispersed as dots in the urban centers, including Harbin, Changchun, and Daqing, which are relatively economically developed regions. The medium-high HAI level overlapped with the medium HAI level in space and was distributed in the shape of a “7”. The areas with a low HAI level were mainly located in the ecological space, including forestland areas in the northeastern and eastern regions as well as wetland and water area in the western part of the region.

3.4. Analysis of Coupling Coordination Degree between HAI and ESV

The coupling degree and coordination degree of HAI and ESV in the Songnen Plain from 1990 to 2020 were calculated based on Formulas (9)–(11) (Figure 9). The coupling degree was relatively stable and within the range of 0.7–0.8, which was in the running-in stage [13], indicating that HAI and ESV interacted with each other closely. Meanwhile, the coordination degree increased from 0.2645 in 1990 to 0.2670 in 2020, and the value of “center of mass” was 0.2651, revealing that the coordination level was in a moderate disorder stage [62], but it clearly showed an upward trend. This study indicates that, although the high HAI level had an adverse effect on the ecosystem service, the two systems are transitioning from conflict to coordination.

3.5. Bivariate Spatial Autocorrelation Analysis of HAI and ESV

Based on the bivariate spatial autocorrelation model, we constructed the spatial weight matrix and plotted a Moran’s I scatterplot by using the GeoDa spatial analysis tool (Figure 10). The global Moran’s I values between HAI and ESV were significant with P = 0.001 at 99.9% confidence level and were −0.276, −0.266, −0.288, −0.264, −0.252, −0.243, and −0.242, respectively, for the seven years. The results indicate that there is a negative spatial correlation between the two in the study area. During 1990–2020, the Moran’s I values were all less than 0 and had the weakest correlation in 2020, indicating that human activity and ecosystems are shifting from conflict to coordination. The results can be corroborated by the study on the degree of coupling coordination.

By using the local spatial autocorrelation analysis, the bivariate LISA clustering map of HAI and ESV in the study area was obtained (Figure 11a). In Figure 11a, H-H stands for high HAI and high ESV; L-L denotes low HAI and low ESV; L-H signifies low HAI and high ESV; and H-L represents high HAI and low ESV. Taking the 2020 LISA clustering map as an example, the region had significant spatial heterogeneity. The areas with the H-L cluster had the widest distribution (20.75% of the total area), concentrating in the majority of Changchun, the western part of Harbin, the eastern part of Qiqihar, and the northern and southern regions of Suihua. The L-H areas mainly occurred in the water area and the eastern forestland. The areas with H-H cluster were scattered in the study area, accounting for only 2.73% of the region (Figure 11b). The L-L areas were mainly located in the parts of Qiqihar, Suihua, Baicheng, and Songyuan. In general, the study area was characterized by the clustering of H-L and L-H, and the scattered distribution of H-H and L-L.

4. Discussion

4.1. Spatial–Temporal Change Characteristics of Land Use and ESV

This study identified the significant spatial and temporal heterogeneity of ESV in the study area. Temporally, ESV in the Songnen Plain exhibited a trend of initially declining and then increasing from 1990 to 2020. Moreover, the rate of decline slowed down distinctly after 2010 and reached the lowest value in 2015 (Table 6). This may be the result of the intensity and direction of human activities on the ecosystem, which changed the land use pattern and regional ESV. From 1990 to 2010, the Songnen Plain, serving as an important national commodity grain base and heavy industry base, faced the dual pressures of ensuring food security and promoting economic development [48]; therefore, the areas of wetland, forestland, grassland, and water area declined considerably, while paddy field, dry land, and built-up land significantly increased (Table 4). The decline in ESV during this period is largely attributed to the conversion of natural ecosystems into artificial (semi-artificial) ones. After 2010, the rate of ESV decline slowed, and ESV bottomed out after 2015, which is related to the downturn in the regional economy, population outflow, the construction of ecological civilization, and the implementation of the second round of policies on returning farmland to forestland and grassland, etc. In detail, the population of the region decreased by 13%, from 4.19 $\times$ 10⁷ people in 2010 to 3.65 $\times$ 10⁷ people in 2020 [63,64], which greatly lessened the pressure on the ecological environment. Meanwhile, China put forward the goal of “ecological civilization construction” in 2007 [65]; subsequently, in 2012, ecological protection and sustainable socio-economic development were placed in equally important positions, deepening the concept of ecological protection [66]. Moreover, in 2014, China launched a new round of returning farmland to forestland (grassland and wetland) [67]. In the year 2015, policy documents such as the Guidelines on Accelerating Ecological Civilization and the Integrated Reform Plan for Promoting Ecological Civilization were successively issued [66]; concurrently, relevant measures were implemented in order to improve the eco-environment, such as the drawing of the ecological protection red line and compensating for ecological protection activities [66]. These policies and measures had a positive impact on stabilizing and increasing the area of natural ecosystems and improving ecosystem service capabilities, which enabled ESV increases. In terms of the degree of change in various service functions, the ESV of water supply experienced the most significant decline, followed by hydrological regulation, esthetic landscape, and biodiversity (Figure 4). In the Songnen Plain, due to climate warming, the decline of precipitation [68], and the expanding of croplands (paddy fields and dry lands) (Table 4), the high water resource demand for agricultural irrigation [69,70] has led to a significant decline in the ESV of water supply. In addition, the long-term transformation of natural ecosystems to artificial (semi-artificial) ones in the study area (Table 5) diminished the functions of biodiversity and hydrological regulation, which in turn weakened the value of the corresponding ecosystem services. Spatially, regions with high ESV values were primarily concentrated in ecological areas such as water areas, wetlands, and forestland, while areas with low ESV values were mainly located in areas with intensive human activities and ecological fragility (dry land, built-up land, and unused land). Overall, this spatial pattern is related to the fact that natural ecosystems have higher ecosystem service value equivalents and provide more ESV than artificial (semi-artificial) ones.

4.2. The Feedback Mechanism between HAI and ESV

The results of the coupling coordination degree between HAI and ESV illustrated that human activities and ecosystem services were closely related and interactive (Figure 9). This is consistent with previous studies [13,16,37]. Zhang et al. [13] analyzed the coupling relationship between human activities and ESV in Xinjiang during 2000–2015 based on the coupling degree model, and the results demonstrated that the two were highly coupled (all higher than 0.797), indicating a strong interaction. Huang et al. [37] observed that there was a non-linear coupling relationship between human footprint and ESV in the Pearl River Delta Urban Agglomeration from 2000 to 2020, with mutual influences and constraints. Eigenbrod et al. [71] concluded that there would be an interaction between urbanization and ecosystem service provision in the UK, meaning that ecosystem services were affected by human activities. The main reason for the mutual feedback between human activities and ecosystem service is that human production and survival are inseparable from the support of ecological products and ecosystem service functions. For example, provision services provide raw materials such as food and wood for human survival, and cultural services fulfil spiritual pleasures [1,72]. Meanwhile, human activities also have an effect on the ecosystem by altering the structure and function of the ecosystem [73,74], which includes negative pressure [74,75,76] and positive effect [77] (Figure 12). It can be seen that human activities and ecosystem services are interrelated and mutually restrained.

Furthermore, the temporal pattern of the coordination degree and the bivariate spatial autocorrelation between ESV and HAI in the Songnen Plain showed that high HAI levels exerted a strong negative impact on ESV, but they are transitioning from conflict to coordination. The results demonstrated that the positive effects of human activities on ecosystems in the study area were gradually increasing, and the negative pressures were lowering steadily (Figure 9 and Figure 10). Relevant studies have shown that human activities had varied bearing on ecosystems in different regions and stages of socio-economic development [13,77,78]. Zhang et al. [13] pointed out that, since the agricultural expansion from 2000 to 2015, the ESV in Xinjiang declined with the exception of the value of food supply, indicating an overall negative influence of human activities on the ESV in the region. Sinchembe and Ellery et al. [78] noticed that the hydrological health of wetlands in South Africa has suffered from road construction and invasive alien plants, resulting in wetland degradation and a reduction in the supply capacity of wetland indirect ecosystem services, suggesting a negative implication of human activities on wetland ecosystems. Conversely, Mao et al. [77] found that ecological conservation policies and ecological projects can increase ecological land and improve regional ecosystem service capacity, which suggested that human activities can aid the betterment of ecosystem services as well. The possible mechanisms for the impact of human activities on ecosystems are as follows: In the stage of rapid urban development, along with the advancement of urbanization, industrialization, and agricultural modernization, human beings intend to achieve urban development and economic growth through the expansion of built-up land and cropland, as well as the compression of ecological space [79]. At the same time, high-intensity human activities also affect the physical, chemical, biological, and hydrological properties of the soil [8,80]. Such changes in the quantity and quality of ecosystems bring about the alterations in ecosystem structure and function and the decline of ecosystem service capacity [37,81]. When ecosystems reach a certain level of degradation and begin to constrain socio-economic development, the government is inclined to enact ecological conservation policies; simultaneously, with the deepening of the concept of ecological conservation and the promotion of ecological conservation technology, humans tend to balance economic development and ecological conservation [77]. For example, measures such as crop rotation and fallow, rest grazing, and grazing prohibition are adopted to improve the ecological service capacity of cropland and grassland [82]; the red line of the urban development boundary is delineated to restrict the expansion of built-up land; Ecological De-farming is promoted to increase the supply area of ecosystem service [46]; and the improvement of ecosystem service capacity is achieved by reducing pollutant emissions and through ecological restoration [77] (Figure 12).

4.3. Policy Implications

Benefiting from ecological conservation policies and its implementation, pressures on ecological systems induced by human activities have been alleviated to some extent (the proportion of the high HAI level in the study area decreased (Figure 7)). Nevertheless, most of these policies were promoted at the national level and fell short of considering the special ecosystem structure and function and the socio-economic development situation of the Songnen Plain. In addition, the research area also undertakes multiple tasks such as revitalizing Northeast China, ensuring national food and ecological security, protecting the black soil, and facing the challenges of climate change. Taking into account the regional characteristics, spatial differences, and climate challenges, the following recommendations can be given.

Firstly, differentiated policies should be formulated to achieve accurate management in heterogeneity areas. In the H-L areas with relatively well developed economic activity, such as Changchun City, Harbin City, and Qiqihar City (Figure 11a), the occupation of ecological land and cropland by build-up land should be strictly controlled, and the existing build-up land should be fully used. Additionally, green industries should be actively developed to reduce pollution emissions. Whereas, in the L-H zones, namely the areas of water area, wetland, and forestland (Figure 11a), it is urgent to identify ecological protection zones, strictly abide by ecological red lines, construct ecological corridors, and enhance the regional ecosystem service capacity. For cropland that needs to harmonize grain production and ecological services, the red line for 1.8 billion acres of cropland should be safeguarded to prevent the uncontrolled expansion of cropland; meanwhile, conservation tillage and other technologies should be promoted to increase soil organic matter, protect the biodiversity of cropland ecosystems, reduce the dosage of chemical fertilizers and pesticides, and improve ecosystem service capacity and health [82].

Secondly, it is necessary to formulate targeted ecological policies for ecologically fragile areas to curb ecosystem deterioration. In the Songnen Plain, there are many rolling hilly regions (Keshan County, Baiquan County, and Kedong County in Qiqihar City) with severe soil and water erosion [36]. In these areas, service functions of cultivated ecosystems, such as soil conservation, water supply, and food production, have been declining under the dual influence of natural and human factors. Moreover, the region also includes ecologically fragile areas such as saline-alkali land [9] and agro-pastoral ecotone [83]. In light of this, the intensity of human disturbance in the fragile ecological environment should be weakened, and the ecological restoration should also be implemented. In addition, it is essential to monitor the carrying capacity and carrying threshold of the ecosystem dynamically, provide early warning of ecological disasters, and restore the degraded ecosystem.

Thirdly, the government should unveil ecological policies to actively respond to the challenges posed by climate warming on the regional ecosystems. Previous studies have shown that climate warming and economic interests have led to a gradual expansion of rice planting areas and a northward shift of the cultivation area [69,77]. This phenomenon, the significant expansion of the paddy field area at the end of the study period in the Songnen Plain, is also proved by this study (Table 4). The expansion of paddy fields, along with the existing dry land planting, will further aggravate the consumption of water resources, lower the capacity of water supply in ecosystem services, increase additional pressure on groundwater, river, and lake resources [13,84], and also bring about soil degradation and a series of related problems [85]. Consequently, policymakers ought to enact relevant policies in response to climate change, controlling the scale of dry land-to-paddy conversion and implementing water-saving irrigation technology and soil and water conservation measures so as to promote the coordinated development of humans and ecosystems.

Finally, it is worth mentioning that, apart from the Songnen Plain, such grain production areas worldwide as the Midwestern United States, Brazil’s Cerrado Region, and the Pampas Region of Argentina [86,87,88] may face multiple challenges in sustainable food production, ecosystem conservation, economic development, and adaptation to global climate change. These regions (and their countries), when formulating sustainable development policies, could conduct research based on the “assessment-correlation-policy” framework of this study to enhance the specificity and scientific rigor of their policies. They could also apply the “modified Equivalent Factor Method” and the “HAI Assessment Model” from this study, which adequately consider local social, economic, and natural conditions, to assess ESV and HAI, offering high applicability. Based on further understanding of the spatiotemporal correlation between regional ESV and HAI and using the research outcomes and policy measures from the Songnen Plain as a reference, targeted strategies could be developed for diverse regions and ecologically fragile areas. These strategies may include monitoring ecosystem dynamics, controlling land conversion, expanding ecological reserves, and implementing ecological technologies, aiming to achieve regional socio-economic synergy and sustainable development.

5. Conclusions

In this research, we analyzed the Songnen Plain land use changes using land use data every five years between 1990 and 2020 and measured ESV and HAI, employing the modified Equivalent Factor Method and HAI Assessment Model (based on land use, night-time light, and population spatial distribution data). Furthermore, we explored the spatiotemporal correlations between the two utilizing the coupled coordination degree and the bivariate spatial autocorrelation model. The conclusions are as follows:

(1)

The dominant land types in the Songnen Plain from 1990 to 2020 were dry land (about 49.2–52.4% of the total area), forestland (about 14.0–14.6% of the total area), grassland (about 7.5–10.7% of the total area), and wetland (about 7.0–8.9% of the total area); overall, the region represented the following characteristics: cropland converted into build-up land, ecological land converted into cropland, and dry land transformed into paddy fields.

(2)

From 1990 to 2020, ESV declined from 950.96 billion yuan in 1990 to 836.31 billion yuan in 2015, and increased to 864.60 billion yuan in 2020. The ESV of natural ecosystems showed a declining tendency, while artificial (semi-artificial) ones showed a rising trend. Except for the significant increase in ESV of food production (5.44%), the ESV of other ecosystem service functions showed a downward trend. Among them, the ESV of water supply decreased dramatically (31.23%). In terms of spatial distribution, the ESV in the western and northeastern regions was relatively high, and the ESV in other regions was the opposite. The high ESV areas were mainly distributed in the ecological areas such as water areas, wetlands, and forestland. The distribution of the low ESV areas coincided with the areas of intensive human activities and ecological fragility, such as dry land, build-up land, and unused land. Moreover, ESV increased significantly in the northeast during the study period.

(3)

During 1990–2020, HAI showed an upward trend, and the spatial distribution characteristic was generally low in the northeast and southwest, as well as high in the middle and southeast. The high HAI levels gradually shrank after reaching their peak in 2000 and were dispersed as dots in the urban centers, including Harbin, Changchun, and Daqing, which are relatively economically developed regions; the areas of medium-high HAI levels overlapped with the areas of medium HAI levels in a ‘7’ shape; the areas of low HAI levels were highly coincident with the ecological space.

(4)

HAI had a strong interaction with ESV (the coupling degree was between 0.7 and 0.8), and the degree of coordination was in a moderate disorder (from 0.2645 to 0.2670). Although human activities had a negative impact on ESV, they were transitioning from conflict to coordination. For the spatial dimension, the two systems were H-L and L-H clustering and H-H and L-L scattered distribution.

(5)

In the future, the decision makers should formulate differentiated policies. In economically developed areas (H-L areas) like Changchun, Harbin, and Qiqihar, strict control over land occupation and the promotion of green industries are crucial. For ecological zones (L-H areas), it is urgent to strictly abide by ecological red lines and construct ecological corridors. For fragile regions, it is essential to weaken the intensity of human disturbance and implement targeted restoration efforts. Additionally, they should control the scale of dry land-to-paddy conversion, promote water-saving irrigation technologies, and implement soil and water conservation measures to address climate warming.

(6)

The “assessment-correlation-policy” research framework proposed in this study is suitable for characterizing and understanding the mechanisms of ecosystem services and human activities. The approaches employed in this study, which fully consider regional uniqueness, data availability, and inter-regional comparability, can be applied to quantitatively studying the feedback relationship between ecosystem services and human activities in other similar regions. Additionally, the policies proposed in this study could be extended to other major grain-producing areas globally, offering Chinese theoretical insights and practical references for these regions.