Active Layer Thickness Retrieval Over the Qinghai-Tibet Plateau Using Sentinel-1 Multitemporal InSAR Monitored Permafrost Subsidence and Temporal-Spatial Multilayer Soil Moisture Data

Increasing near-surface temperature over the Qinghai-Tibet Plateau (QTP) has led to permafrost degradation and increasing active layer thickness (ALT). In this study, the ALT was estimated based on ground subsidence monitored by multitemporal interferometric synthetic aperture radar (MT-InSAR) and temporal-spatial multilayer soil moisture data. For the ground subsidence monitoring, a modified Stefan piecewise elevation change model based on air temperature data was integrated into a new small baseline subset (NSBAS) chain. A total of 33 scenes of Sentinel-1 data (S-1) were collected over one year to build the MT-InSAR analysis network. Moreover, both soil moisture active/passive (SMAP) L4 surface and root zone soil moisture data and ERA-Interim reanalysis data were used to build an ALT retrieval model. In particular, the global-scaled soil moisture data (SMAP and ERA-Interim) fraction was separated based on the Sentinel-1 amplitude-based land cover classification results and in situ soil moisture data. A typical ALT estimation method based on the point scale groundwater information was also performed to evaluate the performance of the proposed method. Based on the validation of the ground-based ALT observations, the proposed method outperformed the traditional point scale groundwater information-based method, with a correlation coefficient of 0.67, RMSE of 0.70 and ubRMSE of 0.51, respectively. The ERA-Interim-based estimation results were underestimated due to the overestimation of the ERA-Interim soil moisture data. Obvious differences were observed between the ALT of the alpine meadow areas and alpine desert areas. Our results demonstrate that the combination of temporal-spatial multilayer soil moisture data and the MT-InSAR method with S-1 images is a promising approach for the large-scale characterization of ALT.

