An advanced coverage estimation method to quantify biological soil crust coverage using Sentinel-2 imagery in desert and sandy land of China

ABSTRACT

Monitoring the distribution and area change of biological soil crusts (BSCs) can enhance our understanding of the interactions between nonvascular plants and the environment in drylands. However, using only pixel-based binary classification methods results in large-area estimation errors at large scales. The lack of available calculation methods for directly measuring BSC coverage using multispectral satellite images makes it challenging to obtain BSC area data for further studies at large scales. To address these issues, this study developed feature space conceptual models for desert and sandy land based on the characteristics of BSC in drylands. The desert feature space comprised the normalized difference vegetation index (NDVI) combined with the brightness index (BI), encompassing moss, lichen, and non-BSC. The sandy land feature space relied on the biological soil crust index (BSCI) and the NDVI, including vegetation, mixed BSCs and sandy soil. Using Sentinel-2 satellite imagery and a spectral unmixing model, the abundance of BSCs was quantified in four BSC growth areas located in the Gurbantunggut Desert and Mu Us Sandy Land of China. Validation of the method indicated that the root mean square error (RMSE) of the BSC coverage estimation results was 10% and 8% in desert and sandy land, respectively (estimation accuracies of 79% and 81%, respectively). This demonstrated that the proposed method can effectively estimate BSC coverage at a subpixel scale. The resulting BSC coverage data can provide the possibility to evaluate the functions of regional ecosystems.

KEYWORDS:

• Biological soil crust
• BSC coverage
• feature space
• spectral unmixing model

1 Introduction

Biological soil crusts (BSCs) are among the most distinctive and pioneering organisms in drylands and are often used as an important indicator of ecosystem stability and the restoration of degraded ecosystems (Zhou et al. 2019). BSCs can be broadly categorized into three types, namely, algal crusts, lichens and mosses. They begin to slowly develop when algae colonize the soil surface (Chen et al. 2009). BSCs possess remarkable adaptability, enabling them to survive in arid environments and withstand extreme conditions such as high temperatures, alkaline pH levels, and high salinity (Chongfeng et al. 2017). Their presence plays a significant role in accelerating species succession and influencing soil properties, particularly in areas with limited water and nutrients (Chamizo et al. 2012). Given their wide distribution, these organisms enhance fertility in terms of nitrogen and carbon, reduce soil erosion, and fix nitrogen and carbon from the atmosphere, allowing them to play critical roles in dryland ecosystems, especially in the face of climate change (Pietrasiak et al. 2013). Consequently, the development of effective BSC monitoring and management strategies is essential.

Remote sensing technology is a necessary tool used to map the distribution of BSCs. Different indices have been developed based on hyperspectral and multispectral data according to the subtle spectral characteristics of BSCs (Beaugendre et al. 2017; Ben-Dor et al. 2003; Chen et al. 2019; Escadafal 1989; Goldshleger et al. 2010; Karnieli 1997; Mathieu et al. 1998; Rodríguez-Caballero, Escribano, and Cantón 2014; Ustin et al. 2009; Weber et al. 2008). High-spatial-resolution hyperspectral sensors are very sensitive to BSC changes and can effectively detect the spectral characteristics of chlorophyll and other components in BSCs. Weber et al. (2008) and Chamizo, Cantón, et al. (2012) developed the continuum removal crust identification algorithm and the crust development index, respectively, based on the absorption characteristics of BSCs and successfully applied these indices to analyze hyperspectral images. However, hyperspectral satellite data have limitations, such as acquisition difficulties and high processing costs, and are thus not applicable for large-scale applications, especially those performed at such large spatial scales as the arid and semiarid areas of China. In contrast, multispectral satellites have a very short revisit period and a large geographical cover range and are thus a feasible choice for computing in large-scale mapping. The use of such data can expand the research scope to a global scale. However, only a few studies have successfully used multispectral sensors to recognize BSCs (Chen et al. 2005; Karnieli 1997; Panigada et al. 2019). Karnieli (1997) and Chen et al. (2005) used the differences between BSC and other surface components to develop the crust index and biological soil crust index (BSCI), respectively, based on multispectral sensors. When mapping BSCs in large-scale areas using multispectral sensor imagery, the reflected signals observed by the sensors contain a mixture of BSCs and non-BSCs due to the limited resolution of sensors, such as a 10 m Sentinel-2 image and a 30 m Landsat-8 image. Binary classification assumes that each pixel is either entirely BSC or non-BSC, which may not be accurate for areas with mixed vegetation cover or partial BSC coverage. Therefore, if the pixel size is larger than the research object in remote sensing images, it is necessary to consider the actual coverage of the research object within each pixel (Rodríguez-Caballero et al. 2017). On the other hand, if the statistical BSC area is derived by pixel counting, great errors will arise in large-scale quantitative monitoring. These errors can affect the estimation of ecological parameters related to BSC area, such as the annual carbon sequestration capacity of drylands, surface erosion, runoff and the calculation of the carbon and nitrogen contents. Ignoring such errors can hinder the accurate assessment of the ecological effects of BSCs in subsequent studies Yong-Sheng et al. (2023); Bullard, Strong, and Aubault (2022); Maestre et al. (2013). Considering the crucial role of BSCs in influencing desert and sandy land ecosystems, it becomes necessary to accurately map the spatial abundance distribution of BSCs at a large scale and calculate the BSC area to enable a quantitative assessment of their ecological impact. However, the index-based detection methods mentioned above are not sufficient for calculating BSC coverage, mainly because image pixels containing surface components at different proportions may have the same spectral index value (Wang et al. 2022).

Machine learning and deep learning, both of which are useful methods for mining information from data, have been applied to estimate BSC coverage. To date, some studies have realized BSC coverage mapping with the help of the random forest model and neural network model (Chen et al. 2019; Liming et al. 2021). These studies have shown that advanced computing algorithms are effective when applied to quantify BSC coverage. However, the calculation of these methods is complex and requires considerable a priori knowledge. Additionally, such methods may be suitable for areas with a single BSC type because these methods do not fully consider spectral characteristic differences among various BSCs. Spectral mixture analysis (SMA) has been thought to be an effective method for calculating the component proportions in mixed pixels (Powell et al. 2007; Wenliang 2020). Its principle is that the spectrum of a pixel is a linear combination of the reflectance of all surface components within the pixel range (Pacheco and McNairn 2010). A spectral unmixing model was first developed for hyperspectral analysis, in which the information of mixed pixels can be extracted from many spectral bands (Tian et al. 2021). Additionally, different spectral mixture methods have been developed based on SMA, such as normalized spectral mixture analysis and multiple endmember spectral mixture analysis (MESMA) (Fernández-Manso, Quintano, and Roberts 2012; Quintano, Fernández-Manso, and Roberts 2013).

In drylands, BSCs are predominantly found in deserts and sandy lands. When conducting large-scale mapping of BSC coverage using remote sensing technology, it is crucial to determine the appropriate BSC estimation scheme according to the dominant BSC types in different land cover areas. This differentiation is necessary because the factors influencing BSCs vary in different land types, leading to variations in growth, succession, and water content, which consequently result in distinct spectral characteristics (Wang et al. 2023). Deserts typically exhibit low vegetation coverage and are primarily in arid regions. The development of BSCs in deserts is primarily influenced by natural factors. Sandy lands are located mainly in semihumid and semiarid areas. Due to population growth, sandy lands in China have experienced significant damage from activities such as excessive grazing and extensive exploitation. In response, the government has implemented ecosystem restoration projects such as the Grain for Green and Three-North Shelterbelt Program (Gao et al. 2022; Lin et al. 2022). These projects have led to a combination of natural and anthropogenic influences on BSCs in sandy lands. As a result, it becomes necessary to conduct studies on these areas separately. On the other hand, in deserts where only a few surface component types exist, differences among BSC types are evident, and the ecological effects of both lichens and mosses are significant. Unfortunately, no study has been performed to implement the classification of different types of BSCs without being based on additional classifiers. Additionally, it is not possible to distinguish among areas of different BSC types within a pixel in multispectral images. Classifying the images first can largely improve the purity of the mixed pixels and reduce the uncertainty of the mixed spectrum. This can be achieved with high-resolution hyperspectral sensors. Rodríguez-Caballero, Escribano, and Cantón (2014) used hyperspectral data to quantify the proportion of algal crust and lichen after determining different types of BSCs using a support vector machine. However, only applying this idea to multispectral sensors can help to accurately calculate the BSC coverage at a large scale. We tested a new BSC coverage estimation model in the Gurbantunggut Desert and Mu Us Sandy Land, where BSCs are widely distributed.

The main objectives of this study are as follows: (1) exploring the spectral characteristics of the main surface component types in desert and sandy land; (2) proposing feature spaces for two land cover types based on separation indices; (3) constructing a BSC coverage estimation model based on the spectral unmixing model; (4) mapping the BSC coverage of the study areas using Sentinel-2 images; and (5) assessing the robustness of the proposed method over areas with different land cover types.

2 Materials and methods

2.1 Study area

2.1.1 Desert

The Gurbantunggut Desert in Xinjiang is the second-largest desert in China (Zhou et al. 2020). The annual precipitation range in this desert is 70–150 mm, and the annual average temperature is 6–10°C (Shuai et al. 2021). The southern and southwestern parts of the desert are close to the oasis area of the economic belt on the northern slope of Tianshan Mountain, which is an active area of land development in Xinjiang (Figure 1a). Reclamation is most obvious in the southwestern desert region, greening construction has been carried out in the desert fringe areas, and a desert forest and sand-control forest have been rebuilt and newly built. In winter and spring, snow and ice melt, replenishing the desert with plenty of water; thus, short-lived plants have developed. Lichens and mosses are common organisms in the Gurbantunggut Desert, and their presence is closely related to the characteristics of the desert environment. Due to the blocking effect of the Himalayas, the moisture flow from the Indian Ocean cannot reach the region, and the desert has been in drought. Despite the water and nutrient deficiencies, lichens and mosses can survive and reproduce in this harsh desert environment due to their specialized water absorption structures and photosynthetic abilities (Zaady, Karnieli, and Shachak 2007). During the non-plant growth season, the main surface component types are dry plant materials, sands, mosses and lichens (Table 1). The northern and central-eastern areas of the desert were selected as the study area. Area A is the area where mosses are the dominant surface component type, and Area B is the area where lichens are developed (Figure 1a). Because these BSCs are in deserts and have developed over a long period, there are significant differences in the physical states of the mosses and lichens, as shown in the photos in Figure 1a.

Figure 1. Study areas and sample points, (a) areas A and B in the Gurbantunggut desert, Xinjiang and (b) areas C and D in the Mu Us Sandy Land. The photos of each study area depict the biological soil crusts present in those regions. Areas A and B represent desert areas with predominantly mosses and lichens as the dominant biological soil crusts. Areas C and D are sandy lands characterized by the presence of mixed biological soil crusts. The satellite images and photos were derived from Sentinel-2 and field surveys. China cover data were from the Chinese Academy of Sciences (Bingfang et al. 2017). (c) Overview of drylands in China.

Table 1. Regional overview of the four study areas.

