A Study of Diffusion Equation-Based Land-Use/Land-Cover Change Simulation

: Simulating and predicting the development and changes in urban land change can provide valuable references for the sustainable development of cities. However, the change process of urban land-use/land-cover is a complex process involving multiple factors and multiple relationships. This dilemma makes it very challenging to accurately simulate the results and to make predictions. In response to this problem, we started with the physical characteristics of the land-use/land-cover change process and constructed a diffusion equation to simulate and predict urban land-use/land-cover changes. The diffusion equation is used to describe the diffusion characteristics of the land-use/land-cover change process, which helps to understand the urban land-use/land-cover change process. The experimental results show that (1) the diffusion equation we constructed can simulate urban land-use/land-cover changes, (2) the simulation process of the model is not limited by the time interval of the time series data itself, and (3) the model only requires one parameter without other constraints.


Introduction
As the primary area of production and living for humans, cities have undergone tremendous changes for human survival and development. More than half of the world's population lives in cities [1,2]. Land provides a basic guarantee for the survival and development of mankind while at the same time restricts human activities. For example, as the population grows, unreasonable urban expansion causes soil degradation, land waste, and damage to the ecological environment. These problems in turn restrict people's related activities and bring about challenges to the sustainable development of cities, especially in housing [3], environment [4,5], climate change [6], infrastructure [7], food security, health, education, decent work, safety, natural resources [8], etc.
In the face of such drastic changes in urban land use, exploring urban land-use/landcover change patterns and studying the factors driving urban land-use/land-cover changes have long-term significance for the sustainable development of cities. In the study of land-use/land-cover change, some scholars have devoted themselves to detecting urban land changes [9,10] and urban boundary expansions [11,12]. Some scholars have focused on exploring the factors driving urban land change, including nonlinear regression and simulations of the characteristics driving urban land change. Simultaneously, some researchers have explored the impact of different factors on urban land-use/land-cover changes in specific regions [13][14][15][16], such as the effect of war on Syria's land change and the impact of Pakistan's urbanization on election politics [17]. In recent years, the sustainable development of cities has also attracted the unanimous attention of scholars [18][19][20][21][22].
As a complex nonlinear spatiotemporal process, many scholars have studied all aspects of the causes of land-use/land-cover type changes, from the overall influencing factors to the single influencing factors. Both are used to solve or determine the mechanism of urban land-use/land-cover type changes. Within this process, different exploration methods have been developed. Using machine learning [23,24] to predict land-use/landcover can describe the nonlinear process better to achieve better simulation results. Still, this method often ignores the spatial relationship in the process of land-use/land-cover changes. Simultaneously, some neural network models, a black-box model, have created difficulties for people to understand and exploit the patterns of land change. In other aspects, the regular pattern of land-use/land-cover change obtained using the statistical learning method is a prediction based on existing data [25]. Although the law of change can be expressed directly with a statistical formula, this formula is weak in describing nonlinearity. The agent-based model (ABM) [26] often defines rules based on the simulation area's specific conditions, and the model is not highly transferable. Meanwhile, the cellular automaton (CA) [27][28][29] model can utilize spatial information by constructing a neighborhood range. It implicitly expresses the decision-making process through spatial transformation rules, and the heterogeneity of decision-making comes only from spatial heterogeneity. It implicitly expresses the decision-making process by constructing spatial transformation rules. Some scholars believe that the most powerful theories are often simplest. Based on this recognition, they start from the spatial form of the city and use the diffusion-limited aggregation (DLA) to simulate the changes in urban space [30]. In recent years, some scholars have combined multi-source data to quantify urban areas on the basis of percolation theory [31].
From some statistical data, we can draw a conclusion. In the process of urban land change, there is a common constant trend, and it is more evident in the shift of impervious urban surfaces; that is, the urban scope presents an overall expansion trend [32]. In the past 30 years, the global urban area has expanded by 80% [2]. From 1985 to 2015, the international urban area increased from 362,700 square kilometers to 653,400 square kilometers. In terms of countries, the United States, China, and India have the largest number of cities, and they have undergone different urbanization tracks [33,34]. Among them, China's urban growth is mainly concentrated in the Beijing-Tianjin-Hebei region, the Yangtze River Delta, and the Pearl River Delta [35]. The growth of these specific urban areas has accounted for more than 60% of the total urban expansion in China. As one of the typical representative cities in China's urbanization, Shenzhen's urbanization process has developed rapidly since the establishment of the Shenzhen Special Economic Zone. In the past 40 years, Shenzhen's permanent resident population has increased from 314,000 to 13.0266 million.
To explore the law of the land-use/land-cover change process, we proposed a landuse/land-cover simulation model based on the diffusion equation which can express this nonlinear relationship in time and space and form a display change relationship formula while using time and space information. Based on Shenzhen city's land-use/land-cover data, a diffusion equation is constructed to simulate the land-use/land-cover change process of Shenzhen city and to predict the future land-use/land-cover change situation of Shenzhen city. Based on the land-use/land-cover data [36], we first compare the simulation effect of the diffusion equation with the CA-Markov model, and the results show the effectiveness and superiority of the diffusion equation. We then further compare the effects of the diffusion coefficient and spatial neighborhood size on the model.

The Study Area and the Data
Shenzhen is located on the southern coast of China and belongs to Guangdong Province, between 113 • 46 -114 • 37 E and 22 • 27 -22 • 52 N. Its jurisdiction is characterized by a narrow east-west length and a north-south diameter. The east-west diameter distance is 157 kilometers, as shown in Figure 1. In 2018, the city's total land area reached 1997.47 square kilometers (excluding islands), and the built-up area is 927.96 square kilometers. Until 2019, Shenzhen had nine administrative districts and one new district. At the end of 2018, Shenzhen had a permanent population of 13 million, with an urbanization rate of 100%. The landform of Shenzhen is dominated by low mountains, flat terraces, and terraced hills. The plain occupies 22.1% of the land area, and the forest coverage rate is 44.6%. The land-use/land-cover classification data of Shenzhen used in this paper is from the 30 m Landsat classification data [36], which is obtained from cloud-free Landsat 5 TM, Landsat 7 ETM+, and Landsat 8 OLI remote-sensing images (WRS path 122, row 44; WRS path 121, row 44) that was downloaded from the US Geological Survey (USGS) website. The data are presented on a land-use/land-cover map of Shenzhen over the past 30 years, with a time interval of 2-4 years, and there are six types of land-use/land-cover involved: forest, cultivated land, water body, grassland, built-up area, and bare land. The AdaBoost classification method based on C4.5 decision tree was used to classify remote sensing images, and then in the verification set with more than 300 random sampling points, the overall classification accuracy reached 90.04% (kappa = 0.89). Thus, it is believed that the information is accurate and reliable. From the data, we can see that, in the past 30 years, the changes in Shenzhen have been highly complex, and different types have shown other trends. Taking the built-up area as an example, the size of the built-up area in Shenzhen City increased in general from 1988 to 2015, which was also mentioned in [36]. It is mentioned that, in the past 27 years, the proportion of Shenzhen's built-up area increased from 5.63% in 1988 to 41.77% in 2015.

The Notations
In actual simulations, we usually use computers to calculate the urban land changes using the diffusion model. It is indispensable first to define the basic unit of action of the diffusion equation. In this article, we adopt a simple method: treating a pixel as a unit of measure both in the CA-Markov model and the diffusion model. Analogous to the definition in the CA-markov model, the basic unit of action is called a cell. Simultaneously, we noticed that no matter what model is used to simulate urban land-use/land-cover changes, the spatial neighborhood's size is essential. In this article, to evaluate the effect of the algorithm under different neighborhood sizes, we considered three different spatial neighborhoods to measure the diffusion model's response. They are the 4 neighborhood, 8 neighborhood, and 25 neighborhood. The relationship between the neighborhoods and the central cell is shown in Figure 2. Finally, to facilitate an understanding of the diffusion model in this article, we used the symbol definitions shown in Table 1 before the specific derivation and provide a detailed explanation of each symbol that is used.

Notations
Definitions The land use category of the grid in row i and column j t n , n ∈ {0, 1, 2, · · · , k} t is the time moment, and k means the biggest year Φ(X i,j , t) The land-use type of the grid in row i and column j at time t p n , n ∈ {0, 1, 2, · · · , k} The domain proportion probability at time n The influence from the neighborhood of the central grid at the moment of t a(X i,j , t) The diffusion coefficient that represents the diffusion capacity of land type X i,j

The Methods
In this section, we explain the algorithm flow in three parts. First, we introduce the idea of the algorithm as a whole and give an overview of the corresponding process. Next, we give the mathematical expression of the algorithm and the related discrete process. Finally, for the key parameters in the diffusion equation, the diffusion coefficient, we give two different solutions.

The Experimental Framework
The general idea of the experiment is to use the historical data of all previous years to predict the land-use/land-cover in the next year. After obtaining the laboratory data, we extract the sequential data blocks from the grid in the spatial neighborhood range and convert them into decomposition coefficients for each unit grid. The next step is to quantify each grid's diffusion effect in the spatial neighborhood on the central grid. Our quantization method predicts the simulation of the diffusion process for each grid in the spatial neighborhood to obtain the influence of the spatial neighborhood grid on the central grid at the next moment. In this way, we obtain the subsequent change in the land-use/land-cover type on each grid in the spatial neighborhood. Based on these changes, we can preliminarily determine the land-use/land-cover type of the central grid at the next moment.
In the diffusion equation, the speed of diffusion is closely related to the substance concentration at the current moment. In our diffusion equation, the proportion of same land-use/land-cover types to the central grid in the current spatial neighborhood is called the "concentration" of the central grid. The specific judgment process is as follows: if the "concentration" of the current center grid increases, then the land-use/land-cover type shows a trend of diffusion and the land-use/land-cover type remains unchanged. If the "concentration" of the current center grid does not change significantly, then the grid has not been affected by the spread of other land-use/land-cover types and the landuse/land-cover type of the center grid remains unchanged. If the "concentration" of the land-use/land-cover type of the current grid decreases, then the land-use/land-cover type of the grid is affected by the spread of other land-use/land-cover types in its neighborhood and there is a tendency to change to different land-use/land-cover types. Therefore, the current grid's land-use/land-cover type becomes the land-use/land-cover type with the highest proportion in the spatial neighborhood or the most apparent land-use/landcover type of the spread. Figure 3 shows the processing flow of our diffusion model for each grid Φ(X i,j , t n ). Under this process, we first simulated Shenzhen's changing pattern from 1988 to 2013 and used the land-use/land-cover classification data of Shenzhen in 2015 as a test. Then, to test the effect of our algorithm in different spatial neighborhoods, we tested the simulation accuracy of each district in Shenzhen in three spatial neighborhoods.

