An Objective Assessment of Hyperspectral Indicators for the Detection of Buried Archaeological Relics

Hyperspectral images can highlight crop marks in vegetated areas, which may indicate the presence of underground buried structures, by exploiting the spectral information conveyed in reflected solar radiation. In recent years, different vegetation indices and several other image features have been used, with varying success, to improve the interpretation of remotely sensed images for archaeological research. However, it is difficult to assess the derived maps quantitatively and select the most meaningful one for a given task, in particular for a non-specialist in image processing. This paper estimates for the first time objectively the suitability of maps derived from spectral features for the detection of buried archaeological structures in vegetated areas based on information theory. This is achieved by computing the statistical dependence between the extracted features and a digital map indicating the presence of buried structures using information theoretical notions. Based on the obtained scores on known targets, the features can be ranked and the most suitable can be chosen to aid in the discovery of previously undetected crop marks in the area under similar conditions. Three case studies are reported: the Roman buried remains of Carnuntum (Austria), the underground structures of Selinunte in the South of Italy, and the buried street relics of Pherai (Velestino) in central Greece.


Introduction
Remote sensing allows performing non-invasive large-scale analysis and monitoring of archaeological sites.The presence of buried archaeological structures can be detected by observing below the surface usually using ground-penetrating radar sensors [1,2].Laser scanning has been also employed in different applications [3] such as mapping areas beneath the vegetation [4], as in the recent case of the discovery of previously undocumented settlement features in Cambodia [5].
For the case of optical sensors, buried structures are usually detected in vegetated areas (cultivated or not) by studying vegetation anomalies on the surface known as crop marks, which are caused by differential growth in the vegetation.In presence of buried structures, the amount of fertile soil decreases causing stress in the vegetation (negative crop marks), while it can increase in presence of a ditch or higher concentration of organic material (positive crop marks).Similarly, differences in humidity or soil filling below the surface may generate anomalies in the soil known as soil marks, which can also be remotely observed.
Hyperspectral sensors measure the reflected solar radiation in hundreds of contiguous and narrow spectral bands, and naturally excel at estimating different parameters in vegetated areas, from health status to average leaf area and chlorophyll content in the plants.These tasks are usually achieved by combining the information of specific spectral bands in what is usually known as a vegetation index or, more generally, a spectral index.Specific vegetation indices can provide important information on plant stress in the areas concealing buried archaeological relics, as shown through the analysis of spectra acquired on ground in [6].Features obtained through other image analysis techniques, such as rotation of the feature space through Principal Components Analysis (PCA), may also provide additional insights regarding the presence of objects of interest [7].PCA may also be employed to reduce the dimensionality of the images, facilitating the data exploration step.
However, a detailed assessment of the qualitative and quantitative performance of spectral indices is still missing in the literature, as the quality comparisons carried out visually are usually estimated subjectively [8].An assessment restricted to PCA and Tasseled Cap features concludes that the former is better suited for such tasks [9].Nevertheless, PCA only represents a potential intermediate step for the analysis of these data, as it creates a new dataset with the same dimensionality of the original one, in which single principal components (PCs) with associated high eigenvalues must be anyway individually assessed to find the one better suited for cropmark detection.As the eigenvalues only reflect how much variance is captured by every PC, statistical variance in the data does not ensure that the same variations are present in the features of interest.In [6], the authors compared more than 70 vegetation indices using both broad and narrow bands, calculated for different phenological stages of barley and cereal crops using ground spectroradiometers: the indices found as relevant were then evaluated on an image acquired by the EO-1 Hyperion spaceborne sensor on the Thessalian plain in Greece.The figure of merit used in the work to evaluate the effectiveness of the extracted features is based on the measured contrast in the area of interest, which may not be due to the actual presence of buried structures, as for the mentioned variance in PCs.Similar approaches computed spectral separability as a criterion to evaluate the efficiency of the different indices, without explicitly quantifying their difference in terms of performance [7,10,11].
The lack of an unbiased and objective assessment for the hundreds of potential spectral features which could be extracted for a given dataset makes it difficult for a non-specialist to select only a few among them as a pre-processing step to be followed up by visual analysis.On the one hand, an in-depth study of only a few of them could leave crucial information undetected; on the other hand, considering all of the possibilities would be time-consuming and very hard to put in practice for an operator.
Taking a first step in tackling this issue, this paper proposes to use information theoretical measures to quantitatively estimate the mutual dependence of a given spectral index and the location of buried archaeological structures.If a map highlighting the presence of buried relics is available for a restricted area of interest with low spectral contrast, the mutual information computed between this and a given spectral index could automatically rank the various available spectral features according to their discrimination power for that region of interest.
Information theory has provided along the years several valuable tools for remote sensing applications, especially for model selection, texture extraction, distortion evaluation and informational content characterization of satellite images [12].For specific applications, ideas derived from Shannon's theory have been applied to Synthetic Aperture Radar (SAR) data [13,14], and to derive similarity measures for spectral signatures acquired by hyperspectral sensors [15].The notion of mutual information (MI) is of particular interest to estimate the statistical dependence of images, regardless of the range of the image values or sensor's characteristics: as an example, in [16] the authors carried out MI-based coregistration of optical and SAR images.
This paper exploits the statistical dependence between maps of known underground archaeological features derived from previous studies and relevant spectral features extracted from a hyperspectral dataset.The goal is evaluating the capability of the retrieved images to detect the buried archaeological remains.The method provides a list of retrieved maps which can be easily computed from a hyperspectral dataset, with high discrimination power for a specific archaeological task for which some targets are known.
The paper is structured as follows.Section 2 presents the method for ranking of spectral features according to their relevance at detecting crop marks based on mutual information, after a brief reminder on spectral indices.Section 3 presents two airborne hyperspectral datasets and one satellite multispectral dataset.The derived buried structures maps are ranked and discussed in Sections 4 and 5.We conclude in Section 6.

