Integration of geophysical datasets by a conjoint probability tomography approach : application to Italian active volcanic areas

We expand the theory of probability tomography to the integration of different geophysical datasets. The aim of the new method is to improve the information quality using a conjoint occurrence probability function addressed to highlight the existence of common sources of anomalies. The new method is tested on gravity, magnetic and self-potential datasets collected in the volcanic area of Mt. Vesuvius (Naples), and on gravity and dipole geoelectrical datasets collected in the volcanic area of Mt. Etna (Sicily). The application demonstrates that, from a probabilistic point of view, the integrated analysis can delineate the signature of some important volcanic targets better than the analysis of the tomographic image of each dataset considered separately. Mailing address: Dr. Paolo Mauriello, Dipartimento di Scienze e Tecnologie per l’Ambiente ed il Territorio, Università degli Studi del Molise, Campobasso, Italy; email: mauriello@unimol.it Vol51,1,2008_DelNegro 16-02-2009 21:28 Pagina 167


Introduction
Interpretation in geophysics is currently done using a single geophysical dataset to obtain an image of a single geophysical parameter.The geologic parameters of interest are then usually estimated by overlaying a series of these geophysical images.Clearly, much can be learned from experiments designed to isolate those data which are thought to be most sensitive to some particular property of a geological structure.Since current estimation techniques are, however, very subjective and the geologic parameters are not derived directly, much effort has been made to develop codes for joint inversion of geophysical data.This approach is considered of benefit to structural assessment, as it provides a common platform for visualizing the results, thus avoiding overlaying disparate images to look for common sources of anomalies (Voorhies, 1991).Curiously, the conceptual foundation for this type of inversion has sometimes been questioned.However, it must once again be stressed that some properties of geological materials can contribute signals to many kinds of data, yet they are apparently not uniquely determined by any single kind of datum.In such cases, more plausible estimates of the properties might be obtained by using more than one kind of datum.In the last two decades, there have been different approaches to joint inversion of disparate datasets with varying degrees of success (e.g.Lines et al., 1988;Sasaki, 1989;Haber and Oldenburg, 1997;Berge et al., 1999;Hoon et al., 2001;Kozlovskaya, 2001;Nunnari et al., 2001;Chen and Hoversten, 2003;Musil et al., 2003;Gallardo and Meju, 2007).
This paper deals with the integration problem using probability tomography as a basic tool to look for common sources of anomalies.Probability tomography is a 3D imaging approach useful to explore the information potential of a geophysical dataset.It was originally formulated to single out the points underground with the highest occurrence probability of polar sources of self-potential anomalies observed on the ground (Patella, 1997a,b).The approach was then extended to the geoelectric (Mauriello et al., 1998;Mauriello andPatella, 1999a, 2005a), natural electromagnetic induction (Mauriello andPatella, 1999b, 2000), gravimetric (Mauriello and Patella, 2001a,b), ground deformation (Iuliano et al., 2001) and magnetic (Iuliano et al., 2001(Iuliano et al., , 2002a,b;,b;Mauriello and Patella, 2005b) methods.A further expansion of this imaging approach was addressed to single out the contribution of dipolar sources of anomalies (Iuliano et al., 2002b;Mauriello and Patella, 2006a,b).
In the following sections, after recalling the general theory of the 3D probability tomography, we formalise the problem of the integration of different geophysical datasets, previously outlined in Iuliano et al. (2002b), and introduce the concept of conjoint occurrence probability as a tool to highlight the existence of common sources of anomalies.The integration method is tested on some geophysical datasets collected in the Italian active volcanic areas of Mt.Vesuvius (Naples) and Mt.Etna (Sicily).

Theory
Consider a reference coordinate system with the (x, y)-plane placed at sea level and the z-axis positive downward, and a 3D datum domain V as drawn in fig. 1.In particular, V is bounded at top by a non-flat ground survey area and at bottom by a surface through the maximum depths at which datum points can be located.Finally, let A(r) be a vector anomaly field at a set of datum points rϵ(x,y,z), with r∈V.
We assume that A(r) can be discretized as (2.1) i.e. as a sum of partial effects due to Q sources.If the generic q-th source, located at rqϵ(xq, yq, zq), represents the effect of a polar source, the multiplier aq gives directly the pole strength, whereas if it represents a dipolar source, then aq is explicated as (pq•∇q), i.e. the inner product of the dipole moment pq by the vector operator ∇q (Mauriello and Patella, 2006a).The polar or dipolar source effect at r is analytically expressed by the vector kernel s(r, rq), which can reduce to a scalar function in many cases of geophysical interest.
We define the power Λ associated with A(r) within V as Indicating with Ã(r) the vector anomaly field normalized by the square root of its power Λ, viz.
(2.3) and recalling eq.(2.1), can be written as , (2.4) where .From eq. (2.2) it follows (2.5) which, using eq.(2.4), can be expanded as . (2.6) Taking the generic q-th integral in eq.(2.6) and applying Schwarz inequality, accounting for the identity (2.5), we obtain (2.7) From eq. (2.7) we define as anomaly source occurrence probability the function η(rq) given by (Mauriello and Patella, 2006a) (2.8) with (2.9) The 3D η(rq) function, which due to inequality (2.7) satisfies the condition (2.10) is assumed to give a measure of the probability with which either a polar or a dipolar source, responsible of the observed A(r) anomaly field, is present at r q.
The concept of probability associated with η(rq) is motivated as follows.In general, a probability measure p is defined as a function assigning to every subset E of a space of states U a real number p(E) such that (Gnedenko, 1979) (2.11c) Considering that the presence of a source at rq is independent from the presence of another source at another point, the function (2.12) can be defined as a probability density, allowing a probability measure to find a source at rq to be deduced in agreement with axioms ((2.11a)-(2.11c)).Actually, eq.(2.8) differs from eq. (2.12)only for an unknown constant factor and also has the advantage of giving the sign of the source.Thus, we can conventionally assume function (2.8) as an occurrence probability measure of a source of anomaly.It is worth recalling the physical meaning of the sign of the η-function.If the anomaly field A(r) is inspected in terms of polar sources, the sign of η indicates either the sign of the source (e.g. the electrical charge in the selfpotential method), or the sign of the departure of the physical parameter describing the observed field from a reference value (e.g. the density contrast in the gravity method or the resistivity contrast in the geoelectrical method).If, instead, the anomaly field A(r) is inspected in terms of dipolar sources, the algebraic sign indicates the direction of the dipole component along the axis to which the η-function refers, i.e. the direction along which either a source sign inversion (as, e.g., in the self-potential and magnetic methods), or a parameter contrast sign inversion (as, e.g., in the gravity and geoelectrical methods) occurs.

Practical procedure
The 3D tomography imaging procedure of a dataset collected on an uneven topography consists in a scanning procedure based on the knowledge of the function s(r,rq) which takes the role of scanner function.
In practice, since we do not know the anomaly source distribution generating the observed field A(r), we use an elementary positive source of unit strength to scan the whole solid space (the tomospace) to search where the actual sources are most probably located.The result of the cross-correlation in eq.(2.8) for any tern xq, yq, zq of the tomospace will give the probability of occurrence of a positive source (η>0) or a negative source (η<0) in that point as responsible for the observed A(r) surface field.By scanning the whole tomospace we can at last obtain a full 3D image reconstruction of the source distribution underground in a probabilistic sense.
In current survey practice, ground-based or airborne geophysical data are collected over a surface, hence the 3D V-domain, drawn in fig. 1, collapses to a 2D S-domain.Volume integrals in previous equations are reduced to surface integrals extended over S, which, for the sake of generality, is assumed to be a non-flat surface.The η(rq) function is now written as (2.13) where (2.14) Assuming 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, as in fig. 1, using a topography surface regularization factor g(z) given by (2.15) eq.(2.13) and eq.(2.14) can be regularized as (2.16) and (2.17) where the integration intervals along the x and y-axis are [−X, X] and [−Y, Y], respectively.
Ar srr

#
For further details about the regularization procedure in the case of other current datum domains, e.g.profile lines or pseudosections, the reader is referred to Mauriello and Patella (2006a,b).

The conjoint probability tomography
Suppose we dispose of N different geophysical datasets, which we wish to integrate to enhance the geophysical prediction power by finding a set of common anomaly sources located at the same points underground.We indicate with (n=1,2,…, N) the n-th vector anomaly field normalized by the square root of its power Λn, which, recalling eq. ( 2.4), is explicated as Of course, since the N fields to be integrated may not necessarily have non-vanishing sources at the same points, Q is now assumed as the maximum number of points where we can hypothesize the existence of sources of all the N geophysical datasets.
Given eq.(2.5), it follows (3.2) which, using eq.(3.1), can be expanded as Taking the generic q-th addendum from eq. (3.3) and applying Schwarz inequality, accounting for identity (3.2), we obtain which is used to define a conjoint anomaly source occurrence probability function ηN(rq) as (3.9) The function ηN(rq) is now interpreted as the probability that all source elements an responsible for the respective An fields (n=1,2,..,N) simultaneously exist at a given point rq.We point out that the ãn,q distributions are unknown.Thus, as discussed in Section 2.2, we use elementary positive sources of unit strength to scan the tomospace.Under this condition, the concept of probability originally introduced with eq.(2.8) also meets the requirement that the conjoint probability that N independent events obtain is equal to the product of the N single probabilities.
From the practical point of view, the computation of the conjoint anomaly source occurrence probability given in eq.(3.5) is a very simple task, thanks to the property expressed by eq.(3.7).It suffices to make, point by point, the algebraic product of the single occurrence probabilities, rewritten as in eq.(3.8), calculated by the procedure given in Section 2.2.
Basically, our integrated approach aims at revealing from multiple datasets the localization of common sources of anomalies, independently of their strengths, in a probabilistic sense.Hence, time variations of the strengths of the sources, not coupled with a displacement or deformation of their geometry, have no influence on the results of the conjoint probability tomography.
It is also worth pointing out that the N datasets must not necessarily correspond with N different geophysical methods.Indeed, repeated surveys with a single method may be required in % # some places either to monitor the persistence of a potentially migrating source of anomaly, or to reveal weak signals of a stable anomaly source in areas with a highly varying noise.

The Vesuvius volcanic area
Vesuvius is considered a high risk active volcanic system, because of its history of recurrent devastating eruptions in the last 2000 years and the high density of villages all around its lower slopes.In the recent period (1631-1944), the volcanic activity showed a cyclic behavior alternating a rather continuous summit activity, ending with explosive eruptions, to short quiet periods never exceeding seven years (Scandone et al., 1993).Since 1944 Mt.Vesuvius has remained dormant.
In the Mt.Vesuvius volcanic area, gravity (GR), magnetic (MG) and Self-Potential (SP) datasets were collected in the second half of the past century.In particular, a detailed GR map was elaborated by Ciani et al. (1960), which has since then remained almost unchanged, as documented by the many integrations made over the years until the last survey performed at the end of the past century (Berrino et al., 1998), all showing full consistency with the original dataset.The MG dataset refers to an aeromagnetic survey carried by AGIP in 1978 (Cassano and la Torre, 1987), which was repeated 21 years later (Supper and Seiberl, 2000), again showing a substantially unvaried field.Finally, the only SP survey available was performed in the mid 90s of the past century (Di Maio et al., 1998).Although the GR, MG and SP datasets span over about four decades of observations, the GR and MG time stability makes the integration we are going to discuss referable to a common period of the volcano's history.Figure 2 shows a sketch map of the Mt.Vesuvius district and the areas which have been investigated by the GR, MG and SP methods.We first describe the GR-MG and then the GR-SP conjoint tomographies.
Figures 3 and 4 show the results of the 3D probability tomography applied separately on the GR and MG datasets, using the procedure recalled above and previously developed in Mauriello and Patella (2001a,b) for the GR method, and in Iuliano et al. (2001Iuliano et al. ( , 2002a,b) ,b) and Mauriello and Patella (2005b) for the MG method.A detailed discussion of the tomographic images in figs.3 and 4 has recently been given by Iuliano et al. (2002a,b).Figure 5 shows the conjoint GR-MG tomography obtained using eq.(3.2) for N=2.The dominant feature appearing in the slice sequence is the same as that visible in the single GR and MG tomographies.In fact, the central, vertically elongated nucleus in fig. 5 is a replica of that appearing both in fig. 3 and in fig.4, in the depth range from 0.4 km a.s.l down to 2 km b.s.l., with comparable high occurrence proba-bilities.We associate this feature with the volcano chimney.The Vesuvius plumbing system is in fact interpreted as made of a single dominant central conduit, whose top portion can by this analysis be interpreted as plugged by magnetized volcanic materials with density lower than the reference crustal density.
In this example, as previously said, the integrated tomography does not provide any additional information which is not already deducible from each separate tomography.It is indeed proposed as a test on the capability of the new method to single out common sources while depressing signals not referable to common sources, as those appearing only in the GR tomography around the central nucleus.
The second application regards the integration of the GR and SP survey data.Figures 6 and  7 show the results of the 3D dipolar probability tomography applied separately on the GR and SP datasets, following the procedure recalled above and amply developed in previous papers (Iuliano et al., 2002b;Mauriello and Patella, 2006a,b).The GR dipolar tomography drawn in fig.6 is extracted from the whole GR tomography and refers to the smaller area surveyed by the SP method (see fig. 2).Also for these single tomographies a detailed discussion can be found in Iuliano et al. (2002b).Figure 8 shows the conjoint GR-SP dipolar tomography obtained again using eq.(3.2) for N=2.A W-E alignment of nuclei with non-vanishing probabilities crossing midway the whole area is well delineated, especially in the deepest slices from 1.2 km to 2 km b.s.l.. Accordingly, a second NE-SW alignment manifests a similar probability increase in the same depth interval.Both alignments appear to be an extension, much better focussed both horizontally and at greater depths, of the polarized patches visible in the SP dipolar tomography of fig. 7. The interpretation is that the two features may represent the GR-SP conjoint signature of two nearly perpendicular fracture systems, at whose intersection the Vesuvius volcano very likely started to grow.

The Mt. Etna volcanic area
Etna is located in northeastern Sicily and is Fig. 3.The Mt. Vesuvius case-history.The gravity pole source probability tomography.The top slice is the reference Bouguer anomaly survey map (redrawn after Cassano and La Torre, 1987).The left-hand tomographic column refers to the depth range from 0.8 km a.s.l. to 2 km b.s.l. with a slicing step of 0.4 km, whereas the right-hand column refers to the depth range from sea level to 7 km b.s.l. with a slicing step of 1 km.Fig. 4. The Mt. Vesuvius case-history.The magnetic dipole source probability tomography.The top slice is the reference aeromagnetic anomaly survey map (redrawn after Cassano and La Torre, 1987).The left-hand tomographic column refers to the depth range from 0.8 km a.s.l. to 2 km b.s.l. with a slicing step of 0.4 km, whereas the right-hand column refers to the depth range from sea level to 7 km b.s.l. with a slicing step of 1 km.  the largest and most active volcano in Europe.
It formed about a zone of complex geodynamics related to the subduction of the African plate beneath the Eurasian plate.Its volcanism is thought to have been activated by trends cutting transversally the main foredeep-foreland system (Barberi et al., 1973).In this volcanic area, gravity (GR) and Dipole Geoelectrical (DG) datasets were collected in the mid eighties of the past century.Figure 9 shows a sketch map of the Mt.Etna district surveyed by the GR method and the centres and expansion axes of the DG soundings.
Figure 10 shows the residual Bouguer GR  map, redrawn after Loddo et al. (1989), who adopted a reference density of 2.67 g/cm 3 .This map has been used to elaborate the polar tomography shown in fig.11 by the procedure reported in Mauriello and Patella (2001a,b).Noteworthy in this tomography is the presence, in the central volcanic area, of two positive nu-clei aligned E-W and separated by a saddle, completely surrounded, inland, by negative spots.
Figure 12 shows the DG sounding diagrams, redrawn after Loddo et al. (1989).This set of curves has been used to perform the 3D probability tomography drawn in fig.13   procedure described in Mauriello and Patella (1999a), assuming as reference resistivity the average apparent resistivity of 500 Ωm.The tomography of fig. 13 shows in the central volcanic area the existence of a single positive nucleus stretching seaward, hosted within a completely negative background.

by the
The comparison between the single GR and DG tomographies does not show a fully conformal collocation of the dominant GR and DG positive nuclei in the central area.In fact, the core of the DG positive nucleus falls exactly over the saddle separating the two GR positive nuclei, and the western GR positive nucleus partially invades the negative DG area west-  ward.As suggested by Mauriello et al. (2004), a plausible explanation for such a density increase, where the resistivity appears to decrease, may be the presence of a dense, hydrothermally induced mineral particles deposition within the top fractured portion of a deeply rooted barrier, which has already been evidenced in the central volcanic area down to 30 km of depth, at least, by a magnetotelluric study, and assimilated to a slowly cooled magma dyke (Mauriello et al., 1997(Mauriello et al., , 2004)).
Figure 14 shows the conjoint GR-DG polar tomography.The central barrier appears now imaged by a pair of nuclei with opposite sign aligned E-W, deepening down to 3 km b.s.l., at most, immersed inland in a positive background.The barrier can thus be interpreted as made of three different flanked blocks, belonging to the summit portion of the deeply rooted dyke.Starting from the east, 1) the weakly positive eastern nucleus may be ascribed to a mostly sound block (density and resistivity greater than the respective reference values); 2) the central area below the Etna cone, characterized by almost null conjoint probabilities, corresponding to the central positive nucleus in fig.13 and the saddle area in fig.11, may be ascribed to a fractured block with empty intersti-  tial spaces (density equal and resistivity greater than the respective reference values); 3) the negative western nucleus, corresponding to the western GR positive nucleus in fig.11, may finally be ascribed to a fractured block with filled interstitial spaces due to mineral paragenesis (density greater and resistivity lower than the respective reference values).

