Multibaseline SAR Interferometry for Automatic DEM Reconstruction

Abstract

Multibaseline SAR interferometry can be successfully exploited for automatic phase unwrapping and high quality Digital Elevation Model (DEM) reconstruction. The information coming from several interferograms with different baselines increases the elevation ambiguity interval and allows automatic phase unwrapping. The height of each pixel in the image is considered as a random variable: a Maximum Likelihood (ML) estimation of the height is carried out by exploiting the probability density function of the interferometric phase, that depends on the local coherence value. Atmospheric effects are taken into account and considered as a low frequency distortion (with correlation length greater than 1-2 km) of the phase to height conversion function for each interferogram. After the phase unwrapping step is possible to combine all the information available getting a combined DEM that is more reliable than each single DEM. Finally, using multiple images is possible to produce high resolution coherence maps (i.e. an ensemble average is used instead of a space average) that show what actually remains unchanged within all images. Results obtained using ERS-1 data of the area of Bonn (3 days repeat cycle) and Tandem data of two italian volcanos (Vesuvius and Etna) are presented.

Keywords: SAR-Interferometry, Phase Unwrapping, DEM reconstruction, Atmospheric Effects, Coherence Estimation

Introduction

Two of the main problems encountered in the generation of precise DEMs from interferometric SAR images are phase unwrapping ambiguity and elevation distortion caused by atmospheric effects.

When more than one interferogram of the same region is available, a multibaseline approach allows to overcome these difficulties, increasing the final product accuracy (Ferretti A. et al., 1996) (Rocca F. et al., 1997). Baseline optimization procedures together with a good estimation of the standard deviation of the phase noise allow, in fact, the use of powerful statistics techniques, such as Maximum Likelihood (ML) and Maximum A Posteriori (MAP), to properly combine all the information available and to accurately estimate the height difference of each pixel with respect to a reference one with known elevation. The resulting DEM is less affected by atmospheric artifacts, since the uncorrelated atmospheric and noise phase contributions coming from single interferograms are averaged, thus reducing the elevation error dispersion.

When the ML DEM is generated it is possible to get the phase difference with respect to each interferogram. These phase residues are proportional to targets motion and/or atmospheric changes.

The same approach can be exploited to estimate the coherence using an ensemble average instead of a space one. The achieved coherence map highlights what remains unchanged during the time interval between the first and the last acquisition and could be exploited for image segmentation and classification.

In conclusion, the proposed technique can provide:

• a multiple interferograms averaged DEM
• an ''atmospheric'' noise map for each interferogram
• a multiple interferograms averaged coherence map, that gives a measure of SNR on a fine spatial resolution.

The Algorithm

The image is divided in small blocks (about 1 x 1 km) such that:

1. The orbits can be considered rectilinear.
2. Phase distortion due to possible atmospheric effects can be considered linear.
3. The phase to height conversion function can be well approximated, for each interferogram, by a linear function:

dh = A * dr + B * dy + C * dphi +D (1)

where dh is the height variation, dr is the range variation, dy is the azimuth variation and dphi is the interferometric phase variation. All these variations are defined with respect to a reference point chosen inside the block. The reference point will be a pixel having high coherence value in all the interferograms.

In order to compensate for possible baseline errors and atmospheric distortion, the parameters A, B, C and D are iteratively optimized as more and more points are unwrapped. A RLMS optimization algorithm is used to minimize the error between the elevation values obtained from the interferograms and a reference DEM. If no a priori information is available (i.e. a rough DEM or GCPs) the baselines are optimized relative to one of the input interferograms, considered the reference for the other ones: the phase to height conversion function is considered fixed and known for this interferogram. Since the impact of phase noise and atmospheric phase artifacts on the final DEM depends upon the mean value of the normal baseline, the interferogram with the highest value is chosen as the reference one: the higher the baseline, the lower the DEM distortion.

The patching of the unwrapped regions inside each block is operated only when all the blocks are processed. The algorithm starts from the reference point (called block centroid), as in a region growing algorithm (Hock et al., 1995), and it looks among the neighbor pixels for the most "reliable" point to be unwrapped. The unwrapping operation is reliable if all the interferograms estimate the same height variation for the running pixel with respect to the reference point. The dispersion "r" of the probability density function (p.d.f.) of the random variable dh(i,j) around the maximum is an indication of the "reliability" of the pixel P(i,j). The narrower the p.d.f. the more reliable the unwrapping.

