**Advanced Multitemporal Phase Unwrapping Techniques for DInSAR Analyses**

Antonio Pepe

*Istituto per il Rilevamento Elettromagnetico dell'Ambiente (IREA), Research National Council (CNR) Italy* 

### **1. Introduction**

56 Recent Interferometry Applications in Topography and Astronomy

Abdelsalam, D. ; Magnusson, R.; & Kim, D. (2011). Single-shot dual wavelength digital holography based on polarizing separation, Applied Optics, 50, pp. 3360-3368. Abdelsalam, D.; & Kim, D. (2011). Coherent noise suppression in digital holography based on flat fielding with apodized apertures, Optics Express, 19, pp. 17951-17959. Abdelsalam, D.; Baek, B.; Cho, Y.; & Kim, D. (2010). Surface form measurement using singleshot off-axis Fizeau interferometer, J. Opt. Soc. Korea, 14, pp. 409-414. Abdelsalam, D.; Shaalan, M.; & Eloker, M. (2010). Surface microtopography measurement of

Abdelsalam, D.; Shaalan, M.; Eloker, M.; & Kim, D. (2010). Radius of curvature

Cuche, E. ; Marquet, P.; & Depeursinge, D. (2000). Spatial filtering for zero-order and twin

Cuche, E.; Bevilacque, F.; & Depeursinge, C. (1999). Digital holography for quantitative

Cuche, E.; Marquet, P.; & Depeursinge, C. (2000). Aperture apodization using cubic spline

David, F.; John, T.; Lopez, R.; & Stahl, H. (1993). Vector formulation for interferogram

Ghiglia, D. ; & Pritt, M. (1998). Two-dimensional Phase Unwrapping : Theory, Algororithm,

Holden, J. (1949). Multiple-beam interferometer-intensity distribution in the reflected

Kinosita, K. (1953). Numerical evaluation of the intensity curve of a multiple-beam Fizeau

Kühn, J.; Colomb, T.; Montfort, F.; Emery, Y.; Cuche, E.; & Depeursinge, C. (2007). Real-time

Kumar, U. ; Bhaduri, B. ; Kothiyal, M. ; & Mohan Krishna. (2009). Two wavelength micro-

Malacara, D. ; Servin, M.; & Malacara, Z. (2005). Interferogram analysis for optical testing,

Pan, F.; Xiao, W.; Liu, S. ; Wang, F., & Li, R. (2011). Coherent noise reduction in digital

Wang, JY.; & Silva, DE. (1980). Wave-front interpretation with Zernike polynomials.

United States of America : Taylor and Francis Group, pp. 384-385.

dual-wavelength digital holographic microscopy with a single hologram

interferometry for 3-D surface profiling, Optics and Lasers in Engineering, 47, pp.

holographic phase contrast microscopy by slightly shifting object, Optics Express,

reflection, Optics and Lasers in Engineering, Vol. 48, pp. 643-649. Born, M.; Wolf, E. (1980). Principles of Optics*,* Cambridge University Press, pp 459-490. Bünnagel, R. (1956). Einfaches verfahren zur topographischen darstellung einer optischen

phase contrast imaging, Optics. Letters., 24, pp. 291-293.

surface fitting, Applied Optics, 32, pp. 4738-43. Gabor, D. (1948). A new microscopic principle, Nature, 161, pp. 777-778.

and Lasers in Engineering, Vol. 48, pp. 543-547.

planflache. Opt Acta, 3, pp. 81-85.

communications, 182, pp. 59-69.

Gehrcke, E. (1906). Interferenzen, pp 39, Germany.

223-229.

19, pp. 3863-3869.

Applied Optics, 19, pp. 1510-8.

and Software, Wiley, New York , USA.

system, Proc. Phys. Soc., 62, pp. 405-417.

fringe, J. Phys. Soc. Japan, 8, pp. 219-225.

acquisition, Optics Express, 15, pp. 7231-42.

Steve, B. (2006). Hand book of CCD astronomy, Cambridge, UK. Tolansky, S. (1960). Surface microtopography, London : Longmans.

a standard flat surface by multiple-beam interference fringes at reflection, Optics

measurement of spherical smooth surfaces by multiple-beam interferometry in

image elimination in digital off-axis holography, Applied Optics, 39, pp. 4070-4075.

interpolation :application in digital holographic microscopy, Optics

**9. References** 

Differential Synthetic Aperture Radar Interferometry (DInSAR) represents nowadays a wellestablished remote sensing technique to generate spatially dense surface deformation maps of large areas on Earth (Massonnet & Feigl, 1999; Bürgmann et. al., 2000). To achieve this task, the phase difference (known as interferogram) between SAR data pairs related to temporally-separated observations is exploited, retrieving a measure of the ground displacement projection along the radar line of sight (LOS) with centimeter to millimeter accuracy (Gabriel et. al., 1989).

Altough DinSAR methodology was originally applied to analyze single deformation episodes ( Massonnet et. al., 1993; Goldstein et. a., 1993; Massonnet et. al., 1995; Peltzer et. al., 1995), it has recently applied with success to investigate the temporal evolution of the detected deformation phenomena through the generation of displacement time-series (Ferretti et. al., 2001; Berardino et. al., 2002) . In this case, the analysis is based on the computation of deformation time-series via the inversion of a properly chosen set of interferograms, produced from a sequence of temporally-separated SAR acquisitions relevant to the investigated area. In this context, two main categories of advanced DInSAR techniques for deformation time-series generation have been proposed in literature, often referred to as Persistent Scatterers (PS) (Ferretti et al., 2001; Hooper et. al., 2004), and Small Baseline (SB) (Berardino et. al., 2002; Mora et. al., 2003) techniques, respectively. The PS algorithms select all the interferometric data pairs with reference to a single common master image, without any constraint on the temporal and spatial separation (baseline) among the orbits. In this case, the analysis is carried out at the full resolution spatial scale, and is focused on the pixels containing a single dominant scatterer thus ensuring very limited temporal and spatial decorrelation phenomena (Zebker & Villasenor, 1992). Instead, for what concerns the SB techniques, the interferograms are generated by considering multiple master images, because this allows having interferometric data pairs with small temporal and spatial baselines. Accordingly, distributed targets can be also investigated, and the analysis may exploit both single-look and multi-look interferograms. Within the SB framework, a popular approach is the one referred to as Small BAseline Subset (SBAS), which was originally developed for analyzing multi-look interferograms (Berardino et. al., 2002) and subsequently adapted to the full resolution ones; in the latter case the low

Advanced Multitemporal Phase Unwrapping Techniques for DInSAR Analyses 59

wherein a and r are the azimuth and range pixel spacing, respectively and az0 and rg0 the corresponding azimuth and range first lines. The discrete counterpart of the partial derivatives are then computed by using the popular wrapped-differences-of-wrapped-

> <sup>ˆ</sup> <sup>ˆ</sup> <sup>ˆ</sup> *x y*

*xy x y xy x y xy*

, 1, , 1, ,

Accordingly, the phase gradient vector is estimated by wrapping possible phase differences

implicitly assuming that, in a properly sampled interferogram, the phase differences of

 ˆ ˆ *rot* 

*xy x y*

 

ˆ ,, , ,

*rxy xy x y x y*

Consequently, the integration of such phase gradients depends on the chosen integration path. This suggests us that a way to unwrap interferograms is to correct the phase gradients in such a way that the unwrapped phase can be recovered independently of the integration direction. Moreover, we can observe that the phase gradient estimate has the advantage that its errors are localized and come in integer multiples of so that its curl (hereafter referred to as residue field) can be profitably used to reconstruct the full phase terms. The residue

*x y x y x y x y* 

*x y y x*

, 1, , 1 ,

Its values are either zero (no residues) or +/- (positive or negative residue, respectively). Therefore, the goal of the phase unwrapping procedure is to eliminate potential integration paths enclosing unequal numbers of positive and negative residues. Residue derives from two sources in the radar measurements. The former is related to the true discontinuities in the data: the fringe spacing may be so fine on certain topographic slopes or, from large interobservation displacement, such to exceed the Nyquist criterion of half-cycle spacing. The latter is the noise present in the data set, whatever from thermal and other noise sources or from decorrelation due to baseline length and temporal changes in the scene. However, residues from whatever source require compensation in the phase unwrapping procedure.

, , , ,

 

*xy xy xy xy xy*

, ,1 , ,1 ,

 

, ,

, ,,

 

, , ,

 ,

 

(4)

 

 

. Accordingly, the phase gradient

 

**<sup>ψ</sup> x y** (3)

interval by adding the correct multiples of and by

, ,

 

> 

 

 

> 

(6)

depends both on the noise level (i.e. the lower is the

0 (5)

 

 

interval. Moreover, the probability

). For this reason, the introduced

whose components with respect to the two spatial axis *x*ˆ and *y*ˆ are:

phases estimator , *arctg* sin cos 

> 

 ,

 

adjacent samples are likely to be restricted to the

field expression (see Figure 1) is the following:

 

coherence value the more likely the derivatives exceed

vector can be defined as follows:

*x*

*y*

greater than +/- in the

that the phase difference exceed

operator is not conservative, that is:

resolution results are properly extended to the full resolution scale, as discussed in (Lanari et. al., 2004).

Despite the differences among the various implementations of the PS and SB algorithms, which are outside the scope of this work, a common problem to be faced is the Phase Unwrapping (PhU) operation representing the retrieval process of the full phase signals from their (measured) modulo-2π restricted components, often referred to as "interferometric fringes". In the context of the advanced DInSAR techniques, the need of jointly analyzing sequences of multi-temporal interferograms has more recently promoted the development of new PhU approaches with improved performances with respect to those focused on processing single interferograms. In this chapter, following the description of the basic rationale of the most widely used PhU techniques, we present a short review of the most recent space-time PhU methodologies with a particular emphasis on the Extended Minimum Cost Flow (EMCF) and its recent improvements. Some experimental results obtained by applying these latter approaches will be shown to demonstrate the validity of the presented algorithms.

### **2. Phase unwrapping**

We start by introducing the key ideas at the base of most of the Phase Unwrapping (PhU) techniques proposed up to now in literature. As stated before, the phase unwrapping problem concerns the retrieval of the full interferometric phase (unwrapped ones) from its computed modulo-2π restricted components. Within DInSAR processing codes, PhU operation actually represents one the most critical task to be successfully accomplished. The problem in reconstructing the true phase from the measured one arises in the presence of aliasing errors mainly caused by the phase noise caused by low coherence and undersampling phenomena due to locally high fringe rates (related to the deformation signals associated to the interferograms).

To clarify the key aspects of the phase unwrapping problem, let us firstly focus on the problem to unwrap single interferograms. To this aim, let *a r <sup>z</sup>* , *<sup>g</sup>* be the interferometric phase measured in correspondence to the pixel of SAR coordinates , *<sup>z</sup> <sup>g</sup> a r* . The unwrapping problem can be stated as the searching of the multiple 2-integer that must be added to the wrapped phases to retrieve the full interferometric phase *a r <sup>z</sup>* , *<sup>g</sup>* . Indeed:

$$
\Psi\left(a\_{z'}r\_{\mathcal{g}}\right) = \phi\left(a\_{z'}r\_{\mathcal{g}}\right) + 2\pi H\left(a\_{z'}r\_{\mathcal{g}}\right) \tag{1}
$$

This represents an ill-posed problem admitting an infinite set of acceptable solutions that can benefit from the knowledge of external information. In the absence of such information, the unwrapped field is typically reconstructed by computing the phase spatial gradients and integrating them over a "consistent" path. Accordingly, PhU algorithms can be grouped in two main categories depending on the fact that the unwrapping solution is path-dependent or not. To efficiently solve the problem, it is convenient to represent it in terms of discrete coordinates *x*,*y* being:

$$\begin{aligned} a\_z &= a\_{z0} + \chi \Delta a \\ r\_\chi &= r\_{\chi^0} + y \Delta r \end{aligned} \tag{2}$$

resolution results are properly extended to the full resolution scale, as discussed in (Lanari

Despite the differences among the various implementations of the PS and SB algorithms, which are outside the scope of this work, a common problem to be faced is the Phase Unwrapping (PhU) operation representing the retrieval process of the full phase signals from their (measured) modulo-2π restricted components, often referred to as "interferometric fringes". In the context of the advanced DInSAR techniques, the need of jointly analyzing sequences of multi-temporal interferograms has more recently promoted the development of new PhU approaches with improved performances with respect to those focused on processing single interferograms. In this chapter, following the description of the basic rationale of the most widely used PhU techniques, we present a short review of the most recent space-time PhU methodologies with a particular emphasis on the Extended Minimum Cost Flow (EMCF) and its recent improvements. Some experimental results obtained by applying these latter approaches will be shown to demonstrate the validity of

We start by introducing the key ideas at the base of most of the Phase Unwrapping (PhU) techniques proposed up to now in literature. As stated before, the phase unwrapping problem concerns the retrieval of the full interferometric phase (unwrapped ones) from its computed modulo-2π restricted components. Within DInSAR processing codes, PhU operation actually represents one the most critical task to be successfully accomplished. The problem in reconstructing the true phase from the measured one arises in the presence of aliasing errors mainly caused by the phase noise caused by low coherence and undersampling phenomena due to locally high fringe rates (related to the deformation

To clarify the key aspects of the phase unwrapping problem, let us firstly focus on the

phase measured in correspondence to the pixel of SAR coordinates , *<sup>z</sup> <sup>g</sup> a r* . The unwrapping problem can be stated as the searching of the multiple 2-integer that must be

This represents an ill-posed problem admitting an infinite set of acceptable solutions that can benefit from the knowledge of external information. In the absence of such information, the unwrapped field is typically reconstructed by computing the phase spatial gradients and integrating them over a "consistent" path. Accordingly, PhU algorithms can be grouped in two main categories depending on the fact that the unwrapping solution is path-dependent or not. To efficiently solve the problem, it is convenient to represent it in terms of discrete

> 0 0 *z z g g*

*a a xa r r yr* 

  *a r a r Ha r <sup>z</sup>* , ,2 , *<sup>g</sup> <sup>z</sup> <sup>g</sup> <sup>z</sup> <sup>g</sup>* (1)

*a r <sup>z</sup>* , *<sup>g</sup>* be the interferometric

*a r <sup>z</sup>* , *<sup>g</sup>* . Indeed:

(2)

et. al., 2004).

the presented algorithms.

**2. Phase unwrapping** 

coordinates *x*,*y* being:

signals associated to the interferograms).

problem to unwrap single interferograms. To this aim, let

added to the wrapped phases to retrieve the full interferometric phase

wherein a and r are the azimuth and range pixel spacing, respectively and az0 and rg0 the corresponding azimuth and range first lines. The discrete counterpart of the partial derivatives are then computed by using the popular wrapped-differences-of-wrappedphases estimator , *arctg* sin cos . Accordingly, the phase gradient vector can be defined as follows:

$$
\hat{\nabla}\Psi = \left< \Delta\_x \boldsymbol{\nu} \right>\_{-\pi,\pi} \hat{\mathbf{x}} + \left< \Delta\_y \boldsymbol{\nu} \right>\_{-\pi,\pi} \hat{\mathbf{y}} \tag{3}
$$

whose components with respect to the two spatial axis *x*ˆ and *y*ˆ are:

$$\begin{aligned} \left< \Delta\_x \boldsymbol{\nu} \left( \mathbf{x}, \boldsymbol{y} \right) \right>\_{-\pi, \pi} &= \left< \boldsymbol{\nu} \left( \mathbf{x} + \mathbf{1}, \boldsymbol{y} \right) - \boldsymbol{\nu} \left( \mathbf{x}, \boldsymbol{y} \right) \right>\_{-\pi, \pi} = \left< \boldsymbol{\phi} \left( \mathbf{x} + \mathbf{1}, \boldsymbol{y} \right) - \boldsymbol{\phi} \left( \mathbf{x}, \boldsymbol{y} \right) \right>\_{-\pi, \pi} \\ \left< \Delta\_y \boldsymbol{\nu} \left( \mathbf{x}, \boldsymbol{y} \right) \right>\_{-\pi, \pi} &= \left< \boldsymbol{\nu} \left( \mathbf{x}, \boldsymbol{y} + \mathbf{1} \right) - \boldsymbol{\nu} \left( \mathbf{x}, \boldsymbol{y} \right) \right>\_{-\pi, \pi} = \left< \boldsymbol{\phi} \left( \mathbf{x}, \boldsymbol{y} + \mathbf{1} \right) - \boldsymbol{\phi} \left( \mathbf{x}, \boldsymbol{y} \right) \right>\_{-\pi, \pi} \end{aligned} \tag{4}$$

Accordingly, the phase gradient vector is estimated by wrapping possible phase differences greater than +/- in the , interval by adding the correct multiples of and by implicitly assuming that, in a properly sampled interferogram, the phase differences of adjacent samples are likely to be restricted to the , interval. Moreover, the probability that the phase difference exceed depends both on the noise level (i.e. the lower is the coherence value the more likely the derivatives exceed ). For this reason, the introduced operator is not conservative, that is:

$$\operatorname{rot}\left(\vec{\nabla}\,\varphi\right) = \nabla \times \left(\vec{\nabla}\,\varphi\right) \neq 0\tag{5}$$

Consequently, the integration of such phase gradients depends on the chosen integration path. This suggests us that a way to unwrap interferograms is to correct the phase gradients in such a way that the unwrapped phase can be recovered independently of the integration direction. Moreover, we can observe that the phase gradient estimate has the advantage that its errors are localized and come in integer multiples of so that its curl (hereafter referred to as residue field) can be profitably used to reconstruct the full phase terms. The residue field expression (see Figure 1) is the following:

$$\begin{split} \boldsymbol{\nu}\left(\mathbf{x},\boldsymbol{y}\right) &= \nabla \times \Big(\hat{\nabla}\,\boldsymbol{\nu}\left(\mathbf{x},\boldsymbol{y}\right)\Big) = \Delta\_{\boldsymbol{x}}\Big[\left<\Delta\_{\boldsymbol{y}}\,\boldsymbol{\nu}\left(\mathbf{x},\boldsymbol{y}\right)\right>\_{-\pi,\pi}\Big] - \Delta\_{\boldsymbol{y}}\Big[\left<\Delta\_{\boldsymbol{x}}\,\boldsymbol{\nu}\left(\mathbf{x},\boldsymbol{y}\right)\right>\_{-\pi,\pi}\Big] = \\ &= \left<\Delta\_{\boldsymbol{x}}\boldsymbol{\nu}\left(\mathbf{x},\boldsymbol{y}\right)\right>\_{-\pi,\pi} + \left<\Delta\_{\boldsymbol{y}}\boldsymbol{\nu}\left(\mathbf{x}+1,\boldsymbol{y}\right)\right>\_{-\pi,\pi} - \left<\Delta\_{\boldsymbol{x}}\boldsymbol{\nu}\left(\mathbf{x},\boldsymbol{y}+1\right)\right>\_{-\pi,\pi} - \left<\Delta\_{\boldsymbol{y}}\boldsymbol{\nu}\left(\mathbf{x},\boldsymbol{y}\right)\right>\_{-\pi,\pi} \end{split} \tag{6}$$

Its values are either zero (no residues) or +/- (positive or negative residue, respectively). Therefore, the goal of the phase unwrapping procedure is to eliminate potential integration paths enclosing unequal numbers of positive and negative residues. Residue derives from two sources in the radar measurements. The former is related to the true discontinuities in the data: the fringe spacing may be so fine on certain topographic slopes or, from large interobservation displacement, such to exceed the Nyquist criterion of half-cycle spacing. The latter is the noise present in the data set, whatever from thermal and other noise sources or from decorrelation due to baseline length and temporal changes in the scene. However, residues from whatever source require compensation in the phase unwrapping procedure.

Advanced Multitemporal Phase Unwrapping Techniques for DInSAR Analyses 61

tree is zero. Networks on small trees are used to prevent any single branch from becoming too long and isolating large sub-areas from the rest of the image. A consequence of the indiscriminate branch growth until charge neutrality is achieved from all trees, and all residues accounted for, is that in residue-rich regions the tree growth is so dense that the region is isolated from the remainder of the image and no unwrapped phase estimate can be obtained. By concluding, this conservative approach nearly eliminates mistakes but, at the

The second major class of Phase Unwrapping algorithms, commonly in use today, was presented by (Ghiglia & Romero, 1994), who applied a mathematical formalism first developed by Hunt [69] to the radar interferometry phase unwrapping problem. Hunt developed a matrix formulation suitable for general phase reconstruction problems; Ghiglia found that a discrete cosine transform technique permits accurate and efficient least-squares inversion, even for the very large matrices encountered in the radar interferometry special case. In particular, the un-weighted LS method performs the following minimization

min ˆ , ,

Equation (7) represents a variational problem, whose Euler equation is the Poisson one:

**<sup>ψ</sup> ψ ψ**

where **ψ** is the unknown unwrapped phase field. We may stress that, with respect to the branch-cut approaches, in these cases the solution will be no longer congruent with the

> <sup>2</sup> , , ,, , , *x x y y xy xy x y x y*

under the Neumann boundary condition, which can be finally solved by using simple cosine

One major difference between the residue-cut and least-squares solutions is that in the residue-cut approach only integral numbers of cycles are added to the measurements to produce the result. Conversely, in the least-squares approach, any value may be added to ensure smoothness and continuity in the solution, thus the spatial error distribution may differ between the approaches, and the relative merits of each method must be determined

Least-squares methods are very computationally efficient when they make use of fast Fourier transform techniques but the resulting unwrapping is not very accurate, because they tend to spread the errors instead of concentrating them on a limited set of points. To overcome this problem a weighting of the wrapped phase can be useful. However, the proposed, weighted least squares algorithms are iterative and not as efficient as the un-weighted ones, and the result accuracy strongly depends on the used

 

*i j*

<sup>2</sup>

 

(8)

(7)

 

> 

*xy xy*

expense of providing an incomplete unwrapping result.

**2.2 Least-squares algorithms** 

original interferometric phase.

depending on the application.

weighting mask.

or Fourier transform filtering (Fornaro et. al., 1996).

problem

One of the major PhU algorithm exploits the fact that residues mark the endpoints of lines in the interferogram along which the true phase gradient exceed /sample, these lines are commonly referred to as "branch-cuts" or "ghost-lines". Most of the algorithm up to now proposed are based on a proper compensation of these residues and, among these, we will concentrate in the following on the residue-cut and the least-square approaches (Zebker & Lu, 1998).

Fig. 1. Representation of the residue field.

### **2.1 Branch-cut algorithms**

Within this class of algorithms we will point out our interest in particular on the so-called "residue-cut tree algorithm". The initial residue-cut phase unwrapping procedure proposed in (Golstein et. al., 1998) is implemented by first of all identifying the locations of all residues in an interferogram, and then connecting them with branch-cuts, so as to prevent the existence of integration paths that can encircle unbalanced numbers of positive and negative residues. The tree algorithm is a relatively conservative algorithm that tends to grow rather dense networks of trees in residue-rich regions. The algorithm initially connects closely spaced, oppositely charged, pairs of residues with cuts that prevent integration paths between them and, if all permitted integration paths enclose equal numbers of positive and negative residues, the consistency is assured. Progressively longer trees are permitted until all residues are connected to, at least one, other residue and until the next charge on each tree is zero. Networks on small trees are used to prevent any single branch from becoming too long and isolating large sub-areas from the rest of the image. A consequence of the indiscriminate branch growth until charge neutrality is achieved from all trees, and all residues accounted for, is that in residue-rich regions the tree growth is so dense that the region is isolated from the remainder of the image and no unwrapped phase estimate can be obtained. By concluding, this conservative approach nearly eliminates mistakes but, at the expense of providing an incomplete unwrapping result.

### **2.2 Least-squares algorithms**

60 Recent Interferometry Applications in Topography and Astronomy

One of the major PhU algorithm exploits the fact that residues mark the endpoints of lines in the interferogram along which the true phase gradient exceed /sample, these lines are commonly referred to as "branch-cuts" or "ghost-lines". Most of the algorithm up to now proposed are based on a proper compensation of these residues and, among these, we will concentrate in the following on the residue-cut and the least-square approaches (Zebker &

Within this class of algorithms we will point out our interest in particular on the so-called "residue-cut tree algorithm". The initial residue-cut phase unwrapping procedure proposed in (Golstein et. al., 1998) is implemented by first of all identifying the locations of all residues in an interferogram, and then connecting them with branch-cuts, so as to prevent the existence of integration paths that can encircle unbalanced numbers of positive and negative residues. The tree algorithm is a relatively conservative algorithm that tends to grow rather dense networks of trees in residue-rich regions. The algorithm initially connects closely spaced, oppositely charged, pairs of residues with cuts that prevent integration paths between them and, if all permitted integration paths enclose equal numbers of positive and negative residues, the consistency is assured. Progressively longer trees are permitted until all residues are connected to, at least one, other residue and until the next charge on each

Lu, 1998).

Fig. 1. Representation of the residue field.

**2.1 Branch-cut algorithms** 

The second major class of Phase Unwrapping algorithms, commonly in use today, was presented by (Ghiglia & Romero, 1994), who applied a mathematical formalism first developed by Hunt [69] to the radar interferometry phase unwrapping problem. Hunt developed a matrix formulation suitable for general phase reconstruction problems; Ghiglia found that a discrete cosine transform technique permits accurate and efficient least-squares inversion, even for the very large matrices encountered in the radar interferometry special case. In particular, the un-weighted LS method performs the following minimization problem

$$\min\_{\tilde{\Psi}} \left\{ \sum\_{i} \sum\_{j} \left| \nabla \tilde{\Psi}(\mathbf{x}\_{i} y) - \hat{\nabla} \Psi(\mathbf{x}\_{i} y) \right|^{2} \right\} \tag{7}$$

where **ψ** is the unknown unwrapped phase field. We may stress that, with respect to the branch-cut approaches, in these cases the solution will be no longer congruent with the original interferometric phase.

Equation (7) represents a variational problem, whose Euler equation is the Poisson one:

$$
\nabla^2 \tilde{\boldsymbol{\varphi}} \left( \mathbf{x}, \boldsymbol{y} \right) = \nabla \cdot \nabla \tilde{\boldsymbol{\varphi}} \left( \mathbf{x}, \boldsymbol{y} \right) = \Delta\_x \left( \left< \Delta\_x \boldsymbol{\varphi} \left( \mathbf{x}, \boldsymbol{y} \right) \right>\_{-\pi, \pi} \right) + \Delta\_y \left( \left< \Delta\_y \boldsymbol{\varphi} \left( \mathbf{x}, \boldsymbol{y} \right) \right>\_{-\pi, \pi} \right) \tag{8}
$$

under the Neumann boundary condition, which can be finally solved by using simple cosine or Fourier transform filtering (Fornaro et. al., 1996).

One major difference between the residue-cut and least-squares solutions is that in the residue-cut approach only integral numbers of cycles are added to the measurements to produce the result. Conversely, in the least-squares approach, any value may be added to ensure smoothness and continuity in the solution, thus the spatial error distribution may differ between the approaches, and the relative merits of each method must be determined depending on the application.

Least-squares methods are very computationally efficient when they make use of fast Fourier transform techniques but the resulting unwrapping is not very accurate, because they tend to spread the errors instead of concentrating them on a limited set of points. To overcome this problem a weighting of the wrapped phase can be useful. However, the proposed, weighted least squares algorithms are iterative and not as efficient as the un-weighted ones, and the result accuracy strongly depends on the used weighting mask.

Advanced Multitemporal Phase Unwrapping Techniques for DInSAR Analyses 63

 , , 1, , 1 , <sup>2</sup> *xy x y rxy K xy K x y K xy K xy*

that relates the integer unknowns to the measurable residues. At this stage, the phase unwrapping problem can be formulated as the searching of the K terms that satisfy the

min , , , ,

being c( ) the so-called cost functions allowing us to individuate areas where the location of branch-cuts is likely or unlikely. It is easy to verify that, whether the costs were chosen constant, the problem (12) would be equivalent to search for the minimum total cut-line length. Cost functions are essentially expressed as a function of the estimated local interferogram quality (by exploiting the spatial coherence, or the phase gradient density or other properly identified quality factors). The problem given in (12) is a non-linear minimization problem with integer variables, and, if the following change of variables is

> 

 

,, ,,

*c xyK xy c xyK xy*

*c xyK xy c xyK xy*

*xx yy*

(14)

(15)

, min 0, , , max 0, , , min 0, , , max 0, ,

*x x x x y y y y*

,, ,, *x y*

*x y xx yy*

It can be seen that the problem stated in (14) can be transformed to define a minimum cost flow problem on a network (see Figure 2), with the new variables representing the net flow running along the network arcs. Once the network has been solved, the solutions in terms of

> ,,, ,,, *xxx yyy*

 

*K x y K x y K x y K x y K x y K x y*

This section is focused on the presentation of the rationale at the base of the so-called Extended Minimum Cost Flow (EMCF) PhU technique that represents a way to efficiently incorporate temporal information in the phase unwrapping problem of sequence of multitemporal differential interferograms. Before discussing the characteristics of the proposed

 

it can be re-formulated via two different linear problem, as follows:

the 2-multiple integer functions will be expressed as follows:

**2.4 Extended Minimum Cost Flow (EMCF) algorithm** 

min

*K K*

*K xy K xy K xy K xy K xy K xy K xy K xy*

*c x y K x y c x y K x y* 

*xx yy K K x y x y*

constraints (11) and solve the following minimization problem:

*x y*

considered

(13)

(11)

(12)

Several other approaches, which can be found in the bibliography at the end of this work, have also been investigated but, in the following, we will address the capability of one of these, known as the minimum cost flow phase unwrapping technique, developed for the two-dimensional case by (Costantini, 1998), and belonging to the branch-cuts PhU algorithm class.

### **2.3 Minimum cost flow algorithms**

Within the branch-cuts phase unwrapping algorithms an easy and fast algorithm is based on a solution of an equivalent minimum cost flow network and will be here addressed.

Branch-cut methods are based on the integration of the estimated neighbouring pixel differences of the unwrapped phase along conservative paths, thus avoiding the regions where these estimated differences are inconsistent. The problem of building cuts delimiting these regions is very difficult and the resulting phase unwrapping algorithm is very computationally expensive. However, we may exploit the fact that the neighbouring pixel differences of the unwrapped phases are estimated with possibly an error that is an integer multiple of . This circumstance allows the formulation of the phase unwrapping problem as the one of minimizing the weighted deviations between the estimated and the unknown neighbouring pixel differences of the unwrapped phases with the constraint that the deviations must be integer multiple of . With this constraint, the unwrapping results will not depend critically on the weighting mask we used, and errors are prevented to spread.

Minimization problems with integer variables are usually computationally very complex. However, recognizing the network structure underlying the phase unwrapping problem, it makes possible to employ very efficient strategies for its solution. In fact, the problem can be equated to the one of finding the minimum cost flow on a network, for the solution of which there are very efficient algorithms. To explain its basic principles and clarify how it can be performed, we refer to the unknown, unwrapped phase field and we impose that the result is consistent, thus requiring the irrotational property of this field:

$$\begin{aligned} \nabla \times \nabla \,\boldsymbol{\nu}\left(\mathbf{x}, \boldsymbol{y}\right) &= \Delta\_{\mathbf{x}}\left(\Delta\_{\mathbf{y}}\boldsymbol{\nu}\left(\mathbf{x}, \boldsymbol{y}\right)\right) - \Delta\_{\mathbf{y}}\left(\Delta\_{\mathbf{x}}\boldsymbol{\nu}\left(\mathbf{x}, \boldsymbol{y}\right)\right) = \\ \nabla \cdot \boldsymbol{\Delta}\_{\mathbf{x}}\boldsymbol{\nu}\left(\mathbf{x}, \boldsymbol{y}\right) + \Delta\_{\mathbf{y}}\boldsymbol{\nu}\left(\mathbf{x} + \mathbf{1}, \boldsymbol{y}\right) - \Delta\_{\mathbf{x}}\boldsymbol{\nu}\left(\mathbf{x}, \boldsymbol{y} + \mathbf{1}\right) - \Delta\_{\mathbf{y}}\boldsymbol{\nu}\left(\mathbf{x}, \boldsymbol{y}\right) = 0 \end{aligned} \tag{9}$$

Obviously, we can also express each term of (9) with respect to the wrapped phase derivates by introducing, for each phase term, a corresponding, unknown -multiple term, as follows:

$$\begin{aligned} \Delta\_x \boldsymbol{\nu} \left( \mathbf{x}, \boldsymbol{y} \right) &= \left< \Delta\_x \boldsymbol{\nu} \left( \mathbf{x}, \boldsymbol{y} \right) \right>\_{-\pi, \pi} + 2\pi K\_x \left( \mathbf{x}, \boldsymbol{y} \right) \\ \Delta\_y \boldsymbol{\nu} \left( \mathbf{x} + 1, \boldsymbol{y} \right) &= \left< \Delta\_y \boldsymbol{\nu} \left( \mathbf{x} + 1, \boldsymbol{y} \right) \right>\_{-\pi, \pi} + 2\pi K\_y \left( \mathbf{x} + 1, \boldsymbol{y} \right) \\ \Delta\_x \boldsymbol{\nu} \left( \mathbf{x}, \boldsymbol{y} + 1 \right) &= \left< \Delta\_x \boldsymbol{\nu} \left( \mathbf{x}, \boldsymbol{y} + 1 \right) \right>\_{-\pi, \pi} + 2\pi K\_x \left( \mathbf{x}, \boldsymbol{y} + 1 \right) \\ \Delta\_y \boldsymbol{\nu} \left( \mathbf{x}, \boldsymbol{y} \right) &= \left< \Delta\_y \boldsymbol{\nu} \left( \mathbf{x}, \boldsymbol{y} \right) \right>\_{-\pi, \pi} + 2\pi K\_y \left( \mathbf{x}, \boldsymbol{y} \right) \end{aligned} \tag{10}$$

These relations eventually lead to the following equation

Several other approaches, which can be found in the bibliography at the end of this work, have also been investigated but, in the following, we will address the capability of one of these, known as the minimum cost flow phase unwrapping technique, developed for the two-dimensional case by (Costantini, 1998), and belonging to the branch-cuts PhU algorithm

Within the branch-cuts phase unwrapping algorithms an easy and fast algorithm is based on

Branch-cut methods are based on the integration of the estimated neighbouring pixel differences of the unwrapped phase along conservative paths, thus avoiding the regions where these estimated differences are inconsistent. The problem of building cuts delimiting these regions is very difficult and the resulting phase unwrapping algorithm is very computationally expensive. However, we may exploit the fact that the neighbouring pixel differences of the unwrapped phases are estimated with possibly an error that is an integer multiple of . This circumstance allows the formulation of the phase unwrapping problem as the one of minimizing the weighted deviations between the estimated and the unknown neighbouring pixel differences of the unwrapped phases with the constraint that the deviations must be integer multiple of . With this constraint, the unwrapping results will not depend critically on the weighting mask we used, and errors are

Minimization problems with integer variables are usually computationally very complex. However, recognizing the network structure underlying the phase unwrapping problem, it makes possible to employ very efficient strategies for its solution. In fact, the problem can be equated to the one of finding the minimum cost flow on a network, for the solution of which there are very efficient algorithms. To explain its basic principles and clarify how it can be performed, we refer to the unknown, unwrapped phase field and we impose that the result

*x y y x*

,, ,

*xy xy x y*

*xy x y*

Obviously, we can also express each term of (9) with respect to the wrapped phase derivates by introducing, for each phase term, a corresponding, unknown -multiple term, as

 

*xy xy K xy*

, , 2,

*xx x*

*yy y*

These relations eventually lead to the following equation

*xx x*

*yy y*

,

 

,

 

, ,

 

> 

 

 

 

 

 

*x y x y Kx y*

1, 1, 2 1,

*x y x y K xy*

,1 ,1 2 ,1

*xy xy K xy*

, , 2,

*xy x y xy xy*

, 1, , 1 , 0

 

> 

(9)

(10)

is consistent, thus requiring the irrotational property of this field:

a solution of an equivalent minimum cost flow network and will be here addressed.

class.

**2.3 Minimum cost flow algorithms** 

prevented to spread.

follows:

$$K\_x\left(\mathbf{x},\mathbf{y}\right) + K\_y\left(\mathbf{x}+\mathbf{1},\mathbf{y}\right) - K\_x\left(\mathbf{x},\mathbf{y}+\mathbf{1}\right) - K\_y\left(\mathbf{x},\mathbf{y}\right) = -\frac{r\left(\mathbf{x},\mathbf{y}\right)}{2\pi} \tag{11}$$

that relates the integer unknowns to the measurable residues. At this stage, the phase unwrapping problem can be formulated as the searching of the K terms that satisfy the constraints (11) and solve the following minimization problem:

$$\min\_{\mathbf{K}\_{\mathbf{x}}\mathbf{K}\_{\mathbf{y}}} \left| \sum\_{\mathbf{x}} \sum\_{\mathbf{y}} c\_{\mathbf{x}}(\mathbf{x}, \mathbf{y}) \Big| \mathbf{K}\_{\mathbf{x}}(\mathbf{x}, \mathbf{y}) \right| + \sum\_{\mathbf{x}} \sum\_{\mathbf{y}} c\_{\mathbf{y}}(\mathbf{x}, \mathbf{y}) \Big| \mathbf{K}\_{\mathbf{y}}(\mathbf{x}, \mathbf{y}) \Big| \right| \tag{12}$$

being c( ) the so-called cost functions allowing us to individuate areas where the location of branch-cuts is likely or unlikely. It is easy to verify that, whether the costs were chosen constant, the problem (12) would be equivalent to search for the minimum total cut-line length. Cost functions are essentially expressed as a function of the estimated local interferogram quality (by exploiting the spatial coherence, or the phase gradient density or other properly identified quality factors). The problem given in (12) is a non-linear minimization problem with integer variables, and, if the following change of variables is considered

$$\begin{aligned} K\_x^-(\mathbf{x}, \boldsymbol{y}) &= \min \left[ \mathbf{0}, K\_x(\mathbf{x}, \boldsymbol{y}) \right] \\ K\_x^+(\mathbf{x}, \boldsymbol{y}) &= \max \left[ \mathbf{0}, K\_x(\mathbf{x}, \boldsymbol{y}) \right] \\ K\_y^-(\mathbf{x}, \boldsymbol{y}) &= \min \left[ \mathbf{0}, K\_y(\mathbf{x}, \boldsymbol{y}) \right] \\ K\_y^+(\mathbf{x}, \boldsymbol{y}) &= \max \left[ \mathbf{0}, K\_y(\mathbf{x}, \boldsymbol{y}) \right] \end{aligned} \tag{13}$$

it can be re-formulated via two different linear problem, as follows:

$$\min\_{\mathbf{K}\_{\boldsymbol{x}},\mathbf{K}\_{\boldsymbol{y}}} \left\{ \sum\_{\mathbf{x}} \sum\_{\mathbf{y}} \frac{\left[c\_{\mathbf{x}}(\boldsymbol{x},\boldsymbol{y})\mathbf{K}\_{\boldsymbol{x}}^{+}(\boldsymbol{x},\boldsymbol{y}) + c\_{\mathbf{y}}(\boldsymbol{x},\boldsymbol{y})\mathbf{K}\_{\boldsymbol{y}}^{+}(\boldsymbol{x},\boldsymbol{y})\right]}{\left[c\_{\mathbf{x}}(\boldsymbol{x},\boldsymbol{y})\mathbf{K}\_{\boldsymbol{x}}^{-}(\boldsymbol{x},\boldsymbol{y}) + c\_{\mathbf{y}}(\boldsymbol{x},\boldsymbol{y})\mathbf{K}\_{\boldsymbol{y}}^{-}(\boldsymbol{x},\boldsymbol{y})\right]} \right\} \tag{14}$$

It can be seen that the problem stated in (14) can be transformed to define a minimum cost flow problem on a network (see Figure 2), with the new variables representing the net flow running along the network arcs. Once the network has been solved, the solutions in terms of the 2-multiple integer functions will be expressed as follows:

$$\begin{aligned} K\_{\mathbf{x}}\left(\mathbf{x},\boldsymbol{y}\right) &= K\_{\mathbf{x}}^{+}\left(\mathbf{x},\boldsymbol{y}\right) - K\_{\mathbf{x}}^{-}\left(\mathbf{x},\boldsymbol{y}\right) \\ K\_{\mathbf{y}}\left(\mathbf{x},\boldsymbol{y}\right) &= K\_{\mathbf{y}}^{+}\left(\mathbf{x},\boldsymbol{y}\right) - K\_{\mathbf{y}}^{-}\left(\mathbf{x},\boldsymbol{y}\right) \end{aligned} \tag{15}$$

### **2.4 Extended Minimum Cost Flow (EMCF) algorithm**

This section is focused on the presentation of the rationale at the base of the so-called Extended Minimum Cost Flow (EMCF) PhU technique that represents a way to efficiently incorporate temporal information in the phase unwrapping problem of sequence of multitemporal differential interferograms. Before discussing the characteristics of the proposed

Advanced Multitemporal Phase Unwrapping Techniques for DInSAR Analyses 65

Fig. 3. SAR data representation in the temporal/perpendicular baseline plane for the ERS SAR data analyzed in the following experiments. (A) SAR image distribution. (B) Delaunay triangulation. (C) Triangulation after removal of triangles with sides characterized by spatial and/or temporal baseline values exceeding the selected thresholds (corresponding in our

The proposed unwrapping solution, in addition to the spatial characteristics of each DInSAR interferogram, also exploits the temporal relationships among a properly selected interferogram sequence, thus allowing the performance improvement of the MCF technique.

experiments to 300 m and 1500 days, respectively).

approach, some considerations about the MCF technique are first in order. The original minimum cost flow algorithm, thought to be applied to a regular spatial grid, was improved in order to deal with sparse data (Costantini, 1998). In particular, the grid of the investigated samples is typically chosen to be relevant to the coherent pixels present into the DInSAR interferograms, while the Delaunay triangulation is used to define the neighbouring points and the elementary cycles in the set of the identified coherent sparse pixels. The possibility to extend this approach to the three-dimensional case, in order to simultaneously unwrap an interferometric sequence, has already been investigated in (Costantini et. al., 2002). However, in this case, the problem cannot be formulated in terms of network minimum cost flow procedures. Accordingly, no computationally efficient codes are available and, therefore, the overall unwrapping process can be extremely time-consuming, particularly for long interferogram sequences.

Fig. 2. A pictorial representation of the equivalent network to be solved for retrieving the integer multiples needed to reconstruct a conservative phase field. The network nodes are associated to each residue-cut and the bidirectional arcs are related to the phase differences arcs.

approach, some considerations about the MCF technique are first in order. The original minimum cost flow algorithm, thought to be applied to a regular spatial grid, was improved in order to deal with sparse data (Costantini, 1998). In particular, the grid of the investigated samples is typically chosen to be relevant to the coherent pixels present into the DInSAR interferograms, while the Delaunay triangulation is used to define the neighbouring points and the elementary cycles in the set of the identified coherent sparse pixels. The possibility to extend this approach to the three-dimensional case, in order to simultaneously unwrap an interferometric sequence, has already been investigated in (Costantini et. al., 2002). However, in this case, the problem cannot be formulated in terms of network minimum cost flow procedures. Accordingly, no computationally efficient codes are available and, therefore, the overall unwrapping process can be extremely time-consuming, particularly

Fig. 2. A pictorial representation of the equivalent network to be solved for retrieving the integer multiples needed to reconstruct a conservative phase field. The network nodes are associated to each residue-cut and the bidirectional arcs are related to the phase differences

for long interferogram sequences.

arcs.

Fig. 3. SAR data representation in the temporal/perpendicular baseline plane for the ERS SAR data analyzed in the following experiments. (A) SAR image distribution. (B) Delaunay triangulation. (C) Triangulation after removal of triangles with sides characterized by spatial and/or temporal baseline values exceeding the selected thresholds (corresponding in our experiments to 300 m and 1500 days, respectively).

The proposed unwrapping solution, in addition to the spatial characteristics of each DInSAR interferogram, also exploits the temporal relationships among a properly selected interferogram sequence, thus allowing the performance improvement of the MCF technique.

Advanced Multitemporal Phase Unwrapping Techniques for DInSAR Analyses 67

phases as a starting point for the subsequent spatial unwrapping operation, performed on each single interferogram; again the standard MCF network programming technique is applied and, in this case, we use the "costs" of the previous solutions to set the weights associated with the single arcs involved in the spatial unwrapping. Overall, the basic rationale of the procedure is quite simple: to exploit the temporal relationships among the computed multi-temporal interferograms to bootstrap the subsequent spatial unwrapping operation. We remark that the presented approach is focused on multilook interferograms. Moreover, the possibility to apply the MCF network programming algorithm to both the temporal and spatial unwrapping steps leads to a computationally efficient procedure. Let us start our analysis by investigating the generation process of the interferograms needed for the algorithm implementation discussed in the following sections. Accordingly, we consider a set of N+1 independent SAR images of the same area, which are co-registered with respect to a reference one, with respect to which we may compute the temporal and spatial (perpendicular) baseline components. Accordingly, each SAR image can be represented by a point in the temporal/perpendicular baseline plane (see Figure 3A) where we may also compute a Delaunay triangulation (see Figure 3B). Despite the constraints on the maximum allowed baseline extensions, we remark that the obtained triangulation may involve arcs relevant to data pairs whose baselines exceed the assumed maxima, thus potentially leading to generate interferograms strongly decorrelated. To avoid these effects without loosing our triangulation representation, we may remove all the triangles involving arcs with too large baselines, as shown in Figure 3C. Equivalently, we may also remove the triangles corresponding to interferograms including data pairs with large Doppler centroid differences; this is often the case for interferograms involving ERS data acquired after 2000, i.e., following the gyroscopes failure events (Miranda et. al., 2008). Notice that this triangles removal step may lead to discarding some acquisitions and/or to the generation of more than one independent subset of triangles, i.e., to a data representation consistent with the one described in the D-InSAR Small Baseline Subset (SBAS) procedure. Therefore, the compatibility between this data organization and the one exploited in the SBAS technique is clearly evident. Following the identification of the final triangulation, the computation of the DInSAR interferogram sequence is performed. At this stage, we can generate the "mask" of the pixels in the spatial plane that are considered coherent within the generated sequence. This mask can be obtained, for instance, by considering those pixels with an estimated coherence value greater than a selected threshold, which are common to a part or even to the entire interferogram sequence. Accordingly, as a final step, we compute a second Delaunay triangulation which involves the arcs connecting the neighbouring pixels of the

The presented PhU procedure is based on a two-step processing approach that benefits of the information available from both the grids. In particular, the key idea is to carry out first, for each arc connecting neighbouring pixels, a "temporal" unwrapping operation which implies the basic MCF approach. The second step relies on the use of these results as a starting point for the "spatial" unwrapping performed on each single interferogram. The key issues of these two processing steps are described in the following analysis which is focused on the use of multilook interferograms. The temporal PhU method benefits from the relationships existing between the phase differences of pixel pairs relevant to the wrapped and unwrapped signals. In particular, if we consider the spatial arc connecting the generic A and B pixel pair in the spatial plane, the unknown, unwrapped phase difference can be

computed mask.

expressed as follows

Fig. 4. Delaunay triangulation in the azimuth/range plane involving the set of spatially coherent pixels (in black).

The starting point of our approach, mostly oriented to deformation time-series generation, is the computation of two Delaunay triangulations. The former is relevant to the SAR data acquisitions distribution in the so-called "Temporal/Perpendicular baseline" plane and allows us to identify the DInSAR interferograms sequence to be computed; the latter is carried out in the Azimuth/Range plane and involves the coherent pixels common to the interferograms within the sequence. The unwrapping operation of the overall data set is then performed via a two-step processing procedure: first of all, we identify all the segments of the computed "spatial" triangles and, for each of these, we carry out a "temporal" PhU step by applying the MCF technique to the grid relevant to the temporal/perpendicular baseline plane. The second step uses, for each arc, the previously computed unwrapped

Fig. 4. Delaunay triangulation in the azimuth/range plane involving the set of spatially

The starting point of our approach, mostly oriented to deformation time-series generation, is the computation of two Delaunay triangulations. The former is relevant to the SAR data acquisitions distribution in the so-called "Temporal/Perpendicular baseline" plane and allows us to identify the DInSAR interferograms sequence to be computed; the latter is carried out in the Azimuth/Range plane and involves the coherent pixels common to the interferograms within the sequence. The unwrapping operation of the overall data set is then performed via a two-step processing procedure: first of all, we identify all the segments of the computed "spatial" triangles and, for each of these, we carry out a "temporal" PhU step by applying the MCF technique to the grid relevant to the temporal/perpendicular baseline plane. The second step uses, for each arc, the previously computed unwrapped

coherent pixels (in black).

phases as a starting point for the subsequent spatial unwrapping operation, performed on each single interferogram; again the standard MCF network programming technique is applied and, in this case, we use the "costs" of the previous solutions to set the weights associated with the single arcs involved in the spatial unwrapping. Overall, the basic rationale of the procedure is quite simple: to exploit the temporal relationships among the computed multi-temporal interferograms to bootstrap the subsequent spatial unwrapping operation. We remark that the presented approach is focused on multilook interferograms. Moreover, the possibility to apply the MCF network programming algorithm to both the temporal and spatial unwrapping steps leads to a computationally efficient procedure. Let us start our analysis by investigating the generation process of the interferograms needed for the algorithm implementation discussed in the following sections. Accordingly, we consider a set of N+1 independent SAR images of the same area, which are co-registered with respect to a reference one, with respect to which we may compute the temporal and spatial (perpendicular) baseline components. Accordingly, each SAR image can be represented by a point in the temporal/perpendicular baseline plane (see Figure 3A) where we may also compute a Delaunay triangulation (see Figure 3B). Despite the constraints on the maximum allowed baseline extensions, we remark that the obtained triangulation may involve arcs relevant to data pairs whose baselines exceed the assumed maxima, thus potentially leading to generate interferograms strongly decorrelated. To avoid these effects without loosing our triangulation representation, we may remove all the triangles involving arcs with too large baselines, as shown in Figure 3C. Equivalently, we may also remove the triangles corresponding to interferograms including data pairs with large Doppler centroid differences; this is often the case for interferograms involving ERS data acquired after 2000, i.e., following the gyroscopes failure events (Miranda et. al., 2008). Notice that this triangles removal step may lead to discarding some acquisitions and/or to the generation of more than one independent subset of triangles, i.e., to a data representation consistent with the one described in the D-InSAR Small Baseline Subset (SBAS) procedure. Therefore, the compatibility between this data organization and the one exploited in the SBAS technique is clearly evident. Following the identification of the final triangulation, the computation of the DInSAR interferogram sequence is performed. At this stage, we can generate the "mask" of the pixels in the spatial plane that are considered coherent within the generated sequence. This mask can be obtained, for instance, by considering those pixels with an estimated coherence value greater than a selected threshold, which are common to a part or even to the entire interferogram sequence. Accordingly, as a final step, we compute a second Delaunay triangulation which involves the arcs connecting the neighbouring pixels of the computed mask.

The presented PhU procedure is based on a two-step processing approach that benefits of the information available from both the grids. In particular, the key idea is to carry out first, for each arc connecting neighbouring pixels, a "temporal" unwrapping operation which implies the basic MCF approach. The second step relies on the use of these results as a starting point for the "spatial" unwrapping performed on each single interferogram. The key issues of these two processing steps are described in the following analysis which is focused on the use of multilook interferograms. The temporal PhU method benefits from the relationships existing between the phase differences of pixel pairs relevant to the wrapped and unwrapped signals. In particular, if we consider the spatial arc connecting the generic A and B pixel pair in the spatial plane, the unknown, unwrapped phase difference can be expressed as follows

Advanced Multitemporal Phase Unwrapping Techniques for DInSAR Analyses 69

Fig. 5. EMCF-based region-growing algorithm. (a) Pictorial representation of the *p-th* iteration of the algorithm: the phase values at the candidate pixel (located at the center) are predicted by the phases of the eight-seed-pixels (gray boxes) located in its neighborhood. (b)

Diagram block of the algorithm.

$$\begin{split} \Delta \tilde{\boldsymbol{\nu}}\_{AB} &= \tilde{\boldsymbol{\nu}}\left(\boldsymbol{A}\right) - \tilde{\boldsymbol{\nu}}\left(\boldsymbol{B}\right) = \left\langle \tilde{\boldsymbol{\rho}}\left(\boldsymbol{A}\right) - \tilde{\boldsymbol{\sigma}}\left(\boldsymbol{B}\right) \right\rangle\_{-\pi,\pi} + 2\pi H\_{AB} = \\ &= \Delta \tilde{\boldsymbol{\rho}}\_{AB} + 2\pi H\_{AB} \end{split} \tag{16}$$

Therefore, we may compute a Delaunay triangulation shown in Figure 4. In order to define a set of elementary cycles relevant to the coherent pixels only, and we may impose the irrotational property for the phase gradient, in a discrete space. At this stage, the PhU problem can be solved, as done for the original case, in a very efficient way by recognizing that a network structure underlines it and, by searching for the relevant network minimum cost flow solution, as fully described in (Pepe & Lanari, 2006). The unwrapped DInSAR phase differences computed via the previous temporal PhU step are then finally used as a starting point for the spatial unwrapping operation. This second step is carried out on the single interferograms through the application of the basic MCF unwrapping technique.

### **3. EMCF-based region growing technique**

A different class of PhU algorithms is based on the general concept of the "region growing" allowing us to considerably increase the spatial density of the areas with correct unwrapped results. In particular, we present a space-time region growing (RG) technique whose core is represented by the EMCF algorithm described in section 2. Basic idea of such approaches is to retrieve the unwrapped phase of a set of investigated pixels that are in the neighbourhood of an already unwrapped region. In particular, the unwrapped phase of the analyzed pixels is computed by exploiting the phase values of neighbouring unwrapped points (Xu & Cumming, 1999). Notice that with respect to classic region growing approaches here we apply a method that allows to handle with non-linear deformation trends: this is beneficial to analyze areas characterized by strongly non-linear displacements behaviours. Some experimental results will be also presented to demonstrate the effectiveness of the proposed PhU methodology. In particular, the implemented RG-EMCF algorithm has been applied to sequences of multilook differential interferograms, relevant to the Central Nevada region, characterized by large non-linear deformation phenomena, as well as to a test site located in the Gardanne area (France), affected by strong decorrelation phenomena.

The method provides a prediction of the (unwrapped) phase values at a pixel in the neighborhood of a region of pixels that have already been unwrapped. Following the preliminary PhU operation, we identify two sets of pixels. The former consists of *Q* high quality points **S** *SS S* 01 1 , ,..., *<sup>Q</sup>* , hereinafter referred to as seed points, that are characterized by temporal coherence values greater than a given threshold (typically set to 0,7 for DInSAR analyses). Notice that the temporal coherence is a measure of the correctness of the obtained unwrapping results originally proposed in (Pepe & Lanari, 2006) that is not fully discussed here for the sake of brevity. The region-growing procedure starts by analyzing the candidate pixels in the proximity of the DInSAR reference point location, and proceeds along a path (for instance, a spiral) including all the previously selected candidate pixels.

For the sake of generalization, let us describe the *p-th* iteration of the algorithm, aimed at analyzing the multitemporal unwrapped phases at the generic candidate pixel *Cp* . This is done by exploiting the knowledge of the unwrapped phase values **Ψ***Spi* at the seed pixels *pi S* , located within a rectangular area centered around the location of the candidate pixel [see Fig. 5(a)]. Thereby, we obtain different predictions of the (unwrapped) phases at the

, <sup>2</sup>

*AB AB H*

 

   

(16)

*AB AB*

Therefore, we may compute a Delaunay triangulation shown in Figure 4. In order to define a set of elementary cycles relevant to the coherent pixels only, and we may impose the irrotational property for the phase gradient, in a discrete space. At this stage, the PhU problem can be solved, as done for the original case, in a very efficient way by recognizing that a network structure underlines it and, by searching for the relevant network minimum cost flow solution, as fully described in (Pepe & Lanari, 2006). The unwrapped DInSAR phase differences computed via the previous temporal PhU step are then finally used as a starting point for the spatial unwrapping operation. This second step is carried out on the single interferograms through the application of the basic MCF unwrapping technique.

A different class of PhU algorithms is based on the general concept of the "region growing" allowing us to considerably increase the spatial density of the areas with correct unwrapped results. In particular, we present a space-time region growing (RG) technique whose core is represented by the EMCF algorithm described in section 2. Basic idea of such approaches is to retrieve the unwrapped phase of a set of investigated pixels that are in the neighbourhood of an already unwrapped region. In particular, the unwrapped phase of the analyzed pixels is computed by exploiting the phase values of neighbouring unwrapped points (Xu & Cumming, 1999). Notice that with respect to classic region growing approaches here we apply a method that allows to handle with non-linear deformation trends: this is beneficial to analyze areas characterized by strongly non-linear displacements behaviours. Some experimental results will be also presented to demonstrate the effectiveness of the proposed PhU methodology. In particular, the implemented RG-EMCF algorithm has been applied to sequences of multilook differential interferograms, relevant to the Central Nevada region, characterized by large non-linear deformation phenomena, as well as to a test site located in the Gardanne area (France), affected by strong decorrelation phenomena. The method provides a prediction of the (unwrapped) phase values at a pixel in the neighborhood of a region of pixels that have already been unwrapped. Following the preliminary PhU operation, we identify two sets of pixels. The former consists of *Q* high quality points **S** *SS S* 01 1 , ,..., *<sup>Q</sup>* , hereinafter referred to as seed points, that are characterized by temporal coherence values greater than a given threshold (typically set to 0,7 for DInSAR analyses). Notice that the temporal coherence is a measure of the correctness of the obtained unwrapping results originally proposed in (Pepe & Lanari, 2006) that is not fully discussed here for the sake of brevity. The region-growing procedure starts by analyzing the candidate pixels in the proximity of the DInSAR reference point location, and proceeds along a path (for

 

2

instance, a spiral) including all the previously selected candidate pixels.

For the sake of generalization, let us describe the *p-th* iteration of the algorithm, aimed at analyzing the multitemporal unwrapped phases at the generic candidate pixel *Cp* . This is done by exploiting the knowledge of the unwrapped phase values **Ψ***Spi* at the seed pixels *pi S* , located within a rectangular area centered around the location of the candidate pixel [see Fig. 5(a)]. Thereby, we obtain different predictions of the (unwrapped) phases at the

 

**3. EMCF-based region growing technique** 

*AB AB*

  *H*

Fig. 5. EMCF-based region-growing algorithm. (a) Pictorial representation of the *p-th* iteration of the algorithm: the phase values at the candidate pixel (located at the center) are predicted by the phases of the eight-seed-pixels (gray boxes) located in its neighborhood. (b) Diagram block of the algorithm.

Advanced Multitemporal Phase Unwrapping Techniques for DInSAR Analyses 71

In this section, we provide some experimental results relevant to two different test site areas (Casu, 2009): a wide area in central Nevada (USA) extending for about 600 x 100 km and the

The Central Nevada SAR data set consists of 264 ERS-1/2 SAR data frames (track 442, frames: from 2781 to 2871), spanning the 1992-2000 time interval, and we generate 148 multilook DInSAR interferograms, with a spatial resolution of about 160 x 160 m. Note that SAR data set and interferogram distribution are the same exploited by (Casu et al., 2008). First of all, we select the set of Seed Points from which the unwrapped phase will be propagated. To do this, we just impose a threshold value to the temporal coherence achieved by applying the EMCF algorithm, in particular we impose a value of 0.8, obtaining

more than 850,000 Seed pixels. In Figure 6a it is shown the selected Seed Point mask.

In order to evaluate the phase unwrapping step performances, we compare the achieved RG-EMCF results with those obtained by using a conventional SBAS processing chain. In particular, Figure 6b and Figure 6c show the temporal coherence masks of the reliable pixels (obtained by imposing previous mentioned thresholds for the RG-EMCF algorithm and for the conventional one) relevant to the conventional and new Region-Growing algorithm, respectively. By a first qualitative analysis, it is clear that the RG-EMCF approach permits to obtain a larger number of points at a given coherence threshold. More systematically, we measured an increase of about 88% for the number of "grown" pixels, while the total image

Moreover, we also compute the temporal coherence histogram, shown in Figure 9, relevant to the common "grown" pixels only. It is clear the coherence improvement achieved by the RG-EMCF results. As mean coherence value we obtain 0.90 for the RG-EMCF approach while 0.85 for the conventional SBAS results. Also in this case, it is clearly visible the

As final remark, we present in Figure 7 the mean deformation velocity maps of the two sets of unwrapped phases processed via the SBAS approach. Figure 8 clearly demonstrates the RG-EMCF capability to increment reliable pixels not only in non deforming areas, but also

In general, Region-Growing is very helpful to unwrap wide areas. Indeed, to manage large set of points can imply PhU infeasibility, due to the complexity of the problem to be solved as well as to the high computational burden. Therefore, a cascade of global and local PhU steps permits to correctly unwrap large areas, by first deliberately reducing the amount of investigated points (Seed Pixels) and, subsequently, reconsidering the removed ones

However, a Region-Growing approach can be applied also in highly decorrelated region, where the low coherence implies to deal with a reduced amount of Seed Pixels. Therefore, a

Accordingly, we applied the proposed RG-EMCF PhU approach to the Gardanne (France) region which is highly vegetated and where an open pit mine is present. Therefore, the zone

Region-Growing step can help increasing the unwrapped pixel spatial density.

**3.1 Real data results** 

increment is of about the 32%.

improvement of the new proposed algorithm.

in zone affected by strong and non linear deformation.

(Candidate Pixels) in a proper Region-Growing step.

highly vegetated region of Gardanne (France).

selected candidate pixel, namely *<sup>i</sup>* **<sup>Ψ</sup>** *Cp* , computed by integrating the phase differences along the paths connecting the candidate and the i-th seed pixel:

$$\mathbf{W}^{(i)}\left(\mathbf{C}\_{p}\right) = \mathbf{W}\left(\mathbf{S}\_{pi}\right) + \Delta\mathbf{W}\left(\mathbf{C}\_{p'}, \mathbf{S}\_{pi}\right) \quad \text{i} = 1, \dots, \mathbf{K}\_{p} \tag{17}$$

where the estimate of the (unwrapped) phase difference **Ψ***C S p p* , *<sup>i</sup>* over the given candidate/seed arc is retrieved by applying the temporal PhU strategy detailed in (Pepe & Lanari, 2006). However, since different integration paths lead to independent phase predictions at the candidate pixel, we carry out a weighted average of such individuals predictions:

$$\overline{\Psi}\left(\mathbf{C}\_{p}\right) = \left(\sum\_{i=1}^{\mathbb{K}\_{p}} w\_{i} \Psi^{(i)}\left(\mathbf{C}\_{p}\right)\right) \bigg/ \left(\sum\_{i=1}^{\mathbb{K}\_{p}} w\_{i}\right) \tag{18}$$

where the weights {w} are set taking into account the (achieved) temporal network minimum costs (Pepe & Lanari, 2006). The average prediction in (18) is eventually used to compute the (unwrapped) phase vector at the selected candidate pixel.

$$\mathbf{\bar{\Psi}}\left(\mathbf{C}\_{p}\right) = \overline{\mathbf{\bar{\Psi}}}\left(\mathbf{C}\_{p}\right) + \mathbf{r}\left(\mathbf{C}\_{p}\right) \quad \mathbf{i} = \mathbf{1}, \dots, \mathbf{K}\_{p} \tag{19}$$

where the , *C CC p pp* **<sup>r</sup> Φ Ψ** term on the right-hand side of (19) guarantees that unwrapped and wrapped phase vectors differ by -integer multiples.

A check whether the prediction in (19) is correct is finally performed: if this is the case, the selected candidate pixel is added to the set of seed points, otherwise it is discarded from the following analyses. To the purpose, two different reliability tests are exploited.

Test 1- We invert (19) by the SBAS strategy and measure the related temporal coherence value (Pepe et al. 2006). The test is passed when the measured temporal coherence is greater than a given threshold.

Test 2- It relies on the assumption the larger the values of the residual phase vector the larger the probability to retrieve incorrect phase estimates. Accordingly, we analyze the residual phase dispersion through:

$$\Lambda\left(\mathbf{C}\_p\right) = \frac{\left|\sum\_{k=0}^{M-1} \exp\left[j\pi\left(\mathbf{C}\_p\right)\right]\right|}{M} \tag{20}$$

This second test is passed when this additional coherence value is greater than a selected threshold.

The block diagram of the algorithm is sketched in Fig. 5b.

### **3.1 Real data results**

70 Recent Interferometry Applications in Topography and Astronomy

where the estimate of the (unwrapped) phase difference **Ψ***C S p p* , *<sup>i</sup>* over the given candidate/seed arc is retrieved by applying the temporal PhU strategy detailed in (Pepe & Lanari, 2006). However, since different integration paths lead to independent phase predictions at the candidate pixel, we carry out a weighted average of such individuals

1 1

**<sup>r</sup> Φ Ψ** term on the right-hand side of (19) guarantees that

where the weights {w} are set taking into account the (achieved) temporal network minimum costs (Pepe & Lanari, 2006). The average prediction in (18) is eventually used to

A check whether the prediction in (19) is correct is finally performed: if this is the case, the selected candidate pixel is added to the set of seed points, otherwise it is discarded from the

Test 1- We invert (19) by the SBAS strategy and measure the related temporal coherence value (Pepe et al. 2006). The test is passed when the measured temporal coherence is greater

Test 2- It relies on the assumption the larger the values of the residual phase vector the larger the probability to retrieve incorrect phase estimates. Accordingly, we analyze the

> <sup>1</sup> 0 exp

This second test is passed when this additional coherence value is greater than a selected

*p*

(20)

*jr C*

*M*

*M*

*k*

*p*

*C*

The block diagram of the algorithm is sketched in Fig. 5b.

*K K p p i p i p i i i C wC w* 

 

compute the (unwrapped) phase vector at the selected candidate pixel.

 

following analyses. To the purpose, two different reliability tests are exploited.

unwrapped and wrapped phase vectors differ by -integer multiples.

*<sup>i</sup>* **<sup>Ψ</sup>** *Cp* , computed by integrating the phase differences

, 1,..., *<sup>i</sup>* **Ψ ΨΨ** *C S CS i K p p <sup>i</sup> p pi <sup>p</sup>* (17)

**<sup>Ψ</sup>** (18)

1,..., **Ψ Ψ** *C C Ci K p pp p* **<sup>r</sup>** (19)

selected candidate pixel, namely

where the , *C CC p pp*

than a given threshold.

threshold.

residual phase dispersion through:

predictions:

along the paths connecting the candidate and the i-th seed pixel:

In this section, we provide some experimental results relevant to two different test site areas (Casu, 2009): a wide area in central Nevada (USA) extending for about 600 x 100 km and the highly vegetated region of Gardanne (France).

The Central Nevada SAR data set consists of 264 ERS-1/2 SAR data frames (track 442, frames: from 2781 to 2871), spanning the 1992-2000 time interval, and we generate 148 multilook DInSAR interferograms, with a spatial resolution of about 160 x 160 m. Note that SAR data set and interferogram distribution are the same exploited by (Casu et al., 2008). First of all, we select the set of Seed Points from which the unwrapped phase will be propagated. To do this, we just impose a threshold value to the temporal coherence achieved by applying the EMCF algorithm, in particular we impose a value of 0.8, obtaining more than 850,000 Seed pixels. In Figure 6a it is shown the selected Seed Point mask.

In order to evaluate the phase unwrapping step performances, we compare the achieved RG-EMCF results with those obtained by using a conventional SBAS processing chain. In particular, Figure 6b and Figure 6c show the temporal coherence masks of the reliable pixels (obtained by imposing previous mentioned thresholds for the RG-EMCF algorithm and for the conventional one) relevant to the conventional and new Region-Growing algorithm, respectively. By a first qualitative analysis, it is clear that the RG-EMCF approach permits to obtain a larger number of points at a given coherence threshold. More systematically, we measured an increase of about 88% for the number of "grown" pixels, while the total image increment is of about the 32%.

Moreover, we also compute the temporal coherence histogram, shown in Figure 9, relevant to the common "grown" pixels only. It is clear the coherence improvement achieved by the RG-EMCF results. As mean coherence value we obtain 0.90 for the RG-EMCF approach while 0.85 for the conventional SBAS results. Also in this case, it is clearly visible the improvement of the new proposed algorithm.

As final remark, we present in Figure 7 the mean deformation velocity maps of the two sets of unwrapped phases processed via the SBAS approach. Figure 8 clearly demonstrates the RG-EMCF capability to increment reliable pixels not only in non deforming areas, but also in zone affected by strong and non linear deformation.

In general, Region-Growing is very helpful to unwrap wide areas. Indeed, to manage large set of points can imply PhU infeasibility, due to the complexity of the problem to be solved as well as to the high computational burden. Therefore, a cascade of global and local PhU steps permits to correctly unwrap large areas, by first deliberately reducing the amount of investigated points (Seed Pixels) and, subsequently, reconsidering the removed ones (Candidate Pixels) in a proper Region-Growing step.

However, a Region-Growing approach can be applied also in highly decorrelated region, where the low coherence implies to deal with a reduced amount of Seed Pixels. Therefore, a Region-Growing step can help increasing the unwrapped pixel spatial density.

Accordingly, we applied the proposed RG-EMCF PhU approach to the Gardanne (France) region which is highly vegetated and where an open pit mine is present. Therefore, the zone

Advanced Multitemporal Phase Unwrapping Techniques for DInSAR Analyses 73

Fig. 7. Mean deformation velocity maps, in SAR coordinates, relevant to Central Nevada area. a) Conventional SBAS results. b) RG-EMCF results. By courtesy of dr. Francesco Casu

(Casu, 2009).

Fig. 6. Temporal coherence masks relevant to Central Nevada area. a) Seed Point mask (~850,000 pixels). b) RG-EMCF results (~1,770,000 pixels). c) Conventional SBAS results (~1,340,000 pixels). Grown points are highlighted in red while in white are represented the Seed Points of Figure 6a. By courtesy of dr. Francesco Casu (Casu, 2009).

