Monitoring Analysis of Urban Subsidence in Northern Henan Province Based on TS-InSAR Technology

,


Introduction
Surface subsidence represents a geological phenomenon resulting from the confluence of diverse natural and anthropogenic factors.Natural contributors to surface subsidence encompass tectonic subsidence, seismic events, volcanic activity, climatic variations, alterations in ground stress, and inherent soil consolidation processes.Conversely, anthropogenic activities such as the exploitation of subterranean fluid resources like groundwater and natural gas, extraction of solid minerals, karst collapse, and consolidation settlement associated with civil engineering endeavors in soft soil regions also contribute to surface subsidence [1][2][3][4].This process typically endures over prolonged periods, posing challenges in its timely detection, while the resultant damages are arduous to rectify [5,6], thereby impeding the urbanization process significantly.Land subsidence, characterized by the horizontal sinking of the Earth's surface, can yield severe ramifications, particularly when uneven subsidence occurs in urban or densely populated localities.Such nonuniform subsidence can inflict substantial damage upon structures, infrastructure, and pose existential threats to human lives and property.Vigilant monitoring endeavors in these areas assume paramount importance in disaster mitigation and safeguarding economic assets [7][8][9][10].Presently, land subsidence poses a threat to over 150 countries and regions worldwide, particularly densely populated urban agglomerations, such as Italy, Japan, and Pakistan, precipitating substantial economic losses [11][12][13][14][15][16][17][18].China, too, grapples with a pronounced impact from land subsidence, notably in regions like the Yangtze River Delta, the North China Plain, the Fenwei Basin, the Pearl River Delta, the Northeast Plain, the Huaihe Plain, the Jianghan Plain, the coastal plains, and mountain fault basins.As of early 2016, over 105 prefecture-level cities across 21 provinces (including municipalities directly under the Central Government) in China have experienced discernible land subsidence [4,[19][20][21][22][23], underscoring the urgent imperative for comprehensive monitoring initiatives to address this burgeoning issue.The protracted duration and widespread prevalence of land subsidence have engendered significant losses in both livelihoods and productivity [24,25].Consequently, the monitoring of the long-term spatiotemporal evolution of land subsidence assumes paramount importance in instituting remedial measures and preemptive warning systems for critical projects [26,27].Traditional surveying methods such as leveling, Global Navigation Satellite System surveying, and discrete point geodesy may be employed for high-precision ground monitoring.However, these conventional methodologies encounter substantial limitations when applied to urban land subsidence monitoring.The expansive research area and intricate terrain configurations in urban locales impose constraints on conventional surveying methods, including low spatial distribution rates, exorbitant costs, time-intensive processes, limited capacity to procure large-scale monitoring results, susceptibility to environmental variables, and inadequate reflection of continuous surface settlement information within the study area [19,[28][29][30][31].In contrast, synthetic aperture radar interferometry (InSAR) technology, conceived in the 1990s and continuously refined since then, proffers a plethora of high-precision monitoring outcomes [32][33][34][35].InSAR technology embodies characteristics, such as stable repetitive cycles, all-weather, and all-day monitoring capabilities, and has evinced efficacy in large-scale surface deformation monitoring [22,27,[36][37][38].Presently, InSAR technology finds application not only in rapid response to geological disasters such as seismic events and volcanic eruptions but also in the long-term monitoring of crustal movements, glaciers, landslides, urban and mining area land subsidence, and the stability of artificial structures [22,[39][40][41][42][43].However, the accuracy of traditional D-InSAR technology is encumbered by environmental variables such as atmospheric effects and space-time baselines [44].Gradually developed sequential InSAR technology effectively enhances monitoring accuracy and reliability by superimposing complex images to mitigate the impact of atmospheric delays and terrain effects.This advancement augments settlement accuracy from the erstwhile range of 10-20 mm to 2-3 mm, effectively surmounting the limitations of traditional D-InSAR technology.
In recent years, temporal InSAR techniques such as PS-InSAR [32,33] and small baseline subsets InSAR (SBAS-InSAR) [45] have emerged as pivotal technical modalities in surface settlement monitoring and have garnered widespread adoption both domestically and internationally.Among these techniques, the PS-InSAR technique discerns stable persistent scatterer points within the study area, particularly in artificial structures within urban environments [2,46,47].This technique evinces remarkable effectiveness in monitoring urban surface subsidence.
Drawing upon these considerations, this investigation employs PS-InSAR technology in tandem with Sentinel-1A radar image data procured from 24 ascending orbits spanning the temporal domain from March 2017 to February 2019.The amassed data undergo processing utilizing SARScape software to delineate the distributional attributes of surface settlement within the designated research domain.Recognizing the topographic disparities between the western and eastern sectors of the study area, characterized by rugged terrain and plains, respectively, it becomes apparent that the dense vegetation cover significantly influences the identification of PS points.In this experimental framework, the amplitude dispersion index methodology is applied to ascertain the PS points within the study area, thereby augmenting their spatial density.Furthermore, the mean settlement velocity across the Northern Henan Plain during the period from 2017 to 2019 is computed.A comparative assessment of secondary data layers within the study locale is executed to furnish valuable empirical support for mitigating land subsidence in this geographical expanse.

