16 January 2018 Interferometric synthetic aperture radar phase unwrapping based on sparse Markov random fields by graph cuts
Author Affiliations +
Abstract
Phase unwrapping (PU) is one of the key processes in reconstructing the digital elevation model of a scene from its interferometric synthetic aperture radar (InSAR) data. It is known that two-dimensional (2-D) PU problems can be formulated as maximum a posteriori estimation of Markov random fields (MRFs). However, considering that the traditional MRF algorithm is usually defined on a rectangular grid, it fails easily if large parts of the wrapped data are dominated by noise caused by large low-coherence area or rapid-topography variation. A PU solution based on sparse MRF is presented to extend the traditional MRF algorithm to deal with sparse data, which allows the unwrapping of InSAR data dominated by high phase noise. To speed up the graph cuts algorithm for sparse MRF, we designed dual elementary graphs and merged them to obtain the Delaunay triangle graph, which is used to minimize the energy function efficiently. The experiments on simulated and real data, compared with other existing algorithms, both confirm the effectiveness of the proposed MRF approach, which suffers less from decorrelation effects caused by large low-coherence area or rapid-topography variation.
Zhou, Chai, Xia, Ma, and Lin: Interferometric synthetic aperture radar phase unwrapping based on sparse Markov random fields by graph cuts

1.

Introduction

Interferometric synthetic aperture radar (InSAR) is a powerful tool to measure the digital elevation model (DEM) of the earth’s surface.1 The InSAR system makes use of two or more SAR images covering the same scene to generate interferograms. The acquired interferometric phase is the principal value of the absolute phase, i.e., modulo-2π observation, which corresponds to topography or deformation and is called the wrapped phase. Under this condition, phase unwrapping (PU) is the process of resolving the absolute phase through the wrapped phase. Unfortunately, from a mathematical viewpoint, PU is an ill-posed problem, whose significant property is that there is no unique solution to it, if no further information is added. In fact, an assumption taken by most PU algorithms is that the absolute value of phase differences between neighboring pixels is less than π, the so-called Itoh condition.2 In essence, this assumption requires that the measured regions have spatial continuity. The Itoh condition can be violated if the absolute phase surface is discontinuous, or if the wrapped phase is dominated by noise. The existence of large low-coherence regions and rapid-topography variations poses challenges to two-dimensional (2-D) PU algorithms.

PU approaches comprise three main classes: path following, minimum norm, and Bayesian approaches. Path following algorithms apply line integration schemes over the wrapped phase image, and basically rely on the assumption that the Itoh condition holds along the integration path. This approach includes the so-called branch cuts,3,4 quality-guided,5 and minimum discontinuity,6 etc. These path following algorithms, although 2π-true, are sensitive to the Itoh conditions, which are often violated in the 2-D case. Minimum norm methods try to find a phase solution for which the Lp norm of the difference between absolute phase differences and wrapped phase differences is minimized. When p=2, the Lp norm-based PU method is transformed to the least squares (LS) PU method,7 and when p=1, it is changed to the minimum-cost flow (MCF) PU method.8 A drawback of the LS algorithm is that it tends to smooth discontinuities, unless they are provided as binary weights.4 The rather popular MCF algorithm converts a 2-D PU problem to an MCF problem on a network, which connects the residues with branch cuts and minimizes the total weighted length. The L1 norm-based criterion performs better than the L2 norm in preserving discontinuities.4 However, the MCF approach, as a 2π-true unwrapper, introduces integer times of 2π error on discrete points.

The Bayesian approach relies on a data observation mechanism model, as well as a priori knowledge of the phase to be modeled.910.11 Based on this approach, Markov random fields (MRFs) constraining global phase variations are used for solving a 2-D PU problem, which has been proposed in several studies.1213.14 Ying et al.12 proposed an MRF approach for a 2-D PU problem. This approach utilizes an efficient algorithm for parameter estimation using a series of dynamic programming connected by the iterated conditional modes. Dias and Valadao13 proposed a PU algorithm based on graph cuts. This algorithm solves integer optimization problems by computing a sequence of binary optimizations, which is converted to the {0, 1} cut of graphs, so the graph cuts approach is utilized. Chen et al.14 proposed an integrated phase denoising and unwrapping algorithm based upon MRF, the energy function of which is defined according to the local independence property inferred from the MRF structure and then minimized to obtain the estimate of the true phase value. Experiments on synthetic and real data prove that this algorithm performs better than other competing algorithms.

Compared to some commonly used PU algorithms, the MRF approach for PU has two major advantages. First, strength of the MRF approach lies in its ability to take contextual information into account, making full use of local properties of interferograms. The pairwise difference of the unwrapped phase over all neighbors is minimized as a global optimization. Second, some newly developed optimization algorithms, e.g., graph cuts13 and loopy belief propagation,15 etc., have proven to be very powerful to solving this minimization for a 2-D PU problem. However, since the traditional MRF algorithm for PU is usually defined on a rectangular grid, such as the well-known four-connected grid, it fails easily if large parts of the wrapped data are influenced by high phase noise caused by a large low-coherence area or terrain with large topographic gradients, which normally makes unwrapping a hard task. Therefore, the case in which data are available on a sparse set of points is considered. Accordingly, the MCF approach has been subsequently changed into the Delaunay-MCF approach in order to deal with sparse data.1617.18 What is more, multitemporal/multibaseline InSAR (i.e., persistent scatterers interferometry, PSInSAR) attracted increasing interests,1920.21 it also deals with only a sparse set of pixels, so in that case, sparse-grid PU algorithms are indispensable.

