Geometrically Accurate Time Series of Archived Aerial Images and Airborne Lidar Data in a Forest Environment

Reconstructing three-dimensional structural changes in the forest over time is possible using archived aerial photographs and photogrammetric techniques, which have recently been introduced to a larger audience with the advent of digital photogrammetry. This paper explores the feasibility of constructing an accurate time-series of archived aerial photographs spanning 42 years using different types of geometric data and estimation methods for image orientation. A recent airborne laser scanning (lidar) data set was combined with the image block and assessed for geometric match. The results suggest that it is possible to establish the multitemporal geometry of an image block to an accuracy that is better than 0.5 m in 3D and constant over time. Even geodetic ground control points can be omitted from the estimation if the most recent images have accurate direct sensor orientation, which is becoming a standard technique in aerial photography. This greatly reduces the costs and facilitates the work. An accurate multitemporal image block combined with recent lidar scanning for the estimation of topography allows accurate monitoring and retrospective analysis of forest vegetation and management operations.


Introduction
Forests exist in a permanent state of change in response to natural and human-driven processes.Monitoring of forests is crucial to applications in forestry and nature protection.Historical records such as results of forest surveys provide information on the historical landscape structure from the stand up to the regional level (Axelsson 2001).Historical records of the state, management operations and natural events in stands would be informative, but they are rarely available apart from those pertaining to research forests.
Remote sensing (RS) methods were studied for the appraisal of changes in forests (Fensham and Fairfax 2002, Rhemtulla et al. 2002, Saksa et al. 2003, Heikkonen and Varjo 2004, Yu et al. 2004).However, RS has been most extensively employed for the assessment of state variables.
Aerial photographs have certain advantages over other data in image-based, long-term multitemporal RS, e.g. they allow three-dimensional (3D) measurements for accurate estimation of trees and canopy (Miller et al. 2000, Fujita et al. 2003, Korpela and Tokola 2005).However, in most retrospective studies of forest vegetation, the analyses have been carried out in 2D using monoscopic, geometrically rectified images.This approach is of limited precision due to the inherent spatial mismatch between images (Taylor et al. 2000, Fensham andFairfax 2002).For example, accurate orthorectification requires that both the sensor orientation and relief used for differential image correction are accurate and detailed.An unambiguous canopy relief cannot be defined.Moreover, geometric changes are a part of the vegetation change detected in multitemporal analysis.Another aspect is the time span, e.g. in Finland, a system for topographic photography was developed in the 1930s.Aerial photography has since been developed, which means that the archived material represents a large variation in the factors that define the usability for vegetation analysis (Hall 2003).
Photogrammetric techniques have until recently been the exclusive domain of the specialists.The concept of digital photogrammetry (Sarjakoski 1981) has only recently emerged.Costly photogrammetric instruments are being replaced by digital photogrammetric workstations in which many of the arduous routines have been automated (Schenk 1999).Direct sensor orientation (DSO) systems have revolutionized photogrammetry (Krauss 1997, Cramer et al. 2000).In DSO, the exterior orientation of the camera is observed during the flight.The advent of digital photogrammetry and new sensor techniques such as airborne laser scanning (ALS) are also widening the opportunities available for new applications in forestry (Leckie et al. 2003, Naesset et al. 2004, Korpela 2004).
A prerequisite for precise, 3D retrospective photogrammetric measurements is accurate solution of the exterior orientation of all images in a 3D coordinate system that remains constant over time.Aerial triangulation (AT) is the technique commonly used for solving exterior orientation.Conventional AT employs ground control points (GCPs) and tie points, which are observed in the images.Old photographs often lack GCPs; thus, the geometry needs to be brought back in time by some means.The hypothesis here is that it is achievable using multitemporal image observations of immobile points and postpositioned GCPs, which can even be completed or replaced by recent DSO observations.The use of DSO is motivated by the possibility of establishing the geometry of a multitemporal image block with a small number of GCPs, or, in the extreme case, with no GCPs, which would reduce the costs.
The objective here was to explore the feasibility of constructing an accurate time series of archived aerial photographs and recent ALS data for the retrospective, precise 3D analysis of forests in the conditions encountered in boreal forests.The geometry of a multiscale, multitemporal image block spanning 42 years was determined using 3 versions of AT, all of which apply multitemporal tie points but differ in the use of GCPs and DSO observations.The different versions were computed to examine the stability of AT, geometric errors in the GCPs and image observations, and the feasibility of solving multitemporal AT without GCPs.
The results of a recent ALS (lidar) run were combined with the photogrammetric time series and examined for the geometric match using natural targets.The idea in combining lidar and photogrammetry is to use lidar for the determination of terrain elevation; a constant element over time, and to use images for vegetation observations (cf.St-Onge et al. 2004).A good geometric match between different types of observation over time is needed.

Area of Multitemporal Interest
The study area (12 km × 17 km) is located in the municipalities of Orivesi, Juupajoki and Ruovesi in Southern Finland (61°50´N, 24°20´E).Scots pine (Pinus sylvestris L.), Norway spruce (Picea abies L. Karst.), and silver and downy birch (Betula pendula Roth, Betula pubescens Ehrh.) are the main dominant tree species.Deciduous trees are in full leaf from late May until late September.Tilly soils with stones and boulders on the surface are frequently encountered.Spruce occurs naturally at these sites and attains the heights of 25-27 m at the age of 100 years.There are also sites with assorted soils on eskers and glacial delta formations with few stones or bedrock on the surface; these barren sites are dominated by pine.Lush grass-herb sites are rare; however, most sites on mineral soil develop dense covers of grass and shrub vegetation 1-5 years after a clear-cut.The elevation varies from 135 m to 200 m above sea level and the topography is hilly in places.The basins are covered with lakes, ponds, open mires, spruce swamps and pine bogs.The peatlands were largely drained in the 1960s and 1970s.Many of the ponds are in the process of terrestrialization, a process presumably speeded up by drainage.Forests have been actively managed and there has been a dense network of forest roads since the 1970s.The few settlements in the area are scattered along cultivated fields.
The geodetic infrastructure is normal for the conditions occurring in Finland.There are 10-12 trigonometric geodetic points per 100 km 2 ; these are located on hilltops and were measured using field triangulation between the early 1950s and late 1970s.Geodetic satellite positioning was initiated in 1988 in the region (Fig. 1); these points are located on roadsides.Since 2002 VRS™-GPS (Virtual reference station) has been possible in the area for accurate and rapid positioning, e.g. for DSO.All geodetic points had planimetric XY coordinates in the Kartastokoordinaatistojärjestelmä KKJ, the national Gauss-Krüger map projection system used in Finland.Elevations were given in the N60 system, the National system of orthometric heights.The combined KKJ/N60 system was considered to be a tangential XYZ coordinate frame, which is required in photogrammetric projects.The effects of Earth's curvature and geoid undulations were neglected (Kraus 1997).The elevations of the geodetic points had been determined alternatively, using trigonometric observations or geodetic satellite positioning.In Finland, the focus has been on planimetric coordinates at the expense of elevation.Only recently has there been a shift in thinking triggered by the adoption of geodetic satellite positioning.The nominal accuracy of the elevation observations varied significantly; the trigonometric elevations could have been off by as much as 1-4 m (Korhonen 1966).The best accuracy was assumed for those GPS points measured and adjusted relative to points with levelled, accurate elevations (e.g.98M3743 in Fig. 1).The point set was densified in 2000 using RTK-GPS (Real-time kinematic).Stones visible in aerial images taken in 1997 were positioned relative to trigonometric point 76M7627, with a relative accuracy of app.0.05 m in 3D (Korpela 2000).The elevation was levelled to 9 points, using the RTK-GPS points as references.

Image Material 1962-2004
Aerial photographs were requested of two archives and 288 photographs between 1962 and 2004 were included (Fig. 2).Large and medium-scale images were favoured.The first photographs for forestry purposes were taken in 1962 and the first topographic colour and colour-infrared (CIR) photographs in 1973.The proportion of forestry photographs increased in the 1980s; these were mostly in CIR film and taken under leaf-on conditions at a scale of 1:30000, which has been customary in Finland.The scale of the topographic photographs (1946−) was between 1:16000 and 1:60000.The area is exceptional, with large-scale photography for research, teaching and road planning.Some of the oldest images had unusable areas degraded by clouds and the early CIR images had suffered from overexposure.The original films were scanned in 1997-2004, using resolutions of 14, 15 or 21 µm.
Archives provided the camera calibration certificates and available flight documentation.The calibrated camera constant, sometimes referred to as focal length or principal distance, location of the principal point in the fiducial coordinate system, coordinates or separation of the fiducial marks, radial and tangential lens distortions, optical or photographic resolving power and information on chromatic aberration were found in the certificates (cf.Krauss 1993, p. 29-79).The exterior orientation parameters i.e. the XYZ coordinates of the projection centre and the 3 camera rotations had been observed for 2 flights in 2004 using a direct sensor orientation system.

Lidar Campaign in 2004
ALS (e.g.Wehr and Lohr 1999) was carried out on August 5, 2004 (Table 1) over an area that matches flight 97251 (Fig. 2).A reference GPS was at point 88M5041 during the campaign.The XYZ geometry of the lidar points was thus established relative to 88M5041 with a trigonometric elevation.The planimetric KKJ coordinates were transformed to the ETRS89 (~WGS84) system of the GPS.The needed ellipsoid height was computed using the FIN2000 geoid model (Ollikainen 2002).N60 elevations are 18.69 m lower than the ellipsoid in the area.The ALS consisted of three 6-km-long strips with 20% overlap and 2 shorter strips flown perpendicularly at the ends of the long strips.Blom Geomatics AS in Norway flew the campaign and processed the data.The lidar points were delivered in the conformal WGS84− UTM35 projection with ellipsoid heights; these were transformed back into the KKJ/N60 coordinate system (JHS 2002).