Conclusions
Adhering to the propensity interpretative approach of modern science, a global probability tomography method has been developed to analyse single and/or integrated geophysical field datasets.The aim of the new method is to try to overcome the limits imposed by current deterministic approaches, which are either based on preconceived assumptions or dictated by more or less idealized models of the geophysical reality, generally much more complex than one may think.The few examples illustrated above demonstrate that, from a probabilistic point of view, acceptable articulated solutions to a complex geophysical interpretation problem can readily be elicited from the whole spectrum of potential solutions displayed by 3D probability tomography, without imposing external constraints.

Fig. 1 .
Fig. 1.The 3D datum domain, characterized by irregular boundary surfaces.The (x,y)-plane is placed at sea level and the z-axis points into the earth.

Fig. 2 .
Fig. 2. A geographical sketch map of the Mt.Vesuvius volcanic area with indication of the gravity, magnetic and self-potential surveyed areas.

Fig. 5 .
Fig. 5.The Mt. Vesuvius case-history.The integrated gravity and magnetic source probability tomography.The top slices are the reference GR and MG survey maps, already shown in fig.3 and fig.4,respectively.The left-hand tomographic column refers to the depth range from 0.8 km a.s.l. to 2 km b.s.l. with a slicing step of 0.4 km, whereas the right-hand column refers to the depth range from sea level to 7 km b.s.l. with a slicing step of 1 km.