Spectral Indices
Spectral indices are derived from optical imaging data, by highlighting the empirical correspondences between an observed physical parameter on ground and a ratio employing two or more specific spectral bands in an image acquired over the same site.In hyperspectral research and applications, spectral indices are usually derived from images in reflectance using narrow spectral bands.Spectral indices may be derived also from multispectral images, if the spectral features of interest are broad enough.This is the case for one of the most popular spectral indices, the Normalized Difference Vegetation Index (NDVI) [17], which estimates the fraction of an image element covered by alive and green vegetation.Vegetation indices aim to take advantage of the unique features of healthy vegetation's spectral profiles, especially in the range exhibiting an abrupt transition between the strong absorption in the visible red frequencies and the strong reflectance in the near infrared.This spectral range, usually known as Red Edge, is linked with the spongy mesophyll cells located in the interior or back of the leaves which reflect near infrared radiation [18].
Regarding archaeological applications, the dominant source of data used for prospection is represented by aerial pictures acquired with a standard camera, and therefore lacking the infrared information.Nevertheless, this portion of the spectrum conveys very important information for such applications [18].Spectral indices linked to the health status of the vegetation may be used to reveal crop marks as local differences in the health of vegetated areas, which in turn may indicate the presence of buried archaeological relics.This detection is also possible for moderately stressed vegetation, well before the effects become visible in RGB pictures by altering the spectral response of the plants in the red region of the spectrum.In this study both broad-band and narrow-band vegetation indices have been ranked according to their correspondence with the known ground truth in archaeological areas.

Mutual Information
The assessment of the spectral information considered in this paper is based on the estimation of the amount of information shared by two objects, more specifically two random variables X and Y.Each random variable represents an image, so X may stand for a map of the known buried relics in a sensitive area, while Y is one of the several spectral indices computed on the same area, as an image can be considered the output of a stochastic process.The two random variables X and Y must have the same size, meaning in our case that they must consist of the same number of pixels.
From Shannon's probabilistic point of view, the quantitative estimation of the information shared between X and Y is done via the mutual information I(X, Y) [19].This quantity measures the amount of information that can be obtained about a random variable by observing another one, and is defined as: Here p(x) and p(y) are the probability distributions of X and Y, while p(x, y) is the joint probability of X and Y.The mutual information is positively defined, with a value I(X, Y) = 0 indicating statistical independence of the two variables.The units of information depend on the base of the logarithm, which is usually 2, 10 or e (in this work base 2 is used), and the maximum value I(X, Y) can assume is related to the entropy values of X and Y.
Mutual information is able to identify non-monotonic and generally non-linear correspondences between discrete variables, therefore it can be successfully applied to estimate the mutual dependence between two variables expressed in different units of measure.For example, X may contain radiance values as recorded from a sensor and Y a labeling of the same size, created manually.
To illustrate intuitively how the mutual information works, we report an example in which the objective is identifying the more suitable spectral band in a hyperspectral image to separate different kinds of crops.We analyze a NASA Airborne Visible-Infrared Imaging Spectrometer (AVIRIS) hyperspectral scene containing agricultural fields acquired over the Salinas Valley in the USA.The image has a size of 512 × 217 pixels with 192 spectral bands in the range from 0.288 to 2.5 µm, and with a Ground Sampling Distance (GSD) of 3.7 m.The ground truth identifying the different crops in the image is reported in Figure 1c, and is labeled from 0 (unclassified, reported in black) to 16, with a number and a color assigned to each different crop.As expected, band 1 centered at 380 nm is severely affected by noise in Figure 1a due among other effects to strong atmospheric absorption, Rayleigh scattering, and lower overall power in sun's emissions in the spectral range between ultraviolet and blue light.Band 42 centered at 750 nm in Figure 1b has opposite characteristics, exhibiting therefore a higher SNR.On a first visual inspection, it looks that the band in Figure 1b would be more effective at discriminating the different crops in the area, and this can be verified by analyzing the joint probability distribution for the two spectral bands in the lower part of Figure 1, as follows.
The scatter plots on the bottom of Figure 1a,b represent gray values in the specific band on the x-axis and the crop they belong to on the y-axis.The higher the amount of pixels with a given gray value x i belonging to a given class y i , the higher the joint probability p(x i , y i ).In the color-coded diagrams, a low value at coordinates (x j , y k ) means that only a few pixels having a radiance value x j belong to class y k , and is shown in dark blue, while the highest values are represented in red.If for a given spectral band or index different ranges of radiance values can be associated to different classes, we expect that band or index to be effective at discriminating the ground covers.On the other hand, if a narrow range of radiance values can belong to several classes such as in Figure 1a, it would be difficult to classify a dataset based on the same spectral information.In layman's terms, if different values in one image correspond to different values in the other one for a high number of pixels, the mutual information goes up, regardless of the domain, scale and ordering of these values.This represents an advantage over other quantities such as correlation, which would fail in this case because it requires both variables to increase or decrease monotonically in a comparable way.Therefore, in the joint probability distributions reported, it is desired that all pixels belonging to a given interval of radiance values are associated to a minimum number of classes on ground (ideally one), in order to minimize confusion.For example, a value larger than 7500 corresponds in Figure 1b to the class labeled as 1, and values around 2000 to the class labeled as 5.Even though some confusion between classes remains, such as the overlap between the distributions of the values for classes 9, 10 and 15, the case of Figure 1a is much more confused.For example, a value slightly smaller than 300 in that band could correspond basically to any class, and observing such radiance value would not give any meaningful information on differentiating the crops in the scene.
The mutual information quantifies the discussed spread in the joint probability distribution: referring to Equation (1), the total value increases meaningfully whenever a pair p(x i , y j ) is considerably larger than 0 and similar in value to ∑ n y=1 p(x i ) or ∑ m x=1 p(y j ) (ideally both).In our cases, n and m represent the total number of classes and gray values assumed by the dataset, respectively.An increase in mutual information would therefore decrease the confusion in recovering a label in Y by observing the radiance value of the same pixel in X.In this paper, we proceed to assess the value of the information extracted from a hyperspectral dataset as follows.First, prior knowledge about the area of interest is exploited in order to derive a map with an interpretation of underground anthropogenic features, which is binary and set to 1 for image elements considered to contain known concealed archaeological relics and to 0 otherwise.Subsequently, a series of spectral indices is extracted from the available hyperspectral image and stored as a set of two-dimensional images.
The mutual information is then computed between the described map and each spectral index: a high value indicates that both vary in similar spatial locations, implying that a change in the map is usually reflected in a change in the value of the computed spectral feature.It must be remarked that such changes can be positive or negative: the spectral index can either increase or decrease for a pixel which really contains a buried structure, but it must vary in the same locations where the relics are located in order to yield a high mutual information value, which is in any case always a positive quantity.
This procedure has two advantages over assessments carried out by measuring the contrast in a sensitive area: firstly, a variation in the data contributes in increasing the mutual information only if it overlaps to real changes as collected in the reference map.Secondly, the mutual information computed for spectral indices with different value ranges can be directly compared, as the amount of variation in absolute value does not affect the mutual dependence of the two variables.
If the mutual information between an image representing a given spectral index and the reference map is 0, this means that buried structures are not visible in the considered image.On the other hand, high values imply that buried structures should be clearly outlined.

Carnuntum
Carnuntum, located on the southern bank of the Danube River, was the capital of the former Roman province known as Pannonia superior.Previous non-invasive studies on this archaeological site revealed that its area spanned more than 6 square kilometers, and it included a school of gladiators [20].An airborne hyperspectral dataset acquired above this archaeological landscape by the sensor AisaEAGLE (Specim Spectral Imaging Oy Ltd., Oulu, Finland) is available for download as a supporting dataset for the ARCTIS MATLAB R toolbox, which is released under the creative commons license [21].The image comprises 65 spectral bands in the range 400-1000 nm, and has a GSD of 0.4 m [8].
The map reported in Figure 2 has been used as a ground check for the buried structures.The map has been manually derived using as reference the results contained in [8], validated by a complete interpretation of aerial archaeological evidence collected in a time span of 50 years [22].In order to focus on a meaningful area and decrease variations in the results caused by other kinds of spectral variations, we restrict our analysis to the yellow rectangle overlaid on the image in Figure 2.

