Spatiotemporal Evolution Analysis and Future Scenario Prediction of Rocky Desertification in a Subtropical Karst Region

: Landscape change is a dynamic feature of landscape structure and function over time which is usually affected by natural and human factors. The evolution of rocky desertiﬁcation is a typical landscape change that directly affects ecological environment governance and sustainable development. Guizhou is one of the most typical subtropical karst landform areas in the world. Its special karst rocky desertiﬁcation phenomenon is an important factor affecting the ecological environment and limiting sustainable development. In this paper, remote sensing imagery and machine learning methods are utilized to model and analyze the spatiotemporal variation of rocky desertiﬁcation in Guizhou. Based on an improved CA-Markov model, rocky desertiﬁcation scenarios in the next 30 years are predicted, providing data support for exploration of the evolution rule of rocky desertiﬁcation in subtropical karst areas and for effective management. The speciﬁc results are as follows: (1) Based on the dynamic degree, transfer matrix, evolution intensity, and speed, the temporal and spatial evolution of rocky desertiﬁcation in Guizhou from 2001 to 2020 was analyzed. It was found that the proportion of no rocky desertiﬁcation (NRD) areas increased from 48.86% to 63.53% over this period. Potential rocky desertiﬁcation (PRD), light rocky desertiﬁcation (LRD), middle rocky desertiﬁcation (MRD), and severe rocky desertiﬁcation (SRD) continued to improve, with the improvement showing an accelerating trend after 2010. (2) An improved CA-Markov model was used to predict the future rocky desertiﬁcation scenario; compared to the traditional CA-Markov model, the Lee–Sallee index increased from 0.681 to 0.723, and ﬁgure of merit (FOM) increased from 0.459 to 0.530. The conclusions of this paper are as follows: (1) From 2001 to 2020, the evolution speed of PRD was the fastest, while that of SRD was the slowest. Rocky desertiﬁcation control should not only focus on areas with serious rocky desertiﬁcation, but also prevent transformation from NRD to PRD. (2) Rocky desertiﬁcation will continue to improve over the next 30 years. Possible deterioration areas are concentrated in high-altitude areas, such as the south of Bijie and the east of Liupanshui.


Introduction
As an important part of global change, the evolution of rocky desertification has a significant impact on the quality of regional ecological environments. Rocky desertification landscape change is a dynamic feature of landscape structure and function over time which is usually affected by terrain, climate, and human factors. Karst rocky desertification is a surface landscape change similar to desertification, characterized by such features as vegetation degradation, soil erosion, and exposed bedrock, against the background of a fragile karst geology and tropical-subtropical humid or semihumid climate [1][2][3]. It comprises a dynamic land degradation process. Guizhou Province in China is one (2) Most of the research periods selected by the existing research institutes have focused on the 1990s through to the beginning of the 21st century, and there are no research results for the past five years, and especially less analysis and simulation research focused on the evolution trend in future scenarios, which makes it difficult to provide effective data support for present and future rocky desertification control. Using the 2001-2020 rocky desertification grade mapping and suitability factors, this study utilizes quantitative statistics and analysis of the temporal and spatial changes of rocky desertification from a dynamic perspective, as well as rocky desertification change intensity and speed, thus revealing the historical evolution law of rocky desertification in Guizhou and improving the CA-Markov transfer matrix to predict the rocky desertification scenario in the next 30 years. Furthermore, we set up three different governance scenarios to analyze the temporal and spatial evolution pattern of rocky desertification, in order to explore the evolution trend of rocky desertification in the future and to provide data support for the accurate management of rocky desertification.

Research Area
Guizhou Province in China is one of the most typical areas featuring subtropical karst landform development. It is located in the hinterland of Southwest China (latitude 24 • 37 -29 • 13 ; longitude 103 • 36 -109 • 35 ). It spans about 595 km from east to west, 509 km from north to south, and has a total area of 176,200 km 2 [44]. The research area was shown as Figure 1.

Data
Using the MODIS data set and National Forest Continuous Inventory data (NFCI) on the Google Earth Engine (GEE) platform, and referring to previous research methods, we constructed rocky desertification maps for different periods [48][49][50]. In the National Forest Continuous Inventory data (NFCI) of Guizhou, rocky desertification is divided into five categories, coded as 00, 10, 21, 22, and 23-24, representing NRD, PRD, LRD, MRD, and Guizhou Province is an area where carbonate is widely distributed and karst is strongly developed. The landform in Guizhou is divided into four basic types: plateau, mountain, hill, and basin. It is mostly mountain, followed by hills, together comprising 92.5% of the province [45]. The soil types in the province are complex and diverse, mainly including red soil, yellow soil, and yellow-brown soil, of which yellow soil has the largest area and is mainly distributed in the central region, representing the main soil type in the province. The annual average temperature in Guizhou is −1 to 25 • C, the annual sunshine is 900-2600 h, the annual precipitation is 500-2500 mm, and the perennial relative humidity is more than 70% [7,46]. Affected by climate, soil, and mountainous terrain, the vegetation types in the province are diverse. The central and northern part is dominated by middle subtropical evergreen broad-leaved forest, the southern part is subtropical evergreen broad-leaved forest with tropical composition, the middle eastern part is humid forest, and the western part is semihumid forest. Cold temperate subalpine coniferous forest is distributed in high-altitude areas, and hidden karst evergreen deciduous broad-leaved mixed forest and secondary deciduous broad-leaved forest are distributed in limestone and dolomite mountains [47].

