Extraction of Rocky Desertiﬁcation Information in the Karst Area Based on the Red-NIR-SWIR Spectral Feature Space

: The complex topography, severe surface fragmentation and landscape heterogeneity of the karst region of southwest China make it extremely difﬁcult to extract information on rocky desertiﬁcation in the region. In order to overcome the disadvantages of the surface parameter-based feature space approach, which is difﬁcult to construct and apply, this study uses the reﬂectance of Landsat 8 Operational Land Imager (OLI) in the red (Red), near-infrared (NIR) and shortwave infrared (SWIR) bands as the feature variables, and establishes a two-dimensional SWIR-NIR, Red-NIR and SWIR-Red reﬂectance spectral feature space. The three models of perpendicular rocky desertiﬁcation index 1 (PRDI1), perpendicular rocky desertiﬁcation index 2 (PRDI2) and perpendicular rocky desertiﬁcation index 3 (PRDI3) were also constructed based on the variation of the degree of rocky desertiﬁcation in each spectral feature space. The accuracy of the rocky desertiﬁcation extracted by these three index models was veriﬁed and compared with the karst rocky desertiﬁcation index (KRDI) and rocky desertiﬁcation difference index (RSDDI), which are constructed based on the surface parameter feature space. The results show that: (1) The waveband reﬂectance-based feature space model provides a new method for large-scale rocky desertiﬁcation information extraction, characterized by easy data acquisition, simple index calculation and good stability, and is conducive to the monitoring and quantitative analysis of rocky desertiﬁcation in karst areas. (2) The overall accuracy and Kappa coefﬁcient of PRDI1 are 0.829 and 0.784, respectively, both higher than other index models, showing the best applicability, accuracy and effectiveness in rocky desertiﬁcation information extraction. (3) According to the results extracted from PRDI1, the total area of rocky desertiﬁcation in Huaxi District of Guizhou province is 320.44 km 2 , with the more serious grades of rocky desertiﬁcation, such as severe and moderate, mainly distributed in the southwestern, western and southeastern areas of Huaxi District. This study provides important information on the total area and spatial distribution of different degrees of rocky desertiﬁcation in the study area, and these results can be used to support the local government’s ecological and environmental management decisions. The method proposed in this study is a scientiﬁc and necessary complement to the characteristic spatial methods based on different surface parameters, and can provide important methodological support for the rapid and efﬁcient monitoring of karstic rocky desertiﬁcation over large areas.


Introduction
The karst region of southwest China is one of the three largest karst regions in the world, covering an area of about 540,000 km 2 , and is the largest continuous exposed carbonate rocky distribution area in the world [1].Karst rocky desertification is an extreme manifestation of land degradation in a fragile subtropical karst environment, where both natural factors and human activities have led to the gradual degradation of vegetation, massive soil loss and the extensive exposure of bedrock [2,3].Its serious impact on the local ecological environment and hydrological conditions has led to frequent geological disasters such as droughts, floods, landslides and ground subsidence [4].Karst rocky desertification also affects local agricultural and forestry production, seriously hindering sustainable local socio-economic development and making local residents live in poverty for a long time [5].Due to the complex topography, severe surface fragmentation and landscape heterogeneity of the karst region of southwest China, it is extremely difficult to extract information on the spatial distribution, area and extent of rocky desertification in the region.
A great deal of research has been carried out to address the difficulty of extracting information on rocky desertification.Traditional research methods usually rely on field surveys of exposed bedrock, vegetation, soil thickness and related indicators, but field surveys are often costly, time-consuming and labor-intensive, especially for large-scale monitoring [6,7].Remote sensing technology has gradually become the main means of rocky desertification monitoring because of its ability to obtain large-scale rocky desertification information in a rapid and timely manner and to achieve large-scale rocky desertification monitoring and assessment [8][9][10][11].Conventional remote sensing image classification methods such as supervised classification and unsupervised classification are highly subjective and have difficulty distinguishing the grade of rocky desertification, cannot directly indicate the development status of rocky desertification, and are not suitable for extracting spatial distribution, area and grade information of rocky desertification [12,13].The mixed image element decomposition method (SMA) is a relatively common method in the study of karst desertification, but due to the discontinuous distribution of topography in karst areas and the influence of human activities, weathering and erosion, end element variation is common and it is very difficult to obtain feature end elements, so the SMA method cannot easily obtain information on rocky desertification [14,15].Traditional vegetation indices are based on the specific spectral characteristics of green vegetation, reflecting the "greenness" of vegetation information, which can reflect the growth of vegetation, but cannot directly indicate the development of rocky desertification and the degree of development [16].In recent years, the characteristic spatial model method based on surface parameters has been gradually applied to the study of karstic desertification; this can distinguish different grades of desertification and directly indicate the degree of development of desertification, and some good research results have been achieved.
A spectral feature space is a spectral space consisting of multi-band spectral information or two or more of the resulting surface ecophysical parameters [17].In addition to the fruitful research results in the fields of desertification and salinization [18,19], the characteristic spatial models constructed based on different surface parameters have also achieved some good research results in the study of rocky desertification.Zhang, Jie et al. [20] constructed a new karst rocky desertification index (KRDI) based on three surface parameters, greenness, humidity and brightness, and used the temporal spectral feature space (TSFS) model to determine the rank of rocky desertification.The results proved that the KRDI is a sensitive indicator of rocky desertification and can be used for the direct extraction of rocky desertification information.Luo Jie et al. [21] constructed the rocky desertification difference index (RSDDI) based on NDVI and Albedo surface parameter feature space, and extracted information on the spatial distribution and extent of rocky desertification in Dougu Township, and the results showed that the RSDDI could extract and classify rocky desertification information more accurately.Guo, Bing et al. [12] constructed two types of feature space models (point-to-point and point-to-line) based on seven surface parameters to monitor the rocky desertification in Dafang County, and the results proved that the feature space models constructed based on different surface parameters could extract rocky desertification in karst areas, providing a new research method for the effective large-scale monitoring of rocky desertification.
The above model of the rocky desertification index, constructed on the basis of the characteristic space of different surface parameters, has the advantages of reflecting changes in the hydrothermal combination of the surface cover and having a clear biophysical significance.However, this method also has some disadvantages, mainly the complicated process of calculating surface parameters and the inclusion of certain errors, which gradually accumulate and increase the uncertainty of the index model calculation [22].In addition, the values of the surface parameters tend to change with the study area, time and season, etc., resulting in changes in the shape of the feature space constructed using them, and it is difficult to form a feature space with a specific shape, such as a triangle or trapezoid, making its application more difficult.In order to solve this problem, this study considers the construction of a rocky desertification index model by establishing a band reflectance spectral feature space, which has the advantages of a simple construction process, easy access to model parameters and good stability.It has been studied extensively in the fields of soil moisture and land drought monitoring [23][24][25], and has achieved good research results, but it has not been reported in the field of karstic rocky desertification research so far, and its applicability to the extraction of spatial distribution and the degree information of karstic rocky desertification deserves further study.
Among the many remote sensing images, the multispectral Landsat series of remote sensing images has been one of the most widely used remote sensing images in recent decades due to its high spatial and temporal resolution, wide coverage, and free availability [26,27].Many studies have used the red (Red: 630-690 nm) band, near infrared (NIR: 775-900 nm) band and shortwave (SWIR: 1550-1750 nm) band of Landsat series remote sensing images to construct different karst rocky indices and vegetation indices, such as normalized rocky indices, specific vegetation indices and normalized vegetation indices [28][29][30], and superimposed them for analysis to invert rocky desertification information.In addition, many studies have shown that the reflectance of surface cover from, for example, rocks, vegetation and buildings varies the most in these three bands, and indices constructed using these three bands can better distinguish between them [14,28].This shows that karst rocky desertification information has a strong correlation with the Red, NIR and SWIR bands of the multispectral imagery.Therefore, this study considers the use of these three bands to construct a feature space model to extract information on the spatial distribution, area and extent of karstic rocky desertification.
In summary, in order to overcome the disadvantages of the feature space method based on different surface parameters such as a complex calculation process, difficult access to model parameters and poor stability, and to improve the applicability of the feature space method, as well as to save extraction costs and achieve the rapid and efficient extraction of karstic rocky desertification information, a new method based on band reflectance spectral feature space was developed in this study.The method uses the Red, NIR and SWIR bands of the multispectral Landsat-8 OLI as feature variables to construct a spectral feature space, and by studying the spatial distribution of different land covers and the variation of different degrees of desertification in the reflectance spectral feature space, an index model can be established that can quickly, accurately and directly indicate the degree of desertification development.The model can be easily replicated on other multispectral images and applied to other karst areas.

Study Area
As shown in Figure 1a, this study area was located in Huaxi District, Guiyang City, Guizhou Province, between 106 • 27 ~106 • 52 E and 26 • 11 ~26 • 34 N, with a total area of 964.32 km 2 , in which karst landscapes are widely distributed, accounting for 94% of the total area [31].The topography of the region is shown in Figure 1b, with a high north and south and low middle, with the highest features in the southeast, ranging from 1014 m to 1684 m above sea level, and the central part being around 1100 m above sea level.The landform type of the whole region is dominated by mountains and hills, with denuded hills interspersed with basins, valleys and depressions, and the landscape is fragmented [32].The geology of the area is complex, with a rich variety of rocky types, including lime rocky, dolomite, common sand shale and purple sand shale [33].The climate is subtropical, humid and monsoonal, with no cold winters or hot summers, an average annual temperature of 14.9 • C and abundant rainfall of 1178.3 mm. Figure 1c shows the distribution of sample points in the study area.
landform type of the whole region is dominated by mountains and hills, with denuded hills interspersed with basins, valleys and depressions, and the landscape is fragmented [32].The geology of the area is complex, with a rich variety of rocky types, including lime rocky, dolomite, common sand shale and purple sand shale [33].The climate is subtropical, humid and monsoonal, with no cold winters or hot summers, an average annual temperature of 14.9 °C and abundant rainfall of 1178.3 mm. Figure 1c shows the distribution of sample points in the study area.

