Feature Selection and Regression Models for Multisource Data-Based Soil Salinity Prediction: A Case Study of Minqin Oasis in Arid China

Abstract

(1) Monitoring salinized soil in saline–alkali land is essential, requiring regional-scale soil salinity inversion. This study aims to identify sensitive variables for predicting electrical conductivity (EC) in soil, focusing on effective feature selection methods. (2) The study systematically selects a feature subset from Sentinel-1 C SAR, Sentinel-2 MSI, and SRTM DEM data. Various feature selection methods (correlation analysis, LASSO, RFE, and GRA) are employed on 79 variables. Regression models using random forest regression (RF) and partial least squares regression (PLSR) algorithms are constructed and compared. (3) The results highlight the effectiveness of the RFE algorithm in reducing model complexity. The model incorporates significant environmental factors like soil moisture, topography, and soil texture, which play an important role in modeling. Combining the method with RF improved soil salinity prediction (R² = 0.71, RMSE = 1.47, RPD = 1.84). Overall, salinization in Minqin oasis soils was evident, especially in the unutilized land at the edge of the oasis. (4) Integrating data from different sources to construct characterization variables overcomes the limitations of a single data source. Variable selection is an effective means to address the redundancy of variable information, providing insights into feature engineering and variable selection for soil salinity estimation in arid and semi-arid regions.

1. Introduction

Widespread soil salinization, a form of desertification, poses a threat to food security and sustainable development, especially in arid and semi-arid regions [1]. Arid regions of Northwest China are severely affected by salinization due to arid climate, scarce precipitation, intense evaporation, and low soil water content [2,3], which hinders modern agricultural development and ecological security [4].

To address this global challenge, timely and cost-effective monitoring of soil salinization is crucial [5]. Traditional methods involving on-site sample collection and laboratory analysis are both time-consuming and expensive due to intensive sampling requirements [6,7,8]. Remote sensing technology, particularly with the Sentinel-2 satellite’s Multi-Spectral Instrument (MSI) launched in 2015, offers a valuable alternative [9,10,11,12,13]. Although it provides new bands (three red-edge bands and two SWIR bands), most scholars still prefer to use the traditional visible and near-infrared bands, and the sensitivity of the new red-edge bands to soil salinity information needs to be explored. However, despite its benefits, relying solely on optical remote sensing has limitations due to weather, vegetation, and soil salt crusts. Microwave remote sensing, sensitive to soil salinity through the soil’s imaginary dielectric constant, provides a solution by overcoming these constraints and enabling soil salinization inversion from radar data [14,15,16]. While there’s potential in combining optical and microwave remote sensing, such studies remain relatively scarce.

The formation and development of salinized soils are closely related to the climate, biology, topography, physical and chemical properties of soils, and other associated factors of the environment in which they are located. The introduction of environmental variables to participate in the establishment of inversion models has become a meaningful way to improve the quality of modeling [17]. Among them, the common vegetation indices, soil texture, soil moisture, slope, slope orientation, elevation, surface roughness, texture, etc., have also been critical environmental variables [1,18,19]. Since most of the environmental covariates are obtained by band math, the information has a certain degree of redundancy. Therefore, it is necessary to select the regionally optimal characterization variables through variable screening [20]. Pearson correlation analysis is the most common method for modeling the optimal screening of environmental variables [17,21]. However, the Pearson method ignores the nonlinear relationship between predictor and target variables. Wrapper methods integrate the feature selection process with the training process, and the predictive ability of the model is used as a measure of feature selection. Commonly, there are gray relational analysis (GRA) [19], Genetic Algorithms (GA) [22], Penalty Term Based Feature Selection Methods [23] (e.g., LASSO, Ridge regression, Elastic Net algorithm, etc.), random forest algorithm (RF), and Best Subset Selection (BSS) [24]. Feature engineering is an essential part of the modeling process. However, there are no established guidelines for the optimal feature set for arid zones. More attention needs to be paid to this area. Few scholars have taken visible spectral features, microwave spectral features, texture features, and terrain features into account simultaneously.

Machine learning and data mining techniques, such as random forest (RF), decision tree (Classification and Regression Trees), support vector machine (SVM), and artificial neural network (ANN), surpass traditional statistical methods in assessing soil salinity. These methods effectively address the nonlinear relationship between soil and environmental factors, offering improved precision [23,25,26,27,28].

Google Earth Engine (GEE) is a cloud-based platform that provides an integrated environment of massive data catalogs (Landsat 4–8, Sentinel 1–3, etc.) and thousands of computers for planetary-level data analysis with a friendly workbench environment that allows for interactive algorithm development via a JavaScript (or Python, R) based application programming interface (GEE API) for interactive algorithm development. A growing number of researchers are now using the platform for large-scale land cover and land use change mapping and monitoring [29,30]. Some scholars such as Ivushkin et al. [31] have also updated global soil salinity maps using soil property maps and thermal infrared images in the GEE cloud platform, thus further improving spatial accuracy. However, few studies have explored the potential of cloud computing platforms (e.g., GEE) in extracting feature variables and constructing soil salinity prediction models.

Consequently, this study employed optical remote sensing data, microwave remote sensing data, and digital terrain data in Minqin Oasis. A comprehensive extraction of 79 environmental variables, including sensitive wavebands, vegetation indices, salinity indices, terrain factors, and texture information, was conducted using the GEE platform. The model was also used to assess and map the spatial distribution of soil salinity in the lower part of the Shiyang River Basin in the arid zone of northwest China. The main objectives are (1) to explore the effectiveness of four feature selection methods to remove redundant information from variables, (2) to determine the optimal set of environmental variables applicable to the study area and explore the importance of the variables, and (3) to develop an optimal prediction model for soil salinity in the study area and map the spatial distribution. This study will provide valuable guidance for predicting soil salinity in arid zones.

2. Materials and Methods

In this study, the overall workflow includes data pre-processing, extraction and selection of feature variables, model construction, and evaluation (Figure 1). The specific research content is as follows.

2.1. Study Area and Field Sampling

The Shiyang River Basin is located east of the Hexi Corridor in Gansu Province, China. It is one of the three major inland river basins of the Hexi Corridor, a traditional agricultural area in China. Minqin Oasis is located in the lower reaches of the Shiyang River (38°03′ N to 39°28′ N, 101°49′ E to 104°12′ E), which is a typical arid inland river basin area [32]. The climate type of this region is arid desertification, with an average annual temperature of 7.8 °C, an average annual precipitation of about 116.5 mm, and an average annual evaporation of 2308 mm [33]. Minqin Oasis is bordered by the Tengger Desert in the east and the Badain Jaran Desert in the west (Figure 2), with severe soil desertification, accounting for 94.5% of the total area [34]. The Minqin oasis relies on the Hongyashan reservoir and irrigation canal system for water supply. Due to the reduction of surface water inflow into the oasis from the upper and middle reaches of the Shiyang River, groundwater has become the main source of irrigation water in order to meet water demand, leading to over-exploitation [35]. This triggered a rise in groundwater levels and an increase in dissolved mineral concentrations. The combination of saline irrigation and inadequate drainage measures led to rapid soil salinization [32]. Soil texture in the study area is mainly of five types: clay loam, loam, sandy loam, loamy sand, and sandy loam. Loamy sand is dominant, and cropland is predominantly loam (Figure 2b).

During the field survey conducted on 3–7 June and 26–30 July 2023, a total of 215 sample points were selected from different land use categories with obvious salinization in the study area; soil conductivity (uS∙cm⁻¹), soil temperature, soil moisture, soil pH, elevation value, and land use types in the soil surface layer (0–10 cm) were collected in each sample point using a soil environment detector (FK-WSYP, Fangke, China), an RTK(Huace, China) device, and an NDVI measuring instrument, and the spatial distribution of sample points is shown in Figure 2c. The conductivity was converted from uS∙cm⁻¹ to the standard unit dS∙m⁻¹.

2.2. Data

2.2.1. Remote Sensing Data and Pre-Processing

In this study, we mainly used Sentinel-1 and Sentinel-2 imaging data, and SRTM DEM data provided by the GEE platform. GEE has been providing atmospherically and geometrically corrected level 2A data since 28 March 2017. In this paper, we selected data with cloudiness within 10% during the two sampling periods (3–7 June and 26–30 July). In order to avoid the influence of clouds and shadows on the model as much as possible, the median synthesis method was utilized to reconstruct the Sentinel-2 minimum cloudiness composite image after removing clouds using s2cloudless. The four bands B2, B3, B4, and B8 with a resolution of 10 m and the three red-edge bands B5, B6, and B7 with a resolution of 20 m were selected for the calculation of the relevant feature variables.

In view of the fact that Sentinel-2 data are affected by cloudy weather, GEE’s Sentinel-1A GRD level 1 data, which has VV and VH polarization modes, is selected for this study. These data were pre-processed by orbit file application, noise removal, radiometric calibration, and terrain correction. The study area image was selected in synchronization with the Sentinel-2 data, and the backscattering coefficients for the VV and VH polarization modes were derived using the median() function and the Refined Lee filter to remove coherent speckle.

Digital elevation model (DEM) is a discrete mathematical representation of surface topography that reflects topographic changes. In this paper, the SRTM DEM 30 m resolution data are used to extract the topographic factors to provide an analytical basis for the study.

2.2.2. Environmental Covariates

• Variables Based on Microwave Data;

Although the study area is located in the arid zone, the surface of the oasis in the region was still covered by a large amount of vegetation during the sampling period, which had a noticeable effect on obtaining the backscattering coefficient of the soil and, thus, the interference of the vegetation layer on the scattering needed to be eliminated. The water cloud model (WCM) simplifies the scattering process of microwave remote sensing by assuming that the vegetation layer is a horizontally homogeneous cloud layer and ignoring the scattering effect between the vegetation layer and the soil layer below the layer [36,37]. The expression is as follows:

$${\mathsf{\sigma }}_{\mathrm{pp}}^0={\mathsf{\sigma }}_{\mathrm{veg}}^0+{\mathsf{\gamma }}^2{\mathsf{\sigma }}_{\mathrm{soil}}^0$$

(1)

$${\mathsf{\sigma }}_{\mathrm{veg}}^0=\mathrm{A}\times {\mathrm{m}}_{\mathrm{veg}}\times \cos \left(\mathsf{\theta }\right)\times \left(1-{\mathsf{\gamma }}^2\right)$$

(2)

$${\mathsf{\gamma }}^2={\mathrm{e}}^{-2\mathrm{B}\times {\mathrm{m}}_{\mathrm{veg}}\times \sec \left(\mathsf{\theta }\right)}$$

(3)