Selinunte
Selinunte is one of the most outstanding cultural heritage sites of Southern Italy, located along the SW coast of Sicily, founded in the 7th century BC by colonizers who came from Megara Hyblaea (ancient Greek colony in Sicily).Selinunte was built over three N-S hills orthogonal to the coastline: moving eastwards, they are the Gaggera, Manuzza and Eastern hills (Figure 3).The public and religious centre of the ancient city (i.e., Acropolis) was located over southern part of Manuzza Hill.The Gaggera and Eastern hills, i.e., two religious areas, were connected to the Acropolis area by a main road, i.e., Plateia, crossing the city area and passing through the urban defence walls.During the two centuries that followed, the city spread eastwards and westwards, along Cottone and Modione stream valleys, which stretch along the borders of Manuzza Hill.The street pattern of Manuzza Hill was first identified on aerophotos by Schmiedt [23] and partly confirmed by successive archaeological investigations [24].[7].MIVIS (Daedalus AA5000) records the incoming radiation into four optical ports covering the Visible (VIS), NIR and SWIR regions from 0.43 to 2.45 µm, plus Thermal Infrared channels from 8 to 13 µm, for a total of 102 channels.The data were calibrated and corrected as described in [7].

Ferai (Velestino)
The ancient settlement of Pherai expands mainly below the modern village of Velestino in Thessaly, central Greece.The site exhibits a continuous habitation from the Neolithic period (magoula Bakalis) until the Early Roman Imperial period (1st ce.CE) [25][26][27], when it was abandoned to be settled again after the Byzantine period (13th ce.AC).The settlement flourished during the Archaic and Classical periods.Fortification walls were also built for the protection of the acropolis and other sectors of the city.Even though the exact boundaries of the ancient settlement remain unknown, the total extent of it has been estimated to range between 80-120 ha [28].
Various geophysical prospection surveys have been conducted in particular by the GeoSat ReSeArch Lab of FORTH, focusing on the Neolithic settlements of Bakalis and Velestino, the Hypereia spring and the Hellenistic stoa.In 2014, a more extensive survey aiming towards the study of the urban development of the site was carried out at the northern sections of the village, where the northern limits of the ancient Greek settlement were expected [29].A high resolution spaceborne satellite image acquired on the 15th of June 2009 on the site by the Quickbird sensor is avalaible.The image comprises four spectral bands (RGB and NIR) with a GSD of 2.4 m.A true color subset indicating the locations of two areas used in next section is reported in Figure 4.

Results
Several spectral indices have been derived from the described datasets and ranked according to their mutual information with a map of known buried structures.A complete list of the indices and further information on them can be found in Appendix A. The reported indicators employ a wide range of spectral bands: although their vast majority is applicable for generic vegetation health monitoring, several of them are tailored for specific vegetation parameters.In some cases, improved versions of already existing indices are tested, such as the Atmospherically Resistant Vegetation Index (ARVI) [30] which can be considered as an improved Normalized Difference Vegetation Index (NDVI) [17] robust to atmospheric effects, or the improved version of the Modified Chlorophyll Absorption Ratio Index (CARI) [31] among several variants of the CARI.Subsequently, the mutual information between different spectral features and the binary maps of buried structures in the areas has been computed for several spectral features.We expect a higher mutual information between the two to denote a better visibility of the archaeological features of interest.In comparison to simple contrast analysis, mutual information implicitly evaluates if an increase of contrast in a given spectral feature is reflected by a change in the reference map.

