Surface Deformation of Expansive Soil at Ankang Airport, China, Revealed by InSAR Observations

Abstract

Ankang Airport is constructed on an expansive soil-fill platform in Shaanxi Province, Central China. Since its completion in 2020, it has suffered surface deformation caused by the consolidation and settlement of the fill layer and instability of the expansive soil slope. Exploring the special deformation law of expansive soil regions by remote sensing and analyzing the deformation characteristics of airports in mountainous areas have always been key issues in related disaster research. Based on the intensity and phase observation data of 37 Sentinel-1 synthetic aperture radar images, this study obtained the spatio-temporal distribution of the deformation of Ankang Airport from May 2020 to October 2021. First, phase optimization was performed on the original interferograms. Second, the persistent scatterer synthetic aperture radar interferometry (PS-InSAR) method was applied to extract the surface deformation information of Ankang Airport, and the accuracy was evaluated. Finally, the singular spectrum analysis method was introduced to jointly analyze the deformation information obtained by the InSAR technology in combination with geological and climatic data. The results show that the excavation area of Ankang Airport was basically stable, the filling area had obvious surface and uneven deformation, and the expansive soil fill slope exhibits deformation characteristics strongly related to slope, rainfall, and fill depth. The deformation was mainly caused by consolidation and settlement, supplemented by the expansion and shrinkage deformation of the expansive soil.

1. Introduction

Ankang Airport is one of the representative projects of Chinese regional aviation network construction, located in Shaanxi Province, China. Expansive soil is mainly distributed in the Ankang Basin where Ankang Airport is located. Therefore, considering the actual cost and transportation situation, the foundation of the airport is built with the expansive soil material. The airport construction lasted 50 months; it is a typical airport located in a mountainous region and, during its construction, excavation, and filling, works with a huge amount of earthwork were carried out. The total filling volume reached 30,000,000 m³, with excavation and filling works staggered. These factors make the uneven land subsidence process likely to occur in the Ankang Airport area, thereby causing related disasters. In contrast to airports located in mountainous regions, expansive soil is widely distributed in the Ankang Airport area. Expansive soil is a disaster-prone geological constituent with obvious expansion and contraction properties and developed fissures [1]. In summary, the expansion and contraction of expansive soil and the staggered distribution of many excavation and filling works have caused potential hazards to the safe and stable operation of the airport. The early research about the stability of Ankang Airport reflect that it has suffered a certain degree of surface subsidence [2]. Therefore, it is necessary to monitor the deformation of Ankang Airport and analyze its stability.

Traditionally, the surface deformation of large-scale artificial structures, such as airports, can be monitored on-site using measurement methods such as GNSS and leveling. Although these traditional geodetic methods have high accuracy, they tend to have poor spatial resolution, omit a large amount of surface information, and are labor-intensive. In recent years, with the rapid development of remote sensing technology with a higher spatio-temporal resolution, various remote sensing technologies have been widely used in deformation monitoring research of large-scale artificial structures. Examples include laser scanning technologies [3], distributed fiber optic sensing (DFOS) [4], and unmanned aerial vehicle measurements [5]. As a new type of Earth observation method, synthetic aperture radar interferometry (InSAR) technology can overcome the shortcomings of traditional observation methods. It can obtain high-resolution surface deformation information and realize large-scale and precise long-term monitoring [6]. Recent studies have shown that time-series InSAR technology is a reliable Remote Sensing method for airport deformation monitoring. For example, Dai and Gao combined InSAR technology with groundwater, active faults, and other factors to interpret the uneven subsidence of Beijing’s Da’xing Airport [7,8]; Marshall combined InSAR technology to study the deformation of the Kuala Lumpur Airport related to the subsidence model and tropical peat intensity [9], and Jiang investigated the long-term reclamation subsidence of Hong Kong’s Chek Lap Kok Airport based on the geological data of reclamation [10]. Furthermore, Wu successfully inverted the 20-year subsidence history of the airport using multi-sensor SAR data [11]; Liu and Zhuo used time-series InSAR technology to analyze the subsidence temporal and spatial characteristics of Xiamen Airport in detail [12,13]. Taking Iqaluit Airport as an example, Short proved that InSAR technology has the ability to obtain the periodic deformation related to special soils such as frozen soil in the airport area [14]. Wu used Sentinel-1 images to study the uneven deformation patterns caused by the construction of Yan’an New Airport’s cut-and-fill project [15]. These studies made full use of SAR data, climate, hydrology, and geological data from various sensors [16] and, based on InSAR technology, interpreted the land subsidence patterns of major airports. The reliability of the time-series InSAR technique for monitoring surface deformations related to soil properties has been demonstrated [16,17].