Fig. 6. Temporal coherence masks relevant to Central Nevada area. a) Seed Point mask (~850,000 pixels). b) RG-EMCF results (~1,770,000 pixels). c) Conventional SBAS results (~1,340,000 pixels). Grown points are highlighted in red while in white are represented the

Seed Points of Figure 6a. By courtesy of dr. Francesco Casu (Casu, 2009).

Fig. 7. Mean deformation velocity maps, in SAR coordinates, relevant to Central Nevada area. a) Conventional SBAS results. b) RG-EMCF results. By courtesy of dr. Francesco Casu (Casu, 2009).

Advanced Multitemporal Phase Unwrapping Techniques for DInSAR Analyses 75

Fig. 9. Mean deformation velocity maps relevant to the Gardanne area (France). a)

In this section, we propose to exploit the highly efficient EMCF algorithm as the core procedure to unwrap sequences of large single-look DInSAR interferograms suitable for the

Conventional SBAS results. b) RG-EMCF results.

**4. On unwrapping large interferograms** 

Fig. 8. Temporal coherence histogram relevant to common "grown" pixels obtained by exploiting conventional (dashed line) and new (continuous line) Region-Growing algorithm. By courtesy of dr. Francesco Casu (Casu, 2009).

is characterized by very low temporal coherence and is affected by strongly non linear deformation, implying to be a very critical test area for multi-temporal DInSAR analysis.

