CN114371477A - InSAR multi-baseline phase unwrapping method - Google Patents

Abstract

The invention discloses an InSAR multi-baseline phase unwrapping method, and belongs to the technical field of interferometric synthetic aperture radars. The invention provides an optimization model for establishing phase gradient ambiguity by using 2 norm constraint, and converting a solved problem into a mixed integer quadratic programming problem; the energy function expression form of the Markov random field model is improved, so that the model is converted into a quadratic programming problem to simplify the solving process. The method provided by the invention can effectively improve the execution efficiency and the solving precision of InSAR multi-baseline phase unwrapping, the solving result has the characteristics of determinacy and uniqueness, and the execution time is controllable under the given precision.

Description

InSAR multi-baseline phase unwrapping method

Technical Field

The invention relates to the technical field of Interferometric Synthetic Aperture radars (InSAR), in particular to an InSAR multi-baseline phase unwrapping method.

Background

InSAR processing enables the extraction of high-precision, large-scale Digital Surface Models (DSMs) of the ground, with irreplaceable roles in mapping applications. Phase unwrapping is a key core step of InSAR processing, and because the interferometric phase acquired through InSAR data is wrapped, namely between-pi and pi, the phase unwrapping operation is required to restore the phase to an absolute phase, and subsequent ground elevation inversion and other processing can be performed.

Phase unwrapping can be divided into single baseline phase unwrapping techniques and multi-baseline phase unwrapping techniques. The single-baseline phase unwrapping technology performs phase unwrapping through an interferogram, and mainly includes a phase unwrapping method based on path tracking, a phase unwrapping method based on norm optimization, a phase unwrapping method based on statistics, and a phase unwrapping method based on network planning. The phase unwrapping method based on path tracking adopts branch tangent lines to connect phase residual points, so that positive and negative residual points are balanced, and then the branch tangent lines are avoided to integrate the phase to obtain an unwrapped phase, and a phase island is easily generated by the method; a norm minimization objective function of the difference between the unwrapping phase gradient and the wrapping phase gradient (such as 0 norm, 1 norm and 2 norms) is established based on a norm optimization phase unwrapping method, the phase gradient is obtained by solving an optimization objective, and an unwrapping phase is obtained by integration; the phase unwrapping method based on statistics restores an absolute phase according to statistical information of an interference phase, and phase unwrapping and filtering processing are often performed at the same time, so that the problem of global diffusion of errors is easily caused; the phase unwrapping method based on network planning converts the residual minimization problem into the minimization problem of network cost flow, and utilizes the minimum cost flow network planning algorithm to perform phase unwrapping, so that the method is the most effective single-baseline phase unwrapping method at present.

The multi-baseline phase unwrapping technology carries out phase unwrapping through a plurality of interferograms, and by utilizing the combination processing of the interferograms corresponding to different long baselines and short baselines, the gradient of a phase jump part can be correctly estimated, and the phase unwrapping precision is improved. The multi-baseline phase unwrapping technology is divided into two steps of planning, wherein in the first step, the relation among baseline size, phase ambiguity and phase gradient is established, the phase gradient is estimated by utilizing Chinese remainder theorem, in the second step, the minimization criterion of unwrapping phase and wrapping unwrapping under gradient constraint is established, and an objective function is solved to obtain an unwrapping phase. However, the current multi-baseline phase unwrapping technology still has the problems of low algorithm execution efficiency, high noise interference susceptibility and the like.

Disclosure of Invention

Aiming at the limitation of the current InSAR multi-baseline phase unwrapping technology, the invention provides an InSAR multi-baseline phase unwrapping method, which uses Mixed Integer Quadratic Programming (MIQP) to improve the execution efficiency of gradient estimation and uses a Markov Random Field (MRF) model to improve the anti-noise capability of the phase unwrapping method, thereby realizing the improvement of InSAR multi-baseline unwrapping calculation efficiency and calculation accuracy.

The technical scheme adopted by the invention is as follows:

an InSAR multi-baseline phase unwrapping method comprises the following steps:

(1) inputting n interferograms

And a corresponding vertical base line B₁,B₂,…,B_nEstablishing a 2-norm optimization model objective function f of the phase gradient ambiguity by utilizing the relation among interference phases, ambiguities, vertical baselines and terrain heights among different interferograms;

(2) converting the target function f into a standard quadratic expression form, and solving the target function by using mixed integer quadratic programming to obtain phase gradient ambiguity k₁,k₂,...,k_n；

(3) In obtaining phase gradient ambiguity k₁,k₂,...,k_nUnder the condition of (1), establishing a Markov random field energy function E which can express both phase smoothness and phase closure;

