1 Introduction

Permafrost change has received much attention in recent years due to its sensitivity to global warming (Liu and Chen, 2000). Increased temperature can cause phase changes between soil ice and water in the active layer, which has strong influences on land-surface processes in cold regions (e.g., thermal conditions; biological, hydrological, and pedological activities) (Hinzman et al., 1991 ; Kane et al., 1991 ; Cheng and Wu, 2007; Qin et al., 2016 ). By contrast, climate warming also has great negative impact on cold-region infrastructures such as the Qinghai–Tibet highway and railway. Therefore, it is necessary to investigate seasonal and long-term surface changes in these permafrost regions.

A variety of traditional methods has been used to figure out the seasonal or long term changes of the active layer, for instance, mechanical probing (Brown et al., 2000 ; Cheng and Wu, 2007; Zorigt et al., 2016 ) and frost or thaw tubes (Harris, 1970; Warren and Brown, 1972; Mackay, 1973a; Mark Nixon and Taylor, 1998). Data acquired from temperature sensors can also be used to assess differences in thaw depth during freezing and thawing periods (Zhang et al., 1997 ; Brown et al., 2000 ; Wu et al., 2010 ). With the development of the technology in geophysical survey in the past decades, time-domain reflectometry (TDR) and GPS were also applied to determine the variation of the active layer (Baker et al., 1982 ; Little et al., 2003 ). However, some of these traditional measurements, such as drilling boreholes in wintertime, are not easy to obtain due to the tough fieldwork conditions on the Tibetan Plateau, with high altitude and cold weather. Furthermore, the greatest disadvantage of all these in-situ methods is that only sparse, point-based observations can be supplied, though these observational results have high accuracy. Therefore, mapping permafrost displacement over a broad area by D-InSAR should be an adaptive way for determining seasonal and long-term surface changes of permafrost.

Advanced InSAR methods such as Permanent Scatters InSAR (PS-InSAR) and SBAS-InSAR have been proposed for detecting surface change using D-InSAR technology (Ferretti et al., 2000 ; Usai, 2001; Berardino et al., 2002 ). SBAS-InSAR has already been applied to retrieve long-term surface-deformation data over cold regions or mine areas, where permanent scatters are difficult to find (Castañeda et al., 2009 ; Liu et al., 2010 ; Liu et al., 2012 ; Chen et al., 2013 ; Short et al., 2014 ). Furthermore, SBAS-InSAR can significantly reduce the de-correlation induced by a long spatial and temporal baseline because the least-square method is incorporated into this approach for improving temporal resolution (Berardino et al., 2002 ). These specific advantages show a potential application of the SBAS-InSAR method for deformation measurement over permafrost regions. Previous studies of deformation measurements with D-InSAR technology have been conducted in some cold regions (Li et al., 2004 ; Liu et al., 2010 , 2012; Daout et al., 2017 ), and C-band and L-band data also have been utilized for detecting the surface changes of permafrost over the Tibetan Plateau (Xie et al., 2010 ; Chen et al., 2013 ; Daout et al., 2017 ). However, baseline or temporal de-correlations still have great adverse impacts on D-InSAR results, due to large baselines between acquisitions and vegetation growth or snowfall; and radar results also are significantly affected by phase delay, which is in turn caused by atmospheric factors such as humidity (Zebker et al., 1997 ). Although some previous studies have inverted seasonal and long-term surface variation of the active layer on the Tibetan Plateau, some of these studies still lack the ground truth. In addition, few studies use Sentinel-1 data for inverting time-series displacements of ground surfaces in permafrost regions of the Tibetan Plateau and validating its applicability.

In this study, using the SBAS-InSAR method, we acquired and analyzed the seasonal and long-term deformations of the TGL and LDH regions during two freeze thaw cycles, from October 2014 to June 2016, with Sentinel-1 images. The objective of this study was to estimate the applicability of Sentinel-1 C-band SAR data over specific permafrost regions under different vegetation coverage and soil-moisture regimes by using in-situ leveling observations.