Methods for the Image
Orientation of a Multitemporal Image Set

Problem Statement
The objective of image orientation is to establish, for each image, an exact mapping from the object space to the xy film plane: (X, Y, Z) ⇒ (x, y).The pinhole camera model (1) was chosen for this process.The so-called collinear equations for a camera ray are: (1) In (1) x and y are the film plane (camera) coordinates of the object point (X, Y, Z), c is the camera constant and x 0 and y 0 are the coordinates of the principal point in the film plane.(X 0 , Y 0 , Z 0 ) are the object coordinates of the projection centre and r ij are the elements (directional cosines between the coordinate axes of the XYZ and xyz systems) of an orthogonal 3 × 3 matrix, which are non linear functions of 3 camera rotations ω, φ and κ (Krauss 1993).

Fiducial Mark Transformation and Film Deformations
In digital photogrammetry, image observations apply to features in the x´y´ or column,row coordinate system of digital images.The computations (1) require metric xy camera coordinates.Affine transformation was used in fiducial mark transformation for mapping between the 2 systems.To study film deformations, the mappings were also estimated with Helmert transformation, which is rigid and does not compensate for affine deformations.
The cameras had 4 or 8 fiducial marks, which were measured in the images using a semiautomatic algorithm in the EspaBlock program (EspaSystems, Espoo, Finland).Some fiducial marks were positioned manually by applying strong contrast enhancement.
In general, film deformations induce errors in the camera model (1).It is possible to extend (1) by additional parameters to compensate for film deformations, lens errors and errors in the calibrated values of c, x 0 , and y 0 , but this was not needed here.Camera calibration has to be carried out in combination with aerial triangulation (Krauss 1997) if a calibration certificate cannot be found for a particular camera or non-metric, uncalibrated amateur cameras have been used.

Hybrid Bundle Block Adjustment
Aerial triangulation (AT) was carried out using hybrid bundle block adjustment by applying weighted least-square (WLS) estimation (Kraus 1997).Hybrid refers to the combined use of image observations of ground control points (GCPs) and tie points and of direct sensor orientation (DSO) observations of the unknown exterior orientation parameters.Bundle block adjustment is equivalent to nonlinear regression and calls for an iterative solution, in which the unknown parameters x k are first given approximate values that are gradually improved.In WLS, corrections are solved into vector x using In (2), A is the matrix of observation equations, P a diagonal weight matrix for the homogenization of observations and y the vector of residuals.The nonzero elements in A correspond to the differential quotients of the linearized error equations of observations calculated using the approximate values of the unknowns.The elements in P are the inverses of the a priori variances of observation errors.
A program was written for constructing different versions of A, P, and y for (2).The partial derivatives of the linearized collinear equations (1) with respect to the unknowns, needed for A, were taken from Krauss (1993, p. 279-280).Each image point yields two rows in A, one for each coordinate.If DSO observations were used in solving a version of (2), the error equations of the DSO observations were added to A (partial derivates of DSO error equations are constant and equal to 1), the inverses of a priori variances to P, and the residual vector y was extended to include the residuals in the DSO observations, e.g ∆X 0 = X 0(DSO) -X 0(adjusted) .
The accuracy was evaluated by computing the mean-square error of unit weight (σ 0 ) and the mean-square errors of the unknown parameters (σ xk ).
In ( 3) n is the number of observation equations and u is the total number of unknowns (Krauss 1993, p. 386).The mean-square errors σ xk of the individual unknown parameters x k were computed using the main diagonal elements q kk of the inverted (dispersion) matrix Q: Matlab was used for numerical solution of ( 2)-( 5).Table 2 lists the a priori standard deviations (SDs) used in the weight matrices P in (2).

