3.1 Datasets

The greatest feature of the railway is its long linear engineering; to monitor and analyze this type of ground object, certain requirements must be met. The first is large area coverage because the road direction is not the same everywhere and, in monitoring, the need for a unified reference standard requires large-scale synchronized measurements. Second, to meet the requirements of InSAR technology, such projects must apply a high-precision InSAR data processing algorithm to ensure high-resolution and fine-scale detection.

The TerraSAR-X data were acquired in strip-map mode with an incidence angle range of 39–40^∘ at horizontal–horizontal (HH) polarization. The potential of the X-band data for detecting higher deformation gradients compared to that of other sensors arises from the high spatial and temporal resolutions of these data. Nevertheless, the coverage of the strip-map mode data is too small to study long, linear engineering projects. Therefore, in this paper, the C-band ASAR data and TerraSAR-X data were used in combination to analyze the stability of the Qinghai–Tibet Railway. The TerraSAR-X data were selected to verify the accuracy of the ASAR T405 data results over the first segment of the railway in Yangbajain, and the ASAR T133 data were used to analyze the deformation of the railway near the Naqu area because the ASAR T405 data could not cover this area completely, and the ASAR T133 data could also verify the accuracy of the ASAR T405 data results over the Naqu area. The data coverage is shown in Fig. 2 with the blue dotted line.

The ASAR T405 data were acquired from August 2003 to September 2010, but there were no data for 2006; therefore, we processed the data in three stages (2003–2005, 2007, and 2008–2010). The ASAR T133 data were acquired from November 2007 to August 2010. The TerraSAR-X data were acquired from December 2011 to November 2012.

3.2 Methodology

The main steps of the FRAM-SBAS (full rank matrix small baseline subset InSAR) method are as follows.

Firstly, the principle of interferogram generation is based on a specific time baseline and space baseline, and the appropriate redundant interferogram is selected to maximize the interferogram coherence. The main constraints are Eq. (1):

where ΔB_⊥ is the vertical baseline of data interference pairs, Δt is the time baseline, and ΔDC is the Doppler frequency difference.

To minimize the spatial and temporal decorrelations, we constructed a baseline network (Fig. 3) using the following criteria: perpendicular baselines shorter than 200 m and a daytime interval baseline of less than 180 d. Each acquisition node in the network has at least two link pairs, meaning that each node has a minimum number of connections with the other nodes (two are used in this paper).

Figure 3Temporal and perpendicular baselines for the interferograms used in this study. Different graphs represent the baselines of different data generated during different time periods. In each figure, the horizontal axis represents the time (year) and the vertical axis represents the perpendicular baseline (m).

Download

Secondly, coherence points are selected. The coherence point is selected based on the principle of full rank matrix, which effectively improves the quality of coherent point selection and provides the basis for subsequent least-squares inversion. By constructing a single set interferogram network, each point can construct the matrix described by Eq. (2), where A is an M∗N dimensional matrix, M is the number of interferograms, and N is the number of images. For any pixel in any interferogram, the coherence is greater than a certain threshold, and the −1 and 1 flags can be set at the corresponding positions of the matrix A of Eq. (2). For example, the first interferogram consists of the first image and the third image; that is, ${\mathit{\delta }}_{\mathit{\varphi }\mathrm{1}}={\mathit{\varphi }}_{\mathrm{3}}-{\mathit{\varphi }}_{\mathrm{1}}$. If the interference of a certain point in the interferogram satisfies the conditions, then the corresponding positions are $A_{\mathrm{11}}=-\mathrm{1}$ and A₁₃=1, and the remaining positions of the first line are 0. Similarly, in the second interferogram, if the coherence of the point in one of the interferograms is less than the coherence threshold, the diversion is set to zero. All interferograms are considered to obtain each point pair matrix, and then the rank of each matrix is determined (Eq. 2). If the matrix is full rank, the point is selected as the coherence point. The method can be used to select points that are coherent in the time series and coherent in the partial time interval but the interference network is connected, thereby increasing the number and precision of the coherent points.

Thirdly, discrete point phase unwrapping is performed. In the FRAM-SBAS method, discrete coherence point data are resampled onto a regular Cartesian grid and phase unwrapping is performed using a network flow method. Then, the phase jump is checked according to the closed ring residual method, and the jump phase is corrected for the jump region.

