Simulation of Dynamic Urban Growth Boundary Combining Urban Vitality and Ecological Networks: A Case Study in Chengdu Metropolitan Area

Abstract

The Chengdu Metropolitan Area, located on the eastern edge of the world’s highest plateau, has experienced a period of integrating urban and rural area development for decades. With rapid urbanization and population growth, the vulnerability and security of the ecological environment have become critical aspects to consider in sustainability. Moreover, the presence of different levels of vitality in the study area has a crucial impact on land-use change. Hence, we propose a growth boundary study based on the theory of urban vitality and ecological networks. We focus on identifying the inefficient urban land and urban development potential land, explore their expansion probabilities to conduct spatial simulations for the next 15 years, and combine the ecological network to form a reasonable spatial pattern. Results showed that the proposed model could simulate the urban growth state more accurately within a certain space scale and integrate different limits and influences to simulate different growth strategies under multiple planning periods. Thus, the proposed model can be an effective decision support tool for the government.

1. Introduction

Urban sprawl is called inefficient because it generates low-density development that is “sprawled” over the landscape [1]. In 2016, the United Nations Conference on Housing and Sustainable Urban Development (Habitat III) was held in Quito, the capital of Ecuador. The “New Urban Agenda” adopted at the conference stated that “By 2050, the world’s urban population is expected to nearly double, making urbanization one of the most transformative trends in the 21st century” and that, “We encourage spatial development strategies to consider the need for urban expansion… to prevent over expansion and marginalization of cites”. Over the last three decades, China has experienced rapid urban development and large-scale urban–rural migration. About 75% of China’s population is expected to live in cities within 20 years, resulting in a stronger demand for transport, energy, water, and other basic life necessities, which will be a huge challenge for resource and urban planners [2].

Urban Growth Boundaries (UGBs) are one of the most used policy tools around the world to control urban sprawl. The practice of UGBs can be traced back to the Greater London Plan in the 1930s, when large-scale “green belts” were established on the edge of cities to limit urbanization to a defined area [3]. The concept of smart growth, introduced to the United States in the 1990s to control uncontrolled urban sprawl, emphasizes environmental, social, and economic sustainability and is a compact, centralized, and efficient program. Under this concept, UGB was proposed as an important policy tool and has since been used in many countries around the world, including the United States, Canada, Japan, and the United Kingdom [3]. Chinese government reformed the spatial planning system and proposed a new spatial planning requirement, including the implementation of the “Three Districts and Three Lines” project [4], where delineating UGB plays a key role in regional spatial planning.

As far as the current research is concerned, cellular automata (CA) are a “bottom-up” dynamic model that is discrete in time and space. Thus, it is capable of simulating the spatiotemporal evolution of complex systems [5]. In recent decades, this model has been widely used by researchers in land-use simulation and urban sprawl studies, thereby becoming an important technical approach for delineating urban growth boundaries [6,7,8,9]. The research methods based on the CA model can be mainly divided into logistic regression (LR) and Markov chain (MC), as well as CA combining artificial intelligence (AI) algorithms [10]. However, given that LR is a linear simulation method that is not able to simulate complex and nonlinear changes in land use. Hence, recent research focuses on MC and AI methods. CA-MARKOV simulations require a dataset of land cover representations of the initial state, a Markov transition matrix, a set of suitability images (one for each land cover class), multiple iterations, and an adjacency filter. Therefore, transition rules were set up to develop suitability maps for each simulated land cover using multi-criteria evaluation (MCE) and fuzzy membership for land cover classes [10], which took into account the influence of various spatial influencing factors on change in land use, but failed to characterize their influence on the changing process of land use. The artificial intelligence algorithm combined with CA can take into account random and nonlinear processes in the process of land-use change, as well as the impact of multiple driving factors on the temporal and spatial changing process of land use. At present, the main methods used in this field include ANN-CA, Random Forest (RF), Support Vector Machine (SVM) [11], and so on—among which ANN-CA has been widely used by previous scholars, as it is based on artificial neural network (ANN). Considering that the ANN algorithm has proven to be an effective method to map the complex and nonlinear relationship between historical land use and various auxiliary data sources [10], it has proven to be highly feasible. For instance, Tayyebi et al. proposed a method based on an artificial neural network (ANN), GIS, and RS-UGB model to construct a cluster with 90 growth lines centred on Tehran to simulate the complex geometric environment of Tehran, Iran [8]; Liang designed a CA simulation based on the system dynamics (SD) model and an ANN to achieve a bottom-up urban growth simulation [12]. However, this model still has various disadvantages despite its advantages in processing complex data. For example, simulation accuracy is limited by the transformation rules and data mining scale, thus, it is necessary to extract the difference between the two periods of land-use data before simulation. Meanwhile, the transformation of land use used increases the computational scale of data and makes the model more complex [12]. Random Forest (RF) is a reliable non-parametric ensemble model that sits at the top of the classifier hierarchy. In addition, it uses a guided sampling strategy to create a “forest” of trees made up of various individual decisions. Each tree is based on a subset of feature variables and randomly selected observations. The final output RF model is a policy generated by averaging the decisions of individual trees or by voting [13]. Moreover, this method has various applications and has been extensively studied. Thus, a genetic algorithm was used to optimize the model, and an RF model to estimate the possibility of urban development. However, these existing models are insufficient in the ability to evolve a certain land-use patch and simulate the spatiotemporal dynamics of multiple land-use types. Therefore, Liang et al. proposed the PLUS model to solve these problems, which combines a new Land Expansion Analysis Strategy (LEAS) and CA model based on multi-type random patch seeds (CARS). This method is superior in use efficiency and simulation accuracy [14].

