the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
Feasibility of threedimensional density tomography using dozens of muon radiographies and filtered back projection for volcanos
Shogo Nagahara
Seigo Miyamoto
This study is the first trial to apply the method of filtered back projection (FBP) to reconstruct threedimensional (3D) bulk density images via cosmicray muons. We also simulated threedimensional reconstruction image with dozens of muon radiographies for a volcano using the FBP method and evaluated its practicality.
The FBP method is widely used in Xray and CT image reconstruction but has not been used in the field of muon radiography. One of the merits of using the FBP method instead of the ordinary inversion method is that it does not require an initial model, while ordinary inversion analysis needs an initial model.
We also added new approximation factors by using data on mountain topography in existing formulas to successfully reduce systematic reconstruction errors. From a volcanic perspective, lidar is commonly used to measure and analyze mountain topography.
We tested the performance and applicability to a model of Omuroyama, a monogenetic scoria cone located in Shizuoka, Japan. As a result, it was revealed that the density difference between the original and reconstructed images depended on the number of observation points and the accidental error caused by muon statistics depended on the multiplication of total effective area and exposure period.
Combining all of the above, we established how to evaluate an observation plan for volcanos using dozens of muon radiographies.
1.1 Muon radiography and its principles
Muon radiography is a method that can be used to make a map of the inner bulk density structures of large objects such as volcanoes, archeological targets, and so on, using secondary cosmicray muons. These muons are generated by the interactions between highenergy primary cosmic rays (the main component is proton) and nuclei in the atmosphere. The flux, energy spectrum, and the zenith angle dependence of secondary cosmicray muons have been well researched (e.g., Dorman, 2004; Honda et al., 2004; Patrignani et al., 2016; Nishiyama et al., 2016). Their behavior including energy loss in various materials has also been investigated (Groom et al., 2001). Therefore, when we assume “density length”, which is the integral of multiplication of density and material thickness, we can evaluate the number of penetrating muons. Muon detectors have been developed in the field of particle physics and cosmicray physics. To make a bulk density map, we need to measure not only the counts of penetrating muons from the target, but also the direction. For example, nuclear emulsion films (Morishima et al., 2017), hodoscope by scintillating plastic bars (Jourde et al., 2013, Ambrosino et al., 2015), glass resistive plate chambers (Ambrosino et al., 2015), the multiwire proportional chambers (Oláh et al., 2018) are capable of doing that. By implementing these muon detectors around the target, we can obtain the penetrating muon flux for each direction from the detector; then by comparing to the initial muon flux, we also obtain the attenuation of muons for each directions. By using the topographic data of the target, it is possible to derive the twodimensional averaged bulk density from the muon attenuation and the path length of the target material.
The principle of Xray radiography and muon radiography is very similar. There are two significant differences between these two methods: the first is the attenuation length. Typical Xray beams can penetrate the material in less than 1 m water equivalent. Conversely, some muons can penetrate 1 km of water equivalent because their kinetic energy is very high. The second difference is the origin of the source. The source of cosmicray muons is completely environmental and we cannot control the flux while Xray beams are generated by accelerating the electron artificially. Typically, the number of observed muons is much smaller than in ordinary Xray radiography.
The first significant result for volcanology was the twodimensional bulk density imaging of the shallow conduit in Mt. Asama by Tanaka et al. (2007a). Several observations have been made since this research (e.g., Tanaka et al., 2007b, 2014; Lesparre et al., 2012).
1.2 Threedimensional bulk density imaging
The internal structure of volcanoes gives important information for volcanology. For example, the shape of a shallow conduit affects the eruption dynamics (Ida, 2007). However, muon radiography by only one direction makes just a 2D image, and this density is the average of material along the muon path direction. Therefore, if we find some contrast in the 2D density image, we cannot distinguish the actual position of this density anomaly along the muon path direction. To observe the real conduit shape, it is necessary to obtain the density image from different directions to reconstruct the threedimensional bulk density image.
Tanaka et al. (2010) attempted to observe the target from two directions in Mt. Asama. Nishiyama et al. (2014, 2017) conducted a 3D density analysis in the ShowaShinzan lava dome, combined with gravity observation data, which is also sensitive to density. Jourde et al. (2015) evaluated this jointinversion method between muon radiography and gravity, and they observed and conducted 3D density analyses by using threepoint muon radiography and gravity data (RosasCarbajal et al., 2017). These previous studies required prior information on internal density distribution because of insufficient observation data, and they were performed using the inversion technique.
In this study, we propose the application of a 3D density reconstruction analysis method using filtered back projection (FBP), which does not require prior information. In the FBP method, it is possible to obtain a 3D density distribution from many projection images only without using the inversion technique (Deans, 2007). This method is applied to Xraycomputed tomography (CT), so it is possible to apply to muon radiography data in principle. However, muon radiography differs from Xray CT in three points. First, there is a constraint on the number of observation points and position. In Xray CT, there are hundreds of observation points, and each position is controllable. However, for muon radiography, we can only use several dozen points, and the positions are limited because of topography. Second, the cosmicray muon attenuation flux is not a simple exponential. Therefore, the influence of muon statistical error depends on the results of 3D density, which is not trivial. Third, in the case of muon radiography, typically the amount of signal is much less than for Xray because the source of cosmicray muons is completely environmental. Therefore, it is important to study the features of the FBP method in the case of realistic observations with various numbers of muon radiographies. So we should consider not only the reconstruction error of the FBP method, but also how the error of muon statistics propagates to the final image.
The radon transform is used to obtain projection images from all directions with respect to a density distribution. In muon radiography, this corresponds to acquiring observation data on density length from all directions. For three dimensions, the radon transform $p(X,Z,\mathit{\beta})$ of an object with density $\mathit{\rho}(x,y,z)$ is given by the following:
where x, y and z are the positions in a 3D volume; X and Z are the tangents of azimuth and elevation angle values, respectively; β is the observation point position at a counterclockwise angle with respect to the y axis; and D is the distance between the observation point and the origin. Figure 1 shows the geometric definition for these parameters.
In a 3D case, if observation data are conebeam data and observation points only exist on the circumference, a complete inverse radon transform does not exist. Therefore, approximation is needed. Feldkamp (1984) proposed one of the best methods to approximate a solution with a small elevation angle in two dimensions. This approximation is written as follows:
where ${Z}_{\mathrm{0}}=z/(Dx\mathrm{sin}\mathit{\beta}y\mathrm{cos}\mathit{\beta})$, ${L}_{\mathrm{2}}=\sqrt{\mathrm{1}+{Z}_{\mathrm{0}}^{\mathrm{2}}}(D+x\mathrm{sin}\mathit{\beta}y\mathrm{cos}\mathit{\beta})$, ${X}_{\mathrm{0}}=(x\mathrm{cos}\mathit{\beta}+y\mathrm{sin}\mathit{\beta})/{L}_{\mathrm{2}}$ and h(X) are a Ram–Lak filter (Ramachandran and Lakshminarayanan, 1971). A feature of this method is that it does not require the shape or initial model of the object. However, when there is a density change in the vertical direction, the accuracy of the approximation decreases. In many examples of volcanic muon radiography, we obtain the shape of the volcano by using other methods; therefore, the influence of changes in the shape can improve the accuracy of the approximation. To estimate the elevation angle, we use the ratio of the path length of the observed muon $q\left(X,{Z}_{\mathrm{0}},\mathit{\beta}\right)$ to the approximation of ${q}_{h}\left(X,{Z}_{\mathrm{0}},\mathit{\beta}\right)$ (see Fig. 2), which can be written as follows:
where $p{}^{\prime}\left(X,Z,\mathit{\beta}\right)$ is the approximation of the density length for the inverse radon transform. Finally, the reconstruction calculation formula can be written as follows:
where m, n is the index of X, β, respectively. We name this the “path length normalization approximation”.
In this section, we describe the specific components of the simulation calculation. The simulation calculation is divided into the following four steps:

parameter setup

simulation calculation of the observed muon counts

reconstruction calculation using data created in step 2

calculations for evaluating the reconstruction results.
3.1 Parameter setup for target and detector
We simulated and reconstructed the density structure of Omuroyama, which is located in Shizuoka, Japan. We chose this volcano for two reasons. First, this volcano is easily observable from all directions because there are no large structures around the surrounding muon shields in a topographical view. Second, there are no occurrences of muon radiography for these large scoria hills. Omuroyama is a large scoria hill. We base the internal structural model of the large scoria hill on observations at the time of its formation (Luhr and Simkin, 1993). However, there are currently no direct examples of these observations.
Figure 3 shows the contour map of the Omuroyama model used in the simulation. We assume that the x axis is in the east–west direction, the y axis is in the north–south direction and the origin is the summit.
We configure the internal density distribution similarly to a checkerboard with a side length of 100 m and a density of 1000 and 2000 kg m^{−3}. We presume that the first internal density distribution is defined as the original image and is expressed as ${\mathit{\rho}}^{\text{ori}}(x,y,z)$.
The field of view was set from −2 to 2 (−63.4 to 63.4 in degrees) horizontally and 0 to 1 (0 to 45 in degrees) vertically, and the angular resolution was set to 0.04 (2.3 in degrees) in tangent. The observed muon statistics affect the density reconstruction error: the number of muons observed is proportional to the effective area of the device S and the exposure period T. The total effective area and exposure period ST of all muon devices was set as 1000 m^{2} ⋅ days. For example, when the number of observation points is 16, each ST per point is $\mathrm{1000}/\mathrm{16}=\mathrm{62.5}$ m^{2} ⋅ days.
All observation points were assumed to be on the circumference of radius D=500 m placed at the center $\left(x,y\right)=(\mathrm{50}$ m, 50 m) of the mountain. The position of the observation points on the circumference is equal to the rotation angle from the reference line. The position β(rad) of the observation point is defined counterclockwise from the straight line parallel to the y axis and passes through the center $\left(x,y\right)=(\mathrm{50}$ m, 50 m) of the mountain. The value of β, on which the observation point is placed, must always be 1 at β=0, with the rest arranged at equal intervals along the circumference. For example, for the 16observationpoint case, the position of the observation point is ${\mathit{\beta}}_{n}=\frac{\mathrm{2}\mathit{\pi}}{\mathrm{16}}n$ (n=0, 1,..., 15). Figure 3 also shows the observation point arrangement when there are 16 observation points.
3.2 Simulation calculation of muon count observation
The simulation calculation of the observed number of muons is mainly performed with the following procedures.

