Wide Area Detection and Distribution Characteristics of Landslides along Sichuan Expressways

: Wide area landslide detection is a major international research hotspot in the ﬁeld of geological hazards, and the integration of multi-temporal optical satellite images and spaceborne interferometric synthetic aperture radar (InSAR) appears to be an effective way to realize this. In this paper, a technical framework is presented for wide area landslide detection: (i) multi-temporal satellite optical images are used to detect landslides with distinguishable geomorphological features; (ii) Generic Atmospheric Correction Online Service (GACOS) assisted InSAR stacking is employed to generate annual surface displacement rate maps in radar line of sight using satellite SAR images from both ascending and descending tracks, which are in turn utilized to automatically detect active landslides from ground motion using hotspot analysis, and (iii) the distribution characteristics of the detected landslides are investigated by examining their relationships with topographic and hydrological factors. Three expressways in Sichuan Province, China—namely the Yakang (Ya’an-Kangding), Yaxi (Ya’an-Xichang), and Lushi (Luding-Shimian) expressways—and their surrounding regions (a total area of approximately 20,000 square kilometers) were chosen as the study area. A total of 413 landslides were detected, among which 320 were detected using multi-temporal satellite optical images, and 109 were detected using GACOS-assisted InSAR stacking. It should be noted that only 16 landslides were detected by both approaches; these landslides all exhibited not only obvious geomorphological features but also ground motion. A statistical analysis of the topographic and hydrological factors shows that of the detected landslides: 81% are distributed at elevations of 1000–2500 m, over 60% lie within the elevation range of 100~400 m, and 90% present with medium and steep slopes (20 ◦ ~45 ◦ ), and 80% are located within areas seeing an annual rainfall of 950~1050 mm. Nine landslides were found to pose potential safety hazards to the expressways. The research ﬁndings in this paper have directly beneﬁtted the Sichuan expressways; equally important, it is believed that the technical framework presented in this paper will provide guidance for hazard mitigation and the prevention of transportation hazards in the future.

Geophysical prospecting technology is one non-destructive geophysical method that can be used to detect landslides, which utilizes both natural and artificially created physical fields, (e.g., gravitational, magnetic, electrical, seismic, and radiometric fields) acquired on the surface of the Earth, in the air and/or underground (in wells and shafts) to determine the location of geological structures, ore bodies, and so forth and their fundamental characteristics [21,22]. However, the method has certain errors and is costly to implement and is highly affected by geological conditions. The transition zone from the Qinghai-Tibet Plateau to the Sichuan Basin is characterized by clouds and fogs, dense vegetation, high elevation differences, and strong water vapor variations. Therefore, wide area landslide detection in this region is challenging. In this paper, a technical framework is presented for wide-area landslide detection. First, multi-temporal satellite optical images are used to detect landslides with distinguishable geomorphological features and qualitatively assess the threat to highways. Secondly, Sentinel-1 images from both ascending and descending tracks are interferometrically processed and then analyzed with GACOS-assisted InSAR stacking to generate annual surface displacement rate maps in the radar line of sight (LOS). Deforming landslides are automatically extracted by the hotspot analysis method. Finally, the main conditions that control and influence the development and distribution of landslides are analyzed through the statistics of features of all the detected landslides, such as elevation, elevation difference, area, slope angle, slope aspect, rainfall, and other characteristics. The results will be directly applied to landslide hazard mitigation and prevention of hazards along Sichuan expressways.

Study Area
The research area of this study encompasses approximately 20,000 km 2 along and around the Yakang, Yaxi, and Lushi expressways in Sichuan Province, China (Figure 1), located in the transitional zone between the Sichuan Basin and the Qinghai-Tibet Plateau. The area comprises a large topographic terrain, complex geomorphological types, and complicated geological conditions. Land cover is dominated by forests, grasslands, and shrubs ( Figure 2a). Because of the strong influence of rainfall, earthquakes, and human engineering activity, landslides often occur along and near the expressways in the study area. For instance, at 03:50 on 21 August 2020, the Zhonghaicun landslide occurred due to heavy rainfall, and road traffic was interrupted for a long period [23], which seriously threatened the construction and safe operation of the expressways. Therefore, it is necessary to carry out wide area landslide detection along and around the expressways in the region as a matter of urgency.

