Delineating the future boundaries of urban development using the FLUS Model: A Case Study of Zhaoyuan City, China

Zhaoyuan City in Shandong Province, China was used as an example to explore the driving factors, constraints, and model parameterization using the FLUS-UGB model. Using observed data from 2018 to determine the model accuracy, the Kappa coefficient was found to be 0.8132, the overall accuracy was 0.8924, and the reliability was high. The forecast results show that, from 2018 to 2030, the urban land area growth rate of Zhaoyuan City was 11.40%, with the central and northwestern regions showing higher growth rates. The urban land area growth rate of Jinling Town was the highest at 42.53%, while Mengzhi and Wenquan Streets and Xinzhuang Town’s urban land area growth rate exceeded 10%. Other towns’ urban land area growth rate was less than 10%.


Introduction
With the development of cities and the acceleration of urbanization and consequent population growth, there has been an expansion in urban land use and a corresponding reduction in cultivated and ecological land uses. In order to protect ecosystem services and resources from cultivated land and improve the efficiency of urban land use, it is necessary to delineate urban development boundaries [1][2].
Recent research on the boundaries of urban development can be divided into three categories. The first is based on evaluation of available land resources and population growth to predict changes in the range of urban area [3][4]. The second is based on cellular automata (CA) and traditional urban growth models such as multi-agent systems (MAS) to simulate land use changes [5]. The third category integrates the above two methods and combines them with constraints from mathematical models such as the CLUE-S, SLEUTH, MCE-CA , and Constraint CA models. [6][7][8][9][10][11].
The FLUS model is a future land use change prediction model based on combining CA and neural networks [12][13]. It mines the mapped relationships between driving factors and urban expansion, and shows high accuracy in its ability to simulate complex land use changes [14][15], making it suitable for the study of land use changes in small-scale cities and towns [16][17]. Recently, research on the boundaries of urban development has mainly focused on large and medium-sized cities [18]. This study uses Zhaoyuan City as an example to study future potential land use and layouts of the space in order to provide a reference for the land use planning work.

Overview of the study area
Zhaoyuan City is located in northeastern Shandong Province, between 120°08′～120°38′east and 37°05′～37°33′north. It is adjacent to Qixia City to the east, Laizhou City to the west, Laiyang and Laixi cities to the south, Longkou City to the north, and the Bohai Sea to the northwest. The city has a total area of 1,432.32km 2 and a coastline of 13.5km. It is rich in mineral resources and lacks freshwater resources [19]. With the rapid development of cities and towns in the region, there is a growing tension between expanding urbanization and protecting the natural environment and associated ecosystem services and cultivated land resources.

Data source and preprocessing
Land use data from Zhaoyuan City used in this study was from 2009-2018 from the land use survey change database, which was then reclassified into seven types of land: cultivated, woodland, grassland, water area and tidal flats, urban land, other construction-related land, and unused land. DEM data at a 30m resolution was downloaded from the Geospatial Data Cloud (http://www.gscloud.cn/) and processed and analyzed by ArcGIS 10.2 to obtain slopes and topographic undulation data. Red line data indicating environmentally protected sea and land area as well as permanent basic farmland data were collected from the Zhaoyuan Municipal Bureau of Natural Resources. Road and river data were extracted from land use data. Point of Interest (POI) data was derived from Baidu API to get the distances to the school, the scenic park, and the key hospital.

FLUS model
Due to its convenience and high accuracy, the FLUS model has been widely used in research on urban expansion and land use change. The FLUS model is mainly divided into two major modules: "Calculation of Probability of Suitability Based on Neural Network (ANN)" and "Cellular Automata Based on Adaptive Inertial Mechanism".
The suitability probability calculation process involves identifying the driving factors and the land use data at the beginning of the research period. The ANN model then fits the complex relationship between them through a large number of calculations and loop iterations. The specific formula for the calculation process is as follows: where: sp (p, k, t) represents the probability that the land type k on grid p appears at time t. ω j,k is the weight value from the hidden layer to the output layer; ω i, j is the weight value from the input layer to the hidden layer; and x i (p, t) is the relationship function between the variable i in the input layer and the p grid at time t.
Adaptive inertia mechanism cellular automata calculates the total probability of each grid square by determining the adaptive inertia coefficient, combining the calculated suitability probability conversion cost matrix and neighborhood factor parameters, and uses the roulette model to determine whether land-type conversion occurs at the next moment, and finally realizes a rationalized land-use simulation through loop iteration. The specific formula is: where: , is the total probability of the grid p being converted to land type k at time is the total number of grids of land use type k after the last iteration of the N×N Moore neighborhood window, in this paper N=3; ωk is the neighborhood effect weight of each land type; is the adaptive inertia coefficient of land type k at time t; → is the conversion cost from land type c to land type k; and is the difference between the number of grids of land use type k and the number of demand at time t-1.

Markov chain
Markov chain theory is used to predict the total number of pixels in the future. Based on the land use data from 2009 to 2018, this method predicts the number of pixels of land demand in each category in 2030. The specific formula is as follows: ∈[0,1) (4) where: S t and S t+1 respectively represent the land use status at time t and t+1; and P ij represents the transition probability matrix at time t.

Driving factors
The factors affecting the expansion of urban constructed land include natural and social factors [20][21]. The following driving factors were selected: elevation, slope, topographic undulation, distance from waters and beaches, distance from a highway, distance from a transportation hub, distance from a school, distance from scenic parks, and distance from large hospitals.

Constraints
Combining the principles of ecological priority and cultivated land protection, the environmental protection red line, permanent basic farmland and tidal flats were selected as areas unavailable for urban expansion (Figure 1). The land use type did not change during the simulation conversion, and the remaining areas were classified as convertible areas.

Cost matrix
The cost matrix reflects the difficulty of conversion between various types of land. It is represented by a value between 0 and 1. The closer the value is to 1, the easier it is to convert, and 0 indicates that it cannot be converted. Combining the 2009-2018 land use type conversion matrix with expert advice, the cost matrix predicting land use change in 2030 was determined. The detailed data is shown in Tables 1 and 2.

Neighborhood factor
The neighborhood weight represents the self-expansion ability of various types of land as driven by external factors. Existing studies used the total area change of the land type as the neighborhood weight parameter [21], however, this method does not consider the area base of each area. For land-types that have a small base but are fast-growing, it is difficult to accurately reflect their expansion ability. This study uses the dimensionless value of the total area change rate of the land type as the neighborhood parameter, with this specific calculation formula: * (6) where: X* is the dimensionless value of the change rate of the total area of the land; X is the change rate of the total area of the land; max is the maximum value of the change rate of each locality; min is the minimum value of the change rate of each locality; and X* is used as the weight parameter for the neighborhood of the land type. Table 3 displays the neighborhood weight X* values. From the value of X*, other constructed land, urban land and other land are shown to have have strong expansion capacity.

Number of target pixels
Markov Chains were used to predict the demand for each land type in 2018 and 2030 by using the forecast results in 2018 as the accuracy test data, with the overall accuracy reaching 99.64%. The specific data are shown in Table 4.

Simulation accuracy test
In the GeoSOS-FLUS software platform, the land use data and driving factors in 2009 were input in order to conduct ANN sample training. The obtained suitability probability data and related parameters were then input into the CA module to simulate land use in 2018. The results met the consistency requirements, with an overall accuracy of 0.8924, a Kappa coefficient of 0.8132, and a Kappa coefficient of 0.81-1 which is considered almost identical.

Boundaries of urban development
The land use data in 2018 was used as the base period data to simulate land use in 2030, and the area of urban land used in 2030 was predicted as the urban development boundary (Figure 2). The forecast results show that from 2018 to 2030, the area of urban land, other constructed land, and other land will increase, the area of arable land, woodland, and grassland will decrease, the conversion rate of waters and beaches will be low, and the probability of new urbanized land around the original urban land will increase. These findings are in line with the law of historical change.  From 2018 to 2030, Zhaoyuan City's urban land area growth rate was 11.40%. The highest growth rate of urban land area was 42.53% in Jinling Town, and the growth rates of urban land area in Mengzhi Street, Wenquan Street, and Xinzhuang Town all exceeded 10%. The growth rate of urban land area in other towns was less than 10%, and Xiadian Town showed negative growth.

Conclusion
Based on the FLUS model, this study predicts the urban development boundary of Zhaoyuan City in 2030. The study uses the principle of "ecology first, cultivated land coordination", integrates spatial constraints, explores the model parameterization method, proposes a new method to quantify neighborhood weights, and optimizes the forecasting model. The results show that the urban land area growth rate of Mengzhi Street, Wenquan Street, Quanshan Street, and Daqinjia Street in the central city is higher than that of other towns, ranging from 7% to 18%. Outside the central city, the growth rate of urban land area in Jinling and Xinzhuang Towns was higher, 42.53% and 13.04%, respectively. The growth rate of urban land area in Jinling Town was the highest among all towns. The central part of Zhaoyuan City is close to the central city area, so daily life for community members may be more convenient, with more obvious areas of urban expansion. The northwest part of Zhaoyuan City is close to the coast, with gentle slopes, less ecological land distribution, and is close to the the main river, which is suitable for urban land expansion. Therefore, the urban development boundary has a larger area in the central and northwestern parts of Zhaoyuan City.