Simulation of Dynamic Urban Expansion under Ecological Constraints Using a Long Short Term Memory Network Model and Cellular Automata

Abstract

Rapid urban expansion has seriously threatened ecological security and the natural environment on a global scale, thus, the simulation of dynamic urban expansion is a hot topic in current research. Existing urban expansion simulation models focus on the mining of spatial neighborhood features among driving factors, however, they ignore the over-fitting, gradient explosion, and vanishing problems caused by the long-term dependence of time series data, which results in limited model accuracy. In this study, we proposed a new dynamic urban expansion simulation model. Considering the long-time dependence issue, long short term memory (LSTM) was employed to automatically extract the transformation rules through memory units and provide the optimal attribute features for cellular automata (CA). This study selected Lanzhou, which is a semi-arid region in Northwest China, as an example to confirm the validity of the model performance using data from 2000 to 2020. The results revealed that the overall accuracy of the model was 91.01%, which was higher than that of the traditional artificial neural network (ANN)-CA and recurrent neural network (RNN)-CA models. The LSTM-CA framework resolved existing problems with the traditional algorithm, while it significantly reduced complexity and improved simulation accuracy. In addition, we predicted urban expansion to 2030 based on natural expansion (NE) and ecological constraint (EC) scenarios, and found that EC was an effective control strategy. This study provides a certain theoretical basis and reference value toward the realization of new urbanization and ecologically sound civil construction, in the context of territorial spatial planning and healthy/sustainable urban development.

1. Introduction

The rapid population growth, social economic development, and the continuous expansion of urban built-up land have led to an array of significant challenges on a global scale, such as the aggravation of climate change, depletion of natural resources, and the degradation of the ecological environment [1,2]. In particular, contradictions between urban expansion and ecological environmental degradation have become increasingly severe [3,4]. The growth rate of urban built-up land in the 21st century is 2.14 times that of the 1990s in China [5]. High-value ecological land, such as cultivated land, forests, and water body, is faced with irreversible risk of being quickly converted to built-up land [6]. Therefore, the government points out that it is necessary to strengthen ecosystem protections by completing the delineation of the three control lines that include ecological protection, and the permanent basic farmland red lines and urban development boundaries [7]. Consequently, the key issues that require urgent resolution encompass how to objectively and comprehensively simulate the spatiotemporal dynamics of urban expansion in the future, scientifically and rationally optimize the spatial land planning system, and improve the quality of the ambient ecological environment [8].

There are two major issues involved in the modeling and prediction of dynamic future urban expansion trends: (1) how to employ the most objective parameterized method to obtain optimally accurate prediction results and (2) whether urban expansion simulations based on variable scenarios lead to substantial differences [9,10]. Various models have been developed and applied to predict the future dynamic changes in urban expansion, such as system dynamics models [11], agent-based models [12], machine learning models [13], cellular automata (CA) models [7], and deep learning models [14]. Different models have their own advantages and disadvantages. System dynamics models possess strong analytical capacities, whereas, it is difficult to address the interactions and feedback relationships between driving factors [10]. Machine learning models typically explore urban expansion processes through the establishment of mathematical relationships between conversion probability and various variables. However, it is difficult to interpret physical mechanisms under urban expansion using these techniques [15]. CA models, which adopt iterative, bottom-up calculations combined with the nearest neighbor structure, boundary constraints, and cell transformation rules, have a strong capacity to simulate the self-organization phenomenon of complex systems. Consequently, they are extensively employed in various simulations to realize the modeling of urban expansion [16]. However, they often rely on other models for the design of transformation rules, such as artificial neural networks (ANNs), support vector machines (SVMs), random forests (RFs), and recursive neural networks (RNNs) [17]. Due to the possible over-fitting of the training data, some scholars through approximation error analysis, it is confirmed that a deep neural network (DNN) can provide better generalization ability than ANN to overcome this problem in the noise environment with deep depth and wide width of hidden layer on top [18,19,20,21]. Most existing studies assume that input and output elements are independent of each other, and focus on the mining of spatial neighborhood features and the improvement of transitional suitability, while ignoring the long-term dependence of time series data and the effects of contextual features on findings [22]. When extracting the historical information of time series data, it is time-consuming to repeatedly conduct back propagation learning; the gradient signal will explode or vanish, which results in the end of the learning process [23,24]. Subsequently, simulation performance is inevitably limited. Time convolutional network is an improved convolutional network in DNN, which can solve the problem of gradient explosion or vanishing. However, it is a unidirectional structure, which requires high integrity of time series data. Additionally, when migrating a model from a problem that requires little memory information to one that requires longer memory, the results are poor [25].

To resolve the above issues, a new CA model that incorporated deep learning technology was proposed to simulate spatiotemporal urban expansion dynamics by exploiting its respective advantages in the extraction of temporal and spatial features. As a time-recursive neural network, long short term memory (LSTM) cannot only effectively process relatively long-time intervals and delay time series, but also accurately predict subsequent development trends. It has been widely utilized the in speech recognition, image classification, data mining, analysis, and prediction domains [26,27,28]. Xue and Xue (2018) [29] applied LSTM to conduct research on traffic flow prediction, and it was demonstrated that LSTM could further improve prediction accuracy based on contextual features. Liu and Sun (2020) [30] proposed a collaborative LSTM-CA computing method to establish a simulation model for the identification of volcanic ash cloud diffusion, and verified that this method can achieve high accuracy (96.1%). For this study, the LSTM-CA coupling model and two traditional models (ANN-CA and RNN-CA) were utilized for the same test data and parameter settings. Lanzhou, which is a semi-arid region in Northwest China was selected as an example, and existing data from 2000 to 2020 were used to compare and verify the model.

It is difficult to accurately predict future urban expansion according to historical trends alone, as newly released policies can easily change urban development patterns. Scenario-based modeling, which combines historical trends and related policies, is attracting academic and governmental attention due to its potential for evaluating policies that support sustainable development goals [31,32]. Thapa and Murayama (2012) [31] designed environment-protecting and resource-saving scenarios from related policies and planning maps to identify future urban growth allocation in Kathmandu Valley, Nepal. In the context of urban expansion, the Chinese government has put forward an “Ecological Protection Red Line” (EPRL) policy, which demarcates areas where construction is restricted within cities, such as ecological and natural resources, and reverses the trend of urban expansion from the perspective of quantifying the suitability of land use [22,33]. Subsequently, many scholars began to simulate urban expansion based on the scenarios of ecological protection and ecological constraint (EC).

As relates to the establishment of EC scenarios, ecological suitability partitioning through the determination of ecological sources (ESs) and resistance surfaces (RSs) is currently considered to be an effective way to ensure ecological security [34]. Peng et al. (2018) [35] identified ES in conjunction with the risk of ecological degradation and the importance of ecological functions, and took impervious surface as an important factor for the correction of RS. However, the ES and RS evaluation systems in existing studies remain inadequate. The resulting low EC expansion mode on the space form showing a continuous open spread outside a city and town may not split the connectivity of the natural ecology and agricultural landscape corridor. Additionally, it might impact the integrity of the regional ecological system function structure, and further weaken the urban hydrology adjustment ability, while increasing the thermal environment and heat island effects [36]. Habitat quality and landscape connectivity are essential toward achieving regional ecological security, and play a key role in the determination of ES [6]. The ambient environment, socioeconomics, and ecological factors also play vital roles in the maintenance of ecological environmental protection and can be employed as significant parameters for the construction of RS [37]. Therefore, improving the ES and RS evaluation system for ecological suitability partitioning and the extraction of constraint factors can provide a new way to simulate urban expansion under EC.

Urban expansion in the economically developed regions of Eastern China has received more attention and research; however, less research has been done in the regions of Northwest China. Akin to many other cities worldwide, conflicts between the rapid expansion of urban built-up and the chronic shortage of ecological resources is vitally prominent in Lanzhou. Lanzhou is a typical semi-arid valley city with significant ecological vulnerability due to meteorological factors of low precipitation and high evaporation [38]. Therefore, this study aimed to: (1) clarify the spatiotemporal evolution of urban expansion in Lanzhou from 2000 to 2020; (2) determine a comprehensive and efficacious evaluation system and method for ecological suitability partitioning to improve the EC scenario; and (3) develop a comprehensive and easy-to-use LSTM-CA framework for the prediction of dynamic urban expansion to 2030 under different scenarios. This study will be useful toward understanding the response relationships between ecosystems and urban expansion in semi-arid valley basins, while providing significant guidance for ecosystem restoration.

The remainder of this paper is structured as follows: Section 2 outlines the study area and data sources; Section 3 introduces the model basics and methodology used in this research; Section 4 presents and analyzes the prediction results based on different models and scenarios; finally, Section 5 presents the discussion, and Section 6 ends with the conclusion.

2. Study Area and Data Sources

2.1. Study Area

Lanzhou (Figure 1) (103°40′ E, 36°03′ N) is located in the Northwest of China. It is one of the typical river valley basin-type cities in semi-arid areas [39]. As of July 2020, the city covers an area of approximate 1.31 × 10⁴ km² and has a population of over 4.13 × 10⁶. The terrain of the study area is high in the southwest and low in the northeast, with an average altitude of 1520 m. The Yellow River flows from the southwest to northeast, forming a beaded valley with alternate canyons and basins, which runs 35 km from east to west, and 2–8 km from north to south [40]. The landforms are complex and diverse, with a staggered distribution of mountains, plateaus, plains, and river valleys. The average annual temperature, precipitation, relative humidity, and evaporation are about 10.3 °C, 327 mm, 56%, and 1448 mm respectively [41]. Lanzhou primarily includes important ecological spaces, such as critical areas for biodiversity and water conservation, special protection areas, vital national forest parks, natural desert grasslands, natural ecosystems, and habitats for wild species. As one of the most dynamic economic development areas in Northwestern China, and an important node city of the “Belt And Road” initiative, Lanzhou currently faces the serious degradation of ecosystem functionality and the overall environment.

2.2. Data Sources

All of the data sets used in this study, and their formats and sources are described in

Table A1. Data Format and Source.

