Comparing Methods of Rocky Desertification 1 Monitoring at the Sub-pixel Scale in a Highly 2 Heterogeneous Karst Region 3

Karst rocky desertification (KRD) is a process where strong anthropogenic disturbances 12 and exposure of carbonate bedrock occurs in fragile karst ecosystems. The fractional cover of rocky 13 outcrops is a key indicator and mechanistic driver of KRD and can be accurately assessed using 14 remote sensing technology. Nevertheless, rugged karst terrain relief can cause shadow effects on 15 satellite imagery and combine with high heterogeneity of karst landscapes to prevent fractional 16 cover retrievals. In this study, we explored the feasibility of applying multispectral high spatial 17 resolution ALOS imagery for fractional cover extraction of rocky outcrops. We selected the 18 dimidiate pixel model (DPM), which has been applied in previous studies, and spectral mixture 19 analysis (SMA; including simple endmember spectral mixture analysis (SESMA) and multiple 20 endmember spectral mixture analysis (MESMA)) to explore the feasibility of using remote images 21 for KRD monitoring and improve accuracy for estimating fractions. Results from MESMA achieved 22 high overall accuracy (76.4%) in monitoring percentage of rocky outcrop fraction in the study area. 23 SESMA appears to underestimate percentage of rocky outcrop likely because the development of 24 KRD was driven by complex factors (soil erosion, dissolution and anthropogenic disturbances). This 25 results in spectral reflectance of rocky outcrop being variable in different settings. Predicted exposed 26 bedrock coverage using SESMA and MESMA was similar in sun-lit and shaded areas although 27 predictions from SESMA were smaller than reference data. DPM underestimated the fractional 28 cover of rocky outcrops on south-facing slopes and overestimated it in shaded areas. Furthermore, 29 SESMA and MESMA effectively reduced topographic effects. We conclude that it is better to extract 30 percentage of rocky outcrop using MESMA in the karst region of southwestern China. Remote 31 sensing is emerging as a feasible method to extract surface condition information in heterogeneous 32 and rugged landscapes. 33


