## *Article* **A Two-Stage Fuzzy Optimization Model for Urban Land Use: A Case Study of Chongzhou City**

**Jinjiang Yao 1 , Bingkui Qiu 2 , Min Zhou 3 , Aiping Deng 3,\* and Siqi Li 3**

**Abstract:** Under the background of New-type Urbanization, with the continuous advancement of urbanization and the all-round development of cities, all kinds of demands are also rising. In the case of demand, it is difficult to quickly adjust from the land supply side and to guide the optimization of the structure and layout of land use is one of the methods to achieve this based on the current situation and shortage of urban land use structure and spatial arrangement. Because of the complexity, uncertainty and dynamics of the land use system, it is necessary to use an uncertain model to accurately describe and propose the approximate optimal solution, so this study analyzes the influencing mechanism of land use and optimize the land use structure under uncertainties by using a Bayesian network and fuzzy mathematical programming. Based on the results of the two stages of analysis, the cellular automata simulation is completed. The framework is applied to Chongzhou city in western China. The results indicated that the optimal land space for cultivated land is in the middle and the south based on the joint influence probability of arable land and urban construction land. The conversion probability of the area near the east is low, and the joint impact probability of construction land in all areas is generally similar except for the western protection area. After the optimization of the fuzzy planning, the optimal construction land scale is 69.42 km2 . Under the condition that the cultivated land's red line is guaranteed, there is still 98.87 km2 of space for the increase in cultivated land. It is found through simulation that the increase in construction land would occur in the central and western parts of Chongzhou, which may be caused by the urban siphon effect. According to Monte Carlo verification, when the conversion probability exceeds 50%, the cultivated land could be turned into urban construction land, with an accuracy of 91.99%. Therefore, this proposed framework is helpful to understand the process of land use and provides a reference for making scientific and reasonable territorial spatial planning and guiding land use practice under uncertainties.

**Keywords:** two-stage land use; Bayesian network; fuzzy mathematical programming; Chongzhou city; uncertainty

## **1. Introduction**

Land-use systems are complex, diverse and large, and the academic community tries to narrow the boundaries of land-use systems from the perspective of scenario setting to facilitate the use of various models and methods for research [1]. In terms of land use simulation, different land use systems, methods and models derive different driving mechanisms [2–4] containing various assumptions in the simulation process, which provide new ideas and angles for exploring the uncertainty in land use system and form our understanding of land use system today. Although there are systematic differences in scenario-based land use simulation and modeling, uncertainty does not change with the model methods and scenarios set [5], namely, changes caused by changing model

**Citation:** Yao, J.; Qiu, B.; Zhou, M.; Deng, A.; Li, S. A Two-Stage Fuzzy Optimization Model for Urban Land Use: A Case Study of Chongzhou City. *Sustainability* **2021**, *13*, 13961. https://doi.org/10.3390/su132413961

Academic Editors: Marc A. Rosen, Miguel Amado and George D. Bathrellos

Received: 16 October 2021 Accepted: 15 December 2021 Published: 17 December 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

parameters and modifying a hypothesis do not change the uncertainty as well. Reducing uncertainty in land use forecasts is desirable, and some technical efforts can be made, but it is problematic to determine an "accurate" model for a specific purpose. Compared with other systems, the complexity and hugeness of land use systems do not allow repeated proofreading and collection of data that occurred in the past, which makes it more difficult to test the validity of simulation. Reducing such known uncertainties is a common problem faced by experts and scholars. What academics need to pay attention to is how to improve their understanding of the laws of land use systems and how to use these uncertainties to guide land use. It is a current problem that needs to be solved urgently, and it is also the way used to guide practice in a timely manner. Moreover, with the continuous advancement of urbanization and the all-round development of cities, all kinds of demands are also rising, so it is not feasible to reduce the demand in a short period of time [6]. Therefore, in the case of demand, it is difficult to quickly adjust from the land supply side, and to guide and play the role of land in the supply side reform elements may be one of the solutions.

With the development of land use research, more and more scholars have been inspired to capture uncertainty. The uncertainty of spatial factors will affect land use decision making, and incorporating uncertainty into the optimization of land use structure has gradually become one of the important modeling directions of the future [7,8]. Land use structure optimization and spatial optimization together provide an important theoretical basis for the formulation of territorial spatial planning. Territorial spatial planning is an important guide for land use practice, and planning often has problems with unclear planning intentions and inaccurate spatial grasp at the driving level of land use changes [9]. This is due to the inappropriate perception by the framers of uncertainties and laws in land use. In addition to emphasizing the situational and social construction in space [10], planning should also objectively describe and measure the actual existence and summarize it, including correctly understanding the law of land use and identifying various uncertainties [5]. Uncertainty is inherent. Although the tools used to study land use change may never know all the information, some uncertain methods are still expected to be applied to express and capture some of the uncertainties. Some representative studies are as follows, Verstegen et al. (2012) constructed a PCraster land use change model, which integrated the functions of visualization, simulation and uncertainty analysis. Van et al. (2012) simulated land supply with known drivers of land use, and then analyzed uncertainty using the Monte Carlo method. Meyer et al. (2012) used Bayesian network to help assess the uncertainties brought about by stakeholder knowledge modeling. Different from the uncertainty of spatial influencing factors, the uncertainty in the optimization of land use structure is often ignored. Most studies used ordinary planning to balance the land use demand of various sectors [11–15]. However, it is necessary to consider the uncertainty in the optimization of the land use structure. In order to perfect the research on this aspect, fuzzy mathematical programming is applied to it. Moreover, the uncertainty in the land use system has not been widely studied, and there are few relevant data. Some studies that take into account this problem are also always described as random distributions [16], which means that some research choses to avoid uncertainty and generally define the sources of uncertainty. The related research only analyzed the uncertainty of the influencing factors, but seldom applied the results obtained by the uncertainty method to the spatial simulation.

Mathematical models are widely used in optimization research in various fields to solve the problems of multi-constraint conditions and multi-objective optimal configuration. In the field of land use, mathematical planning has long been applied. The mathematical programming model is the optimal allocation scheme determined by mathematical programming theory [17–27]. For example, integer programming [28,29] is used to solve the optimization problem of integer variables. Dynamic programming [30] is used to solve the optimal problem of decision-making processes. Linear programming [31,32] is used to solve the one-dimensional optimal allocation problem. Nonlinear programming [33–35] is used to solve multi-dimensional allocation optimization problems. Among them, the use of integer programming needs to limit the variable to integers, so the scope of application is not wide enough. Dynamic programming can realize the optimization of the dynamic decision-making process of linkage, but there may be a multi-dimensional solution dilemma. Generally, the time cost is high. Linear programming is simple to solve and more applicable than other kinds of programming, but only for simple linear programming problems. Fuzzy programming, random programming, chance-constrained programming, parameter programming, etc., are all linear programming methods. Stochastic programming is mainly used for stochastic processes, which can be solved by intelligent models. Parameter planning is applicable to exploring the influence of parameter changes on the final results and has a specific scope of application. Both types of planning are scenario specific. Fuzzy programming and chance-constrained programming have a relatively wide range of applications. They can express the general decision-making process and can be solved even when constraints are violated and the value of variables is not clear. The results can be used as a reference for decision makers. The applicability of nonlinear programming is not as broad as that of linear programming because of its multi-dimensional characteristics, but nonlinear programming has better performance for some complex problems.

In view of the shortcomings of the above models, the purpose of this paper is to fully consider the uncertain factors in the land use system and optimize the land use quantity and the spatial layout by combining cellular automata model so as to provide quantitative support for decision makers. Based on the complexity, uncertainty and dynamic characteristics of land use systems, the modeling method of capturing uncertainty is used to find more spatial information from the perspective of uncertainty, and the function of multiple spatial impact factors was constructed to depict the law of land use impact. To a certain extent, this expands the depth of the analysis of influencing factors, provides certain theoretical significance for the study of land use development and further deepens the cognition of the uncertainty in the process of land use. Under the background of the new urbanization strategy proposed by the Chinese government to pay more attention to the coordinated development of population, economy, society and ecological environment [36], based on the current situation and shortage of urban land use structure and spatial arrangement, combined with regional development policies, starting from the reform of land supply side structure, using the method of Bayesian network and fuzzy mathematical programming, this paper makes an in-depth analysis of the influencing mechanism of land use and the optimization of land use structure in two stages. It is helpful to provide theoretical support for the formulation of territorial space planning. The optimization result aiming at economic carbon efficiency is combined as well to provide a certain reference for the sustainable development of urban land use. China's New-type Urbanization and agriculture-oriented cities will be used as case demonstration.

#### **2. Methodology**

## *2.1. Two-Stage Land Use Optimization Model*

Figure 1 shows a two-stage land use optimization framework under uncertain conditions. This study divides land use analysis into two stages. The first stage is to analyze land use influencing factors through a Bayesian network. It performs well in dealing with uncertainties and can reflect a more real world. The second stage is to optimize land use through fuzzy mathematical programming. In addition to economic and social constraints, this stage also involves environmental constraints marked by carbon emissions. In order to link the two stages, cellular automata are used to finally realize the analysis results of the two stages, and the neighborhood influence is added.

**Figure 1.** Two-stage optimization framework. **Figure 1.** Two-stage optimization framework.

#### 2.1.1. Bayesian Network for Land Use 2.1.1. Bayesian Network for Land Use

For the occurrence or non-occurrence of land use events, there are usually only two probabilities of 1 and 0, that is, either occurrence or non-occurrence. Few studies have considered the occurrence probability of land use events, and Bayesian network has the ability to solve this problem. Bayesian network is particularly flexible, and new information can often be found during modeling to improve the performance and applicability of the network [37]. It is a way to deal with uncertainty [12] and provides more information for interpretation of the evaluation results. Bayesian network clearly shows the complex relationship between land use change, its driving factors and its change process with its good graphical description method. Although the influence of each uncertain fac-For the occurrence or non-occurrence of land use events, there are usually only two probabilities of 1 and 0, that is, either occurrence or non-occurrence. Few studies have considered the occurrence probability of land use events, and Bayesian network has the ability to solve this problem. Bayesian network is particularly flexible, and new information can often be found during modeling to improve the performance and applicability of the network [37]. It is a way to deal with uncertainty [12] and provides more information for interpretation of the evaluation results. Bayesian network clearly shows the complex relationship between land use change, its driving factors and its change process with its good graphical description method. Although the influence of each uncertain factor on the estimation results can be measured by sensitivity analysis, it requires the estimation

of complex partial derivatives. Bayesian network not only does not need complex differentiation, but it is closer to the reduction of the real situation, in which the actual process is probabilistic in nature [38]. A node *X* points to another node *Y*, where *X* is the parent of *Y*. All parent nodes can be used as a set *X* = {*x*1, *x*2, . . . *xk*}, child nodes can also be represented by the set *Y* = {*y*1, *y*2, . . . *yk*}, represented as a directed acyclic graph, which is the advanced form of Markov chain. The analysis of land use suitability can also be completed by Bayesian network. The Bayesian formula is shown in Formula (1):

$$P(X|Y) = \frac{P(Y|X) \cdot P(X)}{P(Y)}\tag{1}$$

In Formula (1), *P*(*X*|*Y*) represents the probability that the parent node *X* is in a certain state under the condition of child node *Y*, that is, the posterior probability. *P*(*Y*|*X*) is the same. *P*(*X*) is the prior probability of *X* without regard to *Y*, and *P*(*Y*) is the same thing. According to Bayesian theory, in land use analysis, it is usually assumed that all factors are independent of each other, so their relationship can be expressed by the joint probability Formula (2):

$$P(\mathbf{x}\_1, \mathbf{x}\_2, \dots, \mathbf{x}\_k) = P(\mathbf{x}\_k | \mathbf{x}\_1, \mathbf{x}\_2, \dots, \mathbf{x}\_{k-1}) \cdots P(\mathbf{x}\_2 | \mathbf{x}\_1) P(\mathbf{x}\_1) \tag{2}$$