The Mathematical Representation of the Diffusion Process
With the land-use/land-cover data, we observed a phenomenon: the city's area gradually expands outward with human activity; for example, the built-up area increases year by year. Based on this phenomenon, we can assume that, without the intervention of government policies, the changes in urban built-up areas satisfy a specific diffusion equation. Simultaneously, without the promotion of corresponding policies, the changing process of land-use/land-cover type in a certain place is also affected by the surrounding land-use/land-cover types. For example, when the land-use/land-cover type near the area is the same, the land-use/land-cover types in the area are more difficult to change to other land-use/land-cover types. In addition, by distributing more of the same type of land-use/land-cover in the entire spatial neighborhood, the spread of this type also becomes more possible, making the changes in other different land-use/land-cover types nearby consistent with the area. For example, for a fixed spatial range, if the proportion of built-up area in the field is more extensive, with human activities, other land-use/landcover types, such as forests and grasslands, become easier to convert to built-up areas. At this time, for this spatial range, the proportion of built-up areas in the region increases and the proportion of other types decreases. Considering these circumstances, we constructed the following diffusion equation: where a(X, t) means the diffusion coefficient that represents the diffusion capacity of different land types. F(X, t) shows the influence of the neighborhood grid to the central grid.
In the diffusion equation, the diffusion coefficient reflects the diffusion characteristics in historical data. The second derivative term simulates the nonlinear phenomenon in the process of land-use/land-cover change. Regarding the influence of the spatial scope, the more grids in the neighborhood with the same land-use/land-cover type as the central grid, the more beneficial it is to this type's next diffusion process. Otherwise, the diffusion process of the central grid suffers more resistance. It is embodied in the equation that, if the current spatial neighborhood is within the range of the same land type as the central grid, the grid's proportion is p k . The spatial neighborhood's influence on the entire diffusion process is −(1 − p k ).
The above diffusion equation must be solved on a computer, so discretizing the diffusion equation is necessary. In our experiment, we used the second-order precision central difference method to discretize the above equation. The discrete form of the first derivative term is as follows: The second derivative can be expressed as follows: Bringing the above two equations into the diffusion equation, we obtain the discretized diffusion equation form as follows: Simplifying this formula, we can obtain the solution iteration formula as follows: Let r = ∆t ∆X 2 , we can further obtain the following: This formula is the iterative formula for our actual simulation.

The Solving of the Diffusion Coefficient
In the diffusion equation, the diffusion coefficient indicates the diffusion ability of a particular diffusion process. When the diffusion media changes, the diffusion coefficients should be different. In this model, the diffusion coefficient is related to the specific spatial location where the diffusion process occurs and when the diffusion process occurs. This two-dimensional function includes time and space. Simultaneously, the diffusion coefficient describes the diffusion ability of the diffusion process, so in actual experiments, the diffusion coefficient must represent an individual variation ability, that is, at least a first-order derivative related to time and space.
Under these premises, our model's diffusion coefficient is a function of time and space. For different diffusion units, their diffusion coefficients should differ so that the diffusion coefficient and spatial correlation can be satisfied. Second, for each unit, its diffusion coefficient at this time should be related to the historical state so that it can satisfy the time correlation. Finally, considering that the diffusion process' occurrence is also affected by the surrounding spatial neighborhood, the neighborhood range must be considered when calculating the diffusion coefficient of a certain grid.
In our time series data grid, we first extract all of the time series data of the current grid in the fixed spatial neighborhood, and then, we obtain the spatial neighborhood state-change data block of the current grid in the time dimension. For this data block, we calculate the number of neighborhood grids with the same land-use/land-cover type as the central grid at each moment and obtain the concentration for the year. After obtaining the concentration for different years, we can predict subsequent years' ratios using the numerical simulation methods. When the actual situation in the forecast year cannot be obtained, we use a simple polynomial fitting to obtain the slope of change at the next moment. In this way, we can further obtain the slope of the change in the proportion of land types between the forecast year and the previous year, and this slope of evolution is regarded as the diffusion coefficient. At this point, we completed the solution of the diffusion coefficient. This solution pays attention to the time change in the process of land change, which is not only related to the moment of time but also related to the neighborhood of the current grid.
In Figure 4, we show the process of solving the diffusion coefficient when the spatial neighborhood is eight neighborhoods, where t n , n ∈ {0, 1, . . . , k} represents the different moments and p n , n ∈ {0, 1, . . . , k} corresponds to the probabilities at different times. Given the probabilities of all of the moments, we can predict the proportion probability p k+1 of the next moment; by simulating the change in the proportion from p 0 to p k+1 , we can obtain the change slope from p k to p k+1 . Moreover, during the calculation process, it may appear that the proportion of the current grid domain in the time dimension has not changed, so the slope of the change appears to be 0. In the experiment, when the slope is 0, this indicates that the land-use/landcover type in the neighborhood has not changed continuously, but it does not mean that there is no diffusion process. The diffusion term of the entire diffusion equation no longer contributes to the whole diffusion process, which is inconsistent with the actual situation. At this time, the role of this land-use/land-cover type in the spreading process should be more vital than useless. Based on this consideration, we conducted a test on the Luohu District of Shenzhen. The reasons for choosing this district are that Luohu District is one of the old urban areas of Shenzhen and that the types and distribution of land use types are representative. Different experimental results were obtained under different diffusion coefficient settings, and finally, we chose 1 as the final diffusion coefficient. In our experiment, not all of them were 1. The diffusion coefficient is artificially set to 1 only when the land-use/land-cover type in the neighborhood has not changed in all historical data, and the slope is zero.

Results
In the experiment, each pixel of the land-use/land-cover data is regarded as a primary operation unit and the diffusion process is measured against a pixel. In each round of experimental results, the kappa coefficient and overall accuracy were used as evaluation indicators. The whole experiment was divided into three parts. First, we compared the experimental results between the CA-Markov model and the diffusion equation with the same land-use/land-cover data in Shenzhen city. The spatial neighborhood is both the 8 neighborhood. Then, we focused on two important factors in the diffusion model: the diffusion coefficient and the spatial neighborhood range. We successively compare the influence of the diffusion coefficient and the size of the spatial neighborhood on the experimental results using the controlled variable method. In all of the experiments, the kappa coefficient and overall accuracy were used as evaluation indicators.

The Land-Use/Land-Cover Prediction Using the Diffusion Equation
In the process of using the diffusion equation to predict land-use/land-cover changes, each pixel was used as a diffusion unit. The land change process for each diffusion unit is related to the changing pattern over time and the land-use/land-cover type in the spatial neighborhood. In this part, the standard 8 neighborhood is adopted as the spatial neighborhood range in the diffusion model. Then, we obtained the diffusion coefficient by calculating the changing trend over time, from 1988 to 2013, and we determined the landuse/land-cover type next by calculating the influence of different land types on the current unit's spatial neighborhood. After determining the diffusion coefficient, by combining the spatial neighborhood and the influence of other grids in the neighborhood on the central cell using the diffusion equation, we obtained the concentration of the central grid next. If the current grid concentration increases or does not change significantly, then the land type of the current grid in the field represents the trend of diffusion and the current type remains unchanged; if the current grid concentration decreases, then the land type of the current grid becomes the land-use/land-cover type with the highest proportion in the field. Under this judgment condition, we tested Shenzhen as a whole.
The range of the spatial neighborhood used in the test process is the 8 neighborhood. Under the same preconditions and judgment conditions, we also obtained the overall forecast of Shenzhen in 2015, which is shown in Figure 5. Compared with the original 2015 land-use/land-cover classification data, the Kappa of the predict result is 0.6760, the overall classification accuracy is 0.7995, the classification confusion matrix is shown in Table 2.  In the comparison experiment, we used the CA-Markov model provided in the IDRISI software as a comparison. The state transition matrix of Shenzhen from 2011 to 2013 and the land-use/land-cover of Shenzhen in 2013 were input in the software, and the spatial neighborhood size was set as the 8 neighborhood. Finally, we created the land-use/landcover situation in Shenzhen in 2015, Figure 6, and the Kappa coefficient of the experimental result is 0.6741. Then, we tested each district of Shenzhen separately. For the district of Shenzhen, the kappa coefficient and accuracy are shown in Table 3.

The Influence of Different Neighborhoods
The change between the land-use/land-cover type is affected by the surrounding areas, and the prediction effect of the different regions on urban land-use/land-cover changes differs as well. We also tested the impact of varying neighborhood sizes on land-use/landcover prediction in three different data sets. In this paper, we selected three districts of Shenzhen as representative: Nanshan District, Luohu District, and Yantian District, which are representatives of the old community of Shenzhen City. When evaluating the impact of different spatial neighborhood sizes on the experiment, we chose three classic neighborhood ranges, and the specific definitions can be seen in Figure 2.
There is a common trend among Nanshan District, Luohu District, and Yantian District under different neighborhood conditions. The more neighborhoods used in the diffusion process, the higher the accuracy of the results. However, in the 8 neighborhood 8 and 25 neighborhood, the difference was not significant. The land-use/land-cover prediction results of different regions are shown in Figures 7-9. To show this difference more intuitively, we list the experimental results of various districts in Shenzhen in different neighborhoods in Table 4.

The Influence of Diffusion Coefficient
In the diffusion equation, the diffusion coefficient determines the size of the diffusion capacity. By using different diffusion coefficients, some phenomena can be better observed. First, we used the change of slope in the time dimension, that is, the change in the proportion of the same land type in the neighborhood, as the central grid to predict the change slope from 2013 to 2015. Taking this slope value as the diffusion coefficient, a set of experimental results was obtained in various districts of Shenzhen. Then, we adopted the state transition matrix used in the CA-Markov model, taking the diagonal elements of the state transition matrix as the diffusion coefficient. Specifically, the diagonal transition probability of the land-use/land-cover type corresponding to the center grid is the diffusion coefficient of the center grid. Thus, we obtained the second set of experimental results. In both experimental cases, the spatial neighborhoods used are 8 neighborhood. The comparison of the experimental results is shown in Table 5 below.
In the experiment, the state transition matrix represents the overall land transformation of the region, and the change slope of the time dimension only represents the transformation within the spatial neighborhood of the central grid. By introducing a state transition matrix with global transition information, there are usually better simulation results in some old urban areas. In some new urban areas or urban areas that undergo surface changes that are more drastic, it is more accurate to use the time change slope in the neighborhood as the result of the diffusion coefficient simulation.

Discussion
In our model, the time interval of time series data does not limit the algorithm, and it can simulate the land-use/land-cover in subsequent years. For example, using the twotime series data of 2000 and 2005, we can handle the land-use/land-cover situation in 2010, the land-use/land-cover situation in 2011, or later years. Second, the model parameters of our method are simple. Only one parameter is needed to simulate changes in urban land-use/land-cover. This brings great convenience to the use of the model.
In this study, we constructed a simple diffusion equation model to simulate landuse/land-cover changes. It achieves the same prediction effect as the classic CA-Markov model and avoids strict requirements on the time interval of time series data, but there are still many areas for improvement. We think that some auxiliary data such as road network data, land policy data, etc., improve the prediction results. These extra data consider different driving factors for future land-use/land-cover changes, improving our prediction results. First, for the diffusion coefficient of the model, we can consider introducing multi-source data to improve the situation where the diffusion coefficient is 0. Second, When calculating the influence of the neighborhood, the overall distribution characteristics can also be added to the model. For example, when forecasting the Nanshan District of Shenzhen city, the overall change of Shenzhen can be added to our model. In this way, the local change information and the overall change information were both considered.

Conclusions
This research proposed a better simulation model of land-use/land-cover change and conducted a simulation result of Shenzhen. First, by selecting the same 8 neighborhood and by comparing it with the CA-Markov model, our proposed method performs better in terms of kappa coefficient and overall accuracy. This also confirms that the diffusion equation model which we proposed can simulate urban land-use/land-cover changes to a certain extent. Then, we tested the differences in land-use/land-cover changes in different areas in different districts of Shenzhen, and the results showed that the results were more accurate when more neighborhoods were considered. However, the accuracy improvement from the 4 neighborhood to the 8 neighborhood is more obvious, while the difference between the 8 neighborhood and the 25 neighborhood is not significant. From this point of view, our algorithm has a certain degree of robustness, and it can resist the adverse effects of different neighborhood ranges. Finally, we started from the state of land-use/land-cover itself, did not introduce other artificial assumptions of transformation conditions, and only simulated land-use/land-cover change based on the influence from surrounding areas during land-use/land-cover change, avoiding some errors caused by human factors. Unlike other algorithms that require continuous time series data with the same time interval when predicting land-use/land-cover change, this method can use time series data of any time span, which provides great convenience for subsequent land-use/land-cover change prediction.