CN110109100B - Multi-baseline least square phase unwrapping method based on quality map weighting - Google Patents

Multi-baseline least square phase unwrapping method based on quality map weighting Download PDF

Info

Publication number
CN110109100B
CN110109100B CN201910271422.5A CN201910271422A CN110109100B CN 110109100 B CN110109100 B CN 110109100B CN 201910271422 A CN201910271422 A CN 201910271422A CN 110109100 B CN110109100 B CN 110109100B
Authority
CN
China
Prior art keywords
phase
winding
gradient
initialized
winding phase
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201910271422.5A
Other languages
Chinese (zh)
Other versions
CN110109100A (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.)
University of Electronic Science and Technology of China
Original Assignee
University of Electronic Science and Technology of China
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 University of Electronic Science and Technology of China filed Critical University of Electronic Science and Technology of China
Priority to CN201910271422.5A priority Critical patent/CN110109100B/en
Publication of CN110109100A publication Critical patent/CN110109100A/en
Application granted granted Critical
Publication of CN110109100B publication Critical patent/CN110109100B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S13/00—Systems using the reflection or reradiation of radio waves, e.g. radar systems; Analogous systems using reflection or reradiation of waves whose nature or wavelength is irrelevant or unspecified
    • G01S13/88—Radar or analogous systems specially adapted for specific applications
    • G01S13/89—Radar or analogous systems specially adapted for specific applications for mapping or imaging
    • G01S13/90—Radar or analogous systems specially adapted for specific applications for mapping or imaging using synthetic aperture techniques, e.g. synthetic aperture radar [SAR] techniques
    • G01S13/9021—SAR image post-processing techniques
    • G01S13/9023—SAR image post-processing techniques combined with interferometric techniques
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S7/00—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00
    • G01S7/02—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S13/00
    • G01S7/41—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S13/00 using analysis of echo signal for target characterisation; Target signature; Target cross-section

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)
  • Magnetic Resonance Imaging Apparatus (AREA)
  • Image Analysis (AREA)

Abstract

The invention discloses a multi-baseline least square phase unwrapping method based on mass diagram weighting, which comprises the steps of firstly calculating a primary phase gradient main value of a multi-baseline InSAR wrapped phase, extending the primary phase gradient main value in a mirror symmetry mode, calculating a second-order phase gradient, judging the quality of the wrapped phase by using a mass diagram, performing phase compensation on a residual error point in the wrapped phase by using the primary phase gradient main value as prior information, performing optimized weighting on a second-order phase gradient of the wrapped phase based on the mass diagram, and finally performing multi-baseline least square phase unwrapping to obtain an unwrapped phase.

Description

Multi-baseline least square phase unwrapping method based on quality map weighting
Technical Field
The invention belongs to the technical field of Radar, and particularly relates to the technical field of phase unwrapping in Interferometric Synthetic Aperture Radar (InSAR) elevation mapping.
Background
The interferometric synthetic aperture radar (InSAR) technology is developed in the early 70 s of the 20 th century, combines a radio interferometric measurement technology with an SAR imaging technology, can perform all-time and all-weather interferometric measurement on large-area terrains, and can extract three-dimensional information of the earth surface by using phase information of an SAR image while acquiring a two-dimensional SAR image. The operating principle of the InSAR is that SAR image pairs of the same scene on the ground are obtained through simultaneous observation of two antennas or different aerial observation of a single antenna, SAR image interference phase information is extracted, and a Digital Elevation Model (DEM) is inverted by combining the geometric relation of an imaging scene, so that high-precision three-dimensional topographic surveying and mapping are realized, and the method is widely applied to the fields of natural disaster monitoring, topographic surveying and mapping, natural resource investigation and the like. The multi-baseline InSAR (Multi-baseline InSAR) technology is an extension of the single-baseline InSAR technology, and can overcome the contradiction that the high measurement precision and the phase unwrapping difficulty of a single-baseline InSAR system cannot be considered to a certain extent. The method adopts a multi-baseline InSAR measurement technology to fuse a plurality of interference images, and can well identify the area with large fluctuation in the actual terrain, so that the research on the multi-baseline InSAR system theory and technology is increasingly important in the field of high-precision three-dimensional information acquisition. See the literature "data processing technology of airborne interferometric synthetic aperture radar", and tombia et al, scientific press.
Phase Unwrapting (PU) refers to a process of recovering a winding Phase to a true Phase corresponding to a terrain elevation, and is a crucial step in an InSAR data processing flow. The two SAR images are subjected to conjugate multiplication to obtain interference phases, the interference phases are wound in an interval of [ -pi, pi), the difference between the interference phases and the real phases is integral multiples of 2 pi, phase unwrapping is a main error source for obtaining digital elevation information, and the accuracy of elevation information extraction is directly related. The least square phase unwrapping algorithm based on Fast Fourier Transform (FFT) is a global optimal algorithm, can obtain an optimal solution in the least square sense, is a common phase unwrapping method in practical application at present, and has the basic idea that the phase main value of the phase gradient of the adjacent pixel points of the real phase of the original scene is the same as the phase main value of the phase gradient of the adjacent pixel points of the wrapped phase obtained from echo data. See in detail The literature "Zhang Z, Guo Z. study on The unweighted least-square phase unwarping algorithm [ J ]. Proceedings of SPIE-The International Society for Optical Engineering,2010,7848(1): 30-36.".
The method aims at solving the problems of overlapping and shading and the like existing in a single-baseline InSAR system and low phase unwrapping precision caused by undulating terrain, and needs to improve the traditional phase unwrapping algorithm.
Disclosure of Invention
The invention discloses a multi-baseline least square phase unwrapping method based on mass diagram weighting, which comprises the steps of firstly calculating a primary phase gradient main value of a multi-baseline InSAR wrapped phase, extending the primary phase gradient main value in a mirror symmetry mode, calculating a second-order phase gradient, judging the quality of the wrapped phase by using a mass diagram, performing phase compensation on a residual error point in the wrapped phase by using the primary phase gradient main value as prior information, performing optimized weighting on a second-order phase gradient of the wrapped phase based on the mass diagram, and finally performing multi-baseline least square phase unwrapping to obtain an unwrapped phase.
For the convenience of describing the present invention, the following terms are first defined:
definitions 1, interferometric synthetic Aperture Radar (InSAR)
The interferometric synthetic aperture radar is a synthetic aperture radar technology which utilizes two or more SAR images obtained at different observation angles in the same observation scene to perform interferometric imaging processing, extracts interferometric phase information, and then inverts terrain height and elevation change information by combining radar system parameters, radar platform geometric position parameters and observation terrain information. See the literature "synthetic aperture radar imaging principle", edited by buzz, electronic technology university press.
Definition 2, Baseline
The base line refers to the distance between two antennas in the InSAR system, and the length of the base line of the InSAR system is marked as B in the invention. See the literature "synthetic aperture radar imaging principle", edited by buzz, electronic technology university press.
Definition 3, winding phase:
the winding phase refers to that in actual interference processing, the interference phase is wound in an interval of [ -pi, pi) and is different from the real phase by integral multiple of 2 pi. For details, see literature, "spaceborne synthetic aperture radar interferometry", edited by wangtao et al, scientific press.
Definition 4, phase unwrapping
The process of recovering the winding phase to the true phase corresponding to the terrain elevation is referred to as phase unwrapping. See the literature "data processing technology of airborne interferometric synthetic aperture radar", and tombia et al, scientific press.
Definition 5, unwrapping phase:
the unwrapping phase refers to a true phase recovered from the wrapping phase after the phase unwrapping process in the interference process. For details, see literature, "spaceborne synthetic aperture radar interferometry", edited by wangtao et al, scientific press.
Definition 6, fast Fourier transform and inverse fast Fourier transform
The basic idea of Fast Fourier Transform (FFT) is to decompose an original N-point sequence into a series of short sequences, and to achieve the purposes of eliminating repeated computation, reducing multiplication and simplifying the structure by properly combining the discrete Fourier transforms corresponding to the short sequences by fully utilizing the symmetry and periodicity of exponential factors in the discrete Fourier Transform. The Inverse of the Fast Fourier Transform is denoted as Inverse Fast Fourier Transform (IFFT). See the literature "digital signal processing course", Chen peiyu, Qinghua university Press.
Definition 7, quality map
The quality map is an evaluation model of the phase data quality and represents the quality of the corresponding phase of each pixel in the phase map. The mass map is generally divided into three categories: a coherence coefficient map, a pseudo coherence coefficient map, and a phase derivative variation map. The coherence coefficient map is derived from radar complex images and shows coherence of different areas; the pseudo-coherence coefficient map is that under the condition that real information of two complex images cannot be obtained, the amplitude is assumed to be consistent with the phase information change of the two complex images; the phase derivative change map characterizes the degree of 'bad' phase data quality, a region with rapidly changed phase is defined as a low-quality region, and the phase derivative change map can be used as a first choice of the quality map when the terrain fluctuation is large. For details, the literature "research on unwrapping algorithm for optimizing the measurement of the end surface shape of the optical fiber" is the master academic thesis of Zhejiang university, rufeng.
Definition 8, residual error point
Residual points, also called singular value points, refer to the fact that phase gradients are wound around a current pixel point 2 x 2 window to perform loop integration, a point with an integration result of 2 pi is defined as a positive residual point, a point with an integration result of-2 pi is defined as a negative residual point, and factors such as phase noise pollution and severe topographic relief can cause the existence of residual points. See the literature "data processing technology of airborne interferometric synthetic aperture radar", and tombia et al, scientific press.
Define 9, Main value function
For trigonometric functions, the inverse function is a multi-valued function, and in order to express a determined single-valued function relationship, the multi-valued inverse trigonometric function is restricted, and a specific value range is taken to meet the single-valued requirement in the function definition, and the restricted relationship is called the main value function of the multi-valued inverse trigonometric function and is denoted as W [. See the literature "mathematical guidelines: a practical mathematics handbook ", editions of ebohade, zeitler, etc., scientific press.
Definition 10, mirror symmetry continuation method
The mirror symmetry continuation method is a method for realizing periodic continuation by turning over an image in a mirror symmetry manner. For details, the document "digital image processing", the editions of Okauri et al, electronic industry Press.
Definition 11, normalization
Normalization, also known as dispersion normalization, is a linear transformation of the raw data, mapping the result to [0, 1%]Interval, transformation function is expressed as
Figure GDA0003526491410000031
Wherein xmaxIs the maximum value of the sample data, xminIs the minimum value of the sample data. See the literature "mathematical guidelines: a practical mathematics handbook ", editions of ebohade, zeitler, etc., scientific press.
The invention provides a multi-baseline least square phase unwrapping method based on quality map weighting, which comprises the following steps:
step 1, initializing parameters required by a multi-baseline InSAR least square phase unwrapping method based on mass diagram weighting:
initializing parameters required by a multi-baseline InSAR least square phase unwrapping method based on mass diagram weighting, comprising the following steps of: the number of data sets of the winding phase of the multi-baseline InSAR is recorded as L, the distance direction point number of the winding phase is recorded as M, the direction point number of the winding phase is recorded as N, the length of the baseline is recorded as BlAnd the I group winding phase of the multi-baseline InSAR is recorded as phil(a, b), wherein L is 1,2, …, L, a is 1,2, …, M, b is 1,2, …, N, a represents the a-th point of the distance direction, b represents the b-th point of the azimuth direction; the winding phase is FFT transformed into the corresponding pixel coordinate in the frequency domain as (j, k), where j is 1,2, …, M, k is 1,2, …, N, j denotes the j-th point from the direction to the frequency domain, k denotes the k-th point from the direction to the frequency domain, and the FFT is the fast fourier transform of the conventional standard in definition 6;
step 2, calculating a primary phase gradient main value of a winding phase distance direction and a primary phase gradient main value of an azimuth direction:
using the formula Δl,x(a,b)＝W[φl(a+1,b)-φl(a,b)]And calculating to obtain the winding phase phi of the first grouplThe principal value of the first-order phase gradient in the (a, b) direction of distance, denoted as Δl,x(a, b); using the formula Δl,y(a,b)＝W[φl(a,b+1)-φl(a,b)]And calculating to obtain the winding phase phi of the first groupl(a, b) principal values of azimuthal first-order phase gradients, denoted as Δl,y(a, b), L is 1,2, …, L, a is 1,2, …, M, b is 1,2, …, N, where x denotes a distance direction of the winding phase, y denotes an azimuth direction of the winding phase, L denotes a step1, M represents the number of winding phase distance direction points initialized in step 1, N represents the number of winding phase direction points initialized in step 1, W [ ·]To define the principal value function of the conventional criterion in 9, phil(a, b) is the l-th set of winding phases defined in step 1, a represents the a-th point of the distance direction, b represents the b-th point of the azimuth direction;
step 3, extending a first-order phase gradient main value matrix in mirror symmetry, and calculating a second-order phase gradient:
using a formula
Figure GDA0003526491410000041
For the first group winding phase philPrincipal value of gradient of (a, b) distance to first order phasel,x(a, b) carrying out mirror symmetry continuation to obtain a principal value of the distance to the first-order phase gradient after the continuation, and recording the principal value as rhol,x(a, b); using a formula
Figure GDA0003526491410000042
For the first group winding phase phil(a, b) principal value of azimuthal first-order phase gradient Δl,y(a, b) carrying out mirror symmetry continuation to obtain a main value of the azimuth first-order phase gradient after continuation, and recording the main value as rhol,y(a, b), L is 1,2, …, L, a is 1,2, …,2M, b is 1,2, …,2N, wherein L represents the number of winding phase data sets initialized in step 1, M represents the number of winding phase distance direction points initialized in step 1, N represents the number of winding phase direction points initialized in step 1, Δl,x(a, b) and. DELTA.l,y(a, b) are respectively the first group winding phase phi calculated in step 2l(a, b) first-order phase gradient main values in the distance direction and the azimuth direction, wherein the mirror symmetry continuation is a mirror symmetry continuation method of the traditional standard defined in 10;
using the formula rhol(a,b)＝[ρl,x(a+1,b)-ρl,x(a,b)]+[ρl,y(a,b+1)-ρl,y(a,b)]And calculating to obtain the winding phase phi of the first grouplSecond order phase gradient of (a, b), denoted as ρl(a, b), L is 1,2, …, L, a is 1,2, …,2M, b is 1,2, …,2N, wherein L represents the number of winding phase data sets initialized in step 1,m represents the number of winding phase distance direction points initialized in step 1, N represents the number of winding phase direction points initialized in step 1, ρl,x(a,b)、ρl,y(a, b) respectively representing the first-order phase gradient main values of the extended winding phase distance direction and the extended azimuth direction;
step 4, identifying residual error points in the winding phase:
using a formula
Figure GDA0003526491410000051
Calculating to obtain a loop integral of the first-order gradient main value of the pixel point (a, b) in the a-th distance direction and the b-th azimuth direction in the winding phase in the 2 multiplied by 2 neighborhood, and marking as R (a, b), wherein rhoi(i is 1,2,3,4) is a winding phase first-order gradient principal value, ρ1＝ρl,x(a,b)、ρ3＝-ρl,x(a, b +1) is the principal value of the first-order phase gradient of the distance direction after continuation in step 3, ρ2＝ρl,y(a+1,b)、ρ4＝-ρl,y(a, b) is a main value of the first-order phase gradient of the azimuth direction after continuation in step 3, L is 1,2, …, L, a is 1,2, …,2M, b is 1,2, …,2N, L represents the number of winding phase data sets initialized in step 1, M represents the number of winding phase distance direction points initialized in step 1, and N represents the number of winding phase azimuth direction points initialized in step 1;
the method for identifying residual error points in winding phases comprises the following steps: when the first-order gradient principal value of the winding phase pixel point (a, b) is in the 2 × 2 neighborhood, the loop integration result | R (a, b) | is 2 pi, (a, b) is a residual point, and when R (a, b) | 0, (a, b) is a continuous point;
and 5, performing phase compensation on the residual error points based on the quality map:
using a formula
Figure GDA0003526491410000052
Calculating to obtain a quality diagram of the l group of winding phases, and recording the phase quality of pixel points (a, b) in the a-th distance direction and the b-th azimuth direction as zl(a, b), L is 1,2, …, L, a is 1,2, …,2M, b is 1,2, …,2N, wherein L represents the number of winding phase data sets initialized in step 1, and M represents the winding phase initialized in step 1The number of azimuth distance points, N represents the number of azimuth direction points of the winding phase initialized in step 1, rhol,x(a,b)、ρl,y(a, b) respectively representing the first-order phase gradient main values of the distance direction and the azimuth direction after continuation in the step 3,
Figure GDA0003526491410000061
respectively representing the mean values of first-order phase gradient main values of the distance direction and the azimuth direction of the current pixel point after extension in a mxm (M is more than or equal to 2 and less than or equal to M) window;
using a formula
Figure GDA0003526491410000062
Calculating to obtain a phase compensation factor, which is recorded as ni(i ═ 1,2,3,4), where z isi(i ═ 1,2,3,4) is the phase mass, z1Denotes the phase quality, z, corresponding to the residual point (a, b) in the winding phase2Represents the phase mass, z, corresponding to point (a +1, b)3Represents the phase mass, z, corresponding to point (a +1, b +1)4Representing the phase quality corresponding to the point (a, b +1), wherein R (a, b) is the loop integral of the first-order gradient main value of the winding phase pixel point (a, b) calculated in the step 4 in the neighborhood of 2 multiplied by 2; using a formula
Figure GDA0003526491410000063
Performing residual point phase compensation, wherein a is 1,2, …,2M, b is 1,2, …,2N, and rho isi(i is 1,2,3,4) is the primary value of the first order phase gradient in step 4,
Figure GDA0003526491410000064
the primary value of the first-order phase gradient after phase compensation, M represents the number of winding phase distance direction points initialized in the step 1, and N represents the number of winding phase direction points initialized in the step 1;
and 6, performing weighted optimization on the second-order phase gradient based on the quality map:
winding phase mass z after compensating residual point phase in step 5l(a, b) normalized to [0,1 ] using the dispersion normalization method of the conventional standard in definition 11]Within the interval, the formula Z is adoptedl,x(a,b)＝min(zl(a-1,b),zl(a, b)), calculating a range gradient weight, denoted as Zl,x(a, b); using the formula Zl,y(a,b)＝min(zl(a,b-1),zl(a, b)), calculating to obtain an azimuth weight gradient, and recording as Zl,y(a, b), L1, 2, …, L, a 1,2, …,2M, b 1,2, …,2N, wherein z isl(a, b) is the phase quality of the pixel point of the ith winding phase, the a-th distance direction and the b-th direction of the L group of winding phases calculated in the step 5, L represents the number of winding phase data sets initialized in the step 1, M represents the number of winding phase distance direction points initialized in the step 1, N represents the number of winding phase direction points initialized in the step 1, and a function min (parameter 1, parameter 2) represents the minimum value in a calculation parameter list;
using a formula
Figure GDA0003526491410000065
Calculating to obtain a second-order phase gradient compensation factor of the mth group of winding phases, and recording the second-order phase gradient compensation factor as
Figure GDA0003526491410000066
L1, 2, …, L, M1, 2, …, L, a 1,2, …,2M, b 1,2, …,2N, wherein L represents the number of winding phase data sets initialized in step 1, M represents the number of winding phase distance points initialized in step 1, N represents the number of winding phase direction points initialized in step 1, Z represents the number of winding phase direction points initialized in step 1, and Z represents the number of winding phase direction points initialized in step 1l,x(a, b) and Zl,y(a, B) are respectively the gradient weight of the distance direction and the gradient weight of the azimuth direction, BlDenotes the length of the i-th base line defined in step 1, BmRepresents the length of the mth baseline defined in step 1;
using a formula
Figure GDA0003526491410000071
Calculating to obtain a second-order phase gradient of the m-th group of winding phases after weighting optimization, and recording the second-order phase gradient as
Figure GDA0003526491410000072
Where L denotes the number of winding phase data sets initialized in step 1, and M denotes the number of initialization in step 1,2, …, L, a 1,2, …,2M, b 1,2, …,2NN represents the number of winding phase azimuth points initialized in step 1,
Figure GDA0003526491410000073
second order phase gradient compensation factor, p, for the mth set of winding phasesm(a, b) is the second order phase gradient of the winding phase calculated in step 3;
and 7, performing phase unwrapping by using the weighted and optimized second-order phase gradient:
the mth group of winding phase second-order phase gradients after the weighting optimization in the step 6
Figure GDA0003526491410000074
Calculating to obtain a winding phase second-order phase gradient by adopting a fast Fourier transform method of the traditional standard in definition 6
Figure GDA0003526491410000075
The two-dimensional FFT result of (1) is denoted as P (j, k), and the formula is adopted
Figure GDA0003526491410000076
Calculating to obtain the extended unwrapping phase
Figure GDA0003526491410000077
The two-dimensional FFT result of (1) is represented by Φ (j, k), j is 1,2, …,2M, k is 1,2, …,2N, j represents the jth point in the distance frequency domain, k represents the kth point in the azimuth frequency domain, M represents the number of winding phase distance points initialized in step 1, and N represents the number of winding phase azimuth points initialized in step 1; calculating the extended unwrapping phase by adopting a fast Fourier inverse transformation method of the traditional standard in definition 6 for phi (j, k), and recording the extended unwrapping phase as
Figure GDA0003526491410000078
When a is 1,2, …,2M, b is 1,2, …,2N, take
Figure GDA0003526491410000079
Wherein a is 1,2, …, M, b is 1,2, …, N, to obtain the final unwrapping phase, which is recorded asIs psi (a, b).
The invention has the innovation point that the multi-baseline least square phase unwrapping method based on the mass diagram weighting is provided, the mass diagram is used as prior information to compensate residual points in an wrapped phase, the problem of residual point phase error transmission in the traditional single-baseline least square phase unwrapping algorithm is solved, the wrapping phase second-order phase gradient is optimally weighted based on the mass diagram to reduce unwrapping errors, and the phase unwrapping precision is improved.
The method has the advantages of simple realization, high precision and good applicability, and can better identify the high-range fluctuating area in the real terrain by fusing a plurality of groups of winding phase information through the multi-baseline InSAR technology; compensating residual error points in the winding phase, and avoiding the error transmission problem caused by the phase residual error points in the traditional algorithm; the second-order phase gradient of the winding phase is optimized in a weighting mode, and the phase unwrapping precision is improved.
Drawings
FIG. 1 is a schematic block flow diagram of a method provided by the present invention;
fig. 2 shows a multi-baseline InSAR phase unwrapping simulation parameter of the method provided by the present invention.
Detailed Description
The invention can be verified by adopting a simulation experiment method, and all steps and conclusions are verified to be correct on MATLAB R2017b software. The specific implementation steps are as follows:
step 1, initializing parameters required by a multi-baseline InSAR least square phase unwrapping method based on mass diagram weighting:
initializing parameters required by a multi-baseline InSAR least square phase unwrapping method based on mass diagram weighting, comprising the following steps of: the number of data sets of the winding phase of the multi-baseline InSAR is recorded as L being 5, the number of direction points of the distance of the winding phase is recorded as M being 500, the number of direction points of the winding phase is recorded as N being 500, the length of the baseline is B1＝3m，B2＝7m，B3＝13m，B4＝19m，B531m, the I group winding phase of the multi-base line InSAR is recorded as phil(a, b), wherein l is 1,2, …,5, a is 1,2, …,500, b is 1,2, …,500, a represents the a-th point of the distance direction, b represents the azimuth directionPoint b; the winding phase is transformed into the corresponding pixel coordinate in the frequency domain by FFT and is marked as (j, k), wherein j is 1,2, …,500, k is 1,2, …,500, j represents the j point in the distance frequency domain, and k represents the k point in the azimuth frequency domain;
step 2, calculating a primary phase gradient main value of a winding phase distance direction and a primary phase gradient main value of an azimuth direction:
using the formula Δl,x(a,b)＝W[φl(a+1,b)-φl(a,b)]And calculating to obtain the winding phase phi of the first grouplThe principal value of the first-order phase gradient in the (a, b) direction of distance, denoted as Δl,x(a, b); using the formula Δl,y(a,b)＝W[φl(a,b+1)-φl(a,b)]And calculating to obtain the winding phase phi of the first groupl(a, b) principal values of azimuthal first-order phase gradients, denoted as Δl,y(a, b), L is 1,2, …, L, a is 1,2, …, M, b is 1,2, …, N, L is 5, M is 500, N is 500, where L is the number of winding phase data sets initialized in step 1, M is the number of winding phase distance direction points initialized in step 1, N is the number of winding phase direction points initialized in step 1, x is the distance direction of the winding phase, y is the direction of the winding phase, W [ · is]To define the principal value function of the conventional criterion in 9, phil(a, b) is the l set of winding phases defined in step 1, a represents the a-th point of the distance direction, b represents the b-th point of the azimuth direction;
step 3, extending a first-order phase gradient main value matrix in mirror symmetry, and calculating a second-order phase gradient:
using a formula
Figure GDA0003526491410000081
For the first group winding phase philPrincipal value of gradient of (a, b) distance to first order phasel,x(a, b) carrying out mirror symmetry continuation to obtain a principal value of the distance to the first-order phase gradient after the continuation, and recording the principal value as rhol,x(a, b); using a formula
Figure GDA0003526491410000091
For the first group winding phase phil(a, b) principal value of azimuthal first-order phase gradient Δl,y(a, b) proceeding mirror symmetry continuation to obtain the extended squarePrincipal value of bit direction first order phase gradient, denoted as ρl,y(a, b), L is 1,2, …, L, a is 1,2, …,2M, b is 1,2, …,2N, L is 5, M is 500, N is 500, where L is the number of winding phase data sets initialized in step 1, M is the number of winding phase distance points initialized in step 1, N is the number of winding phase direction points initialized in step 1, Δ is the number of winding phase direction points initialized in step 1, and L is 1,2, …l,x(a, b) and. DELTA.l,y(a, b) are respectively the first group winding phase phi calculated in step 2l(a, b) first-order phase gradient main values in the distance direction and the azimuth direction, wherein the mirror symmetry continuation is a mirror symmetry continuation method of the traditional standard defined in 10;
using the formula rhol(a,b)＝[ρl,x(a+1,b)-ρl,x(a,b)]+[ρl,y(a,b+1)-ρl,y(a,b)]And calculating to obtain the winding phase phi of the first grouplSecond order phase gradient of (a, b), denoted as ρl(a, b), L is 1,2, …, L, a is 1,2, …,2M, b is 1,2, …,2N, L is 5, M is 500, N is 500, where L represents the number of winding phase data sets initialized in step 1, M represents the number of winding phase distance points initialized in step 1, N represents the number of winding phase azimuth points initialized in step 1, ρ is the number of winding phase azimuth points initialized in step 1, andl,x(a,b)、ρl,y(a, b) respectively representing the first-order phase gradient main values of the extended winding phase distance direction and the extended azimuth direction;
step 4, identifying residual error points in the winding phase:
using a formula
Figure GDA0003526491410000092
Calculating to obtain a loop integral of the first-order gradient main value of the pixel point (a, b) in the a-th distance direction and the b-th azimuth direction in the winding phase in the 2 multiplied by 2 neighborhood, and marking as R (a, b), wherein rhoi(i is 1,2,3,4) is a winding phase first-order gradient principal value, ρ1＝ρl,x(a,b)、ρ3＝-ρl,x(a, b +1) is the principal value of the first-order phase gradient of the distance direction after continuation in step 3, ρ2＝ρl,y(a+1,b)、ρ4＝-ρl,y(a, b) is the principal value of the first-order phase gradient in the azimuth direction after continuation in step 3, L is 1,2, …, L, a is 1,2, …,2M, b is 1,2, …,2N, L is 5, M is 500, N is 500,l represents the number of winding phase data sets initialized in step 1, M represents the number of winding phase distance orientation points initialized in step 1, and N represents the number of winding phase orientation points initialized in step 1;
the method for identifying residual error points in winding phases comprises the following steps: when the first-order gradient principal value of the pixel point (a, b) in the winding phase is in the ring integration result | R (a, b) | 2 pi in the 2 × 2 neighborhood, (a, b) is a residual point, and when R (a, b) | 0, (a, b) is a continuous point;
and 5, performing phase compensation on the residual error points based on the quality map:
using a formula
Figure GDA0003526491410000101
Calculating to obtain a quality diagram of the l group of winding phases, and recording the phase quality of pixel points (a, b) in the a-th distance direction and the b-th azimuth direction as zl(a, b), L is 1,2, …, L, a is 1,2, …,2M, b is 1,2, …,2N, L is 5, M is 500, N is 500, where L represents the number of winding phase data sets initialized in step 1, M represents the number of winding phase distance points initialized in step 1, N represents the number of winding phase azimuth points initialized in step 1, ρ is the number of winding phase azimuth points initialized in step 1, andl,x(a,b)、ρl,y(a, b) respectively representing the first-order phase gradient main values of the distance direction and the azimuth direction after continuation in the step 3,
Figure GDA0003526491410000102
respectively representing the mean values of the first-order phase gradient main values of the distance direction and the azimuth direction after the current pixel point is extended in a 2 multiplied by 2 window;
using a formula
Figure GDA0003526491410000103
Calculating to obtain a phase compensation factor, which is recorded as ni(i ═ 1,2,3,4), where z isi(i ═ 1,2,3,4) is the phase mass, z1Denotes the phase quality, z, corresponding to the residual point (a, b) in the winding phase2Represents the phase mass, z, corresponding to point (a +1, b)3Represents the phase mass, z, corresponding to point (a +1, b +1)4Representing the phase quality corresponding to point (a, b +1), R (a, b) being the winding phase calculated in step 4The loop integral of the first-order gradient main value of the pixel point (a, b) in the neighborhood of 2 multiplied by 2; using a formula
Figure GDA0003526491410000104
Performing residual point phase compensation, wherein a is 1,2, …,2M, b is 1,2, …,2N, M is 500, and N is 500, wherein M represents the number of winding phase distance direction points initialized in step 1, N represents the number of winding phase direction points initialized in step 1, and ρ isi(i is 1,2,3,4) is the primary value of the first order phase gradient in step 4,
Figure GDA0003526491410000105
is the primary value of the first-order phase gradient after phase compensation;
step 6, weighting and optimizing the multi-baseline second-order phase gradient based on the quality diagram:
winding phase mass z after compensating residual point phase in step 5l(a, b) normalized to [0,1 ] using the conventional dispersion normalization method in definition 11]Within the interval, the formula Z is adoptedl,x(a,b)＝min(zl(a-1,b),zl(a, b)), calculating a range gradient weight, denoted as Zl,x(a, b); using the formula Zl,y(a,b)＝min(zl(a,b-1),zl(a, b)), calculating to obtain an azimuth weight gradient, and recording as Zl,y(a, b), L is 1,2, …, L, a is 1,2, …,2M, b is 1,2, …,2N, L is 5, M is 500, N is 500, wherein L represents the number of winding phase data sets initialized in step 1, M represents the number of winding phase distance points initialized in step 1, N represents the number of winding phase position points initialized in step 1, z is zl(a, b) is the phase quality of the pixel point in the ith distance direction and the mth azimuth direction of the l group of winding phases calculated in the step 5, and min (parameter 1, parameter 2) represents the minimum value in the calculation parameter list;
using a formula
Figure GDA0003526491410000111
Calculating to obtain a second-order phase gradient compensation factor of the mth group of winding phases, and recording the second-order phase gradient compensation factor as
Figure GDA00035264914100001111
L1, 2, …, L, M1, 2, …, L, a 1,2, …,2M, b 1,2, …,2N, L5, M500, N500, wherein L represents the number of winding phase data sets initialized in step 1, M represents the number of winding phase distance points initialized in step 1, N represents the number of winding phase direction points initialized in step 1, Z represents the number of winding phase direction points initialized in step 1, and L, a 1,2, …, L, a 1, 2M, 2N, L, M, 500, N, 500l,x(a, b) and Zl,y(a, B) are respectively the gradient weight of distance direction and the gradient weight of azimuth direction, and the length of the base line is B1＝3m、B2＝7m、B3＝13m、B4＝19m、B5＝31m；
Using a formula
Figure GDA0003526491410000112
Calculating to obtain a second-order phase gradient of the m-th group of winding phases after weighting optimization, and recording the second-order phase gradient as
Figure GDA0003526491410000113
Where L denotes the number of winding phase data sets initialized in step 1, M denotes the number of winding phase distance direction points initialized in step 1, N denotes the number of winding phase direction points initialized in step 1, M denotes the number of winding phase direction points initialized in step 1,2, …, L, a denotes 1,2, …,2M, b denotes 1,2, …,2N, L denotes 5, M denotes 500, N denotes 500,
Figure GDA0003526491410000114
second order phase gradient compensation factor, p, for the mth set of winding phasesm(a, b) is the second order phase gradient of the winding phase calculated in step 3;
and 7, performing phase unwrapping by using the weighted and optimized second-order phase gradient:
the mth group of winding phase second-order phase gradients after the weighting optimization in the step 6
Figure GDA0003526491410000115
Calculating to obtain a winding phase second-order phase gradient by adopting a fast Fourier transform method of the traditional standard in definition 6
Figure GDA0003526491410000116
The two-dimensional FFT result of (1) is denoted as P (j, k), and the formula is adopted
Figure GDA0003526491410000117
Calculating to obtain the extended unwrapping phase
Figure GDA0003526491410000118
The two-dimensional FFT result of (1) is represented by Φ (j, k), j is 1,2, …,2M, k is 1,2, …,2N, M is 500, N is 500, j represents the jth point in the distance frequency domain, k represents the kth point in the azimuth frequency domain, M represents the number of winding phase distance direction points initialized in step 1, and N represents the number of winding phase azimuth direction points initialized in step 1; calculating the extended unwrapping phase by adopting a fast Fourier inverse transformation method of the traditional standard in definition 6 for phi (j, k), and recording the extended unwrapping phase as
Figure GDA0003526491410000119
When a is 1,2, …,2M, b is 1,2, …,2N, take
Figure GDA00035264914100001110
The final unwrapped phase is obtained when a is 1,2, …, M, b is 1,2, …, N, and is denoted as ψ (a, b).

Claims (1)

1. A multi-baseline least square phase unwrapping method based on quality map weighting is characterized by comprising the following steps:
step 1, initializing parameters required by a multi-baseline InSAR least square phase unwrapping method based on mass diagram weighting:
initializing parameters required by a multi-baseline InSAR least square phase unwrapping method based on mass diagram weighting, comprising the following steps of: the number of data sets of the winding phase of the multi-baseline InSAR is recorded as L, the distance direction point number of the winding phase is recorded as M, the direction point number of the winding phase is recorded as N, the length of the baseline is recorded as BlAnd the I group winding phase of the multi-baseline InSAR is recorded as phil(a, b), wherein L is 1,2, …, L, a is 1,2, …, M, b is 1,2, …, N, a represents the a-th point of the distance direction, b represents the b-th point of the azimuth direction; the winding phase is FFT transformed to the corresponding pixel coordinate in the frequency domain and is recorded as(j, k), wherein j is 1,2, …, M, k is 1,2, …, N, j represents the j-th point of the distance to frequency domain, k represents the k-th point of the azimuth frequency domain, and FFT is the fast fourier transform of the conventional standard;
step 2, calculating a primary phase gradient main value of a winding phase distance direction and a primary phase gradient main value of an azimuth direction:
using the formula Δl,x(a,b)＝W[φl(a+1,b)-φl(a,b)]And calculating to obtain the winding phase phi of the first grouplThe principal value of the first-order phase gradient in the (a, b) direction of distance, denoted as Δl,x(a, b); using the formula Δl,y(a,b)＝W[φl(a,b+1)-φl(a,b)]And calculating to obtain the winding phase phi of the first groupl(a, b) principal values of azimuthal first-order phase gradients, denoted as Δl,y(a, b), L is 1,2, …, L, a is 1,2, …, M, b is 1,2, …, N, where x represents the distance direction of the winding phase, y represents the azimuth direction of the winding phase, L represents the number of winding phase data sets initialized in step 1, M represents the number of winding phase distance direction points initialized in step 1, N represents the number of winding phase azimuth direction points initialized in step 1, W [ · o]Is a function of the principal value of the conventional standard, phil(a, b) is the l set of winding phases initialized in step 1, wherein a represents the a-th point of the distance direction, and b represents the b-th point of the azimuth direction;
step 3, extending a first-order phase gradient main value matrix in mirror symmetry, and calculating a second-order phase gradient:
using a formula
Figure FDA0003526491400000011
For the first group winding phase philPrincipal value of gradient of (a, b) distance to first order phasel,x(a, b) carrying out mirror symmetry continuation to obtain a principal value of the distance to the first-order phase gradient after the continuation, and recording the principal value as rhol,x(a, b); using a formula
Figure FDA0003526491400000012
For the first group winding phase phil(a, b) principal value of azimuthal first-order phase gradient Δl,y(a, b) carrying out mirror symmetry continuation to obtain a main value of the azimuth first-order phase gradient after continuation, and recording the main value as rhol,y(a, b), L is 1,2, …, L, a is 1,2, …,2M, b is 1,2, …,2N, wherein L represents the number of winding phase data sets initialized in step 1, M represents the number of winding phase distance direction points initialized in step 1, N represents the number of winding phase direction points initialized in step 1, Δl,x(a, b) and. DELTA.l,y(a, b) are respectively the first group winding phase phi calculated in step 2l(a, b) the first-order phase gradient main values in the distance direction and the azimuth direction are extended in a mirror symmetry mode to be a conventional standard mirror symmetry extension method;
using the formula rhol(a,b)＝[ρl,x(a+1,b)-ρl,x(a,b)]+[ρl,y(a,b+1)-ρl,y(a,b)]And calculating to obtain the winding phase phi of the first grouplSecond order phase gradient of (a, b), denoted as ρl(a, b), L is 1,2, …, L, a is 1,2, …,2M, b is 1,2, …,2N, wherein L represents the number of winding phase data sets initialized in step 1, M represents the number of winding phase distance direction points initialized in step 1, N represents the number of winding phase direction points initialized in step 1, ρ isl,x(a,b)、ρl,y(a, b) respectively representing the first-order phase gradient main values of the extended winding phase distance direction and the extended azimuth direction;
step 4, identifying residual error points in the winding phase:
using a formula
Figure FDA0003526491400000021
Calculating to obtain a loop integral of the first-order gradient main value of the pixel point (a, b) in the a-th distance direction and the b-th azimuth direction in the winding phase in the 2 multiplied by 2 neighborhood, and marking as R (a, b), wherein rhoi(i is 1,2,3,4) is a winding phase first-order gradient principal value, ρ1＝ρl,x(a,b)、ρ3＝-ρl,x(a, b +1) is the principal value of the first-order phase gradient of the distance direction after continuation in step 3, ρ2＝ρl,y(a+1,b)、ρ4＝-ρl,y(a, b) is the principal value of the first-order phase gradient of the azimuth direction after continuation in step 3, L is 1,2, …, L, a is 1,2, …,2M, b is 1,2, …,2N, L represents the number of winding phase data sets initialized in step 1, M represents the number of winding phase distance direction points initialized in step 1, and N is shown in tableShowing the number of winding phase azimuth points initialized in the step 1;
the method for identifying residual error points in winding phases comprises the following steps: when the first-order gradient principal value of the winding phase pixel point (a, b) is in the 2 × 2 neighborhood, the loop integration result | R (a, b) | is 2 pi, (a, b) is a residual point, and when R (a, b) | 0, (a, b) is a continuous point;
and 5, performing phase compensation on the residual error points based on the quality map:
using a formula
Figure FDA0003526491400000022
Calculating to obtain a quality diagram of the l group of winding phases, and recording the phase quality of pixel points (a, b) in the a-th distance direction and the b-th azimuth direction as zl(a, b), L is 1,2, …, L, a is 1,2, …,2M, b is 1,2, …,2N, wherein L represents the number of winding phase data sets initialized in step 1, M represents the number of winding phase distance direction points initialized in step 1, N represents the number of winding phase direction points initialized in step 1, ρ isl,x(a,b)、ρl,y(a, b) respectively representing the first-order phase gradient main values of the distance direction and the azimuth direction after continuation in the step 3,
Figure FDA0003526491400000031
respectively representing the mean values of first-order phase gradient main values of the distance direction and the azimuth direction of the current pixel point after extension in a mxm (M is more than or equal to 2 and less than or equal to M) window;
using a formula
Figure FDA0003526491400000032
Calculating to obtain a phase compensation factor, which is recorded as ni(i ═ 1,2,3,4), where z isi(i ═ 1,2,3,4) is the phase mass, z1Denotes the phase quality, z, corresponding to the residual point (a, b) in the winding phase2Represents the phase mass, z, corresponding to point (a +1, b)3Represents the phase mass, z, corresponding to point (a +1, b +1)4Representing the phase quality corresponding to the point (a, b +1), wherein R (a, b) is the loop integral of the first-order gradient main value of the winding phase pixel point (a, b) calculated in the step 4 in the neighborhood of 2 multiplied by 2; using a formula
Figure FDA0003526491400000033
Performing residual point phase compensation, wherein a is 1,2, …,2M, b is 1,2, …,2N, and rho isi(i is 1,2,3,4) is the primary value of the first order phase gradient in step 4,
Figure FDA0003526491400000034
the primary value of the first-order phase gradient after phase compensation, M represents the number of winding phase distance direction points initialized in the step 1, and N represents the number of winding phase direction points initialized in the step 1;
and 6, performing weighted optimization on the second-order phase gradient based on the quality map:
winding phase mass z after compensating residual point phase in step 5l(a, b) normalized to [0,1 ] using the dispersion normalization method of the conventional standard]Within the interval, the formula Z is adoptedl,x(a,b)＝min(zl(a-1,b),zl(a, b)), calculating a range gradient weight, denoted as Zl,x(a, b); using the formula Zl,y(a,b)＝min(zl(a,b-1),zl(a, b)), calculating to obtain an azimuth weight gradient, and recording as Zl,y(a, b), L1, 2, …, L, a 1,2, …,2M, b 1,2, …,2N, wherein z isl(a, b) is the phase quality of the pixel point of the ith winding phase, the a-th distance direction and the b-th direction of the L group of winding phases calculated in the step 5, L represents the number of winding phase data sets initialized in the step 1, M represents the number of winding phase distance direction points initialized in the step 1, N represents the number of winding phase direction points initialized in the step 1, and a function min (parameter 1, parameter 2) represents the minimum value in a calculation parameter list;
using a formula
Figure FDA0003526491400000035
Calculating to obtain a second-order phase gradient compensation factor of the mth group of winding phases, and recording the second-order phase gradient compensation factor as
Figure FDA0003526491400000036
1,2, …, L, M1, 2, …, L, a 1,2, …,2M, b 1,2, …,2N, and combinations thereofWhere L represents the number of winding phase data sets initialized in step 1, M represents the number of winding phase distance direction points initialized in step 1, N represents the number of winding phase direction points initialized in step 1, and Z represents the number of winding phase direction points initialized in step 1l,x(a, b) and Zl,y(a, B) are respectively the gradient weight of the distance direction and the gradient weight of the azimuth direction, BlDenotes the length of the first baseline initialized in step 1, BmRepresents the length of the mth baseline initialized in step 1;
using a formula
Figure FDA0003526491400000041
Calculating to obtain a second-order phase gradient of the m-th group of winding phases after weighting optimization, and recording the second-order phase gradient as
Figure FDA0003526491400000042
M is 1,2, …, L, a is 1,2, …,2M, b is 1,2, …,2N, wherein L represents the number of winding phase data sets initialized in step 1, M represents the number of winding phase distance points initialized in step 1, N represents the number of winding phase direction points initialized in step 1,
Figure FDA0003526491400000043
second order phase gradient compensation factor, p, for the mth set of winding phasesm(a, b) is the second order phase gradient of the winding phase calculated in step 3;
and 7, performing phase unwrapping by using the weighted and optimized second-order phase gradient:
the mth group of winding phase second-order phase gradients after the weighting optimization in the step 6
Figure FDA0003526491400000044
Calculating to obtain a second-order phase gradient of a winding phase by adopting a traditional standard fast Fourier transform method
Figure FDA0003526491400000045
The two-dimensional FFT result of (1) is denoted as P (j, k), and the formula is adopted
Figure FDA0003526491400000046
Calculating to obtain the extended unwrapping phase
Figure FDA0003526491400000047
The two-dimensional FFT result of (1) is represented by Φ (j, k), j is 1,2, …,2M, k is 1,2, …,2N, j represents the jth point in the distance frequency domain, k represents the kth point in the azimuth frequency domain, M represents the number of winding phase distance points initialized in step 1, and N represents the number of winding phase azimuth points initialized in step 1; calculating the extended unwrapping phase by adopting the traditional standard fast Fourier inverse transformation method for phi (j, k), and recording the extended unwrapping phase as
Figure FDA0003526491400000048
When a is 1,2, …,2M, b is 1,2, …,2N, take
Figure FDA0003526491400000049
The final unwrapped phase is obtained when a is 1,2, …, M, b is 1,2, …, N, and is denoted as ψ (a, b).
CN201910271422.5A 2019-04-04 2019-04-04 Multi-baseline least square phase unwrapping method based on quality map weighting Active CN110109100B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201910271422.5A CN110109100B (en) 2019-04-04 2019-04-04 Multi-baseline least square phase unwrapping method based on quality map weighting

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201910271422.5A CN110109100B (en) 2019-04-04 2019-04-04 Multi-baseline least square phase unwrapping method based on quality map weighting

