**1. Introduction**

The repeated alternating action of frost heaving and thawing settlement of frozen soil in permafrost area will cause damage to the local geological environment, and then easily lead to geological disasters such as landslide, debris flow, foundation rupture, and collapse [1,2], especially linear engineering such as railway and highway [3–5]. Uneven deformation of the roadbed is a common and serious disaster in linear engineering, which seriously affects the operation of linear engineering such as high-speed railways and highways. The Qinghai–Tibet Plateau (QTP) is the largest permafrost region except polar regions [6]. Approximately 610 km of Qinghai–Tibet railway is laid in the permafrost area, crossing national nature reserves such as Hoh Xil, Sanjiangyuan, Qiangtang. In view of the particularity and importance of the location of the Qinghai–Tibet railway, it is of great practical significance to monitor the surface deformation of the Qinghai–Tibet corridor in permafrost area [7–17].

Traditional geodetic methods, such as leveling and global positioning system (GPS) measurement, can achieve high-precision measurement of surface deformation. However,

**Citation:** Liu, H.; Huang, S.; Xie, C.; Tian, B.; Chen, M.; Chang, Z. Monitoring Roadbed Stability in Permafrost Area of Qinghai–Tibet Railway by MT-InSAR Technology. *Land* **2023**, *12*, 474. https://doi.org/ 10.3390/land12020474

Academic Editor: Xiaoyong Bai

Received: 5 January 2023 Revised: 10 February 2023 Accepted: 13 February 2023 Published: 14 February 2023

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

due to the special location and harsh environmental conditions of the Qinghai–Tibet plateau, both traditional leveling and GPS measurement require a lot of human and material costs, which is difficult and inefficient. Although it will not cause this problem to predict the stability of Qinghai–Tibet line foundation in combination with historical disaster data, field survey data [18], settlement index or allowable bearing capacity index [19], the accuracy is rough and low, and the prediction of influence factors cannot be better applicable to the existing situation under the external human intervention.

The wide application of synthetic aperture radar interferometry (InSAR) makes up for the shortcomings of the above method. InSAR monitors ground deformation by analyzing the phase information of two aperture radar images [20]. It has been widely used in the ground deformation caused by earthquakes, volcanoes, glaciers, landslides, and land subsidence [21–24]. With the launch of new SAR satellites such as Sentinel-1A/B, researchers can more easily obtain SAR images. Rich SAR data sources and accessibility promote the application of InSAR technology in the field of permafrost deformation monitoring. Due to its large coverage of SAR images and the relatively short return visit time, the interference coherence has been greatly improved, which better meets the accuracy requirements of deformation monitoring such as linear engineering in permafrost area, and it is used to monitor frozen soil deformation. Researchers have already demonstrated the ability to use InSAR technique to detect to this freeze/thaw-related ground motion over permafrost regions. As shown in Table 1, current researches using InSAR technology to monitor frozen soil mainly focus on the Arctic, Qinghai–Tibet Plateau and other regions of the world. It can be seen that the application of InSAR technology to permafrost monitoring has broad prospects.

**Arctic Qinghai–Tibet Plateau Other Regions** Zwieback et al. [25], Bartsch et al. [26], Strozzi et al. [27], Rudy et al. [28], Liu et al. [29–33] Zhou et al. [34], Xu et al. [35], Zou et al. [36], Xiang et al. [37], Zhang et al. [38–44], Wang et al. [45,46], Reinosch et al. [47], Lu et al. [48], Chen et al. [49], Daout et al. [50] Chen et al. [51], Rouyet et al. [52], Antonova et al. [53], Li et al. [54], Liu et al. [29]

**Table 1.** Research on monitoring frozen soil with InSAR technology.

As can be seen from the above table, the research on permafrost in the Qinghai–Tibet plateau in recent years has been fruitful. However, over the years, most research on the frozen soil deformation along the Qinghai–Tibet railway has focused on the section from Wudaoliang to Tuotuohe, especially the Beiluhe area, and less so on the Chaerhan Salt Lake section in the north and the section from Yangbajing to Dangxiong in the south. The scope of the existing research has been very limited, and the overall deformation of the permafrost section of the Qinghai–Tibet railway has not been monitored. In addition, the existing research only involves local disaster characteristics, and there is a lack of systematic investigation and analysis of the overall disaster characteristics and laws. Most of the deformation research on the frozen soil section of the Qinghai–Tibet railway is limited to the results of more than a decade. Due to the long time, it is difficult to reflect the current situation of roadbed deformation, which is not conductive to the discovery of the existence and potential risk of deformation of railway roadbed. At the same time, the previous research results have not been compared and discussed, and the methods and effects of the measures have not been qualitatively and quantitatively evaluated. The relevant research only stays at the level of producing deformation results. Therefore, many years after the completion of the Qinghai–Tibet railway project, it is very necessary to