Calculate the density length $p(X,Z,\mathit{\beta})$ from ${\mathit{\rho}}^{\text{ori}}\left(x,y,z\right)$ for each observation direction viewed from the observation point.

Calculate the theoretical muon flux ${F}_{\mathrm{0}}(X,Z,\mathit{\beta})$ by using a previously prepared relationship among the muon flux, elevation angle and penetration density length. We used the cosmicray muon flux model of Honda et al. (2004) and the muon energy attenuation of Groom et al. (2001) for the calculations made here.

Calculate the theoretical muon count observation ${N}_{\mathrm{0}}(X,Z,\mathit{\beta})$ by multiplying ${F}_{\mathrm{0}}(X,Z,\mathit{\beta})$ the device area S of the observation period T and the solid angle of spatial decomposition in the observation direction.
Figure 4a shows the observation state at observation point A in Fig. 3, and Fig. 4b shows the theoretical muon count observation ${N}_{\mathrm{0}}(X,Z,{\mathit{\beta}}_{\mathrm{1}})$ at that time.
It is not suitable to use the muon flux table in the region of the 10 m water equivalent or less because of small change. To avoid this region, we did not use these data when the path length $q\left(X,{Z}_{\mathrm{0}},\mathit{\beta}\right)$ was 10 m or less.
3.3 Reconstruction calculation
The reconstruction calculation procedure is as follows:

Calculate the muon flux ${F}_{\mathrm{0}}(X,Z,\mathit{\beta})$ from the muon number ${N}_{\mathrm{0}}(X,Z,\mathit{\beta})$, device shape and observation period.

Calculate the observed density length $p(X,Z,\mathit{\beta})$ from ${F}_{\mathrm{0}}(X,Z,\mathit{\beta})$, as well as the relationship among the muon flux, elevation angle and penetration density length.

In the path length normalization approximation, calculate path length $q\left(X,{Z}_{\mathrm{0}},\mathit{\beta}\right)$ and approximation path length
${q}_{h}\left(X,{Z}_{\mathrm{0}},\mathit{\beta}\right)$ from the shape information.

