Analysis of Landsat-8 OLI Imagery for Estimating Exposed Bedrock Fractions in Typical Karst Regions of Southwest China Using a Karst Bare-Rock Index

Abstract

Karst rocky desertification (KRD) has become the primary ecoenvironmental problem in the karst regions of southwest China. The rapid and efficient acquisition of exposed bedrock fractions (EBF) is crucial for the monitoring and assessment of KRD degree and distribution within the highly heterogeneous landscapes. Remote-sensing indices provide a useful method for the quick mapping of the EBF at large scales. The currently available rock indices, however, are faced with insensitivity to bedrock change characteristics, which greatly limits their performances and suitability. To address this problem, we proposed a novel karst bare-rock index (KBRI) that applies shortwave-infrared (SWIR) and near-infrared (NIR) bands from Landsat-8 OLI imagery to maximally distinguish between exposed bedrock and other land cover types in southwest China. A linear regression model was thus established between KBRI and the EBF derived from in situ measurements. The model developed here was then validated with an independent experiment and applied over a large geographic area to produce regional maps of EBF in southwest China. Experimental results showed good performance on root mean square error (5.59%), mean absolute error (4.63%), root mean absolute percentage error (13.59%), and coefficient of determination (0.72), respectively. The advantages of the proposed method are reflected in its simplicity and minimal requirements for auxiliary data while still achieving comparatively better accuracy than existing related indices. Thus, the KBRI has the great potential for the application in other regions around the world with the similar geological backgrounds, thereby helping to address the similar or other related environmental issues. Results of this study provide baseline data for the KRD assessment and karst-ecosystem management in southwest China.

1. Introduction

The karst area in southwest China is located in the center of the East Asia area, one of the three largest continuous distributions of karst development zones in the world, covering about 540,000 km² and supporting more than 220,000,000 people [1,2]. Due to its small carrying capacity of an ecological environment and a large population, the sharp contradiction between humans and land resources has led to a unique type of desertification, namely, karst rocky desertification (KRD). It is a special type of land degradation, in which a karst area covered by vegetation and soil is transformed into a bare rocky landscape with very limited vegetation and soil [2,3]. This has become the source of environmental disasters (e.g., droughts, floods, landslides, etc.) and economic poverty in the karst area of southwest China [4,5]. When KRD develops, the ecosystem is involved in reverse succession, during which the original vegetation is degraded or disappears, leading to bedrock exposure and low vegetation cover as the most obvious land-surface symptoms [6,7]. Therefore, the fractional coverage of exposed bedrock and vegetation is commonly represented as the key ecological indicators of KRD [7,8,9]. To date, extensive research has been conducted on estimation of fractional vegetation cover [10,11,12]. Consequently, the rapid and efficient acquisition of exposed bedrock fractions (EBF) is thus crucial for the monitoring and assessment of KRD degree and distribution at the local and regional scale.

Although traditional ground surveys can generate accurate information related to the EBF, in situ observations are usually expensive, time-consuming, labor-intensive and spatially limited [13,14], and thus fail to conduct multitemporal comparison and multiscale research [15]. By contrast, remote-sensing technology has the great advantages of low cost of data acquisition, large-area coverage, long time-series observations, and spatial continuity [16,17], rendering it an indispensable means for the quantitative extraction of EBF information in a large area. Early work with satellite imagery focused principally on visual interpretation and computer-assisted digital image processing [8,9], which was lack of objectivity and efficiency, and it was susceptible to the interpreter’s bias [18]. The karst ecosystem in southwest China is famous for its high degree of landscape heterogeneity and broken surface [19]. Consequently, land covers are often mixtures of several types in karst regions [20], making it quite difficult to accurately extract the main land-surface symptoms of KRD. Vegetation indices and spectral mixture analysis (SMA) have been widely used to extract land cover fractions applied on a subpixel scale by remote sensing [21,22,23]. Nevertheless, the applications of these two methods in the research of karst feature extraction have some limitations. On the one hand, vegetation indices are mainly based on the spectral reflectance characteristics of green vegetation and can thus reflect the greenness information of vegetation [24]. In other words, vegetation indices are not designed for directly detecting the EBF information. On the other hand, the SMA method makes a basic assumption that mixed pixels in which land surface is composed of a limited number of primary features (i.e., endmembers), and the observed mixed reflectance is the function of an endmember’s spectrum and its corresponding abundance [25,26]. However, its practical use is limited by the difficulty of endmember selection and endmember variability [27]. For the karst regions in southwest China, with high landscape complexity and discontinuous landform distribution, one pixel often contains mixed spectral information of vegetation, bedrock, and bare soil. Therefore, it is very difficult to select the endmembers and thus not easily applicable to employ the SMA method to obtain land-cover information in highly heterogeneous landscapes of karst regions [4,7]. It is thus necessary to develop a relatively simple and efficient approach to quantify the EBF in karst areas.

