Modeling and Locating the Wind Erosion at the Dry Bottom of the Aral Sea Based on an InSAR Temporal Decorrelation Decomposition Model

1. Introduction

The fine particles produced by wind erosion are essential for sand and dust activities. The sand and dust activities, especially the salt dust activities, are a disaster for surrounding residents and ecosystems in the Aral Sea [1]. Wind erosion accelerated the desertification process and promoted the formation of the Aralkum Desert. Sufficient fine particles from this new desert provide favorable conditions for sand and dust activities. The salt storms generated by winds caused various diseases to the surrounding residents and result in the death of a large area of vegetation, especially the crops. Accurate localization of these wind erosion areas is necessary for human intervention in salt and dust storms. Scientists developed numerous models based on traditional experimental physics and remote sensing technology to distinguish different degrees of wind erosion.

Wind erosion research results based on traditional experimental physics provide a solid theoretical foundation for wind erosion research. According to the study of Liu and Zobeck based on experimental physics methods, wind blowing, impact and abrasion of moving particles are the three primary forms of particle moving and separation [2,3]. Bagnold derived the threshold of the starting wind speed of the particle fluid according to the moment balance equation [4]. The research results of Lu and Shao show that when the ground particles are mainly dust, the vertical dust flux caused by the impact is proportional to the 3-4 power of the wind speed [5]. Numerous scientists have studied the forms and laws of the movement of sand and dust and provided many physical models of different forms of movement [3,6–15]. Particles suspended in the air fall to the ground mainly through dry and wet sedimentation. The physical processes involved in dry and wet sedimentation are relatively complex, and there is no effective means to describe the dry and wet sedimentation process [6,16–18]. Clarifying the influencing factors of the wind erosion process is the basis for establishing the wind erosion model, and Chepil’s classification of wind erosion factors plays an essential role in advancing the research on the wind erosion process [19]. Based on Chepil’s research results, subsequent researchers have proposed many well-known wind erosion models [20–25]. The research of wind erosion process based on traditional experimental physics has essential value and significance for human understanding of the fundamental laws of wind erosion activity. However, the uncertainty of the wind erosion model itself and the difficulty of obtaining high spatial and temporal resolution data greatly limit the wind erosion model in the application of describing large-scale wind erosion scenarios.

InSAR (Interferometric Synthetic Aperture Radar) decorrelation can be used to describe the random changes of the surface due to wind erosion, expecting to solve the problem of quantitative characterization of the wind erosion intensity [26]. Coherence is commonly employed to assess the similarity of InSAR echo signals. It quantifies the correlation between two complex InSAR echo signals by calculating their correlation coefficient. Computation of radar echo signal coherence typically involves spatial averaging of the radar echo signals within a moving window. Decorrelation, numerically equal to 1 minus the coherence coefficient, denotes the loss of coherence. Various factors contribute to decorrelation, including temporal, thermal, and spatial. Estimating temporal decorrelation involves eliminating the impacts of thermal and spatial decorrelation from the total decorrelation [27].

In 1992, Zebker and Villasenor studied the relationship between ITD and surface erosion of the Death Valley in California, and the results showed that ITD and wind erosion degrees have a negative linear correlation. In another study in 2000, Wegmuller confirmed this relationship, and he pointed out that the above relationship is still valid when the vegetation coverage is less than 40% [28]. Later, InSAR decorrelation was used to study the desertification due to wind erosion, also showing that this technology has the ability to detect ground changes due to wind erosion [29–31]. In some other studies related to dune stability, it has been confirmed that temporal decorrelation technology can detect changes in the surface of dunes due to wind erosion [32,33].

Although many indicators in the current study were proposed to describe the WEI, few can be used to describe and measure the WEI of a resolution unit accurately. This paper proposed a novel WEI of a resolution unit based on the surface deformation caused by wind erosion. Zebker’s research shows that the ITD in the resolution unit is related to the displacement of scatters in it [26]. Based on the relationship between scatters’ random displacement and the surface random deformation within a resolution unit, we related the WEI with the SITD of a resolution unit. In fact, the resolution unit’ ITD often involves the contribution of soil and vegetation in the bare lands of the Aral Sea. Wegmuller’s study indicates that the ITD of a resolution unit was equal to the sum of the weighted ITD of all scatterers within the resolution unit [28]. However, few studies are involved in decomposing the ITD contributions of all scatterers within a resolution unit. Therefore, this paper focuses on the:

• Model the WEI of a resolution unit and relate the WEI with the SITD of this resolution unit.
• SITD estimation within a resolution unit.
• Mapping the WEI of the dry lakebed of the Aral Sea.

2. Study Area

Part of the Aral Sea is in southern Kazakhstan, and the rest is in northern Uzbekistan. Due to the extended distance from the sea, the Aral Sea has a classical continental dry climate. The study area is shown in Figure 1.