Introduction
Karst landscapes, a special environment formed above carbonate bedrock, are among the most ecologically fragile regions in the world [1,2].Southwest China has one of the largest continuous karst landscapes in the world; an area of 540 000 km 2 [3].Critically, this karst region has high human population density and limited agriculture [4].Consequently, the close human / land relationship has led subsistence farming towards over-exploitation and severe environmental damage including a process known as karst rocky desertification (KRD).KRD is the progressive land degradation in which soil is partially or completely eroded [5].Effects of KRD include widespread rocky outcrops, declining land productivity, and the formation of desert-like landscapes [6].As a result, KRD, following soil erosion in the loess plateau and desertification in northwest China, is considered to be one of the most important ecological and environmental problems in China [7].Therefore, it is imperative to monitor and assess KRD conditions for the sustainable development of this region.
There are various methods available to map KRD.Previous research has largely depended on the field survey of the cover of vegetation and rocky outcrops, slope and soil distribution to map and qualitatively estimate the extent of KRD.These methods are time consuming, expensive and limited by ruggedness and the logistics of undertaking large spatial-scale research [8].Fractional ground cover metrics extracted from remotely sensed images have been widely applied to describe land degradation and human disturbance [9].When KRD occurs, the most obvious land-surface symptoms are low vegetation cover and bedrock exposure.Therefore, the fractional cover of vegetation and exposed rocks are the most useful measures of land-surface changes indicating KRD.
Satellite images were being used to map KRD and its historical changes beginning in late 1990s [10,11].Commonly used moderate resolution images included Landsat Multispectral Scanner (MSS) [12], Landsat Thematic Mapper (TM) [13,14] and Advanced Space-borne Thermal Emission and Reflection Radiometer (ASTER) [15,16].These optical satellite images have resolutions between 15 and 30 meters and are useful for extracting land cover and change conditions at regional levels.
However, large variation in karst landforms (e.g.poljes, valleys, cockpits, towers, and sinkholes) and various degrees of soil erosion create abundant small-scale niches for vegetation growth and cause soil discontinuity and vegetation fragmentation.The dissolution of carbonate rocks produces densely distributed ditches on land surfaces that intermingle with rocky outcrops, bare soil, grass, shrub and forest within a KRD region [12].Consequently, it is challenging to identify a pure, rocky spectrum on a relatively fine-scale (e.g.SPOT 10×10 m) using remotely sensed data [16].Therefore, the high degree of heterogeneity in karst landscapes results in single satellite image pixels including more than a single land object.
One potential solution is to attempt to estimate the proportion of land cover at the sub-pixel scale for heterogeneous landscapes [17].A DPM is commonly used to calculate fractional vegetation cover (FVC) at sub-pixel scales [18].This method has already been applied in the karst region for KRD monitoring [19,20].Another widely used method is spectral mixture analysis (SMA) [21].This approach supposes that reflectance for one pixel is a linear mixture of several endmembers and that each endmember is a unique land cover type with a specific feature.The aim of SMA is to decompose mixed spectra and calculate proportions of each land cover type in a single pixel.The spectral unmixing model has been widely applied to plant species identification [22,23], fire severity [24,25] and urban remote sensing [26,27] with some success.However, there have been few SMA-related studies in the karst region of southwest China.
A critical step for SMA is endmember selection.Unlike the impervious surfaces found in urban areas, bare rocks often intermingle with vegetation and soil.In this region, the degree of rock exposure has changed due to variation in natural conditions and human disturbance [28].
Furthermore, as high albedo endmembers, rocks, cement road surfaces, building roofs and limestone soils can cause spectral confusion in an image because they have similarly high reflectance.All these directly affect the endmember selection in KRD.In rural karst regions, roads and building roofs are between 4 and 10 meters wide.High spatial resolution imagery (e.g., SPOT-5 and ALOS imagery) likely provides more information to separate roads and roofs from KRD in this region.Compared with medium resolution images that are commonly used in the karst region, high resolution imagery could provide greater detail and capture more spatial variation that would be helpful to explore the mechanisms of KRD dynamics and endmember selection.Thus, high spatial resolution imagery is promising for fine-scale karst land applications.
In addition to high lateral heterogeneity, the relative high elevation contrast in the karst region causes significant shadow effects in remote images.The topographic effect in images limits interpretation because weak reflectance from shaded areas commonly complicates information extraction [29].Use of NDVI is able to mitigate variable topographic effects on the spectral properties of land surfaces and so is wildly used in the DPM for FVC estimations in karst regions [30].In contrast, SMA treats shade as a separate image component and eliminates the shadow fraction by area redistribution [31].However, there is no related research that proves its applicability for topographic correction in karst regions.
Our study attempted to take advantage of high resolution optical multispectral ALOS images to extract fractional cover of rocky outcrops at a sub-pixel level in the karst region of southwest China.
We applied and compared the efficacy of DPM and two SMAs methods through 1) accuracy evaluation of extraction of the rocky outcrop fractions, 2) assessment of topographic effects of fractional images from each method, and 3) analyzing the KRD distribution.Results from this study serve as a technical reference for the application of optical remote sensed imagery in heterogeneous and rugged terrain regions.

Data and Preprocessing
ALOS imagery was acquired on June 4th, 2009.These data contained four multispectral bands (blue, green, red and near infrared) with a 10 m spatial resolution and a panchromatic band with 2.5 m spatial resolution.We chose multispectral ALOS imagery for KRD information extraction at subpixel scales for two reasons.Rocky outcrops and a variety of vegetation commonly intermingle in an ALOS 10×10 m pixel.Second, pan-sharpen ALOS imagery has higher spatial resolution that can be used for the quantitative assessment of linear spectral unmixing.The solar elevation and azimuth of this image is 73.3°and 93.8°, respectively.Although relatively large solar elevation generally reduces shadow effects on images [33], this image still had a high percentage (about 25% of total area in this study) of shadow (Figure 1A).The image had level L1G processing (systematic correction).
Atmospheric correction of multispectral ALOS imagery is conducted according to the MODTRAN5 model.We georeferenced all images and projected them to a Universal Transverse Mercator (UTM) map projection.The original multispectral and panchromatic images were fused using the Gram-Schmidt procedure in the ENVI software package [34] to produce four-band, pan-sharpened multispectral ALOS images.

