A cellular automata approach of urban sprawl simulation with Bayesian spatially-varying transformation rules

ABSTRACT

Incorporating spatial nonstationarity in urban models is essential to accurately capture its spatiotemporal dynamics. Spatially-varying coefficient methods, e.g. geographically weighted regression (GWR) and the Bayesian spatially-varying coefficient (BSVC) model, can reflect spatial nonstationarity. However, GWR possess weak ability eliminating the negative effects of non-constant variance because the method is sensitive to data outliers and bandwidth selection. We proposed a new cellular automata (CA) approach based on BSVC for multi-temporal urban sprawl simulation. With case studies in Hefei and Qingdao of China, we calibrated and validated two CA models, i.e. CA_BSVC and CA_GWR, to compare their performance in simulating urban sprawl from 2008 to 2018. Our results demonstrate that CA_BSVC outperformed CA_GWR in terms of FOM by ~2.1% in Hefei and ~3.6% in Qingdao during the calibration stage, and showed more accuracy improvement during the validation stage. The CA_BSVC model simulated urban sprawl more accurately than the CA_GWR model in regions having similar proximity to the existing built-up areas, especially in less developed regions. We applied CA_BSVC to predict urban sprawl at Hefei and Qingdao out to 2028, and the urban scenarios suggest that the proposed model shows better performance and reduced bias in reproducing urban sprawl patterns, and extends urban simulation methods by accounting for spatial nonstationarity.

KEYWORDS:

• Urban sprawl modeling
• cellular automata
• spatially-varying coefficient
• spatial nonstationarity
• figure-of-merit (FOM)

1. Introduction

Urbanization alters the land surface from its natural state to an artificial urban state, often leading to important impacts on the natural environment and human society (Karimi et al. 2019; Liu et al. 2019). Urban sprawl has not only led to substantial changes in land use and greater economic prosperity (Hegazy and Kaloop 2015), but also resulted in local climate change, greenhouse gas emissions, urban heat islands and increased flood risks (Baur, Förster, and Kleinschmit 2015; Chen et al. 2019; Estoque, Murayama, and Myint 2017; Li et al. 2017b; Muis et al. 2015). Urban sprawl varies spatially because of the differences in topography, population, economy and policy (Ramezani, Haddad, and Geroliminis 2015; Xu et al. 2018). Some models and analyses tools of urban sprawl have taken into account spatial nonstationarity, but most existing models are incapable of capturing the spatially nonstationary property of urban development (Mirbagheri and Alimohammadi 2017). To better forecast future urban scenarios, it is necessary to develop new urban simulation methods that accurately address spatial nonstationarity. The excepted new methods should improve our understanding of spatial nonstationarity and in turn, better support urban planning policy and regulation (Chen et al. 2019; Debonne et al. 2018; Xia et al. 2019).

Among urban simulation methods, geographical cellular automata (CA) are bottom-up, self-organizing approaches to the urban simulation that focuses on historical pattern reconstruction and future scenario forecast (Liu et al. 2018; Sakieh et al. 2015; Wu et al. 2019). A typical CA model is characterized by five entities: cells, driver-based transformation rules, neighborhood effects, constraints, and stochastic factors. The overall land transformation rules are usually computed considering the effects of both the global and local factors driving urban sprawl (Ghosh et al. 2017; Mustafa et al. 2018). The integrated transformation rules consider all five entities to determine whether the state of each cell will change at the next time step (Clarke 2018; Shafizadeh-Moghadam et al. 2017a). CA models typically use biophysical, demographic, socioeconomic and infrastructural factors to define driver-based transformation rules and reflect urban dynamics (Kantakumar, Kumar, and Schneider 2016). These rules are conventionally acquired using logistic regression (LR), a technique assumes that each factor has the same effect everywhere in space. Many methods used in CA modeling have the similar assumption, where CA parameters are spatially invariant. Consequently, they cannot reflect spatial nonstationarity.

Spatial nonstationarity in urban development is manifested as spatially different promoting and hindering effects of the same driving factor (Feng and Tong 2019). Many metrics have been used to evaluate this phenomenon and a few CA models have been developed to properly address the spatial nonstationarity when simulating urban development (Chen et al. 2014; Huang et al. 2019). For instance, partition-based CA approaches divide the cell space into sub-regions and then use LR to identify the spatial differences in urban sprawl (Abolhasani et al. 2016; Helbich and Griffith 2016). Ke, Qi, and Chen (2015) employed both K-means and KNN-cluster methods to partition a study area into subregions, using a unique transformation rule in each sub-region to reflect the spatial differences. In CA modeling of urban development, spatial nonstationarity can also be incorporated into the neighborhood effect that varies spatially. To address this, a variable-grid CA model was proposed to aggregate local and regional urban dynamics simultaneously for urban simulation (Vliet, White, and Dragicevic 2009). Since then, better performance has been achieved by incorporating the distance-decaying neighborhood into CA models (Li et al. 2017a). Although these methods could somewhat reflect spatial nonstationarity, they all utilize spatially invariant parameters in the transformation rules. Most recently, Feng and Tong (2019) reported spatially nonstationary maps of land-use gradient, which were incorporated into the neighborhood of CA models that have substantial modeling performance improvement. This indicates that it is more important to explicitly include spatially-varying parameters into the driver-based transformation rules instead of the neighborhood effect.