In Figure 9, we present the mean deformation velocity maps computed by applying the SBAS approach to a set of 75 ERS-1/2 and 8 ENVISAT images acquired in the 1992-2004 time period and coupled in 243 interferograms, the latter being unwrapped via both the procedure implemented in the SBAS processing chain and the RG-EMCF algorithm. Note that, DInSAR data have been multilook, obtaining a final spatial ground resolution of about 80 by 80 meters.

We just remark that, simultaneous exploitation of ERS and ENVISAT data implies that no cross sensor interferograms can be generated, due to the different signal wavelength of the two sensors. Therefore, the Temporal/Perpendicular baseline network will result decomposed in several subsets (at least two), corresponding to the different sensor interferograms. Also in this case, looking at Figure 9, it is clearly visible the strong increase of the correctly unwrapped pixels. Indeed, "grown" pixel number pass from about 2,500 for the conventional case to more than 87,000 for the RG-EMCF one, with a global improvement (for the whole area, i. e., accounting also for the Seed Points) of about the 113%. It is worthy to remark that the RG approach available in the SBAS processing essentially does not produced any results (only 2,500 "grown" pixels starting from 73,000 Seeds): this is because the strong decorrelation affecting the area.

Fig. 8. Temporal coherence histogram relevant to common "grown" pixels obtained by exploiting conventional (dashed line) and new (continuous line) Region-Growing algorithm.

