Altimetry fitting methods for InSAR digital elevation models

Digital elevation models (DEM) are numerical representations of terrain altimetry data largely used in spatial analysis applications. The principles for acquisition, storage, management, update, spatial analysis, visualization and integration with other systems are well known. However, as DEM applications are becoming widespread also increases the concern about the quality of the elevation data and the propagation of DEM errors in derived products. This work is related to Revista Brasileira de Cartografia N 60/03, outubro 2008. (ISSN 1808-0936) 263 DEM generated by synthetic aperture radar interferometry using P and X band collected in Amazon region. The main objective is developing methods for DEM correction that provides flexibility to reduce global and local errors of variable magnitudes commonly present in DEM. The adaptive methodology intends to adjust DEM into terrain landscape by calculating and incorporating error compensation surfaces generated by offset height, triangulation with linear interpolation and inverse squared distance interpolation (IQD) methods. The surfaces calculation and efficiency assessment of methods is provided by analysis of statistical parameters extracted from corrected DEM using a ground control point structure measured by confident survey techniques. Results revealed that application of proposed method, besides improving DEM geometric quality, also provided confident numerical estimation of DEM global and local accuracy.

DEM generated by synthetic aperture radar interferometry using P and X band collected in Amazon region.The main objective is developing methods for DEM correction that provides flexibility to reduce global and local errors of variable magnitudes commonly present in DEM.The adaptive methodology intends to adjust DEM into terrain landscape by calculating and incorporating error compensation surfaces generated by offset height, triangulation with linear interpolation and inverse squared distance interpolation (IQD) methods.The surfaces calculation and efficiency assessment of methods is provided by analysis of statistical parameters extracted from corrected DEM using a ground control point structure measured by confident survey techniques.Results revealed that application of proposed method, besides improving DEM geometric quality, also provided confident numerical estimation of DEM global and local accuracy.

INTRODUCTION
Digital elevation models (DEM) are extensively used in a wide range of earth science fields for different spatial analysis applications.Therefore, nowadays there is a great availability of DEM data representing earth surface altimetry, mainly in developed countries for applications such as orthophotomaps production, flood prevention, control of erosion, agriculture planning, visibility analysis, 3D visualization and so on.Until recently did not exist many ways for DEM data acquisition and generation, but during the last years some effective techniques had appeared that increasingly improved DEM availability expanding potential applications.SAR interferometry (InSAR), for example, is one of new promising techniques that are becoming very popular.The problem with many DEM is that available data do not always present adequate format, resolution and geometric quality for all required applications and while DEM applications are becoming widespread also increase the concern about the quality of the available DEM elevation data and the propagation of DEM errors through the derived analysis.It is well known by researchers that results of many DEM-based quantitative operations are very influenced by magnitude and spatial distribution of elevation errors in DEM.However, currently available DEMs frequently report only the average magnitude of errors as the root mean square error (RMSE) or mean absolute error (MAE), which does not provide information on systematic bias or on the spatial patterns of the DEM errors (HEUVELINK, 1998).
A very important issue to be considered when remote sensing techniques are used in mapping earth resources is the geometric quality and error assessment of data collected and derived that are involved in all processing steps.Countries worldwide generally adopt cartographic standards, some of them enforced by laws, in order to guarantee geometric compatibility with national systematic mapping scales and spatial resolution.In the case of 3D data for DEM applications the most suitable way to evaluate errors and provide adequate correction relies in a set of mapping functions to associate feature positions in DEM to their corresponding in a cartographic reference frame.So far, many efforts have been made to solve registration and geometric correction problems based in this type of approach (AUDETTE et al., 2000;GRUEN & AKCA, 2005).However, relevant research contributions are still required in this study field, mainly for development of effectively adaptive registrations techniques.It is also important to develop procedures that provide geometric control and certify DEM internal quality (POTTMANN et al., 2004).Such quality topics are essential issues for terrain mapping applications.
The present paper describes a developed method for DEM altimetry correction using data generated by InSAR and a control point structure determined by Global Positioning System (GPS) combined with electronic topographic station surveys.The used approach provides statistical estimation of the altimetry correction effectiveness.The DEM data used were generated by SAR P band fully polarimetric system and X band HH polarization mode.
The methodology adopted in this work consists in the main following topics: SAR data acquisition; field topographic control structure determination; geometric quality evaluation of original DEM models; development of strategies and algorithms for DEM error correction and DEM quality evaluation.All computational implementation of strategies and algorithms to reach the objectives were taken in Interactive Data Language (IDL), a computing environment for interactive analysis and visualization of data that integrates an array-oriented language with mathematical analysis and graphical display techniques.
The heart of working methodology purposes optimal DEM fitting to geographic landscape by calculating and incorporating error compensation surfaces generated by offset height, triangulation with linear interpolation and inverse squared distance interpolation (IQD) methods.The compensation surfaces are calculated from comparison of original DEM data positions to corresponding known cartographic land marks.The resulting correction surfaces, calculated as a regular grid with the same cell resolution as original models, are added to DEM in a pixel by pixel basis.The work also included calculation of global and local accuracy parameters that help to characterize the geometric quality of resulting DEM.The final results revealed that altimetry data of study area, like many other currently available DEM, incorporate local and global error patterns and consequently need to be corrected prior to be used in any possible analyses applications.Final conclusions pointed out that application of proposed altimetry correction method, besides improving DEM geometric quality, also provided statistically confident numerical estimation of DEM global and local accuracy parameters that can be used to evaluate error propagation in DEM derived products.

