Phase unwrapping method based on parallel local minimum reliability dual expanding for large-scale data

Abstract. Considering data reliability, the optimized connections of dual residues for phase discontinuity reconstruction provide a more reliable scheme and produce more robust unwrapped results. However, their practical implementation usually involves a time-consuming iterative global operation, which is not suitable for application to the phase unwrapping (PU) of big blocks of interferometric synthetic aperture radar (InSAR) phase data. A parallel PU method based on local minimum reliability dual expanding is proposed. With given quality weight maps, the dual reliability is defined based on residues, and the minimum reliability residue pairs are introduced to represent possible discontinuous boundaries. We provide a dynamic expanding method for the pairs with local minimum reliability searching and dual merging. The final minimum balanced trees obtained are used for path integrating of PU with the aid of reliability maps. The calculation of reliability maps, the residue pair searching, and the dynamic expanding are designed to be carried out in parallel. We employ the interface propagation scheme based on the eikonal equation and flood-filling for the parallel implement. Two big blocks of airborne InSAR data were processed by the proposed method, the experimental results and analysis verified its robustness and efficiency for the large-scale PU problem.


Introduction
Interferometric synthetic aperture radar (InSAR) is an important Earth observation technology that has been developing for decades, and two-dimensional (2-D) phase unwrapping (PU) 1 is a key part of the InSAR interferometry process. PU retrieves absolute phase values from ambiguous modulo 2π phase in order to reconstruct three-dimensional (3-D) information from observation targets. The main ideas of PU are based on an assumption of Nyquist's theorem that the differences between neighbor sampling phases are less than half a cycle. Based on this assumption, the difference values can be obtained precisely using the wrapping function. The difficulty with PU is that the assumption is untenable under certain circumstances, such as where there is highly undulating terrain or systematic noises are present. In such cases, the problem is ill-posed, and the neighbor phase difference can be more than half a cycle, which means that its unique value cannot be specified directly. PU works mainly by focusing on direct or indirect detection and identification of reliable paths for unwrapping integration by bringing in more auxiliary constraints, which are expected to exclude the abnormal neighbor phase differences.
The data quality indices usually play a key part in existing PU methods. They can be used as weights in the least square methods, modulating the optimization criteria. [2][3][4] Moreover, quality map-guided unwrapping methods based on data quality indices frequently appear in the research literature. Recently, improved schemes have been proposed, such as interwoven linked list, 5 priority queuing, 6,7 and quality mapping with edges and shadow areas. 8 Another implemented form of them, tiling strategies 9,10 based on residue clustering characteristics, have recently been provided with a convex hull algorithm for more effective schemes. The quality indices work more indirectly in reconstructing discontinuity compared to others based on residues, such as branchcuts methods. The classic branch-cut method 11 connects residues to form discontinuity boundaries, in line with rules such as minimal distances. The simplicity of its principle brings about high efficiency on performance, but its drawbacks of islands or holes are problems for noisy areas.
More reliable methods are those relying on optimization with graph theory, for instance, the minimum cost flow (MCF) method, 12 the graph-cuts method, 13 and the statistical-cost network flow algorithm (SNAPHU). 14,15 These methods are usually time-consuming for large blocks of data, because of their iterative optimization processing and a large number of residues. Data tiling is usually also considered to help improve the efficiency of the method. Regular data tiling strategies are applied in SNAPHU, 14 the branch-cut, 16 and quality-guided 17 methods; irregular data tiling is also applied in the MCF method, 18,19 which uses subgraphs involving partial residues to reduce the workload and in the extended minimum cost flow method, 20 which splits the MCF network into simpler subnetworks. Another strategy for dealing with large-scale data is parallel computing; for instance, parallel algorithms are designed in PU based on field programmable gate array discrete cosine transform 21 and graphics processing unit (GPU) regiongrowing 22 technologies.
In this paper, we propose a regular-grid PU method based on a parallel local minimum reliability dual expanding scheme, which is used to reconstruct the discontinuity boundaries. First, based on the definition of dual reliability, the minimum reliability residue pairs are introduced to represent possible discontinuous boundaries. We provide a dynamic expanding method for the pairs with local minimum reliability searching and dual extension. The final minimum balanced trees (MBTs) obtained are used for path integrating of PU with the aid of reliability maps. The calculation of reliability maps, the residue pair searching, and the dynamic expanding are all designed to be implemented in parallel. In Sec. 2, we provide the concepts of dual reliability, local minimum reliability dual expanding, and MBTs. In Sec. 3, we give the implement schemes of the structure models, parallel interface propagation, dynamic expanding and connection lines tracing, and the flood-filling. In Sec. 4, two PU experiments with wrapped phase data from airborne InSAR are carried out with the proposed method, and the process and final results are analyzed. The study conclusions are presented in Sec. 5.

