*Article* **Interferometric SAR Observation of Permafrost Status in the Northern Qinghai-Tibet Plateau by ALOS, ALOS-2 and Sentinel-1 between 2007 and 2021**

**Lichuan Zou 1,2,3, Chao Wang 1,2,3,\* , Yixian Tang 1,2,3, Bo Zhang 1,2,3, Hong Zhang 1,2,3 and Longkai Dong 1,2,3**


**Abstract:** With global warming, permafrost is undergoing degradation, which may cause thawing subsidence, collapse, and emission of greenhouse gases preserved in previously frozen permafrost, change the local hydrology and ecology system, and threaten infrastructure and indigenous communities. The Qinghai-Tibet Plateau (QTP) is the world's largest permafrost region in the middle and low latitudes. Permafrost status monitoring in the QTP is of great significance to global change and local economic development. In this study, we used 66 scenes of ALOS data (2007–2009), 73 scenes of ALOS-2 data (2015–2020) and 284 scenes of Sentinel-1 data (2017–2021) to evaluate the spatial and temporal permafrost deformation over the 83,000 km<sup>2</sup> in the northern QTP, passing through the Tuotuohe, Beiluhe, Wudaoliang and Xidatan regions. We use the SBAS-InSAR method and present a coherence weighted least squares estimator without any hypothetical model to calculate long-term deformation velocity (LTDV) and maximum seasonal deformation (MSD) without any prior knowledge. Analysis of the ALOS results shows that the LTDV ranged from −20 to +20 mm/year during 2007–2009. For the ALOS-2 and Sentinel-1 results, the LTDV ranged from −30 to 30 mm/year during 2015–2021. Further study shows that the expansion areas of permafrost subsidence are concentrated on braided stream plains and thermokarst lakes. In these areas, due to glacial erosion, surface runoff and river alluvium, the contents of water and ground ice are sufficient, which could accelerate permafrost subsidence. In addition, by analyzing LTDV and MSD for the different periods, we found that the L-band ALOS-2 is more sensitive to the thermal collapse of permafrost than the C-band sensor and the detected collapse areas (LTDV < −10 mm/year) are consistent with the GF-1/2 thermal collapse dataset. This research indicates that the InSAR technique could be crucial for monitoring the evolution of permafrost and freeze-thaw disasters.

**Keywords:** permafrost; InSAR; Qinghai-Tibet Plateau; ALOS; ALOS-2; Sentinel-1; thermal melting collapse

## **1. Introduction**

The Qinghai-Tibet Plateau (QTP) is known as the Asian water tower, with an average altitude of more than 4000 m [1]. It is bounded by the Pamir Plateau in the west, Hengduan Mountain in the east, the southern end of the Himalayas in the south and Kunlun Altun Mountain, Qilian Mountain in the north [2]. The QTP is a high terrain and thus receives more solar radiation energy than lower elevation areas [3]. The Chinese mainland climate is affected by the South Asian and East Asian monsoons, resulting in a diversity of climates in different regions, such as the rainy climate in China's southern part of the Yangtze River and drought in Northwest China [4]. In addition, the QTP has many glaciers, lakes, groundwater and surface rivers, making the QTP a super water tower in the plateau area,

**Citation:** Zou, L.; Wang, C.; Tang, Y.; Zhang, B.; Zhang, H.; Dong, L. Interferometric SAR Observation of Permafrost Status in the Northern Qinghai-Tibet Plateau by ALOS, ALOS-2 and Sentinel-1 between 2007 and 2021. *Remote Sens.* **2022**, *14*, 1870. https://doi.org/10.3390/rs14081870

Academic Editors: Takeo Tadono and Masato Ohki

Received: 28 February 2022 Accepted: 9 April 2022 Published: 13 April 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

which affects the water system layout of all of Asia [5]. The QTP is a region with a large amount of permafrost at high latitudes [6]. As a key component of the Earth's cryosphere, permafrost plays an important role in the surface energy balance, carbon and water cycles, terrestrial ecosystem, and hydrological system [7]. In recent years, with global warming, permafrost degradation has accelerated [8], and degradation has had an impact on the environment and the energy and material balance. Therefore, it is very important to monitor the permafrost status on a large scale for a long time series [7].