Hyperspectral remote-sensing imagery, consisting of hundreds of continuous narrow bands obtained by imaging spectrometry [8,28], can provide abundant diagnostic spectral channels to enhance the ability to distinguish different land covers in KRD areas. Recently, a KRD synthesis index (KRDSI) was proposed to extract the fractional coverage of nonvegetation land-cover types (i.e., nonphotosynthetic vegetation, exposed bedrock, and bare soil) based on spectral reflectance differences between different land-cover types using the short-wave infrared (SWIR 2000–2400 nm) of hyperspectral images [4,7]. However, the application of hyperspectral remote sensing is limited by its inherent contradiction between the spectral resolution and the signal-to-noise ratio [29]. In addition, the relatively high cost and narrow swath width render its practicality limited to the small-scale land cover extraction. In this regard, the multispectral images obviously have more advantages [18]. Huang and Cai [17] firstly established a normalized difference rock index (NDRI) using a Landsat ETM image and supervised classification method to map karst-exposed bedrock. However, this method could not differentiate well between exposed bedrock and built-up land. Xie et al. [18] proposed two kinds of carbonate rock indices (CRIs) using Landsat-8 OLI imagery, based on the unique spectral response of karst rock and other land cover types. Nevertheless, given that the experimental site was relatively small and the study didn’t involve any analysis of the spectral reflectance comparison between karst bedrock and built-up land, it still remains questionable whether the CRIs can be well applicable to distinguish exposed bedrocks in built-up areas at a large scale. Additionally, Li and Wu [13] constructed a new NDRI by using the shortwave-infrared (SWIR) and near-infrared (NIR) band deployed on Landsat-8 OLI imagery. In their study, a Dimidiate Pixel Model was employed to estimate the bedrock exposure rate, assuming that a pixel was proportionally composed of two components, namely, rock and vegetation. However, the model itself ignored the basic fact that the mixed pixels in karst areas with large heterogeneity consist of multiple typical land features [8]. Moreover, most rock-related indices are site-specific and thus not well applicable when altering the study areas.

To solve these problems, we propose a simple mathematical approach rather than a mechanism model that could be used to extract the EBF information directly and automatically within complex landscapes. The principal purposes of this study are: (1) to develop a new spectral index sensitive to the EBF by investigating the unique spectral features of exposed bedrock and other typical land-cover types, using Landsat-8 OLI imagery; (2) to establish a linear regression model between the proposed index and the EBF derived from in situ measurements to quantify the EBF values for the whole study area; (3) to conduct the comparison analysis between the new index and existing relevant indices in terms of EBF estimation, in order to quantify the effectiveness of the newly proposed index in this study; (4) to validate the robustness of the developed model with independent experimental data.

2. Materials and Methods

2.1. Study Area

The Xiaojiang watershed is located in eastern Yunnan Province, southwest China (24°10′–24°45′N, 103°30′–104°05′E; Figure 1), which has been suffering from serious KRD for decades due to the fragile karst ecosystem and serious conflicts between human beings and land resources. It is a first-level tributary of the Nanpan River, with its main part located in Luxi County, Honghe Hani, and Yi Autonomous Prefecture. The total area of the whole watershed attains 1008.39 km², with 80.36% of the area covered by karst environment (Figure 1c). The Xiaojiang watershed is characterized by a subtropical plateau monsoon climate, with distinct dry and rainy seasons [30]. The mean annual temperature in the area is 15.2 °C, the mean annual rainfall is 966.8 mm, and over 80% of precipitation is concentrated between June and October [30].

The Nandong underground river watershed (hereinafter referred to as the Nandong watershed) was treated as an independent validation area for model performance evaluation in this study. It is located in the southeast of Yunnan (23°10′–23°43′N, 103°10′–103°42′E; Figure 1), which is a typical super-large underground river watershed in southwest China [31], and also known as a serious KRD region [32]. The total area of the Nandong watershed is 1617.13 km², in which 88.74% of the whole area is covered by karst landscape (Figure 1d). The watershed is located on the Tropic of Cancer and characterized by a subtropical monsoon climate, with a mean annual temperature of 19.2 °C and mean annual rainfall of 830 mm [32].

2.2. Landsat-8 Data and Preprocessing