compare the deformation monitoring results in recent years with the previous monitoring results to understand whether the deformation disaster has changed.

Most of the studies have only used permanent scatterers to obtain the deformation characteristics of ground objects along the Qinghai–Tibet railway. The scattering characteristics are directly affected by external environmental factors such as soil moisture and surface water, and almost no ground objects can maintain stable backscattering characteristics for a long time. At the same time, the vegetation coverage in the study area is low, mostly exposed sandy or mucus saline soil, and the interference coherence is high in the short term. This phenomenon makes it more difficult to obtain persistent scattering points. Therefore, it is difficult to obtain high-density stable scattering points only by obtaining traditional persistent scatterer (PS) points. Distributed scatterers (DS) refers to ground objects, mostly bare ground and sparse vegetation that resolve the backscattering coefficients of all scatterers in a unit roughly the same. The spatial density of the measurement points is increased on the region characterized by DS while retaining the high-quality information obtained using the PS technique on deterministic targets. PS and DS were properly combined to increase the density of measurement points and further improve the coherence point density and parameter estimation accuracy [55]. Considering the limitations of traditional persistent scattering points, we obtain a high-quality and high-density time series deformation point set by combining PS points and distributed scattering points (DS) with limited time baselines, which greatly improves the defects of the above methods so as to obtain the existing deformation of railway roadbed and potential risks along the route.

In this study, we used the free obtained Sentinel-1A Radar image to monitor the overall deformation of the permafrost section of the Qinghai–Tibet railway by combining PS points and DS points, compared and discussed the deformation of local areas for many years to monitor the uneven settlement of railway roadbed in space and the seasonal changes in the time domain.

#### **2. Materials and Methods**

#### *2.1. Study Area*

The repeated alternating action of frost heaving and thawing settlement of frozen soil in permafrost area will have a serious impact on linear projects such as Qinghai–Tibet railway. Approximately 610 km of Qinghai–Tibet railway is paved in permafrost area, from Nachitai to Ando station. The study area covers 32◦20 –35◦94 , 91◦42 –94◦62 , with an altitude of 3514–5544 m.

The thickness of frozen soil is approximately 60~120 m, and the thickness of active layer (ALT) is between 0.8 and 4 m, with an average of approximately 2 m It is a typical continental climate. It is cold in winter (to −20.7 ◦C), warm in summer (to +22.7 ◦C), and the annual average temperature is −0.98–1.03 ◦C. The ground is generally frozen for 7 months (October to April of the next year) [56]. The annual precipitation is between 191 and 485 mm, mainly concentrated in summer (May to October). The elevation fluctuates from north to south [57], and its topographic characteristics are shown in Figure 1. There are four monitoring stations close to this section, namely Golmud, Wudaoliang, Tuotuohe and Anduo meteorological stations. Land types mainly include glaciers, snow, exposed rocks and other land types [58]. Due to the melting and freezing of frozen soil, the maximum seasonal deformation of the surface can reach 10 cm. The distribution of frozen soil is shown in the upper left corner of Figure 1. The landslide within 1 km and debris flow within 5 km caused by frost heaving and thawing settlement of frozen soil will pose a threat to the railway. This study will take the permafrost section of a 5 km buffer zone as the research object.

**Figure 1.** The green diamond is the meteorological station. The station name is marked in green font. The white mark is the altitude of each meteorological station. The red pentagram is the railway station along the way, and the name is marked in red font. The terrain distribution has been marked in black font in the Figure. The Qinghai–Tibet line passes through mountain permafrost, plateau discontinuous permafrost and island permafrost from north to south. The study area is a permafrost section passing through the blue area, which has higher risk than other sections. The black rectangular box indicates the coverage of the Sentinel-1 image, the red rectangular box is the swath of each scene used in the study, and the black arrow indicates the navigation direction and line of sight direction.

