Multi-baseline SAR Interferometry using Elaboration of Amplitude and Phase Data

Multi-baseline interferometric Synthetic Aperture Radar (In-SAR) systems can be exploited to estimate the Digital Elevation Model (DEM) of the observed scene without ambiguities and with an increased accuracy, even in the case of high sloped ground regions. The techniques usually used exploit only the interferometric phase information and they are based on Maximum Likelihood (ML) estimation. An important problem to be taken into account is the mutual correlation of the (complex) interferometric images which impedes the closed form evaluation of the interferometric phases likelihood function. Moreover the statistical independence approximation of the phase interferograms is usually adopted. In this paper we present a method exploiting both amplitude and phase of the interferometric images, with the purpose of expressing the multi-baseline likelihood function in closed form, and we show that, when the number of baselines increases, to achieve an higher estimation accuracy the images mutual correlation cannot be neglected. We also show that to obtain a full resolution speckle reduced intensity image from several full resolution multi-baseline interferometric (complex) images, a phase compensation and a whitening operation have to be performed before averaging the data intensities.


Introduction
Synthetic Aperture Radar Interferometric systems (In-SAR) allow the estimate of the Digital Elevation Models (DEM) of the observed scene [1,2], 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). Multi-baseline In-SAR systems exploit more than two interferometric acquisitions, taken along different sensor trajectories with appropriate spatial separation [3,4].
By properly choosing the baselines, different objectives can be reached: increasing the image range resolution [5,6], reconstructing a 3D reflectivity map (SAR tomography) [7,8], and increasing the terrain height estimation accuracy [9][10][11][12][13], overcoming the height ambiguity problem, which is encountered in single baseline systems in the case of steep and/or discontinuous height profiles. As far as the DEM estimation is concerned, multi-baselines interferometric techniques have the potential of providing DEMs with meters accuracy. They usually exploit only the interferometric phases of SAR images, which are directly related to the ground altimetric profile [1]. The height profile estimation is generally performed maximizing the likelihood function (Maximum Likelihood (ML) estimation) of the multiple phase interferograms. The amplitude of the interferometric images, instead, is usually not exploited, since it is not related in a straightforward way to the height profile, and then it does not provide additional information on it, unless an appropriate scattering model, relating the height slopes to the reflectivity intensity [14], is included or other features, evaluated from amplitude and phase interferometric images, are exploited [15,16].
A problem to be considered in the application of multi-baseline techniques is that a closed form of the likelihood function of the multiple phase interferograms is not available when the interferometric phases are mutually correlated, as usually happens in multi-baselines configurations [10]. A closed form has been found only in the case of dual-baseline systems [17] and it has not yet been generalized to the general case of more than two baselines). In the general case of more than two baselines, the interferometric phases are assumed to be statistically independent, so that the multi-baseline likelihood function can be found by computing the product of the single-baseline likelihood functions. Of course, we expect that this approximation will influence the DEM estimation accuracy.
In this paper we present a DEM estimation technique which exploits both amplitude and phase of the multi-baselines interferometric images. The use of the amplitude information has not the objective of improving the height estimation accuracy, but it is aimed to express in a closed-form the "exact" likelihood function of the multi-baseline data, properly taking into account the mutual correlation (coherence) of all the interferometric images. This expression can be easily derived from the well known multivariate Gaussian model of the interferometric (complex) images, obtained in the fully developed speckle assumption [18]. The method presented, besides taking into account the coherence between all couples of images, has the advantage of providing a speckle reduced estimate of the scene reflectivity amplitude. It is shown that, in the assumption of large signal to (thermal) noise ratio, the height estimate is independent from the speckle free reflectivity amplitude, while the reflectivity amplitude estimate requires the knowledge on the scene height profile, since it introduces phase terms which have to be compensated before properly combining the multi-baseline images. Numerical results on simulated data show that, when the number of the baselines increases, to obtain an improvement in the DEM reconstruction accuracy, it is crucial to properly take into account the images coherence. Moreover, they show that to obtain a full resolution speckle reduced amplitude image, both coherence and height profile have to be exploited in the multi-baseline data processing.

