Читать книгу Muography - Группа авторов - Страница 50
3.3 LINEAR JOINT INVERSION
ОглавлениеThe formulation of the joint inversion (gravity and muography) has been provided by several researchers. Here, let us describe one simple (linear) formulation by Davis & Oldenburg (2012) and Nishiyama et al. (2014). Note that this formulation is not the only one.
Suppose the target volume of interest is composed of small rectangular prisms (Fig. 3.4) with density values of ρ j (j = 1, 2, ⋯, n). The vertical component of the gravitational effect produced at the i‐th gravity station can be written as
Figure 3.4 Schematic illustration of linear joint inversion of muography and gravity data.
where G ij is the gravitational contribution of the j‐th voxel to the i‐th gravity station for unit density. Given a geometry of the prism as , the gravity kernel G ij is expressed as
where G = 6.674 × 10−11 m3/kg/s2 is the universal gravitational constant. The analytical formula of the integration in equation 3.7 is introduced in the Supplemental Information. On the other hand, the density‐length X i , derived from muography analysis along the i‐th line of sight, is approximated in the same way with equation 3.6:
where L ij represents the length of the i‐th trajectory confined in the j‐th prism. Since the gravity anomaly and density‐length are both written as a linear combination of unknown densities, concatenating the data vectors and design matrices
leads to the formulation of the linear inverse problem
This simple formulation was first provided by Davis & Oldenburg (2012). Here, the data vector d is a column vector with n muon + n grav elements, where n muon is the number of the muography rays and n grav is the number of the gravity stations. The design matrix A has (n muon + n grav) ×n elements.
One obvious solution to equation 3.10 could be obtained by multiplying the inverse matrix of A from the left. However, such an inverse matrix does not exist because n muon + n grav ≠ n in general. Even when n muon + n grav = n coincidentally, A is rank deficient in most cases and A –1 d does not provide a realistic solution; it tends to place a huge density anomaly in the vicinity of muography detectors or gravity stations. Nishiyama et al. (2014) propose a solution to circumvent the rank‐deficit problem by employing a Bayesian probabilistic approach (see Tarantola, 2005, for details). In this approach, the information we have on the true value of data d is described by a Gaussian probability density function (pdf) with its peak at the observed value d obs in a (n muon + n grav)‐dimensional space. The observation error around the peak is described by a covariance matrix, C d . Besides, the intuition we have on the density distribution, that the density values should be around a certain value, can be introduced as a Gaussian pdf with its peak ρ 0 and covariance matrix C ρ (prior pdf). Bayes theorem then convolutes the two pdfs and provides an updated pdf on ρ , so‐called posterior pdf. As long as the data and prior pdfs are Gaussian, the posterior pdf also takes a Gaussian form and its peak ρ ' and covariance C ρ ' are given as
(3.11)
and
(3.12)
This solution is not robust enough yet. It varies greatly depending on how the target volume is gridded. One strategy is to introduce a smoothing constraint on the prior covariance matrix C ρ , such as
Figure 3.5 Three‐dimensional representation of density distribution inside Showa‐Shinzan lava dome. The red points and blue lines indicate the gravity stations and the lines of sight of the muography detector, respectively.
Redrawing Figure 6 of Nishiyama et al. (2017).
where σ ρ is the allowance of deviation from the initial guess density, d(i,j) is the distance between the i‐th and j‐th prisms, and λ is the correlation length, which controls the correlation of neighboring prisms. When λ is large, the smoothing effect propagates to a longer distance and the inversion tends to be over‐determined, whereas it becomes ill‐posed when λ is small.
Nishiyama et al. (2017) applied the above linear inversion to a lava dome located in Usu volcanic region, Japan (Showa‐Shinzan lava dome). The surveyed target is characterized by an uplifted plateau with a diameter of ~ 1 km and a dacitic lava dome with a diameter of 200 m on the plateau. This target had been once surveyed by muography in the pioneering era of muography (Tanaka et al., 2007). The first muography shows a two‐dimensional density map that suggests that the lava and surrounding plateau have distinct contrast between high and low density. Additional muography and gravity surveys were then performed to constrain the exact shape of the underground lava block (Nishiyama et al., 2017). Gravity surveys were performed at thirty stations on/around the dome and the muography observation was performed from 500 m west of the summit using emulsion films. Fig. 3.5 shows the obtained three‐dimensional density distribution, which was recovered from one muography and thirty gravity data. The result indicates the high‐density lava block extends vertically in a cylindrical shape within 200 m depth from the surface.
Cosburn et al. (2019) applied a similar linear joint inversion to muography data taken from an underground tunnel and gravity data on the surface and inside the tunnel in Los Alamos canyon, where a regional stratigraphy is very well‐studied. Above the tunnel, there is a layer of high‐density beds (deposits from volcanic surge) spreading horizontally. This configuration provides a unique opportunity to verify the resolution of the joint inversion and to tune the smoothing methods and related parameters. To properly reconstruct such a layered structure in the target volume, they introduce an anisotropic function with two correlation lengths, λ xy and λ z , to allow for independent depth and lateral variations. Specifically, instead of a single correlation length (equation 3.13), they introduce
(3.14)
where x ij , y ij , and z ij are horizontal distances and vertical distance of i‐th and j‐th prisms, respectively. By assigning smaller values for vertical correlation length (λ z < λ xy ), stronger constraints are imposed on prisms at the same elevation, and then a layered structure is preferably reconstructed.