The Landsat-8 satellite was launched on 11 February 2013, equipped with two sensor payloads, i.e., Operational Land Imager (OLI) and Thermal Infrared Sensor (TIRS), collecting images of the Earth with a 16-day repeat cycle and a scan swath of 185 km. The Landsat archive was reorganized into a tiered data Collection structure in 2016 [33]. In this study, Landsat-8 Collection 1 Level 1 Precision and Terrain (L1TP)-corrected data were used. The L1TP products are composed of eight multispectral bands with a spatial resolution of 30 m, one panchromatic band with a resolution of 15 m, and two thermal bands collected at 100 m, but resampled to 30 m to match the multispectral bands [34,35]. Given that, during the nongrowing season, especially in January, deciduous trees generally shed leaves and grass withers, exposed bedrock easily appears in karst regions of southwest China [13]. Therefore, in this study, two cloud-free OLI images (with georeference information of UTM zone 48 N, WGS84) covering the Xiaojiang and Nandong watershed were acquired separately on 9 January and 14 March 2017, downloaded from https://earthexplorer.usgs.gov/. The detailed description of Landsat-8 OLI images is shown in

Table 1. Landsat-8 Operational Land Imager (OLI) images used in this study.

Acquisition Date	Scene ID	Path	Row	Sun Azimuth (°)	Sun Elevation (°)	Cloud Amount (%)
9 January 2017	LC81290432017009LGN01	129	43	150.58	37.48	0.61
14 March 2017	LC81290442017073LGN00	129	44	132.15	54.14	0.14

.

The preprocessing of Landsat OLI images was conducted according to the following steps. First, the images were calibrated to at-sensor radiance images. Second, the radiance images were then corrected to surface-reflectance images using the Fast Line-of-sight Atmospheric Analysis of Hypercubes (FLAASH) tools based on a MODTRAN radiative transfer model embedded in the Environment for Visualizing Images (ENVI) software (version 4.8) (Boulder, CO, USA) [35]. Third, all images were cropped to cover the study areas using the watershed boundary vector data. Moreover, due to the aforementioned high landscape heterogeneity in the study areas, in order to guarantee more pure pixels of exposed bedrock, bare soil, vegetation, and so forth in the subsequent spectral analysis, the Gram–Schmidt Pan Sharpening tool in ENVI software was thus employed to elevate the resolution of OLI multispectral image to 15 m using the panchromatic band of L1TP image.

2.3. Auxiliary Data and Accuracy Assessment

Auxiliary data used in this study mainly refer to field-survey data, consisting of 340 quadrats, each measuring a 30 m square. Among them, 264 quadrats, investigated during 20–25 January 2017, were used to provide ground-truth information of the EBF in the Xiaojiang watershed, to establish the model relating the proposed spectral index to the EBF. To validate the model robustness, another 76 quadrats, investigated during 8–10 March 2017, were used to derive the in situ measurements of EBF in the Nandong watershed, which were then compared with the EBF obtained from the developed model. Additionally, a handheld global-positioning satellite (GPS) was used to record the geographical coordinates of these quadrats.

Several indicators were introduced in this paper in order to quantitatively assess the accuracy of the developed model, including the root mean square error (RMSE), mean absolute error (MAE), root mean absolute percentage error (RMAPE), and coefficient of determination (R²). These assessment indicators were widely used in model evaluation studies and their mathematical expressions were described as the following equations [36,37,38]:

$${\delta }_{\mathrm{RMSE}}=\sqrt{\frac{1}{n}{\displaystyle \sum _{i=1}^n{\left(y_i-{\widehat{y}}_i\right)}^2}}$$

(1)

$${\delta }_{\mathrm{MAE}}=\frac{1}{n}{\displaystyle \sum _{i=1}^n\left|y_i-{\widehat{y}}_i\right|}$$

(2)

$${\delta }_{\mathrm{RMAPE}}=\frac{1}{n}\sqrt{{\displaystyle \sum _{i=1}^n{\left[\left(y_i-{\widehat{y}}_i\right)/y_i\right]}^2}}$$

(3)

$${\mathrm{R}}^2={\left({\displaystyle \sum _{i=1}^n\left({\widehat{y}}_i-\overline{\widehat{y}}\right)\times \left(y_i-\overline{y}\right)}/\left(\sqrt{{\displaystyle \sum _{i=1}^n{\left({\widehat{y}}_i-\overline{\widehat{y}}\right)}^2}}\times \sqrt{{\displaystyle \sum _{i=1}^n{\left(y_i-\overline{y}\right)}^2}}\right)\right)}^2$$

(4)

where n is the total amount of samples; $y_i$ and ${\widehat{y}}_i$ represent the observation value and estimated value, respectively; and $\overline{y}$ and $\overline{\widehat{y}}$ denote the average observation value and average estimated value, respectively. When RMSE, MAE, and RMAPE are low and R² is high, then the best-fit model is achieved [39].

2.4. The Karst Bare-Rock Index (KBRI)

Since the vegetation/soil cover and bedrock exposure are basically the classification and assessment standard for the KRD degree [9], we focused more on these three land-cover types. Meanwhile, due to the large terrain differences and broken land surfaces in the study area, there existed a certain amount of crescent and stick shape shadow on the remotely sensed imagery [40], also recognized as the information-missing area. Thus, in this study, land covers within the whole watershed were divided into six typical classes, vegetation, rock, soil, shadow, built-up, and water. It is worth noting that vegetation is actually the mixture of cropland, forestland, and grassland, while the built-up includes urban and rural residential land, roads, and other built-up land. With the assistance of Google Earth imagery with a resolution of 1 m, 50 pure pixels for each class were randomly selected manually on the multispectral image of Xiaojiang watershed to depict the spectral profiles.

As shown in Figure 2, the spectral difference between bare rock and other land-cover types is apparent at SWIR bands, especially SWIR1 (band 6 in OLI image). Some scholars used this feature to propose related indices for extracting karst rock [13,17]. Furthermore, it is noticed that the surface reflectance of rock, soil, and built-up appears a rising trend from band 5 (NIR) to band 6 (SWIR1) whereas that of vegetation, shadow, and water shows a decreasing trend. Meanwhile, the slope of increase of reflectance of rock is obviously larger than that of soil and built-up. Consequently, the band combination of SWIR1 and NIR can be employed to identify the exposed bedrock and other types in the study area. To achieve a higher level of reflectance contrast between rock and other land cover types, a root function was applied on the new spectral index [41]. In order to restrict the index value range (–1~1), the equation was divided by 20. Therefore, the KBRI is calculated from the satellite image using the following formula:

$$\mathrm{KBRI}=\frac{\rho SWIR1-\rho NIR}{20\times \sqrt{\rho SWIR1+\rho NIR}}$$

(5)

where $\rho SWIR1$ and $\rho NIR$ represent the surface reflectance of SWIR1 and NIR band of Landsat-8 OLI image, respectively.

2.5. Comparison with Related Indices

To determine the effectiveness of the new index proposed in this paper, we also developed two simple rock indices (SRI) using the relatively simple band ratios deployed on the same Landsat-8 OLI imagery. For the convenience of differentiation, the two indices were separately named SRI1 and SRI2, with their expressions defined as the following equations:

$$\mathrm{SRI}1=\rho Green/\rho SWIR2$$

(6)

$$\mathrm{SRI}2=2\times \rho SWIR1/(\rho NIR+\rho SWIR2)$$

(7)

where $\rho Green$, $\rho NIR$, $\rho SWIR1$, and $\rho SWIR2$ represent the surface reflectance of green (band 3), NIR (band 5), SWIR1 (band 6), and SWIR2 (band 7) of the Landsat-8 OLI image, respectively.

Additionally, experiments on other related indices were also conducted in a comparison analysis, including the two carbonate-rock indices (i.e., CRI1 and CRI2) [18], NDRI from Landsat ETM imaging [17], and NDRI deployed on Landsat-8 OLI imaging [13]. It is worth noting that although the latter two indices were designated with the same name, they were established with different band combinations. Thus, to avoid confusion, the former index was renamed as NDRI1 and the latter as NDRI2. Furthermore, the normalized-difference built-up index (NDBI) [42], and the normalized-difference vegetation index (NDVI) were also introduced in the comparison analysis. The formulas of these related indices were shown as follows:

$$\mathrm{CRI}1=\rho Blue/\rho NIR$$

(8)

$$\mathrm{CRI}2=(\rho Blue-\rho NIR)/(\rho Blue+\rho NIR)$$

(9)

$$\mathrm{NDRI}1=(\rho SWIR1-\rho Red)/(\rho SWIR1+\rho Red)$$

(10)

$$\mathrm{NDRI}2=(\rho SWIR2-\rho NIR)/(\rho SWIR2+\rho NIR)$$

(11)

$$\mathrm{NDBI}=(\rho SWIR1-\rho NIR)/(\rho SWIR1+\rho NIR)$$

(12)

$$\mathrm{NDVI}=(\rho NIR-\rho Red)/(\rho NIR+\rho Red)$$

(13)

where $\rho Blue$, $\rho NIR$, $\rho SWIR1$, $\rho Red$, and $\rho SWIR2$ represent the surface reflectance of blue (band 2), NIR (band 5), SWIR1 (band 6), red (band 4), and SWIR2 (band 7) of Landsat-8 OLI imaging, respectively.

2.6. A Linear Regression Model