In this paper, we propose a sparse MRF method for PU, which extends the traditional regular-grid MRF algorithm to deal with a sparse data grid. Our proposed algorithm differs from the regular-grid MRF algorithm mainly in two aspects. First, the proposed algorithm uses only the subset of the data that are considered as coherent pixels, which are based on the coherence coefficient of the SAR images. Since the selected data are no longer on a regular grid, the Delaunay triangulation method is then used to divide the entire SAR image into a series of 2-D triangles according to the coherent pixels. Second, to speed up the graph cuts algorithm for sparse MRF, dual elementary graphs are designed and merged to obtain the Delaunay triangle graph, which is used to minimize energy function efficiently. The remainder of this paper is organized as follows. In Sec. 2, we explain the formulation of the sparse MRF based on Delaunay triangulation and describe the implementations of the graph cuts algorithm to optimize. In Sec. 3, we give the two simulated experiments, in which some deeper analysis and comparison between the proposed MRF algorithm and other existing algorithms is provided. In Sec. 4, two experiments using real data are carried out, and the performance of the proposed MRF method is also compared with other existing algorithms. Finally, the study’s conclusions are presented in Sec. 5.

2.

Sparse Markov Random Field Method for Phase Unwrapping

In this section, we describe sparse MRF representation based on Delaunay triangulation for a 2-D PU problem and implementations of the graph cuts optimization algorithm. Our method is inspired from the fact that areas of low interferometric coherence are affected by high phase noise and therefore are excluded during the PU process. Given this situation, a sparse set of coherent points is generated using a coherence coefficient of the SAR images. The image coordinates of the coherent points are then used to create a 2-D Delaunay triangulation. A prior of local phase difference is computed to unwrap the phase through minimization of the local phase difference. In particular, this prior is formed by computing a pairwise energy function induced only by the coherent point and the triangulated grid. Moreover, the minimization of the energy function is achieved via the graph cuts algorithm, in order to escape local minima.

2.1.

Coherent Point Generation and Delaunay Triangulation

We define coherent points to be pixels that can be robustly unwrapped due to their high coherence coefficient of the SAR images. We produce the interferometric coherence coefficient map estimated from the normalized interferogram and the coregistered intensity images, which is defined in Ref. 22. Accordingly, we generate coherent points by computing and applying a threshold to the coherence coefficient map. The main difficulty with defining coherent points is the selection of a threshold. If the threshold is too low, too many low-quality phase values will be able to corrupt the unwrapped solution. On the other hand, if the threshold is too high, then too many pixels will be masked out, and phase values on adjacent pixels will violate the Itoh condition. Hence, for selecting the threshold, we use an automated method proposed in Ref. 4 by examining the distribution of the coherence coefficient map, which helps to reach a balance between quality and quantity of coherent points. Since low-coherent points are masked out, the running time of the subsequent optimization algorithm is reduced. The 2-D triangulation of the coherent points aims to represent the entire image with a set of 2-D triangles, which reflects the pairwise difference of the unwrapped phase between a pixel and its neighboring pixels. Many 2-D triangulation methods exist, and the representative method is Delaunay triangulation, which can be generated by a divide and conquer algorithm for triangulations in 2-D case.23 Accordingly, we compute a Delaunay triangulation, which involves the edges connecting the neighboring pixels of the computed mask, as shown in Fig. 1.

Fig. 1

Pairwise sparse MRF connected through the Delaunay triangulation graph. The circles represent Delaunay triangulation nodes and the solid lines represent Delaunay triangulation edges. p denotes a node in a Delaunay triangulation and q denotes the neighbor node of p on the Delaunay triangulation edge. The local adaptation of weights wpq expresses prior knowledge of phase variability on the edge from node p to q.

JARS_12_1_015006_f001.png

2.2.

Sparse Markov Random Field Model Representation

The wrapped phase φ and the absolute unwrapped phase ϕ are related by multiples of 2π, i.e.,

(1)

ϕ=φ+2πk,
where the wrap count k is an integer. The PU problem determines the unknown wrap count for each pixel. The energy function for PU based on sparse MRF can be defined as follows:

(2)

E(kp|φ)=pP,qNpV(Δϕpq)·wpq,
where P are the set of nodes in the Delaunay triangulation and Np is the neighborhood system of the Delaunay triangulation node p. φp and kp are the wrapped phase and the wrap count at Delaunay triangulation node p; q is the neighbor nodes of p on the Delaunay triangulation edge; Δϕpq is the phase variability, which is the pairwise difference of the unwrapped phase between node p and neighbor q; and V(·) is the edge potential, which is a real-valued function. We use the classical Lp norm, which is the most widely used class of clique potentials in PU; it is given by V(Δϕ)=|ΔϕW(Δϕ)|p, where W(·) is termed the wrapping operator. As defined in Ref. 13, edge potential is convex given that p1, while it is nonconvex given that p<1. It is worth mentioning that convex potential has priority over nonconvex potential both in stability and computational efficiency. Here, we choose p=1, which exactly solves the classical minimum L1 norm PU problem. wpq is the local adaptation of weights, which expresses prior knowledge on phase variability between node p and its neighbor q. Here, we set the value of wpq as the mean of coherence coefficient γp of the node p and γq of the neighbor q, which is given by

(3)

wpq=(γp+γq)2.

This optimization problem can be solved by graph cuts optimization, which is explained in Sec. 2.3. Initially, the labels of all nodes are set to zero, i.e., kpt=0=0, for all p. At each iteration step, every node’s label will either be 1 (phase plus 2π) or 1 (phase minus 2π) or 0 (phase remains unchanged), i.e., kpt+1=kpt+xpt+1, in which the superscript t denotes iteration and xpt+1{0,1,1}. Every iteration aims to decrease the value of the energy function as much as possible. When the energy ceases to decrease, the iteration is terminated. For a pair of neighboring nodes, the edge potential to be minimized is defined as follows:

(4)

Epq(xpt+1,xqt+1)=|ϕpt+1ϕqt+1W(ϕpt+1ϕqt+1)|·wpq,