is characterized by very low temporal coherence and is affected by strongly non linear deformation, implying to be a very critical test area for multi-temporal DInSAR analysis.

In Figure 9, we present the mean deformation velocity maps computed by applying the SBAS approach to a set of 75 ERS-1/2 and 8 ENVISAT images acquired in the 1992-2004 time period and coupled in 243 interferograms, the latter being unwrapped via both the procedure implemented in the SBAS processing chain and the RG-EMCF algorithm. Note that, DInSAR data have been multilook, obtaining a final spatial ground resolution of about

We just remark that, simultaneous exploitation of ERS and ENVISAT data implies that no cross sensor interferograms can be generated, due to the different signal wavelength of the two sensors. Therefore, the Temporal/Perpendicular baseline network will result decomposed in several subsets (at least two), corresponding to the different sensor interferograms. Also in this case, looking at Figure 9, it is clearly visible the strong increase of the correctly unwrapped pixels. Indeed, "grown" pixel number pass from about 2,500 for the conventional case to more than 87,000 for the RG-EMCF one, with a global improvement (for the whole area, i. e., accounting also for the Seed Points) of about the 113%. It is worthy to remark that the RG approach available in the SBAS processing essentially does not produced any results (only 2,500 "grown" pixels starting from 73,000 Seeds): this is because

