CN111398959B - InSAR time sequence earth surface deformation monitoring method based on earth surface stress strain model - Google Patents
InSAR time sequence earth surface deformation monitoring method based on earth surface stress strain model Download PDFInfo
- Publication number
- CN111398959B CN111398959B CN202010265777.6A CN202010265777A CN111398959B CN 111398959 B CN111398959 B CN 111398959B CN 202010265777 A CN202010265777 A CN 202010265777A CN 111398959 B CN111398959 B CN 111398959B
- Authority
- CN
- China
- Prior art keywords
- deformation
- insar
- arc
- time sequence
- earth surface
- 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
- 238000000034 method Methods 0.000 title claims abstract description 44
- 238000012544 monitoring process Methods 0.000 title claims abstract description 33
- 238000010586 diagram Methods 0.000 claims abstract description 5
- 239000011159 matrix material Substances 0.000 claims description 33
- 238000012876 topography Methods 0.000 claims description 16
- 238000012937 correction Methods 0.000 claims description 9
- 238000006243 chemical reaction Methods 0.000 claims description 3
- 230000001419 dependent effect Effects 0.000 claims description 3
- 238000003384 imaging method Methods 0.000 claims description 3
- 230000010354 integration Effects 0.000 claims description 3
- 230000000873 masking effect Effects 0.000 abstract description 3
- 238000012545 processing Methods 0.000 abstract description 2
- 238000004088 simulation Methods 0.000 description 6
- 239000008186 active pharmaceutical agent Substances 0.000 description 5
- 238000011160 research Methods 0.000 description 4
- 230000006978 adaptation Effects 0.000 description 2
- 238000004422 calculation algorithm Methods 0.000 description 2
- 238000007796 conventional method Methods 0.000 description 2
- 238000001914 filtration Methods 0.000 description 2
- 238000012986 modification Methods 0.000 description 2
- 230000004048 modification Effects 0.000 description 2
- 238000010587 phase diagram Methods 0.000 description 2
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000004364 calculation method Methods 0.000 description 1
- 230000001427 coherent effect Effects 0.000 description 1
- 230000001066 destructive effect Effects 0.000 description 1
- 238000001514 detection method Methods 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 230000008030 elimination Effects 0.000 description 1
- 238000003379 elimination reaction Methods 0.000 description 1
- 238000005305 interferometry Methods 0.000 description 1
- 238000005259 measurement Methods 0.000 description 1
- 238000005457 optimization Methods 0.000 description 1
- 230000002123 temporal effect Effects 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
- G01S13/00—Systems using the reflection or reradiation of radio waves, e.g. radar systems; Analogous systems using reflection or reradiation of waves whose nature or wavelength is irrelevant or unspecified
- G01S13/88—Radar or analogous systems specially adapted for specific applications
- G01S13/89—Radar or analogous systems specially adapted for specific applications for mapping or imaging
- G01S13/90—Radar or analogous systems specially adapted for specific applications for mapping or imaging using synthetic aperture techniques, e.g. synthetic aperture radar [SAR] techniques
- G01S13/9021—SAR image post-processing techniques
- G01S13/9023—SAR image post-processing techniques combined with interferometric techniques
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01B—MEASURING LENGTH, THICKNESS OR SIMILAR LINEAR DIMENSIONS; MEASURING ANGLES; MEASURING AREAS; MEASURING IRREGULARITIES OF SURFACES OR CONTOURS
- G01B7/00—Measuring arrangements characterised by the use of electric or magnetic techniques
- G01B7/16—Measuring arrangements characterised by the use of electric or magnetic techniques for measuring the deformation in a solid, e.g. by resistance strain gauge
Landscapes
- Engineering & Computer Science (AREA)
- Remote Sensing (AREA)
- Radar, Positioning & Navigation (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Electromagnetism (AREA)
- Computer Networks & Wireless Communication (AREA)
- Radar Systems Or Details Thereof (AREA)
- Image Analysis (AREA)
Abstract
The invention provides an InSAR time sequence earth surface deformation monitoring method based on an earth surface stress strain model, which comprises the following steps: step 1, collecting a time sequence SAR image covering a region to be monitored, realizing registration and differential interference of the time sequence SAR image based on the existing software, and generating an interference pattern meeting a preset time-space baseline threshold; and 2, constructing a local Dirony triangle network based on all pixel points in the interference diagram. According to the invention, the physical and mechanical relation between the deformation of the adjacent points of the earth surface is described by using the earth surface stress strain model, so that the accuracy of InSAR time sequence earth surface deformation monitoring can be remarkably improved, modeling and parameter solving are performed based on arc segments formed by the adjacent points, therefore, the interference pattern is not required to be unwrapped, the InSAR data processing efficiency can be effectively improved, the physical and mechanical relation of the observed data in space is considered, and therefore, the high-accuracy InSAR deformation monitoring result with large space coverage density can be obtained without masking a low-coherence region in the interference pattern.
Description
Technical Field
The invention relates to the field of earth measurement based on satellite interference synthetic aperture radar, in particular to an InSAR time sequence earth surface deformation monitoring method based on an earth surface stress strain model.
Background
Synthetic aperture radar interferometry (Interferometric Synthetic Aperture Radar, SAR, inSAR) techniques utilize two-shot SAR images to measure surface deformation. The time sequence InSAR technology (Multiple Temporal InSAR, MTInSAR) analyzes a series of SAR images on time sequence, so that time sequence earth surface deformation results of a research area can be obtained, and meanwhile, the influence of factors such as atmospheric delay, topography residual errors, destructive interference errors and the like on the earth surface deformation results can be effectively weakened. The monitoring targets of the MTInSAR technology are largely divided into permanent scatterers (Permanent Scatterer, PS) and distributed scatterers (Distributed Scatterer, DS). PS keeps higher signal to noise ratio in the whole detection period, but the distribution density of PS is extremely low in natural environment, so that the monitoring requirement is difficult to meet. DS keeps better signal-to-noise ratio on some short-time space baseline interferograms, the space-time variation characteristics of adjacent pixels are similar, DS distribution density is higher (such as sparse vegetation areas, deserts and the like) in natural environment, and application is wider.
When acquiring the time sequence earth surface deformation at the DS target, the existing method often needs to perform complex and complicated mathematical optimization estimation (namely phase reconstruction) based on all possible interference pairs, or perform point-by-point estimation based on a least square criterion by utilizing the unwrapped short-time blank baseline interferogram. The former has lower efficiency in estimating the surface deformation, while the latter is easy to be influenced by factors such as unwrapping error due to the need of unwrapping in advance. In addition, the existing method does not consider the physical and mechanical relationship between adjacent DS when estimating time sequence earth surface deformation, and can not obtain more accurate deformation results in a low coherence area.
Disclosure of Invention
The invention provides an InSAR time sequence earth surface deformation monitoring method based on an earth surface stress strain model, which aims to solve the problems that the traditional method only carries out mathematical spatial filtering on high coherence points, or unwraps an interference pattern, masks and then carries out time sequence solving point by point, does not consider the physical and mechanical relationship between adjacent high coherence points, and is easy to cause lower spatial density of monitoring results.
In order to achieve the above object, an embodiment of the present invention provides an InSAR time-series earth surface deformation monitoring method based on an earth surface stress strain model, including:
step 3, establishing a functional relation between time sequence deformation phase gradient of the target arc section and interference phase differences of all arc sections within a surrounding range of 1km multiplied by 1km based on a surface stress strain model;
and 5, repeating the step 3 and the step 4 for each side of the local Dirony triangulation network until the time sequence deformation phase difference on all the arc sections is calculated, and performing space integration on all the arc sections by taking one earth surface stable point or known deformation point in the triangulation network as a reference to obtain a time sequence deformation result.
Wherein, the step 1 specifically includes:
and collecting the T Jing Shixu SAR image covering the area to be monitored, realizing the registration and differential interference of the time sequence SAR image based on the existing software, and finally generating M interferograms meeting the preset time-space baseline threshold.
Wherein, the step 2 specifically includes:
based on all pixel points in the M interferograms, a local Dirony triangle network is constructed, and the local Dirony triangle network contains
The arc segments, the longest arc segment, should be less than the preset threshold of arc segment length, which is set according to the range of the assumed homogeneous region.
Wherein, the step 2 further comprises:
assume a certain arc segment A in the constructed local Dirony triangulation network 0 Is in the form of point P i Sum point P j To be the start and stop point, point P i Sum point P j The corresponding three-dimensional coordinates are respectively
And->
Point P i Sum point P j Three-dimensional earth surface deformation occurring between SAR image acquisition moments corresponding to the mth interferogram is respectively
And->
According to the earth surface stress strain model, the following steps are obtained:
wherein,,
e represents east-west direction, n represents north-south direction, u represents vertical direction,/->
An unknown parameter matrix representing a surface stress strain model can be expressed as:
wherein,,
represents the partial derivative, d m ＝[d m，e d m，n d m，u ] T Represents the three-dimensional earth surface deformation field which occurs between the acquisition moments of the primary SAR image and the auxiliary SAR image corresponding to the mth interferogram, and x= [ x ] e x n x u ] T Representing the three-dimensional directions of east-west, north-south and vertical.
Wherein, the step 2 further comprises:
due to point P i Sum point P j The distance is relatively short, and the geometrical view angle difference of the SAR satellite when observing the two points is negligible, namely at the point P i Sum point P j At this point, three-dimensional deformation is projected to coefficient matrix B of InSAR line of sight deformation geo Is identical, coefficient matrix B geo The following is shown:
B geo ＝[a b c] (3)
wherein a represents a projection coefficient of deformation of east-west deformation in the InSAR sight line direction, b represents a projection coefficient of deformation of north-south deformation in the InSAR sight line direction, c represents a projection coefficient of deformation of vertical deformation in the InSAR sight line direction, and the projection coefficient can be determined according to the imaging geometry of SAR satellites;
multiplying B by the same sign of (1) geo The following formula is obtained:
wherein,,
respectively represent the point P i Sum point P j Projection deformation of the three-dimensional earth deformation along the InSAR line of sight, i.e. the observation at two points in the mth interferogram, additionally +.>
Can be regarded as the one-dimensional deformation observation value of InSAR in the arc section A 0 And the deformation gradient parameters along east-west direction, south-north direction and vertical direction are set.
Wherein, the step 3 specifically includes:
based on the earth surface stress strain model, in the mth interferogram, the arc section A 0 Surface deformation gradient in certain space range
Can be assumed to be constant, thus, based on the equation (4), the target arc segment A in the mth interferogram can be established 0 Upper surface deformation gradient->
All K arc sections A within a certain range of the surrounding k Comprising arc section A 0 Is->
A functional relationship between (k=1, 2,., K), resulting in:
wherein,,
B sm ＝[Δ 1 ，Δ 2 ，...，Δ k ，…，Δ K ] T ，Δ k representing arc segment A k Coordinate increment between two endpoints, +.>
Representing arc segment A k Difference in InSAR deformation observations between the two endpoints.
Wherein, the step 3 further comprises:
taking a first scene SAR image of the M interferograms as a reference, and easily obtaining a conversion matrix B between deformation observed values corresponding to the InSAR interferograms and deformation of the InSAR sight line to time sequence ground surface trs Wherein B is trs The matrix size of the matrix is M x (T-1), each row corresponds to one interferogram, each column corresponds to one SAR image, the reference image is removed, in each row, the columns corresponding to the M interferogram main images are-1, the columns corresponding to the auxiliary images are +1, and other elements are 0;
based on (5) and matrix B trs M interferograms can be establishedThe phase difference DeltaL of deformation on all arc sections defo And arc section A 0 Time sequence earth surface deformation phase gradient
Functional relationship between:
wherein DeltaL defo ＝[(ΔL 1 ) T ，(ΔL 2 ) T ，…，(ΔL m ) T ，…，(ΔL M ) T ] T ，
Wherein t=2, 3, T,
representing the sign of the Cronecker product, < ->
Know->
Respectively representing InSAR line-of-sight deformation phase at arc segment A during the acquisition of the t-th SAR image 0 Gradients in east-west, north-south, and vertical directions;
target arc segment A 0 The actual InSAR observed value comprises the earth surface deformation phase difference
Including topography residual dependent phase->
And noise phase, wherein the target arc segment A 0 Topography residual phase->
Expressed as:
wherein,,
wherein,,
representing arc segment A in the mth interferogram 0 Topography residual phase on ∈>
Vertical baseline representing mth interferogram, dz 0 Representing arc segment A 0 The terrain residual, λ, ρ, θ on the upper surface represent the satellite wavelength, satellite-to-target point distance, satellite incident angle at the target point, respectively.
Wherein, the step 3 further comprises:
assuming that the terrain residual correlation phase is spatially uncorrelated, arc segment A k The residual phase of the topography above can be considered as noise, and is based on the strain of the surface stress by combining the formula (6) and the formula (7)The model can establish the functional relation between the time sequence deformation phase gradient of the target arc section and the interference phase difference delta L of all arc sections in a certain surrounding range:
wherein DeltaL is the interference phase on the corresponding K arc segments in all M interferograms, which is obtained by making the difference between the true observed phases of the corresponding arc segment endpoints, B= [ B ] defo ，B dz ]，
Matrix B dz The size of (M×K) ×1, matrix B dz Line number and B of (2) defo Consistent with DeltaL, matrix B dz And arc segment A in delta L 0 The elements of the positions corresponding to the upper interference phase difference observation values are +.>
The other elements are all 0.
Wherein, the step 3 further comprises:
weighting different arc segments according to the coherence of high coherence points, wherein the coherence is an InSAR observation value of c, and the corresponding variance sigma 2 Can be expressed as:
InSAR observed value variance of two endpoints of easily obtained component arc section
And->
Furthermore, the variance of the difference between the InSAR observations at two points can be obtained>
The method comprises the following steps:
ignoring the covariance of InSAR observations between different arc segments, obtaining a variance matrix D of an observation vector DeltaL based on the coherence values of the equation (9) and the equation (10) and all points ΔL Wherein D is ΔL The diagonal elements respectively correspond to the variances of the InSAR observed values delta L on the arc segments, the rest elements are 0, and the weight array P of the observed value vector delta L ΔL Can be expressed as:
P ΔL ＝(D ΔL ) -1 (11)。
wherein, the step 4 specifically includes:
based on equation (8) and equation (11), an initial solution of the unknown parameter vector can be solved based on a weighted least squares criterion
The method comprises the following steps:
in actual data, the difference between the neighboring point observations may contain 2pi ambiguity, i.e., the observation ΔL may contain gross errors, in order to reduce the gross errors to the unknown parameter vector
According to the statistical characteristics of the observed value, i.e. when the observed value correction is greater than a certain threshold v thr When the corresponding observed value is considered to be a gross error, where v thr Calculated according to the following formula:
wherein,,
observed value correction v 0 Calculated according to the following formula:
removing the correction v from DeltaL 0 ＞v thr Is re-calculated according to equation (12)
I.e. is not influenced by the rough difference, the arc section A can be obtained 0 Difference dz between topography residual errors between upper starting and stopping points 0 And a time-sequential earth deformation gradient in east-west, north-south and vertical directions +.>
Obtaining arc section A 0 Two end points P i Sum point P j Time sequence earth surface deformation phase difference between
The scheme of the invention has the following beneficial effects:
according to the InSAR time sequence earth surface deformation monitoring method based on the earth surface stress strain model, the earth surface stress strain model is utilized to describe the physical mechanical relation between earth surface adjacent point deformation, the accuracy of InSAR time sequence earth surface deformation monitoring can be remarkably improved, modeling and parameter solving are carried out based on arc segments formed by adjacent points, therefore, an interferogram is not required to be unwrapped, inSAR data processing efficiency can be effectively improved, the physical mechanical relation of observation data in space is considered, and therefore, a high-accuracy InSAR deformation monitoring result with large space coverage density can be obtained without masking a low coherence region in the interferogram.
Drawings
FIG. 1 is a flow chart of the present invention;
FIG. 2 is a time-series earth surface deformation graph of the simulated SAR image acquisition moment of the invention;
FIG. 3 is an atmospheric delay phase diagram of the simulated SAR image acquisition time of the present invention;
FIG. 4 is a spatiotemporal baseline distribution diagram of a simulated SAR image of the present invention;
FIG. 5 is a simulated short baseline InSAR interferogram of the present invention;
FIG. 6 is a graph of the topography residual phase in a simulated InSAR interferogram of the present invention;
FIG. 7 is a phase diagram of the phase loss of coherent noise in a simulated InSAR interferogram of the present invention;
FIG. 8 is a graph of the average deformation rate of the original simulated interferograms of the present invention as solved by the methods of the present invention and by the conventional methods of (a), (c) and (b) with (b) unmasked and (c) and (d) masked.
Detailed Description
In order to make the technical problems, technical solutions and advantages to be solved more apparent, the following detailed description will be given with reference to the accompanying drawings and specific embodiments.
The invention provides an InSAR time sequence earth surface deformation monitoring method based on an earth surface stress strain model, which aims at solving the problems that the existing traditional method only carries out mathematical spatial filtering on high coherence points, or unwrapping and masking interference patterns and then carrying out time sequence solving point by point, does not consider the physical and mechanical relationship between adjacent high coherence points and easily causes lower space density of monitoring results.
As shown in fig. 1 to 8, an embodiment of the present invention provides an InSAR time-series earth surface deformation monitoring method based on an earth surface stress strain model, including: step 1, collecting a time sequence SAR image covering a region to be monitored, realizing registration and differential interference of the time sequence SAR image based on the existing software, and generating an interference pattern meeting a preset time-space baseline threshold; step 2, constructing a local Dirony triangle network based on all pixel points in the interference diagram; step 3, establishing a functional relation between time sequence deformation phase gradient of the target arc section and interference phase differences of all arc sections within a surrounding range of 1km multiplied by 1km based on a surface stress strain model; step 4, introducing robust estimation to remove arc segments containing 2 pi ambiguity, and finally realizing the solution of the time sequence deformation phase difference of the target arc segments based on a weighted least square criterion; and 5, repeating the step 3 and the step 4 for each side of the local Dirony triangulation network until the time sequence deformation phase difference on all the arc sections is calculated, and performing space integration on all the arc sections by taking one earth surface stable point or known deformation point in the triangulation network as a reference to obtain a time sequence deformation result.
According to the InSAR time sequence earth surface deformation monitoring method based on the earth surface stress strain model, which is disclosed by the embodiment of the invention, the earth surface stress strain model describes the physical and mechanical relationship between the three-dimensional earth surface deformation of the earth surface adjacent points, and the function relationship between the InSAR one-dimensional time sequence deformation phase gradient of the target arc section and the interference phase difference of all arc sections within the range of 1km multiplied by 1km around can be established based on the earth surface stress strain model through the geometric relationship between the InSAR observation value and the three-dimensional earth surface deformation.
Wherein, the step 1 specifically includes: and collecting the T Jing Shixu SAR image covering the area to be monitored, realizing the registration and differential interference of the time sequence SAR image based on the existing software, and finally generating M interferograms meeting the preset time-space baseline threshold.
According to the InSAR time sequence earth surface deformation monitoring method based on the earth surface stress strain model, the setting of the time-space base line threshold value can be determined according to earth surface change conditions of specific research areas, orbit control of new generation SAR satellites (such as sentinel No. 1 satellites) is good, the space base line of time sequence SAR images of the same research area is often far smaller than the space incoherent base line threshold value, therefore, the influence of the setting of the space base line threshold value on the selection of interferograms is small, the time base line threshold value can be determined empirically according to vegetation density degree of the research area, SAR satellite wavelength and other factors, and generally, the time base line can be set to be 100 days, and the space base line can be set to be 200 meters.
Wherein, the step 2 specifically includes: based on all pixel points in the M interferograms, a local Dirony triangle network is constructed, and the local Dirony triangle network contains
The arc segments, the longest arc segment, should be less than the preset threshold of arc segment length, which is set according to the range of the assumed homogeneous region.
According to the InSAR time sequence earth surface deformation monitoring method based on the earth surface stress strain model, the preset threshold value of the arc length is generally selected according to the range of the homogeneous region, under the general condition, the preset threshold value of the arc length can be set to be 1km by considering that the range of the homogeneous region is 1km multiplied by 1km, at the moment, all two points meeting the threshold value are connected into a line, and the number of the last connecting lines is
Wherein, the step 2 further comprises: assume a certain arc segment A in the constructed local Dirony triangulation network 0 Is in the form of point P i Sum point P j To be the start and stop point, point P i Sum point P j The corresponding three-dimensional coordinates are respectively
And
point P i Sum point P j Three-dimensional surface deformation occurring between SAR image acquisition moments corresponding to the mth interferogram is +.>
And->
According to the earth surface stress strain model, the following steps are obtained:
wherein,,
e represents east-west direction, n represents north-south direction, u represents vertical direction,/->
An unknown parameter matrix representing a surface stress strain model can be expressed as:
wherein,,
represents the partial derivative, d m ＝[d m，e d m，n d m，u ] T Represents the three-dimensional earth surface deformation field which occurs between the acquisition moments of the primary SAR image and the auxiliary SAR image corresponding to the mth interferogram, and x= [ x ] e x n x u ] T Representing the three-dimensional directions of east-west, north-south and vertical.
Wherein, the step 2 further comprises: due to point P i Sum point P j The distance is relatively short, and the geometrical view angle difference of the SAR satellite when observing the two points is negligible, namely at the point P i Sum point P j At this point, three-dimensional deformation is projected to coefficient matrix B of InSAR line of sight deformation geo Is identical, coefficient matrix B geo The following is shown:
B geo ＝[a b c] (3)
wherein a represents a projection coefficient of deformation of east-west deformation in the InSAR sight line direction, b represents a projection coefficient of deformation of north-south deformation in the InSAR sight line direction, c represents a projection coefficient of deformation of vertical deformation in the InSAR sight line direction, and the projection coefficient can be determined according to the imaging geometry of SAR satellites;
multiplying B by the same sign of (1) geo The following formula is obtained:
wherein,,
respectively represent the point P i Sum point P j Projection deformation of the three-dimensional earth deformation along the InSAR line of sight, i.e. the observation at two points in the mth interferogram, additionally +.>
Can be regarded as the one-dimensional deformation observation value of InSAR in the arc section A 0 And the deformation gradient parameters along east-west direction, south-north direction and vertical direction are set.
Wherein, the step 3 specifically includes: based on the earth surface stress strain model, in the mth interferogram, the arc section A 0 Surface deformation gradient in certain space range
Can be assumed to be constant, thus, based on the equation (4), the target arc segment A in the mth interferogram can be established 0 Upper surface deformation gradient->
All K arc sections A within a certain range of the surrounding k Comprising arc section A 0 Is->
Functional relation betweenWherein, (k=1, 2,) K), yields:
wherein,,
B sm ＝[Δ 1 ，Δ 2 ，…，Δ k ，…，Δ K ] T ，Δ k representing arc segment A k Coordinate increment between two endpoints, +.>
Representing arc segment A k Difference in InSAR deformation observations between the two endpoints.
Wherein, the step 3 further comprises: taking a first scene SAR image of the M interferograms as a reference, and easily obtaining a conversion matrix B between deformation observed values corresponding to the InSAR interferograms and deformation of the InSAR sight line to time sequence ground surface trs Wherein B is trs The matrix size of the matrix is M x (T-1), each row corresponds to one interferogram, each column corresponds to one SAR image, the reference image is removed, in each row, the columns corresponding to the M interferogram main images are-1, the columns corresponding to the auxiliary images are +1, and other elements are 0;
based on (5) and matrix B trs The deformation phase difference delta L of all arc sections in the M interferograms can be established defo And arc section A 0 Time sequence earth surface deformation phase gradient
Functional relationship between:
wherein DeltaL defo ＝[(ΔL 1 ) T ，(ΔL 2 ) T ，…，(ΔL m ) T ，…，(ΔL M ) T ] T ，
Wherein t=2, 3, T,
representing the sign of the Cronecker product, < ->
And->
Respectively representing InSAR line-of-sight deformation phase at arc segment A during the acquisition of the t-th SAR image 0 Gradients in east-west, north-south, and vertical directions;
target arc segment A 0 The actual InSAR observed value comprises the earth surface deformation phase difference
Including topography residual dependent phase->
And noise phase, wherein the target arc segment A 0 Topography residual phase->
Expressed as:
wherein the method comprises the steps of，
Wherein,,
representing arc segment A in the mth interferogram 0 Topography residual phase on ∈>
Vertical baseline representing mth interferogram, dz 0 Representing arc segment A 0 The terrain residual, λ, ρ, θ on the upper surface represent the satellite wavelength, satellite-to-target point distance, satellite incident angle at the target point, respectively.
Wherein, the step 3 further comprises: assuming that the terrain residual correlation phase is spatially uncorrelated, arc segment A k The terrain residual phase can be regarded as noise, and the functional relation between the time sequence deformation phase gradient of the target arc section and the interference phase difference delta L of all arc sections in a certain surrounding range can be established based on the earth surface stress strain model by combining the formula (6) and the formula (7):
wherein DeltaL is the interference phase on the corresponding K arc segments in all M interferograms, which is obtained by making the difference between the true observed phases of the corresponding arc segment endpoints, B= [ B ] defo ，B dz ]，
Matrix B dz The size of (M×K) ×1, matrix B dz Line number and B of (2) defo Consistent with DeltaL, matrix B dz And arc segment A in delta L 0 The elements of the positions corresponding to the upper interference phase difference observation values are +.>
The other elements are all 0.
Wherein, the step 3 further comprises: weighting different arc segments according to the coherence of high coherence points, wherein the coherence is an InSAR observation value of c, and the corresponding variance sigma 2 Can be expressed as:
InSAR observed value variance of two endpoints of easily obtained component arc section
And->
Furthermore, the variance of the difference between the InSAR observations at two points can be obtained>
The method comprises the following steps:
ignoring the covariance of InSAR observations between different arc segments, obtaining a variance matrix D of an observation vector DeltaL based on the coherence values of the equation (9) and the equation (10) and all points ΔL Wherein D is ΔL The diagonal elements respectively correspond to the variances of the InSAR observed values delta L on the arc segments, the rest elements are 0, and the weight array P of the observed value vector delta L ΔL Can be expressed as:
P ΔL ＝(D ΔL ) -1 (11)。
wherein, the step 4 specifically includes: base groupIn equations (8) and (11), an initial solution for the unknown parameter vector can be solved based on a weighted least squares criterion
The method comprises the following steps:
in actual data, the difference between the neighboring point observations may contain 2pi ambiguity, i.e., the observation ΔL may contain gross errors, in order to reduce the gross errors to the unknown parameter vector
According to the statistical characteristics of the observed value, i.e. when the observed value correction is greater than a certain threshold v thr When the corresponding observed value is considered to be a gross error, where v thr Calculated according to the following formula:
wherein,,
observed value correction v 0 Calculated according to the following formula:
removing the correction v from DeltaL 0 ＞v thr Is re-calculated according to equation (12)
I.e. is not influenced by the rough difference, the arc section A can be obtained 0 Difference dz0 between topography residual errors between starting and stopping points and time sequence surface shape along east-west direction, north-south direction and vertical directionGradient->
Obtaining arc section A 0 Two end points P i Sum point P j Time sequence earth surface deformation phase difference between
The effect of the InSAR time sequence earth surface deformation monitoring method based on the earth surface stress strain model according to the embodiment of the invention can be further illustrated by the following simulation experiment, and the simulation data is described: (1) the earth surface deformation field of 10 moments of linear change is simulated in a certain area (the image size is 100 multiplied by 100) as shown in figure 2, and the corresponding atmospheric phases of 10 moments are simulated based on a fractal curve surface (the fractal dimension is 2.67) as shown in figure 3, wherein the maximum value of the atmospheric phases is 1 radian; (2) in combination with the sentry-1A/B satellite data space-time baselines as shown in fig. 4, 30 interferograms were generated as shown in fig. 5. Wherein, the interference pattern is added with the topography residual phases which are uniformly distributed in the interval [ -10m,10m ] as shown in figure 6, and the phase loss interference noise based on the real interference pattern coherence simulation as shown in figure 7.
When the time sequence earth surface deformation of the target arc section is solved, firstly, a terrain residual error is solved by using a presumed deformation model (such as a linear deformation model), then the phase of the terrain residual error is subtracted from an interferogram, and the time sequence earth surface deformation is solved by using the residual phase, wherein in the whole process, only the observation phase on the target arc section is considered, and the observation values on other arc sections in a certain range around are not considered.
The simulation experiment utilizes the method and the traditional method to solve the simulation data under the conditions that the original simulation interferogram is not masked and the pixels with average coherence lower than 0.6 are masked. In the solving result, because the method considers the physical and mechanical relation between the deformation of the surface adjacent points and adopts robust estimation to perform robust elimination, the method can obtain accurate time sequence surface deformation results as shown in fig. 8 (a) and (c) even in a low coherence area without a mask. In contrast, the conventional method can obtain more reliable deformation results only in the high coherence region as shown in fig. 8 (b) and (d). On this basis, the root mean square error of the time-series deformation obtained by the two methods is shown in table 1:
TABLE 1 root mean square error for InSAR timing deformation
As can be seen from a combination of table 1 and fig. 8, the algorithm of the present invention can obtain more accurate and complete spatial coverage of the InSAR time sequence earth surface deformation field compared with the conventional algorithm.
According to the InSAR time sequence earth surface deformation monitoring method based on the earth surface stress strain model, the earth surface stress strain model with a physical and mechanical meaning is introduced to calculate time sequence earth surface deformation, so that the calculation precision and the space coverage density of the time sequence earth surface deformation are further improved, meanwhile, interference pattern unwrapping is not needed in the InSAR time sequence earth surface deformation monitoring method based on the earth surface stress strain model, and time sequence earth surface deformation estimation errors caused by interference pattern unwrapping errors are effectively weakened.
While the foregoing is directed to the preferred embodiments of the present invention, it will be appreciated by those skilled in the art that various modifications and adaptations can be made without departing from the principles of the present invention, and such modifications and adaptations are intended to be comprehended within the scope of the present invention.
Claims (10)
1. An InSAR time sequence earth surface deformation monitoring method based on an earth surface stress strain model is characterized by comprising the following steps:
step 1, collecting a time sequence SAR image covering a region to be monitored, realizing registration and differential interference of the time sequence SAR image based on the existing software, and generating an interference pattern meeting a preset time-space baseline threshold;
step 2, constructing a local Dirony triangle network based on all pixel points in the interference diagram;
step 3, establishing a functional relation between time sequence deformation phase gradient of the target arc section and interference phase differences of all arc sections within a surrounding range of 1km multiplied by 1km based on a surface stress strain model;
step 4, introducing robust estimation to remove arc segments containing 2 pi ambiguity, and finally realizing the solution of the time sequence deformation phase difference of the target arc segments based on a weighted least square criterion;
and 5, repeating the step 3 and the step 4 for each side of the local Dirony triangulation network until the time sequence deformation phase difference on all the arc sections is calculated, and performing space integration on all the arc sections by taking one earth surface stable point or known deformation point in the triangulation network as a reference to obtain a time sequence deformation result.
2. The method for monitoring the surface deformation of the InSAR chronology based on the surface stress strain model according to claim 1, wherein the step 1 specifically comprises the following steps:
and collecting the T Jing Shixu SAR image covering the area to be monitored, realizing the registration and differential interference of the time sequence SAR image based on the existing software, and finally generating M interferograms meeting the preset time-space baseline threshold.
3. The method for monitoring the surface deformation of the InSAR chronology based on the surface stress strain model according to claim 2, wherein the step 2 specifically comprises:
based on all pixel points in the M interferograms, a local Dirony triangle network is constructed, and the local Dirony triangle network contains
The arc segments, the longest arc segment should be smaller than the preset threshold of the arc segment length, the preset threshold of the arc segment length is set according to the range of the assumed homogeneous regionAnd is defined around.
4. The method for monitoring the surface deformation of the InSAR chronology based on the surface stress strain model as set forth in claim 3, wherein the step 2 further comprises:
assume a certain arc segment A in the constructed local Dirony triangulation network 0 Is in the form of point P i Sum point P j To be the start and stop point, point P i Sum point P j The corresponding three-dimensional coordinates are respectively
And->
Point P i Sum point P j Three-dimensional surface deformation occurring between SAR image acquisition moments corresponding to the mth interferogram is +.>
And->
According to the earth surface stress strain model, the following steps are obtained:
wherein,,
e represents east-west direction, n represents north-south direction, u represents vertical direction,/->
An unknown parameter matrix representing a surface stress strain model can be expressed as:
wherein,,
represents the partial derivative, d m ＝[d m,e d m,n d m,u ] T Represents the three-dimensional earth surface deformation field which occurs between the acquisition moments of the primary SAR image and the auxiliary SAR image corresponding to the mth interferogram, and x= [ x ] e x n x u ] T Representing the three-dimensional directions of east-west, north-south and vertical.
5. The method for monitoring the surface deformation of the InSAR chronology based on the surface stress strain model as set forth in claim 4, wherein the step 2 further comprises:
due to point P i Sum point P j The geometrical view angle difference when the SAR satellite observes the two points is ignored when the distance is relatively close, namely at the point P i Sum point P j At this point, three-dimensional deformation is projected to coefficient matrix B of InSAR line of sight deformation geo Is identical, coefficient matrix B geo The following is shown:
B geo ＝[a b c] (3)
wherein a represents a projection coefficient of deformation of east-west deformation in the InSAR sight line direction, b represents a projection coefficient of deformation of north-south deformation in the InSAR sight line direction, c represents a projection coefficient of deformation of vertical deformation in the InSAR sight line direction, and the projection coefficient can be determined according to the imaging geometry of SAR satellites;
multiplying B by the same sign of (1) geo The following formula is obtained:
wherein,,
respectively represent the point P i Sum point P j Projection deformation of the three-dimensional earth deformation along the InSAR line of sight, i.e. the observation at two points in the mth interferogram, additionally +.>
Can be regarded as the one-dimensional deformation observation value of InSAR in the arc section A 0 And the deformation gradient parameters along east-west direction, south-north direction and vertical direction are set.
6. The method for monitoring the surface deformation of the InSAR chronology based on the surface stress strain model according to claim 5, wherein the step 3 specifically comprises:
based on the earth surface stress strain model, in the mth interferogram, the arc section A 0 Surface deformation gradient in certain space range
Assuming constant, therefore, based on equation (4), the target arc A in the mth interferogram can be established 0 Upper surface deformation gradient->
All K arc sections A within a certain range of the surrounding k Comprising arc section A 0 Is->
A functional relationship between, where k=1, 2, …, K, yields:
wherein,,
B sm ＝[Δ 1 ,Δ 2 ,…,Δ k ,…,Δ K ] T ，Δ k representing arc segment A k Coordinate increment between two endpoints, +.>
Representing arc segment A k Difference in InSAR deformation observations between the two endpoints.
7. The method for monitoring the surface deformation of the InSAR chronology based on the surface stress strain model as set forth in claim 6, wherein the step 3 further comprises:
taking a first scene SAR image of the M interferograms as a reference, and easily obtaining a conversion matrix B between deformation observed values corresponding to the InSAR interferograms and deformation of the InSAR sight line to time sequence ground surface trs Wherein B is trs The matrix size of the matrix is M x (T-1), each row corresponds to one interferogram, each column corresponds to one SAR image, the reference image is removed, in each row, the columns corresponding to the M interferogram main images are-1, the columns corresponding to the auxiliary images are +1, and other elements are 0;
based on (5) and matrix B trs The deformation phase difference delta L of all arc sections in the M interferograms can be established defo And arc section A 0 Time sequence earth surface deformation phase gradient
Functional relationship between:
wherein DeltaL defo ＝[(ΔL 1 ) T ,(ΔL 2 ) T ,…,(ΔL m ) T ,…,(ΔL M ) T ] T ,
Where t=2, 3, …, T,
representing the sign of the Cronecker product, < ->
And->
Respectively representing InSAR line-of-sight deformation phase at arc segment A during the acquisition of the t-th SAR image 0 Gradients in east-west, north-south, and vertical directions;
target arc segment A 0 The actual InSAR observed value comprises the earth surface deformation phase difference
Including topography residual dependent phase->
And noise phase, wherein the target arc segment A 0 Topography residual phase->
Expressed as:
wherein,,
wherein,,
representing arc segment A in the mth interferogram 0 Topography residual phase on ∈>
Vertical baseline representing mth interferogram, dz 0 Representing arc segment A 0 The terrain residual, λ, ρ, θ on the upper surface represent the satellite wavelength, satellite-to-target point distance, satellite incident angle at the target point, respectively.
8. The method for monitoring the surface deformation of the InSAR chronology based on the surface stress strain model as set forth in claim 7, wherein the step 3 further comprises:
assuming that the terrain residual correlation phase is spatially uncorrelated, arc segment A k The terrain residual phase can be regarded as noise, and the functional relation between the time sequence deformation phase gradient of the target arc section and the interference phase difference delta L of all arc sections in a certain surrounding range can be established based on the earth surface stress strain model by combining the formula (6) and the formula (7):
wherein DeltaL is the interference phase on the corresponding K arc segments in all M interferograms, which is the true view of the corresponding arc segment end pointsThe difference of the measured phases is obtained, and B= [ B ] defo ,B dz ],
Matrix B dz The size of (M×K) ×1, matrix B dz Line number and B of (2) defo Consistent with DeltaL, matrix B dz And arc segment A in delta L 0 The elements of the positions corresponding to the upper interference phase difference observation values are +.>
The other elements are all 0.
9. The method for monitoring the surface deformation of the InSAR chronology based on the surface stress strain model as set forth in claim 8, wherein the step 3 further includes:
weighting different arc segments according to the coherence of high coherence points, wherein the coherence is an InSAR observation value of c, and the corresponding variance sigma 2 Can be expressed as:
InSAR observed value variance of two endpoints of easily obtained component arc section
And->
Furthermore, the variance of the difference between the InSAR observations at two points can be obtained>
The method comprises the following steps:
ignoring the covariance of InSAR observations between different arc segments, obtaining a variance matrix D of an observation vector DeltaL based on the coherence values of the equation (9) and the equation (10) and all points ΔL Wherein D is ΔL The diagonal elements respectively correspond to the variances of the InSAR observed values delta L on the arc segments, the rest elements are 0, and the weight array P of the observed value vector delta L ΔL Can be expressed as:
P ΔL ＝(D ΔL ) -1 (II)。
10. the method for monitoring the surface deformation of the InSAR chronology based on the surface stress strain model according to claim 9, wherein the step 4 specifically comprises:
based on equation (8) and equation (11), an initial solution of the unknown parameter vector can be solved based on a weighted least squares criterion
The method comprises the following steps:
in the actual data, the difference between the observed values of adjacent points contains 2pi ambiguity, i.e. the observed value DeltaL contains coarse difference, in order to reduce the coarse difference to the unknown parameter vector
According to the statistical characteristics of the observed value, i.e. when the observed value correction is greater than a certain threshold v thr When the corresponding observed value is considered to be a gross error, where v thr Calculated according to the following formula:
wherein,,
observed value correction v 0 Calculated according to the following formula:
removing the correction v from DeltaL 0 >v rhr Is re-calculated according to equation (12)
I.e. is not influenced by the rough difference, the arc section A can be obtained 0 Difference dz between topography residual errors between upper starting and stopping points 0 And a time-sequential earth deformation gradient in east-west, north-south and vertical directions +.>
Obtaining arc section A 0 Two end points P i Sum point P j Time sequence earth surface deformation phase difference between
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010265777.6A CN111398959B (en) | 2020-04-07 | 2020-04-07 | InSAR time sequence earth surface deformation monitoring method based on earth surface stress strain model |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010265777.6A CN111398959B (en) | 2020-04-07 | 2020-04-07 | InSAR time sequence earth surface deformation monitoring method based on earth surface stress strain model |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN111398959A CN111398959A (en) | 2020-07-10 |
| CN111398959B true CN111398959B (en) | 2023-07-04 |
Family
ID=71431454
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202010265777.6A Active CN111398959B (en) | 2020-04-07 | 2020-04-07 | InSAR time sequence earth surface deformation monitoring method based on earth surface stress strain model |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN111398959B (en) |
Families Citing this family (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN112835043B (en) * | 2021-01-06 | 2023-03-21 | 中南大学 | Method for monitoring two-dimensional deformation in any direction |
| CN112797886B (en) * | 2021-01-27 | 2022-04-22 | 中南大学 | Winding phase oriented InSAR time sequence three-dimensional deformation monitoring method |
| CN113281744B (en) * | 2021-03-11 | 2023-03-21 | 中南大学 | Time sequence InSAR method based on hypothesis test and self-adaptive deformation model |
| CN114814839B (en) * | 2022-03-22 | 2024-05-03 | 武汉大学 | Wide area earth surface deformation detection method and system based on InSAR phase gradient stacking |
| CN115201823B (en) * | 2022-07-22 | 2023-08-04 | 电子科技大学 | Ground surface deformation monitoring method utilizing BDS-InSAR data fusion |
| CN116338690B (en) * | 2023-03-28 | 2024-01-16 | 中南林业科技大学 | Model constraint-free time sequence InSAR terrain residual error and earth surface deformation estimation method |
Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CA2652639A1 (en) * | 2008-02-06 | 2009-08-06 | Halliburton Energy Services, Inc. | Geodesy via gps and insar integration |
| WO2011154804A1 (en) * | 2010-06-07 | 2011-12-15 | Universitat Politècnica De Catalunya | Method for estimating the topography of the earth's surface in areas with plant cover |
| CN110888130A (en) * | 2019-10-30 | 2020-03-17 | 华东师范大学 | Coal mine area ground surface deformation monitoring method based on lifting rail time sequence InSAR |
Family Cites Families (8)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| GB0316858D0 (en) * | 2003-07-18 | 2003-08-20 | Univ Nottingham | A process for radar measurements of the position and movement of indivual targets and of areas of land and ice |
| US9417323B2 (en) * | 2012-11-07 | 2016-08-16 | Neva Ridge Technologies | SAR point cloud generation system |
| KR101605450B1 (en) * | 2014-08-04 | 2016-03-22 | 서울시립대학교산학협력단 | Method of stacking multi-temporal MAI interferogram and Apparatus Thereof |
| CN108594224B (en) * | 2018-03-30 | 2021-09-03 | 中国电力工程顾问集团中南电力设计院有限公司 | Three-dimensional time sequence deformation monitoring method fusing different platforms and orbit SAR data |
| CN108627833B (en) * | 2018-05-15 | 2021-08-24 | 电子科技大学 | Atmospheric phase compensation method based on GB-InSAR |
| CN109061641B (en) * | 2018-07-06 | 2020-01-17 | 中南大学 | InSAR time sequence earth surface deformation monitoring method based on sequential adjustment |
| CN109884635B (en) * | 2019-03-20 | 2020-08-07 | 中南大学 | Large-range high-precision InSAR deformation monitoring data processing method |
| CN110456345B (en) * | 2019-06-28 | 2020-11-10 | 深圳市水务规划设计院股份有限公司 | Building inclination monitoring method based on InSAR technology |
-
2020
- 2020-04-07 CN CN202010265777.6A patent/CN111398959B/en active Active
Patent Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CA2652639A1 (en) * | 2008-02-06 | 2009-08-06 | Halliburton Energy Services, Inc. | Geodesy via gps and insar integration |
| WO2011154804A1 (en) * | 2010-06-07 | 2011-12-15 | Universitat Politècnica De Catalunya | Method for estimating the topography of the earth's surface in areas with plant cover |
| CN110888130A (en) * | 2019-10-30 | 2020-03-17 | 华东师范大学 | Coal mine area ground surface deformation monitoring method based on lifting rail time sequence InSAR |
Non-Patent Citations (1)
| Title |
|---|
| INSAR矿区地表三维形变监测与预计研究进展;朱建军 等;测绘学报;第48卷(第2期);135-144 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN111398959A (en) | 2020-07-10 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN111398959B (en) | InSAR time sequence earth surface deformation monitoring method based on earth surface stress strain model | |
| CN106772342B (en) | Time sequence differential radar interference method suitable for large-gradient ground surface settlement monitoring | |
| CN104123464B (en) | Method for inversion of ground feature high elevation and number of land subsidence through high resolution InSAR timing sequence analysis | |
| CN110888130A (en) | Coal mine area ground surface deformation monitoring method based on lifting rail time sequence InSAR | |
| KR101490981B1 (en) | Method for Correction of Ionospheric Distortion of Synthetic Aperture Radar Interferogram and Apparatus Thereof | |
| CN111208512B (en) | Interferometric measurement method based on video synthetic aperture radar | |
| CN108663678B (en) | Multi-baseline InSAR phase unwrapping algorithm based on mixed integer optimization model | |
| CN103675790A (en) | Method for improving earth surface shape change monitoring precision of InSAR (Interferometric Synthetic Aperture Radar) technology based on high-precision DEM (Digital Elevation Model) | |
| CN111273293B (en) | InSAR residual motion error estimation method and device considering terrain fluctuation | |
| CN110456346B (en) | Power transmission tower inclination monitoring method based on InSAR technology | |
| CN112711021B (en) | Multi-resolution InSAR (interferometric synthetic Aperture Radar) interactive interference time sequence analysis method | |
| CN110673145A (en) | InSAR (interferometric synthetic Aperture Radar) surface deformation monitoring method and system based on discontinuous coherence | |
| CN111856459B (en) | Improved DEM maximum likelihood constraint multi-baseline InSAR phase unwrapping method | |
| CN112444188B (en) | Multi-view InSAR sea wall high-precision three-dimensional deformation measurement method | |
| CN116381680A (en) | Urban earth surface deformation monitoring method based on time sequence radar interferometry technology | |
| CN113238228B (en) | Three-dimensional earth surface deformation obtaining method, system and device based on level constraint | |
| Wang et al. | Microdeformation monitoring by permanent scatterer GB-SAR interferometry based on image subset series with short temporal baselines: The Geheyan Dam case study | |
| CN115453520A (en) | Surface deformation measurement method and device based on dual-frequency multi-polarization differential interference | |
| CN113341410B (en) | Large-range under-forest terrain estimation method, device, equipment and medium | |
| CN112268517B (en) | Method for monitoring deformation of power transmission tower equipment by using PSInSAR | |
| Rossi et al. | JERS-1 SAR image quality and interferometric potential | |
| Lombardi et al. | Accuracy of high resolution CSK interferometric Digital Elevation Models | |
| CN116363057B (en) | Surface deformation extraction method integrating PCA and time sequence InSAR | |
| CN116642440A (en) | Small baseline set interferometry method based on wavelet filtering | |
| CN117761716A (en) | InSAR landslide deformation time sequence monitoring method and storage medium fusing CR and GNSS |
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 |