In order to compute the p.d.f. of this random variable we use the coherence maps associated to each interferogram. From the absolute value of the coherence it is possible to achieved the expression of the p.d.f. of the interferometric phase ( Lee J.S. et al., 1994) ( Bamler R. and Just D., 1993) and thus of the elevation. The data from the "N" interferograms are "N" independent measures of the same physical variable (elevation), and we can compute its "a posteriori" p.d.f.(f), i.e. the p.d.f. of dh conditioned to the data:

f(dh/dphi1,dphi2,..,dphiN) = K* g(dphi1,dphi2,..,dphiN/dh)*ap(dh) = g(dphi1/dh) * g(dphi2/dh) * .. * g(dphiN/dh) * ap(dh) (2)

where:

K is a normalization constant; g(dphin/dh) is the p.d.f. to observe the value dphin for the interferometric phase of the interferogram "n" when the actual height variation is dh; ap(dh) is the a priori information (this function is a constant if no reference DEM of the region is available).

In Figure 1 it is shown an example of this kind of computation, where three interferograms with three different baseline values are considered. The p.d.f. a posteriori presents a neat maximum corresponding to the ML value.

Figure 1: example of p.d.f. a posteriori computation. Three independent interferograms are considered (the normal baseline values are 106, 146 and 230 m). The coherence value is supposed to be the same and equal to 0.5 (3 looks).

The estimated value of dh (dh_est) maximizes the likelihood function f. When the value of dh_est is determined it is easy to choose the value of 2pi to be added to each interferogram: the correct value of the phase will correspond, in fact, to the height value nearest to the estimated one. In this way is then possible to compute the reliability (r) associated to each transition from the edge of the already unwrapped region to the neighbor pixels.

Many different measures of confidence can be considered for the "a posteriori" p.d.f. ; we have chosen a very simple one:

r = integral[dh_est-M,dh_est+M] f(dh) d(dh) (3)

being M a parameter that fixes the range of allowed altitude variation (e.g. to avoid aliasing, as will be discussed below). The reliability is then the probability that the correct value of the height variation lies inside the interval [dh_est-M,dh_est+M]. It is a positive value less than 1 and can be considered a measure of the multi image "topographic" coherence.

The algorithm will unwrap a point only if

1. This reliability is the maximum value among all the points around the edge of the region already unwrapped.
2. The reliability is greater than a threshold: r(i,j) > THR_REL

When a new point is unwrapped we can carry on the optimization of the geometric parameters and we can compute the reliability values for the neighborhood of this new point. Again the algorithm will look for the more reliable pixel to be unwrapped and so on.

Once the N unwrapped phase maps have been derived, the blocks should be aligned. The centroid Po with the highest value of coherence is chosen as the reference point and the phases of all the points in the image must be unwrapped with respect to it. However, since all the pixels within each block have been already unwrapped with respect to the centroid, the phase alignment is carried out simply unwrapping all the centroids with respect to Po.

The phase alignment is then performed using the same strategy used within the block (ML estimation). Again, the algorithm starts from Po and looks among the neighbor blocks for the most reliable centroid, minimizing the probability to propagate possible unwrapping errors.

The N unwrapped phase maps could now be converted to elevation maps if the non linear relations between phase and elevation were known. That is if the orbital parameters were known with high accuracy. Since in general this is not the case, an optimization procedure has been adopted to maximize the elevation consensus from all the unwrapped interferograms. If no a priori DEM is available, again the largest baseline is assumed as a reference (it is less sensitive to orbital parameters errors and atmospheric artifacts) and its orbital parameters are supposed to be correct. The phase to elevation relation of the remaining N-1 interferograms is first approximated up to the third order, then the parameters are changed to minimize the average (weighted with the coherence) elevation difference with respect to the reference elevation map. Finally all the elevation maps are averaged with a weight dependent on the coherence and the baseline.

The benefits of this algorithm are twofold. In the first place the coherence information is properly exploited: the coherence maps highlight the best path for the unwrapping algorithm and give an estimation of the standard deviation of the error on the final DEM. In the second place there is less risk of aliasing with respect of conventional single interferogram phase unwrapping ( Goldstein R.M et al., 1988 ) (Prati C. et al., 1990). Theoretically it would be enough to have three interferograms with baselines that are prime with respect to each other to remove ambiguities (Chinese remainder theorem). In a practical case, where data are noisy and baselines random, the use of multiple interferograms increases significantly the elevation ambiguity level (Figure 1).

Atmospheric Effects and Multi-Image Coherence Estimation

Once the Digital Elevation Model is available (ML DEM) it is possible to compute N differential interferograms between the original data and the synthetic version obtained from the ML DEM and the optimized baselines values. These phase difference maps can highlight interesting changes in the phase of some targets and/or atmospheric effects due to the change in the refraction index from one acquisition to another. Both effects are clearly visible in the phase error maps we got from the data sets we used to test the algorithm (Rocca F. et al., 1997).

A multi-interferogram approach can be usefully exploited to estimate the coherence using an ensemble average instead of a space one. When a good DEM is available with the same resolution of the SAR images, it is in fact possible to combine the data of all the interferograms available to compute a multi-baseline coherence map of the area of interest on a fine spatial resolution: the increased number of freedom degrees (due to multiple interferograms) allows to get high resolution products (say 20 x 20 m). First the topographic contribution on the phases of each interferogram is compensated for using the DEM. Then the phases are high-pass filtered to eliminate local distortions due to atmospheric effects. Finally the mean phase value is subtracted in each interferogram so that all the data can be considered phase aligned. The final result highlights what actually remains unchanged in all the images.

Experiments

Bonn data set

The topography presents no particular difficulty, but the presence of the river and low coherence forested areas make this area suitable for testing the feasibility of automatic phase unwrapping of not connected pixels. In Figure 2 the multi-image reflectivity map of the area used for the test is shown. This image was simply obtained by averaging the absolute values of all the images available, in order to reduce the speckle noise and to highlight the image features.

Figure 2: Bonn area: multi-image reflectivity map. Dimensions: 600 (range) x 400 (azimuth).

The combined DEM we obtained at the end of the processing (Figure 3) is smooth and presents no relevant phase unwrapping error. Not unwrapped pixels have been interpolated.

Figure 3: Bonn area: combined DEM.

The two banks of the river have been correctly unwrapped even if there is no path between them. Unfortunately, no reference DEM of this areas was available to us, so no comparison was possible.

The Vesuvius data set

The second test area was the region around Mt. Vesuvius. Seven Tandem interferograms were used (Table 1). The baseline values range from 39 to 253 m.

Figure 4: Vesuvius: multi-image reflectivity map. Dimensions: 700 (range) x 540 (azimuth).

No a priori information was exploited during the processing, the final result is presented in Figure5.

Figure 5: Vesuvius: combined DEM.

The estimated error standard deviation map is reported in Figure 6. Black areas correspond to not unwrapped pixels (the reliability was under threshold). The maximum value is 8 m. This map was obtained from the 7 coherence maps and the optimized baselines. The distortion due to possible atmospheric artifacts is not taken into account, so it can be considered a lower bound.

Figure 6: Vesuvius: estimated error standard deviation. For unwrapped areas the range of values is 1..8 meters. Black areas correspond to not unwrapped pixels. Distortion due to atmospheric effects is not taken into account.

Figure 7: Vesuvius: error between the multibaseline DEM and the SPOT DEM of the same area. Black corresponds to -50 m white to +50 m. The standard deviation is 10 m.

Finally, the multi-interferogram coherence map, estimated on a fine spatial resolution (20 x 20 m) highlights interesting features, not visible in a standard (two images) coherence image (Figure 8). Other results concerning this data set are shown in a companion paper (Rocca F. et al., 1997).

Figure 7: Vesuvius: coherence maps comparison. On the left hand side a standard coherence image of the area of interest (April 1996 Tandem interferogram, Bn= 106 m), while on the right hand side the multibaseline coherence map.

The Etna data set

The third test area was the region around Etna volcano. Six Tandem interferograms were used (Table 2). The area selected is about 32 x 27 km. As usual, the images were averaged only by a factor 5 in azimuth direction.

The final DEM obtained at the end of the processing is shown in Figure 9.

Figure 9: Etna: combined DEM.

Figure 10: Etna: error between the Multibaseline SAR DEM and the CNR DEM (courtesy of CNR-Gruppo Nazionale per la Vulcanologia, Reading University and University College London).

Finally, a very interesting comparison is reported in Figure 11. The topographic contribution was subtracted from three different interferograms (August 1995, September 1995, May 1996) to highlight possible atmospheric artifacts, using first the Multibaseline SAR DEM then the reference CNR DEM. A high coherence area (about 5 x 5 km) was selected to reduce the phase noise contribution. The images show that the atmospheric distortions are almost completely removed in the SAR DEM. The phase distortion power ranges from 0.38 to about 1.5 squared radiants. The low-pass trend of these effects is clearly visible.

Figure 11: Atmospheric artifacts detected using the two DEMs (Multibaseline SAR and CNR), considering a high coherence area (about 5 x 5 km). The 3 interferograms used for this comparison are (from left to right) August 1995, September 1995, May 1996. The images show that the atmospheric distortions are almost completely removed in the final DEM.

Figure 12: Power spectrum of the August 1995 data. The added line follows the -8/3 slope to be expected from turbulence phenomena.

Conclusions

This paper describes a multibaseline approach for InSAR DEM generation. It is shown that the combination of more than two SAR images allows to get an automatic technique able to produce high quality results. Coherence maps resolution can be substantially improved due to the increased number of input interferograms. Moreover, the combination of many uncorrelated phase artifacts (mainly due to atmospheric changes) strongly reduces their impact on DEM accuracy. Even if some aspects of the processing chain must be still improved and optimized, the results on real data are good. Comparison with two reference DEMs available of two different test sites showed a good quality of the final products: the error standard deviation is about 10 and 13 m. The next step will be to integrate this software with a geocoding algorithm to obtain a multibaseline DEM on the UTM grid and possibly to combine ascending and descending data to further improve the accuracy.

Acknowledgements

We'd like to thank Ms. Roberta Marchini for the processing of the differential interferograms and the Istituto Internazionale di Vulcanologia - CNR (Catania), the Reading University and the University College London for the reference Etna DEM.

References

Ferretti A., Monti Guarnieri A., Prati C., Rocca F., 1996

Multi Baseline Interferometric Techniques and Applications. Proc. FRINGE '96, http://www.geo.unizh.ch/rsl/fringe96/papers/ferretti-et-al/.

Bamler R. and Just D., 1993

Phase Statistics and Decorrelation in SAR Interferograms. Proc. IGARSS '93, pp. 980-984.

Richard M. Goldstein, Howard A. Zebker, 1988

Satellite Radar Interferometry: Twodimensional Phase Unwrapping. Radio Science, July 1988, 23, pp.713-720.

Prati C.,Giani M. and Leuratti N., 1990

SAR interferometry: a 2D Phase Unwrapping Technique Based on Phase and Absolute Value Informations. Proc. IGARSS 90, 1990, Washigton D.C., USA, pp. 2043-2046.

Hock Lim, Wei Xu, Xiaojing Huang, 1995

Two New Practical Methods for Phase Unwrapping. Proc. IGARSS 95, Florence - Vol. I, pp.196-198.

Lee J.S., Hoppel K.W., Mango S.A., Miller A.R., 1994

Intensity and Phase Statistics of Multilook Polarimetric and Interferometric SAR Imagery. IEEE Trans. Geosci. Remore Sens., 1994, 32, (5).

Rocca F., Prati C., Ferretti A., 1997

An Overview of SAR Interferometry. Proc. 3rd ERS Symposium, Florence 1997, http://florence97.ers-symposium.org

Goldstein R., 1995

Atmospheric limitations to repeat-track radar interferometry. Geophysical Research Letters, 1995, 2,(18), pp. 2517-2520.

Keywords: ESA European Space Agency - Agence spatiale europeenne, observation de la terre, earth observation, satellite remote sensing, teledetection, geophysique, altimetrie, radar, chimique atmospherique, geophysics, altimetry, radar, atmospheric chemistry