| Mike Seymour and Ian Cumming | University of British Columbia Dept. of Electrical and Computer Engineering 2356 Main Mall, Vancouver, B.C.,Canada V6T 1Z4 http://www.ee.ubc.ca/sar/ |
|
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.
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 
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 [2]. 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 [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 ( 
) 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 [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 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 [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.
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 [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: 
-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 [9],
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 [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.
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.
