Ecological security pattern-based simulation for land use structure change: a case study in Ezhou City, China

Abstract

The ecological environment quality is an important constraint and an optimization objective for land resource allocation. Integrating ecological service value (ESV) accounting and ecological security pattern (ESP) delineation, and combining with the land use structure of 2004/2010/2016 in Ezhou City, this research laid out the urban ESP based on ESV with Net Primary Productivity (NPP), and made it as the main influence factor to simulate land use structure in 2022. The results indicated that: 1) The water body has the biggest contribution to ESV, while the construction land has the minimum; 2) 91 ecological corridors are extracted, of which 28 were important ecological corridors; there were 36 ecological nodes extracted, including 17 important nodes; 3) According to ESV, Ezhou City was divided into four security zones. The area of ecological restoration zone was the largest, and human activity core zone area was the smallest; 4) In the no ESP protection scenarios and ESP protection scenarios separately, the net increase area of construction land is from 868.5 hm² to 52.74 hm² in the ecological core protection area; the construction land in the human activity core area has been increased by 2342.31 hm² in protected scene, 766.23 hm² more than that of the unprotected scene. The results show that the division of security zones promoted the relocation of construction land from ecological protection core zone to human activity core zone, which can protect the ecological environment effectively, and the ESP-based simulation can provide the decision-making reference to coordinate the relationship of regional land resource allocation and the ecological environment protection.

Keywords

simulation of land use structure / Net Primary Productivity / ecological security pattern / ecosystem service value / Ezhou City

1 Introduction

Land use structure is the spatial distribution and quantitative combination of various types of land resources under certain input and constraint conditions. As a product of natural factors and human activities, land use structure change is with complicated process and evolution mechanism, which is also the competition result of different land use types under various driving forces (Li et al., 2003). Current research is focusing on establishing the land use structure optimization theory, simulating the land use structure with appropriate model, and optimizing the future land use pattern (Yan and Ji, 2001). A number of models and algorithms, such as CA-Markov model (Kityuttachai et al., 2013; Xu et al., 2017), probit regression model (Liu et al., 2015), cellular automata model (Wang et al., 2012) have been conducted on the land use structure simulation and predication in different regional scales. The land carrying capacity (Zhu et al., 2014) and land resource optimization allocation (Wang et al., 2012; Yuan and Liu, 2014) are analyzed to reveal the relationship between land use change and the driving factors in terms of nature, resource, population and economy (Pan et al., 2010; Munteanu et al., 2014; van Vliet et al., 2015). Ecological factor is also concerned and used as an important constraint and an optimization objective for the land resource allocation (Feng et al., 2009; Zhao et al., 2011). Under different ecological environment conditions, the framework for land use optimization models is proposed (Lin and Liu, 2002), and optimization method of land use structure is explored effectively (Xiang and Meng, 2013). In these studies, some ecological security constraint elements are considered in the process of land use change (Yu et al., 2014), and ecological source factors (Wei et al., 2016) such as water and woodlands were used for the land use simulation. The spatial-temporal pattern of land ecological security was evaluated with some math models (Meng et al., 2014; Yu et al., 2017). The interaction mechanism between land use change and ecological environment quality is explored, in which the importance and restriction of ecological environment in land use was illustrated (Song et al., 2015; Du et al., 2017; Wang et al., 2018).

Ecological security is a state in which the ecological and environmental resources in a country or a region are not restricted or threatened by resources and the ecological environment (Ma et al., 2013). More attention is paid to the ecosystems service functions on human social systems (Li et al., 2017; Verhagen et al., 2018; Xu et al., 2019). Based on multi-source data, including statistical data, remote sensing data and monitoring data (Raudsepp-Hearne et al., 2010; Luo et al., 2013), index system is built to assess ecological security. A lot of different models are invented, such as pressure-state-response model (Bai and Tang, 2010), dynamic systems model (Li et al., 2015), gray relational model (Yu et al., 2012; Wang et al., 2013; Yang et al., 2013; Zhang et al., 2014), to analyze the spatial-temporal evolution (You et al., 2011; Yang et al., 2013; Ye et al., 2014) of ecological security and divide ESP (Cong et al., 2018). Then, the land ecological security and the quality of agricultural ecological environment is analyzed (Gao et al., 2014; Zhang et al., 2015; Su et al., 2016). These researches also discussed the relationship between regional ecological environment and urban development, integrated ecological security into urban development planning decision, and combined the ecosystem pattern and ecological evolution to analyze ecological early warning, urban expansion and land space planning (Zhao and Ma, 2014; Liu et al., 2018). Examples include the study that analyzed the impact of ecological barriers in urban expansion with the minimum cumulative resistance model (MCR) (Zhang et al., 2017); the study that analyzed the urban form and spatial pattern in the process of urban growth under the various kinds of influence factors, such as water resources, geological disasters and biodiversity (Tang et al., 2009; Ren et al., 2013); the study that defined the ecological bottom line of urban expansion based on ecosystem services and security patterns, and the study that delineated the urban growth boundary through the ESP and urban expansion simulation (Zhou et al., 2014a). These researches focused on the harmonious development between urbanization and ecological security, and emphasized the restraint effect of ecological environment on urban development (Hodson and Marvin, 2010; Zhou et al., 2014b; Yu et al., 2017).

A lot of researches in the land use structure change simulation and ecological security have been carried out based on the identification of external conflicts and pressures in the ecological safety system. These studies reveal a close relationship between land use change and ecosystem services triggered by urbanization. For land resources allocation and utilization, as well as understanding and handling the forming process and impact mechanism of land use structure, the systematic study of the relationship between ecological services and land use change is of great theoretical and practical significance. The relationship of the spatiotemporal variation of land use types with the ESV accounting and the ESP construction has been discussed in many literatures. However, the ESV that is used as the key constraint condition to make ESP for land use structure change simulation lacks sufficient practical discussion. The methods that mainly use panel data supplemented by remote sensing data and the administrative divisions as the evaluation unit can reflect the overall situation. Since lack of field monitoring data and the mismatch of data scale, the spatial characteristics of ESP and its internal diversity of administrative units can’t be accurately reflected.

Ezhou City is an important open port city in the middle and lower reaches of the Yangtze River Economic Belt in China. It is also an important part of Wuhan City Circle and a demonstration area of ecological civilization in Hubei Province, China. In the process of fast urbanization, dramatic economic and spatial restructuring result in tremendous land use/land cover change. These impacts lead to continuous deterioration and losses in the ecological values of fragile environment. Therefore, the contradiction and conflicts between the urban development and the ecological environmental protection appears very prominent, and it is urgent to explore some solution. By setting the ESP as a constraint to simulate the land use pattern in different scenarios, we put forward the best land use scheme, hoping to provide a discussion for the land use research under the background of ecological civilization construction, and provide a corresponding method reference for urban ecological environment protection and sustainable use of land resources. In this study, we adopted a variety of methods to perform a spatiotemporal analysis of land use structure during 18 years (2004–2022), and the research is designed with the objectives of (i) mapping NPP by CASA model Carnegie-Ames-Stanford Approach model) and estimating ESV based on NPP; (ii) dividing ESP based on ESV; (iii) mapping spatiotemporal land use status in 2004/2010/2016, and simulating land use structure in 2022 under different kinds of scenarios.

2 Study area and materials

2.1 Study area

