Imaging method and imaging system of muon scattering

Abstract

The application provides a muon scattering imaging method, comprising the following steps: dividing a three-dimensional space of an imaging area passed by muons into voxels; in each scattering event, (1) obtaining a position of a nearest point passed by muons in the imaging area, (2) calculating a scattering angle initial value according to an incident track and an exit track of the muons passing through the imaging area and comparing the scattering angle initial value with a scattering angle threshold value, (3) calculating and counting each path length of each voxel passed by the muons from the incident point to the nearest point and each path length of each voxel passed by the muons from the nearest point to the exit point; performing voxelization and merging on the path length and the scattering angle value; calculating imaging density and corresponding imaging values, and performing image reconstruction. The present imaging method improves the clarity of muon imaging by performing voxelization and merging on the path length and the scattering angle and combining the angle cap method to reconstruct the imaging. The application also provides an imaging system using the above imaging method.

Description

Technical Field

[0001] This invention relates to the technical field of cosmic ray imaging, and in particular to an imaging method and system for muon scattering. Background Technology

[0002] Primary cosmic rays entering the atmosphere interact with atmospheric atomic nuclei to produce a large number of secondary particles, most of which are π mesons of similar numbers but different charges. Charged π mesons then decay into muons.

[0003] As a highly penetrating high-energy charged particle, the muon can easily penetrate thick shielding layers. Its rest mass is 105.7 MeV, approximately 207 times the mass of an electron, and it can carry either a positive or negative charge. As a charged particle with a mass between that of an electron and a proton, when a muon interacts with matter, μ... + More like a light proton or a heavy positron, μ - It is more like a heavy electron. Muons interact with matter primarily through three mechanisms: energy loss, complete absorption, and multiple Coulomb scattering.

[0004] Since Anderson discovered muons in 1936 by measuring the energy loss of cosmic ray particles in a cloud chamber, muon imaging can be applied to some unique fields: (1) it can image ultra-large structures such as volcanoes and pyramids; (2) it can detect whether containers contain high-Z substances; (3) it can perform non-destructive imaging flaw detection on hard-to-penetrate facilities such as nuclear reactors.

[0005] Inversion methods for muon imaging include transmission methods, the Point of Closest Approach (PoCA) method and its simplified versions, statistical methods, and rapid detection methods. Currently, when using muons for imaging, muon scattering events are related to factors such as the size and area of ​​the detector panel, the zenith angle, and the total detection time, resulting in often limited imaging data. Under these conditions, these imaging methods often face problems such as low image quality and unstable imaging results. Summary of the Invention

[0006] To address the aforementioned problems in existing technologies, this invention provides an imaging method for muon scattering. Based on voxel partitioning and image value reconstruction, this method performs voxel partitioning followed by large-voxel merging of parameters such as path length and scattering angle, and combines this with the angle capping method to reconstruct the image values ​​within voxels, thereby improving the clarity of muon imaging. This invention also provides an imaging system employing the above imaging method.

[0007] To achieve the above objectives, the present invention provides the following technical solution:

[0008] Imaging methods for muon scattering include the following steps:

[0009] (a) Divide the three-dimensional space of the imaging region traversed by the muon into voxels;

[0010] (b) In each scattering event,

[0011] (1) Obtain the position of the nearest point in the imaging region through which the muon passes.

[0012] (2) Calculate the initial scattering angle based on the incident and exit tracks of the muon through the imaging region; compare the initial scattering angle with the scattering angle threshold:

[0013] When the initial value of the scattering angle is greater than the scattering angle threshold, the scattering angle value is statistically calculated as the square of the scattering angle threshold.

[0014] When the initial scattering angle is not greater than the scattering angle threshold, the scattering angle value is statistically calculated as the square of the initial scattering angle itself.

[0015] (3) Calculate and count the path length of each voxel from the incident point to the nearest point, and calculate and count the path length of each voxel from the nearest point to the exit point.

[0016] (c) Perform large voxelization merging based on the statistical path length and scattering angle values. The large voxelization merging is to merge the voxel data in the grid adjacent to the central small voxel.