By courtesy of dr. Francesco Casu (Casu, 2009).

the strong decorrelation affecting the area.

80 by 80 meters.

Fig. 9. Mean deformation velocity maps relevant to the Gardanne area (France). a) Conventional SBAS results. b) RG-EMCF results.

### **4. On unwrapping large interferograms**

In this section, we propose to exploit the highly efficient EMCF algorithm as the core procedure to unwrap sequences of large single-look DInSAR interferograms suitable for the

Advanced Multitemporal Phase Unwrapping Techniques for DInSAR Analyses 77

Fig. 10. Examples of Triangulations. (a) Dataset of 96 points with highlighted 8 of them, labeled to as A, B, C, D, E, F, G, H, respectively. (b) The 8-points Delaunay triangulation drawn with red lines (c) Constrained Delaunay Triangulation generated from the set of points of (a), and by using as constraints the triangulation of (b). (d) Delaunay triangulation

carried out on each single interferogram via the basic MCF approach but, in order to preserve the unwrapped phases relevant to the primary network, the weights used for the spatial MCF minimization must be properly set. If we refer to the generic PQ arc, this is

> 100 1

drastically decreasing the overall computational burden.

*PQ Constrained*

*L PQ G w PQ G Cst*

*Constrained*

min min (21)

 

*Constrained*