Because of the change in the water vapor content in the atmosphere, phase artifacts in InSAR images caused by path delays, such as radar signal propagations through the stratified and turbulent atmosphere and ionosphere, frequently degrade the interpretability of the phase and correlation signatures of the terrain. The effect of atmospheric delay consists of three parts: (1) the long wavelength effect of the atmosphere, which is similar to the orbit error effect; (2) the short wavelength effect of the atmosphere (i.e., turbulent atmospheric artifacts); and (3) the vertical stratification of the atmosphere, which causes height-dependent refractivity variations. In this paper, the three phase delays are calculated using a network methodology. The methodology estimates the phase delay for each SAR acquisition; then, each atmospheric artifact is simulated. The proposed method can effectively eliminate the atmospheric phase delay in the interferograms.

Fourthly, orbital and atmospheric error removal is performed. The orbit error removal is performed using the network method proposed by Biggs et al. (2007). The atmospheric error is divided into long-wavelength atmospheric delay error and turbulent atmospheric delay error and terrain-related atmospheric delay error. The three errors are removed using the network methods.

Traditional atmospheric delay phase (APS) estimates are based on a single interferogram (Ferretti, 2001). The atmospheric phase in the interferogram is the difference in atmospheric phase delay between the sub-image and the main image. If one of the two images is used to generate other interferograms, the phase delay signal on the image is also passed to the other interferograms, which also makes a correlation between the two interferograms. In this paper, we will use the network method to estimate the atmospheric delay error of each image acquisition time and then use these estimates to obtain the delay error of a single moment to reconstruct the atmospheric delay error of each interferogram.

After removal of the DEM error and the deformation phase, it can be assumed that the residual phase is mainly caused by the atmosphere. Suppose δφ_j(x,y) represents the residual phase value at (x,y) on the jth interferogram and that $\mathit{\phi }(t_A,x,y)$ and $\mathit{\phi }(t_B,x,y)$ represent the phase values of the imaging moments t_A and t_B at (x,y), respectively. Each interferogram can be expressed by Eq. (3):

Based on a short baseline set network, we can construct equations such as Eq. (4):

where A represents the M⋅N matrix. The element A_kl of the matrix A is defined according to the following rules: if l=t_B, then A_kl=1; if l=t_A, then $A_{\mathrm{kl}}=-\mathrm{1}$; otherwise, A_kl=0. δφ is a known vector of M dimension, representing the number of interferograms M; φ is an N-dimensional unknown vector representing the atmospheric phase values of N imaging moments. Equation (4) can be rewritten as follows:

where δφ_k(x,y) represents the residual phase of interferogram k and the corresponding position is (x, y).

Since the matrix A is the rank-deficient matrix, a unique solution cannot be obtained. Generally, the singular value decomposition (SVD) method can be used to solve the solution and the atmospheric delay at each moment is obtained; then, the phase value of each interferogram is simulated by using Eq. (5). In the calculation of the variance of the residual phase of each interferogram, if the interferogram has the lowest atmospheric variance, the atmospheric phase of the interferogram is assumed to be zero. This constraint is added to Eq. (5) to calculate the atmospheric delay phase of all other image acquisition moments (Li, 2014).

Fifthly, the deformation result is obtained. The interference pattern is settled using the least-squares method to obtain the deformation results of the study area.

The specific procedure is illustrated in Fig. 4.

Figure 4Process flow of the full rank matrix small baseline subset InSAR (FRAM-SBAS) time-series analysis method.

Download

4.1 InSAR results

The deformation information obtained by ASAR T405 during the construction of the Lhasa–Naqu section of the Qinghai–Tibet Railway from 2003 to 2005 is shown in Fig. 5a. The deformation of the study area is very small during this period, and the maximum deformation is approximately 5 mm yr⁻¹.

Figure 5The InSAR results. Panel (a) represents the ASAR T405 results for the Lhasa–Naqu section of the Qinghai–Tibet Railway between 2003 and 2005. Panel (b) represents the ASAR T405 results for the Lhasa–Naqu section of the Qinghai–Tibet Railway in 2007, and the elliptical red dotted line represents the region of large deformation. Panel (c) represents the ASAR T405 data from 2008 to 2010. Panel (d) represents the ASAR T133 data from 2008 to 2010. The red line represents the railway. The black circle represents the main stations along the railway.