(5)

ϕpt+1=φp+2π(kpt+xpt+1),

(6)

ϕqt+1=φq+2π(kqt+xqt+1),
where ϕpt+1 and ϕqt+1 are the unwrapped phase values of node p and neighbor q during iteration, respectively. Considering all pairs of nodes, the energy function of each iteration is given by

(7)

minxE(x)=pqEpq(xpt+1,xqt+1),
where the summation is over all edges. The overall process of PU based on sparse MRF is shown in Fig. 2.

Fig. 2

Flowchart of the proposed sparse MRF method for a 2-D PU problem.

JARS_12_1_015006_f002.png

2.3.

Graph Cuts Optimization for Sparse Markov Random Fields

Because the edge potentials of Epq(xpt+1,xqt+1) are convex, this optimization problem can be solved by the standard graph cuts method. The graph cuts algorithm minimizes Eq. (7) through a sequence of binary optimizations, with each binary problem mapped onto a certain graph ς=υ,ξ (υ and ξ denote the set of nodes and edges, respectively) with two special terminals α and β. Each node p is connected to the terminals α and β by edges tpα and tpβ (referred to as a t-links), respectively. Each pair of nodes {p,q} is connected by an edge epq (referred to as an n-link). The set of edges consists of tpα, tpβ, and epq. The binary minimization is obtained by computing a minimum cut on the graph. Any cut leaves each node p with exactly only one t-link. If the cut separates p from terminal α, then node p is assigned label α. Similarly, node p is assigned label β when the cut separates node p from terminal β.

To make the representation of the 2-D PU problem, we design dual elementary graphs shown in Figs. 3(a) and 3(b), respectively. The first one with terminal α set as 1 (phase plus 2π) and β set as 0 (phase remains unchanged) is illustrated in Fig. 3(a), and the weights assigned to the edges of the elementary graph are defined in Eq. (8). The second one with terminal α as 1 (phase minus 2π) and β as 0 (phase remains unchanged) is illustrated in Fig. 3(b), and the weights assigned to the edges of the elementary graph are defined in Eq. (9). We merge the two elementary graphs to obtain the Delaunay triangle graph, which is illustrated in Fig. 3(c). The Delaunay triangle graph is used to minimize Eq. (7) by choosing two of the possible swap, 0 and 1 or 0 and 1, which makes the process more efficient. Using the min-cut/max-flow formulation, the optimal swap for the entire graph can be computed:

(8)

