1 Introduction

On the Qinghai-Tibet Plateau (QTP), permafrost underlies a total area of about 1.06 × 106 km2 and is shallow as well as warm (with temperature higher than −2 °C) (Ran et al. 2018). Therefore, it is particularly sensitive and vulnerable to climate change and anthropogenic disturbances (Wu et al. 2013). The QTP is warming at a rate that is double that of the global average (0.32 °C per decade compared with the global average of 0.16 °C per decade) (Yao et al. 2019). In response to climate warming and precipitation extremes, the permafrost regions of the central QTP have undergone intensifying rates of degradation (Mu et al. 2020). The degradation is mainly manifested as a reduction in the permafrost distribution area, the thinning of permafrost thickness, the thickening of the active layer, and the melting of ground ice (Yang et al. 2010).

Rapid thawing of ice-rich permafrost increases hillslope susceptibility to sliding by altering the geotechnical properties of hillslope materials, causing the release and mobilization of the stored sediment (Patton et al. 2019; Zhang, Li et al. 2022, Zhang, Lin et al. 2022), leading to the occurrence of thermokarst landslides. This type of landslide refers to landscapes with a variety of topographic depressions or collapses of unstable ground surface, including active layer detachment (ALD), thermal erosion gullies (TEG), and retrogressive thaw slumps (RTS) (Kokelj and Jorgenson 2013).

Retrogressive thaw slumps (RTSs) are prevalent slope failures related to the thawing of ice-rich permafrost and the attendant soil slumping in permafrost regions (Van Everdingen 1998; Balser 2015; Huang et al. 2020). A RTS typically consists of a sheer ice-rich headwall and a gentle slump floor occupied by mudflows (Lewkowicz 2007). Thaw slumps are widely distributed in the Beilu River Basin, Fenghuo Mountain, and Qilian Mountains and show an accelerated development trend (Luo et al. 2019). The thickness and area of individual RTSs are small, but RTSs usually occur in regional clusters and retreat at a rapid rate, posing risks to nearby vegetation, infrastructure, and communities by affecting sediment transport and topographic gradients (Runge et al. 2022). More importantly, abrupt permafrost thaw could mobilize large amounts of organic carbon through RTSs (Abbott and Jones 2015), intensifying global climate warming (Turetsky et al. 2019). For the QTP, the research about the impact of thermokarst landforms on the carbon cycle is mainly focused on RTSs (Mu et al. 2017). The comprehensive impact on landscapes and ecosystem, and the abruptness of their development make RTSs a particularly important permafrost disturbance feature that requires better investigation and prediction capabilities. Such works can provide scientific references for permafrost environmental management, ecological protection, and engineering construction.

Previous research mainly focused on the distributions, dynamic behaviors, and driving factors of RTSs in subregions of the QTP. In a study on RTS dynamic behaviors, light detection and ranging (LiDAR) and terrestrial laser scanning (TLS) were used to observe the seasonal deformation in a thermokarst gully in the northeastern QTP between 2016 and 2018 (Zhong et al. 2021). The study found that the surface subsidence and headwall retreat of the thermokarst gully reached 3.364 m and 10.66 m during this period (Zhong et al. 2021). Although these conventional methods have been proven to be suitable for monitoring the deformation of individual RTS, they are difficult to be adopted in inaccessible or remote areas with harsh natural environment and complex geological conditions (Hu et al. 2016). The high-resolution optical satellite images have the characteristics of wide coverage and rich image information. Based on manual interpretation, automatic mapping, and field investigation, a total of 438 RTSs were delineated around the Beiluhe region from optical remote sensing images in 2017 (Luo et al. 2019). Advancing these studies, more than 800 RTSs were mapped by applying a deep learning algorithm along the Qinghai-Tibet Engineering Corridor in 2019, and the number of RTSs in the Beilu River Basin has increased to 474 (Xia et al. 2022). However, when the collapse area is similar to the surrounding environment and the RTS features are not obvious, there are limitations for identifying RTSs only from optical images (Huang et al. 2020). Therefore, thermokarst landforms on the QTP are still poorly understood.

