CN116642440A - Small baseline set interferometry method based on wavelet filtering - Google Patents
Small baseline set interferometry method based on wavelet filtering Download PDFInfo
- Publication number
- CN116642440A CN116642440A CN202310494895.8A CN202310494895A CN116642440A CN 116642440 A CN116642440 A CN 116642440A CN 202310494895 A CN202310494895 A CN 202310494895A CN 116642440 A CN116642440 A CN 116642440A
- Authority
- CN
- China
- Prior art keywords
- phase
- interference
- wavelet
- time sequence
- surface deformation
- 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.)
- Granted
Links
- 238000000034 method Methods 0.000 title claims abstract description 63
- 238000001914 filtration Methods 0.000 title claims abstract description 30
- 238000005305 interferometry Methods 0.000 title claims abstract description 22
- 238000012545 processing Methods 0.000 claims abstract description 13
- 238000012937 correction Methods 0.000 claims abstract description 3
- 239000011159 matrix material Substances 0.000 claims description 14
- 238000004458 analytical method Methods 0.000 claims description 6
- 238000006073 displacement reaction Methods 0.000 claims description 6
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Chemical compound O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 claims description 6
- 230000001427 coherent effect Effects 0.000 claims description 5
- 230000001133 acceleration Effects 0.000 claims description 3
- 238000004364 calculation method Methods 0.000 claims description 3
- CLOMYZFHNHFSIQ-UHFFFAOYSA-N clonixin Chemical compound CC1=C(Cl)C=CC=C1NC1=NC=CC=C1C(O)=O CLOMYZFHNHFSIQ-UHFFFAOYSA-N 0.000 claims description 3
- 238000013461 design Methods 0.000 claims description 3
- 238000010586 diagram Methods 0.000 claims description 3
- 238000007689 inspection Methods 0.000 claims description 3
- 238000010276 construction Methods 0.000 claims description 2
- 239000005436 troposphere Substances 0.000 claims description 2
- 238000001514 detection method Methods 0.000 abstract description 10
- 230000006855 networking Effects 0.000 abstract 1
- 238000005516 engineering process Methods 0.000 description 5
- 230000000694 effects Effects 0.000 description 4
- 101100295776 Drosophila melanogaster onecut gene Proteins 0.000 description 3
- 238000007796 conventional method Methods 0.000 description 2
- 230000001419 dependent effect Effects 0.000 description 2
- 238000005315 distribution function Methods 0.000 description 2
- 238000012804 iterative process Methods 0.000 description 2
- 238000005259 measurement Methods 0.000 description 2
- NAWXUBYGYWOOIX-SFHVURJKSA-N (2s)-2-[[4-[2-(2,4-diaminoquinazolin-6-yl)ethyl]benzoyl]amino]-4-methylidenepentanedioic acid Chemical compound C1=CC2=NC(N)=NC(N)=C2C=C1CCC1=CC=C(C(=O)N[C@@H](CC(=C)C(O)=O)C(O)=O)C=C1 NAWXUBYGYWOOIX-SFHVURJKSA-N 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 230000009191 jumping Effects 0.000 description 1
- 230000007774 longterm Effects 0.000 description 1
- 238000000691 measurement method Methods 0.000 description 1
- 238000005457 optimization Methods 0.000 description 1
- 238000012950 reanalysis Methods 0.000 description 1
- 238000005070 sampling Methods 0.000 description 1
- 238000012795 verification Methods 0.000 description 1
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01B—MEASURING LENGTH, THICKNESS OR SIMILAR LINEAR DIMENSIONS; MEASURING ANGLES; MEASURING AREAS; MEASURING IRREGULARITIES OF SURFACES OR CONTOURS
- G01B15/00—Measuring arrangements characterised by the use of electromagnetic waves or particle radiation, e.g. by the use of microwaves, X-rays, gamma rays or electrons
- G01B15/06—Measuring arrangements characterised by the use of electromagnetic waves or particle radiation, e.g. by the use of microwaves, X-rays, gamma rays or electrons for measuring the deformation in a solid
-
- 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/885—Radar or analogous systems specially adapted for specific applications for ground probing
-
- 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
Landscapes
- Engineering & Computer Science (AREA)
- Remote Sensing (AREA)
- Physics & Mathematics (AREA)
- Radar, Positioning & Navigation (AREA)
- Electromagnetism (AREA)
- General Physics & Mathematics (AREA)
- Computer Networks & Wireless Communication (AREA)
- Radar Systems Or Details Thereof (AREA)
Abstract
The invention discloses a small baseline set interferometry method based on wavelet filtering, and relates to the technical field of earth surface deformation detection; the method comprises the following steps: acquiring SAR images for processing; adopting a soft threshold rule to identify and select stable points, calculating interference phase wavelet coefficient variance, and identifying stable pixels in an interference image by taking the interference phase wavelet coefficient variance as a standard; performing stable point networking and phase unwrapping; performing phase filtering and iterative correction; removing the DEM error of the interferogram, performing Legendre expansion on the phase of the DEM error, and estimating and eliminating the phase of the DEM error; generating a surface deformation time sequence; removing the atmospheric delay component from the surface deformation time sequence to obtain a final surface deformation time sequence; and estimating the linear deformation rate of the surface deformation by using a least square fitting method. The method improves the density and deformation detection precision of the target points identified by the traditional SBAS-InSAR method, and realizes a better small baseline set interferometry method.
Description
Technical Field
The invention belongs to the technical field of earth surface deformation detection, and particularly relates to a small baseline set interferometry method based on wavelet filtering.
Background
The small baseline set interferometry technique (Small Baseline Subset Interferometric Synthetic Aperture Radar, SBAS-InSAR) is a measurement technique that combines interference pairs of short, empty baselines and relatively small doppler center frequency differences, identifies phase information of surface stable high coherence points over long time intervals, and establishes a time-series phase model to estimate surface displacement, and is widely used for long-term surface deformation detection. How to accurately identify the stable high coherence points on the earth surface, calculate the coherence point displacement information, and effectively remove the influence of errors at the same time, thereby being a key factor for ensuring the deformation detection precision.
When the conventional SBAS-InSAR technology is used for identifying the surface stable points, a hard threshold method of 'one-cut' is adopted, so that data loss is caused, and the number and the density of the identified surface stable points are limited; when the minimum cost flow method is used for phase unwrapping, iteration and optimization are not performed, the unwrapping effect of a phase discontinuous region is not ideal, and unwrapping errors are easily caused; when the DEM error is estimated and the influence of atmospheric delay is eliminated, the mode of time domain high-pass filtering and space domain low-pass filtering is directly adopted, the loss of a low-frequency part and unstable components (such as jump caused by a hard threshold value) of data can be possibly caused, the real atmospheric condition is not considered, and the deformation measurement error of a complex area of the meteorological condition is larger. The above problems tend to increase the earth displacement estimation error, thereby reducing the deformation detection accuracy.
Aiming at the problem of low earth surface deformation detection precision caused by the defects in the traditional SBAS-InSAR technology, the invention provides a small baseline set interferometry method based on wavelet filtering.
Disclosure of Invention
The invention aims to provide a small baseline set interferometry method based on wavelet filtering, which aims to solve the problems of limited number and density of identified earth surface stable points, larger deformation measurement error of a complex meteorological condition area, non-ideal unwrapping effect of a phase discontinuous area and the like in the traditional SBAS-InSAR technology in the background technology.
In order to achieve the above purpose, the invention adopts the following technical scheme:
a small baseline set interferometry method based on wavelet filtering, comprising the steps of:
s1, acquiring SAR images for processing to obtain an interferogram set;
s2, identifying and selecting stable points by adopting a soft threshold rule; calculating the wavelet coefficient variance of the interference phase, and identifying stable pixels in the interference pattern by taking the wavelet coefficient variance as a standard;
s3, stable point network construction and phase unwrapping are carried out; constructing a delaunay triangle network by using the stable pixels identified in the S2, and unwrapping the sparse stable point network by using a minimum cost flow unwrapping algorithm;
s4, performing phase filtering and iterative correction; filtering and estimating high-frequency interference by using Legendre wavelet to the phase unwrapping result, correcting the unwrapped phase error and removing the spatial phase loss interference noise, and repeating S3 until the estimated high-frequency interference is smaller than a given threshold value;
s5, removing DEM errors of the interferograms processed in the step S4; performing Legend expansion on the DEM error phase, and estimating and eliminating the DEM error phase;
s6, generating a surface deformation time sequence; adopting a weighted iterative least square method, re-weighting the observed value as a function (such as double square) of the residual error obtained in the previous step, solving an equation by using a method for solving a generalized inverse matrix, and estimating a surface deformation time sequence;
s7, removing the atmospheric delay component from the surface deformation time sequence to obtain a final surface deformation time sequence;
s8, estimating the linear deformation rate of the surface deformation by using a least square fitting method.
Preferably, the step S1 of acquiring the SAR image for processing specifically includes:
acquiring N+1 SAR images with the same coverage range, and acquiring a time sequence t 0 ，t 1 ，…，t N N is a natural number;
and carrying out processing steps of image registration, small baseline group pairing, interference processing, orbit error estimation and terrain phase removal by using the SAR images with similar observation geometry, so as to obtain K interferogram sets with proper space baselines.
Preferably, the identifying of the stable pixel in the interferogram in S2 is specifically as follows:
s21, estimating phase noise in a complex domain by a wavelet analysis method;
s22, removing the original phase to reduce noise; when noise reduction is carried out on the interference phase, a soft threshold rule is adopted to carry out threshold control;
s23, establishing a Gaussian scattering model for the target scatterers distributed on the ground surface to obtain an interference phase density distribution function.
Preferably, the wavelet analysis method in S21 obtains a noise component based on a wavelet coefficient variance, specifically as follows:
the estimated interferometric phase noise phase components in the ith interferogram for pixel x are:
wherein the method comprises the steps ofAnd->Representing the real and imaginary parts of the noise phase, respectively;
the associated time domain variance is:
wherein the method comprises the steps ofc and s represent real and imaginary parts, respectively;
the relationship between the interference phase and complex phase observations satisfies:
wherein R and I are each a noise component comprising phase noiseAnd->Real and imaginary parts of the initial phase observations of (a);
standard deviation of interference phaseThe method comprises the following steps:
preferably, in S22, a soft threshold rule is used to perform threshold control, which specifically includes:
taking sigma as a reference variance of the interference phase, and according to a calculation formula of the variance
Wherein χ is 2 Is a chi-square probability density function;
σ is estimated by the expected variance of the deformation based on the InSAR gaze direction:
wherein DeltaR is radar side-looking skew difference and sigma l For its variance;
for a given significant level a of significance,the confidence interval is:
the pixels passing the inspection standard are stable pixels, otherwise are non-stable pixels.
Preferably, in the step S5, the DEM error phase is subjected to legendre unwrapping, which specifically includes:
for any pixel x (r, c) in the unwrapped interferogram, r and c are the coordinates of the point in the azimuth-oblique coordinate system, and the phase contribution of the DEM error is expressed as follows:
wherein lambda is radar wavelength, R represents the slant distance between the satellite and the target point, B ⊥ A vertical baseline distance representing two observations of a satellite, θ being the radar side view angle;
the legendre polynomial expansion form is:
wherein P is m Is an m-order Legend polynomial expansion;
the Legend spread of the DEM error phase is:
the first four components (m=0, 1,2, 3) in the above formula contain DEM errors of 85% or more and are removed in the interferogram.
Preferably, the generating of the surface deformation time sequence in S6 is specifically as follows:
definition ψ= [ ψ ] 1 ,K,ψ N ] T Vector composed of N unknown deformation phases, δφ= [ 1 δφ,K, q δφ] T Phase observations after K known unwraps;
the phase relation of the coherent point and the surface deformation in the interference diagram is expressed as:
wherein A is a design matrix,for observing errors +.>The estimated value of the displacement of the time sequence;
the least squares solution of the above formula is:wherein P is a weight matrix;
calculating new weights and updating coefficients by iterative processing until a preset threshold is met, e.g.Delta is a preset threshold.
Preferably, in S7, the removing of the atmospheric delay component is performed on the time series of the surface deformation, specifically as follows:
estimating the atmospheric delay phase of the radar line-of-sight direction using ERA-5, LOS single-path troposphere delay at elevation zIs the surface elevation z and the reference plane elevation z ref The integral of the refractive index of air in between, expressed as:
wherein θ is the angle of incidence;
R d ＝287.05J·kg -1 ·K -1 ，R v ＝461.495J·kg -1 ·K -1 the dry air and water vapor specific gas constants, respectively;
g m is a weighted average of z and z ref Gravitational acceleration between;
p is the partial pressure of dry air in Pa;
e is the partial pressure of water vapor in Pa;
t is the temperature in K;
constant k 1 ＝0.776K·Pa -1 ，k 2 ＝0.716K·Pa -1 ，k 3 ＝3.75·10 3 K 2 ·Pa -1 ；
z ref Setting to 10000m;
the first term in the above equation corresponds to the dry air component of the delay path, and the second term relates to air humidity; and estimating and eliminating the atmospheric delay component in the time sequence through the formula to obtain the final surface deformation time sequence.
Compared with the prior art, the invention has the beneficial effects that:
(1) The invention provides a small baseline set interferometry method based on wavelet filtering, which improves the density and deformation detection precision of target points identified by the traditional SBAS-InSAR method and realizes a better small baseline set interferometry method.
(2) The method is controlled by adopting a scale-related soft threshold rule when the earth surface stable point is identified, so that the data loss caused by the traditional 'one-cut' hard threshold method is avoided, and the density of target points identified by the traditional SBAS-InSAR method is improved.
(3) The method utilizes Legend wavelet to filter and estimate the phase unwrapping result to estimate high-frequency interference, and iteratively corrects the unwrapped phase error and removes spatial phase loss interference noise, thereby reducing the unwrapping error and improving the unwrapping precision.
(4) The Legend expansion type estimation and DEM error removal are adopted in the method, so that the accuracy of ground deformation detection is improved.
(5) According to the method, the ERA-5 atmospheric model is used for replacing a time domain high-pass filtering mode and a space domain low-pass filtering mode to estimate and remove the influence of atmospheric delay, so that the accuracy of ground deformation detection is improved.
Drawings
FIG. 1 is a flow chart of a small baseline set interferometry method based on wavelet filtering in the present invention;
FIG. 2 is a graph showing the comparison of the identification of stable points according to the present invention and the identification of stable points according to the conventional method in example 2;
FIG. 3 is a graph showing the unwrapping result of the time series in example 2;
FIG. 4 is a time series estimation chart of the surface deformation in example 2;
fig. 5 is a graph of verification of the accuracy of the results in example 2.
Detailed Description
The following description of the embodiments of the present invention will be made clearly and completely with reference to the accompanying drawings, in which it is apparent that the embodiments described are only some embodiments of the present invention, but not all embodiments. All other embodiments, which can be made by those skilled in the art based on the embodiments of the invention without making any inventive effort, are intended to be within the scope of the invention.
Example 1:
a small baseline set interferometry method based on wavelet filtering, comprising the steps of:
step one: assuming that there are n+1 SAR images with the same coverage, acquiring a time sequence t 0 ，t 1 ，…，t N N is a natural number. The SAR images have similar observation geometry, and K sets of interferograms with proper spatial baselines are obtained through the steps of image registration, small baseline pairing, interference processing, orbit error estimation, terrain phase removal and the like.
Step two: and calculating the interference phase wavelet coefficient variance, and identifying the stable pixel by taking the interference phase wavelet coefficient variance as a standard. Firstly, estimating phase noise in a complex domain by a wavelet analysis method, and removing (namely noise reduction) from an original phase; then, a Gaussian scattering model (Gaussian scattering model) is built for the target scatterers of the surface distribution to obtain an interference phase density distribution function.
In order to avoid the hard threshold function 'one-cut' effect from jumping or dithering in the wavelet domain when the interference phase is denoised, a Scale dependent soft threshold rule (Scale-dependent soft thresholding scheme) is used for threshold control. The method acquires noise components based on wavelet coefficient variances, so that prior information is not needed, and data loss is avoided.
Assume that the estimated interferometric phase noise phase component in the ith interferogram for pixel x isWherein->And->Representing the real and imaginary parts of the noise phase, respectively; its associated time domain variance is +.>Wherein->c and s represent real and imaginary parts, respectively. The relationship between the interference phase and complex phase observations satisfies: />Wherein R and I are each +.>And->Real and imaginary parts of the initial phase observations of (a).
Standard deviation of interference phaseThe method comprises the following steps:
assuming σ is the reference variance of the interference phase, according to the calculation formula of the variance, thenWherein χ is 2 As a chi-square probability density function. For a given significant level α, ->The confidence interval is:
the pixels passing the inspection standard are stable pixels, otherwise are non-stable pixels. It should be noted that σ may be estimated based on the expected variance of the InSAR line-of-sight deformation:
wherein DeltaR is radar side-looking skew difference and sigma l Is its variance.
Studies have shown that for SAR data, the phase noise can be expressed as a function of the number of views (as shown in equation (4), while the interference coherence of SAR images is independent of the number of views (as shown in equation (5)). Therefore, the method based on the time sequence interference phase noise statistical characteristics can identify stable pixels under the full resolution, and is superior to the traditional small baseline set technology.
Wherein N is r And N az Respectively representing the distance direction and azimuth direction vision number, wherein SNR is radar image signal-to-noise ratio, mu dB Representing phase noise. The coherence of the interferogram can be expressed as:
wherein E x]Representing the expectation of x, |u 1 I and u 1 The i indicates the amplitudes of the two SAR images forming the interference, respectively, and the γ indicates the interference coherence thereof.
Step three: and (3) constructing a delaunay triangle network (Delaunay triangulation network) by utilizing the stable pixels identified in the step two, and unwrapping the sparse stable point network by adopting a minimum cost flow unwrapping algorithm.
Step four: and (3) filtering the phase unwrapping result by using Legend wavelet to estimate high-frequency interference, correcting the phase of the corresponding target point, and repeating the third step until the estimated high-frequency interference is smaller than a given threshold value. After this step, the phase noise mainly contains DEM errors and atmospheric delay components, and is removed in a subsequent step.
Step five: and carrying out Legendre expansion on the DEM error phase, and estimating and eliminating the DEM error phase.
For any pixel x (r, c) in the disentangled interferogram, r and c are the coordinates of the point in an azimuth-oblique coordinate system, and Δh (x) is used for representing the error between the DEM and the real terrain, then the phase contribution of the DEM error can be expressed as:
wherein lambda is radar wavelength, R represents the slant distance between the satellite and the target point, B ⊥ The vertical baseline distance of two observations of the satellite is represented, θ being the radar side view angle.
The legendre polynomial expansion form is:
wherein P is m Is an expansion of an m-order Legend polynomial. The Legend spread of the DEM error phase is:
the first four components (m=0, 1,2, 3) in the above formula contain DEM errors of 85% or more and are removed in the interferogram. After this step, the phase noise mainly contains an atmospheric delay component, which is removed in a subsequent step.
Step six: and (3) adopting a weighted iterative least squares (Iteratively Reweighted Least Squares, IRLS) method, re-weighting the observed value as a function (such as double square) of the residual error obtained in the previous step, and solving an equation by using a method for solving a generalized inverse matrix to estimate the earth surface deformation time sequence.
The radar satellite images at the time a and the time b twice respectively, and in the phase interference diagram generated by interference processing, the relation between the interference phase value of the coherent point and the earth surface deformation phase can be expressed as follows:
wherein ψ= [ ψ ] 1 ,K,ψ N ] T Vector composed of N unknown deformation phases, δφ= [ 1 δφ,K, q δφ] T Is K known unwrapped phase observations. The matrix form of formula (10) is expressed as: the relationship may be expressed as:
wherein A is a design matrix,for observing errors +.>Is a time series displacement amount estimation value. In order to solve the equation (11) and also to reduce the effect of outliers (e.g., phase unwrapping errors), a weighted iterative least squares (Iteratively Reweighted Least Squares, IRLS) method is used to re-weight the observations as a function of the residual from the previous step (e.g., double square). To use the IRLS method, consider the least squares solution of formula (11) as:
wherein P is a weight matrix. Since the observations are mutually independent subsets that overlap in time, the solution of the equation must be performed using a method that solves the generalized inverse matrix when performing the inversion operation. This allows for a more realistic assessment of the observed noise and results in a more robust result. The iterative process of IRLS starts with estimating new weights of observations, and the iterative process is shown in formula (13):
wherein t is a harmonic factor, and the solving method is Holland and Wetsch, df is the degree of freedom (the number of interferograms minus the number of SAR images plus 1, i.e. k-N + 1),the process of calculating the new weights and updating the coefficients is iterated until a preset threshold is met (e.g +.>δ is a preset threshold). After n iterations, the completion variance covariance matrix of the surface deformation time sequence can be obtained at the same time>Wherein (1)>And Q n The final variance factor and the weight matrix, respectively.
Step seven: the atmospheric delay phase of the radar line of sight (LOS) was estimated using ERA-5 (European Centre for Medium-Range Weather Forecasts Reanalysis v, fifth generation ECMWF atmospheric analysis dataset). LOS single path tropospheric delay at elevation zIs the surface elevation z and the reference plane elevation z ref The integral of the refractive index of air in between, expressed as:
wherein θ is the incident angle, R d ＝287.05J·kg -1 ·K -1 ，R v ＝461.495J·kg -1 ·K -1 The specific gas constants of dry air and water vapor, g m Is a weighted average of z and z ref The gravitational acceleration between the two is that P is the partial pressure of dry air in Pa, e is the partial pressure of water vapor in Pa, and T is the temperature in K. Constant k 1 ＝0.776K·Pa -1 ，k 2 ＝0.716K·Pa -1 ，k 3 ＝3.75·10 3 K 2 ·Pa -1 ，z ref Typically set to 10000m. (14) Where the first term corresponds to the dry air component of the delay path and the second term relates to the air humidity. And (3) estimating and eliminating the atmospheric delay component in the time sequence through the formula (14) to obtain the final surface deformation time sequence.
Step eight: and finally, estimating the linear deformation rate by using a least square method.
Example 2:
an embodiment using the wavelet filtering-based small baseline set interferometry method of example 1 was performed.
Step one: 120 Sentinel-1 (C band) SLC products covering the Beijing area are collected, and the steps of small baseline pairing, registration, interference processing, orbit error estimation, terrain phase removal and the like are carried out, so that 235 interferogram sets with proper space baselines are obtained.
Step two: and (3) adopting a (2) type soft threshold rule to identify and select stable points. By contrast, the density of stable coherent points obtained by the treatment of the traditional hard threshold method is 299/km 2 The coherence point density identified by the soft threshold method provided by the invention is improved to 532/km 2 . The density and number of stable coherent points obtained by recognition are obviously improved, see fig. 2. The left picture in fig. 2 is the stable point identified by the present invention, and the right picture is the stable point identified by the conventional method.
Step three: and (3) constructing a Delaunay triangle network, and disentangling the sparse stable point network in all the interferograms by adopting a minimum cost flow disentanglement algorithm, wherein the disentanglement result is shown in fig. 3.
Step four: and (3) filtering and estimating high-frequency interference by utilizing Legend wavelets to the phase unwrapping result, setting a high-frequency interference part threshold value as 12 according to the sampling period of the data, repeating the step three for iteration, correcting the unwrapped phase error, removing the spatial phase loss interference noise and reducing the unwrapping error.
Step five: and (3) carrying out Legendre expansion on the phase unwrapping result in the step four according to the formula (9), calculating the first four terms, estimating the DEM error, and removing the DEM error from the phase.
Step six: and (3) iterating according to formulas (10) - (13), re-weighting the observed value as a function (such as double square) of the residual error obtained in the previous step, solving an equation by using a method for solving a generalized inverse matrix, and estimating the earth surface deformation time sequence. See fig. 4.
Step seven: the atmospheric delay error is estimated pixel by pixel according to equation (14) and removed from the time series.
Step eight: and estimating the linear deformation rate by using a least square fitting method. Referring to fig. 5, the correlation between the two data (R 2 ) Reaching 0.82, the error is better than 6mm/year.
The foregoing is only for aiding in understanding the method and the core of the invention, but the scope of the invention is not limited thereto, and it should be understood that the technical scheme and the inventive concept according to the invention are equivalent or changed within the scope of the invention by those skilled in the art. In view of the foregoing, this description should not be construed as limiting the invention.
Claims (8)
1. A small baseline set interferometry method based on wavelet filtering, comprising the steps of:
s1, acquiring SAR images for processing to obtain an interferogram set;
s2, identifying and selecting stable points by adopting a soft threshold rule; calculating the wavelet coefficient variance of the interference phase, and identifying stable pixels in the interference pattern by taking the wavelet coefficient variance as a standard;
s3, stable point network construction and phase unwrapping are carried out; constructing a delaunay triangle network by using the stable pixels identified in the S2, and unwrapping the sparse stable point network by using a minimum cost flow unwrapping algorithm;
s4, performing phase filtering and iterative correction; filtering and estimating high-frequency interference by using Legendre wavelet to the phase unwrapping result, correcting the unwrapped phase error and removing the spatial phase loss interference noise, and repeating S3 until the estimated high-frequency interference is smaller than a given threshold value;
s5, removing DEM errors of the interferograms processed in the step S4; performing Legend expansion on the DEM error phase, and estimating and eliminating the DEM error phase;
s6, generating a surface deformation time sequence; adopting a weighted iterative least square method, re-weighting the observed value as a function of the residual error obtained in the previous step, and solving an equation by using a method for solving a generalized inverse matrix to estimate the earth surface deformation time sequence;
s7, removing the atmospheric delay component from the surface deformation time sequence to obtain a final surface deformation time sequence;
s8, estimating the linear deformation rate of the surface deformation by using a least square fitting method.
2. The wavelet filtering-based small baseline set interferometry method according to claim 1, wherein the step of acquiring the SAR image in S1 is performed as follows:
acquiring N+1 SAR images with the same coverage range, and acquiring a time sequence t 0 ，t 1 ，…，t N N is a natural number;
and carrying out processing steps of image registration, small baseline group pairing, interference processing, orbit error estimation and terrain phase removal by using the SAR images with similar observation geometry, so as to obtain K interferogram sets with proper space baselines.
3. The small baseline set interferometry method based on wavelet filtering according to claim 1 or 2, wherein the identifying stable pixels in the interferogram in S2 is specifically as follows:
s21, estimating phase noise in a complex domain by a wavelet analysis method;
s22, removing the original phase to reduce noise; when noise reduction is performed on the interference phase, a soft threshold rule is adopted for threshold control.
4. A wavelet-filtering-based small baseline set interferometry method according to claim 3, wherein the wavelet analysis method in S21 obtains noise components based on wavelet coefficient variance, specifically as follows:
the estimated interferometric phase noise phase components in the ith interferogram for pixel x are:
wherein the method comprises the steps ofAnd->Representing the real and imaginary parts of the noise phase, respectively;
the associated time domain variance is:
wherein the method comprises the steps ofw=c, s, c and s represent real and imaginary parts, respectively;
the relationship between the interference phase and complex phase observations satisfies:
wherein R and I are each a noise component comprising phase noiseAnd->Real and imaginary parts of the initial phase observations of (a);
standard deviation of interference phaseThe method comprises the following steps:
5. the wavelet filtering-based small baseline set interferometry method according to claim 4, wherein the threshold control is performed by adopting a soft threshold rule in S22, specifically as follows:
taking sigma as a reference variance of the interference phase, and according to a calculation formula of the variance
Wherein χ is 2 Is a chi-square probability density function;
σ is estimated by the expected variance of the deformation based on the InSAR gaze direction:
wherein DeltaR is radar side-looking skew difference and sigma l For its variance;
for a given significant level a of significance,the confidence interval is:
the pixels passing the inspection standard are stable pixels, otherwise are non-stable pixels.
6. The wavelet filtering-based small baseline set interferometry method according to claim 1, wherein the step S5 of performing legendre expansion on DEM error phases is as follows:
for any pixel x (r, c) in the unwrapped interferogram, r and c are the coordinates of the point in the azimuth-oblique coordinate system, and the phase contribution of the DEM error is expressed as follows:
wherein lambda is radar wavelength, R represents the slant distance between the satellite and the target point, B ⊥ A vertical baseline distance representing two observations of a satellite, θ being the radar side view angle;
the legendre polynomial expansion form is:
wherein P is m Is an m-order Legend polynomial expansion;
the Legend spread of the DEM error phase is:
the first four components (m=0, 1,2, 3) in the above formula contain DEM errors of 85% or more and are removed in the interferogram.
7. The wavelet filtering-based small baseline set interferometry method according to claim 1, wherein the generating of the surface deformation time sequence in S6 is specifically as follows:
definition ψ= [ ψ ] 1 ,K,ψ N ] T Vector composed of N unknown deformation phases, δφ= [ 1 δφ,K, q δφ] T Phase observations after K known unwraps;
the phase relation of the coherent point and the surface deformation in the interference diagram is expressed as:
wherein A is a design matrix,as an observation error, ψ) is a time series displacement amount estimation value;
the least squares solution of the above formula is:wherein P is a weight matrix;
calculating new weights and updating coefficients by iterative processing until a preset threshold is met, e.g.Delta is a preset threshold.
8. The wavelet filtering-based small baseline set interferometry method according to claim 1, wherein the step S7 of performing the atmospheric delay component removal on the surface deformation time sequence is specifically as follows:
estimating the atmospheric delay phase of the radar line-of-sight direction using ERA-5, LOS single-path troposphere delay at elevation zIs the surface elevation z and the reference plane elevation z ref The integral of the refractive index of air in between, expressed as:
wherein θ is the angle of incidence;
R d ＝287.05J·kg -1 ·K -1 ，R v ＝461.495J·kg -1 ·K -1 the dry air and water vapor specific gas constants, respectively;
g m is a weighted average of z and z ref Gravitational acceleration between;
p is the partial pressure of dry air in Pa;
e is the partial pressure of water vapor in Pa;
t is the temperature in K;
constant k 1 ＝0.776K·Pa -1 ，k 2 ＝0.716K·Pa -1 ，k 3 ＝3.75·10 3 K 2 ·Pa -1 ；
z ref Setting to 10000m;
the first term in the above equation corresponds to the dry air component of the delay path, and the second term relates to air humidity; and estimating and eliminating the atmospheric delay component in the time sequence through the formula to obtain the final surface deformation time sequence.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310494895.8A CN116642440B (en) | 2023-05-05 | 2023-05-05 | Small baseline set interferometry method based on wavelet filtering |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310494895.8A CN116642440B (en) | 2023-05-05 | 2023-05-05 | Small baseline set interferometry method based on wavelet filtering |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN116642440A true CN116642440A (en) | 2023-08-25 |
| CN116642440B CN116642440B (en) | 2024-07-23 |
Family
ID=87639086
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202310494895.8A Active CN116642440B (en) | 2023-05-05 | 2023-05-05 | Small baseline set interferometry method based on wavelet filtering |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN116642440B (en) |
Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP2004191053A (en) * | 2002-12-06 | 2004-07-08 | Mitsubishi Electric Corp | Synthetic aperture radar device and numerical altitude model creation method |
| US20060034545A1 (en) * | 2001-03-08 | 2006-02-16 | Universite Joseph Fourier | Quantitative analysis, visualization and movement correction in dynamic processes |
| CN107132539A (en) * | 2017-05-03 | 2017-09-05 | 中国地质科学院探矿工艺研究所 | Landslide early-stage identification method of time sequence InSAR (interferometric synthetic Aperture Radar) based on small baseline set |
| CN110187346A (en) * | 2019-06-21 | 2019-08-30 | 湖南科技大学 | A kind of ground SAR Detection of Gross Errors method under complex working condition |
| CN115906951A (en) * | 2022-11-24 | 2023-04-04 | 重庆理工大学 | Learnable activation function method based on orthogonal Legendre wavelets |
-
2023
- 2023-05-05 CN CN202310494895.8A patent/CN116642440B/en active Active
Patent Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20060034545A1 (en) * | 2001-03-08 | 2006-02-16 | Universite Joseph Fourier | Quantitative analysis, visualization and movement correction in dynamic processes |
| JP2004191053A (en) * | 2002-12-06 | 2004-07-08 | Mitsubishi Electric Corp | Synthetic aperture radar device and numerical altitude model creation method |
| CN107132539A (en) * | 2017-05-03 | 2017-09-05 | 中国地质科学院探矿工艺研究所 | Landslide early-stage identification method of time sequence InSAR (interferometric synthetic Aperture Radar) based on small baseline set |
| CN110187346A (en) * | 2019-06-21 | 2019-08-30 | 湖南科技大学 | A kind of ground SAR Detection of Gross Errors method under complex working condition |
| CN115906951A (en) * | 2022-11-24 | 2023-04-04 | 重庆理工大学 | Learnable activation function method based on orthogonal Legendre wavelets |
Non-Patent Citations (2)
| Title |
|---|
| MINGLIANG GAO等: "Regional Land Subsidence Analysis in Eastern Beijing Plain by InSAR Time Series and Wavelet Transforms", REMOTE SENSING, 26 February 2018 (2018-02-26) * |
| 朱锋等: "基于ＩｎＳＡＲ和小波变换的不均匀沉降段识别 ——以京津高铁北京段为例", 地理与地理信息科学, 31 January 2014 (2014-01-31) * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN116642440B (en) | 2024-07-23 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN111273293B (en) | InSAR residual motion error estimation method and device considering terrain fluctuation | |
| Shirzaei | A wavelet-based multitemporal DInSAR algorithm for monitoring ground surface motion | |
| CN106772342B (en) | Time sequence differential radar interference method suitable for large-gradient ground surface settlement monitoring | |
| Sica et al. | $-Net: Deep Residual Learning for InSAR Parameters Estimation | |
| Loew et al. | Generation of geometrically and radiometrically terrain corrected SAR image products | |
| Shimada | Verification processor for SAR calibration and interferometry | |
| CN108007476B (en) | Interference calibration method and system for space-based interference imaging radar altimeter | |
| CN111398959B (en) | InSAR time sequence earth surface deformation monitoring method based on earth surface stress strain model | |
| Liao et al. | Improved topographic mapping through high-resolution SAR interferometry with atmospheric effect removal | |
| CN111856459B (en) | Improved DEM maximum likelihood constraint multi-baseline InSAR phase unwrapping method | |
| CN110703220B (en) | Multi-baseline PolInSAR vegetation parameter inversion method considering time decoherence factors | |
| CN112711021B (en) | Multi-resolution InSAR (interferometric synthetic Aperture Radar) interactive interference time sequence analysis method | |
| CN111239736B (en) | Single-baseline-based surface elevation correction method, device, equipment and storage medium | |
| CN114114181B (en) | Satellite-borne SAR interference baseline correction method based on orbit error phase basis | |
| Pinheiro et al. | Generation of highly accurate DEMs over flat areas by means of dual-frequency and dual-baseline airborne SAR interferometry | |
| CN110018476A (en) | Time difference baseline set timing interference SAR processing method | |
| CN117406221A (en) | InSAR deformation monitoring method for high-resolution DEM of corner reflector | |
| CN113341410B (en) | Large-range under-forest terrain estimation method, device, equipment and medium | |
| CN116642440B (en) | Small baseline set interferometry method based on wavelet filtering | |
| CN112946647A (en) | Atmospheric error correction InSAR interferogram stacking geological disaster general investigation method and device | |
| Monaldo | Improvement in the estimation of dominate wavenumber and direction from spaceborne SAR image spectra when corrected for ocean surface movement | |
| CN115453533A (en) | Time sequence InSAR distributed target phase optimization method | |
| Zhao et al. | A combined multi-interferogram algorithm for high resolution DEM reconstruction over deformed regions with TerraSAR-X data | |
| Sabater et al. | Physical analysis of atmospheric delay signal observed in stacked radar interferometric data | |
| CN118376993B (en) | Time sequence interference synthetic aperture radar terrain residual error estimation method |
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 |