InSAR ground deformation multi-source vertical contribution analysis method based on deformation signal spectrum decomposition

Abstract

The application discloses an InSAR surface deformation multi-source vertical contribution analytical method based on deformation signal spectrum decomposition and belongs to the technical field of surface deformation monitoring. The steps are as follows: first, long-time SAR images are acquired, vertical deformation data and related rates and deformation variables are obtained through SBAS-InSAR technology processing; then, a plurality of auxiliary conditions are fused, regional spatial partitioning is realized through K-Means clustering; subsequently, CEEMDAN decomposition and mean threshold multi-scale reconstruction are performed on deformation signals of each partition, high-frequency items, low-frequency items and trend items are obtained; finally, ICA blind source separation and variance contribution rate analysis are performed, and shallow, medium and deep surface deformation information is inversed through directional reorganization; the K-Means clustering algorithm is adopted to perform spatial partitioning on the regional deformation field, the homogeneity and precision of subsequent signal processing are improved; and the CEEMDAN decomposition and multi-scale reconstruction method is introduced, adaptive and fine scale separation of deformation time series signals is realized; the blind source separation technology is applied to multi-scale component reorganization, and the physical attribution of the vertical position of the subsidence source is realized.

Description

Technical Field

[0001] This invention belongs to the field of surface deformation monitoring technology, specifically, it relates to an InSAR method for analyzing the multi-source vertical contribution of surface deformation based on deformation signal spectral decomposition. Background Technology

[0002] Land subsidence, a slowly evolving geological hazard caused by human activities such as groundwater over-extraction, underground resource exploitation, and large-scale engineering construction, has become a critical issue threatening the safety of urban buildings, the stable operation of major infrastructure, and the sustainable development of the regional ecological environment, driven by accelerated urbanization and industrialization. Its core hazards include: long-term subsidence leading to continuous loss of surface elevation, causing road cracking, building tilting, damage to underground pipelines, and even exacerbating flood risks, resulting in irreversible negative impacts on economic and social development. Therefore, achieving high-precision monitoring of the land subsidence process and clarifying the deformation contribution mechanism at different depths are core prerequisites for formulating scientific prevention and control strategies and predicting subsidence development trends.

[0003] Traditional ground subsidence monitoring technologies mainly rely on leveling and the Global Positioning System (GPS). While leveling can achieve single-point monitoring with millimeter-level accuracy, it is limited by its "point-distribution" characteristic, making it unable to cover large areas. It also suffers from drawbacks such as long retesting cycles and high costs in terms of manpower and resources, making it difficult to capture the spatiotemporal dynamics of subsidence. Although GPS technology can achieve regional monitoring, it is affected by satellite signal obstruction (such as urban building clusters and vegetation-covered areas), resulting in a low density of monitoring points. Furthermore, its vertical monitoring accuracy is lower than that of leveling, failing to meet the requirements for refined analysis.

[0004] The emergence of time-series synthetic aperture radar interferometry has provided a breakthrough solution for land subsidence monitoring. This technology, relying on the all-weather, wide-area coverage capability of synthetic aperture radar and combined with a "small baseline set" data processing strategy, can achieve millimeter-level precision in monitoring surface deformation and acquire high-density, long-term time-series deformation data. It effectively overcomes the shortcomings of traditional technologies, such as "sparse data points, long cycles, and high costs," and has become the mainstream technology for regional land subsidence monitoring.

[0005] However, SBAS-InSAR technology still has certain shortcomings: the surface deformation signals it acquires are essentially the superimposed and coupled results of deformations at different depths and with different mechanical mechanisms underground. Specifically, ground subsidence in a region usually includes three types of vertical contributions: first, rapid, short-period deformation in the shallow layer (e.g., 0-100m) caused by seasonal rainfall and short-term engineering activities (e.g., foundation pit excavation, temporary loads); second, elastic-inelastic compressive deformation in the aquifer system in the middle layer (e.g., 100-300m) due to the interannual extraction and recharge cycle of groundwater; and third, secondary consolidation (creep) deformation in the thick cohesive soil in the deep layer (e.g., below 300m) due to long-term stress. These deformation signals, originating from different depths and possessing elastic, plastic, and viscoelastic characteristics, are mixed together in the spatiotemporal domain and are ultimately observed by InSAR technology as a single surface line-of-sight deformation sequence, forming the "mixed signal" problem.

[0006] While existing technologies attempt to decompose InSAR time-series signals using signal processing methods (such as Empirical Mode Decomposition (EMD) and Independent Component Analysis (ICA), most remain at the level of "planarization." In actual processing, they can only separate periodic terms (such as seasonal fluctuations) and trend terms (such as long-term subsidence) in the time series, failing to establish an objective and quantitative correlation between the mathematically decomposed components and specific geological compression layers (shallow, intermediate, and deep). For example, traditional EMD decomposition is susceptible to mode mixing, making it difficult to accurately distinguish deformation signals of different periods; and single ICA technology lacks spatial partitioning and multi-scale capabilities. The reconstruction support cannot achieve vertical attribution of deformation signals. This current situation of "emphasizing signal decomposition and neglecting geological correlation" may lead to the inability of existing technologies to effectively isolate and quantify the contribution ratio of each vertical layer to the total settlement. This may restrict the accurate identification of the main control factors of settlement (such as determining whether it is the influence of shallow engineering or the over-extraction of deep groundwater), the scientific prediction of future trends (such as distinguishing between reversible elastic deformation and irreversible plastic deformation), and the formulation of targeted treatment measures (such as targeted control of mid-level groundwater extraction or reinforcement of shallow soil). The refined research and prevention technology of ground settlement need to be further improved. Summary of the Invention

[0007] To address the problems raised in the existing background technology, this invention provides an analytical method for multi-source vertical contribution of InSAR surface deformation based on deformation signal spectral decomposition.

[0008] To solve the above problems, the present invention adopts the following technical solution.

[0009] The method for analyzing the multi-source vertical contribution of InSAR surface deformation based on deformation signal spectral decomposition is as follows:

[0010] S1. Acquire long-term synthetic aperture radar (SAR) image data covering the study area, with time points t0, t1, ..., tn; process the image data using SBAS-InSAR technology, extract the time-series deformation of the ground target points along the radar line of sight, combine it with the radar incident angle to convert it into vertical deformation, and obtain the annual average deformation rate V' of each target point and the cumulative deformation D0, D1, ..., Dn relative to the stable reference point at each observation time;

[0011] S2, at the aforementioned average annual deformation rate Using the core feature as the basis, and integrating meteorological and hydrogeological auxiliary conditions, the K-Means clustering algorithm is used to spatially partition the set of deformation points in the region into K categories with similar deformation features and spatiotemporal patterns.

[0012] S3. For each spatial partition, perform adaptive noise-complete ensemble empirical mode decomposition on the temporal deformation signals of all deformation points within the partition, decomposing the non-stationary deformation sequence of each point into a series of intrinsic mode functions and a residual term; employ a multi-scale reconstruction method based on mean thresholding to calculate the pre-... indivual The mean amplitude of the components is used as a threshold. IMFs with amplitudes greater than this threshold are merged into high-frequency terms, while those with amplitudes less than or equal to this threshold are merged into high-frequency terms. The terms are merged into low-frequency terms, and the residual terms are used as trend terms.

[0013] S4. Apply independent component analysis to the three component datasets of high-frequency term, low-frequency term and trend term to perform blind source separation and obtain multiple independent components; analyze the variance contribution rate of each independent component, and perform targeted recombination of specific independent components from high-frequency term, low-frequency term and trend term to respectively retrieve shallow surface deformation, middle surface deformation and deep surface deformation information.

[0014] Preferably, .

[0015] Furthermore, the SAR image data acquired in S1 is a continuous observation image covering the same study area with a time span of not less than one year.

[0016] Preferably, the meteorological auxiliary conditions in S2 include the average annual precipitation and average annual evapotranspiration of the study area during the same period, and the hydrogeological auxiliary conditions include data such as the thickness of the first compressible layer, the thickness of the second compressible layer, the thickness of the third compressible layer, the variation of the first confined groundwater level, the variation of the second confined groundwater level, and the variation of the third confined groundwater level in the study area.

[0017] Furthermore, the multi-scale reconstruction method based on the mean threshold in S3 specifically includes the following steps:

[0018] S301. Data Matrix Construction: Represent the InSAR deformation time series data within the partitions after K-Means clustering as a matrix. ,

[0019]

[0020] in This indicates the number of deformation points within the partition. The length of the time series, matrix elements Indicates the first The deformation point at the first Deformation at each time point;

[0021] S302, CEEMDAN Adaptive Decomposition: For the matrix Each row of timing signals in Performing CEEMDAN decomposition decomposes the original signal into a set of intrinsic mode functions and a sum of residual terms:

[0022]

[0023] in Indicates the first One eigenmode function Represents the residual term;

[0024] S303. Multi-scale reconstruction based on amplitude mean: For the decomposition result of each deformation point, calculate its preceding value. indivual The average amplitude of the component As an adaptive threshold:

[0025]

[0026] The signal is then reconstructed based on this threshold, including the high-frequency term. From amplitude greater than of Composed of multiple components:

[0027]

[0028] low frequency items From amplitude less than or equal to of Composed of multiple components:

[0029]

[0030] Trend Item For residual terms :Right now .

[0031] Furthermore, the blind source separation and vertical layering information extraction in S4 specifically include the following steps:

[0032] S401, ICA Input Matrix Construction: Set the high-frequency terms of all deformation points within the partition. Low-frequency term set Trend Item Set As three independent input data matrices ,in This represents the total number of deformation points within the partition.

[0033] S402, Signal Centering and Whitening: Centering is performed on each input matrix, i.e. ,in For matrix The mean vector over time; for the centered data Perform whitening treatment, that is ,in and They are covariance matrix The eigenvector matrix and the diagonal matrix of eigenvalues;

[0034] S403. Independent Component Extraction: The FastICA algorithm is used to solve the separation matrix by maximizing the negative entropy. From the whitening signal Extracting independent components, namely ,in The separated independent component matrix; updated iteratively. Until convergence, to maximize the negative entropy objective function. :

[0035]

[0036] in Since it is not a quadratic function, take , It is a constant, and , It is a standard Gaussian random variable;

[0037] S404. Vertical Stratified Quantitative Reorganization Based on Variance Contribution: Calculating the variance contribution rate of each independent component obtained after ICA separation of the three component datasets (high-frequency, low-frequency, and trend term). :

[0038]

[0039] in This is the variance calculation function. This represents the total number of independent components in the dataset.

[0040] Based on the typical Quaternary stratigraphic structure, the geological significance of each component is analyzed, and quantitative recombination is performed according to preset rules to obtain shallow, middle and deep surface deformation information.

[0041] Furthermore, the typical Quaternary stratigraphic structure, from top to bottom, includes: a shallow unconfined aquifer mainly composed of loose silt and fine sand; a middle layer mainly composed of alternating layers of medium and fine sand and silty clay; and a deep confined aquifer system mainly composed of thick layers of cohesive soil and dense gravel.

[0042] Furthermore, the preset recombination rule is as follows:

[0043] Shallow surface deformation information = the component with the largest variance contribution rate in the high-frequency term + the component with the smallest variance contribution rate in the low-frequency term + the component with the smallest variance contribution rate in the trend term;

[0044] Deep surface deformation information = the component with the smallest variance contribution rate in the high-frequency term + the component with the largest variance contribution rate in the low-frequency term + the component with the largest variance contribution rate in the trend term;

[0045] Information on mid-level surface deformation = the component with the second-highest variance contribution rate in the high-frequency terms + the component with the second-highest variance contribution rate in the low-frequency terms + the component with the second-highest variance contribution rate in the trend terms.

[0046] Preferably, when processing SAR image data using SBAS-InSAR technology in S1, it is implemented using Sarproz, GAMMA, or MintPy synthetic aperture radar interferometry software, and the processing includes image stitching, image registration, interferogram generation, phase unwrapping, and geocoding steps.

[0047] Furthermore, the calculations of CEEMDAN decomposition, multi-scale reconstruction based on mean threshold, and independent component analysis are implemented using Matlab or Python programming tools, and the numerical precision used in the calculation process is retained to three decimal places.

[0048] Preferably, the method further includes step S5: using stratified monitoring data and leveling monitoring data to verify the accuracy of the shallow, middle and deep surface deformation information retrieved. When the verification error is less than 0.5 mm / year, the analytical results of the method are deemed valid.

[0049] Beneficial effects

[0050] Compared with the prior art, the beneficial effects of the present invention are as follows:

[0051] (1) The present invention uses the K-Means clustering algorithm to spatially partition the regional deformation field, which effectively improves the homogeneity and accuracy of subsequent signal processing.

[0052] (2) This invention introduces CEEMDAN decomposition and a multi-scale reconstruction method based on mean threshold, which realizes adaptive and refined scale separation of deformed time series signals.

[0053] (3) This invention applies blind source separation technology to multi-scale component reconstruction, realizing the physical attribution of the vertical position of sedimentation source.

[0054] Figure 1 This is a flowchart illustrating the InSAR surface deformation multi-source vertical contribution analysis method based on deformation signal spectral decomposition proposed in this invention.

[0055] Figure 2 This is a schematic diagram of the SBAS-InSAR process for the InSAR multi-source vertical contribution analysis method for surface deformation based on deformation signal spectral decomposition proposed in this invention.

[0056] Figure 3 The cumulative deformation time series of surface target points obtained by SBAS-InSAR monitoring;

[0057] Figure 4 The sequence of high-frequency, low-frequency, and trend terms reconstructed by MTR after CEEMDAN decomposition;

[0058] Figure 5 The shallow deformation time series, intermediate deformation time series, and deep deformation time series are obtained by separating high-frequency terms, low-frequency terms, and trend terms using ICA and then recombining them according to physical mechanisms. Detailed Implementation

[0059] To make the objectives, technical solutions, and advantages of the embodiments of this application clearer, the technical solutions of the embodiments of this application will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are some embodiments of this application, but not all embodiments. Generally, the components of the embodiments of this application described and shown in the accompanying drawings can be arranged and designed in various different configurations.

[0060] Therefore, the following detailed description of the embodiments of this application provided in the accompanying drawings is not intended to limit the scope of the claimed application, but merely to illustrate selected embodiments of the application. All other embodiments obtained by those skilled in the art based on the embodiments of this application without inventive effort are within the scope of protection of this application.

[0061] Example 1:

[0062] like Figures 1-5 As shown, the implementation steps of the InSAR multi-source vertical contribution analysis method for surface deformation based on deformation signal spectral decomposition are as follows:

[0063] S1. Acquire long-term synthetic aperture radar (SAR) image data covering the study area, with time points t0, t1, ..., tn; process the image data using SBAS-InSAR technology, extract the time-series deformation of surface target points along the radar line of sight, combine it with the radar incident angle to convert it into vertical deformation, and obtain the annual average deformation rate V' of each target point and the cumulative deformation D0, D1, ..., Dn relative to the stable reference point at each observation time, where D0 = 0;

[0064] S2, with an average annual deformation rate Using the core feature as the basis, and integrating meteorological and hydrogeological auxiliary conditions, the K-Means clustering algorithm is used to spatially partition the set of deformation points in the region into K categories with similar deformation features and spatiotemporal patterns.

[0065] S3. For each spatial partition, perform adaptive noise-complete ensemble empirical mode decomposition on the temporal deformation signals of all deformation points within the partition, decomposing the non-stationary deformation sequence of each point into a series of intrinsic mode functions and a residual term; employ a multi-scale reconstruction method based on mean thresholding to calculate the pre-... indivual The mean amplitude of the components is used as a threshold. IMFs with amplitudes greater than this threshold are merged into high-frequency terms, while those with amplitudes less than or equal to this threshold are merged into high-frequency terms. The terms are merged into low-frequency terms, and the residual terms are used as trend terms, where m is the total number of IMFs obtained from the decomposition.

[0066] S4. Apply independent component analysis to the three component datasets of high frequency, low frequency and trend terms to perform blind source separation and obtain multiple independent components. Analyze the variance contribution rate of each independent component and perform targeted recombination of specific independent components from high frequency, low frequency and trend terms to retrieve shallow surface deformation, middle surface deformation and deep surface deformation information respectively.

[0067] S5: The accuracy of the shallow, middle and deep surface deformation information obtained by using stratified benchmark monitoring data and leveling monitoring data is verified. When the verification error is less than 0.5 mm / year, the analytical results of this method are deemed valid.

[0068] The SAR image data acquired in S1 are continuous observation images covering the same study area with a time span of no less than one year.

[0069] The meteorological auxiliary conditions in S2 include the average annual precipitation and average annual evapotranspiration in the study area during the same period. The hydrogeological auxiliary conditions include the total thickness of the first compressible layer, the thickness of the second compressible layer, the thickness of the third compressible layer, the dynamic change amplitude of the first confined groundwater level, the change amplitude of the second confined groundwater level, and the change amplitude of the third confined groundwater level in the study area.

[0070] The multi-scale reconstruction method based on mean thresholding in S3 specifically includes the following steps:

[0071] S301. Data Matrix Construction: Represent the InSAR deformation time series data within the partitions after K-Means clustering as a matrix. ,

[0072]

[0073] in This indicates the number of deformation points within the partition. The length of the time series, matrix elements Indicates the first The deformation point at the first Deformation at each time point;

[0074] S302, CEEMDAN Adaptive Decomposition: For the matrix Each row of timing signals in Performing CEEMDAN decomposition decomposes the original signal into a set of intrinsic mode functions and a sum of residual terms:

[0075]

[0076] in Indicates the first One eigenmode function Represents the residual term;

[0077] S303. Multi-scale reconstruction based on amplitude mean: For the decomposition result of each deformation point, calculate its preceding value. indivual The average amplitude of the component As an adaptive threshold:

[0078]

[0079] The signal is then reconstructed based on this threshold, including the high-frequency term. From amplitude greater than of Composed of multiple components:

[0080]

[0081] low frequency items From amplitude less than or equal to of Composed of multiple components:

[0082]

[0083] Trend Item For residual terms :Right now .

[0084] The blind source separation and vertical layering information extraction in S4 specifically include the following steps:

[0085] S401, ICA Input Matrix Construction: Set the high-frequency terms of all deformation points within the partition. Low-frequency term set Trend Item Set As three independent input data matrices ,in This represents the total number of deformation points within the partition.

[0086] S402, Signal Centering and Whitening: Centering is performed on each input matrix, i.e. ,in For matrix The mean vector over time; for the centered data Perform whitening treatment, that is ,in and They are covariance matrix The eigenvector matrix and the diagonal matrix of eigenvalues;

[0087] S403. Independent Component Extraction: The FastICA algorithm is used to solve the separation matrix by maximizing the negative entropy. From the whitening signal Extracting independent components, namely ,in The separated independent component matrix; updated iteratively. Until convergence, to maximize the negative entropy objective function. :

[0088]

[0089] in Since it is not a quadratic function, take , It is a constant, and , It is a standard Gaussian random variable;

[0090] S404. Vertical Stratified Quantitative Reorganization Based on Variance Contribution: Calculating the variance contribution rate of each independent component obtained after ICA separation of the three component datasets (high-frequency, low-frequency, and trend term). :

[0091]

[0092] in This is the variance calculation function. This represents the total number of independent components in the dataset.

[0093] Based on the typical Quaternary stratigraphic structure, the geological significance of each component is analyzed, and quantitative recombination is performed according to preset rules to obtain shallow, middle and deep surface deformation information.

[0094] A typical Quaternary stratigraphic structure, from top to bottom, includes: a shallow unconfined aquifer mainly composed of loose silt and fine sand; a middle layer of slightly confined to confined aquifer mainly composed of alternating layers of medium and fine sand and silty clay; and a deep confined aquifer mainly composed of thick layers of cohesive soil and dense gravel.

[0095] The default recombination rules are:

[0096] Shallow surface deformation information = the component with the largest variance contribution rate in the high-frequency term + the component with the smallest variance contribution rate in the low-frequency term + the component with the smallest variance contribution rate in the trend term;

[0097] Deep surface deformation information = the component with the smallest variance contribution rate in the high-frequency term + the component with the largest variance contribution rate in the low-frequency term + the component with the largest variance contribution rate in the trend term;

[0098] Information on mid-level surface deformation = the component with the second-highest variance contribution rate in the high-frequency terms + the component with the second-highest variance contribution rate in the low-frequency terms + the component with the second-highest variance contribution rate in the trend terms.

[0099] When using SBAS-InSAR technology to process SAR image data in S1, it is achieved through Sarproz, GAMMA, or MintPy synthetic aperture radar interferometry software, and the processing includes image stitching, image registration, interferogram generation, phase unwrapping, and geocoding steps.

[0100] The calculations for CEEMDAN decomposition, multi-scale reconstruction based on mean threshold, and independent component analysis were performed using Matlab or Python programming tools, and the numerical precision used in the calculations was retained to three decimal places.

[0101] Example 2:

[0102] This embodiment takes the typical land subsidence area of ​​the North China Plain—the southern Hebei Plain—as the research object. The method of this invention is used to realize the analysis of the multi-source vertical contribution of surface deformation in this area. The specific implementation process is as follows. All calculation steps are completed by Matlab R2023a or Python 3.9, and the numerical precision is uniformly retained to three decimal places.

[0103] I. Pre-implementation preparations:

[0104] (I) Overview of the study area:

[0105] The study area is located in the southern plain of Hebei Province, covering a total area of ​​approximately 2800 km². It has a temperate monsoon climate with an average annual precipitation of 520-650 mm and an average annual evapotranspiration of 1200-1400 mm. The Quaternary strata are well-developed and can be divided into three layers from top to bottom: the shallow layer (0-100 m) is a shallow aquifer composed of loose silt and fine sand; the middle layer (100-300 m) is a slightly confined to confined aquifer composed of interbedded medium and fine sand and silty clay; and the deep layer (>300 m) is a deep confined aquifer system composed of thick cohesive soil and dense gravel. Groundwater over-extraction exists in the region, with an average annual ground subsidence of 5-30 mm. The method of this invention is needed to analyze the subsidence contribution of different vertical strata.

[0106] (II) Data Preparation:

[0107] SAR image data: Sentinel-1A rising orbit SAR data covering the study area from January 3, 2018 to December 20, 2022 were acquired, totaling 148 scenes, with a spatial resolution of 5m (azimuth) × 20m (range).

[0108] Supporting data:

[0109] Meteorological data: Monthly precipitation and evapotranspiration data for the study area from 2018 to 2022 (sourced from the China Meteorological Administration's surface meteorological observation stations), from which the five-year average annual precipitation and five-year average cumulative evapotranspiration were calculated;

[0110] Hydrogeological data: Based on data from 52 boreholes in the study area, spatial distribution maps of the thickness of the first compressible layer (14-67m), the second compressible layer (58-250m), and the third compressible layer (128-611m) were obtained by ArcGIS Kriging interpolation, along with monthly average water level data from 133 shallow, medium, and deep groundwater monitoring wells;

[0111] Validation data: Monthly monitoring data from 2018 to 2022 from 3 stratified monitoring stations (monitoring shallow, middle and deep layer deformation respectively, with a monitoring accuracy of 0.1 mm) and 15 leveling monitoring points (monitoring accuracy of 0.3 mm) in the study area.

[0112] II. Step S1: InSAR Vertical Deformation Data Acquisition:

[0113] The SBAS-InSAR technology is used to process SAR image data. The specific workflow is implemented using Gamma software, and the steps are as follows:

[0114] 1. Image registration: Radiometric calibration and geometric correction were performed on 148 SAR images. Based on the digital elevation model of the study area, the topographic phase was removed. In GAMMA software, using the precise orbit data provided by ESA, all auxiliary images were accurately registered to the main image and unified into the same geometric reference frame.

[0115] 2. Construction of Short Baseline Interferometric Network and Differential Interferometry: For the registered SAR dataset, a time baseline threshold of 36 days was set to construct a short baseline interferometric network. Differential interferometry was performed on each interferometric pair in the network to generate differential interferograms. During this process, SRTM DEM data was used to simulate and remove terrain phase and flat terrain phase. To improve the signal-to-noise ratio and computational efficiency, all interferograms were subjected to 20×4 (azimuth × range) multi-look processing.

[0116] 3. Coherent point target identification and phase unwrapping: Initial candidate points are selected based on the amplitude deviation index threshold (0.75). Then, based on the phase stability between adjacent points, reliable permanent and distributed scatterers are selected as the final observation points by setting parameters such as weed_standard_dev to 1.0 and weed_max_noise to 2.3. Subsequently, the "3D-quick" algorithm and a 200-meter unwrapping window are used to unwrap the phase of the selected point targets and restore their absolute phase.

[0117] 4. Temporal Deformation Inversion and Residual Phase Correction: For the unwrapped phase, an observation equation is established in the spatiotemporal dimension, and the deformation rate and elevation error of each point are calculated using the singular value decomposition algorithm. To suppress residual errors in the phase, a spatiotemporal filter (time window: 60 days, spatial window: 300 meters) is applied to separate and remove the atmospheric delay phase, which mainly manifests as a random signal.

[0118] 5. Deformation Transformation and Extraction: Extract the time-series deformation of surface target points along the radar line of sight, combine it with the radar incident angle, and then use the formula... ( (where the angle of incidence is the angle of incidence), which is then converted into vertical deformation;

[0119] 6. Results Output: Obtain the annual average deformation rate of 2,088,089 coherent target points in the study area. And the cumulative deformation variables of each target point from January 3, 2018 to December 20, 2022, for a total of 148 observation times. The time series of cumulative deformation variables for some target points are as follows: Figure 3 As shown.

[0120] III. Step S2: Spatial partitioning based on K-Means:

[0121] Feature data construction: Using the "annual average deformation rate V', thickness of the first compressible layer, thickness of the second compressible layer, thickness of the third compressible layer, annual average precipitation, annual average evapotranspiration, first confined water level variation, second confined water level variation, and third confined water level variation" of each coherent target point as clustering features, a multidimensional feature matrix is ​​formed, and the feature data is standardized.

[0122] K value determination: The "elbow rule" and DB index are used to evaluate the clustering effect. When K=4, the DB index is the smallest and the clustering effect is the best.

[0123] Spatial partitioning: The K-Means clustering algorithm was used to cluster the feature matrix, dividing the study area into four partitions (R1-R4) with similar deformation characteristics and spatiotemporal patterns.

[0124] R1 zone (740,985 points): average annual settlement rate -39.28~3.31 mm / year, average annual precipitation 508.81~535.67 mm, thickness of the first compressible layer 20.08~61.54 m, thickness of the second compressible layer 73.42~250.19 m, thickness of the third compressible layer 162.74~565.85 m, variation of the first confined water level 3.71~18.40 m, variation of the second confined water level 3.14~25.42 m, variation of the third confined water level 2.87~28.59 m;

[0125] Zone R2 (149,956 points): Annual average settlement rate -69.32~6.17 mm / year, annual average precipitation 473.93~540.63 mm, thickness of the first compressible layer 19.01~66.52 m, thickness of the second compressible layer 108.05~247.23 m, thickness of the third compressible layer 341.01~610.87 m, variation of the first confined water level 7.5~20.86 m, variation of the second confined water level 3.93~51.19 m, variation of the third confined water level 2.56~52.96 m;

[0126] R3 Zone (901,650 points): Annual average settlement rate -48.81~-9.51 mm / year, annual average precipitation 479.36~533.53 mm, thickness of the first compressible layer 13.54~56.74 m, thickness of the second compressible layer 106.35~223.38 m, thickness of the third compressible layer 246.08~560.13 m, variation of the first confined water level 1.82~14.47 m, variation of the second confined water level 4.16~45.48 m, variation of the third confined water level 3.49~38.43 m;

[0127] R4 area (295,498 points): average annual settlement rate -28.02~11.28 mm / year, average annual precipitation 424.17~526.40 mm, thickness of the first compressible layer 26.25~45.14 m, thickness of the second compressible layer 143.71~190.75 m, thickness of the third compressible layer 262.66~542.76 m, variation of the first confined water level 0.95~5.61 m, variation of the second confined water level 2.97~44.71 m, variation of the third confined water level 3.59~24.96 m. [1]

[0128] IV. Step S3: CEEMDAN decomposition and MTR multi-scale reconstruction:

[0129] Taking region R1 (where deformation is most significant) as an example, the time-series deformation signals of 2,088,089 deformation points within the region are selected from the end-of-month data to form a time series, which is then processed as follows:

[0130] S301: Data Matrix Construction:

[0131] The vertical deformation time series data of region R1 is represented as a matrix. , where matrix elements Indicates the first The deformation point at the first Vertical deformation of deformation at each time point;

[0132] S302: CEEMDAN adaptive decomposition of time-series signals in each row of matrix X CEEMDAN decomposition was performed, with a noise level of 0.2 and an ensemble number of 100. The decomposition yielded m = 8 intrinsic mode functions (IMF1~IMF8) and 1 residual term R(t), satisfying the formula: Among them, IMF1~IMF3 have higher frequencies, while IMF4~IMF8 have lower frequencies, and the residual term R(t) shows a monotonic trend.

[0133] S303: Multi-scale reconstruction threshold calculation based on mean threshold: For the decomposition results of each deformation point, calculate the average amplitude of the first m-1=7 IMF components. Amplitude calculation uses the formula The average threshold of region R1 is obtained. =0.827mm;

[0134] Signal reconstruction: High-frequency terms Combined amplitude greater than IMF1~IMF3 represent short-cycle deformations with a period of 1-3 months;

[0135] low frequency items Combined amplitude less than or equal to IMF4 to IMF8 represent long-cycle deformations with a period of 6 to 18 months;

[0136] Trend Item Directly extract the residual term This represents the long-term subsidence trend from 2018 to 2022;

[0137] The reconstructed sequences of high-frequency terms, low-frequency terms, and trend terms are as follows: Figure 4 As shown.

[0138] V. Step S4: ICA blind source separation and vertical stratification information extraction:

[0139] S401: Construction of the ICA input matrix:

[0140] Set the high-frequency terms of the 30821 deformation points in region R1 respectively Low-frequency term set Trend Item Set Construct as input data matrix .

[0141] S402: Signal centering and whitening centering processing: for Centralize them separately, the formula is as follows ,in For matrix The mean vector in the time dimension;

[0142] Whitening process: Calculate the centered matrix covariance matrix ,right Eigenvalue decomposition yields the eigenvector matrix. and eigenvalue diagonal matrix Through formula To achieve whitening, the correlation between observed signals is eliminated, and the whitened signal... The covariance matrix is ​​the identity matrix.

[0143] S403: Independent component extraction uses the FastICA algorithm to extract independent components, and sets a non-quadratic function. The separation matrix is ​​solved by the following steps. :

[0144] Initialize the separation matrix ;

[0145] Iterative update : ,in , They are respectively The first and second derivatives;

[0146] Orthogonalization: for the updated Orthogonalization is performed to avoid extracting redundant components;

[0147] Convergence criterion: When When the iteration converges, the separation matrix is ​​output. ;

[0148] Independent component calculation: using the formula ,from Eight, eight, and four independent components were extracted from the samples, respectively.

[0149] S404: Vertical Stratified Quantitative Reorganization Based on Variance Contribution:

[0150] Using formula Calculate the variance contribution rate of each independent component, where Total number of independent components (high-frequency items) =8, low-frequency term =8, Trend Item =4), the calculation results are as follows:

[0151] High-frequency independent components: =42.3% (maximum) =28.7% (Second) =1.2% (minimum);

[0152] Low-frequency independent components: =3.1% (minimum) =25.6% (Second) =38.9% (maximum);

[0153] Independent components of the trend: =2.5% (minimum) =22.8% (Second) =58.7% (maximum).

[0154] Geological significance analysis: Combining the typical Quaternary stratigraphic structure of "shallow unconfined aquifer - intermediate slightly confined aquifer - confined aquifer - deep confined aquifer system", the high-frequency term with high variance contribution rate (VC1=42.3%) corresponds to the seasonal fluctuations of shallow groundwater and the influence of rainfall; the low-frequency term with high variance contribution rate (VC8=38.9%) corresponds to the influence of interannual water level changes in the intermediate aquifer; and the trend term with high variance contribution rate (VC4=58.7%) corresponds to the influence of plastic compression of deep cohesive soil.

[0155] Quantitative recombination: Obtaining vertical layered deformation information according to preset rules:

[0156] Shallow surface deformation = high frequency term (42.3%) + Low-frequency items (3.1%) + Trend Item (2.5%)

[0157] Middle-level surface deformation = high-frequency term (28.7%) + Low-frequency items (25.6%) + Trend Item (22.8%)

[0158] Deep surface deformation = high frequency term (1.2%) + Low-frequency items (38.9%) + Trend Item (58.7%)

[0159] The recombined shallow, medium and deep deformation time series are as follows: Figure 5 As shown.

[0160] VI. Step S5: Accuracy Verification

[0161] The vertical stratified deformation analysis results of zones R1 to R4 were verified using stratified marker monitoring data and leveling monitoring data.

[0162] Stratified marker verification: The measured deformation of shallow, middle and deep layers from three stratified marker monitoring stations was compared with the analytical results, and the average annual error from 2018 to 2022 was calculated. The results were 0.21 mm / year, 0.28 mm / year and 0.35 mm / year, respectively, all less than 0.5 mm / year.

[0163] Leveling verification: The cumulative measured settlement values ​​of 15 leveling monitoring points were compared with the cumulative values ​​of the analytical results of "shallow + middle + deep layers". The average absolute error was 0.42 mm and the relative error was 1.8%.

[0164] Validity determination: All verification errors were less than 0.5 mm / year, indicating that the method of the present invention is valid for the analysis results of the multi-source vertical contribution of surface deformation in the southern plain of Hebei.

[0165] In summary, using the analytical method of this invention, the contribution ratio of vertical stratification deformation in the southern plain area of ​​Hebei Province from 2018 to 2022 was obtained as follows:

[0166] Shallow deformation contributes 15.2% to 22.7% (mainly affected by rainfall and seasonal groundwater extraction);

[0167] The contribution of mid-layer deformation is 38.5% to 45.3% (mainly affected by interannual groundwater over-extraction).

[0168] The contribution of deep deformation is 32.0% to 46.3% (mainly affected by the secondary consolidation of deep cohesive soil).

[0169] This result provides a basis for "targeted treatment" of ground subsidence in the study area.

[0170] The above embodiments only illustrate preferred embodiments of the present invention, and their descriptions are relatively specific and detailed, but they should not be construed as limiting the scope of the present invention. It should be noted that those skilled in the art can make various modifications, improvements, and substitutions without departing from the concept of the present invention, and these all fall within the protection scope of the present invention.

Claims

1. A method for analyzing the vertical contribution of multi-source surface deformation in InSAR based on deformation signal spectral decomposition, characterized in that, The steps are as follows: S1. Acquire long-term synthetic aperture radar (SAR) image data covering the study area, with time points t0, t1, ..., tn; process the image data using SBAS-InSAR technology, extract the time-series deformation of the ground target points along the radar line of sight, combine it with the radar incident angle to convert it into vertical deformation, and obtain the annual average deformation rate V' of each target point and the cumulative deformation D0, D1, ..., Dn relative to the stable reference point at each observation time; S2, at the aforementioned average annual deformation rate Using the core feature as the basis, and integrating meteorological and hydrogeological auxiliary conditions, the K-Means clustering algorithm is used to spatially partition the set of deformation points in the region into K categories with similar deformation features and spatiotemporal patterns. S3. For each spatial partition, perform adaptive noise-complete ensemble empirical mode decomposition on the temporal deformation signals of all deformation points within the partition, decomposing the non-stationary deformation sequence of each point into... One intrinsic mode function and one residual term; using a multi-scale reconstruction method based on mean threshold, the pre-calculation... indivual The mean amplitude of the components is used as a threshold. IMFs with amplitudes greater than this threshold are merged into high-frequency terms, while those with amplitudes less than or equal to this threshold are merged into high-frequency terms. The terms are merged into low-frequency terms, and the residual terms are used as trend terms. S4. Apply independent component analysis to the three component datasets of high-frequency term, low-frequency term and trend term respectively to perform blind source separation and obtain multiple independent components; The variance contribution rate of each independent component is analyzed, and specific independent components from high-frequency, low-frequency and trend terms are recombined in a targeted manner according to preset rules to respectively retrieve information on shallow surface deformation, middle surface deformation and deep surface deformation. The preset recombination rule is as follows: Shallow surface deformation information = the component with the largest variance contribution rate in the high-frequency term + the component with the smallest variance contribution rate in the low-frequency term + the component with the smallest variance contribution rate in the trend term; Deep surface deformation information = the component with the smallest variance contribution rate in the high-frequency term + the component with the largest variance contribution rate in the low-frequency term + the component with the largest variance contribution rate in the trend term; Information on mid-level surface deformation = the component with the second-highest variance contribution rate in the high-frequency terms + the component with the second-highest variance contribution rate in the low-frequency terms + the component with the second-highest variance contribution rate in the trend terms.

2. The InSAR surface deformation multi-source vertical contribution analysis method based on deformation signal spectral decomposition according to claim 1, characterized in that: The SAR image data acquired in S1 are continuous observation images covering the same study area with a time span of not less than one year.

3. The InSAR surface deformation multi-source vertical contribution analysis method based on deformation signal spectral decomposition according to claim 1, characterized in that: The meteorological auxiliary conditions in S2 include the average annual precipitation and average annual evapotranspiration of the study area during the same period. The hydrogeological auxiliary conditions include the thickness of the first compressible layer, the thickness of the second compressible layer, the thickness of the third compressible layer, the variation of the first confined groundwater level, the variation of the second confined groundwater level, and the variation of the third confined groundwater level in the study area.

4. The InSAR surface deformation multi-source vertical contribution analysis method based on deformation signal spectral decomposition according to claim 1, characterized in that: The multi-scale reconstruction method based on mean threshold in S3 specifically includes the following steps: S301. Data Matrix Construction: Represent the InSAR deformation time series data within the partitions after K-Means clustering as a matrix. , in This indicates the number of deformation points within the partition. The length of the time series, matrix elements Indicates the first The deformation point at the first Deformation at each time point; S302, CEEMDAN Adaptive Decomposition: For the matrix Each row of timing signals in Performing CEEMDAN decomposition decomposes the original signal into a set of intrinsic mode functions and a sum of residual terms: in Indicates the first One eigenmode function Represents the residual term; S303. Multi-scale reconstruction based on amplitude mean: For the decomposition result of each deformation point, calculate its preceding value. indivual The average amplitude of the component As an adaptive threshold: The signal is then reconstructed based on this threshold, including the high-frequency term. From amplitude greater than of Composed of multiple components: low frequency items From amplitude less than or equal to of Composed of multiple components: Trend Item For residual terms :Right now .

5. The InSAR surface deformation multi-source vertical contribution analysis method based on deformation signal spectral decomposition according to claim 4, characterized in that, The blind source separation and vertical layering information extraction in S4 specifically include the following steps: S401, ICA Input Matrix Construction: Set the high-frequency terms of all deformation points within the partition. Low-frequency term set Trend Item Set As three independent input data matrices ,in This represents the total number of deformation points within the partition. S402, Signal Centering and Whitening: Centering is performed on each input matrix, i.e. ,in For matrix The mean vector over time; for the centered data Perform whitening treatment, that is ,in and They are covariance matrix The eigenvector matrix and the diagonal matrix of eigenvalues; S403. Independent Component Extraction: The FastICA algorithm is used to solve the separation matrix by maximizing the negative entropy. From the whitening signal Extracting independent components, namely ,in The separated independent component matrix; updated iteratively. Until convergence, to maximize the negative entropy objective function. : in Since it is not a quadratic function, take , It is a constant, and , It is a standard Gaussian random variable; S404. Vertical Stratified Quantitative Reorganization Based on Variance Contribution: Calculating the variance contribution rate of each independent component obtained after ICA separation of the three component datasets (high-frequency, low-frequency, and trend term). : in This is the variance calculation function. This represents the total number of independent components in the dataset. Based on the typical Quaternary stratigraphic structure, the geological significance of each component is analyzed, and quantitative recombination is performed according to preset rules to obtain shallow, middle and deep surface deformation information.

6. The InSAR surface deformation multi-source vertical contribution analysis method based on deformation signal spectral decomposition according to claim 5, characterized in that: The typical Quaternary stratigraphic structure, from top to bottom, includes: a shallow unconfined aquifer mainly composed of loose silt and fine sand; a middle layer of slightly confined to confined aquifer mainly composed of alternating layers of medium and fine sand and silty clay; and a deep confined aquifer mainly composed of thick layers of cohesive soil and dense gravel.

7. The InSAR surface deformation multi-source vertical contribution analysis method based on deformation signal spectral decomposition according to claim 1, characterized in that: When processing SAR image data using SBAS-InSAR technology in S1, it is achieved through Sarproz, GAMMA, or MintPy synthetic aperture radar interferometry software, and the processing includes image stitching, image registration, interferogram generation, phase unwrapping, and geocoding steps.

8. The InSAR surface deformation multi-source vertical contribution analysis method based on deformation signal spectral decomposition according to claim 5, characterized in that: The calculations for CEEMDAN decomposition, multi-scale reconstruction based on mean threshold, and independent component analysis are implemented using Matlab or Python programming tools, and the numerical precision used in the calculation process is retained to three decimal places.

9. The InSAR surface deformation multi-source vertical contribution analysis method based on deformation signal spectral decomposition according to claim 1, characterized in that... The method also includes step S5: using stratified monitoring data and leveling monitoring data to verify the accuracy of the shallow, middle and deep surface deformation information retrieved. When the verification error is less than 0.5 mm / year, the analytical results of the method are deemed valid.

Images