[0017] (d) Based on the scattering angle and path length after bulk merging, the imaging density and corresponding imaging values ​​are calculated, and image reconstruction is performed.

[0018] In one embodiment of the present invention, after obtaining the position of the nearest point, the voxel coordinates of the nearest point are calculated.

[0019] In one embodiment of the present invention, when dividing the three-dimensional space of the imaging region into voxels, the region is uniformly divided in the x, y, and z directions, thereby dividing it into voxels in the form of small cuboids.

[0020] In one embodiment of the present invention, the typical value of the scattering angle threshold is the statistical standard deviation of the initial scattering angle values ​​collected and calculated.

[0021] In one embodiment of the present invention, the calculation of the initial value of the scattering angle and the comparison with the scattering angle threshold includes the following steps:

[0022] (b1) In each scattering event, detect the incident line corresponding to the incident trajectory when the muon passes through the voxel, and the exit line corresponding to the exit trajectory.

[0023] (b2) Calculate the initial value of each scattering angle based on the incident line and the exit line;

[0024] (b3) Calculate the initial value of the scattering angle and the statistical standard deviation of the set of initial scattering angle values ​​as the scattering angle threshold;

[0025] (b4) Iterate through all initial values ​​of scattering angle and compare them with the scattering angle threshold;

[0026] (b5) After imaging with the typical value of this scattering angle threshold, increase or decrease the scattering angle threshold by 30% and image again. Repeat this several times to select the final imaging result.

[0027] In one embodiment of the present invention, when acquiring and statistically analyzing path length, scattering events, and scattering angle values:

[0028] Allocate a path array L to count and store the path length of each voxel from the incident point to the nearest point, and the path length of each voxel from the nearest point to the exit point when the muon passes through each voxel in the imaging region.

[0029] Assign event array N poca This is used to count and store the number of scattering events of muons passing through each voxel in the imaging region;

[0030] Assign an angle array Θ to count and store scattering angle values.

[0031] In one embodiment of the present invention, during the large voxelization merging process, the path array L and the event array N... poca And the angle array Θ is respectively quantized and merged;

[0032] The first major merging of the path array L is as follows:

[0033] For the entire N x ·N y ·N z The imaging region is of a certain size, covering all M values ​​near each central voxel. x ·M y ·M z Path information L within a small voxel l,m,n The summation is performed according to the following formula and stored in the central voxel L. i,j,k Inside:

[0034]

[0035] and use L i,j,k Replace the original array L l,m,n ,

[0036] After the first major merging of the path array L, a second major merging is performed:

[0037] Separate each dimension and perform large-scale element merging operations sequentially in the x, y, and z directions:

[0038]

[0039] and use Replace the original array L i,j,k .

[0040] In one embodiment of the present invention, the imaging density is the angle array Θ divided by the path array L, and the division of the array is performed element by element.

[0041] In one embodiment of the present invention, based on the event array N poca The number of events recorded within a large voxel is used to remove imaging values ​​from the edges of imaging regions with low event counts, or to set the imaging value to zero.

[0042] The present invention also provides an imaging system for muon scattering, comprising:

[0043] The first set of detectors, located on one side of the imaging object, is used to measure the incident track data of the muon incident imaging region.

[0044] The second set of detectors is used on the other side of the imaged object to measure the exit trajectory data of muons emitted from the imaging region.

[0045] The memory is used to store incident trajectory data, exit trajectory data, incident path array L, and event array N. poca and the data of the angle array Θ;

[0046] The processor receives incident trajectory information and outgoing trajectory information, calculates and statistically analyzes the path length, scattering events, and scattering angle information of muons passing through voxels, and obtains imaging values ​​based on the above information to perform image reconstruction.

[0047] Based on the above technical solution, the technical effects achieved by the present invention are as follows:

[0048] The muon scattering imaging method provided by this invention first subdivides voxels according to rules, calculates and statistically analyzes data such as path length, scattering events, and scattering angle values, and then merges the data in the adjacent grids of the central small voxel into a larger voxel as the imaging value of the central small voxel, thereby increasing the amount of data, reducing data fluctuations, and improving imaging quality.

[0049] The imaging method of the present invention, when reconstructing intravoxel imaging values, combines a scattering angle capping method on the basis of large voxelization and merging to remove initial scattering angle values ​​greater than the scattering angle threshold; that is, in scattering events greater than the scattering angle threshold, their contribution to the imaging value is calculated at most as the contribution of a scattering angle of a certain capping value (scattering angle threshold) to the imaging value, thus reducing the destructive effect of large-angle events on the stability of imaging values. Attached Figure Description

[0050] Figure 1 This is a schematic diagram of the detection system used in the muon scattering imaging method of the present invention.

[0051] Figure 2 This is a flowchart of the muon scattering imaging method of the present invention.

[0052] Figure 3 This is a schematic diagram of the muon passing through the imaging region according to the present invention.

[0053] Figure 4 This is a block diagram of the muon scattering imaging system of the present invention.

[0054] Figure 5 This is a comparison image of U-shaped tungsten blocks imaged using the muon scattering imaging method of this invention and the traditional PoCA method. Detailed Implementation

[0055] To facilitate understanding of the present invention, a more comprehensive description will be given below in conjunction with the accompanying drawings and specific embodiments. The drawings illustrate preferred embodiments of the invention. However, the invention can be implemented in many different forms and is not limited to the embodiments described herein. Rather, these embodiments are provided to provide a thorough and complete understanding of the disclosure of the invention.

[0056] It should be noted that when a component is said to be "fixed to" another component, it can be directly attached to the other component or there may be an intervening component. When a component is said to be "connected to" another component, it can be directly connected to the other component or there may be an intervening component.

[0057] Unless otherwise defined, all technical and scientific terms used herein have the same meaning as commonly understood by one of ordinary skill in the art to which this invention pertains. The terminology used herein in the description of the invention is for the purpose of describing particular embodiments only and is not intended to be limiting of the invention.

[0058] Example 1

[0059] Figure 1 This is a schematic diagram of the detection system used in the muon scattering imaging method of this embodiment. Figure 2This is an overall flowchart of the muon scattering imaging method in this embodiment, in conjunction with reference. Figure 1 and Figure 2 This embodiment provides an imaging method for muon scattering, including the following steps:

[0060] (a) Divide the three-dimensional space of the imaging region traversed by the muon into voxels;

[0061] (b) In each scattering event,

[0062] (1) Obtain the position of the nearest point in the imaging region through which the muon passes.

[0063] (2) Calculate the initial scattering angle based on the incident and exit tracks of the muon through the imaging region; compare the initial scattering angle with the scattering angle threshold:

[0064] When the initial value of the scattering angle is greater than the scattering angle threshold, the scattering angle value is statistically calculated as the square of the scattering angle threshold.

[0065] When the initial scattering angle is not greater than the scattering angle threshold, the scattering angle value is statistically calculated as the square of the initial scattering angle itself.

[0066] (3) Calculate and count the path length of each voxel from the incident point to the nearest point, and calculate and count the path length of each voxel from the nearest point to the exit point.

[0067] (c) Perform large voxelization merging based on the statistical path length and scattering angle values. The large voxelization merging is to merge the voxel data in the grid adjacent to the central small voxel.

[0068] (d) Based on the scattering angle and path length after bulk merging, the imaging density and corresponding imaging values ​​are calculated, and image reconstruction is performed. Specifically, a set of detectors is set on both sides of the imaging object, each set of detectors having two detector units, such as... Figure 1 As shown, two points p0p1 on the two detector units on one side of the imaging object determine the incident line of the muon, and two points q0q1 on the two detector units of the other set of detectors determine the exit line of the muon. This embodiment assumes that the scattering point (least nearest neighbor PoCA) is located at the intersection of the incident and exit lines when they are coplanar, or at the midpoint of the common perpendicular of the two lines when they are skew lines.

