A novel remote sensing monitoring index of salinization based on three-dimensional feature space model and its application in the Yellow River Delta of China

Abstract

Previous studies were mostly conducted based on two-dimensional feature space to monitor salinization, while studies on dense long-term salinization monitoring based on three-dimensional feature space have not been reported. Based on Landsat TM/ETM+/OLI images and three-dimensional feature space method, this study introduced six typical salinization surface parameters, including NDVI, salinity index, MSAVI, surface albedo, iron oxide index, wetness index to construct eight different feature space monitoring index. The optimal soil salinization monitoring index model was proposed base on field observed data and then the evolution process of salinization in Yellow River Delta (YRD) were analyzed and revealed during 1984–2022. The salinization monitoring index model of MSAVI-Albedo-I_Fe2O3 feature space had the highest accuracy with R² = 0.93 and RMSE = 0.678g/kg. The spatial distribution of salinization in YRD showed an increasing trend from inland southwest to coastal northeast and the salinization intensity showed an increasing trend during 1984–2022 due to the implements of agricultural measures such as planting salt-tolerant crops, microbial remediation and fertility improvement. The rate of salinization deterioration in the northeast part was greater than others. Zones of salinization improvement were mainly located in cultivated land of the southwest parts.

Keywords:

• Soil salinization
• 3-D feature space
• Landsat images
• spatial distribution
• vegetation index

1. Introduction

Soil salinization, one of the typical land degradation types, can be attributed to the comprehensive actions of natural factors and human activities (Chang et al. 2015 Luo 2015,; Wang et al. 2018; Abdennour et al. 2021). It seriously affects the stability and sustainable development of the ecological environment. Saline soil can be divided into inland saline soil and coastal saline soil according to the geo-graphical location. Coastal saline soil is the main type of saline soil in Shandong, which is mostly distributed in YRD and the southern coast of Laizhou Bay (Liu 2016). The modern YRD with only 160 years has the youngest land in China since the Yellow River diverted into the Bohai Sea in 1855, which was bred by the sediment accumulation of the Yellow River (Liu 2016; Wang et al. 2018). The problem of salinization in YRD is particularly prominent, which greatly restricts the local economic development (Li et al. 2020; Wang et al. 2020; Duan et al. 2021), seriously threatening the stability of wetland ecosystem and the layout of regional ecological protection red line (Yao et al. 2006). Therefore, accurate monitoring of soil salinization, and then grasping its spatial and temporal evolution patterns can provide approach and decision supports for the accurate remediation and treatment for the salinization in YRD.

With the launch of more satellites, remote sensing has become one of the important approaches to realize large-scale salinization monitoring (Wang et al. 2018; Guo et al. 2022a). Cao et al. (2014) explored the feasibility of extracting salinization information by remote sensing in YRD based on Landsat satellites. In recent years, many scholars have utilized typical parameters that inversed from remote sensing images to construct two-dimensional feature space for monitoring salinization information and achieved good results with the overall accuracy of 80–90% (Jia et al. 2020; He et al. 2021). Ha et al. (2009) selected the salinity index and surface reflectance to construct feature space model, and then obtained the spatial distribution of soil salinization in Keriya oasis. Wang et al. (2010) utilized the NDVI and salt index to propose the NDVI-SI feature space model, and then quantitatively analyzed the salinization of Yutian Oasis in Xinjiang. Ding et al. (2013) proposed the MSAVI-WI feature space model to monitor soil salinization in arid areas. Ding et al. (2013) analyzed and explored the effects of three different land use types on the physical and chemical properties of saline-alkali soils in the YRD using a combination of field investigation and indoor analysis. Wang et al. (2016) analyzed the spatial variability of soil water and salinity in Kenli District of the YRD and the micro-domain characteristics of soil salinity in the study area were explored in summer. Zou et al. (2017) investigated the temporal and spatial evolution characteristics of soil salinity in the YRD under reclamation and natural conditions, and found that artificial reclamation could reduce the degree of soil salinization to a certain extent. Guo et al. (2019a) developed a novel feature space monitoring index model of salinization based on vegetation index and salt index in YRD. However, soil degradation caused by salinization is often affected by both man-made and natural factors, known as primary or secondary salinization (Abdennour et al. 2019). The above studies are mostly based on the two-dimensional feature space model to monitor soil salinization information, and cannot consider the influence of the interaction of multiple types of factors on the inversion of salinization information. In mathematics, the number of variables in three-dimensional space is more than that in two-dimensional space, and three-dimensional space can consider the influence of more interaction factors on salinization. In addition, the salinization process is always influences by many natural-artificial factors. Therefore, three-dimensional space is more conductive to improve the inversion accuracy of salinization than two-dimensional space. In addition, previous studies are mostly conducted to analyze the evolution process and pattern of salinization from the perspective of sparse time series. Li and Gao (2018) utilized the spatial overlay analysis and geostatistical simulation methods to analyze the spatial and temporal evolution of soil salinization from 2013 to 2016 based on ArcGIS 10.2. Li et al. (2020) adopted the long-term sequence Landsat images from 1976 to 2018 to investigate the relationship between the evolution of coastal erosion and deposition in the YRD since 1976. Although it can reveal the temporal and spatial variation characteristics of regional salinization to a certain extent, it is difficult to monitor and capture the mutation process of salinization and its critical threshold information.