Figure 1. Study area and the corresponding land types. The 2017 land surface coverage map with a resolution of 10 m is from Tsinghua University (http://data.ess.tsinghua.edu.cn/).

Figure 1. Study area and the corresponding land types. The 2017 land surface coverage map with a resolution of 10 m is from Tsinghua University (http://data.ess.tsinghua.edu.cn/).

In the Aral Sea, the maximum temperature difference between spring and summer can reach 60 °C, and the average annual precipitation is about 100mm, which is very rare. The significant temperature difference and scarce rainfall are very conducive to wind erosion. The significant temperature difference and scarce rainfall are very conducive to wind erosion. Furthermore, as shown in Figure 1, this area is dominated by bare land. Besides, the desertification of the dry lakebed is badly severe, and the loose soil is vulnerable to erosion. Thus, the dry climate, sparse vegetation coverage, and loose soil properties make this area a potential wind erosion area. Sampling data was from a field survey about the Aral Sea’s desertification in November 2018, and parts of the sampling sites are listed in Figure 1.

3. Method and Data

3.1. Soil Sampling and Volumetric Soil Moisture Data

We obtained soil data at the sampling point in this field survey, including soil photos and soil salinity data. The landscape photos can visually indicate whether the soil is vulnerable to wind erosion. Therefore, we can use these data to qualitatively assess the ITDDM model and the estimation results of SITD. The volumetric soil moisture (topsoil from 0 to 7 cm) data comes from the EAR5 data set on the Google Earth Engine platform and will be used to extract potential dust emission areas.

3.2. Vegetation Fraction Coverage Data

Vegetation fraction coverage (VFC) is related to the vegetation and soil’s backscattering weights under the assumption that only these two land types are in this resolution unit [34]. We can use VFC to estimate the soil microwave backscattring coefficient (SMBC) and soil microwave backscattring coefficient (VMBC) related to the ITD weights of vegetation and soil. In this paper, we used NDVI of Landsat 8 to compute VFC based on the Pixel dichotomy [35], and the spatial resolution of NDVI data was 30 m.

3.3. Microwave Backscattering Data

Sentinel-1 dataset C-band Synthetic Aperture Radar (SAR) dataset will be used to estimate SMBC of VMBC of a resolution unit. The spatial resolution of this data is about 10 m. Usually, the normalized microwave backscattring coefficient (NMBC) is between 0 and 1. The NMBC and the microwave backscattring coefficient (MBC) described in dB have a relationship like this [28]:

$${\sigma }_{dB}=10{l\mathrm{o}\mathrm{g}}_{10}{\sigma }_N,$$

(1)

Where ${\sigma }_{dB}$ denotes the MBC in dB, and ${\sigma }_N$ is NMBC. Due to the very sparse vegetation coverage in the arid and semi-arid area, the volume scattering is very weak, and thus we use the ITD of VV polarization to characterize wind erosion of topsoil. The ITD weights of vegetation and soil are related to the SMBC and VMBC in the resolution unit. Therefore, the SMBC and VMBC in a resolution unit will be estimated firstly by the VFC and the total MBC. The backscattering images with VV polarization taken on June 23, 2020, and July 05, 2020 (https://code.earthengine.google.com/), were obtained, and this time interval is the same as that for calculating the VFC.

3.4. Soil Sampling Data

We used the Sentinel-1A satellite C-band single look complex (SLC) SAR images to calculate the ITD of a resolution unit. In this study, the satellite mode was right-looking, and the orbit cycle was 12 days. We acquired two SLC image pairs with a descending strip-map pattern on June 23, 2020, and July 05, 2020. The incident angle was about 34.23°, and the spatial resolution was 20 m after multi look processing. Other information of two SLC image pairs is shown in Table 1.

Table 1. Information and interferometry pattern of two SLC image pairs.

Footprint	Acquisition date	Obit number	Combination mode	Time baseline	Normal baseline
North Aral Sea	2020/06/23	33138	master	12 days	-33.064 m
	2020/07/05	33313	slaver
Aral Sea South	2020/06/23	33138	master	12 days	-35.154 m
	2020/07/05	33313	slave

* The Sentinel-1A data is from the copernicus dataspace (https://dataspace.copernicus.eu/).

Table 1. Information and interferometry pattern of two SLC image pairs.

Footprint	Acquisition date	Obit number	Combination mode	Time baseline	Normal baseline
North Aral Sea	2020/06/23	33138	master	12 days	-33.064 m
	2020/07/05	33313	slaver
Aral Sea South	2020/06/23	33138	master	12 days	-35.154 m
	2020/07/05	33313	slave

* The Sentinel-1A data is from the copernicus dataspace (https://dataspace.copernicus.eu/).

As is shown in Table. 1, the SLC images acquired on June 23, 2020, were set as the master images, and the others were set as slave images. The absolute time baseline of these two SLC image pairs is 12 days, and the normal baseline is about -33.064 m and -33.154 m for the SLC image pairs of the study area, respectively. The relatively short time baseline setup mainly meets the assumption of identical backscattering levels for acquisitions 1 and 2 [28]. The short spatial baseline is much smaller than the critical baseline, indicating that the influence of the spatial baseline on ITD is negligible [26].

3.5. Method

3.5.1. Wind Erosion Intensity Modeling

Generally, the areas where wind erosion occurs within a resolution unit are randomly distributed. Therefore, the degree of wind erosion of the 1 - D soil profile can help us analyze and define a resolution unit’s degree of wind erosion. Figure 2 shows the deformation of the soil surface in the resolution unit after wind erosion.

The ground surface will sink slightly and randomly after wind erosion for most arid bare lands, except for desert regions. Due to particles’ deposition or surface movement, the ground surface may slightly uplift after wind erosion in desert areas. Fortunately, there are few active deserts in the study area. The deformation of different positions in the resolution unit can be regarded as a random signal $\delta \left(x,y\right)$ at the point $\left(x,y\right)$. Since the root mean square can describe the intensity of the random signal, we can use the root mean square of the erosion depth at different locations within the resolution unit to define the WEI.

$$WEI=\sqrt{\frac{1}{S}\iint {\delta \left(x,y\right)}^{2}dxdy},$$

(1)

In order to simplify this process, we only consider the discrete case. Assuming that there are N differential elements in a resolution unit, the wind erosion intensity of the resolution unit can be simplified as:

$$WEI=\sqrt{\frac{1}{Nds}{\sum }_{i=1}^N{\delta }_i^{2}ds}=\sqrt{\frac{1}{N}{\sum }_{i=1}^N{\delta }_i^2},$$

(2)

3.5.2. InSAR Temporal Decorrelation

Coherence is used to measure the similarity of the InSAR echoes and is often used to describe different degrees of terrain changes. The coherence between two complex echo signals ($s_1$ and $s_2$) is defined as their correlation coefficient $\gamma$.

$$\gamma =E\left\{s_1·s_2^{*}\right\}/\sqrt{E\left\{{\left|s_1\right|}^2\right\}·E\left\{{\left|s_2\right|}^2\right\}},$$

(3)

Decorrelation $\rho$, which is equal to $1-\gamma$, is usually caused by thermal decorrelation ${\rho }_{thermal}$, spatial decorrelation ${\rho }_{spatial}$, and the temporal decorrelation ${\rho }_{temporal}$. Temporal decorrelation due to environmental changes over time can be described by the following formula [26].

$${\rho }_{temporal}=\rho /{\rho }_{thermal}/{\rho }_{spatial},$$

(4)

3.5.3. The Ralationship between InSAR Temporal Decorrelation and Wind Erosion Intensity

According to the microwave backscattering theory, the differential element’s random subsidence in the resolution unit can be equivalent to the random displacement of the scatterers (or scattering differential bins) in the vertical direction (Figure 3).

Zebker derived the relationship between the root mean square (RMS) displacement of scatterers and ITD [26]:

$${\rho }_{temporal}=exp\left\{-\frac{1}{2}{\left(\frac{4\pi }{\lambda }\right)}^2\left({\delta }_y^2{sin}^2\theta +{\delta }_z^2{cos}^2\theta \right)\right\},$$

(5)

Where ${\rho }_{temporal}$ is the ITD for scatterers within a resolution unit, ${\delta }_y$ denotes the horizontal displacement of the scatterer, ${\delta }_z$ denotes the vertical displacement of the scatterer, and $\theta$ is the incident angle. Since wind erosion in a resolution unit can be equivalent to the vertical displacement of the scatterers within the resolution unit, Equation (5) can be simplified to:

$${\rho }_{temporal}=exp\left\{-\frac{1}{2}{\left(\frac{4\pi }{\lambda }\right)}^2{\delta }_z^2{cos}^2\theta \right\},$$

(6)

However, it must be noted that the vertical displacement of the scatterers may be positive or negative because of the existence of horizontal displacement of the scatterers or the deposition of the soil particles. In fact, except for the more mobile deserts, most of the bare land in the Aral Sea exhibited random subsidence rather than uplift under wind erosion. Furthermore, the ground uplift caused by the particles’ horizontal displacement and the particles’ deposition is so tiny that it can be ignored. We plotted the relationship between the ITD of the C-band and the wind erosion intensity when the incident angle was 34°, and the results are shown in Figure 4.

According to Figure 4 and formula (6), ITD and WEI have a non-linear negative correlation, which coincides well with Wegmuller’s field survey [28]. However, in bare land of the Aral Sea, a resolution unit usually contains soil and vegetation. Therefore, we must decompose these contributions and estimate the SITD, which can be used to describe the wind erosion degree of soils.

3.5.4. InSAR temporal Decorrelation Decomposition Model

Usually, there is sparse vegetation on the arid bare lands. According to Wegmuller, the ITD of a resolution unit can be decomposed into the contributions of soil and vegetation when there is only soil and vegetation within it [28].

$$\mathsf{\gamma }=\frac{{\sigma }_v}{\sigma }{\gamma }_v+\frac{{\sigma }_s}{\sigma }{\gamma }_s,$$

(7)

${\sigma }_v$, ${\sigma }_s$, and $\mathsf{\sigma }$ are the VMBC, SMBC and total MBC of this resolution unit, respectively. ${\gamma }_s$, ${\gamma }_v$ and $\mathsf{\gamma }$ denote the SITD, vegetation ITD (VITD) and total ITD of this resolution unit, respectively. However, before formula 7 can be used to estimate the SITD, we must know the backscattering coefficient of soils and vegetation of a resolution unit. Therefore, a backscattering contribution decomposition (MBCD) model within a resolution unit proposed in our previous work will be used to unmix the backscattering contributions of vegetation and soil and estimate the VMBC and SMBC within a resolution unit [34].

3.5.5. Backscattering Contribution Decomposition and Estimation Within a Resolution Unit

In our previous work about the backscattering contribution decomposition within a resolution unit, we have reached an important conclusion: the SMBC, VMBC, and total MBC of a resolution unit satisfy a simple linear relationship. This relationship can be expressed as [34]:

$${\sigma }_{veg}{f}_{veg}+{\sigma }_{soil}\left(1-f_{veg}\right)=\sigma ,$$

(8)

${\sigma }_{veg}$, ${\sigma }_{soil}$ and $\sigma$ are the VMBC, SMBC and total MBC of a resolution unit, respectively. $f_{veg}$ is the VFC of this resolution unit. Our previous work also provided the estimation method of VMBC and SMBC, which will be used to estimate the SITD [34].

3.5.6. SITD Estimation Based on LSM-SVD Method

The continuity of the spatial distribution of soil and vegetation allows them to have nearly the same temporal decorrelation in adjacent resolution units. Thus, the temporal decorrelation of a resolution unit can be estimated by the sample points within the buffer of this resolution unit. If n sample points around the resolution unit (resolution P) were used to estimate the ITD of vegetation and soil of this resolution unit, n linear equations can be written as:

$$\begin{gathered} w_v^1{\gamma }_v+w_s^1{\gamma }_s={\gamma }_1 \\ {w}_v^2{\gamma }_v+w_s^2{\gamma }_s={\gamma }_2 \\ ⋮ \\ {w}_v^i{\gamma }_v+w_s^i{\gamma }_s={\gamma }_i \\ ⋮ \\ {w}_v^n{\gamma }_v+w_s^n{\gamma }_s={\gamma }_n, \end{gathered}$$

(9)

Where $w_v^i$ and $w_s^i$ are the ITD weights of vegetation and soil of the $ith$ sampling point, ${\gamma }_v$ and ${\gamma }_s$ denote VITD and SITD of the resolution unit P, respectively, and ${\gamma }_i$ is the composite ITD of the $i\mathrm{t}\mathrm{h}$ sampling point. The system equations (formula (9)) can also be written in the form of a matrix:

$$\left[\begin{array}{cc}\begin{array}{c}{w}_v^1\\ \begin{array}{c}{w}_v^2\\ ⋮\\ \begin{array}{c}{w}_v^i\\ \begin{array}{c}⋮\\ w_v^n\end{array}\end{array}\end{array}\end{array}& \begin{array}{c}{w}_s^1\\ \begin{array}{c}{w}_s^2\\ \begin{array}{c}⋮\\ \begin{array}{c}{w}_s^i\\ \begin{array}{c}⋮\\ w_s^n\end{array}\end{array}\end{array}\end{array}\end{array}\end{array}\right]\left[\begin{array}{c}{\gamma }_v\\ {\gamma }_s\end{array}\right]=\left[\begin{array}{c}{\gamma }_1\\ \begin{array}{c}{\gamma }_2\\ ⋮\end{array}\\ \begin{array}{c}{\gamma }_i\\ ⋮\\ {\gamma }_n\end{array}\end{array}\right],$$

(10)

Let $\mathbf{w}$ be the ITD weight matrix, ${\mathsf{\gamma }}_{\mathbf{p}}$ be the InSAR temporal decorrelation decomposition matrix of the resolution unit P, $\mathsf{\gamma }$ be the InSAR temporal decorrelation composite matrix, and make the following conventions:

$$\mathbf{w}=\left[\begin{array}{cc}\begin{array}{c}{w}_v^1\\ \begin{array}{c}{w}_v^2\\ ⋮\\ \begin{array}{c}{w}_v^i\\ \begin{array}{c}⋮\\ w_v^n\end{array}\end{array}\end{array}\end{array}& \begin{array}{c}{w}_s^1\\ \begin{array}{c}{w}_s^2\\ \begin{array}{c}⋮\\ \begin{array}{c}{w}_s^i\\ \begin{array}{c}⋮\\ w_s^n\end{array}\end{array}\end{array}\end{array}\end{array}\end{array}\right],{\mathsf{\gamma }}_{\mathbf{p}}=\left[\begin{array}{c}{\gamma }_v\\ {\gamma }_s\end{array}\right],\mathsf{\gamma }=\left[\begin{array}{c}{\gamma }_1\\ \begin{array}{c}{\gamma }_2\\ ⋮\end{array}\\ \begin{array}{c}{\gamma }_i\\ ⋮\\ {\gamma }_n\end{array}\end{array}\right],$$

(11)

Then formula (9) can be can be simplified as:

$$\mathbf{w}{\mathsf{\gamma }}_{\mathbf{p}}=\mathsf{\gamma },$$

(12)

Then the least squares estimation of ${\mathsf{\gamma }}_{\mathbf{p}}$ can be expressed as:

$$\widehat{{\mathsf{\gamma }}_{\mathbf{p}}}={\left({\mathbf{w}}^{\mathrm{T}}\mathbf{w}\right)}^{-1}{\mathbf{w}}^{\mathrm{T}}\mathsf{\gamma },$$

(13)

Considering the huge computational load, we can use the SVD method to estimate SITD. Assume that $\mathbf{w}$ can be decomposed as follows:

$$\mathbf{w}=\mathbf{M}\mathbf{S}{\mathsf{{\rm E}}}^{\mathrm{T}},$$

(14)

The SVD-based least square estimation of ${\mathsf{\gamma }}_{\mathbf{p}}$ can be expressed as:

$$\widehat{{\mathsf{\gamma }}_{\mathbf{p}}}={\mathbf{E}}_{2\times r}{\left({\mathbf{S}}_{r\times r}\right)}^{-1}{\left({\mathbf{M}}_{n\times r}\right)}^{\mathrm{T}}{\mathsf{\gamma }}_{n\times 1},$$

(15)

When the first singular value of $\mathbf{w}$ is greater than or equal to 90% of the sum of all its singular values, the value of $r$ is 1, and $\mathbf{S}$ has the following form:

$$\mathbf{S}={ϵ}_1,$$

(16)

Where ${ϵ}_1$ is the first singular value of $\mathbf{w}$. Otherwise, S has the following form:

$$\mathbf{S}=\left[\begin{array}{cc}{ϵ}_1& 0\\ 0& {ϵ}_2\end{array}\right]$$

(17)

Where ${ϵ}_2$ is the second singular value of $\mathbf{w}$.

4. Results

4.1. VFC and MBC of the Study Area

The VFC with a resolution of 30 m calculated by Landsat 8 NDVI and the Sentinel-1A C-band microwave backscattering coefficient of VV polarization with a resolution of 10 m from the Google Earth Engine platform were used to estimate the VMBC and SMBC. To evaluate the vegetation’s influence on the SITD estimation, we also calculated the SMBC to VMBC (SDV) ratio. Suppose the total ITD of the resolution unit replaces the SITD of a resolution unit. In that case, the error mainly depends on the VMBC and SMBC in the resolution unit. The VFC, total MBC, SMBC, VMBC, and SDV are shown in Figure 5. As shown in Figure 5(a), Vegetation coverage in the Aral Sea and its surrounding areas is very sparse. The vegetation coverage of most of the study areas is between 0 and 0.1608. In comparison, the vegetation coverage is relatively high in the northwest and northeast (VFC ranges from 0.1608 to 0.3137) and the southern part of the study area (VFC is between 0.3137 and 1).

It can be seen from Figure 5(b) that, except for the area around the east branch of the South Aral Sea, the backscattering coefficients in other regions are relatively low. Figure 5(c) shows that the soil’s backscattering coefficient is very close to the total backscattering coefficient of the resolution unit. However, there are obvious errors in estimating soil backscattering coefficients in the northwest and northeast of the study area, most probably due to the uniform spatial distribution of vegetation [34]. In addition, the results of backscattering coefficient estimation have been fully verified in our previous work, so this paper will not refer to the verification of backscattering coefficient estimation results [34]. Furthermore, the areas with incorrect backscattering coefficient estimation results often have relatively high soil water content or vegetation coverage, so we can remove these places from the study area, which will be discussed in Section 4.2. Figure 5(c) shows that the backscattering coefficient of soil is much higher for most of the study area than that of vegetation. According to formula (10) and SDV shown in Figure 5, the influence of vegetation on SITD in the resolution unit may be negligible. However, when the SMBC in a resolution unit is close to the VMBC (such as in severely desertified areas), the ITD weights of vegetation and soil are almost the same. In this case, the impact of vegetation on SITD is not negligible. The process of MBC and ITD estimation is badly time-consuming. Extracting the potential wind erosion area can significantly reduce the target area and effectively reduce the amount of calculation. Numerous studies have shown that wind erosion is almost impossible in these areas when the soil moisture exceeds 10%, or the VFC is higher than 40% [36,37,38,39,40,41,42]. When the vegetation coverage exceeds 40%, the volumetric soil moisture is usually higher than 10%. Our survey results also show that when the vegetation coverage exceeds 40%, the soil volumetric water content is usually higher than 10% [34]. Therefore, these areas can be removed from the study area to reduce the amount of calculation effectively. In the following section, we obtained soil moisture data. Soil moisture and vegetation coverage data were used to extract potential wind erosion areas.

4.2. Potential Wind Erosion Areas

We used GIS software to label the areas with vegetation coverage higher than 0.4 or volumetric soil moisture higher than 0.1 in the study area. Then we got the potential wind erosion area, and the result is shown in Figure 6. It should be noted that since the volumetric soil moisture (VSM) of only one pixel in area B exceeds 10%, we have ignored this area. Figure 6 shows that except for areas B, E, and F, the rest of the lakebed is vulnerable to wind erosion.

4.3. InSAR Temporal Decorrelation of Soil

We calculated the SITD of all regions except regions A, B, D, E, and F based on Equation 15. The process of SITD estimation was like that of MBC. The upper and lower thresholds of the VFC difference are set to 0.2 and 0.05, respectively, and the buffer size is set to 100 meters in this paper. Even so, the calculation process is still quite time-consuming. Therefore, we resample the ITD weights of soil and vegetation and the total ITD to 50 meters to further reduce the amount of calculation. The estimated result of SITD is shown in Figure 7. Figure 7 shows that the areas with severe temporal decorrelation are mainly distributed in the bare lands of the middle, northeast, and southeast of the South Aral Sea. In addition, the coastal areas of the North Aral Sea are also areas with severe soil InSAR temporal decorrelation. Beyond the boundaries of the water body of the Aral Sea in 1973, the soil InSAR temporal decorrelation in most areas is very slight. This significant contrast indicates that the dry lakebed is highly susceptible to wind erosion.

The pixel-scale SITD and the corresponding landscape photos are shown in Figure 8.

The landscape photos in Figure 8 clearly show the soil’s actual condition near the sampling points, which can be used to visually judge the degree of soil susceptibility to wind erosion. The SITD will be used to describe the spatial distribution of wind erosion of different degrees. Like the description in part 3.3, because the waters and wetland cannot be the potential place where wind erosion will happen, the SITD was not be estimated for these two land types. In the next section, SITD will be used to describe the WEI of the study area.

4.4. Wind Erosion at the Dry Bottom of the Aral Sea

Do an inverse transform to Equation (6), and then we can convert SITD to WEI. Furthermore, we divided WEI into eight levels, and the results are shown in Figure 9. Figure 9 shows that the wind erosion in the study area mainly occurs on the lakebed. Besides, our previous field survey showed that the dry lakebed is rich in salt and toxic substances [34]. Therefore, the fine particles produced in the wind erosion process should also be rich in these materials.

Salt storms that carry fine particles rich in salt and toxic substances will badly threaten the surrounding ecosystems and human health. Therefore, we simply count the wind erosion area of the dry lakebed. We used GIS software to remove areas with VFC greater than or equal to 4, areas with VSM greater than or equal to 0.1, wetlands, and water bodies from the Aral Sea water surface in 1973. Then, we used this vector to extract the WEI map of the potential wind erosion area on the dry bottom of the Aral Sea. At last, we performed simple statistics on the wind erosion intensity, and the results are shown in Table 2.

5. Discussion

Due to the harsh natural environmental conditions of the Aral Sea, it is a challenge to verify the accuracy of dust activity intensity described by soil temporal decorrelation through the measurement of dust activity intensity. The sampling data for validation is from a joint desertification survey referring to the Aral Sea in 2018 by China and Uzbekistan. However, the survey was not designed specifically for this research, so we can only provide limited validation by analyzing the partial sampling data. Nevertheless, through rigorous theoretical derivation, little sampling points, and supporting literature, we can still prove the accuracy of the spatiotemporal distribution results of dust activity intensity described in this paper. The accuracy of the spatiotemporal distribution results of dust activity intensity in the Aral Sea depends on three aspects:

1. Can soil temporal decorrelation accurately describe dust activity intensity?

2. Whether soil temporal decorrelation only depends on the mathematical expectation of phase random variation within the resolution unit and whether other factors, such as variation of soil dielectric constant and soil roughness, will affect soil temporal decorrelation.

3. Are the estimation results of soil temporal decorrelation accurate?

For the first question, we have provided a comprehensive argumentation through the proposal of dust activity intensity in Section 3.5.2 and the derivation of the relationship between dust activity intensity and soil temporal decorrelation in Section 3.5.3, which demonstrates the feasibility of using soil temporal decorrelation to describe dust activity intensity. Next, we will discuss the impact of the other two factors.

5.1. The Impact of Non-Phase Factors on Soil Temporal Decorrelation

According to the definition of soil temporal decorrelation, it mainly depends on the variations in the backscattering coefficient and phase during the interferometric period. Significant changes in soil roughness, moisture content, and salinity can all lead to notable variations in backscattering levels. Through an investigation of precipitation during the interferometric period, no significant rainfall was observed in these regions within 1-2 days before the second imaging (in June and July, the temperature in the Aral Sea is high, and any small amount of rainfall occurring a few days before the second imaging would quickly evaporate). This survey indicates that the soil moisture content in most areas is unlikely to undergo significant changes, thus not significantly influencing the backscattering coefficient. The soil salinity of the Aral Sea primarily comes from the evaporated seawater, and after the seawater dries up, the salt deposits in the soil at the bottom of the dried-up lake. Therefore, the soil salinity is unlikely to undergo significant changes over a relatively short time interval (12 days), and thus, it cannot cause drastic variations in the backscattering coefficient. The spatial continuity of soil types within the resolution unit and similar wind conditions within the unit ensure slight variation in erosion levels among different regions within the unit. Consequently, the surface roughness variation is relatively minor. Therefore, the soil temporal decorrelation of a resolution unit primarily depends on the mathematical expectation of phase random variations at different locations within this resolution unit, representing dust activity intensity, and these phase variations are mainly caused by wind erosion.

5.2. The Results of SITD Estimation Assessment

The accuracy of soil temporal decorrelation estimation primarily depends on the assumptions made during the decomposition and estimation. According to Wegmuller’s research, two assumptions regarding the soil temporal decorrelation decomposition model are often valid for bare soil and sparse low vegetation regions in arid and semi-arid areas. The estimation was based on the assumption that "spatial adjacent resolution units have the same soil and vegetation temporal decorrelation." Soil and vegetation’s temporal decorrelation primarily depend on all scatterers’ root mean square displacement within a resolution unit. Spatial adjacent resolution units have similar vegetation types and structures, soil types and structures, and highly similar wind conditions. Therefore, the root mean square of soil random erosion of a resolution unit is almost the same as that of an adjacent resolution unit. Similarly, a resolution unit’s vegetation temporal decorrelation is nearly the same as that of an adjoining resolution unit. Hence, the assumption that a resolution unit’s soil or vegetation temporal decorrelation is almost the same as that of a spatial adjacent resolution unit is reasonable.

The differences in the backscattering coefficients of soil and vegetation among adjacent resolution units are the prerequisite for estimating the soil and vegetation temporal decorrelation. Because of the slight differences in soil dielectric constant (due to spatial variations in soil moisture and salinity), soil roughness (due to the spatial difference of wind erosion depth within a resolution units), vegetation dielectric constant (due to spatial variations in vegetation moisture content, vegetation type, vegetation density, health, and condition), and vegetation roughness (due to the spatial distribution variations in leaf structure, branching patterns, vegetation canopy height and density, vegetation growth stage, and environmental factors.), the backscattering coefficients of soil and vegetation in adjacent resolution units often exhibit differences. Consequently, it is possible to accurately estimate the soil temporal decorrelation of a resolution unit.

As described in Section 4.3, we use the SITD maps (SITDMs) and the corresponding landscape photos to verify the estimation results of SITD. As is shown in landscape photos of L5, L6 and L15, some parts of arid soil are covered by relatively dense vegetation, while some are covered by very sparse vegetation. Besides, these two landscape photos also showed that the soil is relatively loose here. According to the study of Bagnold, the ground surface wind speed will significantly decrease when the vegetation exists, and the degree of wind erosion will also be declined [4]. So, the wind erosion of soils around L5, L6 and L15 with few vegetation coverages is severe but is relatively slight for soils around L5, L6 and L15 with relatively dense vegetation coverage, and the wind erosion of soils shown in SITDMs coincides well with the actual soil property shown in the corresponding landscape photos. Although there is few vegetation at site L7, the soil was very stable according to the ITDM of L7. Precipitation around L7 between June 22, 2020, and July 05, 2020, is a possible cause of the stable status of soils. The relatively stable status of L8, L9 and L11 shown in their corresponding SITDMs is most probably the result of dense vegetation coverage showed in their corresponding landscape photos. The relatively slight wind erosion of L10 (SITDM of L10) may be due to soil’s high-water content, which is most probably the result of the short distance between L10 and water. As shown in the landscape photos of site L12, the desertification trend is severe. However, the SITDM of L12 indicates the soil here is very stable, and the relatively stable status of the soil is most probably due to the rain between June 23, 2020, and July 05, 2020, which can reduce the wind erosion significantly [33]. In Figure 8, the soil property of L13 is just the same as that of L12, and the corresponding SITDM indicates the wind erosion here is very severe, so the actual soil property shown in the landscape photo is also in good agreement with the wind erosion degree shown in the corresponding SITDM. The landscape photos of L14 and U7 shows that there is few vegetation at these places, but this place is rich in relatively large stones. The large stones are hard to be moved by wind, and meanwhile, they can reduce the wind speed on the ground surface, which made this place remarkably stable. For L16, as is shown in the landscape photos of L16, there are many artificial buildings and roads here, the SITD and backscattering coefficient of which are much higher than that of vegetation [32]. Therefore, the SITD should be very high for the resolution unit with buildings in it according to formula (10). The high SITD of L16 showed in its SITDM coincides nicely with the soil property shown in the landscape photo of L16. Large areas of bare soil were found around U5 and U6, and there are few large stones here, so U5 and U6 are very vulnerable to soil erosion. The analysis above is consistent with the severity of wind erosion shown in the SITDMs of U5 and U6. In conclusion, the SITDMs, except for some influenced by precipitation, are all consistent with the actual soil property shown in the corresponding landscape photos, and this result indicates that SITD can be used to describe the severity of wind erosion.

The 2008 NOAA AVHRR images (Figure 10) show that the dust emission areas are mainly distributed in the bare land between the two branches of the South Aral Sea, the northeast and southeast of the South Aral Sea [43]. These areas spatially coincide well with the results of this study.

Besides, investigation regarding the dust storm events in the Aral Sea (2005-2008) also indicates that the northeastern and southeastern parts of the South Aral Sea are the right places where dust storms often occur [43]. Although sandstorm events also happened in the bare land between the two branches of the South Aral Sea, the frequency is much lower than that of the northeast and southeast of the South Aral Sea. Therefore, the locations of dust emission identified in this study, except for the bare land between the two branches of the South Aral Sea, are consistent with the investigation results. A plausible explanation is that desertification had already occurred in the bare land between the two branches of the South Aral Sea in 2008 or even before, but the desertification was mild. After about 12 years of desertification, this region’s desertification becomes severe.

Furthermore, the dust activity intensity of this region described in this paper is confirmed by the landscape photographs of the three sampling points falling in this area (Figure 11). Therefore, it is highly likely that the bare land between the two branches of the South Aral Sea has developed into a new outbreak area for salt and dust storms.

5.3. Wind Erosion at the Dry Bottom of the Aral Sea

(1) Overview of the wind erosion in the Aral Sea

We conducted a simple statistic of different degrees of the WEI according to Figure 9, and the results showed that the dry lakebed is suffering from different degrees of wind erosion. The area with wind erosion intensity greater than 1.5 cm accounts for 0.2790% of the lakebed, which is about 176 square kilometers. The area with wind erosion intensity ranging from 1 to 1.5 cm accounts for 2.2560%, approximately 1,426 square kilometers. The site with a wind erosion intensity of 0.5 to 1 cm accounts for 35.5420%, approximately 22,471 square kilometers. The area with a wind erosion intensity of 0.4 to 0.5 cm accounts for 15.3211%, approximately 9,687 square kilometers. The site with wind erosion intensity between 0.3 and 0.4 cm accounts for 15.4784%, about 9,786 square kilometers. The area with a wind erosion intensity of 0.2 to 0.3 cm accounts for 27.2261%, approximately 17,214 square kilometers. The site with wind erosion intensity ranging from 0.1 to 0.2 cm accounts for 3.723%, about 2,354 square kilometers. The area with wind erosion intensity between 0 and 0.1 cm occupies 0.1736%, approximately 110 square kilometers.

(2) Spatial distribution of wind erosion

According to the WEI of the study area, the severe wind erosion areas ($\mathrm{W}\mathrm{E}\mathrm{I}\ge 1 \mathrm{c}\mathrm{m}$) are mainly distributed in bare lands of the middle, the northeast, and southeast of the South Aral Sea. The spatial distribution of the moderate wind erosion area ($0.3 \mathrm{c}\mathrm{m}\le \mathrm{W}\mathrm{E}\mathrm{I}\le 1 \mathrm{c}\mathrm{m}$) is the same as that of the severe wind erosion area. The moderate and the severe wind erosion area coincide well with the area of sand and dust activities, and these areas can be treated as the focus areas for people to intervene in sand and dust activities [43]. The severe wind erosion areas within the dry lakebed are the most probably where salt dust storms occur.