Compared with field-based or optical sensors methods, the Interferometric Synthetic Aperture Radar (InSAR) method has the characteristics of all-weather, high-resolution, high-precision, and wide coverage with low cost (Massonnet and Feigl 1998). In the last two decades, D-InSAR and multi-temporal InSAR methods, such as Persistent Scatterers interferometry (Ferretti et al. 2001) and Small Baseline Subset (Berardino et al. 2002) have been successfully used to monitor surface deformation in the permafrost regions of the QTP (Zhang, Li et al. 2022, Zhang, Lin et al. 2022). Chen et al. (2018) used the Persistent Scatterer InSAR to measure the surface subsidence of the thermokarst landforms on Eboling, and found that the subsidence trends ranged from 0.3 to 3 cm/yr between 2006 and 2011 in the thermal gully regions. Moreover, InSAR technology has been increasingly explored and successfully applied to different types of slope deformation identification and monitoring (Dini et al. 2019). However, the deformation pattern affected by decorrelation is often inconsistent with the RTS outline, which hinders the potential application of InSAR techniques (Liu et al. 2015). Compared with the counterparts in the circum-Arctic, there is still a lack of systematic studies of the characteristics of regional RTSs deformation on the QTP. Previous studies of RTS displacement measurements focused on relatively small spatial scales (Jiao et al. 2022). To achieve a better understanding of RTS dynamic variation and a more comprehensive inventory, a combination of multi-source images and different techniques is necessary.

This study aimed to characterize the spatiotemporal variations of surface deformation of retrogressive thaw slumps using a multi-temporal InSAR technique, and to assess the potential of deformation properties in mapping the retrogressive thaw slumps. We selected the Mati Hill as the study area, which is located in the Beluhe Basin on central Tibetan Plateau, because a RTS inventory in this region has been compiled from optical remote sensing images in a previous investigation (Xia et al. 2022). Based on the Sentinel-1A descending SAR images acquired between 2017 and 2021, we derived a five-year-long average deformation velocity and displacement time series of the study area, and analyzed the deformation properties of the RTS inventory. Furthermore, we compiled a new inventory of RTSs by interpreting the InSAR-identified surface displacements and high-resolution optical images. The spatial distribution characteristics of the RTSs were also revealed. Finally, the latent RTSs were detected based on multiple elements of deformation, terrain factor, and vegetation. This study conducted a novel InSAR application to detect the RTSs in the permafrost region, which can help better monitor RTS developments.

2 The Study Area

The Beilu River Basin is located in the hinterland of the Qinghai-Tibet Plateau (Fig. 1), including the terraces of the Beiluhe and Xiushui Rivers, Hoh Xil Hill to the north, and Fenghuo Mountain to the south (Luo et al. 2019). The study area is located in the Mati Hill (center location 92.71° E, 35.01° N) where thaw slumps are densely distributed. The elevation ranges from 4690 to 4910 m, with a mean of 4790 m. The region has a flat terrain with an average slope of less than 10°. Fluvial and deluvial sediments formed the upland plain landforms. The surficial sediments are dominated by fine to gravelly sands (Yin et al. 2017). The underlying substrate consists of silty clay and mudstone with sandstone interlayers (Lin et al. 2011). Vegetation in the study area is characterized by the alpine meadow and alpine grassland (Luo et al. 2015).

Fig. 1
figure 1

Source Thaw slumps in 2019 are from Xia et al. (2022). The permafrost map is from the Qinghai-Tibet Plateau Science Data Center for China, National Science and Technology Infrastructure at http://data.tpdc.ac.cn

Topography and permafrost distribution of the study area. The red polygons represent thaw slumps in 2019.

Full size image

The mean annual air temperature and annual precipitation from a nearby meteorological station are − 3.8 °C and about 300 mm, respectively. Most of this region is underlain by ice-rich permafrost: 70% of the ground has an ice content > 30%, and 20% has an ice content > 50%, with permafrost thicknesses ranging from 20 to 80 m (Zhou et al. 2000). The average annual ground temperature from boreholes is − 2.0 to − 0.5°C and the active layer thickness is 1.5–2.0 m (Luo et al. 2015; Wu et al. 2015). Based on the optical remote sensing images, visual interpretation and deep learning methods have been applied to identify RTSs in this area (Huang et al. 2020). The results show that due to abnormally high temperatures and sufficient precipitation during the thawing season, a notable increase in the number of thaw slumps has been observed in the Beiluhe region over the past few decades.