Spatially-varying coefficient models (SVCMs) are a new class of methods to analyze spatial phenomena with the modeling parameters changing across space (Gelfand et al. 2003; LeSage 2004), explicitly delineating spatial nonstationarity (Duncan et al. 2019; Helbich and Griffith 2016; Hu et al. 2016; Krehl 2018). GWR is a typical SVCMs that incorporates spatial nonstationarity into rule parameterization without predefined regional categorization (partition) (Nakaya et al. 2005). It uses distance-weighted subsets of samples to produce regression coefficients for every point; these coefficients are sensitive to sample locations (Brunsdon, Fotheringham, and Charlton 1996). Recent studies show that, as compared to conventional CA models, the integrated model of GWR and CA (named CA_GWR) can better simulate urban sprawl and accurately evaluate the consequent losses of ecosystem service (Mirbagheri and Alimohammadi 2017; Feng and Tong 2018; Chen et al. 2020). While this method can portray spatial nonstationarity, it still has some features that may lead to modeling failure. First, outliers may affect the regression coefficients via spatial enclave effects, because GWR depends on the attributes of the surrounding samples (LeSage 2004). These outliers are often linked to unpredictable urban sprawl (caused, for example, by changing policies) that does not follow the general rules of urban development. Second, nonstationarity variance exists widely across the observations and samples, which may fail to satisfy the assumption of normality in modeling errors (Shing 2008). Moreover, related methods such as geographically weighted logistic regression (GWLR), multiscale geographically weighted regression (MGWR) and geographically weighted temporally correlated logistic regression (GWTCLR) could also be used to derive land transformation rules; however, these are variates of standard GWR and share similar features. In contrast, a Bayesian spatially-varying coefficient (BSVC) model can examine the influence of an outlier and detect non-constant variance to alleviate deviation in coefficient estimates (LeSage 2004). Therefore, it should be useful to couple BSVC with CA for modeling urban sprawl simulation. Here, we compare GWR and BSVC to test which method can better reproduce historical urban patterns and forecast future urban scenarios.

In this study, we chose Chinese cities to test the performance of CA models because they have experienced unprecedented economic growth and rapid extent expansion since 1978 (Teets 2012). From 1978 to 2018, the Chinese urban population increased from 172 million to 831 million, and the urbanization rate increased from 18% to 60% (National Bureau of Statistics of China 2019). To model this rapid urban land-use change, we integrated BSVC and CA to create CA_BSVC, a novel urban sprawl simulation model that attempts to reduce the influence of non-constant variance and outliers in the model estimates. Using Hefei and Qingdao of China as the test cases, we aim to (1) validate the proposed CA_BSVC model considering its ability to describe spatial nonstationarity, (2) compare the simulation accuracies and spatial differences of CA_BSVC with CA_GWR, a GWR-based CA model; and (3) analyze the spatiotemporal urban sprawl in Hefei and Qingdao, and forecast their future scenarios. This study can help to understand and portray spatial nonstationarity in urban sprawl, and improve CA methods for better urban simulation and forecast.

2. The study area and datasets

2.1. The study area

We selected two cites (i.e. Hefei and Qingdao; see Figure 1) to examine the proposed model. Hefei is the capital of Anhui Province in the western Yangtze River Delta (Figure 1(b)). It covers 11,445 km² and incorporates nine jurisdictions: four urban districts that define the Hefei City center (marked by red boundary lines) and five rural counties. The topography of Hefei is dominated by low-lying plains and hills, with more than 8 million people living mainly in the plain area. The city is located in the mid-latitude zone with distinct monsoon seasons. The average annual temperature is 15.7°C and the precipitation about 1000 mm. Another city, Qingdao, is a bustling and important port city in north China (Figure 1(c)), with seven urban districts and three counties under its jurisdiction, and the population was more than 9 million as of 2018 (Qingdao Municipal Bureau of Statistics 2018). Qingdao is a coastal hilly city with high terrain in the east and low terrain in the west. Among them, mountains account for 15.5% of the total area, hills 2.1%, plains 37.7%, and depressions 21.7%. The city has a temperate monsoon climate, with an average temperature of 25.3°C and an average annual precipitation of 662 mm. Both cities experienced great urban sprawl in the last decade, which mainly concentrated in their urban districts (Hao, Zhu, and Zhong 2015; Zhao and Zou 2018). We therefore selected the urban districts to test the performance of our model (CA_BSVC) by comparing it with a GWR-based CA method (CA_GWR).

Figure 1. Location of the study areas: Hefei (b) and Qingdao (c)

2.2. Input data and variables

The maps of land-use change and its driving factors are core input to the CA_BSVC model. We used cloud-free Landsat TM and Landsat OLI images to produce 2008, 2013 and 2018 land-use maps (Table 1). After geo-correction and radiometric-rectification, we applied a supervised Maximum-Likelihood classifier to these images and interpret land-use as urban, non-urban or water land categories. Six widely used factors were selected as the land-use change drivers in this study considering their ability to reflect the urban sprawl potential (Chen, Liu, and Li 2017; Kamusoko and Gamba 2015; Liao et al. 2014; Yao et al. 2017b). We also classified land-use change drivers as biophysical, socioeconomic and infrastructural (Table 1). These factors were used to construct CA transformation rules and produce transformation potential reflecting the probability of a cell converting its state from non-urban to urban. The biophysical factor was reflected by the digital elevation model (DEM) acquired from the Geospatial Data Cloud (www.giscloud.com). The socioeconomic factor denotes the number of persons per pixel (PPP) retrieved from the WorldPop project (www.worldpop.org.uk). The infrastructural factors include the proximity to subways, railways, highways and primary roads, all produced using Euclidean distances based on datasets derived from OpenStreetMap (www.openstreetmap.org). All factor maps were projected to UTM-WGS 1984 Zone 50 N and resampled at the same 30 m resolution (Figure 2).

Figure 2. Factors that drive urban sprawl at Hefei and Qingdao from 2008 to 2013

Table 1. Urban land-use maps and selected land-use change driving factors

Category	Name	Data source	Definition	Type
Land-use map	Land-2008	Landsat TM images from Geospatial Data Cloud	Land-use map for 2008	Dependent
Land-use map	Land-2013	Landsat TM images from Geospatial Data Cloud	Land-use map for 2013	Dependent
Land-use map	Land-2018	Landsat OLI images from Geospatial Data Cloud	Land-use map for 2018	Dependent
Driving factor	DEM	DEM product collected at Geospatial Data Cloud	Terrain condition reflected by elevation	Biophysical
Driving factor	PPP	WorldPop project website	Population density per pixel	Socioeconomic
Driving factor	DisSubway	Free editable world OpenStreetMap	Proximity to nearest subway	Infrastructural
Driving factor	DisRail	Free editable world OpenStreetMap	Proximity to nearest railway	Infrastructural
Driving factor	DisHighway	Free editable world OpenStreetMap	Proximity to nearest highway	Infrastructural
Driving factor	DisRoad	Free editable world OpenStreetMap	Proximity to nearest primary road	Infrastructural

3. Methods

3.1. An overview of our method