(4) estimating a weight coefficient of a Markov random field energy function E by using a numerical method, wherein a smoothness weight coefficient s is obtained by maximizing a neighborhood phase, and a closure weight coefficient c is expressed by using a phase coherence rho of an interferogram;

(5) substituting the phase gradient ambiguity k of the first planning solution₁,k₂,...,k_nConverting the energy function E into a standard quadratic expression form, minimizing the energy function by using a quadratic programming method, and solving an unwrapping phase;

and completing InSAR multi-baseline phase unwrapping.

Further, the specific mode of the step (1) is as follows:

(101) inputting n interferograms, wherein the ith interferogram corresponds to a winding phase

Degree of ambiguity k_iPerpendicular base line B_iThe relationship to the terrain height h is as follows:

wherein, λ represents wavelength, R represents distance from antenna phase center to ground target point, and θ represents antenna incident angle;

with equal terrain heights h, the equation is obtained:

(102) using a permutation and combination mode, combining two by two to obtain a difference, taking a 2 norm and then summing:

further, the markov random field energy function E is:

wherein E (phi)_i,j) Indicating the phase phi at the pixel point (i, j)_i,jThe energy function of (2) includes two parts respectively representing smoothness and closure of the unwrapping phase, C_i,jRepresenting a set of phase values, s, adjacent to the pixel (i, j)_i,jWeight representing smoothness, c_i,jRepresenting the weight of the closure.

Further, the specific mode of the step (5) is as follows:

(501) substituting calculated phase gradient ambiguity k₁,k₂,...,k_nThe complex exponential part of the energy function E is removed, and the expression form after expansion and simplification is as follows:

wherein s is^rAnd s^cRepresenting the weight s in the row and column directions, respectively, c^rAnd c^cRepresenting the weights of the weights c in the row and column directions, k, respectively^rAnd k^cRespectively representing the fuzziness of the fuzziness k in the row direction and the column direction;

(502)E(φ_i,j) Is to a phase point phi in the interference pattern_i,jAssuming that the size of the interference pattern is M multiplied by N pixels, modeling is carried out on all points in the interference pattern, and the model is made

Then after summing i, j, it is converted and written in matrix form as follows:

e (Ψ) is an energy function for the whole interferogram, wherein Ψ represents a one-dimensional vector obtained by straightening a two-dimensional interferogram φ in columns;

φ_i,jthe ith row and jth column element, Ψ, of the matrix φ_tThe t-th element representing the vector Ψ;

φ_i,jto Ψ_tThe mapping of (a) is: t ═ i + N · j, i.e.:

Ψ＝[φ_0,0,φ_1,0,…,φ_1,N-1,…,φ_M-1,N-1]^T

f represents a primary coefficient vector, F_tDenotes the t-th element, phi, of the vector F_i,jCoefficient of (D) to F_tThe mapping of (a) is: t ═ i + N · j, i.e.:

h represents a quadratic coefficient matrix, H_i,jRepresents the ith row and jth column element, phi, of the matrix H_l,mφ_p,qCoefficient of (b) to H_i,jThe mapping of (a) is: i is l + m.m, j is p + m.q;

(503) for energy function

And performing minimum iterative solution to obtain a smoothed unwrapping phase psi, and finally readjusting psi into a two-dimensional matrix form phi, wherein phi is the obtained unwrapping phase.

Compared with the background technology, the invention has the following advantages:

1. aiming at the problem that the ambiguity of the phase gradient is low in efficiency when the conventional InSAR multi-baseline phase unwrapping technology uses 1 norm optimization to solve the phase gradient, the invention provides the method for solving the phase gradient by using 2 norm optimization and the time complexity is O (K)^n(n-1)/2) Reduced to O (n)³) And the calculation efficiency is improved.

2. Aiming at the problems of low efficiency and error diffusion caused by a heuristic global optimization method used for solving the Markov random field problem conventionally, the invention uses a Markov model to establish the relation between a unwrapping phase and an wrapping phase and the ambiguity of the phase, improves the expression form of the Markov random field model, and enables the Markov random field model to be solved by using a quadratic programming mode, thereby improving the efficiency on one hand and improving the noise coping capability of the algorithm on the other hand.

3. The solving algorithms used by the invention are all realized by adopting an optimization algorithm, the solving result has the characteristics of determinacy and uniqueness, and the execution time is controllable under the given precision.

Drawings

Fig. 1 is an InSAR multiple baseline phase unwrapping flow diagram.

FIG. 2 is a diagram of a Markov random field model showing the interaction between the true phase, the wrapped phase, and the unwrapped phase of an interferogram.