Region	Surface component type	Scope (Lon, Lat)	Climatic zone
Area A	Lichen, moss, sand, and plant litter	Lon: [88.33, 88.57] Lat: [44.64,44.72]	Arid area
Area B	Lichen, moss, sand, and plant litter	Lon: [88.72, 88.86] Lat: [44.44,44.49]	Arid area
Area C	MBSCs, vegetation, and sandy soil	Lon: [106.12, 106.31] Lat: [37.57, 37.69]	Semi-arid area
Area D	MBSCs, vegetation, and sandy soil	Lon: [107.13, 107.35] Lat: [37.54, 37.67]	Semi-arid area

2.1.2 Sandy land

The sandy study areas are in the Mu Us Sandy Land, including Litong District, Yanchi County and Dingbian County (Figure 1c). These regions are typical sandy land management areas influenced by policy and have a semiarid continental climate with precipitation concentrated in the summer season. BSCs develop in the study areas due to the presence of grazing bans and effective wind and sand management practices. Field photographs depict BSCs in an advanced stage of succession, as evidenced by Figure 1b. There is a significant ecological niche overlap between mosses and lichens. In this study, we categorized this stage of BSC as mixed BSCs (MBSCs). For representative area selection, we identified two prevalent types commonly observed in sandy land: shrub-dominated areas (Area C in Figure 1b) and grassland-dominated areas (Area D in Figure 1b). Due to the accumulation of vegetation humus and calcium, the main types of soil in this region are chestnut calcareous soil, brown calcareous soil, and black loamy soil. These soils are mixed with sand, forming a soil with a higher proportion of sand particles. We referred to this type of soil as sandy soil (Table 1).

2.2 Data

2.2.1 Hyperspectral data

A FieldSpec 4 HR NG spectrometer (Analytical Spectral Devices Inc.) was used to measure the spectra of BSCs with different coverages and other types of surface components from September 26-15 October 2021. Its spectrum measurement range is 350–2500 nm and the field of view is 25°. We marked the measuring area of the spectrometer with a ring with a diameter of 37 cm, and the spectrometer was 83 cm above the ground. The spectral data of BSC sample sites with different coverages were measured under dry and clear-sky conditions, including 42 moss sample sites in the Gurbantunggut Desert and 23 MBSCs sample sites in the Mu Us Sandy Land. Five measurements were made at each site, and a white board was used to correct the spectrometer before each measurement. The average of the five measurements was then taken as the final value of the reflectance at that site. After each spectral measurement, the area within the spectral measurement range was photographed using a digital camera. Based on eCognition and ENVI 5.3, the supervised classification method was used to extract BSC areas from the photographs and calculate their coverage. In addition, spectral samples of plant litter, vegetation, sandy soil, sand and lichen were obtained.

2.2.2 Satellite data

The Sentinel-2 satellite observation system is a land monitoring constellation composed of two satellites. It is a group of remote sensing earth observation satellites in the Copernicus project series of the European Space Agency (ESA). The Sentinel-2 multispectral instrument (MSI) has 13 spectral bands, including the visible (VIS), red-edge (RE), near-infrared (NIR) and shortwave infrared (SWIR) bands, and has the advantages of a high spatial resolution and 5-day revisit period (Bannari et al. 2018). The VIS and NIR bands both have a spatial resolution of 10 m, and the RE and SWIR bands have a spatial resolution of 20 m. In the Gurbantunggut Desert, four cloud-free Sentinel-2 surface reflectance images from October 15 and 17, 2021, were selected. In the Mu Us Sandy Land, ten cloud-free Sentinel-2 surface reflectance images taken on September 1, 3, 6, and 8, 2021, were chosen. Subsequently, these images were fused using median compositing to create the final image for the four study areas. The satellite image time was synchronized with the sampling times. Moreover, the images were masked by applying desert and sandy land boundaries derived from the land cover data.

2.2.3 Unmanned aerial vehicle data

From September 26 to 15 October 2021, we conducted field surveys in the Mu Us Sandy Land and Gurbantunggut Desert. During this survey period, the impact of short-lived vegetation on the spectral data was minimized as they had mostly withered or died in the desert. Moreover, as it was a non-precipitation season, the weather conditions were generally clear and sunny. Multiple flight sample plots were carefully selected within each sample area, and each sample area covered an area of 3 × 3 km². Within each flight sample plot, a DJI Phantom 4 unmanned aerial vehicle (UAV) with a five-channel multispectral camera was used to shoot the sample images with a high spatial resolution of 0.02 m, and the flight paths of the UAV are shown in Figure 2a. The UAV was flown at a height of 30 meters, and shooting was performed under clear-sky, windless and optimal field-of-vision conditions. The flight duration ranged from approximately 15 to 20 minutes, covering an area of 0.03 km² in each flight. The UAV was equipped with four bands, including the blue band (450 ± 16 nm), green band (560 ± 16 nm), red band (650 ± 16 nm) and NIR band (840 ± 26 nm) (Wang et al. 2023). DJI Terra software was used for radiometric correction, image mosaicing and fusion, format conversion, and qualified orthophoto band image production based on the real-time kinematic technology of UAV to enhance the accuracy of geolocation matching with satellite imagery, as shown in Figure 2. Visual interpretation and supervised classification were used to extract the BSC areas from each orthophoto. We used object-based nearest neighbor classification, which is specifically divided into segmentation and classification. The UAV images were segmented in eCognition software, and the segmentation scale was set to 10. The mean, brightness, and standard deviation were selected as features for the supervised classification. In addition, the classification results were verified by randomly generating points and labeling the sample points by visual interpretation. We repeated these steps three times and took the average of the three results as the final BSC coverage result (Figure 2c). To correspond the satellite pixel area to the UAV image area, we superimposed the UAV images and the satellite images in the same coordinate system and segmented the UAV images along the boundaries of the satellite pixels to calculate the true BSC coverage of the satellite pixels (Figure 2d). The above steps were also applied to sandy land. Finally, each field area corresponding to a satellite pixel was taken as a sample. A total of 233 samples were obtained to assess the accuracy of BSC coverage estimation results, including 133 desert samples and 110 sandy land samples. All the elements in the UAV orthophoto image were divided into different categories, including BSC regions and non-BSC regions in the desert, as well as vegetation, MBSCs, and sandy soil in the sandy land. After performing visual interpretation, certain mixed pixels that closely resembled BSCs were excluded. Subsequently, a post-classification accuracy assessment was conducted using randomly selected validation points. Through meticulous visual inspection and analysis of the confusion matrix, it was determined that the average accuracy of BSC extraction reached 96%.