Data Collection and Processing
The remote sensing image used in this study was a Landsat-8 OLI image downloaded from the Geospatial Data Cloud (http://www.Gscloud.cn,accessed on 21 June 2022).The imaging time was 18 March 2020, with a spatial resolution of 30 m and a temporal resolution of 16 days, possessing a high spatial and temporal resolution.The pre-processing process consisted mainly of radiation correction and atmospheric correction, the purpose of which was to remove various aberrations from the image data [34], and the purpose of atmospheric correction is the process of removing radiometric errors caused by atmospheric effects and inverting the true surface reflectance of the features [35].The processing process was as follows: first, radiometric calibration using the Radiation toolbox in Environment for Visualizing Images 5.3 (ENVI5.3)was carried out to calculate top-of-atmosphere reflectance, converting the image grey-scale values into physically meaningful irradiance values [36].Secondly, the Fast Line-of-sight Atmospheric Analysis of Hypercubes (FLAASH) atmospheric correction model in ENVI5.3 was used to invert the surface true reflectance, and after atmospheric correction, the apparent reflectance to the earth was obtained for each band [37].Finally, the corrected remote sensing image was cropped, using the administrative division vector map of Huaxi District to obtain a remote sensing image map covering the whole Huaxi District.

Data Collection and Processing
The remote sensing image used in this study was a Landsat-8 OLI image downloaded from the Geospatial Data Cloud (http://www.Gscloud.cn,accessed on 21 June 2022).The imaging time was 18 March 2020, with a spatial resolution of 30 m and a temporal resolution of 16 days, possessing a high spatial and temporal resolution.The pre-processing process consisted mainly of radiation correction and atmospheric correction, the purpose of which was to remove various aberrations from the image data [34], and the purpose of atmospheric correction is the process of removing radiometric errors caused by atmospheric effects and inverting the true surface reflectance of the features [35].The processing process was as follows: first, radiometric calibration using the Radiation toolbox in Environment for Visualizing Images 5.3 (ENVI5.3)was carried out to calculate top-of-atmosphere reflectance, converting the image grey-scale values into physically meaningful irradiance values [36].Secondly, the Fast Line-of-sight Atmospheric Analysis of Hypercubes (FLAASH) atmospheric correction model in ENVI5.3 was used to invert the surface true reflectance, and after atmospheric correction, the apparent reflectance to the earth was obtained for each band [37].Finally, the corrected remote sensing image was cropped, using the administrative division vector map of Huaxi District to obtain a remote sensing image map covering the whole Huaxi District.
Since karstic desertification occurs only in karst areas, and there were some non-karst areas in the study area, this study used 1:50,000 lithological vector data (http://www.karstdata.cn/,accessed on 22 June 2022) from Huaxi District to mask non-karst areas before constructing the rocky desertification index model [12].As the reflectance of rocky areas is similar to that of built-up areas and bare land, it is difficult to distinguish them using the spectral features of the images, which can easily lead to misclassification of rocky desertification [38], so built-up areas and bare land need to be removed before the construction of the rocky desertification supervisory information extraction model; as water bodies are also unlikely to undergo rocky desertification, they also need to be removed.Therefore, this study used the 2020 land use vector data of Huaxi District to mask built-Remote Sens. 2023, 15, 3056 5 of 20 up areas, bare land and water bodies.This was obtained from European Space Agency (ESA) World Cover 10 m land cover raster data (http://scihub.copernicus.eu/,accessed on 30 August 2022) in ArcGIS10.2using raster-to-vector conversion.The land cover raster data were produced from Sentinel-1 and Sentinel-2 data, with an overall accuracy of 85%.The masking process was completed in the Mask toolbox of ENVI5.3, and the required data sources are shown in Table 1.

Spectral Characteristics of Each Land Cover Type
As the most direct characterization factor of rocky desertification, the change in rocky exposure rate in remote sensing monitoring directly determines the change pattern of the degree of rocky desertification.Although there are many factors affecting rocky reflectance, such as the degree of rocky exposure, mineral composition and surface smoothness, the main factor affecting rocky reflectance for the same rocky type is the degree of rocky exposure, and its reflectance increases with the increase in rocky exposure.Based on the land cover characteristics of the study area, this study classified the land cover throughout the study area into five typical categories: vegetation, rocks, soils, buildings and water bodies.In this study, with the help of 1 m resolution Google Earth images, 50 pure image elements of each feature were manually and randomly selected on multispectral images of the Huaxi District to plot the mean Popper curves, as shown in Figure 2. As can be seen in Figure 2, the reflectance of the rocks increased from the red band (band 4: 630-690 nm) and the near-infrared band (band 5: 775-900 nm) to the short-wave infrared 1 band (band 6: 1550-1750 nm), and the reflectance of other features varied considerably in these three bands, so the combination of these three bands could be used to differentiate the features and further extract information on rocky desertification.In this study, the red, near-infrared and short-wave infrared bands, which have significant reflectance variations, were chosen as the bands for monitoring rocky desertification.

Methods
The Landsat 8 OLI remote sensing images of the study area were pre-processed with radiometric calibration, atmospheric correction and image cropping to obtain the surface reflectance of each band, and the reflectance of the red band (Red), near infrared band (NIR) and shortwave infrared band (SWIR) were used as feature variables.SWIR-NIR, Red-NIR and SWIR-Red reflectance spectral feature spaces were supplemented by lithology data, land use data, DEM elevation data, etc.The corresponding rocky desertification index models, PRDI1, PRDI2 and PRDI3, were constructed according to the changing characteristics of the degree of rocky desertification in each reflectance spectral

Methods
The Landsat 8 OLI remote sensing images of the study area were pre-processed with radiometric calibration, atmospheric correction and image cropping to obtain the surface reflectance of each band, and the reflectance of the red band (Red), near infrared band (NIR) and shortwave infrared band (SWIR) were used as feature variables.SWIR-NIR, Red-NIR and SWIR-Red reflectance spectral feature spaces were supplemented by lithology data, land use data, DEM elevation data, etc.The corresponding rocky desertification index models, PRDI1, PRDI2 and PRDI3, were constructed according to the changing characteristics of the degree of rocky desertification in each reflectance spectral feature space.The accuracy of these three index models extracted from rocky desertification was verified and compared, and the index model with the highest accuracy was compared with the previously proposed index models KRDI and RSDDI, which are constructed based on the surface parameter feature space.The overall flowchart of the study is shown in Figure 3, and some of the main components of the flowchart are described in the following sections.

Methods
The Landsat 8 OLI remote sensing images of the study area were pre-processed with radiometric calibration, atmospheric correction and image cropping to obtain the surface reflectance of each band, and the reflectance of the red band (Red), near infrared band (NIR) and shortwave infrared band (SWIR) were used as feature variables.SWIR-NIR, Red-NIR and SWIR-Red reflectance spectral feature spaces were supplemented by lithology data, land use data, DEM elevation data, etc.The corresponding rocky desertification index models, PRDI1, PRDI2 and PRDI3, were constructed according to the changing characteristics of the degree of rocky desertification in each reflectance spectral feature space.The accuracy of these three index models extracted from rocky desertification was verified and compared, and the index model with the highest accuracy was compared with the previously proposed index models KRDI and RSDDI, which are constructed based on the surface parameter feature space.The overall flowchart of the study is shown in Figure 3, and some of the main components of the flowchart are described in the following sections.OLI remote sensing image of the study area, SWIR-NIR was constructed, and its feature space is shown in Figure 4a.In the constructed SWIR-NIR spectral feature space, the pixels in the two-dimensional space were selected and the distribution of features corresponding to the pixels can be clearly seen on the remote sensing image map.As can be seen in Figure 4a, the pixels of different surface coverings were distributed in different locations in the feature space, with pixels of water bodies mainly in the blue oval area and pixels of buildings, roads, etc., mainly in the yellow oval area.Due to the different materials of the buildings, some of the pixels are darker while others are brighter, and can be observed that the pixels of the darker buildings are mainly distributed in the dark yellow ovals in Figure 4a, while the pixels of the brighter buildings are mainly distributed in the bright yellow ovals.In addition to these surface covers, some pixels of bare soil were found in the SWIR-NIR feature spatial scatter plot, presumably because the image was taken in late March, a season when it was under cultivation and there was a large amount of bare soil exposed.As areas such as water bodies, buildings, roads and bare ground are unlikely to be rocky desertified, they need to be removed to avoid the misclassification of rocky areas and masked to obtain the spectral feature space shown in Figure 4b.This spectral feature space was in an approximate triangular shape, and the pixels in the Figure were mainly vegetation, rocky or a mixture of vegetation and rocky.It can be observed that the bare rocky areas were mainly distributed in the purple oval area and the vegetation was mainly distributed in the lower border of the triangle, where the dark green oval area was mainly dark vegetation cover, such as coniferous forest, and the bright green oval area was mainly bright vegetation cover such as grassland.The vegetation from the dark vegetation cover to the bright vegetation cover went through the process of change from coniferous forest, broad-leaved forest, scrub and grassland, forming the concept of a vegetation line, denoted by L 1 , which was proposed by Xia Xueqi [39]  where ρNIR and ρSWIR represent the surface reflectance of the near-infrared band (band 5) and shortwave infrared band 1 (band 6) of the Landsat-8 OLI image, respectively; M is the slope of the vegetation line L1 and a is the intercept of the vegetation line on the perpendicular axis.The scatter plot shown in Figure 5a was obtained by constructing the feature space The scatter of image elements in geometric space can be abstracted into a single monomorph (trihedron or polyhedron), with pure elements (those covered by only one type of ground cover) corresponding to the vertices of the single monomorph, and mixed elements (those made up of two or more types of ground cover) distributed within or on the edges of the single monomorph [40].Figure 4c is a simulation of the SWIR-NIR feature space scatter plot.The pure image elements formed by bare rocky, light vegetation and dark vegetation are distributed at the three vertices of A, B and C, respectively, and the mixed image elements between them are distributed inside and on the edge of the triangle.It was found that there was a certain distribution pattern of rocky desertification image elements in the triangle: the closer the distance from point A, the more pure bare rock distribution, the less vegetation cover and the more serious rocky desertification; the closer the distance from the vegetation line L 1 , the more vegetation cover, the less pure bare rocky distribution and the less rocky desertification.It can be seen that more and more rocky areas were exposed along the direction perpendicular to the vegetation line, which means that the degree of rocky desertification is increasing, and that different degrees of rocky desertification can be better distinguished by an ellipse parallel to the vegetation line L 1 .Therefore, this study constructed a rocky desertification index based on the purity of the rocky desertification image element to the vegetation line, named the perpendicular rocky desertification index 1 (PRDI1).The following rocky desertification index was established based on the formula for the distance from the point to the line.
where ρN IR and ρSW IR represent the surface reflectance of the near-infrared band (band 5) and shortwave infrared band 1 (band 6) of the Landsat-8 OLI image, respectively; M is the slope of the vegetation line L 1 and a is the intercept of the vegetation line on the perpendicular axis.
Based on the SWIR-NIR reflectance spectral feature space established by the 2D Plot tool in ENVI 5.3, the minimum shortwave infrared reflectance value corresponding to each NIR reflectance value in the feature space was extracted based on the unique distribution of vegetation image elements in the 2D spectral feature space, and the trend line was obtained as the vegetation line; see Figure 4d.Therefore, the value of M in this study was 0.4338 and the value of a was 0.0047.Substituting the values of M and a into Equation (1) yields: 3.1.2.Model Construction Based on the Red-NIR Spectral Feature Space The scatter plot shown in Figure 5a was obtained by constructing the feature space for the red band (Red) and the near infrared band (NIR).It was found that the water bodies were mainly distributed in the blue ellipse area, the dark buildings were mainly distributed in the dark yellow ellipse area and the bright buildings were mainly distributed in the bright yellow ellipse area.The two-dimensional feature space scatter plot shown in Figure 5b was obtained after masking the areas where rocky desertification is unlikely to occur.The feature space scatter plot has a flat and long triangular shape due to the low reflectance of rocky areas and vegetation in the red wavelength band (Red).The purple ellipse in Figure 4b is mainly the bare rocky area, the lower line of the triangle is mainly the vegetation cover area, the dark green ellipse is mainly the dark vegetation cover area, such as coniferous forest, and the bright green ellipse is mainly the bright vegetation cover area, such as grassland.The vegetation from the dark vegetation cover area to the bright vegetation cover area goes through the process of changing from coniferous forest, broadleaved forest, scrub and grassland, forming a vegetation line, which is represented by L 2 .Figure 5c shows a simulation of the scatter plot of the Red-NIR feature space.Similarly, the spatial distribution of rocky desertification pixels in this feature space also has a certain regularity, with bare rocky distribution areas near point A, vegetation distribution areas near L 2 , and mixed vegetation and rocky pixels inside and at the edges of the triangle; an ellipse parallel to the vegetation line L 2 can better distinguish different degrees.The ellipse parallel to the vegetation line L 2 can better distinguish the different degrees of rocky desertification.Therefore, the perpendicular rocky desertification index 2 (PRDI2) can be constructed based on the degree of purity of the rocky desertification image element to the vegetation line L 2 , and the formula is established as follows.
where ρN IR and ρRed represent the surface reflectance of the near-infrared band (band 5) and red band (band 4) of the Landsat-8 OLI image, respectively; N is the slope of the vegetation line L 2 and b is the intercept of the vegetation line on the perpendicular axis.
distinguish different degrees.The ellipse parallel to the vegetation line L2 can better distinguish the different degrees of rocky desertification.Therefore, the perpendicular rocky desertification index 2 (PRDI2) can be constructed based on the degree of purity of the rocky desertification image element to the vegetation line L2, and the formula is established as follows.

PRDI2=
N×ρRed-ρNIR+b where  and  represent the surface reflectance of the near-infrared band (band 5) and red band (band 4) of the Landsat-8 OLI image, respectively; N is the slope of the vegetation line L2 and b is the intercept of the vegetation line on the perpendicular axis.Unlike the previous two spectral feature spaces, the spectral feature spaces constructed by SWIR and Red are narrower, but also approximate a triangular shape.As shown in Figure 6a, water bodies are mainly distributed in the blue oval region, and dark and light buildings are mainly distributed in the yellow oval region.The spectral feature space map shown in Figure 6b was obtained after masking those areas where rocky desertification does not occur.The lower left part of the figure near the origin was mainly the vegetation cover area, which is easily saturated due to the low reflectance of vegetation in the red band, making the distribution of vegetation in the feature space more concentrated [41].
Due to artificial mining, weathering and erosion, the rocky surface color changes to varying degrees [42], with darker colorful rocky areas and brighter rocky areas located in the dark purple oval and bright purple oval areas, respectively.It can be observed that there is also a certain distribution pattern of rocky desertification image elements in this spectral feature space (as shown in Figure 6c): the closer to the AB side, the more rock distribution and the more serious rocky desertification; the closer to the origin O, the more vegetation coverage and the less rocky desertification.The interior of the triangle is mainly a mixture of rock and vegetation image elements, and the ellipse perpendicular to the OP line can better distinguish different degrees of rocky desertification.Therefore, the distance from any point in the triangular feature space to the origin O (0, 0) is used to represent the process of rocky desertification.Based on the formula for the distance between two points, the expression for the rocky desertification index is as follows.
where ρRed and ρSW IR represent the surface reflectance of the red band (band 4) and shortwave infrared band 1 (band 6) of the Landsat-8 OLI image, respectively.

Calculation of Feature Space Model Based on Surface Parameters
To determine the reliability of the rocky desertification information extraction index model proposed in this paper, this study compares it with two previously constructed rocky desertification information extraction index models based on the feature space of surface parameters, namely, the karst rocky desertification index (KRDI) and the rocky desertification difference index (RSDDI).KRDI is a rocky desertification information extraction index model based on the three surface parameters of brightness, greenness and humidity of the tassel hat transformation, which has good sensitivity to different karst rocky desertification classes [20].The RSDDI is a model for extracting rocky desertification information based on two surface parameters: NDVI and Albedo.This index model

Calculation of Feature Space Model Based on Surface Parameters
To determine the reliability of the rocky desertification information extraction index model proposed in this paper, this study compares it with two previously constructed rocky desertification information extraction index models based on the feature space of surface parameters, namely, the karst rocky desertification index (KRDI) and the rocky desertification difference index (RSDDI).KRDI is a rocky desertification information extraction index model based on the three surface parameters of brightness, greenness and humidity of the tassel hat transformation, which has good sensitivity to different karst rocky desertification classes [20].The RSDDI is a model for extracting rocky desertification information based on two surface parameters: NDVI and Albedo.This index model can extract and classify rocky desertification information more accurately and conveniently [21], calculated using the following equations.
where B, G and W represent the brightness, greenness and moisture components of the tassel hat shift, respectively.Their calculation is referenced in the literature [43].
where NDV I is the normalized vegetation index and Albedo is the albedo of Landsat-8 OLI image.Their calculation is referenced in the literature [44].

Rocky Desertification Classification
According to Equations ( 2), ( 4) and ( 5)-( 7), this study used the Band Math tool of ENVI5.3 to obtain five rocky desertification indices, namely PRDI1, PRDI2, PRDI3, KRDI and RSDDI.Referring to the grading method of Guo, Bing and Wei, Haishuo et al. [18,45], these five rocky desertification indices were classified into five grades based on the natural breakpoint method in ArcGIS 10.2, namely, no rocky desertification (No), slight rocky desertification (Slight), moderate rocky desertification (Moderate), intense rocky desertification (Intense) and severe rocky desertification (Severe).The natural breakpoint method is a grading method obtained based on Jenks' optimization algorithm, which takes full account of the histogram distribution of the rocky desertification index and will set boundaries at locations with relatively large differences in data values, so that the differences in rocky desertification within each grade of the classification are minimized and the differences between grades of rocky desertification are maximized [46,47].

Accuracy Assessment
To compare the performance of the models in distinguishing different degrees of karstic desertification, a confusion matrix approach was used to evaluate the accuracy of the models.Four main accuracy evaluation elements were calculated and observed based on the confusion matrix, namely, production accuracy (P sm ), user accuracy (P ni ), overall classification accuracy (P 0 ) and kappa coefficient [48].Firstly, random sample points were created in ArcGIS 10.2, and 205 random sample points were obtained after deleting those falling in non-karst areas; the distribution of sample points is shown in Figure 1c.The coordinates of these 205 sample points were then imported into high-resolution Google Earth software, and a 30 × 30 m sized sample square was delineated at each sample point.Combined with the field survey, the degree of rocky desertification of each sample square was judged as the true value of rocky desertification, and finally the confusion matrix of model estimates and observations was established for each sample point.The overall precision, kappa coefficient, user precision and producer precision of each rocky desertification information mention model were calculated.

Comparison of the Accuracy of the Models
In this study, the confusion matrix for different degrees of rocky desertification the Huaxi District of China was studied using the perpendicular rocky desertificati index 1 (PRDI1) model as an example, as shown in Table 3.The accuracy of PRD PRDI2, PRDI3, KRDI and RSDDI models were compared by means of creating confusi

Comparison of the Accuracy of the Models
In this study, the confusion matrix for different degrees of rocky desertification in the Huaxi District of China was studied using the perpendicular rocky desertification index 1 (PRDI1) model as an example, as shown in Table 3.The accuracy of PRDI1, PRDI2, PRDI3, KRDI and RSDDI models were compared by means of creating confusion matrices, and the results are shown in Figure 8.According to the calculated statistics, the overall accuracy of PRDI1 was 0.829 and the Kappa coefficient was 0.784, which were higher than the other models.Moreover, the no rocky desertification, intensive rocky desertification and heavy rocky desertification values extracted by PRDI1 all had better user accuracy and producer accuracy, of above 84%, indicating that PRDI1 performed the best in the classification of rocky desertification in each model.The overall accuracy and Kappa coefficient of PRDI2 were only 0.254 and 0.073, respectively, and the producer and user accuracies for all five classes of rocky desertification were less than 60%, indicating that the index model was poor in both individual classes of rocky desertification and overall rocky desertification extraction, and that the model is not suitable for rocky desertification information extraction in the study area.The overall precision was 0.615, the Kappa coefficient was 0.512, and the only categories with producer precision greater than 80% were no rocky desertification and heavy rocky desertification.The producer precision and user precision of the other three grades of rocky desertification were lower, indicating that the model is only applicable to the extraction of no rocky desertification and heavy rocky desertification, and not to the extraction of slight rocky desertification, moderate no rocky desertification and severe rocky desertification.The overall precision of KRDI is 0.820 and the Kappa coefficient is 0.770, with good user precision and producer precision of 84% or more for no rocky desertification, intensive rocky desertification and severe rocky desertification, indicating that KRDI performs relatively well in rocky desertification classification and is more accurate for the extraction of this type of rocky desertification information in the study area.The overall precision of RSDDI extraction was 0.818, the Kappa coefficient was 0.767, and the user precision was higher than 84% for no rocky desertification, intensive rocky desertification and severe rocky desertification, and the producer precision was also higher than 78%, indicating that RSDDI performs better in the rocky desertification classification and is more accurate for the extraction of this type of rocky desertification information in the study area.The above comparison shows that the overall classification accuracy, Kappa coefficient, producer accuracy and user accuracy of PRDI1 are better than other models, and the extraction ability of rocky desertification is stronger for the overall and each class.It shows that PRDI1, an exponential model based on the SWIR-NIR feature space, has good performance and is feasible and applicable to the extraction of karstic rocky desertification information.

Spatial Distribution Characteristics of Rocky Desertification Extracted by PDRI1
The above study shows that the index model PRDI1, constructed based on the SWIR-NIR reflectance spectral feature space, had a higher accuracy than other models, and the spatial distribution of rocky desertification extracted by the index model was consistent with the actual distribution of rocky desertification in the field survey of the study area, which indicates that the method based on the SWIR-NIR reflectance spectral feature space is feasible and applicable for extracting the grade and spatial distribution of rocky desertification.The method based on the SWIR-NIR reflectance spectral feature space is feasible and applicable for extracting the grade and spatial distribution of rocky desertification.Therefore, in this study, PDRI1 was used to invert the spatial distribution of rocky desertification in the Huaxi District, and the results are shown in Figure 9.As can be seen in Figure 9, the serious rocky desertification areas are mainly concentrated in the western, southwestern, northwestern and southeastern parts of Huaxi District, mainly due to the high altitude, steep terrain, thin soil thickness and severe soil erosion in the above-mentioned areas, as well as the influence of human activities, resulting in more serious rocky desertification in these areas [49].In the eastern and central regions, the low altitude, gentle terrain, thick soil layer and lighter soil erosion, as well as the implementation of the project to return farmland to forest and grass and the comprehensive rocky desertification management project measures in place since 2005 and 2008, respectively, have allowed the vegetation ecosystem in the region to recover and soil erosion to be reduced, resulting in a lesser degree of rocky desertification [50].

Spatial Distribution Characteristics of Rocky Desertification Extracted by PDRI1
The above study shows that the index model PRDI1, constructed based on the SWIR-NIR reflectance spectral feature space, had a higher accuracy than other models, and the spatial distribution of rocky desertification extracted by the index model was consistent with the actual distribution of rocky desertification in the field survey of the study area, which indicates that the method based on the SWIR-NIR reflectance spectral feature space is feasible and applicable for extracting the grade and spatial distribution of rocky desertification.The method based on the SWIR-NIR reflectance spectral feature space is feasible and applicable for extracting the grade and spatial distribution of rocky desertification.Therefore, in this study, PDRI1 was used to invert the spatial distribution of rocky desertification in the Huaxi District, and the results are shown in Figure 9.As can be seen in Figure 9, the serious rocky desertification areas are mainly concentrated in the western, southwestern, northwestern and southeastern parts of Huaxi District, mainly due to the high altitude, steep terrain, thin soil thickness and severe soil erosion in the above-mentioned areas, as well as the influence of human activities, resulting in more serious rocky desertification in these areas [49].In the eastern and central regions, the low altitude, gentle terrain, thick soil layer and lighter soil erosion, as well as the implementation of the project to return farmland to forest and grass and the comprehensive rocky desertification management project measures in place since 2005 and 2008, respectively, have allowed the vegetation ecosystem in the region to recover and soil erosion to be reduced, resulting in a lesser degree of rocky desertification [50].

Sources of Error and Applicability of Each Model
The desirability of a constructed spectral index model is determined by its sensitivity to information about the target feature of interest [46].In a karst environment, the spectral reflectance of exposed bedrock, vegetation, exposed lime-rocky soils and construction sites show significantly different characteristics in SWIR of OLI images (Figure 3), with water, hydroxyl and carbonate rock being the main determinants of the spectral absorption characteristics in the SWIR band [51].Therefore, SWIR is one of the suitable bands for characterizing typical land cover types in karst areas [52].In addition, the combination of the Red, NIR and SWIR bands was finally identified as the basis for the new proposed index model due to the significant differences in spectral reflectance trends for different classes of features in the red to shortwave infrared bands.
In general, establishing a relatively optimal index is the key to extracting land cover [53].According to Feng, Haixia and Pei, Jie et al. [14,41], it can be seen that the NIR band is highly sensitive to vegetation, the SWIR band is highly sensitive to rocks, and the red band is relatively less sensitive to vegetation and rocks.The comparison in Figure 7 shows that PRDI1 was able to extract all levels of rocky desertification well, probably because the index was constructed from the NIR and SWIR bands, which are sensitive to

Sources of Error and Applicability of Each Model
The desirability of a constructed spectral index model is determined by its sensitivity to information about the target feature of interest [46].In a karst environment, the spectral reflectance of exposed bedrock, vegetation, exposed lime-rocky soils and construction sites show significantly different characteristics in SWIR of OLI images (Figure 3), with water, hydroxyl and carbonate rock being the main determinants of the spectral absorption characteristics in the SWIR band [51].Therefore, SWIR is one of the suitable bands for characterizing typical land cover types in karst areas [52].In addition, the combination of the Red, NIR and SWIR bands was finally identified as the basis for the new proposed index model due to the significant differences in spectral reflectance trends for different classes of features in the red to shortwave infrared bands.
In general, establishing a relatively optimal index is the key to extracting land cover [53].According to Feng, Haixia and Pei, Jie et al. [14,41], it can be seen that the NIR band is highly sensitive to vegetation, the SWIR band is highly sensitive to rocks, and the red band is relatively less sensitive to vegetation and rocks.The comparison in Figure 7 shows that PRDI1 was able to extract all levels of rocky desertification well, probably because the index was constructed from the NIR and SWIR bands, which are sensitive to vegetation and rocks, and the distribution of each scatter in the spectral feature space constructed based on them is more uniform, not concentrated in one place or scattered elsewhere, so the model is able to better extract PRDI2, which has better extraction ability for no, slight and moderate rocky desertification, but poorer extraction ability for intensive and severe rocky desertification, probably because PRDI2 is constructed based on the reflectance of red and near-infrared bands, which are sensitive to soil and vegetation, but less sensitive to rocks, resulting in the underestimation of intensive and heavy rocky desertification.The poor extraction ability of PRDI3 for heavy rocky desertification and the high extraction ability for other classes of rocky desertification may be due to the fact that PRDI3 is constructed from the short-wave infrared band and the red band, and although the short-wave infrared band is sensitive to rocks, the red band is less sensitive to both rocks and vegetation, making the extraction ability of heavy rocky desertification weaker.For KRDI and RSDDI, constructed based on surface parameter feature space, they are index models with high accuracy that have been verified by previous authors; the size and spatial distribution of rocky desertification extracted by them and PRDI1 were relatively similar, but the accuracy verification showed that the accuracy of PRDI1 proposed in this study was slightly higher than that of KRDI and RSDDI, thus indicating that the exponential model PRDI1 constructed based on the SWIR-NIR reflectance spectral feature space is feasible for karst rocky desertification information extraction.In addition, it can be seen from Figure 8 that the accuracy of each index model was relatively low in both mild and moderate rocky-deserted areas, due to the fact that there were more mixed pixels consisting of rocks, vegetation and bare ground in these areas, which leads to misclassification and the omission of pixels [54].
In addition to these factors, as the mathematical expressions and related parameters of PRDI1 and PRDI2 are determined based on a fixed vegetation baseline, which is only an ideal assumption, in reality the shape, width, thickness and other characteristics of the vegetation line are related to the type of vegetation, leaf water content, etc. Strictly speaking, the vegetation line is not a fixed line, and this may also cause errors in the model [39].In the feature space, there are a small number of images outside the triangle, which may be a mixture of uncensored buildings, bare ground, etc.The images of these features participate in the construction of the model, which brings some influence on the accuracy of the model, and the model should be further optimized in future research.
Although the feature space method based on surface parameters has the advantages of easy access to surface parameters, the ability to reflect changes in the hydrothermal combination of the surface cover and clear biophysical significance, we found in the course of our research that because the values of the surface parameters tend to change with the study area and the time season, the shape of the feature space constructed using them also changes, sometimes even in irregular polygons.It is difficult to form a feature space with a specific shape, such as a triangle or trapezoid, so this method is more difficult to apply.The method of constructing feature space based on wave reflectance can overcome these shortcomings, and the shape of the feature space constructed by it is relatively fixed, in the shape of an approximate triangle, which is not easily influenced by other factors.

Research Limitations
The study area was located in the karst mountains area of southwest China, with the most severe land degradation, uneven ground surface and complex topographic conditions, all of which increase the difficulty of accurately extracting information on karstic rocky desertification [28].At present, there are three main difficulties in extracting information on rocky desertification: firstly, as the spectral characteristics of soil and rocky are very similar, it is very easy to confuse rocky and soil on remote sensing images, and the presence of soil will directly affect the accuracy of rocky desertification information extraction.This affects the accuracy of the results.Secondly, due to the large variation in topography, there are some backlit areas which are also affected by reflected radiation between adjacent objects, resulting in a certain number of shadow areas on the satellite images, and these shadows may also lead to uncertainty in the estimation of rocky desertification [16]; future studies should consider the effect of shadows.Thirdly, as remote sensing images are flat images after projection, and the topography of karstic rocky desertification areas is undulating, it is difficult to obtain the results of the field survey of rocky desertification after projection.Coupled with the irregularity of the distribution of rocks and vegetation in the field, it is difficult to calculate quantitatively the proportion of rocks and vegetation within the sampling units, which brings great difficulties to the accuracy verification of the results of rocky desertification information extraction.
In summary, the index PRDI1, constructed based on the SWIR-NIR reflectance spectral feature space, takes advantage of the two-dimensional spatial information to effectively characterize the dynamic changes of ground cover types and rocky desertification with clear biophysical significance, and the method is simple, effective and easy to access and operate.In addition, the question of how to improve the performance of the index PRDI1 by considering the influence of factors such as shading and soil and rocky confusion on the accuracy of PRDI1 will be the focus of future work.

Conclusions
Rocky desertification has caused serious ecological and socio-economic problems in the karst region of southwest China, with far-reaching effects, while the rugged terrain and complex topography of the karst region make it extremely difficult to extract rocky desertification information over large areas.To this end, this study proposes a simple method to effectively and automatically extract information on rocky desertification in the karst region of southwest China.The method firstly constructed two-dimensional SWIR-NIR, Red-NIR and SWIR-Red reflectance spectral feature spaces using the SWIR, NIR and Red reflectances of Landsat-8 OLI remote sensing imagery based on the analysis of the spectral response characteristics of the main land cover types.Secondly, three new rocky desertification information extraction index models, PRDI1, PRDI2 and PRDI3, were developed based on the point-to-point and point-to-line principles.Finally, the three index models were compared with the index models KRDI and RSDDI, which were constructed based on the surface parameter feature space.The results show that the index model constructed by this study based on the wave reflectance feature space method had the advantages of easy access to data, simple calculation and good stability, which were convenient for the monitoring and quantitative evaluation of rocky desertification in karst areas.In addition, the index model PRDI1 established in this study had an overall accuracy of 0.829 and a Kappa coefficient of 0.784, which was more accurate than other index models, and its performance was good.The overall accuracy of PRDI1 was 0.9% and 1.1% higher than that of KRDI1 and RSDDI, and was more accurate and effective in extracting information on the spatial distribution and extent of rocky desertification.Therefore, this study provides important information on the total area of rocky desertification and the spatial distribution of each grade of rocky desertification in Huaxi District.These findings can provide basic research data for the monitoring and assessment of karst rocky desertification, which can be used to support the local government's ecological and environmental management decisions.

Figure 1 .
Figure 1.(a) Location map of the study area.(b) Elevation map of the Huaxi District.(c) Distribution map of sample points in Huaxi District.

Figure 1 .
Figure 1.(a) Location map of the study area.(b) Elevation map of the Huaxi District.(c) Distribution map of sample points in Huaxi District.

Figure 2 .
Figure 2. Spectral profiles of the five main land cover types in the study area.

Figure 2 .
Figure 2. Spectral profiles of the five main land cover types in the study area.

Figure 2 .
Figure 2. Spectral profiles of the five main land cover types in the study area.

Figure 3 .
Figure 3.The flow chart of research method.Figure 3. The flow chart of research method.

Figure 3 .
Figure 3.The flow chart of research method.Figure 3. The flow chart of research method.

3. 1 .
Model Construction Based on Reflectance Spectral Feature Space 3.1.1.Model Construction Based on the SWIR-NIR Spectral Feature Space Using band 5 (NIR) and band 6 (SWIR) of the atmospherically corrected Landsat 8

Figure 4 .
Figure 4. (a) SWIR-NIR spectral feature space before being masked.(b) SWIR-NIR spectral feature space after being masked.(c) Distribution of rocky desertification classes in the SWIR-NIR feature space: No-no rocky desertification; Slight-slight rocky desertification; Moderate-moderate rocky desertification, Intensive-intensive rocky desertification; Severe-severe rocky desertification.(d) Fitted curve of vegetation line L1 in the SWIR-NIR spectral feature space.

Figure 4 .
Figure 4. (a) SWIR-NIR spectral feature space before being masked.(b) SWIR-NIR spectral feature space after being masked.(c) Distribution of rocky desertification classes in the SWIR-NIR feature space: No-no rocky desertification; Slight-slight rocky desertification; Moderate-moderate rocky desertification, Intensive-intensive rocky desertification; Severe-severe rocky desertification.(d) Fitted curve of vegetation line L 1 in the SWIR-NIR spectral feature space.

Figure 5 .Figure 5 . 3 .
Figure 5. (a) Red-NIR spectral feature space before being masked.(b) Red-NIR spectral feature space after being masked.(c) Distribution of rocky desertification classes in Red-NIR spectral fea-Figure 5. (a) Red-NIR spectral feature space before being masked.(b) Red-NIR spectral feature space after being masked.(c) Distribution of rocky desertification classes in Red-NIR spectral feature space: No-no rocky desertification; Slight-slight rocky desertification; Moderate-moderate rocky desertification, Intensive-intensive rocky desertification; Severe-severe rocky desertification.(d) Fitted curve of vegetation lines L 2 in Red-NIR spectral feature space.Based on the Red-NIR reflectance spectral feature space established by the 2D Plot tool in ENVI5.3, the minimum red reflectance value corresponding to each NIR reflectance value in the feature space was extracted based on the unique distribution of vegetation image elements in the 2D spectral feature space, and the trend line was obtained as the vegetation line; see Figure 5d.Therefore, the value of N in this study is 0.1667 and the value of a is 0.0015.Substituting the values of N and b into Equation (3) yields: PRDI2 = |0.1667× ρRed − ρN IR + 0.0015| √ 0.1667 2 + 1 (4) 3.1.3.Model Construction Based on the SWIR-Red Spectral Feature Space

Figure 6 .
Figure 6.(a) SWIR-Red spectral feature space before being masked.(b) SWIR-Red spectral feature space after being masked.(c) Distribution of rocky desertification classes in the SWIR-Red feature space: No-no rocky desertification; Slight-slight rocky desertification; Moderate-moderate rocky desertification, Intensive-intensive rocky desertification; Severe-severe rocky desertification.

Figure 6 .
Figure 6.(a) SWIR-Red spectral feature space before being masked.(b) SWIR-Red spectral feature space after being masked.(c) Distribution of rocky desertification classes in the SWIR-Red feature space: No-no rocky desertification; Slight-slight rocky desertification; Moderate-moderate rocky desertification, Intensive-intensive rocky desertification; Severe-severe rocky desertification.

Figure 9 .
Figure 9. (1) Results of rocky desertification extracted using PDRI1 (the red squares represent the different levels of rocky desertification).(2) Google Earth remote sensing images (the red dots represent the location of the different levels of rocky desertification on the Google Earth image).(3) Field photographs: (a) no rocky desertification; (b) slight rocky desertification; (c) moderate rocky desertification; (d) intensive rocky desertification; and (e) severe rocky desertification.

Figure 9 .
Figure 9. (1) Results of rocky desertification extracted using PDRI1 (the red squares represent the different levels of rocky desertification).(2) Google Earth remote sensing images (the red dots represent the location of the different levels of rocky desertification on the Google Earth image).(3) Field photographs: (a) no rocky desertification; (b) slight rocky desertification; (c) moderate rocky desertification; (d) intensive rocky desertification; and (e) severe rocky desertification.

Table 3 .
Confusion matrix, using PRDI1 as an example.