where σ_pp⁰ is the total backscattering coefficient; σ_veg⁰ is the backscattering coefficient of the vegetation; σ_soil⁰ is the backscattering coefficient of the soil surface; γ² is the attenuation factor for the penetration into the vegetation layer; m_veg is the moisture content of vegetation; θ is the angle of incidence of radar waves; and A and B are the vegetation parameters.

The results of the backscattering coefficients processed by the WCM model are shown in Figure 3, and it can be seen that the backscattering coefficients of both polarization modes are reduced. In this study, the backscattering coefficients of the two polarization modes (VV, VH) provided by Sentinel-1A were extracted by the WCM model after removing the influence of vegetation as the base data. Six polarization combinations were constructed as the eigenvariables: WCMVV, WCMVH, (WCMVV − WCMVH)/(WCMVV + WCMVH)(Com1), WCMVV + WCMVH(Com2), WCMVV − WCMVH(Com3), and WCMVV/WCMVH(Com4) [37,38].

2.

Variables based on optical remote sensing data;

The correlation between soil physicochemical properties and spectral indices is significant. Therefore, multiple indexes are constructed through band operation, and indexes better related to soil salinity are selected to evaluate the degree of soil salinization. This study set 33 spectral indices and transformed bands based on existing research results in arid regions, as shown in

Table 1. Spectral indices and transformed band lists.

	Index	Formula	References
Salinity Index	Salinity Index (SI-T)	$\left(\mathrm{R}∕\mathrm{NIR}\right)\times 100$	[39]
	Brightness Index (BI)	$\sqrt{{\mathrm{R}}^2+\mathrm{NI}{\mathrm{R}}^2}$	[39]
	Salinity Index (SI)	$\sqrt{\mathrm{B}\times \mathrm{R}}$	[40]
	Salinity Index 1 (SI1)	$\sqrt{\mathrm{R}\times \mathrm{G}}$	[40]
	Salinity Index 2 (SI2)	$\sqrt{{\mathrm{G}}^2+{\mathrm{R}}^2+{\mathrm{NIR}}^2}$	[41]
	Salinity Index 3 (SI3)	$\sqrt{{\mathrm{G}}^2+{\mathrm{R}}^2}$	[41]
	Salinity Index (S1)	B/R	[41]
	Salinity Index (S2)	(B − R)/(B + R)	[42]
	Salinity Index (S3)	(G$\times$R)/B	[42]
	Salinity Index (S5)	(B$\times$R)/G	[42]
	Salinity Index (S6)	(R$\times$NIR)/G	[42]
	Normalized Difference Salinity Index (NDSI)	(R − NIR)/(R + NIR)	[40]
Vegetation Index	Normalized Difference Vegetation Index (NDVI)	(NIR − R)/(NIR + R)	[43]
	Canopy Response Salinity Index (CRSI)	$\sqrt{\frac{\left(\mathrm{NIR}\times \mathrm{R}\right)-\left(\mathrm{G}\times \mathrm{R}\right)}{\left(\mathrm{NIR}\times \mathrm{R}\right)+\left(\mathrm{G}\times \mathrm{R}\right)}}$	[44]
	Enhanced Vegetation Index (EVI)	$2.5\times \left(\frac{\mathrm{NIR}-\mathrm{R}}{\mathrm{NIR}+6\mathrm{R}-7.5\mathrm{B}+1}\right)$	[45]
	Difference Vegetation Index (DVI)	$\mathrm{N}\mathrm{I}\mathrm{R}-\mathrm{R}$	[43]
	Soil Adjusted Vegetation Index (SAVI)	$(\mathrm{NIR}-\mathrm{R})/(\mathrm{NIR}+\mathrm{R}+\mathrm{L})$	[46]
	Modified Soil Adjusted Vegetation Index (MSAVI)	$(2\mathrm{NIR}+1-\sqrt{{\left(2\mathrm{NIR}+1\right)}^2-8\left(\mathrm{NIR}-\mathrm{R}\right) }$)/2	[47]
	Salinity Remote Sensing Index (SRSI)	$\sqrt{{\left(\mathrm{NDVI}-1\right)}^2+{\mathrm{SI}}^2}$	[48]
	Ratio Vegetation Index (RVI)	$\mathrm{NIR}/\mathrm{R}$	[49]
	Generalized Difference Vegetation Index (GDVI)	$\frac{{\mathrm{NIR}}^2-{\mathrm{R}}^2}{{\mathrm{NIR}}^2+{\mathrm{R}}^2}$	[50]
Spectral Transformation	Principal Component Analysis (PCA)	PC1, PC2, PC3	[51]
	Tasseled Cap (TC)	TC1, TC2, TC3	[52]
	Sentinel-2 MSI Bands	B4, B3, B2, B8
Red-edge Spectral Index	Red-edge Salinity Index S1 (RES1)	B/$\mathrm{RedEdge}$	[19]
	Red-edge Salinity Index S2 (RES2)	(B − $\mathrm{RedEdge}$)/(B + $\mathrm{RedEdge}$)	[19]
	Red-edge Salinity Index S3 (RES3)	(G$\times \mathrm{RedEdge}$)/B	[19]
	Red-edge Salinity Index S5 (RES5)	(B$\times \mathrm{RedEdge}$)/G	[19]
	Red-edge Salinity Index S6 (RES6)	(RedEdge$\times$NIR)/G	[19]
	Red-edge Salinity Index SI (RESI)	$\sqrt{(\mathrm{B}\times \mathrm{RedEdge})}$	[19]
	Red-edge Salinity Index S1 (RESI1)	$\sqrt{(\mathrm{G}\times \mathrm{RedEdge})}$	[19]
	Red-edge Salinity Index SI2 (RESI2)	$\sqrt{{\mathrm{G}}^2+{\mathrm{RedEdge} }^2+{\mathrm{NIR}}^2}$	[19]
	Red-edge Salinity Index SI3 (RESI3)	$\sqrt{{\mathrm{G}}^2+{\mathrm{RedEdge} }^2}$	[19]

