Landslide-Hazard-Avoiding Highway Alignment Selection in Mountainous Regions Based on SAR Images and High-Spatial-Resolution Precipitation Datasets: A Case Study in Southwestern China

: Landslides recurrently cause severe damage and, in some cases, the full disruption of many highways in mountainous areas, which can last from a few days to even months. Thus, there is a high demand for monitoring tools and precipitation data to support highway alignment selections before construction. In this study, we proposed a new system highway alignment selection method based on coherent scatter InSAR (CSI) and ~1 km high-spatial-resolution precipitation (HSRP) analysis. Prior to the CSI, we calculated and analyzed the feasibility of Sentinel-1A ascending and descending data. To illustrate the performance of the CSI, CSI and SBAS–InSAR were both utilized to monitor 80 slow-moving landslides, which were identified by optical remote-sensing interpretation and field investigation, along the Barkam–Kangting Highway Corridor (BKHC) in southwestern China, relying on 56 Sentinel-1A descending images from September 2019 to September 2021. The results reveal that CSI has clearer deformation signals and more measurement points (MPs) than SBAS-InSAR. And the maximum cumulative displacements and rates of the landslides reach − 75 mm and − 64 mm/year within the monitoring period (CSI results), respectively. Furthermore, the rates of the landslides near the Jinchuan River are higher than those of the landslides far from the river. Subsequently, to optimize the highway alignment selection, we analyzed the spatiotemporal evolution characteristics of feature points on a typical landslide by combining the − 1 km HSRP, which was calculated from the 30 ′ Climatic Research Unit (CRU) time-series datasets, with the climatology datasets of WorldClim using delta spatial downscaling. The analysis shows that the sliding rates of landslides augment from the back edge to the tongue because of fluvial erosion and that accelerated sliding is highly related to the intense precipitation between April and September each year (ASP). Consequently, three solution types were established in our method by setting thresholds for the deformation rates and ASPs of every landslide. Afterward, the risk-optimal alignment selection of the BKHC was finalized according to the solution types and consideration of the construction’s possible impacts. Ultimately, the major problems and challenges for our method were discussed, and conclusions were given.


Introduction
Throughout history, reliable road networks have been vital local assets, connecting communities and unlocking economic growth [1,2].Since China's Belt and Road (B&R) of governments and researchers around the world, observation networks still suffer from low station density and spatial resolution in mountainous regions [53].Thus, several interpolation methods, such as inverse distance weighting, kriging methods, and regression analysis, are used to generate meteorological data for ungauged areas.However, the station density decides the result accuracy of these methods [53][54][55][56][57][58][59].In general, precipitation data products are released by several climate research organizations, including general circulation models (GCMs), the Climatic Research Unit (CRU), the Global Precipitation Climatology Centre (GPCC), and Willmott and Matsuura (W&M) [60][61][62].Particularly, CRU products include the monthly mean precipitation, which is generated from data obtained from observational stations but has a low spatial resolution (~55 km).Nevertheless, the climate changes drastically per kilometer in the mountainous areas of southwestern China, especially the precipitation, which induces landslide instability.In practice, the ~1 km HSRP needs to be combined with the results of the TS-InSAR to optimize the highway alignment selection.Consequently, it is necessary to spatially downscale and correct CRU monthly mean precipitation products.Previous studies have proved that the delta downscaling framework is suitable for monthly precipitation data downscaling using CRU products [60,[63][64][65][66].
Field investigation and optical remote-sensing interpretation are employed in the conventional highway alignment selection method; sometimes, drilling wells are needed too.The drawback of this method is the high human and economic costs.As a major traffic connection to the Bangladesh-China-India-Myanmar International Economic Corridor and the Silk Road Economic Belt in the Tibetan area of Sichuan Province, the alignment selection problem of the BKHC is severe because of the complicated geological background, strong tectonic movement, rock body rupture, as well as frequent geohazards.To solve the problem of the conventional method, in this paper, we propose a new systematic highway alignment selection method based on CSI and ~1 km monthly HSRP analysis.
First, the FSs of the Sentinel-1A ascending and descending data were calculated.Consequently, based on 56 Sentinel-1A satellite images from September 2019 to September 2021 (before the design was finalized), both CSI and SBAS-InSAR were utilized in a typical section of the BKHC (from Jinchuan to Danba), which has 80 slow-moving landslides identified by optical remote-sensing interpretation and field investigation.The comparison results of the two methods illustrate that CSI has better performance.Subsequently, three solution types were established in our method by setting thresholds for deformation rates (CSI results) and ASPMPs of every landslide.Finally, the risk-optimal alignment selection of the BKHC was finalized according to the solution types and consideration of the construction's possible impacts.

