Growth Simulations of Urban Underground Space with Ecological Constraints Using a Patch-Based Cellular Automaton

Abstract

The growth simulation of urban underground space (UUS) under the consideration of ecological constraints can effectively reveal the characteristics and trends of UUS changes, and provide a basis for planning the construction of sustainable and livable ecological cities. Therefore, this study considers urban ecological space as a constraint mechanism for UUS development and conducts a simulation study of the dynamic and complex UUS growth process, with a view toward guiding UUS planning under a long-term overall vision. In this study, a patch-based cellular automaton (CA) model is constructed to simulate the dynamic and complex growth process of UUS, subject to the ecological constraints generated by the agent-based land allocation optimization model. The spatial drivers of UUS growth simulation are determined based on the Random Forest (RF) algorithm. The results of the research case in Tianfu New District, Chengdu City, demonstrate that UUS expansion with ecological constraints exhibits sustainable characteristics. However, the growth rate of the UUS development scale is significantly lower when ecological constraints are present compared to when they are not. This study’s results contribute to urban management by finding a balance between UUS development and ecological space conservation, and providing theoretical support for rational UUS planning and decision making in the construction of low-carbon cities.

1. Introduction

The three-dimensional development of cities using urban underground space (UUS) is relevant to the world’s leading metropolitan areas. Cities in developed and developing countries are expanding and growing much faster than their populations [1,2]. Sustainable cities are pursuing compact city policies to exclude new constructions in undeveloped and green areas [3,4]. UUS is one component of the three-dimensional urban space (including above-ground, surface, and underground space). It is a unique space with a vast potential for exploration in achieving the goal of low-carbon cities, which has not been fully utilized [5,6,7]. The results of many studies have demonstrated that the rational and effective use of UUS is beneficial for solving urban development problems, especially for achieving urban sustainability [5,8,9].

UUS planning should be guided by a far-reaching overall vision to ensure that these spaces can be developed and utilized in an orderly and sustainable manner, as it is difficult to demolish or reconstruct UUSs once they have been developed and utilized in a short period [10,11,12,13,14]. In cities worldwide, there is usually a plot of land (either large or small) with a high density that needs to be developed underground, and there is still much space that can be used underground compared with the surface space in most cities [15,16,17,18]. Currently, utilizing UUS usually follows the principle of “first-come-first-served”, whether for underground resource exploitation or physical structures [5,19,20]. All planners and policymakers need to be aware of the opportunities and threats of using UUS to improve urban livability and establish well-thought-out UUS planning systems [6,7,11,21].

The results of modeling the growth simulation of UUS provide a reference for scientific UUS planning and sustainable UUS infrastructure construction. The growth simulation of UUS is an essential tool and method to effectively guarantee the coordinated planning of urban aboveground–underground integrated space [9,12,22,23], the scientific and organized development of underground space and resources [24,25], and the environmental and sustainable construction of UUS infrastructure [26,27,28]. The growth simulation of UUS takes the historical UUS utilized distribution as the critical input and combines traffic conditions [29,30,31], the geographical environment [10,24,32], and socioeconomic [16,33] and other driving factors, and adopts an empirical approach to simulate the spatial and temporal development of UUS use.

The mechanism of cellular automatons (CAs) makes them appropriate for UUS growth simulation because they can describe the spatial interactions in neighboring or layered areas. A growth simulation of UUS, similar to urban space growth simulation, focuses on the self-organization and spatial evolution of UUS or urban space development, where CA models have been widely used in urban space growth simulation [34,35,36] to explore the urban expansion mechanisms of various cities [37,38,39]. CA models can also effectively aggregate various socioeconomic and geographic factors, including demographics, traffic, and location, to simulate spatial utilization dynamics by defining reasonable transition rules based on spatial interactions between cell states and regions [40,41]. In brief, CA models contain a set of cells, each of which maintains a state that may be influenced by neighboring cells [42,43,44,45]. In particular, constructing and developing UUS should meet human needs for urban living, working, and traveling activities and protect the diversity of urban ecological space by ensuring that its area is not reduced [8,11,22]. Modeling growth for UUS simulation not only excludes ecological and natural protected areas, such as water bodies and mountains, but should also include blue–green spaces (BGSs) that exist more widely in the city (e.g., urban forests, grasslands, and farmlands) [46,47].

However, no existing simulation growth model of UUS has explicitly incorporated ecological constraints into the dynamic feedback of its growing spatial configuration as it evolves. The research question in this study is two-fold: (1) How do we construct a model that can accurately simulate the process of UUS growth with ecological constraints? (2) What are the differences between the simulation growth results of urban subspace growth with and without ecological constraints? Therefore, this study aims to construct a growth simulation method for UUS under ecological constraints. The growth simulation of this method is established based on the patch-based CA (patch-CA) model, in which the ecological conservation zones (urban BGS protection areas) are based on the agent-based land allocation optimization (AgentLA) model. The spatial drivers of UUS growth simulation are determined based on the Random Forest (RF) algorithm. The novel contributions of this study, which differ from previous works, are as follows:

• To prioritize the importance of ecological conservation zoning in the future development of UUS use, the AgentLA model was used for ecological conservation zones (urban BGS protection areas) as a spatial constraint in the UUS growth simulation;
• A method based on the RF algorithm to determine the spatial drivers of UUS growth simulation was constructed based on the process and underlying mechanism of UUS development and the spatial driving factors considered in previous urban simulation studies;
• A method was constructed to simulate UUS growth based on the patch-CA model, in which the three spatial pattern types of UUS growth (filling, edge extension, and outlying) were simulated using two complementary procedures of organic and spontaneous growth in the proposed patch-CA model;
• UUS growth simulations under ecological constraints are evaluated using the Tianfu New District in Chengdu, China, as an example, and the validity and reproducibility of the model proposed in this study are verified.