. In this study, principal component analysis (PCA) and tasseled cap transformation (TC) were utilized to reduce the dimensionality of the dataset and extract relevant features. PCA is a widely used statistical technique that transforms the original variables into a new set of uncorrelated variables called principal components, which capture the maximum variance in the data. This helps in simplifying the dataset while retaining essential information. Similarly, the TC transformation is specifically designed for multispectral remote sensing data and generates new indices (Brightness, Greenness, Wetness) that summarize the information content related to soil and vegetation properties. Similar to the use of salinity indices, vegetation indices, and microwave backscattering coefficients, PCA and TC transformations provide crucial insights into the soil and environmental conditions [21,24].

3.

Texture feature extraction;

The texture is the spatial variation and repetition of the grayscale of the remote sensing image. The texture feature parameter is a vital feature parameter reflecting the spatial distribution of the grayscale pattern of remote sensing images. The grayscale matrix method is the most mature texture extraction method, which can describe the spatial distribution and structural characteristics of each grayscale pixel of the image and is widely used to extract texture feature values of remote sensing images [23]. Although 18 GLCM outputs are provided in GEE, in order to avoid dimensional disaster, five indicators with high correlation with soil salinity, namely Correlation (COR), Contrast (CON), Dissimilarity (DIS), Entropy (ENT), and Angular Second Moment (ASM), were selected as texture features based on preliminary experiments and existing research references [53,54]. And the B2, B8, and gray value bands of the Sentinel-2 were used as the computational bands.

Table 2. Texture characteristics chosen for the article and the formula for calculating them.

Feature Name	Formula	Clarification
Contrast (CON)	${\displaystyle \sum _{\mathrm{i},\mathrm{j}=0}^{\mathrm{n}-1}}\frac{{\mathrm{P}}_{\mathrm{i}-\mathrm{j}}}{1+{(\mathrm{i}-\mathrm{j})}^2}$	Describes the clarity of the remotely sensed image and how it produces changes in a small area.
Dissimilarity (DIS)	${\displaystyle \sum _{\mathrm{i},\mathrm{j}=0}^{\mathrm{n}-1}}{\mathrm{P}}_{\mathrm{i},\mathrm{j}}\left|\mathrm{i}-\mathrm{j}\right|$	Describes the localized contrast of a remotely sensed image.
Entropy (ENT)	$\sqrt{{\displaystyle \sum _{\mathrm{i},\mathrm{j}=0}^{\mathrm{n}-1}}{\mathrm{P}}_{\mathrm{i}.\mathrm{j}}(-\ln {\mathrm{P}}_{\mathrm{i},\mathrm{j}})}$	Indicators of texture feature clutter in remotely sensed images.
Angular Second Moment (ASM)	${\displaystyle \sum _{\mathrm{i},\mathrm{j}=0}^{\mathrm{n}-1}}{{\mathrm{P}}^2}_{\mathrm{i},\mathrm{j}}$	Homogeneity and consistency of grayscale distribution of remotely sensed images.
Correlation (COR)	${\displaystyle \sum _{\mathrm{i},\mathrm{j}=0}^{\mathrm{n}-1}}\frac{{\mathrm{P}}_{\mathrm{i},\mathrm{j}}(\mathrm{i}-{\mathrm{u}}_{\mathrm{i}}\left)\right(\mathrm{j}-{\mathrm{u}}_{\mathrm{i}})}{\sqrt{\left({\mathsf{\sigma }}_{\mathrm{i}}^2\left)\right({\mathsf{\sigma }}_{\mathrm{j}}^2\right)}}$	Describes the directionality of linear targets in remotely sensed images.

shows the formula and descriptions of the texture feature computation selected in this paper. The texture feature calculation function glcmTexture() is provided in GEE, in which the default 3 × 3 window size is used for calculating GLCM.

4.

Topographic factor extraction.

DEM-derived topographic variables directly or indirectly determine the direction of soil water movement, altering the pattern and location of soil salt accumulation [55]. Studies have confirmed the relationship between topography and soil salinity [43]. In this paper, four feature variables, elevation, slope, aspect, and surface roughness, were extracted based on DEM and resampled to the exact resolution as the Sentinel data.

2.3. Feature Selection Algorithms

Variable screening reduces prediction error by removing potentially irrelevant predictor variables. In this study, Pearson correlation analysis, Spearman correlation analysis, and Kendall correlation analysis were used for initial variable selection. In statistics, all three correlation coefficients are used to measure the correlation between variables and take values between −1 and 1. Pearson’s correlation coefficient is used for continuous variables that are typically distributed. The Spearman’s correlation coefficient is similar, except that the Spearman’s calculation requires that the two variables be transformed into ordinal numbers and are independent of outliers, which makes it suitable for nonlinear situations. Kendall’s correlation coefficient measures the strength of monotonicity between the two ordinal variables. Eliminating highly correlated and redundant variables at an early stage can prevent the subsequent feature selection methods from being biased by multicollinearity [56,57,58].

