Synergistic Coupling of Multi-Source Remote Sensing Data for Sandy Land Detection and Multi-Indicator Integrated Evaluation

Abstract

Accurate and timely extraction and evaluation of sandy land are essential for ecological environmental protection; it is urgent to do the research to support the sustainable development goals (SDGs) of Land Degradation Neutrality. This study used Sentinel-1 Synthetic Aperture Radar (SAR) data and Landsat 8 OLI multispectral data as the main data sources. Combining the rich spectral information from optical data and the penetrating advantages of radar data, a feature-level fusion method was employed to unveil the intrinsic nature of vegetative cover and accurately identify sandy land. Simultaneously, leveraging the results obtained from training with measured data, a comprehensive desertification assessment model was proposed, which combines multiple indicators to achieve a thorough evaluation of sandy land. The results showed that the method based on feature-level fusion achieved an overall accuracy of 86.31% in sandy land detection in Gansu Province, China. The integrated multi-indicator model C22_C/FVC is the ratio of correlation texture features of VH to vegetation cover based on which sandy land can be classified into three categories. When C22_C/FVC is less than 2.2, the pixel is classified as fixed sandy land. Pixels of semi-fixed sandy land have an indicator value between 2.2 and 5.2. Shifting sandy land has values greater than 5.2. Results showed that shifting sandy land and semi-fixed sandy land are the predominant types in Gansu Province, with 85,100 square kilometers and 87,100 square kilometers, respectively. The acreage of fixed sandy land was the least, 51,800 square kilometers. The method presented in this paper is robust for the detection and evaluation of sandy land from satellite imageries, which can potentially be applied for conducting high-resolution and large-scale detection and evaluation of sandy land.

1. Introduction

The natural ecological environment is increasingly being disturbed by human activities such as economic development and resource exploitation. Areas of sandy land are expanding, posing a serious threat to the sustainable development of the ecological environment and human life [1,2]. Additionally, it can lead to global issues such as biodiversity loss and soil degradation [3,4,5]. Land desertification in arid and semi-arid areas is one of the most prominent ecological environmental issues. It is an important focus of attention and active management [6,7,8]. Sandy land refers to areas where the surface is covered by sand dunes or sand in semi-arid or semi-humid regions, characterized by strong winds, limited water flow, and sparse vegetation. Sandy land is widespread in northern China. In order to assess the effectiveness of ecological restoration projects, it is necessary to timely evaluate the current status and dynamics of sandy land, which will enable more efficient and targeted implementation of policies and plans for desertification control [9,10]. Indeed, accurately determining whether a region has undergone desertification and effectively extracting and evaluating information about sandy land is crucial for regional, ecological, environmental protection, and sustainable development, as the implementation of various measures to prevent and control desertification is essential for achieving these goals.

In recent years, remote sensing technology has provided a more objective and accurate data foundation for detecting and evaluating sandy land, thanks to its wide observation coverage, real-time capability, and strong dynamics [11,12,13,14]. It has become an indispensable tool for detecting and evaluating sandy land information at regional and even global scales [15]. Various remote sensing monitoring methods have been proposed. In the field of detecting sandy land, methods based on soil characteristic parameters [16,17,18,19,20], single sensors [21,22], and fused multi-sensor approaches have been proposed. In the context of characterizing sandy land, methods based on characteristic parameters yield high accuracy. However, there are challenges when it comes to large-scale extraction of desertification information. This approach requires on-site measurement of soil characteristic parameters, which can be time-consuming and resource-intensive. Additionally, data acquisition takes a longer time, resulting in poor timeliness. The use of a Normalized Difference Snow Index (NDSI) based on the red and SWIR2 bands of Landsat 5 was proposed [23], which effectively highlights and identifies the shifting sand and sand dunes. However, it faces challenges in identifying mixed pixels, especially those that consist of a combination of sand and vegetation. Therewith, the Normalized Difference Sand and Infrared Index (NDSAI) was proposed for the rapid identification of sandy and non-sandy areas [24], but it shows relatively poor performance in distinguishing between sandy areas and buildings. On the other hand, the Normalized Difference Exposed Sand Index (NDESI) was then introduced [25], which is based on the blue, red, and two shortwave infrared bands. This index significantly improves the detection accuracy for arid bare sand areas. Afterward, Wu et al. conducted research on the decomposition of mixed pixels using multispectral remote sensing data to detect information about sandy land. However, the problem of vegetation obscuring the underlying sandy soil in desertified areas has not been adequately resolved [21]. Li et al. focused on the advantages of multiple sensor fusion and proposed a feature-level fusion method that can effectively identify the underlying sandy soil through vegetation. However, obtaining fully polarimetric SAR data and hyperspectral data at a high resolution has proven to be challenging, making it difficult to achieve large-scale detection of sandy land information [26]. The study conducted by Zhao et al. focused on the polarization characteristics of various polarization decompositions and selected the optimal planar features for achieving object classification, which exhibits strong regional specificity [27]. The research conducted by Di et al. employed both dual-polarization and fully polarimetric radar data for land cover classification, but in environments with diverse and complex land cover types, there may be instances of mixed or misclassified results [28]. Wang et al. proposed a new classification scheme for mud and sand flats on intertidal flats using fully polarimetric SAR data. It can reduce the misclassification problem between mudflats and sand flats, but their poor performance in distinguishing targets with similar scattering characteristics [29]. Zerrouki et al. provided an effective approach to detect deserted regions based on Landsat imagery and Variational AutoEncoder, which makes it suitable for resolving complex recognition problems like separating desertification cases from like-desertification pixels [30]. Wei et al. propose a refined desertification information extraction method based on multisource feature spaces and geographical zoning modeling. It can be used to realize the transformation from the surface reference variables with biophysical meanings into remote sensing spectral information and complete the quantitative calculation and automatic classification of desertification information [31]. Tan et al. proposed an unsupervised classification method based on fully polarimetric radar data, called ND-LSC (Unsupervised Classification based on Polarimetric Radar Data). It combined the advantages of the ND strategy, ASLIC algorithm, and LSC algorithm to effectively obtain polarimetric parameters in areas with diverse land cover types and complex structures, such as the Otingdag Sandy Land [32].

Sandy land evaluation involves the classification of sandy land into fixed sand, semi-fixed sand, and shifting sand based on information such as vegetation cover. The area of each type of sandy land is evaluated to assess the desertification status of the area under investigation. The introduction of vegetation indices such as the Normalized Difference Vegetation Index (NDVI) by Yucker et al. marked a new era in the establishment of remote sensing monitoring indicator systems [33]. Vegetation indices are considered important indicators in the monitoring and assessment of desertification as they effectively distinguish between surfaces with vegetation cover and bare soil. Therefore, they serve as crucial factors in the evaluation of desertification and its progression. However, in sandy areas, vegetation is very sparse, and the range of vegetation cover change thresholds is limited. Therefore, relying solely on NDVI cannot objectively reflect the development process of desertification in arid and semi-arid regions. Based on the changes in the physical properties of topsoil, Xiao et al. employed the Topsoil Grain Size Index (TGSI) to indicate soil texture for the monitoring of desertification in arid regions [34]. By combining remote sensing data with ground measurements, a linear regression model was established by Huo et al. based on NDVI and Fractional Vegetation Cover (FVC). The sandy land assessment was accomplished through the classification of vegetation cover levels [35,36]. Subsequently, Lamchin et al. constructed a land desertification assessment system based on three indicators: NDVI, Albedo, and TGSI. They analyzed the changes in land desertification in Mongolia as well as the influencing factors [37]. Changes in sandy land in Mexico were analyzed by Becerril-Pina et al. through the introduction of the Bare Soil Index (BSI) and NDVI as analytical tools [38]. A method proposed by Wang et al. incorporated FVC, Enhanced Vegetation Index (EVI), and Modified Soil Adjusted Vegetation Index (MSAVI) to establish a decision tree for land desertification classification using appropriate thresholds [39]. However, this method selected vegetation indices, which have high internal correlations and lack distinct representativeness.