Furthermore, in most studies of UGBs, land-use maps derived from remote sensing are generally used as the only criterion for urban land identification. However, unfortunately, it is difficult to use them to reflect the use of urban land entities, and the specific conditions of human activities result in a certain degree of the unsuitability of the simulation of the urban growth boundary [15]. Herein, urban vitality is considered to be a necessary condition for urban success and one of the necessary characteristics to identify a city, with features such as diversity, high density, and hybridity [16]. In past studies, the amount of urban vitality was regarded as an important criterion for distinguishing urban and rural areas [17,18]. In addition, with the development of big data collection and processing technology, human activity data have become a major data source for identifying urban built-up areas in recent years, such as point of interest (POI) data, mobile phone signaling data, and water and electricity data. Meanwhile, Long et al. used mobile phone signaling to evaluate the implementation effects of urban growth boundary [15,17,18,19]. Moreover, inventory land growth has become an important model of urban spatial growth in many areas of the world, particularly for countries and regions in the middle and late stages of urbanization. Therefore, we need more accurate ways to identify areas that need to be developed or have potential growth in future urban expansion. We also regard the evolution of urban growth boundaries as a dynamic process to achieve more accurate and reasonable guidance for urban spatial growth.

The research focuses on the Chengdu Metropolitan Area (CMA) in southwest China. Over the past decade, it has been one of the most rapidly urbanizing regions in the world. In addition, its location on the eastern edge of the world’s highest plateau, the Tibetan Plateau, makes it the closest metropolitan area to the roof of the world; it has a considerable degree of complexity and fragility in the ecological environment. Hence, we must pay attention to the influence of the ecological environment on urban spatial expansion.

In practice, this study proposes an urban boundary delineation method based on urban human activities and urban vitality. The objectives are to identify areas with future development potential and existing urban low-efficiency land for the study area and incorporate them into the simulation of urban spatial expansion. Combined with the ecological security system of the ecological network, the random forest method is used in the PLUS model to calculate the future development probability of each land use. Moreover, the random seed CA prediction method is used under this constraint to simulate land use in three periods in the next 15 years. In this way, the urban and ecological spatial pattern of the study area is optimized, thereby delineating dynamic urban growth boundaries.

2. Study Area and Data

2.1. Research Area

The CMA is located in southwest China and on the east side of Longmen Mountain, on the eastern edge of the Qinghai–Tibet Plateau, which is the transition zone between the Sichuan Basin and the Qinghai–Tibet Plateau, as well as the world’s closest metropolitan area to the roof of the world. Sensitive areas, with more complex geological features and many seismic fault zones, make it more reasonable to formulate urban expansion strategies to avoid risky areas. The CMA consists of four cities, namely, Chengdu, Deyang, Meishan, and Ziyang, with an area of 33,136 km² (Figure 1).

In recent years, with the implementation of China’s Western Development Strategy and the Belt and Road Initiative, the region’s urbanization rate has increased from 50% in 2010 to 60% in 2017. In the following 15 years, with the implementation of the “Chengdu–Chongqing Double Cities Economic Area” strategy, regional development will become an important step in China’s inland economic growth. Hu presented ideas for Chengdu’s planning, such as balancing urban and rural development, focusing on efficiency and equity, paying more attention to citizens’ interests, and coordinating short- and long-term development [20] Therefore, the CMA needs a UGB policy to define the scope of UGB and rationally guide the development and utilization of its urban land. The area of built-up urban areas in China has nearly tripled in the past 15 years, reflecting the significant effect of CMA’s urbanization. However, the expansion of urban land into metropolitan areas has occupied agricultural land and ecological space and adversely affected the sustainability of regional development.

2.2. Sources of Data