Data and Methods
2.1.Study Area and Datasets.The North Henan Plain, strategically located in the northern sector of Henan Province and proximal to the Yellow River, assumes substantial geographical significance as an integral segment of the North China Plain.It is demarcated by the Yellow River to the south, the Taihang Mountains to the west, Hebei Province to the north, and Shandong Province to the east.The topography of the western sector is distinguished by undulating hilly terrains, whereas the eastern sector is predominantly characterized by expansive plains.The overall topographical gradient exhibits a gradual elevation from west to east, with a subtle orientation toward the northeast.The focal area of this study encompasses the northern part of Henan Province, specifically incorporating the urban centers of Anyang, Puyang, and Hebi.Geographically positioned between latitudes 35°and 36.5°N and longitudes 113.5°and 115.5°E, this region epitomizes a quintessential model for resource exploitation within the plains of Henan Province.This area is delineated by a warm temperate semihumid and semiarid continental monsoon climate, manifesting cold, dry winters, and springs, alongside hot, rainy summers, and autumns, thereby engendering pronounced seasonal distinctions.Over an extended temporal average, temperatures fluctuate between 11 and 14.2°C, with an annual mean precipitation ranging from 537 to 622 mm.

2
Journal of Sensors Nonetheless, notable interannual variability is observed in the temporal and spatial precipitation patterns.Temporally, the bulk of precipitation is recorded between June and September, constituting over 70% of the annual aggregate.Spatially, precipitation distribution is primarily concentrated within the eastern plains.The geographical coordinates of the study area are elucidated in Figure 1.
The dataset employed in this study was procured from the Sentinel-1 satellite, which is administered by the European Space Agency.The Sentinel-1 constellation, consisting of two spacecraft, specifically Sentinel-1A and Sentinel-1B, was inaugurated in 2014.These satellites facilitate terrestrial monitoring by offering a 6-day revisit cycle, thus ensuring consistent and systematic data collection for analytical research purposes.Sentinel-1A, a sophisticated high-resolution synthetic aperture radar (SAR) satellite, orbits the Earth in a sun-synchronous trajectory at an elevation of 693 km.It boasts an orbital inclination of 98.18°, a period of 99 min, and a 12-day revisit interval.Outfitted with a C-band SAR apparatus, Sentinel-1A is capable of all-weather, versatile imaging across multiple polarization modes, including single polarization, dual polarization, and quad polarization, enabling the acquisition of detailed surface imagery under various atmospheric conditions.The precision of this satellite's detection capabilities extends to millimetric or even submillimetric accuracy.Sentinel-1A is operational in four principal modes: SM (strip mode), IW (Interferometric wide-field mode), EW (Extra wide-file model), and WV (Wave Mode), each designed to cater to specific observational requirements.Comprehensive details regarding the essential parameters of Sentinel-1A are enumerated in Table 1.
For this experiment, a comprehensive dataset comprising 24 radar images procured from the ascending orbit of the Sentinel-1A satellite was utilized to encompass the designated study region.This dataset extended over a period from March 14, 2017 to February 8, 2019.The chosen configuration for this analysis included the VV polarization mode and the IW operational mode, with the radar signal wavelength set at 5.56 cm.For topographical analysis within the study area, the digital elevation model (DEM) was sourced from the ASTER GDEM dataset, characterized by a 30 m resolution.This dataset is a product of the Terra satellite, an advanced Earth observation platform administered by NASA.Detailed  attributes of the Sentinel-1A radar images employed in this experimental framework are delineated in Table 2, providing a basis for further empirical analysis.The dataset acquired from the designated research area was incorporated as secondary leveling data within the "Qingfeng-Nanle Mengzhou-Wenxian Land Subsidence Monitoring Technical Service Level Monitoring Results Report."This analysis encompassed a total of 150 leveling points, of which seven were flagged for exhibiting anomalous values.The precision of the leveling network's adjustment, particularly within the Qingfeng-Nanle groundwater depression cone in the study area, is documented in Table 3.This provides a quantitative basis for assessing the network's calibration accuracy in response to subsurface water level fluctuations.

