Detection of Neolithic Settlements in Thessaly (Greece) Through Multispectral and Hyperspectral Satellite Imagery

Thessaly is a low relief region in Greece where hundreds of Neolithic settlements/tells called magoules were established from the Early Neolithic period until the Bronze Age (6,000 – 3,000 BC). Multi-sensor remote sensing was applied to the study area in order to evaluate its potential to detect Neolithic settlements. Hundreds of sites were geo-referenced through systematic GPS surveying throughout the region. Data from four primary sensors were used, namely Landsat ETM, ASTER, EO1 - HYPERION and IKONOS. A range of image processing techniques were originally applied to the hyperspectral imagery in order to detect the settlements and validate the results of GPS surveying. Although specific difficulties were encountered in the automatic classification of archaeological features composed by a similar parent material with the surrounding landscape, the results of the research suggested a different response of each sensor to the detection of the Neolithic settlements, according to their spectral and spatial resolution.


Introduction
The spectral capability of early satellite sensors opened new perspectives in the field of archaeological research. The recent availability of hyperspectral and multispectral satellite imageries has established a valid and low cost alternative to aerial imagery in the field of archaeological remote sensing. The high spatial resolution and spectral capability can make the VHR satellite images a valuable data source for archaeological investigation, ranging from synoptic views to small details [1].
Since the beginning of the 20th century, aerial photography has been used in archaeology primarily to view features on the earth's surface, which are difficult if not impossible to visualize from the ground level [2][3][4]. Archaeology is a recent application area of satellite remote sensing and features such as ancient settlements can be detected with remote sensing procedures, provided that the spatial resolution of the sensor is adequate enough to detect the features [5].
A number of different satellite sensors have been employed in a variety of archaeological applications to the mapping of subsurface remains and the management and protection of archaeological sites [6][7][8][9][10][11]. The advantage of satellite imagery over aerial photography is the greater spectral range, due to the capabilities of the various on-board sensors. Most satellite multi-spectral sensors have the ability to capture data within the visible and non-visible spectrum, encompassing a portion of the ultraviolet region, the visible, and the IR region, enabling a more comprehensive analysis [12]. Multispectral imagery such as Landsat or ASTER is considered to be a standard means for the classification of ground cover and soil types [13]. Concerning the detection of settlement mounds the above sensors have been proven to be helpful for the identification of un-vegetated and eroded sites [5]. In recent years the high spatial resolution imageries of IKONOS and Quickbird have been used for the detection of settlements and shallow depth monuments [14][15][16]. Hyperspectral imagery (both airborne and satellite) has been also applied in archaeological investigations on an experimental basis and need further investigation [2,17,18].
The goal of this particular project was the application of different methods and means of satellite remote sensing for the detection of Neolithic settlements. In this study four satellite remote sensing images with different spatial resolutions (ASTER, Landsat, HYPERION, IKONOS) were examined in order to search their potential for automatic extraction of Neolithic settlements, by means of pixelbased and object -based methods. This paper seeks to address these issues through a multi -sensor case study in Thessaly, Greece, where different satellite image processing techniques contributed to the detection of the so called 'magoules' that are found in the Thessalian plains. The satellite data were statistically analyzed, together with other environmental parameters, to examine any kind of correlation between environmental, archaeological and satellite data. Moreover, different methods were compared and integrated methodologies for the detection of Neolithic settlements were extracted. The results of the study suggested that the complementary use of different imagery can provide more satisfactory results.

