Download

A- A+
Alt. Display

# Detecting Looting Activity through Earth Observation Multi-Temporal Analysis over the Archaeological Site of Apamea (Syria) during 2011–2012

## Abstract

This study aims to investigate the potentials of open-access, freely distributed Earth Observation images for detecting large-scale looted areas. The analysis was conducted using medium-resolution Landsat 7 ETM+ images over the archaeological site of Apamea, at Syria. The site has been systematically and extensively looted in the recent past, attracting the interest of scholars.

We propose a multi-temporal analysis of cloud-free multispectral Landsat 7 ETM+ images throughout the period between January 2011 to April 2012, just at the beginning of the Syrian war. The analysis was completed through the interpretation of pseudo-color temporal composites, investigation of the multi-temporal spectral profiles, correlations between the spectral bands, and the application of principal component analysis (PCA). The overall analysis was limited within the spectral range of 450–1750 nm. This wavelength range corresponds to the first five spectral bands of the Landsat images. Furthermore, we explored the big data cloud platform Google Earth Engine to detected looted areas. A supervised classification strategy was designed and performed on this cloud platform employing the Random Forest classifier. Finally, a time-stamp change detection approach was implemented.

The overall findings were compared with available images from Google Earth Digital Globe and published articles and reports related to the Apamea archaeological site. It was found that the high revisit temporal resolution of the Landsat sensor was able to detect and map the looting activity in the area as a result of the spectral change in the archaeological landscape, despite its 30 m spatial resolution. At the same time, however, the analysis has provided other false-true detections in other areas in the landscape.

Keywords:
How to Cite: Agapiou, A., 2020. Detecting Looting Activity through Earth Observation Multi-Temporal Analysis over the Archaeological Site of Apamea (Syria) during 2011–2012. Journal of Computer Applications in Archaeology, 3(1), pp.219–237. DOI: http://doi.org/10.5334/jcaa.56
Published on 11 Aug 2020
Accepted on 27 May 2020            Submitted on 15 Feb 2020

## 1. Introduction

Increased interest from the scientific community for monitoring archaeological areas through satellite imagery has been reported in the recent past (Luo et al. 2019; Stewart, Oren & Cohen-Sasson 2018; Abate & Lasaponara 2019; Danti, Branting & Penacho 2017; Tapete & Cigna 2017; Chen, Lasaponara & Masini 2017; Agapiou & Lysandrou 2015). These studies have shown the benefits of using Earth Observation as a systematic tool for monitoring natural and man-made hazards (Tzouvaras et al. 2019), including also looted areas (Cigna et al. 2014; Tapete et al. 2013).

A recent review study by Tapete and Cigna (2019) properly summarizes the current status of remote sensing methods for detecting looting activity around the world, using either simple (visual) or more advanced processing techniques. These findings showcase that researchers, scholars and other interested parties have turned their interest to explore spaceborne sensors for monitoring looting activities. Despite the fact that satellite observations cannot prevent the illegal actions on the ground, the identification of new looted areas, probably unknown to local stakeholders, is considered as a critical step towards the increase of awareness for potential illegal trafficking (Parcak 2017; Parcak et al. 2016). Nevertheless, the overall research using satellite sensors can be restricted by their accessibility within the specific time-window of the looted activity.

The detection of looting marks using high resolution optical and radar images has been widely discussed in the recent past (Caspari 2018; Lauricella et al. 2017; Agapiou, Lysandrou & Hadjimitsis 2017; Casana 2015; Casana & Panahipour 2014), using different types of sensors such as the WorldView, GeoEye and the IKONOS products. Nevertheless, the exploitation of open-access, freely distributed medium resolution datasets is still limited considered in the literature, mainly due to pixel resolution. Indeed, from the relevant literature, it is evident that the use of medium resolution optical satellite datasets – like those of Sentinel-2 and Landsat series – is still not sufficiently investigated. As Tapete and Cigna (2019) highlight in their work, only one study has demonstrated the potential use of Sentinel-2 image (i.e. the work of Tapete & Cigna 2018) for detecting looted areas. As the authors argue (Tapete & Cigna 2018, 3.1.2), “the use of ….data — such as free-of-charge Sentinel-2 and Landsat imagery — is still at the very early stage. Nevertheless, it is to explore at what extent these datasets can be utilized for regional mapping, owing to their extensive spatial coverage per single frame and availability for the entire landmass, as well as their short revisit time” and we here adopt this position.

The short revisit time of open-access, freely distributed satellite images can be a critical factor for detecting illegal ground activities. Consequently, their capabilities need to be further studied in order to understand their potentials for monitoring large scale archaeological sites. As Li & Roy (2017) argue, Landsat-8, Sentinel-2A, and Sentinel-2B optical sensors together can provide a global median average revisit of 2.9 days, and, over a year, a global median minimum revisits 14 min (±1 min) and maximum revisit interval of 7.0 days. Of course, the main barrier using these medium resolution datasets is their coarse resolution compared to other existing higher spatial resolution satellite data.

In this study, we focus on the archaeological site of Apamea (Syria) during the period 2011–2012, aiming to capture (known) looting activity from its beginning. For this purpose, simple and more advanced image processing methodologies are tested, presented and discussed evaluating thus potentials and limitations of medium resolution images for monitoring sites under threat.

The paper is organized as follow: at first, a short introduction of the Apamea archaeological site and the looting activity carried out in the last years is presented. Then the methodology adopted for the study as well as the datasets used are given, followed by the results section. The study ends with the conclusion section; summarizing highlights of the current investigation.

## 2. Looting activities in the Apamea region

The archaeological site of Apamea, in northwest Syria (Lat: 35.42°; Long: 36.40°) is a city founded in 300 BC that holds impressive archaeological remains. The site includes the Great Colonnade which ran for nearly 2 km, and the Roman Theatre, one of the largest surviving theatres of the Roman Empire. Recent excavation reports for the site can be found in Finlayson (2012) while further archaeological literature can be found in Balty (1969).