where *GConstrained* is the set of constrained edges, and *Cst*min is the temporal minimum network cost relevant to the given spatial arc. Moreover, L is a very large integer number, and is a threshold value that is typically set not greater than 5% of the total number of interferograms. Based on (21), the flow into the constrained MCF network is automatically forced not to cross the constrained arcs and, as a consequence, the estimates of the primary network unwrapped phases are fully preserved. In other words, the presented approach allows us to effectively "propagate" the unwrapped solution from the primary network to the connected sub-networks, largely improving the phase unwrapping performances, and

*PQ G Cst*

computed from the 96 points in (a).

easily achieved by imposing:

generation of deformation time-series through the SBAS approach. This PhU approach allows us to directly apply the SBAS inversion to the unwrapped full resolution DInSAR phase sequences, with no need to pass through the analysis of the corresponding sequences of multi-look DInSAR interferograms (Lanari et. al., 2004). To properly solve this PhU problem, we suggest to apply an effective divide-and-conquer approach to the space-time phase unwrapping problem. The key idea is to divide the complex minimum cost flow network problems, implementing the whole PhU step, into that of simpler sub-networks, which are solved by applying the EMCF approach. More precisely, we start by identifying, and solving, a primary network that involves a selected set of very coherent pixels in our interferograms. The results of this primary network minimization, representing the backbone structure of the overall network, are subsequently used to constrain the solution of the remaining sub-networks, including the entire set of coherent pixels. To achieve this task, the second EMCF PhU step relies on the generation of a Constrained Delaunay Triangulation (CDT) (Chew, 1989), whose constrained edges are relevant to the set of successfully unwrapped pixels analyzed during the first PhU operation. We remark that our approach has some similarities with [6] where a two-scale strategy was also suggested to unwrap large interferograms; however our strategy is similar but inverted because in our scheme the primary PhU step is used to figure out a (global) PhU solution, and the secondary one to locally "propagate" the PhU solution in low coherent areas. The key idea of the proposed approach is to split our complex MCF network problem into that of simpler sub-networks. This solution can be efficiently implemented, as shown in the following, through two subsequent processing steps that are both carried out by using the EMCF technique.

Basically, the first PhU step is carried out on a set of very coherent pixels, used to compute a Delaunay triangulation in the Azimuth/Range plane, and the achieved PhU results are eventually exploited to successfully unwrap the remaining pixels. To achieve this task we solve a "constrained optimization problem" based on the computation of a Constrained Delaunay Triangulation (CDT) in the plane from the grid of the overall coherent pixels. To clarify this issue, let us provide some basic information about the CDT, which is a triangulation of a given set of vertices with the following properties: 1) a pre-specified set of non-crossing edges (referred to as constraints, or constrained edges) is included in the triangulation, and (2) the triangulation is as close as possible to a Delaunay one. As an example, in Figure 10, it is shown a simple CDT relevant to a set of 96 points, see Fig. 10(a); in this case the selected constraints are represented by the 14 edges of the Delaunay triangulation generated from an eight-points subset of the originally 96 ones, see Fig. 10(b); the computed CDT is shown in Fig. 10(c), whereas the corresponding Delaunay triangulation is presented in Fig. 10(d).

Similarly to the case of Figure 10, we compute a CDT for the spatial grid of the coherent pixels , whose constrained edges are the arcs of the previously identified primary network. Since the EMCF Phase Unwrapping algorithm can work with generic triangular irregular grids, not necessarily Delaunay triangulations, it can be also applied to the irregular spatial grid obtained via our CDT. However, in this case we must solve a constrained optimization problem because we want to preserve, for each interferogram, the unwrapped phase values already obtained by solving the primary network. Accordingly, we perform the second unwrapping step again through the EMCF approach but applying the temporal PhU step only to the unconstrained arcs of the generated CDT. Moreover, the spatial PhU step is

generation of deformation time-series through the SBAS approach. This PhU approach allows us to directly apply the SBAS inversion to the unwrapped full resolution DInSAR phase sequences, with no need to pass through the analysis of the corresponding sequences of multi-look DInSAR interferograms (Lanari et. al., 2004). To properly solve this PhU problem, we suggest to apply an effective divide-and-conquer approach to the space-time phase unwrapping problem. The key idea is to divide the complex minimum cost flow network problems, implementing the whole PhU step, into that of simpler sub-networks, which are solved by applying the EMCF approach. More precisely, we start by identifying, and solving, a primary network that involves a selected set of very coherent pixels in our interferograms. The results of this primary network minimization, representing the backbone structure of the overall network, are subsequently used to constrain the solution of the remaining sub-networks, including the entire set of coherent pixels. To achieve this task, the second EMCF PhU step relies on the generation of a Constrained Delaunay Triangulation (CDT) (Chew, 1989), whose constrained edges are relevant to the set of successfully unwrapped pixels analyzed during the first PhU operation. We remark that our approach has some similarities with [6] where a two-scale strategy was also suggested to unwrap large interferograms; however our strategy is similar but inverted because in our scheme the primary PhU step is used to figure out a (global) PhU solution, and the secondary one to locally "propagate" the PhU solution in low coherent areas. The key idea of the proposed approach is to split our complex MCF network problem into that of simpler sub-networks. This solution can be efficiently implemented, as shown in the following, through two subsequent processing steps that are both carried out by using the EMCF

Basically, the first PhU step is carried out on a set of very coherent pixels, used to compute a Delaunay triangulation in the Azimuth/Range plane, and the achieved PhU results are eventually exploited to successfully unwrap the remaining pixels. To achieve this task we solve a "constrained optimization problem" based on the computation of a Constrained Delaunay Triangulation (CDT) in the plane from the grid of the overall coherent pixels. To clarify this issue, let us provide some basic information about the CDT, which is a triangulation of a given set of vertices with the following properties: 1) a pre-specified set of non-crossing edges (referred to as constraints, or constrained edges) is included in the triangulation, and (2) the triangulation is as close as possible to a Delaunay one. As an example, in Figure 10, it is shown a simple CDT relevant to a set of 96 points, see Fig. 10(a); in this case the selected constraints are represented by the 14 edges of the Delaunay triangulation generated from an eight-points subset of the originally 96 ones, see Fig. 10(b); the computed CDT is shown in Fig. 10(c), whereas the corresponding Delaunay

Similarly to the case of Figure 10, we compute a CDT for the spatial grid of the coherent pixels , whose constrained edges are the arcs of the previously identified primary network. Since the EMCF Phase Unwrapping algorithm can work with generic triangular irregular grids, not necessarily Delaunay triangulations, it can be also applied to the irregular spatial grid obtained via our CDT. However, in this case we must solve a constrained optimization problem because we want to preserve, for each interferogram, the unwrapped phase values already obtained by solving the primary network. Accordingly, we perform the second unwrapping step again through the EMCF approach but applying the temporal PhU step only to the unconstrained arcs of the generated CDT. Moreover, the spatial PhU step is

technique.

triangulation is presented in Fig. 10(d).

Fig. 10. Examples of Triangulations. (a) Dataset of 96 points with highlighted 8 of them, labeled to as A, B, C, D, E, F, G, H, respectively. (b) The 8-points Delaunay triangulation drawn with red lines (c) Constrained Delaunay Triangulation generated from the set of points of (a), and by using as constraints the triangulation of (b). (d) Delaunay triangulation computed from the 96 points in (a).

carried out on each single interferogram via the basic MCF approach but, in order to preserve the unwrapped phases relevant to the primary network, the weights used for the spatial MCF minimization must be properly set. If we refer to the generic PQ arc, this is easily achieved by imposing:

$$\begin{array}{l} \text{for} \\ \text{for}\_{PQ} = \begin{cases} L & PQ \in \{\text{G}\_{\text{Construction}}\} \\ \text{100} & \left\{ PQ \notin \{\text{G}\_{\text{Construction}}\} \right\} \cap \left\{ \text{Cst}\_{\text{min}} < \rho \right\} \\ \text{1} & \left\{ PQ \notin \{\text{G}\_{\text{Construction}}\} \right\} \cap \left\{ \text{Cst}\_{\text{min}} > \rho \right\} \end{array} \tag{21}$$

where *GConstrained* is the set of constrained edges, and *Cst*min is the temporal minimum network cost relevant to the given spatial arc. Moreover, L is a very large integer number, and is a threshold value that is typically set not greater than 5% of the total number of interferograms. Based on (21), the flow into the constrained MCF network is automatically forced not to cross the constrained arcs and, as a consequence, the estimates of the primary network unwrapped phases are fully preserved. In other words, the presented approach allows us to effectively "propagate" the unwrapped solution from the primary network to the connected sub-networks, largely improving the phase unwrapping performances, and drastically decreasing the overall computational burden.

Advanced Multitemporal Phase Unwrapping Techniques for DInSAR Analyses 79

Fig. 12. ERS-1/2 DInSAR results. (a) Geocoded map of the mean deformation velocity of the investigated area (b) Zoomed view of the area of Napoli city and Campi Flegrei caldera (Italy) highlighted by the black box in (a). The plots show the DInSAR/leveling comparison of the deformation time-series corresponding to the pixels in (a) labeled as b (25A leveling

benchmark) and c (236 leveling benchmark), respectively.

Fig. 11. SAR data representation in the temporal/perpendicular baseline plane for the ERS-1/2 SAR data analyzed in the following experiments, relevant to the Napoli (Italy) bay area. (a) SAR image distribution. (b) Delaunay Triangulation after the removal of triangles with sides characterized by spatial and temporal baseline values, as well as doppler centroid differences, exceeding the selected thresholds (corresponding in our experiments to 400 m, 1500 days and 1000 Hz, respectively).

The proposed approach was validated by analyzing a dataset of 86 ERS-1/2 SAR images acquired on descending orbits (Track 36, Frame 2781) between June 8, 1992 and August 23, 2007, over the Napoli (Italy) Bay area (see Figure 11). From these data, we identified a set of 234 data pairs characterized by perpendicular baseline values smaller than 400 m, and a maximum time interval of 1500 days. Precise satellite orbital information and a 3 arcsec SRTM DEM of the area were used to generate a sequence of single-look DInSAR interferograms. The computed interferograms were then unwrapped by applying the EMCF approach. To achieve this task, we first identified the spatial grid of all the coherent pixels to be unwrapped, composed of about 530 000 pixels, where 50 000 of them are very coherent. From the latter pixels, we generated, in the spatial plane, a Delaunay triangulation. The arcs of this triangulation represent the constrained edges of the implemented CDT: this structure connects the overall set of coherent pixels , and is essential for the second EMCF PhU step, leading to the final estimate of the unwrapped interferograms on the chosen spatial grid.

Figure 12 shows a false color map of the detected mean deformation velocity, where only points with high data quality are included, superimposed on a multi-look SAR amplitude image of the area. Moreover, in order to further investigate the achieved accuracy of the proposed approach, we focused on the Napoli city and surroundings, including Campi Flegrei caldera [see Fig. 3(a)] where independent geodetic information (leveling and GPS measurements) was available. In particular, for our analysis, we considered pixels located in correspondence to continuous GPS stations and leveling benchmarks and, for each of these points, we compared the retrieved DInSAR time-series with those obtained from the geodetic measurements, projected on the radar LOS. Note that, although we did not perform any filtering of the atmospheric phase artifacts affecting the DInSAR time-series, there is a good agreement between the SAR and the geodetic measurements. These results f confirm the effectiveness of the proposed phase unwrapping approach.

