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

Jian Gao; Zhongchang Sun

doi:10.1117/1.JRS.13.038506

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

Jian Gao, Zhongchang Sun

Author Affiliations +

Journal of Applied Remote Sensing, Vol. 13, Issue 3, 038506 (September 2019). https://doi.org/10.1117/1.JRS.13.038506

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.

1. 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)¹ 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.²^–⁴ 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,⁵ priority queuing,⁶^,⁷ and quality mapping with edges and shadow areas.⁸ Another implemented form of them, tiling strategies⁹^,¹⁰ 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 branch-cuts methods. The classic branch-cut method¹¹ 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,¹² the graph-cuts method,¹³ and the statistical-cost network flow algorithm (SNAPHU).¹⁴^,¹⁵ 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,¹⁴ the branch-cut,¹⁶ and quality-guided¹⁷ methods; irregular data tiling is also applied in the MCF method,¹⁸^,¹⁹ which uses subgraphs involving partial residues to reduce the workload and in the extended minimum cost flow method,²⁰ 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²¹ and graphics processing unit (GPU) region-growing²² 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.

2. Principle and Method

In the 2-D definition domain, the unwrapped phase $Φ$ at $(x, y)$ can be given as

Φ (x, y) = Φ_{0} + \int_{l} W [\nabla Ψ (u, v)] d s,

where

\nabla Ψ

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).

Fig. 1

(a) Discontinuity boundaries and residues, and (b) they are enveloped in low-reliability areas of the constructed dual reliability map.

2.1.

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

Eq. (1)

P_{S} (x, y) = \int_{c} Q (u, v) d s,

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

P^{+} (x, y) = \min {P_{S} (x, y) | S \in R^{+}} P^{-} (x, y) = \min {P_{S} (x, y) | S \in R^{-}},

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:

Eq. (2)

P (x, y) = P^{+} (x, y) + P^{-} (x, y) .

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.

2.2.

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

Eq. (3)

c (S_{a}, S_{b}) = \underset{c}{\arg \min} {\int_{c} Q (u, v) d s | c : S_{a} \to S_{b}} .

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.

Fig. 2

(a) A minimum reliability residue pair, (b) a disturbed residue pair, and (c) their dual reliability, where the areas encircled by solid lines expand from the inside positive residues and the ones by dashed lines from negatives.

If all residue pairs are detected, and the connecting lines are determined according to Eq. (3), the obtained phase discontinuity boundaries would determine the final unwrapping result. A PU optimization problem based on minimum reliability residue pairs can be written as

Eq. (4)

\min \int_{D} Q (u, v) d s = \min \sum_{D} P,

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}$ .

Fig. 3

(a) A single residue pair, (b) two overlapped residue pairs with the same direction, and (c) two overlapped residue pairs with different directions.

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.

2.3.

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

P^{+} (x, y) = \min {P_{S} (x, y) | S \in R^{+} \cup Ω} P^{-} (x, y) = \min {P_{S} (x, y) | S \in R^{-} \cup Ω} .

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.

Fig. 4

(a) Connections between residues and border and (b), (c) reliability chart for two cases.

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

H (s) = \min {P (u, v) | (u, v) \in s} .

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

Eq. (5)

l (x, y) = \underset{s}{\arg \max} {H (s) | s : (x_{r}, y_{r}) \to (x, y)} .

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

Eq. (6)

F (x, y) = H (l) = \min {P (u, v) | (u, v) \in l} .

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.

3. Parallel Numerical Implementation

3.1.

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$ .

Fig. 5

Propagation between (a) points and (b) blocks, (c) centre points are separated by a connecting line of corner points, and (d) the block interface moving as status changes. The dark area refers to processed points or blocks, the gray area to interface points or blocks, and the white area to uncovered points or blocks.

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

Eq. (7)

| Δ p^{+} (x, y) | = | Δ p^{-} (x, y) | = q (x, y),

where

| • |

is the absolute value operator and

Δ

is the difference operator. This expression is in the form of a classic eikonal equation. ²³ 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:

Eq. (8)

q_{x} = \frac{1}{σ_{x}^{2} + α}, q_{y} = \frac{1}{σ_{y}^{2} + α},

where

σ_{x}^{2}

and