In this study, a three-dimensional feature space salinization monitoring model has been constructed and proposed by introducing multiple types of surface parameters, and then the long time series data sets of salinization from 1984 to 2022 have been retrieved. Finally, the evolution mode and characteristics of salinization in the modern YRD have been analyzed from the perspectives of gradual change and mutation.

2. Materials and methods

2.1. Study area

The Yellow River Delta(YRD) (Figure 1), located in northern Shandong Province of China (118cha′–119located in northern Shandong Province of China (118c0 km². The YRD is a warm temperate semi-humid continental monsoon climate with clear four seasons and flat terrain (Li et al. 2011; Guo et al. 2019a). The temperature is higher in summer, while that of winter climate is cold, with a large temperature difference. The average annual precipitation in YRD is 594 mm, and its spatial distribution is very uneven (Liu et al. 2015; Chen et al. 2017; Guo et al. 2019b). The maximum evaporation is in summer, followed by that in spring and autumn. The study region is prone to drought with great evaporations and relative less precipitation. The terrain is high in southwest and low in northeast. Due to its unique geographical location, the YRD has high groundwater level and high groundwater mineralization (Guo et al. 2022b; Guo et al. 2023). At the same time, frequent storm surges and seawater recharge further exacerbated the salinization condition. The YRD is located in the transition zone of atmosphere, ocean, river and land. It has typical multiple ecological interfaces and fragmented surface landscape patterns (especially wetlands), which leads to severe environmental gradients and rapid salinization succession.

Figure 1. Overview of the study area.

2.2. Data collection and preprocessing

The datasets of Landsat5 TM, Landsat7 ETM + and Landsat8 OLI images (Path/Row, 121/34, Leavel-1) during 1984–2022 were obtained from the geographic spatial data cloud (http://www.gscloud.cn/). The spatial resolution of the images was 30 m, and the data acquisition time was mainly concentrated in March, April and May (Table 1). The cloud coverage of most images was less than 5%. Geometric correction module of ENVI5.3 was applied to calibrate the geometric deformation and error of the images and Fast Line-of-sight Atmospheric Analysis of Spectral Hypercube (FLAASH) module of ENVI5.3 was utilized to eliminate the influence of atmospheric scattering and absorption (Hu et al. 2014;Liu et al. 2022). In addition, the Landsat Gapfill tool was applied to repair the bands of Landsat7 image to eliminate the periodic scanning bands (Tang et al. 2005).

Table 1. List of Landsat images and its detail information.

Id	Time	Sensor	Resolution/m	Id	Acquisition date	Data type	Resolution/m
1	1984.05.30	TM	30	21	2004.04.19	TM	30
2	1985.03.14	TM	30	22	2005.04.22	TM	30
3	1986.04.18	TM	30	23	2006.04.09	TM	30
4	1987.04.21	TM	30	24	2007.03.27	TM	30
5	1988.04.07	TM	30	25	2008.04.14	TM	30
6	1989.04.10	TM	30	26	2009.04.01	TM	30
7	1990.03.12	TM	30	27	2010.05.06	TM	30
8	1991.04.16	TM	30	28	2011.03.22	TM	30
9	1992.04.02	TM	30	29	2012.04.01	ETM+	30
10	1993.05.23	TM	30	30	2013.05.30	OLI	30
11	1994.04.24	TM	30	31	2014.05.01	OLI	30
12	1995.03.26	TM	30	32	2015.05.04	OLI	30
13	1996.05.31	TM	30	33	2016.03.03	OLI	30
14	1997.02.11	TM	30	34	2017.04.23	OLI	30
15	1998.04.19	TM	30	35	2018.03.25	OLI	30
16	1999.04.06	TM	30	36	2019.04.13	OLI	30
17	2000.04.08	TM	30	37	2020.04.15	OLI	30
18	2001.03.10	TM	30	38	2021.04.18	OLI	30
19	2002.04.14	TM	30	39	2022.04.05	OLI	30
20	2003.04.01	TM	30

2.3. Methods

2.3.1. Principle of feature space

Feature space refers to the plane space constructed by typical parameters as horizontal and vertical axes in the coordinate system (Ding et al. 2013; Guo et al. 2019b). This model has been widely used in monitoring the process of land degradation, including desertification, rocky desertification and salinization (Ha et al. 2009; Guo et al. 2020; Guo et al. 2023). Based on the feature space, this study selected the optimal parameter combination to establish the feature space index model for salinization monitoring. Studies of salinization on field observed data can obtain accurate results, which are not suitable for large-scale salinization monitoring (Li et al. 2017; Wu et al. 2018). With the development of earth observation technology in recent years, remote sensing monitoring technology has been widely used in salinization monitoring (Li et al. 2012). Many studies are conducted based on fuzzy analytic hierarchy process, comprehensive index method and image classification method, which would amplify the effects of some factors on salinization (Wu et al. 2018; Wang et al. 2020). Moreover, these methods cannot explain the interaction between surface parameters and ignore the spatial characteristics of each parameter (Gu et al. 2021). In recent years, many scholars have introduced some surface parameters to construct feature space models for salinization monitoring, and have achieved good results with the overall accuracy of 80–90% (Ding et al. 2013; Wang et al. 2013; Yao et al. 2021). The feature space model can better consider the interactions of typical surface parameters on the process of salinization. Therefore, the feature space model has been applied to monitor the salinization in this paper. With the aggravation of salinization, the surface soil salt accumulation is not conducive to the growth of vegetation, and the vegetation coverage shows a decreasing trend (Figure 2) (Meng et al. 2016;Yao et al. 2021). With the accumulation of saline in surface soil, Albedo and I_Fe2O3 both show an increasing trend. In this study, MSAVI-Albedo-I_Fe2O3 feature space has been taken as an example, A represents the severe salinization zone (Low MSAVI, high Albedo, and high I_Fe2O3), C represents the non-salinization zone (High MSAVI, low Albedo, and low I_Fe2O3). The greater the distance from any point M to point C (M, N, P), the more severe soil salinization. On the contrary, the soil salinization degree is slighter.

Figure 2. Principle of Feature space.

2.3.2. Calculation of typical parameters

Considering the special natural ecological landscape and soil-forming background of saline soil in YRD, six typical surface parameters of salinization were selected in this study: normalized differences vegetation index (NDVI), modified soil adjusted vegetation index(MSAVI), salinity index (SI), surface albedo (Albedo), iron oxide index (I_Fe2O3), wetness index (WI). NDVI and MSAVI can accurately reflect the grown condition and changes of vegetation (Sun et al. 2018; Ma et al. 2018). MSAVI can eliminate the influence of bare soil in vegetation index. SI is an important indicator to characterize and indicate the salinization information, and salt in seawater is one of the important sources for salinization in YRD (Lu et al. 2020; Guo et al. 2022c). Surface albedo refers to the ratio of solar shortwave radiation to the total solar radiation in all directions, which is an important parameter for numerical climate model and surface energy balance (Betts et al. 1997; Xiao et al. 2011). With the exacerbation of soil salinization, the vegetation coverage of non-salt-tolerant plants showed a decreasing trend. The salt accumulation process in the surface soil would increase the surface albedo (Stroeve et al. 2005; Alejandro et al. 2007; Wang et al. 2010; Wang et al. 2017). There are significant differences in humidity between coastal and inland areas due to the special location of YRD. When drought occurs, the soil dryness increases (decreasing wetness index), and the underground mineralized water molecules move upwards, and a large amount of salt accumulates in the surface soil, resulting in occurrence of salinization. Therefore, the WI can better reflect the soil salinization information (Zhang et al. 2008; Hu et al. 2019). Soil salinization can not only influence the growth and types of vegetation in salinization areas, but also affects the composition and changes of chemical substances in the soil. Thus, I_Fe2O3 is one of the important indicators of soil salinization in YRD (Guo et al. 2020; Peng et al. 2013). The calculation formulas for the above six typical parameters are shown in Table 2:

Table 2. Formula of surface parameters.

Surface parameters	Formula
NDVI (Sun et al. 2018)	$\text{NDVI}=(\text{Nir}-\text{Red})/(\text{Nir}+\text{Red})$
MSAVI (Ma et al. 2018)	$\text{MSAVI}=(2\text{Nir}+1-\sqrt{{(2\text{Nir}+1)}^2-8(\text{Nir}-\text{Red})})/2$
SI (Lu et al. 2020)	$\text{SI}=\sqrt{\text{Blue}\times \text{Red}}$
Albedo (Guo et al. 2022c)	$\text{Albedo}=(0.356\times \text{Blue}+0.130\times \text{Red}+0.373\times \text{Nir}+0.085\times {\text{Swir}}_1+0.072\times {\text{Swir}}_2-0.0018)/1.016$
WI (Hu et al. 2019)	$\text{WI}=0.1446\times \text{Blue}+0.1761\times \text{Green}+0.3322\times \text{Red}+0.3396\times \text{Nir}\mathrm{ }-0.6210\times {\text{Swir}}_1-0.4186\times {\text{Swir}}_2$
IFe2O3 (Peng et al. 2013)	${\mathrm{I}}_{{\text{Fe}}_{{}^2}{\mathrm{O}}_{{}^3}}=\frac{\text{Red}}{\text{Nir}}$

Note: Blue, Green, Red, Nir, Swir₁ and Swir₂ refer to the reflectance of blue band, green band, red band, near infrared band, and short infrared bands, respectively.

2.3.3. Index standardization

Since the dimensions of surface parameters are different, the indices should be standardized to enhance the inversion accuracy of salinization in this study (Eq. (1)(1) $$P_i=(S_i-S_{i,\min })/(S_{i,\max }-S_{i,\min })$$(1) ): (1) $$P_i=(S_i-S_{i,\min })/(S_{i,\max }-S_{i,\min })$$(1) where $P_i$ refers to the standardized indicator of i; $S_i$ is the original index of i; $S_{i,\min }$ is the minimum value of the original index of i; $S_{i,\max }$ is the maximum value of the original index of i.

2.3.4. Trend analysis method

The average salinization monitoring index during certain periods is obtained by the Raster Calculator of ArcGIS 10.7 (Hong et al. 2016; Eq. (2)(2) $$\overline{\text{SMI}}=\frac{\sum _{i=1}^n{\text{SMI}}_{\mathrm{i}}}{n}$$(2) ). (2) $$\overline{\text{SMI}}=\frac{\sum _{i=1}^n{\text{SMI}}_{\mathrm{i}}}{n}$$(2) where $\overline{\text{SMI}}$ is the average value of a pixel in n years, SMI_i is the salinization monitoring index value for the first year of a pixel.

The linear regression analysis method can eliminate the influence of salinization mutation in specific years to some extent. The equation is as follows (Eq. (3)(3) $$C=\frac{\sum _{j=1}^n{\text{SMI}}_j{t}_j-\frac{1}{n}\sum _{j=1}^n{\text{SMI}}_j\sum _{j=1}^n{t}_j}{\sum _{j=1}^n{t}_j^2-\frac{1}{n}{\left(\sum _{j=1}^n{t}_j\right)}^2}$$(3) ). (3) $$C=\frac{\sum _{j=1}^n{\text{SMI}}_j{t}_j-\frac{1}{n}\sum _{j=1}^n{\text{SMI}}_j\sum _{j=1}^n{t}_j}{\sum _{j=1}^n{t}_j^2-\frac{1}{n}{\left(\sum _{j=1}^n{t}_j\right)}^2}$$(3) where C represents the gradient of the trend line, n is the total number of the years (n = 39). If C > 0, it indicates that SMI shows an increasing trend, vice versa.

2.3.5. Gravity center model

The gravity center can reflect the non-uniformity and bias of spatial and temporal distribution of geographical elements, which has been widely used in the fields of land use, economy, population and ecology (Liang et al. 2021; Wang et al. 2022; Guo et al. 2023). In this paper, the gravity center can indicate the bias of the salinization. In addition, the migration trajectory of gravity center can reflect the change trend of salinization in spatial distribution (Yang et al. 2014). The gravity center of a region is defined as (x, y), and ${\text{MAI}}_i$ is set as the attribute value of the ith plane space unit (grid). Its Cartesian coordinates are set as $(x_i,y_i)$ (Eqs. (4)-(5)). The gravity center of n plane space units (grid) is defined as a point $(\overline{x},\overline{\mathrm{y}}).$ (4) $$\overline{x}=\frac{\sum _{i=1}^n{\text{MAI}}_i{x}_i}{\sum _{i=1}^n{\text{MAI}}_i}$$(4) (5) $$\overline{y}=\frac{\sum _{i=1}^n{\text{MAI}}_i{y}_i}{\sum _{i=1}^n{\text{MAI}}_i}$$(5)

3. Results

3.1. Construction of feature space

Based on the standardized surface parameters (NDVI, MSAVI, SI, Albedo, I_Fe2O3, WI) that derived from Landsat images in April 2020, the three-dimensional feature spaces that composed of different surface parameters were constructed by Matlab2022. In order to eliminate the interference of urban impervious surface and artificial building area on salinization monitoring, these zones were eliminated combined with the national land use data. Eight 3-D feature spaces were shown in Figure 3.

Figure 3. Three-dimensional feature spaces that composed of different surface parameters (a) MSAVI-Albedo-IFe2O3; (b) NDVI-WI- IFe2O3; (c) NDVI-Albedo- IFe2O3; (d) MSAVI-Albedo-SI; (e) NDVI-Albedo-SI; (f) NDVI-WI-SI; (g) MSAVI-WI-IFe2O3; (h) MSAVI-WI-SI.

3.2. Establishment of salinization monitoring index model (MSAVI-Albedo-I_Fe2O3 feature space as an example)

Taking MSAVI-Albedo-I_Fe2O3 feature space as an example, the distribution laws of different degrees of salinization soil in feature space were significantly different. Based on the distance from any point to point (M, N, P), this study set four point groups that distributed in different regions of the feature space. Through the 19 field observed samples in each point sets, the corresponding relations between different soil salinization degrees and four point sets were determined. As shown in Figure 4, red point sets were mainly distributed in non-salinization zones, while orange point sets were mainly concentrated in slight salinization zones. The yellow and green point sets were mainly located in mild and severe salinization areas, respectively. Therefore, different point sets for levels of soil salinization were better identified and distinguished in MSAVI-Albedo-I_Fe2O3 feature space, which was consistent with the results of field observed data.

Figure 4. Spatial distribution of levels of soil salinization in three-dimensional feature space. (a) Non salinization; (b) Slight salinization; (c) Mild salinization; (d) Severe salinization.

As shown in Figure 5, MSAVI, Albedo and I_Fe2O3 had obvious nonlinear relationships in 3-D feature space. With the decrease of MSAVI, Albedo and I_Fe2O3 both increased, and the salinization showed an exacerbating trend. The larger the distance from any point M to point O (M, N, P) in the feature space, the more serious the corresponding salinization. The distance L1 could be utilized to distinguish different levels of salinization (Guo et al. 2019b; Guo et al. 2022c), so MSAVI-Albedo-I_Fe2O3 feature space monitoring index (MAI) was constructed as follows (Eq. (6)(6) $$\text{MAI}={\mathrm{L}}_1=\sqrt{{\mathrm{(}\text{MSAVI}-\mathrm{M}\mathrm{)}}^2{+\mathrm{(}\text{Albedo}-\mathrm{N}\mathrm{)}}^2{+\mathrm{(}\mathrm{I}}_{{\text{Fe}}_{{}^2}{\mathrm{O}}_{{}^3}}-\mathrm{P}{)}^2}$$(6) ): (6) $$\text{MAI}={\mathrm{L}}_1=\sqrt{{\mathrm{(}\text{MSAVI}-\mathrm{M}\mathrm{)}}^2{+\mathrm{(}\text{Albedo}-\mathrm{N}\mathrm{)}}^2{+\mathrm{(}\mathrm{I}}_{{\text{Fe}}_{{}^2}{\mathrm{O}}_{{}^3}}-\mathrm{P}{)}^2}$$(6)

Figure 5. Establishment of 3-D feature space monitoring index model of salinization.

3.3 Accuracy verification and comparisons

Figure 6 showed that different salinization monitoring index were obtained based on the raster calculator of ArcGIS 10.7 and the established monitoring models were shown in Table 3.

Figure 6. Spatial distribution of different salinization models (a) MSAVI-Albedo-IFe2O3; (b) NDVI-WI- IFe2O3; (c) NDVI-Albedo- IFe2O3; (d) MSAVI-Albedo-SI; (e) NDVI-Albedo-SI; (f) NDVI-WI-SI; (g) MSAVI-WI- IFe2O3; (h) MSAVI-WI-SI.

Table 3. Accuracy comparison of different salinization monitoring index models.

Feature space	Formula	R^2	RMSE (g/kg)
MSAVI-Albedo-IFe2O3	$\text{MAI}=\sqrt{{\mathrm{(}\text{MSAVI}-1\mathrm{)}}^2{+\mathrm{(}\text{Albedo}-1\mathrm{)}}^2{+\mathrm{(}\mathrm{I}}_{{\text{Fe}}_{{}^2}{\mathrm{O}}_{{}^3}}{)}^2}$	0.93	0.678
NDVI-WI- IFe2O3	$\text{NWI}=\sqrt{{\mathrm{(}\text{NDVI}-1\mathrm{)}}^2{+\text{WI}}^2{+\mathrm{(}\mathrm{I}}_{{\text{Fe}}_{{}^2}{\mathrm{O}}_{{}^3}}{)}^2}$	0.83	0.799
NDVI-Albedo- IFe2O3	$\text{NAI}=\sqrt{{\mathrm{(}\text{NDVI}-1\mathrm{)}}^2{+\mathrm{(}\text{Albedo}-1\mathrm{)}}^2{+\mathrm{(}\mathrm{I}}_{{\text{Fe}}_{{}^2}{\mathrm{O}}_{{}^3}}{)}^2}$	0.81	0.987
MSAVI-Albedo-SI	$\text{MAS}=\sqrt{{\mathrm{(}\text{MSAVI}-1\mathrm{)}}^2{+\mathrm{(}\text{Albedo}-1\mathrm{)}}^2{+\text{SI}}^2}$	0.79	1.124
NDVI-Albedo-SI	$\text{NAS}=\sqrt{{\mathrm{(}\text{NDVI}-1\mathrm{)}}^2{+\mathrm{(}\text{Albedo}-1\mathrm{)}}^2{+\text{SI}}^2}$	0.76	1.356
NDVI-WI-SI	$\text{NWS}=\sqrt{{\mathrm{(}\text{NDVI}-1\mathrm{)}}^2{+\text{WI}}^2{+\text{SI}}^2}$	0.74	1.647
MSAVI-WI- IFe2O3	$\text{MWI}=\sqrt{{\mathrm{(}\text{MSAVI}-1\mathrm{)}}^2{+\text{WI}}^2{+\mathrm{(}\mathrm{I}}_{{\text{Fe}}_{{}^2}{\mathrm{O}}_{{}^3}}{)}^2}$	0.67	1.966
MSAVI-WI-SI	$\text{MWS}=\sqrt{{\mathrm{(}\text{MSAVI}-1\mathrm{)}}^2{+\text{WI}}^2{+\text{SI}}^2}$	0.61	2.354

In order to verify the inversion accuracy of the proposed models, in this study, the univariate linear regression analysis (correlation coefficient and root mean square error) was performed based on the 64 field observations (soil salinity of 0–20 cm) and the monitoring index (p < 0.01). The results showed that (Table 3) there were significant differences in the inversion accuracy of different three-dimensional feature space monitoring index models. The salinization monitoring index model based on MSAVI-Albedo-I_Fe2O3 feature space had the best applicability in YRD with R² = 0.93 and RMSE space hg/kg, which was higher than previous studies with the R² ranging from 0.8–0.9 (Wang et al. 2010; Ding et al. 2013; Guo et al. 2019a; Guo et al. 2020; Yao et al. 2021). Taking the MSAVI-Albedo-I_Fe2O3 feature space model as an example, the accuracy verifications based on univariate linear regression analysis model was shown in Figure 7.

Figure 7. Accuracy verification of the feature space model.

3.4. Temporal variation of salinization during 1984–2022

The salinization remote sensing monitoring indexes (SMI) during 1984–2022 were obtained by the optimal model of MSAVI-Albedo-I_Fe2O3 feature space and ArcGIS 10.7. The average SMI in the YRD during 1984–2022 showed obvious temporal change laws (Figure 8). The average remote sensing monitoring index of salinization in recent 39 years was 0.97. During 1984–2022, the SMI showed an overall increasing trend. The maximum SMI in 2016 was 1.06, while that of 1996 was the smallest of 0.85. The average SMI in 2006–2012 was larger than the average value of 1984–2022. The average SMI values in 1984–1986, 1988–2000, 2013–2015 and 2019–2022 were lower than that in recent 39 years.

Figure 8. Interannual variation of average salinization monitoring index.

3.5. Spatial distribution pattern of salinization during 1984–2022

In order to explore the spatial differentiation patterns of salinization in YRD, the salinization monitoring index(SMI) was divided into four grades by using the Natural Breaks method of ArcGis10.7 combined with the field observed data, 《Outline of Planning for Improvement and Utilization of Saline-alkali Land》, and the vegetation coverage: SMI < 0.8 as non salinization, 0.8 < SMI < 1.0 as slight salinization; 1.0 < SMI < 11.2 as mild salinization; SMI > 1.2 as severe salinization (Wang et al. 2020).

There were significant differences in spatial distribution for levels of salinization during 1984–2022 (Figure 9). The slight salinization zone had the largest area of 3448.75 km², accounting for 49.02%, mainly concentrated in the central inland areas. The area of mild salinization was 1504.63 km², accounting for 21.39%. The area of non-salinization was 1368.32 km², accounting for 19.45%, which was mainly distributed in the central and southwestern inland areas. Zone of severe salinization had the smallest area of 709.57 km², accounting for 10.14%. Mild and severe salinization zones were mostly located in the northeast and coastal areas. Because in these regions, the terrain was lower, prone to seawater recharge and other phenomena. In addition, the depth of groundwater was shallow.

Figure 9. Spatial distribution of salinization in Yellow River Delta.

3.6. Change intensity of salinization during 1984–2022

In order to further analyze the spatial and temporal changes of salinization in YRD in recent 39 years, the change trend coefficient was calculated by ArcGIS 10.7 and Eq.3(3) $$C=\frac{\sum _{j=1}^n{\text{SMI}}_j{t}_j-\frac{1}{n}\sum _{j=1}^n{\text{SMI}}_j\sum _{j=1}^n{t}_j}{\sum _{j=1}^n{t}_j^2-\frac{1}{n}{\left(\sum _{j=1}^n{t}_j\right)}^2}$$(3) (Figure 10). The trend coefficient (−0.04 to 0.04) were divided into five different levels, including severely reduced (−0.04 to 0.02), slightly reduced (−0.02 to 0.005), stable (−0.005 to 0.005), slightly increased (0.005–0.02), and severely increased (0.02–0.04).

Figure 10. Change trend and intensity of salinization during 1984-2022.

The area of salinization improvement (severely reduced and slightly reduced) was 698.76 km², accounting for 10.16%, which was mainly located in the southwestern YRD, mainly distributed in Lijin County, the west and middle of Dongying District, and the west of Guangrao County. The area of stable zones was 4407.30 km², accounting for 64.09%, mainly distributed in the middle of Kenli District and the south of Hekou District. The area of salinization deterioration(slightly increased and severely increased) was 1770.26 km², accounting for 25.75%, mainly concentrated in the northern and eastern coastal parts of the YRD, mainly distributed in the north of Hekou District, the east of Kenli District, the east of Dongying District and the east of Guangrao County.

3.7. Migration trajectory of salinization gravity center during 1984–2022

In order to better reveal the uneven and bias of temporal and spatial changes of salinization, the gravity center was utilized to analyze the migration trajectory of soil salinization gravity center during 1984–2022.The gravity center is generally expressed in longitude and latitude of the map (Lv et al. 2020;Wei et al. 2022). In this study, the gravity centers of salinization in the YRD during 1984–2022 were calculated with the tool of Spatial Analyst of ArcGIS 10.7 (Figure 11). The results showed that the gravity centers of salinization were mainly concentrated in the junction of Kenli District and Lijin County, indicating that the gravity center of soil salinization shifts from southwest to northeast, and the condition of soil salinization in the northeast was more significant than that of the southwest parts. The migration trajectories of salinization gravity center at different time scales (5-year, 10-year, 20-year and 40-year) were consistent with slighter differences. During 1984–2022, the gravity center of salinization in the YRD moved northeastward on the whole. At 5-year and 10-year scales, the gravity center of salinization shifted slightly in the southwest of the area where the gravity center occurred (before 2013), and then quickly moved northward (after 2013), but the gravity center of salinization at different time scales was located in the northern junction of Kenli District and Lijin County. This indicated that the degree of salinization in the northeast of the YRD was more intensive than that in the southwest from 1984 to 2013, while the salinization in the north was significantly aggravated from 2013 to 2022, and the deviation of salinization center was larger than that before 2013. Because after 2013, a large area of reservoir ponds was built in Hekou District and Kenli District for aquaculture. Although the degree of salinization increased, the rational development and utilization of salinized land, and the establishment of farmland industry demonstration areas in the southwest region contributed greatly to the improvements of the soil salinization. Therefore, the gravity center of salinization moved northeast after 2013, which indicated that the quality of eco-environment of the northeast parts was greatly enhanced due to the rational development and utilization of saline soil, such as screening salt-tolerant plant, irrigation and drainage engineering measures improvement, and chemical and biological measures improvement (Yu et al. 2022).

Figure 11. Migration trajectory of salinization gravity center during 1984-2022.

4. Discussion

4.1. Advantages of the proposed three-dimensional feature space model

Current researches were mostly conducted based on two-dimensional feature space, while the influence of the interaction between multiple types of factors on the salinization process was ignored. The salinization process in the YRD was comprehensively affected by the interaction of natural and artificial factors. Based on the construction of three-dimensional feature space, this study could take into account the influence of interactive actions among different types of factors on the salinization process. Previous studies utilized the two-dimensional feature space models to monitor the salinization process, which obtained better results. Ding et al. (2013) proposed a concept of MSAVI-WI feature space and established the soil salinity monitoring index (MWI) model and the results indicated that there was significant correlation between the MWI and surface soil salinity with R² = 0.844. Guo et al. (2019a) proposed a rapid monitoring of soil salinization in the YRD based on VI-SI feature space and found that the MSAVI-SI₁ feature space model index had the highest inversion accuracy with R² = 0.912. Guo et al. (2020) utilized five typical surface parameters to construct ten different feature spaces and found that the inversion accuracy of the I_Fe2O3 feature space detection index based on the point to point model was the highest with R² = 0.86. The formation and evolution of salinization were always comprehensively influenced by natural factors and human activity, such as precipitation, vegetation, evaporations, topography, and underground water (Guo et al. 2019b). The two-dimensional feature space model could only consider two surface parameters, which would reflect the influences of some factors (Guo et al. 2019a). The three-dimensional feature space model would take more influencing factors into considerations to comprehensively reflect the salinization evolutions (Guo et al. 2020). In this paper, The salinization monitoring index model of MSAVI-Albedo-I_Fe2O3 feature space had the highest accuracy with R² = 0.93 and RMSE = 0.678g/kg. The reason was that the MSAVI could better reflect the vegetation condition in zones with sparse vegetation region. In addition, with the aggravation of salinization, more salt would accumulate in the surface layer of the soil, the surface was more sensitive to this change process (Wang et al. 2010). Soil solute transport process was the direct cause of salinization, so that the Fe⁺³ was a better indicator to reflect and monitor the salinization evolution (Ha et al. 2009; He et al. 2021; Wu et al. 2021).

4.2. Causes of spatial and temporal distribution of salinization in the YRD

The salinization of the YRD showed a decreasing trend from northeast to southwest and from coastal to inland in spatial distribution (Zhang et al. 2008; Zhang et al. 2011). The annual gravity centers of salinization were mainly concentrated in the junction of Kenli District and Lijin County. Moreover, the gravity center of salinization moved to northeast of the study area during 1982–2022. The salinization intensity in the coastal area of the northeast of the YRD was higher than that in the inland area of the northwest (Wang et al. 2020; Li et al. 2011). Affected by land-sea interaction, the coast area of the YRD had dramatically dynamic changes. The salt of saline soil in YRD that developed from saline silt in the coastal regions mainly came from sea water (Duan et al. 2021). The reclamation of coastal areas before 2000 was mostly unreasonable, which led to severe salinization. In addition, irrational irrigation and extensive management led to extremely serious secondary salinization in the region (Li et al. 2011; Chen et al. 2021). In addition, the tidal flat was one of the typical land use types in the YRD. The tidal flat referred to the zonal that distributed parallel to the coastline. Due to the influence of seawater, the groundwater level was high and the soil contains more salt (An et al. 2017). Salinization improvement zones were mainly located in cultivated land of the southwest parts, while stable zones were distributed in the central inland area. During the past decades, on the basis of perfecting irrigation and drainage water conservancy facilities, comprehensive improvement technical measures are integrated with agricultural measures such as rational farming, planting salt-tolerant crops, microbial remediation and fertility improvement (Zhang et al. 2015; Cheng et al. 2022).

5. Conclusion

Fully considering the causes and environmental background of formation and evolution of salinization in the YRD, six typical surface indexes were selected and extracted from Landsat images, and eight three-dimensional feature space monitoring index models of salinization were construed. Based on filed observed data, the optimal monitoring index model was determined, and then the temporal and spatial variation characteristics of salinization in the YRD were explored. The main conclusions were as follows: (1) The salinization monitoring index model of MSAVI-Albedo-I_Fe2O3 feature space had the highest accuracy with R² = 0.93 and RMSE = 0.678g/kg, which indicated that the new proposed model had better applicability in the salinization monitoring for the YRD;(2) The salinization in the northeast parts in the YRD was more serious than that of southwest parts. In addition, there was an improvement trend of salinization during 1984–2022;(3) The aggravation degree of salinization in the northeastern parts was higher than that in other parts; (4) Salinization improvement zones were mainly located in cultivated land of the southwest parts, while the severe salinization zones mainly occurred in the wetland in the northeast regions.

In recent years, more surface parameters can be utilized to indicate and analyze the process of salinization, the optimization of feature variables subset has become an urgent problem to be solved in this research field. It is necessary to use the random forest algorithm to select the optimal surface parameters to construct the feature space model in future studies. Meanwhile, more surface parameters should be utilized to construct the n-dimensional (n > 3) feature space model index to improve the inversion accuracy of salinization, which could consider the interactions between multi-type factors.

Acknowledgements

This work was supported by Natural Science Foundation of Shandong Province (grant no. ZR2021MD047); National Natural Science Foundation of China (grant no.42101306); Open fund of Key Laboratory of National Geographic Census and Monitoring, MNR (grant no.2020NGCM02); the Open Fund of Key Laboratory of Urban Land Resources Monitoring and Simulation, Ministry of Natural Resources (grant no. KF-2020-05-001), and Agricultural Science and Technology Innovation Program (grant no. CAAS-ZDRW202201).

Data availability statement

The data that support the findings of this study are available from the corresponding author, Guo B, upon reasonable request.

Disclosure statement

No potential conflict of interest was reported by the authors.