Abstract

A digital elevation model (DEM) created by automatic image matching or laser scanning – also named as LIDAR, includes a not negligible number of points, not located on the terrain surface but on buildings or vegetation or even miscorrelation. The manual refinement of such a DEM is time consuming. In general, points located not on a continuous surface, which can be differentiated, have to be identified and removed. A Digital Filtering Technique to be applied to the automatically acquired DEM-data is presented. The strategy is based on Linear Prediction of stationary random function after trend removal. The filtering was applied to Photogrammetrically acquired DEM-data via automatic digital correlation techniques. Different terrain compositions like buildings and canopy density; terrain roughness; grid sizes; negative scales have been investigated

The system can also be used in connection with DEM acquired through LIDAR. For both applications very acceptable results have been achieved.

Introduction

Digital Elevation Models (DEM) acquired using photogrammetric correlation techniques or LIDAR have the disadvantage that the resulting modelation surface may not represent the bare terrain but the visible surface including vegetation and buildings, the so called Digital Surface Model (DSM).

The first step to be done after data acquisition is the removal of the points not belonging to the terrain. An automated method using linear prediction has been applied. The method and results are explained.

Digital Elevation Model Filtering

The elimination of points not belonging to the terrain surface is known as Filtering. There are several methods or procedures for interpolation and filtering. Among them are:

a) Splines approximation

b) Shift Invariant Filters

c) Linear Prediction

d) Morphological Filters

Although Morphological Filters are the most frequently been applied, Linear Prediction is a very robust methodology for the filtering of Digital Surface Models.

Let us consider the following three random functions l(u), s(u), and r(u), such that

l(u) = s(u) + r(u)

The observable function is l(u) and r(u) represents the noise; s(u) is called the signal. Interpolation and filtering are therefore the problems of finding an estimate ŝ (u_o) of the random function s(u) at u = u_o, when a discrete set of a function values l(u₁), l(u₂),…. l(u_n) from a given realization l(u) are given.

The favored estimate is such that its results from a linear combination of l(u), or:

ŝ_o = ŝ(u_o) = a^t l

with: a^t = [a₁, a₂,…a_n] is a vector of coefficients and

l = [l(u₁),l(u₂),..,l(u_n)] is the vector of given data values,

that means, the estimated value is a linear combination of the given data values. This is particularly important for those cases where a function cannot be or it can be extremely difficult to be represented in an analytical form, i.e., Digital Terrain Models.

It can be proven that filtered signal value ŝ_o = ŝ(u_o) is:

ŝ_o = s_sol S_ll^-1 l (1)

with:

ŝ_o= predicted or filtered value of the signal at u=u_o

s_sol = is a vector containing the cross-covariance between signal s at u=u_o (s_o)

and observations l_i. (equals to covariance between signals s_iand s_o)

S_ll = Covariance matrix

l = Vector of centered measurements

The covariance matrix is constructed from the covariance function that has general form:

COV = A. e ^-(2)

This expression states that the covariance between two points P_i and P_k is dependent on their reciprocal distance. If the points are close to each other, then the covariance is high. The covariance tends to zero with growing distance between points. A is the vertex value of the function and consequently is the covariance for zero distance. B represents the slope of the covariance function. These parameters are known or they are to be determined empirically in each case.

A fundamental prerequisite for the application of a covariance function is the removal of the trend from the given observations or signals. Only in this case the covariance between points will only be dependent on their reciprocal distances and in such a case we will be handling stationary random functions. The elimination of the trend is accomplished by using a very low degree polynomial or a moving plane. The result of this trend separation is the vector l_i that contains the centered points of measurements. These values l_i describes the deviations of the sample-measured points from the moving pane or the low degree polynomial.

Finally, we can conclude that the elements of the covariance matrix S_ll are in fact the covariance values between the point measurements. The main diagonal elements are the variances of the signals (centered point measurements after trend removal) and the off-diagonal elements are the covariance values between same above signals. That is:

C = (3)