The construction of a Bayesian network involves qualitative and quantitative steps. The first step is to formulate the node and its network structure, and the second is to calculate the parameters of each node. It is found that policy and socio-economic factors show complex and close links under the influence of natural factors such as resource endowment and spatial geographical characteristics [4]. At the same time, the design and implementation of land planning will be affected by local socio-economic factors and external processes [6]. This study extracts candidate variables from spatial data, constructs the network structure according to common sense, expert knowledge and statistical correlation between variables and avoids creating parameter sets with excessive data as much as possible [39]. Therefore, in terms of spatial factor analysis, historical state, elevation, slope, accessibility, natural endowment and other factors are taken into consideration in the construction of the Bayesian network, as shown in Figures 2 and 3. First is the linear features of the impact of land use, railway effects on the city of the population, material and information exchanges with the outside world, belonging to the clear on the larger influence on urban development tool, so it leads to the improvement of the region and all aspects of development, both in cultivated land activities and urban economic activities, and can produce a good liquidity. In addition to assuming the connection with the city, the highway also serves as a bridge for communication within the city, facilitating the effective connection and interaction of internal elements. So linear features can lead to strip development around them. Secondly, point features have a certain radiating effect on urban construction, such as scenic spots and commercial centers. These point features are active centers of population and capital in cities, which attract population and investment, are easier to form intensive land use and have more potential to lead development. Moreover, natural endowments are vital to both urban and agricultural development. For agriculture, water supply will directly have an important impact on crop quality and yield. Rice, the main food crop in Chongzhou City, will have its yield affected by the lack of water resources. At the same time, water is an important resource guarantee for urban construction. Sufficient water supply will support urban life and production and, vice versa, will limit urban sustainable development. In order to improve the quality of water resources, provide a good resource endowment for urban life and production and improve the ecological and landscape conditions of the city, Chongzhou City has implemented measures such as the 'water control ten' river (lake) chief system, black and odorous water treatment and comprehensive water environment management. Therefore, this shows that the current agricultural development and urban construction are both limited by water and dependent on water.

**Figure 2.** Bayesian networks for cultivated land. **Figure 2.** Bayesian networks for cultivated land. **Figure 2.** Bayesian networks for cultivated land.

**Figure 3.** Bayesian network of urban construction land. **Figure 3.** Bayesian network of urban construction land. **Figure 3.** Bayesian network of urban construction land.

Bayesian network can help to search various laws in land use systems, and Gaussian radial basis functions can help to express these laws. The various influencing factors of land use change have different forces, so the attraction of each factor to different land types is also different. For example, the closer the area is to the city center, the easier it is to develop into urban construction land. The attraction of the city center is higher than that of the county center and the town center. The land close to the highway is easier to form strips of urban construction land on both sides of the highway, and the higher the grade of the highway, the easier it is. Similarly, rivers have a huge impact on urban development, and urban construction land is often gathered along the banks. Therefore, there are differences in the influence ranges of different influencing factors. It is generally believed that there is a certain rule with distance. In order to express the difference of forces Bayesian network can help to search various laws in land use systems, and Gaussian radial basis functions can help to express these laws. The various influencing factors of land use change have different forces, so the attraction of each factor to different land types is also different. For example, the closer the area is to the city center, the easier it is to develop into urban construction land. The attraction of the city center is higher than that of the county center and the town center. The land close to the highway is easier to form strips of urban construction land on both sides of the highway, and the higher the grade of the highway, the easier it is. Similarly, rivers have a huge impact on urban development, and urban construction land is often gathered along the banks. Therefore, there are differences in the influence ranges of different influencing factors. It is generally believed that there is a certain rule with distance. In order to express the difference of forces Bayesian network can help to search various laws in land use systems, and Gaussian radial basis functions can help to express these laws. The various influencing factors of land use change have different forces, so the attraction of each factor to different land types is also different. For example, the closer the area is to the city center, the easier it is to develop into urban construction land. The attraction of the city center is higher than that of the county center and the town center. The land close to the highway is easier to form strips of urban construction land on both sides of the highway, and the higher the grade of the highway, the easier it is. Similarly, rivers have a huge impact on urban development, and urban construction land is often gathered along the banks. Therefore, there are differences in the influence ranges of different influencing factors. It is generally believed that there is a certain rule with distance. In order to express the difference of forces generated by distance, a Gaussian radial basis function is adopted, as shown in Formula (3):

$$P\_k(\mathbf{x}\_k|m) = \exp\left(-\frac{1}{2\sigma\_k^2}(\mathbf{x}\_k - \overline{\mathbf{x}}\_k)^2\right) \tag{3}$$

$$\text{where } \mathbf{P}(\mathbf{x}\_k|m) = \mathbf{x}\_k \text{ is the conditional coefficient of the random variables } \mathbf{x}\_k \text{ and } \mathbf{x}\_k \text{ is the marginal distribution.}$$

2 ଶ ( − ̅ )ଶቇ (3) where () represents the probability of developing into a certain land type, and represents the standard deviation of the influencing factor variable. The larger the value of , the wider the influence range of the variable, while the smaller the value, the more 2 where () represents the probability of developing into a certain land type, and represents the standard deviation of the influencing factor variable. The larger the value of , the wider the influence range of the variable, while the smaller the value, the more where *Pk*(*xk* |*m*) represents the probability of developing into a certain land type, and *σ* represents the standard deviation of the influencing factor variable. The larger the value of *σ*, the wider the influence range of the variable, while the smaller the value, the more concentrated the influence range. *xk* represents the value of the variable, and *xk* represents the central value, that is, the variable value most likely to develop into the target land type, at which time the probability of *Pk*(*xk* |*m*) reaches the maximum.

(3):

(3):

#### 2.1.2. Land Use Optimization Mathematical Programming Model

Next, we estimated the uncertainty in the land use system is crucial for investigating land use policies, evaluated urban risk response strategies and quantified the impact of land use change on climate and the environment. Uncertainty involves a wide range, and this study uses a small entry point to express one aspect of the uncertainty, namely, the uncertainty of the target conditions for economic carbon efficiency, as shown in Formula (4):

$$\text{Max } Z = \frac{\sum\_{i=1}^{n} \mathbf{x}\_{i} \cdot \text{UB}\_{i}}{\sum\_{i=1}^{n} \mathbf{x}\_{i} \cdot \text{WC}\_{i} + \sum\_{i=1}^{n} \mathbf{x}\_{i} \cdot \text{PC}\_{i} + \sum\_{i=1}^{n} \mathbf{x}\_{i} \cdot \text{EC}\_{i} + \sum\_{i=1}^{n} \mathbf{x}\_{i} \cdot \text{SC}\_{i} + A \mathbf{C}}} \tag{4}$$

In Formula (4), *xi* represents the area of each type of land; *Max* means the maximum economic benefit; and *UBi* represents the level of economic benefit generated per unit land area (kg/km2 ). In this study, the gross output of primary industry, secondary industry and tertiary industry represents the gross output of cultivated land and urban construction land. *UBi* is the ratio of the gross output of each type of land to the area of each type of land before optimization. *Min* represents the minimization of carbon emission. *WCi* represents the carbon emissions (kg/km2 ) generated by wastewater quantified on all kinds of land, which can be obtained by the ratio of the amount of cultivated land wastewater and urban wastewater to the area of all kinds of land obtained from statistical data. *PCi* is the carbon emission from human body on various types of land (kg/km2 ). *ECi* represents the carbon emission from energy consumption of various types of land (kg/km2 ). *SCi* represents the carbon emissions from solid waste generated from various types of land (kg/km2 ). *AC* stands for carbon emissions from livestock and poultry farming. In order to more clearly understand the relationship between economic efficiency and carbon emissions, this research integrates economic benefit targets and carbon emission targets to explore the economic carbon efficiency, that is, how much economic benefit each unit of carbon emission produces and whether it is achieved between the sacrifice of the climate and the economic benefit balance. Therefore, the objective function is optimized.

#### (1) Environmental constraints

Since both agricultural and urban development have an impact on the soil and water environment, it is necessary to consider the status of soil. In the Formula (5), *r* represents soil and water loss rate (%). *Ar* represents the area of loss (km2 ). Both farmland operation and urban production and living will produce sewage. If sewage is discharged into the environment without treatment, it will not only pollute the water environment but also cause excessive heavy metals and poisons in the soil, which will increase the burden of the ecological environment and make it unbearable, affecting the environmental sustainability. In the Formula (5), *Gi* represents the sewage discharge coefficient (ton/km2 ) generated on a certain type of land, which can be obtained by the ratio of sewage discharge to the area of each type of land before optimization. *PG* denotes sewage treatment plant capacity (tons). Both agricultural production and urban production and life need to consume a large number of solid materials and produce a large amount of solid waste. Considering that the rural solid waste is treated on-site, which is not easy to count, only the solid waste generated in the city is considered. However, the processing capacity of urban solid waste is limited, so certain constraints must be met. In Formula (5), *USWi* represents the discharge coefficient of solid waste generated on a certain type of land (ton/km2 ), which can be obtained by the ratio of the amount of solid waste generated by the city to the area of urban construction land before optimization. *PSW* stands for solid waste disposal capacity (tons).

