Simulating urban growth by coupling macro-processes and micro-dynamics: a case study on Wuhan, China

ABSTRACT

The urban form influences the quality of urban functions and is strongly correlated with the sustaining capabilities of urban development. However, in the context of rapid urbanization, unreasonable land expansion as a universal phenomenon poses a great challenge for urban management. Notably, the urban expansion process is self-organizing, and the evolving macroscopic pattern can be used to predict microscopic behavioral characteristics. Therefore, the analysis of macro- and micro-interactions can provide new ideas for urban modeling. Traditional geographic cellular automata (CA) models often have poor morphological reproducibility, and the few models that combine top-down and bottom-up CA use strict coupling constraints, resulting in inadequate self-organizing natural expressions and poor precision performances. In this study, we proposed a new land growth simulation model based on a soft constraint mechanism that couples micro-dynamics with macro-processes. Specifically, a geographic micro-process model (GMP) based on the meta-process accumulation concept was applied to capture the evolution characteristics of the macro-urban form and spatially deduce the future urban intensity gradient. The soft coupling between the macro and micro levels of the model was supported by a punishment mechanism that was developed for this study. A specially designed index, the morphology similarity (MS) index, was developed to evaluate and understand the heterogeneity of the simulated and real urban forms from a micro-perspective. The model was applied to Wuhan, the largest city in central China, to demonstrate that the proposed model has a high simulation accuracy [with a Kappa value of 0.8506 and a figure-of-merit (FoM) value of 0.3034 in the optimal parameter combination] and imitative ability [maximum sensitivity (MS) value of 0.01341 in the optimal parameter combination vs. MS value of 0.01336 in the true scenario]. The evaluation system developed in this study also demonstrated the high robustness and reliability of the future multi-scenario simulation conducted in this work.

KEYWORDS:

• Land growth simulation model
• geographic micro-process model
• urban intensity gradient
• urban modeling

1. Introduction

Urban form is a comprehensive result of natural and human activities (Glaeser and Kahn 2003). Evolving from the original settlements, each city has formed a unique spatial pattern by relying on the economic, social, political, and other incentives of human civilization to expand outward (Abaya Gomez 2015; Kropf 2018). This type of dynamic development carries urban functions and reshapes the urban landscape (Zhang and Wang 2006). Urban spatial expansion also dominates urban forms (Zhang, Ghosh, and Park 2023). However, the unprecedented speed of urbanization is causing serious urban problems, such as reduced biodiversity, damaged arable land, declining air quality, threatened water resources, polluted soil, land subsidence, and heavy traffic (Chung, Frank, and Pokhrel 2021; Jaeger and Schwick 2014; Li et al. 2016; Mahtta et al. 2022; Xu et al. 2019). Many contemporary studies have confirmed that a compact urban design is capable of reducing carbon emissions (Ye et al. 2015), enhancing facility accessibility (Hidalgo and Gutiérrez 2013), shortening travel distance (Xu et al. 2020), and mitigating the urban heat island effect (Li et al. 2020), so as to support the well-being of urban inhabitants. Hence, urban planners should realize the significance of the urban form and develop reasonable plans for the future development of the city accordingly.

At present, urban modeling is primarily based on the self-organization aspect of urban evolution. Based on this self-organizing nature, the relationship between macroscopic and microscopic environments is that of affiliation and interaction (Artime and Domenico 2022). Notably, macroscopic patterning emerges from dynamic interactions at the microscopic level, and deconstructing the emergent macroscopic structure can reconstruct microscopic dynamic processes (Barnett and Seth 2023). Thus, it is possible to use top-down and bottom-up interaction modes to realize the self-organization of urban evolution. A cellular automata (CA) model is commonly used to realize the spatiotemporal evolution of complex systems, based on spatial interaction and temporal causality. Combined with a geographic information system (GIS), this underpins our understanding of the evolution of spatial and temporal patterns in urban areas (Stevens, Dragicevic, and Rothley 2007). This spatially explicit modeling method has been extensively applied in the field of global urban expansion (Feng et al. 2020; Isinkaralar and Varol 2023; Isinkaralar, Varol, and Yilmaz 2022), urban growth boundary identification (Wang et al. 2020), urban development zone design (Liang et al. 2021), land use pattern optimization (Li et al. 2022), and mitigation of human-environment conflicts (Gong et al. 2009), thus, providing strong technical support for land use policy and land management decision-making. At present, innovation has made substantial progress at the method level (Feng et al. 2011; Liang et al. 2021; Liu et al. 2017); this is especially true for the precision improvement in transformation mechanisms, configuration of neighborhood relations (Liang et al. 2020), multi-scale strategy (Chen et al. 2020) and dimension enhancement (Chen and Feng 2022). These studies strove hard to reproduce urban land change and predict future landscape scenarios; however, they ignored the heterogeneity of land demand within cities (Lin et al. 2023).

The heterogeneity of urban spatial demand stems from the game theory behind spatial supply and demand (Bednar-Friedl, Koland, and Steininger 2011). Agglomeration and dispersion forces within each city (Xu et al. 2019) assemble the physical building space elements (Lei et al. 2021), such as impervious surfaces, traffic roads, and social space elements (food service and medical service facilities and shopping services), into a more compact structure near the core urban area. Thus, these forces generally cause resource scarcity in suburban areas. The relationship between distance decay and urban form is presented by way of a mutual entwine of sophistication, which further drives the crop layout, travel cost, and decision of enterprises in location selection, along with other aspects within cities (Glaeser 2007; Tobler 1970). Therefore, the phenomenon of spatial distribution imbalance has continuously pushed scholars to develop more intuitive and accurate approaches for quantifying land demand during urban growth. In general, land density decay is characterized by a series of mathematical solutions, including the negative exponential, inverse power, Tanner, Inverse-S-shaped, and Gaussian-based functions and geographic micro-process model (Batty and Sik Kim 1992; Chen and Feng 2022; Clark 1951; Jiao 2015; Jiao et al. 2021; Yang et al. 2022), providing insights into the law of urban land growth in the macro process. The geographic micro-process (GMP) model regards urban spatial expansion as the accumulation of geographical meta-processes, and the probability of all newer urban developments in vacant buildable land refers to the distance attenuation law. With few adaptable parameters and ease of implementation, it has great potential for calculating macroscopic land density gradients.

