Validation of Height Models from ERS Interferometry
The initial processing steps within the InSAR processing chain have become increasingly better documented and well understood in the years since the widespread availability of ERS-1 data first opened up the field of SAR interferometry. Image coregistration (Small D., et. al., 1993a), interferogram generation (Gatelli F., et. al., 1994), and baseline estimation have matured while phase unwrapping is undergoing refinement (Davidson G. and Bamler R., 1996), (Ghiglia D.C. and Romero L.A., 1994). The steps following phase unwrapping are now the subject of increasing research: absolute phase determination, phase to height conversion, forward geocoding, and height model validation (Small D., et. al., 1995b). This paper summarizes some of the latest results at the University of Zürich in the field of height model validation using ERS-1/2 tandem data.
Obstacles to successful height model generation are investigated: coherence maintenance, suboptimal phase unwrapping. In order to test the accuracy of ERS repeat-pass interferometry, test sites with high quality reference DEMs have been chosen. The InSAR processing algorithm and validation methodologies are briefly described, and applied to data from the test sites. Results are presented and evaluated, and conclusions are drawn.
2. InSAR Processing
After SLC image coregistration (Small D., et. al., 1993a), interferogram generation, and flattening (Small D., et. al., 1995b), the image geometry must be refined to ensure successful transfer from slant range geometry into the chosen map reference system.
The accuracy of the precise ERS orbits delivered by the D-PAF is approximately 30cm along and cross-track with a slightly higher radial accuracy (8cm) (Massmann F.-H., 1995). However, successful phase to height conversion requires better accuracy (Li F.K. and Goldstein R.M., 1990).
The geometry must therefore be refined using tiepoints. Two refinements are performed using two types of tiepoints. Height tiepoints are measured in areas of uniform unwrapped phase, while position tiepoints (Northing, Easting) are measured in parts of the image that can be localized with high geometric accuracy (e.g. bridges, road crossings).
Height tiepoints distributed well across range and azimuth are used to refine the baseline model (Small D., et. al., 1993a). An iterative non-linear least squares fit algorithm is used to estimate a set of parameters including the phase constant, the cross-track baseline component and its trend along the azimuth dimension. The radial orbit component is usually left at its nominal value, as it is known more accurately.
Before such refinement (open loop) height estimates have errors on the order of kilometres. Once the phase constant (only) has been localized, height errors are typically reduced by an order of magnitude. Further refinement of the cross-track baseline component (and trend in azimuth) reduces the height estimation error to still lower values.
Position tiepoints are used to refine the SAR image acquisition geometry (as during conventional GTC production (terrain-geocoding)). As with the baseline geometry above, several approaches are also possible here. In conventional ERS GTC production, the azimuth starting point, pixel spacing, near range value, and range pixel spacing are all refined based on tiepoints well distributed across the scene. One may also use the tiepoints to refine the orbit itself. Given a polynomial description of the orbit, one typically refines the constant and linear terms, leaving higher order terms at their nominal values.
Based on the refined baseline and the imaging geometries, the interferometric height of each point in the image with a valid unwrapped phase may then be calculated (Small D. et. al., 1996).
The third equation above is the simplified version appropriate for a 3-parameter datum shift transformation (translation only: ). The more general 7-parameter datum shift transformation is somewhat more complicated, as it adds a scale and three rotation parameters.
Backward geocoding is employed in conventional ERS GTC production (Meier E., et. al., 1993). Height model positions are transformed into WGS84 Cartesian coordinates (Frei U., et. al., 1993), providing a reference frame common with the orbit data. Given precise knowledge of the satellite orbit and each DEM position, one may solve for the satellite position that satisfies the Doppler equation above. The associated range value is then easily calculated as the distance from the satellite to the DEM position, and the (range, azimuth) position within the slant range image is now localized. Resampling into map coordinates follows.
Forward geocoding presumes knowledge of the local height values in a slant range matrix (see the ellipsoid equation above). An initial starting point within the scene is selected and then an iterative method (e.g. Newton-Raphson) is used to solve the set of simultaneous non-linear equations. That solution then serves as the starting point for the solution of the equations defining the next point. A set of irregularly-gridded points in WGS84 geometry are the result. As the points do not coincide with a regular grid in any map geometry, a regridding step follows the solution of the equations.
Both backward and forward geocoding methodologies have been implemented at RSL that directly support full 7-parameter datum shift transformations. Both implementations geocode an ERS quarter scene in less than an hour on a modern UNIX workstation. Provided that the reference and interferometric height models are of comparable accuracies, the quality of the geocoding from the two methods is nearly identical (Small D. et. al., 1996).
The forward method is necessary when one has no reference elevation model, i.e. in the DEM generation application. Other applications are more satisfactorily served by the backward method, as one can generate a DEM of an area once with the forward method, and then use the backward method later successively to terrain geocode all types of slant-range information (e.g. backscatter values, coherence, differential phase). Using the backward method, one is not restricted to height information from a single InSAR pair: the reference elevation model used could be an amalgam formed from several acquisitions, even a mosaic formed through combination of DEMs from different InSAR sensors.
Given a high quality reference elevation model, backward geocoding enables such comparison, but it is not a "blind" test of the end-to-end InSAR elevation model generation, as the quality of the geocoding is non-representative. It does however supply quick feedback on coherence, local incidence angle, height model error interdependencies. It also allows use of a refined height model assembled from heterogeneous sources. Since it uses the reference height values during geocoding, it cannot be part of a "blind" test of InSAR-generated DEMs.
Only a "blind" forward geocoding from slant range into map geometry, with reference heights introduced exclusively at the final height comparison, enables a true end-to-end validation of the InSAR height model generation process. Forward geolocation of areas with poor phase information will be erroneous, and may not register well with ground truth information.
The relative success of the forward geocoding can be tested through overlay of the ERS images in map geometry, geocoded using the forward and backward methodologies. Areas where features do not overlap indicate disagreement in the height models (resulting in planimetric shifts during geocoding).
Validation of the height models is carried out principally through calculation of the difference between the forward geocoded height model and the reference model. RMS difference values together with histograms of the distribution of the height differences provide quality measures.
Note that foreshortened and layover areas undergo extended interpolation when transformed into map geometry - the local slope (and also coherence) therefore influences the locally achievable height accuracy. Combination of at least one ascending and descending pair can be used to mitigate this influence. Consistency checks can be performed between several forward-geocoded height models to increase the height model accuracy.
4.1 Test Site Bonn
The Bonn test site provides an example of a vegetated scene with moderate topography. A variety of baselines and repeat pass intervals are available, as the data was acquired during the first ERS-1 ice phase.
Figure 1 shows the location of the Bonn test site. The Rhine river crosses through the northeast corner of the scene. One finds arable land, mixed deciduous and coniferous forests, open-pit mines, and densely populated areas within the scene.
Two interferograms: 14.03.92 / 17.03.92 and 14.03.92 / 29.03.92 provide a good combination of baseline variety and repeat-pass temporal interval. Precise orbit products were ordered from the D-PAF for use in geocoding and baseline estimation. A reference DEM (originally produced by digitizing map sheet contours) was provided by the GEOS group at the D-PAF for use in geocoding research.
Figure 2 shows a 15-day repeat DEM-flattened interferogram. With the exception of areas where the reference elevation model was no longer up-to-date (some mining pits, gravel pits), topographic fringes are no longer visible. The fringes that remain are likely of atmospheric origin. Note that this interferogram has a time interval of 15 days - the relatively long time interval increases its susceptibility to such differential effects.
The 14-17.03.92 interferogram has a shorter 3-day repeat and a much more sensitive baseline. Regional phase anomalies such as those seen in Figure 2 are not present. The interferogram was filtered, the phase was unwrapped, and the baseline geometry refined. Finally, the digital elevation model was calculated and geocoded. An area within the scene "Bonn-Nörvenich" with terrain variation of approximately 80m, and relatively high coherence values was selected for more detailed study. Geocoded interferometric heights were compared point-by-point with the reference elevations, with an elevation-difference map being produced (Small D., et. al., 1995b). The difference map showed that the interferometric height estimates agree remarkably well with those from the reference DEM. A quantitative representation of the height difference distribution is presented in Figure 3. The histogram shows the frequency of each height difference value within the area. Note that all pixels within the region (of all coherences) were used to compute the histogram, and that no significant systematic height bias is visible. A global RMS error of 2.7m was calculated over the 12x13km area for the interferogram calculated from the March 14/17 1992 pair. Note that the reference DEM was quantized to 1m intervals - the high accuracy was achievable due to the slowly rolling topography together with the large baseline (>400m).
Results were not so satisfying for the March 14/29 1992 pair. The longer time interval both reduced mean coherence over the scene and introduced more coherent phase shifts, distorting the height estimates, and resulting in an RMS error of about 12m over the same region.
The interferometric height model from the Bonn scene was geocoded using both the forward and backward methodologies. Comparison of the geocoded height models showed that planimetric geocoding differences occur where the interferometric height estimates differ (Small D. et. al., 1996). Examples are gravel pits (where the reference height model was no longer current), and forested areas (where the interferometric height estimate was noisy).
4.2 Test Site Bern
Bern was chosen as a test site to investigate the influence of a variety of slopes as well as vegetation states on the DEM generation process. Within the tandem data inventory the Oct. 95 and Nov. 95 pairs were selected as having the best combination of simultaneously high height sensitivity and coherence (Stebler O., et. al., 1996). The more extreme topography and numerous hilly forested areas make phase unwrapping more challenging in comparison to the Bonn scene. Larger height errors are the result. Many areas are rendered inaccessible to phase unwrapping.
DEM flattening was performed on the interferograms to investigate the consistency of their phase behaviours, and the potential accuracy achievable (assuming an optimal phase unwrapping algorithm). No noticeable global phase trends remained after this step. Aside from lakes, rivers, and layover regions, most areas showed remarkably flat phase values. An exception was found in forest stands.
Figure 5 shows an extract from the Nov. 95 DEM-flattened interferogram juxtaposed with a scanned extract from a Swiss cartographic map, courtesy of the Swiss Federal Office of Topography.
Note how the forest stands are clearly distinguishable in the phase image. The behaviour is consistent for all three interferograms studied (Nov. 91, Oct. 95, and Nov. 95). The phase difference can be explained by the fact that the reference height model (originally derived through aerial stereo photography) is based on ground height measurements, while the C-band interferometric height estimates are based upon radar returns echoed from the tree canopy (Small D., et. al., 1995a).
A region-growing phase unwrapping algorithm was used to recover absolute phase difference values. Inspection of height difference maps calculated from the resulting single-interferogram height models showed that phase unwrapping errors had corrupted some calculated height values. Although more robust phase unwrapping methods (Davidson G. and Bamler R., 1996) might have yield improved results, in the presence of variable temporal decorrelation it is inadvisable to generally assume a 100% successful phase unwrapping. We therefore employ a scheme to minimize the impact of such errors.
The accuracy of the final forward geocoded InSAR height model product can be improved through combination of the results from individual interferograms. Areas with inconsistent height values are marked as "no data" areas, joining the areas that weren't within reach of the phase unwrapping algorithm.
The distributions of the single-scene forward-geocoded height differences are shown in Figure 6. Note that phase unwrapping errors produce "islands" of incorrect height values, corrupting the height model. Consistency checking substantially reduces this error source, resulting in the height difference distributions seen in Figure 7. Figure 7(a) shows the height difference histogram for all unwrapped areas. If one restricts the area under study to parts of the image with a local incidence angle in the range from 14-30°, the distribution of height differences improves slightly - see Figure 7(b).
Figure 8 shows the dependence of the height accuracy achieved in the final end-product (forward geocoded) height model on local incidence angle. For local incidence angles ranging from 0 to 45°, the figure shows the mean height error, together with its standard deviation at each LIA. Improvement is significant in comparison to results from a single interferogram. The final height accuracies are respectable: the remaining problem is their lack of omnipresence. Efforts to increase the robustness of the phase unwrapping step are aimed at remedying this problem.
The combined forward-geocoded InSAR height model, with the reference DHM25 visible where no consistent InSAR height estimate was available, is shown in Figure 9. The juxtaposition demonstrates that the InSAR height estimates (where they are available) agree with the reference, but that many regions remain off limits to successful phase unwrapping. Given additional data sets that met the criterion of a large enough baseline (for good height sensitivity) further refinement would be possible (Massonnet D., et. al., 1996). The addition of at least one ascending tandem pair (unfortunately unavailable for this study) would increase the accuracy of the end result by increasing the local resolution in areas foreshortened in the descending Oct. 1995 / Nov. 1995 geometries.
Areal validation of height maps generated by repeat-pass ERS InSAR provides confidence in the InSAR technique. For a 12x13 km area near Bonn, Germany, RMS accuracies of 2.7 m were achieved, with no observable systematic biases over the 40x50 km quarter scene.
Height accuracy increases with longer baselines, and decreases drastically where slopes become extreme, as well as in areas of low coherence (e.g. forest). Spectral-shift filtering dramatically decreases phase variance and is of critical importance for the large baselines that are optimal for the extraction of topography.
Phase unwrapping errors must be either manually corrected, or mitigated through consistency checks with data from other interferograms. Improvement of the height estimate through combination of multiple tandem ERS-1/2 pairs is advisable - combination of ascending/descending pairs is necessary to offer a consistent ground resolution across the scene.
Experience with our Bern scene observing the seasonal variation of 1-day repeat coherence suggests that the winter season appears to be optimum. The minimum coherence useful for mapping is dependent on the slopes within the scene, the required accuracy, and the scene's baseline. However, as a rule of thumb, requiring that the mean coherence be above 0.5 can be a useful discriminator.
The optimum baseline for mapping purposes depends on the slopes that can be expected within the scene. Given a relatively flat scene, with gently rolling topography, a baseline of 300-400m provides the best height accuracy while not excessively sacrificing spatial resolution during spectral shift filtering. Scenes with more significant slopes must make do with smaller baselines and less accuracy, lest they fall prey to unsuccessful phase unwrapping. Consistency checks using multiple interferograms can be used to improve accuracy.
Given the right conditions, repeat-pass ERS interferometry can produce height models with respectable accuracy. The weakness of ERS-1 InSAR height derivation lies in hilly forested areas, where low coherences combine with topography to render height estimation problematic. However, new techniques continue to emerge that work to increase the robustness of InSAR-based height estimation in the face of all known limitations.
This work was supported by subcontract 11318/95/I-HGE with Matra-Cap Systemes, France within a larger contract to ESA-ESRIN. Development of the InSAR processing chain was supported by a contract with ESA under the supervision of the German Remote Sensing Data Centre (DLR-DFD).
ERS SAR data and precise orbit ephemeris information were provided courtesy of ESA-ESRIN. The Swiss reference digital elevation model was provided courtesy of the Swiss Federal Office of Topography. The reference elevation model for the Bonn area was provided courtesy of the D-PAF, for use in SAR geocoding research.
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
|Copyright 2000 - European Space Agency. All rights reserved.|