Joint inversion of muon tomography and gravimetry - a resolving kernel approach

Both muon tomography and gravimetry are geophysical methods that provide information on the density structure of the Earth's subsurface. Muon tomography measures the natural flux of cosmic muons and its attenuation produced by the screening effect of the rock mass to image. Gravimetry generally consists in measurements of the vertical component of the local gravity field. Both methods are linearly linked to density, but their spatial sensitivity is very different. Muon tomography essentially works like medical X-ray scan and integrates density information along elongated narrow conical volumes while gravimetry measurements are linked to density by a 3-dimensional integral encompassing the whole studied domain. We develop the mathematical expressions of these integration formulas -- called acquisition kernels -- to express resolving kernels that act as spatial filters relating the true unknown density structure to the density distribution actually recoverable from the available data. The resolving kernels provide a tool to quantitatively describe the resolution of the density models and to evaluate the resolution improvement expected by adding new data in the inversion. The resolving kernels derived in the joined muon/gravimetry case indicate that gravity data are almost useless to constrain the density structure in regions sampled by more than two muon tomography acquisitions. Interestingly the resolution in deeper regions not sampled by muon tomography is significantly improved by joining the two techniques. Examples taken from field experiments performed on La Soufri\`ere of Guadeloupe volcano are discussed.


Introduction
Muon tomography is a geophysical method that offers a new way to determine the density of large rock volumes by measuring their screening effect on the natural cosmic muons flux crossing the rock volume to probe (e.g. Nagamine 2003, see Tanaka 2013 for a brief review). The small cross-section of muons in ordinary matter (Barrett et al. 1952) allows the hard component of the muon spectrum (Tang et al., 2006;Gaisser & Stanev 2008) to cross hectometers, and even kilometers, of rock. Most muons crossing the rock volume have a negligible scattering relative to the instrument angular resolution and travel along straight trajectories ranking muon tomography among the class of straight-ray scanning imaging methods (Marteau et al., 2011). In practice, muon tomography is performed by using a series of pixelated particle detectors that allow to determine the trajectories of the muons passing through the rock body. Portable field telescopes presently used sample hundredths of directions and allow to scan an entire volcano from a single view-point in a couple of weeks (Fig. 1). By counting the number of muons passing through the target, the attenuation onto the incident muon flux is determined for each sampled direction and used to produce a radiography of the object opacity (expressed in g.cm −2 ) or of average density along ray-paths if the object geometry is known.
Since the pioneering works by Nagamine et al. (1995a,b) and Tanaka et al. (2001), recent studies illustrate the interest of the method to image spatial and temporal variations of the density inside volcanoes (Tanaka et al. 2005(Tanaka et al. , 2007a(Tanaka et al. ,b, 2008(Tanaka et al. , 2009a Portal et al., 2013). When compared with the relatively large number of publications devoted to qualitative applications of muon tomography, only a small number of studies address methodological issues and quantitative assessments of the method. Lesparre et al. (2010) establish a feasibility formula where the achievable density resolution is related to the measurement duration (i.e. time resolution), the total apparent rock thickness (i.e. total opacity) and the telescope acceptance (i.e. the detection capacity of the matrices). The feasibility formula writes as an inequality and gives practical hints to design field experiments and evaluate which density heterogeneities can be resolved inside a given geological target, for a given amount of time and a given telescope. In a more recent study, Jourde et al. (2013) present experimental evidences of a flux of upward-going particles that occurs in certain field conditions. These particles have trajectories parallel to those of the muons emerging from the rock body to radiography but they travel through the telescope from rear to front. These upward-going particles may constitute a huge Poissonian noise that could strongly alter the radiographies. Jourde et al. (2013) give practical recommendations for choosing experimental sites likely to give an as best as possible signal-to-noise ratio, and they also puts strong constrains on the time resolution of the electronic detection chain necessary to statistically recognize particles coming from the rear face of the telescopes.
In the present study, we extend our methodological work by quantitatively examining how the combination of muon tomography and gravity data may improve the reconstruction of the density distribution inside geological bodies. Studies combining muon data and gravity measurements remain scarce, and we emphasize the early study by Caffau et al. (1997) who compared muon tomography with gravity measurements. More recently, Davis & Oldenburg (2012) and Nishiyama et al. (2014) presented joined inversions of gravity data and muon tomography using a straightforward linear regularized inversion based on block models. Combining density radiographies and gravity measurements is a quite natural intention since both methods provide information directly related to the density distribution. However, both methods sample the density structure in very different ways, and the scanner-like principle of muon tomography confers to this method a high-resolution that could question the interest to perform a costly complementary gravimetry survey. Conversely, muon tomography has limitations due to difficulties to install muon telescopes at a large number of places around the geological body to image. Consequently the number of opacity radiographies available to perform the tomography reconstruction is always small with, eventually, wide uncovered angular sectors that make the inverse Radon transform ill-posed to perform the tomography reconstruction. Another strong limitation of muon tomography comes from the fact that it at best only samples the density structures located above the horizontal plane passing through the telescope location.
Considering these issues, gravity data may bring useful information to better constrain the recovery of the density distribution, and it is the purpose of the present study to quantitatively document the benefits of a joint inversion of gravity and muon data to recover the density structure. In order to obtain results independent of any particular parametrization (e.g. block discretization), we adopt a formulation in terms of resolving kernels whose expression only depend on the geometrical properties of the data acquisition (i.e. locations of measurement points and telescope acceptance functions). We begin by establishing the relationships between the density structure and both muon tomography and gravity data. Next, we derive the resolving kernels respectively corresponding to individual gravity and muon inversions and to joint gravity-muon inversion. The resolving kernels translate the information contained in the data into information concerning the density structure. We conclude the article with examples taken from real field experiments conducted on La Soufrière of Guadeloupe to illustrate the practical interest to combine muon radiographies and gravity measurements. The muon tomographies experiments were already described in various articles (Lespare et al., 2012 and Jourde et al. 2013). Three sites called Ravine Sud, Rocher Fendu and Savaneà Mulets were explored and are represented on Fig. 1(a). The gravimetry survey is currently running, for the purpose of this article we simulated one hundred measurements regularly spaced on a grid that covers the dome ( Fig. 1(b)).

The sampling of the density distribution by muon tomography and gravimetry
Here, we recall the main formula relating the density distribution to the data, i.e. fluxes of muons and gravity measurements. In the inverse problem framework, these formula describe the forward problem for each method. In the remaining, we suppose that the muon data have been cleaned from perturbing effects such as upward going fluxes of particles as described in

Muon tomography
The primary information used in muon tomography consists in cosmic muons flux attenuation measurements resulting from the screening produced by the geological volume to scan. Attenuation is measured by counting the number of muons emerging from the volume for each observation axis: s m = (r m , P m (ϕ, θ)), of the telescope ( Fig. 1(a)). r m represents the position of the telescope, P m the observation axis acceptance pattern which depends on (ϕ, θ)  the azimuth and zenith angles referenced at r m (see Fig. 2). Note that r m is the same for all the observation axes on a given site. P m depends on the telescope geometry and angular orientation on the site. Our standard field telescopes count 31 × 31 observation axes and, in a field experiment where the telescope successively occupies several places around the target, the number, M , of data may easily reach several hundredths. For example if we use the muon tomography data from the 3 Soufrière sites M = 3 × 31 × 31. In practice M is lower as many axes point downward or above the volcano. P m (cm 2 .sr.rad −2 ) shape depends on the detection matrices structure (see Lesparre et al. 2012a,b for further details). It has a steep peak centered on a small solid angle region Ω m . It is identically null outside Ω m ( Fig. 2 and Fig. 3). Observe that P m must not be confounded with the integrated pixel acceptance T m (cm 2 .sr) used for instance in Lesparre et al. (2010), The number of muons attributed to a given line of sight actually corresponds to all muons detected in Ω m . Inside the geological volume the trajectories of these muons describe a conical volume whose apex is located at the telescope, r m .
The attenuation of the muon flux caused by the rock screen depends on the amount of matter encountered by the particles along their trajectories. For a given straight trajectory t = (r, ϕ, θ) (r is a telescope site and (ϕ, θ) the azimuth and zenith angles referenced at r) it is quantified by the density line integral along t, the opacity where ρ is the density, L the particle path length, andρ is the average density along t. The differential flux associated with t may be expressed as a function δφ t = R A R T Figure 2. Muon tomography reference frame and notations. R A and R T respectively are the absolute orthonormal and the instrument reference frames. An observation axis s m = (r m , P m (ϕ, θ)) is represented with r m the vector that localizes the telescope position and P m its acceptance pattern (restrained to the solid angle Ω m ). The spherical coordinates (ϕ, θ) are here localizing the steep acceptance peak mentioned in section 2.1.
Note that the integration is restricted to the small solid angle Ω m because of the compact support of the observation axis acceptance P m (Fig. 3). δφ t is not linearly related to , however for small opacity fluctuations we assume that δφ t may be approximated by its first order development around the local average density, ρ 0 (r). ρ 0 (r) is the prior density model of the geological structure. For a given path t it reads, where 0 = t ρ 0 (ξ) dξ. Rearranging the terms and letting α t = dδφ t d ( 0 ), we obtain, Inserting eq. 5 into eq. 3 we get the approximate equation where φ 0 = φ m (ρ 0 ) is the flux corresponding to the prior density model ρ 0 (r) and t = (r m , ϕ, θ). In the remaining of the present paper, we shall use the centred and normalized flux, whereρ min andρ max are expected extreme values of the density.

Gravimetry
Gravimetry aims to estimate the gravity field generated by surrounding objects measuring locally the vertical acceleration they produce. The vertical acceleration g is directly related to the density spatial distribution through the Newton law: where the vector r n represents the location of the n th measurement point (in our example n runs from 1 to 100). As for muon tomography we use the normalized gravity anomalyg n defined asg 3. Resolving kernel approach

The acquisition kernels
We define X, the space that contains the set of continuous L 2 functions going from R 3 into R.
The 3D density distribution ρ belongs to X and it is related to the muon flux measurements, φ m , and to the gravity data,g n , through the action of acquisition kernels G and M which also belong to X. This reads, where ·, · X is X inner scalar product, and M and N are respectively the number of muon tomography and gravimetry data. From eq. (6) (7) (8) and (9) we obtain explicit expressions for M and G,  Observe that the 1/ξ 2 term in eq. (12) comes from the spherical coordinates elementary volume expression, ξ 2 dξdΩ, inserted in eq. (6). Examples of acquisition kernels are plotted in Fig. 4 and discussed in sections 4.1 and 4.2. The X structure allows to introduce prior information into the problem. For instance, the classical inner product of L 2 continuous functions, can be replaced by the weighted inner product, where the weight function w (w(r , r ) > 0, w(r , r ) = w(r , r )), plays the role of a covariance function that may be used to neglect the impact of the free air zone around the studied structure for gravimetry and muon tomography (see eq. (25) and its comment below). It may also serve to introduce a correlation length for the density variations.

The resolving kernel
The 3D density distribution,ρ(r), obtained by solving the set of linear equations (10,11) is a degraded version of the true density distribution, ρ(r), both because of the limited number of data available and of the filtering (i.e. blurring) effect of the acquisition kernels. In the remaining, we shall use the set of undifferentiated acquisition kernels and the set of undifferentiated normalized data, We now formulate the inverse problem in the framework of functional spaces where the family of acquisition kernels constitutes a set of generating functions of a subspace X K of X of dimension K (Tarantola & Nercessian, 1984;Bertero et al., 1985). This implicitly assumes that the ζ k are linearly independent with respect to the retained inner product, i.e. no acquisition kernel can be written as a linear combination of the other kernels. The noticeable instances where this important assumption is not satisfied correspond to situations where several data have been acquired identically, i.e. either at the same location for gravity measurements or with the same position and orientation of the telescope for muon tomography. In such cases the dimension of X K is reduced since the redundant data may be merged (i.e. averaged) into a single one.
The best density distribution that can be recovered through the inversion process (it is the best because it takes all the information contained in the data and makes the less hypotheses about X K complementary subspace) is a linear combination of the generating functions, The components a k of eq. (18) are obtained by minimizing the quadratic distance Y between the data and the corresponding values given by the density model, Y is the weighted Euclidean space that contains the measurements. The weights W k permits to introduce prior information about the measurements quality. It is possible to introduce crossed terms W ij if the measurements are not independent, but it is not the case here. We get where S k,j is the (k, j) component of the Gram matrix inverse defined as S k,j = W j × ζ k , ζ j X . Using eq. (20), eq. (18) becomes, The presence of ζ j , ρ X in the right hand part of this equation indicates that the density distribution actually recovered,ρ(r), is assembled from projections of the true unknown density, ρ(r), onto the acquisition kernels. The recovered density is a filtered version of the true density distribution, and the filter (i.e. the resolving kernel) depends on the data. This can be made more explicit by rewriting eq. (21) as (Bertero et al., 1985), where we introduce the resolving kernel, ζ j is the acquisition kernel ζ j modulated by the prior information represented by the w function. For instance, w may be an indicator function used to limit the support of the acquisition kernels to the volume of interest.

Characterisation of the resolving kernel
A resolving kernel, ∆(r, r ), is a function defined in the whole space that plays the role of a spatial filter. When applied to the true density distribution, it gives the reconstructed density.
The amplitude and the shape of ∆ render the achievable resolution of the reconstructed density structure. According to eq. (23) it is a linear combination of the acquisition kernels. ∆ may be characterized in different ways by using several properties to quantify its resolution and anisotropy. These properties should be the least possible dependant to a specific resolving kernel and allow the user to easily appreciate the resolution and its eventual bias. In the present study, we simply compare the resolving kernels against the ideal kernel represented by a Dirac distribution δ(r − r ). This is achieved through the projection γ of ∆(r, r ) onto a δ(r − r ), where the X scalar product is defined by eq. (15). Here, the Hamming function is used as the weight function where K w is a normalising constant, H L is a rectangular pulse that restricts w to r − r ∈ [−L/2; L/2], and L = 25 m. If r or r is located outside the volcano we take w(r, r ) = 0. This choice is explained in section 4.4. We now display resolving kernels corresponding to the data acquisition shown in Fig. 1(a) for muon tomography and in Fig. 1(b) for gravimetry. Muon radiographies are taken from three sites equidistantly located along the Southern edge of the volcano. Gravity measurements are assumed to be done on a regular grid over the entire lava dome.
Accounting for the fact that the acquisition kernels ζ k are either for gravity or for muon tomography (eq. (16)), we successively consider the case of resolving kernels obtained for muon tomography alone, for gravity data alone, and for a combination of muon tomography and gravity data. We compute ∆(r, r ) for two positions r = {r 1 ; r 2 } located along a vertical line that goes through the center of the dome (Fig. 5). Points r 1 and r 2 are respectively inside and below the volume of the lava dome spanned by the lines of sight of the telescopes ( Fig. 1(a)). The parameter γ is computed and plotted on the four characteristic slices represented on Fig. 5.   Fig. 1(a)) and point 2 is located at Z2 below the ray coverage of the telescope. Fig. 4(b) shows a gravimetry acquisition kernel G. Remind that the data used in the present study are normalized relatively to a reference model with density ρ 0 (eq. (9)). The gravimetry acquisition kernels are very sensitive to density fluctuations close to the measurement point because of its 1/r 2 term. The gravity data are actually the component of the gravity field anomaly taken along the local vertical, and the acquisition kernel becomes less and less sensitive as we get closer to the horizontal plane that contains the measurement point. The gravimetry inverse problem is systematically ill-posed (e.g. Al-Chalabi, 1971) because no matter the number of measurements the resolving kernel mostly integrates information around the measurement positions, i.e. near the surface. An illustration of this problem is given for the resolving kernels of r 1 and r 2 (Fig. 6(a),7(a)). For gravimetry inversions it is more realistic to model ρ(r) by a function that depends on a few discrete parameters (even if it means losing the linearity between the data and the measurements) rather than trying a continuous inversion.

Gravimetry kernels
Observe that G integrates the density over the entire volume and provide information for point r 2 located below the lines of sight of the telescope. Fig. 4(b) shows a typical muon tomography acquisition kernel M (eq. (12)). It has a conical shape whose aperture angle depends on the distance between the front and the rear detection matrices of the telescope. The apex of the kernel is located at the telescope, and the kernel widens as we move away from the telescope thus the local sensitivity is decreasing. Moreover the triangular shape of the intra-pixel acceptance P m (see Fig. 2 and Fig. 3) makes the sensitivity maximum along the main line of sight s m .

Muon tomography kernels
The Fig. 4(b) shows we are as sensitive to a density change occurring on a few tenth of meters in front of the telescope as to the same change happening on a few hundred of meters beside the volcano. It reveals how deterministic is the telescope position. If one desires to image or monitor a specific region belonging to a bigger structure, the measurement will be much more sensitive if the telescope is in front of it. The important heterogeneities inside the muon tomography acquisition kernels forbid us to use the Radon transform mathematical corpus. For an equivalent resolution and scanning the kernels can be regularized taking the telescope away from the volcano and reducing the angular aperture. We then get into the typical experimental conditions of a medical X-ray tomography. But the consequences are a weaker particle flux (a longer acquisition time) and a greater sensitivity to potential noises. So a compromise has to be found, but the actual lack of understanding of the noises and the already very long acquisition times we are facing lead us to take the telescope the closest we can to the volcano.
We draw the reader's attention to the fact that, despite their compact support, the acquisition kernels overlap each others for neighbour main lines of sight (see Fig. 3). As will be seen below, this characteristic is fundamental to understand the shape of the resolving kernels.
A muon tomography resolving kernel is a linear combination of muon tomography acquisition kernels M (eq. (23)), and the inversion process optimizes the b i coefficients to obtain the best density model (eq. (24)). Fig. 6(b) and Fig. 6(c) show the resolving kernel for point r 1 for different combinations of tomography datasets.
When using data acquired from a single place located at the Southernmost edge of the volcano (Ravine Sud, see Fig. 5(a)), the resolving kernel ( Fig. 6(b)) encompasses lines of sight spanning a limited range of azimuths. Consequently, the filtering effect of ∆(r 1 , r ) integrates ρ along a long narrow cone to give the estimated densityρ(r 1 ). The fact that ∆(r 1 , r ) has not a compact support like the M kernels comes from the overlapping of neighbour acquisition kernels that produces a transfer of information among lines of sight.
When simultaneously using all three muon tomography sites, the resolving kernel includes acquisition kernels that span a wider range of sight azimuths. Consequently, the resolving kernel is more localized onto point r 1 (Fig. 6(c)). However, the number of radiographies remains small and the kernel has a spider shape visible in the horizontal slice at the far right of Fig. 6(c).
Observe that the resolving kernel ∆(r 2 , r) = 0 since all acquisition kernels are null in this part of the volcano.

Joined muon tomography and gravimetry kernels
We now consider resolving kernels computed by using both muon M and gravity G acquisition kernels.
Gravimetry does not improve significantly the inversion process at r 1 , and the resolving kernel ∆(r 1 , r ) (Fig. 6(d)) looks very similar to the one obtained for the muon radiographies alone (top right of Fig. 6(c)). The information provided by muon tomography is dominant relative to gravimetry excepted at the immediate vicinity of the gravity measurement points.
The situation is very different for point r 2 where the resolving kernel ( Fig. 7(b)) obtained by joining muon radiographies and gravity data appear very different from the gravity kernel ( Fig. 7(a)). The most conspicuous effect is that muon data compensates the great sensitivity of gravimetry at near-surface locations by shifting the center of mass of the resolving kernel downward. This considerably improves the vertical resolution achievable in the deepest parts of the volcano.
The conclusions are different if only one tomography acquisition is available. In that case the gravimetry measurements have an impact on the upper part of the dome because they contribute to resolve the ambiguity about the anomaly spatial depth relatively to the acquisition position. But then the zone below the dome will lack data to be properly constrained.

Impact of prior information
The choice made for the X and Y weight functions w(r , r ) (eq. 15) and W i=1...K (eq. (19)) has an important influence on the obtained resolving kernels ∆(r, r ). (d) gravimetry and tomography data Resolving kernel value normalized at point 1 in logarithmic scale (e) Figure 6. Resolving kernel at point r 1 (Fig. 5) obtained for the gravity data alone, (a), the Ravine Sud muon tomography alone, (b), all three muon radiographies , (c), and the joined muon and gravity datasets (d). See Fig. 1(b) for the locations of gravity measurements and the three sites for muon radiographies. The resolving kernel absolute value is normalized with reference to the value computed at point r 1 and represented with a log 10 scale.
In the X space, the diagonal term w(r , r = r ) permits to adjust the local degree of prior knowledge on ρ(r). For instance, w(r, r) = 0 in regions where ρ(r) is assumed sufficiently well-known to have no impact on our measurements. This corresponds to situations where ρ 0 (r) = ρ(r) and where the concerned regions have not to be accounted for in the inversion process. In our case we use it to cancel the free-air impact on muon tomography and gravimetry, but we can also constrain it to incorporate direct field measurements of the density.
The non-diagonal part w(r , r = r ) may be used to introduce a spatial correlation in ρ. This can be done throughζ which the convolution of ζ with w (eq. (25)). Here, we use a simple Hamming function with a 25 m correlation length everywhere in the dome (eq. (27)), andζ is a smoothed version of ζ which attenuates the 1/r 2 effect previously mentioned (for muon tomography it permits to get closer to the X-ray tomography experimental conditions previously detailed). The correlation introduced by the Hamming function increases the acquisition kernels sensitivity further from the measurement point toward the central and the Northern parts of the dome. This produces a better localization of the resolving kernel at r 2 as can be checked by comparing Fig. 7(c) with Fig. 7(b) where no spatial correlation was applied. The counterpart of this effect is a de-sharpening of the kernel at point r 2 . w is a regularizing low-pass filter that removes spurious short-wavelength fluctuations in the density (c) gravimetry and tomography data + prior information Resolving kernel value normalized at point 2 in logarithmic scale (d) Figure 7. Resolving kernel at point r 2 (Fig. 5) obtained for the gravity data alone (a) and for the joined muon and gravity datasets (b). Figure (c) is obtained with both muon and gravity datasets and some prior information about the density spatial correlation. See Fig. 1(b) for the locations of gravity measurements and the three sites for muon radiographies. The resolving kernel absolute value is normalized with reference to the value computed at point r 2 and represented with a log 10 scale. model and reduces the ill-conditioning of the inverse problem (e.g. Bertero et al., 1988).
The choice of w is problem-dependent and must be sustained by prior knowledge. The Hamming function acts as a low pass filter with a limited support compatible with the large homogeneous zones observed on the field: massive andesite, hydrothermally altered material and possibly large cavities.
In the Y space, the weights W i=1...K allow to assign different quality factors to the available data at one inversion location. For instance, in muon tomography, the W s permit to account for the fact that all observation axes have not the same integrated acceptance T m (Fig. 3). The quality of the gravity data strongly depends on the ground stability (i.e. tilt stability during measurement sequences) and the presence of wind (i.e. vibrations of the gravity-meter). The non-diagonal terms W ij , i =j are null as the measurements are independent the ones from the others. (b) gravimetry and tomography data

γ maps
The γ(r) index defined in eq. (26) may be used to estimate the resolution achievable everywhere in the volcano. Fig. 8 shows slices of the γ function obtained for the gravity data (left part of Fig. 8) and by joining the three muon tomography data sets together with the gravity measurements ( Fig. 8(b)). The gravimetry γ slices clearly reveal the important sensitivity of the data to density variations located in the immediate vicinity of the measurement points and the very low sensitivity to density structure located deeper in the lava dome.
The muon-gravimetry γ slices confirm the results obtained for points {r 1 ; r 2 } and show the considerable improvement of the resolution obtained when jointly using the muon and gravity datasets. They also reproduce the asymmetric resolution due to the conical shape of the muon acquisition kernels M. Since the places occupied by the telescope are located along the Southern edge of the volcano, a finer resolution is obtained for the Southern part of the lava dome. This corresponds to the dark-red circular sector visible in the upper horizontal slice in Fig. 8(b). As expected, the resolution is coarser in the central and Northern parts of the dome. The same slice shows that the North-Eastern region of the volcano is resolved by the gravity data alone since no lines of sight of the telescope cross this part of the dome.
The γ map is a useful tool to plan an acquisition survey. One can easily compute how γ is changed with different possible measurement campaigns and select the most pertinent one depending on the region of interest, the available time on the field and the accessibility of the site. This choice is critical as muon tomography acquisitions are long (a few weeks) and gravity measurements delicate. It can also be used to design a mesh for the problem. The meshing elements density can roughly follows the γ map fluctuations.
We emphasize that other definitions may be used for γ and that a single index may prove insufficient to characterize the shape of resolving kernels. Consequently, we recommend to perform a 3D examination of individual resolving kernels at locations of particular importance (i.e. like detecting places where density changes occur).

Conclusion
The resolving kernel analysis discussed in the present study allows to quantitatively assess the way by which gravity data and density muon radiographies may be joined to reconstruct the density distribution inside geological bodies. A main result concerns the improvement of the resolution obtained in the deep regions of the density model when joining muon and gravity data, and despite the fact that these regions are not sampled by muon tomography. Part of the information brought by the muon data is transferred to the deep regions of the model through the long-range coupling of the gravity acquisition kernels ( Fig. 4(a)) used to construct the joined resolving kernels (eq. (23)). The compact support of the muon acquisition kernels ( Fig. 4(b)) allows muon tomography to obtain high-resolution models of the density structure of La Soufrière but restricted to the upper part of the dome. Gravimetry cannot be used on its own to make a proper continuous inversion because it is too sensitive to the very close subsurface. Therefore coupling the two techniques does not lead to significant improvement of the density model in regions where muon data are available (compare Fig. 6(c) and Fig. 6(d)).
The muon tomography acquisition kernel has a conic shape. and noise considerations obligate us to put the cone apex just in front the studied structure. It results in a large decrease of the sensitivity between the front and the rear of the volcano. This problem adds to the heterogeneous tomography sampling and forbid us to use the standard Radon transform usually adopted in X-ray tomography medical experiments to inverse the density. It also shows how deterministic is the telescope position if one desires to image or monitor a specific region belonging to a bigger body.
The positive weight function w of the inner product (eq. (15)) may be used to introduce prior information both by limiting the support of the density distribution to reconstruct and by introducing a spatial correlation smoothing the 1/r 2 effect of muon tomography and gravimetry acquisition kernels. This extends the range of sensitivity of the measurements and results in a solution with a more homogeneous quality.
Finally the γ maps give an overview of the resolving kernel geometry everywhere in the dome which can easily be used to optimally plan future acquisition surveys with respect to data already available.