Study Area and Data Sources 2.1. Study Area
As depicted in Figure 1, the BKHC goes from northeast to southwest, accompanied by the G248 national highway.This corridor is located in the cascade zone between the Qinghai-Tibet Plateau and Sichuan Basin, which goes through the tectonic erosion of the high-middle mountain landform, tectonic erosion of the alpine landform, river erosion and accumulation ladder landform, tectonic erosion of the alpine canyon landform, and tectonic erosion of the alpine mountain-plateau landform.Topographically, the BKHC is characterized by deep valleys and rugged mountains, with elevations ranging from 1800 m to 5300 m.As depicted in Figure 2, more than 18 active faults intersect the route's lines.Strong tectonic movements have produced deep, narrow river valleys.Meanwhile, the rainy period in this region usually lasts from April to September.The monthly precipitation in this period far exceeds that in any other months.The complex geological, hydrological, and geotechnical processes provide dynamic conditions, which promote the occurrence of landslides.Along this highway corridor, a typical section from Jinchuan to Danba lies on the river scouring and accumulation ladder landform.Generally, such landforms consist of shorter cliffs and steep slopes.With the river scouring and precipitation, the stability of the vertical sliding surfaces and rock-layer structure becomes weak.Consequently, fully developed landslides are distributed along this section.A 10 km buffer has been taken along this section as the study area.

SAR Data Sources
A set of 56 descending Sentinel-1A images (from September 2019 to September 2021) and corresponding precise orbit datasets (POD) are provided by the European Space Agency (ESA) (Table 1).In addition, the 30 m resolution digital elevation model (DEM) provided by the National Aeronautics and Space Administration (NASA) is employed to estimate and remove the topographic phase.The coverage of the SAR images is shown in Figure 1 with the blue-dashed box.3) [66,67].
As previous studies have proven, CRU datasets exhibit better performance than other similar gridded products.In addition, 323 weather stations across China were employed by the CRU group to generate CRU time-series data (Figure 3) [64][65][66][67].And reference WorldClim datasets comprised monthly mean precipitations at 2 m, which were generated based on 9000-60,000 weather stations located globally using the thin- Along this highway corridor, a typical section from Jinchuan to Danba lies on the river scouring and accumulation ladder landform.Generally, such landforms consist of shorter cliffs and steep slopes.With the river scouring and precipitation, the stability of the vertical sliding surfaces and rock-layer structure becomes weak.Consequently, fully developed landslides are distributed along this section.A 10 km buffer has been taken along this section as the study area.

SAR Data Sources
A set of 56 descending Sentinel-1A images (from September 2019 to September 2021) and corresponding precise orbit datasets (POD) are provided by the European Space Agency (ESA) (Table 1).In addition, the 30 m resolution digital elevation model (DEM) provided by the National Aeronautics and Space Administration (NASA) is employed to estimate and remove the topographic phase.The coverage of the SAR images is shown in Figure 1 with the blue-dashed box.3) [65,66].
As previous studies have proven, CRU datasets exhibit better performance than other similar gridded products.In addition, 323 weather stations across China were employed by the CRU group to generate CRU time-series data (Figure 3) [63][64][65][66].And reference WorldClim datasets comprised monthly mean precipitations at 2 m, which were generated based on 9000-60,000 weather stations located globally using the thin-plate spline interpolation method.Remarkably, cross-validation correlations indicated that these datasets exhibited good performance globally because of the introduction of the satellitederived covariates and distance to the nearest coastal covariates and reflected orographic effects well.
plate spline interpolation method.Remarkably, cross-validation correlations indicated that these datasets exhibited good performance globally because of the introduction of the satellite-derived covariates and distance to the nearest coastal covariates and reflected orographic effects well.