The Sentinel-1 series of SAR satellite data, provided worldwide for free by the European Space Agency (ESA), has greatly promoted the development and popularization of InSAR technology. However, using Sentinel-1 SAR data alone to monitor airport surface deformation still poses challenges, for example, the long observation period, surface scattering characteristics, and other factors affect the coherence, and stable persistent scatterer (PS) points can only be acquired in a limited area [18,19]. In this study, PS-InSAR technology considering the phase optimization of homogenous points is used, combined with geological and climatic data, to study the surface deformation of Ankang Airport due to filling engineering and the distribution of expansive soil. This study has two goals. First, we combine the Fast Statistically Homogeneous Pixel Selection-synthetic aperture radar interferometry (FaSHPS-InSAR) method with the conventional StaMPS data processing framework [20,21] to obtain the spatial distribution and deformation information of homogeneous ground scattering points at Ankang Airport, and compare the deformation results with the leveling data measurement results to verify the potential of this data method to study the surface deformation information. Second, multi-source data, such as construction data, geotechnical characteristics, regional rainfall data, and SAR interferometric data, were integrated to explore the acquired temporal and spatial characteristics of deformation and its possible evolution process, the possible causes of deformation were analyzed, and stability analysis was conducted.

The remainder of this paper is organized as follows. First, the geological characteristics of the expansive soil distribution area where Ankang Airport is located and the surface deformation caused by high fill are described. Second, the InSAR data used and key technical details of the data processing strategy are introduced. In the next section, the spatio-temporal distribution of the deformation information from Ankang Airport obtained by InSAR technology is introduced, and the accuracy checked using leveling data. The newly compiled results are interpreted and analyzed in combination with the construction and rainfall data. Finally, the conclusions are presented in the final section and recent work is summarized.

2. Study Area

Ankang Airport is located in Ankang City, Shaanxi Province, China. This area belongs to the intersection of the southern foothills of Qinling Mountains and the northern foothills of Bashan Mountains. There is a regionally active fault here: the Yuehe fault. The Yuehe fault was affected by the Dabashan torsional tectonic system in the late Yanshan period. At present, it is a fault that is mainly tensile, has long-term activity, and is still active today. Therefore, the study area itself is at risk of earthquakes (Figure 1). The airport is 15 km northwest of the city center and designed to have an annual passenger throughput of 300,000 passengers and a cargo throughput of 750 tons. The total planning area of the airport is about 128,000 m², the total area of the airport terminal is 5500 m², and the airport runway is 2600 m long and 45 m wide. It is a 4C-level civil airport and was built in 2016 and completed in May 2020. It represents China’s gradual improvement of the civil airports in recent years [22].

The original topography of the airport is dominated by hills, high in the north and low in the south, spanning the Luojia River. Therefore, a large amount of filling work has been carried out in the construction of the airport; the spatial distribution of cut and fill is shown in Figure 2. The volume filled in the study area reached 30,000,000 m³ and the filling area reached 52,000 m², accounting for 41% of the total area of the airport. The maximum filling depth was 48.08 m, and the maximum excavation depth was 45.86 m. The excavation and filling works were staggered, and the maximum vertical height difference in the excavation and filling area reached 93.96 m. The airport fill originated from the undisturbed soil in the adjacent excavation area. The soil was silty clay, dominated by expansive soil. Expansive soil is a naturally formed multi-fissured geological body with significant expansion and contraction characteristics, over-consolidation, multi-fracture, and other undesirable properties [23]. During the construction process of the airport’s interior area, expansive soil improvement and soil compaction took place, the area was concreted, and the strength of the structure was significantly improved. However, there was an area of exposed soil on the expansive high-fill slopes on the south, north, and east sides of the airport. The geological body was dominated by medium-strength expansive soil, and obvious surface fissures developed in this area.

A comprehensive analysis of the hydrological conditions and geological environment in the study area was conducted which showed that it will increase the possibility of disasters, including the concentration of high-intensity rainfall and a large quantity of expansive soils. During the heavy rainfall period, expansive soil absorbs water and expands; during the drought period, it loses water and shrinks. The alternate wet and dry periods cause repeated expansion and contraction, resulting in surface soil fissures. Second, carrying out high-fill engineering under complex landform conditions leads to potential geological problems such as the settlement and deformation of the high-fill foundation caused by consolidation of the original foundation and compression of the structure itself. Therefore, it is necessary to continuously monitor the ground deformation of the airport in combination with precipitation data, geological information, and InSAR observations.

3. Dataset and Methodology

3.1. Dataset

A set of 40 ascending Sentinel-1 images acquired between May 2020 and December 2021 were used to characterize and monitor the ground surface deformation at Ankang Airport. We chose 20210308 as the master image for co-registration with other images. The thresholds of the temporal and perpendicular baselines were set to 60 days and 160 m, respectively. In total, 36 interferograms were generated (Figure 3).

An unmanned aerial vehicle digital elevation model (UAV DEM) with a grid spacing of 3 m and vertical accuracy of 1 m was provided by Shaanxi Jiawei Spatial Geographic Information Technology Co., Ltd., Ankang, China. This DEM is also available to simulate and remove the topographic phase contribution within the interferograms process.

Table 1. Specific information of Sentinel-1A images and DEM data.

Data	Parameters	Description
	Type	SLC
	Track number	84
Sentinel-1	Orbit number	Ascending
	Azimuth resolution	1
	Range resolution	1
	Source	UAV
	Acquired Date	2020/6/11