Considering the heterogeneity of spatial demand as an imperative clue to solve the puzzle of urban self-organization and integrating this approach into the land growth modeling framework can better demonstrate the characteristics of urban self-organization, guide the simulation toward a more similar land form, and provide an efficient reference for the optimization of land-use spatial patterns and the selection of future urban development paths. Wang et al. (2019) combined an inverse S-shaped function with the CA model to simulate urban land in a top-down and bottom-up coupled framework using a series of divided buffer zones. By incorporating spatiotemporal Gaussian-based functions into the CA model, Yang et al. (2023) designed an innovative method for achieving a spatiotemporally explicit representation of land demand. However, in the modeling process, the connection between the macro- and micro-levels often requires a tough coupling strategy, which does not reflect the self-organization of the urban system well and may result in an unexpected decrease in the model accuracy. Owing to the distance attenuation effect, there is an unavoidable phenomenon wherein it becomes more difficult to preserve precision as the simulation progresses outward from the urban center and sub-centers. Forcibly restricting the number of new land cells, combined with the gradient analysis, might have led to a decline in the accuracy of the results. Therefore, this study examines a new method that can naturally capture the complexity generated in simple and logically coherent rules of city planning, so that spatial form can emerge from presupposition models (Li et al. 2017; Liu et al. 2022).

In the related accuracy evaluation system, the consistency of the classification results of urban and non-urban land receives excessive attention, and the degree of morphological similarity tends to be forgotten. This paradigm implies that the morphological reproducibility of the proposed model, as well as the correctness of future landscape scenarios, cannot be determined effectively. At present, the landform measurement method has shifted from simply measuring macro features, such as land increment, growth rate, and density, to calculating more detailed microfields, including the degrees of spread and fragmentation and complexity (Jiao et al. 2018; Liu et al. 2021). Therefore, the use of both macro-quantitative characteristics and emerging micro-measurement methods can be expected to evolve into a general trend. Angel et al. (2020) investigated the geometric attributes of two-dimensional urban footprint. In particular, the five subcategories of shape compactness of the urban footprint in their view include contiguity, proximity, cohesion, exchange, and fullness; this perspective provided the scholars in the field novel ideas for the further exploration of morphological evaluation systems.

To address the aforementioned problems, in this study, we focused on the following questions: (1) Can we tailor the GMP model constructed from the perspective of the urban natural growth process to the geographical CA model effectively? (2) Can the model be a breakthrough in the traditional threshold termination mode in the gradient constraint during the coupling process? (3) Is the model capable of verifying the form-imitation ability of the proposed model and measuring the correctness of the landform in the simulation results of multiple future scenarios. (4) Is the morphological evaluation system established by combining macro- and micro-features robust and effective?

Therefore, this study proposes a new explicit spatiotemporal urban growth model. First, the GMP model was used to characterize the urban intensity gradients. Second, based on the relationship between the expected land growth and the demand in the gradient constraint, a new soft-constraint method was developed to bridge the macro-processes and micro-dynamics. Furthermore, based on the land evolution trend in its macro-processes, we designed different scenarios of future land growth, i.e. compact, stable, and sprawl growths. Finally, we presented an index to measure the morphological reproduction ability of the model and offered a comprehensive morphological evaluation index system. Our study framework was applied to Wuhan (China), to explore and verify the accuracy and validity of the model; notably, our study model can provide technical support for land decision-making and management in a larger study area in the future.

2. Materials

2.1. Study area

Wuhan, the capital of Hubei Province, is the most flourishing city in inland China. It serves as a key point of access to the upper and lower reaches of the Yangtze River, which is one of the most important river basins in China. Despite the impact of the COVID-19 epidemic, the total gross domestic product (GDP) of the country has remained in the top 10 countries for years. The seventh census (National Bureau of Statistics of China 2021b) states that Wuhan has a permanent population of 12.327 million, indicating that it has become China’s eighth megacity. To meet the demands of rapid and vigorous development, the built-up area of Wuhan has more than quadrupled over the past two decades, from 210 km² in 2000 to 885 km² in 2020 (National Bureau of Statistics of China 2021a). However, urban problems, such as farmland loss, ecological deterioration, and disorderly expansion, are impacting the well-being of the dwellers, hindering the sustainable development of cities, and threatening the stability of the Yangtze River Economic Belt. The local government has launched several plans to curb the unhealthy development of urban land in the region; to promote the effective implementation of these countermeasures in the region, we selected Wuhan as the study area and attempted to capture the characteristics of land dynamic change and reproduce the growth of its urban land cover, while preparing the model for future large-scale simulations (Figure 1).

Figure 1. Geographic location of the study area (Wuhan, Hubei Province, China).

2.2. Data and processing

