TECHNIQUES FOR COMPARING MULTI-TEMPORAL ARCHIVE AERIAL IMAGERY FOR GLACIER MONITORING WITH POOR GROUND CONTROL

: The measurement of geometric changes in Alpine glaciers is an essential aspect to assess their reaction to climate change effects. Archive aerial images may integrate valuable information to this purpose at times when other types of remotely sensed data are not available. The application of photogrammetric techniques such as Structure-from-Motion (SfM) and Multi-View-Stereo matching allows the extraction of dense point clouds to model the glacier environment. The comparison of multiple datasets requires to setup a stable reference system, a task that in archive photos is generally carried out by means of ground control points (GCPs). This paper would like to propose and assess some techniques to cope with the lack of ground control. Multitemporal SfM (MSfM) is presented and tested on a dataset including six different aerial blocks collected by means of analogue (PAN/RGB) and digital airborne cameras from 1967 to 2006. These images have been downloaded from the IGNF online repository and cover the area of the Val Veny (Brenva and Miage glaciers) in the Mount Blanc massif, at the border between Italy and France. Coupled with other solutions (i.e., extraction of GCPs from maps and DTMs and ICP co-registration of point clouds), MSfM has revealed as a suitable technique for coregistration of multiple photogrammetric blocks of aerial photos with minimum ground control. Some tests carried out in the case study area demonstrated that the integration of MSfM and ICP coregistration for refinement may significantly improve the comparison between multiple point clouds, which is a fundamental pre-requisite for the analysis of glacier changes over time.


INTRODUCTION
Monitoring glaciers to assess their reaction to climate change effects includes several aspects and requires multiple technologies and sensors.The evaluation of geometric changes plays a fundamental role to detect and to measure ice-mass variations.In recent years tailored UAV (Unmanned Aerial Vehicles) missions have been organized to investigate specific glaciers (see, e.g., Fugazza et al., 2018), while high-resolution satellite imagery provides information to reconstruct the geometry in the recent decades, to be continued in the future.Historical archive aerial photos (Poli et al., 2020) may provide valuable information to this purpose combined with up-to-date photogrammetric techniques, such as Structure-from-Motion (SfM) and Multi-View-Stereo (MVS) dense matching (Barazzetti et al., 2009).Point clouds of reconstructed 3D surface of glaciers may be used to detect thickness changes and height differences over time.This task can be accomplished by applying suitable techniques for computing distances between each pair of point clouds (Fey et al., 2017;Malekian et al., 2022).Under the technical point-of-view, one of the crucial points to cope with the alignment of point clouds derived from archive photos is how to deal with poor ground control (Tian et al., 2022).Ground control is fundamental to setup a stable reference system for the computation of coordinates to be compared within time.In addition, a good ground control also helps compensate for possible bias in the exterior orientation (EO).This task is usually implemented by introducing enough well distributed ground control points (GCPs), whose coordinates are measured by using differential GNSS techniques.On the other hand, the development of IMU/GNSS sensors may directly provide the EO of a photogrammetric block, but this solution cannot be used with archive photos dating back a few decades ago.Unfortunately, GCPs are often unavailable to be used with archive photos.Due to the lower resolution and data quality of these datasets, on one side, and the changes on the terrain (new constructions, vegetation growth, soil erosion, landslides, etc.), it is often difficult to measure new GCPs to be recognizable in photos collected in the past decades.During some research projects aiming at studying glacier changes within time, the authors have realized the need for alternative techniques to replace or to integrate the use of GCPs when comparing multitemporal datasets including photos from the archives.These techniques have been tried during a study regarding a group of Alpine glaciers in the Mount Blanc massif in the Italian Alps at the border with France.Six datasets of digitized analogue aerial photos provided by the National Geographic and Forestry Institute of France (IGNF) have been selected, downloaded, and used for photogrammetric processing.These datasets cover a time span of approximately 40 years from 1967 to 2006 and include images acquired with either analogue (panchromatic and RGB) and digital airborne cameras.While the change in ice thickness of these glaciers was relatively small until the mid-1990s, this study revealed an increasing reduction rate at the beginning of 21 st century.The proposed paper would like to describe the methodological approach for dealing with poor ground control.More scientific results on the case study can be found in Malekian et al. (2022;2023).

Overview
On the basis of the preliminary issues illustrated in Section 1, we consider a generic set of aerial blocks, which may come either from the archives and from more recent aerial missions.In general, historical datasets of archive photos have not been collected for investigation/monitoring of glaciers but consisted in typical photogrammetric blocks at scales used for topographic mapping projects (1:7,000 -1:30,000) and standard overlap (60%-70% along strip and 20%-30% side lap).Exception to this is the case of blocks flown over Antarctic and Arctic regions, which were collected with the purpose of scientific investigation (see, e.g., Fieber et al., 2018;Tian et al., 2022).While for block collected in recent years the ground control can be established by a set of GCPs or by IMU/GNSS data (if available), in the past the datasets may consist in the images only.
In the following, three different techniques will be illustrated to overcome the lack of ground control: 1. Multitemporal Structure-from-Motion (MSfM); 2. Extraction of GCPs from maps and DTMs; and 3. Surface coregistration based on ICP (Iterative Closest Point).
These techniques are not intended to be only use independently from one another, but they may be also cascaded in a structured procedure.

Multitemporal Structure-from-Motion
Structure-from-Motion Photogrammetry (SfM -see Granshaw, 2018) is commonly referred to as the entire automatic procedure for computing the EO of a block of images and to output a dense point cloud derived from the application of Multi-View-Stereo (MVS) matching techniques (see Eltner et al., 2016;James et al., 2019).While at the beginning SfM could be only applied to close-range blocks, successive development in processing and computing techniques extended its application to any types of images (Barazzetti et al., 2020).Today SfM can be used to process a photogrammetric block by using low-cost and open-source software packages running on a standard laptop.Cloud-computing resources are also popular to process large blocks without needing local hard resources.
When photogrammetric blocks must be compared to detect changes in the surveyed area, the standard approach would be to process each dataset independently from others.Each of them is georeferenced by means of a suitable set of GCPs to be included in the bundle block adjustment BBA as weighted pseudoobservations or to apply a similarity transformation to the derived point cloud.The adopted approach depends on the options offered by each software package.After image EO and MVS matching are completed, a dense point cloud is obtained to be compared with respect to other datasets.In this approach, GCPs may not be the same at different epochs, but they should all refer to the same reference system (see Barbarella et al., 2015).Direct orientation based on INS/GNSS integrated sensors is another solution to cope with georeferencing of an aerial photogrammetric block.This solution, anyway, cannot be available when using datasets recorded before this technique entered in the state-of-the-art, i.e., before late '90s.In many cases, GCPs cannot be used to align multitemporal datasets obtained from archive photos.This may be due to different reasons: the institution distributing the images did not collect GCPs, different reference systems and techniques adopted for their measurements, low quality.In such a case, an alternative procedure is to bundle together multitemporal blocks and to compute the image EO/camera calibration at the same time.This solution (referred to hereafter as multitemporal SfM -MSfM) presents the advantage that images from all epochs will be aligned together, as proposed by Feurer and Vinatier (2018) in their so-called "Time-SIFT" method.The knowledge of accurate GCPs becomes less relevant, because: (1) GCPs are not necessary for the alignment, since blocks are already aligned after SfM; (2) the availability of GCPs is useful only for global georeferencing, requiring a quality that may also be at lower accuracy (see Subsect. 2.3).Of course, if one of the adopted photogrammetric blocks is provided with accurate GCPs, this could be used for georeferencing all datasets.
Disadvantages and problems may be related to the different quality, resolution, and radiometric content (panchromatic vs RGB) of the aerial images, and changes in the content.This last aspect is very crucial when dealing with glaciers and, in general, the high mountain environment (Scaioni et al., 2018).For the sake of completeness, Di Rita et al. ( 2020) already suggested to apply MSfM to some datasets of UAV images recorded on an Alpine glacier for monitoring purpose.

Extraction of GCPs from maps and DTMs
If GCPs are not available, a possible solution is to use other types of geographic information to derive their coordinates.To do this, the planimetric component is usually split from the altimetric component.The former may be obtained from existing maps produced close to the acquisition time of the photogrammetric block to be georeferenced.In such a case, the accuracy of GCP planimetric position is usually in the order of ±10 m.The latter may be derived from digital terrain models (DTMs), that in general are more stable over time (at least in some areas) due to the minor changes in terrain topography (see Okolie and Smit, 2022).In such a case, the accuracy of GCP height is in the order of ±5 m.The estimated accuracy of GCPs derived from other geographic data may seem quite poor and thus not useful.On the other hand, in many cases this accuracy is enough to provide a first alignment between different photogrammetric datasets, to be refined later (see Subsect. 2.4).An example of this application is reported in Tian et al. (2022), where data from satellite photos collected in the '60s have been compared to recent satellite data.

Coregistration of point clouds using ICP
This approach allows to align a pair of point clouds based on algorithms for point sets registration.If more point clouds need to be aligned, as in the case study presented in Section 3, one is chosen as reference ("master") and all others ("slaves") are registered to it.The criteria for "master" selection may rely on the one with its own georeferencing, or the one with the best data quality, or the one with content less dissimilar from others.
The most popular technique to apply is Iterative Closest Point (ICP), see Pomerleau et al. (2013) for a review about this subject.ICP allows to compute a rigid-body transformation between a pair of point sets by starting from an initial approximate alignment.This can be obtained from manual picking of at least three corresponding points, or by using independent georeferencing of both datasets.In fact, ICP is generally also applied to refine the coregistration already obtained by using other techniques described in Subsections 2.2 and 2.3.
Since in this study we would like to apply change-detection techniques (Lindenbergh and Pietrzyk, 2015) to understand the glacier evolution from multiple point clouds derived from photogrammetric processing of historical aerial photos, we should consider that some areas might undergo important modifications over time.The ICP-based coregistration should only include some manually stable areas to avoid biased outputs.This task is usually accomplished by manually extracting these stable regions.Otherwise, a robust registration technique may be applied to automatically discriminate between areas with and without changes (see Wujanz et al., 2016).Stable regions will be used to compute the rigid-body transformation to align the whole "slave" point cloud to the "master".After coregistration is computed, it is important to assess the quality of this result by checking residuals.If this test is positive, the same coregistration parameters obtained for stable areas can be applied to the remaining part of the "slave" point cloud.
After point cloud coregistration, the comparison can be operated (Tavakoli et al., 2021).Several techniques exist to this purpose, either on the direct basis of the point clouds, and the use of derived models such as meshes (e.g., Triangulated Irregular Networks) or grid raster digital elevation models (DEMs).
In the final analysis of changes, it is important to consider the of uncertainty in point cloud registration.This option is possible in some advance algorithms, such as M3C2 (Lague et al., 2013).