The Elastic Net algorithm, the GRA algorithm, and the RFE algorithm were used to minimize the correlation between variables further. To reduce the risk of overfitting and multicollinearity and to provide interpretable core models, shrinkage methods that regularize the estimated coefficients by reducing variance and residuals have been widely used [59]. Elastic Net combines Lasso and Ridge’s advantages and controls the penalty term’s size by two parameters, λ and ρ. It retains the nature of the feature selection of Lasso while considering the stability of Ridge regression.

The Recursive Feature Elimination (RFE) algorithm works by adding or removing specific feature variables to obtain the optimal combination of variables that maximizes model performance [60]. In this study, the RFE algorithm is implemented using the caret package of the R language(4.3.0), in which four models, namely, random forests (rfFuncs), linear regression (lmFuncs), support vector machines (SvmFuncs), and bagging decision trees (treebagFuncs), are selected for variable ranking, and the optimal method is chosen based on RSME.

Gray correlation analysis (GRA) is a multifactorial statistical analysis method. The basic idea is to determine whether the serial curves are closely linked based on the degree of similarity in their geometric shapes [61]. The closer the curves are, the greater the correlation between the corresponding sequences is, and vice versa, the smaller it is [62]. The formula is as follows:

$${\mathrm{r}}_{\mathrm{j}}=\frac{1}{\mathrm{n}}{\displaystyle \sum _{\mathrm{k}=1}^{\mathrm{n}}}{\mathsf{\xi }}_{\mathrm{j}}\left(\mathrm{k}\right)={\mathrm{z}}_{\mathrm{k},\mathrm{j}}=\frac{\mathrm{a}+\mathsf{\rho }\mathrm{b}}{\left|{\mathrm{x}}_{\mathrm{kj}}-{\mathrm{y}}_{\mathrm{j}}\right|+\mathsf{\rho }\mathrm{b}}$$

(4)

where ρ is the resolution coefficient, generally taken as 0.5, and a and b are the bipolar minimum and maximum differences, respectively.

2.4. Modeling Strategy

2.4.1. Modeling Algorithm

Partial least squares method (PLSR) is a classical multivariate statistical data analysis method, which integrates the advantages of modeling methods such as MLR and PCR, and can well extract and reflect the information on data variability, and at the same time improve the explanatory power of the best composite variables by decomposing the information about the data, and it has been widely used in salinity inversion in arid zones [11,63,64].

The random forest model (RF) is an integrated learning method based on decision trees, where multiple random samples are obtained through bootstrap sampling, and corresponding decision trees are built from these samples to form a random forest [65]. RF-based regression modeling has achieved better inversion accuracy in soil salinity inversion and has been the most used algorithm in recent years [66].

2.4.2. Model Evaluation

Hyperparameters for each algorithm were determined using a grid search algorithm with 3-fold cross-validation. Model accuracy was assessed by calculating the coefficient of determination (R²), root mean square error (RMSE), and relative analysis error (RPD) [67]. Reliable models are usually characterized by high R² and RPD values, and low RMSE values. The specific formulas are as follows:

$${\mathrm{R}}^2=\frac{\mathrm{S}\mathrm{S}\mathrm{R}}{\mathrm{S}\mathrm{S}\mathrm{T}}=\frac{{\displaystyle {\sum }_{\mathrm{i}=1}^{\mathrm{n}}}{(\widehat{{\mathrm{y}}_{\mathrm{i}}}-\overline{\mathrm{y}})}^2}{{\displaystyle {\sum }_{\mathrm{i}=1}^{\mathrm{n}}}{({\mathrm{y}}_{\mathrm{i}}-\overline{\mathrm{y}})}^2}$$

(5)

$$\mathrm{MSE}=\sqrt{\frac{1}{\mathrm{m}}{\displaystyle \sum _{\mathrm{i}=1}^{\mathrm{m}}}{({\mathrm{P}}_{\mathrm{i}}-{\mathrm{M}}_{\mathrm{i}})}^2}$$

(6)

$$\mathrm{RPD}=\frac{{\mathrm{S}}_{\mathrm{y}}}{\mathrm{RMSE}}$$

(7)

The overall error of each algorithm was first evaluated, and then the predictions were further validated using salinity classes. Salinity classes were categorized into five classes based on Richards et al. [68]: non-saline (0–2 dS∙m⁻¹), slightly saline (2–4 dS∙m⁻¹), moderately saline (4–8 dS∙m⁻¹), strongly saline (8–16 dS∙m⁻¹), and extremely saline soils (>16 dS∙m⁻¹).

3. Results and Discussion

3.1. Sample Descriptive Statistics

We collected 215 samples with conductivity values ranging 0–22.19 dS∙m⁻¹, mean 2.1026 dS∙m⁻¹, and standard deviation 3.83538 dS∙m⁻¹ (coefficient of variation 1.824). These were divided 8:2 into training and validation sets. Figure 4a shows datasets with high coefficients of variation, indicating significant soil salinity variability. Figure 4b shows the box line diagrams of the distribution of sample points under different land use types. Overall, the sample points with EC values above 8 dS∙m⁻¹ are mainly distributed on unutilized land, which belongs to the heavily and extremely salinized soils. The cultivated land type has the second largest number of sample points, which belongs to non-salinized or mildly salinized soil. However, there are also a certain number of sample points with high EC values, indicating the existence of secondary salinization problems. Woodland, mainly pike and poplar forests, reached 4.56 dS∙m⁻¹, suggesting surface salt accumulation due to evapotranspiration, low precipitation, and high groundwater mineralization. Grasslands in the study area are mainly low-cover grasslands lacking moisture and with severe desertification, and the points with high EC values are near the salt ponds, aligning with findings in the Minqin oasis, where soil salinization is widespread, with saline wastelands having the highest salinity at oasis edges [32].