Statistical Model of Multi-baseline Images
Let us consider a multi-baseline interferometric SAR system (see Fig. 1 for the case of two baselines B 1 and B 2 ). where L is the number of the interferometric images, N×M is the images dimension, a(n,m) and φ a (n,m) are the amplitude and phase of the radar reflectivity in the (n,m)-th image pixel, X l (n,m) is the complex envelope of the multiplicative speckle noise, v l is additive thermal noise, h(n,m) denotes the height value at the image pixel (n,m), and l α is given by [19]: where λ is the wavelength, B ⊥l is the orthogonal baseline between the l-th antenna and the reference antenna (master antenna), measured along the direction orthogonal to the range direction, θ is the view angle, c is the speed of light, and R 0 is the distance between the master antenna and the center of the scene.
Note that the channel diversity depends on the change of l α with the orthogonal baseline.
The problem consists in estimating the height values h(n,m) starting from the L complex SAR images Y l (n,m), taking into account that the speckle free reflectivity amplitude a(n,m) is also unknown.
Since we consider each image pixel independently from the others, in (1) the dependence on the pixel coordinates can be understood, and in each pixel we can define the data (complex valued) random vector Y=[Y 1 ,Y 2 ,…Y L ] T and the speckle (complex valued) random vector X=[X 1 ,X 2 ,… ,X L ] T , where the T denotes the transpose. Moreover, since after SAR image focusing the signal to thermal noise ratio is tipically sufficiently high [20], in (1) we can neglect the additive thermal noise v l . Then, (1) can be rewritten in a vector form as: In (2) a, j a and Φ are deterministic, while in the assumption of fully developed speckle [18], the complex elements of the speckle vector X have real and imaginary parts which are zero mean mutually uncorrelated Gaussian variables, with variance 2/π (this is a normalization required for obtaining a unit average speckle intensity). Then, X can be modeled as a proper (or circularly symmetric) complex Gaussian random vector [21] with zero mean and covariance matrix C X of size L×L, whose generic element of index (p,q) is given by: where E[·] denotes the expectation operation and γ pq is the coherence between the complex elements p X and q X of the vector X, Note that in the considered assumptions, it is easy to show thar the coherence γ pq assumes real values.
The probability density function (pdf) of the vector X can, then, be written as: From (3) we observe that the data vector Y is related to the speckle vector X by a linear relationship. Then Y is still an L-dimensional proper complex Gaussian random vector with zero mean and covariance matrix C Y , whose generic element of index (p,q) is given by: where α pq =α p -α q .
From (7), we note that Y C depends on a and h, but does not depend on the phase φ a . The pdf of the vector Y is, then, given by: with y=[y 1 , y 2 , …… ,y L ] T Now, exploiting (7), we can write: depends only on h and it is not dependent on the amplitude a.

Maximum Likelihood Estimation of Height Profile and Reflectivity Amplitude
Since the pdf (11) depends only on the parameters a and h, the reflectivity phase value φ a cannot be estimated from the interferometric images and it does not influence the ML estimation of the parameters a and h.
The ML estimate of a and h can be obtained by maximizing the log-likelihood function, obtained by substituting in (11) the measured complex radar images y : The ML solution for h and a is then given by: Note that the log-likelihood function (12) will exhibit an unique maximum respect to a, while it presents many local maxima respect to the parameter h, since h appears in the exponent of the complex exponential terms in (4), each of which is periodic with a period depending on the baseline value. Larger periods correspond to smaller baselines. Then, (12) it will be periodic respect to h with a period equal to the least common multiple of the periods of each exponential term. If the least common multiple of these periods is outside the interval range of the admissible height values for the observed scene, or if at least two baselines are in irrational ratio, the likelihood function will exhibit an unique global maximum in the search interval. Anyway, several local maxima will be present, so that a global maximization procedure will be required to find the solution of (13) respect to h.
To show that we can solve (13) in two separate steps, first searching the solution for h and then for a, we compute the derivative of (12) with respect to h and a, and equate it to zero.
Let us consider first the derivative of (12) with respect to h. Since the first term in (12) is independent from h, we have: Therefore, the ML estimate for h is the value ĥ satisfying the following nonlinear equation: Since, as stated before, the matrix does not depend on a, the values of the maxima of (13) respect to h can be found independently on the estimate of the reflectivity amplitude a, using (15) or, equivalently, through the minimization procedure: Moreover, to avoid possible ambiguities, the solution search should be restricted to a specific interval. This can be done only if some a priori information on the range of variation for h is available. For instance, one can exploit the recent availability of low spatial resolution DEMs of a large part of the Earth surface computed by the interferometric data acquired by the SRTM C-Band SAR system [23] and free distributed by NASA [24].
As far as the estimate of a is concerned, with simple computations, we get: The ML solution for a is the value â which equals to zero the expression (17) and is given by: From (18), we note that the ML solution â for a depends on the height value h, through the matrix Φ , and it is expressed in closed form. In particular in (18) amounts to perform the square root of the average intensity of the data are obtained after whitening the phase corrected images, in this way the phase coherence is restored by compensating the phase terms introduced by the height profile and the data obtained are decorrelated.

Numerical Results
The proposed technique has been tested on synthetic data, which have been simulated using the Envisat-ASAR parameters given in Table 2. Four different sets of three, six, nine and fifteen interferometric images, of size 500×500 pixels, have been simulated following the model given by (3) and using the orthogonal baseline values given in Table 2.
The multiplicative (complex) speckle vector X in (3), for each image pixel, has been generated as a realization of a proper complex Gaussian random vector with zero mean and covariance matrix C X , whose generic elements Xpq c , given by (5)  corresponds to the complete decorrelation of the interferometric images (γ=0), and in this case the two interferometric images cannot be used to form the phase interferogram. Such baseline is usually referred as the critical baseline. Note that the baselines of Table 2 are all smaller than the critical baseline, which in the case considered is 1081 = ⊥c B m. In Figs. 2(a) and 2(b) the amplitude of the radar reflectivity map and the height profile (measured in meters) used to simulate SAR images are shown. Fig. 2(c), instead, shows the inaccurate a priori DEM exploited in the estimation, exhibiting an height error with maximum values of m 20 ± respect to the true DEM, and with a standard deviation of 3.5 m. Note that the height error values producing a phase of π/4 in correspondence of the minimum and maximum baselines are 3.5 m and 1.2 m, respectively. Then, height error error values between 1.2 m and 3.5 m produce a not negligible unknown phase factor at least in one of the images considered, while height error values larger than 3.5 m produce an unknown phase factor in all images of the considered set. Since the height estimate is independent on the value of the reflectivity amplitude a, as shown in the previous Section, we can first estimate h using (16). Note that in (16) the knowledge of the speckle covariance matrix C X is required. As it is usually done in interferometic applications, it has been estimated from the data by using a sliding window whose dimensions have to be small enough to assume that the reflectivity amplitude and the height values in the window pixels are constant, and large enough to have a sufficiently accurate estimate of the statistical average [25]. The size generally adopted is 5×5 pixels. Moreover, we note that when the number of baselines increases and the baselines assume values which are close each other, the matrix C X becomes ill conditioned, and special care has to be devoted to its inversion, which can be performed by uding singular value decomposition (svd) [22].
The height estimation accuracy can be evaluated by computing the root mean square error (RMSE h ). We remind that the RMSE h is the result of two contributions: statistical errors due to the presence of noise, and errors due to solution ambiguity [10]. The latter contribution can be minimized by restricting the solution search interval around the height values of the low resolution inaccurate a priori DEM, which for instance cam be provided by . In such case, the RMSE h can be assumed dependent only from the phase noise, and can be consistently compared with the square root of the Cramer Rao Lower Bound for h (CRLB h ), which represents the maximum accuracy attainable with an efficient estimator [26] and it does not take into account the possible occurrence of ambiguous solutions. Then, to reduce the occurrence of ambiguous solutions,we search for the ML solution in an interval 40 m wide around the a priori height profile of the scene. The RMSE h values in case of three, six, nine and fifteen images are presented in Table 3, while the height reconstruction obtained from six images and fifteen images are shown in Figs. 3(a) and 3(b) respectively. Table 3 also shows the square root of CRLB h, to assess the maximum accuracy that can be achieved by processing the considered data sets. As expected, RMSE h decreases and approaches the corresponding CRLB h 1/2 values when the number of images increases. The better quality of the height profile obtained with 15 baselines is evident also in Fig. 3.  To show the better performance of the presented method, it is useful to make a comparison with the method of estimating the altimetric profile starting from only phase data. Since in this case a closed form of the joint pdf of the phase interferograms is not available, it is commonly assumed that the interferometric phases are statistically independent, so that the joint pdf can be easily computed by performing the product of the marginal pdfs of each phase interferogram [9]. This assumption is not usually verified in the case of multi-baselines interferograms [10], since the spectra of the different images are overlapping each other. The adoption of this approximation is as more critical as more the number of baselines increases, since in this case the image spectra overlapping becomes more relevant. Fig. 4 shows the results obtained from only phase data in the case of six and fifteen images and it shows a worse quality respect to the reconstruction of Fig. 3. The RMSE h values are presented in Table 3. It is noted that the estimation exploiting only phase data exhibits a worst performance. The difference in performance becomes more evident when the number of images increases, since the adopted approximation of statistical independence among the multiple phase interferograms becomes more critical when the baselines step up within the critical baseline extent.  As far as the estimation of the amplitude reflectivity a is concerned, it can be evaluated using (18). Note that the estimate (18) depends on h through the matrix Φ=Φ(h) (see (4)) and on the speckle covariance matrix C X . To analyze the influence of h and of the mutual correlation of the multi-baseline images on the estimation of a , we can consider the following three cases: The matrix is approximated with h 0 is the a priori inaccurate DEM. This amounts to use a rough approximation of the phase factors introduced by the height profile. The image correlation, instead, is taken into account through the matrix X C . In this case (18) becomes: Eq. (22) amounts to perform the average of the intensities of the images obtained after whitening and phase compensating the L measured images for the available a priori DEM h 0 . Note that phase compensation partially restores the phase coherence, while the whitening operation uncorrelates the data.