Ezhou City is located in the southeastern part of Hubei Province, with east longitude 114°32'E‒115°05′E, north latitude 30°00'N‒30°06′N, and the city covers 1596 km² in 2017 (Fig. 1). Its geomorphic unit includes four types: Yangtze River alluvial terrace, hilly landform, hilly plain and stagnant alluvial plain. The average annual rainfall is 1282.8mm, the average annual evaporation is 1520 mm, and the average temperature is 17°C as it has a subtropical monsoon climate. Since it is known as “The city of hundreds of lakes” and “The land of fish and rice” , Ezhou City has a good natural geographical condition, and the ecological environment is suitable for the growth of various plants since its fertile soil and rich water resource. The soil types include red soil, purple soil, tidal soil and paddy soil. The soil organic matter content range of the whole city is 12.8‒55.7 g/kg, and the pH range is 5.28–8.1. The natural vegetation is mainly trees, shrubs, and artificial cultivation of economic forest is mainly tung tree, rape, chestnut, and walnut.

2.2 Data sources and processing

The data sets used in this paper mainly include land cover data, NDVI data, gridded meteorological data, and field soil data.

The basic data for the land use structure dynamics analyses was Landsat (2004, 2010, 2016) images from “Data Center for Resources and Environment Sciences, Chinese Academy of Sciences” (available at Resource and Environment Science and Data Center website). The preliminary results of land use data were obtained by visual interpretation on remote sensing images, and then the final land use data were obtained by modifying them based on field survey data. Then the land use structure was divided into 6 types, namely: cultivated land, garden land, woodland, water area, construction land and others land use. The types description and its correspondence with current land classification standards in China are shown in the table below (Table 1).

The Normalized Difference Vegetation Index (NDVI) is downloaded from the MODIS vegetation index product MOD13Q1 data set from USGS website, and its spatial resolution is 250 m with a temporal resolution of 16 days. The 46 month-to-month images of Ezhou City in 2016 were selected and the images were processed by MRT (MODIS Reprojection Tool). After data mosaic and combination, re-projection, the maximum value of the month data was verified, and the NDVI of January to December was obtained. Landsat TM data to get a monthly NDVI data with a spatial resolution of 30m are integrated.

The meteorological data are from China National Meteorological Information Center. It mainly included monthly mean precipitation, monthly mean temperature data, monthly solar total radiation data, and the geographical location information of each monitoring site (coordinates and elevation) of 11 monitoring sites around the study area in 2016. These data were interpolated by Kriging method with ArcGIS^®. The grid size and coordinate system of the meteorological data are consistent with remote sensing data.

Soil data was got from the agricultural bureau of Ezhou City. Some data were supplemented with field sampling in July 2016 and September 2017, and the sampling depth is 0–30 cm. Test items included the contents of total nitrogen (N), total phosphorus (P) and total potassium (K) using soil agrochemical routine method.

3 Research methods

The framework of the paper is shown in Fig. 2. First, the land use structure of the study area is extracted, of which spatiotemporal characteristics are analyzed. Secondly, CASA model is used to calculate NPP, and ESV is obtained from nine services types based on NPP. Thirdly, the ESP is given with MCR model according to ESV. Finally, CA-Markov model is used to simulate the land use structure in the near future with two scenarios.

3.1 Methods of ESP Delineation

3.1.1 NPP estimate—CASA model

NPP is the total amount of organic matter accumulated by green vegetation per unit area of a unit of time, which is defined as the difference between organic matter produced by photosynthesis and organic matter consumed by autotrophic respiration. There are three kinds of model for NPP estimation: statistical model, parametric model and process model. The paper use CASA parameter model to simulate it (Yu et al., 2009). CASA model is a model of light use efficiency in which remote sensing data, meteorological data and vegetation types are used as the input parameters. The NPP is mainly determined by photosynthetic active radiation (APAR) absorbed by plants and light energy utilization (ε) (Potter et al., 1993). The NPP can be calculated by the equations followed:

(1)\mathrm{N}\mathrm{P}\mathrm{P}(x,t)=\mathrm{A}\mathrm{P}\mathrm{A}\mathrm{R}(x,t)\times \epsilon (x,t),

(2)\mathrm{A}\mathrm{P}\mathrm{A}\mathrm{R}(x,t)=\mathrm{S}\mathrm{O}\mathrm{L}(x,t)\times \mathrm{F}\mathrm{P}\mathrm{A}\mathrm{R}(x,t)\times \mathrm{0.5},

(3)\epsilon (x,t)=T_{\mathrm{l}\mathrm{o}\mathrm{w}}(x,t)\times T_{\mathrm{l}\mathrm{i}\mathrm{g}\mathrm{h}\mathrm{t}}(x,t)\times W_{\epsilon }(x,t)\times {\epsilon }_{\max },

where x is the pixel location and t is the time (in months); \mathrm{A}\mathrm{P}\mathrm{A}\mathrm{R}(x,t) is the photosynthetic active radiation absorbed by the canopy of vegetation at pixel x in month t (MJ·m⁻²); \epsilon (x,t) is the light use efficiency of the vegetation at pixel x at month t (gC ·MJ⁻¹), and is obtained by correcting the maximum utilization of light energy with the temperature stress factor and the water stress factor; \mathrm{S}\mathrm{O}\mathrm{L}(x,t)is the total solar radiation received by pixel x in month t (MJ·m⁻²); \mathrm{F}\mathrm{P}\mathrm{A}\mathrm{R}(x,t) is the proportion of photosynthetic active radiation absorbed by vegetation at pixel x in month t, and it is estimated based on NDVI; T_{\mathrm{l}\mathrm{o}\mathrm{w}}(x,t) and T_{\mathrm{l}\mathrm{i}\mathrm{g}\mathrm{h}\mathrm{t}}(x,t) are respectively the stress coefficients of light energy utilization at low temperature and high temperature; W_{\epsilon }(x,t)is the water stress coefficient; {\epsilon }_{\max }is the maximum light energy utilization of different vegetation in the ideal states (gC·MJ⁻¹).

3.1.2 Estimation of ESV

Referring to the classification and estimation methods of ESV by Costanza et al. (1997), Xie et al. (2015), Ouyang et al. (1999), etc., and considering the topography and landform characteristics, vegetation and soil conditions, and the status quo of social and economic development of the research area, the ESV, in this study, is classified into the following nine services types: raw material production function, nutrient cycling function, gas regulation function, climate regulation function, hydrological regulation function, ecological environmental purification function, soil formation and protection function, biodiversity maintenance function, and recreational and culture function.

1) Value estimation of raw material production function

The most basic and important function of the ecosystem is that green vegetation can transform CO₂ and H₂O into organic matter through photosynthesis, which can provide the most primitive energy for human beings and other organisms. The quantity of the organic matter produced by the ecosystem is estimated by energy substitution method and the equation below:

(4)V_{\mathrm{o}\mathrm{m}}={\displaystyle \sum _{x=\mathrm{1}}^{n}NPP(x)\times P_{\mathrm{o}\mathrm{m}}},

where V_{\mathrm{o}\mathrm{m}} is the value of material production, yuan·(m²·yr) ⁻¹; NPP(x) is the quantity of organic matter produced annually by pixel x, gC·(m²·yr)⁻¹; P_{\mathrm{o}\mathrm{m}} is the price of organic matter, yuan.

2) Value estimation of nutrient cycling function

N, P, K, and other nutrients in organic environment could be transformed into organic matter by photosynthesis, which can provide nutrients to maintain the ecological balance. The value of nutrient cycling function can be estimated by

