CN115856889A - InSAR time sequence deformation monitoring method capable of automatically correcting errors - Google Patents
InSAR time sequence deformation monitoring method capable of automatically correcting errors Download PDFInfo
- Publication number
- CN115856889A CN115856889A CN202310054817.6A CN202310054817A CN115856889A CN 115856889 A CN115856889 A CN 115856889A CN 202310054817 A CN202310054817 A CN 202310054817A CN 115856889 A CN115856889 A CN 115856889A
- Authority
- CN
- China
- Prior art keywords
- phase
- interferogram
- threshold
- unwrapping
- 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
Images
Abstract
The invention provides an InSAR time sequence deformation monitoring method with automatic error correction, which relates to the technical field of surface deformation monitoring and comprises the following steps: s1, acquiring a target area, selecting a time interval and selecting an SAR image; s2, setting a space-time baseline threshold, constructing a differential interferogram set of a short baseline set, and performing phase unwrapping and geocoding to obtain an unwrapped interferogram set; s3, identifying stable points in the target area, and performing combined correction of orbit errors and atmospheric phases on the unwrapped interferogram set based on the stable points to obtain a corrected interferogram set; s4, screening the correction interference pattern set according to screening conditions to obtain a high-quality interference pattern set; and S5, resolving and processing the high-quality interference pattern set based on a deformation algorithm to obtain ground deformation information of the target area. The method can automatically identify stable points in a research area, is used for performing combined correction processing on the orbit error and the atmospheric phase by interference, and intelligently performs quality screening and deformation calculation on the interference pattern.
Description
Technical Field
The invention relates to the technical field of terrain monitoring, in particular to an InSAR time sequence deformation monitoring method capable of automatically correcting errors.
Background
The interferometric radar refers to a synthetic aperture radar (InSAR) adopting an interferometric technique, is a newly developed space-to-ground observation technique, and is a product of combining a traditional SAR remote sensing technique and a radio astronomical interference technique. The method comprises the steps of transmitting microwaves to a target area by utilizing a radar, receiving echoes reflected by the target, obtaining an SAR complex image pair imaged in the same target area, obtaining an interference diagram by conjugate multiplication of the SAR complex image pair if a coherence condition exists between the complex image pair, and obtaining a path difference of the microwaves in two imaging processes according to a phase value of the interference diagram, so that the landform, the landform and the surface micro-change of the target area are calculated, and the method can be used for digital elevation model establishment, crust deformation detection and the like. In recent years, with the continuous emission of SAR satellites and the improvement of InSAR data processing level, PSI and SBAS technologies are taken as representatives, and the time sequence InSAR technology can reach the deformation monitoring precision of millimeter level.
However, in the time sequence InSAR engineering application process, since the differential interferogram is doped with many error items such as orbit error, unwrapping error, atmospheric delay and the like, correction of these errors or elimination of poor interference pairs requires a great deal of subjective experience. For example, subjective selection of stable points is involved in track error and laminar atmospheric phase correction, and errors in differential interference phases are difficult to correct accurately when stable point selection is inappropriate.
Disclosure of Invention
Aiming at the technical problem, the invention provides an InSAR time sequence deformation monitoring method capable of automatically correcting errors, which can automatically identify stable points in a research area, is used for joint correction processing of interference on orbit errors and atmospheric phases, and intelligently screens the quality of an interferogram and solves the deformation of the interferogram.
In order to achieve the technical purpose, the invention provides an InSAR time sequence deformation monitoring method with automatic error correction, which comprises the following steps:
s1, acquiring a target area, selecting a time interval, and selecting an SAR image of the time interval in the target area;
s2, setting a space-time baseline threshold, constructing a differential interferogram set of a short baseline set based on the space-time baseline threshold, and performing phase unwrapping and geocoding to obtain an unwrapped interferogram set;
s3, identifying stable points in the target area, and performing combined correction of orbit errors and atmospheric phases on the unwrapping interference pattern set based on the stable points to obtain a corrected interference pattern set;
s4, screening the correction interference pattern set according to screening conditions to obtain a high-quality interference pattern set;
and S5, resolving and processing the high-quality interference pattern set based on a deformation algorithm to obtain ground deformation information of the target area.
In an embodiment of the present invention, the step S3 includes:
s31, acquiring the average coherence of the unwrapping interferogram, and setting a first coherence threshold;
s32, calculating the phase standard deviation of the unwrapping interferogram, and setting a first phase standard deviation threshold value;
s33, performing quality screening on the unwrapping interferogram set according to the first coherence threshold and the first phase standard deviation threshold, and if the average coherence of the unwrapping interferogram is smaller than the first coherence threshold or the phase standard deviation of the unwrapping interferogram is larger than the first phase standard deviation threshold, screening the unwrapping interferogram to obtain the unwrapping interferogram set after quality screening;
s34, stacking time sequence InSAR processing is carried out on the unwrapped interferogram set after the quality screening, and the phase change rate is obtained according to a stacking time sequence algorithm;
s35, setting a threshold interval of a stable point, obtaining deformation quantities of all points in the target area according to the phase change rate, and if the deformation quantities of a certain point are distributed in the threshold interval of the stable point, the point is the stable point;
s36, calculating the joint phase of the orbit error and the atmosphere in the layer according to the position information of the stable point, and performing phase correction on the unwrapped interferogram set by taking the joint phase as an error signal to obtain a corrected interferogram set.
In an embodiment of the present invention, the step S36 includes:
establishing a functional relation among the unwrapping phase, the elevation and the position according to the position information of the stable point, wherein the formula is as follows:
wherein x and y are respectively the row number and the column of the pixel in the geographic coordinate systemThe number of the mobile station is,
is the unwrapped phase at point (x, y), h (x, y) is the elevation of point (x, y), and->
Is the random phase error of point (x, y), a i =(a 0 ,a 1 ,...a 6 ) Is the correlation coefficient;
and solving the joint phase of the orbit error and the layered atmosphere by a least square method, and removing the joint phase from the unwrapping phase of the above formula to obtain a third interferogram.
In an embodiment of the present invention, in step S34, the stacking timing algorithm is: setting a minimum stack interferogram number threshold
Stacking time sequence deformation calculation is carried out on the unwrapping interferogram set subjected to quality screening according to an intermittent mode, and when the number of the stacked interferograms is greater than or equal to>
Then, the stack yields the phase change rate:
wherein the content of the first and second substances,
representing the rate of change of phase, N being the number of unwrapped interferograms>
For the time base of the jth unwrapped interferogram, the->
Is the phase value of the jth unwrapped interferogram.
In one embodiment of the present invention, in step S35, the stable point threshold interval is near the median of the deformation of all the points in the target region.
In an embodiment of the present invention, step S4 includes:
s41, acquiring average coherence of the corrected interferogram, setting a second coherence threshold, calculating a phase standard deviation of the corrected interferogram, and setting a second phase standard deviation threshold;
s42, quality screening is carried out on the correction interferogram set according to a second coherence threshold and a second phase standard deviation threshold, if the average coherence of the correction interferogram is smaller than the second coherence threshold or the phase standard deviation of the correction interferogram is larger than the second phase standard deviation threshold, the correction interferogram is screened out, and the correction interferogram set after quality screening is obtained;
s43, screening the corrected interference pattern set after quality screening by adopting a triangular phase closed-loop method, forming the corrected interference pattern set after quality screening into a triangular phase closed-loop network formed by one or more triangles, and generating three unwrapping interference differential phases by each triangular phase closed-loop;
s44, setting a closed loop phase threshold, and calculating the root mean square of the corresponding closed loop phase according to the three unwrapping interference differential phases of each group of triangular phase closed loops:
when the root mean square of the closed loop phase is less than the closed loop phase threshold, the three unwrapped interfering differential phases are retained,
when the root mean square of the closed loop phase is larger than the closed loop phase threshold, eliminating the three unwrapping interference differential phases;
s45, repeating the step S44 to obtain a high-quality interferogram set.
In a specific embodiment of the present invention, the deformation algorithm is a weighted least square deformation calculation algorithm with an additional smooth constraint condition, and the step S5 includes:
s51, partitioning a target area, acquiring high coherence points in each area block, and extracting unwrapping phases and coherence information of each high coherence point;
s52, establishing a functional relation among the unwrapping phase, the deformation rate and the elevation residual error in each area block, solving the deformation rate of each high-coherence point according to the coherence information, solving according to the functional relation to obtain the elevation residual error, and splicing the results of each area block to obtain the accumulated deformation at each moment;
and S53, carrying out turbulence atmospheric correction and noise removal processing on the accumulated deformation at each moment according to a filtering method to obtain the earth surface deformation information of the target area.
In an embodiment of the present invention, step S52 includes:
for a single high coherence point to be solved, the observation equation is:
wherein the content of the first and second substances,
is a time interval, is->
Is changed at the change rate->
Is an elevation residual->
For unwrapping phase, M is the number of high-quality interferograms, N +1 is the number of SAR images, θ is the angle of incidence, B is the vertical baseline of the high-quality interferograms, R is the distance between the SAR image radar antenna and the ground target, π is a constant, λ is the wavelength;
introducing virtual observations
Adding a smooth constraint condition with equal acceleration for solving as initial data, wherein the smooth constraint condition is as follows:
in the formula
To indicate the ground orderMarking the deformation rate at the kth moment, introducing a smoothing parameter alpha to control the constraint strength, and combining an observation equation and a smoothing constraint condition to obtain:
wherein
For smooth constraint coefficients>
Is a rate factor, is based on the rate of the signal>
Is an elevation residual term coefficient>
Is a smooth constraint matrix;
solving the rate item and the elevation residual error, and splicing the result of each area block to obtain the accumulated deformation at each moment.
In an embodiment of the present invention, step S53 includes:
s531 processing the accumulated deformation of each moment according to spatial low filtering to remove a noise phase in the accumulated deformation;
s532, removing turbulence atmospheric components in the accumulated deformation by adopting time low-pass filtering;
s533 combines the results of step S531 and step S532 to obtain the information of the surface deformation of the target area.
In an embodiment of the present invention, the step S2 includes:
s21, registering the SAR images, setting a space-time baseline threshold, forming a short baseline set according to the space-time baseline threshold, and forming an interferogram set according to the principle of the short baseline set;
s22, converting DEM data of the SAR image into a terrain phase value, and differentiating the interferogram set and the terrain phase generated by the DEM to obtain a differential interferogram set;
s23, filtering, phase unwrapping and geocoding the differential interferogram to obtain an unwrapped interferogram set.
The beneficial effects of the invention are as follows:
(1) The phase change rate is obtained through a stacking time sequence algorithm, stable points of a target area are obtained, a functional relation is established based on the stable points, the functional relation is solved through a least square method, and the automatic correction of differential interference phase orbit errors and layered atmospheric phases of InSAR unwrapping can be achieved.
(2) The method can automatically judge the quality of the interferograms, screen high-quality interferograms and perform deformation calculation to ensure the ground deformation monitoring precision.
(3) By means of block processing, the invention can realize large-range and automatic ground surface deformation monitoring.
Drawings
Various other advantages and benefits will become apparent to those of ordinary skill in the art upon reading the following detailed description of the preferred embodiments. The drawings are only for purposes of illustrating the preferred embodiments and are not to be construed as limiting the invention. Also, like reference numerals are used to refer to like parts throughout the drawings. In the drawings:
FIG. 1 is a flow chart of a method of an embodiment of the present invention;
FIG. 2 is a technical roadmap for an embodiment of the present invention;
FIG. 3 is a flowchart of a joint calibration of orbital error and stratified atmospheric phase according to an embodiment of the present invention;
FIG. 4 is a flow chart of automated screening of interferogram quality according to an embodiment of the present invention;
FIG. 5 is a flowchart of short baseline set time-series deformation calculation and filtering according to an embodiment of the present invention;
FIG. 6 is a diagram illustrating the joint calibration effect of orbital error and laminar atmospheric phase according to an embodiment of the present invention;
FIG. 7 is a diagram illustrating an effect of an automatic screening of interferograms according to an embodiment of the present invention;
FIG. 8 is a graph showing deformation rate results obtained in accordance with an embodiment of the present invention.
Detailed Description
The following detailed description of embodiments of the present invention is provided in connection with the accompanying drawings and examples. The following examples are intended to illustrate the invention but are not intended to limit the scope of the invention.
Referring to fig. 1, the present invention provides an InSAR time sequence deformation monitoring method with automatic error correction, including:
s1, a target area is obtained, a time interval is selected, and SAR images of the time interval in the target area are selected.
S2, setting a space-time baseline threshold, constructing a differential interferogram set of the short baseline set based on the space-time baseline threshold, and performing phase unwrapping and geocoding to obtain an unwrapped interferogram set.
And S3, identifying stable points in the target area, and performing combined correction of the orbit error and the atmospheric phase on the unwrapping interference pattern set based on the stable points to obtain a corrected interference pattern set.
And S4, screening the correction interference pattern set according to the screening conditions to obtain a high-quality interference pattern set.
And S5, resolving and processing the high-quality interference pattern set based on a deformation algorithm to obtain ground deformation information of the target area.
Referring to fig. 2, the technical process of the present invention is: 1) And acquiring an SAR image in the target region, wherein the SAR image comprises an SLC head file, a track file and DEM data. 2) And carrying out registration operation on the SAR image, constructing a short baseline set, carrying out difference based on DEM data, and carrying out primary phase unwrapping and geocoding to obtain an unwrapped interferogram set. 3) And performing combined correction of the orbit error and the layered atmospheric phase on the unwrapping interference pattern set according to the stable point in the target area to obtain a corrected interference pattern set. 4) And performing quality detection and screening based on the unwrapping phase and coherence of the corrected interferogram set, obtaining interferograms with better quality according to screening conditions, taking the interferograms as a high-quality interferogram set, and screening out the interferograms which do not meet the screening conditions by putting the interferograms into a poor interferogram list. 5) And performing deformation calculation and filtering on the high-quality interference image set, and obtaining deformation rate and cumulative deformation information of the target area.
Specifically, an area requiring terrain monitoring is determined as a target area, and an SAR image in an appropriate time interval is acquired, which is also referred to as a main image.
Specifically, step S2 includes:
s21, registering the SAR images, setting a space-time baseline threshold, forming a short baseline set according to the space-time baseline threshold, and forming an interferogram set according to the principle of the short baseline set;
s22, converting DEM data of the SAR image into a terrain phase value, and differentiating the interferogram set and the terrain phase generated by the DEM to obtain a differential interferogram set;
s23, filtering, phase unwrapping and geocoding the differential interferogram to obtain an unwrapped interferogram set.
In an embodiment of the present invention, the step S2 includes the following steps:
1. the method comprises the steps of carrying out registration processing on an original SAR image, wherein the registration of the SAR image is also the registration of an interference complex image pair, and comprises coarse registration, pixel level registration and sub-pixel level registration.
Coarse registration
The first step is as follows: calculating the pixel of the center point of the main image
Coordinate values (x, y, z) in the orbital coordinate system. The azimuth time and the ground distance time are respectively calculated by utilizing the line number, the PRE value, the pixel number and the RSR value, and the vector of the satellite position is calculated by utilizing a cubicstation method according to the azimuth, the ground distance time, the speed and the position vector in an SLC head file. And then calculating coordinate values (x, y, z) of the central point pixel by utilizing multiple iterations to set a convergence threshold value for partial derivatives of the equation.
The second step: based on the generalLe equation, calculating the row and column values of the point on the secondary image
。
The doppler equation is as follows:
wherein the content of the first and second substances,
a position vector representing a satellite, based on the location of the satellite, and>
a position vector representing a pixel point, and>
。
the third step: according to
And &>
To calculate the offset between the main image and the auxiliary image>
。
。
Pixel level registration
The precision of pixel-level registration is 1 pixel, and the present embodiment adopts a window-based automatic registration method, and the registration operation is performed in the frequency domain of the primary and secondary images. The matching window of the window matching method is set according to the size of the main image, in the embodiment, 256 pixels by 256 pixels can be set, the search window is selected on the auxiliary image, and the quality evaluation indexes of the two windows are calculated in different rows and columns according to the integral pixel offset in the operation process of the search window.
Sub-pixel level registration
The sub-pixel level registration is also realized in the frequency domain, the registration precision reaches 1/8 pixel, and the steps are as follows:
1. firstly, performing oversample processing on a main image and an auxiliary image, and then interpolating the image pair by using an interpolation method, wherein the interpolation interval is 0.01 pixel.
2. A maximum coherence estimate is performed, which can also be implemented based on a window search method to calculate the offset.
3. And fitting the offset and resampling the auxiliary image, performing data fitting of the offset by adopting a least square method, calculating a coordinate conversion relation of the complex image pair by using a second-order polynomial, and resampling the auxiliary image according to the conversion relation.
2. Setting a space-time baseline threshold, freely combining the SAR images into a plurality of subsets according to the space-time baseline threshold to be used as a short baseline set, wherein the baseline distance between the SAR images is smaller and the image baseline distance between the subsets is larger in each subset, and generating an interferogram of the subsets, namely an interferogram set, according to interference conditions.
3. Obtaining DEM data in a target area, reading the DEM data, calculating DEM point intervals and interference fringe image element point intervals, calculating a ratio of the DEM point intervals and the interference fringe image element point intervals, then resampling the DEM data, obtaining an image element value of the DEM data by using a nearest neighbor interpolation method according to the obtained ratio, converting the DEM data into a radar coordinate system, assigning values to the converted DEM data by using a nearest neighbor interpolation method, and finally converting the DEM data into a phase value.
4. And (4) differentiating the interferogram set and the DEM generated terrain phase, namely subtracting the phase value of the DEM from the interferogram set to obtain a differential interferogram set.
5. The difference interferogram is filtered, and the Goldstein filtering algorithm is adopted in the embodiment, which specifically includes the following steps:
firstly, dividing a differential interference pattern into a plurality of mutually overlapped phase blocks, then converting the phase blocks into a frequency domain according to Fourier transform, and carrying out corresponding smoothing treatment on a power spectrum of the phase blocks:
wherein S {. Star } is a smoothing operator, F (u, v) and H (u, v) are frequency domain interferograms before and after filtering respectively, (u, v) are coordinate values of a certain pixel point in the interferogram in the frequency domain, alpha is a filtering parameter,
。
when a =0 is set as the control parameter,
indicates no filtering effect and when α =1, is present>
A low-pass filter consisting of a rectangular smoothing window is shown, indicating that the filtering effect is very strong. In addition, a larger phase block size and a higher value of the filter parameter α are advantageous for enhancing coherence of the interferogram, and overlapping phase blocks (overlap less than 75%) helps to reduce discontinuities at the interferogram boundaries.
Considering the nonuniformity of phase noise distribution in the interferogram, if the whole interferogram is filtered by a fixed alpha value, the adaptivity of filtering is easily reduced, and there may be over-filtering in high-coherence regions and under-filtering in low-coherence regions, so it is difficult to select a reasonable global filtering parameter alpha value, and for this, the average value of coherence values may be used instead of the alpha value, so that the filtering degree of the non-coherence region of the interferogram is stronger than that of the coherence region:
in the formula, parameter
Representing the block of valid phases in the coherence map, i.e. the sliding window minus the average of the absolute values of the coherence coefficients over the overlap region.
6. Phase unwrapping is performed on the differential interferogram, a network-based phase unwrapping algorithm is adopted in the embodiment, and a wrapping function is defined as:
wherein, (i, j) is a pixel point,
for winding phase>
Representing the unwrapped phase, the residual is defined as: />
The objective function is then expressed as:
wherein the content of the first and second substances,
is a priori the weight of the residual. The following variable substitutions were made:
substituting equation (7) into equation (6), the non-linear minimization problem of equation (6) can be converted into a linear minimization problem:
according to the above, the differential interferograms can be phase unwrapped and geocoded to obtain a set of unwrapped interferograms.
Referring to fig. 3, specifically, step S3 includes:
s31 obtains the average coherence of the unwrapped interferogram and sets a first coherence threshold.
S32, calculating the phase standard deviation of the unwrapped interferogram and setting a first phase standard deviation threshold.
S33, quality screening is carried out on the unwrapping interferogram set according to the first coherence threshold and the first phase standard deviation threshold, and if the average coherence of the unwrapping interferogram is smaller than the first coherence threshold or the phase standard deviation of the unwrapping interferogram is larger than the first phase standard deviation threshold, the unwrapping interferogram is screened out, and the unwrapping interferogram set after quality screening is obtained.
In an embodiment of the invention, when the unwrapping interferogram is obtained, the coherence of the unwrapping interferogram can be obtained at the same time, and the average coherence Coh of the unwrapping interferogram is calculated 1 Setting a first coherence threshold
Simultaneously calculating the phase standard deviation Std of the unwrapped interferogram 1 And sets a first phase standard deviation threshold->
If yes, is>
Or->
That is, if the average coherence of the unwrapping interferogram is smaller than the coherence threshold and the phase standard deviation of the unwrapping interferogram is larger than the phase standard deviation threshold, the quality of the unwrapping interferogram is poor, the unwrapping interferogram is screened out, and the poor unwrapping interferogram does not participate in the deformation calculation of the subsequent stacking time sequence.
S34, stacking time sequence InSAR processing is carried out on the unwrapped interferogram set after the quality screening, and the phase change rate is obtained according to a stacking time sequence algorithm.
S35, setting a threshold interval of the stable point, obtaining the deformation quantity of all points in the target area according to the phase change rate, and if the deformation quantity of a certain point is distributed in the threshold interval of the stable point, the point is the stable point.
S36, calculating the joint phase of the orbit error and the layered atmosphere according to the position information of the stable point, and performing phase correction on the unwrapping interferogram set by taking the joint phase as an error signal to obtain a corrected interferogram set.
Specifically, the stacked time sequence InSAR processing is performed on the set of unwrapped interferograms after the quality screening, and in order to enable most points to be subjected to deformation calculation, the minimum number of stacked interferograms is set in the embodiment
And (4) threshold, carrying out stacking time sequence deformation calculation by using intermittent coherent points. When the number of stacked interferograms->
Then, the phase change rate obtained by stacking is:
wherein the content of the first and second substances,
representing the rate of change of phase, N being the number of unwrapped interferograms>
For the time base of the jth unwrapped interferogram, the->
Is the phase value of the jth unwrapped interferogram. The standard deviation is:
after the phase change rate is obtained, deformation quantities of stable points in the region are considered to be distributed near the deformation median of all the points, so that when a stable point threshold interval is set, the threshold interval is the median of the deformation quantities of all the points in the target region. And, a stable point proportion threshold value is set
As a stable point selection criterion. When the deformation quantity at a certain point is distributed over all point deformation quantities->
When in the interval, this point is considered as a stable point.
After the stable point is obtained, establishing a functional relation among the unwrapping phase, the elevation and the position according to the position information of the stable point, wherein the functional relation is as follows:
wherein x and y are the row number and the column number of the picture element in the geographic coordinate system respectively,
is the unwrapped phase at point (x, y), h (x, y) is the elevation of point (x, y), and->
Is the random phase error of point (x, y), a i =(a 0 ,a 1 ,...a 6 ) Is the correlation coefficient.
And solving the joint phase of the orbit error and the layered atmosphere by a least square method, and removing the joint phase from the unwrapping phase of the formula to obtain a corrected interferogram set.
Specifically, step S4 includes:
s41, average coherence of the correction interferogram is obtained, a second coherence threshold value is set, phase standard deviation of the correction interferogram is calculated, and a second phase standard deviation threshold value is set.
S42, quality screening is carried out on the corrected interferogram set according to the second coherence threshold and the second phase standard deviation threshold, if the average coherence of the corrected interferogram is smaller than the second coherence threshold or the phase standard deviation of the corrected interferogram is larger than the second phase standard deviation threshold, the corrected interferogram is screened out, and the corrected interferogram set after quality screening is obtained.
S43, screening the corrected interference pattern set subjected to quality screening by adopting a triangular phase closed-loop method, forming the corrected interference pattern set subjected to quality screening into a triangular phase closed-loop network formed by one or more triangles, and generating three unwrapped interference differential phases by each triangular phase closed-loop;
s44, setting a closed loop phase threshold, and calculating the root mean square of the corresponding closed loop phase according to the three unwrapping interference differential phases of each group of triangular phase closed loops:
when the root mean square of the closed loop phase is less than the closed loop phase threshold, the three unwrapped interfering differential phases are retained,
when the root mean square of the closed loop phase is larger than the closed loop phase threshold, rejecting the three unwrapped interference differential phases;
s45, repeating the step S44 to obtain a high-quality interferogram set.
Referring to FIG. 4, in an embodiment of the present invention, when obtaining the corrected interferogram, the coherence can be obtained at the same time, and the average coherence Coh of the corrected interferogram is calculated 2 Setting a second coherence threshold
Simultaneously calculating the phase standard deviation Std of the corrected interferogram 2 And setting a second phase standard deviation threshold->
If yes, is>
Or
If the average coherence of the corrected interferograms is smaller than the second coherence threshold and the phase standard deviation of the corrected interferograms is larger than the second phase standard deviation threshold, the quality of the corrected interferograms is poor, and the corrected interferograms are placed into the poor corrected interferogram set 1 to be screened out, so that the corrected interferogram set after quality screening is obtained. The first coherence threshold and the second coherence threshold may be the same, or the first phase standard deviation threshold and the second phase standard deviation threshold may be the same.
After the correction interference pattern set is subjected to quality screening, whether a triangular closed loop exists in the correction interference patterns subjected to the quality screening is judged, the correction interference pattern set with the triangular closed loop is screened according to a triangular phase closed loop method, the correction interference patterns which do not meet the conditions are placed into the poor-quality correction interference pattern set 2 to be screened, the correction interference patterns which meet the conditions have high-quality unwrapping phases, and the correction interference patterns are used as the good-quality interference pattern set.
In an embodiment of the present invention, if the original SAR image is 3, then the SAR images 1, 2, 3 are used as a set of triangular phase closed loops to generate three unwrapping interference differential phases
: setting a closed loop phase threshold->
And calculating the root mean square of the closed loop phase as a judgment standard of the quality of the unwrapping phase. When it satisfies
Temporarily retaining all unwinding phases of the closed loop; when/is>
Then related to rejecting closed loop
And unwinding the phase.
In another embodiment of the present invention, the original SAR image is 4, and the set of the corrected interferograms forms a triangular phase closed-loop network formed by two triangles, and the unwrapping phases thereof are respectively:
when the reference for unwinding is in agreement, there is theoretically->
，
. In practical cases, due to interference pattern filteringThe right side of the equation is often not 0 due to windup errors. In order to select an interferogram with serious unwrapping error, a closed loop phase threshold is set
And calculating the root mean square of the closed loop phase as a judgment standard of the quality of the unwrapping phase. When it is satisfied with
Temporarily retaining all the unwrapping phases of the closed loop; when in use
If so, the part associated with the closing ring is rejected>
And unwinding the phase.
Through the steps, a high-quality interferogram set is obtained.
Referring to fig. 5, in the embodiment of the present invention, the deformation algorithm is a weighted least squares deformation calculation algorithm with an additional smooth constraint condition, specifically, step S5 includes:
s51, the target area is partitioned, high coherence points in each area block are obtained, and unwrapping phase and coherence information of each high coherence point are extracted.
In particular, the high coherence points have stable and time invariant radar scattering properties, have reliable phase information, and are not affected by the orbital spatial separation of the interferometric system.
In an embodiment of the present invention, the high coherence Point (PS) is selected as follows:
and screening high coherence points according to the time sequence correlation coefficient, and according to the following steps:
wherein the content of the first and second substances,
is a coherence factor and->
M and S respectively represent local pixel information blocks of the two SAR images constituting the interference pair, a denotes a complex conjugate operator, and M and n denote the size of the local moving window. Coherence factor->
The larger the noise of the interference phase is. />
And calculating the coherence coefficient of each pixel in the time interval. If N +1 SAR images are provided, N interferograms can be formed, and each pixel forms a coherence coefficient sequence
。
Calculating the average value of the coherent coefficient of the pixel in the time interval:
then setting a suitable threshold
And extracting the image elements with the average coherence coefficient larger than or equal to the threshold value as high coherence points.
Through the selection process, the pixel points with serious decorrelation can be filtered.
S52, establishing a functional relation among the unwrapping phase, the deformation rate and the elevation residual error in each area block, solving the deformation rate of each high-coherence point according to the coherence information, solving according to the functional relation to obtain the elevation residual error, and splicing the results of each area block to obtain the accumulated deformation at each moment.
Specifically, step S52 includes:
for a single high coherence point to be solved, the observation equation is:
wherein, the first and the second end of the pipe are connected with each other,
is a time interval, is->
At a rate of deformation>
Is an elevation residual->
For unwrapping phase, M is the number of good quality interferograms, N +1 is the number of SAR images, θ is the angle of incidence, B is the vertical baseline of the good quality interferograms, R is the distance between the SAR image radar antenna and the ground target, π is a constant, λ is the wavelength.
Introducing virtual observations
Adding a smooth constraint condition with equal acceleration for solving as initial data, wherein the smooth constraint condition is as follows:
in the formula
Representing the deformation rate of the ground object at the kth time instant. And (3) introducing a smoothing parameter alpha to control the constraint strength, and combining an observation equation and two equations of smoothing constraint conditions, namely an equation (14) and an equation (15):
wherein
For smoothing the constraint coefficients, the larger the value, the smoother the solution time series cumulative deformation. />
Is a rate term coefficient>
For elevation residual term coefficients>
For the smooth constraint matrix, i.e.:
solving the rate item and the elevation residual error, and splicing the result of each area block to obtain the accumulated deformation at each moment.
And S53, carrying out turbulence atmospheric correction and noise removal processing on the accumulated deformation at each moment according to a filtering method to obtain the earth surface deformation information of the target area.
Specifically, step S53 includes:
s531, processing the accumulated deformation at each moment according to spatial low-pass filtering to remove a noise phase in the accumulated deformation;
s532, removing turbulence atmospheric components in the accumulated deformation by adopting time low-pass filtering;
s533 combines the results of step S531 and step S532 to obtain the information of the surface deformation of the target area.
One embodiment is illustrated in the following:
taking Beijing as an example, the target area is used as the target area, and SAR data of 23 sentinels in total are collected in the area from 3 months in 2019 to 12 months in 2019. The attitude of the satellite is in the ascending orbit direction, the heading angle of the satellite is-13.47 degrees, the incident angle of the radar is 43.75 degrees, and the time interval is uniform.
And registering the 23-scene SAR image, constructing an interferogram network by taking a space-time baseline threshold as a standard, acquiring DEM data of the interferogram network, performing pretreatment such as differential interference, phase unwrapping, geocoding and the like, and producing an unwrapped differential interference phase in Beijing area to obtain an unwrapped interferogram set.
The stable point in the target area is identified by using the stacking time sequence InSAR method, and the orbit error and atmospheric phase joint correction is performed on all the unwrapping phases based on the stable point information, and the obtained effect graph is shown in FIG. 6.
And (3) performing quality evaluation on the interferograms based on the coherence, the phase standard deviation and the phase closed-loop threshold, and screening to obtain a high-quality interferogram set with better quality, wherein the automatic screening result is shown in fig. 7.
And performing weighted least square deformation calculation and filtering processing with additional smooth constraint conditions on the high-quality unwrapping phase to obtain reliable ground deformation information, wherein the obtained deformation result is shown in figure 8.
The above embodiments are only for illustrating the invention and are not to be construed as limiting the invention, and those skilled in the art can make various changes and modifications without departing from the spirit and scope of the invention, therefore, all equivalent technical solutions also belong to the scope of the invention, and the scope of the invention is defined by the claims.
Claims (10)
1. An InSAR time sequence deformation monitoring method with automatic error correction is characterized by comprising the following steps:
s1, acquiring a target area, selecting a time interval, and selecting an SAR image of the time interval in the target area;
s2, setting a space-time baseline threshold, constructing a differential interferogram set of a short baseline set based on the space-time baseline threshold, and performing phase unwrapping and geocoding to obtain an unwrapped interferogram set;
s3, identifying stable points in the target area, and performing combined correction of orbit errors and atmospheric phases on the unwrapping interference pattern set based on the stable points to obtain a corrected interference pattern set;
s4, screening the correction interference pattern set according to screening conditions to obtain a high-quality interference pattern set;
and S5, resolving and processing the high-quality interference pattern set based on a deformation algorithm to obtain ground deformation information of the target area.
2. The method according to claim 1, wherein step S3 comprises:
s31, acquiring the average coherence of the unwrapping interferogram, and setting a first coherence threshold;
s32, calculating the phase standard deviation of the unwrapping interferogram, and setting a first phase standard deviation threshold value;
s33, performing quality screening on the unwrapping interferogram set according to the first coherence threshold and the first phase standard deviation threshold, and screening the unwrapping interferogram if the average coherence of the unwrapping interferogram is smaller than the first coherence threshold or the phase standard deviation of the unwrapping interferogram is larger than the first phase standard deviation threshold to obtain the unwrapping interferogram set after quality screening;
s34, stacking time sequence InSAR processing is carried out on the unwrapped interferogram set after quality screening, and phase change rate is obtained according to a stacking time sequence algorithm;
s35, setting a threshold interval of a stable point, obtaining deformation quantities of all points in the target area according to the phase change rate, and if the deformation quantities of a certain point are distributed in the threshold interval of the stable point, the point is the stable point;
s36, calculating the joint phase of the orbit error and the layered atmosphere according to the position information of the stable point, and performing phase correction on the unwrapping interferogram set by taking the joint phase as an error signal to obtain a corrected interferogram set.
3. The method of claim 2, wherein step S36 comprises:
establishing a functional relation among the unwrapping phase, the elevation and the position according to the position information of the stable point, wherein the formula is as follows:
wherein x and y are the row number and the column number of the picture element in the geographic coordinate system respectively,
is the unwrapped phase at point (x, y), h (x, y) is the elevation of point (x, y), and->
Is the random phase error of point (x, y), a i =(a 0 ,a 1 ,...a 6 ) Is the correlation coefficient;
and solving the joint phase of the orbit error and the layered atmosphere by a least square method, and removing the joint phase from the unwrapping phase of the above formula to obtain a third interferogram.
4. The method of claim 2, wherein in step S34, the stacking timing algorithm is: setting a minimum stacked interferogram number threshold
Stacking time sequence deformation calculation is carried out on the unwrapping interferogram set subjected to quality screening according to an intermittent mode, and when the number of the stacked interferograms is greater than or equal to>
Then, the stack yields the phase change rate:
wherein the content of the first and second substances,
representing the rate of change of phase, N being the number of unwrapped interferograms>
For the time base of the jth unwrapped interferogram, the->
Is the phase value of the jth unwrapped interferogram.
5. The method according to claim 2, wherein in step S35, the stable point threshold interval is around the median of the deformation quantities of all points in the target region.
6. The method according to claim 2, wherein step S4 comprises:
s41, acquiring the average coherence of the correction interferogram, setting a second coherence threshold, calculating the phase standard deviation of the correction interferogram, and setting a second phase standard deviation threshold;
s42, quality screening is carried out on the corrected interferogram set according to the second coherence threshold and the second phase standard deviation threshold, if the average coherence of the corrected interferogram is smaller than the second coherence threshold or the phase standard deviation of the corrected interferogram is larger than the second phase standard deviation threshold, the corrected interferogram is screened out, and a corrected interferogram set after quality screening is obtained;
s43, screening the corrected interference pattern set after quality screening by adopting a triangular phase closed-loop method, forming the corrected interference pattern set after quality screening into a triangular phase closed-loop network formed by one or more triangles, and generating three unwrapping interference differential phases by each triangular phase closed-loop;
s44, setting a closed loop phase threshold, and calculating the root mean square of the corresponding closed loop phase according to the three unwrapping interference differential phases of each group of triangular phase closed loops:
when the root mean square of the closed loop phase is less than the closed loop phase threshold, the three unwrapped interfering differential phases are retained,
when the root mean square of the closed loop phase is larger than the closed loop phase threshold, rejecting the three unwrapped interference differential phases;
s45, repeating the step S44 to obtain a high-quality interferogram set.
7. The method according to claim 1, wherein the deformation algorithm is a weighted least squares deformation solving algorithm with additional smooth constraint, and the step S5 comprises:
s51, partitioning a target area, acquiring high coherence points in each area block, and extracting unwrapping phases and coherence information of each high coherence point;
s52, establishing a functional relation among the unwrapping phase, the deformation rate and the elevation residual error in each area block, solving the deformation rate of each high-coherence point according to the coherence information, solving according to the functional relation to obtain the elevation residual error, and splicing the results of each area block to obtain the accumulated deformation at each moment;
and S53, carrying out turbulence atmospheric correction and noise removal processing on the accumulated deformation at each moment according to a filtering method to obtain the earth surface deformation information of the target area.
8. The method of claim 7, wherein step S52 comprises:
for a single high coherence point to be solved, the observation equation is:
wherein the content of the first and second substances,
is a time interval, is->
Is changed at the change rate->
In the form of an elevation residual,
for unwrapping phase, M is the number of high-quality interferograms, N +1 is the number of SAR images, θ is the angle of incidence, B is the vertical baseline of the high-quality interferograms, R is the distance between the SAR image radar antenna and the ground target, π is a constant, λ is the wavelength;
introducing virtual observations
Adding a smooth constraint condition with equal acceleration for solving as initial data, wherein the smooth constraint condition is as follows:
in the formula
Expressing the deformation rate of the ground target at the kth moment, introducing a smoothing parameter alpha to control the constraint strength, and combining an observation equation with a smoothing constraint condition:
wherein
For smooth constraint coefficients>
Is a rate factor, is based on the rate of the signal>
Is an elevation residual term coefficient>
Is a smooth constraint matrix;
solving the rate item and the elevation residual error, and splicing the result of each area block to obtain the accumulated deformation at each moment.
9. The method according to claim 7, wherein step S53 comprises:
s531 processing the accumulated deformation of each moment according to spatial low filtering to remove a noise phase in the accumulated deformation;
s532, removing turbulence atmospheric components in the accumulated deformation by adopting time low-pass filtering;
s533 obtains the surface deformation information of the target area by combining the results of step S531 and step S532.
10. The method according to claim 1, wherein step S2 comprises:
s21, registering the SAR images, setting a space-time baseline threshold, forming a short baseline set according to the space-time baseline threshold, and forming an interferogram set according to the principle of the short baseline set;
s22, converting DEM data of the SAR image into a terrain phase value, and differentiating the interferogram set and the terrain phase generated by the DEM to obtain a differential interferogram set;
s23, filtering, phase unwrapping and geocoding the differential interferogram to obtain an unwrapped interferogram set.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310054817.6A CN115856889B (en) | 2023-02-03 | 2023-02-03 | InSAR time sequence deformation monitoring method capable of automatically correcting errors |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310054817.6A CN115856889B (en) | 2023-02-03 | 2023-02-03 | InSAR time sequence deformation monitoring method capable of automatically correcting errors |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN115856889A true CN115856889A (en) | 2023-03-28 |
| CN115856889B CN115856889B (en) | 2023-05-19 |
Family
ID=85657483
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202310054817.6A Active CN115856889B (en) | 2023-02-03 | 2023-02-03 | InSAR time sequence deformation monitoring method capable of automatically correcting errors |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN115856889B (en) |
Cited By (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN116148855A (en) * | 2023-04-04 | 2023-05-23 | 之江实验室 | Method and system for removing atmospheric phase and resolving deformation of time sequence InSAR |
| CN116299245A (en) * | 2023-05-11 | 2023-06-23 | 中山大学 | Time sequence InSAR deformation rate result self-adaptive mosaic correction method |
| CN116403116A (en) * | 2023-06-08 | 2023-07-07 | 江苏省地质调查研究院 | D-InSAR earth surface deformation space-time characteristic fine sensing method for comprehensively planning multiple observation indexes |
| CN116559866A (en) * | 2023-07-11 | 2023-08-08 | 南京天辰礼达电子科技有限公司 | Ground-based synthetic aperture radar atmosphere compensation method |
| CN116736306A (en) * | 2023-08-15 | 2023-09-12 | 成都理工大学 | Time sequence radar interference monitoring method based on third high-resolution |
| CN116736305A (en) * | 2023-08-14 | 2023-09-12 | 北京观微科技有限公司 | Ground object deformation determination method and device, electronic equipment and storage medium |
Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103675790A (en) * | 2013-12-23 | 2014-03-26 | 中国国土资源航空物探遥感中心 | Method for improving earth surface shape change monitoring precision of InSAR (Interferometric Synthetic Aperture Radar) technology based on high-precision DEM (Digital Elevation Model) |
| CN106950556A (en) * | 2017-05-03 | 2017-07-14 | 三亚中科遥感研究所 | Heritage area deformation monitoring method based on distributed diffusion body sequential interference SAR technology |
| CN111474544A (en) * | 2020-03-04 | 2020-07-31 | 广东明源勘测设计有限公司 | Landslide deformation monitoring and early warning method based on SAR data |
| US20200394780A1 (en) * | 2017-06-15 | 2020-12-17 | The University Of Nottingham | Land deformation measurement |
| CN112213723A (en) * | 2020-09-11 | 2021-01-12 | 速度时空信息科技股份有限公司 | Method for monitoring landslide in real time by using SBAS technology |
| CN115113203A (en) * | 2022-07-20 | 2022-09-27 | 江苏省水利科学研究院 | Method for removing InSAR atmospheric phase |
-
2023
- 2023-02-03 CN CN202310054817.6A patent/CN115856889B/en active Active
Patent Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103675790A (en) * | 2013-12-23 | 2014-03-26 | 中国国土资源航空物探遥感中心 | Method for improving earth surface shape change monitoring precision of InSAR (Interferometric Synthetic Aperture Radar) technology based on high-precision DEM (Digital Elevation Model) |
| CN106950556A (en) * | 2017-05-03 | 2017-07-14 | 三亚中科遥感研究所 | Heritage area deformation monitoring method based on distributed diffusion body sequential interference SAR technology |
| US20200394780A1 (en) * | 2017-06-15 | 2020-12-17 | The University Of Nottingham | Land deformation measurement |
| CN111474544A (en) * | 2020-03-04 | 2020-07-31 | 广东明源勘测设计有限公司 | Landslide deformation monitoring and early warning method based on SAR data |
| CN112213723A (en) * | 2020-09-11 | 2021-01-12 | 速度时空信息科技股份有限公司 | Method for monitoring landslide in real time by using SBAS technology |
| CN115113203A (en) * | 2022-07-20 | 2022-09-27 | 江苏省水利科学研究院 | Method for removing InSAR atmospheric phase |
Non-Patent Citations (1)
| Title |
|---|
| 陆世安;: "珠三角地区短基线地表形变估计――以珠海市为例" * |
Cited By (12)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN116148855A (en) * | 2023-04-04 | 2023-05-23 | 之江实验室 | Method and system for removing atmospheric phase and resolving deformation of time sequence InSAR |
| CN116148855B (en) * | 2023-04-04 | 2023-07-21 | 之江实验室 | Method and system for removing atmospheric phase and resolving deformation of time sequence InSAR |
| CN116299245A (en) * | 2023-05-11 | 2023-06-23 | 中山大学 | Time sequence InSAR deformation rate result self-adaptive mosaic correction method |
| CN116299245B (en) * | 2023-05-11 | 2023-07-28 | 中山大学 | Time sequence InSAR deformation rate result self-adaptive mosaic correction method |
| CN116403116A (en) * | 2023-06-08 | 2023-07-07 | 江苏省地质调查研究院 | D-InSAR earth surface deformation space-time characteristic fine sensing method for comprehensively planning multiple observation indexes |
| CN116403116B (en) * | 2023-06-08 | 2023-08-29 | 江苏省地质调查研究院 | D-InSAR earth surface deformation space-time characteristic fine sensing method for comprehensively planning multiple observation indexes |
| CN116559866A (en) * | 2023-07-11 | 2023-08-08 | 南京天辰礼达电子科技有限公司 | Ground-based synthetic aperture radar atmosphere compensation method |
| CN116559866B (en) * | 2023-07-11 | 2023-09-29 | 南京天辰礼达电子科技有限公司 | Ground-based synthetic aperture radar atmosphere compensation method |
| CN116736305A (en) * | 2023-08-14 | 2023-09-12 | 北京观微科技有限公司 | Ground object deformation determination method and device, electronic equipment and storage medium |
| CN116736305B (en) * | 2023-08-14 | 2023-10-27 | 北京观微科技有限公司 | Ground object deformation determination method and device, electronic equipment and storage medium |
| CN116736306A (en) * | 2023-08-15 | 2023-09-12 | 成都理工大学 | Time sequence radar interference monitoring method based on third high-resolution |
| CN116736306B (en) * | 2023-08-15 | 2023-10-24 | 成都理工大学 | Time sequence radar interference monitoring method based on third high-resolution |
Also Published As
| Publication number | Publication date |
|---|---|
| CN115856889B (en) | 2023-05-19 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN115856889A (en) | InSAR time sequence deformation monitoring method capable of automatically correcting errors | |
| CN111208512B (en) | Interferometric measurement method based on video synthetic aperture radar | |
| CN109782282B (en) | Time series InSAR analysis method integrating troposphere atmospheric delay correction | |
| CN112986993B (en) | InSAR deformation monitoring method based on space constraint | |
| CN109633648B (en) | Multi-baseline phase estimation device and method based on likelihood estimation | |
| CN103487809B (en) | A kind of based on BP algorithm and time become the airborne InSAR data disposal route of baseline | |
| CN111273293B (en) | InSAR residual motion error estimation method and device considering terrain fluctuation | |
| CN113624122A (en) | Bridge deformation monitoring method fusing GNSS data and InSAR technology | |
| Li et al. | A new analytical method for estimating Antarctic ice flow in the 1960s from historical optical satellite imagery | |
| Crosetto et al. | Radargrammetry and SAR interferometry for DEM generation: validation and data fusion. | |
| CN113093184B (en) | Interferometric measurement method based on video synthetic aperture radar | |
| CN107945216A (en) | More images joint method for registering based on least-squares estimation | |
| Knöpfle et al. | Mosaicking of digital elevation models derived by SAR interferometry | |
| CN109633639A (en) | The high-precision rapid registering method of TOPSAR interference data | |
| CN110780327B (en) | Marine target cooperative positioning method based on satellite-borne AIS and infrared camera | |
| CN117523404B (en) | Method, device, terminal and storage medium for monitoring dynamic change of storage yard | |
| CN113552565B (en) | Phase unwrapping method for SAR data high noise and large gradient change area | |
| CN102944308B (en) | Attitude error correcting method of time-space joint modulation interference imaging spectrometer | |
| CN108008382B (en) | A kind of method of more base spaceborne interferometric SAR systematic survey mountain terrains | |
| CN109886910A (en) | External digital elevation model DEM modification method and device | |
| CN111598929B (en) | Two-dimensional unwrapping method based on time sequence differential interference synthetic aperture radar data | |
| CN109946682B (en) | GF3 data baseline estimation method based on ICESat/GLAS | |
| Charrier et al. | Analysis of dense coregistration methods applied to optical and SAR time-series for ice flow estimations | |
| JP4937791B2 (en) | Radar image processing device | |
| CN116299453A (en) | Satellite-borne SAR non-trace-along mode interference elevation measurement 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 |