TEST SITE AND DATASETS
The study area that provided data for testing the experiments is a parcel of Amazonian Tapajós National Forest, stated in southern region of Santarém County, inside Tapajós River basin, next to São Jorge village (Brazilian state of Pará).The region is covered by primary forest, re-growth vegetation of various stages, pastures and uncovered ground soil (TIMBÓ ELMIRO et al., 2006).The data acquisition was conducted in different dates.In September 2000 a radar system from AeroSensing RadarSysteme GmbH that operates X (HH polarization) and P band full polarimetric mode was flown over sections of Brazilian Amazon in the framework of cooperation between the Brazilian Army and the National Institute for Space Research (INPE).One-pass X band and two-pass band InSAR data was acquired over the Tapajós National Forest and surroundings.Digital elevation models for both frequencies were generated by using the company proprietary software (MOREIRA, 1996).SAR radiation P band penetrates into forest returning after ground interaction, so producing a DEM related to forest ground (HOFMANN, et al., 1998).Because X band radiation does not penetrate the forest canopy, the DEM produced by X band data is related to the top of forest and is commonly called Digital Surface Model (DSM).If both models generated from X and P band are submitted to adequate geometric correction that provide effective calibration then resulting difference between them is related to the forest height and is commonly known as Digital Height Model (DHM).DHM is an important DEM derived product used in many natural resources applications.During the mission flight corner reflectors were positioned on several points inside test area to ensure absolute altimetry calibration.However, a problem with the airplane inertial system navigation prevented the proper DEM calibration, so the method proposed in this paper intends to improve InSAR DEM calibration (DUTRA et al., 2002).The points of ground control structure were acquired by GPS and topography surveying techniques.GPS technique used dual frequency receivers operating in static relative positioning method.The base station was installed near to the central part of the DEM test area, so the maximum base-rover distance reached no more than seven kilometers, assuring the same satellite constellation to both GPS receivers.Topography survey used electronic station instruments for measuring closed traverses, all starting from two GPS measured points.These survey procedures provided accuracy for control points of about five centimeters.
Some preliminary studies involving the study area can be found in MURA et al. (2001), DUTRA et al. (2002) and TIMBÓ ELMIRO et al. (2006).Figure 1 presents a general view of study area.In the bottom a tile of LANDSAT-TM5 image shows a rectangle that delimits the DEM area which is represented at right side.

Fig. 1 -General view of study area in Tapajós National
Forest -PA -Brazil.

WORKING METHODOLOGY
The general approach adopted to implement the DEM altimetry correction in this work encompasses several steps that can be summarized in the following list of topics: 1) P band and X band InSAR data acquisition; 2) survey of field control points with GPS satellite tracking and electronic topographic stations; 3) conversion and registration of all data to unique terrestrial geodetic frame (WGS84) and preliminary assessment of InSAR DEM geometric quality; 4) InSAR DEM outliers detection and filtering; 5) altimetry correction of InSAR DEM from field control points using adaptive techniques; 6) quality assessment of processed DEM; and 7) information extraction from corrected DEM.Some of the above topics are briefly described in the next sections with special emphasis to topics 5 and 6.

Preliminary assessments of DEM
SAR remote sensing general principles and characteristics of the used survey instruments and methods indicates that P band and X band DEM is supposed to produce similar altimetry values to that obtained by field surveys in cleared areas such as roads, agricultural parcels and pasture fields (HOFMANN et al., 1998;BALTZER, 2001).Preliminary numerical analysis, visual interpretation using available optical sensor images and field observations permitted to accept the hypothesis that all data being studied are very well matched horizontally and so, errors to be corrected are mainly in DEM altimetry component.Thus, corrections calculations, tests and evaluations had been carried out only on the DEM altimetry component.At the same time the preliminary evaluations confirmed good horizontal alignment, they also revealed the presence of outliers and errors in vertical geometry that needs altimetry calibration.Prior to correction surfaces calculation interactive efforts were applied to detect and eliminate outliers found in InSAR data.This step was done by using convolution filters developed in IDL environment that uses statistical criterion to identify suspicious values treating them with local median.