3.2. Feature Selection

3.2.1. Correlation Analysis

Figure 5a represents three correlation coefficients for variables with EC at p < 0.01 and p < 0.05. Figure 5b shows a heat map of Pearson’s correlation coefficients after excluding the variables with no significant correlations. Coefficients ranged from −0.39 to 0.35, indicating significant correlations (p < 0.05) between most environmental covariates and soil conductivity, except for salinity indices S1, S2, S3, S4, S5, SI, SI1, and SI3. Notably, salinity indices using red-edge bands exhibited superior correlation, likely due to enhanced spectral information and a higher signal-to-noise ratio. Among all 79 variables, pc3, the third principal component of the principal component analysis, had the highest Pearson correlation coefficient (r = 0.39), Spearman correlation coefficient of RE3SI1 (r = −0.26), and Kendall correlation coefficient of B8_asm (r = −0.18), whereas RE3SI1 and B8_asm did not have significant linear correlation, and there may be nonlinear correlation. Based on the results of three correlation analyses, 34 significant and relatively highly correlated variables were initially selected from 79 environmental covariates. They are B2, B2_contrast, B8, BI, Brightness, CRSI, Com3, Com4, DVI, EVI, GDVI, Greenness, MSAVI, NDSI, NDVI, RE1S6, RE2S1, RE2S2, RE3S3, RE3S5, RE3SI3, RVI SAVI, SI2, SI_T, WCMVH, Wetness, elevation, gray_contrast, gray_diss, pc2, pc3, RE3SI1, and RE2SI. Figure 5b revealed a high degree of correlation (r > 0.9) among the selected variables, necessitating further screening to mitigate information redundancy and multicollinearity.

3.2.2. Feature Variable Combination Selection

The 34 variables underwent GRA, Elastic Net, and RFE screening. Prior to algorithm application, all variable values underwent feature scaling for magnitude and value range standardization. For RFE, four models were tested, and Figure 6a,b display the RMSE and R² for each model under various feature subsets. In both RMSE and R², the RF algorithm outperformed, reaching its best at 11 features. Consequently, the RFE algorithm of the RF model was selected for feature selection. The Elastic Net algorithm employed ρ = 0.5 for L1 and L2 regularization, with grid search determining the optimal λ value of 0.07. In the GRA algorithm, rho = 0.5 regulated correlation coefficient differences. The GRA algorithm ranked variable importance; features with importance over 0.7 were selected.

Figure 6c–e shows the feature and importance ordering after feature selection using the Elastic Net, GRA, and RFE-RF algorithms, respectively. It can be seen that both Elastic Net and RFE-RF are able to distinguish the importance of variables better, and there is little difference in the importance of most of the variables under the GRA algorithm. Figure 6f shows a schematic diagram of the variables selected by the three algorithms. Notably, Com4, SI2, WCMVH, and Wetness appear in all three results, with Wetness, SI2, and Com4 exhibiting considerable importance. These three variables represent exactly three data sources: the tassel-capped transformed Wetness component, the salinity index, and the microwave backscattering coefficient polarization combination, respectively. It can be seen that, in addition to the traditional optical remote sensing variables, environmental covariates as well as microwave data also have a high sensitivity to soil salinity and are a factor to be considered when modeling soil salinity prediction.

3.3. Modeling Assessment

As the GRA algorithm lacks the capability to provide the optimal subset of feature variables with known importance ordering, this paper employs a circular, iterative approach. Initially, all variables enter PLSR or RF modeling, undergoing accuracy tests (R², RMSE). Subsequently, the last variable in the importance ordering is eliminated, and the remaining variables re-enter the modeling operation. This process continues iteratively, with each cycle involving the removal of the last variable until only two variables remain in the operation. R² and RMSE values are generated for each step, allowing the determination of the optimal number of variables. For optimal model results, a grid search with three-fold cross-validation is employed to select parameters. PLSR focuses on determining the principal components (n_components), while RF involves parameters such as n_estimators, max_depth, min_samples_leaf, min_samples_split, and max_features. Default values are set for other parameters.

Table 3. The optimal parameters for each algorithm and the number of feature variables used in the model.

Algorithms	Parameters	Number of Features
GRA-PLSR	n_components = 3	10
Elastic-PLSR	n_components = 8	10
RFE-PLSR	n_components = 2	11
GRA-RF	n_estimators = 27, max_depth = 4,min_samples_split = 10,	10
Elastic-RF	min_samples_leaf = 1, max_features = 3	10
RFE-RF	n_estimators = 10, max_depth = 6,min_samples_split = 3,	11

presents the optimal parameters and the number of feature variables for each algorithm.