(5)V_{\mathrm{c}}={\displaystyle \sum _{x=\mathrm{1}}^n{V}_{\mathrm{c}i}(x)}={\displaystyle \sum _{x=\mathrm{1}}^{n}NPP(x)R_{i\mathrm{1}}\times R_{i\mathrm{2}}\times P_i},

where V_{\mathrm{c}} is the value of nutrient cycling, yuan·(m²·yr)⁻¹; V_{\mathrm{c}i}(x)is the value of nutrient elements i accumulated on the unit area of the pixel x, yuan·(m²·yr)⁻¹; i represents three kinds of element, N, P, and K; R_i1 is the distribution rate of elements i in organic matter produced by different ecosystems, %; R_i2 is the proportion of elements converted into chemical fertilizer, %; and P_i is the price of fertilizer i, yuan·t⁻¹.

3) Value estimation of gas regulation function

Gas regulation function mainly refers to plants emit O₂ during photosynthesis and adjust the air quality. The estimation formula is as follows:

(6)V_{\mathrm{s}\mathrm{f}\mathrm{o}}={\displaystyle \sum _{x=\mathrm{1}}^n\mathrm{1.2}}\times NPP(x)\times V_{{\mathrm{o}}_{\mathrm{2}}},

where V_{\mathrm{s}\mathrm{f}\mathrm{o}} is the value of gas regulation, yuan ·(m²·yr)⁻¹; V_{{\mathrm{o}}_{\mathrm{2}}}is the unit price of industrial oxygen, yuan ·gC⁻¹.

4) Value estimation of climate regulation function

The value of climate regulation is reflected in the amount of CO₂ absorbed by plants. The formula is as follow:

(7)V_{\mathrm{x}\mathrm{s}\mathrm{c}}={\displaystyle \sum _{x=\mathrm{1}}^n\mathrm{1.62}\times NPP(x)\times V_{\mathrm{c}{\mathrm{o}}_{\mathrm{2}}}},

where V_{\mathrm{x}\mathrm{s}\mathrm{c}} is the value of climate regulation, yuan·(m²·yr) ⁻¹ and V_{\mathrm{c}{\mathrm{o}}_{\mathrm{2}}} is the price per unit of CO₂, yuan ·gC⁻¹.

5) Value estimation of hydrological regulation function

The hydrological regulation function includes water regulation and water supply. The former means the ecosystem regulation conservation function for regional water (including the water flow in the lakes and rivers) when the underlying surface is water. The latter means when the underlying surface is soil, the ecosystem provides filtration, maintenance and storage of water for vegetation and various types of organisms in the soil. Since the regulating function of the ecosystem to the water is similar to the type of reservoir water storage, the alternative engineering method to estimate the ecosystem service value of water conservation is used based on the average storage cost of the reservoir construction. The formula is as follow:

(8)V_{\mathrm{w}\mathrm{c}}={\displaystyle \sum _{x=\mathrm{1}}^{n}V(x)\times P_{\mathrm{w}}},

where V_{\mathrm{w}\mathrm{c}} is the water conservation value, yuan·(m²·yr)⁻¹; V(x) is the annual water conservation capacity per unit area of pixel x, m³·(m²·yr)⁻¹; P_w is unit capacity cost of reservoir construction, yuan·m⁻³.

When the underlying surface is soil, the formula is as follows:

(9)V_{\mathrm{s}}(x)={\displaystyle \sum _{x=\mathrm{1}}^n{P}_{\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{n}}(x)\times K_{\mathrm{w}}\times R_{\mathrm{w}}},

where V_{\mathrm{s}}(x) is the amount of water conservation capacity per unit area of pixel x when the underlying surface is soil; P_{\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{n}}(x) is the monthly precipitation in pixel x, mm·m⁻²; K_{\mathrm{w}} is the ratio of runoff rainfall to total rainfall; R_{\mathrm{w}} is the runoff reduction coefficient of surface vegetation.

When the underlying surface is water, the formula is as follows:

(10)V_{\mathrm{w}\mathrm{c}}(x)={\displaystyle \sum _{x=\mathrm{1}}^n{P}_{\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{n}}(x)-E{T}_a(x)},

where V_{\mathrm{w}\mathrm{c}}(x) is the amount of water conservation capacity per unit area of pixel x when the underlying surface is water; and E{T}_a(x)is the actual monthly evaporation at pixel x, mm·m⁻².

6) Value estimation of ecological environmental purification function

Ecological environmental purification refers to the reduction of pollutants and toxicity in the environment through the metabolism of biological groups, and it mainly includes absorption of pollutants, blocking of dust, killing of bacteria and reduction of noise. In this paper, the absorption and decomposition of SO₂、HF、NO_x, and adsorption of dust by various land use types were evaluated. The formula is as follows:

(11)V_{\mathrm{c}\mathrm{e}}={\displaystyle \sum _{a=\mathrm{1}}^{\mathrm{9}}{\displaystyle \sum _{b=\mathrm{1}}^{\mathrm{4}}{Q}_{ab}\times C_{ab}\times S_a}},

where, V_{\mathrm{c}\mathrm{e}} is the value of ecological environmental purification, yuan·(m²·yr)⁻¹; Q_{ab}is the amount of harmful substances b absorbed by the land use type a, g·(m²·yr) ⁻¹; C_{ab} is the unit governance costs of land use type a to harmful substances b, yuan·g⁻¹; S_ais the area of land use type a, m²; a is different land use types, the value is 1‒9; b is the different harmful substances, the value is 1‒4.

7) Value estimation of soil formation and protection function

Remote sensing data are used to estimate soil conservation, and soil formation and protection value is calculated by opportunity cost method. Considering the seasonal dynamics of precipitation and vegetation and their different coupling effects with seasonal changes, USLE model was used to estimate soil protection amount, and the monthly model equation modified that is proposed by Qi et al. (2011):

(12)A={\displaystyle \sum _{m=\mathrm{1}}^{\mathrm{12}}[K\times LS\times R_m\times (\mathrm{1}-P\times C_m）]},

where A is soil conservation amount, t·hm⁻²; R is the factor of rainfall-off erosivity factor, based on the annual average precipitation and the monthly average precipitation, Weschmeier’s empirical formula (Mccool et al., 1982) was used to calculate the precipitation, MJ·mm·(hm²·h)⁻¹; K is the soil erodibility factor, Williams’ empirical formula (Williams et al., 1983) was used to calculate it, (t·hm²·h)·(hm²·MJ·mm)⁻¹; L is the slope length, S is the slope direction; C is cover-management factor, it was assigned according to vegetation coverage (Cai et al., 2000), the value range is 0–1, dimensionless; P is soil and water conservation factor, including engineering measures and tillage measures, based on the empirical P values adopted by the United States and China under different soil and water conservation measures, it was assigned by the slope percentage calculated by DEM and the superposition analysis of land use data, dimensionless; m is the month.

Then, the formula for soil formation and protection value is as follows:

(13)V_{\mathrm{g}}={\displaystyle \frac{A}{p_{\mathrm{b}}\times h}}\times v,

where, V_g is the value of soil formation and retention function; P_b is the soil bulk density, 1.35 g/cm³ according to “Soil erosion classification standard SL190-2007”; h is the thickness of soil layer, v is the annual yield of soil per unit area.

8) Value estimation of biodiversity maintenance function