Figure 2. Rmote sensing images and classification results of biological soil crust areas in the desert. (a) a high-resolution true-color image of biological soil crusts taken by an unmanned aerial vehicle (spatial resolution: 0.02 m). (b) a Sentinel-2 red band image of biological soil crust areas. (c) classification results of biological soil crust areas based on eCognition software; the colored area represents biological soil crusts. (d) biological soil crust classification results for a single satellite pixel.

2.3 Feature spaces

2.3.1 NDVI-BI feature space for desert

In each desert pixel, we limited the number of endmember types to 3, namely, non-BSC, lichen and moss, because the vegetation coverage of this area was less than 4%. Non-BSC only includes sand and a small amount of plant litter. The ability to successfully determine the proportion of spectral contributions from different surface component types using multispectral data depends on the contrast of the spectra between these components. Therefore, the NDVI and BI were selected as the separation indices, which were defined by the following equations:

(1) $NDVI=\frac{NIR-\mathrm{R}\mathrm{e}\mathrm{d}}{\mathrm{N}\mathrm{I}\mathrm{R}+\mathrm{R}\mathrm{e}\mathrm{d}}$(1)

(2) $BI=\sqrt{Gree{n}^2+Re{d}^2+NI{R}^2}$(2)

where Green, Red and NIR represent the reflectances of the green, red and NIR bands, respectively.

The separation indices were redistributed to a plane space and used as the x-axis and y-axis to construct a desert feature space, as shown in Figure 3. Moss has a low BI value and high NDVI value, lichen has low NDVI and BI values, and non-BSC has a high BI value and low NDVI value. The reflectances of pure lichen, moss and non-BSC serve as three vertices in the feature space, forming a triangle, and the mixed spectrum is distributed inside the triangle.

Figure 3. The feature space conceptual model for deserts based on the normalized difference vegetation index (NDVI) and brightness index (BI). The photos taken in the field show different surface component types in the desert.

2.3.2 BSCI-NDVI feature space for sandy land

There were three types of surface components mainly distributed in the sandy land, including green vegetation, sandy soil and MBSCs. The NDVI and BSCI were selected as the separation indices, and the BSCI was defined by Eq.3 (Chen et al. 2005).

(3) $BSCI=\frac{1-L\ast \left|R_{red}-R_{green}\right|}{R_{GRNIR}^{mean}}$(3)

where $R_{green}$ and $R_{red}$ represent the reflectances of the green and red bands of the multispectral satellites, respectively. $R_{GRNIR}^{mean}$ is the average reflectance of the green, red and NIR bands. L represents the coefficient.

In the sandy land feature space constructed using the NDVI and BSCI, MBSCs exhibit high BSCI values and moderate NDVI values. Green vegetation is characterized by high NDVI values, while sandy soil has low NDVI values. Both green vegetation and sandy soil have low BSCI values (as shown in Figure 4). As a result, the scattered points in the feature space form triangles. Any mixed pixel containing all three surface component types is in the triangle.

Figure 4. The feature space conceptual model for sandy land based on the biological soil crust index (BSCI) and normalized difference vegetation index (NDVI). Field photos captured different surface component types in the sandy land.

2.3.3 In situ measured spectra assessment

The NDVI versus BI and BSCI versus NDVI scatter plots were plotted for deserts and sandy lands, respectively, based on the field-measured spectra of each surface component type. Subsequently, we conducted a comparison between the patterns observed in these scatter plots and the conceptual models constructed in subsections 2.3.1 and 2.3.2.

2.4 BSC coverage estimation and validation

2.4.1 Spectral unmixing model

The reflectance of a mixed pixel is contributed only by three surface component types with different abundances and is a linear combination of the reflectances of different endmembers, and the abundances of the different endmembers in a mixed pixel are taken as the endmember weights (Deng and Changshan 2013). Therefore, the pixel reflectance of desert and sandy land can be expressed with Eq. 4 and Eq. 5, respectively.

(4) $R_{desert}=f_{moss}\ast R{R}_{moss}+f_{sand}\ast R{R}_{sand}+f_{lichen}\ast R{R}_{lichen}$(4)

(5) $R_{sandy\hspace{0.17em}land}\hspace{0.17em}=\hspace{0.17em}{f}_{MBSCs}\ast R{R}_{MBSCs}+f_{sandy\hspace{0.17em}soil}\ast R{R}_{sandy\hspace{0.17em}soil}+f_{Green\hspace{0.17em}Vegetation}\ast R{R}_{Green\hspace{0.17em}Vegetation}$(5)

where $R_{desert}$ and $R_{sandyland}$ are the reflectances of desert and sandy land, respectively. $R{R}_{moss}$, $R{R}_{lichen}$ and$R{R}_{non-BSC}$ are the reflectances of pure moss, lichen and non-BSC, respectively. $f_{moss}$, $f_{lichen}$ and $f_{non-BSC}$ are the endmember coverages and are greater than or equal to 0.