Data
Multi-source satellite optical images and Sentinel-1 radar images were used in this study, as shown in Table 1. The selected satellite optical images included Gaofen-2 satellite images and Google Earth archived images. The Gaofen-2 satellite was successfully launched on 19 August 2014. The spatial resolution of the panchromatic band is 0.8 m, the multi-spectral resolution is 3.2 m, the width is 45 km, and the revisiting time is five days, meeting the requirements of landslide detection. The archived images of Google Earth are mainly satellite images with a spatial resolution below 1 m, such as WorldView, Quickbird, and GeoEye, which provide abundant data resources for landslide detection. The European Space Agency's Sentinel-1A/B satellites operate day and night, performing C-band microwave SAR imaging and providing radar imagery with wide coverage, (e.g., 250 km × 250 km) and a short repeat cycle (12 days for a single satellite and six days for both). In this study, a total of 521 Sentinel-1A/B images in four tracks (Table 1) were selected for the single-look complex (SLC) images of VV polarization in the interferometric wide swath (IW).       [24]. In the case of landslide identification, Li et al. (2022) proposed to divide landslides into four types in terms of their histories and movement: actively deforming slopes (I), reactivated historically deformed slopes (II), stabilized historically deformed slopes (III), and undeformed but potentially unstable slopes (IV) [25] (Figure 2). The actively deforming slopes (I) refer to the slopes without historical failures which are actively deforming but there has been no historical slope failure; the reactivated historically deformed slopes (II) refer to the slopes with historical failures that have been reactivated, (i.e., in motion); the stabilized historically deformed slopes (III) refer to the slopes with historical failures that are in a stable state, and the undeformed but potentially unstable slopes (IV) refer to the slopes without historical failures that are in a stable state but with a slide-prone structure [25]. Types I and II can often be identified through mapping their deformation signals since they are in motion, and geomorphologic features are usually utilized to identify Types II and III since they exhibit historical failures. In contrast, it is difficult (if not impossible) to identify Type IV since there is neither a deformation signal nor a geomorphologic feature with such slopes.

Multi-Temporal Satellite Optical Image Interpretation
The multi-temporal satellite optical image interpretation method has the characteristics of wide range, non-contact, low cost, and high precision, and is suitable for wide area landslide detection [18]. Utilizing spectral and textural differences in satellite optical images acquired at different times, combined with topographic and geomorphologic features, it is straightforward to identify actively deforming slopes with obvious deformation signals ( Figure 3a) and historically deformed slopes (Figure 3b,c) [10]. The actively deforming slopes with obvious deformation signals can exhibit cracks at tops (heads) of the slope and small-scale collapses on the lower parts of the slope [2,26]. Because of historical failures, the historically deformed slopes often exhibit certain geomorphological features, such as scarps, flanks, cracks, ridges, etc. [25,26]. Hence, such geomorphological features can be used to identify the abovementioned landslides. Note that satellite optical images can be seriously affected by the presence of clouds, fogs, and their corresponding shadows on the ground, which in turn might result in false identification of landslides. To address this, multi-temporal satellite optical images are often employed to ensure there are sufficient cloud-free images to examine the evolution of the slopes. failures, the historically deformed slopes often exhibit certain geomorphological features, such as scarps, flanks, cracks, ridges, etc. [25,26]. Hence, such geomorphological features can be used to identify the abovementioned landslides. Note that satellite optical images can be seriously affected by the presence of clouds, fogs, and their corresponding shadows on the ground, which in turn might result in false identification of landslides. To address this, multi-temporal satellite optical images are often employed to ensure there are sufficient cloud-free images to examine the evolution of the slopes.