The remainder of this paper is organized as follows: Section 2 provides a literature review, including a review of simulation models and constraint mechanisms of UUS growth from previous research. The framework, using the solution and proposed models, is presented in Section 3. The information and data of the study area are described in Section 4. Section 5 describes the case and corresponding calculation results, including the model parameters and UUS growth simulation results. The discussion and analysis are presented in Section 6. Lastly, Section 7 summarizes the study and presents its limitations.

2. Literature Review

2.1. Simulation Model of UUS Growth

The development of UUS utilization is always a dynamic and complex process, resulting from the constant temporal and spatial changes of natural processes and human activities [10,13,16]. In terms of the history of UUS utilization, the trend has shifted from a simple single space to a comprehensive space, such as a functional compound, network systematization, and deep layout [25,27,48]. Simultaneously, research on UUS is focusing more on its complexity and dynamic development [5,16,48,49,50,51,52].

Humans have intervened and regulated UUS development to various degrees in different periods and regions. However, the self-organization of the UUS system cannot be changed due to the uncertainty of UUS development, and because all the information on UUS development cannot be fully grasped [16,53,54]. Therefore, the UUS system can be considered a self-organization system, wherein the degree of self-organization works differently at different times and under different policy contexts.

The CA model has been proven to be an effective tool for simulating UUS development because of its high degree of consistency in dynamics and self-organization compared with the UUS system. The CA model is a discrete-time dynamic system defined in a discrete lattice state space that can simulate the spatio-temporal evolution of complex systems [55,56,57,58]. CA models are effective because of their flexible structure and support for bottom-up simulations [42,59,60]. CA models [34,35,36,37,38,39,40,41] have been widely used in urban growth simulation, land-use change (LUC), and analysis over the past two decades because they have the advantage of representing the drivers of urban development through transition rules and reflecting the spatial externalities of urban sprawl through neighborhood rules, as illustrated in Table 1.

Simultaneously, UUS development is a complex and dynamic system driven by various natural, social, and economic factors. Therefore, the characteristics of CA models and the above studies (see Table 1) demonstrate their similarity to the simulation of UUS development processes. Additionally, patch-based CA models [61,62,63] can simulate patch or parcel-based changes to reflect realistic UUS use change patterns, unlike traditional CA models that simulate cell-based changes. The UUS growth simulation in this study can be divided into three basic spatial modes (as illustrated in Figure 1): infill, edge expansion, and outlying. These three modes are simulated in the proposed CA model using two complementary procedures: organic and spontaneous growth [56,61,62,64,65,66]. The infilling and edge expansion that occur during UUS growth (which can be simulated by organic growth) refer to the expansion of the existing UUS patches that convert adjacent non-UUS cells by combining the neighborhood effect (Equation (9)) with the suitability of UUS occurrence in the spatial configuration module. Outlying UUS growth (which can be simulated by spontaneous growth) is intended to simulate leapfrog development by growing new patches that are not connected to existing UUS areas. New patches are generated as organic extensions in future UUS development, and existing and newly generated patches are combined over time to produce a diverse UUS development pattern.

Table 1. A representative list of literature on the urban growth simulation of CA models.

Number	Researcher(s)	Year(s)	Model(s)	Field of Study	Case Study
1	Mitsova et al. [44]	2011	CA-Markov chain model	Urban growth projections at a regional scale	Twenty-seven counties in Ohio, Indiana, and Kentucky, USA
2	Sakieh et al. [34]	2014	CA model	Simulating urban expansion and scenario prediction	Karaj City, Iran
3	Deep et al. [43]	2014	CA-Markov chain model	Measure the urban sprawl	Dehradun City, India
4	Guan et al. [35]	2015	Hybrid parallel CA model	Urban growth simulation	California, USA
5	Liao et al. [36]	2015	Logistic-CA model	Land-use simulations	Xiamen City, China
6	Dahal et al. [42]	2015	An irregular CA model	Urban land-use dynamics	San Marcos, Texas, USA
7	Gounaridis et al. [37]	2019	Random Forest—CA model	Explore future land use/cover change	Attica, Greece
8	Liao et al. [56]	2019	Logistic regression-ordered weighted averaging-CA model	Urban sprawl scenario simulations	Xiamen City, China
9	Shu et al. [38]	2020	Traditional logistic CA model	Simulate the urban growth	Tongshan County, China
10	Chen et al. [39]	2020	Geographically weighted regression—CA model	Simulate the urban expansion	Chongqing City, China
11	Roodposhti et al. [45]	2020	CA land-use models	Land-use change (LUC)	Ahvaz, Iran
12	Cao et al. [55]	2020	Logistic-CA model	Urban spatial growth	Hangzhou City, China
13	Zhou et al. [58]	2020	Random Forest and CA-Markov model	Multi-scenario simulation of urban land change	Shanghai City, China
14	Dinda et al. [40]	2021	CA-Markov Chain model	Urban growth pattern and loss in urban green space	Kolkata, India
15	Wang et al. [59]	2021	SM-Logistic-CA model	Simulating urban land growth	Beijing City, China
16	Liang et al. [62]	2021	Patch-generating land-use simulation model	Drivers of land expansion and landscape dynamics	Wuhan City, China
17	Liu et al. [41]	2022	ESP-FLUS Model	Urban growth boundaries	Min Delta Region, China
18	Jin et al. [57]	2023	CA-based Future Land Use Simulation (FLUS) model	Land-use change (LUC)	Myanmar
19	Lin et al. [63]	2023	Novel landscape-driven patch-based CA model	Urban land-use changes	Guangzhou City, China

Table 1. A representative list of literature on the urban growth simulation of CA models.

