Location via proxy:   [ UP ]  
[Report a bug]   [Manage cookies]                
Next Article in Journal
Land Surface Greening and CO2 Fertilization More than Offset the Gross Carbon Sequestration Decline Caused by Land Cover Change and the Enhanced Vapour Pressure Deficit in Europe
Previous Article in Journal
Identification and Analysis of Heatwave Events Considering Temporal Continuity and Spatial Dynamics
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

A Parallel InSAR Phase Unwrapping Method Based on Separated Continuous Regions

1
School of Geographic and Biologic Information, Nanjing University of Posts and Telecommunications, Nanjing 210023, China
2
Department of Geomatics Engineering, University of Calgary, Calgary, AL T2N 1N4, Canada
3
School of Civil Engineering, Anhui Jianzhu University, Hefei 230601, China
4
Key Laboratory of Digital Earth Science, Aerospace Information Research Institute, Chinese Academy of Sciences, Beijing 100094, China
5
School of Marine Technology and Geomatics, Jiangsu Ocean University, Lianyungang 222005, China
*
Author to whom correspondence should be addressed.
Remote Sens. 2023, 15(5), 1370; https://doi.org/10.3390/rs15051370
Submission received: 18 January 2023 / Revised: 24 February 2023 / Accepted: 27 February 2023 / Published: 28 February 2023

Abstract

:
Phase unwrapping is an imperative step in interferometry processing that has a significant influence on the quality of subsequent products. Many existing phase unwrapping algorithms have been designed to solve for the unwrapped phase under the assumption that noisy areas with discontinuities are small or that reliable continuity can be recovered there. They attempt to restore the unwrapped phase by using continuity and data quality measures, such as residues. However, when the observing field is divided into separate zones of continuous phase due to a large range of noise, such as those caused by rivers or mountains, it is difficult to use traditional phase unwrapping techniques to recover global continuity in these noisy areas. To address this challenge, we present a two-dimensional parallel phase unwrapping method that is designed to handle cases where the continuity of the phase is separated by closed noisy loops. Based on continuity distances, this method aims to identify continuous regions that are free of hidden phase discontinuities and restore phase continuity between the separated regions. A heterogeneous residual diffusion scheme is used to restore the unwrapped phase outside continuous regions. The parallel algorithm for extracting continuous regions, restoring continuity between the regions, and diffusing residuals was implemented on a GPU device to increase the processing efficiency. We applied our method to typical TanDEM-X data covering rivers, islands, and mountains and demonstrated that it is a promising solution for large-scale, heavily noisy phase unwrapping problems.

1. Introduction

The phase unwrapping (PU) problem is a fundamental challenge in the processing of interferometric synthetic aperture radar (InSAR) data. A large body of research has been devoted to developing various solutions to this problem, employing a range of strategies and techniques. The central issue addressed by these efforts is how to accurately restore the unwrapped phase data, while minimizing the impact of phase discontinuities. A key requirement for most PU methods is that the phase difference between adjacent samples should not exceed half a cycle within an interferogram [1,2]. However, achieving this goal can be challenging. Phase discontinuities can arise for a variety of reasons, including noise, errors in the interferometric phase measurement, and the presence of topographic features that introduce phase jumps. Therefore, finding ways to effectively cope with these phase discontinuities is a key focus of PU research.
One class of methods for unwrapping phase data in interferometry is to design an index map that represents the quality of the wrapped phase data. This index map is expected to approximate the probability distribution of discontinuities in the data. A high-quality unwrapped seed is then set in the interferogram, and the unwrapped area is gradually expanded using a certain growing algorithm, where the expansion is guided by the quality index. To handle cases where multiple high-quality seeds are present, a region-growing scheme can be utilized, where neighboring regions are merged once they have been grown [3,4]. The continuity between regions can be determined by comparing the phase difference between connecting boundaries. If the difference is less than half a cycle, the two regions are considered to have reliable neighboring continuity and can be easily merged. However, if neighboring continuity cannot be guaranteed, additional strategies may be necessary. One such strategy is to use discrete reference points to form an irregular grid [5,6]; another strategy is to use a particle filter to select expansion paths that are guided by the quality of the data [7].
One of the key strategies for quality-guided PU methods is the explicit detection and determination of discontinuity boundaries. These boundaries are defined between the locations where the phase residue, which is the difference between the wrapped and unwrapped phase, changes abruptly. To construct a graph representing these discontinuity boundaries, researchers typically use principles such as minimizing the total length of the boundaries. Once the graph is constructed, the solution for unwrapping can be obtained by difference integrating away from the discontinuity boundaries. There are several different approaches to constructing this graph, including the branch-cut method [8,9], network flow methods such as SNAPHU [10], and the use of minimum spanning trees [11] and graph-cut [12] algorithms.
Hybrid unwrapping strategies have been proposed to improve the unwrapping of wrapped interferograms by dividing them into high- and low-quality areas. These methods utilize a combination of traditional PU techniques and more robust methods, such as curve fitting [13,14,15], the convex hull algorithm [16], or minimum norms [17], to handle areas of low-quality phase data. These hierarchical unwrapping schemes are dependent on the overall quality of the data. When the high-quality areas are completely separated by noise and the neighboring continuity cannot be guaranteed, their unwrapping performance is significantly hindered. In addition, there are methods to achieve PU processing with the help of a continuum filtering framework [18] or auxiliary information of observing terrain [19].
Recently, there has been an increasing trend to incorporate deep learning techniques into PU research. The principle of PU methods using deep learning is based mainly on the learning to solve the large number of parameters attached in neural networks from sample datasets and extract useful information in the wrapped phase data. For example, some researchers have used deep learning to generate quality maps [20], estimate the phase gradient [21,22,23,24], predict discontinuity boundaries [25,26], and perform image segmentation [27,28,29,30]. This has been shown to be very promising to increase the accuracy and efficiency of unwrapping.
However, many techniques for addressing phase discontinuities primarily rely on the assumption that the regions affected by noise are localized and small or that continuity information can be inferred from the wrapped phases present in the noisy data. When the interferometric phase contains large areas of noise that separate contiguous regions or when it is challenging to obtain reliable continuity information from the noisy areas, the PU methods based on this assumption often fail to produce accurate results. In this study, we propose a novel two-dimensional unwrapping method that utilizes continuous regions to address the issue of continuity separation.
The quality of data outside the separated continuous regions is inadequate, and ensuring continuity in these areas is uncertain. Based on the assumption of a certain level of continuity between adjacent regions, our method constructs a global continuity framework from the separated regions. This framework is then used to unwrap the phases located outside the continuous regions. This approach is particularly well suited for situations where high-quality data are obscured by obstacles such as radar shadows, water bodies, and forests. To improve the computational efficiency, we also implemented a GPU-based algorithm to parallelize the proposed method, making it well-suited to large-scale PU problems.
In Section 2, we provide a detailed explanation of the principle behind the proposed method. Its numerical implementation is discussed in Section 3, followed by an examination of the performance of the developed method through unwrapping experiments and a comparative analysis of three real interferometric phase data in Section 4. In Section 5, we engage in a thorough discussion of the proposed method. We conclude with a summary of our findings and suggestions for future research in Section 6.

2. Principle and Method

2.1. Defining and Extracting Continuous Regions