The Shannon-wiener index in within-habitat diversity was used to calculate the biodiversity maintenance value:

(14)H=-{\displaystyle \sum _{t=\mathrm{1}}^m{P}_t\times {\log }_{\mathrm{2}}{P}_t},

where H is within-habitat diversity, P_t is the proportion of individuals belonging to species t in all individuals, and m is the number of species. According to the survey data provided by the agricultural bureau and forestry bureau of the research area, the biodiversity index is calculated and interpolated in the study area. Then the biodiversity index is converted into currency value by referring to the graded value of species conservation index of “Forest Ecosystem Service Function Evaluation Standard of China”.

9) Value estimation of recreational and culture function

In this study, the recreational and culture value are reflected through tourism income, since the research area have the superior geographical conditions and rich ecological tourism resources, such as Xishan scenic spot, Lianhua mountain scenic spot, Liangzi lake eco-tourism resort, Honglian lake resort, Geshan scenic spot, White pheasant mountain scenic spot, Longpan ji, Wuwang city, Taiping mountain tourist area and so on. In 2016, the total income of cultural tourism in Ezhou City was 5.267 billion yuan. The ratio of recreational cultural value of each ecological type cultivated land: garden land: forest land: water area: others should be 1: 3.06: 12.2: 26.1:3.3, the value of other types is not included in this study (He et al., 2016).

The sum of the above nine values was the ESV determined in this study. According to each value evaluation model, the functional value of each grid is calculated in ArcGIS^®.

3.1.3 Method of modeling ESP‒MCR model

MCR is a derivative of the cost distance model. It is originally applied to the study of species diffusion process, and was tried in land use, tourism planning in recent years. It represents the work done against the resistance from the “source” to the destination. Three main factors of source, distance and base resistance is considered in the model. The model is expressed as

(15)\mathrm{M}\mathrm{C}\mathrm{R}=f_{\min }{\displaystyle \sum _{j=n}^{i=m}{D}_{ij}\times R_i},

where f is the unknown positive function and represents the positive correlation between minimum cumulative resistance and ecological process; D_{ij} represents the space distance of the species from source j to the spot cell i, and R_i represents the resistance coefficient.

In this study, the source selection is based on the demarcation of ecological red line in Ezhou City, and the patch area larger than 5 hm² in the ecological red line area serves as the ecological source area. Such places included the nature reserves, forest, wetland, scenic spots, drinking water sources, important waters, etc., and the total area is 461.68 km² (Table 2).

The process of land use change is the process of “source overcoming obstacle.” In a way, the resistance surface can reflect the trend of the types change, and the magnitude of resistance is the key to judge this trend. In this study, the base resistance surface is established based on the distribution of ESV. As a whole, the higher the value of the unit ecological service is, the less resistance is, and vice versa.

This research identified the ESP based on the land use map, ecological source map, ESV distribution map, and the accumulated resistance value map.

3.1.4 Method of building ecological corridor

The migration of organisms in different ecological sources needs to overcome the landscape resistance. Landscape resistance means the difficulty of migration in different source. The MCR model can be used to simulate the migration of organisms in different sources and calculate the optimal migration route. Five specific steps are included when using MCR model: a) determine the ecological sources and ecological resistance factors, b) assign the weight value foreach resistance factor, c) form the resistance surface, d) generate the cost distance grid data, e) construct ecological corridor (Wang et al., 2019).

Several data sources collected in 2016, including DEM, slope, distance to road, distance to township center, distance to water area and land use data, as ecological sources and ecological resistance factor. According to the scoring method presented by Huang and Chen (2014) and Huang et al. (2019), we obtain the resistance value and weight coefficient of six different resistance surfaces (Table 3). The weights of slope, altitude, distance to water area, distance to construction land, distance to town center and distance to road are 0.15, 0.11, 0.14, 0.17, 0.21, and 0.22 respectively. Among the six resistance factors, weight of the distance to road is the highest, which has the greatest impact on the cumulative resistance value. It is due to the frequent human activity near the road. In the most area, altitude difference is no more than 100 m. Thus, the slope ranges from 0 to 10 degrees. Also, note that water occupies a large area. Hence, the weights of these three factors are relatively low, and have small impact on the cumulative resistance value. Six resistance surfaces are superimposed to obtain the comprehensive resistance surface, then the potential ecological corridor is obtained by using the cost distance and cost path tools in ArcGIS.

Gravity model can be used to calculate the interaction intensity between two different sources. (Wu and Wang, 2015; Shi et al., 2018; Liang and Zhao, 2020) The model is expressed as

(16)G_{\mathrm{xy}}={\displaystyle \frac{W_x{W}_y}{{{\displaystyle N}}_{xy}^{\mathrm{2}}}}={\displaystyle \frac{\left[{\displaystyle \frac{\ln (S_x)}{R_x}}\right]\left[{\displaystyle \frac{\ln (S_y)}{R_y}}\right]}{{{\displaystyle \left({\displaystyle \frac{T_{xy}}{T_{\max }}}\right)}}^{\mathrm{2}}}}={\displaystyle \frac{{{\displaystyle T}}_{\max }^{\mathrm{2}}\times \ln (S_x)\times \ln (S_y)}{{{\displaystyle T}}_{xy}^{\mathrm{2}}\times R_x\times R_y}}.

We attributed the interaction force greater than 1200 to the important ecological corridor, and the rest to the general ecological corridor. Where G_xy is the force between source x and source y; W_x and W_y are the weight values of the two sources; N_xy is the normalized value of potential corridor resistance; S_x and S_y are the area of source x and y, respectively; R_x is the resistance value of source x; T_xy is the accumulated resistance value between source x and y; while T_max is the maximum accumulated resistance value among all corridors. When source x and y have strong interaction force, the G_xy value will be large. Also, it means that the corridor has high importance.

3.2 Simulation method of land use structure

Before the simulation, LULC about its spatial pattern, the conversion matrix in the future situation, and the spatial variation in the forecast year should be analyzed. CA-Markov model is used to carry out the spatial and temporal simulation of land use change, since it has the advantages of both cellular automata theory and Markov model theory on time series and space prediction.

3.2.1 Influence factor selection and weight determination

Based on the principles of accessibility, representativeness, spatio-temporal consistency, six factors are selected as the key influencing factors to reflect the structure change: elevation, slope, distance from the water bodies, distance from the road traffic, distance from the administrative center, distance from itself.

Here, elevation factor and slope factor are important geographical driving factors for land use change. Elevation often had a significant impact on the land use patterns, and is divided into five ranges according to the topographic and geomorphological features of the study area: 10–50, 50–85, 85–140, 140–242, 242–485. Based on the superposition analysis results with the current land use status, the areas with large changes are mainly concentrated in the range of 10–140. Slope was extracted from elevation, and it is divided into five grades: 1–5. By superposition with the current land use situation map, it is shown that the land use types in slopes 4 and 5 are stable. The land use changes mainly took place in the grades of 1 to 3.

The study area is close to the Yangtze River, and there are many lakes, among which Liangzi Lake, Wutong Lake, and Red Lotus Lake are famous important aquatic products areas at home and abroad. Therefore, water factor is an important driving factor influencing the spatial layout. The paper takes 500 m as the interval to make buffer zone for water body, and the result showed that land type changed frequently within the range of 4500 m around the water area.

For local economic development, roads traffic factor is a fundamental driving force. Four high-speed ways run through in the study area. The land use types that are close to these roads are affected by the distance from roads. By establishing a buffer zone with an interval of 500 m for the roads, the analysis shows that the land use change within the range of 3000 m around the road was greatly affected by road traffic factor driving.

