Identification of Surface Deformation-Sensitive Features under Extreme Rainfall Conditions in Zhengzhou City Based on Multi-Source Remote Sensing Data

Abstract

Extreme precipitation is one of the most prevalent meteorological disasters occurring today. Its occurrence not only causes significant social and economic losses but also indirectly affects surface deformation, creating safety hazards for diverse ground features. Although there are presently high-precision, comprehensive tools such as continuous scattering interferometry to observe surface deformation, it takes a long time to locate potentially vulnerable objects. A monitoring scheme for surface deformation anomalies was devised to address the timeliness issue of identifying sensitive surface features under extreme rainfall conditions. An SAR image of Sentinel-1A is used to derive the surface deformation in three years before and after a rainstorm in the main urban area of Zhengzhou, and the anomaly surface deformation objects after extreme precipitation are screened to determine the surface deformation-sensitive objects. The results indicate that, in the past three years, a 22.14 km² area in Zhengzhou City has experienced a settlement speed greater than 10 mm/yr. Under the influence of the “7–20” rainstorm in the main urban area of Zhengzhou City, among them, the area of highly sensitive agricultural land for deformation is 2,581,215 m², and there are 955 highly sensitive houses for deformation, with an excellent recognition effect. This method is effective in rapidly locating surface deformation-sensitive or potentially damaged features; it can provide a reference for the vulnerability and risk assessment of buildings.

1. Introduction

In recent years, the tension between climate extremes and urban development has become increasingly prominent [1]. Flooding caused by extreme precipitation not only threatens the safety of people’s lives and properties but the large amount of precipitation may also lead to soil saturation and foundation settlement, which poses a safety hazard to buildings [2]. Accurate and timely updating of the hazardous area or susceptibility maps is an important guide to master disaster data and enhance the ability of emergency response to prevent and mitigate disasters.

Remote sensing technology has enabled the implementation of multi-scale and multi-frequency monitoring methods: from small-scale detection at high update frequencies to large-scale monitoring at low update frequencies and providing high update frequencies when necessary. Satellite-based synthetic aperture radar interferometry (InSAR) is an effective satellite-based remote sensing method for monitoring surface deformations over large areas with millimeter- to centimeter-scale vertical accuracy and bi-weekly or monthly temporal resolution [3]. However, the application of D-InSAR in monitoring surface deformation is usually affected by incoherence and atmospheric effects. In this context, multi-temporal InSAR processing methods have been proposed, such as the permanent scatterer interferometry (PS-InSAR) method [4] and small baseline set (SBAS-InSAR) methods [5]. For urban surface deformation, PS-InSAR is one of the most widely used and mature techniques in the field [6]. For geologically stable cities, surface deformation is a slow process of change, and potential areas of deformation hazards are difficult to identify in a short period of time using existing hazard assessment indicators. For example, the deformation of a slowly subsiding site tends to show a linear downward trend with a certain periodicity in a short inter-annual time period, but when the surface deformation values of a site repeatedly show anomalies beyond the original evolutionary characteristics, the site can be considered as a surface deformation-sensitive area. For example, Li et al. [7] identified deformation anomalies by counting the residual terms decomposed from a dam deformation time series. Shi et al. [8] used a deep learning method to examine the deviation of the expected value of surface deformation from the actual value and judged values larger than the standard deviation of all deviations from the statistics as anomalies. Although surface deformation anomalies may not vary significantly in numerical terms, increasing the number of anomalies can lead to structural fatigue and pose safety hazards for buildings. Therefore, for a city like Zhengzhou, which has been geologically structurally safe and persistently stabilized by human impacts in recent years, surface deformation anomalies can be used to identify hazardous or susceptible areas [9].

In statistics, an outlier is a data point that does not belong to a particular group; it is an anomalous observation that is very different from the other values and departs from a well-formed data set. There are various methods for outlier testing, and the mainstream ones are as follows: 3sigma, Z-score, boxplot, Grubbs hypothesis testing, KNN, LoF, COFsos, DBSCAN, iForest, PCAA, utoEncoder and so on. Each of these methods has its own advantages and disadvantages, among which the boxplot method has been widely used in the study of anomaly testing by virtue of its computational simplicity, high efficiency and the reduced influence of anomalies (excellent detection of anomalies on sea surface temperature, spatio-temporal precipitation and general circulation model data) [10]. The robust skewed boxplot method can be improved to be more accurate in detecting anomalies in skewed distributions [11]. However, if statistical methods are to be used to identify outliers in time series, it is necessary to ensure that the time series is stationary. However, the evolution process of surface deformation values in cities often exhibits periodicity and trends. Therefore, the time series obtained by removing trend time-series elements and periodic time-series elements from the original surface deformation time-series data is not in line with the evolution characteristics of the original surface deformation time-series data. We call it residual, and the number of values that exceed the threshold specified by the boxplot is the outlier.

In order to obtain the residuals, it is necessary to decompose the time series. A number of time-series analysis methods for signal decomposition have been proposed for successful application in geoscience research [12]. The main methods that can be used for time-series analysis include least squares (LS), Kalman filtering (KLF), wavelet decomposition (WD), empirical modal decomposition (EMD), and seasonal trend decomposition process (STL) with locally weighted regression [13]. In addition to this, singular spectrum analysis (SSA) is also widely used in geoscience time-series analysis [14], including trend detection at different temporal resolutions, extraction of seasonal term components, simultaneous extraction of residuals from small and large cycles, identification of periodicities with different magnitudes, extraction of complex trends and periodicities, analysis of short time series with structural patterns, and detection of change points. Unlike other decomposition methods, SSA does not rely on assumptions about the time series. It can, therefore, be used to disentangle anomalous components in time series. For example, Sahoo et al. used SSA to discover that 220Rn time-series data exhibit anomaly behavior before earthquakes occur, which can provide assistance in understanding the seismic events that occurred during the study period [15]. In addition, another study demonstrated the potential of using SSA and Fisher–Shannon statistical methods to effectively capture and estimate internal vegetation anomalies in forest cover [16]. Similarly, in meteorology, the decomposition of the residuals of non-hydrostatic delay time series using SSA can be used to predict and understand dust storm occurrence [17]. Thus, SSA is an auxiliary method to effectively identify time-series anomalies.