In the karst geological environment, exposed bedrock, bare soil, and vegetation have their own unique spectral-absorption characteristics [4]. Here, we hypothesized that, with the increase of fractions of exposed bedrock, bare soil, and vegetation, their corresponding spectral responses will be more prominent [4]. In other words, the proposed index could be used to effectively capture the spectral-absorption feature, and thus be adopted to retrieve the EBF information. The observed EBF derived from field surveys was related to the KBRI values, using a linear regression model, after verifying that a linear relationship existed. Hence, this linear-regression model was then constructed using the Statistical Product and Service Solutions (SPSS) software package (version 16.0) (Armonk, NY, USA), and used to estimate the EBF over the entire Xiaojiang watershed, southwest China, using satellite information (KBRI-derived values). The model can be mathematically defined as the following equation:

F_rock = a + b × KBRI

(14)

where F_rock represents the EBF, with a and b as constant coefficients. Statistically significant differences were set with p values = 0.05. After the coefficients of a and b is determined, the linear regression model can be applied to the whole OLI image. A full map of EBF throughout the Xiaojiang watershed can thus be generated. Here, 264 in situ measurements of EBF across the Xiaojiang watershed were then compared with the EBF values obtained from the linear-regression model using the assessment indicators listed in Section 2.3.

To illustrate the effectiveness of the proposed KBRI index, other related indices listed in Section 2.5 also had their linear relationships with EBF explored using the same method described above, except using the corresponding index instead of KBRI deployed on the same satellite image.

In addition, to validate the robustness of the developed models, an independent validation experiment was conducted in the Nandong watershed, located 60 km southwest of the Xiaojiang watershed, with similar landscape characteristics. In this test experiment, 76 additional samples were collected throughout Nandong to be compared with the EBF information obtained from the original model, using the quantitative assessment indicators.

3. Results

3.1. Index Images

All indices applied in this study, including KBRI, SRI1, SRI2, CRI1, CRI2, NDRI1, NDRI2, NDBI, and NDVI, were extracted from the Landsat-8 OLI image. The differences between these index images can be easily demonstrated by visual judgment. Figure 3 presents their visual differences in the Xiaojiang watershed on 9 January 2017.

As shown in Figure 3b–d,i, the index images had a relatively small and concentrated range of pixel values, intuitively manifested as no significant difference between the pixel colors, suggesting that CRI1, CRI2, NDRI1, and SRI2 are not sensitive to the EBF changes and thus not conducive to establishing a well-fit regression model between the EBF and the corresponding index. In the NDRI2 image (Figure 3e), the exposed bedrock was confused with bare soil. Likewise, in the NDBI image (Figure 3f), the difference between exposed bedrock and built-up land was not large enough to achieve a satisfactory classification. However, in the KBRI image (Figure 3a), the difference between the exposed bedrock and all other land-cover types was so clear that the exposed bedrock could be easily interpreted.

3.2. Mapping of EBF in Xiaojiang Watershed

In order to estimate the EBF in a karst-dominated environment, a linear-regression model was constructed. Four indicators, RMSE, MAE, RMAPE, and R², were employed to evaluate the performance of the developed models. Figure 4 shows the linear relationships between EBF and the corresponding index in the Xiaojiang watershed.

As shown in Figure 4, most indices, including KBRI, CRI1, CRI2, NDRI2, NDBI, and SRI2, had a significant positive correlation with the EBF, with different regression coefficients. Nevertheless, NDRI1, NDVI, and SRI1 showed a significant negative correlation with the EBF. It is noticed that the linear correlation between the EBF and the KBRI results showed the highest degree of correlation (Figure 4a). Meanwhile, it is easily illustrated that the KBRI values ranged more widely compared with other spectral indices, which conformed to the actual conditions in the study area. In contrast, except for KBRI, SRI1, and SRI2, the scatterplots of six other indices clearly showed that the phenomenon of the same abscissa corresponding to multiple ordinates has occurred many times (Figure 4b–g), indicating that these spectral indices are not sensitive to the EBF variation characteristics.

Table 2. Accuracy assessment of the linear regression models in the Xiaojiang and Nandong watershed.

Indices	Xiaojiang				Nandong
	Root Mean Square Error (RMSE) (%)	Mean Absolute Error (MAE) (%)	Root Mean Absolute Percentage Error (RMAPE) (%)	Coefficient of Determination (R2)	RMSE (%)	MAE (%)	RMAPE (%)	R2
KBRI	5.59	4.63	13.59	0.72	8.62	7.01	5.93	0.70
CRI1	7.92	6.12	26.04	0.43	19.55	15.42	52.19	0.13
CRI2	7.94	6.20	24.99	0.43	18.21	14.46	50.89	0.15
NDRI1	10.05	7.45	33.60	0.09	14.74	10.85	37.20	0.00
NDRI2	5.90	4.77	14.87	0.68	11.44	8.57	33.03	0.45
NDBI	5.86	4.85	12.31	0.69	11.11	8.16	32.66	0.43
NDVI	7.82	5.99	26.15	0.45	15.59	11.87	38.13	0.16
SRI1	11.24	8.57	35.51	0.18	13.91	11.04	8.81	0.26
SRI2	39.44	38.46	59.06	0.51	10.62	8.47	6.43	0.40