3 Data and Method

In this section, we describe in detail the various data used, including SAR images, optical images, and so on. We also explain the methods used to obtain deformation and to detect RTSs.

3.1 Data Sources

To obtain the deformation characteristics of RTSs, we collected the published datasets of RTS inventory in this study. Xia et al. (2022) used an iteratively semi-automatic method built on deep learning to delineate thaw slumps in the 2019 PlanetScope CubeSat images. A total of 875 RTSs were mapped over a ~54,000 km2 area along the Qinghai Tibet Engineering Corridor. In the same region as this study, 268 slumps were clustered in the Mati Hill region (Fig. 1). This open-source dataset is available at https://doi.org/10.5281/zenodo.6397029 (Xia et al. 2021).

The Sentinel-1A satellite was launched by the European Space Agency (ESA) in April 2014. The descending Sentinel-1A SAR images from October 2017 to July 2021 were downloaded from the ASF (Alaska Satellite Facility) and we used the Level 1 single-look complex (SLC) images in interferometric wide-swath (IW) mode with VV polarization in this study. The original data width is about 250 km and consists of three sub-swaths with overlaps, and swath 1 and 2 were selected as the region of interest for processing (Torres et al. 2012). Precision orbit data (POD) published by the ESA was used to improve the quality of the registration and correct orbit error of the SAR data.

The Shuttle Radar Topography Mission (SRTM) DEM data in 30 m resolution were also adopted to remove the topographic phase, geocode the deformation result, and calculate terrain factors, which include elevation, slope, and aspect. For optical images, we relied on the Sentinel-2A images from 2017 to 2021 and used images in July−August, the peak growing season months of vegetation (Runge et al. 2022). The normalized difference vegetation index (NDVI) was calculated based on the B4 (red) and B8 (near-infrared) bands of the 10 m resolution Sentinel-2 images. The World Imagery with a 0.5 m resolution provided by ESRI can supply us with extra reference for better manual identification of RTSs. Table 1 describes the datasets used in this study.

Table 1 Datasets used in this study
Full size table

3.2 Methods

The flowchart for deformation analysis and latent RTS detection by integrating the InSAR-derived deformation map with optical remote sensing images is shown in Fig. 2. We initially generated the InSAR-derived deformation map. Then, the visual interpretation method with a deformation layer and optical images was employed to delineate the RTSs. Thereafter, spatial masking based on deformation analysis was used to identify latent slumps.

Fig. 2
figure 2

The processing flow chart in this study. RTS retrogressive thaw slumps, SBAS small baseline subsets, SRTM shuttle radar topography mission, DEM digital elevation model, NDVI normalized difference vegetation index

Full size image

3.2.1 Small Baseline Subset InSAR

The Small Baseline Subsets (SBAS)-InSAR method was applied to obtain the surface deformation. SBAS-InSAR is a typical multi-temporal InSAR method proposed by Berardino and Lanari (Berardino et al. 2002). The basic principle of SBAS technology is to generate a differential interferogram based on several SAR image pairs with short space-time baselines. Phase noise is reduced by multi-viewing and filtering the interferogram. We used the spatial complex multilooking to mitigate the phase noise. The range and azimuth look numbers are 5 and 1, respectively, generating the ground pixels of approximately 20 × 20 m. We chose the SAR image acquired on 19 September 2019 as the common reference image, whereas all other images were coregistered to the reference geometry. We constructed interferograms with temporal baseline not exceeding 50 days (Chen et al. 2022). A total of 463 interference pairs were generated (Fig. 3). Then singular value decomposition (SVD) was used to solve multiple small baseline sets to obtain the least square solution under the smallest norm and to obtain the time series deformation sequence in the study area.

Fig. 3
figure 3

Connections between the Sentinel-1 SAR acquisitions with the perpendicular and temporal baselines. The circles denote the Sentinel-1 SAR acquisitions. The black lines denote the interferogram pairs. The red circle refers to the reference image

Full size image

3.2.2 Manually Mapping Thaw Slumps

Retrogressive thaw slumps can be mapped from multi-spectral data as there are color differences between the RTS floor and the surrounding undisturbed vegetated landscape (Witharana et al. 2022). Therefore, visual interpretation of aerial images is a classical approach for the identification and mapping of RTSs, and it is particularly convenient when investigating inaccessible regions (Guzzetti et al. 2012). In previous studies, the deep learning technique has enabled the automatic delineation of RTSs at a relatively high precision. However, a limitation of this method is that it cannot distinguish between two or more RTSs if they are close to each other (Huang et al. 2020). Moreover, some RTSs that have re-vegetated on the surface cannot be identified using remote sensing images alone (Xia et al. 2022).

In this study, we visually interpreted RTSs by integrating deformation results and optical images after obtaining the ground deformation rate map from SBAS processing. First, an annual deformation rate of 1 cm was set as thresholds for the detection and mapping of possible RTSs. We also used the time series information to exclude the deformation unrelated to permafrost thawing. The suspicious RTSs were then superimposed onto the multi-temporal optical remote sensing images and NDVI map and were further screened according to the surface geomorphological features (for example, headwall, recently-exposed sediments, poorly-developed vegetation) (Luo et al. 2019) to eliminate misjudged RTSs and to finalize the RTS mapping. It should be noted that the mapped boundary based on deformation features may result in underestimation because obvious deformation may only occur in some parts of the RTS region during a specific observation period (Kang et al. 2023). Therefore, the optical images were used to depict the outline of the RTSs, guided by the InSAR measurements.

4 Results and Analysis

In this section, we present the results of ground deformation derived from InSAR. Then, we analyze the spatial and temporal deformation characteristics of RTS by combining deformation results and optical images. Finally, we apply the method presented in Sect. 2.2.2 to detect RTSs in the study region and analyze their spatial distribution.

4.1 Surface Deformation Characteristics of Retrogressive Thaw Slumps

According to the deformation results (Fig. 4), most of the permafrost in the study area is in a relatively stable state, with a deformation rate ranging from − 10 to 10 mm/year. The box plots in Fig. 5 highlight a correlation between linear deformation rate and slope angle and vegetation change. The surface deformation in the gentle slope area is faster, that is, decreasing linear velocity with increasing slope angle between 17° and 22°. The rapid deformation area is accompanied by vegetation degradation simultaneously.

Fig. 4
figure 4

Source Xia et al. (2022)

Map of the pluriannual ground surface linear deformation rate in the line-of-sight (LOS) directions for the descending Sentinel-1 track 150 between 2017 and 2021. The black polygons represent the outline of retrogressive thaw slumps (RTSs) in 2019.

Full size image
Fig. 5
figure 5

Box plots showing the distribution of InSAR-derived line-of-sight (LOS) deformation rate under a different terrain slope angles, and b different normalized difference vegetation index (NDVI) values

Full size image

Combining the optical images and distribution of RTSs, we found that the majority of the greater subsidence is located in and around the RTSs, with deformation rates ranging from − 35 to 20 mm/year. In the 268 existing RTSs, the average deformation rate of 145 slumps was more than 5 mm/year, and the extreme deformation rate of 228 collapses was more than 10 mm/year, with only a few collapses in a stable state. Because RTSs occur in the gentle slope area, the variation characteristics of RTSs are more diverse than permafrost thawing. A negative value (red) may be the surface subsidence caused by the melting of ground ice; a positive signal (blue) indicates a displacement close to the LOS direction, which may be the accumulation of sediment at the foot of the slope or the downslope movements toward the satellite. In addition, it can be seen that many deformations happened outside of the RTSs. On the one hand, the RTS inventory used in this study was obtained in 2019 (Xia et al. 2022), while the deformation rate was calculated based on the SAR images from 2017 to 2021. On the other hand, the deformation area does not represent the occurrence of RTSs, and it is more likely to be a prone area of RTSs. As a result, some regions with significant deformation have no RTS.

It can be seen from Fig. 6 that the negative values mainly occur at the headwall of RTSs, while the positive values are mainly distributed in the interior of the scar region and at the foot of the slope. Through the distribution of deformation, the development trend of RTSs can be further assessed. In Fig. 6a, the positive deformation is mainly distributed in the internal area, which shows that the RTS has not been further developed in the past five years, but the meltwater mobilized debris and soil from the headwall, resulting in sediment mobilization and transport (Paquette et al. 2020). In Fig. 6c, the greater slant-range decreasing is mainly located at the headwall of the RTS, which indicates that the thaw slump is still active during the study period, and the mobilized material continues to be transported downslope and exposed ice-rich headwalls (Kokelj et al. 2021). Finally, we found that 42 RTSs remain active during the monitoring period, and it was verified by comparison with multi-temporal optical images. The NDVI inside the RTSs is significantly lower than that in the surrounding areas, which can help us better map the RTS.

Fig. 6
figure 6

Zoom-in view of three typical retrogressive thaw slumps (RTSs). The black polygons represent the outline of RTSs in 2019. a–c World Imagery optical images of typical RTSs in 2021. d–f The average deformation rate of typical RTSs acquired from the Small Baseline Subsets (SBAS)-InSAR. g–i The normalized difference vegetation index (NDVI) of typical RTSs in 2021

Full size image

In combination with InSAR deformation results, optical images, and vegetation information, a typical RTS was specifically analyzed (Fig. 7). The length of this RTS is almost 600 m, and the maximum elevation difference between the front and rear edges reaches 55 m. It can be seen from the optical image and profile that the angular headwall is located at a high altitude, and the NDVI in the undisturbed area is significantly higher than that in the scar area, which is maintained at about 0.1. The debris formed by the ablation of material on the steep headwall has a high water content and therefore flows easily even at a shallow angle. The RTS was formed in 2017, and the exposed headwall continues to develop and expand toward the upslope. The area of this RTS changed from 0.11 to 0.244 km2 from 2017 to 2021, and the largest growth occurred between 2019 and 2020, while the front and both sides of the boundary remain unchanged. We have observed subsidence both inside and near the RTS. The area with rapid deformation (rate > 1 cm/yr) is consistent with the range of retreating headwall.

Fig. 7
figure 7

Analysis of typical retrogressive thaw slump (RTS) variation. a Dynamic change in a thaw slump during the active developing stage. The red arrow represents the direction of RTS expansion. The black dotted line AB is a section line passing through the interior of the RTS. P1–P3 are representative points located at the headwall of the RTS. b The line-of-sight (LOS) deformation rate and normalized difference vegetation index (NDVI) in different years along the profile AB in a

Full size image

In order to further analyze the temporal and spatial variation trend of RTS, we selected three representative points to generate time series results. P1–P3 are located at the headwall of this RTS. The deformation sequence results (Fig. 8) show that all points have a periodic change, and it is composed of three motion stages: stable, abrupt thaw, and linear subsidence. The reason is that the study area is located in the permafrost region. In the “no thaw” case, the permafrost remains unchanged and the active layer undergoes exactly repeating annual freeze-thaw cycles (Liu et al. 2015). In winter, as the soil water in the pores of the active layer freezes into ice, the ground surface uplifts; in summer, the pore ice in the active layer melts into water, and the surface subsides.

Fig. 8
figure 8

The time series change and motion stages of a retrogressive thaw slump (RTS) headwall. See Fig. 7a for the location of P1–P3

Full size image

A typical RTS is caused by the melting of ground ice in sloping terrain, hence the InSAR-derived results include at least two types of ground surface deformation: (1) periodic, seasonal subsidence and uplift; and (2) irreversible subsidence caused by the thawing of ice-rich permafrost. Before the collapse, the surface vegetation coverage is undisturbed, so the region is stable in the early stage. Retrogressive thaw slumps show seasonal cycles and initiate in summer by melting of exposed ground ice in a headwall (Lewkowicz and Way 2019). As a result, the headwall shows an abrupt thaw trend in the deformation sequence. Due to lower air temperature dampening ice melting, the RTS stabilizes in autumn.

After the initiation of thaw slumps, the amplified permafrost disturbances and loss of protective vegetation increase the heat flux to the subsurface and contribute to accelerated permafrost thaw (Patton et al. 2021). In addition, the topographic depression formed by the RTS results in the accumulation of snow in winter. The thicker snowpack subsequently persists into spring. It acts as an insulation that prevents the soil heat to the cold atmosphere (Zhang 2005). In spring, meltwater from the snowpack increases pore water pressure, further decreasing hillslope stability (Patton et al. 2021). Therefore, the result is a long-term linear subsidence trend and increasing seasonal amplitudes (Liu et al. 2015).

4.2 Updated Retrogressive Thaw Slump Inventory and Their Spatial Distribution

Section 3.1 indicated that the RTS region is characterized by rapid deformation and disturbed vegetation. Based on manual interpretation of the optical images, deformation maps, and vegetation, which were taken between 2017 and 2021, a total of 375 RTSs were identified in the permafrost areas of the Mati Hill. Of these RTSs, 268 are consistent with the previous inventory, and the additional 97 slumps are new findings of this study (Fig. 9). It can be seen that most of the RTSs are concentrated in the hilly area in the north and east of the Mati Hill, and a small number of slumps are distributed in the southwest.

Fig. 9
figure 9

All detected retrogressive thaw slumps (RTSs). The size of the circle indicates the area of a RTS

Full size image

Comparing with multi-temporal optical images, we found that 21 RTSs were missing in the previous inventory of the Mati Hill area. Figure 10 shows the detailed deformation and NDVI of two typical RTSs that formed in 2019. In the absence of obvious headwall and slump floor, we cannot confirm it as a RTS simply through color difference and lower NDVI. However, the spatial distribution of deformation is similar to that of a RTS. As a result, RTSs can be mapped with high precision based on multi-source remote sensing data.

Fig. 10
figure 10

Zoom-in view of missing retrogressive thaw slumps (RTSs). a and b Sentinel-2 optical images of typical RTSs in 2019. c and d The deformation rate of typical RTSs. e and f The normalized difference vegetation index (NDVI) of typical RTSs in 2019

Full size image

Figure 11 shows the geometric characteristics of mapped RTSs. The areas of the RTSs range from 0.8 to 20 ha, and 95% of the RTSs are smaller than 8 ha. The minimum perimeter is 117 m, and the maximum perimeter can reach 3.7 km. Of slumps, 84% have a perimeter less than 1000 m. Compared with the landslides in unfrozen areas, the scale of RTS is generally small.

Fig. 11
figure 11

Geometric statistics of the retrogressive thaw slumps (RTSs) based on mapped RTSs. a Histogram showing the perimeter sizes of all RTSs. b Histogram showing the area sizes of all RTSs

Full size image

Figure 12 presents statistics on terrain factors and vegetation characteristics. The RTSs are found at medium elevations (4718–4879 m), with an average elevation of 4,809 m. The slope ranges from 0.93° to 16.5°, with an average of 7.88°. Figures 12a and b respectively show that the RTSs preferentially exist at elevations from 4750 to 4850 m, and on slopes from 4° to 10°. The RTSs are predominantly located on north-facing slopes, which commonly receive less solar radiation and have a thinner active layer than those oriented southward (Fig. 12c). In 2017, the NDVI of RTS areas ranged from 0 to 0.2, so the original surface vegetation coverage was low; by 2021, the NDVI of most RTSs has further reduced. Figure 12d indicates that the development of RTSs will damage the original vegetation and lead to failure to recover in a short term. Through comparison with the terrain of the whole study area (Table 2), it was found that the average elevation and average slope of RTSs are higher than that of the whole study area, but the vegetation index does not show this pattern, which reflects that the spatial distribution characteristics of mapped RTSs are slightly different from the overall regional characteristics. The RTS distribution is affected by elevation, slope, aspect, and vegetation, which is similar to the highlights of previous studies (Yin et al. 2017; Huang et al. 2020; Xia et al. 2022).

Fig. 12
figure 12

Statistics of retrogressive thaw slump (RTS) terrain and vegetation factors. Landscape refers to the entire study area. a The elevation frequencies of the landscape and RTSs. b The slope frequencies of the landscape and RTSs. c The aspect frequencies of RTSs. d The normalized difference vegetation index (NDVI) distribution of RTSs in 2017 and 2021

Full size image
Table 2 Characteristics of the retrogressive thaw slump (RTS) area and the entire study area
Full size table

5 Discussions

In this section, we first compare our results with previous studies. Based on statistical results in Sect. 3, we use the spatial mask method to detect latent RTSs. Finally, the application potential of our method is discussed.

5.1 Comparison with Previous Studies

We used the SBAS-InSAR method to measure the deformation of large-scale RTSs in the Mati Hill area. In the permafrost region of the QTP, many studies have used multi-temporal InSAR method to monitor ground surface deformation. Chen et al. (2013) used ALOS-1 data to estimate surface displacements in Beiluhe from 2007 to 2010. The surface deformation rate was primarily in the range of – 20 mm/year to 20 mm/year. Wang et al. (2020) used TerraSAR-X and Sentinel-1 data to obtain deformation results for the Beilu River Basin from 2018 to 2019 and found that the linear deformation rate ranged from − 24.5 to 5 mm/year. We compared our results with previous studies in the same region (Table 3). In general, the linear deformation rate estimated by this study is higher than that in previous studies that focused on gradual thaw. The disparity may be related to the differences in the observation periods and the density of RTSs. To put this work in a spatial context, we also compared our study with several other studies on thermokarst landforms in other permafrost regions. Chen et al. (2018) used ALOS-1 data to measure the subsidence of thermal erosion gullies in the Eboling Mountain from 2006 to 2011. The subsidence was up to 3 cm/yr inside and around the thermal erosion gullies, which is consistent with the results of this study.

Table 3 Comparison between the results of this study and previous studies
Full size table

In addition, we identified more RTSs by integrating deformation results, optical images, and NDVI. Some of the RTSs were newly developed and others were unrecognized in previous studies. Luo’s inventory was obtained based on images from 2008 to 2017 (Luo et al. 2019), Xia’s results were derived from data in 2019 (Xia et al. 2022), while our results were acquired in 2021. The spatiotemporal variation characteristics derived from InSAR method were captured, which play an important role in detecting latent RTSs.

5.2 Latent Retrogressive Thaw Slump Identification

In the early stage of development, latent RTSs have not yet formed at a large scale, which brings challenges to recognition based on optical images. However, our analysis showed that the latent RTS regions often have a greater deformation rate than the surrounding areas. As a consequence, we mainly detected latent RTSs through deformation masking in this study. It is worth noting that the majority of deformation areas (|velocity| > 10 mm/year) captured by SBAS-InSAR may not be RTS but reflect errors and other landforms such as permafrost thawing or gelifluction. To remove deformation not associated with RTSs, we applied several environmental spatial masks, including elevation, slope, aspect, and NDVI (Huang et al. 2020; Runge et al. 2022; Xia et al. 2022). The specific thresholds are shown in Table 4. After applying multiple masks and removing previous RTSs, a total of 69 latent RTSs were detected in the study area, and a significant increase occurred in the north of the Mati Hill (Fig. 13).

Table 4 Threshold of spatial masks applied
Full size table
Fig. 13
figure 13

Distribution of all retrogressive thaw slumps (RTSs) in the Mati Hill region. The pink circles represent the RTSs in 2019 identified by Xia et al. (2022). The green circles represent newly occurred RTSs, while the blue circles denote RTSs missed by the previous study. The purple circles represent the region with latent RTSs detected by applying the spatial masks

Full size image

Figure 14 shows two typical latent RTSs. Compared with active RTSs (Fig. 6), we found that the latent RTSs have no obvious shape and color characteristics in remote sensing images (Fig. 14a, b). R1 and R2 are two points with greater subsidence in the potential RTS region. Figure 14c shows the time series deformation results for the location of point R1 with a linear subsidence rate of − 13 mm/year and an accumulative deformation of 65 mm. Figure 14d also illustrates the LOS displacement time series of point R2. The LOS displacement rate of R2 was – 11 mm/year. In addition, we found that accelerated subsidence occurred in the summer of 2020, which is an evident feature of ground ice melting.