In order to quantitatively monitor the surface deformation susceptible areas along the extreme precipitation conditions, this paper firstly utilizes the InSAR technique to obtain the spatial and temporal transformations of surface deformation in Zhengzhou city for a total of three years before and after 20 July 2021; after that, the singular spectrum analysis (SSA) method is used to decompose the time series of the grid points in each area, and the residual term is used to recognize the residuals under extreme precipitation conditions by using the statistical test method. The residuals were used to identify the residuals under extreme precipitation conditions, and the locations with multiple occurrences of anomalous information were recognized as highly sensitive areas of surface deformation. This study can help us quickly, and with low cost, identify areas of abnormal surface deformation, especially for some buildings, and provide a reference for relevant departments to identify areas of surface deformation with safety hazards.

2. Study Area and Data

2.1. Study Area

Zhengzhou City is located in the central north of Henan Province (Figure 1a), at 113°27′~113°51′ E, 34°36′~35°00′ N between the urban area north of the Yellow River, west of the Mount Song, the topographic trend from southwest to northeast gradually low, a step decline, an area of about 1010 km², and the strata are mainly composed of Quaternary loose silt, silty clay, and gravel layers. Adverse rock and soil conditions such as collapsible loess are mainly distributed in the western part of Zhengzhou, as shown in Figure 1b; the climate conditions belong to a temperate continental monsoon climate, with an average annual precipitation of approximately 632 mm, and the precipitation time in the year is mainly concentrated in June–August [18]. The Jalu River, Dongfeng Canal, Jinshui River, Xionger River and Soxu River in the city undertake the task of precipitation and wastewater discharge. Further, 77% of the city’s features is covered by buildings and 13% is cultivated land.

From 18:00 on 18 July to 0:00 on 21 July 2021, there was a rare continuous heavy rainfall weather process in Zhengzhou. The maximum hourly precipitation observed by Zhengzhou Station from 08:00 on 20 July to 08:00 on 21 July reached 201.9 mm/h (16:00–17:00), breaking through the historical extreme value on the Chinese Mainland. During the 720 extremely heavy rainstorm disaster in Zhengzhou (hereinafter referred to as the disaster), the city’s cumulative average precipitation was 449 mm, but a single day of precipitation broke through after 1951. In Zhengzhou, since the establishment of the station, with 60 years of historical records, the city’s drainage system has been seriously overloaded, the city flooding has been serious, and 380 people died due to the disaster, accounting for 95.5% of the province; the direct economic loss was CNY 40.9 billion, accounting for 34.1% of the province, and multiple urban roads and buildings were damaged to varying degrees. In this paper, the main urban areas of Zhengzhou City (Jinshui District, Huiji District, Zhongyuan District, Erqi District, Guancheng District) were selected for monitoring.

2.2. Data Sets

In this study, the PS-InSAR analysis utilizes Sentinel-1 SAR data provided by the European Space Agency (ESA). The Sentinel-1 SAR operates at C-band frequencies, and Sentinel-1A consists of four modes of operation: SM, IW, EW, and WV, including interferometric wide-field (IWF) modes with a resolution of 5 m × 20 m and an amplitude of 250 km. mode with a revisit period of 12 days. A total of 90 single-view complex (SLC) images acquired from ascending orbits in interferometric wide-field (IW) mode were used in this study. The study area has an incidence angle of approximately 38.96°, a distance-directional resolution of 2.32 m, and an azimuthal resolution of 13.91 m. These images cover a time span from 14 July 2019 to 12 June 2022, with a temporal resolution of 12 days, and the data were obtained from the Alaska Aeronautical Communications Agency. The DEM is the U.S. Shuttle Endeavour’s radar topography mapping SRTM (Shuttle Radar Topography Mission (SRTM)) data, providing elevation information at a spatial resolution of 30 m. Precision orbit files were obtained from the ESA Copernicus Program data distribution site.

In order to analyze the causes of urban deformation anomalies, we reviewed the “Environmental Geological Survey and Evaluation of Major Cities in Henan Province-Zhengzhou City Summary Report” and drew a sketch map of the distribution of undesirable geotechnical bodies in Zhengzhou City (Figure 1b). The rainfall data were obtained from the NOAA Global Environment website, the subway construction information was obtained from the Zhengzhou Municipal Bureau of Statistics, and the land cover type data were obtained from the global land use cover data produced by ESRI based on the Sentinel-2 data using the deep learning method. The 2019 housing vector boundaries of Zhengzhou City were obtained from web sharing and checked for accuracy to be used for the research in this paper. Optical images obtained using Google Earth and street photos of June 2021 provided by Baidu Maps were used to study the changes in urban construction.

3. Methods

This article obtained the cumulative surface deformation data of Zhengzhou City using VV polarized Sentinel-1A images over a time span from 14 July 2019 to 12 June 2022, using the PS InSAR method. Then, we interpolated into grid form and applied SSA method to extract residuals on each grid time series. Afterwards, the boxplot method was used to calculate all residuals on each grid and filter out outliers.

3.1. InSAR Processing

In this study, the PS-InSAR method in the ENVI SARScape software (Version: 5.6) was used to retrieve temporal deformation of the study area based on Sentinel-1A images. The PS-InSAR method was first proposed by Ferretti et al. [4]. Although various algorithms for PS-InSAR techniques have been developed in recent years, all of them aim at retrieving localized deformations on highly coherent scatterers (i.e., persistent scattering (PS) pixels, e.g., buildings or bare rocky objects) from the contained differential interference phase.

