CN113204869B - A Phase Unwrapping Method Based on Rank Information Filtering - Google Patents

A Phase Unwrapping Method Based on Rank Information Filtering Download PDF

Info

Publication number
CN113204869B
CN113204869B CN202110465645.2A CN202110465645A CN113204869B CN 113204869 B CN113204869 B CN 113204869B CN 202110465645 A CN202110465645 A CN 202110465645A CN 113204869 B CN113204869 B CN 113204869B
Authority
CN
China
Prior art keywords
interferogram
pixel
phase
state
state variable
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
CN202110465645.2A
Other languages
Chinese (zh)
Other versions
CN113204869A (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.)
Guilin University of Electronic Technology
Original Assignee
Guilin University of Electronic Technology
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 Guilin University of Electronic Technology filed Critical Guilin University of Electronic Technology
Priority to CN202110465645.2A priority Critical patent/CN113204869B/en
Publication of CN113204869A publication Critical patent/CN113204869A/en
Application granted granted Critical
Publication of CN113204869B publication Critical patent/CN113204869B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/10—Noise analysis or noise optimisation

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Mathematical Physics (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Optimization (AREA)
  • Mathematical Analysis (AREA)
  • Computational Mathematics (AREA)
  • General Engineering & Computer Science (AREA)
  • Data Mining & Analysis (AREA)
  • Software Systems (AREA)
  • Databases & Information Systems (AREA)
  • Algebra (AREA)
  • Computing Systems (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Geometry (AREA)
  • Image Processing (AREA)

Abstract

本发明公开了一种基于秩信息滤波的相位解缠方法，本方法把基于AMPM的局部相位梯度估计技术与秩信息滤波器结合起来，建立基于秩信息滤波的相位解缠程序，把干涉图相位解缠问题转化为秩信息滤波框架下的状态估计问题；引入H‑∞算子对状态变量信息矩阵进行优化，提高状态变量估计精度；利用基于堆排序的快速路径跟踪策略来指导相位解缠路径，保证基于秩信息滤波的相位解缠程序沿高质量像元到低质量像元路径解缠干涉图。模拟数据以及实测数据实验结果表明本文算法的有效性，能从噪声缠绕干涉图获得较为稳健的结果。

Figure 202110465645

The invention discloses a phase unwrapping method based on rank information filtering. The method combines the local phase gradient estimation technology based on AMPM with the rank information filter, establishes a phase unwrapping program based on rank information filtering, and converts the interferogram phase The unwrapping problem is transformed into a state estimation problem under the rank information filtering framework; the H-∞ operator is introduced to optimize the state variable information matrix to improve the estimation accuracy of the state variable; the fast path tracking strategy based on heap sorting is used to guide the phase unwrapping path , to ensure that the rank-information filtering-based phase unwrapping procedure unwraps the interferogram along the path from high-quality pixels to low-quality pixels. The experimental results of simulated data and measured data show that the algorithm in this paper is effective, and it can obtain more robust results from noise-entangled interferograms.

Figure 202110465645

Description

Phase unwrapping method based on rank information filtering
Technical Field
The invention relates to the field of phase unwrapping, in particular to a phase unwrapping method based on rank information filtering.
Background
Phase unwrapping is one of the important steps in data processing such as interferometric synthetic aperture radar measurement (InSAR), optical interferometry, synthetic aperture sonar (InSAS), and magnetic resonance imaging. Because of the periodicity of the trigonometric function, the interference phase of the characterization target parameter obtained from interferometry is limited to the phase main value (-pi, pi) interval, commonly known as the winding phase, and the true interference phase of the reaction target parameter is recovered from the winding phase, namely the so-called phase unwrapping.
Traditional phase unwrapping methods, such as branch-cut method, least square method, quality guiding method, etc., have good unwrapping effect under the condition of less noise and less residual points, and have large unwrapping result error once the noise is large and the residual points are densely distributed. Subsequently, a class of unwrapping algorithms based on nonlinear filtering and state estimation is proposed successively, including extended kalman filter phase unwrapping algorithm (EKFPU), unscented kalman filter phase unwrapping algorithm (kfpu), volumetric kalman filter phase unwrapping algorithm (CKFPU), unscented information filter phase unwrapping algorithm (uipu), etc. These algorithms compensate to some extent for the deficiencies of conventional phase unwrapping algorithms, but efficient and accurate unwrapping of noise interferograms remains a very challenging problem.
Disclosure of Invention
The present invention addresses the shortcomings of the prior art described above by providing a method for phase unwrapping based on Rank Information Filtering (RIF) that can obtain more robust results from noise-wrapped interferograms.
The technical scheme for realizing the aim of the invention is as follows:
a phase unwrapping method based on rank information filtering comprises the following steps:
1) According to the existing nonlinear filtering phase unwrapping model, combining a local phase gradient estimation technology based on a correction matrix beam model (AMPM) with a rank information filter, establishing a phase unwrapping model based on rank information filtering, and converting an interferogram phase unwrapping problem into a state estimation problem under a rank information filtering frame;
2) According to the rank information filtering phase unwrapping model obtained in the step 1), a one-dimensional rank information filtering phase unwrapping algorithm is established;
3) And (3) according to the one-dimensional rank information filtering phase unwrapping algorithm obtained in the step (2), replacing the one-dimensional coordinates with the two-dimensional coordinates, and establishing the two-dimensional rank information filtering phase unwrapping algorithm.
Further, in step 1), using the relationship between the unwrapping phases of adjacent pixels of the interferogram, the interferogram phase unwrapping system equation can be expressed as follows:
Figure BDA0003043259350000021
Figure BDA0003043259350000022
wherein x is k Representing the unwrapping phase of the k pixels of the interferogram as a state variable to be evaluated; mu (mu) k-1 Representing the true phase gradient of the k-1 pixel of the interferogram, the estimated value can be obtained by using the local gradient estimation technology based on AMPM
Figure BDA0003043259350000023
w k-1 Representing phase gradient estimation errors; zeta type toy k Observation noise vector, ζ, for interference image k-element state variable 1,k And xi 2,k Measurement noise added to the imaginary and real parts of the complex interference signal, h x k ]Noiseless observation vectors which are the state variables of k pixels of the interference image; z k Is the observation vector of the state variable of the k pixel of the interference image.
Further, in step 2), the interference pattern k-1 pixel state estimation value and the estimation error variance are respectively set as
Figure BDA0003043259350000024
And->
Figure BDA0003043259350000025
The phase unwrapping based on the rank information filtering comprises the steps of:
2-1) generating rank information sampling points according to a rank sampling principle;
2-2) carrying out state prediction on the interference image k pixels to obtain one-step prediction values and one-step prediction error variances of the state variables of the interference image k pixels;
2-3) converting the state space of the interference image k pixels into the information space to obtain an information matrix Y of one-step prediction of the state variables of the interference image k pixels k Sum information vector
Figure BDA0003043259350000026
2-4) carrying out state update on the interference image k pixels to obtain the unwrapped phase of the interference image k pixels and the variance of the phase estimation error.
Further, in step 2-1), the method includes the steps of generating rank information sampling points according to the rank sampling principle:
Figure BDA0003043259350000027
Figure BDA0003043259350000028
in the method, in the process of the invention,
Figure BDA0003043259350000029
for the initial estimation of the state variable of the k pixels of the interference image, χ i,k Represents a rank information sample point generated according to the rank sampling principle, n represents the dimension of the state vector, n=1, u 1 And u 2 Represents a standard state offset, where u 1 ＝0.2,u 2 =0.5; i represents the number of sampling points; />
Figure BDA00030432593500000210
Representing an interference pattern k-1 pixel state estimation value; />
Figure BDA00030432593500000211
Representing an estimation error variance; />
Figure BDA00030432593500000212
Representing the value of the state variable of the k-1 picture element.
Further, in step 2-2), the one-step predicted value of the interferogram k pel state variable:
Figure BDA0003043259350000031
one-step prediction error variance of interferogram k pel state variables:
Figure BDA0003043259350000032
in the method, in the process of the invention,
Figure BDA0003043259350000033
as the weight coefficient, Q k Process error variance caused by the local gradient estimation technology based on AMPM; />
Figure BDA0003043259350000034
Representing one-step predicted values of state variables of k pixels of the interference image; />
Figure BDA0003043259350000035
Representing one-step prediction error variance of the state variable of the k pixels of the interference image; t represents a transpose operation.
Further, in step 2-3), the information matrix Y of one-step prediction of the interferogram k-pel state variables k Sum information vector
Figure BDA0003043259350000036
Figure BDA0003043259350000037
Figure BDA0003043259350000038
Information matrix Y using H-infinity operator k When the optimization is performed, the information matrix in the formula (5) becomes:
Figure BDA0003043259350000039
wherein, gamma s As the attenuation factor, 0.8.ltoreq.gamma.is usually taken s Less than or equal to 2, I is the same as
Figure BDA00030432593500000310
Identity matrix of the same dimension.
Further, the method comprises the steps of,in step 2-4), the cross covariance
Figure BDA00030432593500000311
Measurement prediction value +.>
Figure BDA00030432593500000312
The method comprises the following steps:
z i,k ＝h[χ i,k ] (8)
Figure BDA00030432593500000313
Figure BDA00030432593500000314
information status distribution i k And corresponding information matrix distribution I k The method comprises the following steps:
Figure BDA00030432593500000315
Figure BDA00030432593500000316
Figure BDA00030432593500000317
wherein eta is k Observing vector residual error for interference image k pixel state variable, R k Observing noise variance for the interferogram k pel state variable,
updated information matrix Y k And information vector y k The method comprises the following steps of:
Figure BDA0003043259350000041
Figure BDA0003043259350000042
the updated state estimation values and the state estimation error variances are respectively as follows:
Figure BDA0003043259350000043
Figure BDA0003043259350000044
wherein,,
Figure BDA0003043259350000045
representing the state estimation of the k pixels of the interference image, namely the unwrapping phase of the k pixels of the interference image; />
Figure BDA0003043259350000046
Representing the phase estimation error variance of the k pixels of the interference image; the one-dimensional RIF unwrapping algorithm can unwrap the interferogram wrapping phases in a row-by-row or column-by-column manner until all wrapping phases in the interferogram are unwrapped.
Further, in step 3), two-dimensional coordinates (m, n) are used to replace one-dimensional k, and if the interference image element (m, n) is the pixel to be unwound, the initial state predicted value of the pixel is obtained
Figure BDA0003043259350000047
And its prediction error variance->
Figure BDA0003043259350000048
The calculation can be as follows:
Figure BDA0003043259350000049
Figure BDA00030432593500000410
Figure BDA00030432593500000411
wherein the pixels of the interference pattern (a, s) are unwrapped pixels in 8 pixels adjacent to the pixel to be unwrapped,
Figure BDA00030432593500000412
and->
Figure BDA00030432593500000413
Representing state estimates of the pixels of the interferograms (a, s) and their estimated error variances; />
Figure BDA00030432593500000414
Phase gradient estimation values representing interference pixels (m, n) and (a, s) can be obtained by using an AMPM-based local gradient estimation technique; d, d (a,s) Weights representing the state estimates of the pixels of the interferograms (a, s);
rank information sampling point χ of interferogram (m, n) pixel i,(m,n) The calculation can be as follows:
Figure BDA00030432593500000415
one-step predicted value of state variable
Figure BDA00030432593500000416
Figure BDA0003043259350000051
State variable one-step prediction error variance
Figure BDA0003043259350000052
Figure BDA0003043259350000053
Wherein Q is [(m,n)|(a,s)] For the variance of process errors caused by the local gradient estimation technology based on AMPM, a fast path tracking strategy based on heap ordering is utilized to guide a phase unwrapping path, and a RIF phase unwrapping program is guided to complete recursion estimation of interferogram wrapping pixels along a path from high-quality pixels to low-quality pixels.
The method applies a high-efficiency and steady rank information filter to interferogram phase unwrapping, and combines an AMPM-based local phase gradient estimation technology with a fast path tracking strategy based on heap ordering, so as to provide a phase unwrapping algorithm based on rank information filtering. The method comprises the steps of establishing a phase unwrapping program based on rank information filtering, and acquiring phase gradient information of the phase unwrapping program based on rank information filtering by using a local phase gradient estimation technology based on AMPM; and a fast path tracking strategy based on heap ordering is utilized to guide a phase unwrapping path, so that a phase unwrapping program based on rank information filtering is ensured to unwrap an interferogram along a path from a high-quality pixel to a low-quality pixel. The simulation data and the actual measurement data experimental results show that the effectiveness of the algorithm can obtain a more robust result from the noise winding interference diagram.
Drawings
FIGS. 1 a-1 f are simulated interferometry, wherein FIGS. 1 a-1 c are true unwrapping phases of three interferometry, FIG. 1d is a noise wrapping phase diagram of the true interferometry phase of FIG. 1a, FIG. 1e is a noise wrapping phase diagram of the true interferometry phase of FIG. 1b, and FIG. 1f is a noise wrapping phase diagram of the true interferometry phase of FIG. 1 c;
FIGS. 2 a-2 c are the results of unwrapping FIG. 1d with the method of the present invention, where FIG. 2a shows unwrapped phases, FIG. 2b shows phase unwrapped errors, and FIG. 2c shows phase unwrapped error histograms;
FIGS. 3 a-3 c are the results of unwrapping FIG. 1e with the method of the present invention, where FIG. 3a shows unwrapped phases, FIG. 3b shows phase unwrapped errors, and FIG. 3c shows phase unwrapped error histograms;
FIGS. 4 a-4 c are the results of unwrapping FIG. 1f with the method of the present invention, where FIG. 4a shows unwrapped phases, FIG. 4b shows phase unwrapped errors, and FIG. 4c shows phase unwrapped error histograms;
fig. 5 a-5 c show the results of unwrapping measured data using the method of the present invention, wherein fig. 5a shows a partial Etna volcanic interference pattern, fig. 5b shows unwrapping phase, and fig. 5c shows unwrapping phase re-wrapping results.
Detailed Description
The following description of the technical solutions according to the embodiments of the present invention will be provided fully with reference to the accompanying drawings in which it is apparent that the described embodiments are only some embodiments, but not all embodiments of the present invention. All other embodiments, which can be made by those skilled in the art based on the embodiments of the invention without making any inventive effort, are intended to be within the scope of the invention.
Examples:
a phase unwrapping method based on rank information filtering comprises the following steps:
1) According to the existing nonlinear filtering phase unwrapping model, combining an AMPM-based local phase gradient estimation technology with a rank information filter, establishing a phase unwrapping model based on rank information filtering, and converting an interferogram phase unwrapping problem into a state estimation problem under a rank information filtering framework;
using the relationship between the unwrapped phases of adjacent pixels of the interferogram, the interferogram phase unwrapping system equation can be expressed as follows:
Figure BDA0003043259350000061
Figure BDA0003043259350000062
wherein x is k Representing the unwrapping phase of the k pixels of the interferogram as a state variable to be evaluated; mu (mu) k-1 Representing the true phase gradient of the k-1 pixel of the interferogram, the estimated value can be obtained by using the local gradient estimation technology based on AMPM
Figure BDA0003043259350000063
w k-1 Representing phase gradient estimation errors; zeta type toy k Observation noise vector, ζ, for interference image k-element state variable 1,k And xi 2,k Measurement noise added to the imaginary and real parts of the complex interference signal, h x k ]Noiseless observation vectors which are the state variables of k pixels of the interference image; z k An observation vector which is an interference image k pixel state variable;
2) According to the rank information filtering phase unwrapping model obtained in the step 1), a one-dimensional rank information filtering phase unwrapping algorithm is established;
let the state estimation value of the interference pattern k-1 pixel and the estimation error variance respectively be
Figure BDA0003043259350000064
And->
Figure BDA0003043259350000065
The phase unwrapping based on the rank information filtering comprises the steps of:
2-1) rank information sampling points generated according to the rank sampling principle:
Figure BDA0003043259350000066
Figure BDA0003043259350000067
in the method, in the process of the invention,
Figure BDA0003043259350000068
for the initial estimation of the state variable of the k pixels of the interference image, χ i,k Represents a rank information sample point generated according to the rank sampling principle, n represents the dimension of the state vector, n=1, u 1 And u 2 Represents a standard state offset, where u 1 ＝0.2,u 2 =0.5; i represents the number of sampling points; />
Figure BDA0003043259350000069
And->
Figure BDA00030432593500000610
Respectively representing the state estimation value and the estimation error variance of the interference image k-1 pixel;
Figure BDA00030432593500000611
representing a k-1 pel state variable value;
2-2) carrying out state prediction on the interference image k pixels to obtain one-step prediction values and one-step prediction error variances of the state variables of the interference image k pixels;
one-step predicted value of state variable of k pixels of an interferogram:
Figure BDA0003043259350000071
one-step prediction error variance of interferogram k pel state variables:
Figure BDA0003043259350000072
in the method, in the process of the invention,
Figure BDA0003043259350000073
as the weight coefficient, Q k Process error variance caused by the local gradient estimation technology based on AMPM; />
Figure BDA0003043259350000074
Representing one-step predicted values of state variables of k pixels of the interference image; />
Figure BDA0003043259350000075
Representing one-step prediction error variance of the state variable of the k pixels of the interference image; t represents a transpose operation;
2-3) converting the state space of the interference image k pixels into the information space to obtain an information matrix Y of one-step prediction of the state variables of the interference image k pixels k Sum information vector
Figure BDA0003043259350000076
Information matrix Y for one-step prediction of interference image k pixel state variable k And information vector->
Figure BDA0003043259350000077
Figure BDA0003043259350000078
Figure BDA0003043259350000079
Information matrix Y using H-infinity operator k When the optimization is performed, the information matrix in the formula (5) becomes:
Figure BDA00030432593500000710
wherein, gamma s As the attenuation factor, 0.8.ltoreq.gamma.is usually taken s Less than or equal to 2, I is the same as
Figure BDA00030432593500000711
A co-dimensional identity matrix;
2-4) carrying out state update on the interference image k pixels to obtain unwrapped phases of the interference image k pixels and a variance of phase estimation errors; cross covariance
Figure BDA00030432593500000712
Measurement prediction value +.>
Figure BDA00030432593500000713
The method comprises the following steps:
z i,k ＝h[χ i,k ] (8)
Figure BDA00030432593500000714
Figure BDA0003043259350000081
information status distribution i k And corresponding information matrix distribution I k The method comprises the following steps:
Figure BDA0003043259350000082
Figure BDA0003043259350000083
Figure BDA0003043259350000084
wherein eta is k Observing vector residual error for interference image k pixel state variable, R k Observing noise variance for the interferogram k pel state variable,
updated information matrix Y k And information vector y k The method comprises the following steps of:
Figure BDA0003043259350000085
Figure BDA0003043259350000086
the updated state estimation values and the state estimation error variances are respectively as follows:
Figure BDA0003043259350000087
Figure BDA0003043259350000088
wherein,,
Figure BDA0003043259350000089
representing the state estimation of the k pixels of the interference image, namely the unwrapping phase of the k pixels of the interference image; />
Figure BDA00030432593500000810
Representing the phase estimation error variance of the k pixels of the interference image; the one-dimensional RIF unwrapping algorithm can unwrap the interferogram wrapping phases in a row-by-row or column-by-column manner until all wrapping phases in the interferogram are unwrapped;
3) And (3) according to the one-dimensional rank information filtering phase unwrapping algorithm obtained in the step (2), replacing the one-dimensional coordinates with the two-dimensional coordinates, and establishing the two-dimensional rank information filtering phase unwrapping algorithm.
Using two-dimensional coordinates (m, n) to replace one-dimensional k, and setting an interference image element (m, n) as an element to be unwound, and predicting the initial state of the element
Figure BDA00030432593500000811
And its prediction error variance->
Figure BDA00030432593500000812
The calculation can be as follows:
Figure BDA00030432593500000813
Figure BDA00030432593500000814
Figure BDA00030432593500000815
wherein the pixels of the interference pattern (a, s) are unwrapped pixels in 8 pixels adjacent to the pixel to be unwrapped,
Figure BDA00030432593500000816
and->
Figure BDA00030432593500000817
Representing state estimates of the pixels of the interferograms (a, s) and their estimated error variances; />
Figure BDA0003043259350000091
Representing phase gradient estimates between interference pixels (m, n) and (a, s) can be obtained using AMPM-based local gradient estimation techniques [13 ]]；d (a,s) Weights representing the state estimates of the pixels of the interferograms (a, s);
rank information sampling point χ of interferogram (m, n) pixel i,(m,n) The calculation can be as follows:
Figure BDA0003043259350000092
one-step predicted value of state variable
Figure BDA0003043259350000093
Figure BDA0003043259350000094
State variable one-step prediction error variance
Figure BDA0003043259350000095
Figure BDA0003043259350000096
Wherein Q is [(m,n)|(a,s)] For the variance of process errors caused by the local gradient estimation technology based on AMPM, a fast path tracking strategy based on heap ordering is utilized to guide a phase unwrapping path, and a RIF phase unwrapping program is guided to complete recursion estimation of interferogram wrapping pixels along a path from high-quality pixels to low-quality pixels.
In order to verify the performance of each method, different algorithms comprise an iterative least squares method (ILS), a quality guide method (QGPU) and the method, the simulation and the actual measurement interferograms are unwrapped under the same MATLAB software environment (Intel i5-8265U@1.60G CPU+8GB RAM), and the unwrapped results of each algorithm are compared and analyzed.
Fig. 1 shows three different simulated interferograms of 256×256 pixels, wherein fig. 1 a-1 c show the real interference phases, fig. 1 d-1 f show the noisy winding phase diagrams, respectively, with signal to noise ratios of 7.44dB, 2.18dB, 0.73dB, and the three interferograms are unwrapped using a Rank Information Filter Phase Unwrapping (RIFPU) algorithm.
FIGS. 2 a-2 c are the results of unwrapping FIG. 1d with the method of the present invention, where FIG. 2a shows unwrapped phases, FIG. 2b shows phase unwrapped errors, and FIG. 2c shows phase unwrapped error histograms;
FIGS. 3 a-3 c are the results of unwrapping FIG. 1e with the method of the present invention, where FIG. 3a shows unwrapped phases, FIG. 3b shows phase unwrapped errors, and FIG. 3c shows phase unwrapped error histograms;
FIGS. 4 a-4 c are the results of unwrapping FIG. 1f with the method of the present invention, where FIG. 4a shows unwrapped phases, FIG. 4b shows phase unwrapped errors, and FIG. 4c shows phase unwrapped error histograms;
fig. 5 a-5 c show the results of unwrapping measured data using the method of the present invention, wherein fig. 5a shows a partial Etna volcanic interference pattern, fig. 5b shows unwrapping phase, and fig. 5c shows unwrapping phase re-wrapping results.
The result of unwrapping the noise-wrapped phase diagrams shown in fig. 1 d-1 f by the method is shown in fig. 2 a-4 c, and it can be seen that the unwrapping phase obtained by the method has better consistency with the real interference diagram, and the phase unwrapping error is smaller. The first list lists the root mean square error of Iterative Least Squares (ILS), quality guided methods (QGPU) and the method unwrapping different signal-to-noise interference patterns, it can be seen that the root mean square error of the method is much smaller than the QGPU and ILS methods.
Table one phase unwrapping error for each algorithm
Figure BDA0003043259350000101
In the experimental data, fig. 5a is a partial Etna volcanic interference diagram, the unwrapping phase of the method is shown in fig. 5b, and the unwrapping phase is shown in fig. 5 c. The unwrapping phase obtained by the method is continuous and consistent, and the rewinding phase diagram stripes are consistent with the original interference diagram stripes, which shows that the method obtains more effective unwrapping results.
The preferred embodiments of the invention disclosed above are merely to aid in the description of the invention and are not intended to limit the invention to the specific embodiments described. Obviously, many modifications and variations are possible in light of the above teaching. The embodiments were chosen and described in order to best explain the principles of the invention and the practical application, to thereby enable others skilled in the art to best understand and utilize the invention.