presented the result of accuracy assessment of the developed models. It can be easily seen that KBRI owned the lowest RMSE and MAE, the second-lowest RMAPE, and the highest R² (5.59%, 4.63%, 13.59%, and 0.72, respectively) among all nine indices. Although NDBI produced a lower RMAPE value than KBRI, the RMSE and MAE values of NDBI were higher than that of KBRI, while the R² value of NDBI was lower than that of KBRI. Generally, the result of assessment accuracy indicated that KBRI, the new spectral index proposed in our study, has the best performance relative to other remote sensing indices. Conversely, SRI2 achieved the highest RMSE, MAE, and RMAPE values, with a relatively low R² (39.44%, 38.46%, 59.06%, and 0.51, respectively). Meanwhile, SRI1 owned the second-highest RMSE, MAE, and RMAPE values, and the second-lowest R² value (11.24%, 8.57%, 35.51%, and 0.18, respectively). In addition, NDRI1 presented the third-highest RMSE, MAE, and RMAPE values, following SRI2 and SRI1, and the lowest R² value (10.05%, 7.45%, 33.60%, and 0.09, respectively). Therefore, all facts suggest that the SRI1, SRI2, and NDRI1 have the worst performances and thus are not fit for the EBF estimation. Hence, estimating the EBF using a linear-regression model developed with KBRI, as proposed in this study, is feasible and reliable. Figure 5 illustrates the full maps of EBF throughout the Xiaojiang watershed produced via the developed models with the corresponding index.

According to the KRD classification standard [8], the EBF was categorized as five different levels, corresponding to no KRD, slight KRD, moderate KRD, severe KRD, and extreme KRD, respectively.

Table 3. Statistical results of EBF in the Xiaojiang and Nandong watershed.

Region	EBF Level	<5%	5–15%	15–30%	30–50%	>50%
Xiaojiang	Proportion (%)	29.73	14.91	26.99	26.77	1.59
	Area (km2)	299.84	150.35	272.21	269.92	16.06
Nandong	Proportion (%)	10.82	9.84	24.18	44.42	10.75
	Area (km2)	174.94	159.06	390.94	718.34	173.84

showed the statistical results of EBF estimated by the linear-regression model developed with KBRI. For each EBF level, the area proportion was calculated by dividing the area of this EBF level by the total area of the watershed. The results showed that in the Xiaojiang watershed, the proportion of pixels with EBF less than 5% attained 29.73%, higher than that of the other four levels. Although pixels with EBF larger than 50% accounted for merely 1.56% of the total area, the combined proportion of pixels with EBF ranging from 15% to 50% was 53.76%, indicating that the Xiaojiang watershed was dominated by moderate and severe KRD. Therefore, Ecological Restoration Projects (ERPs) initiated by the national and local Chinese governments, aiming at alleviating severe rocky desertification and improving vegetation recovery [43], should be further promoted in this area.

3.3. An Independent Validation

An independent experiment was conducted in the Nandong watershed to validate the robustness of the developed models. Figure 6 presents the full maps of EBF throughout Nandong using the developed models with relevant spectral indices. Additional 76 validation samples were collected across Nandong by field survey to quantitatively assess the model performance when applied in other karst environments.

Table 2 presents that, among all nine indices, the KBRI had the lowest value of RMSE, MAE, and RMAPE while owning the highest value of R² (8.62%, 7.01%, 5.93%, and 0.70, respectively). For KBRI, compared with the estimation accuracy in Xiaojiang, the RMSE and MAE increased, while the indicator of RMAPE and R² decreased, indicating that the overall accuracy of the model developed with KBRI has slightly declined in Nandong, but could still meet the requirements of the application. In contrast, the estimation accuracy of most other indices experienced a substantial decline except SRI1 and SRI2. For instance, the NDRI2’s RMSE, MAE, and RMAPE separately rose from 5.90% to 11.44%, from 4.77% to 8.57%, and 14.87% to 33.03%, respectively, while its R² decreased from 0.68 to 0.45. Although SRI1 and SRI2 didn’t undergo an apparent decline in estimation accuracy, even performing better in some indicators, their overall accuracy remained comparatively low. This reflects the poor applicability of the developed models with eight other remote-sensing indices except KBRI when altering the study area. In summary, the proposed KBRI is able to effectively distinguish between exposed bedrock and bare soil (vegetation), and the developed regression model using KBRI can be well employed to estimate the EBF in other karst regions of southwest China.