DEM	Resolution Positioning accuracy	3 (m) 1 (m)
	Elevation accuracy	1 (m)

lists the primary parameters of the SAR dataset and DEM used in this study.

The available and independently monitored ground deformation ground-truth points we obtained are two leveling points (L1 and L2) located in the central area of the airport (in Section 4.2). In addition, we used the rainfall data provided by the Ankang Meteorological Station to analyze the features of deformation time series, the observation data of this station are provided by The China Meteorological Data Service Center.

3.2. Methodology

The intensity and phase information of the SAR images were applied to study the spatio-temporal characteristics of the surface deformation at Ankang Airport. The data processing methods used included original image preprocessing, time-series phase optimization, and PS-InSAR deformation calculation.

First, Sentinel-1A image data were preprocessed using GAMMA software, including image co-registration, topographic phase contribution removal, image resampling, and data de-ramping [24]. As Sentinel-1 adopts the progressive scanning (TOPS) observation imaging mode [25], the Doppler centroid rapidly varies along the track, which leads to differences between different bursts. It is prone to obvious phase ramps [26]. Hence, intensity matching was adopted for the burst overlap areas to ensure that all SAR images were accurately co-registered [27].

In the case of varying surface conditions (e.g., planting vegetation or paving runways), the stability of the scattering properties of the target may be compromised, which reduces the number of PSs being detected [18,28]. Therefore, we designed a data-processing strategy to optimize the interference phase by combining intensity information. First, the intensity images were stacked in a time series to select the homogeneous points. The selection of homogeneous points is an algorithm that measures the similarity between the domain and central pixels by means of statistical inference [20]. Assuming that similar objects have the same backscattering characteristics and phase center, the parameter estimation of similar pixels can be used to improve the signal-to-noise ratio of echo signals in complex scenes [29]. In this study, the FaSHPS-InSAR open-source toolbox proposed by Jiang et al. was mainly used for homogeneous point selection [30]. The theoretical basis of FaSHPS was to transform the hypothesis testing problem into a parameter estimation method for confidence interval estimation under the condition that a Gaussian hypothesis was established. Compared with traditional nonparametric hypothesis testing, it has higher computational efficiency and less uncertainty.

The intensity samples of each scatter point $p$ as the central pixel on the N intensity images were $\left\{A_1,A_2,\dots ,A_N\right\}$, and the sample points for the mathematical expectation $\mathsf{\mu }\left(p\right)$ of $p$ were estimated as $\overline{\mathrm{A}}\left(p\right)=\left(A_1\left(p\right)+A_2\left(p\right)+\dots +A_N\left(p\right)\right)/\mathrm{N}$. To obtain the expected interval estimate, we determined the distribution of $\overline{\mathrm{A}}\left(p\right)$ according to the central limit theorem. As the number of samples N increased, $\overline{\mathrm{A}}\left(p\right)$ gradually approached a Gaussian distribution. Assuming that N was sufficiently large to hold the Gaussian hypothesis, the interval estimate of $\overline{\mathrm{A}}\left(p\right)$ was obtained as shown in Equation (1).

$$\mathrm{P}\{\mathsf{\mu }\left(p\right)-{\mathrm{z}}_{1-\frac{\alpha }{2}}·\sqrt{\frac{Var\left(A\left(p\right)\right)}{N}}<\overline{\mathrm{A}}\left(p\right)<\mathsf{\mu }\left(p\right)+{\mathrm{z}}_{1-\frac{\alpha }{2}}·\sqrt{\frac{Var\left(A\left(p\right)\right)}{N}}\}=1-\alpha$$

(1)

where $\mathrm{P}\left\{.\right\}$ indicates the probability of the interval, ${\mathrm{z}}_{1-\frac{\alpha }{2}}$ is the $1-\frac{\alpha }{2}$ quantile of the standard normal distribution, and $Var A\left(p\right)$ is the variance of $p$. It is assumed that when the image is a single look, the intensity image of the homogeneous area obeys the Rayleigh distribution [31]. Its coefficient of variation $\mathrm{CV}$ is quantitative and can be expressed as

$$\mathrm{CV}=\sqrt{\frac{VarA\left(p\right)}{E^2}}=\sqrt{\frac{4}{\pi }-1}$$

(2)

Therefore, when the dispersion of the backscattering coefficients in the N-intensity image stacks was not large, the scattered intensity samples were considered stable and homogeneous. At this time, Equation (1) is rewritten as Equation (3), which includes only $\mathsf{\mu }\left(p\right)$ in the form of an interval.

$$\mathrm{P}\left\{\mathsf{\mu }\left(p\right)-{\mathrm{z}}_{1-\frac{\alpha }{2}}·\mathsf{\mu }\left(p\right)·\sqrt{\frac{\frac{4}{\pi }-1}{N}}<\overline{\mathrm{A}}\left(p\right)<\mathsf{\mu }\left(p\right)+{\mathrm{z}}_{1-\frac{\alpha }{2}}·\mathsf{\mu }\left(p\right)·\sqrt{\frac{\frac{4}{\pi }-1}{N}}\right\}=1-\alpha$$