FIG. 3 is a diagram of a neighborhood pixel set for a target pixel.

Detailed Description

The invention is further described below with reference to the accompanying drawings.

An InSAR multi-baseline phase unwrapping method is based on Markov random fields and mixed integer quadratic programming and can recover from an interferogram wrapped phase to an absolute phase.

Suppose there are n interferograms:

the size of the interference pattern is M multiplied by N, and the corresponding vertical baseline lengths are respectively as follows: b is₁,B₂,...B_nSolving the unwrapping phase: phi is a₁,φ₂,...,φ_n。

As shown in fig. 1, the method comprises the following steps:

(1) according to the InSAR observation geometric model, the winding phase of the ith interferogram

Degree of ambiguity k_iPerpendicular base line B_iThe relationship to the terrain height h is as follows:

where λ represents the wavelength, R represents the distance from the antenna phase center to the ground target point, θ represents the antenna incident angle, and using the equal terrain heights h, the equation can be obtained:

a new 2-norm objective function is provided, a permutation and combination mode is used, pairwise combination is carried out, the difference is obtained, the 2 norms are taken, and then the sum is obtained:

because of n interference patterns, two-by-two combination can be obtained

And (4) a combination mode.

(2) The expression form of converting the objective function f into a standard quadratic form is as follows:

wherein:

X＝[k₁,k₂,...,k_n]^T,X∈Zⁿ

representing the ambiguity vector to be solved, the vector length being n,

representing a diagonal matrix of observed wrapped phase values,

defining:

P＝4π·sum_col(ΦQ)

sum_col(-) denotes the column-wise summation of the matrices, resulting in a row vector.

Is a constant independent of the unknowns and can be disregarded, so that the objective function f (x) can be simplified to:

solving a quadratic programming problem

Obtaining the ambiguity X, and further obtaining the phase gradient.

(3) Under the condition of obtaining phase gradient ambiguity, establishing an energy function model E of the Markov random field, as shown in FIG. 2, wherein the energy function is composed of two parts of representing phase smoothness and phase closure, and the phase phi at the pixel point (i, j)_i,jThe energy function of (a) is:

C_i,jrepresenting the set of phase values adjacent to pixel (i, j), s, as shown in figure 3_i,jIs the weight of smoothness, c_i,jIs the weight of the closeness.

(4) Estimating the weight coefficient of the Markov random field energy function, the smoothness weight coefficient s_i,jBy maximizing the neighborhood phase acquisition, the closure weight coefficient c_i,jPhase coherence representation using interferograms:

smoothness weight s_i,jThe estimation method comprises the following steps:

closure weight c_i,jUse ofThe coherence coefficient is calculated:

S₁and S₂Two complex SAR images corresponding to the interferograms are generated respectively, E (-) represents expectation, and superscript denotes conjugation operation.

(5) Substituting the phase gradient ambiguity k of the first planning solution₁,k₂,...,k_nThe complex exponential part of the energy function E can be removed, and the expression form after expansion and simplification is:

respectively by s^rAnd s^cWeight representing the weight s in the row and column directions, denoted by k^rAnd k^cRepresenting the ambiguity of the ambiguity k in the row and column directions.

Upper face E (phi)_i,j) Is for a phase point phi_i,jIf all points are considered, let

Then i, j after summing can be converted and written as a matrix as follows:

the above equation is the overall energy function.

Ψ represents a one-dimensional vector for straightening a two-dimensional interferogram, the vector having a length of MN, φ_i,jTo Ψ_tMapping of (2): t ═ i + N · j, i.e.:

Ψ＝[φ_0,0,φ_1,0,…,φ_1,N-1,…,φ_M-1,N-1]^T

f represents a primary coefficient vector, phi_i,jCoefficient of (D) to F_tMapping of (2): t ═ i + N · jNamely:

h represents a quadratic coefficient matrix with the size of MN multiplied by MN, phi_l,mφ_p,qCoefficient of (b) to H_i,jThe mapping of (a) is: i is l + m.m, and j is p + m.q.

For energy function

After the minimization iterative solution is performed, the smoothed unwrapped phase Ψ can be obtained.

(6) Finally, the psi is readjusted to be in a two-dimensional matrix form phi, and InSAR multi-baseline phase unwrapping is completed.

In summary, the present invention uses a mixed integer quadratic programming method to estimate the phase gradient ambiguity of the interferogram, and uses a Markov random field model to solve the unwrapped phase. The invention provides an optimization model for establishing phase gradient ambiguity by using 2 norm constraint, and converting a solved problem into a mixed integer quadratic programming problem; the invention improves the energy function expression form of the Markov random field model, so that the model is converted into a quadratic programming problem to simplify the solving process. The method provided by the invention can effectively improve the execution efficiency and the solving precision of InSAR multi-baseline phase unwrapping, the solving result has the characteristics of determinacy and uniqueness, and the execution time is controllable under the given precision.

Claims (4)

1. An InSAR multi-baseline phase unwrapping method is characterized by comprising the following steps:

(1) inputting n interferograms

And a corresponding vertical base line B₁，B₂，...，B_nEstablishing 2 norm optimal of phase gradient ambiguity by using the relation among interference phase, ambiguity, vertical base line and terrain height among different interferogramsModeling an objective function f;

(2) converting the target function f into a standard quadratic expression form, and solving the target function by using mixed integer quadratic programming to obtain phase gradient ambiguity k₁，k₂，...，K_n；

(3) In obtaining phase gradient ambiguity k₁，k₂，...，k_nUnder the condition of (1), establishing a Markov random field energy function E which can express both phase smoothness and phase closure;

(4) estimating a weight coefficient of a Markov random field energy function E by using a numerical method, wherein a smoothness weight coefficient s is obtained by maximizing a neighborhood phase, and a closure weight coefficient c is expressed by using a phase coherence rho of an interferogram;

(5) substituting the phase gradient ambiguity k of the first planning solution₁，k₂，...，k_nConverting the energy function E into a standard quadratic expression form, minimizing the energy function by using a quadratic programming method, and solving an unwrapping phase;

and completing InSAR multi-baseline phase unwrapping.

2. The InSAR multi-baseline phase unwrapping method as recited in claim 1, wherein the step (1) is specifically performed by:

(101) inputting n interferograms, wherein the ith interferogram corresponds to a winding phase

Degree of ambiguity k_iPerpendicular base line B_iThe relationship to the terrain height h is as follows:

wherein, λ represents wavelength, R represents distance from antenna phase center to ground target point, and θ represents antenna incident angle;

with equal terrain heights h, the equation is obtained:

(102) using a permutation and combination mode, combining two by two to obtain a difference, taking a 2 norm and then summing:

3. the InSAR multi-baseline phase unwrapping method as recited in claim 2, wherein the Markov random field energy function E is:

wherein E (phi)_i，j) Indicating the phase phi at the pixel point (i, j)_i，jThe energy function of (2) includes two parts respectively representing smoothness and closure of the unwrapping phase, C_i，jRepresenting a set of phase values, s, adjacent to the pixel (i, j)_i，jWeight representing smoothness, c_i，jRepresenting the weight of the closure.

4. The InSAR multi-baseline phase unwrapping method as recited in claim 3, wherein the step (5) is specifically performed by:

(501) substituting calculated phase gradient ambiguity k₁，k₂，...，K_nThe complex exponential part of the energy function E is removed, and the expression form after expansion and simplification is as follows:

wherein s is^rAnd s^cRespectively representWeight s weight in row and column directions, c^rAnd c^cRepresenting the weights of the weights c in the row and column directions, k, respectively^rAnd k^cRespectively representing the fuzziness of the fuzziness k in the row direction and the column direction;

(502)E(φ_i，j) Is to a phase point phi in the interference pattern_i，jAssuming that the size of the interference pattern is M multiplied by N pixels, modeling is carried out on all points in the interference pattern, and the model is made

Then after summing i, j, it is converted and written in matrix form as follows:

e (Ψ) is an energy function for the whole interferogram, wherein Ψ represents a one-dimensional vector obtained by straightening a two-dimensional interferogram φ in columns;

φ_i，jthe ith row and jth column element, Ψ, of the matrix φ_tThe t-th element representing the vector Ψ;

φ_i，jto Ψ_tThe mapping of (a) is: t ═ i + N · j, i.e.:

Ψ＝[φ_0，0，φ_1，0，...，φ_1，N-1，...，φ_M-1，N-1]^T

f represents a primary coefficient vector, F_tDenotes the t-th element, phi, of the vector F_i，jCoefficient of (D) to F_tThe mapping of (a) is: t ═ i + N · j, i.e.:

h represents a quadratic coefficient matrix, H_i，jRow i of the representation matrix HColumn j element, phi_l，mφ_p，qCoefficient of (b) to H_i，jThe mapping of (a) is: i is l + m.m, j is p + m.q;

(503) for energy function

And performing minimum iterative solution to obtain a smoothed unwrapping phase psi, and finally readjusting psi into a two-dimensional matrix form phi, wherein phi is the obtained unwrapping phase.

Images