In the PS-InSAR interferometric model, the phase at any point can be expressed as:

ϕ = ϕ_defo + ϕ_atmo + ϕ_topo + ϕ_geo + ϕ_noise(1)

where ϕ represents the interferometric phase, ϕ_defo represents the deformation phase of the surface line of sight, ϕ_atmo represents the atmospheric disturbance phase, ϕ_topo represents the terrain phase, ϕ_geo represents the flat phase, and ϕ_noise represents the thermal noise phase of the system.

By using the PS-InSAR technique in the ENVI SARScape software, the SAR image acquired on 4 January 2021 was selected as the master image, and the other images were co-aligned with the master image. All images were aligned to the master image coordinate system using a DEM to remove the terrain leveling effect. Then, we used a linear model to estimate the deformation rate and residual elevation information in all differential interferograms and calculated the pixels with high coherence. Then, the estimated target elevation was estimated, and the deformation rate in the LOS direction and the atmospheric compensation was generated again with the deformation results. Finally, all the steady point scattering points (PS points) were geocoded and transformed to the corresponding positions. PS-InSAR technique is quite mature in urban surface monitoring, and its data confidence meets the experimental requirements.

3.2. PS Point to Raster Conversion

The generated 1.25 million PS points are interpolated into 20 m image elements, and since the number and density of PS points are sufficiently high, inverse distance weight (IDW) interpolation is used; the basic principle is:

Firstly, calculate the distance from the unknown point to all known points.

(2)$d_i=\sqrt[2]{(x-x_i{)}^2+(y-y_i{)}^2}$

where d_i is the distance from the unknown point to the i-th known point, x_i and y_i are the coordinates of the i-th known point, and x and y are the coordinates of the unknown point.

Then, calculate the weight of unknown points based on distance.

(3)$w_i=\frac{1/d_i}{\sum _1^n 1/d_i}$

The weight from the unknown point to the i-th known point is then calculated.

Finally, calculate the value of the interpolation point.

(4)$Z_0=\sum _{i=1}^{i=n} w_i\ast Z(x_i,y_i)$

where Z (x_i, y_i) is the value of the i-th known point.

Each interpolated point can be obtained using the above method.

Then, the 90-period raster layers are converted into netCDF files, which improves the computational efficiency and enables clearer results. In addition, the surface deformation hazard grading index is formulated with reference to the “Specifications for risk assessment of geological hazard”, as shown in Table 1 below. This specification was proposed by the Ministry of Natural Resources of the People’s Republic of China and is applicable to the risk assessment of geological disasters in China. It stipulates that when the local surface deformation rate is less than 10 mm/yr, it is a low-risk situation of weak deformation development; when the surface deformation rate is 10–30 mm/yr, there is a moderate risk of moderate deformation and weak development; when the surface deformation rate is greater than 30 mm/yr, it is a high-risk situation for strong deformation development. The following analysis focuses on the areas where the surface deformation rate is greater than 10 mm/yr.

3.3. Time-Series Decomposition

Singular spectrum analysis (SSA) is a signal processing technique based on singular value decomposition for analyzing time-series data, which can decompose time-series data into trend, period, and residual terms. Firstly, for the one-dimensional time series, x_t (1 ≤ t ≤ N) is embedded into the delayed sequence through a sliding window to form a trajectory matrix, where L is the window length and K = N − L + 1.

(5)$\mathrm{X}=({\mathrm{x}}_{\mathrm{ij}}{)}_{\mathrm{i}, \mathrm{j}=1}^{\mathrm{L},\mathrm{K}}=\left(\begin{array}{ccccc}{\mathrm{x}}_1& {\mathrm{x}}_2& {\mathrm{x}}_3& \cdots & {\mathrm{x}}_{\mathrm{K}}\\ {\mathrm{x}}_2& {\mathrm{x}}_3& {\mathrm{x}}_4& \cdots & {\mathrm{x}}_{\mathrm{K}+1}\\ ⋮& ⋮& ⋮& \cdots & ⋮\\ {\mathrm{x}}_{\mathrm{L}}& {\mathrm{x}}_{\mathrm{L}+1}& {\mathrm{x}}_{\mathrm{L}+2}& \cdots & {\mathrm{x}}_{\mathrm{N}}\end{array}\right)$

The singular value decomposition (SVD) is then performed on the trajectory matrix to decompose the trajectory matrix into the following form:

(6)$\mathrm{X}=\mathrm{U}\mathsf{\Sigma }{\mathrm{V}}^{\mathrm{T}}$

where U ∈ R^L×L, Σ ∈ DIAG^L×K and V ∈ R^K×K, U and V are called You matrices and are unitary orthogonal matrices, and Σ is called singular value and is a diagonal matrix.

The singular value decomposition decomposes the trajectory matrix into a linear combination for the unitary matrix, the diagonal array, and the unitary matrix; then, X can be transformed into:

(7)$\mathrm{X}=\sum _{\mathrm{i}=1}^{\mathrm{r}} {\mathsf{\sigma }}_{\mathrm{i}}{\mathrm{U}}_{\mathrm{i}}{\mathrm{V}}_{\mathrm{i}}^{\mathrm{T}}=\sum _{\mathrm{i}=1}^{\mathrm{r}} {\mathrm{X}}_{\mathrm{i}}$

where r is the rank of the matrix, σ is the square root of the eigenvalues of the covariance matrix of the trajectory matrix from smallest to largest, and d = max(i), also known as the singular spectrum of the original sequence.

Divide the set of subscripts into m disjoint subsets and sum the matrices contained within each group, and the other i₁…i_p are the matrices contained in group I. Then