Carnuntum
The mutual information values for single spectral bands and all the Principal Components reported in Figure 5 were retrieved for the Carnuntum dataset.A visual inspection of these spectral features in Figure 6 confirms that a higher mutual information indicates that buried structures are more clearly visible, as follows.In Figure 6(1) a true color combination of the hyperspectral image is reported, in which crop marks are hard to observe.The average mutual information value for the three bands used is 0.017, while it increases up to 0.114 for band 43, centered at 787 nm in the infrared range, where the presence of features in the field corresponding to the buried structures beneath is obvious (Figure 6(2)).The four Principal Components related to the eigenvectors with the highest associated eigenvalues are reported in Figure 6(3-6).It is remarkable that the computed mutual information for all PCs is largest for the first and the fourth one, which happen to be the same ones in which crop marks are visible in Figure 6, with values of 0.124 and 0.055, respectively.
As expected the average values related to the single spectral bands are higher with respect to the ones from the PCs, as the former are strongly correlated and the information they contain is highly redundant.On the other hand, the first PCs capture the relevant information in the image, as they are ideally decorrelated and ordered according to their variance: therefore, low-ranked PCs are mostly composed by noise and have negligible utility for further analysis.The ranking for the spectral indices collected in Appendix A is reported in Figure 7.In the figure, a performance expressed in percentage is derived for each index by comparing the related mutual information value with the one obtained by the best index.In order to display clearly the differences between the results, only a rotated relevant detail of the image is shown in the subfigures.For each subfigure a meaningful histogram stretch has been derived by applying a linear-with-2%-saturation stretch restricted to the area of interest.A visual analysis makes obvious that the top-ranked indicators clearly enhance the differences between crop marks and the surrounding areas: as the mutual information drops, concealed archaeological structures get harder to spot, up to becoming undetectable for an observer for the indices with associated lowest ranks.
The values reported in Figure 7 are relative quantities: to better understand them, we must consider that by definition I(X, Y) ≤ min (H(X), H(Y)), where H(X) and H(Y) is the Shannon entropy of X and Y, respectively, defined as Therefore, the maximum possible value which I(X, Y) can assume is related to the entropy for the binary map of the buried structures interpretation as in Figure 2, which being equal to 0.7 is much smaller than the entropy of any spectral band or index which can assume a wide range of values.The normalized value I norm = I(X, Y)/ min (H(X), H(Y)), comprised between 0 and 1, represents the decrease in uncertainty in the analyzed spectral index with respect to the maximum possible value.In our experiments, I norm assumes a maximum value of 0.213 and a minimum value of 0.043, for BAI and Photochemical Reflectance Index (PRI), respectively.These correspond to mutual information values of 0.15 and 0.03, respectively.All other values are accordingly proportional to the ones reported in Figure 7.The indices are sorted from best to worst according to the mutual information with the map reported in Figure 2. A score is reported as the relative mutual information with respect to the one computed for the best index.The buried structures become less evident as the mutual information decreases.
The index with the associated highest score is the Burned Area Index (BAI) [32].Although BAI is primarily used to assess damages caused by fire, previous studies have reported this index to exhibit high discrimination capabilities for the task of detecting crop marks in the Mediterranean region [6].The performance of BAI shows that indices originally defined for different applications can be employed successfully in archaeological research, and highlights that simple (i.e., linear) modifications of existing indices can even improve the final results and assist the interpretation of the image.Other interesting results include the Simple Ratio (SR) index [33] which, in spite of its simplicity, was ranked third.The narrow spectral ranges centered around 740 nm and 720 nm used to calculate the Vogelmann Red Edge (VRE) index [34] seem to be among the most promising spectral windows for detecting crop marks in areas with similar conditions, as found also in other studies [18,21].The NDVI index, widely used in several remote sensing archaeological applications, was able to enhance the crop marks but less effectively than other -seldom used-indices: nevertheless, it must be remarked that it has been applied using two specific narrow bands (for further details see Appendix A); results obtained using different bands or simulating typical wider bands which characterize multispectral sensors could yield very different results.Indeed, different sensors are characterized by different spectral responses and are affected by specific uncertainties in their measurements.An essential part of the latter is related to different behaviours of the Relative Spectral Response filters which describe the relative sensitivity of the sensor to energy at different wavelengths [35].Attempts to quantify these variations have been made by analyzing Full Width at Half Maximum characteristics [36].Even though the total relative uncertainties due to calibration for vegetation indices should be within 5% for spaceborne sensors [37], recent studies performed by the authors have shown that in some cases (e.g., IKONOS) this could reach differences up to 8%, well beyond the pre-calibration errors of the sensor [35].Therefore, specific characteristics of each system should be taken into consideration, prior to their application for archaeological research.Global Environmental Monitoring Index (GEMI) [38] and ARVI were also ranked among the most promising indices, just below SR.The close ranking of the ARVI and NDVI indices is justified by the accurate atmospheric correction carried out previously on the image [8].If this pre-processing step is not carried out, the two indices are expected to perform similarly only if atmospheric scattering and absorption effects are limited during the airborne flight.Although this can be true is some cases, ARVI is usually the best choice whenever scattering and absorption phenomena are strong or unknown, which is often the case whenever researchers have to deal with archive datasets [39].
Generally, it is evident that the ranking of the indices is case-sensitive: nevertheless, some interesting general conclusions can be drawn from the described experiment.Firstly, the ranking shows that indices originally defined within different fields of study such as BAI can be applied to archaeological research; secondly, the most sensitive spectral range to the variations between healthy and stressed vegetation, where the latter leads to the formation of crop marks, ranges from the red to the NIR portions of the spectrum.It is interesting to remark that any single band in this spectral region would underperform if compared to these indices (see Figure 5), as it would be ranked around position 22 among the results considered in Figure 7. Finally, improvements of existing indices may enhance the final outcomes: see the performance for different versions of the Chlorophyll Absorption Ratio Index (CARI).

Selinunte
The proposed method was applied to spectral indices obtained from the Selinunte dataset.The area of Manuzza Hill (i.e., yellow rectangle in Figure 8) was selected because the street network, which was identified in [23], was partly detected from the analysis of MIVIS data [7,40,41] (Figure 8) and was characterized by homogenous land cover.The street network of the Manuzza hill includes a primary road (i.e., Plateia) and a series of parallel secondary street axes (i.e., Stenophoi).The Plateia are as wide as about 6m (i.e., about two pixels size), whereas the Stenophoi are, at maximum, 3 m wide.In Figure 8, pixels related to Stenophoi street axes may be mixed (i.e., only partially belonging to the class of interest) and thus unreliable: therefore, the mutual information is computed between the main crossroad and the same set of indices derived in Figure 7.The ranking of the indices from best to worst once again well agrees with the visual interpretation of the data, as in general the visibility of the main crossroad decreases along with the relative Mutual Information score (Figure 9).The absolute values range from 0.1 to 0.03, except for the SIPI index, which was corrupted and appears almost completely saturated, having a value of 0 correctly representing absolute statistical independence from the features of interest.The most interesting aspects are the ones confirming the applicability of the described method for a specific case of study only.Firstly, it is remarkable that the parallel vertical roads reported along the main crossroad in Figure 8 are better visible roughly in the indices ranked among the top 10: this suggests that unknown buried structures may be detected once the best spectral features are identified.Secondly, indices such as BAI and VRE, which were very good indicators in Carnuntum, perform in this case poorly, and exhibit barely visible features, confirming that the performance of a specific spectral index is case-sensitive and also depends on the spectral configuration of the imaging sensor.
We verified the ranking of the indices as computed from Manuzza hill on the west Acropolis area (i.e.blue rectangle in Figure 8), which features a series of parallel secondary street axes (i.e., Stenophoi).The indices are reported in Figure 10 with the same contrast stretch used in Figure 9.The archaeological features are known also in this area, as reported in Figure 8, but the reference map was not used in this case to rank the spectral indices.The ranking well agrees with the visibility of the buried relics in the north of this image subset, as both Manuzza Hill and the west Acropolis areas were characterized by the same land cover (i.e., grass and xeritic shrubs) and same lithological outcrops (i.e., weathered calcarenite soil).On the one hand, the vegetation cover in the southern area is more sparse with respect to the area containing Manuzza Hill; on the other hand, the west Acropolis area is located in the valley of the Modione stream and its vegetation is less dry [7].The indices are sorted from best to worst, according to the mutual information with the main crossroad in the area reported in the yellow rectangle in Figure 8.A score is derived as the relative mutual information with respect to the one computed for the best index.As in the case of study reported in Figure 7, the buried structures become less evident as the mutual information decreases.8).The indices are sorted according to the mutual information with the main crossroad in the area reported in the yellow rectangle of the same figure.

Ferai (Velestino)
We analysed the available Quickbird image of the site, employing the original multispectral dataset instead of the pansharpened one, as in the pansharpening process the archaeological features had been visibly corrupted, becoming much less evident in spite of the increase in resolution.Therefore, the following experiments suggest that better results may be obtained by using original multispectral acquisitions at lower resolutions for archaeological applications.
Two different sections of the site have been selected for the particular study: the first area is located where a preceding geophysical study identified the road network at the northern limits of the ancient Greek settlement (green rectangle in Figure 4).The second area was located to the south-east, where visual inspection of the panchromatic image has given some suggestion for the existence of additional ancient streets (red rectangle in Figure 4).
The studies carried in situ combined the employment of magnetic (multi-sensor SENSYS magnetometer MX system) and electromagnetic techniques (GEM-2 multifrequency EM and CMD multi-oils spacing EM), which were effective at revealing part of the ancient street system on the site (Figure 11).About 6-7 NW-SE streets were identified, having a width of about 6-7 m (Figure 12).Sections of the particular features are also visible in the high-resolution Quickbird and GeoEye-1 satellite image (the latter is not reported in this study), and in historical aerial photographs acquired from the 1960s and 1970s [29].As the image is multispectral, we derived a set of broadband indices which differ from the ones analysed in the previous cases.In common with the previous cases of study, we derived GARI, GNDVI, NDVI, SAVI, SR, Iron Oxide, MCARI, NLI, and MTVI.We also computed other indices such as Red Green Ratio, Modified NLI and Optimized SAVI [6].
The MI between the indices and the vector map for the buried features existing in the field reported in Figure 13 have been computed and the resulting top 6 and bottom 6 ranked spectral indices are reported in Figure 14.It is evident how the features of interest become hard to spot in the low-ranked indices.In all of these images a linear stretch with 2 % saturation has been applied, and the same stretch has been used to analyse the rest of the image to look for previously unknown archaeological structures.An interesting subset of the image with similar soil conditions is reported in Figure 15.In this figure, the green lines represent the extension of the street network of the ancient settlement to the south.Even though there have been no excavations in the particular area, the visual inspection of the true color subset of the Quickbird multispectral image indicates clearly the residues of the streets that have a parallel alignment and they are consistent with the spacing of other roads to the north.White and red circles represent respectively known and previously unknown archaeological features which are observable in the top ranked indices under the same conditions and histogram stretch.The subsets of the spectral indices of interest are reported in Figure 16.
In general we notice that the top 6 indices (GARI, GNDVI, SR, NDVI, OSAVI and SAVI) are able to highlight the street network in both sections of the site, as they have similar type of crops that were at the same phenological stage.The results of the study are of importance not only from the image processing aspect, but also from the archaeological perspective since they provide an additional evidence of the extent of the city which is in accordance to the hypotheses made in the past [28].On the other hand, the above indices are different from top-ranking ones in the previously reported case studies (Carnuntum and Selinunte), indicating that there is a variation of the discrimination power of each index according to the type of crops and their phenological status, among other parameters.4).Among the broadband indices computed, the top 6 and bottom 6 are reported according to the mutual information with the reference map of archaeological relics.The buried structures become hard to spot in the indices with an associated low score.