In 2007, the Lhasa–Naqu section of the Qinghai–Tibet Railway was functionally completed and opened to traffic. Figure 5b shows the deformation information of the line obtained by ASAR T405 in 2007. The area of the line is obviously variable compared with that before the opening of the railway. In the area circled by the elliptical red dotted line, the deformation is large, with a maximum value of 20 mm yr⁻¹. The operation of the train has a certain impact on the railway, which may be because the roadbed bears the weight of the train, causing the roadbed to become compacted and sink.

From 2008 to 2010, 3 years after the smooth operation of the Qinghai–Tibet Railway, the deformation of the line slowed compared with that in 2007 (Fig. 5c–d). The main deformation area is consistent with the deformation zone in 2007, mainly in the area of Damxung, which may be due to the geothermal exploitation in the area. A comparison of ASAR T405 and ASAR T133 reveals that the maximum deformation in this area is 15 mm yr⁻¹. However, the deformation in the northern part of the railway is relatively stable.

Figure 6The deformation maps of the overlap area of the TerraSAR-X and ASAR data. (a) The ASAR T405 data range from 2008 to 2010. (b) TerraSAR-X data range from 2011 to 2012.

The overlap area between TerraSAR-X and ASAR is located above Yangbajain, and the deformation field obtained from the two datasets is analyzed. Figure 6a indicates that ASAR acquired the deformation field in the region from 2008 to 2010. There is an obvious uplift area at the corner of the railway from 2008 to 2010, and the maximum cumulative uplift is up to 7 mm yr⁻¹. Figure 6b shows the deformation field acquired by TerraSAR-X from 2011 to 2013, and the deformation of the uplift area at the corner is obviously smaller and tends to be stable, but there are two subsidence areas in the lower left corner and upper right corner, where the maximum subsidence is 10 mm yr⁻¹. The reasons for the subsidence area are analyzed by superposing geological hazards.

Figure 7The deformation maps of the railway sequence changes from 2003 to 2010 in the Damxung section. (a) 2003–2005; (b) 2007; (c) 2008–2010.

For the present analysis, a section of the railway in Damxung is selected (Fig. 7). In 2003–2005 (Fig. 7a), the section of the road surface was basically stable. In 2007 (Fig. 7b), the construction of the railway resulted in a region with a large area of settlement, where the maximum settlement rate was 20 mm yr⁻¹. In 2008–2010 (Fig. 7c), the subsidence area was reduced, and the maximum settlement rate was 10 mm yr⁻¹; the surface tended to be stable.

Based on the results of the deformation field obtained for 2009–2010, some high-voltage towers along the Lhasa–Naqu railway are analyzed (Fig. 8). The results show that different sections of the towers have different drops or lifts. The maximum lifting capacity of D306 is 17 mm, and the maximum settlement of D269 is 18 mm. Most of the larger variables are located at the corner of the railway and within a section of the line.

Figure 8Deformation of a high-voltage power tower. (a) The location of the power tower. (b) The time series of the deformation characteristics of the power tower.

In this paper, the SAR data from 2003 to 2012 are analyzed. It is concluded that during the construction of the Qinghai–Tibet Railway, the linear variable along the railway was approximately 10 mm yr⁻¹. After the completion of the traffic, the linear variable along the railway was 4–8 mm yr⁻¹. Li (2012) used SBAS technology to analyze the ENVISAT ASAR data from 1997 to 2010 in the vicinity of Yangbajain Damxung of the Qinghai–Tibet Railway. It was found that the settlement rate near the railway was 2 mm yr⁻¹, and the impact of frozen soil was approximately 10 mm yr⁻¹. Zhang et al. (2017) used Sentinel-1 data to analyze the deformation variables of the Qinghai–Tibet Railway during the period of 2014–2016 in the Qinghai–Tibet Plateau. It is concluded that the settlement rate of the Qinghai–Tibet Railway is approximately −10 mm yr⁻¹ and the settlement rate of the rail-stabilized area is approximately −5 mm yr⁻¹. Ma et al. (2011) and Dong et al. (2013) found that the overall settlement rate of the Qinghai–Tibet Railway subgrade is < 10 mm yr⁻¹. At the junction of the fracture, we verified this finding with GPS, and the GPS result was highly consistent with the deformation field acquired by InSAR. Chen et al. (2012) used C- and L-band small baseline SAR interferometry to analyze the interaction between permafrost and infrastructure along the Qinghai–Tibet Railway, and the results showed surface motions along the embankment primarily in the range of −20 to +20 mm yr⁻¹.

