surface. This term can be evaluated by assigning the optimal density to the topographic material and integrating the gravitational effect (see Supplemental Information). The third term, Δg trend, is a so‐called regional trend, which originates from the deeper parts of the Earth, such as lithology, tectonics, and mantle dynamics. The latter two terms must be properly understood and removed to extract the first term, the gravitational effect of interest.
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
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)