has become widespread due to continuous climatic warming [4]. The degradation of permafrost is associated with local and widespread collapse, instability of the ground surface, and surface subsidence [5]. The ALT is a unique indicator of climate change because this factor is affected by air temperature changes, soil moisture variations and underground heat transfer [6]; the ALT is a significant factor for monitoring the permafrost degradation of the QTP [7], [8]. However, the observations of the ALT of the QTP are limited because of the cost of obtaining direct measurements and maintaining active monitoring programs in permafrost regions [9]. Therefore, it is crucial to develop practical methods to retrieve the ALT over large spatial areas.
Traditional field-based analytical methods for estimating ALT are usually based on ground temperature monitoring, frost or thaw tubes and mechanical probing [5], [10]- [12]. However, these conventional methods are extremely laborintensive, time-consuming and limited to sparse point densities [8]. Ground penetrating radar (GPR) is an efficient and powerful tool for analyzing permafrost distribution characteristics, ALT and evolution processes. Compared with the sparse point densities of the traditional field methods, GPR provides line profile explorations of the underground geophysical status to a maximum exploration depth of 4 m [13]. Many semiempirical model-based analytical methods have been developed to estimate ALT using temperature, soil water content, soil property data and heat transfer equations [5]. Among all these methods, the Stefan [14] and Kudryavtsev models [15] are the most widely used for the QTP and other high latitude permafrost regions [16]- [18]. The Stefan or simplified Stefan model combining air temperature and soil parameters as input data has been applied to point and regional scales [19]- [23]. Certain studies have combined physically based parameters with the Kudryavtsev solution to improve the simulation precision [24]. Additionally, stochastic approaches, such as Bayesian models, are effective methods for ALT and maximum thickness of seasonally frozen ground (MTSFG) estimations. Qin et al. developed a Bayesian model to evaluate the MTSFG based on temperature and precipitation data [9]. However, these point and regional scale-based methods are limited by low spatial resolution and have difficulty meeting the application requirements at a fine scale. In addition, numerical model, which incorporating phase change effects has been regarded as effective method for simulating and forecasting the thermal regime of permafrost over a short time interval [25]. These models require detailed surface meteorological variables, and soil thermal parameters [26]. For example, geophysical institute permafrost laboratory 2 (GIPL2) model is a numerical transient model that has been widely used to simulate thermal conditions of the permafrost including ALT and MTSFG [26], [27].
Interferometric synthetic aperture radar (InSAR) and multitemporal InSAR (MT-InSAR) techniques have been widely applied in monitoring surface displacements to identify terrain stability issues in QTP permafrost regions [28]- [34]; these studies showed that InSAR and MT-InSAR were suitable for documenting the slope movement processes and studying the seasonal displacement dynamics related to ground freezing and thawing, and the seasonal displacement was successfully used to estimate the ALT induced by the freeze-thaw process [8], [35]- [37]. The resolution of the ALT mapped by a synthetic aperture radar (SAR) image is far superior to global scale results, which make this technique promising for monitoring changes in ALT [37]. The key factors of the ALT retrieval were to accurately estimate the seasonal displacement and to distinguish a reasonable relationship between the seasonal deformation process and the ALT. Regarding seasonal deformation estimation, MT-InSAR methods perform better than conventional InSAR methods [38], and choosing an appropriate seasonal deformation model is crucial to MT-InSAR analysis over permafrost regions [39]. Currently, the representative seasonal deformation models used in permafrost regions are simplified or modified Stefan model-based seasonal deformation models [4], [29], [35], sinusoidal seasonal models [34], [39]- [42] and climatic factor-based season deformation models [31]; the relationship between the MT-InSAR-derived surface subsidence and ALT are mainly distinguished by the vertical distribution of the pore water within the active layer [8], [36], [41] because during seasonal freezing and thawing, the water-ice phase changes can induce cm-scale uplift and subsidence [43]. Liu et al. [44] proposed a vertical distribution model for the water content within the active layer to convert the seasonal displacement to ALT under the assumption that the active layer is fully saturated. Wang et al. [8] built an ALT retrieval model that combined the ground subsidence using the InSAR technique, types of land cover and ground water information. A multilayered groundwater model for different types of land cover was proposed. Li et al. proposed a new method based on the heat transfer process of soil and the time lags between the InSAR measured amplitude of the seasonal oscillations and meteorologically recorded temperatures [37]. However, these studies did not consider the temporal-spatial distribution of the soil water content but assumed that the temporal soil moisture is homogeneous and the spatial variation of the soil moisture is only discriminated at the point scale.
In this paper, an ALT estimation model based on the ground subsidence derived from MT-InSAR and temporal-spatial multilayer soil moisture data is proposed. First, time series multilayer ERA-Interim reanalysis and SMAP L4 soil moisture data were acquired from the European Centre for Medium-Range Weather Forecasts (ECMWF) and National Aeronautics and Space Administration (NASA); due to the limitations in the resolution of these data, the soil moisture data of different types of land cover were separated by the ratio of the in situ soil moisture statistics. The land cover classification was performed using an S-1 amplitude map. Then, the temporal-spatial profile of the soil moisture distribution was established. Second, the new SBAS (NSBAS) technique with a modified Stefan piecewise elevation change model was developed to monitor the freeze-thaw seasonal displacement. A total of 33 S-1 images taken at intervals as short as 12 days from 28 November 2017 to 29 December 2018 were collected to build the MT-InSAR analysis network. Finally, combining the seasonal deformation estimated from the MT-InSAR and the temporal-spatial profile of the soil moisture distribution, the ALT was retrieved. Compared with previous studies, our method considers the variation in the temporal-spatial soil moisture at different soil depths.

A. STUDY AREA
The study area is located in the south of the Qinghai province (longitude: 90.715−93.751, latitude: 34.204−35.836) from Tuotuohe to Wudaoliang (Figure 1a). This area is a warm continuous permafrost region with widespread seasonal frozen soil. The elevation ranges from 4400 to 5455 m ( Figure 1b). Recorded by the Beiluhe Weather Station, the mean annual precipitation is approximately 369.8 mm/yr. Over 90% of the precipitation falls between May and September [45].
A previous study indicated that the thickness of the permafrost ranges from 50 to 80 m and the upper interface of the permafrost ranges from −1.8 to −2.2 m [46]. The active layer thickness (ALT) in this area varies from 1.6 to 3.4 m [10], [45]. The ALT has increased over the past 10 years, with notable trends in desert grassland (5.3 cm/year) and degrading alpine meadow (4.4 cm/year) areas [45]. The overall mean ALT increased at a rate of 6.3 cm/yr from 2006-2010 [41]. The main types of vegetation in the permafrost area are alpine meadow, alpine steppe, and alpine desert steppe [8], [47]. The study area is same as that of our previous study; more detailed information about the study area can be found in [39].