Since 2010, the site drew the attention of scholars due to the looting activities started after the beginning of the Syrian world, in March 2011. In parallel, resulting from the extensive and systematic looting activities, the site has been popularized in the media as well. The destruction of the site, visible in high-resolution satellite images, revealed ongoing illegal excavation activities taking place on the eastern, northeastern, and western regions of the city. As Casana & Panahipour (2014) argue, the looting trenches of the site typically measure up to 3m. These illegal excavations were taking place probably by using bulldozers and other heavy machinery, based on the size, pattern, as well as the number of the looting holes (Casana & Panahipour 2014). This observation was also reported earlier by UNESCO (Observatory of Syrian Cultural Heritage 2010).

The illegal activities in the Apamea site were investigated by several scholars using high-resolution satellite images, like the radar sensor TerraSAR-X Staring Spotlight (ST) covering the period October 22nd 2014 to June 21st 2015, with an azimuth resolution of up to 0.24 m (Tapete, Cigna & Donoghue 2016). Other researchers, such as Casana & Panahipour (2014), have reported looting activities at the site using a limited number of very high-resolution images with sub-meter analysis. In specific, Casana & Panahipour (2014) used WorldView imageries acquired on June 2nd 2008, November 28th 2013 and March 6th 2014. Their results have demonstrated that most of the looting activity could be dated to the first phase of the conflict, having occurred by April 2012, although it continued into November 2012. They also found that in April 2012 only a few holes were dug into the privately-owned land, while a considerable number of holes were observed outside the walls. Unfortunately, the site was repeatedly looted after that period, and this was confirmed throughout the interpretation and analysis of recent imageries (Tapete & Cigna 2019, Table 2).

The visual inspection of the looting activities of the Apamea archaeological site during the years 2011–2012 can be performed through the compressed red-green-blue (RGB) images of the Google Earth (Figure 1). Figure 1 (top), indicates no looting activity at the archaeological site on July 20th 2011, but extensive looting activities in the area can be depicted in Figure 1 (bottom) taken on April 04th 2012.

Figure 1

Apamea archaeological site on July 20th 2011 (top) and April 04th 2012 (bottom). Looted areas are highlighted on the figure April 04th 2012 (bottom) with a red polygon while the densest looting areas are visible with a blue color in the Figure (images from Google Earth Digital Globe).

## 3. Datasets and Methodology

### 3.1. Datasets

#### 3.1.1. Landsat 7 ETM+ datasets

For the needs of the study, 14 Landsat 7 ETM+ images were used covering the period from January 2011 until April 2012. The images were provided by the Landsat NASA/USGS space program. The acquisition dates of the images are shown in Table 1. The images were a-priory selected to have minimum cloud-coverage (less than 10%).

Table 1

Details of the Landsat 7 ETM+ datasets used for the aims of the study.

no date Period* Notes**

1. 2011-01-12 1st period (T0) Before looting events
2. 2011-03-01 Before looting events
3. 2011-06-21 Before looting events
4. 2011-07-07 Before looting events
5. 2011-07-23 2nd period (T1) Looting period
6. 2011-08-08 Looting period
7. 2011-08-24 Looting period
8. 2011-09-09 Looting period
9. 2011-09-25 Looting period
10. 2011-10-11 Looting period
11. 2011-10-27 Looting period
12. 2011-11-28 Looting period
13. 2012-03-19 Looting period
14. 2012-04-04 Looting period

The images, after the interpretation of the Google Earth software’s data (Google Earth 2020), were categorized into two periods based on the presence of looting activity. The first period, T0, refers to the phase before looting activity, while the second period, T1, refers to the looting period were illegal excavations are visible. As mentioned before, after the T1 period, other looting events were also reported from other researchers (see introduction). Four Landsat images were used for the T0 period, while for the T1 period, ten medium Landsat resolution images were used as indicated in Table 1.

#### 3.1.2. Other available data

Compressed RGB high-resolution images from the Google Earth platform were used as ground truth (see Figure 1) to confirm the overall results from the analysis of the Landsat datasets. Two high-resolution images were retrieved for the time-window of our analysis (2011–2012), namely an image with an acquisition date of July 20th 2011 (used for the T0 period) and another image on April 04th 2012 (used for the T1 period).

The analysis was also validated with existing published material from colleagues that worked in the Apamea area, as well as other online material from different organizations and groups. Such kind of data can be found in the works of Tapete & Cigna (2016, 2019).

### 3.2. Methodology

For the needs of the study, a multi-temporal analysis of the satellite images was followed. Despite the medium spatial resolution of the Landsat series, this study aims to help us to understand their use for monitoring purposes of archaeological sites and landscapes (Agapiou, Alexakis & Hadjimitsis 2019), as is the case of illegal looting activity of the Apamea site. The Landsat 7 ETM+ sensor was selected since other open-access satellite datasets like the Sentinel 2A and 2B sensors of the Copernicus Programme having better spatial resolution (10m and 20m resolution products) are available only after 2015. Also, the newest Landsat sensor, namely the Landsat 8 LDCM, was launched in 2013. Therefore, for the specific needs and the period of observation of the current study, no alternative option in terms of freely distributed datasets was available.

Initially, the images were selected and downloaded, based on their cloud-coverage characteristics over the area of interest. All images were already pre-processed by the provider, in terms of radiometric and geometric corrections. The geo-registration of all Landsat scenes was consistent, and within prescribed image-to-image tolerances of ≦ 12-meter radial, root means square error (RMSE) (Landsat Collections 2020). Regarding the radiometric corrections, the atmospheric effects were taken into consideration. The images were downloaded at bottom-of-atmosphere reflectance values minimizing thus radiometric errors due to the sun elevation angle, sensor bias and atmospheric scattering.

Subsequently, the individual spectral bands were stacked together, creating a n-dimensional spectral space of the Apamea region, including Landsat images before and during the reported looting period (T0 and T1 periods, respectively). Here we process only the first five bands of the sensor, namely the blue (band 1: 450–520 nm), green (band 2: 450–520 nm), red (band 3: 630–690 nm), near-infrared (band 4: 770–900 nm, NIR) and short-wave infrared (band 5: 1550–1750nm, SWIR) bands. The sensor records the target-leaving radiance from the visible to the short-wave infrared parts of the spectrum, which are considered very important for monitoring soil and vegetated areas as well as for land use/land cover changes. In specific, the NIR spectrum can be used to detect vegetated areas due to the response of the healthy vegetation to reflect a large amount of the NIR incoming radiance. At the same time, it absorbs much of the energy in the red part of the spectrum. The SWIR part of the spectrum is used to discriminate soil types and geological formations, while it is also used to detect soil moisture.

Initially, the spectral profile of the looted area was extracted for the first five bands of the Landsat 7 ETM+ images using the Google Earth Engine environment. In addition to the spectral profiles, the Normalized Difference Vegetation Index – NDVI – (equation 1) was estimated to monitor any vegetation dynamics of the area.

(Eq. 1)

Where the ρNIR refers to the reflectance at the NIR part of the spectrum (band 4) and ρRED the reflectance at the red band (band 3). Additionally, various pseudo-color composites were created. This first-step investigation allowed us to observe significant spectral changes of the area, but not the detection of the individual looting holes of the area (due to the spatial constraints of the sensor). In addition, utilizing the integrated Landsat 7 ETM+ dataset, various pseudo color composites were visualized, and interpreted.

Then a Principal Component Analysis (PCA) was implemented. The PCA is a statistical analysis process of the image that accounts the variation of the pixel values within the image (Campbell 2007). This orthogonal transformation re-projects the spectral bands into a new n-dimensional space of linearly uncorrelated variables. While PCA is applied in single image approaches for target detection, it can also be applied in a temporal dataset. Once the PCA is applied in the temporal n-space, the analysis takes into consideration the overall temporal variance of the images. This type of analysis can be performed for images over the same area and minimum radiometric noises. In our case, this is partially true since the geometric rectification of the Landsat 7 ETM+ archive was less than 0.5 of pixel size, while at the same time, the radiometric variation of the incoming radiance and sun elevation angles have been used to normalize the overall results. Therefore, the first principal components (PCs) can better explain any significant seasonal changes of the area. The following PCs can be used to interpret additional temporal variations of the area, such as those of the looting activity.

Furthermore, the Pearson Correlation coefficient (equation 2) and the Mahalanobis spectral distance (equation 3) were used to calculate the spectral similarities and spectral distance of the Landsat’s bands.

(Eq. 2)
$Rx,y=\frac{cov\left(X,Y\right)}{{\sigma }_{x}{\sigma }_{y}}$

Whereas R is the Pearson’s correlation coefficient, cov is the covariance, σx is the standard deviation of x and σy is the standard deviation of y.

(Eq. 3)

Where Y is the Mahalanobis spectral distance, Yi is the data value vector at row I, is the mean vector and S–1 is inverse of the covariance matrix. Based on the findings of the Mahalanobis spectral distance and Pearson’s correlation values, the normalized ratios for some interesting Landsat pairs were calculated. The normalized ratio was estimated using the ‘Max Spectral Distance’ equation (equation 4).

(Eq. 4)

Moreover, supervised pixel-based classification algorithms were tested in the area using the Google Earth Engine big data cloud platform. The specific environment allows us to explore the capabilities of the Google server infrastructure, minimizing thus the time and computer resources for the image processing (Agapiou 2017; Orengo & Petrie 2017). Several supervised classifiers were fitted with specific training areas, namely the ‘looted areas’, ‘vegetation’, ‘soil’, ‘water’, ‘urban’ and ‘roads’. Using these training areas, the classification results were obtained, evaluated, and compared with Google Earth images. A statistical analysis of the classification process was performed based on the classification confusion matrix.

Finally, we designed a time-stamp change detection approach using the Landsat images. This approach uses consecutive pairs of the Landsat images of Table 1. In specific, we explore the rate of spectral change based on the following equation (Eq. 5):

(Eq. 5)

Whereas ΔR/R is the rate of spectral change, Ri,j is the reflectance value of the band (i) at a given date of satellite overpass (j) while Ri,j–1 refer to the reflectance value of the band (i) of the previous satellite overpass. The overall results are presented in the section below.

## 4. Results

### 4.1. Spectral profiles and pseudo-color composites

The temporal changes of the spectral bands during the T0 to T1 periods are shown in Figure 2. Figure 2a to 2e indicate the spectral reflectance of the bands 1 to 5, respectively. The NDVI is also plotted, indicating the vegetation dynamics of the area. Precipitation values (in mm per day) are plotted in Figure 2 with a blue color (data from NASA 2020).

Figure 2

Spectral profile over the looted area for the first five bands (band 1 at Figure 2a; band 2 at Figure 2b; band 3 at Figure 2c; band 4 at Figure 2d and band 5 at Figure 2e).

While a remarkable increase of the reflectance value at the near-infrared part of the spectrum (band 4, Figure 2d, see arrow) is observed at the beginning of March 2011 (image no. 2 of Table 1), this should be linked with vegetation growth (see NDVI increase during this period). Beyond this period, the reflectance values remain stable for all bands until the end of October 2011 (image no. 11 of Table 1). From this point (image no. 12 of Table 1, 2011-11-28) a decrease of reflectance is observed, and especially for the SWIR part of the spectrum (see the arrow at Figure 2e). A decrease of reflectance is likely to be noticed after a rainfall, which is not, however, our case. Therefore, this period between the 2011-10-27 and 2011-11-28 indicates a spectral change of the soil which cannot be explained either from the vegetation dynamics (NDVI) or the precipitation of the area.