Work Flow of Multitemporal Aerial Triangulation
AT was carried out in phases.Subsets of images were attached to the image block, after which (2) was computed.The number of blunders in observations and initial approximations was kept small and could be found by examining y.
(2) converged after 2-3 iterations and elements in x became negligible.Serious errors in observations or in initial approximations resulted in the divergence of x.The work advanced backwards in time starting from the 2002 images (Fig. 2).The 50 images from 2004 with both DSO and image observations completed the block.GCPs and tie points were observed monoscopically using the EspaBlock program and an experimental photogrammetric program.Latter was used for completing the observations when a satisfactory solution was available and (1) could be used for verifying if a point lacked an observation in an image.The initial approximations of the unknown coordinates of tie points and elevation points were determined from maps.Camera angles ω and φ were set to zero, assuming perfect vertical photography, and the initial approximations for κ, X 0 and Y 0 were determined from the photo-indexes and flying height was used for Z 0 .
There were full GCPs with known XYZ coordinates and elevation points with a known Z coor-dinate that had been established in 1964-2000 (Fig. 1).For example, the RTK-GPS-measured GCPs were postpositioned only in 2000 and had signals in 2002 and 2004.Image observations of these GCPs apply to natural objects before 2002, of which the earliest was from 1966.Similarly, elevation points were observed only in the images of 2002 and 2004.The earliest image observation of a signalled geodetic GCP was from 1999.To make use of the signalled GCPs (1999)(2000)(2001)(2002)(2003)(2004) and the DSO observations of 2004, it was necessary to find multitemporal tie points (Fig. 3).In the full block of 288 images there were 279   tie points, 9 elevation points and 22 full GCPs, which gives 2583 unknown parameters.There were 3599 observed image points (Fig. 4).47% of the tie points were multitemporal with a time span ranging from 1 to 42 years.

Fiducial Mark Transformations
Helmert transformations had larger RMSEs for photographs stored for a lengthy period before scanning (Table 3) indicating the presence of film deformations.For example, 8 frames of flight 97248 were scanned in 1999 and had a mean RMSE of 10.9 µm (affine 6.9 µm).Remaining 8 were scanned 4 years later and had a mean RMSE of 63 µm (affine 6.7 µm; Fig. 5).The differences in the RMSEs between the affine and Helmert transformations were negligible for photographs that had been scanned shortly after development.
A difference was observed between images from the 2 archives; the other archive uses machinery for maintaining constant temperature and humidity, while the photographs in the other archive are stored in a normal office room.

Aerial Triangulation with Ground Control Points and Direct Sensor Orientation Observations
The DSO observations of flights 04402 and 04403 were used in solving (2) with the weights in Table 2.The maximum number of observations was used in this AT, referred to as the full AT.

General Accuracy of the Image Block
The value for σ 0 (3) was 0.9965; σ 0 ≈ 1 for wellperformed measurements, a good model, and   properly chosen weights.The largest standard errors (SEs) of the exterior orientation parameters (Table 4) were obtained for a frame of flight 00106, which covers an area outside the convex hull formed by the GCPs and had only 6 tie points located in the other half of the frame.
Similarly, the largest SEs in the tie point coordinates were observed for isolated points on the borders of the block (Table 5, Fig. 6) and the SEs in the X, Y and Z coordinates were strongly correlated.Thus, the spatial pattern for X in Fig. 6 represents the Y and Z coordinates.
Photogrammetric XYZ coordinates were computed for the 22 GCPs, using least square ray intersection (Korpela 2000, p. 101-102) with from 2 to 67 image observations per GCP.The solutions were compared against true values (Fig. 7, Fig. 8).The largest differences in X, Y and Z (−0.22 m, +0.34 m and +0.60 m) were observed for the detached geodetic point 98M3743 seen in two 1:30000 photographs.Point 88M5041, the GPS reference in the lidar campaign, was observed in 42 images and its position was 0.19 m below the true value.The mean differences between the DSO observations of the projection centre coordinates (X 0 , Y 0 , Z 0 ) and solutions of AT for flights 04402 and 04403 were 0.02 m, −0.03 m and 0.04 m (SDs 0.023-0.029m).These are within the a priori SDs of 0.1 m, indicating reasonable congruence between the GCP and DSO observations.

Accuracy of Image Pairs with Stereo Overlap
The accuracy of the image block was analysed by forming 313 image pairs with over 50% overlap (normalization/rectification for stereo viewing was not carried out), in which the points were computed photogrammetric coordinates in 2533 image pair × point combinations using least square ray intersection.These were compared with the solutions for AT (tie points) or fixed, true values (GCPs).
For 1872 tie point observations, the SDs of the differences in X, Y and Z were 0.16 m, 0.16 m and 0.37 m, respectively; these are larger than the mean SEs in AT (Table 5).For GCPs, the SDs of the differences were 0.13 m, 0.14 m and 0.36 m. 70% of all tie point and GCP observations refer to a signalled point.The average obtained accuracy of about 0.15 m for X and Y, and 0.35 m for Z, is to be interpreted as the achievable accuracy of an observation made in an arbitrary image pair assuming further that the target is as clearly detectable in the images as the GCPs and tie points were.For example, the top of a deciduous tree is not as clear a target as the crest of a roof; hence, the accuracy for treetops will be less.The types of points used in AT are given in Table 7.
Images that represent different scales and image pairs with varying overlap are present in the block (Fig. 2).These factors affect the accuracy of photogrammetric measurements and were examined by selecting image pair × point combinations of the signalled GCPs from 2002 and 2004 (Table 6).The same image observation was used in all combinations of image pairs in computing the photogrammetric XYZ position, which are therefore not independent.The effect of scale is seen in the increasing inaccuracy shown when the scale decreases from 1:6000 to 1:30000.The accuracy obtained for 1:30000 is likely too optimistic, due to the small number of observations.Large image overlaps and focal lengths decrease the estimation accuracy of Z, which can be seen by comparing a scale of 1:14000 with scales of 1:12000 and 1:16000.

