Topological Characterization and Advanced Noise-Filtering Techniques for Phase Unwrapping of Interferometric Data Stacks

Abstract

This chapter addresses the problem of phase unwrapping interferometric data stacks, obtained by multiple SAR acquisitions over the same area on the ground, with a twofold objective. First, a rigorous gradient-based formulation for the multichannel phase unwrapping (MCh-PhU) problem is systematically established, thus capturing the intrinsic topological character of the problem. The presented mathematical formulation is consistent with the theoretical foundation of the discrete calculus. Then within the considered theoretical framework, we formally describe an innovative procedure for the noise filtering of time-redundant multichannel multilook interferograms. The strategy underlying the adopted multichannel noise filtering (MCh-NF) procedure arises from the key observation that multilook interferograms are not fully time consistent due to multilook operations independently applied on each single interferogram. Accordingly, the presented MCh-NF procedure suitably exploits the temporal mutual relationships of the interferograms. Finally, we present some experimental results on real data and show the effectiveness of our approach applied within the well-known small baseline subset (SBAS) processing chain, thus finally retrieving the relevant Earth’s surface deformation time series for geospatial phenomena analysis and understanding.

1. Introduction

Multichannel (or multitemporal) InSAR techniques address the processing of interferometric data stacks obtained by combining multiple SAR acquisitions over the same area. These approaches can be essentially categorized in two main classes: persistent scatterers (PS) and small baseline (SB)-based techniques. The solution of the multichannel phase unwrapping (MCh-PhU) problem is generally required in multichannel InSAR techniques, whenever multidimensional SAR data set, conveying information about complex Earth’s crust events, have to be systematically investigated on suitable space-time scales for geospatial phenomena understanding [1–29]. In this chapter, we focus on two different related main issues.

Primarily, we present a rigorous gradient-based formulation of the MCh-PhU problem that is consistent with the theoretical foundation of the discrete calculus [30–34]. Emphasis is placed on the topological characterization of the underlying discrete setting provided by the differential operators of the discrete calculus, which are formally associated with matrix operators and represent the discrete counterparts of the classical differential operators of the vector calculus. Accordingly, MCh-PhU problem formulation is systematically established in terms of discrete differential operators, which are defined by the topology of the intrinsically discrete spaces upon which they act, thus capturing the essential topological character of the problem within a systematic matrix formalism [35]. It is worth highlighting that our approach provides an unambiguous and theoretical-consistent formalism for the MCh-PhU problem, overcoming the conceptual inconsistencies of the existing gradient-based formulations [1, 17, 29]. Indeed, the existing approaches pose some conceptual limitations from a mathematical viewpoint since they rely on an intrinsically discrete setting based description and, at the same time, resort to the concepts of gradient and curl of the conventional vectorial calculus, which inherently imply a reference to an underlying continuum space and the notion of the infinitesimal [30]. In addition, the proposed formal framework enables meaningful analytical investigations on a mathematical consistent playground, also providing interesting implications and permitting to express previous obtained results in a more general way.

Then we present an innovative procedure to filter out the noise affecting the phase components of a redundant set of (multitemporal) multilook small-baseline interferograms. This is achieved by independently solving, for each pixel of the scene, a nonlinear optimization problem based on computing the wrapped phase vector that minimizes the (weighted) circular variance of the difference between the original and noise-filtered interferograms [43]. This noise-filtering procedure arises from the key observation that multilook interferograms are not fully time consistent because they are generated through multilook operations that are independently carried out on each single interferogram. Indeed, the wrapped discrete curl of the interferometric phases defined on a graph whose nodes and edges describe SAR acquisitions (in the time/perpendicular baseline domain) and inherent interferograms, respectively, is different from zero. This modulo-2π cyclic inconsistency of multichannel interferometric phases is properly handled by the presented multichannel noise-filtering (MCh-NF) procedure. The presented technique is very easy to implement because it does not imply a preliminary time-consuming selection of statistically homogenous pixels (SHP), as for instance required by the SqueeSAR technique [44], and it has no need of any a priori information on the statistics of complex-valued SAR images. The effectiveness of the presented noise-filtering approach as well as its impact on the quality of multichannel phase unwrapping procedures are also fully investigated.

2. Multichannel phase unwrapping problem

In this Section, we review the mathematical formulation of the multichannel phase unwrapping (MCh-PhU) problem within the purview of the discrete calculus. As a matter of fact, a graph-based description naturally arises in formulating the MCh-PhU problem due to the underlying discrete irregular data structure. Indeed, as far as discrete settings (e.g., graphs) are concerned, resorting to the conventional vectorial calculus might not be adequate since it inherently implies a reference to an underlying continuum space. On the contrary, discrete calculus offers a rigorous methodological framework since it treats a discrete domain as entirely its own entity. In particular, discrete calculus provides proper differential operators that make it possible to purely operate onto a finite, discrete structure without referring to the continuous space and notion of the infinitesimal [30]. More specifically, the introduction of some well-known algebraic structures [30–33] capturing the essential topological character of the underlying graphs permits to phrase pertinent differential operators as matrices. Therefore, one of the most important consequences of this approach is that the purely topological nature of the discrete differential operators is made more apparent and concrete. Accordingly, by systematically adopting the relevant key concepts and tools available within this theoretical framework, we here provide a description of the MCh-PhU problem on a suitable mathematical playground. For such a purpose, we first establish the notation and terminology used throughout the subsequent Sections 2.1, and then the problem at hand is reviewed within this formalism (Sections 2.2 and 2.3). We remark that the focus here is on presenting key concepts that are useful for the following analyses; however, a comprehensive treatment of the discrete calculus and related huge fields of mathematics (e.g., algebraic topology, exterior calculus, and differential forms) is clearly beyond the scope of this work but can be found in refs. [30–33].

2.1. The theoretical framework of discrete calculus

A graph $$G\text{(}V\text{;}ℰ\text{)}$$ is defined by two sets: $$V$$ and $$ℰ$$. The former is the set of nodes (or vertices) of the graph, and the latter represents the corresponding set of edges. Let Q and M be the cardinality of $$V$$ and $$ℰ$$, respectively. The vector space $${ℝ}^M$$ is referred to as the edge space, and the vector space $${ℝ}^Q$$ is referred to as the vertex space, with $$ℝ$$ denoting the field of real numbers. Without loss of generality, we here assume that the graph $$G$$ is connected (i.e., every pair of vertices in the graph is connected [33]). Moreover, an orientation establishes a default direction on an edge that is considered positive or negative, thus yielding an oriented graph. The M × Q incidence matrix $$\Pi =[{\Pi }_{mq}]$$ of an oriented graph $$G$$ specifies its edge–node connectivity relations, whose entries are defined as follows [30–33]:

$${\Pi }_{mq}\equiv \{\begin{array}{c}-1,\text{if }m\text{th edge starts at }q\text{th node}\\ +1,\text{if }m\text{th edge ends at }q\text{th node}\\ 0,\text{otherwise}\end{array}$$E1

with m = 1, 2,..., M and q = 1, 2,..., Q. It is important to note that the column rank of $$\Pi$$ is Q − 1. The incidence matrix $$\Pi$$ generates an orthogonal decomposition $$ℝ\text{(}\Pi \text{)}\oplus N\text{(}{\Pi }^{\text{T}}\text{)}={ℝ}^M$$, where $$ℝ\text{(}\Pi \text{)}$$ is the column space of the incidence matrix $$\Pi$$, and $$N\text{(}{\Pi }^{\text{T}}\text{)}$$ denotes the kernel (or null-space) of the matrix $${\Pi }^{\text{T}}$$. The notion of cycle space is also fundamental in graph theory. The cycle space of the graph $$G$$, namely, $$C=C(G)$$, is the subspace of edge space $${ℝ}^M$$ spanned by all the cycles in $$G$$. The dimension of $$C(G)$$ is also referred to as the cyclomatic number of $$G$$ [33]. It is also well known that, for every connected graph $$G$$ with Q nodes and M edges, the dimension of the cycle space is given by R = dim($$C$$ ($$G$$)) = M − Q + 1 [30]. Each basis for $$C(G)$$ (i.e., the cycle basis) is therefore uniquely specified by an M × R matrix $$\Omega$$, called cycle matrix. Thus, the column vectors of $$\Omega =\text{[}{\omega }^1,\dots ,{\omega }^R]$$ form a basis for an R-dimensional vector subspace (the cycle space of $$C(G)$$ of $${ℝ}^M$$. $$C(G)$$ is indeed the null-space of $${\Pi }^{\text{T}}$$ so that a cycle basis provides a basis for $$N\text{(}{\Pi }^{\text{T}}\text{)}$$ [30]. Accordingly, a fundamental property of a linear graph is expressed by the remarkable relations:

$${\Pi }^{\text{T}}\Omega =0$$E2

$${\Omega }^{\text{T}}\Pi =0$$E3

Indeed, several methods for defining a cycle set have been studied, and they can be used to define incidence relations between edges and cycles. Specifically, the definition of a cycle set from the edge set can be obtained algebraically and geometrically (i.e., from an embedding). Algebraic methods find a suitable M × R matrix $$\Omega$$ whose columns provide a basis for the null-space of $${\Pi }^{\text{T}}$$, with $$R=\dim (N({\Pi }^{\text{T}}))$$. Geometric methods for defining a cycle set (i.e., from an embedding) permit to identify algorithmically a cycle set (representing the faces) in this embedding and may be used to produce the edge–face incidence matrix $$\Omega$$ (as illustrated in Figure 1). In particular, it is possible to consider a normalized irreducible cycle basis forming elementary (or irreducible) cycles [30–33], i.e., cycles that contain no other cycles, so that we can associate to each elementary cycle an elementary cycle vector $${\omega }^r={[{\omega }_1,\mathrm{...},{\omega }_M]}^{\text{T}}$$, whose entries are defined as follows:

$${\omega }_i\equiv \{\begin{array}{c}+1\text{if }r\text{th cycle traverses }i\text{th edge forward }\\ -1\text{if }r\text{th cycle}\text{traverses }i\text{th edge backward}\\ 0\text{otherwise }\end{array}$$E4

Accordingly, the so defined $$\Omega =\text{[}{\omega }^1,\dots ,{\omega }^R]$$ provides a particularly attractive basis for $$N({\Pi }^{\text{T}})$$, i.e., the cycle basis formed by all the elementary cycle vectors associated with the elementary cycles in $$G$$. We also note that $$\Omega$$ defines the incidence connectivity relations between edges and cycles (see Figure 1).

It is instructive to highlight that the topological operators $$\Pi$$, $${\Pi }^{\text{T}}$$, and $$\Omega$$ provide the discrete counterparts of the classical gradient ($$\nabla$$), divergence ($$\nabla \cdot$$), and curl ($$\nabla \times$$) operators of the vector calculus for continuous space, respectively. Accordingly, they can be regarded as differential operators on the discrete setting [30]. In addition, it is worth emphasizing that identities (2) and (3) mimic the properties of their classical vector calculus counterparts $$\nabla \cdot \nabla \times =0$$ (div curl = 0) and $$\nabla \times \nabla =0$$ (curl grad = 0), respectively. It should be pointed out that $$\Pi$$ yields differences along edges of nodal “potentials” represented by $$\Pi x$$. Conversely, given an arbitrary $$f\in {ℝ}^M$$, a solution of the equation $$\Pi x=f$$ (if it exists) is called the potential of $$f$$. Note also that $$x$$ (if it exists) is not unique since the constant column $$\mathbb{I}={[1,\mathrm{...},1]}^{\text{T}}\in {ℝ}^Q$$ is an element of the kernel of $$\Pi$$. Of course, not every $$f\in {ℝ}^M$$ is the discrete gradient of some $$x$$ since $$f$$ may contain a curl component. Indeed, a prescribed $$f\in {ℝ}^M$$ can be written as a nodal difference ($$f=\Pi x$$) if it is cyclically consistent, i.e., if it satisfies $${\Omega }^{\text{T}}f=0$$ (i.e., there is no component of the flow in the cycle space). Note also that $$\Omega$$ is the (cross) differential operator of the graph whose expression can be given in terms of a normalized cycle basis; $$N({\Omega }^{\text{T}})$$ denotes the subspace of $${ℝ}^M$$ with zero flow circulation (curl-free) around cycles. Moreover, $${\Pi }^{\text{T}}f$$ yields nodal accumulations from flows along edges. As a result, the differential operators, as basic tools of the discrete calculus, have been established and properly phrased on the discrete space. This mathematical abstraction meaningfully captures the topological structure of the underlying discrete setting. Note that the topological characterization of the graph is essentially embodied in the algebraic structure of the associated discrete (matrix) operators and their interrelations. We also stress the significant distinction between the discrete operators and the commonly adopted discretized versions of the continuous differential operations obtained via the method of finite differences in numerical analysis; the latter generally lack the desirable topological behaviour [42].

Figure 1.

An example graph shown along with its edge–node incidence matrix, $$\Pi$$, and cycle (edge-face incidence) matrix, $$\Omega$$. Note that M=5, Q=4 and R=2.

As a final remark, some considerations on the total unimodularity (TU) property inherent to the matrix operators, which is extremely important in integer linear programming, are in order. We recall that a matrix is TU if the determinant of every submatrix is either zero or ±1. For any graph, the edge–node incidence matrix is TU. On the contrary, the face–edge incidence matrix, in general, is not TU [30–34]. However, the edge–face incidence matrix is TU when each edge is included in exactly two faces that traverse the edge in opposite directions (e.g., a planar graph with a minimum cycle basis [30]). In this circumstance, total unimodularity of the edge–face incidence matrix stems from the fact that the face–edge incidence matrix is the edge–node incidence matrix of the dual graph [30]. Indeed, a TU constraint matrix (and integer constraints) guarantees that the solution of the related optimization problem (see also optimization problems (21) and (22) in the following) will be integer. Nonetheless, TU property has a further practical significance since the relaxed problem, obtained by neglecting the integer constraints, can also be solved using generic linear (not integer) programming solvers.

2.2. Rigorous gradient-based formulation of the MCh-PhU problem

Once the basic concepts of the discrete calculus and graph theory are presented, we are in the position to frame the formulation of the MCh-PhU problem on an appropriate mathematical playground. Let us consider a set of Q single-look-complex (SLC) SAR data acquired over a certain area of interest. One of them is assumed as the reference (master) image, with respect to which the images are properly coregistered. This set is characterized by the corresponding acquisition times $$\{t_1,\dots ,t_Q\}$$ and perpendicular baselines $$\{b_{\perp 1},\dots ,b_{\perp Q}\}$$. Accordingly, for each coregistered SLC pair, a multilook phase interferogram (suitably depurated by the flat-Earth and topography contributions, by using a priori information about the topography and acquisition geometry) can be produced [37]; however, a common practice (within the multitemporal SB InSAR class) is first to identify a suitable small-baseline subset of the relevant multibaseline (temporal and spatial-perpendicular baselines) interferometric-pair set [6]. This is done to confine the effect of decorrelation phenomena associated with inherent angular and temporal electromagnetic backscattering variations [38]. Furthermore, a subset of common pixels of the M interferograms are then usually identified via the estimated coherence [37,38], so that only P pixels characterized by relatively high coherence values are singled out. In other words, the coherence index, providing a quantitative estimation of the decorrelation effects, permits to discriminate in favor of the “reliable” pixels.

The final aim is to reconstruct the absolute (i.e., not restricted in the principal [–π,π) interval) interferometric phases values from the wrapped (i.e., observed only in the principal [–π,π) interval) interferometric phase pertinent to M multichannel interferograms.

The problem we are interested in can be naturally framed on a discrete setting. Indeed, one possibility is to regard the discrete set of P selected (typically sparsely distributed) coherent pixels as a set of nodes, $V_{\text{B}}$, in the Euclidean (azimuth, range) plane, and the set of Q SAR acquisitions representing a set of nodes, $V_{\text{A}}$, in the Euclidean (temporal-baseline, perpendicular-baseline) plane [26]. Accordingly, a formal description of the problem at hand can be given in terms of a couple of abstract graphs: $G_{\text{A}}=(V_{\text{A}};{ℰ}_{\text{A}})$ and $G_{\text{B}}=(V_{\text{B}};{ℰ}_{\text{B}})$ where the corresponding edge sets ${ℰ}_{\text{A}}$ and ${ℰ}_{\text{B}}$ have to be properly defined. Accordingly, with Q and M, we denote the cardinality of $V_{\text{A}}$ and ${ℰ}_{\text{A}}$, and with P and N the cardinality of $V_{\text{B}}$ and ${ℰ}_{\text{B}}$, respectively. Note also that defining ${ℰ}_{\text{A}}$ pertains to the M interferometric pairs selection. Defining a meaningful edge set for a collection of nodes is now concerned since different criteria can be adopted to achieve it. The dimensionality of the ambient space in which the graph is embedded deserves some considerations. In this regard, we recall that a graph is called planar if it can be embedded in the plane [30–33]. Note also that a graph is not generally guaranteed to be planar, even if the nodes are embedded in two dimensions. Since planarity is important for the workability of the implicated optimization procedure with powerful numerical solvers, a typical option to preserve graph planarity is resorting to the Delaunay triangulation in the Euclidean plane for establishing the edge set from nodes embedded in two dimensions [26]. Note that such an option specifically pertains to the solution strategy [26, 35]; nonetheless, our general formulation applies as well when different edge structures are adopted. Accordingly, once $G_{\text{A}}$ and $G_{\text{B}}$ have been somehow defined, the topological properties inherent to each graph are algebraically captured by the related differential operators, which are summarized in Table 1.

Symbol	Quantity	Related meaning
Q	Nodes number of $G_{\text{A}}$	Number of SAR acquisitions
M	Edges Number of $G_{\text{A}}$	Number of interferometric pairs
R	Dimension of the cycle space of $G_{\text{A}}$
$${\Pi }_{\text{A}}$$	M × Q incidence matrix of $G_{\text{A}}$	Discrete-gradient operator
$${\Pi }_{\text{A}}^{\text{T}}$$		Discrete-divergence operator
$${\Omega }_{\text{A}}^{}$$	M × R cycle matrix of $G_{\text{A}}$	Discrete-curl operator
P	Nodes number of $G_{\text{B}}$	Number of selected pixels
N	Edges number of $G_{\text{B}}$
L	Dimension of the cycle space of $G_{\text{B}}$
$${\Pi }_{\text{B}}$$	N × P incidence matrix of $G_{\text{B}}$	Discrete-gradient operator
$${\Pi }_{\text{B}}^{\text{T}}$$		Discrete-divergence operator
$${\Omega }_{\text{B}}^{}$$	N × L cycle matrix of $G_{\text{B}}$	Discrete-curl operator

Table 1.

Adopted notation

First of all, we consider the absolute phase relevant to the multichannel SAR acquisition as a node variable pertinent to both the graphs $G_{\text{A}}$ and $G_{\text{B}}$ ; by using a matrix representation, this information can be conveniently arranged in a Q × P matrix $$\Phi$$ as follows:

$$\Phi =[{\phi }^1,,\dots ,{\phi }^P]=\left[\begin{array}{c}\begin{array}{c}{\phi }_1\\ ⋮\\ {\phi }_Q\end{array}\end{array}\right]$$E5

where $$\forall p\in \left\{1,2,\mathrm{...},P\right\}$$ $${\phi }^p\in {ℝ}^Q$$ encodes in a vectorial manner the pth node variable relevant to the graphs $G_{\text{A}}$ ; similarly, $$\forall q\in \left\{1,2,\mathrm{...},Q\right\}$$ $${\phi }_q\in {ℝ}^P$$ encodes the qth node variable relevant to $G_{\text{B}}$.

Widely adopted global gradient-based PhU approaches, which have historically been developed for the single-channel case, generally consist in three processing steps [1, 29]. First, an estimation of the (wrapped) phase gradient is obtained; the estimated phase gradient is then suitably corrected (in terms of 2π multiples),and subsequently integrated to attain the unwrapped (absolute) phase.

Within the formulation of the MCh-PhU problem we concern [26], a twofold estimation of the discrete gradient field is carried out onto the considered two graphs $G_{\text{A}}$ and $G_{\text{B}}$, as discussed in the following. The stack of the absolute interferometric phases relevant to the M (vectorized) interferograms can be formally represented through a P × M matrix denoted by

$$\Psi =[{\psi }^1,\dots ,{\psi }^M]=\left[\begin{array}{c}\begin{array}{c}{\psi }_1\\ ⋮\\ {\psi }_P\end{array}\end{array}\right]$$E6

wherein the P-dimensional vector $${\psi }^m$$ refers to the absolute phase field pertinent to the mth interferogram. Accordingly, $$\Psi$$ is formally related to the absolute (unwrapped) phase matrix $$\Phi$$ via the discrete gradient operator $${\Pi }_{\text{A}}$$ :

$${\Psi }^{\text{Τ}}={\Pi }_{\text{A}}\Phi$$E7

Note also that

$${\Psi }^{\text{Τ}}=[{\psi }_{\text{1}}^{\text{T}},\dots ,{\psi }_P^{\text{T}}]=[{\Pi }_{\text{A}}{\phi }^1,\dots ,{\Pi }_{\text{A}}{\phi }^P]$$E8

By applying the discrete gradient operator $${\Pi }_{\text{B}}$$ to the absolute phase of each interferometric pair, we obtain

$${\Pi }_{\text{B}}\Psi =[{\Pi }_{\text{B}}{\psi }^1,\dots ,{\Pi }_{\text{B}}{\psi }^M]$$E9

As a result, by using Eq. (7), we get

$${\Pi }_{\text{B}}\Psi ={\Pi }_{\text{B}}{\Phi }^{\text{T}}{\Pi }_{\text{A}}^{\text{T}}$$E10

Second, we consider the (wrapped) interferometric phase that is uniquely defined only in the principal value range since it is obtained as the phase of a complex function. Hence, it is convenient to formally introduce the non-injective (modulo-2π) wrapping operator $$W:\phi \in ℝ\to \mathrm{mod}\left(\phi +\pi ,2\pi \right)-\pi \in [-\pi ,\pi )$$. It should be noted that the following trivial identities hold:

$$\begin{array}{l}W(W(A)\pm W(B)\text{)}=W(A\pm B\text{)}=W(W(A)\pm B\text{)}=W(A\pm W(B)\text{)}\text{a}\\ W\text{(}A\text{)}=A+\text{2}\pi Z\text{b}\end{array}$$E11

where $$A$$ and $$B$$ represent two generic matrices and $$Z$$ is a suitable integer matrix. Given the stack of the unknown absolute (unwrapped) interferometric-phases $$\Psi$$, the corresponding stack of the wrapped phases $$\tilde{\Psi }$$ can be conveniently expressed in terms of the wrapped discrete gradient of (wrapped) observed phase as follows:

$${\tilde{\Psi }}^{\text{Τ}}=[W({\Pi }_{\text{A}}W\text{(}{\phi }^1)),\dots ,W({\Pi }_{\text{A}}W\text{(}{\phi }^P))]=W({\Pi }_{\text{A}}W\text{(}\Phi ))=W({\Pi }_{\text{A}}\Phi )$$E12

where we have exploited Eq. (11b). Note also that the pth column of $${\tilde{\Psi }}^{\text{Τ}}$$, i.e., $${\tilde{\psi }}_p^{\text{T}}=W({\Pi }_{\text{A}}W\text{(}{\phi }^p))$$, can be regarded as an estimate of the absolute-phase discrete gradient on the graph $G_{\text{A}}$. It is worth remarking that the observed (multilook) interferometric phase can however be corrupted by noise [39–41], which is taken into account by considering an additive phase noise term $$D$$. Accordingly, by using Eq. (11a), we get

$${\tilde{\Psi }}^{\text{Τ}}=W(W({\Pi }_{\text{A}}\Phi )+D\text{)}=W({\Pi }_{\text{A}}\Phi +D\text{)}$$E13

More specifically, whenever a possible spatial filtering (e.g., conventional multilooking followed by a noise-filtering step [42]) is independently applied to each SAR interferometric data pair, the resulting term $$D$$ in Eq. (13) implies that the phase interferograms $${\tilde{\psi }}^m$$, with $$m\in \left\{1,2,\mathrm{...},M\right\},$$ are no more fully time consistent (in the sense of [26, 43, 44]). To clarify this point, we observe that by using Eq. (13), it turns out that

$$W\left({\Omega }_A^T{\Psi }^T\right)=W\left({\Omega }_A^T\left({\Pi }_A\mathrm{+D+}2\pi Z\right)\right)=W\left({\Omega }_A^T{\Pi }_A\Phi +{\Omega }_A^{T}D\right)=W\left({\Omega }_A^{T}D\right)\ne 0$$E14

where we have used Eq. (11b) with $$A={\Pi }_{\text{A}}\Phi +D$$ and noted that $${\Omega }_{\text{A}}^{\text{T}}Z$$ is an integer matrix and $${\Omega }_{\text{A}}^{\text{T}}{\Pi }_{\text{A}}=0$$ (according to Eq. (3)). Eq. (14) reads as “the wrapped discrete curl of the interferometric phase on $G_{\text{A}}$ (i.e., pertinent to the ‘temporal’ domain) is generally different from zero”; it formally expresses the (modulo-2π) cyclic inconsistency of the multichannel interferometric phase inherent to independently filtered SAR interferograms. Note also that Eq. (14) represents, within our framework, the generalization (to a wider class of discrete settings) of the “phase triangularity” condition in ref. [44], capturing the underlying structure of the problem within a suitable matrix formalism. With reference to the mth interferometric pair, the estimated absolute interferometric-phase gradient on the graph $G_{\text{B}}$ is then obtained by wrapping the discrete gradient of (wrapped) interferometric-phase field: $$g^m=W\text{(}{\Pi }_{\text{B}}{\tilde{\psi }}^m\text{)}$$. Thus, by stacking the so-obtained absolute phase gradient estimations, we get the N × M matrix $$G=[g^1,g^2,\dots ,g^M\text{]}$$, where

$$G=W\text{(}{\Pi }_{\text{B}}\tilde{\Psi }\text{)}$$E15

Finally, by substituting Eq. (13) in Eq. (15), we obtain

$$G=W\left({\Pi }_{\text{B}}\tilde{\Psi }\right)=W\left({\Pi }_{\text{B}}W\left({\left[{\Pi }_{\text{A}}\Phi \right]}^{\text{T}}+D^{\text{T}}\right)\right)$$E16

From Eq. (16), by using Eq. (11b), we get

$$G=W\left({\Pi }_{\text{B}}{\Phi }^{\text{T}}{\Pi }_{\text{A}}^{\text{T}}+{\Pi }_{\text{B}}{D}^{\text{T}}+\text{2}\pi {\Pi }_{\text{B}}Z\right)=W\left({\Pi }_{\text{B}}{\Phi }^{\text{T}}{\Pi }_{\text{A}}^{\text{T}}+{\Pi }_{\text{B}}{D}^{\text{T}}\right)$$E17

where in the last equality we have noted that $${\Pi }_{\text{B}}Z$$ is also an integer matrix. It should be emphasized that, under the assumption $$D=o$$, the equality between Eqs. (10) and (17) holds only up to an integer matrix multiplied by $$2\pi$$.

2.3. The MCh-PhU problem as constrained optimization

In this Section, the nonlinear inversion MCh-PhU problem is reformulated as a (nonlinear) constrained optimization problem. According to the presented general formulation, we introduce in the following the MCh-PhU problem as the solution of the following matrix equation:

$${\Pi }_{\text{B}}{\Phi }^{\text{T}}{\Pi }_{\text{A}}^{\text{T}}+{\Pi }_{\text{B}}{D}^{\text{T}}+2\pi K=G$$E18

where the columns of $$G$$ represent the interferometric-phase pseudo-gradients estimated from the observed phase, and $$K$$ is an (unknown) N × M integer matrix, whose columns represent the corresponding (2π-normalized) corrections to be added to the (wrapped) interferometric-phase pseudo-gradients in order to recover the absolute interferometric-phase discrete gradients. It is worth noting that the term pseudo-gradient is used here to emphasize that the integration of the estimated gradient is path-dependent (non-conservative behavior); the term residues [1] is also typically used to connote the inconsistency of the estimated phase gradient. As a matter of fact, matrix equation (Eq. (18)) describes an ill-posed problem, in which the data $$G$$ generally do not constrain sufficiently the problem to get a unique solution. Additional suitable constraints and a priori assumption have, thus, to be introduced to solve the problem. First, for restoring the cyclic consistency (see Section 2.1) of the estimated pseudo-gradients pertinent to the graphs $G_{\text{B}}$ and $G_{\text{A}}$, two corresponding sets of (equality) constraints have to be enforced, respectively. More specifically, pre-multiplying both sides of Eq. (18) by $${\Omega }_{\text{B}}^{\text{T}}$$ and taking into account Eq. (3), we obtain

$${\Omega }_{\text{B}}^{\text{T}}\left(G-2\pi K\right)=0$$E19

Similarly, by premultiplying both sides of the transposed version of Eq. (18) by $${\Omega }_{\text{A}}^{\text{T}}$$ and taking also into account Eq. (3), we obtain

$${\Omega }_{\text{A}}^{\text{T}}\left(G-{\Pi }_{\text{B}}{D}^{\text{T}}-2\pi K\right){}^{\text{T}}=0$$E20

Constraints stated by Eq. (19) imply that the columns of $$G-2\pi K$$ must lie in the null-space of $${\Omega }_{\text{B}}^{\text{T}}$$. Since the matrix $${\Pi }_{\text{B}}$$ represents a basis to span the null-space of $${\Omega }_{\text{B}}^{\text{T}}$$ (see Eq. (3)), we may then write $$G-2\pi K={\Pi }_{\text{B}}X$$, where $$X$$ is a new variable. Accordingly, the corrected pseudo-gradients stack $$G-2\pi K$$ is enforced to be a stack of discrete gradients, which can thus be unambiguously integrated. Similarly, Eq. (20) implies $${\text{[}G-{\Pi }_{\text{B}}{D}^{\text{T}}-2\pi K\text{]}}^{\text{T}}={\Pi }_{\text{A}}^{}Y$$. As a result, the two sets of constraints, stated by Eqs. (19) and (20), guarantee that the solution of the problem is effective in preserving the cyclic consistency (curl-free) property of the corrected gradients pertaining to the graphs $G_{\text{B}}$ and $G_{\text{A}}$, respectively. As a matter of fact, the solution of Eq. (18) cannot be determined by using the two sets of constraints (Eqs. (19) and (20)) only; thus, the inverse problem must be first regularized [45]. The minimum-norm methods search for a global solution that minimizes a generalized error-norm associated with an optimality criterion, so incorporating prior information about the behavior of the solution [1]. Accordingly, we resort to a regularization approach using l₁-norm minimization in weighted version, as a specific case of l_p-norm general formulation. Formally, the MCh-PhU problem may be then formulated as a constrained optimization problem for the field of integer corrections:

$$\widehat{K}=\underset{K\in {\mathbb{Z}}^{N\times M}}{\arg \min }{\Vert K\Vert }_{1,C}$$E21

subject to

$$\{\begin{array}{c}{\Omega }_{\text{A}}^{\text{T}}{K}^{\text{T}}={\Omega }_{\text{A}}^{\text{T}}{\left[G-{\Pi }_{\text{B}}{D}^{\text{T}}\right]}^{\text{T}}{\left(2\pi \right)}^{-1}\\ {\Omega }_{\text{B}}^{\text{T}}K={\Omega }_{\text{B}}^{\text{T}}G{\left(2\pi \right)}^{-1}\end{array}$$E22

wherein

$${\Vert K\Vert }_{1,C}={\displaystyle \sum _{m=1}^M{\displaystyle \sum _{n=1}^N{c}_{nm}}\left|k_{nm}\right|}$$E23

represents the weighted l₁-norm [46] of the matrix $$K$$, $$C={[c_{nm}]}_{N\times M}$$ denotes a suitable weighting matrix, and $$\mathbb{Z}$$ indicates the field of integer numbers. As far as the existence of an integer solution for Eqs. (21) and (22) is concerned, it should be noted that the considerations at the end of Section 2.1 apply. Since the first matrix equation in Eq. (22) includes a generally not null (unwanted) term $${\Omega }_{\text{A}}^{\text{T}}D$$, its fulfillment deserves further discussion. Although the evaluation of $$W({\Omega }_{\text{A}}^{\text{T}}D)$$ can be obtained according to Eq. (14), however, a full estimation for $${\Omega }_{\text{A}}^{\text{T}}D$$ is generally not a simple task. Further discussion is provided in Section 3. The solution of the optimization problem (Eqs. (21) and (22)) is also referred to as the minimum weighted discontinuity solution (in a weighted l₁-norm sense) [1, 23]. As a matter of fact, finding the global minimum point of the problem stated by Eqs. (21) and (22) for an arbitrary pair of graphs is, in general, a difficult task. A suboptimal strategy aimed at solving Eqs. (21) and (22) consists in adopting a two-stage approach. This is, in particular, the solution strategy implemented through the extended minimum cost flow (EMCF) technique [26], in which the edge structure of each considered graph is usually defined via a Delaunay triangulation in the Euclidean plane, to take advantage from efficient solvers for minimum cost flow (MCF) problems [47, 49, 50]. We remark that the distinctive characteristic of the EMCF approach is the extensive use of the computationally efficient MCF method. Moreover, a dual-level parallel model for EMCF has also been proposed in refs. [35] and [36]. Moreover, different approaches toward full 3D phase unwrapping have recently been proposed in refs. [63] and [64].

3. Noise-filtering of multichannel SAR interferograms

In this Section, we review the basic concepts concerning the filtering of noise that corrupts a stack of multitemporal SAR interferograms. First, the noise-filtering operation for single-channel multilook interferograms is discussed; subsequently, the general framework of the multichannel noise-filtering (MC-NF) approach, which is intimately connected with the problem of multichannel phase unwrapping, is described.

3.1. Decorrelation noise in SAR interferograms

In order to introduce the problem at hand, let us first consider one single-channel SAR interferogram obtained starting from two SAR (synthetic aperture radar) images, namely, $$i_1$$ and $$i_2$$, acquired (over the same scene on Earth) at two different times, namely, $$t_1$$ and $$t_2$$, respectively. The two SAR images can be represented via two complex-valued signals, say $$i(x,r)$$ and $$i_2(x,r)$$, with $$x$$ and $$r$$ denoting the two independent spatial variables (with respect to azimuth and range direction, respectively) in the radar geometry. The two complex signals can be expressed as follows [29, 51]:

$$\begin{array}{l}{i}_1(x,r)={\gamma }_1(x,r)e^{-j\frac{4\pi }{\lambda }r}+n_1(x,r)\\ i_2(x,r)={\gamma }_2(x,r)e^{-j\frac{4\pi }{\lambda }\left(r+\delta r\right)}+n_2(x,r)\end{array}$$E24

where $$\delta r$$ is the sensor-to-target slant range difference at time $$t_2$$ with respect to time $$t_1$$, and $${\gamma }_1(x,r)$$ and $${\gamma }_2(x,r)$$ are the corresponding (complex-valued) reflectivity functions of the illuminated scene at time $$t_1$$ and $$t_2$$, respectively. Furthermore, the two additive (noise) contributions $$n_1(x,r)$$ and $$n_2(x,r)$$ describe random quantities that are included in Eq. (24). As a result of these noise terms and of the intrinsic random nature of the two images reflectivity functions, when the two SAR images are interfered to form a so-called interferogram, i.e., when their phase difference $$\psi$$ is extracted, the interferometric phase will be noisy. An important parameter influencing the quality of the retrieved interferometric phase is the (complex) cross-correlation factor between the two involved SAR images, which is typically defined as

$$\chi =\frac{E[i_1\cdot i_2^{*}]}{\sqrt{E\left[{\left|i_1\right|}^2\right]\cdot E\left[{\left|i_2\right|}^2\right]}}=\rho e^{j\psi }$$E25

where $$\rho \in [0,1]$$, $$\psi \in [-\pi ,\pi )$$, and the asterisk denotes the conjugate complex value. Noteworthy, the cross-correlation factor (Eq. (25)) is a complex-valued term that can be decomposed in terms of amplitude $$\rho$$ (i.e., $$\rho =\left|\chi \right|$$) and phase $$\psi$$. For interferometric SAR images, $$\chi$$ can be evaluated by performing spatial averaging (known as multilooking) operations on a statistically homogeneous area. Indeed, the symbol $$E[]$$ in Eq. (25), which is representative of the statistical expectation operation [52], can then be replaced by the spatial averaging operation. The amplitude factor $$\rho$$, which is known to as coherence, accounts for the similarity between the two SAR images, whereas $$\psi$$ is the multilook interferometric phase. A value for the coherence that approaches zero is representative of an uncorrelated scene, whereas coherence value that is close to unity corresponds to a noise-free interferogram.

There are several causes that are responsible for coherence decrease. As matter of fact, the cross-correlation factor in Eq. (25) depends on different noise sources, and it can be conveniently factorized as follows [51]:

$\chi ={\chi }_{\text{the}}\cdot {\chi }_{\text{tem}}\cdot {\chi }_{\text{spa}}\cdot {\chi }_{\text{dop}}\cdot {\chi }_{\text{mis}}\cdot {\chi }_{\text{vol}}\cdot e^{j\psi }$E26

where

• ${\chi }_{\text{the}}$ is the contribution of the thermal noise.

• ${\chi }_{\text{temp}}$ accounts for the effects due to (temporal) changes in the complex-valued reflectivity function between the two passages of the radar sensor over the illuminated area. The so-called temporal decorrelation is very difficult to be statistically modeled being associated to complex modifications of the electromagnetic response of the scene: They can be induced by human activities and/or natural causes.

• ${\chi }_{\text{spa}}$ is the term that takes into account the fact that from one SAR image to another the same ground resolution cell is imaged from two slightly different looking angles. The change change of the looking angle, in turn, leads to a shift between the range spectra of the two SAR images, and accordingly, it causes decorrelation since the range spectra of the two interfering SAR images are only partly overlapped. It can be shown that range spectra shift depend on the perpendicular baseline of the considered SAR data pair, and there is a limit value for the perpendicular baseline (known to as critical baseline) for which the two range spectra are completely non overlapped (i.e., the images are definitely uncorrelated one another) [29, 51].

• The term ${\chi }_{\text{dop}}$ takes into account of the so-called Doppler decorrelation effects due to the fact that SAR azimuth spectra are centered on a specific frequency (Doppler Centroid). When two SAR images with considerably different Doppler Centroid values interfere, a decorrelation noise contribution arises from the imperfect overlapping of the two related azimuth spectra. Hence, in the case that the two azimuth spectra are not overlapped at all, we have ${\chi }_{\text{dop}}=0$.

• ${\chi }_{\text{mis}}$ accounts for possible misregistration between two SAR images.

• ${\chi }_{\text{vol}}$ accounts for volumetric decorrelation effects [53].

The multilook operation, leading to the multilook phase $$\psi$$ in Eq. (25), reduces the level of noise corrupting interferograms, although this is paid in terms of a reduction of spatial resolution of interferograms. Multilook interferometric phase can be described by using a random quantity and, accordingly, it can be characterized via the knowledge of its probability density function. It has been shown in literature [40, 41, 54] that the probability density function (pdf) of an L-multilook interferometric phase (with L being the number of averaging samples in the averaging window used for the estimation of the statistical average operation involved in the calculation of Eq. (25)) can be given in terms of a Gauss hypergeometric function:

$$p_L\left(\psi \right)=\frac{\text{Γ}\left(L+\frac{1}{2}\right){(1-{\rho }^2)}^L\beta }{2\sqrt{\pi }\mathrm{\Gamma }\left(L\right){\left(1-{\beta }^2\right)}^{L+\frac{1}{2}}}+\frac{{(1-{\rho }^2)}^L}{2\pi }{}_{2}F{}_1\left[L,1;\frac{1}{2};{\beta }^2\right]$$E27

where $$\psi \in [-\pi ,\pi )$$, $$\beta =\rho cos(\psi -{\psi }_0)$$, $$\rho$$ represents the coherence, $${}_{2}F{}_1$$ denotes the Gauss hypergeometric function, $$\text{Γ}(\cdot )$$ is the gamma function, and $${\psi }_0$$ represents the expected “true” value of the interferometric phase. The peak of the distribution is located at $$\psi ={\psi }_0$$.

The pdf in Eq. (27) is sketched in Figure 2a for different values of L and in Figure 2b for different values of the $$\rho$$, with $${\psi }_0=0$$. By observing Figure 2a, it is clear that pdfs become narrower as the number of looks L increases (as expected). This finding is extremely important because it demonstrates that the interferometric phase may be thought to be corrupted by an additive noise random signal, namely, $$v$$, that has the same pdf as in Eq. (27) but with a zero-mean expected value, i.e., we may assume as valid the following additive model for the interferometric noise [54]: $$\psi ={\psi }_0+v$$. To further investigate about the statistics of multilook interferograms, we can observe that the validity of Eq. (27) is only restricted to the [–π,π) interval. However, this restriction does not apply when the phase signal is directly derived in the complex plane instead of the real plane. In the works of Lopez (2003) [55] and Lopez and Pottier (2007) [40], a comprehensive analytical derivation of the noise statistics in the complex plane is derived. Nonetheless, the Cramér–Rao bound for the standard deviation of multilook phase is given by [59]

$${\sigma }_v=\frac{1}{\sqrt{2L}}\frac{\sqrt{1-{\rho }^2}}{\rho }$$E28

that shows that standard deviation depends on the coherence $$\rho$$ and multilook factor L. Note that the phase standard deviation approaches the limit (Eq. (28)) asymptotically as the number of looks increases.

Figure 2.

Probability density function of the interferometric phase $$\psi$$ [rad]: (a) for different values of the number of looks L (1—black line, 2—red line, 5—blue line, 10—green line, 20—orange line), with $$\rho =0.7$$ ; (b) for different values of the correlation coefficient $$\rho$$ (0.1—black line, 0.2—red line, 0.5—blue line, 0.7—green line, 0.8—orange line) with L = 4.

3.2. Single-channel noise-filtering approaches for multilook interferograms

In order to mitigate the effects of decorrelation noise artifacts affecting SAR interferograms, several noise-filtering techniques, mostly working on single-channel data, have been proposed in literature over the years [42, 54, 56, 57]. As shown in previous Section 3.1, the statistics of multilook interferograms can be characterized via a probability density function expressible in closed form (Eq. (27)), and the noise standard deviation generally depends upon the coherence $$\rho$$ and the number of looks L [see also Eq. (28)]. Three different multilook interferograms, which are characterized by the same perpendicular baseline (of about 100 m), have been obtained by using three SAR sensors working at the different (C, X, and L) bands of the microwave region and are depicted in Figure 3. As it is evident from Figure 3, L-band interferograms are less affected by noise than the ones pertinent to C and/or X-band.

Figure 3.

Multilook interferogram computed by using different SAR data pairs: (a) July 11, 2011–August 16, 2011, X-band Cosmo-SkyMed (CSK); (b) September 15, 2004–October 29, 2004, C-band ASAR/ENVISAT; (c) July 30, 2007–September 14, 2007, L-band ALOS/PALSAR-1.

It should be emphasized that coherence and noise levels can also significantly change from one SAR system to another depending on the operational wavelength.Multilook processing has been proved to be effective for noise reduction, but this is paid in terms of a decrease of the image spatial resolution. Noise filtering constitutes a preliminary step before phase unwrapping. Indeed, the multilook operation usually involves an averaging on neighboring SAR pixels, hence reducing the spatial resolution of the interferograms. Several algorithms have been proposed in literature. The most commonly used noise filter is the boxcar filter applied in the complex plane. Another frequently used option is provided by the Golstein’s frequency-domain algorithm [42], which is an empirical approach proposed for geophysical applications. In this case, a complex interferogram (amplitude and phase) is segmented into overlapping rectangular patches and for each patch the relevant power spectrum Z is computed.

Figure 4.

Multilook interferogram relevant to the SAR data pair July 11, 2011–August 16, 2011, X-band Cosmo-SkyMed (CSK): (a) Original, (b–d) Goldstein filtering, with (b) $$\alpha =0.25$$, (c) $$\alpha =0.5$$, and (d) $$\alpha =1.0$$.

The response of the Golstein’s filter is then computed from the power spectrum as follows:

$$H(\xi ,\eta )={\left|Z(\xi ,\eta )\right|}^{\alpha }$$E29

where $$\xi$$ and $$\eta$$ denote the relevant spectral variables, respectively. Notice that the filtering effect vanishes when $$\alpha =0$$ ; conversely, the filtering effect is more pronounced as the parameter $$\alpha$$ approaches unity. We show in Figure 4 the result of the application of the Goldstein’s filter to a multilook interferogram, relevant to the Mt. Etna (Italy) volcano, obtained by using the Cosmo-Skymed sensor of the Italian Space Agency (ASI). Specifically, different values of the filtering parameter $$\alpha$$ have been considered in Figure 4. The limited effectiveness of the filtering capabilities of Goldstein approach is evident from the result depicted in Figure 4. A modification of the Goldstein filter that relies on an adaptive choice of the filtering factor $$\alpha$$ (which depends on the spatial coherence $$\rho$$) has also (more recently) been proposed by Baran in 2003 [58]. Other filters, such as the median filter [59] and the two-dimensional Gaussian filter, are also used to reduce noise while performing multilooking operations. It is worth noting that boxcar and Goldstein filters do not adapt to the direction of the fringes because these filters are operated in a square window. In order to overcome such a limitation, Lee et al. 1998 [54] then proposed a self-adaptive filter based on local gradient slope estimation that is able to improve noise-filtering performance by exploiting directional characteristics of an InSAR interferogram. Several adaptations and relevant improvements of the Lee filter have subsequently proposed in literature over the recent years [56, 57], most of them based on the exploitation of the intrinsic directional behavior of InSAR interferograms. In fact, compared with the fringe phase and gradient, the fringe direction variation is gently, thus making it possible to use fringe direction to guide interferogram filtering.

3.3. The multichannel noise-filtering (MCh-NF) algorithm

The noise-filtering methods discussed in the previous Section have historically been developed to analyse and filter out the noise affecting single interferograms, separately, thus without taking into account their mutual temporal relationships. A multichannel noise-filtering problem arises when a stack of SAR interferograms need to be jointly exploited. In this case, it is profitable to develop/use noise-filtering approaches that not only exploit spatial/frequency information but can also take into account temporal relationships among available multichannel interferograms, in order to identify and filter out the noise affecting the whole interferometric data stack in the more reliable way as possible. A specific multichannel noise-filtering (MCh-NF) method [43], which is based on using a stack of time-redundant multilook interferograms, is described in this Section. The MCh-NF method is here described by adopting the same rigorous formalism and terminology used for the topological description of multichannel phase unwrapping problem presented in Section 2. According to the adopted symbolism, let us consider Q SAR images and let M be the number of multilook interferograms characterized by small perpendicular and temporal baselines.

The resulting interferometric data stack of the M (wrapped) small-baseline multilook interferograms can be expressed as $${\tilde{\Psi }}^{\text{Τ}}=[{\tilde{\psi }}_{\text{1}}^{\text{T}},\dots ,{\tilde{\psi }}_P^{\text{T}}]$$ ; thus, the M-dimensional vector $${\tilde{\psi }}_p^{\text{T}}$$ described the (vectorized) multichannel interferometic-phase pertinent to the pth pixel, with $$p\in \{1,\dots ,P\}$$ and P denoting the number of coherent pixels common to all interferograms. In particular, $$\tilde{\Psi }$$ can be expressed in terms of discrete gradient $${\Pi }_{\text{A}}$$, according to Eq. (13), as:

$${\tilde{\Psi }}^{\text{Τ}}=W\left({\Pi }_{\text{A}}\Phi +D\right)$$E30

wherein $$\Phi$$ represents the (unknown) phases associated with the available SAR images, and the matrix $$D$$ describes the additive noise-term that corrupts the stack of interferograms. The noise term should be estimated and properly filter out from the generated interferograms. As discussed in Section 3.2, the term $$D$$ arises since both a multilook operation and a noise-filtering procedure are typically applied to each single interferogram, separately. Both these operations are independently carried out on each single interferogram; hence, they are not necessarily time consistent. The fact that the interferograms are not fully time consistent can be formally expressed, according to Eq. (14), in terms of discrete curl $${\Omega }_{\text{A}}^{\text{T}}$$, in the form:

$$\mathrm{W}\left({\Omega }_A^T{\tilde{\Psi }}^T\right)=W\left({\Omega }_A^{T}D\right)\ne 0$$E31

which represents the topological generalization of the phase-triangularity condition exploited by the SqueeSAR technique [44]. Therefore, the multichannel noise-filtering (MCh-NF) approach suitably addresses the temporal inconsistencies inherent in the time-redundant multilook interferograms, which can be mathematically described in terms of the (modulo-2π) cyclic inconsistency of the multichannel interferometric phases [see Eq. (31)]. More specifically, MCh-NF is based on the solution of the a nonlinear optimization problem, as detailed in the following. First, $$\forall p\in \{1,\dots ,P\}$$, the Q-dimensional vector (Q is the number of SAR acquisitions) representing the (unknown) wrapped phases $${\tilde{\mathrm{\Phi }}}^p=W({\mathrm{\Phi }}^p)$$ is estimated as follows:

$${\widehat{\tilde{\mathrm{\Phi }}}}^p=\underset{{\tilde{\mathrm{\Phi }}}^p\in {ℝ}^Q}{\arg \max }\left|{\overline{\zeta }}_p\circ e^{j({\tilde{\mathrm{\Psi }}}_p^T-W({\Pi }_{\text{A}}{\tilde{\mathrm{\Phi }}}^p\text{)})}\right|$$E32

where $$j=\sqrt{-1}$$ denotes the imaginary unit, P denotes the number of coherent pixels to which the noise-filtering procedure is applied, the symbol $$\circ$$ represents the Hadamard product, and $$\overline{{\zeta }_p}={[{\overline{\zeta }}_p^1,\cdots ,{\overline{\zeta }}_p^M]}^T$$ is an M-dimensional normalized weighting vector representing our confidence on the quality of the exploited M (small-baseline) interferometric phases pertinent to the pth pixel, with

$${\overline{\zeta }}_p^m=\frac{{\zeta }_p^m}{{\displaystyle \sum _{h=1}^M{\zeta }_p^h}}$$E33

wherein the generic elements $${\zeta }_p^m$$ can be related to the spatial coherence as detailed after. Subsequently, these estimated vectors $${\widehat{\tilde{\mathrm{\Phi }}}}^p$$ are used to reconstruct a new (noise-filtered) stack of time-consistent interferograms $${\widehat{\tilde{\Psi }}}^{\text{Τ}}=[{\widehat{\tilde{\psi }}}_{\text{1}}^{\text{T}},\dots ,{\widehat{\tilde{\psi }}}_P^{\text{T}}]$$, where $${\widehat{\tilde{\mathrm{\Psi }}}}_p^T=W({\Pi }_{\text{A}}{\widehat{\tilde{\mathrm{\Phi }}}}^p)$$ and $$p\in \{1,\dots ,P\}$$. We emphasize that, according to Eq. (32), the MCh-NF technique is based on searching for the (unknown) wrapped-phase vector $${\tilde{\mathrm{\Phi }}}^p\in {ℝ}^Q$$ that minimizes the (weighted) circular variance of the random (phase) vector representative of the phase difference, $${\tilde{\mathrm{\Psi }}}_p^T-W({\Pi }_{\text{A}}{\tilde{\mathrm{\Phi }}}^p\text{)}$$, between the “original” and the “reconstructed” interferograms.

The evaluation of the weights for the optimization problem in Eq. (32) is now addressed. Let $${\stackrel{⌣}{\Theta }}^{}=[{\stackrel{⌣}{\mathrm{\Theta }}}_{i,j}^{}]$$ be a matrix description for a generic 2-D phase map, whose corresponding vectorized representation is provided by the P-dimensional vector $$\mathrm{\Theta }$$. Each pixel of the phase map is identified by discrete range and azimuth coordinates, denoted by i and j, respectively. Accordingly, each pair (i, j) is uniquely associated with a monodimension index $$p\in \{1,\dots ,P\}$$ identifying an element of the vector $$\Theta$$. The spatial coherence relevant to $$\mathrm{\Theta }$$ (i.e., $$\stackrel{⌣}{\mathrm{\Theta }}$$) evaluated around the pixel (i, j) (associated with the index p) is defined as

$${\zeta }_p\left(\mathrm{\Theta }\right)=\frac{1}{\left(2L_R+1\right)\left(2L_A+1\right)}\left|{\displaystyle \sum _{l=-L_R}^{L_R}{\displaystyle \sum _{h=-L_A}^{L_A}\exp \left[j{\stackrel{⌣}{\mathrm{\Theta }}}_{i+l,j+h}\right]}}\right|$$E34

where $$2L_R+1$$ and $$2L_A+1$$ are the number of azimuth and range pixels within the used boxcar averaging window, which is centred around the generic pixel identified by the discrete range and azimuth coordinates, i and j, respectively. In particular, the mth weight $${\zeta }_p^m$$ is expressed, according to Eq. (34), in terms of spatial coherence relevant to the (vectorized) interferograms $${\tilde{\psi }}^m$$ and evaluated around the pixel associated with the index p, in the functional form $${\zeta }_p^m={\zeta }_p^{}({\tilde{\psi }}^m)$$, $$\forall m\in \{1,\dots ,M\}$$. Therefore, $${\zeta }_p={[{\zeta }_p^1,\cdots ,{\zeta }_p^M]}^{\text{T}}$$ can be evaluated in terms of the spatial coherence directly from the stack of M multilook interferograms $$\tilde{\Psi }=[{\tilde{\psi }}^1,\dots ,{\tilde{\psi }}^M]$$. As experimentally demonstrated in ref. [43], the “reconstructed” interferograms with MCh-NF are significantly less affected by noise than the original ones. However, a group of the reconstructed interferograms, although limited, can exhibit spatial coherence values smaller than the ones relevant to the corresponding original interferograms, thus implying that a partial corruption of the spatial coherence occurs during the minimization procedure. In particular, it happens in correspondence to interferograms that were originally significantly coherent, and this is due to the fact that the strategy in Eq. (32) tends to “inject” coherence from very coherent to incoherent interferograms, by exploiting the time redundancy of the small baseline data pairs. Accordingly, in order to also preserve the spatial coherence of the very coherent interferograms, a simple nonlinear combination between the original and the reconstructed interferograms is carried out, thus further increasing the phase quality of the whole set of M reconstructed interferograms. In particular, the two sets of interferograms are combined through the following (wrapped) weighted averaging operation:

$${\overline{\psi }}^m=\arctan \left(\frac{{\zeta }^m\circ \sin \left({\tilde{\psi }}^m\right)+{\widehat{\zeta }}^m\circ \sin \left({\widehat{\tilde{\psi }}}^m\right)}{{\zeta }^m\circ \cos \left({\tilde{\psi }}^m\right)+{\widehat{\zeta }}^m\circ \cos \left({\widehat{\tilde{\psi }}}^m\right)}\right)\forall m\in \left\{1,\dots ,M\right\}$$E35

where the symbol $$\circ$$ represents the Hadamard product, and $${\zeta }^m$$ and $${\widehat{\zeta }}^m$$ are two P-dimensional vectors. In particular, $${\zeta }^m={[{\zeta }_1^m,\cdots ,{\zeta }_P^m]}^{\text{T}}={[{\zeta }_1^{}({\tilde{\psi }}^m),\cdots ,{\zeta }_P^{}({\tilde{\psi }}^m)]}^{\text{T}}$$ is expressed in terms of the spatial coherence relevant to the original multilook interferogram $${\tilde{\psi }}^m$$. Similarly, $${\widehat{\zeta }}^m={[{\widehat{\zeta }}_1^m,\cdots ,{\widehat{\zeta }}_P^m]}^{\text{T}}={[{\zeta }_1^{}({\widehat{\tilde{\psi }}}^m),\cdots ,{\zeta }_P^{}({\widehat{\tilde{\psi }}}^m)]}^{\text{T}}$$ is expressed in terms of the spatial coherence relevant to the reconstructed multilook interferogram $${\widehat{\tilde{\psi }}}^m$$. The block diagram of the MCh-NF algorithm is depicted in Figure 5. A pertinent pseudo-code for computing the filtered interferometric data stack is also presented (Figure 6).

Note that the exploitation of “conventional” small baseline multilook interferograms is the distinctive characteristic of MCh-NF approach with respect to other previous solutions, such as the SqueeSAR [44] and Phase Linking [62] methods and other recently proposed multitemporal-filtering techniques [60, 61] based on constraining the analysis to distributed scatterers [29], which are identified through a pixel-by-pixel selection procedure performed at the full resolution complex SAR image spatial grid. Such a selection permits to rely on the distributed scattering hypothesis, under which the probability density function (pdf) of the complex-valued SAR image may be regarded as being a zero-mean multivariate circular normal distribution, and an appropriate maximum likelihood (ML) estimation step of the filtered phase values associated to each SAR acquisition is implemented. On the contrary, the presented MCh-NF approach focuses on conventional multilook interferograms generated without any a priori pixel selection step. Accordingly, in this case, it is not possible to rely on the validity of the above-mentioned distributed scattering hypothesis. Therefore, both the phase linking [62] and the phase triangulation of the SqueeSAR [44] algorithms require a preliminary identification of the statistically homogeneous pixels (SHPs) on the full-resolution range-azimuth grid. In particular, in ref. [44], the selection strategy of these pixels is based on the application of the Kolmogorov–Smirnov test to carefully select a homogeneous statistical population. Clearly, this requires working at the full resolution spatial scale and implies the analysis of the amplitude values of the complex SAR image pixels. A more detailed comparison among the presented MCh-NF method and the ones provided in ref. [44] can be found in ref. [43].

Figure 5.

MCh-NF block diagram

4. Experimental results

We present in this section some results we obtained by processing a data set consisting of 39 SAR images (Track 308, Frame 2754), collected by the ENVISAT sensor between December 2002 and August 2010 over the Abruzzi region (Italy). The test-site area includes the city of L’Aquila and its surroundings, which were struck on 6 April 2009, by an Mw 6.3 earthquake that caused more than three hundred fatalities, thousands of evacuees, as well as severe industrial and residential building damages. Starting from the available SAR images, we retrieved a stack of 338 small baseline differential SAR interferograms with a maximum perpendicular baseline of 400 m and a maximum time span of 2000 days [43]. The interferograms have been computed by performing a complex multilook operation with 4 and 20 looks in the range and azimuth directions, respectively. For the interferogram generation, we used precise satellite orbit information and a three-arcsecond shuttle radar topography mission (SRTM) digital elevation model (DEM) of the region to remove the topographic phase contributions. Finally, the multilook interferograms have been prefiltered by applying the Goldstein’s filter [42].

Figure 6.

Pseudo code of the MCh‐NF algorithm.

Figure 7.

Comparison between the original (left column) and noise-filtered (right column) multilook interferograms relevant to the area of the Abruzzi region (Italy). (a–c) October 30, 2005, to November 8, 2009, August 21, 2005, to June 6, 2010, and August 1, 2004, to April 12, 2009, interferograms, characterized by perpendicular baseline values of 192, 145, and 395 m, respectively. (d–f) Noise-filtered multilook interferograms corresponding to the ones in panels a–c, respectively.

To investigate the performance of the presented noise-filtering approach, we applied the nonlinear minimization procedure in Eq. (32) to the stack of the generated (original) multilook small baseline interferograms. As a result, we retrieved a new set of noise-filtered interferograms that are characterized by generally improved coherence values. This is clearly visible in Figure 7a–f, where we compare three unfiltered (left side) interferograms with the corresponding (right side) noise-filtered interferograms. It is evident how the fringes due to the earthquake are well recovered. Such interferometric data stacks can then be used for the generation of surface deformation time series using the small baseline subset (SBAS) [6] processing chain, whose parallel version (P-SBAS) has been proposed in refs. [13, 14, 35, 36]. This is matter for the analysis presented in the next subsection.

Figure 8.

Block diagram of the advanced EMCF-SBAS processing chain.

4.1. The use of MCh-NF algorithm within MCh-PhU framework

We present in this subsection how the MCh-NF algorithm can be efficiently used within the SBAS processing chain, where phase unwrapping procedures are implemented through the MCh-PhU technique known as extended minimum cost flow (EMCF), also discussed in the first part of the chapter. Figure 8 shows the diagram block of the advanced EMCF-based SBAS processing chain [26, 35, 43], which integrates the conventional SBAS codes, exploiting the EMCF MCh-PhU procedure to perform phase unwrapping operations, with the presented MCh-NF noise-filtering technique. In addition, an effective procedure for the selection of time-redundant interferograms is also included; interested readers can find additional details in ref. [43]. To provide an example of the potential of the advanced processing chain incorporating both the discussed MCh-PhU and MC-NF techniques, we here focus on the area of Yellowstone caldera, representing one of the largest and most active volcanic systems in the world. We analyze the temporal evolution of the surface deformation occurring in this area by applying the implemented EMCF-SBAS processing chain to a set of 22 ENVISAT images (Track 41, Frame 2709), acquired from May 2005 to September 2010, from which we have retrieved a corresponding set of 122 small baseline interferograms [43]. As in the previous case study (relevant to Abruzzi area), the prescribed maximum values of 400 m and 2000 days for perpendicular and temporal baseline, respectively, have also been considered. The retrieved mean deformation velocity map depicted in Figure 9, which has been obtained by applying the processing chain (MCh-NF + EMCF-SBAS) of Figure 8, allows us to recognize the complex deformation scenario affecting the Yellowstone Caldera region and its surroundings, where uplift effects and broad subsidence patterns characterize the detected deformation field. In addition, the deformation time series relevant to four pixels, whose locations are identified by the black stars labeled as P1, P2, P3, and P4, are also illustrated in Figure 9.

Figure 9.

(a) Mean deformation velocity map (in color) of Yellowstone Caldera, computed in coherent pixels only and superimposed on the SAR amplitude image (gray-scale representation) of the zone, retrieved by applying the advanced EMCF-SBAS processing chain. The black square marks the location of the reference SAR pixel. (b–e) Deformation time series relevant to the four pixels identified via black stars in Fig. 8(a).

5. Conclusion

Within the context of SAR interferometry, two main issues related to the multichannel phase unwrapping and noise filtering for interferometric data stacks processing have been addressed. First, a rigorous gradient-based formulation for the multichannel phase unwrapping (MCh-PhU) problem has been systematically established, thus providing a topological characterization of the problem within the purview of the theoretical foundation of the discrete calculus. Then the innovative MCh-NF procedure for the noise filtering of time-redundant multichannel multilook interferograms has been properly presented within the considered topological framework, by adopting a consistent formalism. Finally, some experimental results obtained with real data have been shown, thus demonstrating the effectiveness of our approaches and their relevance for geospatial phenomena analysis and understanding.

Acknowledgments

This work was supported by the Italian Ministry of University and Research (MIUR) under the project “Progetto Bandiera RITMARE.” We would like to thank the European Space Agency for providing the ENVISAT ASAR data and the University of Delft, Delft, The Netherlands, for the related precise orbits. We would also like to thank Italian Space Agency (ASI), which has provided us the Cosmo-SkyMed SAR images under the framework of the European Union’s Seventh Program for research, technological development, and demonstration MED-SUV project (grant no. 308665). Finally, the authors thank the Japanese Space Agency (JAXA), which has provided the used ALOS-1 data through the project entitled “Advanced Interferometric SAR Techniques for Earth Observation at L-band” (ID project 1149) in the framework of “The 4-th ALOS Research Announcement for ALOS-2” call.