Evaluating expressway stability using interferometric synthetic aperture radar and measuring its impact on the occurrence of geohazards: a case study of Shanxi Province, China

ABSTRACT

This study designs a coherence-based multi-master-image stacking InSAR workflow as an efficient large-scale expressway stability monitoring method, which can improve the overall coherence and number of interferometric pairs via multi-master-image stacking and an average coherence threshold. The average annual ground deformation rate of Shanxi Province, China, was determined using Sentinel-1A data from 2017 to 2021 with a parallel processing strategy to improve efficiency. Compared to global navigation satellite system data, the InSAR results have a root mean square error of 2.90 mm/year. To evaluate the relationship between different factors and the occurrence of geohazards along the expressway network, we conducted a binary logistic regression analysis. The result shows that the deformation rate, annual precipitation, bedrock hardness, distance to stream, and traffic load are significantly related to the occurrence of geohazards along the expressway network in Shanxi Province, among which the deformation rate has the highest significant value, while the traffic load has the lowest, indicating that using InSAR alone to determine geohazards along expressway network is biased and it takes more than InSAR to conduct a comprehensive yet accurate expressway geohazard vulnerability evaluation. This study provides guidance for the road maintenance to ensure the safe operation of expressway systems.

KEYWORDS:

• Expressway stability
• geohazards
• interferometric synthetic aperture radar
• coherence-based multi-master-image stacking

1. Introduction

The expansion of expressways has greatly promoted socio-economic development in China (Faber 2014). However, large volumes and high speed of vehicular traffic destabilize expressways, resulting in severe consequences (Abdelwahab and Abdel-Aty 2002; Qu et al. 2014). Particularly in Shanxi Province, which is located in the eastern part of the Loess Plateau in north China, the wide distribution of collapsible loess has caused serious water and soil loss (Zhang et al. 2015; Li et al. 2020), also contributing to destabilizing expressways. Furthermore, underground mining activities have intensified ground subsidence in goaf areas, thereby jeopardizing the safety of overlying road infrastructure (Xu et al. 2020; Liu et al. 2021; Yi et al. 2021). In addition, vehicular traffic affects the stability of expressways, making them extremely vulnerable to geological disasters (Guo et al. 2019). Therefore, the stability of expressway is in urgent need of routinely monitoring in Shanxi Province.

The surface deformation of an expressway directly reflects its stability (Dai et al. 2018). Currently, using conventional methods (e.g. Global Positioning System and leveling) to evaluate the stability of expressway networks exhibits certain difficulties, including low sampling density, low coverage, and high cost in terms of time and labor (Shi et al. 2014; Liu et al. 2020), which make it unsuitable for large-scale expressway stability evaluation. As an advanced spaceborne remote sensing earth observation approach (Wright, Parsons, and Lu 2004; Lyu et al. 2020; Chen et al. 2021), Interferometric synthetic aperture radar (InSAR) has been applied to ground deformation recognition due to its high accuracy, wide coverage, and capability to operate in all weather conditions (Chaussard et al. 2013; Ge et al. 2014; Wang et al. 2022). The application of multi-temporal InSAR to surface deformation monitoring has achieved particularly good results. Different multi-temporal techniques, such as permanent scattering (Ferretti, Prati, and Rocca 1999), small baseline subsets (Berardino et al. 2002), and stacking InSAR (Sandwell and Price 1998), have not only reduced errors (e.g. atmospheric delay and decorrelation noise) but also enhanced accuracy and further promoted the application of InSAR technology.

In 1998, Usai and Klees (1998) made the first attempt to measure road deformation using InSAR by extracting the phase values of high-coherence pixels along the road and assessing their stability during the monitoring period. Subsequently, InSAR has been extensively used in the deformation monitoring and geohazards vulnerability along road network applied InSAR to the monitoring of highway stability in Tibet and Foshan, respectively, and demonstrated the capability of InSAR in deformation extraction of linear infrastructures. Later, Zhao et al. (2019), Hu et al. (2021) and Hussain et al. (2021) integrated InSAR results with a series of geological factors to study the relationship between geohazard and these factors. However, these studies were typically carried out in relatively small areas (i.e. within a single SAR data frame) and seldom included attribution analysis regarding the geohazard-causing factors.

Due to low efficiency of conventional InSAR approach in large-scale deformation extraction with regard to expressway network, as well as the absence of comprehensive large-scale expressway deformation monitoring research in the Loess Plateau region, in this study, we designed a multi-master-image stacking InSAR workflow specially for large-scale surface deformation extraction and applied it to the surface deformation extraction of the expressway network of Shanxi Province. The InSAR results, along with a variety of disaster-causing factors, were subsequently input into a binary logistic regression model to conduct an attribution analysis of geohazards along the expressway network. In doing so, this study fills the niche of large-scale expressway deformation monitoring and geohazard attribution research and provides a reference for related research in other regions.

2. Study area and SAR data

Shanxi Province is located on the eastern side of the Loess Plateau, North China, and has a total area of 156,700 km² (Figure 1). Shanxi Province features a complex topography, with mountainous terrain accounting for 80% of its total area. Most areas in the province are located at an altitude of over 1,000 m (Li et al. 2013) and are geologically composed of well-developed strata (Guo et al. 2018). Annual precipitation throughout the province is less than 620 mm (Fan and Wang 2011), concentrated in the summer months. The cumulative distance of currently operating expressway in the province is 5,745 km (Figure 1). A total of 236 geological hazards (including ground subsidence, fissures, and landslides) occurred along the expressway network (within 1 km on both sides) from 2017 to 2021, seriously threatening the safety of motorists, vehicles, and infrastructure.

Figure 1. Map of Shanxi Province and data coverage. Sentinel-1 frames are delineated with blue boxes; expressway network is shown in grey line; orange and red dots are the location of subsidence/fissure and landslide, respectively; and the purple triangle denotes the location of GNSS station.

C-band IW (interferometric wide swath)-mode SAR data of Sentinel-1A satellite were used for InSAR processing (Yagüe-Martínez et al. 2016; Radman, Akhoondzadeh, and Hosseiny 2021), and required 10 frames to cover the study area (Figure 1). For each frame, a total of 58 VV-polarized ascending single-look complex images taken from March 2017 to December 2021 were collected at a frequency of one scene per month, the detailed data information is shown in Table 1. The resolution of SAR data was 20 m after multi-looking, with a coefficient of 5 × 1. An external digital elevation model (DEM), namely ALOS PALSAR DEM, was obtained by correcting SRTM using PALSAR data from 2006 to 2011 and up-sampling to 12.5 m resolution was used for topographic phase removal and geocoding. For the purposes of this study, data from 2017 to 2021 were collected from 65 global navigation satellite system stations for the subsequent validation of InSAR results.

Table 1. Detailed data information of SAR acquisitions.

Tile identifier	Acquisition dates*	Total acquisitions
P11F121	20170312 ~ 20211204	58
P11F116	20170312 ~ 20211204	58
P11F111	20170312 ~ 20211204	58
P113F126	20170319 ~ 20211223	58
P113F121	20170319 ~ 20211223	58
P113F116	20170319 ~ 20211223	58
P113F111	20170319 ~ 20211223	58
P40F127	20170314 ~ 20211230	58
P40F122	20170314 ~ 20211230	58
P40F117	20170314 ~ 20211230	58

(Acquisition dates are represented in yyyymmdd format.).

3. Methods

3.1 Coherence-based multi-master-image stacking workflow

The stacking InSAR algorithm was initially proposed by Sandwell and Price (1998), whose workflow is mainly performing DInSAR to time-series SAR images in chronological order, and the surface deformation can be determined by calculating the weighted average of all DInSAR unwrapped phases. And later it has been applied to multiple scenarios including groundwater extraction (Bui et al. 2021), landslide survey (Zhang et al. 2021), permafrost deformation (Dai et al. 2018), due to its simplicity. The algorithm can be demonstrated using equation 1:

(1) $V_{def}=-\frac{\lambda }{4\pi }\ast \frac{{\sum }_1^i\mathrm{\Delta }{t}_i{\phi }_i}{{\sum }_1^i\mathrm{\Delta }{t_i}^2}$(1)

where $V_{def}$ is the deformation rate; the unwrapped phase and the time interval of the $i$-th interferogram are denoted by ${\phi }_i$ and $\mathrm{\Delta }{t}_i$, respectively; and $\lambda$ is the wavelength of the radar echo.

Stacking InSAR is more efficient than other time-series algorithms, and more suitable for deformation extraction over a large area or from a large amount of data. However, the stacking algorithm requires a large number of unwrapped images. Generally, it requires good coherence, which the conventional stacking workflow lacks, to correctly inverse surface deformation (Zebker and Villasenor 1992; Gatelli et al. 1994). Therefore, in this study, we used a coherence-based multi-master-image-stacking workflow, which increased the data sampling rate and improved the overall coherence of interferograms, to extract surface deformation more accurately.

The phase information is more reliable when the overall coherence of the interferogram is higher, leading to a more accurate deformation interpretation (Touzi et al. 1999). The coherence of a pixel $\left(x,y\right)$ is calculated using equation 2:

(2) ${\gamma }_{\left(i,j\right)}=\frac{\left|{\sum }_{x=1}^m\left.{\sum }_{y=1}^{n}M\left(x,y\right)S^{\ast }\left(x,y\right)\right|\right.}{\sqrt{{\sum }_{x=1}^m{\sum }_{y=1}^n|M\left(x,y\right){|}^2{\sum }_{x=1}^m{\sum }_{y=1}^n|S\left(x,y\right){|}^2}}$(2)

where $M$ and $S$ represent the master and slave images, respectively; * indicates the complex conjugate; and $M$ and $n$ denote the local window size. The average coherence of the entire scene is calculated using equation 3:

(3) ${\gamma }_{avg}=\frac{{\sum }_{i=1}^l{\sum }_{j=1}^w{\gamma }_{\left(x,y\right)}}{l\ast w}$(3)

where $l$ and $w$ are the length and width of the interferogram, respectively.

Interferometric pairs with low coherence (below the average coherence threshold) were ruled out, and those with high coherence were used in the subsequent stacking process (Figure 2). First, the Sentinel-1 SAR images, corresponding orbit state data, and the 12.5 m ALOS PRISM DEM of the study area were prepared and uploaded to the Supercomputing Center of Wuhan University for parallel processing to improve efficiency. Each Sentinel-1 image frame underwent pre-processing, including data decompression and co-registration, to obtain single-look complex image stacks. Second, baseline matching was performed using a multi-master-image strategy. One hundred and sixty interferometric pairs (with minor variations between different tracks) were obtained using a spatial baseline of 10 m, and 80 pairs with better coherence were selected for further processing with regard to the average coherence index. Third, interferometry was performed on these image pairs to obtain interferograms, and the minimum cost flow method was used for phase unwrapping with a threshold of 0.5 to obtain the initial unwrapped interferograms. Fourth, phase correction, including atmospheric fitting and spatio-temporal filtering, was performed to obtain refined unwrapped interferograms. Finally, stacking and geocoding were performed to obtain the frame-wise deformation rate map. These maps were downloaded to a local workstation for mosaicking and adjustment using the mosaic module in ArcGIS software with an averaging strategy applied to the overlapping areas to reduce differences between adjacent frames, and the comparison of the deformation rate before and after the mosaicking and adjustment are shown in Figure 3. After that, the global deformation rate map of the entire study area was obtained.

Figure 2. Flow chart of the stacking workflow used in this study. SLC: single-look complex; MCF: minimum cost flow; SCWU: Supercomputing center of Wuhan University.

Figure 3. Comparison of deformation rate before and after mosaicking and adjustment. Left column is the location indicator of the zoom-ins; A.Before and B.Before are the deformation rate before mosaicking and adjustment; A.After and B.After are the deformation rate after mosaicking and adjustment.

3.2. Binary logistic regression

The impact of different factors on the occurrence of geological disasters along the expressway network in Shanxi Province was evaluated using the binary logistic regression model, which is commonly used to predict the likelihood of a specific event occurring in a region based on the multiple regression relationship between a dependent variable (e.g. a geohazard) and multiple independent variables (e.g. causal factors). It can also be used to determine the extent to which different independent variables influence the dependent variable (Peng, Lee, and Ingersoll 2002). In our case study, 236 random non-hazard points (same as the hazard points) were generated within the buffer zone (Figure 4a). As for the factor layer, after thorough review of existing literatures (Devkota et al. 2013; Banerjee, Ghose, and Pradhan 2018; Zhao et al. 2019; Fatih, Faruk, and Ekrem 2021; Ye et al. 2022) and evaluation of accessibility of data, 7 factors were adopted in addition to the InSAR deformation rate results obtained following the steps in Section 3.1 (Figure 4b), elevation (Figure 4c), gradient (Figure 4d), precipitation (Figure 4e), bedrock hardness (Figure 4F), stream distribution (Figure 4g), fault distribution (Figure 4h), and traffic load (Figure 4i) were collected. Detailed data information and preparation strategy as listed in Table 2. The point attribute table containing all the above-mentioned layer information was then exported to Microsoft Excel. Afterward, the Excel data was imported into SPSS software to conduct binary logistic regression analysis. The function of variables (f) influencing the conditional probability (p) of the occurrence of a geohazard was expressed using a logit transformation (Srimaneekarn et al. 2022):

Figure 4. Factors of geohazards assessed in the binary logistic regression. a. Random hazard points. b. Deformation rate. c. Elevation. d. Gradient. e. Annual precipitation. f. Bedrock hardness. g. Stream distribution. h. Fault distribution. i. Traffic load.

Table 2. Detailed data information and preparation strategy.

Factor layer	Data type	Source	Data preparation
Elevation	Raster	ALOS DEM	The elevation is the same as the external DEM used in the InSAR process, with a resolution of 12.5m.
Gradient	Raster	ALOS DEM	The gradient was generated using ALOS DEM with the help of ArcGIS software.
Precipitation	Vector	National Tibetan Plateau Data Center	The precipitation is generated by averaging the precipitation of Shanxi Province from 2000 to 2020.
Bedrock hardness	Raster	National Geological Archives of China	The bedrock hardness map is generated based on the lithology extracted from the 1:200000 geological map released in 2017.
Stream distribution	Vector	ALOS DEM	The stream is extracted from ALOS DEM using ArcGIS software, and used to determine the distance between point and nearest stream.
Fault distribution	Vector	National Geological Archives of China	The fault is extracted from the 1:200000 geological map released in 2017, and used to determine the distance between point and nearest fault.
Traffic load	Vector	Department of Transportation of Shanxi	The traffic load data is calculated by averaging the annual data from 2017 to 2020, and assigned to the expressway network shape file using ArcGIS.

(4) $f\left(p\right)=In\frac{p}{1-p}$(4)

The logistic regression model was as follows:

(5) $In\frac{p}{1-p}={\beta }_0+{\beta }_1{X}_1+{\beta }_2{X}_2+\cdots +{\beta }_n{X}_n\hspace{0.17em}$(5)

The relationship between independent and dependent variables was evaluated using the significance index obtained from the Wald test, which is the ratio of the estimating parameter ${\beta }_n$ to its standard error (Kleinbaum and Klein 2010; Hosmer, Lemeshow, and Sturdivant 2013). The lower the significance value of a factor, the greater its influence on the occurrence of geohazards.

4. Results

The resulting surface deformation rate map of Shanxi Province is depicted in Figure 5. The average annual surface deformation rate in the study area ranged from −33 to 22 mm/year. Decorrelated areas of the map occurred because of the influence of dense vegetation, layover, or the shadowing effect. The deformation of the expressway network in Shanxi Province was determined using a 1 km buffer zone on both sides of the expressway, as shown in Figure 6. Most of the areas along the expressway are relatively stable, yet major deformation is observable in three areas (Figure 6, a – c), and possible reasons are the disturbance of passing vehicles, underground mining activity in adjacent regions, and early deformation of geological hazards or water and soil loss.

Figure 5. Deformation rate map of Shanxi Province. Negative value represents subsiding, positive value represents uplifting.

Figure 6. Deformation rate map of the expressway network in Shanxi Province. Negative value represents subsiding, positive value represents uplifting.

Although the proposed InSAR workflow takes into account various errors (i.e. topographic and atmospheric errors), the results may still suffer from a small amount of residual terrain error, atmospheric delay error, and noise error. Therefore, to evaluate the performance of the proposed workflow, the average annual vertical deformation rate at 65 GNSS stations was calculated as the ground truth, and the InSAR results were projected vertically for comparison (Figure 7). The root mean square error (RMSE) was 2.90 mm/year.

Figure 7. Comparison between InSAR results and GNSS data.

The results of the binary regression analysis are shown in Figure 8. The variance inflation factor of all independent variables is lower than 2. This indicates no multicollinearity (Lavery et al. 2019), meaning that the independent variables in the regression analysis do not affect one another. Deformation rate, annual precipitation, bedrock hardness, distance to stream, and traffic load all had a significant impact on the occurrence of geohazards, based on a standard of significance value <0.05 (Ou et al. 2021), whereas elevation, slope, and distance to fault did not.

Figure 8. Results of the binary logistic regression analysis for the eight factors. a. Deformation rate. b. Elevation. c. Gradient. d. Annual precipitation. e. Bedrock hardness. f. Distance to stream. g. Distance to fault. h. Traffic load. VIF: variance inflation factor; sig: significance.

5. Discussion

Based on the expressway deformation information obtained from InSAR shown in Figure 6, we calculated the proportion of pixels in each deformation rate interval, as shown in Table 3. It can be seen from the table that the annual average deformation rate within ±10 mm/year accounts for 93.47%, which indicates that most of the expressway network in Shanxi is relatively stable without large deformation. However, it is worth noting that more than half of the expressway network has a negative deformation rate, reaching 54.60%, showing a subsiding trend. In addition, the area with the absolute value of deformation rate greater than 20 mm/year accounts for 0.04%. These areas are suffering from serious deformation, as shown in the zoom-ins in Figure 6. Therefore, although their proportion is relatively small, they cannot be ignored. And the expressway maintenance department needs to carry out deformation treatment in time to ensure road stability and the safety of passing vehicles.

Table 3. Statistical table of pixel count and percentage in each deformation rate interval.

Interval (mm/year)	Pixel count	Percentage
<-20	7725	0.03%
[−20, −15]	212331	0.72%
[−15, −10]	825201	2.80%
[−10, −5]	3547684	12.03%
[−5, 0]	11508127	39.03%
[0, 5]	9607595	32.58%
[5, 10]	2898753	9.83%
[10, 15]	708745	2.40%
[15, 20]	167897	0.57%
>20	3184	0.01%

Among these factors, traffic load had the smallest significance value of 0.002, indicating that this factor has the greatest significance on the occurrence of geohazards. Large traffic loads and high-speed vehicular movement greatly disturb the subgrade, resulting in pavement and slope deformation on both sides of the expressway. The distance to the stream was the second-largest contributor to the occurrence of geohazards (significance value = 0.008). Water erodes the slope toe and softens the soil, which reduces the shear strength of the soil and destabilizes the slope. Bedrock hardness was the key internal factor controlling the occurrence of geohazards as well as their scale and type (significance value = 0.017).

The deformation rate obtained via InSAR had a significance value of 0.036, higher than that of the traffic load, distance to stream, and bedrock hardness, which means that the deformation rate is not as significant to the occurrence of geohazard as the other three factors. There are three possible reasons for this. First, InSAR captures ground deformation in the direction of a radar’s line of sight. Although these data can be converted into a vertical direction, this representation still differs from the sliding direction of geohazards in some cases (Dai et al. 2022). Second, some geohazards occur instantaneously without the slow accumulation of subsidence. In such cases, severe subsidence occurs over a relatively short time, far exceeding the maximum deformation gradient that InSAR can detect. Third, densely vegetated areas along the expressway lead to decorrelation, such that ground deformation cannot be interpreted using InSAR. In future studies, both ascending and descending tracks of SAR data with longer wavelengths should be adopted to address these problems. The significance value of precipitation on geohazard occurrence is equivalent to that of the deformation rate. During a precipitation event, a large volume of surface water infiltrates the slope, increasing its mass and sliding power. Rainwater also infiltrates rock fractures, reducing the cohesion and internal friction angles of the rock and thereby reducing slope stability.

The significance values of elevation and gradient on geohazard occurrence are larger than 0.005, indicating that they are not significant to the occurrence of geohazard. This is because of land leveling, slope cutting, and ditch filling activities during expressway construction, all of which considerably decrease the topographic relief along the expressway. In addition, the distribution of faults in the study area is scarce and concentrated, resulting in a higher significance value.

6. Conclusion

This study designed a coherence-based multi-master-image stacking workflow as an efficient large-scale expressway stability monitoring method to improve the overall coherence and the number of interferometric pairs via a multi-master-image strategy and average coherence threshold. Using Sentinel-1 data obtained from 2017 to 2021, the annual average ground deformation rate was determined for Shanxi Province. Comparing the GNSS data with the InSAR results, the RMSE was 2.90 mm/year. To evaluate the impact of different factors, including the deformation rate, on geohazard occurrence along the expressway network, this study included binary logistic regression analysis. Among the factors considered, the deformation rate, annual precipitation, bedrock hardness, distance to stream, and traffic load all had a significant impact on the occurrence of geohazards along the expressway network in Shanxi Province, indicating that an accurate and comprehensive evaluation of expressway network to geohazard requires the participation of multiple datasets other than InSAR results. This study fills the niche of large-scale expressway deformation monitoring and geohazard attribution analysis in Shanxi Province and provides a reference for future studies in other regions.

Acknowledgements

The authors would like to thank the European Space Agency for providing Sentinel-1 data and Qianxun SI for GNSS data. The authors would also like to thank the anonymous reviewers whose comments helped significantly improve this article.

Disclosure statement

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Data availability statement

Sentinel-1 SAR data can be downloaded from https://scihub.copernicus.eu/; Orbit state data can be downloaded from https://s1qc.asf.alaska.edu/aux_poeorb/; And ALOS PRISM DEM data can be downloaded from https://asf.alaska.edu/; The precipitation data is obtained from National Tibetan Plateau Data Center (http://www.ngac.org.cn/); The bedrock hardness and fault map are obtained from the 1:200000 geological map from National geological archives of China (http://www.tpdc.ac.cn/).

Additional information

Funding

This research was funded by the Natural Science Foundation of China projects (NSFC) (No: 41801397), the National Key R&D Program of China (2018YFC0825803).