Principles of muon radiography

When passing through an object, muons interact with the matter and lose energy through ionization, bremsstrahlung, pair production, and photonuclear reactions. The energy loss behavior can generally be described as(1)where a represents the ionization term, and b is the sum of the other three radiation loss terms. Both a and b depend on the muon energy and composition of the material. Their values can be obtained from the tables provided by Groom et al. [35]. is the integral of the density of the object along the muon path, which is referred to as the density length of the object. If a muon with energy E_min loses its energy after penetrating an object with density length X₁, X₁ can be expressed as(2)where Eμ is the remaining muon mass. Combining Eqs. (1) and (2), the value of X₁ can be numerically calculated. Alternatively, Monte Carlo simulations can be used to determine the maximum densities at which muons of different energies can penetrate. The simulation results are fitted to derive a formula describing the relationship between X₁ and E_min. E_min can be determined from the measured cosmic-ray muon flux I, which can be expressed as an integral of the differential energy spectrum ϕ(E):(3)Unlike traditional imaging techniques that use artificial radiation sources, transmission muography uses naturally existing radiation sources whose energy spectra ϕ(E) cannot be manually modulated. Moreover, because the targets in transmission muography experiments are usually large, using detectors to measure the cosmic-ray muon energy spectrum ϕ(E) in real-time, before the muons penetrate the target, is impractical. Nonetheless, theoretical studies and experimental measurements show that the sea-level cosmic-ray muon energy spectrum is relatively stable and can be reliably modeled. The angular distribution of sea-level cosmic-ray muons is approximately , where θ is the zenith angle of the cosmic-ray muons. Their energy spectra follow a power law with a negative exponent. Therefore, the transmitted muon flux mainly consisted of muons with energies near E_min. Classic formulas such as Gaisser/Tang and Reyna are commonly used to describe the cosmic-ray muon energy spectrum at the sea level because they provide reliable descriptions in the energy range most relevant to transmission muography [36].

Combining Eqs. (1)–(3) and the sea-level cosmic-ray muon energy spectrum model, the density length of the object in different directions can be obtained, resulting in a 2D projection map of the object’s density length X.

3D image reconstruction of muon tomography

Based on a 2D projection map of X obtained from Sect. 2.1, the 3D density distribution can be solved for. Let m be the number of directions. By uniformly dividing the 3D region of interest (ROI) into n voxels, the integral can be discretized as(4)where xi represents the density length along the ith ray, ρj is the density of the jth voxel, and is the intersection length of the ith ray within the jth voxel. The relationship between the density lengths along different directions and density distribution can be expressed as a system of linear equations:(5)where represents the density lengths of the object along different directions, represents the intersection length of rays within the voxels along different directions, and represents the densities of the voxels.

Because the density length data are obtained from cosmic-ray muon flux data, which are usually noisy and incomplete, the linear equations in Eq. (5) faces issues such as being underdetermined and inconsistent and having a large sparse system matrix. Algebraic reconstruction techniques are commonly employed to address these problems. The most basic algebraic reconstruction algorithm is the ART algorithm [37], which is based on the Kaczmarz optimization algorithm. The ART principle can be understood as the gradual projection of a trial solution onto the hyperplanes determined by each equation to obtain a new solution. If the system of equations has a unique solution, then the trial solution converges iteratively. In every iteration of the ART algorithm, each equation of Eq. (5) updates ρ once. This equation-by-equation update method is sensitive to noise and can cause divergence if the noise is severe. The SIRT algorithm was developed based on the basic ART algorithm [38]. In each iteration of the SIRT algorithm, the contribution of each equation to ρ is first calculated without updating ρ. After going through all the equations, their contributions are averaged and ρ is updated accordingly. The updated equation for SIRT is(6)where k represents the number of iterations and λ is used to control the iteration step size, which is usually set to 1.

Both ART and SIRT solve the reconstruction problem based solely on the projection data X, without introducing additional constraints. However, muon tomography often lacks sufficient views, resulting in the ill-posedness of the image reconstruction. To improve the quality of the reconstructed image further, additional constraints must be introduced during the reconstruction procedure. We observed that the density distribution of the targets in the muon tomography analysis typically exhibited piecewise constancy, indicating that its gradient-magnitude image was sparse. This feature allowed using TV minimization as a constraint.