Accuracy and Value of Different Types of Ground Points
In multitemporal AT, it is preferable that the image points are unambiguously measurable and remain immobile over time.Image residuals of the RTK-GPS-measured GCPs (stones) had a mean L2 norm (length) of 16.4 µm in 221 image observations from 1966 to 2000, in which the GCPs were observed without a signal.There were 398 observations with a signal in the images from 2002 and 2004, which had a 9.2-µm mean L2 norm for the residuals.The difference may also be explained by the presence of film deformations and differences in scale.The GCPs had signals only in the recent images, which were generally taken from a lower flying altitude.It is also possible that the points  The GCPs and tie points were determined per object type by visual interpretation of aerial images and the image residuals were analysed by object type (Table 7).The smallest mean L2 norm was observed for a junction of a ditch.It cannot be concluded that such points are particularly advantageous, because of the small number of image observations.This tie point was situated on the edge of the image block and out of the coverage of the flights with DSO.The XYZ coordinates of this tie point were adjusted in AT to the image observations because there were no nearby GCPs or tie points observed in the images with DSO that would have forced the geometry of the bundle block.A low mean L2 norm for the residuals of 3.8 µm indicates good accuracy, but it is misleading.
The signalled points had a mean L2 norm of 8.7 µm, which is less than the a priori measurement accuracy (Table 2).The largest residuals were observed for the GCP 88M5041, which is a small, almost symmetric patch of bedrock in a cultivated field.This GCP was observed without a signal in 25 images from 1973 to 2000.It is likely that the image observations did not apply to the same object point over time, because the borders and extent, i.e. the symmetry, of the bare rock with respect to the exact GCP position may have changed.
In general, finding multitemporal tie points was difficult in large-scale images with forest only, if there were no clearly observed areas with stones; in such cases even treetops of dead trees had to be selected.Three types of object did not qualify as multitemporal tie points.Road markings could not be considered because they are repainted frequently.Shadows move with the apparent movement of the sun, and small tussocks of low vegetation change their appearance rapidly.Insulators on electric poles, although smaller than the image pixels, were often clearly detectable and were considered immobile.A hummock or hollow in mire vegetation was used if the vegetation appeared to be unvarying.This does not guarantee that the elevation remained unchanged, since the waterlevel may have also varied.Similarly, buildings were used only if there were no signs of alterations in the structures over time.Residential buildings were considered immobile and thus favoured.Stones along forest roads were avoided.Corners of roof structures were used only if the crest of the roof could not be identified.In general, it was considered that any motion or deformation of multitemporal points would induce geometric, random or systematic inaccuracy in the image block.

Aerial Triangulation with Ground Control Points Only
AT was solved without DSO observations such that the image observations of the 22 GCPs and 9 elevation points were the only support for exterior orientation.The DSO observations of flights 04402 and 04403 were omitted in solving (2).At the solution, the 279 tie points and 9 elevation points shifted with respect to the positions of the full AT (Fig. 9) by −0.01 m, −0.01 m and 0.15 m in X, Y and Z, respectively; the SDs of the shifts were 0.03 m, 0.02 m and 0.22 m.The shifts were spatially correlated and large in the southern part of the block (Fig. 9, Fig. 10), which is outside the convex hull of the GCPs and thus prone to errors.The area is covered by a large gravel pit with ongoing activity since the 1960s and by the Lakkasuo mire complex and finding multitemporal tie points was particularly difficult there.The mean L2 norm of the image residuals changed from 10.8 µm for full AT to 10.5 µm.
The XY shifts were small implying congruence between GCPs and DSO observations.