Methods
In this paper, we propose a new systematic highway alignment selection method to finalize the landslide-hazard-avoiding highway alignment of the BKHC, which consists of four steps, i.e., data selection (feasibility calculation of SAR data), CSI and monthly HSRP analysis, and alignment selection (based on the locations, three solution types, and the construction's possible impacts, as depicted in Figure 4).Each of the four steps will be illustrated in detail.

Methods
In this paper, we propose a new systematic highway alignment selection method to finalize the landslide-hazard-avoiding highway alignment of the BKHC, which consists of four steps, i.e., data selection (feasibility calculation of SAR data), CSI and monthly HSRP analysis, and alignment selection (based on the locations, three solution types, and the construction's possible impacts, as depicted in Figure 4).Each of the four steps will be illustrated in detail.

Data Selection
The relationships between the line of sight (LOS) and slope displacements are summarized in Figure 5 [68].The slope orientations and features, which are reflected in the DEM, determine the geometric distortion occurrence.The feasibility of the SAR

Data Selection
The relationships between the line of sight (LOS) and slope displacements are summarized in Figure 5 [67].The slope orientations and features, which are reflected in the DEM, determine the geometric distortion occurrence.The feasibility of the SAR data (FS) is calculated based on the DEM (slope angles and aspects) and the satellite orbit parameters (slant range) [68].Generally, the shadow effect occurs when the beam cannot illuminate a parcel of terrain because of a physical barrier.Likewise, the layover effect indicates the pixels affected by distortion occurring when the front of the radar beam reaches the top of a slope before it reaches its base.And both shadow and layover effects are more likely to occur in far ranges and on slopes facing the same direction as the beam or in near ranges for small incidence angles.Herein, we take the local incidence angle, which is, namely, the angle between the incident radar beam and a line that is normal to that surface, as a basis calculation factor.The local incidence angle was calculated using Equation ( 2) Local incicence angle = arccos( where   denotes the height from the satellite to the center of the Earth,  ℎ denote the geodetic height of the sub-satellite point on the Earth,  1 denotes the near slan range, and   represents the resolution of the slant range. The FS calculation results are classified into 4 classes according to the effects o Hence, the FS can be calculated using Equation (1) [68]; all the calculation factors must be in radians.The slope and aspect can be derived using the DEM, and the Flat_area can be extracted by classifying the slope.The Lay_Shad_mask needed a reclassification from the original data, which was created by assigning a value equal to 0 for the pixels affected by the layover and/or shadow and 1 for the others.The local incidence angle was calculated using Equation (2), where R H denotes the height from the satellite to the center of the Earth, R h denotes the geodetic height of the sub-satellite point on the Earth, L 1 denotes the near slant range, and P r represents the resolution of the slant range.FS = −sin((Slope (rad) × sin(Aspect (rad) + LOS azimuth angle (rad)) The local incidence angle was calculated using Equation ( 2) where R H denotes the height from the satellite to the center of the Earth, R h denotes the geodetic height of the sub-satellite point on the Earth, L 1 denotes the near slant range, and P r represents the resolution of the slant range.
The FS calculation results are classified into 4 classes according to the effects of the terrain's geometry (Table 2).The SAR data composition selection can be determined based on these 4 classes [68].

Value
Type

CSI
As TS-InSAR methods, PS and DS algorithms have been widely used in landslide monitoring.The PS algorithm focuses on the selection of PS targets, while the DS algorithm identifies DS targets and their optimal phases.CSI combines these two targets for further analysis using a standard PSI tool, which can reflect the details of landslide deformations better with enough points [42].

PS Target Selection
As the first generation of the TS-InSAR, the PS algorithm is a mature algorithm that integrates into a variety of commercial or open-source tools.In CSI, StaMPS is employed to select PS points based on amplitude and phase information because of the close relationship between the amplitude stability and phase stability when the coherence is high [42,69].Hence, the amplitude dispersion index (ADI) can first be used to select PS targets and is calculated using Equation (3) as follows: where σ A and µ A denote the mean and standard deviation of a series of amplitude values for one pixel, respectively.When the ADI value of a certain pixel is under 0.4, it can be considered as a PS target.Subsequently, the deformation signal will be assumed to be spatially correlated for each PS target, and the temporal coherence (γ) can be calculated using Equation (4) as follows [69]: where N and φ noise,i denote the number of images and the estimated residual phase noise of ith SLC image, respectively.Certain PS targets will be kept with the original phase values for further time-series analysis when the their γ values are high enough.

DS Target Selection
DS target selection has two steps, including SHP identification and optimal phase estimation, as follows: (1) SHP identification Generally, it is necessary to obtain sufficient and high-temporal-resolution SAR images (at least 2 years of datasets with a 12~48-day revisiting cycle) to help with highway alignment selection, especially in southwestern China, because of the spatiotemporal discorrelation.Thus, a KS test can be used to identify the SHPs, which needs a large number of samples [70].The sample's complex coherence matrix can be estimated using Equation (5) as follows: where N denotes the number of images, x ′ represents the pixels in a fixed window (Λ 25 × 25) centered on pixel x, Y(x ′ ) is the normalized complex scattering vectors, and H indicates a Hermitian transpose.Particularly, only if the number of SHPs with high weights (ω(x, x ′ ) > 0.5) in window Λ is larger than 20, pixel x can be selected as a DS target in the CSI.
(2) Estimation of the optimal phase For each pixel, the optimal phase series Ψ = [Ψ 1 , Ψ 2 , . . . ,Ψ N ] T can be estimated accord- ing to the previous section.Ψ 1 is set as 0, and the maximum likelihood estimation (MLE) of Ψ can be calculated using Equation ( 6) as follows: where η = [0, e jΨ 1 , . .., e jΨ N ] T , and the symbol "o" denotes the mathematical operator of the Hadamard product between two matrices obtained using Equation (5).Then, a phaselinking approach can be utilized to solve this equation, which can be expressed as a closed form in Equation (7) as follows [71]: where k is the iteration step.
To evaluate the quality of the estimated optimal phases, the goodness-of-fit index (GF) can be calculated using Equation (8) as follows: The GF, illustrated as the extension of the temporal coherence, ϕ mn , is the phase value of item (m, n) in the coherence matrix, and Ψ m and Ψ n are the estimated optimal phases.Then, those DS targets with high GF values (higher than the predefined threshold) will be selected as the final DS targets.

Combination of PS and DS Targets
PS and DS pixels can be connected to form a Delaunay triangulation network for phase analysis after first dropping the common pixels.Subsequently, the phase is first corrected for the spatially uncorrelated part of the look-angle error.The final time-series deformation and rate can be retrieved after the 3D phase unwrapping and spatiotemporal filtering [69,72].All these procedures can be implemented in the StaMPS.

HSRP Analysis
To obtain the monthly ~1 km HSRP datasets, delta downscaling is a well-suited method, including four steps [53].First, a climatology dataset is constructed for each month based on 30 ′ CRU and WorldClim time-series data.Second, the 30 ′ anomaly time-series data are derived for precipitation based on the 30 ′ CRU time-series data and the constructed precipitation dataset, which can be expressed using Equation (9) as follows: where Ano PRE (yr, m) are the anomalies for the precipitation, PRE(yr, m) are the absolute precipitation values, CRUClim_PRE(m) is the 30 ′ climatology for the precipitation, and m and yr correspond to month (from January to December) and year, respectively.Third, the bilinear interpolation method is employed for the 30 ′ anomaly grids at each time.The 0.5 ′ spatial resolution datasets (~1 km) will be retrieved to match the WorldClim data.Finally, the high-spatial-resolution anomaly time-series dataset is transformed to an absolute climatic time-series dataset based on the WorldClim data using Equation (10) as follows: where PRE(yr, m, 0.5 ′ ) are the absolute precipitation values at a 0.5 ′ spatial resolution (~1 km HSRP), Ano PRE (yr, m, 0.5 ′ ) denotes anomalies at the 0.5 ′ spatial resolution for the precipitation, and WorldClim_PRE(m, 0.5 ′ ) are climatology datasets from WorldClim at a 0.5 ′ spatial resolution for the precipitation.

Alignment Selection
As the final step of our method, all the accurate results above should be considered comprehensively.The annual deformation rate is one of the most direct factors reflecting the activities of landslides.Usually, when the rate's absolute value is higher than 10 mm/yr, the landslide is considered as an active event [73].And when the rate's absolute value is higher than 30 mm/yr, the landslide is considered as a vigorous active event (The annual deformation rate of every landslide's geometric center point is selected in our method.).
As for the HSRP, because our study area is in a typical Tibetan Plateau monsoon climate zone, which only has a low mean annual rainfall [42], the rainfall intensities in the rainy season will decide the stability of landslides in such areas [74].Hence, the total precipitation from April to September of every landslide within the monitoring period (ASPMP) was extracted as our judgment factor based on the ~1 km HSRP.Furthermore, for highways at the local scale, the HSRP is close to a normal distribution, and the ASPMP's mean value (ASPMPM) for all the landslides within the monitoring period has a higher reference price than any others in the statistics.

Results and Analysis
As depicted in Section 3, our proposed method was utilized to typically section the BKHC (from Jinchuan to Danba), which has 80 slow-moving landslides identified by optical remote-sensing interpretation and field investigation, based on 56 Sentinel-1A satellite images from September 2019 to September 2021 (before the design was finalized).As a result, the risk-optimal alignment selection of the BKHC was finalized.

Results of the FS and Data Selection
As shown in Figure 6a,b, the FSs of the ascending and descending Sentinel-1A images were calculated.Almost half of these landslides were located in the low-effect area in the descending SAR data, while almost half of them were located in the shadow or high-effect area in the ascending SAR data.Thus, we selected 56 descending images of Sentinel-1A (2019-2021) as the only data source for the CSI. Figure 6c reveals that these landslides are distributed in the river valley; consequently, these landslides have developed well because of the river's perennial scouring and dense rainfall.
the BKHC (from Jinchuan to Danba), which has 80 slow-moving landslides identified by optical remote-sensing interpretation and field investigation, based on 56 Sentinel-1A satellite images from September 2019 to September 2021 (before the design was finalized).As a result, the risk-optimal alignment selection of the BKHC was finalized.

Results of the FS and Data Selection
As shown in Figure 6a,b, the FSs of the ascending and descending Sentinel-1A images were calculated.Almost half of these landslides were located in the low-effect area in the descending SAR data, while almost half of them were located in the shadow or high-effect area in the ascending SAR data.Thus, we selected 56 descending images of Sentinel-1A (2019-2021) as the only data source for the CSI. Figure 6c reveals that these landslides are distributed in the river valley; consequently, these landslides have developed well because of the river's perennial scouring and dense rainfall.

Results of CSI and HSRP
To illustrate the improved performance of the CSI, we employed SBAS-InSAR for the same study area based on the same SAR data source [14].Figure 7a,b shows the annual rates along the LOS direction, as measured using SBAS-InSAR and CSI, respectively.The CSI's deformation signal and total number of MPs are clearer and higher than those of the SBAS-InSAR, respectively.The maximum annual rate reaches −64 mm/yr within the monitoring period.
The monthly ~1 km HSRP datasets can be retrieved using the method in Section 3.3; Figure 7c shows an example of the HSRP in June 2020.On this basis, the HSRP and ASPMP of every landslide can extracted one by one.

Results of CSI and HSRP
To illustrate the improved performance of the CSI, we employed SBAS-InSAR for the same study area based on the same SAR data source [14].Figure 7a,b shows the annual rates along the LOS direction, as measured using SBAS-InSAR and CSI, respectively.The CSI's deformation signal and total number of MPs are clearer and higher than those of the SBAS-InSAR, respectively.The maximum annual rate reaches −64 mm/yr within the monitoring period.
The monthly ~1 km HSRP datasets can be retrieved using the method in Section 3.3; Figure 7c shows an example of the HSRP in June 2020.On this basis, the HSRP and ASPMP of every landslide can extracted one by one.

Analysis of a Typical Landslide
To demonstrate the applicability of our method to every landslide, we selected the largest landslide as a typical case.This landslide (102.876°E, 30.963°N) is located in Niela Village in Danba County, which is next to the Jinchuan River.The low vegetation coverage makes the interferograms have a high coherence.Meanwhile, the orientation and slope angle satisfy the observations of the descending Sentinel-1A data.
As can be seen from Figure 8, this landslide is about 1150 m long and 1200 m wide, which covers an area of up to 1,275,267 m 2 .CSI can obtain clearer deformation signals and more MPs than SBAS-InSAR, and the maximum annual deformation's LOS velocity reached 59 mm/yr.To further reveal the spatiotemporal evolution, the cumulative deformation curve of three points (P1, P2, and P3) from the slope's back edge to the tongue was combined with the monthly HSRP.To avoid any misunderstanding, we inverted the result's value to reflect the slip down.As shown in Figure 9, the cumulative displacement in the period of the monitoring reached almost 75 mm, and accelerated sliding is highly related to intense rainfall events between April and September each year (the temporal window for our ASPMP).Therefore, the gathered precipitation may have contributed to the formation of gullies on the landslide body.Water from the gullies infiltrates deposits composed of silty clay, gravelly soil, and cracks; consequently, a slip surface formed.With increasing water content and slope weight, the slip accelerates significantly.

Analysis of a Typical Landslide
To demonstrate the applicability of our method to every landslide, we selected the largest landslide as a typical case.This landslide (102.876• E, 30.963 • N) is located in Niela Village in Danba County, which is next to the Jinchuan River.The low vegetation coverage makes the interferograms have a high coherence.Meanwhile, the orientation and slope angle satisfy the observations of the descending Sentinel-1A data.
As can be seen from Figure 8, this landslide is about 1150 m long and 1200 m wide, which covers an area of up to 1,275,267 m 2 .CSI can obtain clearer deformation signals and more MPs than SBAS-InSAR, and the maximum annual deformation's LOS velocity reached 59 mm/yr.To further reveal the spatiotemporal evolution, the cumulative deformation curve of three points (P1, P2, and P3) from the slope's back edge to the tongue was combined with the monthly HSRP.To avoid any misunderstanding, we inverted the result's value to reflect the slip down.As shown in Figure 9, the cumulative displacement in the period of the monitoring reached almost 75 mm, and accelerated sliding is highly related to intense rainfall events between April and September each year (the temporal window for our ASPMP).Therefore, the gathered precipitation may have contributed to the formation of gullies on the landslide body.Water from the gullies infiltrates deposits composed of silty clay, gravelly soil, and cracks; consequently, a slip surface formed.With increasing water content and slope weight, the slip accelerates significantly.Moreover, we draw a profile map of the landslide and extract a profile line from A1 to A2 in the annual deformation map (Figure 10).As seen in Figures 9 and 10, the slope tongue's rate of and accumulated deformation are the highest, those of the slope's middle part are the second highest, and those of the slope's back edge are the lowest.According to our field investigation, fluvial erosion by the Dajinchuan River is another triggering factor.Especially in the flood season, the river water's level may increase by nearly 5 m with respect to its normal status, with the peak flux exceeding 650 m 3 /s.The hydrodynamic pressure variations induced by the rapid river-waterlevel's changes may lead to instability and movement at the front edge of the landslide [42].
In the first alignment design version, the recommended line's mainline (K) and alternative line (B) are tunneled through the slope's back edge and middle part, respectively.The annual movement rate of the geometric center point on the landslide is approximately 36 mm/yr, which is higher than "30 mm/yr", and the corresponding ASPMP (598.75 mm) is higher than ASPMPM (598.57mm).According to our method's rule, we chose tunnel and reinforce.And because construction along line B will intensify the deformation according to the spatiotemporal evolution, the final design version selects the mainline (K).Moreover, we draw a profile map of the landslide and extract a profile line from A1 to A2 in the annual deformation map (Figure 10).As seen in Figures 9 and 10, the slope tongue's rate of and accumulated deformation are the highest, those of the slope's middle part are the second highest, and those of the slope's back edge are the lowest.According to our field investigation, fluvial erosion by the Dajinchuan River is another triggering factor.Especially in the flood season, the river water's level may increase by nearly 5 m with respect to its normal status, with the peak flux exceeding 650 m 3 /s.The hydrodynamic pressure variations induced by the rapid river-water-level's changes may lead to instability and movement at the front edge of the landslide [42].

Landslide-Hazard-Avoiding Alignment Selection
From the above selection case, the 79 remaining landslides can be reconsidered legitimately for the first alignment design version.As described in Table 3, 33 landslides in total had a high impact (including the above case), 1 landslide had a moderate impact, and 3 landslides in total had a low impact (The other landslides far away from the line had no impact and, therefore, are not mentioned in the table.).The annual rates and ASPMPs of the landslides were extracted to obtain solutions for the alignment selection.As Table 3 shows, 13 type A, 16 type B, and 8 type C were selected for every landslide.And the corresponding final alignment design version was applied for the BKHC construction.

Landslide-Hazard-Avoiding Alignment Selection
From the above selection case, the 79 remaining landslides can be reconsidered legitimately for the first alignment design version.As described in Table 3, 33 landslides in total had a high impact (including the above case), 1 landslide had a moderate impact, and 3 landslides in total had a low impact (The other landslides far away from the line had no impact and, therefore, are not mentioned in the table .).The annual rates and ASPMPs of the landslides were extracted to obtain solutions for the alignment selection.As Table 3 shows, 13 type A, 16 type B, and 8 type C were selected for every landslide.And the corresponding final alignment design version was applied for the BKHC construction.

Advantages
Because of the tremendous costs caused by design alterations, a risk-optimal highway alignment selection is essential before the construction.Thus, monitoring and stability analyses of landslides are vital through the whole alignment-design workflow.Compared with traditional methods, the new systematic highway alignment selection approach we proposed combines the CSI and ~1 km monthly HSRP analysis.It can give a more accurate quantitative reference to help to finalize the alignment.
As the fundamental tool in landslide monitoring, TS-InSAR has to face spatiotemporal decorrelations, atmospheric delay anomalies, geometric distortions, etc., especially in the mountainous area in southwestern China [14][15][16][17][18][19][20].Consequently, SqueeSAR TM has been widely employed to solve a part of these problems, in some period of time, as a representative next-generation TS-InSAR technique.After reviewing the existing research about SqueeSAR TM , we found that it has evolved to series of algorithms [35][36][37][38][39][40][41][42], from which we selected the CSI as the key step of our method to provide a much denser observation network for performing more reliable phase unwrapping to obtain more accurate MPs.
However, no matter which algorithm is utilized, geometric distortion still restricts the accuracy and effectiveness for monitoring from the root.In most studies or actual projects, both the ascending and descending SAR images are used without any FS analysis.The tremendous data redundancy reduces the TS-InSAR's efficiency, especially in wide ranges.Therefore, in our method, data selection is regarded as being the first key step based on FS calculations, which can reduce the data downloading and processing times.
Rainfall is the most important factor inducing landslides, especially in the rainy season.After obtaining enough deformation signals and accurate rates of landslides, the accuracy of the precipitation data should be improved as well.In the mountainous area of southwestern China, high-spatial-resolution precipitation data cannot be obtained because of the low station density, even with the establishment of additional stations in recent years [53][54][55][56][57][58][59].On reviewing the related research, we found that the station density decides the result accuracy of conventional interpolation methods, including inverse distance weighting, kriging methods, and regression analysis.Thus, delta spatial downscaling is used to obtain the ~1 km HSRP at the local scale based on CRU time-series datasets and WorldClim climatology datasets in our method.Compared with the data from the climate stations far from the BKHC, the HSRP exhibits smaller deviations and is more accurate [60][61][62][63][64][65][66].Moreover, our method extracted the ASPMP of every landslide and ASPMPM based on the HSRP, which would result in a reasonable reference index for the alignment selection.
Herein, our method incorporated the advantages of these new techniques, attempted to combine all the accurate results, and considered the locations and deformation rates of landslides and ASPMPs as the major factors, while construction possibly impacts the secondary factors.On these bases, the risk-optimal highway alignment can be finalized reasonably and efficiently based on our three solution types.

Limitations and Further Directions
Regardless of the TS-InSAR method's advantages that we incorporated, the computational expense of the CSI is still unacceptable for a normal computer.In our case study, we used a computer with an i9-11900 CPU, a 128 GB RAM, and a 4 TB hard disk to conduct the experiments.Hence, how to improve the computational efficiency and reduce the consumption of the storage space have become critical problems to be solved.Moreover, L-band images have more penetration and X-band images have a higher resolution than C-band images, respectively.We only processed C-band SAR images, but the other two types of band images were not tested.It is expected that our method can not only benefit from them but also face computing and FS-analyzing challenges.In addition, the Generic Atmospheric Correction Online Service's (GACOS's) atmospheric delay products were not used in our method to improve the accuracy of the monitoring, which can provide a new approach for the atmospheric correction of repeat-pass InSAR [75].
As for the HSRP analysis, the data-downscaling method was employed for the 30 ′ CRU and WorldClim time-series data in our method after reviewing previous studies [65,66].However, on one hand, the averaged 30 ′ elevation was used for the original CRU data, which weakened the representation of the precipitation on the actual land surface, especially in the mountainous area.On the other hand, the weather stations used for the CRU evaluation are located in valleys near counties or cities, which cannot reflect the real conditions of wilderness areas.Meanwhile, as the reference climatology dataset, the WorldClim dataset performed well.But a new and better reference climatology dataset needs to be generated for southwestern China, and the collection of public and private climate data can be a good solution.
There are many other environmental factors that affect landslide stability, such as earthquakes, reservoir impoundments, fluvial erosion, and weathering.But in the final step of our method, we only selected the ASPMP as the dominant environmental factor to help to select the alignment.In practice, other factors should be considered comprehensively to establish a judgement model.For instance, landslides are typical post-disasters of earthquakes, especially in the eastern margin of the Qinghai-Tibetan Plateau.This region is a very active seismic area in which three major earthquakes have been recorded since the 1930s, namely, the 1933 Mw 7.3 Diexi earthquake, the 1976 Mw 6.7 Songpan-Pingwu earthquake, and the 2008 Mw 7.9 Wenchuan earthquake.During these strong seismic activities, the rock mass was damaged by the discordant deformation from different existing lithologies.Hence, the locations of earthquake zones should be considered synthetically for better highway alignment selections in future work.And the possible impacts of construction on landslides need a quantitative analysis rather than just a qualitative deduction.
In the future, our research emphasis will focus on parallel computing, automatic processing based on deep-learning models, more landslide-event-training samples, multiplatform SAR images, higher-spatial-resolution DEMs, etc.

Conclusions
In this study, a section of the BKHC was taken as the study area, which has a severe alignment design problem due to the complicated geological background, strong tectonic movement, rock body rupture, as well as frequent geohazards.Hence, this paper proposed a new systematic highway alignment selection method for a typical section of the BKHC, which has 80 slow-moving landslides identified by optical remote-sensing interpretation and field investigation, based on CSI and ~1 km monthly HSRP analysis.
Specifically, the FSs of Sentinel-1A ascending and descending images were first calculated to select suitable SAR datasets according to the locations of identified landslides.Consequently, 56 Sentinel-1A descending images from September 2019 to September 2021 were processed using CSI and SBAS-InSAR to obtain the time-series deformation and rates of landslides and illustrate the strengths of the CSI.The results reveal that CSI has clearer deformation signals and more MPs than SBAS-InSAR.And the maximum cumulative displacements and rates of the landslides reach −75 mm and −64 mm/yr within the monitoring period (CSI results), respectively.To demonstrate the applicability of our method to every landslide, we selected the largest landslide as a typical case.Through further spatiotemporal evolution analysis of this case, we found that accelerated sliding is highly related to intense rainfall events and that fluvial erosion impacts the slope's stability perennially.Under this circumstance, we utilized delta downscaling to obtain the ~1 km HSRP based on CRU and WorldClim time-series data, and the ASPMPs of every landslide were extracted afterward.Subsequently, three solution types were established in our method by setting thresholds for deformation rates and ASPMPs.Finally, the risk-optimal alignment selection of the BKHC was reasonably and efficiently finalized based on these solution types and the construction's possible impacts.
This proposed method expands the application range of the TS-InSAR in actual transport engineering and reduces the time and economic and labor costs associated with recurrent landslides.Though some limitations need to be resolved in the future, these research results offer a valuable consultation for the alignment selection of highways in mountainous areas and, consequently, for landslide management.

Figure 2 .
Figure 2. Faults in the study area.

Figure 2 .
Figure 2. Faults in the study area.

Figure 3 .
Figure 3. Distribution of national weather stations across China (modified from figure 1 from in [62].

Figure 3 .
Figure 3. Distribution of national weather stations across China (modified from Figure 1 from in [62]).

Figure 5 .
Figure 5. Relationship between the LOS and the downslope displacements for different slope orientations (scenarios 1, 2, and 3 are facing the sensor and scenarios 4, 5, and 6 are facing away from the sensor) and slopes (α and α') (adapted with permission from Ref. [67].2016, Keren Dai).

Figure 6 .
Figure 6.The FSs of the Sentinel-1A data: (a) ascending and (b) descending; (c) the terrain of the study area.

Figure 6 .
Figure 6.The FSs of the Sentinel-1A data: (a) ascending and (b) descending; (c) the terrain of the study area.

Figure 8 .
Figure 8.The annual deformation rates derived from descending Sentinel-1a data using two TS-InSAR methods: (a) SBAS-InSAR and (b) CSI; (c) is the optical remote-sensing interpretation; (d) is the field investigation.

Figure 9 .
Figure 9.The cumulative displacement curves of typical points and monthly HSRPs of the landslide.

Figure 10 .
Figure 10.Profile map of the landslide and the profile line along A1-A2.

Figure 10 .
Figure 10.Profile map of the landslide and the profile line along A1-A2.

Table 2 .
Feasibility of SAR data.

Table 3 .
The impacts of the landslides on the study area and the corresponding advice for the alignment design.