‘ Looting marks ’ in space-borne SAR imagery: Measuring rates of archaeological looting in Apamea (Syria) with TerraSAR-X Staring Spotlight Remote Sensing of Environment

In archaeological remote sensing, space-borne Synthetic Aperture Radar (SAR) has not been used so far to mon- itor ‘ looting ’ (i.e.illegalexcavationsinheritagesites)mainlybecauseofthespatialresolution ofSARimages,typ-ically not comparable to the ground dimensions of looting features. This paper explores the potential of the new TerraSAR-X beam mode Staring Spotlight (ST) to investigate looting within a work ﬂ ow of radar backscattering change detection. A bespoke time series of ﬁ ve single polarisation, ascending mode, ST scenes with an unprece- dented azimuth resolution of 0.24m was acquiredoverthe archaeological site of Apamea in westernSyria,from October 2014 to June 2015 with a regular sampling of one image every two months. Formerly included in the Tentative List of UNESCO, the site has been heavily looted from at least early 2012 to May 2014, as con ﬁ rmed by Google Earth Very High Resolution (VHR) optical imagery. Building upon the theory of SAR imaging, we de- velop a novel conceptual model of ‘ looting marks ’ , identify marks due to occurrence of new looting and discrim-inate them from alteration (e.g. ﬁ lling) of pre-existing looting holes. ‘ Looting marks ’ appear as distinctive patterns of shadow and layover which are visible in the ground-range reprojected ST image and generated by the morphology of the holes. The recognition of looting marks within ratio maps of radar backscatter ( σ 0 ) be- tween consecutive ST scenes allows quanti ﬁ cation of the magnitude, spatial distribution and rates of looting ac-tivities. Inagreement with the estimatesbased on Google Earth imagery, the STacquiredin October 2014 shows that ~45% of the site was looted. In the following eight months new looting happened locally, with holes mainly dugalongthemarginsofthealreadylootedareas.Texturevaluesof~0.31clearlydistinguishtheseholesfromthe unaltered,baregroundnearby.Hotspotsofchangeareidenti ﬁ edbasedonthetemporalvariabilityof σ 0 ,andcol-ourcompositesindicatewhererepeatedlootingandalterationofexistingholesoccurred.Mostlootingmarksare observednorthofthetwomainRomandecumani.Lootingintensi ﬁ edalmoststeadilyfromDecember2014,with over 1500 new marks in February – April 2015. The estimated rates of looting increased from 214 looting marks/ month in October – December 2014 to over 780 marks/month in April – June 2015, and numerically express the dynamic nature of the phenomenon to which Apamea is still exposed. The method of identifying looting marksinVHRradarimagesthereforeprovesareliableopportunityforarchaeologistsandimageanalyststomea-sure remotely the scale of looting and monitor its temporal evolution.