Aerial Triangulation Based on Direct Sensor Orientation Only
AT was solved such that the 22 full GCPs and 9 elevations points were treated as tie points with unknown coordinates in (2).Thus, the only exterior geometric support came from the DSO of the 50 images of flights 04402 and 04403.At the solution, the points had shifted with respect to the positions of the full AT (Fig. 11, Fig. 12).The 22 full GCPs, now free to move in AT, shifted an average of −0.07 m, 0.15 m and 0.02 m in X, Y and Z, respectively; the SDs of the shifts were 0.15 m, 0.12 m and 0.29 m.For example, the solution of point 88M5041 was 0.12 m above the true value and the RTK-GPS-positioned GCPs were an average of −0.10 m down from the true positions.These shifts indicated the presence of inconsistencies and correlated errors in the true elevations of the GCPs, assuming that the DSO observations and the calibrated focal length of the DSO images were reliable.The 279 tie points and 9 elevation points of full AT had shifted −0.06 m, 0.11 m and −0.07 m in X, Y and Z, respectively; the SDs of the shifts were 0.08 m, 0.07 m and 0.13 m.
The smallest shifts were observed in the coverage of flights 04402 and 04403 with DSO, where the results of this AT were most reliable.The points had shifted down in the area of the RTK-GPS points, which also coincides with the coverage of DSO flights 04402 and 04403 (Fig. 12).
As in full AT, the points with signals were computed photogrammetric coordinates using all image pairs with over 50% overlap (Section 4.2.2).The mean differences between the photogrammetric and true coordinates were less than 0.02 m and the SDs were slightly smaller than in full AT in Table 6.Thus, the geometry of the image pairs had not deteriorated, although no GCPs were used.

Geometric Match Between Lidar and Photogrammetric Observations
Manufactured targets with accurate geodetic XYZ observations were not available for assessment of the absolute accuracy of the lidar (cf.Csanyi et al. 2005).Thus, the relative mismatch between the photogrammetric observations and lidar points could only be estimated.The match in XY was examined in 4 target locations.One was a flat roof (Fig. 13), while the others were shores of ponds having open, low mire vegetation (Fig. 14).The XY match was estimated to be 0.5 m or better, based on the agreement of the superimposed lidar points and features in the images.
Vertical differences were checked using the images from 1997 (flight 97251) and 2004 (flight 04403) each with 2 different exterior orientations.In all, 43 presumably flat and smooth surfaces, mostly gravel-paved forest roads, were meas-ured 4 different photogrammetric elevations that were compared against lidar points in a circle with a radius of 1 m (Table 8).The measurement accuracy of the lidar was assumed to be app.0.1 m.The lidar points had echoed from 0.23 to 0.34 m below the photogrammetric points, whose  The lidar data provider reported that a height offset from 4 cm to 15 cm is possible because of the 'blind' processing, i.e. processing without control targets.The observed offsets are likely explained by the inaccurate elevations of the GCPs, especially that of point 88M5041, which fixed the geometry of the lidar points.The 43 target areas were located in the coverage of several image quadruplets, triplets and pairs, and photogrammetric measurement accuracy varies between and within such image sets because of orientation errors.Similarly, the lidar data can have errors that cannot be removed in the processing: strip adjustment in this case.

Examples on the Potential Use of the Oriented Images
The accuracy of the image orientation enabled stereointerpretation and multiple image matching, using epipolar constraints (Korpela 2004, p. 30).Fig. 15 illustrates how multitemporal images can be used for the determination of individual tree height growth and Fig. 16 illustrates monitoring of terrestrialisation.GCPs or DSO observations provided the exterior XYZ information needed in multitemporal AT in this study.DSO can consist of observations of projection centre coordinates and, in some cases, of the 3 camera rotations.DSO observations are insensitive to flying altitude and imaging geometry, e.g.image overlap, since no triangulation is needed.Here the projection center coordinates were observed with 0.1-m accuracy using DSO, which was presumably better than the accuracy of many of the GCPs, especially for the elevation.In AT a GCP observed with a 10-µm error in a 1:30000 image is off by 0.3 m on the ground, which exemplifies the advantage of DSO for high flying altitudes.DSO observations can have systematic errors that can only be controlled using GCPs (Cramer et al. 2000).Concerning DSO, the most interesting finding here was that the exterior geometry could be brought back to that of 1962 satisfactorily using only the DSO observations of 2004 images and the multitemporal tie points.
Tangential XYZ-systems are needed in photogrammetric projects, and it was assumed here that the terrestrial KKJ/N60 system with orthometric elevations approximates one in a small area.Satellite positioning, which is currently used in DSO and in geodesy, relies on accurate global 3D coordinate systems.If several coordinate systems are needed, accurate transformation between them is required to ensure high accuracy.Similarly, attention should be focused on the accuracy of the GCPs.Here, an error was made when the GCP set was densified using RTK-GPS measurements relative to an old geodetic point with an unreliable elevation.Similarly, the lidar points were fixed to a point with an unreliable elevation.I suggest that GCPs be postestablished and measured using accurate satellite positioning, and do not recommend the use of trigonometric geodetic points unless documentation reliably suggests their high accuracy.
Multitemporal tie points played a central role in the AT; they should be selected and measured carefully.Immobility over time is as essential as measurability in the images.The multitemporal tie points served in place of GCPs in the old images.Finding good multitemporal tie points is likely easier in an urban environment.Using very small-scale images over a long time interval can make it difficult to reliably identify immobile targets; which was experienced here in parts of the study area.Similarly, finding multitemporal tie points in large areas of intact forests can be difficult.In such difficult cases it is possible use 2D or 1D immobility constraints between tie points.For example, the same treetop may be observed at different time points as a tie point and the resulting XY coordinates are constrained to be to have a small separation in XY (with additional observation equations in A and observation weights in P in (2)).Similarly, tie points observed in planar areas may be constrained to have equal Z coordinates.3D immobility of the multitemporal tie points that was applied here is not necessarily required.
The image material in this study does not represent an average case in Finland.Large-scale images are rarely available except for in urban areas.The availability of large-scale images with large image overlaps undoubtedly aided in finding reliable multitemporal tie points.Otherwise, the study area represents the normal overall landscape structure in southern Finland.Repeating the procedures described here and extending the temporal coverage should, however, be possible in other areas.Automatic photogrammetric procedures were not used, apart from measurement of the fiducial marks in this study.Commercial photogrammetric workstations have algorithms for the automatic measurement of tie points.In a multitemporal image set, these could be utilized at least for finding the 'transient' tie points between consecutive frames.How well different commercial software packages are suited for and can support multiscale, multitemporal AT, remained largely unknown.Bundle block adjustment in this study was carried out using an experimental program written in Basic and Matlab was used for the numerical calculations.I attempted to use EspaBundle (EspaSystems) for AT, but its algorithms failed in solving the needed initial approximations.EspaBundle was instead used to verify the results of AT using a workflow with external files.
The project was based on the feasibilities of using multitemporal, archived image data together with modern ALS for precise retrospective analysis of forests.If lidar points are used for constructing an accurate digital terrain model, it is important that the photogrammetric observations and lidar-based terrain models remain in the same coordinate system over time.Preferably, the absolute geometric accuracy should be correct for all data, not only the relative match.The present study exemplified how the orientation accuracy of airborne lidar data can be evaluated against the photogrammetric data using natural and artificial targets.Comparison of the results of the three different versions of multi-temporal AT exemplified the difficulties inherent in the concept of absolute accuracy.The low, unknown accuracy of the GCP elevations, and their associated errors made it difficult to draw conclusions about the absolute geometric accuracy of the image or lidar data.The geometry of the multitemporal image block was also complex, because there were images at multiple scales, and some of the tie points were observed in as many as 87 images.The DSO observations were assumed to be accurate throughout the analysis, but this could not be verified in detail due to the errors in the GCPs.However, accuracy levels of below 0.5 m in 3D were observed for both the photogrammetric data in 1962-2004 and the lidar observations, which constituted satisfactory use of the material for precise retrospective analysis of trees and other objects.For comparison, the image registration method used in Miller et al. (2000) resulted in ± 1-5 m horizontal and vertical errors.In that study multitemporal aerial photographs were used for monitoring changes in coniferous canopies.
The possible applications in which integrated image and lidar data could be used are many.In photogrammetric single-tree forest inventory, it is possible to measure treetop positions in 3D, derive tree heights and height growth, measure crown dimensions, interpret tree species, and monitor removal and mortality of trees, all within the limits of discernibility of trees in aerial images (Talts 1977, Korpela 2004, Korpela and Tokola 2005).Other phenomena that could be monitored include: -Retrospective identification and timing of management operations and studying stand dynamics, -Changes in the crown cover of virgin or restored peatlands, -Rate of terrestrialization in paludified ponds, -Local, stand-specific terrain models using, automatic photogrammetric tools and old images in which the stand is seen as a clear-cut.
Archived photographs present massive amounts of data that are still largely stored on film.Archives can consist of hundreds of thousands of photographs.Digital storage and transfer of such vast amounts of data are challenging.It may be notably less expensive to store the films than to maintain a digital archive.However, it may be advisable to begin digitizing the oldest photographs to save the data stored therein.Retrospective analysis of forests is one form of usage for the data.

Fig. 1 .
Fig. 1.Different types of ground point used in the study area.The code is given for the geodetic points (triangles) by the National Land Survey of Finland.

Fig. 2 .
Fig. 2. Coverage of 288 aerial photographs from 1962 to 2004 per flight.The Indeces are drawn in the KKJ coordinate system with grid {2509 × 2521E, 6849 × 6866N}.The flight code, date of photography, film type, nominal scale, and camera and lens type are given in each index.NAG refers to a normal-angle camera with a 210-mm focal length.All other photographs were taken with a wide-angle lens having a 150-mm focal length.The first 2 digits in the flight code indicate the year of the photograph.

Fig. 3 .
Fig. 3.An illustration of the principle of multitemporal AT using bundle block adjustment.The position of the projection centre (X 0 , Y 0 , Z 0 ) and 3 rotations of the camera ω, φ and κ are the unknown exterior orientation parameters solved simultaneously for each image.There are tie points with unknown XYZ coordinates, and GCPs with known XYZ, XY or Z coordinates that are observed in the images.A tie point is either multitemporal or a 'transient point' observed in 2 consecutive images.DSO observations of (X 0 , Y 0 , Z 0 ) and (ω, φ, κ) may be available for the more recent images and can be utilized in solving the adjustment.Initial approximations are needed for the unknown camera parameters, coordinates of tie points and partially known GCPs.
Type of observation SD Image coordinates (x,y) 11 µm Projection centre (X 0 , Y 0 , Z 0 ) 0.1 m a) Attitude parameters ω and φ 0.000045 rad a) Attitude parameter κ 0.000072 rad a) a) Nominal RMS values of Applanix POS AV 510 DSO system used by the National Land Survey of Finland in flights 04402 and 04403.

Fig. 5 .
Fig. 5. Residual vectors of the Helmert transformation of 16 images of flight 97248 seen in xy film plane [m].The vectors are scaled by a factor of 500.Eight longer residual vectors are notably correlated between images and indicate stretch and/or shrinkage along the film diagonals.See text for more explanation.

Fig. 6 .
Fig. 6.Standard errors of the X coordinate of the tie points and elevation points [m].

Fig. 7 .
Fig. 7. Differences in the Z coordinate of GCPs and elevation points between true positions and photogrammetric positions [m].

Fig. 8 .
Fig. 8. Differences in the XY coordinates of GCPs between the true positions and photogrammetric positions [m].

Fig. 9 .
Fig. 9. XY shifts of the 279 tie points and 9 elevation points from the solutions for full AT to the solutions for AT based on the use of GCPs [m].Fig. 10.Z shifts of the 279 tie points from the solutions for full AT to the solutions for AT based on the use of GCPs [m].

Fig. 11 .
Fig. 11.XY shifts of the GCPs and tie points from the solutions for full AT (Section 4.2) to the solutions for AT based on the DSO observations only [m].Fig. 12. Z shifts of the object points from the solutions for full AT to the solutions for AT based on the DSO observations only [m].The rectangle depicts the coverage of flights with DSO.

Fig. 13 .
Fig. 13.An almost flat roof of 29.5 m × 17 m as seen in a near-nadir aerial image.Lidar points with Z = 165.8± 1.0 m are superimposed in white.Ten chimneys rise 60-70 cm from the roof and 7 lidar points with Z above 166 m are outlined.Their XY coordinates are within 0.5 m from the centres of the chimneys.The point spacing is about 1.2 m along the scan lines and 1.0 m and 1.4 m in the flight path direction (left-right).

Fig. 14 .
Fig. 14.A paludified pond in an aerial image of 43 m × 43 m.A sample of lidar points with Z = 168.04± 1.0 m are superimposed in white.The photography was done 2 weeks prior to the lidar.The area is under 2 overlapping lidar strips with a point spacing of about 1.4 m along the scan lines (top-down) and a more irregular spacing in the direction of the flight paths.The lidar did not register echoes from open water (dark areas) due to mirror reflection.Submergent and floating aquatic vegetation may have been registered by the lidar.The match in XY is better than 0.5 m, based on visual agreement of the shorelines in areas with high point density.

Fig. 15 .
Fig. 15.An example of tree height growth determination, 1985-2004.The treetop is positioned in 3D using the 1985 image pair and is at Z = 176.7 m.The XYZ position is moved up in 1-m steps for a total of 8 m, which is the total height growth between 1985 and 2004.The 8-m vertical movement is superimposed in the images, all of which represent oblique views, thus enabling the Z estimation from monocular images.The arrows depict treetop positions in 1989-2002.

Fig. 16 .
Fig. 16.An example of changes in shoreline extent of a pond under paludification in 1973-2004.Multitemporal image pairs were used for determining the elevation of the water level, which had varied from 141.6 m to 142.5 m.Open water is marked with W. A small inlet in the lower right corner was blocked up by vegetation between 1989 and 1995.The shoreline has shifted more than 5 m in places.

Table 1 .
Parameters of the airborne laser scanning run.

Table 2 .
A priori standard deviations of different types of observation in the hybrid bundle block adjustment.

Table 3 .
Mean statistics of the residuals (RMSEs) of the fiducial mark transformations[µm]. df = degrees of freedom.

Table 4 .
Mean statistics of standard errors (σ xk ) of the exterior orientation parameters [mrad] and [m].

Table 5 .
Mean statistics of standard errors (σ xk ) of tie point coordinates [m].

Table 6 .
SDs of the differences in the photogrammetric and true XYZ coordinates of the signalled GCPs observed in different image pair combinations with varying imaging geometry [m].

Table 7 .
Image residuals and maximal span for different types of tie point and GCP.