## *2.2. Data*

In this paper, the descending image of Sentinel-1 satellite in 2020 is selected to obtain the risk situation of the section of Qinghai–Tibet line, involving 122 images in four maps, as shown in Table 2.


**Table 2.** Relevant parameters of four scenes of radar image.

The pixel sizes in the central azimuth and range directions of the image are 2.32 m and 13.98 m, respectively, and VV polarization mode is selected. The central incidence angle is between 33 and 34◦, and the spatial range covered by each SAR image is approximately <sup>252</sup> × 190 km2. The digital elevation model (DEM) of the study area adopts SRTM-1 with a resolution of 30 m. In addition, the Precision Orbit Data (POD) of Sentinel-1A is provided by European Space Agency (ESA). The data of glacier and frozen soil distribution originate from the "Heihe project data management center" [59].

#### *2.3. Methods*

In this study, permanent scatterers and distributed scatterers (PS, DS) were combined to select highly coherent targets and increase the coverage of coherent points. Compared with PS method, this method has a better effect and higher efficiency in the analysis of foundation deformation along the Qinghai–Tibet corridor. The processing flow chart is demonstrated in Figure 2. The deformation results largely depend on the error correction during interference processing, the selection of coherent points and phase unwrapping. Therefore, this section is expanded with the principles and methods of each part to finally obtain the desired results.

**Figure 2.** Process of joint point selection of PS and DS to obtain high-density point targets and roadbed time series deformation.

#### 2.3.1. Error Correction

High coherence of target means high coherence in space domain and stable phase in time domain. For MT-InSAR, the image will be registered, and the terrain and orbit have generated phase compensation. The residual phase *ϕ<sup>n</sup>* of unwrapped interferogram obtained after differential interference is usually composed of residual terrain, deformation, atmosphere and noise phases [60], which are wound between deformation phases.

$$
\boldsymbol{\phi}^{\rm u} = \boldsymbol{\phi}\_{\rm topo}^{\rm u} + \boldsymbol{\phi}\_{\rm defo}^{\rm u} + \boldsymbol{\phi}\_{\rm aps}^{\rm u} + \boldsymbol{\phi}\_{\rm noise}^{\rm u} + k \cdot 2\pi\mathsf{\tau},\tag{1}
$$

where *ϕ<sup>n</sup> topo* is the residual topographic phase; *ϕ<sup>n</sup> defo* is the deformation phase; *<sup>ϕ</sup><sup>n</sup> aps* is the atmospheric phase screen, indicating the signal delay caused by weather conditions; *k* is the integer fuzzy number; *ϕ<sup>n</sup> noise* is the phase noise caused by time decorrelation, error registration, uncompensated spectral offset decorrelation, orbit error, soil moisture and thermal noise.

The residual topographic phase of point p is as follows:

$$
\!\!\!\!p\_{p,topo}^{\mu} = k\_p^z \cdot \epsilon\_{p\prime}^z \tag{2}
$$

where *k<sup>z</sup> <sup>p</sup>* = **<sup>4</sup>***πbn <sup>λ</sup>sinθRn* is height-phase coefficient and *<sup>z</sup> <sup>p</sup>* is elevation error; *bn* is the baseline of the nth image relative to the reference (main) image; *θ* is the local angle of incidence; *λ* is the carrier wavelength; *Rn* is the (zero Doppler) distance between the target and the nth orbit acquisition.

The deformation phase generated by the displacement at point p can be effectively divided into two components:

$$
\boldsymbol{\upvarphi}\_{p, \text{defo}}^{\mu} = \boldsymbol{k}^{t} \boldsymbol{v}\_{p} + \mu\_{\text{NL}}.\tag{3}
$$

where *v* is constant speed (line of sight); *k<sup>t</sup>* = **<sup>4</sup>***πtn <sup>λ</sup>* is the time phase factor; *tn* is the time baseline; *μNL* is the phase term due to possible nonlinear motion.

Atmospheric phase screen (APS) can greatly reduce the influence by considering the phase difference between two nearby points. Therefore, APS calibration can also be carried out by establishing Delaunay triangulation network.

For the two PS candidate points connected, the measured tomography signals *y***1**(*m*) and *y***2**(*m*) can be expressed as [61]

