InSAR (Interferometric Synthetic Aperture Radar) land subsidence periodicity and trend feature extraction method considering spatial heterogeneity

Abstract

The invention discloses an InSAR (Interferometric Synthetic Aperture Radar) land subsidence periodicity and trend feature extraction method considering spatial heterogeneity. According to the method, long-time-sequence synthetic aperture radar image data is acquired, high-precision monitoring is carried out on land subsidence by using a PS-InSAR technology, on the basis, dimension reduction processing is carried out on deformation information by using principal component analysis (PCA), a combinatorial algorithm integrating spatio-temporal spectral clustering and blind source signal separation ICA is provided, different types of subsidence characteristics are efficiently identified, and the accuracy of land subsidence monitoring is improved. According to the method, the periodic characteristics, the tendency characteristics and the randomness characteristics of different types of land subsidence are extracted, the land subsidence can be analyzed more comprehensively, meticulously and accurately, a more reliable and scientific basis is provided for subsidence monitoring, risk assessment and prevention and treatment decision, and the method has a wide application prospect.

Description

Technical Field

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

[0002] Land subsidence refers to the geological phenomenon that the surface elevation gradually decreases due to 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 urbanization, land subsidence has had a serious impact on infrastructure, ecological environment and social security. Traditional land subsidence monitoring methods, such as leveling and GPS monitoring, are limited by the sparsity of the monitoring network and the complexity of data processing, making it difficult to achieve large-scale and timely 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 has become one of the research hotspots in the field of land subsidence research. With the widespread application of InSAR technology, data monitoring methods based on radar images can achieve large-scale subsidence monitoring. However, the subsidence data obtained by InSAR technology usually have problems such as data redundancy and complex spatiotemporal characteristics, and different subsidence characteristics are often intertwined, making it difficult to extract effective subsidence information simply through traditional methods. Therefore, the introduction of data dimensionality reduction and signal decomposition technology in subsidence data analysis has become the key to improving monitoring accuracy and revealing subsidence mechanisms.

[0003] In order to solve these problems, the present invention innovatively combines PCA dimension reduction, spatiotemporal spectral clustering and blind source signal separation (ICA) to construct a new method for monitoring and analyzing land subsidence. First, the principal component analysis technology is used to reduce the dimension of multidimensional time series deformation data; secondly, based on the data after PCA dimension reduction, combined with the spatiotemporal correlation of underground aquifer boundaries, groundwater level fluctuations, and precipitation, the spatiotemporal spectral clustering algorithm is used to cluster the land subsidence data to identify the spatial distribution and temporal variation of different subsidence categories, providing an important basis for subsequent signal analysis. Finally, for each category of land subsidence data, 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 settlement caused by seasonal factors, climate change, etc.; the trend component reflects the long-term settlement trend caused by factors such as groundwater level changes and geological settlement; the random component reflects the settlement fluctuation caused by extreme weather time or local influence. By combining spatiotemporal spectral clustering and independent component analysis methods, it is possible to more accurately identify and separate different sources of subsidence signals in different regions, providing a reliable basis for the prediction, monitoring and prevention of ground subsidence. Summary of the invention

[0004] The purpose of this invention is to organically combine PCA dimensionality reduction, spatiotemporal spectral clustering and blind source signal separation (ICA) and other technologies to construct an InSAR ground subsidence periodicity and trend feature extraction method that takes into account spatial heterogeneity, and innovatively solve the problems of data redundancy and complex spatiotemporal features in ground subsidence. Through dimensionality reduction, pattern recognition and signal decomposition of large-scale, high-dimensional subsidence data, the key spatiotemporal features in the subsidence data can be effectively extracted, and different subsidence characteristics and their changing laws can be revealed, thereby improving monitoring accuracy, optimizing 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] The InSAR land subsidence periodicity and trend feature extraction method considering spatial heterogeneity includes the following steps:

[0007] S1. First, obtain long-time synthetic aperture radar SAR image data covering the study area; assume that the image acquisition time is t0, t1, …, t n The PS-InSAR technology is used to process these images, extract the deformation information of each stable permanent scatterer PS point on the surface in the direction of the radar line of sight, and convert it into vertical deformation, and finally obtain the annual average settlement rate V of each PS point and the cumulative settlement amount S0, S1, …, S at each observation time. n ; These original comparison data are needed to build the model;