Traditional measurement methods of permafrost deformation include GPS [9], leveling surveys [10], and drilling [11]. However, due to the harsh environment of the QTP, these methods cannot monitor permafrost on a large scale [8]. The multitemporal interferometric synthetic aperture radar (MT-InSAR) technique is a useful tool to map ground deformation [12]. MT-InSAR has been used to monitor the freeze-thaw cycle of permafrost [13–30], and to retrieve the thickness of the active layer [31–36] and permafrost degradation [37–40]. In these studies, some researchers have been committed to monitoring permafrost for a long time. Zhang [8] used Sentinel-1, ENVISAT and ERS-1 data to evaluate the ground deformation of permafrost and the risk along the Qinghai-Tibet Railway (QTR) from 1997 to 2018. The results show that the estimated deformation rate ranged from −20 to +10 mm/year and most of the QTR appeared to be stable. Daout [41] used ENVISAT and Sentinel-1 data to construct the spatial and temporal dynamics of permafrost deformation in the northeastern QTP from 2003 to 2019. The results show the pervasive subsidence of the permafrost of up to ~2 cm/year, increasing by a factor of 2 to 5 from 2003 to 2019. However, because the C-Band SAR data are easily affected by the region's vegetation and the atmosphere, the results may be affected by spatial and temporal decorrelation. The ALOS Phased Array type L-band Synthetic Aperture Rada (PALSAR) is preferred for ground subsidence monitoring in areas covered by vegetation and where there is a high rate of ground deformation [42]. Therefore, in order to improve the coherence of targets, we used L-band datasets to monitor the ground deformation of permafrost from 2007 to 2021.

The ground deformation process of permafrost is complex. With tectonic activity, erosion, and sedimentation all interacting in the QTP [43], it is difficult to accurately describe the freezing and thawing cycle of permafrost. Therefore, research has attempted to understand the deformation characteristics of permafrost. The sinusoidal model [44,45] and degree-day model [8,46] were used to describe the seasonal variation in the ground surface due to up-down deformation cycles of permafrost. However, it remains controversial which type of model is better at describing seasonal deformation [47]. To extract the temporal characteristics of permafrost directly from the SAR data, Wang [47] directly converted the network of interferograms into a deformation time series without a preset deformation model. Then, the long-term deformation velocity and seasonal deformation were extracted. However, for seasonal deformation, Wang assumed that the highest terrain elevation occurred from January to February, and the lowest elevation occurred from August to October. Wang also averaged the intra-annual deformation value. The average intra-annual deformation may smooth the features of the permafrost deformation. In addition, using prior knowledge may not be suitable for application to the QTP with spatial heterogeneity. In this study, we proposed a long-term deformation velocity and maximum seasonal deformation model without any prior knowledge to directly extract the deformation features of permafrost.

To reveal the status of the permafrost, we extracted time series deformation directly. First, we used 66 scenes of ALOS data (2007–2009), 73 scenes of ALOS-2 data (2015–2020) and 284 scenes of Sentinel-1 data (2017–2021) to reveal the spatial and temporal permafrost deformation in the northern QTP. Second, thermal collapse of permafrost was detected. Finally, we revealed the relationship between the maximum seasonal deformation and the long-term deformation velocity.

#### **2. Study Area and Dataset**

#### *2.1. Study Area*

The study area is located in the Hoh Xil region of the QTP. Hoh Xil is located in Qinghai Province, China, between Kunlun Mountain and Hoh Xil Mountain [48]. The Hoh Xil area, located at a high altitude, covers an area of approximately 2.35×10<sup>5</sup> square kilometers, with an average altitude of 5000 m [49]. The cold climate conditions make Hoh Xil a natural permafrost field. More than 90% of the land is covered by permafrost, with a thickness of 80~120 m [50]. There are a large number of glaciers and lakes. Glaciers have created a large-scale glacial erosion of mountains and produced a large number of deposits [51]. These deposits freeze and thaw repeatedly in the extremely cold environment, and are decomposed into sand of different particle sizes, to form different landscapes [52]. In addition, the ice and snow on the high mountains continue to flow into the Hoh Xil basin and finally converge in the lake [53]. Eventually, the Hoh Xil region evolved many types of landforms. The study area (Figure 1) is located in the northern part of the QTP, passing through the Tuouohe, Beiluhe, Wudaoliang and Xidatan (90.713~95.828◦E, 33.981~36.197◦N). The mean annual ground temperature of the study area is close to 0 ◦C, with elevations from 2600 to 7000 m. The annual precipitation is between 50 mm and 400 mm [54]. The typical geomorphic types include glaciers, lakes, hot melt lakes, alpine meadows, and alpine grasslands. There is a large amount permafrost and seasonally frozen soil which provides abundant opportunities for us to study the permafrost deformation.

**Figure 1.** The study area with topography and surface traces of four overlapping ALOS ascending tracks, three ALOS-2 ascending tracks and two Sentinel-1 descending tracks in the northern QTP.

#### *2.2. ALOS and ALOS-2 Data*