Let the voxel density at the coordinates be . Then, the TV in the density image is.(7)Using the gradient descent algorithm to solve the TV-minimization problem, we let:(8)(9)(10)and(11)for , and perform TV gradient descent:(12)where α controls the TV gradient descent step size, and the total number of gradient descent iterations is controlled by N_TV. In Ref. [32] TV minimization is combined with ART; however, considering that data noise is more significant in muon tomography, we combined SIRT with TV minimization for image reconstruction. The SIRT-TV procedure terminates when the number of iterations reaches N_max. The overall process of the muon tomography image reconstruction is illustrated in Fig. 1. Notably, after the SIRT-step in each iteration k, a non-negative constraint is applied to . Any additional prior information related to ρ can also be incorporated in this step.

Fig. 1Full-screen View

(Color online) Flowchart of muon tomography image reconstruction using the SIRT-TV algorithm

Muography simulation

The simulation is performed using Geant4, which is a Monte Carlo toolkit developed by CERN based on C++ [39]. Geant4 can simulate the transport and interactions of particles in matter and is a commonly used and reliable simulation platform for transmission muography. The simulation model can be divided into three parts: the QinShiHuang tomb phantom, muon track detector, and cosmic-ray muon source.

The QinShiHuang tomb was studied using 20 geophysical approaches as part of the National 863 Hi-tech Project. Archaeological results confirmed the depth of the underground palace, as well as the size and location of the tomb chamber and surrounding walls [40]. Based on these archaeological data, we built a phantom model of the QinShiHuang tomb for our simulation, which included the mound, rammed earth, loam wall, stone wall, burial chamber, and surrounding land, as shown in Fig. 2.

Fig. 2Full-screen View

(Color online) Schematic of the QinShiHuang tomb phantom. (a) 3D view of the phantom. (b) Cross-sectional view of the phantom at z=0, with red dots indicating the projections of the detector positions on the cross-section. (c) Cross-sectional view of the phantom at x=0. (d) Cross-sectional view of the phantom at y=0. The rectangles in (c) and (d) represent the ROI for 3D image reconstruction

The cosmic-ray muon source is sampled using the Reyna formula [41]:(13)where θ is the zenith angle of cosmic-ray muons, p is the momentum of cosmic-ray muons in GeV/c, , c₁=0.00253, c₂=0.2455, c₃=1.288, c₄=-0.2555, and c₅=0.0209. This formula is suggested for a momentum range of 1 Gev/c<p<2000/cosθ Gev/c and zenith angle range of 0°≤θ≤90°.

The coordinate system is defined as shown in Fig. 2(b)–(d). The muon track detector is a 1 m × 1 m ideal detector placed horizontally, recording cosmic-ray muons’ direction. The detector positions are 5 m away from the underground palace boundaries. Notably, 12 detector locations are evenly spaced at 30° surrounding the palace center, numbered 1 to 12, as shown in Fig. 2(b). Because cosmic-ray muons travel downward, detectors must be positioned below the palace level to ensure that the field of view covers the ROI. The detector depth affects the angular size of the ROI, overburden thickness, and muon zenith angles. We set the detectors’ burial depth to 90 m, as shown in Fig. 2(c), where maximum zenith angle of muons crossing ROI is 75°. Following Ref. [42], we determined the minimum muon energy required to reach the detector through the phantom, sampling only muons above this threshold. The simulated cosmic-ray muon statistics per view equal a 180-d physical measurement.

Image reconstruction

Since the tomb phantom is much larger than the detectors, they are treated as points for calculating muon flux and density length maps. Simulated cosmic-ray muon fluxes after penetrating the QinShiHuang tomb phantom, detected by each detector, are divided into 75 × 180 directions by zenith and azimuth angles, as shown in Fig. 3. Each pixel’s direction in Fig. 3 is represented by its center coordinate. The minimum energy E_min of cosmic-ray muons penetrating the phantom is calculated using Eq. (3). For density length conversion, Geant4 simulates the average muon penetration range in loess, with results in Fig. 4. A second-order polynomial is fitted to the curve in Fig. 4.