DPM
We used DPM to calculate fractional vegetation cover (FVC) [35].The percentage of rocky outcrops equals one minus the FVC value as soil background was covered by crops.The Normalized Difference Vegetation Index (NDVI) was applied: (1) where and are the surface reflectance of near infrared and red bands respectively.It has been demonstrated that FVC follows a linear relationship with the NDVI: where NDVIs and NDVIv are representative values of NDVI for bare rock pixels (where FDPM = 0) and a 100% vegetated pixel (FDPM=1) respectively.

Spectral Unmixing Models
Simple endmember liner spectral mixture analysis (SESMA) and multi-endmember spectral mixture analysis (MESMA) were applied to extract bedrock cover.SESMA extracts only one optimal endmember for each independent class.MESMA expands the values and focuses on selecting one or several representative endmembers within each land cover type.Both SESMA and MESMA assume that the measured spectrum is a linear combination of the spectra from all components within the pixel [36].The mathematical model for this assumption is expressed as: ( where i is the number of spectral bands used; k values are 1, . . ., n (number of endmembers); Ri is the spectral reflectance of band i in a pixel that contains one or more endmembers; fk is the proportion of endmember k within the pixel; Rik is the spectral reflectance of endmember k within the pixel on band i, and is the error for band i.To solve fk, the following conditions must be satisfied: (1) selected endmembers should be independent, (2) the number of endmembers should be less than or equal to the number of spectral bands used, and (3) selected spectral bands should not be highly correlated.
Theoretically, therefore, there were a maximum of three endmembers (land cover types) in four-band ALOS data.A common approach for obtaining fk is to use a least-squares estimation process by minimizing the residual error with the sum of fk of all optical endmembers equal to one.The spectral mixture analysis consists of three primary steps: 1) selection of candidate endmembers to build a spectral library, 2) optimal endmember selection, and 3) decomposing the mixed pixels to extract fractional images.

Spectral Library Building
We applied the vegetation-high albedo-shadow model that was found to be most suitable for non-urban areas to extract remotely sensed information [37].Vegetation endmembers included grass, shrub, forest and crops in the karst region.Ground features with high albedo include rocky outcrops, road surfaces, and rooftops.Low albedo features include shadow and water.As only three endmembers could be selected for model input, road surfaces and rooftops were masked according to their reflectance in the red band.Their values are generally higher than that of rocky outcrops based on field surveys.Open water areas were also masked using the ratio of green and near infrared bands [38].Consequently, other features, including vegetation, rocky outcrops and shadow, were kept for SMA modeling.We used a spectral scatter plot to select candidate endmembers.This scatter plot triangulated with the three vertices representing vegetation, rocky outcrop and shadow.Finally, 291, 185 and 159 candidate endmembers for three land types (vegetation, rocky outcrop and shadow) were identified for further processing.

Optimal Endmember Selection
Determining high quality image endmembers has been identified as a critical stage of spectral mixture modeling [39].We initially selected optimal endmembers from images according to field survey sites to compare with two other methods.The SESMA method is a fully-constrained linear spectral unmixing approach based on a vegetation-high albedo-low albedo model [27].The spectral scatter plots were produced as a triangle (here defined as △OAB).The best three endmembers could be identified by the largest area (S) among all triangles formed by any three pixels: Where and denote the endmember vectors of OA and OB, respectively.
For the MESMA method, we used an average RMSE (root mean squared error) (EAR) approach to select the most appropriate endmembers.The endmembers are selected by producing the lowest RMSE within a class [31].EAR was calculated according to Eq. (4) using: ( where i is an endmember, j is the modeled spectrum, N is the number of endmembers, and n is the number of modeled spectral.The "−1" corrects for the zero-error resulting from an endmember model itself.Here, eight endmembers of vegetation, rocky outcrop and shadow were selected according to the lowest RMSE, respectively.

Fraction Extraction
Fractional maps of endmember land cover components were generated by optimal endmember models using least-squares solutions.To reduce shadow effects, we performed shade normalization of the fraction images by dividing the cover of each endmember by the total percent cover of all nonshade endmembers (1 minus shade fraction) in each pixel.This suppresses the shade fraction so that we extract more information from the other two fractions (vegetation and rock) [40].

Accuracy Assessment
Accuracy assessments for rocky outcrop coverage were conducted using error matrices.A Kappa coefficient was used to measure the accuracy of the predicted rocky outcrop coverage [41].
Field validation sites were sampled from 2009 to 2011 in Du'an county.Each site was 30 × 30 m.
Transect sampling [43] and visual observations were combined to estimate the fractional cover of vegetation and rocky outcrop at each site.Each site took 5-8 pictures as recorded data.The transect sample sites were mainly located near roads or paths as high mountain areas are difficult to access.
Most validation sites for accuracy assessments were determined by visual estimates of rocky outcrop cover from a distance during field visits.Pan-sharpened multispectral ALOS (2.5 m) images were used as auxiliary data to identify land features for some sites where the rocky cover is hard to estimate.
Images from Google Earth were also used to identify features in shaded areas for validation.As indicated above, we built three spectral libraries using ALOS multispectral images based on separate methods.Eight representative spectral curves from each land type were collected from ALOS images according to field sample sites and MESMA (Figure 2a, b).Spectral curves in vegetation and shadow classes had similar shapes.Although the spectral data extracted based on field sample sites were deemed pure endmembers, spectral curves of rocky outcrops had large variation.One of the reasons for this is that differences in aspect and slope at sample sites caused variation in image spectral signature.In addition, surface characteristic changes (e.g., color) of rocky outcrops affect spectral reflectance.Therefore, there is uncertainty in how rocky outcrop characteristics are represented using spectral data.Based on MESMA, eight representative spectral curves of each land type were collected according to the lowest RMSE within each class.Spectral curves within each class are similar, though the curves of rocky outcrops had some variation (Figure 2b).This indicates that the spectral signature of rocky outcrops might also vary.
One representative spectral curve from each class was collected from SESMA (Figure 2c) and spectral curves of three land types presented separately.The spectral curve that has high reflectance in the 520-600 nm and 610-690 nm bands represents the characteristics of rocky outcrops.Comparing these three methods, the spectra of representative endmembers determined at field sample sites was not as consistent and was difficult to use to represent the rocky outcrop features.Therefore, the spectral selections from MESMA and SESMA were used to extract rocky outcrop coverage.Error matrixes for sunlit areas were used to explore estimate error in the percentage of rocky outcrops.In sunlit areas, MESMA successfully predicted the cover of bedrock outcrops in each level.
Producer's accuracies (PA) were between 0.79 and 0.93.SESMA achieved a higher accuracy in the areas with lower bedrock outcrop cover between 10% and 50%.However, accuracies were lower when cover was large than 50% (the producer's accuracies were less than 0.19).Many areas with high bedrock outcrop cover had been mis-classified as low bedrock outcrop area.This means that the coverage of bedrock outcrops was underestimated based on SESMA in sunlit areas.DPM obtained high accuracies in low and high bedrock outcrop cover areas.For the bedrock outcrop cover between 10% and 70%, the producer's accuracies were less than 60%, mainly because of the underestimation of the coverage of bedrock outcrops.
The accuracy of percentage of rocky outcrops was lower in the shaded areas than in sunlit areas.
Using MESMA, the highest producer's accuracy was 1.00 and the lowest was 0.34.A confusion matrix showed that the rocky outcrop cover was slightly overestimated for each cover class.Although the results from SESMA reached high producer's accuracies in areas where the rocky outcrop cover was between 10% and 50%, the predicted coverage was underestimated in totally shaded areas.
Alternatively, DPM clearly overestimated rocky outcrop cover in shaded areas suggesting that it is greatly affected by shadow.The mean cover of rocky outcrops estimated from MESMA had a similar value in sunlit and shaded areas in each rocky outcrop category (Figure 3).Although the cover of the shaded area is smaller than that in the sunlit area, the largest difference in predicted values in these two areas was no greater than 6.4%.The predicted value from MESMA also matched the reference data well.Each prediction in sunlit and shaded areas fell into the range of reference data.The estimated percentages of bedrock outcrop estimated by SESMA also are similar in sunlit and shaded areas.However, the predicted cover values are underestimated, especially in areas with large areas of rocky outcrop.For example, where cover ranged from 51% to 70% the prediction was just less than 46.1%.Using DPM, we found large differences in predicted values for sunlit and shaded areas in each rocky outcrop category.The predicted cover in sunlit areas was between 5% and 15% smaller than reference values in areas with rocky outcrop cover less than 90%.However, we found the opposite results in shaded areas where predicted values were 15%-54% larger than reference data.Therefore, the percent of rocky outcrop bedrock was underestimated in sunlit areas and overestimated in shaded areas as predicted by DPM.Due to predicted differences in the three methods, the rocky outcrop cover of 40.8%, 28.4% and 36.6% in the whole study area were obtained from MESMA, SESMA and DPM respectively.
The standard error (SE) of predicted cover was distributed between 1% and 20% (Figure 3).The SE from MESMA and SESMA had similar trends in each class of predicted cover.They had low SE values (5%-12%) in the most rock covered areas except for areas with the highest rocky cover.Low values indicated that the predicted cover of rocky outcrop is more accurate.In total, the SE in the sunlit areas was lower than that in shaded areas when we used these two methods.Consequently, predicted percentage was more accurate in the sunlit areas than in shaded areas.DPM had higher SE values compared with these two SMA methods suggesting that the value of predicted cover was distributed discretely compared with reference data.Shadow effects made similar land features look different and varied spectra in sunlit and shaded areas on satellite images (Figure 4a).Compared with the original satellite images, fractions from MESMA and SESMA had few spectral changes and therefore had similar grey values for similar land types in sunlit areas (white border) and shaded areas (yellow border) (Figure 4b, c).This indicated that the rocky fractions predicted by MESMA and SESMA might reduce shadow effects.Compared with that in Figures 4b and 4c, the rocky fraction from DPM was also affected by shadow illustrated by a large difference in white and yellow bounded areas (Figure 4d).
We used the cosine of incident angle to represent the shadow effect of satellite images.The lower the correlation coefficient between the rocky fraction and the cosine of incident angle, the smaller the impact of the shadow [29].This result showed that the percentage of rocky fractions from MESMA and SESMA had small correlation coefficients with the cosine of incident angle (0.02 and 0.01, respectively; P < 0.01).Contrast this with the rocky fraction of DPM and cosine of incident angle that had a relatively high correlation coefficient (0.33; P < 0.01).A higher regression coefficient means the rocky fraction from DPM was affected by terrain relief.This indicates that shade normalization by MESMA and SESMA is better at reducing topographic effects to predict the percent of rocky outcrops.With the exception of non-KRD and potential KRD, the summed areas of light KRD, moderate KRD and severe KRD occupied respectively 62.7%, 39.5% and 56.5% of the total area based on MESMA, SESMA and DPM (Figure 5).Based on MESMA, the non-KRD, potential KRD, light KRD, moderate KRD and severe KRD occupied respectively 6.8%, 30.5%, 34.0%, 20.6% and 8.1% of the study area.Potential, light and moderate KRD dominated the region.The main predicted KRD type from SESMA is potential KRD that occupied 53.8% of the study area.As the study area includes typical KRD landscapes where vegetation degradation, severe soil erosion and exposed rock is common [13], KRD areas appeared to be underestimated by SESMA.Although the total area of KRD from DPM appeared reasonable, its distribution was affected by terrain relief.This was because the underestimation of rocky outcrop cover in sunlit areas corresponded to the distribution of potential and light KRD.In contrast, the overestimation of rocky coverage in shaded areas largely correlated to the distribution of moderate and severe KRD areas.