The administrative center is an important institution in China, which have a significant impact on the regional pattern. For the factor of the administrative center, buffer zones are built with the center of 25 township government locations. The results showed that the areas with frequent land-class conversion were concentrated in the 3000 m area around the administrative center.

The land use situation itself also has impact on the change of land use type. For example, it is easier to reclaim the land around the cultivated land and expand the construction land on the edge of the city. Therefore, the buffer analysis of the current land use situation is carried out. It is believed that the area with frequent land conversion is within the scope of 100 m, while the frequency of spontaneous conversion is significantly reduced over 300 m range.

Since each influence factor has different impact on the land use change, one of the key questions for simulation is to allocate the weight value of influence factors reasonably. Assuming that the influence factor and land use structure change are not related, the structure change per unit condition should be roughly the same at different periods, which was close to the average of the area variation. But, in fact, the variation of each type of land is different with its the average variation value. The larger the difference is, the greater the effect of the influence factors will be. Therefore, the variation coefficient to descript the correlation between different influence factors and structure changes is used as a basic weight for simulation. The variation can be calculated as

(17){\sigma }_i={\displaystyle \frac{\mathrm{1}}{\overline{x_{ij}}}}\sqrt{{\displaystyle \frac{{\displaystyle \sum _{j=\mathrm{1}}^n{(x_{ij}-\overline{x_{ij}})}^{\mathrm{2}}}}{n_i}}},

where {\sigma }_i is the correlation, i is influence factor, j is an interval of the influence factor, n is the number of intervals, x_i{}_j is variable in the interval, \overline{x_{ij}} is the average variation value.

3.2.2 Creating the land suitability atlas

The land suitability atlas includes two parts: limiting factor atlas and influence factor atlas. The limiting factor atlas is the source map, which is the unoccupied land use type during the study period. The influence factor atlas includes ESP and six influence factors mentioned above. Fuzzy membership function in Multi-criteria evaluation (MCE) is used to determine the land suitability atlas. With normalizing the suitability factors, the MCE model is used to generate land use maps. The land suitability atlas is synthesized by collection editor model in IDRISI^®.

3.2.3 Land use pattern simulation

After the stack of the land use maps of 2004 and 2010, the area and probability of the cell transfer to other land types are obtained. With Markov transition matrix, the land use structure in 2016 is simulated based on the structure in 2010. In the simulation, the model iteration number is 6, 5×5 contiguity filter is used as the cellular neighborhood regular in CA. The real land use structure in 2016 is used to testify the simulation results.

3.2.4 The simulation

The scenarios of ecological protection scenarios (ESP protection scenarios) and no ecological protection scenarios (no ESP protection scenarios) are simulated, and the land use structure in 2022 is simulated based on the initial state of land use structure in 2016. In ESP protection scenarios, the simulation focused on the ecological protection of land status; the ecological source is considered as a key obstacle factor, and is normalized to a Boolean image where the ecological source was 0 and the other was 1. In no ESP protection scenarios, only the six factors are considered.

4 Results and analysis

4.1 Soil data analysis

Based on data from 202 soil samples, we got the content of total nitrogen (N), total phosphorus (P), and total potassium (K) in the soil. The content of total N in soil of the study area is between 374.17 and 376.59 (unit: mg/kg), and the maximum value appears in the southern region. The content value of total P increases gradually from 412.19 to 454.64 in the direction of NW‒SE. The content of total K ranges from 1.99% to 2.26%. No significant difference of the content of N appears in the study area. The sampling points are shown in the Fig. 3.

4.2 The land use change in 2004‒2016

Supplemented by field survey and change survey data of land bureau, we draw the land use status map in 2004/2010/2016 from remote sensing data by object-oriented method in ENVI5.1 (Fig. 4), The land use structure and its changes are obtained in ArcGIS^®, as shown in Table 4. We verified the accuracy of interpretation results based real land use data, and the kappa coefficients of the three periods of data were 89.64%, 92.19%, and 87.38%.

The area of construction land grew fastest from 2004 to 2016. Especially in 2010 to 2016, the increasement reaches to 58.26%. At the same time, cultivated land and forest land showed an opposite trend, among which the cultivated land decreased the most by 23.31% and 18.59% respectively in the two periods. Note that the water area increased by 7.69% from 2004 to 2010, and decreased by 20.22% from 2010 to 2016. Then the decrease of woodland and other land is slight, and the area of garden plots remains basically unchanged. In sum, the expansion of construction land and the contraction of cultivated land and water area are the main characteristics of land use change in these two periods, which means that significant urban expansion in Ezhou City occurred in the past two decades, and the pressure of cultivated land protection and ecological environment protection is increasing.

4.3 The estimates of ESV

According to the Eqs. (1)–(3), the NPP of Ezhou City in 2016 is calculated as shown in Fig. 5. The total NPP in 2016 amounted to 4173.65 billion gC, and the NPP per unit area is between 0–741 gC·(m²·yr) ⁻¹, with an average NPP of 261.49 gC·(m²·yr) ⁻¹. The areas with lower NPP values are concentrated in the water area and urban areas, such as Liangzi Lake, Yanglan Lake and the Yangtze River zone. The area with high NPP value tends to concentrate in the forest land and surface vegetation area. Except for the above part, the NPP distribution is relatively uniform. Comparing the NPP results with the land use status and other data, we can find that the areas with high NPP are mostly with high vegetation coverage, and has high precipitation, high accumulated temperature, and high solar radiation.

Based on the nine indexes calculated from NPP, the ESV per unit area of Ezhou City in 2016 is traced (Fig. 6), and the ESV contribution rate of land use types are measured (Table 5).

Water area is the biggest contributor to the ESV per unit area reaching 33.31%, and its hydrological regulation value and environmental purification value are prominent than any other categories. Forest land is the second for its contribution rate of ESV per unit area reaching 21.89%, and its hydrological regulation value and biodiversity maintenance value are remarkable, which is in consistent with the good ecological environment. Construction land has the lowest value with 0.5573 yuan·(m²·yr)⁻¹ and its contribution rate is only 2.21%. That means construction land makes little contribution to the environment in those nine ESV types. The contradiction between the regional ecological environment protection and the construction land expansion is more prominent. The ESV per unit area of cultivated land, garden land and other land use decrease in turn. Most of the areas with high ESV were important water areas and woodlands, and the ESV results were highly coincident with the important ecological source date, indicating that the ESV results had high credibility.

4.4 Regional assignment of ESP

To simplify the corridor network structure and better show the connection characteristics between the source areas, we combine the small patches adjacent to the large ecological source areas and obtain 14 ecological source areas. The MCR model is used to obtain the minimum cost path between different sources, and 91 ecological corridors are obtained (Fig. 7). The interaction matrix is calculated by the Gravity model, and the result of the interaction matrix can be used to judge the relationship between different sources (Table 6).

As shown in Table 6, the maximum value of the interaction force appears between the source 7 and 8, which is as high as 19104.3. It means that the ecological cost from source 7 to source 8 is the lowest and the corridor between the two sources is easy to build. So, measures need to be taken to protect the corridor and strengthen the links between the two sources. The minimum force value is 50.221, which appears between source 10 and source 12. It means that the source connection between source 10 and 12 is weakest. In addition to the influence of distance between two sources, it is also related to the larger ecological resistance caused by the large scope of construction land in the eastern region.