Dataset	Format	Resolution	Time	Data Source
Land use/cover	Vector (Polygon)	30 m	2000, 2010, 2020	Global Geographic Information Product Platform (http://www.globallandcover.com) (accessed on 15 November 2020)
Nighttime light	Vector (Polygon)	130 m	2018	Application Platform of High-resolution Earth Observation System for Hubei Province, China (http://59.175.109.173:8888) (accessed on 16 November 2020)
Faults	Vector (Polyline)			“Gansu Hydrogeological Map” (Gansu Bureau of Geology and Mineral Hydrogeology engineering Geological Exploration Institute) (http://www.gssgy.com) (accessed on 15 November 2019)
DEM	Raster	30 m	2015	USGS Earth Explorer (https://earthexplorer.usgs.gov) (accessed on 26 December 2019)
Slope	Raster	30 m	2015	USGS Earth Explorer (https://earthexplorer.usgs.gov) (accessed on 26 December 2019)
NDVI	Raster	30 m	2015	MODIS (https://modis.gsfc.nasa.gov) (accessed on 20 December 2020)
Settlements	Vector (Point)		2019	Crawled POI data Official website of Gold Maps (https://lbs.amap.com) (accessed on 14 December 2020)
Roads	Vector (Polygon)		2020	Open Street Map (https://webmap.cn) (accessed on 13 December 2020)
Temperature	Raster		2015	China Meteorological Data Service Centre (http://data.cma.cn) (accessed on 14 November 2020)
Precipitation	Raster		2015	China Meteorological Data Service Centre (http://data.cma.cn) (accessed on 2 November 2019)
Natural protected areas	Vector (Polygon)		2015	Portal website of Gansu Forestry and Grass Bureau and its administrative departments (http://lycy.gansu.gov.cn) (accessed on 25 December 2019)
Population density	Excel		2017	Data Center for Resources and Environmental Sciencesof the Chinese Academy of Sciences (http://www.resdc.cn) (accessed on 28 December 2019)
GDP	Excel		2017	Data Center for Resources and Environmental Sciencesof the Chinese Academy of Sciences (http://www.resdc.cn) (accessed on 28 December 2019)
Vegetation type	Raster	1 km	2001	Data Center for Resources and Environmental Sciencesof the Chinese Academy of Sciences (http://www.resdc.cn) (accessed on 28 December 2019)
Surface evapotranspiration	Raster	1 km	2014	MOD16 (https://modis.gsfc.nasa.gov) (accessed on 20 December 2019)
Airports	Vector (Point)		2019	Crawled POI data Official website of Gold Maps (https://lbs.amap.com) (accessed on 14 December 2020)
NPP	Raster		2010	Data Center for Resources and Environmental Sciencesof the Chinese Academy of Sciences (http://www.resdc.cn) (accessed on 28 December 2019)

in Appendix A. Time series land use data from 2000, 2010, and 2020 were derived from the Global Geographic Information Product Platform (Figure 2), with a resolution of 30 m, including eight land use types (cultivated land, forest land, grass land, shrubland, wetland, water, built-up land, and unused land). Urban expansion is affected by traffic, location, population density, public facilities, government planning policies, and other factors [42]. On this basis, the 14 data points were identified as the driving factors of the urban expansion simulation (Figure 3). Natural environmental factors included: DEM, slope, distance to fault, and NDVI; climate factors included: precipitation and temperature; social and economic factors included: nighttime lights, GDP, population, NPP, and nature protected areas; locational factors included: distance to roads, distance to the settlements, and distance to airports. All data were resampled to 30 m, and normalized to [0, 1] to reduce computation costs and accelerate the simulation process. Mean annual temperature and precipitation data were obtained through spatial interpolation of meteorological station data in Gansu and all surrounding counties. Firstly, the meteorological point and altitude, latitude, and longitude data were conducted and multiple regression, and interpolation was generated by adopting the inverse distance weighting method. Then, the data of 130 meteorological stations (70%) were used as the training set, and the remaining 56 (30%) were used as the verification set, with a total of 186 meteorological stations. The overall accuracy of the interpolation was 89%.

3. Methods

The general structure of the proposed framework is illustrated in Figure 4 and involved two main technical steps. (1) LSTM-CA prediction model (a) and (2) simulation prediction based on natural expansion (NE) and EC scenarios. Among them, the NE scenario was an unconstrained setting. The EC scenario was set as follows: ES were identified by constructing a comprehensive evaluation system of habitat quality, ecosystem services, and landscape connectivity; RS evaluation system was constructed from the natural environment, socioeconomics, and ecology; minimum cumulative resistance (MCR) model was implemented to partition ecological suitability, and the partition results were used as EC factors (b).

3.1. LSTM-CA Model

During the process of building the integration model, CA was used to control the spatial dimension features as the upper framework, and the input data was required to have attribute information corresponding to the neighborhood in CA. As the underlying structure, LSTM was employed to automatically extract the transformation rules in CA and look for the rules on time dimension during the course of processing serialized data [30]. Consequently, the input data format of the LSTM-CA integration model was a three-dimensional tensor (sample, time step, and attribute), which was controlled by the neighborhood in CA. The time step and attribute corresponded to the features of time and space, respectively. Since the training process of the integrated model was automatic, it only needed to construct different inputs in terms of the spatiotemporal features of the original data, whereafter the corresponding results could be directly obtained.

On the basis of RNN, LSTM introduced a “processor”, namely a new memory unit structure for determining the usefulness of data. A memory unit was primarily comprised of four elements (input gate, self-circulating connecting neuron, forgetting gate, and output gate) [43]. The most important component of the LSTM model was the state of cells. We employed an activation function to input the memory cell information of the network into the LSTM for calculation. The output results that conformed to the rules of the algorithm continued to be the input of the subsequent layer, while those that did not conform were dispensed with via the forgetting gate. During testing, the learned local model assigned weights to update the relationship between the input factors and predictions for each step of the LSTM.

The calculation procedure for the LSTM was as follows (Figure 5):

Step 1: Connect the hidden state of the previous time step with the input of the current time step. The hidden state controls the information in the memory cell through the output gate. Each time step contains a hidden state, usually using the previous hidden state to update the memory unit at this moment.

Step 2: The connection result is transferred to the forgetting gate, and the data is selectively deleted following remapping by the sigmoid activation function. By processing the previous hidden state and the current input, the range of values after the output is [0, 1], and those close to 0 are deleted, while those close to 1 are retained.

Step 3: The connection result is transferred to the input gate, where data from candidate memory cells will be added to the memory cells after remapping by the activation function;

Step 4: The connection result is delivered to the tanh activation function, and the candidate memory cells will be output. The calculation method of candidate memory cell is similar to the three gates, the output results are all [0, 1], and those close to 1 are output.

Step 5: Using the outputs from Step 2 to Step 4 and combining with the memory cells of the previous time step, new memory cells are calculated and generated;

Step 6: The output gate determines the next hidden state, which contains valid information from the previous time step. Similar to the first two gates, the previous time step hidden state and the current input (Step 1) are passed to the output gate, which is mapped by the tanh activation function to determine which data in the memory cell will be passed to the hidden state.

The updating calculation formula of the t layer of LSTM was as follows:

$$\left\{\begin{array}{l}{i}_t=\sigma (W_i\cdot \left[h_{t-1},x_t\right]+b_i)\\ f_t=\sigma (W_f\cdot \left[h_{t-1},x_t\right]+b_f)\\ o_t=\sigma (W_o\cdot \left[h_{t-1},x_t\right]+b_o)\\ {\tilde{c}}_t=\mathit{Tanh}(W_c\cdot \left[h_{t-1},x_t\right]+b_c)\\ c_t=f_t{c}_{t-1}+i_t\odot {\tilde{c}}_t\\ h_t=o_t\odot \mathit{Tanh}(c_t)\end{array}\right.$$

(1)

where, $i_t,f_t,o_t$ represent input gate, forgetting gate, and output gate, respectively; ${\tilde{c}}_t$ represents candidate memory cells; $c_{t-1}$ and $c_t$ represent the memory unit of the previous time step and the current time step; $h_{t-1}$ and $h_t$ represent the hidden state of the previous time step and current time step;$x_t$ represents the cell state of the current timestep input; $\sigma$ is the activation function (sigmoid, etc.); $W_i$ represents the weight matrix of the input gate; and $b_i$ represents the offset of an input gate. Unlike traditional RNN, the transfer from the state $c_{t-1}$ of the previous unit of memory to the current state $c_t$ does not necessarily completely depend on the state calculated by activation function $\sigma$, also to be controlled by $i_t$ and $f_t$. This makes it easier for the entire network to learn about long-term dependencies between sequences.

The model was divided into training and prediction, where given the input and output during the training phase, the calculation and storage of the gradient in the LSTM were propagated back through time to determine the dependency between the model variables and parameters. The land use data and driving factors were normalized for preprocessing and then iterated, and the loss value can be quickly reduced to a certain value. At this time, the iterative process basically converges and a training model with high accuracy is obtained. During the prediction phase, the set parameters can be applied to a given input and the output can be obtained. Based on the urban expansion simulation data features, a typical four-layer LSTM neural network structure was constructed, where the first two layers were LSTM layers, and the last two layers were fully connected layers. The first LSTM layer contained 80 memories, which output each time step, while the second layer contained 100 memories, which only accepted the output of the last time step. To prevent over-fitting, in the full connection the rejection probability of the first layer was set at 20%. Softmax was activated at the second layer and the rejection probability was set at 50%. The learning rate was set to 0.001, the loss function was the classical mean square error, the activation function was the tahn function, and the optimizer selects the commonly used adaptive moment estimation Adam algorithm, as the optimizer has the advantages of lower resource requirements, a fast convergence speed, and high accuracy.

3.2. Partitioning of Ecological Suitability

3.2.1. Identification of Ecological Sources

The weights of habitat quality, ecosystem services, and landscape connectivity obtained by the Delphi method were 0.25, 0.50, and 0.25, respectively. The ES result was obtained by the superposition of protected nature areas, and all data were normalized to [0, 1].

Habitat quality adopted the INVEST model [44], and in terms of the sensitivity of the habitat to threat sources, the degree of degradation of the habitat under different land use types can be identified so as to obtain the habitat quality of the study area (Equation (2)) [45].

$$Q_{xj}=H_j\left[1-\left(\frac{V_{xj}^z}{V_{xj}^z+k^z}\right)\right]$$

(2)

where $Q_{xj}$ is the habitat suitability; $k$ is the subsaturation constant, generally half of the maximum value of $V_{xj}$; and $z$ is a constant. In this study, the model parameters were determined by referring to previous results and collinearity test results (

Table 1. Habitat quality assessment parameters.

Threat Source	Weight	Sensitivity						Maximum Influence Distance (km)	Attenuation Type
		Cultivated	Forest	Grass	Shrubland	Wetland	Water
Build-up	1	0.4	0.75	0.45	0.65	0.8	0.8	10	Exponential
Human activity intensity	1	0.75	0.75	0.5	0.65	0.8	0.8	10	Exponential
Highway	0.8	0.8	0.65	0.3	0.55	0.7	0.65	3	linear
Railway	1	0.5	0.55	0.2	0.45	0.6	0.55	5	linear
Main road	1	0.7	0.58	0.6	0.6	0.6	0.55	5	linear

) [6,46]. The parameter settings in the model are presented in Table 1. The nighttime light index is often used to characterize the intensity of human activities.

Water yield refers to the process and capacity of water conservation in the ecosystem [47]. Soil retention is based on the potential and actual amount of soil erosion caused by rainfall, soil, slope, vegetation, and land management [48]. Biodiversity is an important component of ecosystem services [49,50]. For this study, water yield, soil conservation, and biodiversity were employed to characterize ecosystem services.

Landscape connectivity is defined as the continuity of corridors, networks, or spatial matrices that quantitatively characterize the relationships between ecosystems [51]. Morphological spatial pattern analysis was adopted to identify and extract the core area, and the probability of connectivity (PC) index (Equation (3)) was selected to evaluate the landscape connectivity of the core area [52]. Repeated experiments show that 500 m is a more suitable fixed distance for connectivity analysis using the PC index.

$$PC=\frac{{\displaystyle \sum _{i=1}^n{\displaystyle \sum _{j=1}^n{a}_i\times a_j\times p_{ij}^{\ast }}}}{A_L^2}$$

(3)

where $A_L$ is the total area of the landscape; $a_i$ and $a_j$ are the areas of ecological patches $i$ and $j$, respectively; $p_{ij}$ is the maximum connectivity value of all paths between patches $i$ and $j$; and $n$ is the number of patches.

3.2.2. Identification of Resistance Surfaces

During the process of urban expansion, location conditions, traffic networks, and other positive factors of economic and social development, terrain conditions, hydrological environments, and other adverse constraints will increase the economic costs and ecological risks of grid unit development [53]. This research integrated the characteristics and data availability of the study area, DEM, slope as a natural environment factor, the net primary productivity (NPP), nighttime lights index, and road network as social and economic factors, ecosystem service, and the surface of the ecological vulnerability as ecological factors, to establish a cumulative RS evaluation system, and the weight of each factor was obtained through an analytic hierarchy process.

3.2.3. Minimum Cumulative Resistance Model

MCR can effectively simulate the hindrance of landscape in space movement to better express the interactive relationships of landscape ecological processes. Considering the three factors of source, distance, and resistance, ecological protection land can be identified to maximize the benefits of built-up and ecological land (Equation (4)) [53]. The EP scenario aims to protect important ecological areas with specific location and ecosystem functions to maintain regional ecological security. On the basis of existing data (Table A1 in Appendix A), MCR model was implemented to partition ecological suitability, and the normalization of partition results were incorporated in the algorithm as constraint factors.

$$MCR=f\min {\displaystyle \sum _{j=n}^{i=m}{D}_{ij}\times R_i}$$

(4)

where, $MCR$ is the value of minimum cumulative resistance, and $D_{ij}$ is the distance of ecological land from source unit $y$ to landscape unit $i$. $R_i$ represents the resistance coefficient of landscape unit $i$ to the movement process; $f$ is the relationship function between cumulative resistance and ecological suitability.

4. Results

4.1. Spatiotemporal Pattern Evolution of Urban Expansion in Lanzhou

In terms of area statistics, Lanzhou has experienced considerable urban expansion over the last 20 years (Figure 2). Built-up land increased from 229.27 km² in 2000 to 757.76 km² in 2020, which translates to a 3.31 fold increase in urbanized area at an annual growth rate of 23.05%. However, the change trend varied for different periods. From 2000 to 2010, the study area was at a stage of slow expansion, whereas from 2010 to 2020, it experienced a stage of rapid expansion with an annual growth rate of 28.23%, where the total built-up area increased by 559.54 km². From the perspective of spatial distribution, each portion of the study area shrank slightly from 2000 to 2010, and the changing patches were relatively broken. From 2010 to 2020, not only was the built-up land in the main urban area extended to the surrounding areas, but the Lanzhou New District at the junction of Yongdeng and Gaolan also became the main space for urban expansion. Various counties urban expansion of number: Gaolan (35.18%) > Yong deng (28.74%) > Yuzhong (9.07%) > Chengguan (6.15%) > Qilihe (5.65%) > Xigu (5.32%) > Honggu (5.16%) > Anning (4.73%).

Furthermore, other land use types were analyzed. From 2000 to 2010, the forest land area decreased from 156.55 to 122.69 km², at an annual reduction rate of 2.16%. From 2010 to 2020, the forest land increased in the northwest of Yongdeng and the southeast of Yuzhong, which were approved as nationally protected nature areas in 2001 and 1988, respectively. Meanwhile, the cultivated lands, grass lands, and shrublands changed little during 2000–2010, but revealed a decreasing trend year by year from 2010 to 2020, with a net loss of 460.65 km², 191.49 km², and 163.35 km², at annual reduction rates of 0.83%, 0.27%, and 9.34%, respectively. Wetlands decreased year by year, and the annual reduction rate was higher from 2000 to 2010 than from 2010 to 2020, whereas changes in water and unused land were not obvious (

Table 2. Area and proportion of different land use types from 2000 to 2020.

Land Use Type	2000		2010		2020
	Area (km2)	Percent (%)	Area (km2)	Percent (%)	Area (km2)	Percent (%)
Cultivated	5224.50	39.92	5258.99	40.19	4820.34	36.84
Forest	156.55	1.2	122.69	0.94	361.69	2.76
Grass	7187.99	54.93	7223.70	55.2	7032.21	53.74
Shrubland	175.20	1.34	174.81	1.34	11.46	0.09
Wetland	1.20	0.01	0.60	0	0.52	0
Water	42.89	0.33	38.58	0.29	43.27	0.33
Build-up	229.27	1.75	198.22	1.51	757.76	5.79
Unused	68.70	0.52	68.71	0.53	59.06	0.45

). Comparatively, it was demonstrated that cultivated land was the largest land resource encroached upon via urban expansion in Lanzhou, as it was the land type most intimately associated with anthropogenic activities. The net loss of forest lands was much smaller than that of cultivated land, as the forested portions of the study area accounted for less land, were further away from built-up land and at a higher altitude, possessed stronger ecosystem services, and were within nationally protected nature areas.

4.2. Verification of Urban Expansion Simulation Model

We adopted historical land use data and 14 driving factors to estimate the regression coefficient and correct the model parameters. The proposed model was encoded in Python to obtain the simulated land use map in 2020. The accuracy evaluation criteria of the OA, Kappa coefficient, PA, UA, and FOM were used to compare the simulation results of LSTM-CA, ANN-CA, and RNN-CA models in 2020. To obtain effective comparative results, the parameter settings of the three models in the experiment were identical, with the accuracy evaluation results shown in

Table 3. Accuracy compared with traditional models.

Model	OA	Kappa	FOM	PA	UA	LOSST	LOSSV
RNN-CA	0.8326	0.7827	0.2605	0.3846	0.5063	0.0265	0.0261
ANN-CA	0.8837	0.8268	0.3899	0.5358	0.6582	0.0198	0.0186
LSTM-CA	0.9101	0.8662	0.4439	0.6122	0.7325	0.0153	0.0142

. The OA of the RNN-CA was 88.37%, and its simulation ability was the poorest, whereas the ANN-CA was only a little higher than RNN-CA (83.26%), both of which were less than 91.01% of the LSTM-CA. The accuracy of the other four LSTM-CA items were also better than the two traditional models, which indicated that long-term dependence was a key factor that could not be ignored in the study of urban expansion. The simulation performance could be improved by extracting conversion rules from memory units, which verified the effectiveness and superiority of this model.

Figure 6 illustrates the simulated and actual urban expansion in 2020 based on the three models above. The simulation results of LSTM-CA were the most consistent with the actual situation, followed by ANN-CA, whereas the difference of RNN-CA was obvious. We also recorded the LOSS value curves of each model (Figure 7), LOSS value of the test set (LOSS_T), and LOSS value of the verification set (LOSS_V). The lower the LOSS value, the higher the precision of model training. From the perspective of the LOSS curve and simulation map, the simulation map of these three models was basically consistent with the actual distribution of land use types. However, there were slight difference patterns. The identification of the built-up areas simulated by RNN-CA was quite different from the actual map. The prediction results of land use types with small built-up and forest lands are less than the actual map, on the contrary, the results of land use types with large areas of grass and cultivated land are more than the actual map. The LOSS value was the highest when it converged, and the LOSS curve was opposite to the trend of the other two. With the increased number of iterations, LOSS presented an increasing trend, which demonstrated that the prediction results had over-fitting and gradient rise phenomena. The proportion of built-up areas of ANN-CA was more dispersed than that of the actual land use map, and the LOSS value was also higher than that of the LSTM-CA. By comparing the LOSS curve, this model was weaker than the LSTM-CA in simulating the dynamic change of spatiotemporal urban expansion. The LSTM-CA simulation effect was the best, and the LOSS value fitted by LSTM-CA was the smallest, which indicated that the model could not only improve accuracy, but also avoided the gradient explosion or vanishing caused by the gradient rise.

4.3. Partitioning of Ecological Suitability

The high quality habitats were distributed evenly in the areas with high forest and grass land coverage in Yongdeng, Gaolan, and Yuzhong. The overall level of habitat quality in these areas was higher; thus, the risk of ecological degradation was higher. The threat degree of water was high, and the corresponding habitat quality was low (Figure 8a). The higher ecosystem service areas were primarily distributed north of the Yellow River, with an overall mean value of about 0.5 (Figure 8b). Areas with high landscape connectivity were distributed mainly in the Liancheng National Nature Reserve of Yongdeng County and the Xinglong Mountain National Nature Protected Areas of Yuzhong County. With low and relatively low levels accounting for more than 50%, this indicated that the landscape connectivity of the ecological land was poor (Figure 8c). The results of habitat quality, ecosystem services, and landscape connectivity were superimposed to obtain the ES, which accounted for 27.46% of the total study area (Figure 8d). These were mainly distributed across the west of Yongdeng, Honggu, and National Nature Protected Areas. From the perspective of land use type, the ES consisted primarily of forest, grass, and cultivated land. These areas not only had high vegetation cover and ecosystem services such as water yields, soil conservation, and biodiversity conservation, but were also are located at high elevations.

The expansion resistance of ES was mainly distributed in the areas with strong economic capacities, close to urban boundaries and high anthropogenic disturbances, and the corresponding possibility of urban built-up land expansion was high (Figure 9a). Conversely, the expansion resistance of urban built-up land was primarily distributed in regions with complex terrains far from road networks, with good ecological environments, and strong ecosystem services, where the corresponding ES was likely to expand (Figure 9b). The MCR model was implemented to calculate and generate the minimum cumulative resistance surface, and the partition results of ecological suitability were obtained (Figure 9c). Suitable land for urban built-up has the common characteristics of a gentle terrain and contiguous agglomeration, which is far from ES. Urban expansion in these regions poses a negligible threat to ecological security and conveys a low degree of damage to the ecological environment. Land that is suitable for ecological protection is divided into key protected areas and general protected areas, where the key protected areas are mainly forests, high-cover grasslands, and water resource areas, which account for 32.70% of the total area. As core areas for the preservation of regional ecological security, key protected areas maintain ecological integrity, and are not suitable for urban construction; thus, ecological protections need to be strengthened. General protected areas are primarily distributed around the ES and serve as ecological buffers. The findings of this partitioning identify EC for urban expansion simulations and predictions, from the perspective of coordinating ecological protection and urban construction, while providing a scientific basis for sustainable urban development.

4.4. Simulation and Prediction of Dynamic Urban Expansion to 2030

The strategy of ecological control and sustainable development is a policy of the Chinese government to protect natural ecosystems. In this strategy, land within the EPRL is not allowed to be developed into urban land. Therefore, we designed two NE and EC scenarios to simulate the urban expansion of Lanzhou in 2030. Figure 10 illustrates the prediction results of LSTM-CA under these two scenarios. According to the spatiotemporal evolution of urban expansion from 2000 to 2020, cultivated land was the largest land resource encroached upon by urban expansion in Lanzhou. Consequently, this paper mainly analyzed the spatial changes in built-up land representing urban expansion, and cultivated land occupied by urban expansion. Built-up land was an increasing trend in the NE scenario, which was more obvious than the EC scenario. Figure 10b3 revealed that the area of built-up land expanded more under the NE scenario than the EC scenario, with a total area of 19.72 km², which was primarily distributed in the suitable built-up areas of Yongdeng, Gaolan, and Yuzhong. This indicated that without policy interventions, the three counties of Lanzhou would experience rapid urbanization over the next decade. We also found that the expansion rate of built-up land in both scenarios was lower than that in 2010–2020, which may have been mainly due to the shortage of land in the study area and the lack of developable land in the central urban area in recent years. Moreover, the Lanzhou New District has not been properly planned since its establishment in 2012, and its development has been relatively slow. The spatial variation trend of cultivated land in the two scenarios was completely opposite to that of built-up land, where the occupation of cultivated land in EC scenario was inhibited. Figure 10b5 illustrates that the area was occupied more by cultivated land in the NE scenario than in the EC scenario with a total of 6.49 km², which was evenly distributed across all districts and counties.

Furthermore, it was observed that the unused land showed a decreasing trend to different degrees under the two (NE and EC) scenarios, with 5.73 km² and 21.55 km² being replaced by other types of ecological land, respectively. Forest land and water increased slightly under the EC scenario, which indicated that ecological protection should be paid more attention to in the study area. Further, we extracted new built-up land under the two scenarios from 2020 to 2030, and superimposed them on the ecological suitability partitioning result map. It was found that the new built-up land under the EC scenario was located outside the key ecological protected areas, and mainly expanded mainly to the boundary of the central urban area and Lanzhou New District. These results showed that the urban expansion pattern based on EC had a positive impact on regional ecological environment protection, and was conducive to the decision-making of territorial spatial planning.

5. Discussion

5.1. LSTM-CA: A Reliable Coupled Model for Urban Expansion Simulation

The extraction of transformation rules is an important step in dynamic urban expansion simulation [54]. Researchers have employed a variety of methods to extract transition rules, including ANN, SVM, RF, and RNN [14]. These traditional coupled CA models are essentially based on the Markov chain, and assume that the state of a given unit in the next step is related only to the state of the previous step. Urban expansion is a dynamic geographical process with obvious long-term time dependence, where the simulation accuracy will be significantly reduced if the long-term time dependence is ignored [55]. Considering long-term dependence, LSTM-CA can capture spatial and temporal neighborhood features by automatically extracting transformation rules from memory cells. The simulation results revealed that the LSTM-CA model could significantly improve overall accuracy, and was better than the traditional ANN-CA and RNN-CA models, which further confirmed the importance of modeling transformation rules in urban expansion.

5.2. Necessity of Ecological Constraints on Urban Expansion

Ecological security protection is a critical link that prevents urban development from degrading ecological environments. Urban expansion may lead to changes in the composition and configuration of ecosystem functions for various land use types. Previous studies have explored the effects of urban expansion on ecosystems. However, EC is more important than post-destruction management [56,57]. The study area was located in a northwest semi-arid region comprised mainly of grass lands, which was different from the ecological service system of the coastal areas; making it difficult to set EC scenarios. Consequently, as much as was possible, this study collected the latest and most comprehensive data using the most applicable methods. From the perspective of protecting the regional ecological environment, an optimal ES and RS evaluation system was developed, and the ecological suitability was partitioned using MCR technology. Taking the partitioning results as the EC condition, the urban expansion simulation results were obtained.

By comparing the urban expansion of Lanzhou to 2030 under the two scenarios, it was indicated that under the NE scenario, the expansion of urban built-up land area was significantly reduced. The encroachment of cultivated land was also effectively controlled, and the forest land area and water showed an increasing trend. This signified that it was essential to consider habitat quality, ecosystem services, and landscape connectivity to obtain regional ES. Meanwhile, the natural environment, socioeconomic, and ecological factors played a vital role in the maintenance of the ecological environment. Finally, based on the results of ecological suitability partitioning, the following suggestions were made for the territorial spatial planning and healthy urban development of Lanzhou. Firstly, based on the important ecological protection areas, the boundary of EPRL should be quickly demarcated. Combined with the simulation results of urban expansion, it can serve as a reference for the demarcation of permanent basic farmland red line and urban development boundary, so as to avoid the encroachment of a large amount of cultivated land. Secondly, to coordinate the economic management and ecological protection of Lanzhou New District. To understand the influencing mechanisms of different land uses on ecosystem services at different stages of economic development, realize the integrated development of the Lanzhou-Baiyin metropolitan area driven by economic growth, and conduct land reclamation to prevent ecosystem degradation. Thirdly, certain forest systems with relatively low ecosystem services should carry out water yield, soil erosion, and biodiversity restoration, improve the overall landscape connectivity of the region, and increase the urban green space.

5.3. Limitations and Future Work Prospects

Although the LSTM-CA model proposed in this study achieved optimal results in the simulation of urban expansion, there remain some deficiencies due to limitations in data collection and research methods.

Firstly, there were significant neighborhood externalities between different driving factors and adjacent plots. The integrated model also considered both long-term dependence and neighborhood externalities to further improve accuracy. Consequently, future technological breakthroughs can be made by collecting additional driving factors, extracting neighborhood features with the assistance of a three-dimensional convolutional neural network (CNN), solving long-term dependency problems with LSTM and other methods, and obtaining an optimal integration framework by comparing CNN-LSTM and LSTM-CNN models in examples to improve simulation performance.

Secondly, toward the establishment of the partitioning model of ecological suitability, the role of carbon storage in the evaluation of ecosystem services on climate and environment was ignored. It has been confirmed by many scholars that global climate change has a significant impact on the provision of ecosystem services. Meanwhile, carbon storage is also affected by the cumulative effects of extreme climate and land use changes [10,58]. Further study will be required to determine whether there are significant disparities in land use and carbon storage in different regions.

Ultimately, the expansion of built-up land will lead to a decrease in ecosystem service value, which is not fully utilized in studying the optimization of ecological security patterns from the perspective of economic development. Ecosystem service value (ESV) is a key comprehensive index that represents ecological change and an effective means to define the ecological level of urbanization [33]. There are few studies on the quantitative analysis of ESV loss and its serious consequences based on ecological protection policies. The investigation of different types of ESV may assist us with understanding whether sustainable urbanization policies might effectively improve localized ecological conditions.

6. Conclusions

This study undertook to analyze the spatiotemporal evolution of urban expansion in Lanzhou over the last 20 years, and simulated future urban expansion scenarios under different development modes using LSTM-CA models. As an important node city of the “Belt and Road” initiative, Lanzhou has experienced extreme urban expansion, which caused serious ecological degradation. We constructed a complete ES and RS evaluation system using MCR technology to partition ecological suitability. The findings of partitioning were employed as EC to simulate urban expansion, mitigate the threat of urbanization to the ecological environment, and assist with spatial territorial planning. Fast-developing areas need to obtain longer time dependence. Therefore, this study combined LSTM and CA to establish a simulation model of future dynamic urban expansion for semi-arid regions from the perspective of time and space, which can improve scientific prediction results, while providing references for urban planning.

The contributions of this study were as follows:

The LSTM-CA model cannot only effectively deal with relatively long-term intervals and delayed time series, but also accurately predict subsequent developmental trends in the time series, with strong robustness. In view of the nonlinear and highly complex features of urban expansion, the importance of EC and the context features of time series data were fully considered in the LSTM network model. Subsequently, the urban expansion prediction model coupled with LSTM-CA was constructed. The model automatically extracted transformation rules, which significantly reduced the complexity of modeling while improving simulation accuracy.

It was demonstrated that ecological constraints were effective controls for limiting urban expansion in Lanzhou. The ES was established through the consideration of preventing ecological degradation and maintaining the sustainability of ecosystem services and the integrity of the landscape. The RS was constructed by integrating the natural environment, socioeconomic, and ecological factors. The key ecological protected area within the study site covered 4279.00 km², which accounted for 32.70% of the total area and provided a reliable basis for the delineation of the EPRL. By comparing the urban expansion projections under the two scenarios to 2030, it was revealed that the growth of built-up land under the EC scenario decelerated, and all the newly built-up land, according to the ecological suitability partitioning results, was located outside the key protected areas. Therefore, it will be necessary for the government to formulate and implement stringent ecological protection policies.