Principle and Method
In the 2-D definition domain, the unwrapped phase Φ at ðx; yÞ can be given as E Q -T A R G E T ; t e m p : i n t r a l i n k -; s e c 2 ; 1 1 6 ; 2 5 3 Φðx; yÞ where ∇Ψ refers to the gradient of the wrapped phase Ψ on path l and W to the wrapping function, and l is an integrating path from the reference unwrapped phase Φ 0 to the target. The core idea of most PU methods focuses on directly or indirectly detecting and shielding phase discontinuity boundaries so that integration of path l can be reliably designed. We try to find a dual reliability measure for the dual residues on discontinuity boundaries. For example, for the two discontinuity boundaries c 1 and c 2 from four coupled residues in Fig. 1(a), the constructed dual reliability map makes the discontinuity boundaries enveloped in low-reliability areas as shown in (b).

The Dual Reliability Definition Based on Phase Data Quality
For a given phase quality weight Q in the defined domain, if there is a reference residue S and a target at ðx; yÞ, the reliability P of the target is defined as an integration of the quality field following a path c, namely (1) where Q is a positive increasing function of the goodness of the data, Qðu; vÞ refers to the weight at ðu; vÞ on c, s refers to an arc length measure of c, and c refers to a special curved path from S to the target, which minimizes the integration. There are many potential candidates for this path that are based on popular quality factors, such as the coherence coefficient and the derivative variance.
If there are many reference residues, there will be multiple values of the integration reliability at one target position, and we select the smallest integration reliability as the final definition of the minimum risk. At ðx; yÞ reliability, P þ based on positive reference residues and P − based on negative ones are defined as E Q -T A R G E T ; t e m p : i n t r a l i n k -; s e c 2 . 1 ; 1 1 6 ; 4 4 7 where R þ and R − refer to positive and negative reference residue sets and the selected references for the two reliability values are labeled as S þ and S − . The sum of P þ and P − gives the dual reliability: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 2 ; 1 1 6 ; 3 8 0 Pðx; yÞ ¼ P þ ðx; yÞ þ P − ðx; yÞ: (2) The phase reliability P þ , P − , and P provide comprehensive reliability descriptions of phase data. If a position is with a low-reliability value, it should be near to a positive or negative residue or in a low-quality area, which suggests that it is tending to located around a discontinuous boundary.

Minimum Reliability Residue Pairs and Local Minimum Reliability Dual Expanding
According to definitions of phase reliability, we can get reliability P þ and P − and corresponding reference residues S þ and S − at any position, including residues. For a positive residue S a and a negative residue S b , if the selected reference negative residue at the position of S a is S b , and the selected reference positive residue at the position of S b is S a , we call the residues S a and S b a minimum reliability residue pair. The residues in the pair both have minimum reliability values derived from each other, and the direct connecting line, cðS a ; S b Þ, from the positive residue to the negative, is with the minimum reliability values; that is The typical value of dual reliability P on the pair connecting line c is constant according to the definitions. However, residues with the same sign may interfere with each other in reliability calculating, the dual reliability on disturbed connecting lines may be not constant. In any case, the connecting lines should be "valley" lines of the dual reliability map. Examples are shown in Fig. 2. (4) where D is the union of the sets of the optimal connecting lines of all the residue pairs, and the formula meets a weighted minimum L 0 norm of the PU. It is NP-hard for the global optimization, and we try to solve the problem with a scheme of local minimum reliability expanding.
The minimum reliability residue pairs may appear in different forms. The simplest form is a single pair; however, many residue pairs have different overlapping forms, with the same or opposite directions. The pair combinations are with some important features. First, one residue belongs to only one dual minimum residue pair, so the detected pairs should be removed from the reference residue sets, or they may disturb other pairs' detection. The quality weight on detected connecting lines should be zero according to the L 0 norm, namely, the boundaries act as superconductors in the reliability calculation. For the same direction, overlaps will enforce the connection. However, overlaps with opposite directions may weaken the original connections. Examples of overlapping pairs are shown in Fig. 3, where the connection between TD in (b) is doubled, and connection c 4 in (c) is eliminated with neutralization from c 5 .
As mentioned above, not all residues can be coupled in pairs in once reliability calculation according to the definitions. Therefore, the detection of all pairs is performed stepwise in the form of dual expanding. When coupled residues are removed from the reference sets, reliability values need to be updated for new dual pairs' detection. The procedure runs iteratively until there is no more detectable pairs and the reference sets do not change. The newly detected pairs may change existing ones; in any case, they always keep residue balanced. We call such balanced, connected residues as MBTs.

Path Priority for Path Integration
Under ideal conditions, all residues are included in detected dual pairs for a PU problem. However, there are some situations in which this is not the case; for instance, the residues are outside the defined domain. We regard all points on the border as reliability references, positive and negative. The definitions of positive reliability and negative reliability are renewed as E Q -T A R G E T ; t e m p : i n t r a l i n k -; s e c 2 . 3 ; 1 1 6 ; 6 1 5 According to the new reliability definitions, the dual reliability of a connecting line between an inside residue and the border may be not constant, and the connections cannot be determined directly by pair detection. In this case, we use a reliability comparison strategy to determine the final connection of dual residues. A typical example of the connection between residues and a border is shown in Fig. 4, where residues A and B cannot make up a detectable dual pair. It is clear that connecting line c 2 in (b) and c 1 would be broken first.
The reliability comparison is carried out during path integrating for unwrapped phase. We take the minimal reliability on given path s as its reliability index HðsÞ, namely E Q -T A R G E T ; t e m p : i n t r a l i n k -; s e c 2 . 3 ; 1 1 6 ; 4 8 6 HðsÞ ¼ minfPðu; vÞjðu; vÞ ∈ sg: Too low-path reliability suggests that the path passes through an unreliable area that may contain a phase discontinuity. Among the different paths, the most reliable one l from reference ðx r ; y r Þ to target ðx; yÞ should be E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 5 ; 1 1 6 ; 4 1 8 lðx; yÞ ¼ arg max s fHðsÞjs∶ðx r ; y r Þ → ðx; yÞg: According to this definition, we take the reliability index of the most reliable integration path as priority F at ðx; yÞ, which determines the integration order during unwrapping; namely E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 6 ; 1 1 6 ; 3 5 6 Fðx; yÞ ¼ HðlÞ ¼ minfPðu; vÞjðu; vÞ ∈ lg: Correspondingly, the priority based on reliability provides a quantized description of the most reliable integration path, which should guide the unwrapping to obtain the most reliable result.

Interface Propagation Method for Reliability and Priority
The most important thing in PU is to find the most reliable integration paths, which are based on reliable discrete sampling points. A demonstration of regular grid frame used in the dense 2-D PU is shown in Fig. 5(c), where the circles, termed centre points, denote phase sampling points and the squares, termed corner points, lie between them. The edges connecting the corners are potential discontinuity boundaries. Each corner point is with reliability values p þ and p − , dual reliability value p, reference residue labels s þ and s − for dual pair detecting, forerunner neighbor labels n þ and n − for connection line tracing, and priority f. Each corner edge has two properties the quality weight q and difference d between the forward and reverse connecting lines of the residue pairs. If d is nonzero, the edge will be a superconductor with zero-valued q during forward propagation and a barrier during path integration. Central connecting edges are potential integration paths. Each centre point has wrapped phase ψ, integration priority f i and unwrapped phase φ, and each centre edge has wrapped phase difference Δψ, which is used in calculating φ and q.
According to the definition, positive reliability value p þ and negative reliability value p − at the reference residues' position are both zero, but it is not convenient to get other values directly from the definitions. The relation between p þ or p − and quality weight q at ðx; yÞ can be rewritten as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 7 ; 1 1 6 ; 2 0 5 jΔp þ ðx; yÞj ¼ jΔp − ðx; yÞj ¼ qðx; yÞ; where j • j is the absolute value operator and Δ is the difference operator. This expression is in the form of a classic eikonal equation. 23 In other words, the reliability calculation is similar to wave-front propagation from sources with a given speed field, and the reliability corresponds to the travel time from the nearest reference residue. The procedure of solving eikonal equations is an implementation of interface propagation from reference residues before detecting and expanding of dual minimum reliability pairs, and flood-filling is employed for the backward interface propagation on the inverse dual reliability map. If all reference residues are correctly labeled and the definition of quality weight is specified, positive and negative reliability values at any position can be calculated using an eikonal equation solver. There are many candidates for defining the positive quality description function, including the coherence coefficient and the derivative variance. Here we use definitions that include the variance of the phase gradient in the horizontal and vertical directions: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 8 ; 1 1 6 ; 6 9 9 where σ 2 x and σ 2 y are the variances calculated in the x and y directions within a window and α is a constant positive weight.
According to Eq. (7), the reliability value p þ or p − at grid position ði; jÞ is calculated for its definite neighboring reliability values using an upwind eikonal equation solver as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 9 ; 1 1 6 ; 6 0 7 where q xAE and q yAE denote the weight values in the propagation directions, and the distance between neighbors is taken as a unit value in the case of grid sampling. During backward propagation of the flood-filling, the integration priority value f at ði; jÞ is calculated by comparison on neighboring priority values and its reliability value, as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 0 ; 1 1 6 ; 5 2 5 f i;j ¼ min½p i;j ; maxðf iAE1;j ; f i;jAE1 Þ:

Interface Propagating and MBTs Constructing in Parallel
In our proposed parallel scheme, a large data block is divided into small blocks, and each block or pixel has three designed statuses "unvisited," "active," and "converged." Status unvisited corresponds to unprocessed blocks or pixels, active is for blocks during propagation, and converged refers to finished calculation. Initially, reference residues for forward propagation and reference unwrapped phase pixel for backward propagation are labeled converged; the blocks and neighbor pixels of reference residues for forward propagation and blocks and neighbor pixels of the reference unwrapped phase pixel for backward propagation are labeled active; the others are labeled unvisited. Operations executed on active pixels of active blocks may change their status, and they run in loops until all blocks are converged. An example of the propagation between pixels and blocks is shown in Figs. 5(a), 5(b), and 5(d), where units in the dark area are converged with known reliability/priority values, those in the white area are unvisited, and those in gray are active. The calculation of the reliability/priority of active units is performed according to Eq. (9) or Eq. (10). Once the calculating finished, active units will be converted to be converged, and new active units will be checked and labeled. When calculation on pixel a in (a) is finished, the active pixel interface g-a-f may move to g-d-e-f. Similarly, when the operation on block F in (d) is finished, active block interface A-B-C-D-E-F-G may move to A-B-C-D-E-J-G. This process is similar to the framework of the fast iterative method, 24,25 a well-known parallel version of the fast marching method, 26 which is usually employed as an eikonal equation solver.
After the positive and negative forward propagations are complete, the obtained positive and negative reference residue maps s þ and s − are used to detect dual pairs for constructing MBTs. If a residue pair is detected, the positive or negative forerunner maps n þ ∕n − are implemented to find the corresponding connecting edges. The value of d for this should be renewed at the same time. When all the detectable pairs have been found, the pixels whose reliability, reference residues, and forerunners derived from them are renewed by forward propagation. The forward propagation and residue pair detection run alternately until no detectable residue pairs remain. The final dual reliability map p is the sum of p þ and p − , namely E Q -T A R G E T ; t e m p : i n t r a l i n k -; s e c 3 . 2 ; 1 1 6 ; 1 3 8 where ðm; nÞ refers to the unwrapped neighbor with the greatest p, and Δψ refers to their wrapped difference.
Algorithm 1 describes the parallel iterative algorithm for forward propagation, whose time complexity is approximately Oðm × tÞ, where t is the number of active times.
The parallel iterative algorithm for backward propagation is described in Algorithm 2, whose time complexity is approximately Oðm × tÞ.
From the backward propagation procedure, we obtain an integration path reliability map, and the unwrapping process runs according to the integration paths with the highest priority. An example of MBT construction and unwrapping using a priority map is shown in Fig. 6. The original wrapped phase and residues are marked in part (a), where reference residues are marked Algorithm 1 Parallel forward propagation. 1: Detect residues and initialize s þ , s − , p þ , p − , n þ , n − , d .

Experiments and Analysis
The machine used for the unwrapping process in the experiments was equipped with 2 CPUs (12 cores, 1.9 GHz), a GPU (15 multiprocessors, 2880 cores, 0.75 GHz) with 12 GB of device memory (3 004 MHz), a 32 GB host memory (dual-channel, DDR4, 2 133 MHz), and a 64-bit operating system. The experimental data are from interferograms of real airborne SAR, and they are processed in small blocks in parallel.

PU Experiment on Data A
The first real wrapped phase data block was with a size of 16;384 × 10;928 pixels. The data contained a gentle variation overall with serious noise at many locations. The wrapped phase was shown in Fig. 7, which is not a full-resolution visualization for its huge size. The data were decomposed into blocks of 32 × 32 pixels, giving 175,446 blocks in total. First, we detected the residues in the data, and the result showed that there are 259,053 positive and 259,065 negative residues. Then the weight maps in the horizontal and vertical directions were calculated according to Eq. (8). The quality weight maps for the horizontal direction were shown as logarithmic values in Fig. 8, wherein the noisy areas are labeled with low-weight values. With detected residues and weight maps, the reliability can be calculated. Here we provide the initial positive reliability map in Fig. 9 and the corresponding positive reference map in Fig. 10. It is clear that high-reliability values surround noisy areas in the reliability map. Then the residue pair detection and reliability map updating were performed iteratively for 14 times, and after that, there were no detectable residue pairs. The final dual reliability map is shown in Fig. 11, and the final positive reference map in Fig. 12. Compared to initial maps, it is obvious that there are only fewer residues left, and a lot of them are located near borders and bring about strips of reference regions as shown. A priority map obtained after backward propagation is shown in Fig. 13. The difference between the reliability map and the priority map can be easily found in the marked area after flood-filling. With the priority map, a final result of unwrapping is shown in Fig. 14. This is an ideal unwrapping result. All discontinuous disturbances and residues   did not have bad impacts on the result. They only affect their locale unwrapping, which guarantees overall quality.

PU Experiment on Data B
For further comparison, we took another real but more challenging airborne data. The second data block is the same size of 16;384 × 10;928 pixels and contains more noises areas. As shown in Fig. 15, there are lots of noise at the left and the right, and streak noises at the bottom, which bring about 476,291 positive and 476,456 negative residues found in residue detection. According to the processing flow, we carried out the weight calculation, the iterative reliability map calculation and the dual pair detection, flood-filling for the priority map, and final unwrapping with the priority map obtained. The initial positive reliability map is shown in Fig. 16 and the final dual reliability map in Fig. 17, and it can be found that the reliability values in the central area are significantly higher than in other areas. The final unwrapped phase is shown in Fig. 18(a). Unwrapping in most of the visually recognizable continuous areas performed well and provided a good result of phase recovery.
For comparison, we also employed several other existing PU tools to process data B. The methods involved include the quality-guided flood method, 27 MCF method, 14,28 and SNAPHU. 14 For the quality-guided flood method, the proposed reliability map is used as a weight map, and  Fig. 16 The initial positive reliability map of data B. the corresponding result is shown in Fig. 18(b). It is obvious that the flood method cannot guarantee a reliable unwrapped result in the case of heavy noise, where wrong unwrapped phase appear in the large areas marked with ellipses. The same problem exists in MCF's results shown in Fig. 18(c), and sharper edges around noise indicate more frequent phase discontinuities, although the graph-based PU method did better to some extent, compared to the quality-guided method. SNAPHU with iterative optimization is very time-consuming for large-scale data with heavy noise, so the tool ran in a tiling parallel mode. Data B are divided into 16 × 11 tiles, and each tile is of about 1000 × 1000 pixels with 50-pixel wide overlap bands on its borders. The SNAPHU running costed 19 min 9.616 s, and its unwrapped result is shown in Fig. 18(d). The main problem is that the consistency between tiles cannot be guaranteed. Also we carried out some quantitative analysis of the experiments, including the number of dual pairs, the length of phase discontinuity boundaries and the running. There were 14 iterations in experiment A and 16 iterations in experiment B for dual pair detecting, and the number of pairs found of each iteration is shown in Fig. 19. We can see that the number of detectable pairs in each iteration is exponentially decreasing; namely, the MBT constructing on million-level residues can be finished after a dozen iterations.
As for the length of phase discontinuity boundaries, the L 0 norm and L 1 norm of the unwrapped phase φ can be expressed as x < 0.5 1; x ≥ 0.5 for L 0 ½x; x ≥ 0.5 for L 1 ; and j • j is an absolute value operator and ½• is a rounding operator. It should be noted that the phase cycle is normalized to 1 here for convenience. We made statistics on discontinuous lengths and processing times, and the results were listed in Table 1. It is apparent according to this table that the MBTs contributed most of the discontinuous boundaries, and the percentages are 98.75% and 98.69% for experiment A, 86.83% and 85.52% for experiment B. Considering the increasing need for high efficiency, the run time is also a key factor. Experiment A takes 4 min 48.860 s, and experiment B takes 7 min 38.648 s with more residues. Although the running time is related to the number of residues and the number of iterations, the statistical data show that it is not worse than a linear relationship.

Conclusions
A parallel PU method for large-scale data based on local minimum reliability dual expanding is proposed in this paper. First, we define dual reliability for phase data with reference residues and quality weight maps, a definition of dual minimum reliability residue pairs based on dual reliability for discontinuity boundaries. The detection and expansion of the dual pairs construct the MBTs, and other discontinuity boundaries are implied in the priority map derived from floodfilling on the inverse dual reliability map. The final unwrapped result is produced by path integration based on the priority map and MBTs. We also provide a parallel interface propagation scheme for reliability calculation and priority calculation, which bring about processing efficiency. The PU method thus provides a robust and efficient choice for large-scale data.

Acknowledgments
This research was funded by the Key R&D Program Projects in Hainan Province (No. ZDYF2019008) and the Strategic Priority Research Program of the Chinese Academy  Fig. 19 The number of dual residue pairs found in each iteration.