Claims (5)

1.一种基于秩信息滤波的相位解缠方法，其特征在于，包括如下步骤：1. A phase unwrapping method based on rank information filtering, is characterized in that, comprises the steps: 1)根据已有的非线性滤波相位解缠模型，把基于AMPM的局部相位梯度估计技术与秩信息滤波器结合起来，建立基于秩信息滤波的相位解缠模型，把干涉图相位解缠问题转化为秩信息滤波框架下的状态估计问题；1) According to the existing nonlinear filtering phase unwrapping model, the local phase gradient estimation technology based on AMPM is combined with the rank information filter to establish a phase unwrapping model based on the rank information filtering, and the interferogram phase unwrapping problem is transformed into is the state estimation problem under the framework of rank information filtering; 2)根据步骤1)中得到的秩信息滤波相位解缠模型，建立一维秩信息滤波相位解缠算法；2) According to the rank information filtering phase unwrapping model obtained in step 1), a one-dimensional rank information filtering phase unwrapping algorithm is established; 3)根据步骤2)中得到的一维秩信息滤波相位解缠算法，用二维坐标代替一维坐标，建立二维秩信息滤波相位解缠算法；3) According to the one-dimensional rank information filtering phase unwrapping algorithm obtained in step 2), replace the one-dimensional coordinates with two-dimensional coordinates, and establish a two-dimensional rank information filtering phase unwrapping algorithm; 步骤1)中，利用干涉图相邻像元解缠相位之间的关系，干涉图相位展开系统方程可表示为如下：In step 1), the relationship between the phases of the adjacent pixels of the interferogram is unwrapped, and the system equation of the interferogram phase unwrapping can be expressed as follows:
Figure FDA0004241825440000011
Figure FDA0004241825440000011
Figure FDA0004241825440000012
Figure FDA0004241825440000012
 式中，xk表示干涉图k像元解缠相位，作为待评估的状态变量；μk-1表示干涉图k-1像元真实相位梯度，可利用基于AMPM的局部梯度估计技术来获取其估计值
Figure FDA0004241825440000013
wk-1表示相位梯度估计误差；ξk为干涉图k像元状态变量的观测噪声矢量，ξ1,k和ξ2,k分别为附加在复干涉信号虚部和实部的量测噪声，h[xk]为干涉图k像元状态变量的无噪声观测矢量；zk为干涉图k像元状态变量的观测矢量；In the formula, x k represents the unwrapped phase of pixel k in the interferogram, which is used as the state variable to be evaluated; μ k-1 represents the real phase gradient of pixel k-1 in the interferogram, and the AMPM-based local gradient estimation technique can be used to obtain its estimated value
Figure FDA0004241825440000013
w k-1 represents the phase gradient estimation error; ξ k is the observation noise vector of the interferogram k pixel state variable, ξ 1,k and ξ 2,k are the measurement noise added to the imaginary part and real part of the complex interference signal, respectively , h[x k ] is the noise-free observation vector of the interferogram k pixel state variable; z k is the observation vector of the interferogram k pixel state variable; 步骤2)中，设干涉图k-1像元状态估计值及其估计误差方差分别为
Figure FDA0004241825440000014
和/>
Figure FDA0004241825440000015
则基于秩信息滤波的相位解缠包括如下步骤：In step 2), it is assumed that the estimated value of the interferogram k-1 pixel state and its estimated error variance are respectively
Figure FDA0004241825440000014
and />
Figure FDA0004241825440000015
Then the phase unwrapping based on rank information filtering includes the following steps: 2-1)根据秩采样原理生成秩信息采样点；2-1) Generate rank information sampling points according to the rank sampling principle; 2-2)对干涉图k像元进行状态预测，得到干涉图k像元状态变量一步预测值以及一步预测误差方差；2-2) Carry out state prediction to interferogram k pixel, obtain interferogram k pixel state variable one-step prediction value and one-step prediction error variance; 2-3)对干涉图k像元进行状态空间向信息空间转换，得到干涉图k像元状态变量一步预测的信息矩阵Yk和信息向量
Figure FDA0004241825440000016
2-3) Convert the interferogram k pixel from the state space to the information space, and obtain the information matrix Y k and information vector of the one-step prediction of the interferogram k pixel state variable
Figure FDA0004241825440000016
 2-4)对干涉图k像元进行状态更新，得到干涉图k像元解缠相位以及相位估计误差方差；2-4) Update the state of the interferogram k pixel to obtain the unwrapped phase of the interferogram k pixel and the phase estimation error variance; 步骤3)中，用二维坐标(m,n)代替一维k，设干涉图像元(m,n)为待解缠像元，则该像元初始状态预测值