The Program DTMCOR

Developed at the Institute for Photogrammetry and Engineering Survey of the University of Hannover, the software analyses and filters a Digital Elevation Models. The programs locate firstly possible Spikes or Blunders by the introduction of a tolerable minimum and maximum height.

The elimination of the Trend is carried out via the use of a moving plane.

Z_i = a₀+ a₁X_i+ a₂Y_i (4)

The area covered by the DEM is divided into a mesh of equal size. The dimensions of the square grid are suggested by the program based on the distribution of the points. For the processing of the points of a particular mesh (Processing Unit = 1 mesh), the points located in the 8 surrounding meshes are considered. In this way, the moving plane coefficients are computed using the points located in the 9 contiguous meshes (Area of Consideration), via least squares. See Fig. 1

The trend removal is based on the tilted plane in the Area of Consideration resulting in the centered measurement values l_i.

Assuming a normal distribution of those l_i values their standard deviation (s_z) is computed by the program and a multiplication factor (fac) is entered into the program for the computation of a threshold or tolerance factor (T_z).

T_z = fac s_z (5)

All those points within the corresponding processing unit whose deviations (i.e., centered measurement values li) are exceeding the above tolerance (T_z) are excluded. The trend separation is repeated in a loop until no more defective heights are recognized by the system. The erroneous heights of the processing mesh are deleted from the records. The erroneous heights of the neighboring patches (i.e., area of consideration) remains for the computations of the next patch. They only remain unconsidered for the current patch.

Figure 2 shows the above-explained iterative process of trend removal and elimination of erroneous height values.

The figure 2 shows a typical terrain height-profile. It consists of 10 points, three of which are blunders. The inclination of the moving plane represented by the colored straight lines, varies considerably with each iteration. Four iterations are needed in the example, until the moving plane stabilizes and fits the terrain surface and no more blunders are identified. The standard deviation decreases drastically from iteration to iteration with the removal of the largest blunders to the smallest.

In relation with the Linear Prediction DTMCOR makes use of a Covariance Function with the following form:

C(P_iP_k) = A . e ^-1.30103 (6)

A and B are parameters of the function. A is the vertex-value of the signal-covariance function. It is a filter factor for the normed covariance function. Its specifies the relationship between random and systematic components of the height discrepancies (i.e., centered measured heights l_i). A value A = 1.0 (in the program limited to 0.99), means no random errors are available. The parameter B is the slope of the function and represents the distance at which the influence of points is reduced to 5%. On the other hand its value also limits the width of the mesh for the local prediction.

The interpolated surface is defined as in equation (1) where the main diagonal elements of the covariance matrix contains variances V_ll = 1, meaning all measurements are regarded as being of the same accuracy. As the vertex-value A have been found to be appropriate at 0.7, interpolation and filtering are possible. The higher the vertex-value of the function is, the smaller the variance s_l² is, and smaller is the filtering effect (See Figure 3).

Figure 3.

Covariance Function

The differences between the centered measured values at each DEM point and the predicted value according to formulae (1), are computed. Once again assuming normal distribution of the above discrepancies, their corresponding standard deviation (s_zp) is computed. A multiplication factor is introduced in the program for the calculation of a threshold or tolerance value (T_zp)

T_zp = fac s_zp (7)

Text Box: If the computed discrepancies are exceeding the threshold, the corresponding DEM points are also eliminated in a loop fashion.

Figure 4. Surface of prediction Corresponding real height values

Text Box: Figure 5. Profile through area including a point not belonging to the surface

(left: left hand profile, center: center profile, right: right hand profile

red line = surface of prediction in these profiles based on the area,

points = real height point)

Experimental Tests

DEM's automatically acquired using digital correlation techniques have been filtered using DTMCOR. Three different terrain areas have been used, namely:

Area Type A: Urban flat area with heavy building density, buildings of

different heights and dimensions, moderate canopy.

Area Type B: Urban Flat area with moderate building density, regular building

size and height, moderate to heavy canopy