KRD Extraction
Monitoring rocky desertification is an important task for environmental management in southwest China.Previous research on remote sensing of KRD estimated the percentage of rocky outcrop cover using DPM [11].According to our results, this method underestimates fractional rocky cover by 5% to 15% in sunlit areas.This could be caused by the nonlinearity of NDVI; NDVI tends to be sensitive to vegetation cover where it is sparse [45,46].Therefore, the fractional vegetation cover might be overestimated in sunlit areas and the coverage of rocky outcrop underestimated accordingly.
The high heterogeneity of the karst landscape posed a within-pixel mixing problem for remote sensing information extraction and SMA is considered to be a feasible solution [25].Although there are no studies that use MESMA-based fraction images to estimate KRD levels, SMA has proven to be efficient in detecting a signal from impervious surfaces in urban areas [27].However, unlike the regularity of impervious urban surfaces, the features of rocky outcrops have considerable surface variability.Because of differences in human disturbance intensity, degrees of weathering and erosion, rocky outcrop surface colors may also vary.Field investigations have shown that bedrock outcrops often had white-gray surface color after the soil was ploughed, mined or the vegetation clipped with soil runoff.Especially in mined areas, rocky outcrops had high spectral reflectance on cut surfaces (Figure 2a; a spectral curve of bedrock with high values).When natural forces, like weathering or dissolution, dominated, surface bedrock outcrops appeared gray.Color differences of rocky outcrop surfaces might cause variation in spectral reflectance (Figure 2b).
Using SESMA, we selected one optimal endmember for three land cover types from a vertex of spectral scatter plots [27].This approach only allowed one spectrum as a pure pixel for each endmember.In our study, the spectra of rocky outcrops from SESMA had higher reflectance compared with those from MESMA (Figure 2).If this one optimal endmember represents rocky outcrops well, the fractional cover would be estimated more accurately.Otherwise, the spectrum might represent part of a bedrock feature (like a cut rock surface) and the fractional cover would be easily underestimated.According to the error matrix of accuracy assessment, results from SESMA underestimated between 1% and 27% of the fractional cover of rocky outcrops.Many high rocky outcrop areas were classified into low or moderate rocky covered area.This may have been caused by SESMA's optimal spectral curve not representing the general features of rocky outcrops in karst areas.For that reason, an optimal endmember collection based on the SESMA model is difficult to incorporate the spectral variation of feature classes in the karst region.The underestimation of SESMA resulted in us predicting that the rocky outcrop coverage of the whole study area was 28.4%.This didn't match the observation that the study area is a typical KRD area where the major KRD land types are light and moderate KRD [13,16,44].MESMA extends SESMA by allowing the number and type of endmembers to vary on a perpixel basis [22].It overcomes SESMA's limitations by requiring a model to meet minimum fit, fraction and residual constraints while testing multiple models for each image pixel.The optimal endmembers from MESMA are selected by producing the lowest RMSE within a class [31].MESMA allows a certain number of optimal endmembers to be selected from images.When one land cover has variable spectral reflectance, it is better to have many representative spectral curves to represent its features.Indeed, both natural factors and human activities have contributed to the occurrence of KRD, adding to image variability [11,16].This diversity of factors has resulted in some differences in the exposed state of the bedrock.The results of accuracy assessment showed that much of rocky outcrop cover could be accurately predicted.Therefore, multiple endmembers are more reliable than single endmembers to detect and monitor karst exposed rock.
A possible shortcoming in our optimal endmember selection is the number of endmembers.
Although our study masked built-up areas and water, and selected remotely sensed images in the crop-growing season so that bare soil was covered by crops, images still had bare soil.KRD is, by definition, the landscape of exposed bedrock after the soil is lost [5].Rocky outcrops appear as ragged on the images and are accompanied by countless rock fissures after weathering and soil erosion.
There is often some soil, deciduous plants or bryophytes in those rock cracks [16].Because of the restricted number of bands of ALOS image, three types of endmembers were suitable for selecting a model fraction image.Field measurements of the spectra on the karst land surface show that nonphotosynthetic vegetation, bare soil and rocky outcrops had similar spectral responses at 350-850 nm.However, differences were mainly focused at 1350-2350 nm [47,48].Therefore, bare soil might be confused with exposed rock in SMA modeling in our study.In future studies, imagery with more spectral bands (such as WorldView-3 that has 8 bands) should be used to enable both the selection of more endmembers and improve the accuracy of remote sensing information extraction in highly heterogeneous regions.