Fig. 3Full-screen View

(Color online) Simulated cosmic-ray muon fluxes penetrating the QinShiHuang tomb phantom, measured by 12 detectors using Geant4. The numbers above the figures correspond to the detector IDs in Fig. 2. θ and φ represent the zenith and azimuth angles of the incident muons in the horizontal coordinate system

Fig. 4Full-screen View

(Color online) The average penetration range of muons in loess as a function of muon energy obtained from Geant4 simulations

Since we are interested in the underground palace structure, we define the ROI as shown in Fig. 2(c) and (d) and perform 3D reconstruction for density distribution within the ROI. The density length data from different directions includes contributions from inside and outside the ROI. For the tomb phantom, knowing the density distribution outside the ROI allows us to calculate its density length. Subtracting this from the obtained data gives density length within the ROI, as shown in Fig. 5. Some directions have negative density length values due to cosmic-ray muon flux fluctuations. Only non-negative values are used in reconstruction.

Fig. 5Full-screen View

(Color online) Density length maps of the tomb phantom ROI converted from the 12-view simulated muon fluxes

In 3D image reconstruction, the ROI is divided into 60 × 50 × 20 voxels, with a voxel size of 3 m × 3.2 m × 4 m. This voxelization ensures appropriate voxel numbers for solving linear equations in Eq. (5) while maintaining resolution for fine structure reconstruction. SIRT and SIRT-TV reconstruct images from the same density length data for comparison. Following Ref. [32], we set λ=1, α=0.2, and N_TV=20. Initial density values are set to ρ₀=1.6 g/cm³, matching the loess density in the phantom. Both algorithms iterate N_max=50 times. An additional constraint in each iteration restricts voxel density outside the burial mound to 0.00129 g/cm³, matching the phantom’s air density. Slices at three representative locations are selected for comparison with phantom slices, as shown in Fig. 6.

Fig. 6Full-screen View

(Color online) Comparison of the slices of 3D density distribution between the phantom and the reconstruction results of SIRT and SIRT-TV. The dashed lines indicate the locations of the slices

A comparison of the reconstructed image slices, obtained by the two algorithms, with those of the phantom reveals that the shapes and positions of the tomb chamber and walls in the SIRT-TV reconstruction results match those of the phantom. By contrast, the slices obtained using the SIRT algorithm exhibit significant salt-and-pepper artifacts, particularly at the bottom of the reconstructed image, around the tomb chamber, and at the top of the burial mound. These artifacts can potentially lead to misinterpretation of the underground palace’s structure. The reconstruction results obtained using SIRT-TV do not show such significant artifacts, preliminarily verifying that the inclusion of TV minimization can effectively suppress artifacts.

Reconstructed image quality assessment

Figure 6 shows a visual comparison of the reconstructed images obtained using SIRT and SIRT-TV. However, visual comparison lacks quantitative information and is noncomprehensive because the reconstructed images are 3D. To quantitatively compare the reconstruction results of SIRT and SIRT-TV, using appropriate image quality assessment metrics is crucial. Classic metrics include mean squared error (MSE), peak signal-to-noise ratio (PSNR), and structural similarity (SSIM) [43]. These metrics can evaluate the overall similarity between the reconstructed image and phantom. However, this study is primarily focused on the feasibility of accurately reconstructing the structures of the underground palace, specifically the tomb chamber and walls. Referring to the image quality assessment metrics used in Refs. [28], we use the Jaccard similarity index J to assess the geometric similarity of the reconstructed images to the phantom for the tomb chamber and walls. J is defined as follows:(14)where P represents the set of voxels belonging to the tomb chamber or walls in the phantom and R represents the set of voxels belonging to the same structure in the reconstructed image. The term denotes the number of elements at the intersection of P and R, and denotes the number of elements in the union of P and R. A larger indicates a higher geometric similarity, and when , the geometric structures of the phantom and reconstructed image are identical. When calculating for the images reconstructed from Sect. 3, we encounter the issue of inconsistent voxel segmentation between the phantom and reconstructed images because the phantom established in Geant4 is not voxelized. To address this problem, we modified Eq. (14). We define a vector p, which represents the proportion of overlap between the voxels defined in the reconstructed images and the corresponding structure in the phantom. Consequently, element pi corresponding to the ith voxel can take three types of values:(15)We define a vector r to determine whether each voxel belongs to a structure in a reconstructed image. Element ri corresponding to the ith voxel can take two values:(16)Therefore, Eq. (14) can be redefined as:(17)In this phantom study, partial overlaps exist between voxels and structures; therefore, J(p, r) cannot reach 1. When calculating r, image segmentation extracts tomb chamber and wall structures from reconstructed images. Since tomb chamber and wall densities correspond to minimum and maximum densities respectively, segmentation uses density thresholds to extract these structures based on voxel density values. For the tomb chamber:(18)And for the walls:(19)Ideally, thresholds should be set according to structure densities in the phantom. Due to differences between reconstructed image densities and the phantom—especially higher tomb chamber densities compared to that of air in the phantom—an optimization approach finds suitable thresholds for segmentation. Because chamber density is lower than loess, while wall density is higher, the chamber threshold is searched within (0, 1.6 g/cm³] and wall threshold within (1.6 g/cm³, 2.7 g/cm³], using 0.1 g/cm³ steps. For each threshold, J is calculated, with highest J taken as geometric similarity, denoted as J_chamber for the chamber and J_wall for walls. Thresholds are searched for reconstructed images from Sect. 3. Voxel sets for chamber and walls are extracted using thresholds with highest J, and positions from segmentation are shown in Fig. 7.