4.2 Deformation and hazards

The Qinghai–Tibet Railway runs along the Qinghai–Tibet Plateau and its eastern margin. Under the influence of the uplift of the Qinghai–Tibet Plateau, the topography is generally high in the west and low in the east. The Qinghai–Tibet Railway passes through structural units, such as the Naqu orogenic belt and the Lhasa block. The climatic region of the Qinghai–Tibet Plateau has obvious vertical zoning characteristics, large temperature differences between winter and summer, and strong freeze–thaw weathering, which makes the geological environment along the Qinghai–Tibet Railway and adjacent areas sensitive and conducive to the development and occurrence of geological disasters.

Figure 9Deformation and hazards. The magenta rhombi show the geological hazard points.

In this paper, the upper Yangbajain railway section is selected and the geological hazard and deformation fields are superimposed (Fig. 9). The geological hazards occur mostly in areas with large deformation fields, which are also in the corner along the railway. There is basically no distribution of geological hazard points in the stable area of the deformation field, which shows that the stability of the surrounding geological structure is affected by railway construction and other factors, resulting in frequent geological hazards. At the same time, to ensure the safe and normal operation of the railway, it is necessary to perform long-term monitoring of the geological and geomorphological conditions around the railway to avoid casualties and economic losses.

4.3 Deformation and the Bengco fault

Garthwaite et al. (2013) obtained the deformation field within the Qinghai–Tibet Plateau by using the large-scale InSAR technique. The present slip rate of the Bengco fault is 1–4 + 1 mm yr⁻¹. At 90.30^∘ E, the slip rate is 1 mm yr⁻¹, which is 40 km west of the rupture zone of a 1951 earthquake. At 90.75^∘ E, which is located in the western part of the rupture zone, the sliding rate reaches 4 mm yr⁻¹ (Garthwaite et al., 2013). Ryder also used the InSAR method to estimate the deformation rate of the Bengco fault. The right-lateral slip deformation at 90.4^∘ E is approximately 1±1 mm yr⁻¹, and that at 90.9^∘ E is approximately 4±1 mm yr⁻¹. It is believed that this fracture deformation is mainly caused by the post-earthquake stick-slip effect of two earthquakes in 1951 and 1952 (Ryder et al., 2014). Taylor and Peltzer (2006) studied a series of V-shaped conjugate tectonic belts using the InSAR technique. The present slip rates of several tectonic belts in the central and western regions are consistent with GPS results (Taylor and Peltzer, 2006).

Figure 10Bengco fault and the distribution of the SAR imagery. The blue dotted line shows the extents of the ASAR images. The red line shows the Bengco fault. The purple line shows the railway.

Figure 11Bengco fault slip rate. The red line shows the Bengco fault. The purple line shows the railway. The red circle shows the GPS location. The black dotted line shows the profiles of the Bengco fault. The purple line shows the highway. The left figure shows the Bengco fault slip rate. The right figure shows the detailed information on profiles.

Using the 32 ASAR image scenes between 2004 and 2010 (two adjacent tracks, T176 and T405) (Fig. 10), the same treatment method as described in Sect. 3.2 was used to obtain the deformation rate of the Bengco fault zone. Three profiles, A-A', B-B', and C-C', are selected in the deformation rate map to analyze the slip rate of the Bengco fault (Fig. 11). It can be seen from the results that the formation rate of the Bengco area is between 1 and 3 mm yr⁻¹. The rate in the eastern segment is approximately 2–3 mm yr⁻¹, and that in the western segment is approximately 1–2 mm yr⁻¹. The intersection of the A-A' section line and the fault essentially coincides with the intersection of the Qinghai–Tibet Railway and the fault-breaking fault zone. From the results of the section deformation rate, it can be seen that the fault slip has little effect on the overall deformation of the railway.