The first SAR dataset was acquired by the Japanese Advanced Land Observing Satellite (L-band Synthetic Aperture Radar ALOS PALSAR). The wavelength of ALOS PALSAR is 23.62 cm with a pixel resolution of 12.5 m and a swath of 70 km. The incidence angles range between 36◦ and 38◦ . The images correspond to the ascending orbit that is northeastlooking [55]. In this study, a total of 66 scenes of ascending stripmap mode ALOS PALSAR data along four adjacent tracks covering seven frames from 2007 to 2009, were collected in the northern QTP. The acquired ALOS data are shown in Table 1. The second SAR dataset was acquired by the Advanced Land Observing Satellite-2 (L-band ALOS-2 PALSAR-2) launched in May 2014 by the Japan Aerospace Exploration Agency (JAXA). The revisit cycle of the ALOS-2 PALSAR-2 is 14 days [56]. In this study, a total of 73 scenes of ascending ALOS-2 PALSAR-2 data along three tracks covering four frames from 2015 to 2020 were collected in the northern QTP. The acquired ALOS-2 data are shown in Table 2.


**Table 1.** ALOS/PALSAR Data.

#### **Table 2.** ALOS-2/PALSAR-2 Data.


## *2.3. Sentinel-1 Data*

Sentinel-1 is a two satellites constellation that carries C-band SAR sensors (~5.6 cm wavelength) under the Copernicus Programme coordinated and managed by the European Space Agency (ESA). Sentinel-1 data are provided to users free of charge. The revisit cycle of Sentinel-1 single satellite A or B is 12 days in all weather conditions, day and night [57]. In this study, a total of 284 scenes of descending Sentinel-1 TOPS data along two tracks covering two frames from 2017 to 2021, were collected in the northern QTP. The Sentinel-1 datasets were acquired in interferometric wide swath (IW) mode. The acquired Sentinel-1 data are shown in Table 3.

**Table 3.** Sentinel-1 Data.


#### **3. Methodology**

A coherence weighted least square estimator without any hypothetical model was used to calculate the long-term deformation velocity (LTDV) and maximum seasonal deformation (MSD) without any prior knowledge. The main steps are shown in Figure 2. The main steps are: (1) importing ALOS, ALOS-2 and Sentinel-1 images and images coregistration; (2) Optimal interference pair selection and interferogram generation; (3) Multilook filtering, phase unwrapping and coherent target selection; (4) Least squares inversion; (5) Atmospheric phase correction and topographic residual correction; (6) Geocoding the average deformation map and extracting the time series deformation; (7) Extracting the long-term deformation velocity and maximum seasonal deformation; and (8) Analyzing the deformation of permafrost with auxiliary data.

**Figure 2.** Flowchart of the processing approach.

#### *3.1. Interferogram and Phase Unwrapping*

First, we converted raw SAR data (ALOS, ALOS-2, Sentinel-1) to SLC format and used DEM data to register the slave image with the master image. Then we selected the optimal interferometric pair by calculating the temporal and spatial baselines of the primary and secondary images. After optimal interference pair selection, the next step is to generate interferograms. In the interferogram, each point has a phase value, which contains much important information. Let us assume point P in the master and slave image. Each SAR image contains phase and amplitude information. For point P, the phase can be modeled as:

$$\phi(P) = \varphi \, + \frac{4\pi}{\lambda} \mathbf{r} + \mathfrak{a} + \mathfrak{n} \tag{1}$$

where ϕ is the original phase controlled by target scatter attributes, such as roughness and moisture, λ is the wavelength of the SAR sensor, r is the distance between the sensor and the target, α is the atmospheric phase, and n represents the noise phase. In the interferogram, for point P, the phase can be modeled as:

$$
\Delta\Phi(P) = \Delta\varphi + \frac{4\pi}{\lambda}\Delta\mathbf{r} + \Delta\mathbf{x} + \Delta\mathbf{v} \tag{2}
$$

where ∆ϕ is the phase difference generated by the change in the target's properties, <sup>4</sup><sup>π</sup> λ ∆r is the distance difference between the two observations to point P, ∆α is the different atmospheric phase, and ∆ν is the noise phase.

Through Equation (2), the interferogram can be generated from two registered SAR images. The range of the interferogram is −π to π. After generating the interferogram, the next step is phase unwrapping. Phase unwrapping is the process of estimating the whole phase difference between a reference pixel and other pixels. The phase unwrapping algorithm calculates the phase difference between adjacent pixels and then integrates these phase differences, which is equivalent to calculating the number of color fringes starting from the reference pixel. For most unwrapping algorithms, minimizing the residual weighting function between the estimated unwrapping phase difference and the original phase difference is used. In this study, we used a minimum cost flow (MCF) unwrapping algorithm. The phase unwrapping algorithm target is to minimize Equation (3):