Principles of PS-InSAR Technology.
Certain locations on the terrestrial surface manifest enduring physical attributes and consistent radar backscatter, meriting their classification as "PS points."These PS points are primarily identified within man-made structures, including buildings, concrete pavements, bridges, and naturally exposed rock formations, characterized by their high degree of coherence.The PS-InSAR methodology capitalizes on the stability of these PS points across both temporal and spatial dimensions, applying differential analytical techniques to adjacent point targets.This strategy facilitates the acquisition of highly accurate ground deformation data within the specified research zone, concurrently alleviating the influence of spatial-temporal coherence disruptions and atmospheric disturbances that typically challenge traditional interferometric methods.
The operational framework of the PS-InSAR technique is predicated on the collection of N + 1 SAR images of a targeted area from a satellite.Among this collection, a singular image is selected as the referential or primary image, a decision influenced by criteria such as temporal baseline, spatial baseline, and Doppler centroid frequency.The remaining N images are then meticulously registered and resampled to conform to the pixel dimensions of the primary image.Through the execution of differential interferometric processing on the N image pairs, an interferogram is constructed.The interferometric phase for each pixel within the image pairs is delineated as follows: The interference phase of a pixel comprises several components, namely the topographic phase φ top , flat phase φ ref (also known as reference ellipsoid phase), terrain phase φ def , atmospheric delay phase φ atm , and noise phase φ noi .The flat phase arises as a systematic error phase due to the ground not serving as an ideal reference datum.This error can be mitigated by utilizing precision orbit data.Following the flattening process, the interference phase of the pixel is primarily influenced by the terrain phase.To account for this, DEM data encompassing the study area is incorporated to simulate the terrain phase.However, it should be noted that the accuracy of the DEM data introduces residual elevation error phase during this simulation.The differential interference phase associated with ground point deformation can be mathematically expressed as follows: In the given passage, the following variables are defined: φ diff represents the pixel's differential interference phase, φ res Ultimately, the surface deformation phase can be mathematically expressed as follows:

PS-InSAR Processing Flow
(1) The formulation of the connectivity graph commences with an analysis of the temporal baseline, spatial baseline, and Doppler centroid frequency.Upon meticulous examination, the SAR imagery captured on March 9, 2018 is designated as the principal master image.This selection paves the way for the establishment of master-slave data pairings between the chosen master image and additional datasets, strictly adhering to a predefined critical baseline threshold.This meticulous procedure engenders the creation of SAR data pairings and connectivity graphs, which are indispensable for executing differential interferometry.The spatiotemporal baseline of the PS-InSAR is shown in Figure 2.
(2) The procedure continues with the sequential registration of the 23 supplementary images to the primary image, followed by the generation of interferograms for each pair comprising an auxiliary image and the main image via interferometric processing.Within the purview of this experimental investigation, the amplitude dispersion index-defined as the quotient of a pixel's temporal mean amplitude and its standard deviation-was employed to pinpoint candidate PS points.This approach is particularly effective in urban environments where the amplitude standard deviation is markedly low, thus enabling the efficient identification of points of interest without necessitating the analysis of phase information.To amend the distortions introduced by flat phase and topographic phase effects, the interferogram undergoes a process of deflattening, facilitated by the integration of external DEM data.It is critical to acknowledge that the success in mitigating topographic phase effects is heavily dependent on the precision and resolution of the utilized external reference DEM.Consequently, a higher degree of accuracy and resolution in the external reference DEM significantly enhances the outcome of topographic phase correction.
(3) The initial phase of Persistent Scatterer Inversion primarily hinges on the identification of a specified quantity of PS points, with the historical phase information of these dependable individual targets (where each image pixel is representative of a singular target) being scrutinized.This examination facilitates the calculation of the velocity of displacement at the research site and the residual topographic phase, which are subsequently utilized to rectify the synthesized interferogram.A salient feature of PS-InSAR technology is its adeptness in leveraging densely arrayed scatterers to counteract the effects of signal propagation delay and correct the interferogram.(4) In the subsequent phase of Persistent Scatterer Inversion, following the correction of the interferogram, it is employed to ascertain deformation-related metrics.
After the acquisition of surface deformation parameters for each SAR image, the atmospheric phase is deduced by applying the results derived from the preliminary linear model inversion.This step is followed by atmospheric adjustments.The ultimate rate of surface deformation is ascertained through a secondary inversion process.( 5) The encoding of geographical information involves the meticulous selection of points with high coherence Journal of Sensors as PS points, given that the deformation data obtained from these points precisely mirrors the aggregate deformation experienced within the study zone.To guarantee the selection of points of exceptional quality, a coherence coefficient threshold of 0.75 was established, leading to the subsequent geocoding process that culminates in the generation of a vector file encapsulating the PS points.After an initial screening, a total of 405,427 PS points were identified within the area of study.The methodological sequence employed by the PS-InSAR technology in this investigation is illustrated in Figure 3.

PS-InSAR Land Subsidence Results
. Utilizing the PS-InSAR methodology, the temporal and spatial dynamics of surface deformation within the designated study region were delineated for the timeframe spanning March 2017 to February 2019.Nonetheless, the concentration of PS points predominantly within man-made structures, notably urban centers, as opposed to their sparse distribution in rural locales, hampers their efficacy in mirroring the settlement distribution traits of the area under study comprehensively.The distribution of PS points is shown in Figure 4.
To surmount this challenge and garner a more holistic comprehension of surface deformation, the outcomes of the data processing were subjected to further analysis employing the inverse distance weighting interpolation technique.This approach facilitated the interpolation among data points, thereby generating a continuum of deformation outcomes across the study region, as illustrated in Figure 5. Significantly, Figure 5 elucidates a discernible subsidence trend pervading the entire study area, with the phenomenon of subsidence being preeminently manifest in the central and eastern sectors, where subsidence rates fluctuate between 10 and 50 mm/year.Conversely, the western segment witnesses relatively insubstantial subsidence, with the majority of locales registering subsidence rates below 15 mm/year, and certain areas even showcasing uplift phenomena.
To probe into the settlement interrelations among the cities of Anyang, Hebi, and Puyang, two transect lines were delineated: one extending between Anyang and Puyang, and the other spanning from Hebi to Puyang, as demonstrated in Figure 5.Each transect was segmented into 100 equidistant points, thereby facilitating the derivation of settlement metrics along these vectors.Figures 6 and 7 categorically reveal that subsidence is chiefly prevalent in rural sectors, whereas the urban conurbations of Anyang, Hebi, and Puyang manifest pronounced bulging in comparison to the adjacent rural areas, engendering an inverted "U" profile.This observation is corroborated in Figure 5, which denotes the minimal subsidence rate within the entire study area to be situated in the western precinct of Anyang City.Given the prominence of Anyang, Hebi, and Puyang as principal urban centers within the study locale, accruing precise subsidence data for these metropolitan areas holds substantial importance.An exhaustive analysis of both the temporal and spatial distribution characteristics, in conjunction with the temporal evolution traits of land subsidence within the study domain, is imperative for ameliorating the adverse impacts engendered by urban land subsidence.
Throughout the observation timeframe, the surface terrain within the western sector of Anyang District exhibited notable stability, with certain locales even undergoing an elevation increase.Conversely, the subsidence anomaly predominantly manifested in the central and eastern sectors of Anyang, culminating in the emergence of three distinct subsidence epicenters within this vicinity.The rate of ground subsidence within the central depression zone of Anyang City was quantified to range between approximately 11 and 19.84 mm/year, spanning an expanse of 100.14 km 2 .A temporal analysis of settlement patterns in Anyang City

Journal of Sensors
February of the subsequent year.Meanwhile, subsidence within the eastern precincts of Anyang City exhibits a tendency toward mitigation from March to October, whereas the western section evidences a progressive elevation trend.Hebi City demonstrates a pronounced trend of overall ground settlement, with this phenomenon being markedly more acute within the tricity area.Notably, two significant areas of subsidence are discernible: one in the northeast, at the confluence of Hebi City and Xunxian County, and another in the southeast, at the boundary between Hebi City and Qi County.The subsidence rates within these locales span from 20 to 27 mm/year, enveloping an estimated territory of 30.74 km 2 .Temporal analysis of the surface subsidence sequence in Hebi City, as delineated in Figure 9, indicates that from the inception of the monitoring phase until September 2017, subsidence levels in Hebi City remained minimal, exhibiting a general tendency toward gradual elevation gain.However, from October 2017 through January 2018, the    In Puyang City, the phenomenon of land settlement predominantly manifests within the western and eastern sectors, with a settlement velocity estimated between 18 and 25.6 mm/ year, encompassing an area of roughly 162.18 km 2 .An examination of the temporal dynamics of land subsidence in Puyang City, as illustrated in Figure 10, discloses a progressive escalation in settlement rates extending from the urban core to the periphery, with a significant concentration around the urban fringe.During the interval from April to August 2017, the settlement activity within these regions exhibited notable stability.Nonetheless, commencing in January 2018, the onset of settlement was observed in the western portion of Puyang City, with the affected area broadening progressively until a cumulative settlement of 30 mm was recorded.The period between February and October 2018 saw a gradual return to stability in the settlement patterns within Puyang City, accompanied by a minor elevation gain in the central locality.Approaching the conclusion of the monitoring timeline, there was a resurgence in the rate of settlement on the outskirts of Puyang City.

PS-InSAR Accuracy
Assessment.The PS-InSAR monitoring outcomes postexperimental intervention were juxtaposed with terrestrial subsidence level data derived from groundwater funnel observations in the North China Plain (Henan segment), specifically along the Qingfeng-Nanle leveling route.Both datasets demonstrated a point-based distribution; however, there was a notable density difference with PS points being densely clustered in spatial arrangement, as opposed to the more dispersed distribution of leveling points, leading to a spatial incongruity between the two sets of data.This divergence suggests a variance in resolution.To mitigate this discrepancy, an averaging procedure was applied to the PS-InSAR data, using the leveling points as referential centroids.This entailed computing the mean value of points with high coherence within a designated proximity as the PS-InSAR monitoring outcome.In instances where PS points were absent within the stipulated range, the precision validation for the corresponding water point was foregone.Consequently, a total of 96 PS-InSAR monitoring outcomes were compiled.
Subsequently, a scatter diagram was constructed to facilitate a comparison between the leveling data and the PS-InSAR findings.Figure 11 delineates the trend line deduced via SSA, revealing that both the leveling and PS-InSAR datasets exhibit parallel deformation trends, with their settlement metrics showing substantial congruence.Nonetheless, the settlement estimations derived from PS-InSAR monitoring were observed to surpass those ascertained through leveling.The root mean square error (RMSE) and standard deviation associated with the PS-InSAR results were recorded at AE12.9 and AE13.31 mm, respectively, underscoring the reliability of the PS-InSAR monitoring methodology.

Discussion
The interrelation between groundwater levels and surface deformation constitutes a focal point of considerable scholarly interest.Prior investigations have robustly corroborated the linkage between the phenomena of land subsidence in the North China Plain and the overextraction of groundwater resources [48,49].The practice of excessive groundwater  withdrawal precipitates a sustained diminution in the groundwater table.Invoking the effective stress principle as proposed by Terzaghi, it is observed that a reduction in the pore water pressure ensues from the descending water table.This reduction in pore pressure negates its contribution to the structural integrity and deformational capacity of the soil matrix, thereby amplifying the effective stress and leading to enhanced ground settlement.
The research and evaluation pertaining to the sustainable exploitation of groundwater within the North China Plain, with a particular focus on Henan Province, reveal that the areas under study are predominantly located within zones of intensive groundwater extraction.Table 4 illustrates the alterations in groundwater levels within the research locale.Within this context, Anyang and Hebi are identified as regions of superficial groundwater overextraction, whereas Puyang is distinguished as an area of profound confined water overextraction.Significantly, the geographic positioning of Puyang coincides with the Qingfeng-Nanle depression cone, indicative of a decrease in superficial groundwater levels.Crucially, the geographical distribution of these locales is consistent with the incidences of land subsidence identified within the scope of the study.Second, the ramifications of project construction on urban land subsidence emerge as a pivotal consideration in the ongoing urbanization narrative.The implications of such construction activities can be bifurcated into two distinct categories: localized land subsidence, attributable to the concentrated construction of edifices on the surface, and regional land subsidence, engendered by expansive engineering projects executed subterraneously.A juxtaposition of this study's outcomes with remote sensing imagery of the designated area reveals that, throughout the monitoring timeframe, urban engineering endeavors significantly influenced Anyang and Hebi, as illustrated in Figure 12.Conversely, the impact of engineering activities in the Puyang region was markedly less pronounced, verging on negligible.To a certain degree, these construction projects play a contributory role in the initiation or amplification of terrestrial subsidence.

Conclusions
The current investigation harnesses the capabilities of PS-InSAR technology, leveraging Sentinel-1A radar imagery across 24 research zones, to scrutinize land subsidence in urban locales across the northern plains of Henan Province.The objective is to delineate the spatiotemporal distribution patterns of surface deformation within the designated area from March 2017 to February 2019.This endeavor involves a comparative analysis with secondary leveling data collated over the identical timeframe.The precision of the PS-InSAR methodology applied in this study is rigorously validated, facilitating an in-depth examination of deformation causatives, integrating geological records, time series deformation analyses, groundwater usage statistics, and the impact of urban development, among other variables.The study's outcomes elucidate that: (1) The results derived from the application of PS-InSAR technology align closely with the leveling data, revealing a pronounced propensity for land subsidence across the study domain, particularly within its eastern sectors, where the maximum annual subsidence rate reaches up to 50 mm.Seasonal analysis indicates a significant incidence of subsidence during the autumn and winter months.Comparative validation of leveling and experimental data unveils the RMSE and standard deviation of AE12.9 mm and AE13.31 mm, respectively, underscoring the robustness and reliability of the detection results achieved through PS-InSAR technology.(2) Further analysis of geological data from the study area establishes a tangible link between land subsidence and two predominant factors: the overextraction of groundwater and extensive urban construction activities.The spatial congruence of subsidence epicenters within areas of groundwater overexploitation mirrors the general distribution of the research zone.A consistent pattern of groundwater overdraft is observed, precipitating a reduction in aquifer levels and, consequently, diminishing the load-bearing capacity of the foundational strata.Moreover, the influence of engineering construction on land subsidence is significant.The widespread execution of urban engineering initiatives modifies the foundational stress conditions, thereby exerting a contributory effect on ground settlement.

FIGURE 1 :
FIGURE 1: Map of geographical location of the study area.

FIGURE 2 :
FIGURE 2: PS-InSAR spatiotemporal baseline distribution: (a) indicates the time baseline and (b) indicates the space baseline.

FIGURE 6 :
FIGURE 6: Distribution characteristics of land subsidence on the profile line from Anyang to Puyang.

FIGURE 7 :
FIGURE 7: Distribution characteristics of land subsidence on the profile line from Hebi to Puyang.

FIGURE 11 :
FIGURE 11: Comparison diagram of PS-InSAR monitoring results and leveling results based on SSA smoothing.

TABLE 3 :
Adjustment accuracy of leveling net in study area.Accidental mean error Mean error of height difference Mean elevation error Actual closing error/limit error Length of leveling net phase, φ l def represents the linear deformation phase, φ ref error represents the residual phase after correcting for the reference ellipsoid phase, φ nl def represents the nonlinear deformation phase, and φ δh represents the elevation error phase.The residual phase is subjected to spatial filtering through low-pass filtering to remove high-frequency noise phase.The nonlinear deformation phase, denoted by φ nl def , is identified using time-domain filtering.

TABLE 4 :
Changes of urban groundwater level in northern plain of Henan Province (unit (m)).