Figure FDA0004241825440000017
及其预测误差方差/>
Figure FDA0004241825440000018
可按如下计算：In step 3), use two-dimensional coordinates (m, n) to replace one-dimensional k, and set the interference image element (m, n) as the pixel to be unwrapped, then the predicted value of the initial state of the pixel
Figure FDA0004241825440000017
and its forecast error variance/>
Figure FDA0004241825440000018
Can be calculated as follows:
Figure FDA0004241825440000021
Figure FDA0004241825440000021
Figure FDA0004241825440000022
Figure FDA0004241825440000022
Figure FDA0004241825440000023
Figure FDA0004241825440000023
 式中，干涉图(a,s)像元为待解缠像元相邻8个像元中的已解缠像元，
Figure FDA0004241825440000024
和/>
Figure FDA0004241825440000025
表示干涉图(a,s)像元的状态估计值及其估计误差方差；/>
Figure FDA0004241825440000026
表示干涉像元(m,n)与(a,s)之间的相位梯度估计值，可利用基于AMPM的局部梯度估计技术获取；d(a,s)表示干涉图(a,s)像元状态估计的权值；In the formula, the interferogram (a, s) pixel is the unwrapped pixel among the 8 adjacent pixels of the pixel to be unwrapped,
Figure FDA0004241825440000024
and />
Figure FDA0004241825440000025
Represents the state estimation value of the interferogram (a, s) pixel and its estimation error variance; />
Figure FDA0004241825440000026
Indicates the phase gradient estimation value between the interferometric pixel (m,n) and (a,s), which can be obtained using AMPM-based local gradient estimation technology; d (a,s) represents the interferogram (a,s) pixel The weight of the state estimate; 干涉图(m,n)像元的秩信息采样点χi,(m,n)，可按如下计算：The rank information sampling point χ i,(m,n) of the interferogram (m,n) pixel can be calculated as follows:
Figure FDA0004241825440000027
Figure FDA0004241825440000027
 状态变量一步预测值
Figure FDA0004241825440000028
State variable one-step predictor
Figure FDA0004241825440000028
Figure FDA0004241825440000029
Figure FDA0004241825440000029
 状态变量一步预测误差方差
Figure FDA00042418254400000210
State variable one-step forecast error variance
Figure FDA00042418254400000210
Figure FDA00042418254400000211
Figure FDA00042418254400000211
 其中，w为权重系数；n表示状态向量的维数；
Figure FDA00042418254400000212
表示干涉图像元(m,n)初始状态预测值；u1和u2表示标准状态偏量；/>
Figure FDA00042418254400000213
表示干涉图像元(m,n)预测误差方差；Q[(m,n)|(a,s)]为基于AMPM的局部梯度估计技术导致的过程误差方差，利用基于堆排序的快速路径跟踪策略来指导相位解缠路径，引导RIF相位解缠程序沿高质量像元到低质量像元的路径完成对干涉图缠绕像元递推估计。Among them, w is the weight coefficient; n represents the dimension of the state vector;
Figure FDA00042418254400000212
Indicates the initial state prediction value of the interference image element (m,n); u 1 and u 2 represent the standard state deviation; />
Figure FDA00042418254400000213
Represents the prediction error variance of the interference image element (m,n); Q [(m,n)|(a,s)] is the process error variance caused by the AMPM-based local gradient estimation technique, using a fast path tracking strategy based on heap sorting To guide the phase unwrapping path, guide the RIF phase unwrapping program to complete the recursive estimation of the interferogram wrapped pixels along the path from high-quality pixels to low-quality pixels. 2.根据权利要求1所述的基于秩信息滤波的相位解缠方法，其特征在于，步骤2-1)中，包括根据秩采样原理生成的秩信息采样点：2. the phase unwrapping method based on rank information filtering according to claim 1, is characterized in that, in step 2-1), comprises the rank information sampling point that generates according to rank sampling principle:
Figure FDA0004241825440000031
Figure FDA0004241825440000031
Figure FDA0004241825440000032
Figure FDA0004241825440000032
 式中，