$$\begin{cases} y\_1(m) = \mathcal{g}\_1(m, s\_1) \exp[j \mathcal{Q}\_{APS}(m)] \\ y\_2(m) = \mathcal{g}\_2(m, s\_2) \exp[j \mathcal{Q}\_{APS}(m)]' \end{cases} \tag{4}$$

where <sup>∅</sup>*APS***1**(*m*) and <sup>∅</sup>*APS***2**(*m*) represent APS present in *<sup>y</sup>***1**(*m*) and *<sup>y</sup>***2**(*m*), respectively; *s***<sup>1</sup>** and *s***<sup>2</sup>** are the inclined elevations of the two PS points connected, respectively; *g***1**(*m*, *s***1**) and *g***2**(*m*, *s***2**) represent the math ideal measurement value of the two connection points without APS interference.

Because the spatial frequency of APS is low, the candidate points with long spatial distance usually have great differences in APS. At this time, it is difficult to deal with connection points with large APS differences. Therefore, the long connection arc is rejected by setting the distance threshold. When the arc length is short, two adjacent scatterer candidate points have similar APS:

$$
\mathcal{Q}\_{\text{APS}}(\mathfrak{m}) \approx \mathcal{Q}\_{\text{APS}}(\mathfrak{m}).\tag{5}
$$

Therefore, the APS can be easily calibrated by subtracting the phase of the reference point, and the relative tomographic signal **Δ***y* (*m*) of the connecting arc is as follows:

$$
\Delta y(m) = y\_2(m) \exp(-j \cdot \lceil y\_1(m) \rceil) = \mathbf{g}\_2(m, \mathbf{s\_2}) \cdot \exp(-j \cdot \lceil \mathbf{g\_1}(m, \mathbf{s\_1}) \rceil), \tag{6}
$$

where is the operation of phase retention.

### 2.3.2. Selection of Coherent Points

Permanent and distributed scatterers are combined to increase the coverage of coherent points. Ferretti et al. proposed the Dispersion of Amplitude (DA) based on the definition of PS [62]. When the main scatterer exists in the pixel, its phase is mainly determined by the phase of the main scatterer, which is less affected by noise, and the phase standard deviation has the following relationship with the amplitude:

$$
\sigma\_{\varphi} \approx \frac{\sigma\_{A}}{m\_{A}} \equiv D\_{A\prime} \tag{7}
$$

where *σϕ* is the standard deviation of phase; *σ<sup>A</sup>* is the standard deviation of amplitude A; *mA* is the mean amplitude of N SAR images in time dimension; *DA* is the dispersion index.

When extracting DS, considering its statistical distribution characteristics [63], the extracted DS is transformed into extracting several pixels subject to the same statistical distribution, which can be called homogeneous pixels. Then, the recognition of planar targets is transformed into identifying homogeneous pixels first, and then estimating their phase. The Kolmogorov–Smirnov algorithm is used to identify homogeneous pixels. By determining whether their cumulative distribution function (CDF) is the same, we can judge whether they are homogeneous pixels. K–S algorithm defines *DN* as the maximum absolute value of the difference between the two cumulative distribution functions:

$$D\_N = \max\_{-\infty < \mathbf{x} < \infty} |S\_N(\mathbf{x}) - P\_N(\mathbf{x})|\_{\prime} \tag{8}$$

where *SN* and *PN* are the cumulative distribution functions of two different SAR image pixels, respectively. By setting a certain threshold (*Dthr*), when *DN* ≤ *Dthr*, it is considered that the two samples obey the same statistical distribution, that is, they are judged to be homogeneous pixels so as to ensure that the identified homogeneous pixel is directly or indirectly connected with the central pixel P, and there is no independent region. In addition, the algorithm takes the identified homogeneous pixel as DS, which makes the subsequent phase estimation more accurate, so the results are more reliable.

#### 2.3.3. Phase Unwrapping

The accuracy of surface deformation information acquisition mainly depends on phase unwrapping. The main idea of phase unwrapping algorithm based on network flow is to minimize the difference between the derivative of unwrapping phase and the derivative of winding phase. This method can not only greatly reduce the time and space complexity of phase unwrapping algorithm and improve the calculation speed, but also limit the whole error to a small range and prevent the retransmission of error, so as to improve the accuracy of unwrapping results. In this study, the minimum cost flow algorithm is based on irregular networks: first, the phase with high coherence coefficient is extracted as a high-quality phase data set. Next, a Delaunay triangulation is established according to the position of these phase points. Then, the residual points in the triangulation are identified, and the minimum cost flow algorithm is used to connect the positive and negative residual point pairs to establish the branch tangent. Finally, the unwrapping phase value is obtained by the method of adding and subtracting 2nπ around or through the branch tangent.