For a given interferogram, we can extract regions with a continuous phase based on specific quality indexes, for example the weighted distance to phase residues. Such a quality index is a quantitative description of the tendency to phase discontinuity, as shown in Figure 1, and has been used to guide path integration for coupling dual residues during phase unwrapping. Extraction of continuous regions is closely related to a quality measure [31], for convenience, which here is called the continuity distance D and is defined as
D ± ( x , y ) = min { c s Q ( u , v ) d ρ | s S ± }
D ( x , y ) = D + ( x , y ) + D ( x , y ) ,
where Q is the local phase quality descriptor, s refers to the reference residue point, c s refers to the shortest path to s, d ρ refers to the differential length on c s , and S ± refers to the positive or negative phase residue set.
In the case that the phase retains a certain continuity among neighboring regions, the appropriate use of the continuity may make it easier to connect and merge them. The definition of the continuous region Θ can be expressed as follows:
Θ = { ( x , y ) | D ( x , y ) γ , ( x , y ) Ω } ,
where Θ refers to a continuous region and connected domain, γ refers to the corresponding continuity distance threshold, and Ω refers to areas masked for very low coherence or others. The function of the threshold γ is to guarantee that no hidden phase discontinuity appears in continuous regions. It is worth noting that different regions might have different thresholds for better area coverage.
The continuous regions should be as large as possible once the thresholds have been set. The thresholds can be described as
( γ 1 , , γ n ) = arg max γ Ξ , Ξ = 1 n Θ i ,
where Ξ refers to the region set. The threshold values should be as small as possible, and they could be specified by searching from small to large in the continuity distance range.

2.2. Region Residue Pairs’ Detection and Region Phase Alignment

Setting local phase references for the separated continuous regions, we can obtain the continuous unwrapped phase in each region, as cycle jumps only exist outside. Taking into account cycle jumps between the phase references, it is necessary to align the phase of all regions for a global continuity frame. Phase connections must be established between them.

2.2.1. Definitions of Region Phase Difference and Region Residue

To establish phase connections between separated continuous regions, we duplicated and propagated unwrapped phase Φ within the continuous regions outward according to the principle of the minimum Euclidean distance. The mapping from a given point to its corresponding unwrapped phases within the continuous region can be described as
( x , y ) arg min ( u , v ) ( x u ) 2 + ( y v ) 2 , ( x , y ) Ξ , ( u , v ) Ξ ,
where ( x , y ) refers to the location of a point outside and ( u , v ) refers to a point within regions. With the rule, the phase outside of regions can be set, such as Φ ( x , y ) = Φ ( u , v ) , and the areas outside (white areas in Figure 2) are divided and attached to their nearest regions. Phase connections between continuous regions depend on the phase at expanded boundaries. We denote the expanded boundary as l, shown as dashed lines in Figure 2.
The mean and variance of the expanded phase difference are calculated between continuous regions on the boundaries. For the boundary l between regions M and N, the average phase difference Δ Φ l from M to N is calculated by path integrating as follows:
Δ Φ l = 1 ρ l ( Φ N Φ M ) d ρ ,
where d ρ is the differential length and Φ N and Φ M refer to the expanded unwrapped phase. Δ Φ l refers to the phase shift of the two regions on boundary l, whose quality can be described by the variance σ 2 . Considering the boundary length and variance, we define a discontinuity cost index Λ as follows:
Λ l = ρ σ Δ Φ 2 = ρ 2 l ( Φ N Φ M Δ Φ M N ) 2 d ρ .
Λ describes the continuity quality of expanding boundaries and is taken as the cost of discontinuity on the boundary. However, the mean of region phase difference Δ Φ l is unreliable at region boundaries with poor continuity.
Likewise, a region discontinuity index, the region residue, is defined similarly at the intersections of different expanded boundaries referred to as region corners (solid circles in Figure 2). The intersections of region boundaries and borders are called border corners (open circles in Figure 2). Assuming three boundaries or more meet at a certain point, the region residue r at the region corner is calculated as follows:
r = i = 1 n Δ Φ l i 2 π ,
where · is the rounding operator and the sum calculation is performed clockwise (or counterclockwise). Considering the case of raster data, the value of integer n is generally 3 or 4, and the residue values are generally ± 1 , ± 2 , , the larger the absolute value and the smaller the probability of occurrence according to our experiments.
For example, we calculated the residue value at region corner o 4 in the bottom right of Figure 2, and the phase difference values on the boundaries are noted in Figure 3. The calculation is as follows:
r o 4 = 1.6 + 0.8 + 0.2 = 2 1 + 0 = + 1 .
Non-zero residue values indicate that one or more phase discontinuity boundaries exist, which might cause regional inconsistencies. It is necessary to detect possible discontinuities before using phase differences to align the region phases.
We took the region corners and border corners as vertices and the region boundaries as connecting edges, which can form a graph structure (Figure 4a). A graph structure is used to illustrate the processing of region residues. Note that the closed boundaries without region corners, such as that between continuous regions C and D, are not shown in the graph, for no residue could be detected there.

2.2.2. Balancing Region Residues and Aligning Region Phases

The corners farther from residues are expected to have a higher continuity measure. Similarly, we define the continuity distance at region corner t as
R t ± = min Λ l i | l i L s t , s O ± Z
where O ± is the positive or negative residue set, Z is the border corner set, and L s t is a boundary subset connecting the reference residue s and the target corner t. The residues and border corners are taken as the source of the continuity distance, and Λ refers to the distance gain of the path sections. Taking the corners in Figure 4a as an example and assuming Λ values of each boundary and the positive residue o 2 ( r = + 1 ) , o 4 ( r = + 1 ) and negative residue o 1 ( r = 2 ) , we calculated the positive and negative continuity distances noted as ( R + , R ) in Figure 4a.
In order to balance the residues, it should find the path with the minimum continuity distance ( Λ ). The focus is on how to balance the dual residues, based on the distance map that is generated at all vertices. The ideal coupling edges W of the dual residues is expected to minimize the objective function as follows:
J = i = 1 n Λ l i , l i W ,
where l i visits all elements of W and n is the size of W . The minimization problem is N-P hard [32], and we provide an approximate solution scheme based on the corner continuity distance map.
Figure 4. A graph (a) derived from Figure 2 and its continuity distance map after the first coupling (b) and the second coupling (c). Distance values noted in the pair as ( R + , R ) near the vertices and Λ values noted beside the edges. Solid and open arrows indicate the positive and negative propagation, respectively.
Figure 4. A graph (a) derived from Figure 2 and its continuity distance map after the first coupling (b) and the second coupling (c). Distance values noted in the pair as ( R + , R ) near the vertices and Λ values noted beside the edges. Solid and open arrows indicate the positive and negative propagation, respectively.
Remotesensing 15 01370 g004
A pair of residues can be defined as a positive and a negative residue that are sources of the continuity distance for each other. The connecting edges between the residues are marked with phase discontinuities. With a positive residue point s A and a negative residue point s B , the dual pair Γ s B s A is expected to satisfy the following condition:
Γ s B s A : Υ s B + s A Υ s A s B , s A O + , s B O ,
where Υ s A is the source of the negative continuity distance of s A and Υ s B + is the source of the positive continuity distance of s B , while ↦ refers to the mapping operator. For residue pairs that need to take border corner s Z into account, the corresponding pairs Γ s Z s A and Γ s B s Z are defined as follows:
Γ s Z s A : Υ s A s Z , s A O + , s Z Z
Γ s B s Z : Υ s B + s Z , s B O , s Z Z .
According to the above definitions, it is convenient in Figure 4 to find region residue pairs Γ o 1 o 2 and Γ z 4 o 4 (Figure 4b). However, the residue balancing is not finished after the one time of pair detection, for residue o 1 ( r = 2 ) needs more operations.
Before the next pair detecting, we need to update the vertices’ continuity distance values affected by the balanced residues. Since the Λ value of each edge for coupling in Equation (9) is counted only once, Λ on the coupling edges was set as 0 when updating the distance values. A border corner is expected to contain both positive and negative residues indefinitely. The distance values of o 1 , o 2 , o 3 , and o 4 all were updated as shown in Figure 4b.
Because the newly detected connecting edge may overlap the previous detection, it is dealt with by strengthening in the same direction and weakening in the opposite direction. The merging result of the old and new pairs is marked in Figure 4c, showing that the edges o 2 o 3 and o 3 o 4 become coupled, and the edge o 4 z 4 becomes uncoupled, namely
Γ z 4 o 4 Γ o 1 z 4 = Γ o 1 o 4 ,
where ⊕ is the merging operator of the residue pairs. If all the residues are balanced after the detection of the pairs, the corresponding expanded boundaries form a discontinuity edge set W among the regions.
The operation of aligning unwrapped phases between continuous regions is to correct the integer-cycle shift. For example, the aligned phases of the region N are as follows:
Φ ¯ ( x , y ) = Φ ( x , y ) + 2 π K N , ( x , y ) N
where Φ refers to the unaligned unwrapped phase, Φ ¯ refers the phase after aligning, and K N is the integer-cycle offset. To align all the separate regions, one of them is set as the reference. For other regions, given the offset value of a neighboring region M and the connecting edge l, the offset of N is calculated as follows:
K N = K M + Δ Φ l 2 π , l W ,
where l cannot be in the discontinuity boundary set.