2) MULTILAYER SOIL MOISTURE DATA
In this study, two types of soil moisture data were used to build the ALT estimate model. One type is SMAP level 4 surface and root zone soil moisture data from NASA, distributed over a 9 km (0.1 • ) grid. The SMAP surface measurements provide direct sensing of the soil moisture in the top 5 cm of the soil column and produce model-derived estimates of the root zone soil moisture at a depth of 1 m. Such estimates are obtained by merging the SMAP observations with estimates from a land surface model in a data assimilation system [48]. These data are available at https://nsidc.org/. The other type of data used is the ERA-Interim reanalysis data, produced by the ECMWF integrated forecasting system model and providing a four-layer representation of soil (0-7 cm, 7-21 cm, 21-72 cm, 72-189 cm) distributed over a 13.187 km (0.125 • ) grid [49]. These data are available from https://apps.ecmwf.int.

3) IN SITU SOIL MOISTURE AND ALT DATA
The in situ soil moisture data were used to separate the soil moisture fractions of the SMAP and ERA-Interim soil moisture grid by the different types of land cover. Five field campaigns in Beiluhe were performed on 12 August 2015, 29 July 2016, 9 August 2016, 28 August 2018 and 31 August 2018, and a total of 87 samples were collected. The sampling locations are displayed in Figure 2. The time domain reflectometer (TDR) and gravimetric soil moisture measurements were used to measure the volumetric soil moisture with 5 cm depth. The correlation between the gravimetrically determined soil moisture and TDR measured in this area were compared, and the confidence was suggested [50]. Meanwhile, the groundwater content of the alpine meadow and alpine desert with depth of 5,20,50,80,120,150,180 cm from May 2014 to August 2016 were collected (in situ ALT A2 in Figure 2b) [8].
The in situ ALT data were measured by GPR for 28 August 2018 (in situ ALT A1 in Figure 2b). A detailed introduction to the method can be found in [39]. Two in situ monitoring sites with different types of land cover measured by Wang et al. [8] (in situ ALT A2 in Figure 2b), the filed investigation of ALT measured by Lu et al. in June 2015 [46] (in situ ALT A3 in Figure 2b) and 10 instrumented borehole ALT measured by Yin et al. [45] (in situ ALT A4 in Figure 2b) from 2014 to 2016 were also used in this study.

4) ANCILLARY DATA
In this study, certain auxiliary data were used to build the seasonal deformation model and analyze the ground subsidence and ALT results. The modified Stefan seasonal deformation model was built based on the air temperature data. We obtained the 2-m air temperature ERA-Interim reanalysis data from ECMWF. This parameter has been proven to resolve the near-surface air temperature for various mountain regions well [4]. In addition, the influencing factor analysis of the surface feature corresponding to the multilayer soil temperature data was evaluated. The multilayer soil temperature data have four layers of soil (0-7 cm, 7-21 cm, 21-72 cm, 72-189 cm), which are produced by the ECMWF integrated forecasting system model.

III. METHODOLOGY
The objective of this study was to estimate the ALT based on the multilayer soil moisture data and freeze-thaw seasonal subsidence results monitored by MT-InSAR. The processing flow chart consists of four main steps: 1) InSAR preprocessing of the S-1 images, 2) seasonal subsidence modeling, 3) NSBAS deformation estimation, and 4) ALT conversion modeling. An overview of the proposed method is shown in Figure 3.

A. SENTINEL-1 INSAR PROCESSING
The InSAR preprocessing included interferogram network selection, coregistration, differential interferogram phase generation, phase unwrapping and phase correction. First, the interferogram network was selected based on the perpendicular-temporal baseline, and then, all the SAR images were coregistered; after the coregistration, the interferogram network was refined based on the coherence, and a total of 77 interferometric pairs were used to generate the interferograms. The topographic phase was removed using a 1-arcsecond grid (30 m) shuttle radar topography mission (SRTM) digital elevation model (DEM) [51]. Multilooking with 2 pixels in range and 10 pixels in azimuth was performed to remove the phase noise. A minimum cost flow (MCF) phase unwrapping technique was used for unwrapping the phase [52]. These processes were implemented based on the generic mapping tools (GMTSAR) [53]. Then, the interferometric atmospheric phase delay and orbital error corrections were implemented by PyAPS [54] and network deramping methods [55]. The GIAnT [56], [57] was used to perform the phase correction. A detailed description of the InSAR processing can be found in [39].

B. SEASONAL AND LONG-TERM SUBSIDENCE MODEL
In permafrost regions, subsidence and uplift are dominantly associated with the freezing and thawing of the active layer VOLUME 8, 2020 and uppermost permafrost [4], [8]. The deformation process cannot be monitored using a linear model. In this paper, a piecewise deformation change model [4] containing long-term and seasonal deformation was applied as follows: where D is the total ground displacement on any day t since the onset of thawing, V is the long-term displacement rate (mm/yr), t T and t F are the onsets of thawing and freezing, respectively; according to the air temperature data, the thawing and freezing onsets were 9 April 2018 and October 6, 2018, respectively; T is the number of consecutive permafrost years, S is the seasonal subsidence rate (mm), and I is the composite index to combine the thaw and freeze indices as given by the following: where √ A T (t) and √ A F (t) are the square roots of the accumulated degree days of thaw (ADDT) and freeze (ADDF) (units: • C days), respectively, calculated using the air temperatures. In this study, we calculated ADDT and ADDF using the daily averaged 2-m air temperatures acquired from ECMWF based on the Stefan equation [58]. The scaling factor α can be expressed as follows: where k F = 1.4 W · M −1 · K −1 and k T = 0.6 W · M −1 · K −1 are the thermal conductivities estimated by Romanovsky and Osterkamp [59] and n F and n T are the n-factors (n F = 0.61, n T = 0.62, calculated by Cao et al. [3]) indicating the ratio between the ground surface temperature and air temperature [4].

C. NSBAS METHOD
In this paper, considering the effects of the thawing and freezing of the permafrost and the time span of the observations, the seasonal model expressed in (1) was used as follows: where λ is the wavelength, B ⊥ denotes the perpendicular baseline, θ and R denote the incidence angle and slant range distance, respectively, and t represents the time interval between two SAR scenes.
In recent studies, The seasonal subsidence with a short time interval was successfully monitored based on the PSI (StaMPS-InSAR) [60] or DInSAR [61] method using Sentinel-1 data. The linear model with a short time interval used in seasonal subsidence monitoring is acceptable. However, in the StaMPS-InSAR method, the PS points were insufficient especially over the alpine meadow areas, while in the traditional DInSAR method, the absolute deformation has to be estimated using seasonal function absolute deformation model. The applicability of the seasonal deformation model based NSBAS method was well demonstrated in our pervious published study [62], to further improve the accuracy of seasonal deformation monitoring, the modified Stefan model was proposed and integrated into the NSBAS chain in this study. The NSBAS method was proposed and described by López-Quiroz et al. [63], and the method chain described in [64] which can be applied over pixels where critical interferometric links are missing and is effective for conditions where strong spatial-temporal decorrelation occurs in an interferometry network [64]. In traditional SBAS, the linear relationship between the set of interferometric phase observations and individual SAR scene phase values is as follows: where M is the number of interferograms, i and j denote the acquisition time, l is a specific pixel, d l is the vector including the interferometric phase ( I ), m l is a vector containing the phase delay increments between two successive images, and G l is a matrix of zeros and ones directly related to the set of interferograms generated from the available data. The singular value decomposition (SVD) method was used to solve (5). To overcome the SVD biases when independent images do not overlap, NSBAS uses additional constraints as follows [63]: where N denotes the acquisition dates, B k ⊥ is the perpendicular baseline between the satellite paths at acquisitions 1 and k, and f(·) is a parametric representation of the temporal form of the deformation, used as a regularization function. For a seasonal deformation surface, the parameter f ( t k ) in (6) can be replaced by D ( t k ) from (1). Then, the constrained system can be solved by a least-squares inversion: where − T line is the time interval of the linear deformation term.
In this study, the GIAnT [56], [57] was used to integrate the NSBAS and the seasonal deformation model. After the coefficients V and S are obtained, the time-series surface elevation changes can be deduced from these parameters.