Figure FDA0004241825440000033
为干涉图k像元状态变量初始估计值，χi,k表示根据秩采样原理生成的秩信息采样点，n表示状态向量的维数，n＝1，u1和u2表示标准状态偏量，其中u1＝0.2,u2＝0.5；i表示采样点个数；/>
Figure FDA0004241825440000034
和/>
Figure FDA0004241825440000035
分别表示干涉图k-1像元状态估计值及其估计误差方差；/>
Figure FDA0004241825440000036
表示k-1像元状态变量值。In the formula,
Figure FDA0004241825440000033
is the initial estimated value of the state variable of interferogram k pixel, χ i, k represent the rank information sampling points generated according to the rank sampling principle, n represents the dimension of the state vector, n=1, u 1 and u 2 represent the standard state deviation , where u 1 =0.2, u 2 =0.5; i represents the number of sampling points;/>
Figure FDA0004241825440000034
and />
Figure FDA0004241825440000035
Respectively represent the estimated value of the k-1 pixel state of the interferogram and its estimated error variance; />
Figure FDA0004241825440000036
Indicates the k-1 pixel state variable value. 3.根据权利要求1所述的基于秩信息滤波的相位解缠方法，其特征在于，步骤2-2)中，干涉图k像元状态变量一步预测值：3. the phase unwrapping method based on rank information filtering according to claim 1, is characterized in that, in step 2-2), interferogram k pixel state variable one-step predictive value:
Figure FDA0004241825440000037
Figure FDA0004241825440000037
 干涉图k像元状态变量一步预测误差方差：Interferogram k-pixel state variable one-step prediction error variance:
Figure FDA0004241825440000038
Figure FDA0004241825440000038
 式中，
Figure FDA0004241825440000039
为权重系数，Qk为基于AMPM的局部梯度估计技术导致的过程误差方差；/>
Figure FDA00042418254400000310
表示干涉图k像元状态变量一步预测值；/>
Figure FDA00042418254400000311
表示干涉图k像元状态变量一步预测误差方差；T表示转置运算；χi,k表示根据秩采样原理生成的秩信息采样点；u1和u2表示标准状态偏量；n表示状态向量的维数。In the formula,
Figure FDA0004241825440000039
is the weight coefficient, Q k is the process error variance caused by the AMPM-based local gradient estimation technique; />
Figure FDA00042418254400000310
Indicates the one-step predictive value of the interferogram k pixel state variable; />
Figure FDA00042418254400000311
Indicates the one-step prediction error variance of the interferogram k pixel state variable; T indicates the transposition operation; χ i,k indicates the rank information sampling points generated according to the rank sampling principle; u 1 and u 2 indicate the standard state deviation; n indicates the state vector of dimensions. 4.根据权利要求1所述的基于秩信息滤波的相位解缠方法，其特征在于，步骤2-3)中，干涉图k像元状态变量一步预测的信息矩阵Yk和信息向量
Figure FDA00042418254400000312
4. the phase unwrapping method based on rank information filtering according to claim 1, is characterized in that, in step 2-3), the information matrix Y k and information vector of interferogram k pixel state variable one-step prediction
Figure FDA00042418254400000312
Figure FDA00042418254400000313
Figure FDA00042418254400000313
Figure FDA00042418254400000314
Figure FDA00042418254400000314
 利用H-∞算子对信息矩阵Yk进行优化，则式(5)中的信息矩阵变为：Using the H-∞ operator to optimize the information matrix Y k , the information matrix in formula (5) becomes:
Figure FDA00042418254400000315
Figure FDA00042418254400000315
 式中，γs为衰减因子，取0.8≤γs≤2，I是与
Figure FDA00042418254400000316
同维的单位矩阵；/>
Figure FDA00042418254400000317
表示干涉图k像元状态变量一步预测值；/>
Figure FDA0004241825440000041
表示干涉图k像元状态变量一步预测误差方差。In the formula, γ s is the attenuation factor, taking 0.8≤γ s ≤2, I is the same as
Figure FDA00042418254400000316
identity matrix of the same dimension; />
Figure FDA00042418254400000317
Indicates the one-step predictive value of the interferogram k pixel state variable; />
Figure FDA0004241825440000041
Indicates the one-step prediction error variance of the interferogram k-pixel state variable. 5.根据权利要求1所述的基于秩信息滤波的相位解缠方法，其特征在于，步骤2-4)中，互协方差
Figure FDA0004241825440000042
及量测预测值/>
Figure FDA0004241825440000043
为：5. the phase unwrapping method based on rank information filtering according to claim 1, is characterized in that, in step 2-4), mutual covariance
Figure FDA0004241825440000042
and measured predicted value/>
Figure FDA0004241825440000043
for: zi,k＝h[χi,k] (8)z i,k ＝h[χ i,k ] (8)
Figure FDA0004241825440000044
Figure FDA0004241825440000044
Figure FDA0004241825440000045
Figure FDA0004241825440000045
 信息状态分布ik和相应的信息矩阵分布Ik为：The information state distribution i k and the corresponding information matrix distribution I k are:
Figure FDA0004241825440000046
Figure FDA0004241825440000046
Figure FDA0004241825440000047
Figure FDA0004241825440000047
Figure FDA0004241825440000048
Figure FDA0004241825440000048
 式中，ηk为干涉图k像元状态变量观测矢量残差，Rk为干涉图k像元状态变量观测噪声方差，In the formula, η k is the observation vector residual of the interferogram k pixel state variable, R k is the observation noise variance of the interferogram k pixel state variable, 更新的信息矩阵Yk和信息向量yk分别为：The updated information matrix Y k and information vector y k are respectively:
Figure FDA0004241825440000049
Figure FDA0004241825440000049
Figure FDA00042418254400000410
Figure FDA00042418254400000410
 更新后的状态估计值及其状态估计误差方差分别为：The updated state estimates and their state estimation error variances are:
Figure FDA00042418254400000411
Figure FDA00042418254400000411
Figure FDA00042418254400000412
Figure FDA00042418254400000412
 其中，