When using indices instead of reflectance, the relationship between two separation indices and BSC coverage can be expressed using Eq. 6 as follows:

(6) $Index=f_1\ast Inde{x}_1+f_2\ast Inde{x}_2+f_3\ast Inde{x}_3$(6)

where $Index$ is the separation index value of mixed pixels. $Inde{x}_1$, $Inde{x}_2$ and $Inde{x}_3$ are index values of different endmembers. $Accessisdenied$ $f_2$ and $f_3$ are the coverages of different endmembers and are greater than or equal to 0.

Based on the above methods, the abundance of each surface component type can be determined using two separation indices of any pixel and the coordinates of the three vertices of the triangle. In deserts, the BSC coverage represents the combined coverage of moss and lichen in mixed pixels. By quantifying the coverage of each surface component type in each mixed pixel and averaging it across the study area, the area proportions can be determined.

The selection of endmembers was a parallel step to the above steps. Minimum noise fraction (MNF) rotation can determine the intrinsic dimensionality of the image data and separate the noise in the data to reduce the amount of computation required in subsequent processing. In this study, MNF rotation was applied to the Sentinel-2 images of the study areas using ENVI. Then, the pixel purity index (PPI) algorithm was used to determine the purity of the pixels based on the first four MNF bands. High PPI values represent the high purity of image pixels (Quintano, Fernández-Manso, and Roberts 2013); thus, points with higher PPI values were selected to calculate the average value of the separation indices. After visually inspecting each point on Google Maps, these points were input into the spectral unmixing model to calculate the BSC coverage.

2.4.2 Accuracy assessment

The accuracy assessment relies on comparing field-measured coverage results with estimated coverage results from validation points. The strength of their correlation is determined by calculating the correlation coefficient squared (R²). When R² approaches 1, it signifies a strong linear relationship between the variables, allowing for a simple linear model to effectively describe their connection. In addition, the root mean square error (RMSE) and normalized mean square error (NMSE) were used to evaluate the accuracy of the proposed method. These indices can measure the deviation between the estimated value and the observed value, and a low value indicates the high accuracy and predictive capability of the model. The RMSE and NMSE were calculated and expressed as percentages using the following equations:

(7) $RMSE=\sqrt{\frac{{\sum }_{i=1}^n{\left(p_i-q_i\right)}^2}{N}}\times 100\mathrm{\%}$(7)

(8) $NMSE=\frac{mean\left({\left(p-q\right)}^2\right)}{mean\left({\left(mean\left(q\right)-q\right)}^2\right)}\times 100\mathrm{\%}$(8)

where $p_i$ represents the estimated BSC coverage of sample i, $q_i$ represents the measured value of sample i, and N represents the total number of samples. $Mean$ represents the average values of the samples.

In addition, estimation accuracy (EA) was used to assess the prediction ability of a model. A high EA value indicates a higher estimation accuracy of the model.

(9) $EA=\left(1-\frac{RMSE}{Mean}\right)\times 100\mathrm{\%}$(9)

where $Mean$ represents the average values of the measured samples.

In addition, the validation samples were categorized based on their land cover types. Then, statistical data illustrating the deviation between the observed and predicted values were visualized through a confusion matrix and histogram.

3 Results

3.1 Ground unit spectrum

Figure 5a shows that the spectral slope and reflectance of BSCs are lower than those of sand and plant litter in the VIS-NIR wavelengths. This is due to BSCs being characterized by dark crusts, which contribute to the reduction in spectral reflectance. Additionally, BSC-covered areas tend to exhibit increased surface roughness, further lowering their spectral properties. When lichen and moss have equal coverage, the spectral reflectance of moss is significantly lower than that of lichen. Plant litter exhibits a spectral curve similar to sand but slightly higher in the VIS to NIR wavelengths. Moreover, moss exhibits significant absorption in the red band, primarily attributed to its higher chlorophyll content compared to lichen, particularly under high-coverage conditions. Conversely, lichen and sand display subtle spectral absorption at 680 nm (Figure 5b). The spectral curve of vegetation is different from those of other surface components, exhibiting absorption and reflection characteristics in the red and green bands, respectively, and high reflectance in the NIR band. When compared to sandy soil, the spectral reflectance of MBSCs and their slope between the green and red bands are lower because of the presence of moss in MBSCs.

Figure 5. Spectral curves based on a spectrometer for different surface components from 400 nm to 900 nm in the (a) desert and (b) sandy land. Marked points indicate the Sentinel-2 band reflectance for different surface component types.

The results of the above spectral analysis show that the selected indices are useful for distinguishing among different main surface component types in deserts and sandy lands. Indices were calculated based on the measured spectra derived from a UAV to construct the feature space, and the results are shown in Figure 6.

Figure 6. Spectral feature spaces in different land cover types, including (a) desert and (b) sandy land. The spectral data were derived from the field investigation.

In deserts, the BI can distinguish between BSC and non-BSC areas, and the NDVI can separate mosses from lichens. In the NDVI-BI space, mosses appear in the uppermost part of the triangle, with an NDVI value of 0.27 and a BI value range of 0.1–0.3 due to its dark black surface (Figure 6a). Non-BSC areas, located at the lower right, have a very high BI but low NDVI. Lichens, located at the lower left, have similar NDVI values to non-BSC. These three types of surface components show a triangular shape, which is also consistent with the conceptual model.

In sandy lands, further calculation of the index values for the three surface components reveals that MBSCs, vegetation and sandy soil also show a distinct triangle in the BSCI-NDVI feature space (Figure 6b). Vegetation is situated in the lower-right corner, MBSCs are located at the top, and sandy soil is located at the lower left.

