A scatter ray calculation method and system based on Monte Carlo perturbation calculation

Abstract

The present invention discloses a scattered ray calculation method and system based on Monte Carlo perturbation calculation. The method includes the following steps: (1) establishing a Monte Carlo radiation transport calculation model; (2) sampling source particles; (3) calculating the offset parameters of source particles; (4) performing transport simulation on the source particles; (5) calculating the scattered ray count contributions of the current simulated source particles and their offsets to different sub-beams to a two-dimensional plane measurement device; (6) calculating the scattered ray contributions of a unit intensity radiation source. The system includes the following modules: a source sampling module, a transport simulation module, a data processing module, and a storage medium of the system. The present invention can effectively improve the efficiency of calculating scattered rays by the Monte Carlo method and achieve high-efficiency and accurate scattered ray calculation.

Description

Technical Field

[0001] The present invention relates to the application fields of nuclear energy and nuclear technology such as radiotherapy, industrial nondestructive testing, etc., and particularly to a scattered ray calculation method and system based on Monte Carlo perturbation calculation. Background Art

[0002] In the fields of radiotherapy dose monitoring, medical image reconstruction, industrial nondestructive monitoring, reactor monitoring, etc., generally, the source strength distribution, the flux / dose distribution in the medium, or the medium composition is reconstructed based on the transmission distribution after penetrating the medium. The transmission distribution formed after the ray penetrates the medium includes the "primary ray" that does not react with the medium and the "scattered ray" that reacts with the medium. If the "scattered ray" can be accurately removed, then the source strength, the flux / dose distribution in the medium, or the medium composition can be quickly and accurately obtained through the "primary ray" distribution. Therefore, studying an accurate scattered ray calculation method is of great significance to the above fields.

[0003] Currently, due to speed requirements, the scattered ray calculation mainly uses analytical methods. By constructing an analytical scattering kernel, the scattering distribution after the incident ray passes through the medium is obtained, or a simplified Monte Carlo method is used to calculate the scattering distribution. These methods have a relatively fast calculation speed. However, due to many approximate assumptions, it is difficult to give an accurate scattered ray distribution when the ray has complex scattering effects with the medium, resulting in insufficient final reconstruction accuracy. Although the Monte Carlo method has high accuracy, in order to obtain an accurate scattered ray distribution, it often requires a large amount of calculation time. Existing research has used the method of track reuse to accelerate the Monte Carlo simulation speed. However, due to the need to store a large number of tracks, the memory consumption is huge. Therefore, the Monte Carlo method is rarely used in scattered ray calculation. Summary of the Invention

[0004] Object of the Invention: The object of the present invention is to provide an accurate and efficient scattered ray calculation method and system based on Monte Carlo perturbation calculation.

[0005] Technical Solution: A scattered ray calculation method based on Monte Carlo perturbation calculation according to the present invention includes the following steps:

[0006] (1) Establishing a Monte Carlo radiation transport calculation model: Establish a Monte Carlo radiation transport calculation model according to the source distribution information, the medium image information, and the two-dimensional plane measurement device information;

[0007] (2) Sampling source particles: Divide the entire irradiation field into n sub-beams according to sub-fields of a certain size. Sample the source particle parameters (r i , Ω i ) according to the source distribution information, and confirm the sub-beam i to which the source particle belongs according to the source particle position;

[0008] (3) Source particle offset parameter calculation: Shift the source particle to another sub-beam j, and calculate the source offset correction parameter Bis_S based on the current sub-beam source distribution information and the source distribution information of the sub-beam after offset. i,j = P j (r j , Ω j ) / P i (r i , Ω i ), where P i (r i , Ω i ) is the generation probability of the current source particle in sub-beam i, and P j (r j , Ω j ) is the generation probability of the offset source particle after being shifted to sub-beam j;

[0009] (4) Conduct transport simulation on the source particle: When the particle takes each step forward, shift the current particle to other sub-beams, and calculate the offset correction parameter for the current step based on the perturbation calculation method;

[0010] (5) Calculate the scattered ray count contribution of the current simulated source particle and its shifted to different sub-beams to the two-dimensional plane measurement device: When the source particle transport simulation reaches the point (x, y) of the two-dimensional measurement plane device, count the scattered ray contribution Flux i (x, y) of the particle to the two-dimensional plane measurement device. At the same time, obtain the total offset correction parameter of the particle shifted to different sub-beams according to the offset parameters calculated in steps (3) and (4) where is the offset correction parameter when the particle is transported in the k-th step from sub-beam i to sub-beam j. Finally, statistically obtain the scattered ray contribution Flux j (x, y) = Flux i (x, y) * Bis_Tot i,j ;

[0011] (6) Calculate the scattered ray contribution of the unit intensity radiation source: After simulating all source particles, statistically obtain and normalize the scattered ray contributions of each sub-beam to the two-dimensional plane device, and then multiply by the number of particles emitted by the unit intensity radiation source to calculate the scattered ray contribution value of all sub-beams under unit intensity.

[0012] The specific content of step (3) is as follows:

[0013] (3.1) According to the source particle parameters (r i , Ω i ) obtained by sampling and the sub-beam to which the source particle belongs, calculate the generation probability P of the source particle in the corresponding sub-beami (r i , Ω i ), where i is the sub - beam number to which the source particle belongs;

[0014] (3.2) Shift the source particle space to other sub - beams to obtain new source particle parameters (r j , Ω j ). Calculate the production probability P j (r j , Ω j ) of the source particle in the current sub - beam, where j is the sub - beam number to which the source particle is shifted;

[0015] (3.3) Based on the above two sets of probabilities, calculate the offset correction parameter Bis_S when the source particle in the i - th sub - beam is shifted to the j - th sub - beam i,j = P j (r j , Ω j ) / P i (r i , Ω i ).

[0016] The specific steps of step (4) are as follows:

[0017] (4.1) Transport the source particle in the i - th sub - beam to the k - th step. The starting point parameters are (r i,k,start ), and the ending point parameters are (r i,k,end ). Among them, i is the sub - beam number where the current particle is located, and k is the number of the transport step of the current particle. Shift the starting point of the k - th step to other sub - beams to obtain the starting point coordinates (r j,k,start ) and ((r j,k,end )) of the k - th step under different sub - beams, where j is the sub - beam number to which the particle is shifted;

[0018] (4.2) According to the starting point coordinates of the particle transport in the k - th step under different sub - beams, combined with the transport medium image data, calculate the equivalent paths and that the particle in the current sub - beam i and the particle shifted to different sub - beams j pass through in the k - th step of transport. At the same time, obtain the reaction cross - section and of the particle at the k - th step according to the material of the grid element where the starting point coordinates are located. Calculate the transport offset correction parameter when the source particle in the i - th sub - beam is transported to the k - th step and shifted to the j - th sub - beam

[0019] (4.3) According to the ending point coordinates of the particle transport in the k - th step under different sub - beams, combined with the transport medium image data, calculate the collision offset correction parameter when the reaction of the current particle after the k - th step of transport is shifted from sub - beam i to sub - beam j Among them, m is the nuclide number of the collision that occurs after the k-th step of the actual transport of the particle, is the nucleon density of the nuclide of the collision that occurs after the k-th step of the actual transport of the particle, is the nucleon density of the m-th nuclide at the position where the corresponding end coordinate is located after the particle transport deviates to the j-th sub-beam in the k-th step;

[0020] (4.4) Calculate the offset correction parameter of the current step according to the transport offset correction parameter and the collision offset correction parameter calculated in steps (4.2) and (4.3) The source distribution information described in step (1) includes the open field size, source-skin distance, source-axis distance, particle type, and accelerator energy.

[0021] The medium image information described in step (1) includes CT image resolution, size, and CT value.

[0022] The two-dimensional plane measurement device information described in step (1) includes the distance from the source to the two-dimensional plane measurement device.

[0023] A scattered ray calculation system based on Monte Carlo perturbation calculation, characterized in that the system includes the following modules:

[0024] Source sampling module: This module samples source particle parameters according to source distribution information and calculates source particle offset parameters;

[0025] Transport simulation module: This module realizes the transport process of particles in the medium and counts the contributions of scattered rays of the current particle and the offset particle to the counting area;

[0026] Data processing module: This module reads in source distribution information, medium image information, and two-dimensional plane measurement device information and converts them into a Monte Carlo radiation transport calculation model, and at the same time reads in reaction cross-section data involved in the particle transport process;

[0027] System storage medium: This module stores the Monte Carlo radiation transport calculation model containing source distribution information, medium image information, and two-dimensional plane measurement device information, reaction cross-section data involved in the particle transport process, and transport track information, offset parameters, and counting information during the particle transport process.

[0028] A computer storage medium, on which a computer program is stored, and when the computer program is executed by a processor, it implements the above-mentioned scattered ray calculation method based on Monte Carlo perturbation calculation.

[0029] A computer device, including a storage, a processor, and a computer program stored on the storage and executable on the processor, and when the processor executes the computer program, it implements the above-mentioned scattered ray calculation method based on Monte Carlo perturbation calculation.

[0030] Beneficial effects: Compared with the prior art, the present invention has the following advantages:

[0031] 1. Based on the Monte Carlo method for scatter ray calculation, accurate scatter ray contributions can be given, thereby obtaining an accurate primary ray distribution and achieving high-precision reconstruction of source strength / medium dose / medium composition;

[0032] 2. During particle transport, the simulated contribution of each particle is offset to different source positions through perturbation calculation, thereby effectively improving the Monte Carlo simulation efficiency and achieving efficient scatter ray calculation. At the same time, the memory consumption will not increase significantly during the transport process;

[0033] 3. It can effectively improve the reconstruction accuracy of source strength / medium composition / dose flux in the medium, etc. in applications such as industrial inspection, radiotherapy dose monitoring, and medical image reconstruction. Description of the Drawings

[0034] Figure 1 It is a flowchart of a scatter ray calculation method based on Monte Carlo perturbation calculation. Detailed Embodiments

[0035] The technical solution of the present invention will be further described below in conjunction with the drawings.

[0036] Example 1: As Figure 1 shown, dose monitoring and verification are performed on a typical nasopharyngeal carcinoma patient, and the scatter ray distribution on a two-dimensional planar device is calculated. The current field size of this patient is 10 cm × 8 cm, the particle type is photon, the accelerator energy is 6 MV, the source-skin distance is 95 cm, the source-axis distance is 100 cm, and the distance from the source to the two-dimensional planar measurement device is 160 cm. The specific steps are as follows:

[0037] (1) Establishment of Monte Carlo radiation transport calculation model

[0038] According to the source distribution information of this example, including a field size of 10 cm × 8 cm, a source-skin distance of 95 cm, a source-axis distance of 100 cm, a particle type of photon, an accelerator energy of 6 MV, CT image information, including patient CT image resolution, size, CT value, etc., and two-dimensional planar measurement device information, including a distance of 160 cm from the source to the two-dimensional planar measurement device, etc., a Monte Carlo radiation transport calculation model is established.

[0039] (2) Sampling of source particles

[0040] The entire irradiation field is divided into 320 sub-beams according to sub-fields of a certain 0.5 cm × 0.5 cm, and source particle parameters (r i ,Ω i), confirm the sub-beam to which the source particle belongs according to the position of the source particle, and assign the number i to the sub-beam to which the source particle belongs.

[0041] (3) Calculation of source particle offset parameters

[0042] (3.1) According to the source particle parameters (r i , Ω i ) obtained by sampling and the sub-beam to which the source particle belongs, calculate the generation probability P i (r i , Ω i ) of the source particle in the corresponding sub-beam, where i is the number of the sub-beam to which the source particle belongs;

[0043] (3.2) Spatially offset the source particle to other sub-beams to obtain new source particle parameters (r j , Ω j ), and calculate the generation probability P j (r j , Ω j ) of the source particle in the current sub-beam, where j is the number of the sub-beam to which the source particle is offset;

[0044] (3.3) Based on the above two sets of probabilities, calculate the offset correction parameter Bis_S when the source particle in the i-th sub-beam is offset to the j-th sub-beam i,j = P j (r j , Ω j ) / P i (r i , Ω i );

[0045] (4) Perform transport simulation on the source particle

[0046] (4.1) The source particle in the i-th sub-beam is transported to the k-th step, the starting point parameters are (r i,k,start ), and the ending point parameters are (r i,k,end ), where i is the number of the sub-beam where the current particle is located, and k is the number of the transport step of the current particle. Offset the starting point of the k-th step to other sub-beams to obtain the starting point coordinates (r j,k,start ) and ((r j,k,end )) of the k-th step under different sub-beams, where j is the number of the sub-beam to which the particle is offset;

[0047] (4.2) According to the starting point coordinates of the particle transport at the k-th step under different sub-beams, combined with the image data of the transport medium, calculate the equivalent paths and that the particle in the current sub-beam i and the particle offset to different sub-beams j pass through at the k-th step of transport. At the same time, obtain the reaction cross-section of the particle at the k-th step according to the material of the grid element where the starting point coordinates are located and Calculate the transport offset correction parameter when the particle of the i-th sub-beam source is transported to the j-th sub-beam at the k-th step

[0048] (4.3) According to the end coordinates of the particle transport at the k-th step under different sub-beams, combined with the image data of the transport medium, calculate the collision offset correction parameter when the reaction after the k-th step of the current particle is transported from sub-beam i to sub-beam j Where m is the nuclide number of the collision that occurs after the k-th step of the particle's real transport, is the nucleon density of the nuclide that undergoes a collision after the k-th step of the particle's real transport, is the nucleon density of the m-th nuclide at the position corresponding to the end coordinate after the particle is transported to the j-th sub-beam at the k-th step;

[0049] (4.4) Calculate the offset correction parameter of the current step according to the transport offset correction parameter and the collision offset correction parameter calculated in steps (4.2) and (4.3)

[0050] (5) Calculate the scattered ray count contribution of the current simulated source particle and its offset to different sub-beams to the two-dimensional plane measurement device

[0051] (5.1) Calculate the total correction factor of the i-th sub-beam source particle transported to the j-th sub-beam according to steps (3) and (4)

[0052] (5.2) Calculate the scattered ray contribution Flux of the j-th sub-beam to the two-dimensional plane device when the i-th sub-beam source particle is transported to the j-th sub-beam j (x,y) = Flux i (x,y) * Bis_Tot i,j , where Flux i (x,y) is the scattered ray contribution of the i-th sub-beam source particle to the two-dimensional plane device after transport, and (x,y) are the coordinates of the detection point on the two-dimensional plane measurement device

[0053] (6) Calculate the scattered ray contribution of the unit intensity radiation source

[0054] After simulating all source particles, count and normalize the scattered ray contributions of each sub-beam to the two-dimensional plane device, and then multiply by the number of particles emitted by the unit intensity radiation source to calculate the scattered ray contribution values of all sub-beams under unit intensity

[0055] Example 2: A scattered ray calculation system based on Monte Carlo perturbation calculation, the system includes the following modules:

[0056] Source sampling module: This module samples source particle parameters according to source distribution information and calculates source particle offset parameters.

[0057] Transport simulation module: This module implements the transport process of particles in the medium and statistically calculates the scattered ray contributions of the current particles and the offset particles to the counting region.

[0058] Data processing module: This module reads in source distribution information, medium image information, and two-dimensional plane measurement device information and converts them into a Monte Carlo radiation transport calculation model. At the same time, it reads in reaction cross-section data involved in the particle transport process.

[0059] System storage medium: This module stores the Monte Carlo radiation transport calculation model containing source distribution information, medium image information, and two-dimensional plane measurement device information, reaction cross-section data involved in the particle transport process, as well as transport track information, offset parameters, and counting information during the particle transport process.

[0060] Embodiment 3: A computer storage medium having a computer program stored thereon, which when executed by a processor implements the above-described scattered ray calculation method based on Monte Carlo perturbation calculation.

[0061] Embodiment 4: A computer device includes a storage, a processor, and a computer program stored on the storage and executable on the processor. When the processor executes the computer program, it implements the above-described scattered ray calculation method based on Monte Carlo perturbation calculation.

Claims

1. A scatter ray calculation method based on Monte Carlo perturbation calculation, characterized in that It includes the following steps: (1) Establishing a Monte Carlo radiation transport calculation model: Establish a Monte Carlo radiation transport calculation model based on source distribution information, medium image information, and two-dimensional plane measurement device information; (2) Source particle sampling: Divide the entire irradiation field into n sub-beams according to sub-fields of a certain size, and sample the source particle parameters (r i0 , Ω i0 ) according to the source distribution information, and confirm the sub-beam i0 to which the source particle belongs based on the source particle position; (3) Source particle offset parameter calculation: Offset the source particle to another sub-beam j0, and calculate the source offset correction parameter Bis_S according to the current sub-beam source distribution information and the source distribution information of the sub-beam after offset i0,j0 = P j0 (r j0 , Ω j0 ) / P i0 (r i0 , Ω i0 ), where P i0 (r i0 , Ω i0 ) is the generation probability of the current source particle in the sub-beam i0, and P j0 (r j0 , Ω j0 ) is the generation probability of the offset source particle after offset to the sub-beam j0; (4) Performing transport simulation on the source particles: When the particle takes each step forward, shift the current particle to other sub-beams, and calculate the offset correction parameter for the current step based on the perturbation calculation method; (4.1) The source particles of the i-th sub-beam are transported to the k-th step, with the starting point parameters being (r i,k,start ), and the ending point parameters being (r i,k,end ). Among them, i is the sub-beam number where the current particle is located, and k is the number of the current particle transport step. The starting point of the k-th step is offset to other sub-beams to obtain the starting point coordinates (r j,k,start ) and (r j,k,end ) at the k-th step under different sub-beams, where j is the sub-beam number to which the particle is offset; (4.2) Calculate the equivalent path that the particle passes through in the k-th step of particle transport under different sub-beams based on the starting point coordinates of the k-th step of particle transport under different sub-beams and in combination with the image information of the transport medium, and calculate the current sub-beam i where the particle is located and the equivalent path that the particle is deflected to different sub-beams j in the k-th step of particle transport. and Meanwhile, obtain the reaction cross-section of the particle at the k-th step according to the material of the cell where the starting point coordinates are located. and Calculate the transport offset correction parameter when the particle in the i-th sub-beam source is deflected to the j-th sub-beam in the k-th step of particle transport: (4.3) Calculate the collision offset correction parameter when the reaction of the current particle after the k-th step of transport deviates from sub-beam i to sub-beam j according to the end coordinates of the k-th step of particle transport under different sub-beams and combining the image information of the transport medium. Where m is the nuclide number of the collision that occurs after the k-th step of the actual particle transport. is the nucleon density of the nuclide where the collision occurs after the k-th step of the actual particle transport. is the nucleon density of the m-th nuclide at the position corresponding to the end coordinates after the k-th step of the particle transport is offset to the j-th sub-beam. (4.4) Calculate the offset correction parameter for the current step based on the transport offset correction parameter and the collision offset correction parameter obtained according to steps (4.2) and (4.3). (5) Calculate the scattered ray contribution of the current simulated source particle and its offsets to different sub-beams to the two-dimensional plane measurement device: When the source particle transport simulation reaches the point (x, y) of the two-dimensional measurement plane device, count the scattered ray contribution Flux of this particle to the two-dimensional plane measurement device i (x, y), and obtain the total offset correction parameter of the particle offset to different sub-beams according to the offset parameters calculated in steps (3) and (4) where is the offset correction parameter for the particle to offset from sub-beam i to sub-beam j at the k-th step of particle transport. Finally, count the scattered ray contribution Flux j (x, y) = Flux i (x, y) * Bis_Tot i,j ; (6) Calculating the scattered ray contribution of a unit-intensity radiation source: After simulating all source particles, statistically analyze and normalize the scattered ray contributions of each sub-beam to the two-dimensional plane device, and then multiply by the number of particles emitted by the unit-intensity radiation source to calculate the scattered ray contribution values of all sub-beams under unit intensity.

2. The method according to claim 1, wherein The specific content of step (3) is as follows: (3.1) Calculate the generation probability \(P_{i0}(r, \Omega)\) of the source particle in the corresponding sub-beam based on the source particle parameters \((r, \Omega)\) obtained by sampling and the sub-beam to which the source particle belongs, where \(i0\) is the sub-beam number to which the source particle belongs; i0 , \(\Omega\) i0 ) and the sub-beam to which the source particle belongs, calculate the generation probability \(P_{i0}(r, \Omega)\) of the source particle in the corresponding sub-beam, where \(i0\) is the sub-beam number to which the source particle belongs; i0 (r i0 , \(\Omega\) i0 ), where \(i0\) is the sub-beam number to which the source particle belongs; (3.2) Offset the source particle space into other sub-beams to obtain new source particle parameters (r j0 , Ω j0 ), and calculate the generation probability P j0 (r j0 , Ω j0 ) of the source particle in the current sub-beam, where j0 is the number of the sub-beam to which the source particle is offset; (3.3) Based on the above two sets of probabilities, the offset correction parameter Bis_S for the particles of the i0-th sub-beam source offset to the j0-th sub-beam is calculated i0,j0 = P j0 (r j0 , Ω j0 ) / P i0 (r i0 , Ω i0 ).

3. The method according to claim 1, wherein The source distribution information described in step (1) includes field size, source-skin distance, source-axis distance, particle type, and accelerator energy.

4. The method according to claim 1, characterized in that, The medium image information described in step (1) includes CT image resolution, size, and CT value.

5. The method according to claim 1, wherein The two-dimensional plane measurement device information described in step (1) includes the distance from the source to the two-dimensional plane measurement device.

6. A scatter ray calculation system based on Monte Carlo perturbation calculation, characterized in that, Implementing the system based on the scattered ray calculation method based on Monte Carlo perturbation calculation as described in claim 1 includes the following modules: Source sampling module: This module samples source particle parameters according to source distribution information and calculates the offset correction parameters for source particles; Transport simulation module: This module realizes the transport process of particles in the medium and statistically analyzes the scattered ray contributions of the current particle and the offset particle to the counting area; Data processing module: This module reads in source distribution information, medium image information, and two-dimensional plane measurement device information and converts them into a Monte Carlo radiation transport calculation model, and at the same time reads in the reaction cross-section data involved in the particle transport process; Storage medium of the system: This module stores the Monte Carlo radiation transport calculation model containing source distribution information, medium image information, and two-dimensional plane measurement device information, the reaction cross-section data involved in the particle transport process, as well as the transport track information, offset parameters, and counting information during the particle transport process.

7. A computer storage medium, on which a computer program is stored, characterized in that, When the computer program is executed by a processor, it realizes the scattered ray calculation method based on Monte Carlo perturbation calculation as described in any one of claims 1-5.

8. A computer device, comprising a memory, a processor, and a computer program stored on the memory and executable on the processor, characterized in that, When the processor executes the computer program, it realizes the scattered ray calculation method based on Monte Carlo perturbation calculation as described in any one of claims 1-5.

Images