Publications (2)

Publication Number Publication Date
CN110109100A CN110109100A (en) 2019-08-09
CN110109100B true CN110109100B (en) 2022-05-03

Family

ID=67485256

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201910271422.5A Active CN110109100B (en) 2019-04-04 2019-04-04 Multi-baseline least square phase unwrapping method based on quality map weighting

Country Status (1)

Country Link
CN (1) CN110109100B (en)

Families Citing this family (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110793463B (en) * 2019-09-25 2020-11-10 西安交通大学 Unwrapped phase error detection and correction method based on phase distribution
CN110907932B (en) * 2019-11-26 2022-03-18 上海卫星工程研究所 Distributed InSAR satellite height measurement precision influence factor analysis method and system
CN113238227B (en) * 2021-05-10 2022-09-30 电子科技大学 Improved least square phase unwrapping method and system combined with deep learning
CN114265062B (en) * 2021-11-11 2023-11-10 电子科技大学 InSAR phase unwrapping method based on phase gradient estimation network

Citations (3)

* 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
CN107730491A (en) * 2017-10-13 2018-02-23 北京工业大学 A kind of phase unwrapping package method based on Quality Map
CN108008383A (en) * 2017-11-09 2018-05-08 电子科技大学 A kind of four FFT phase unwrapping methods of more baseline high accuracy of iteration

Family Cites Families (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US7054504B2 (en) * 1999-02-25 2006-05-30 Ludwig Lester F Relative optical path phase reconstruction in the correction of misfocused images using fractional powers of the fourier transform
JP3498734B2 (en) * 2000-08-28 2004-02-16 セイコーエプソン株式会社 Image processing circuit, image data processing method, electro-optical device, and electronic apparatus
US7409293B2 (en) * 2004-06-03 2008-08-05 Honeywell International Inc. Methods and systems for enhancing accuracy of terrain aided navigation systems

Patent Citations (3)

* 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
CN107730491A (en) * 2017-10-13 2018-02-23 北京工业大学 A kind of phase unwrapping package method based on Quality Map
CN108008383A (en) * 2017-11-09 2018-05-08 电子科技大学 A kind of four FFT phase unwrapping methods of more baseline high accuracy of iteration

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
"A phase unwrapping method based on minimum cost flows method in irregular network";Yu Yong 等;《IEEE International Geoscience and Remote Sensing Symposium 》;20021231;第1726-1728页 *
"Least-squares two-dimensional phase unwrapping using FFT"s";M.D. Pritt 等;《IEEE Transactions on Geoscience and Remote Sensing》;19941231;第706-708页 *
蒋留兵 等."基于误差迭代补偿的干涉合成孔雷达四向加权最小二乘相位解缠法".《科学技术与工程 》.2017, *

Also Published As

Publication number Publication date
CN110109100A (en) 2019-08-09

Similar Documents

Publication Publication Date Title
CN110109100B (en) Multi-baseline least square phase unwrapping method based on quality map weighting
CN109633648B (en) Multi-baseline phase estimation device and method based on likelihood estimation
Pascazio et al. Multifrequency InSAR height reconstruction through maximum likelihood estimation of local planes parameters
CN104459633B (en) Wavelet field InSAR interferometric phase filtering method in conjunction with local frequency estimation
CN103487809B (en) A kind of based on BP algorithm and time become the airborne InSAR data disposal route of baseline
US5774089A (en) Method to resolve ambiguities in a phase measurement
CN103439708B (en) Polarized InSAR interferogram estimation method based on generalized scattering vector
CN102955157B (en) Fast correlation coefficient method for interferometric synthetic aperture radar image precise registration
CN103698764A (en) Interferometric synthetic aperture radar imaging method under sparse sampling condition
CN105954750A (en) Strip-map synthetic aperture radar non-sparse scene imaging method based on compressed sensing
CN104730519A (en) High-precision phase unwrapping method adopting error iteration compensation
Girard et al. Sparse representations and convex optimization as tools for LOFAR radio interferometric imaging
CN103630898B (en) To the method that multi-baseline interference SAR phase bias is estimated
Li et al. An interferometric phase noise reduction method based on modified denoising convolutional neural network
Baselice et al. Multibaseline SAR interferometry from complex data
Li et al. Least squares-based filter for remote sensingimage noise reduction
Hu et al. Inverse synthetic aperture radar imaging exploiting dictionary learning
Fornaro et al. Inversion of wrapped differential interferometric SAR data for fault dislocation modeling
Hu et al. Widely-distributed radar imaging based on consensus ADMM
CN115015929A (en) Efficient high-precision InSAR phase filtering network based on sparse model drive
CN112946644B (en) Based on minimizing the convolution weight l1Norm sparse aperture ISAR imaging method
Zhu et al. Beyond the 12m TanDEM-X dem
Chaubey et al. Improvement in InSAR phase unwrapping using external DEM
Liu et al. Efficient nonuniform Fourier reconstruction for spaceborne/airborne bistatic SAR
CN110426708A (en) A kind of satellite-borne SAR guarantor's phase imaging method based on orientation multichannel

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