Temporal-Spatial Simulation of Landscape Variations Using Combined Model of Markov Chain and Automated Cell

Classification and preparation of land use map are the most commonly used methods in remote sensing data. Numerous advanced classification methods have been developed in a recent year. Among those methods, one can name support vector machine. Identifying land cover changes can play an important role in future decisions of regional managers. Therefore, the present study aimed at employing CA–Markov model to present a powerful but simple model for simulation and prediction of Yazd city. Accordingly, first, satellite images of 2000, 2005, 2010, and 2016 were classified by SVM method, and then land use maps of Yazd city were extracted. The results of the detection of changes indicated the expansion of residential areas and the reduction of dry land and vegetation during a 16-year period. Prediction of land use development in 2040 using CA–Markov model with a kappa coefficient of 84.43 shows the high accuracy of this model. Besides, the urban land use area shows 18,432 hectares increase in 2040 compared to 2010, while the area of dry land uses (− 15,570) rocky areas (− 186) and vegetation (− 2658) hectares decrease compared to 2010. Moreover, the results of CA–Markov model indicated the development of Yazd city by degradation of vegetation.


Introduction
The nature of ground is not constant and always changes. Detecting land use changes and reviewing its trends are one of the basic requirements for managing and evaluating natural resources (Fathizad et al. 2017;Hamidy et al. 2016;Alipur et al. 2016). Monitoring land use changes in time intervals is achieved more accurately through a remote sensing technique in a shorter time period and a lower cost Aliabad and Shojaei 2019;Achu et al. 2019;Kermani and Rohrbach 2018). Remote sensing is of central importance because of certain features such as wide visibility, integration, the use of various electromagnetic energy components to record phenomena, short return periods, the possibility of using hardware and software, low cost and faster inspection (Hamidy et al. 2016;Kermani and Rohrbach 2018;Shojaei et al. 2018;Nasiri et al. 2017;Reinsma 2019). Today, evaluation of ground changes (the conversion of forest, pasture, agricultural lands, and gardens to residential areas and urban facilities) through change detection process using multi-temporal remote sensing imagery is an alternative method for costly, time-consuming and often low-priced methods (Arabi Aliabad et al. 2018). Iranmehr (2014) employed TM, ETM + , and OLI Landsat images related to years of 1985, 2003, and 2013 to detect and predict the land use and cover changes around Zayanderud. First, all pre-processed images were reconstructed and then hybrid classification method was applied. Finally, land use and cover maps were prepared for each year. After estimating the accuracy of the resulting maps using error matrix, the comparison method after classification was used to detect changes. The results showed that the highest urban development occurred between 1985 and 2003 with an average of 952 hectares per year and the least growth occurred between 2003 and 2013 with an area of about 552 hectares each year. Arkhi (2015) employed an object-oriented approach to study land use/cover changes in Abdanan. The results indicated the decreasing trend of medium and good rangeland in the area which shows the degradation of the region through replacing medium and good rangeland by poor and dry land. Makroni and Soltaniyan (2016) prepared the land use map in years of 2003 and 2014 using maximum likelihood method. The results of the monitoring changes showed that in the course of the study, the area of residential and rangeland use has been increased and the area of wetlands and agricultural lands has been reduced. Haack and Rafter (2006) examined changes in urban development by preparation and evaluation of land use maps using satellite imagery and detection techniques between the years 1978 and 2000 in the Katmandu Nepal valley. Petropoulos et al. (2012) developed a land use map using both object-oriented and support vector machine (SVM) methods and concluded the suitability of both methods in the preparation of land use map; however, SVM method has higher accuracy and kappa coefficient than object-oriented method. Rizk Hegazy and Rashed (2015) evaluated the urban growth and land use changes of Dakahlia province in Egypt using remote sensing method and GIS between the years of 1985 and 2010. The results showed 30% increase in urban growth and 33% decrease in agricultural lands. Fan et al. (2008) employed TM images of 1977 and ETM + images of 2003 to evaluate the land use and Fig. 1 Location of the study area in Province and Iran cover changes of Pearl River between 1977 and 2003 and to predict the changes for 10 years. They classified images using supervised classification with maximum likelihood method. Then, they were evaluated the accuracy of obtained maps by random sampling and selecting 690 samples in field observation. The results showed 91.5 and 86.8% accuracy for classified maps of 1977 and 2003, respectively. After classification of land use and cover changes between 1977 and 2003, a comparison method was performed to identify the urban expansion trend and degradation of agricultural lands between 1990 and 2000. Finally, urban expansion was predicted for 2008 and 2013 using Markov chain and CA-Markov. In the present study, land use changes of Yazd city was performed using SVM method and prediction of these changes for 2040 was done by CA-Markov model.

