Characteristics of Freeze–Thaw Cycles in an Endorheic Basin on the Qinghai-Tibet Plateau Based on SBAS-InSAR Technology

: The freeze–thaw (F-T) cycle of the active layer (AL) causes the “frost heave and thaw settlement” deformation of the terrain surface. Accurately identifying its amplitude and time characteristics is important for climate, hydrology, and ecology research in permafrost regions. We used Sentinel-1 SAR data and small baseline subset-interferometric synthetic aperture radar (SBAS-InSAR) technology to obtain the characteristics of F-T cycles in the Zonag Lake-Yanhu Lake permafrost-affected endorheic basin on the Qinghai-Tibet Plateau from 2017 to 2019. The results show that the seasonal deformation amplitude (SDA) in the study area mainly ranges from 0 to 60 mm, with an average value of 19 mm. The date of maximum frost heave (MFH) occurred between November 27th and March 21st of the following year, averaged in date of the year (DOY) 37. The maximum thaw settlement (MTS) occurred between July 25th and September 21st, averaged in DOY 225. The thawing duration is the thawing process lasting about 193 days. The spatial distribution differences in SDA, the date of MFH, and the date of MTS are relatively signiﬁcant, but there is no apparent spatial difference in thawing duration. Although the SDA in the study area is mainly affected by the thermal state of permafrost, it still has the most apparent relationship with vegetation cover, the soil water content in AL, and active layer thickness. SDA has an apparent negative and positive correlation with the date of MFH and the date of MTS. In addition, due to the inﬂuence of soil texture and seasonal rivers, the seasonal deformation characteristics of the alluvial-diluvial area are different from those of the surrounding areas. This study provides a method for analyzing the F-T cycle of the AL using multi-temporal InSAR technology.


