Preliminary Study on InSAR-Based Uplift or Subsidence Monitoring and Stability Evaluation of Ground Surface in the Permafrost Zone of the Qinghai–Tibet Engineering Corridor, China

Abstract

Against the background of global warming, permafrost areas are facing increasing thawing, and the threat to the surface of the Qinghai–Tibet Engineering Corridor (QTEC) is serious. It is imperative to understand the current surface deformation and analyze the changes spatiotemporal characteristics for future warnings. At present, observation of a long time series and overall coverage of vertical ground deformation in QTEC are lacking. This paper takes the permafrost deformation of the QTEC as its research object. It uses the pretreated LiCSAR product and combines it with the LiCSBAS package to obtain monitoring results of the long time series deformation of the engineering corridor’s surface. The SAR image acquisition date is taken as the constraint, the results covering the whole processing area are selected, and then the vertical deformation information covering the entire engineering corridor area by ignoring the north–south displacement is calculated. The results show that the surface of the study area, as a whole, slightly subsided between May 2017 and March 2022, and the vertical deformation rate was mostly distributed at −27.068 mm/yr − 18.586 mm/yr, with an average of −1.06 mm/yr. Vertical deformation dominated at 52.84 percent of the study area, of which settlement accounted for 27.57 percent and uplift accounted for 25.27 percent. According to the statistics of the normal distribution of deformation velocity per pixel, a total of 77% of the engineering corridor was stable, with a vertical deformation rate between −6.964 mm/yr and −4.844 mm/yr, and 17.7% of the region was sub-stable, with a settling rate of −12.868 mm/yr − –6.964 mm/yr. The unstable regions included areas with settlement rates greater than 12.868 mm/yr and uplift rates greater than 10.748 mm/yr, representing 4.4 percent and 0.9 percent of the total area, respectively, for a total of 5.3 percent. The results of this paper can be used as the theoretical basis and as basic data for decision making and scientific research in various departments, and they are of great significance for surface stability assessment and early warnings along engineering corridors and traffic projects.

1. Introduction

Frozen ground, also known as frozen soil, is a geological formation where the temperature remains at or below 0 °C in winter, while permafrost is a type of rock formation where the temperature remains at or below 0 °C for at least two consecutive years [1]. Permafrost is a significant component of the Earth’s surface, covering approximately 20% of the total land area, and it is found on all continents, except Australia [2]. In China, seasonal permafrost (excluding transient permafrost) is widely distributed over an area of approximately 5.36 × 106 km², and its permafrost is the third largest in the world, covering an area of approximately 1.59 × 106 km², excluding glaciers and lakes [3]. Approximately 80% of the permafrost in China is located on the Qinghai–Tibet Plateau (QTP) [4]. As global warming continues to accelerate, permafrost areas are facing increasingly severe thawing problems, which have resulted in a series of disasters and significant economic losses [5]. These issues have attracted significant attention from scholars and policymakers around the world. Studying the dynamics of permafrost thawing and its impacts on the environment and society is crucial for developing effective strategies to mitigate these problems.

The Qinghai–Tibet Engineering Corridor (QTEC) is situated in the QTP hinterland, where unprecedented threats to its surface have emerged due to the dual influences of climate warming and human activities. The area is characterized by strong and persistent permafrost degradation, increased active layer thickness, and severe surface subsidence, which pose significant challenges to the normal operation of linear infrastructure projects, such as the Qinghai–Tibet Highway (QTH), Qinghai–Tibet Railway (QTR), DC (direct current) transmission towers, and the Golmud–Lhasa pipeline [6,7]. Therefore, comprehending the current surface deformation and analyzing the temporal and spatial changes is essential for future early warnings. Traditional field monitoring methods rely on single-point, linear, or small-scale areas, such as level surveys, buried instruments (positive inverse plumbs, multipoint displacement meters, etc.), differential GPS (global positioning systems), borehole and level surveys, pit explorations, physical explorations (ground-penetrating radar and high-density electrical methods), optical fibers, ground 3D laser scanning, and unmanned aerial surveys [8,9,10,11]. However, these methods have limited temporal and spatial resolutions, and they require significant human and material resources to carry out site monitoring in the field [12,13]. Therefore, monitoring methods or monitoring results with a large scope, long time series, high accuracy, low consumption, and high spatial and temporal resolutions are highly sought after by operating companies and decision makers.

Synthetic aperture radar (SAR) is an active microwave remote sensing technique that possesses the unique characteristics of being able to operate in all weather conditions and independent of cloud cover. As a result, it has found wide-ranging applications in various fields, such as digital elevation model (DEM) generation [12,14], surface deformation measurement [15,16], urban deformation monitoring [17,18], infrastructure deformation monitoring [19,20], mine deformation monitoring [21,22,23], seismic deformation monitoring [24,25], volcanic deformation monitoring [26,27], landslide monitoring [28,29], glacier movement monitoring [30,31,32], permafrost process monitoring [33,34,35,36,37], and many others. However, traditional InSAR (interferometric synthetic aperture radar) processing methods often rely on single ascending or descending orbit data, which can only provide deformation results along the line-of-sight (LOS) based on the relative position of the satellite. Such methods have limitations, as they are highly dependent on the SAR satellite’s geometric parameters and may not accurately reflect the actual deformation of the ground surface. Moreover, most processing methods start with single look complex (SLC) SAR images, which consume significant hard disk and computer computational resources and have low computational efficiency. Additionally, incorporating new data becomes difficult.

In recent years, advancements in computer technology, big data, and cloud computing have led to the development and application of cloud platforms in the field of traditional optical remote sensing. Platforms such as the Google Earth Engine (GEE) (https://developers.google.cn/earth-engine/, accessed on 29 March 2023), Pixel Information Expert Engine (PIE-Engine) (https://engine.piesat.cn/, accessed on 29 March 2023), and AI Earth of Alibaba Damo Academy (https://engine-aiearth.aliyun.com/#/, accessed on 29 March 2023) have emerged and been vigorously developed and applied [38,39,40]. By combining InSAR technology with cloud platforms, the pressure of traditional desktop-side data processing can be reduced, the entry platform for InSAR can be lowered, and the technology can be combined with practical applications to realize industrialization. One such example of a big data processing platform is LiCSAR (Looking Into Continents from Space with Synthetic Aperture Radar), which processes free and open-source Sentinel-1 raw data through supercomputing for large-scale interferometry to automatically generate products such as geocoded interferograms (filtered and unfiltered), coherence maps, and unwrapped phase images (after filtering) for large regions [41]. Another tool, LiCSBAS [42], takes the LiCSAR product as input and applies the new small baseline subsets (NSBAS) algorithm for timing analysis and checking loop closures before performing deformation timing inversion to improve accuracy [42,43,44]. In addition, LiCSBAS supports the use of the Generic Atmospheric Correction Online Service (GACOS) for InSAR for atmospheric phase removal [45,46]. Considerable results have been achieved based on the research and application of LiCSAR products and LiCSBAS tools [47,48,49,50,51], including work in permafrost regions [52]. By leveraging these technological advancements, the pressure of traditional desktop-side data processing can be reduced, and the efficiency and accuracy of InSAR processing can be improved, paving the way for future research and practical applications.

Despite the growing popularity of SAR in monitoring surface deformation in permafrost areas along the QTH and QTR (i.e., the QTEC), there exist certain limitations that hinder its effectiveness. These limitations pertain to the computing performance and theoretical knowledge of InSAR. As a result, most of the InSAR-based studies conducted in these areas are based on small areas and a single orbit (ascending or descending), with LOS satellite SAR being the primary focus. Research exploring the use of multi-orbital SAR data for capturing vertical deformation over large areas and long time series remains relatively unexplored in this field.

This paper aimed to address the current limitations of InSAR-based studies on surface deformation in permafrost areas of the QTEC. To achieve this, this paper utilized the preprocessed LiCSAR product and combined it with the LiCSBAS package to obtain the LOS deformation monitoring results of the engineering corridor surface over a long time series. The results were screened to cover the entire processing range and to neglect the north–south displacement for decomposing the vertical (UD, up and down) and east–west (EW) deformation results. The long time series UD deformation results covering the entire project corridor were analyzed using various methods, such as parametric statistics, Sen trend analysis (Theil–Sen median method), M–K (Mann–Kendall) significance test, and Hurst index. This analysis aimed to identify the differences between the deformation’s characteristics and development in order to achieve a comprehensive understanding of the surface deformation. Furthermore, this paper aimed to provide theoretical support and basic data for the decision making and scientific research of various departments by completing an overall understanding of surface deformation and achieving the goal of the evaluation and future early warnings of surface stability along engineering corridors and transportation projects.

This paper is organized as follows. Section 2 provides a detailed introduction to the geographical location and basic geographical features of the QTEC. Additionally, it includes information on the regional divisions of the processing area of the InSAR. Section 3 focuses on the methods and data used in this study, including the acquisition of the InSAR deformation results and post-processing of the data. In Section 4, the deformation results obtained in this paper are analyzed, and the stability of the ground surface is evaluated. The relevant results obtained by this research are thoroughly examined, including the accuracy of the results, their characterization, and their practical application. Based on these results, the methods used in this study, the analysis of the results, and the potential follow-up work are discussed in Section 5. Finally, Section 6 presents the conclusions.

2. Study Area Overview

The definition of the QTEC lacks clarity and is often understood as the area along the QTH or QTR within a certain range, or a combination of both within a specific boundary. The boundary distance is a subject of debate, with values such as 5 km, 10 km, and 20 km being commonly used [6,7,53]. This area is typically classified as a buffer zone. In this study, we used the 20 km buffer zones of the QTH and QTR (from Xidatan to Amdo) and combined the results of their respective 20 km buffer zones to obtain the engineering corridor boundaries, as illustrated in Figure 1. It is essential to note that during the creation of the buffer zones, the use of the projection coordinate system was necessary to ensure accuracy, as the buffer distance unit was in kilometers. In this paper, the WGS 1984 geographic coordinate system was utilized, with the Albers equal-area projection (orthogonal isoparametric conic projection) applied to ensure computational efficiency and data processing convenience. To further enhance the efficiency of the data processing, we divided the entire engineering corridor into four small areas (areas 1, 2, 3, and 4 in Figure 1b,c), with overlapping parts between the subareas. The mosaic ranges of the subareas encompassed the entire engineering corridor.

The geographical coordinates of the engineering corridor were 91°15′46.836″E–95°06′35.0532″E, 32°0′26.442″N–36°34′6.2076″N. The toponymic information involved in the scope of the engineering corridor was derived from the toponymic and notation datasets in the 1:1 million public version of Basic Geographic Information Vector Data (Version 2021) (https://www.webmap.cn/commres.do?method=result100W, accessed on 29 March 2023) released by the National Catalogue Service for Geographic Information (see Tables S2 and S3 for details).

The QTEC is a unique and challenging engineering project located in the QTP, with severe climatic conditions and high elevations ranging from 2725 m to 6175 m, as depicted in Figure 1b. The average elevation is 4618 m, with a standard deviation of 556 m. The terrain is characterized by alternating mountain systems (high plains) and river valley bottoms, with valleys interspersed between basins and mountain systems. The area is primarily situated in the hinterland of the QTP, with the exception of the southern part of the Tanggula Mountains, which is influenced by an oceanic climate, and the northern part, which experiences the arid climate of the Qaidam inland. The region experiences a dry and cold climate, with significant vertical zoning as the altitude increases [4]. The permafrost zone along the route spans three larger natural climatic zones, including the arid climatic zone north of the Kunlun Mountains, the plateau dry cold climatic zone between the Kunlun Mountains and Tanggula Mountains, and the plateau subarid climatic zone south of the Tanggula Mountains. These zones differ significantly in various climatic elements [54].

The permafrost zone in the QTEC exhibits a pronounced seasonal temperature variation, with the annual average temperature ranging from −6.9 °C to 7.4 °C. The lowest temperature, which usually occurs in January or December, ranges from −17.4 °C to −14.5 °C, while the highest temperature is recorded in July, averaging between 6.5 °C and 8.1 °C. The annual average temperature variation ranges from 15 °C to 26 °C, with the extreme variation not exceeding 50 °C. The average daily variation throughout the year ranges from 10 °C to 19 °C, with an extreme daily variation of 35 °C [55]. In terms of precipitation, the range is mainly between 200 mm and 500 mm, with a maximum average annual precipitation variation rate of 5.1 mm/yr, increasing with decreasing latitude. Overall, there is a trend of warming and increasing humidity within the engineering corridor [56].

The study area exhibits large spatial variation and vertical and latitudinal zonality in the distribution of permafrost, with a dry zonality spatial distribution pattern [4]. The thickness of the permafrost increases as the altitude rises due to the decreasing temperatures [6,57,58]. Permafrost along the Qinghai–Tibet Highway is distributed continuously over a large area, as depicted in Figure 1c. The distribution of permafrost is influenced by regional geomorphic units, topography, loose sediments, and rock layers [58,59]. Generally, low temperature (with the mean annual ground temperature (MAGT) ≤ −1 °C), stable permafrost is found in high mountain areas, while high temperature (−1 °C < MAGT ≤ 0 °C) [60], unstable permafrost is mainly distributed in high plains and river valleys [57]. Thaw zones are mainly distributed in river valleys at lower elevations.

3. Materials and Methods

3.1. Surface Uplift or Settlement Calculation

LiCSAR, an acronym for looking into continents from space with synthetic aperture radar, is a sophisticated operating system specially designed for the large-scale interferometric processing of Sentinel-1 data [41]. It is designed to automatically generate unwrapped and geocoded interferograms, as well as coherence coefficient maps, with a spatial resolution of 0.001° (about 110 m at the equator) in the WGS-84 coordinate system [47,52]. The LiCSAR processing chain uses a variety of custom tools and algorithms source codes for the core processing functions, which are based on the advanced commercial SAR processing software, GAMMA (https://www.gamma-rs.ch/software, accessed on 26 July 2023), and consist of four main processes, i.e., Preparation of Frame Epoch SLC, Resampling to RSLC, Formation of Differential Interferograms, and Unwrapping Interferograms [42]. These data products cover large areas and are freely and openly accessible for downloading through a dedicated online portal (https://comet.nerc.ac.uk/COMET-LiCS-portal/, accessed on 29 March 2023). All LiCSAR products’ Frame ID information used in this paper is attached in the Supplementary Materials, i.e., the Txt: Frame ID file.

LiCSBAS is an open-source software package developed by the Collaborative Management of Environmental Targets (COMET) for InSAR time series analysis [41,42]. It is integrated with the LiCSAR automated Sentinel-1 data InSAR processor [42,47]. The software package can be accessed at https://github.com/yumorishita/LiCSBAS/ (accessed on 29 March 2023, version 1.5.11). By utilizing the freely available LiCSAR products, users can obtain the results of InSAR time series analysis while reducing processing time and disk space [47,61]. Within the LiCSBAS processing scheme (Figure 2), interferograms with multiple unwrapping errors are automatically identified and removed by phase-closing loops (steps 1-2 in the time series analysis process). Multiple noise indices are employed to mask the results and obtain reliable LOS time series deformation accumulation maps and deformation velocity maps (steps 1-5 and 1-6 in the time series analysis process). Moreover, the Generic Atmospheric Correction Online Service for InSAR (GACOS, http://www.gacos.net/, accessed on 29 March 2023) [46] can be easily implemented to remove atmospheric delay noise (steps 0-3 in preparing the stack of UNW data).

The LiCSBAS toolkit was utilized to obtain the surface deformation data of the study area from 2015 to 2022 (see Table S1 for details), including the LOS cumulative deformation and deformation rate of the ascending and descending tracks. Subsequently, a two-dimensional deformation decomposition on the same or similar (with a time interval of no more than 6 days) results was conducted to acquire the surface vertical deformation results using the acquisition date of the processed images as the filtering condition (UD deformation calculation process in Figure 2); detailed information on this can be found in

Table 1. Results of the date nodes with line-of-sight ascending and descending deformation for obtaining vertical deformation.

Area One	Areas Two, Three, and Four	Area One	Areas Two, Three, and Four	Area One	Areas Two, Three, and Four
20 February 2017	25 February 2017	21 August 2019	26 August 2019	2 October 2020	7 October 2020
4 March 2017	9 March 2017	2 September 2019	7 September 2019	14 October 2020	19 October 2020
9 September 2017	14 April 2017	13 December 2019	5 January 2020	26 October 2020	31 October 2020
21 April 2017	26 April 2017	12 January 2020	17 January 2020	11 February 2021	16 February 2021
24 April 2017	29 September 2017	24 January 2020	29 January 2020	22 August 2021	27 August 2021
11 March 2018	16 March 2018	5 February 2020	10 February 2020	3 September 2021	8 September 2021
28 April 2018	3 May 2018	17 February 2020	22 February 2020	9 October 2021	14 October 2021
27 June 2018	2 July 2018	24 March 2020	29 March 2020	21 October 2021	26 October 2021
21 July 2018	26 July 2018	5 April 2020	10 April 2020	2 November 2021	7 November 2021
7 September 2018	12 September 2018	11 May 2020	16 May 2020	26 November 2021	1 December 2021
13 October 2018	18 October 2018	23 May 2020	28 May 2020	8 December 2021	13 December 2021
12 December 2018	17 December 2018	16 June 2020	21 June 2020	1 January 2022	6 January 2022
29 January 2019	3 February 2019	28 June 2020	3 July 2020	25 January 2022	30 January 2022
30 March 2019	4 April 2019	15 August 2020	20 August 2020	6 February 2022	11 February 2022
11 April 2019	16 April 2019	27 August 2020	1 September 2020	18 February 2022	23 February 2022
16 July 2019	21 July 2019	8 September 2020	13 September 2020	14 March 2022	19 March 2022
28 July 2019	2 August 2019	20 September 2020	25 September 2020	26 March 2022	31 March 2022

.

Based on the LOS deformation results obtained from the ascending and descending orbits, the surface vertical deformation values can be effectively calculated by utilizing the SAR satellite flight geometry parameters while disregarding the north–south deformation, which is generally considered to be less sensitive in InSAR deformation monitoring [20,62]. The method to calculate these values involves the integration of the obtained SAR data with the corresponding flight geometry parameters, resulting in the efficient calculation of the desired deformation values. The specific formula for the calculations is as follows:

$$\left[\begin{array}{c}{d}_{insar}^A\\ d_{insar}^D\end{array}\right]=\left[\begin{array}{cc}\cos {\theta }_{inc}^A& -\mathrm{c}\mathrm{o}\mathrm{s}{\alpha }_{azi}^A\sin {\theta }_{inc}^A\\ \cos {\theta }_{inc}^D& -\mathrm{c}\mathrm{o}\mathrm{s}{\alpha }_{azi}^D\sin {\theta }_{inc}^D\end{array}\right]\left[\begin{array}{c}{d}_u\\ d_e\end{array}\right]$$

(1)

where $d_{insar}^A$, $d_{insar}^D$, ${\theta }_{inc}^A$, ${\theta }_{inc}^D$, ${\alpha }_{azi}^A$, and ${\alpha }_{azi}^D$, $d_u$, and $d_e$ represent the LOS direction results of the ascending InSAR monitoring, the LOS direction results of the descending InSAR monitoring, the incidence angle of the ascending SAR radar, the incidence angle of the descending SAR radar, the azimuth angle of the ascending SAR radar, the azimuth angle of the descending SAR radar, the vertical deformation information, and the east–west deformation information, respectively. Equation (1) can be used to solve the east–west and vertical deformation rates and the cumulative deformation information for the same time period (i.e., nodes). This process can be quickly implemented using a Python script (steps 2-5 in UD deformation calculation in Figure 2), i.e., LiCSBAS_decomposeLOS.py of the LiCSBAS package tools (https://github.com/yumorishita/LiCSBAS/wiki/4_other_tools#licsbas_decomposelospy, accessed on 29 March 2023).

For a comprehensive understanding of the principles and advantages of LiCSBAS, please refer to the literature [42,47,61]. The method and the corresponding parameter settings employed in this paper are elaborated in

Table 2. Detailed information and key parameters of the deformation result processing.

Workflow	Step	Description	Parameters
Prepare stack of UNW Data	0-1. Download	Retrieving GeoTIFF files of UNW from the COMET-LiCS web portal based on the frame ID	-f: Frame ID, -s: 20141001, -e: 20220331, --get_gacos: y, --n_para: 12
	0-2. Convert (and Downsample)	Converting the GeoTIFF files of UNW and COH to float32 and uint8 formats, respectively	-n: 1, --n_para: 12
	0-3. GACOS	Applying a tropospheric correction to the UNW data using GACOS data	-g: Path to the dir containing all GACOS data, --n_para: 12
	0-4. Mask UNW	Masking specified areas or areas with low coherence in the UNW data	-c: 0.2, --n_para: 12
	0-5. Clip UNW	Clipping a specified rectangular area of interest from the unw and cc data	-g: 92.87/95.23/34.45/36.23 (area one), --n_para: 12
Time Series Analysis	1-1. Quality Check	Assessing the quality of the UNW data and identifying bad interferograms based on average coherence and coverage	-c: 0.05, -u: 0.3
	1-2. Loop Closure	Identifying bad UNW by checking loop closure and determining a preliminary reference point that contains all valid UNW data and exhibits the smallest RMS of loop phases	-l: 1.5 rad, --n_para: 12
	1-3. SB Inversion	Inverting the SB network of UNW to obtain the time series cumulative deformation and velocity using the NSBAS approach	--n_unw_r_thre: 1, --gpu: y
	1-4. Bootstrap	Calculating the standard deviation of the velocity using the bootstrap method and STC	--mem_size: 8000, --gpu: y
	1-5. Mask TS	Creating a mask for the time series deformation using several noise indices	-c: 0.05, -u: 1.5, -v: 100, -T: 1, -g: 10, -s: 5, -i: 50, -l: 5, -r: 2
	1-6. Filter (and Deramp) TS	Applying a spatio-temporal filter (high-pass in time and low-pass in space) with a Gaussian kernel, similar to StaMPS	-s: 1, -r: 2, --hgt_linea: y
UD deformation calculation	2-1. Date filtering	Filtering the measurement area with both ascending and descending orbital deformation time points	Implementation through R language conditional functions
	2-2. Cum2vel	Calculating the velocity and its standard deviation from the cumulative deformation of the time series	-s: 20170225, -e: 20220331, --vstd: y
	2-3. Cum2flt	Generating a float32 file that represents the cumulative displacement over a specified date period derived from the original time series cumulative deformation	-d: each of the date notes, -m: 20170225
	2-4. Flt2geotiff	Converting the filtered velocity and time series cumulative deformation results from a float32 format image file to a GeoTIFF file	--a_nodata: −9999
	2-5. Decompose LOS	Decomposing 2 (or more) LOS displacement data into EW and UD components, assuming no displacement in the NS direction	-f: Text file containing input GeoTIFF file paths of LOS displacement (or velocity), E component, and N component (Format: dispfile1 Efile1 Nfile1 dispfile2 Efile2 Nfile2…), -r: cubic
	2-6. Mosaic	Consolidating multiple raster datasets into a new raster dataset, such as the velocity and cumulative deformation results	mosaic operator: mean
	2-7. Clip	Extracting a portion of the mosaic velocity and cumulative deformation based on the boundary (*.shp data) of the Qinghai Tibet Engineering Corridor	use input features for clipping geometry: yes

Abbreviations: NUW, unwrapped interferometric phases; COH, coherence map; GACOS, generic atmospheric correction online service; RMS, root mean square; SB, small baseline; NSBAS, new small baseline subset; STC, spatio-temporal consistency [42]; TS, time series cumulative deformation; LOS, line-of-sight; EW (E), east and south; UD, up and down; NS (N), north and south; *.shp, vector file in ESRI Shapefile format; StaMPS, a software package to extract ground displacements from time series of synthetic aperture radar (SAR) acquisitions (https://homepages.see.leeds.ac.uk/~earahoo/stamps/, accessed on 26 July 2023).

. The source code utilized in the processing is provided in the Supplementary Materials as the Code: Batch_LiCSBAS file.

3.2. Post-Processing Analysis of Results

3.2.1. Deformation Velocity Partitioning

The Gaussian distribution, also known as the normal distribution, is a widely observed probability distribution in various natural phenomena [63,64], as well as in deformation. The normal distribution function is mathematically defined as follows:

$$f(x)=\frac{1}{\sqrt{2\pi }\delta }{e}^{-\frac{1}{2}{\left(\frac{x-\mu }{\delta }\right)}^2}$$

(2)

where $x$ is the observed variable, and in this paper, it is the vertical deformation rate; $\mu$ and $\delta$ are the mean and standard deviation of the observations, respectively. The probability density function of the normal distribution is shown in Figure 3.

According to the characteristics of the normal distribution [63,65], the vertical deformation rate within the engineering corridor can be classified as (min, μ − 3σ), (μ − 3σ, μ − 2σ), (μ − 2σ, μ − σ), (μ − σ, μ + σ), (μ + σ, μ + 2σ), (μ + 2σ, μ + 3σ), and (μ + 3σ, max) in order to indicate strong settlement, subsidence, slight settlement, stability, slight uplift, uplift, and strong uplift. In this way, the surface subsidence or uplift is classified into three classes: stable zone, sub-stable zone (including slight settlement or slight uplift), and unstable zone (including uplift or settlement and strong settlement or uplift).

3.2.2. Cumulative Deformation Trend Detection and Significance Test

To test for trends in the time series variables, several methods can be employed, such as linear slope trend analysis [66,67] and Sen slope estimation [68,69]. The former involves a regression analysis based on long time series data, while the latter employs nonparametric statistics for robust trend estimation, providing advantages such as high computational efficiency and insensitivity to statistical errors and niche data [65]. Therefore, in this study, the Sen slope estimation method was selected to explore trends in the cumulative deformation variables over time. Typically, Sen slope estimation is combined with the Mann–Kendall test to evaluate the significance of changes. The calculation of Sen’s slope estimation is as follows:

$$\beta =Median\left(\frac{d_j-d_i}{j-i}\right)\forall j>i$$

(3)

where $d_i$ and $d_j$ are time series morphological variables (i.e., raster data and morphological variables of the corresponding individual grid points); $Median\left(d\right)$ is the median operation of the observed cumulative morphological variables; $\beta$ is the trend estimation parameter, where greater than 0 indicates an increasing trend of surface morphological variables, less than 0 indicates a decreasing trend of surface morphological variables, and equal to 0 indicates no change.

The Mann–Kendall (M-K) test is a nonparametric statistical test that was originally proposed by Mann in 1945 and further refined by Kendall and Sneyers, and it is widely used to predict long-term trends in the time series data of various elements (e.g., hydrometeorology, including temperature, precipitation, and evapotranspiration) [65]. The method has low statistical requirements for observed data and is suitable not only for normally distributed variables but also for non-normally distributed variables [68,69]. In addition, not only is it suitable for variables with linear trends but also for variables with nonlinear changes; it is not affected by missing values and outliers, and it has been very widely used in trend significance tests for long time series data [65,70]. In addition, it is also applicable in the significance analysis of surface deformation [71]. For the time series deformation variable monitoring results of $D_i$, $i=\mathrm{1,2},3,...,n$, define the annotated test statistic $Z$ for the original hypotheses: H0, the data in the series are randomly aligned, i.e., there is no significant trend; H1, there is an upward or downward trend in the series. The M-K statistical test is calculated as follows:

$$Z=\left\{\begin{array}{c}\frac{S-1}{\sqrt{Var\left(S\right)}},S>0\\ 0,S=0\\ \frac{S+1}{\sqrt{Var\left(S\right)}},S<0\end{array}\right.$$

(4)

$${\sum }_{i=1}^{n+1}{\sum }_{j=i+1}^{n}sign\left(d_j-d_i\right)$$

(5)

$$sign\left(d_j-d_i\right)=\left\{\begin{array}{c}1,\left(d_j-d_i\right)>0\\ 0,\left(d_j-d_i\right)=0\\ -1,\left(d_j-d_i\right)<0\end{array}\right.$$

(6)

$$Var(S)=\frac{1}{18}[n(n-1)(2n+5)-{\sum }_{i=1}^m{t}_i(t_i-1)(2t_i+5)]$$

(7)

where $d_i$ and $d_j$ are the time series morphological variable data (i.e., morphological variables of individual cells corresponding to the raster data $D_i$ and $D_j$); $n$ is the number of cumulative morphological variable periods; $m$ is the number of knots (recurring data sets) in the series; $t_i$ denotes the number of repeated data in the $i$th set of repeated data sets; and $Z$ is the M–K test statistic.

By finding the critical value in the normal distribution table, at a given significance level α, the original hypothesis H0 is rejected when $\left|Z\right|>Z_{1-\frac{\alpha }{2}}$, which means that the trend is considered significant; when $\left|Z\right|<Z_{1-\frac{\alpha }{2}}$, the original hypothesis H0 is accepted, which means that the trend is not significant. At a given significance level α, the critical value $Z_{1-\frac{\alpha }{2}}$ is $\pm$1.96; moreover, when $\left|Z\right|$ is greater than 1.65, 1.96, or 2.58, this means that the trend passes the significance test with a 90%, 95%, and 99% confidence level, respectively. The categories for determining the significance of the trend are shown in

Table 3. Significance categories for trends in cumulative surface form variables.

β	Z	Trend Characteristic
β > 0	Z > 2.58	Extremely significant increase
	1.96 < Z ≤ 2.58	Significant increase
	1.65 < Z ≤ 1.96	Microsignificant increase
	Z ≤ 1.65	No significant increase
β = 0	Z	No change
β < 0	Z ≤ 1.65	No significant reduction
	1.65 < Z ≤ 1.96	Slightly significant reduction
	1.96 < Z ≤ 2.58	Significant reduction
	Z > 2.58	Extremely significant reduction

.

3.2.3. Trend Prediction of Surface Deformation Change

Autosimilarity and long-range dependence are commonly observed phenomena in nature and are frequently used in research fields such as meteorology, hydrology, geology, and seismology. Analyzing the autosimilarity and long-range dependence of time series variables can facilitate the effective prediction of future changes [72,73,74]. Hurst’s index is a quantitative method that effectively describes the autosimilarity and long-range dependence of time series variables. It was first proposed and utilized by the British hydrologist Hurst in 1951 to quantify the relationship between the flow and storage capacity of the Nile reservoir [65,75,76]. Since then, Hurst’s index has been used in various computational methods, such as the absolute value method, aggregated variance method, wavelet analysis, periodogram method, residual analysis, and rescaled polar variance analysis (R/S) [77]. Previous studies have suggested that Hurst indices estimated by R/S analysis and wavelet analysis are more reliable than those estimated using other methods [65,73]. Thus, this paper used the R/S analysis method, which is also the earliest method used by Hurst, to estimate the Hurst index. The definition and formula for this method with R/S analysis are as follows:

The cumulative form variable time series $D_i$ corresponds to the raster endomorphic variables $d_i,i=\mathrm{1,2},3,\dots ,n$. For any positive integer $m$, define the time series as follows:

$$∆d_i=d_i-d_{i-1}$$

(8)

$$\overline{∆d_{\left(m\right)}}=\frac{1}{m}{\sum }_{i=1}^m∆d_i$$

(9)

$$X\left(t\right)={\sum }_{i=1}^m(∆d_i-\overline{∆d_{\left(m\right)}})$$

(10)

$$R\left(m\right)=maxX(t,m)-minX(t,m)$$

(11)

$$S(m)={\left[\frac{1}{m}{\sum }_{i=1}^m{(∆d_i-\overline{∆d_{\left(m\right)}})}^2\right]}^{\frac{1}{2}}$$

(12)

$$\frac{R}{S}={\left(cm\right)}^H$$

(13)

In the formulas, $R$ is the polar deviation; $S$ is the standard deviation; and $c$ is a constant. The monitored time series cumulative variables are divided into $n$ subseries, where $i=\mathrm{1,2},3,\dots ,n$; $m$ is any positive integer and $0<m<n$; $H$ is the Hurst index; $X(t,m)$ is the cumulative deviation of $m$ subseries, where $1\le t\le m$; and $\overline{∆d_{\left(m\right)}}$ and $∆d_i$ are the mean and difference series, respectively.

For Equation (13), $R\left(m\right)/S\left(m\right)\cong R/S$, and the presence of $R/S\propto m^H$ indicates the existence of the Hurst phenomenon for the analyzed time series cumulative deformation variables, which can be obtained by taking the logarithm of both sides of Equation (13) and fitting them in a double logarithmic coordinate system using the least squares method [65,78].

The Hurst index (H) is a statistical parameter that quantifies the autosimilarity and long-range dependence of time series data, where the value range is [0, 1]. A value of H = 0.5 indicates that the time series cumulative quantity is a random sequence with independent distribution and limited variance, lacking long-term regularity. Conversely, when 0< H < 0.5, the time series cumulative-shaped variable exhibits inverse persistence, meaning that future changes will be opposite to past changes. When 0.5 < H < 1, the time series cumulative-shaped variable exhibits sustainability, where future changes will be consistent with past changes. The magnitude of the Hurst index indicates the strength of the inverse persistence (closer to 0) or sustainability (closer to 1) of the time series. The R/S analysis method is a reliable way to estimate Hurst indices, and it is widely used in various research fields, such as meteorology, hydrology, geology, and seismology, for the effective prediction of future changes based on the analysis of the autosimilarity and long-range dependence of time series observed variables [65,73,77].

4. Results Analysis

4.1. Results of the Accuracy Verification

In order to validate the accuracy of the vertical deformation measurements acquired by InSAR, a comparison was conducted between the InSAR results and the field monitoring data on vertical deformation. Additionally, a comparison was conducted between the one-dimensional line-of-sight deformation results based on ascending and descending orbits. The field monitoring data utilized consisted of observations collected from 16 observation points (OPs) located along the QTH, as presented in

Table 4. Geographical location and deformation rate (mm/yr) of the site monitoring points.

Observation Sites	Longitude (E)	Latitude (N)	UD Deformation Rate	Asc. LOS Deformation Rate	Des. LOS Deformation Rate
OP1	94°03.081′	35°37.020′	−0.769	−3.331	1.795
OP2	93°57.795′	35°33.109′	1.046	−2.941	4.818
OP3	93°43.561′	35°30.132′	−7.807	−8.141	−1.896
OP4	93°34.098′	35°24.548′	−3.664	−7.398	2.993
OP5	93°26.776′	35°21.839′	−8.299	−10.236	−1.760
OP6	93°26.678′	35°21.819′	−8.516	−10.558	−1.690
OP7	93°06.678′	35°12.258′	−2.244	−3.813	1.442
OP8	93°02.521′	35°08.303′	−3.599	−4.759	−1.312
OP9	92°53.914′	34°40.346′	−0.419	−0.476	2.700
OP10	92°44.608′	34°34.532′	−6.525	−5.204	−4.624
OP11	92°43.568′	34°28.656′	−2.019	−2.704	0.163
OP12	92°25.838′	34°12.968′	−0.999	2.308	−2.960
OP13	92°20.386′	34°00.675′	−6.992	−4.581	−6.343
OP14	92°14.064′	33°46.399′	−6.570	−4.333	−5.852
OP15	91°56.752′	33°23.874′	−1.828	−2.628	9.093
OP16	91°45.164′	33°04.292′	−2.527	1.334	−4.150

Note: In this study, we used the abbreviation OP to refer to the observation point, UD to represent the vertical direction (up and down), Asc. to indicate the ascending orbit, and Des. to represent the descending orbit. The term LOS stands for the line-of-sight direction, and the deformation rate is expressed in mm/yr. It is important to note that the time span of each deformation rate analyzed in this study ranged from 25 February 2017, to 31 March 2022.

. The field monitoring method involved measuring the relative elevation difference of the monitoring points with respect to a reference point using an infrared level. The data were collected once a month, with each measurement being taken around the 21st of the month, and the deformation monitoring results from 16 points were utilized for comparative purposes in this study. The selected time span for the analysis was from 22/23 May 2017 to 22/23 December 2019. These results were used to compare and validate the InSAR measurements.

Table 4 presents the vertical deformation rates of the InSAR-based measurements at 16 monitoring points. The analysis indicates that, with the exception of OP2, where the surface showed a slight uplift at an uplift rate of 1.046 mm/yr, the remaining 15 points exhibited subsidence. Among these points, the subsidence rates at OP1, OP9, and OP12 were less than or equal to 1 mm/yr. In contrast, the subsidence rates at OP4, OP7, OP8, OP11, OP15, and OP16 ranged from 1.8 mm/yr to 3.7 mm/yr. Moreover, the subsidence rates at the remaining six points exceed 6.5 mm/yr, with OP5 and OP6 showing subsidence rates exceeding 8 mm/yr. Notably, the deformation results based on single ascending or descending LOS-oriented measurements were more variable and, in several instances, inconsistent, with positive and negative signs of deformation rates observed. This inconsistency suggests that the LOS-oriented deformation results incorporated geometric parameters, such as the radar incidence and azimuth angles, that impacted the strength of the SAR signal. Relying solely on one-dimensional ascending or descending LOS-oriented deformation information may, therefore, yield inaccurate results. In summary, the findings suggest that subsidence is occurring at most of the monitoring points, with varying rates of deformation observed. The results based on single LOS-oriented measurements, however, should be interpreted with caution, as they may be affected by various factors, including the geometric parameters of the SAR system.

Figure 4 compares the cumulative results of the time series deformation based on the vertical direction with the field monitoring data, revealing good agreement between the two datasets, although some differences are evident at certain time periods. Several factors account for these differences. First, the time ranges of the data monitoring are not entirely consistent, particularly regarding the starting dates. Second, the reference points used for the two deformation monitoring results are not in the same location. Specifically, the reference points for the InSAR measurements were determined automatically by the LiCSBAS toolkit processing steps 1-3, while the results of the field monitoring are related to the set reference points. Third, the InSAR vertical deformation results were determined based on the ascending and descending orbital data, which require temporal consistency, thus reducing the temporal resolution. Fourth, if there is a large horizontal displacement of the ground surface, especially north–south deformation, this will lead to inconsistencies in the observed results. Conversely, on-site monitoring occurs once a month at a relatively fixed time, with the data collection having a higher temporal resolution compared to the InSAR measurement results. Finally, the monitoring data used were from the railroad roadbed annex, which is significantly influenced by train operations and manual ballast filling, leading to differences in the results. In summary, while the cumulative results of the time series deformation based on the vertical direction agree well with the field monitoring data, certain differences exist due to variations in the monitoring time ranges, reference points, temporal resolution, and external influences on the railroad roadbed annex. These factors should be considered when interpreting the results of deformation monitoring.

The in situ monitoring data obtained at various observation points (OPs) along the QTR indicated significant settlement at most points (Figure 4), namely, OP3, OP4, OP5, OP6, OP10, OP13, OP14, OP15, and OP16, which is consistent with the vertical deformation results obtained using InSAR technology. The field monitoring data also revealed slight settlement or stability with a fluctuating lift phenomenon at OP1, OP2, OP7, OP8, and OP12. These findings are well corroborated by the InSAR vertical results, with discrepancies not exceeding 30 mm, which is smaller than the errors between the onsite and single LOS results, thus demonstrating the accuracy of the vertical deformation derived by combing the ascending and descending orbit results. Furthermore, the comparison of the one-dimensional LOS deformation results of the ascending and descending tracks with field monitoring data indicated high agreement in certain regions, such as OP1, OP3, OP8, OP10, OP13, and OP14, where the trends of the ascending and descending tracks were consistent. However, the error with the field monitoring was larger than the vertical error, exceeding 30 mm. In certain regions, such as OP4 and OP15, the LOS deformation trends of the ascending and descending tracks were inconsistent, indicating that reliance solely on one-dimensional LOS-oriented deformation information could lead to incorrect conclusions. This discrepancy was attributed to the different radar beam-ground geometry relationship during satellite imaging and is discussed in greater depth in subsequent sections.

Upon examining Table 4 and Figure 4, it is apparent that while some areas indicate surface uplift in terms of the deformation rate, the time series cumulative deformation results from both in situ monitoring and InSAR measurements exhibit little or stable uplift, or even a slight decrease in certain cases, such as OP2 and OP9. Nevertheless, in most instances, the deformation rate can accurately depict the long-term variation of a site. Conversely, the in situ monitoring data utilized for comparison are subject to several perceived factors, including inconsistent observation time periods and varying reference point selection. Despite these uncertainties, the vertical deformation and field monitoring trend exhibit high consistency, with errors predominantly within 10 mm and a maximum not exceeding 30 mm, which is of higher accuracy than that of a single LOS direction. Therefore, by combining the ascending and descending track data, the obtained vertical deformation can accurately reflect the actual settlement or uplift of the ground surface.

4.2. Characteristics of the Deformation Results

For the engineering corridor, according to the results (Figure 5), the vertical deformation rate was concentrated between −27.068 mm/yr and 18.586 mm/yr, and the deformation values in this part of the area were mostly between −7.502 mm/yr and 5.542 mm/yr, followed by −14.024 mm/yr and −7.502 mm/yr and 5.524 mm/yr and 12.064 mm/yr. The mean value of the vertical deformation rate results was −1.06 mm/yr, and the median and standard deviation were 0.136 mm/yr. The mean value of the vertical deformation rate was −1.06 mm/yr, and the median and standard deviation were 0.136 mm/yr and 5.904 mm/yr, respectively, which indicates that the surface deformation condition of the Qinghai–Tibet Engineering Corridor is approximately similar to the data processing range, with most areas in a slight uplift or stable state and a few areas in a severe uplift or severe settlement state, but the settlement rate in the settlement area is larger and needs more attention, especially in the area close to the road project. This is because excessive vertical settlement can have a huge impact on the normal operation of the road and even cause traffic accidents.

The distribution of the deformation rate for the engineering corridor is presented in Figure 6, with the exception of the surrounding area due to the low coherence of the descending track data in the annex of the Tongtian River. The areas along the highway and railroad exhibiting significant settlement include the Budongquan–Sonandaji Protection Station–Qumar River section, located in close proximity to Salt Lake in Hoh Xili. The breach of the Zonag Lake dike in 2012 resulted in a sharp increase in the water volume of the downstream Salt Lake area, altering the original thermal condition of the ground surface, thus leading to serious degradation of underground ice and permafrost and a noticeable surface settlement. Furthermore, the Tanggula Mountain section exhibited strong subsidence, with the subsidence rate exceeding −27.068 mm/yr in some areas. Later sections provide a detailed analysis of the areas where the subsidence rate exceeded −20 mm/yr. Conversely, the section from Tongtian River to Yanshiping Town to Wenquan exhibits obvious uplift, with the uplift rate ranging from 0 to 12.064 mm/yr.

By combining the results in Figure 7 and Figure 8, it can be seen that the cumulative deformation variables in the engineered corridor from 25 February 2015 to 31 March 2022 ranged from −725.739 mm to 313.117 mm, with most areas ranging from −11.526 mm to 20.938 mm. The mean, median, and standard deviation of the cumulative deformation variables were −3.683 mm, 1.170 mm, and 27.974 mm, respectively. These results suggest that most of the areas in the engineered corridor are undergoing slight uplift, but in the areas that experienced subsidence, the amount of subsidence was greater, resulting in an overall mean cumulative deformation that was negative. The time series cumulative deformation map shows a high level of consistency with the deformation rate. Specifically, subsidence was most pronounced in the areas around the Suonandajie Protection Station, Tongtian River, and Tanggula Mountain, while uplift was obvious in the section from Yanshiping Town to Wenquan.

4.3. Surface Stability Evaluation and Early Warning

The vertical deformation rate of the engineering corridor was analyzed and graded based on the normal distribution characteristics and significance, as shown in Figure 9a. The majority of the area, accounting for 77% of the whole engineering corridor, is classified as a stable area with a deformation rate ranging between −6.964 mm/yr and 4.844 mm/yr (Figure 10). The sub-stable area occupies 17.7% of the total area, including a slightly uplifted area of 7.7% and a slightly subsided area of 10%. The unstable area comprises strongly subsided (1.2%), subsided (3.2%), strongly uplifted (0.1%), and uplifted areas (0.8%), accounting for 5.3% of the total area. The subsided area (4.4%) is predominant, and the uplifted area only accounts for less than 1% of the total area. Overall, the project corridor is evidently experiencing subsidence, which poses several potential disaster risks.

The spatial distribution analysis reveals the presence of five prominent areas of subsidence (see Figure A1). Area one (Figure A1a) stretches from north to south in an orderly fashion along the Aqinggangqian Rijiu–Kunlun Mountain range. This area displays a clear subsidence strip in an almost east–west direction, with numerous glaciers distributed throughout the region, particularly near the location of the Kunlun tunnel. Area two (Figure A1b) is situated around the Salt Lake in Hoh Xili, with a subsidence belt running in a northwest–southeast direction, encompassing the Salt Lake, Suonandajie protection station, Qingshui Lake, and Qumar River. There is also an apparent subsidence of the river terrace from Baladacai Qu to the Qumar River confluence. The erosion of flowing water and serious weathering along the slope of the surface rock layer are the primary causes of subsidence in this region. Area three (Figure A1c) comprises the three small lake regions, namely, Goulucuo Lake, Zhaliwacuo Lake, and Chacuo Lake. Numerous runoffs, such as Luorigou Qu, Xiacang Nige Qu, and Naoduozhuorou Qu, are present in the subsidence region, which also contains a series of subsidence funnels, similar to those observed in area two. Area four (Figure A1d) is located on the southern side of the Tongtian River (Yanshiping section) and is divided into two subregions by the highway and railroad routes. The north side of the highway extends from the center at Qiamutelang Qu to the southwest at Qielou Qu, while the south side stretches from the center of Naruo Qu to Simeilongrenma village. This region is characterized by a series of subsidence areas and subsidence zones, with several subsidence funnels primarily caused by the weathering of rocks along the river terrace slope. Area five (Figure A1d) is located from Cangzhan village to the Toujoura Mountains, stretching from the northwest to the Touerjiu Mountains (Koulong). This region has a substantial distribution of runoff from Lema Layong, Riana Cangbu, Nadongdaer, Doxiong Qu, Touliju Qu, Tangrong Qu, etc. The subsidence is particularly evident at the bottoms of mountain river valleys flowing through the region, and this may be related to the erosive action of flowing water and the degradation of permafrost due to changes in the thermal state of the local area, resulting in subsidence [37,58]. Apart from these five prominent areas, some smaller subsidence regions exist, such as those along Wudaolang (Porao Congqiong–Porao Congqie–Chage Jialong–Ajia Chage–Ajia Longya), the Xiushui River, near the test section of the Beiluhe Station (34°49.747′N, 92°56.109′E, http://beiluhe.ac.cn/guanyuwomen/, accessed on 26 July 2023), and the Tanggula Mountains near the town of Kaixinling.

In comparison to the large area of subsidence, the area of surface uplift is relatively smaller and primarily concentrated in specific regions. These regions include Xidatan– Wangkun–Kunlun Mountain Tunnel, Kunlun River–Sanchahe Big Bridge, near Hongliang River, Tongtian River (Gari section)–near Gaizeng Mengjiang Macuo Qinma, close to Nimaquan, Yanshiping–Wenquan–highway and railway division, the north side of Tanggula Pass, and the area south of Amdo. Our analysis of Google images reveals that the uplifted areas are mostly located in the river terraces, at the foot of slopes, and in some small areas around lakes. These regions include Kunlun River, Qingshui Lake–Chumar River corner, Chumar River, Hongliang River, Riachi Qu, east side of Ria Chiqu, Zhaliwacuo Lake, Chacuo Lake, Yaxicuo Lake, Tuotuo River, and Tongtian River. We observed that the more curved the river in these regions, the more serious the uplift. These observations suggest that the surface of the project corridor area is affected by erosion and deposit caused by flowing water, leading to surface subsidence in some areas (mainly erosion) and sedimentation (mainly deposit) in others, which results in the observed uplift [7,33].

The findings of the Sen trend analysis and M–K significance test, based on the time series cumulative deformation, are presented in Figure 9b. A comparison of Figure 9a,b and Figure 10a reveals a high degree of agreement, whereby regions with negative deformation rates exhibit an extreme significant decrease and regions with positive deformation rates demonstrate an extreme significant increase. Both results passed the significance test, indicating that in the unstable region, both subsidence and uplift reveal a rise in cumulative deformation. Based on the statistical results in Figure 10b, the trend of highly significant deformation in the entire region encompassed 78.8% of the total area: a 39.8% highly significant decrease and a 39% highly significant increase. Significant areas occupied 5.4% of the total area, with the decreasing and increasing areas each accounting for 2.7%. The slightly significant increase and decrease areas occupied 1.3% of the total area each, amounting to 2.6%. Finally, the insignificant area constituted 13.2% of the total area, with the increasing and decreasing areas each accounting for 6.6%.

The Hurst index results are presented in Figure 9c and Figure 10c, indicating that 64.7% of the areas within the engineering corridor will exhibit an opposite trend in the future, while 35.3% of the areas will continue to exhibit the same trend. Combining the spatial distribution maps in Figure 9a,c, it can be inferred that for the five large subsidence regions (see Figure A1), the future subsidence trend in the location of area one will stabilize, leading to a reduction in subsidence. The Salt Lake of Hoh Xili in area two will continue to subside in the future, which is in line with the current trend. Area 3 will continue to experience subsidence in the future. Areas four and five, located in the vicinity of the Tongtian River and Doxiong Qu, respectively, will experience alleviation and stabilization of subsidence in the future.

Currently, several regions are experiencing strong uplift, including the Tongtian River annex–Gaizeng Mengjiang Macuo Qinma, the vicinity of YanShiPing–Nianyan Qu–Bumade–Quezai Qu–Wenquan–highway and railroad division, and the southern regions of Amdo Mailong, Dopuer Qu, and Jiangqing Qu, among others. The analysis predicts that these regions will continue to experience uplift in the future, with a strong tendency for surface elevation.

4.4. Linear Traffic Engineering Stability Evaluation and Early Warning

To investigate the settlement or uplift patterns along the linear infrastructure of the Qinghai–Tibet Highway and the Qinghai–Tibet Railway, as well as their future trends, we selected the highway–railway vector boundary from Xidatan (the starting point of the highway) to the Amdo section. We then extracted the corresponding vertical deformation rate and Hurst index data with a 1 km interval. The results of this analysis are presented in Figure 11.

In Figure 11, the y-axis on the left indicates the vertical deformation rate, the y-axis on the right indicates the Hurst index, and the x-axis shows the mileage number. The starting and ending points of the highway are Xidatan and Amdo, respectively, while for the railroad, they are Yuzhufeng Station and Amdo Station, respectively. The horizontal lines in the figure correspond to the deformation rates, ranging from −18.772 mm/yr to 16.652 mm/yr, as indicated by the different colors and line styles (solid red, dashed red, dashed green, solid black, solid green, solid light-blue, and solid blue). Additionally, the blue dashed line denotes a location with a Hurst index of 0.5.

Based on an analysis of Figure 11, it can be observed that for both the Qinghai–Tibet Highway and the Qinghai–Tibet Railway, a majority of the sections show a vertical deformation rate of less than 0 mm/yr, indicating that they are subsiding (below the solid, black line). The areas that exhibit instability are primarily in the subsiding regions (above the dashed, green line and below the solid, green line), with most areas showing slight subsidence (between the solid, green line and the solid, blue line). Similarly, the uplifted regions are mostly slightly uplifted (between the dashed, green line and the solid, red line), with only a few areas displaying strong subsidence (below the solid, blue line). The locations of strong uplift and strong subsidence are concentrated near K3264 and K1133 on the highway and near K1335 on the railway, respectively. The settlement in the surrounding area was found to better characterize the deformation of the traffic engineering corridor, except for the influence of its own load. Overall, the stability along the Qinghai–Tibet Railway is better than that of the Qinghai–Tibet Highway.

Taking into account the Hurst index, it can be inferred that the region where the highway is experiencing strong subsidence, specifically the K2951–K2952 section, will continue to subside in the future. However, the K3369–K3372, K3377–K3380, and K3398 areas are expected to exhibit a decreasing trend in the subsidence rate, leading to the stabilization of the surface in the future. As for the region where the railroad is undergoing strong subsidence, which is located near K1335, it is expected to stabilize again in the future. In the area of strong uplift, both roads and railroads are anticipated to stabilize in the future.

Regarding the stable and sub-stable areas, for the Qinghai–Tibet Highway, sections in or near K2927–K2941, K2949–K2972, K3021, K3048–K3058, and K3099 will maintain a significant subsidence rate in the future, rendering them unstable. Conversely, in the vicinity of K3045 and K3289, an upward trend will continue in the future, resulting in an unstable state. With respect to the Qinghai–Tibet Railway, settlement rates near K1010–K1016, K1027–K1047, and K1288 will remain high in the future, leading to an unstable surface. Meanwhile, an upward trend will continue in the future near K1362, also rendering the surface unstable.

5. Discussion

5.1. Necessity of Surface 2D Deformation InSAR Monitoring

The characterization of surface deformation based solely on one-dimensional ascending or descending line-of-sight deformation information is limited due to the actual complexity of the surface and the characteristics of SAR satellite imaging (Figure 12). The LOS direction does not necessarily reflect the actual situation of surface deformation. Thus, it is essential to consider the actual situation of the deformation. In Figure 12a, where the direction of the surface deformation is exactly the same as the LOS direction, the InSAR measurement results can precisely reflect the real deformation of the surface. However, when the surface deformation and direction are perpendicular to the LOS direction measured by the SAR satellite (Figure 12b), the InSAR measurement cannot capture the surface deformation information at all. Generally, the surface deformation direction is at an arbitrary angle to the LOS direction (Figure 12c). In such cases, the LOS direction deformation can represent the surface deformation to some extent and is a vector superposition of the real 3D surface deformation in the LOS direction [79,80].

In order to obtain accurate deformation information of the ground surface, it is necessary to reconstruct the three-dimensional deformation field of the surface by using at least three different directions of the InSAR monitoring results or equivalent a priori information [79,80,81,82,83]. Although the launch of more SAR satellites into space has provided the possibility of InSAR-based reconstruction of the real 3D deformation field of the ground surface, in practical applications it is still challenging to reconstruct the real 3D field due to the differences in the wavelength, geometric parameters, and spatial and temporal resolutions of the SAR data provided by different satellites [79]. This is particularly true for the reconstruction of 3D deformation fields with long time series. Furthermore, the inaccessibility of SAR data further limits the reconstruction of real 3D deformation fields.

Although the results of the surface deformation presented in this paper are not strictly three-dimensional, they represent a significant improvement over traditional one-dimensional LOS deformation results. Specifically, the vertical deformation results obtained for the entire permafrost zone in the Qinghai–Tibet Engineering Corridor, while ignoring the north–south deformation, represent a substantial advancement. Moreover, the validation results indicate that the vertical deformation results obtained through the approach presented here are more accurate and reliable than 1D LOS deformation results and, therefore, more practical for real-world applications. Notably, the obtained results did not incorporate the geometric characteristics of the SAR satellites, which facilitates a comparison across the results obtained from the different satellite platforms. Furthermore, the use of ascending and descending orbit LOS direction solutions is a more scientifically robust approach than simply projecting 1D LOS direction deformation results onto the vertical direction. Ultimately, the InSAR processing approach presented in this paper has the potential to reduce the level of expertise required for InSAR deformation monitoring, thus improving and expanding the application of InSAR monitoring in various industries and contexts.

In the future, it is necessary to introduce more heterogeneous data, such as GPS data, and other methods to process the deformation data obtained from the (additional) SAR data and carry out the reconstruction of the real 3D deformation field of the ground surface in combination with the results of the InSAR deformation of ascending and descending orbits.

5.2. Study of Surface Lift Driving Mechanism and Deformation Pattern

Previous studies have demonstrated that the gradual ground subsidence and uplift observed in the permafrost zone of the Qinghai–Tibet Engineering Corridor can be attributed to the gradual freezing and thawing processes that occur within the permafrost environment. Specifically, ground subsidence is seasonal in nature and closely associated with the thawing and subsequent volume contraction of the active layer due to the phase transition from solid, frozen water to liquid during high summer temperatures. Conversely, during freezing periods characterized by low winter temperatures, ground uplift occurs due to the reverse phase change of water from a liquid to a solid state [10,11,35]. In situations of climate warming or local disturbances, excess ice or ice-rich sediment conduction during the thawing of the uppermost permafrost layer can lead to additional long-term ground subsidence, which can serve as an indicator of permafrost degradation [16,33,35].

Further research is required to investigate how environmental factors affect deformation. Based on the results obtained in this study, it is necessary to conduct an analysis of the correlations between large-scale surface deformation information and environmental factors in the study area [84]. For example, Zhao et al. utilized InSAR technology to obtain surface deformation information for the permafrost section from Xidatan to Amdo in the Qinghai–Tibet Engineering Corridor. They analyzed the correlations between surface deformation and three influencing factors, namely, volumetric ice content, average annual ground temperature, and active layer thickness. The results indicated a strong positive correlation between the surface deformation and volumetric ice content, as well as the average annual ground temperature, with the bias correlation coefficients exceeding 0.7. Additionally, a moderate positive correlation was observed with active layer thickness, with an average bias correlation coefficient of 0.42 [85].

The literature [33] analyzed the relationship between LOS-oriented deformation rates, seasonal deformation magnitudes, and 14 potential environmental factors in the permafrost zone of the Tibetan Plateau (central Tibetan Plateau, along the Qinghai–Tibet Highway from the northern Kunlun Mountains to the southern boundary of the permafrost zone and seasonal permafrost) obtained based on InSAR using the Geodetector method, and the results showed that the slope and latitude are the main drivers of the spatial pattern of seasonal deformation within the study area. Latitude is the most important factor driving the spatial pattern of linear deformation, which is similar to the latitudinal pattern (the strong subsidence or uplift signals along the latitudinal direction in the linear deformation) observed from linear deformation maps [33]. In addition, based on the results of vertical deformation, the development of suitable surface deformation patterns in permafrost areas needs to be further developed to improve the theory of surface deformation information inversion models in permafrost areas.

In the future, it will be necessary to enhance and develop the use of vertical surface deformation information in permafrost areas obtained using the methods provided in this paper. It is necessary to study the deformation characteristics by means of zoning [13]; to use different models (linear and nonlinear) to invert the surface deformation information and to analyze their differences and accuracy [10,11]; and to carry out modeling of permafrost deformation in different periods (freeze-up and thaw-down periods), especially the relationship between deformation amplitude, deformation rate, and time series cumulative deformation variables [33].

5.3. LiCSBAS Processing LiCSAR Results Applicability

The InSAR data processing procedure is intricate, and its outcomes are contingent upon the operators’ professional expertise. Furthermore, discrepancies arise in the results obtained by operators with diverse knowledge backgrounds [15,86]. Therefore, it is imperative to investigate the reliability of the results presented in this paper since the approach employed for InSAR processing differs from the traditional methodology.

Adequate input data serves as the fundamental assurance for achieving high-quality and accurate output. In this study, the vertical deformation results are derived from the solution obtained through the utilization of ascending and descending orbit data. Therefore, ensuring the input of high-quality (i.e., high coherence) ascending and descending orbit data is crucial to obtaining precise vertical deformation output [22]. Figure 13 illustrates the average coherence coefficient maps for all the ascending and descending orbit data employed in this study, spanning from August 2015 to March 2022.

The distribution map of mean correlation coefficients (Figure 13) reveals that both ascending and descending data exhibit coherence coefficients greater than 0.3, except for regions containing lakes and smaller areas. Notably, extensive contiguous areas in the northern part of the engineering corridor exhibit high coherence (coherence coefficient exceeding 0.7), aligning with the severe permafrost degradation observed in the northeastern region of the Qinghai-Tibet Plateau in recent years. The absence or degradation of permafrost leads to predominantly bare ground surfaces with less pronounced freeze-thaw cycles, resulting in a relatively stable surface with high coherence.

The statistical analysis conducted on an image-by-image basis illustrates that for the ascending track data, the average coherence coefficients range from 0.025 to 0.917, with 90% of the area exhibiting coefficients between 0.357 and 0.525. Similarly, for the descending track data, the average coherence coefficients range from 0.025 to 0.917, with 90% of the regions displaying coefficients ranging from 0.343 to 0.506. These findings indicate the substantial coherence of the input interferogram data utilized in this study, thereby ensuring accurate output results.

In the LiCSBAS toolkit’s InSAR time series analysis, the Bootstrap method can be employed in steps 1-4 (time series analysis process in Figure 2) to calculate the standard deviation of deformation velocity (referred to as VSTD) [42,87]. This calculation relies on the linear cumulative displacement results obtained from steps 1-3 (time series analysis process in Figure 2). The VSTD serves as an indicator of the deformation results’ reliability and the suitability of the linear deformation model [42]. A higher VSTD estimation suggests the presence of noise or nonlinearity within the displacement time series. This indicates that the time series inversion using the linear deformation model is not appropriate, emphasizing the need to enhance the accuracy of the obtained deformation results. Conversely, a lower VSTD implies the inverted deformation results’ high accuracy and the strong applicability of the linear model employed for the inversion of time series deformation. Figure 14 presents the VSTD values corresponding to the ascending and descending orbit data used for InSAR monitoring spanning from 25 February 2017, to 31 March 2022.

Figure 14 illustrates that the VSTD remains below 0.4 mm/yr in the majority of regions for both ascending and descending orbits, with only a few regions exhibiting a VSTD exceeding 0.6 mm/yr. Remarkably, there is a significant spatial agreement between the two. The statistical analysis conducted on an image-by-image basis reveals that 95% of the Line-of-Sight (LOS)-oriented deformation exhibits a VSTD below 0.375 mm/yr. Furthermore, for the descending orbit data, 95% of the areas display a VSTD below 0.360 mm/yr. These findings demonstrate the high accuracy of the deformation rate results derived solely from the calculation of ascending and descending orbit data, thus underscoring the reliability of the vertical time series deformation results. Consequently, it can be inferred that both the vertical time-series deformation calculation and the deformation rate calculation yield dependable outcomes.

Furthermore, to ensure consistency with the deformation rate results from Chen et al. [33] (available for download at http://www.ncdc.ac.cn/portal/metadata/f17d720f-b765-4342-b11d-f8d81a77f730, accessed on 29 March 2023), the date (from 2014 to 2019) and data processing ranges were constrained using the same parameter method employed in this study. Consequently, the parametric method utilized in this paper produced deformation results aligning with the study range of Chen et al., enabling a direct comparison with their published results (refer to Figure 15).

The findings reveal a significant spatial agreement between the deformation rate distributions presented in this study (Figure 15a) and those reported by Chen et al. (Figure 15b). The regions exhibiting high and low values correspond closely. Specifically, the spatial distribution locations and range sizes of regions displaying LOS-oriented deformation results on the ground surface, whether close to or distant from the SAR satellite, are essentially the same. This correspondence not only attests to the reliability of the ascending and descending LOS deformation results obtained using the methods and parameters employed in this study but also reinforces the dependability of the vertical deformation results.

6. Conclusions

In this paper, the LiCSAR product combined with the LiCSBAS toolkit was used to obtain surface deformation information in the permafrost zone of the Qinghai–Tibet Engineering Corridor for the period from May 2017 to March 2022, and the results were used as the basis for a stability evaluation study of the engineering corridor surface. The following conclusions were obtained:

(1)

Based on the LiCSAR product and LiCSBAS package, we can quickly obtain the surface deformation monitoring results of large-scale and long time series, and the calculations exhibit low consumption of computational resources, high computational efficiency, and the capability to be conducted in batch automation. It can also be used with other toolkits to quickly crop, mosaic, and select the deformation results for time periods of interest. This provides new methods and options for InSAR monitoring in the context of big data, cloud platforms, and cloud computing, and it lays a solid foundation for the development of large-scale surface deformation monitoring in the future.

(2)

The surface of the study area was in a slight settlement state from May 2017 to March 2022, and the vertical deformation rate was mostly distributed in the range of −27.068–18.586 mm/yr, with an average of −1.06 mm/yr. The results of the field monitoring show that the error of the vertical time series’ cumulative deformation was mostly less than 10 mm and the maximum was not more than 30 mm; while the error of the single ascending and descending track monitoring results was mostly more than 50 mm, and there are multiple deformation trend discrepancies. This shows that the vertical deformation results obtained using the same date or similar dates to obtain the deformation results for the ascending and descending tracks can better reflect the real settlement or uplift of the ground surface in the permafrost area.

(3)

A total of 77% of the engineering corridor was in a stable state, with vertical deformation rates between −6.964 mm/yr and 4.844 mm/yr, while 17.7% of the area was in a sub-stable state, wherein, 7.7% of the total area was considered unstable, including settlement rates between −12.868 mm/yr and −6.964 mm/yr, accounting for 10% of the total area, and the slightly uplifted area (uplift rate between 4.844 mm/yr and 10.748 mm/yr) accounting for 7.7% of the total area. The unstable area included an area with a settling rate greater than 12.868 mm/yr and an uplift rate greater than 10.748 mm/yr, accounting for 4.4% and 0.9% of the total area, respectively, totaling 5.3%. There were five large subsidence areas within the project corridor, containing numerous subsidence funnels, while the uplift areas were much smaller and sporadically distributed compared to the subsidence areas.

(4)

The stability of the areas along the Qinghai–Tibet Railway is significantly higher than that of the Qinghai–Tibet Highway, and there are fewer sections located in unstable areas. Four areas with serious settlement and one area with obvious uplift were found along the highway, while only two areas were found to be unstable along the railway, one each for settlement and uplift. The areas that need to be focused on in the future for the Qinghai–Tibet Highway are the five areas of subsidence and two areas of uplift, while the areas of subsidence and uplift along the railroad are areas two and one, respectively. The results obtained based on the method outlined in this paper can provide effective data support and the specific locations of high-risk areas for the safe operation of highways and railroads, as well as effective reference solutions for long-term monitoring and future early warning.