The area statistics indicated (see Table 3) that the proportion of pixels with EBF ranging 30–50% reached 44.42%, much higher than that of other four levels, which verified that severe KRD was distributed widely across Nandong watershed. In comparison with Xiaojiang, the proportion of pixels with EBF less than 5% in Nandong attained 10.82%, less than half of the counterpart of Xiaojiang. Meanwhile, pixels with EBF higher than 30% accounted for 55.17% of the total area in Nandong, nearly twice as much as the counterpart of Xiaojiang. Overall, the KRD degree in Nandong was much more severe than that of the Xiaojiang watershed. Therefore, urgent countermeasures need to be taken in this area.

4. Discussion

Fractional ground cover extracted from remotely-sensed imagery has been widely used to estimate land degradation and human disturbances [44,45,46]. As the key ecological indicator of KRD monitoring and assessment, EBF quantification has been the focus of previous KRD-related research [4,7,8,47]. A commonly adopted approach is based on the regression model between vegetation indices and fractional cover. Nevertheless, vegetation indices are designed for vegetation greenness information detection [48], and thus not suitable for exposed-bedrock extraction, especially in rugged karst areas characterized by severe KRD, which is the trigger to build a more applicable spectral index for the rapid and accurate extraction of exposed bedrock.

Therefore, the main contribution of this study is the proposal of the novel approach for identifying the exposed bedrock in a highly heterogeneous karst environment by using Landsat-8 OLI imagery. In the first step, a novel index KBRI was presented aiming to enhance the difference between the exposed bedrock and other land-cover types dramatically. Then, a linear regression model relating satellite-derived information (i.e., KBRI values) to field-surveyed measurements of EBF was constructed to estimate the EBF over the whole study area. Finally, model robustness was validated through its application to another karst area with similar landform characteristics. Experiments showed comparatively good and satisfactory mapping results.

Generally speaking, the establishment of a relatively optimal index was the key for the specific land cover extraction [49]. Meanwhile, the performance and suitability of a particular index is determined by whether it is sensitive to the characteristics of interest [50]. In the karst environment, the spectral-reflectance features of exposed bedrock, vegetation, bare limestone soil, and built-up land show significant differences at the SWIR band of OLI imagery (Figure 2). Groups such as water, hydroxyl, carbonate, and sulfate are the main determinants of the spectral-absorption characteristics in the SWIR region [51]. Therefore, the SWIR is one of the appropriate bands for characterization of typical land-cover types in the karst area [4,7]. In addition, due to the obvious differences of the variation trend in spectral reflectance of different typical classes from the NIR to SWIR band, the combination of the SWIR and NIR bands was ultimately determined as the basis of the newly proposed index. In the meantime, the improvement of the mathematical operations in the applied root function further enhanced the spectral response differences between exposed bedrock and bare limestone soil and built-up land, which greatly helped to overcome the problem of differentiating the exposed bedrock and built-up land in spite of their relatively similar spectral characteristics [17].

In this study, we hypothesized that, with the increase of the fractions of exposed bedrock, the corresponding spectral-absorption feature will become more prominent [4]. Hence, the proposed index here was able to capture the spectral absorption and could thus provide quantitative information of fractions of exposed bedrock. The regression model contributed to expanding the EBF information on the sampling points scale to the regional scale. In a comparison analysis, the KBRI was found to be superior to other rock-related indices in the linear fitting relationship with in situ measurements of EBF, with relatively low RMSE, MAE, RMAPE, but high R² (Table 2). Although the overall accuracy of the model developed with NDBI was pretty close to that of KBRI, it underwent considerable decline when changing the study area. In contrast, the model developed with KBRI remained robust despite a slight decrease of estimation accuracy. Moreover, although Huang and Cai [17] proposed that the joint use of NDRI1 and NDVI was able to effectively map the karst rock in southwest China, their correlations with the EBF was comparatively small in this study, especially NDRI1. Similarly, Li and Wu [13] established NDRI2 from Landsat-8 OLI remotely sensed images to calculate the bedrock exposure in southwest China; however, it didn’t perform well in our study due to the difficulty of differentiating between exposed bedrock and bare soil. Additionally, although the two bedrock indices proposed by Xie et al. [18] (i.e., CRI1 and CRI2) were considered to be useful for estimating the fractions of exposed carbonate rock, the estimation models developed by them didn’t present high assessment accuracy, even lower when applied to other karst areas. Furthermore, CRI1 and CRI2 were not sensitive to the EBF changes, presented as concentrated and short-range value distribution (Figure 4b,c). In this study, we also built two simple rock indices (i.e., SRI1 and SRI2) based on simple-band ratios, which, however, didn’t show strong relationships with the EBF and were thus not fit for EBF retrieval. Besides these multispectral indices, Yue et al. [7] developed a hyperspectral index KRDSI to estimate the fractional coverage of exposed bedrock, with R² and RMSE as 0.53 and 0.11, respectively. Their accuracy was lower compared with our accuracy. The main reason lay in the severe weathering process of carbonate rocks in karst areas and the subsequent variability of absorption characteristics of carbonate ion [52]. In general, all of the findings of this study suggested that the KBRI could be used as an alternative and more accurate remote-sensing index to distinguish exposed bedrock from other land covers and to determine the coverage percentage of exposed bedrock.