The sources of research data are divided into three categories: land-use map derived from remote sensing, spatial big data, and other spatial impact factor data. Land-use remote sensing data come from LANDSAT remote-sensing image data as an information source and are obtained through manual visual interpretation. Considering the large scope of the study area, the resolution of the data was adjusted to 300 m × 300 m. As for the spatial big data, POI data from the AutoNavi map open data platform, involving business, catering, scenic spots, commercial services, public services, and other points of interest, are geospatial big data representing real geographic entities. Each piece of data contains attribute information, such as the name, address, type, latitude and longitude, and administrative area of the geographic entity. The POI data used in this article come from the Gaode Map Open Platform (https://map.gaode.com/ (accessed on 10 July 2022) in 2020, and the AMappoi tool is used to collect the POI data. Each piece of data contains information, such as name, address, and category, as well as latitude and longitude coordinates. The heat map comes from the Baidu map open platform, which reflects the population density clustered in each region. Finally, the calculation results are displayed in different colors and brightness. Considering the spatial difference in human flow, the Baidu heat map has strong timeliness and is updated every 15 min. Given the cost of data acquisition, the thermal data on 1 July 2020, were used. In addition, the night vision remote sensing data came from the NPP-VIIRS satellite, and the data were carried by the Suomi National Polar-orbiting Partnership Satellite. Provided by the Visible Infrared Imaging Radiometer (VIIRS), this satellite has a higher spatial resolution than others, it can detect a weaker light radiation [21] and its spatial resolution is 500 × 500 m. The population density data are derived from the WorldPop population dataset.

Other spatial image factor data mainly include spatial data, such as elevation, slope, land erosion degree, NDVI vegetation coverage, location of seismic fault zone, geological disaster occurrence, rivers, roads, ecological protection areas, and planned green corridor data sets. These data are provided by Data Center for Resources and Environmental Sciences, Chinese Academy of Sciences (RESDC) (http://www.resdc.cn (accessed on 12 July 2022)).

3. Methods

The research is divided into three steps: constructing an ecological network based on ecological sensitivity evaluation, identifying the scope of built-up areas, and proposing land-use classification, and simulating the evolution of land use based on PLUS, thereby delineating boundaries based on the simulation results (Figure 2).

3.1. Building an Ecological Network

Owing to the special geographical location of the study area and the sensitive ecological environment, the protection of ecologically sensitive areas is included in the urban growth forecast within the region. The ecological network concept is a suitable basis for inserting biodiversity conservation into sustainable landscape development, which can bridge the paradox between reserve conservation and urban development. The networks that provide connectivity offer habitats and corridors that help conserve biodiversity [22].

3.1.1. Ecological Sensitivity Evaluation

The ecological sensitivity of the research area must be evaluated to preliminarily analyze the natural conditions, social economy, and policy factors of the research area and provide suitable land for future urban land expansion. The land not suitable for urban expansion was screened out through evaluation. Combined with the actual situation of the research area, 10 indicators (e.g., topography, ecology, hydrology, land use, and land coverage) involved in selecting the evaluation system were determined [23,24,25]. Notably, the study area is located on the eastern edge of the Tibet Plateau. In addition, geo-logical disasters, such as earthquakes, occur frequently [26], which have a significant impact on the local ecological sensitivity and urban development. Therefore, the frequency and spatial distribution of geological disasters and the distance from the fault zone are also important indicators in the evaluation. Finally, the indicators are divided into three to five grades by the Natural Break method, and the weights are determined by the analytic hierarchy process (

Table 1. Ecological sensitivity evaluation.

Evaluation Items	Weight	Grade
		9	7	5	3	1
Elevation (m)	0.15	2864–6507	1859–2864	1061–1859	589–1061	225–598
Slope (%)	0.15	41.89–82.16	24.1–41.89	13.21–24.1	4.51–13.21	0–4.51
Land Relief (m)	0.1	1555–4980	578–1555	173–578	56–173	0–56
Distance to rivers and lakes (m)	0.05	0–12075	12,075–28,980	28,980–48,784	48,784–73,901	73,901–123,168
Distance to seismic fault zone (m)	0.1	0–64332	64,332–96,783	96,783–120,972	120,792–18,993	18,993—264,422
Frequency of geological disasters in 2020 (times)	0.15	4785–5982	3589–4785	2392–3589	1196–2392	0–1196
NDVI vegetation coverage rate (%)	0.1	0.63–0.99	0.52–0.63	0.40–0.52	0.23–0.40	0–0.23
Soil Average erosion modulus (t/km2 ·a)	0.1	8000–15,000	5000–8000	2500–5000	1000–2500	<200
Land-use types	0.1	Bare land, wetland	--	Farmland, woodland, grassland	--	Urban land
Identified or planned ecological spaces	Natural reserves, forest parks, planed green corridor

).

Based on the ArcGIS, the superposition analysis was carried out according to the above-mentioned determination factor superposition method, and the results of the comprehensive ecological sensitivity evaluation of the study area (

Table 2. Classification of ecological sensitivity evaluation results.

Grade	Value
Very low	25–59
Low	60–75
Medium	76–95
High	96–119
Very high	119–181

) were obtained. Five grades were used, and the Aggregate Polygons tool was used to aggregate the highly sensitive patches with a distance of less than 1000 m to form the ecological bottom-line area of the study area.

3.1.2. Construction of Ecological Network Based on LCP methods

Recently, the least cost path (LCP) method has been widely used in the process of ecological network construction. This model considers the geographic information of the landscape and the behavior of organisms. By referencing relevant literature [23,24,25,26] to the identification of ecological sources, evaluation of landscape resistance, extraction of ecological corridors, and setting of the corridor width, we propose the following ecological network construction methods:

• Identification of ecological sources. The patches of the high and extremely sensitive areas were sorted from large to small. In addition, the patches with an area ranking in the top 10% and area not less than 100 ha were selected as ecological sources. The remaining patches were defined as ecological targets.
• Evaluation of landscape resistance. Landscape resistance refers to the difficulty of species migrating among different landscape units. The higher the ecological suitability of the patch, the lower the migration cost. When assigning the landscape resistance, the landscape resistance value of each ecological unit in this study was determined with reference to the estimation of the ecological diversity value land of various ecological land in China by Xie et al., and other related studies [27] (

  Table 3. Value of migration cost for various land-use types.

  Land-Use Type	Cost
  Urban land	1000
  Woodland	1
  Grassland	3
  Farmland	20
  Wetland & rivers	300
  Roads	1000

  ).
• Extraction of ecological corridors. The cost distance tool in ArcGIS10.5 is used to calculate the path cost from all ecological sources to the ecological target, and then, the cost path tool is used to extract ecological corridors.
• Ecological corridor planning. The green corridor is adjusted according to the actual situation and divided into three types: landscape leisure green corridor, ecological protection green corridor, and ecological greenway along the road.

3.2. Identify Urban Built-Up Areas

3.2.1. Generating Kernel Density

The accurate extraction of urban built-up areas has great practicability and is of great significance for measuring the urbanization process and judging the urban environment. However, measuring urban built-up areas with a single point of the dataset is difficult [7]. Studies in recent years proved that urban vitality measurement indicators represented by POI data could better explain the insufficiency of extracting night-time light data (NTL) and could make the extraction of urban built-up areas more objective and accurate [28]. Heat map data and population density data can extend the perspective to specific human activities, enhance the accuracy of identification, and make up for the problems of low spatial resolution of NTL, including the inability of POI data to reflect human activities. After referring to relevant studies, the four data were given equal weights [18,28]. Additionally, the influence of abnormal pixels on the results was excluded through geometric correction of NTL data and population density data. Compared with other data, the preprocessing of POI point data is more complicated. To obtain the density distribution of POI data, the kernel density estimation method is used, which is a high-quality density estimation method by estimating the density of point features in the surrounding neighborhood through the kernel density function [29]. The POI kernel density was generated by ArcGIS, and the same method was used to obtain the heat map data. Finally, the superposition analysis was used to synthesize the four data to obtain the comprehensive density data, that is, the city’s vitality value distribution (Figure 3), whose range is between 0 and 10,822.

3.2.2. Extraction of Urban Vitality Area

The generated comprehensive density data are a smooth surface. To obtain more reflections on the relationship between the density and the contour area, the ArcGIS extraction contour tool is used to set contour lines with an interval of 100 m (Figure 4a). After referring to the relevant literature [18,28,29], the Density-Graph (D-G) method was used to identify the vitality boundary between urban and rural areas. In urban areas, the comprehensive density is high, and the contour lines are relatively tight, whereas the density in rural areas is relatively low, and the contour lines are sparser. Given this feature, the boundary between urban and rural areas is identified by observing the changing law of density and contour growth. Among them, the area of the closed curve enclosed by the contour line is $S_d$, its theoretical radius is $\sqrt{S_d}$, the density value is the abscissa, and the radius increment is the ordinate to fit the D-G curve [26]. Ultimately, the goodness of fit is R² > 0.9 (Figure 4b). As the density value decreases, the curve keeps rising until a certain point where the rising rate accelerates. The graph shows that the point lies between the density values of 7000–7500. Compared with the urban area identified by the remote sensing image, the range with a contour value of 7500 has a high degree of coincidence with the recognition range. Thus, it is defined as the dynamic boundary between the main urban area and the rural area.

In addition, we observed that for some independent towns outside the main urban area, we could not extract from them with a unified D-G method because of the large gap between their density values and the main urban area. Therefore, for such towns, the Density Counter Tree method is used to extract the city boundaries, which is a method developed based on local contour trees. This method uses one or more trees to quantitatively represent a contour map [19].

3.2.3. Identify Land-Use Types

We compare the identification results with the remote sensing images of urban land and combine the differences between the two identification results and the actual situation. Then, we propose the following five new definitions of land-use types:

• Urban vitality land (UVL): The identified vitality areas overlaps with the remote sensing image of urban land.
• Urban development-potential land (UDPL): UDPL is the land that belongs to the urban vitality area but not defined as urban land in the remote sensing image and does not belong to the ecological network. This land-use type is the main source of the next phase of urban expansion.
• Urban inefficient land (UIL): UIL belongs to the urban land in the remote sensing image but the vitality value does not meet the urban standard. This area may be a new urban area under development or a decaying industrial area. We advocate the renewal of this area to promote the restoration and improvement of vitality.
• Ecological network land (ENL): ENL including the identified ecological network land, rivers, and lakes prohibit urban development, which includes bare land, forest land, cultivated land, grassland, wetland, and other land-use types.
• Non-urban land (NUL): Except for all the land-use types mentioned above, NUL, including forest land, allows urban development, village construction land and others. As this type of land use does not involve restrictive factors, such as ecological protection, for the convenience of simulation, it is unified into the same type of land use.

3.3. Urban Expansion Simulation Based on PLUS Mode

To simulate the expansion of various types of land more accurately, we adopted the Patch-generating Land-Use Simulation (PLUS) model developed by Liang et al. [9].

The PLUS model contains two modules: a random forest framework based on a land expansion analysis strategy (LEAS); and a CA based on multi-type random patch seeds (CARS) [14]. The model was available to simulate the change of land-use patches and to analyze the underlying drivers of land-use dynamics by applying a combination of the LEAS and CARS. The resulting output by the PLUS model has proved more accurate, had more reliable spatial patterns, and allowed for important insights concerning the drivers of land expansion [14,30,31,32].

In this study, the PLUS model will be used to calculate the expansion probabilities of the five land use types and to generate the patches with high development probability for each land-use type to simulate the land-use changes. In addition, we take five years as a period and control the number of patches of different land-use types in each period based on different growth strategies in different periods to achieve the goal of dynamic simulation.

3.3.1. Land-Use Expansion Analysis Strategy (LEAS)

The proposed LEAS requires two dates of land-use data as its TAS, thus, in this study, the land-use maps of 2015 and 2020 were used. We overlaid the two periods of land-use data and extracted the cells with changed states from the later date of the land-use data, which represented the change regions for each land-use type. Then, sampling points were randomly selected and divided into subsets based on their land-use types, which were separately analyzed using a data mining method [14].

The random forest algorithm is used to generate the factors of each type of land-use expansion and driving force one by one to simulate the appearance of each land-use type in the same geographic space probability [30]. Referring to the previous literature [9,14,30,31,32], we selected five spatial driving factors: the distance to the road, the distance to the railway, the distance to the rail traffic, the elevation, and the slope. The random forest parameters in the LEAS module are set as follows: using the uniform sampling method, the number of decision trees is 20, the sampling rate is 0.01, and the number of features used to train the random forest is 5, which is the same as the number of driving factors. The constraint map is a binary image, where 1 indicates that the land-use type can be converted to other land-use types, and 0 indicates that the land-use type cannot be converted to other land-use types (Figure 5f). The ENL and river wetland land were selected as the prohibited construction zone. Finally, we obtained the development probabilities of five types of land use (Figure 5a–e). As we can see, the development of urban vitality areas is more likely to be closer to the periphery of the city, whereas the development of ecological land is far from urban areas.

3.3.2. Predicting Future Demand for Land Expansion

For the future demand simulation of urban land, we forecast in three stages: 2020–2025, 2025–2030, and 2030–2035. We use the Markov method to forecast for each land-use type and adjust for different expansion strategies in the forecast stage (

Table 4. Land expansion strategies and demands at each stage.

Period	Expansion Strategies	Prediction of the Number of Grids for Each Land-Use Type
		UVL	UDPL	UIL	ENL	NUL
2020–2025	Incremental expansion first	13,136	4932	4477	193,497	154,597
2025–2030	Incremental and inventory expansion equally valued	18,293	3565	3352	190,823	154,597
2030–2035	Focus on inventory development	21,704	2194	2302	190,362	154,597

). In the past five years, Chengdu Metropolitan Area’s (CMA) expansion strategy has been based on inventory development, focusing on the development of potential urban land [33]. In the 2025–2030 time period, urban land redevelopment and urban expansion are equally important. Taking a strategy, such as promoting the diversification of land use and enhancing the vitality of the city, is encouraged. In the time of 2030–2035, the urban expansion will slow down, and more emphasis will be placed on the development of UIL.

The CARS module is a CA model that includes a patch-generation mechanism based on multi-type random seeds of land use [14], which integrates the impacts of macro “top-down” land-use demand and the “bottom-up” simulation on the land system. It incorporates an innovative multi-type random seeds-generating mechanism to simulate micro-land-use competition to drive the current land-use amounts to meet the macro-demand under the neighborhood effect and development probability (the result of LESA) [32].

Table 5. Neighborhood effect weights of various types of land-use.

Land-Use Type	Weights
UVL	0.2
UDPL	0.25
UIL	0.2
ENL	0.15
NUL	0.2

shows the neighborhood effect weights of various types of land use in CMA, represents the neighborhood effect of a pixel which is determined by the proportion of land use of type in the neighborhood of a pixel and the neighborhood weights [32]. When the neighborhood effect of type land use is zero, the multi-type random seeds-generating mechanism will generate random seeds of each land-use type through the Monte Carlo method. Additionally,

Table 6. Transformation matrix of various types of land-use.

	UVL	UDPL	UIL	ENL	NUL
UVL	1	1	0	0	0
UDPL	1	1	0	0	0
UIL	1	1	1	0	0
ENL	1	1	1	1	0
NUL	0	0	0	0	1

shows the transformation matrix of various types of land use the value is 0 or 1, the value 0 means it cannot be converted with a specific land-use type, and 1 means can be converted. Owing to the characteristics of ecological land, it does not participate in any land conversion.

3.3.3. CA Prediction Based on Multi-Class Random Patch Seeds (CARS)

CARS module is a CA model that includes a patch-generation mechanism based on multi-type random seeds of land use [14], which integrates the impacts of macro “top-down” land-use demand and the “bottom-up” simulation on the land system. It incorporates an innovative multi-type random seeds generating mechanism to simulate micro-land-use competition to drive the current land-use amounts to meet the macro-demand under the neighborhood effect and development probability(the result of LESA) [14].Table 5 shows the neighborhood effect weights of various types of land use in CMA, represents the neighborhood effect of a pixel which is determined by the proportion of land use of type in the neighborhood of a pixel and the neighborhood weights [31]. When the neighborhood effect of type land-use is zero, the multi-type random seeds generating mechanism will generate random seeds of each land-use type through the Monte Carlo method. Additionally, Table 6 shows the transformation matrix of various types of land use, the value is 0 or 1, the value 0 means it cannot be converted with a specific land-use type, and 1 means can be converted. Owing to the characteristics of ecological land, it does not participate in any land conversion.

3.3.4. Model Validation

To verify the prediction accuracy, Kappa and FOM coefficients are used to check the simulation accuracy. Among them, the Kappa coefficient is a test method used to evaluate the consistency of the impact classification results. The larger the value of the Kappa coefficient, the more accurate the result of the land-use simulation [30]. The FOM coefficient is a performance evaluation method, which is the ratio of the accurately predicted sample to other samples. A larger FoM value means higher simulation accuracy [30].

3.4. UGB Delineation

The ecological network is regarded as a restrictive growth boundary, that is, a rigid growth boundary. Although the use of rigid boundaries brings uncertainty, it is a useful addition to environmental protection [34]. For the simulated urban vitality land (UVL) at each stage is the extent of the flexible urban growth boundary. The UIL and UDPL in each stage are the main space for future urban expansion and redevelopment, which both participate in the next stage of the space expansion process. The role of the flexible boundary is to guide urban expansion to take place within the boundary to form a compact urban space. To form a smooth border, the edge filter is used to remove broken patches and sharpen the edges of urban cells and convert them to polygons to form smooth boundaries.

4. Results and Analysis

4.1. Ecological Network Construction Results

Through the ecological assessment, 5588.9 km² of high and extremely sensitive areas were identified (Figure 6a). The ecological source patches are concentrated in the Longmen Mountain and Longquan Mountain, and 15 ecological corridors were established for the other six ecological targets. The width of the ecological corridor varies according to its different functions. The purpose of the ecological protection corridor is to protect ecological diversity and reserve enough space for species migration. According to the corridor width proposed by Roling J, the width of such corridors was 60 m, and the width of the landscape leisure green corridor and the ecological greenway along the road is 5–15 m according to the relevant planning [35] (Figure 6b). The reason is that such land is concentrated in the city center, and the land supply is relatively tight. The ecological network is composed of ecological sources, ecological goals, and ecological corridors, which are also the rigid boundaries of urban expansion. They are distributed around the study area and rely on Longmen Mountain and Longquan Mountain. Each patch is large and concentrated. The Mt. Longmen fault zone in the transition zone between the Chengdu Plain and the Qinghai-Tibet Plateau is the largest patch, providing a large area of ecological resources and ecological security barriers for the urban area. Its network structure is relatively simple, and Longquan Mountain assumes the role of the ecological center, forming a spatial pattern connecting the east and west with the Longquan Mountain ecological area as the core, ensuring the connectivity between various ecological patches. The main components of the ecological network are cultivated land, forest land, grassland, bare land, and wetlands. Among them, cultivated land accounts for 52.4%, which is the largest source of land for the ecological network, followed by forest land, accounting for 38.6%.

4.2. Land-Use Identification Results

A total of 672 km² of land were identified as UVL, and the urban low-utility land area was 412 km² (

Table 7. Area of various land-use types after identification.

Land-Use Type	UNL	UIL	UDPL	ENL	NUL
Area (km2)	672	412	850	17,509	13,914

), accounting for 38% of the urban area identified by remote sensing. This area is located on the outer edge of the central urban area of Chengdu, particularly in the east and northeast of Chengdu, and the southwest. There are more UIL in the central and western regions. At present, the types of land used in these areas are low-density residential and industrial areas (Figure 7). Problems exist, such as low development density, single function, disorderly expansion, and poor appeal to the crowd, and they have great urban vitality and increased demand. UDPL is concentrated on the outer edge of the central urban area of Chengdu. This area is close to the main urban area, has close economic ties, and has good traffic conditions, which is the main direction for future urban expansion. Moreover, the Jianyang area (Figure 7) located on the south side of Longquan Mountain is close to Tianfu International Airport. Its comprehensive urban vitality is higher than in other areas, and it will carry more urban functions in the future. Therefore, this area has more potential development land.

In general, the current CMA urban land-use efficiency is not high, and there are a large number of low-efficiency lands that need to be improved. For a long time, CMA’s urban development has ignored human activities and urban vitality, particularly the lack of reasonable land-use policies, which is mainly reflected in the random conversion of another land into urban land. For these reasons, urban functions cannot be realized, which enhances suburban sprawl.

4.3. Delineation Results of Urban Growth Boundaries

The tree-stage dynamic urban growth boundaries are drawn combined with the PLUS simulation results (Figure 8 and

Table 8. Area of various land-use types after simulation.

Land-Use Type/Year	UNL	UIL	UDPL	UNL	ENL
2025 (Km2)	1182	443	403	17,415	13,914
2030 (Km2)	1646	321	302	17,174	13,914
2035 (Km2)	2043	287	297	16,816	13,914

). The kappa coefficient of the simulation is 0.948, the overall accuracy is 94.75%, and the FoM index is 0.122, which is quite efficient for further analysis. The resulting UGB raster format was converted into a vector format, patches smaller than 5 km² were removed, and patches within 2 km were aggregated to form a continuous growth boundary. Figure 8a–c show the specific conditions of UGB in three different locations in three time periods. Different spatial expansion strategies were adopted in each time-period. The low-vitality areas of towns decreased significantly over time, demonstrating the effectiveness of the simulation with a process of inventory development behaviors. Moreover, most of the urban-potential development land is converted into UGB. From the specific perspective, the central urban area of Chengdu (Figure 8a) will be mainly used for urban development potential in the space expansion in 2020–2025, leaving a large number of urban low-utility land waiting to be improved. This situation will be improved in 2030–2035, as inefficient land continues to decrease and urban land further expands in the northeast direction, gradually integrating with the Deyang area. The Dujiangyan area (Figure 8b) located in the northwest of CMA has virtually not experienced a spatial expansion in the 15-year expansion plan. The reason is the stagnation of spatial expansion caused by too many restrictive factors in the surrounding area, which is inconsistent with the actual situation and needs to be reconnected. Thus, further improvements will be made in the following research of the Jianyang area which relies on the international airport; its (Figure 8c) space expansion is fast and has been in the process of acceleration.

4.4. Analysis

4.4.1. Analysis of Land-Use Conversion

Figure 9a shows the flow of land-use conversion during the simulation period from 2020 to 2035. Owing to the new land-use classification, most cultivated land has flowed to NUL and continued to exist as cultivated land. A considerable part of the cultivated land has been converted into ENL, reaching 7083 km². Nearly 80% of the forest land has been converted into ENL, which is an important part of the CMA landscape pattern. The land within the UGB range mainly comes from cultivated land, non-urban construction land, and original urban land. Figure 9b shows the sources of transformation between dynamic urban land, UIL, and UDPL in three stages. The growth of UVL land mainly comes from NUL and UDPL. With the implement of planning, the dependence of new UVL on NUL is gradually reduced. By the end of the planning period, the number of UIL and UDPL has decreased significantly, which means that most of the UDPL has been successfully converted into UGB. UIL has also been successfully converted into UGB in this process, reflecting the success of land inventory development strategy. Through the analysis of land-use transformation, the transformation of urbanized areas with different vitality values in the simulation can be seen. Thus, our proposed UGB delineation method regards human activities as an important factor affecting urban expansion, follows complex mathematical rules, and regards the simulation of urban growth boundaries as a dynamic process. Moreover, advocates that low-utility land should be transformed into urban vitality, which increases the practicality of UGB during implementation.

4.4.2. Expansion Direction Analysis

As shown in Figure 10, an eight-quadrant map is drawn to compare the expansion directions of the 2025–2035 simulations for four important cities. For the downtown area of Chengdu (Figure 10a), the first to eighth quadrants have almost the same expansion trend. In the Deyang urban area (Figure 10b), the expansion patterns at different stages differ. In 2025–2020, there will be developed in all directions. In 2030–2035, with the development of north and east, there will be evident development in quadrants II, V, and VI. For the Jianyang area (Figure 10c), 2030–2035 year-on-year growth dominates, particularly in quadrants V, VI, VII, and VIII. The development from 2025 to 2030 is mainly concentrated in the I, VII, and VIII quadrants. The growth in the south and north directions is more evident. For the Meishan area (Figure 10c), the growth is more concentrated on the west side, whereas quadrants V, VI, and others in the east direction show little increase.

5. Discussion and Conclusions

This study uses a set of research methods to predict the spatial expansion of CMA and delineate dynamic UGBs for it. Combined with the identification of UIL and UDPL, the two are included in the urban spatial land. In the prediction simulation of dynamic expansion and evolution, the ecological network construction technology method is used. Moreover, the ecological barrier of the CMA is constructed as a rigid boundary to control the growth. The technical methods include the PLUS model, LCP method, kernel density method, and D-G model. In short, the model can be applied to cities or metropolitan areas that are in the middle and later stages of urbanization and have certain spatial continuity characteristics.

Compared with the current UGBs delineation technical methods, this study has made certain progress: The urban vitality and land coverage are combined to consider the urban land expansion probability, thereby highlighting the dynamic transformation process of urban inefficient land and development potential land, which responses to the development strategy of the inventory expansion in the study area; It highlights the constraints derive by ecological pattern, ensures the ecological security and provides a considerable amount of green infrastructure and open space for cities and towns in the metropolitan area. However, the specific implementation results need further practice verification.

Owing to the data collection and processing technology of big data, through processing, the situation of human activities in the geographical space over a period of time can be clearly identified. Moreover, through a comparison with remote sensing images, inefficiencies and vitality in towns and cities can be distinguished. Insufficient areas can also be identified as NUL with development potential. As a dynamic policy tool, UGB should not blindly pursue urban extensional growth in the implementation process but should select appropriate land for appropriate extensional expansion and redevelopment of inefficient land. This process should also reflect in the process of UGB delineation. This study selects the CMA as a case, identifies a large number of urban low-efficiency land and development potential land, and incorporates them into expansion simulation and UGB delineation. The study also combines ecological sensitivity evaluation to construct an ecological network as a rigid growth boundary has formed the urban land development and control strategy and ecological landscape pattern in the next 15 years. Over the course of a 15-year dynamic simulation, inefficient lands were continuously reduced, meeting a considerable portion of the urban expansion needs.

When the local government is using the results of this research, it can be combined with the implementation of territorial spatial planning. Based on China’s newly revised “Land Management Law”, territorial spatial planning should coordinate the layout of agricultural, ecological, urban, and other functional spaces, delineate, and implement permanent basic farmland red lines, ecological protection red lines, and urban development boundaries. Our study result can serve as a reference for the delineation of urban development boundaries of the Chengdu metropolitan area. In the process of delimiting the urban development boundary, we advocate that in the process of delineating urban development boundaries, five years are used as a planning period, and different development strategies are adopted for each stage to delineate growth boundaries based on the actual conditions (Figure 10). Herein, the balance between the two development strategies of inventory development and incremental development in different stages was highlighted: In the first five years (2020–2025), the incremental development strategy will continue. The area around the central urban area of Chengdu and the Jiayang area has a higher expansive potential, and the government should give priority to the development potential areas of related areas. Meanwhile, in the second stage (2025–2030), the government should focus on guiding the slowdown of incremental development, and pay attention to the renewal and development of areas in the Northeast, West, and South of Chengdu, and take measures to improve traffic accessibility, enhance the diversity of urban land functions and land-use Intensity, and so on, to enhance the urban vitality of some regions, hence, to improve the carrying capacity of the urban population. Moreover, in the third stage (2030–2035), it is expected that the urbanization process of the study area will enter the later stage, and the incremental development will further slowdown, particularly in the Central urban area of Chengdu. Except for some areas (such as Jiayang, Measham, Deyang, etc.) that still have the potential for extensional growth, the development mode of other areas should be adjusted to the inventory development mode. However, areas with low density will be redeveloped to further improve the urban structure and functional layout. In addition, the delineation of the ecological network is from the perspective of the ecological function of the study area, and further extends the functionality of the ecological protection boundaries for ecological function assurance, environmental quality safety, and resource utilization, as well as further guides urban expansion in the appropriate range. Then, it can form a complete network with both ecological protection and recreational functions.

In the research, some problems still need to be solved. First, the identification methods for low-utility land are relatively simple. Owing to the large scope of the study area, this study lacks detailed research and observation on microscopic land-use and lacks investigations on the land-use properties, plot ratio, building density, and traffic conditions of relevant land parcels. Hence, some land use in special circumstances is included in the low-efficiency land. Therefore, when making a development plan, the land within the development scope should be re-identified first. Second, the identification and simulation of small urban patches will have some deviations. Considering their large differences in scale and vitality from the surrounding central urban areas, using the same parameter settings may reduce the accuracy of the simulation. Therefore, this research framework is more suitable for simulating a single continuous, highly concentrated urban area. Moreover, separate simulations are required for individual, relatively independent patches.