2.3. Region Growing Based on Heterogeneous Residual Diffusion

The proposed region growing method is based on heterogeneous residual diffusion outside of separate continuous regions, which ensures that the unwrapped phase outside regions is continuous and congruent with that inside.

2.3.1. Initializing Unwrapped Phase Outside Regions

To initialize the unwrapped phase outside of separate regions, we used the Poisson equation and the Dirichlet boundary condition. The Poisson equation is represented by
2 Φ ^ ( u , v ) = 2 Ψ ( u , v ) , ( u , v ) Ξ Φ ^ ( x , y ) = Φ ¯ , ( x , y ) Ξ ,
where 2 is the Laplacian operator, Φ ^ refers to the initial solution, 2 Ψ is the Laplacian of the wrapped phase Ψ , and Φ ¯ refers to the aligned unwrapped phases, which can guarantee a unique unwrapping phase solution, but it may not be congruent with the wrapped phases [33].
We then define a linear discontinuity indicator, which can be used to judge the existence of discontinuous boundaries in grid data. Given a horizontal or vertical line segment l, where unwrapped phases Φ ¯ s and Φ ¯ e are at the start and end, respectively, and wrapped phase samples, Ψ i , ( i = 1 , 2 , , n ) , are distributed in sequence, we define the indicator κ as
κ c = Φ ¯ s Φ ¯ e + i = 1 n + 1 W ( Ψ i Ψ i 1 ) ,
where W ( · ) is the wrapping operator, Ψ 0 = W ( Φ ¯ s ) , and Ψ n + 1 = W ( Φ ¯ e ) . When κ is non-zero, it indicates that l crosses one or more phase discontinuous boundaries, but note that opposite does not hold.
As shown in Figure 5, the linear discontinuity indicator can be used to identify areas where discontinuous boundaries exist, such as the region marked in green. Additionally, the indicator can be used to develop a decision scheme on solving for the unwrapped phase horizontally or vertically, by comparing the value of κ for horizontal and vertical line segments. Generally, solutions obtained using the Poisson equation are preferred when κ = 0 , compared to those of κ 0 or undefined.
In a more concise manner, the initial solution outside the continuous regions could be described as
Φ ^ 0 = Φ ^ x κ x = 0 , κ y 0 or undefined Φ ^ y κ y = 0 , κ x 0 or undefined Φ ^ x + Φ ^ y 2 otherwise .

2.3.2. Heterogeneous Residual Diffusion

The proposed method utilizes the concept of residual phase ϑ , defined as
ϑ = W ( Φ ^ Ψ ) ,
to guide the unwrapping process. Nonlinear diffusion of the residual phase is based on heterogeneous smoothing, which heavily smooths the residual phase in the direction where the wrapped phase is smoothest. This is accomplished by employing a differential equation and a boundary condition, represented by
ϑ t ( u , v , t ) = g ( | Ψ | 2 ) ϑ + h ( | Ψ | 2 ) ϑ
ϑ ( x , y ) = 0 , ( u , v ) Ξ , ( x , y ) Ξ ,
where ϑ refers to the second derivative of ϑ along the tangent direction of the contour of Ψ and ϑ refers to the second derivative in the gradient direction that is normal to the tangent direction (as shown in Figure 6).
ϑ = 1 | Ψ | 2 Ψ y Ψ x H ϑ Ψ y Ψ x T ϑ = 1 | Ψ | 2 Ψ x Ψ y H ϑ Ψ x Ψ y T
is used, where H ϑ is the Hessian matrix of ϑ , g ( · ) is a hypersurface minimal function defined as g ( s ) = ( 1 + s ) 1 2 , and h ( s ) = g ( s ) + 2 s g ( s ) .
In nonlinear diffusion, the iterative process is represented by
Φ ^ n + 1 = Φ ^ n + ( g ( | Ψ | 2 ) ϑ + h ( | Ψ | 2 ) ϑ ) Δ t
where Δ t is the time interval. Iterative nonlinear diffusion with known boundaries gradually adjusts the initial unwrapped phase to fit the wrapped phase.

2.3.3. Continuous Region Growing

After obtaining the iterative solution, we made the fitted solution congruent with its original wrapped phase in regions of high quality. Specifically, we select Φ ^ ( u , v ) if it satisfies the following condition:
( x u ) 2 + ( y v ) 2 δ 2 , | W ( Φ ^ ( u , v ) Ψ ( u , v ) ) | τ , ( x , y ) Ξ , ( u , v ) Ξ ,
where δ represents the sampling interval, | · | denotes the absolute operator, and τ is a threshold of deviation, with τ ( 0 , π ) .
The resulting congruent unwrapped phase Φ ˜ is represented by
Φ ˜ ( u , v ) = Φ ^ ( u , v ) + W ( Ψ ( u , v ) Φ ^ ( u , v ) ) ,
where ( u , v ) is considered a newly added element of the continuous regions. Through the iterative process, the continuous regions continue to expand and the outside ones shrink, until a globally consistent unwrapped phase is achieved.

3. Numerical Implementation