The urban land data for 1995, 2000, 2005, 2010, 2015, and 2020 were extracted from the annual China Land Cover Dataset (CLCD) released by the Wuhan University (Jie and Xin 2021; URL: https://essd.copernicus.org/articles/13/3907/2021/). This dataset, which had an overall accuracy of 79.31% (based on 5463 visually interpreted samples), originally contained the land cover information for China (resolution of 30 m) for the timeseries of 1990–2019; the timeseries was updated to 2021 and is widely applied in the field of land science (Wang 2022; Xu et al. 2022).

The point-of-interest (POI) data and road network vector datasets, which represent the economic and social vitality of cities, were derived from Amap (https://ditu.amap.com/) and Open Street Map (OSM) (https://www.openstreetmap.org) respectively; notably, for this study, we used the data for 2019. For the analysis of traffic data, we used ArcGIS 10.6 to modify and eliminate inaccuracies. The physiogeographical components and digital elevation model (DEM) data for 2015 were derived from the Shuttle Radar Topography Mission (SRTM) dataset (http://srtm.csi.cgiar.org/srtmdata/), and the slope information was constructed using ArcGIS 10.6. The geographical constraint factor, water area, and scenic area were derived from the CLCD dataset and the area-of-interest (AOI) vector as derived from web crawling.

Previous studies have shown that the timing of data collection should be as close as possible to that of the land-use data, even though referring to different time periods is allowed. Although the driving factors were collected from different time periods, previous studies confirmed that this kind of variance can be permitted as long as the time periods are close to that of the land-use data (Liang et al. 2021).

The city centers and sub-centers (shown in Figure 1) were identified using the contour tree method and POI data (Chen et al. 2017). A widely adopted gradient analysis method was used to cover the study area in the donut-shaped buffer zone, with an interval of 1 km; in total, 47 ring layers were generated (Dong et al. 2019; Jiao 2015; Xu et al. 2019). All the raster data used in this study were resampled to 30-m resolution and the geographic and projection coordinate systems in all the collected data were equated precisely.

3. Methodology

3.1. Method overview

This study uses an innovative bottom-up and top-down coupled soft restraint by capturing the demand relationship of macro land growth in its microevolution, to reproduce and predict land expansion patterns accurately (Figure 2). Based on the transition analysis strategy (TAS) (Liang et al. 2021), first, we applied a machine learning method to identify the relationship between land change and the driving factors for 2000 and 2015; we considered several socioeconomic and physical geographical factors to explore the suitability of land development in the study area. Second, the GMP model was utilized to characterize the land density distribution and corresponding urban intensity gradients for 1995, 2000, 2005, 2010, 2015, and 2020. Third, the expected value of micro-growth, combined with the actual growth value of each iteration, was fused into the geographical CA to realize a punishing mechanism through the activation of the nonlinear function, and the accuracy was verified by comparing the simulation results and land data for 2020. Furthermore, beginning from the historical floating trend of the GMP parameter, we defined three scenarios of future land growth: compact, sprawl, and stable development. Finally, the solved GMP function was introduced into the model to complete the multi-scenario simulation tasks for 2025 and 2030, and the reliability of the study was proved according to the designed urban-form evaluation framework.

Figure 2. Framework of the proposed urban growth model.

3.2. Micro-dynamics modeling

Previous studies demonstrate that the CA model is one of the most powerful tools for simulating recurrent evolutionary processes of various elements within a certain system (Sante et al. 2010; Wang et al. 2020). The process starts with the Game of Life and encompasses several fields, such as basic mathematics, medicine, transportation, finance, and economics.

Geoscience plays an indispensable role in environmental monitoring, land simulation, ecological protection, and disaster warning (Feng et al. 2011; Li and Yeh 2020). In this context, the cells can be envisioned as units of a dense structure, squares on a chessboard, or even hexagons in a beehive. The cells and its neighbors can transmit information about the state at the next time point; using this approach, various dynamic simulation studies have been realized based on specified evolutionary rules, objective constraints, and the number of neighbors. Utilizing this theoretical basis and considering its four components, namely, development potential, neighborhood effect, constraint, and random jamming, we introduced the following modifications.

3.2.1. Development suitability

Development suitability, which is the foundation of evolutionary rules, is a key step in CA realization. During this study, 3000 newly developed and 3000 non-urban land units were randomly sampled to identify the differences in the annual land data. Furthermore, several experimental analyses have shown that land use change is related to geographical conditions, land development history, and socioeconomic activities (Chen et al. 2020; Glaeser and Kahn 2003). Based on this conclusion and the driving factors selected in previous studies (Liang et al. 2021; Wang et al. 2019; Yang et al. 2023), we constructed 10 raster images that considered the economic, social, and geographical natural elements of the study area as well as the natural factors (elevation, slope, distance to water bodies, and social and economic factors). Notably, we also considered other factors, such as the Euclidean distance to primary, secondary, and tertiary roads, motorway, center and sub-centers, highway toll gates, train stations, and airport access entrances.

Furthermore, we used a random forest (RF) algorithm, which is one of the most popular machine learning methods. The RF model has gained a high reputation, owing to its fast training speed, ability to carry out parallelization and timely feature-correlation detection, and capacity to maintain efficient accuracy. In particular, RF can perform quite well when dealing with time-demanding tasks characterized by high dimensionality. To determine the most suitable hyperparameters for the samples, we applied genetic algorithms (GA) to optimize the model. Using GeatPy 2.7.0, a Genetic and Evolutionary Algorithms toolbox for Python, the objective function was set as the average F1-score (used in statistics to measure the accuracy of a binary model; specifically, the weighted average of accuracy and recall rates) in a 10-fold-cross-check of the training data; our optimization goal was to maximize the result. The ratio of the training set to the test set was 8:2, pre-optimization test-set accuracy was 0.80, and post-optimization accuracy was 0.81; most importantly, we observed an effect that prevented overfitting in the training-set accuracy results, which decreased from 1.00 to 0.98. The development suitability (P_potential) was derived using the equation in Feng et al. (2020) and Wu (2002) [see Equation (1)(1) $P_{potential}=f_{ml}\left(x_1,x_2,x_3\dots x_n\right)$(1) ]:

(1) $P_{potential}=f_{ml}\left(x_1,x_2,x_3\dots x_n\right)$(1)

where $f_{ml}$ is the selected machine-learning method, the variables of x represent the previously mentioned economic, social, and geographical natural factors, and n depends on the total number of factors.

3.2.2. Setting the neighborhood effect

The influence of each cell’s neighbors cannot be ignored; a smaller number of neighbors implies that the transferred instruction tends to be more local, whereas a large number of neighbors implies that the relevant information spreads more widely. In general, Moore neighborhoods with eight surroundings and von Neumann neighborhoods with four surroundings are the prevalent set statuses, although other forms, such as 5 × 5, are often used in published reports. We used the classical 3 × 3 neighborhood, to implement logical signals that could carry duties (Wu 2002).

(2) $P_{neighbor}=\frac{f_{n\times n}}{\left(n^2-1\right)}$(2)

Where P_neighbor refers to the influence of each cell on other cells in its $n\times n$ surrounding area, $f_{n\times n}$ represents the relationship of the transformation condition, and n is strongly connected to the predefined number of neighbors.

3.2.3. Introducing constraints

In reality, various natural conditions and planning decisions reduce the growth of built-up urban areas; this implies that the corresponding numerical space is not eligible to accept cells. Thus, the influence of these constraints should be eliminated, to ensure the rationality of the modeling goal. The equation used to analyze the constraints (P_constrain) is shown below (Wang et al. 2020; White and Engelen 2000):

(3) $P_{constrain}=f_{con}\left(C_1,C_2,C_3\dots C_n\right)$(3)

where $f_{con}$ is the specified constraining rule; $C_1,C_2,C_3\dots C_n$ are the constraining attributes; and n is the total number of constraints. In particular, we selected water and scenic parks as the constrained spaces.

3.2.4. Create random disturbance

We aimed to reproduce a complex urban evolution process; therefore, it was difficult to exclude the uncertain spatial-distribution state that was hidden under definite development laws. In virtual modeling, these uncertainties can be characterized by the addition of stochastic disturbances (P_disturb). We used Equation (4)(4) $P_{disturb}=1+{\left(-\mathrm{l}\mathrm{n}\gamma \right)}^{\theta }$(4) to guarantee the repeatability of our model (Feng et al. 2011; White and Engelen 1993):

(4) $P_{disturb}=1+{\left(-\mathrm{l}\mathrm{n}\gamma \right)}^{\theta }$(4)

where γ and θ denotes the relevant parameters of randomness; the value of γ fluctuates between 0 and 1, and the value of θ was set to 2 by trial and error in this study.

Combining these four aspects, the proposed CA model was constructed using the following paradigm (Wu 2002):

(5) $P=P_{potential}\times P_{neighbor}\times P_{constrain}\times P_{disturb}$(5)

where P is the probability of each cell becoming an urban land cell.

3.3. Macro-process modeling

An extensive body of literature on the spatial distribution patterns of urban geographical elements has confirmed that social interaction, economic activity, and the density of population and land obey the distance attenuation law, which can be fitted well using a power function (Batty and Sik Kim 1992; Clark 1951). The power-shaped falloff curves, which grasp the significance of attenuation characteristics from a concise perspective, are powerful quantitative approaches that can help us unfold the mystery of urban land. Among the various attenuation functions, the GMP (Jiao et al. 2021) model, which can simulate the land-moving currents in a dynamic accumulation process, provides a fresh perspective for weighing urban development units. More specifically, the principle of this model is to view urban spatial growth as a sequential process, wherein urban building blocks gradually sprawl outward with the assistance of policies and economic and social resources. Thus, in this scenario, urban growth can be divided into substantial meta-processes over a fixed time period (Figure 3(a)), wherein urban seeds are allowed to bloom in vacant spots in a power-function manner. In other words, the probability of newly growing urban land in developable vacant areas decreases outward from the center and follows the form of a power function at each time node. According to this idea, at the initial time node, $t_0=\delta t$, the urban land distribution function can be defined as $y_0={\mathrm{r}}^{-\alpha }$, where $\mathrm{r}$ stands for the outward distance from the city center and $\alpha$ is the relative parameter of the distribution function. In the following time phases, i.e. $t=2\delta t,3\delta t,4\delta t\dots n\delta t\left(\mathrm{n}\to \mathrm{\infty }\right)$, the overall land density distributions can be decomposed into two sections: developed and newly-added; this can be expressed as $y_1=y_0+{\mathrm{r}}^{-\alpha }\left(1-y_0\right)$,$\hspace{0.17em}{y}_2=y_1+{\mathrm{r}}^{-\alpha }\left(1-y_1\right)\cdots$ $y_n=y_{\text{break}}n-1+{\mathrm{r}}^{-\alpha }\left(1-y_{n-1}\right)\left(\mathrm{n}\to 0\right)$. Hence, by mathematical induction, the urban land density function can be denoted as follows (Jiao et al. 2021):

Figure 3. Conceptual diagram of the geographic micro-process (GMP) model developed in this study. (a) portrays the accumulation of geographic meta-processes, and (b) portrays the influence of parameter regulation.

(6) $y=1-{\left(1-r^{-\alpha }\right)}^n$(6)

where $\alpha$ is the relevant gradient parameter in the meta-process of urban growth, and n represents the cumulative number.

The manner in which parameter $\alpha$ links with urban morphology is notable (Jiao et al. 2021). To be more specific, when we keep n constant, as shown in (Figure 3(b)) (fixed n = 6), a steeper curve will form, and the value of $\alpha$ will increase (the slope of the curve increases, while adjusting the value of $\alpha$ from 1 to 2.5). This indicates that the spatial agglomeration degree of the study object is getting increasingly stronger (see the green area in Figure 3(b)). Similarly, a lower value indicates a relatively loose shape of the spatial entity. As for the n parameter, an invariant state of $\alpha$ leads to the reverse direction between n and the curve rate. For instance, if the $\alpha$ value is set at 2 and the n value is adjusted from 150 to 500, the curve will portray a more complete form (see the orange area in Figure 3(b)). This indicates that the n parameter can reflect the accumulation effect in the process of land growth.

3.4. Bridging macro and micro levels with soft constraints

We already knew that the growth amount in each of the buffer zones was revealed after applying the macroscopic GMP equation. Thus, whenever the urban seeds exceeded a set threshold in a certain ring, a mask within its range would arise to deter further increments. However, this method, intended to maintain the growth number (with rigid control), is often not an appropriate approach. Thus, we endeavored to find a smoother coupling method, while considering two concerns: (1) How can the model satisfy the need to better reflect the self-organization of urban growth? (2) How can we ensure and improve the model’s convergence ability? Therefore, we proposed an innovative solution: we added a penalty term to bridge the gap between the macro and micro components.

Specifically, after each iteration of a certain number of urban seeds, our model recalculated the quantity that increased at the on-going stage of each buffer zone and determined the demand of each buffer zone through a fractional relationship between the expected growth and current amounts. To avoid deficiencies in the traditional linear function, e.g. the inability to respond flexibly and failure to converge, we calculated the demand coefficient value in the sigmoid function. The classical nonlinear function, which featured a unique shape, returned smoothly for the given independent variables. The penalty term can be defined as follows:

(7) $\left\{\begin{array}{c}{Q}_{\left(i,current\right)}=Q_{\left(i,initial\right)}+N\\ x_i=\frac{Q_{\left(i,\mathrm{e}\mathrm{x}\mathrm{p}ect\right)}-Q_{\left(i,current\right)}}{Q_{\left(i,\mathrm{e}\mathrm{x}\mathrm{p}ect\right)}}\\ y_i=Co{n}_{(DT>D)}\left(\frac{1-c}{1+e^{-2k{x}_i+k}}+c\right)\end{array}\right.$(7)

where $Q_{\left(i,expect\right)},Q_{\left(i,current\right)}$ and $Q_{\left(i,initial\right)}$ are the expected, current, and initial amounts of urban-land seeds, respectively, in each ring. The step size of each iteration was set to N; N urban seeds were randomly discarded each time. x_i denotes the demand coefficient for each ring. y_i represents the i-th buffer penalty term. Notably, altering the c and k values enabled us to modify the morphology of the sigmoid-shaped function. D and D_T in the equation represent the parameters related to the action time of the penalty term. If the multiple relationship, D_T, between the current growth number of cells and the expected number reached the critical value D, the penalty term began working.

3.5. Accuracy assessment

An accuracy assessment was conducted using traditional classification consistency evaluation indices, including the Cohen’s Kappa coefficient (Kappa) and figure-of-merit (FoM) indices. To fill the research gap regarding the verification of model form reproduction ability, we constructed a morphological similarity (MS) index to evaluate morphological consistency.

Using a confusion matrix constructed from real-world images with predicted ones, Kappa was able to achieve the validation of the stimulated results for a global perspective. The equation for this calculation is shown in Equation (8)(8) $Kappa=\frac{P_{diagonal}-P_{product}}{1-P_{product}}$(8) (Liu et al. 2017):

(8) $Kappa=\frac{P_{diagonal}-P_{product}}{1-P_{product}}$(8)

where $P_{diagonal}$ refers to the fractional relationship between the diagonal figures of the confusion matrix constructed by all types and the total number of samples, i.e. the gross number of correctly classified samples divided by the total number of samples. $P_{product}$ denotes the actual quantity of each type multiplied by the predicted quantity of that type and divided by the square of the total sample number for the conditions of all types of sum figures. In this study, as the simulation target was urban land, we calculated the values for two types of objects: urban land and non-urban land.

The two indices, FoM and MS, were used to perform accuracy testing for a local perspective. In general, evaluation operations require the cooperation of newly-born points and their counterparts in the real situation. The equation used to calculate the FoM index is shown in Equation (9)(9) $FoM=\frac{N_{hit}}{N_{miss}+N_{hit}+N_{false}}$(9) (Pontius and Millones 2011; Pontius et al. 2008):

(9) $FoM=\frac{N_{hit}}{N_{miss}+N_{hit}+N_{false}}$(9)

$\mathrm{w}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{e}{N}_{miss},N_{hit},N_{false}$ are the number of cells whose actual state changes but the simulation result remains unchanged, the number of cells whose actual state changes and is consistent with the simulation result, and the number of cells whose actual state remains unchanged but the simulation result changes, respectively.

When we compared the images for the FoM and Kappa indices, it did not provide a morphological understanding; however, the MS index provided new insights into the matching of pixels. In general, it is difficult to visually check morphological differences in simulations. If we rely on the similarity within a certain range, the MS index can capture micro morphological features by window scanning. We followed the approach adopted by Yang et al. (2023) and used the following functions:

(10) $\{\begin{array}{c}{\mathrm{f}}_{\mathrm{i}}=a^{\left(-\frac{d_i}{\epsilon }\right)}\\ MS=\frac{SUM\left(f_m\right)}{N}\end{array}$(10)

The exponential variant $f_i$ was used to quantify the distance attenuation of each city cell from the center cell in the selected neighborhood window i×i. In Equation (10)(10) $\{\begin{array}{c}{\mathrm{f}}_{\mathrm{i}}=a^{\left(-\frac{d_i}{\epsilon }\right)}\\ MS=\frac{SUM\left(f_m\right)}{N}\end{array}$(10) , d_i is the distance between the cells in a predefined window and the center cell. $\alpha$ and $\epsilon$ are related to the speed of distance-based decay. $f_m$ is the average value of $f_i$ after each scan. The average sum value of $f_m$ was determined as the final MS value. Theoretically, the MS index ranged from 0 to 1. Notably, both 0 and 1 represented extreme cases. The former indicated that there was no urban land in the study area, and the latter denoted that the study area was completely covered by urban land. Therefore, the value of the MS index varied according to the actual situation. In this experiment, substantial attempts indicated that the MS index ranged from 0 to 0.05. Applying this index to the results of different models and then, conducting a comparison with the value in a real-world situation can help us better determine the simulation accuracy.

4. Results

4.1. Performance of geographic micro-process (GMP) fitting

We applied the Levenberg – Marquardt algorithm in the GMP model; the GMP-fitting results for the period of 1995–2020 are exhibited in (Table 1 and Figure 4(a)). The six-year average fitting R² value was 0.9507, which proved that the GMP was an excellent model for depicting the phenomenon of urban-land distance attenuation. The accumulation process of land cells was observed from the upward trend of the curve in (Figure 4(a)). This implies that urban landscapes are undergoing complex changes (due to the growing buoyancy). According to the $\alpha$ parameter (Table 1 and Figure 4(b)), we divided the urban land development into two phases, i.e. the increasing (1995–2005) and decreasing (2010–2020) phases. In the first stage, there was a slight increase of 10.35%, from 1.9512 to 2.1531, indicating that the city followed a compact development strategy during this period. However, after reaching the peak, this figure continued to fall during the later phase (from 2.1483 to 1.9273), with an amplification of 10.29%, implying that urban growth was more sprawled during this period. This discrepancy can be observed in (Figure 4(c)). Newly born cells tended to place their roots in the middle and rear buffer zones, demonstrating a backward morphological migration. For a similar $\alpha$ value, a higher n value represented a cumulative effect of quantity. As shown in Table 1, from 1995 to 2020, n increases drastically from 35.5309 to 151.2663, even though the value of $\alpha$ was similar, indicating the rapid growth of land coverage in the city.

Figure 4. Illustration of the geographic micro-process (GMP) model used in this study. (a) GMP-fitting results; (b) portrays the variations in the$\alpha$ parameter. The graphical shapes represent different time nodes. (c) increment proportion of land cells in each buffer ring. Green-based color denotes the data for 1995–2005, yellow-based color denotes the data for 2010–2020.

Table 1. Geographic micro-process (GMP)-fitting results and calculated cell number in different scenarios.

Year	Cells number	$\alpha$	n	R^2
1995	342969	1.9512	35.5309	0.9527
2000	457734	2.0099	61.6351	0.9320
2005	566709	2.1531	178.7755	0.9346
2010	774160	2.1483	180.0099	0.9529
2015	939769	2.0763	190.0693	0.9657
2020	1089207	1.9273	151.2663	0.9661

4.2. Analysis of modeling results

When the optimal combination of hyperparameters was found using the GA, without the addition of macro-level control, our baseline model (Figure 5(a)) was achieved. Then, by debugging the penalty term parameters, we dynamically allocated the GMP-fitting results to the microsimulation module in multiple-parameter combination modes (Figure 5(c)). In contrast to the information shown in (Figure 5(b)), without using GMPs as additional information, the stimulated map tended to gather in the front area, leaving the rear part with a lower number of cells, which demonstrates the intensity anomaly shown in (Figure 5(d)). By comparing the various combinations of parameters, we noted that the high consistency of the intensity value indicated that the performance of the penalty term was robust (Figure 5(d)). We further examined the fitting diagram (Figures 3 and 4(a)); keeping the k parameter constant, the curve became flatter with increasing parameter value, and it shifted in the opposite direction when we exchanged the setting. A steeper curve indicated that more variable space was provided for the input value and that the penalty term was more sensitive. Thus, we observed the different values of these parameters to explore the differences in the accuracy between the GMP-based and the baseline normal models in all the buffers. As shown in (Figure 6), all GMP-based models were the best performers, in terms of the Kappa value comparison for the tail-part buffers (Figure 6(b)). Similarly, (Figure 6(c)) portrays that in the GMP-based models, the FoM value was more effective for distance attenuation; however, the baseline model was the best performer for the front buffers. As shown in (Figure 6(d)), we noted complex back-and-forth fluctuations in each curve; however, the GMP-based model was more similar to the real-scenario curve.

Figure 5. Comparison of simulation results and demand intensity of each buffer for 2020. (a–c) simulation results of non- geographic micro-process (GMP)-based model, true situation, and GMP-based models respectively. (d) graph portraying the standardized distribution of the newly-born cells.

Figure 6. Fitting diagram of the penalty term and the accuracy evaluation results of each model. (a) penalty-fitting results; (b) Kappa value results; (c) figure-of-merit (FoM) index results; (d) morphological similarity (MS) index results; geographic micro-process (GMP).

5. Future scenario simulations

5.1. Urban intensity gradient forecasting

Based on the number of cells in the original multi-period urban data, the urban intensity gradients for 2025 and 2030 were predicted using the gray prediction algorithm (at strict 5-yr intervals). In general, this gray-based model extrapolates urban cell development trends by searching for heterogeneity among the input factors. Even if an unexpected situation arises, e.g. weak historical sequence traits or an insufficient sampling amount, the difficulties can be managed (Hsu and Chen 2003; Kayacan, Ulutas, and Kaynak 2010).

According to the principle of the GMP model, $\alpha$ is strongly linked to the land form; thus, it can serve as an excellent proxy for modeling in a variety of morphological situations. Owing to the fluctuation of history, we designed three different scenarios: (1) Stable scenario:$\alpha$ was set to the value in initial status; (2) Sprawl scenario: The maximum decrement in our record was given to the initial $\alpha$ value; (3) Compact scenario: The maximum increment was assigned to the initial $\alpha$ value. Through nonlinear programming analysis, we obtained an analytical solution for future urban intensity gradients.

5.2. Redefining accuracy in morphological perspective

The accuracy of scenario forecasting cannot be evaluated using Kappa, FoM, MS, or other related indicators. Therefore, we rebuilt a precision framework for assessing the scenario results, based on Angel et al. (2020). (1) Concentration: We utilized the urban compactness index (UCI), which is defined as the ratio of the radius of 75-% land density area (R^0.75) to the radius of 25-% land density area (R^0.25). (2) Connectivity: We calculated the ratio of the average distance from all the urban cells to the urban center and sub-centers (AD) to the average distance from all growable cells to the urban center and sub-centers (SD). (3) Exchange: We estimated the share of urban land occupied in a circle of equal area centered on the global centroid (SO).

5.3. Spatial simulation and analysis

The urban intensity gradients and GMP-related parameters for future land cells, results of GMP-fitting under various scenario in different future years, and accuracy verification results are shown in (Table 2, Figures 7(d), 8(d)), respectively. The statistical fluctuations (Table 2) and the overall rise in the curve shape from 2020 to 2030 (Figure 7(d)) both reveal that the n value can be a good proxy of the cumulative effect of urban cells that have equivalent $\alpha$ values. In addition, in our experiment, the adaptive gray model for future prediction was proven to be effective, based on the relative error between the two curves obtained by backtracking the fitting trend (Figure 8(d)). Furthermore, the land distribution in each buffer zone was difficult to discern with the naked eye. Therefore, we conducted observations on the utilization of the indicator system, to emphasize the assessment of morphological characteristics.

Figure 7. Simulations for various scenarios for 2025 and geographic micro-process (GMP)-fitting results. (a–c) results for compact, stable, and sprawl development, respectively; (d) GMP-fitting results for various scenarios for 2025 and 2030.

Figure 8. Simulation of various scenarios for 2030 and Grey model fitting results. (a–c) results for compact, stable, and sprawl development, respectively; (d) graph portraying the accuracy of the gray prediction algorithm.

Table 2. Geographic micro-process (GMP)-fitting results and calculated cell number for different scenarios.

Year	Cells number	Scenario	$\alpha$	n
2025	1377634	Compact	2.1531	441.0967
		Sprawl	1.9273	219.5153
		Stable	2.0444	315.2021
2030	1701490	Compact	2.1531	675.9306
		Sprawl	1.9273	331.3155
		Stable	2.0444	479.4774

The results in Table 3 demonstrate that the variation tendency of the three indicators conforms to the objective law of the predefined scenarios, i.e. the more compact the urban land, the greater the concentration and exchange value, and the smaller the connectivity degree. If more cells grow into the urban center and sub-centers, the size of the heartland area will increase; this synchronizes with the shrinking R^0.25 value and explains the decline in the AD value. Simultaneously, the radius of the equal-area circle (R) increased from 17.6690 (in 2025) to 19.7007 (in 2030), and the number of cells in the equal-area circle (N) of the looser growth mode was much smaller for the given R value, proving that our simulation results have high reliability.

Table 3. Evaluation results for different scenarios.

Year	Scenario	Type of indicator
		Concentration			Connectivity		Exchange
		R^0.75（km）	R^0.25（km）	UCI	AD（km）	AD/SD	R（km）	N	SO
2025	Compact	14.5438	30.1729	0.482	16.0011	0.7115	17.669	694854	0.6379
	Sprawl	13.8676	31.3177	0.4428	16.172	0.7191	17.669	677587	0.6221
	Stable	14.2321	30.6876	0.4638	16.3622	0.7276	17.669	686787	0.6305
2030	Compact	17.7281	36.7864	0.4819	12.7113	0.5652	19.7007	922399	0.6812
	Sprawl	17.1602	38.7704	0.4426	12.9104	0.5741	19.7007	904470	0.668
	Stable	17.4672	37.6742	0.4636	12.8081	0.5695	19.7007	914008	0.675

6. Discussion

6.1. Understanding the role of the penalty-term parameters

As shown in (Figure 6(a)), the sensitivity of the penalty term can be changed by adjusting the parameters of the functions c and k. The D parameter is strongly related to the duration of the penalty term; the larger the value of D, the later the action of the penalty term, which abates the correction force. The N parameter can control the step size of the iteration; the larger the N value, the faster the convergence speed of the model, and the more the probability of each buffer zone to deviate from a predetermined threshold. Thus, the sensitivity, duration, and step size play important roles in the composition, operation, and precision of the proposed penalty term. In Section 4.2, the aspects of land quantity variations and model precision fluctuations in each buffer have not been demonstrated effectively; however, these aspects are worth exploring.

In the absence of macro-constraints, each new urban cell tends to preferentially occupy a place, with a high global development probability; this leads to a distortion of the urban form (MS index column in Table 4). Conversely, after imposing the penalty term, the expected growth value and the original development probability compete with each other. First, the cells will seize the place with the greatest local development probability, with the dominance of a punishing force, i.e. central buffers with relatively high growth demand; this could explain the crests in all the curves (Figure 9). Then, the cells turn to former central buffers, which usually have a high global development probability, as a sign of waning growth demand in previous areas. When the growth requirements of these buffers were met, the advantage transitioned to the rear part of the buffers, reflected by the significant fluctuations in the corresponding positions (Figure 9).

Figure 9. Effects of each parameter of the penalty term: (a–d) effects of the c, k, D, and N parameters on the experimental results, respectively; geographic micro-process (GMP).

Table 4. Comparison of precision changes of different combinations of c, k, D, and N parameters.

Type of model					Kappa	FoM	MS
Non-GMP-based CA					0.8389	0.2862	0.01402
GMP-based CA (no punishment)					0.8459	0.2998	0.01382
GMP-based CA	c = 0.2	k = 6	D = 0.8	N = 300–1000	0.8516	0.3024	0.01340
	c = 0.8	k = 6	D = 0.8	N = 300–1000	0.8518	0.3024	0.01338
	c = 0.8	k = 6	D = 0	N = 300–1000	0.8516	0.3019	0.01337
	c = 0.8	k = 16	D = 0.8	N = 300–1000	0.8516	0.3018	0.01337
	c = 0.8	k = 6	D = 0	N = 8000–10000	0.8519	0.3034	0.01341

*MS of the true scenario is 0.01336.

Further empirical comparisons provide evidence that when the modifications i.e. larger c and k parameters, were set to reduce sensitivity (Figure 6(a)), the punishing effects weakened, and the deviation from the predefined expectations in the buffer zones was smaller (Figures 9(a, b)). In addition, when a longer-lasting penalty term appears, e.g. a smaller D parameter, the punishing effect is suppressed (Figure 9(c)). However, the step size has the most significant impact; the larger the N value, the greater the effect of the penalty term (Figure 9(d)). Table 4 portrays the strength of the punishment term, which was usually proportional to the Kappa or FoM values; however, an inverse relationship was observed in the morphological performance. Specifically, soft constraints rendered urban cells more likely to search for locations that had a relatively high development probability, especially as the number of buffers gradually increased and the simulation capability gradually declined. In such circumstances, hardline control of cellular growth at positions of low developmental probability is often counterproductive.

6.2. Implication of simulating urban land form

The urban form is very important to the development of a city, and controlling it is imperative when simulating the growth process of the city. In our study, we learned that CA-based urban simulation models largely relied on the influence of the development of the suitability of probability and the inter-neighborhood attraction at the micro-level; this kind of classical simulation paradigm has been ingrained in several existing studies. However, each city has complex systems. The agglomeration effect of its center and sub-centers leads to an allometric relationship between various urban elements and the urban land; notably, urban core areas have a much denser economic and social factor distribution. Unless the model is provided with extra information about distant areas, such as the planning of development zones, considering the total growth number of the study city without implementing more macro-level controls only allows urban cells to gather in places near the core area; thus, the suburban region lose the opportunities for growth. However, tough constraint methods often fall short when the probability is generally low in certain buffers. In this case, by adding detailed macro-requirements combined with the “soft control” strategy, cells can intelligently turn to other buffers that have a high growth probability, without significantly comprising the simulated form.

The good fitting performance of the GMP model demonstrates its excellent ability to characterize urban landforms. From the perspective of land microcosmic growth, the geographic entity can be abstracted into a series of interrelated and continuous units in time and space, to synchronously connect with the morphological evolution of the entire city. Therefore, such a function can help us understand the formation mechanism of urban form and have a better grasp of the future direction of landform prediction. Altering the α parameter of the GMP model allows urban land to decline more sharply or more gently at the macro level; this can help the model to control the compactness of the urban form. In addition, the fluctuation range of $\alpha$ does not deviate greatly under the observed trend in a long timeseries. Thus, it is feasible to integrate historical experiences into scenarios for different forms of future urban development.

Most previous studies aimed to test the accuracy of the results from the perspective of classification consistency. However, a higher Kappa or FoM value may not represent a better morphological consistency in the simulated results. Notably, when the interval between each simulation was short, the quantitative change in land growth did not reach the morphological qualitative change standard that could be observed directly. Such morphological errors continue to accumulate over time, and a series of simulation results under this paradigm will eventually leave the real development track of urban land. In this study, beyond the use of traditional indicators, such as Kappa and FoM, we constructed an index based on a microscopic MS test, which could help scholars discover the deviations in the simulated morphology over time. Note that the MS test is based on the distance attenuation feature to the specified center, to realize morphological difference discrimination, while focusing on the similarity of the neighborhood and without requiring the target map one-for-one to source records. However, the method offers all cells a relatively looser matching space. Thus, this type of fuzzy testing can ensure that land morphology maintains a high degree of similarity, without compromising accuracy. The main questions that arose from this analysis were assessing how to: 1) quantify the morphological characteristics of land growth under different scenarios and 2) design detailed index systems to verify the accuracy of multi-scenario results. The responses to these two questions have seldom been reported.

Scholars must realize the importance of morphological tests. The macro-level forewarned the model of the expected demand of all buffers; however, each dynamic evolution at the micro-level may not be able to transport cells to their real positions. Therefore, a comparison of the fluctuations in land increments for different scenarios cannot guarantee microscopic correctness.

The evaluation framework proposed in this study overcomes these limitations to some extent. The concentration index abstracts urban land, quantifies the macro increment distribution for multiple scenarios from the perspective of the land density ratio, and then, determines whether each result meets the required goal. The latter two indicators, namely, connectivity and exchange, give more consideration to the differences in the cellular distribution of each scenario in the microspace. Thus, the criteria that can be based on a combination of macro- and micro-levels can facilitate a better exploration of future urban growth scenarios.

6.3. Limitations and future directions

In this study, we demonstrated a top-down and bottom-up coupling model based on soft constraints to simulate the dynamics of urban land growth. This approach has several limitations.

First, on the macro scale, the urban intensity gradient derived from trend extrapolation may lack the interference from exogenous factors, and the effect of spacing between buffer rings needs detailed studies. Specifically, these effects can cause future forecasting values to fluctuate severely outside historical changes; therefore, a lack of such considerations may affect the precision of the model. Similarly, multi-scenario settings also rely on past experience to realize the allocation of macrogradients; more intricate factors need to be studied to address this issue. Second, a shorter simulation period generally led to higher accuracy. As the simulation continues, it is necessary to constantly update the impact factors related to development potential; otherwise, the future prediction accuracy of a longer timeseries will be limited. Third, the combination of penalty-term parameters has a subtle effect on the accuracy and form of the results; an optimal combination is required to maximize the model performance.

Our future study will be devoted to the construction of systematic macro-gradient acquisition rules and the development of a more robust micro-dynamic evolution mechanism and a comprehensive urban-form structure evaluation framework, while moving toward a higher dimension and carrying out simulations at larger scales.

7. Conclusion

Due to complicated urban forms, the reconstruction of the fundamental expression of cities in chaos can benefit the studies on urban cognition, management order construction, and sustainable urban development. Untangling the vast web of cities also emphasizes the search for order in complexity. To address the shortcomings of the top-down and bottom-up coupling models proposed in the existing studies, we present a soft constraint method and a more comprehensive urban morphological evaluation system to achieve an urban land dynamic simulation.

Notably, the spatiotemporal evolution of cities can be considered as the continuous geographical accumulation of a meta-process. Applying this understanding to describe and forecast incremental land changes offers valuable insights into the allocation rules of macrogradients. Furthermore, creating breakthrough soft constraint strategies to couple the macroscopic and microscopic modules is conducive to improving the accuracy of the model and progressing toward a more realistic land form. We are also committed to developing morphological evaluation indicators from macro- and micro-perspectives, to understand the accuracy of land morphology in the simulation process. In our study, the Kappa, FoM, and MS values reached 0.8506, 0.3034, and 0.01341, respectively, for the optimal parameter combination, indicating that the proposed model had a high simulation accuracy and a good imitative ability. Finally, we addressed the fact that the evolution of cities is closely related to their form. Thus, sorting the universal law in this staggered relationship requires us to jump out of the original understanding paradigm and embrace more challenging simulation strategies.

Acknowledgments

All the authors would like to express their gratitude for the support provided by the National Natural Science Foundation of China (Grant No. 41971368), as well as for the constructive comments and suggestions from the editors and three anonymous reviewers. Their meticulous reviews have greatly contributed to the enhancement of this manuscript.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability statement

The data that support the findings of this study are available from the first author upon reasonable request.

Additional information

Funding

This work was supported by the National Natural Science Foundation of China [41971368].