Assessing Land Resource Carrying Capacity in China’s Main Grain-Producing Areas: Spatial–Temporal Evolution, Coupling Coordination, and Obstacle Factors

Abstract

The land resources in the main grain-producing areas (MGPAs) provide a solid foundation for grain production, and promotion of the sustainable utilization of land resources in these areas is crucial for ensuring national food security. To comprehensively assess the land resource carrying capacity (LRCC) in China’s MGPAs, we utilized the driver-pressure-state-impact-response (DPSIR) framework and applied the analytic hierarchy process (AHP) and entropy weight (EW) method to analyze the spatial–temporal evolution of LRCC in China’s MGPAs from 2000 to 2020. By establishing a coupling coordination model, we explored the coupling coordination relationships among LRCC subsystems and identified key factors hindering the balanced development of LRCC using an obstacle degree model. The research results indicate that the LRCC in China’s MGPAs generally shows an increasing trend with a continuously growing rate, with the LRCC in the southern regions surpassing that in the northern regions. The overall coupling coordination of LRCC demonstrates an increasing trend, although the growth rate is decreasing. The coupling coordination level of LRCC in the southern regions is generally higher than that in the northern regions, and the gap in coordination levels between various regions is narrowing. The state and pressure subsystems significantly influence the balanced development of LRCC. Indicators such as arable land area per capita, grain production yield per unit area of arable land, grain production yield per capita, forest land area per capita, and grassland area per capita play vital roles in the development of LRCC. Based on these findings, we have put forward targeted recommendations.

1. Introduction

Land resources are one of the most important finite resources on earth [1,2], supporting human production and the development of all kinds of constructive activities [3,4]. According to statistics, land surface alone accounts for less than one third of the Earth’s surface area [5]. As Malthus (1798) pointed out, land cannot infinitely produce the necessary consumer goods for human beings and, if the pace of population growth is faster than that of land production, the maximum population will be determined by the extreme productive capacity of land resources [6]. With the rapid development of industries and the explosive growth of population in the process of urbanization, the contradiction between limited land resources and the larger scale and higher intensity of social and economic activities is becoming increasingly severe [3,7]. Data from the Food and Agriculture Organization of the United Nations show that the land use area per capita has declined annually from 4.25 hectares/person in 1961 to 1.65 hectares/person in 2021 at an average annual rate of 1.61% [8]. Research has pointed out that once the pressure exerted by human activities exceeds the carrying capacity of land resources, it will result in irreversible damage to the land resource system [9]. Therefore, ensuring the sustainable use of land resources is an inevitable choice for achieving sustainable development [10,11], which has become a global concern [12,13].

Carrying capacity, a term used in the field of engineering mechanics, refers to the relationship between a load and the carrier object that carries the particular load [3,14]. Subsequently, it has been widely introduced and applied in various research fields, including biology, ecology, and many other research fields related to human activities and the resource environment, and its meaning has also been expanded [15]. Land resource carrying capacity (LRCC) has been defined by the United Nations Educational, Scientific and Cultural Organization as “the intensity of human activities that can be carried by a region while maintaining an acceptable standard of living”, which appears to be the most widely accepted definition [14,16]. As global land development intensifies, the human–land relationship has become increasingly strained, leading to a comprehensive study of LRCC worldwide [17]. Previous research has focused on agricultural land [18] and urban land [19,20], with a scope ranging from a single city [21] to larger spatial regions such as major cities [22,23] and urban agglomerations [24]. In terms of research content, spatial and temporal changes [25,26], improvement measures [25,27], and prediction of evolutionary trends [3] have been studied extensively. With regard to methodology, the evaluation of LRCC has gradually evolved from using a single grain index focusing solely on population [23,25] to a comprehensive index assessed by an evaluation indicator system that considers various human activities [3,13].

Although these studies have provided valuable information to enhance the understanding of LRCC, there are also some limitations to these existing studies. First, the lack of a robust theoretical framework to support the selection of evaluation indicators, leading to a relatively vague understanding of the interrelationships among the subsystems within land resources. Second, previous studies have examined the coupling and coordination between LRCC and external factors, such as industrial development [7] and poverty incidence [27], but insufficient attention has been paid to the internal coupling and coordination of LRCC itself. Third, little research has identified the factors that impede balanced development of LRCC, making it difficult to propose specific strategies for improving the carrying capacity of land resources. Lastly and most importantly, most of the previous studies on LRCC in China have been conducted in the coastal cities or economically developed areas, while ignoring the specificity and importance of the main grain-producing areas (MGPAs). Ensuring food security in China is highly dependent on stable grain production in the MGPAs, which account for about 78% of the country’s total grain production [28]. Land resources in the MGPAs are used not only for grain production but also for regional socioeconomic activities, exserting enormous pressure on the land resources in the MGPAs [29,30]. This fact is demonstrated by the significant decline in the area of arable land, with a decrease of 1.98 million hectares observed from 2012 to 2021 [31]. Although China has prioritized the protection of arable land resources as a fundamental national policy, this decreasing trend presents a formidable challenge to both China’s and the world’s food security.

Therefore, the present study took the MGPAs of China as the study areas, and based on the driver-pressure-state-impact-response (DPSIR) framework, the indicator system for evaluating LRCC was determined, and the analytic hierarchy process (AHP) and entropy weight (EW) method was used to more accurately and rationally determine the weights of the evaluation indicators of the DPSIR framework. Based on the above, we analyzed the characteristics and patterns of the LRCC of China’s MGPAs in both spatial and temporal dimensions during 2000–2020, established a degree of coupling coordination model to explore the coupling coordination relationships among the subsystems of LRCC under the DPSIR framework, and identified the main factors hindering the balanced development of LRCC through the obstacle degree model. The research results could provide a scientific basis for the Chinese government to formulate measures to ensure that the demand for land resources from human activities in the MGPAs is under reasonable control. This is of great significance for achieving the sustainable use of regional land resources in the MGPAs and ensuring the national food security of China in the future.

2. DPSIR Framework

The DPSIR framework was initially proposed by the Organization for Economic Co-operation and Development, subsequently developed by the European Environment Agency, and has been adopted by the United Nations [32]. It systematizes the cause–effect relationships among interacting components of social, economic, and environmental systems [33]. Due to its potential and usefulness in providing clear and meaningful explanations [34], the DPSIR framework has been widely used in policymaking and research related to the management and protection of water [35], land [36], and marine [37] resources.

Key performance indicators (KPIs) quantify and simplify data that is not easily observable in a clear manner, serving as the primary tool for measuring and quantifying the DPSIR framework [38,39]. The DPSIR framework consists of five subsystems with different implications [40]:

1.

Driver

The driver subsystem refers to human social and economic activities that may cause changes in regional land resources. Population growth and economic development are the most fundamental drivers [41], because a growing population and economy lead to increasing demand and consumption of land resources [33,42]. The natural population growth rate and urbanization rate are commonly used to measure population growth [43,44,45]. GDP growth rate, GDP per capita, and share of tertiary industry could be used to measure economic development [43,44,45,46]. In addition, according to Gao et al. (2022) and Fan et al. (2023), the value of agricultural output per capita, livestock and poultry population per capita, and value of industrial output per capita could reflect the unique economic characteristics of an MGPA [44,47].

2.

Pressure

The pressure subsystem represents a direct or indirect cause of the regional land resources changes under the human social and economic activities. According to Wang et al. (2018), Tan et al. (2021), and Zhang et al. (2022), population density, annual disposable income per capita, and GDP density could measure the pressure from human social and economic activities [41,43,48]. Agricultural activities in an MGPA impose pressure on land resources, the application of chemical fertilizers, pesticides, plastic film, and emissions of livestock and poultry manure per unit area of arable land could be used to assess agricultural pollution [43,44,46,49,50]. According to Luo et al. (2020) and Gao et al. (2022), pressure on land resources caused by industry in the MGPAs could be assessed through chemical oxygen demand emissions in industrial wastewater per unit area of land and industrial solid waste emissions per unit area of land, reflecting industrial pollution [20,47].

3.

State

The state subsystem denotes the situation or conditions of the regional land resources. According to China’s current land use classification standards, land resources in the MGPAs are generally represented by eight types, including arable land, garden land, forest land, grassland, wetland, construction land, transportation land, and water bodies [31,51]. Therefore, the per capita area of these eight types of land resources was used as an evaluation indicator to evaluate the quantity of land resources [44,45,48]. Considering the critical importance of arable land in the MGPAs, grain production yield per unit area of arable land was chosen to measure the quality of land resources [43,44,48].

4.

Impact

The impact subsystem refers to feedback results of regional land resources on the economy, society, and ecosystem. Excessive exploitation of land resources can exacerbate the risk of natural disasters by increasing social and economic vulnerability and reducing the resilience of livelihoods [52]. Geological hazards (i.e., landslides, rock falls, debris flows, ground fissures, and subsidence) are natural disasters that are directly related to land resources [53]. According to Zhang et al. (2022), Shan et al. (2022), and Fan et al. (2023), the effect of land resources on social well-being could be measured by the urban registered unemployment rate and Engel coefficient [43,44,54]. In addition, according to Zhou et al. (2013), Cheng et al. (2017), and Shi et al. (2021), grain production yield per capita, arable land area growth rate, rate of crops affected by natural disasters, and rate of crops damaged by natural disasters could reflect the function of grain production of an MGPA [55,56,57].

5.

Response

The response subsystem stands for positive and effective measures and countermeasures taken by humans to adapt to or improve regional land resources [43,55,58,59,60]. According to Tan et al. (2021), He and Wang (2022), and Fan et al. (2023), the proportion of afforested area, control rate of forest diseases, pests and rodents, rate of street and road cleaning in built-up areas, harmless treatment and disposal rate of domestic waste, completed industrial pollution control investments per unit area of land, and environmental protection expenditure per unit area of land could be evaluation indicators [44,48,61]. In addition, with regard to grain production in the MGPAs, the degree of agricultural mechanization per unit area of arable land and proportion of effective irrigated area were chosen as evaluation indicators to reflect the function of grain production of the MGPAs [57].

The essence of the DPSIR framework is a causal chain of these five subsystems. Specifically, drivers exert pressures on regional land resources, resulting in changes in states of regional land resources, and subsequently leading to a series of impacts on human and regional land resources, which promote responses as feedback to the drivers [40,59,62,63,64]. The interrelationship between the five subsystems of the DPSIR framework is illustrated in Figure 1.

The DPSIR framework can provide substantial answers to the following questions concerning the MGPAs of China [34]: What is happening to the land resources (state)? What are the causes of these changes (driver and pressure)? What are the consequences of these changes on human and natural resources (impact)? What measures have been taken to address these changes (response)? Therefore, the DPSIR framework was employed for the evaluation of the LRCC of the MGPAs in China in the present study.

3. Materials and Methods

3.1. Study Area

The division of China’s MGPAs originated from the reform of the grain circulation system in 2001 [65]. Based on the general characteristics of grain production and consumption, as well as historical traditions and resource disparities, 31 provinces in China (excluding Hong Kong, Macau, and Taiwan) are categorized into three major grain functional areas: the MGPAs (13 provinces), the balanced production and marketing areas (11 provinces), and the major grain-marketing areas (7 provinces). The study areas in the current study are 13 provinces in the MGPAs of China, consisting of Heilongjiang, Liaoning, Jilin, Inner Mongolia, Hebei, Jiangsu, Anhui, Jiangxi, Shandong, Henan, Hubei, Hunan, and Sichuan. The green areas in Figure 2 depict the geographical distribution of the MGPAs in China.

Geographically, the study areas can be further divided into five regions: Northeast China (Heilongjiang, Liaoning, and Jilin), North China (Inner Mongolia and Hebei), East China (Jiangsu, Anhui, Jiangxi, and Shandong), Central China (Henan, Hubei, Hunan), and Southwest China (Sichuan). These areas are predominantly characterized by plains or low hills with warm and humid climates and fertile soils that are conducive to grain production [25]. In 2022, the grain sown area of the MGPAs reached 89,029 thousand hectares, accounting for about 75.24% of the total grain sown area in China. In the same year, the MGPAs achieved grain production of 53,718 million tons, representing about 78.25% of China’s total grain production [28].

3.2. Data Source

Evaluating LRCC requires a large amount of data, and existing official statistical yearbooks and databases provide comprehensive and authoritative information related to LRCC. In addition, official data also have the advantages of comparability, which facilitate cross-regional and comparative temporal analysis. Many previous studies have utilized official data as the foundation for their research, so this study also chooses to use official statistical data as the primary source of data [22,66]. The study uses panel data from 2000 to 2020 for the 13 MGPAs in China. The social and economic data were collected from the National Bureau of Statistics of China (http://data.stats.gov.cn, accessed on 27 May 2023), while land use data (with a resolution better than 1 m) were primarily collected from the National Land Survey Achievements Sharing Application Service Platform (https://gtdc.mnr.gov.cn, accessed on 27 May 2023). Agricultural data were obtained from the China Rural Statistical Yearbook, environmental data were sourced from the China Environmental Statistical Yearbook, and urban construction-related data were obtained from the China Urban and Rural Construction Statistical Yearbook. Missing data were complemented using either the linear interpolation method or the linear extrapolation method. The data collected from the above steps need to undergo the necessary transformation to provide the final assessment indicators, and the specific calculation methods are listed in Table S1.

3.3. Methods

3.3.1. Construction of the Evaluation Indicator System

The evaluation indicators for each subsystem within the DPSIR framework were identified by referring to the KPIs used in previous relevant studies and taking into account the characteristics of land resources in the MGPAs of China. A total of 42 variables were selected to construct the LRCC evaluation index system for the MGPAs in China, as depicted in

Table 1. Details of evaluation indicators for LRCC based on the DPSIR framework.

Subsystems	Indicators	Units	Attribute
Driver	D1 Natural population growth rate	‰	Negative
	D2 Urbanization rate	%	Negative
	D3 GDP growth rate	%	Positive
	D4 GDP per capita	CNY/person	Positive
	D5 Share of tertiary industry	%	Positive
	D6 Value of agricultural output per capita	CNY/person	Positive
	D7 Livestock and poultry population per capita	heads/person	Positive
	D8 Value of industrial output per capita	CNY/person	Positive
Pressure	P1 Population density	people/km2	Negative
	P2 Annual disposable income per capita	CNY/person	Positive
	P3 GDP density	CNY/m2	Positive
	P4 Chemical fertilizer usage per unit area of arable land	t/km2	Negative
	P5 Pesticide usage per unit area of arable land	t/km2	Negative
	P6 Plastic film usage per unit area of arable land	t/km2	Negative
	P7 Emissions of livestock and poultry manure per unit area of arable land	kg/km2	Negative
	P8 COD emissions in industrial wastewater per unit area of land	t/km2	Negative
	P9 Industrial solid waste emissions per unit area of land	t/km2	Negative
State	S1 Arable land area per capita	m2/person	Positive
	S2 Garden land area per capita	m2/person	Positive
	S3 Forest land area per capita	m2/person	Positive
	S4 Grassland area per capita	m2/person	Positive
	S5 Wetland area per capita	m2/person	Positive
	S6 Construction land area per capita	m2/person	Positive
	S7 Transportation land area per capita	m2/person	Positive
	S8 Water body area per capita	m2/person	Positive
	S9 Grain production yield per unit area of arable land	kg/m2	Positive
Impact	I1 Frequency of occurrence of geological hazards	times	Negative
	I2 Direct economic losses caused by geological hazards	million CNY	Negative
	I3 Urban registered unemployment rate	%	Negative
	I4 Engel coefficient	%	Negative
	I5 Grain production yield per capita	kg/person	Positive
	I6 Arable land area growth rate	%	Positive
	I7 Rate of crops affected by natural disasters	%	Negative
	I8 Rate of crops damaged by natural disasters	%	Negative
Response	R1 Proportion of afforested area	%	Positive
	R2 Control rate of forest diseases, pests, and rodents	%	Positive
	R3 Rate of street and road cleaning in built-up areas	%	Positive
	R4 Harmless treatment disposal rate of domestic waste	%	Positive
	R5 Completed industrial pollution control investments per unit area of land	CNY/km2	Positive
	R6 Environmental protection expenditure per unit area of land	CNY/km2	Positive
	R7 Degree of agricultural mechanization per unit area of arable land	kW/km2	Positive
	R8 Proportion of effective irrigated area	%	Positive

.

3.3.2. Determination of Evaluation Indicator Weights

A combination of subjective and objective weighting methods has commonly been used in previous studies to determine the weights of the evaluation indicators of the DPSIR framework [34,54,61]. This can avoid the subjective arbitrariness of subjective weighting method and address the limitation of the objective weighting method that it tends to ignore the subjective information of the decision maker [54]. Therefore, the AHP-EW method was used to determine the weights of the evaluation indicators of the DPSIR framework in the current study. By combining subjective judgment and objective data, the weights of each evaluation indicator are more accurate and reasonable.

1.

AHP method

The AHP method was used to obtain subjective weights based on expert opinions, providing logical and credible results. The basic steps of the AHP method are as follows [67,68]:

Step 1: The judgment matrix is constructed by pairwise comparison among the evaluation indicators with each layer. The matrix is determined by an expert group consisting of three associate professors and two professors who have long been engaged in land resource management research. The comparisons are made using Saaty scaling [39], as shown in

Table 2. Saaty’s scale.

Rating	Definition	Interpretation
1	Equal importance	Two indicators have equal importance
3	Somewhat more important	One indicator has slightly more importance than the other
5	Much more important	One indicator has more importance than the other
7	Very much more important	One indicator has very more importance than the other
9	Absolute importance	One indicator has absolute importance over the other
2, 4, 6, 8	Intermediate values	Intermediate values indicate intermediate rating. They are used when compromise is needed

.

The judgment matrix $X$ is expressed as follows:

$$X=x_{mm}=\left[\begin{array}{cc}\begin{array}{cc}{x}_{11}& x_{12}\\ x_{21}& x_{22}\end{array}& \begin{array}{cc}\cdots & x_{1m}\\ \cdots & x_{2m}\end{array}\\ \begin{array}{cc}⋮& ⋮\\ x_{m1}& x_{m2}\end{array}& \begin{array}{cc}\ddots & ⋮\\ \cdots & x_{mm}\end{array}\end{array}\right]$$

(1)

where the comparison between evaluation indicators $i$ and $j$ is $x_{ij}$, and the comparison between evaluation indicators $j$ and $i$ is $x_{ji}=\frac{1}{x_{ij}}$. When $i=j$, ${x_{ij}=x}_{ji}=1$. $i=\mathrm{1,2},\cdots ,m$; $j=\mathrm{1,2},\cdots ,m$.

Step 2: Normalized to obtain a vector $Z_j$, and the calculation formula is as follows:

$$Z_j=\sqrt[m]{x_{j1}{x}_{j2}\dots x_{jm}}$$

(2)

Step 3: Calculate the weight coefficient (i.e., the subjective weight) $S_j$ of each evaluation indicator according to the following formula:

$$S_j=\frac{Z_j}{{\sum }_{j=1}^m{Z}_j}$$

(3)

Step 4: Consistency index ($CI$) and consistency ratio ($CR$) are used to examine the consistency. The calculation formula of $CI$ is as follows:

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

(4)

where ${\lambda }_{max}$ is the maximum eigenvalue of the judgment matrix, and the calculation formula is as follows:

$${\lambda }_{max}={\sum }_{i=1}^m\frac{{\left(Xw\right)}_i}{n{w}_i}$$

(5)

The larger the $CI$ value, the worse the degree of consistency. When $CI=0$, the judgment matrix has completion consistency.

The calculation formula of $CR$ is as follows:

$$CR=\frac{CI}{RI}$$

(6)

where $RI$ is the average random consistency index of the matrix, and the $RI$ values of the 1–9 order judgment matrix are shown in

Table 3. Value of average random consistency index.

Matrix Order	1	2	3	4	5	6	7	8	9
$RI$	0	0	0.58	0.90	1.12	1.24	1.32	1.41	1.45

. The judgment matrix can be considered to have satisfactory consistency only when $CR<0.1$.

2.

EW method

The EW method was used to derive the objective weights according to the data itself, and its evaluation results are systematic and explanatory. The specific process of the entropy weight method is as follows [59,69]:

Step 1: Assuming that there are $n$ objects to be evaluated and $m$ evaluation indicators, the original matrix can be constructed:

$$X=\left[\begin{array}{cc}\begin{array}{cc}{x}_{11}& x_{12}\\ x_{21}& x_{22}\end{array}& \begin{array}{cc}\cdots & x_{1m}\\ \cdots & x_{2m}\end{array}\\ \begin{array}{cc}⋮& ⋮\\ x_{n1}& x_{n2}\end{array}& \begin{array}{cc}\ddots & ⋮\\ \cdots & x_{nm}\end{array}\end{array}\right]$$

(7)

where $x_{ij}$ is the data of the $i$-th objects to be evaluated under the $j$-th evaluation indicator in the original data, $i=\mathrm{1,2},\cdots ,n,j=\mathrm{1,2},\cdots ,m$.

Step 2: In order to avoid the impact of evaluation indicator differences, the original data set is normalized by using the formula:

Positive indicators:

$$Y_{ij}=\frac{x_{ij}-\mathit{min}\left(x_j\right)}{\mathit{max}\left(x_j\right)-\mathit{min}\left(x_j\right)}$$

(8)

Negative indicators:

$$Y_{ij}=\frac{{\mathit{max}\left(x_j\right)-x}_{ij}}{\mathit{max}\left(x_j\right)-\mathit{min}\left(x_j\right)}$$

(9)

where $j$ represents the $j$-th evaluation indicator, and $i$ represents the $i$-th province among all 13 provinces, $i=\mathrm{1,2},\cdots ,n$; $j=\mathrm{1,2},\cdots ,m$. $Y_{ij}$ is the normalized value of the $j$-th evaluation indicator of the $i$-th objects, and the value range is [0,1]. $x_{ij}$ is the original data of the $j$-th evaluation indicator of the $i$-th objects. $\min \left(x_j\right)$= the minimum value in the original data of the $j$-th evaluation indicator, and $\max \left(x_j\right)$= the maximum value in the original data of the $j$-th evaluation indicator.

Step 3: Calculate the entropy value $E_j$ of the $j$-th evaluation indicator using the following formula:

$$E_j=-\frac{1}{ln\left(n\right)}{\sum }_{i=1}^n{p}_{ij}\mathit{ln}\left(p_{ij}\right)$$

(10)

where $ln$ is the natural $log$, $p_{ij}$ refers to the proportion of the $i$-th objects to be evaluated in the $j$-th evaluation indicator, and the calculation formula is:

$$p_{ij}=-\frac{Y_{ij}}{{\sum }_{i=1}^n{Y}_{ij}}$$

(11)

Step 4: Calculate the entropy weight (i.e., the objective weight) $O_j$ of each evaluation indicator according to the formula as follows:

$$O_j=\frac{1-E_j}{{\sum }_{j=1}^m(1-E_j)}=\frac{1-E_j}{m-{\sum }_{j=1}^m{E}_j}$$

(12)

3.

Comprehensive evaluation model

The AHP method and the EW method are combined to calculate the comprehensive weight $C_j$, and the formula for the calculation is as follows [70,71]:

$$C_j=\frac{S_j{O}_j}{{\sum }_{j=1}^m{S}_j{O}_j}$$

(13)

After standardizing the data and calculating the comprehensive weights of the evaluation indicators, the score for the LRCC in a given year of a province in the MGPAs is obtained using the following comprehensive evaluation model:

$${LRCC}_i={\sum }_{j=1}^m{C}_j{Y}_{ij}$$

(14)

where ${LRCC}_i$ represents the LRCC score of the $i$-th province in a given year, $i=\mathrm{1,2},\cdots ,n$. $C_j$ represents the comprehensive weight of the $j$-th evaluation indicator, $j=\mathrm{1,2},\cdots ,m$. $Y_{ij}$ is the standardized value of the $j$-th evaluation indicator of the $i$-th province.

3.3.3. Determination of LRCC Score Classification

Both a higher LRCC score and a lower LRCC score do not represent a suitable carrying capacity of land resources. Instead, an ideal LRCC score that ensures a balance between the carrying capacity of land resources and the pressure exerted by human activities is deemed appropriate for realizing the sustainable use of land resources [25,72]. Specifically, in regions with relatively better socioeconomic development, LRCC scores tend to be higher due to the relatively larger consumption of land resources. However, it is better to descend the LRCC scores to a relatively lower level within these regions to avoid irreversible damage to land resources. Similarly, in regions with relatively lagging socio-economic development, LRCC scores are generally lower. It is better to elevate the LRCC scores to a relatively higher level within these regions to avoid ineffective use and wastage of land resources. The analysis above is depicted in Figure 3.

Since the LRCC system is complex, it is difficult to directly determine the classification of the LRCC scores. Therefore, with reference to Luo et al. (2022) [22], the LRCC scores were classified into five categories by calculating the average quintiles of the LRCC scores ($Q_1–Q_4$) of the 13 provinces in the MGPAs in a certain year, as shown in Figure 2. Provinces with LRCC scores lower than or equal to the first quintile ($Q_1$) are considered to be great surplus provinces, indicating that their socioeconomic development is substantially below the carrying capacity of land resources. Provinces with LRCC scores greater than or equal to the first quintile ($Q_1$) but lower than the second quintile ($Q_2$) are defined as surplus provinces, meaning that the socioeconomic development in these provinces is below the carrying capacity of land resources. Provinces with LRCC scores greater than or equal to the second quintile ($Q_2$) but lower than the third quintile ($Q_3$) are classified as balanced provinces, indicating the socioeconomic development approach to the carrying capacity of land resources within these provinces. Provinces with LRCC scores greater than or equal to the third quintile ($Q_3$) but lower than the fourth quintile ($Q_4$) are classified as overloaded provinces, which indicates that the socioeconomic development is above the carrying capacity of land resources within these provinces. Provinces with the socioeconomic development substantially above the carrying capacity of land resources are defined as seriously overloaded, with LRCC scores greater than or equal to the fourth quintile ($Q_4$).

3.3.4. Coupling Coordination Analysis

The degree of coupling coordination model was used to reflect the interaction between subsystems of LRCC in the MGPAs of China, evaluated based on the DPSIR framework. The coupling degree refers to the level of interaction between the five subsystems within the LRCC [73]. The specific calculation formula is as follows:

$$CP={\left[\frac{D\times P\times S\times I\times R}{{\left(\frac{D+P+S+I+R}{5}\right)}^5}\right]}^{\frac{1}{5}}$$

(15)

where $CP$ is the coupling degree of the five subsystems within the LRCC, which ranges from 0 to 1. The closer the $CP$ is to 1, the greater the coupling degree, which represents a stronger interaction between the five subsystems. The high coupling degree of a subsystem means that a small change in the subsystem can have significant effects on other subsystems, leading to a tightly coupled system. $D$, $P$, $S$, $I$, and $R$ are the respective comprehensive evaluation results of the driver, pressure, stress, impact, and response subsystems.

Although the coupling degree can indicate the degree of mutual influence, it does not reflect the extent of coordinated development degree between various subsystems [74]. Therefore, the following degree of coupling coordination model is constructed:

$$CC=\sqrt{CP\times CD}$$

(16)

$$CD=\alpha D+\beta P+\gamma S+\delta I+\epsilon R$$

(17)

where $CC$ is the coupling coordinated degree of the five subsystems within the LRCC, ranging from 0 to 1. The closer the $CC$ is to 1, the greater the coupling coordinated degree, which represents a higher level of synergy and cooperation between the five subsystems. $CD$ is the comprehensive coordination index, which is used to reflect the orderly or disorderly coupling extent of each subsystem. $\alpha$, $\beta$, $\gamma$, $\delta$, and $\epsilon$ are undetermined contribution coefficients of the comprehensive coordination index, and $\alpha +\beta +\gamma +\delta +\epsilon =1$. In the present study, each subsystem was considered equally important for the comprehensive coordination index, which is the current practice in most studies [75,76,77]. Therefore, $\alpha$, $\beta$, $\gamma$, $\delta$, and $\epsilon$ all take the value 0.2. Referring to previous studies [75,76,78], the division for the coupling coordinated degree is shown in

Table 4. Classification of the coupling coordinated degree.

Value Range	Coordination Level	Coupling Coordination Type
$0\le CC<0.1$	1	Extreme disordered recession
$0.1\le CC<0.2$	2	Severe disordered recession
$0.2\le CC<0.3$	3	Moderate disordered recession
$0.3\le CC<0.4$	4	Mild disordered recession
$0.4\le CC<0.5$	5	On the verge of disordered recession
$0.5\le CC<0.6$	6	Barely coordinated development
$0.6\le CC<0.7$	7	Primary coordinated development
$0.7\le CC<0.8$	8	Intermediate coordinated development
$0.8\le CC<0.9$	9	Well-coordinated development
$0.9\le CC<1$	10	High-quality coordinated development

.

3.3.5. Obstacle Factor Analysis

The obstacle degree model was introduced to reveal specific hindering factors for the balanced development of LRCC [79,80]. The specific calculation formula is as follows:

$$O_i=\frac{(1-S_i)\times w_i}{{\sum }_{i=1}^n(1-S_i)\times w_i}$$

(18)

where $O_i$ is the obstacle degree of the evaluation indicator $i$ in a total of $n$ evaluation indicators, $i=\mathrm{1,2},\cdots ,n$, indicating the degree of influence of the single evaluation indicator on the LRCC. $S_i$ is the standardized value of the evaluation indicator $i$, $w_i$ is the comprehensive weight of the indicator $i$.

$$U_j=\sum O_{ij}$$

(19)

where $U_j$ is the obstacle degree of the $j$ subsystem in a total of $m$ subsystems, $j=\mathrm{1,2},\cdots ,m$. The larger the value, the greater the obstacle effect of the subsystem on the improvement of the LRCC.

4. Results

4.1. Weights of Evaluation Indicators

The subjective, objective, and comprehensive weights of the evaluation indicators determined sequentially by the AHP, the EW, and the AHP-EW methods are shown in

Table 5. Weights of the evaluation indicators for LRCC.

Subsystems	Codes	Subjective Weights	Objective Weights	Comprehensive Weights
Driver	D1 Natural population growth rate	0.0158	0.0174	0.0113
	D2 Urbanization rate	0.0135	0.0118	0.0065
	D3 GDP growth rate	0.0162	0.0062	0.0041
	D4 GDP per capita	0.0086	0.0318	0.0112
	D5 Share of tertiary industry	0.0071	0.0129	0.0038
	D6 Value of agricultural output per capita	0.0257	0.0467	0.0492
	D7 Livestock and poultry population per capita	0.0257	0.0171	0.0180
	D8 Value of industrial output per capita	0.0066	0.0278	0.0075
Pressure	P1 Population density	0.0085	0.0096	0.0033
	P2 Annual disposable income per capita	0.0060	0.0310	0.0076
	P3 GDP density	0.0049	0.0618	0.0124
	P4 Chemical fertilizer usage per unit area of arable land	0.0137	0.0116	0.0065
	P5 Pesticide usage per unit area of arable land	0.0160	0.0080	0.0052
	P6 Plastic film usage per unit area of arable land	0.0160	0.0069	0.0045
	P7 Emissions of livestock and poultry manure per unit area of arable land	0.0217	0.0054	0.0048
	P8 COD emissions in industrial wastewater per unit area of land	0.0163	0.0063	0.0042
	P9 Industrial solid waste emissions per unit area of land	0.0163	0.0052	0.0035
State	S1 Arable land area per capita	0.0581	0.0591	0.1407
	S2 Garden land area per capita	0.0339	0.0267	0.0371
	S3 Forest land area per capita	0.0220	0.0581	0.0524
	S4 Grassland area per capita	0.0126	0.0822	0.0424
	S5 Wetland area per capita	0.0156	0.0611	0.0391
	S6 Construction land area per capita	0.0364	0.0320	0.0477
	S7 Transportation land area per capita	0.0251	0.0362	0.0372
	S8 Water body area per capita	0.0300	0.0310	0.0381
	S9 Grain production yield per unit area of arable land	0.1212	0.0296	0.1468
Impact	I1 Frequency of occurrence of geological hazards	0.0414	0.0004	0.0007
	I2 Direct economic losses caused by geological hazards	0.0573	0.0016	0.0038
	I3 Urban registered unemployment rate	0.0120	0.0028	0.0014
	I4 Engel coefficient	0.0083	0.0092	0.0031
	I5 Grain production yield per capita	0.0625	0.0408	0.1045
	I6 Arable land area growth rate	0.0625	0.0020	0.0051
	I7 Rate of crops affected by natural disasters	0.0216	0.0045	0.0040
	I8 Rate of crops damaged by natural disasters	0.0216	0.0073	0.0065
Response	R1 Proportion of afforested area	0.0111	0.0207	0.0094
	R2 Control rate of forest diseases, pests, and rodents	0.0069	0.0064	0.0018
	R3 Rate of street and road cleaning in built-up areas	0.0075	0.0242	0.0074
	R4 Harmless treatment disposal rate of domestic waste	0.0069	0.0115	0.0033
	R5 Completed industrial pollution control investments per unit area of land	0.0121	0.0522	0.0259
	R6 Environmental protection expenditure per unit area of land	0.0193	0.0480	0.0380
	R7 Degree of agricultural mechanization per unit area of arable land	0.0283	0.0236	0.0274
	R8 Proportion of effective irrigated area	0.0272	0.0113	0.0126

. The $CR$ value was 0.06, which was less than 0.1, indicating that the consistency of the matrix was considered to have passed the test. Regarding the weights of the subsystems, the subsystems were ranked in descending order of weight values as state (0.5815), impact (0.1291), response (0.1258), driver (0.1116), and pressure (0.0520). In terms of the weights of individual evaluation indicators, the weight value of grain production yield per unit area of arable land (S9, 0.1468), arable land area per capita (S1, 0.1407), and grain production yield per capita (I5, 0.1045) were higher than the other evaluation indicators. This indicates that the state subsystem is the major subsystem that determines changes in the LRCC in the MGPAs, grain production yield and the area of arable land are the main indicators affecting the LRCC in the MGPAs.

4.2. Evaluation Results of LRCC

4.2.1. Results of LRCC Scores

Table 6. LRCC scores of the MGPAs.

Regions	Provinces	2000	2005	2010	2015	2020	Growth (%)
Northeast China	Heilongjiang	0.1755	0.1909	0.2457	0.3058	0.3399	3.3603
	Liaoning	0.1989	0.2058	0.2434	0.2866	0.3208	2.4189
	Jilin	0.1999	0.2185	0.2648	0.2895	0.3430	2.7363
	Average	0.1914	0.2051	0.2513	0.2940	0.3346	2.8322
North China	Inner Mongolia	0.2358	0.2446	0.2802	0.3103	0.3627	2.1763
	Hebei	0.2317	0.2907	0.3242	0.3961	0.4268	3.1015
	Average	0.2338	0.2677	0.3022	0.3532	0.3948	2.6542
East China	Jiangsu	0.2050	0.2216	0.2568	0.2992	0.3360	2.5013
	Anhui	0.3888	0.4316	0.4859	0.5056	0.5421	1.6758
	Jiangxi	0.1803	0.2091	0.2616	0.3059	0.3561	3.4615
	Shandong	0.1890	0.2438	0.2616	0.3096	0.3270	2.7790
	Average	0.2408	0.2765	0.3165	0.3551	0.3903	2.4441
Central China	Henan	0.1748	0.1842	0.2337	0.2537	0.3231	3.1193
	Hubei	0.2098	0.2364	0.2815	0.3432	0.3949	3.2129
	Hunan	0.2238	0.2344	0.2904	0.3610	0.4026	2.9795
	Average	0.2028	0.2183	0.2685	0.3193	0.3735	3.1006
Southwest China	Sichuan	0.3082	0.3371	0.4272	0.5168	0.5595	3.0264
MGPAs	Average	0.2247	0.2499	0.2967	0.3449	0.3873	2.7596

describes the LRCC of 13 provinces in China’s MGPAs in 2000, 2005, 2010, 2015, and 2020. The average LRCC score of the MGPAs increased by 72.36% from 2000 to 2020, with an average annual growth rate of 2.7596%. All 13 provinces showed an upward trend in their LRCC scores over the years, indicating that the LRCC of the MGPAs was increasingly higher. Among the 13 provinces, Jiangxi (3.4615%) had the greatest average annual growth rate in LRCC scores, followed as Heilongjiang (3.3603%), Hubei (3.2129%), Henan (3.1193%), Hebei (3.1015%), Sichuan (3.0264%), Hunan (2.9795%), and Shandong (2.7790%). The average annual growth rate in LRCC scores of these eight provinces exceeded those of the MGPAs (2.7596%). Sichuan (0.5595) had the highest score for LRCC in 2020, followed by Anhui (0.5421), Hebei (0.4268), Hunan (0.4026), and Hubei (0.3949). The LRCC scores of these five provinces were above the average LRCC score of the MGPAs (0.3873). In addition, the LRCC scores of Anhui have consistently ranked first since 2000, but it was surpassed by Sichuan in 2012 for the first time due to its smallest average annual growth rate (1.6758%). In 2020, the LRCC scores of Henan surpassed those of Liaoning for the first time, even though Henan ranked last among the 13 provinces from 2000 to 2019.

Regarding the regions, they were ranked in descending order according to their annual growth rate of LRCC scores as Central China (3.1006%), Southwest China (3.0264%), Northeast China (2.8322%), North China (2.6542%), and East China (2.4441%), with the first three regions having the annual growth rate of LRCC scores above the average of MGPAs (2.7596%). The five regions were ranked in descending order of their LRCC scores in 2020 as Southwest China (0.5595), North China (0.3948), East China (0.3903), Central China (0.3735), and Northeast China (0.3346), with the first three regions having LRCC scores above the average of MGPAs (0.3873).

4.2.2. Evolution of LRCC

The evolution of the number of provinces in the five different LRCC classifications is shown in Figure 4. The number of provinces with a great surplus LRCC remained constant at two in the period from 2000 to 2020. The number of provinces with a surplus LRCC remained relatively stable over the years. It remained at three from 2000 to 2009, decreased to two in 2011, after increasing from three to four in 2010, and then remained at three after 2012. The provinces with a balanced or overloaded LRCC experienced the most drastic changes before 2015. After 2015, however, there was no change in the number of provinces in both categories, remaining at two and four, respectively. The number of provinces with seriously overloaded LRCC decreased from three to two in 2003 and increased to three in 2004. During the period from 2005 to 2020, two provinces consistently experienced a serious overload of LRCC.

The LRCC of the 13 provinces statistically plotted over time is shown in Figure 5. The LRCC is relatively low in three provinces in Northeast China, as well as in Jiangsu, Jiangxi, and Henan. Among them, the LRCC of Heilongjiang and Jiangxi showed an overall increasing trend, shifting from the great surplus or surplus categories to the balanced category. Henan and Liaoning experienced relatively less change and were consistently in the great surplus category for an extended period of time. In Jilin and Jiangsu, the LRCC frequently shifted between the balanced and surplus categories. Compared to the six provinces mentioned above, Inner Mongolia, Hebei, Shandong, Hubei, and Hunan had relatively high LRCC. Among these, the LRCC of Hebei shifted from the severely overloaded category to the overloaded category, while the LRCC of Hubei and Hunan instead shifted from the overloaded category to the severely overloaded category. The LRCC of Inner Mongolia and Shandong showed frequent fluctuations, with the former fluctuating back and forth between the balanced, overloaded, and seriously overloaded categories, and the latter fluctuating back and forth between the surplus, balanced, and overloaded categories. In addition, the LRCC in both Anhui and Sichuan remained the highest and unchanged from 2000 to 2020, with both provinces remaining in the seriously overloaded category.

In 2020, the socioeconomic development of Liaoning and Henan was substantially below the carrying capacity of land resources, resulting in both provinces having the LRCC in the great surplus category. The LRCC of Heilongjiang and Shandong were in the surplus category, indicating that the socioeconomic development in both provinces was below the carrying capacity of land resources. The LRCC in the balanced category in Jilin and Jiangsu reflected a balanced approach to socioeconomic development in relation to the carrying capacity of land resources within these two provinces. In Inner Mongolia, Hebei, Hubei, and Hunan, the LRCC fell into the overloaded category, suggesting that their socioeconomic development exceeded the carrying capacity of their land resources. Socioeconomic development far exceeded the carrying capacity of land resources within Anhui and Sichuan, as indicated by the fact that LRCC fell into the seriously overloaded category.

4.3. Coupling Coordination Relationships between the Subsystems of LRCC

4.3.1. Results of the Degree of Coupling Coordination

Table 7. Degree of coupling coordination of LRCC of the MGPAs.

Regions	Provinces	2000	2005	2010	2015	2020	Growth (%)
Northeast China	Heilongjiang	0.3203	0.3738	0.4826	0.5679	0.6312	3.4500
	Liaoning	0.3464	0.3938	0.4568	0.5509	0.6012	2.7950
	Jilin	0.3149	0.3773	0.4329	0.4959	0.5928	3.2136
	Average	0.3272	0.3816	0.4574	0.5382	0.6084	3.1499
North China	Inner Mongolia	0.4008	0.4144	0.4410	0.5192	0.5915	1.9651
	Hebei	0.4163	0.5002	0.5740	0.6705	0.7199	2.7764
	Average	0.4086	0.4573	0.5075	0.5949	0.6557	2.3930
East China	Jiangsu	0.4057	0.4421	0.4950	0.5476	0.6210	2.1514
	Anhui	0.5005	0.5946	0.6653	0.7213	0.7551	2.0775
	Jiangxi	0.3498	0.4061	0.4654	0.5124	0.6096	2.8161
	Shandong	0.2994	0.4573	0.4828	0.5573	0.5939	3.4840
	Average	0.3889	0.4750	0.5271	0.5847	0.6449	2.5611
Central China	Henan	0.3419	0.3823	0.4507	0.5067	0.5667	2.5587
	Hubei	0.3867	0.4175	0.5052	0.5850	0.6534	2.6574
	Hunan	0.4036	0.4394	0.5411	0.6401	0.7026	2.8106
	Average	0.3774	0.4131	0.4990	0.5773	0.6409	2.6832
Southwest China	Sichuan	0.3546	0.4576	0.6045	0.7384	0.8052	4.1857
MGPAs	Average	0.3724	0.4351	0.5075	0.5856	0.6495	2.8202

describes the degree of coupling coordination of the LRCC of the 13 provinces in China’s MGPAs in 2000, 2005, 2010, 2015, and 2020. The average degree of coupling coordination of LRCC of the MGPAs increased by 74.41% from 2000 to 2020, with an average annual growth rate of 2.8202%. All 13 provinces showed an upward trend in their degree of coupling coordination of LRCC over the years, indicating the coupling coordination relationships among the subsystems of the LRCC of the MGPAs was improving. Among the 13 provinces, Sichuan (4.1857%) had the greatest average annual growth rate in degree of coupling coordination of its LRCC, followed by Shandong (3.4840%), Heilongjiang (3.4500%), and Jilin (3.2136%). The average annual growth rate in the degree of coupling coordination of LRCC of these four provinces exceeded those of the average of MGPAs (2.8202%). Sichuan (0.8052) also had the highest degree of coupling coordination of its LRCC in 2020, followed by Anhui (0.7551), Hebei (0.7199), Hunan (0.7026), and Hubei (0.6534). The degree of coupling coordination of these five provinces was above the average of LRCC score of the MGPAs (0.6495). Similar to the LRCC score, the degree of coupling coordination of Anhui has consistently ranked first since 2000, but it was surpassed by Sichuan in 2014 for the first time, due to its relatively low average annual growth rate (2.0775%). The province with the lowest degree of coupling coordination of LRCC changed 12 times, and Jilin was the province with the lowest degree of coupling coordination for the longest period of time (six years). From 2018 to 2020, Henan had the lowest degree of coupling coordination among the 13 provinces.

Regarding the regions, they were ranked in descending order of the annual growth rate of degree of coupling coordination as Southwest China (4.1857%), Northeast China (3.1499%), Central China (2.6832%), East China (2.5611%), and North China (2.3930%), with the first two regions having the annual growth rate of their degree of coupling coordination above the average of MGPAs (2.8202%). The five regions were ranked in descending order of their degree of coupling coordination in 2020 as Southwest China (0.8052), North China (0.6557), East China (0.6449), Central China (0.6409), and Northeast China (0.6084), with the first two regions having a degree of coupling coordination above the average of MGPAs (0.6459).

4.3.2. Evolution of Coupling Coordination Relationships

The evolution of the number of provinces in seven different LRCC coupling coordination relationship classifications is shown in Figure 6. In general, the number of provinces in the disordered recession categories was reduced, and an increasing number of provinces were in the coordinated development categories. Between the years 2000 and 2020, no province experienced an LRCC that fell into the categories of extreme disordered recession, severe disordered recession, or high-quality coordinated development. During the same period, the number of provinces with a LRCC in the categories of mild disordered recession, on the verge of disordered recession, barely coordinated development, primary coordinated development, and intermediate coordinated development showed great fluctuations. In 2004, the last province (Jiangxi) that had previously been categorized as experiencing moderate disordered recession in LRCC shifted to the mild disordered recession category. The LRCC of one province (Anhui) entered the primary coordinated development category for the first time in 2007. In 2012, the final province (Jiangxi) with LRCC in the mild disordered recession category, transitioned to being on the verge of the disordered recession category. Two provinces (Anhui and Sichuan) entered the intermediate coordinated development category based on their LRCC coupling coordination degree in 2013. In 2016, no provinces had LRCC in the disordered recession category. The first province (Sichuan) with LRCC in the well-coordinated development category appeared in 2020.

Figure 7 illustrates the LRCC coupling coordination relationship of the 13 provinces statistically plotted from 2000 to 2020. At the beginning of the 21st century, Anhui was the only province where the degree of coupling coordination of LRCC subsystems was in the coordinated development category. In contrast, the degree of coupling coordination of LRCC subsystems in the other 12 provinces remained in the disordered recession category. From 2000 to 2020, continued improvement in the coupling coordination relationship of LRCC subsystems was observed in three provinces in Northeast China, as well as in Hebei, Jiangsu, Anhui, Henan, and Sichuan. During the same period, fluctuations in the coupling coordination relationship of LRCC subsystems in Inner Mongolia, Jiangxi, Shandong, Hubei, and Hunan were relatively more frequent, with several setbacks within a positive process overall.

In 2020, the LRCC subsystems of all 13 provinces in the MGPAs were in the coordinated development category. The LRCC subsystems of Sichuan exhibited the highest degree of coupling coordination and was categorized as well-coordinated development. One province each in North China, East China, and Central China, namely Hebei, Anhui, and Hunan, displayed a relatively high degree of coupling coordination in their LRCC subsystems, classifying them under the category of intermediate coordinated development. The degree of coupling coordination of LRCC subsystems in Heilongjiang, Liaoning, Jiangsu, Jiangxi, and Hubei was relatively low, and they were categorized in the primary coordinated development category. One province each in Northeast China, North China, East China, and Central China, namely Jilin, Inner Mongolia, Shandong, and Henan exhibited the lowest degree of coupling coordination in their LRCC subsystems, which were classified under the barely coordinated development category.

4.4. Obstacle Factors for LRCC

4.4.1. Obstacle Subsystems and Obstacle Degree of LRCC

Table 8. Obstacle subsystems and obstacle degree of LRCC of the MGPAs from 2000 to 2020.

Year	Driver		Pressure		State		Impact		Response
	Ranking	Obstacle Degree (%)	Ranking	Obstacle Degree (%)	Ranking	Obstacle Degree (%)	Ranking	Obstacle Degree (%)	Ranking	Obstacle Degree (%)
2000	4	12.22	5	3.38	1	56.35	2	14.06	3	14.00
2001	4	12.15	5	3.41	1	56.43	3	13.99	2	14.01
2002	4	12.10	5	3.46	1	56.73	3	13.85	2	13.86
2003	4	11.84	5	3.50	1	56.93	2	14.12	3	13.61
2004	4	11.71	5	3.65	1	56.94	3	13.79	2	13.91
2005	4	11.71	5	3.73	1	56.86	3	13.80	2	13.91
2006	4	11.74	5	3.85	1	56.67	3	13.77	2	13.97
2007	4	11.56	5	4.01	1	56.61	2	13.98	3	13.84
2008	4	11.49	5	4.08	1	56.86	2	13.85	3	13.72
2009	4	11.58	5	4.12	1	56.92	2	13.86	3	13.53
2010	4	11.25	5	4.21	1	57.11	2	13.82	3	13.61
2011	4	10.99	5	4.48	1	57.34	3	13.56	2	13.63
2012	4	10.88	5	4.50	1	57.62	2	13.57	3	13.44
2013	4	10.82	5	4.56	1	57.87	2	13.48	3	13.26
2014	4	10.59	5	4.32	1	58.81	2	13.25	3	13.03
2015	4	10.59	5	4.34	1	58.84	2	13.26	3	12.97
2016	4	10.66	5	4.09	1	58.79	2	13.27	3	13.19
2017	4	10.59	5	4.00	1	58.85	2	13.35	3	13.22
2018	4	10.42	5	4.00	1	59.14	2	13.58	3	12.86
2019	4	10.25	5	3.98	1	59.69	2	13.59	3	12.50
2020	4	9.39	5	4.24	1	59.61	2	13.60	3	13.16

shows the obstacle subsystems and the obstacle degree of LRCC. From 2000 to 2020, an overall increase in the obstacle degree of the state and pressure subsystems was observed, indicating a gradual rise in the influence of these two subsystems on the balanced development of LRCC. Conversely, there was a general reduction in the obstacle degree within the driver, impact, and response subsystems, indicating their diminishing constraints on LRCC.

The state subsystem exhibited the highest obstacle degree among the five subsystems, emerging as the primary factor limiting the balanced development of the LRCC in the MGPAs. The response and impact subsystems had closely comparable obstacle degrees, securing the second and third positions since 2012. In specific years like 2001, 2002, 2004 to 2006, and 2011, the obstacle degree of the response subsystem exceeded that of the impact subsystem, but this trend was reversed in other years. The driver and pressure subsystems displayed the lowest obstacle degrees, ranking fourth and fifth, respectively, indicating a relatively minor impact of these two subsystems on the LRCC.

4.4.2. Obstacle Indicators and Obstacle Degree of LRCC

The top five obstacle indicators and the obstacle degree of LRCC (cumulative obstacle degree over 50%) were selected as the main obstacle indicators, which are listed in

Table 9. Top five obstacle indicators and obstacle degree of LRCC of the MGPAs from 2000 to 2020.

Year	First Place		Second Place		Third Place		Fourth Place		Fifth Place
	Obstacle Indicator	Obstacle Degree (%)	Obstacle Indicator	Obstacle Degree (%)	Obstacle Indicator	Obstacle Degree (%)	Obstacle Indicator	Obstacle Degree (%)	Obstacle Indicator	Obstacle Degree (%)
2000	S1	14.27	S9	12.70	I5	12.52	D6	6.23	S3	5.41
2001	S1	14.30	S9	12.77	I5	12.51	D6	6.21	S3	5.42
2002	S1	14.37	S9	12.77	I5	12.46	D6	6.20	S3	5.44
2003	S1	14.35	S9	13.38	I5	12.62	D6	6.14	S3	5.43
2004	S1	14.65	I5	12.56	S9	12.55	D6	6.06	S3	5.54
2005	S1	14.80	I5	12.54	S9	12.28	D6	6.01	S3	5.59
2006	S1	14.92	I5	12.39	S9	11.75	D6	5.97	S3	5.64
2007	S1	15.23	I5	12.64	S9	10.93	D6	5.82	S3	5.75
2008	S1	15.52	I5	12.53	S9	10.43	S3	5.85	D6	5.68
2009	S1	15.62	I5	12.64	S9	10.44	S3	5.80	D6	5.61
2010	S1	15.82	I5	12.56	S9	10.22	S3	5.89	D6	5.54
2011	S1	16.16	I5	12.50	S9	9.92	S3	6.06	S4	5.48
2012	S1	16.27	I5	12.43	S9	9.67	S3	6.25	S4	5.56
2013	S1	16.60	I5	12.30	S9	8.97	S3	6.43	S4	5.79
2014	S1	16.46	I5	12.11	S9	10.10	S3	6.12	S4	5.22
2015	S1	16.73	I5	12.05	S9	9.56	S3	6.23	S4	5.76
2016	S1	16.69	I5	12.12	S9	9.63	S3	6.29	S4	5.83
2017	S1	16.57	I5	12.15	S9	9.64	S3	6.32	S4	5.88
2018	S1	16.65	I5	12.37	S9	9.47	S3	6.41	S4	5.98
2019	S1	17.03	I5	12.52	S9	8.64	S3	6.55	S4	6.10
2020	S1	17.06	I5	12.50	S9	8.69	S3	6.57	S4	6.15

Note: D6 = value of agricultural output per capita; S1 = arable land area per capita; S3 = GDP growth rate; S4 = GDP per capita; S9 = grain production yield per unit area of arable land; I5 = grain production yield per capita.

. The main obstacle indicator was S1 (arable land area per capita), with the largest obstacle degree of the 42 indicators, which greatly impacted the balanced development of LRCC. And the overall rise in the obstacle degree of S1 indicates its increasing impact on LRCC. The obstacle degrees of S9 (grain production yield per unit area of arable land) and I5 (grain production yield per capita) were close to each other, with S9 holding the second position and I5 the third position until 2003. Subsequently, the decrease in the obstacle degree of S9 led to a reversal in their rankings, with I5 surpassing S9 after 2005. The declining obstacle degree of D6 (value of agricultural output per capita) led it to drop from fourth place among 42 indicators to fifth place in 2008, and it ceased to be among the top five obstacle indicators in 2011. The growing obstacle degree of S3 (forest land area per capita) elevated its ranking from the fifth position among obstacle indicators to the fourth in 2008. S4 (grassland area per capita) emerged as the fifth-ranked obstacle indicator since 2011, with a general increase in obstacle degree.

5. Discussion

5.1. Evaluation of Indicator Weights

Weights are measures of indicator importance, and reasonable weights can improve the accuracy of LRCC evaluation results [81]. For MGPAs, there is an objective conflict between economic development needs and food production. Population growth and urbanization have led to a continuous increase in total food demand, as well as the expansion of urban construction land, the conversion of arable land into construction land, and the exploitation of forests for agriculture and horticulture [36]. The rapid changes in land resource status and utilization structure affect the regional resource carrying capacity status and constrain sustainable social and economic development. The state subsystem covers important information such as arable land area per capita, construction land area per capita, and grain production yield per unit area of arable land, indicating the current land resource status and representing the level of land health and sustainable development [44,45,48]. It should be the primary subsystem determining the changes in LRCC in the MGPAs. As the foundation of human survival, food should be the primary constraint on LRCC [82]. As the most fundamental element of food production, arable land should be a key indicator reflecting LRCC [83].

5.2. Evaluation of LRCC

During the period from 2000 to 2020, there was an overall positive trend in the coordinated development of the economy, ecology, and society in China’s MGPAs. The condition of LRCC improved progressively, and there was a gradual increase in grain production levels. Land serves as the primary carrier for grain production and is a crucial pillar of the national food security strategy. In light of this, the government has introduced relevant policies to promote the development of LRCC in the MGPAs. In 2003, the Central Committee explicitly proposed the “strictest farmland protection system”. In 2015, the State Council issued the “General Plan for Ecological Civilization System Reform”, which explicitly emphasized the need to embrace the concept of spatial equilibrium and to maintain a balance between population, economy, the resource environment, and development. It stressed that population size, industrial structure, and growth rate should not exceed the carrying capacity of local land resources and environmental capacity. Due to the comprehensive impact of government-led ecological conservation and restoration projects, advancements in agricultural technology, increased environmental awareness, and optimized land utilization, the LRCC in the MGPAs has been steadily increasing year by year [81,84]. However, rapid population growth and urbanization have led to a tense relationship between human resources and land resources, and the MGPAs still face significant pressures on land resources [81]. According to research by Luo, from 2001 to 2017, the overall pressure on arable land in the MGPAs showed a downward trend, but in 2017, eight provinces in the MGPAs still had arable land pressure at or above the critical pressure zone [85]. Currently, the overall LRCC in the MGPAs is in an unbalanced development state, with overloading as the main feature.

When viewed from a spatial dimension, the overall spatial distribution of LRCC in the MGPAs displays a characteristic pattern of “higher in the south, lower in the north”. The evaluation of LRCC requires a comprehensive examination that takes into account the supply capacity of land resources as well as the pressure generated by human activities [47,72,86]. Provinces with a surplus of LRCC tend to be concentrated in the less densely populated northern regions, such as Heilongjiang and Liaoning. These areas have large land areas, lower population densities, relatively higher arable land areas per capita, and thus exert less pressure on land resources. Furthermore, these regions tend to have relatively lower chemical fertilizer usage per unit area of arable land and pesticide usage per unit area of arable land, leading to reduced pollution pressure on land resources caused by agricultural and pastoral activities. Conversely, provinces with LRCC overload tend to be concentrated in the southern regions with large populations or rapid economic development, such as Anhui and Sichuan. These areas have limited exploitable land resources, experience fast population growth, resulting in higher population densities and relatively lower arable land areas per capita. Moreover, these regions exhibit relatively higher levels of urbanization, leading to strong demand for non-agricultural land use, which has caused their socioeconomic development to exceed the carrying capacity of land resources.

5.3. Evaluation of Coupling Coordination Relationships

Given that land in the MGPAs not only serves as a vital underpinning for economic growth but also functions as a primary vehicle for ecological governance, uncovering the degree of coupling coordination between LRCC and environment is a focal point for achieving high-quality, sustainable development in the MGPAs [45]. The degree of coupling coordination can be utilized to analyze the level of coordinated development in the ecological economy of the MGPAs, revealing the extent of mutual influence among subsystems of LRCC [45]. From 2000 to 2020, the level of coupling coordination in China’s MGPAs has been steadily increasing, with the rate of growth gradually on the rise. Notably, in 2020, there was a reduction in the gap of coupling coordination levels. This may be attributed to the formidable task of achieving further upgrades once a certain threshold of coordination development has been reached [87,88]. Over the next few years, there may be a deceleration in the rate of increase of coupling coordination levels, potentially resulting in a further reduction in the coordination gap. In terms of the quantity of coupling coordination types, with the passage of time, provinces falling under the category of uncoordinated decline have gradually transitioned toward coordinated development. It was not until 2016 that all MGPAs entered the realm of coordinated development. The interactions among subsystems of LRCC are steadily intensifying, progressing towards a more organized direction. This may be attributed to increased attention from local governments on promoting coordinated development [88,89].

From a spatial perspective, the benign interactions among subsystems in the MGPAs exhibit an overall “southern high and northern low” spatial distribution, with the spatial distribution of coupling coordination types shifting from dispersion to aggregation. This may be attributed to the rapid economic development in the southern regions, despite issues such as continuously increasing population density, expanding construction land scale, and the concentration of urban populations. These areas have made significant investments in public resources, infrastructure, and environmental protection, which have facilitated the formation of positive human–environment interactions and ultimately led to a continuous promotion of the coordination between land resources and the environment, resulting in high-quality coordinated development [45,88,90]. In contrast, the northern regions, despite their favorable ecological environment, experience slower economic development, relatively weak infrastructure construction, and a trend of population migration towards more advantageous areas [88,89,90]. Furthermore, in these regions, most industries are undergoing a slow transition, and the demand for construction land from the industrial structure is robust. This has led to a mismatch between population and land urbanization development, resulting in poor overall coordination within the LRCC subsystems [89].

5.4. Diagnosis of Obstacle Factors

Identifying obstacle factors and key obstacle indicators aids in analyzing the critical factors influencing LRCC and provides targeted policies and recommendations to promote the future development of MGPAs [72]. With the rapid development of the economy and society, MGPAs continuously adjust and improve relevant land-use policies, while enhancing the protection of the ecological environment. This has led to a year-by-year increase in the alignment between the driving force subsystem, impact subsystem, response subsystem, and the demand for LRCC, resulting in a decreasing trend in obstacle degree. Indicators of the state subsystem and pressure subsystem, as direct reflections of the ecological environment of land resources, have become the primary obstacles affecting LRCC, especially with the accelerated urbanization process. Therefore, in the context of the new economic normal, if MGPAs aim to enhance LRCC, they should focus on the pressure and state subsystems, adopt more scientifically rational land planning, increase efforts in environmental protection, alleviate pressure on land resources, and improve the state of land resources [47,72,91].

The key obstacle indicators influencing the LRCC in MGPAs include arable land area per capita, grain production yield per unit area of arable land, grain production yield per capita, forest land area per capita, and grassland area per capita. Among these five indicators, four belong to the state subsystem, providing further evidence of the state subsystem’s impact on LRCC. In various years, the primary obstacles are often arable land area per capita, grain production yield per unit area of arable land, and grain production yield per capita, indicating the significant influence of human activities on land-based grain production on land ecology. The mismatch between population and arable land area results in a pressing food demand [72,91]. Additionally, indicators that directly affect LRCC, such as forest land area per capita and grassland area per capita, have also become significant obstacles, necessitating further efforts to protect forests and grasslands.

6. Conclusions and Suggestions

Based on the DPSIR framework, this study comprehensively evaluated the LRCC of China’s MGPAs from 2000 to 2020. The main conclusions are summarized as follows:

(1) The LRCC shows an upward trend, but the land still faces significant pressure, reflected in the predominance of overloaded LRCC. Attributed to factors such as population density and levels of economic development, the overall LRCC exhibits a spatial differentiation feature of “higher in the south and lower in the north”.

(2) The degree of coupling coordination shows an upward trend, but the magnitude of increase is decreasing, mainly categorized as primary coordination and intermediate coordination. The benign interactions within the subsystems manifest an overall spatial distribution of “higher in the south and lower in the north”, with a narrowing gap in coordination levels among provinces.

(3) The state subsystem and the pressure subsystem are the main obstacles influencing the LRCC in the MGPAs. Arable land area per capita, grain production yield per unit area of arable land, grain production yield per capita, forest land area per capita, and grassland area per capita are the main indicators hindering LRCC.

Based on the above research findings, several policy implications can be provided for the high-quality development of MGPAs:

(1) The LRCC in the MGPAs still has significant potential for improvement. For provinces with surplus LRCC, population growth can be appropriately promoted, while ensuring the safety of land resources, moderate development of land resources can be carried out, the current land utilization level can be improved, and the land carrying potential can be enhanced by increasing the efficiency of land resource utilization. For provinces with overloaded LRCC, strict land use control and an urban land scale auditing system should be implemented. Arable land resources should be strictly protected, and the conversion of agricultural land to construction land should be controlled. Making full use of the potential of existing construction land is the main direction of urban land use.

(2) The level of coupling coordinated of LRCC in the MGPAs needs to be further improved. Based on regional resource advantages, efforts should be made to optimize industrial structure, promote economic development, and increase investment in infrastructure and environmental protection.

(3) Optimizing agricultural production technology and management to increase the grain yield per unit of arable land, thus boosting overall grain production and alleviating pressure on arable land resources. It is necessary to establish a dynamic monitoring system for the protection of forests and grasslands to prevent excessive deforestation and overgrazing.

This study employed the DPSIR framework to comprehensively assess the LRCC in China’s MGPAs and provided targeted recommendations. However, there are still two limitations that need to be addressed. First, while this study utilized the AHP-EW method to ensure both practical significance and objectivity of the weightings, quantitative validation of this method’s reliability was not conducted. Nevertheless, the AHP-EW method has been widely applied in previous research and is considered relatively mature. Therefore, it remains a useful tool. Second, this study used province-level data, which may overlook the heterogeneity of county-level LRCC, degree of coupling coordination, and obstacle factors. Future research could consider incorporating more detailed county-level data for a more refined analysis.