(8)${\mathrm{X}}_{\mathrm{I}}={\mathrm{X}}_{{\mathrm{i}}_1}+\cdots +{\mathrm{X}}_{{\mathrm{i}}_{\mathrm{p}}}$

(9)$\mathrm{X}={\mathrm{X}}_{{\mathrm{I}}_1}+\cdots +{\mathrm{X}}_{{\mathrm{I}}_{\mathrm{m}}}$

Finally in the reconstruction step, we transform each matrix X_Ij, in Equation (5) into a new sequence of length N, i.e., we obtain the decomposed sequence. Define Y as an L*K matrix with elements y_ij, 1 ≤ I ≤ L, 1 ≤ j ≤ K. Let L* = min (L, K), K* = max (L, K), and N = L + K − 1. If L < K, y^*_ij =y_ij; otherwise, y^*_ij = y_ji. Then, the formula for transforming the matrix into a one-dimensional sequence of diagonal averages can be expressed as:

(10)${\mathrm{y}}_{\mathrm{k}}=\left\{\begin{array}{l}\frac{1}{\mathrm{k}}{\displaystyle \sum _{\mathrm{m}=1}^{\mathrm{k}}} {\mathrm{y}}_{\mathrm{m},\mathrm{k}-\mathrm{m}+1}^{*}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}1\le \mathrm{k}<{\mathrm{L}}^{*}\\ \frac{1}{{\mathrm{L}}^{*}}{\displaystyle \sum _{\mathrm{m}=1}^{{\mathrm{L}}^{*}}} {\mathrm{y}}_{\mathrm{m},\mathrm{k}-\mathrm{m}+1}^{*}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}{\mathrm{L}}^{*}\le \mathrm{k}\le {\mathrm{K}}^{*}\\ \frac{1}{\mathrm{N}-\mathrm{k}+1}{\displaystyle \sum _{\mathrm{m}=\mathrm{k}-{\mathrm{K}}^{*}+1}^{\mathrm{N}-{\mathrm{K}}^{*}+1}} {\mathrm{y}}_{\mathrm{m},\mathrm{k}-\mathrm{m}+1}^{*}\hspace{1em}{\mathrm{K}}^{*}<\mathrm{k}\le \mathrm{N}\end{array}\right.$

The final one-dimensional sequence is obtained as:

$({\mathrm{y}}_1,{\mathrm{y}}_2,{\mathrm{y}}_3,...,{\mathrm{y}}_{\mathrm{k}-1},{\mathrm{y}}_{\mathrm{k}},{\mathrm{y}}_{\mathrm{k}+1}...,{\mathrm{y}}_{\mathrm{L}-1},{\mathrm{y}}_{\mathrm{L}},{\mathrm{y}}_{\mathrm{L}+1},...,{\mathrm{y}}_{\mathrm{K}+\mathrm{L}-3},{\mathrm{y}}_{\mathrm{K}+\mathrm{L}-2},{\mathrm{y}}_{\mathrm{K}+\mathrm{L}-1})$

The remaining principal component sequences can be reconstructed similarly. In this way, the signal decomposition of the time series is completed.

3.4. Box-and-Line Diagram Guidelines

Box-and-line plots are a simple method for identifying outliers using outlier cutoff point judgments. This graph-based approach to outlier identification is appealing not only because of its simplicity but also, and more importantly, because it does not use extreme potential outliers that could distort the calculation of outlier metrics and reduce sensitivity to outliers, so data outside the inner bounds are usually defined as outlier data.

According to previous studies, Zhengzhou City surface deformation has a certain trend and periodicity [19,20]; the residual is not in line with the development of surface deformation law of anomalies. We allow for the existence of anomalous values in the process of surface deformation development, but the fluctuation of such anomalous values can only be permitted within a certain range; however, according to the second law of geography, the rule of change of surface deformation in each region exists in a heterogeneous manner, so that every location has its own independent range of permissible residual fluctuations. For this study, the surface deformation in the study area is interpolated into a raster (Figure 2a), and each pixel in the raster represents the change in surface deformation values at that location (Figure 2b). The fluctuation range of outliers for each pixel is the inner limit corresponding to the residual statistics in the time series of the corresponding pixel, and values exceeding the allowable fluctuation range (inner limit) are labeled as outliers (Figure 2c).

4. Results

4.1. Spatial and Temporal Changes on the Surface of Zhengzhou City

During the period from July 2019 to June 2022, the main areas of surface uplift in Zhengzhou City were distributed in the central and western parts of the city, while the subsidence areas were mainly concentrated in the northern and eastern parts (Figure 3a). The area of weakly developed areas with deformation rates within 10 mm/yr was 1504 km², accounting for 98% of the total area and belonging to the low-risk area. The area of moderately developed areas within the range of 10–30 mm/yr. is 22.47 km². These areas with moderate risk of deformation are mainly distributed in the eastern and northeastern parts of Zhengzhou City (Figure 3b). They are distributed in a block-like manner and exhibit subsidence phenomena. In addition, the area of hazardous areas with a settlement rate greater than 30 mm/yr. only accounts for 0.1% of the total area (as shown in Figure 3c,d). According to the risk assessment indicators for surface deformation, the surface situation in Zhengzhou City is characterized by low risk in most areas, but there are still some areas with moderate risk of settlement, and there are almost no high-risk surface deformation areas. Therefore, the following will only discuss the situation of moderate-risk settlement areas in detail. In order to visually analyze the development of deformation in the area of moderately hazardous subsidence, the PS points in the area with larger patches were sampled, and three sampling points (P01, P02, and P03) were selected, the distribution of which is shown in Figure 3b, and as a comparison, P04 was selected in the city center with weakly developed ground uplift.

P01 is located in the east of Guxing Town (Figure 4a), surrounded by mostly agricultural land and low-rise buildings, which are also intersections of multiple expressways. The cumulative deformation curves and average cumulative deformation curves of all sampling points in this area gradually decrease over time (as shown in Figure 4b), and the annual settlement rate is concentrated at 12 mm/yr. P02 is located to the east of Jialu River Station (as shown in Figure 4c), surrounded by newly developed high-rise residential buildings with low population density. Similarly, the sampling point curve in this area shows a downward trend (Figure 4d), with an annual settlement rate concentrated at 11 mm/yr. P03 is located in Yangqiao Village (Figure 4e), surrounded by newly developed high-rise residential buildings with a concentrated population density. There is a high demand for agricultural water in this area, and the main source of agricultural water comes from groundwater extraction. The sampling point curve in this area shows a downward trend (Figure 4f), and the annual settlement rate is concentrated at 11 mm/yr. P04 is located at Henan University of Education (Figure 4g). Although it is surrounded by high-rise office buildings and has a high population density, the sampling point curve in this area shows an upward trend (Figure 4h), with an annual uplift rate concentrated at 6 mm/yr. It is a typical area of surface uplift after the recovery of groundwater levels in recent years [21]. Therefore, the moderately dangerous subsidence areas in Zhengzhou City are mainly located near agricultural land with high demand for groundwater and the development of new residential areas.

The precipitation of Zhengzhou City from May 2019 to June 2022 is shown in Figure 3e, where the gray background is the rainy season of that year, and, similarly, in Figure 4b,d,f,h, the gray background is the rainy season of the corresponding year. By observing the surface deformation of each sampling area, even around 20 July 2021, it is difficult to detect anomaly surface deformation fluctuations based on existing indicators.

4.2. Identification of Sensitive Features

The change in the number of surface deformation anomaly grids in Zhengzhou City is shown in Figure 5a, and its fluctuation is generally the same as the change in precipitation; the number of deformation anomaly grids in the rainy season is more than before, and the larger the precipitation, the more anomalous grids. However, the number of anomalous grids appears to be abnormally high at the beginning and the end of the curve, especially in the first period, which is the starting period where the values inside the grids are all zero, but there are still about 4.5 × 10⁵ grids determined as anomalies, which is illogical. This error occurs because when SSA decomposition is performed, ideally, the sum of the trend and period terms at the beginning of the curve would be equal to the actual values, and the trend and period terms would be unbiased over the entire time series, but these two terms at the beginning and the end of the curve cannot be modeled completely and accurately, and the un-modeled factors are incorporated into the residual terms (Figure 5b), which are not the residuals from the real surface deformation, and, thus, these large numbers of false residuals are categorized as outliers in the boxplot statistics, leading to this error. However, as the sequence advances, this false residual rapidly becomes smaller until the problem described above reappears at the tail end of the approach. Although the false residuals are not conducive to counting the outliers in the first and last parts of the sequence, this is not entirely a bad thing because the boxplot counts all the residuals of the entire time series, and the false residuals can be introduced as noise, which makes the judgment indexes in the counting of outliers stricter and improves the accuracy of the identification of outliers in the middle of the sequence.

Therefore, the number of anomalous grids for the two rainy seasons in the middle of the sequence is judged with high accuracy, and the following analysis focuses on the change in the number of anomalous grids for the rainy season period in 2021.

Since the surface cover of Zhengzhou City is mainly covered by buildings (76.5%), crops (11.96%), and grasslands (6.5%), and the remaining surface cover types are only 5% of the total area (of which the water area accounts for 4.48%), only the anomalies of the anomalous building area, the crop area, and the grassland area will be analyzed and, in addition, due to the availability of the house vector boundaries (definition: if the position of a house intersects with the position of an anomaly pixel, then the house is an anomaly house). In addition, due to the deviation between the vector boundary of the house and the actual situation, the anomalous area of building area will be replaced by the number of anomalous buildings. In order to judge the surface deformation anomalies in a more refined way, it is defined that two or fewer deformation anomalies occurring at the same time at the same coordinates of the pixels in the same grid in the rainy season of 2021 are recognized as low-sensitive pixels of surface deformation, and more than two are recognized as high-sensitive pixels; at the same time, the anomalies occurring are classified into subsidence anomalies and upliftment anomalies (e.g., Figure 5d). The area and number of anomalies of various features are shown in Figure 5c and Table 2. Among them, “High sensitive areas” represent the number of times anomalies occur in pixels where different ground objects are located, and 05, 04, and 03 represent the number of times abnormal subsidence occurs, which is 5, 4, and 3. Thus, 5, 4, and 3 indicate that the number of times abnormal uplift occurs is 5, 4, and 3. “Total area” represents the total area of pixels with different abnormal occurrences. “Crops” and “Rangeland” represent the areas where different abnormal deformations occur in crop and grassland areas, respectively. “Number of buildings” represents the number of houses with different abnormal deformations.

4.2.1. Crop and Grassland Areas

In Zhengzhou City, for 99.7% of the total area for the low-sensitive area, the occurrence of low-sensitive features to buildings and most of the anomalies are mostly subsidence and crops, and the grassland area accounted for only 7%. However, in the highly sensitive area, there was 72% of the area for the crop area, which occurred three times in the anomaly settlement of the crop area, which accounted for 65% of the total sensitive area. The reason for this result is that a strata house surrounded by farmland in the northern part of Zhengzhou City experienced three obvious anomaly deformations during the rainy season of 2021 (Figure 6b), but since the surrounding features are farmland with poor backscattering, it is not possible to compute its PS points, and, thus, a large number of blank zones of PS points are generated around the building, but since the inverse distance weight interpolation method assigns a distance based on its positional weights, it spreads these three anomalies to a large number of blank zones around the building (Figure 6a), resulting in an error. Similarly, similar agricultural land patches are distributed in a large number of sparsely built-up areas in Zhengzhou City, but the detection of their anomalies may be subject to error, and, therefore, is only used as a reference. After measurement, the area of this blue area is about 1,520,000 m². After correction, the highly sensitive area is dominated by anomalous subsidence phenomena, and the area of highly sensitive agricultural land and grassland is slightly higher than that of the built-up area, which is the main contributor to anomalous subsidence features.