{|tpα|=E(1,0)E(0,0),if  [E(1,0)E(0,0)]>0|tpβ|=E(0,0)E(1,0),if  [E(1,0)E(0,0)]0|tqα|=E(1,1)E(1,0),if  [E(1,1)E(1,0)]>0|tqβ|=E(1,0)E(1,1),if  [E(1,1)E(1,0)]0|epq|=E(0,1)+E(1,0)E(0,0)E(1,1),

(9)

{|tpα|=E(1,0)E(0,0),if  [E(1,0)E(0,0)]>0|tpβ|=E(0,0)E(1,0),if  [E(1,0)E(0,0)]0|tqα|=E(1,1)E(1,0),if  [E(1,1)E(1,0)]>0|tqβ|=E(1,0)E(1,1),if  [E(1,1)E(1,0)]0|epq|=E(0,1)+E(1,0)E(0,0)E(1,1).

Fig. 3

The dual elementary graphs are designed to minimize Eq. (7). (a) The first elementary graph for energy function is constructed, where 1 and 0 represent two terminals. In this case, E(1,0)E(0,0)>0 and E(1,1)E(1,0)>0. (b) The second elementary graph for energy function is constructed, where 1 and 0 represent two terminals. In this case, E(1,0)E(0,0)>0 and E(1,1)E(1,0)>0. (c) The Delaunay triangle graph is obtained in the end results from the dual elementary graphs.

JARS_12_1_015006_f003.png

3.

Simulated Data Experiment

In order to analyze the accuracy of the proposed MRF algorithm, simulated data are used in the experiment, which are obtained from an interferogram simulated by the MATLAB toolbox for InSAR.24 The function simulates an interferogram by radarcoding a fractal DEM and uses the specified percentage to create water areas of height 0 (approximately). The phase of water areas is uniform noise, and the coherence is Gaussian below 0.2. The geometric decorrelation noise is modeled from the terrain slopes according to the fractal dimension of the simulated DEM and added to the simulated interferogram. Smaller fractal dimension implies smoother surface while higher fractal dimension means steeper topography. In order to evaluate the performance of the proposed MRF method against two kinds of decorrelation effects (one is caused by large low-coherence region and the other is rapid-topography variation induce geometric decorrelation), two simulated data (A and B) are selected for analysis. Major parameters of the simulated data A and B are listed in Table 1.

Table 1

Major parameters of the simulated data A and B.

DatasetSatellite height (km)Looking angleWavelength (m)Perpendicular baseline (m)Water area percentage (%)Fractal dimension
Simulated data A785190.057100502.1
Simulated data B785190.057100102.4

In order to compare the performance of the proposed MRF algorithm to other existing ones, the traditional regular-grid MRF (RG-MRF) algorithm is used for comparison to verify the improvement of the proposed algorithm, and the Delaunay-MCF algorithm, which also deals with sparse data, serves as the competing algorithms for comparison. The RG-MRF algorithm used here is the same as in Ref. 13, in which the energy function is defined as follows:

(10)

E(k|φ)=i,jPV(Δϕijh)·vij+V(Δϕijv)·hij,
where (·)h and (·)v denote pixel horizontal and vertical phase differences, and the L1 norm potential is used, which is the same as the proposed algorithm. hij and vij are weights, which represent the local coherence coefficient at horizontal and vertical edge. Moreover, some improvements on unwrapped accuracy can also be made by using the RG-MRF method with mask (RGM-MRF), if we set hij and vij=0 for low-coherence areas, which is also used for comparison. The Delaunay-MCF method is used to compare as suggested in Ref. 17, in which the optimization objective is as follows:

(11)

ϵ=i=1Nwi·|Δϕα+Δϕβ+Δϕγ2π|,
where Δϕα, Δϕβ, and Δϕγ represent the phase gradient at each edge of a generic triangle, respectively (whose edges are labeled α, β, and γ, respectively). N represents the overall number of edges relevant to the Delaunay triangulation, and the weight wi is the local coherence coefficient at each edge. To quantitively analyze the performance of each algorithm, the root mean square (RMS) error of the unwrapping accuracy is defined as follows:

(12)

RMS=i=1N(ϕ^ϕ)2N,
where ϕ is the vector collecting from the reference true unwrapped phase, ϕ^ is the vector collecting from the estimated unwrapped phase generated by each algorithm. N represents the overall number of high-coherence pixels. It is also worth mentioning that, to fairly compare the PU result, RMS error of each algorithm is calculated based on the same sparse data. We implemented the algorithms in C or C++ and ran the experiments on a modern Pentium 5. All of the experiments were run on the same machine (2.5 GHz and 8 Gbyte RAM).

3.1.

Experiment on the Simulated Data A

Simulated data A has the characteristic that the area of low-coherence region is large but the topography is comparatively flat. Figures 4(a)4(c) show the unwrapped phase image, the wrapped phase image, and the coherence coefficient map, respectively. In Fig. 4(b), it can be noticed that the wrapped phase image is dominated by uniform phase noise induced by the water areas. In Fig. 4(c), we can see that the coherence coefficient of the water areas is below 0.2, which normally makes unwrapping the rest of the data unfeasible. Figure 4(d) shows that 147,823 spatially coherent pixels are generated, where pixels with a coherence coefficient lower than 0.2 are masked out. The Delaunay triangulation is computed, which involves 442,157 edges connecting the neighboring pixels of the computed mask, as shown in Fig. 4(e).

Fig. 4

Test on the simulated data A with a size of 512×512  pixels. (a) True interferometric phase image, (b) flattened interferogram, (c) coherence coefficient map, (d) mask map of the spatially coherent pixels where pixels with coherence higher than 0.2 (in red), and (e) Delaunay triangulation involving the set of spatially coherent pixels.

JARS_12_1_015006_f004.png

Figure 5 illustrates the unwrapped results of simulated data A by each algorithm. In Fig. 5, we can visually notice that among them, the proposed MRF method has the deviation error in a smaller scale and very close to it is the Delaunay-MCF method. Otherwise, the deviation errors of the two regular-grid MRF algorithms are higher than those of the former two, especially the RGM-MRF algorithm, which generates an unreasonable unwrapped result. To quantitively analyze the accuracy and efficiency of each algorithm, the RMS error and running time are calculated and listed in Table 2. Furthermore, it can also be noticed that the proposed MRF algorithm generates the lowest RMS error, followed by the Delaunay-MCF algorithm. On the contrary, the result of the two regular-grid MRF methods has higher RMS errors compared to the other two methods. This is because they are both based on the four-connected grid graph, different isolated regions surrounded by low-coherence areas cannot be connected to each other, which causes phase jumps between them easily. In addition to this, in Table 2, it is also noticed that although the proposed MRF algorithm is more time consuming than the Delaunay-MCF algorithm, it is computationally efficient compared with the other two MRF algorithms, and the reason is that it is based on a set of 2-D Delaunay triangles instead of the entire image, which drastically reduces the running time of the optimization process.

Fig. 5

Unwrapped results on the simulated data A using the (a) proposed MRF algorithm, (b) RG-MRF algorithm, (c) RGM-MRF algorithm, (d) Delaunay-MCF algorithm. Deviation errors between the unwrapped result generated by the (e) proposed MRF algorithm, (f) RG-MRF algorithm, (g) RGM-MRF algorithm, (h) Delaunay-MCF algorithm, and the reference topography phase.

JARS_12_1_015006_f005.png

Table 2

Rms error and running time of each algorithm on the simulated data A.

AlgorithmsTime (s)RMS error (rad)
Proposed MRF14.192.15
RG-MRF24.922.8
RGM-MRF18.1414.36
Delaunay-MCF10.242.25

In order to research into the performance of each algorithm more deeply, we calculate the RMS error against the water percentage of the simulated DEM. The higher water percentage means that more area in the simulated interferogram has the low coherence coefficient level. Figure 6 shows the relationship between the RMS error and water percentage (from 10% to 70%). In Fig. 6, it can be seen that the proposed MRF algorithm (red line) generates lower RMS error compared with the competing algorithms in the majority of the water percentage scale and very close to it is the Delaunay-MCF algorithm (green line). The RGM-MRF algorithm (cyan line), on the other hand, suffers the highest RMS error, followed by the RG-MRF algorithm (blue line). Moreover, it can still be noticed that the error curves of the four algorithms are similar with a low-level RMS error in the area with low water percentage (below 40%). However, in the area with high water percentage (above 40%), the proposed MRF algorithm generates lower RMS error than the competing algorithms, which means that it suffers less from decorrelation effects caused by expanding of the low-coherence area.

Fig. 6

Relationship between the unwrapped RMS error (rad) of each algorithm and the water percentage of the simulated DEM (%).

JARS_12_1_015006_f006.png

3.2.

Experiment on the Simulated Data B

Simulated data B has the characteristic that the topography varies more drastically, which suffers serious geometric decorrelation effect. Figures 7(a) and 7(b) show the unwrapped phase image and the wrapped phase image, respectively. In Fig. 7(b), it can be noticed that the wrapped phase is rather noisy, thus inducing a huge number of phase jumps (i.e., residues), making the unwrapping a hard task. Figure 7(c) shows the coherence coefficient map, which is computed from the terrain slopes (geometric decorrelation). Figure 7(d) shows that 239,870 spatially coherent pixels are generated, where pixels with a coherence coefficient lower than 0.5 are masked out. The Delaunay triangulation is computed, which involves 717,804 edges connecting the neighboring pixels of the computed mask, as shown in Fig. 7(e).

Fig. 7

Test on the simulated data B with a size of 512×512  pixels. (a) True interferometric phase image, (b) flattened interferogram, (c) coherence coefficient map, (d) mask map of the spatially coherent pixels where pixels with coherence higher than 0.5 (in red), and (e) Delaunay triangulation involving the set of spatially coherent pixels.

JARS_12_1_015006_f007.png

Figure 8 shows the unwrapped results of simulated data B by each algorithm. In Fig. 8, it is illustrated visually that the unwrapped result using the proposed MRF and RGM-MRF algorithms has lower deviation error than that of the Delaunay-MCF and RG-MRF algorithms. Table 3 shows the corresponding RMS errors and running time of each algorithm. Among them, the proposed MRF algorithm suffers from the lowest RMS error. For the Delaunay-MCF algorithm, on the contrary, it abnormally generates the highest RMS error. Moreover, the performance of the two regular-grid MRF algorithms is very distinctive from each other. The RG-MRF algorithm generates very high RMS error, while the RGM-MRF algorithm has much lower RMS error and the reason may be that low-coherence pixels are excluded during the PU process. In addition, Table 3 also demonstrates that the proposed MRF algorithm is more efficient than other two MRF algorithms.

Fig. 8

Unwrapped results on the simulated data B using the (a) proposed MRF algorithm, (b) RG-MRF algorithm, (c) RGM-MRF algorithm, and (d) Delaunay-MCF algorithm. Deviation errors between the unwrapped result generated by the (e) proposed MRF algorithm, (f) RG-MRF algorithm, (g) RGM-MRF algorithm, and (h) Delaunay-MCF algorithm and the reference topography phase.

JARS_12_1_015006_f008.png

Table 3

RMS error and running time of each algorithm on the simulated data B.

AlgorithmsTime (s)RMS error (rad)
Proposed MRF15.201.49
RG-MRF24.365.72
RGM-MRF17.981.57
Delaunay-MCF11.845.80

Here, we calculate the RMS error against the fractal dimension of the simulated DEM for further analysis. The higher fractal dimension means local topography varies more drastically, which results in low-coherence induced by geometric decorrelation. The relationship between the RMS error and fractal dimension (between 2 and 3) is shown in Fig. 9. It can be observed that, in the area with low fractal dimension (below 2.3), meaning that the topography is comparatively flat, the four algorithms have identical error curves and have a low-level RMS error. However, in the area with high fractal dimension (above 2.3), meaning that the topography varies drastically, the behavior of these curves is very distinctive from each other. The Delaunay-MCF algorithm (green line) generates the highest RMS error, followed by the RG-MRF algorithm (blue line). On the opposite, the proposed MRF algorithm (red line) introduces the lowest RMS error in the majority of the fractal dimension scale and very close to it is the RGM-MRF algorithm (cyan line), which demonstrates that the performance of the proposed MRF algorithm suffers less from decorrelation effects caused by topography variation.

Fig. 9

Relationship between the unwrapped RMS error (rad) of each algorithm and the fractal dimension (between 2 and 3) of the simulated DEM.

JARS_12_1_015006_f009.png

4.

Real Data Experiment

In this section, the performance of the proposed method is tested on two real applications, which both intend to reconstruct the real DEM. The two experiments are applied to two interferometric pairs of TerraSAR-X, respectively, which are both downloaded from the Infoterra website. The major parameters of the two interferometric pairs are listed in Table 4. Moreover, the Advanced Spaceborne Thermal Emission and Reflection Radiometer Global (ASTER G) DEM covering the same area is used to simulate the topography phase, which is, as the reference, compared with every unwrapped result generated by different algorithms.

Table 4

Parameters of two interferometric pairs of TerraSAR-X.

DatasetWavelength (m)Resolution (m)PolarizationOrbitLooking angleDate
Uluru0.0311HHDescending45.8February 12, 2009
February 23, 2009
Grand0.0313HHDescending39.2March 10, 2008
CanyonMarch 21, 2008

4.1.

Experiment on the Real Data of the Uluru

In the first real data experiment, the performances of the four methods are tested on two TerraSAR-X spotlight images covering the Uluru, in Australia. Figure 10(a) shows the SAR intensity image of this data, which is geographically featured by an isolated hill that rises abruptly from and is surrounded by relatively flat lowlands. Figure 10(b) shows the topography phase simulated by ASTER G DEM, which is used as the reference topography phase. We use the two SLC images of the Uluru to generate a flattened interferogram [see Fig. 10(c)] with eight looks in the range and azimuth direction, having a pixel spacing of approximately 8×8  m. The Goldstein interferogram filter is applied to reduce the phase noise before the competing unwrapper. In Fig. 10(c), we can notice that it presents large phase rates to produce aliasing between isolated hill and flat lowlands; in addition, the back of the isolated hill has the characteristic that the coherence coefficient is low due to the shadow, such that the unwrapping becomes a hard task. Figure 10(d) shows the coherence coefficient map using the adaptive estimation window size of 7×7. We generate the mask map of the pixels by considering those pixels with an estimated coherence value greater than a selected threshold 0.3, as shown in Fig. 10(e). The number of selected coherent sparse pixels is 542,392. The Delaunay triangulation is then generated to define the neighbor nodes in the set of the identified coherent sparse pixels, which is shown in Fig. 10(f). The number of edges relevant to the Delaunay triangulation is 1,624,132.

Fig. 10

Test on the real data of the Uluru, which is a TerraSAR-X phase image with a size of 1274×680  pixels. (a) SAR intensity image, (b) topography phase simulated by ASTER G DEM, which is used as the reference topography phase, (c) flattened filtered interferogram, (d) coherence coefficient map, (e) mask map of the spatially coherent pixels where pixels with coherence higher than 0.3 (in red), and (f) Delaunay triangulation involving the set of spatially coherent pixels.

JARS_12_1_015006_f010.png

Figure 11 shows the unwrapped results obtained by the four algorithms, and Table 5 shows the corresponding RMS errors and running time of each algorithm. In Fig. 11, we can notice that one obvious vertical discontinuity is seen in the unwrapped result generated by the RG-MRF algorithm, and the unwrapped result obtained by RGM-MRF algorithm shows that the discontinuities problem becomes much more serious. In Fig. 11, it is also shown visually that the proposed MRF and Delaunay-MFC algorithms both alleviate the discontinuities successfully and keep the continuity of the unwrapped result in the whole image. The relative performance of the four algorithms can also be seen in the metrics in Table 5, where the proposed MRF and Delaunay-MCF algorithms suffer less from RMS error in comparison with the reference, while the two regular-grid MRF algorithms have a higher RMS error level than that of the former two, and the reason is given in the end of Sec. 3.1. In Table 5, we can see that the proposed MRF algorithm’s execution time is smaller than that needed in the other two MRF algorithms.

Fig. 11

Unwrapped results on the real data of the Uluru using the (a) proposed MRF algorithm, (b) RG-MRF algorithm, (c) RGM-MRF algorithm, (d) Delaunay-MCF algorithm. Deviation errors between the unwrapped result generated by the (e) proposed MRF algorithm, (f) RG-MRF algorithm, (g) RGM-MRF algorithm, and (h) Delaunay-MCF algorithm and the reference topography phase.

JARS_12_1_015006_f011.png

Table 5

RMS error and runtime of each algorithm on the real data of the Uluru.

AlgorithmsTime (s)RMS error (rad)
Proposed MRF86.673.14
RG-MRF235.465.31
RGM-MRF120.5411.91
Delaunay-MCF52.453.88

4.2.

Experiment on the Real Data of the Grand Canyon

In the second real data experiment, the performances of the four methods are tested on two TerraSAR-X stripmap images covering the Grand Canyon, in the United States. We remark that the selection of this Grand Canyon test-site area is related to the difficulty in the unwrapping operation, where some regions have the characteristic that the coherence coefficient is very low and the topography varies drastically since it is located across the canyon. The SAR intensity image of this data is illustrated in Fig. 12(a). Figure 12(b) shows the topography phase simulated by ASTER G DEM, which is used as the reference topography phase. Figure 12(c) shows the flattened interferogram with 10 looks in the range and azimuth direction, with a pixel spacing of approximately 30×30  m. Before processing two algorithms, a Goldstein interferogram filter is used for phase denoising. Figure 12(d) shows the coherence coefficient map using the adaptive estimation window size of 7×7. Figure 12(e) shows that 969,868 spatially coherent pixels are generated, where pixels with a coherence coefficient lower than 0.3 are masked out. The Delaunay triangulation is computed, which involves 5,813,674 edges connecting the neighboring pixels of the computed mask, as shown in Fig. 12(f).

Fig. 12

Test on the real data of the Grand Canyon, which is a TerraSAR-X phase image with a size of 1171×1283  pixels. (a) SAR intensity image, (b) topography phase simulated by ASTER G DEM, which is used as the reference topography phase, (c) flattened filtered interferogram, (d) coherence coefficient map, (e) mask map of the spatially coherent pixels where pixels with coherence higher than 0.3 (in red), and (f) Delaunay triangulation involving the set of spatially coherent pixels.

JARS_12_1_015006_f012.png

The unwrapped results obtained by the four algorithms are illustrated in Fig. 13, and Table 6 shows the corresponding RMS errors and running time of each algorithm. In Fig. 13, the unwrapped result by the proposed MRF and RGM-MRF algorithms seems more credible and seamless than that of the other two algorithms, while the unwrapped result by the RG-MRF and Delaunay-MCF algorithm has the 2nπ unwrapping errors accumulated throughout the whole image. As indicated by the metrics in Table 6, the proposed MRF and RGM-MRF algorithms perform better in generating the lower RMS error, while the Delaunay-MCF algorithm abnormally generates the highest RMS error, and close to it is the RG-MRF algorithm. This is partly because topography varies drastically around the canyon while the phase fringes break and generate a lot of residues, which degrades the accuracy of the Delaunay-MCF and RG-MRF algorithms. On the contrary, the proposed MRF algorithm suffers less from the decorrelation effects caused by topography variation. Table 6 also demonstrates that the proposed MRF algorithm is more efficient in comparison with the other two MRF algorithms.

Fig. 13

Unwrapped results on the real data of the Grand Canyon using the (a) proposed MRF algorithm, (b) RG-MRF algorithm, (c) RGM-MRF algorithm, (d) Delaunay-MCF algorithm. Deviation errors between the unwrapped result generated by the (e) proposed MRF algorithm, (f) RG-MRF algorithm, (g) RGM-MRF algorithm, and (h) Delaunay-MCF algorithm and the reference topography phase.

JARS_12_1_015006_f013.png

Table 6

RMS error and runtime of each algorithm on the real data of the Grand Canyon.

AlgorithmsTime (s)Rms error (rad)
Proposed MRF316.3110.37
RG-MRF658.9720.01
RGM-MRF430.4310.48
Delaunay-MCF250.2530.66

Finally, we give the comparison between simulated and real data experiments. In the experiment on the simulated data A, interferograms are simulated with lots of low-coherence regions caused by water areas. In this case, it can be noticed that the proposed algorithm generates the lowest error level, with the Delaunay-MCF algorithm resulting in a slightly larger error. However, the two regular-grid algorithms have the high error level. The relative performance of the four algorithms is also found in the real-data experimental result of the Uluru, where the area of the decorrelation region is also large due to shadow. The two experimental results both demonstrate that the proposed algorithm suffers less from the influence of expanding of the low-coherence area. In the experiment on the simulated data B, interferograms with lots of low coherence caused by rapid-topography variation are simulated. In this condition, it generates the similarity between the proposed and RGM-MRF algorithms, and they both have lower error level compared with other two algorithms. The similar results also can be noticed in the real-data experimental result of the Grand Canyon, where the terrain changes fiercely across the canyon. The two experimental results both confirm that the proposed algorithm suffers less decorrelation effects caused by topography variation.

5.

Conclusions

We proposed a PU solution based on sparse MRF for extending the traditional regular-grid MRF PU algorithm dealing with a sparse data to process InSAR interferograms for the generation of DEM. Our method is inspired by the fact that areas of low interferometric correlation are affected by high phase noise and therefore are excluded during the PU process. The proposed algorithm involves three main steps. Initially, Delaunay triangulation is computed based on a sparse set of coherent points using a coherence coefficient of the SAR images. Following this step, a sparse MRF representation based on the Delaunay triangulation is defined. Finally, to speed up the graph cut algorithm for sparse MRF, we design dual elementary graphs and merge them to obtain the Delaunay triangle graph, which is used to minimize energy function efficiently. The experiments are carried out both on simulated and real data, compared with the other existing algorithms. All the experimental results confirm that the proposed MRF method is an effective PU method, which suffers less from decorrelation effects caused by a large low-coherence area and rapid-topography variation. We underline that the proposed solution is based on a set of 2-D Delaunay triangles instead of the entire image, which reduces the running time of the optimization process drastically. Accordingly, the overall approach is computationally efficient.

It is worth mentioning that we focus only on DEM extraction using the proposed algorithm, but InSAR is also useful to measure ground deformation, for instance, due to subsidence (volcanic or fluid extraction/injection) or coseismic/postseismic deformation. In the future, we will research the application of the proposed algorithm for those problems, especially in coseismic deformation, where the interferograms suffer from large gradients in ground deformation near the fault rupture and discontinuities in phase across the rupture. Furthermore, although the algorithm is developed mainly for the single-pass interferometry or short-period singlebaseline interferometry, which is used for the generation of topography or deformation map, the proposed algorithm can be extended into the three dimensions (3-D) to solve the multibaseline interferometry PU problem. In a future study, we will utilize 3-D Delaunay triangulation to extend this algorithm for the interferogram series reconstruction problem.

Acknowledgments

This research was jointly supported by the National Natural Science Foundation of China (Grant No. 41501461), the sponsorship of Jiangsu Overseas Research and Training Program for University Prominent Young and Middle-aged Teachers and Presidents, the Natural Science Foundation of Jiangsu Province of China (Grant No. BK20140419). The authors would like to thank the anonymous reviewers for their constructive comments, which were taken into account in revising this paper.

References

1. H. Zebker and R. Goldstein, “Topographic mapping from interferometric SAR observations,” J. Geophys. Res. 91, 4993–4999 (1986). http://dx.doi.org/10.1029/JB091iB05p04993 Google Scholar

2. K. Itoh, “Analysis of the phase unwrapping problem,” Appl. Opt. 21(14), 2470 (1982).APOPAI0003-6935 http://dx.doi.org/10.1364/AO.21.002470 Google Scholar

3. S. Madesn, H. Zebker and J. Martin, “Topographic mapping using radar interferometry: processing techniques,” IEEE Trans. Geosci. Remote Sens. 31(1), 246–256 (1993).IGRSD20196-2892 http://dx.doi.org/10.1109/36.210464 Google Scholar

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

5. W. Xu and I. Cunning, “A region growing algorithm for InSAR phase unwrapping,” IEEE Trans. Geosci. Remote Sens. 37(1), 124–134 (1991).IGRSD20196-2892 http://dx.doi.org/10.1109/36.739143 Google Scholar

6. T. Flynn, “Two-dimensional phase unwrapping with minimum weighted discontinuity,” J. Opt. Soc. Am. A 14(10), 2692–2701 (1997).JOAOD60740-3232 http://dx.doi.org/10.1364/JOSAA.14.002692 Google Scholar

7. D. Fried, “Least-squares fitting a wave-front distortion estimate to an array of phase difference measurements,” J. Opt. Soc. Am. A 67(3), 370–375 (1977).JOAOD60740-3232 http://dx.doi.org/10.1364/JOSA.67.000370 Google Scholar

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

9. J. Marroquin and M. Rivera, “Quadratic regularization functionals for phase unwrapping,” J. Opt. Soc. Am. A 12(11), 2393–2400 (1995).JOAOD60740-3232 http://dx.doi.org/10.1364/JOSAA.12.002393 Google Scholar

10. G. Nico, G. Palubinskas and M. Datcu, “Bayesian approach to phase unwrapping: theoretical study,” IEEE Trans. Image Process. 48(9), 2545–2556 (2000).IIPRE41057-7149 http://dx.doi.org/10.1109/78.863057 Google Scholar

11. J. Dias and J. Leitao, “The ZM algorithm for interferometric image reconstruction in SAR/SAS,” IEEE Trans. Image Process. 12(4), 408–422 (2002).IIPRE41057-7149 http://dx.doi.org/10.1109/TIP.2002.999675 Google Scholar

12. L. Ying et al., “Unwrapping of MR phase images using a Markov random field model,” IEEE Trans. Med. Imaging 25(1), 128–136 (2006).ITMID40278-0062 http://dx.doi.org/10.1109/TMI.2005.861021 Google Scholar

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

14. R. Chen et al., “Integrated denoising and unwrapping of InSAR phase based on Markov random fields,” IEEE Trans. Geosci. Remote Sens. 51(8), 4473–4485 (2013).IGRSD20196-2892 http://dx.doi.org/10.1109/TGRS.2013.2268969 Google Scholar

15. J. Frey, R. Koetter and N. Petrovic, “Very loopy belief propagation for unwrapping phase images,” in Advances in Neural Information Processing Systems (NIPS), pp. 737–743 (2001). Google Scholar

16. M. Costantini and P. A. Rosen, “A generalized phase unwrapping approach for sparse data,” in Proc. IEEE Int. Geoscience and Remote Sensing Symp., Hamburg (1999). http://dx.doi.org/10.1109/IGARSS.1999.773467 Google Scholar

17. A. Pepe and R. Lanari, “On the extension of the minimum cost flow algorithm for phase unwrapping of multitemporal differential SAR interferograms,” IEEE Trans. Geosci. Remote Sens. 44(9), 2374–2383 (2006).IGRSD20196-2892 http://dx.doi.org/10.1109/TGRS.2006.873207 Google Scholar

18. A. Shanker and H. Zebker, “Edgelist phase unwrapping algorithm for time series InSAR analysis,” J. Opt. Soc. Am. A 27(3), 605–612 (2010).JOAOD60740-3232 http://dx.doi.org/10.1364/JOSAA.27.000605 Google Scholar

19. A. Ferretti, C. Prati and F. Rocca, “Permanent scatterers in SAR interferometry,” IEEE Trans. Geosci. Remote Sens. 39(1), 8–20 (2001).IGRSD20196-2892 http://dx.doi.org/10.1109/36.898661 Google Scholar

20. A. Hooper et al., “A new method for measuring deformation on volcanoes and other natural terrains using InSAR persistent scatters,” Geophys. Res. Lett. 31, L23611 (2004).GPRLAJ0094-8276 http://dx.doi.org/10.1029/2004GL021737 Google Scholar

21. B. Kampes and R. Hanssen, “Ambiguity resolution for permanent scatterer interferometry,” IEEE Trans. Geosci. Remote Sens. 42(11), 2446–2453 (2004). http://dx.doi.org/10.1109/TGRS.2004.835222 Google Scholar

22. R. Hanssen, Radar Interferometry—Data Interpretation and Error Analysis, Kluwer, Norwell (2002). Google Scholar

23. L. Guibas and J. Stolfi, “Primitives for the manipulation of general subdivisions and the computation of Voronoi diagrams,” ACM Trans. Graphics 4(2), 74–123 (1985).ATGRDF0730-0301 http://dx.doi.org/10.1145/282918.282923 Google Scholar

24. B. Kampes and S. Usai, “Doris: the Delft object-oriented radar interferometric software,” in 2nd Int. Symp. on Operationalization of Remote Sensing, Enschede (1999). Google Scholar

Biography

Lifan Zhou is a lecturer at the Changshu Institute of Technology. He received his BS degree in geographic information systems from Wuhan University in 2006 and his PhD degree in geographic information systems from Zhejiang University in 2014. He is the author of more than 10 journal papers. His main research interests are in the fields of phase unwrapping and synthetic aperture radar interferometry signal processing and applications.

Dengfeng Chai received his BS degree in surveying engineering from Wuhan University, his MS degree in photogrammetry and remote sensing from Wuhan University, and his PhD in applied mathematics from Zhejiang University in 1997, 2000, and 2006, respectively. He is an associate professor at Zhejiang University. He is the author of more than 20 journal papers. His research interests include developing statistical approaches for object recognition and extraction from remotely sensed images.

Yu Xia is an associate professor at the Changshu Institute of Technology. He received his MS degree in computer science and PhD in computer vision both from Jiangnan University in 2009 and 2013, respectively. He is the author of more than 20 journal papers. His research interests include computer vision, digital image processing, and pattern recognition.

Peifeng Ma is a postdoctor at the Chinese University of Hong Kong. He received his MS degree in signal and information processing from the Chinese Academy of Sciences in 2012, and the PhD in earth system and geoinformation science from the Chinese University of Hong Kong in 2016. He is the author of more than 10 journal papers. His research interests include SAR tomography and persistent scatterer interferometry and their application to infrastructural health monitoring.

Hui Lin is a professor at the Chinese University of Hong Kong. He received his MS degree in cartography and remote sensing from the Chinese Academy of Sciences in 1983, and his MA and PhD degrees in geographic information systems from the University at Buffalo in 1987 and 1992, respectively. He is the author of more than 200 journal papers. His research interests are in satellite remote sensing for cloudy and rainy environments and GIS applications.

© The Authors. Published by SPIE under a Creative Commons Attribution 3.0 Unported License. Distribution or reproduction of this work in whole or in part requires full attribution of the original publication, including its DOI.
Lifan Zhou, Lifan Zhou, Dengfeng Chai, Dengfeng Chai, Yu Xia, Yu Xia, Peifeng Ma, Peifeng Ma, Hui Lin, Hui Lin, } "Interferometric synthetic aperture radar phase unwrapping based on sparse Markov random fields by graph cuts," Journal of Applied Remote Sensing 12(1), 015006 (16 January 2018). https://doi.org/10.1117/1.JRS.12.015006 . Submission: Received: 31 July 2017; Accepted: 20 December 2017
Received: 31 July 2017; Accepted: 20 December 2017; Published: 16 January 2018
JOURNAL ARTICLE
19 PAGES


SHARE
RELATED CONTENT


Back to Top