Number	Researcher(s)	Year(s)	Model(s)	Field of Study	Case Study
1	Mitsova et al. [44]	2011	CA-Markov chain model	Urban growth projections at a regional scale	Twenty-seven counties in Ohio, Indiana, and Kentucky, USA
2	Sakieh et al. [34]	2014	CA model	Simulating urban expansion and scenario prediction	Karaj City, Iran
3	Deep et al. [43]	2014	CA-Markov chain model	Measure the urban sprawl	Dehradun City, India
4	Guan et al. [35]	2015	Hybrid parallel CA model	Urban growth simulation	California, USA
5	Liao et al. [36]	2015	Logistic-CA model	Land-use simulations	Xiamen City, China
6	Dahal et al. [42]	2015	An irregular CA model	Urban land-use dynamics	San Marcos, Texas, USA
7	Gounaridis et al. [37]	2019	Random Forest—CA model	Explore future land use/cover change	Attica, Greece
8	Liao et al. [56]	2019	Logistic regression-ordered weighted averaging-CA model	Urban sprawl scenario simulations	Xiamen City, China
9	Shu et al. [38]	2020	Traditional logistic CA model	Simulate the urban growth	Tongshan County, China
10	Chen et al. [39]	2020	Geographically weighted regression—CA model	Simulate the urban expansion	Chongqing City, China
11	Roodposhti et al. [45]	2020	CA land-use models	Land-use change (LUC)	Ahvaz, Iran
12	Cao et al. [55]	2020	Logistic-CA model	Urban spatial growth	Hangzhou City, China
13	Zhou et al. [58]	2020	Random Forest and CA-Markov model	Multi-scenario simulation of urban land change	Shanghai City, China
14	Dinda et al. [40]	2021	CA-Markov Chain model	Urban growth pattern and loss in urban green space	Kolkata, India
15	Wang et al. [59]	2021	SM-Logistic-CA model	Simulating urban land growth	Beijing City, China
16	Liang et al. [62]	2021	Patch-generating land-use simulation model	Drivers of land expansion and landscape dynamics	Wuhan City, China
17	Liu et al. [41]	2022	ESP-FLUS Model	Urban growth boundaries	Min Delta Region, China
18	Jin et al. [57]	2023	CA-based Future Land Use Simulation (FLUS) model	Land-use change (LUC)	Myanmar
19	Lin et al. [63]	2023	Novel landscape-driven patch-based CA model	Urban land-use changes	Guangzhou City, China

2.2. UUS Growth Constraint Mechanisms

Constraint mechanisms are a critical element of urban CAs [67,68], and cells belonging to protected or planned UUS areas are subject to strict management development restrictions. Many scholars have carried out simulations of urban CAs based on constraint mechanisms to reflect the different constraints (e.g., water bodies, mountains, lakes, forests, and government policies) on UUS or land-use changes [56,69,70,71,72,73,74]. The cited studies demonstrate that the incorporation of ecological constraints into scenario simulations of UUS CAs produces significant results.

Ecological constraints have a significant role in simulating UUS growth, and a remarkable compensatory effect exists between these ecological landscapes composed of urban BGSs and UUS development. Several recent studies have confirmed that urban BGSs, comprising waterbodies, urban parks, forests, grasslands, and other green spaces, are the main sources of the urban ecological carbon sink that has strong carbon sequestration effects and capacity [75,76,77]. The plants in urban BGSs absorb CO₂ from the air through photosynthesis, which is then stored in the form of biomass, thus accumulating organic carbon [46,76,78,79]. The storage of carbon in vegetation and soil creates sustainable ecosystems in cities’ BGSs [80,81].

Three constraints—urban construction land boundaries, UUS restricted development boundaries (mainly considering geological constraints such as mountains and water systems), and ecological protection areas—hinder UUS development. In particular, the restricted boundaries of urban BGSs are considered to be constraints for ecological protection areas. In this study, an AgentLA model was applied to solve the ecological zoning problem because of the effectiveness and computational efficiency of AgentLA models compared with traditional methods [82,83]. The AgentLA model was utilized to generate a layout of the protected areas of urban BGSs, which were further incorporated into a patch-CA model as spatial constraints to simulate future UUS growth.

3. Methodology

This analysis of the simulation process of UUS growth can be summarized using three main stages (see Figure 2). The first stage is the generation of UUS development constraints: urban construction land boundaries, restricted UUS development boundaries, and ecological protection areas (urban BGSs). The restricted boundaries of urban BGSs are generated using an AgentLA model (as detailed in Section 3.1). The second stage (see Section 3.2 for details) is the basic simulation for constructing the CA model of UUS growth. The main task of this stage is to retrieve the transition rules that determine whether the cell state changes. The transition rules include spatial drivers as well as neighborhood effects. The third stage (as detailed in Section 3.3) is the construction of the UUS growth simulation using the patch-based CA model, which includes the following processes: suitability-based sorting, stochastic seeding, and the moving-window scanning process.

3.1. Generating UUS Development Constraints

Urban construction land boundaries, restricted UUS development boundaries (mainly geological constraints such as mountains and water systems), and ecological protection areas, are considered to be constraints on UUS development. Urban construction land and UUS restricted development boundaries are specifically identified by urban planning and geological survey departments; thus, in this study, these two UUS restricted areas are assumed to be clearly and specifically located in urban space. The AgentLA model was utilized to generate ecological reserves as a spatial constraint for the UUS growth simulation; that is, no development is allowed within the protected urban BGSs.

The AgentLA model optimized the spatial pattern of ecological reserves (urban BGS protection) to maximize the total ecological compactness and suitability. To implement the AgentLA model, the amount of available land to be protected within urban BGSs (i.e., the number of plots) should be initially specified. The model then randomly distributes individual agents (i.e., mobile units) in space. The agent’s population is equal to the given amount of land protected within the urban BGSs. Subsequently, each agent makes the decision to iteratively relocate based on the fitness function [83]:

$$f_{ij}=w_v{v}_{ij}+w_c{c}_{ij}$$

(1)

$$c_{ij}=\frac{{\sum }_{\Omega }{x}_{mn}{e}^{-\raisebox{1ex}{$d_{mn}$}\!\left/ \!\raisebox{-1ex}{$\gamma $}\right.}}{{\sum }_{\Omega }{e}^{-\raisebox{1ex}{$d_{mn}$}\!\left/ \!\raisebox{-1ex}{$\gamma $}\right.}}$$