Topographic Effects
In karst regions, rugged terrain is stimulated by the strong development of carbonate rock.Steep elevations are common but very localized.Rugged terrain not only affects vegetation growth but also affects the extraction of vegetation information [30].This is not unusual.For medium and high spatial resolution satellite imagery, the spectral reflectance of land cover is often affected by terrain relief.In our study this caused variation in DN values in sunlit and shaded areas.Previous studies have shown that a topographic correction model does not effectively improve the classification accuracy of ground objects [11,33,49].Furthermore, there are few studies that compare spectral differences of the same land objects between sunlit and shaded areas in a karst region after topographic correction.The results of our study have shown that rocky outcrop coverage predicted by SESMA and MESMA was similar in sunlit and shaded areas, largely because the shadow fraction, as an independent component, was extracted by least squares in the SMA model and removed through shade normalization [40].In fact, shadows from mountains, tall trees, and even protruding stones, commonly occur in karst regions.However, previous studies paid less attention to shadow and its impacts on remote sensing applications in karst areas [29].One of the advantages of SMA model is that it accounts for shade, demonstrating that SMA provides a new perspective of topographic correction for remote sensing applications in a mountainous region.
DPM, as a simple and effective method for managing topographic effects, is widely used to monitor vegetation cover [28,50].However, DPM predicted vegetation cover was very different between sunlit and shaded areas in our study.The reason appears to be that, when an aspect changed from sun to shade, the reflectance of the near infrared band decreased more rapidly than the red band due to the shadow effect (Figure 2).Based on the DPM formula, the predicted value of rocky outcrop cover becomes large in shaded areas.The difference between sunlit and shaded areas caused by topographic effects influenced the accurate interpretation of land surface information from KRD and further affected KRD mapping.Therefore, topographic effects based on DPM interpretation could not be ignored in the optical remote sensing monitoring in this mountainous region.