Study Area and Data
Thessaly is a relatively closed geographical unit, with definite mountainous borders (Mt. Antichasia and Olympus in the north, Mt. Ossa, Mavrovouni and Pelion in the east, Mt. Othris in the south, and Mt. Pindus in the west, reaching heights of 2,000 m) and two accesses to the sea, one through the Tempe gorge (NE) and another between the Othrys and Maurovouni mountains to the gulf of Volos.. All Thessalian basins show continuous habitation during all phases of the Neolithic period. As a matter of fact, Thessaly is famous for its long-lasting sites on its extensive fertile soils ( Figure 1). The Neolithic settlement mounds are typically low hills of 1-5 meters height and a mean diameter of 300 meters, mainly consisting of loam and mud based materials. There are hundreds of Neolithic settlements/tells called magoules all over Thessaly, with different kind of vegetation now above them. Due to the intensive cultivation of the land in the past, not all of them are visible. Past field archaeological surveys were able to identify a number of them based mainly on the surface concentration of sherds and lithic material [19][20][21][22]. However most of the magoules (137) are mainly at East Thessaly (Larisa Plain) and less (63) in west Thessaly (Karditsa Plain). These two plains consist of Quaternary alluvial deposits. The study involved satellite image detection of Neolithic Settlements in Thessaly by incorporating the following satellite and digital spatial data (Table 1): -4 ASTER images.
-1 Landsat ETM image. -1 HYPERION image: Only 137 of the 242 total HYPERION bands were used in the analysis, because many of the bands exhibited low signal to noise ratio or other problems.
-4 IKONOS images: For each image, the multispectral bands were fused with the high resolution panchromatic band in order to exploit the spectral information of the four multispectral bands (blue, green, red, near infrared) and the effective spatial resolution of the panchromatic band. -18 Air photos acquired from the Geographic Service of the Hellenic Army -GYS.
-The results of topographic mapping through systematic GPS surveying of more than 342 Neolithic settlements of Thessaly. -A DEM of 20 m pixel size of the study area. The DEM was constructed after digitizing in GIS environment 24 topographic maps scale 1:50.000 from the Geographic Service of the Hellenic Army. It has to be mentioned that ASTER DEM was also exploited in the particular study but it did not cover sufficiently the whole area of interest, and second, the specific images have different area coverage and only the ASTER mosaic was able to cover the whole region of Thessaly.

Research Methodology and Results
The image processing of satellite data was carried out in two steps starting with the basic preprocessing procedures followed by more sophisticated image processing steps.

Preprocessing of Satellite Images
The construction of image mosaics ( Figure 2) followed the masking of the sea, the clouds and the snow areas using Erdas Imagine 9.1 software package. The next step had to do with the transformation of the projection systems of all images to the Hellenic Geodetic System of Reference (EGSA87/HGSR87) so that they can all be fused to the same projection system. The final step of image preprocessing was the conversion of DN (Digital Number) values of images to reflectance. Different equations to convert the DN values to radiance were employed. The conversion of the DN values of ASTER images was achieved through the equation: [23]. For the IKONOS images the equation: Lrad = DN/Unit Conversion Coefficient was used [24]. The conversion of DN values of Landsat images to radiance was accomplished through the equation: Lrad = DN * Grescale + Brescale where Grescale and Brescale are band specific rescaling factors [25]. For the case of HYPERION images "signal to noise" ratio was used to select 137 bands from the total of 242. Then DN values were converted to radiance values according to the equations: VNIRL = DN /40, SWIRL = DN / 80 (USGS, 2007). The last conversion had to do with the conversion of the radiance of all images to reflectance through the general algorithm by [26] (1): where : Pp unitless planetary reflectance L λ spectral radiance at the sensor's aperture d 2 earth-sun distance in astronomical units ESUN λ mean solar exoatmospheric irradiances Θ s solar zenith angle in degrees.