Plane surface method altimetry correction
All DEM error calculations were processed by direct comparison to the field control structure, considered as ground trust positions due to the high precision of ground acquisition methods.The control points were separated in two sets.One set used to calculate the error surfaces and the other set used to provide the quality evaluation tests.The first level of correction applied to DEM consisted of a simple vertical shift, that is, the determination of an altimetry offset calculated by averaging differences between certified field control points and their corresponding in InSAR models.
This procedure results in correction surface corresponding to a horizontal leveled plane (TIMBÓ ELMIRO et al., 2006).The simple mathematical formulation to determine plane correction surface is given in Equation ( 1) where e m is the offset corresponding to the plane surface correction, N is the total number of field control points used, h T is the field control point altitude and h S is the corresponding InSAR altitude.The offset precision can be given by a spread estimator like root mean squared error (RMSE) given in Equation 2. (2) Alternatively precision estimator can be expressed by the mean absolute error (MAE) given in Equation 3.
Figure 2 shows in gray level image representation the correction surfaces for offset method, P band at left and X band at right.The constant grey level for each band represents the value that should be added to original DEM.Dark grey features in P band at left side represent positions of sample ground points used in calculation of correction surfaces.Table 1 shows global quality parameters of DEM offset correction surface for P and X band evaluated in terms of minimum, maximum, mean, absolute deviation (MAE) and standard deviation (RMSE).The results were obtained by comparing the independent set of control points for tests with corresponding points in corrected DEM. Figure 3 shows offset method local quality using precision maps representation for P, at left side, and X, at right side.Precision maps are continuous surface that represents DEM precision point by point associating an estimated RMSE to each DEM individual cell.The precision maps were generated using fifteen sets of ground test points spread over the study area.The differences between test points and the corresponding points in DEM corrected by offset method produced residuals and one RMSE for each set of test points resulting fifteen individuals RMSE.Interpolating the fifteen individuals RMSE related to each set of ground test points results a regular grid of RMSE that represents the precision map of fitted DEM.

Fig. 3 -Precision maps for DEM corrected by offset
surface method, P left side and X right side.
The selected interpolation method used to generate the precision maps was triangulation with linear interpolation that is described in the next section.

Triangulation surface altimetry correction
The offset correction method discussed in the previous step applies a correction with the same value to all points over the entire extent of DEM, so it does not provide an efficient way to treat the local errors of variable magnitudes that are common in DEM.An alternate adaptive approach that fits better in these cases is the calculation and incorporation of error compensation surfaces generated by triangulation with linear interpolation methods.
The triangulation with linear interpolation method adopted in this work uses the optimal Delaunay triangulation.Delaunay triangulation algorithm creates triangles by connecting lines between sample data points.The sample points are connected in such a way that no triangle edges are intersected by other triangles.The result is a patchwork of triangular faces over the entire extent of the grid (SHEWCHUK, 1999).Delaunay triangulation algorithm implemented maximizes the small angles of triangles producing the most equilateral possible triangles.Each triangle defines a plane over the grid nodes that lies within the triangle, with the tilt and elevation of the triangle determined by the three original data points defining the triangle.All grid nodes within a given triangle are defined by the triangular surface (BURROUGH & MCDONNELL, 1998).This method is considered an exact interpolator and it works best when vertices corresponds to significant data points on the area, that is, samples represents peaks and valleys in errors (ERXLEBEN et al., 2002).Figure 4 shows the correction surfaces for DEM generated by triangular with linear interpolation in gray level representation for P band and X band at the left and right side respectively.Dark grey features in P band at left side represent positions of sample ground points used in calculation of correction surfaces.Table 2 summarizes global quality parameters of DEM correction surface for P and X band in terms of minimum, maximum, mean, absolute deviation (MAE) and standard deviation (RMSE).The results were obtained comparing the independent set of control points for tests with corresponding points in corrected DEM. Figure 5 shows triangulation method local quality using precision maps for P, left side, and X, right side.The precision maps were generated in a similar way to that of previous offset method, i.e. comparing the test points with their corresponding in DEM corrected by triangulation method.

Inverse distance surface altimetry correction
The third DEM correction approach applied to solve problems of global and local errors of different magnitudes is the inverse squared distance method (IQD) that is a weighted average interpolator, which can be either exact or smoothing.Using IQD data are weighted during interpolation, so that the influence of one point, relative to another declines with distance from the grid node to be calculated (BURROUGH & MCDONNELL, 1998).This type of method combines the idea of proximity with the gradual change and explicitly supposes that things close together are more related than things apart.
The core of mathematical formulation to IQD interpolation surface is given in Equation 4where j refers to points to be interpolated, i refers to sample data points, d is the distance between both and n is the number of sample data points used to interpolate.

∑ ∑
Normally, inverse squared distance behaves as an exact interpolator because it produces infinities when Σ(d -2 ij) = 0 (i.e. at the sample data points).When calculating a grid node, the weights assigned to the data points are fractions, the sum of all the weights being equal to 1.When a particular observation is coincident with a grid node, the distance between that observation and the grid node is 0 and that observation is given a weight of 1, all other observations are given weights of 0. Thus, the grid node is assigned the value of the coincident observation.One of the drawbacks of inverse distance is the commonly generation of "bull's-eyes" surrounding the solitary observation positions with values that differ greatly from their surroundings within the grid area.A smoothing parameter can be assigned during inverse distance interpolation to reduce this effect by smoothing the interpolated grid (ERXLEBEN et al., 2002).
Figure 6 shows in gray level representation P and X band correction surfaces respectively at the left and right side for DEM generated by inverse squared distance interpolation method.Dark grey features in P band at left side represent positions of ground sample points used in calculation of correction surfaces.Table 3 summarizes the global quality parameters of surface correction generated by inverse squared distance interpolation method for P and X DEM.The meanings of general terms in Table 3 are the same as in Tables 1 and 2. Figure 7 shows IQD method local quality using precision maps for P at left side and X at right side.The precision maps were generated in a similar way to that of previous offset and triangulation method.These results were obtained comparing data from the fifteen independent control point test sets with their corresponding in the DEM corrected by IQD method.

RESULTS AND DISCUSSIONS
The implementation of offset method is straightforward and represents the least general cost of all tested procedures, but it is not the best way to adjust altimetry errors due to its rigid way for compensating variable errors.The remaining two methods are more adaptive than offset to treat global and local errors.Triangulation with linear interpolation, as can be observed in Figure 4, provides a patchwork of triangular faces over the entire extent of DEM.The overall cost of triangulation is greater than offset method but it provides reasonably better results.As can be seen in Figures 4 and 5 errors in triangulation are corrected according to their extent and not with constant values.The general comparison of global parameters in Tables 1, 2 and 3 and local results of precision maps in Figures 3, 5 and 7 shows that numeric parameters favored triangulation and IQD over offset method.IQD provides a much more smooth variation appearance than triangulation but it presents some "bull's eyes" artifacts as can be seen in Figure 6.The general cost of IQD is greater than offset and triangulation methods.Errors in IQD and triangulation methods are corrected according to their extent and not in a fixed way as happen in the offset method.By analysis of general results in Tables 1, 2 and 3 and Figures 2 to 6 it is possible to conclude that IQD produced the best local results, followed by triangulation and offset.IQD and triangulation demonstrated to be equally adequate to DEM correction.The analysis of equations, tables, figures and general results pointed that all three methods confirmed the presence of global and local errors in DEM and that measures of only global errors are not sufficient for characterization of DEM geometric quality.It is also important to characterize the DEM local quality.Figure 8 shows DEM P at left side corrected by IQD method and DEM X at right side corrected by triangulation with linear interpolation method.The results in Figure 8 were obtained by adding the calculated correction surface with original DEM pixel by pixel.

CONCLUSIONS AND FINAL REMARKS
For several possible DEM analyses applications the researchers can produce their own data and, of course, control de quality of the entire process of error propagation.However, production of DEM is not, in general, the major objective of scientific researchers and, so, they prefer to use available DEMs generated by third party.The major problem is that available DEMs do not report clearly its own errors, sometimes hiding information on spatial patterns of the DEM errors.So, the methods presented in this work can offer a potential contribution for DEM evaluation and DEM correction that relies on adaptive ways to treat global and local errors.Main conclusions pointed out that proposed altimetry correction method, besides improving DEM general geometric quality, also provided confident numerical estimation of DEM global and local accuracy parameters.These points are very important issues in DEM analysis due to the possibility of error propagation assessments in DEM derived products.
Future research topics complementing this study, including development of advanced surface correction methods based in splines and kriging are now going on and results will soon be submitted for publication.

Fig. 2 -
Fig. 2 -P band and X band offset correction surfaces in gray level image representation.

Fig. 4 -
Fig. 4 -P band and X band triangulation correction surfaces in gray level image representation.

Fig. 6 -
Fig. 6 -Gray level image representation of P and X band correction surfaces generated by IQD method.

Fig. 8 -
Fig. 8 -DEM P at left side corrected by IQD and DEM X at right side corrected by triangulation.

TABLE 2 -
TRIANGULATION CORRECTION QUALITY PARAMETERS FOR P AND X DEM.