Our CA_BSVC model has three steps: data preprocessing, model construction, and model execution (Figure 3). We first classified the Landsat imagery to produce urban patterns at the start time and end time, and extracted factors driving urban sprawl from the raw data. We next used systematic sampling to select training data to build the transformation rules using GWR and BSVC. We then calibrated and validated the CA_GWR and CA_BSVC models to simulate the historical urban patterns in the end year and compared these two models in detail. Finally, we applied the CA_BSVC model to forecast a future urban scenario for Hefei and Qingdao.

Figure 3. CA_BSVC workflow for modeling historical and future urban sprawl

3.2. The fundamental CA model

CA are computationally intensive models with self-organizing and bottom-up characteristics. They can be applied to reproducing complex land-use patterns and urban sprawl with a set of transformation rules and a series of model-defined interactions among neighboring cells (Ghosh et al. 2017; Grinblat, Gilichinsky, and Benenson 2016; Kim et al. 2017). One core feature of CA models is the use of specific rules to determine whether a cell will change from one state to another at the next time step. The cell states ($S_i^{t+1}$) at time step t + 1 are defined as (Barredo et al. 2003; de Almeida et al. 2003; Gao et al. 2020):

(1) $S_i^{t+1}=Fun\left(S_i^t,P_{{\mathrm{\Omega }}_i^t},P_{driver},Rest\left(S_i^t\right),Rnd()\right)$(1)

where Fun denotes the rules defining the total transformation probability $P_{tt}$; $S_i^t$ denotes the cell state at the present time step (t); $P_{{\mathrm{\Omega }}_i^t}$ denotes the neighborhood-based transformation probability; $P_{driver}$ denotes the driver-based transformation probability; $Rest\left(S_i^t\right)$ denotes the regions unavailable for urban development, including broad water bodies, natural reserves, and basic farmlands; and Rnd() denotes the unknown errors that can be given by (García et al. 2011):

(2) $Rnd()=1+{\left(-lnr\right)}^{\alpha }$(2)

where r is a random number in the range of [0,1]; and $\alpha$ is an exponential factor to scale the effects of randomness.

The neighborhood-based probability reflected by a neighborhood’s configuration and size is a key component of CA modeling (Dahal and Chow 2015; Wu et al. 2012). It is often defined using a regular symmetric window:

(3) $P_{{\mathrm{\Omega }}_i^t}\left(\mathrm{n}\mathrm{o}\mathrm{n}\mathrm{u}\mathrm{r}\mathrm{b}\mathrm{a}\mathrm{n}\mathrm{t}\mathrm{o}\mathrm{u}\mathrm{r}\mathrm{b}\mathrm{a}\mathrm{n}\right)=\frac{{\sum }_{i=1}^{m\times m-1}con\left(S_i^t=dev\right)}{m\times m-1}$(3)

where m × m represents the neighborhood size, and $con\left(S_i^t=dev\right)$ denotes the number of developed cells in the neighborhood. In this study, we applied a 5 × 5 square window to reflect neighborhood impacts.

The driver-based transformation probability is calculated using a set of drivers of urban sprawl defined through LR, a method that has been widely used in the literature (Wu 2002). To improve model robustness, Feng and Tong (2020) proposed a local adjustment parameter (LAP) to offset the increase in the neighborhood probability and a time-increment parameter (TIP) to offset the decay of the driver-based probability. As a result, the driver-based transformation probability ($P_{driver}$) and the total transformation probability ($P_{tt}$) are given by:

(4) $\left\{\begin{array}{c}{P}_{driver}\left(nonurbantourban\right)=\frac{exp\left(-\left({\sum }_{j=0}^k{a}_j{x}_j\right)+\epsilon \right)}{1+exp\left(-\left({\sum }_{j=0}^k{a}_j{x}_j\right)+\epsilon \right)}\\ \begin{array}{c}{P}_{tt}\left(nonurbantourban\right)=\frac{\left(P_{{\mathrm{\Omega }}_i^t}\text{cdotLAP}+P_{driver}\cdot {\left(1+TIP\right)}^{t-1}\right)\text{cdotRest}\left(S_i^t\right)}{2}\\ \end{array}\end{array}\right.$(4)

where LAP ranges from 0.5 to 1 and TIP ranges from 0.0 to 0.1; k is the number of driving factors; $a_j$ denotes the weight of factor x_j; and $\epsilon$ are the fitting residuals.

3.3. The CA_BSVC model

3.3.1. Spatial varying coefficient models

Compared with LR-based CA models, the most distinct differences of the CA_GWR and CA_BSVC models are their land transformation rules. These rules were trained using sample points (a total of 4,805) selected from the maps of urban sprawl (dependent variable) and driving factors (independent variables). The sample selection method in this study was the systematic sampling. We used SVCMs (e.g. GWR and BSVC) to define the weight of each driver, that is, we applied these two methods to train the land transformation rules. SVCMs examine spatial nonstationarity in natural phenomena based on the assumption that the regression coefficients are functions of their spatial locations (Xu and Huang 2015). GWR is a representative SVCM built by extending classical linear regression (Brunsdon, Fotheringham, and Charlton 1996). The GWR is defined as (LeSage 2004):

(5) $W_{\left(lat,lon\right)}Y=W_{\left(lat,lon\right)}X{a}_{\left(lat,lon\right)}+{\epsilon }_{\left(lat,lon\right)}$(5)

where $W_{\left(lat,lon\right)}$ is an n × n spatial weight matrix at location $\left(lat,lon\right)$; $Y$ is the dependent variable vector (n × 1); $X$ is the independent variable matrix (n × k); $a_{\left(lat,lon\right)}$ denotes the model coefficient vector at location $\left(lat,lon\right)$; and ${\epsilon }_{\left(lat,lon\right)}$ represents the random error vector that follows the $N\left(0,{\sigma }^2\right)$ distribution.

The definition of spatial weights $W_{\left(lat,lon\right)}$ are critical for the GWR coefficient estimation because they substantially affect the regression results. Commonly used spatial weight definitions are distance threshold, inverse distance and Gaussian function. For artificial land surface and suburban areas with highly clustered patterns, it is challenging to capture spatial variations using a small bandwidth, given the few samples included in local regression (Mirbagheri and Alimohammadi 2017). However, the Gaussian method (Mcmillen 2004) given by

(6) $w_{\left(lat\mathrm{\_}i,lon\mathrm{\_}i\right),\left(lat\mathrm{\_}j,lon\mathrm{\_}j\right)}=exp\left(\frac{-{d_{\left(lat\mathrm{\_}i,lon\mathrm{\_}i\right),\left(lat\mathrm{\_}j,lon\mathrm{\_}j\right)}}^2}{b^2}\right)$(6)

applies an exponential distance-decay function to generate spatial weights that sufficiently capture the spatial variation. Here, $w_{\left(lat\mathrm{\_}i,lon\mathrm{\_}i\right),\left(lat\mathrm{\_}j,lon\mathrm{\_}j\right)}$ represents the weight given to the processing location $\left(lat\mathrm{\_}j,lon\mathrm{\_}j\right)$ on the regression location $\left(lat\mathrm{\_}i,lon\mathrm{\_}i\right)$; $d_{\left(lat\mathrm{\_}i,lon\mathrm{\_}i\right),\left(lat\mathrm{\_}j,lon\mathrm{\_}j\right)}$ is the distance between a given location $\left(lat\mathrm{\_}i,lon\mathrm{\_}i\right)$ and a sample location $\left(lat\mathrm{\_}j,lon\mathrm{\_}j\right)$; and b is the bandwidth. The CA_GWR model defined the bandwidth using a cross-validation method where a best bandwidth of 0.3163 was related to 33 neighboring points. The same bandwidth was applied to CA_BSVC, making sure the two models were compared under the same search range. The GWR method for training land transformation rules can be performed using the “GWR tool” in the UrbanCA software (Feng and Tong 2020).

3.3.2. A Bayesian spatially-varying coefficient model

As opposing to GWR, the BSVC modeling errors (${\epsilon }_i$) are assumed to vary across space. The errors were formulated as ${\epsilon }_{\left(lat,lon\right)}\text{~}N\left(0,{\sigma }^2{V}_{\left(lat,lon\right)}\right)$, where $V_{\left(lat,lon\right)}$ is a diagonal matrix ($n\times n$) representing the non-constant features (Subedi, Zhang, and Zhen 2018). BSVC applies the contiguity-based smoothing to reflect the relationships among the retrieved spatially-varying coefficients. The BSVC model coefficients are (LeSage 2004):

(7) $\begin{array}{rl}& \left\{\hspace{0.17em}\right.a_{\left(lat\mathrm{\_}i,lon\mathrm{\_}i\right)}=\left(W_{\left(lat\mathrm{\_}i,lon\mathrm{\_}i\right),\left(lat\mathrm{\_}1,lon\mathrm{\_}1\right)}{I}_k,\cdot \cdot \cdot ,\text{break}{W}_{\left(lat\mathrm{\_}i,lon\mathrm{\_}i\right),\left(lat\mathrm{\_}n,lon\mathrm{\_}n\right)}{I}_k\right)\left(\begin{array}{c}{a}_{\left(lat\mathrm{\_}1,lon\mathrm{\_}1\right)}\\ \cdot \cdot \cdot \\ a_{\left(lat\mathrm{\_}n,lon\mathrm{\_}n\right)}\end{array}\right)+u_{\left(lat\mathrm{\_}i,lon\mathrm{\_}i\right)}{u}_{\left(lat\mathrm{\_}i,lon\mathrm{\_}i\right)}\text{break}\text{~}N\left(0,{\sigma }^2{\delta }^2{\left(X^T{W_{\left(lat\mathrm{\_}i,lon\mathrm{\_}i\right)}}^{2}X\right)}^{-1}\right)\end{array}$(7)

where $W_{\left(lat\mathrm{\_}i,lon\mathrm{\_}i\right),\left(lat\mathrm{\_}j,lon\mathrm{\_}j\right)}$ is a row-normalized matrix; $I_k$ is the $k$ identity matrix; $u_{\left(lat\mathrm{\_}i,lon\mathrm{\_}j\right)}$ is a stochastic term; and ${\delta }^2$ is a scale factor controlling the smoothing influence. When ${\delta }^2$ is small, the smoothing method has more considerable impacts on the regression coefficients; when ${\delta }^2$ increases, the impacts gradually reduce to those given by GWR.

The calculation of regression coefficients using Equation (6)(6) $w_{\left(lat\mathrm{\_}i,lon\mathrm{\_}i\right),\left(lat\mathrm{\_}j,lon\mathrm{\_}j\right)}=exp\left(\frac{-{d_{\left(lat\mathrm{\_}i,lon\mathrm{\_}i\right),\left(lat\mathrm{\_}j,lon\mathrm{\_}j\right)}}^2}{b^2}\right)$(6) is much more complicated. Here, a Bayesian method is introduced to assist in solving the equation (LeSage 2004). Similar to other Bayesian methods, the core role of BSVC is calculating the posterior probability density for all model parameters. Regarding this, Markov Chain Monte Carlo (MCMC) and Gibbs sampling are typical Bayesian methods that can simulate probability density $P(\theta |Y)$ with large random samples (LeSage 2004). When using MCMC for sampling, the marginal distribution cannot be considered a stable state before the converge, and the process of removing the previous unstable value is called burn-in in the Bayesian methods.

To implement MCMC, we retrieved the parameters $a_{\left(lat,lon\right)}$,$\sigma$,$\delta$, V and y from the conditional posterior distribution $P(a_{\left(lat,lon\right)}|\sigma ,\delta ,V,y)$ where $y$ denotes the model coefficients for other observations. Figure 4 presents the Gibbs sampling procedure to define the driver-based transformation rules. Based on BSVC, the transformation probability in Equation (2)(2) $Rnd()=1+{\left(-lnr\right)}^{\alpha }$(2) can be rewritten as:

(8) $P_{driver}\left(\mathrm{n}\mathrm{o}\mathrm{n}\mathrm{u}\mathrm{r}\mathrm{b}\mathrm{a}\mathrm{n}\mathrm{t}\mathrm{o}\mathrm{u}\mathrm{r}\mathrm{b}\mathrm{a}\mathrm{n}\right)=\frac{\mathrm{e}\mathrm{x}\mathrm{p}\left(W_{\left(lat,lon\right)}X{a}_{\left(lat,lon\right)}+{\epsilon }_{\left(lat,lon\right)}\right)}{1+exp\left(W_{\left(lat,lon\right)}X{a}_{\left(lat,lon\right)}+{\epsilon }_{\left(lat,lon\right)}\right)}$(8)

Figure 4. Workflow for the definition of driver-based transformation rules using BSVC

3.3.3 Model assessment

We used figure-of-merit (FOM) and total operating characteristic (TOC) to evaluate model performance. FOM is based on the overlay of the actual starting map (IniUrban), the actual ending map and the simulated ending map. Through map overlaying, we identified output cells as Initial Urban (IniUrban), Correct Simulation (CorrectSim), Missed Simulation (MissSim), False Simulation (FalseSim), Nonurban Persistence (PerCell) and Water Bodies (Pontius et al. 2008b). TOC is a ratio-free metric that can demonstrate model performance in land-use change modeling and judge the spatial distribution relationship between the existing urban area and growth urban area (Pontius Jr and Si 2014). The FOM and TOC metrics can be given by (Tong and Feng 2020):

(9) $\begin{array}{c}\\ FOM=\frac{\mathrm{C}\mathrm{o}\mathrm{r}\mathrm{r}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{S}\mathrm{i}\mathrm{m}}{\mathrm{C}\mathrm{o}\mathrm{r}\mathrm{r}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{S}\mathrm{i}\mathrm{m}+\mathrm{M}\mathrm{i}\mathrm{s}\mathrm{s}\mathrm{S}\mathrm{i}\mathrm{m}+\mathrm{F}\mathrm{a}\mathrm{l}\mathrm{s}\mathrm{e}\mathrm{S}\mathrm{i}\mathrm{m}}\end{array}$(9)

(10) $\left\{\begin{array}{c}{\mathrm{X}}_{\mathrm{t}\mathrm{o}\mathrm{c}}=\frac{\mathrm{F}\mathrm{a}\mathrm{l}\mathrm{s}\mathrm{e}\mathrm{S}\mathrm{i}\mathrm{m}}{\mathrm{F}\mathrm{a}\mathrm{l}\mathrm{s}\mathrm{e}\mathrm{S}\mathrm{i}\mathrm{m}+\mathrm{P}\mathrm{e}\mathrm{r}\mathrm{C}\mathrm{e}\mathrm{l}\mathrm{l}}\\ {\mathrm{Y}}_{\mathrm{t}\mathrm{o}\mathrm{c}}=\frac{\mathrm{C}\mathrm{o}\mathrm{r}\mathrm{r}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{S}\mathrm{i}\mathrm{m}}{\mathrm{C}\mathrm{o}\mathrm{r}\mathrm{r}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{S}\mathrm{i}\mathrm{m}+\mathrm{M}\mathrm{i}\mathrm{s}\mathrm{s}\mathrm{C}\mathrm{e}\mathrm{l}\mathrm{l}}\end{array}\right.$(10)

where ${\mathrm{X}}_{\mathrm{t}\mathrm{o}\mathrm{c}}$ and ${\mathrm{Y}}_{\mathrm{t}\mathrm{o}\mathrm{c}}$ jointly define each point in TOC curves.

4. Results

4.1. Land transformation probability maps

We used the methods in Equation (7)(7) $\begin{array}{rl}& \left\{\hspace{0.17em}\right.a_{\left(lat\mathrm{\_}i,lon\mathrm{\_}i\right)}=\left(W_{\left(lat\mathrm{\_}i,lon\mathrm{\_}i\right),\left(lat\mathrm{\_}1,lon\mathrm{\_}1\right)}{I}_k,\cdot \cdot \cdot ,\text{break}{W}_{\left(lat\mathrm{\_}i,lon\mathrm{\_}i\right),\left(lat\mathrm{\_}n,lon\mathrm{\_}n\right)}{I}_k\right)\left(\begin{array}{c}{a}_{\left(lat\mathrm{\_}1,lon\mathrm{\_}1\right)}\\ \cdot \cdot \cdot \\ a_{\left(lat\mathrm{\_}n,lon\mathrm{\_}n\right)}\end{array}\right)+u_{\left(lat\mathrm{\_}i,lon\mathrm{\_}i\right)}{u}_{\left(lat\mathrm{\_}i,lon\mathrm{\_}i\right)}\text{break}\text{~}N\left(0,{\sigma }^2{\delta }^2{\left(X^T{W_{\left(lat\mathrm{\_}i,lon\mathrm{\_}i\right)}}^{2}X\right)}^{-1}\right)\end{array}$(7) to generate GWR and BSVC transformation probability maps in Hefei and Qingdao by considering the land-use and driving factor maps. Figure 5 displays the initial urban area and urban sprawl during 2008–2013 in Hefei (Figure 5(a)) and Qingdao (Figure 5(d)), respectively. The urban sprawl in Hefei occurred mainly at the periphery of the existing urban area. Limited by the spatial extent of the existing urban area and the hills to the east, Qingdao’s major urban sprawl increased northward along the edge of the hills. We also noticed evident urban sprawl in Huangdao district (the region 11 in Figure 1(c)). The transformation probability maps show that large values are related to clustered urban sprawl occurring around urban districts. The regions with high and middle transformation probability (in red and yellow) in BSVC are substantially fewer than GWR, where BSVC generated transformation probability more consistent with the actual urban sprawl between 2008 and 2013. Furthermore, the BSVC transformation probability map in Qingdao accurately reflected the urban sprawl potential in the northern part. While in GWR, transformation probabilities of this part are quite similar to those in the surrounding areas.

Figure 5. 2008–2013 urban sprawl and the GWR and BSVC transformation probability maps: (a) Hefei’s urban sprawl; (b) GWR transformation probability map in Hefei; (c) BSVC transformation probability map in Hefei; (d) Qingdao’s urban sprawl; (e) GWR transformation probability map in Qingdao; and (f) BSVC transformation probability map in Qingdao

4.2. Model calibration

The calibration of the CA_GWR and CA_BSVC model was done with the input data and the retrieved CA transformation rules that capture urban sprawl at Hefei and Qingdao from the 2008–2013 interval. Visual inspection reveals that both the simulated urban distributions cluster in the city center, with more isolated urban sprawl patches in actual pattern (Figures 5 and 6). Only minor differences are found between the CA_GWR and CA_BSVC simulations, where urban sprawl modeled by GWR appears in the urban fringes. In contrast, the urban sprawl modeled by BSVC gives high transformation probability to small cities, and visual inspection shows that urban sprawl in small cities by BSVC is greater than by GWR.

Figure 6. The 2013 simulated urban patterns based on the CA_GWR and CA_BSVC models: (a) Hefei’s urban sprawl simulated by CA_GWR; (b) Hefei’s urban sprawl simulated by CA_BSVC; (c) Qingdao’s urban sprawl simulated by CA_GWR; and (d) Qingdao’s urban sprawl simulated by CA_BSVC.

Table 2 shows that CA_BSVC with FOM 37.2% for Hefei and 56.8% for Qingdao performed better than GWR with FOM 35.1% for Hefei and 53.2% for Qingdao. These indicate that the FOMs of CA_BSVC in both cities are about 3% higher than CA_GWR, indicating better modeling accuracy of CA_BSVC. Compared to CA_GWR, CA_BSVC correctly simulated more urban sprawl (0.5% and 0.4%, respectively), missed less urban sprawl (0.5% and 0.4%, respectively), and made less false urban sprawl (0.1% and 1.3%, respectively). These metrics all show that CA_BSVC is better than CA_GWR in simulating the 2013 urban patterns for both Hefei and Qingdao.

Table 2. Model assessment for CA_GWR and CA_BSVC in 2013

City	Model	Assessment component
		IniUrban	CorrectSim	MissSim	FalseSim	PerCell	Water	FOM
Hefei	GWR (%)	29.8	9.1	3.8	13.0	34.6	9.6	35.1
	BSVC (%)	29.8	9.6	3.3	12.9	34.7	9.6	37.2
	BSVC minus GWR (%)	0.0	0.5	−0.5	−0.1	0.1	0.0	2.1
Qingdao	GWR (%)	9.9	17.0	3.9	11.1	52.8	5.3	53.2
	BSVC (%)	9.9	17.4	3.5	9.8	54.1	5.3	56.8
	BSVC minus GWR (%)	0.0	0.4	−0.4	−1.3	1.3	0.0	3.6

While metrics show the statistical aspect of the model performance, the spatial patterns can display more about the performance. We therefore overlaid the 2008 actual pattern, the 2013 actual pattern and the 2013 simulated pattern to produce new maps for evaluating the spatial successes and errors (Figure 7). Although the two pair of evaluation maps are similar in most areas, there are distinct differences in specific regions. For example, in the black ellipses, CA_BSVC for Hefei has more clustered CorrectSim and fewer clustered FalseSim than by CA_GWR (Figure 7(a,b)). Similar results of fewer clustered FalseSim have been observed in Qingdao, especially in the black ellipses (Figure 7(c,d)). For both models, CorrectSim and MissSim are dispersedly distributed while more clustered FalseSim can be detected in the CA_GWR simulations. These indicate that clustered PerCell are more likely misclassified as urban by the CA_GWR model. Overall, CA_BSVC outperformed CA_GWR on both successes and errors in 2013, especially for CorrectSim (more) and MissSim (fewer).

Figure 7. The 2013 simulation successes and errors for CA_GWR and CA_BSVC with selected regions: (a) Hefei’s urban sprawl simulated by CA_GWR; (b) Hefei’s urban sprawl simulated by CA_BSVC; (c) Qingdao’s urban sprawl simulated by CA_GWR; and (d) Qingdao’s urban sprawl simulated by CA_BSVC.

4.3. Model validation and assessment

To validate the CA_BSVC model, we predicted the 2018 urban pattern in Hefei and Qingdao by using the 2013 actual map as the initial map (Figure 8). Minor spatial differences can be detected by the comparison between BSVC and GWR in Hefei, which could be attributed to the fact that most regions have been urbanized in 2013 and there are fewer regions left in 2018 suitable for further urban sprawl (Figure 5(a)). However, the simulation results of the two models varied considerably in the northern part of Qingdao, especially in the regions A and B (Figure 8(c,d)). Table 3 reports sharp increases in FOM and CorrectSim by the comparison between the BSVC and GWR simulation results in Hefei, indicating that BSVC has higher accuracy when the urban sprawl area is relatively small. A higher FOM change in 2018 than that in 2013 also took place in Qingdao, where the decrease of False alarms played a key role.

Figure 8. The 2018 simulated urban patterns based on the CA_GWR and CA_BSVC models: (a) Hefei’s urban sprawl simulated by CA_GWR; (b) Hefei’s urban sprawl simulated by CA_BSVC; (c) Qingdao’s urban sprawl simulated by CA_GWR; and (d) Qingdao’s urban sprawl simulated by CA_BSVC.

Table 3. Model assessment for CA_GWR and CA_BSVC in 2018

City	Model	Assessment component (%)
		IniUrban	CorrectSim	MissSim	FalseSim	PerCell	Water	FOM
Hefei	GWR	42.8	0.7	3.9	5.3	37.7	9.6	6.7
	BSVC	42.8	1.6	3.0	8.2	34.8	9.6	12.4
	BSVC minus GWR	0.0	0.9	−0.9	2.9	2.9	0.0	5.7
Qingdao	GWR	30.8	14.5	2.4	15.3	31.8	5.3	48.7
	BSVC	30.8	12.9	4.0	9.7	37.2	5.3	52.9
	BSVC minus GWR	0.0	−1.6	1.6	−5.6	5.6	0.0	4.2

The stacked map of the 2013 actual map, 2018 actual map and 2018 simulated map yields the distribution of successes and errors of CA_BSVA and CA_GWR in the validation stage. Due to less urban sprawl in Hefei, the overall distribution of the stacked map is relatively consistent. Visually, in the west and north of Hefei, CA_BSVA presents higher CorrectSim than CA_GWR; in Qingdao, the clustered FalseSim is observed in the CA_GWR simulation (Figure 9(d)), which is similar to the CA_GWR simulation results from 2008 to 2013 and seriously affected the model accuracy.

Figure 9. The 2018 simulation successes and errors for CA_GWR and CA_BSVC with selected regions: (a) Hefei’ urban sprawl simulated by CA_GWR; (b) Hefei’ urban sprawl simulated by CA_BSVC; (c) Qingdao’ urban sprawl simulated by CA_GWR; and (d) Qingdao’ urban sprawl simulated by CA_BSVC.

Figure 10 illustrates the TOCs of BSVC and GWR in Hefei and Qingdao. The areas under curve (AUC) of BSVC and GWR for Hefei are about 0.7, and those for Qingdao are about 0.9, which are all greater than that of the uniform baseline (0.5). These metrics reveal that urban sprawl in Hefei and Qingdao tended to occur in the places near the existing urban areas. The curves of BSVC and GWR for Qingdao are consistently higher than the Uniform, while the curves for Hefei are lower than the Uniform line within 200 km² and then exceed the line. Based on the FOM, the success and error maps, and the TOC, we concluded that BSVC is superior to GWR in urban sprawl simulation. We therefore applied CA_BSVC to predict future urban sprawl in both Hefei and Qingdao.

Figure 10. Using TOC to evaluate the simulation performance of CA_GWR and CA_BSVC: (a) Hefei; and (b) Qingdao

4.4. Future predictions and urban dynamics

We forecasted the urban scenarios for 2028 using CA_BSVC to better understand urban land-use change (Figure 11), and then delineated the hot-zones of urban development with 5 km-wide concentric buffers centered to the urban centroid. Our remote sensing classification shows that Hefei’s urban areas increased from 192 km² in 2008 to 305 km² in 2018 and are expected to reach 438 km² by 2028 based on the CA-Markov estimation, with an annual urban sprawl of 13 km² during 2018–2028. Qingdao’s urban areas increased from 135 km² in 2008 to 647 km² in 2018 and are expected to reach 933 km² by 2028, corresponding to a total urban sprawl rate of 44% during 2018–2028 with annual urban sprawl of 29 km². Most urban sprawl in Hefei is encompassed by the three innermost concentric rings (buffer ring with 15 km wide), and urban sprawl in Qingdao is encompassed by the four innermost concentric ring Qingdao (buffer ring with 20 km wide). Table 4 presents that the largest urban patches in Hefei will occur in the second and third rings, which are 5 ~ 15 km away from the urban centroid. Qingdao’s urban sprawl is farther from the urban centroid, mainly concentrates between the third to fifth rings, with the most occurring in the fourth ring (15–20 km).

Figure 11. The 2018–2028 urban sprawl and the concentric buffers for Hefei (a) and Qingdao (b)

Table 4. The urban sprawl in each concentric buffer during 2018–2028

Distance (km)	Hefei			Qingdao
	Area (km^2)	Percent (%)	Cumulative (%)	Area (km^2)	Percent (%)	Cumulative (%)
0–5	9.4	7.1	7.1	1.6	0.6	0.6
5–10	57.8	43.3	50.4	13.6	4.7	5.3
10–15	63.1	47.2	97.6	61.0	21.3	26.6
15–20	3.2	2.4	100.0	101.3	35.4	62.0
20–25				68.0	23.8	85.8
25–30				24.4	8.5	94.3
30–35				16.3	5.7	100.0

To further exhibit the spatial differences in urban sprawl in 2008–2028, we partitioned the study areas in four directions based on the urban centroids (Figure 12). The urban sprawl areas in Hefei in the four directions expected to exceed those in 2008, with 131.4 km², 106.4 km², 102.3 km² and 74.6 km² emerged in the east, west, north and south directions, respectively. Different from Hefei, a single central city located in flat terrain, Qingdao’s urban development is obviously blocked by mountains, hills and oceans. Urban areas in the south and north were dominant in 2008, accounting for 71.6% of the total. Benefit from a shorter distance to the existing urban areas, a quick increase in urban sprawl has been depicted in the east during 2008–2028 with 263.4 km² (Figure 12(c)). A similar sprawl occurs in the north, where urban area increases from 46.3 km² in 2008 to 285.0 km². Urban areas in the south and west have also grown considerably, reaching 151.0 km² and 144.8 km² in the two decades, respectively.

Figure 12. The 2008–2028 urban sprawl of Hefei and Qingdao in the four predefined directions

5. Discussion

We developed a new CA_BSVC model in this study and compared its simulation results with an early reported CA_GWR model in simulating urban sprawl in Hefei and Qingdao. Compared to GWR, the proposed BSVC method can not only detect the locally varying relationships between urban sprawl and its driving factors, but also avoid the negative influences of outliers and non-constant variance (Hutabarat et al. 2013; LeSage 2004). This study shows that BSVC can improve the simulation accuracy and enhance the modeling successes in less developed regions.

Compared to the simulation results modeled by CA_GWR, the proposed CA_BSVC model achieved higher simulation accuracy, whose FOM increased by more than 2.1% in the calibration and 4.2% in the validation. This indicates that CA_BSVC can effectively simulate urban sprawl and capture its spatial nonstationarity, with improvement over CA_GWR. The transformation probability maps in Figure 5 present substantial spatial differences between CA_GWR and CA_BSVC. The distribution of the transformation probability by BSVC is relatively reasonable (high probability covers fewer areas while lower probability covers more), while the medium-value probability by GWR covers vast of the study area, which suggests the priority of urban sprawl in these areas. In addition, the BSVC transformation probability map in Qingdao more properly described the urban sprawl potential in region A, where the dispersed urban areas were identified by BSVC but not by GWR. We therefore concluded that BSVC has advantages in modeling multi-temporal urban sprawl and capturing the dispersed urban areas.

The simulation results of urban sprawl in 2013 and 2018 can roughly be inferred from the transformation probability maps. The result assessment (c.f. Table 2) indicates that CA_BSVC outperformed CA_GWR in the CorrectSim, MissSim and FalseSim metrics. This could be attributed to the ability of CA_BSVC to accurately identify the spatial variations in the CA parameters that can reflect the effects of driving factors. Although BSVC produced slightly lower FalseSim in Hefei (−0.1%) and Qingdao (−1.3%) compared to GWR, the aggregation of FalseSim can be observed in the black ellipse in Figure 7, indicating that GWR tends to produce clustered FalseSim. The CA_BSVC and CA_GWR models have also been evaluated by TOC (c.f. Figure 10), where the BSVC curve in Hefei is consistently higher than the GWR curve, indicating that CA_BSVC simulated urban sprawl more accurately than CA_GWR in areas with the same proximity to the existing urban areas. This feature has also been proved by FOM in 2018. FOM by BSVC nearly double that of GWR when urban sprawl mainly occurring near the existing urban areas in Hefei.

FOM is considered an important index to delineate the simulation ability of CA models (Pontius et al. 2008a; Tong and Feng 2020). An experiment of 13 modeling applications showed that FOM ranges from 1% to 59% (Pontius et al. 2008a), where only half of the applications have FOMs greater than 20% and 4 of these exceed 30%. The simulation applications in other Chinese cities (Liu et al. 2018) showed that FOM ranged from 19% to 36% (most of them smaller than 25%). A new CA model based on the dynamic land parcel subdivision vector (Yao et al. 2017a) showed an improvement in FOM by 3.7%, compared to a CA model based on the random forest algorithm. A variable weighted LR-based CA model (Shu et al. 2017) showed FOM improvement of 4.44% as compared to an non-weighted LR-based CA model. FOM differences of 1.5% have been reported by the comparison between CA models based on artificial neural networks and support vector regression (Shafizadeh-Moghadam et al. 2017b). Compared to the literature, our CA_BSVC model has good performance not only in the FOM itself (37.2% in Hefei and 56.8% in Qingdao) but also in the accuracy improvement (up to 3.6% in the calibration and 5.7% in the validation).

Rapid urban sprawl has caused negative impacts on environments, energy consumption and sustainability worldwide (Arifwidodo 2012; Artmann, Inostroza, and Fan 2019). In recent years, many countries and local governments advocated different urban development strategies to promote compact, sustainable cities, implying increasing population density and infrastructure utilization through reasonable spatial arrangement (Xu et al. 2020; Yin, Mizokami, and Aikawa 2015). Given the regulations the Chinese government have introduced to restrict the increasing urban areas and control urban sprawl (Wang, Krstikj, and Koura 2017), the future urban scenarios may be dominated by compact patterns in limited areas. During 2013–2018, Hefei has confronted the situation of restrict the city’s extensional expansion. The actual urban sprawl in 2013–2018 was 29 km², which was only 34.9% of that in 2008–2013. Table 4 indicates that the model accuracy of CA_BSVC in 2013–2018 is almost twice that of CA_GWR. This result suggests that CA_BSVC can better simulate urban sprawl in densely urbanized areas and should be suitable for scenario forecast of cities with compact patterns.

The two models were performed using a computer workstation equipped with i7-9750 h CPU, 16 G memory and operating system Windows 10. The computational time was 90 seconds for training the GWR-based transformation rules and was 5026 seconds for BSVC. This shows that GWR’s computational efficiency is much higher than BSVC. However, the BSVC’s running time is acceptable since the performance of modern computers is getting better. During the model implementation for simulating urban sprawl, the computational time of the two models is the same because they used similar input maps: the initial urban map, the final urban map and the transformation probability map. Although both methods can be used to retrieve CA transformation rules, GWR is relatively simple because it is an extension of linear regression (Brunsdon, Fotheringham, and Charlton 1996); whereas, more knowledge is needed to understand the mathematical concepts of BSVC because it relates to Markov Chain Monte Carlo and Gibbs sampling (Wolf, Oshan, and Fotheringham 2018). GWR is widely included in spatial analysis software with a friendly graphical user interface (e.g. ArcGIS and SAM) while BSVC can be performed using R packages and personal code shared on GitHub (Bakar, Kokic, and Jin 2016; Finley and Banerjee 2020).

Limitations exist in our study because the conventional SVCMs usually apply a single scale to the regression subsets, compromising the modeling accuracy of spatial phenomena (Murakami et al. 2019). To overcome the limitation, future work could apply MGWR with various bandwidths for each explanatory variable (Oshan et al. 2019a, 2019b), where the correlation between the MGWR estimates and the true values is stronger than that produced by GWR (Fotheringham, Yang, and Kang 2017). We should also consider the applicability of GWR variates to build CA models, where the variates include GWLR, MGWR and GWTCLR. Future work could also strengthen our CA_BSVC model by using multiscale SVC approaches (Wolf, Oshan, and Fotheringham 2018) to further assess the scale sensitivity of CA models in simulating urban sprawl. Overall, this study improves our approach to CA-based modeling of urban sprawl and provides a useful resource for assessing future scenarios of other rapidly urbanizing cities using CA_GWR and CA_BSVC via the free UrbanCA software package (Feng and Tong 2020).

6. Conclusions

Spatial nonstationarity is an inherent characteristic of land-use change and urban sprawl. Incorporating nonstationarity into the spatially explicit models is a powerful way to simulate and forecast urban sprawl more accurately. We proposed a novel model based on BSVC to simulate the spatially-varying urban dynamics in a rapidly urbanizing area of China. We compared CA_BSVC with CA_GWR to assess its relative performance by simulating urban sprawl at Hefei and Qingdao. Our results demonstrate that: 1) CA_BSVC outperformed CA_GWR in terms of FOM by ~3% during the model calibration and by ~5% during the model validation; 2) CA_BSVC simulated urban sprawl more accurately than CA_GWR in regions with the similar desistance to the existing urban areas, especially in less developed areas; and 3) fewer FalseSim clusters were observed in the BSVC results comparing with GWR. The CA_BSVC model overcomes the outlier and nonstationarity variance effect of the CA_GWR model, hence reducing the clustering of modeling errors and more accurately simulating urban sprawl adjacent to the existing urban areas. Compared to conventional SVCMs, BSVC can substantially improve the modeling accuracy, especially when simulating urban sprawl in densely urbanized areas and forecasting urban scenario with compact patterns.

The proposed CA_BSVC model has stronger simulation ability compared with CA_GWR and strengthens the CA modeling approaches by incorporating spatially-varying transformation rules. The new model improves our understanding of urban dynamics by including spatial nonstationarity. Furthermore, through the combination of SVCM, the Bayesian method and CA, the proposed CA_BSVC model was proved more accurate in predicting urban sprawl in densely developed areas. It can certainly offer better insights into the examination and simulation of spatially nonstationary urban sprawl. The model should also promote the assessment of compact cities excepted in various planning strategies, and provide appropriate scenarios based on the current urban development policies.

Data and codes availability statement

Data and codes used in the research are available at https://figshare.com/s/5436916184e60f533899. Data are shared as an ArcGIS package file (mpk). Data list is available at the figshare site.

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

This study was supported by the National Key R&D Program of China [2018YFB0505000 and 2018YFB0505400], and the National Natural Science Foundation of China [42071371 and 41631178].