3) Fully coherent interferometric multi-look (FCML): in this case the matrix Φ=Φ(h) is approximated with
, where ĥ is the estimated height value. In this case (18) becomes: Eq. (23) amounts to perform the average of the intensities of the images obtained after whitening and phase compensating the L measured images using the estimated height profile. Note that in this case phase compensation produces an higher phase coherence, since the used height values are more accurate, while whitening uncorrelated the data.
Figs. 5, and 6 show the IML, the PCML and the FCML reconstructed image amplitudes, in the case of six and fifteen images, respectively. As expected, comparing by a visual inspection Figs. 5 and 6 with Fig. 2(a), it is quite evident that FCML outperforms the other methods. where a is the true amplitude reflectivity value and â is its estimated value, and the equivalent numbers of looks (ENL), defined as:. For a single image, the ENL is about equal to one, since the variance and the square statistical average assumes the same value. Note that the higher the ENL, the better is the speckle reduction and, then, the capability of distinguishing different intensity levels. The values of the parameters defined above, evaluated on the amplitude images obtained starting from three, six, nine and fifteen images, are shown in Table 4.
The ENL calculated on an uniform reflectivity patch of the initial image amounts to 1.05, as expected. The ENL obtained with more images in the FCML case increases with the number of images. This behaviour is not met in IML and PCML cases. Note that in all cases, the ENL obtained starting from L images is smaller than L, due to the mutual images correlation. The results in Table 4 show that when the number of images increases, performance of FCML becomes noticeably better than IML because the effect of the mutual correlation between the images becomes more significant. In fact, we note that, due to the data correlation, the intensity average performed by IML produces only a small effect for the speckle reduction on the reflectivity amplitude, which does not increases when the number of the images increases. Moreover, comparing the accuracy parameters obtained with PCML and FCML processing, we observe that the effect of the improvement of the height profile accuracy is always evident and it increases when the number of the images increases. Table 4 also shows the corresponding square root of the CRLB for a, to assess the maximum accuracy achievable from the considered data sets.

Conclusions
In this paper we presented a method to improve the DEM estimation accuracy, starting from correlated multi-baseline interferometric images and a low resolution inaccurate a priori DEM of the observed scene. The method takes into account the correlation among all the interferometric images. This is very important, since multi-baseline acquisitions are almost always mutually correlated. The method exploits both amplitude and phase of the interferometric data and it is based on an ML estimation technique. It allows to estimate both the height profile and the speckle reduced reflectivity intensity. It has been shown that the estimation procedure can be split in two steps: first the height estimation is performed independently from the reflectivity amplitude, then the speckle reduced reflectivity amplitude is evaluated taking into account both the images correlation and the phase factors introduced by the height profile.
Numerical results on simulated data show that the obtained height estimation accuracy is better than that achievable with the methods exploiting only phase data, which adopt the approximation of statistical independence of the phase interferograms, due to the difficulty of expressing their "exact" likelihood function in closed form. Performance improvement is as larger as the number of baselines considered becomes higher. Moreover, as expected, it is shown that the height accuracy improves with the number L of the processed images with a low slower than the linear one, due to the images correlation.
As far as the estimation of the reflectivity amplitude is concerned, numerical experiments show that both the image correlation and the phase factors introduced by the height profile of the observed scene have to be considered. In particular, it is shown that the exploitation of the improved accuracy height estimation obtained in the first step, allows to a noticeably improvement of the reflectivity amplitude quality.