|Mike Seymour and Ian Cumming||University of British Columbia
Dept. of Electrical and Computer Engineering
2356 Main Mall, Vancouver, B.C.,Canada V6T 1Z4
Keywords: SAR interferometry, DEMs, terrain height estimation,
Published in the Proceedings of the 3rd ERS Scientific Symposium, pages 421-426, Florence, Italy, March 17-21, 1997.
SAR interferometry exploits the travel time phase information in the complex-valued SAR images to derive terrain heights or terrain height changes . 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 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 ( ) or equivalently, the normal ( ) and parallel ( ) 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 . In the following, we develop an algorithm for estimating topography which makes use of low-quality sparse DEM data.
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 ), we can show that these DEMs can aid in InSAR processing.
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 . 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 ( ) and the perturbation of the terrain height ( ) due to errors in the baseline magnitude ( ) and the baseline orientation ( ) respectively: . Similarly, the DEM terrain height estimate can be written as the true topography plus an error term ( ): . 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 ( ) 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 yields
For 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.
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 . 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 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 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 . 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.
To test the efficacy of the algorithm, a simulation of a noisy interferogram ( ) was generated with the following parameters (see Figure 1). The ``true'' height model Figure 2 was generated by gridding scattered data samples of a Terrain Resource Information Map . 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: -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 percentile is approximately 10 meters while the original DCW data had a value of 122 m at the 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.
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 , meters while 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  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.
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 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.
Table 1: Preliminary algorithm for estimation of topography using coarse low-quality DEMs.