2 Study area

TGL (91.8°E–92.0°E, 33.0°N–33.2°N) is located in the interior of a continuous permafrost zone on the Tibetan Plateau (Figure 1), with a high altitude of around 5,100 m a.s.l.. This region is covered by degraded alpine meadow, with plant stature less than 10 cm and coverage of 20%–30%. The mean annual air temperature is approximately −4.9 °C, and the extreme maximum and the minimum values of daily average temperature are 17.6 °C and −29.6 °C, respectively (Jiao et al., 2014 ). Precipitation primarily occurs in summer, between May and September, with an annual precipitation of about 400 mm; and the maximum precipitation occurs in July and August. The active-layer thickness of the in-situ leveling observational field in the TGL region is about 3.15 m. The soils are composed of gravelly sandy loam, including around 5%–30% gravel, 65%–95% sand, and 3%–22% silt, as well as 2%–13% clay (Hu et al., 2015 ).

LDH (91.5°E–91.85°E, 31.6°N–32.1°N) is located on the south margin of a continuous permafrost zone of the Tibetan Plateau (Figure 1). Vegetation coverage, soil moisture, and topography in this region significantly differ from those in the TGL region. Vegetation type at LDH is alpine paludal meadow, with a vegetation coverage of more than 80%. The average altitude of this area is above 4,500 m a.s.l., the mean annual air temperature is around −1.3 °C, and the mean annual precipitation (from 1971 to 2004) is about 411 mm (Yang et al., 2012 ). The active-layer thickness is around 1 m. Soil in this area is much wetter than that at TGL, and a pure-ice layer exists at the bottom of the active layer.

3 Data sources 3.1 Sentinel-1 C-band data and DEM

The Sentinel-1 program was established by the European Space Agency to replace the older Earth-observation mission of the European remote-sensing satellite (ERS) and ensure continuous data acquisition. The Sentinel-1 mission includes two polar-orbiting satellites with an altitude of 693 km, which means they can comprise a constellation and improve the observation efficiency. These two satellites both carry a C-band synthetic aperture radar (SAR) device that can operate day and night, and it also can acquire all-weather images (European Space Agency, 2013). Two sets of 27 and 42 Sentine-1A/B data acquired from October 2014 to June 2017 and covering TGL and LDH were utilized in this study; see the list in Tables 1 and 2, respectively. All acquisition images used in this study were collected in Interferometric Wide Swath (IW) mode, with a swath width of around 250 km and spatial resolution of 5m×20m.

Precise orbit files of remotely sensed data, which are necessary in interferometry processing were obtained from the Payload Data Ground Segment (PDGS). The Shuttle Radar Topography Mission DEM (SRTM DEM), with a spatial resolution of 30 m and absolute height error of less than 16 m (Farr et al., 2007 ), was also used in this study because the indispensable DEM is used for topography subtraction during differential interferometry processes.

3.2 In-situ leveling data

A leveling observational field was arranged at TGL in 2010, consisting of one leveling datum point and 25 benchmarks evenly distributed over an area of 16 m². The leveling-measurement datum point is an iron tube, inserted into a 60-m-deep borehole, with about 15 cm protruding aboveground. The leveling observing site was also established at the LDH region in October 2015, including 16 benchmarks with an area of around 500 m² and one datum point (same as the datum point in TGL). We hypothesized that the leveling datum points of the two observation sites kept stable during each freeze–thaw cycle because the iron tubes were deeply inserted into the ground (Mackay, 1973b, 1977). These leveling marks were measured every month or two with a SOKKIA SDL30 digital level, which performs well in adverse environments, with a height-observing accuracy of 1.0 mm (standard deviation for 1-km double-run leveling) using random bidirectional (RAB) code staves.

4 Methods

The SBAS method assumes an unwrapping phase of each interferogram as observation and then acquires the cumulative phase by solving a linear least-squares problem.

Considering that N+1 images at the same orbit are collected in chronological sequence of t₀, t₁, …, t_n, and hypothesizing that each image could serve as a master image and produce at least one interferogram with other acquisitions, total M interferograms are produced during the interfering process; and they satisfy the following equation:

$\frac{{N + 1}}{2} \le M \le N\left( {\frac{{N + 1}}{2}} \right)$

(1)

Assuming that the j-th interferogram is calculated by acquisition at the time of t_A and t_B, so the azimuth and range coordinates of x and r can be described as

$\begin{align} \text{δ} {\emptyset _j}\left( {x,r} \right) =& {\emptyset _B}\left( {x,r} \right) - {\emptyset _A}\left( {x,r} \right) \\ \approx & \frac{{4{\text{π}}}}{\lambda }\left[ {d\left( {{t_B},x,r} \right) - d\left( {{t_A},x,r} \right)} \right] \\ &+ \Delta \emptyset _{{\text{topo}}}^j\left( {x,r} \right) + \Delta \emptyset _{\rm {APS}}^j\left( {{t_B},{t_A},x,r} \right) + \Delta \emptyset _{\rm {noise}}^j\left( {x,r} \right)\end{align} $

(2)

where λ is the transmitted-signal central wavelength and j∈(1, …, M); d(t_B, x, r) and d(t_A, x, r) mean the cumulative deformation in the direction of line of sight (LOS) at times t_A and t_B, relating to the value of d(t₀, x, r)=0.

$\Delta \emptyset _{{\text{topo}}}^j$ (x, r) represents the topography phase, most of which can be removed by an available external DEM; therefore, the residual topography phase can be neglected during the interfering process because of its small magnitude. $ \Delta \emptyset _{{\text{APS}}}^j$ (x, r) means the phase component caused by atmospheric disturbances. $\Delta \emptyset _{{\text{noise}}}^j $ (x, r) indicates de-correlation noise, such as system noise (caused by system thermal and quantisation). Equation (2) can be simplified as below without taking $\Delta \emptyset _{\rm {APS}}^j $ (x, r) and $\Delta \emptyset _{\rm {noise}}^j$ (x, r) into consideration:

$\begin{align} \text{δ} {\emptyset _j}\left( {x,r} \right) &= {\emptyset _B}\left( {x,r} \right) - {\emptyset _A}\left( {x,r} \right) \hfill \\ & \approx \frac{{4 \text{π}}}{\lambda }\left[ {d\left( {{t_B},x,r} \right) - d\left( {{t_A},x,r} \right)} \right]\end{align} $

(3)

For the purpose of obtaining the cumulative deformation with physical meaning, the phase in Equation (3) can be expressed as the mean phase velocity multiplied by the acquisition-time difference:

${v_j} = \frac{{{\emptyset _j} - {\emptyset _{j - 1}}}}{{{t_j} - {t_{j - 1}}}}$

(4)

Therefore, the phase value of the j-th interferogram can be given by

$\mathop \sum \limits_{k = {{t_{A,j+ 1}} }}^{{t_{B,j}}} \left( {{t_k} - {t_{k - 1}}} \right){v_k} = \text{δ} {\emptyset _j}$

(5)

That is the integration of phase velocity in each period between master and slave images. With this equation organized in matrix form, it can eventually be expressed as:

${\rm{B\nu }} = {\text{δ}}\emptyset$

(6)

where B is always a rank-deficiency matrix with a size of M×N because many different images are chosen as the master image during the computing process. By contrast, a generalized matrix of B can be obtained after applying the singular-value decomposition (SVD) method, and then the velocity (v) can be calculated by employing the minimum-norm solution. Finally, deformation of various periods can be achieved according to the integration of phase velocity at the same time.

In this study, InSAR Scientific Computing Environment (ISCE) software was implemented in D-InSAR processing because it is freely attainable and widely used in the InSAR scientific community (Rosen et al., 2012 ). An adaptive filter (Goldstein-Werner filter) was applied in this software to smooth interferograms (Goldstein and Werner, 1998), and statistical-cost, network-flow algorithm for phase unwrapping (SNAPHU) was employed in phase unwrapping in the final stage of interfering. Genetic InSAR Analysis Toolbox (GIAnT) software was applied in the time-series analyzing and deformation computing process (Agram et al., 2012 , 2013). A high-coherence area in most interfering results was assumed to be stale during the monitoring period, and this area was selected as the phase reference area during the time-series processing, which means that all inverted deformations referred to this area. Deformations of a specific pixel (from SBAS-InSAR time-series results) that corresponded with an in-situ measurement area was extracted and compared with ground truth, to validate the D-InSAR results. A flow chart of data processing is shown in Figure 2.

Before time-series InSAR analysis, interferograms should first be produced. In the TGL region, all total 103 interferograms (Figure 3) were produced when setting the spatial perpendicular baseline less than 100 m and the interval less than 200 days between master and slave acquisitions. All the imagery data was cropped as 0.2°×0.2° grids to save computing time. A total of 324 interferograms (Figure 4) were produced in the LDH region. The largest perpendicular baseline and time separation were set as 100 m and 300 days, respectively. Similarly with the TGL region, only a small region in LDH, around 0.50°×0.35° in terms of latitude and longitude, was calculated in the interfering process.

During SBAS processing, the threshold of coherence and minimum number of interferograms were set as 0.23 and four separately in TGL, which means that pixels would be included in the time-series processing with a coherence higher than 0.23 in more than four interferograms. In the LDH region, pixels' coherence were to be no lower than 0.2 in no fewer than 50 interferograms during the time-series calculating process. The lower coherence threshold was employed for avoiding seriously incoherent results due to LDH's high vegetation coverage and abundant soil-water content; more interfering results restraint was used because a large number of interferograms were produced during D-InSAR process.

5 Results

In Figure 5, we compared the results observed from the SBAS-InSAR method with terrestrial leveling measurements at TGL. The deformation obtained from the two methods showed a consistent trend and significant seasonal variations during the study period. As observed by D-InSAR and terrestrial leveling, subsidence of around 5 mm occurred from late October to December 2014, though the soils should have been uplifted due to the effects of frost heave during that period. Similar phenomenon can also be identified during the same period of 2015 and 2016, which means that land subsidence did occur during the autumn freezing stage in the study area.

As InSAR results have shown, the land surface in TGL rose by about 13 mm from mid-December 2014 to mid-July 2015; and it increased by 4 mm, with remarkable fluctuations from late-January to August 2016 (Figure 5). Meanwhile, the active-layer thickness increased 6 mm, according to the leveling measurement. The terrestrial leveling observation indicated no subsidence during summer 2016. On the contrary, the ground rose around 15 mm by October 2016, as compared with June 2016. The in-situ observations are sparser than the satellite acquisitions in temporal distribution. Displacement results of InSAR measurements were larger than the leveling observations before July 2016, but these values were less than the ground truth after that time.

Figure 6 demonstrates the cumulative deformation of TGL during the two-year measurement period, i.e., the total ground displacement on November 3, 2016 (Figure 6b) relative to the surface height on October 27, 2014 (Figure 6a). The study area of TGL is included in the black rectangle, and deformation of these two maps is in the satellite's LOS direction, as measured in mm. A positive value means upwards from the land surface moving towards the satellite; negative represents land surface moving away from the satellite (subsidence). Some pixels appear as white in the figures because the coherence of these areas in the interferograms is lower than the coherence threshold. Therefore, these pixels were not included in calculating the deformation during the SBAS process. Regions that surround the study area appear as either subsidence or uplift, which means ground-deformation change is heterogeneous on the Tibetan Plateau. As evident from Figure 6, the region of interest underwent a small subsidence of around 5 mm during the monitoring period from late October 2014 to early November 2016.

As shown in Figure 7, the trend and magnitude of deformation obtained from D-InSAR in the LDH region were inconsistent with in-situ observations. The leveling results were much higher than the measurements of radar technology, indicating an insensitivity of the deformation measured by Sentinel-1 C-band data to the freeze–thaw cycles of active-layer soils, indicating the Sentinel-1 data set performs poorly in the LDH region.

6 Discussion

The deformation process over cold regions is relatively intricate. On a local scale, permafrost displacement is determined by many factors, such as water/ice content of the active layer, its thickness, vegetation coverage, ground-surface roughness, and geomorphological processes, among many others (Liu et al., 2010 ; Chen et al., 2013 ). In permafrost regions, seasonal deformation is mainly affected by freeze thaw processes of the active layer, due to water–ice phase change, along with temporal variation, i.e., ground-surface heaving in cold seasons and subsiding in warm seasons. In addition, long-term surface variation is caused by ice–water phase changes in soils, which might be a useful proxy for permafrost agradation or degradation.

6.1 Seasonal and long-term changes of ground-surface deformations in the TGL region

It is difficult to analyze the temporal variations in ground-surface deformation in detail due to the sparse in-situ measurements during the monitoring period; and thus the displacement observed by D-InSAR was analyzed in the following sections. Generally, the freeze–thaw processes over permafrost regions on the Tibetan Plateau can be divided into four stages: thawing settlement in summer, frost heave in autumn, cooling process in winter, and warming process in spring (Zhao et al., 2000 ). However, we found some anomalous deformations inconsistent with those stages were presented during the monitoring period (Figure 6). For example, the uplift of land surface occurred from May to mid-July 2015 and between June and October 2016, which can be observed from InSAR and the terrestrial leveling method, respectively. These anomalies can be explained by the rainfall or meltwater infiltration in summer. Seasonal ground thaw commonly occurs from early May to late October each year (lasting around 110 days), and the greatest concentration of precipitation falls at this region during the same period, accounting for about 83% of the annual precipitation (Jiao et al., 2014 ). The soil type of the TGL region is gravelly sandy loam with large porosity, and thus the percolating rainwater migrates downwards to the freezing active-layer soils beneath the thawing surfaces during the thawing period, subsequently increasing the ice content of the frozen active layer and hence the frost heaving of the ground surface. These phenomena have been verified in laboratory and field experiments at various permafrost regions, such as on the Tibetan Plateau and along the western Arctic coast (Cheng, 1982; Mackay, 1983). Water infiltration during the thawing period may have great impact on hydrological processes and on latent and sensible changes over permafrost regions. Therefore, these implications should be considered in the studies on hydrological and energy exchanges, by relating these processes with land-surface displacement.

In the TGL region, simple linear regressions were applied on InSAR-derived and ground-based measurement results for detecting long-term variation in the active-layer thickness. Both the InSAR results and leveling measurements showed a decreasing trend during the study period, with a larger slope of D-InSAR (Figure 5). The main reason for this trend should be the low sampling frequency of the terrestrial leveling observations, which may cause the loss of some larger values of deformation during the observing period. For instance, surface height was not recorded between 18 April 2015 and 30 October 2015 or between 31 October 2016 or 8 March 2017 by leveling measurement. Another possible explanation is that D-InSAR results included some residual phase during the interfering process due to DEM errors or atmospheric delay, and this question should be solved in future studies.

Land-surface displacement over cold regions might be used to indicate permafrost changes, but it is difficult to monitor slight changes of active-layer thickness. For example, it is impossible to detect a surface settlement of 5 mm using ground temperature or drilling methods. In the TGL region, the decreased displacement trend might signify permafrost degradation because the melting of ice in the top of permafrost can induce surface settlement. Some studies have shown that significant variation of active-layer thickness occurred in the past three decades in the eastern and inland Tibetan Plateau (including the TGL region) because of remarkable climate warming (Zhao et al., 2004 ). A similar trend of permafrost degradation was also obtained on the Alaska North Slope between 1992 and 2000, with a magnitude around 10–40 mm per decade (Liu et al., 2012 ), while the subsidence of about 1–3 mm per year occurred, according to terrestrial leveling measurement and D-InSAR results in the TGL region, separately. This phenomenon demonstrates a permafrost degradation in these regions under climate warming.

6.2 Applicability estimation of D-InSAR method over the TGL and LDH regions

Although it is difficult to estimate the accuracy of D-InSAR observation with statistical methods due to limited ground measurements, the differences between D-InSAR-derived displacements and leveling observations at TGL were less than 10 mm. However, irregular and incorrect deformations were obtained by the D-InSAR method over LDH. Such a large difference between TGL and LDH might be attributed to the differences in vegetation coverage, soil-water content, or a combination of the two.

Table 3 shows the porosity, vegetation coverage, and water content of all layers sampled at the two study sites. Details from this table illustrate that water content and porosity of soil at LDH is larger than at TGL in various sampled layers. We have discussed the frost heave in TGL during the ground-thawing period due to large soil porosity and concentrated precipitation in summer. However, the same phenomenon did not display in LDH according to in-situ measurements, though porosity in LDH is larger than that in TGL; and precipitation there concentrated at the same period. Maybe normal deformation was related to exorbitant soil-water content in the LDH region during the thawing period, and longer leveling measuring periods would be necessary to validate this hypothesis.

Soil-water content and vegetation coverage may be the control factors for deformation-estimating in LDH. Soil with higher moisture content may result in specular scattering and coherence loss, which leads to incorrect deformation records because the phase values might be smoothed by the adaptive phase filter, and the unwarping phase may be produced incorrectly in these areas (Short et al., 2014 ). By contrast, soil-water content changed substantially from April to July (Table 3) in the LDH region. This change could cause complex phase changes, which has been demonstrated in laboratory and field experiments (Castañeda et al., 2009 ; Morrison et al., 2011 ). Another reason for the difference is high vegetation coverage in the LDH region, with some research indicating that vegetational variations show a significant consistency with the phase diversity; and it could adversely influence the final observational results of the ground-surface deformations (Zwieback and Hajnsek, 2016).

7 Conclusions

The D-InSAR method can supply estimates of seasonal and long-term displacements of land surface over cold regions with broader coverage, and this approach has even higher temporal resolution than our ground-based leveling method. Compared with in-situ leveling measurements, Sentinel-1 C-band data performed well in the sparsely vegetated and drier TGL region. Some detailed variation, such as frost heave in summer time, can be identified by the D-InSAR method; and a long-term surface-subsidence tendency also can be detected by this approach. Seasonal change of land surface at the TGL site was less than 20 mm, which was mainly caused by the freeze–thaw cycle of active-layer soils; and long-term subsidence was around 1 to 3 mm per year, due to the thawing of ground ice. However, Sentinel-1 data showed poor results in the LDH region; and deformation was greatly undervalued by the D-InSAR method, mainly due to higher vegetation coverage and greater soil-water content.

Ground-based results are obtained by measuring limited points on the ground, but InSAR-observed deformations displayed changes over a broad area. Therefore, it is impossible to analyze accuracy between InSAR and leveling measurements based on statistical methods. In the future, some corner reflectors should be installed to enable acquiring more detailed comparative information between InSAR and ground-based measurements. Radar data with different bands also should be tested in the LDH region for the collection of surface deformations. For example, L-band data may produce better results in this region.

Acknowledgments:

This work was supported by the Innovation Groups of the National Natural Science Foundation of China (41421061), the Chinese Academy of Sciences (KJZD-EW-G03-02), the project of the State Key Laboratory of Cryosphere Science (SKLCS-ZZ-2017), and CUHK Direct Grant (4053206). The Sentinel-1 SAR data were provided by the European Space Agency (ESA) through Sentinels Scientific Data Hub. The authors would like to thank the two anonymous reviewers, as well as the editor for their helpful comments and suggestions.