(3)

When the value of $\mathsf{\mu }\left(p\right)$ was calculated, the interval represented by Equation (3) was also uniquely determined, such that the confidence interval of $p$ can be estimated. When processing the intensity image stack, it was assumed that the reference pixel is $p_K$, and the number of pixels to be measured is $K$. First, take $\overline{\mathrm{A}}\left(p_K\right)$ as the true value of $p_K$, then estimate the sample mean $\overline{\mathrm{A}}\left(p_K\right)$ of K-1 pixels and compare them with the confidence interval in Equation (3) one by one. The points that fall within the interval are homogeneous points. To reduce the estimation bias, α = 50%.

After the homogeneous point sample selection was completed, it was used as a parameter, and a time sequence phase optimizer was introduced to optimize the phase of the original interferograms. The core theory of the time series phase optimizer was: in the SAR stack of N SLC images, the time vector of each point p was usually considered to obey a complex multivariate Gaussian model with a covariance matrix [32], as shown in Equation (4).

$$\mathrm{f}\left(\mathrm{g}\right)=\frac{1}{{\pi }^n\det \left(\sum \right)}\exp \left(-{\mathrm{g}}^H{\sum }^{-1}\mathrm{g}\right)$$

(4)

where superscript $H$ is Hermitian transposition, the temporal vector is $\mathrm{g}={[{\mathrm{g}}_1,{\mathrm{g}}_2,\dots ,{\mathrm{g}}_n]}^T$, and $\sum$ represents the covariance matrix. One of the main objectives of phase optimization is to estimate $\sum$ using observation matrix $G=\left[{\mathrm{g}}_1,{\mathrm{g}}_2,\dots ,{\mathrm{g}}_l\right]$. The maximum likelihood estimator of $\sum$ is:

$$\widehat{\sum }=\frac{G{G}^H}{l}=\hat{\mathsf{\Gamma }}\circ \widehat{\varphi }$$

(5)

where $\circ$ is the Hermitian transposition,$\widehat{\varphi }$ consists of the phase difference between two acquisitions, and $\hat{\mathsf{\Gamma }}$ indicates a real symmetric matrix in which each element is the coherence between different acquisitions, also known as spatial coherence, can be expressed as

$${\hat{\mathsf{\Gamma }}}_{i,k}=\frac{\left|{{\displaystyle \sum }}_{p\in \mathsf{\Omega }}\exp \left(j\left(\angle g_{i,p}-\angle g_{k,p}-{\varphi }_{fft,p}\right)\right)\right|}{l}$$

(6)

where ${\hat{\mathsf{\Gamma }}}_{i,k}$ is the coherence between $i\mathrm{th}$ and $k\mathrm{th}$ acquisition. $\mathsf{\Omega }$ denotes a set that includes $l$ homogeneous spatial pixels that satisfy the stationary assumption. $\angle g_{i,p}$ is the observed value of $i\mathrm{th}$ SLC images at pixel $p$. ${\varphi }_{fft,p}$ is a phase compensation used to remove the spatial phase ramp in set $\mathsf{\Omega }$ [33]. Next, we use Equation (4) to perform a maximum joint likelihood function on the samples in the set Ω to achieve timing phase optimization [34].

Considering the superiority of the PS-InSAR algorithm, we introduce the phase-optimized interferograms into the traditional StaMPS data processing framework. StaMPS method select PS pixels based on amplitude dispersion and estimated phase stability. Subsequently, the three-dimensional phase unwrapping method was applied to unwrapped the wrapped phase. After the atmospheric phase removal and spatio-temporal filtering, we can perform the extraction of deformation phase [35,36,37]. Finally, we obtain the velocity and deformation time series in LOS direction. The overall data processing flow is shown in Figure 4.

To verify the regional stability of the deformation results, the standard deviation of the mean velocity derived from the PS-InSAR technique was calculated using Equation (7) [38,39]:

$${\mathsf{\sigma }}_{\Delta v}^2\approx \left(\frac{\lambda }{4\pi }\right)\frac{{\mathsf{\sigma }}_{\phi }^2}{M_{{\mathsf{\sigma }}_{Bt}^2}}$$

(7)

where $\lambda$ is the wavelength of the radar wave, ${\mathsf{\sigma }}_{\phi }^2$ indicates the phase dispersion, $M$ represents the number of interferograms, and the vertical baseline deviation of the interferograms is ${\mathsf{\sigma }}_{Bt}$.

In order to explore the relationship between deformation time series and rainfall, this paper introduces the singular spectrum analysis (SSA) method to analyze the deformation time series of Ankang Airport.

Singular spectrum analysis is an orthogonal analysis method that can decompose a set of relatively complex time series into a few, independent, and easy-to-interpret parts. The main principle is shown in Equation (8) [40].

$$R=X{X}^T=US{V}^T$$

(8)

where $X$ is the Hankel matrix constructed from the original discrete time series, $R$ is the autocorrelation matrix constructed from $X$, and U and V are orthogonal matrices belonging to R.$S=\left[\mathrm{diag}\left({\lambda }_1,{\lambda }_2,\dots .\right),\dots ,0\right]$, and $0$ is a zero matrix. ${\lambda }_i$ is the singular value of the autocorrelation matrix R, ${\lambda }_1\ge {\lambda }_2\ge \dots \ge 0$, and the larger the corresponding eigenvalue ${\lambda }_i$ of the non-zero singular value, the larger the proportion of the component. The specific process is as follows: first, remove the noise part of the deformation time series; second, calculate the eigenvalue matrix and decompose the trend deformation and periodic deformation [41]; and finally, analyze the periodic deformation in combination with the regional daily rainfall.

4. Results

4.1. Interferograms Phase Optimization

For the i interferograms obtained by differentially obtaining SLC images, there was a common phenomenon of spatio-temporal decoherence, and it was difficult to obtain sufficient effective PS points. To solve the above problems, this study used a sequential phase optimizer based on the multidimensional complex circular Gaussian distribution described in Section 3.1, to perform phase optimization on the original interferograms stack.

First, spatial samples with similar backscattering properties in space were obtained based on the intensity image set and backscattering coefficient. In this study, the parameter statistics FaSHPS homogenous point selection algorithm was used, which had a high-test effect on the samples. The homogeneous point selection samples are shown in Figure 5.

Using the timing phase optimizer described in Section 3.2, timing phase optimization was performed on all interferograms. A comparison between the timing-phase optimization results and the original phase is shown in Figure 6. The original phase depicted in Figure 6a has a certain spatial coherence in the high-coherence region, but the signal-to-noise ratio was low in the medium- and low-coherence regions. Notably, the phase optimization method used in this paper not only had good results in the high coherence area but also restored part of the signal in the noise area, and the coherence and integrity of the spatial phase were significantly improved.

The coherence distribution of the original interferograms is shown in Figure 7a. The spatial coherence in the airport area was poor, and especially in the high fill slope area, the inner vegetation coverage area, and the runway area, there was a poor coherence distribution. Figure 7b depicts the coherence distribution of the interferograms after the phase optimization. The coherence of the airport area was significantly improved, and the spatial noise had been suppressed to a certain extent; in particular, the spatial continuity of coherence was better in the high fill slope, airport runway area, and terminal area. This illustrates the suitability of the adopted phase optimizer for the phase optimization of the Ankang Airport interferometry.

Using the StaMPS software, we set consistent phase mean and amplitude dispersion threshold parameters and selected stable PS points from the two sets of interferograms. The distribution of PS points identified using the original interferograms set was shown in Figure 7c, and the distribution of PS points identified using the phase-optimized interferograms set is shown in Figure 7d. The identification results of stable PS points showed that, compared with the original interferograms, the phase-optimized interferograms shows higher phase and intensity stability in the study area, and more PS points were identified, especially in low shrubs. Compared to the original interferograms, the high-fill slope area covered by the main vegetation had obvious advantages. Using the data processing strategy in this study, under the same parameter settings, a total of 2129 PS points were identified within the airport range, whereas using the original interferograms, only 1372 PS points were identified. More importantly, in the airport runway and expansive soil high-fill slope area with historical deformation, the number of PS points had increased significantly. This is conducive to better monitoring of the safe operation of the airport functional area and a more comprehensive stability analysis of high-fill slopes.

4.2. InSAR-Derived Results and Reliability Analysis

In this section, the set of interferograms processed by the phase optimization step described in Section 4.1 was solved by the StaMPS method, and good deformation monitoring results were obtained at Ankang Airport. The annual average deformation rate map of the airport from May 2020 to October 2021 is shown in Figure 8, where the deformation rate was in the direction of LOS, the positive value in blue indicates movement towards the satellite, and the negative value in red indicates movement away from the satellite. The InSAR results showed that most of the airport area was relatively stable, and the deformation rate was between (−15–19 mm/year). However, four distinct deformation areas were identified in the high-fill slope area of the airport. Moderate (−27.5–−35.9 mm/year) deformation areas and severe (−35.9–−44.4 mm/year) deformation areas were concentrated in the slope area of the airport expansive soil filling (I, II, III, and IV). The maximum deformation rate appeared in the I area, and the annual maximum deformation rate was −44.4 mm/year. This area was the slope formed by the area with the largest fill layer in the airport, and expansive soil was used as the main filling soil. Therefore, it was necessary to analyze the time-series deformation in combination with the special hydraulic laws of expansive soil. In addition, combined with the distribution map of the airport cut-and-fill shown in Figure 2, the magnitude of the deformation rate was proportional to the depth of the fill layer. The aforementioned deformation laws are discussed in detail in Section 5.

The standard deviation of the calculated deformation rate is shown in Figure 9a. It can be seen from Figure 9a that the standard deviation of the slope area is lower than 5.3 mm/year, the standard deviation of the runway area is 3.5 mm/year, and the overall standard deviation is lower than 5.5 mm/year, showing the reliability of the data processing strategies.

In addition, the leveling observation data of two ground reference points located on the airport slope and north side of the runway at Ankang Airport (L1 and L2 in Figure 9a) were collected. As the observation value of the leveling point reflects the vertical deformation and the InSAR deformation result mainly reflects the deformation along the LOS direction, according to the incident angle of the pixel corresponding to the leveling point, the leveling deformation result was uniformly projected in the LOS direction for accuracy verification. The comparison results of the deformation time series of points L1 and L2 observed by leveling observation and InSAR technology are shown in Figure 9b,c. The difference in the deformation results of the L1 point in the slope area is less than 3 mm, and the difference in the deformation results of the L2 point in the runway area is less than 5 mm. Although the two observation methods are not strictly registered in time, the deformation trend is consistent, further illustrating the reliability of the InSAR monitoring results.

4.3. Stability Analysis of High Fill Slope

The high-fill expansive soil slope stability influence area was distributed around the airport flight area, and the filling foundation was mainly expansive soil, which may have caused deformation related to the hydraulic properties of the expansive soil. The stability of an expansive soil slope area was also a key issue in airport safety operations.

According to the excavation and filling construction conditions (Figure 10a), Figure 10b shows the spatial distribution of the deformation rate of the high-fill expansive soil slope at the airport. The results showed that the deformation rate of most slope areas was relatively low, but InSAR technology identified three obvious deformation areas (Figure 10b in I, II, III, and IV). Area I and II are located on the fill slope on the north and south sides of the airport runway. According to the excavation and filling construction conditions (Figure 2), area I was the area with the largest fill layer depth (48.08 m) in the entire airport with a history of slope failure. The maximum deformation rate of area I was −42.19 mm/year, which was the largest deformation in the entire airport area; area III was located on the fill slope on the east side of the airport. The fill depth of this area was 18.9 m, and the maximum deformation rate reached −37.1 mm/year; area IV was located south of the terminal area, and the engineering slope at the airport expressway had a maximum deformation rate of −37.4 mm/year. The filling depth here was only 14.206 m, but according to previous geological data, the filling soil in this area was mainly strong expansive soil. The above deformation areas appeared in the fill slope, indicating that the consolidation and compression of the fill layer may be one of the potential causes of slope instability.

From the site survey data, we found that there was a certain area of exposed expansive soil in the fill slope area of the airport. Considering the expansion, contraction, and crack-prone characteristics of expansive soils that are highly correlated with soil moisture, this study collected regional rainfall data in the study area from May 2020 to July 2021, and, based on the maximum daily rainfall and rainfall time intensity, this was divided into two rainy seasons from May 2020 to July 2020 and April 2021 to August 2021. Four feature points (P1, P2, P3, and P4 in Figure 10a) were extracted from three regions, and the deformation time series and rainy season time distribution obtained by InSAR technology were combined for analysis. A comparison between the deformation time series of the P1–P4 feature points and the distribution during the rainy season is shown in Figure 11.

P1 and P2 were located on the fill slopes on the north and south sides of the runway area in deformation area A, respectively. As shown in Figure 11a,b, the deformation rates reached −42.19 mm/year and −31.59 mm/year, respectively. The overall deformation trend was away from the sensor. However, in the two rainy seasons, the deformation rate slowed down significantly, and an expansion phenomenon appeared, which was consistent with the soil mechanical properties of expansive soil to a certain extent. P3 showed a similar deformation time series to P1 and P2: During the rainy seasons, there was obvious rebound deformation phenomenon along the LOS direction, which was related to the strong expansive soil distributed in this area. The deformation time series of P4 is slightly different, although the deformation phenomenon away from the SAR sensor also occurred in this area, no obvious rebound phenome-non occurred in the rainy season. The maximum cumulative deformation of P4 reached −53.12 mm, and the deformation rate of P4 was −34.89 mm/year.

To describe the spatial distribution pattern of the high-fill slope deformation in area A in detail, a 3D model was constructed using the UAV to obtain the DEM and DSM information of Ankang Airport, as shown in Figure 12a. The annual average deformation rate of the slope in area A obtained using InSAR technology was registered to the 3D model of the airport, as shown in Figure 12b. The fill slope on both sides of the runway was divided into two parts, north and south. In Figure 12b, area N represents the north-side slope, and area S represents the south-side slope. To analyze the correlation between the spatial distribution of the slope and the slope rate, the deformation rate and 3D model registration results of the slope in the north and south of area A obtained by InSAR were enlarged, as shown in Figure 12c,d, respectively.

As shown in Figure 12c, the spatial distribution characteristics of the deformation of the northern slope were as follows: the deformation rate gradually increases from east to west along the slope; the deformation rate gradually increased with the increase in elevation; and the larger deformation points gathered in the top area at the west end of the slope. As shown in Figure 12d, the deformation rate of the southern slope of area A increased gradually with increasing elevation but did not show obvious spatial distribution characteristics in the east–west direction, and the larger deformation points were gathered in the top area of the slope. The spatial distribution of slope deformation at N and S presented characteristics related to the slope and elevation.

4.4. Stability Analysis of Runway

Ankang Airport currently has a runway with a total length of 2700 m and a total width of 45 m. The ground of the runway is relatively weak and prone to consolidation and compression, leading to uneven ground settlement. From the InSAR deformation monitoring results of the runway area shown in Figure 13b, the airport runway is in a stable state, and only the maximum fill area of the airport has a moderate deformation area (−19 mm/year to −27.5 mm/year).

According to the distribution of fill and cut described in Section 2, P5–P8 were extracted along the runway of Ankang Airport, and four feature points were used for time-series deformation analysis. Based on a comparative analysis of the deformation time series of the feature points in the filling area, the deformation time series of the four feature points and the distribution of the regionally intensive rainfall period are shown in Figure 14.

The deformation rate of P5 is −4.27 mm/year, and the maximum cumulative deformation reaches −18.7 mm. The deformation time series is shown in Figure 14a. In both rainy seasons, P5 had been uplifted to different degrees along the LOS direction, the surface was the largest, the lift amount reached 6.4 mm, and the overall deformation trend tended to be stable. P6 was in the deepest filling area of the airport, the filling depth reached 48.08 m, the deformation rate was −15.05 mm/year, and the maximum cumulative deformation reached −14.13 mm. The deformation time series is shown in Figure 14b. In the rainy season of 2020, there was obvious uplift along the LOS direction, whereas in the rainy season of 2021, the deformation rate slowed down significantly, and the overall deformation trend tended toward subsidence along the LOS direction. Consolidation compression correlation; P7 was in the excavation area of the airport, the deformation rate was 0.19 mm/year, and the maximum cumulative deformation reached 2.3 mm. The deformation time series is shown in Figure 14c. Compared with points P5, 6, and 8, the deformation time series of point P7 is relatively stable and has no obvious correlation with rainfall distribution, and the overall deformation trend fluctuates between ±3 mm. The deformation rate of P8 was −4.07 mm/year, and the maximum cumulative deformation reached 15.4 mm. The deformation time series is shown in Figure 14d. It was similar to P5 and P6, and both have a certain degree of uplift along the LOS direction during the rainy season; however, the fluctuation range of P8 was large, and the cumulative deformation was relatively small (0.79 mm).

We created a section line RR1 (marked in Figure 15a) along the airport runway, obtained the elevation change information of the runway along the section line of RR1 (Figure 15b), and extracted the deformation rate in the 25 m area on both sides of RR1, as shown in Figure 15c. There was obvious uneven deformation along the airport runway. The runway deformation in area A was significantly larger than that in other areas; multiple settlement funnels were observed along the runway, and the maximum deformation rate in the vertical direction was −24 mm/year.

Based on the above analysis, the following conclusions were drawn: (1) the expansive soil in the excavation area is more stable than the expansive soil in the filling area; and (2) the overall deformation of the runway area is sufficiently small to meet functional requirements, such as general navigation.

5. Discussion

5.1. Driving Factors of Ankang Airport Deformation

The common deformation types of large-scale high-fill buildings were post-construction settlement, caused by consolidation and compression of filled foundation, and uneven settlement, caused by discontinuous construction [42].

The deformation rate of the PS point in the study area and the depth of the fill volume at its location were extracted to describe the relationship between the distribution of the fill volume in the Ankang expansive soil airport and the surface deformation, as shown in Figure 16.

As shown in Figure 16, with an increase in the thickness of the fill body, the number of settlement observation points and the degree of deformation increased significantly, and the distribution and thickness of the fill soil affected the distribution and size of the surface deformation. The areas with the most obvious deformation were all distributed in the area where the depth of the fill layer was >30 m, and the deformation monitoring points in this area accounted for 62% of all deformation and settlement monitoring points. Monitoring points with a deformation rate greater than −24 mm/year appear only in the area where the depth of the fill layer is >15 m. By contrast, in the area where the depth of the fill layer is less than 15 m, there are fewer settlement points and lower deformation rates. As the compressive deformation of the fill soil at different depths and the pressure on the undisturbed soil increases with the depth of the fill layer [43], the consolidation and compression of the fill layer is one of the reasons for the surface deformation of Ankang Airport.

In addition, it is necessary to consider the ground deformation caused by the dry–wet cycle of expansive soil [44,45]. The area has a subtropical monsoon climate, with rainfall concentrated between late June and late August. The analysis in Section 4.3 and Section 4.4 shows that the deformation trend of the filling area in the entire airport area is relatively similar and has certain regularity; the deformation rate slows down in the rainy season and then springs back, and there are periodic deformations to different degrees, which is consistent with the hydraulic effect of swelling and shrinking of expansive soil. However, the trend deformation caused by post-consolidation compression is the main component of the deformation time series, and errors are caused by the noise phase in the InSAR settlement process, which hinders the exploration of the relationship between the deformation time series and rainfall. Therefore, this paper introduced the singular spectrum analysis (SSA) method to analyze the deformation time series of Ankang Airport. The P1, P3, P5, and P7 deformation time series for SSA decomposition, and the trend deformation and periodic deformation time series are shown in Figure 17.

P1 and P3 are in the slope fill area, P5 is in the runway fill area, and P7 in the runway excavation area, serving as the control group. The first column (a,c,e,g) of Figure 17 describes the deformation trend of the four points. The deformation trends of P1 and P3 located on the filled slope were similar, showing an obvious settlement trend. The trend deformation of P5 in the fill area of the runway is small, and the trend deformation of P7 does not appear obvious and tends to be stable. The periodic deformation was obtained after removing the trend deformation and high-frequency noise (shown in Figure 17b,d,f,h). When the rainfall increased significantly, P1, P3, and P5 showed an obvious uplift phenomenon; P1 and P3 located on the expansive soil fill slope show obvious periodic deformation, which is similar to the monitoring results of similar cohesive soils using InSAR technology [46]. However, the P7 point located in the excavation area remained stable during the observation period, with no periodic deformation clearly related to rainfall. This is because the expansive soil in the airport runway, terminal building, and other areas has been mixed with lime, cement, or other chemical additives to inhibit the expansion and contraction characteristics of the expansive soil. In addition, the foundation of the runway area has been compacted several times and covered with concrete as a water barrier. However, the expansive soil slope around the airport still has some areas of bare leakage area, which are also the focus of this research [47].

5.2. Potional of InSAR in Reaveling Surface Deformation in Expansive Soil Airport

The above analysis confirms that the deformation information obtained by the data processing method adopted in this paper is consistent with the characteristics of the geological environment.

Monitoring the surface deformation of the airport in expansive soil region by InSAR Technology, the low-coherence effect in runway and vegetation area is the first problem to be solved in data processing [19,48]. This paper adopted the timing phase optimization method described in Section 3.2 to improve the coherence of interferograms. As the result shown in Figure 7, Section 4.1, our method has much denser PS points in the whole airport area, especially in the low-coherence area such as slope and runway. The results are similar to the experimental results of using the likelihood function to estimate the phase at Shenzhen Airport [49]. This illustrates the great potential of adopting methods for airport deformation monitoring. The second problem is that airport deformation is often a complex manifestation of multiple deformations [50]. In order to reasonably interpret the monitoring results of InSAR, the SSA method is introduced in this paper to separate the periodic and trend deformations of the deformation time series. As shown in Figure 17, the relationship between the periodic term deformation and rainfall is more pronounced. This method can help us understand the development law of expansive soil in engineering construction.

However, the proposed method still has its limitations. First, in the selection of homogeneous point samples, the window size needs to be set carefully. A window size that is too small will lead to biased estimates due to few samples, while a window size that is too large will cause the interference signal to be non-stationary and underestimate coherence. Second, the determination of eigenvalues in the SSA method still requires subjective judgment, and we will establish an adaptive SSA method for analysis in future research. In addition, this paper does not consider the topographic error and the effect of temperature on the deformation of concrete structures [11,51]. Therefore, the monitoring of surface deformation of expansive soil at Ankang Airport still needs further research.

6. Conclusions

In this study, interferogram phase optimization and PS-InSAR technologies were used to monitor the ground deformation of Ankang Airport from its completion to the initial stage of operation. Furthermore, SSA technology was introduced to decompose the trend deformation of the fill area and the periodic deformation of the expansive soil from the InSAR time series and analyzed with rainfall information. A total of 39 interferograms were used, and the monitoring period was from May 2020 to October 2021. It fully reveals the temporal and spatial distribution characteristics of the surface deformation after the overall completion of Ankang Airport and provides important deformation information for it. The research results have important guiding significance for airport safety navigation, stability analysis, and future construction. The main conclusions drawn from this study are as follows:

In this study, a method for selecting homogenous points under the assumption of a Gaussian complex circle is introduced, and phase optimization of the interferograms is carried out according to the samples of homogenous points. Compared to the original interferograms, this method identifies more persistent scatterers in the airport. At Ankang Airport, the method used in this paper identified 2129 persistent scatterer pixels, while the original Stamps only obtained 1372 persistent scatterer pixels. The results show that this method can significantly improve the coherence of interferograms and verify the reliability of the monitoring results according to the level monitoring data and standard deviation distribution.

From May 2020 to October 2021, three relatively obvious deformation areas were identified at Ankang Airport, all of which were located on the high-fill slopes of the airport. The deformation trend of these filled slopes is positively correlated with the slope and depth of the fill layer. The runway area of the airport is generally stable, but small-scale land subsidence occurs in the fill area of the runway, which is mainly caused by consolidation and compression of the fill layer in this area. In addition, the analysis results of the deformation time series in different areas show that the expansive soil located in the fill area is more sensitive to rainfall, particularly in the high-fill slope area without integral concrete pouring. There is bare soil in this area, which is conducive to expansive soil disasters. The deformation results showed an obvious periodic deformation.

In general, the monitoring results of the time-series InSAR technology, based on phase optimization, in the Ankang Airport area show that the airport area is dominated by the trend deformation caused by the consolidation and compression of the fill. In the exposed slope and high-fill area, the deformation time series is consistent with the rainfall distribution, and there are characteristics of expansive soil disaster development. At present, the Ankang Airport is still in a stage of rapid deformation. Therefore, it is necessary to monitor its long-term deformation.