Fig. 14
figure 14

The deformation characteristics of latent retrogressive thaw slumps (RTSs) in the study area. a World Imagery optical images of R1 in 2021. b World Imagery optical images of R2 in 2018. c InSAR displacement time series of point R1. d InSAR displacement time series of point R2

Full size image

We evaluated the identification results through changes in NVDI, which can act as an indicator of surface disturbance (Runge et al. 2022). We calculated the NDVI values based on the Sentinel-2 images from 2017 to 2022. The results show that there was a significant decrease in NDVI in latent RTS areas after 2021 (Fig. 15), which represented surface disturbance caused by subsequent thaw slumps.

Fig. 15
figure 15

Time series of normalized difference vegetation index (NDVI) from 2017 to 2022 in two latent thaw slump regions

Full size image

5.3 Application Potential of Multi-temporal InSAR in Retrogressive Thaw Slump Investigations

Against the background of climate change, the total RTS-affected area has increased (Luo et al. 2019). Multi-temporal InSAR can be used to map the spatiotemporal variation characteristics of RTSs, which can provide benchmark information for studying the relationship between regional and local forcing conditions. In addition, the InSAR-identified deformation in this study allowed us to obtain a more complete inventory than optical images that are usually limited by cloud covers. The spatial distribution of RTSs is frequently used as a proxy for the past or present occurrence of permafrost, and the kinematics of these thermokarst landforms can be used to derive indirect information about permafrost state.

Furthermore, latent RTSs can be accurately recognized by combining the InSAR-identified deformation, distribution of ground ice, and other geomorphologic elements. On the other hand, the relationship between the melting subsidence of ground ice and the surface deformation of InSAR can be established (Wang et al. 2023), which could provide a reference for predicting the change of RTSs. Therefore, the large-scale InSAR deformation results for the permafrost region of the QTP have important application potential for monitoring the changes of RTSs and other thermokarst landforms.

6 Conclusion

In this study, we used the Sentinel-1A SAR data and the SBAS method to obtain deformation results in the Beiluhe region and acquired the deformation characteristics of RTSs. We applied manual delineation on high-resolution remote sensing images and displacement velocity map derived from InSAR to map RTSs. Further, the latent RTSs were detected based on the distribution characteristics. The following conclusions can be drawn:

  1. (1)

    In the Mati Hill region, the permafrost areas exhibited significant deformation, with linear subsidence rates generally exceeding 10 mm/year. This is mainly due to the widespread distribution of RTS. Compared to normal permafrost degradation, the displacement velocity of RTSs is more rapid (− 35 to 20 mm/year), and the greater deformation is mainly concentrated in the headwall and the central part.

  2. (2)

    For the active RTSs, the greater subsidence area is consistent with the expanded scar zone or the retrogressive headwall. The deformation of headwall is composed of three motion stages: stable, abrupt thaw, and linear subsidence.

  3. (3)

    Compared to the inventory based on optical images in 2019, our inventory overall shows better completeness and novelty. We newly identified 76 RTSs in the study region, which is 26% more than the inventory in 2019. The RTSs preferentially occur in regions with gentle and north-facing slope and tend to appear at medium elevation with poor vegetation coverage.

This study verified the effectiveness of the time series InSAR technique in depicting the spatiotemporal variations of retrogressive thaw slumps, which can provide a vital reference for mapping RTSs on the QTP. However, affected by the resolution of the raw data, deformation signals of some of the smaller-scale RTSs are difficult to extract. Since the dynamic change of RTSs is more complex than conventional permafrost degradation, it is necessary to further improve the deformation analysis of RTSs by combining SAR data from different tracks and satellites. With the successful launch of the LuTan-1 satellite with L-band and the public availability of RADARSAT-1 images, abundant data will be more favorable for carrying out a study about long-time series deformation of RTSs. Combined with deep learning methods, it is expected to obtain highly accurate distribution and deformation properties of RTSs on the Tibetan Plateau.