In the proposed method, the data were organized into three primary structures: Regular Grid, Regions, and Irregular Grid. These structures were designed to efficiently store and manipulate the data related to the phase samples, corners (region corners and border corners), and continuous regions.
The Regular Grid is the most-fundamental structure, which stores information about the wrapped phase grid, quality grid, distance grid, flag grid, source grid, discontinuity grid, unwrapped phase grid, and other elements (see Figure 7, left). The Regions structure stores information about continuous regions and includes elements such as the Regular Grid instance, region id array, region flag array, and region offset array (see Figure 7, middle). The Irregular Grid stores information related to corners and relies on both the Regular Grid and the Regions structures. It includes elements such as the Regions instance, corner id array, boundary id array, corner flag array, residue array, corner source array, corner distance array, boundary flag array, boundary cost array, boundary couple, and others (see Figure 7, right).
The proposed method was implemented by extracting continuous regions, unwrapping within these regions, aligning the regions, and growing the regions. To achieve computational efficiency, a GPU-based parallel processing scheme was designed for the main processing. In this parallel mode, massive local independent operations are carried out for each local grid point, corner, boundary, and surrounding neighborhood. The main process is briefly described as follows:
  • Construction of the phase continuity distance grid: the first step in constructing the phase continuity distance grid is the detection and marking of phase residues in the wrapped phase. This is followed by the calculation of the quality grid, which is defined as the inverse of the wrapped gradient variance, for example
    Q x = 1 σ x 2 , Q y = 1 σ y 2 ,
    where σ x 2 and σ y 2 refer to the variance of the phase gradient in the horizontal and vertical directions. The positive and negative continuity distance grids of the wrapped phase are then calculated according to Equation (2).
    As an optional operation, the procedure constructs a preliminary minimum balanced tree (MBT) for coupling dual-phase residues and marking simple discontinuity boundaries. This step can reduce the number of unbalanced residues and simplify the distribution of the following continuous regions. The parallel calculation of the distance grid and implementation of the dual residue pair coupling in the sampling phase can be found in our previous work [31]:
  • Unwrapping within the continuous regions: With the given unwrapped phase references, local direct unwrapping is performed independently in each region. The unwrapping process is based on the relations between Φ i , j and its unwrapped neighbors Φ u , v such that
    Φ i , j = Φ u , v + W ( Ψ i , j Ψ u , v ) , ( i , j ) Ξ , ( u , v ) N ( i , j ) Ξ ,
    where no known discontinuity boundary exists between ( i , j ) and ( u , v ) . This process is applied to all pixels in the regions.
    If there are any neighboring phase pixels Φ i and Φ i + 1 in the continuous regions, away from known discontinuity boundaries, that meet the condition:
    | Φ i , j Φ u , v | > π , ( i , j ) Ξ , ( u , v ) N ( i , j ) Ξ ,
    a higher threshold is required as in the previous step, and the unwrapped phases are re-evaluated again.
  • Continuous region boundary expansion: The local unwrapped phase in the continuous regions is expanded outwards as follows with known Φ u , v , if d i , j is greater than d u , v + 1 or undefined:
    Φ i , j = Φ u , v , d i , j = d u , v + 1 , ( i , j ) Ξ , ( u , v ) N ( i , j ) Ξ ,
    where d refers to the distance to the nearest grid region. The phase difference and discontinuity cost are calculated on each boundary between the regions according to Equations (5) and (6). Region residue detection is performed on each region corner according to Equation (7).
  • Corner continuity distance calculation: In order to improve the efficiency of corner continuity distance calculation, we utilized a parallel implementation of the fast march method (FMM) [34,35]. The positive and negative continuity distance of each corner point is determined through the examination of the discontinuity cost of boundaries and the distribution of region residues. The calculation is performed using known R u ± , in accordance with the following equation:
    s = arg min u { R u ± + Λ l | l : i u , u N ( i ) } R i ± = R s ± + Λ l s i , Υ i ± = Υ s ± ,
    where Λ l represents the discontinuity cost on boundary l and N ( i ) refers to the corner neighborhood of corner i.
  • Balancing region residues: We utilized the procedure defined in Equations (10)–(12) to identify pairs of detectable residues and assigned them the appropriate continuity distance values. Subsequently, the residues and boundary indexes are adjusted in accordance with the process outlined in Algorithm 1. Here, the boundary couple index, C, takes a value of zero in the case of an uncoupled pair, and the direction index of the pair, b, is defined as ± 1 , depending on whether the direction of the pair and the direction of the boundary are the same.
    Algorithm 1 Balancing region residues.
      1:
    repeat
      2:
    for all unbalanced residue ri in parallel do
      3:
      test coupling conditions in Equations (10)–(12)
      4:
      if found a pair Γ v u  then
      5:
        r u + = r u + 1 for Equations (10) and (11);
      6:
        r v = r v + 1 for Equations (10) and (12);
      7:
       trace pair boundaries L v u ;
      8:
       for  l j in L v u  do
      9:
         C j = C j + b ;
    10:
       end for
    11:
      end if
    12:
    end for
    13:
     update distance R ± and sources Υ ± ;
    14:
    until residues are all balanced.
  • Aligning region phase: A reference region was chosen, and the unwrapped phase of all other regions was adjusted accordingly, using Equation (13), on boundaries without a couple for C. This resulted in the establishment of a global phase continuity framework.
  • Initializing solution outside regions: The initial unwrapped phase outside the regions is obtained after detecting phase discontinuities in the horizontal and vertical directions. This was accomplished by solving the Poisson equations according to the linear discontinuities defined in Equation (14). The local Poisson equation in the horizontal or vertical direction outside the continuous regions can be written as
    1 2 1 Φ ^ i 1 Φ ^ i Φ ^ i + 1 T = 2 Ψ i ,
    and the equation changes at the boundaries. For example, if ( i + 1 ) is in a continuous region, the equation is rewritten as
    1 2 Φ ^ i 1 Φ ^ i T = 2 Ψ i + Φ ¯ i + 1 , ( i + 1 ) Ξ .
    In our GPU implementation, we employed the parallel cyclic reduction (PCR) method [36] to solve the Poisson equation.
  • Heterogeneous diffusion: Heterogeneous diffusion is applied to ensure that the unwrapped phase remains congruent with the wrapped phase, in accordance with Equations (17) and (18). This process is repeated to expand the regions until the phase stabilizes and converges. The final unwrapped result is obtained once convergence is achieved. It is important to note that the time interval Δ t should not break the stability of the diffusion process, and we set Δ t < 1 4 in our implementation.
The above-described main steps of the algorithms were implemented in a parallel manner utilizing the capabilities of GPU devices. A simplified representation of the entire process is presented in Figure 8.

4. Experiments and Analysis

