Method for Extracting Periodic and Trend Characteristics of InSAR Ground Subsidence Considering Spatial Heterogeneity

Abstract

The present invention discloses a method for extracting the periodic and trend characteristics of InSAR ground subsidence considering spatial heterogeneity. This method obtains long-time series synthetic aperture radar image data, uses the PS-InSAR technology to conduct high-precision monitoring of ground subsidence. On this basis, principal component analysis (PCA) is adopted to reduce the dimension of deformation information, and a combined algorithm integrating spatio-temporal spectral clustering and independent component analysis (ICA) of blind source signals is proposed to efficiently identify different types of subsidence characteristics, extract the periodic characteristics, trend characteristics, and random characteristics of different categories of ground subsidence, and can conduct a more comprehensive, detailed, and accurate analysis of ground subsidence, providing a more reliable and scientific basis for subsidence monitoring, risk assessment, and prevention and control decision-making, and having broad application prospects.

Description

Technical Field

[0001] The present invention relates to the fields of land subsidence monitoring, groundwater resource management, urban infrastructure monitoring, geological disaster warning, etc., and particularly to a method for extracting the periodic characteristics and trend characteristics of land subsidence by using InSAR technology. Technical Background

[0002] Land subsidence refers to the geological phenomenon in which the surface elevation gradually decreases due to the consolidation and compression of loose underground strata under the influence of natural or human factors. It is usually closely related to groundwater extraction, underground resource exploitation, geological structure, etc. With the acceleration of the urbanization process, land subsidence has caused serious impacts on infrastructure, ecological environment and social security. Traditional land subsidence monitoring methods, such as leveling and GPS monitoring, are limited by problems such as the sparsity of the monitoring network and complex data processing, and it is difficult to achieve large-scale and high-timeliness subsidence monitoring. In recent years, with the development of remote sensing technology and big data analysis technology, the land subsidence monitoring method based on InSAR has emerged and become one of the research hotspots in the field of land subsidence research. With the wide application of InSAR technology, the data monitoring method based on radar images can achieve large-scale subsidence monitoring. However, the settlement data obtained by InSAR technology usually has problems such as data redundancy and complex spatio-temporal characteristics, and different settlement characteristics are often intertwined, making it difficult to extract effective settlement information simply through traditional methods. Therefore, introducing data dimensionality reduction and signal decomposition technology into settlement data analysis has become the key to improving monitoring accuracy and revealing the settlement mechanism.

[0003] To solve these problems, the present invention innovatively combines technologies such as PCA dimensionality reduction, spatio-temporal spectral clustering, and blind source signal separation (ICA) to construct a new land subsidence monitoring and analysis method. First, the principal component analysis technology is used to reduce the dimensionality of multi-dimensional time-series deformation data; secondly, based on the data after PCA dimensionality reduction, combined with the spatio-temporal correlation of the underground aquifer boundary, the variation range of the groundwater level, and precipitation, the spatio-temporal spectral clustering algorithm is used to cluster and analyze the land subsidence data to identify the spatial distribution and time variation law of different subsidence categories, providing an important basis for subsequent signal analysis. Finally, for the land subsidence data of each category, the blind source signal separation technology is applied to separate the mixed deformation signal into periodic components, trend components and random components through the maximum independence criterion. In land subsidence monitoring, the periodic component mainly reflects the periodic subsidence caused by seasonal factors, climate change, etc.; the trend component reflects the long-term subsidence trend caused by factors such as groundwater level change and geological subsidence; the random component reflects the subsidence fluctuation caused by extreme weather time or local influence. By combining spatio-temporal spectral clustering and independent component analysis methods, different sources in the settlement signals of different regions can be more accurately identified and separated, providing a reliable basis for the prediction, monitoring and prevention of land subsidence. Summary of the Invention

[0004] The object of the present invention is to organically combine technologies such as PCA dimensionality reduction, spatio-temporal spectral clustering, and blind source signal separation (ICA), and construct a method for extracting the periodic and trend characteristics of InSAR ground subsidence considering spatial heterogeneity, innovatively solving problems such as data redundancy and complex spatio-temporal characteristics in ground subsidence. Through dimensionality reduction, pattern recognition, and signal decomposition of large-scale and high-dimensional subsidence data, the key spatio-temporal characteristics in the subsidence data can be effectively extracted, different subsidence characteristics and their variation laws can be revealed, thereby improving the monitoring accuracy, optimizing the risk assessment, and providing a scientific basis for prevention and control measures.