Calculate the density reconstructed image ${\mathit{\rho}}^{\text{rec}}\left(x,y,z\right)$ from the density length $p(X,Z,\mathit{\beta})$ by using Eq. (4).
4.1 Systematic error evaluation
We evaluated the systematic error, which is defined as the density difference between the original and reconstructed images at two points. First, we compared the differences between the methods for approximating the elevation angle (i.e., Feldkamp approximation and path length normalization approximation). Second, we quantified the relationship between the observation points and systematic errors.
4.1.1 The relationship between the observation points and systematic errors
We modeled scenarios with 4, 8, 16, 32, 64, 128 and 256 observation points. The reconstruction results are shown in Fig. 5. The systematic error $\mathit{\delta}{\mathit{\rho}}^{\text{sys}}\left(x,y,z\right)$ was defined as $\mathit{\delta}{\mathit{\rho}}^{\text{sys}}\left(x,y,z\right)={\mathit{\rho}}^{\text{rec}}\left(x,y,z\right){\mathit{\rho}}^{\text{ori}}(x,y,z)$. To evaluate the systematic error of all the reconstruction results, we calculated the average of $\mathit{\delta}{\mathit{\rho}}^{\text{sys}}\left(x,y,z\right)$ over the entire object area as the average value of systematic error μ^{sys}, and the sample standard deviation $\mathit{\delta}{\mathit{\rho}}^{\text{sys}}\left(x,y,z\right)$ was defined as the systematic error distribution σ^{sys}.
The relationship between the number of observation points and systematic error deviation σ^{sys} for the entire mountain body is shown in Fig. 6. As the number of observation points increases, the systematic error decreases. At an angular resolution of 0.04, there is almost no change at 64 or more points. At a resolution of 0.02, there is no change with more than 128 points. Therefore, when paying attention to the method of approximating the elevation angle, there are a number of implications when the number of observation points is 64 or more.
4.1.2 Comparison of Feldkamp approximation and path length normalization approximation
We simulated both the Feldkamp approximation and path length normalization approximation. Figure 7 shows the reconstruction results of both approximations. In the Feldkamp approximation, the average value of the systematic error μ^{sys} was $\mathrm{0.22}\times {\mathrm{10}}^{\mathrm{3}}$ kg m^{−3}, whereas it was $\mathrm{0.01}\times {\mathrm{10}}^{\mathrm{3}}$ kg m^{−3} for the path length normalization approximation.
4.2 Evaluation of accidental errors
We also evaluated the accidental error $\mathit{\delta}{\mathit{\rho}}^{\text{acc}}\left(x,y,z\right)$ in the reconstruction results. This value is the density error caused by the statistical error of muon count $N(X,Z,{\mathit{\beta}}_{\mathrm{1}})$. We assumed that ${N}_{\mathrm{0}}(X,Z,\mathit{\beta})$ follows a Poisson distribution. We generated 500 types of values with errors assigned, according to the Poisson distribution (in the following, referred to as “muon statistical error”) of ${N}_{\mathrm{0}}(X,Z,\mathit{\beta})$. This is referred to as muon count with statistical error ${N}_{j}^{\text{stat}}(X,Z,\mathit{\beta})$ (j=1 to 500). Here, the index “j” represents the trial of different seeds of random numbers set to ${N}_{j}^{\text{stat}}(X,Z,\mathit{\beta})$ for every X,Z and β.
The accidental error $\mathit{\delta}{\mathit{\rho}}^{\text{acc}}\left(x,y,z\right)$ was defined as follows:
Figure 8a, b and c show the spatial distribution of the accidental errors. The accidental error did not depend on the location in the plane. The accidental error was smaller in a section with higher altitude, i.e., a section with a large elevation angle at observation. Moreover, we saw this trend regardless of the number of observation points.
We defined the average of $\mathit{\delta}{\mathit{\rho}}^{\text{acc}}\left(x,y,z\right)$ over the entire object area as the average systematic error value μ^{acc}, and the sample standard deviation of $\mathit{\delta}{\mathit{\rho}}^{\text{acc}}\left(x,y,z\right)$ was taken as the accidental error distribution σ^{acc}. Even if the number of observation points increased, no significant changes were observed in the accidental error.
5.1 Limit of systematic errors
In Fig. 6, the systematic error does not converge to zero even if the number of observation points increases to more than 200. The observation point position β is represented by a counterclockwise rotation (see Fig. 1 definition of parameters). The interval of β is the angular resolution of the observation point. Increasing the number of observation points is equivalent to increasing the angular resolution of β. When comparing the resolution of X with the resolution of β for the 64point observation, the resolution of β is $\mathrm{360}/\mathrm{64}=\mathrm{5.6}$^{∘}, the angular resolution is 2.3^{∘} and the resolution of β is lower than X. However, for 256 points, the angular resolution of β is 1.4^{∘}, which is higher than the angular resolution of X. Figure 6 shows that the systematic error converges near the number of observation points when the resolution of β exceeds the resolution of X. These results indicate that the systematic error depends on the poor resolutions of both X and β.
5.2 Influence of path length on approximation
Why is the average systematic error value different between the Feldkamp approximation and path length normalization approximation? For a volcano with a structure similar to Omuroyama, which is cone shaped with a crater on the summit, the length of the muon path and the elevation angle tend to be shorter than the path length estimated in the horizontal plane (see Fig. 2). In the path length normalization approximation, given that the approximation is made with the path length as a reference, the difference in path length is not important in the Feldkamp approximation; however, the difference in path length is not taken into consideration and is influenced by the change in the path length. As a result, in the Feldkamp approximation, the average value of the systematic error is negative because of the presence of results with short path lengths.
5.3 Elevation angle dependence of accidental errors
Why does the accidental error δρ^{acc} become smaller as the elevation section increases? The accidental error δρ^{acc} occurs as a result of muon statistical error. The muon statistical error follows a Poisson distribution. As the number of observed muons increases, the muon statistical error becomes relatively small. Conversely, the muon flux increases as the elevation angle increases. In a section with high altitude, the reconstruction calculation uses both data with a large elevation angle and data with a large number of observed muons, thus reducing the accidental error.
5.4 Relation between muon counts and accidental errors
We performed these simulations under the condition that the total effective area of the observation device is equal. For a 16point observation, ST per point is 62.5 m^{2} ⋅ days; for a 32point observation, the device area S per point is 2 times greater at 31.25 m^{2} ⋅ days. Nevertheless, the results for the final accidental error values did not depend on the number of observation points (see Table 1). In Eq. (4), the operator is $\sum _{n=\mathrm{1}}^{N}$, where N is the number of observation points, and the factor ${p}_{c}\left({X}_{m},{Z}_{\mathrm{0}n},{\mathit{\beta}}_{n}\right)$ corresponds to the number of observed muons. p doubles if N is divided by 2 because the effective area also doubles. As a result of the calculation, $\mathit{\rho}\left(x,y,z\right)$ in Eq. (4) remains the same for every x, y and z value (i.e., each voxel). This is why the accidental error is nearly identical between the fourpoint observation and 64point observation.
This discussion is able to apply for actual observation with any muon detector type. In the case of an emulsiontype detector, it is easy to divide the effective area S. In the case of hodoscopetype detectors, we can divide the exposure period T by moving the detector to another observation point (e.g., Tanaka, 2016).
5.5 Evaluation of observation plan
We summarized systematic error and accidental error for Omuroyama and ST = 1000 m^{2} ⋅ days in Table 1. We can consider the better conditions of observation from this table. In this table, systematic error is larger than accidental error, excluding the case of 64 points and the 450 m cross section. When the number of observation points is 4 to 32, ST = 1000 m^{2} ⋅ days is sufficient, but in the case of 64 points, it is better to use more STs.
5.6 Limit of this simulation
In this evaluation, the observation points were arranged in a circular orbit. In the future, it is necessary to study more realistic observation point placements. For example, it is difficult to put the observation points on the same plane or in the same interval of β because of topography. We should also work these cases as a next step.
We simulated the systematic error of the 3D density structure of Omuroyama volcano by using several muon detectors via the FBP method with and without information on mountain topography.
 (i)