Composition of RGB Composites
Several RGB composites were constructed in an effort to examine their efficiency in the detection of the Neolithic settlements. For the ASTER image with acquisition date 19-03-2003, where most of the magoules are registered, the RGB→1,2,3, RGB→3,2,5 and RGB→2,3,7 composites ( Figure 3) were the most successful for the visual detection of the Neolithic settlements (Out of 239 settlements, 39 of them were highly visible, 49 average visible and 151 poorly visible). Those composites appeared to have the highest Optimum Index Factor. High OIF values indicate bands that contain much "information" with little correlation. By using the OIF method, three band components of an RGB can be evaluated on their effectiveness for display [27]. OIF is defined by equation (2).
where s i is the standard deviation of band i and r(ij) is the correlation coefficient of band i and band j. Similarly, RGB composites of IKONOS images were able to detect 27 in a total of 48 settlements. It has to be noted that 19 of the detectable magoules, namely the highest of all corresponding to an average altitude of 4.6 m, were highly visible in all RGB composites. On the other hand, RGB composites of Landsat and HYPERION images were not very promising (for HYPERION composites only five settlements were detected in a total of 21). Finally, average altitude aerial images contributed to an excellent detection of all the five settlements that were inside the spatial limits of the airphoto mosaic. As a general conclusion however, the most crucial factors for the detection of magoules proved to be the acquisition date of the image due to the fact that the land around the majority of the settlements is cultivated (Figure 4). Visual interpretation is commonly used for visual extraction of obvious and large or medium scale archaeological structures like settlement mounds [28][29][30]. For IKONOS images it was possible to detect most of the settlements with just a simple visual interpretation of any kind of RGB composite due to the high spatial resolutions of the specific image. The visual detection of them was achieved based on shape, linearity, tone, and texture size between different patterns around them [14]. The same task was accomplished for the airphotos. However, the lack of airphoto data and their small spatial coverage of the study area turned air photos to have ancillary role in the whole study ( Figure 5).

Spectral Profile Comparison and Classification
The identification of spectral signatures was considered to be a crucial task for the detection of Neolithic settlements especially for the classification process. That task was accomplished in order to exploit any potential distinct spectral characteristics of surface and subsurface settlements patterns compared with the surrounding material [2]. Signatures were collected from all tells and were divided into two categories: those collected from plain areas and those collected from mountainous areas due to different soil cover ( Figure 6). The basic statistics for each band for all satellite images have been evaluated. Each band was reclassified in two categories: a) for all pixels within the range of <reflectance>+/-σ and b) for all the pixels outside the specific range. As a result, binary files were created and Boolean addition in GIS environment was followed to produce a final classification map (Figure 7). After the creation of the spectral signature modeling map, 64 settlements in a total of 120 (56.6%) were established in areas of very high possibility.

Principal Component Analysis
Principal Component Analysis involves a mathematical procedure that transforms a number of correlated variables in a smaller number of uncorrelated variables called principal components. The method was applied to ASTER, Landsat and HYPERION images to decorrelate the data and to reduce the dimension of the study [31]. PCA of ASTER images concluded to the best results where 39 settlements were highly discriminated and 47 medium discriminated in a total of 247. Furthermore, 14 magoules that were not visible in the original images were clearly visible after applying PCA to ASTER images (Figure 8).

Data Fusion
Image fusion is a standard satellite image procedure of combining images of different spatial resolution to obtain a single final composite image. Image fusion is applied to digital imagery for different reasons such as to enhance certain features that are not visible in either of the single data alone [1] and to sharpen the images [32]. The images that can be used can be from different sensors and resolutions. By using ERDAS imagine software various fusion combinations and techniques were tried, such as ASTER (15 m) visible channels with the PCA product (PC1) of HYPERION (30 m) or the high resolution (1 m) bands (datafusion products) of IKONOS with the PCA product (PC1) of the HYPERION. PC1 of HYPERION image was selected in order to exploit the best radiometric resolution available compared to the rest high spatial resolution images. The results were highly promising for the cases of fusion (using PCA technique, namely re-scaling the high resolution image to fit the data range of PC1 following an inverse PC transformation, and cubic convolution interpolation) between high spatial resolution and high spectral resolution images (Figure 9).

Spectral Mixer Utility
In our effort to exploit the high spectral resolution of HYPERION images, a spectral mixer application through the use of Erdas Imagine 9.1 software was also applied. Spectral Mixer produces three bands to be assigned to the red, green, and blue color guns, but in this case instead of just assigning each band to a color gun one can select a weighted average of spectral bands to be assigned to a color gun [33]. For HYPERION images only the bands that had reflectance values above 0.3 were chosen and a weighting coefficient of 0.14 was applied for each band. The new RGB that was created (RGB1) employed the mixing of the bands (

Radiometric Enhancement
Due to the variable quality of the original images, the radiometric enhancement was vital for the appearance of the images and the better recognition of the terrain features. After applying radiometric enhancement to ASTER images (acquisition date of 19-03-2003) 57 settlements were detected. A non-linear radiometric enhancement of the HYPERION PCA image, followed by an inversion of brightness was able to highlight eight settlements from a total of nine. (Melia 1, Melia 2, Anagennisi 2, Moshohori 3, Kipseli 2, Prodromos 1 of Larisa, Nikaia 17 and Kuparissia 2). Similar type of non-linear radiometric enhancement of the high resolution IKONOS images through the modification of the histogram was able to outline the round shape of known magoules, as well as to identify 10 more targets of similar geometry that need to be verified by the ground truthing activities that will follow ( Figure 11).  The fact that different kinds of marks, such as crop, soil and shadow marks, are generally associated with the presence of buried archaeological remains [34][35][36] was exploited at the IKONOS images so as to detect some completely flat magoules such as Anagennisi 2 (Figure 12). Soil and moisture differences within near-surface archaeological deposits can influence surface vegetation patterns creating crop marks of various kinds. In addition soil marks can appear as changes in color or texture in freshly ploughed fields before the growing crops mask the surface of the soil [37].

Land Classification and Vegetation Indices
In most cases, difficulties in the detection of archaeological sites originate due to the fact that the spectral response of archaeological sites and surrounding areas is almost the same [2,14,38]. However, in the domain of predictive modeling, the specification of the environmental attributes that correlate to the location of the archaeological sites is of importance. For this reason, in order to investigate the regime of the land use surrounding the magoules, several methods of supervised classification were applied to Landsat and ASTER images. For the classification procedure five classes were defined: Uncovered land, Uncultivated land, Cultivated land, Urban area and Water reservoirs. Mahalanobis fuzzy classification proved to be the most efficient one in terms of the overall accuracy assessment (based on the error matrix) compared to all the classification algorithms that were applied (Maximum Likelihood, Minimum Distance, Mahalanobis Distance, Parallelepiped, Spectral Angle Mapper, Maximum Likelihood (fuzzy), Minimum Distance (fuzzy), Mahalanobis Distance (fuzzy)) ( Table 2 and Figure 13).  Due to the small agreement between the land use classification results that produced between Landsat and ASTER sensors, the Normalised Difference Vegetation Index (NDVI) was computed to analyse the difference of vegetation during various acquisition dates. Vegetation indices are mainly extracted from reflectance data from the red and near infrared (NIR) bands [39]. The NDVI was obtained by the following equation (3): As expected, the NDVI of the spring ASTER image was higher than the summer Landsat image (Figure 14).

De-correlation Stretch
The de-correlation stretch is a process that is used to enhance (stretch) the color differences found in the input pixels. The principal component transformation is similar, except the fact that the transformation vectors are derived from the correlation matrix rather than the covariance matrix. Decorrelation stretch to the ASTER images managed not only to detect easily 36 Neolithic settlements ( Figure 15), but also to estimate the area of each settlement in GIS environment.

Spatial Enhancement
Spatial enhancement of images is considered to be a standard satellite image enhancement. In order to emphasize the marks arising from the presence of magoules to Thessaly plain various spatial filters were applied to all the images. Of the several types of filters that were applied in the specific study, only two of them, Sobel Right Diagonal 3x3 and Laplace 3x3, proved to be very useful for the detection of Neolithic settlements. Although the spatial filters were applied to all bands of ASTER and the first three principal components of HYPERION images, they were especially satisfactory when they were applied at the first band of ASTER image ( Figure 16).  The values of the matrices can be seen in Table 3. The extraction of statistics about the number of settlements that were detected by each filter indicated that Sobel right diagonal filter was the most reliable one achieving a discrimination of almost 150 sites (Table 4).

Object Based Remote Sensing
The ASTER image was segmented and classified based on an object based approach through the use of e-Cognition software. The object based technique is considered as very useful for heterogeneous land covers [14]. Segmentation is the most important phase in object based classification. The image is subdivided to homogeneous areas based on their spatial characteristics, shape, scale and object hierarchy level [40]. The second phase includes the classification of image, where training objects are selected to train the classification in a similar way to the pixel based classification but instead of using pixels as training samples, geometric objects are used. Subsequently, classification parameters are defined [14]. The application of object -oriented methodology to the ASTER images managed to detect easily only 15 settlements in a total of 234, whereas 185 settlements were not discriminated at all and 34 were medium discriminated from the neighbor pattern. For the application of segmentation to ASTER image we used a scale factor of 5. However, the fact that the settlements don't have uniform shape and spatial characteristics was the main reason for the poor results of this methodology ( Figure  17).

Predictive Modeling
After applying all the above enhancement processes a predictive model was designed to locate potential magoules in the wider region of the Thessaly plain. The results of land use classification, NDVI estimates and those from the spectral signatures and classification of the ASTER image (acquisition date 19-03-2003) were combined together with a DEM constructed by digitization of 1:50.000 scale topographic maps. All these data were reclassified and a certain weight factor was applied to each cell of the raster layers. The weighting and rating factors were specified based on the statistical analysis of the specific parameters in relation to the correlation of them with the known magoules and their importance in terms of the location of the magoules. All the raster layers were rated ( Table 5) and equation (4)   The final predictive map consisted of pixel areas with different probability for the existence of Neolithic settlements. It was estimated that 92 of the already known settlements are laid on areas of high probability and 23 in areas of medium probability.

Conclusions
The various approaches applied on different satellite images for the detection of Neolithic settlements in Thessaly illustrate the benefits that satellite remote sensing can provide in archaeological investigations. It was proven that an integration of images from different satellite sensors can contribute to a faster and more accurate and qualitative detection of archaeological sites.
Specifically, ASTER images proved to be the most reliable and efficient for the detection of Neolithic settlements, being able to combine a medium spatial resolution with high spectral resolution. In contrast, Landsat images concluded to quite poor results, mainly due to the acquisition date of the imagery, which produced low signal to noise ratio for the archaeological targets. The high spectral abilities of HYPERION especially after merging it with the high resolution images of IKONOS seem to have an increased potential not only in detecting but also in outlining the particular features. The image processes that proved to be more effective were the spatial filtering, the process of decorrelation stretch and the radiometric enhancement. The integration of land use classification data with NDVI and spectral signatures resulted to very promising modeling maps. On the other hand, the object based classification method proved that most of Neolithic mounds lack uniform shape characteristics that can be easily distinguished from the surrounding vegetation patterns. Although most of them have a circular or oval shape, they belong to the same land use type of the wider region that makes them almost impossible to separate from the other features of the terrain. Furthermore, although the use of conventional aerial photos can often pinpoint particular features based on crop marks, the intensive use of space and the relatively leveling of the ground at the magoules has masked their particular features. This makes them impossible to be identified without the exploitation of the enhanced radiometric resolution of satellite imagery. Thus, satellite remote sensing may offer further advantages with other type of archaeological targets, and it offers potential for further investigation. The vegetational regime at the mounds proved to be a crucial factor for their detection. In case there had been different kinds of vegetation on the settlements and the surrounding areas, the automatic extraction by means of remote sensing would have been easier.
The above processes were limited to the satellite imagery. The particular methods can be also employed for the detection and mapping of similar archaeological targets such as Bronze Age mounds and settlements, monumental tholos tombs and others. The results of this study can be further enhanced through manipulation of the above conclusions with the spatial tools of GIS applied to the distribution of the magoules on the geomorphologic attributes of the terrain. In this way, a more integrated and synthetic tool for the detection of the magoules and the study of the Neolithic settlement patterns can be produced.