Automated DEM extraction in digital aerial photogrammetry : precisions and validation for mass movement monitoring

Automated procedures for photogrammetric image processing and Digital Elevation Models (DEM) extraction yield high precision terrain models in a short time, reducing manual editing; their accuracy is strictly related to image quality and terrain features. After an analysis of the performance of the Digital Photogrammetric Workstation (DPW) 770 Helava, the paper compares DEMs derived from different surveys and registered in the same reference system. In the case of stable area, the distribution of height residuals, their mean and standard deviation values, indicate that the theoretical accuracy is achievable automatically when terrain is characterized by regular morphology. Steep slopes, corrugated surfaces, vegetation and shadows can degrade results even if manual editing procedures are applied. The comparison of multi-temporal DEMs on unstable areas allows the monitoring of surface deformation and morphological changes. Mailing address: Dr. Arianna Pesci, Istituto Nazionale di Geofisica e Vulcanologia, Sede di Bologna,Via D. Creti 12, 40128 Bologna, Italy; e-mail: pesci@bo.ingv.it


Introduction
Digital Elevation Models (DEMs) play an important role in Earth surface deformation studies: eruptive events, gravitative instabilities, landslides and glacier evolution, geomorphological variations and crustal deformations can be detected and evaluated by means of multi-temporal DEM comparisons (Kaab and Funk, 1999;Baldi et al., 2002;Honda and Nagai, 2002;Mora et al., 2003;van Westen and Lulie Getahun, 2003).
Models registered in a common reference system can be directly compared defining height residuals to identify deforming areas, estimate mass movement and physical surface changes.
Co-registration is a crucial operation which depends on the correct choice and measurement of control points (the final results could be degraded by systematic errors).
Several techniques are suitable for DEM extraction like digital aerial and terrestrial photogrammetry, airborne and terrestrial laser scanning, GPS methodology with its different measurement approaches and active and passive remote sensing, with optical satellite imagery systems (Fraser et al., 2002).
In particular, digital aerial photogrammetry is a powerful tool in surface model generation extracting high resolution DEMs by means of automated image matching procedures (Farrow and Murray, 1992;Heipke, 1995;Kraus, 1998;Mitchell and Chadwick, 1999;Schenk, 1999).
Several institutes and agencies planned and executed photogrammetric flights over Italian regions providing a wide historical data set: in particular, in 1954, 1976 and 1998 an almost total coverage of the peninsula territory was achieved.
At present many acquisitions are performed using digital cameras providing high resolution images directly, while in the frame of digital processing, historical images (film photographs) acquired by means of analogic cameras have to be converted into a digital format using high quality photogrammetric scanners.
DEM internal precision and grid dimension are strictly related to the quality of metric camera, photographs scale, scanner efficiency, resolution at ground, view angles, surface morphology, vegetation, shadows and atmospheric conditions.
In this paper, different automated procedures are analysed and compared to provide accurate DEMs in the shortest time, reducing manual interactions in the frame of the Digital Photogrammetric Workstation DPW 770 Helava (Dowman, 1990;Dowman et al., 1992).
Several applications were performed to define the realistic DEM precision, estimating threshold values for movement detection on different area typologies.