D. ALT RETRIEVAL MODEL
The ALT retrieval over the permafrost area is based on the assumption that the surface subsidence is caused by the change in the volume of the pore ice transitioning to water in the active layer. The relationship between the displacement of the permafrost and ALT can be expressed as follows [44]: where Mv is the soil moisture, ρ w and ρ i are the densities of pure water and ice, respectively, H is the ALT, and z is the seasonal displacement. The soil moisture is not saturated in the thawing season and varies temporally and with different soil depths and different types of soil and ecosystems. Therefore, the temporal-spatial multilayer soil moisture data were introduced for the ALT retrieval. The soil moisture data used in this study are distributed over 9 km (SMAP L4) or 13 km (ERA-Interim) resolution; it is hard to separate the soil moistures of different types of land cover in a single pixel. In this study, the in situ soil moisture data were measured in different periods (Figure 2). The soil moisture ratio between the alpine meadow and alpine desert areas was computed according to the classification results of an S-1 amplitude image (Tables 1 and 2). Then, the fractions of the soil moisture contributions of different land covers in the ERA-Interim or SMAP pixel were determined based on the ratio. The mixed pixels were separated to different endmembers at a coarse resolution.
The time-series ERA-Interim soil moisture data contains four layers: 0-7 cm, 7-21 cm, 21-72 cm, 72-189 cm, and the SMAP L4 soil moisture comprises the surface (0-5 cm) and root zone (1 m) layers. The piecewise multilayer retrieval model can be built based on these data and the MT-InSAR observed surface subsidence as follows: VOLUME 8, 2020  where H is the ALT throughout thawing season, N is the number of intervals, H i is the thawing depth of the permafrost in a specific time interval, z i is the difference in the seasonal displacement amplitude in successive S-1 acquisitions, j is the specific soil moisture layer; the selection of j is based on the accumulation of H i ; and Mv i,j denotes the mean soil moisture between successive time intervals (daily) in a specific layer. Considering the ERA-Interim soil moisture data, the layer of 72-189 cm was chosen as surrogate data when the accumulated ALT was greater than 189 cm. In terms of the SMAP L4 soil moisture data, the surface soil moisture data were selected as the first layer when the accumulation of H i was less than 1 m, while the root zone soil mois-ture layer was chosen when the melting depth was greater than 1 m.

A. SPATIAL-TEMPORAL SEASONAL SUBSIDENCE
The estimated seasonal subsidence of the study area is presented in Figure 4a. The seasonal settlement ranged from -79.98 to 5 mm, while in the profile zones A1-A2 and B1-B2 (Figure 4a), the seasonal settlement ranged from -60 to 0 mm ( Figure 5a). Figure 4b depicts the map of the linear settlement, where the linear subsidence trends were as high as 29.9 mm/yr inside the study area, and the linear subsidence values in most areas were approximately 10 mm/yr, which is similar to previous results measured by [8]. The DEM errors inferred from Figure 4c ranged from −50 to 50 m in most of the areas studied, and the mean and standard deviation acquired were similar to those of the sinusoidal seasonal model-based results [39].
To evaluate the accuracy of the seasonal subsidence results, a correlation analysis was performed between the sinusoidal seasonal model-based results (Figure 4d) and the seasonal subsidence model used in this study (Figure 4a), with a correlation coefficient of 0.7 (Figure 5a). In our previous study, the spatial distribution and variability characteristics 84342 VOLUME 8, 2020 of the seasonal amplitude based on the sinusoidal seasonal model were observed, and opposing fluctuations were found between the elevation and seasonal amplitude [39]. The correlation coefficient between the slope of the topography and the seasonal amplitude over the profile zones A1-A2 and B1-B2 was 0.54 ( Figure 5b). In this study, an improved correlation was found with a coefficient of 0.6 ( Figure 5c). This improved coefficient demonstrated that the seasonal subsidence model used in this study is more suitable for characterizing the surface displacement. Figure 6 shows the time-series displacement of selected statistics of the sample regions (Figure 4a) over the alpine meadow (P1) and alpine desert areas (P2), where a negative value indicates ground subsidence. Both of the sample regions are characterized by seasonal displacement. As shown in Figure 6, an obvious difference was observed between the time-series subsidence amplitudes of the alpine meadow and alpine desert areas. The alpine meadow regions had peak-to-peak seasonal displacement with an amplitude of 32 mm, while the displacement was relatively stable over P2 (Figure 6b) compared with that of P1 (Figure 6a), with a peak-to-peak seasonal displacement amplitude of 20 mm (Figure 6a and b, the box plots of median values (red lines)). Similar to the sinusoidal seasonal model results, a larger spatial variation was observed for the alpine meadow areas, whereas larger variation ranges (whiskers of box plot for P2) and more outliers (red points in box plot for P2) occurred over certain areas of the desert, especially during the winter season due to the dry conditions and deep ALT. As shown in Figure 6, the permafrost table was stable at the end of winter and obviously subsided from 28 March 2018 to 24 September 2018. Rapid freezing uplift occurred over the alpine meadow areas in the initial stage of winter from October to December with an uplift of 20 mm and subsequently slowed, while the freezing uplift trends were stable throughout the whole winter over the alpine desert areas. These effects occurred because the time series permafrost displacement was controlled by the thermal state of the soil at different depths, while a number of environmental factors, such as peat layers and water availability, affected the thermal conductivity of the soil [65]. At the surface and subsurface of the alpine meadow areas, different types of soil and soil water content directly influenced the soil thermal conductivity, which led to the slowing of the uplift. In the alpine desert regions, the differences in the types of soil and soil water content were relatively small compared with those of the alpine meadow areas [8], and as a result, the difference in the time series deformation was insignificant over the whole winter period. Subsidence oscillations were present during the thawing period due to precipitation and variation of the soil moisture [37], [39].

B. ALT ESTIMATION RESULTS
In this study, two types of multilayer soil moisture data (ERA-Interim and SMAP L4) were used to establish the spatial-temporal distribution of the ground soil water content, and then, the seasonal subsidence during the melting period was converted to the ALT based on Equation (9). Figure 7 shows the estimated ALT based on the ERA-Interim (Figure 7a) and SMAP L4 (Figure 7b). The results show that the estimated ALT values are lower than those based on SAMP L4 data. It is because that the ERA-Interim data had higher soil moisture values than the SMAP L4 data ( Figure 8). As shown from the in situ soil moisture measurements from 2015 to 2018 (Table 1), the soil moisture VOLUME 8, 2020  ranged from 0.06 to 0.24 m 3 /m 3 . Comparing the in situ soil moisture to the SMAP L4 and ERA-Interim data, SMAP L4 had higher consistency, and the ERA-Interim data was overestimated (Figure 8), which means that the ALT results obtained based on ERA-Interim data were underestimated. The spatial distribution of the estimated ALT (Figure 7) was similar to the seasonal subsidence (Figure 4a and d). The ALT results showed strong spatial variability over the study area.
To assess the spatial pattern of the estimated ALT over different ground features in detail, the ALT was extracted along profile lines P1-P2 and separated based on the classification image (Figure 9a and b). The classification image was acquired using the S-1 amplitude map. Figure 9a shows the classification results using the support vector machine (SVM) method [66], and fine correlation was achieved compared with high resolution TerraSAR-X image based classification results [39]. To validate the classification results of Sentinel-1 amplitude image quantificationally, the Sentinel-2 optical image with a typical sample region was classified and treated as reference results (Figure 10a and b). The accuracy statistics were calculated with the confusion matrices from the Sentinel-2 classification results ( Table 3 ). The overall accuracy (OA), Kappa, producer's accuracy (PA), and user's accuracy (UA) are commonly used to assess the classification accuracies. The statistics results show that the overall accuracy was 0.91 (Table 3 ), as shown from the OA, PA, and Kappa in Table 3, higher classification accuracy in the alpine meadow area compared with the alpine desert area was achieved. The surface feature in the alpine desert area was more complex (Figure 10a), which led to relatively poor classification accuracy.
Obvious differences were observed between the alpine meadow and alpine desert areas in both the ERA-Interimbased ( Figure 9b) and SMAP L4-based (Figure 9c) estimation results. For the alpine desert regions, the mean ALT retrieved using the ERA-Interim and SMAP L4 data were 2.37 and 2.90 m, respectively, while the mean ALT retrieved were 1.19 and 1.62 m over the alpine meadow regions, respectively. There was a difference of approximately 0.5 between the ERA-Interim and SMAP L4-based results. This result was consistent with the relationship between the ALT and subsidence under different water content conditions built by Wang et al. [8]. The results indicated that changes in the soil moisture content greatly influenced the determination of the ALT. Large variation ranges in the ALT were found over profile P1-P2, especially in the alpine desert areas. Comparing the Figure 9b with Figure 9c, similar spatial variations were acquired from P1-P2. This result indicated that the spatial distribution of the ALT estimated was mainly controlled by the subsidence results (Figure 4a), and the range of the ALT was controlled by the soil moisture.

