Interferometric synthetic aperture radar phase unwrapping based on sparse Markov random fields by graph cuts

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.

Interferometric synthetic aperture radar phase unwrapping based on sparse Markov random fields by graph cuts

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 L p norm of the difference between absolute phase differences and wrapped phase differences is minimized. When p ¼ 2, the L p norm-based PU method is transformed to the least squares (LS) PU method, 7 and when p ¼ 1, it is changed to the minimumcost 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 L 1 norm-based criterion performs better than the L 2 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. [9][10][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. [12][13][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 Valadao 13 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 cuts 13 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 wellknown 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. [16][17][18] What is more, multitemporal/multibaseline InSAR (i.e., persistent scatterers interferometry, PSInSAR) attracted increasing interests, [19][20][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 regulargrid 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.

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.

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 w pq expresses prior knowledge of phase variability on the edge from node p to q.

Sparse Markov Random Field Model Representation
The wrapped phase φ and the absolute unwrapped phase ϕ are related by multiples of 2π, i.e., E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 1 ; 1 1 6 ; 7 0 4 (1) 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: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 2 ; 1 1 6 ; 6 5 1 where P are the set of nodes in the Delaunay triangulation and N p is the neighborhood system of the Delaunay triangulation node p. φ p and k p 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 L p norm, which is the most widely used class of clique potentials in PU; it is given by VðΔϕÞ ¼ jΔϕ − WðΔϕÞj p , where Wð·Þ is termed the wrapping operator. As defined in Ref. 13, edge potential is convex given that p ≥ 1, 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 L 1 norm PU problem. w pq 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 w pq as the mean of coherence coefficient γ p of the node p and γ q of the neighbor q, which is given by E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 3 ; 1 1 6 ; 4 5 3 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., k t¼0 p ¼ 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., k tþ1 p ¼ k t p þ x tþ1 p , in which the superscript t denotes iteration and x tþ1 p ∈ f0; 1; −1g. 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: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 4 ; 1 1 6 ; 3 2 5 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 5 ; 1 1 6 ; 2 9 3 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 6 ; 1 1 6 ; 2 6 5 where ϕ tþ1 p and ϕ tþ1 q 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 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 7 ; 1 1 6 ; 2 1 4 where the summation is over all edges. The overall process of PU based on sparse MRF is shown in Fig. 2.

Graph Cuts Optimization for Sparse Markov Random Fields
Because the edge potentials of E pq ðx tþ1 p ; x tþ1 q Þ 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 ς ¼ hυ; ξi (υ 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 t α p and t β p (referred to as a t-links), respectively. Each pair of nodes fp; qg is connected by an edge e pq (referred to as an n-link). The set of edges consists of t α p , t β p , and e pq . 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: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 8 ; 1 1 6 ; 6 8 7 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 9 ; 1 1 6 ; 5 7 6

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 lowcoherence 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.
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: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 0 ; 1 1 6 ; 2 3 4 where ð·Þ h and ð·Þ v denote pixel horizontal and vertical phase differences, and the L 1 norm potential is used, which is the same as the proposed algorithm. h ij and v ij 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 h ij and v ij ¼ 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: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 1 ; 1 1 6 ; 6 7 5 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 w i 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: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 2 ; 1 1 6 ; 5 6 8 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).

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

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

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

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

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