Figure 3a presents the R-G-B pseudo color composite for the T0 period (date of acquisition: 2011-03-01, no.2. of Table 1). In contrast, Figure 3b and 3c present the NIR-R-G and the SWIR-NIR-R pseudo color composites respectively for the same date. Similarly, Figure 3d to 3f present the pseudo color composites for the Landsat image acquired on 2011-11-28 (see no.12. of Table 1). While the spatial resolution of the Landsat image inhibits a detail interpretation of both images, some interesting observations can be made for these pseudo-color composites.

Figure 3

(a) R-G-B pseudo color composite of T0 period (date of acquisition: 2011-03-01, no.2. of Table 1); (b) NIR-R-G pseudo color composite of T0 period (date of acquisition: 2011-03-01, no.2. of Table 1); (c) SWIR-NIR-R pseudo color composite of T0 period (date of acquisition: 2011-03-01, no.2. of Table 1); (d) R-G-B pseudo color composite of T1 period (date of acquisition: 2011-11-28, no.12. of Table 1); (e) NIR-R-G pseudo color composite of T1 period (date of acquisition: 2011-11-28, no.12. of Table 1); (f) SWIR-NIR-R pseudo color composite of T1 period (date of acquisition: 2011-11-28, no.12. of Table 1).

The most interesting outcome is the presence of vegetation within the archaeological site, which is represented in red colour in Figure 3e, indicating high reflectance values in the near-infrared part of the spectrum. The shape and the size of this vegetated area are also alike with the one highlighted as a looted area in Figure 1 (see red polygon Figure 1, bottom). Indeed, while vegetation is evident within the archaeological site and its surrounding during other periods (see Figure 2b), both the size and the shape of the vegetation observed at the image taken at 2011-11-28 (Figure 2e) correspond to the boundaries of the looting activity of the site (Figure 1, bottom). New vegetation is normally grown in areas whereas the top-soil has been disturbed, as the case of a looting activity. However, as expected, the detection of individual looting holes is not feasible as this is the primary constraint of the resolution of the Landsat 7 ETM+ images (30 m pixel resolution).

Further image statistics were elaborated using Pearson’s correlation value (R2). Each spectral band was correlated with the rest of the bands of all Landsat images, creating thus a 14 × 14 correlation matrix, and in total 4,900 Pearson’s correlation coefficient values (= 5 bands × 14 × 14 images). The results are shown and visualized in Tables 2, 3, 4, 5, 6, for bands 1 to 5, respectively.

Table 2

Pearson’s correlation value (R2) for the spectral band 1. The T0 period is indicated with a red polygon (no.1–no.4).

Table 3

Pearson’s correlation value (R2) for the spectral band 2. The T0 period is indicated with a red polygon (no.1–no.4).

Table 4

Pearson’s correlation value (R2) for the spectral band 3. The T0 period is indicated with a red polygon (no.1–no.4).

Table 5

Pearson’s correlation value (R2) for the spectral band 4. The T0 period is indicated with a red polygon (no.1–no.4).

Table 6

Pearson’s correlation value (R2) for the spectral band 5. The T0 period is indicated with a red polygon (no.1–no.4).

The range value of the Pearson’s correlation value is between –1 to +1. Values close to 0, indicated with light yellow colour in Tables 2, 3, 4, 5, 6, highlight spectral bands with no correlation (→0 correlation coefficient), while green colour shows a strong (negative or positive respectively) correlation (→1 correlation coefficient). The diagonal values are equal to 1 (the same band correlated to itself), while the lowest Pearson’s correlation value (R2) is highlighted with the blue font at each table. For the Apamea case study, our interest is focused on combination of bands from both T0 and T1 periods, with regression value close to 0 (light yellow zones).

As it was found, the most interesting combination is the image pair no. 3 with no. 12 (2011-06-21 and 2011-11-28 respectively) at the blue, green and SWIR bands, giving an R2 value score equal to 0.25, –0.04 and 0.42. For the red and NIR bands, the most noteworthy pairs are no. 3 and no. 10 (2011-06-21 and 2011-10-11), and no.3 and no. 5 (2011-06-21 and 2011-07-23), respectively. In addition, an interesting pair is the one of images no. 1 and no. 5 (2011-01-12 and 2011-07-23).

Further to the correlation analysis, the Mahalanobis spectral distance of the bands was estimated. Tables 7, 8, 9, 10, 11 show the results of the Mahalanobis separability index for the spectral bands 1 to 5, respectively. Higher values, indicated with green colour, emphasize higher spectral distance between the pair of bands. In contrast, lower values, indicated with light blue color, suggest no significant spectral distance. Our interest is focused at pairs with high spectral separability, especially those of periods T0 and T1. The highest spectral value for each band is highlighted with a blue font in Tables 7, 8, 9, 10, 11. For both bands 1 and 2, the highest spectral distance was found for image no.1 and no. 6 (see Table 1), while for bands 3 to 5, the best combination was found for the image no. 1 and no. 5; no. 1 and no. 12, and no. 4 with no. 8, respectively.

Table 7

Mahalanobis distance for the band 1. Higher values indicate higher separability between the pair-wise bands. The T0 period is indicated with a red polygon (no.1–no.4).

Table 8

Mahalanobis distance for the band 2. Higher values indicate higher separability between the pair-wise bands. The T0 period is indicated with a red polygon (no.1–no.4).

Table 9

Mahalanobis distance for the band 3. Higher values indicate higher separability between the pair-wise bands. The T0 period is indicated with a red polygon (no.1–no.4).

Table 10

Mahalanobis distance for the band 4. Higher values indicate higher separability between the pair-wise bands. The T0 period is indicated with a red polygon (no.1–no.4).

Table 11

Mahalanobis distance for the band 5. Higher values indicate higher separability between the pair-wise bands. The T0 period is indicated with a red polygon (no.1–no.4).

Following the Mahalanobis spectral distance and the Pearson’s correlation values, the ‘Max Spectral Distance’ (see equation 4) was estimated for the images: (i) no. 3 with no. 12 for band 1, (ii) no.1 with no. 6 for band 1, (iii) no. 3 with no. 12 for band 5 and (iv) no. 4 with no. 8 for band 5.

The results after the application of the ‘Max Spectral Distance’ in these four pairs of images are shown in Figure 4. A Google Earth image over the Apamea area (before and after the looting event) is displayed in Figure 4a and 4d, respectively. Figure 4b shows the ‘Max Spectral Distance’ of images no. 3 and no. 12 of Table 1 (2011-06-21 and 2011-11-28). Similarly, Figure 4c shows the results from the normalized ratio of no.1 with no. 6 (2011-01-12 and 2011-08-08) using the blue band. The last two figures (Figure 4e and 4f) demonstrate the results for the SWIR band for images no. 3 with no. 12 (2011-06-21 and 2011-11-28), and no. 4 with no. 8 (2011-07-07 and 2011-09-09). Substantial spectral changes are highlighted with red color. While spectral changes are observed over the looted areas (Figure 4d), the use of a single pair of spectral bands cannot distinguish the exact boundaries of the looted area. Indeed, substantial spectral differences are evident in the whole area, which can be linked with seasonal or land-use changes. The normalized ratio of no.1 with no. 6 (2011-01-12 and 2011-08-08) (Figure 4c) and no. 3 with no. 12 (2011-06-21 and 2011-11-28) (Figure 4e), have similar errors and cannot support the reliable extraction of the area characterized with looted activity.

Figure 4

High-resolution image from Google Earth of the area of Apamea before (a) and after the looting event (d). The ‘Max Spectral Distance’ normalized ratio of the blue bands of images no. 3 and no. 12 (b) and no.1 and no. 6 (c) are shown on the first row. (e) and (f) show the results for the SWIR band of the pairs no. 3-no. 12 no. 4-no. 8. Changes are highlighted with red color while no significant changes with blue color.

### 4.2. PCA analysis

To further process the spectral cube of the Landsat datasets (14 Landsat 7 ETM+ datasets, 70 spectral bands), PCA was applied. The analysis was performed to map any temporal changes of the area based on the spectral variance of the images. The latest can be linked with the spectral heterogeneity of the overall spectral cube in the given area. In our study, this spectral variance can be linked with the looting activity over the archaeological landscape of Apamea, once however the seasonal spectral variations (e.g. vegetation dynamics) and other land-use changes are removed. The results are shown in Figure 5.

Figure 5

High-resolution image from Google Earth of the area of Apamea before (a) and after the looting event (d). (b) first principal component – PC1; (c) second principal component – PC2; (e) third principal component – PC3 and (f) fourth principal component – PC4. The looted area is highlighted with a red polygon in Figure 5f.

Figure 5a and 5d present the archaeological site of Apamea before and after the looting activity, while the Figure 5b, 5c, 5e and 5f present the first to fourth principal components (PC1-PC4) outcomes, respectively. The light grayscale tones of the Figure 5b, 5c, 5e and 5d indicate areas with high spectral variance scores. In Figure 5b, this variance is evident in the agricultural areas in the northeastern and north-western part of the archaeological site (see red circles in Figure 5b). In comparison, high variation is recorded in the water body at the eastern part of the site (see red arrows in Figure 5b). This observation can be justified by the seasonal changes of the water body and its quality, as well as the various phenological stages of the cultivated areas. These areas are also highlighted in the second principal component (see arrows in Figure 5c). The agricultural area also dominates in the third principal component (see the circle in Figure 5d); however, some variations begin to be visible in the archaeological site per se. Therefore, the first three PCs (PC1–PC3) do not seem to be able to capture the looting activity of the area, mainly due to the seasonal changes of the landscape. This, as mentioned earlier, was something expectable since temporal changes can influence the overall spectral cube significantly and give heterogenous spectral profiles through time. In contrary, the fourth PC (PC4, Figure 5f) was able to enhance better the looting activity. The red polygon shown in Figure 5f indicates the looted area (see also Figure 5d) which fits very well with the outlines of the high variance (of the fourth PC). At the same time, we can also notice the dark tones (low spectral variance) within the archaeological site, whereas no looting was reported (see red arrows at Figure 5f).

### 4.3. Supervised classification

Various supervised classifiers have then been implemented through the Google Earth Engine. The Landsat 7 ETM+ image taken at 2011-11-28 was used for the classification process, based on the outcomes of section 4.1. Training areas for the supervised classification process have been selected, after a careful interpretation of the image. As mentioned earlier, the following categories have been used for this classification process: ‘looted areas’, ‘vegetation’, ‘soil’, ‘water’ and ‘urban’ areas. Similarly, other polygons have been selected for the same classes for validation purposes. The exact number of pixels selected for each class can be found in Table 12. The classification results using the Random Forest (RF) classifier are depicted in Figure 6. The overall accuracy of the RF classifier was estimated to be more than 90%, which was, however, reached after several repetitions of the training model to minimize misclassification errors. The confusion matrix of the RF classifier is provided in Table 12. The classification results per class can be seen under the columns class (%) of Table 12.

Figure 6

(a) R-G-B pseudo colour composite of the Landsat 7 ETM+ taken at 2011-11-28; (b) NIR-R-G pseudo color composite of the Landsat 7 ETM+ taken at 2011-11-28; (c) training areas used for the classification purposes on top of a high-resolution satellite base-map provided by ESRI ArcGIS Online service (d) RF classification results and (e) RF classification results after the application of a majority filter.

Table 12

Confusion matrix of the Random Forest classification algorithm implemented at the Landsat 7 ETM+ image taken at 2011-11-28 (results refer to the classification process (Class.) and post-classification analysis (Test).

Class Name # Points Looted areas Vegetation Water Soil Urban

Class. Test Class. (%) Test (%) Class. (%) Test (%) Class. (%) Test (%) Class. (%) Test (%) Class. (%) Test (%)

Looted areas 143 42 99.3 86 0.7 7 0 0 0 2 0 5
Vegetation 211 45 0 4 99.53 96 0 0 0.47 0 0 0
Water 32 11 0 0 0 9 100 91 0 0 0 0
Soil 779 131 0.92 4 0 0 0 0 99.08 89 0 7
Urban 291 17 0 0 0 0 0 0 3.29 8 96.71 92

Figure 6a and 6b show the Landsat 7 ETM+ image at the R-G-B and NIR-R-G pseudo colour composites. The archaeological site of Apamea is also indicated with an arrow. Figure 6c shows the various training areas used for the classification purposes on top of a high-resolution satellite image, provided by ESRI ArcGIS Online base-maps service. The RF classification results at a 30-meter resolution using are depicted in Figure 6d. Figure 6e shows the same classification results after the implementation of a 3 × 3 majority filter for minimizing the ‘salt and pepper’ classification error.

As depicted in Figure 6d and 6e, the RF classifier was able to classify with relatively good accuracy the looted area (indicated with red color) with limited false positives. Although other looted areas have been reported in the literature, this was not able to be confirmed by the Google Earth images, which were used as the ground truth datasets in our case study. The validation accuracy of this classification is indicated in Table 12 under the columns named test (%). The validation results are similar to the one of the classification analysis, indicating an overall high classification accuracy. Classification results were also generated with the application of the SVM classifier and other pixel-based classifiers, using different sampling strategies, which did not perform useful classification outputs (results not shown here).

### 4.4. Time-stamp change detection of the archaeological site

While the classification findings are encouraging in terms of the detection of looted, this was feasible only to our prior knowledge of the event (looting). Indeed, the supervised classification analysis at the Google Earth Engine was performed once we trained the RF classifier with specific areas of interest (i.e. looted areas).

For this reason, the time-stamp change detection approach, as mentioned before (see equation 5), was implemented. The following Figures (Figures 7, 8, 9, 10, 11) illustrate the rate of spectral change for the five spectral bands of the Landsat 7 ETM+ sensor, namely the blue (Figure 7), green (Figure 8), red (Figure 9), near-infrared (Figure 10) and middle infrared bands (Figure 11). The two different periods discussed in our study, namely the T0 and, T1, are indicated in the figures, while the transition period from T0 to T1 is illustrated with a blue rectangle.

Figure 7

ΔR/R rate of change over the archaeological site of Apamea (red line) at the blue band of Landsat 7 ETM+, in comparison with other areas covered with soil (yellow line), urban areas (purple line), vegetated regions (green line) and water bodies (blue line).

Figure 8

ΔR/R rate of change over the archaeological site of Apamea (red line) at the green band of Landsat 7 ETM+, in comparison with other areas covered with soil (yellow line), urban areas (purple line), vegetated regions (green line) and water bodies (blue line).

Figure 9

ΔR/R rate of change over the archaeological site of Apamea (red line) at the red band of Landsat 7 ETM+, in comparison with other areas covered with soil (yellow line), urban areas (purple line), vegetated regions (green line) and water bodies (blue line).

Figure 10

ΔR/R rate of change over the archaeological site of Apamea (red line) at the NIR band of Landsat 7 ETM+, in comparison with other areas covered with soil (yellow line), urban areas (purple line), vegetated regions (green line) and water bodies (blue line).

Figure 11

ΔR/R rate of change over the archaeological site of Apamea (red line) at the SWIR band of Landsat 7 ETM+, in comparison with other areas covered with soil (yellow line), urban areas (purple line), vegetated regions (green line) and water bodies (blue line).

This time-stamp approach can indicate reflectance changes between two consecutive Landsat images of the same part of the spectrum. Temporal spectral changes over ‘looted areas’ are shown with a red line in all graphs at Figures 7, 8, 9, 10, 11. Yellow, purple and green lines indicate the spectral changes over ‘soil’, ‘urban’ and ‘vegetation’ areas respectively.

Small differences are evident in all the five spectral bands. However, significant temporal spectral differences over the archaeological site are evident only in specific periods. For instance, at Figure 7 (band 1), we can observe a noticeable spectral change during the first months of 2011 for the ‘soil’ and ‘looted’ areas with an increase of the reflectance values. The change for both classes is the same, indicating that the archaeological area was still not disturbed. Both classes had a similar dramatic drop of their reflectance at the summer months of 2011, which the first significant spectral differences between these classes are reported after November of 2011 until the end of the monitoring period (April 2012). This observation is aligned with the previous findings which indicate that from 2011-11-28 (image no. 12 of Table 1) a decrease of reflectance can be observed (see also arrow at Figure 2e). The temporal change of the spectral signature of the urban areas remains similar within a spectral range of approximately 10%.

Comparable temporal changes for the ‘soil’ and ‘looted’ areas are also reported for the green (Figure 8) and red bands (Figure 9). At the same time, spectral changes are also noticeable in the SWIR part of the spectrum (Figure 11). Therefore, key spectral regions can be considered the blue (Figure 7) and SWIR (Figure 11) part of the spectrum, since they can better distinguish changes over soil areas. In contrast, the NIR (Figure 10), which is usually linked with the temporal changes of vegetation cannot provide any meaningful information for the specific case study.

## 5. Conclusion

While monitoring of looting cannot prevent illegal activities, the use of Earth Observation sensors has attracted the interest of scholars to explore new ways for monitoring archaeological sites, especially in conflict zones. In the recent past, researchers were focused on the exploitation of high-resolution optical satellite datasets, through various processing analysis, showing promising results. However, as it was found, the use of medium resolution images has been limitedly discussed due to the constraints of the spatial resolution of these sensors.

In this study, we focused on the application of various processing analysis techniques at a multi-temporal medium resolution dataset of Landsat 7 ETM+ images over the site of Apamea in Syria, during the period 2011–2012. It was shown that the looted area was able to be detected, in the pseudo-colour composites, PCA and classification analysis. Despite some false positives results, we could further narrow-down specific periods for future elaboration using high-resolution datasets. In specific, an interesting period was found between 2011-10-11 to 2011-11-28 (image no. 10 to no. 12 of Table 1). During this period, we reported spectral changes due to the looting activity.

To avoid confusion with seasonal spectral changes due to the phenological cycle of vegetation, or land use/land cover changes of the area, we have also examined other areas in the surrounding area of Apamea. By monitoring other areas of interest – beyond the Apamea archaeological site per se -, we can minimize errors from false positives. For instance, in Figure 9, during the beginning period of our monitoring, we can observe a high increase of the blue band reflectance value over the looted and soil areas. In this case, both these regions seem to have approximately an increase of 40%. This change can be however linked with the increase/decrease of the moisture of the soil (e.g. after rainfall) and not directly to the looting activity.

From the findings of this study, two critical aspects can be highlighted. Firstly, the use of medium resolution images – as those of Landsat – was only feasible due to the extent of the looted area, so large that it could be detected in a 30m pixel resolution; secondly, the exploitation of temporal data cube was feasible due to the systematic observation of the spaceborne sensors’ data which can be nowadays found in both NASA/USGS (Landsat) and the European Copernicus (Sentinel) space programs: an increased spatial resolution dataset – compared to the 30m resolution Landsat images – can be found through the exploitation of the Sentinel 2A and 2B sensors providing 10m and 20m spatial resolution images. These sensors with a high revisit time (5 days) can act as a systematic observatory tool for local stakeholders and other interested parties. A time-stamp approach based on the analysis of pairs of images can spot areas of ‘change’ and scholars can then proceed with further high-resolution or ground investigations.

While the spatial resolution of open-access datasets like the Landsat is not expected to change in the forthcoming years, the exploitation of multi-temporal datasets has recently been available (2018) in high-resolution commercial sensors like the Dove Satellite Constellation, which combines a fleet of nanosatellites acting as a satellite constellation. These technological trends are expected to impact future space-based remotely sensed applications as the case of the systematic monitoring of archaeological landscapes and sites.

## Acknowledgements

The results are part of the project ‘Synergistic Use of Optical and Radar data for cultural heritage applications’, (PLACES), under the Research and Innovation Foundation grant agreement CULTURE/AWARD-YR/0418/0007 funded by the Republic of Cyprus.

## Competing Interests

The author has no competing interests to declare.

## References

1. Abate, N and Lasaponara, R. 2019. Preventive Archaeology Based on Open Remote Sensing Data and Tools: The Cases of Sant’Arsenio (SA) and Foggia (FG), Italy. Sustainability, 11(15): 4145. DOI: https://doi.org/10.3390/su11154145

2. Agapiou, A. 2017. Remote sensing heritage in a petabyte-scale: satellite data and heritage Earth Engine© applications. International Journal of Digital Earth, 10(1): 85–102. DOI: https://doi.org/10.1080/17538947.2016.1250829

3. Agapiou, A, Alexakis, DD and Hadjimitsis, DG. 2019. Potential of Virtual Earth Observation Constellations in Archaeological Research. Sensors, 19: 4066. DOI: https://doi.org/10.3390/s19194066

4. Agapiou, A and Lysandrou, V. 2015. Remote sensing archaeology: Tracking and mapping evolution in European scientific literature from 1999 to 2015. Journal of Archaeological Science: Reports, 4: 192–200. DOI: https://doi.org/10.1016/j.jasrep.2015.09.010

5. Agapiou, A, Lysandrou, V and Hadjimitsis, DG. 2017. Optical Remote Sensing Potentials for Looting Detection. Geosciences, 7: 98. DOI: https://doi.org/10.3390/geosciences7040098

6. Balty, J. 1969. Apamée de Syrie. Bilan des recherches archéologiques 1965–1968. Actes du colloque tenu à Bruxelles les 29 et 30 avril 1969. In: Balty, J, Dulière, C and Theunissen, M (eds.), Colloque Apamée de Syrie (1st: 1969: Brussels, Belgium). Bruxelles: Centre belge de recherches archéologiques à Apamée de Syrie. pp. 147.

7. Campbell, JB. 2007. Introduction to Remote Sensing. New York: The Guilford Press.

8. Casana, J. 2015. Satellite Imagery-Based Analysis of Archaeological Looting in Syria. Near Eastern Archaeology, 78: 142–152. DOI: https://doi.org/10.5615/neareastarch.78.3.0142

9. Casana, J and Panahipour, M. 2014. Satellite-Based Monitoring of Looting and Damage to Archaeological Sites in Syria. Journal of Eastern Mediterranean Archaeology and Heritage Studies, 2: 128–151. DOI: https://doi.org/10.5325/jeasmedarcherstu.2.2.0128

10. Caspari, G. 2018. Assessing Looting from Space: The Destruction of Early Iron Age Burials in Northern Xinjiang. Heritage, 1: 320–327. DOI: https://doi.org/10.3390/heritage1020021

11. Chen, F, Lasaponara, R and Masini, N. 2017. An overview of satellite synthetic aperture radar remote sensing in archaeology: From site detection to monitoring. Journal of Cultural Heritage, 23: 5–11. DOI: https://doi.org/10.1016/j.culher.2015.05.003

12. Cigna, F, Lasaponara, R, Masini, N, Milillo, P and Tapete, D. 2014. Persistent Scatterer Interferometry Processing of COSMO-SkyMed StripMap HIMAGE Time Series to Depict Deformation of the Historic Centre of Rome, Italy. Remote Sensing, 6: 12593–12618. DOI: https://doi.org/10.3390/rs61212593

13. Danti, M, Branting, S and Penacho, S. 2017. The American Schools of Oriental Research Cultural Heritage Initiatives: Monitoring Cultural Heritage in Syria and Northern Iraq by Geospatial Imagery. Geosciences, 7: 95. DOI: https://doi.org/10.3390/geosciences7040095

14. Finlayson, C. 2012. New Excavations and a Reexamination of the Great Roman Theater at Apamea, Syria, Seasons 1–3 (2008–2010). American Journal of Archaeology, 116(2): 277–319. DOI: https://doi.org/10.3764/aja.116.2.0277

15. Google Earth. Available at https://www.google.com/earth/ [Last accessed 09 February 2020].

16. Landsat Collections. Available at https://www.usgs.gov/land-resources/nli/landsat/landsat-collection-1?qt-science_support_page_related_con=1#qt-science_support_page_related_con [Last accessed 09 February 2020].

17. Lauricella, A, Cannon, J, Branting, S and Hammer, E. 2017. Semi-automated detection of looting in Afghanistan using multispectral imagery and principal component analysis. Antiquity, 91: 1344–1355. DOI: https://doi.org/10.15184/aqy.2017.90

18. Li, J and Roy, DP. 2017. A Global Analysis of Sentinel-2A, Sentinel-2B and Landsat-8 Data Revisit Intervals and Implications for Terrestrial Monitoring. Remote Sensing. 2017, 9: 902. DOI: https://doi.org/10.3390/rs9090902

19. Luo, L, Wang, X, Guo, H, Lasaponara, R, Zong, X, Masini, N, Wang, G, Shi, P, Khatteli, H, Chen, F, Tariq, S, Shao, J, Bachagha, N, Yang, R and Yao, Y. 2019. Airborne and spaceborne remote sensing for archaeological and cultural heritage applications: A review of the century (1907–2017). Remote Sensing of Environment, 232: 111280. DOI: https://doi.org/10.1016/j.rse.2019.111280

20. NASA. 2020. POWER Data Access Viewer, Prediction Of Worldwide Energy Resources. Available at https://power.larc.nasa.gov/data-access-viewer/ [Last accessed 09 June 2020].

21. Orengo, HA and Petrie, CA. 2017. Large-Scale, Multi-Temporal Remote Sensing of Palaeo-River Networks: A Case Study from Northwest India and its Implications for the Indus Civilisation. Remote Sensing, 9: 735. DOI: https://doi.org/10.3390/rs9070735

22. Parcak, S. 2017. Moving from Space-Based to Ground-Based Solutions in Remote Sensing for Archaeological Heritage: A Case Study from Egypt. Remote Sensing, 9: 1297. DOI: https://doi.org/10.3390/rs9121297

23. Parcak, S, Gathings, D, Childs, C, Mumford, G and Cline, E. 2016. Satellite evidence of archaeological site looting in Egypt: 2002–2013. Antiquity, 90(349): 188–205. DOI: https://doi.org/10.15184/aqy.2016.1

24. Stewart, C, Oren, ED and Cohen-Sasson, E. 2018. Satellite Remote Sensing Analysis of the Qasrawet Archaeological Site in North Sinai. Remote Sensing, 10: 1090. DOI: https://doi.org/10.3390/rs10071090

25. Tapete, D, Casagli, N, Luzi, G, Fanti, R, Gigli, G and Leva, D. 2013. Integrating radar and laser-based remote sensing techniques for monitoring structural deformation of archaeological monuments. Journal of Archaeological Science, 40: 176–189. DOI: https://doi.org/10.1016/j.jas.2012.07.024

26. Tapete, D and Cigna, F. 2017. Trends and perspectives of spaceborne SAR remote sensing for archaeological landscape and cultural heritage applications. Journal of Archaeological Science: Reports, 14: 716–726. DOI: https://doi.org/10.1016/j.jasrep.2016.07.017

27. Tapete, D and Cigna, F. 2018. Appraisal of Opportunities and Perspectives for the Systematic Condition Assessment of Heritage Sites with Copernicus Sentinel-2 High-Resolution Multispectral Imagery. Remote Sensing, 10: 561. DOI: https://doi.org/10.3390/rs10040561

28. Tapete, D and Cigna, F. 2019. Detection of Archaeological Looting from Space: Methods, Achievements and Challenges. Remote Sensing, 11: 2389. DOI: https://doi.org/10.3390/rs11202389

29. Tapete, D, Cigna, F and Donoghue, DNM. 2016. “Looting marks” in spaceborne SAR imagery: Measuring rates of archaeological looting in Apamea (Syria) with TerraSAR-X Staring Spotlight. Remote Sensing of Environment, 178: 42–58. DOI: https://doi.org/10.1016/j.rse.2016.02.055

30. Tzouvaras, M, Kouhartsiouk, D, Agapiou, A, Danezis, C and Hadjimitsis, DG. 2019. The Use of Sentinel-1 Synthetic Aperture Radar (SAR) Images and Open-Source Software for Cultural Heritage: An Example from Paphos Area in Cyprus for Mapping Landscape Changes after a 5.6 Magnitude Earthquake. Remote Sensing, 11: 1766. DOI: https://doi.org/10.3390/rs11151766

31. UNESCO, Observatory of Syrian Cultural Heritage. 2010. Available at https://en.unesco.org/syrian-observatory/news/apamea-afamia [Last accessed 09 February 2020].