|
1.IntroductionInterferometric 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- 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 -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 norm of the difference between absolute phase differences and wrapped phase differences is minimized. When , the norm-based PU method is transformed to the least squares (LS) PU method,7 and when , 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 norm-based criterion performs better than the norm in preserving discontinuities.4 However, the MCF approach, as a -true unwrapper, introduces integer times of 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–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–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.16–18 What is more, multitemporal/multibaseline InSAR (i.e., persistent scatterers interferometry, PSInSAR) attracted increasing interests,19–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 UnwrappingIn 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 TriangulationWe 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. 2.2.Sparse Markov Random Field Model RepresentationThe wrapped phase and the absolute unwrapped phase are related by multiples of , i.e., where the wrap count 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: where are the set of nodes in the Delaunay triangulation and is the neighborhood system of the Delaunay triangulation node . and are the wrapped phase and the wrap count at Delaunay triangulation node ; is the neighbor nodes of on the Delaunay triangulation edge; is the phase variability, which is the pairwise difference of the unwrapped phase between node and neighbor ; and is the edge potential, which is a real-valued function. We use the classical norm, which is the most widely used class of clique potentials in PU; it is given by , where is termed the wrapping operator. As defined in Ref. 13, edge potential is convex given that , while it is nonconvex given that . It is worth mentioning that convex potential has priority over nonconvex potential both in stability and computational efficiency. Here, we choose , which exactly solves the classical minimum norm PU problem. is the local adaptation of weights, which expresses prior knowledge on phase variability between node and its neighbor . Here, we set the value of as the mean of coherence coefficient of the node and of the neighbor , which is given byThis 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., , for all . At each iteration step, every node’s label will either be 1 (phase plus ) or (phase minus ) or 0 (phase remains unchanged), i.e., , in which the superscript denotes iteration and . 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: where and are the unwrapped phase values of node and neighbor during iteration, respectively. Considering all pairs of nodes, the energy function of each iteration is given by where the summation is over all edges. The overall process of PU based on sparse MRF is shown in Fig. 2.2.3.Graph Cuts Optimization for Sparse Markov Random FieldsBecause the edge potentials of 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 is connected to the terminals and by edges and (referred to as a t-links), respectively. Each pair of nodes is connected by an edge (referred to as an n-link). The set of edges consists of , , and . The binary minimization is obtained by computing a minimum cut on the graph. Any cut leaves each node with exactly only one t-link. If the cut separates from terminal , then node is assigned label . Similarly, node is assigned label when the cut separates node 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 ) 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 (phase minus 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 , which makes the process more efficient. Using the min-cut/max-flow formulation, the optimal swap for the entire graph can be computed: 3.Simulated Data ExperimentIn 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 1Major parameters of the simulated data A and B.
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: where and denote pixel horizontal and vertical phase differences, and the norm potential is used, which is the same as the proposed algorithm. and 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 and 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: where , , and represent the phase gradient at each edge of a generic triangle, respectively (whose edges are labeled , , and , respectively). represents the overall number of edges relevant to the Delaunay triangulation, and the weight 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: where is the vector collecting from the reference true unwrapped phase, is the vector collecting from the estimated unwrapped phase generated by each algorithm. 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 ASimulated 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. Table 2Rms error and running time of each algorithm on the simulated data A.
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. 3.2.Experiment on the Simulated Data BSimulated 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. Table 3RMS error and running time of each algorithm on the simulated data B.
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. 4.Real Data ExperimentIn 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 4Parameters of two interferometric pairs of TerraSAR-X.
4.1.Experiment on the Real Data of the UluruIn 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 . 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 . 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. Table 5RMS error and runtime of each algorithm on the real data of the Uluru.
4.2.Experiment on the Real Data of the Grand CanyonIn 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 . 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 . 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 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. Table 6RMS error and runtime of each algorithm on the real data of the Grand Canyon.
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.ConclusionsWe 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. AcknowledgmentsThis 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. ReferencesH. 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
K. Itoh,
“Analysis of the phase unwrapping problem,”
Appl. Opt., 21
(14), 2470
(1982). http://dx.doi.org/10.1364/AO.21.002470 APOPAI 0003-6935 Google Scholar
S. Madesn, H. Zebker and J. Martin,
“Topographic mapping using radar interferometry: processing techniques,”
IEEE Trans. Geosci. Remote Sens., 31
(1), 246
–256
(1993). http://dx.doi.org/10.1109/36.210464 IGRSD2 0196-2892 Google Scholar
D. C. Ghiglia and M. D. Pritt, Two-Dimensional Phase Unwrapping: Theory, Algorithms and Software, Wiley-Interscience, New York
(1998). Google Scholar
W. Xu and I. Cunning,
“A region growing algorithm for InSAR phase unwrapping,”
IEEE Trans. Geosci. Remote Sens., 37
(1), 124
–134
(1991). http://dx.doi.org/10.1109/36.739143 IGRSD2 0196-2892 Google Scholar
T. Flynn,
“Two-dimensional phase unwrapping with minimum weighted discontinuity,”
J. Opt. Soc. Am. A, 14
(10), 2692
–2701
(1997). http://dx.doi.org/10.1364/JOSAA.14.002692 JOAOD6 0740-3232 Google Scholar
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). http://dx.doi.org/10.1364/JOSA.67.000370 JOAOD6 0740-3232 Google Scholar
M. Costantini,
“A novel phase unwrapping method based on network programming,”
IEEE Trans. Geosci. Remote Sens., 36
(3), 813
–821
(1998). http://dx.doi.org/10.1109/36.673674 IGRSD2 0196-2892 Google Scholar
J. Marroquin and M. Rivera,
“Quadratic regularization functionals for phase unwrapping,”
J. Opt. Soc. Am. A, 12
(11), 2393
–2400
(1995). http://dx.doi.org/10.1364/JOSAA.12.002393 JOAOD6 0740-3232 Google Scholar
G. Nico, G. Palubinskas and M. Datcu,
“Bayesian approach to phase unwrapping: theoretical study,”
IEEE Trans. Image Process., 48
(9), 2545
–2556
(2000). http://dx.doi.org/10.1109/78.863057 IIPRE4 1057-7149 Google Scholar
J. Dias and J. Leitao,
“The ZM algorithm for interferometric image reconstruction in SAR/SAS,”
IEEE Trans. Image Process., 12
(4), 408
–422
(2002). http://dx.doi.org/10.1109/TIP.2002.999675 IIPRE4 1057-7149 Google Scholar
L. Ying et al.,
“Unwrapping of MR phase images using a Markov random field model,”
IEEE Trans. Med. Imaging, 25
(1), 128
–136
(2006). http://dx.doi.org/10.1109/TMI.2005.861021 ITMID4 0278-0062 Google Scholar
J. Dias and G. Valadao,
“Phase unwrapping via graph cuts,”
IEEE Trans. Image Process., 16
(3), 698
–709
(2007). http://dx.doi.org/10.1109/TIP.2006.888351 IIPRE4 1057-7149 Google Scholar
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). http://dx.doi.org/10.1109/TGRS.2013.2268969 IGRSD2 0196-2892 Google Scholar
J. Frey, R. Koetter and N. Petrovic,
“Very loopy belief propagation for unwrapping phase images,”
in Advances in Neural Information Processing Systems (NIPS),
737
–743
(2001). Google Scholar
M. Costantini and P. A. Rosen,
“A generalized phase unwrapping approach for sparse data,”
in Proc. IEEE Int. Geoscience and Remote Sensing Symp.,
(1999). http://dx.doi.org/10.1109/IGARSS.1999.773467 Google Scholar
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). http://dx.doi.org/10.1109/TGRS.2006.873207 IGRSD2 0196-2892 Google Scholar
A. Shanker and H. Zebker,
“Edgelist phase unwrapping algorithm for time series InSAR analysis,”
J. Opt. Soc. Am. A, 27
(3), 605
–612
(2010). http://dx.doi.org/10.1364/JOSAA.27.000605 JOAOD6 0740-3232 Google Scholar
A. Ferretti, C. Prati and F. Rocca,
“Permanent scatterers in SAR interferometry,”
IEEE Trans. Geosci. Remote Sens., 39
(1), 8
–20
(2001). http://dx.doi.org/10.1109/36.898661 IGRSD2 0196-2892 Google Scholar
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). http://dx.doi.org/10.1029/2004GL021737 GPRLAJ 0094-8276 Google Scholar
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
R. Hanssen, Radar Interferometry—Data Interpretation and Error Analysis, Kluwer, Norwell
(2002). Google Scholar
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). http://dx.doi.org/10.1145/282918.282923 ATGRDF 0730-0301 Google Scholar
B. Kampes and S. Usai,
“Doris: the Delft object-oriented radar interferometric software,”
in 2nd Int. Symp. on Operationalization of Remote Sensing,
(1999). Google Scholar
BiographyLifan 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. |