ARCHAEOLOGICAL FEATURE DETECTION FROM ARCHIVE AERIAL PHOTOGRAPHY WITH A SFM-MVS AND IMAGE ENHANCEMENT PIPELINE

Understanding and protecting cultural heritage involves the detection and long-term documentation of archaeological remains alongside the spatio-temporal analysis of their landscape evolution. Archive aerial photography can illuminate traces of ancient features which typically appear with different brightness values from their surrounding environment, but are not always well defined. This research investigates the implementation of the Structure-from-Motion Multi-View Stereo image matching approach with an image enhancement algorithm to derive three epochs of orthomosaics and digital surface models from visible and near infrared historic aerial photography. The enhancement algorithm uses decorrelation stretching to improve the contrast of the orthomosaics so as archaeological features are better detected. Results include 2D / 3D locations of detected archaeological traces stored into a geodatabase for further archaeological interpretation and correlation with benchmark observations. The study also discusses the merits and difficulties of the process involved. This research is based on a European-wide project, entitled “Cultural Heritage Through Time”, and the case study research was carried out as a component of the project in the UK.


INTRODUCTION 1.1 Background
Detection, interpretation, documentation and monitoring of archaeological features are fundamental for understanding the historic environment, as well as maintaining and protecting cultural heritage.Archaeological features commonly survive as subsurface remains (e.g.ancient roads, ditches, ruined buildings) which can be evidenced as crop/soil marks and/or earthworks on historical photographs (Evans and Jones, 1977).Such features can be vulnerable to landscape dynamics due to natural and cultural processes such as weather and climate change, pollution, agriculture, urbanization and other human activities ranging from tourism to war.
Excavations and geophysical surveys have long been utilised for the identification and interpretation of archaeological remains (Gheyle et al., 2016;Kvamme, 2003).In addition, a wide range of geoinformatics technologies, including photogrammetry, remote sensing, laser scanning, Global Navigation Satellite Systems (GNSS), topographic surveys etc., are well established techniques for archaeological documentation and monitoring (Xiao et al., 2018).The advance of sensor miniaturisation alongside unmanned aerial vehicle (UAV) technology has also facilitated multi-temporal 3D cultural landscape modelling.For example, recent studies have investigated the derivation of vegetation indices from multispectral imagery (Agapiou et al., 2017) and surface morphological attributes from laser scanned digital surface models (DSMs) (Toumazet et al., 2017) to automatically detect and record archaeological structures.
Modern technology aside, archaeologists have long used historic aerial photography for investigations (Evans and Jones, 1977;Kaimaris et al., 2012;Verhoeven, 2012).Archaeological features on photographs usually appear with different chromatic indicators and textures.Their visibility (ground response) mainly depends upon soil moisture conditions, crop regime, vegetation type, weather conditions, light and time of the day as well as the state of the features themselves (e.g.excavated, buried etc.).To distinguish archaeological features on historical photographs can be a particularly challenging task as the images often lack accompanying information about conditions prevailing at the time they were captured.
The study presented here proposes an image enhancement pipeline to detect and record the 3D/2D locations of possible archaeological traces from visible/near infrared (NIR) archive aerial photography.This pipeline, implemented in conjunction with the Structure-from-Motion (SfM) and Multi-View Stereo (MVS) image matching approach (James et al., 2017), can ultimately document a time-series of 3D archaeological features in support of understanding cultural heritage loss.The research has been conducted as part of the "Cultural Heritage Through Time" project (CHT2, 2017), supported by the European Joint Programming Initiative for Cultural Heritage, which aims to merge heterogeneous information and expertise to deliver enhanced four-dimensional (4D) digital products for heritage sites (Fieber et al., 2017).The proposed pipeline is tested at Corbridge Roman Town (Historic England, 2017), one of the three CHT2 UK case study sites located on the Hadrian's Wall landscape (Hadrian's Wall, 2017).

STUDY AREA
The study area of Corbridge Roman Town (54° 58′ 41.11″ N, 2° 1′ 45.36″ W) is located westwards from the edge of the modern village of Corbridge, Northumberland, UK (Figure 1).It was the most northerly urban settlement in the Roman Empire.Archaeological investigations have revealed that the Roman presence in the area evolved over the first four centuries AD, with settlement centring on what is now the current English Heritage Guardianship Area c. AD 85 (Bishop and Dore, 1989;English Heritage, 2018b;Haynes and Turner, 2017).Both military and civilian activities took place across the settlement.Remains of legionary compounds, civilian buildings, gravel quarries, bridge abutment, roads, cemetery spaces and funerary sites have been observed (English Heritage, 2018b).Of particular importance to the history of Roman settlement in this area is the intersection of Roman roads running E-W, the 'Stanegate', and N-S, the 'Dere Street', as seen in Figure 1.The 'Stanegate' road extended the Roman Town to the west and the 'Dere Street' provided a crossing over the River Tyne to the south (English Heritage, 2018a;Haynes and Turner, 2017).Figure 1 also shows the possible location of a second Century AD Roman mausoleum base structure (Region C) which was unearthed by Gillam and Daniels (1961).
Today the only part of the study area where archaeological features are visible and maintained is Region A, as shown in blue in Figure 1.Region A is open to the public under English Heritage management.Archaeological excavations in this area commenced in the early 1900s (English Heritage, 2018b;Haynes and Turner, 2017).It is important to note that the ancient Roman settlement extends well beyond the managed site in all directions.
Current research (Haynes and Turner, 2017) is attempting to characterise the pattern of broader settlement, but for this to be successfully achieved a better understanding of landscape evolution is essential.Reaching this understanding is a significant challenge, partly because the study area is subject to fluvial flood hazard from the River Tyne to the south as well as extensive ploughing activities.Unexcavated archaeological features have been disturbed by these activities, as already reported in the late nineteenth Century (Forster, 1881;Haynes and Turner, 2017).

METHODOLOGY
The methodological workflow consists of two main stages.The first stage involves multi-epoch, archive aerial photography collation alongside 3D surface model and orthomosaic reconstruction using the SfM-MVS approach, as described in Fieber et al. (2017).The second stage comprises an image enhancement pipeline to improve the recognition of archaeological features.
Multi-temporal SfM-MVS orthomosaics, generated in stage one, together with rectified aerial photographs, constitute the input datasets for this pipeline.These datasets are then subject to two processing steps, namely decorrelation contrast stretching and noise filtering.
Decorrelation stretching improves the colour contrast by minimising the correlation from the image bands in a similar way to principal component analysis (Campbell, 1996).The decorrelation stretching technique has long been successfully applied to visible (Lo Curzio and Magliulo, 2010) and multispectral imagery (Gillespie, 1992).Here, it is applied to three-band visible and NIR imagery.A normalised 5x5 pixel average convolution filter smooths the three decorrelated bands to remove speckle noise.Then the decorrelated band, the one which reveals the most significant information related to archaeological traces (Green and Red bands for visible and NIR imagery, respectively), is utilised to support their 2D digitisation.The 2D digitised location and elevation of the archaeological features, extracted from the SfM-MVS 3D surface models, can be recorded into a geodatabase for archaeological interpretation and cross-validation with benchmark geophysical observations.

Datasets
Three epochs of aerial photographic datasets, acquired in 1984, 2006 and 2016, were utilised (Table 1).Aerial photographs from 1984 (captured using a Hasselblad H5D-60) and 2006 (captured using a Canon EOS-1Ds Mark II) were digitised at 2000 dpi resolution (Fieber et al., 2017) and supplied by the Historic England Archive.It should be noted that digitisation was not performed using a photogrammetric scanner, and the resultant imagery is therefore likely subject to distortions -e.g.Thomas et al. (1995).The modern aerial dataset from 2016 was provided by Historic England in digital format, and included oblique imagery captured using a NIKON D3X and a NIKON D800E.Of all regions shown in Figure 1, only Region C was not visible in the 1984 black and white (bw) datasets (Table 1 and Figure 2c  The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLII-2, 2018 ISPRS TC II Mid-term Symposium "Towards Photogrammetry 2020", 4-7 June 2018, Riva del Garda, Italy

SfM-MVS and georeferencing
The first stage of the methodological workflow, SfM-MVS processing, was undertaken using Agisoft PhotoScan (version 1.4.1)(PhotoScan, 2016).To ensure that SfM-MVS products from all collected datasets were accurately georeferenced into a common fixed reference frame (Ordnance Survey (OS) Great Britain 1936, OSGB36), 40 ground control points (GCPs) were surveyed in October 2017 using Real Time Kinematic (RTK) GNSS.The number of GCPs used per epoch as control information within the SfM-MVS workflow is reported in Table 1.GCPs could not always be identified precisely in archival and recent imagery (e.g.corner of a traffic island).In addition, the study area includes some inaccessible private land which did not allow an optimal distribution of GCPs.For these main reasons poor georeferencing results were delivered.The 1984 NIR orthomosaic could not be reconstructed since fewer than three aerial photographs existed for that epoch (Table 1).Therefore, the two single NIR images (one over Region B and another one over Region C as seen in Figure 1) were georeferenced with respect to the 2016 orthomosaic using five common points, delivering a 0.70 m average planimetric RMSE.

Data co-registration
In order to generate a spatially consistent 3D time-series of archaeological datasets, a rigorous co-registration from epoch-toepoch was performed, prior to implementing the image enhancement pipeline.The co-registration was applied to the three dense point clouds only over Region A (Figure 1), since adjacent fields were subjected to environmental change (e.g.crop growth) that would adversely affect co-registration results.The 1984 and 2006 dense point clouds were co-registered with respect to the reference 2016 dense point cloud.Buildings and vegetation were also excluded from all datasets.A sevenparameter Helmert transformation was applied using the iterative closest point (ICP) routine as implemented in the OPALS software (Pfeifer et al., 2014).To verify the ICP performance, statistics of the cloud-to-cloud differences over Region A were calculated with the aid of the M3C2 algorithm (Lague et al., 2013), as listed in  Not only did the implementation of ICP clearly minimise the point-to-point differences (Table 2) but, especially for the 2016-1984 co-registration, it also removed an apparent systematic tilt error, as evidenced in Figure 3.The Helmert transformation parameters were then applied to the dense point clouds initially reconstructed from the 2006 and 1984 datasets only over Region B (Figure 1).The co-registration was not applied to the dense point clouds over Region C, and as such it did not further improve the georeferencing result.The transformed dense point clouds were reimported to Agisoft PhotoScan and the MVS workflow was repeated without the necessity of the SfM and georeferencing steps.This procedure finally resulted in the reconstruction of three consistent dense point clouds, DSMs and orthomosaics.From an average 0.10 m ground sampling distance (GSD), as reported in Table 1, DSMs were generated at each epoch with an average 0.20 m spatial resolution.As is clearly evidenced from both sub-regions in Figure 4 and Figure 5, the 2016 imagery does not reveal any archaeological traces.Whilst is possible that the archaeological structures might have been more deeply buried or damaged agriculture and/or erosion, key points to consider are the day/time and climate conditions of the year the image was acquired as these can affect the crop height and growth.

RESULTS AND DISCUSSION
The 2016-2006 negative elevation differences over Region B (Figure 4c) and Region C (Figure 5c) indicate crop change.The most likely explanation for this is that the 2006 photographs were acquired in July during the growing season (Table 1).Specifically, a negative elevation change of approximately -1.0 m to -0.6 m is observed west and north of Region B (Figure 4c) as well as along the field tramlines over Region C (Figure 5c).By contrast, a positive elevation change of +0.6 m to +1.0 m is apparent south of Region A towards the River Tyne (Figure 4c), while the ±0.2 m minimal change over Region A lies within the estimated maximum RMSE after ICP implementation (Table 2).As a result, it can be assumed that the 2016-2006 positive elevation differences represent changes in crop height.However, it cannot be directly deducted whether the negative changes have been solely caused by geomorphological/fluvial dynamics or anthropogenic activities.
To estimate the 2006-1984 actual landscape change is also a significant challenge.After ICP co-registration, a ±0.2 m elevation difference is observed over Region A with the exception of greater positive change over buildings and dense vegetation (see (i) in Figure 4).The ±0.2 m complies with the maximum RMSE in Table 2.However, an apparent systematic error still remains and is partly indicated as a positive elevation change exceeding 1.0 m NW of Region B. It is likely that a combination of different factors have caused this error, such as poor imaging network geometry, unreliable modelling of the camera lens distortion, suboptimal GCP distribution, distortion caused by the scanned photography, non-photogrammetric nature of archive photography and noise from vegetation variations.Most of these are well known error sources in SfM-MVS derived products, as demonstrated in recent studies (Eltner et al., 2016;James et al., 2017).According to these studies the inclusion of oblique imagery in the SfM-MVS workflow could significantly minimise possible distortions in the derived DSMs.To that end, the 1984 vertical datasets, as seen in Figure 2c, proved to be insufficient for the quantification of real 3D landscape change.
To reliably estimate surface change and alleviate erroneous deformations in elevation differences, observations other than aerial photography could also be utilised, such as laser scanning observations.Airborne laser scanning outperforms photogrammetric products from aerial imagery in estimating exact changes in ground level/soil depth, as laser scanning can capture underlying terrain.Such multi-epoch observations could further support archaeological investigations by simultaneously separating the twofold cause of landscape change and identifying buried ancient features.
Besides the difficulties in 3D landscape evolution modelling, the historical aerial datasets provided useful semantic information advantageous for the 2D/3D archaeological feature documentation.Both visible and NIR imagery illuminated linear features which are related to ancient Roman roads and the mausoleum base.However, the 1984 NIR dataset provided finer details of linear structures across the Roman roads and south of Region A, compared to the 2006 visible dataset (Figure 4a).This can be also attributed to the possibility that the cropland was possibly set aside from production in 1984 as no field tramlines are observed in 1984 imagery.Therefore, archaeological features were not as obscured by farming activities, as in the case of 2006 (Figure 4b).Additional historical evidences should be considered to back up the aforementioned hypothesis.
In the case of the 2006 visible orthomosaic over the mausoleum base in Region C, linear traces were blended with the crop characteristics (Figure 5a).Analysis of that year's orthomosaic illustrated high band-to-band correlation, as seen in Figure 6a.This explains the difficulty in identifying the structural traces on the 2006 visible, as opposed to the 1984 NIR, imagery (Figure 5a).After the band-to-band correlation reduction and noise filtering, structural traces were better distinguished from the background cropland in both the 1984 and 2006 datasets (Figure 5b).Image enhancement mostly improved the archaeological feature identification over Region C as features over Region A were already quite well defined in both 2006 visible and 1984 NIR orthomosaics (Figure 4a).Moreover, the image enhancement aided in detecting additional features in the 2006 orthomosaic (see (ii) in Figure 4), which are not directly related to archaeological remains.It is speculated that these features, seen as linear traces with NS direction, may represent sediment deposition, caused by fluvial flood hazard from the River Tyne.This can also be explained by the observed positive elevation change (Figure 4c).Further analysis of flow accumulation and direction with the use of GIS tools could potentially support this statement.
Vegetation indices can constitute an alternative to decorrelation stretching for crop mark detection (Agapiou et al., 2017).A Red-Green difference vegetation index (DVI) was derived from the 1984 and 2006 orthomosaics as it is the clearest index to employ.Moreover, it was not possible to derive the 1984 normalised difference index (NDVI) as the pixel brightness values of the aerial orthomosaic could not be converted into spectral reflectance without ground spectral reference datasets.DVI results over Region C are shown in Figure 7a and b.The DVI outputs are very similar to the enhanced orthomosaics in Figure 5b.However, the decorrelated 1984 orthomosaic provided an output with higher contrast compared to the 1984 DVI (Figure 7a versus Figure 5b).Surface morphological attributes (e.g.shaded relief, local relief, slope, curvature, openness etc.) can distinguish certain geomorphological characteristics of archaeological structures separating them from their surrounding environment (Toumazet et al., 2017).For example, Toumazet et al. (2017) modelled archaeological elements based on their specific ring-shape as detected on a local relief.Here, Figure 7c shows the morphological attribute of shaded relied as derived from the 2006 DSM.The shaded relief image does not illustrate any linear traces of the mausoleum base structure and neither do the aforementioned morphological attributes.This verifies that the features detected on the 2006 orthomosaic are completely buried under the crop/grass land, which also applies to the detected features of Region B. The depth of the buried features can only be estimated with subsurface observations.Long-term monitoring of farming activities across the study area with the use of advanced UAV imagery would potentially aid the geodatabase maintenance of unexcavated 3D archaeological features.
To finally assess the accuracy of the research outcome and derive quantitative metrics, benchmark observations are required.Such observations can be obtained with geophysical surveys (e.g.ground penetrating radar and gradiometry).Such observations are being pursued as part of further ongoing research, as described in Haynes and Turner (2017).Geophysics surveying results in an independent set of 2D digitised archaeological features that can be used to verify the position of the digitised features of the presented workflow.Regarding the third dimension, multi-epoch DSMs generated with airborne laser scanning can also complement the SfM-MVS derived DSMs from archival imagery.

CONCLUSIONS AND FUTURE WORK
This paper has investigated the detection and documentation of 2D/3D archaeological features from archival aerial datasets in both visible and NIR wavelengths.Tests have demonstrated that SfM-MVS, in conjunction with an image enhancement, pipeline can potentially reveal cultural heritage lost assets.The paper has outlined some of the difficulties in post-processing archival photographs for deriving DSMs.Distortion was observed mostly when only vertical imagery was included in the SfM-MVS process.Co-registration, as implemented with ICP, partly reduced the distortion and improved the epoch-to-epoch consistency of the SfM-MVS outputs.The research has also presented the potential that archival photography can provide in identifying significant information related to unexcavated archaeology.The decorrelation stretching algorithm enhanced the contrast of visible imagery and enabled linear buried structures to be better illuminated.In order to comprehensively understand the landscape evolution of Corbridge Roman settlement, in addition to information derived from aerial historical evidence, other factors should be considered.These include the ongoing geomorphological dynamic processes and farming activities in the surrounding environment.
Ongoing research is investigating different algorithmic approaches to process black and white imagery.This will extend the presented time-series of historic datasets.Detection of archaeological features from geophysical sub-surface observations, together with relevant archaeological interpretation, are scheduled for 2018.The digitised features derived from all heterogeneous observations will also be compared against those features depicted on the historical map in Figure 1.Tests will be conducted at all three CHT2 UK case study sites and the results will ultimately support further archaeological investigations and cultural heritage conservation across the Hadrian's Wall landscape.

Figure 2
Figure2shows three perspective views of the reconstructed 3D surface models together with the estimated positions of camera exposure stations and GCPs.Average root mean square errors (RMSEs) at ground control points of 0.40 m in planimetry and 0.10 m in elevation were delivered for the 2016 dataset.For the 2006 dataset average RMSEs at GCPs of 0.40 m in planimetry and 0.14 m in elevation were delivered.RMSEs at GCPs of 1.07 m in planimetry and 1.92 m in elevation were computed for the 1984 dataset.GCPs could not always be identified precisely in archival and recent imagery (e.g.corner of a traffic island).In addition, the study area includes some inaccessible private land which did not allow an optimal distribution of GCPs.For these main reasons poor georeferencing results were delivered.The 1984 NIR orthomosaic could not be reconstructed since fewer than three aerial photographs existed for that epoch (Table1).Therefore, the two single NIR images (one over Region B and another one over Region C as seen in Figure1) were georeferenced with respect to the 2016 orthomosaic using five common points, delivering a 0.70 m average planimetric RMSE.

Figure 4a and
Figure 4a and Figure 5a show the reconstructed orthomosaics post SfM-MVS implementation and co-registration.The orthomosaics, enhanced after decorrelation stretching and noise filtering, are illustrated in Figure 4b and Figure 5b.Archaeological features of the Roman roads over Region B and C were then digitised, as indicated in Figure 4c and Figure 5c, respectively.To quantify the landscape change from 1984 onwards, inter-epoch elevation differences were generated by subtracting each DSM from the immediately more recent DSM on a pixel-by-pixel basis (Figure 4c and Figure 5c).The 1984 and 2006 digitised features are overlaid upon the 2016 elevation differences to indicate their possible location.

Figure 4 .
Figure 4. Results of the applied image enhancement pipeline for the 1984, 2006 and 2016 datasets alongside inter-epoch elevation differences and digitised archaeological features over Region B (© Ordnance Survey (OS) mapping copyright 2016).

Figure 5 .Figure 5
Figure 5. Results of the applied image enhancement pipeline for the 1984, 2006 and 2016 datasets alongside inter-epoch elevation differences and digitised archaeological features over Region C.

Figure 6 .
Figure 6.3D Scatterplots of pixel brightness value per band for the 2006 orthomosaic over Region C (a) before and (b) after decorrelation stretching.

Figure 7 .
Figure 7.The (a) 1984 and (b) 2006 Red-Green difference vegetation index with the 2006 surface morphological attribute of shaded relief over Region C.
. Aerial photographic datasets with details after SfM-MVS implementation.

Table 2 .
The M3C2 distances of the 2016-1984 epoch pair, both before and after the ICP implementation, are illustrated in Figure3.

Table 2 .
Statistics of point-to-point cloud differences.