Introduction
In archaeological heritage science, the term 'looting' refers to unauthorised excavations without any scientific purpose that aim to remove goods of historical or cultural value, frequently used to feed the clandestine market trafficking antiquities. Although the true scale of this phenomenon has been a matter for conjecture (Brodie, Doole, & Renfrew, 2001), there is no doubt across the international community that looting is a plague across the world. Looting even affects some western countries with long standing tradition of archaeological conservation (Proulx, 2013).
Since the 1950s many efforts to encourage protection of World Heritage were done by the United Nations Educational, Scientific and Cultural Organization (UNESCO) via recommendations, conventions and promotion of international cooperation agreements (UNESCO, 1956(UNESCO, , 1970. Nevertheless, the last decades saw looting spread especially in remote areas of developing countries, with limited or absent surveillance, also taking advantage of situations such as conflicts when authoritative efforts are focused elsewhere. In such circumstances, looting tends to worsen, if not even to accelerate. Actions to protect the local Remote Sensing of Environment 178 (2016) [42][43][44][45][46][47][48][49][50][51][52][53][54][55][56][57][58] heritage are then further constrained by the inaccessibility of the site due to security considerations. The methods (e.g., hand-digging or using machines) and rates (e.g., systematic and steady or occasional and sudden) with which looting occurs contribute to determine the degree of destruction of the archaeological context. It is to be acknowledged that looting also causes irreversible damage to the landscape and the anthropogenic environment.
Building upon a well-established tradition of aerial photography and discrimination of spectral signatures, the above mentioned research exploited optical imagery exclusively. In particular, image analysts aimed to identify the circular holes left by looters, in most cases undertaking visual estimates and manual zoning (e.g., Cunliffe, 2014). In other studies, semi-automated data processing was developed based on spatial autocorrelation (Lasaponara et al., 2012(Lasaponara et al., , 2014. On the contrary, there is an extreme paucity of studies exploring the capabilities of Synthetic Aperture Radar (SAR) imagery for this specific type of application. The well-known advantage of SAR sensors to image, in theory, under any weather conditions, should overcome the frequently claimed limitations of optical imagery being occasionally hampered from clear visibility due to clouds (UNITAR, 2014) and the variability of illumination from scene to scene (Stone, 2008). An exception is the research published by , which demonstrated that, in agreement with other satellitebased observations, even SAR time series of medium spatial resolution (~25-30 m) enabled the detection of changes in areas with no record of authorised archaeological excavations. But so far this study stands alone.
Practitioners long argued that SAR was not of sufficient resolution for successful implementation in this branch of archaeological remote sensing, for which sub-metre level of detail is required. Nonetheless, in the current context of SAR mission development, the new Staring Spotlight Mode (ST), which has been recently released by the German Aerospace Center (DLR), opens interesting perspectives. This mode enhances the TerraSAR-X azimuth resolution by means of ST imaging in combination with an extended azimuth pattern steering (Mittermayer, Wollstadt, Prats-Iraola, & Scheiber, 2014). With an azimuth resolution of up to 0.24 m over the scene extent varying between 2.5 to 2.8 km in azimuth and 4.6 to 7.5 km in range, ST can be rightly regarded as revolutionary compared to other SAR beam modes as it brings, for the first time, SAR to a resolution level closer to that of Very High Resolution (VHR) optical imagery (e.g., WorldView-2/3 offer~0.3-0.4 m panchromatic resolution). In the radar domain, at present no better resolution can be obtained with the existing space missions. SAR images acquired by COSMO-SkyMed in X-band (i.e. similar to TerraSAR-X) can reach 1 m resolution in Spotlight mode, while the corresponding mode of L-band ALOS PALSAR 2 provides 3 m range by 1 m azimuth resolution.
This paper aims to investigate how to exploit the unprecedented imaging capability of ST images to detect the presence or modification of looting holes or pits, and measure their rates of occurrence across an area of interest.
To this scope, we investigate the archaeological site of Apamea in western Syria. In this site, the occurrence of looting during the Syrian civil war which started in March 2011 is ascertained, and also well documented in official reports (UNITAR, 2014), scientific publications (Casana & Panahipour, 2014;Cunliffe, 2014) and international media (Lawler, 2014) based on optical imagery from commercial satellites.
We first apply principles of radar imaging to a range of geometric parameters (width, depth and orientation) that reflect those characterising looting holes or pits, and simulate how these holes appear in an ST image. We then discuss the results obtained by analysing a bespoke ST time series acquired over Apamea from October 2014 to June 2015, with sampling frequency of one image every two months. In particular, we combine texture extraction, analysis of the radar backscatter and amplitude-based change detection to: (1) identify looting marks, map their spatial distribution and temporal evolution; and (2) estimate rates of archaeological looting.
The intended impact of this paper is twofold. As a technology-led research, it demonstrates the capabilities of ST beam mode for identifying and mapping looting marks in radar imagery, and monitoring looting as a dynamic phenomenon and measuring its rates. It is envisaged that the proposed method will be considered as a new remote sensing opportunity to monitor looting for use in similar contexts and in sites for which other types of assessment are not possible. Our results also update the knowledge about the condition of the archaeological heritage in Apamea which, as of mid-2015, could only be inferred from commercial imagery updated until mid-2014.

Study area
The archaeological site of Apamea is located in western Syria (centre coordinates: 35.41992°N, 36.401220°E),~50 km north-west of the town of Hama, and is situated on the right bank of the Orontes, at the top of a high relief overlooking the Ghab plain (Fig. 1a). Within the site, elevation ranges between 212 and 273 m a.s.l. and topographic gradients are generally very gentle, 4.5°on average (Fig. 1b).
This strategic topographic position was exploited to control the trade routes by the founder Seleucus Nicator in 300-299 BC, who named the site after his Bactrian wife Apama (or Afamia). 7 km-long ramparts refortified by the emperor Justinian during the 6th century AD enclose the site (Fig. 2a), and are its architectural boundary, marking the natural ridge of the relief along the eastern, southern and western sides, where slopes are as steep as 23° (Fig. 1b). As a convention for this paper, we conform to the archaeological literature about Apamea (e.g., Balty, 1969;Vannesse, Haut, Debaste, & Viviers, 2014) and refer to this boundary to indicate the extent of the archaeological site and region of interest of our spatial analysis.
The site covers an area of~2.44 km 2 , of which almost 40% was uncovered since the 1930s by international archaeological missions (Fig. 2a). Before the Syrian civil war, the excavated sectors were mainly concentrated east of the modern road that crosses Apamea from north to south (number 4 in Fig. 2a) and included the Roman theatre in the western head towards Qalaat Al-Madiq (number 5 in Fig. 2a). Privatelyowned cultivated fields laid to the west of the modern road, and agricultural land covered most of the southern portion of the site, thereby resulting in~1.47 km 2 of unexcavated areas (grey polygons in Fig. 2a). Numerous photographs available online (e.g., Rdefrankrijker, 2009) show the characteristics of the landscape prior to looting.
The high archaeological potential of Apamea was acknowledged in 1999 among the elements to support its candidacy for inscription in the World Heritage List (WHL) of UNESCO. The most remarkable features of the site include the monumental columns and porticoes of the Cardo Maximus (number 1 in Fig. 2a). This 20 m-wide monumental road runs across the site from north to south for about 1.8 km and was formerly flanked by some of the most noticeable public places of the ancient town, including the Roman temple of Tycheion, the Agorà and the Roman Market. Luxurious villas such as the House of Consoles and the House of Pilasters (number 6 in Fig. 2a) were located along the Decumanus Maximus (number 2 in Fig. 2a), alongside the socalled Eastern Cathedral built during the 6th century AD (number 7 in Fig. 2a).
After residing in the Tentative List of UNESCO for several years after the WHL submission, unfortunately Apamea can be now considered among the world heritage at risk due to the recent history of looting. Illegal excavations spread across the site during the Syrian conflict, with looting holes extensively pock-marking the landscape. The international heritage community largely agrees that illegal excavations in Apamea were the result of a systematic and organised looting campaign, carried out on a massive scale (Westcotte, 2014), during which looters exploited bulldozers and/or machinery like backhoes or excavators (Casana & Panahipour, 2014;Lawler, 2014; see Section 3). Apart from explaining the rapid progression of looting in a relatively short period of time, this background information also helps to understand the morphology and distribution with which looting holes appear across Apamea.
Based on optical imagery available on Google Earth we can infer that, as of 4th April 2012, 0.93 km 2 (i.e.~38%) of the archaeological site was looted. In particular, looting devastated~75% of the excavated sectors plus~12% of the unexcavated areas. In this regard, Casana & Panahipour (2014) hypothesised a major trend of the archaeological looting by observing that the phenomenon was mostly concentrated east of the modern road (number 4 in Fig. 2a), i.e. in the governmentowned heritage site. More recent satellite images available in Google Earth, such as that acquired on 28th September 2012 by IKONOS, reveal the need to revise their initial hypothesis. Looting holes increasingly started to appear in the private-owned land west of the modern road, from a few clusters of looting holes covering 0.015 km 2 in September 2012 to 0.105 km 2 in March 2014, thereby heavily covering the northwestern sector of the formerly cultivated land (Fig. 2b). As recently pointed out by UNITAR (2014), there is evidence that looting continued across 2013 and still represents a threat to site preservation. Media and published reports provided rough numbers of looting holes (Danti & Prescott, 2014 refer to over 5000 holes) but, to the best of our knowledge, an assessment of the rates at which looting is occurring in Apamea has not as yet been attempted. Fig. 2b displays the map of looting areas that we have drawn directly on the Pléiades satellite image acquired on 6th March 2014 and made available through Google Earth. Looting holes are found across more than 44% of the total extent of the site in March 2014. Of these, nearly 80% were distributed across the former zones excavated by archaeologists, while the remainder were dug in the unexcavated sectors lying west of the modern road. Although a new Pléiades scene acquired on 2nd May 2014 has recently appeared on Google Earth, this covers only 2/3 of the site and thus the March image is still the most up-to-date and complete of the Google Earth Apamea site timeline as of mid-2015. The looting map derived from this image, and reported in Fig.  2b, is therefore used in this study as a temporal reference for the SAR analysis (see Section 5).
Our choice of referring to Google Earth images is justified not only because these data were freely accessible and had been already used by previous studies of looting in Apamea cited above, but also in light of the evidence that no other more updated optical image of comparable resolution was available as of mid-2015 when our study was completed. The DigitalGlobe catalogue unfortunately provides a confirmation, as optical image acquisition over the geographical area including Apamea was discontinuous in 2014 to mid-2015, also with variable spatial resolution and site coverage and, sometimes, extensive cloud coverage. Although this represented a constraint for this research, it did not prevent us to retrieve the baseline scenario of looting in Apamea prior to our monitoring campaign, at a spatial resolution similar to our TerraSAR-X Staring Spotlight data stack (see Sections 4.1 and 5). Anyway, it is worth mentioning that it was beyond the scope of this research to undertake a dual monitoring (optical and radar) or to perform a comparative or contrasting assessment of the two remote sensing technologies.

Looting marks in VHR SAR imagery
As documented by UNITAR (2014) with regard to cultural heritage across the Syrian Arab Republic, looting manifests in different ways, including holes or pits, footpaths dug to reach the entrance of tombs, removal and accumulation of soil, debris from collapse or vandalism of ruins, trenches and illegal excavations. In this work we focus on looting holes and identify how these are imaged by the side-looking viewing geometry of the TerraSAR-X Staring Spotlight (ST) acquisition mode. This accounts for the fact that holes are worldwide among the most common forms of looting and predominate in Apamea.
According to photographic and video documentation captured by the Directorate General of Antiquities and Museums, DGAM (UNESCO, 2015;Westcotte, 2014) and the Association for the Protection of Syrian Archaeology -APSA2011 (2014a, 2014b), the typical morphology of looting holes in Apamea mainly reflects the use of excavators and bulldozers. Most of the observed holes are characterised by square, rectangular or pseudo-circular openings, sometimes L-shaped and, in some cases, perfectly circular or more irregular shapes. The planimetric dimensions of the openings range between tens of centimetres to a few metres. Similarly, excavation depths are from less than one metre to a couple of metres, and even more in exceptional cases. This is a consequence not only of the method used by the looter to dig the soil, but also of an opportunistic rationale by which the looter's perception of the likelihood of finding goods is the driver to progress with excavation. In most cases, however, the depth is not the predominant dimension. Moreover, given the width of the openings, it is reasonable to assume that the holes are not as deep as they are wide. VHR optical imagery depicting the site in 2012-2014 and accessed from Google Earth confirm the planimetric characteristics of the holes. These are characterised by regularly shaped footprints and rectangular and squared openings that are mainly oriented with sides parallel to the north-south and east-west directions, following the orthogonal city plan of the archaeological site. Only very few exceptions from this pattern are observed. Individual holes and looting clusters are often surrounded by deposits of the excavated material, generally extending few tens of centimetres tõ 1-1.5 m around each hole, and the photographic documentation (UNESCO, 2015) suggests that these are up to a few tens of centimetres thick.
Based on this evidence, we now introduce the conceptualisation of how a looting hole is imaged in VHR SAR imagery. Fig. 3a shows a schematic representation of a regularly shaped looting hole, as seen by drawing a cross section along the range direction. The dimensions of the hole are defined by its width l and depth h, with l predominating on h. This conceptualisation is analogous to those developed by several authors working on building reconstruction and width and height retrieval based on shadow and layover signatures in VHR SAR imagery (e.g., Franceschetti, Guida, Iodice, Riccio, & Ruello, 2005;Brunner, Lemoine, Bruzzone, & Greidanus, 2010). In this respect, looting holes are conceptualised as depressions in the digital terrain models rather than elements sticking out from the ground, and their simplified scattering model is similar to that of a road between two buildings, where the latter constitute the walls of the looting hole.
Building upon the methodology implemented for the simulation of SAR distortions by Cigna, Bateson, Jordan, & Dashwood (2014), we modelled the TerraSAR-X satellite sensor location as illumination source of a hill-shading and shadowing model. We defined the azimuth and altitude angles of the radar sensor according to the orientation parameters of the Staring Spotlight Line-Of-Sight (LOS), and used heading angle γ of 350.2°(i.e. −9.8°from the north-south direction) and incidence angle θ of 39.7°at scene centre (see Section 4.1).
We considered only single backscattering from the walls and base of the looting hole. This is believed to predominate with respect to doublebouncing caused by the dihedral corner reflector arising from the walls and base. According to the Rayleigh roughness criterion (Rayleigh, 1945; formerly published in 1877), the excavated surfaces are rough relative to the 3.1 cm radar wavelength and 39.7°incidence angle of the ST imagery (according to which surfaces with features higher than 5 mm are considered rough), and hence produce negligible specular reflection.
The resulting portions of the hole affected by radar shadow and layover are indicated in Fig. 3. Whilst radar layover produces strong radar returns and thus very bright pixels in the radar image (which are represented in cyan in Fig. 3), areas of shadow are those that cannot be illuminated by the radar LOS thus producing no radar returns to the sensor, and are therefore imaged in the SAR scene as darker pixels ( Fig. 5 and Fig. 10). Areas affected by shadow extend l S = h ⋅ tan θ on the ground, from the edge of the hole producing the shadow (i.e. the closest to the sensor) and away from the sensor, whilst layover extends l L ¼ h tanθ from the farthest edge and towards the sensor. This means, for instance, that a~0.83 m wide portion of a 1 m deep hole will be in shadow and a~1.20 m wide portion will be affected by layover. It becomes apparent that the shadow and layover produced by each hole can be either in close proximity, adjacent or overlapped, depending on the relationship between l and h. In particular, these are in close proximity when l S +l L b l, i.e. when the depth h is smaller than lcos θsinθ (i.e.~0.5⋅ l in this case). They are adjacent and therefore interfere when h equals 0.5⋅l, or overlap if h exceeds the latter. Shadow occludes the looting hole when l S ≥l, hence when h reaches or exceeds l tanθ (i.e.~1.2⋅l in this case). When the holes are imaged by the SAR sensor, the backscattered signal is captured in the slant range geometry, as shown in Fig. 3a. The dimensions along the Slant Range (SR) of the shadow and layover produced by the presence of the hole are SR S ¼ h cosθ and SR L =hcosθ respectively. Under the assumption of a flat-Earth model and absence of the hole, the simple reprojection of the looting marks from the slant range onto the Ground Range (GR) geometry generates marks equal to G R S ¼ h sinθ cosθ and GR L ¼ h tanθ for the shadow and layover respectively (Fig. 3b), where the relation between the slant range and ground range is defined by SR = GR sin θ. It is to be noted that, whilst the ground range shadow mark falls within the hole producing it and is (sinθ) −2 times wider than l S (i.e.~2.45 wider in this case), the layover mark extends from the farthest edge and away from the sensor, and its extension equals l L .
To conceptualise the theoretical radar signatures of the holes, analyse an assorted sample of looting hole configurations, and account for the low topographic gradients across the site (see Section 2), we generated a flatterrain Digital Elevation Model (DEM) with a number of looting holes reflecting a range of shapes, widths, depths and orientations. This simulation did not account for the random effects of the radar speckle and was performed at 1 cm spatial resolution, to derive simplified, typical looting marks expected from the different configurations. In particular, we defined the following five sets of looting holes (Fig. 4): • Set A: seven square holes, characterised by fixed depth of 1 m, and increasing width from 0.5 to 3.5 m, in 0.5 m increments; Fig. 3. Cross section of a regularly shaped looting hole along the range direction: (a) radar shadow and layover within the hole resulting from the side-looking geometry of the radar, alongside the respective 'looting mark' in the slant-range direction, i.e. SR S and SR L ; and (b) reprojected 'looting mark' onto the ground range geometry, i.e. GR S and GR L . • Set B: seven square holes, characterised by fixed width of 2 m, and increasing depth from 0.25 to 2 m, in 0.25 m increments; • Set C: six square holes, characterised by fixed width of 2 m and depth of 1 m, and increasing tilt α, from 12.5°to 75°, in 12.5°increments; • Set D: six rectangular holes, with fixed width of 2 by 3.5 m, depth of 1 m, and increasing tilt α, from 0°to 150°, in 30°increments; • Set E: four circular holes, with 1 and 2 m diameter and depths of 0.5 and 1 m, and three L-shaped holes with 1 m depth and various orientation and width.
It is acknowledged that the above set cannot be exhaustive of all possible combinations that can be found in real-world situations. Nonetheless, this selection does not preclude the reader to extend the conceptualisation to predict the radar layover and shadow marks for longer, deeper, differently tilted or shaped holes, following the incremental rationale used in each set we have analysed in this study. Fig. 4a shows the results of the simulation indicating which portions of each looting hole are affected by radar shadow and layover, similarly to that represented in Fig. 3a. The simulation reveals the characteristic size and shape of the areas of each looting hole that are affected by layover and shadow with varying morphology and orientation. Common to all looting holes is the combined presence of shadow in their western portion, generated by the walls facing away from the sensor (i.e. eastward), and layover in their eastern portions, generated by walls facing towards the satellite sensor (i.e. westward). It can also be observed that the extension of shadow and layover areas decreases from the theoretical values l S and l L when the edges of the hole tilt away from the range and azimuth directions.
We define as 'looting marks' the combined patterns of radar shadow and layover produced by looting holes as depicted in the slant-range geometry, i.e. along the LOS (i.e. SR S and SR L in Fig. 3a-b), and then reprojected along the ground range using a flat-Earth model (i.e. GR S and GR L in Fig. 3b). Fig. 4b reproduces the ground range reprojection of looting marks for our five sets of simulated holes, highlighting the characteristic shapes and dimensions of the marks when the imagery is geocoded back from the radar to the map geometry using a flat-Earth model. This figure therefore illustrates how looting marks appear in map geometry, without radar speckle. Random effects of speckle and resolution degradation are clearly expected in the real data with respect to the simulated model. The usefulness of such a conceptual model for the identification and count of looting holes across the site of Apamea is explained in Section 4.2 and discussed in Section 5.
As mentioned at the beginning of this section, looting holes are frequently surrounded by mounds of debris accumulated from the excavation. According to our conceptual model in Fig. 3a, such deposits generate layover only if they are steeper than θ (i.e. 39.7°) and face the sensor, or shadow only if they are steeper than 90°-θ (i.e. 50.3°) and face away from the sensor. In such a case, this signal would add onto the layover and shadow components of the looting mark, respectively, thus making its footprint larger. Nevertheless, this situation is very unlikely, because the loose terrain composing the mounds would collapse at such steep slopes. Furthermore, the photographic documentation published by UNESCO (2015) confirms that our assumption reflects the reality observed in Apamea.

Staring Spotlight data stack
Our stack of TerraSAR-X Staring Spotlight (ST) data consists of the following five scenes acquired bi-monthly between October 2014 and June 2015, with temporal baselines of 55 or 66 days according to the TerraSAR-X nominal repeat cycle of 11 days: These images were made available to the authors at no cost according to an experimental acquisition campaign that was agreed with the German Aerospace Center (DLR) in the framework of the TSX-New-Modes-2013 LAN2377 grant.
The acquisition geometry of the stack was designed to cover the entire archaeological site with a single image frame (Fig. 5a), therefore the scene centre was set at the location 35.4°N, 36.50°E. Ascending orbits with 350.2°heading angle and 39.7°incidence angle at scene centre (beam spot_051) were exploited, with a right-looking antenna orientation. A minimum incidence angle of 39.4°in the near range, and a maximum of 40.0°in the far range were used.
This configuration of acquisition parameters resulted in ground coverage of 3.0 km along the azimuth direction and 5.6 km along the range direction, imaging 16.8 km 2 in total, including the entire archaeological site within the ancient city walls, Tell Al-Madiq and the modern town to the west, and the Apamea reservoir to the south-east of the walls and the three dams (Fig. 5a). Although from the literature it is known that looting also affected the neighbouring Tell Jifar (Casana & Panahipour, 2014), this site could unfortunately not be included within the same image swath.
The five scenes were acquired with HH polarisation, 9.65 GHz frequency (corresponding to 3.1 cm wavelength) and 300 MHz range bandwidth variant. Azimuth pixel spacing of 0.17 m and slant range spacing of 0.45 m were achieved, the latter corresponding with 0.70 m ground range spacing at scene centre. These indicate a spatial resolution of 0.24 m in azimuth and approximately 0.92 m along the ground range direction, as per the ST mode technical specifications (Fritz & Eineder, 2013).
The zooms in Fig. 5b-g emphasise the extraordinary level of detail provided by the ST data, and its resolution similarity with optical imagery available through Google Earth. Accounting for the effect of the ascending acquisition geometry, we use the strong backscatter and associated shadows to identify the standing marble columns of the Cardo Maximus, the portico and the monumental column (Fig. 5c-e). The ST also clearly images the major military garrison, complete with bunkers and artillery emplacements ( Fig. 5f-g), which was installed by the Syrian government where the former tourist cafe stood (Lawler, 2014;UNITAR, 2014).
With regard to the temporal resolution, the ST image acquisition campaign was designed to cover an observation period of suitable duration and to sample looting in Apamea at fixed intervals, thus conforming to the rationale of regular monitoring of the condition of the site. For the case of Apamea, the temporal regularity and technical consistency achieved with the TerraSAR-X mission were considered as an unprecedented opportunity to retrieve reliable and comparable estimates of looting rates, thus overcoming the technical constraint of relying on single images or irregular time series, as happened in previous studies (Casana & Panahipour, 2014;Cunliffe, 2014;UNITAR, 2014).
The analysed ST data stack for Apamea covers over half a year, from late autumn to summer. The local climate can be classified as hot semiarid according to the Köppen-Geiger system (Kottek, Grieser, Beck, Rudolf, & Rubel, 2006), with scarce monthly precipitation, generally not exceeding 70 mm in the winter season and mostly due to few isolated events. Fig. 6 reports the meteorological data made available by the National Oceanic and Atmospheric Administration (NOAA) for the weather station located in Hama (35.117°N, 36.750°E, elevation 303 m a.s.l.). The data cover the period October 2014-April 2015, which is of interest for this research. TerraSAR-X ST images were acquired in dry months, with total monthly precipitation not exceeding 7 mm, therefore no major concerns arise in terms of interference to the radar signal due to local condition of the soil or vegetation (see Section 5). The general absence of dense vegetation is also confirmed by the land cover of the site as seen from Google Earth imagery since 2003 and a plethora of images accessible online.

Image processing for change detection
We processed our ST stack using a SAR amplitude change detection workflow similar to that employed by ; Tapete et al. (2013); and Tapete, Cigna, Donoghue, & Philip (2015).
TerraSAR-X Single Look Slant Range Complex (SSC) Level 1 products including amplitude and phase information were first imported into the GAMMA SAR and Interferometry software and then co-registered to the first acquisition (22nd October 2014) using the cross-correlation method and polynomial offset models in range and azimuth, obtaining subpixel co-registration precisions of less than~0.1 pixel along the azimuth direction (hence smaller than~0.016 m) for the four slave scenes.
Radar intensity was computed in terms of radar backscattering coefficient σ 0 (sigma nought), indicating the radar signal backscattered from the imaged targets to the sensor, and normalised to the horizontal surface of the WGS84 ellipsoid (Fig. 5a, b, c, f). This coefficient depends on the properties of both the imaged surface (i.e. dielectric constant, roughness and local incidence angle), and the radar signal (i.e. wavelength and polarisation). To reduce the radar speckle, Multi-Look Intensity (MLI) images in slant range coordinates were also generated using a multi-look factor of 5 in the azimuth direction, thus bringing the azimuth pixel spacing to 0.83 m.
Each looting hole generates a combined pattern of radar shadow and layover (namely 'looting mark'; Section 3). Following the conceptual model (Fig. 4b), in the radar image this appears as a pair of abrupt spatial changes in the backscatter σ 0 compared to more uniform radar signatures of the surrounding un-looted ground. Our conceptual model has consequently the added value of explaining not only how a looting hole is imaged by the radar sensor, but also its effect on local texture. By extracting the textural properties of the scene, we could therefore retrieve an improved site-scale understanding of the extent of looted areas.
Texture images are generally computed using first and second order statistics such as the coefficient of variation, contrast, inverse moment and uniformity (e.g., Kurvonen & Hallikainen, 1999;Ulaby, Kouyate, Brisco, & Williams, 1986;Paudyal, Eiumnoh, & Aschbacher, 1995). These are used to replicate the capability of the human eye to delineate common spatial patterns, easing the identification of the boundaries of objects, surface features and, more generally, radar backscatter discontinuities across the scene.
In this study, texture values at each location i were derived by employing a moving kernel of m by m pixels centred at pixel i, and by computing the difference between the logarithm of the average radar backscatter and the average of the logarithms of the backscatter of the m 2 pixels j within the kernel (GAMMA RS, 2014): To implement the above, we exploited the multi-looked SAR amplitude images with σ 0 in the linear scale, and used kernels of 9 by 9 pixels, meaning that the size of the filtering area on the ground was of~8.3 m, as 0.92 m was the nominal ground resolution. A Gaussian weighting was also applied within the kernels, with weighting coefficients w j decreasing from 1.0 to 0.02 with increasing distance from the centre i to the margins of the window. This was performed to reduce further the effect of radar speckle, to enhance over 50 times the value of the backscattering at the centre of the kernel with respect to that at the margins, and to highlight edges between spatially variable radar signatures ( Fig.  7 and Fig. 8). This allowed us to detect the sharp spatial variations in backscatter occurring in the presence of looting holes and derive a generalised delineation of looted vs. non-excavated zones of the site in order to estimate the spatial extent of looting as at the end of 2014 with respect to the scenario provided by the currently available, most up-to-date and complete Google Earth image for March 2014.
Multi-temporal RGB (Red-Green-Blue) and RC (Red-Cyan) colour composition with groups of three scenes or two respectively was used to identify the location and extent of sectors that underwent significant radar backscatter changes during the eight monitored months. In the derived RGB and RC products (Fig. 9a), areas in white to dark grey indicate unchanged backscatter properties, whereas the tints of coloured pixels reveal the scene dates that recorded temporal changes in the backscatter properties.
To analyse the temporal variability of the radar backscatter across the site between October 2014 and June 2015, we computed the mean σ 0 ðiÞ  and standard deviation Stdv(i) of the five multi-looked images as follows: where n is the number of averaged scenes (i.e. t 1 to t n ), σ t 0 (i) is the radar backscatter of pixel i at time t, and σ 0 values are all expressed using the dB scale.
Whilst the computation of the mean backscatter of the 5 images allowed us to enhance the radar signature of the different surfaces, their standard deviation with respect to the mean revealed sectors of the site that underwent backscatter changes over the eight monitored months and allowed their quantification. The map of temporal variability (Fig. 9b) preserved the initial pixel spacing of the multi-looked scenes, hence 0.45 m and 0.83 m in slant range and azimuth respectively, corresponding to~0.92 m ground resolution pixels.
To quantify the backscatter changes that occurred across the site with bi-monthly frequency, we computed ratios between consecutive image pairs formed based on the ST stack, and between the first and last acquisition. These allowed us to enhance morphological changes that occurred at the scale of the image pixels in the multi-looked scenes (~0.92 m), alongside soil moisture variations between wetter (i.e. October 2014 and June 2015) and drier (i.e. December 2014, February 2015 and April 2015) acquisition days (see Section 4.1). On the other hand, the effects of local topography on the radar backscatter that were common to each pair were compensated, and the variations due to local incidence angles cancelled out (e.g., Nico, Pappalepore, Pasquariello, Refice, & Samarelli, 2000;Scheuchl, Ullmann, & Koudogbo, 2009;Boldt & Schulz, 2012), due to the homogeneity in the acquisition parameters and geometry of the scenes. Prior to ratioing, we spatially filtered the multi-looked scenes using a 3 by 3 pixel window and Gaussian weighting. For each pair of scenes, k and z, the ratioing operator was implemented by using the formula: where t k and t z are the acquisition times of scenes k and z respectively, and σ 0 values are expressed using the linear scale. R tz/tk (i) takes on values between 0 and 1 for pixels brighter at time t k with respect to t z , and greater than 1 for pixels brighter at t z with respect to t k . To account for data skewness, the resulting ratios were finally expressed using the dB scale, so negative values indicate pixels where the backscatter was greater at t k than t z , whereas positive values indicate pixels with greater backscatter at t z rather than t k . Co-registered MLIs and all derived products were finally geocoded and ellipsoid-corrected to the WGS84 datum by exploiting the orbital information of the master image, as provided by DLR with nominal pixel localisation accuracy below 20 cm ('SCIENCE' orbit types; Fritz & Eineder, 2013). We refined the geocoding accuracy of ellipsoidcorrected scenes and products using a set of more than 30 GCPs to obtain the best overlap with the available optical reference. We used Virtual Earth and WorldView images accessed via ArcGIS Online that depicted the situation prior to the civil war started in 2011, alongside Google Earth optical imagery covering the years 2003, 2004, 2007 and 2011, the latter in July where no evidence of looting is observed.
It is worth noting that we did not implement any terrain correction during the geocoding of the ST scenes to the map geometry. This choice was driven, primarily, by the very gentle topography of the archaeological site (see Section 2; Fig. 1). No major topographic distortions such as layover or shadow are present across the study area, except for those induced by the presence of the looting holes that are the focus of this study. Moreover, a DEM with at least sub-metre resolution and accuracy (e.g., a LiDAR DTM) would have been necessary to terrain-correct the ST scenes without compromising their information content. To the best of our knowledge, no models with such a resolution were available for the study area. Nevertheless, even if suitable DEMs were available, they would not be helpful to detect looting marks. Indeed, terrain correction of the radar data with an accurate, VHR and contemporaneous DEM would cancel out the topographic distortions produced by the morphology and micro-topography due to the presence of the looting holes. In other words, this operation would remove the radar shadow and layover marks that are detected in this study and that are used to retrieve evidence of looting and estimate the rates of this phenomenon. Moreover, if such a model was provided for future work, the use of elevation data non-contemporaneous with the ST imagery would be detrimental for such an analysis. As proved in Section 5, looting is a dynamic phenomenon which causes a rapid evolution of local micro-topography of the study area. In a time-series approach this would mean repeated acquisitions of VHR DEMs, for instance with drones or airborne sensors, which is not cost-effective and is prohibited for security reasons in areas of conflict.
For these reasons, although terrain correction of the radar scenes is generally recommended for other applications, technical, scientific and practical aspects made a simple ellipsoid correction and flat terrain assumption the ideal and most appropriate processing choices to follow.
As a last note, from a practical point of view, the methodology of change detection described above presents the whole workflow from SSC Level 1data to derived products. This is the best option as it allows every step to be supervised by the operator from processing to postprocessing. Nevertheless, this workflow is flexible to apply to other SAR product formats, such as those referred to as Multi Look Ground Range Detected (MGD), Geocoded Ellipsoid Corrected (GEC) or Enhanced Ellipsoid Corrected (EEC), as per the TerraSAR-X technical specifications (Fritz & Eineder, 2013). These three SAR products are already pre-processed with multi-looking and speckle reduction operations, projection and re-sampling to the WGS84 reference ellipsoid and terrain correction. These may therefore be a preferred option for those operators who are interested in an image format easy to handle in GIS software, for them to follow our method of texture extraction, RGB and RC colour composition and looting mark mapping.

Results and discussion
To the best of our knowledge, as of mid-2015, the Pléiades images acquired on 6th March and 2nd May 2014 were the most recent VHR optical satellite scenes covering Apamea that were available online through Google Earth or reported in the scientific and grey literature ( Fig. 2b; see also Section 2). Therefore our ST scene collected on 22nd October 2014 is the first sub-metre resolution satellite image acquired since then, providing new evidence about the condition of the archaeological site.
In terms of looting areal extent, the map retrieved from the ST acquisition in October 2014 (Fig. 7a) is consistent with the situation recorded in Google Earth (see looting polygons in Fig. 2b). Looting within the ancient walls extended up to~45% of the archaeological site, thereby showing an increase of~1% with respect to that observed in March 2014. Of this,~80% is found over the previously excavated archaeological areas, while the remainder falls over the agricultural fields that were formerly unexcavated by archaeologists.
At a first glance, the radar backscatter (Fig. 5b) shows differential distribution of the radar signatures across the site. This reflects the different surface roughness and local microtopography that distinguish looted areas from smoother bare soil of the un-looted agricultural fields in the western and southern sectors of the archaeological site. This difference is further enhanced by extracting the textural properties of the ST scenes. Fig. 7a displays the texture map obtained from the ST image acquired on 22nd October 2014, by applying a kernel of 9 by 9 pixels and Gaussian weighting (see Section 4.2).
At site scale, the texture map provides an immediate picture of the extent of looting, due to a noticeable contrast with the areas in Apamea that are still untouched by illegal excavations. The latter are indeed characterised by spatial similarity in the backscatter values of adjacent pixels, and result in much lower texture values than looted surfaces. This proves one of the advantages of working with SAR, when compared with optical imagery. In the optical domain, feature detection or classification can be constrained by insufficient colour contrast or unfavourable or inconsistent light conditions, whereas in SAR imaging the alteration of surface roughness caused by looting is used as a morphological marker.
In the texture map, each looted area is delineated by square, rectangular and/or circular features. These have main dimension or diameter in the range of 5-7 pixels, corresponding with 4.5-6.5 m on the ground. The structures of the monumental colonnade are also well delineated. Strips with texture values of~1.3 on average and length of 730, 310 and 220 m (Fig. 7a) coincide with the three main sections of the porticoes (Fig. 5c-e). Similarly, the western walls are clearly demarcated, with texture values of~1.4 on average and always exceeding 0.85. The eastern walls are less demarcated due to the different visibility to the satellite LOS and their texture values can be as low as 0.4 and~0.7 on average.
Throughout the monitoring period, texture values are always higher for looted areas, where small-scale variations in the backscatter occur due to the local difference in the morphology of excavated areas and deposits, and their typical radar marks of shadow and layover, with respect to un-looted zones. The latter are characterised by mean texture values of 0.15 ± 0.05, while nearly double values are observed for looted areas, thereby resulting in a clear separation between these two classes. In particular, mean values of 0.31 ± 0.16 are found both where looting marks are particularly dense and in sectors which were unexcavated prior to the beginning of the Syrian civil war. On the other hand, we observe 0.27 ± 0.17 texture values where looting holes overlapped onto excavated archaeological features (e.g., the House of the Pilasters; number 7 in Fig. 2a). Although the separation between these two sub-classes of looted areas is not sharp, the discrepancy between the values observed can be attributed to the different complexity of the looting textures. Conversely, smoother, non-excavated fields in the western and southern sectors of the site are not highly textured (~0.12) due to the above mentioned spatial similarity in the backscatter values of adjacent pixels. Fig. 7b shows the result of the pixel reclassification based on the numerical analysis of the observed texture values across 18 sample areas, 5 of which are un-looted and 13 densely looted (Fig. 7a). This map proves valuable at both large and small scales. More recent looting features show texture values up to about 0.8 to 1.3, and are concentrated along the margins of the already looted areas. Conversely older looting marks in the southern sector of the archaeological site are characterised by texture values ranging between 0.31 and 0.50, and seem to reflect a progressive degradation of the looting holes, as also suggested by the VHR optical imagery available through Google Earth.
In this regard, Fig. 8 shows the texture variation induced by new looting (i.e. looting occurring in areas not previously excavated or looted) and how texture and its reclassification can be effectively used to map areas of new excavations, such as the new looting marks appearing in December 2014 and February 2015. The example provided reflects the trend observed across the site throughout the monitoring period. From October 2014 to June 2015 new looting happened at the local scale, frequently in the form of individual to clusters of holes. Given the scale of such a phenomenon, some of the new looting might not have been detected so well by using other SAR beam modes, such as TerraSAR-X StripMap mode at 3 m resolution.
The second main evidence we retrieved from the analysis of the TerraSAR-X ST time series is that, across the site, significant changes of radar backscattering and its texture were also found in areas where looting had already occurred. Colour composition of any combination of the five ST backscattering images highlights clearly the occurrence of this type of alteration at all temporal scales (Fig. 9a-e). These colour composite maps complement the textural information, especially where repeated looting and reworking of previous looting could have generated a much more complex morphology and mixture of texture.
In this regard, the map of temporal variability provides a quantitative assessment of how much these areas changed throughout the whole monitoring period. As shown in Fig. 9f-j, these patterns of alteration are distinctly marked as hot spots compared with the surrounding pixels. Values of temporal variability are typically N5 dB and generally up to 10 dB in the eastern and northern sector of the archaeological site. Some variability in the radar backscattering is also observed in the southern and western sectors, where some differences in the soil moisture conditions between the wetter (October 2014 and June 2015) and drier (December 2104, February and April 2015) months are detected across the agricultural fields.
Repeated looting and reworking of areas previously looted are a phenomenon that was already observed in the Google Earth imagery from March to May 2014, for instance in the northern-western part of the site and west of the upper section of the colonnade (Fig. 5c-e). In these areas, colour composites and temporal variability maps highlight a sequence of changes and high temporal variability ( Fig. 9c and h). This is a further proof-of-evidence of the continuity and agreement between what was already known until May 2014 based on Google Earth, and what our analysis using ST scenes provides to characterise and map the varied forms with which looting occurs within the archaeological site.
In this context, we exploit the radar backscattering ratios (see Section 4.2) and the conceptual model of 'looting marks' (see section 3) to record and map both new looting and alteration to existing looted areas. Recalling from Section 3, the typical pattern of 'looting mark' consists of an area of radar shadow coupled with layover, with a shape that follows the morphology of the looting hole as imaged in the radar scene along the slant range and then re-projected onto the ground range using a flat-Earth model. By analogy, the appearance, disappearance or alteration of this pattern in one of the radar scenes is recognisable distinctively in the derived ratios R tz/tk (i), based on its dissimilarity with respect to values for unaltered ground.
To this aim, Fig. 10 provides the interpretation keys of looting marks, by comparing appearance and disappearance of looting holes against circumstances where no alteration is recorded in the radar backscattering and their derived ratio maps. This matrix of interpretation keys accounts for the ascending orbit geometry and acquisition parameters of the TerraSAR-X ST images used in this research (see Section 4.1). A similar matrix can be generated for a descending geometry, but in that case the patterns will be reversed.
Bare ground not affected by looting or other types of morphological alteration does not generate any clear pattern in either the input radar scenes or their ratios (Fig. 10a), and R tz/tk (i) takes on values of~0 dB as σ 0 at time t z almost equals that at time t k . As extensively discussed by Cigna et al. (2013) with regard to SAR change detection in semiarid environments, different R tz/tk (i) values could reflect variations of soil properties such as an increase or decrease of moisture content. Nonetheless, this type of change does not cause any morphological alteration and therefore the combined pattern of shadow and layover is not generated as it would be in case of excavation due to looting.
In this regard, Fig. 10b provides an example of an evident looting mark (mark Type 1) as a result of new looting where the ground was previously untouched. On the other hand, Fig. 10c illustrates how a looting mark appears in the ratio map if pre-existing looting hole has been modified morphologically to the extent of having being filled in from time t k to time t z . The effect of such reworking is a 'reverse looting mark' or 'filling mark' with an area of increased radar backscattering coupled with a decrease, due to the disappearance of the hole and its shadow and layover mark (mark Type 2).
On the contrary, if a looting hole has not been modified and kept the same shape from t k to t z , no pattern of looting mark is detected in the ratio map (Fig. 10d). This acts as a further proof that the patterns as per Fig. 10b-c are reliable markers of looting activities. In particular, the unchanged looting hole reported in Fig. 10d refers to a hole that was already visible in May 2014 from the Google Earth imagery. This example thus demonstrates how looting assessment based on the TerraSAR-X ST time series applies not only to looting occurred after the first SAR acquisition, but also retrospectively to check whether previous looting holes were altered or their condition is unchanged. Fig. 10e illustrates the occurrence of multiple looting marks as a result of systematic digging of an aligned set or cluster of looting holes. Shape, orientation and distribution of the looting marks reflect the type of looting that occurred from the first to the second image of the pair. Comparison of the ratio map with the corresponding radar backscattering images provides confirmation of the causes for the observed marks. In this specific example, it is apparent that the repeated pattern of looting holes reflects systematic excavation activity.
By implementing these interpretation keys, the occurrence of looting marks were mapped across Apamea for each of the monitoring intervals between consecutive TerraSAR-X ST images i.  Fig. 11 compares the resulting four maps of looting marks that appeared or disappeared during each monitoring interval. These maps allow the spatio-temporal analysis of the looting dynamics across the archaeological site, discriminating between the two types of marks (i.e. 'looting mark' or Type 1, and 'filling mark' or Type 2). It is worth noting that each map provides the number of marks appearing or disappearing during the respective monitoring interval, and not the total number of marks (and therefore holes or pits) present across the site.
The clearest evidence provided by these maps is that, during all monitoring intervals, new looting and filling marks are mostly found in the portion of the archaeological site lying north of the Decumanus Maximus (Fig. 11a, area B). More precisely, the majority of the marks observed during the four monitoring intervals falls in the upper part of the Roman city, north of the second main decumanus (Fig. 11a, area C).
In terms of dynamic evolution within the Justinian walls, 470 occurrences of new looting and filling marks are mapped in the interval October 2014-December 2014 (Fig. 11a). These are scattered north of the second main decumanus and also pock-mark the north-western sector of the site. The latter was formerly outside the official boundaries of the excavated archaeological site, which was looted between April 2012 and March 2014 (see Section 2).
From December 2014 to February 2015 we observe intensification in the appearance and disappearance of marks across the site, with total number of occurrences of Types 1 and 2 equal to 1200. Their spatial distribution also suggests a migration of the activity towards south-east (Fig.  11b). Further intensification of new looting and filling marks is found in February 2015-April 2015, with more than 1500 occurrences distributed across the whole upper portion north of the second main decumanus (Fig. 11c). A similar amount of occurrences and coverage are retrieved in the last interval April 2015-June 2015, but in this case the total number of new marks is~1430 and these are distributed more sparsely, with some of them scattered in the area between the two decumani (Fig. 11d).
The above figures allow us to estimate the observed rates of new looting and filling marks occurring during each interval, and how these compare with the preceding months: • October 2014-December 2014: 132 new looting marks/month (Type 1) and 82 filling marks/month (Type 2), hence a total rate of 214 marks/month; • December 2014-February 2015: 420 new looting marks/month (Type 1) and 235 filling marks/month (Type 2), hence a total rate of 655 marks/month; • February 2015-April 2015: 432 new looting marks/month (Type 1) and 255 filling marks/month (Type 2), hence a total rate of 687 marks/month; • April 2015-June 2015: 535 new looting marks/month (Type 1) and 245 filling marks/month (Type 2), hence a total rate of 780 marks/ month.
For a correct interpretation of the above figures, it is crucial to note that these numbers and rates refer to occurrences of new looting and filling marks in each pair of consecutive ST scenes, but do not reflect the total number of looting holes across the site. Rates are calculated with specific regard to each time interval and thereby the number of occurrences that we mapped does not follow an incremental or cumulated rationale. As demonstrated in Fig. 10, looting marks can be the result of repeated looting on previously looted areas. A dynamic analysis therefore proves suitable to address the challenge of monitoring sites such as Apamea where looting manifests in varied forms and at rates varying in time. The identification and mapping of looting marks therefore provide an exceptional added value to the analysis and monitoring of the conditions of the archaeological site, by revealing location and rates of looting occurrence even within already looted areas.

Conclusions
Our experiment over the archaeological site of Apamea allows us to draw three main conclusions.
Firstly, by exploiting the very high resolution imaging capability provided by the novel TerraSAR-X Staring Spotlight (ST) mode, we demonstrated that satellite SAR imagery can now offer a valid and reliable opportunity to monitor looting from space. With such a level of spatial detail, SAR can complement well-established methods based on the analysis of aerial photography and optical imagery, and even fill the gap when the latter are unavailable or are discontinuous, such as in remote sites and areas of conflict, or hindered by cloud cover. Similar to their typical signature in optical imagery, the morphology of looting holes in the satellite SAR scenes is the key element to retrieve patterns of occurrence and change of looting activities. The matrix of interpretation keys for such Fig. 10. Interpretation keys for looting marks in radar backscattering change detection map: (a) un-looted bare ground that is still unaltered by looting; (b) looting mark associated to new looting activity (mark Type 1); (c) looting mark generated by reworking and filling of existing looting area (mark Type 2); (d) unchanged looting hole; (e) cluster of looting holes. TerraSAR-X ST data © DLR 2014-2015; Google Earth Images © CNES/Astrium. marks is, for the first time, defined in this paper and implemented to map new looting features occurring across the area of interest.
Secondly, the regular revisiting time of the TerraSAR-X ST stack (not always available for studies based on optical images) allowed precise estimation of looting rates in order to help to understand the dynamic evolution of this phenomenon across the site investigated. Knowledge of the type, magnitude and spatial distribution of looting occurrences is essential to assess the exposure of local heritage and landscape, as well as the residual risk. As reported in recent papers in the specialist literature (e.g., AAAS, 2014), it is apparent that this is the kind of information sought by archaeologists, conservators, heritage bodies and international organisations to quantify the phenomenon and the resulting damage to archaeology and environment.
Last but not least, the multi-temporal analysis of the TerraSAR-X ST time series presented in this paper contributes to updating the information and records on the condition of Apamea between the end of 2014 and mid-2015 that no other remote sensing platform captured at such resolution level as of mid-2015. The evidence gained from ST image texture, temporal variability and change detection maps is consistent with that observed from Google Earth until mid-2014. From the patterns and observed rates of looting, it is apparent that in the period October 2014-June 2015 the site was still under threat, with major concentration of looting activities north of the two main decumani. The bi-monthly sampling of ST images proved to be sufficiently frequent to observe appreciable differences in time. This reflects the scale with which looting occurred in these months in Apamea, although it is envisaged that acceleration or deceleration of the phenomenon in the subsequent months might require different monitoring timescales.
There are, anyway, some constraints and limitations that need to be accounted for. With regard to the above first conclusion, the demonstration that TerraSAR-X ST (and more generally VHR SAR images acquired from space) can be used to detect looting and measure its rates of occurrence opens the question about the use of these high resolution images in the field of archaeological remote sensing. Data availability can be a constraint because, to the best of our knowledge, at present TerraSAR-X ST images are the only source of VHR data from civil space missions and can be accessed from either announcements of opportunity or under commercial routes (for detail the reader should refer to Airbus Defence & Space and the TerraSAR-X mission websites). Certainly, consistency of acquisition parameters and temporal regularity of ST are advantageous properties that, coupled with very high resolution, enhance the contribution that these SAR data can offer to complement optical imagery towards a quantitative assessment of looting occurrence and rates. The availability of VHR optical images nearly simultaneous to our TerraSAR-X ST data stack would have been ideal to further analyse the complementarity of radar and optical data in this framework. Future research could focus on a synchronised monitoring experiment with both technologies to compare their performance in estimating and measuring looting dynamics.
In the perspective of future implementation, our change detection method is suitable to highlight not only new looting marks in areas previously untouched by looters, but also already looted areas that are further excavated, altered or filled (see Section 5 and matrix of interpretation keys in Fig. 10). A potentiality of this method therefore relates to both the discovery of new sites or areas affected by looting and the updating of looting scenarios in already looted sites.