Fig. 6 .
Fig. 6.The Mt. Vesuvius case-history.The gravity dipole source probability tomography.The top slice is the scaled reference Bouguer map, extracted from the larger map already shown in fig.3.The following tomographic column refers to the depth range from 0.8 km a.s.l. to 2 km b.s.l. with a slicing step of 0.4 km.

Fig. 7 .
Fig. 7.The Mt. Vesuvius case-history.The self-potential dipole source probability tomography (redrawn after Di Maio et al., 1998).The top slice is the reference SP survey map.The following tomographic column refers to the depth range from 0.8 km a.s.l. to 2 km b.s.l. with a slicing step of 0.4 km.

Fig. 8 .
Fig. 8.The Mt. Vesuvius case-history.The integrated gravity and self-potential dipolar source probability tomography.The top slices are the reference GR and SP anomaly survey maps, already shown in fig.6 and fig.7, respectively.The following tomographic column refers to the depth range from 0.8 km a.s.l. to 2 km b.s.l. with a slicing step of 0.4 km.

Fig. 9 .
Fig. 9.A geographical sketch map of the Mt.Etna volcanic area corresponding to the gravity survey area, with indication of the sounding centers (dots) and expansion axes (arrows) of the dipole geoelectrical survey.

Fig. 11 .
Fig. 11.The Mt. Etna case-history.The gravity source probability tomography.The top slice is the reference topographic map.The following tomographic column refers to the depth range from 1 km b.s.l. to 8 km b.s.l. with a slicing step of 1 km.

Fig. 13 .
Fig. 13.The Mt. Etna case-history.The resistivity source probability tomography.The top slice is the reference topographic map.The following tomographic column refers to the depth range from 0.5 km a.s.l. to 4 km b.s.l. with a slicing step of 0.5 km.

Fig. 14 .
Fig. 14.The Mt. Etna case-history.The integrated gravity and resistivity source probability tomography.The top slice is the reference topographic map.The following tomographic column refers to the depth range from 1 km b.s.l. to 4 km b.s.l. with a slicing step of 1 km.