To address the challenge of identifying sandy land with vegetation cover and the challenge of lacking unified evaluation indicators for desertification monitoring, this study focuses on the following: (1) combining both spectral and radar data by taking advantage of rich spectral information from multispectral data and the penetrating characteristics of radar data to employ a feature-level fusion method to unveil the intrinsic nature of vegetative cover and accurately identify sandy land; (2) analyzing the distribution patterns of different levels of desertification samples based on multiple indices derived from optical and radar data. The type of desertification for each sample point is determined based on the measured vegetation cover data and indicators for desertification classification; (3) to propose an integrated multi-indicator model for evaluating sandy land, enabling the completion of sandy land assessment.

2. Study Area and Data Description

2.1. Study Area

The study area is in the entire Gansu Province, northwest China, bounded between 32°11′N and 42°57′N latitude and 92°13′E and 108°46′E longitude. The province is adjacent to the Tengger Desert and the Badain Jaran Desert. The northern part of Gansu Province is mainly covered by shifting, semi-fixed, and fixed sand, forming a natural landscape of sandy grasslands. The southern part has higher altitude, abundant rainfall, and higher pasture productivity and serves as a major livestock farming base. It is divided into 12 prefecture-level municipalities and 2 autonomous prefectures. Gansu has a diverse range of climate types, including subtropical monsoon climate, temperate monsoon climate, temperate continental (arid) climate, and plateau alpine climate, from south to north. The annual precipitation in different areas of Gansu ranges from 36.6 to 734.9 mm, and the average annual temperature ranges from 0 to 15°C. Most areas in the province have a dry climate, and arid and semi-arid regions account for 75% of the total area. The major land cover types in Gansu Province include arable land, orchards, forests, grasslands, wetlands, urban and industrial land, transportation land, water, and water conservancy facilities. The location of the study area is shown in Figure 1.

2.2. Data Acquisition and Preprocessing

2.2.1. Remote Sensing Data

The remote sensing data used in this study include Sentinel-1 SAR data and Landsat 8 Operational Land Imager (OLI) multispectral data. The specific image information is shown in

Table 1. Basic information of data sources.

Data Type	Imaging Time	Spatial Resolution	Number of Scenes
Sentinel-1	2021.07–2021.08	10 m	34
Landsat 8 OLI	2021.07–2021.08	30 m	50

. Sentinel-1 is composed of two polar-orbiting satellites, namely Sentinel-1A and Sentinel-1B. The sensors onboard are active microwave remote sensing sensors. The Sentinel-1 data used in this study was acquired in Interferometric Wide (IW) mode. In this mode, the radar beam can be electronically steered from aft to fore in each burst’s azimuthal direction, in addition to controlling the beam’s range like in a scanning radar. This avoids the fan effect and results in uniform image quality across the entire swath. The data are acquired in dual-polarization (VV+VH), and the data product is in Single Look Complex (SLC) format, which is a Level-1 product that provides both phase and amplitude information. The Landsat 8 OLI sensor consists of 9 spectral bands with an imaging swath of 185 km. The OLI sensor includes a narrow panchromatic band (Band 8) that allows for better discrimination of vegetation and non-vegetation features in the pan-sharpened image. The rich spectral information provided by Landsat 8 OLI data supports the inversion of multiple indicator factors in this study. The acquisition time of the remote sensing data in this study is coincident with the time of field data collection, ensuring the accuracy validation of the results.

2.2.2. Data Preprocessing

In this study, Landsat 8 OLI data was obtained via the Google Earth Engine (GEE) platform. Due to the large extent of the study area, there was a significant amount of data to be processed. The GEE platform eliminates the need to download and preprocess a large amount of data, as it allows for batch processing of satellite imagery on the Google Cloud. This capability enables processing without being constrained by spatial or temporal limitations. Therefore, the preprocessing, mosaicking, clipping, and generation of multiple index products were carried out on the GEE platform, significantly reducing the time and cost associated with data processing.

Processing of Sentinel-1 data was conducted using the SNAP software developed by the European Space Agency (ESA). The specific preprocessing workflow is as follows: (1) Applying Orbit File: The orbit vectors contained in the SAR image metadata are typically not highly accurate. Therefore, precise orbit data acquired from subsequent publications are used for optimization during this step. (2) Calibration: The purpose of radiometric calibration is to ensure that pixel values in the SAR image directly reflect the backscattering coefficient of the scatterers. (3) Deburst: During the deburst process, adjacent sub-swaths with the same time tag are merged in the range direction. In the azimuth direction, sub-swaths are merged based on their zero Doppler time. (4) Generation of a Polarimetric Matrix C2: Since the Sentinel-1 satellite has only dual-polarization channels, only the C2 covariance matrix can be generated. The C2 matrix is a complex conjugate symmetric matrix, and it only requires the recording of elements C11 (real), C12 (complex with separate storage for real and imaginary parts), and C22 (real) for subsequent polarimetric decomposition. (5) Multi-looking: Multi-look processing helps reduce the speckle noise in SAR images. (6) Range Doppler Terrain Correction: Due to the actual changes in terrain and the tilt of the SAR sensor, SAR data often exhibit range deformations. The purpose of terrain correction is to eliminate these deformations and to align the geometric features of SAR images more consistent with the real world.

2.2.3. Field Data

The selection of sampling sites was of crucial importance. The sampling sites were determined following three principles: representativeness, inclusiveness of all grassland types, and inclusiveness of sampling sites with high, medium, and low vegetation cover levels. The size of each sampling plot was 1 m × 1 m, and the survey included information on grass species, soil type, vegetation cover, etc. Vegetation cover was mainly measured by the “visual estimation method” and/or the “acupuncture method”. The acupuncture method measures vegetation cover from the statistics of the presence/absence of vegetation cover at the same intervals according to a certain distance rule along two transect lines.

A total of 314 field measurements were obtained in Gansu Province (Figure 2), sufficient for validating the results, among which 104 field measurements were sandy lands and the other 210 field measurements were non-sandy lands.

3. Methodology

3.1. Overall Approach