We chose to validate the efficiency and reliability of the proposed PU method using three real interferograms acquired by TanDEM-X satellites. The choice to test the method on TanDEM-X data was motivated by the unique characteristics of the satellite’s data acquisition, which include the ability to collect high-resolution data from a variety of geographical regions and the ability to capture phase information with a high degree of accuracy.
The first interferogram covers the city of Wuhan, China, in 2011, and the observed area includes rivers and lakes, which result in separated regions with a high-quality local phase. Observations of Singapore, which cover the sea and several separate islands, comprised the second set of experimental data. The third experimental data came from observations of Songshan, China, and the observing area included undulating hills, which resulted in shadows and overlapping in the data. The experiments were carried out on a machine equipped with a CPU (8 cores, 2.70 GHz), a GPU (NVIDIA Quadro M2000M, 640 cores, 1.03 GHz) with 4 GB of device memory (1.25 GHz), a 16 GB host memory (dual-channel, DDR4, 2133 MHz), and a 64-bit operating system.
The three experimental interferograms were of 2944 × 5014 pixels, 1752 × 2984 pixels, and 2073 × 3529 pixels, respectively. There were 72,751 positive and 73,183 negative, 13,178 positive and 13,151 negative, and 86,080 positive and 86,093 negative phase residues detected in the three interferograms, respectively (Figure 9a–c). The initial distance maps were constructed based on these positive and negative residues (Figure 9d–f), and the coupling of residue pairs was performed four times for each of the three datasets. This process resulted in 10,094 pairs being found and 79 positive and 11 negative phase residues remaining uncoupled in the Wuhan data, 14,895 pairs and 193 positive and 153 negative residues remaining uncoupled in the Singapore data, and 82,719 pairs and 3718 positive and 3690 negative residues remaining uncoupled in the Songshan data.
The maximum values in the distance maps, D m a x , were 192,475, 3,013,342, and 128,245, respectively, and the initial ratio thresholds for region extraction were set at 0.1%, 0.1%, and 1.0%. The step interval was set at 0.1%, and the stop thresholds were 0.8%, 0.3%, and 1.9%, correspondingly. The extraction yielded 7504 continuous regions for the interferogram of Wuhan, 552 continuous regions for the interferogram of Singapore, and 7246 continuous regions for the interferogram of Songshan. The coverage percentage of extracted regions in the valid areas was 84.10, 87.64, and 53.87, respectively. These values are listed in Table 1.
A detailed analysis of the iterative extraction process of the Wuhan data is presented in Figure 10. As the threshold increased, it was observed that the corresponding detection area gradually decreased, and the extracted continuous region gradually became larger. It is noted that a higher threshold may decompose a large region into more small regions, resulting in an increase in the number of regions.
Using local phase references, direct unwrapping was performed in each region. Next, the unwrapped phase of the reliability regions was expanded to generate boundaries between regions. The numbers of boundaries with and without corners, as well as the number of corners were calculated for each of the three data, as shown in Table 2.
Region residue detection was then performed at the intersection of these boundaries, with the results being recorded for each of the three data. After iteratively constructing corner reliability maps and performing residue coupling (4, 6, and 6 times), the number of pairs of unbalanced residues was determined for each interferogram, as shown in Table 2. The number of unbalanced residues in each coupling iteration is presented in Figure 11.
In the detail of Figure 12, the phase in Regions A and B varied greatly, resulting in Positive Residues 1 and 2 and Negative Residues 3 and 4. After residue coupling, the marked discontinuity boundaries changed the phase cycle of Areas B, C, and D, which is a typical process of residue coupling and phase alignment. Region E, with an extreme phase change, had different effects on the upper and lower regions.
After aligning the regions for the global frame, Poisson solving with the Dirichlet boundary condition provided the initial unwrapped phase, considering the discontinuities defined in Equation (14). The initial solutions were obtained from four times of PCR iterative solving with wrapped phase differences along the horizontal and vertical directions. The three initial solutions were then subjected to 1000 iterations of nonlinear diffusion, respectively. Congruent areas were updated every 20 iterations with τ = 0.05 for the data of Wuhan, τ = 0.12 for the data of Singapore, and τ = 0.05 for the data of Songshan. The initial unwrapped phase from PCR solving and the final results from nonlinear diffusing are shown in Figure 13.
In order to provide further insight into the processing, two windowed areas of 250 × 250 pixels and 350 × 350 pixels were analyzed in detail. Initial solutions filled the areas outside the continuous regions, and the continuity of the phase between outside and inside the continuous regions was preserved. However, it is noted that stripes appeared in areas with severe noise, which was caused by the integration path passing through phase discontinuity boundaries. The final results from nonlinear diffusion made effective corrections in the areas with stripes. Some small pits appeared in the areas with severe noise, which was attributed to the attenuation of the smoothing effect of nonlinear diffusion in these noise-affected areas, where the effect of guidance from the wrapped gradient was not reliable.
To further evaluate the effectiveness of the proposed method, we calculated the wrapped phase differences between the unwrapped and original wrapped phases. This quantity provides insight into the impact of the unwrapping process, specifically the removal of cycle jumps. As shown in Figure 14, the residual phase maps revealed that the proposed method effectively completed the unwrapping process, as evidenced by the minimal and scattered residual phases after masking out water bodies in the Wuhan and Singapore data. The roots of variance of the residuals in the results were 0.0343, 0.0470, and 0.0928, respectively. The area coverage of residual absolute values larger than 0.012 were 1.951%, 2.545%, and 9.672%. However, it is noted that the relatively high number of residual phases in the Songshan region was a result of strip-shaped shadows or layovers caused by steep mountainous terrain. Despite the use of a smoothing procedure using nonlinear diffusion, the traces of the discontinuity boundaries could not be fully eliminated.
We compared the performance of our proposed method with several existing unwrapping methods, including SNAPHU [10], region growing [6], GAMMA minimum cost flow (MCF), and PHASS [37,38] by processing the same dataset (as illustrated in Figure 15). The SNAPHU method and GAMMA method employ similar optimization schemes based on network flow, while the region growing method utilizes a growing scheme that relies on reliability sorting. For SNAPHU, we utilized the DEFO mode with coherence as the weight and used the final solution. For the region growing method, we used the pixel-level reliability to connect different growing regions. For the GAMMA method, we used the minimum cost flow algorithm to minimize the cycle errors. For the PHASS method, graph optimization algorithms based on MCF with regions were employed.
The results of the unwrapping process indicated that SNAPHU performed well in the Wuhan and Singapore cases within land areas. However, SNAPHU’s connection processing between regions was poor, resulting in cycle jumps between regions. The region growing method, on the other hand, focuses on connecting different growing regions at pixel-level continuity, but the instability of the noisy areas resulted in different cycle deviation patches. The GAMMA method showed similar results with cycle errors among separated regions. Despite the small area of error, the PHASS method gave better results on the Wuhan and Singapore data, but it performed poorly on the Songshan data. For the Songshan case, SNAPHU, region growing, GAMMA, and PHASS all failed to produce a global continuity solution across noisy areas.
We conducted an analysis of the runtime costs in our PU experiments (as shown in Table 3). On the data experiments of Wuhan and Singapore, the PHASS method spent the least time, mainly due to its low coherence region masking and CPU parallel acceleration. The times of our method were close to the PHASS method on the two data. On the Songshan data, our method significantly took the shortest time.
A further examination of the time costs of each processing step in our proposed method revealed that, in the Wuhan data experiment, 70.695 s were spent, with 15.8% of that time allocated to constructing the distance map, 15.4% to extracting the continuous region, 3.2% to aligning regions, and 65.6% to solving outside of regions (Table 4). In the Singapore data experiment, 3.809 s, 0.838 s, 0.555 s, and 18.122 s were spent on these steps, respectively. Similarly, in the Songshan experiment, 7.858 s, 6.456 s, 1.453 s, and 20.230 s were allocated to these steps, respectively. We observed that more than half of the total runtime was spent solving outside regions in all three experiments.

5. Discussion