Study Area
Yazd city with 2397 km 2 is located between the latitude of 29° 52′-33° 27′ N and longitude 52° 55′-56° 37′ E. Due to the location of Yazd city in desert margin, it has dry climate. Also, due to the lack of raining, no vegetation can be seen in this area. Figure 1 shows the location of study area.

Land Use Map
In the present study, sensor data of TM (2000 and2010), ETM + (2005), and OLI (2016) were used to extract land use map (a 16-year period). The satellite images of May and June was used because of vegetation cover in spring. These images had a proper condition in terms of the lack of cloud cover and rainfall. Detailed data are presented in Table 1.
To prepare the data for processing, after preparation of satellite data, geometric correction operations and dereferencing were performed using road network vector layer and satellite imagery of the area. Spectral correction of images also was done to enhance the phenomenon and increase the quality of images and remove the adverse effects of light and atmosphere on them. Then, the band 4-3-2 false-color composite (RGB) was formed using the correlation between bands and classification was conducted by SVM method.
According to the objective of the present research and vegetation type in the area, four categories including residential, dry, vegetation, and cliff areas were identified and classified. To ensure the classification, a comparison was made between land use map and field observations. In the present study, random sampling method was employed to classify  data. Samples were selected based on land use map and field observation randomly.

Estimation of Produced Map Accuracy
To evaluate the accuracy of the classification, the error matrix, the statistical parameters of the overall accuracy, producer accuracy, user accuracy, kappa coefficient, total kappa, user error, producer error, and total error were calculated, and then the majority filter was applied to obtain a uniform image and the removal of dispersed pixels over the images of the classification (Fig. 2).

Support Vector Machine (SVM)
Statistical learning theory-based SVM method is a non-parametric supervised statistical method (Pao 1989), in which the samples that form the classes' boundaries is obtained using all bands and an optimization algorithm and then an optimal linear decision-making boundary is calculated to separate the classes. These samples are support vectors (Keshavarz and Yazdi 2005) (Fig. 3). As shown in Fig. 2, binary samples are called support vectors, and the optimal margin method is used to calculate the decision boundary of two distinct classes. To calculate the margin, all class + 1 samples are placed on one side and all class − 1 samples on the other side of the boundary. Besides, to calculate decision-making boundary, the distance between the closest training samples of both classes in the vertical direction of the decision boundary should be maximized (Keshavarz and Yazdi 2005). By parallel expansion of decision boundary from both sides through two lines to pass through the closest samples of two classes, optimal hyperplane boundary is formed, which a boundary with the maximum is spacing between two classes. Those two parallel lines are called marginal hyperplanes. For pixels outside the marginal hyperplanes, Eq. 1 is applied.
For class pixels, where x is a decision boundary, W T threshold, W is a perpendicular n-dimensional vector to the decision boundary (Goodarzi Mehr et al. 2012).

Modeling Using CA-Markov model
Markov chain and CA are both discrete dynamic models in time and position. Producing no geographic perceptions is the inherent problem of the Markov chain. Although the probability of conversion on each basic node may be accurate, it lacks stochastic space distribution within each land use group, i.e., there is no spatial component in the modeling output. As a result, CA-Markov is a proper method for spatial-temporal dynamic modeling of use and cover changes and remote sensing method and GIS data can effectively contribute to this model (Li and Reynolds 2000). In CA-Markov model, Markov chain converts the temporal changes between land cover/use class based on probabilities, while spatial changes are controlled by local rules determined by CA (Eastman 2006). To monitor the spatial patterns of land cover and land use in Yazd, Markov chain model was first used to create the transmission probability matrix and the area of the transmission of various land cover classes, which provides more detailed information on the transfer between classes in different types of land cover. Subsequently, the prediction of land cover changes using the CA-Markov model and based on the probability of a transmission derived from the Markov chain analysis, which could give us a simulation map for land cover in the future. In this research, first, the necessary changes and corrections of images were performed in the ArcGIS 10.2 software, and then models of IDRISI software were used to simulate land use changes in Yazd area.

Image Classification
In the present study, according to a review of previous studies, SVM method was employed to classification. Based on the familiarity with the study area also field observation four class including vegetation, dry, cliff, and residential land use were identified. The images of 2000, 2005, 2010, and 2016 can be seen in Figs. 4, 5, 6, and 7, respectively.

Validation
The necessity of using any kind of satellite image processedland use map is its accuracy. An error matrix was used to assess the accuracy of land use/land cover maps in 2016 by ground control points (100 points). Besides, false-color composite image interpretation, land use map provided by the Natural Resources Administration of Yazd and Google Earth (areas that have not changed over time) (255 truth ground points for each year) and error matrix were used for accuracy of images classification (maps related to years 2000, 2005, and 2010). The results showed the highest accuracy for a map of 2016 derived from the processing of OLI sensors images with 0.9126 kappa coefficient, and 93.57 total accuracy, which followed by a map of 2005 derived from the processing of ETM + sensors images with 0.8877 kappa coefficient, and 91.57 total accuracy. Table 2 shows producer and user's accuracy.
The results of area measurement of land use in the study area and mentioned time intervals are presented in Table 3.
The results of area measurement indicated the increased area of vegetation, residential, and cliff area; and decreased area of drylands in the first 5 years. After 2005, except for   Then, in this research, the horizontal alignment method was used to detect changes using the comparison method after classification. An assessment of land use change between 2000 and 2016 showed that an increase in the area occurred in the built areas (13,275 ha), as well as the decrease in the area in the dry (9965 ha), vegetation (3294 ha), and rocky (16 hectares) areas.
As shown in Table 4, the largest increase in the area occurred in the built areas. During 16 years, 13,275 hectares of land uses (10,310 hectares of drylands and 2654.55 hectares of vegetation) have been converted to residential areas. The results of land use changed during the studied years are presented as the following diagram.
The results showed significant development in Yazd city because of an increase in immigration. Therefore,   fundamental change occurred in the city's inner texture and structure. Also, an increase in immigration resulted in changes in vegetation boundaries (agricultural and dry land). In addition, vegetation was reduced during the period, and most of it was captured by the developed regions.

Modeling the Land Use Changes Using CA-Markov model
The land cover maps of 2000 and 2005 were used to predict the land cover of 2010. This is done to determine the validity and efficiency of the model in prediction and modeling by mapping from the model and the map of the actual land cover. The result of the implementation phase of the Markov chain analysis, a transition likelihood matrix is a transmitted area matrix. Table 5 shows the likelihood transition matrix of land use between 2000 and 2005.
The results of evaluating the prediction accuracy with the Markov chain model using the 2010 land use map are presented in Table 6.
To validate the Markov chain model, simulated land use map by the model and real map resulting from 2010 classification were compared. Kappa coefficient derived from error matrix between the mentioned maps was obtained to be 84.43%. Besides, the difference between the estimated and processed area of 2010 is shown in Table 6.
Then, data including land use/cover map of 2005 as the base map and 2000-2005 transformation matrix of areas were combined in IDRISI Selva software using the spatial function of CA to predict land use map of 2010. Kappa coefficient derived from error matrix between the modeling-based map and land use map of 2010 was 84.43%. Besides, the difference between the predicted and processed area of 2010 is shown in Table 6. The map obtained from validation of CA-Markov model for 2010 is presented in Fig. 8.

Modeling the Land Use Changes Using CA-Markov model for 2040
The land cover maps of 2000 and 2016 were used to predict the land cover of 2040 (Fig. 9). Table 7 shows the matrix of transition probability for land use changes during the statistical years of 2010-2040 using Markov chain method.

Discussion and Conclusion
Classification of land cover using satellite images is one of the most important applications of remote sensing and most of the algorithms have been developed for this purpose. In the present study, SVM algorithm was employed as a classification method of satellite images. Besides, in this research, CA-Markov model was used to evaluate the trend of land use changes and predict this trend for the arid area of Yazd city. Based on the results of evaluating the land use changes between 2000 and 2016, a significant increase in the area of residential area (13,275 ha) and a decrease in the area of drylands (9965 ha 1989-1975 and 1990-2001. Finally, with the help of simulated models, it was concluded that some types of trees whose developmental trend is expected with significant decline requires a proper forestry management. Arabi Aliabad et al. (2018) investigated the land use changes and the way of using lands by remote sensing data of 1973 and 2010. Besides, they simulated the urban growth of Mumbai city for the 10-year period from 2020 to 2030. The results of detecting changes showed that the highest growth rates occurred between 1973 and 1990, with city growth rates of 142 percent. In contrast, the city growth rate between 1990 and 2001 and between 2001 and 2010 were 40 and 38%. The most affected lands were dry lands. The results of the model showed that urban growth will continue in the next 2 decades, so that in 2020 and 2030, the city will have 26% and 12% growth rate (Most of these changes will be along the north). Iranmehr (2014) used Landsat images of TM, ETM + , and OLI related to years of 1985, 2003, and 2013 to identify and predict the land use and cover around Zayanderud. First, all pre-processed images were reconstructed and then hybrid classification method was applied. Finally, land use and cover maps were prepared for each year. After estimating the accuracy of the resulting maps using error matrix, the comparison method after classification was used to detect changes. The results showed that the most urban development occurred between 1985 and 2003 with an average of 952 hectares per year and the least growth occurred between 2003 and 2013 with an area of about 552 hectares each year. Moreover, a combined model of CA-Markov was employed to predict the city development of 2013. Then, the obtained map was compared with other maps related to 2013 to evaluate the capability of a model for predicting changes. Finally, urban areas map was predicted for 2040. The results showed 85,247 hectares growth in the urban area.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat iveco mmons .org/licen ses/by/4.0/.