Dataset presentation
The glaciers under investigation belong to the Mount Blanc massif in the western sector of Italian Alps at the border with France.These glaciers flow in the Val Veny, as can be seen in Figure 1, which also reports the locations of one out of six photogrammetric blocks from aerial missions operated in 1967-2006.As shown in Table 1, the oldest datasets were collected by means of analogue cameras, while the latest could already exploit digital aerial cameras (DAC) that are today the state-ofthe-art technology in topographic aerial mapping.Datasets are available as panchromatic or RGB colour images (see Table 1).
Images have been downloaded in digital form through the geoportal "IGN -Remontez le temps" of the IGNF (https://remonterletemps.ign.fr/).This archive also includes other older aerial missions over the same area (i.e., in 1952 and 1958) that would make the analysis even more interesting.The image quality and the cameras adopted to collect these datasets was more difficult to be processed.For this motivation, we decided not to use them in this study, but to postpone their use in the future (see Sect. 4).Photos recorded in analogue format have been scanned with photogrammetric scanners at pixel size of 20 µm.Some datasets are also provided with a calibration certificate containing information on the camera and lens (Kraus, 2008), and with camera locations recorded by onboard GNSS sensors (1996, 2000 and 2001).No GCPs are available.In Malekian et al. ( 2023) some considerations about image quality are reported.

Application of multitemporal SfM
Multitemporal SfM (MSfM) has been applied to the image blocks listed in Table 1 to understand the effectiveness of this procedure to align aerial datasets featuring multiple characteristics (e.g., type of camera, radiometric content, flight geometry, photo scale) and different acquisition time.Both groups of differences are particularly challenging, but the changes on the ground are significant, since they account for the glacier retreat, different snow cover, soil erosion, anthropogenic modification over time.
In some previous publications (Malekian et al., 2022;2023), all datasets listed in Table 1 were independently oriented using SfM.Some sets of GCPs were obtained as described in Subsection 2.3.The ICP alignment between point clouds derived from MVS dense matching was operated at a later stage by using the approach presented in Subsection 2.4.
Here we have followed a different approach, which is based on georeferencing only one dataset (2000) by using 7 GCPs, as presented in previous publications (Fig. 1).Once georeferenced, Dataset 2000 has been used as geographic reference to align all other photogrammetric blocks to it and to evaluate changes.Photogrammetric processing has been carried out in Agisoft Metashape Professional ® ver.1.7.5 (in the following AMP), which represents a standard package available at low-cost but with the capability to process datasets of large aerial images as the ones employed here.
The process for the orientation of selected aerial blocks is illustrated in the following paragraphs.Before applying SfM, images have been pre-processed to improve the image quality as illustrated in Malekian et al. (2023).et al. (2023).Their values have been retained sufficient considering the type and quality of adopted images (see Aguilar et al., 2013).Another important consideration concerns the inner orientation and camera calibration.Due to the lack of calibration certificates for some datasets, self-calibration based on standard Brown model has been included in the BBA (Luhmann et al., 2014).This option has been selected with the purpose of defining a common processing workflow for all types of aerial blocks.In addition, this solution is supported by the presence of a sufficient number of images per each block (see Luhmann et al., 2016) and a suitable block geometry in terms of strip overlap.Approximate focal lens and pixel size have been provided.The unique block acquired by DAC (Dataset 2006) has not been processed with self-calibration, since images have been delivered as distortion free and a camera calibration certificate was retrieved upon request to the IGNF.The most datasets were collected using different high-end analogue cameras (Table 1), thus they reported the outer frame with fiducial marks and instrument recordings.According to the experiments already reported in Malekian et al. (2023), we preferred to mask the outer camera frame to exclude the corresponding area from the extraction of key-points and tie points to input in the BBA.This also avoided the recognition and measurement of fiducial marks (Forlani et al., 1999).A simple rectangular mask has been applied to all the images of the same dataset.Figure 2 shows an example of a such mask applied to an image from Dataset 1967.This image also reported tie points after (blue points) and before (blue and white points) internal blunder rejection.

Multitemporal SfM:
In a second stage, multiple blocks have been progressively included in a multitemporal SfM (MSfM) project in AMP.Dataset 2000 has been used to establish the datum thanks to the presence of 7 GCPs.The measurement of these points has been also extended to some images from other datasets where they could be clearly identified and manually measured (Fig. 3).This solution has led to an average number of 10 rays per GCP, which increased the reliability of corresponding observations in the BBA.
To setup the MSfM, a common block including photos from all datasets has been first created by aligning and merging individual projects in AMP.This software allows to import them as "chunks" (i.e., subprojects) in a new project.By using the extracted tie points in each block and by looking at intrablock correspondences, all "chunks" could be aligned with respect to one of them selected as reference (in this case: Dataset 2000).After merging all blocks in a unique project, they have been all aligned w.r.t. the same datum.At this point, a BBA has been run again to improve the image orientation, which has received as input 103,259 tie points (average 2.92 rays per point).Also in this case, processing has been computed at "Highest" quality level.According to the poor quality of GCP measurement on the images (see Fig. 3) and the source of 3D coordinates (from existing maps and a DTM), we decided to assign a low accuracy (GCP=10 m) not to influence the intrinsic accuracy of estimated EO in the BBA, but only to guarantee a general absolute common georeferencing to all blocks.This decision is different from the one adopted in previous publications, where a higher accuracy was selected for GCP object coordinates, resulting in smaller residuals.

Dense cloud generation:
After the EO estimation based on MSfM, a dense point cloud has been generated per each original dataset in AMP.The quality level "High" has been selected, as trade-off between final quality and processing time.
Point clouds have been cleaned and segmented by using a common Region-of-Interest (RoI).The size of dense points clouds has ranged from 26 million points (Dataset 1967) to 56 million points (Datasets 1996 and2006).Dense point clouds have been exported from AMP in LAS interchange format and imported in the open-source software CloudCompare (CC) Ver.2.12.4 "Kyiv" (www.Cloudcompare.org)for successive processing operations.In CC each point cloud has been processed first to remove duplicated points by selecting a minimum distance between closest points larger than 20 cm.Secondly, a filter to remove residual isolated points was applied.The quality assessment has been performed in CC base on five stable areas, as presented in next paragraph.After that, the same areas have been used for computing a rigid body transformation to minimize distances between different epochs based on ICP registration.Then each set of transformation parameters has been applied to the remaining part of each "slave" point clouds to refine the alignment with respect to the "master" obtained from MSfM.In Figure 4, point clouds obtained from Datasets 1967 and 1988 are compared by using M3C2.Differences are shown to highlight glacier changes.

Quality assessment:
The final empirical quality assessment has been based on the comparison of dense point clouds inside the RoI.Due to the presence of several and large changes between different epochs, the analysis has been limited to five stable areas, as in previous publications.Points within five stable windows have been segmented and compared by using M3C2 (Multi-scale Model to Model -see Lague et al., 2013) technique, that was applied to detect changes between point cloud pairs (see also Tavakoli et al., 2021).M3C2 is based on the computation of local normals, which are used to tailor in each location the direction to detect changes.In addition, M3C2 is resistant to changes in point density and point cloud noise and allows to exploit information on the point clouds' accuracy (if available) to distinguish between real displacements and noise (James et al., 2017).All points in stable areas per each epoch have been merged together and compared to other epochs.Table 3 reports the Root Mean Square Errors (RMSE) between different epochs, where the older point cloud has been used as "master" and the other as "slave" in the M3C2 comparison.In Table 4 we reported more statistics about the comparisons between consecutive epochs only.

Discussion
The core of the methodology for the alignment of multiple photogrammetric blocks proposed here consists in the MSfM as alternative to the independent orientation of single blocks.MSfM has demonstrated to be more convenient from the operational point-of-view, requiring a unique set of GCPs with the only purpose of global georeferencing.
As far as the quality of block co-registration is concerned, the analysis of residuals on stable areas between successive epochs (Table 4) shows that, in terms of RMSE, the application of MSfM has resulted in significant improvements w.r.t. the independent SfM on single datasets.In the case of older datasets, the improvement in terms of RMSE has ranged from approx.30% to slightly less than 50%.With more recent datasets, the improvement has reached approx.80%.These results are motivated by the good performance of MSfM in terms of extracted tie points, as shown in Table 2.
It is important to notice that MSfM works well also for the integration of photogrammetric blocks featuring big differences in terms of sensing technology (analogue vs digital), radiometric content (panchromatic vs RGB), and relevant topographic changes.
By looking at the components contributing to RMSE in Table 4, standard deviations play the most relevant role.This demonstrates that data dispersion counts more than bias after MSfM.This is confirmed by results after the ICP refinement reported in the rightmost column of Table 4. Mean differences are reduced in terms of a few centimetres (max mean difference is -22 cm between Datasets 1996 and 2000).
Table 3 shows the RMSE after M3C2 comparison of points in stable areas for all possible combinations of epochs.General results are like the ones obtained from consecutive datasets reported in Table 4. Also in this case, larger values of RMSEs have been obtained when older datasets have been compared.This is motivated by the improved quality of aerial photos that could be obtained in the 25 years.Dataset 1996 has revealed to be more problematic, probably due to the lower number of images and the block geometry.Results shows that this technique is highly suitable for the alignment of selected datasets, which also included quite big challenges in terms of camera technology, radiometric content, and relevant topographic changes.This outcome suggests that MSfM can be generally applied to other case studies with similar performances.In addition, a further improvement could be achieved from the integration between MSfM and the refinement of point-cloud alignment based on ICP computed in stable areas.Here stable areas, which can be also exploited for quality assessment, have been manually selected.One of the possible future developments should be their automatic selection (see Wujanz et al., 2016).

CONCLUSIONS AND FUTURE DEVELOPMENTS
Other interesting evolutions of the presented methodology is the application of new image matching techniques based on Artificial Intelligence (e.g., Deep Learning), see Morelli et al. (2022) andMildenhall et al. (2022).
Another important aspect to develop is about the comparison of point clouds based on M3C2 technique.Recently new developments of this method were published (Winiwarter et al., 2021;Zahs et al., 2022), also considering the possibility of the contemporary analysis of multiple datasets at the same time, instead of pairwise comparison (Anders et al., 2021).

Figure 2 .
Figure 2. Satellite view of the area under investigation including: glacier outlines (see legend), ground control points (GCPs), and locations of aerial photos from the Dataset 2000 (acknowledgements to Google Earth ® for the background satellite image).
3.2.1.Independent SfM of single blocks: In a first stage each block has been independently processed without considering overlap with other datasets.GCPs have been measured only in Dataset 2000 and used to setup the Italian geodetic datum (RDN 2008.0 / UTM32N) in the BBA.In other datasets, a minimum constraint datum has been used (see Luhmann et al., 2014) to compute the EO in a first stage.Image orientation has been computed based on SfM implemented in AMP: processing has always been computed at "Highest" quality level setup, that corresponds on the use of the full-resolution imagery.The image size in the considered blocks ranged from approx.16 Mpixels to 136 Mpixels.Almost all photos in each block could be oriented apart from a few images (less than 3% of the total number of available images).Residuals in image and object space (the latter only for Dataset 2000) are reported in Malekian

Figure 2 .
Figure 2. Example of an image from Dataset 1967 (panchromatic film recorded by analogue aerial camera Leica RC10) masked to avoid the interference of the outer frame during Structure-from-Motion (SfM).Overlaid points represent extracted tie points, while blue points are the ones which passed the internal blunder rejection (acknowledgements to IGNF ® for the aerial photo).

Figure 3 .
Figure 3. Images of a GCP tracked in multiple photos across different epochs.

Figure 4 .
Figure 4. Result of the comparison based on M3C2 technique between Datasets 1967 and 1988 (panchromatic) over the terminus of Miage Glacier in Val Veny (distances in [m]).In the upper part the point clouds obtained from the application of MSfM are shown.After ICP refinement of coregistration, point clouds have been compared by means of M3C2 (see the map of changes in the lower part).

Table 2 .
Root Mean Square Errors (RMSE) of GCPs and other statistical parameters after bundle block adjustment of Multitemporal SfM (MSfM); results refer to sub-blocks (one per epoch) and to the multitemporal block as a whole.

Table 3 .
RMSE after M3C2 comparison of point clouds in five stable areas from different epochs (results in [m]).

Table 4 .
(Malekian et al., 2023)e application of Multitemporal Structure-from-Motion (MSfM) as alternative technique for the alignment of multiple datasets of archive aerial photos collected at different times.The target application consisted in monitoring changes in Alpine glaciers as consequence of climate change.MSfM was tested in a case study (Val Veny in the Mount Blanc massif, Italy) where multiple datasets of aerial photos could be retrieved from the online archive repository of National Geographic and Forestry Institute of France (IGNF).Statistics after comparison of points in stable areas according to three different approaches: single SfM based on independent GCPs per each dataset(Malekian et al., 2023); MSfM with GCPs from Dataset 2000; MSfM followed by ICP refinement.