I. INTRODUCTION
Synthetic aperture radar interferometric systems (In-SAR) allow the estimation of digital elevation models (DEM) of the observed scene [Reference Massonet and Feigl1, Reference Itoh2], starting from two SAR images of the same scene acquired with slightly different view angles along two different trajectories, spatially separated by a distance called baseline (single-baseline configuration).
The phase image directly extracted from a complex SAR interferogram is wrapped since it represents a measure, modulo-2π mapped, of the actual phase difference between two co-registered SAR images of the same scene; furthermore, the interferogram is subject to phase noise. In order to build the elevation map of the scene, the phase image must be first unwrapped so as to form a continuous phase map [Reference Itoh2]. In this perspective, the development of phase unwrapping (PU) algorithms has attracted growing interest during the last years and has rapidly become one of the key topics for SAR interferometry applications [Reference Costantini3–Reference Yu, Li and Bao11]. The main difficulties involved in the PU problem are intrinsic nonlinearity of the SAR interferogram and the nature of the phase noise; consequently, the problem can be regarded as a classic nonlinear estimation problem. Besides, the PU problem becomes an ill-posed problem when the absolute value of phase difference between neighboring pixels exceeds (Itoh's condition is violated) [Reference Pascazio and Schirinzi5, Reference Budillon, Ferraiuolo, Pascazio and Schirinzi6] due to the presence of abrupt terrain discontinuities and/or interferometric noise. Moreover, many PU techniques adopt a noise-filtering stage to reduce phase noise, followed by the PU algorithm itself. It should be apparent that this approach is largely suboptimal because a portion of the information contained in the noisy pixels vanishes during the noise-filtering stage and cannot be recovered by subsequent processing.
In the last years, multichannel techniques [Reference Pascazio and Schirinzi5–Reference Dai and Zha12], which require multiple SAR acquisitions, have gained growing attention because of their ability to overcome, or at least mitigate, the ambiguity problem and to improve PU accuracy by exploiting baseline and/or frequency diversity. Some of them are based on statistical estimation techniques, such as maximum likelihood (ML) [Reference Pascazio and Schirinzi5–Reference Eineder and Adam8] and maximum a posteriori (MAP) Bayesian [Reference Ferraiuolo, Pascazio and Schirinzi9, Reference Shabou, Baselice and Ferraioli10] estimation, and have proved their effectiveness even in the presence of strong height discontinuities and for high sloped and low-coherence areas.
Recently, the usefulness of the Kalman filter [Reference Loffeld, Nies, Knedlik and Wang4] for PU and noise reduction has been demonstrated. The PU algorithm presented in this paper is based on the extended Kalman filter that performs simultaneously PU, noise reduction, and data fusion. Our solution differs from that proposed in [Reference Loffeld, Nies, Knedlik and Wang4] because it makes use of a maximum likelihood estimator (MLE), instead of local frequency estimators, and it is suitable for multichannel configuration. Furthermore, a new state model [Reference Chirico and Schirinzi13] incorporating the information of the pixels in a small window around the pixel of interest has been introduced. The main consequence of the new state-space model is to lessen the overall impact of the pixels affected by unwrapping errors on the estimate of the pixels yet to be unwrapped; in such a way, the propagation effect of the unwrapping errors has been mitigated. In addition, our algorithm allows the estimation of the height profile, producing at the same time an estimate of the height slope along the range and azimuth directions. One of the key points of our algorithm is its formulation in terms of fixed-interval smoothing that makes it appropriate for the batch nature of SAR interferogram data. The paper is organized as follows: in Section II, the mathematical statement of the PU problem for the multichannel case is introduced, while in Section III, the proposed PU algorithm is described. Then, Section IV reports the results obtained in the simulations and Section V is devoted to the conclusions.
II. PU PROBLEM STATEMENT
Given C ≥ 1 complex SAR interferograms of size M × N obtained in multifrequency and/or multibaseline configuration and corresponding to the same imaged scene, let I = {(m, n)|m = 1,2, … , M, n = 1,2, … , N} be a rectangular grid of points (or pixels) in each of the interferograms. With reference to the cth acquisition channel, where c assumes its value in the range 1,2, … , C, the complex value of the SAR interferogram at pixel (m, n) ∈ I has the following expression:
$$\eqalign{u_c \lpar m\comma \; n\rpar &= a_c \lpar m\comma \; n\rpar \exp\! \lcub j\tilde \varphi _c \lpar m\comma \; n\rpar \rcub \cr & =a_c \lpar m\comma \; n\rpar \lcub \cos\! \lpar \varphi _c \lpar m\comma \; n\rpar \rpar +j\sin \!\lpar \varphi _c \lpar m\comma \; n\rpar \rpar \rcub \cr & = a_c \lpar m\comma \; n\rpar \lcub \cos \!\lpar \alpha _c h\lpar m\comma \; n\rpar \rpar +j\sin \!\lpar \alpha _c h\lpar m\comma \; n\rpar \rpar \rcub \comma \; }$$
where we exploited the 2π periodicity of the sin(·) and cos(·) functions. In the previous expression, 
$\tilde \varphi _c \lpar m\comma \; n\rpar $ is the measured interferometric phase, φc(m, n) is the true unambiguous interferometric phase, a c(m, n) is the interferometric amplitude, h(m, n) is the height profile and αc = 4πB c⊥/(λcR 0sinθ), λc being the wavelength, B c⊥ the orthogonal baseline measured along the direction orthogonal to the range, θ the view angle and R 0 the distance from the centre of the scene. Since 
$\tilde \varphi _c \lpar m\comma \; n\rpar $ is ambiguous to within-integer multiples of 2π, it is said to be wrapped. Now let W{·} indicate the nonlinear wrapping operator given by:
$$\eqalign {\tilde \varphi _c \lpar m\comma \; n\rpar &=W\lcub \varphi _c \lpar m\comma \; n\rpar \rcub \cr & ={\rm atan2}\left({{\mathop{\rm Im}\nolimits} \lpar e^{j\varphi _c \lpar m\comma n\rpar } \rpar \comma \; {\mathop{\rm Re}\nolimits} \lpar e^{j\varphi _c \lpar m\comma n\rpar } \rpar } \right)\in \left]{ - \pi\comma \; \pi } \right]\comma \; }$$
which wraps the absolute phase φc(m, n) into the principal interval ] − π,π] and outputs the wrapped phase 
$\tilde \varphi _c \lpar m\comma \; n\rpar $. Owing to decorrelation effects and to the presence of speckle [Reference Bamler and Hartl14], the observed interferometric phase is also affected by noise. The multichannel PU problem consists of inferring the value of the ground height h(m, n) over the entire imaged scene, using the M × N × C noisy observations of the modulo 2π mapped interferometric phase (or wrapped phase) 
$\tilde \varphi_c$:
$$\tilde \varphi _c \lpar m\comma \; n\rpar =W\lcub \varphi _c \lpar m\comma \; n\rpar +v_c \lpar m\comma \; n\rpar \rcub \comma \;$$where v c(m, n) describes the phase noise that enters the system through the nonlinear wrapping operator, which is assumed to have the well-known probability density function (pdf) depending on the interferogram coherence [Reference Bamler and Hartl14]. The inversion of (3), which is the core element of the PU process, is an ill-posed problem where many solutions may correspond to a given interferogram. To overcome the obstacle, some additional a priori information must be introduced; this is accomplished hereinafter through the introduction of a state model describing spatial variation of the phase within the SAR interferogram.
III. MULTICHANNEL PU ALGORITHM
In this section, the PU problem is reformulated in state-space form; to this end, we provide a state-space representation of spatial variation of phase signal and we also provide a nonlinear observation model that expresses the relation between the true absolute phase and the observable (noisy corrupted) wrapped phase.
A) State-space model
Let x(m, n) = [h(m, n) h x(m, n) h y(m, n)]T be the state vector at pixel (m, n), where h(m, n) is the height profile, h x(m, n) is the height slope along the x direction, and h y(m, n) is the height slope along the y direction. The PU problem is inherently an ill-posed problem because different solutions may correspond to a single observed SAR interferogram; this effect becomes even more severe due to the presence of height discontinuities and phase noise. Consequently, the PU problem has to be constrained introducing some smoothness or regularity assumptions on the surface to be reconstructed. The Kalman filter fulfills this requirement offering an effective way to express the relationship between a pixel and its neighborhood in terms of state-space models. A simple way to do this is to consider the height at a given pixel as the contribution of the height and the height slope of the L unwrapped pixels around the pixel of interest; thus, the state-space model adopted in this paper, which can be applied both in the forward and backward recursion of the extended Kalman filter, is:
$$\eqalign{{\bi x}\lpar m\comma \; n\rpar =&\displaystyle{1 \over L}\sum\limits_{\lpar k\comma j\rpar \in W_d \cap I} {\left[{\matrix{ 1 & {m - k} & {n - j} \cr 0 & 1 & 0 \cr 0 & 0 & 1 \cr } } \right]} \cr &\times{\bi x}\lpar k\comma \; j\rpar g\lpar k\comma \; j\rpar +{\bi v}\lpar m\comma \; n\rpar \comma \; }$$where I is the rectangular grid previously defined, W d denotes the pixels within a d × d window around pixel (m, n) and the function g(k, j) equals one if the corresponding pixel has already been unwrapped and zero otherwise. The noise vector v accounts for the process noise and is modeled as a zero mean vector with covariance matrix Q. The predicted covariance matrix P− to be used in the Kalman filter can be calculated from (4):
$$\eqalign{{\bi P}^ - \lpar m\comma \; n\rpar &=\displaystyle{1 \over {L^2 }}\sum\limits_{\lpar k\comma \,j\rpar \in W_d \cap I} {\left[{\matrix{ 1 & {m - k} & {n - j} \cr 0 & 1 & 0 \cr 0 & 0 & 1 \cr } } \right]} \cr & \quad\times{\bi P}\lpar k\comma \; j\rpar \left[{\matrix{ 1 & {m - k} & {n - j} \cr 0 & 1 & 0 \cr 0 & 0 & 1 \cr } } \right]^T g\lpar k\comma \; j\rpar +{\bi Q}\lpar m\comma \; n\rpar .}$$B) Multichannel observation model
The measurement equation (2) presents a high degree of nonlinearity; in addition, the phase noise is subject to a wrap-around effect that circularly shifts the phase noise samples in the interval ]−π,π]. For the given reasons, (2) is not well suited as a basis for the design of an extended Kalman filter. Instead, we adopt the nonlinear observation model proposed in [Reference Loffeld, Nies, Knedlik and Wang4] which renders the PU problem tractable by the extended Kalman filter, providing a good filter performance as well. In order to formulate the measurement model for the multichannel configuration, we first have to express the measurement equation at pixel (m, n) for the cth channel. With no loss of generality, and without altering the phase statistic of the observations, the sought after nonlinear observation model is achieved normalizing the amplitude of the complex SAR interferogram and then forming the vector made by the normalized in-phase and quadrature components of the complex interferogram (1), as shown:
$$\eqalign{{\bi z}_c \lpar m\comma \; n\rpar =&\left[{\matrix{ {{\mathop{\rm Re}\nolimits} \left\{{\displaystyle{{u_c \lpar m\comma \; n\rpar } \over {a_c \lpar m\comma \; n\rpar }}} \right\}} \cr {{\mathop{\rm Im}\nolimits} \left\{{\displaystyle{{u_c \lpar m\comma \; n\rpar } \over {a_c \lpar m\comma \; n\rpar }}} \right\}} \cr } } \right]=\left[{\matrix{ {\cos \!\lpar \alpha _c h\lpar m\comma \; n\rpar \rpar } \cr {\sin \!\lpar \alpha _c h\lpar m\comma \; n\rpar \rpar } \cr } } \right] \cr & +{\bi w}_c \lpar m\comma \; n\rpar {\rm\comma \; }\quad {\rm }c=1\comma \; 2\comma \; \ldots\comma \; C\comma \; }$$where zc(m, n) is the measurement vector at pixel (m, n), and the measurement noise samples wc(m, n) are modeled as independent random vectors with zero mean and covariance matrix Rc(m, n). The linearized measurement matrix for each single channel is given by:
$$\eqalign{{\bi H}_c \lpar m\comma \; n\rpar &=\left. {\displaystyle{d \over {d{\bi x}}}\left[{\matrix{ {{\mathop{\rm Re}\nolimits} \left\{{\displaystyle{{u_c \lpar m\comma \; n\rpar } \over {a_c \lpar m\comma \; n\rpar }}} \right\}} \cr {{\mathop{\rm Im}\nolimits} \left\{{\displaystyle{{u_c \lpar m\comma \; n\rpar } \over {a_c \lpar m\comma \; n\rpar }}} \right\}} \cr } } \right]} \right\vert _{{\bf \hat x}^. \lpar m\comma n\rpar } \cr &=\alpha _c \left[{\matrix{ { - \sin \!\lpar \alpha _c \hat h^ - \lpar m\comma \; n\rpar \rpar } & 0 & 0 \cr {\cos \!\lpar \alpha _c \hat h^ - \lpar m\comma \; n\rpar \rpar } & 0 & 0 \cr } } \right]{\rm\comma \; } \cr c&=1\comma \; 2\comma \; \ldots\comma \; C\comma \; }$$
where 
$\hat h - \lpar m\comma \; n\rpar $ is the predicted ground height computed using equation (4). In principle, due to the periodic nature of the sine and cosine functions in (6), the sole information coming from the in phase and quadrature components of the complex interferogram together with the state model in (4), does not suffice to completely resolve the ambiguity problem except when Itoh's condition is only slightly violated or when it is not violated at all. In order to overcome this limitation, the approach proposed in this paper tries to combine the advantages of the Bayesian estimation, guaranteed by the Kalman filter, to the already proven effectiveness of the maximum likelihood height estimator (MLHE) described in [Reference Pascazio and Schirinzi5] for resolving the phase ambiguity due to highly sloped height profiles eventually characterized by strong discontinuities. At each pixel (m, n), the ML height estimate is the value of h that maximizes the multichannel likelihood function, and can therefore be expressed as:
$$\hat h_{ML} \lpar m\comma \; n\rpar =\mathop {\arg \max }\limits_h \prod\limits_{i=1}^C {\Lambda \left({\varphi _i \lpar m\comma \; n\rpar \vert h\lpar m\comma \; n\rpar } \right)}\comma \;$$where h(m, n) is the height profile to be estimated and
$$\eqalign{&\Lambda \lpar \varphi _i \lpar m\comma \; n\rpar \vert h\lpar m\comma \; n\rpar \rpar =\displaystyle{{1 - \vert \gamma \vert ^2 } \over {2\pi \lpar 1 - \vert \gamma \vert ^2 \cos ^2 \lpar \varphi _i - \alpha _i h\rpar \rpar }} \cr & \times\left\{{1+\displaystyle{{\vert \gamma \vert \cos \lpar \varphi _i - \alpha _i h\rpar \cos ^{ - 1} \lsqb - \vert \gamma \vert \cos \lpar \varphi _i - \alpha _i h\rpar \rsqb } \over {\lsqb 1 - \vert \gamma \vert ^2 \cos ^2 \lpar \varphi _i - \alpha _i h\rpar \rsqb ^{1/2} }}} \right\}}$$is the pdf [Reference Bamler and Hartl14] (i.e. the single-channel likelihood function for the ith channel) of the measured wrapped φi given the actual phase αih. In the expression of the single-channel likelihood function (9), γ is the coherence coefficient whose value is determined by the decorrelation effects responsible for the differences between the couple of SAR images that originated the interferogram.
In our approach, the ML estimate (8) is considered as a supplementary pseudo-measurement to be used in the Kalman filter update step. Equations (6) and (8) describe the set of all of the available measurements for the nonlinear multisensor data fusion problem defined by the expressions in (4) and (7). Now, we move to the expression of the global observation model for the multibaseline configuration, in which we bring together all of the available observations coming from the C acquisition channels and the ML height estimate (the latter eventually filtered by means of a median filter). The proposed approach to deal with the multisensor data fusion problem is to combine the observations into a single composite observation model. To this end, using (6), (7), and (8), we introduce the (2C + 1) × 1 composite observation vector:
$${\bi z}\lpar m\comma \; n\rpar =\lsqb \hat h_{{\rm ML}} \lpar m\comma \; n\rpar\;{\bi z}_1^T \lpar m\comma \; n\rpar {\bi z}_2^T \lpar m\comma \; n\rpar \cdots {\bi z}_C^T \lpar m\comma \; n\rpar \rsqb ^T$$and the composite linearized observation matrix of size (2C + 1) × 3:
$$\eqalign{&{\bi H}\lpar m\comma \; n\rpar \cr & \quad=\lsqb {\bi H}_{ML}^T \lpar m\comma \; n\rpar \quad {\bi H}_1^T \lpar m\comma \; n\rpar \quad {\bi H}_2^T \lpar m\comma \; n\rpar \quad \cdots \quad {\bi H}_C^T \lpar m\comma \; n\rpar \rsqb ^T\comma \; }$$where HML(m, n) = [1 0 0]. Finally, the composite noise vector w of size (2C + 1) × 1 and its covariance matrix R of size (2C + 1) × (2C + 1) are given respectively by:
$${\bi w}\lpar m\comma \; n\rpar =\lsqb {\bi w}_{ML} \lpar m\comma \; n\rpar \quad {\bi w}_1^T \lpar m\comma \; n\rpar \quad {\bi w}_2^T \lpar m\comma \; n\rpar \quad {\bi w}_C^T \lpar m\comma \; n\rpar \rsqb ^T$$and
$$\eqalign{&{\bi R}\lpar m\comma \; n\rpar =E\lsqb {\bi w}\lpar m\comma \; n\rpar {\bi w}^T \lpar m\comma \; n\rpar \rsqb \cr & =\left[{\matrix{ {{\bi R}_{ML} \lpar m\comma \; n\rpar } & {{\bf 0}_{1 \times 2} } & {{\bf 0}_{1 \times 2} } & {{\bf 0}_{1 \times 2} } & {{\bf 0}_{1 \times 2} } \cr {{\bf 0}_{2 \times 1} } & {{\bi R}_1 \lpar m\comma \; n\rpar } & {{\bf 0}_{2 \times 2} } & {{\bf 0}_{2 \times 2} } & {{\bf 0}_{2 \times 2} } \cr {{\bf 0}_{2 \times 1} } & {{\bf 0}_{2 \times 2} } & {{\bi R}_2 \lpar m\comma \; n\rpar } & {{\bf 0}_{2 \times 2} } & {{\bf 0}_{2 \times 2} } \cr {{\bf 0}_{2 \times 1} } & {{\bf 0}_{2 \times 2} } & {{\bf 0}_{2 \times 2} } & \ddots & {{\bf 0}_{2 \times 2} } \cr {{\bf 0}_{2 \times 1} } & {{\bf 0}_{2 \times 2} } & {{\bf 0}_{2 \times 2} } & {{\bf 0}_{2 \times 2} } & {{\bi R}_C \lpar m\comma \; n\rpar } \cr } } \right]\comma \; }$$where RML(m, n) is the variance of ML estimate of height h.
C) Data fusion algorithm based on extended Kalman Smoother (EKS)
The multichannel PU algorithm described here is an EKS, which is a two-step technique for nonlinear fixed-interval smoothing consisting of a forward recursion and a backward recursion both followed by a smoothing step which yields the desired smoothed estimate. On the other hand, the fusion of the noisy interferometric data is automatically obtained, thanks to the composite observation vector in (10). The extension of the Kalman filter algorithm to the two-dimensional case requires the introduction of a scanning pattern for the M × N SAR interferograms, for example, we may consider this scanning pattern: left to right, advance one line, then right to left. We can formally model an arbitrary scanning pattern by rearranging the pixels in each of the interferograms with the reordering function f : (m, n) ∈ I → k ∈ {1,2, … , M × N}. Now, adopting a 1D notation, the two-dimensional nonlinear fixed-interval smoother can be summarized as follows:
Summary of nonlinear fixed-interval smoother for multichannel InSAR phase unwrapping.
where 
${\bi \hat x}_f \lpar k\vert k - 1\rpar $ and 
${\bi \hat x}_b \lpar k\vert k+1\rpar $, denoting the forward and backward estimate of the state vector, respectively, are evaluated using (4), whereas their covariance matrices Pf (k|k − 1), Pb (k|k + 1) are evaluated using (5), 
${\bi \hat z}_f \lpar k\vert k - 1\rpar $ and 
${\bi \hat z}_b \lpar k\vert k+1\rpar $ are the forward and backward predicted measurements calculated using (6) and (10). The notation 
${\bi \hat x}_f \lpar k\vert j\rpar $ indicates the estimate of x(k) given the measurements z at pixels 1,2, … j, whereas the notation 
${\bi \hat x}_b \lpar k\vert j\rpar $ indicates the estimate of x(k) given the measurements at pixels MN, MN-1, … , j. In Table 1, 
${\bi \hat x}\lpar k\rpar $ is the smoothed estimate and P(k) is the error covariance matrix of the smoothed estimate. The initial conditions for backward recursion are derived from the last two equations of the smoothing step imposing the following conditions 
${\bi \hat x}\lpar MN\rpar ={\bi \hat x}_f \lpar MN\vert MN\rpar $ and P(MN) = Pf (MN|MN).
D) Filter initialization
Correct filter initialization should guarantee that the covariance associated with the initial estimate reflects realistically its accuracy; in other words, the error in the initial state estimate has to be consistent with the initial state covariance. According to the Bayesian model, the initial state estimate is assumed to be a normally distributed random variable with a known mean and covariance matrix. A practical way to implement consistent filter initialization is by means of the two-point differencing method. Consider the phase measured at the two adjacent pixels around the first pixel (i.e. k = 1) of the scanning pattern defined by the function f, and supposing that the first pixel is 
${\bi \hat x}_f \lpar 1\rpar ={\bi \hat x}_f \lpar 1\comma \; 1\rpar $ then:
$${\bi \hat x}_f \lpar 0\vert 0\rpar =\displaystyle{1 \over {\alpha _l }}\left[{\matrix{ {\varphi _l \lpar 1\comma \; 1\rpar } \cr {\varphi _l \lpar 2\comma \; 1\rpar - \varphi _l \lpar 1\comma \; 1\rpar } \cr {\varphi _l \lpar 1\comma \; 2\rpar - \varphi _l \lpar 1\comma \; 1\rpar } \cr } } \right]$$and the corresponding covariance matrix is:
$${\bi P}_f \lpar 0\vert 0\rpar =\displaystyle{{\sigma _l^2 } \over {\alpha _l^2 }}\left[{\matrix{ 1 & 1 & 1 \cr 1 & 2 & 1 \cr 1 & 1 & 2 \cr } } \right]\comma \;$$where it is assumed that channel l is the channel with the lowest variance. The forward filter starts updating the state estimate at k = 2.
IV. NUMERICAL RESULTS
In this section, we report the experimental results carried out on simulated interferometric data. Three different elevation profiles have been set-up for the purpose of testing different characteristics of the PU algorithm. In the first test case, we use a Gaussian surface in order to check the capability of the proposed method to phase unwrap smooth areas.
The second test case reproduces discontinuities that are typical of the urban environment and is henceforth meant to investigate whether the proposed multichannel PU technique is able to resolve ambiguities due to steep phase jumps (e.g. buildings). The concluding test case evaluates the performance of the proposed method on more realistic interferometric SAR data, generated from the real DEM of the mountainous area around Isolation Peak, CO, USA, available in [Reference Ghiglia and Pritt15], which has been often considered as a test case for evaluating SAR PU algorithms performance (see [Reference Pascazio and Schirinzi5, Reference Ferraiuolo, Pascazio and Schirinzi9, Reference Yu, Li and Bao11, Reference Dai and Zha12]). The availability of true height profile allows us to precisely quantify the accuracy of the proposed technique. All the phase patterns, simulated using the ENVISAT-ASAR parameters of Table 2 in model (1), are corrupted by a decorrelation phase noise with the statistical distribution (9) described in [Reference Bamler and Hartl14] corresponding to a coherence value of γ = 0.5 for all the acquisitions and according to models (2) and (3). Owing to the low coherence value, the phase profiles are quite difficult to be correctly unwrapped using the single-channel approach. The results obtained with the proposed method using the 20 baselines of Table 2 are compared to the single-channel version of our algorithm described in [Reference Chirico and Schirinzi13] and to the MLHE described in [Reference Pascazio and Schirinzi5]; for the sake of comparison the following results will be referred to the reference channel (i.e. channel 1). In the EKS algorithm, a 3 × 3 (i.e. d = 3) pixels window W d and the scanning pattern mentioned earlier (i.e. left to right, advance one line, then right to left) are selected. The amount of process noise Q in expression of the state prediction covariance matrix (5) can be set according to the maximum expected variation of the true absolute phase slope in the interferogram; we observed that Q = diag(0.2, 0.5, 0.5) gives satisfactory results.
Envisat-ASAR parameters and baselines.
In relation to the first test case, the main characteristics of the Gaussian phase surface for the first acquisition channel (this will be referred to as the reference channel) corresponding to orthogonal baseline 1, are shown in Fig. 1(a)–1(c). The noisy interferogram acquired through the reference channel is shown in Fig. 1(d), whereas the comparison between the single-channel EKS PU algorithm [Reference Chirico and Schirinzi13] and the multichannel algorithm proposed here is presented in Figs 1(e) and 1(f). Owing to the presence of noise and the phase wraparound effect, the single-channel approach fails the unwrapping process, on the other hand the success of the proposed method becomes apparent from Fig 1(f). The mean square error of the reconstructed height is 3.78 m.
Gaussian phase surface. (a) True (noiseless) phase referred to ch. 1, (b) true horizontal phase slope referred to ch. 1, (c) true vertical phase slope referred to ch. 1, (d) noisy wrapped phase acquired through ch. 1, (e) estimated phase referred to ch. 1 (single-channel EKS using baseline 1), (f) estimated phase referred to ch. 1 (multichannel EKS using 20 SAR interferograms, baselines 1–20), (g) estimated phase referred to ch. 1 (MLHE using 20 SAR Interferograms, baselines 1–20).
The second test case is related to a discontinuous height profile. The true unwrapped phase profile, the true horizontal, and the vertical slopes of the second test case are depicted in Figs 2(a)–2(c), respectively. A number of pixels that abundantly violate the unambiguous phase difference condition are presented in this image (Figs 2(b) and 2(c)), due to this the single-channel EKS PU algorithm [Reference Chirico and Schirinzi13] is not able to perform PU of the current height profile (Fig. 2(e)), whereas the new multichannel approach can overcome this problem. In Figs 2(e) and 2(f), we compare the results obtained with the single channel of [Reference Shabou, Baselice and Ferraioli10] and the multichannel approach (Fig. 2(f)) proposed here, respectively. It can be seen that the reconstructed phase profile appears very similar to the reference one of Fig. 2(a); the resulting mean square error of the height estimate is 3.57 m.
Urban phase surface. (a) True (noiseless) phase referred to ch. 1, (b) true horizontal phase slope referred to ch. 1, (c) true vertical phase slope referred to ch. 1, (d) noisy wrapped phase acquired through ch. 1, (e) estimated phase referred to ch.1 (single-channel EKS using baseline 1), (f) estimated phase referred to ch. 1 (multichannel EKS using 20 SAR interferograms, baselines 1–20), (g) estimated phase referred to ch. 1 (MLHE using 20 SAR Interferograms, baselines 1–20).
The third test case is related to the real scene around isolation peak [Reference Ghiglia and Pritt15]. The true phase profile of the third dataset, referred as usual to the reference channel, is shown in Fig. 3(a). In addition, the true horizontal and the vertical slopes are presented in Figs 3(b) and 3(c), respectively. Figure 3(d) shows the noisy corrupted wrapped phase acquired by the reference channel. The reconstructed phase profile, again referred to the first acquisition channel and obtained with the single-channel version of the EKS PU algorithm [Reference Chirico and Schirinzi13], is shown in Fig. 3(e) and it presents some artifacts due to the strong noise affecting the phase together with the phase wraparound effect. Figure 3(f) presents the reconstructed phase obtained using the newly proposed multichannel PU algorithm with seven SAR interferograms; it is apparent that the introduction of multibaseline observations permits reliable estimation of the phase profile even in the presence of low coherence data. The mean square error of the reconstructed height is 3.51 m.
Realistic phase surface available in [Reference Yu, Li and Bao11]. (a) True (noiseless) phase referred to ch.1, (b) true horizontal phase slope referred to ch. 1, (c) true vertical phase slope referred to ch. 1, (d) noisy wrapped phase acquired through ch. 1, (e) estimated phase referred to ch. 1 (single-channel EKS using baseline 1), (f) estimated phase referred to ch. 1 (multichannel EKS using 20 SAR Interferograms, baselines 1–20), (g) estimated phase referred to ch. 1 (MLHE using 20 SAR Interferograms, baselines 1–20).
V. CONCLUSIONS
In this paper, the potentiality of the Kalman filter to perform DEM estimation of the observed scene from multibaseline interferometric SAR signals has been investigated. The proposed method optimally fuses the information coming from the MLE with the interferometric data. Also exploiting the spatial information of the neighboring pixels, the proposed method brings together the advantages of the MLHE and those of the Bayesian estimation provided by the Kalman filter. Numerical results have proved the capability of the presented algorithm of performing both PU and phase noise reduction in a single step. Moreover, its capability of reconstructing a larger class of height profiles, overcoming the limitations induced by the presence of large phase discontinuities, has been shown. From the numerical point of view, the presented method is efficient enough and certainly acceptable in terms of computing speed, since it does not require an accurate numerical maximization of non-convex objective functions. In the case of high sloped or urban areas, only a rough ML estimate of height profile is required.
Davide Chirico received the Laurea degree (cum laude) in Telecommunication Engineering from the University of Cassino, Cassino, Italy, in 2009. He is currently working toward the Ph.D. degree with the University of Naples “Parthenope”, Naples, Italy. Since 2010, he has been with MBDA-Missile Systems where he is a systems engineer at the C2&NEC Solutions Department. His current research interests include the fields of estimation theory and radar data processing, with emphasis on nonlinear filtering techniques applied to tracking and navigation problems.
Gilda Schirinzi Full professor of Telecommunications. Gilda Schirinzi graduated in Electronic Engineering in 1983 at University of Naples “Federico II”. In the same year, she joined the Electronic Engineering Department as a research fellow. From 1985 to 1986, she was at European Space Agency, ESTEC, The Netherlands. In 1988, she joined the Istituto di Ricerca per l'Elettromagnetismo e i Componenti Elettronici (IRECE), of the Italian National Council of Researches (CNR), in Naples. In 1992, she was appointed as the Head of Electromagnetics Division of IRECE, and in 1997 she became a Senior Researcher. In November 1998, she joined the University of Cassino as an Associate Professor of Telecommunications. In 2005, she became a Full Professor. Since November 2008, she has been at University of Naples “Parthenope”, in the Telecommunication group. She has taught Signal Theory, Electrical Communications, Microwave Remote Sensing Systems, and Image Processing and Coding. Her main scientific interests are in the field of Synthetic Aperture Radar (SAR) signal processing and coding, SAR interferometry, microwave imaging techniques, and image and signal processing for remote sensing applications.