Fig. 11. SAR data representation in the temporal/perpendicular baseline plane for the ERS-1/2 SAR data analyzed in the following experiments, relevant to the Napoli (Italy) bay area. (a) SAR image distribution. (b) Delaunay Triangulation after the removal of triangles with sides characterized by spatial and temporal baseline values, as well as doppler centroid differences, exceeding the selected thresholds (corresponding in our experiments to 400 m,

The proposed approach was validated by analyzing a dataset of 86 ERS-1/2 SAR images acquired on descending orbits (Track 36, Frame 2781) between June 8, 1992 and August 23, 2007, over the Napoli (Italy) Bay area (see Figure 11). From these data, we identified a set of 234 data pairs characterized by perpendicular baseline values smaller than 400 m, and a maximum time interval of 1500 days. Precise satellite orbital information and a 3 arcsec SRTM DEM of the area were used to generate a sequence of single-look DInSAR interferograms. The computed interferograms were then unwrapped by applying the EMCF approach. To achieve this task, we first identified the spatial grid of all the coherent pixels to be unwrapped, composed of about 530 000 pixels, where 50 000 of them are very coherent. From the latter pixels, we generated, in the spatial plane, a Delaunay triangulation. The arcs of this triangulation represent the constrained edges of the implemented CDT: this structure connects the overall set of coherent pixels , and is essential for the second EMCF PhU step, leading to the final estimate of the unwrapped interferograms on the chosen spatial grid.

Figure 12 shows a false color map of the detected mean deformation velocity, where only points with high data quality are included, superimposed on a multi-look SAR amplitude image of the area. Moreover, in order to further investigate the achieved accuracy of the proposed approach, we focused on the Napoli city and surroundings, including Campi Flegrei caldera [see Fig. 3(a)] where independent geodetic information (leveling and GPS measurements) was available. In particular, for our analysis, we considered pixels located in correspondence to continuous GPS stations and leveling benchmarks and, for each of these points, we compared the retrieved DInSAR time-series with those obtained from the geodetic measurements, projected on the radar LOS. Note that, although we did not perform any filtering of the atmospheric phase artifacts affecting the DInSAR time-series, there is a good agreement between the SAR and the geodetic measurements. These results f

confirm the effectiveness of the proposed phase unwrapping approach.

1500 days and 1000 Hz, respectively).

Fig. 12. ERS-1/2 DInSAR results. (a) Geocoded map of the mean deformation velocity of the investigated area (b) Zoomed view of the area of Napoli city and Campi Flegrei caldera (Italy) highlighted by the black box in (a). The plots show the DInSAR/leveling comparison of the deformation time-series corresponding to the pixels in (a) labeled as b (25A leveling benchmark) and c (236 leveling benchmark), respectively.

Advanced Multitemporal Phase Unwrapping Techniques for DInSAR Analyses 81

P. Berardino, G. Fornaro, R. Lanari, and E. Sansosti, "A new algorithm for surface

IEEE Trans. Geosci. Remote Sens., vol. 40, no. 11, pp. 2375–2383, Nov. 2002. A. Hooper, H. Zebker, P. Segall, and B. Kampes, "A new method for measuring deformation

Geophys. Res. Lett., vol. 31, L23611, doi:10.1029/2004GL021737, Dec. 2004. O. Mora, J. J. Mallorquí, and A. Broquetas, "Linear and nonlinear terrain deformation maps

H. A. Zebker and J. Villasenor, "Decorrelation in interferometric radar echoes," IEEE Trans.

R. Lanari, O. Mora, M. Manunta, J. J. Mallorquì, P. Berardino, and E. Sansosti, "A Small

H. A. Zebker and Y. Lu "Phase Unwrapping Algorithms for radar interferometry: residue-

R. M. Goldstein, H. A. Zebker and C. L. Werner: "Satellite radar interferometry: two-

D. C. Ghiglia and L. A. Romero: "Robust two-dimensional weighted and unweighted phase

G. Fornaro, G. Franceschetti, R. Lanari, E. Sansosti: "Robust phaseunwrapping techniques: a

M. Costantini: "A novel phase unwrapping method based on network programming", IEEE

M. Costantini, F. Malvarosa, F. Minati, L. Pietranera and G. Milillo, "A Three-dimensional

N. Miranda, B. Rosich, C. Santella, and M. Grion: "Review of the impact of ERS-2 piloting

A. Pepe, and R. Lanari R., "On the extension of the minimum cost flow algorithm for phase

W. Ku, and I. Cumming, "A Region-Growing Algorithm for InSAR Phase Unwrapping, "

F. Casu, "The Small BAseline Subset technique: performance assessment and new

Geosci. Remote Sens., vol. 44, no. 9, pp. 2374-2383, Sept. 2006.

IEEE Trans. Geosci. Remote Sens., vol. 37, pp. 124–134, Jan. 1999

Geosci. Remote Sens., vol. 30, no. 5, pp. 950–959, Sep. 1992.

dimensional phase unwrapping", Radio Sci., 23, 1988;

comparison", J. Opt. Soc. Am A., 13, 2355,1996;

vol 36, 3 May 1998, pp 813-821;

2002, vol. 3, pp. 1741-1743.

Sens., vol. 41, pp. 2243−2253, Oct. 2003.

2004.

1998;

1994;

ROM;

thesis, 2009

deformation monitoring based on small baseline differential SAR interferograms,"

on volcanoes and other natural terrains using InSAR persistent scatterers,"

from a reduced set of interferometric SAR images," IEEE Trans. Geosci. Remote

Baseline Approach for Investigating Deformation on Full resolution Differential SAR Interferograms," IEEE Trans. Geosci Remote Sens., vol. 42, no. 7, Jul.

cut, least-squares, and synthesis algorithm J. Opt Soc. Am A. 15, 586-598,

unwrapping that uses fast transform and iterative methods", J. Opt. Soc. Am. A, 11,

Phase Unwrapping Algorithm for Processing of Multitemporal SAR Interferometric Measurements," in Proc. Geosc, and Rem. Sensing Symposium, Toronto (Canada),

modes on the SAR Doppler stability", in Proc. Fringe,Frascati, Italy, Dec. 2003, CD-

unwrapping of multitemporal differential SAR interferograms," IEEE Trans.

developments for surface deformation analysis of very extended areas", PhD

### **5. Conclusion**

This chapter has presented a short review of some advanced multi-temporal phase unwrapping approaches for the generation of surface deformation time series through the application of the Small Baseline Subset (SBAS) DInSAR technique. Following a description of the basic rationale of the PhU problem, we have described in details the space-time PhU algorithm known as Extended Minimum Cost-Flow (EMCF). Finally, we focused on the recent improvements of this algorithm to analyze large interferograms affected by significant decorrelation effects and/or with severe non-linear deformation signals. The presented results clearly demonstrate the validity of these approaches and their valuable applicability in real cases.

### **6. Acknowledgment**

The ERS-1/2 SAR data were provided by the European Space Agency, and the DEM of the investigated zone has been acquired through the SRTM archive. The exploited precise orbital information was supplied by the Technical University of Delft (The Netherlands). I would like to thank all my colleagues at IREA; a special thank goes in particular to dr. Lanari, dr. Manunta, and dr. Casu who have actively contributed to the development of some of the presented algorithms.

### **7. References**


This chapter has presented a short review of some advanced multi-temporal phase unwrapping approaches for the generation of surface deformation time series through the application of the Small Baseline Subset (SBAS) DInSAR technique. Following a description of the basic rationale of the PhU problem, we have described in details the space-time PhU algorithm known as Extended Minimum Cost-Flow (EMCF). Finally, we focused on the recent improvements of this algorithm to analyze large interferograms affected by significant decorrelation effects and/or with severe non-linear deformation signals. The presented results clearly demonstrate the validity of these approaches and their valuable

The ERS-1/2 SAR data were provided by the European Space Agency, and the DEM of the investigated zone has been acquired through the SRTM archive. The exploited precise orbital information was supplied by the Technical University of Delft (The Netherlands). I would like to thank all my colleagues at IREA; a special thank goes in particular to dr. Lanari, dr. Manunta, and dr. Casu who have actively contributed to the development of

D. Massonnet, and K. L. Feigl, "Radar Interferometry and its application to changes in the

R. Bürgmann, P. A. Rosen, and E. J. Fielding, "Synthetic aperture radar interferometry to

A. K. Gabriel, R. M. Goldstein, and H. A. Zebker, "Mapping small elevation changes over

D. Massonnet, M. Rossi, C. Carmona, F. Adragna, G. Peltzer, K. Feigl, and T. Rabaute, "The

R. M. Goldstein, H. Engelhardt, B. Kamb, and R. M. Frolich, "Satellite radar interferometry

D. Massonnet, P. Briole, and A. Arnaud, "Deflation of Mount Etna monitored by spaceborne radar interferometry," Nature, vol. 375, no. 6532, pp. 567–570, Jun. 1995. G. Peltzer, and P. A. Rosen, "Surface displacement of the 17 May 1993 Eureka Valley

A. Ferretti, C. Prati, and F. Rocca, "Permanent scatterers in SAR interferometry," IEEE

Trans. Geosci. Remote Sens., vol. 39, no.1, pp. 8-20, Jan. 2001.

earth's surface," Rev. of Geophys., vol. 36, no. 4, pp. 441–500,

measure Earth's surface topography and its deformation," Annual Review Earth

large areas: Differential interferometry," J. of Geophys. Res., vol. 94, no. B7, pp.

displacement field of the Landers earthquake mapped by radar interferometry,"

for monitoring ice sheet motion: Application to an antarctic ice stream," Science,

earhtquake observed by SAR interferometry," Science, vol. 268, no. 5215, pp. 1333-

**5. Conclusion** 

applicability in real cases.

**6. Acknowledgment** 

**7. References** 

some of the presented algorithms.

doi:10.1029/97RG03139, 1998.

9183–9191, March 1989.

1336, Jun. 1995.

Planet Science, vol. 28, pp. 169–209, May 2000.

Nature, vol. 364, no. 6433, pp. 138–142, Jul. 1993.

vol. 262, no. 5139, pp. 1525–1530, Dec. 1993.


**4** 

**Simulation of 3-D Coastal Spit Geomorphology** 

Interferometric synthetic aperture radar (InSAR or IfSAR), is a geodetic technique uses two or more single look complex synthetic aperture radar (SAR) images to produce maps of surface deformation or digital elevation (Massonnet, and Feigl 1998; Burgmann et al., 2000; Hanssen 2001). It has applications as well, for monitoring of geophysical natural hazards, for instance earthquakes, volcanoes and landslides, also in engineering, in particular recording of subsidence and structural stability. Over time-spans of days to years, InSAR can detect the centimetre-scale of deformation changes (Zebker et al.,1997). Further, the precision DEMs with of a couple of ten meters can produce from InSAR technique compared to conventional remote sensing methods. Nevertheless, the availability of the precision DEMs may a cause of two-pass InSAR; regularly 90 m SRTM data may be accessible for numerous territories (Askne et al.,2003). InSAR, consequently, provides DEMs with 1-10 cm accuracy, which can be improved to millimetre level by DInSAR. Even so, alternative datasets must acquire at high latitudes or in areas of rundown coverage (Nizalapur et al., 2011). However, the baseline decorrelation and temporal decorrelation make InSAR measurements unfeasible (Lee 2001; Luo et al., 2006; Yang et al., 2007; Rao and Jassar 2010). In this regard, Gens (2000) reported the length of the baseline designates the sensitivity to height changes and sum of baseline decorrelation. Further, Gens (2000) stated the time difference for two data acquisitions is a second source of decorrelation. Indeed, the time differences while comparing data sets with a similar baseline length acquired one and 35 days a part suggests only the temporal component of the decorrelation. Therefore, the loss of coherence in the same repeat cycle in data acquisition are most likely because of baseline decorrelation. According to Roa et al. (2006), uncertainties could arise in DEM because of limitation InSAR repeat passes. In addition, the interaction of the radar signal with troposphere can also induce decorrelation. This is explained in several studies of

Hanssen (2001); Marghany and Hashim (2009); and Rao and Jassar (2010).

Generally, the propagation of the waves through the atmosphere can be a source of error exist in most interferogram productions. When the SAR signal propagated through a

**1. Introduction** 

**Using Differential Synthetic Aperture** 

*Institute for Science and Technology Geospatial (INSTEG) Universiti Teknologi Malaysia, Skudai, Johor Bahru,* 

**Interferometry (DInSAR)** 

Maged Marghany

*Malaysia* 