Figure FDA00042418254400000413
表示干涉图k像元状态估计，亦即干涉图k像元解缠相位；/>
Figure FDA00042418254400000414
表示干涉图k像元相位估计误差方差；w为权重系数；zk为干涉图k像元状态变量的观测矢量；/>
Figure FDA00042418254400000415
表示干涉图k像元状态变量一步预测值；/>
Figure FDA00042418254400000416
表示用H-∞算子对信息矩阵Yk进行优化后的信息矩阵；n表示状态向量的维数；上述一维RIF解缠算法按逐行或逐列的方式对干涉图缠绕相位进行解缠，直至干涉图中所有缠绕相位解缠完毕。in,
Figure FDA00042418254400000413
Indicates the interferogram k pixel state estimation, that is, the interferogram k pixel unwrapping phase; />
Figure FDA00042418254400000414
Represents the phase estimation error variance of pixel k in the interferogram; w is the weight coefficient; z k is the observation vector of the state variable of pixel k in the interferogram; />
Figure FDA00042418254400000415
Indicates the one-step predictive value of the interferogram k pixel state variable; />
Figure FDA00042418254400000416
Indicates the information matrix optimized by the H-∞ operator on the information matrix Y k ; n indicates the dimension of the state vector; the above one-dimensional RIF unwrapping algorithm unwraps the interferogram winding phase in a row-by-row or column-by-column manner , until all the wrapped phases in the interferogram are unwrapped.
CN202110465645.2A 2021-04-28 2021-04-28 A Phase Unwrapping Method Based on Rank Information Filtering Active CN113204869B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110465645.2A CN113204869B (en) 2021-04-28 2021-04-28 A Phase Unwrapping Method Based on Rank Information Filtering

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110465645.2A CN113204869B (en) 2021-04-28 2021-04-28 A Phase Unwrapping Method Based on Rank Information Filtering

Publications (2)

Publication Number Publication Date
CN113204869A CN113204869A (en) 2021-08-03
CN113204869B true CN113204869B (en) 2023-06-23

Family

ID=77029221

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110465645.2A Active CN113204869B (en) 2021-04-28 2021-04-28 A Phase Unwrapping Method Based on Rank Information Filtering

Country Status (1)

Country Link
CN (1) CN113204869B (en)

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110440935A (en) * 2019-08-12 2019-11-12 桂林电子科技大学 A kind of phase developing method based on Extended information filter
CN111724307A (en) * 2020-06-19 2020-09-29 山东财经大学 An image super-resolution reconstruction method based on maximum posterior probability and non-local low-rank prior, terminal and readable storage medium

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP4758572B2 (en) * 2001-07-27 2011-08-31 富士フイルム株式会社 Phase unwrapping method for fringe image analysis

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110440935A (en) * 2019-08-12 2019-11-12 桂林电子科技大学 A kind of phase developing method based on Extended information filter
CN111724307A (en) * 2020-06-19 2020-09-29 山东财经大学 An image super-resolution reconstruction method based on maximum posterior probability and non-local low-rank prior, terminal and readable storage medium

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
Efficient phase unwrapping algorithm based on cubature information particle filter applied to unwrap noisy continuous phase maps;Xianming Xie等;《Optics Express》;20190325;第27卷(第7期);9906-9924 *
Phase unwrapping algorithm based on a rank information filter;Xianming Xie等;《Applied Optics》;20210729;第60卷(第22期);6648-6658 *
基于高级InSAR时序分析方法的高速公路沉降分析;张庆云等;《科学技术与工程》;20180718;第18卷(第20期);20-26 *

Also Published As

Publication number Publication date
CN113204869A (en) 2021-08-03

Similar Documents

Publication Publication Date Title
JP2018195069A (en) Image processing apparatus and image processing method
CN112099007B (en) Azimuth multi-channel SAR fuzzy suppression method suitable for non-ideal antenna directional diagram
JP7027365B2 (en) Signal processing equipment, signal processing methods and programs
CN113589286B (en) Unscented Kalman filtering phase unwrapping method based on D-LinkNet
CN110018476A (en) Time difference baseline set timing interference SAR processing method
CN113506212B (en) Improved hyperspectral image super-resolution reconstruction method based on POCS
Zhou et al. Nonlocal means filtering based speckle removal utilizing the maximum a posteriori estimation and the total variation image prior
CN112597433A (en) Plug and play neural network-based Fourier phase recovery method and system
CN119471684A (en) Signal processing optimization method and system for spaceborne multi-band synthetic aperture radar
CN115453533B (en) A Phase Optimization Method for Distributed Targets in Time-Series InSAR
CN113204869B (en) A Phase Unwrapping Method Based on Rank Information Filtering
CN113920255B (en) High-efficient mapping system based on point cloud data
CN118505557B (en) Construction of point cloud denoising model, and point cloud denoising method and device
CN111079893A (en) Method and apparatus for obtaining generator network for interference fringe pattern filtering
JP7512150B2 (en) Information processing device, information processing method, and program
Genser et al. Spectral constrained frequency selective extrapolation for rapid image error concealment
CN110856014A (en) Moving image generation method, moving image generation device, electronic device, and storage medium
CN117315437A (en) Phase shift value estimation method and system of phase shift interferogram based on convolutional neural network
CN114814728B (en) A sound source localization method, system, electronic device and medium
CN117168444A (en) Kalman filtering method and system based on iteration error state
Li et al. 3D shape measurement based on Res-Attention-Unet for deep learning
CN113163201A (en) Video multi-frame reconstruction method and device based on single-pixel camera
CN110108307A (en) A kind of terrible imaging method of adjustable orthogonalization depth
CN119251091B (en) Image denoising method and system based on adversarial frequency mixing technology
CN115482892A (en) Composite material component multi-source noise filtering method based on digital speckles

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