A. EVALUATION OF ALT RESULTS
To validate the precision of the ALT retrieved results, both regional scale estimation of the ALT based on the temperature-depth profile and in situ ALT data were introduced. At the global scale, the permafrost thickness can be obtained through linear interpolation from temperaturedepth profiles [3]. As shown in Figure 11, the melting time delay was determined in the ERA-Interim multilayer temperature data, and then, the melting process could be acquired. First, the time lags (e.g., T2-T1 in Figure 11) between two neighboring layers (e.g., Layer 1 and Layer 2) were calculated. Then, the melting rates of the different layer intervals were acquired using the time lags and design depth of the ERA-Interim data. Finally, the ALT was computed with the melting rate and melting period. The estimated ALT for the ERA-Interim data scale was 2.45 m. This result was similar to the total mean ALT of the SMAP L4 data-based results.   At the point scale, the in situ ALT data (A1-A4 in Figure 2b, Figure 12) were used to evaluate the reliability of the estimated ALT results. Especially, the in situ ALT of A4 (Figure 2b) including 10 sampling points with a land cover of Swamp meadow (SM), undisturbed Alpine meadow (AM), Degradation alpine meadow (DM) and Desert grassland (DG), the definitions of those samples identifications (ID) were SM1-SM3, AM1-AM3, DM1-DM2 and DG1-DG2, respectively. The detailed location of the specific sample can be seen in [45] and Figure 2b.
The ALT increased with notable trends. The in situ ALT acquired in the previous study was converted to the dates of our study by increasing the values at rates of 5.3 cm/yr and  4.2 cm/yr over the alpine desert and alpine meadow areas, respectively [45]. Figure 13b shows that the ALT estimation based on the SMAP L4 soil moisture data had reliable relative accuracy, with a correlation coefficient of 0.67. As shown in the scatterplot of Figure 13b and the linear polynomial fit expression, the ALT was underestimated. There was a relatively high absolute difference compared with the in situ results, with an RMSE of 0.7 m. Apart from the estimated error, this difference occurred because the in situ measurement time was inconsistent with period studied, and the actual change in the ALT may not strictly satisfy the rates of increase used. In addition, as shown in Figure 12, the ALT showed significant spatial variations even for the same type of land cover, which led to a higher RMSE between the estimation results and in situ measurements at the point scale. To get a more reliable estimation of the RMSE, the bias can be removed by defining the ubRMSE; the characterized random error was 0.5 m over the study area. As shown in Table 4, the underestimation and ubRMSE were mainly caused by the alpine meadow areas. Considering the GPR data, the in situ measurements located in the alpine meadow and alpine desert areas (Figure 2b A1 and Figure 12) show ALT of 2.5 and 3 m, respectively. The estimated ALT based on the SMAP L4 data were 1.75 and 3.31 m, respectively ( Table 4). The ALT was underestimated over the alpine meadow areas in A4 (Table 4). These results were similar to the ALT estimation results based on the TerraSAR-X data of Wang et al. [8]. For comparison with the in situ GPR of Lu et al. [46] from 5 June 2015 over A3, the ALT of 8 June 2018 was computed. Fine accordance was acquired for the alpine meadow areas with a mean difference of 0.1 m (Table 4).
Certain issues may have led to the errors in the estimations. First, the soil moisture of the alpine meadow areas had more spatial variation than alpine desert areas, and the spatial-temporal soil moisture at the global scale was inadequate for the retrieval model (Equation (9)), although the soil moisture ratios between the alpine meadow and alpine desert areas were separated by the in situ soil moisture measurements (Table 1). Second, as shown in Figure 12, the ALT showed higher spatial variations over the alpine meadow areas, which led to more discordance with the estimated results. In addition, the seasonal deformation that was used in the determination of the ALT was underestimated due to the uniform starting date defined; however, the thawing dates were different between the alpine meadow and alpine desert areas [8]. Additionally, the ALT was affected by the lithology and topography aspects. The solar radiation intensity differed for different aspects and may have led to large differences in the deformation [29]. The modified Stefan subsidence model was built according to the air temperature data, thus leading to subsidence estimation errors. As discussed in our previous study [39], climate anomalies such as rainfall were key factors impacted the phase measured and thereby masked the actual surface subsidence [67].
In addition, to evaluate the ALT estimation performance of our proposed method more comprehensively, a typical ALT estimation method proposed by wang et al. [8], was also performed to acquire the ALT of the study area. The retrieval precision was validated based on in situ ALT (Figure 13b and Table 4). Figure 13 shows that the ALT estimation based on point scale soil moisture data [8] was higher underestimated compared with our proposed method, with a correlation coefficient of 0.67, RMSE of 0.92 and ubRMSE of 0.54. As also shown in Table 4, compared with the in situ results, the absolute difference error of the traditional method is relatively high. In some samples the absolute difference error greater than 1m, such as desert grassland area (DG1 and DG2 in Table 4). The main reason is that only two surface classification categories were defined in this study, while, the transition region between Meadow and desert exist in the study area, the misclassification of those regions leads to the ALT estimation error. In the following study, more categories of the land surface should be defined and classified.

B. EVALUATION OF ALT RESULTSADVANTAGES AND LIMITATIONS
Previous studies have improved the InSAR-based ALT retrieval model using the vertical distribution of the soil water content [8], [44]. In these studies, the ground water distribution for different types of land cover was modeled. However, these studies were based on the assumption that the active layer was close to saturation [36], [41], [44], or only single soil moistures were input for different land covers for the surface soil moisture model [8]. As mentioned in [36], remotely sensed soil moisture data such as the soil moisture and ocean salinity (SMOS) and SMAP can be directly used in the ALT retrieval algorithm. To bypass these limitations, both the ERA-Interim and SMAP temporal-spatial multilayer data were introduced into the retrieval model (Equation (9)). First, VOLUME 8, 2020 the multidimensional soil moisture with temporal, spatial and depth information was considered. Second, the soil moisture proportions for different types of land cover were separated based on in situ measurements. In terms of the SMAP L4 soil moisture data, only the surface and root zone were provided for the estimation. Wu et al. [22] investigated the change in the ALT from 2002 to 2012 and concluded that the ALT was most closely related to the soil temperature at a depth of 0.5 m. Therefore, the SMAP L4 soil moisture data of the surface and root zone layers were reasonable for estimating the ALT.
However, there are still certain limitations to be considered that could be improved in future work. 1) The ERA-Interim soil moisture data provided more layers than the SMAP L4 data; however, these data were overestimated compared with in situ soil moisture observations ( Figure 8 and Table 1) and should be calibrated according to the in situ borehole data.
2) The precision of the global scale soil moisture data was not validated; when validation analysis of these data was performed using ground-based measurements, the optimal soil moisture data could be selected. 3) Both the ERA-Interim and SMAP L4 data are at a global scale; although improvements were acquired compared with the single point selections of previous studies, it is not enough to separate the soil moisture proportions based on the in situ measurements. In a future study, certain downscaling methods using the backscatter coefficient of multitemporal single or dual-polarized data could be developed to downscale the soil moisture data to a fine resolution [68].

VI. CONCLUSION
In this paper, we applied multilayer soil moisture distribution information and the MT-InSAR technique using S-1 data to estimate the ALT over a permafrost region in the QTP.
To obtain the groundwater distribution information in detail, the multilayer temporal-spatial soil moisture data were used, and thus, the multidimensional distribution of the groundwater in study area was acquired. Based on the experimental results obtained, the following conclusions can be drawn.
First, the seasonal subsidence in the area studied ranged from −60 to 0 mm. Higher slope terrain property responses were found in the modified Stefan model-based seasonal subsidence estimation results than in the sinusoidal model-based results, which demonstrated the applicability of the modified Stefan model to subsidence modeling. Spatial variations in the subsidence of different ecosystems were discovered in this study. The freezing induced uplift over alpine meadow areas was rapid in the initial stage of winter from October to December and subsequently slowed.
In addition, the experimental results showed that the retrieved mean ALT was 1.62 m in the alpine meadow areas and 2.9 m in the alpine desert areas. A relatively fine consistency was shown in the SMAP L4 soil moisture data-based results compared with that of field investigations. The ALT was underestimated in the ERA-Interim soil moisture-based results due to the overestimation of the data.
Finally, the ALT estimation based on MT-InSAR and temporal-spatial soil moisture data has higher performance than traditional methods. Thus, the surface and root zone soil moisture information characterized by large scale temporalspatial soil moisture product is a promising method that can be used for the ALT monitoring over permafrost regions.
In future work, ERA-Interim soil moisture data calibration and global scale remote sensing soil moisture data downscaling should be applied to improve the precision of the soil moisture information.  ZHENGJIA ZHANG received the B.Sc. degree in geoscience information system from the Hefei University of Technology, Hefei, China, in 2012, and the Ph.D. degree in microwave remote sensing from the Institute of Remote Sensing and Digital Earth, Chinese Academy of Sciences, Beijing, China, in 2017. He is currently a Lecturer with the Faculty of Information Engineering, China University of Geosciences, Wuhan, China. His research interests include synthetic aperture radar (SAR) interferometric technique, time-series interferometric SAR, and their applications on nature hazards and permafrost regions. VOLUME 8, 2020