Data
Using the MODIS data set and National Forest Continuous Inventory data (NFCI) on the Google Earth Engine (GEE) platform, and referring to previous research methods, we constructed rocky desertification maps for different periods [48][49][50]. In the National Forest Continuous Inventory data (NFCI) of Guizhou, rocky desertification is divided into five categories, coded as 00, 10, 21, 22, and 23-24, representing NRD, PRD, LRD, MRD, and SRD, respectively [51]

Data
Using the MODIS data set and National Forest Continuous Inventory data (NFCI) on the Google Earth Engine (GEE) platform, and referring to previous research methods, we constructed rocky desertification maps for different periods [48][49][50]. In the National Forest Continuous Inventory data (NFCI) of Guizhou, rocky desertification is divided into five categories, coded as 00, 10, 21, 22, and 23-24, representing NRD, PRD, LRD, MRD, and SRD, respectively [51]. We obtained rocky desertification level maps for Guizhou Prov-   The state of rocky desertification in Guizhou presents a distribution pattern of "heavy in the west, light in the east, heavy in the south, and light in the north". The areas with severe rocky desertification are mainly distributed in Bijie in the northwest, Liupanshui in the west, and southwest Guizhou and Anshun in the southwest, while the areas with- The state of rocky desertification in Guizhou presents a distribution pattern of "heavy in the west, light in the east, heavy in the south, and light in the north". The areas with severe rocky desertification are mainly distributed in Bijie in the northwest, Liupanshui in the west, and southwest Guizhou and Anshun in the southwest, while the areas without rocky desertification are mainly distributed in southeast Guizhou, Zunyi, and other places.
DEM data were obtained from Geospatial Cloud (http://www.gscloud.cn/, accessed on 13 May 2021). The resolution of the DEM is 30 m. In order to be consistent with the rocky desertification mapping result, we resampled the DEM to 250 m. The suitability factors used included temperature, humidity, light, precipitation, elevation, slope, gross domestic product (GDP), lighting, lithology, and soil, which were obtained from the resource and environmental science data center of the Chinese Academy of Sciences (https://www. resdc.cn/, accessed on 13 May 2021) and the Karst Science Research Data Center (https: //www.resdc.cn/, accessed on 13 May 2021). The suitability factor set is detailed in Figure 3.

Processing Workflow
First, the rocky desertification level (NRD, PRD, LRD, MRD, SRD) map was obtained from the previous research [51], where the multiple line model (ML), random forest model (RF), and support vector machine model (SVM) were compared to the constructed rocky desertification classification model. The overall accuracies (OAs) of these models were 73.8%, 78.2%, and 80.6%, respectively. The SVM model had the best performance. After combining vegetation types and vegetation seasonal phases, the SVM model accuracy reached 91.1%. SVM model was suitable for remote sensing data classification [52]. Second, we analyzed the temporal and spatial evolution of rock desertification from three aspects: dynamic degree, evolution intensity, and evolution velocity. Third, the suitability factor set was designed for future scenario prediction of rocky desertification. We predicted the rocky desertification scenarios in the next 30 years based on an improved CA-Markov model and set three different rocky desertification control scenarios to analyze the future rocky desertification evolution results. The results provide detailed data support for the rocky desertification ecological environment restoration. The processing workflow is shown in Figure 4.

Processing Workflow
First, the rocky desertification level (NRD, PRD, LRD, MRD, SRD) map was obtained from the previous research [51], where the multiple line model (ML), random forest model (RF), and support vector machine model (SVM) were compared to the constructed rocky desertification classification model. The overall accuracies (OAs) of these models were 73.8%, 78.2%, and 80.6%, respectively. The SVM model had the best performance. After combining vegetation types and vegetation seasonal phases, the SVM model accuracy reached 91.1%. SVM model was suitable for remote sensing data classification [52]. Second, we analyzed the temporal and spatial evolution of rock desertification from three aspects: dynamic degree, evolution intensity, and evolution velocity. Third, the suitability factor set was designed for future scenario prediction of rocky desertification. We predicted the rocky desertification scenarios in the next 30 years based on an improved CA-Markov model and set three different rocky desertification control scenarios to analyze the future rocky desertification evolution results. The results provide detailed data support for the rocky desertification ecological environment restoration. The processing workflow is shown in Figure 4.

Transition Matrix
The evolution direction and scale of rocky desertification are typically complicated within a certain period and often vary from place to place. Therefore, the evolution transition matrix was constructed to clarify the general direction and trend of evolution [13,53], shown as Equation (1)

Transition Matrix
The evolution direction and scale of rocky desertification are typically complicated within a certain period and often vary from place to place. Therefore, the evolution transition matrix was constructed to clarify the general direction and trend of evolution [13,53], shown as Equation (1): s 1,1 s 1,2 · · · s 1,k s 2,1 s 2,2 · · · s 2,k · · · · · · · · · s k,1 s k,1 · · · s k,k where S is the area (unit, 10 3 km 2 ); m and n are the start and end time of a period, respectively; and s i,j is the area of rock desertification grade i transformed into rock desertification grade k during this period (unit, 10 3 km 2 ).
The transfer speed is calculated as shown in Equation (2): where p i is the evolution speed of rocky desertification (unit, 10 3 km 2 /a), S 1 is the area at the start of a period, S 2 is the area at the end, and T is the length of the period.

Dynamic Degree
The dynamic degree of each level of rocky desertification evolution can be quantitatively interpreted from the change rate of rocky desertification. This not only allows for comparison of the difference of rocky desertification levels among regions, but also plays an important role in the simulation of future rocky desertification level change. The dynamic degree provides a description of the change rate of rocky desertification. It refers to the quantitative change of a certain level of rocky desertification within a certain period [54,55]. The calculation method is as shown: where K represents the dynamic degree of a certain level of rocky desertification change in the time domain, U a represents the area of a certain level of rocky desertification in the beginning year a, U b represents that in the end year, and T represents the time span.

Traditional CA-Markov Model
The Markov model is a method to predict the probability of an event occurring. It allows for the forecasting of what will happen in the future, using a Markov chain based on the past and known conditions of the event. This method has been widely used in land cover-change modeling. However, when using the traditional Markov model, it is difficult to predict spatial pattern changes in time. The cellular automata (CA) model has strong spatial operation ability and can effectively simulate the spatial changes of the system. The CA-Markov model combines the ability of the CA model to simulate the spatial variation of a complex system with the long-term prediction ability of the Markov model. The model is applied to the future scenario prediction of rocky desertification, which can guarantee the prediction accuracy of rocky desertification level evolution and effectively simulate the spatial change in the rocky desertification pattern. Therefore, this model is scientific and practical [56][57][58].
The CA-Markov model can simulate each grid as a cell in the spatial distribution pattern of rocky desertification, where the rocky desertification level is the cell state. A Markov process is used to match the possible states of rocky desertification, where the area or proportion of phase transformation between rocky desertification types is the state transition probability. Equation (4) can be used to predict the change state of the rocky desertification structure: where S(T) and S(T 0 ) are the rocky desertification structure states at times T and T 0 , respectively, and S m,n is the Markov transition matrix.

Factor Weighting Based on Analytical Hierarchy Process (AHP)
AHP, a useful multicriteria decision-making method, was used to assign weights to each established factor for reasonable assessment [59]. We use the following steps shown in Figure 5 to calculate the weights of the factors based on the AHP method. The Saaty scale is shown in Table 1 [60]. The normalized weights of all factors were examined for the consistency ratio (CR) [61].
Remote Sens. 2022, 14, x FOR PEER REVIEW 9 Figure 5. The processing of factor weighting based on AHP. Since different change accelerations are different, we set up an acceleration to b express and simulate the rocky desertification level matrix. Therefore, we optimize  Moderately The judgment slightly to moderately favors one parameter 4 Intermediate Preference between 3 and 5 5 Strongly The judgment strongly or essentially favors one parameter 6 Intermediate Preference between 5 and 7 7 Very strongly Very strong preference or importance 8 Intermediate Preference between 7 and 9 9 Extremely Quite preferred or quite important

Improved CA-Markov Model
Since different change accelerations are different, we set up an acceleration to better express and simulate the rocky desertification level matrix. Therefore, we optimized the Markov transition matrix based on acceleration. The formula is as follows: where a is the acceleration matrix and the matrix b guarantees that the sum of each row of the Markov transition matrix is 1. The flowchart of the improved CA-Markov prediction model for rocky desertification is shown in Figure 6. where a is the acceleration matrix and the matrix b guarantees that the sum of each row of the Markov transition matrix is 1. The flowchart of the improved CA-Markov prediction model for rocky desertification is shown in Figure 6.

Lee-Sallee Index
In this paper, the modified Lee-Sallee shape index is used to judge the accuracy [43]. The expression is given in Equation (6): where L is the revised Sallee index (with value ranging from 0 to 1), A is the real distri-

Lee-Sallee Index
In this paper, the modified Lee-Sallee shape index is used to judge the accuracy [43]. The expression is given in Equation (6): where L is the revised Sallee index (with value ranging from 0 to 1), A 0 is the real distribution map of rocky desertification in a certain year, and A 1 is the rocky desertification distribution map predicted by the CA-Markov model in this year.

Figure of Merit (FOM) Index
We also perform the comparison quantitatively by computing the components of the FOM. The FOM is a ratio, where the numerator is the intersection of simulated and reference change, while the denominator is the union of simulated and reference change [62]. Equation (7) shows how the FOM is derived from its four components: Misses, Hits, Wrong Hits, and False Alarms [63]. FOM = (Hits)100% Misses + Hits + Wrong Hits + False Alarms (7) where Misses = error owing to observed change predicted as persistence; Hits = correct, observed change predicted as change with the same; Wrong Hits = error owing to observed change predicted as change but with a wrong gaining category; False Alarms = error due to observed persistence predicted as change.

Dynamic Degree of Rocky Desertification Evolution
We calculated the annual dynamic degree of rocky desertification in four periods, as detailed in Table 2. A positive number indicates an increase in the rocky desertification area, while a negative value indicates its decrease. A positive value of the dynamic degree indicates growth in the rate, while a negative value indicates a decrease in the rate. Table 2. Evolution area and annual average dynamic degree of rocky desertification (NRD: no rocky desertification, PRD: potential rocky desertification, LRD: light rocky desertification, MRD: medium rocky desertification, SRD: severe rocky desertification).  The change in NRD area presented a continuous upward trend. The growth of this type of area increased from 5.346 × 10 3  The average annual dynamic degree decreased from 0.288% to −5.337%. Therefore, the NRD area increased year by year, while the area of the other rocky desertification levels continued to decline, indicating that rocky desertification has been well-controlled and governed. The average annual dynamic degree was small during the first period and then gradually increased during the next three periods. This demonstrates that the speed of rocky desertification control has been accelerating annually.

Spatiotemporal Rocky Desertification Evolution Intensity
We set up five rocky desertification grades. The evolution intensity of rocky desertification was the span distance of rocky desertification level transformation, with value ranging from −4 to 4. The spatial and temporal distribution of rocky desertification evolution intensity is shown in Figure 5, where positive values indicate that the rocky desertification situation in the area has been ameliorated and the numerical value indicates the amelioration strength. Meanwhile, a negative value indicates the deteriorating state of rocky desertification, where the numerical value indicates the deterioration strength. Overall, the transformation of rocky desertification mainly occurred between adjacent levels; that is, the evolution intensity of rocky desertification was generally 1 or −1. In particular, there was a large proportion of change between NRD and PRD. The larger the separation level, the smaller the probability of rocky desertification transformation. This is because the development of rocky desertification is a gradual process-in the natural state, it is impossible for the evolution of rocky desertification with great intensity to occur over a short period of time. Of course, there were still some cross-level rocky desertification evolution areas, which need to be paid special attention to in rocky desertification control.
As shown in Figure 7, the amelioration area of rock desertification was larger than that of rock desertification deterioration. This demonstrates that the implementation of slope and ladder, water conservancy project construction, afforestation, and other measures in the rocky desertification control project has been effective and the rocky desertification control effect is good. The amelioration area of rocky desertification was not only from PRD to NRD, but also from SRD to MRD, MRD to LRD, and LRD to PRD.

Spatiotemporal Rocky Desertification Evolution Speed
Dividing the area of rocky desertification evolution by the year span, we obtained the speed of rocky desertification evolution, as shown in Figure 10. Rocky desertification evolution areas rarely appeared in Qiandongnan, which is dominated by nonlithologic areas. In Qiandongnan, the water and temperature conditions are good, the vegetation is flourishing, and rock desertification is rare. The regions with rapid rocky desertification evolution were mainly distributed in areas with medium lithology, but not in Bijie and Liupanshui, where the lithology of rocky desertification is serious. This means that areas with low rocky desertification level are more likely to undergo evolution. In regions with high levels of rocky desertification, it is relatively difficult to control the desertification and its evolution is slower. From 2001 to 2005, the rapid evolution regions were mainly distributed in Guiyang and northern Zunyi. From 2005 to 2010, the rapid evolution regions were mainly distributed in Guiyang and Anshun. From 2010 to 2015, Zunyi and Tongren were the main regions where rocky desertification evolved rapidly. From 2015 to 2020, the rapid evolution regions of rocky desertification were widespread, radiating from the capital city Guiyang to the surrounding areas, mainly distributed along the Wujiang and Liujiang rivers. This means that the change of rocky desertification was affected by both natural conditions and human disturbances.
As shown in Figure 11, the evolution speeds of PRD, LRD, MRD, and SRD were all positive. This means that the four levels of rocky desertification all showed a trend of alleviation. The evolution speeds of NRD and PRD were relatively high. Considering Table 1, which shows the area and proportion of rocky desertification at each level, due to the obvious control effect of LRD, MRD, and SRD, their proportions decreased annually, with most of them eventually evolving to PRD. The evolution speeds of NRD and PRD were higher than those of other rocky desertification levels. The evolution speeds of NRD were 2.591 × 10 3 km 2 /a, 3.155 × 10 3 km 2 /a, 2.263 × 10 3 km 2 /a, and 1.544 × 10 3 km 2 /a in the four periods, respectively. The evolution speeds of PRD were 1.551 × 10 3 km 2 /a, 1.847 × 10 3 km 2 /a, 2.193 × 10 3 km 2 /a, and 2.884 × 10 3 km 2 /a in the four periods, respectively. This indicates that the evolution between NRD and PRD was frequent, widely distributed, and that they can easily alternate. The tendency of rock desertification can occur if little attention is paid to vegetation and soil protection in NRD areas. If the PRD regions are governed strongly, they can easily be transformed into NRD. Therefore, not only should the control of MRD and SRD be strengthened, but attention should also be paid to prevention and control in PRD regions. The slow evolution of rocky desertification occurred mostly in the regions where the evolution intensity was larger than two levels. This means that the evolution of rocky desertification has hierarchical inertia, and the higher the level of rocky desertification, the slower its evolution and the more difficult its governance.

Spatiotemporal Rocky Desertification Evolution Speed
Dividing the area of rocky desertification evolution by the year span, we obtained the speed of rocky desertification evolution, as shown in Figure 10. Rocky desertification evolution areas rarely appeared in Qiandongnan, which is dominated by nonlithologic distributed in Guiyang and northern Zunyi. From 2005 to 2010, the rapid evolution regions were mainly distributed in Guiyang and Anshun. From 2010 to 2015, Zunyi and Tongren were the main regions where rocky desertification evolved rapidly. From 2015 to 2020, the rapid evolution regions of rocky desertification were widespread, radiating from the capital city Guiyang to the surrounding areas, mainly distributed along the Wujiang and Liujiang rivers. This means that the change of rocky desertification was affected by both natural conditions and human disturbances. As shown in Figure 11, the evolution speeds of PRD, LRD, MRD, and SRD were all positive. This means that the four levels of rocky desertification all showed a trend of alleviation. The evolution speeds of NRD and PRD were relatively high. Considering Table 1, which shows the area and proportion of rocky desertification at each level, due to the obvious control effect of LRD, MRD, and SRD, their proportions decreased annually, with most of them eventually evolving to PRD. The evolution speeds of NRD and PRD were higher than those of other rocky desertification levels. The evolution speeds of NRD were 2.591 × 10 3 km 2 /a, 3.155 × 10 3 km 2 /a, 2.263 × 10 3 km 2 /a, and 1.544 × 10 3 km 2 /a in the 3 2 3 The net profit speed of rocky desertification evolution refers to the sum of the speeds of the deteriorating and ameliorating parts of a certain level of rocky desertification. The annual statistics are shown in Table 3, which indicates that the average annual evolving net profit speed was positive. The evolution speed increased annually, from 1.704 × 10 3 km 2 /a in 2001-2005 to 4.937 × 10 3 km 2 /a in 2015-2020, demonstrating that the effect of rocky desertification control is becoming more and more obvious.

Improved CA-Markov Prediction Model Performance
Markov models predict the future based on the probability of events occurring in the past. From the above analysis in this paper, the evolution of rocky desertification in Guizhou showed a trend of accelerating amelioration over the past two decades, especially after 2010. Using the traditional CA-Markov model to predict rocky desertification in Guizhou cannot reflect this accelerating trend.
For our experiment, we used the remote sensing inversion results of rocky desertification from 2001 to 2015 to construct a Markov transition matrix with an interval of 5 years. The set of suitability factors was composed of 10 factors, namely temperature, humidity, sun, precipitation, elevation, slope, GDP, light, lithology, and soil, as shown in Figure 4. We used the AHP method to calculate the factor weight, and the result is shown in Table 4. The CR value is 0.04, which is lower than 0.1, meeting the threshold value of the AHP method. In this paper, 5 rocky desertification levels (as the dependent variables) and the corresponding 10 suitability factors (as independent variables) of the image were calculated by logistic regression. As shown in Table 5, the ROC regression results for the selected 10 suitability factors reflected the spatial distribution of rocky desertification well. The reflection ability for NRD was the strongest, with ROC of 0.9209. The ROC values for MRD, LRD, SRD, and PRD were 0.8637, 0.8481, 0.8361, and 0.8169, respectively. The ROC values of every rocky desertification level were all above 0.8, which indicates that there was good consistency between the predicted and real distributions. The net profit speed of rocky desertification evolution refers to the sum of the speeds of the deteriorating and ameliorating parts of a certain level of rocky desertification. The annual statistics are shown in Table 3, which indicates that the average annual evolving net profit speed was positive. The evolution speed increased annually, from 1.704 × 10 3 km 2 /a in 2001-2005 to 4.937 × 10 3 km 2 /a in 2015-2020, demonstrating that the effect of rocky desertification control is becoming more and more obvious. Markov models predict the future based on the probability of events occurring in the past. From the above analysis in this paper, the evolution of rocky desertification in Guizhou showed a trend of accelerating amelioration over the past two decades, especially after 2010. Using the traditional CA-Markov model to predict rocky desertification in Guizhou cannot reflect this accelerating trend.
For our experiment, we used the remote sensing inversion results of rocky desertification from 2001 to 2015 to construct a Markov transition matrix with an interval of 5 years. The set of suitability factors was composed of 10 factors, namely temperature, humidity, sun, precipitation, elevation, slope, GDP, light, lithology, and soil, as shown in Figure 4. We used the AHP method to calculate the factor weight, and the result is shown in Table 4. The CR value is 0.04, which is lower than 0.1, meeting the threshold value of the AHP method. In this paper, 5 rocky desertification levels (as the dependent variables) and the corresponding 10 suitability factors (as independent variables) of the image were calculated by logistic regression. As shown in Table 5, the ROC regression results for the selected 10 suitability factors reflected the spatial distribution of rocky desertification well. The reflection ability for NRD was the strongest, with ROC of 0.9209. The ROC values for MRD, LRD, SRD, and PRD were 0.8637, 0.8481, 0.8361, and 0.8169, respectively. The ROC values of every rocky desertification level were all above 0.8, which indicates that there was good consistency between the predicted and real distributions. After constructing the suitability factor set, the CA-Markov model was used to predict the rocky desertification in 2020, as shown in Figure 12b. Then, based on the remote sensing inversion result of rocky desertification in 2020, shown in Figure 12a, the Markov transition matrix from 2001 to 2015 was optimized, as shown in Table A1. Using an iterative operation, where the difference between the results of the two iterations is less than 0.1, the acceleration parameter a was calculated. The optimized rocky desertification result for 2020 is shown in Figure 12c. Finally, a was used for future scenario prediction of rocky desertification after 2020.  After constructing the suitability factor set, the CA-Markov model was used to predict the rocky desertification in 2020, as shown in Figure 12b. Then, based on the remote sensing inversion result of rocky desertification in 2020, shown in Figure 12a, the Markov transition matrix from 2001 to 2015 was optimized, as shown in Table A1. Using an iterative operation, where the difference between the results of the two iterations is less than 0.1, the acceleration parameter a was calculated. The optimized rocky desertification result for 2020 is shown in Figure 12c. Finally, a was used for future scenario prediction of rocky desertification after 2020.  Combined with Figure 12 and Table 6, the improved CA-Markov prediction model of rocky desertification was used to calculate the results in 2020. Compared with the traditional CA-Markov model, the method proposed in this paper improved the consistency of the proportion of each rocky desertification level. The RMSE of each rocky desertification level decreased from 2.016 to 0.056. The proportions of the confusion matrix between simulation result and reality result are shown in Table 7. As shown in Table 8, the improved CA-Markov model was used to simulate the rocky desertification data in 2020, and the Lee-Sallee shape value was 0.723, higher than that obtained with the traditional CA-Markov model (which was 0.681). Compared with the Lee-Salle value of 0.629 obtained by Ma et al. [37], the improved CA-Markov model had higher accuracy. The improved CA-Markov model also had a higher overall accuracy and Kstandard value. Thus, the improved CA-Markov model can accurately reflect the evolution trajectory of rocky desertification.  [37] 0.629 - The FOM result is shown in Figure 13.

Future Scenario Prediction Based on Improved CA-Markov Prediction Model
Based on the improved CA-Markov model, taking the remote sensing inversion rocky desertification result in 2020 as a baseline, and using the set of suitability factors in 2020, we set the iteration interval to 5 years and predicted the future scenario of rocky desertification from 2025 to 2050. The results are shown in Figure 14.

Future Scenario Prediction Based on Improved CA-Markov Prediction Model
Based on the improved CA-Markov model, taking the remote sensing inversion rocky desertification result in 2020 as a baseline, and using the set of suitability factors in 2020, we set the iteration interval to 5 years and predicted the future scenario of rocky desertification from 2025 to 2050. The results are shown in Figure 14.

Future Scenario Prediction Based on Improved CA-Markov Prediction Model
Based on the improved CA-Markov model, taking the remote sensing inversion rocky desertification result in 2020 as a baseline, and using the set of suitability factors in 2020, we set the iteration interval to 5 years and predicted the future scenario of rocky desertification from 2025 to 2050. The results are shown in Figure 14. In this paper, the amelioration part of the future rocky desertification scenario is defined as the region where the evolution intensity of rocky desertification is greater than or equal to 2, while the deterioration part of rocky desertification is the region is where the evolution intensity of rocky desertification is less than or equal to −2. The amelioration and deterioration parts in the future scenarios of rocky desertification from 2020 to 2030 were accurately located, as shown in Figure 15. From 2020 to 2030, the amelioration regions of rocky desertification were significantly larger than the deterioration regions. The amelioration regions in future scenarios of rocky desertification were widely distributed, in areas such as Bijie, Liupanshui, Qianxinan, Guiyang, and Qiannan. The deterioration regions in future scenarios were mainly distributed in Bijie and Liupanshui. Bijie and Liupanshui are the two regions featuring both amelioration and deterioration areas at the same time, which should be taken into consideration when planning rocky desertification governance.   In this paper, the amelioration part of the future rocky desertification scenario is defined as the region where the evolution intensity of rocky desertification is greater than or equal to 2, while the deterioration part of rocky desertification is the region is where the evolution intensity of rocky desertification is less than or equal to −2. The amelioration and deterioration parts in the future scenarios of rocky desertification from 2020 to 2030 were accurately located, as shown in Figure 15. From 2020 to 2030, the amelioration regions of rocky desertification were significantly larger than the deterioration regions. The amelioration regions in future scenarios of rocky desertification were widely distributed, in areas such as Bijie, Liupanshui, Qianxinan, Guiyang, and Qiannan. The deterioration regions in future scenarios were mainly distributed in Bijie and Liupanshui. Bijie and Liupanshui are the two regions featuring both amelioration and deterioration areas at the same time, which should be taken into consideration when planning rocky desertification governance. In this paper, the amelioration part of the future rocky desertification scenario is defined as the region where the evolution intensity of rocky desertification is greater than or equal to 2, while the deterioration part of rocky desertification is the region is where the evolution intensity of rocky desertification is less than or equal to −2. The amelioration and deterioration parts in the future scenarios of rocky desertification from 2020 to 2030 were accurately located, as shown in Figure 15. From 2020 to 2030, the amelioration regions of rocky desertification were significantly larger than the deterioration regions. The amelioration regions in future scenarios of rocky desertification were widely distributed, in areas such as Bijie, Liupanshui, Qianxinan, Guiyang, and Qiannan. The deterioration regions in future scenarios were mainly distributed in Bijie and Liupanshui. Bijie and Liupanshui are the two regions featuring both amelioration and deterioration areas at the same time, which should be taken into consideration when planning rocky desertification governance.

Future Scenario Prediction of Rocky Desertification under Different Governance Scenarios
Three different control scenarios were considered: historical evolution, major governance, and complete governance. In the major governance scenario, three levels of rocky desertification-PRD, LRD, and MRD-were taken as the key control aims of the restoration of karst rocky desertification. These three levels of karst rocky desertification are more easily transformed into other levels in the middle stage of the succession process of karst rocky desertification. Therefore, they can be more effectively governed and restored [38]. Previous studies have shown that, due to the favorable rainfall and heat conditions for plant growth, through the implementation of vegetation restoration and management projects, the vegetation coverage can be increased by 10-30% within 10 years [37,64]. Considering the existing measures for controlling and restoring rocky desertification in the historical evolution situation, in the major governance scenario, the mean value was taken as 20%. In this study, the transition matrix of karst rocky desertification under historical evolution scenarios was adjusted, in order to calculate the transition matrix for the major governance scenario. Here, the probability of transformation of PRD, LRD, and MRD to more serious karst rocky desertification level was decreased by 20% and that to minor karst rocky desertification was increased by 20%.
The complete governance scenario emphasizes ecological restoration and management strategies that restrict regional development, with an aim to realize the comprehensive restoration of karst rocky desertification. Human activities in karst areas are restricted to relieve the pressure on the land. Based on the rock desertification transfer matrix of complete governance scenarios, the adjustment of transformation probability of karst rock desertification was extended to two grades: NRD and SRD. The rules of the transformation matrix modification were as follows: the probability of all levels of karst rocky desertification being converted to a more serious level was reduced by 20%, while that to minor karst rocky desertification was increased by 20%.
The prediction results of rocky desertification under different control scenarios in 2030 are shown in Figure 16 and Table 9.

Future Scenario Prediction of Rocky Desertification under Different Governance Scenarios
Three different control scenarios were considered: historical evolution, major governance, and complete governance. In the major governance scenario, three levels of rocky desertification-PRD, LRD, and MRD-were taken as the key control aims of the restoration of karst rocky desertification. These three levels of karst rocky desertification are more easily transformed into other levels in the middle stage of the succession process of karst rocky desertification. Therefore, they can be more effectively governed and restored [38]. Previous studies have shown that, due to the favorable rainfall and heat conditions for plant growth, through the implementation of vegetation restoration and management projects, the vegetation coverage can be increased by 10-30% within 10 years [37,64]. Considering the existing measures for controlling and restoring rocky desertification in the historical evolution situation, in the major governance scenario, the mean value was taken as 20%. In this study, the transition matrix of karst rocky desertification under historical evolution scenarios was adjusted, in order to calculate the transition matrix for the major governance scenario. Here, the probability of transformation of PRD, LRD, and MRD to more serious karst rocky desertification level was decreased by 20% and that to minor karst rocky desertification was increased by 20%.
The complete governance scenario emphasizes ecological restoration and management strategies that restrict regional development, with an aim to realize the comprehensive restoration of karst rocky desertification. Human activities in karst areas are restricted to relieve the pressure on the land. Based on the rock desertification transfer matrix of complete governance scenarios, the adjustment of transformation probability of karst rock desertification was extended to two grades: NRD and SRD. The rules of the transformation matrix modification were as follows: the probability of all levels of karst rocky desertification being converted to a more serious level was reduced by 20%, while that to minor karst rocky desertification was increased by 20%.
The prediction results of rocky desertification under different control scenarios in 2030 are shown in Figure 16 and Table 9.

Improvement of Future Prediction Accuracy by Modifying Markov Transition Matrix
On one hand, the proposed prediction model of rocky desertification based on the CA-Markov method depends on the reliability and stability of the transition probability matrix. On the other hand, it depends on the interpretation of a suitability factor set for rocky desertification. In this paper, the transfer matrix from 2001 to 2015 was taken as the probability matrix and the rocky desertification in 2015 served as a base image. The rocky desertification in 2020 was predicted by the proposed CA-Markov model and compared with the real results in 2020. Then, the acceleration coefficient matrix was used to adjust the Markov transition probability matrix, thus improving the prediction accuracy. This is similar to previous studies, such as that of Xu et al. [64] who quantified the impact of rocky desertification expansion on a designated area. They analyzed the neighborhood effect of karst rocky desertification of different grades and corrected the probability matrix of rocky desertification. This strategy is suitable for medium-or high-resolution image analysis. However, this study was based on MODIS 250 m resolution data products, which does not quite fit the scale of neighborhood effects. Therefore, we used the error term between the predicted results and the real results to iteratively optimize the probability matrix, thus improving the matrix optimization and prediction accuracy.
Thanks to the support of sufficient sample data, the method of [51], in this paper, obtained remote sensing inversion results of rocky desertification level in 2020 bearing similarity to the real results, providing a reliable data basis for CA-Markov prediction of future prospects. Although there were errors within a certain allowable range, the results can still be considered reliable. The suitability factor set in this paper was derived from the average values from 2001 to 2020, including such factors as temperature, precipitation, and light. The suitability factor set can reflect the level difference of rocky desertification. However, there were still time differences between these suitability factors in different

Improvement of Future Prediction Accuracy by Modifying Markov Transition Matrix
On one hand, the proposed prediction model of rocky desertification based on the CA-Markov method depends on the reliability and stability of the transition probability matrix. On the other hand, it depends on the interpretation of a suitability factor set for rocky desertification. In this paper, the transfer matrix from 2001 to 2015 was taken as the probability matrix and the rocky desertification in 2015 served as a base image. The rocky desertification in 2020 was predicted by the proposed CA-Markov model and compared with the real results in 2020. Then, the acceleration coefficient matrix was used to adjust the Markov transition probability matrix, thus improving the prediction accuracy. This is similar to previous studies, such as that of Xu et al. [64] who quantified the impact of rocky desertification expansion on a designated area. They analyzed the neighborhood effect of karst rocky desertification of different grades and corrected the probability matrix of rocky desertification. This strategy is suitable for medium-or high-resolution image analysis. However, this study was based on MODIS 250 m resolution data products, which does not quite fit the scale of neighborhood effects. Therefore, we used the error term between the predicted results and the real results to iteratively optimize the probability matrix, thus improving the matrix optimization and prediction accuracy.
Thanks to the support of sufficient sample data, the method of [51], in this paper, obtained remote sensing inversion results of rocky desertification level in 2020 bearing similarity to the real results, providing a reliable data basis for CA-Markov prediction of future prospects. Although there were errors within a certain allowable range, the results can still be considered reliable. The suitability factor set in this paper was derived from the average values from 2001 to 2020, including such factors as temperature, precipitation, and light. The suitability factor set can reflect the level difference of rocky desertification. However, there were still time differences between these suitability factors in different years. In future research, we will consider varying factors in different years to compare and optimize the suitability factor set, in order to improve the performance of the prediction model.

Temporal and Spatial Evolution of Rocky Desertification and Its Influence
From the spatial perspective, the distribution pattern of rocky desertification in Guizhou was "heavy in the west and south, light in the east and north", consistent with the results of previous studies [65,66]. The spatiotemporal change of rocky desertification has a great relationship with the natural environment and geological background. In this paper, we found that the severe rocky desertification areas were mainly distributed in three regions: Bijie in northwest Guizhou, Liypanshui in west Guizhou, and Qinxinan and Anshun in southwest Guizhou. Pure carbonate is widely distributed in these regions, and the terrain has a larger slope and high soil erosion. However, non-rocky desertification was mainly distributed in Qindongnan and Zunyi in Guizhou, which are nonkarst areas with good geothermal and superior forest site conditions.
Considering the time dimension, the evolution of rocky desertification in Guizhou was obvious from 2001 to 2020, predominantly showing a gradual amelioration trend. The evolution of rocky desertification was continuous, mostly occurring in adjacent levels: as expected, the greater the separation level, the smaller the evolution probability of rocky desertification. The net profit speed of each rocky desertification level transfer showed a trend of improvement. Among them, the transfer speed of potential rocky desertification was the fastest. Especially from 2015 to 2020, the transfer net profit speed reached 4.937 × 10 3 km 2 /a. The dynamic degrees from 2015 to 2020 were also the most dramatic. The area of PRD was relatively large and easy to transfer, and the transfer direction was alternating. Therefore, it would be ideal to pay attention to the prevention and control of PRD deterioration. Meanwhile, SRD was the slowest to evolve, with net profit speed of transfer being only 0.423 × 10 3 km 2 /a from 2015 to 2020. The change rate of rocky desertification is inversely proportional to its level, which is consistent with previous research [67]. The spatial and temporal evolution pattern of rocky desertification is not only related to the natural environment, including geology and climate, but also closely related to human activities and the background of China's social and economic development and ecological construction during this period. Since 1999, China has implemented the "six major" forestry projects. In 2003, the decision on "Accelerating the Development of Forestry" was issued, and the forest resources in Guizhou were further improved. After 2010, with the further implementation of the greening and beautification project and the construction of beautiful countryside, forest site conditions have been greatly improved and forest resources are growing. Guizhou has been practicing the scientific judgment that "lucid waters and lush mountains are invaluable assets". Thus, it can be stated that rocky desertification has been ameliorated based on scientific management and the returning farmland to forests project.

Discussion on Suitability Governance of Future Scenarios of Rocky Desertification
Based on the prediction of the improved CA-Markov model, the rocky desertification in Guizhou will continue to be ameliorated from 2025 to 2050. In particular, for Qianxinan in southwest Guizhou, rocky desertification will be strongly ameliorated. With the development of economy and society, the national investment in ecological environment management projects is also increasing. The implementation of such projects will promote the amelioration of the rocky desertification situation. The ecological environment is becoming better and better. The regions with severe rocky desertification in 2025 were mainly distributed in Bijie, Liupanshui, Qindongnan, and Anshun, and the areas with severe rocky desertification will be smaller from 2025 to 2030 than those from 2020 to 2025. However, Bijie and Liupanshui are still relatively serious regions. Therefore, these regions should be given priority for the control of rocky desertification.
It was found that the occurrence of rocky desertification is lower in the southeast and south of Guizhou; furthermore, these areas are dominated by forest. This demonstrates that vegetation can prevent rocky desertification. Forestation is an effective method to prevent rocky desertification. The main vegetation types in rocky desertification regions and non-rocky desertification regions can be further studied in future research. Studying what kind of forest species and tree species are suitable for planting in the rocky desertification areas can provide effective data and planning support for future rocky desertification governance.

Conclusions
In this study, we analyzed the temporal and spatial evolution law of rocky desertification in Guizhou, China, from 2001 to 2020, considering three aspects: dynamic degree, intensity, and speed of rocky desertification evolution. The CA-Markov model was improved to better predict the rocky desertification level. The Lee-Sallee index reached 0.723, overall accuracy reached 0.839, and FOM reached 0.530, indicating a significant improvement over the traditional CA-Markov model. The future scenario of rocky desertification in Guizhou, China, from 2025 to 2050 was predicted, and three different governance scenarios were set to analyze the possible improvement and deterioration areas of rocky desertification in the future. The conclusions of this paper are as follows: (1) Considering the dynamic degree of rocky desertification, the annual average dynamic degree gradually increased, showing that the speed of rocky desertification control has accelerated year by year. The government has vigorously carried out a pilot project of comprehensive rocky desertification control, which has curbed the deterioration from rocky desertification in Guizhou Province. From the perspective of the intensity and speed of rock desertification evolution, the transformation rates of NRD and PRD were the highest, while that of SRD was the lowest. Therefore, in the process of rocky desertification control, relevant actors should pay attention to the effect of MRD and SRD control, as well as avoiding the deterioration of NRD and PRD. (2) By improving the CA-Markov model, we improved the accuracy of prediction for rocky desertification future scenarios and could better explore the evolution rule of rocky desertification. Based on the improved CA-Markov model results, we found that the serious regions of rocky desertification in the future are mainly in Bijie, Liupanshui, Qindongnan, and Anshun. Thus, more attention should be paid to these areas in the control of rocky desertification. The improved Markov model improves the prediction accuracy, can predict the future rocky desertification scenario more accurately, and is conducive to the accurate analysis of the future rocky desertification landscape pattern. At the same time, through the simulation of three different control levels, it provides a basis for scientific control of rocky desertification. The accurate positioning of areas that may deteriorate in the next 10 years can well guide the accurate control of local rocky desertification in Guizhou. In future research, more landscape prediction models should be tested in the future scenario prediction of rocky desertification, and better prediction methods should be explored.

Conflicts of Interest:
The authors declare no conflict of interest. The founding sponsors had no role in the design of the study; in the collection, analyses or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.