4.2.2. Building Areas

During the rainy season of 2021, the above methods identified a total of 875 houses that showed highly sensitive surface deformation, with almost equal numbers of houses with anomaly settlement and anomaly uplift. These houses are mainly located in the old and densely populated urban areas of the Guancheng and Erqi districts, and most of them are abandoned and old houses (Figure 7).

In order to observe the identification of highly sensitive buildings more intuitively, another four highly sensitive buildings are randomly selected to observe the cumulative surface deformation changes at their respective nearest PS points, and the results are shown in Figure 8. House a is located at the intersection of Longhai Expressway and Yingxiang Road (Figure 8a), the building is a colored steel tile house (Figure 8e), and it is nearby bare ground under construction. This kind of house has no foundation and is easily subject to the settlement phenomenon by rainwater in the rainy season, and three anomaly settlements occurred during the rainy season in 2021 (Figure 8i). House b is the Zhengzhou Book Purchase Center (Figure 8b), a commercial building built in 1999, located in the busy downtown area, with many urban underground space projects and building construction projects nearby, and its own underground structural condition is easily affected (Figure 8f); similarly, three anomaly settlements occurred during the rainy season in 2021 (Figure 8i), and three anomaly settlements occurred during the rainy season in 2021 (Figure 8i). House c is a dwarf layer building next to Building 1 of the Family House of the Chinese Medicine College (Figure 8c), which is surrounded by dense old buildings, and it can be seen from Figure 8g that the structure of this dilapidated brick house is not very stable anymore; with the poor infrastructure and drainage capacity of these old urban areas, the house experienced three anomaly uplift phenomena in the rainy season. House d is similar to house c (Figure 8d), which is also a dilapidated building in the old urban areas, and it already has obvious cracks on its walls (Figure 8h), indicating that its foundation structure has been unstable for a long time. According to the actual situation of these sampled buildings, we verified the identification ability of potential damaged buildings, and the scheme proposed in this paper works well in terms of identifying surface deformation-sensitive buildings.

5. Discussion and Limitations

Although the coupled results of SSA decomposition are consistent with the actual data, the first and last parts of the individual decomposition results do not fit perfectly with the real situation, although the effect of the error on the time period of surface deformation in this paper is relatively small; when the time series is shorter, this error will inevitably be spread to the corresponding time period of the study, and the method of optimizing away this imperfect fit is a fundamental means of solving the problem of the large error. In addition to this, extending the time series as much as possible is a direct means of minimizing the error, as well as more accurately identifying the data for the trend term, thus further improving the accuracy of the residual term.

Due to the limitation of the features in the study area, the PS-InSAR method is not good for the deformation detection in the non-building area. Of course, it can be combined with the SBAS-InSAR method, and the SBAS method can be utilized to fill in the area with poor spatial continuity of the PS points, but the poor correlation of the SAR data at this location will still lead to the lack of data or distortion in these areas. It is necessary to further obtain the surface deformation information from other observation means to improve the monitoring accuracy.

Most areas of Zhengzhou City are covered by soft soil, loess, and general soil, which mainly contain clay minerals, such as illite, kaolinite, and montmorillonite. For some houses with unstable foundation structures, such as ground cracks, the presence of these minerals may lead to abnormal uplift of the house due to the infiltration of precipitation into the underground of the house, which may cause the soil to absorb water and expand. For some houses with old or damaged structures, precipitation may often exacerbate structural damage, leading to abnormal settlement of the houses.

Although the scheme provided in this paper identifies areas with anomalous values of surface deformation and utilizes a more stringent metric (anomalies of three or more times) to identify buildings with highly sensitive surface deformation, assessing the potential hazards of a building requires further access to its geologic movement and displacement characteristics of the site; mechanical parameters of the building materials; load-bearing properties; and distribution and interrelationships between the structural (columns, beams, and load-bearing walls) and non-structural (cladding, distribution) elements, as well as other building data. Therefore, the program can only be used as an auxiliary tool in the search for potentially damaging buildings.

6. Conclusions

In this paper, a semi-automated scheme is proposed for extracting valuable information from wide-area surface deformation data to identify features that are prone to deformation on urban surfaces. This scheme successfully realizes the dynamic monitoring and information extraction of surface deformation-prone features during the rainy season during the 7–20 extraordinary rainstorm by studying the Sentinel-1A SAR data and the Sentinel-2 10 m resolution land cover products transiting through Zhengzhou City, using singular spectrum analysis and the box-and-line diagram criterion. At the same time, remote sensing was fused with building contour data in Zhengzhou City to identify sensitive buildings with urban surface deformation under the influence of flooding, demonstrating the capability and application of multi-source SAR data to acquire disaster spatio-temporal information from multiple angles. The program found that (1) the change in the area of sensitive features in Zhengzhou City is similar to the change in precipitation, and the larger the precipitation, the larger the sensitive area. (2) In the rainy season with extreme precipitation, the area of highly sensitive areas accounted for a very small proportion of Zhengzhou City, and, in addition, the highly sensitive areas are more often characterized by the phenomenon of anomalous subsidence and are dominated by agricultural land and grassland, but the authenticity of the area of these two features is for reference only. (3) In the rainy season where 7–20 is located, the number of houses with highly sensitive surface deformation is 955, and the number of buildings with anomaly settlement is slightly more than that of anomaly uplift; these highly sensitive houses are mainly distributed in the old urban areas, and after sampling and comparing with the actual situation, it was found that the main characteristics of these highly sensitive houses are the age of the houses or their own unstable structure, and the credibility of this program to identify sensitive buildings with surface deformation is relatively high.

