The content of this article comes from the 4th issue of "Journal of Surveying and Mapping" in 2024 (drawing review number: GS Jing (2024) No. 0714)
InSAR phase dewrapping maximum/minimum cut value improvement algorithm
Gao Yandong
,1,2, Jia Yikun,1,2, Li Shijin1,2, Chen Yu1,2, Li Huaizhan1,2, Zheng Nanshan1,2, Zhang Shubi1,2
1. Key Laboratory of Land Environment and Disaster Monitoring, Ministry of Natural Resources, China University of Mining and Technology, Xuzhou 221116, Jiangsu, China
2. School of Environment and Geomatics, China University of Mining and Technology, Xuzhou 221116, Jiangsu, China
Abstract:InSAR has been widely used in high-precision DEM inversion, and phase unwrapping technology is one of the key steps affecting the accuracy of DEM acquisition, but the large gradient region has always been the core problem affecting the accuracy of unwrapping results. In order to solve this problem, this paper proposes a maximum flow/minimum secant phase dewrapping algorithm based on the improved weight of potential function. Firstly, in order to solve the problem of unreasonable weight setting of PUMA model, the prior information of phase gradient change is obtained by using the existing external DEM, and the maximum absolute phase gradient value of the window is substituted into the corresponding potential function formula to obtain the weight. Then, by adjusting the threshold value of the potential function weight, the problem of unwrapping error caused by the failure of the PUMA potential function due to the unreasonable setting of the potential function weight was solved, and the phase unwrapping accuracy of the large gradient change region was improved. Finally, the proposed algorithm is verified by simulation data and real TanDEM-X InSAR data, and compared with the existing methods. The results show that the proposed algorithm can improve the unwrapping accuracy by at least 44.93% in the simulation data, and in the large gradient change region in the real data, the proposed algorithm can obtain a larger range of effective unwrapping results than the existing algorithms. Key words: InSAR ; potential function ; Maximum flow/minimum cut; Phase unwrapping
Fund:This work was financially supported by the National Natural Science Foundation of China (42001409)
Author: Gao Yan-dong (1988-), Male, Ph.D., Associate Professor, Research Direction为单/多基线相位解缠, InSAR High Precision DEM Counter-Performance. E-mail:[email protected]
This article cites format
Gao Yandong, Jia Yikun, Li Shijin, Chen Yu, Li Huaizhan, Zheng Nanshan, Zhang Shubi. Improved algorithm for the maximum stream/minimum cut value of InSAR phase dewrapping. Journal of Surveying and Mapping[J], 2024, 53(4): 644-652 doi:10.11947/j.AGCS.2024.20220633
GAO Yandong. The improved max-flow/min-cut weight algorithm for InSAR phase unwrapping. Acta Geodaetica et Cartographica Sinica[J], 2024, 53(4): 644-652 doi:10.11947/j.AGCS.2024.20220633
Read more
InSAR has been widely used in the fields of DEM acquisition, mine safety monitoring, and landslide disaster early warning [1-3]. Phase unwrapping is one of the key steps in InSAR data processing, and the accuracy of the results will directly affect the accuracy of the final product [4]. Phase unwrapping aims to obtain the true phase from the phase entangled between (-π,π] [5-6]. However, the phase unwrapping accuracy in the large gradient region is not ideal, which is one of the main factors affecting the accuracy of the final product obtained by InSAR [7].
Phase unwrapping has always been a technical difficulty and research hotspot in the field of InSAR data processing, which can be divided into two main categories [8]. The first type is the path following method represented by the branch cutting method [9], the mass graph method [10] and the minimum discontinuity method [11], which can quickly and accurately obtain the unwrapping results by calculating the residual points, setting the branch tangent lines, and optimizing the unwrapping integration path. However, this kind of method will produce error transmission and even unwrapping voids in low-quality areas, which will seriously affect the accuracy of unwrapping results. The other type is the gradient estimation optimization method represented by the Lp norm method [12], the minimum-cost flow (MCF) method [13], and the statistical cost network flow (SNAPHU) [14], which aims to minimize the difference between the phase gradient estimate and the real phase gradient to obtain the final unwrapping result, and this type of method has good stability. However, such methods are prone to global error transmission, which affects the final unwrapping result [15]. In addition, the phase continuity assumption is a common criterion for unwrapping methods, which assumes that in the absence of noise, the true phase gradient is equal to the winding phase gradient [16], i.e., the phase gradient is limited to (-π,π], however, in the real world, especially in the region of large gradient variation, this assumption is not applicable, so most unwrapping methods cannot obtain ideal results in the region of large gradient variation [17]. In order to solve the above problems, a multi-baseline phase unwrapping method was proposed. Among them, the representative multi-baseline phase unwrapping methods mainly include the Chinese remainder theorem [18] and the two-stage programming method [19], which can get rid of the limitation of phase continuity assumption by combining the interference phases of different baseline lengths, and can also obtain more ideal unwrapping results in the region of large gradient changes. However, the multi-baseline phase unwrapping method has the problem of poor noise robustness, and not only that, not all regions have the interferometric phase data required by the multi-baseline method at this stage, which is also one of the main factors limiting the application of such methods. Therefore, it is particularly important to obtain higher precision unwrapping results in large gradient regions through single-amplitude interferometric phase.
In recent years, the study of phase unwrapping in the region with large gradient variation of single-amplitude interference phase has received extensive attention. Ref. [20] proposes a phase unwinding strategy based on the assistance of the mining settlement prediction model, which uses the mining settlement prediction model to predict the deformation of the mining area in combination with the subsidence parameters of the mining area, and then iteratively processes it to solve the unwinding problem of the large gradient change area. Ref. [21] proposes a phase dewrapping method with the assistance of prior information, which can solve the phase unwrapping problem in large gradient regions by combining prior information. In Ref. [22], a slice phase unwrapping algorithm based on external DEM assistance was proposed, which simulates the solvable phase of a single baseline through an external DEM, so as to slice the unwrapped phase, so as to reduce the phase gradient change and solve the phase unwrapping problem in the large gradient change area. Based on the Markov random field (MRF), the phase unwrapping max-flow/min-cut algorithm (PUMA) was proposed for the first time in Ref. [23], and it has been proved that the convex potential function has a good effect on dealing with continuous phases, and the non-convex potential function is more suitable for dealing with large gradient phases. Therefore, PUMA has obvious advantages for unwinding in large gradient regions. In Ref. [24], a linear process MRF (LP-MRF) method was proposed, which combined with the non-convex potential function and set the weights by using an 8-neighborhood linear window to detect phase continuity, which improved the phase unwrapping ability of PUMA in the large gradient region. An improved PUMA (IM-PUMA) method was proposed in Ref. [25], which uses the quasi-convex potential function to adaptively identify the large gradient phase, and uses the external DEM to obtain the maximum gradient prior information to assist in setting the weights, which improves the accuracy of the unwrapping results in the large gradient change region. Although IM-PUMA can obtain ideal results, however, the setting of the weight of the model is polarized in the large gradient change region, and the weight setting in the flat terrain area does not reach the ideal large value, and in the case of extreme large gradient change, the weight setting will lead to a certain unwrapping error due to the small weight setting, so the weight setting of the IM-PUMA model still needs to be further improved.
Based on the IM-PUMA method, this paper proposes a PUMA based on potential function to improve the weight. Firstly, the external auxiliary DEM is used to obtain the prior phase gradient information, and the maximum phase gradient of the window is substituted into the corresponding potential function formula, and the weight is obtained by normalization. The method of setting the threshold value of the weight is used to solve the problem of unwrapping error caused by the potential function of the weight setting too small in the area of extreme large gradient change, and the overall phase unwrapping accuracy is improved, so as to obtain high-precision DEM products. After the simulation data and real TanDEM-X InSAR data experiments, the results show that the proposed algorithm can obtain higher accuracy results than the existing algorithms in the large gradient change region.
1 Phase unwrapping theory
1.1 Basic principles of phase unwrapping
The phase information obtained by the actual measurement of InSAR technology is the principal value of the true phase, that is, the winding phase, and its value is between (-π, π). An integer multiple of 2π between the winding phase and the true phase, and the relationship between the true phase and the winding phase is
(1)
where s is the pixel; ψ(s) is the true aspect; φ(s) is the winding phase; k(s) is the gradient ambiguity coefficient. For single-baseline phase unwrapping, the phase gradient is constrained by the phase continuity assumption (-π,π), which can be expressed as
(2)
where Δψ(s) is the true phase gradient between pixel s and pixel s-1. This assumption is a necessary premise for the single-baseline phase unwrapping algorithm at this stage, but the true phase gradient change in the large gradient change region is often greater than the winding phase gradient change, which is also one of the main factors for the phase unwrapping error.
1.2 IM-PUMA方法
PUMA is a phase unwrapping method based on Bayesian theory, the accuracy of PUMA unwrapping results depends on two prior information, potential function and weight setting, traditional PUMA sets the weight required by the potential function through the coherence coherence, but the coherence coefficient is not enough to represent the phase continuity between adjacent pixels, so it will cause some interference to the potential function, resulting in the increase of the unwrapping error. Therefore, in view of the shortcomings of the traditional PUMA model, the IM-PUMA method was proposed in Ref. [25], which adopts a more reasonable quasi-convex potential function to adaptively identify the large gradient phase, and uses the gradient information of the external DEM to assist in setting the weights, which greatly increases the reliability of the potential function weights. Based on PUMA, the IM-PUMA method abstracts the interferogram into an MRF, where the nodes of the MRF correspond to the interferogram pixels, and the edges between the nodes correspond to the winding phase gradient Δφ between each pair of pixels in the interferogram. Based on MRF, the energy function model E(k|ψ) of the IM-PUMA method is constructed
(3)
where p is the prime point of the interferometric image; P is the set of interferometric image pixels; NP is a neighborhood system of pixel p; Δψpq is the true phase gradient between pixel p and adjacent pixel q; Δkpq is the difference between the gradient ambiguity coefficient between pixel p and adjacent pixel q. Kp and Kq are the gradient blurring coefficients of pixel P and adjacent pixel Q, respectively. WPQ stands for weight, which is an indicator of the phase continuity between pixel p and adjacent pixel q. According to the calculation of the maximum gradient prior information of the externally assisted DEM, the pixel-symmetric cluster Cpq, V(·) composed of pixel p and adjacent pixel q is a potential function
(4)
where a is the adjustment parameter. According to the IM-PUMA energy function model, the value of the energy function of the t-iteration is defined as
(5)
Define δ∈{0,1} and minimize the energy function of the t-th iteration, then the value of the energy function of the t+1 iteration is expressed as
(6)
During the ceremony,
and
The phase ambiguity coefficient estimates of pixel P at the T and T+1 iterations are described, respectively. while
，
and
The phase ambiguity coefficient estimates of the adjacent pixel q in the t and t+1 iterations are deputies, respectively. while
, at the initial moment
，
and
The phase gradient estimates for the t and t+1 iterations are described, respectively. From Eq. (5) and Eq. (6), it can be seen that the energy minimization problem can be mapped to a binary optimization problem, and the energy function value of the group Cpq in Eq. (6).
can be expressed as
、
、
、
These four energy expressions are mapped to a "directed basic graph" with assigned values according to the method of literature [26], and a "directed graph" with assigned values is obtained by combining the directed basic graphs of all groups, and the energy function value of the t+1st iteration can be obtained by performing the maximum flow/minimum cut operation on the complete directed graph [27], until the energy function value is no longer reduced, and the final unwrapping phase can be obtained. Although IM-PUMA can obtain better results than traditional PUMA in the large gradient region, the potential function weight of IM-PUMA is only obtained by gradient value, and there is a polarization problem in the weight setting of the large gradient change region, especially in the more extreme large gradient change region, the potential function will not function normally in the process of energy minimization, resulting in some defects of IM-PUMA in the large gradient change region.
2 PUMA algorithm with improved weights based on potential function
In order to solve the above problems, this paper proposes a PUMA algorithm based on the improved weight of potential function based on the IM-PUMA algorithm. Using the definition of the maximum gradient similar to the IM-PUMA algorithm, C is the maximum weight, G is the maximum absolute phase gradient of the window, and G(i,j) is calculated in a sliding window with a size of 2n+1 centered on pixels (i,j), taking the 3×3 window as an example, as shown in Figure 1.
Figure 1
Figure 1 Example of a 3×3 window
Fig. 1 3×3 window example
The absolute phase gradients of the window in the horizontal, vertical, 45° and 135° directions within the window are denoted as Gh, Gv, Gd1 and Gd2, respectively. The phase gradient in the 4 directions can be expressed as
(7)
According to Eq. (7), the maximum absolute phase gradient of the window at pixels (i,j) is
(8)
It is worth noting that the difference between the maximum absolute phase gradient of the window proposed in this paper and the maximum gradient in IM-PUMA is that the maximum gradient of the window in IM-PUMA can be negative, while the maximum absolute phase gradient of the window proposed in this paper is based on the absolute value of the gradient values in 4 directions, that is, there is a constraint
(9)
Due to the existence of Eq. (9) constraint, the polarization phenomenon in the weight setting of the large gradient change region after the weight normalization of IM-PUMA can be completely avoided, and the weight is always small in the region where the large gradient rises or falls, which conforms to the distribution principle that the weight of the potential function should be small at the large gradient change.
In IM-PUMA, only the maximum gradient of all pixels (with negative values) is normalized to obtain the weights, which is not ideally combined with the quasi-convex potential function model in the IM-PUMA energy function model. Similarly, if only the maximum absolute phase gradient G(i,j) of the window is normalized to obtain the weight for energy minimization, the effect of polarization in the large gradient change region can only be eliminated, and the phase continuity information can be roughly expressed. Therefore, in order to make better use of the phase gradient to reasonably represent the phase continuity information, this paper proposes to further optimize the weights of the maximum absolute phase gradient G(i,j) of the window combined with the quasi-convex potential function formula, and the maximum weights at pixels(i,j) can be expressed as
(10)
The maximum weight C of all pixels is normalized so that the value of C is within the interval [0,1].
(11)
where Cmin is the minimum value in C; Cmax is the maximum value in C. In this case, the normalized maximum weight C′ is not the final potential function weight. C′(i,j) close to 1 indicates that there is a large gradient at the pixel (i,j), and a small weight should be assigned to the corresponding potential function. C′(i,j) close to 0 indicates that the terrain at pixel (i,j) is gentle, and a large weight should be given to the corresponding potential function. Therefore, the potential function weight w(i,j) of the energy function model at pixels (i,j) is given by Eq. (12).
(12)
The potential function weight expression of the interferogram energy function model is
(13)
It is worth noting that in Eq. (10), the square term is given to the maximum absolute phase gradient G(i,j) of the window, which makes the originally small C′(i,j) smaller, and combined with Eq. (12), it can be seen that the value of w(i,j) will increase, that is, the operation of assigning the square term to the maximum absolute gradient G(i,j) of the window increases the weight of the flat terrain area, so that the weight of the flat terrain area may be close to 1, which is in line with the principle that the phase of the flat terrain area should have high confidence. In addition, the weights calculated by using the formula of the quasi-convex function itself can be more reasonably matched with the quasi-convex function model, which further improves the credibility of the phase continuity information in the large gradient region.
In the region of extreme large gradient variation, G(i,j) will appear to be very large compared to other regions, that is, the weight w(i,j) in this case will be very close to 0, which will cause the potential function model to not function normally. This is because the weight considers this to be an almost vertical fault, and the effect of the potential function on the fault should be erased as much as possible. However, there is also a certain gradient in the more extreme large gradient change region of the practical problem, so it is obviously unreasonable for the weight w(i,j) to be close to 0 in this case. In order to solve the above problems, this paper proposes a way to solve the problem of too small weights by setting a threshold for the weights.
Assuming that the weight threshold is M, the obtained potential function weight w is limited by the threshold, and the weight expression W of the potential function is limited to
(14)
If w(i,j)=M, then the maximum absolute phase gradient of the window at the pixel (i,j) is said to be the critical gradient GM, which can be determined according to the external auxiliary DEM data, and the weight threshold can be calculated according to the critical gradient GM and Eq. (12). The specific method of obtaining the critical gradient GM is as follows: calculate the phase gradient of the external DEM azimuth and distance respectively, count the distribution law of all phase gradients, and determine the GM according to the phase gradient distribution, generally GM corresponds to the first 10%~30% interval of all phase gradient values, and the phase gradient exceeding the GM value can be regarded as the more extreme large gradient change in the interferogram, and the weight of the large gradient region where the maximum absolute gradient value exceeds the GM needs to be compensated. In this way, the problem that the potential function cannot play a role due to the weight being too small can be solved. The value threshold M is in the range of [0.1, 0.3].
3 Test results and analysis
In this paper, the effectiveness of the proposed algorithm is verified and analyzed by using simulation data and real TanDEM-X InSAR data, and compared with MCF, SNAPHU and IM-PUMA methods. In this paper, the root mean square error is used as a quantitative evaluation index of phase unwrapping accuracy
(15)
where N is the total number of pixels;
is the unwrapping phase;
is the reference unwrapping phase.
3.1 Simulation data test
Figure 2 shows the simulation data used in this paper, and the data size is 458×157. Figure 2(a) shows the simulation-aided unwinding phase; Fig. 2(b) shows the corresponding winding interference phase; Fig. 2(c) shows the weights calculated by the IM-PUMA algorithm based on the simulated unwinding phase. Fig. 2(d) shows the weight calculated by the algorithm based on the simulated unwinding phase, and the weight threshold of the simulation data in this paper is 0.3. As can be seen from Fig. 2(a) and Fig. 2(b), there is an obvious large gradient change region in the simulation data, which creates a good condition for verifying the performance of different phase unwrapping algorithms in the large gradient change region. In order to identify large gradient regions more finely, the window size for calculating the weights of the simulation data is set to 3×3.
Figure 2
Fig.2 Experimental dataset of simulation data
Fig. 2 Simulation data experimental dataset
As shown in Figure 3, Figure 3 (a)-(d) are the unwrapping results of MCF, SNAPHU, IM-PUMA and the proposed algorithm, respectively. Fig. 3(e)-(h) are the errors of the unwrapping results of different methods and Fig. 2(a), respectively, and it can be seen that MCF and SNAPHU have obvious unwrapping errors, and IM-PUMA also has obvious unwrapping patches in the more extreme large gradient region due to the unreasonable weight setting. As can be seen from Figure 3(h), although the proposed algorithm has also produced some unwrapping errors, it has eliminated most of the unwrapping patches. In this paper, RMSE is used to quantitatively evaluate the unwrapping error of different methods, as shown in Table 1, due to the influence of the more extreme large gradient changes, the RMSE results of MCF, SNAPHU and IM-PUMA unwrapping results are not ideal, and the proposed algorithm can better solve the unwrapping problem in the large gradient change region, so its RMSE is the highest accuracy among several unwrapping methods, which also verifies that the proposed algorithm can obtain better unwrapping results than the existing methods in the more extreme large gradient change region.
Figure 3
Fig.3. The results and errors of different phase unwrapping methods are tested on simulation data
Fig. 3 Unwrapping results and error plots for different phase unwrapping methods with simulated data
表1 模拟数据的评估结果Tab. 1 Evaluation results of simulated data
| Phase unwrapping method | RMSE |
|---|---|
| DBM | 1.633 4 |
| SNAPHU | 1.385 3 |
| IM-PUMA | 1.999 5 |
| This article is an algorithm | 0.762 9 |
A new window opens| Download the CSV
3.2 Real data test
In order to further verify the effectiveness of the proposed algorithm, the real TanDEM-X data in Weinan area were selected for processing, as shown in Figure 4. Fig. 4(a) is the Goolgle image of the real data, Fig. 4(b) is the unwrapping phase simulated by the external auxiliary DEM, Fig. 4(c) is the real interference phase of TanDEM-X, Fig. 4(d) is the weight of the IM-PUMA, and Fig. 4(e) is the weight of the algorithm in this paper, and the weight threshold of the real data in this paper is 0.3. As can be seen from Figure 4(a), there is a large area of large gradient change in the study area, which will seriously affect the phase unwrapping accuracy. According to the situation of the real data study area, the window size for calculating the real data weight is set to 3×3.
Figure 4
Figure 4 Experimental dataset of real data
Fig. 4 Real data experimental dataset
Fig. 5 shows the unwrapping results and error diagrams of different phase unwrapping methods, and Fig. 5 (a) and (e) are the unwrapping results and corresponding unwrapping errors of the MCF method, respectively, and there is an obvious unwrapping jump error in the unwrapping results of the MCF method, which also shows that the MCF method is prone to produce large unwrapping errors in the large gradient change region. Fig. 5 (b) and (f) are the unwrapping results and errors of the SNAPHU method, respectively, although the SNAPHU method can obtain better unwrapping results than the MCF method, but there is still a large area of unwrapping and jumping errors. Fig. 5 (c) and (g) are the unwrapping results and errors of IM-PUMA, respectively, although there is still a certain area of unwrapping and jumping error, compared with the first two methods, the unwrapping results of IM-PUMA are significantly improved. However, due to the problem of weight setting, there is still a partial unwrapping error. Fig. 5 (d) and (h) are the unwrapping results and errors of the proposed algorithm, respectively, compared with the first three methods, the proposed algorithm obtains the most ideal unwrapping results, and there is the least area of unwrapping jump error. In order to further demonstrate the performance of the proposed unwrapping method, this paper uses RMSE to compare and analyze different methods (Table 2), and the RMSE of both MCF and SNAPHU methods is larger due to the influence of unwrapping error in large gradient change region. The RMSE accuracy of IM-PUMA is improved compared with the first two methods, but it is still not very ideal. The RMSE of the proposed algorithm is the best among the four methods, which further proves that the proposed algorithm can achieve better results than the existing methods in the large gradient change region.
Figure 5
Fig.5 The unwrapping results and errors of different phase unwrapping methods are tested on real data
Fig. 5 Unwrapping results and error plots for different phase unwrapping methods with real data
表2 真实TanDEM-X数据的评估结果Tab. 2 Evaluation results of real TanDEM-X data
| Phase unwrapping method | RMSE |
|---|---|
| DBM | 5.387 0 |
| SNAPHU | 4.952 0 |
| IM-PUMA | 2.866 3 |
| This article is an algorithm | 2.442 8 |
A new window opens| Download the CSV
4 Summary and outlook
Phase unwrapping in large gradient regions is one of the key problems affecting high-precision DEM inversion. Based on the IM-PUMA method, this paper mainly studies the problem of unreasonable weight setting of the algorithm model, and proposes a maximum flow/minimum secant phase unwrapping algorithm based on the improved weight of the potential function. Combined with the prior phase gradient information of the external DEM, a weight setting method based on the maximum absolute phase gradient estimation of the window is proposed, and the weights are obtained by using the quasi-convex potential function formula, which makes the setting of the weights more reasonable to match with the potential function. In addition, in order to solve the problem of unwrapping error in the more extreme large gradient change region, this paper uses the method of setting a threshold for the weight to solve the problem of unwrapping error caused by the weight setting too small potential function in the extreme large gradient change region, so as to improve the phase unwrapping accuracy of the large gradient change region. The experimental results show that in the simulation data, compared with the existing phase unwrapping methods, the proposed algorithm can obtain higher accuracy results in the large gradient change region, and the accuracy of the unwrapping results can be increased by 44.93%, which better solves the problem of unwrapping error caused by unreasonable IM-PUMA weight setting. In the real TanDEM-X InSAR data processing, due to the problem of large gradient change, the existing unwrapping algorithms cannot obtain the ideal unwrapping results, especially the MCF produces a wide range of unwrapping errors, the proposed algorithm adopts a more scientific weight setting method to effectively solve the problem of unwrapping error in the large gradient change region, and the experimental results show that the proposed algorithm can obtain better unwrapping results than the existing methods.
First trial: Zhang Yanling review: Song Qifan
Final Judge: Jin Jun
information