Target scanning and 3-d forward and inversion method based on cosmic ray muons

Abstract

The application discloses a target scanning and three-dimensional forward and inverse method based on cosmic muons, which comprises the following steps: (1) fitting cosmic muon flux and energy spectrum outside the target according to experimental data; (2) determining the interaction model of muons and matter; (3) collecting the shape structure outside the target; (4) assuming that the target has a certain internal structure; (5) preparing the detector measurement response curve data; (6) constructing a three-dimensional muon interaction forward model; (7) generating the muon counting rate or survival rate in each direction received by the detector at the preset position through calculation; (8) converting the counting rate or survival rate into the equivalent medium length based on experimental measurement; (9) performing three-dimensional inversion calculation to obtain the real internal structure of the target; (10) then performing post-processing on the inversion result based on zero hypothesis, alternative hypothesis and other statistics.

Description

Technical Field

[0001] This invention relates to the field of cosmic ray muon imaging, and more specifically to a target object scanning and three-dimensional forward and inverse modeling method based on cosmic ray muons. Background Technology

[0002] Cosmic ray muons are secondary ray particles generated by the interaction of high-energy particles from interstellar space with atomic nuclei in Earth's atmosphere. They possess a wide energy range and high energy, allowing them to penetrate large targets and offering unparalleled advantages over other physical detection methods, thus promising broad application prospects. Cosmic ray muon imaging technology is a crucial area for the development of the green nuclear technology industry in the 21st century, often referred to as a "geo-earth fault scanning imaging CT." Research on intensity-attenuated cosmic ray muon imaging technology has been applied to mineral exploration. Developed three-dimensional inversion imaging techniques have been used to measure and reconstruct the spatial distribution of underground igneous copper, lead, zinc, and other non-ferrous metal ore bodies, achieving excellent mineral exploration results. In imaging targets of a certain scale, it offers unparalleled advantages over artificial ray devices. Currently, conventional mineral exploration techniques, such as magnetics, electromagnetism, gravity, and seismic methods, are selective in identifying the lithology of geological bodies and offer multiple interpretations and judgments of anomalies, resulting in various application limitations and drawbacks in the field of mineral exploration. Mun imaging technology can perform two-dimensional and three-dimensional imaging of geological bodies with different densities and lithologies underground, thereby accurately delineating different geological bodies and overcoming the shortcomings of other geophysical exploration methods. This represents a significant innovation in mineral exploration technology. However, current mun imaging algorithms cannot effectively identify non-target information and cannot accurately obtain the uncertainty of the anomaly distribution of the inversion results' cell density. To address these issues, this invention proposes a method for detection, forward and inverse modeling, and statistical data processing using a track mun detection system. Summary of the Invention

[0003] The purpose of this invention is to propose a method for detection, forward and inverse modeling, and statistical data processing using a track muon detection system. This detection method can be used not only on the ground or in existing tunnels, but also in small-diameter boreholes for cosmic ray muon scanning to identify geological bodies below the surface.

[0004] The present invention specifically adopts the following technical solution:

[0005] The target object scanning and 3D forward and inverse modeling method based on cosmic ray muons includes the following steps:

[0006] (1) The cosmic ray muon flux and energy spectrum outside the target object were obtained by fitting the experimental data;

[0007] (2) Determine the interaction model between muons and matter;

[0008] (3) Select a target object shape or use a target object shape from a real application scenario, and obtain the outline of the target cultural relic;

[0009] (4) Assume the internal properties of a target object, including its shape, density distribution, and location;

[0010] (5) Place the detector at certain specific locations around or inside the target object; (so that the direction of the muon rays covers the area of ​​the target object with the unknown density distribution);

[0011] (5) Measure the detector response curve for later use;

[0012] (6) Using (1)-(4) as known input conditions, construct the Muon forward model;

[0013] (7) Calculate the muon count rate or survival rate of the detector at the preset position in each direction by using the detector response curve (5) and the forward model (6), which is the calculation data based on the forward model.

[0014] (8) Convert the measured count rate or survival rate into the experimentally measured equivalent medium length;

[0015] (9) Compare the equivalent medium length obtained by forward modeling and actual measurement, and combine the geological preset conditions to continuously iterate the forward geological model and adjust the inversion parameters until a satisfactory geological model iteration result is obtained, that is, the target density distribution imaging result.

[0016] (10) Based on the null hypothesis and alternative hypothesis, determine the non-target factors and the inversion uncertainty.

[0017] Preferably, muons lose energy during their interaction with matter, and after losing all their kinetic energy, they decay into electrons and neutrinos.

[0018] The energy loss of muons mainly occurs through four mechanisms: ionization excitation, bremsstrahlung, electron-positron pair production, and photosynthetic inelastic scattering. For muons with energies below 500 GeV, ionization excitation is the primary mode of energy loss. The rate satisfies equation (1).

[0019]

[0020] Where Z and A are the atomic number and atomic mass of the substance being acted upon, respectively, ρ is the density of the substance, the relativistic coefficient β = v / c, v is the muon velocity, c is the speed of light, and K = 2m. e / I 2 I is the average excitation potential (on the order of eV), m e For electron mass (m) e=0.511MeV), since Z / A = 1 / 2 for most elements, K and Z are roughly linearly related, and K contributes less to the energy loss rate in the logarithmic term than in the proportional term. Therefore, according to the above analysis, the energy loss rate of muons is almost only related to the matter density ρ.

[0021] Preferably, the cosmic ray muon flux and energy spectrum at the Earth's surface can be obtained either through quasi-empirical formulas fitted from experimental data or through Monte Carlo simulation. Preferably, the region to be determined beneath the surface is divided into several discrete small grids, and the inversion process involves solving for the density value in each grid. Assuming the density value at a certain grid point (x,y) is ρ(x,y), the equivalent stratum length L traversed by the muon ray in the ray direction j at detector position i is calculated. ij Represented as:

[0022] L ij =∑ k ρl ijk (2)

[0023] Among them, l ijk To determine the track length of a muon within a given grid cell along its path, and taking into account all ray directions and detector positions, we can formulate a matrix equation:

[0024]

[0025] Where (L1,L2,…,L) N ) represents the equivalent length corresponding to the muon decay information, denoted as L; (ρ1, ρ2, ..., ρ M ) represents the density value in each cell after dividing the region, i.e., the density distribution, and ρ is the variable to be solved. Let l be the length of the track of the muon passing through each cell for each data point. Since the above equation is often an underdetermined equation, solving the equation generally requires minimizing the objective function in an iterative manner. The objective function is a linear superposition of "Misfit" and "Norm", that is, Objective function = Misfit + β × Norm. "Misfit" represents the difference between the measured data and the predicted model calculation results, and "Norm" is the scale of whether the model conforms to the prior information. The inversion achieves the minimization of the objective function through iteration, adjusts the parameter values, and obtains the optimal density distribution model.

[0026] Preferably, the inversion methods involved include two types: Method 1 is the null hypothesis equivalent length, that is, for directions where the difference between experimental and predicted values ​​is large, the predicted value is used to replace the experimental value; Method 2 is the equivalent length of experimental measurement. Based on these two lengths, the inversion methods include two types. For Method 1, the difference between the two types of equivalent lengths is used as input data, i.e., (L1, L2, ..., L...) in the inversion formula. N Using the topological data of the imaging area, the density of the surrounding rock, and other constraints as known conditions, inversion is performed to obtain the density distribution based on the difference in equivalent lengths. For method 2, three-dimensional inversion is performed on the experimental data and the null hypothesis data respectively to obtain two density distributions based on the experimental data and the null hypothesis. The difference between their density distributions is the result of method 2.

[0027] Preferably, the determination of non-target factors and inversion uncertainty is based on the null hypothesis. The null hypothesis is that the geological model assumes uniform density and that there are no density anomalies in the strata. The inversion result based on the null hypothesis means that any density non-uniformity in the inversion result is not caused by geological factors, but by other background factors such as the algorithm. Therefore, the result of the null hypothesis is considered as the background result.

[0028] Preferably, the inversion process based on alternative hypotheses and experimental data involves a sampling process and multiple inversions to obtain an inversion ensemble based on the alternative hypotheses. This inversion ensemble approximates multiple measurements of the target object and is a collection of measurement results. Based on this ensemble, we can obtain the probability of a similar density distribution pattern appearing after multiple measurements, and obtain the density uncertainty by calculating the variance of the density distribution in each grid cell.

[0029] The present invention has the following beneficial effects:

[0030] The target object scanning and three-dimensional forward and inverse modeling method based on cosmic ray muons described in this application can be used not only on the ground or in existing tunnels, but also in small-diameter boreholes for cosmic ray muon scanning to identify geological bodies (including ore bodies) below the surface, effectively solving the problems of deep geological scientific research and blind ore body delineation. Attached Figure Description

[0031] Figure 1 Framework for cosmic ray muon target scanning and 3D forward and inverse modeling algorithms;

[0032] Figure 2 A schematic diagram of the measurement of geological anomaly targets based on cosmic ray muons;

[0033] Figure 3 A schematic diagram showing the differences in counting at the locations of three detectors placed underground and at location D0 due to the presence or absence of density anomalies.

[0034] Figure 4 This is the typical forward and inverse modeling process for muon imaging;

[0035] Figure 5 A schematic diagram of key elements in the inversion of the Earth's surface, such as grids, detector locations, and rays;

[0036] Figure 6 Schematic diagram of two inversion algorithms for data processing of muon 3D imaging;

[0037] Figure 7 This is an inversion process based on experimental data;

[0038] Figure 8 This is a single inversion process based on the null hypothesis;

[0039] Figure 9 Generates a sampling-based null hypothesis ensemble (multiple inversions);

[0040] Figure 10 This is a sampling-based method for generating alternative hypothesis ensembles and calculating density uncertainty.

[0041] Figure 11 The results of Method 1 (solid blocks) and Method 2 (blocks outlined by lines) based on experimental data of a city wall rampart are shown.

[0042] Figure 12 The correlation between the density distributions of Method 1 and Method 2 based on experimental data from a city wall rampart. Figure 6 (The two methods shown);

[0043] Figure 13 Based on experimental data from a city wall embankment, denoted as , where 'a' is the average density distribution obtained from the ensemble of the alternative hypothesis, 'b' is the density uncertainty distribution obtained from the ensemble of the alternative hypothesis, 'c' is the average density distribution obtained from the ensemble of the null hypothesis, and 'd' is the correlation between the average density of the ensemble of the alternative hypothesis and the average density of the ensemble of the null hypothesis. Detailed Implementation

[0044] The specific embodiments of the present invention will be further described below with reference to the accompanying drawings and specific examples:

[0045] The target object scanning and 3D forward and inverse modeling method based on cosmic ray muons includes the following steps:

[0046] (1) Calculate the cosmic ray muon flux and energy spectrum at the Earth's surface;

[0047] (2) Determine the interaction model between muons and matter; calculate the cosmic ray muon flux and energy spectrum at the Earth's surface. This can be obtained either by a quasi-empirical formula fitted to experimental data or by Monte Carlo simulation.

[0048] (3) Assume a surface shape or use a surface shape from a real-world application scenario;

[0049] (4) Assume a target object, including its shape, density distribution and location;

[0050] (5) Collect detector response curves for later use;

[0051] (6) Using (1)-(4) as known input conditions, construct a forward model;

[0052] (7) Calculate the muon count rate of the detector at the preset position in each direction through the detector response curve (5) and the forward model (6), which is the “experimental” data obtained through simulation.

[0053] (8) Convert the count rate of (7) into the equivalent medium length of the “measurement”;

[0054] (9) Perform inversion calculation by combining the equivalent medium length calculated using surface topography data and the "measured" equivalent medium length;

[0055] (10) Compare the results of the inversion calculation with the assumed known target. If the degree of agreement does not meet the preset conditions, continue to adjust the inversion parameters and model equations until a satisfactory imaging result is obtained.

[0056] (11) Based on the null hypothesis and alternative hypothesis, determine the non-target factors and the inversion uncertainty.

[0057] In geophysics, the "forward problem" is a method of predicting measurement results (predicted data) based on geological models (geological structures or properties), principles (such as interactions), and known conditions. The "inversion problem," on the other hand, is the reverse process of determining the parameters of a geological model based on relevant principles or physical models, given measurement data. Generally, because the forward model is deterministic and its parameters are complete and noise-free, the solution to the forward problem is unique. In the inversion problem, however, the actual data always contains statistical and systematic errors, and the information to be solved is far more complex than the known data, making it an underdetermined problem with multiple solutions. Typically, during cosmic ray muon geological scanning, a muon track detector system is placed below the Earth's surface to receive muon rays from different directions above.

[0058] like Figure 2 As shown, muons lose energy during their interaction with matter, and after losing all their kinetic energy, they decay into electrons and neutrinos.

[0059] Muons lose energy primarily through four mechanisms: ionization excitation, bremsstrahlung, electron-positron pair production, and photosynthetic inelastic scattering. For muons with energies below 500 GeV, ionization excitation is the main mode of energy loss during their interaction with matter. The rate satisfies equation (1).

[0060]

[0061] Where Z and A are the atomic number and atomic mass of the substance being acted upon, respectively, ρ is the density of the substance, the relativistic coefficient β = v / c, v is the muon velocity, c is the speed of light, and K = 2m. e / I2, where I is the average excitation potential (on the order of eV), m e For electron mass (m) e =0.511MeV). Since most elements Z / A≈1 / 2, K and Z are roughly linearly related. Furthermore, K contributes less to the energy loss rate in the logarithmic term than in the proportional term. Therefore, according to the above analysis, the energy loss rate of muons is almost only related to the matter density ρ.

[0062] For example, the minimum energy required for a muon to penetrate 8m of standard rock or 20m of water is approximately 4 GeV. When a muon passes through an imaging region, energy is lost, and flux attenuates. The amount of flux attenuation is related to the density distribution of the imaging region. Therefore, if the density distribution in the imaging region is known, the attenuation information of the muon flux (count) in each direction can be calculated based on the muon energy loss law and the energy spectrum of the incident muons. This process maps the density distribution model to the probe data and is called forward modeling. Conversely, if the attenuation information of muons in each direction is known, such as the muon azimuth count or survival rate, the density distribution in the imaging region can be inferred from the muon energy loss law and the energy spectrum of the incident muons. This process maps the probe data to the density distribution model and is called inversion. Azimuth counts are important in scenarios with or without density anomalies, and the forward and inverse modeling methods are similar. Figure 3 , Figure 4 As shown.

[0063] The area to be explored beneath the Earth's surface is divided into several discrete small grids, such as... Figure 5 As shown, the inversion process involves solving for the density value in each grid cell. Assuming the density value at a certain (x,y) grid point is ρ(x,y), the equivalent formation length L through which the muon ray passes in the ray direction j at detector position i is... ij It can be represented as:

[0064] L ij =∑ k ρl ijk (2)

[0065] Among them, l ijkTo determine the geometric length within a grid cell along the ray path, taking into account all ray directions and detector positions, we can formulate a matrix equation:

[0066]

[0067] Where (L1,L2,…,L) N ) represents the equivalent length corresponding to the muon decay information, denoted as L; (ρ1, ρ2, ..., ρ M ) represents the density value in each cell after dividing the region, i.e., the density distribution, and ρ is the variable to be solved. Let l be the length of the track of the muon passing through each cell for each data point. Since the above equation is often an underdetermined equation, solving the equation generally requires minimizing the objective function in an iterative manner. The objective function is a linear superposition of "Misfit" and "Norm", that is, Objective function = Misfit + β × Norm. "Misfit" represents the difference between the measured data and the predicted model calculation results, and "Norm" is the scale of whether the model conforms to the prior information. The inversion achieves the minimization of the objective function through iteration, adjusts the parameter values, and obtains the optimal density distribution model.

[0068] Due to the ambiguity of the solution, the density distribution obtained through inversion needs to address the following issues:

[0069] 1. Does the measured density distribution or density anomaly match the actual geological conditions? What is the degree of positive correlation? If similar patterns still appear after multiple measurements, it indicates that it is not caused by statistical noise. Using the alternative hypothesis ensemble method, calculate the central value and standard deviation of the density distribution at each grid point to evaluate the experimental inversion results.

[0070] 2. What morphologies in the solved density distribution are caused by non-geological factors? Could systematic errors caused by inversion algorithms, measurements, etc., lead to the appearance of non-geological morphologies in the density distribution? How can non-geological morphologies be subtracted from the density distribution?

[0071] 3. Assuming multiple measurements are performed, will a similar density distribution pattern still appear?

[0072] 4. How to assess the uncertainty of density in each grid cell?

[0073] 5. How to determine the threshold of background morphological density? That is, at what threshold is the density considered background (or uniform density), and above the threshold is the density anomaly?

[0074] By proposing, for example Figure 6The data post-processing architecture shown addresses the above problem. Methods 1 and 2 in the diagram primarily solve problem 2. For Method 1, it starts with two types of equivalent lengths: one is the null hypothesis equivalent length, where the predicted value is replaced by a sampled value for directions where the difference between the experimental and predicted values ​​is large; the other is the equivalent length of the experimental measurement. The difference between these two types of equivalent lengths is used as input data, i.e., (L1, L2, ..., L...) in the inversion formula. N Using the topological data, surface data, surrounding rock density, and other constraints of the imaging area as known conditions, an inversion is performed (see...). Figure 7 This yields a density distribution based on the difference in equivalent lengths. For Method 2, three-dimensional inversion is performed on the experimental data and the null hypothesis data respectively (see...). Figure 8 Method 2 obtains two density distributions: one based on experimental data and the other based on the null hypothesis. The difference between these density distributions is the result of Method 2. The null hypothesis assumes a uniform density geological model and that there are no density anomalies in the strata. The inversion result based on the null hypothesis means that any density non-uniformity in the inversion result is not caused by geological factors, but by other background factors such as the algorithm. Therefore, the result based on the null hypothesis can be considered the background result. We therefore also introduce an ensemble of the null hypothesis, which is a collection of multiple inversions of the background density distribution. This ensemble can solve problem 5; the specific process is detailed in [link to documentation]. Figure 9 As shown. This process and the single null hypothesis inversion process (see...) Figure 7 The difference lies in the addition of a full sampling step in each process. Similarly, based on the inversion process using experimental data, we introduce a sampling process and perform multiple inversions to obtain an inversion ensemble based on alternative hypotheses. This ensemble can well answer questions 1, 3, and 4. The inversion ensemble based on alternative hypotheses can be considered as a collection of multiple measurements of the target object. Based on this ensemble, we can obtain the probability of a similar density distribution pattern appearing after multiple measurements. By analyzing the variance of the density distribution in each grid, we can obtain the density uncertainty, which is very helpful in understanding the inversion results based on experimental data.

[0075] Based on the method mentioned in this invention, the experimental data of the Muzi (a type of ancient Chinese insect) collected on-site from the city wall rampart of a certain archaeological target were analyzed and processed. Figure 11 , Figure 12 This demonstrates the results of Method 1 and Method 2 and their correlation. It can be seen that the two methods have a good correlation, and the results are very similar, both effectively suppressing morphological problems caused by background non-geological factors. Figure 13 (a) Mean density distribution of the alternative hypothesis ensemble. This distribution is very similar to the inversion results based on experimental data, meaning that similar results can be obtained even with multiple measurements; Figure 13(b) The density uncertainty distribution obtained from the alternative hypothesis ensemble, which can be used to analyze the confidence level of each grid density; Figure 13 (c) The mean density distribution of the null hypothesis ensemble, which can be used to extract the background threshold of the inversion density; Figure 13 (d) The correlation between the average density of the alternative hypothesis ensemble and the average density of the null hypothesis ensemble shows a certain correlation and difference. Areas with good correlation indicate the uniformity of density, while areas with poor correlation correspond to density anomalous areas in the internal structure of the horse face.

[0076] Of course, the above description is not intended to limit the present invention, and the present invention is not limited to the examples given above. Any changes, modifications, additions or substitutions made by those skilled in the art within the scope of the present invention should also fall within the protection scope of the present invention.

Claims

1. A method for object scanning and 3D forward and inverse problem based on cosmic muons, characterized in that, The method comprises the following steps: (1) fitting the cosmic muon flux and energy spectrum outside the target object according to experimental data; (2) determining the interaction model of muons and matter; (3) selecting a target object shape or using a target object shape in a real application scenario, and obtaining the contour of the target object; (4) assuming the internal properties of the target object, including its shape, density distribution and position; (5) placing the detector at certain positions around or inside the target object, so that the direction of the muon rays covers the area of the target object with the to-be-determined density distribution; (6) measuring the response curve of the detector; (7) taking (1)-(4) as input known conditions to construct a muon forward model; (8) calculating the muon count rate or survival rate received by the detector at the preset position in each direction based on the forward model, that is, the calculation data based on the forward model, through the detector response curve (6) and the forward model (7); (9) converting the measured count rate or survival rate into an experimental measurement equivalent medium length; (10) comparing the equivalent medium lengths obtained by forward calculation and measurement, combining the preset geological conditions, and iteratively adjusting the inversion parameters until a satisfactory geological model iteration result, i.e., the target object density distribution imaging result, is obtained; (11) determining the non-target object factors and inversion uncertainty based on the null hypothesis and alternative hypothesis. The inversion methods include two methods. Method 1 is zero hypothesis equivalent length, that is, for the direction with large difference between experimental value and predicted value, the experimental value is replaced by the predicted value. Method 2 is the equivalent length of experimental measurement. The inversion methods based on the above two lengths include two methods. For method 1, the difference between the above two equivalent lengths is taken as input data, that is, the difference between the two equivalent lengths is taken as input data in the inversion formula Based on the imaging area topological data, the surrounding rock density and other constraint conditions, the inversion is carried out to obtain the density distribution based on the difference between the equivalent lengths. For method 2, three-dimensional inversion is carried out on the experimental data and zero hypothesis data respectively to obtain two density distributions based on the experimental data and the zero hypothesis data. Then the difference between the density distributions is the result of method 2.

2. The cosmic-ray muon-based object scanning and 3D forward-inversion method of claim 1, wherein, Muons will lose energy in the process of interaction with matter, and after all their kinetic energy is lost, they will decay into electrons and neutrinos. The energy loss modes of muons are mainly divided into four kinds: ionization excitation, bremsstrahlung, pair production, and photo-nuclear inelastic scattering. For muons with energy less than 500 GeV, ionization excitation is the main mode of energy loss, and the energy loss rate of muons satisfies formula (1) （1） where Z and A are atomic number and atomic mass of the substance being acted upon, is the density of the substance, the relativistic coefficient , v is the muon velocity, c is the speed of light, , I is the average excitation potential, eV order, m e is the electron mass m e = 0.511 MeV, since most elements , K and Z are roughly linear, and K in the logarithmic term, the contribution to the energy loss rate is less than the proportional term, so according to the analysis of the above formula, the energy loss rate of muon is almost only related to the density of the substance .

3. The cosmic-ray muon based object scanning and 3D forward and inverse method according to claim 1, wherein, The cosmic muon flux and energy spectrum at the ground surface are obtained by a quasi-empirical formula fitted by experimental data or by a Monte Carlo program simulation.

4. The cosmic-ray muon based object scanning and 3D forward and inverse method according to claim 1, wherein, The subsurface region to be solved is divided into several discrete small grids, and the inversion process is to solve the density value in each grid. It is assumed that the density value in a certain place is represented as: The position of the detector The direction of the ray The equivalent formation length through which the muon ray passes is represented as: （2） where, is the path length of the muon in the i-th bin of its path, Considering all the ray directions and detector positions, a matrix equation can be written: （3） wherein is the equivalent length corresponding to the muon attenuation information, denoted as ; is the density value in each unit after the division of the region, i.e., the density distribution, which is a quantity to be solved, denoted as , the matrix is the track length of the muon passing through each unit corresponding to each data, denoted as ; since the above equation is often an underdetermined equation, the solution of the equation generally needs to minimize the objective function in an iterative manner, and the objective function Objective function is a linear superposition of "Misfit" and "Norm", i.e. "Misfit" represents the difference between the measured data and the estimated model calculation result, and "Norm" is a scale of whether the model conforms to the prior information, and the inversion realizes the minimization of the objective function through iteration, adjusts the parameter value, and obtains the best density distribution model.

Images