Then we extracted the intersection points of ecological corridor as ecological nodes, and total 36 ecological nodes were extracted. We chose the points on the important ecological corridors as the important nodes, while others were general ecological nodes, including 17 important nodes and 19 general ecological nodes. In the process of ecological civilization city, we should focus on strengthening the construction of important ecological corridors and the protection of important ecological nodes.

According to the MCR model, the cumulative resistance in the research area is calculated (Fig. 8). Based on cumulative minimum resistance results, considering the distribution of source areas, land use structure and ecological service value, the ecological security pattern of Ezhou city is divided into four grades according to the principle of ecological priority.

Figure 9 shows the ecological security level distribution in 2016. The core area of ecological protection covers an area of 413.62 km², accounting for 25.91% of the total area of Ezhou City. It mainly includes Liangzi Lake Nature Reserve, important water sources and other areas with high ecological service value. As the core area for maintaining ecosystem balance, it has the greatest impact on the regional ecological environment. The ecological buffer zone, covered an area of 332.18 km², accounting for 20.81% of the total area, is located at the periphery of the core protected area. Compared with the core area of ecological protection, the ecological service value of the core area of human activity is slightly lower. The structure in these areas is relatively stable. It is an important component of the ecological environment and can be used as a backup area for ecological protection. The ecological restoration area that has little impact on the city’s ecological environment, covers an area of 602.85 km², and it is the reserve resource area of the anthropic zone and ecological buffer zone reserve area. Human activity core area, an area of 247.82 km², accounting for 15.52% of the total area, mainly exists in city center and development zone. Since there are high-intensity human activities in this area, it is an important supporting area for the economic development. Those areas also have the lowest ecosystem service value, with frequent changes in land use types. And the changes of ecosystem service are most volatile in those areas.

4.5 Simulation and accuracy verification

In the Markov module of IDRISI^®, the land use data in 2004 and 2010 is added. The time interval and prediction interval are 6. The probability of one land use type conversion to other land use types during 2004–2010 is obtained as shown in Table 7. Specifically, in the table, the row represents the probability of one class turning to other classes. The column represents the probability of other classes to be converted to the class.

The land use structure in 2010 is used as the base year to simulate the land use structure in 2016. The difference of quantity and spatial distribution between the simulation map and the status map are evaluated using Crosstab in CA-Markov model.

The evaluation result shows that value of kappa is 0.9038. It is higher than the conventional value (kappa= 0.75). Therefore, this indicates that there is a high consistency and small differences between the simulated map and the actual map. The analysis result shows that the model has high reliability and could be applied to the simulation. Besides, we randomly selected some plots with different land use types in the simulation results and conducted field investigation in June 2021. The results show that 84% of the current land use types are consistent with the simulation results.

Taking the land use structure map in 2016 as the initial state, the simulation time interval is 6, and the land use structure maps in 2022 under two scenarios are simulated, as shown in Fig. 10. The net growth area of construction land with the two scenarios are shown in Table 8.

In both scenarios, the spatial distribution of construction land is different although they have the same total area with 39530.79 hm². In the ecological core protected zone, the net increment of construction land is 868.5 hm² in scenario without ESP protection, and −52.74 hm² in the scenario with ESP protection. The quantity of construction land remains stable in ESP protection scenario, which indicates that the original construction land is moved out and there was no new construction land in this area, and the ecological source area is fully protected. In the ecological buffer zone, the net increase of construction land shows obvious increase in ESP protection scenario than NO ESP protection scenario, which is 2860.29 hm² and 3378.21 hm² respectively. The quantity difference between the two simulation results is up to 517.92 hm². There are two main possible sources, one is that the area within its own region have increased under the influence of various factors, the other is that part of the land originally distributed in the ecological protected core zone is also transferred to the ecological buffer zone. In the human activity core zone, construction land in NO ESP protection scenario is an increase of 1576.08 hm², conversely, 2342.31 hm² in ESP protection scenario, an increase of 766.23 hm² compared with the former. It means economic activity in human activities core zone is more frequent. Therefore, the constraint of ESP played a significant role in guiding the transfer of construction land from ecological protected core zone to human activity core zone, which could protect the regional ecological environment, and guarantee the quantity demand of urban development for construction land. While under the scenario of ESP protection, the construction land in the ecological core area was completely moved out, and construction land area in the other three types increased a lot. Meanwhile, the efficient protection of ecological corridors and nodes helps to the construction of urban ecological network structure, and it also helps to the construction of eco city structure.

5 Conclusions

In this study, we selected CASA model to calculate NPP in 2016, calculated nine index data related to ESV, added them together to obtain ESV in Ezhou city, gained the ecological corridors using MCR and Gravity model, extracted the ecological nodes from ecological corridors, and then evaluated the land use structure in 2022 under two kinds of scenarios.

Among the ESV per unit area in Ezhou City in 2016, the value of water area is 8.4087 yuan·(m²·yr)⁻¹, accounting for 33.31% of the total, and its hydrological regulation/ecological environmental purification function are prominent. It shows that water resources are the most important part of the ecological environment of Ezhou City, and play an extremely important coordinating role in ecological environmental system. Therefore, it is necessary to strengthen the protection and supervision of water resources. Meanwhile, we should take measures to make rational and efficient use of water resources. The value of construction land was the lowest, only 0.5573 yuan·(m²·yr)⁻¹ and the contribution rate of ESV was only 2.21%, which also verified the intense conflict between urban development and ecological environment protection. It indicates that the regional development should find a balance between the protection of ecological environment and the expansion of construction land.

Based on ecological source data and comprehensive ecological resistance data, the ecological corridors in the study area were obtained. 28 important ecological corridors and 63 general ecological corridors were selected using Gravity model. The corridor between source 7 and source 8 is very important, and methods need to be taken to strengthen the connection between them. Then we extracted 17 important ecological nodes and 19 general ecological nodes from the intersection of ecological corridor. In the future, it is necessary to pay attention to the protection of ecological nodes.

According to the ESV, Ezhou City is divided into four safety-level zones. The ecological core protected area is the core area to maintain the regional ecosystem balance, with the powerful ecological function and an area of 413.62 km². The structure of the ecological buffer zone which is located at the periphery of the core protected area, is stable relatively, and has not a significant environmental impact on Ezhou City, with an area of 332.18 km². The ecological restoration area has a small impact relatively on the ecological environment. It is a buffer area for human activities and ecological buffer areas, with an area of 602.85 km². With intensive development of human activities, the core area of human activities that is mainly concentrated in the downtown is an important supporting area for the urban economic development, and has the lowest ecosystem service value, with an area of 247.82 km².

From the simulation results of ESP protection scenario, the expansion of construction land is totally constrained in the ecological core protected areas. Conversely, the construction land within the core area of human activities, which has more intensive and frequent human activities, has a tendency of continuous expansion. Therefore, it can be considered that the division of ESP has a good restrictive effect on the spatial distribution and development direction of the construction land expansion. Compared the simulation results of the two scenarios, the result under ESP protection scenario had a more reasonable land use distribution. It means that our result under ESP protection scenario has high credibility. This paper constructed ecological security pattern in Ezhou City, trying to create a reasonable layout of ecological space and urban construction space. Under the scenario of ESP protection, the construction land was completely removed from the ecological core area and the ecological area was totally protected. Ecological protection is particularly important under the background of sustainable development. Simultaneously, economic development should focus on improving land use efficiency rather than land use scale. How to efficiently layout and use the existing land should be the focus of attention. This study provides reference of land use change simulation research with ecological environment constraint, and a decision-making support for coordinating the relation between urbanization process and ecological environment protection in other developing cities.

References

[1] Bai X R, Tang J (2010). Ecological security assessment of Tianjin by PSR model. Procedia Environ Sci, 2(6): 881–887

[2] Cai C F, Ding S W, Shi Z H, Huang L, Zhang G Y (2000). Study of applying USLE and geographical information system IDRISI to predict soil erosion in small watershed. J Soil Water Conserv, 14(2): 19–24

[3] Cong D M, Zhao S H, Yu T, Chen C, Wang X B (2018). Urban growth boundary delimitation method integrating comprehensive ecological security pattern and urban expansion simulation—a case study of planning areas in Tianshui City (2015‒2030). J Nat Resourc, 33(1): 14–26

[4] Costanza R, d’Arge R, de Groot R, Farber S, Grasso M, Hannon B, Limburg K, Naeem S, O’Neill R V, Paruelo J, Raskin R G, Sutton P, van den Belt M (1997). The value of the world’s ecosystem services and natural capital. Nature, 387(6630): 253–260

[5] Du Y Y, Hu Y N, Yang Y, Peng J (2017). Building ecological security patterns in southwestern mountainous areas based on ecological importance and ecological sensitivity: a case study of Dali Bai Autonomous Prefecture, Yunnan Province. Acta Ecol Sin, 37(24): 8241–8253

[6] Feng Y X, Luo G P, Yin C P, Zhou D C, Lu L, Xu W Q, Dai L (2009). Change of land-use degree and ecological security assessment of inland river basins in the arid area—a case study on Manas River Basin in Xinjiang. J Nat Resourc, 24(11): 1921–1932

[7] Gao Q, Shi X Y, Zhang C, Zhang M R, Ma Y W (2014). Dynamic assessment and prediction on quality of agricultural eco-environment in county area. Transact Chinese Soc Agr Eng, 30(5): 228–237

[8] He L, Jia Q J, Li C, Zhang L, Xu H (2016). Land use pattern simulation based on ecosystem service value and ecological security pattern. Transact Chinese Soc Agr Eng, 32(03): 275–284

[9] Hodson M, Marvin S (2009). ‘Urban ecological security’: a new urban paradigm? Int J Urban Reg Res, 33(1): 193–215

[10] Huang L M, Chen J F (2014). Suitability evaluation of urban construction land based on features extraction of a MCR surface. Resourc Sci, 36(7): 1347–1355

[11] Huang M Y, Yue W Z, Feng S R, Cai J J (2019). Analysis of spatial heterogeneity of ecological security based on MCR model and ecological pattern optimization in the Yuexi County of the Dabie Mountain Area. Resourc Sci, 34(4): 771–784

[12] Kityuttachai K, Tripathi N, Tipdecho T, Shrestha R (2013). Ca-Markov analysis of constrained coastal urban, growth modeling: Hua Hin seaside city, Thailand. Sustainability, 5(4): 1480–1500

[13] Li C, Zhang F G, Song N P, Kong X B, Chen H W (2003). Discussion on theory and method of optimal regional allocation of land use structure. Geogr Geo-Inform Sci, 19(2): 52–55

[14] Li F Z, Lu S S, Sun Y N, Li X, Xi B Y, Liu W Q (2015). Integrated evaluation and scenario simulation for forest ecological security of beijing based on system dynamics model. Sustainability, 7(10): 13631–13659

[15] Li Y, Feng Y, Guo X R, Peng F (2017). Changes in coastal city ecosystem service values based on land use—a case study of Yingkou, China. Land Use Policy, 65: 287–293

[16] Lin Z P, Liu X N (2002). A case study on land use pattern under ecological security in Ecotone between agriculture and animal husbandry in northeastern China. Chinese J Eco, 21(6): 15–19

[17] Liu J H, Li W F, Zhou W Q, Han L J, Qian Y G (2018). Scenario simulation balancing multiple objectives for optimal land use allocation of the Beijing-Tianjin-Hebei megaregion. Acta Ecol Sin, 38(12): 250–259

[18] Liu K, Li Y E, Wu Q, Shen J F (2015). Driving force analysis of land use change in the developed area based on Probit regression model: a case study of Nanjing City, China. Chinese J App Eco, 26(7): 2131–2138 in Chinese)

[19] Luo N S, Li J J, Luo F Z (2013). Empirical analysis on the relationship between the China urbanization and regional eco-efficiency. China Populat, Resourc Environ, 23(11): 53–60

[20] Liang Y Y, Zhao Y D (2020). Construction and optimization of ecological network in Xi’an based on landscape analysis. Chinese J App Ecol, 31(11): 3767–3776in Chinese)

[21] Ma F J, Liu J T, Eneji A E (2013). A review of ecosystem services and research perspectives. Acta Ecol Sin, 33(19): 5963–5972

[22] Mccool D K, Wischmeier W H, Johnson L C (1982). Adapting the universal soil loss equation to the Pacific Northwest. Transact ASAE, 25(4): 0928–0934

[23] Meng J J, Yan Q, Xiang Y Y (2014). The optimization of ecological security pattern based on land use and assessment of schemes in Ordos. J Desert Res, 34(2): 590–596

[24] Munteanu C, Kuemmerle T, Boltiziar T, Butsic V, Gimmi U, Lúboš Halada, Kaim D, Király G, Éva Konkoly-Gyuró , Kozak J, Lieskovský J, Mojses M, Müller D, Ostafin K, Ostapowicz K, Shandra O, Štych P, Walker S, Radeloff V C (2014). Forest and agricultural land change in the Carpathian region—a meta-analysis of long-term patterns and drivers of change. Land Use Policy, 38: 685–697

[25] Ouyang Z Y, Wang X K, Miao H (1999). A primary study on Chinese terrestrial ecosystem service and their ecological-economic values. Acta Ecol Sin, 10(5): 606–613

[26] Pan J H, Shi P J, Zhao R F (2010). Research on optimal allocation model of land use structure based on LP-MCDM-CA model: the case of Tianshui. J Mt Sci, 28(4): 407–414

[27] Potter C S, Randerson J T, Field C B, Matson P A, Vitousek P M, Mooney H A, Klooster S A (1993). Terrestrial ecosystem production: a process model based on global satellite and surface data. Global Biogeochem Cycles, 7(4): 811–841

[28] Qi S H, Jiang M X, Yu X B (2011). Evaluating soil erosion in Jiangxi Province with USLE model and remote sensing technology during 1995‒2005. China Environ Sci, 31(07): 1197–1203

[29] Raudsepp-Hearne C, Peterson G D, Bennett E M (2010). Ecosystem service bundles for analyzing tradeoffs in diverse landscapes. Proc Natl Acad Sci USA, 107(11): 5242–5247

[30] Ren P, Hong B T, Zhou J M (2013). Study on ecological security evaluation and spatial characteristics of cultivated land of agriculture production area in the upper Yangtze River. China Populat, Resourc Environ, 23(12): 65–69

[31] Song W, Deng X Z, Yuan Y W, Wang Z, Li Z H (2015). Impacts of land-use change on valued ecosystem service in rapidly urbanized North China Plain. Ecol Model, 318: 245–253