Systematic error, which is defined as the density difference between the original and reconstructed images in each voxel internal mountain, depends on the angular resolution of the muon detectors and the number of observation points.
 (ii)
By comparing the systematic error with and without information on mountain topography, the systematic error deviations are nearly identical. However, the mean value of systematic error becomes more precise in the former case; i.e., the value is more precise when a new method of approximation of path length normalization is used.
In addition, we studied the propagation of muon statistics to the final reconstruction results. By assuming that the multiplication of total effective area and exposure period is fixed and by changing only the number of observation points, the accidental error caused by muon statistics does not change. This accidental error depends only on the total muon statistics for all observation points.
Considering the above, we established how to evaluate an observation plan of dozens of muon radiographies.
We assumed that there are tens observation points in this study. The actual observations, which involve many nuclear emulsion muon detectors, were executed by Morishima et al. (2017). Furthermore, Oláh et al. (2018) succeeded in developing a highquality and inexpensive multiwire proportional chamber system. Considering such recent advances, the CT volcanic observation of volcanoes by using numerous muon detectors will be possible in the near future.
The muon flux data we used have been published and are accessible through Honda et al. (2004) and Groom et al. (2001). Omuroyama DEM data can download from the Geospatial Information Authority of Japan.
SN led the research, wrote the paper, and made the simulation code. SM supervised the work of SN, contributed to discussions, and edited the paper.
The authors declare that they have no conflict of interest.
The authors appreciate Masato Koyama of Shizuoka University
and Yusuke Suzuki of Izu Peninsula Geopark Promotion Council for advice and help with the geological structure in Omuroyama. The support of The
University of Tokyo is acknowledged.
Edited by: Lev Eppelbaum
Reviewed by: two anonymous referees
Ambrosino, F., Anastasio,A., Bross, A., Béné, S., Boivin, P., Bonechi, L., Cârloganu, C., Ciaranfi, R., Cimmino, L., Combaret, Ch., Alessandro, R. D'., Durand, S., Fehr, F., Français, V., Garufi, F., Gailler, L., Labazuy, Ph., Laktineh, I., Lénat, J.F., Masone, V., Miallier, D., Mirabito, L., Morel, L., Mori, N., Niess, V., Noli, P., PlaDalmau, A., Portal, A., Rubinov, P., Saracino, G., Scarlini, E., Strolin, P. and Vulpescu, B. : Joint measurement of the atmospheric muon flux through the Puy de Dôme volcano with plastic scintillators and Resistive Plate Chambers detectors, J. Geophys. Res.Solid Earth, 120, 7290–7307, https://doi.org/10.1002/2015JB011969, 2015.
Deans, S, R.: The Radon Transform and Some of Its Applications, Courier Corporation, 2007.
Dorman, L. I.: Cosmic Rays in the Earth's Atmosphere and Underground, Kluwer Academic Publishers, Dordrecht, Boston, London, 2004.
Feldkamp, L. A., Davis, L. C., and Kress, J. W.: Practical conebeam algorithm, J. Opt. Soc. Am, A1, 612–619, 1984.
Groom, D. E., Mokhov, N. V., and Striganov, S. I.: Muon stopping power and range tables 10 MeV–100 TeV, Atom. Data Nuc. Data Tab., 78, 183–356, 2001.
Honda, M., Kajita, T., Kasahara, K., and Midorikawa, S.: New calculation of the atmospheric neutrino flux in a threedimensional scheme, Phys. Rev. D, 70, 043008, https://doi.org/10.1103/PhysRevD.70.043008, 2004.
Ida, Y.: Driving force of lateral permeable gas flow in magma and the criterion of explosive and effusive eruptions, J. Volcanol. Geotherm. Res., 162, 172–184, 2007.
Jourde, K., Gibert, D., Marteau, J., Jean de Bremond d'Ars, Gardien, S., Girerd, C., Ianigro, J. C., and Carbone, D.: Effects of upwardgoing cosmic muons on density radiography of volcanoes, Geophys. J. Int., 40, 6334–6339, https://doi.org/10.1002/2013GL058357, 2013.
Jourde, K., Gibert, D., and Marteau, J.: Improvement of density models of geological structures by fusion of gravity data and cosmic muon radiographies, Geosci. Instrum. Method. Data Syst., 4, 177–188, https://doi.org/10.5194/gi41772015, 2015.
Lesparre, N., Gibert, D., Marteau, J., Komorowski, J. C., Nicollin, F., and Coutant, O.: Density muon radiography of La Soufri`ere of Guadeloupe volcano: comparison with geological, electrical resistivity and gravity data, Geophys. J. Int., 190, 1008–1019, 2012.
Luhr, J. F. and Simkin, T.: Paricutin: The Volcano Born in a Mexican Cornfield, Geoscience Press Inc. 1993.
Morishima, K., Kuno, M., Nishio, A., Kitagawa, N., Manabe, Y., Moto, M., Takasaki, F., Fujii, H., Satoh, K., Kodama, H., Hayashi, K., Odaka, S., Procureur, S., Attie, D., Bouteille, S., Calvet, D., Filosa, C., Magnier, P., Mandjavidze, I., Riallot, M., Marini, B., Gable, P., Date, Y., Sugiura, M., Elshayeb, Y., Elnady, T., Ezzy, M., Guerriero, E., Steiger, V., Serikoff, N., Mouret, J., Charles, B., Helal, H., and Tayoubi, M., Discovery of a big void in Khufu's Pyramid by observation of cosmicray muons, Nature, 552, 386, 2017.
Nishiyama, R., Tanaka, Y., Okubo, S., Oshima, H., Tanaka, H. K. M., and Maekawa, T,: Integrated processing of muon radiography and gravity anomaly data toward the realization of highresolution 3D density structural analysis of volcanoes: Case study of ShowaShinzan lava dome, Usu, Japan, J. Geophys. Res.Solid Earth, 119, 699–710, https://doi.org/10.1002/2013JB010234, 2014.
Nishiyama, R., Taketa, A., Miyamoto, S., and Kasahara, K.: Monte Carlo simulation for background study of geophysical inspection with cosmicray muons, Geophys. J. Int., 206, 1039–1050, 2016.
Nishiyama, R., Miyamoto, S., Okubo, S., Oshima, H., and Maekawa, T.,: 3D Density Modeling with Gravity and MuonRadiographic Observations in ShowaShinzan Lava Dome, Usu, Japan, Pure Appl. Geophys., 174, 1061–1070, 2017.
Oláh, L., Tanaka, H. K. M., Ohminato, T., and Varga, D.: Highdefinition and lownoise muography of the Sakurajima volcano with gaseous tracking detectors, Sci. Rep., 8, 3207, https://doi.org/10.1038/s41598018214239, 2018.
Patrignani, C., et al. (Particle Data Group): The Review of Particle Physics (2017), Chin. Phys. C, 40, 100001 (2016) and 2017 update, 2016.
Ramachandran, G. N. and Lakshminarayanan, A. V.: Threedimensional Reconstruction from radiographs and electron Micrographs, P. Natl. Acad. Sci. USA, 68, 2236–2240, 1971.
RosasCarbajal, M., Jourde, K., Marteau, J., Deroussi, S., Komorowski, J.C., and Gibert, D.: Threedimensional density structure of La Soufrière de Guadeloupe lava dome from simultaneous muon radiographies and gravity data, Geophys. Res. Lett., 44, 6743–6751, https://doi.org/10.1002/2017GL074285, 2017.
Tanaka, H. K. M.: Instant snapshot of the internal structure of Unzen lava dome, Japan with airborne muography, Sci. Rep., 6, 39741, https://doi.org/10.1038/srep39741, 2016.
Tanaka, H. K. M., Nakano, T., Takahashi, S., Yoshida, J., Takeo, M., Oikawa, J., Ohminato, T., Aoki, Y., Koyama, E., Tsuji, H., and Niwa, K.: High resolution imaging in the inhomogeneous crust with cosmicray muon radiography: The density structure below the volcanic crater floor of Mt. Asama, Japan, Earth Planet. Sci. Lett., 263, 104–113, https://doi.org/10.1016/j.epsl.2007.09.001, 2007a.
Tanaka, H. K. M., Nakano, T., Takahashi, S., Yoshida, J., Ohshima, H., Maekawa, T., Watanabe, H., and Niwa, K.: Imaging the conduit size of the dome with cosmicray muons: The structure beneath ShowaShinzan Lava Dome, Japan, Geophys. Res. Lett., 34, L22311, https://doi.org/10.1029/2007GL031389, 2007b.
Tanaka, H. K. M., Taira, H., Uchida, T., Tanaka, M., Takeo, M., Ohiminato, T., and Tsuji, H.: Threedimensional computational axial tomography scan of a volcano with cosmic ray muon radiography, J. Geophys. Res., 115, B12332, https://doi.org/10.1029/2010JB007677, 2010.
Tanaka, H. K. M., Kusagaya, T., and Shinohara, H.: Radiographic visualization of magma dynamics in an erupting volcano, Nat. Commun., 5, 3381, https://doi.org/10.1038/ncomms4381, 2014.