[0008] S2, after obtaining the annual average sedimentation rate V and the cumulative sedimentation at each observation time S0, S1, ..., S n After the multi-dimensional time series deformation data is obtained, the principal component analysis (PCA) method is used to reduce the dimension of the data; the principal component with the largest variance contribution is screened out by performing eigenvalue decomposition on the covariance matrix of the deformation data, thereby effectively reducing the data dimension and retaining the main spatiotemporal information of the sedimentation change. Based on this, the number of independent components for subsequent blind source signal separation (ICA) is determined;

[0009] S3. Based on the data after PCA dimension reduction, the spatiotemporal spectral clustering algorithm is used to perform cluster analysis on the ground subsidence data. The algorithm identifies the spatial distribution and temporal variation of the subsidence pattern by analyzing the spatiotemporal correlation of the underground aquifer boundary, groundwater level fluctuation, and precipitation, and divides the ground subsidence signal into multiple different categories, each of which represents different subsidence characteristics.

[0010] S4. For each category of deformation data, blind source signal separation (ICA) technology is used to separate independent components. The mixed signal is decomposed into periodic components, trend components and random components through the maximum independence criterion. The periodic component data, trend component data and random component data are then classified and sorted.

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

[0012] A1. Implementation of spatiotemporal spectral clustering algorithm:

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

[0014]

[0015] Where 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 each two deformation points Xi and Xj, their similarity is calculated:

[0018]

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

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

[0021]

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

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

[0024]

[0025] A1-4. Calculation of eigenvectors:

[0026] By calculating L sym The eigenvectors corresponding to the first k largest eigenvalues ​​of construct the eigenvector space U∈R n×k ; For each eigenvector Uij Normalize by row and get the matrix T∈R n×k , where each row represents a feature vector;

[0027] A1-5. K-means cluster analysis:

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

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

[0030] A2-1. Constructing signal matrix:

[0031] After completing the spatiotemporal spectral clustering, each cluster corresponds to a set of InSAR data. These data sets form a matrix X, which contains the observed signal data at multiple time points; the signal set after spatiotemporal spectral clustering can be taken as X, whose size is m×t, where m (1≤m≤n) is the number of InSAR observation points in the cluster, and t is the length of the time series;

[0032]

[0033] Among them, x i,j is the observed value of the i-th deformation point at the j-th moment;

[0034] A2-2. Centralized signal processing:

[0035] Before performing blind source separation on the signal, the data needs to be centered; that is, the mean of each signal point is subtracted to make the mean of the data zero:

[0036] X′=XE(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 eigenvalue matrix Λ:

[0040]

[0041] CV=VΛ (8)

[0042] Z=V·X (9)

[0043] Among them, Z is the signal matrix after whitening;

[0044] A2-4. Objective function definition and negative entropy maximization:

[0045] The non-Gaussianity of the signal is a key indicator of ICA separation. In order to maximize the non-Gaussianity, negative entropy 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] Among them, G(·) is a nonlinear function, and the logarithmic cosine function G(x)=log(cosh(x)) is often selected. By optimizing the negative entropy, the maximum non-Gaussianity can be obtained, thereby achieving signal separation.

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

[0049] In order to achieve signal separation, it is necessary to find an inverse transformation matrix W so 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] Among them, S is the mixing matrix and D is the source signal matrix.

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

[0053] In order 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] Where G′ and G″ are the first-order derivative and second-order derivative of G, respectively; k is the number of iterations. W is updated through continuous iterations until convergence.

[0056] A2-7. Normalized inverse transformation matrix:

[0057] In order 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. Signal separation result output:

[0060] After 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 S1 acquires radar observation image data, radar image processing software used includes: Sarproz and MintPy.

[0063] Preferably, the calculations of formulas (1) to (14) are implemented using Matlab or Python computing tools.

[0064] Compared with the prior art, the invention has the following advantages:

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

[0066] 1. The present invention uses principal component analysis (PCA) to reduce the dimensionality of multidimensional land subsidence data, removes redundant information, reduces computational complexity, and improves the effect of subsequent clustering.

[0067] 2. Combining spatiotemporal spectral clustering with the association between underground aquifer boundaries, groundwater levels, precipitation, etc., it is possible to accurately identify the sedimentation characteristics of different regions, breaking through the problem of traditional methods ignoring spatial heterogeneity.

[0068] 3. Blind signal separation (ICA) technology decomposes ground subsidence signals into periodic, trend and random components, helping to reveal the complex causes of subsidence and provide effective monitoring of long-term changes and emergencies. BRIEF DESCRIPTION OF THE DRAWINGS

[0069] Figure 1 A schematic diagram of the process of the InSAR land subsidence periodicity and trend feature extraction method taking into account spatial heterogeneity proposed by the present invention;

[0070] Figure 2A schematic diagram of the PS-InSAR process of the InSAR land subsidence periodicity and trend feature extraction method proposed by the present invention taking into account spatial heterogeneity;

[0071] Figure 3 The InSAR rebound point time series and blind source signal separation result diagram selected by the present invention;

[0072] Figure 4 This is the result diagram of the InSAR settlement point time series and blind source signal separation selected by the present invention. DETAILED DESCRIPTION

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

[0074] Embodiment 1:

[0075] See also Figure 1 ,The InSAR ground subsidence periodicity and trend feature extraction method considering spatial heterogeneity includes the following steps;

[0076] S1. First, obtain long-time synthetic aperture radar (SAR) image data covering the study area. Assume that the image acquisition time is t0, t1, …, t n ; The PS-InSAR technology is used to process these images, extract the deformation information of each stable permanent scatterer (PS point) on the surface in the direction of the radar line of sight, and convert it into vertical deformation, and finally obtain the annual average settlement rate V of each PS point and the cumulative settlement amount S0, S1, …, S at each observation time. n ;

[0077] S2, after obtaining the annual average sedimentation rate V and the cumulative sedimentation at each observation time S0, S1, ..., S n After the multi-dimensional time series deformation data is collected, the principal component with the largest variance contribution is screened out by performing eigenvalue decomposition on the covariance matrix of the deformation data, thereby effectively reducing the data dimension and retaining the main spatiotemporal information of the settlement change.

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

[0079] S3. Based on the data after PCA dimension reduction, the spatiotemporal spectral clustering algorithm is used to cluster the ground subsidence data. The algorithm analyzes the spatiotemporal correlation of factors such as the boundaries of underground aquifers, groundwater level fluctuations, and precipitation in different regions, identifies the spatial distribution and temporal variation of subsidence, and divides the ground subsidence signals into multiple different categories, each of which represents different subsidence characteristics. The calculation process is as follows:

[0080] A1. Represent the InSAR deformation time series data after dimensionality reduction and the groundwater level variation, underground aquifer thickness, 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] According to the difference between each two deformation points Xi and Xj, their similarity is calculated:

[0085]

[0086] Among them, δ is a parameter that controls the similarity measurement and determines the distance threshold between points. Small δ values ​​will lead to higher similarity and thus improve clustering accuracy.

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

[0088]

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

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

[0091]

[0092] A4. Calculation of eigenvectors:

[0093] By calculating L sym The eigenvectors corresponding to the first k largest eigenvalues ​​of construct the eigenvector space U∈R n×k For each eigenvector U ij Normalize by row and get the matrix T∈R n×k , where each row represents a feature vector.

[0094] A5. K-means clustering analysis:

[0095] The K-means algorithm is used to cluster each row vector of the matrix T and discretize it into C categories. Finally, the cluster label of each row vector is matched with the data point in the original data to determine the cluster category of each deformation point.

[0096] S4. For each category of deformation data, blind source signal separation (ICA) technology is used to separate independent components. The ICA method is based on the maximum independence criterion and decomposes the mixed signal into independent periodic components, trend components and random components; its calculation process is as follows:

[0097] A1. Construct signal matrix:

[0098] After completing the spatiotemporal spectral clustering, each cluster corresponds to a set of InSAR data. These data sets form a matrix X, which contains the observed signal data at multiple time points. The signal set after spatiotemporal spectral clustering can be taken as X, whose size is m×t, where m (1≤m≤n) is the number of InSAR observation points in the cluster, and t is the length of the time series.

[0099]

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