InSAR Automatic Landslide Detection
InSAR technology can be employed to generate information on surface elevation and surface deformation from the interference signals of two or more SAR images covering the same area [27,28]. When the microwave signal from the SAR satellite sensor passes through the atmosphere, the propagation path will be affected to varying degrees with time. Any differences in the atmosphere at the time of SAR image acquisition during different time periods will directly affect the InSAR observations. This influence is usually known as the atmospheric effect, which is also one of the most significant error sources in InSAR deformation monitoring [29,30]. Most landslides in the study area are located in mountainous areas, where topography-related tropospheric delays are particularly significant, and atmospheric disturbances are not negligible where there are rivers, lakes, and other bodies of water nearby, or strong convective weather is present during satellite transits. Therefore, prior to landslide detection, it is necessary to carry out a reasonable and effective atmospheric correction on interference images. The GACOS (http://www.gacos. net/ (accessed on 25 October 2021)) system developed by Prof. Zhenhong Li's team is based on HERS ECMWF data with a spatial resolution of 0.125 • and a temporal resolution of six hours. It uses the Iterative Tropospheric Decomposition (ITD) model to separate the stratification and turbulence signals from the total tropospheric delay, and can provide near-real-time, global coverage of InSAR atmospheric delay error correction mapping for atmospheric delay error correction in InSAR measurements [31][32][33].
Considering the large topographic relief in the study area, the presence of strong convective weather, and the obvious atmospheric influence on InSAR interferometric images, the GACOS-assisted InSAR stacking method is an efficient means to generate annual surface displacement rate maps [34,35]. In this study, a multi-looking factor of 4 in range, (i.e., a resolution of approximately 17 m in the ground direction) and 1 in azimuth, (i.e., a resolution of approximately 14 m in the ground direction) were used to process the Sentinel-1 SAR images. The SAR images were first pre-processed to obtain single-looking complex SLC images, and the topographic influence was attenuated using the 30 m resolution ALOS Global Digital Surface Model (AW3D30) released by the Japan Aerospace Exploration Agency (JAXA). A time baseline of 60 days and a spatial baseline of 150 m were set to generate interferograms, and a series of unwrapping maps were obtained by adaptive filtering and phase unwrapping processing, after which atmospheric correction products of GACOS were introduced to attenuate atmospheric effects ( Figure 4). It should be pointed out that the average coherence greater than 0.4 can be observed in 80% of the study area, suggesting the feasibility of InSAR for this area (Figure 2b).  The InSAR stacking technique is based on the assumption that surface deformation varies approximately linearly and the average rate per pixel is estimated from N interferometric images [36]. This has the advantages of low computational resource requirements, straightforward automated processing, simple, effective, and fast application, and the ef- The InSAR stacking technique is based on the assumption that surface deformation varies approximately linearly and the average rate per pixel is estimated from N interferometric images [36]. This has the advantages of low computational resource requirements, straightforward automated processing, simple, effective, and fast application, and the effective suppression of random errors. The least square method was used to solve the average unwrapping rate following GACOS correction, the interference phases were weighted by time intervals and, finally, the annual surface deformation rate map of LOS direction in the study area was generated [35], as in Equation (1).
where V mean is the mean velocity, ϕ i demonstrates the unwrapped phase and ∆T i is the temporal baseline of the ith interferogram.
The Getis-Ord G * i statistic (G * i )-also known as hot-spot analysis [37,38]-is a method commonly used to analyze location-related tendencies (clustering) among the attributes of spatial data (points or areas). The method is an adaptation of the General G-statistic [37], a global method for quantifying the degree of spatial autocorrelation over an area. It is first necessary to determine whether a feature point has significant spatial clustering characteristics by looking at each pixel in the environment of adjacent elements within a certain distance. Pixels with significant spatial clustering characteristics have high value attributes and are surrounded by other elements that also have high values [33,34]. It is time-consuming and laborious work to visually interpret the annual average surface deformation rates of a large area in order to ascertain locations of landslide risk. Therefore, the Getis-Ord hotspot analysis method was introduced to automatically extract effective deforming areas from the annual surface deformation rate map. The calculation formula is as follows: is computed for feature (i) at a distance (d) standardized as a z-score. x j is the annual surface deformation value of each neighbor pixel. ω ij is the spatial weight for the target neighbor i and j pair and n is the total number of pixels. The Euclidean method was used to calculate the geographical distances from each feature to its neighboring features.
The degree of clustering of deformation points is assessed based on the z-score and p-Value. The higher the z value, the greater the likelihood of clustering at a particular location. In this study, deformation points satisfying p < 0.01 (99% confidence level) and z > ±|2.58| were selected as having significant spatial clustering properties [39].

Landslide Identification with Integrated Remote Sensing Techniques
After obtaining geomorphological and deformational information as demonstrated in the previous two sections, landslides can be identified (Table 2). Type I landslides (actively deforming slopes) are detected mainly on the basis of deformation signals; since some of them may exhibit geomorphological features such as cracks, scarps and collapses, multitemporal satellite optical image interpretation may provide assistance in some cases. Type II landslides (reactivated historically deformed slopes) should be detected on the basis of deformation signals and geomorphological features, whilst Type III landslides (stabilized historically deformed slopes) should be detected by interpreting geomorphological features with optical remote sensing images. Type IV landslides (undeformed but potentially unstable slopes) have no deformation signal or geomorphological feature, and neither InSAR nor optical remote sensing imagery can be utilized to detect Type IV landslides. Such landslides can be identified according to their geological structures, and geophysical prospecting technology should be used for this purpose [21,22]. Geophysical prospecting technology is costly to implement, and its performance is highly affected by geological conditions. Hence, it is still difficult to detect Type IV landslides over wide regions at present.

Annual Surface Deformation Rate
The GACOS-assisted InSAR stacking method was used to process 521 SAR images acquired during the period from 2017 to 2021 to generate LOS annual deformation rate maps from both ascending and descending tracks encompassing appropriately 40,000 km 2 of three expressways in Sichuan, as shown in Figure 5a,b, respectively. The annual surface deformation rate maps show that most of the landslides threatening the expressway are concentrated in the Luding-Shimian, Shimian-Hanyuan and Hanyuan-Yingjing sections, as shown in the three red areas in Figure 5, with the landslides in the area 1 ( Figure 6a The large capacity of the Pubugou Hydropower Station is likely to be a particularly significant factor in the landslides in Zone 1. There are 16 landslides in Zones 1, 2 and 3 with an annual surface deformation rate greater than 20 mm/a, representing great safety hazards to the expressway ( Figure 6). From Figure 6a,d, it can be seen that more landslides are detectable in the ascending track images than in the descending track images in Zone 1, whereas landslides A and H are detectable in both the ascending and descending track images. However, landslides B, C, D, E, F, G and I are only detectable in the ascending track images and landslide J is only detectable in the descending track images. The situation in Zones 2 and 3 is similar to that in Zone 1, where there are landslides that could be detected by both ascending and descending track images and by descending or ascending track images only.

Annual Surface Deformation Rate
The GACOS-assisted InSAR stacking method was used to process 521 SAR images acquired during the period from 2017 to 2021 to generate LOS annual deformation rate maps from both ascending and descending tracks encompassing appropriately 40,000 km 2 of three expressways in Sichuan, as shown in Figure 5a,b, respectively. The annual surface deformation rate maps show that most of the landslides threatening the expressway are concentrated in the Luding-Shimian, Shimian-Hanyuan and Hanyuan-Yingjing sections, as shown in the three red areas in Figure 5, with the landslides in the area 1 ( Figure 6a The large capacity of the Pubugou Hydropower Station is likely to be a particularly significant factor in the landslides in Zone 1. There are 16 landslides in Zones 1, 2 and 3 with an annual surface deformation rate greater than 20 mm/a, representing great safety hazards to the expressway ( Figure 6). From Figure 6a,d, it can be seen that more landslides are detectable in the ascending track images than in the descending track images in Zone 1, whereas landslides A and H are detectable in both the ascending and descending track images. However, landslides B, C, D, E, F, G and I are only detectable in the ascending track images and landslide J is only detectable in the descending track images. The situation in Zones 2 and 3 is similar to that in Zone 1, where there are landslides that could be detected by both ascending and descending track images and by descending or ascending track images only.

Annual Surface Deformation Rate
The GACOS-assisted InSAR stacking method was used to process 521 SAR images acquired during the period from 2017 to 2021 to generate LOS annual deformation rate maps from both ascending and descending tracks encompassing appropriately 40,000 km 2 of three expressways in Sichuan, as shown in Figure 5a,b, respectively. The annual surface deformation rate maps show that most of the landslides threatening the expressway are concentrated in the Luding-Shimian, Shimian-Hanyuan and Hanyuan-Yingjing sections, as shown in the three red areas in Figure 5, with the landslides in the area 1 (Figure 6a The large capacity of the Pubugou Hydropower Station is likely to be a particularly significant factor in the landslides in Zone 1. There are 16 landslides in Zones 1, 2 and 3 with an annual surface deformation rate greater than 20 mm/a, representing great safety hazards to the expressway ( Figure 6). From Figure 6a,d, it can be seen that more landslides are detectable in the ascending track images than in the descending track images in Zone 1, whereas landslides A and H are detectable in both the ascending and descending track images. However, landslides B, C, D, E, F, G and I are only detectable in the ascending track images and landslide J is only detectable in the descending track images. The situation in Zones 2 and 3 is similar to that in Zone 1, where there are landslides that could be detected by both ascending and descending track images and by descending or ascending track images only.

Annual Surface Deformation Rate
The GACOS-assisted InSAR stacking method was used to process 521 SAR images acquired during the period from 2017 to 2021 to generate LOS annual deformation rate maps from both ascending and descending tracks encompassing appropriately 40,000 km 2 of three expressways in Sichuan, as shown in Figure 5a,b, respectively. The annual surface deformation rate maps show that most of the landslides threatening the expressway are concentrated in the Luding-Shimian, Shimian-Hanyuan and Hanyuan-Yingjing sections, as shown in the three red areas in Figure 5, with the landslides in the area 1 (Figure 6a The large capacity of the Pubugou Hydropower Station is likely to be a particularly significant factor in the landslides in Zone 1. There are 16 landslides in Zones 1, 2 and 3 with an annual surface deformation rate greater than 20 mm/a, representing great safety hazards to the expressway ( Figure 6). From Figure 6a,d, it can be seen that more landslides are detectable in the ascending track images than in the descending track images in Zone 1, whereas landslides A and H are detectable in both the ascending and descending track images. However, landslides B, C, D, E, F, G and I are only detectable in the ascending track images and landslide J is only detectable in the descending track images. The situation in Zones 2 and 3 is similar to that in Zone 1, where there are landslides that could be detected by both ascending and descending track images and by descending or ascending track images only.
No deformation signal or geomorphological feature Geophysical prospecting technology

Annual Surface Deformation Rate
The GACOS-assisted InSAR stacking method was used to process 521 SAR images acquired during the period from 2017 to 2021 to generate LOS annual deformation rate maps from both ascending and descending tracks encompassing appropriately 40,000 km 2 of three expressways in Sichuan, as shown in Figure 5a,b, respectively. The annual surface deformation rate maps show that most of the landslides threatening the expressway are concentrated in the Luding-Shimian, Shimian-Hanyuan and Hanyuan-Yingjing sections, as shown in the three red areas in Figure 5, with the landslides in the area 1 (Figure 6a The large capacity of the Pubugou Hydropower Station is likely to be a particularly significant factor in the landslides in Zone 1. There are 16 landslides in Zones 1, 2 and 3 with an annual surface deformation rate greater than 20 mm/a, representing great safety hazards to the expressway ( Figure 6). From Figure 6a,d, it can be seen that more landslides are detectable in the ascending track images than in the descending track images in Zone 1, whereas landslides A and H are detectable in both the ascending and descending track images. However, landslides B, C, D, E, F, G and I are only detectable in the ascending track images and landslide J is only detectable in the descending track images. The situation in Zones 2 and 3 is similar to that in Zone 1, where there are landslides that could be detected by both ascending and descending track images and by descending or ascending track images only.

Landslide Detection and Mapping
The combined multi-temporal satellite optical image visual interpretation and GACOSassisted InSAR stacking methods detected a total of 413 landslides. As shown in Figure 7a, 320 landslides were detected through satellite optical image visual interpretation. A total of 109 landslides were automatically detected by the hotspot analysis method based on InSAR-derived LOS annual surface displacement rate maps, and the hotspot analysis results for Zones 1 and 2 are shown in Figure 8. It is worth noting that only 16 landslides with obvious deformation signals and geomorphological features were detected by both approaches, mainly because the landslide types that can be detected by the two methods are different [40]. Then, according to the classification of landslides in this study, the 413 detected landslides were divided into 371 actively deforming slopes, 23 reactivated historically deformed slopes and 19 stabilized historically deformed slopes (Figure 7b).
Where the landslides that can be detected using high-resolution satellite optical images include landslides that have occurred and landslides that are incubating with obvious deformation signals, and since the spatial resolution of the satellite optical images used in this paper was 0.8 m, landslides without geomorphological features cannot be detected by this method. In contrast, the GACOS-assisted InSAR stacking method can only be used to detect landslides in motion, (i.e., Types I and II). Although the two approaches can detect different types of landslides, they are nevertheless both landslides and, therefore, when it comes to landslide detection, it is necessary to combine the two technologies in order that they may complement and verify each other.    Landslide cataloging must determine not only the location of any landslide but also its extent. Both the geomorphological features of the satellite optical images and the annual surface deformation rate generated by the GACOS-assisted InSAR stacking method can be used as a basis for mapping the extent of landslides, as shown in Figure 9. Four cases are included in this research, combining the two methods. In case 1, landslides can be detected by both approaches, as shown in Figure 9a, and mapping the extent of landslides is based on the range of annual surface deformation rate map, and geomorphological features taken from satellite optical images. In case 2, landslides can only be identified by multi-temporal satellite optical image interpretation, as shown in Figure 9b, and the extent of the mapped landslide is determined based on high-resolution satellite images over multiple periods. In case 3, landslides can only be detected by InSAR, as shown in Figure 9c, and mapping the extent of landslides requires a combination of deformation extent, multi-temporal satellite optical images, DEM, slope and other multi-source data. In case 4, as shown in Figure 9d, neither of these two methods can detect landslides; this is also the main area in which landslides are missed due to the extent of dense vegetation. In these areas, the geographical features are affected by vegetation, meaning that satellite optical images cannot be interpreted, and InSAR results are incoherent or geometrically distorted. In such areas, airborne LiDAR is believed to be an effective tool to detect old landslides covered by dense vegetation [40].  (ii) LOS annual surface deformation rate maps from ascending tracks, and (iii) LOS annual surface deformation rate maps from descending tracks. Solid red line means landslide is detectable and dashed red line means landslide is not detectable.

Inventory of Landslides That Pose a High Risk to Expressways
By scrutinizing the multi-temporal satellite optical images, annual surface displacement rate map and distance from the expressway, it was concluded that nine of the 413 landslides detected were potential threats to the expressway, of which three, four and two were safety hazards to Lushi, Yaxi and Yakang expressways, respectively. The length, width, area, distance from the expressway, satellite optical image texture characteristics, and annual surface displacement rate of the landslides were combined to qualitatively evaluate these nine landslides, as detailed in Table 3. (ii) LOS annual surface deformation rate maps from ascending tracks, and (iii) LOS annual surface deformation rate maps from descending tracks. Solid red line means landslide is detectable and dashed red line means landslide is not detectable.

Inventory of Landslides That Pose a High Risk to Expressways
By scrutinizing the multi-temporal satellite optical images, annual surface displacement rate map and distance from the expressway, it was concluded that nine of the 413 landslides detected were potential threats to the expressway, of which three, four and two were safety hazards to Lushi, Yaxi and Yakang expressways, respectively. The length, width, area, distance from the expressway, satellite optical image texture characteristics, and annual surface displacement rate of the landslides were combined to qualitatively evaluate these nine landslides, as detailed in Table 3. This is a deforming landslide, which was locally deformed and eroded by the river; it has the possibility of overall instability.

Effects of Spaceborne SAR Imaging Geometry
Due to the inherent oblique observation geometry of all imaging radar systems, surface slopes and similar terrain features lead to geometric distortions of the data acquired by SAR systems [12], the most relevant of which are layovers and shadows. These phenomena are particularly evident in high mountain valley areas and will result in parts of the detection area being left without a valid signal, ultimately leading to missed landslide detection [41]. In order to better compensate for the geometric distortions caused by SAR imaging geometry, this study combined the ascending and descending tracks of SAR images to obtain the deformation characteristics of detection targets under different imaging geometries, thus effectively reducing the areas that would not be effectively observed due to geometric distortions.
Quantitatively evaluating the combined ascending and descending tracks of SAR images can effectively reduce geometric distortions. In this paper, the geometric distortion and effective observable area were counted, as shown in Figure 10. The effective observation area made up 89.4% of the total area of the ascending track, the layover area represented 10.3% and the shadow area was less than 0.3%. The observable area in the descending track comprised 83.7% of the total area, the layover area made up 15.8% and the shadow area was just 0.5%. In the combined ascending and descending imagery, up to 97.3% of the region could be observed effectively, and the geometric distortion region was reduced to 2.7%. This result was consistent with the analysis that combined observation of ascending and descending tracks along the Sichuan-Tibet Railway, effectively reducing the single-track geometric distortion area by 30% to 1.5% [16]. These results show that the area of geometric distortion can be effectively reduced by using the joint observation of the ascending and descending tracks, whereupon the rate of missed detections of landslides can be reduced.

Identifying Deformation without Landslides
The area of deformation in the LOS annual surface deformation rate map includes not only landslides that are deforming but also alluvial fans formed when tributaries enter the main river channel, subsidence caused by agricultural irrigation, urban ground subsidence, and deformation caused by human engineering activities, as shown in Figure 11a-d, respectively. Certain cases, such as Zone I in Figure 8c and Zones I, S, and T in Figure 8d, will also be considered valid deformation areas when hotspot analysis is performed on the LOS annual surface deformation rate map. To effectively detect landslides, it is necessary to combine DEM and satellite optical images for comparison and validation, to exclude any misclassification of automatic landslide detection by hotspot analysis. the main river channel, subsidence caused by agricultural irrigation, urban ground subsidence, and deformation caused by human engineering activities, as shown in Figure  11a-d, respectively. Certain cases, such as Zone I in Figure 8c and Zones I, S, and T in Figure 8d, will also be considered valid deformation areas when hotspot analysis is performed on the LOS annual surface deformation rate map. To effectively detect landslides, it is necessary to combine DEM and satellite optical images for comparison and validation, to exclude any misclassification of automatic landslide detection by hotspot analysis.

Development Characteristics of Landslides along Sichuan Expressways
According to the landslide mapping method proposed in this paper, the surface area of the 413 landslides detected amounted to 47.32 km 2 (the total area of the study area is about 20,000 km 2 ), with minimum and maximum areas of 1176 m 2 and 6.980 km 2 , respectively.

Distribution Pattern of Detected Landslides with Respect to Topographic Factors
Rugged mountain terrain is more prone to landslides than flat terrain [42]. In previous studies, topographic parameters, such as elevation, relative altitude, slope angle, and aspect, were considered essential factors for the formation, development, and occurrence of landslides [42], and were widely used to evaluate the distribution of landslides [42,43]. The AW3D30 data at 30 m spatial resolution was employed to generate topographic factors, such as elevation, slope angle, slope aspect, and relative altitude, as shown in Figure  12.  Spatial relationships between landslide distribution and slope angle, slope aspect, precipitation and excess precipitation: (a) spatial relationship between slope angle and detected landslides; (b) spatial relationship between slope aspect and detected landslides; (c) spatial relationship between precipitation and detected landslides, and (d) spatial relationship between excess precipitation and detected landslides.
Elevation plays an important role in landslide distribution because it is related to land cover, roads and river valleys, slope material, and the seismic wave amplification effect [42]. From Figures 7 and 13a, it can be seen that approximately 80.5% of landslides occur between 1000 and 2500 m above sea level, and the maximum landslide density is distributed between 1000 and 2000 m. Below 1500 m, the number of landslides increases with elevation, while above 1500 m, the number decreases with elevation. Because slopes above 1500 m are less affected by rainfall and river erosion, the number and development of landslides are relatively low. As shown in Figure 13f, more than 77% of the landslides are located in the area with a relative altitude below 400 m, and the densest area of landslides is located in the area with relative altitudes of 100-300 m. When the relative altitude is below 200 m, the number of landslides increases with increased relative altitude. However, when the relative altitude is higher than 200 m, the number of landslides decreases with increased relative altitude. Figure 12. Spatial relationships between landslide distribution and slope angle, slope aspect, precipitation and excess precipitation: (a) spatial relationship between slope angle and detected landslides; (b) spatial relationship between slope aspect and detected landslides; (c) spatial relationship between precipitation and detected landslides, and (d) spatial relationship between excess precipitation and detected landslides.
The slope aspect may be the only topographic parameter with orientation in common topographic distribution analysis, and so it is usually combined with other directional parameters for further examination [42]. Previous studies have shown that the slope aspect influences weathering conditions, temperature, precipitation, deposition conditions, and land cover, and so plays a role in the susceptibility of an area to landslides [45,46]. As shown in Figure 12b, the slope aspect was divided into the following 16 Figure 13. Relationships between distribution of detected landslides and elevation, slope angle, slope aspect, precipitation, excess precipitation, and relative altitude, respectively: (a) relationship between distribution of detected landslides and elevation; (b) relationship between distribution of detected landslides and slope angle; (c) relationship between distribution of detected landslides and slope aspect; (d) relationship between distribution of detected landslides and precipitation; (e) relationship between distribution of detected landslides and excess precipitation, and (f) relationship between distribution of detected landslides and relative altitude.
The slope aspect may be the only topographic parameter with orientation in common topographic distribution analysis, and so it is usually combined with other directional parameters for further examination [42]. Previous studies have shown that the slope aspect influences weathering conditions, temperature, precipitation, deposition conditions, and land cover, and so plays a role in the susceptibility of an area to landslides [45,46]. As shown in Figure 12b Figure 13c, the slope aspects ENE, E, ESE, SE, SSE, SSW, SW, WSW, and W together account for 82.3% of the Figure 13. Relationships between distribution of detected landslides and elevation, slope angle, slope aspect, precipitation, excess precipitation, and relative altitude, respectively: (a) relationship between distribution of detected landslides and elevation; (b) relationship between distribution of detected landslides and slope angle; (c) relationship between distribution of detected landslides and slope aspect; (d) relationship between distribution of detected landslides and precipitation; (e) relationship between distribution of detected landslides and excess precipitation, and (f) relationship between distribution of detected landslides and relative altitude.

Distribution Pattern of Detected Landslides with Respect to Hydrological Factors
On an annual scale, precipitation is widely recognized as the driver of the seasonal acceleration and deceleration of slow-moving landslides [47][48][49]. However, precipitation may not be the only factor to initiate a landslide that has been constantly active for hundreds of years [50]. A comparison was made of 20-year average precipitation (2001-2020) relative to the landslide distribution ( Figure 12c). Overall, 80% of the identified landslides are located in areas with an annual rainfall of 950 to 1050 mm/a. When rainfall is less than 960 mm/a, the number of landslides increases with rainfall, and when rainfall is greater than 960 mm/a, the number of landslides decreases with rainfall. In addition, the identified landslides were compared with the average excess precipitation between 2017 and 2021 ( Figure 12d). Here, excess precipitation is defined as the difference between annual precipitation and the 20-year average, where a positive value indicates excess. As shown in Figure 13e, it can be found that the landslides being deformed (as identified by InSAR from 2017 to 2021) increased with excess rainfall, indicating that rainfall is one of the driving factors of landslide occurrence. However, regions with less rainfall may experience more landslide activity, partly depending on other factors such as bedrock formation, land uplift, and tectonic activity [51].

Four General Rules of Landslide Distribution
By summarizing the spatial distribution patterns of landslides in this study and other areas (Table 4), four general rules have been determined. The first is the general rule of the distribution of elevation and landslides. There is an elevation boundary that makes the number of landslides reach the peak at that elevation boundary, whereupon the phenomenon of decreasing above that elevation boundary will occur. For example, the elevation boundary of the China-Pakistan Karakoram Highway is 2100 m [52]; along the Jinsha River Basin, the elevation boundary is 2000 m [13]. With increased altitude, due to human activity, river erosion, and rock weathering, the number of landslides is small and the degree of development is slow. The second rule determines that landslides mainly occur in areas with a relative altitude below 500 m, but they are also related to the research area. For instance, in the Bailong River Basin, landslides are most densely distributed in areas with a relative elevation difference of 50 to 150 m [53]. In the Jinsha River Basin, landslides are most densely distributed in areas with a relative elevation difference of 500 to 700 m [13]. The third rule determines that landslides are most frequently developed in the slope range 10 • -45 • because, as the slope angle increases, the rock weathering layer subsequently becomes thinner or non-existent [54] and shear stress increases [55], resulting in a lower number of landslides on steep slopes. The fourth rule determines that landslides tend to occur in mountainous areas with abundant rainfall. For example, on the west coast of the US, 75% of landslides are located in mountainous areas with relatively abundant rainfall (≥2000 mm/a) [56]. However, there are many exceptions. For example, along the China-Pakistan Karakoram Highway area, the rainfall range is 162-941 mm/a, with 84.3% of landslides occurring in areas in the semi-arid region (<400 mm/a) [52]; along the Sichuan-Tibet Railway Ya'an-Linzhi). No significant statistical relationship emerged between rainfall and landslide distribution [39].

Conclusions
In this paper, a technical framework for wide area landslide detection is proposed and successfully applied to an area of approximately 20,000 km 2 along and near three expressways in Sichuan, leading to the following conclusions: