Ecological Security Assessment and Warning of Cultivated Land Quality in the Black Soil Region of Northeast China

Abstract

The ecological security of cultivated land critically depends on maintaining the quality of the land under cultivation. For the security of the nation’s grain supply, the evaluation and early warning of cultivated land quality (CLQ) are essential. However, previous studies on the assessment of the ecological safety of CLQ only rigidly standardized the assessment indicators and failed to investigate the positive and negative trends and spatiotemporal driving factors of the indicators. The main objective of this study was to develop a drive–pressure–state–response (DPSR) model to identify the hierarchical structure of indicators, using an improved matter–element model to assess the CLQ in the black soil region of northeastern China from 2001 to 2020. A panel data model was employed to explore the crucial drivers of CLQ warnings. The findings reveal that socioeconomic development has a potential impact on the improvement of CLQ. CLQ is generally in a secure state, with 69.71% of cities with no warnings and only 3.46% and 0.13% of cities under serious and extreme warnings, respectively. Compared with 2001, the CLQ in 2020 effectively improved by socioeconomic development and the conservation and reasonable utilization of arable land. According to the early warning results, the cultivated land in the northern regions was of higher quality than that in the southern regions. Moreover, the CLQ was significantly positively correlated with the agricultural GDP growth rate, grain yield per unit of cultivated land area, annual precipitation, and the habitat quality index, and was significantly negatively correlated with land carrying capacity. The findings of this study can provide a scientific and targeted basis for black soil conservation and utilization.

1. Introduction

Cultivated land is a valuable and non-renewable resource that contributes to ecological security, food security, and social development, in addition to offering various critical services and products [1]. In addition to providing food and agricultural products for human consumption, cultivated land supplies valuable agricultural raw materials for industrial development, enhancing industrialization in the country [2]. Changes in cultivated land quality (CLQ) play an essential role in regional grain security [3]. Owing to its high levels of organic matter and nutritional potential, black soil is one of the most precious resources with respect to cultivated land [4]. The black soil region in northeast China, one of China’s major grain-producing areas, has been negatively affected by long-term overexploitation and utilization, which poses a potential risk to national grain security [5]. The black soil region of northeast China contributes 10.8% of the nation’s overall grain production and is crucial to national ecological security and grain security [6,7,8]. Therefore, ecological security assessments and CLQ warnings in the black soil region are imperative in terms of maintaining grain security and sustainable development.

CLQ is a comprehensive representation of cultivated land conditions, soil fertility, and soil quality with the aim of meeting human’s well-being requirements [9]. The ecological security of CLQ refers to a healthy and unthreatened state of the natural conditions of arable land and the quantity of arable land and soil [10]. In addition to the influence of natural factors, such as soil degradation and extreme precipitation, CLQ is also influenced by social factors, such as economic development and human activities [11,12]. Moreover, social factors can be used to quantitatively assess the ecological security of CLQ through time-series data [13,14]. Urban expansion has been driven by the rapid urbanization process, which is gradually endangering arable land in peri-urban regions and reducing the amount of arable land [15,16]. The abuse of pesticides and chemical fertilizers to achieve high levels of food production can directly result in the overuse of cultivated land resources [17]. Ecological environmental variables, such as precipitation, vegetation cover, and habitat quality, also play an essential role in the ecological security assessment of CLQ. Precipitation and vegetation cover are important determinants of soil erosion rates and natural drivers of declines in CLQ [18,19,20]. Furthermore, these socioeconomic and ecological environmental variables have both temporal and spatial attributes, which contribute to the long-term dynamic monitoring of the ecological security of CLQ [21,22,23]. Therefore, in contrast to soil attribute variables, the ecological security of CLQ can be assessed based on socioeconomic development and the ecological environment to explore the spatiotemporal dynamics of CLQ.

The matter–element model [24] is a multicriteria decision analysis (MCDA) model proposed by a Chinese scholar that transforms the incompatibility problem into a multiparameter compatibility problem through the use of a correlation function. Furthermore, the matter–element model primarily uses correlation degrees to measure how closely indicators and results conform to each other when evaluating the study’s objectives in qualitative and quantitative ways. Matter–element models have been used to solve a variety of practical issues in many fields, including groundwater health monitoring, urban risk assessment, and river health monitoring [25,26]. Such models were initially applied in evaluation studies to assess, e.g., land suitability [27], river health [28], and dynamic water quality [29]. However, there is a lack of studies incorporating matter–element models that are applicable to CLQ assessment and early warning. Additionally, in the existing matter–element models, the correlation function is insufficiently developed and optimized, and the influence of different indicator types on the correlation function is ignored, leading to significant variations between the theoretical and practical results.

Existing studies on the ecological security assessment of CLQ have focused on the effects of soil attributes such as soil quality, soil texture, and soil nutrient information on CLQ [30,31]. However, the influences of monitoring technology and the environment make it difficult to obtain time-series soil profile data, especially in large-scale areas, and it is impossible to accurately explore the spatiotemporal dynamics of CLQ using soil profile data only. Additionally, the method used to assess the ecological security of CLQ is too homogeneous in its treatment of indicators and it only standardizes the indicators in a fixed way, without fully considering the role of positive indicators in promoting the assessment objectives and the role of negative indicators in hindering the assessment objectives.

In this context, the objective of this study was to explore the spatiotemporal dynamics of the ecological security of CLQ based on socioeconomic development and to focus on the effects of positive and negative indicators on CLQ. We developed a drive–press–state–response model to determine the levels of fourteen indicators covering the four dimensions of social economy, human activities, cultivated land attributes, and ecological environment. Then, we applied the enhanced matter–element model to the assessment and early CLQ warning in the black soil region of northeast China and analyzed the spatial effect and critical driving factors of CLQ from 2001 to 2020. The panel data model was utilized to analyze the spatiotemporal driving factors affecting CLQ. The aim of this study was to assess and provide early warnings of CLQ trends resulting from socioeconomic and ecological factors and to analyze the critical drivers in order to address potential security threats with respect to cropland ecological security and black soil conservation.

2. Materials and Methods

2.1. Study Area

The black soil region of northeast China (118°35′ E~135°05′ E, 38°43′ N~53°33′ N) includes 40 cities in the provinces of Jilin, Liaoning, and Heilongjiang, as well as 4 cities in the eastern Inner Mongolia Autonomous Region (Figure 1 and Appendix A). The average annual temperature of the whole region is −7 to 11 °C; the winter is cold and long; and the annual precipitation in most areas is 450 to 850 mm. Due to the influence of the monsoon climate, the rainy season in the northeastern black soil area is mainly concentrated in the hot summer. This characteristic of precipitation and long-term sunshine has an essential impact on crop growth and yield.

2.2. Indicator System

We selected various indicators from existing studies to assess the ecological safety of CLQ and developed a multilevel indicator model that includes the four levels of socioeconomic drivers, human activity pressure, ecological state, and cultivated land response (DPSR), with a total of fourteen indicators [32,33]. These indicators have been widely used for ecological security assessments and they have strong logical relationships [34,35]. Moreover, we selected indicators with time-series properties to explore the spatiotemporal dynamics of CLQ from 2001 to 2020. The weights of each indicator were calculated using the analytic hierarchy process and the entropy method, and the calculation results passed the test (

Table 1. Indicator information and weights.

Criterion Layer	Index Layer	Indicator Type	Weights
Socioeconomic Drivers	Urbanization intensity (A1)	−	0.0521
	Greenery coverage of built-up areas (A2)	−	0.0172
	Percentage of secondary industry GDP (A3)	−	0.0326
	Agricultural GDP growth rate (A4)	+	0.0691
Human Activity Pressure	Population growth rate (B1)	−	0.0306
	Land carrying capacity (B2)	−	0.0539
	Amount of chemical fertilizer per unit of cultivated land area (B3)	−	0.0661
Ecological State	Annual precipitation (C1)	+	0.0726
	Normalized difference vegetation index (C2)	+	0.1066
	Habitat quality index (C3)	+	0.0496
Cultivated Land Response	Cultivated land area growth rate (D1)	+	0.0867
	Replanting index (D2)	+	0.1307
	Land reclamation rate (D3)	+	0.0890
	Grain yield per unit of cultivated land area (D4)	+	0.1432

Note: “+”, positive indicator; “−”, negative indicator.

).

2.3. Data Sources

The economic data used in this study were mainly retrieved from provincial statistical yearbooks and the China Statistical Yearbook (accessed on 1 March 2022). Land-use data were obtained from the time-series dataset published by Yang and Huang [36] (https://zenodo.org/record/5210928 (accessed on 1 March 2022)). Meteorological data were retrieved from the precipitation dataset of the Resource and Environment Science and Data Center (https://www.resdc.cn/ (accessed on 1 March 2022)). The maximum synthesis method was used to fuse the normalized difference vegetation index (NDVI) derived from Landsat images from May to September of every year in the study period.

2.4. Improved Matter–Element Model for Early Warning

The matter–element model can describe the quantitative and qualitative changes in indicators and is suitable for multiple indicator assessments [37]. However, in real-world scenarios, the correlation function of the matter–element model needs to be improved and adjusted according to the type of indicator. The CLQ early warning results obtained using the original matter–element model can be significantly unrelated to the actual situation when the effects of positive and negative indicators on early warning are ignored. To address this issue, in this study, we introduced fuzzy sets into the original model, which were used to establish a correlation function of the CLQ early warning model that matches practical applications.

(1)

Establish the early warning matter–element of CLQ.

The early warning matter–element of CLQ is represented by an ordered ternary array $R=(N,c,v)$, where $N$ is the level of the early warning, $c_i$ is the indicator feature of the early warning, and $v_i$ is the value of $c_i$.

$$R=[N,c,v]=\left[\begin{array}{ccc}N& c_1& v_1\\ & c_2& v_2\\ & ⋮& ⋮\\ & c_n& v_n\end{array}\right]$$

(1)

(2)

Determine the classical and limited fields.

The classical field ($R_j$) indicates the range of values for the city indicator ($c_i$) under different early warning levels. The limited field ($R_p$) indicates the range of values for the city indicator ($c_i$) under all early warning levels.

$$R_j=(N_j,c_i,X_{ji})=\left[\begin{array}{ccc}{N}_j& c_1& (a_{j1},b_{j1})\\ & c_2& (a_{j2},b_{j2})\\ & ⋮& ⋮\\ & c_n& (a_{jn},b_{jn})\end{array}\right]$$

(2)

Here, $N_j$ is the $j$th early warning level, and $X_{ji}$ is the range of values of the $i$th early warning indicator at the $j$th early warning level.

$$R_p=(p,c_i,X_{pi})=\left[\begin{array}{ccc}p& c_1& (a_{p1},b_{p1})\\ & c_2& (a_{p2},b_{p2})\\ & ⋮& ⋮\\ & c_n& (a_{pn},b_{pn})\end{array}\right]$$

(3)

Here, $p$ is the total number of early warning levels, and $X_{pi}$ is the range of values of the $i$th indicator.

(3)

Establish the correlation function based on fuzzy side distance.

The original matter–element model uses a primary correlation function to calculate the correlation degree (Equation (4)). However, this correlation function ignores the effect of positive and negative indicators of the city’s early warning. For example, the range of the second early warning level of the positive indicator $c_1$ is $[a_{21},b_{21}]$. Whereas the primary correlation function is used to calculate the correlation degree, it is assumed that when ${value}_{v_{21}}=\frac{a_{21}+b_{21}}{2}$, $K\left(x\right)$ obtains the maximum value. However, this is inconsistent with the fact that not all indicators achieve the optimal solution at the center point.

$$K(x)=\left\{\begin{array}{cc}-\frac{\rho (x,X_{ji})}{\left|X_{ji}\right|}& x\in X_{ji}\\ \frac{\rho (x,X_{ji})}{\rho (x,X_{pi})-\rho (x,X_{ji})}& x\notin X_{ji}\end{array}\right.$$

(4)

$$\rho (x,X_i)=\left|x-\frac{a_i+b_i}{2}\right|-\frac{b_i-a_i}{2}$$

(5)

Here, $\rho (x,X_i)$ is the extensible distance between value $x$ of the $i$th indicator and the interval $X_i$, and $a_i$ and $b_i$ are the range of interval $X_i$.

To address this issue, in this study, we introduced the fuzzy set $\widehat{A}$ to represent the optimal point ${\widehat{x}}_{0ji}$, with different correlation functions established for positive and negative indicators. The fuzzy side distance correlation function for positive indicators is as follows:

$$\widehat{A}=\left\{K({\widehat{x}}_{0ji})=\max (K(x))\left|x={\widehat{x}}_{0ji}\right.\right\}$$

(6)

$$K(x)=\left\{\begin{array}{cc}-\frac{{\rho }_R(x,{\widehat{x}}_{0ji},X_{ji})}{\left|X_{ji}\right|}& x\in X_{ji}\\ \frac{{\rho }_R(x,{\widehat{x}}_{0ji},X_{ji})}{{\rho }_R(x,{\widehat{x}}_{0ji},X_{pi})-{\rho }_R(x,{\widehat{x}}_{0ji},X_{ji})}& x\notin X_{ji}\end{array}\right.$$

(7)

Here, ${\rho }_R(x,x_{0ji},X_{ji})$ is the right-side distance between $x$ and the interval $X_{ji}$ when $x$ achieves the optimal solution at ${\widehat{x}}_{0ji}$. The right-side distance is calculated according to Equation (8).

$${\rho }_R(x,x_0,X)=\left\{\begin{array}{ll}x-b\hfill & x\ge b\hfill \\ a-x\hfill & x\le x_0\hfill \\ \frac{a-x_0}{b-x_0}(b-x)\hfill & x\in 〈x_0,b〉\hfill \end{array}\right.$$

(8)

Here, $x_0$ is the optimal point in interval $X$, and interval $X$ ranges from $a$ to $b$. In Equation (7), ${\widehat{x}}_{0ji}$ is calculated using the linearly increasing trapezoidal affiliation function (Equation (9)).

$$F(x,X)=\left\{\begin{array}{cc}0& x\le a\\ \frac{x-a}{b-a}& a<x<b\\ 1& x\ge b\end{array}\right.$$

(9)

Similarly, for the negative indicators, we established a fuzzy side distance correlation function and a calculation function of the left-side distance (${\rho }_L(x,x_{0ji},X_{ji})$) and calculated the optimal point (${\widehat{x}}_{0ji}$) using the linearly decreasing trapezoidal affiliation function.

$$K(x)=\left\{\begin{array}{cc}-\frac{{\rho }_L(x,{\widehat{x}}_{0ji},X_{ji})}{\left|X_{ji}\right|}& x\in X_{ji}\\ \frac{{\rho }_L(x,{\widehat{x}}_{0ji},X_{ji})}{{\rho }_L(x,{\widehat{x}}_{0ji},X_{pi})-{\rho }_L(x,{\widehat{x}}_{0ji},X_{ji})}& x\notin X_{ji}\end{array}\right.$$

(10)

$${\rho }_L(x,x_0,X)=\left\{\begin{array}{ll}a-x\hfill & x\le a\hfill \\ x-b\hfill & x\ge x_0\hfill \\ \frac{b-x_0}{a-x_0}(x-a)\hfill & x\in 〈a,x_0〉\hfill \end{array}\right.$$

(11)

$$F(x,X)=\left\{\begin{array}{cc}1& x\le a\\ \frac{b-x}{b-a}& a<x<b\\ 0& x\ge b\end{array}\right.$$

(12)

(4)

Calculate the single dominance index.

After defining the fuzzy side distance correlation function, we can calculate the single dominance index ($k_s\left(x\right)$) based on the correlation degree ($K\left(x\right)$). The dominance index is obtained by normalizing the correlation degree, which makes the matter–element model more objective in representing the early warning results.

$$k_s(x)=\left\{\begin{array}{cc}\frac{K(x)}{\underset{x\in X_0}{\max }K(x)}& K(x)\ge 0\\ \frac{K(x)}{\underset{x\notin X_0}{\max }\left|K(x)\right|}& K(x)<0\end{array}\right.$$

(13)

(5)

Calculate the integrated dominance index.

The indicator weights (${\alpha }_i$) are calculated by the entropy method and the analytic hierarchy process. According to the single dominance index ($k_s\left(x\right)=\left[k_s\right(x_1),k_s(x_2),k_s(x_3),\dots ,k_s(x_N\left)\right]$), the integrated dominance index of the $j$th warning level ($C_j\left(x\right)$) is calculated according to Equation (14).

$$C_j(x)={\displaystyle \sum _{i=1}^N{\alpha }_i{k}_{sj}(x_i)}$$

(14)

2.5. Panel Data Model for Driving Factors of CLQ

Panel data are a type of data with multiple attributes that contain both cross-sectional and time-series information. There are three main types of panel data models—random effects (RE) models, fixed effects (FE) models, and mixed effects (ME) models—that are suitable for analyzing different individual influences [38,39,40]. In this study, we developed a panel data model to investigate the relationship between the dominance index (DI) and different driving factors by considering the time effect and individual effect of early warning targets. Unit root and co-integration tests were performed on the original index data to verify the stationary state and cointegration of the data. Then, we constructed a panel data model to analyze the driving factors of CLQ early warnings.

$$D{I}_{iy}={\displaystyle \sum _{j=1}^{14}{\beta }_j{x}_{jiy}}+{\mu }_{iy}+{\xi }_{iy}+{\epsilon }_{iy}$$

(15)

Here, $i(i=1,2,\dots ,39)$ is the code for warning city $i$; $y(y=2001,2002,\dots ,2020)$ is time information $y$; $D{I}_{iy}$ is the integrated dominance index of CLQ; $x_j$ represents the 14 explanatory variables; ${\beta }_j$ is the parameter of variable $x_j$; ${\mu }_{iy}$ and ${\xi }_{iy}$ are the intercept terms of the individual effect and the time effect, respectively; and ${\epsilon }_{iy}$ is the random disturbance term.

3. Experiments and Results

3.1. CLQ Early Warning Results

3.1.1. The Classical Field and Limited Field of CLQ

The classical domain ($R_j$) and the nodal domain ($R_p$) are calculated according to Equations (1)–(3) using the integrated weight method and the natural breakpoint method (

Table 2. The classical field and limited field.

Indicator	Classical Field					Limited Field
	I	II	III	IV	V
	No Warning	Light Warning	Moderate Warning	Severe Warning	Extreme Warning
A1	(23.24, 50.32)	(50.32, 59.18)	(59.18, 68.1)	(68.1, 76.13)	(76.13, 93.71)	(23.24, 97.71)
A2	(8.78, 33.96)	(33.96, 37.86)	(37.86, 42)	(42, 50.62)	(50.62, 78.33)	(8.78, 78.33)
A3	(11.32, 40.84)	(40.84, 48.05)	(48.05, 56.3)	(56.3, 70.66)	(70.66, 88.69)	(11.32, 88.69)
A4	(10.64, 54.37)	(3.38, 10.64)	(−4.32, 3.38)	(−17.85, −4.32)	(−38.62, −17.85)	(−38.62, 54.37)
B1	(−5.78, −0.22)	(−0.22, 0.19)	(0.19, 0.66)	(0.66, 1.52)	(1.52, 3.14)	(−5.78, 3.14)
B2	(8.84, 136.41)	(136.41, 200.67)	(200.67, 294.88)	(294.88, 414.39)	(414.39, 592.3)	(8.84, 592.3)
B3	(1.79, 31.07)	(31.07, 42.68)	(42.68, 53.75)	(53.75, 65.86)	(65.86, 85.4)	(1.79, 85.4)
C1	(19.95, 43.57)	(17.75, 19.95)	(15.24, 17.75)	(12.22, 15.24)	(7.3, 12.22)	(7.3, 43.57)
C2	(86.57, 129.18)	(70.24, 86.57)	(63.24, 70.24)	(54.13, 63.24)	(47.1, 54.13)	(47.1, 129.18)
C3	(0.66, 0.8)	(0.6, 0.66)	(0.54, 0.6)	(0.48, 0.54)	(0.38, 0.48)	(0.38, 0.8)
D1	(0.02, 15.06)	(−1.04, 0.02)	(−3.25, −1.04)	(−9.52, −3.25)	(−11.47, −9.52)	(−11.47, 15.06)
D2	(57.33, 94)	(50.21, 57.33)	(41.84, 50.21)	(30.99, 41.84)	(16.23, 30.99)	(16.23, 94)
D3	(53.15, 95.79)	(43.91, 53.15)	(32.21, 43.91)	(17.65, 32.21)	(4.59, 17.65)	(4.59, 95.79)
D4	(289.93, 670.96)	(233.76, 289.93)	(177.27, 233.76)	(112.86, 177.27)	(28.32, 112.86)	(28.32, 670.96)

). Specifically, in this study, we used all data from 2001 to 2020 to determine the warning thresholds for each indicator based on the median, the mean, and the maximization principle (two-thirds of the total), dividing the warning and no-warning zones (I). We used the natural breakpoint method for more detailed grading, delineating light warnings (II), moderate warnings (III), severe warnings (IV), and extreme warnings (V) for zones with warnings. We took the year 2001 as the baseline to establish a relative risk system, which was used to analyze the secure state of CLQ in the study area. This threshold determination method considers the influence of economic development, land reclamation, human activities, and ecological conditions on the CLQ in the whole region.

3.1.2. Indicator Levels and Warning Result Levels

As for the positive and negative indicators, the corresponding fuzzy side distance correlation functions were determined, and the dominance index of each indicator was calculated according to Equations (4)–(13). Additionally, each indicator’s early warning level was established according to the maximum affiliation theory, and the unit with the highest dominance index value was determined as the early warning level of the indicator. Then, the integrated dominance index was calculated according to Equation (14) to determine the integrated warning level. The results for the cities in Liaoning Province in 2010 are presented as an example (

Table 3. Early warning indicators and results for Liaoning Province in 2010.

	A1	A2	A3	A4	B1	B2	B3	C1	C2	C3	D1	D2	D3	D4	Integrated
Shenyang	V	V	III	I	III	V	IV	I	V	IV	I	I	II	I	I
Dalian	IV	IV	III	II	III	V	V	I	III	IV	III	I	I	III	III
Anshan	III	III	III	I	II	IV	IV	I	I	II	I	III	I	II	I
Fushun	IV	II	IV	I	I	II	III	I	I	I	II	IV	I	II	I
Benxi	IV	IV	IV	I	I	II	II	I	I	I	II	V	I	III	I
Dandong	III	II	III	I	I	II	IV	I	I	II	II	IV	I	III	I
Jinzhou	I	II	II	I	I	IV	IV	I	IV	IV	I	I	II	II	I
Yingkou	II	III	III	I	III	V	V	I	II	III	II	III	I	I	I
Fuxin	II	II	II	I	II	II	II	I	IV	IV	II	I	I	II	II
Liaoyang	II	II	IV	II	II	IV	V	I	II	II	I	II	II	I	II
Panjin	III	II	IV	II	IV	IV	III	I	V	IV	II	I	II	I	I
Tieling	I	II	III	I	I	III	V	I	II	III	I	I	I	I	I
Chaoyang	I	I	III	V	I	II	II	I	II	II	III	II	II	II	II
Huludao	I	II	II	I	II	III	III	I	II	II	III	II	II	III	II

) to illustrate the correlation between the indicator level and the integrated warning result level. Table 3 shows that most cities have a no-warning level, but each indicator has a distinct warning level. The integrated warning level in Shenyang is a no warning, but indicators A1, A2, B2, and C2 are extreme warnings. Indicators C1, D2, and D3 in Dalian are no warnings, but the integrated warning level is a moderate warning. Indicators A1, A2, A3, B1, C2, and D1 in Fuxin are light warnings, and the integrated warning level is also a light warning. As a result, each indicator’s warning level has a unique impact on the integrated warning results. The superiority of the performance of the matter–element model lies in its ability to transform the incompatible relationship between multiple indicators into a compatible relationship, showing the correlation between the individual indicator levels and the integrated warning results.

3.1.3. CLQ Early Warning Results Based on the Improved Matter–Element Model

An early warning indicator system for CLQ was constructed using 2001 as the base year to analyze changes in CLQ in the black soil region of northeast China from 2001 to 2020; Figure 2 shows the spatial variation of the CLQ early warning. There are now significantly fewer areas in which warnings are present, indicating that CLQ has changed over time to become warning-free. As shown in Figure 2a, the CLQ in the northern part of the black soil region is significantly better than that in the southern part; however, there are still CLQ warnings in nearly half of the cities. In 2005, the warning level clearly decreased in the eastern part of the black soil region, but significant warnings remained in the western part. As shown in Figure 2c,d, most cities have no warnings, and only a few cities have light or moderate warnings. There are no dramatic warnings for the CLQ in 2020, except for a few cities with light or moderate warnings.

3.2. Spatiotemporal Driving Factors of CLQ

A panel data model was constructed to analyze the driving factors of spatiotemporal variability affecting CLQ (Section 2.5). The early warning integrated dominance index is the explained variable, and the 14 indicators are the explanatory variables. We tested the raw data and the panel data model sequentially by a unit root test, co-integration test, and model test, and then we estimated the parameters. The estimated parameters reflect the correlation between each early warning indicator and the integrated dominance index and indicate whether each indicator has a positive or negative impact on the CLQ early warning results.

3.2.1. Unit Root Test

Stationarity tests of the 14 indicators were performed using Stata software. We performed unit root tests for the explained variable y (the integrated dominance index) and the explanatory variables (A1, A2, A3, A4, B1, B2, B3, C1, C2, C3, D1, D2, D3, and D4) using IPS and HT tests, respectively. As shown in

Table 4. The results of the unit root test.

Variable	IPS	HT	Result
y	−4.075 *** (0.000)	0.233 *** (0.000)	Stationary
A1	−1.031 (0.151)	0.604 (0.051)	Non-stationary
A1(D1)	−14.663 *** (0.000)	−0.040 *** (0.000)	Stationary
A2	−1.338 (0.091)	0.406 *** (0.000)	Non-stationary
A2(D1)	−22.218 *** (0.000)	−0.283 *** (0.000)	Stationary
A3	3.910 (1.000)	0.781 (1.000)	Non-stationary
A3(D1)	−12.056 *** (0.000)	0.135 *** (0.000)	Stationary
A4	−15.696 *** (0.000)	−0.149 *** (0.000)	Stationary
B1	−6.052 *** (0.000)	0.018 *** (0.000)	Stationary
B2	8.415 (1.000)	0.752 (0.997)	Non-stationary
B2(D1)	−7.901 *** (0.000)	0.139 *** (0.000)	Stationary
B3	3.594 (1.000)	0.617 (0.106)	Non-stationary
B3(D1)	−15.06 *** (0.000)	−0.182 *** (0.000)	Stationary
C1	−11.643 *** (0.000)	−0.118 *** (0.000)	Stationary
C2	2.404 (0.992)	0.834 (1.000)	Non-stationary
C2(D1)	−2.718 *** (0.000)	0.482 *** (0.000)	Stationary
C3	−11.216 *** (0.000)	−0.089 *** (0.000)	Stationary
D1	−0.667 (0.252)	0.429 *** (0.000)	Non-stationary
D1(D1)	−18.393 *** (0.000)	−0.184 *** (0.000)	Stationary
D2	3.300 (1.000)	0.730 (0.983)	Non-stationary
D2(D1)	−15.374 *** (0.000)	0.032 *** (0.000)	Stationary
D3	2.370 (0.991)	0.772 (1.000)	Non-stationary
D3(D1)	−10.638 *** (0.000)	0.301 *** (0.000)	Stationary
D4	−1.804 ** (0.036)	0.556 *** (0.001)	Stationary

Note: A1(D₁) is the first-order difference for indicator A1. ***, ** represent significant at statistical level of 1%, 5% and 10% respectively.

, variables y, A4, B1, C1, C3, and D4 are stationary. Although variables A1, A2, A3, B2, B3, C2, D1, D2, and D3 are non-stationary, their first-order difference results are stationary.

3.2.2. Cointegration Test

Because a portion of the raw data is non-stationary, which can lead to pseudo-regression errors in the regression results, it was necessary to perform cointegration tests to eliminate the negative effects of the data. Therefore, we took the raw data of all variables and the raw data of some variables (stationary state) with the first-order difference data of other variables (non-stationary state) as the input data to perform two co-integration tests. As shown in

Table 5. The result of the cointegration test.

Data Type	Statistical Values	T	P
Raw data	Modified Dickey–Fuller t	−4.5754	0
	Dickey–Fuller t	−8.0923	0
	Augmented Dickey–Fuller t	−3.158	0.008
Raw data and first-order differential data	Modified Dickey–Fuller t	−8.4458	0
	Dickey–Fuller t	−10.1816	0
	Augmented Dickey–Fuller t	−4.8681	0

, the raw and first-order difference data of the explanatory variables were all co-integrated with the explained variable in both experiments.

3.2.3. Panel Data Model Test

In this study, the F test, FM test, and Hauman test were used to test different types of panel models (

Table 6. The results of the panel data model test.

Data Type	Test Method	Statistical Values	P	Hypothesis		Test Result
				H1	H0
Raw data	F test	17.35	0	FE	ME	H1: FE
	LM test	869.58	0	RE	ME	H1: RE
	Hauman test	62.33	0	FE	RE	H1: FE
Raw data and first-order differential data	F test	17.84	0	FE	ME	H1: FE
	LM test	1000.94	0	RE	ME	H1: RE
	Hauman test	76.16	0	FE	RE	H1: FE

). The test results indicate that all tests are significant at the 1% level, leading to the rejection of the null hypothesis. In the two experiments, we obtained the consistent conclusion that the FE model outperforms the ME model, the RE model outperforms the ME model, and the FE model outperforms the RE model after the F test, FM test, and Hauman test. Therefore, the FE model was used to analyze the significant relationship between the explanatory and explained variables.

3.2.4. Regression Results of the Panel Data Models

Two sets of experiments were performed to confirm the impacts of data stationarity and co-integration on the outcomes of regression analysis. In one set of experiments (

Table 7. Regression results of the different panel data models.

Part 1				Part 2
Variable	Models			Variables	Models
	ME	RE	FE		ME	RE	FE
A1	−0.0023 **	−0.0017 ***	−0.0008	A1(D1)	0.0001	0.0001	0.0000
	(−2.6157)	(−3.2917)	(−1.3969)		(0.0456)	(0.1095)	(−0.0190)
A2	0.0005	−0.0020 ***	−0.0026 ***	A2(D1)	−0.0006	−0.0004	−0.0004
	(0.3752)	(−3.3386)	(−4.1581)		(−0.7673)	(−0.5393)	(−0.5367)
A3	−0.0011	−0.0010 **	−0.0009 *	A3(D1)	−0.0020	−0.0005	−0.0002
	(−1.1614)	(−2.0746)	(−1.7488)		(−1.6364)	(−0.4950)	(−0.1982)
A4	0.0020 ***	0.0017 ***	0.0017 ***	A4	0.0021 ***	0.0016 ***	0.0015 ***
	(5.3208)	(5.3520)	(5.5239)		(4.8386)	(4.8551)	(4.4034)
B1	−0.0190 **	−0.0123 ***	−0.0084 **	B1	−0.0189 **	−0.0086 *	−0.0050
	(−2.4192)	(−2.8537)	(−1.9654)		(−2.3687)	(−1.7511)	(−1.0120)
B2	−0.0002	−0.0001	−0.0020 ***	B2(D1)	−0.0059 ***	−0.0033 ***	−0.0030 ***
	(−0.9182)	(−0.9502)	(−4.0388)		(−3.0559)	(−3.1382)	(−2.9527)
B3	−0.0002	−0.0022 **	−0.0022 **	B3(D1)	0.0021	−0.0007	−0.0012
	(−0.1777)	(−2.4536)	(−2.0449)		(1.3997)	(−0.4067)	(−0.7329)
C1	0.0065 ***	0.0076 ***	0.0073 ***	C1	0.0021	0.0069 ***	0.0073 ***
	(3.0641)	(8.3536)	(8.0117)		(1.1083)	(7.0800)	(7.3627)
C2	0.0009	0.0022	0.0209 ***	C2(D1)	0.0016	−0.0299	−0.0325 *
	(0.3985)	(1.2477)	(3.2691)		(0.0662)	(−1.5886)	(−1.7454)
C3	0.5170 ***	0.4168 ***	0.4239 ***	C3	0.5382 ***	0.4573 ***	0.3960 ***
	(3.7052)	(4.9059)	(4.8351)		(4.9382)	(5.4798)	(4.3455)
D1	0.0105 ***	0.0108 ***	0.0113 ***	D1(D1)	0.0033	0.0032	0.0032
	(2.9915)	(4.8716)	(5.2041)		(1.0419)	(1.2150)	(1.2449)
D2	0.0032 ***	0.0049 ***	0.0052 ***	D2(D1)	0.0034 ***	0.0020 **	0.0015
	(3.3978)	(7.1506)	(7.1650)		(3.5421)	(2.0031)	(1.5319)
D3	0.0017	0.0030 *	0.0167 ***	D3(D1)	0.0442 ***	0.0135	0.0093
	(0.9046)	(1.8694)	(4.5051)		(3.4390)	(1.2861)	(0.8972)
D4	0.0004 **	0.0008 ***	0.0009 ***	D4	0.0007 ***	0.0010 ***	0.0011 ***
	(2.2789)	(9.3454)	(10.6688)		(5.8150)	(15.9082)	(15.7183)
R2	0.562		0.654	R2	0.471		0.582

Note: ***, ** and * represent significant at statistical level of 1%, 5% and 10% respectively.

, Part 1), we performed a regression analysis using raw data for all variables. In another set of experiments (Part 2), we used raw data for some variables (stationary state) and first-order difference variables for others (non-stationary state) as the input data. The FE model in the Part 2 experiment was determined to be the final panel data model suitable for the analysis of the spatiotemporal driving factors of CLQ. As shown in Table 7, because the unit root and cointegration tests were not performed in the Part 1 experiment, more explanatory variables were significantly correlated with the explained variables because more explanatory factors were misfit with the explained variables due to the pseudo-regression phenomena. The regression results of the FE model in the Part 2 experiment show that five explaining variables (A4, B2, C1, C3, and D4) were significantly correlated with the explained variables. One explanatory variable (B2) was significantly and negatively correlated with the CLQ warning results. Variable B2 (land carrying capacity) indicates the population density of each city. The more people there are in a specific land area, the more pressure there is on the environment, and the more the CLQ is affected. Four explanatory variables (A4, C1, C3, and D4) were significantly and positively correlated with the CLQ warning results. Variable A4 (agricultural GDP growth rate) indicates the state of agricultural economic development, which has a potential impact on improvements in the CLQ. Variable C1 (annual precipitation) primarily functions in the grain-growing period. A high value of variable C3 (habitat quality index) implies a healthy ecological environment, which benefits the CLQ. In addition to being essential for food security, variable D4 (grain yield per unit of cultivated land area) serves as a critical indicator of the CLQ.

4. Discussion

The black soil region of northeast China is a crucial grain-producing region, and its CLQ trend directly affects national food security [41]. Socioeconomic development in the last 20 years has contributed to improvements in the ecological environment of the arable land, which have a potential influence on the ecological safety of CLQ [42,43]. In 2001, the northern portion of the black soil region had better CLQ than the southern portion, with Chifeng, Daqing, Mudanjiang, and Yanbian under severe warnings. This phenomenon can be attributed to the fact that these cities are in a state of high urbanization, with a concentrated urban population, underutilization of arable land, and a decline in rural population. There is a correlation between the urbanization rate and the CLQ [44]. With the development of society and the implementation of policies related to cultivated land protection, arable land has been effectively protected and appropriately utilized in the black soil region, which plays an essential role in improving food production, protecting cultivated land security, and maintaining ecological health [45].

Table 8. Statistical analysis of CLQ warning results.

Early Warning Level	Security Level	Number	Percentage (%)
No warning	Security	544	69.74
Light warning	Proximity security	148	18.97
Moderate warning	Critical Insecurity	60	7.69
Severe warning	Insecurity	27	3.46
Extreme warning	Extremely insecurity	1	0.13

shows the statistical results for the different warning levels between 2001 and 2020. Among them, extreme warnings only account for 0.13%. The results of this study demonstrate that most cities in northeast China were in a secure state in 2020 (Table 8), indicating that the improvements in the CLQ and the ecological preservation of cultivated land both significantly benefited from socioeconomic development.

In this study, we divided CLQ changes into three types—stable, gradual, and abrupt—to analyze the trends in CLQ changes over the previous 20 years. We developed classification guidelines by statistically analyzing the CLQ in 39 cities, as shown in

Table 9. Principles of classification of the types of changes in CLQ.

Type of Change	Number	Principle
Abrupt change	6	$\left|{\alpha }_i-{\alpha }_{i-1}\right|\ge 2\mathrm{and}f\left(i\right)\ge 2$
Gradual change	6	$(\frac{d{\alpha }_i}{di}>0\mathrm{or}\frac{d{\alpha }_i}{di}<0$$)\mathrm{and}f\left(i\right)\ge 2$
Stable change	27	Others

Note: ${\alpha }_i$ denotes the CLQ warning level in year $i$, and $f\left(i\right)$ denotes the number of year $i$ appearances; order of importance: abrupt change > gradual change > stable change.

. An abrupt change indicates that a city’s CLQ early warning level experienced a dramatic increase or decrease in the past 20 years. This type of change requires the attention of government departments, including an analysis of the current status of CLQ in the city and prompt identifications of possible safety hazards. A gradual change indicates a continuous increase or decrease in a city’s CLQ early warning level over time. It is essential to focus on the ongoing decline and the factors that may contribute to the degradation of CLQ. A stable change indicates that a city’s CLQ has been in a stationary state; cities of this type are in an almost warning-free state. This classification method has an essential guiding role in CLQ monitoring and early warning, providing scientific decision-making support for black soil conservation and utilization by indicating the degree to which attention should be paid to areas with varying change tendencies.

The results of the panel data model show (Table 7) that the primary driving factors of CLQ are agricultural GDP growth rate, land carrying capacity, grain yield per unit of cultivated land area, annual precipitation, and habitat quality index. These results indicate that the index system established in this study can scientifically provide early warning results for CLQ, taking into account four dimensions: socioeconomic development, human–land relationships, ecological environment, and cultivated land characteristics [46,47]. The development of appropriate political, technical, and management strategies can mobilize farmers’ willingness to produce to attain higher incomes. However, these strategies must be based on the appropriate exploitation of cultivated land [48]. A high-quality tillage layer and management strategy can contribute significantly to improvements in CLQ [49,50]. Continued population growth and human activity put varying degrees of external pressure on cultivated land, increasing the demand for multiple agricultural products [51,52,53]. The overutilization of cropland can also damage its health and lead to a decline in its quality [54,55]. An increase in grain yield indicates an improvement in the cropland’s productivity, which is also influenced by other factors, such as cultivated land conditions, cultivated land fragmentation, the application of chemical fertilizers, and conservation policies [56,57,58]. As a result, an increase in grain yield does not entirely reflect an improvement in the CLQ, but it is undeniable that grain yield is an essential sign of the CLQ state [59]. The rainy season in northeast China is concentrated in the summer; these conditions of high temperatures, rainfall, and long periods of sunshine are favorable for crop growth [60,61]. Therefore, rainfall showed a significant positive effect on the CLQ. The habitat quality index considers the influence of nearby land-use types on arable land, serving as a characterization indicator of habitat quality. In addition to reflecting regional biodiversity, it also illustrates the impact of various habitat threats and disturbances on CLQ.

5. Conclusions

Timely monitoring and accurate early CLQ warnings are essential for cultivated land safety. In this study, we constructed a DPSR model to select indicators covering four dimensions: socioeconomic drivers, human activity pressure, ecological state, and cultivated land response. We introduced fuzzy sets into the matter–element model to construct a fuzzy side distance-based matter–element model for CLQ assessment and early warning. The experimental results in the black soil region of northeast China from 2001 to 2020 validate the applicability of the optimized matter–element model for ecological safety assessments and CLQ early warning. In the past two decades, socioeconomic development has effectively improved the ecological safety of CLQ. Compared to 2001, Dalian, Yingkou, Fuxin, Panjin, Chaoyang, and Huludao had an early warning situation for CLQ in 2020. Specifically, Yingkou, Fuxin, Panjin, Chaoyang, and Huludao are under light warnings, and Dalian is under a moderate warning. In the whole region, the CLQ is the best in Heilongjiang Province and the worst in Liaoning Province. This result implies that rapid socioeconomic development positively affects CLQ. To better understand the key variables influencing CLQ, a panel data model was constructed to analyze the critical spatiotemporal drivers of CLQ. The experimental results demonstrate that the developments of the agricultural economy, human activities, grain yield, precipitation, and habitat quality were significant driving factors in changes in CLQ. The analysis results reported in this study provide a scientific basis for regional CLQ assessment and early warning.

However, due to the limitations of experimental conditions and difficulties in obtaining data, we did not use natural soil conditions and profile data in this study to assess the ecological safety of cropland quality. Time-series soil profile data, especially for soil natural conditions, physical properties, and chemical properties, are difficult to obtain. This is a challenging task, because such data require annual monitoring and specialized equipment. In region-scale areas, it is even more difficult to obtain such data accurately.