(2)

where $f_{ij}$ denotes the fitness at cell $(i,j)$; $v_{ij}$ denotes the value of the ecological suitability in cell $(i,j)$; $c_{ij}$ denotes the neighborhood condition; $w_v$ and $w_c$ denote the weights of $v_{ij}$ and $c_{ij}$, respectively; $\Omega$ denotes the Moore neighborhood of cell $(i,j)$; $x_{mn}$ = 1 if cell $(i,j)$ is already occupied by an agent; otherwise, it equals 0; $d_{mn}$ denotes the distance between cell $(i,j)$ and $(m,n)$; and $\gamma$ denotes a compensation parameter that can be taken between $\left[\mathrm{0,10}\right]$.

The fitness of each agent first evaluates the current position ($f_{ij}$). Subsequently, it searches for the vacant position with the highest fitness among the neighboring positions adjacent to cell $(i,j)$ and can be denoted as $f_L$. Meanwhile, there is a randomly selected empty space from the whole cell space, which has at least one agent in its 3 × 3 neighborhood-adjacent position. The fitness of this position is denoted as $f_G$. The agents stay at or move to one of these positions with the highest fitness by comparing the values of $f_{ij}$, $f_L$, and $f_G$. During each iteration, the spatial pattern formed by the decisions of all agents is evaluated using the $F$ function:

$$F=SV-SL$$

(3)

$$SV=\frac{\sum v_{ij}}{v_{max}}$$

(4)

$$SL=\frac{L-L_{min}}{{L_{max}-L}_{min}}$$

(5)

where $SV$ and $SL$ denote the agent’s scores of the ecological compactness and suitability of the currently formed pattern; $v_{max}$ denotes the potential maximum total ecological suitability value given the land quantity; $L$ denotes the total perimeter of patches formed by the agents; and $L_{max}$ and $L_{min}$ denote the perimeters of a cell and individual agents (i.e., perimeter per cell × cell count), respectively, given the land quantity. If the $F$ value converges to a stable level, the model iteration is terminated to obtain the best model for the environmental protection area.

3.2. Basic Simulations of the CA Model for UUS Growth

The urban CA, in general, can be described using the following equation [36,40,55,56,84], and the meta-CA for UUS growth simulation also obeys such a definition:

$$S_{ij}^{t+1}=f(S_{ij}^t,{\Omega }_{ij}^t,Con,N)$$

(6)

where $S_{ij}^{t+1}$ and $S_{ij}^t$ denote the states of cell $(i,j)$ at the moment of $t+1$ and $t$, respectively; ${\Omega }_{ij}^t$ denotes the neighborhood effect; $Con$ denotes a constraint function of urban expansion; $N$ denotes the total number of cells in the whole cell space; $f$ denotes the transition rules; and $t$ denotes the number of iterations of the CA simulation.

The transition rules determine the state changes of cells within urban CAs, and many studies have applied transition probabilities to represent the chance of cellular state changes and established corresponding transition rules [36,85,86]. The variation in cell state is thought to be determined by the local transition probabilities influenced by spatial drivers [56,83]. According to the process and underlying mechanism of UUS development [5,6,16,26] and the spatial driving factors considered in previous urban simulation studies [63,83,87,88], distance to train station (a), distance to metro station (b), distance to express road (c), distance to main road (d), distance to subsidiary road (e), distance to by-pass road (f), slope (g), elevation (h), average rainfall (i), average temperature (j), population (k), and GDP (l) were used as spatial driving factors of the UUS simulations.

The contribution of various spatial drivers to the local transition probability is determined using several methods, such as expert experience, spatial statistics, and artificial intelligence [56,87,89]. As a tree-based algorithm, an RF algorithm demonstrates advantages in dealing with interactions among predictors, overcoming overfitting problems, and handling high-dimensional data sets [90,91]. This study uses the importance of the RF as well as its associated Gini feature to select the factors that can influence the changes in UUS. By comparing the UUS utilization patterns, the constant pixels are 0, and the changed pixels are 1. As the binary values of the pixels are set as the classification criteria, the 12 spatial driving factors are normalized and arranged in a set of random trees to evaluate the classification results. Thus, the variable with higher Gini feature importance has a greater impact on UUS utilization changes and vice versa. The Gini feature importance was calculated as follows [41]:

$$p\left(X\right)=\frac{1}{N_T}{\sum }_{t=1(t\in T)}^{N_T}(E{R}_{tx}-E_{tx})$$

(7)

where $p\left(X\right)$ denotes the Gini importance of variable $X$; $N_T$ denotes the number of trees; $T$ denotes all trees in the forest; $E_{tx}$ denotes the out-of-bag error of tree $t$ before variable $X$ permutation; and $E{R}_{tx}$ denotes the out-of-bag error of tree $t$ after variable $X$ permutation. The factorial selection procedure was implemented in R studio using an RF package.

Based on the above analysis, the conversion probability calculation formula is established as follows:

$$P_{l,ij}={\sum }_{k=1}^{12}{p}_k\left(X\right)x_{k,ij}$$

(8)

where $P_{l,ij}$ denotes the local development potential of cell $(i,j)$; $x_k$ denotes the $k$th spatial variable; $k$ denotes the spatial driving factors; $i$= a, b, c, d, e, f, g, h, i, j, k, and l correspond to distance to train station (a), distance to metro station (b), distance to express road (c), distance to main road (d), distance to subsidiary road (e), distance to by-pass road (f), slope (g), elevation (h), average rainfall (i), average temperature (j), population (k), and GDP (l), respectively; and $p_k\left(X\right)$ denotes the corresponding weight coefficient.

The neighborhood effect is a central component of the urban CA. It embodies the essence of “bottom-up” self-organized evolution in CAs. Generally, the more developed the cells in a neighborhood, the more likely the central cell will be transformed into UUS use. Most studies use Von Neumann or Moore neighborhoods to construct neighborhood rules. A 3 × 3 Moore neighborhood was used in this study to calculate the effect of adjacent developed cells on the development probability of the central cell [43,45,92]:

$$P_{\Omega ,ij}^t=\frac{{\sum }_{3\times 3}con(S_{ij}^t=developed)}{3\times 3}$$

(9)

where $P_{\Omega ,ij}^t$ denotes the transition probability estimated using the neighborhood function and $con(S_{ij}^t=UUSdeveloped)$ is calculated based on the number of urban UUS cells in the 3 × 3 kernel neighborhood.

Besides the influence of various drivers and neighborhood effects, the urbanization process is usually affected by some stochastic contingencies that make the expansion of UUS uncertain. Consequently, stochastic factor mechanisms and development constraints are introduced in the UUS CA model, and the final transition probability that determines the change in the central cellular state is expressed as follows [34,35,56,84]:

$$P_{ij}^t=P_{l,ij}\times Con(\stackrel{̑}{A}\cdot )\times \left[1+(-\mathit{ln}\gamma {)}^{\alpha }\right]$$

(10)

where $P_{ij}^t$ denotes the probability under which $cell(i,j)$ may change its state from non-developed to developed in the $t$-th iteration of the UUS cellular automaton; $Con(\stackrel{̑}{A}\cdot )$ denotes a constraint function; $\left[1+(-\mathit{ln}\gamma {)}^{\alpha }\right]$ denotes a random disturbance factor, wherein $\gamma$ denotes a random value that ranges from 0 to 1; and $\alpha$ denotes a parameter for controlling the degree of random interference that generally ranges from 1 to 10.

3.3. Simulating UUS Growth Using the Patch-Based CA Model

Two complementary procedures (organic and spontaneous growth) of the patch-based CA model were used in this study to simulate three UUS growth types: infilling, edge expansion, and outlying (as illustrated in Figure 1). The AgentLA model was used to generate ecological reserves as spatial constraints for the UUS growth simulation (as illustrated in Figure 3a) before using the patch-based CA model to simulate UUS growth. The patch growing algorithm [61,63,65,66,83,87] used to simulate three UUS growth types was mainly achieved as illustrated in Figure 3b. An operational code of the patch-based CA Model can be downloaded from the following webpage: https://github.com/HPSCIL/Patch-generating_Land_Use_Simulation_Model (accessed on 5 September 2023).

• Suitability-based sorting: The patch growth algorithm ranks the suitability of all cells according to Equation (10). The corresponding subset of cells with the highest suitability is retained according to the expected area of UUS, and the cells with very low suitability are excluded.
• Stochastic seeding: Stochastic seeding is performed on a subset of the retained cells of UUS. A stochastically selected cell is tested for stochastic seeding; the selected cell is used as the seed cell if its suitability is greater than a random number uniformly distributed in the interval of [0, 1], and the cell is accepted as a seed to start a new patch. Otherwise, the stochastic seeding test continues to iterate until a seed cell appears. For the two complementary growth procedures of UUS, we treat the seed cells connected to the built-up area as organically grown cells and the seed cells not connected to it as spontaneously grown cells.
• Moving-window scanning process: A 3 × 3 cell is chosen as the scanning window size, with the seed cell as the center and the other cells within the window as candidates for the self-growing of new patches. The candidate cells determine the next seed cell according to the stochastic seeding principle, while the previous seed cell transforms into a newly grown patch. Then, the scanning window is moved to the center of the new seed cell and continues stochastic seeding within the window to grow patches until the area reaches the expected area of UUS and the self-growing of patches ends. Notably, when the scanning window moves to the new seed cell while the candidate cells are already within the window, then the patch growth is made more compact by multiplying the isometric parameter in the interval of [0, 2] to narrow or amplify its adaptation.

3.4. Validation of Model Accuracy

To validate the accuracy of the data of the UUS growth simulation of this case study, and to enable their use for UUS growth prediction under different scenarios for future years, three parameters—the OA, Kappa coefficient, and figure of merit (FoM)—were selected for the modeling accuracy assessment [41,56,93].

OA is evaluated based on whether each cell has been correctly assigned a value of 0 or 1 [56,88,94], as calculated using the following formula:

$$OA=\frac{TP+TN}{TP+TN+FP+FN}$$

(11)

where True Positive (TP) denotes that the cell state is correctly classified as a positive (1) sample—that is, the prediction is assigned a value of 1 and the cell is actually assigned a value of 1 as well; True Negative (TN) denotes that the cell state is correctly classified as a negative (0) sample—that is, the prediction is assigned a value of 0 and the cell is actually assigned a value of 0 as well; False Positive (FP) denotes that the cell state is incorrectly classified as a positive (1) sample—that is, the prediction is assigned a value of 1 but the cell is actually assigned a value of 0; False Negative (FN) denotes that the cell state is incorrectly classified as a negative (0) sample—that is, the prediction is assigned a value of 0 but the cell is actually assigned a value of 1.

The Kappa coefficient is used as a metric of the CA model to measure the agreement between the simulation and the reference. The Kappa coefficient varies from 0 to 1, with larger Kappa coefficients representing higher overall simulation accuracy; 0.75 is often considered to be the threshold for high simulation accuracy [93,95] and is calculated as follows:

$$K=\frac{P_0+P_e}{1-P_e}$$

(12)

$$P_e=\frac{(TP+FN)\times (TP+FP)+(FN+TN)TN\times (TN+FP)}{N^2}$$

(13)

where $P_0$ denotes the OA as the sum of the number of correctly classified samples in each category divided by the total number of samples; $P_e$ denotes the sum of the product of the number of true samples of each type and the number of predicted samples, divided by the square of the total number of samples (see Equation (13) for the calculation); and $N$ denotes the number of metacell samples.

Besides the OA and the Kappa coefficient, the FoM was used to assess the accuracy of the simulated changes. Larger FoM values demonstrate more accurate simulations. In theory, there is no fixed range of FoM values, and most FoMs reported in previous studies are below 0.3 [41,61,96]. The FoM is calculated as follows:

$$FoM=B/(A+B+C+D)$$

(14)

where $B$ denotes the pixel numbers with an accurate prediction of UUS growth; $A$ denotes the numbers that had an observed change but were predicted as a persistence of UUS growth; $C$ denotes the numbers that had an observed persistence but were predicted as a change in UUS growth; and $D$ denotes the numbers that had an observed change but were predicted as changing to an incorrect category of UUS growth.

4. Information and Data of the Study Area

4.1. Study Area

The study area of this paper is Tianfu New District, which is a national-level new area in Chengdu City, Sichuan Province, China. Among them, Tianfu New District includes the Tianfu New District Directly Administered Area, the Chengdu High-tech District, Shuangliu District, Longquanyi County, Xinjin County, Jianyang City (part of the area), Pengshan County, and Renshou County (part of the area). The specific location of the area is illustrated in Figure 4. Chengdu City is located in the central part of Sichuan Province and is its capital. It is an important city in western China. Chengdu City was the first demonstration area for park city construction in China. The government of Tianfu New District, which is a part of Chengdu City, endeavored to build a park city that fully practiced the new development concept. In Tianfu New District, the total planning area is 1578 square kilometers, the urban population will be controlled at 5.8–6.3 million people by 2030, and the various ecological land is about 928 square kilometers, including about 90 square kilometers of green space; the green area per capita can reach 15 square meters.

4.2. Experimental Data

The data of Tianfu New District in this study were mainly obtained from different local government departments’ reports in Chengdu, including the “Territorial and Spatial Master Plan of Chengdu Directly Administered Area of Tianfu New District (2019–2035)”, the “Master Plan of Sichuan Tianfu New District (2010–2030)”, the “Chengdu Statistical Yearbook 2021”, “Park City Construction in Chengdu Directly Administered Area of Tianfu New District of Sichuan (2021–2025)”, and other reports and related atlases. In addition, web crawler and remote sensing data were included, which are categorized and organized as follows:

• Data on urban planning land, such as study area boundaries, the area of building land and its spatial distribution, and land-use types and their spatial distribution;
• Data on UUS planning, such as the area and spatial distribution of UUS development and the visionary planning situation of UUS;
• Geographic information-related data, such as slope, elevation, aspect, average rainfall, and average temperature data;
• Transportation-related data information, such as the spatial distribution of subway stations and lines, the location of train stations, the spatial distribution of highway lines, and the spatial layout of city roads.

Land-use data for the study area were obtained from the Resource and Environment Science and Data Center in China (http://www.resdc.cn (accessed on 18 July 2023)), which is classified into six main categories: arable land, forest land, grassland, water bodies, built-up areas, and unused land. The vector road network data were obtained from OpenStreetMap (OSM). Population and GDP data were obtained from the Chengdu Statistical Yearbook 2021. Elevation data were obtained from Google Earth using BIGEMAP software, elevation maps (dem files) were downloaded, and slope data were obtained from a DEM using ArcGIS software. The climate-related data were obtained from the Geographic Remote Sensing Ecology Network (http://www.gisrs.cn (accessed on 18 July 2023)). To facilitate subsequent data processing, we uniformly converted the various types of data in

Table 2. List of data sources for UUS growth simulation.

Categories	Data Name	Data Source	Year	Spatial Resolution
Land use	Land-use data	Resource and Environment Science and Data Center in China	2010–2020	300 m × 300 m
Topographical factors	Elevation	BIGEMAP Map Downloader	2010–2020	100 m × 100 m
	Slope	BIGEMAP Map Downloader	2010	100 m × 100 m
Socioeconomic factors	Population	Chengdu Statistical Yearbook 2021	2010–2020	1 km × 1 km
	GDP	Chengdu Statistical Yearbook 2021	2010–2020	1 km × 1 km
Traffic network factors	Distance to express road	OpenStreetMap (OSM)	2010–2020	300 m × 300 m
	Distance to main road	OpenStreetMap (OSM)	2010–2020	300 m × 300 m
	Distance to subsidiary road	OpenStreetMap (OSM)	2010–2020	300 m × 300 m
	Distance to by-pass road	OpenStreetMap (OSM)	2010–2020	300 m × 300 m
Transportation infrastructure factors	Distance to metro station	OpenStreetMap (OSM)	2010–2020	300 m × 300 m
	Distance to train station	OpenStreetMap (OSM)	2010–2020	300 m × 300 m
Constraint data of UUS development	Urban construction land boundary	Master Plan of Sichuan Tianfu New District (2010–2030)	2010–2030	Vector polygons
	Mountains and water systems	Master Plan of Sichuan Tianfu New District (2010–2030)	2010–2030	Vector polygons
	Urban BGS protection areas	Master Plan of Sichuan Tianfu New District (2010–2030), Field survey and map street view	2010–2030	Vector polygons
Climate factors	Average temperature	Geographic Remote Sensing Ecological Network	2010–2020	100 m × 100 m
	Average rainfall

into vector data. The spatial coordinate system and resolution (300 m), and the corresponding Euclidean distance and point density, were calculated from the road network and point-of-interest data, and the spatial driving factors for the study of regional UUS growth simulation were constructed as illustrated in Figure 5.

5. Implementations and Results

5.1. Setting of Model Parameters

The sample data set was divided into training data (from 2010 to 2015) and validation data (from 2016 to 2020), as the validation data were used to evaluate the prediction accuracy to ensure the accuracy and stability of the model. The overall accuracy (OA) and Kappa coefficient values (as illustrated in Equations (11) and (12)) increase slightly as the neighborhood size increases from 3 × 3 to 7 × 7. Moreover, using too large of a neighborhood (e.g., 9 × 9) may reduce the sensitivity of the CA model neighborhood interaction (and, thus, the accuracy of the simulation). Therefore, based on better model performance, this CA model set the neighborhood size to 7 × 7. The number of sub-trees of the trained RF model is 100, the maximum depth of the sub-tree is 10, and the maximum number of features considered for each split of the sub-tree node is 5. In addition, the results of calculating the importance of the RF model for the spatial driver are illustrated in Figure 6, which demonstrates that the GDP, population, distance to a metro station, distance to a main road, and slope are significant for UUS development. This is mainly because socioeconomic factors (S1) are the most important factors in UUS development.

Constructing and utilizing UUS requires a large amount of capital, which is closely related to the city’s GDP; however, the larger the population, the larger the scale of UUS development required. The transportation network distance (S2) and transportation infrastructure (S3) are influential because the transportation network plays an important role in determining the shape of UUS, and the intensity of UUS development is higher in areas with convenient transportation networks. Topographic factors (S4) also play an important role in UUS development in the Tianfu New District, which is rich in topographic features, including mountains, lakes, hills, and plains. Therefore, owing to the large topographic variability, the average temperature and average rainfall factors (S5) play an influential role when considering the potential floodplains caused by climate change. The performance of the AgentLA model is sensitive to the $w_v$ setting. The model F value is highest when the $w_v$ value is 0.7–0.8 [83]; thus, the $w_v$ value is set to 0.75 for the AgentLA model.

5.2. Results of UUS Growth Simulation

The average values of the OA, FoM, and Kappa coefficient were calculated for the validation data from 2016 to 2020, according to the above formulae. The OA was calculated as 90.42%, the Kappa coefficient as 0.8424, and the FoM as 0.638426. The calculated results of the Kappa coefficient and FoM were higher than those reported in the previous related literature [41,61,93,95,96]; therefore, the constructed UUS growth simulation model has good simulation accuracy and can be used to simulate the growth of UUS for future years.

Using the UUS utilization in the study area in 2010 as the initial state and applying the growth simulation model of UUS constructed in this study, the growth simulation results of the UUS utilization in the study area in 2020 were obtained, as illustrated in Figure 7b, and compared with the UUS planning in the study area in 2020. The data sources were the Master Plan of Chengdu Tianfu New District, Sichuan Province, 2010–2030; the Logistic-CA model [55]; the artificial neural network (ANN)-CA model [97]; the and RF-patch-CA (unconstrained) model for comparison, as illustrated in Figure 7a,c–e, respectively.

Figure 7 demonstrates that the UUS growth simulation results of the model constructed in this study are closer to the set planning results than are those of the other three CA models; locally, the other three CA models demonstrate the obvious phenomenon of UUS central metacell clustering, while the model constructed in this study effectively combats this phenomenon. To further verify the accuracy of the UUS growth simulation model constructed in this study, the Logistic-CA model, ANN-CA model, and RF-patch-CA (unconstrained) model were calculated at the metacell level, and the OA, FoM, and Kappa coefficient results of the model constructed in this study are illustrated in

Table 3. Precision comparison of UUS utilization growth simulation results of different models.

Validation of Model Accuracy	The Model Constructed in This Study	RF-Patch-CA (Unconstrained)	ANN-CA	Logistic-CA
OA	0.9042	0.7647	0.7743	0.7658
Kappa coefficient	0.7424	0.6150	0.6303	0.6161
FoM	0.3384	0.1164	0.1060	0.1178

.

Table 3 demonstrates that the OA, FoM, and Kappa coefficient values of the model constructed in this paper are higher than those of the other three models, in which the Kappa coefficient improved by 0.1273, 0.1121, and 0.1263, and the FoM increased by 0.2220, 0.2325, and 0.2206, respectively. Therefore, these results further show that the model constructed in this study simulates the growth change pattern of UUS with greater accuracy.

6. Discussion and Analysis

6.1. Unconstrained by Ecological Constraints of UUS Growth Simulation

Using the UUS growth simulation model proposed in this study, the simulation of UUS growth without ecological constraints in 2020, 2025, and 2030 considering only the prohibited construction areas, such as mountains and water systems, was calculated (as illustrated in Figure 8), and the change in the area during the UUS growth simulation in each region was calculated (as illustrated in

Table 4. Comparison of UUS use area changes without ecological constraints.

Name of Area	Area Change in UUS Use (km2)		Unconstrained by Ecological Constraints of Area Change in UUS Use (km2)
	2010–2015	2015–2020	2020–2025	2025–2030
Pengshan County	6.6132	9.6533	6.7464	1.3467
Renshou County	8.2665	14.0530	7.4398	6.6132
Xinjin County	19.8395	19.0128	9.0931	3.3066
Shuangliu District	86.3844	83.3909	61.4105	19.2829
Tianfu New District Directly Administered Area	16.9462	67.9414	72.6906	87.9161
Chengdu High-tech District	18.1862	20.4799	3.7994	1.4447
Longquanyi County	41.3323	26.4526	9.5714	5.8122
Jianyang City	1.6533	1.3467	1.1599	1.8796

).

Table 4 and Figure 8 demonstrate that the area of UUS utilization in 2020, 2025, and 2030 will be 637.47 km², 809.38 km², and 936.98 km², respectively. The area of UUS utilization will increase by 299.51 km² from 2020 to 2030, with an average annual growth rate of 4.70%. Further analysis of the area changes for UUS utilization in each region demonstrates that different regions have different characteristics of change. Their changes can be roughly divided into three categories: (1) the scale of UUS development remains unchanged, with Jianyang City as a representative region; (2) the scale of UUS development grows rapidly and continuously, representing Shuangliu District, the Tianfu New District Directly Administered Area, and Longquanyi County; (3) the scale of UUS development appears to grow at the beginning, but the growth rate of these areas gradually decreases after 2020, representing Pengshan County, Renshou County, Chengdu High-tech District, and Xinjin County.

6.2. Constrained by Ecological Constraints of UUS Growth Simulation

The simulation of UUS growth with ecological constraints considers the UUS utilization constraints of urban BGS protection areas. Using the UUS growth simulation model proposed in this study, the simulation of UUS growth with ecological constraints in 2020, 2025, and 2030, considering the urban BGS protection areas, was calculated (as illustrated in Figure 9), and the area change during the UUS growth simulation in each region was calculated (as illustrated in

Table 5. Comparison of UUS use area changes under ecological constraints.

Name of Area	Area Change in UUS Use (km2)		Constrained by Ecological Constraints of Area Change in UUS Use (km2)
	2010–2015	2015–2020	2020–2025	2025–2030
Pengshan County	6.6132	9.6533	3.2665	0.7533
Renshou County	8.2665	14.0530	3.0931	2.6934
Xinjin County	19.8395	19.0128	5.3066	1.9197
Shuangliu District	86.3844	83.3909	16.4708	5.6101
Tianfu New District Directly Administered Area	16.9462	67.9414	23.2354	27.2920
Chengdu High-tech District	18.1862	20.4799	3.1332	1.0128
Longquanyi County	41.3323	26.4526	9.0931	4.7994
Jianyang City	1.6533	1.3467	0.9599	1.0796

).

Table 5 and Figure 9 demonstrate that the area of UUS utilization in 2020, 2025, and 2030 will be 637.47 km², 702.03 km², and 747.19 km², respectively. Compared to the case without ecological constraints, the scale of UUS development continues to grow, but its growth rate is significantly lower. The area of UUS utilization will increase by 109.72 km² from 2020 to 2030, with an average annual growth rate of 1.72%, which is significantly less than that in the unconstrained case. Further analysis of the area changes for UUS utilization in each region demonstrates that the area of UUS utilization in Shuangliu District, the Tianfu New District Directly Administered Area, and Longquanyi County will continue to grow from 2015 to 2030. However, owing to the development constraints of urban BGSs’ protection areas on UUS, its growth rate will gradually decrease and be significantly smaller than the growth rate under the “no ecological constraints” scenario.

6.3. Comparative Analysis with or without Ecological Constraints

To compare the difference of the UUS growth simulation process with and without ecological constraint scenarios, we compare and analyze the spatial distribution of ecological constraints on urban space in the study area (as illustrated in Figure 10b), the UUS growth simulation without ecological constraints (as illustrated in Figure 10a), and the UUS growth simulation with ecological constraints (as illustrated in Figure 10c), and the changes in the UUS utilization area of UUS with and without ecological constraints are compared (as illustrated in Figure 11).

By comparing the results of the UUS growth simulations with and without ecological constraint scenarios in Figure 10 and Figure 11, we find that urban ecological space constraints are important in influencing the UUS development scale. In both cases, the scale of future UUS development is increasing, but the magnitude of the increase is significantly different. Under ecological constraints, the scale of future UUS development is much smaller than that under unconstrained development. Most of the urban BGS protection areas are located in the Tianfu New District Directly Administered Area, wherein the growth rate of UUS development is significantly slower than that in the unconstrained scenario, owing to the constraint of urban BGS protection areas. Owing to the fragmented distribution of urban BGS protection areas, many fragmented area patches (utilization areas of UUS) have been formed in the Tianfu New District Directly Administered Area. Without ecological constraints, the growth of the UUS development scale in the region will demonstrate a disorderly expansion due to a lack of constraints in the form of urban BGS protection areas. Therefore, UUS development under ecological constraints will be more sustainable and can protect urban BGSs on the ground more effectively, thereby promoting urban parks in Tianfu New District and constructing a low-carbon urban area.

7. Conclusions

This study combined the AgentLA model with the patch-CA model to construct a UUS growth simulation model. By taking Chengdu Tianfu New District as an example, the development process of the UUS scale from 2015 to 2030 was simulated, and two scenarios with and without ecological constraints were provided to predict the simulation results of UUS growth from 2020 to 2030. The simulation results with Chengdu Tianfu New District demonstrate that the model constructed in this study is suitable for simulating the process of changing the scale of UUS development. According to the relative importance measure of the RF algorithm, the GDP, population, distance to the subway station, distance to the main road, and slope have a large impact on the probability of conversion to UUS development. In the case of UUS development without ecological restrictions, the scale of UUS development in Tianfu New District will demonstrate a rapid and large-scale growth in the future; specifically, the Tianfu New District Directly Administered Area is growing faster than other areas. Simultaneously, under the consideration of ecological space constraints, although the scale of UUS development will continue to grow, UUS expansion is characterized by a slow and sustainable growth rate.

Overall, this study explains how the sustainable elements of ecological space can be incorporated into the proposed model, which combines the AgentLA model with the patch-CA model and helps to better simulate the complex dynamics of UUS over time. By comparing the results of different scenarios, we found that UUS planning, considering the urban ecological space constraint, can better implement the concept of a “park city” in the Tianfu New District, preserve more urban BGSs, and optimize the urban spatial structure while the UUS growth simulation combined with spatial drivers can also provide a theoretical basis for future UUS planning and urban planners. More importantly, this study identifies the spatial drivers of UUS growth simulations and their weights based on the RF algorithm, which provides a reference for the study of growth simulations of other UUS developments. On this basis, the patch-based CA model was constructed under ecological constraints and has a high prediction accuracy. The model also demonstrates a high portability for studies on the simulation of the growth of UUS in other cities.

This work represents an ecologically constrained approach to UUS growth simulation; however, limitations are inevitable. This study attempts to construct a UUS growth simulation process in which urban ecological space is effectively protected, but only urban BGSs, mainly large urban parks (including those with water bodies), are considered in the specific case and analysis, while some smaller community-based parks are not considered. A further and more comprehensive refinement of urban ecological space could better reflect the model’s orientation to UUS development. Another limitation is that the development and utilization of UUS can increase the area of ecological space in urban areas by re-constructing it on saved land; however, this aspect has not been considered in constructing our model. Simultaneously, the process of UUS development and utilization can create three-dimensional layout forms, such as urban atriums or sunken spaces, and the creation of such three-dimensional layout forms will include green vegetation planting, which positively impacts urban ecological space and the construction of livable cities.