3.2 Spectra unmixing of satellite images

Figure 7 shows that the scatter plots derived from Sentinel-2 satellite images all show the distribution characteristics of triangles, consistent with the conceptual model constructed in section 2.3 and the field spectral results provided in section 3.1.

Figure 7. The normalized difference vegetation index (NDVI)-brightness index (BI) and NDVI-biological soil crust index (BSCI) feature spaces for different land cover types based on the Sentinel-2 images. The scatter distribution results were presented using density distribution plots. Figure 7a shows the combined results for areas A (b) and B (c) in the desert. Areas C (d) and D (e) are located in the sandy land. The average index values of the endmembers in different regions were listed.

Area A in the desert is a moss-dominated area. There are only a few points on the oblique edge of the feature space (Figure 7b). Moss and lichen develop together in areas with high BSC coverage, but pure moss and lichen pixels are relatively few. In Area B of the desert, lichen is the main BSC type and has a very low NDVI value; thus, the scatter points obtained for these areas are located mainly below the triangle (Figure 7c). When the lichen coverage is high, the development of moss in this area causes the NDVI to increase. Based on the PPI results, the NDVI and BI values for pure lichen are 0.05 and 0.21, respectively. For pure moss, the values are 0.28 for NDVI and 0.19 for BI, while for pure sand, the NDVI and BI values are 0.09 and 0.83, respectively. In sandy land, scatter plots in the BSCI-NDVI feature space and the index values of the pure pixels show similarities among different areas (Figure 7d,e). By using the PPI and conducting field surveys, we determined that the pure MBSCs have an NDVI value of 0.22 and a BSCI value of 9.3 in Area C (Figure 7d). This very high BSCI value is mainly due to the low reflectance of such areas. Furthermore, the index values of the pure pixels in Area C closely resemble those in Area D (Figure 7e). However, slight deviations in the index values of the pure pixels between the two areas were observed. Sandy soil areas have high reflectance with very low BSC coverage and vegetation coverage, resulting in low BSCI values. In general, the distribution of all scatter plots aligns with our field investigation results and corresponds to the succession trend of BSCs.

3.3 Quantification of BSC coverage

The BSC coverage distributions in the desert and sandy land are shown in Figure 8. In the desert region, denoted as Area A, the BSC coverage reached 70.4%, while the sand coverage accounted for 28.1%; in Area B, the BSC coverage was 65.0%, with an area of 56.6 km² (Table 2). The sand coverage in Area B was relatively higher at 33.5%, indicating a larger expanse of exposed sand compared to Area A. The BSC distribution in both areas showcased a fragmented pattern with multipoint development and striped characteristics. Particularly, well-developed BSCs were observed in the lowlands between hills with low and gentle terrain. In the sandy land, the dominant surface component type in Area C was sandy soil, accounting for 49.0% coverage. Following sandy soil, the coverage of MBSCs was measured at 34.4%. Vegetation coverage was relatively low, amounting to only 15.8% in this area. In Area D, the BSC covered an area of 92.4 km², accounting for 31.8% of the total area (Table 2). In contrast, vegetation coverage in Area D was measured at 14.1%, indicating sparse vegetation growth. Sandy soil coverage was more prominent in this area, occupying 52.1% of the total area.

Figure 8. Biological soil crust coverage distribution map in the four study areas.

Table 2. Biological soil crust coverage and area statistics in different regions.

Region	BSC coverage (%)	BSC area (km^2)
Area A	70.4%	162.1
Area B	65.0%	56.6
Area C	34.4%	73.7
Area D	31.8%	92.4

3.4 Accuracy assessment

Figure 9 displays the comparison between the predicted BSC coverage classes and the actual BSC coverage classes based on the field survey data. The high values along the diagonal and the low values in the off-diagonal cells highlight the reliability and accuracy of the proposed method, as well as the consistency of the results with field coverage data. These findings demonstrate that multispectral data could be used directly to estimate BSC coverage. As shown in Figure 10, the histograms showed that there was a difference in the estimation abilities of the prediction models between desert and sandy land. In the desert, the BSC coverage estimation result has an R² value of 0.86, with an RMSE of 10% and an NMSE of 18%. The statistical results of the validation points showed that the method had a good ability to estimate the BSC coverage in deserts (EA value is 79%). Among the 133 validation points, the deviation between the predicted and observed values was within 10% for 101 points (Figure 10a). In the sandy land, an R² value of 0.89 indicated a strong correlation between the predicted and observed values. The RMSE and NMSE of the BSC coverage estimation result were 8% and 13%, respectively, with an EA of 81%. Figure 10b shows that the BSC coverage of 71 points was underestimated using the proposed method. Compared to the validation results of the desert, the deviation range between the estimated coverage and observed coverage at the validation points of the sandy land was smaller.

Figure 9. Confusion matrices between observed and estimated biological soil crust coverage for the (a) desert and (b) sandy land.

Figure 10. Histograms of the estimated errors of the proposed biological soil crust coverage models for the (a) desert and (b) sandy land.

4 Discussion

The distribution and composition of BSCs are crucial for understanding desertification and climate change because BSCs are highly related to ecosystem processes such as the carbon, nitrogen and water cycles. However, mapping the distribution of BSC coverage is challenging. One of the main difficulties faced in this study is to develop appropriate models for different land cover types and large-scale rapid coverage estimations. Based on the spectral differences between BSCs and other surface components in the desert and sandy land, two such models have been developed: the feature space model and the spectral unmixing model. The feature space model takes into account the unique spectral characteristics of BSCs, allowing it to differentiate them accurately. On the other hand, the spectral unmixing model calculates the proportions of different components present in mixed pixels. By combining the feature space model with the spectral unmixing model, high accuracy can be achieved in estimating BSC coverage at the subpixel scale for both land cover types. This approach has the potential to enable large-scale estimations of BSC coverage and provide accurate ecological assessments.