The phase difference between two adjacent phase points (such as targets p and q) can be written as [64]

$$
\Delta \boldsymbol{\sigma}\_{pq}^{\boldsymbol{\mu}} = \mathcal{W} \Big\{ k\_p^z \Delta \boldsymbol{\epsilon}\_{pq}^z + k^t \Delta \boldsymbol{\sigma}\_{pq} + \Delta \boldsymbol{\omega}\_{pq}^{\boldsymbol{\mu}} \Big\}, \tag{9}
$$

where *W* is the winding operator, Δ*<sup>z</sup> pq* is the relative elevation error, Δ*vpq* is the relatively constant velocity, and Δ*w<sup>n</sup> pq* = *μ<sup>n</sup> pq*,*NL* +*ϕ<sup>n</sup> pq*,*aps* +*ϕ<sup>n</sup> pq*,*noise* is the phase difference between the model and measured values between p and q points in the nth interferogram.

Therefore, the Delaunay triangulation connection network is constructed in the spatial domain as follows. Assuming that the largest connected network is searched after K iterations, and there are *<sup>P</sup>*(*n*) arcs and *<sup>Q</sup>*(*n*) PS points (n = 1, ··· , K) in the nth search result of the connected network, then the maximum connected network contains *P*(*k*) connected arcs and *Q*(*k*) points. The integration of the nth search connection network can be modeled as

$$
\Delta \mathbf{S}^{(n)} = \mathbf{G}^{(n)} \cdot \mathbf{S}^{(n)},\tag{10}
$$

where Δ*S*(*n*) is the relative height of the *P*(*n*) point arc of the nth search connection network; *G*(*n*) is the transformation matrix from point arc to point, which is composed of 0, <sup>−</sup> 1, 1; *<sup>S</sup>*(*n*) is the absolute height of the *<sup>Q</sup>*(*n*) point in the nth search connection network:

$$
\Delta S^{(n)} = \begin{bmatrix}
\Delta S\_1^{(n)} \\
\Delta S\_2^{(n)} \\
\vdots \\
\Delta S\_p^{(n)} \\
p \\
p
\end{bmatrix}\_p \times 1, \mathsf{G}^{(n)} = \begin{bmatrix}
1 & \cdots & -1 & \cdots & 0 & \cdots \\
\vdots \\
\cdots & 1 & 0 & \cdots & -1 & \cdots \\
\vdots \\
\cdots & 0 & \cdots & 1 & -1 & \cdots \\
\vdots \\
\end{bmatrix}\_{p^{(n)}}, \mathsf{S}^{(n)} = \begin{bmatrix}
\Delta S\_1^{(n)} \\
\Delta S\_2^{(n)} \\
\vdots \\
\vdots \\
\Delta S\_{Q(n)}^{(n)}
\end{bmatrix}\_{Q^{(n)}}.\tag{11}
$$

#### **3. Results**

In order to improve the phase solution accuracy, it is necessary to check whether the spatio-temporal baseline meets the requirements, and reduce the impact on the deformation phase. In the study area, 21 swaths were analyzed by MT-InSAR, and a total of 118 pairs of interference pairs were generated. The temporal and spatial baseline distribution and combination mode are shown in Figure 3. The re-entry period of SAR satellite is 12 days, the vertical baseline is less than 130 m, and the coherence is relatively high. The finally obtained surface time series deformation information along the Qinghai–Tibet line is shown in Figure 4. The areas with serious deformation are mainly distributed from WangKun station to Budongquan station (WB) in the north section of permafrost area and from Tanggula station to Za'gya Zangbo station (TZ) in the south section.

**Figure 3.** The center point is the reference main image, each line represents an interference pair, and the X and Y axes represent its spatio-temporal baseline distribution.

**Figure 4.** Overall deformation results obtained by time series InSAR processing. The surface deformation of frozen soil section of Qinghai–Tibet corridor is mainly distributed in WB in the north and TZ in the south, which is represented by red rectangular box in the Figure. The black arrow indicates the Los direction.

#### **4. Discussion**