While the proposed unwrapping method was able to achieve globally consistent results in our experiments, there are still some aspects that need further discussion. One of these is the use of phase residue MBT coupling during the construction of the continuity distance map as a preprocessing step. This step effectively reduces the number of phase residues and generates a distance map with larger continuous regions. However, it is imperative to note that more coupling iterations may not always produce better results. In some cases, excessive coupling may lead to incorrect results, which highlights the importance of careful consideration when determining the optimal number of iterations.
Another issue that needs to be addressed is the relatively high processing time required for unwrapping outside regions. Improving the efficiency of this aspect of the method should be a focus of future research. This may involve modifying or replacing the proposed nonlinear diffusion technique, particularly if strict adherence to discontinuity boundaries is necessary. Additionally, the setting of the congruity parameter τ plays a significant role in both the speed of region growing and the location of discontinuity boundaries. A larger value of τ relaxes the congruity residuals, which improves efficiency, but also increases the likelihood of the convergent regions expanding across the true discontinuity boundary, resulting in an incorrect boundary location. On the other hand, a smaller value of τ results in more accurate discontinuity boundaries, but may prevent the expansion of regions, requiring a longer convergence time or even preventing convergence altogether. Through our experiments, we found that the value of τ should be around 0.1 cycles according to our empirical analysis.

6. Conclusions

In this paper, a novel unwrapping method based on separate continuous regions was proposed to address the PU problem of regions separated by noise while neighboring continuity cannot be guaranteed. The proposed method was based on a continuity distance measure, which is computed as the minimum weighted distance from the dual-phase residues. Distance thresholds are used to identify regions with good phase continuity, which enables correct unwrapping in local regions. By assuming phase continuity between regions, the unwrapped phase of the region is expanded to calculate cycle numbers and continuity quality between regions. The region residues are detected using shift numbers. The coupling of region residues is designed to identify region discontinuities, which are considered during constructing the global continuity framework. The proposed method also employs Poisson solving and nonlinear diffusion to solve outside the unwrapped regions.
The proposed method effectively addresses the issue of continuous region isolation caused by noise, by using the coupling of region residues to ensure the reliability of the phase alignment of separated regions. Additionally, the method is easily implemented in parallel, which can significantly improve the efficiency of unwrapping processing. The results of the experiments conducted on real InSAR data demonstrated the reliability and efficiency of the developed method. This method was able to overcome the limitations of prior unwrapping methods, which were seen as less reliable when dealing with regions separated by noise. The proposed method is a step forward in the field of InSAR data unwrapping. It is expected to have a positive impact on a wide range of applications that depend on InSAR data.

Author Contributions

Conceptualization, H.J. and J.G.; methodology, J.G. and H.J.; software, J.G. and H.J.; validation, H.J., Z.S. and J.G.; formal analysis, J.G. and H.J.; investigation, H.J.; resources, H.J., Z.S. and R.W.; data curation, J.G. and H.J.; writing—original draft preparation, J.G. and H.J.; writing—review and editing, J.G. and Y.H.; visualization, J.G. and H.J.; supervision, H.J. and R.W.; project administration, H.J. and J.G.; funding acquisition, Z.S., H.J. and J.G. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the Key R&D Program Projects in Hainan Province (Nos. ZDYF2020192 and ZDYF2019008), the National Natural Science Foundation of China (No. 42274028), the Doctoral Research Foundation of Anhui Jianzhu University (No. 2021QDZ07), and the CSC Scholarship (No. 202008320034).

Data Availability Statement

