1. Introduction
The archaeological site of Egnazia is located on the Adriatic coast of Puglia region (Southern Italy) [
1] (
Figure 1). The coastal landscape is characterized by a sequence of sub-horizontal surfaces sloping from the Murge Plateau (locally about 400 m high) towards sea from an elevation of about 120 m. These surfaces are the result of marine abrasion during the Middle Pleistocene because of the superimposition of regional uplift and the glacio-eustatic sea level changes [
2]. In addition, a relict drainage network characterized by deep valleys, locally named gravine, parallel to each other and perpendicular to the coastline, were shaped for sapping processes [
3,
4]. The stratigraphy of the Egnazia area is characterized at the base by the presence of the Calcare di Bari unit (Cretaceous), which represents the local basement, marked by karstic processes. Biocalcareous sandstone related to the Calcarenite di Gravina unit (Lower Pleistocene) overlies in transgression the limestone and outcrops extensively along the cost [
5].
Mentioned by authors such as Pliny, Strabone, Horace, the city of Egnazia had great importance in the ancient world due to its geographical position. Thanks to the presence of the port and the Via Traiana, it was, in fact, an active center of trade and commerce.
The history of ancient Gnathia has unfolded over many centuries [
6]. The first settlement, consisting of a village of huts on the hill of the acropolis, was built in the 15th century BC (Bronze Age). The village itself was rebuilt in the 14th century BC following a large fire [
7]. Clayey beaten, post holes and a refractory cooking surface carved directly into the rock, attest to the presence of nuclei in different areas of the area dating back to the late Bronze Age [
8].
In the 11th century BC (Iron Age), the invasion of populations from the Balkan area, the Iàpigi, is recorded, while with the 8th century BC begins the Messapian phase. From this period, besides the necropolis along the coast on the sides of the acropolis and in the southern area [
9], the mighty defense walls remain, 7 m high, 2 km long, delimiting the three inland sides of a nearly trapezoidal urban area of about 40 hectares, being the fourth side coinciding with the shoreline (
Figure 1). It enclosed a discontinuous residential fabric with non-straight road layouts, alternating with agricultural spaces, pastures, areas of worship and funerary ones, a typical arrangement of the Messapian cities of Salento [
7]. The fortification dates back to the beginning of the 4th century BC and consists of a double retaining wall, an embankment of stones, a wall in square blocks and a ditch. At the end of the 4th century BC or at the beginning of the 3rd century BC the ditch was filled, a double row wall was built and a new and wider ditch was dug, as proved by an archaeological excavation in 1978 in the northernmost edge of the walls [
10]. The new defensive structure had therefore covered and incorporated the old one, thus increasing the safety of the city.
Starting from the 3rd century BC, Egnazia will cease with the Roman occupation and the most important urban development occurs. The period of growth is accentuated with the assumption of the status of municipia after the end of the social war in 89 BC, a condition that suggests the need for an arrangement of administrative centers adapted to the new needs [
6]. The main signs of this change concern the civil basilica, the forum, the enhancement of the port (starting with Agrippa), the creation of two piers and, probably, also the construction of the so-called cryptoporticus [
7]. In the late Republican age and the early imperial age (2nd century BC–1st century AD), a more systematic organization took place. In fact, the arrangement of the road network, in particular from via Minucia, allows for a more defined subdivision of residential districts and public areas (Del Monte & Randino, 2014): the new insulae thus include houses, artisanal structures and commercials; the burial areas undergo transformations, with the expansion of the existing necropolis and the creation of the western necropolis that extends outside the walls; in the Augustan age there was the organization of the core of the forum with the construction of the public baths on the southern side [
6,
7].
A new functional renovation is promoted by Trajan between 108 and 110 AD; of all the initiatives, the most significant is the reorganization of the main artery which will take its name from the emperor and resumes the oldest road route of the Via Minucia.
Only towards the end of the 4th century AD, the city undergoes a period of pause in settlement following the destruction of numerous sectors caused by traumatic events of natural origin (earthquake with a tsunami of 365 AD). Already from the first years of the 5th century, new developments took place, including those of artisanal and commercial systems, the monumental strengthening of Christian buildings, the promotion of economic activities and commercial exchanges, architectural interventions [
6,
7].
In a similar context, thanks to the absence of modern overlays, geophysical prospection can be a useful tool for the subsurface mapping of buried features enhancing the knowledge of the archaeological sites. Thus, in order to ascertain the presence of the ditch along a portion of the southern side of the fortification, a geoelectrical survey has been performed, consisting of nine parallel 2D Electrical Resistivity Tomography (ERT) lines, located as in
Figure 1.
In this paper, the Extended data-adaptive Probability-based Electrical Resistivity Tomography Inversion (E-PERTI) [
11] has been applied for imaging the buried features. E-PERTI is the latest development within the probability tomography imaging approach as adapted to geoelectrical methods [
12,
13,
14,
15].
The probability tomography imaging was at the beginning used to characterize the occurrence probability of resistivity departures from a reference background resistivity. Although not capable of providing estimates of the true resistivity of the buried bodies, it proved useful to highlight the shape and position of anomalous resistivity highs and lows in different geological contexts [
16,
17,
18].
The next step was the launch of the data-adaptive Probability-based ERT Inversion (PERTI) method, which was formulated with the aim of giving an estimation of the most probable resistivity values within a surveyed volume as weighted averages of the whole apparent resistivity dataset [
19]. Using this approach, different near-surface prospections were successfully processed and interpreted [
20,
21,
22,
23,
24].
In general, the PERTI method cannot be considered a standard tomographic inversion scheme. It does not make any reference to the geometry rank of the anomaly sources as an a priori constraint to the imaging algorithm, unlike other methods that were developed to work on this a priori characteristic of the expected geometry of the model such as the minimum gradient support inversion [
25,
26], the blocky inversion (L1 norm) [
27], or the covariance-based inversion [
28]. Conceptually, the PERTI method relates only to the pure physical aspects of the electric field interaction with the subsoil’s structures as it is sensed on the ground. More precisely, the most probable resistivity values are estimated point by point in the surveyed volume, based exclusively on the field apparent resistivity dataset. The approach is thought independent of the true geometry of the bodies, the full resolution of which stands in principle well beyond the intrinsic illumination power of the geoelectrical method.
The E-PERTI method has at last been developed in order to enhance robustness to noise with respect to the original PERTI, as well as to optimize resistivity estimates. It consists of processing, with the PERTI algorithm, many different subsets of data extracted from the apparent resistivity dataset both randomly and by sequential vertical scanning and horizontal windowing across the datum space [
11]. The result is a more or less dense cluster of resistivity values in each point of the surveyed volume, used at last to predict the best estimation of the most probable resistivity in the same point by an intrinsic linear regression model with ordinary least squares techniques.
In this paper, the first application of the E-PERTI scheme to a large apparent resistivity field dataset is illustrated. In the following, at first, a description of the ERT survey technique and a summary of the E-PERTI tool starting from the original PERTI theory are provided, then the explanation of the ERT plan designed in the Egnazia survey site and the presentation and discussion of the results reached from the ERT data inversion are shown.
2. The ERT Survey Technique
Geoelectrical surveys are among the most suitable geophysical tools for the study of the subsoil. The aim is to get information regarding the geometry and localization of buried structures in the light of their intrinsic electrical resistivities. A geoelectrical field survey consists of measuring on the ground the potential difference across a pair of electrodes, generated by a direct current injected in the subsoil through a second pair of electrodes. By changing the positions and spacing of the electrodes according to specific rules, horizontal and vertical resistivity changes in the subsoil can be investigated. Presently, the ERT approach is used routinely, thanks to the latest technological developments that make it a fast, versatile and cost-effective method for depicting resistivity structures. The dipole-dipole (DD) electrode configuration moved on the ground along a profile is the most suited array for gathering ERT data. A dipole is a paired electrode set with the two electrodes located relatively close to one another. The convention for a DD array is to maintain an equal amplitude for both the current (active) dipole and the potential (passive) dipole along the profile, with the distance between the inner electrodes being an integer multiple of the dipole amplitude (see
Figure 2).
A primary advantage of the DD array is the ease of deployment in the field due to the short wire lengths needed to connect the active and passive dipoles to the measuring ammeter and voltmeter, respectively. Furthermore, it allows fast data acquisition as it supports multiple receiver channel measurements simultaneously. The DD layout is also largely more sensitive to lateral resistivity contrasts than other conventional arrays. For this reason, it can be particularly effective in archaeological contexts whose goal is often to highlight anthropogenic structures (walls, cisterns, tombs, trenches, etc.) that cause strong lateral resistivity contrasts at the boundary with the hosting subsoil. Nonetheless, it should be pointed out that the dipole-dipole array is very poor in mapping horizontal structures. Furthermore, the main disadvantage of this array is the decrease in signal strength as the distance between the dipole pair increases. However, in ERT prospecting for archaeology, both for the typical small size and closed shape of the target bodies and for the modest depths involved, these drawbacks appear largely minimized.
The first-instance result from a DD ERT survey is the so-called pseudosection, which is a straightforward representation of the well-known apparent resistivity parameter. On a not uniform subsoil, the apparent resistivity results to be a sort of complex weighted average of the true resistivities of the buried bodies, varying as the dipoles disposition along the profile. Referring to
Figure 2, for a generic
nth measurement with the current electrodes in the position A
n and B
n, and the potential electrodes in the position M
n and N
n, the apparent resistivity
ra,n is calculated as
ra,n =
πk(
k + 1)(
k + 2)
a(
Vn/
In), where
In is the intensity of the direct current injected into the ground from the active dipole and
Vn is the potential drop across the passive dipole. Each
ra,n value is attributed at the midpoint between the two dipoles at a depth half their mutual distance, as in
Figure 2. The plotted
ra,n values are then contoured and colored by any conventional interpolation software. The resulting pseudosection represents a rough and sometimes severely distorted image of the true resistivity distribution of the subsurface. It thus needs to be converted by an inversion software into a true resistivity section in order to create a realistic image of the subsoil.
3. Outline of the E-PERTI Method
The E-PERTI method [
11] is an extension of the PERTI method [
19] that in turn was directly derived from the resistivity probability tomography approach [
13]. The latter consisted in the determination of the resistivity anomaly occurrence probability function
calculated in the space below the surface composed of
M contiguous cubic cells with equal volume Δ
V, each identified by the location of its center (
) (
m = 1, 2, …,
M) and its intrinsic resistivity
. In addition, the dataset was supposed to be composed of
N apparent resistivity measurements,
, each identified by a running index
n (
n = 1, 2, …,
N) determined by the position of the electrode device on the ground. Thus, the
was defined as:
where
In Equation (1), is the difference between and the reference apparent resistivity , computed at the same node as for using a reference model, and is the difference between the intrinsic resistivity and the resistivity of the reference model in the mth cell.
At the base of the PERTI method, there is the assumption that the reference uniform resistivity is supposed to be the unknown value
that corresponds to a generic
mth cell centered at (
), allowing to rewrite
as reported in [
14]:
The rationale for the PERTI is that if at a point (
) it results
, then in the cell centered at (
) the probability to find an increase or a decrease of the resistivity with respect to
is zero. In other words, in that cell, the intrinsic resistivity does not differ from
. Thus, referring to Equation (2), since
is always different from zero, the
condition leads to
which represents the required solution for the application of the PERTI method. We only need to change repeatedly the coordinates (
) to retrieve, point by point, the resistivity pattern within
V [
19].
Unlike deterministic data inversion tools, starting from the pioneering work by Loke and Barker [
29] onward, the PERTI method follows a probabilistic paradigm, totally independent of a priori information. Therefore, in principle, the two approaches are not comparable, despite the fact that in common practice, decisions are taken on the basis of the ERT reconstructions anyway produced. It is just to comply with this practical reason, that a comparison between PERTI and the well-known ERTLAB and RES2D/3DINV commercial software was made on synthetic models, showing that PERTI is generally as efficacious as the two previous tools in detecting, distinguishing and shaping the sources of the observed apparent resistivity anomalies [
19].
Shortly, the main features of the PERTI method are [
19]: (i) unnecessity of a priori reference models; (ii) full, unconstrained adaptability to any kind of dataset, acquisition technique and spatial regularity; (iii) drastic reduction of computing time, thus allowing for in-field real-time inversion.
It is now worth expressing a few words about the first of the PERTI peculiarities listed above. Specifically, we refer to the concept of depth of investigation (DOI), which is associated with an empirical index, based on at least two classic inversions of the dataset, attributing different values to the resistivities of the reference model [
30,
31,
32]. The DOI is considered essential to provide information on how well the subsurface volumes are covered by the data and to discern between similar laterally correlated resistivity values. With PERTI, however, since it is a purely data-adaptive method not needing any reference model to start with, such a DOI index cannot be calculated.
Therefore, in order to improve the accuracy of the method, the extension of PERTI (E-PERTI) [
11] was developed considering the Bernoulli conceptual scheme in probability theory [
33]. According to this scheme, each single geoelectrical measurement is considered as a “test”, that is, the realization of a specific set of conditions that can give rise to an elementary event of a space
U of elementary events (test space). Grouping a set of
N tests such as a resistivity dataset, each measurement is part of a new space of tests, called
UN, composed by arbitrary points of the space
U. In addition, in ideal conditions, as every
nth geoelectric measurement is an individual process, we can consider the
N apparent resistivity tests as automatically satisfying the condition of mutual independence. Hence, if
N given tests are independent, then any
Q tests arbitrarily extracted from the given ones are still independent [
33].
Assuming
Nq ≤
N for
q = 1, 2, …,
Q and
Q are arbitrarily large, a set of
Nq data (tests) within the
N available data is extracted. Applying the PERTI approach, for each
Nq set Equation (3) becomes
Thus, letting
we obtain
Q estimates of
, say
from which to derive the best estimator of
, using, for example, the standard least squares procedure for the determination of the slope of a linear equation of the type
where
Some results are now worth being recalled from [
11] about the main advantages of extending the PERTI method, that is, about the role of the E-PERTI linear least squares best fit expressed by Equation (8), which is directly based on the PERTI formula (4).
Using a DD ERT simulation, it was shown that by extracting many different sets of Nq data from the N-data original pseudosection, with Nq < N, more or less large spreads of the true resistivity values were observed in all of the common points of the Nq-data PERTI sections. In principle, such a spread was assumed as clear evidence of how much an inversion process may depend in an unpredictable way on the profile extent, array setup and data density. At the same time, it was proof of how important is to dispose of a ductile, reliable statistical tool, such as the E-PERTI, straight resulting from the adopted inversion process.
It was also shown that by corrupting the same N-data pseudosection with a random error in one sector, the E-PERTI section resulting from combining a set of Nq-data PERTI sections, was visibly a more realistic image than the original N-data PERTI section, as to position and true resistivity of the buried bodies. Such a result has been assumed as proof of the robustness of the proposed scheme against noise.
In conclusion, E-PERTI has been revealed to be a dynamic and objective tool, whose results can eventually be used as reference models in any other inversion algorithm.
5. Results and Discussion
In
Appendix A, for each profile, the following pictures are shown:
Final E-PERTI sections resulting from the estimation of
ρm with a linear best-fit procedure applied point by point to the resistivity values belonging to the PERTI sections obtained by vertical scanning. In fact, particular emphasis was given to the vertical scanning approach with the ultimate goal of obtaining horizontal two-dimensional sections at gradually increasing depths compatibly with the acquisition parameters of the tomographies. For this type of objective, the procedure of the vertical scanning of the sections seemed to us the most suitable in order to show the horizontal trend of the ditch (
Figure A3,
Figure A6,
Figure A9,
Figure A12,
Figure A15,
Figure A18,
Figure A21,
Figure A24 and
Figure A27).
A common logarithmic scale has been used for visualization of both apparent resistivity and modeled real resistivity maps and, in order to enhance the resistivity contrasts, the same color sequence within the same [ρa,min, ρa,max] interval has been assigned to each set.
Regarding the E-PERTI application, the rationale at the base of the new approach is to inspect in greater detail the electrical nature of the subsoil analyzing different subsets of an ERT dataset. To this respect, it is also worth recalling that with the original PERTI tool the resistivity values that result from using the whole ERT dataset were initially assumed as the most probable values, according to Equation (4). With the E-PERTI, instead, the best estimate of the most probable resistivity values, following Equation (8), is the result from a linear least-squares best-fit of a more or less spread distribution of resistivity values obtained in each single point by the PERTI tool applied to all selected subsets, the whole dataset included.
As regards the modality of selection of the subsets from the original ERT dataset, the horizontal stepwise windowing has the meaning of repeatedly applying the PERTI algorithm to different limited portions of a pseudosection [
11]. The longer the pseudosection, the higher can be the number of windows opened along it, which can be contiguous or even overlapping. Of course, in each window, the relative subset vertically includes all the data from the surface down to the latest pseudodepth. The purpose is to get, from PERTI, resistivity estimates inside any window not influenced by points outside. This modality applied to the Egnazia case study has allowed us to draw the left-side images in
Figure A2,
Figure A5,
Figure A8,
Figure A11,
Figure A14,
Figure A17,
Figure A20,
Figure A23 and
Figure A26.
The second conceived approach is the vertical scanning, which means that the PERTI algorithm is repeatedly applied to horizontal stripes of a pseudosection with progressively increasing thickness starting from top-down to the latest pseudodepth, that is, the whole pseudosection [
11]. Therefore, the aim of the vertical scanning is just to observe how deeper points as stripes become thicker influence resistivity at shallow levels. For the case under study, the results of this approach are shown in the middle sections in
Figure A2,
Figure A5,
Figure A8,
Figure A11,
Figure A14,
Figure A17,
Figure A20,
Figure A23 and
Figure A26. In the images, the last section obtained by a vertical scanning of the pseudosection has not always been processed using the entire data set so the result differs slightly from the final E-PERTI sections.
The last considered is the random extraction approach, consisting in applying the PERTI algorithm repeatedly to different subsets made of a varying number of
Nq data (tests) extracted from the whole pseudosection dataset (
Nq <
N) [
11]. The main purpose is to check for the stability of the resistivity nuclei in the various sections. Wandering or sporadic nuclei are likely to represent artifacts because of their variability in shape and intensity. At Egnazia, this last approach has provided the right-hand sequence of images in
Figure A2,
Figure A5,
Figure A8,
Figure A11,
Figure A14,
Figure A17,
Figure A20,
Figure A23 and
Figure A26.
By means of these entire PERTI results, the last disposes of a cluster of resistivity values in each point of the section, which can even manifest some large dispersion, as amply shown in [
11]. The spreading of the values is a consequence of the different roles that each partial PERTI can have in providing the most probable value of the true resistivity in a point. It follows that a statistical estimation of it is needed, using the best-fit linear regression tool, according to Equation (8). A more refined approach would be the use of a weighted best-fit regression, by assigning a damping or enhancing multiplier to all those resistivity values coming from visibly noise-corrupted PERTI, and/or from strongly depth-dependent PERTI sections if the research interest is, for example, addressed to the shallowest subsoil portion, as in this case study. Although a weighted best-fit may be seen as introducing some subjectivity into the process, in principle it does not conflict with the axioms of probability [
11]. In field practice, a careful choice of different weights can help achieve a sounder interpretation of the results. Under this aspect, E-PERTI reveals to be a dynamic tool, which can also be of support in other inversion algorithms.
As the thickest stripes that have been processed by the vertical scanning approach coincide with the full pseudosections, the relative PERTI sections (bottom section of the central sequence in
Figure A2,
Figure A5,
Figure A8,
Figure A11,
Figure A14,
Figure A17,
Figure A20,
Figure A23 and
Figure A26) can be used for comparison to evaluate the quality of the E-PERTI sections. At glance, the E-PERTI sections appear rather similar to the PERTI ones. However, a careful inspection allows one to recognize the ability of the new method to isolate and better define the limits of the buried conductive mass, as mainly in ERTs 1, 5, 6, 7, 9 (
Figure A3,
Figure A15,
Figure A18,
Figure A21 and
Figure A27). Given that the 3D reconstruction of this shallow confined conductive body represents the final goal of the ERT survey at Egnazia, it can reasonably be assumed to correspond with the loose sandy-clayey sediments filling the ancient ditch.
In
Figure 3, a 3D view of the E-PERTI sections is reported. As depicted by a dotted line, the high-conductivity volumes are distributed along a route that is parallel to the walls and that can be associated with the filling of the ancient ditch of the city, which had been supposed by archaeologists to lie close to that line.
Finally, by assembling all DD profiles, the E-PERTI solution of the entire apparent resistivity dataset has provided a 3D ERT model that in
Figure 4 is shown as a sequence of horizontal slices at increasing depth from 2 m down to 5 m beneath the ground level, while in
Figure 5 is depicted in full compact form. In both these last two figures, in order to better highlight the resistivity lows attributable to the presence of the ditch, the color scale has been modified, by compressing the left-hand half of the log scale below 1.8 and leaving the right-hand half colorless.
6. Conclusions
In this paper, for the first time, an ERT case study in archaeology in which the E-PERTI method was fully adopted to process and interpret the field survey dataset was shown. The purpose was to detect and model the course of the buried ancient ditch hypothesized to exist just outside a southern portion of the walls of the 4th century BC, built in defense of the ancient urban area of Egnazia (Puglia, Italy).
With the E-PERTI method, an interpretative model derived from the statistical determination of the most probable resistivity values within the surveyed space was proposed. This approach is theoretically justified by the probabilistic paradigm upon which the original PERTI method is based, which allows the extraction of a wealth of independent data subsets from the field dataset, following Bernoulli’s conceptual scheme in probability theory.
In practice, the most probable model resistivity at each point of the surveyed space has been determined by a least-squares approach aimed at finding the slope of a linear equation that relates the multiplicity of the provisional resistivity estimates to the partial data subsets. Needless to say that the term “most probable” used to qualify the resistivity values of the model is a relative concept, as it intrinsically depends on amount and accuracy of the measures forming the entire initial database.
At last, the present field case history allows the confirmation of what was stated in the previous methodological study [
11]. That is to say that E-PERTI is a robust approach to noise dejection and that the technique of extraction of subsets of the data is an easy viable technique for optimizing the modeling process in ERT applications. Ultimately, if necessary, the result deriving from the application of the E-PERTI method can be used as a valuable reference model in any standard iterative inversion process.