TY - JOUR AU - Mauriello,, Paolo AB - Abstract We present a generalized theory of the probability tomography applied to the gravity method, assuming that any Bouguer anomaly data set can be caused by a discrete number of monopoles, dipoles, quadrupoles and octopoles. These elementary sources are used to characterize, in an as detailed as possible way and without any a priori assumption, the shape and position of the most probable minimum structure of the gravity sources compatible with the observed data set, by picking out the location of their centres and peculiar points of their boundaries related to faces, edges and vertices. A few synthetic examples using simple geometries are discussed in order to demonstrate the notably enhanced resolution power of the new approach, compared with a previous formulation that used only monopoles and dipoles. A field example related to a gravity survey carried out in the volcanic area of Mount Etna (Sicily, Italy) is presented, aimed at imaging the geometry of the minimum gravity structure down to 8 km of depth bsl. 3D probability tomography, gravity method, multipole analysis, Mt Etna volcanic area 1 Introduction Geophysical probability tomography (GPT) is an approach to virtually explore the subsoil in the search for the most probable localization of the sources of anomalies appearing in any given datum domain. GPT was originally formulated for the self-potential method (Patella 1997a, 1997b) and then extended to the geoelectrical, EM induction, magnetic and gravity methods (Mauriello and Patella 1999a, 1999b, 2001a, 2001b, 2008a, 2008b). In all these formulations, the causative bodies were assimilated to aggregates of either monopoles or dipoles, or both. The GPT sensitivity and resolution have been largely tested on experimental data sets in many application fields, generally with a good performance (Di Maio et al1998, Lapenna et al2000, Saracco et al2004, Takács et al2005, Matias et al2006, Bhattacharya et al2007, Chianese and Lapenna 2007, Alaia et al2008a). An extension of the GPT approach has recently been proposed for the geoelectrical method (Alaia et al2008b), by assuming a data set to be the conjoint response of monopoles, dipoles and quadrupoles. The geoelectric GPT algorithm has been tested on synthetic examples and then successfully applied to a field survey performed in the historical park of Pompei (Naples, Italy) for archaeological purposes. A further extension of the GPT approach that includes also octopoles has been elaborated for the self-potential method (Alaia et al2009). These simplest sources have been used to find the most probable location of the centres and some peculiar points of the boundaries of the source bodies. In synthesis, the motivation for the GPT method was to find an imaging solution that only has characteristics that are directly traced to the data, using the most elemental physical sources of the geophysical signals and without introducing a priori information. This solution, which can be interpreted as a minimum structure model (Kitanidis 1997), although not necessarily providing an accurate image of the real source bodies since it neglects to utilize structural information that may be available, is of fundamental importance and may be useful as a benchmark. The aim of this paper is to adapt the new formalism to the specific properties of the gravity field. For an easier understanding of the following abstract mathematical formulation, it is worth giving soon a simple explanation of the structural meaning we intend to associate with each of the elementary sources that will be considered. A monopole will be introduced to characterize a mass centre. A dipole will be used to identify a point where a sharp density discontinuity occurs in the direction of only one of the coordinate axes of a rectangular reference system, parallel to the dipole axis. A square quadrupole will be used to identify a point where a sharp density discontinuity occurs in the directions of only two of the coordinate axes, parallel to the sides of the quadrupole. A cubic octopole will at last be used to identify a point where a sharp discontinuity occurs in the directions of all the three coordinate axes, parallel to the edges of the cube. All these multipole elements in 0D, 1D, 2D and 3D can be included within the special class of mass elements formalized by Jacoby and Smilde (2009). Such an extended multipole analysis can thus provide, in principle, a complete set of critical points necessary to better characterize the geometry of the most probable minimum structure source, compatible with the measured data set. The new approach is thus expected to be useful in characterizing discrete, confined real bodies, but less apt, of course, to directly conform to unconfined structures. The further limitation must be recalled as to the inability of the GPT approach to provide estimates of the density contrasts between anomalous bodies and hosting material (Mauriello and Patella 2001a, 2001b). A few tests on simple synthetic models and the analysis of a field case in the Mt Etna volcanic area (Sicily, Italy) will be illustrated, in order to evaluate the feasibility and accuracy of the new approach. 2 The GPT generalized formalism for gravity 2.1 The Bouguer anomaly function and its information power Let us consider a reference system with a horizontal (x,y)-plane and the z-axis positive downwards, and a 2D datum domain S as in figure 1. The S-domain is generally a non-flat ground survey area described by a topographic function z(x,y). Indicating with Ba(r) the Bouguer anomaly at a set of datum points r ≡ (x,y,z), with r∈S, we assume that Ba(r) can be discretized as 1 i.e. as a sum of effects due to a set of M monopoles, the mth element of which is located at rm≡ (xm, ym, zm) and has strength pm ⋅Pm, where pm and Pm are the monopole moment and a point operator zero-order tensor, respectively; a set of N dipoles, whose nth element is located at rn≡ (xn,yn,zn) with strength (dun ⋅ Lun) (u = x,y,z), where dun and Lun are the dipole moment and a line operator first-order tensor, respectively; a set of G quadrupoles, whose gth element is located at rg≡ (xg,yg,zg) with strength quvg ⋅ Suvg (u,v = x,y,z), where quvg and Suvg are the quadrupole moment and a square operator second-order tensor, respectively; and a set of H octopoles, whose hth element is located at rh≡ (xh,yh,zh) and has strength ouvwh ⋅ Cuvwh (u,v,w = x,y,z), where ouvwh and Cuvwh are the octopole moment and a cube operator third-order tensor, respectively. Figure 1 Open in new tabDownload slide The datum domain (S-domain) generating the Ba(r) map on top. The (x,y)-plane is placed at sea level and the z-axis points into the earth. Figure 1 Open in new tabDownload slide The datum domain (S-domain) generating the Ba(r) map on top. The (x,y)-plane is placed at sea level and the z-axis points into the earth. The dot in the definition of the source strength tensors indicates the inner product. The operator tensors P, Lu, Suv and Cuvw (u,v,w = x,y,z) are made explicit in figure 2. Figure 2 Open in new tabDownload slide Explicit representation of the symbolic tensor operators appearing in the definition of the strengths of the monopole, dipole, quadrupole and octopole source elements defined in equation (1). Figure 2 Open in new tabDownload slide Explicit representation of the symbolic tensor operators appearing in the definition of the strengths of the monopole, dipole, quadrupole and octopole source elements defined in equation (1). The effect of the M, N, G and H source elements at a point r ∈ S is regulated by the base function s(r, ri) = (zi − z)/|ri − r|3 (i = m,n,g,h), which represents the vertical component of the gravitational acceleration due to a point mass of unitary strength (Mauriello and Patella 2001a, 2001b). We define the information power Λ, associated with Ba(r), over the surface S as 2 which, using equation (1), is expanded as 3 2.2 The source monopole occurrence probability We consider a generic mth integral of the first sum in equation (3) and apply Schwarz inequality, thus obtaining 4 Inequality (4) is used to define a source monopole occurrence probability (SMOP) function (Mauriello and Patella 2001a, Alaia et al2008b) as 5 where 6 and s(r, rm) has the role of a source monopole scanner. The 3D SMOP function, which satisfies the condition −1 ⩽ η(p)m ⩽ 1, is given as a measure of the probability of a source monopole placed at rm, being responsible for the observed Ba(r) field. At each rm, a positive, null or negative value of η(p)m indicates the occurrence probability of a positive, null or negative differential mass (ΔM), respectively. This is equivalent to say the occurrence probability of a density increase, invariance or decrease, respectively, with respect to the reference density. For computational purposes, we assume that the projection of S onto the (x,y)-plane can be fitted to a rectangle R of sides 2X and 2Y along the x- and y-axis, respectively. Using the topography surface regularization factor g(z) given by (Mauriello and Patella 2001a, 2001b) 7 equation (5) is definitely written as follows: 8 with 9 The s(r,rn) function will take the role of a source monopole scanner in the computation procedure that will be described later. 2.3 The source dipole occurrence probability We take a generic nth integral of the second sum in equation (3) and apply, as previously, Schwarz's inequality to each u-component. We can thus define a source dipole occurrence probability (SDOP) function (Iuliano et al2002, Alaia et al2008b) as 10 with 11 where the surface regularization has been accounted for. Also η(d)n,u falls in the range [−1,1]. Thus, at each rn, three values of η(d)n,u can be computed. They are interpreted as a measure of the probability of a single source dipole located at rn, being responsible of the whole Ba(r) field. Each first derivative of s(r, rn) will take the role of a source dipole scanner in the forthcoming computation procedure. The first derivatives of s(r, rn) are given in expanded form in the appendix. 2.4 The source quadrupole occurrence probability Accordingly, we consider now a generic gth integral of the third sum in equation (3) and apply again Schwarz's inequality to each uv-element (u,v = x,y,z), which is arranged to define the source quadrupole occurrence probability (SQOP) function (Alaia et al2008b) as 12 with 13 As previously, also the 3D SQOP function falls in the range [−1,1]. Thus, at each rg, nine values of η(q)g,uv are taken as a measure of the probability for a quadrupole source located at rg, to be responsible of the Ba(r) data set. Since Suvg is a symmetric square tensor, it follows that η(q)g,uv = η(q)g,vu. Hence, at each rg, the three diagonal plus the three right-up or left-down off-diagonal terms of η(q)g,uv are sufficient. However, as we are interested in finding the position of square quadrupoles, we will finally consider only the three off-diagonal terms (u ≠ v) (Alaia et al2008b). Each of the three off-diagonal second derivative of s(r,rg) will have the role of a source quadrupole scanner in the forthcoming computation procedure. The three off-diagonal second derivatives of s(r, rg) are reported in the appendix. 2.5 The source octopole occurrence probability Finally, we consider a generic hth integral of the fourth sum in equation (3) and apply again Schwarz's inequality to each uvw-term (u,v,w = x,y,z), allowing a source octopole occurrence probability (SOOP) function to be defined as 14 with 15 As previously, the 3D SOOP function falls in the range [−1,1]. At each rh, 27 values may be calculated, which are interpreted as a probability measure of a single octopole source located at rh, being responsible of the whole Ba(r) data set. However, as we are interested in finding only the position of cubic octopoles, we will limit our analysis only to the SOOP functions with u≠v, u≠w, v≠w, which reduce to a single SOOP function, since a change of the order of derivation of s(r, rm) does not alter the final result. Thus, we shall consider the SOOP function with u = x, v = y and w = z. The above-defined third derivative of s(r, rh) will take the role of a source octopole scanner in the forthcoming computation procedure. Its expression is reported in the appendix. 2.6 The η-functions as probabilities It is now worth clarifying the role of probability assigned to the η-functions. In general, a probability measure ϕ is defined as a function assigning to every subset E of a space of states U a real number ϕ(E) such that (Gnedenko 1979) 16 if E ∩ F = 0 with E,F ⊂ U, then 17 18 Considering now that the presence of a source at ri (i = m,n,g,h depending on whether the source is a monopole, a dipole, a quadrupole or an octopole) is independent of the presence of another similar source at another point, the following function 19 where V is the total volume containing all the sources, can be defined as a probability density function, since we can deduce from it a probability measure to find a source at ri satisfying axioms (16), (17) and (18). In fact, equation (19) differs from each of the definitions of η(ri) given in equations (8), (10), (12) and (14) for an unknown constant factor and also has the advantage of giving information about the sign of the source. Thus, we assume each of the definitions of η(ri) as our conventional source occurrence probability function. 2.7 The 3D probability tomography calculations The 3D multipole tomography imaging of a Ba(r) data set will consist of a repeated scanning procedure in the lower half-space, using one after the other the source scanner functions previously defined, appearing inside the integrals in equation (8), (10), (12) or (14), depending on whether monopoles, dipoles, quadrupoles or octopoles are to be detected. In practice, using a discretized version of the integrals in equations (8), (10), (12) and (14), the values of the SMOP, SDOP, SQOP and SOOP functions, respectively, are calculated for any assigned ri in the solid half-space. Each value will give the probability which a corresponding positive source (η > 0) or negative source (η < 0) located at that point takes as responsible for the observed Ba(r) data set. Using a regular grid of points ri in the lower half-space, for each source type we can obtain a 3D representation of the buried source distribution in a probabilistic sense. In particular, picking out from each plot the points where the corresponding SMOP, SDOP, SQOP or SOOP function reaches its local maxima and minima, the most probable minimum set of monopoles, dipoles, quadrupoles or octopoles can at last be extracted. Finally, assembling all these plots into a single 3D diagram, the most probable minimum structure of the sources of the observed Ba(r) data set can be imaged. Ideally, this procedure simulates a scanning process of the solid half-space in the search for the points at which to place a tuned set of sources of various polar orders, capable altogether of giving the best approximation to the measured Ba(r) field. 3 Synthetic examples We show a few synthetic examples, in order to highlight the main aspects of this GPT multipole generalization. 3.1 The cube model At first, we consider the cube model, which is assigned a density contrast Δσ = 0.5 g cm−3, sides 6 m long parallel to the coordinate axes and the ΔM centre at x = 0, y = 0 and z = 6 m. The Ba data set has been calculated at the nodes of a square grid with a step of 1 m from −18 m to 18 m along both the x- and y-axis, using the formula for the rectangular parallelepiped reported by Parasnis (1997, p 81, equation 3.27). Figure 3 shows the calculated Ba map on the x,y plane. We recall that an identical Ba image results if a change of the model distance scale is made, provided that the Ba values are all multiplied by the same scale factor. Figure 3 Open in new tabDownload slide The Ba map for the cube model with Δσ = 0.5 g cm−3, sides 6 m long and ΔM centre at x = 0, y = 0 and z = 6 m. Figure 3 Open in new tabDownload slide The Ba map for the cube model with Δσ = 0.5 g cm−3, sides 6 m long and ΔM centre at x = 0, y = 0 and z = 6 m. Figure 4 shows the results from the application of the multipole GPT algorithms to the Ba map in figure 3. For the sake of clarity, in this and all the following 3D examples, sufficiently small SMOP, SDOP, SQOP and SOOP nuclei will be drawn, each enclosing a local maximum absolute value (MAV) of the corresponding η-function. Only the sign of the SMOP nuclei will explicitly be indicated in order to easily get the sign of the density contrast of the sources to which they refer. The algebraic sign of each of the SDOP, SQOP and SOOP nuclei is, instead, readily deducible from that of the corresponding SMOP nucleus, knowing the positive directions of the coordinate axes. Figure 4 Open in new tabDownload slide The SMOP (Aa), x-SDOP (Ab), y-SDOP (Ac), z-SDOP (Ad), xy-SQOP (Ba), xz-SQOP (Bb), yz-SQOP (Bc) and xyz-SOOP (Bd) tomographies derived from the Ba synthetic map in figure 3. The body with light blue lines is the cube model. The SMOP nucleus is positive. Figure 4 Open in new tabDownload slide The SMOP (Aa), x-SDOP (Ab), y-SDOP (Ac), z-SDOP (Ad), xy-SQOP (Ba), xz-SQOP (Bb), yz-SQOP (Bc) and xyz-SOOP (Bd) tomographies derived from the Ba synthetic map in figure 3. The body with light blue lines is the cube model. The SMOP nucleus is positive. The SMOP image shows a positive nucleus around the ΔM source centre. The SDOP image shows, instead, three distinct doublets of nuclei with opposite sign very close to the centres of the corresponding opposite faces of the cube. Three distinct quadruplets appear around the centres of the cube edges in the SQOP tomographies of the off-diagonal terms, and an octoplet located at the vertices of the cube is the peculiar result from the SOOP image. The parameters of the nuclei in figure 4 are listed in table 1. A shift of 0.1 m along the z-axis is estimated for the cube centre from its true position. Furthermore, an average error of about 3% affects the estimate of the edge length of the cube, from the distance between the MAV points of two opposite nuclei in each multiplet. Table 1 Characterization of the SMOP, SDOP SQOP and SOOP nuclei in figure 4. A: nucleus type; B: bounding isosurface; C: maximum absolute value (MAV); D: (x,y,z) of the MAV point. A . B . C . D . SMOP (+) (0.0, 0.0, 5.9) x-SDOP (+) (−3.1, 0.1, 6.0) x-SDOP (−) (3.1, 0.0, 6.0) y-SDOP (+) (−0.1, −3.0, 6.0) y-SDOP (−) (−0.1, 3.0, 6.0) z-SDOP (+) (0.0, 0.1, 3.0) z-SDOP (−) (0.0, 0.0, 9.0) xy-SQOP (+) (3.0, 3.0, 5.9) (−3.0, −3.0, 5.9) xy-SQOP (−) (3.1, −3.0, 6.0) (−3.0, 3.0, 6.0) xz-SQOP (+) (−3.3, 0.0, 3.1) (3.3, 0.0, 9.0) xz-SQOP (−) (3.3, 0.0, 3.1) (−3.3, 0.0, 9.0) yz-SQOP (+) (0.0, −3.2, 3.1) (0.0, 3.3, 9.0) yz-SQOP (−) (0.0, 3.2, 3.1) (0.0, −3.2, 9.0) xyz-SOOP (+) (3.1, 3.1, 3.0) (−3.1, −3.1, 3.0) (3.1, −3.1, 9.0) (−3.1, 3.1, 9.0) xyz-SOOP (−) (3.1, −3.1, 3.1) (−3.1, 3.1, 3.1) (−3.1, −3.1, 9.1) (3.1, 3.1, 9.0) A . B . C . D . SMOP (+) (0.0, 0.0, 5.9) x-SDOP (+) (−3.1, 0.1, 6.0) x-SDOP (−) (3.1, 0.0, 6.0) y-SDOP (+) (−0.1, −3.0, 6.0) y-SDOP (−) (−0.1, 3.0, 6.0) z-SDOP (+) (0.0, 0.1, 3.0) z-SDOP (−) (0.0, 0.0, 9.0) xy-SQOP (+) (3.0, 3.0, 5.9) (−3.0, −3.0, 5.9) xy-SQOP (−) (3.1, −3.0, 6.0) (−3.0, 3.0, 6.0) xz-SQOP (+) (−3.3, 0.0, 3.1) (3.3, 0.0, 9.0) xz-SQOP (−) (3.3, 0.0, 3.1) (−3.3, 0.0, 9.0) yz-SQOP (+) (0.0, −3.2, 3.1) (0.0, 3.3, 9.0) yz-SQOP (−) (0.0, 3.2, 3.1) (0.0, −3.2, 9.0) xyz-SOOP (+) (3.1, 3.1, 3.0) (−3.1, −3.1, 3.0) (3.1, −3.1, 9.0) (−3.1, 3.1, 9.0) xyz-SOOP (−) (3.1, −3.1, 3.1) (−3.1, 3.1, 3.1) (−3.1, −3.1, 9.1) (3.1, 3.1, 9.0) Table 1 Characterization of the SMOP, SDOP SQOP and SOOP nuclei in figure 4. A: nucleus type; B: bounding isosurface; C: maximum absolute value (MAV); D: (x,y,z) of the MAV point. A . B . C . D . SMOP (+) (0.0, 0.0, 5.9) x-SDOP (+) (−3.1, 0.1, 6.0) x-SDOP (−) (3.1, 0.0, 6.0) y-SDOP (+) (−0.1, −3.0, 6.0) y-SDOP (−) (−0.1, 3.0, 6.0) z-SDOP (+) (0.0, 0.1, 3.0) z-SDOP (−) (0.0, 0.0, 9.0) xy-SQOP (+) (3.0, 3.0, 5.9) (−3.0, −3.0, 5.9) xy-SQOP (−) (3.1, −3.0, 6.0) (−3.0, 3.0, 6.0) xz-SQOP (+) (−3.3, 0.0, 3.1) (3.3, 0.0, 9.0) xz-SQOP (−) (3.3, 0.0, 3.1) (−3.3, 0.0, 9.0) yz-SQOP (+) (0.0, −3.2, 3.1) (0.0, 3.3, 9.0) yz-SQOP (−) (0.0, 3.2, 3.1) (0.0, −3.2, 9.0) xyz-SOOP (+) (3.1, 3.1, 3.0) (−3.1, −3.1, 3.0) (3.1, −3.1, 9.0) (−3.1, 3.1, 9.0) xyz-SOOP (−) (3.1, −3.1, 3.1) (−3.1, 3.1, 3.1) (−3.1, −3.1, 9.1) (3.1, 3.1, 9.0) A . B . C . D . SMOP (+) (0.0, 0.0, 5.9) x-SDOP (+) (−3.1, 0.1, 6.0) x-SDOP (−) (3.1, 0.0, 6.0) y-SDOP (+) (−0.1, −3.0, 6.0) y-SDOP (−) (−0.1, 3.0, 6.0) z-SDOP (+) (0.0, 0.1, 3.0) z-SDOP (−) (0.0, 0.0, 9.0) xy-SQOP (+) (3.0, 3.0, 5.9) (−3.0, −3.0, 5.9) xy-SQOP (−) (3.1, −3.0, 6.0) (−3.0, 3.0, 6.0) xz-SQOP (+) (−3.3, 0.0, 3.1) (3.3, 0.0, 9.0) xz-SQOP (−) (3.3, 0.0, 3.1) (−3.3, 0.0, 9.0) yz-SQOP (+) (0.0, −3.2, 3.1) (0.0, 3.3, 9.0) yz-SQOP (−) (0.0, 3.2, 3.1) (0.0, −3.2, 9.0) xyz-SOOP (+) (3.1, 3.1, 3.0) (−3.1, −3.1, 3.0) (3.1, −3.1, 9.0) (−3.1, 3.1, 9.0) xyz-SOOP (−) (3.1, −3.1, 3.1) (−3.1, 3.1, 3.1) (−3.1, −3.1, 9.1) (3.1, 3.1, 9.0) Of practical interest is to retrieve the shape and position of the source body. Figure 5 suggests that a quick modelling can be done by plotting into a single image all the nuclei drawn in figure 4. Figure 5 Open in new tabDownload slide A joint representation of the SMOP (red), SDOP (light blue), SQOP (green) and SOOP (purple) nuclei, viewed from top (a) and laterally (b), useful to retrieve the source body of the Ba map in figure 3. The body with light blue lines is the cube model. The SMOP nucleus is positive. Figure 5 Open in new tabDownload slide A joint representation of the SMOP (red), SDOP (light blue), SQOP (green) and SOOP (purple) nuclei, viewed from top (a) and laterally (b), useful to retrieve the source body of the Ba map in figure 3. The body with light blue lines is the cube model. The SMOP nucleus is positive. 3.2 The sphere model The Ba map in figure 3 has a very close resemblance with the map generated by a spherical body. It is thus instructive to analyse the GPT response of a sphere. For this purpose, since the gravity effects of all the concentric spheres, which have the same mass surplus or deficit with respect to the hosting material, are equivalent to the effect of a positive or negative differential mass (ΔM) concentrated at the common centre, we consider a point source at x = 0, y = 0 and z = 6 m with ΔM = 56.52 × 106 g. The Ba map has been computed at the nodes of a square grid with the same characteristics as in the previous case, using the formula for the sphere reported by Parasnis (1997, p 76, equation 3.20). Figure 6 depicts the Ba map thus obtained. Figure 6 Open in new tabDownload slide The Ba map generated by a point differential mass ΔM = 56.52 × 106 g placed at x = 0, y = 0 and z = 6 m. Figure 6 Open in new tabDownload slide The Ba map generated by a point differential mass ΔM = 56.52 × 106 g placed at x = 0, y = 0 and z = 6 m. Figure 7 shows the results from the application of the multipole GPT imaging. As previously, in the SMOP image a nucleus appears about the ΔM centre. The SDOP images show again three doublets of nuclei with opposite sign aligned along the coordinate axes. The nuclei of the off-diagonal SQOP quadruplets appear again as three distinct groups disposed on planes parallel to the reference planes. Finally, also the SOOP octoplet appears located in a cubic fashion quite similar to the previous case. Considered singularly, the SMOP, SDOP, SQOP and SOOP nuclei are so regularly located that no difference can be detected with respect to the previous cube model. The situation changes notably if we plot the nuclei altogether into a multipole image as in figure 8. It is no longer possible, now, to combine a set of SDOP, SQOP and SOOP nuclei crossed by a single plane as in the previous case. In other words, a cube's face can no longer be traced. The multipole GPT seems thus capable to differentiate the response of a cube from that of a single point source located at the centre of the cube. Figure 7 Open in new tabDownload slide The SMOP (Aa), x-SDOP (Ab), y-SDOP (Ac), z-SDOP (Ad), xy-SQOP (Ba), xz-SQOP (Bb), yz-SQOP (Bc) and xyz-SOOP (Bd) tomographies derived from the Ba map in figure 6. The SMOP nucleus is positive. Figure 7 Open in new tabDownload slide The SMOP (Aa), x-SDOP (Ab), y-SDOP (Ac), z-SDOP (Ad), xy-SQOP (Ba), xz-SQOP (Bb), yz-SQOP (Bc) and xyz-SOOP (Bd) tomographies derived from the Ba map in figure 6. The SMOP nucleus is positive. Figure 8 Open in new tabDownload slide A joint representation of the SMOP (red), SDOP (light blue), SQOP (green) and SOOP (purple) nuclei under two different angles of view, useful to retrieve the source body of the Ba map in figure 6. The SMOP nucleus is positive. Figure 8 Open in new tabDownload slide A joint representation of the SMOP (red), SDOP (light blue), SQOP (green) and SOOP (purple) nuclei under two different angles of view, useful to retrieve the source body of the Ba map in figure 6. The SMOP nucleus is positive. The fact that SDOP, SQOP and SOOP nuclei appear also for the ΔM single point model is fully consistent with the definition of the multipoles of different orders given in the introduction. Qualitatively speaking, this result may be interpreted as another aspect of the non-uniqueness of the inverse problem, i.e. of the model responsible of a given Ba data set. In other words, we may assert that the polyhedral geometry in figure 8 likely represents the most probable multipolar model of minimum structure approximating the Ba map in figure 6. The parameters of the SMOP, SDOP, SQOP and SOOP nuclei are listed in table 2. Table 2 Characterization of the SMOP, SDOP SQOP and SOOP nuclei in figure 7. A: nucleus type; B: bounding isosurface; C: maximum absolute value (MAV); D: (x,y,z) of the MAV point. A . B . C . D . SMOP (+) (0.0, 0.0, 6.0) x-SDOP (+) (−3.0, 0.0, 6.0) x-SDOP (−) (3.1, 0.0, 6.0) y-SDOP (+) (−0.1, −3.0, 6.0) y-SDOP (−) (0.0, 3.0, 6.1) z-SDOP (+) (0.1, 0.1, 3.1) z-SDOP (−) (0.1, 0.1, 9.1) xy-SQOP (+) (2.2, 2.2, 5.9) (−2.2, −2.3, 5.9) xy-SQOP (−) (2.2, −2.3, 5.9) (−2.2, 2.2, 5.9) xz-SQOP (+) (−2.5, −0.2, 4.5) (2.5, −0.1, 7.5) xz-SQOP (−) (2.5, −0.2, 4.5) (−2.5, −0.1, 7.5) yz-SQOP (+) (−0.2, −2.5, 4.6) (−0.1, 2.5, 7.5) yz-SQOP (−) (−0.2, 2.5, 4.5) (−0.1, −2.5, 7.5) xyz-SOOP (+) (2.0, 2.1, 4.5) (−2.0, −2.1, 4.4) (1.9, −2.1, 7.6) (−2.0, 2.1, 7.5) xyz-SOOP (−) (2.1, −2.0, 4.4) (−2.0, 2.0, 4.4) (−2.0, −2.0, 7.5) (2.0, 2.0, 7.5) A . B . C . D . SMOP (+) (0.0, 0.0, 6.0) x-SDOP (+) (−3.0, 0.0, 6.0) x-SDOP (−) (3.1, 0.0, 6.0) y-SDOP (+) (−0.1, −3.0, 6.0) y-SDOP (−) (0.0, 3.0, 6.1) z-SDOP (+) (0.1, 0.1, 3.1) z-SDOP (−) (0.1, 0.1, 9.1) xy-SQOP (+) (2.2, 2.2, 5.9) (−2.2, −2.3, 5.9) xy-SQOP (−) (2.2, −2.3, 5.9) (−2.2, 2.2, 5.9) xz-SQOP (+) (−2.5, −0.2, 4.5) (2.5, −0.1, 7.5) xz-SQOP (−) (2.5, −0.2, 4.5) (−2.5, −0.1, 7.5) yz-SQOP (+) (−0.2, −2.5, 4.6) (−0.1, 2.5, 7.5) yz-SQOP (−) (−0.2, 2.5, 4.5) (−0.1, −2.5, 7.5) xyz-SOOP (+) (2.0, 2.1, 4.5) (−2.0, −2.1, 4.4) (1.9, −2.1, 7.6) (−2.0, 2.1, 7.5) xyz-SOOP (−) (2.1, −2.0, 4.4) (−2.0, 2.0, 4.4) (−2.0, −2.0, 7.5) (2.0, 2.0, 7.5) Table 2 Characterization of the SMOP, SDOP SQOP and SOOP nuclei in figure 7. A: nucleus type; B: bounding isosurface; C: maximum absolute value (MAV); D: (x,y,z) of the MAV point. A . B . C . D . SMOP (+) (0.0, 0.0, 6.0) x-SDOP (+) (−3.0, 0.0, 6.0) x-SDOP (−) (3.1, 0.0, 6.0) y-SDOP (+) (−0.1, −3.0, 6.0) y-SDOP (−) (0.0, 3.0, 6.1) z-SDOP (+) (0.1, 0.1, 3.1) z-SDOP (−) (0.1, 0.1, 9.1) xy-SQOP (+) (2.2, 2.2, 5.9) (−2.2, −2.3, 5.9) xy-SQOP (−) (2.2, −2.3, 5.9) (−2.2, 2.2, 5.9) xz-SQOP (+) (−2.5, −0.2, 4.5) (2.5, −0.1, 7.5) xz-SQOP (−) (2.5, −0.2, 4.5) (−2.5, −0.1, 7.5) yz-SQOP (+) (−0.2, −2.5, 4.6) (−0.1, 2.5, 7.5) yz-SQOP (−) (−0.2, 2.5, 4.5) (−0.1, −2.5, 7.5) xyz-SOOP (+) (2.0, 2.1, 4.5) (−2.0, −2.1, 4.4) (1.9, −2.1, 7.6) (−2.0, 2.1, 7.5) xyz-SOOP (−) (2.1, −2.0, 4.4) (−2.0, 2.0, 4.4) (−2.0, −2.0, 7.5) (2.0, 2.0, 7.5) A . B . C . D . SMOP (+) (0.0, 0.0, 6.0) x-SDOP (+) (−3.0, 0.0, 6.0) x-SDOP (−) (3.1, 0.0, 6.0) y-SDOP (+) (−0.1, −3.0, 6.0) y-SDOP (−) (0.0, 3.0, 6.1) z-SDOP (+) (0.1, 0.1, 3.1) z-SDOP (−) (0.1, 0.1, 9.1) xy-SQOP (+) (2.2, 2.2, 5.9) (−2.2, −2.3, 5.9) xy-SQOP (−) (2.2, −2.3, 5.9) (−2.2, 2.2, 5.9) xz-SQOP (+) (−2.5, −0.2, 4.5) (2.5, −0.1, 7.5) xz-SQOP (−) (2.5, −0.2, 4.5) (−2.5, −0.1, 7.5) yz-SQOP (+) (−0.2, −2.5, 4.6) (−0.1, 2.5, 7.5) yz-SQOP (−) (−0.2, 2.5, 4.5) (−0.1, −2.5, 7.5) xyz-SOOP (+) (2.0, 2.1, 4.5) (−2.0, −2.1, 4.4) (1.9, −2.1, 7.6) (−2.0, 2.1, 7.5) xyz-SOOP (−) (2.1, −2.0, 4.4) (−2.0, 2.0, 4.4) (−2.0, −2.0, 7.5) (2.0, 2.0, 7.5) 3.3 The rotated and tilted cube model We show now what happens when the edges of the cube are no longer parallel to the reference coordinate axes. A new model is thus analysed by rotating the cube previously dealt with by 45° around both the vertical axis and the y-oriented horizontal axis through the ΔM centre. Figure 9 shows the Ba map of this new source body configuration. Figure 9 Open in new tabDownload slide The Ba map for the tilted cube model with the same parameters as in figure 6, rotated by 45° around the vertical and horizontal axes through the ΔM centre. Figure 9 Open in new tabDownload slide The Ba map for the tilted cube model with the same parameters as in figure 6, rotated by 45° around the vertical and horizontal axes through the ΔM centre. Figure 10 illustrates the results from the application of the multipole GPT imaging. The SMOP image still shows a nucleus located around the ΔM centre. In contrast, the SDOP, SQOP and SOOP nuclei exhibit a mixed behaviour compared with that of the coaxial cube model. While in the former case they distinctly represent the faces, edges and vertices of the cube, respectively, now the same multiplets seem to indifferently simulate any of such features. The only apparent regularity is that the SDOP, SQOP and SOOP nuclei are still arranged in three doublets, three quadruplets and one octoplet, respectively. However, it must be stressed that this behaviour is not arbitrary or a numerical artefact. Since our procedure simply implies the search for the MAV points of the first derivatives and second- and third-order cross-derivatives of the base function with respect to the reference axes, it can easily be checked how effectively the positions of all the MAV points exactly comply with the meaning associated with the multipoles of various order, reported in the introduction. Figure 10 Open in new tabDownload slide The SMOP (Aa), x-SDOP (Ab), y-SDOP (Ac), z-SDOP (Ad), xy-SQOP (Ba), xz-SQOP (Bb), yz-SQOP (Bc) and xyz-SOOP (Bd) tomographies derived from the Ba synthetic map drawn in figure 9. The body with light blue lines is the inclined cube model. The SMOP nucleus is positive. Figure 10 Open in new tabDownload slide The SMOP (Aa), x-SDOP (Ab), y-SDOP (Ac), z-SDOP (Ad), xy-SQOP (Ba), xz-SQOP (Bb), yz-SQOP (Bc) and xyz-SOOP (Bd) tomographies derived from the Ba synthetic map drawn in figure 9. The body with light blue lines is the inclined cube model. The SMOP nucleus is positive. In order to evaluate the goodness of our approach, figure 11 shows an assembled version of all the MAV nuclei from two different angles of view. It can easily be appreciated how the position and shape of the GPT-reconstructed minimum structure fully conforms to the original model. Of course, for this simple case of only one source body, the minimum structure coincides with the body itself. Figure 11 Open in new tabDownload slide A joint representation of the SMOP (red), SDOP (light blue), SQOP (green) and SOOP (purple) nuclei, under two different angles of view, useful to retrieve the source body of the Ba map in figure 9. The SMOP nucleus is positive. Figure 11 Open in new tabDownload slide A joint representation of the SMOP (red), SDOP (light blue), SQOP (green) and SOOP (purple) nuclei, under two different angles of view, useful to retrieve the source body of the Ba map in figure 9. The SMOP nucleus is positive. 3.4 The two-prism model The fourth example is the coaxial two-prism model, whose aim is to test the resolution power of the new GPT method. The first prism is simply a cube with Δσ = 0.5 g cm−3 and the second one is a parallelepiped with Δσ = −1.0 g cm−3. Three cases are shown with three different distances between the ΔM centres of the two prisms. Positions and edge lengths of the two bodies are detailed in the caption of figure 12. The Ba data sets have been computed at the nodes of a square grid by a 1 m step in the rectangle [−60,60] × [−30,30] m2. Figure 12 Open in new tabDownload slide The Ba map for the two-prism model made of a cube with Δσ = 0.5 g cm−3, edge 12 m long and centre at x = 20 m (a), x = 14.5 m (b) and x = 11 m (c), y = 4 m, z = 15 m, and a parallelepiped with Δσ = −1.0 g cm−3, x- and z-oriented edges 13 m long and y-oriented edges 13.5 m long, and centre at x = −20 m (a), x = −15 m (b) and x = −11.5 m (c), y = −4.75 m, z = 15 m. Figure 12 Open in new tabDownload slide The Ba map for the two-prism model made of a cube with Δσ = 0.5 g cm−3, edge 12 m long and centre at x = 20 m (a), x = 14.5 m (b) and x = 11 m (c), y = 4 m, z = 15 m, and a parallelepiped with Δσ = −1.0 g cm−3, x- and z-oriented edges 13 m long and y-oriented edges 13.5 m long, and centre at x = −20 m (a), x = −15 m (b) and x = −11.5 m (c), y = −4.75 m, z = 15 m. Figure 12 shows the Ba maps for the three cases in order of decreasing distance between the centres from the top (a) to the bottom plot (c). It is quite evident that the decreasing distance is the cause of an increasing compression of the Ba contour lines in the region of highest mutual interference. Figure 13 displays the GPT results for the three cases, where, for brevity only, the combined multipole images are reported. In the top one, which refers to the distance between the centres greater than three times the average edge length of the bodies, the interaction between the two prisms is rather negligible and their true shape can still be recognized. In the middle image, which refers to the distance between the centres of about 2.5 times the average edge length, all the facing SDOP, SQOP and SOOP nuclei depart from their initial places to converge to the centre of the two bodies' system. Finally, in the bottom picture, which refers to a distance a little greater than two times the average edge length, the detached facing nuclei of the same type and sign, indicated in the figure, are wholly melted midway between the two prisms. The facing faces, edges and vertices of the two nearby bodies have therefore completely lacked resolution. Figure 13 Open in new tabDownload slide A joint representation of the SMOP (red and yellow), SDOP (light blue), SQOP (green) and SOOP (purple) nuclei for the two-cube model with decreasing distance between the centres of the two cubes. The sequence of the images is the same as that of the Ba maps in figure 12. Figure 13 Open in new tabDownload slide A joint representation of the SMOP (red and yellow), SDOP (light blue), SQOP (green) and SOOP (purple) nuclei for the two-cube model with decreasing distance between the centres of the two cubes. The sequence of the images is the same as that of the Ba maps in figure 12. The melting effect is an intrinsic limitation of the method. Equation (3) and following definitions of the η-functions suggest that mutual dependences exist on sources inside and outside the volume probed. Interferences are likely to vanish as long as the source bodies are so largely separated as to produce well distinct anomalies. It is difficult at this stage of the research to establish a general rule to fix the threshold below which melting starts to dominate. Accounting for the example illustrated here, we might tentatively say that a threshold may be estimated at a horizontal spacing between the two facing faces of the prisms of the same order of the mean depth of their centres. 4 A field example We show that the application of the 3D GPT to a gravity survey of Mt Etna (Eastern Sicily, Italy), carried out in the frame of a multi-method geophysical project, aimed at delineating the structural setting of the whole volcanic apparatus. Etna is the biggest and most active volcano in Southern Europe, which formed within a large extension zone related to the subduction of the African under the Eurasian Plate. Figure 14 shows the survey area and the relative Ba residual map, obtained from the Ba field map after the application of a 2D high-pass filter of 50 km cut-off wavelength. The reference mean crustal density of 2.67 g cm−3 was taken for slab and terrain corrections. Further details on field data acquisition and processing are given in Loddo et al (1989). The multipole GPT scanning process has been elaborated down to 8 km of depth bsl. Figure 14 Open in new tabDownload slide The top map shows the Mt Etna survey area with indication of the magnetotelluric survey area (red rectangle) and sounding stations (green triangles) and the passive seismic tomography profiles (blue lines). The (bottom map) shows the residual Ba map (redrawn after Loddo et al (1989)). Numbers in the Ba map indicate the dominant anomalies. Figure 14 Open in new tabDownload slide The top map shows the Mt Etna survey area with indication of the magnetotelluric survey area (red rectangle) and sounding stations (green triangles) and the passive seismic tomography profiles (blue lines). The (bottom map) shows the residual Ba map (redrawn after Loddo et al (1989)). Numbers in the Ba map indicate the dominant anomalies. The SMOP image in figure 15 displays a set of distinct nuclei containing the monopoles, which are interpreted as the ΔM centres responsible of the Ba closed anomalies in figure 14. Combining in pairs and altogether the SMOP, SDOP, SQOP and SOOP nuclei into single pictures, regardless of their algebraic sign except for the SMOP nuclei, the images in figure 16 are obtained. Figure 15 Open in new tabDownload slide The Mt Etna 3D SMOP tomography of the residual gravity map reported in figure 14. Figure 15 Open in new tabDownload slide The Mt Etna 3D SMOP tomography of the residual gravity map reported in figure 14. Figure 16 Open in new tabDownload slide A joint representation of the SMOP and SDOP (a), SMOP and SQOP (b), SMOP and SOOP (c), SMOP, SDOP, SQOP and SOOP (d) nuclei, resulting from the application of the GPT method to the Mt Etna gravity map in figure 14. Figure 16 Open in new tabDownload slide A joint representation of the SMOP and SDOP (a), SMOP and SQOP (b), SMOP and SOOP (c), SMOP, SDOP, SQOP and SOOP (d) nuclei, resulting from the application of the GPT method to the Mt Etna gravity map in figure 14. Figure 16 shows a complex assemblage of the SDOP, SQOP and SOOP multiplets. All the related nuclei appear clustered around the eight monopoles of figure 15, tending thus to configure eight distinct blocks. Some remarkable features can now be noted in this 3D block pattern. Compared with the results from the previous synthetic examples with simple confined bodies, the first feature is a rather frequent incompleteness of the multiplets around the monopoles. For instance, the three doublets of the dipole sources never total six nuclei, and in one case, around the monopole n.1 in figure 15, they even drop to only three, one for each doublet. A similar situation occurs for the 3 quadruplets, which drop to only 5 nuclei from a total of 12, around the monopole n.3 in figure 15, and for the octoplet, which drops to only 2 nuclei from a total of 8, around the monopole n.4 in figure 15. Also, we do not happen to see in each cluster a number of multiplets, belonging to any of the source types, greater than the corresponding maximum number as detected by the GPT of the synthetic cases previously considered. Another feature is the presence of the detachment and melting effects due to closeness of the causative bodies, as explained along with the discussion of the two-prism model. This feature can particularly be viewed in figure 16, midway between the monopoles n.6 and n.7 and the monopoles n.2 and n.5, marked in figure 15. No tilts can be appreciated around any horizontal axis, while moderate horizontal rotations can be observed for the sources n.3, n.4 and n.8, reported in figure 15. Summarizing, the biggest central blocks n.5, n.6 and n.7 appear as confined bodies with a roughly quadrangular section. All the other peripheral blocks are reasonably to be interpreted as partially unconfined bodies, because of either an insufficient number of data or lack of density contrast at least along one of the reference directions. In order to get an initial evaluation of both feasibility and accuracy of the multipole GPT method in real cases, a comparison is briefly illustrated between the above results and those from magnetotelluric soundings (Mauriello et al2004) and passive seismic tomography (Chiarabba et al2000), located as in the top map of figure 14. To facilitate the comparison, table 3 lists only the estimated depths of the monopole and bottom vertical dipole for each of the eight isolated sources plotted in figures 15 and 16, because they cannot easily be deduced from figure 16. Table 3 Estimated depths of the monopole and bottom vertical dipole for each of the eight isolated sources plotted in figures 15 and 16. Source number . Monopole depth in km bsl . Bottom vertical dipole depth in km bsl . 1 4.0 >8.0 2 3.6 >8.0 3 2.6  5.2 4 2.7  6.0 5 2.8  5.8 6 3.1  6.1 7 4.8 >8.0 8 3.0  5.8 Source number . Monopole depth in km bsl . Bottom vertical dipole depth in km bsl . 1 4.0 >8.0 2 3.6 >8.0 3 2.6  5.2 4 2.7  6.0 5 2.8  5.8 6 3.1  6.1 7 4.8 >8.0 8 3.0  5.8 Table 3 Estimated depths of the monopole and bottom vertical dipole for each of the eight isolated sources plotted in figures 15 and 16. Source number . Monopole depth in km bsl . Bottom vertical dipole depth in km bsl . 1 4.0 >8.0 2 3.6 >8.0 3 2.6  5.2 4 2.7  6.0 5 2.8  5.8 6 3.1  6.1 7 4.8 >8.0 8 3.0  5.8 Source number . Monopole depth in km bsl . Bottom vertical dipole depth in km bsl . 1 4.0 >8.0 2 3.6 >8.0 3 2.6  5.2 4 2.7  6.0 5 2.8  5.8 6 3.1  6.1 7 4.8 >8.0 8 3.0  5.8 The denser blocks n.5 and n.6 in figure 16 appear to conform to two weakly differentiated highly resistive nuclei with resistivity ⩾4000 Ω m, modelled by the magnetotelluric soundings about in the same area and within the same depth interval, as shown in figure 17. A correspondence can also be noticed with the results of the W–E and N–S passive seismic tomography sections. The N–S section, which laps the western face of the multipole source n.6, discloses a high P-wave velocity bulge with velocity as high as 6.2 km s−1, raising from more than 15 km of depth up to about 4 km bsl, as clearly evident in the top section of figure 18. It shows a lateral extent of about 10 km, which is nearly the same as that of the gravity block n.6 along the x-axis. A similar P-wave bulge-shaped anomaly with same depth features appears in the W–E tomography section, which laps the northern face of the gravity source n.6, as depicted in the bottom section of figure 18. Also in this case, the lateral extent of the bulge-shaped seismic anomaly of about 10 km results to be nearly the same as that of the block n.6 along the y-axis. Figure 17 Open in new tabDownload slide Resistivity pattern at different depths bsl as deduced by the magnetotelluric soundings in the red rectangular area sketched in the top map of figure 14 (redrawn after Mauriello et al (2004)). Figure 17 Open in new tabDownload slide Resistivity pattern at different depths bsl as deduced by the magnetotelluric soundings in the red rectangular area sketched in the top map of figure 14 (redrawn after Mauriello et al (2004)). Figure 18 Open in new tabDownload slide P-wave velocity pattern across the sections through the blue profiles sketched in the top map of figure 14 (redrawn after Chiarabba et al (2000)). The white dots indicate the hypocentres of the volcanic earthquakes. Figure 18 Open in new tabDownload slide P-wave velocity pattern across the sections through the blue profiles sketched in the top map of figure 14 (redrawn after Chiarabba et al (2000)). The white dots indicate the hypocentres of the volcanic earthquakes. The magnetotelluric slices in figure 17 show also a hint of resistivity decrease down to about 20 Ω m at the westmost corner of the surveyed area. This tendency to a resistivity low may be associated with the lower density block n.4, occurring in the same sector. Finally, the N–S P-wave velocity tomography in the top section of figure 18 shows the signature of a bowl-shaped velocity low with a gradual decrease from 5.8 km s−1 to 4.2 km s−1, immediately to the north of the previous velocity high in the depth range 3–10 km bsl and laterally extending for about 6 km. This low may be correlated with the lower density block n.7, whose eastern face is lapped by the N–S seismic profile. A similar feature occurs in the W–E section, immediately to the west of its velocity high, which may again be correlated with the lower density block n.7, whose southern face is now lapped by the W–E profile. The above analysis shows, in conclusion, the probable existence of a dense, rigid and resistive structure placed in the central part of the volcano, which appears to split into two distinct blocks in the upper 8 km of depth. From the volcanological point of view, each of them may reasonably describe a thick set of slowly cooled vertical magma dikes, subject to stress charge and sudden rupture, as documented by the quake hypocentre cluster appearing in the seismic sections of figure 18. The above gravity interpretation is of course far from a decisive modelling. Further refinements would, however, go beyond the main scope of the paper. Nevertheless, it has served the purpose of showing the degree of reliability and accuracy that the new multipole gravity analysis can reach when applied to a real field case. 5 Conclusion We have exposed the theory of the probability tomography for the gravity method by a generalized approach including a combination of monopoles, dipoles, quadrupoles and octopoles as elemental sources of Bouguer anomalies. These physical sources are used to detect the position of the centres of mass of the sources and critical points of their boundaries, by means of a set of ad hoc formulated occurrence probability functions. The points where these functions reach local maxima and minima are claimed to characterize the most probable minimum structure of multipole sources compatible with the anomalies observed in the datum space. We can thus state, in short, that the proposed method is a data adaptive anomaly source imaging tool that optimizes the problem of identifying the positions of critical points subject to the constraint to represent the local maxima and minima of a set of probability functions. Noteworthy is that these functions are directly derived in the frame of the gravity theory. Basically, the proposed method is nothing but a tool for solving an intrinsically non-deterministic decision problem. It must also be pointed out that this probability tomography approach cannot give estimates of the density contrasts which characterize the sources of anomalies. It is thus more appropriate in those circumstances in which the density contrast of the bodies to be detected is either known in advance or not a strictly necessary parameter, and the problem is only to check the existence of the bodies and to reconstruct their geometries. This is often the case in target-oriented near-surface applications to archaeology, mining industry, cavity detection, engineering geology and others, or when the integration with other powerful geophysical methods does not make essential to know the density contrasts. Otherwise, the method can be considered as one of the most valuable supports to the classical interpretation. The individually calculated local maximum absolute values of the SMOP, SDOP, SQOP and SOOP functions can in fact be used as benchmarks in the search for a solution of the canonical gravity inverse problem (Jacoby and Smilde 2009). Improving the geometrical definition of the sources of gravimetric anomalies has thus been the true purpose of the new approach. A few synthetic examples using cubes, spheres and prisms have been illustrated in order to demonstrate the full capabilities of the multipole probability tomography in the search for the most probable geometry and location of the buried sources. Finally, a field case related to a gravity survey carried out in the volcanic area of Mt Etna (Sicily, Italy) has been presented in order to delineate the geometry of the Bouguer anomaly sources in the central volcanic area within the first 8 km of depth below sea level. Both the synthetic and field examples have shown that the multipole probability tomography approach can be very useful in characterizing the geometry of discrete, confined bodies. Of course, unconfined sources of anomalies, such as those caused by e.g. smooth changes in the depth of a crystalline basement, or any other set of layers with greater coverage, will not suitably be imaged due to lack of some of the critical points which can be dealt with as multipole sources. Appendix A.1 A.2 A.3 A.4 A.5 A.6 A.7 A.8 where A.9 Acknowledgments The authors wish to thank an anonymous referee and Professor Dr Wolfgang Jacoby, whose detailed comments and suggestions have notably helped to improve the manuscript. References Alaia R , Patella D , Mauriello P . , 2008a Application of the geoelectrical 3D probability tomography in a test-site of the archaeological park of Pompei (Naples, Italy) , J. Geophys. Eng. , vol. 5 (pg. 67 - 76 )http://dx.doi.org/10.1088/1742-2132/5/1/0071742213217422140 Google Scholar Crossref Search ADS WorldCat Alaia R , Patella D , Mauriello P . , 2008b Imaging quadrupolar geophysical anomaly sources by 3D probability tomography. Application to near surface geoelectrical surveys , J. Geophys. Eng. , vol. 5 (pg. 359 - 70 )http://dx.doi.org/10.1088/1742-2132/5/4/0011742213217422140 Google Scholar Crossref Search ADS WorldCat Alaia R , Patella D , Mauriello P . , 2009 Imaging multipole self-potential sources by 3D probability tomography , Prog. Electromagn. Res. B , vol. 14 (pg. 311 - 39 )http://dx.doi.org/10.2528/PIERB09021614 Google Scholar Crossref Search ADS WorldCat Bhattacharya B , Shalivahan B , Jardani A , Bera A . , 2007 Three-dimensional probability tomography of self-potential anomalies of graphite and sulphide mineralization in Orissa and Rajasthan, India , Near Surf. Geophys. , vol. 5 (pg. 223 - 30 ) OpenURL Placeholder Text WorldCat Chianese D , Lapenna V . , 2007 Magnetic probability tomography for environmental purposes: test measurements and field applications , J. Geophys. Eng. , vol. 4 (pg. 63 - 74 )http://dx.doi.org/10.1088/1742-2132/4/1/0081742213217422140 Google Scholar Crossref Search ADS WorldCat Chiarabba C , Amato A , Boschi E , Barberi F . , 2000 Recent seismicity and tomographic modelling of the Mount Etna plumbing system , J. Geophys. Res. , vol. 105 (pg. B5 10923 - 38 )http://dx.doi.org/10.1029/1999JB90042701480227 Google Scholar Crossref Search ADS WorldCat Di Maio R , Mauriello P , Patella D , Petrillo Z , Piscitelli S , Siniscalchi A . , 1998 Electric and electromagnetic outline of the Mount Somma-Vesuvius structural setting , J. Volcanol. Geotherm. Res. , vol. 82 (pg. 219 - 38 )http://dx.doi.org/10.1016/S0377-0273(97)00066-803770273 Google Scholar Crossref Search ADS WorldCat Gnedenko B V . , 1979 , Kurs Teorii Verojatnostej Moscow Mir published in Italian as Teoria Della Probabilità, Rome: Editori Riuniti Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Iuliano T , Mauriello P , Patella D . , 2002 Looking inside Mount Vesuvius by potential fields integrated geophysical tomographies , J. Volcanol. Geotherm. Res. , vol. 113 (pg. 363 - 78 )http://dx.doi.org/10.1016/S0377-0273(01)00271-203770273 Google Scholar Crossref Search ADS WorldCat Jacoby W , Smilde P . , 2009 , Gravity Interpretation: Fundamentals and Application of Gravity Inversion and Geological Interpretation [with CDROM] Berlin Springer Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Kitanidis P K . , 1997 The minimum structure solution to the inverse problem , Water Resour. Res. , vol. 33 (pg. 2263 - 72 )http://dx.doi.org/10.1029/97WR0161900431397 Google Scholar Crossref Search ADS WorldCat Lapenna V , Patella D , Piscitelli S . , 2000 Tomographic analysis of self-potential data in a seismic area of southern Italy , Ann. Geofis. , vol. 43 (pg. 361 - 74 ) OpenURL Placeholder Text WorldCat Loddo M , Patella D , Quarto R , Ruina G , Tramacere A , Zito G . , 1989 Application of gravity and deep dipole geoelectrics in the volcanic area of Mt. Etna (Sicily) , J. Volcanol. Geotherm. Res. , vol. 39 (pg. 17 - 39 )http://dx.doi.org/10.1016/0377-0273(89)90018-803770273 Google Scholar Crossref Search ADS WorldCat Matias H C , Monteiro Santos F A , Rodrigues Ferreira F E , Machado C , Luzio R . , 2006 Detection of graves using the micro-resistivity method , Ann. Geophys. , vol. 49 (pg. 1235 - 44 ) OpenURL Placeholder Text WorldCat Mauriello P , Patella D . , 1999a Resistivity anomaly imaging by probability tomography , Geophys. Prospect. , vol. 47 (pg. 411 - 29 )http://dx.doi.org/10.1046/j.1365-2478.1999.00137.x0016802513652478 Google Scholar Crossref Search ADS WorldCat Mauriello P , Patella D . , 1999b Principles of probability tomography for natural-source electromagnetic induction fields , Geophysics , vol. 64 (pg. 1403 - 17 )http://dx.doi.org/10.1190/1.14446451070485X Google Scholar Crossref Search ADS WorldCat Mauriello P , Patella D . , 2001a Gravity probability tomography: a new tool for buried mass distribution imaging , Geophys. Prospect. , vol. 49 (pg. 1 - 12 )http://dx.doi.org/10.1046/j.1365-2478.2001.00234.x0016802513652478 Google Scholar Crossref Search ADS WorldCat Mauriello P , Patella D . , 2001b Localization of maximum-depth gravity anomaly sources by a distribution of equivalent point masses , Geophysics , vol. 66 (pg. 1431 - 37 )http://dx.doi.org/10.1190/1.14870881070485X Google Scholar Crossref Search ADS WorldCat Mauriello P , Patella D . , 2008a Localization of magnetic sources underground by a probability tomography approach , Prog. Electromagn. Res. M , vol. 3 (pg. 27 - 56 )http://dx.doi.org/10.2528/PIERM08050504 Google Scholar Crossref Search ADS WorldCat Mauriello P , Patella D . , 2008b Resistivity tensor probability tomography , Prog. Electromagn. Res. B , vol. 8 (pg. 129 - 46 )http://dx.doi.org/10.2528/PIERB08051604 Google Scholar Crossref Search ADS WorldCat Mauriello P , Patella D , Petrillo Z , Siniscalchi A , Iuliano T , Del Negro C . , 2004 A geophysical study of the Mount Etna volcanic area , Mt Etna: Volcano Laboratory Bonaccorso A , Calvari S , Coltelli M , Del Negro C , Falsaperla S . (AGU Geophys. Mon. Series 143) pp 273–91 Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Parasnis D S . , 1997 , Principles of Applied Geophysics London Chapman and Hall Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Patella D . , 1997a Introduction to ground surface self-potential tomography , Geophys. Prospect. , vol. 45 (pg. 653 - 81 )http://dx.doi.org/10.1046/j.1365-2478.1997.430277.x0016802513652478 Google Scholar Crossref Search ADS WorldCat Patella D . , 1997b Self-potential global tomography including topographic effects , Geophys. Prospect. , vol. 45 (pg. 843 - 63 )http://dx.doi.org/10.1046/j.1365-2478.1997.570296.x0016802513652478 Google Scholar Crossref Search ADS WorldCat Saracco G , Labazuy P , Moreau F . , 2004 Localization of self-potential sources in volcano-electric effect with complex continuous wavelet transform and electrical tomography methods for an active volcano , Geophys. Res. Lett. , vol. 31 (pg. 1 - 5 )http://dx.doi.org/10.1029/2004GL01955400948276 Google Scholar Crossref Search ADS WorldCat Takács E , Pethö G , Szabó I . , 2005 Comparative investigations about the applicability of current density pseudosections in the interpretation of 2D VLF vertical magnetic anomalies , Acta Geod. Geophys. Hung. , vol. 40 (pg. 127 - 46 )http://dx.doi.org/10.1556/AGeod.40.2005.2.21217897715871037 Google Scholar Crossref Search ADS WorldCat © 2009 Nanjing Institute of Geophysical Prospecting TI - Imaging multipole gravity anomaly sources by 3D probability tomography JO - Journal of Geophysics and Engineering DO - 10.1088/1742-2132/6/3/009 DA - 2009-09-01 UR - https://www.deepdyve.com/lp/oxford-university-press/imaging-multipole-gravity-anomaly-sources-by-3d-probability-tomography-CATm6CcNoy SP - 298 VL - 6 IS - 3 DP - DeepDyve ER -