σ_{y}^{2}

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

Eq. (9)

p_{i, j}^{+} = \min (p_{i \pm 1, j}^{+} + q_{x \pm}, p_{i, j \pm 1}^{+} + q_{y \pm}), p_{i, j}^{-} = \min (p_{i \pm 1, j}^{-} + q_{x \pm}, p_{i, j \pm 1}^{-} + q_{y \pm}),

where

q_{x \pm}

and

q_{y \pm}

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

Eq. (10)

f_{i, j} = \min [p_{i, j}, \max (f_{i \pm 1, j}, f_{i, j \pm 1})] .

3.2.

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,²⁴^,²⁵ a well-known parallel version of the fast marching method,²⁶ 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

φ_{i, j} = φ_{m, n} + Δ ψ,

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 \times t)$ , where $t$ is the number of active times.

Algorithm 1

Parallel forward propagation.

1: Detect residues and initialize

s^{+}

,

s^{-}

,

p^{+}

,

p^{-}

,

n^{+}

,

n^{-}

,

d

.

2: Searching active blocks and make list

L

3: Repeat

4: While active block list

L

is not empty do

5: For all block

b \in L

parallel do

6: For iteration number

i = 0

to

m

do

7: Calculate and update

(p, d, s, n)

on

b

8: End for

9: If

b

converges then

10: Label b converged

11: For all neighbors

b_{n}

of

b

do

12: Calculate and update

(p, d, s, n)

on borders of

b_{n}

13: If

b_{n}

changes then

14: Label

b_{n}

active

15: Insert

b_{n}

into

L

16: End if

17: End for

18: End if

19: End for

20: For all

b \in L

parallel do

21: If

b

is converged then

22: Remove

b

from

L

23: End if

24: End for

25: End while

26: Search for dual pairs with

s^{+}

,

s^{-}

and make list

T

27: For all block

s \in T

parallel do

28: Search connecting lines with

n^{+}

or

n^{-}

, and update

d

29: End for

30: Make active block list

L

with newly coupled residues

31: Until list

T

is empty

The parallel iterative algorithm for backward propagation is described in Algorithm 2, whose time complexity is approximately $O (m \times t)$ .

Algorithm 2

Parallel backward propagation.

1: Appoint the reference point and label active block into list

L

2: While active block list

L

is not empty do

3: For all block

b \in L

parallel do

4: For iteration number

i = 0

to

m

do

5: Calculating and updating f on

b

6: End for

7: If

b

converges then

8: Label b converged

9: For all neighbor

b_{n}

of

b

do

10: Calculating and updating

b_{n}

11: If

b_{n}

changes then

12: Label

b_{n}

active

13: Insert

b_{n}

into

L

14: End if

15: End for

16: End if

17: End for

18: For all

b \in L

parallel do

19: If

b

is converged then

20: Delete

b

from

L

21: End if

22: End for

23: End while

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 with circles for the positives and with squares for the negatives. The final reliability map and detected MBTs are shown in (f), and the priority map is shown in part (g), where a unwrapped reference phase is located at the upper left. The unwrapped phase based on a certain priority level is shown in part (h), and the final unwrapped phase is shown in part (i).

Fig. 6

(a) Wrapped spiral phase and residues, weight map on (b) the horizontal direction and (c) the vertical direction, initial (d) positive reliability map and (e) negative reliability map, (f) final dual reliability map and with MBTs detected, (g) priority map, (h) unwrapped phase above a certain priority level, and (i) the final unwrapped phase.

4. 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.

4.1.

PU Experiment on Data A

The first real wrapped phase data block was with a size of $16,384 \times 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 \times 32 pixels$ , giving 175,446 blocks in total.

Fig. 7

The wrapped phase with 518,118 residues from data A.

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.

Fig. 8

The horizontal quality weight map of data A.

Fig. 9

The initial positive reliability map of data A.

Fig. 10

The initial positive reference map of data A.

Fig. 11

The final dual reliability map of data A.

Fig. 12

The final reference map of data A.

Fig. 13

The final priority map of data A.

Fig. 14

The final unwrapped phase of data A.

4.2.

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 \times 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.

Fig. 15

The wrapped phase with 952,747 residues of data B.

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.

Fig. 16

The initial positive reliability map of data B.

Fig. 17

The final dual reliability map of data B.

Fig. 18

The unwrapped phase of data B from (a) the proposed method, (b) quality-guided flood method, (c) MCF method, and (d) SNAPHU.

For comparison, we also employed several other existing PU tools to process data B. The methods involved include the quality-guided flood method,²⁷ MCF method,¹⁴^,²⁸ and SNAPHU.¹⁴ For the quality-guided flood method, the proposed reliability map is used as a weight map, and 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 \times 11$ tiles, and each tile is of about $1000 \times 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.

Fig. 19

The number of dual residue pairs found in each iteration.

As for the length of phase discontinuity boundaries, the $L^{0}$ norm and $L^{1}$ norm of the unwrapped phase $φ$ can be expressed as

\sum_{i = 1}^{M - 1} \sum_{j = 1}^{N - 1} k (| φ_{i + 1, j} - φ_{i, j} |) + k (| φ_{i, j + 1} - φ_{i, j} |),

where

k (x) = {\begin{cases} 0, & x < 0.5 \\ 1, & x \geq 0.5 for L^{0} \\ [x], & x \geq 0.5 for L^{1} \end{cases},

and

| • |

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.

Table 1

Discontinuity lengths and running times from the unwrapping experiments.

	L0 nrom			L1 nrom			Time
	MBT (%)	F-F (%)	Total (pix)	MBT (%)	F-F (%)	Total (pix)	Time
Exp. A	98.75	1.25	1,031,054	98.69	1.31	1,070,106	4 min 48.860 s
Exp. B	86.83	13.17	3,153,762	85.52	14.48	3,398,380	7 min 38.648 s

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.

5. 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 flood-filling 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 of Sciences (Grant No. XDA19030104), Scientific Research Foundation of Nanjing University of Posts and Telecommunications (Grant No. NY215109), and the National Natural Science Foundation of China (Grant Nos. 41501378 and 41501497). The authors declare no conflicts of interest.

References

1.

D. C. Ghiglia and M. D. Pritt, Two-Dimensional Phase Unwrapping: Theory, Algorithms and Software, Wiley, Blackwell, New York (1998). Google Scholar

2.

J. Gao, J. Li and L. Shi, “Two-dimensional phase unwrapping method using cost function of L0 norm,” IOP Conf. Ser.: Earth Environ. Sci., 57 (1), 12 –41 (2017). https://doi.org/10.1088/1755-1315/57/1/012041 Google Scholar

3.

M. Wang et al., “Robust wrapping-free phase retrieval method based on weighted least squares method,” Opt. Lasers Eng., 97 34 –40 (2017). https://doi.org/10.1016/j.optlaseng.2017.05.008 Google Scholar

4.

X. Wang, S. Fang and X. Zhu, “Weighted least-squares phase unwrapping algorithm based on a non-interfering image of an object,” Appl. Opt., 56 (15), 4543 –4550 (2017). https://doi.org/10.1364/AO.56.004543 APOPAI 0003-6935 Google Scholar

5.

M. Zhao and Q. Kemao, “Quality-guided phase unwrapping implementation: an improved indexed interwoven linked list,” Appl. Opt., 53 (16), 3492 –3500 (2014). https://doi.org/10.1364/AO.53.003492 APOPAI 0003-6935 Google Scholar

6.

H. Zhong, S. Zhang and J. Tang, “Path following algorithm for phase unwrapping based on priority queue and quantized quality map,” in Int. Conf. Computat. Intell. and Software Eng., 1 –4 (2009). https://doi.org/10.1109/CISE.2009.5364436 Google Scholar

7.

H. Zhong et al., “An improved quality-guided phase-unwrapping algorithm based on priority queue,” IEEE Geosci. Remote Sens. Lett., 8 364 –368 (2011). https://doi.org/10.1109/LGRS.2010.2076362 Google Scholar

8.

K. Chen et al., “An object image edge detection based quality-guided phase unwrapping approach for fast three-dimensional measurement,” in IEEE/ASME Int. Conf. Adv. Intell. Mechatron., 571 –576 (2013). https://doi.org/10.1109/AIM.2013.6584153 Google Scholar

9.

H. Yu and H. Lee, “A convex hull algorithm based fast large-scale two-dimensional phase unwrapping method,” in IEEE Int. Geosci. and Remote Sens. Symp. (IGARSS), 3824 –3827 (2017). https://doi.org/10.1109/IGARSS.2017.8127834 Google Scholar

10.

H. Yu et al., “Large-scale L0-norm and L1-norm 2-D phase unwrapping,” IEEE Trans. Geosci. Remote Sens., 55 (8), 4712 –4728 (2017). https://doi.org/10.1109/TGRS.2017.2698452 IGRSD2 0196-2892 Google Scholar

11.

R. M. Goldstein, H. A. Zebker and C. L. Werner, “Satellite radar interferometry: two-dimensional phase unwrapping,” Radio Sci., 23 (4), 713 –720 (1988). https://doi.org/10.1029/RS023i004p00713 Google Scholar

12.

M. Costantini, “A novel phase unwrapping method based on network programming,” IEEE Trans. Geosci. Remote Sens., 36 (3), 813 –821 (1998). https://doi.org/10.1109/36.673674 IGRSD2 0196-2892 Google Scholar

13.

J. M. Bioucas-Dias and G. Valadao, “Phase unwrapping via graph cuts,” IEEE Trans. Image Process., 16 (3), 698 –709 (2007). https://doi.org/10.1109/TIP.2006.888351 IIPRE4 1057-7149 Google Scholar

14.

C. W. Chen and H. A. Zebker, “Two-dimensional phase unwrapping with statistical models for nonlinear optimization,” in IEEE Int. Geosci. and Remote Sens. Symp., 3213 –3215 (2000). https://doi.org/10.1109/IGARSS.2000.860386 Google Scholar

15.

X. Feng et al., “A new method about placement of the branch cut in two-dimensional phase unwrapping,” in 1st Asian and Pac. Conf. Synth. Aperture Radar, 755 –759 (2007). https://doi.org/10.1109/APSAR.2007.4418721 Google Scholar

16.

Q. Huang et al., “Parallel branch-cut algorithm based on simulated annealing for large-scale phase unwrapping,” IEEE Trans. Geosci. Remote Sens., 53 (7), 3833 –3846 (2015). https://doi.org/10.1109/TGRS.2014.2385482 IGRSD2 0196-2892 Google Scholar

17.

H. Zhong et al., “Parallel quality-guided phase unwrapping algorithm in shared memory environment,” Geomatics Inf. Sci. Wuhan Univ., 42 (1), 130 (2017). https://doi.org/10.13203/j.whugis20140728 Google Scholar

18.

M. Even, “Accelerating phase unwrapping based on integer linear programming by processing of subgraphs,” in IEEE Int. Geosci. and Remote Sens. Symp. (IGARSS), 6441 –6444 (2016). https://doi.org/10.1109/IGARSS.2016.7730683 Google Scholar

19.

Z. Wang et al., “A novel fast phase unwrapping method for large interferometric datasets,” in IEEE Int. Geosci. and Remote Sens. Symp. (IGARSS), 6445 –6448 (2016). Google Scholar

20.

A. Pepe et al., “New advances of the extended minimum cost flow phase unwrapping algorithm for SBAS-DInSAR analysis at full spatial resolution,” IEEE Trans. Geosci. Remote Sens., 49 (10), 4062 –4079 (2011). https://doi.org/10.1109/TGRS.2011.2135371 IGRSD2 0196-2892 Google Scholar

21.

H. Chen et al., “An efficient FPGA-based parallel phase unwrapping hardware architecture,” IEEE Trans. Comput. Imaging, 3 (4), 996 –1007 (2017). https://doi.org/10.1109/TCI.2017.2663767 Google Scholar

22.

S. E. Popov, “Improved phase unwrapping algorithm based on NVIDIA CUDA,” Programm. Comput. Software, 43 (1), 24 –36 (2017). https://doi.org/10.1134/S0361768817010054 PCSODA 0361-7688 Google Scholar

23.

V. I. Arnold, Lectures on Partial Differential Equations, Springer-Verlag, New York (2004). Google Scholar

24.

A. Capozzoli et al., “A comparison of fast marching, fast sweeping and fast iterative methods for the solution of the Eikonal equation,” in 21st Telecommun. Forum Telfor (TELFOR), 685 –688 (2013). https://doi.org/10.1109/TELFOR.2013.6716321 Google Scholar

25.

W. Jeong and R. T. Whitaker, “A fast iterative method for Eikonal equations,” SIAM J. Sci. Comput., 30 (5), 2512 –2534 (2008). https://doi.org/10.1137/060670298 SJOCE3 1064-8275 Google Scholar

26.

J. A. Sethian, Level Set Methods and Fast Marching Methods. Evolving Interfaces in Computational Geometry, Fluid Mechanics, Computer Vision, and Materials Science, Cambridge University Press, Cambridge (1999). Google Scholar

27.

K. Itoh, “Analysis of the phase unwrapping algorithm,” Appl. Opt., 21 (21), 2470 (1982). https://doi.org/10.1364/AO.21.002470 APOPAI 0003-6935 Google Scholar

28.

A. V. Goldberg, “An efficient implementation of a scaling minimum-cost flow algorithm,” J. Algorithms, 22 1 –29 (1997). https://doi.org/10.1006/jagm.1995.0805 JOALDV 0196-6774 Google Scholar

Biography

Jian Gao received his PhD in photogrammetry and remote sensing from the State Key Laboratory of Information Engineering in Surveying, Mapping, and Remote Sensing, Wuhan University, Wuhan, China, in 2012. He is currently a lecturer in the Department of Geographic Information, Nanjing University of Posts and Telecommunications, Nanjing, China. His current research interests include radar interferometry, phase unwrapping, and remote sensing image processing.

Zhongchang Sun received his PhD in cartography and geographic information systems from the Centre for Earth Observation and Digital Earth, Chinese Academy of Sciences, Beijing, China, in 2011. He is currently an associate professor at the Institute of Remote Sensing and Digital Earth of the Chinese Academy of Sciences. He has authored or co-authored more than 20 SCI papers. His current research interests include urban environment remote sensing, land surface dynamics remote sensing, and radar remote sensing.

CC BY: © The Authors. Published by SPIE under a Creative Commons Attribution 4.0 Unported License. Distribution or reproduction of this work in whole or in part requires full attribution of the original publication, including its DOI.

Citation Download Citation

Jian Gao and Zhongchang Sun "Phase unwrapping method based on parallel local minimum reliability dual expanding for large-scale data," Journal of Applied Remote Sensing 13(3), 038506 (12 September 2019). https://doi.org/10.1117/1.JRS.13.038506

Received: 14 February 2019; Accepted: 13 August 2019; Published: 12 September 2019

Access the abstract

JOURNAL ARTICLE
18 PAGES

DOWNLOAD PAPER SAVE TO MY LIBRARY

GET CITATION

CITATIONS

Cited by 5 scholarly publications and 1 patent.

Explore citations on Lens.org

KEYWORDS

Reliability

Sun

Interfaces

Interferometric synthetic aperture radar

Floods

Data processing

Superconductors

1.

Introduction

2.

Principle and Method

Fig. 1

2.1.

The Dual Reliability Definition Based on Phase Data Quality

Eq. (1)

Eq. (2)

2.2.

Minimum Reliability Residue Pairs and Local Minimum Reliability Dual Expanding

Eq. (3)

Fig. 2

Eq. (4)

Fig. 3

2.3.

Path Priority for Path Integration

Fig. 4

Eq. (5)

Eq. (6)

3.

Parallel Numerical Implementation

3.1.

Interface Propagation Method for Reliability and Priority

Fig. 5

Eq. (7)

Eq. (8)

Eq. (9)

Eq. (10)

3.2.

Interface Propagating and MBTs Constructing in Parallel

Algorithm 1

Algorithm 2

Fig. 6

4.

Experiments and Analysis

4.1.

PU Experiment on Data A

Fig. 7

Fig. 8

Fig. 9

Fig. 10

Fig. 11

Fig. 12

Fig. 13

Fig. 14

4.2.

PU Experiment on Data B

Fig. 15

Fig. 16

Fig. 17

Fig. 18

Fig. 19

Table 1

5.

Conclusions

Acknowledgments

References

Biography

Show All Keywords

Keywords/Phrases

Search In:

Publication Years