Author Contributions

Conceptualization, L.H., L.C., Q.W., J.H. and B.Y.; methodology, L.H., Q.W. and L.C.; software, L.H., B.Y. and Q.W.; validation, L.H. and J.H.; formal analysis, L.H.; investigation, L.H.; resources, J.H, L.H. and L.C.; data curation, L.H. and J.H.; writing—original draft preparation, L.H.; writing—review and editing, L.H. and L.C.; visualization, L.H.; supervision, L.C.; project administration, L.C.; funding acquisition, L.C. All authors have read and agreed to the published version of the manuscript.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Sentinel-1A (https://search.asf.alaska.edu/#/, accessed on 10 October 2023), SRTM DEM (http://srtm.csi.cgiar.org/srtmdata, accessed on 10 October 2023). Satellite precise orbit data (https://dataspace.copernicus.eu/explore-data/data-collections/sentinel-data/sentinel-1, accessed on 10 October 2023). <Environmental Geological Survey and Evaluation of Major Cities in Henan Province-Zhengzhou City Summary Report> (https://www.ngac.cn/dzzlfw_sjgl/d2d/dse/category/fileView.do?path=/dzzlfwsjgl/admin_file/core/file/toView.do?id=4028808563f8a21c01643a6b217436c4, accessed on 30 November 2009). Landcover (https://www.arcgis.com/apps/instant/media/index.html?appid=fc92d38533d440078f17678ebc20e8e2, accessed on 30 December 2021). Baidu Street View Map, (https://map.baidu.com/@12651557.08,4106835.09,21z,87t,74.85h#panoid=09028900121707121435495398I&panotype=street&heading=140.85&pitch=19.26&l=13&tn=B_NORMAL_MAP&sc=0&newmap=1&shareurl=1&pid=09028900121707121435495398I, accessed on 30 June 2021).

Acknowledgments

The authors would like to thank the editor and reviewers for their contributions to the paper.

Conflicts of Interest

The authors declare no conflict of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1. (a) Location of the study area, (b) shows the distribution of the study soil types, (c) distribution of feature distribution type distribution.

Figure 2. Principle of the outlier test method for each image element, the image time-series data. Raw data in the lower left corner in (a) are decomposed by (b) to the residual term (c), and the unreasonable value in (c) is judged as an outlier outlier.

Figure 3. The annual average distribution map of surface deformation rate in Zhengzhou City (a), as well as the distribution maps of different levels of dangerous deformation rate development (b–d). In the legend, dark green represents the high risk rate of surface deformation, while light green and red represent the medium risk rate of surface deformation. (e) The daily cumulative precipitation of 57,083 meteorological station in Zhengzhou City from 20:00 to 20:59 (with a gray background indicating the rainy season of that year).

Figure 4. Area maps (a,c,e,g) of the sampling areas of the four PS points and their corresponding cumulative changes of all PS points in the sampling frame (b,d,f,h), different colors represent the cumulative surface deformation of all PS points at different positions within the sampling area, and the box-and-line plots are the deformation rate distributions of all all PS points in the sampling frame.

Figure 5. (a) Shows the variation in the number of anomalous pixels, (b) shows the result of time-series decomposition of a PS point, which is used to analyze the phenomenon of thick heads and thick tails of the number of anomalous pixels in (a), and (c) shows the statistics of different anomalies in different features, with the red background representing the anomalous subsidence phenomenon and the blue background representing the anomalous uplift phenomenon. (d) shows the number of times surface deformation anomalies occur at different pixel positions, with 0 before the number indicating the number of abnormal subsidence times and 0 without indicating the number of abnormal uplift times.

Figure 6. (a) Shows the apparent error anomaly detection image region, (b) shows the cumulative shape variable change for all PS points in the sampling box in (a), the gray background is the rainy season period of the year, and the black boxes in (b) mark the points with apparent anomalous shape variable values during the rainy season, different colors represent the cumulative surface deformation of all PS points at different positions within the sampling area, the red line indicates the occurrence of surface deformation anomalies at the corresponding time point.

Figure 7. Distribution of highly sensitive housing.

Figure 8. (a–d) Shows the zone maps of four randomly selected highly sensitive buildings, the white dots are the PS points, the cumulative deformation values of the PS points on each building are indicated in (i,j), the gray background is the rainy season, the box marks out the obvious anomalous values, and the actual houses corresponding to (a–d) are indicated in (e–h). Situation.

Table 1

Indicators for evaluating the risk of surface deformation.

Hazard Level	Deformation Evaluation	Deformation Rate (mm/yr)
Lower risk	Weak development	≤10
Moderate risk	Medium development	10~30
High risk	intensive development	≥30

Table 2

Area and number of physical anomalies at each site.

Highly Sensitive Areas	Total Area (m)2	Crops (m)2	Rangeland (m)2	Number of Buildings
05	8800	319	80	0
04	45,600	2543.45	11,764	21
03	2,764,400	2,326,908	86,070	468
3	702,800	228,980	33,773	376
4	46,800	21,001	1476	10
5	3600	1464	0	0

References

1. Lesk, C.; Rowhani, P.; Ramankutty, N. Influence of extreme weather disasters on global crop production. Nature; 2016; 529, pp. 84-87. [DOI: https://dx.doi.org/10.1038/nature16467] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/26738594]

2. Shahid, S.; Pour, S.H.; Wang, X.; Shourav, S.A.; Minhans, A.; bin Ismail, T. Impacts and adaptation to climate change in Malaysian real estate. Int. J. Clim. Chang. Strateg. Manag.; 2017; 9, pp. 87-103. [DOI: https://dx.doi.org/10.1108/IJCCSM-01-2016-0001]

3. Zhou, M.; Guo, J.; Liu, X.; Hou, R.; Jin, X. Analysis of GNSS-Derived Tropospheric Zenith Non-Hydrostatic Delay Anomaly during Sandstorms in Northern China on 15th March 2021. Remote Sens.; 2022; 14, 4678. [DOI: https://dx.doi.org/10.3390/rs14184678]

4. Ferretti, A.; Prati, C.; Rocca, F. Permanent scatterers in SAR interferometry. IEEE Trans. Geosci. Remote Sens.; 2001; 39, pp. 8-20. [DOI: https://dx.doi.org/10.1109/36.898661]

5. Berardino, P.; Fornaro, G.; Lanari, R.; Sansosti, E. A new algorithm for surface deformation monitoring based on small baseline differential SAR interferograms. IEEE Trans. Geosci. Remote Sens.; 2002; 40, pp. 2375-2383. [DOI: https://dx.doi.org/10.1109/TGRS.2002.803792]

6. Ramirez, R.A.A.; Lee, G.-J.; Choi, S.-K.; Kwon, T.-H.; Kim, Y.-C.; Ryu, H.-H.; Kim, S.; Bae, B.; Hyun, C. Monitoring of construction-induced urban ground deformations using Sentinel-1 PS-InSAR: The case study of tunneling in Dangjin, Korea. Int. J. Appl. Earth Obs. Geoinf.; 2022; 108, 102721. [DOI: https://dx.doi.org/10.1016/j.jag.2022.102721]

7. Li, M.; Shen, Y.; Ren, Q.; Li, H. A new distributed time series evolution prediction model for dam deformation based on constituent elements. Adv. Eng. Inform.; 2019; 39, pp. 41-52. [DOI: https://dx.doi.org/10.1016/j.aei.2018.11.006]

8. Shi, X.; Pan, B. LADSDIn: LiCSAR-Based Anomaly Detector of Seismic Deformation in InSAR. IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens.; 2023; 16, pp. 4400-4413. [DOI: https://dx.doi.org/10.1109/JSTARS.2023.3272026]

9. Wang, A. Monitoring Data Fusion and Hydrogeological Interpretation of Urban Land Subsidence in Zhengzhou City. Ph.D. Thesis; Wuhan University: Wuhan, China, 2018.

10. Sun, Y.; Genton, M.G. Adjusted functional boxplots for spatio-temporal data visualization and outlier detection. Environmetrics; 2012; 23, pp. 54-64. [DOI: https://dx.doi.org/10.1002/env.1136]

11. Zhao, C.; Yang, J. A Robust Skewed Boxplot for Detecting Outliers in Rainfall Observations in Real-Time Flood Forecasting. Adv. Meteorol.; 2019; 2019, 1795673. [DOI: https://dx.doi.org/10.1155/2019/1795673]

12. Bedford, J.; Bevis, M. Greedy Automatic Signal Decomposition and Its Application to Daily GPS Time Series. J. Geophys. Res.-Solid Earth; 2018; 123, pp. 6992-7003. [DOI: https://dx.doi.org/10.1029/2017JB014765]

13. Zhao, Z.; Meng, Y. Review of remotely sensed time series data for change detection. J. Remote Sens.; 2016; 20, pp. 1110-1125.

14. Li, C.; Yang, P.; Zhang, T.; Guo, J. Periodic signal extraction of GNSS height time series based on adaptive singular spectrum analysis. Geod. Geodyn.; 2023; [DOI: https://dx.doi.org/10.1016/j.geog.2023.04.003]

15. Sahoo, S.K.; Katlamudi, M.; Gakka, U.L. Singular spectrum analysis on soil radon time series (222Rn) in Kachchh, Gujarat, India: Detection of periodic oscillations and earthquake precursors. Arab. J. Geosci.; 2020; 13, 973. [DOI: https://dx.doi.org/10.1007/s12517-020-05946-y]

16. Telesca, L.; Lovallo, M.; Cardettini, G.; Aromando, A.; Abate, N.; Proto, M.; Loperte, A.; Masini, N.; Lasaponara, R. Urban and Peri-Urban Vegetation Monitoring Using Satellite MODIS NDVI Time Series, Singular Spectrum Analysis, and Fisher–Shannon Statistical Method. Sustainability; 2023; 15, 11039. [DOI: https://dx.doi.org/10.3390/su151411039]

17. Zhou, C.; Lan, H.; Burgmann, R.; Warner, T.A.; Clague, J.J.; Li, L.; Wu, Y.; Zhao, X.; Zhang, Y.; Yao, J. Application of an improved multi-temporal InSAR method and forward geophysical model to document subsidence and rebound of the Chinese Loess Plateau following land reclamation in the Yan’an New District. Remote Sens. Environ.; 2022; 279, 113102. [DOI: https://dx.doi.org/10.1016/j.rse.2022.113102]

18. Li, X.; Zhu, Y. Precipitation Diurnal Variation in Summer in Zhengzhou from 1961 to 2017. Henan Sci.; 2019; 37, pp. 924-932.

19. Gao, C. Analysis of Ground Subsidence Patterns Andcauses in Zhengzhou City Based onInSAR Technology. Ph.D. Thesis; Henan Polytechnic University: Jiaozuo, China, 2022.

20. He, Q.; Yang, H. Correlation analysis of surface deformation and groundwater level in Zhengzhou city. Prog. Geophys.; 2017; 32, pp. 2636-2641.

21. Zhang, J. Research on the Genetic Mechanism of Land Subsidence in Zhengzhou. Ph.D. Thesis; North China University of Water Resources and Electric Power: Zhengzhou, China, 2021.