$$\sum\_{\mathbf{i,j}} \mathbf{w}\_{\mathbf{i,j}}^{(\mathbf{x})} \Big| \Delta \boldsymbol{\phi}\_{\mathbf{i,j}}^{(\mathbf{x})} - \Delta \boldsymbol{\psi}\_{\mathbf{i,j}}^{(\mathbf{x})} \Big| + \sum\_{\mathbf{i,j}} \mathbf{w}\_{\mathbf{i,j}}^{(\mathbf{y})} \Big| \Delta \boldsymbol{\phi}\_{\mathbf{i,j}}^{(\mathbf{y})} - \Delta \boldsymbol{\psi}\_{\mathbf{i,j}}^{(\mathbf{y})} \Big| \tag{3}$$

∆ψ (y) is the unwrapping phase difference in the y direction and ∆ψ (x) is the wrapping phase difference in the x direction.

We used ISCE software (https://github.com/isce-framework/isce2, accessed on 20 March 2022) for this p;rocessing. To increase the signal-to-noise ratio, during processing, we set 14:4, 16:8 and 3:9 multilook in the azimuth and range direction for ALOS, ALOS-2 and Sentinel-1 interferograms, respectively. After the multilook processing, the spatial resolution (azimuth pixel resolution × range pixel resolution) of ALOS, ALOS-2 and Sentinel-1 was 50 m × 10 m, 5 m × 35 m and 23 m × 3 m, respectively. To accelerate the processing rate, we set 15 threads based on an Intel(R) Xeon(R) Gold 6129 CPU with 16 CPU cores for parallel processing.

#### *3.2. Time Series Deformation Processing*

After the phase unwrapping, we selected highly coherent targets (coherence > 0.7) to extract the time series deformation. Let us assume that there are N + 1 scene SAR images coregistered to the master image, and the time order is (*t*<sup>0</sup> . . . *tN*) which can generate M interferograms. The original phase of each SAR image is φ = h 0, φ 1 . . . φ *N* i , where we assume φ <sup>0</sup> = 0. φ contains the ground deformation phase, atmospheric phase, topographic phase and noise phase. For each pixel, the model is described as:

$$
\Delta \phi = \mathbf{A} \phi \tag{4}
$$

where ∆φ represents the interferometric phase and A is the M × N − 1 coefficient matrix. In matrix A, 1 represents the master image, and −1 represents the slave image. Then, the coherence weighted least square method is used to obtain raw phase φ [58,59]:

$$\hat{\boldsymbol{\phi}} = \arg\min \left\| \left| \boldsymbol{W}^{\frac{1}{2}} (\Delta \boldsymbol{\phi} - \boldsymbol{A} \boldsymbol{\phi}) \right| \right\|\_{2} = \left( \boldsymbol{A}^{T} \boldsymbol{W} \boldsymbol{A} \right)^{-1} \boldsymbol{A}^{T} \boldsymbol{W} \boldsymbol{\Delta} \boldsymbol{\phi} \tag{5}$$

The estimated original phase *φ*ˆ includes the following items:

$$
\delta = \delta\_{\rm dis} + \delta\_{\rm tropo} + \delta\_{\rm gecon} + \delta\_{\rm resid} \tag{6}
$$

where *φ*ˆ *dis* represents the ground deformation phase, *φ*ˆ *tropo* represents the atmospheric phase, *φ*ˆ *geom* represents the topographic phase, and *φ*ˆ *resid* represents the residual phase. W is a M × M diagonal weight matrix:

$$\mathcal{W} = \text{diag}\left\{ \frac{2L\gamma^1}{1-\gamma^{1^2}}, \dots, \frac{2L\gamma^{M^2}}{1-\gamma^{M^2}} \right\} \tag{7}$$

$$\gamma^{\hat{\jmath}^2} = \frac{1}{\sigma\_{\Delta \Phi^{\hat{\jmath}}}^2} \tag{8}$$

where *σ* 2 ∆φ*<sup>j</sup>* is the phase variance of the *j*th interferogram calculated through the integration of the phase probability distribution function and L is the number of looks. After the raw phase is calculated by the least squares method, the displacement time series can be obtained:

$$
\delta \phi\_{\rm dis} = \delta - \delta\_{\rm tropo} - \delta\_{\rm gecon} - \delta\_{\rm resid} \tag{9}
$$

For atmospheric phase removal, global atmospheric models (GAMs), particularly ERA5, have great potential in atmospheric phase screen (APS) correction for InSAR applications on the Tibetan plateau [60]. The ERA5 model has been successfully applied in the APS correction in the QTP [41,47,54]. Therefore, in this study, we used ERA-5 reanalysis data to estimate the tropospheric delay [61]. For the topographic residual correction, the topographic phase residual caused by a DEM error was estimated based on the proportionality with the perpendicular baseline time series [62]:

$$\boldsymbol{\phi}^{i} - \boldsymbol{\phi}\_{tropo}^{i} = \frac{-4\pi}{\lambda} \left( \frac{B\_{\perp}^{i}}{r \sin(\theta)} \boldsymbol{z}\_{\varepsilon} + \sum\_{k=0}^{2} \frac{c\_{k}(t\_{i} - t\_{1})^{k}}{k!} + \sum\_{l \in I\_{\mathbb{S}}} s\_{l} H(t\_{i} - t\_{1}) \right) + \boldsymbol{\phi}\_{resid}^{i} \tag{10}$$

where *B i* ⊥ is the perpendicular baseline, *θ* is the incidence angle, *H*(*t<sup>i</sup>* − *t*1) is a Heaviside step function and *z<sup>ε</sup>* , *c<sup>k</sup>* and *s<sup>l</sup>* are the unknown parameters, which can be estimated by minimizing the *L* 2 -norm of the residual phase. The time series deformation processing was implemented by MintPy (https://github.com/insarlab/MintPy, accessed on 20 March 2022) [63].

#### *3.3. Long-Term Deformation Velocity and Maximum Seasonal Deformation*

We can extract two import indices from the time series deformation: long-term deformation trend and maximum seasonal deformation. The long-term deformation may be used to describe the water storage in the active layer. The maximum seasonal deformation describes the active layer thickness. These two indices were obtained by direct calculation without any assumed model. Wang [47] used the intra-annual highest-lowest terrain elevation difference to represent the intra-annual seasonal deformation after extracting and separating the linear trend from the deformation time series. Wang used prior knowledge and assumed that the highest terrain elevation occurred from January–February, and the lowest elevation occurred from August to October. In addition, the averaged value was used for the final intra-annual deformation value. The average intra-annual deformation may smooth the features of the permafrost deformation. In addition, using prior knowledge may not be suitable for application to the QTP with spatial heterogeneity. The long-term linear trend velocity can be modeled as:

$$D'(t) = D(t) - v \cdot t \tag{11}$$

where *D*(*t*) is the time series deformation, *t* is the temporal span of image acquisition dates away from the first acquisition date, *v* is the long-term linear trend velocity, and *D*0 (*t*) is the deformation time series minus the long-term linear trend. The maximum seasonal deformation can be modeled as:

$$\text{msd} = \max\_{i \in T} \left( D'\_{\text{max}} - D'\_{\text{min}} \right) \tag{12}$$

where *D*0 *max* is the highest terrain elevation in year *i* and *D*<sup>0</sup> *min* is the lowest terrain elevation in year *i*. T is the number of years spanned by the observation period.

#### **4. Results and Analysis**

#### *4.1. Long-Term Deformation and Time Series Deformation in the Northern QTP*

The coherence weighted least squares method with 66 scenes of ALOS data (2007– 2009), 73 scenes of ALOS-2 data (2015–2020) and 284 scenes of Sentinel-1 data (2017–2021) was used to retrieve the permafrost deformation in the northern QTP. The results illuminate on a large scale the permafrost deformation in the northern QTP from 2007 to 2021. The ground deformation results are shown in Figure 3.

**Figure 3.** Long-term deformation results of the study area. Background is the DEM. (**a**) ALOS result (2007–2009). (**b**) ALOS-2 result (2015–2020). (**c**) Sentinel-1 result (2017–2021).

In Figure 3, the red color indicates the subsidence in the LOS direction, and the blue color is the uplift toward the satellite. There are a large number of lakes (deep blue color) in the western part of the study area. From the three different deformation maps obtained by using three different SAR sensors, we find that the long-term deformation velocity in the northern QTP ranges from −20 to +20 mm/year (ALOS, 2007–2009), −30 to 15 mm/year (ALOS-2, 2015–2020), and −25 to 30 mm/year (Sentinel-1, 2017–2021). Further study shows that the expansion areas of permafrost subsidence are concentrated on braided stream plains and thermokarst lakes. In these areas, due to glacial erosion, surface runoff and river alluvium, the contents of water and ground ice are sufficient, which could accelerate permafrost subsidence.

In Figure 4a, we selected two typical areas, Salt Lake and a braided stream plain. Salt Lake (35.532◦N, 93.409◦E) is located in the northeastern section of Hoh Xil (Kekexili) in the QTP. The water of Salt Lake mainly comes from seasonal rivers. On 14 September 2011, Zonag Lake burst its banks, which caused the area of Salt Lake to expand significantly. The outburst of Zonag Lake could affect the freeze-thaw behavior of the surrounding permafrost. We used Google Earth images to show the historical expansion of Salt Lake (Figure 4b,d). From the different period images, we find that the Salt Lake experienced significant expansion and that the distance from Salt Lake to the QTR is decreasing. Figure 4c shows the result processed from stacks of ascending ALOS datasets from 2007 to 2009. The red color indicates ground deformation away from the satellite. We find that the degradation of the permafrost mainly occurred in the northeast direction. This may be due to the existence of glaciers and braided stream plains in northeast of Salt Lake. Glaciers have created a large-scale glacial erosion of mountains, resulting in a large number of deposits. In addition, with sufficient moisture, the effect of thawing processes is stronger than that of freezing processes. This could increase the thickness of the active layer. Figure 4d shows the result processed from stacks of descending Sentinel-1 datasets from 2017 to 2021. Similarly, the red color represents the subsidence along the LOS direction. Compared with the 2007– 2009 result and 2017–2021 result, the area of permafrost degradation is expanding. We used ArcMap software to analyze the temporal and spatial changes in the permafrost within a 30-km rectangle around Salt Lake. The calculation result shows an increase of 51.7 km<sup>2</sup> areas (deformation rate > 10 mm/year) within the analysis rectangle after the outburst of Zonag Lake in 2011. The second landscape example is a braided stream plain (34.8791◦N, 93.1845◦E). We used Google Earth image to show the spatial properties of the braided stream plain (Figure 4f). This braided stream plain is surrounded by glaciers in the north, south and east. The terrain in the west is relatively low. There are many lake tributaries in this area. The surrounding glaciers provide many water resources for the braided stream plain. With sufficient moisture, the permafrost may undergo degradation. To illuminate the freeze-thaw cycle of permafrost, we used the ALOS, ALOS-2 and Sentinel-1 datasets to map the ground deformation of permafrost around the braided stream plain from 2007 to 2021. The results are shown in Figure 4g–i.

**Figure 4.** Permafrost evolution in two typical areas, Salt Lake and braided stream plain. (**a**) Average deformation rate during 2017–2021. (**b**) Salt Lake in December 2008 from Google Earth. (**c**) Average deformation rate during 2007–2009. (**d**) Salt Lake in December 2018 from Google Earth. (**e**) Average deformation rate during 2017–2021. (**f**) Braided stream plain from Google Earth image. (**g**) Average deformation rate during 2007–2009. (**h**) Average deformation rate during 2015–2020. (**i**) Average deformation rate during 2017–2021.

To further study the time series deformation of permafrost around Salt Lake, we selected four points (Figure 4c). The details for these points are provided in Table 4. A 25 km profile (Figure 4g) was selected to reveal the relationship between the permafrost deformation and the elevation.


**Table 4.** The details of the four selection points around Salt Lake.

Figure 5a,b shows the long-term deformation trend deformation of selected points around Salt Lake. From March 2007 to April 2008, the four points could balance the freezethaw process. However, from July 2009 to March 2010, point A breaks the cycle, and the ability to thaw is stronger than the ability to freeze, which may increase the thickness of the active layer. From Figure 5b and Table 3, at point C, the average uplift deformation from 2007 to 2008 changed into the average subsidence deformation from 2017 to 2021. At point D, the results reveal that the subsidence of the ground increased 10 times from 2007 to 2021. From Figure 5c, we find that there is a negative correlation between the subsidence and the elevation. This may be due to the melting of glaciers around the braided stream plain, resulting in a large amount of water, accumulating on the slope. Finally, the thickness of the permafrost active layer is increased, which causes permafrost degradation.

**Figure 5.** The time series deformation of the study points and the deformation of the selected profile. (**a**) Time series deformation of study points during 2007–2009. (**b**) Time series deformation of study points during 2017 to 2021. (**c**) The deformation of selected profile.

In the freezing and thawing process of permafrost, the moisture in permafrost can move and change, which may cause rock damage, ground deformation and thawing disasters. The sudden thawing of permafrost may lead to environmental disasters, such as ground collapse and rapid erosion. The thawing mud flow of the permafrost occurs on the slope of the permafrost. In summer, the thawing of the active layer could not flow downslope due to the existence of permafrost, resulting in abundant water in the active layer. Under the condition of excessive moisture, the permafrost creeps downward along the slope, resulting in the thermal melting collapse. We combined optical data (GF-1, GF-2) [64] to reveal the thermal melting collapse of permafrost.

In Figure 6a, the red color indicates the positions of thermal melting collapse, which were generated from 2019–2020 GF-1/2 data. In Figure 6bc, the red color indicates the subsidence in the LOS direction. In Figure 6b, we used the ALOS-2 dataset to reveal the thermal melting collapses of permafrost from 2015 to 2020. The result shows that L-band ALOS-2 can effectively detect the location of thermal melting collapses. The result shows spatial consistency with the GF-1/2 datasets. Compared with the ALOS-2 results, C-band Sentinel-1 rarely detected the location of thermal melting collapses. The results show that the ALOS-2 data can detect displacement that cannot be detected by C-band sensors. To reveal the relationship between the time series deformation and the thermal melting collapses. We selected four points and plotted the time series deformation. The details of the four points are provided in Table 5.

**Figure 6.** The thermal melting collapse, average deformation and time series deformation of the study points. (**a**) Thermal melting collapse from GF-1/2 during 2019–2020. (**b**) Average deformation rate during 2015–2020. (**c**) Average deformation rate during 2017–2021. (**d**) Time series deformation of points during 2015–2020. (**e**) Time series deformation of points during 2017–2021.



From Table 5, from 2015 to 2020, the L-band ALOS-2 result is different from the Cband Sentinel-1 result. For the ALOS-2 result, the permafrost under points A, B and C is undergoing subsidence. The Sentinel-1 result shows that the permafrost is still stable. Due

to the different observation periods of the two sensors, ALOS-2 has a longer observation period. In addition, results are also affected by the internal characteristics of the sensor, such as wavelength and incident angle.

#### *4.2. Long-Term Deformation Velocity and Maximum Seasonal Deformation in the Northern QTP*

After determining the time series deformation derived from the raw phase, the longterm deformation trend and maximum seasonal deformation magnitude can be calculated. The spatial distributions of the maximum seasonal deformation magnitude are shown in Figure 7.

**Figure 7.** The long-term deformation velocity and maximum seasonal deformation. (**a**) Longterm deformation velocity during 2015–2020 from ALOS-2. (**b**) Long-term deformation velocity during 2017–2021 from Sentinel-1. (**c**) ALOS-2 maximum seasonal deformation. (**d**) Sentinel-1 maximum seasonal deformation. (**e**) Maximum seasonal deformation during 2015–2020 from ALOS-2. (**f**) Maximum seasonal deformation during 2017–2021 from Sentinel-1. (**g**) ALOS-2 LTD and MSD. (**h**) Sentinel-1 LTD and MSD.

Due to the limitation of the revisit period of ALOS data, we calculated the long-term deformation velocity (Figure 7a,b) and the maximum seasonal deformation (Figure 7e,f) of ALOS-2 and Sentinel-1 data. The experimental results show that the average maximum seasonal deformation of the two sensors is 66 mm. As shown in Figure 7a,b,e,f, the area with large maximum seasonal deformation may tend to undergo strong subsidence and uplift. From the ALOS-2 result, the average maximum seasonal deformation during 2015–2020 is 68.74 mm. According to the Sentinel-1 result, the average MSD during 2017–2021 is 64.85 mm. The areas with high MSD are mainly concentrated in glacier outwash plains and braided stream plains. We selected two typical areas, A and B. The upper left of study area A is located in Duoergai Co Lake, surrounded by glaciers to the north. The MSD is large, which may be due to glacial erosion. A large amount of ice deposits is brought to low-lying areas, affecting the water storage capacity and ice melting capacity of permafrost. In addition, a large MSD occurred in thermal karst landforms, such as area B, which is characterized by the formation of hot melt lakes and proglacial lakes. Due to the thermal influence of underground ice in the permafrost, the land shrinks and settles, and finally melted water flows along the groove of the slope. In addition, we reveal the relationship between the maximum seasonal deformation and the long-term deformation velocity. From Figure 7g,h, there is a weak correlation between MSD and LTD, with coefficients of

determination (*R* 2 ) [65] of 0.208 for ALOS-2 and 0.301 for Sentinel-1. This shows that the relation between MSD and LTD is not simply linear and is related to geomorphic types, moisture, ice and other factors.

#### **5. Discussion**

To evaluate the permafrost deformation in the northern QTP, we compared other studies that are shown in Table 6.


**Table 6.** Previous studies on permafrost deformation in the QTP.

Xie et al. [66] used the PSI technique with Envisat datasets to reveal the ground deformation in Beiluhe. The result shows that the average deformation from 2003 to 2007 was −20~3 mm/year. Chen [67] used multisource datasets (ALOS and Envisat) based on the IPTA and SABS techniques to illuminate the displacement of permafrost in Beiluhe. The result shows that the average displacement is −20~20 mm/year. Chen et al. [28] found that the average deformation of permafrost in Beiluhe is −25~10 mm/year by using the SBAS method with ALOS data. From 2014 to 2020, Li et al. and Wang et al. [68] used the MT-InSAR method to research the deformation of permafrost in the QTP. Li et al. [69] found that the average deformation in Yangbaijin ranged from −30 to 10 mm/year from 2014 to 2015. Wang et al. [54] used Sentinel-1 data based on NSBAS to reveal the deformation of the whole QTP. The result shows that the average deformation rate is −20~20 mm/year. Zhang et al. [8] used three different sensor datasets (ERS-1, Envisat and Sentinel-1) to reveal the freeze-thaw cycle of permafrost in Wudaoliang and Tuouohe from 1997 to 2008. The result shows that the average deformation rate ranges from −20 to 10 mm/year. From the above related research, our results are consistent with the previous results.

In addition, to evaluate the accuracy of the SBAS-InSAR results in this paper, we used leveling data from other researchers in the QTR [70]. The locations of the three leveling points are shown in Figure 8. The elevations of the three leveling points were collected on 9 January 2018, and 11 February 2018. The deformation results of 3 January 2018 and 8 February 2018, were extracted from 150-475 (Path-Frame) Sentinel-1 InSAR results. In order to ensure the accuracy, we convert the leveling results to the displacement of the LOS direction based on the radar incident angle. The comparison between the InSAR experiment results and the leveling results of A, B and C is shown in Table 7. Table 7 shows that the absolute errors of the leveling data and InSAR measurements at points A, B and C were 3.1, 5.7 and 6.8 mm, respectively, which were relatively small.

Through a comparison with previous studies and leveling data, this paper provides a theoretical basis and data support for the study of permafrost evolution. However, there are some limitations in our study. First, there are some errors in the InSAR processing, such as orbit errors, unwrapping phase errors, residual topographic errors, and ERA5 APS correction errors. These errors may lead to incorrect permafrost deformation. Second, from 2010 to 2015, due to the lack of data, we have no ability to monitor the deformation of permafrost in the QTP during this period. In the future, we will combine more SAR data

and auxiliary data to help us better understand the freeze-thaw cycle and degradation of permafrost in the QTP.

**Figure 8.** Geographic location map of leveling points along the QTR.

**Table 7.** The absolute error between the measured values of leveling points A, B, and C and the InSAR observed values.


#### **6. Conclusions**

In this study, we used 66 scenes of ALOS data (2007–2009), 73 scenes of ALOS-2 data (2015–2020) and 284 scenes of Sentinel-1 data (2017–2021) based on SBAS-InSAR to reveal the spatial and temporal deformation of permafrost in the northern QTP, passing through the Tuotuohe, Beiluhe, Wudaoliang and Xidatan regions between 2007 and 2021. In addition, a coherence weighted least squares estimator without any hypothetical model was used to calculate the long-term deformation velocity (LTDV) and maximum seasonal deformation (MSD) without any prior knowledge. The conclusions are summarized as follows:

1. Analysis of the ALOS results shows that the LTDV ranged from −20 to +20 mm/year during 2007–2009. For the ALOS-2 and Sentinel-1 results, the LTDV ranged from −30 to 30 mm/year during 2015–2021. Over a 15-year observation period, permafrost degradation is accelerating. The expansion areas of permafrost degradation are concentrated on braided stream plains, thermokarst lakes and glacier outwash plains. In addition, the high MSDs are mainly concentrated in glacier outwash plains and

braided stream plains, which may be due to glacial erosion and a large amount of ice deposits brought to low-lying areas, affecting the water storage capacity and ice melting capacity of permafrost.


In the future, more SAR datasets and auxiliary datasets will be used to illuminate the freeze-thaw cycle and retrieve the active layer thickness on the QTP on a large scale.

**Author Contributions:** Conceptualization, L.Z., C.W. and Y.T.; formal analysis, L.Z.; methodology, L.Z. and C.W.; software, L.Z., H.Z., B.Z. and L.D.; validation, L.Z., C.W. and Y.T.; writing—original draft, L.Z.; writing—review and editing, C.W.; Y.T. and H.Z.; visualization, L.Z.; supervision, C.W.; project administration, H.Z.; funding acquisition, C.W. All authors have read and agreed to the published version of the manuscript.

**Funding:** This paper was funded by the National Natural Science Foundation of China, Grant No. 41930110.

**Data Availability Statement:** ALOS is provide by JAXA and available from the Alaska Satellite Facility (ASF) (https://vertex.daac.asf.alaska.edu, accessed on 14 February 2022). ALOS-2 data is provided and available from JAXA EORA2 program. Sentinel-1 data is provided by the European Space Agency (ESA) and available from the Alaska Satellite Facility (ASF) (https://vertex.daac.asf. alaska.edu, accessed on 14 February 2022). GF-1/2 is provided by the National Tibetan Plateau Data Center and available from National Tibetan Plateau Data Center (http://data.tpdc.ac.cn, accessed on 14 February 2022). SRTM DEM data is available at https://srtm.csi.cgiar.org/srtmdata/, accessed on 14 February 2022.

**Acknowledgments:** The authors would like to thank JAXA EORA2 program (project id: ER2A2N030) and its ALOS, ALOS-2 data, NASA, CGIAR-CSI for SRTM DEM data, ESA for Sentinel-1 data and National Tibetan Plateau Data Center for GF-1/2 data. Profs. Qingbai Wu, Ji Chen, Jing Luo and Beiluhe Permafrost Station of the Northwest Institute of Eco-Environment and Resources of CAS, and Drs. Fan Wu and Lu Xu of our team are acknowledged for helpful discussion and field support.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