[32] Su Y, Chen X, Liao J, Zhang H, Wang C, Ye Y, Wang Y (2016). Modeling the optimal ecological security pattern for guiding the urban constructed land expansions. Urban for Urban Green, 19: 35–46

[33] Shi N N, Han Y, Wang Q, Quan Z J, Luo Z L, Ge J S, Han R Y, Xiao N W (2018). Construction and optimization of ecological network for protected areas in Qinghai Province. Chinese J Ecol, 37(6): 1910–1916

[34] Tang J, Mao Z L, Wang C Y, Xu X M, Han W Z (2009). Regional land use structure optimization based on carbon balance: a case study in Tongyu County, Jilin Province. Resourc Sci, 31(1): 130–135

[35] van Vliet J, De Groot H L F, Rietveld P, Verburg P H (2015). Manifestations and underlying drivers of agricultural land use change in europe. Landsc Urban Plan, 133: 24–36

[36] Verhagen W, Van Teeffelen A J A, Verburg P H (2018). Shifting spatial priorities for ecosystem services in europe following land use change. Ecol Indic, 89: 397–410

[37] Wang C, Peng Q, Tang N, Li H Y (2018). Spatio-temporal evolution and the synergy and trade-off relationship of cultivated land multi-function in 2005‒2015: a case of Shapingba District, Chongqing City. Scientia Geographica Sinica, 38(4): 590–599

[38] Wang G, Wang J L, Gong L Y, Su B L, Liu Q B (2013). Spatial-temporal evolution of regional eco-security based on GIS-Markov model—a case study of Ganjingzi District in Dalian, Liaoning Province. Scientia Geographica Sinica, 33(8): 957–964

[39] Wang L P, Jin X B, Du X D, Zhou Y K (2012). Land use scenarios simulation of Foshan City based on gray model and cellular automata model. Transact Chinese Soc Agr Eng, 28(3): 237–242

[40] Wang Y Y, Shen C Z, Jin X B, Bao G Y, Liu J, Zhou Y K (2019). Developing and optimizing ecological networks based on MSPA and MCR model. Ecologic Sci, 38(2): 138–145

[41] Wei W, Shi P J, Zhou J J, Xie B B, Li C H, Lei L (2016). Configuration partition of land use optimization in arid inland river basin based on ecological security pattern. Transact Chinese Soc Agr Eng, 32(18): 9–18

[42] Williams J R, Renard K G, Dyke P T (1983). EPIC: a new method for assessing erosion’s effect on soil productivity. J Soil Water Conserv, 38(5): 381–383

[43] Wu Z, Wang H (2015). Establishment and optimization of green ecological networks in Yangzhou City. Chinese J Ecol, 34(7): 1976–1985

[44] Xiang Y Y, Meng J J (2013). Research on optimization of land use structure in Wuhan urban agglomeration based on ecological benefit. Resourc Environ Yangtze Basin, 22(10): 1297–1304

[45] Xie G D, Zhang C S, Xiao Y, Lu C X (2015). The value of ecosystem services in China. Resourc Sci, 37(9): 1740–1746

[46] Xu J, Xiao Y, Xie G D, Jiang Y (2019). Ecosystem service flow insights into horizontal ecological compensation standards for water resource: a case study in Dongjiang Lake basin, China. China Geogr Sci, 29(2): 214–230

[47] Xu X J, Liu H Y, Lin Z S, Liu J X, Li L H (2017). Scenario analysis of land use change in Jiangsu coast based on CA-Markov model. Res Soil Water Conserv, 24(1): 213–218

[48] Yan C, Ji L (2001). An index of equilibrium of urban land-use structure and information dimension of urban form. Geogr Res, 20(2): 146–152

[49] Yang Q S, Qiao J G, Ai B (2013). Simulation of urban ecological security pattern based on cellular automata: a case of Dongguan City, Guangdong Province of South China. Chinese J Ecol, 24(9): 2599–2607In Chinese)

[50] Ye Y Y, Su Y X, Zhang H O, Liu K, Wu Q T (2014). Ecological resistance surface model and its application in urban expansion simulations. Acta Ecol Sin, 69(4): 485–496

[51] You W B, He D J, Wu L Y, Hong W, Zhan S H, Qin D H, You H M (2011). Temporal-spatial differentiation and its change in the landscape ecological security of Wuyishan Scenery District. Acta Ecol Sin, 31(21): 6317–6327

[52] Yu D, Shi P, Shao H, Zhu W, Pan Y (2009). Modelling net primary productivity of terrestrial ecosystems in east asia based on an improved casa ecosystem model. Int J Remote Sens, 30(18): 4851–4866

[53] Yu F, Li X B, Wang H (2014). Optimization of land use pattern based on eco-security: a case study in the huangfuchuan watershed. Acta Ecol Sin, 34(12): 3198–3210

[54] Yu H Y, Zhang F, Cao L, Wang J, Yang S T (2017). Spatial-temporal pattern of land ecological security at a township scale in the Bortala Mongolian Autonomous Prefecture. Acta Ecol Sin, 37(19): 6355–6369

[55] Yu J, Fang L, Cang D B, Zhu L, Bian Z F (2012). Evaluation of land eco-security in Wanjiang district base on entropy weight and matter element model. Transact Chinese Soc Agr Eng, 28(5): 260–266

[56] Yuan M, Liu Y L (2014). Land use optimization allocation based on multi-agent genetic algorithm. Chinese Soc Agri Eng, 30(1): 191–199

[57] Zeng Y, Qiu X F, Liu C M, Aoda P (2010). Ecological safety evaluation of three gorges reservoir area in Chongqing with the pressure-state-response model. Prog Geogr, 29(9): 1095–1099

[58] Zhang J P, Qiao Q, Liu C L, Wang H H, Pei X (2017). Ecological land use planning for Beijing City based on the minimum cumulative resistance model. Acta Ecol Sin, 37(19): 6313–6321

[59] Zhang J Q, Wu Y J, Ge Y, Wang C H, Kung H (2014). Eco-security assessments of poor areas based on gray correlation model: a case study in Enshi. Geogr Res, 33(8): 1457–1466

[60] Zhang L, Chen Y, Wang S T, Men M X, Xu H (2015). Assessment and early warning of land ecological security in rapidly urbanizing coastal area: a case study of Caofeidian new district, Hebei, China. Chinese J App Ecol, 26(8): 2445–2454 in Chinese)

[61] Zhao D, Li F, Wang R S (2011). Optimization of urban land structure based on ecological green equivalent: a case study in Ningguo City, China. Acta Ecol Sin, 31(20): 6242–6250

[62] Zhao H B, Ma Y J (2014). Spatial-temporal pattern and obstacle factors of cultivated land ecological security in major grain producing areas of northeast China: a case study in Jilin Province. Chinese J App Ecol, 25(2): 515–524 in Chinese)

[63] Zhou K H, Liu Y L, Tan R H, Song Y (2014a). Urban dynamics, landscape ecological security, and policy implications: a case study from the Wuhan area of central China. Cities, 41: 141–153

[64] Zhou R, Wang X J, Su H L, Qian X, Sun B (2014b). Delimitation of urban growth boundary based on ecological security pattern. Urban Planning Forum, (4), 57–63

[65] Zhu X Z, Li X W, Jia K J, Qi F (2014). A study on system dynamics of land comprehensive carrying capacity in Shanghai City. China Land Sci, 2: 90–96

RIGHTS & PERMISSIONS

Higher Education Press