Discussion
The detection of crop marks depends on several factors which are usually only partially known, such as climate, characteristics and phenological stage of the crops [42], agricultural practices, material and depth of relevant buried structures, and soil characteristics.Therefore, some studies have concluded that archaeological sites do not exhibit specific spectral characteristics that can be exploited to efficiently detect crop marks in any remotely sensed image, as each site may exhibit different geomorphological characteristics [7,40,41].Given that no universal algorithm can be efficiently applied to all cases studies, several spectral indices might perform similarly at revealing the presence of archaeological relics on a given area of interest [6].
This study aims at filling the existing gap in the literature between the extraction of spectral information and the preliminary evaluation of their efficiency for a given task, providing archaeologist with a case-sensitive workflow which can be tailored for each archaeological site of interest to implicitly account for variations due to the phenological status of the vegetation and climate in the area of interest.Another advantage of the proposed technique is that mutual information is robust to changes in the representation of the data from DN to radiance to reflectance, as the relations on which it is computed are non-linear.The application of the proposed methodology confirms that the capability of the retrieved images to detect the buried archaeological structures is case sensitive, implying that results obtained on one site cannot be directly applied to others, as the targets may exhibit different spectral features, due to the mentioned different characteristics.
The reported experiments show that the obtained ranking of features well agrees with a visual interpretation of the area for a series of underground buried structures, and may increase the chances of revealing new details on known buried structures, or to discover new underground relics in the same area.
The main drawback of the method is that an accurate (even if incomplete) map of the underground structures in an area of interest must be given, and a large enough number of image elements must be related to the features of interest, in order for the computed statistics to be meaningful.Another aspect of interest is that results are limited to the type of feature contained in the reference map.For example, if only buried walls are contained, the indices may be ranked according to their discrimination power to find walls, but not necessarily other structures, as for example pits which are usually related to more vigorous vegetation, or relics composed of other materials which altered the soil in different ways.If a detailed map containing different classes exists, though, mutual information would excel in ranking the features based on their overall discrimination power: going beyond the binary map would increase the entropy of the reference map, and along with it the range of values for the mutual information and its discrimination power [19].
As mentioned, no generic definition of a spectral signature having intrinsic archaeological interest exists: in [43] spectral bands particularly sensitive to vegetation stress have been identified and their use is suggested to maximize the probability of detecting crop marks if no a priori knowledge of the site is available.On the other hand, several additional agents could play a role in the formation of crop marks, and therefore for the efficiency of a given spectral indicator in a specific archaeological site or typology of sites.The proposed case-sensitive approach is ideally independent from the set of variables which usually need to be taken into account for this task, being based solely on the statistical dependence between a spectral index and the buried archaeological features.Such approach could be used on a pre-flight campaign, in which a limited amount of hyperspectral data on known buried sites is acquired in order to estimate the best indicators.The integration of geophysical prospection can improve final results by exploiting additional hidden targets for the evaluation.
Current space-borne hyperspectral sensors provide medium resolution datasets, which are sometimes difficult to efficiently employ in archaeological research.However, future hyperspectral missions characterized by a short revisit time will enable for the first time the acquisition of reliable hyperspectral image time series: future work will therefore include the study of phenological changes in the crops and how these can be captured via mutual information.

Conclusions
Remote sensing provides an efficient non-invasive technology for the exploration of archaeological sites.This paper focuses on the evaluation of information extracted from hyperspectral dataset for the detection of concealed archaeological structures, usually performed indirectly by observing stressed or very vigorous vegetation through the appearance of the so-called crop marks.
So far, few efforts have been made to assess objectively and quantitatively the discrimination power of the hundreds of spectral indicators linked to the health status of vegetation which could be potentially extracted from a hyperspectral image for this purpose.In order to overcome this limitation, this paper proposes to use information theoretical measures to quantitatively estimate the mutual dependence of a given spectral index and the location of buried archaeological structures for a specific case of study.If a map containing the location of a few buried relics is known prior to an ongoing excavation for a restricted area of interest with low spectral contrast, the mutual information computed between such map and a given spectral index could automatically rank the available spectral features according to their discrimination power.This semi-automatic ranking of the indices-for each specific case study-according to their 'archaeological' informational content may help experts in focusing on other areas in close proximities, based on the most meaningful indices for further investigation activities.Such researches are usually carried out via visual analysis, and the pre-selection of the most suitable spectral features could maximize the success rate of the interpretation of crop marks.