In the process of sandy land monitoring, focusing on the advantageous features of optical and radar data highlights both the multispectral thematic information and effectively eliminates or suppresses irrelevant details. This approach also preserves the target’s backscattering information in synthetic aperture radar (SAR) data, allowing for the detection of vertical structural information of land features. Consequently, it enhances the accuracy of sandy land monitoring. By analyzing the distribution of measured data on these features and utilizing machine learning techniques to extract important features, the evaluation of sandy land was accomplished. Based on the pre-processing of multi-source remote sensing data, we first detected sandy land using Landsat 8 OLI images acquired from a single sensor. The abundance of sandy land was obtained using a spectral unmixing method, and an appropriate threshold was selected for quantitative detection of sandy land. Secondly, Sentinel-1 data was used to detect sandy land by employing polarization decomposition techniques. Polarization decomposition features were obtained, and a support vector machine (SVM) classification method was used to detect sandy land. Subsequently, the Landsat 8 OLI images and Sentinel-1 radar data were fused, at pixel-level and feature-level, separately used for sandy land detection. The pixel-level fusion was conducted using three algorithms: Hue Saturation Value (HSV), Gram Schmidt (GS), and Principal Component Analysis (PCA). Finally, the accuracy of the six sandy land detection methods was assessed using field sample data. The best method was then selected to detect sandy land in Gansu Province. Next, textural features and index inversions were extracted based on Sentinel-1 radar data and Landsat 8 optical data. The distribution of sample points with different types of sandy land in different indicators was analyzed. The random forest algorithm was used to determine the effective index factors that accurately reflect the type of land desertification. Based on the distribution of these effective factors in two-thirds of the field-sampled data, an optimal threshold was determined. An integrated multi-indicator model was constructed to complete the evaluation of sandy land in Gansu Province. Finally, the results were validated using the remaining one-third of the field sample data. The research workflow is depicted in Figure 3.

3.2. Sandy Land Detection

3.2.1. Sandy Land Detection Based on Optical Data

Firstly, the pure pixel index (PPI) method was employed to extract pure pixels and acquire endmember spectral curves to prepare for spectral unmixing. Secondly, the linear spectral unmixing (LSU) algorithm was applied to obtain the abundance information of land cover based on the spectral characteristics of different materials in the Landsat images. Lastly, an appropriate threshold was selected to detect sandy land.

The extraction of endmembers is a critical step in the process of mixed pixel decomposition [40], as it directly affects the accuracy of the final unmixing results. In this study, the process of PPI was implemented using the ENVI software package, involving several steps. The image was first subjected to Minimum Noise Fraction (MNF) transformation to whiten the noise while reducing the dimensionality of the remote sensing data. Then, the PPI algorithm was applied to calculate the purest spectrum for each pixel by iteratively mapping the n-dimensional scatterplot to a random unit vector. The extreme pixels obtained from each mapping were recorded, and the Digital Number (DN) value of each pixel represented the frequency of being marked as an extreme pixel. Therefore, pixels with higher DN values indicated higher purity, allowing for the extraction of pure endmembers from the image. The LSU method can extract the abundance information of land cover types from remote sensing images based on their spectral characteristics. The result of the unmixing process is a series of grayscale images representing the spectral signatures of the end members. The pixel values in the output images indicate the proportion of each endmember spectrum within that pixel.

3.2.2. Sandy Land Detection Based on Radar Data

To reveal the major scattering characteristics of the targets in radar data, polarimetric decomposition can be performed to separate the polarization features caused by different scattering mechanisms of the objects. In this study, polarimetric decomposition was performed on preprocessed Sentinel-1 data. Since Sentinel-1 is dual-polarization radar data, only H/A/α decomposition (

Table 2. Characteristic parameters corresponding to each polarization decomposition.

Decomposition Method	Extract Features	Feature Meaning
H/A/α	α	α is the average polarization scattering angle of H/A/α decomposition, identifying the main scattering mechanism.
	H	H is the polarization entropy of H/A/α decomposition, which measures the degree of polarization.
	A	A is the anisotropy of H/A/α decomposition, which measures the relative magnitude of non-dominant scattering.

) could be implemented. The dual-polarization SAR system uses X polarization for transmission and simultaneously receives the backscattered echoes in both X and Y polarizations, forming an XX/XY dual-polarization working mode. In this mode, X can be either horizontal (H) or vertical (V) polarization, while Y can also be either H or V. Therefore, common dual-polarization combinations include HH/HV, HH/VV, and VV/VH. In this study, the dual-polarization SAR data used was in VV/VH mode. Therefore, taking VV/VH mode as an example, the polarization scattering matrix of the dual-polarization SAR data is as follows (1) [41]:

$$S_2=\left[\begin{array}{cc}{S}_{VV}& S_{VH}\\ 0& 0\end{array}\right]$$

(1)

The target vector k can be expressed as follows (2):

$$k={\left[\begin{array}{cc}{S}_{VV}& S_{VH}\end{array}\right]}^T$$

(2)

Therefore, the covariance matrix for this mode is as follows (3):

$$C_2=〈{kk}^{*T}〉=\left[\begin{array}{cc}〈{\left|S_{VV}\right|}^2〉& 〈S_{VV}{S^{*}}_{VH}〉\\ 〈S_{VH}{S^{*}}_{VV}〉& 〈{\left|S_{VH}\right|}^2〉\end{array}\right]$$

(3)

The matrix C2 is a complex conjugate symmetric matrix. The diagonal elements C11 and C22 represent the backscatter coefficients of the VV and VH polarization channels, respectively, in the dual-polarization SAR data. The real parts of C11 and C22 in the covariance matrix are used in sandy land detection. Five feature parameters were selected for extracting sandy land information based on radar data: the H/A/α polarization features of Sentinel-1 and the C11 and C22 features of the covariance matrix C2. Multiple sample data points were uniformly selected for four land cover classes, i.e., water, vegetation, buildings, and sandy land. The SVM classification method was used to classify the land cover and detect sandy land. The SVM classification method has been widely used in land cover classification studies due to its advantages of requiring fewer samples and achieving high classification accuracies. The greatest advantage of this method is that it does not require data dimensionality reduction during classification and still offers high performance in terms of algorithm convergence, training speed, and classification accuracy.

3.2.3. Synergistic Coupling of Multi-Source Data for Sandy Land Detection

The fusion of Landsat 8 OLI multispectral images and Sentinel-1 data allows for the highlighting of multispectral thematic information, the elimination or suppression of irrelevant information, and the retention of the backscatter information of targets in SAR data, enabling the detection of the vertical structural information of objects and thus improving the classification results. In this study, firstly, the spectral information of land cover was extracted from the Landsat 8 data, based on the bands sensitive to sandy areas. These spectral features were then combined with the radar data through pixel-level fusion and feature-level fusion techniques. Next, suitable training samples were selected, and the SVM classification method was employed to obtain the information about sandy areas.

In order to achieve accurate extraction of sandy areas, three methods, namely HSV, PCA, and GS, were chosen for pixel-level fusion (

Table 3. Description of fusion methods.

Fusion Methods	Fusion Effect
HSV Fusion	The edge information of the multispectral image, the target spectrum information, and the high-resolution features of the panchromatic image were retained, and the texture details of the image were enhanced.
PCA Fusion	PCA is efficient in data compression and the first component contains the majority of information. When the panchromatic image was used to replace the first principal component for inverse transformation, the phenomenon of spectral distortion appeared to a certain extent.
GS Fusion	The spectral information of the original multi-spectral image can be maintained, the spatial information was also significantly enhanced, and the spectral fidelity effect was better.

).

3.3. Integrated Multi-Indicator Model for Evaluating Sandy Land

An integrated multi-indicator model was constructed to evaluate sandy land based on the detected sandy land information in Gansu Province. Firstly, evaluation indicators for sandy land were derived from Landsat 8 OLI multispectral remote sensing data, including FVC, Modified Soil-Adjusted Vegetation Index (MSAVI), Enhanced Vegetation Index (EVI), Albedo, and Bare Soil Index (BSI). Secondly, textural features based on gray-level cooccurrence, including Mean, Homogeneity, Entropy, Energy, Dissimilarity, Contrast, and Correlation, were extracted from Sentinel-1 data in VV and VH polarizations. The distribution patterns of different sandy land types on these 25 indicators were analyzed, and the importance of multiple indicator features was evaluated using Random Forest. Based on the distribution patterns of measured data in the effective factors, an integrated multi-indicator model was constructed.

The evaluation of different types of sandy land was determined based on measured data of vegetation cover and a sandy land type classification index (

Table 4. Index for sandy land evaluation.

Type of Sandy Land	Vegetation Cover Index
Shifting sand	Sandy terrain or dunes with vegetation cover less than 10%.
Semi-fixed sand	Vegetation cover is 10–30% on the uniformly distributed sandy terrain (dunes), inhibiting the movement of wind-blown sand.
Fixed sand	Sand dunes or sandy terrain with vegetation cover greater than 30%, where wind erosion is insignificant, and the surface is relatively stable.

). By using optical and radar data, multiple index products and texture features were generated. The numerical ranges of the samples among different product features were statistically analyzed to search for appropriate threshold values for classifying different types of sandy land.

The type of sandy land for the sample points was determined based on the measured data of vegetation cover and the index for sandy land evaluation. A total of 104 sample points were used, two-thirds (69 samples) for modeling and the remaining one-third (35 samples) for validation (

Table 5. Sample points of different types of sandy land.

Sample Points	Shifting Sandy Land	Semi-Fixed Sandy Land	Fixed Sandy Land	Total
Total	17	35	52	104
Modeling	11	23	35	69
Verification	6	12	17	35

).

3.3.1. Evaluation Parameter Indicator of Sandy Land Based on Optical Data

The parameters retrieved based on Landsat 8 OLI included Fraction of Vegetation Cover (FVC), Modified Soil-Adjusted Vegetation Index (MSAVI), Enhanced Vegetation Index (EVI), Surface Albedo, and Bare Soil Index (BSI). The numerical ranges of these parameters across different product datasets were statistically analyzed using a set of modeling sample points (69 in total). Threshold values were explored to differentiate various types of sandy land.

NDVI quantifies vegetation by measuring the difference between near-infrared and red light (Equation (4)):

$$NDVI={\displaystyle \frac{{\rho }_{NIR}-{\rho }_R}{{\rho }_{NIR}+{\rho }_R}}$$

(4)

where ${\rho }_R$ represents the reflectance in the red band and ${\rho }_{NIR}$ represents the reflectance in the near-infrared band. A higher NDVI value indicates better vegetation growth and lower degrees of sandy land.

FVC is estimated using a pixel-based binary model [42,43,44]. In areas with high degrees of sandy land, the vegetation cover is relatively low. FVC is calculated according to the following (5):

$$FVC={\displaystyle \frac{NDVI-{NDVI}_{min}}{{NDVI}_{max}-{NDVI}_{min}}}$$

(5)

FVC refers to Fraction of Vegetation Cover, ${NDVI}_{min}$ represents the minimum NDVI value for bare soil pixels, and ${NDVI}_{max}$ represents the maximum NDVI value for vegetation pixels. Due to various factors such as land surface conditions and vegetation types, the maximum and minimum values of NDVI can vary over time and space. In this study, the cumulative frequency of 5% was used as the minimum NDVI value, and 95% was used as the maximum NDVI value.

The MSAVI is capable of reducing the influence of soil background on vegetation pixels, thereby enhancing the sensitivity of vegetation classification in arid and semi-arid regions [45]. MSAVI is calculated as following (6):

$$MSAVI={\displaystyle \frac{\left(2{\rho }_{NIR}+1-\sqrt{{\left(2{\rho }_{NIR}+1\right)}^2-8\left({\rho }_{NIR}-{\rho }_R\right)}\right)}{2}}$$

(6)

MSAVI is a vegetation index that takes into account the influence of soil background, where ${\rho }_R$ represents the reflectance in the red band and ${\rho }_{NIR}$ represents the reflectance in the near-infrared band.

EVI is an enhanced vegetation index that incorporates the blue band, which helps to mitigate the atmospheric effects. It is designed to address the limitations of NDVI in being less sensitive to high vegetation cover and being influenced by soil background in areas with low vegetation cover [46]. EVI is calculated as following (7):

$$EVI=G\left[{\displaystyle \frac{{\rho }_{NIR}-{\rho }_R}{{\rho }_{NIR}+C_1{\rho }_R+C_2{\rho }_B+L}}\right]$$

(7)

where ${\rho }_B$ represents the reflectance in the blue band, ${\rho }_R$ represents the reflectance in the red band, ${\rho }_{NIR}$ represents the reflectance in the near-infrared band, C1 represents the atmospheric correction factor for the red band, C2 represents the atmospheric correction factor for the blue band, and L represents the canopy background correction factor. The empirical values for these parameters are typically set as follows: G = 2.5, C1 = 6, C2 = 7.5, and L = 1.0.

BSI is an index that combines soil background information and vegetation information [47,48,49]. It is well-suited for classifying areas ranging from low vegetation cover to high vegetation cover. BSI has strong applicability in arid and desertification-prone regions and is calculated as following (8):

$$BSI={\displaystyle \frac{\left({\rho }_{SWIR}+{\rho }_R\right)-\left({\rho }_{NIR}+{\rho }_B\right)}{\left({\rho }_{SWIR}+{\rho }_R\right)+\left({\rho }_{NIR}+{\rho }_B\right)}}$$

(8)

where ${\rho }_B$ represents the reflectance in the blue band, ${\rho }_R$ represents the reflectance in the red band, ${\rho }_{NIR}$ represents the reflectance in the near-infrared band, and ${\rho }_{SWIR}$ represents the reflectance in the shortwave infrared 1 band.

Albedo is used in the monitoring of sandy land to indicate surface temperature as well as changes in dryness and humidity [50,51,52]. A higher albedo value indicates a larger area coverage of bare soil, indicating a higher degree of sandy land and more severe vegetation degradation. Albedo is calculated as following (9):

$$Albedo=0.356{\rho }_B+0.13{\rho }_R+0.373{\rho }_{NIR}+0.085{\rho }_{SWIR1}+0.072{\rho }_{SWIR2}-0.0018$$

(9)

where ${\rho }_B$ represents the reflectance in the blue band, ${\rho }_R$ represents the reflectance in the red band, ${\rho }_{NIR}$ represents the reflectance in the near-infrared band, ${\rho }_{SWIR1}$ represents the reflectance in the shortwave infrared band 1, and ${\rho }_{SWIR2}$ represents the reflectance in the shortwave infrared band 2.

LST (Land Surface Temperature) is the temperature of the Earth’s surface. Due to varying vegetation cover, different types of land desertification exhibit variations in surface temperature. Higher temperatures indicate lower vegetation cover in sandy areas. There are typically three algorithms used to estimate LST using thermal infrared information: atmospheric correction, single-window algorithm, and split-window algorithm. In this study, the single-window algorithm was employed to estimate LST. The specific calculation formula is as following (10):

$$T_s=\left[a\left(1-C-D\right)+\left(b\left(1-C-D\right)+C+D\right)T_6-D{T}_a\right]/C$$

(10)