Fig. 7Full-screen View

(Color online) Comparison between the structures extracted from reconstructed images obtained using two algorithms and the structures in the phantom, based on the best density thresholds. The blue crosses and red dots represent the voxel positions of the chamber and the walls in the reconstructed images, respectively. The blue box outlines the boundaries of the tomb chamber in the phantom, and the red box outlines those of the walls

Both the wall and tomb chamber geometries in the SIRT reconstruction result show noticeable overflow compared to the phantom. In contrast, the counterparts extracted from the SIRT-TV reconstructed image align better with the phantom. Table 1 lists the J-calculation results for the structures in the reconstructed images obtained using the two algorithms. Both J_chamber and J_wall in the SIRT-TV reconstruction results are superior to those obtained using SIRT. In particular, J_wall is 2.8 times that of SIRT.

Iterative convergence

Convergence speed is an important metric for evaluating iterative image reconstruction algorithms. In Sect. 4.1, J_chamber and J_wall of the reconstruction results are compared after the same number of iterations for both algorithms. However, the number of iterations required for each algorithm to converge may differ. To compare the convergence speeds of the two algorithms, they were iterated 150 times, and the curves showing the changes in J_chamber and J_wall with the number of iterations are plotted in Fig. 8. J_chamber and J_wall of SIRT-TV generally increased with the number of iterations until convergence, whereas those for SIRT initially increased but then decreased with additional iterations before convergence.

Fig. 8Full-screen View

(Color online) J and ΔJ/J of the reconstructed images obtained by SIRT and SIRT-TV as functions of iterations

To quantitively evaluate the convergence speed, we defined a convergence metric as follows:(20)where k represents the number of iterations and ε is the threshold for convergence. The ΔJ/J values for both algorithms are plotted in Fig. 8 as functions of the number of iterations. ΔJ/J for both algorithms fluctuate while converging due to structure thresholds and voxelization discreteness. With ε set to 0.002, J_chamber and J_wall for SIRT-TV converge after 36 and 24 iterations, while SIRT has not yet converged. SIRT-TV reaches convergence faster due to TV-minimization constraint. In SIRT, J_chamber and J_wall converge at different speeds because the tomb chamber is a large cubic volume, while the walls are complex and thinner.

Impact of TV parameters on reconstruction

The parameter values in the TV-minimization algorithm may affect the quality of the reconstructed images. To assess this impact, we varied the maximum iteration number N_TV and the TV gradient descent step size α. J_chamber and J_wall are shown as functions of the number of iterations for different parameter values in Fig. 9.

Fig. 9Full-screen View

(Color online) J_chamber and J_wall of reconstructed images obtained under different TV parameters, as functions of iterations

For most combinations, convergence occurs within 50 iterations, indicating that choosing a maximum iteration number N_max=50 as the termination criterion for SIRT-TV is reasonable. For J_wall, when N_TV=10 and α=0.1, the trend resembles that of SIRT, and a continuous downward trend is visible after 150 iterations, showing the ineffectiveness of TV-minimization. With different α values, J_wall is always the lowest when N_TV=10. Unlike J_wall, for different α values (except α=0.1), J_chamber is always the highest when N_TV=10, and remains above 0.5 irrespective of the parameter values. This result suggests that compared to the walls, the reconstruction of the tomb chamber does not require a stringent TV constraint.

Overall, the different parameter values in Fig. 9 (except N_TV=10 or α=0.1) do not significantly affect the reconstruction results, indicating that the SIRT-TV algorithm can reliably produce high-quality images.

Impact of measurement configuration and duration

As mentioned in Sect. 1, the main challenges in 3D reconstruction of muon tomography are the presence of noise and incomplete data. In Sect. 3, we simulated relatively ideal conditions to test the feasibility of the algorithm. This section discusses the impact of less-ideal measurement configurations and durations (i.e., fewer views and lower muon statistics) on reconstruction with SIRT-TV.

4.4.1

Measurement configuration

Using the same simulation data from Sect. 3, the number of views is reduced to four, and three different detector configurations are selected for reconstruction. The parameters of the image reconstruction algorithm were the same as those used in Sect. 3. Figure 10 shows slices of the phantom and the reconstructed images for the 12-view and three 4-view configurations. Reducing the number of views leads to more severe artifacts. Among these 4-view configurations, the symmetrical 4-view configuration shows better geometric symmetry in its reconstruction results, whereas the reconstruction results for the other two 4-view configurations exhibit distortions in the shapes of the tomb chamber and walls. The calculation results of J_chamber and J_wall in Table 2 also indicate that, among the three different 4-view configurations, the reconstructed image quality is the highest for the symmetrical 4-view configuration. These artifacts can be explained by the ray densities in the voxels. Given that the voxelization for each reconstruction is the same, more views lead to more rays, higher ray densities in the voxels, and more information regarding voxel density. When voxelization and the number of views were the same for each reconstruction, asymmetrical measurement configurations caused greater variation in ray densities across voxels. This implies that some voxels are crossed by only a few rays, making it difficult to reconstruct their densities. A symmetrical measurement configuration results in better reconstruction quality.

Fig. 10Full-screen View

(Color online) Slices of the phantom, and reconstructed images for the 12-view configuration and three 4-view configurations

Table 2Full-screen View

J_chamber and J_wall of the images obtained from the 12-view reconstruction and three 4-view reconstructions

Measurement configuration	Jchamber	Jwall
12-view	0.5929	0.5498
4-view configuration 1	0.5161	0.4552
4-view configuration 2	0.4818	0.4452
4-view configuration 3	0.4576	0.3777

Show more

4.4.2

Measurement duration

The number of cosmic-ray muons detected by the detectors is proportional to the measurement duration. This affects the statistical fluctuations, which, in turn, affects the quality of the reconstructed image. To investigate the reconstructed image quality for different measurement durations, we select datasets from the beginning of the simulation data in Sect. 3, corresponding to 7, 14, 30, 90, and 180 d of equivalent measurement time. The algorithm parameters are the same as those used in Sect. 3. Slices of the reconstructed images are shown in Fig. 11.

Fig. 11Full-screen View

(Color online) Slices of the phantom and SIRT-TV reconstructed images for different equivalent measurement durations

In the reconstructed image slices for 7-d and 14-d measurement durations, severe artifacts appear between the chamber and walls, especially at the bottom. For measurement durations beyond 30 d, artifacts are suppressed, and the chamber-earth boundary becomes clearer. The 90-d and 180-d reconstructed image slices show similar results, with J_chamber and J_wall both , as shown in Table 3.

Table 3Full-screen View

J_chamber and J_wall of reconstructed images obtained using SIRT-TV under different equivalent measurement durations

Time (day)	Jchamber	Jwall
7	0.2615	0.2050
14	0.4514	0.2465
30	0.5123	0.3383
90	0.5525	0.4994
180	0.5929	0.5498

Show more