[0005] In order to achieve the above object, the present invention adopts the following technical solutions:

[0006] A method for extracting the periodic and trend characteristics of InSAR ground subsidence considering spatial heterogeneity, comprising the following steps;

[0007] S1. First, obtain long-term synthetic aperture radar (SAR) image data covering the study area; assume the image acquisition times are t0, t1, …, t n . Use the PS-InSAR technology to process these images, extract the deformation information of each stable permanent scatterer (PS) point on the ground surface in the radar line-of-sight direction, and convert it into vertical deformation. Finally, obtain the annual average subsidence rate V of each PS point and the cumulative subsidence amounts S0, S1, …, S at each observation time n ; These original comparison data are required for building the model;

[0008] S2. After obtaining the multi-dimensional time-series deformation data of the annual average subsidence rate V and the cumulative subsidence amounts S0, S1, …, S n , use the principal component analysis (PCA) method to perform dimensionality reduction processing on the data; by performing eigenvalue decomposition on the covariance matrix of the deformation data, select the principal components with the largest variance contribution, thereby effectively reducing the data dimension, retaining the main spatio-temporal information of the subsidence change, and determining the number of independent components for subsequent independent component analysis (ICA) of blind source signals based on this;

[0009] S3. Based on the data after PCA dimensionality reduction, apply the spatio-temporal spectral clustering algorithm to perform clustering analysis on the ground subsidence data; this algorithm identifies the spatial distribution and temporal variation law of the subsidence pattern by analyzing the spatio-temporal correlation of the underground aquifer boundary, the amplitude of the underground water level change, and precipitation, and divides the ground subsidence signal into multiple different categories, each category representing different subsidence characteristics;

[0010] S4. For the deformation data of each category, the independent component separation is carried out by using the blind source separation ICA technology; through the maximum independence criterion, the mixed signal is decomposed into periodic components, trend components and random components; then the periodic component data, trend component data and random component data are classified and sorted out.

[0011] Preferably, the spatio-temporal spectral clustering combined with the blind source separation ICA method mentioned in S3 and S4 specifically includes the following steps:

[0012] A1. Implementation of the spatio-temporal spectral clustering algorithm:

[0013] A1-1. Represent the InSAR deformation time series data after dimensionality reduction processing, the groundwater level amplitude, the thickness of the underground aquifer, and precipitation as a matrix X ∈ R 4×n×t :

[0014]

[0015] In the formula, n represents the number of deformation points; t is the total length of the deformation time series;

[0016] A1-2. Calculate the similarity between deformation points and construct the similarity matrix W:

[0017] According to the difference between every two deformation points Xi and Xj, calculate their similarity:

[0018]

[0019] Among them, δ is a parameter that controls the similarity measure and determines the distance threshold between points; a small δ value will result in a higher similarity, thereby improving the clustering accuracy;

[0020] Then, construct the degree matrix D, and its diagonal elements are the sum of the similarities in each row:

[0021]

[0022] Among them, is the sum of the similarities in the i-th row (1 ≤ i ≤ n);

[0023] A1-3. Calculate the Laplacian matrix L and the symmetric Laplacian matrix Lsym:

[0024]

[0025] A1-4. Calculation of eigenvectors:

[0026] By calculating the eigenvectors corresponding to the first k largest eigenvalues of L sym , construct the eigenvector space U ∈ R n×k ; for each eigenvector Uij Normalize by row to obtain matrix \(T\in R\) n×k , where each row represents an eigenvector;

[0027] A1-5, K-means clustering analysis:

[0028] Use the K-means algorithm to cluster each row vector of matrix \(T\) and discretize it into \(C\) classes; finally, the clustering label of each row vector corresponds to the data point in the original data to determine the clustering category of each deformation point;

[0029] A2. Implementation steps of blind source separation:

[0030] A2-1. Construct a signal matrix:

[0031] After completing spatio-temporal spectral clustering, each clustering cluster corresponds to a set of InSAR data. These data sets are formed into a matrix \(X\), which contains the observed signal data at multiple time points; the signal set after spatio-temporal spectral clustering can be used as \(X\), and its size is \(m\times t\), where \(m(1\leq m\leq n)\) is the number of InSAR observation points within the cluster, and \(t\) is the length of the time series;

[0032]

[0033] where \(x\) i,j is the observed value of the \(i\)-th deformation point at the \(j\)-th moment;

[0034] A2-2. Signal centering processing:

[0035] Before performing blind source separation on the signal, it is necessary to center the data; that is, subtract the mean of each signal point so that the mean of the data is zero:

[0036] \(X' = X - E(X)\) (6)

[0037] where \(E(X)\) is the mean of each column, and \(X'\) is the centered data matrix;

[0038] A2-3. Signal whitening processing:

[0039] Calculate the covariance matrix of the signal and perform eigenvalue decomposition on it to obtain the eigenvector matrix \(V\) and the eigenvalue matrix \(\Lambda\):

[0040]

[0041] \(C_V = V\Lambda\) (8)

[0042] \(Z = V\cdot X\) (9)

[0043] where \(Z\) is the whitened signal matrix;

[0044] A2-4. Definition of the objective function and maximization of negentropy:

[0045] The non-Gaussianity of the signal is a key index for ICA separation. To maximize the non-Gaussianity, the negentropy is selected as the objective function to quantify the non-Gaussianity of the signal, and its expression is:

[0046] J(b) = (E{G(b)} - E{G(C)}) 2 (10)

[0047] where G(·) is a non-linear function, and the logarithmic cosine function G(x) = log(cosh(x)) is often selected. By optimizing the negentropy, the maximum non-Gaussianity can be obtained, thus realizing signal separation.

[0048] A2-5. Estimation of the inverse transformation matrix:

[0049] To achieve signal separation, an inverse transformation matrix W needs to be found such that the signal matrix Y after the action of W is close to the source signal matrix D, that is:

[0050] Y = WX = WSD ≈ D (11)

[0051] where S is the mixing matrix and D is the source signal matrix.

[0052] A2-6. Iterative update of the inverse transformation matrix:

[0053] To obtain the optimal inverse transformation matrix W, the fixed-point iteration method is used to update W:

[0054] W(k + 1) = E[ZG′(W T (k)Z)] - E[ZG′′(W T (k)Z)]W(k) (12)

[0055] In the formula, G′ and G″ are the first derivative and the second derivative of G respectively; k is the iteration number. By continuously iteratively updating W until convergence.

[0056] A2-7. Normalization of the inverse transformation matrix:

[0057] To ensure the stability of the algorithm, it is necessary to normalize W after each update so that the norm of the inverse transformation matrix is 1:

[0058]

[0059] A2-8. Output of the signal separation result:

[0060] Through multiple iterations, the inverse transformation matrix W is finally obtained, and then the separated signal matrix Y is obtained; each column in this matrix is the independent source signal separated by ICA, which is close to the original source signal matrix D:

[0061] Y≈D(14)

[0062] Preferably, when obtaining the radar observation image data in S1, the radar image processing software used includes: Sarproz and MintPy.

[0063] Preferably, the calculations of the formulas (1) to (14) are implemented by calculation tools such as Matlab and Python.

[0064] Compared with the prior art, the innovation and beneficial effects of the present invention

[0065] The present invention combines principal component analysis (PCA), spatio-temporal spectral clustering and blind source signal separation (ICA) technologies for the first time, and conducts refined extraction of ground settlement characteristics in different regions, breaking through the limitation that traditional ground settlement technologies cannot fully consider spatial heterogeneity and the diversity of settlement characteristics. By integrating principal component analysis (PCA), spatio-temporal spectral clustering and blind source signal separation (ICA), not only the multi-dimensional complexity of settlement signals is considered, but also the periodic characteristics, trend characteristics and random characteristics of ground settlement can be effectively extracted, enabling a more comprehensive, detailed and accurate analysis of ground settlement characteristics, providing a more reliable and scientific basis for settlement monitoring, risk assessment and prevention and control decisions, and having a wide application prospect. The specific beneficial effects are as follows:

[0066] 1. The present invention uses principal component analysis (PCA) to perform dimensionality reduction processing on multi-dimensional ground settlement data, removing redundant information, reducing the computational complexity, and at the same time improving the effect of subsequent clustering.

[0067] 2. By combining spatio-temporal spectral clustering with the correlations of underground aquifer boundaries, groundwater levels, precipitation, etc., the settlement characteristics of different regions can be accurately identified, breaking through the problem that traditional methods ignore spatial heterogeneity.

[0068] 3. The blind source signal separation (ICA) technology decomposes the ground settlement signal into periodic, trend and random components, helps to reveal the complex causes of settlement, and provides effective monitoring of long-term changes and emergencies. BRIEF DESCRIPTION OF THE DRAWINGS

[0069] Figure 1 It is a schematic flow chart of the InSAR ground settlement periodicity and trend feature extraction method considering spatial heterogeneity proposed by the present invention;

[0070] Figure 2Schematic diagram of the PS-InSAR process for the method of extracting the periodic and trend characteristics of InSAR ground subsidence considering spatial heterogeneity proposed by the present invention;

[0071] Figure 3 Time series of InSAR rebound points and blind source signal separation results selected for the present invention;

[0072] Figure 4 Time series of InSAR subsidence points and blind source signal separation results selected for the present invention. Specific implementation manner

[0073] Next, the technical solutions in the embodiments of the present invention will be clearly and completely described in conjunction with the accompanying drawings in the embodiments of the present invention. Obviously, the described embodiments are only a part of the embodiments of the present invention, rather than all the embodiments.

[0074] Embodiment 1:

[0075] Please refer to Figure 1 , the method for extracting the periodic and trend characteristics of InSAR ground subsidence considering spatial heterogeneity includes the following steps;

[0076] S1. First, obtain long-term synthetic aperture radar (SAR) image data covering the study area. Assume the obtained image times are t0, t1,..., t n ; Use the PS-InSAR technology to process these images, extract the deformation information of each stable permanent scatterer (PS point) on the ground surface in the radar line-of-sight direction, and convert it into vertical deformation. Finally, obtain the annual average subsidence rate V of each PS point and the cumulative subsidence amounts S0, S1,..., S at each observation time n ;

[0077] S2. After obtaining the multi-dimensional time series deformation data of the annual average subsidence rate V and the cumulative subsidence amounts S0, S1,..., S n ; By performing eigenvalue decomposition on the covariance matrix of the deformation data, select the principal component with the largest variance contribution, thereby effectively reducing the data dimension and retaining the main spatio-temporal information of the subsidence change;

[0078] The specific operations include: calculating the covariance matrix of the data matrix, performing eigenvalue decomposition, and selecting the first few principal components. These principal components represent the main trends of the subsidence change and can reflect the main spatio-temporal characteristics of the ground subsidence. Based on the data after PCA dimensionality reduction, determine the number of independent components required for subsequent blind source signal separation (ICA).

[0079] S3. Based on the data after PCA dimensionality reduction, apply the spatio-temporal spectral clustering algorithm to cluster the land subsidence data. This algorithm identifies the spatial distribution and temporal variation patterns of land subsidence by analyzing the spatio-temporal correlations of factors such as the boundaries of underground aquifers, the amplitude of groundwater levels, and precipitation in different regions, and classifies the land subsidence signals into multiple different categories, each category representing different subsidence characteristics. The calculation process is as follows:

[0080] A1. Represent the InSAR deformation time series data after dimensionality reduction, the amplitude of groundwater level, the thickness of the underground aquifer, and precipitation as a matrix X ∈ R 4×n×t :

[0081]

[0082] where n represents the number of deformation points; t is the total length of the deformation time series;

[0083] A2. Calculate the similarity between deformation points and construct the similarity matrix W:

[0084] Calculate the similarity between each pair of deformation points Xi and Xj according to their differences:

[0085]

[0086] where δ is a parameter controlling the similarity measure, which determines the distance threshold between points. A small δ value will result in a higher similarity and thus improve the clustering accuracy.

[0087] Then, construct the degree matrix D, whose diagonal elements are the sum of similarities in each row:

[0088]

[0089] where is the sum of similarities in the i-th row (1 ≤ i ≤ n); ]>

[0090] A3. Calculate the Laplacian matrix L and the symmetric Laplacian matrix Lsym:

[0091]

[0092] A4. Calculation of eigenvectors:

[0093] By calculating the eigenvectors corresponding to the top k largest eigenvalues of L sym construct the eigenvector space U ∈ R n×k . Normalize each eigenvector U ij row by row to obtain the matrix T ∈ R n×k , where each row represents an eigenvector.

[0094] A5, K-means clustering analysis:

[0095] Use the K-means algorithm to cluster each row vector of matrix T and discretize it into C classes. Finally, the clustering labels of each row vector are corresponded to the data points in the original data to determine the clustering category of each deformation point.

[0096] S4. For the deformation data of each category, use the independent component analysis (ICA) technique to perform independent component separation. The ICA method is based on the maximum independence criterion and decomposes the mixed signal into independent periodic components, trend components, and random components. The calculation process is as follows:

[0097] A1. Construct a signal matrix:

[0098] After completing the spatio-temporal spectral clustering, each clustering cluster corresponds to a set of InSAR data. These data sets are formed into a matrix X, which contains the observed signal data at multiple time points. The signal set after spatio-temporal spectral clustering can be used as X, and its size is m×t, where m (1≤m≤n) is the number of InSAR observation points within the cluster, and t is the length of the time series.

[0099]

[0100] where, x i,j is the observed value of the i-th deformation point at the j-th moment.

[0101] A2. Signal centering processing: Before performing blind source separation, first perform centering processing on the signal, that is, subtract its mean from each signal point, so that the mean of the processed signal matrix is zero.

[0102] X′ = X - E(X) (6) [[ID=Z8]]

[0103] where, E(X) is the mean of each column, and X′ is the centered data matrix.

[0104] A3. Signal whitening processing: Calculate the covariance matrix of the signal and perform eigenvalue decomposition on it to obtain the eigenvector matrix V and the eigenvalue matrix Λ:

[0105]

[0106] CV = VΛ (8)

[0107] Z = V·X (9)

[0108] where, Z is the whitened signal matrix.

[0109] A4. Definition of the objective function and maximization of negentropy:

[0110] The non-Gaussianity of the signal is a key index for ICA separation. To maximize the non-Gaussianity, negentropy is selected as the objective function to quantify the non-Gaussianity of the signal, and its expression is:

[0111] J(b) = (E{G(b)} - E{G(C)}) 2 (10)

[0112] where G(·) is a non-linear function, and the logarithmic cosine function G(x) = log(cosh(x)) is often selected. By optimizing the negentropy, the maximum non-Gaussianity can be obtained, thus realizing signal separation.

[0113] A5. Estimation of the inverse transformation matrix:

[0114] To separate the signal, an inverse transformation matrix W needs to be found such that the signal matrix Y after the action of W is close to the source signal matrix D, that is:

[0115] Y = WX = WSD ≈ D (11)

[0116] where S is the mixing matrix and D is the source signal matrix.

[0117] A6. Iteratively update the inverse transformation matrix:

[0118] To obtain the optimal inverse transformation matrix W, the fixed-point iteration method is used to update W:

[0119] W(k + 1) = E[ZG′(W T (k)Z)] - E[ZG′′(W T (k)Z)]W(k) (12)

[0120] In the formula, G′ and G″ are the first derivative and the second derivative of G respectively; k is the number of iterations. By continuously iteratively updating W until convergence.

[0121] A7. Normalize the inverse transformation matrix:

[0122] To ensure the stability of the algorithm, it is necessary to normalize W after each update so that the norm of the inverse transformation matrix is 1:

[0123]

[0124] A8. Output the signal separation result:

[0125] Through multiple iterations, the inverse transformation matrix W is finally obtained, and then the separated signal matrix Y is obtained. Each column in this matrix is the independent source signal separated by ICA, which is close to the original source signal matrix D:

[0126] Y ≈ D (14).

[0127] Example 2:

[0128] Based on Example 1, the differences are as follows:

[0129] Refer to Figure 2 , in the present invention, the Sentinel-1A radar remote sensing images from 2018 to 2022 are first processed by Sarproz and MintPy software. The precision of the PS-InSAR monitoring results is verified by using leveling point data and layered marker monitoring data. The deformation information of each stable permanent scatterer (PS point) on the ground surface in the radar line-of-sight direction is extracted and converted into vertical deformation. The Hebei Plain area is used as the research area, and the change of the cumulative ground settlement deformation amount is obtained by using the monitoring method. And according to the position of each PS point, the groundwater level amplitude, aquifer boundary, precipitation, etc. are extracted.

[0130] The main steps are as follows:

[0131] Step 1: Using 136 scenes of Sentinel-1A ascending orbit SAR data covering the Hebei Plain from 2018 to 2022, the images are processed by Sarproz and MintPy software, including image registration, interferometric processing, phase unwrapping, etc., to obtain the deformation information in the line-of-sight direction, and it is converted into vertical deformation information through geometric correction. In order to verify the precision of the PS-InSAR results, error analysis is carried out using layered marker monitoring data, leveling point data, etc., so as to obtain high-precision ground settlement data.

[0132] Step 2: Based on the administrative regions and groundwater level data of the Hebei Plain, combined with the underground aquifer boundary information and precipitation raster data, through spatial analysis tools on the ArcGIS platform, the spatial characteristic data such as the groundwater level amplitude, underground aquifer thickness, and precipitation at the location of each PS point are extracted.

[0133] Step 3: Import the time series deformation data of the PS points into the matlab platform, and perform principal component analysis (PCA) for dimensionality reduction to extract the main deformation characteristics in the data, reduce the data dimension, and provide the simplified data input for the subsequent spatio-temporal spectral clustering analysis.

[0134] Step 4: Import the time series deformation data after PCA dimensionality reduction and the characteristics such as groundwater level amplitude, underground aquifer thickness, and precipitation into the Python platform, and perform clustering analysis through the spatio-temporal spectral clustering model. This model classifies each PS point into different categories according to the spatio-temporal correlation between the deformation information and the environmental characteristics.

[0135] Step 5: Apply the blind source signal separation technology to the deformation time series signals of each category, and decompose them into periodic components, trend components, and random components to further analyze the driving factors and evolution laws of various deformations.

[0136] The above are only the preferred specific embodiments of the present invention, but the protection scope of the present invention is not limited thereto. Any person skilled in the art within the technical scope disclosed by the present invention, according to the technical solution and inventive concept of the present invention, making equivalent substitutions or changes should be covered within the protection scope of the present invention.

Claims

1. Method for Extracting Periodic and Trend Characteristics of InSAR Ground Subsidence Considering Spatial Heterogeneity including the following steps; S1. First, obtain long-term synthetic aperture radar (SAR) image data covering the study area. Assume the image acquisition times are t0, t1, …, t n ; Use the PS-InSAR technique to process these images, extract the deformation information of each stable permanent scatterer (PS) point on the ground surface in the radar line-of-sight direction, and convert it into vertical deformation. Finally, obtain the annual average settlement rate V of each PS point and the cumulative settlement amounts S0, S1, …, S at each observation time n ; S2. After obtaining the annual average settlement rate V and the cumulative settlement amounts S0, S1, …, S of multi-dimensional time series deformation data at each observation time, the principal component analysis (PCA) method is used to perform dimensionality reduction on the data; by performing eigenvalue decomposition on the covariance matrix of the deformation data, the principal components with the largest variance contribution are selected, thereby effectively reducing the data dimension, retaining the main spatio-temporal information of the settlement change, and determining the number of independent components for subsequent independent component analysis (ICA) of blind source signals based on this; n After obtaining the multi-dimensional time series deformation data of n , the principal component analysis (PCA) method is used to reduce the dimension of the data; by performing eigenvalue decomposition on the covariance matrix of the deformation data, the principal components with the largest variance contribution are selected, so as to effectively reduce the data dimension, retain the main spatio-temporal information of the settlement change, and determine the number of independent components for subsequent independent component analysis (ICA) of blind source signals based on this; S3. Based on the data after PCA dimensionality reduction, apply the spatio-temporal spectral clustering algorithm to perform clustering analysis on the land subsidence data; this algorithm identifies the spatial distribution and temporal variation law of subsidence by analyzing the spatio-temporal correlation of the underground aquifer boundary, the amplitude change of the groundwater level, and precipitation, and divides the land subsidence signals into multiple different categories, with each category representing different subsidence characteristics; S4. For the deformation data of each category, adopt the blind source signal separation ICA technology to perform independent component separation; through the maximum independence criterion, decompose the mixed signal into periodic components, trend components, and random components; then classify and organize the periodic component data, trend component data, and random component data.

2. The method for extracting the periodic and trend characteristics of InSAR ground subsidence considering spatial heterogeneity according to claim 1, wherein The spatio-temporal spectral clustering combined with the blind source signal separation method in S3 and S4 specifically includes the following steps: A1. Implementation of the spatio-temporal spectral clustering algorithm: A1-1. Represent the InSAR deformation time series data after dimensionality reduction, the groundwater level amplitude, the thickness of the underground aquifer, and precipitation as a matrix \(X\in R\) 4×n×t :[[]]END]] In the formula, n represents the number of deformation points; t is the total length of the deformation time series; A1-2. Calculate the similarity between deformation points and construct the similarity matrix W: According to the difference between each two deformation points Xi and Xj, calculate their similarity: where δ is a parameter that controls the similarity measure and determines the distance threshold between points; a small δ value will result in a higher similarity, thereby improving the clustering accuracy; Then, construct the degree matrix D, and its diagonal elements are the sum of the similarities in each row: Among them, is the total similarity of the i-th row, where 1 ≤ i ≤ n; A1-3. Calculate the Laplacian matrix L and the symmetric Laplacian matrix Lsym: A1-4. Calculation of eigenvectors: By calculating the eigenvectors corresponding to the top k largest eigenvalues of L sym , a feature vector space U ∈ R n×k is constructed; for each eigenvector U ij , row normalization is performed to obtain a matrix T ∈ R n×k , where each row represents an eigenvector; A1-5. K-means clustering analysis: Use the K-means algorithm to cluster each row vector of the matrix T and discretize it into C categories; finally, the clustering label of each row vector corresponds to the data point in the original data to determine the clustering category of each deformation point; A2. Implementation steps of blind source signal separation: A2-1. Construct the signal matrix: After completing the spatio-temporal spectral clustering, each clustering cluster corresponds to a set of InSAR data; form these data sets into a matrix X, which contains the observed signal data at multiple time points; use the signal set after spatio-temporal spectral clustering as X, and its size is m×t, where m is the number of InSAR observation points within the cluster, 1≤m≤n; t is the length of the time series; where x i,j is the observed value of the i-th deformation point at the j-th moment; A2-2. Signal centering processing: Before performing blind source separation on the signal, it is necessary to perform centering processing on the data, that is, subtract the mean of each signal point so that the mean of the data is zero: X′ = X - E(X)(6) where E(X) is the mean of each column, and X′ is the centered data matrix; A2-3. Signal whitening processing: Calculate the covariance matrix of the signal and perform eigenvalue decomposition on it to obtain the eigenvector matrix V and the eigenvalue matrix Λ: CV = VΛ (8) Z = V·X (9) where Z is the whitened signal matrix; A2-4. Definition of the objective function and maximization of negentropy: The non-Gaussianity of the signal is the key index for ICA separation. To maximize the non-Gaussianity, the negentropy is selected as the objective function to quantify the non-Gaussianity of the signal, and its expression is: J(b) = (E{G(b)} - E{G(C)}) 2 (10) Among them, G(·) is a non-linear function, and the logarithmic cosine function G(x) = log(cosh(x)) is often selected to optimize the negative entropy to obtain the maximum non-Gaussianity, thereby realizing signal separation; A2-5. Estimation of the inverse transformation matrix: Find an inverse transformation matrix W to achieve signal separation, such that the signal matrix Y after the action of W is close to the source signal matrix D, that is: Y = WX = WSD ≈ D (11) where S is the mixing matrix and D is the source signal matrix; A2-6. Iteratively update the inverse transformation matrix: To obtain the optimal inverse transformation matrix W, the fixed-point iteration method is used to update W: W(k + 1) = E[ZG′(W T (k)Z)] - E[ZG″(W T (k)Z)]W(k) (12) In the formula, G′ and G″ are the first derivative and the second derivative of G respectively; k is the number of iterations; W is continuously iteratively updated until convergence; A2-7. Normalize the inverse transformation matrix: To ensure the stability of the algorithm, it is necessary to normalize W after each update so that the norm of the inverse transformation matrix is 1: A2-8. Output of the signal separation result: Through multiple iterations, the inverse transformation matrix W is finally obtained, and then the separated signal matrix Y is obtained; each column in this matrix is the independent source signal separated by ICA, which is close to the original source signal matrix D: Y ≈ D (14).

3. The method for extracting the periodic and trend characteristics of InSAR ground subsidence considering spatial heterogeneity according to claim 1, wherein: When acquiring radar observation image data, the radar image processing software used includes: Sarproz, GAMMA, and MintPy.

4. The method for extracting the periodic and trend characteristics of InSAR ground subsidence considering spatial heterogeneity according to claim 2, wherein: The calculations of formulas (1) to (14) are implemented by calculation tools such as Matlab and Python.

Images