A cumulative error enhancement descent method for muon transmission tomography

Abstract

The application provides a cumulative error enhanced descent method for muon transmission imaging, comprising the following steps: S1, formulating a muon transmission imaging target, designing a muon transmission imaging model based on block coordinate descent and a block update strategy, and constructing an i-th sub-problem for the k-th iteration in muon transmission imaging; S2, for the i-th sub-problem for the k-th iteration in muon transmission imaging, using an external point penalty function method to constrain the density value of each voxel in each sub-problem, and designing a cumulative error enhanced descent algorithm; S3, under the driving of the muon transmission imaging model and the cumulative error enhanced descent algorithm, completing density structure imaging of the object to be measured. Through the design of the muon transmission imaging model based on block coordinate descent and the block update strategy, the application can convert the solution of a large-scale optimization problem into the solution of a series of small-scale optimization problems, can significantly reduce the occupation of memory resources, and can ensure the accuracy of the imaging result.

Description

Technical Field

[0001] This invention belongs to the field of muon imaging technology, and relates to a muon transmission imaging scenario for large objects under test, and in particular a method for enhancing and reducing cumulative errors in muon transmission imaging. Background Technology

[0002] Muon imaging is a novel imaging technique that uses muons to visualize the internal density and structure of an object under test. Based on its imaging principle, it can be broadly categorized into transmission imaging and scattering imaging. Cosmic muons can penetrate objects several kilometers in size without causing damage. Therefore, compared to traditional X-ray imaging techniques, muon transmission imaging has the advantage of requiring no artificial radiation source and no protective shielding. Consequently, muon imaging has become a research hotspot in the field of geophysical tomography in recent years.

[0003] With the successful development of high-performance muon detectors, muon imaging technology has advanced rapidly. However, current research mainly focuses on nuclear safety and muon scattering imaging, while its application in earth sciences and muon transmission imaging algorithms remains underdeveloped. Geophysical imaging techniques are typically used for large objects. To improve imaging resolution, more X-ray data needs to be acquired and the object divided into the smallest possible voxels, resulting in a large inversion problem. It is worth noting that current research in earth exploration largely focuses on avoiding the calculation and storage of the Hessian matrix to accelerate inversion and reduce memory usage, but it does not consider the interference of environmental noise in X-ray imaging. Therefore, to connect muon transmission imaging technology with real-world applications, researchers in this field should design inversion methods for muon transmission imaging using an engineering-oriented approach. Summary of the Invention

[0004] The purpose of this invention is to address the aforementioned problems in existing technologies by proposing a cumulative error enhancement and reduction method for muon transmission imaging. Based on a muon transmission imaging model and block update strategy using block coordinate descent, the solution of a large-scale optimization problem can be transformed into the solution of a series of small-scale optimization problems, further improving the accuracy of imaging results and reducing memory resource consumption.

[0005] The objective of this invention can be achieved through the following technical solution: a method for enhancing and reducing cumulative errors in muon transmission imaging, comprising the following steps:

[0006] S1. Define the target form for muon transmission imaging as follows:

[0007]

[0008] Where ρ∈R n It is the vector of the prediction model to be solved; G∈R m×n And G ij δ represents the length of the path of the i-th ray within the j-th voxel. i It is an estimate of the measurement error of the i-th ray. Represents the square of the L2 norm of a matrix or vector; It is the model constraint function, matrix W m W x W y and W z It is a discrete representation of the model's smoothness; υ is the regularization factor; ρ ref It is a reference model that can be a zero vector, and is usually set according to the known density information of the object to be measured;

[0009] S2. For the i-th (i = 1, 2, ..., s) subproblem in the k-th iteration of muon transmission imaging, the muon transmission imaging model based on block coordinate descent takes the form of:

[0010]

[0011] Where, ρ i It is the model vector to be solved in the i-th subproblem of the k-th iteration. This represents the solution to the i-th subproblem in the k-th iteration, and is an approximate term. This represents the difference between the current iteration and the previous iteration; the diagonal matrix Λ d This represents the distance in three-dimensional space between the voxels in the current block and the first group of voxels; r i ∈R q It is ρ i Density reference values ​​corresponding to mid-voxels; This represents the difference in voxel density within an abnormal density body, m i ∈R q It is ρ i The average density of the density anomaly containing the mid-voxel; r(ρ i ) is used to compute the separable terms in the muon imaging problem; α s and α o These are weighting coefficients;

[0012] S3. Design a block update strategy based on block coordinate descent to construct the i-th (i=1,2,...,s) subproblem for the k-th iteration in muon transmission imaging;

[0013] S4. For the i-th subproblem in the k-th iteration of muon transmission imaging, the external point penalty function method is adopted to constrain the voxel density value in each subproblem, and the cumulative error enhancement descent algorithm is designed as follows;

[0014] 1) Initialize the i-th subproblem of the k-th iteration of the muon transmission imaging problem:

[0015] First, let v be the maximum number of iterations of the algorithm, q be the number of voxels in the i-th group, p be the number of rays, and l be the upper and lower bounds of the voxels in this group. i u i Secondly, the problem is transformed into an unconstrained optimization problem using an exterior penalty function; finally, the i-th subproblem of the k-th iteration of muon transmission imaging is expressed as a time-varying linear equation system:

[0016] A(t)x(t)-b(t)=0,

[0017] in, T denotes the transpose of the matrix. It is the relationship between voxels and rays in the current subproblem, I∈R q×q It is the identity matrix. It is the smoothing weight of the i-th voxel, Λ s (t)∈R q×q It is a diagonal matrix, when ρ ia When ∈[l,u], Λ s (t) aa =0, otherwise Λ s (t) aa =(λ / (ul)) 0.5 , where ρ ia This represents the a-th voxel in the i-th group; It is the equivalent length of the ray in the current block, ν. s (t)∈R q When ρ ia <l time ν s (t) a =l, when ρ ia >u time ν s (t) a =u, otherwise ν s (t) a =0; x(t)∈R (p +5q)×q Here are the model variables to be solved; t = τΔ, where τ is the update order and τ = 0, 1, 2..., Δ > 0 is the sampling interval, and the system of equations can be discretized as A τ x τ =b τ ;

[0018] 2) Calculate the solution to the i-th muon transmission imaging subproblem in the k-th iteration in the (τ+1)-th iteration:

[0019] Define x τLet be the solution obtained in the τth iteration; therefore, starting from the 1st iteration, based on the solutions obtained in each iteration and the system of linear equations, calculate the solution of the i-th subproblem in the τth iteration; the solution of the i-th subproblem in the k-th iteration in the (τ+1)th iteration is defined as...

[0020]

[0021] in, This represents the pseudo-inverse operation of a matrix, where 0 < η < 1;

[0022] S5. Driven by the muon transmission imaging model and the cumulative error enhancement and descent algorithm, the density structure imaging of the object under test is completed.

[0023] In the above-described cumulative error enhancement and reduction method for muon transmission imaging, in step S3, the block update strategy stipulates that the voxels of the (i+1)th muon transmission imaging subproblem in the k-th iteration satisfy the following conditions: the voxel does not appear in other blocks; the voxel has at least one neighboring voxel in the ith block in three-dimensional space and the neighboring voxel is marked as a density anomalous voxel; the voxel is traversed by at least two rays and these rays do not originate from the same detector.

[0024] In the aforementioned cumulative error enhancement and reduction method for muon transmission imaging, the block update strategy stipulates that the iteration of the current round is terminated when the number of voxels in the i-th muon transmission imaging subproblem is 0.

[0025] In the aforementioned method for enhancing and reducing cumulative errors in muon transmission imaging, in step S4, the external point penalty function is defined as...

[0026]

[0027] Wherein, the penalty factor λ > 0; ρ represents the voxel density, and l and u are the minimum and maximum values ​​that it is allowed to take, respectively.

[0028] In the above-described method for enhancing and reducing cumulative errors in muon transmission imaging, step S4 determines whether the algorithm has terminated and outputs the solution to the i-th subproblem of muon transmission imaging:

[0029] The termination threshold of the algorithm is defined as μ∈R, which is used to determine the convergence of the cumulative error enhancement descent algorithm for the i-th subproblem in the k-th iteration; in the i-th muon transmission imaging subproblem in the k-th iteration, the L2 norms of the errors in the τ-th and τ+1-th iterations are respectively e τ =||A τ x τ -b τ ||2 and e τ+1 =||A τ+1 x τ+1-b τ+1 ||2; The termination criterion for the cumulative error augmentation descent algorithm is designed to determine |e τ -e τ+1 |Is it less than or equal to the termination threshold μ? If the condition is met, terminate the iteration and return the result x after the (τ+1)th iteration. τ+1 If the condition is not met, return to the solution calculated in iteration τ+1, until the condition is met or the algorithm exceeds the maximum number of iterations v, at which point x is returned accordingly when the condition is met. τ+1 Or x obtained after the vth iteration v Finally, update based on the returned results.

[0030] In the aforementioned method for enhancing and reducing cumulative errors in muon transmission imaging, step S5 involves obtaining the solution to the subproblem of muon transmission imaging. The problem is transformed into a solution to the original problem of muon transmission imaging, thereby completing the density structure imaging of the object under test.

[0031] Compared with existing technologies, the cumulative error enhancement and reduction method for muon transmission imaging presented here has the following advantages:

[0032] 1. This invention aims to reduce memory usage and mitigate noise interference with imaging results by designing a muon transmission imaging system. Through the design of a muon transmission imaging model based on block coordinate descent and a block update strategy, the solution of a large-scale optimization problem can be transformed into the solution of a series of small-scale optimization problems, significantly reducing memory resource consumption.

[0033] 2. Furthermore, for imaging scenarios with upper and lower bound constraints, an external point penalty function is used to constrain the density values ​​of voxels in each subproblem. The specific implementation method and the setting of its solution algorithm effectively realize the density structure imaging of the object under test, ensuring the accuracy of the imaging results. Attached Figure Description

[0034] Figure 1 This is a schematic diagram showing the position and size of a standard rock block in three-dimensional imaging space in an embodiment of the present invention;

[0035] Figure 2 This is a schematic diagram showing the position and size of the brick, clay, and plexiglass blocks in the three-dimensional imaging space in an embodiment of the present invention;

[0036] Figure 3 This is a schematic diagram showing the position of the muon detector in three-dimensional space in an embodiment of the present invention;

[0037] Figure 4 This is a slice image of the imaging result of the object under test in the east-west direction in an embodiment of the present invention;

[0038] Figure 5 This is a line graph showing the average density of the imaging results and theoretical results of the object under test in the east-west direction in an embodiment of the present invention.

[0039] Figure 6 The above figures show the simulation results of precision and recall for each iteration in this embodiment of the invention.

[0040] Figure 7 This is a slice image of the imaging result of the object under test in the vertical direction under noise interference in an embodiment of the present invention. Detailed Implementation

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

[0042] like Figures 1 to 7 As shown, a specific embodiment of the present invention is as follows:

[0043] Suppose a muon transmission imaging scene has a density of 2.65 × 10⁻⁶. 3 kg / m 3 And its volume is 15×20×10m 3 Standard rock blocks, such as Figure 2 As shown; it contains three volumes, each 3×4×2m. 3 And the densities are 1.80×10 3 kg / m 3 The bricks, 1.50×10 3 kg / m 3 clay and 1.20×10 3 kg / m 3 Pieces of plexiglass, such as Figure 2 As shown. The simulation data was obtained through forward modeling, which simulates the distribution of incident energy and incident direction of natural muons at sea level. The noise level is comparable to the theoretical noise level in the real environment. This implementation used four muon detectors, such as... Figure 3 As shown, a total of 2.35 × 10⁻⁶ samples were collected. 4 The effective X-ray information was obtained, and the imaging area was non-uniformly meshed into a 6.63×10 grid. 5 Individual elements. The coefficient matrix G and vector d, as well as the error σ of the measurement data, are derived from forward modeling; parameter α s =1; parameter α o =0.05; Based on the characteristics of the imaging space, a reference model ρ is designed. ref The corresponding constraints l and u are as follows: When voxel j is air, l j =-0.10, u j =0.10, when voxel j is the object to be measured. l j =0, u j =2.70.

[0044] Considering that current research in the field of Earth exploration mostly focuses on avoiding the calculation and storage of the Hessian matrix to accelerate inversion and reduce memory usage, but does not consider the interference of environmental noise in X-ray imaging technology on the imaging results, this embodiment provides "a method for enhancing and reducing cumulative errors in muon transmission imaging" to reduce memory resource consumption and mitigate the impact of environmental noise on imaging results. This method includes the following steps:

[0045] Step 1: Based on the known density information of the object to be measured, select several voxels as the first block, and then design the solution model for each block and formulate the block update strategy.

[0046] Specifically:

[0047] Set the initial density ρ of the voxels. 0 =ρ ref The i-th group has q voxels, and its initial density is expressed as... The number of rays is p; the model for the i-th subproblem in the k-th iteration is defined as follows:

[0048]

[0049] Where, ρ i It is the model vector to be solved in the i-th subproblem of the k-th iteration. This represents the solution to the i-th subproblem in the k-th iteration, and is an approximate term. This represents the difference between the current iteration and the previous iteration; the diagonal matrix Λ d This represents the distance in three-dimensional space between the voxels in the current block and the first group of voxels; r i ∈R q It is ρ i Density reference values ​​corresponding to mid-voxels; This represents the difference in voxel density within an abnormal density body, m i ∈R q It is ρ i The average density of the density anomaly containing the mid-voxel; r(ρ i ) is used to compute the separable terms in the muon imaging problem; α s =1 and α o=0.05 is the weighting coefficient. Furthermore, the block update strategy stipulates that the voxels of the (i+1)th muon transmission imaging subproblem in the k-th iteration must satisfy the following conditions: the voxel has not appeared in any other block; the voxel has at least one neighboring voxel in the 3D space located in the ith block, and this neighboring voxel is marked as a density anomalous voxel; the voxel is traversed by at least two rays, and these rays do not originate from the same detector. Simultaneously, this strategy stipulates that the iteration terminates when the number of voxels in the ith muon transmission imaging subproblem is 0.

[0050] Step 2: For the i-th subproblem in the k-th iteration of muon transmission imaging, an exterior penalty function is used to constrain the density values ​​of voxels in each subproblem, and a solution algorithm based on block coordinate descent-error accumulation is designed.

[0051] Specifically:

[0052] S201: Initializing the i-th subproblem of the k-th iteration of the muon transmission imaging problem:

[0053] First, the maximum number of iterations of the algorithm is set to v = 30, the number of voxels in the i-th group is q, and a total of p rays pass through the voxels in this group. The upper and lower bounds of the voxels in this group are l and l, respectively. i with u i Secondly, the problem is transformed into an unconstrained optimization problem by incorporating an exterior penalty function; finally, the i-th subproblem of the k-th iteration of muon transmission imaging is expressed as a time-varying linear equation system:

[0054] A(t)x(t)-b(t)=0,

[0055] in, T denotes the transpose of the matrix. It is the relationship between voxels and rays in the current subproblem, I∈R q×q It is the identity matrix. It is the smoothing weight of the i-th voxel, Λ s (t)∈R q×q It is a diagonal matrix, when ρ ia When ∈[l,u], Λ s (t) aa =0, otherwise Λ s (t) aa =(λ / (ul)) 0.5 , where ρ ia This represents the a-th voxel in the i-th group; It is the equivalent length of the ray in the current block, ν. s (t)∈R q When ρ ia <l time ν s (t) a=l, when ρ ia >u time ν s (t) a =u, otherwise ν s (t) a =0; x(t)∈R (p +5q)×q Here are the model variables to be solved; t = τΔ, where τ is the update order and τ = 0, 1, 2..., Δ > 0 is the sampling interval, and the system of equations can be discretized as A τ x τ =b τ .

[0056] S202: Calculate the solution to the i-th muon transmission imaging subproblem in the k-th iteration in the (τ+1)-th iteration:

[0057] Define x τ Let be the solution obtained in the τth iteration; therefore, starting from the 1st iteration, based on the solutions obtained in each iteration and the system of linear equations, the solution of the i-th subproblem in the τth iteration is calculated; specifically, the solution of the i-th subproblem in the k-th iteration in the τ+1th iteration is defined as

[0058]

[0059] in, This represents the pseudo-inverse operation of a matrix, where η = 0.05. This represents the noise component to simulate additional noise interference. In this implementation... Simulated random noise or Simulate linear noise.

[0060] S203: Determine if the algorithm has terminated, and output the solution to the i-th muon transmission imaging subproblem:

[0061] The termination threshold of the algorithm is defined as μ = 0.001, which can be used to determine the convergence of the cumulative error enhancement descent algorithm for the i-th subproblem in the k-th iteration. In the i-th muon transmission imaging subproblem in the k-th iteration, the L2 norms of the errors in the τ-th and τ+1-th iterations are respectively e τ =||A τ x τ -b τ ||2 and e τ+1 =||A τ+1 x τ+1 -b τ+1 ||2; The algorithm's termination criterion is designed to determine |e τ -e τ+1 |Is it less than or equal to the termination threshold μ? If the condition is met, terminate the iteration and return the result x after the (τ+1)th iteration. τ+1If the condition is not met, return to step S202 until the condition is met or until the algorithm exceeds the maximum number of iterations v, at which point return x when the condition is met. τ+1 Or x obtained after the vth iteration v Finally, update based on the returned results.

[0062] In this embodiment, Python was used for execution. After a simulation of 5 minutes and 15 seconds, the simulation results are as follows: Figure 4 , Figure 5 , Figure 6 , Figure 7 As shown. Figure 4 The east-west slice of the prediction model shown clearly distinguishes the boundaries and locations of the three density anomalies; Figure 5 The average density of the east-west slice of the object under test shown reflects that the density distribution of the predicted model is very close to that of the theoretical model; Figure 6 In the final imaging results shown, the recall and precision of the prediction model are stable at around 98.96% and 90.74%, respectively. Of course, the accuracy of this invention can be further improved by adjusting the design parameters, etc. Figure 7 The vertical slices shown in the imaging results under noise interference are still very clear, indicating that the invention has good robustness.

[0063] In summary, this invention effectively achieves density structure imaging of a target object using muons, ensuring the accuracy of the imaging results and reducing the consumption of memory resources.

[0064] Finally, it should be noted that the above preferred embodiments are only used to illustrate the technical solutions of the present invention and are not intended to limit it. Although the present invention has been described in detail through the above preferred embodiments, those skilled in the art should understand that various changes can be made to it in form and detail without departing from the scope defined by the claims of the present invention.

Claims

1. A method for enhancing and reducing cumulative errors in muon transmission imaging, characterized in that, Includes the following steps: S1. Define the target form for muon transmission imaging as follows: Where ρ∈R n It is the vector of the prediction model to be solved; G∈R m×n And G ij δ represents the length of the path of the i-th ray within the j-th voxel. i It is an estimate of the measurement error of the i-th ray. Represents the square of the L2 norm of a matrix or vector; It is the model constraint function, matrix W m W x W y and W z It is a discrete representation of the model's smoothness; υ is the regularization factor; ρ ref It is a reference model that can be a zero vector, and is usually set according to the known density information of the object to be measured; S2. For the i-th subproblem in the k-th iteration of muon transmission imaging, the muon transmission imaging model based on block coordinate descent takes the form of: Where, ρ i ρ is the model vector to be solved in the i-th subproblem of the k-th iteration. i k This represents the solution to the i-th subproblem in the k-th iteration, and is an approximate term. This represents the difference between the current iteration and the previous iteration; the diagonal matrix Λ d This represents the distance in three-dimensional space between the voxels in the current block and the first group of voxels; r i ∈R q It is ρ i Density reference values ​​corresponding to mid-voxels; This represents the difference in voxel density within an abnormal density body, m i ∈R q It is ρ i The average density of the density anomaly containing the mid-voxel; r(ρ i ) is used to compute the separable terms in the muon imaging problem; α s and α o These are weighting coefficients; S3. Design a block update strategy based on block coordinate descent to construct the i-th subproblem for the k-th iteration in muon transmission imaging; S4. For the i-th subproblem in the k-th iteration of muon transmission imaging, the external point penalty function method is adopted to constrain the voxel density value in each subproblem, and the cumulative error enhancement descent algorithm is designed as follows; 1) Initialize the i-th subproblem of the k-th iteration of the muon transmission imaging problem: First, let v be the maximum number of iterations of the algorithm, q be the number of voxels in the i-th group, p be the number of rays, and l be the upper and lower bounds of the voxels in this group. i u i Secondly, the problem is transformed into an unconstrained optimization problem using an exterior penalty function; finally, the i-th subproblem of the k-th iteration of muon transmission imaging is expressed as a time-varying linear equation system: A(t)x(t)-b(t)=0, in, T denotes the transpose of the matrix. It is the relationship between voxels and rays in the current subproblem, I∈R q×q It is the identity matrix. It is the smoothing weight of the i-th voxel, Λ s (t)∈R q×q It is a diagonal matrix, when ρ ia When ∈[l,u], Λ s (t) aa =0, otherwise Λ s (t) aa =(λ / (ul)) 0.5 , where ρ ia This represents the a-th voxel in the i-th group; It is the equivalent length of the ray in the current block, ν. s (t)∈R q When ρ ia <l time ν s (t) a =l, when ρ ia >u time ν s (t) a =u, otherwise ν s (t) a =0; x(t)∈R (p +5q)×q Here are the model variables to be solved; t = τΔ, where τ is the update order and τ = 0, 1, 2..., Δ > 0 is the sampling interval, and the system of equations can be discretized as A τ x τ =b τ ; 2) Calculate the solution to the i-th muon transmission imaging subproblem in the k-th iteration in the (τ+1)-th iteration: Define x τ Let be the solution obtained in the τth iteration; therefore, starting from the 1st iteration, based on the solutions obtained in each iteration and the system of linear equations, calculate the solution of the i-th subproblem in the τth iteration; the solution of the i-th subproblem in the k-th iteration in the (τ+1)th iteration is defined as... in, This represents the pseudo-inverse operation of a matrix, where 0 < η < 1; S5. Driven by the muon transmission imaging model and the cumulative error enhancement and descent algorithm, the density structure imaging of the object under test is completed.

2. The cumulative error enhancement and reduction method for muon transmission imaging as described in claim 1, characterized in that, In step S3, the block update strategy stipulates that the voxels of the (i+1)th muon transmission imaging subproblem in the kth iteration must satisfy the following: the voxel does not appear in other blocks; the voxel has at least one neighboring voxel in the 3D space in the i-th block and the neighboring voxel is marked as a density anomalous voxel; the voxel is traversed by at least two rays and these rays do not originate from the same detector.

3. The cumulative error enhancement and reduction method for muon transmission imaging as described in claim 2, characterized in that, The update strategy for the block stipulates that the iteration of the current round is terminated when the number of voxels in the i-th muon transmission imaging subproblem is 0.

4. The cumulative error enhancement and reduction method for muon transmission imaging as described in claim 1, characterized in that, In step S4, the exterior penalty function is defined as Wherein, the penalty factor λ > 0; ρ represents the voxel density, and l and u are the minimum and maximum values ​​that it is allowed to take, respectively.

5. The cumulative error enhancement and reduction method for muon transmission imaging as described in claim 4, characterized in that, In step S4, determine whether the cumulative error enhancement descent algorithm has terminated, and output the solution to the i-th muon transmission imaging subproblem: Define the termination threshold of the algorithm as μ∈R, which is used to determine the convergence of the cumulative error enhancement descent algorithm for the i-th subproblem in the k-th iteration; in the i-th muon transmission imaging subproblem in the k-th iteration, the L2 norms of the errors in the τ-th and τ+1-th iterations are respectively e τ =||A τ x τ -b τ ||2 and e τ+1 =||A τ+1 x τ+1 -b τ+1 ||2; The termination criterion for the cumulative error augmentation descent algorithm is designed to determine |e τ -e τ+1 |Is it less than or equal to the termination threshold μ? If the condition is met, terminate the iteration and return the result x after the (τ+1)th iteration. τ+1 If the condition is not met, return to the solution calculated in iteration τ+1, until the condition is met or the algorithm exceeds the maximum number of iterations v, at which point x is returned accordingly when the condition is met. τ+1 Or x obtained after the vth iteration v Finally, update based on the returned results.

6. The cumulative error enhancement and reduction method for muon transmission imaging as described in claim 5, characterized in that, In step S5, the solution to the obtained muon transmission imaging subproblem is... The problem is transformed into a solution to the original problem of muon transmission imaging, thereby completing the density structure imaging of the object under test.

Images