The data underlying this study are available on request.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Ghiglia, D.C.; Pritt, M.D. Two-Dimensional Phase Unwrapping: Theory, Algorithms, and Software; Wiley: Hoboken, NJ, USA, 1998. [Google Scholar]
  2. Yu, H.; Lan, Y.; Yuan, Z.; Xu, J.; Lee, H. Phase Unwrapping in InSAR: A Review. IEEE Geosci. Remote Sens. Mag. 2019, 7, 40–58. [Google Scholar] [CrossRef]
  3. Xu, W.; Cumming, I. A region-growing algorithm for InSAR phase unwrapping. IEEE Trans. Geosci. Remote Sens. 1999, 37, 124–134. [Google Scholar] [CrossRef] [Green Version]
  4. Karsa, A.; Shmueli, K. SEGUE: A Speedy rEgion-Growing Algorithm for Unwrapping Estimated Phase. IEEE Trans. Med Imaging 2019, 38, 1347–1357. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  5. Ojha, C.; Manunta, M.; Pepe, A.; Paglia, L.; Lanari, R. An innovative region growing algorithm based on Minimum Cost Flow approach for Phase Unwrapping of full-resolution differential interferograms. In Proceedings of the 2012 IEEE International Geoscience and Remote Sensing Symposium, Munich, Germany, 22–27 July 2012; pp. 5582–5585. [Google Scholar] [CrossRef]
  6. Herráez, M.A.; Burton, D.R.; Lalor, M.J.; Gdeisat, M.A. Fast two-dimensional phase-unwrapping algorithm based on sorting by reliability following a noncontinuous path. Appl. Opt. 2002, 41, 7437–7444. [Google Scholar] [CrossRef] [PubMed]
  7. Xie, X.; Zeng, Q. Multi-baseline extended particle filtering phase unwrapping algorithm based on amended matrix pencil model and quantized path-following strategy. J. Syst. Eng. Electron. 2019, 30, 78–84. [Google Scholar] [CrossRef]
  8. Dudczyk, J.; Kawalec, A. Optimizing the Minimum Cost Flow Algorithm for the Phase Unwrapping Process in SAR Radar. Bull. Pol. Acad. Sci. Tech. Sci. 2014, 62, 511–516. [Google Scholar] [CrossRef] [Green Version]
  9. Huang, Q.; Zhou, H.; Dong, S.; Xu, S. Parallel Branch-Cut Algorithm Based on Simulated Annealing for Large-Scale Phase Unwrapping. IEEE Trans. Geosci. Remote Sens. 2015, 53, 3833–3846. [Google Scholar] [CrossRef]
  10. Chen, C.; Zebker, H. Phase unwrapping for large SAR interferograms: Statistical segmentation and generalized network models. IEEE Trans. Geosci. Remote Sens. 2002, 40, 1709–1719. [Google Scholar] [CrossRef] [Green Version]
  11. Sawaf, F.; Tatam, R.P. Finding minimum spanning trees more efficiently for tile-based phase unwrapping. Meas. Sci. Technol. 2006, 17, 1428. [Google Scholar] [CrossRef]
  12. Dong, J.; Zhuo, Z.; Li, J.; He, Y. A new phase image reconstruction method using Markov random fields. In Proceedings of the 2017 IEEE/ACIS 16th International Conference on Computer and Information Science (ICIS), Wuhan, China, 24–26 May 2017; pp. 1–5. [Google Scholar] [CrossRef]
  13. Xu, J.; An, D.; Huang, X.; Wang, G. Phase Unwrapping for Large-Scale P-Band UWB SAR Interferometry. IEEE Geosci. Remote Sens. Lett. 2015, 12, 2120–2124. [Google Scholar] [CrossRef]
  14. Zhang, Y.; Xing, Z. A Hybrid Phase Unwrapping Algorithm Based on Quality-Guided and Surface-Fitting. In Proceedings of the 2018 IEEE International Workshop on Electromagnetics: Applications and Student Innovation Competition (iWEM), Nagoya, Japan, 29–31 August 2018; pp. 1–2. [Google Scholar] [CrossRef]
  15. Zhang, Y.; Xing, Z. A Phase Unwrapping Algorithm Based on Branch-cut and B-Spline Fitting in InSAR. In Proceedings of the 2018 IEEE International Conference on Signal Processing, Communications and Computing (ICSPCC), Hong Kong, China, 14–16 September 2018; pp. 1–3. [Google Scholar] [CrossRef]
  16. Yu, H.; Lee, H. A convex hull algorithm based fast large-scale two-dimensional phase unwrapping method. In Proceedings of the 2017 IEEE International Geoscience and Remote Sensing Symposium (IGARSS), Fort Worth, TX, USA, 23–28 July 2017; pp. 3824–3827. [Google Scholar] [CrossRef]
  17. Yu, H.; Lan, Y.; Xu, J.; An, D.; Lee, H. Large-Scale L0-Norm and L1-Norm 2-D Phase Unwrapping. IEEE Trans. Geosci. Remote Sens. 2017, 55, 4712–4728. [Google Scholar] [CrossRef]
  18. Zhang, Y.; Zhang, S.; Gao, Y.; Li, S.; Jia, Y.; Li, M. Adaptive Square-Root Unscented Kalman Filter Phase Unwrapping with Modified Phase Gradient Estimation. Remote Sens. 2022, 14, 1229. [Google Scholar] [CrossRef]
  19. Gao, Y.; Tang, X.; Li, T.; Lu, J.; Li, S.; Chen, Q.; Zhang, X. A Phase Slicing 2-D Phase Unwrapping Method Using the L 1 -Norm. IEEE Geosci. Remote Sens. Lett. 2022, 19, 1–5. [Google Scholar] [CrossRef]
  20. Li, H.; Zhong, H.; Ning, M.; Zhang, P.; Tang, J. Using neural networks to create a reliable phase quality map for phase unwrapping. Appl. Opt. 2023, 62, 1206. [Google Scholar] [CrossRef] [PubMed]
  21. Zhou, L.; Yu, H.; Lan, Y. Deep Convolutional Neural Network-Based Robust Phase Gradient Estimation for Two-Dimensional Phase Unwrapping Using SAR Interferograms. IEEE Trans. Geosci. Remote Sens. 2020, 58, 4653–4665. [Google Scholar] [CrossRef]
  22. Li, L.; Zhang, H.; Tang, Y.; Wang, C.; Gu, F. InSAR Phase Unwrapping by Deep Learning Based on Gradient Information Fusion. IEEE Geosci. Remote Sens. Lett. 2022, 19, 1–5. [Google Scholar] [CrossRef]
  23. Zhao, J.; Liu, L.; Wang, T.; Wang, X.; Du, X.; Hao, R.; Liu, J.; Liu, Y.; Zhang, J. VDE-Net: A two-stage deep learning method for phase unwrapping. Opt. Express 2022, 30, 39794. [Google Scholar] [CrossRef] [PubMed]
  24. Pu, L.; Zhang, X.; Zhou, Z.; Li, L.; Zhou, L.; Shi, J.; Wei, S. A Robust InSAR Phase Unwrapping Method via Phase Gradient Estimation Network. Remote Sens. 2021, 13, 4564. [Google Scholar] [CrossRef]
  25. Sawaf, F.; Groves, R.M. Phase discontinuity predictions using a machine-learning trained kernel. Appl. Opt. 2014, 53, 5439–5447. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  26. Wu, Z.; Wang, T.; Wang, Y.; Wang, R.; Ge, D. Deep-Learning-Based Phase Discontinuity Prediction for 2-D Phase Unwrapping of SAR Interferograms. IEEE Trans. Geosci. Remote Sens. 2022, 60, 1–16. [Google Scholar] [CrossRef]
  27. Spoorthi, G.E.; Gorthi, S.; Gorthi, R.K.S.S. PhaseNet: A Deep Convolutional Neural Network for Two-Dimensional Phase Unwrapping. IEEE Signal Process. Lett. 2019, 26, 54–58. [Google Scholar] [CrossRef]
  28. Sica, F.; Calvanese, F.; Scarpa, G.; Rizzoli, P. A CNN-Based Coherence-Driven Approach for InSAR Phase Unwrapping. IEEE Geosci. Remote Sens. Lett. 2022, 19, 1–5. [Google Scholar] [CrossRef]
  29. Zhou, L.; Yu, H.; Lan, Y.; Xing, M. Deep Learning-Based Branch-Cut Method for InSAR Two-Dimensional Phase Unwrapping. IEEE Trans. Geosci. Remote Sens. 2022, 60, 1–15. [Google Scholar] [CrossRef]
  30. Li, H.; Zhong, H.; Zhang, P.; Ning, M.; Tang, J. Two-Dimensional Phase Unwrapping Based on Residual Prediction Neural Network. IEEE Geosci. Remote Sens. Lett. 2023, 20, 1–5. [Google Scholar] [CrossRef]
  31. Gao, J.; Sun, Z. Phase unwrapping method based on parallel local minimum reliability dual expanding for large-scale data. J. Appl. Remote Sens. 2019, 13, 038506. [Google Scholar] [CrossRef] [Green Version]
  32. Chen, C.W.; Zebker, H.A. Two-dimensional phase unwrapping with use of statistical models for cost functions in nonlinear optimization. JOSA A 2001, 18, 338–351. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  33. Lyuboshenko, I.V.; Maitre, H.; Maruani, A. Least-mean-squares phase unwrapping by use of an incomplete set of residue branch cuts. Appl. Opt. 2002, 41, 2129–2148. [Google Scholar] [CrossRef]
  34. Capozzoli, A.; Curcio, C.; Liseno, A.; Savarese, S. A comparison of Fast Marching, Fast Sweeping and Fast Iterative Methods for the solution of the eikonal equation. In Proceedings of the 2013 21st Telecommunications Forum Telfor (TELFOR), Belgrade, Serbia, 26–28 November 2013; pp. 685–688. [Google Scholar]
  35. Jeong, W.K.; Whitaker, R.T. A Fast Iterative Method for Eikonal Equations. SIAM J. Sci. Comput. 2008, 30, 2512–2534. [Google Scholar] [CrossRef] [Green Version]
  36. Swarztrauber, P.N. The Methods of Cyclic Reduction, Fourier Analysis and the FACR Algorithm for the Discrete Solution of Poisson’s Equation on a Rectangle. SIAM Rev. 1977, 19, 490–501. [Google Scholar] [CrossRef]
  37. Tophu. Available online: https://github.com/opera-adt/tophu (accessed on 16 January 2023).
  38. isce-framework/isce3: InSAR Scientific Computing Environment. Available online: https://github.com/isce-framework/isce3 (accessed on 16 January 2023).