[0069] (1) During the imaging process, it is necessary to acquire and statistically analyze parameters such as path length, scattering events, and scattering angle values. Therefore, this embodiment allocates a path array L and an event array N. poca And the angle array Θ, used to store, calculate and statistically analyze the corresponding path length, scattering events and scattering angle values.

[0070] Where the path array L and the angle array Θ are floating-point arrays, and the event array N poca It is an integer array. The path array L is used to store, calculate, and statistically analyze the corresponding path length data, as well as the voxel density information data; the event array N... poca Used to store, calculate, and count the number of scattering events of muons passing through voxels; angle array Θ, used to store, calculate, and count scattering angle values.

[0071] like Figure 2 As shown in the overall flowchart, the three-dimensional space of the imaging object region is first meticulously divided into multiple voxels; specifically, the imaging range in the x-direction (x min ,x max The innermost part is divided into N. x Each part, for example:

[0072]

[0073] The same operation is performed in the y and z directions, dividing the imaging region equally in the x direction, thereby dividing the region into small cuboids, which are called voxels. In some embodiments, the sides of the small cuboids may be equal.

[0074] Within the three-dimensional space of the imaging region, it can be divided into N x ·N y ·N z Individual factors.

[0075] It should be noted that the imaging area of ​​the object being imaged is slightly larger than the set of all voxels. In other words, the imaging area must include the object being measured (imaged) and also leave a certain margin.

[0076] For N in the three-dimensional space of the imaging region x ·N y ·N z For voxels, density information within each voxel can be stored and statistically analyzed using a path array L. This path array L is a floating-point array, and can be either single-precision or double-precision. Specifically, a three-dimensional array L[N] in the programming language can be used directly. x N y N z It can be used for storage and statistics; it can also allocate one-dimensional data, for example, in i+N... x ·j+N x ·N y • k is used as the coordinates of a one-dimensional array to store information about the three-dimensional space voxel (i,j,k).

[0077] This embodiment uses a three-dimensional array as an example, with the coordinates of the path array L starting from 0. After distributing the array, all elements of the path array L are initialized to zero; the event array N... poca And all elements in the angle array Θ are also initialized to zero.

[0078] (2) Next, the Point of Closest Approach (PoCA) when the muon passes through each voxel in the imaging region is obtained and calculated using the following method:

[0079] The incident trajectory of the muon is the incident line containing p0p1, and the exit trajectory is the exit line containing q0q1. Based on the detected incident and exit trajectories, the auxiliary vector is first calculated. and

[0080]

[0081]

[0082]

[0083] in, Let p0 be the vector pointing to p1. Let q0 be the vector pointing to q1. Let p0 be the vector pointing to q0.

[0084] Next, find the point p on the line p0p1 that is closest to q0q1. d , ( express and (inner product),

[0085]

[0086] The point q on q0q1 that is closest to p0p1 on the line is... d ,

[0087]

[0088] Take p d and q d Midpoint between That is, the PoCA point.

[0089]

[0090] Calculate the voxel position coordinates (i,j,k) to which it belongs from the position coordinates (x,y,z) of the PoCA.

[0091]

[0092]

[0093]

[0094] The above is formula group (1).

[0095] The voxel position coordinates (i,j,k) are stored in the path array L.

[0096] (3) After obtaining the nearest neighbor point (PoCA) coordinates of each voxel when the muon passes through the imaging region, and the corresponding voxel coordinates, it is necessary to calculate and count the path length of each voxel segment from the muon incident point to the nearest neighbor point, and to calculate and count the path length of each voxel segment from the nearest neighbor point (PoCA) to the muon exit point. In this embodiment, the muon path through the imaging region is assumed to be a broken line passing through the PoCA.

[0097] Specifically, such as Figure 3 The diagram shows the muons passing through the imaging region:

[0098] (3.1) Obtain the Muon incident line where p0p1 is located. First, it enters the imaging region (x). min ,x max ,y min ,y max ,z min ,z max Incident point p in If the incident line does not pass through the above imaging area, this step and the following steps can be ignored.

[0099] Calculate p in p c The path length of the line segment passing through each voxel in the imaging region is accumulated in group L. in p c The line segment is the point of incidence p of the muon. in to the nearest point p c The line segment in which it is located.

[0100] The specific calculation method is as follows:

[0101] First, label the line p using parametric equations. in p c The point on,

[0102]

[0103] When t = 0, the corresponding starting point is p. in ;

[0104] When t=1, the corresponding endpoint p c.

[0105] Find i = 0, 1, 2, ..., N plane and line p in p c The parameter set of the intersection points And obtain in the same way and And we obtain:

[0106]

[0107] For {t i Sort by parameter t i With parameter t i+1 The line segment formed by the two points represented is the path segment of the muon through each voxel.

[0108] Among them, point The coordinates (x, y, z) are combined with the above formula group (1) to obtain the coordinates (i, j, k) of the corresponding voxel, and the path length is accumulated in the path array L[i, j, k].

[0109] After the above processing, the path length of each voxel from the Muon incident point to the nearest point is calculated and statistically analyzed.

[0110] (3.2) Obtain the Muon exit line where q0q1 is located. First, it passes through the imaging region (x). min ,x max ,y min ,y max ,z min ,z max The exit point p) out If the outgoing line does not pass through the above imaging area, this step and the following steps can be ignored.

[0111] Using (3.1) to calculate p in p c The path length of each segment passing through each voxel in the imaging region is accumulated in group L using the same method to calculate p. c q out The path length of the line segment passing through each voxel in the imaging region is accumulated in group L. c q out The line segment is from the nearest point p c to Muzi's launch point p out The line segment in which it is located.

[0112] (4) Please refer to the description in (2) above regarding "obtaining and calculating the Point of Closest Approach (PoCA) when the muon passes through each voxel in the imaging region":

[0113] The incident trajectory of the muon is the incident line containing p0p1, and the exit trajectory is the exit line containing q0q1. Based on the detected incident and exit trajectories, the auxiliary vector is first calculated. and

[0114]

[0115]

[0116]

[0117] in, Let p0 be the vector pointing to p1. Let q0 be the vector pointing to q1. Let p0 be the vector pointing to q0.

[0118] Calculate the initial value θ of the scattering angle formed by the incident line p0p1 and the outgoing line q0q1, specifically as follows:

[0119]

[0120] The initial scattering angle θ is calculated statistically, and the statistical standard deviation of the set of initial scattering angle values ​​is calculated as a typical value of the scattering angle threshold c (also known as the angle cap value).

[0121] In actual processing, an important reference value for c is a typical value of the initial scattering angle, such as the statistical standard deviation, as processed in this embodiment. If the scattering angle threshold c is too large, the effect of limiting fluctuations will be poor; if the scattering angle threshold c is too small, the distinction between image values ​​will be poor. Therefore, in some embodiments, the selection of the scattering angle threshold c can start from the statistical standard deviation of the initial scattering angle, and several values ​​can be measured around it, such as within a quarter to four times, to obtain the best imaging effect.

[0122] In the angle array Θ[i,j,k], iterate through all initial scattering angle values ​​θ and compare them with the scattering angle threshold c.

[0123] When the initial scattering angle θ is greater than the scattering angle threshold c, the final scattering angle value α entered into the angle array Θ is the square of the scattering angle threshold c, that is, the scattering angle value α = c. 2 The summation and statistics are then stored in the angle array Θ;

[0124] When the initial scattering angle θ is not greater than the scattering angle threshold c, the final scattering angle value α entered into the angle array Θ is the square of the initial scattering angle θ, that is, the scattering angle value α = θ 2 The summation is entered into the angle array Θ.

[0125] The above is the angle capping method used in this embodiment, which can limit the destructive effect of an excessively large θ on the angle array Θ[i,j,k].

[0126] It should be noted that in the above (2) process of "obtaining and calculating the nearest point (PoCA) when the muon passes through each voxel in the imaging region", after calculating the voxel position coordinates (i,j,k) to which the muon belongs from the position coordinates (x,y,z) of the PoCA, the event array N that records the number of events is then... poca [i,j,k] needs to be incremented by 1.

[0127] (5) It should also be noted that the processing of all events and array accumulation can be carried out in parallel. During the processing, this embodiment records the path array L and event array N for each event information in each voxel. poca Both the path array L and the angle array Θ are summed, so the task can be divided into multiple parts, with each subtask assigned an independent path array L and event array N. poca The event is stored in the angle array Θ, and the portion of the event allocated to this subtask is processed accordingly.

[0128] After all events have been processed, set the L and N values ​​for all subtasks. poca The arrays Θ and Θ are merged separately to obtain the final array. This parallel processing method is suitable for parallel processing using multi-threading or multi-processing. For example, in the OpenMP parallel model, there is a reduction statement specifically for handling this type of problem.

[0129] (6) After processing all events, the resulting arrays are then subjected to maximization to obtain the final image.

[0130] (6.1) Perform the first large-scale merging on the path array L. Throughout N... x ·N y ·N z Within the imaging region of this size, all M values ​​near each central voxel are considered. x ·M y ·M z Each small voxel, and the path information L for these small voxels l,m,n The summation is performed according to the following formula and stored in the central voxel L. i,j,k Inside:

[0131]

[0132] and use L i,j,k Replace the original array L l,m,n ,

[0133] (6.2) In this embodiment, after the first large-scale merging of the path array L, a second large-scale merging (optimization) is performed: During the first large-scale merging of the path array L, a total of (N) x ·N y ·N z There are 10 central voxels, and each central voxel needs to interact with its surrounding (M) voxels. x ·M y ·M z The total computational cost for summing the positions within an individual element is as follows: (N) x ·N y ·N z )·(M x ·M y ·M z ).

[0134] In the second large-scale element merging, each dimension is separated, and the large-scale element merging operation is performed sequentially in the x, y, and z directions (the order does not need to be fixed):

[0135]

[0136]

[0137]

[0138] and use Replace the original array L i,j,k .

[0139] (6.3) During the second major merging of the path array, the following is obtained: During the process, the large voxels containing two adjacent small voxels in the x-direction will largely overlap. Therefore, when performing cyclic index increments in the x-direction, in addition to accessing M to calculate the value of the first central voxel, x For each central voxel i+1, the subsequent evaluation only requires adding one element to the tail of the large voxel containing the i-th value and removing one element from the head, for a total of 2 operations.

[0140] Therefore, the number of element visits in this step can be reduced from N. x ·N y ·N z ·M x Reduced to: N y ·N z ·(M x+2·(N x -1)).

[0141] Adding the computational costs of the three steps together, the total computational cost after merging is:

[0142] 6N x N y N z -x2(N x N y +N x N z +N y N z )+M x N y N z +M y N x N z +M z N x N y ,

[0143] In N x >>1,N x >>M x Under the condition that the y and z directions are similar, the computational cost is approximately 6N. x N y N z Compared to the original implementation in 6.2, the computational cost M x M y M z N x N y N z The calculation speed is approximately fast. This represents the computational cost of the large-scale element merging operation after correct optimization following the steps outlined above.

[0144] When performing one-dimensional merging for each dimension as described above, such as merging the x-dimensional element, the elements in the y and z dimensions are independent and should be computed in parallel to improve speed. The large voxel merging process described above is suitable for GPU parallel computing and will provide excellent speed-up when applied to real-time imaging.

[0145] (6.4) It should be noted that for the event array N poca Given the angle array Θ, the first large-scale merging operation (6.1) on the path array L can be performed, and the selected M... x M y M z The size should also be the same as the parameter in (6.1).

[0146] (6.5) Based on the path array and angle array after large-scale merging, calculate the imaging density and corresponding imaging values, and perform image reconstruction.

[0147] Specifically, the imaging density λ = Θ / L, and the angle array and the path array are divided element by element.

[0148] (7) Combine the event array N poca And the angle array Θ, based on the event array N poca The number of events recorded within the current voxel can be removed by discarding the imaging values ​​of edge regions with low event counts, i.e., setting the imaging value to zero.

[0149] The basis for removing the above imaging values ​​is that the structure of the muon scattering imaging device results in a lower number of events detected at the edge of the imaging area than in the center area, and the operator places the object to be imaged in the center area.

[0150] If a few large-angle scattering events occur at the edge (the initial scattering angle is greater than the scattering angle threshold), the resulting image value may be very large, exceeding the maximum value of the image value in the central region. This reduces the contrast of the image in the central region and affects the image quality.

[0151] Therefore, depending on factors such as the total number of events, the size of the imaging area, and the size of the large voxels, the minimum number of events ranges from several to tens or more. A suitable value must meet two criteria: (i) the value should be large enough that there are no obvious bright spots in the imaged area at the edges where there is no material; (ii) the value should not be too large, causing the central imaging area to be removed. Meeting these conditions is sufficient; excessive precision is unnecessary. In some embodiments, for 3D imaging with thousands of events, a value of 5 can be started from and adjusted to meet criterion (i) to see if the imaging is suitable.

[0152] Example 2

[0153] Figure 4 This is a block diagram of the muon scattering imaging system in this embodiment, as shown below. Figure 4 As shown, this embodiment provides an imaging system for muon scattering, including: a first set of detectors located on one side of the imaging object, used to measure incident trajectory data of muons in the imaging region; a second set of detectors located on the other side of the imaging object, used to measure exit trajectory data of muons exiting the imaging region; and a memory for storing incident trajectory data, exit trajectory data, incident path array L, and event array N. poca The processor receives incident trajectory information and outgoing trajectory information, calculates and statistically analyzes the path length, scattering events and scattering angle information of muons passing through voxels, and obtains imaging values ​​based on the above information to perform image reconstruction.

[0154] Example 3

[0155] Using the muon scattering imaging method of Example 1 and the muon scattering imaging system of Example 2, in N x ·N y ·N z To obtain the xy, xz, and yz cross-sectional images within an imaging region of a certain size, the sizes of the large voxels in the z, y, and x directions can be set to be the same as the size of the imaging region, i.e., M. z =N z , or M y =N y M x =N x .

[0156] Then, the calculation is performed according to the method provided in Example 1. After merging, all voxel information in the z (or y, x) direction is the same, which is actually a two-dimensional cross-sectional view.

[0157] Alternatively, the z (or y, x) directions of each array can be summed and merged first to generate a two-dimensional array, and then the large element merging operation can be performed on the two-dimensional array to further save computing resources.

[0158] Figure 5 The images show a comparison between the muon scattering imaging method of Example 1 and the traditional PoCA method for imaging U-shaped tungsten blocks, as shown in the image. Figure 5 As shown, the left image is an image of a U-shaped tungsten block using the traditional PoCA method. Figure 5 A, The right side shows the imaging of a U-shaped tungsten block using the large volumetric image combined with the angle capping method of Example 1. Figure 5 B. It can be seen that the different imaging values ​​of the traditional PoCA method fluctuate greatly with space. Some ultra-high brightness points affect the overall imaging effect. There are also low brightness points in areas with objects and high brightness points in areas without objects, which ultimately leads to unclear object outlines and poor image quality.

[0159] The muon scattering imaging method used in this embodiment can effectively improve the image clarity and reduce the fluctuation of the image value within each pixel and at different spatial locations. At the same time, it reduces the accuracy requirements for the original incident trajectory and the outgoing trajectory, and the image quality is significantly better than that of the traditional PoCA method.

[0160] The above description is merely an example and illustration of the structure of this invention, and while the description is specific and detailed, it should not be construed as limiting the scope of this invention. It should be noted that those skilled in the art can make various modifications and improvements without departing from the concept of this invention, and these obvious substitutions all fall within the protection scope of this invention.

Claims

1. An imaging method for muon scattering, characterized in that, Includes the following steps: (a) Divide the three-dimensional space of the imaging region traversed by the muon into voxels; (b) In each scattering event, (1) Obtain the position of the nearest point through the imaging region where the muon passes. (2) Calculate the initial value of the scattering angle based on the incident and exit tracks of the muon through the imaging region; compare the initial value of the scattering angle with the scattering angle threshold: When the initial value of the scattering angle is greater than the scattering angle threshold, the scattering angle value is statistically calculated as the square of the scattering angle threshold. When the initial scattering angle is not greater than the scattering angle threshold, the scattering angle value is statistically calculated as the square of the initial scattering angle itself. (3) Calculate and count the path length of each voxel from the incident point to the nearest point, and calculate and count the path length of each voxel from the nearest point to the exit point. (c) Perform large voxelization merging based on the statistical path length and scattering angle values, wherein the large voxelization merging is to merge the voxel data in the grid adjacent to the central small voxel; (d) Based on the scattering angle and path length after bulk merging, the imaging density and corresponding imaging values ​​are calculated, and image reconstruction is performed.

2. The imaging method according to claim 1, characterized in that, After obtaining the location of the nearest point, the voxel coordinates of the nearest point are calculated.

3. The imaging method according to claim 1, characterized in that, When dividing the three-dimensional space of the imaging area into voxels, the x, y, and z directions are uniformly divided to form voxels in the form of small cuboids.

4. The imaging method according to claim 1, characterized in that, The typical value of the scattering angle threshold is the statistical standard deviation of the initial scattering angle values ​​collected and calculated.

5. The imaging method according to claim 4, characterized in that, The calculation of the initial value of the scattering angle and the comparison with the scattering angle threshold include the following steps: (b1) In each scattering event, detect the incident line corresponding to the incident trajectory when the muon passes through the voxel, and the exit line corresponding to the exit trajectory; (b2) Calculate the initial value of each scattering angle based on the incident line and the exit line; (b3) Calculate the initial value of the scattering angle and the statistical standard deviation of the set of initial scattering angle values ​​as the scattering angle threshold; (b4) Iterate through all initial values ​​of scattering angle and compare them with the scattering angle threshold; (b5) After imaging with the typical value of this scattering angle threshold, increase or decrease the scattering angle threshold by 30% and image again. Repeat this several times to select the final imaging result.

6. The imaging method according to claim 1, characterized in that, When acquiring and calculating path length, scattering events, and scattering angle values: Allocate a path array L to count and store the path length of each voxel from the incident point to the nearest point, and the path length of each voxel from the nearest point to the exit point when the muon passes through each voxel in the imaging region. Assign event array This is used to count and store the number of scattering events of muons passing through each voxel in the imaging region; Assign an angle array Θ to count and store scattering angle values.

7. The imaging method according to claim 6, characterized in that, During the large-scale merging process, the path array L and the event array And the angle array Θ is respectively quantized and merged; The first major merging of the path array L is as follows: For the whole The imaging region is of a certain size, covering all the areas near each central voxel. Path information L within a small voxel l,m,n The summation is performed according to the following formula and stored in the central voxel L. i,j,k Inside: ， and use L i,j,k Replace the original array L l,m,n , After the first major merging of the path array L, a second major merging is performed: Separate each dimension, in Perform the bulk element merging operation sequentially in the following directions: and use Replace the original array L i,j,k .

8. The imaging method according to claim 7, characterized in that, The imaging density is the angle array Θ divided by the path array L, and the division of the array is performed element by element.

9. The imaging method according to claim 8, characterized in that, Based on the event array The number of events recorded within a large voxel is used to remove imaging values ​​from the edges of imaging regions with low event counts, or to set the imaging value to zero.

10. An imaging system for muon scattering, employing the imaging method for muon scattering as described in any one of claims 1-9, characterized in that, include: The first set of detectors, located on one side of the imaging object, is used to measure the incident track data of the muon incident imaging region. The second set of detectors, located on the other side of the imaged object, is used to measure the exit trajectory data of muons emitted from the imaging region. The memory is used to store incident trajectory data, exit trajectory data, incident path array L, and event array. and the data of the angle array Θ; The processor receives incident trajectory information and outgoing trajectory information, calculates and statistically analyzes the path length, scattering events, and scattering angle information of muons passing through voxels, and obtains imaging values ​​based on the above information to perform image reconstruction.

Images