Aerial photogrammetry overview
The photogrammetric technique defines shape, size and position of objects using images taken from different points of view.Images can be acquired by analogic or digital cameras; because digital photogrammetry processes numeric images, the available film frames must first be digitized and translated from a continuous to a discrete data set.At the same time the intensity of the signal is given by a grey or colour scale assigned to each pixel.
In general, using the metric properties of an image, together with the coordinates of several Ground Control Points (GCP), the geometry and the position of an object can be described in a defined reference frame.The processing of stereo digital images can be automated using image matching procedures, based on well defined techniques relative to shape or grey/colour intensity for the same zones (Kraus, 1998).Generally, matching algorithms are classified into three principal groups: the ABM (Area Based Matching), the FBM (Feature Based Matching) and the RM (Relational Matching) (Forstner, 1986;Shapiro and Haralick, 1987;Hannah, 1989;Cho, 1995;Heipke, 1996;Kraus, 1997;Wang, 1998).
In ABM small windows composed of grey values are used as matching primitives to perform comparisons over homologous zones of images until the best correspondence is reached: cross correlation or Least Square Matching (LSM) methods can be used.
In FBM, features (points, lines, edge elements, small regions, or other elements) are extracted in each image individually prior to matching them: the search is performed in terms of coordinates relative to the same objects characterized by a set of attributes (length, size, curvature, average brightness for regions, etc.).
The RM concerns global features that are composed of different local features: the relations can be geometric (angle, distance, etc.) or radiometric (grey level between two adjacent regions, etc.).
Digital photogrammetry is able to acquire many 3D points for high spatial resolution DEM generation through semi-automatic procedures, overcoming the problems of processing time and costs of the analytical approach.
The accuracy of digital photogrammetry mainly depends on the camera-object distance, image quality and resolution; the latter is related to the digital camera characteristics or to the scanner resolution when the 'analogic' one is used.Dealing with terrestrial applications, where the camera-object distances are short (from few to some hundred meters), the precision of extracted DEMs ranges from millimetric to few centimetres.In the case of aerial surveys, the image scale may range from 1:1000 to 1:50 000 or more; the precisions decrease from few centimetres to several meters (Kraus, 1998).Obviously, the image correlation method adopted for the DEM automatic extraction plays a key role in the minimization of errors of the photogrammetric working process.
The standard procedures to generate DEMs are based on fundamental steps which consist in internal orientation, external orientation (registration into a defined reference system) and point extraction.
The internal orientation is performed to define the frame position inside the camera, using the camera calibration certificate (provided by the factory) to correct for distorsion and to assign known coordinate values at specific points, that is fiducial markers (fig.1).
The external orientation is obtained by two subsequent steps: relative and absolute orienta-tion.In the first case the aim is to create the stereoscopic model in an arbitrary relative system using points common to both images (tie points).
The absolute orientation is the transformation of the generated stereoscopic model into an external reference frame (for instance UTM) defined by the coordinates of several ground control points, recognized on the images (fig.1); an s-transformation provides the registration of the two systems.
The aim of this work is to analyse the automated procedures for DEM extraction using the Digital Photogrammetric Workstation DPW 770 Helava, and to validate their internal and absolute accuracy.