Study Limitations

Rather than a new mechanism model, our study proposed a simple mathematical approach, which could be utilized to assess the fractions of exposed bedrock directly and automatically over large karst regions characterized by complex landscapes using Landsat-8 multispectral imagery. However, there still exist some limitations regarding our method. First, the inherent high spatial heterogeneity within the complex landscape is a major reason for the uncertainty in EBF retrieval by remote sensing [9]. Moreover, unreasonable land-use practices and ERPs implementation since the 1990s may also contribute to spatial heterogeneity [53]. Image segmentation may be one alternative way to reduce the impact of heterogeneity on bedrock coverage extraction, in which a whole image would be segmented into spatially continuous, disjunctive, and relatively homogeneous subsets [7,54]. Second, because of large terrain differences, there are a large number of backlight areas that are also affected by the reflected radiation between adjacent objects, resulting in a certain number of shaded areas on the satellite image [24]. The topographic shadows may also cause uncertainty in the EBF estimation by inducing changes in slope/aspect, illumination and shading [7,9]. In addition, southwest China is located in a subtropical monsoon climate zone, characterized by frequent rainy and cloudy weather [43,53]. Meanwhile, at present, optical remote-sensing imagery is the mainly used data source for KRD-related information extraction, such as Landsat [55], MODIS [56], and SPOT data [57]. However, due to the influence of the rainy and cloudy weather, it is difficult to obtain high-quality remote-sensing images under clear weather conditions in karst areas [19,58], which greatly limits long time-series dynamic research. Google Earth Engine (GEE), the most advanced cloud-based geospatial data-analysis platform in the world, provides free access to a large catalogue of Earth observation data, algorithms for analyzing data, and a programming interface for creating and running these algorithms using Google’s computation infrastructure [59]. In future studies, GEE will be a very useful tool to effectively mask the clouds/shadows on satellite images and derive long time series of annual/seasonal composite images. Finally, due to the fact that the high-quality sampling data are highly critical for the model establishment [60,61], steep slopes and rugged landscapes lead to very difficult access to the collection of sufficient ground truth training and validation data for modeling [20,62,63]. Therefore, how to develop a more optimized estimation model using limited training data remains a challenge in the next study.

Accurate and rapid estimation of EBF is of great significance both for quantifying past changes and for predicting future changes, which can also provide valuable basic data for KRD mapping and assessment. Future research will mainly focus on spatiotemporal dynamics of the EBF at the regional scale, and reducing the effects of high heterogeneity on the estimation accuracy of EBF, using a GEE platform.

5. Conclusions

This study presents a methodologically simple approach for estimating the EBF in typical regions of southwest China, effectively and automatically. A novel KBRI was firstly proposed, using the band combination of SWIR1 and NIR of Landsat-8 OLI remotely sensed imagery, based on the unique spectral responses of main land-cover types. The KBRI proved to be effective to distinguish exposed bedrock from other typical classes. Then, a linear-regression model was built to investigate the relationship between the proposed KBRI and field-surveyed EBF, and thus to derive the EBF results over the whole study area. The results demonstrated that the model developed with KBRI was superior to the models with other existing rock-related indices, based on the performance of the root mean square error (5.59%), mean absolute error (4.63%), root mean absolute percentage error (13.59%), and coefficient of determination (0.72). Moreover, KBRI demonstrates the potential of the robust application within a larger karst area.

Experimental results showed that pixels with EBF ranging from 15% to 50% accounted for 53.76% of the total area in the Xiaojiang watershed (542.13 km²), suggesting that the Xiaojiang watershed was dominated by moderate and severe KRD. Meanwhile, the proportion of pixels with the same EBF range (i.e., 15–50%) reached 68.60% of the whole watershed of Nandong (1109.28 km²), which indicated that the KRD degree in Nandong was severer in comparison with that of Xiaojiang. The countermeasures were thus in urgent need to alleviate severe rocky desertification and improve the environmental conditions in this area.

Consequently, results of this study could provide basic research data for KRD assessment and monitoring, and be used to enhance regional karst-ecosystem management. Despite the uncertainties of the EBF retrieval, the advantages of our proposed method using the KBRI are that it is easily implemented, robust, and needs minimal requirements for auxiliary data while still achieving comparatively better accuracy, which allows it to be applied in other developing countries around the world which have been confronted with the similar environmental problems.