[0101] A2. Signal centering: Before blind source separation, the signal is firstly centered, that is, the mean value of each signal point is subtracted from the signal point, so that the mean value of the processed signal matrix is ​​zero.

[0102] X′=XE(X) (6)

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

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

[0105]

[0106] CV=VΛ (8)

[0107] Z=V·X (9)

[0108] Among them, Z is the signal matrix after whitening.

[0109] A4. Objective function definition and negentropy maximization:

[0110] The non-Gaussianity of the signal is a key indicator of ICA separation. In order to maximize the non-Gaussianity, negative entropy 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] Among them, G(·) is a nonlinear function, and the logarithmic cosine function G(x)=log(cosh(x)) is often selected. By optimizing the negative entropy, the maximum non-Gaussianity can be obtained, thereby achieving signal separation.

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

[0114] In order to achieve signal separation, it is necessary to find an inverse transformation matrix W so 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] Among them, S is the mixing matrix and D is the source signal matrix.

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

[0118] In order 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] Where G′ and G″ are the first-order derivative and second-order derivative of G, respectively; k is the number of iterations. W is updated through continuous iterations until convergence.

[0121] A7. Normalized inverse transformation matrix:

[0122] In order 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. Signal separation result output:

[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] Embodiment 2:

[0128] Based on Example 1, but with the following differences:

[0129] Reference Figure 2 The present invention first processes the Sentinel-1A radar remote sensing images from 2018 to 2022 through Sarproz and MintPy software, verifies the accuracy of the PS-InSAR monitoring results by using leveling point data and layered mark monitoring data, extracts the deformation information of each stable permanent scatterer (PS point) on the surface in the direction of the radar line of sight, and converts it into vertical deformation. The Hebei Plain area is taken as the study area, and the changes in the accumulated deformation of ground subsidence are obtained by using the monitoring method. The position of each PS point is determined to extract the groundwater level variation, aquifer boundary, precipitation, etc.

[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 were processed by Sarproz and MintPy software, including image registration, interferometric processing, phase unwrapping, etc., to obtain deformation information in the line of sight, and then converted into vertical deformation information through geometric correction. In order to verify the accuracy of the PS-InSAR results, error analysis was performed using layered standard monitoring data and leveling point data to obtain high-precision ground subsidence data.

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

[0133] Step 3: Import the PS point time series deformation data into the MATLAB platform and perform principal component analysis (PCA) dimensionality reduction to extract the main deformation features in the data, reduce the data dimension, and provide simplified data input for subsequent spatiotemporal spectral clustering analysis.

[0134] Step 4: Import the time series deformation data after PCA dimension reduction and the characteristic inputs such as groundwater level fluctuation, underground aquifer thickness, and precipitation into the Python platform, and perform cluster analysis through the spatiotemporal spectral clustering model. This model classifies each PS point into different categories based on the spatiotemporal correlation between deformation information and environmental characteristics.

[0135] Step 5: Apply blind source signal separation technology to each category of deformation time series signal to decompose it 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 preferred specific implementation modes of the present invention, but the protection scope of the present invention is not limited thereto. Any technician familiar with the technical field can make equivalent replacements or changes according to the technical solutions and inventive concepts of the present invention within the technical scope disclosed by the present invention, which should be covered by the protection scope of the present invention.

Claims

1. InSAR land subsidence periodicity and trend feature extraction method considering spatial heterogeneity, The steps include: S1. First, obtain long-term synthetic aperture radar SAR image data covering the study area. Assume that the image acquisition time is t0, t1, …, t n ; The PS-InSAR technology is used to process these images, extract the deformation information of each stable permanent scatterer PS point on the surface in the direction of the radar line of sight, and convert it into vertical deformation, and finally obtain the annual average settlement rate V of each PS point and the cumulative settlement amount S0, S1, …, S at each observation time. n ; S2, after obtaining the annual average sedimentation rate V and the cumulative sedimentation at each observation time S0, S1, ..., S n After the multi-dimensional time series deformation data is obtained, the principal component analysis (PCA) method is used to reduce the dimension of the data; the principal component with the largest variance contribution is screened out by performing eigenvalue decomposition on the covariance matrix of the deformation data, thereby effectively reducing the data dimension and retaining the main spatiotemporal information of the sedimentation change. Based on this, the number of independent components for subsequent blind source signal separation (ICA) is determined; S3. Based on the data after PCA dimension reduction, the spatiotemporal spectral clustering algorithm is used to perform cluster analysis on the ground subsidence data. The algorithm identifies the spatial distribution and temporal variation of subsidence by analyzing the spatiotemporal correlation of the underground aquifer boundary, groundwater level fluctuation, and precipitation, and divides the ground subsidence signals into multiple different categories, each of which represents different subsidence characteristics. S4. For each category of deformation data, blind source signal separation (ICA) technology is used to separate independent components. The mixed signal is decomposed into periodic components, trend components and random components through the maximum independence criterion. The periodic component data, trend component data and random component data are then classified and sorted.

2. The InSAR land subsidence periodicity and trend feature extraction method considering spatial heterogeneity according to claim 1, characterized in that: The spatiotemporal spectral clustering combined with the blind source signal separation method in S3 and S4 specifically includes the following steps: A1. Implementation of spatiotemporal spectral clustering algorithm: A1-1. The InSAR deformation time series data after dimensionality reduction and the groundwater level fluctuation, underground aquifer thickness, and precipitation are expressed as a matrix X∈R 4×n×t : Where 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, their similarity is calculated: Among them, δ is a parameter that controls the similarity measurement and determines the distance threshold between points; a small δ value will lead to a higher similarity, thereby improving clustering accuracy; Then, construct the degree matrix D, whose diagonal elements are the sum of the similarities of each row: in, is the sum of similarities of the i-th row (1≤i≤n); A1-3. Calculate the Laplace matrix L and the symmetric Laplace matrix Lsym: A1-4. Calculation of eigenvectors: By calculating L sym The eigenvectors corresponding to the first k largest eigenvalues ​​of construct the eigenvector space U∈R n×k ; For each eigenvector U ij Normalize by row and get the matrix T∈R n×k , where each row represents a feature vector; A1-5. K-means cluster analysis: Use the K-means algorithm to cluster each row vector of the matrix T and discretize it into C categories; finally, the cluster label of each row vector corresponds to the data point in the original data to determine the cluster category of each deformation point; A2. Implementation steps of blind source signal separation: A2-1. Constructing signal matrix: After completing the spatiotemporal spectral clustering, each cluster corresponds to a set of InSAR data; these data sets form a matrix X, which contains the observed signal data at multiple time points. The signal set after spatiotemporal spectral clustering can be taken as X, whose size is m×t, where m (1≤m≤n) is the number of InSAR observation points in the cluster, and t is the length of the time series; Among them, x i,j is the observed value of the i-th deformation point at the j-th moment; A2-2. Centralized signal processing: Before performing blind source separation on the signal, the data needs to be centered. That is, the mean of each signal point is subtracted to make the mean of the data zero: X′=XE(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 eigenvalue matrix Λ: CV=VΛ (8) Z=V·X (9) Among them, Z is the signal matrix after whitening; A2-4. Objective function definition and negative entropy maximization: The non-Gaussianity of the signal is a key indicator of ICA separation. In order to maximize the non-Gaussianity, negative entropy 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 nonlinear function, and the logarithmic cosine function G(x) = log(cosh(x)) is often selected. By optimizing the negative entropy, the maximum non-Gaussianity can be obtained, thereby achieving signal separation; A2-5. Estimation of the inverse transformation matrix: In order to achieve signal separation, it is necessary to find an inverse transformation matrix W so 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) Among them, S is the mixing matrix, D is the source signal matrix; A2-6. Iteratively update the inverse transformation matrix: In order 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-order derivative and the second-order derivative of G respectively; k is the number of iterations; W is updated by continuous iteration until convergence; A2-7. Normalized inverse transformation matrix: In order 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. Signal separation result output: After 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 InSAR land subsidence periodicity and trend feature extraction method considering spatial heterogeneity according to claim 1, characterized in that : When acquiring radar observation image data, the radar image processing software used includes: Sarproz, GAMMA and MintPy.

4. The algorithm according to claim 2, characterized in that: The calculation of formulas (1) to (14) can be implemented using Matlab and Python computing tools.

Images