In the equation, $T_s$ represents the actual land surface temperature (K), a and b are constants with values of −67.355351 and 0.458606, respectively. C and D are intermediate variables; C = ετ and D = (1 $-$ τ) ([1 + (1 $-$ ε) τ]. Here, ε represents the surface emissivity and τ represents the atmospheric transmittance. $T_6$ represents the brightness temperature (K) measured by the satellite sensor at a certain height, and $T_a$ represents the average atmospheric temperature.

3.3.2. Evaluation Indicators of Sandy Land Based on Radar Data

Radar data possess a certain level of penetration capability, and the extraction of its texture features is of great significance in the identification of sandy land under vegetation cover. In this study, statistical properties were computed based on the GLCM (Gray-Level Co-occurrence Matrix) to quantitatively describe the texture features. These properties include mean, homogeneity, entropy, energy, dissimilarity, contrast, and correlation. The specific calculation formulas are shown in

Table 6. Texture feature calculation formula.

Texture Feature	Formula
Mean	$Mean={\sum }_i{\sum }_{j}p\left(i,j\right)\ast i$
Homogeneity	$Homogeneity={\sum }_i{\sum }_{j}p\left(i,j\right)\ast {\displaystyle \frac{1}{1+{(i-j)}^2}}$
Entropy	$Entropy={\sum }_i{\sum }_{j}p\left(i,j\right)\ast \mathrm{l}\mathrm{n}\mathrm{p}(\mathrm{i},\mathrm{j})$
Energy	$Energy={\sum }_i{\sum }_j{p\left(i,j\right)}^2$
Dissimilarity	$Dissimilarity={\sum }_i{\sum }_{j}p\left(i,j\right)\ast \left|i-j\right|$
Contrast	$Contrast={\sum }_i{\sum }_{j}p\left(i,j\right)\ast {(i-j)}^2$
Correlation	$Correlation={\sum }_i{\sum }_j{\displaystyle \frac{\left(i-Mean\right)\ast \left(j-Mean\right)\ast {p(i,j)}^2}{p\left(i,j\right)\ast {(i-Mean)}^2}}$

.

A sliding window with a size of 9 × 9 was used. The probability, $p\left(i,j\right)$ represents the probability of occurrence of a gray level j (column) given a starting point (row) i.

3.3.3. Construction of an Integrated Multi-Indicator Model for Sandy Land Evaluation

The importance of the 25 extracted indicator features mentioned above was assessed using random forest. Statistical analysis was conducted on the significant features to obtain reasonable thresholds for the evaluation of sandy land.

Random forest (RF) is known for its high predictive accuracy and strong tolerance towards outliers and noise. It effectively analyzes nonlinear data with collinearity and interaction effects. Additionally, it provides variable importance measures (VIM) while analyzing the data [53].

There are variables $X_1,X_2,\dots ,X_M$, for which VIM score statistics need to be calculated. VIM scores are computed based on the Gini index, and the score for variable $X_j$ is denoted as ${{VIM}_j}^{\left(Gini\right)}$.

The statistic ${{VIM}_j}^{\left(Gini\right)}$ represents the average change in node impurity for the jth variable across all trees in the random forest (RF). The calculation formula for the Gini index is as following(11):

$${GI}_m=\sum _{k=1}^k{\widehat{p}}_{mk}(1-{\widehat{p}}_{mk})$$

(11)

K is the number of classes in the bootstrap sample set, and ${\widehat{p}}_{mk}$ is the estimated probability of node m belonging to class k. The Gini index for node m in the case of binary classification data can be calculated as following(12):

$${GI}_m=2{\widehat{p}}_m(1-{\widehat{p}}_m)$$

(12)

${\widehat{p}}_m$ represents the estimated probability of a sample in node m belonging to any class. The importance of variable $X_j$ in node m, which is the change in the Gini index before and after the branching of node m, can be expressed as following(13):

$${{VIM}_j}^{\left(Gini\right)}={GI}_m-{GI}_l-{GI}_r$$

(13)

${GI}_l$ and ${GI}_r$ represent the Gini index of the two new nodes formed by the splitting of node m; the importance of variable $X_j$ in the ith tree, given that variable $X_j$ appears M times in that tree, can be calculated as following(14):

$${{VIM}_{ij}}^{\left(Gini\right)}=\sum _{m=1}^M{{VIM}_{jm}}^{\left(Gini\right)}$$

(14)

The Gini importance of variable $X_j$ in an RF is defined as following(15):

$${{VIM}_j}^{\left(Gini\right)}={\displaystyle \frac{1}{n}}\sum _{i=1}^n{{VIM}_{ij}}^{\left(Gini\right)}$$

(15)

where n represents the number of classification trees in the RF.

3.3.4. Accuracy Assessment Based on Field Data

The accuracy assessment method for the detection of sandy land information is as follows: Based on the location of the field survey plots, the detected results were classified as either sandy land or non-sandy land. Then, the extracted results were compared with the actual types determined through field surveys. The number of correctly classified plots was recorded, and the ratio of correctly classified plots to the total number of plots was calculated to evaluate the classification accuracy of this method. The accuracy assessment method for desertification land evaluation is as follows: Based on the vegetation cover data of the sandy land sampling plots obtained from field surveys, the sample plots were categorized into shifting sand, semi-fixed sand, and fixed sand. Then, these categories were compared with the type of sandy land assigned to the corresponding sampling plots in the sandy land evaluation results to determine if there were any misclassifications. The classification accuracy was calculated as the ratio between the number of correctly classified sampling plots and the total number of sampling plots.

4. Results

4.1. Detection of Sandy Land

To efficiently obtain the optimal method for detecting sandy land information, an experiment was conducted using multiple detection methods in selected representative areas. The experimental area was located in the southern part of Jinchang City, which had good image quality and abundant validation points. The land cover in this area was diverse, including water, vegetation, buildings, and various types of sandy land. Vegetation cover was analyzed as auxiliary information to analyze the results of sandy land detection, as shown in Figure 4. Landsat 8 OLI data imaged in August were chosen for vegetation cover inversion. There were a total of 37 sample points evenly distributed in the experimental area.

4.1.1. Methods of Sandy Land Detection

The six methods of sandy land detection in some areas of Gansu Province were: (1). based on Landsat 8 OLI data mixed pixel decomposition on Landsat 8 OLI data; (2). SVM classification upon polarization decomposition of Sentinel-1 data; (3). SVM classification using pixel-level fusion of Landsat 8 OLI and Sentinel-1 data (GS); (4). SVM classification using pixel-level fusion of Landsat 8 OLI and Sentinel-1 data (PCA); (5). SVM classification using pixel-level fusion of Landsat 8 OLI and Sentinel-1 data (HSV); (6). SVM classification using a feature-level fusion of Landsat 8 OLI and Sentinel-1 data. The results are shown in Figure 5. The SVM method was used to perform the classification of four land cover types: water, vegetation, buildings, and sandy land. Multiple sample data were uniformly selected for each land cover class. The separability of the region of interest was calculated using the Jeffries-Matusita (JM) distance. The JM distance was quantitatively determined to be greater than 1.9 to ensure good separability between the samples. This approach was employed to achieve land cover classification and detect sandy land. In the process of synergistically coupling multiple data sources, we analyzed the spectral curves of land cover obtained from Landsat 8 multispectral images (Figure 5b). It was found that the spectral information of sandy land exhibited a significant inflection point at the SWIR1 wavelength, showing a distinct difference compared to other typical land cover types. In the green and red bands, the reflectance of sandy areas showed an increasing trend and did not intersect with the spectral curves of other land cover types. In the blue and SWIR2 bands, the reflectance of sandy land was similar to that of buildings, which could lead to confusion between the two. In the NIR band, the spectral curve of vegetation exhibited a significant inflection point and was similar to the reflectance of sandy land, making it difficult to classify sandy areas from vegetation. Based on these observations, we selected the SWIR1, green, and red bands from optical data and synergistically coupled them with radar data to detect sandy land. Finally, the accuracy of the six detected sandy land information methods was assessed using field-measured sample data (

Table 7. The accuracy of sandy land detection based on field data.

Methods	Total Accuracy	Sandy Land Accuracy	Non-Sandy Land Accuracy
Sandy land detection based on Landsat 8	75.68%	61.90%	93.75%
Sandy land detection based on Sentinel-1	83.78%	95.24%	68.75%
Sandy land detection based on GS fusion	81.08%	80.95%	81.25%
Sandy land detection based on PCA fusion	75.68%	66.67%	87.50%
Sandy land detection based on HSV fusion	75.68%	71.43%	81.25%
Sandy land detection based on feature-level fusion	89.19%	90.47%	87.5%

). The method with the highest accuracy was selected to perform the detection of sandy land information in Gansu Province.

The overall accuracy of sandy land detection based on Landsat 8 multispectral images was relatively low, at 75.68%. The accuracy for sandy land specifically was as low as 61.90% (Table 7). It was found that by combining vegetation cover with single Landsat 8 OLI multispectral data (5-a), the detection performance was better for the shifting sandy land in the northern part of the study area. Sandy land with vegetation cover was not identified, resulting in lower accuracy, especially in the southern part of the study area, where identification of semi-fixed and fixed sand sandy was challenging. However, for non-sandy land, the identification accuracy was high at 93.75%. The multispectral data contains abundant spectral information, which allows for good recognition of surface features such as vegetation, buildings, and water. However, in the northern part of the study area, there was a misidentification of the road surface as sandy land. As a result, the method of detecting sandy land based on Landsat 8 multispectral image was found to perform well in identifying areas of bare sandy land. However, the underlying sandy land obscured by vegetation cover was difficult to be recognized. The overall accuracy of detecting sandy land based on single Sentinel-1 radar data was 83.78%. The accuracy for sandy land specifically was as high as 95.24%. This method not only detected bare sandy land but also identified sandy land under vegetation cover in the highly vegetated areas in the southern part of the study area. However, it also misclassified a large portion of non-sandy land as sandy land. For example, water in the eastern part and some vegetation areas in the southern part were identified as sandy land, resulting in a low accuracy of 68.75% for non-sandy land (Table 7). Compared to the extraction results of sandy land using Landsat 8 multispectral data, the lack of point-to-point neighborhood information and the irregular arrangement of points in radar data resulted in a weaker ability to recognize and classify land features. However, due to its advantage of canopy penetration, the accuracy of detecting sandy land was improved by 10.81%. This proved the importance of radar data in identifying sandy land under vegetation cover.

The accuracy of identifying non-sandy land information in the detection results using three pixel-level fusion methods was above 81%. This indicates that these methods correctly classified a higher proportion of non-sandy samples compared to the measured sandy samples. However, the accuracy of sandy land and the overall accuracy were relatively low. In particular, the accuracy of sandy land using PCA pixel-level fusion was as low as 66.67%. This means that the results of detecting sandy land based on pixel-level fusion images deviate from the ground truth information and cannot effectively identify sandy land. Additionally, there were image quality issues in the northwest corner of the pixel-level fusion data. The detection results of sandy land based on feature-level fusion were found to be consistent with the original image results. Compared to the results obtained from single multispectral data, the sandy soil obscured by vegetation cover could be effectively identified, and the road surface in the northeast part of the study area could be detected. In comparison to single radar data, it not only identified sandy land but also recognized non-sandy land. For instance, water in the eastern part was accurately monitored, and vegetation in the southern part was retained. When compared to the results of pixel-level fusion, the feature-level fused data exhibited better quality and accurately identified the boundaries of sandy land, facilitating the recognition of semi-fixed and fixed sand. Moreover, the extraction accuracy was improved, resulting in greater precision. It was observed from Table 7 that the feature-level fusion method, which was intended to take advantage of rich spectral features in optical data and the penetration capabilities of radar data, achieved an overall accuracy of 89.19%. The accuracy of sandy land reached as high as 90.47%, allowing for the recognition of both shifting sandy land and sandy land obscured by vegetation cover. The accuracy of non-sandy land was 87.50%, indicating that this method effectively highlighted the relevant multispectral information while eliminating or suppressing irrelevant information. Furthermore, it retained the backscatter information of targets in the SAR data and enabled the detection of vertical structural information of objects, thereby improving the classification results.

4.1.2. Detection and Verification of Sandy Land in Gansu Province

Based on the above research results, it can be concluded that the feature-level fusion method was the most effective in detecting sandy land. Therefore, this method was applied to the extraction of desertification information across Gansu Province (Figure 6). Gansu Province is situated in parts of the Tengger, Badain Jaran, and Kumtag Deserts, particularly in the northern region, including the cities of Jiuquan, Jiayuguan, Zhangye, Jinchang, and Wuwei, where the distribution of sandy land is more extensive. The southern part of Gansu Province is mainly characterized by grasslands and farmland, with relatively favorable ecological conditions. In particular, Gannan Tibetan Autonomous Prefecture is dominated by grasslands and forests. However, there were also areas with desert land in Lanzhou, Baiyin, and the northern regions of Qingyang city.

The accuracy of sandy land detection based on the proposed method was further validated using field survey data from 2021. Based on the locations of the field survey sample sites, the detected results were compared with the actual types determined during the field survey to see if a sample site was misclassified as sandy or non-sandy land. In 2021, a total of 314 sample sites were surveyed in Gansu Province, including 104 sandy land sample sites and 210 non-sandy land sample sites. The classification accuracy is shown in

Table 8. The accuracy of sandy land detection based on field data (Gansu).

Projects	Total Number	Sandy Land	Non-Sandy Land
Sample number	314	104	210
Correct number	271	93	178
Accuracy	86.31%	89.42%	84.76%

.

4.2. Integrated Multi-Indicator Model for Evaluating Sandy Land

The GEE platform was utilized to generate NDVI, MSAVI, FVC, EVI, Albedo, BSI, and LST products using Landsat 8 OLI data in Gansu Province. The SNAP software was used to process Sentinel-1 data and obtain the covariance matrix C₂ data. Texture features were extracted from C₁₁ and C₂₂, resulting in a set of 25 indicators as shown in Figure 7.

4.2.1. Multi-Indicator Analysis of Sand Characteristics Based on Optical Data

The numerical distribution patterns of the indicators obtained from the optical data inversion were analyzed for different types of desertification at 69 sample points. Various threshold values were explored to classify different types of desertification, and the results are shown in Figure 8.

The samples of different types of sandy land were found to have different distributions across various optical indicators. NDVI, MSAVI, FVC, and EVI, which are all vegetation indices, were found to have a higher degree of internal correlation. Therefore, samples of different types of sandy land exhibited a relatively consistent distribution across these indices. The sample points of shifting sand and fixed sand were observed to be clustered in the NDVI, MSAVI, FVC, and EVI indices, with a predominantly distributed range of around 0 to 0.3. Fixed sand, characterized by higher vegetation cover, exhibited higher values for the sample points with a broader distribution range. Due to the influence of vegetation cover in the surrounding areas of the sampling points, the vegetation cover at the satellite retrieval scale can be either underestimated or overestimated, depending on the sparsity or density of the surrounding vegetation. Additionally, the field measurements were conducted at times that did not coincide with the satellite image acquisition and variations in vegetation growth were observed in July and August due to precipitation. These factors introduce deviations when using ground-based measurements to validate satellite retrieval results. Consequently, it becomes challenging to select appropriate threshold values for sandy land assessment based on vegetation index products. According to the distribution of samples at different types of sandy land on the Albedo index, it can be observed that as the type of sandy land decreases, the values of the samples become lower, and the samples tend to cluster. However, the differences among samples at different types of sandy land were relatively small. Therefore, relying solely on this indicator was insufficient for sandy land assessment. The division of different types of sandy land based on the BSI index yielded poor results, as the sample points were clustered within different types with small variations. Analysis of the median, mean, and maximum temperature values within the study area indicated that as the type of sandy land decreased, the temperature decreased. This was due to higher vegetation cover associated with lower degrees of sandy land, resulting in lower temperatures. Although there was relatively high differentiation among different types of sandy land in the distribution of sample points based on the median indicator, it was still insufficient for evaluating the type of sandy land.

4.2.2. Multi-Indicator Analysis of Sand Characteristics Based on Radar Data

Statistical analysis was conducted on the numerical distribution patterns of different types of sandy land samples (69 in total) on radar data texture features. The aim was to explore the threshold values for classifying different types of sandy land. The results are shown in Figure 9 and Figure 10.

According to Figure 9, the distribution of different types of sandy land samples on the C11_Correlation texture feature was relatively observable. It showed a negative correlation with sandy land type, as the data points of different sandy land types converge near the mean value. In particular, the data points of the shifting sand were distributed around 0.97 without any outliers. However, their values were relatively close to the data points of semi-fixed sand, making it difficult to select a reasonable threshold to separate the two. The values of fixed sand were mostly distributed around 0.9, allowing some differentiation from shifting sand, but they cannot be separated from semi-fixed sand. Therefore, relying solely on the C11_Correlation feature was insufficient for sandy land assessment. The samples of different types of sandy land were found to exhibit a positive correlation in the C11_Dissimilarity texture feature. The values of shifting sand converged near the mean value and were primarily distributed between 0 and 2. Although the values increased with the weakening of desertification, both semi-fixed sand and fixed sand were distributed between 0 and 2. Therefore, an appropriate threshold for sandy land assessment could not be determined using this feature. There was no clear pattern or consistent distribution of different sandy land type samples in other texture features of C11. The distribution was more scattered, and these features could not be effectively used as decisive indicators for sandy land assessment.

Based on the distribution of samples from shifting sand, semi-fixed sand, and fixed sand in different texture features of C22 (as shown in Figure 10), it was observed that the texture features of C22 provided a greater type of differentiation among samples of different sandy land types compared to the texture features of C11. In particular, the C22_Correlation index exhibited a noticeable differentiation in the distribution range of samples from different sandy land types. The values of shifting sand tended to converge around the mean and were primarily distributed between 0.9 and 1. However, the values of semi-fixed sand were mainly distributed between 0.7 and 1, encompassing the distribution range of shifting sand, making it difficult to establish a suitable threshold to distinguish between the two. The values of fixed sand were primarily distributed between 0.4 and 0.8, and a threshold near 0.85 could be used to differentiate between shifting sand and fixed sand, but it was unable to accurately classify semi-fixed sand. Therefore, relying solely on this index was insufficient for effectively evaluating sandy land. The distribution of samples from different sandy land types in the C22_Dissimilarity texture feature showed a relatively improved separability compared to the distribution in C11_Dissimilarity. As the sandy land type decreased, the values increased. However, the values of both shifting sand and semi-fixed sand did not fall within the range of 0 to 5, making it difficult to select a suitable threshold to differentiate between these two types of sandy land. Therefore, using this index alone was insufficient for sandy land evaluation. The values of samples from different sandy land types in the C22_Homogeneity texture feature showed a negative correlation with sandy land type. Shifting sand tended to converge around the mean and was primarily distributed between 0.9 and 1.2, while fixed sand was mainly distributed between 0.2 and 0.9. Therefore, a threshold of 0.9 could be used to distinguish between shifting sand and fixed sand. However, distinguishing semi-fixed sand was challenging due to the high degree of dispersion in their distribution. Thus, using this index alone was insufficient for sandy land evaluation. The distribution of samples from different sandy land types in other texture features of C22 did not exhibit any clear patterns and showed a high degree of dispersion. Therefore, these features could not be effectively used as discriminative indicators for sandy land evaluation.

4.2.3. Construction of an Integrated Multi-Indicator Model for Evaluating Sandy Land

Based on the above results, it can be concluded that none of the single optical/radar data-based extracted indicator features can achieve sandy land evaluation effectively. Therefore, random forest was used to assess the importance of the 25 extracted indicators mentioned earlier. The 69 samples representing different types of sandy land were used as training samples to determine the importance of these 25 indicators in classifying different types of sandy land. The results are presented in

Table 9. The importance of each indicator.

Indicator	Weights	Indicator	Weights
C22_Correlation	0.062984	C22_Entropy	0.035900
FVC	0.058598	C22_Energy	0.035632
MSAVI	0.057848	C22	0.035098
NDVI	0 056949	C22_contrast	0.031950
C11	0.055700	C11_Correlation	0.030640
EVI	0.052310	C22_Dissimilarity	0.028150
C11_Mean	0.049878	C11_Energy	0.027598
Albedo	0.047819	C11_Entropy	0.026110
LST_Median	0.047159	C22_Homogeneity	0.025187
C22_Mean	0046510	C11_Homogeneity	0.023315
LST_Max	0.042202	C11_Dissimilarity	0.021639
LST_Mean	0.041436	C11_contrast	0.020442
BSI	0.038945

.

According to Table 9, among the 25 indicators, the correlation feature of C22 and the vegetation cover indicator showed a high level of importance. Therefore, a joint analysis and evaluation of these two indicators were conducted. Based on the distribution patterns of different sandy land type samples in C22_Correlation and FVC, it was found that using the ratio of these two indicators can facilitate the differentiation of sandy land types (as shown in Figure 11).

Based on Figure 11, it was observed that the distribution of samples from different sandy land types in the combined indicator C22_C/FVC showed a noticeable improvement in differentiation compared to using individual indices. The values of samples from shifting sand were primarily distributed within the range of 5 to 7, while samples from semi-fixed sand were mainly distributed between 3 and 5. Samples from fixed sand exhibited values mainly ranging from 0 to 2. However, there were still some outliers in the distribution of certain sample points, which can be attributed to the following reasons: (1) The ground truth data were collected within 1 m × 1 m sampling plots, while the remote sensing data used in this study had a spatial resolution of 30 m, leading to deviations in the validation results.

(2) The classification of different sandy land types of sample points was based on the assessment of vegetation cover data and sandy land indices obtained from field measurements. Upon individually examining the outliers, it was found that most of them had ground truth vegetation cover values near the threshold of the sandy land indices. The observation of ground truth vegetation cover was conducted using visual estimation or pin-pricking methods, which were prone to human errors during measurement and recording. Based on the above analysis, the region where the combined index C22_C/FVC was less than 2.2 was classified as fixed sand; the range between 2.2 and 5.2 represents semi-fixed sand, and values greater than 5.2 indicate shifting sand. This classification scheme was applied to the entire province of Gansu, and the results are presented in the following section.

4.2.4. Sandy Land Evaluation in Gansu Province

The spatial distribution of sandy land types is presented in Figure 12. Shifting sand accounts for 37.97% of the total sandy land, covering 85,100 square kilometers, mainly distributed in the northern parts of Jiuquan City and Wuwei City. Semi-fixed sand accounts for 38.91% of the total sandy land, covering 87,100 square kilometers, mainly distributed in the southern parts of Jiuquan City, Zhangye City, Lanzhou City, and Baiyin City. The fixed sand area accounts for 23.12% of the total sandy land area, covering 51,800 square kilometers, and is mainly distributed in the northern parts of Jiuquan City, the southern parts of Zhangye City, Jinchang City, and the southern parts of Wuwei City.

Based on the evaluation results of sandy land in Gansu Province, the accuracy assessment was conducted using the remaining one-third of the sample points (35 in total), and the results are shown in

Table 10. Accuracy assessment of sandy land evaluation.

Sample Points	Shifting Sand	Semi-Fixed Sand	Fixed Sand	Total
Number of plots	6	12	17	35
Number of correctly identified plots	5	9	14	28
Accuracy	83.3%	75.0%	82.4%	80.0%

. The overall accuracy was 80.0%, with the highest accuracy achieved for shifting sand at 83.3% and the lowest accuracy observed for semi-fixed sand at 75%. Most misclassified semi-fixed sand sample points had actual vegetation cover around 10% to 30%, being misidentified as either shifting sand or fixed sand.

5. Discussion

In response to the challenge of the true nature of sandy land obscured by vegetation cover, this study utilized methods based on both single-sensor and multi-source data fusion to detect sandy land. Based on the detection of sandy land using single Sentinel-1 radar data and feature-level fusion, the results are relatively good. However, radar data lacks neighborhood information between points, and the arrangement of points exhibits significant irregularity. Therefore, sandy land detection solely based on single Sentinel-1 radar data presents challenges in land feature identification and classification, leading to misclassification of non-sandy land as sandy land. The results revealed that the method based on feature-level fusion was found to be considerably more effective in terms of accuracy for identifying sandy land and non-sandy land. The distribution of bare sandy land was accurately detected, and for non-bare sandy land with vegetation cover, the method was able to decompose and quantify the impact of vegetation, thereby revealing the underlying sandy soil and achieving accurate detection. The Normalized Differential Sandy Areas Index (NDSAI) was proposed by Sahar et al. to separate the sandy areas from the non-sandy areas. Because SWIR1 is highly sensitive to vegetation, it was used to develop the new NDSAI. Therefore, NDSAI can better distinguish sandy land from vegetation [24]. In this article, when selecting the fusion bands based on optical data, the SWIR1 band was also chosen. The feature-level fusion algorithm focuses on the advantage of data to reveal the nature of sandy land in sparse vegetation scenes. Its accuracy was better than that of NDSAI, but the data processing was more complicated. Chen et al. proposed a remote sensing sand index called the Sand Differential Emissivity Index (SDEI), which can be conveniently used for large-scale detection of sandy areas. It addresses challenges related to extensive and efficient regional extraction as well as long-term sequence change analysis [54]. However, the results may not achieve superior accuracy in local areas. The algorithm proposed in this paper integrates radar data, making the process relatively complex. Nevertheless, it still attains considerable accuracy, particularly in local regions.

The issue of sparse vegetation in sandy land makes it difficult to objectively assess the type of sandy land using only the vegetation cover index. In this study, multiple optical features and radar texture indices were extracted to address this problem. The importance of these features was evaluated, and the distribution of field-measured data across multiple indicators was analyzed. A high-accuracy integrated multi-indicator model for sandy land evaluation was proposed and applied to Gansu Province. Compared to the land desertification assessment system established by Lamchin et al. using three indicators—NDVI, Albedo, and TGSI—and the method proposed by Wang, which integrates FVC, EVI, and MSAVI indices [37,39], the approach presented in this paper takes into account both terrain texture and vegetation characteristics. This consideration reduces the inherent correlations among multiple features, allowing the expression of information from various perspectives and thereby enhancing the classification accuracy for different types of sandy land. Furthermore, the model has been extensively validated with a fairly large number of field-measured samples, ensuring its reliability and authenticity. However, the discrepancies between estimated vegetation cover and ground measurements were somewhat significant. Ground measurement was made in a 1 m × 1 m sampling area, while the spatial resolution of the data in this study is 30 m. Due to the influence of vegetation cover in the surrounding areas of a sampling point (changes in vegetation density within the adjacent areas may result in changes in vegetation cover inversion at the satellite scale), bias manifested in the verification of the satellite inversion results based on ground-measured data. In addition, the sampling methods of ground-measured data, such as the visual estimation method or acupuncture method, cannot obtain information regarding vegetation structural distribution and type of spatial aggregation. The vegetation cover retrieved by remote sensing contained these two pieces of information, which was one of the main reasons for the discrepancies between ground measurement and satellite inversion results.

In the past, quantitative detection and evaluation of sandy land mainly relied on traditional methods such as on-site mapping based on topographic maps. Although these methods can achieve high accuracy in surveys, they require a significant amount of manpower and resources and are inefficient for large-scale land identification and classification. The feature-level fusion method coupled with an integrated multi-indicator model for sandy land evaluation based on remote sensing was proposed in this study. It efficiently achieved large-scale detection and evaluation of sandy land and obtained objective accuracy verification results by integrating ground survey data. Compared to manually designed methods for sandy land evaluation, the approach proposed in this paper, which utilizes measured data to automatically learn important features, is more scientifically sound and accurate. This research may provide a reference for the detection and evaluation of sandy land in arid and semi-arid regions of northern China in support of effective management and control of sandy land by local authorities.

This study has achieved satisfactory results on the regional scale, but it still needs a lot of effort to promote the national and even global scale. The method is limited to a certain extent, and the degree of automation is insufficient. In the future, we will further optimize the methods and procedures to extend the method to larger spatial scales, contributing to the SDG goal of Land Degradation Neutrality.

6. Conclusions

The research was conducted using optical and radar data as data sources, combined with ground truth sample point data, to investigate methods for sandy land detection. An integrated multi-indicator model for evaluating sandy land was constructed to evaluate sandy land over a large area.

A variety of methods were conducted to detect sandy land, and it was verified based on field sample data that the feature-level fusion method performed the best. This method was applied in Gansu Province, resulting in the large-scale detection of sandy land with an overall accuracy of 86.31%. The classification accuracy for sandy land reached 89.42%, while the classification accuracy for non-sandy land was 84.76%. Based on the detection results, it was found that Gansu Province, with a total area of 436,100 square kilometers, was predominantly affected by sandy land in its northern region, accounting for a high proportion of 51.36%, totaling 224,000 square kilometers. The assessment of the importance of the extracted optical data features and radar data texture information was conducted. It was found that the correlation feature of C22 and the vegetation cover indicator had higher importance. Based on the analysis of the distribution of the measured data in the joint index C22_C/FVC, a set of thresholds was determined. Using these thresholds, the sandy land distribution of Gansu Province was mapped, and statistics of each sandy land type were obtained.

The present study employed a feature-level fusion approach to detect sandy land in Gansu Province. This method effectively addressed the challenge of vegetation cover hindering the identification of sandy soils. Furthermore, an integrated multi-indicator model for evaluating sandy land was constructed to evaluate sandy land. Based on the results, it is evident that Gansu Province needs to further strengthen ecological restoration and improve the quality of the ecological environment to achieve the goal of sandy land prevention and control. In subsequent research, it is planned to adopt a technical approach that combines large-scale and high-resolution remote sensing data with ground surveys. The goal is to automate the entire process and achieve more accurate and real-time extraction and evaluation of sandy land over a large area.