Introduction
The upper layer of permafrost (active layer) experiences a freeze-thaw (F-T) cycle seasonally. The F-T cycle affects the climate, ecology, and water cycle process in the permafrost region [1][2][3][4] by influencing the soil thermal conductivity, surface radiation balance, soil moisture evaporation, and infiltration [5][6][7]. Caused by the difference in physical properties between ice and water, a seasonal deformation of "frost heave and thaw subsidence" was accompanied by the F-T cycle of the active layer (AL). The seasonal deformation amplitude (SDA) can indicate the soil water content of the AL [8] and directly affects the engineering construction stability in the cold region [9][10][11][12]. The timing of the F-T cycle also determines the growing season of vegetation [13,14], affects the energy exchange between land and atmosphere [7,15], and changes the carbon release in permafrost [16,17].
The Qinghai-Tibet Plateau (QTP) is the largest frozen soil region in middle and low latitudes [18], and the permafrost area accounts for approximately 41% of the QTP [19]. Permafrost degradation is characterized by ground temperature increase [20][21][22][23], AL thickening [24][25][26], and permafrost area reduction [27,28] during the past 30 years. Meanwhile, the timing of the F-T cycle has also exhibited changes due to climate warming and permafrost degradation, which is manifested as the delayed onset of freezing and the advanced onset of thawing, and the extended thawing duration on the QTP [5,[29][30][31][32]. Accurate identification of the timing, duration, and amplitude of seasonal deformation is essential for revealing the hydrothermal state of the AL and is important for the studies of hydrology and ecology in the permafrost regions.
The monitoring stations on QTP are sparse due to the extremely harsh climate and complex terrain. An accurate F-T state is one of the primary goals of remote sensing applications on the QTP [33,34]. Microwave radiometers have proven the ability to detect near-surface (about 0-20 cm) ground F-T state [35][36][37][38][39][40] and been applied successfully on the QTP [41][42][43][44][45]. However, due to its coarse spatial resolution (~25 km), its application in fine-scale research is limited. Active microwave sensors, including synthetic aperture radar (SAR), mainly obtain the F-T state through the backscattering coefficient [46][47][48][49][50], generally have a fine spatial resolution (1 km or~10 m for SAR) but are highly susceptible to scatterings from vegetation, snow, and surface roughness changes other than ground F-T state [51][52][53]. The combination of active and passive microwave technology has become a reliable method for monitoring the F-T state [38,[53][54][55][56][57][58][59]. The studies mentioned above focused on identifying the F-T state of the near-surface. Over permafrost terrain, in addition to the F-T state of the near-surface, the F-T state and F-T processes in the AL [44,53] are also very important but still lack investigation.
In this study, we took the Zonag Lake (ZL)-Yanhu Lake (YL) basin as the study area and conducted a detailed analysis of seasonal deformation, and firstly the time characteristics of F-T in the AL were involved in the analysis. The ZL-YL basin was an endorheic region in the hinterland of permafrost on the QTP. However, due to precipitation increase and permafrost degradation, the lake expanded continuously from 2011 to 2019 [91], making the basin a new component of the Yangtze River source. Monitoring the permafrost deformation in the ZL-YL basin is significant for assessing permafrost status and relevance to the hydrological cycle research of the Yangtze River. Previous studies have found the large deformation of both long-term subsidence and seasonal amplitude due to permafrost degradation in the ZL-YL basin [71,75,77,85,92], but the detailed information on the time characteristic and SDA of the F-T cycle is relatively limited.
This study used SBAS-InSAR technology to analyze AL's SDA, freeze/thaw transition, and duration time in the middle latitude QTP permafrost region. The relationships between seasonal deformation characteristics and their influencing factors are also discussed. We first monitored ground surface deformation by SBAS-InSAR technology and then extracted seasonal deformation (including amplitude, the timing of F-T transition, and duration) from the deformation time series. After that, we obtained SDA characteristics and F-T time of AL in the ZL-YL basin. Finally, the seasonal deformation properties were analyzed using multiple influence factor data. The AL's F-T characteristics were first revealed in the middle latitude permafrost environment, and the proposed investigation method can be transferred to other permafrost regions.

Study Area
The ZL-YL basin is located in the Hoh Xil inflow area in the permafrost hinterland of QTP [19]. The average elevation and area of the basin are approximately 4905 m (a.s.l.) and 8.52 × 10 3 km 2 , respectively. The basin terrain is high in the west and low in the east, and small glaciers exist above 5000 m (a.s.l.). The four lakes in the ZL-YL basin are Zonag Lake, Kusai Lake, Heidin Neur Lake, and Yanhu Lake from west to east. The elevation of the four lakes decreases from west to east, among which the height difference between Zonag Lake and Kusai Lake is the largest, about 278 m; the height difference between the three lakes in the east is 9~18 m. The primary recharge sources include the thawing water of ice and snow on the south slope of the Kunlun Mountains and precipitation [93]. According to the long-term monitoring at Wudaoliang station from 1961 to 2017, the average annual precipitation is 294.2 mm, and the average annual air temperature is −4.17 • C, which belongs to a semi-arid climate.
We conducted a field survey in the study area by drilling, pit exploration, and groundpenetrating radar in 2012 and established a comprehensive meteorological monitoring station near the ZL05 borehole on the south bank of the Zonag lake. The temperature observation from 2013 to 2020 shows that the thawing duration is from May to September, and the freezing duration is from October to April. Except for some seasonally frozen ground around the lake due to the thermal influence of water, the types of frozen soils in most areas are permafrost. In the study area, vegetation mainly develops in intermountain basins with relatively low altitudes, such as areas around lakes, while in mountain areas with relatively high altitudes, the surface cover is mainly bare bedrock [94].
After the outburst of Zonag lake in September 2011, its water flowed through Kusai Lake and Heidin Neur Lake and finally poured into Yanhu lake. As a result, the four lakes are connected through surface waterways, forming a typical inflow area. From 2011 to 2019, the area of Yanhu lake continued to increase, and the total catchment area of the study area reached 8.48 × 10 3 km 2 [91]. After a continuous expansion, Yanhu lake is currently only 8 km away from the Qinghai-Tibet Railway. A connecting channel dug on the south bank of the Yanhu lake in 2019 made the basin flow into the Chumar River in the south of the study area, becoming a source of the Yangtze River.

Field Observation Data
The soil information data (soil texture, mass water content, and soil bulk density), types and coverage of vegetation, permafrost thickness, and ALT were obtained during the field investigation. In addition to the comprehensive observation site located in ZL05, we also selected the data of four boreholes for analysis. The four boreholes, HL01, HL02, HL03, and HL04, are all located in river valleys with flat terrain (slope ≤ 1 • ), with vegetation coverage above 85%. Except for the barren land, water bodies, and glaciers, the vegetation-covered ground is classified into four groups: alpine swamp meadow (ASM), alpine meadow (AM), and alpine steppe (AS), and alpine desert (AD) [94]. The monitoring data used in this study include vegetation form, ALT, volume water content (calculated by mass water content and soil bulk density) in AL, and soil texture data (Table 1). In addition, the leveling data observed at Wudaoliang ( Figure 1) were used to verify the SBAS-InSAR observations, and the time range of leveling data used for verification is from 2017 to 2019. The leveling data used for each verification is the average of 25 points in the observation field, covering an area of about 400 m 2 . For detailed information about the devices and methods of the deformation observation field, please refer to other studies from our team [95,96].
Remote Sens. 2022, 14, x FOR PEER REVIEW 4 of 23 vegetation coverage above 85%. Except for the barren land, water bodies, and glaciers, the vegetation-covered ground is classified into four groups: alpine swamp meadow (ASM), alpine meadow (AM), and alpine steppe (AS), and alpine desert (AD) [94]. The monitoring data used in this study include vegetation form, ALT, volume water content (calculated by mass water content and soil bulk density) in AL, and soil texture data (Table 1). In addition, the leveling data observed at Wudaoliang ( Figure 1) were used to verify the SBAS-InSAR observations, and the time range of leveling data used for verification is from 2017 to 2019. The leveling data used for each verification is the average of 25 points in the observation field, covering an area of about 400 m 2 . For detailed information about the devices and methods of the deformation observation field, please refer to other studies from our team [95,96].

Remote Sensing Data
C-band Sentinel-1 SAR data were used to monitor the surface deformation and downloaded from Alaska Satellite Facility (https://vertex.daac.asf.alaska.edu, accessed on 1 April 2020). Level 1 single look complex (SLC) images of interferometric wide-swath (IW) mode with VV polarization in the descending orbit were used in the study, and the orbit number is 150. Finally, we chose 112 acquisition dates from October 2014 to December 2019.
The daily temperature determines the freeze and thaw onsets of the land surface. It is obtained by fitting 8-days MODIS LSTs based on an empirical multi-linear regression model [97], and it has been well used in the study of the distribution permafrost on the QTP [19]. For a better interpretation of the seasonal deformation, the potential influencing factors, including underlying surface data and permafrost-related data, were selected for analysis ( Table 2). The DEM data was used to calculate the slope and aspect. The vegetation form map was produced from field investigation [94]. The NDVI is the growing season average from 2017 to 2019. Soil organic carbon data within 3 m depth was simulated based on field investigation [98]. The soil texture data is the clay, sand, and silt percentage within 3 m depth [99]. The mean annual air temperature and mean annual precipitation are reanalysis data [100]. The permafrost-related data mainly include ALT [101], mean annual ground temperature (MAGT) [101], and ground ice content data [18]. The maximum degree days of freezing (DDF) and maximum degree days of thawing (DDT) were calculated based on the daily average temperature data calculated by MODIS LST. Figure 2 shows the flowchart of SBAS-InSAR processing in this study. We used InSAR Scientific Computing Environment (ISCE) software to generate interferograms [102]. Except for the last image in the time series, all master images interfered with two adjacent images (slave images) in the subsequent time series [103]. The multi-look parameters are 3 and 13 in distance direction and azimuth direction, respectively, and the final pixel is a rectangle of about 40 m. The precise orbit files for correcting track errors are obtained from the Payload Data Ground Segment and downloaded from (https://qc.sentinel1.eo.esa.int, accessed on 1 April 2020). We used an adaptive filtering method to filter the interferograms [104], the filtering strength is 0.6, and the window is 5 × 5. The 1-arcsecond SRTM DEM was used to remove terrain errors, and then the statistical-cost network-flow algorithm for phase unwrapping (SNAPHU) [105] was used to unwrap the interferograms. Finally, we obtained 324 unwrapped interferograms through manual inspection. accessed on 1 April 2020). We used an adaptive filtering method to filter grams [104], the filtering strength is 0.6, and the window is 5 × 5. The 1-arcs DEM was used to remove terrain errors, and then the statistical-cost netwo rithm for phase unwrapping (SNAPHU) [105] was used to unwrap the int Finally, we obtained 324 unwrapped interferograms through manual inspec After acquiring unwrapped phases, we used the Miami InSAR Time-se (MintPy) [106] for time series analysis, and the coherence threshold of time sion is 0.7. The Generic Atmospheric Correction Online Service for InSAR ( lizes the Iterative Tropospheric Decomposition model [107] to separate strat bulent signals from total tropospheric delays, and we used it to correct the error of SBAS-InSAR measurements. Because the lake is not permafrost, it is sible to study the deformation characteristics of the lake, so we masked all lak mask data came from the latest lake data [108] and was used to remove the the water body. Finally, we obtained the time series deformation of the sa sight direction.

Seasonal Deformation Characteristics Extraction
In order to separate the seasonal deformation from the long-term deform a permafrost terrain deformation model that describes linear long-term trends and intra-annual deformation dynamics was used in this study [64,69 where di is the deformation of the ith day in every year relative to the refere time of the first SAR image, that is, 14 January 2017). The v is the long-term rate, A is deformation amplitude, T is the period of seasonal deformation c After acquiring unwrapped phases, we used the Miami InSAR Time-series software (MintPy) [106] for time series analysis, and the coherence threshold of time series inversion is 0.7. The Generic Atmospheric Correction Online Service for InSAR (GACOS) utilizes the Iterative Tropospheric Decomposition model [107] to separate stratified and turbulent signals from total tropospheric delays, and we used it to correct the atmospheric error of SBAS-InSAR measurements. Because the lake is not permafrost, it is not very sensible to study the deformation characteristics of the lake, so we masked all lakes. The water mask data came from the latest lake data [108] and was used to remove the influence of the water body. Finally, we obtained the time series deformation of the satellite line of sight direction.

Seasonal Deformation Characteristics Extraction
In order to separate the seasonal deformation from the long-term deformation trend, a permafrost terrain deformation model that describes linear long-term deformation trends and intra-annual deformation dynamics was used in this study [64,69].
where d i is the deformation of the ith day in every year relative to the reference time (the time of the first SAR image, that is, 14 January 2017). The v is the long-term displacement rate, A is deformation amplitude, T is the period of seasonal deformation caused by the F-T cycle of AL (here set to 365 days), ϕ is the timing of the seasonal deformation, and ε is the noise and deformation residuals caused by other processes. SDA is the maximum deformation magnitude caused by the ice-water phase transition in AL, which can be expressed as the peak-to-peak value (2|A|). By fitting the SBAS-InSAR-derived deformation time series, the unknown parameters in Equation (1) can be acquired. The maximum frost heave (MFH) and maximum thaw settlement (MTS) are the seasonal deformation of the AL, and the F-T time characteristics of the AL based on the deformation in this study mainly come from these two indexes ( Figure 3).
Remote Sens. 2022, 14, x FOR PEER REVIEW 7 of 23 F-T cycle of AL (here set to 365 days), φ is the timing of the seasonal deformation, and ε is the noise and deformation residuals caused by other processes. SDA is the maximum deformation magnitude caused by the ice-water phase transition in AL, which can be expressed as the peak-to-peak value (2|A|). By fitting the SBAS-InSAR-derived deformation time series, the unknown parameters in Equation (1) can be acquired. The maximum frost heave (MFH) and maximum thaw settlement (MTS) are the seasonal deformation of the AL, and the F-T time characteristics of the AL based on the deformation in this study mainly come from these two indexes ( Figure 3). Since the AL's thawing duration and freezing duration form an F-T cycle (365 days), the time characteristics calculated based on the surface deformation include the date of MFH, the date of MTS, and the thawing duration. The calculation method of the thawing duration is as follows: The tMTS and tMFH are the date of MTS and the date of MFH, respectively. The Δt is the thawing duration, representing the duration from the date of MFH to the date of MTS.

Influencing Factor Analysis
Previous studies have shown that seasonal deformation has apparent spatial heterogeneity [8,71,78]. This study used geographical detection to analyze the influencing factor of SDA and the correlation analysis method to analyze the relationship between SDA and time characteristics. The core idea behind the geographical detector is that the similarity of spatial patterns between ground deformation and relevant environmental factors is determined by the q-statistic, determined by the local and global variances [109]. The value of the q-statistic is between 0 and 1, where 1 means the factor controls the spatial patterns of ground deformation and vice versa. Due to AL's lack of water content data on a spatial scale, we only analyzed the influence of volume water content of AL on SDA at borehole sites. Since the AL's thawing duration and freezing duration form an F-T cycle (365 days), the time characteristics calculated based on the surface deformation include the date of MFH, the date of MTS, and the thawing duration. The calculation method of the thawing duration is as follows: The t MTS and t MFH are the date of MTS and the date of MFH, respectively. The ∆t is the thawing duration, representing the duration from the date of MFH to the date of MTS.

Influencing Factor Analysis
Previous studies have shown that seasonal deformation has apparent spatial heterogeneity [8,71,78]. This study used geographical detection to analyze the influencing factor of SDA and the correlation analysis method to analyze the relationship between SDA and time characteristics. The core idea behind the geographical detector is that the similarity of spatial patterns between ground deformation and relevant environmental factors is determined by the q-statistic, determined by the local and global variances [109]. The value of the q-statistic is between 0 and 1, where 1 means the factor controls the spatial patterns of ground deformation and vice versa. Due to AL's lack of water content data on a spatial scale, we only analyzed the influence of volume water content of AL on SDA at borehole sites.

Validation with Leveling Data
In the Wudaoliang leveling observation site, the SBAS-InSAR-derived deformation time series was extracted and converted from the line of sight direction to the vertical ground direction compared with leveling observation. Those observations with SBAS-InSAR and leveling near time, smaller than 10 days, were used for comparison (Table 3).

Validation with Leveling Data
In the Wudaoliang leveling observation site, the SBAS-InSAR-derived deformation time series was extracted and converted from the line of sight direction to the vertical ground direction compared with leveling observation. Those observations with SBAS-In-SAR and leveling near time, smaller than 10 days, were used for comparison (Table 3).    Figure 5 shows the SBAS-InSAR-derived deformation time series and the corresponding model-fitted lines at four boreholes (HL01, HL02, HL03, and HL04 marked in Figure 1). There is a good correlation between the model fitting and the SBAS-InSAR observation, with R 2 above 0.7 at four sites. The SDAs of HL01, HL02, HL03, and HL04 are 42, 44, 41, and 26 mm, respectively. Borehole HL04 has smaller seasonal deformation than the other three sites.  Figure 5 shows the SBAS-InSAR-derived deformation time series and the corresponding model-fitted lines at four boreholes (HL01, HL02, HL03, and HL04 marked in Figure 1). There is a good correlation between the model fitting and the SBAS-InSAR observation, with R 2 above 0.7 at four sites. The SDAs of HL01, HL02, HL03, and HL04 are 42, 44, 41, and 26 mm, respectively. Borehole HL04 has smaller seasonal deformation than the other three sites. According to the field investigation, the four boreholes are in river valleys with flat terrain, vegetation coverage above 85%, and ALT around 1.8 m, and the soil water content in AL of HL01, HL02, and HL03 are similar and significantly higher than HL04 (Table 1). We can find that when the topography, vegetation, and the ALT are similar, the SDA has an obvious positive correlation with the water content in AL. Figure 6 shows that the SDA is about 0~60 mm, with an average value of 19 mm in the study area. According to the statistical results, two apparent SDA classifications exist in the ZL-YL basin. In mountainous areas, the SDA is 0~20 mm with concentrated distribution, and its central value is about 12 mm, while in intermountain basins, the SDA is 20~60 mm with wide distribution, and its central value is about 32 mm. In the exposed area of bedrock at relatively high altitudes, the range of SDA is 0~10 mm. There is a significant difference in the SDA between the upstream and downstream of the basin. The SDA around the basin of Zonag lake is mainly 20~50 mm, while around the basin of Kusai Lake, Heidin Neur Lake, and Yanhu lake downstream is mainly 20~60 mm, among which the maximum SDA around the basin of the Yanhu lake can exceed 60 mm. According to the field investigation, the four boreholes are in river valleys with flat terrain, vegetation coverage above 85%, and ALT around 1.8 m, and the soil water content in AL of HL01, HL02, and HL03 are similar and significantly higher than HL04 (Table 1). We can find that when the topography, vegetation, and the ALT are similar, the SDA has an obvious positive correlation with the water content in AL. Figure 6 shows that the SDA is about 0~60 mm, with an average value of 19 mm in the study area. According to the statistical results, two apparent SDA classifications exist in the ZL-YL basin. In mountainous areas, the SDA is 0~20 mm with concentrated distribution, and its central value is about 12 mm, while in intermountain basins, the SDA is 20~60 mm with wide distribution, and its central value is about 32 mm. In the exposed area of bedrock at relatively high altitudes, the range of SDA is 0~10 mm. There is a significant difference in the SDA between the upstream and downstream of the basin. The SDA around the basin of Zonag lake is mainly 20~50 mm, while around the basin of Kusai Lake, Heidin Neur Lake, and Yanhu lake downstream is mainly 20~60 mm, among which the maximum SDA around the basin of the Yanhu lake can exceed 60 mm.

Spatial Distribution of Seasonal Deformation
ote Sens. 2022, 14, x FOR PEER REVIEW 10 of Figure 6. The spatial distribution (a) and statistical characteristics (b) of seasonal deformation ( time-series deformation of the boreholes is shown in Figures 3 and 5). Figure 3 shows that, in the ZL05, the date of MFH lags behind the lowest surfa temperature by 35 days, while the date of MTS lags behind the highest surface tempe ture by 29 days. The thawing duration of ZL05 is slightly longer than the freezing du tion, which indicates that its permafrost has a certain degree of degradation. Figure  shows that, on the regional scale, the date that MFH occurred was mainly distributed fro 27 November to 21 March of the following year, with an annual average of date of ye (DOY) 37, which is about one month earlier than the onset time of thawing. The date th MTS occurred was mainly distributed from 25 July to 21 September of the following ye with an annual average of DOY 225, about 39 days earlier than the onset time of freezin According to the classification results of the date of MFH and the date of MTS, the int mountain basin and lake basin in relatively low altitude areas, the date of MFH is som what earlier than the piedmont plain, while the date of MTS is relatively later than t piedmont plain.  Figure 3 shows that, in the ZL05, the date of MFH lags behind the lowest surface temperature by 35 days, while the date of MTS lags behind the highest surface temperature by 29 days. The thawing duration of ZL05 is slightly longer than the freezing duration, which indicates that its permafrost has a certain degree of degradation. Figure 7 shows that, on the regional scale, the date that MFH occurred was mainly distributed from 27 November to 21 March of the following year, with an annual average of date of year (DOY) 37, which is about one month earlier than the onset time of thawing. The date that MTS occurred was mainly distributed from 25 July to 21 September of the following year, with an annual average of DOY 225, about 39 days earlier than the onset time of freezing. According to the classification results of the date of MFH and the date of MTS, the intermountain basin and lake basin in relatively low altitude areas, the date of MFH is somewhat earlier than the piedmont plain, while the date of MTS is relatively later than the piedmont plain.

Time Characteristics of Seasonal Deformation Processes
The average thawing duration in the study area is 195 days, larger than the thawing duration calculated based on MODIS LST (average 156 days). In most areas (within three standard deviations), the thawing duration is 175~225 days. The shortest thawing duration is 104~160 days in the high-altitude mountain areas at the edge of the study area and the bedrock at the lower altitude in the middle. In contrast, the longest thawing duration is more than 260 days at the glacier terminus. The thawing duration is longer in the piedmont alluvial-pluvial fan, the middle of the river channel, and the seasonal thawing period around the relatively flat mountain basin, river valley, and lake basin. The average thawing duration in the study area is 195 days, larger than the thawing duration calculated based on MODIS LST (average 156 days). In most areas (within three standard deviations), the thawing duration is 175~225 days. The shortest thawing duration is 104~160 days in the high-altitude mountain areas at the edge of the study area and the bedrock at the lower altitude in the middle. In contrast, the longest thawing duration is more than 260 days at the glacier terminus. The thawing duration is longer in the piedmont alluvial-pluvial fan, the middle of the river channel, and the seasonal thawing period around the relatively flat mountain basin, river valley, and lake basin.  Figure 8 shows a specific negative correlation between SDA and altitude, similar to previous research results [71,74]. Moreover, there is an apparent correlation between seasonal deformation amplitude and the time characteristics of the F-T process in AL.

Relationship between Amplitude and Time of Seasonal Deformation
According to the correlation analysis (Table 4), the result shows a strong negative correlation between SDA and altitude, and there is also an obvious correlation with seasonal deformation time characteristics. Therefore, the correlation of seasonal deformation characteristics of all pixels in space is analyzed. The correlation of SDA, date of MFH, date of MTS, and thawing duration all passed the significance test of 0.001. There is a strong correlation coefficient between thawing duration and the date of MFH and MTS, mainly determined by its calculation method. SDA is negatively correlated with the date of MFH and thawing duration while positively correlated with the date of MTS.  Figure 8 shows a specific negative correlation between SDA and altitude, similar to previous research results [71,74]. Moreover, there is an apparent correlation between seasonal deformation amplitude and the time characteristics of the F-T process in AL. According to the correlation analysis (Table 4), the result shows a strong negative correlation between SDA and altitude, and there is also an obvious correlation with seasonal deformation time characteristics. Therefore, the correlation of seasonal deformation characteristics of all pixels in space is analyzed. The correlation of SDA, date of MFH, date of MTS, and thawing duration all passed the significance test of 0.001. There is a strong correlation coefficient between thawing duration and the date of MFH and MTS, mainly determined by its calculation method. SDA is negatively correlated with the date of MFH and thawing duration while positively correlated with the date of MTS.  Figure 9 shows that except for latitude and longitude, the results reveal that all factors passed the significance test (p < 0.005). The explanation ratios of factors are different. DDF and MAGT contribute more than 60% to the seasonal deformation in the study area. The factors that contribute more than 40% to the seasonal deformation include mean annual air temperature, the ratio of DDF to DDT, DDT, and ALT. The soil organic carbon, precipitation, NDVI, and ground ice content contribute more than 20% to the seasonal deformation. The contribution of soil texture, slope direction, slope, and vegetation type to seasonal deformation is relatively low (less than 10%), among which the vegetation form is only 2%, mainly because the study area is an inland lake area with a relatively simple landform.   Figure 9 shows that except for latitude and longitude, the results reveal that all factors passed the significance test (p < 0.005). The explanation ratios of factors are different. DDF and MAGT contribute more than 60% to the seasonal deformation in the study area. The factors that contribute more than 40% to the seasonal deformation include mean annual air temperature, the ratio of DDF to DDT, DDT, and ALT. The soil organic carbon, precipitation, NDVI, and ground ice content contribute more than 20% to the seasonal deformation. The contribution of soil texture, slope direction, slope, and vegetation type to seasonal deformation is relatively low (less than 10%), among which the vegetation form is only 2%, mainly because the study area is an inland lake area with a relatively simple landform. Figure 5 demonstrates the seasonal deformation variance at four borehole sites, and Figure 6 displays the spatial heterogeneity of seasonal deformation. According to the field investigation, these four boreholes are all located in river valleys, in flat terrain (slope ≤ 1 • ), vegetation coverage is above 85%, and ALT is around 1.8 m ( Table 2). When other conditions are similar, the soil moisture of boreholes HL01, HL02, and HL03 is identical and significantly higher than that of borehole HL04. The vegetation development is closely connected with water supplies in the study area and they are good indicators of soil water content [94], e.g., the swamp meadow environment has much higher soil water content than the meadow or desert environment. Since the soil water content in the AL is not available in the entire study area, we used NDVI as a substitute index of soil water content to analyze the relationship between soil water content and SDA. Moreover, the ALT data also was used to analyze the relationship between ALT and SDA. We analyzed the relationship between various factors and SDA and found that only NDVI and ALT had the most apparent effect on the SDA. Figure 10 shows that the SDA increases obviously with the increase in NDVI and ALT.  Figure 5 demonstrates the seasonal deformation variance at four borehole sites, and Figure 6 displays the spatial heterogeneity of seasonal deformation. According to the field investigation, these four boreholes are all located in river valleys, in flat terrain (slope ≤ 1°), vegetation coverage is above 85%, and ALT is around 1.8 m (Table 2). When other conditions are similar, the soil moisture of boreholes HL01, HL02, and HL03 is identical and significantly higher than that of borehole HL04. The vegetation development is closely connected with water supplies in the study area and they are good indicators of soil water content [94], e.g., the swamp meadow environment has much higher soil water content than the meadow or desert environment. Since the soil water content in the AL is not available in the entire study area, we used NDVI as a substitute index of soil water content to analyze the relationship between soil water content and SDA. Moreover, the ALT data also was used to analyze the relationship between ALT and SDA. We analyzed the relationship between various factors and SDA and found that only NDVI and ALT had the most apparent effect on the SDA. Figure 10 shows that the SDA increases obviously with the increase in NDVI and ALT.   Figure 5 demonstrates the seasonal deformation variance at four borehole sites, and Figure 6 displays the spatial heterogeneity of seasonal deformation. According to the field investigation, these four boreholes are all located in river valleys, in flat terrain (slope ≤ 1°), vegetation coverage is above 85%, and ALT is around 1.8 m ( Table 2). When other conditions are similar, the soil moisture of boreholes HL01, HL02, and HL03 is identical and significantly higher than that of borehole HL04. The vegetation development is closely connected with water supplies in the study area and they are good indicators of soil water content [94], e.g., the swamp meadow environment has much higher soil water content than the meadow or desert environment. Since the soil water content in the AL is not available in the entire study area, we used NDVI as a substitute index of soil water content to analyze the relationship between soil water content and SDA. Moreover, the ALT data also was used to analyze the relationship between ALT and SDA. We analyzed the relationship between various factors and SDA and found that only NDVI and ALT had the most apparent effect on the SDA. Figure 10 shows that the SDA increases obviously with the increase in NDVI and ALT.    Table 5 lists the SDA retrieved in this study and values in the previous studies. Before 2010, the SDA was in the range of 0~20 mm, as revealed both by the L-band data of ALOS PALSAR [110] and the C-band data of ENVISAT ASAR [85]. The SDA values derived by Sentinel-1 after the 2010s are larger than before the 2010s. Our results are similar to Chen's [71]. The seasonal deformation of permafrost acquired by MT-InSAR technology is influenced by many factors such as radar data (mainly wavelength), MT-InSAR algorithm, and model (which is used to separate the seasonal deformation from surface deformation), and monitoring period.

Influencing Factors of Seasonal Deformation Amplitude
(1) Wavelength. The longer the wavelength of the radar, the stronger the penetrating ability, and the less susceptible it is to vegetation [61,113]. The vegetation in permafrost regions of QTP is mainly alpine grassland and alpine desert [94,114]. Previous research results show that, compared with permafrost regions around the Arctic, the wavelength of radar signal has relatively little influence on obtaining the surface deformation of permafrost regions of the QTP [62,115]. (2) MT-InSAR algorithm. The Stamps-InSAR and SBAS-InSAR algorithms are based on permanent scatterers [116] and distributed scatterers [117]. They have been widely used to obtain ground deformation in permafrost regions [60,78,118,119], and their reliability and accuracy have been proved. Because of the requirements for the scattering characteristics in the permafrost region, the SDA obtained by the two algorithms may be different. (3) Model. Different models may lead to differences in the calculated SDA. The seasonal deformation is caused by the F-T cycle in AL and is determined by many factors, such as topography, soil texture, and soil moisture [62,66,68,71,74]. Both the physical model and statistical model have some shortcomings in obtaining seasonal deformation of AL in permafrost areas [95,119,120]. (4) Monitoring period. For Sentinel-1 SAR data with the same wavelength and MT-InSAR algorithm, the SDA obtained by SBAS-InSAR [71,75,112,121] is somewhat different, and the SDA increases with period generally. Previous studies have shown that permafrost in the study area shows an evident degradation trend [93], which will lead to the thawing of ground ice at the bottom of the AL and the increase in the ALT [18]. An increase in the soil water content in the AL leads to an increase in the SDA.
In this study area, the SDA tends to increase with time in the study area, which may be related to permafrost degradation and ground ice thawing. Our research and Chen [71] obtained the SDA of AL based on Sentinel-1 SAR and SBAS-InSAR. Although the model is different, the results are consistent with the observed data. The ALT change in the QTP was 19.5 cm per decade from 1981 to 2018 [122], and the precipitation was about 300 mm [123] in the study area. The interannual variation in the ALT and the soil moisture change in AL may be the undetectable difference between the SDA in this study and Chen's [71].

Influencing Factors of Seasonal Deformation Amplitude
Regarding SDA in the northwest of QTP, unsaturated soil water is mainly concentrated at the bottom of AL and is relatively small, which was 2.5~12 mm from 2003 to 2011 [74]. However, the SDA range from 2007 to 2011 was 0.5~28 mm in the southern boundary of permafrost because the AL is thick and soil moisture content is higher [64]. Previous studies have shown that the SDA positively correlates with the ALT in the same vegetation form [70,71,90]. In the alpine swamp area located in the south of the study area, the SDA ranged from 0 to 74 mm from 2018 to 2019 [77], while in the alpine desert area located in Wudaoliang, the SDA ranged from 2014 to 2019 and was only 9~26 mm [111], which indicated that the influence of vegetation on the SDA was also very prominent. The increased NDVI may be related to the increase in soil water content in the uppermost AL, leading to a larger SDA. Consistent with previous studies [62,111], our study also confirms that SDA at the foot of the slope is also larger than that of the slope, mainly related to the soil water content in AL.
The research in the central of QTP shows that slope and latitude are the main influencing factors of SDA on a large scale [71]. However, our study reveals that the SDA's main influencing factors are DDF and MAGT. It reflects that the main influencing factors are relevant to the studying scales. In addition, similar to other studies, our results also show that the SDA is related to the vegetation cover and the soil water content in AL and ALT. In addition, similar to other studies, our results also show that the SDA is positively correlated with the vegetation cover [70,71,79] and soil water content of AL [8,74] and the ALT [71].

Time Characteristics of Seasonal Deformation
For the first time, the time characteristics of seasonal deformation of AL were obtained by SBAS-InSAR and analyzed in detail. It is different from the characteristics of surface F-T time calculated on MODIS LST.
Compared with the F-T time calculated by MODIS LST, which represents the land surface F-T state, the F-T cycle process of the AL is more complicated. In addition to the factors, e.g., climate, altitude, latitude, and surface conditions that affect the surface F-T cycles [5,29,124], the F-T cycles in the AL are also affected by soil texture and soil moisture [5,78]. The F-T cycle of the AL includes four stages: spring warming, summer thawing, autumn freezing, and winter cooling [125]. It is observed that the date of MTS and the date of MFH is earlier than the onset of surface freezing and surface thawing, respectively, and the thawing duration is smaller than the thawing period. Therefore, the time characteristics of the F-T cycle obtained based on seasonal deformation may better reflect the F-T process of AL.
Although existing studies have found that the seasonal deformation time of AL is obviously different from air temperature or surface temperature [74,84], the research on obtaining the time characteristics of the F-T cycle of AL using InSAR technology is still relatively scarce. In Svalbard in the Arctic, Rouyet et al. [78] used MT-InSAR technology and Sentinel-1 data for the first time to obtain the DOY subsidence maxima of AL. Therefore, by comparing their research, we discussed the DOY difference between the QTP and the circum-Arctic region when the AL reached the subsidence maxima. It should be noted that, since the research results of Rouyet et al. are mainly concentrated in flat areas (the slope is less than 2 • ), we only compared the DOY of MTS in flat areas, and the results are shown in Table 6. The DOY of MTS time of the AL in the study area is earlier than that of Svalbard, located in the circum-Arctic region, but it is closest to the Adventdalen area in Svalbard (the difference is about 10 days). In the Adventdalen area and the flat area of this study area (slope less than 2 • ), the ALT is mainly distributed in 0.9~1.8 and 1.4~3.8 m, while the annual average surface temperature is −0.5~−3.2 and −3.6~−4.75 • C, respectively. The annual average surface temperature and ALT are the critical factors of climate influence on permafrost. The higher the annual average surface temperature, the thinner the ALT and the earlier DOY of MTS of AL. In the thawing season, the migration speed of the thawing front in the AL is also affected by the soil texture and soil water content in the AL [126,127]. The higher the soil water content of AL, the slower the migration speed of the thawing front, while the soil texture affects the migration speed of the thawing front through thermal conductivity, and the higher the thermal conductivity, the faster the migration speed of the thawing front. Due to the lack of data on soil water content and soil texture of AL in the Adventdalen area, we can only speculate that the reason why the DOY of MTS in QTP is earlier than that in the circum-Arctic area is mainly due to the difference in active layer thickness caused by climatic conditions.

Relationship between Amplitude and Time of Seasonal Deformation
Seasonal deformation characteristics exhibit large spatial heterogeneity in the northern area of Kusai Lake and Yanhu Lake ( Figure 11). Four regions of interest (delineated with red lines in Figure 11) were explored to analyze their causes. According to the visual interpretation of remote sensing optical images and field investigation, these four regions are all scoured by seasonal rivers. Region A is mainly a sandy slate with an average slope of 5 • . The slope of region B is relatively small (the average slope is 2 • ), which belongs to the alluvial accumulation area. More than 20 spring groups are in the south of B. Regions C and D are mainly sandy loams located in the upper reaches of the alluvial delta of seasonal rivers and the lower reaches of the alluvial/diluvial delta near lakes. In regions B, C, and D, some seasonal rivers seep into the undercurrent in the alluvial accumulation area from the outlet of the mountain pass to the front section of the inclined plane of the lake basin. Moreover, some seep into groundwater again when they reach the outlet of the exposed spring, overflow, return to surface runoff, and then seep into the lakeside section from the exposed spring.  Figure 12 shows that the date of MFH and MTS in the four regions is earlier than in the surrounding areas, but the thawing duration and SDA are smaller than in the surrounding areas. Due to the continuous supply of seasonal rivers in the thawing period, the average soil moisture in four typical areas will be relatively high [128]. However, in the freezing period, when there is no continuous water supply shortage, the soil moisture in four typical areas is lower than in the surrounding areas due to the poor soil water holding capacity. The above two reasons lead to the earlier date of MFH and MTS in four regions to their surrounding areas and the shorter thawing duration and smaller SDA.  Figure 12 shows that the date of MFH and MTS in the four regions is earlier than in the surrounding areas, but the thawing duration and SDA are smaller than in the surrounding areas. Due to the continuous supply of seasonal rivers in the thawing period, the average soil moisture in four typical areas will be relatively high [128]. However, in the freezing period, when there is no continuous water supply shortage, the soil moisture in four typical areas is lower than in the surrounding areas due to the poor soil water holding capacity.
The above two reasons lead to the earlier date of MFH and MTS in four regions to their surrounding areas and the shorter thawing duration and smaller SDA. and thawing duration in the north of Kusai lake and Yanhu lake. Figure 12 shows that the date of MFH and MTS in the four regions is earlier than in the surrounding areas, but the thawing duration and SDA are smaller than in the surrounding areas. Due to the continuous supply of seasonal rivers in the thawing period, the average soil moisture in four typical areas will be relatively high [128]. However, in the freezing period, when there is no continuous water supply shortage, the soil moisture in four typical areas is lower than in the surrounding areas due to the poor soil water holding capacity. The above two reasons lead to the earlier date of MFH and MTS in four regions to their surrounding areas and the shorter thawing duration and smaller SDA.

Conclusions
We used Sentinel-1 SAR data from 2014 to 2019 to obtain the characteristics of the F-T process of AL from 2017 to 2019 in the ZL-YL basin of the QTP by SBAS-InSAR technology.
The SDA of the AL is about 0~60 mm, and the average SDA is 19 mm. In the study area, many factors jointly control SDA, and the thermal condition of permafrost is the most critical factor affecting the SDA, in which MAGT and DDF have more than 60% influence on the seasonal deformation. In the study area, the thermal condition of permafrost is the main factor in the spatial distribution pattern of SDA, and the SDA has the most apparent relationship with the soil moisture, NDVI, and ALT.
The time characteristics of the F-T cycle show that the date that MFH occurred was mainly distributed from 27 November to 21 March of the following year, with an average time of DOY 37, while the date that MTS occurred was mainly distributed from 25 July to 21 September, with an average time of DOY 225. The average thawing duration of the AL is 193 days, which is about one month longer than the thawing duration of the land surface. SDA is negatively correlated with the date of MFH and thawing duration, while it is positively correlated with the date of MTS. The time characteristics of the F-T cycle obtained based on seasonal deformation may better reflect the F-T process of AL.
In addition, in the alluvial and flood accumulation areas with seasonal river scouring, the soil texture significantly influences the SDA through soil water holding capacity. It mainly shows that the SDA is smaller than its surrounding areas, the date of MFH and MTS is earlier than its surrounding areas, and the thawing duration is shorter than its surrounding areas. Table 7 shows the abbreviations frequently used in this study and their corresponding full names.

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