4.1 Evaluation of the feature space

In this study, we mapped satellite multidimensional spectral data into a two-dimensional feature space. This transformation not only reduces the dimensionality of the data but also highlights the most important patterns of variation and key features within the data. The clear triangular relationships observed in the feature space indicate the correlation and differences in spectral characteristics among different surface component types. By observing the shape and position of triangles in the feature space, the identification and separation of BSC and other surface component types can be achieved.

With the development of BSCs, the absorption in red wavelengths increases and reflectance decreases (Chamizo et al. 2012). Previous studies have shown that the NDVI and BI are both significantly correlated with the growth and succession of BSCs, and these indices exhibit pronounced differences among moss, lichen and non-BSC (Escadafal and Bacha 1996; Rodríguez-Caballero, Knerr, and Weber 2015; Zaady, Karnieli, and Shachak 2007). As the biomass of BSCs increases and subsequent successional species emerge, the increase in chlorophyll a and colored pigments changes the surface color. An important prerequisite for applying the spectral unmixing model is that the indices must be linear with endmember coverage. Moss coverage is positively correlated with the NDVI (Figure 11a). This finding has been supported by earlier studies (Fang and Zhang 2011). On the other hand, the BI reflects changes in BSCs with dark crusts and is negatively correlated with BSC coverage. Consequently, as BSCs cover the surface, the surface color decreases. However, the correlation between reflectance and BSC succession is low (Escadafal and Bacha 1996). The NDVI is an important indicator used to characterize BSC succession, which gradually increases with the development of BSCs (Zaady, Karnieli, and Shachak 2007). Zaady, Karnieli, and Shachak (2007) found that BSC succession can be expressed as a linear regression of the NDVI. Thus, it is reasonable to use the triangular feature space of the NDVI and BI to distinguish between moss, lichen and non-BSC.

Figure 11. Scatter plots of spectral indices and biological soil crust coverage: (a) the normalized difference vegetation index (NDVI) and moss coverage in desert and (b) the biological soil crust index (BSCI) and mixed biological soil crusts coverage in sandy land. Spectral indices were derived from measured spectrometer data. The red lines represent the linear relationships between the indices and biological soil crust coverage.

The sandy land feature space was constructed using the BSCI and NDVI to estimate the BSC coverage of sandy land. Chen et al. (2005) considered that the slope of the spectral curve of BSC in the green and red bands is flatter than that of sand and green vegetation and that BSCs have much lower reflectance in the VIS-NIR bands than sand and green vegetation; they demonstrated that the BSCI can effectively extract BSC areas. Figure 11b shows that the MBSCs coverage exhibits a significant linear relationship (P < 0.01) with the BSCI in sandy land. The BSCI was developed to monitor BSCs using remote sensing and can effectively distinguish sand from BSCs (Chen et al. 2005). The NDVI is still an appropriate indicator for monitoring green vegetation, and a linear relationship between vegetation coverage and the NDVI has been demonstrated (Daughtry, Hunt, and McMurtrey Iii 2004).

4.2 Advantages and limitations of the proposed method

When predicting BSC coverage using multiple environmental variables related to BSCs, obtaining a high prediction accuracy requires high-quality and high-resolution environmental data (Beaugendre et al. 2017). In contrast, our method of directly observing BSCs using satellite images is more reliable and simpler. When establishing the feature space, we considered using the band of the Multispectral Scanner and TM sensors. The advantage of this method is that any satellite sensor containing the VIS and NIR bands can be used to estimate BSC coverage. This is particularly beneficial for calculating long-time series of BSC coverage since it allows for the inclusion of historical satellite imagery before the availability of more advanced satellite sensors (e.g. before 2000). Considering the long-term span of Landsat images, the study of dynamic changes in BSC coverage over extended periods can be achieved. Furthermore, using VIS-NIR bands can help mitigate the impact of plant litter on BSC coverage accuracy. This is because sand and plant litter often exhibit spectral similarities.

Although the spatial resolution of multispectral sensors is lower than that of hyperspectral sensors, the study has successfully demonstrated that multispectral satellite images can be effectively used to map BSC coverage in desert and sandy land regions with satisfactory accuracy. Since the reflectance of surface components and the spectral signals received by the satellite can be affected by other factors, it is recommended to choose the average index value of multiple pure pixels for each image that requires unmixing. Fortunately, the availability of higher-resolution satellite images, such as Landsat (30 m) and Sentinel-2 (10 m), significantly increases the possibility of pure pixels being present in a satellite image. This enables the quick identification of triangles on feature spaces based on two indices and scatter plots, thereby reducing the computational complexity of the spectral unmixing model. Compared to machine learning or deep learning methods, the proposed method in this study requires only a small amount of prior knowledge. Therefore, it is convenient for large-scale applications and long-time series mapping. However, the proposed method in this study was recommended to be used in sparse vegetation areas with vegetation coverage lower than 20%, such as deserts and sandy lands, where vegetation does not dominate. In other land covers, the effect of vegetation on the spectrum of mixed pixels may lead to erroneous estimation of BSC coverage. Therefore, further research is necessary to investigate the applicability and limitations of this method in such areas.

High temperature and intense solar radiation accelerate evaporation from the surface, resulting in a reduction in water in the BSC and a slowdown or even cessation of physiological metabolic activity. These changes will result in a decrease in the spectral absorption of BSC in the red band and affect the determination of BSC coverage. Thus, we chose autumn images instead of summer images for monitoring BSC coverage. Moreover, selecting non-growing season images in deserts can further mitigate the influence of small amounts of short-lived vegetation on the spectrum.

Based on the validation results, we have identified potential reasons for the overestimation or underestimation of BSC coverage in certain areas. These discrepancies are likely attributed to the following factors. The presence of a few shadow areas caused by plant litter and rough surfaces results in low BI values and high BSCI values, respectively, thus causing the estimated BSC coverage to be higher than the actual value. In contrast, a small amount of saline alkali land in the desert will lead to the opposite effect, causing the coverage to be underestimated. We believe that the constructed model is not well suited for estimating BSC coverage in areas with shadows and saline alkali land. In sandy lands, there may be a small amount of sand and plant litter that exhibit higher brightness compared to sandy soil. This higher surface reflectance contributes to lower BSCI values, leading to an underestimation of BSC coverage.

Although this study achieved fine accuracy using linear unmixing models, it has some limitations. The linear spectral unmixing method used in this study only allows one fixed spectrum for each type of endmember. Consequently, it does not account for natural variations in surface properties under different scene conditions, meaning that the same type of surface component may exhibit different spectra. This limitation could lead to decreased accuracy when applying the proposed method on a larger scale. Therefore, it is recommended to use the MESMA model to replace the spectral unmixing model used in this study. The advantage of MESMA is its ability to address spectral variability by allowing the number and types of endmembers to vary on a per-pixel basis, overcoming the limitations of linear SMA (Quintano, Fernández-Manso, and Roberts 2013; Roberts et al. 1998). However, this method also greatly increases the computational workload.

4.3 Ecological significance of this work

Due to the lack of large-scale spatiotemporal data on nonvascular plants, the understanding of the impact of these communities is still limited. A novel BSC coverage estimation method proposed in this study allows us to accurately calculate BSC coverage and observe the spatial heterogeneity of BSC distribution. The proposed method can refine and enhance the results reported by Wang et al. (2022), providing more accurate and reliable BSC coverage data in the study area. Furthermore, the abundance map of BSCs offers additional valuable information that may not be fully captured through land cover maps. Such monitoring methods can enable decision makers to accurately understand the situation of natural resources in the region they manage. Moreover, BSC abundance is an important environmental evaluation indicator, and changes in BSC abundance are related to the environment in which the BSCs are located. A decrease in BSC abundance without human interference may signal the degradation of BSCs, leading to a transition from a continuous distribution to a fragmented one, ultimately culminating in their disappearance. Such changes indicate shifts in environmental factors. Thus, it is also necessary to calculate the previous BSC coverage and dynamic changes in the BSC coverage in combination with long-time series data, which will help to consider the scale and complexity of the impact of natural and human factors on BSCs and to assess the impact of climate change and environmental protection policies on the BSC, further allowing the potential future scale of BSCs to be predicted and simulated. In this process, the proposed coverage estimation method combined with short-revisit-period and wide-coverage multispectral data can solve this demand. When analyzing the changes in BSCs, remote sensing data derived from many multidimensional remote sensing sensors can be used to assess natural environmental changes, e.g. land cover products, primary net productivity products, and ice and snow products (Morel et al. 2019; Romanov 2017; Sulla-Menashe et al. 2019). These data can help researchers quickly simulate the interaction mechanism between the environment and BSCs in drylands.

4.4 Uncertainties and future development

The hydration status of BSCs can be influenced by various factors, introducing uncertainties in their spectral assessment. These uncertainties may result in some deviations in the measurement results, particularly in certain regions. When conducting large-scale monitoring, it becomes crucial to consider the hydration status of BSCs at the time of image acquisition to ensure accurate results. In addition, the proposed feature spaces have certain limitations and may not be directly applicable to regions with different types of surface components beyond the three types mentioned in this study. Due to the limitation of road conditions, the sample distribution was mainly along main roads; thus, the uniformity of the samples was limited to some extent. If other convenient tools, such as long-distance UAVs, could be used for sampling or if long-term continuous sampling can be carried out, the obtained BSC spectral data would likely be more uniform. In some areas where field surveys are difficult, the availability of field data may be limited, which affects the accuracy validation of the quantitative mapping results. Therefore, further research and refinement of the methodology are necessary to address these limitations effectively.

In the future, we plan to apply this method to large-scale drylands to understand the distribution characteristics of BSCs in drylands. Future research directions include using long time series images to obtain the dynamic changes in BSC coverage and using multi-source data to evaluate the relationships between the environmental factors and BSC changes.

5 Conclusion

In this study, we conducted an assessment of the accuracy achieved by combining spectral feature spaces with a spectral unmixing model to quantify BSC coverage in the desert and sandy land of China using multispectral remote sensing imagery. We introduced two feature spaces, namely NDVI-BI and BSCI-NDVI, which demonstrated the triangular relationship between the spectral mapping results of different surface components in mixed pixels of desert and sandy land, respectively. Combining a spectral unmixing model based on the feature space is an effective method to quantify the BSC coverage. All these tasks can be performed on the multispectral sensor to facilitate the implementation of BSC coverage mapping at large scales and to accurately calculate the BSC area in deserts or sandy lands. These findings provide valuable data for quantitatively assessing the improvement abilities and ecological roles of BSCs in drylands. Moreover, they contribute to enhancing our understanding of the significance of BSCs in desert and sandy land ecosystems.

Acknowledgments

This research was financially supported by the National Natural Scientific Foundation of China (Grant No: 41991232), the Strategic Priority Research Program of the Chinese Academy of Sciences (XDA19050504-01) and the National Key Research and Development Plan of China (Grant No: 2016YFC0500201).

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability statement

Data available on request from the authors.

Additional information

Funding

This work was supported by the National Natural Scientific Foundation of China [41991232].