Conclusions
Aimed at overcoming obstacles to monitoring heterogeneous, rugged terrain in rocky deserts using remote sensing, our study applied high spatial resolution ALOS images and compared DPM, SESMA and MESMA to extract key indicators of KRD at sub-pixel scale.The best accuracy assessment results were acquired using MESMA with overall accuracy in the sunlit and shaded areas of 83.7% and 60.4%, respectively.The SESMA approach achieved lower accuracy because it underestimated 1% to 257% of the fractional cover of rocky outcrops.However, the mean coverage of the same objects was similar in sunlit and shaded areas although the accuracy in the shaded area was lower based on these two SMA methods.Correlation analysis and spatial statistics demonstrated that SMA methods could reduce shadow effects on fractional cover extraction in karst regions.Shadows that came from mountains, trees or raised rocks, were widespread and shadow effects were one of the critical factors affecting classification accuracy of land use and land cover in remote sensing applications in the karst region of southwestern China.One of the advantages of the SMA model is that it accounts for the shadow fraction and can remove its influence.
We conclude that SMA provides a new perspective on topographic correction for remote sensing applications in mountainous regions.KRD distribution can be mapped accurately with the fractional estimation by MESMA.The fractional cover of rocky outcrops may be underestimated in the sunlit areas and overestimated in shaded areas with DPM methods.Therefore, the topographic impacts on NDVI should not be ignored in karst vegetation monitoring.

Figure 1 .
Figure 1.Study area and characteristics of KRD.Panel A is a pan-sharpened ALOS image (RGB: NIR, red, and green).This image is located within study area (B).The red, yellow, green and white squares on panels A and B represent reference sites with severe, moderate, light and potential KRD respectively.

Figure 2 .
Figure 2. Spectral curves of rocky outcrop, vegetation and shadow in the study area.a) collected from ALOS image-based field survey sample sites; b) obtained based on MESMA; and c) from SESMA.

Figure 3 .
Figure 3. Predicted percent of rocky outcrop in sunlit and shaded areas.

Figure 4 .
Figure 4. Contrast of rocky outcrop cover in sunlit and shaded areas.Panels a, b, c and d are ALOS sharpened image (RGB: NIR, red, and green bands) fraction rock cover predicted using MESMA, SESMA and DPM, respectively.Circle, square, triangle, diamond and ellipse symbols represent the rock coverage less than 10%, 20%, 30-40%, 55-70% and 75%, respectively.Symbols in white are in sunlit areas and those in yellow are in shaded areas.Bounded areas have identical land cover.

Figure 5 :
Figure 5: Comparison of KRD mapping using MESMA, SESMA and DPM based on ALOS imagery.

Table 3 :
Percent area in each KRD level based on different methods.