CN110658521A - Winding phase-based GBInSAR atmospheric correction method and system - Google Patents

Winding phase-based GBInSAR atmospheric correction method and system Download PDF

Info

Publication number
CN110658521A
CN110658521A CN201910982029.7A CN201910982029A CN110658521A CN 110658521 A CN110658521 A CN 110658521A CN 201910982029 A CN201910982029 A CN 201910982029A CN 110658521 A CN110658521 A CN 110658521A
Authority
CN
China
Prior art keywords
difference
adjacent
phase
coherence
atmospheric
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
Application number
CN201910982029.7A
Other languages
Chinese (zh)
Other versions
CN110658521B (en
Inventor
杨红磊
刘杰
彭军还
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
China University of Geosciences Beijing
Original Assignee
China University of Geosciences Beijing
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by China University of Geosciences Beijing filed Critical China University of Geosciences Beijing
Priority to CN201910982029.7A priority Critical patent/CN110658521B/en
Publication of CN110658521A publication Critical patent/CN110658521A/en
Application granted granted Critical
Publication of CN110658521B publication Critical patent/CN110658521B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • GPHYSICS
    • G01MEASURING; TESTING
    • G01SRADIO 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/00Systems 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/88Radar or analogous systems specially adapted for specific applications
    • G01S13/89Radar or analogous systems specially adapted for specific applications for mapping or imaging
    • G01S13/90Radar or analogous systems specially adapted for specific applications for mapping or imaging using synthetic aperture techniques, e.g. synthetic aperture radar [SAR] techniques
    • G01S13/9021SAR image post-processing techniques
    • G01S13/9023SAR image post-processing techniques combined with interferometric techniques
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01SRADIO 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
    • G01S7/00Details of systems according to groups G01S13/00, G01S15/00, G01S17/00
    • G01S7/02Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S13/00
    • G01S7/40Means for monitoring or calibrating

Landscapes

  • Engineering & Computer Science (AREA)
  • Remote Sensing (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Physics & Mathematics (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • General Physics & Mathematics (AREA)
  • Electromagnetism (AREA)
  • Radar Systems Or Details Thereof (AREA)

Abstract

The invention discloses a winding phase-based GBInSAR atmospheric correction method and system. The method comprises the following steps: acquiring two single-view complex images and a differential interference image of the two single-view complex images in the same region; selecting high coherence points in the interference pattern; obtaining a phase difference value and a distance difference value which are actually measured between adjacent high-coherence points to obtain an adjacent phase difference and an adjacent distance difference; calculating an atmospheric delay coefficient according to the adjacent phase difference and the adjacent distance difference; calculating an atmospheric delay phase according to the atmospheric delay coefficient; according to the difference value of the atmospheric delay phase and the adjacent phase difference, eliminating the adjacent phase difference and the adjacent distance difference corresponding to the difference value of the super threshold to obtain the adjacent phase difference and the adjacent distance difference after rough difference removal; updating the atmospheric delay coefficient according to the adjacent phase difference and the adjacent distance difference after the coarse difference is removed to obtain the atmospheric delay coefficient after the coarse difference is removed; and calculating the atmosphere delay phase after coarse difference correction by using the atmosphere delay coefficient after coarse difference removal. The invention can get rid of the dependence on phase unwrapping.

Description

Winding phase-based GBInSAR atmospheric correction method and system
Technical Field
The invention relates to the technical field of microwave remote sensing measurement, in particular to a GBInSAR atmospheric correction method and system based on winding phase.
Background
Based on the ground-based synthetic aperture radar interferometry (GB-InSAR), is an active microwave radar detection technique, originating from the synthetic aperture radar interferometry (InSAR) technique. The method has the technical advantages of small scale, high flexibility, high precision and the like. In recent years, GB-InSAR technology shows huge deformation monitoring potential in many fields, such as landslide, glacier movement, mining area settlement and the like. However, there are many factors that affect the monitoring accuracy of the GB-InSAR technique, such as noise, system frequency instability, and atmospheric delay, which is a critical factor in interferograms over long time intervals. Studies have shown that a change in relative humidity of 1% at a temperature of 20 c, 1km from the radar, can result in a measurement error of 2 mm. In order to improve the monitoring accuracy of GB-InSAR, the atmospheric delay phase must be corrected well.
To date, there are two main approaches to atmospheric delay correction. The first is a meteorological data compensation method. A mathematical model of meteorological parameters and atmospheric delay coefficients is established by obtaining temperature, humidity and barometric pressure data in a scene. Then, the atmospheric phase is calculated from the atmospheric delay coefficient. The meteorological data compensation method is influenced by the acquired meteorological parameter acquisition position and precision, and the precision is low. The second method is to model the distance as a function of atmospheric phase using high coherence points or artificial control points. Where atmospheric changes are simple, a linear or quadratic model is built between the distance and the phase of the atmospheric control point. The second method can improve accuracy relative to the first method. However, in this method, if there is a phase jump in the interferogram, a process dependent on phase unwrapping is required. Phase unwrapping is a difficult point in the InSAR technology, and different unwrapping methods have different influences on the result; furthermore, in an interferogram with low coherence, it is difficult to obtain a good phase unwrapping result, thereby affecting the accuracy of atmospheric delay correction.
Disclosure of Invention
The invention aims to provide a winding phase-based GBInSAR atmospheric correction method and a winding phase-based GBInSAR atmospheric correction system, which are free from dependence on phase unwrapping, thereby avoiding the influence of the phase unwrapping on atmospheric delay correction and improving the accuracy of the atmospheric delay correction.
In order to achieve the purpose, the invention provides the following scheme:
a winding phase-based GBInSAR atmospheric correction method comprises the following steps:
acquiring two single-view complex images obtained by detecting the same area at different moments and an interference pattern formed by differential interference of the two single-view complex images;
selecting high coherence points in the interference pattern according to coherence of each pixel point in the interference pattern;
acquiring actually measured phase difference values and distance difference values between adjacent high-coherence points in the same azimuth direction to obtain adjacent phase differences and adjacent distance differences;
calculating an atmospheric delay coefficient by using a least square regression method according to the adjacent phase difference and the adjacent distance difference;
calculating an atmospheric delay phase according to the atmospheric delay coefficient;
according to the difference value of the atmospheric delay phase and the adjacent phase difference, eliminating the adjacent phase difference and the adjacent distance difference corresponding to the difference value exceeding a preset threshold value to obtain the adjacent phase difference and the adjacent distance difference after rough difference removal;
updating the atmospheric delay coefficient by using a least square regression method according to the adjacent phase difference and the adjacent distance difference after the rough difference is removed to obtain the atmospheric delay coefficient after the rough difference is removed;
and calculating the atmosphere delay phase after the gross error correction by using the atmosphere delay coefficient after the gross error removal.
Optionally, after obtaining two single-view complex images obtained by detecting the same area at different times and an interferogram formed by differential interference of the two single-view complex images, before selecting a high coherence point in the interferogram according to coherence of each pixel point in the interferogram, the method further includes:
and performing Goldstein filtering processing on the interferogram.
Optionally, the selecting a high coherence point in the interference pattern according to coherence of each pixel point in the interference pattern specifically includes:
calculating the coherence of each pixel in the interference pattern, wherein the coherence represents the similarity of echo phases between two pixel points forming a certain pixel point in the interference pattern through differential interference;
and selecting the pixel points with the coherence exceeding a preset threshold value as high coherence points.
Optionally, the performing Goldstein filtering processing on the interferogram specifically includes:
dividing the interference pattern into a plurality of interference pattern blocks by using a sliding window, wherein the overlapping rate of any two adjacent interference pattern blocks is greater than or equal to a preset rate;
fourier transform is carried out on each interference pattern block to obtain a spectrogram of each interference pattern block;
smoothing each spectrogram to obtain a smoothed amplitude;
processing the frequency spectrum of each spectrogram by using the smoothed amplitude value to obtain a processed spectrogram;
and performing inverse Fourier transform on the processed spectrogram to obtain a filtered interferogram.
Optionally, after obtaining the actually measured phase difference value and distance difference value between the adjacent high coherence points in the same azimuth direction to obtain the adjacent phase difference and adjacent distance difference, before calculating the atmospheric delay coefficient by using a least square regression method according to the adjacent phase difference and adjacent distance difference, the method further includes:
and eliminating the adjacent distance difference higher than the preset distance difference threshold value and the corresponding adjacent phase difference data.
A winding phase based GBInSAR atmospheric correction system comprising:
the image acquisition module is used for acquiring two single-view complex images obtained by detecting the same area by the ground SAR at different moments and an interference image formed by differential interference of the two single-view complex images;
the coherence calculation module is used for selecting high coherence points in the interference pattern according to the coherence of each pixel point in the interference pattern;
the actual measurement difference value acquisition module is used for acquiring actual measurement phase difference values and distance difference values between adjacent high-coherence points in the same azimuth direction to obtain adjacent phase differences and adjacent distance differences;
the atmosphere delay coefficient calculation module is used for calculating an atmosphere delay coefficient by utilizing a least square regression method according to the adjacent phase difference and the adjacent distance difference;
the atmosphere delay phase calculation module is used for calculating an atmosphere delay phase according to the atmosphere delay coefficient;
the rough difference removing module is used for removing the adjacent phase difference and the adjacent distance difference corresponding to the difference value exceeding a preset threshold according to the difference value of the atmospheric delay phase and the adjacent phase difference to obtain the adjacent phase difference and the adjacent distance difference after rough difference removal;
the atmosphere delay coefficient updating module is used for updating the atmosphere delay coefficient by using a least square regression method according to the adjacent phase difference and the adjacent distance difference after the coarse difference is removed to obtain the atmosphere delay coefficient after the coarse difference is removed;
and the atmosphere delay correction module is used for calculating the atmosphere delay phase after the gross error correction by using the atmosphere delay coefficient after the gross error removal.
Optionally, the GBInSAR atmospheric correction system based on winding phase further includes:
and the filtering module is used for carrying out Goldstein filtering processing on the interference pattern.
Optionally, the coherency calculation module includes:
the coherence calculation unit is used for calculating the coherence of each pixel in the interference pattern, and the coherence represents the similarity degree of echo phases between two pixel points forming a certain pixel point in the interference pattern through differential interference;
and the high-coherence point screening unit is used for selecting the pixel points with coherence exceeding a preset threshold value as high-coherence points.
Optionally, the filtering module includes:
the block dividing unit is used for dividing the interference pattern into a plurality of interference pattern blocks by using a sliding window, and the overlapping rate of any two adjacent interference pattern blocks is greater than or equal to a preset rate;
the Fourier transform unit is used for carrying out Fourier transform on each interference pattern block to obtain a spectrogram of each interference pattern block;
the smoothing unit is used for smoothing each spectrogram to obtain a smoothed amplitude;
the frequency spectrum processing unit is used for processing the frequency spectrum of each frequency spectrum graph by using the amplitude of the smoothing processing to obtain a processed frequency spectrum graph;
and the inverse Fourier transform unit is used for performing inverse Fourier transform on the processed spectrogram to obtain a filtered interferogram.
Optionally, the GBInSAR atmospheric correction system based on winding phase further includes:
and the data eliminating module is used for eliminating the adjacent distance difference higher than the preset distance difference threshold value and the corresponding adjacent phase difference data.
According to the specific embodiment provided by the invention, the invention discloses the following technical effects: according to the winding phase-based GBInSAR atmospheric correction method and system, after the high coherence point is calculated, the atmospheric delay coefficient is calculated by using a least square regression method according to the adjacent phase difference and the adjacent distance difference, the gross error is removed according to the atmospheric delay phase and the adjacent phase difference, the atmospheric delay coefficient is updated, and the atmospheric delay coefficient after the gross error is removed is obtained, so that atmospheric delay correction is realized under the condition of not performing phase unwrapping, dependence on the phase unwrapping is avoided, the influence of the phase unwrapping on the atmospheric delay correction is avoided, and the accuracy of the atmospheric delay correction is improved.
Drawings
In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings needed to be used in the embodiments will be briefly described below, and it is obvious that the drawings in the following description are only some embodiments of the present invention, and it is obvious for those skilled in the art to obtain other drawings without inventive exercise.
Fig. 1 is a flowchart of a method of GBInSAR atmospheric correction method based on winding phase in embodiment 1 of the present invention;
FIG. 2 is a difference interference diagram of SAR in the Marantan village mine area foundation in the embodiment of the invention 1;
FIG. 3 is a Malan village mine area ground-based SAR differential interferogram after Goldstein filtering in the embodiment of example 1 of the present invention;
FIG. 4 is a simulated wrapped atmospheric delay phase diagram in an embodiment of example 1 of the present invention;
FIG. 5 is a differential interferogram after atmospheric delay correction in the embodiment of example 1 of the present invention;
FIG. 6 is a comparison graph of the phase of the high coherence point in a certain direction of the Malandkian field differential interferogram, the phase after filtering and the simulated atmospheric phase in the embodiment of example 1 of the present invention.
Detailed Description
The technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the drawings in the embodiments of the present invention, and it is obvious that the described embodiments are only a part of the embodiments of the present invention, and not all of the embodiments. All other embodiments, which can be derived by a person skilled in the art from the embodiments given herein without making any creative effort, shall fall within the protection scope of the present invention.
The invention aims to provide a winding phase-based GBInSAR atmospheric correction method and a winding phase-based GBInSAR atmospheric correction system, which are free from dependence on phase unwrapping, thereby avoiding the influence of the phase unwrapping on atmospheric delay correction and improving the accuracy of the atmospheric delay correction.
In order to make the aforementioned objects, features and advantages of the present invention comprehensible, embodiments accompanied with figures are described in further detail below.
Example 1:
fig. 1 is a flowchart of a method of GBInSAR atmospheric correction method based on winding phase in embodiment 1 of the present invention.
Referring to fig. 1, the GBInSAR atmospheric correction method based on winding phase is applied to a ground-based SAR image. The method comprises the following steps:
step 101: acquiring two single-view complex images obtained by detecting the same area by the ground SAR at different moments and an interferogram formed by differential interference of the two single-view complex images. The interference image and the two single-view complex images are all directed at the same target, and pixel points at the same position in the three images are at the same position of the target.
Step 102: and performing Goldstein filtering processing on the interferogram.
The step 102 specifically includes:
A. dividing the interference pattern into a plurality of interference pattern blocks by using a sliding window, wherein the overlapping rate of any two adjacent interference pattern blocks is greater than or equal to a preset rate; preferably, the predetermined ratio is 75%.
B. Fourier transform is carried out on each interference pattern block to obtain a spectrogram of each interference pattern block; the spectrum of the spectrogram is Z (u, v).
C. Smoothing each spectrogram to obtain a smoothed amplitude S { Z (u, v) };
D. and processing the frequency spectrum of each spectrogram by using the amplitude of the smoothing treatment to obtain the processed spectrogram. The processing formula is as follows:
H(u,v)=S{|Z(u,v)|}α·Z(u,v)
where u and v are spatial frequencies. Z (u, v) is a frequency spectrum after Fourier transform, S { Z (u, v) } is an amplitude obtained after smoothing processing, H (u, v) is a frequency spectrum obtained after Z (u, v) is processed by using the amplitude, and alpha is a smoothing degree.
E. And performing inverse Fourier transform on the processed spectrogram to obtain a filtered interferogram.
Step 103: and selecting high coherence points in the interference pattern according to the coherence of each pixel point in the interference pattern. The step 103 specifically includes:
A. the coherence of each pixel in the interferogram is calculated. The coherence represents the degree of similarity of echo phases between two pixel points forming a certain pixel point in the interferogram through differential interference.
The coherence calculation formula is:
Figure BDA0002235502640000071
wherein the content of the first and second substances,
Figure BDA0002235502640000073
for the coherence, M and n respectively represent the sizes of two sides of the sliding window, and M and S respectively represent two monoscopic complex images. i and j respectively represent two coordinates of the pixel points in the interference image block, then M (i, j) represents the pixel point with the coordinate (i, j) in the single-view complex image M, S (i, j) represents the pixel point with the coordinate (i, j) in the single-view complex image S, S*Is a matrix transposition of S, S*(i, j) represents S*And (5) the pixel point with the middle coordinate of (i, j).
B. And selecting the pixel points with the coherence exceeding a preset threshold value as high coherence points.
Usually, the predetermined threshold is a value between 0.9 and 0.95.
Step 104: and acquiring the actually measured phase difference value and distance difference value between the adjacent high-coherence points in the same azimuth direction to obtain the adjacent phase difference and adjacent distance difference. The interferogram includes two directions. The direction along the track is generally considered to be the azimuth direction, and the direction along which the radar waves are emitted is the range direction.
According to the electromagnetic wave propagation theory, the electromagnetic wave with the wavelength lambda is emitted at the time t, and the electromagnetic wave reaches the target point pixel (i, j) from the emission point through the propagation distance r (i, j) and returns, so that the echo phase in the single-view complex image can be expressed as the following formula:
Figure BDA0002235502640000072
wherein phit(i, j) is a single-view complex image acquired at time t, and n (r, t) is an atmospheric refractive index.
For the ground-based SAR, the spatial distribution of the atmospheric medium can be considered to be uniform due to the small observation area, and the phase related to the atmosphere in the single-view complex image can be simply expressed as
In the same way, another single-view complex image can be obtained. When two single-view complex images are subjected to differential interference, the atmospheric phase in the obtained differential interference pattern can be represented as:
Figure BDA0002235502640000082
wherein n (t)2) And n (t)1) Respectively representing the atmospheric refractive index of two single-view complex images.
Is represented by aThen
Figure BDA0002235502640000084
In the conventional method, a is solved by selecting several control points in a way of establishing a function model of distance and interference phase, but the interference pattern needs to be unwrapped due to the existence of k. However, in the embodiment of the present invention, a is obtained by making a difference between adjacent high coherence points.
The atmospheric phase difference between two adjacent high coherence points on the interferogram can be expressed as
Wherein
Figure BDA0002235502640000086
The measured data is the atmospheric phase difference between two adjacent high coherence points, and the measured data is the distance difference between two adjacent high coherence points. This formula
Figure BDA0002235502640000087
A theoretical relationship between measured data of the atmospheric phase difference and measured data of the distance difference is expressed.
The present invention considers that the phase difference is mainly the atmospheric phase difference for two points having no distortion on the interference pattern for a long time, but such an assumption is difficult to satisfy, and a part of the phase difference deviates from the atmospheric phase difference due to the presence of distortion and noise. But when the distance difference is small enough to guarantee the formula
Figure BDA0002235502640000088
When k in (b) is 0, the phase difference and distance difference of the high coherence points in all interferograms can be calculated, and a can be solved by the least square method.
Step 105: and eliminating the adjacent distance difference higher than the preset distance difference threshold value and the corresponding adjacent phase difference data.
The preset distance difference threshold is set according to the degree of the interference pattern influenced by the atmosphere, and when the interference pattern is seriously influenced by the atmosphere, the preset distance difference threshold is set to be smaller, so that the phase jump of most phase differences is avoided.
The distance threshold may also be calculated by the following formula:
Figure BDA0002235502640000089
λ is the electromagnetic wave wavelength, n (t)2) And n (t)1) Respectively representing the atmospheric refractive index of two single-view complex images. R is the distance threshold.
Step 106: and calculating the atmospheric delay coefficient by using a least square regression method according to the adjacent phase difference and the adjacent distance difference.
Considering the phase and distance differences of all the neighboring high coherence points in the interferogram, the following matrix equation can be obtained:
Δψ=ΔR·a+ε
wherein ε ═ ε12,...εn),ΔR=(Δr1,Δr2,...Δrn),
Figure BDA0002235502640000091
(Δr1,Δr2,...Δrn) Respectively the distance difference from the 1 st group of adjacent high coherence points to the nth group of adjacent high coherence points,
Figure BDA0002235502640000092
phase differences from the 1 st group of adjacent high coherence points to the nth group of adjacent high coherence points, n is the total number of adjacent high coherence points, epsilon is a random error, epsilon12,...εn) Random errors corresponding to the phase difference from the 1 st group of adjacent high coherence points to the phase difference from the nth group of adjacent high coherence points.
Each azimuth contains a plurality of sets of adjacent high coherence points. The matrix is formed by sequentially arranging all the adjacent high-coherence points of a plurality of azimuth groups.
The unknown parameter a can be obtained by the following least squares regression:
a=(ΔR*ΔR)-1ΔR*Δψ
ΔR*is a matrix transpose of Δ R.
Step 107: and calculating the atmospheric delay phase according to the atmospheric delay coefficient.
The expression for the atmospheric retardation phase is:
Figure BDA0002235502640000093
Figure BDA0002235502640000094
the atmospheric retardation phase is a theoretical value of the phase difference.
Step 108: and eliminating the adjacent phase difference and the adjacent distance difference corresponding to the difference value exceeding a preset threshold according to the difference value of the atmospheric delay phase and the adjacent phase difference to obtain the adjacent phase difference and the adjacent distance difference after the rough difference is eliminated.
Since the above process does not phase-unwrapp the interference phase, there may be a full-cycle jump in Δ ψ; when the selected high-coherence point is located in the morphable region, then a morphable phase exists in addition to the atmospheric retardation phase in Δ ψ, and therefore, the above formula a is currently used ═ Δ R*ΔR)-1ΔR*The atmospheric delay factor calculated by Δ ψ is inaccurate. For a specific working area, the deformation area is a small area, and the phase jump possibility of the phase difference at a short distance is small. According to the majority theorem, the phase difference of the deformation area and the phase difference with jumping can be regarded as gross error, and the method provided by the invention carries out elimination according to the set standard deviation threshold value. The standard deviation threshold is set by the following formula:
Figure BDA0002235502640000101
where σ is a preset threshold, i.e., a standard deviation threshold.
And after the preset threshold value is obtained through calculation, coarse difference screening is carried out by using the following formula.
Figure BDA0002235502640000102
I.e. not satisfying the formula
Figure BDA0002235502640000103
And (4) eliminating the phase difference data and the distance difference data.
Step 109: and updating the atmospheric delay coefficient by using a least square regression method according to the adjacent phase difference and the adjacent distance difference after the rough difference is removed to obtain the atmospheric delay coefficient after the rough difference is removed.
Step 110: and calculating the atmosphere delay phase after the gross error correction by using the atmosphere delay coefficient after the gross error removal.
The atmospheric retardation phase after the gross error correction is calculated by the following equation:
Figure BDA0002235502640000104
whereinAnd a' is the atmospheric delay coefficient after coarse difference removal.
The interference image after atmospheric phase correction is shown as
Figure BDA0002235502640000106
Wherein the content of the first and second substances,for the corrected phase, wrap represents the winding phase.
Figure BDA0002235502640000108
Is the measured phase
The following describes the GBInSAR atmospheric correction method based on winding phase according to the present invention, with a concrete embodiment of a ground-based SAR interferogram acquired in the marlange mine area at Ku intervals of 56 hours.
The method comprises the following specific steps:
1. the image is filtered with a Goldstein filter with a window size of 5 x 5 to improve the quality of the interferogram.
2. And calculating coherence according to a coherence formula, and selecting points with coherence larger than 0.95 as high-coherence points to obtain 7724 high-coherence points.
3. Since the interferogram is affected more seriously by the atmosphere, 15m is set as a predetermined distance difference threshold.
4. And (3) along each azimuth direction, calculating the distance difference and the phase difference of adjacent high-coherence points, and rejecting data with the distance difference larger than 200 m. And calculating by using a least square method according to the distance difference and the phase difference to obtain the atmospheric delay coefficient.
5. And simulating the atmospheric phase according to the atmospheric delay coefficient, and performing atmospheric delay correction on the interference pattern.
Fig. 2 is a saran area ground-based SAR differential interference map in the embodiment of example 1 of the present invention.
As can be seen from fig. 2, the image front half phase value becomes larger as the distance becomes larger, typically due to the atmospheric phase.
Fig. 3 is a marander zone ground-based SAR differential interferogram after Goldstein filtering in the embodiment of example 1 of the present invention.
As can be seen from fig. 3, the filtered image is less affected by noise, and the phase quality is improved.
FIG. 4 is a simulated wrapped atmospheric delay phase diagram in an embodiment of example 1 of the present invention.
As can be seen from FIG. 4, the simulated atmospheric phase is more consistent with the phase jump of the differential interferogram in FIG. 2, which indicates that the simulation result is more correct.
Fig. 5 is a differential interferogram after atmospheric delay correction in the embodiment of example 1 of the present invention.
As can be seen from fig. 5, the interference pattern after the atmospheric delay correction is removed, the atmospheric delay removal effect is obvious, and the deformation region can be well identified.
FIG. 6 is a comparison graph of the phase of the high coherence point in a certain direction of the Malandkian field differential interferogram, the phase after filtering and the simulated atmospheric phase in the embodiment of example 1 of the present invention.
As can be seen from fig. 6, the simulated atmospheric phase coincides with the phase of the interferogram.
Example 2:
this embodiment 2 provides a GBInSAR atmospheric correction system based on a winding phase, including:
the image acquisition module is used for acquiring two single-view complex images obtained by detecting the same area by the ground SAR at different moments and an interference pattern formed by differential interference of the two single-view complex images;
the coherence calculation module is used for selecting high coherence points in the interference pattern according to the coherence of each pixel point in the interference pattern;
the actual measurement difference value acquisition module is used for acquiring actual measurement phase difference values and distance difference values between adjacent high-coherence points in the same azimuth direction to obtain adjacent phase differences and adjacent distance differences;
the atmosphere delay coefficient calculation module is used for calculating an atmosphere delay coefficient by utilizing a least square regression method according to the adjacent phase difference and the adjacent distance difference;
the atmosphere delay phase calculation module is used for calculating an atmosphere delay phase according to the atmosphere delay coefficient;
the rough difference removing module is used for removing the adjacent phase difference and the adjacent distance difference corresponding to the difference value exceeding a preset threshold according to the difference value of the atmospheric delay phase and the adjacent phase difference to obtain the adjacent phase difference and the adjacent distance difference after rough difference removal;
the atmosphere delay coefficient updating module is used for updating the atmosphere delay coefficient by using a least square regression method according to the adjacent phase difference and the adjacent distance difference after the coarse difference is removed to obtain the atmosphere delay coefficient after the coarse difference is removed;
and the atmosphere delay correction module is used for calculating the atmosphere delay phase after the gross error correction by using the atmosphere delay coefficient after the gross error removal.
Optionally, the GBInSAR atmospheric correction system based on winding phase further includes:
and the filtering module is used for carrying out Goldstein filtering processing on the interference pattern.
Optionally, the coherency calculation module includes:
the coherence calculation unit is used for calculating the coherence of each pixel in the interference pattern, and the coherence represents the similarity degree of echo phases between two pixel points forming a certain pixel point in the interference pattern through differential interference;
and the high coherence screening unit is used for selecting the pixel points with coherence exceeding a preset threshold value as high coherence points.
Optionally, the filtering module includes:
the block dividing unit is used for dividing the interference pattern into a plurality of interference pattern blocks by using a sliding window, and the overlapping rate of any two adjacent interference pattern blocks is greater than or equal to a preset rate;
the Fourier transform unit is used for carrying out Fourier transform on each interference pattern block to obtain a spectrogram of each interference pattern block;
the smoothing unit is used for smoothing each spectrogram to obtain a smoothed amplitude;
the frequency spectrum processing unit is used for processing the frequency spectrum of each frequency spectrum graph by using the amplitude of the smoothing processing to obtain a processed frequency spectrum graph;
and the inverse Fourier transform unit is used for performing inverse Fourier transform on the processed spectrogram to obtain a filtered interferogram.
Optionally, the GBInSAR atmospheric correction system based on winding phase further includes:
and the data eliminating module is used for eliminating the adjacent distance difference higher than the preset distance difference threshold value and the corresponding adjacent phase difference data.
According to the specific embodiment provided by the invention, the invention discloses the following technical effects: according to the winding phase-based GBInSAR atmospheric correction method and system, after the high coherence point is calculated, the atmospheric delay coefficient is calculated by using a least square regression method according to the adjacent phase difference and the adjacent distance difference, the gross error is removed according to the atmospheric delay phase and the adjacent phase difference, the atmospheric delay coefficient is updated, and the atmospheric delay coefficient after the gross error is removed is obtained, so that atmospheric delay correction is realized under the condition of not performing phase unwrapping, dependence on the phase unwrapping is avoided, the influence of the phase unwrapping on the atmospheric delay correction is avoided, and the accuracy of the atmospheric delay correction is improved.
The embodiments in the present description are described in a progressive manner, each embodiment focuses on differences from other embodiments, and the same and similar parts among the embodiments are referred to each other. For the system disclosed by the embodiment, the description is relatively simple because the system corresponds to the method disclosed by the embodiment, and the relevant points can be referred to the method part for description.
The principles and embodiments of the present invention have been described herein using specific examples, which are provided only to help understand the method and the core concept of the present invention; meanwhile, for a person skilled in the art, according to the idea of the present invention, the specific embodiments and the application range may be changed. In view of the above, the present disclosure should not be construed as limiting the invention.

Claims (10)

1. A GBInSAR atmospheric correction method based on winding phase is characterized by comprising the following steps:
acquiring two single-view complex images obtained by detecting the same area at different moments and an interference pattern formed by differential interference of the two single-view complex images;
selecting high coherence points in the interference pattern according to coherence of each pixel point in the interference pattern;
acquiring actually measured phase difference values and distance difference values between adjacent high-coherence points in the same azimuth direction to obtain adjacent phase differences and adjacent distance differences;
calculating an atmospheric delay coefficient by using a least square regression method according to the adjacent phase difference and the adjacent distance difference;
calculating an atmospheric delay phase according to the atmospheric delay coefficient;
according to the difference value of the atmospheric delay phase and the adjacent phase difference, eliminating the adjacent phase difference and the adjacent distance difference corresponding to the difference value exceeding a preset threshold value to obtain the adjacent phase difference and the adjacent distance difference after rough difference removal;
updating the atmospheric delay coefficient by using a least square regression method according to the adjacent phase difference and the adjacent distance difference after the rough difference is removed to obtain the atmospheric delay coefficient after the rough difference is removed;
and calculating the atmosphere delay phase after the gross error correction by using the atmosphere delay coefficient after the gross error removal.
2. The GBInSAR atmospheric correction method according to claim 1, wherein after obtaining two single-view complex images detected from the same area at different time and an interferogram formed by differential interference of the two single-view complex images, before selecting a high coherence point in the interferogram according to coherence of each pixel point in the interferogram, the method further comprises:
and performing Goldstein filtering processing on the interferogram.
3. The winding phase-based GBInSAR atmosphere correction method according to claim 1, wherein the selecting the high coherence point in the interferogram according to the coherence of each pixel point in the interferogram specifically comprises:
calculating the coherence of each pixel in the interference pattern, wherein the coherence represents the similarity of echo phases between two pixel points forming a certain pixel point in the interference pattern through differential interference;
and selecting the pixel points with the coherence exceeding a preset threshold value as high coherence points.
4. The winding phase-based GBInSAR atmospheric correction method of claim 2, wherein the performing Goldstein filtering processing on the interferogram specifically includes:
dividing the interference pattern into a plurality of interference pattern blocks by using a sliding window, wherein the overlapping rate of any two adjacent interference pattern blocks is greater than or equal to a preset rate;
fourier transform is carried out on each interference pattern block to obtain a spectrogram of each interference pattern block;
smoothing each spectrogram to obtain a smoothed amplitude;
processing the frequency spectrum of each spectrogram by using the smoothed amplitude value to obtain a processed spectrogram;
and performing inverse Fourier transform on the processed spectrogram to obtain a filtered interferogram.
5. The winding phase-based GBInSAR atmospheric correction method of claim 1, wherein after obtaining the measured phase difference value and distance difference value between the adjacent high coherence points in the same azimuth direction to obtain the adjacent phase difference and adjacent distance difference, before calculating the atmospheric delay coefficient by using the least square regression method according to the adjacent phase difference and adjacent distance difference, the method further comprises:
and eliminating the adjacent distance difference higher than the preset distance difference threshold value and the corresponding adjacent phase difference data.
6. A winding phase based GBInSAR atmospheric correction system, comprising:
the image acquisition module is used for acquiring two single-view complex images obtained by detecting the same area by the ground SAR at different moments and an interference image formed by differential interference of the two single-view complex images;
the coherence calculation module is used for selecting high coherence points in the interference pattern according to the coherence of each pixel point in the interference pattern;
the actual measurement difference value acquisition module is used for acquiring actual measurement phase difference values and distance difference values between adjacent high-coherence points in the same azimuth direction to obtain adjacent phase differences and adjacent distance differences;
the atmosphere delay coefficient calculation module is used for calculating an atmosphere delay coefficient by utilizing a least square regression method according to the adjacent phase difference and the adjacent distance difference;
the atmosphere delay phase calculation module is used for calculating an atmosphere delay phase according to the atmosphere delay coefficient;
the rough difference removing module is used for removing the adjacent phase difference and the adjacent distance difference corresponding to the difference value exceeding a preset threshold according to the difference value of the atmospheric delay phase and the adjacent phase difference to obtain the adjacent phase difference and the adjacent distance difference after rough difference removal;
the atmosphere delay coefficient updating module is used for updating the atmosphere delay coefficient by using a least square regression method according to the adjacent phase difference and the adjacent distance difference after the coarse difference is removed to obtain the atmosphere delay coefficient after the coarse difference is removed;
and the atmosphere delay correction module is used for calculating the atmosphere delay phase after the gross error correction by using the atmosphere delay coefficient after the gross error removal.
7. The winding phase based GBInSAR atmospheric correction system of claim 6, further comprising:
and the filtering module is used for carrying out Goldstein filtering processing on the interference pattern.
8. The winding phase based GBInSAR atmospheric correction system of claim 6, wherein the coherence computation module comprises:
the coherence calculation unit is used for calculating the coherence of each pixel in the interference pattern, and the coherence represents the similarity degree of echo phases between two pixel points forming a certain pixel point in the interference pattern through differential interference;
and the high-coherence point screening unit is used for selecting the pixel points with coherence exceeding a preset threshold value as high-coherence points.
9. The winding phase based GBInSAR atmospheric correction system of claim 7, wherein the filtering module comprises:
the block dividing unit is used for dividing the interference pattern into a plurality of interference pattern blocks by using a sliding window, and the overlapping rate of any two adjacent interference pattern blocks is greater than or equal to a preset rate;
the Fourier transform unit is used for carrying out Fourier transform on each interference pattern block to obtain a spectrogram of each interference pattern block;
the smoothing unit is used for smoothing each spectrogram to obtain a smoothed amplitude;
the frequency spectrum processing unit is used for processing the frequency spectrum of each frequency spectrum graph by using the amplitude of the smoothing processing to obtain a processed frequency spectrum graph;
and the inverse Fourier transform unit is used for performing inverse Fourier transform on the processed spectrogram to obtain a filtered interferogram.
10. The winding phase based GBInSAR atmospheric correction system of claim 6, further comprising:
and the data eliminating module is used for eliminating the adjacent distance difference higher than the preset distance difference threshold value and the corresponding adjacent phase difference data.
CN201910982029.7A 2019-10-16 2019-10-16 Winding phase-based GBInSAR atmospheric correction method and system Active CN110658521B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201910982029.7A CN110658521B (en) 2019-10-16 2019-10-16 Winding phase-based GBInSAR atmospheric correction method and system

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201910982029.7A CN110658521B (en) 2019-10-16 2019-10-16 Winding phase-based GBInSAR atmospheric correction method and system

Publications (2)

Publication Number Publication Date
CN110658521A true CN110658521A (en) 2020-01-07
CN110658521B CN110658521B (en) 2021-07-27

Family

ID=69041002

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201910982029.7A Active CN110658521B (en) 2019-10-16 2019-10-16 Winding phase-based GBInSAR atmospheric correction method and system

Country Status (1)

Country Link
CN (1) CN110658521B (en)

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111239736A (en) * 2020-03-19 2020-06-05 中南大学 Single-baseline-based surface elevation correction method, device, equipment and storage medium
CN111650570A (en) * 2020-07-24 2020-09-11 中国地质大学(北京) Three-dimensional atmospheric correction method and system for ground-based interference radar

Citations (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP2413158A1 (en) * 2010-07-26 2012-02-01 Consorci Institut de Geomatica A method for monitoring terrain and man-made feature displacements using ground-based synthetic aperture radar (GBSAR) data
CN104007439A (en) * 2014-05-30 2014-08-27 中国科学院电子学研究所 Interferential circular SAR elevation estimation processing method
CN104407332A (en) * 2014-11-25 2015-03-11 沈阳建筑大学 Correction method for updating DEM (digital elevation model) by aid of ground based SAR (synthetic aperture radar)
CN108459322A (en) * 2018-02-09 2018-08-28 长安大学 A kind of InSAR interference patterns batch filtering and preferred method
CN108627833A (en) * 2018-05-15 2018-10-09 电子科技大学 A kind of atmospheric phase compensation method based on GB-InSAR
CN108627834A (en) * 2018-06-07 2018-10-09 北京城建勘测设计研究院有限责任公司 A kind of subway road structure monitoring method and device based on ground InSAR
CN109782282A (en) * 2019-03-13 2019-05-21 武汉大学 A kind of time series InSAR analysis method of integrated troposphere atmospheric delay correction
CN110007302A (en) * 2019-03-27 2019-07-12 北京空间飞行器总体设计部 A kind of spaceborne double antenna strabismus straight rail interference SAR ocean current speed measurement method
CN110109106A (en) * 2019-04-23 2019-08-09 中国电力科学研究院有限公司 A kind of InSAR interferometric phase unwrapping method in region with a varied topography

Patent Citations (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP2413158A1 (en) * 2010-07-26 2012-02-01 Consorci Institut de Geomatica A method for monitoring terrain and man-made feature displacements using ground-based synthetic aperture radar (GBSAR) data
CN104007439A (en) * 2014-05-30 2014-08-27 中国科学院电子学研究所 Interferential circular SAR elevation estimation processing method
CN104407332A (en) * 2014-11-25 2015-03-11 沈阳建筑大学 Correction method for updating DEM (digital elevation model) by aid of ground based SAR (synthetic aperture radar)
CN108459322A (en) * 2018-02-09 2018-08-28 长安大学 A kind of InSAR interference patterns batch filtering and preferred method
CN108627833A (en) * 2018-05-15 2018-10-09 电子科技大学 A kind of atmospheric phase compensation method based on GB-InSAR
CN108627834A (en) * 2018-06-07 2018-10-09 北京城建勘测设计研究院有限责任公司 A kind of subway road structure monitoring method and device based on ground InSAR
CN109782282A (en) * 2019-03-13 2019-05-21 武汉大学 A kind of time series InSAR analysis method of integrated troposphere atmospheric delay correction
CN110007302A (en) * 2019-03-27 2019-07-12 北京空间飞行器总体设计部 A kind of spaceborne double antenna strabismus straight rail interference SAR ocean current speed measurement method
CN110109106A (en) * 2019-04-23 2019-08-09 中国电力科学研究院有限公司 A kind of InSAR interferometric phase unwrapping method in region with a varied topography

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
HOONYOL LEE ET AL.: "AN EXPERIMENT OF GB-SAR INTERPEROMETRIC MEASUREMENT OF TARGET DISPLACEMENT AND ATMOSPHERIC CORRECTION", 《IGARSS 2008》 *
徐亚明 等: "GB-SAR构建永久散射体网改正气象扰动方法", 《武汉大学学报·信息科学版》 *

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111239736A (en) * 2020-03-19 2020-06-05 中南大学 Single-baseline-based surface elevation correction method, device, equipment and storage medium
CN111239736B (en) * 2020-03-19 2022-02-11 中南大学 Single-baseline-based surface elevation correction method, device, equipment and storage medium
CN111650570A (en) * 2020-07-24 2020-09-11 中国地质大学(北京) Three-dimensional atmospheric correction method and system for ground-based interference radar
CN111650570B (en) * 2020-07-24 2023-01-31 中国地质大学(北京) Three-dimensional atmospheric correction method and system for ground-based interference radar

Also Published As

Publication number Publication date
CN110658521B (en) 2021-07-27

Similar Documents

Publication Publication Date Title
CN106772342B (en) Time sequence differential radar interference method suitable for large-gradient ground surface settlement monitoring
CN103675790B (en) A kind of method improving InSAR technical monitoring Ground Deformation precision based on high accuracy DEM
CN104123464B (en) Method for inversion of ground feature high elevation and number of land subsidence through high resolution InSAR timing sequence analysis
CN109541596B (en) InSAR image processing method and device based on deep learning algorithm
CN111273293B (en) InSAR residual motion error estimation method and device considering terrain fluctuation
CN104111456A (en) Ground surface deformation high-resolution interferometric synthetic aperture radar (InSAR) monitoring method along high speed railway
CN110658521B (en) Winding phase-based GBInSAR atmospheric correction method and system
CN109752715B (en) SAR data total-dispersion body detection method and device
WO2010000870A1 (en) Identification and analysis of persistent scatterers in series of sar images
CN110018476B (en) Time difference baseline set time sequence interference SAR processing method
CN111650570B (en) Three-dimensional atmospheric correction method and system for ground-based interference radar
CN113552565B (en) Phase unwrapping method for SAR data high noise and large gradient change area
CN116338607B (en) Time domain and space domain two-step InSAR troposphere delay correction method
CN111239736B (en) Single-baseline-based surface elevation correction method, device, equipment and storage medium
CN108761458B (en) Morphological refinement-based interference SAR water body digital elevation model correction method
CN112685819A (en) Data post-processing method and system for monitoring dam and landslide deformation GB-SAR
Xiang et al. PS Selection Method for and Application to GB‐SAR Monitoring of Dam Deformation
CN115616575B (en) Interference phase diagram winding method assisted by satellite-borne SAR (synthetic aperture radar) stereo measurement
Haque et al. Time series analysis of subsidence in Dhaka City, Bangladesh using Insar
CN112505686B (en) Elevation error estimation method and device, electronic equipment and storage medium
CN115113202A (en) Interference phase iteration unwrapping method based on two-dimensional Gaussian model
CN113341410A (en) Large-range under-forest terrain estimation method, device, equipment and medium
CN109143188A (en) - 1 data ionosphere corrections method of TOPS sentry
CN116736306B (en) Time sequence radar interference monitoring method based on third high-resolution
Eldhuset High resolution spaceborne INSAR simulation with extended scenes

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