CN116363057B - Surface deformation extraction method integrating PCA and time sequence InSAR - Google Patents
Surface deformation extraction method integrating PCA and time sequence InSAR Download PDFInfo
- Publication number
- CN116363057B CN116363057B CN202310066005.3A CN202310066005A CN116363057B CN 116363057 B CN116363057 B CN 116363057B CN 202310066005 A CN202310066005 A CN 202310066005A CN 116363057 B CN116363057 B CN 116363057B
- Authority
- CN
- China
- Prior art keywords
- time
- matrix
- insar
- surface deformation
- time sequence
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000000605 extraction Methods 0.000 title claims abstract description 8
- 238000000034 method Methods 0.000 claims abstract description 35
- 238000000354 decomposition reaction Methods 0.000 claims abstract description 16
- 230000009467 reduction Effects 0.000 claims abstract description 12
- 238000005516 engineering process Methods 0.000 claims abstract description 8
- 239000000203 mixture Substances 0.000 claims abstract description 7
- 238000010586 diagram Methods 0.000 claims abstract description 5
- 239000011159 matrix material Substances 0.000 claims description 76
- 230000001932 seasonal effect Effects 0.000 claims description 11
- 239000006185 dispersion Substances 0.000 claims description 8
- 230000008569 process Effects 0.000 claims description 7
- 238000003384 imaging method Methods 0.000 claims description 4
- 238000000513 principal component analysis Methods 0.000 claims description 3
- 238000012876 topography Methods 0.000 claims description 3
- 230000008859 change Effects 0.000 claims description 2
- 230000002194 synthesizing effect Effects 0.000 claims description 2
- 230000003247 decreasing effect Effects 0.000 claims 1
- 238000012544 monitoring process Methods 0.000 abstract description 9
- 238000013101 initial test Methods 0.000 abstract description 2
- 238000007619 statistical method Methods 0.000 abstract 1
- 230000000875 corresponding effect Effects 0.000 description 8
- 238000004458 analytical method Methods 0.000 description 5
- 230000009286 beneficial effect Effects 0.000 description 3
- 238000011160 research Methods 0.000 description 3
- 239000000284 extract Substances 0.000 description 2
- 238000005259 measurement Methods 0.000 description 2
- 230000009471 action Effects 0.000 description 1
- 230000001174 ascending effect Effects 0.000 description 1
- 238000006243 chemical reaction Methods 0.000 description 1
- 230000002596 correlated effect Effects 0.000 description 1
- 230000008878 coupling Effects 0.000 description 1
- 238000010168 coupling process Methods 0.000 description 1
- 238000005859 coupling reaction Methods 0.000 description 1
- 238000005336 cracking Methods 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 238000001914 filtration Methods 0.000 description 1
- 230000004927 fusion Effects 0.000 description 1
- 238000005305 interferometry Methods 0.000 description 1
- 238000013507 mapping Methods 0.000 description 1
- 238000005065 mining Methods 0.000 description 1
- 238000012545 processing Methods 0.000 description 1
- 238000012216 screening Methods 0.000 description 1
- 230000002123 temporal effect Effects 0.000 description 1
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/0002—Inspection of images, e.g. flaw detection
- G06T7/0004—Industrial image inspection
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/60—Analysis of geometric attributes
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06V—IMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
- G06V10/00—Arrangements for image or video recognition or understanding
- G06V10/20—Image preprocessing
- G06V10/26—Segmentation of patterns in the image field; Cutting or merging of image elements to establish the pattern region, e.g. clustering-based techniques; Detection of occlusion
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06V—IMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
- G06V10/00—Arrangements for image or video recognition or understanding
- G06V10/70—Arrangements for image or video recognition or understanding using pattern recognition or machine learning
- G06V10/77—Processing image or video features in feature spaces; using data integration or data reduction, e.g. principal component analysis [PCA] or independent component analysis [ICA] or self-organising maps [SOM]; Blind source separation
- G06V10/774—Generating sets of training patterns; Bootstrap methods, e.g. bagging or boosting
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06V—IMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
- G06V10/00—Arrangements for image or video recognition or understanding
- G06V10/70—Arrangements for image or video recognition or understanding using pattern recognition or machine learning
- G06V10/82—Arrangements for image or video recognition or understanding using pattern recognition or machine learning using neural networks
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/20—Special algorithmic details
- G06T2207/20081—Training; Learning
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/20—Special algorithmic details
- G06T2207/20084—Artificial neural networks [ANN]
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/30—Subject of image; Context of image processing
- G06T2207/30108—Industrial image inspection
Landscapes
- Engineering & Computer Science (AREA)
- Theoretical Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Evolutionary Computation (AREA)
- Multimedia (AREA)
- Health & Medical Sciences (AREA)
- Artificial Intelligence (AREA)
- Computing Systems (AREA)
- Databases & Information Systems (AREA)
- General Health & Medical Sciences (AREA)
- Medical Informatics (AREA)
- Software Systems (AREA)
- Geometry (AREA)
- Quality & Reliability (AREA)
- Image Processing (AREA)
Abstract
The invention provides a surface deformation extraction method integrating PCA and time sequence InSAR, and relates to the field of surface deformation monitoring. Firstly, SAR image data are acquired, a differential interference diagram is generated based on a time sequence InSAR technology, high coherence points are extracted through statistical analysis of amplitude and phase information, then phase unwrapping is carried out, and a time sequence InSAR earth surface deformation initial result is obtained. The initial test result is centered, the centered result is decomposed through a singular value decomposition algorithm, data dimension reduction is implemented by integrating data quantity and space-time rules of each component, the signal composition of each main component is interpreted through analyzing the space distribution rule and the time evolution rule of each main component, a surface deformation model is constructed, and error signals are fitted and removed through a polynomial model, so that a final time sequence surface deformation result is obtained. The method only utilizes SAR data, does not introduce external errors, and provides scientific theoretical basis for high-precision surface deformation monitoring.
Description
Technical Field
The invention relates to the field of earth surface deformation monitoring, in particular to an earth surface deformation extraction method integrating PCA and time sequence InSAR.
Background
The deformation of the earth surface is a series of geological disasters such as inclined cracking of a building, deformation of a traffic infrastructure, subsidence of the earth surface and the like which are caused by the phenomena of lifting, tilting, dislocation and the like of the earth surface due to uneven stress under the action of natural or artificial activities. Therefore, it is necessary to monitor the surface deformation of the region of the geological disaster periodically and accurately. Traditional surveying and mapping means, such as triangle elevation measurement, leveling measurement, global satellite navigation system (Global Navigation Satellite System, GNSS) technology and the like, are used for acquiring the earth surface deformation information, and have the advantages of high precision, long working period, low spatial resolution and difficulty in accurately describing the earth surface deformation space-time characteristics of the region. The synthetic aperture radar interferometry (Interferometric Synthetic Aperture Radar, inSAR) technology can acquire high-precision and high-spatial-resolution earth surface deformation information all the time and all the weather, and overcomes the limitation of discrete point monitoring. The time sequence InSAR technology developed in recent years can acquire long-time sequence earth surface deformation information on the basis of ensuring spatial resolution, and is beneficial to describing and revealing earth surface deformation time-space evolution rules.
The earth surface deformation information processed by the time sequence InSAR is often accompanied with noise phases such as atmospheric phase delay, orbit error, residual DEM error and the like. In the existing earth surface deformation time sequence InSAR monitoring application research, a space-time filtering method is mostly adopted to remove the atmospheric delay error and noise signals, and the method has the problems that new noise is introduced or deformation information is misjudged to be noise to be removed and the like. The PCA method is a common data reduction and analysis method, a group of variable data possibly with correlation is converted into a group of linear uncorrelated variables through positive-negative conversion, the PCA method is mainly used for extracting main characteristic components of the data, and the PCA method is rarely applied to the InSAR field at present. Therefore, the method introduces a principal component analysis (Principal Component Analysis, PCA) method into the time sequence InSAR signal analysis mining, extracts a time sequence earth surface deformation model through the space-time analysis of the principal component, and further constructs a high-precision time sequence InSAR earth surface deformation information extraction strategy. The method utilizes SAR self data, does not introduce external errors, and provides scientific theoretical basis for high-precision ground surface deformation monitoring.
Disclosure of Invention
Aiming at the problems existing in the prior art, the invention provides the earth surface deformation extraction method integrating PCA and the time sequence InSAR, which can effectively overcome the defect that deformation information is misjudged as noise and the like when the time sequence InSAR method extracts earth surface deformation signals.
In order to achieve the technical purpose, the invention discloses a surface deformation extraction method for fusing PCA and time sequence InSAR, which comprises the following steps:
generating a differential interference diagram from the acquired SAR image data by using a time sequence InSAR technology;
extracting high coherence points in the differential interferogram based on amplitude dispersion and phase stability principles, and then performing phase unwrapping to generate an initial InSAR time sequence result;
carrying out centering treatment on the initial InSAR time sequence result, and decomposing the centered result through a Singular Value Decomposition (SVD) algorithm;
data dimension reduction is implemented by integrating the data quantity and the space-time law of each component, and the signal composition of each main component is interpreted by analyzing the space distribution law and the time evolution law of each main component; each principal component space comprises a space matrix R, a weight matrix S and a time matrix T, wherein the space matrix represents the distribution rule of the principal components in space, the time matrix represents the distribution rule of the components in time, and the weight matrix represents the importance of each component;
and constructing a surface deformation model by utilizing the interpretation result of each principal component signal, constructing a quadratic polynomial model, fitting error signals through error components in principal component components, and removing to obtain a final time sequence surface deformation result.
Further, the process of generating the differential interferogram is as follows:
q SAR images of a certain area on different dates are obtained, and each SAR image is marked with imaging time attribute, backscattering information and satellite imaging orbit information;
performing geographic coding on the q images by utilizing SRTM DEM data; combining time base line and space vertical base line factors, selecting one piece of SAR image data as an optimal public main image based on a space-time base line threshold value and a Doppler mass center frequency threshold value, and then taking other q-1 pieces of SAR images as auxiliary images; q-1 auxiliary images after registration are obtained through registering each auxiliary image with the main image;
inputting the q SAR images into GAMMA software, setting a short baseline space-time baseline threshold value in the GAMMA software to form a z pair short baseline space-time interference pair, and performing differential interference treatment on the generated short baseline space-time interference pair to generate a z short baseline differential interference diagram; wherein, the short baseline time-space baseline threshold time baseline is set in the GAMMA software and is less than or equal to 60 days, and the space baseline is less than or equal to 400 meters.
Further, the process of generating the initial InSAR time series result is as follows:
acquiring high-coherence points in q SAR images on different dates based on amplitude dispersion and phase stability principles by using StaMPS software, and arranging all the high-coherence points according to a time sequence to form a data set, wherein the data set comprises p high-coherence points;
carrying out phase unwrapping on the z differential interferograms by adopting a 3D unwrapping method, and estimating and removing phase errors caused by DEM and space-time base lines to obtain an initial InSAR time sequence result W p×q 。
Further, PCA is analyzed by adopting fusion principal component to process initial InSAR time sequence result W p×q Generating a plurality of main components, wherein the specific steps are as follows:
3) Centralizing an initial InSAR time sequence result:
initial InSAR time series result W p×q ＝{w 1 ，w 2 ，w 3 ，......，w q -w is i High-coherence point deformation information, w, of the ith image relative to the first image i Mean value of mu i Wherein (i=1, 2,) is @, q, w @. i ＝w i -μ i I.e. each bit of the features minus the respective average value, the centred matrix is W' p×q 。
4) Centering the matrix W 'through SVD algorithm' p×q Preliminary decomposition:
wherein R is p×p Andis a standard orthogonal matrix, sigma p×q Is a diagonal matrix, where each term is a singular value of the matrix;
and carrying out data dimension reduction on the preliminary decomposition result, trying different dimension reduction dimensions v to generate a main component decomposition scheme, comparing the results of different schemes, synthesizing data quantity and the space-time rule recognition characteristics of each component, and selecting the optimal scheme.
Further, the specific steps of selecting the optimal scheme by integrating the data quantity and the recognition characteristics of the space-time law of each component are as follows:
1a) For W 'and W' T Matrix multiplication is carried out to obtain a p multiplied by p dimensional square matrix W' T Wherein W 'is the matrix W' p×q For the initial InSAR time series result W p×q Results after the centering, W' T For the transpose of W', the eigenvalues and eigenvectors of the square matrix are solved using:
(W′W′ T )r j ＝λ j r j
where (j=1, 2,) p, p is a high coherence point number, r j As a feature vector lambda j The characteristic value is the characteristic value corresponding to the characteristic vector;
at this time, sigma is singular value, r is left singular vector, singular value sigma is arranged from big to small in matrix, sigma is reduced quickly, so that the matrix is approximately described by using the first v big singular values, then
1b) For W' T And W 'to obtain a q multiplied by q matrix W' T W', solving eigenvalues and eigenvectors of the square matrix:
(W′ T W′)t i ＝η i t i
wherein (i=1, 2,) q, q being the number of image frames, t i As a feature vector, eta i The characteristic value is the characteristic value corresponding to the characteristic vector;
at this time, μ is a singular value, t is a right singular vector, and the singular values μ are arranged from large to small in the matrix and decrease rapidly, so that the matrix is approximately described by using the first v large singular values, then
1c) W' T And W' T The singular value of W' is square root to form S v×v ；
The SVD algorithm data dimension-reducing matrix W 'is obtained through the steps' p×q The following results:
wherein v is the number of principal components in the principal component space, W' p×q R is the initial result of the time sequence InSAR after centering p×v And T v×q Respectively a space matrix and a time matrix, S v×v For a diagonal matrix, the value of which represents the importance of each component,is T v×q Is a transpose of (a).
Further, interpreting the signal composition of each principal component by analyzing the spatial distribution rule and the time evolution rule of each principal component, wherein the spatial distribution rule is analyzed by an RxS image, and the time evolution rule is analyzed by a T function; r is used for representing deformation distribution, S represents importance of principal component components, T represents time law of deformation, and the specific steps are as follows:
if the complex principal components in all principal components show time correlation, and the space pattern shown by R multiplied by S is consistent with the prior knowledge of the surface deformation space distribution, judging that the principal components showing time correlation are linear surface deformation related components, and if the deformation components show seasonal correlation, estimating as seasonal deformation;
if the component values of the error items T in all the main components fluctuate up and down at zero values, namely the components have no obvious time correlation; wherein if RxS has spatial correlation, the component can be presumed to be a residual orbit or long wavelength error, and if RxS exhibits seasonal variation, the error term can be presumed to be caused by seasonal atmospheric disturbance;
if R S does not exhibit significant deformation characteristics, the error term may be presumed to be a noise-related component.
Further, the error signal is fitted and removed through the error component in the main component, and the strategy for extracting the surface deformation information is as follows:
2a) Extracting the main components related to the surface deformation, and obtaining an optimal time function Q 'through function fitting time evolution rules' t Spatial component R of the corresponding component t And an importance index S t Reconstructing the time sequence earth surface deformation model X based on the following formula d And X is taken as d Pre-removing in the InSAR signal;
X d ＝R t ·S t ·Q t ′ T
2b) Extracting error principal components, reconstructing an error term model by using the function fitting method in the step 2 a), and removing the error term model from the InSAR signal;
2c) Using a quadratic polynomial to fit the principal component signals to remove the residual signals after the models of 2 a) and 2 b), estimating and removing residual atmosphere, orbit, topography error and noise phase which may exist in each interference pair on the time sequence;
wherein m is the high coherence point number in the interferogram, n is the timing number, delta m，n For InSAR initial interferometric phase, x m And y m Representing pixel coordinates, z m The elevation of the pixel points is a, b, … … and h are fitting parameters;
2d) Time sequence earth surface deformation model X to be removed in advance d And adding back the residual time sequence InSAR signal containing nonlinear deformation to obtain a final time sequence earth surface deformation result.
The beneficial effects are that:
the traditional time sequence InSAR error removal method can lead to signal misjudgment caused by coupling of deformation phase and error phase when the earth surface deformation is equivalent to the error signal magnitude, even the error signal magnitude is larger. The initial test result is centered, the centered result is decomposed through singular value decomposition (Singular Value Decomposition, SVD) algorithm, data dimension reduction is implemented by integrating data quantity and space-time law of each component, and the signal composition of each main component is interpreted through analyzing the space distribution law and time evolution law of each main component. And constructing a surface deformation model, fitting an error signal by using a polynomial model, and removing the error signal to obtain a final time sequence surface deformation result. The method only utilizes SAR data, does not introduce external errors, and provides scientific theoretical basis for high-precision surface deformation monitoring.
Drawings
Fig. 1 is a flowchart of a high-precision earth surface deformation monitoring method integrating PCA and time sequence InSAR technology according to an embodiment of the present invention.
Fig. 2 is a main component decomposition result of an initial timing InSAR signal according to an embodiment of the present invention.
FIG. 3 is a graph of the final time-series InSAR earth surface deformation rate according to one embodiment of the present invention.
Detailed Description
In order to make the technical problems, technical schemes and beneficial effects solved by the invention more clear, the invention is further described in detail below with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are for purposes of illustration only and are not intended to limit the scope of the invention.
1. In the embodiment, xuzhou area is selected as a monitoring object; the research area covers a Xuzhou urban area range (33.96-34.46 DEG, 116.97-117.56 DEG E), 66-view C-band ascending Sentinel 1A star (Sentinel-1A) SLC SAR image data of the research area is selected, and the image acquisition time is 2018, 8, 19, 3, 31, and 2022. Acquiring SRTM DEM data with 30m resolution ratio, simulating and removing a terrain phase, and geocoding an interferogram;
embodiments select to use GAMMA and StaMPS (Stanford Method for Persistant Scatterers Algorithm) for SBAS-InSAR technology (Smal l Basel ine Subset InSAR, SBAS-InSAR) process analysis, the method flow chart is shown in FIG. 1. Selecting images acquired on 29 days of 1 month in 2020 as public main images, and registering;
use 5:1, performing multi-view processing on each interference pair as multi-view comparison, setting a time base line threshold value to be less than or equal to 60 days, and a space base line threshold value to be less than or equal to 300m, so as to generate 150 pairs of differential interference patterns;
2. and selecting a high-coherence point for each image in the interference graph based on the amplitude dispersion and phase stability principles, wherein the method comprises the following specific steps of:
1) Setting an amplitude dispersion index to select a high-coherence point, namely removing pixels which do not meet the standard based on the amplitude dispersion, and primarily acquiring the high-coherence point;
2) Estimating high-coherence points after amplitude dispersion screening through phase stability, optimizing phases in iteration, and finally obtaining 216739 high-coherence points;
phase unwrapping is carried out on each pair of interferograms through 3D unwrapping to obtain 150 pairs of unwrapped interferograms, then phase errors caused by DEM and space-time baselines are estimated and removed, and an initial InSAR time sequence result W is obtained 216739×66 Wherein 216739 is a high coherence point number, 66 is a time-series image amplitude number;
3. analysis of initial InSAR time series results W by PCA method 216739×66 Generating a plurality of main components, wherein the specific steps are as follows:
1) Centralizing an initial InSAR time sequence result:
initial InSAR time series result W 216739×66 ＝{w i ，w 2 ，w 3 ，......，w 66 I=1, 2, &..66), wherein w i High-coherence point deformation information, w, of the ith image relative to the first image i Mean value of mu i ，w′ i ＝w i -μ i I.e. each bit of the features minus the respective average value, the centred matrix is W' 216739×66 。
2) W 'is processed by SVD algorithm' 216739×66 The preliminary decomposition is as follows:
wherein R is 216739×216739 Andis a standard orthogonal matrix, sigma 216739×66 Is a diagonal matrix, where each term is a singular value of the matrix;
4. and (3) carrying out data dimension reduction on the preliminary decomposition result, trying different dimension reduction dimensions (v) to generate a main component decomposition scheme, wherein when the number of main component components is 6, each main component space-time rule has obvious identification characteristics, and the dimension reduction comprises the following specific steps:
1a) For W 'and W' T Matrix multiplication is carried out to obtain a 216739 multiplied by 216739 dimensional square matrix W' T Wherein W' T Is the transpose of W' and,
solving eigenvalues and eigenvectors of the square matrix:
(W′W′ T )r j ＝λ j r j
wherein (j=1, 2,) is 216739, r j As a feature vector lambda j The characteristic value is the characteristic value corresponding to the characteristic vector;
at this time, sigma is singular value, r is left singular vector, singular value sigma is arranged from big to small in matrix, sigma is reduced quickly, so we can approximate description matrix with the first v big singular values, then
1b) For W' T And W 'to obtain 66X 66 dimensional square matrix W' T W', solving eigenvalues and eigenvectors of the square matrix:
(W′ T W′)t i ＝η i t i
wherein the method comprises the steps ofWherein (i=1, 2,) 66, t i As a feature vector, eta i The characteristic value is the characteristic value corresponding to the characteristic vector;
at this time, mu is a singular value, t is a right singular vector, and the singular value mu is arranged from large to small in the matrix and is reduced quickly, so that we can describe the matrix approximately by using the first v large singular values, then
1c) W' T And W' T The singular value of W' is square root to form S 6×6 ；
The SVD algorithm data dimension-reducing matrix W 'is obtained through the steps' p×q The following results:
5. each principal component space consists of a space matrix (R), a weight matrix (S) and a time matrix (T), wherein the space matrix represents the distribution rule of the principal components in space, the time matrix represents the distribution rule of the components in time, and the weight matrix represents the importance of each component. Wherein W 'is' 216739×66 R is the initial result of the time sequence InSAR after centering 216739×6 And T 6×66 Respectively a space matrix and a time matrix, S 6×6 For a diagonal matrix, the value of which represents the importance of each component,is T 6×66 The main component decomposition result of the initial time sequence InSAR signal is obtained and is shown in figure 2; fig. 2 (a) - (f) are spatial distribution results corresponding to the first to sixth principal components, respectively, and fig. 2T 1-T6 are time-function variation trends corresponding to the first to sixth principal components, respectively.
6. The signal composition of each principal component is interpreted through the spatial distribution rule and the time evolution characteristic of each principal component in fig. 2, and the specific steps are as follows:
the first Principal Component (PCI) and the second principal component (PC 2) are in certain time correlation (T1 and T2 are in linear change trend along with time), and the spatial pattern of the first principal component and the second principal component is matched with the prior knowledge of the surface deformation spatial distribution, so that the PCI and the PC2 can be judged to be linear surface deformation related components; the sixth principal component (PC 6) exhibits a certain seasonal variation with time, and deformation characteristics exist in the spatial distribution map, which can be determined to be related to seasonal nonlinear deformation;
the third principal component (PC 3) does not show a significant temporal correlation, the T3 value of the entire study period fluctuates up and down at zero values, and its r3×s3 value has a certain spatial correlation, but no significant correlation with elevation, which can be presumed to be mainly a residual orbit or long wavelength error. A fourth principal component (PC 4) in which the r4×s4 value of PC4 is correlated with Gao Chenglve, but does not exhibit a significant deformation characteristic in space, and thus may be caused by seasonal atmospheric disturbances; the time function of the fifth principal component (PC 5) randomly fluctuates around a zero value, the spatial distribution does not show obvious deformation characteristics, and the noise-related component is judged;
7. based on the interpretation mode, PCA results are obtained, and the following strategies are adopted to extract the surface deformation information:
2a) Extracting principal components PCI and PC2 related to the surface deformation, and obtaining an optimal time function Q 'through function fitting time evolution rules' t The spatial component Rt and the importance index St of the corresponding component reconstruct the time-series surface deformation model X based on the following formula d And X is taken as d Pre-removing in the InSAR signal;
X d ＝R t ·S t ·Q t ′ T
2b) Extracting error components PC3 and PC5, reconstructing an error term model by using the function fitting method in the step 2 a), and removing the error term model from the InSAR signal;
2c) Using a quadratic polynomial to fit the principal component signals to remove the residual signals after the models of 2 a) and 2 b), estimating and removing residual atmosphere, orbit, topography errors and noise phases possibly existing in the residual signals of each interference pair on a time sequence;
wherein m is the high coherence point number in the interferogram, n is the timing number, delta m，n For InSAR initial interferometric phase, x m And y m Representing pixel coordinates, z m The elevation of the pixel points is a, b, … … and h are fitting parameters;
2d) X pre-removed d And adding the model back to the residual time sequence InSAR signal containing the nonlinear deformation to obtain a final time sequence earth surface deformation result, as shown in figure 3.
Claims (5)
1. A surface deformation extraction method integrating PCA and time sequence InSAR is characterized by comprising the following steps:
generating a differential interference diagram from the acquired SAR image data by using a time sequence InSAR technology;
extracting high coherence points in the differential interferogram based on amplitude dispersion and phase stability principles, and then performing phase unwrapping to generate an initial InSAR time sequence result;
carrying out centering treatment on the initial InSAR time sequence result, and decomposing the centered result through a Singular Value Decomposition (SVD) algorithm;
data dimension reduction is implemented by integrating the data quantity and the space-time law of each component, and the signal composition of each main component is interpreted by analyzing the space distribution law and the time evolution law of each main component; each principal component space comprises a space matrix R, a weight matrix S and a time matrix T, wherein the space matrix represents the distribution rule of the principal components in space, the time matrix represents the distribution rule of the components in time, and the weight matrix represents the importance of each component;
constructing a surface deformation model by using the interpretation result of each principal component signal, constructing a quadratic polynomial model, fitting error signals by error components in principal component components, and removing to obtain a final time sequence surface deformation result;
interpreting the signal composition of each principal component by analyzing the spatial distribution rule and the time evolution rule of each principal component, wherein the spatial distribution rule is analyzed by an R multiplied by S image, and the time evolution rule is analyzed by a T function; the spatial distribution condition of the main component is represented by a spatial matrix R, the importance of the main component is represented by a weight matrix S, and the time change rule of the main component is represented by a time matrix T, and the specific steps are as follows:
if the complex principal components in all principal components show time correlation, and the spatial pattern shown by R multiplied by S is consistent with the prior knowledge of the surface deformation spatial distribution, judging that the principal components showing time correlation are linear surface deformation related components, and if the linear surface deformation related components show seasonal correlation, estimating as seasonal deformation;
if the error term T component values in all the main components fluctuate up and down at zero values, namely the linear surface deformation related components have no obvious time correlation; wherein if RxS has spatial correlation, the principal component can be presumed to be a residual orbit or long wavelength error, and if RxS exhibits seasonal variation, the error term can be presumed to be caused by seasonal atmospheric disturbance;
if R x S does not exhibit significant deformation characteristics, the error term may be presumed to be a noise-related component;
fitting and removing error signals through error components in the main component components, wherein the strategy for extracting the surface deformation information is as follows:
2a) Extracting the main components related to the surface deformation, and obtaining an optimal time function Q 'through function fitting time evolution rules' t Spatial component R of the corresponding component t And an importance index S t Reconstructing the time sequence earth surface deformation model X based on the following formula d And X is taken as d Pre-removing in the InSAR signal;
X d ＝R t ·S t ·Q t ' T
2b) Extracting error principal components, reconstructing an error term model by using the function fitting method in the step 2 a), and removing the error term model from the InSAR signal;
2c) Using a quadratic polynomial to fit the principal component signals to remove the residual signals after the models of 2 a) and 2 b), estimating and removing residual atmosphere, orbit, topography error and noise phase existing in each interference pair on a time sequence;
wherein m is the high coherence point number in the interferogram, n is the timing number, delta m,n For InSAR initial interferometric phase, x m And y m Representing pixel coordinates, z m The elevation of the pixel points is a, b, c, d, e, f, g and h which are fitting parameters;
2d) Time sequence earth surface deformation model X to be removed in advance d And adding back the residual time sequence InSAR signal containing nonlinear deformation to obtain a final time sequence earth surface deformation result.
2. The method for extracting the surface deformation by fusing PCA and time sequence InSAR as set forth in claim 1, wherein the process of generating the differential interferogram is as follows:
q SAR images of a certain area on different dates are obtained, and each SAR image is marked with imaging time attribute, backscattering information and satellite imaging orbit information;
performing geographic coding on the q images by utilizing SRTM DEM data; combining time base line and space vertical base line factors, selecting one piece of SAR image data as an optimal public main image based on a space-time base line threshold value and a Doppler mass center frequency threshold value, and then taking other q-1 pieces of SAR images as auxiliary images; q-1 auxiliary images after registration are obtained through registering each auxiliary image with the main image;
inputting q SAR images into GAMMA software, setting a short baseline space-time baseline threshold value in the GAMMA software to form a z pair short baseline space-time interference pair, and performing differential interference treatment on the generated short baseline space-time interference pair to generate a z short baseline differential interference diagram; wherein, the short baseline time-space baseline threshold time baseline is set in the GAMMA software and is less than or equal to 60 days, and the space baseline is less than or equal to 400 meters.
3. The method for extracting the surface deformation by fusing PCA and time sequence InSAR as set forth in claim 2, wherein the process of generating the initial InSAR time sequence result is as follows:
acquiring high-coherence points in q SAR images on different dates based on amplitude dispersion and phase stability principles by using StaMPS software, and arranging all the high-coherence points according to a time sequence to form a data set, wherein the data set comprises p high-coherence points;
carrying out phase unwrapping on the z differential interferograms by adopting a 3D unwrapping method, and estimating and removing phase errors caused by DEM and space-time base lines to obtain an initial InSAR time sequence result W p×q 。
4. The method for extracting the surface deformation by fusing PCA and time sequence InSAR as claimed in claim 3, wherein the PCA is adopted to process the initial InSAR time sequence result W by principal component analysis p×q Generating a plurality of main components, wherein the specific steps are as follows:
1) Centralizing an initial InSAR time sequence result:
initial InSAR time series result W p×q ＝{w 1 ,w 2 ,w 3 ,……,w q -w is i High-coherence point deformation information, w, of the ith image relative to the first image i Mean value of mu i Wherein i=1, 2, … …, q, w' i ＝w i -μ i I.e. each bit of the features minus the respective average value, the centred matrix is W' p×q ；
2) Centering the matrix W 'through SVD algorithm' p×q Preliminary decomposition:
wherein R is p×p Andis a standard orthogonal matrix, sigma p×q Is a diagonal matrix, where each term is a singular value of the matrix;
and carrying out data dimension reduction on the preliminary decomposition result, trying different dimension reduction dimensions v to generate a main component decomposition scheme, comparing the results of different schemes, synthesizing data quantity and the space-time rule recognition characteristics of each component, and selecting the optimal scheme.
5. The method for extracting the surface deformation by fusing PCA and time sequence InSAR as set forth in claim 4, wherein the method for selecting the optimal scheme by the identification features of the space-time law of the comprehensive data volume is characterized by comprising the following specific steps:
1a) For W 'and W' T Matrix multiplication is carried out to obtain a p multiplied by p dimensional square matrix W' T Wherein W 'is the matrix W' p×q For the initial InSAR time series result W p×q Results after the centering, W' T For the transpose of W', the eigenvalues and eigenvectors of the square matrix are solved using:
(W′W′ T )r j ＝λ j r j
where j=1, 2, &.. p is the number of high coherence points, r j As a feature vector lambda j The characteristic value is the characteristic value corresponding to the characteristic vector;
at this time sigma j Is singular value, r j As left singular vector, singular value sigma j Arranged from large to small in a matrix, σ j The decrease is fast, so the matrix is described approximately by the first v large singular values, then
1b) For W' T And W 'to obtain a q multiplied by q matrix W' T W', solving eigenvalues and eigenvectors of the square matrix:
(W′ T W′)t i ＝η i t i
wherein (i=1, 2,) q, q being the number of image frames, t i As a feature vector, eta i The characteristic value is the characteristic value corresponding to the characteristic vector;
mu at this time i Is singular value, t i Singular value μ as right singular vector i Arranged from large to small in the matrix and decreasing rapidly, so the matrix is described approximately by the first v large singular values, then
1c) W' T And W' T The singular value of W' is square root to form S v×v ；
The SVD algorithm data dimension-reducing matrix W 'is obtained through the steps' p×q The following results:
wherein v is the number of principal components in the principal component space, W' p×q R is the initial result of the time sequence InSAR after centering p×v And T v×q Respectively a space matrix and a time matrix, S v×v For a diagonal matrix, the value of which represents the importance of each component,is T v×q Is a transpose of (a).
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310066005.3A CN116363057B (en) | 2023-01-16 | 2023-01-16 | Surface deformation extraction method integrating PCA and time sequence InSAR |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310066005.3A CN116363057B (en) | 2023-01-16 | 2023-01-16 | Surface deformation extraction method integrating PCA and time sequence InSAR |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN116363057A CN116363057A (en) | 2023-06-30 |
| CN116363057B true CN116363057B (en) | 2023-11-10 |
Family
ID=86910230
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202310066005.3A Active CN116363057B (en) | 2023-01-16 | 2023-01-16 | Surface deformation extraction method integrating PCA and time sequence InSAR |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN116363057B (en) |
Citations (11)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| ITRM20130426A1 (en) * | 2013-07-19 | 2015-01-20 | Consiglio Nazionale Ricerche | METHOD FOR FILTERING INTERFEROMETRIC DATA ACQUIRED BY SYNTHETIC OPENING RADAR (SAR). |
| CN110058237A (en) * | 2019-05-22 | 2019-07-26 | 中南大学 | InSAR point Yun Ronghe and three-dimensional deformation monitoring method towards High-resolution SAR Images |
| CN112986993A (en) * | 2021-02-07 | 2021-06-18 | 同济大学 | InSAR deformation monitoring method based on space constraint |
| CN113175949A (en) * | 2021-04-23 | 2021-07-27 | 首都师范大学 | Method and system for inverting water release coefficient by combining surface deformation and water level information |
| CN113204023A (en) * | 2021-05-10 | 2021-08-03 | 中国地质大学（武汉） | Dual-polarization phase optimization earth surface deformation monitoring method combining PS target and DS target |
| WO2021227423A1 (en) * | 2020-05-13 | 2021-11-18 | 深圳大学 | Insar digital elevation model construction method and system based on dynamic baseline |
| CN113687353A (en) * | 2021-08-09 | 2021-11-23 | 中国矿业大学 | DS target phase optimization method based on homogeneous pixel time sequence phase matrix decomposition |
| CN114415131A (en) * | 2022-01-14 | 2022-04-29 | 长安大学 | Satellite-borne InSAR (interferometric synthetic Aperture Radar) atmospheric correction method combined with GACOS (gas-assisted operating System) |
| CN114440758A (en) * | 2022-01-09 | 2022-05-06 | 西北大学 | Analysis method for response of landslide to rainfall on regional scale |
| WO2022118101A1 (en) * | 2020-12-04 | 2022-06-09 | Rezatec Limited | System and method for remote dam monitoring |
| CN114693819A (en) * | 2022-04-15 | 2022-07-01 | 昆明理工大学 | Time sequence InSAR image time-space dimension reduction compression method based on tensor decomposition |
-
2023
- 2023-01-16 CN CN202310066005.3A patent/CN116363057B/en active Active
Patent Citations (11)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| ITRM20130426A1 (en) * | 2013-07-19 | 2015-01-20 | Consiglio Nazionale Ricerche | METHOD FOR FILTERING INTERFEROMETRIC DATA ACQUIRED BY SYNTHETIC OPENING RADAR (SAR). |
| CN110058237A (en) * | 2019-05-22 | 2019-07-26 | 中南大学 | InSAR point Yun Ronghe and three-dimensional deformation monitoring method towards High-resolution SAR Images |
| WO2021227423A1 (en) * | 2020-05-13 | 2021-11-18 | 深圳大学 | Insar digital elevation model construction method and system based on dynamic baseline |
| WO2022118101A1 (en) * | 2020-12-04 | 2022-06-09 | Rezatec Limited | System and method for remote dam monitoring |
| CN112986993A (en) * | 2021-02-07 | 2021-06-18 | 同济大学 | InSAR deformation monitoring method based on space constraint |
| CN113175949A (en) * | 2021-04-23 | 2021-07-27 | 首都师范大学 | Method and system for inverting water release coefficient by combining surface deformation and water level information |
| CN113204023A (en) * | 2021-05-10 | 2021-08-03 | 中国地质大学（武汉） | Dual-polarization phase optimization earth surface deformation monitoring method combining PS target and DS target |
| CN113687353A (en) * | 2021-08-09 | 2021-11-23 | 中国矿业大学 | DS target phase optimization method based on homogeneous pixel time sequence phase matrix decomposition |
| CN114440758A (en) * | 2022-01-09 | 2022-05-06 | 西北大学 | Analysis method for response of landslide to rainfall on regional scale |
| CN114415131A (en) * | 2022-01-14 | 2022-04-29 | 长安大学 | Satellite-borne InSAR (interferometric synthetic Aperture Radar) atmospheric correction method combined with GACOS (gas-assisted operating System) |
| CN114693819A (en) * | 2022-04-15 | 2022-07-01 | 昆明理工大学 | Time sequence InSAR image time-space dimension reduction compression method based on tensor decomposition |
Non-Patent Citations (4)
| Title |
|---|
| CAESAR: An Approach Based on Covariance Matrix Decomposition to Improve Multibaseline–Multitemporal Interferometric SAR Processing;Gianfranco Fornaro 等;《IEEE Transactions on Geoscience and Remote Sensing》;第53卷(第4期);全文 * |
| 基于地表变形数据的潜在滑坡识别研究;李萌;彭思佳;白艳萍;黄兆欢;;科技创新与应用(第06期);全文 * |
| 新疆阜康煤田火区DS-InSAR地表形变监测研究;彭锴;《中国优秀硕士学位论文全文数据库 (基础科学辑)》(第3期);全文 * |
| 面向对象高可信SAR数据精确处理;张继贤;黄国满;程春泉;;武汉大学学报(信息科学版)(第12期);全文 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN116363057A (en) | 2023-06-30 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Crippen et al. | NASADEM global elevation model: Methods and progress | |
| Pfeffer et al. | The key role of vertical land motions in coastal sea level variations: A global synthesis of multisatellite altimetry, tide gauge data and GPS measurements | |
| Liu et al. | Estimating Spatiotemporal Ground Deformation With Improved Persistent-Scatterer Radar Interferometry $^\ast$ | |
| Maubant et al. | Independent component analysis and parametric approach for source separation in InSAR time series at regional scale: application to the 2017–2018 Slow Slip Event in Guerrero (Mexico) | |
| Samsonov et al. | Multidimensional Small Baseline Subset (MSBAS) for volcano monitoring in two dimensions: Opportunities and challenges. Case study Piton de la Fournaise volcano | |
| CN112556660B (en) | Sea area gravity anomaly inversion method and system based on satellite height measurement data | |
| CN111398959A (en) | InSAR time sequence earth surface deformation monitoring method based on earth surface stress strain model | |
| CN111856459B (en) | Improved DEM maximum likelihood constraint multi-baseline InSAR phase unwrapping method | |
| CN108427741B (en) | DEM relative error evaluation method based on large number of high-precision control points | |
| Zhao et al. | Assessment of ASTER GDEM performance by comparing with SRTM and ICESat/GLAS data in Central China | |
| CN113281749A (en) | Time sequence InSAR high-coherence point selection method considering homogeneity | |
| CN112051572A (en) | Method for monitoring three-dimensional surface deformation by fusing multi-source SAR data | |
| CN115079172A (en) | MTInSAR landslide monitoring method, equipment and storage medium | |
| Zhang et al. | Reduction of atmospheric effects on InSAR observations through incorporation of GACOS and PCA into small baseline subset InSAR | |
| Barrile et al. | Analysis of hydraulic risk territories: comparison between LIDAR and other different techniques for 3D modeling | |
| Duan et al. | Adaptively selecting interferograms for SBAS-InSAR based on graph theory and turbulence atmosphere | |
| CN114812491A (en) | Power transmission line earth surface deformation early warning method and device based on long-time sequence analysis | |
| Bacalhau et al. | Bathymetry of reservoirs using altimetric data associated to optical images | |
| Yaseen et al. | Local interpolation of coseismic displacements measured by InSAR | |
| CN116363057B (en) | Surface deformation extraction method integrating PCA and time sequence InSAR | |
| CN116068511B (en) | Deep learning-based InSAR large-scale system error correction method | |
| Hu et al. | Isolating orbital error from multitemporal InSAR derived tectonic deformation based on wavelet and independent component analysis | |
| Abili | Comparison of vertical accuracy of open-source global digital elevation models: a case study of Adama City, Ethiopia | |
| CN114046774B (en) | Ground deformation continuous monitoring method integrating CORS network and multi-source data | |
| CN113341410B (en) | Large-range under-forest terrain estimation method, device, equipment and medium |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |