InSAR Terrain Height Estimation Using Low-Quality Sparse DEMs

Mike Seymour and Ian Cumming University of British Columbia
Dept. of Electrical and Computer Engineering
2356 Main Mall, Vancouver, B.C.,Canada V6T 1Z4


SAR interferometry (InSAR) is a technique for estimating topography using the phase difference between coincident SAR images. Many regions of the Earth have existing topographic data that can be used to assist interferometric SAR processing. In the following, we show how coarse low-quality DEMs can be used in the InSAR processing chain. In particular, coarse low-quality DEMs facilitate phase unwrapping by removing as much of the low resolution topographic phase as possible. The problem of unwrapping the residual phase is easier than unwrapping the flat earth corrected interferogram phase because the local phase bandwidth is reduced in addition to the demodulation of the flat earth correction. Lower local fringe frequency reduces the probability of residues and hence the probability of phase unwrapping difficulties. In addition, the coarse DEM can also be used to calibrate the interferogram phase by providing a height template for a fitting process.

Keywords: SAR interferometry, DEMs, terrain height estimation, baseline estimation.

Published in the Proceedings of the 3rd ERS Scientific Symposium, pages 421-426, Florence, Italy, March 17-21, 1997.

1 Introduction

SAR interferometry exploits the travel time phase information in the complex-valued SAR images to derive terrain heights or terrain height changes [1]. In a single ERS SAR image, the travel time phase information is lost because of the disparity between the size of resolution element and the wavelength of the SAR. However, if two well-correlated SAR images of the same ground scene are available, one can extract a noisy estimate of the travel time phase difference from the phase difference between the images. Through careful processing, trigonometry, and some external parameters, the off-nadir angle of a pixel in an ERS scene can be estimated to accurately place the pixel in the ground range plane (< 10 m error). In effect, SAR interferometry can be used to derive a digital elevation model (DEM) of a SAR scene.

Estimation of topography using SAR interferometry is usually seen as a four step process: firstly, produce the interferogram; secondly, unwrap interferogram phase; thirdly, calibrate the interferometer fourthly, convert unwrapped phase to terrain height. For InSAR processing, specific signal processing steps are required to produce an interferogram phase estimate which has the highest SNR possible (e. g. registration of the SAR images and filtering of the interferogram phase). Phase unwrapping is the difficult non-linear problem of estimating how many tex2html_wrap_inline504 phase ambiguities are required to reconstruct a (possibly discontinuous) interferogram phase surface. The difficulty of this task is evidenced by the plethora of published algorithms for phase unwrapping. Calibration of the interferometer geometry refers to estimating the baseline magnitude (B) and the baseline orientation ( tex2html_wrap_inline508 ) or equivalently, the normal ( tex2html_wrap_inline510 ) and parallel ( tex2html_wrap_inline512 ) baselines.

Clearly, there are two ways that a low-quality low resolution DEMs could possibly simplify InSAR processing: by facilitating more sophisticated flattening to aid in the phase unwrapping process, and by providing a template of terrain heights to aid in the calibration of the interferometer. Both stages seem interdependent: one needs unwrapped phase and the geometry of the interferometer to perform height estimation and one needs height estimates which require unwrapped phase to perform the calibration of the interferometer. However, one can decouple these two operations by using the image registration relation [2]. In the following, we develop an algorithm for estimating topography which makes use of low-quality sparse DEM data.

2 DEMs and InSAR

While the main purpose of InSAR processing is to estimate an accurate DEM, if a coarse or inaccurate DEM is available, it can be used to aid in several steps in the InSAR processing chain. There are a number of low-quality low resolution DEMs publicly available. In particular, the United States Geological Survey (USGS) has made a 30 arc-second DEM of most of the world available. While these data are certainly low resolution (approximately 1 km posting) and low accuracy (90 % confidence interval of tex2html_wrap_inline514 ), we can show that these DEMs can aid in InSAR processing.

2.1 Using Low Quality DEMs for Calibration

The accuracy requirements for the baseline components (normal and perpendicular baseline) are approximately 5 cm to provide RMS height errors on the order of 10 m for ERS data (assuming no other error sources). Although precision orbit modeling could give this accuracy, the orbit data for ERS SAR products does not quite meet these exacting requirements [3, 4] and alternative methods of calibration must be used. One method is terrain height tie-pointing [5]. In this case, registration between the known locations of height values and the SAR image must be made. The method pre-supposes locally accurate measures of topography which are easily registered to the reference SAR image.

Consider estimating the interferometer geometry by minimizing sum of squared error (SSE) between the InSAR and the DEM height estimates:


Assuming no other errors, one can write that the InSAR estimated terrain height is a combination of the true terrain height ( tex2html_wrap_inline516 ) and the perturbation of the terrain height ( tex2html_wrap_inline518 ) due to errors in the baseline magnitude ( tex2html_wrap_inline520 ) and the baseline orientation ( tex2html_wrap_inline522 ) respectively: tex2html_wrap_inline524 . Similarly, the DEM terrain height estimate can be written as the true topography plus an error term ( tex2html_wrap_inline526 ): tex2html_wrap_inline528 . The restated calibration problem is then to minimize


Using differentials, the error due to a calibration error at a particular data point p, can be approximated by


InSAR height errors due to mis-calibration of the interferometer are proportional to the ground range ( tex2html_wrap_inline532 ) of the true data. In the case of no shadow and no layover, the ground range of radar data is a monotonically increasing function which can be approximated by a line of constant slope with an offset. Combining the error due to baseline orientation and baseline magnitude into one variable tex2html_wrap_inline534 yields


For tex2html_wrap_inline536 to be zero, e.g. no error in the estimated geometry of the interferometer, it must be true that


If there are no trends in the DEM errors that are correlated with the ground range variable, one can theoretically get accurate estimates of the InSAR geometry. Note that this relation also means that the mean error in the input coarse DEM is translated directly to the output error in the InSAR DEM. The solution is to process data sets which are as large as possible to minimize the possibility of biases in the coarse low-quality DEMs.

2.2 Using Low-Quality DEMs for Phase Unwrapping

One common step in InSAR processing is referred to as phase flattening where the interferogram phase component that is due to a flat earth is removed. Flattening facilitates interpretation of the topography in the interferogram phase and reduces phase wrapping complexity. The residual phase after flattening is directly proportional to terrain height so the residual interferogram looks much like a contour map of terrain height. The complexity of phase unwrapping varies directly with the number of residues and the number of phase residues is a function of local frequency [6]. Therefore, reducing the local frequency of the interferogram simplifies the phase unwrapping problem. If a low-quality DEM is available, flattening which is more sophisticated then flat earth flattening could be used to remove frequency content from the interferogram.

The real difficulty here is to find some way of deriving an appropriate geometry to jump-start the flattening process. It is a circular problem: one needs the interferometer geometry to derive the topography but one needs some measure of topography to derive the geometry. However, the registration relation between the SAR images gives us a crude estimate of the unwrapped phase. This relationship can be used to derive a rough estimate of the interferometer baseline to allow initial flattening using the coarse DEM data.

In practice, the slope error of a reconstructed DEM depends on the method of reconstruction, the fidelity of the data, and the variability of the terrain under consideration. In addition, for radar applications, terrain height errors translate into errors in slant range position. Therefore it is difficult to predict in advance whether or not a DEM will actually help in the InSAR processing algorithm. However, one can examine the influence of terrain height errors on the phase slope by considering the relationship between slant range and ground range and the error in terrain heights. For the 30'' DEM from USGS, the accuracy statistic cited for the terrain height data is tex2html_wrap_inline514 meters at 90 % confidence interval. Since we reconstruct the gridded DEM in the slant range coordinates using triangulation, the averaging reduces the 90% confidence interval to tex2html_wrap_inline540 meters. For height differences, the 80% confidence interval is therefore about 750 meters. Each of the points in the is separated on the ground by approximately 900 meters so the average slope error at 80% confidence interval is about 0.84 meters of height per horizontal meter on ground. Nominally, for flat earth and a range sampling interval of 4 meters in range, the 80 % confidence interval in slope error is about 8 meters which is not enough to cause unwrapping errors for normal baselines on the order of 300 meters.

The algorithm is summarized in Table 1. It consists of three general parts which are iterated to increase the SNR of the output: preliminary flattening, residual phase unwrapping, and calibration. Currently, we unwrap the residual phase using a weighted least squares algorithm [7]. We take the current estimate of the topography and use it again in the preliminary flattening stage. The iteration is terminated when the number of unwrapped pixels stabilizes. The iterations allow the algorithm to cope with isolated residues due to noise.

3 Simulations

To test the efficacy of the algorithm, a simulation of a noisy interferogram ( tex2html_wrap_inline542 ) was generated with the following parameters tex2html_wrap_inline611 (see Figure 1). The ``true'' height model Figure 2 was generated by gridding scattered data samples of a Terrain Resource Information Map [8]. These are publicly available digital mapping products generated by the Surveys and Resource Mapping Branch of the government of British Columbia. The 90% confidence interval of the DEM is 5 m for unobscured ground areas. The nominal sampling interval is 100 meters for flat areas and 75 meters for areas of steep terrain. 30'' DCW data was gridded for the overlapping area of the TRIM elevation data (see Figure 2) and the algorithm was applied. After 3 iterations, the algorithm successfully unwrapped 96 % of valid data. The success of pre-flattening using the coarse DEM is shown in Figure 1 where the filtered residual interferogram phase was generated by simply lowpass filtering the flattened noisy interferogram. The mean error of the InSAR DEM is the same as the mean error of the DCW DEM: tex2html_wrap_inline613 -35 meters. This is as expected because minimizing the sum of squared height errors will conserve the mean of the input data model. We are using the equivalent of only about 100 DCW data points to reconstruct the input topography model. Since the 90% confidence interval of the data is 650 meters, we would expect that the standard deviation of the mean height estimate to be about 40 meters which agrees well with our result.

The global error performance of the algorithm is summarized by the percentile plot of the residual error after the mean error was removed is shown in Figure 3. The value of InSAR deviation at the tex2html_wrap_inline615 percentile is approximately 10 meters while the original DCW data had a value of 122 m at the tex2html_wrap_inline615 percentile. This represents an improvement of approximately 20 dB. Note also the large tail in the percentile distribution of height errors from the InSAR DEM due to phase unwrapping errors.

4 Sardinia - Preliminary Results

The ERS FRINGE group distributed ERS-1 scenes from frame 801, orbit 241, August 2, 1991 and orbit 327, August 8, 1991 of Sardinia as test data. These data are well known for their high coherence and their mountainous character which causes phase unwrapping difficulties. We extracted the overlapping 30'' DCW DEM data for the scene and processed a subsection of the data. The interferogram and reconstructed DCW DEM were coarsely registered by hand.

The sub-scene processed is 400 azimuth pixels by 1000 range pixels large and is approximately 4 km in slant range by 3.2 km in azimuth. According to the FRINGE baseline listings [9], tex2html_wrap_inline619 meters while tex2html_wrap_inline621 meters. The phase and coherence magnitude of the interferogram are shown in Figure 4. The raw interferogram phase was smoothed by a factor of four in azimuth and half-band filtered [10] but not downsampled in range. The coherence magnitude estimate was made over 11 by 11 pixel blocks using a simple linear phase compensation in range.

The modeled interferogram phase and residual phase after unwrapping are shown in Figure 5. The estimated terrain height model is shown in Figure 6. The range of InSAR estimated terrain is between 0 and 1000 meters which is approximately the range of topography in the region. Detailed comparison of the estimated terrain heights with the true topography is difficult because of lack of a detailed DEM. However, we do have access to some topographical maps. A line of topographical data was extracted from the map sheets and converted to slant range representation. The line of data corresponds to azimuth line 130 in the processed InSAR data.

A comparison of two output height estimates, the input 30'' DCW topography and the map topography is shown in Figure 6. The ``InSAR w/DEM'' topography was generated using the DCW DEM for calibration. The ``FRINGE Data'' topography was generated using the FRINGE baseline listings as calibrating the interferometer using the orbit data and the flat earth assumption yielded results which had very large errors. It is clear that the ``FRINGE Data'' has a very large slope bias in comparison with the ``InSAR w/DEM'' approach. Using the DCW DEM to calibrate the interferometer is a substantial improvement over the calibration using the ``FRINGE Data''.

However, it is clear that the DCW topography is not well registered to the interferogram/map data. The offset is approximately 500 m in slant range or about 1.25 kilometers in ground range. Considering that the original spacing of the DCW sample points is on the order of 900 meters and that no special processing was done for the registration, the results are reasonable. Despite the misregistration, the valid portion of InSAR data from 5500 meters to 7000 meters in slant range agrees well with the map data. However the first part of the scene does not fair as well with an obvious ambiguity error starting at approximately 4000 meters which gets worse as one reaches approximately 5000 meters slant range. Although the range of terrain height is reasonable, there are still ambiguity errors present in the data which could possibly be eliminated by iterating the estimated terrain height with the DCW model to get better registration between the DCW data and the interferogram.

5 Conclusion

We have presented an algorithm for processing InSAR data using coarse low-quality DEMs. The algorithm takes advantage of the prior knowledge in the low-quality DEM to ease the difficulty of phase unwrapping. In addition, the input coarse low-quality DEM is used to calibrate the interferogram. We showed that the primary cause of error in the output terrain height estimates is the correlation between the true ground range of the data and the DEM height errors. In practice, this means one should process as large an amount of DEM data as possible to minimize bias in the InSAR estimated topography. We showed some simulations which yielded an improvement of approximately 20 dB in the tex2html_wrap_inline615 percentile value of error. Finally, some preliminary results for the ERS-FRINGE Sardinia data set were shown where the improvement in output terrain height accuracy using the coarse DEM as input was clearly demonstrated. The topography estimated using the FRINGE listing had a substantial range slope in comparison to the true topography of the region.


L. C. Graham, ``Satellite interferometer radar for topographic mapping,'' Proceedings of the IEEE, vol. 62, no. 6, pp. 763-768, 1974.

M. Seymour and I. Cumming, ``An iterative method for estimating baseline geometry,'' in PIERS'96, (Innsbruck), July 1996.

D. Small and D. Neusch, ``Validation of height models from ERS interferometry,'' in FRINGE'96, (Zurich, Switzerland), ESA, 1997.

C. Reigber et al, ``Impact of precise orbits on SAR interferometry,'' in FRINGE'96, (Zurich, Switzerland), ESA, 1997.

D. Small, C. Werner, and D. Nüesch, ``Baseline modelling for ERS-1 SAR interferometry,'' in IGARSS '93, (Tokyo), pp. 1204-1206, Aug 1993.

U. Spagnolini, ``2-D phase unwrapping and instantaneous frequency estimation,'' IEEE Trans. Geo. Remote Sensing, vol. 33, pp. 579-589, May 1995.

D. Ghighlia and L. Romero, ``Robust two-dimensional weighted and unweighted phase wrapping that uses fast transforms and iterative methods,'' J. Opt. Soc. Am. A, vol. 11, pp. 107-117, January 1994.

Surveys and Resources Mapping Branch, British Columbia Specifications and Guidelines for Geomatics. Ministry of Environment, Lands and Parks, Province of British Columbia, release 2.0 ed., January 1992.

G. A. Solaas, ``ERS-1 interferometric baseline algorithm verification,'' Tech. Rep. ES-TN-DPE-OM-GS02 version 1.0, ESA, 1994.

F. Gatelli, A. Monti Guarnieri, F. Parizzi, P. Pasquali, C. Prati, and F. Rocca, ``The wavenumber shift in SAR interferometry,'' IEEE Trans. Geoscience and Remote Sensing, vol. 32, pp. 855-865, July 1994.

Table 1: Preliminary algorithm for estimation of topography using coarse low-quality DEMs.


Figure 1: Simulation results: example of pre-flattening using coarse DEM data.


Figure 2: TRIM data based simulated topography and DCW DEM used for simulation.


Figure 3: Simulation results.


Figure 4: Interferometric SAR data for Sardinia dataset.


Figure 5: Interferogram phase for Sardinia dataset.


Figure 6: InSAR topography estimate for Sardinia dataset.