2.1 Study area
The Dongzhiyuan plateau spans between 106°14’-108°42’E,34°50’-37°19’N, which corresponds to the southeastern edge of the Gansu Province, China (Fig.1a), an arid and semi-arid zone. The topography of the study area is dominated by gully-incised tablelands, on terrain that displays an elevation ranging between 882-1540 m.a.s.l with an average of 1300 m.a.s.l (Fig.1b). The entire Dongzhiyuan covers an area of 2835 km2, while the tableland area is 960 km2. Its average annual temperature is 8.5°C, whilst its average annual precipitation ranges between 480-660 mm. Nevertheless, its average annual evaporation is 1600 mm (Xu et al. 2018). However, trends are not uniform either spatially or over time. Rainfall in the south of the plateau is larger than that in the north, and at least 70% of the yearly value occurs between July and September. Furthermore, storms tend to be high-intensity and short-duration (Jiang et al. 2021a). As expected, vegetation cover mirrors geographical trends on rainfall, spices adapted to arid terrain are predominant in the north, whilst lush, more perennial vegetation is widespread in the south (Ling et al. 2022). The water table exhibits extensive variation due to the high time-clustering of rainfall. Particularly, it is usually deep, as evaporation exceeds ground percolation during most of the year, leading to low water content in the surface soil (Li et al. 2014). Top layers (150 to 200m thick) in the plateau are comprised mostly of Quaternary Upper Pleistocene loess that displays an on-site density of 1.79-1.97g⋅cm-3. Loess depth generally exceeds 100m. Moreover, due to their loose nature, they are susceptible to collapse, leading to the formation of cavities and voids within them. The loess surficial strata are underlaid by a black kiln clay soil. This geomorphological arrangement is observed all over the plateau (Zhang et al. 2021). It must be stressed that these geomorphological units showcase different erosion-resistant traits. The black kiln layer can retain moisture and soil particles, while the loess tends to disaggregate, thus making it susceptible to being washed away by rainwater, as shown in Figs.1e, d.
2.2 Data source
2.2.1 Generation of basic data with the aid of GIS
Diverse data sources were tapped for collecting the information required by this study. Firstly, Google Earth Pro was employed to collect images in the visual spectrum (0.4 to 0.7 mm) to get an overview of the study site and select features of interest. Furthermore, a Digital Elevation Model (DEM) with a 90m resolution was sourced for U.S National Aeronautics and Space Administration (NASA) (http://www.gscloud.cn/search). It was useful for performing InSAR analyses for orbital interference de-weighting and to provide ancillary data for geocoding. Both sources were integrated using ArcGIS 10.7, which was instrumental for mapping the study area, identifying gullies that were vectorized as lines, and allowing for the computation of geomorphological parameters using elevation data. Then, SAR images from the Sentinel satellite constellation were collected through the Copernicus hub (https://scihub.copernicus.eu/) to estimate displacement time histories using SBAS-InSAR. Further details about these datasets are presented in the following section. Obtained displacement data were interpolated to match the resolution of the DEM, allowing for exploration of dependency traits among surficial deformation, slope, and absolute elevation. Finally, rainfall information was sourced from the National Weather Science Data Center of China (http://data.cma.cn/), to explore the role of precipitation in gully erosion. Particularly, records were taken from the Qingyang City rainfall station.
2.2.2 Sentinel-1A data
As stated previously, SAR images were downloaded from the Copernicus hub, supported by ESA (https://scihub.copernicus.eu/). The analysis period ranged from April 2019 to October of the same year. Images selected for the analyses were taken each 12 days, after carefully verifying all shared the same incidence angle (37°). All were taken during the descending pass or the satellite and they correspond to the Vertical transmitting, Vertical receiving (VV) polarization mode, as shown in Table 1. Moreover, the accuracy was improved by considering orbital precision files (Liu et al. 2021).
Table 1. detailed information of sentinel-1A satellite images and Shuttle Radar Topography Mission (SRTM) elevation data.
Data
|
Parameters
|
Description
|
Sentinel-1A
|
Type
|
SLC
|
Imaging mode
|
IW
|
Band and wavelength (cm)
|
C, 5.5
|
Orbit direction
|
Descending
|
Azimuth resolution (m)
|
20
|
Range resolution (m)
|
5
|
Polarization
|
VV
|
SRTM
|
Resolution (m)
|
30
|
Positioning accuracy (m)
|
20
|
Elevation accuracy (m)
|
16
|
SLC stands for single-view complex, IW stands for strip scan mode, and VV stands for polarization mode as vertical polarization.
2.3 SBAS-InSAR
2.3.1 SBAS-InSAR Monitoring Principle
SBAS-InSAR is a procedure that uses SAR images to construct surface change time histories by grouping them into master and slave sets. This allows for effective data extraction from noisy sets (with low coherence) even when using a low number of images. Particularly, the following relationship holds: (Chen et al. 2020).
$$\frac{N+1}{2}\le M\le \frac{N\left(N+1\right)}{2}$$
1
Where N is the number of remote sensing images whilst M is generated considering N + 1 SAR images \ in such a way that M takes the value of 102 possible pairs in this study.
$${\Delta }{\varphi }_{j}\left(x,r\right)=\varphi \left({t}_{B},x,r\right)-\varphi \left({t}_{A},x,r\right)$$
$$\approx {\varphi }_{def,j}\left(x,r\right)+ {\varphi }_{topo,j}\left(x,r\right)+ {\varphi }_{atm,j}\left(x,r\right)+ {\varphi }_{orb,j}\left(x,r\right)+ {\varphi }_{noise,j}\left(x,r\right)$$
2
Equation (2) represents the interferometric phase composition of image j (generated by two SAR images at times tA and tB, given that tA < tB) within the pixel element at azimuth x and r distance from the satellite x, r. φ(tA,x,r) and φ(tB,x,r) are the phase values of the interferograms at time tA and tB at pixel point (x,r). Φdef,j(x,r) represents the phase difference between tA and tB whilst Φtopo,j(x,r) denotes the topographic residual phase after extracting the reference DEM elevation data. Φatm,j(x,r) is the error caused by the atmospheric phase, whilst Φorb,j(x,r) represents the orbital error, and Φnoise,j(x,r) represents random noise. Clearly, Φatm,j(x,r) can be obtained from the interferometric phase by performing time domain high-pass filtering and spatial low-pass filtering. After removing the phase error caused by the atmosphere, the interferometric phase mainly consists of deformation and residual elevation.
After that, it is assumed that the velocity vector in a certain direction can be expressed by a linear model. To obtain the surface deformation information and the residual elevation DEM component, we can transform the Φdef,j(x,j) term in Eq. (2) as follows:
$$\text{v}= {[{v}_{1},{v}_{2},\dots ,{v}_{N}]}^{T}= {[\frac{{\varphi }_{1}}{{t}_{1}-{t}_{0}},\frac{{\varphi }_{2}-{\varphi }_{1}}{{t}_{2}-{t}_{1}}, \dots ,\frac{{\varphi }_{N}-{\varphi }_{N-1}}{{t}_{N}-{t}_{N-1}}]}^{T}$$
3
$${\varphi }_{topo,j}\left(x,r\right)= \frac{4\pi }{\lambda }\bullet \frac{{B}_{{V}_{j}}{\Delta }z}{R\bullet sin\theta }$$
4
Where refers to the radar wavelength, BVj refers to the vertical baseline of the interference relative to j, θrefers to the satellite incidence angle, and △z refers to the altitude difference. Also, it is possible to formulate the C vector as:
$$\text{C}\left(\text{M}\times 1\right)=\left[\frac{4\pi }{\lambda }\bullet \frac{{B}_{{V}_{1}}}{R\bullet sin\theta },\frac{4\pi }{\lambda }\bullet \frac{{B}_{{V}_{2}}}{R\bullet sin\theta }, \dots ,\frac{4\pi }{\lambda }\bullet \frac{{B}_{{V}_{M}}}{R\bullet sin\theta }\right]$$
5
Finally, the displacement can be computed in this manner:
$${C}_{v}+D\bullet \varDelta z=\delta \varphi \left(6\right)$$
6
Where D is an M × (N-1) size matrix of known coefficients. Thus, the deformation rate can be obtained by performing singular value decomposition (SVD) on the result of expression 6.
2.3.2 Data processing workflow
Figure 2 illustrates the flow chart of the data processing done in this study. The main steps are described as follows:
Step 1: Connection map generation: This step will pair the input 16 sentinel-1A SAR images into a large master-slave set considering a 93m spatial baseline and 220 days temporal baseline, leading to 102 derived samples, Outcomes are shown in Fig. 3.
Step 2: The set collected in step one is subjected to basic interferometric processing considering DEM data as a baseline to characterize the geometric baseline. This requires the following tasks: regional coherence generation, de-leveling, filtering, coherence coefficient generation, and phase de-entanglement (minimum cost flow (MCF) is chosen for this experiment), leading to the following output: coherence coefficient maps, interferograms, and de-entanglement maps. This step will align all images to the super master image collecting pixels where coherence is larger than 0.3, then extracts values for points in were the following criteria are satisfied: (1) low geometric deformation and temporal coherence due to the longitudinal gullies are observed, there are large slopes, significant seasonal vegetation cover, and there is the rapid erosion of the gullies (2) Ensure that enough data is kept for further analysis.
Step 3: Use ground control points (GCPs) to estimate and remove phase ramps that remain after step 2 is concluded. the GCPs should be selected at locations where there are no significant topographic features with large altitude variation, no phase transitions, and ideally are far from the study area (in this experiment, a plateau surface with relatively little deformation was considered).
Step 4: This step is the core of the SBAS rate inversion, where the minimum cost flow (MCF) is used to deconvolve data again, while a linear interpolation is performed to obtain a line of sight displacements at locations of interest. Then, atmosphere effects are filtered from the previously obtained deformation rates. Finally, obtained outcomes are projected in the vertical direction. Given the extensive collapses and landslides on both sides of the gully in the study area, the subsidence rate of the gully area calculated by us is almost equal to the erosion rate.