Area Type C: Flat urban / Industrial area. Large open spaces, mixture of low

altitude, small and big size buildings, low to moderate density of

canopies.

Three different image scales, 3 DEM spacing and a DEM in a TIN fashion have been analyzed. The results are shown in Table 1:

AREA TYPE	IMAGE SCALE	DEM Spacing (ft)	TIN	<min Z >max Z (%)	Tilted Plane (%)	Linear Predict. (%)	Total Excluded (%)
B	1:6000	30		16.7	29.6	6.1	52.4
		50		16.46	28.7	5.74	50.91
		100		14.67	26.4	4.88	45.95
B	1:9600	30		16.6	28.9	6.3	51.8
		50		16.5	27.8	5.73	50.03
		100		14.8	25.9	7.99	48.70
B	1:9600		*	16.6	25.14	7.32	49.42
A	1:7920	30		7.32	38.38	15.0	60.7
		50		7.89	39.01	14.9	61.8
		100		9.0	34.9	15.4	59.3
C	1:7920	30		3.14	23.0	11.3	37.44
		50		3.14	22.71	10.84	36.69
		100		2.89	23.24	17.26	43.38

Table 1 Percentage of excluded points

A close look of Table 1 reveals:

a. For given terrain coverage there is a very small increase of the eliminated points with smaller DEM spacing. This is more noticeable in Area B with regular building size and heights.

b. Regardless the type of terrain coverage, the largest percentage of eliminated points is carried out by the iterative trend removal procedure. This corresponds to non-surface terrain points such as correlated points on building roofs. A key parameter for the good performance of the iterative trend removal procedure is the chose of the weight factor for lower points. In all cases a numerical value of 3s_h was used.

c. In general the performance of the method does not change considerably with the use of Triangular Irregular Networks (TIN's).

d. The Image Scale does not influence the performance of the procedure.

e. Linear prediction works better in cases of very abrupt changes of heights. This can be observed in terrain cover type A.

f. The performance of the minimum - maximum height filter depends on the actual knowledge the minimum and maximum terrain heights, which in the majority of cases is not known, but it is supported by a frequency distribution.

Text Box: Figure 6. Area and contours before filtering

Text Box: Figure 7. Area and contours after filtering

Figures 6 and 7 are showing a piece of the terrain cover type C with overlaid contour lines. Figure 6 shows the contour lines before filtering the Digital Elevation Model and Figure 7 after filtering. Contour lines of Figure 6 honors correlated points on rooftop of buildings and tree canopies. Contour lines of Figure 7 are representing the terrain.

Conclusions

The above digital filtering strategy based on iterative trend elimination plus linear prediction has proven to be a very effective tool for removing non-terrain points of a Digital Elevation Model based on digital correlation techniques in a soft copy photogrammetry environment and also by LIDAR.

The filtered Digital Elevation Models can be used without excitation for production of digital orthophotos. Contour lines shall be derived with the aid of breaklines and digital graphical editing.

Further investigation is required using other terrain coverage and other DEM geometry including breaklines.

Acknowledgments

Special thanks are given to Mr. David Betzner for the DEM data acquisition and testing of the filtered data, and to Mr. Julius Burkus for the preparation of some of the above figures.

References

1 JACOBSEN, K. (1980):.Vorschläge zur Konzeption und zur Bearbeitung von Bündelblockausgleichungen, Doctor thesis, Institute for Photogrammetry and Engineering Survey. University of Hannover.

2 KRAUS, K. (1997): Eine neue Method zur Interpolation und Filterung von Daten mit schiefer Fehlerverteilung. Österreichische Zeitschrift für Vermessung & Geoinformation 1/1997, Seite 25-30

3 MIKHAEL, E. (1976): Observations and Least Squares.

4 PASSINI, R. (1997): Transformation of the Geodetic Network of the Emirate of Dubai from UTM Zone 40 Projection System to DLTM (Dubai Local Transverse Mercator). GIM 2/1997, 25-28