The Helava system
The Digital Photogrammetric Workstation DPW 770 Helava processes aerial photographs or other digital images to obtain maps, orthophotomaps, photomosaics, digital terrain models and 3D images.The system is based on image automatic correlation and the software implemented is the Socet Set (SoftCopy Exploitation Tool Set, LH Systems, LLC, 1999).
Concerning the automatic measurement of DEM from original images, it is possible to operate in a short time and to extract precise three-dimensional terrain models: the measuring speed is strictly related to the terrain configuration and to the quality of the digital aerial images: as stated, the three fundamental processing steps (internal orientation, external orientation and point extraction) are executed and DEMs are extracted using algorithms especially designed for different morphological characteristics of the study areas.
The relative orientation procedure can be automated defining a geometry point distribution (tie point pattern), characterized by a specific point number and their position, and requiring the program to check for each point (fig.2).
Naturally, the search for corresponding points is not always satisfied due to the different prospective acquisition of images: the pattern shown in fig. 2 is defined to correctly measure at least 9 points uniformly distributed on the overlapping area, according Von Gruber (Ackermann and Shade, 1993).The method was tested over a single stereo pair, over a strip and over a whole block (formed by more strips).Results indicate that, in the frame of a single stereo pair or dealing with a strip (images sequences), about 50% of points are correctly measured providing a quite uniform distribution redundant data set (the percentage refers to very high point numbers).A manual interaction can be required when a complete block has to be processed because of a low transversal overlapping between strips but, in general, the automated procedure provides a time effective performance.
The Socet Set software allows the automatic extraction of models using the Hierarchical Relaxation Correlation (HRC) algorithm that consists in a Least Square Correlation method (Helava, 1988a,b).The search for the best correlation is performed following 'pyramidal levels': the image resolution is reduced (by a 2 factor, 1:2, 1:4, 1:8, etc.) for each level.The image matching starts from the lower resolution level to the original resolution, gradually improving  the correlation (Miller and De Venecia, 1992).Different approaches are available for different terrain characteristics (table I), (LH Systems, LLC, 1999) and the program is able to use a strategy sequence well suited for each morphological condition (Fabris, 2004).
The whole stereoscopic model is composed of several stereo pairs, and so each natural point belongs to more than two images according to the high overlapping percentage usually adopted (60-70% between subsequent images of a strip, 20-30% between adjacent strips).Correlation algorithms, working on image grey and colour levels, are very sensitive to prospective differences; for this reason the operator may choose the best data for each area.In particular, three different approaches were used and compared based on different overlapping strategies: using 'all the possible overlapping areas', 'each subsequent single pair' and 'the central part on-ly of each single pair'.When the program works in a completely automatic mode, all the possible overlapping areas can be considered, and homologous areas characterized by very different prospectives may produce low correlation values.A manual interaction of the operator is therefore suitable to choose stereo pairs characterized by a correct overlap: generally, the subsequent images in the strip.
A further discrimination concerns the choice of areas internal to the overlap to reduce errors over the image borders: the last operations take a long time.In the next, results relative to DEM extraction following these criteria are represented, concerning a data set composed of 52 images of Stromboli Island distributed along 5 strips at a mean scale of 1:5000 and scanned by the Wehrli Raster Master RM2 Photogrammetric Scanner at 24 µm (1016 dpi) resolution (fig.3).In this case, and in the following examples, all the images used were kept together with their camera calibration certificate.
Stromboli Island is mainly constituted by very steep vegetated zones: buildings and trees are usually distributed over slopes.An adaptive method was used to extract models over a 5 ×5 m grid.As stated, the presence of smoke cause correlation failure so to estimate the realistic accuracies of automatic extracted models, these areas were removed.When a complete images set composed of several strips is analysed, automatic registration can be complicated by a poor or anomalous prospective situation relative to the models used by the stations.
Figure 4 shows the results of the first completely automated process where 'all the possible overlapping areas' are considered and the processing for 'each subsequent single pair', respectively.In the first case, the result is satisfactory only in the lower part of the island (steepest zones are badly represented) while, in the last, performed leading to a set of DEMs which were merged averaging the heights of points in overlapping zones, problems were evidenced only over borders of stereo pairs.The processing for the 'central part of each single pair' was used to reduce these systematic errors: this procedure requires a manual selection of areas to be processed and increases interactive times.
Figure 5 shows the results together with the reference model obtained after an accurate editing.
Each result, that is each extracted DEM, was compared to the reference model: residuals were analysed calculating their mean and the standard deviation values (table II).
The standard deviation values listed in table II are not necessary representative of areas where manual editing is required, but indicate how well the computed model differs from the real one.
The greater the error of the model is, the longer the time required to correct it; in fact the operator has to perform interactive adjustment of the DEM to obtain a visual correspondence between the stereoscopic images and the measured points.In the latter case (C) examined above, the area percentage where editing was necessary, was about 30% of the total surface.