Appendix A
This Appendix contains relevant information on the list of spectral indices considered in the experimental section.The fields Rank C and Rank S represent the performance of each index for the experiment reported in Figures 7 and 9, respectively.In the reported equations the following notation is used: ρ x = reflectance value at x nanometers.2.
Blue = reflectance value at 470 nm.

Figure 1 .
Figure 1.Graphical analysis of the joint probability for different spectral bands of the Salinas dataset.(a) Band 1 of the dataset (top) and joint probability distribution between radiance values in band 1 and the ground truth reported in (c) (bottom).(b) Band 42 of the dataset (top) and joint probability distribution between radiance values in band 42 and the ground truth reported in (c) (bottom).(c) Available ground truth representing the distribution of crop classes (top, black means no information) and sample spectral library containing one full spectrum per class (bottom).

Figure 2 .
Figure 2. Binary reference map derived from known maps of buried structures in the Carnuntum archaeological site (in red), overlaid on spectral band 43 of the hyperspectral dataset, centered around 787 nm (Data source: LBIArchPro and ABT GmbH).The area in the yellow rectangle has been considered to quantify the discrimination power of the extracted spectral features.

Figure 3 .
Figure 3.The study area of Selinunte.

Figure 4 .
Figure 4. True color subset of Quickbird multispectral image acquired over the study area of Ferai.The green rectangle shows the area used to rank the spectral indices, which were then examined in the area denoted by a red rectangle with the same histogram stretch to look for buried archaeological features.

Figure 5 .
Figure 5. Mutual information between known buried sites and single spectral bands and principal components.A higher mutual information value indicates a better discrimination of the archaeological features of interest.For single spectral bands, this increases in the NIR region with respect to the visible spectral range.The top x-axis reports the wavelength for each band and refers to the blue plot only.See some examples in Figure 6.

Figure 6 .
Figure 6.Mutual Information: spectral bands and first four PCs for the Carnuntum dataset (Data source: LBIArchPro and ABT GmbH).

Figure 7 .
Figure 7. Detail of spectral indices described in Appendix A computed on the Carnuntum dataset.The indices are sorted from best to worst according to the mutual information with the map reported in Figure2.A score is reported as the relative mutual information with respect to the one computed for the best index.The buried structures become less evident as the mutual information decreases.

Figure 8 .
Figure 8. Map of buried archaeological structures in the Selinunte area which were identified from MIVIS data [7,40,41], overlaid on spectral band 15 of multispectral infrared and visible imagingspectrometer (MIVIS) dataset, centred around 720 nm.The area in the yellow rectangle has been considered to quantify the discrimination power of the indices, while the area in the blue rectangle has been used to verify the ranking on a nearby site.

Figure 9 .
Figure 9. Detail of spectral indices described in Appendix A computed on the Selinunte dataset.The indices are sorted from best to worst, according to the mutual information with the main crossroad in the area reported in the yellow rectangle in Figure8.A score is derived as the relative mutual information with respect to the one computed for the best index.As in the case of study reported in Figure7, the buried structures become less evident as the mutual information decreases.

Figure 10 .
Figure 10.Detail of spectral indices described in Appendix A computed on the Selinunte dataset in the West Acropolis area (blue rectangle in Figure8).The indices are sorted according to the mutual information with the main crossroad in the area reported in the yellow rectangle of the same figure.

Figure 11 .
Figure 11.Results of the magnetic survey from the north section of the village of Velestino.

Figure 12 .
Figure 12.Diagrammatic interpretation of signals registered by the geophysical (magnetic and electromagnetic surveys) and the satellite imagery from the north section of the village of Velestino.

Figure 13 .
Figure 13.True color subset of Quickbird multispectral image (no pansharpening) in the northern section of Ferai.The overlaid vector shows evident buried structures which were used to rank the spectral indices.

Figure 14 .
Figure 14.Detail of part of the spectral indices described in Appendix A computed on the Ferai image subset in the green rectangle of Figure 4.Among the broadband indices computed, the top 6 and bottom 6 are reported according to the mutual information with the reference map of archaeological relics.The buried structures become hard to spot in the indices with an associated low score.

Figure 15 .
Figure 15.True color subset of Quickbird multispectral image of the Velestino dataset.The overlaid vectors in green show known archaeological features.The areas circled in white contain observable crop marks for the top ranked spectral indices reported in Figure 16.The area circled in red contains a previously unobserved feature which becomes visible in the top ranked indices in the same figure, where it resembles a prolongation of one of the known buried structures.

Figure 16 .
Figure 16.Detail of part of the spectral indices described in Appendix A computed on the Ferai subset in the red rectangle of Figure4).Among the broadband indices computed, the top 6 and bottom 6 are reported according to the mutual information with the reference map of archaeological relics.The buried structures become hard to spot in the indices with an associated low score.