$$\begin{array}{ll}\sum\_{\substack{i=1\\n\end{l}}^{n}r \cdot \mathbf{x}\_{i} \lesssim A\_{r} \\ \sum\_{\substack{i=1\\n\end{l}}^{n} \mathbf{G}\_{i} \cdot \mathbf{x}\_{i} \lesssim P\_{\mathbf{G}} \\ \sum\_{i=1}^{n} \mathbf{G}\_{i} \cdot \mathbf{x}\_{i} \lesssim P\_{\mathbf{G}}\end{array} \tag{5}$$

#### (2) Economic constraints

The fixed asset investment required for agricultural production includes mechanical farming, electromechanical irrigation, mechanical seeding, mechanical harvesting, transportation roads and production facilities, etc. For Sichuan Province, with many mountains and hills, it belongs to an area with a high level of mechanization, so the investment ratio of fixed assets is also large. The fixed asset investment in urban construction has been the main part of the total investment. In Formula (6), *DSFAI* represents the total investment (CNY 100 million); *SFAIi* represents the investment per unit land area (CNY 100 million). In this study, the fixed asset investment in primary industry, secondary industry and tertiary industry represents the fixed asset investment in cultivated land and urban construction land respectively. *ISFAI* represents import investment (CNY 100 million). *OSFAI* represents export investment (CNY 100 million). According to Formula (6), the per capita annual ration consumption is about 135 kg:

$$\sum\_{i=1}^{n} (\text{SFAI}\_{i} \cdot \text{x}\_{i}) + \text{ISFAI} - \text{OSFAI} \gtrsim \text{DSFAI} \tag{6}$$

### (3) Social constraints

Agricultural products are one of its main production targets in the key county-level city for agriculture and ecological protection, so the balance of the food supply is particularly important. In Formula (7), *USi* represents the amount of grain produced per unit land area (kg/km2 ), which can be obtained by the ratio of grain yield to cultivated land area before optimization. *D* stands for grain demand (kg). Agricultural production and urban secondary and tertiary industries both need labor. In Formula (7), *ULi* represents the labor force per unit land area (person/km2 ). In this study, the employees in the primary industry, the secondary industry and the tertiary industry represent the labor force on cultivated land and urban construction land, and the ratio of them to each type of land area before optimization is taken as the value of *ULi* . *TL* stands for total labor force (people).

$$\begin{aligned} \sum\_{i=1}^{n} \mathbf{x}\_{i} \cdot \mathbf{U} \mathbf{S}\_{i} &\gtrsim D\\ \sum\_{i=1}^{n} (\mathbf{U} L\_{i} \cdot \mathbf{x}\_{i}) &\lesssim \mathbf{T} L \end{aligned} \tag{7}$$

#### (4) Programming constraint

The planning conditions mainly include the restriction of cultivated land protection red line and the restriction of construction land area. In Formula (8), *FC* represents the minimum cultivated land protected area (km2 ). *FB* refers to the maximum area of construction land (km2 ) A 338 km2 farmland protection red line has been set, accounting for 31.62% of the total area.

$$\begin{aligned} \alpha\_i &\ge FC, \ i = 1\\ \alpha\_i &\le FB, \ i = 2 \end{aligned} \tag{8}$$

(5) Other Constraints

$$x\_i \ge 0, \ i = 1, 2 \tag{9}$$

In Formula (9), *x*1 represents cultivated land. *x*2 represents urban construction land.

(6) Carbon effect estimation

Before estimating the carbon effect, it is necessary to define the organizational boundary. The organizational boundary set in this study is the carbon effect within the range of cultivated land production activities and urban consumption, which mainly includes carbon sources, carbon emissions generated by mechanical farming and urban production and living energy consumption in cultivated land activities, carbon emissions from

farmland and municipal wastewater treatment, carbon emissions from municipal solid waste treatment, carbon emissions produced by the human body and carbon emissions from livestock and poultry farming. In terms of carbon sinks, cultivated land, shrub land, woodland and water are all included. treatment, carbon emissions produced by the human body and carbon emissions from livestock and poultry farming. In terms of carbon sinks, cultivated land, shrub land, woodland and water are all included. 2.1.3. Cellular Automata Simulation

, = 2

In Formula (9), ଵ represents cultivated land. ଶ represents urban construction

Before estimating the carbon effect, it is necessary to define the organizational boundary. The organizational boundary set in this study is the carbon effect within the range of cultivated land production activities and urban consumption, which mainly includes carbon sources, carbon emissions generated by mechanical farming and urban production and living energy consumption in cultivated land activities, carbon emissions from farmland and municipal wastewater treatment, carbon emissions from municipal solid waste

#### 2.1.3. Cellular Automata Simulation Cellular automata determines cellular state through transformation rules [40]. The

(5) Other Constraints

(6) Carbon effect estimation

land.

*Sustainability* **2021**, *13*, x FOR PEER REVIEW 9 of 25

Cellular automata determines cellular state through transformation rules [40]. The improved transformation rules not only include the influence of cellular state and neighborhood action at the previous moment but also include quantified social, economic, spatial and other factors: improved transformation rules not only include the influence of cellular state and neighborhood action at the previous moment but also include quantified social, economic, spatial and other factors: ௧ାଵ = (௧, , … ) (10)

$$S\_{t+1} = f(S\_{t\prime}N\_{\prime}\dots) \tag{10}$$

≥ 0, = 1,2 (9)

In Formula (10), *S* refers to the set of all possible states, *f* is the transition rule, *N* is the neighborhood of the cell, *t* moment is the current moment and *t* + 1 represents the next moment. In fact, cellular automata take into account spatial explicit weighting factors in suitability, rather than completely random allocation. Each cell's land type is considered the primary type, and that cell is masked, while all other land types are considered the "claim" type and compete for that cell. The "claim" type of cell that influences the possibility of the future land type of the principal cell must be both inherently suitable and close to the principal cell. This study adopts a neighborhood range of 3 × 3, with a kernel as shown in Figure 4. In Formula (10), refers to the set of all possible states, is the transition rule, is the neighborhood of the cell, moment is the current moment and +1 represents the next moment. In fact, cellular automata take into account spatial explicit weighting factors in suitability, rather than completely random allocation. Each cell's land type is considered the primary type, and that cell is masked, while all other land types are considered the "claim" type and compete for that cell. The "claim" type of cell that influences the possibility of the future land type of the principal cell must be both inherently suitable and close to the principal cell. This study adopts a neighborhood range of 3 × 3, with a kernel as shown in Figure 4.

**Figure 4.** Cellular automata kernel.

**Figure 4.** Cellular automata kernel. This study defines three basic scenarios affecting land use change and the basic settings for each scenario. The first situation is the role of neighborhood. If there is an urbanized area around, the central cell is also inclined to develop into an urbanized cell under the influence of the neighborhood state, which is manifested as the concentration of the urbanized area. In the second case, the protected area is regarded as the forbidden area This study defines three basic scenarios affecting land use change and the basic settings for each scenario. The first situation is the role of neighborhood. If there is an urbanized area around, the central cell is also inclined to develop into an urbanized cell under the influence of the neighborhood state, which is manifested as the concentration of the urbanized area. In the second case, the protected area is regarded as the forbidden area for development, and it will not be developed into an urban area no matter whether it has the conditions for the development of urban area. Third case: In the case of limited city size, in order to ensure the balanced improvement of various functions of the city, it is not feasible to consider only a single goal. Therefore, the consideration of the red line of cultivated land protection and the control of construction land scale is added. This paper mainly studies two types of land with fierce competition, namely, cultivated land and urban construction land.

#### **3. Case Study**

Chongzhou is located in southwestern China, covering an area of 1089 km2 , located between 30◦300~30◦530 N and 103◦070~103◦490 E (Figure 5.). In terms of geographical conditions, the terrain of Chongzhou is high in the northwest and low in the southeast. The west is mostly high mountains, and the terrain is steep. High mountains account for 38.4% of the total area. It has mild climatic conditions, abundant agricultural conditions, more rainfall and less sunshine. In terms of history and culture, Chongzhou is one of the cradles of ancient Shu civilization and the birthplace of agricultural civilization. In terms of economy and society, the GDP of Chongzhou in 2017 was CNY 30.04 billion,

and the proportion of primary, secondary and tertiary industries was 12.1:49.5:38.4. The permanent resident population is 664,800, and the number of people employed in the primary, secondary and tertiary industries is 110,500, 282,800 and 149,600, respectively. Chongzhou is located in the west of Chengdu City, in the middle and upper reaches of Minjiang River, and is an important water source and ecological conservation area of Chengdu Plain. According to the strategy of "advancing to the east, expanding to the south, controlling to the west, transforming to the north, and providing excellent" in the General Plan of Chengdu City (2016–2035), Chongzhou is mainly responsible for the development of green industries and providing ecological services. Meanwhile, it strictly implements the ecological red line control, with cultivated land and woodland occupying the main position. According to existing data, generally speaking, cultivated land and urban construction land are mainly distributed in the central and eastern regions. Agriculture and urban construction will form a competitive relationship in terms of terrain. Thus, elevation and slope both constitute the main influencing factors of agricultural development and urban construction. In addition to the above linear features, namely, point features, elevation and slope conditions, for agricultural development, there are also natural conditions such as rainfall, soil productivity and heat that have a certain impact on the characteristics and productivity of agricultural products in the region. Taking into account the accessibility of data, this study selects rainfall as a representative to participate in the construction of the Bayesian network. In order to test the performance of the model, the Bayesian network modeling was replicated in Chongzhou ecological service zones, agricultural protection zones and major urban areas within Chengdu. On the one hand, natural endowment drives the stepwise development of land use. On the other hand, as one of the key pilot cities of New-type Urbanization in western China, it is still in the stage of urban expansion with highly developed social and economic activities, so it is suitable for this model. The general situation of Chongzhou is shown in Figure 6. *Sustainability* **2021**, *13*, x FOR PEER REVIEW 11 of 25

**Figure 5.** A location map of Chongzhou. **Figure 5.** A location map of Chongzhou.

**Figure 6.** Overview of the study area.

The data mainly came from Chengdu Yearbook 2018, Chongzhou Yearbook 2018 and Sichuan Statistical Yearbook 2018, which recorded the economic and social data in 2017. The elevation and spatial land use data are from the geographical science and resources institute of the Chinese Academy of Sciences and Tsinghua University; China's present situation of land use remote sensing monitoring database 30 m global land use data (http://data.ess.tsinghua.edu.cn/fromglc2017v1.html accessed on 6 March 2019) and spatial rainfall data were collected from National Tibetan Plateau Scientific Data Center [41,42]; and spatial data such as roads, water systems, attractions and commercial centers are derived from OpenStreetMap database. According to China's industrial structure and the status quo of land use, combined with the classification standard of China's land use status in 2017, this study divides land types into cultivated land, woodland, shrubby land, water area, protected area and urban construction land. The data on livestock and poultry rearing quantities, feeding cycles, electricity use, fuel use and management costs are from Chongzhou Yearbook 2018, Sichuan Statistical Yearbook 2018, Sichuan Survey Yearbook 2018, National Cost-Benefit Compilation of Agricultural Products 2018 and statistical data of Ministry of Agriculture and Rural Affairs, PRC (http://zdscxx.moa.gov.cn:8080/nyb/pc/search.jsp accessed on 16 December 2021). Nitrogen emissions, CO2, CH4 and N2O emissions factor from the 2006 national greenhouse gas inventory guide, the provincial greenhouse gas inventory preparation guide (2011) and

**Figure 6.** Overview of the study area.

**Figure 5.** A location map of Chongzhou.

**Figure 6.** Overview of the study area. The data mainly came from Chengdu Yearbook 2018, Chongzhou Yearbook 2018 and Sichuan Statistical Yearbook 2018, which recorded the economic and social data in 2017. The elevation and spatial land use data are from the geographical science and resources institute of the Chinese Academy of Sciences and Tsinghua University; China's present situation of land use remote sensing monitoring database 30 m global land use data (http://data.ess.tsinghua.edu.cn/fromglc2017v1.html accessed on 6 March 2019) and spatial rainfall data were collected from National Tibetan Plateau Scientific Data Center [41,42]; and spatial data such as roads, water systems, attractions and commercial centers are derived from OpenStreetMap database. According to China's industrial structure and the status quo of land use, combined with the classification standard of China's land use status in 2017, this study divides land types into cultivated land, woodland, shrubby land, water area, protected area and urban construction land. The data on livestock and poultry rearing quantities, feeding cycles, electricity use, fuel use and management costs are from Chongzhou Yearbook 2018, Sichuan Statistical Yearbook 2018, Sichuan Survey Yearbook 2018, National Cost-Benefit Compilation of Agricultural Products 2018 and statistical data of Ministry of Agriculture and Rural Affairs, PRC (http://zdscxx.moa.gov.cn:8080/nyb/pc/search.jsp accessed on 16 December 2021). Nitrogen emissions, CO2, CH4 and N2O emissions factor from the 2006 national greenhouse gas inventory guide, the provincial greenhouse gas inventory preparation guide (2011) and The data mainly came from Chengdu Yearbook 2018, Chongzhou Yearbook 2018 and Sichuan Statistical Yearbook 2018, which recorded the economic and social data in 2017. The elevation and spatial land use data are from the geographical science and resources institute of the Chinese Academy of Sciences and Tsinghua University; China's present situation of land use remote sensing monitoring database 30 m global land use data (http: //data.ess.tsinghua.edu.cn/fromglc2017v1.html accessed on 6 March 2019) and spatial rainfall data were collected from National Tibetan Plateau Scientific Data Center [41,42]; and spatial data such as roads, water systems, attractions and commercial centers are derived from OpenStreetMap database. According to China's industrial structure and the status quo of land use, combined with the classification standard of China's land use status in 2017, this study divides land types into cultivated land, woodland, shrubby land, water area, protected area and urban construction land. The data on livestock and poultry rearing quantities, feeding cycles, electricity use, fuel use and management costs are from Chongzhou Yearbook 2018, Sichuan Statistical Yearbook 2018, Sichuan Survey Yearbook 2018, National Cost-Benefit Compilation of Agricultural Products 2018 and statistical data of Ministry of Agriculture and Rural Affairs, PRC (http://zdscxx.moa.gov. cn:8080/nyb/pc/search.jsp accessed on 14 December 2021). Nitrogen emissions, CO2, CH4 and N2O emissions factor from the 2006 national greenhouse gas inventory guide, the provincial greenhouse gas inventory preparation guide (2011) and the United Nations Food and Agriculture Organization (http://www.fao.org/faostat/zh/#data/GM accessed on 16 July 2021). Data on the CO2-eq emission factors of standard coal are from China Energy Statistical Yearbook 2018, and the CO2-eq emission factor of electric power data are from literature. [43,44]. To avoid double counting, the organization boundaries of this study are beef cattle, pigs, mutton sheep and meat poultry, as well as the production area, feed area and management area of breeding in Chongzhou [45,46], as shown in Figure 7.

The classification of land use, the sign of land use conversion and the carbon emission form of the study are explained as follows: (1) Biomass storage and carbon sequestration of agricultural products are not considered in this study. (2) Because of the forest and herbs, there are many different crops and the waters of the biomass, and this paper summarizes land use cover change by the change of land type. The vegetation on the land's net ecosystem productivity (NEP) is also changed, and combined with the geographical position and climate characteristics, they determine the net ecosystem productivity of different land types. The average carbon sequestration capacity of different land types is also calculated. (3) The carbon emission form studied in this paper is CO2, that is, CO2 generated by various land use activities.

**Figure 7.** Organizational boundaries of this study. **Figure 7.** Organizational boundaries of this study.

#### **4. Results and Discussion**

[45,46], as shown in Figure 7.

#### The classification of land use, the sign of land use conversion and the carbon emission form of the study are explained as follows: (1) Biomass storage and carbon sequestration *4.1. Influencing Factors of Land Use: Bayesian Network*

of agricultural products are not considered in this study. (2) Because of the forest and herbs, there are many different crops and the waters of the biomass, and this paper summarizes land use cover change by the change of land type. The vegetation on the land's net ecosystem productivity (NEP) is also changed, and combined with the geographical position and climate characteristics, they determine the net ecosystem productivity of different land types. The average carbon sequestration capacity of different land types is also calculated. (3) The carbon emission form studied in this paper is CO2, that is, CO2 generated by various land use activities. **4. Results and Discussion**  *4.1. Influencing Factors of Land Use: Bayesian Network*  A Bayesian network is constructed and completed in Netica software (www.norsys.com accessed on 16 July 2021), using data from railways, highways, scenic spots, commercial centers, slopes, water area, rainfall, elevation and other data processed by ArcGIS. Considering that the area of the research area is the size of a county, this re-A Bayesian network is constructed and completed in Netica software (www.norsys.com, accessed on 16 July 2021), using data from railways, highways, scenic spots, commercial centers, slopes, water area, rainfall, elevation and other data processed by ArcGIS. Considering that the area of the research area is the size of a county, this research created 500 random sample points to train the Bayesian network with the sample data. Since the creation of random points results in the absence of values in the attributes of some points, 4 points were eliminated, and the total sample size became 496. The coordinate system of all kinds of spatial data was unified as GCS\_WGS\_1984, and all data processing was carried out in this coordinate system. The vector spatial data were converted to raster data through the element-to-raster operation, and the pixel size was 250 × 250 m. By generating the nearest neighbor table and extracting the point of arrival, all kinds of spatial attributes were linked with land use, and the attribute characteristics of each grid cell were obtained. Statistical spatial attribute data were the input type of the Netica software. Modules are constructed in Netica, and each module was divided into 10 continuous interval states. The initial probability of each module is equal so as to conduct initial exploration, determine parent nodes and child nodes according to the model design and form a directed acyclic graph.

the United Nations Food and Agriculture Organization (http://www.fao.org/faostat/zh/#data/GM accessed on 16 July 2021). Data on the CO2-eq emission factors of standard coal are from China Energy Statistical Yearbook 2018, and the CO2-eq emission factor of electric power data are from literature. [43,44].To avoid double counting, the organization boundaries of this study are beef cattle, pigs, mutton sheep and meat poultry, as well as the production area, feed area and management area of breeding in Chongzhou

search created 500 random sample points to train the Bayesian network with the sample data. Since the creation of random points results in the absence of values in the attributes of some points, 4 points were eliminated, and the total sample size became 496. The coordinate system of all kinds of spatial data was unified as GCS\_WGS\_1984, and all data processing was carried out in this coordinate system. The vector spatial data were converted to raster data through the element-to-raster operation, and the pixel size was 250 × 250 m. By generating the nearest neighbor table and extracting the point of arrival, all kinds of spatial attributes were linked with land use, and the attribute characteristics of each grid cell were obtained. Statistical spatial attribute data were the input type of the Figure 8 shows the probability map of influencing factors of cultivated land, and Figure 8 shows the probability map of influencing factors of urban construction land. DEM is the height, Slope is the slope, Accessibility is the accessibility to all kinds of lines and points, and so on. A state refers to protected areas, B state refers to cultivated land, C state refers to woodland, D state represents shrub land, E state represents water area and F state represents urban construction land. Among them, the data of accessibility, natural endowment and farming are only statistical data without spatial data, so the preliminary determination is based on the expert knowledge, so that the follow-up research can be carried out smoothly.

Netica software. Modules are constructed in Netica, and each module was divided into 10 continuous interval states. The initial probability of each module is equal so as to conduct

and form a directed acyclic graph.

out smoothly.

initial exploration, determine parent nodes and child nodes according to the model design

Figure 8 shows the probability map of influencing factors of cultivated land, and Figure 8 shows the probability map of influencing factors of urban construction land. DEM is the height, Slope is the slope, Accessibility is the accessibility to all kinds of lines and points, and so on. A state refers to protected areas, B state refers to cultivated land, C state refers to woodland, D state represents shrub land, E state represents water area and F state represents urban construction land. Among them, the data of accessibility, natural endowment and farming are only statistical data without spatial data, so the preliminary determination is based on the expert knowledge, so that the follow-up research can be carried

**Figure 8.** Probability diagram of factors affecting cultivated land: "DistanceToRailway" means the distance to the railway, "DistanceToRoad" means the distance to the road, "DistanceToScenic" means the distance to the scenic spot, "DistanceTo-Center" means the distance to the commercial center,"NaturalEndowment" refers to natural abilities or qualities,"Tractorploughing" refers to mechanized farming or not, "BasicStates" refers to initial land use status and "States" refers to simulated land use status. **Figure 8.** Probability diagram of factors affecting cultivated land: "DistanceToRailway" means the distance to the railway, "DistanceToRoad" means the distance to the road, "DistanceToScenic" means the distance to the scenic spot, "Distance-ToCenter" means the distance to the commercial center, "NaturalEndowment" refers to natural abilities or qualities, "Tractorploughing" refers to mechanized farming or not, "BasicStates" refers to initial land use status and "States" refers to simulated land use status.

According to the results of the Bayesian network in Figure 8, the linear ground features have the same characteristics for cultivated land and urban construction land, that is, the closer the land is to railways, highways, scenic spots and commercial centers, the easier it is to develop into cultivated land, and the trend of decline is from near to far. Due to the large amount of data, all the data of the 496 sample points were not included in this paper, so the sample data were explained in words. The first is the influence brought by the railway. As an important means of transportation to other places, the railway has a wider scope of influence than other linear ground objects. In Figure 8, the states of the "Distancetorailway" node are divided into 10, each state represents a continuous interval, and the value from A to J is increasing. The value of the sample can correspond to one of the intervals, so as to achieve the purpose of grading. The same applies to other state nodes. The distance to the railway has an interval length of 4740 m, as shown in Figure 8, and the state of intervals A and B, about 9 km within range, has the greatest attraction of arable land. With the increase in distance, the influence is more and more weak. The domestic railway's furthest distance influence is also not reduced to 0, which shows that the influence of the railway is wide, but the impact is only large in a certain neighborhood. The distance to the highway is 570 m, the distance to the scenic spot is 2690 m and the According to the results of the Bayesian network in Figure 8, the linear ground features have the same characteristics for cultivated land and urban construction land, that is, the closer the land is to railways, highways, scenic spots and commercial centers, the easier it is to develop into cultivated land, and the trend of decline is from near to far. Due to the large amount of data, all the data of the 496 sample points were not included in this paper, so the sample data were explained in words. The first is the influence brought by the railway. As an important means of transportation to other places, the railway has a wider scope ofinfluence than other linear ground objects. In Figure 8, the states of the "Distancetorailway" node are divided into 10, each state represents a continuous interval, and the value from A to J is increasing. The value of the sample can correspond to one of the intervals, so as to achieve the purpose of grading. The same applies to other state nodes. The distance to the railway has an interval length of 4740 m, as shown in Figure 8, and the state of intervals A and B, about 9 km within range, has the greatest attraction of arable land. With the increase in distance, the influence is more and more weak. The domestic railway's furthest distance influence is also not reduced to 0, which shows that the influence of the railway is wide, but the impact is only large in a certain neighborhood. The distance to the highway is 570 m, the distance to the scenic spot is 2690 m and the distance to the commercial center is 1710 m. The reason why the length of the interval is different is because of the limitation of computing speed and memory. As can be seen from the figure, the influence rules of the three are roughly the same. The optimal impact distance of highway is interval A, about 500 m; the best impact distance of the scenic spot is between A and B, which is about 5 km. The optimal impact distance of the commercial center is between A and B, which is about 3.5 km. In terms of natural endowment, in view of the accessibility of data, elevation, slope and rainfall are selected as state nodes. Among them, rainfall can be seen as a general trend, and regions with greater rainfall have better conditions for agricultural development. Compared with the data integrity of other modules, the anomaly in interval F in the figure may be due to insufficient rainfall data sampling. The probability of elevation's influence on land use state transition presents a "cliff" distribution, that is, more than 67.3% of the cultivated land is distributed below 800 m, which indicates that greater elevation will have

a great limitation on the current agricultural development in Chongzhou. The slope is the same. Moreover, 70.1% of the cultivated land is below 5◦ , and more than 5◦ will have great restrictions on mechanized farming. The slope is divided into five grades according to the relevant provisions of topographic slope classification technology in the second national land survey: I (0◦~2), II (2◦~6), III (6~15), IV (15~25) and V (>25), and above 25 degrees will forbid land reclamation. In fact, the cultivated land in Chongzhou is below the degree, so the slope has a great limitation on agricultural development, which is a factor that must be taken into account in the study. The influence of water sources is similar to that of roads, and it has a greater influence on cultivated land within the range of about 1.7 km between A and B and decreases with the increase in distance. In terms of land use evaluation, mechanized tillage does not play a decisive role in land state. Areas with good natural endowment are more likely to be converted to arable land.

As can be seen from Figure 9, 59.5% of the construction land elevation is mainly kept below 800 m, and 57.3% of the construction land slope is kept below 5◦ . The optimal impact distance of the water area is within the range of about 1.7 km between A and B. In the two state intervals of A and B, within about 9 km, the railway has the greatest attraction to urban construction land. The optimal impact distance of highway is interval A, about 500 m; the best impact distance of the scenic spot is between A and B, which is about 5 km. The optimal impact distance of the commercial center is between A and B, which is about 3.5 km. *Sustainability* **2021**, *13*, x FOR PEER REVIEW 15 of 25

**Figure 9.** Probability diagram of influencing factors of urban construction land: "DistanceToRailway" means the distance to the railway, "DistanceToRoad" means the distance to the road, "DistanceToScenic" means the distance to the scenic spot, "DistanceToCenter" means the distance to the commercial center, "NaturalEndowment" refers to natural abilities or qualities, "Tractorploughing" refers to mechanized farming or not, "BasicStates" refers to initial land use status and "States" refers to simulated land use status. **Figure 9.** Probability diagram of influencing factors of urban construction land: "DistanceToRailway" means the distance to the railway, "DistanceToRoad" means the distance to the road, "DistanceToScenic" means the distance to the scenicspot, "DistanceToCenter" means the distance to the commercial center, "NaturalEndowment" refers to natural abilities or qualities, "Tractorploughing" refers to mechanized farming or not, "BasicStates" refers to initial land use status and "States" refers to simulated land use status.

It is generally believed that all influencing factors, except elevation, slope and rainfall, should present the distribution as shown in Figure 10a if their influences on land use are equidistant attenuations. However, this is not the case in fact. In this study, the influence of various factors on land use under uncertain conditions was obtained through a

Secondly, the influence rules of various factors are not consistent. For example, the attenuation rate of railway influence is slow, while that of highway influence is faster. In order to express the influence law of each factor on land use, Gaussian radial basis function is used to express the non-inversely proportional and different influences. According to the Bayesian network constructed in Netica, the mean value and standard deviation of the influence probability of each factor can be obtained to determine the Gaussian radial basis function to express the influence probability of each factor. Figure 11 shows the Gaussian radial basis function curve of the influence factors of cultivated land drawn in MATLAB R2019a. Figure 12 shows the Gaussian radial basis function curve of influencing factors of urban construction land. Figures 11 and 12a from left to right in turn show: road (green), water (magenta), business center (yellow) and attractions (red) and railway (blue). It can be seen that under the same coordinate measurement, the influence of railway is one of the widest, the highway has the narrowest scope of influence but rail impacts few and far between. Although the road influence area is narrow, the road network is complex, and its influence cannot be underestimated. Second, the influence of the water range is relatively narrow, and the publicizing states of drainage systems are mainly the Golden Horse river, Minjiang river is flow and Black Rock stream (into Dujiangyan). The condition

*4.2. Comparative Analysis of Influence Factors*

#### *4.2. Comparative Analysis of Influence Factors*

It is generally believed that all influencing factors, except elevation, slope and rainfall, should present the distribution as shown in Figure 10a if their influences on land use are equidistant attenuations. However, this is not the case in fact. In this study, the influence of various factors on land use under uncertain conditions was obtained through a Bayesian network, and different results were obtained, as shown in Figure 10b. Firstly, the influence of various factors on land use is not inversely proportional to the distance. Secondly, the influence rules of various factors are not consistent. For example, the attenuation rate of railway influence is slow, while that of highway influence is faster. In order to express the influence law of each factor on land use, Gaussian radial basis function is used to express the non-inversely proportional and different influences. According to the Bayesian network constructed in Netica, the mean value and standard deviation of the influence probability of each factor can be obtained to determine the Gaussian radial basis function to express the influence probability of each factor. Figure 11 shows the Gaussian radial basis function curve of the influence factors of cultivated land drawn in MATLAB R2019a. Figure 12 shows the Gaussian radial basis function curve of influencing factors of urban construction land. Figures 11 and 12a from left to right in turn show: road (green), water (magenta), business center (yellow) and attractions (red) and railway (blue). It can be seen that under the same coordinate measurement, the influence of railway is one of the widest, the highway has the narrowest scope of influence but rail impacts few and far between. Although the road influence area is narrow, the road network is complex, and its influence cannot be underestimated. Second, the influence of the water range is relatively narrow, and the publicizing states of drainage systems are mainly the Golden Horse river, Minjiang river is flow and Black Rock stream (into Dujiangyan). The condition (within the territory of the longest) and the size of the tributary of about 180. Figure 11b shows rainfall (magenta, in millimeters), Figure 11c shows slope (red, in degrees) and Figure 11d shows elevation (blue, in meters). Slope (red, in degrees) is shown in Figure 12b and elevation (blue, in meters) is shown in Figure 12c. The analysis is described in Section 4.1.

#### *4.3. Affect Probability Distribution*

According to the characteristics of the Bayesian network, the probabilities of each node are independent from each other, so the joint influence probability of each factor can be calculated, which is used to express the relationship between each factor and is conducive to the transformation of independent influence probability to comprehensive influence probability. Figure 13 shows the spatial distribution of the joint influence probability of cultivated land. Figure 14 shows the spatial distribution of probability of joint impact of urban construction land. As you can see, in terms of the probability rule, arable land and urban construction land are in competition, but from the point of spatial distribution, cultivated land is the best space in central and southern areas, and the near eastern area transition probability is low and water systems are mainly distributed in the central, southern and western states. The water supply in the east is not only relatively weak, but it is also located in the direction leading to the center of Chengdu and a cluster of industrial parks. The suitability probability of urban construction land is relatively smooth, except for the low probability of the western protection area (brown-red part), the other parts of the suitability probability have little difference. The low suitability probability of a small piece of cultivated land and urban construction land in the south is caused by water area factors. Compared with other areas, this area is far from water area, so the overall probability will be slightly affected.

1

**Figure 10.** Spatial expression of various influencing factors (**a**) and influence probability distribution (**b**).

**Figure 11.** (**a**–**d**) Gaussian radial basis function curves of influencing factors of cultivated land. **Figure 11.** (**a**–**d**) Gaussian radial basis function curves of influencing factors of cultivated land. **Figure 11.** (**a**–**d**) Gaussian radial basis function curves of influencing factors of cultivated land.

**Figure 12.** (**a**–**c**) Gaussian radial basis function curves of influencing factors of urban construction land. **Figure 12.** (**a**–**c**) Gaussian radial basis function curves of influencing factors of urban construction land. **Figure 12.** (**a**–**c**) Gaussian radial basis function curves of influencing factors of urban construction land.

According to the characteristics of the Bayesian network, the probabilities of each node are independent from each other, so the joint influence probability of each factor can be calculated, which is used to express the relationship between each factor and is conducive to the transformation of independent influence probability to comprehensive influence probability. Figure 13 shows the spatial distribution of the joint influence probability of cultivated land. Figure 14 shows the spatial distribution of probability of joint impact of urban construction land. As you can see, in terms of the probability rule, arable land and urban construction land are in competition, but from the point of spatial distribution, cultivated land is the best space in central and southern areas, and the near eastern area transition probability is low and water systems are mainly distributed in the central,

According to the characteristics of the Bayesian network, the probabilities of each node are independent from each other, so the joint influence probability of each factor can be calculated, which is used to express the relationship between each factor and is conducive to the transformation of independent influence probability to comprehensive influence probability. Figure 13 shows the spatial distribution of the joint influence probability of cultivated land. Figure 14 shows the spatial distribution of probability of joint impact of urban construction land. As you can see, in terms of the probability rule, arable land and urban construction land are in competition, but from the point of spatial distribution, cultivated land is the best space in central and southern areas, and the near eastern area transition probability is low and water systems are mainly distributed in the central,

*4.3. Affect Probability Distribution*

*4.3. Affect Probability Distribution*

overall probability will be slightly affected.

overall probability will be slightly affected.

*Sustainability* **2021**, *13*, x FOR PEER REVIEW 18 of 25

southern and western states. The water supply in the east is not only relatively weak, but it is also located in the direction leading to the center of Chengdu and a cluster of industrial parks. The suitability probability of urban construction land is relatively smooth, except for the low probability of the western protection area (brown-red part), the other parts of the suitability probability have little difference. The low suitability probability of a small piece of cultivated land and urban construction land in the south is caused by water area factors. Compared with other areas, this area is far from water area, so the

southern and western states. The water supply in the east is not only relatively weak, but it is also located in the direction leading to the center of Chengdu and a cluster of industrial parks. The suitability probability of urban construction land is relatively smooth, except for the low probability of the western protection area (brown-red part), the other parts of the suitability probability have little difference. The low suitability probability of a small piece of cultivated land and urban construction land in the south is caused by water area factors. Compared with other areas, this area is far from water area, so the

**Figure 13.** Probability distribution of joint influence of cultivated land.

**Figure 14.** Probability distribution of joint influence of urban construction land.

#### **Figure 14.** Probability distribution of joint influence of urban construction land. *4.4. Optimization Analysis of Land Use Structure*

**Figure 14.** Probability distribution of joint influence of urban construction land. *4.4. Optimization Analysis of Land Use Structure*  To solve the fuzzy mathematical programming problem in MATLAB R2019a, it is necessary to turn the fuzzy mathematical programming problem into the general programming problem. The parameters of constraint conditions include three types, including environmental parameters, economic parameters, social parameters, planning parameters and parameters of objective function, as shown in Tables 1 and 2. Finally, the value of is 2 × 10−8; ଵ = 436.87 km2, that is, the area of cultivated land is 436.87 km2; ଶ = 69.42 km2, that is, the area of construction land is 69.42 km2; = 98.38 CNY/kg, that is, the economic benefit brought by 1 kg carbon emission is CNY 98.38. The original construction land scale of Chongzhou was 68.5 km2, while the construction land scale is controlled to *4.4. Optimization Analysis of Land Use Structure*  To solve the fuzzy mathematical programming problem in MATLAB R2019a, it is necessary to turn the fuzzy mathematical programming problem into the general programming problem. The parameters of constraint conditions include three types, including environmental parameters, economic parameters, social parameters, planning parameters and parameters of objective function, as shown in Tables 1 and 2. Finally, the value of is 2 × 10−8; ଵ = 436.87 km2, that is, the area of cultivated land is 436.87 km2; ଶ = 69.42 km2, that is, the area of construction land is 69.42 km2; = 98.38 CNY/kg, that is, the economic benefit brought by 1 kg carbon emission is CNY 98.38. The original construction land scale of Chongzhou was 68.5 km2, while the construction land scale is controlled to To solve the fuzzy mathematical programming problem in MATLAB R2019a, it is necessary to turn the fuzzy mathematical programming problem into the general programming problem. The parameters of constraint conditions include three types, including environmental parameters, economic parameters, social parameters, planning parameters and parameters of objective function, as shown in Tables 1 and 2. Finally, the value of *λ* is 2 × 10−8 ; *x*1 = 436.87 km2 , that is, the area of cultivated land is 436.87 km2 ; *x*2 = 69.42 km2 , that is, the area of construction land is 69.42 km2 ; *z* = 98.38 CNY/kg, that is, the economic benefit brought by 1 kg carbon emission is CNY 98.38. The original construction land scale of Chongzhou was 68.5 km2 , while the construction land scale is controlled to 70 km2 , and the optimized result is 69.42 km2 . Therefore, it can be seen from the solution results that the construction land scale of Chongzhou in 2017 is within a reasonable range and can be increased appropriately but should not be too much. The cultivated land protection red line is 338 km2 , and the optimized results show that Chongzhou can have 436.87 km2 of cultivated land. It is suggested that in addition to the reclamation of wasteland, greater efforts can be made in land consolidation and reclamation to make full use of high-quality agricultural resources. The cultivated land protection red line is 338 km2 . According to the calculation, the economic benefit of 1 kg carbon emission in Chongzhou is CNY 98.38, which belongs to the high economic carbon efficiency. In addition, Chongzhou is not a heavy industrial city with key carbon emission departments as its pillar industry. Its agriculture and tertiary industries are developed, and the secondary industry is dominated by green and intelligent industries. The annual carbon sink of Chongzhou is 2.3 × 109 kg, with a strong natural carbon neutralization field, which to a certain extent ensures that the carbon emissions generated in Chongzhou will not cause climate pressure in Wenchuan, Dujiangyan, Dayi and other surrounding areas and can be neutralized by itself.

**Table 1.** Constraint parameters.

**Table 2.** Objective function parameters.

Finally, considering the application of two stage results, based on the parameters obtained from the Bayesian networks and fuzzy mathematical programming, we build the cellular automata (CA) by considering the impact of land use change of three kinds of basic situations, that is, the role of the neighborhood, the areas of conservation should be prevented from developing into urban areas as well as planning restrictions, as mentioned in Section 2.1.3. The spatial simulation of the optimized results is obtained. As shown in Figure 15, the gray area represents the protected areas and water areas that are prohibited from development and utilization, the green area represents the cultivated land, and the red area represents the construction land. By comparing the land use changes before and after simulation, it can be clearly seen that the construction land increment circled in yellow is located in the central and western parts of Chongzhou, as shown in Figure 15a. Correspondingly, the amount of cultivated land in this region will decrease, but the total area of cultivated land has space to increase. From the point of urban agglomeration, the region in the radiation zone between the state and the Chengdu core are more likely to develop into urban construction land, but the zone influenced by publicizing internal stated may also be influenced by external conditions, such as the urban siphon effect. Small cities such as Chongzhou are more susceptible to the siphon effect of large cities such as Chengdu, which is reflected in space.

such as Chengdu, which is reflected in space.

**Figure 15.** (**a**,**b**) Simulation results of land use spatial optimization of Chongzhou in 2017. **Figure 15.** (**a**,**b**) Simulation results of land use spatial optimization of Chongzhou in 2017.

#### *4.5. Behavior Rule Validation 4.5. Behavior Rule Validation Sustainability* **2021**, *13*, x FOR PEER REVIEW 21 of 25

**Table 2.** Objective function parameters.

**Objection Parameters Unit Value**

yuan/km2 4.41 × 108

Economic benefits per unit of land area yuan/km2 8.60 × 106

Carbon emissions from farmland wastewater treatment kg/km2 6.80 × 103 Carbon emissions from municipal wastewater treatment kg/km2 19.85 × 103 Carbon emissions from solid waste disposal kg/km2 62.28 × 103

Carbon emissions from livestock and poultry production kg 9.59 × 107 Carbon emissions from people ploughing land kg/km2 19.85 × 103 Carbon emissions from urban construction land kg/km2 498.68 × 103 Carbon emissions from energy use on cultivated land kg/km2 1.07 × 103 Carbon emissions from urban energy use kg/km2 2.89 × 106

Finally, considering the application of two stage results, based on the parameters obtained from the Bayesian networks and fuzzy mathematical programming, we build the cellular automata (CA) by considering the impact of land use change of three kinds of basic situations, that is, the role of the neighborhood, the areas of conservation should be prevented from developing into urban areas as well as planning restrictions, as mentioned in Section 2.1.3. The spatial simulation of the optimized results is obtained. As shown in Figure 15, the gray area represents the protected areas and water areas that are prohibited from development and utilization, the green area represents the cultivated land, and the red area represents the construction land. By comparing the land use changes before and after simulation, it can be clearly seen that the construction land increment circled in yellow is located in the central and western parts of Chongzhou, as shown in Figure 15a. Correspondingly, the amount of cultivated land in this region will decrease, but the total area of cultivated land has space to increase. From the point of urban agglomeration, the region in the radiation zone between the state and the Chengdu core are more likely to develop into urban construction land, but the zone influenced by publicizing internal stated may also be influenced by external conditions, such as the urban siphon effect. Small cities such as Chongzhou are more susceptible to the siphon effect of large cities

In the 1940s, the Monte Carlo method was proposed. Monte Carlo is an experimental probabilistic method, which is based on a large number of randomized trials to describe In the 1940s, the Monte Carlo method was proposed. Monte Carlo is an experimental probabilistic method, which is based on a large number of randomized trials to describe the complex uncertainty problems encountered in practice [47]. It is a method that can be systematically and comprehensively explored, which is superior to the parameters of isolation, effectiveness analysis and single-class models in performance, especially in the ability to show the coupling error caused by coupling through the Monte Carlo test. Therefore, this study uses Monte Carlo thought to verify the transformation rules and completes simulation and calculation with ArcGIS and Excel. The Monte Carlo method determines the transition state of land use for each unit based on probability, is a more scientific method to deal with uncertain factors of land use. Considering the size of the study area, this study sampled 100,000 random points in the urban area, which can basically cover the whole area closely. Figure 16 shows the distribution of random points in a part of the area. In addition to the reserve, other non-urban cells into urban construction land cells have certain regularity and transformation rules. In order to verify the rationality of the conversion rules obtained in this study, according to the Monte Carlo idea, when the number of random points is large enough, the situation close to the actual situation can be calculated, that is, when the joint influence probability reaches that value, the cell transformation will occur. The verification results show that if policy restrictions are not taken into account, when the probability of the joint impact formulated in this study exceeds 50%, the state of cell can be changed under the condition of natural development and converted into a cell of urban construction land with an accuracy of 91.99%. the complex uncertainty problems encountered in practice [47]. It is a method that can be systematically and comprehensively explored, which is superior to the parameters of isolation, effectiveness analysis and single-class models in performance, especially in the ability to show the coupling error caused by coupling through the Monte Carlo test. Therefore, this study uses Monte Carlo thought to verify the transformation rules and completes simulation and calculation with ArcGIS and Excel. The Monte Carlo method determines the transition state of land use for each unit based on probability, is a more scientific method to deal with uncertain factors of land use. Considering the size of the study area, this study sampled 100,000 random points in the urban area, which can basically cover the whole area closely. Figure 16 shows the distribution of random points in a part of the area. In addition to the reserve, other non-urban cells into urban construction land cells have certain regularity and transformation rules. In order to verify the rationality of the conversion rules obtained in this study, according to the Monte Carlo idea, when the number of random points is large enough, the situation close to the actual situation can be calculated, that is, when the joint influence probability reaches that value, the cell transformation will occur. The verification results show that if policy restrictions are not taken into account, when the probability of the joint impact formulated in this study exceeds 50%, the state of cell can be changed under the condition of natural development and converted into a cell of urban construction land with an accuracy of 91.99%.

**5. Discussion**

**Figure 16.** Schematic diagram of Monte Carlo random sampling. be turned into urban construction land, with an accuracy of 91.99%. Therefore, this pro- **Figure 16.** Schematic diagram of Monte Carlo random sampling.

Based on the current situation and shortage of urban land use structure and spatial

land use structure under uncertainties by using a Bayesian network and fuzzy mathematical programming. Based on the results of the two stages of analysis, the cellular automata simulation is completed. The results indicated that, according to the joint influence probability of arable land and urban construction land, the best land space for cultivated land is in the middle and the south, while the conversion probability of the area near the east is low. Except for the western protection area, the joint impact probability of construction land in the other parts is generally similar. From the perspective of structure, the construction land increment of 1.5 km2 can achieve the control of the construction land scale. After the optimization of fuzzy planning, the optimal construction land scale is 69.42 km2. Under the condition that the cultivated land red line is guaranteed, there is still 98.87 km2 of space for the increase in cultivated land. In order to make full use of high-quality agricultural resources conditions, the obtained parameters were applied to the simulation, which appeared that the increase in construction land would occur in the central and western part of Chongzhou, which may be caused by the urban siphon effect. According to Monte Carlo verification, when the conversion probability exceeds 50%, the cultivated land could

### **5. Discussion**

Based on the current situation and shortage of urban land use structure and spatial arrangement, this study analyzes the influencing mechanism of land use and optimize the land use structure under uncertainties by using a Bayesian network and fuzzy mathematical programming. Based on the results of the two stages of analysis, the cellular automata simulation is completed. The results indicated that, according to the joint influence probability of arable land and urban construction land, the best land space for cultivated land is in the middle and the south, while the conversion probability of the area near the east is low. Except for the western protection area, the joint impact probability of construction land in the other parts is generally similar. From the perspective of structure, the construction land increment of 1.5 km2 can achieve the control of the construction land scale. After the optimization of fuzzy planning, the optimal construction land scale is 69.42 km2 . Under the condition that the cultivated land red line is guaranteed, there is still 98.87 km2 of space for the increase in cultivated land. In order to make full use of high-quality agricultural resources conditions, the obtained parameters were applied to the simulation, which appeared that the increase in construction land would occur in the central and western part of Chongzhou, which may be caused by the urban siphon effect. According to Monte Carlo verification, when the conversion probability exceeds 50%, the cultivated land could be turned into urban construction land, with an accuracy of 91.99%. Therefore, this proposed framework is helpful to understand the process of land use and provides a reference for making scientific and reasonable territorial spatial planning and guiding land use practice under uncertainties.

Under conventional constraints, this study uses a small entry point to express one aspect of the uncertainty, namely, the uncertainty of the target conditions for economic carbon efficiency. The data were compared with Van et al.'s (2012) simulated land supply with known drivers of land use, and then the uncertainty was analyzed using the Monte Carlo method. Meyer et al. (2012) used a Bayesian network to help assess the uncertainties brought about by stakeholder knowledge modeling. This study considers the uncertainty in the optimization of land use structure and applied the results obtained by the uncertainty method to the spatial simulation.

As shown in the above results, the two-stage land use analysis framework constructed in this study has a good optimization effect. According to the analysis of spatial influencing factors, in the spatial pattern, cultivated land and urban construction land constitute a competitive relationship, and cultivated land is at a competitive disadvantage. If the construction land is not controlled, the agricultural space is bound to cause pressure. From the perspective of land use structure, Chongzhou can realize the growth of construction land and cultivated land area under this optimization framework, but the spatial distribution may change dynamically. Taking construction land as an example, this study simulated its spatial distribution, and the results shows that the increase could occur in the central and western part of Chongzhou, which may be caused by the radiation force of various conditions inside Chongzhou and the radiation force of the external core area of Chengdu. Although the development of the protected area is restricted, it does not constitute the conditions for the development of the urban area. Based on the analysis results of land use in Chongzhou, it can be seen that the spatial influence factors can objectively reflect the law of land use in Chongzhou, and the driving effect of all kinds of spatial influence factors on land use in Chongzhou is significant. The constructed Bayesian networks can capture the uncertainty of spatial influencing factors, especially in the case that all land use information cannot be captured at present, which highlights its advantages in applicability. In addition, fuzzy mathematical planning can supplement the role of land use planning on land use to a certain extent, quantify the impact including economy, environment and policy, help to understand the development of land use in an all-round way and provide some support for optimizing land use so as to promote the orderly, controllable and diversified development of cities. Based on the results of this study, this paper will provide some thoughts for urban land use from the following aspects: Firstly, based on the concept

of sustainable development, more consideration should be given to the environmental benefits and carbon emissions of land use under uncertain and flexible conditions, so as to formulate more promising land plans for cities. Second, based on the influencing factors of land use, planning should correctly grasp the law of land use development, strictly protect high-quality land and guard against the expansion of inefficient use of construction land. The third aspect is to rationally plan functional areas according to different geographical locations of urban space, give full play to regional resource endowments and advantages and formulate more targeted land schemes.

Due to the limitation of current technology and foundations, there is still room for improvement and optimization in this study. The goal of this study is to provide support for the decision-making system, which can reflect the information and elements needed in the decision-making process to a certain extent. The Bayesian network and fuzzy mathematical programming method can capture more uncertainty, which is verified to be correct to in some sense, but the construction of the model is limited to the primary stage and cannot reflect all information and components. The real system is often much more complex than the idea, so it needs to be further explored to clarify more details inside the model, which have been rarely probed so far. The current challenges also include the quantification of external influencing factors, the accessibility of spatial data and the reasonable setting of fuzzy parameters. The uncertainty of land use will continue to exist in the future. Although, considering specific problems, the choice of models and methods may be different, it is still necessary to ensure that models and methods do not have a dominant impact on the research results, and more attention should be paid to practical problems.

#### **6. Conclusions**

From the perspective of land supply-side structural reform, this study took cultivated land and urban construction land in Chongzhou in 2017 as the main research objects, carried out a two-stage study on land use optimization and drew the following conclusions. The spatial distribution of joint impact probability of cultivated land and urban construction land shows that the best land space of cultivated land is in the middle and south, while the regional conversion probability near the east is low, which may be caused by the distribution of water sources and rainfall. The joint impact probability of urban construction land is relatively smooth. Except for the probability of the western reserve being low, the joint impact probability of other parts is generally similar. Compared with other areas, this area is far from the water area, so the overall probability is slightly affected. From the perspective of spatial performance, the form of urban construction land is relatively concentrated, while cultivated land is located in the periphery of urban construction land. In the long run, urban construction land has the dynamic characteristics of "catching up" with cultivated land. If the scale of urban construction land is not controlled and high-quality cultivated land resources are not protected, there will be a risk of imbalance in urban spatial structure, and the space of forest land, irrigated grassland and protected areas will be more tense. In addition to the reclamation of wasteland, Chongzhou can also make greater efforts in land consolidation and reclamation to make full use of high-quality agricultural resource conditions. Construction land can be appropriately increased, but it should not exceed the scale control of construction land. The carbon emissions from production, construction and living in Chongzhou can be fully absorbed by the local carbon sink, which can be self-sufficient, and the carbon emissions are well controlled. Therefore, Chongzhou plays an important role as an ecological conservation area in Chengdu. The results show that the increase in construction land can occur in the central and western parts of Chongzhou, which may be caused by the radiation force of various conditions inside Chongzhou and the external radiation force of the core area of Chengdu. The Monte Carlo verification results show that when the conversion probability is more than 50%, the cultivated land can be turned into urban construction land, with an accuracy of 91.99%, passing the test. In addition, the study provides a methodological reference for the formulation of more sustainable and ecologically friendly land schemes, explores the

land use impact mechanism and land use impact under uncertain conditions, expands the perspective and method of land use impact mechanism research, further deepens the understanding of uncertainty in land use process and provides certain theoretical significance for land use development research.

**Author Contributions:** J.Y. and B.Q. conceived and designed the research; A.D. and B.Q. collected, managed, A.D. verified the data; M.Z. and S.L. calculated and analyzed the data and the results; J.Y. and M.Z. wrote the manuscript. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was supported by National Natural Science Foundation of China No. 41401631.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare that they have no conflict of interest.

## **References**

## *Article* **Multi-Scenario Simulations of Land Use and Habitat Quality Based on a PLUS-InVEST Model: A Case Study of Baoding, China**

**Nan Hu 1 , Dong Xu 1 , Ning Zou 1 , Shuxin Fan 1 , Peiyan Wang 2 and Yunyuan Li 1,\***

**\*** Correspondence<EMAIL_ADDRESS>
**Abstract:** Habitat quality and ecosystem service value (ESV) are important foundations for sustainable development. Baoding, as the strategic hinterland of Beijing–Tianjin–Hebei, is of great significance to regional ecological conservation and sustainable urban development. Based on land-use data from 2000 to 2020, the land-use scenarios of natural development (ND), water protection (WP), forest rehabilitation (FR), and cultivated land protection (CP) in 2030 were predicted by the PLUS model and adopt the InVEST model and equivalent ESV table to assess ecological sustainability. The results show that: (1) From 2000 to 2020, the construction land in Baoding has increased by 812 km2 , and the cultivated land and forest land decreased by 708 km2 and 154 km2 . Habitat quality is obviously deteriorating in 4.66% of the city. (2) Under different scenarios, the order of habitat quality is CP > FR > WP > ND. The habitat quality under each scenario is dominated by medium habitat quality. (3) Under different scenarios, the order of ESV is FR > CP> WP > ND. The fluctuation of forest land and cultivated land scale is affecting the ESV. (4) CP and FR will form a land-use pattern that has "high ecological quality and value", which better balances the economic development and ecological protection of Baoding. This research study will provide a reference for the effective allocation of land resources and will guide the formulation of urban land space planning policy in Baoding.

**Keywords:** land use; habitat quality; ecosystem service value; PLUS-InVEST model; multi-scenario simulations

## **1. Introduction**

To achieve urban–rural integration in the new era, harmonizing regional ecological quality and ecological values is key to promoting regional economic sustainability at the macro level [1,2]. Thus, it is important to accurately assess the future evolution of regional habitat quality trends and calculate the future value of regional ecosystem services. Landuse change represents an important link between human socio-economic activities and natural environment evolution, habitat quality, and ecosystem service values (ESV) [3,4]. Currently, the analysis and prediction of land-use change can be achieved using various models; for example, Logistic-CA [5], CLUES [6,7], and FLUS [8,9], which are based on meta-cellular automata (CA) [10], have been widely used for the spatial optimization of land-use patterns as well as ecological red line delineation. The PLUS model, which is a novel land-use simulation model, exhibits higher simulation accuracy than that of the above-mentioned models and can be better integrated with planning policies to support regional sustainable development studies [11]. For instance, Li et al. coupled the grey multi-objective optimization model (GMOP) and the PLUS model to analyze the spatial and temporal evolution of ESV in the Sichuan–Yunnan region for 2026 under three scenarios: the business-as-usual scenario, the ecological development priority scenario, and the ecological and economic balance scenario [9]. Moreover, Lin et al. obtained the spatial and temporal distributions of carbon stocks in Guangdong Province from 1990 to 2020 by using the InVEST and PLUS models [12], while Lin et al. used the PLUS model for simulating

**Citation:** Hu, N.; Xu, D.; Zou, N.; Fan, S.; Wang, P.; Li, Y. Multi-Scenario Simulations of Land Use and Habitat Quality Based on a PLUS-InVEST Model: A Case Study of Baoding, China. *Sustainability* **2023**, *15*, 557. https://doi.org/10.3390/su15010557

Academic Editors: Qingsong He, Jiayu Wu, Chen Zeng and Linzi Zheng

Received: 18 November 2022 Revised: 22 December 2022 Accepted: 26 December 2022 Published: 28 December 2022

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

land-use and land-cover change in the Fuxian Lake basin [13]. The above-mentioned studies demonstrated the applicability of the PLUS model for the analysis and prediction of land-use change.

Land-use change tends to affect the regional habitat quality, which is the ability of an ecosystem to provide space for sustainable survival and the development of suitable individuals and populations within a certain spatial and temporal scale [14,15]. Stable habitats are the prerequisite and basis for ecosystem services and functions, and they are essential for maintaining and enhancing biodiversity [16]. Habitat quality assessment methods have evolved from single-indicator analyses to multi-indicator integrated measurements and then to distributed model dynamic assessments. With the increasing concern for urban ecological environments and the quality of life, researchers constructed mathematical models such as ARIES [17], SolVES [18], and InVEST [19,20] to carry out a long time-series dynamic assessment, among which the InVEST model based on the distributed algorithm of 3S technology has been widely used. Furthermore, changes in land-use type can affect the structures, processes, and functions of ecosystems, which in turn significantly affect the ESV [21]. Accordingly, scholars studied spatial and temporal changes in ESV in different study areas [22,23], watersheds [24], and ecological zones [25]. Previous studies used the InVEST model to measure habitat quality in different regions, which effectively assesses the impact of changes in natural and social conditions on ecosystem service quality. Consequently, combining the InVEST and ESV will explore the synchronicity of ecological "quality" and "value" to characterize the relationship between ecological development and ecological economy.

With economic development and changes in urban spatial structure, Baoding, as an important city in the Beijing–Tianjin–Hebei city cluster, is facing ecological and environmental problems such as air pollution, water shortage, degradation, and the shrinkage of wetland and loss of biodiversity [26,27], which have become bottlenecks limiting the development of the region. In 2017, it was listed by the Ministry of Housing and Construction as the third batch of pilot cities for "Double Urban Repairs" [28]. To relieve the pressure on the ecological environment of Baoding and promote the development of ecological industries and the ecological integration of Beijing, Tianjin, and Hebei, this paper checks the ecological development-related policy documents of Baoding between 2000 and 2020 and innovatively extracts four policy focuses as the basis for setting scenarios, while combining the results of the scenarios with the highest ecological quality and value for re-validation to provide a reference for the effective allocation of land resources and the spatial planning of land under ecological protection.

### **2. Materials and Methods**

This paper uses the PLUS model and multi-period high-precision land-use data to identify the land-use change patterns in Baoding from 2000 to 2020 and to simulate the landuse situation in Baoding in 2030 under the guidance of four important urban development policies. On this basis, the InVEST model was used to calculate the habitat quality of Baoding in different periods and situations, and the ecological service value equivalence table was introduced to estimate the value of ecosystem services in different situations in Baoding to identify the land-use pattern of "high ecological quality-high ecological value", and to further validate this land-use pattern. The land-use model was further validated to provide reference for the effective allocation of land resources and land space planning under ecological protection (Figure 1).

**Figure 1.** Research Framework.

#### *2.1. Overview of the Study Area*

Baoding is in the western part of the central Hebei province and between the eastern foothills of Taihang Mountains and the western part of the Jizhong Plain (113◦4503200– 116◦1904100 E and 38◦1402900–39◦5703 00 N) (Figure 2). It is adjacent to Beijing and Zhangjiakou in the north, Langfang and Cangzhou in the east, Shijiazhuang and Hengshui in the south, and Shanxi in the west. The topography of the city is high in the northwest and low in the southeast, while the landscape consists of mountains and plains and an uneven regional distribution of river resources. Furthermore, 5 municipal districts, 4 county-level cities, and 15 counties are under its jurisdiction, thereby exhibiting a total area of 22,190 km2 , a resident population of 11,546,000 in 2020, a GDP of CNY 335.33 billion (accounting for 9.38% of the total GDP of Hebei Province) and one of the best economic growth rates in the province [29].

**Figure 2.** Location of the study area.

#### *2.2. Data Source*

The main data sources used in this study include the following: (1) Baoding grain production data and acreage data from the Baoding Statistical Yearbook and grain price data from the National Compilation of Cost and Benefit Information on Agricultural Products [30]; (2) land-use data from the Globeland30 global surface use database (http://www. globallandcover.com, accessed on 10 September 2022), in which the images for the land cover classification of development and update of GlobeLand30 are mainly 30 m multispectral images, including Landsat TM5 ETM+ [31]; (3) elevation data from the geospatial data cloud platform of the Computer Network Information Center of the Chinese Academy of Sciences (http://www.gscloud.cn, accessed on 10 September 2022), with a spatial resolution of 30 × 30 m; (4) temperature and precipitation data from the National Science and Technology Infrastructure Platform at the National Earth System Science Data Centre—Loess Plateau Sub-Center (http://loess.geodata.cn, accessed on 12 September 2022) [32]; (5) GDP and population data from the Resource and Environment Science Data Registration and Publication System; (6) roads, water, and railway data are cited from the National Geographic Information Resources Catalogue Service (www.webmap.cn). To unify the spatial accuracy of the data, the above data were processed by using the cropping and resampling tools of ArcGis 10.5 to convert them into 30 m × 30 m raster files.

#### *2.3. Research Methodology*

#### 2.3.1. PLUS Model

Based on the metacellular automata model, the PLUS model integrates a land expansion strategy analysis module and a simulation model of future land-use changes via multi-class random patch seeding [11]. First, it analyzes the spatial characteristics of various land-use expansion forces and elucidates the driving factors between the two phases of land-use data. Second, it employs the random forest algorithm to sample and calculate the land-use expansion forces one by one to obtain the development probability of each land-use type. Finally, based on the roulette wheel strategy, it combines random patch generation, the transition transfer matrix, and the decreasing threshold mechanism to achieve optimization and to determine the final land-use pattern.

• Selection of land-use change factors

Land-use change factors are affected by multiple factors. On the basis of considering the accuracy and reality of the model, the factors of the natural environment, social economy, and traffic accessibility are comprehensively considered. According to the principles of availability, quantification, and the consistency of driving factors, 10 driving factors were selected (Table 1).

**Table 1.** Land-use change factors.

• Setting neighborhood effects

The neighborhood effects can reflect the expansion capacity of different land types [33]. When applying the PLUS model, it is necessary to determine the ease of interconversion between different land-use types by using neighborhood weights as follows [11]:

$$
\Omega\_{\mathbf{i},\mathbf{k}}^{\mathbf{t}} = \frac{\text{con}\left(\mathbf{c}\_{\mathbf{i}}^{\mathbf{t}-1} = \mathbf{k}\right)}{\mathbf{n} \times \mathbf{n} - 1} \times \mathbf{w}\_{\mathbf{k}} \tag{1}
$$

where con c t−1 i = k represents the total number of grid cells occupied by land-use type k at the last iteration within the n × n window (n is generally taken as 3) [34]; and wk is the weight among the different land-use types between [0, 1]. The neighborhood weights are set based on existing research results and the actual conditions of Baoding (Table 2).

**Table 2.** Neighborhood weights parameters.

• Verification of model accuracy

The PLUS model is often tested using point-by-point comparisons or random validations. Random validations are generally suitable for large-scale simulations, whereas point-by-point comparisons are suitable for relatively small areas [35]. In this study, based on the actual land-use data in the Baoding, this study uses a Kappa coefficient and overall accuracy to test the simulation accuracy of the model. The overall accuracy (OA) refers to the ratio of the number of correctly classified class cells to the total number of classes. The Kappa coefficient is often used to compare the similarity between two images [13]. The expression of the kappa coefficient is as follows:

$$\text{kappa} = \frac{\text{P}\_0 - \text{P}\_\mathbf{c}}{\text{P}\_\mathbf{P} - \text{P}\_\mathbf{c}} \tag{2}$$

where P0 is the observed agreement rate between the reference map and the simulation results; Pc is the expected proportion of correct simulations in the random case; and Pp is the proportion of correct simulations in the ideal classification case and is generally taken as 1.

The results show that the Kappa coefficient of the predicted images in 2020 is 0.869 and the overall accuracy is 0.925, indicating that the PLUS simulation results are spatially consistent and exhibit good applicability. Thus, this approach can be used for future land-use simulation models based on Baoding.

#### 2.3.2. Scenario Setting

Baoding is the ecological environment support area of Beijing–Tianjin–Hebei. According to the characteristics of land-use change in Baoding and the guidance of urban development policies, this study set up 4 scenarios, natural development, water protection, forest rehabilitation, and cultivated land protection, to predict the land-use change in Baoding in 2030. Simulation parameters include the land-use demand, conversion constraint, transition matrix, and neighborhood weights.

1. Natural development: This scenario simulates habitat quality and EVS development until 2030 when the land is not affected by any policies and fully complies with the current situational changes from 2000 to 2020, so no restricted area and land-use change restrictions are set in this scenario. Based on the high-precision land-use data released by the government in 2000, 2010, and 2020, the Markov chain of the PLUS model was used to identify the changing rules for various types of land in the past 20 years, and the identification results were used as the basis for land-use changes in this scenario. It is also the control group.

### 2.3.3. InVEST Model

The habitat quality module of the InVEST model assesses habitat quality, while reflecting the genetic variation and species reproduction potential; this is achieved by considering the sensitivity of each land cover type relative to threat factors as well as the intensity

of external threats [42]. Consequently, the statuses of different resources and conditions in the environment are determined for the survival and development of individuals or populations. The relative impacts of each threat and the distance between the habitat grid and the threat need to be calculated before habitat quality can be calculated using the InVEST model. The formula is as follows:

$$i\_{rxy} = 1 - \left(\frac{d\_{xy}}{d\_{r\,\max}}\right) \text{if } \text{linear} \tag{3}$$

$$d\_{\rm rxy} = \exp\left(-\left(\frac{2.99}{d\_{r\,\max}}\right)d\_{\rm xy}\right) \text{ if exponential} \tag{4}$$

where *dxy* is the linear distance between grid cells *x* and *y*, and *dr max* is the maximum effective distance of threat *r*'s reach across space.

Habitat degradation is calculated as follows:

$$\mathbf{D}\_{\mathbf{x}\mathbf{j}} = \sum\_{\mathbf{r}=1}^{\overline{\mathbf{R}}} \sum\_{\mathbf{y}=1}^{\mathbf{Y}\_{\mathbf{r}}} \left( \frac{\mathbf{w}\_{\mathbf{r}}}{\sum\_{\mathbf{r}=1}^{\overline{\mathbf{R}}} \mathbf{w}\_{\mathbf{r}}} \right) \mathbf{r}\_{\mathbf{y}} \mathbf{i}\_{\mathbf{r}\mathbf{y}} \beta\_{\mathbf{x}} \mathbf{S}\_{\mathbf{j}\mathbf{r}} \tag{5}$$

where y indexes all grid cells on r's raster map, and Yr indicates the set of grid cells on r's raster map.

The final measure of biodiversity and ecosystem function is based on a habitat quality index, which is calculated as follows:

$$\mathbf{Q}\_{\mathbf{x}\mathbf{j}} = \mathbf{H}\_{\mathbf{j}} \left( 1 - \frac{\mathbf{D}\_{\mathbf{x}\mathbf{j}}^{\mathbf{z}}}{\mathbf{D}\_{\mathbf{x}\mathbf{j}}^{\mathbf{z}} + \mathbf{k}^{\mathbf{z}}} \right) \tag{6}$$

where Qxj is the habitat quality index for land-use type j in grid x, and Hj is the habitat suitability for land-use type j, Dxj is the habitat degradation for land-use type j in grid x, k is the half-saturation constant, and z is the default parameter of the model (usually taken as 2.5).

Since the sensitivity of different land-use types varies toward habitat quality threats, only the impact of human activities on the habitat is considered; the more sensitive the landuse unit, the lower its ability to resist disturbances and the more severe the degradation. Habitat suitability and its sensitivity to different threat factors were determined based on relevant information regarding Baoding and previous studies in the neighboring areas [43] (Table 3). Four types of areas with high population activity were selected as stressors, namely: cultivated land, construction land, roads, and railways (Table 4).

**Table 3.** Sensitivity parameters of different land types to habitat threat factors.

**Table 4.** Threat factor parameters.

#### 2.3.4. Calculating the ESV

The ESV in Baoding was determined by referring to the model proposed by Costanza et al. [21]. Then, it was combined with the ESV equivalent factor proposed by Xie et al. [22,23]. The standard equivalent is the economic value of the natural food produced by one hectare of cropland in a year at the national average yield level. It is generally considered to be equivalent to one-seventh of the economic value of food produced on a hectare of cropland at market prices [44,45]. Table 5 shows the ESVs of grain yield per unit area of the ecosystem in Baoding (1567.76 CNY/hm2/a).

$$\text{ESV} = \sum\_{\mathbf{f}=1}^{\mathbf{m}} \sum\_{\mathbf{i}=1}^{\mathbf{n}} (\mathbf{C}\_{\mathbf{f}} \times \mathbf{E}\_{\mathbf{f}\mathbf{i}} \times \mathbf{A}\_{\mathbf{i}}) \tag{7}$$

where ESV is the ecosystem service value in CNY/yr, i is the land-use type and n is the number of land-use types (n = 7 in this study), f is the ecosystem service type and Ai is the area of land-use type i in hm2 . Efi is the equivalent value of ecosystem service f for land-use type i, m is the number of ecosystem service types (m = 11 in this study), and Cf is the economic value of 1 unit of ecosystem services in CNY hm−2yr−1 . This equation can also be represented as follows:

$$\mathbf{E} = \frac{1}{7} \sum\_{\mathbf{i}=1}^{n} \frac{\mathbf{m}\_{\mathbf{i}} \mathbf{p}\_{\mathbf{i}} \mathbf{q}\_{\mathbf{i}}}{\text{TotalArea}} \tag{8}$$

where E is the economic value of 1 unit of ecosystem services, i is the grain crop type, mi is the average price of grain crop i in Baoding, pi is the yield of grain crop i in kg/hm2 , qi is the planting area of grain crop i in hm2 , and TotalArea is the total area planted in grain crops (hm2 ).

**Table 5.** ESVs per unit area of the ecosystem (CNY/hm2 ).

**3. Results**

#### **3. Results** During 2000–2020, the land-use of Baoding was dominated by cultivated land,

*3.1. Analysis of Land-Use Change* forestland, and construction land, with the sum of the areas of these land types accounting

3.1.1. Characteristics of Land-Use Change from 2000 to 2020

*Sustainability* **2023**, *15*, x FOR PEER REVIEW 10 of 19

3.1.1. Characteristics of Land-Use Change from 2000 to 2020 for over 90%, cultivated land and forestland accounted for 49% and 34% of the total area

During 2000–2020, the land-use of Baoding was dominated by cultivated land, forestland, and construction land, with the sum of the areas of these land types accounting for over 90%, cultivated land and forestland accounted for 49% and 34% of the total area on average per year, respectively (Table 6). Over the past 20 years, forestland in the middle- and low-lying mountainous areas of Baoding was gradually converted to cultivated land and construction land (Figure 3). Meanwhile, the grassland area in the hilly region, which acts as a border between the mountainous area and the plain area, slightly increased. Overall, the expansion of construction land in the plain area is prominent, and the main urban area shows a clear trend of expansion toward the north. on average per year, respectively (Table 6). Over the past 20 years, forestland in the middle- and low-lying mountainous areas of Baoding was gradually converted to cultivated land and construction land (Figure 3). Meanwhile, the grassland area in the hilly region, which acts as a border between the mountainous area and the plain area, slightly increased. Overall, the expansion of construction land in the plain area is prominent, and the main urban area shows a clear trend of expansion toward the north. **Table 6.** Proportion of land-use types in Baoding from 2000 to 2020 (%).

> **Table 6.** Proportion of land-use types in Baoding from 2000 to 2020 (%). **Type of Land-Use**

**Figure 3.** Spatial distribution of land-use from 2000 to 2020.

**Figure 3.** Spatial distribution of land-use from 2000 to 2020. Cultivated land and forestland account for the highest share of the total area of Baoding. However, the area of cultivated land and forestland slightly increased during 2000–2020, which was followed by a significant reduction; that is, cultivated land decreased from 11,090 km² in 2000 to 10,380 km² in 2020, while forestland decreased from 7510 km² in 2000 to 7360 km² in 2020. Accordingly, there is a dire need to conserve cultivated land and forestland. In 20 years, the area of construction land grew by 812 km², which represents the highest increase in the area of any land type. Meanwhile, the area of wetland has decreased by 0.53 km², and the area of water areas has increased by 5 km²; this contrasting trend shows that some wetland gradually transformed into water lands, Cultivated land and forestland account for the highest share of the total area of Baoding. However, the area of cultivated land and forestland slightly increased during 2000–2020, which was followed by a significant reduction; that is, cultivated land decreased from 11,090 km2 in 2000 to 10,380 km2 in 2020, while forestland decreased from 7510 km2 in 2000 to 7360 km2 in 2020. Accordingly, there is a dire need to conserve cultivated land and forestland. In 20 years, the area of construction land grew by 812 km2 , which represents the highest increase in the area of any land type. Meanwhile, the area of wetland has decreased by 0.53 km2 , and the area of water areas has increased by 5 km2 ; this contrasting trend shows that some wetland gradually transformed into water lands, and that the ecological protection of the wetland is necessary. Grassland and shrubland slowly increased over these 20 years by 41 km2 and 4 km2 , respectively. Overall, the land-use patterns in Baoding show a significant increase in construction land, a significant decrease in cultivated land and forestland, and a slight shift in other land types.

#### and that the ecological protection of the wetland is necessary. Grassland and shrubland 3.1.2. Land-Use Simulation under Different Scenarios

slowly increased over these 20 years by 41 km² and 4 km², respectively. Overall, the landuse patterns in Baoding show a significant increase in construction land, a significant decrease in cultivated land and forestland, and a slight shift in other land types. 3.1.2. Land-Use Simulation under Different Scenarios Land-use simulations for 2030 were carried out based on the neighborhood weights in Table 2 and the transfer cost matrix under different scenarios. Compared with the 2020 land-use data, the change in the land-use areas in Baoding in 2030 varied under different scenarios (Table 7); the specific scenario settings and corresponding results have been mentioned below.

Land-use simulations for 2030 were carried out based on the neighborhood weights

land-use data, the change in the land-use areas in Baoding in 2030 varied under different

**Table 7.** Proportion (%) and decadal change (km2 ) of land-use types in Baoding under different scenarios in 2030.

Natural development scenario: The increase in construction land area is higher than the other three scenarios, and cultivated land and forestland are still the main contributors to the increase in construction land area. Under this scenario, the forestland area in Baoding decreased by 130 km2 , and the cultivated land area decreased by 541 km2 (becoming the largest land type to be transformed). Meanwhile, the area of construction land increased by 636 km2 , whereas the area of other land types only changed slightly. Based on the spatial distribution of land-use changes (Figure 4), there is a tendency for cultivated land to transform into forestland in the northwestern part of the mountainous area and at the edge of the nature reserve in the southwestern part of Baoding. Nevertheless, construction land in the northern part of the hilly area, the northeastern part of the plain area, and the central city will continue to expand outward along the urban edges, whereas construction land in the southern part of the plain area is likely to transform into cultivated land.

**Figure 4.** Spatial distribution of land-use change under different scenarios in 2030.

Water protection scenario: The land-use change trend in this scenario is similar to ND, the reduction in forestland and cultivated land has not been effectively curbed, but the transfer of ecological land is lower and the range of change is smaller. Meanwhile, the water body and wetland will, respectively, increase by 13 km2 and 2 km2 , because a series of ecological restoration and comprehensive management measures have been carried out in the upstream rivers of Baiyangdian Lake, the water area along the rivers and lakes in the mountainous and hilly areas has increased.

Forest rehabilitation: The area of cultivated land in Baoding was further compressed compared with the ND, with a reduction of 838 km2 , which was the largest reduction in the four scenarios. The expansion of construction land was still obvious, but its expansion rate was effectively controlled from 3% in the ND to 2% because the development of construction land will also be restricted to a certain extent when returning cultivated land to forest. Meanwhile, cultivated land in mountainous and hilly areas (mainly concentrated in the southwest with complex terrain and developed water systems) transferred to woodland.

Cultivated land protection: Under the CP, land use was maintained correspondingly with the land-use situation in 2020. The area of cultivated land remained unchanged, the area of forestland increased by 156 km2 , the area of grassland decreased by 164 km2 , and the area of construction land only increased by 8 km2 . With respect to the areas of land-use change concentrated along the southwest waters of mountainous areas, cultivated land protection manifested as the transfer of grassland to forestland at the edge of the water body and the transfer of grassland to construction land at the junction of hilly area and plain area, and the expansion of construction land in plain areas no longer has the phenomenon of large-scale diffusion due to the strict restriction of cultivated land protection systems and policy guidance with respect to national stock development.

#### *3.2. Analysis of Habitat Quality Change*

#### 3.2.1. Characteristics of Habitat Quality Change from 2000 to 2020

Based on previous studies on habitat quality and the actual situation in Baoding [46], the natural breakpoint method was used to classify habitat quality into five categories: very low (0.0–0.1), low (0.1–0.4), medium (0.40–0.68), high (0.68–0.90), and very high (0.90–1.0). Table 8 shows that the habitat quality in Baoding is relatively stable between 2000 and 2010, exhibiting an average change of 67 km2 or <1% for all land types. However, during 2010–2020, the changes in habitat quality increased sharply, with an average annual change of 413 km2 being observed for all land types. In particular, the low habitat quality area increased by 827 km2 , the medium habitat quality area decreased by −874 km2 , and the habitat quality in other areas changed relatively slightly (1% on average). Therefore, a continuous transformation trend of medium habitat quality areas to low habitat quality areas was observed.

**Table 8.** Proportion of each habitat quality category in Baoding from 2000 to 2020 (%).

To further analyze the changes in habitat quality, the spatial distribution of habitat quality in Baoding was calculated (Figure 5). The land-use type is mainly woodland and the ecological environment faces decreased human disturbance, thereby exhibiting strong stability. Areas with low and medium habitat quality are mainly concentrated in mountainous counties and towns because these contain large construction, agricultural, and forestlands; since they exhibit a strong heterogeneity, habitat degradation is higher. Habitat quality in hilly areas is mainly influenced by the land-use types controlled by

natural factors within the area; for example, water and grasslands dominate the land-use types along the irrigation areas and reservoirs in the southwestern part of Baoding, with patches of woodland being fragmented and interspersed with grasslands, thereby resulting in a low overall habitat quality. In plain areas, due to the significant increase in urban land-use and human activity, habitat quality is degrading and has been mainly categorized as low and medium.

**Figure 5.** Spatial distribution of habitat quality from 2000 to 2020.

The ESV in Baoding from 2000 to 2020 showed a decreasing trend year by year, which was CNY 340, 339 and 332 billion. Forestland has the highest ESV, with an average value of CNY 230 billion (accounting for 68% of the total value); this indicates that although forestland is not the largest land type in Baoding, it plays a more important role in ecological protection and development. Cultivated land, which is the largest land type in terms of area, has an average ESV of CNY 68 million (accounting for 20% of the total value). Because of economic, social, and regional development, cultivated land is more prone to encroachment. Furthermore, it plays the role of restraining the uncontrolled expansion of construction land, which has a significant positive impact on ecological development. Variations in the ecological value of water, and grassland are low, indicating that changes in these regions have less influence on the ESV in Baoding. Therefore, when formulating relevant policies, government departments can consider improving or protecting forestland and cultivated land to improve the ecological quality of Baoding's regional environment to maintain its ecosystem.

## 3.2.2. Habitat Quality Simulation under Different Scenarios