Validation, precision, repeability of photogrammetric DEMs and applications
As stated, the internal precision of extracted DEMs is strictly related to the mean scale of photographs, image quality, pixel dimension and, obviously, morphology of the area: empirically it is estimated about 0.01-0.02% of the mean relative flying height.The real accuracy of data can be verified performing direct comparisons between different models (geo-referenced models) relative to the same 'stable' areas, that is undeformed areas.
Concerning models relative to different epochs, the height differences over the grid points can be calculated leading to residuals distribution characterized by a mean and a standard deviation value.The first indicates the goodness of registration in the same reference system and the second the internal precision.The morphology of the area, as previous described, influences the achievable accuracies and so it will be considered in the following.
DEMs generated at different epochs have to be oriented in the same reference frame when a direct comparison is performed to obtain height differences; this may be accomplished by identifying and measuring the same ground control points whose coordinates are known in an external and stable reference system or using a sufficient number of natural control points recognized a posteriori over all images and measured by means of GPS surveys.Dealing with old images, the absence of a sufficient number of common points, may cause problem in registration introducing systematic effects.An example of the use of aerial photogrammetry for the study of unstable areas, concerns the Vulcano Island (Aeolian Islands, Italy) (fig.6) and in particular the «Forgia Vecchia» slope; in this region many surveys were performed from 1971 up to now.
The main technical characteristics of the last four photogrammetric surveys performed on 1993, 1996 and 2001 are listed in table III, providing information about used calibrated camera and images resolution.
Using ground control points, images of 1996 surveys were oriented and the models were extracted into the same reference system.These points, that are the centres of targets (especially designed to be well identified a posteriori on the images) were rapidly measured during a GPS campaign carried out at the same time as photogrammetric flight (September 1996).
The images relative to different years (1993 and 2001), were oriented using 12 natural points belonging to and measured from the 1996 stereoscopic models, recognized on all the other   images and localized in stable areas, providing the common reference system (table IV).For the oldest surveys (1971 and 1983) only DEMs derived from stereo pairs at about 1:10000 scale, produced by other operators, were at our disposal; in this case the registration was obtained by a least square surface matching procedure with respect to the 1996 model.
The comparisons of DEMs were also performed over four areas characterized by different vegetation and morphological features for the models of 1993, 1996 (at mean scale 1:5000) and 2001 surveys (fig.8a-d).Table V lists the mean and standard deviation values of residuals.
When the presence of vegetation is relevant, manual editing is required to adapt the model to  the real surface: this operation is performed using several visible points on the ground to define breaklines.
Results are in agreement with the expected theoretical values concerning flat areas; obviously the standard deviation increases over zones where vegetation is quite relevant.A further comparison was performed using DEMs extracted from images acquired during the 1996 survey, and characterized by different flying height (mean scales 1:10000 and 1:5000) (table VI).In this case the attention was focalized on image scale dependence.Precision are still in agreement with expected values.
A further check of the accuracy of DEMs was obtained by a kinematic GPS campaign performed over the top of volcano cone, which provided absolute coordinates of points, distributed along independent tracks (fig.9), whose precision was estimated by a cross-over analysis to be about 4 cm.
GPS surveying was executed (September 1996) by three independent operators walking over the edge of the crater, keeping the antenna manually vertical over the ground and collecting data at 1 s sampling rate using trimble 4000SSI receivers.Baselines of a few kilometres were formed with reference GPS stations, fixed on geodetic vertices, and processed by means of the Trimble Geomatics Office (TGO) software.Geographic coordinates were obtained for each measurement epoch using especially designed procedures for kinematic processing to solve ambiguities on-the-fly (Beutler et al., 1995).
Results obtained from kinematic applications and traditional static GPS measurements of control points were given in the WGS84 reference frame, and converted into the UTM cartographic system (32N).The validation of the 1 m×1 m DEM of 1996 survey (1:5000 scale) was assured comparing the GPS data to the photogrammetric model: in fact, height residuals distribution was characterized by a mean value of 0.01 m and a standard deviation of 0.18 m (Baldi et al., 2000).

at mean
The area presents a complex morphology, vegetated zones and a small portion affected by active landslide, not included in the comparisons (fig.10).
Following the same approach described above, results were obtained relative to the whole area minimizing the effect of vegetation by editing procedures (table IX, fig.11).There is a systematic effect in the mean values, prob-ably due to the registration strategy which consisted in an a posteriori recognition on all images of natural points measured with a rapid static GPS survey; considering the photographs scale, the systematic differences are acceptable.
Also in this case the residuals map presents anomalous values in correspondence of dip slopes and vegetated area.In the landslide, which reactivated in 1994 (Autorità di Bacino del Reno, 1998), only a significant but restricted mass movement is detectable in the upper part of the active area and in the lower part of the slope (fig.12).ric surveys (1993, 1976).

Discussion and conclusions
The aim of this work was to analyse the digital photogrammetric procedures for automatic DEM extraction by means of the DPW 770 Helava station (Software Socet Set) and to evaluate the accuracy of the models obtainable in areas characterized by different morphology.
The images used were acquired in recent years over Stromboli Island ( 2001), Vulcano Island (1983÷2001) and Rocca Pitigliana landslide (1976)(1977)(1978)(1979)(1980)(1981)(1982)(1983)(1984)(1985)(1986)(1987)(1988)(1989)(1990)(1991)(1992)(1993); scales range from 1:5000 to 1:37 000.The first step was to check the relative images orientation using a defined tie point pattern: results indicate that, in the frame of a single stereo pair or dealing with a strip (images sequences), about 50% of points are correctly measured providing a redundant data set (the percentage refers to a very high number of points).Manual interaction may be required when a complete block (more strips) has to be processed; this may be due to a low transversal overlapping between strips over steep slopes.The efficiency of default strategies used by the Socet Set software (based on correlation algorithms especially designed for different terrain characteristics) has been demonstrated.Moreover the choice of areas internal to the overlap ('central part of each single pair') reduced errors in automated extraction, reaching higher precision with minimum operator interaction.However the 'editing' is necessary following the automatic procedures to adapt the digital elevation model to the real terrain morphology in about 30% of the whole area.
The real internal precision of models was checked by means of direct comparisons over same areas providing the reliability of photogrammetric technique and the achievable accuracy depending on image quality, mean scale, pixel dimension and terrain characteristics.In particular, it was noticed that these values are quite similar to theoretical values (0.01-0.02% of the flying height) when the morphology is regular and the vegetation is sparse.Obviously the mean image scales play an important role in the final precisions.The photogrammetric models are repeatable (evidenced in table V and table IX) providing the validation and reliability of such technique.
For each application it is necessary to estimate a realistic precision value, selecting areas where there is no evidence of significant mass movements: in this way, the standard deviation of residuals can be representative of model accuracy and can be used as a threshold to detect high deformed zones and estimate morphological changes.
Comparing high precision and co-registered models, surface monitoring is allowed and provided; the described applications, relative to Rocca Pitigliana landslide and Vulcano flank, confirm the efficiency of the methodology.The 3D topographic observations, acquired periodically over unstable areas like volcanoes and landslides, can be used for deformation monitoring and detection of morphological changes.Different studies for modelling processes can be performed based on the quality and resolution of available DEMs.Digital photogrammetry is currently an efficient method for collecting data useful for DEM extraction and recent developments, both in recording techniques and data processing, allow an improvement of resolution and a faster result production.Moreover, the existing historical images used in this work provided reference models, preceding landslide collapse and volcanic activity and allowing mass movement and surface changes estimation.

Fig. 1 .
Fig. 1.Aerial digital image at 1:17 000 scale (Northern Apennine area): four fiducial markers are visible in the corners; the ground control points are uniformly distributed over the overlapping zone.

Fig. 2 .
Fig. 2. Tie point pattern manually defined and used for the processing.The geometry is well suited to cover the stereo pair and to facilitate the measurement of at least one point over each point groups (formed by 4 points).

Fig. 4 .
Fig. 4. DEM automatically extracted using the global elaboration of 'all the possible overlapping areas' and 'each subsequent single pair' respectively (UTM system).

Fig. 5 .
Fig. 5. DEM automatically extracted using the 'central part of each single pair' and the reference DEM obtained by means of editing operations, respectively (UTM system).

Fig
Fig. 8a-d.Area characterized by different morphologies and/or presence of vegetation: a) flat area without vegetation; b) steep area without vegetation; c) steep area and sparse vegetation; d) steep area and dense vegetation.

Fig. 7 .
Fig. 7. Residuals computed over Vulcano Island from 1983 to 1996 (UTM system): it is possible to observe the effects produced by the 1988 activation.

Fig. 9 .
Fig. 9. GPS kinematic: operator tracks are evidenced using different colours for different surveying.

Fig. 11 .
Fig. 11.Residuals relative to the 1976-1993 time span.The area inside the blue line indicate the active slide body.

Table I .
Class of the principal correlation algorithms (based on terrain slope).

Table II .
Comparisons of automated model with reference one: residuals distributions indicate that the case C lead to the best value of standard deviation.

Table III .
Vulcano photogrammetric surveys performed from 1993 to 2001: photogrammetric camera, calibrated focal length, mean image scale, photogram and strip numbers, resolution, pixel and ground pixel dimensions are listed.

Table IV .
Images orientation report: ground control points number, tie points number, residual of external orientation over the 3 cartesian components are listed.Moreover the chosen grid dimension and DTM points number are reported.

Table V .
Results from direct comparisons are listed: mean and standard deviation are estimated over four different areas, comparing about 500 000 data points.

Table VI .
Direct comparisons between models relative to the same area and to the same epoch but obtained from different surveys at different image scales.

Table VII .
Rocca Pitigliana photogrammetric surveys performed from 1976 to 1993: photogrammetric camera, calibrated focal length, mean image scale, photogram and strip numbers, resolution, pixel and ground pixel dimensions are listed.

Table VIII .
Images orientation report: ground control points number, tie points number, residual of external orientation over the three cartesian components are listed.Moreover the chosen grid dimension and DTM points number are reported.

Table IX .
Mean and standard deviation of residuals relative to the 1976-1993 time span in the external landslide area.
Fig. 12.Comparison of DEMs of the Rocca Pitigliana landslide area obtained from two aerial photogrammet-