Model accuracy was measured by the R², RMSE, and RPD of the training and test sets (

Table 4. R-squared and RMSE of cross-validation for six models using the calibration set, and R-squared, RMSE, and RPD results of independent validation using the validation set.

Model	Calibration Set		Validation Set
	R2	$\mathbf{RMSE} (\mathbf{dS}·$m−1)	R2	$\mathbf{RMSE} (\mathbf{dS}·$m−1)	RPD
GRA-PLSR	0.29	3.41	0.16	2.49	1.09
Elastic-PLSR	0.3	3.38	0.31	2.24	1.21
RFE-PLSR	0.28	3.43	0.29	2.28	1.19
GRA-RF	0.72	2.12	0.42	2.06	1.31
Elastic-RF	0.83	1.66	0.62	1.66	1.63
RFE-RF	0.62	2.49	0.72	1.47	1.84

, Figure 7). The R² ranges between 0.16 and 0.72, the RMSE ranges between 1.47 and 2.49, and the RPD ranges between 1.09 and 1.84. Notably, the PLSR algorithm exhibits lower accuracy than RF, with all three PLSR models showing RPDs below 1.4, casting doubt on their reliability. In addition, the predicted values of all three PLSR models appeared to be less tchan 0. Lower accuracy in models using the GRA algorithm for variable selection may be linked to the GCD threshold setting. Repeated experiments may help identify an optimal GCD threshold [69]. Notably, the GRA_RF model exhibits overfitting, as training set R² surpasses validation set R². The REF_RF model, with the best performance (R² = 0.71, RMSE = 1.47, RPD = 1.84), demonstrates robustness. In Figure 7, REF-RF excels in predicting non-saline soils but generally underestimates salinized soils, attributed to a limited number of high EC points. The Elastic-RF model ranks second in accuracy, excelling in predicting moderately to highly salinized soils compared to REF-RF.

Although the R² accuracy of the RFE-RF model in the training and validation sets performed moderately well, there are still some problems with the accuracy and credibility of the model in terms of RPD and RMSE values. Some of the possible factors are the problem of data distribution in the measured dataset, cultivated land film problem, image resolution problem, etc. Compared to empirical/semi-empirical models where the collection of measured data requires consideration of the spatial distribution of data sampling points, machine learning models are more sensitive to the statistical distribution of data values and, for the few ultra-high EC value points in the original dataset, the RF model smooths them by treating them directly as outliers, which results in the construction of a model where all of the medium and high values are underestimated and the low values are overestimated. Therefore, it is necessary to perform data cleaning and normalization steps before proceeding with model construction. In this study, a significant portion of the data, situated on cropland, featured the use of plastic mulch (white or black) to curb evapotranspiration and reduce irrigations, potentially altering soil and vegetation spectral reflectance characteristics, impacting vegetation index accuracy [70,71].

3.4. Importance of Variables

To validate the REF algorithm’s efficacy in reducing information redundancy, the REF-RF model was compared to a model directly constructed using the RF algorithm (Figure 8). The REF-RF model demonstrated significantly higher accuracy than the general RF model, with a 0.18 improvement in R², a 0.39 reduction in RMSE, and a 0.39 improvement in RPD. Features in the REF-RF model fell into five categories: salinity index (SI2, RE3S3, RE1S6), vegetation index (RVI), microwave backscattering coefficient (WCMVH, VV-VH), image spectral enhancement (Wetness, pc3), and topography factor (gray_diss, B2_contrast, elevation). This aligns with this article’s assumptions that, apart from the classical salinity index, vegetation index, spectral enhancement, microwave polarization index, and terrain factor are crucial modeling features. Terrain factors include elevation from DEM and texture features from remotely sensed imagery, reflected in all three feature selection methods.

Sentinel-1 variables rank prominently in the results (Figure 6c–e), with the REF algorithm identifying Com4 and WCMVH as pivotal (third and fifth, respectively). Excluding these features and re-modeling with RF resulted in reduced R² but improved RMSE. An RPD < 1.4 indicated unreliable predictions (Figure 8c), highlighting the lowest accuracy and poorest fit when microwave data was omitted. Samples with soil salinity of 0–2 dS∙m⁻¹ were mainly in cultivated land; vegetation impact was evident in June and July. Microwave data, processed with the water cloud model, effectively compensated for optical remote sensing signal saturation [44]. This underscores microwave data’s significant role in enhancing soil salinity prediction accuracy.

Satellite imagery’s spatial resolution significantly impacts soil salinity mapping [17]. In the medium-to-high salinized soil range, the vegetation type transitions from a mixture of man-made shelterbelts and natural vegetation to various salt- and drought-tolerant vegetation. At a resolution of 10 m, there may be mixed pixels, which adds some noise to models using vegetation indices. Figure 5 reveals that all variable selection algorithms excluded most vegetation indices. RFE_RF picked only one index, RVI, despite a modest 0.18 Pearson’s correlation with soil EC values. Thus, when utilizing low-resolution satellite data for arid zone salinity prediction, careful selection of vegetation indices is essential.

Soil salinity and soil surface water content interact [72] because salt is transported in the same direction as water and salt is transported with water. In addition, varying levels of soil moisture content result in reduced reflectance in the visible (VIS) and near-infrared (NIR) ranges [73], which further affects the response of various spectral indices to soil salinity [74]. The Wetness component in the RFE-RF model constructed in this article mainly reflects the humidity of soil and vegetation, and has a high importance (Figure 6e). Therefore, the water content of the surface soil is an important factor to be considered when conducting soil salinity estimation studies, and the establishment of relevant characteristic variables that respond to soil water content needs to be considered when building predictive models.

There are certain limitations in the selection of environmental covariates in this study. Salinity indices, vegetation indices, texture features, humidity features, and principal component factors are all extracted via Sentinel-2, making it challenging to avoid multicollinearity among variables. This is also a significant reason for the substantial reduction in the number of remaining variables after feature selection. According to Ivushkin et al. [31], thermal infrared imagery also performs well in distinguishing between different levels of soil salinity. Additionally, data such as meteorological data, soil texture, depth to groundwater, crop type, etc., should also be taken into account [55,57,75,76].

3.5. Characterization of the Spatial Distribution of Soil Salinity

The RFE-RF algorithm was used to map the spatial distribution of soil salinity in Minqin oasis (Figure 9a) and categorized according to different salinity classes (Figure 9c). More than 40% of the entire study area is threatened by salinization, which is similar to the findings of Ngabire et al. [77]. The soils in the study area that are seriously affected by salt are mainly located in the desert zone at the edge of the oasis, at the edge of the Badain Jaran Desert, and in the Tengger Desert. The more serious areas in the center are around Hongyashan Reservoir and Qingtuhu Lake (Figure 9d), one of which is the starting point and the other the end point of ecological water transfer from the oasis, and both of which have a high level of groundwater table [78]. Non-saline and mildly saline soils are mainly distributed in the central part of the oasis, where there are mainly residential areas and cultivated land. The degree of salinization gradually increases from southwest to northeast, which is consistent with the flow direction of the Shiyang River and reaches its maximum at Qingtu Lake, the terminal lake of the Shiyang River. Overall, salt-affected arable land reached 47.5% (Figure 9e), with salinized arable land in the central region, mainly around Qingtu Lake, and a large amount of salinized arable land occurring in the western fringe. The areas with the most severe salinization are the intermountain basins and salt ponds. The salt brought by the rising groundwater here and the accumulation of soil from the surrounding surface and runoff will form solonchak on the surface during the dry season. Coupled with the action of wind, the surface salt will accumulate [77]. Comparison of the resultant maps with the soil profiles of the study area recorded in the World Soil Database (HWSD v2.0) (Figure 9b) shows that the spatial distributions of salinized soils in the two are roughly the same. However, there are extremely high values in the HWSD data (EC = 32 dS∙m⁻¹), and some of the high-value areas do not correspond to the actual situation, which may be related to the spatial and temporal resolution of the soil profile data in the HWSD database. For example, the high-value area located in the central Shoucheng Township (rectangular area 1) and the southern high-value area located in Nanhu Township (rectangular area 2), which are mainly cultivated land, do not actually have soils with such high conductivity values, yet the low spatial resolution of the HWSD (1 KM) and the outdated time of data collection resulted in these two locations showing extremely heavy salinization. In contrast, the results obtained by the RFE_RF model were more convincing.

4. Conclusions

The problem of soil salinization in the Minqin oasis in the arid region of northwest China is evident, with more than 40% of the area salinized. Therefore, the use of advanced remote sensing techniques to map soil salinity distribution is essential for food security and rural development in the region.

This study employs Sentinel-1, Sentinel-2, SRTM images, and field-measured data from 215 locations to create 79 salinity indicators using vegetation indices, salinity indices, terrain factors, and microwave backscattering coefficient polarization combinations. Feature selection is performed through the GRA, RFE, and Elastic Net algorithms, and salt prediction models are constructed using RF and PLSR. By comparing the six models, the best model was determined to be the RFE-RF model. Key findings include the following:

• Variable selection reduces model complexity, preventing information redundancy and feature multicollinearity. Relying solely on correlation calculations may not adequately reduce dimensionality due to potential nonlinear correlations. The most effective method for feature selection is the RFE algorithm, which significantly streamlines model construction, reducing the feature variables from 34 to 11. These include six spectral features (SI2, RVI, RE3S3, RE3S6, pc3, and Wetness), two microwave backscattering coefficient features (WCMVH and Com4), two texture features (gray_diss and B2_contrast), and a topographic feature (elevation). In arid areas with low vegetation cover, featuring unutilized land (Gobi Desert) and irrigated cropland, model feature selection should consider not only vegetation and salinity indices but also other soil properties and environmental factors.
• PLSR models lack accuracy for salinity predictions; RF models, adept at capturing nonlinear relationships, are preferable for such predictions. The RFE-RF and Elastic-RF inversion models outperform basic RF models. Among them, REF-RF stands out as the superior salinity inversion model. It excels in low and medium saline soils but tends to underestimate strongly saline soils and above. Model-based feature combination offers effective soil salinity monitoring in the study area’s irrigation zone, aiding decision-making in soil management. Limiting factors include sample point distribution, cropland mulching effects on reflectance, and model parameter optimization.
• Introducing the red-edge band, the water cloud model, and the tassel-cap transformed wetness component produced satisfactory estimates. The red-edge spectral indices and the backscattering coefficients processed by the water cloud model are significantly more sensitive to soil salinity, and these effects will be more pronounced in areas with higher vegetation cover. Therefore, when choosing Sentinel data as the data source for model construction, the red-edge band of Sentinel-2 and using water cloud model pre-processing for Sentinel-1 data should be considered.
• There is a need to study the effects of other environmental variables on soil salinity levels in the future, and the GEE platform makes the whole modeling process of this study online, which facilitates teamwork and provides a convenient method for monitoring soil salinity in a large area. In the future, we will rely on the GEE platform to estimate soil salinity in a larger area and try to improve the method’s reliability further.