Figure 1. The definition of the positive continuity distance ( D + ) and the negative continuity distance ( D ) at point t. The red lines connecting dual residues represent the discontinuity boundaries, and the black lines refer to the the shortest paths to the residues. The full (sum) continuity distance (D) in the gray areas is higher than a given threshold.
Figure 1. The definition of the positive continuity distance ( D + ) and the negative continuity distance ( D ) at point t. The red lines connecting dual residues represent the discontinuity boundaries, and the black lines refer to the the shortest paths to the residues. The full (sum) continuity distance (D) in the gray areas is higher than a given threshold.
Remotesensing 15 01370 g001
Figure 2. An example of continuous regions A, B, C, D, E, F, and their expanding boundaries (dashed lines), region corner o i ( i = 1 , 2 , 3 , 4 ) (solid circles), and border corner z i ( i = 1 , 2 , 3 , 4 ) (open circles).
Figure 2. An example of continuous regions A, B, C, D, E, F, and their expanding boundaries (dashed lines), region corner o i ( i = 1 , 2 , 3 , 4 ) (solid circles), and border corner z i ( i = 1 , 2 , 3 , 4 ) (open circles).
Remotesensing 15 01370 g002
Figure 3. A region corner among continuous regions and region phase differences calculated clockwise.
Figure 3. A region corner among continuous regions and region phase differences calculated clockwise.
Remotesensing 15 01370 g003
Figure 5. The linear discontinuity indicator in different cases. Horizontal κ x = 0 on point a, vertical κ y = 0 on point b, κ y 0 on point c, and κ x and κ y both undefined on point d. The lines in red mark two discontinuity boundaries, and the area in green defines a dangerous region.
Figure 5. The linear discontinuity indicator in different cases. Horizontal κ x = 0 on point a, vertical κ y = 0 on point b, κ y 0 on point c, and κ x and κ y both undefined on point d. The lines in red mark two discontinuity boundaries, and the area in green defines a dangerous region.
Remotesensing 15 01370 g005
Figure 6. Nonlinear residual diffusion and region growing. The diffusion in the tangent direction of the contour is designed to be stronger than that in the normal direction, which drives regions to grow faster in flatter directions.
Figure 6. Nonlinear residual diffusion and region growing. The diffusion in the tangent direction of the contour is designed to be stronger than that in the normal direction, which drives regions to grow faster in flatter directions.
Remotesensing 15 01370 g006
Figure 7. The main data and processing of structures Regular Grid, Regions, and Irregular Grid in the implementation.
Figure 7. The main data and processing of structures Regular Grid, Regions, and Irregular Grid in the implementation.
Remotesensing 15 01370 g007
Figure 8. A simple demonstration of the processing flow for the proposed PU method based on separated continuous regions.
Figure 8. A simple demonstration of the processing flow for the proposed PU method based on separated continuous regions.
Remotesensing 15 01370 g008
Figure 9. The wrapped phases (number of 2 π ) from Wuhan (a), Singapore (b), and Songshan (c) and their phase distance maps (df) after four times of residue coupling, respectively. The water areas in (d,e) are masked.
Figure 9. The wrapped phases (number of 2 π ) from Wuhan (a), Singapore (b), and Songshan (c) and their phase distance maps (df) after four times of residue coupling, respectively. The water areas in (d,e) are masked.
Remotesensing 15 01370 g009
Figure 10. Area coverage and number of changes of regions of the Wuhan data in the iterative extracting process.
Figure 10. Area coverage and number of changes of regions of the Wuhan data in the iterative extracting process.
Remotesensing 15 01370 g010
Figure 11. Numbers of unbalanced residues in the coupling iterations in the experimental data.
Figure 11. Numbers of unbalanced residues in the coupling iterations in the experimental data.
Remotesensing 15 01370 g011
Figure 12. The local unwrapped phase in the continuous regions for Wuhan (a), Singapore (b), and Songshan (c) and their global aligned versions in (df), respectively. Details of windowed areas with boundaries (black), residues (red for positive, green for negative), and coupling edges (red).
Figure 12. The local unwrapped phase in the continuous regions for Wuhan (a), Singapore (b), and Songshan (c) and their global aligned versions in (df), respectively. Details of windowed areas with boundaries (black), residues (red for positive, green for negative), and coupling edges (red).
Remotesensing 15 01370 g012
Figure 13. The initial unwrapped phase outside the continuous regions from PCR for Wuhan (a), Singapore (b), and Songshan (c) and their final version (df) with 1000 times of nonlinear diffusion, respectively. The detailed local unwrapped phase in the regions shows the unwrapped version and diffusion version of the areas marked.
Figure 13. The initial unwrapped phase outside the continuous regions from PCR for Wuhan (a), Singapore (b), and Songshan (c) and their final version (df) with 1000 times of nonlinear diffusion, respectively. The detailed local unwrapped phase in the regions shows the unwrapped version and diffusion version of the areas marked.
Remotesensing 15 01370 g013
Figure 14. The wrapped phase difference between our diffused results and the wrapped phase for the data of Wuhan (a), Singapore (b), and Songshan (c).
Figure 14. The wrapped phase difference between our diffused results and the wrapped phase for the data of Wuhan (a), Singapore (b), and Songshan (c).
Remotesensing 15 01370 g014
Figure 15. The wrapped phase for the experimental data of Wuhan (ad), Singapore (eh), and Songshan (il) from the SNAPHU, region growing, GAMMA, and PHASS PU methods.
Figure 15. The wrapped phase for the experimental data of Wuhan (ad), Singapore (eh), and Songshan (il) from the SNAPHU, region growing, GAMMA, and PHASS PU methods.
Remotesensing 15 01370 g015
Table 1. Statistics of continuous region extraction on the experimental data.
Table 1. Statistics of continuous region extraction on the experimental data.
Data D max γ 0  (%) δ γ  (%) γ 1  (%)RegionsCoverage (%)
Wuhan192,4750.10.10.8750484.10
Singapore3,013,3420.10.10.355287.64
Songshan128,2451.00.11.9724653.87
Table 2. Number statistics of boundaries, corners, residues, and pairs in the experimental data.
Table 2. Number statistics of boundaries, corners, residues, and pairs in the experimental data.
DataBound a Bound b Corner r + 1 r 1 r + 2 r 2 r other Pair
Wuhan19,25352612,890188182250192
Singapore1402359739695210105
Songshan20,29820013,762100499791101054
Bounda is for the boundaries with corners, Boundb for the boundaries without corners, and rother for the residues greater than 2 or less than −2.
Table 3. Runtime (sec) cost of the five different PU methods.
Table 3. Runtime (sec) cost of the five different PU methods.
DataSNAPHURegion GrowingGAMMMAPHASSProposed
Wuhan98.07285.91973.51067.04970.695
Singapore86.99227.80964.58023.04623.324
Songshan93.01456.511161.61043.06835.997
Table 4. Runtime (%) of four steps in the proposed method.
Table 4. Runtime (%) of four steps in the proposed method.
DataConstructing MapExtracting RegionsAligning RegionsSolving Outside
Wuhan15.8%15.4%3.2%65.6%
Singapore16.3%3.6%2.4%77.7%
Songshan21.8%17.9%4.1%56.2%
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Gao, J.; Jiang, H.; Sun, Z.; Wang, R.; Han, Y. A Parallel InSAR Phase Unwrapping Method Based on Separated Continuous Regions. Remote Sens. 2023, 15, 1370. https://doi.org/10.3390/rs15051370

AMA Style

Gao J, Jiang H, Sun Z, Wang R, Han Y. A Parallel InSAR Phase Unwrapping Method Based on Separated Continuous Regions. Remote Sensing. 2023; 15(5):1370. https://doi.org/10.3390/rs15051370

Chicago/Turabian Style

Gao, Jian, Houjun Jiang, Zhongchang Sun, Ruisheng Wang, and Youmei Han. 2023. "A Parallel InSAR Phase Unwrapping Method Based on Separated Continuous Regions" Remote Sensing 15, no. 5: 1370. https://doi.org/10.3390/rs15051370

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop