Monte carlo simulation method and system for accelerator target volume dose distribution

Abstract

This invention discloses a Monte Carlo simulation method and system for dose distribution in an accelerator target area, belonging to the field of Monte Carlo algorithm simulation technology. The method includes: Step 1, before the Monte Carlo simulation begins, spurious particles are emitted from the target area in reverse to perform accompanying Monte Carlo transport, generating a dose sensitivity field across the entire space; Step 2, an adaptive response surface mesh is generated based on the dose sensitivity field, this mesh containing a high-resolution voxel mesh in the high-sensitivity region and macroscopic response surface elements in the low-sensitivity region; Step 3, a small-sample preliminary simulation is performed on the adaptive response surface mesh, and a response surface model is fitted within the macroscopic response surface elements. This invention can construct a global dose sensitivity field through reverse accompanying transport in the target area, and construct an adaptive response surface mesh with high and low precision based on sensitivity differences. For the low-sensitivity region, the response surface model is used to analyze the dose instead of traditional full-particle transport calculations, significantly reducing the simulation computation time.

Description

Technical Field

[0001] This invention relates to the field of Monte Carlo algorithm simulation technology, and more specifically, to a Monte Carlo simulation method and system for dose distribution in accelerator target areas. Background Technology

[0002] Monte Carlo simulation of accelerator target dose distribution is a core technology for accurate radiotherapy dose calculation. With its high-precision particle transport simulation capability, it can accurately reproduce the dose deposition pattern of radiotherapy field and human tissue. It is a key technology to ensure the accuracy of radiotherapy treatment, improve the effect of tumor treatment, and reduce radiation damage to normal tissues. At present, this technology has been widely used in various radiotherapy dose verification and calculation scenarios, especially suitable for dose solving of complex fields, irregular target areas and non-uniform human tissues. However, in cutting-edge clinical scenarios such as online adaptive radiotherapy, the existing Monte Carlo simulation technology is difficult to meet the rigid needs of clinical diagnosis and treatment.

[0003] The most prominent technical shortcoming of current high-precision Monte Carlo simulation is the excessively long computation time per simulation. Taking a standard 20×20cm² radiation field simulation as an example, a single high-precision calculation can take up to one hour, which cannot meet the real-time calculation requirements of online adaptive radiotherapy at the minute or even second level, severely limiting the clinical application of this high-precision simulation technology. Delving into its underlying technical roots, Monte Carlo simulation has inherent physical constraints. Its statistical error follows an inverse square relationship with the total number of simulated particles. That is, improving statistical accuracy requires a quadratic increase in the number of simulated particles. When the simulation statistical uncertainty is optimized from 3% to the clinically required high-precision standard of 1%, the required number of particles and the overall computational load will increase by about 9 times. Even a small optimization of accuracy will lead to a sharp deterioration in computation time.

[0004] Current mainstream simulation acceleration methods, including variance reduction algorithms and multi-threaded parallel computing, can only achieve linear or sublinear performance improvements and cannot overcome the inherent inverse square constraint between statistical error and computational load in Monte Carlo simulations. In scenarios involving online adaptive radiotherapy dose recalculation, patient anatomy and treatment beam parameters are prone to dynamic changes during surgery, requiring repeated dose recalculations in clinical practice. However, existing technologies require performing massive particle simulations with uniform precision across the entire space for each parameter change. To meet the requirement of <1% high statistical accuracy, the time consumed in a single calculation far exceeds the clinical real-time threshold, completely negating the timeliness of online adaptive radiotherapy and hindering the widespread application of high-precision Monte Carlo simulation technology in precision radiotherapy clinical scenarios. Summary of the Invention

[0005] To address the problems existing in the prior art, the present invention aims to provide a Monte Carlo simulation method and system for accelerator target area dose distribution. It can construct a global dose sensitivity field through reverse adjoint transport in the target area, and construct an adaptive response surface grid with high and low precision based on sensitivity differences. For low-sensitivity regions, the response surface model is used to analyze the dose instead of traditional whole-particle transport calculation, which greatly reduces the simulation calculation time.

[0006] To solve the above problems, the present invention adopts the following technical solution:

[0007] The first aspect is the Monte Carlo simulation method for accelerator target dose distribution, including:

[0008] Step 1: Before the Monte Carlo simulation begins, pseudo-particles are emitted from the target region in the reverse direction to perform accompanying Monte Carlo transport, generating a dose sensitivity field across the entire space.

[0009] Step 2: Generate an adaptive response surface mesh based on the dose sensitivity field. This mesh contains a high-resolution voxel mesh in the high-sensitivity region and macroscopic response surface elements in the low-sensitivity region.

[0010] Step 3: Perform small-sample preliminary simulation on the adaptive response surface grid, fit the response surface model within the macroscopic response surface cells, and obtain the initial flux distribution in the high-sensitivity region;

[0011] Step 4: Based on the response surface model and the initial flux distribution, perform the main simulation with dynamic weight allocation. The weight of each particle is determined in real time by the dynamic weight function. The dynamic weight function takes the response surface model and the local statistical uncertainty of the currently simulated particles as input, terminates the particles in the low sensitivity region, and analyzes the contribution dose through the response surface model.

[0012] Step 5: During the main simulation, monitor the statistical uncertainty within the high-sensitivity voxel. When the statistical error in a local region exceeds a predetermined threshold, emit local source particles from that region and repeat steps 3 to 4 for sub-loop resimulation for that local region.

[0013] Step 6: Output the expected dose value and its confidence interval for each voxel.

[0014] Further, step 1 includes:

[0015] Step 11: Convert the three-dimensional target area voxel matrix into a reverse-emission accompanying pseudo-particle source term. The number of pseudo-particles emitted is allocated according to the prescription dose weight of the corresponding voxel. Randomly generate an initial position, initial direction, initial energy and initial weight for each pseudo-particle. Perform accompanying Monte Carlo transport from the target area to the beam exit plane. During the transport process, record the three-dimensional spatial track point sequence of each pseudo-particle and the remaining energy and direction weight corresponding to each track point until the pseudo-particle reaches the beam exit plane or the remaining energy is lower than the cutoff threshold, and generate a reverse track database.

[0016] Step 12: Construct a three-dimensional empty counting matrix with the same resolution as the dose calculation grid. Traverse each track point in the reverse track database, determine the grid voxel to which the track point belongs based on its spatial coordinates, and accumulate the current weight corresponding to the track point into the counting matrix element of the voxel. After accumulating all track points, normalize the accumulated count value of each voxel based on the global maximum value of the counting matrix, and output the normalized three-dimensional matrix as the full-space dose sensitivity field.

[0017] Further, step 2 includes:

[0018] Step 21: Compare the standardized sensitivity coefficient of each voxel in the full-space dose sensitivity field with the preset high sensitivity threshold, and mark voxels with sensitivity coefficients greater than or equal to the threshold as high sensitivity voxels, and mark the rest of the voxels as low sensitivity voxels.

[0019] Step 22: Perform connected component labeling on all high-sensitivity voxels to form a spatially continuous set of high-sensitivity regions, each of which contains its voxel list and spatial bounding box.

[0020] Furthermore, step 2 also includes:

[0021] Step 23: Based on the set of high-sensitivity regions, generate a high-resolution voxel mesh for the voxels within each high-sensitivity region according to a preset high-resolution mesh step size.

[0022] Step 24: Based on the labeling results of low-sensitivity voxels, spatially adjacent low-sensitivity voxels are merged into macroscopic response surface units. Each macroscopic response surface unit stores the average sensitivity coefficient of its geometric boundary and internal voxels, and outputs an adaptive response surface mesh composed of high-resolution voxel mesh and macroscopic response surface units.

[0023] Further, step 3 includes:

[0024] Step 31: A small sample set of particles, accounting for a predetermined proportion of the total planned number of particles, is emitted on the adaptive response surface grid. The dose contribution value and its spatial coordinates at multiple sampling points within each macroscopic response surface unit are recorded, and the particle flux count within each voxel in the high-resolution voxel grid is recorded.

[0025] Step 32: Based on the sampling point data inside each macroscopic response surface unit, perform polynomial least squares fitting on the unit to generate response surface model coefficients, and calculate the initial flux value of each voxel based on the particle flux count in each high-resolution voxel, and output the response surface model of each macroscopic response surface unit and the initial flux distribution of the high-sensitivity region.

[0026] Further, step 4 includes:

[0027] Step 41: Based on the polynomial coefficients and geometric boundaries of each macroscopic response surface unit, and the initial flux distribution of each high-sensitivity voxel, construct the initial state dataset of the main simulation. This dataset contains the polynomial coefficients and geometric boundaries of each macroscopic response surface unit, as well as the initial statistical uncertainty estimate and the zeroed cumulative flux counter of each high-sensitivity voxel.

[0028] Step 42: Emit the main simulated particle from the beam source plane and determine the type of the grid element where the particle is currently located in real time;

[0029] Step 43: When the particle is located within a high-sensitivity voxel, update the flux counter and statistical uncertainty of that voxel, and use the updated local statistical uncertainty as feedback input.

[0030] Furthermore, step 4 also includes:

[0031] Step 44: When the particle is located within a macroscopic response surface unit, determine whether the particle should terminate based on the sensitivity coefficient of that unit and the current local statistical uncertainty using the Russian Roulette algorithm.

[0032] Step 45: If the particle is determined to have terminated, the polynomial response surface model of the macroscopic response surface unit where the particle is located is called according to the particle's current spatial coordinates to calculate the dose contribution value to each voxel in the target area and accumulate it to the dose accumulator of the corresponding voxel, and then the particle is terminated; if the particle is determined not to have terminated, the particle weight is scaled according to the survival probability of Russian roulette and then transport continues.

[0033] Further, step 5 includes:

[0034] Step 51: During the main simulation, at predetermined time intervals or after completing a predetermined number of particles, all high-sensitivity voxels are traversed, the current statistical uncertainty of each voxel is read, and it is compared with a preset local statistical error threshold. High-sensitivity voxels whose statistical uncertainty exceeds the threshold are marked.

[0035] Step 52: Perform spatial connectivity region labeling on the marked high-sensitivity voxels to form one or more local regions to be corrected;

[0036] Step 53: For each local region to be corrected, determine the number of additional local source particles to be emitted based on the difference between the average statistical error of the region and the threshold. Project the initial state of the local source particles from the region back onto the beam source plane. Repeat steps 31 to 43 for the local region and its spatial neighborhood. Add the dose contribution obtained in the sub-cycle to the global dose accumulator and update the statistical uncertainty of the high-sensitivity voxels in the region.

[0037] Further, step 6 includes:

[0038] Step 61: Read the total cumulative dose and effective particle count of each voxel in the global dose accumulator, calculate the expected dose value of each voxel, and calculate the upper and lower boundaries of the confidence interval based on the statistical uncertainty of each voxel and the preset confidence level, and output the expected dose value and its confidence interval of each voxel.

[0039] In a second aspect, the present invention also provides a Monte Carlo simulation system for dose distribution in an accelerator target area, comprising: a reverse accompaniment module for emitting pseudo-particles in the reverse direction from the target area before the start of the Monte Carlo simulation to perform accompaniment Monte Carlo transport and generate a dose sensitivity field in the entire space.

[0040] The adaptive mesh module generates an adaptive response surface mesh based on the dose sensitivity field. This mesh contains a high-resolution voxel mesh in the high-sensitivity region and macroscopic response surface elements in the low-sensitivity region.

[0041] The sample fitting module is used to perform small-sample advance simulations on the adaptive response surface grid, fit the response surface model within the macroscopic response surface cells, and obtain the initial flux distribution in the high-sensitivity region.

[0042] The dynamic weighting module performs the main simulation of dynamic weight allocation based on the response surface model and the initial flux distribution. The weight of each particle is determined in real time by the dynamic weighting function, which takes the response surface model and the local statistical uncertainty of the currently simulated particles as inputs. It terminates particles in the low-sensitivity region and analyzes the contribution dose through the response surface model.

[0043] The local correction module is used to monitor the statistical uncertainty within the high-sensitivity voxel during the main simulation. When the statistical error in a local region exceeds a predetermined threshold, local source particles are emitted from that region, and the sample fitting module to the dynamic weighting module is repeatedly executed for the local region to perform a sub-loop resimulation.

[0044] The dose output module is used to output the expected dose value and its confidence interval for each voxel.

[0045] Compared with the prior art, the beneficial effects of the present invention are as follows:

[0046] (1) This scheme constructs a global dose sensitivity field by reverse adjoint transport in the target area, and constructs an adaptive response surface grid with high and low precision based on the sensitivity difference. For low sensitivity regions, the response surface model is used to analyze the dose instead of the traditional whole particle transport calculation, thus getting rid of the computational redundancy problem of uniform sampling of the entire domain in traditional Monte Carlo simulation and greatly reducing the simulation calculation time.

[0047] (2) This scheme adopts a piecewise dynamic weight function that integrates regional sensitivity coefficient and local real-time statistical uncertainty, and combines Russian roulette algorithm to realize adaptive dynamic control of particle survival and weight. While simplifying the invalid particle transport in low-sensitivity areas, it continuously retains the complete particle flux statistical iteration process in high-sensitivity target areas, avoiding the accuracy loss caused by traditional linear acceleration algorithms, and achieving a balance between computational efficiency and core area dose simulation accuracy.

[0048] (3) This scheme adds a periodic error inspection and local sub-cycle re-simulation mechanism for the entire main simulation process, which can accurately locate the local high-sensitivity area where the statistical error exceeds the standard, and adaptively match the number of local supplementary particles according to the degree of error exceeding the standard, and specifically complete the local dose sample supplementation and error correction, thus solving the defects of uneven local convergence and substandard local accuracy in traditional global simulation.

[0049] (4) This scheme can fully quantify the statistical fluctuations and the range of reliable results brought about by random particle sampling, providing complete and quantifiable data support for radiotherapy dose result verification, clinical error tracing, and adaptive radiotherapy dynamic dose assessment, greatly improving the applicability of high-precision Monte Carlo simulation technology in clinical precision radiotherapy scenarios. Attached Figure Description

[0050] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the accompanying drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are merely some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without any creative effort.

[0051] Figure 1 This is a flowchart of the overall method of the present invention;

[0052] Figure 2 This is a flowchart illustrating the generation of a dose sensitivity field via reverse adjoint transport according to the present invention.

[0053] Figure 3 A flowchart illustrating the generation of the adaptive response surface mesh for this invention;

[0054] Figure 4 This is the main simulation flowchart for the dynamic weight allocation of this invention;

[0055] Figure 5 This is a flowchart of the local error monitoring, sub-loop re-simulation, and final output of the present invention;

[0056] Figure 6 This is a data flow diagram between the various modules of the present invention. Detailed Implementation

[0057] The technical solutions of the present invention will be clearly and completely described below with reference to the accompanying drawings of the embodiments of the present invention. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. All other embodiments obtained by those skilled in the art based on the embodiments of the present invention without creative effort are within the scope of protection of the present invention.

[0058] Example 1:

[0059] Please see Figures 1 to 5 Monte Carlo simulation methods for accelerator target dose distribution include:

[0060] Step 1: Before the Monte Carlo simulation begins, pseudo-particles are emitted from the target region in the reverse direction to perform accompanying Monte Carlo transport, generating a dose sensitivity field across the entire space. The specific operation is as follows:

[0061] Using the radiotherapy target area as the particle emission source, pseudo-particles are emitted in reverse and Monte Carlo transport calculations are performed. The contribution characteristics of particles at each location in the three-dimensional space to the dose deposition of the target area are collected in the whole domain. Finally, a standardized dose sensitivity field covering the entire dose calculation space is constructed. This sensitivity field can accurately quantify the influence of each grid voxel in the three-dimensional space on the final dose result of the target area and distinguish between high dose contribution sensitive areas and low dose contribution flat areas.

[0062] Step 1 specifically includes the following steps:

[0063] Step 11: Convert the three-dimensional target area voxel matrix into a reverse-emission accompanying pseudo-particle source term. The emission quantity of each pseudo-particle is allocated according to the prescription dose weight of the corresponding voxel. Randomly generate an initial position, initial direction, initial energy, and initial weight for each pseudo-particle. Perform accompanying Monte Carlo transport from the target area to the beam exit plane. During the transport process, record the three-dimensional spatial track point sequence of each pseudo-particle and the remaining energy and direction weight corresponding to each track point until the pseudo-particle reaches the beam exit plane or the remaining energy is lower than the cutoff threshold, generating a reverse track database. The specific operations are as follows:

[0064] First, the source term transformation is completed based on the three-dimensional target area voxel matrix preset by the radiotherapy dose calculation system. This voxel matrix is ​​a standard numerical model after discretizing the three-dimensional space of the human target area. Each voxel in the matrix corresponds to a fixed three-dimensional spatial coordinate and clinical target area attributes. Each voxel is also configured with a prescription dose weight that matches the clinical radiotherapy plan. This weight corresponds to the treatment priority and dose importance of different target areas. The system allocates the number of pseudo-particles to each voxel according to the prescription dose weight. The higher the prescription dose weight of the target voxel, the more reverse pseudo-particles are allocated. This ensures that the reverse transport data of high clinically important target areas have a sufficient statistical sample basis and avoids the distortion of sensitivity data in high-value areas.

[0065] After the particle number allocation is completed, each accompanying pseudo-particle to be launched is randomly assigned initial parameters across all dimensions. These parameters include the particle's initial three-dimensional spatial position, initial direction of motion, initial energy value, and initial weight value. The range of all parameter assignments matches the actual physical conditions of the medical radiotherapy accelerator and the spatial range of the human target area. After parameter initialization, the accompanying Monte Carlo reverse transport calculation process is initiated. All pseudo-particles start from their initial positions inside the target area and continuously perform physical transport tracking towards the accelerator beam exit plane, simulating the particle's motion, energy attenuation, and directional shift in the equivalent space of human tissue. Within each particle transport step calculation cycle, the current three-dimensional spatial track point coordinates, the remaining energy value after transport, and the real-time direction weight parameters of the pseudo-particles are collected and recorded in real time, and integrated in temporal and spatial order to form a complete track point sequence for a single particle.

[0066] The reverse transport process of a single pseudo-particle continues iteratively until either of the following termination conditions is met: first, the particle's spatial position moves to the accelerator beam exit plane, and the particle leaves the human dose calculation space; second, the particle's energy continuously decays during transport, and the remaining energy value falls below a preset energy cutoff threshold, causing the particle to lose its dose deposition contribution capability. After the transport of a single particle terminates, all its track data are saved. After all assigned pseudo-particles have completed reverse transport tracking, the track point sequences of all particles, corresponding remaining energy, and directional weight data are summarized and integrated to construct a complete reverse track database. This database completely preserves the transport characteristics and dose contribution basic data of reverse particles across the entire space.

[0067] Step 12: Construct a three-dimensional empty counting matrix with the same resolution as the dose calculation grid. Traverse each track point in the reverse track database, determine the corresponding grid voxel based on the spatial coordinates of the track point, and accumulate the current weight corresponding to the track point into the counting matrix element of that voxel. After accumulating all track points, normalize the accumulated count value of each voxel based on the global maximum value of the counting matrix, and output the normalized three-dimensional matrix as the full-space dose sensitivity field. The specific operations are as follows:

[0068] Based on the generated reverse track database, the numerical calculation logic of spatial weight accumulation and global normalization is used to accurately solve and standardize the dose sensitivity field of the entire space, realizing the quantitative visualization of the dose contribution sensitivity of each region in space. First, a set of three-dimensional empty counting matrices is initialized. The spatial size, voxel resolution, and three-dimensional coordinate mapping relationship of this matrix are completely consistent with the dose calculation grid used in the subsequent forward dose simulation, which can realize accurate one-to-one matching of spatial voxel data. In the matrix initialization state, the count values ​​corresponding to all voxels are uniformly set to zero to ensure that there is no initial data interference in the subsequent weight accumulation calculation.

[0069] After initialization, all track point data stored in the reverse track database are traversed one by one. For each independent track point, its corresponding three-dimensional spatial coordinate parameters are accurately extracted. The position of the three-dimensional counting matrix voxel to which the track point belongs is located through a spatial coordinate matching algorithm. Then, the real-time weight value corresponding to the track point is accumulated into the corresponding voxel's value unit in the counting matrix. This traversal and accumulation logic covers all track points in the database, completing the statistical summary of the reverse particle transport weight in the entire space without omission. After accumulation, the value of each voxel in the three-dimensional counting matrix directly corresponds to the original cumulative contribution of the spatial location to the target area dose deposition, which can intuitively reflect the difference in the basic sensitivity of the dose response at different spatial locations.

[0070] To eliminate spatial sensitivity discrimination errors caused by differences in absolute numerical magnitudes and to achieve standardized and unified comparison of global sensitivity data, a global normalization process is performed on the accumulated three-dimensional counting matrix. First, the accumulated values ​​of all voxels in the counting matrix are traversed, and the global maximum count value is extracted and set as the normalization benchmark value. Using this benchmark value as the denominator and the accumulated count value of each voxel as the numerator, a division operation is performed on each voxel value in the matrix one by one, so that the sensitivity values ​​of all voxels are uniformly mapped to the standard range of zero to one. After the normalization calculation of all voxels is completed, the final standardized three-dimensional numerical matrix is ​​output. This matrix is ​​a full-space dose sensitivity field that can be directly applied to adaptive mesh generation. The value of each voxel in the matrix can accurately characterize the dose contribution sensitivity of the corresponding spatial location.

[0071] In a preferred embodiment of the present invention, step 2 is further included: generating an adaptive response surface mesh based on the dose sensitivity field. This mesh includes a high-resolution voxel mesh in the high-sensitivity region and macroscopic response surface elements in the low-sensitivity region. The specific operation is as follows:

[0072] Based on the global standardized dose sensitivity field output in step 1, an adaptive non-uniform mesh reconstruction of the three-dimensional computational space is completed. According to the difference in the dose contribution sensitivity of each region in the space, the mesh accuracy is differentiated and adapted. The overall mesh construction logic uses the voxel sensitivity value as the sole criterion, dividing the entire dose calculation space into two types of differentiated computational unit structures. High-density, high-precision voxel meshes are constructed for regions with a high degree of influence on the target dose result, while large-size, macroscopic response surface units are constructed for regions with a low degree of influence on the target dose result. The differentiated adaptive response surface mesh can significantly reduce the mesh computation amount in low-sensitivity regions while ensuring the computational accuracy in high-sensitivity dose regions.

[0073] Step 2 specifically includes the following steps:

[0074] Step 21: Compare the normalized sensitivity coefficient of each voxel in the full-space dose sensitivity field with a preset high sensitivity threshold. Voxels with sensitivity coefficients greater than or equal to the threshold are marked as high-sensitivity voxels, and the remaining voxels are marked as low-sensitivity voxels. The specific operation is as follows:

[0075] In the full-space dose sensitivity field, each grid voxel corresponds to a standardized sensitivity coefficient in the range of 0 to 1. The magnitude of this coefficient directly represents the influence of particle transport behavior at the corresponding spatial location on the final dose deposition result of the target area. The larger the value, the higher the dose contribution weight of that spatial location. A pre-configured high sensitivity threshold with fixed parameters is retrieved. This threshold is a fixed judgment value preset based on the accuracy requirements of clinical radiotherapy dose calculation and does not require real-time dynamic adjustment. All three-dimensional voxels in the dose sensitivity field are traversed, and the standardized sensitivity coefficient of each voxel is compared with the preset high sensitivity threshold. When the standardized sensitivity coefficient of a voxel is greater than or equal to the preset threshold, the voxel is assigned a high sensitivity voxel label. When the standardized sensitivity coefficient of a voxel is less than the preset threshold, the voxel is assigned a low sensitivity voxel label. After all voxels in the entire domain have been compared and labeled one by one, the voxel attributes of the entire computational space are accurately divided, providing a clear voxel classification basis for subsequent connected region extraction and differential grid reconstruction.

[0076] Step 22: Perform connected component labeling on all high-sensitivity voxels to form a spatially contiguous set of high-sensitivity regions. Each high-sensitivity region contains its voxel list and spatial bounding box. The specific operations are as follows:

[0077] After completing the global voxel high and low sensitivity labeling, the high-sensitivity voxels in space are mostly distributed in a scattered and irregular form in the target area and surrounding key dose regions, which cannot be directly used for the construction of a regular high-precision mesh. A three-dimensional spatial connected region traversal algorithm is used to determine the neighborhood association of all labeled high-sensitivity voxels. High-sensitivity voxels that are adjacent vertically, horizontally, or front-back in three-dimensional space are grouped into the same connected whole, while clusters of high-sensitivity voxels that are spatially isolated and have no adjacent voxel associations are independently divided into separate connected regions. After the global traversal and aggregation, a set of high-sensitivity regions consisting of multiple independent spatial regions is formed. For each high-sensitivity connected region in the set, the coordinates and sensitivity coefficients of all voxels contained in the region are automatically summarized and recorded to form a voxel list. Simultaneously, the smallest three-dimensional spatial bounding box, i.e., the spatial bounding box, that encloses all voxels in the connected region is fitted and generated to accurately define the spatial range of each high-sensitivity region.

[0078] Step 23: Based on the set of high-sensitivity regions, generate a high-resolution voxel mesh for each voxel within the high-sensitivity region according to a preset high-resolution mesh step size. The specific operation is as follows:

[0079] The pre-set high-resolution grid step size parameter is retrieved. This step size is smaller than the step size of the original base grid of the dose sensitivity field and is a fixed spatial size parameter adapted to the needs of high-precision dose statistics. For each independent high-sensitivity region in the set of high-sensitivity regions, the original base voxel grid inside the region is uniformly subdivided and split using the spatial bounding box corresponding to the region as the boundary range and the pre-set high-resolution grid step size as the uniform subdivision scale. The subdivision process strictly follows the actual spatial contour of the high-sensitivity region, without expanding outward into the low-sensitivity spatial region, and without omitting any spatial position inside the region. Finally, a higher resolution voxel grid with smaller size, higher spatial resolution, and denser number is generated inside each high-sensitivity region. This refined grid structure can accurately match the dose change characteristics of the high-sensitivity region, ensuring the refinement of subsequent particle flux statistics, dose calculation, and error monitoring, and maintaining the bottom line of dose calculation accuracy of the overall simulation.

[0080] Step 24: Based on the labeling results of low-sensitivity voxels, spatially adjacent low-sensitivity voxels are merged into macroscopic response surface elements. Each macroscopic response surface element stores the average sensitivity coefficient of its geometric boundary and internal voxels. The output is an adaptive response surface mesh composed of a high-resolution voxel mesh and macroscopic response surface elements. The specific operations are as follows:

[0081] Based on the global low-sensitivity voxel distribution results marked in step 21, spatial correlation aggregation is performed on all original mesh voxels with low-sensitivity attributes. Adjacent, continuously distributed low-sensitivity voxels in 3D space are integrated into an independent macroscopic response surface unit. Each merged macroscopic response surface unit possesses an independent and complete 3D geometric boundary. The boundary space coordinate parameters of each unit are recorded, and the standardized sensitivity coefficients of all original low-sensitivity voxels within the unit are collected. The average sensitivity coefficient of the macroscopic unit is calculated using an arithmetic mean, and the geometric boundary parameters and average sensitivity coefficient are permanently stored as inherent properties of the unit. After merging and storing the parameters of all low-sensitivity region units, the system globally integrates the high-resolution voxel mesh generated from the subdivision of the high-sensitivity region with the macroscopic response surface unit generated from the merging of the low-sensitivity region. The two mesh structures seamlessly cover the entire 3D dose calculation space, ultimately forming a complete adaptive response surface mesh, completing the output of this step.

[0082] In a preferred embodiment of the present invention, step 3 is further included: performing a small-sample preliminary simulation on the adaptive response surface grid, fitting the response surface model within the macroscopic response surface cells, and obtaining the initial flux distribution in the high-sensitivity region. The specific operations are as follows:

[0083] For macroscopic response surface units with low computational accuracy requirements and gradual spatial dose changes, this step completes the fitting of a polynomial response surface model using a small amount of sampled data. This numerical model replaces the method of solving for dose using massive particle transport data, enabling rapid analytical calculation of dose in low-sensitivity regions. For high-resolution voxel grid regions with high computational accuracy requirements and significant dose gradient changes, this step obtains the initial flux distribution of each fine voxel by statistically analyzing the transport flux of a small sample of particles. This provides initial baseline data for the dynamic adjustment of particle weights and the calculation of statistical uncertainty in the main simulation stage. This small-sample preliminary simulation initializes the parameters of the entire computational space with extremely low computational cost. This avoids the computational waste of traditional Monte Carlo simulations that require calculating all particles from scratch, and provides a complete model foundation and initial data support for the subsequent high-precision and high-efficiency dynamic weight main simulation.

[0084] Step 3 specifically includes the following steps:

[0085] Step 31: A small sample particle set, representing a predetermined proportion of the total planned particle count, is emitted onto the adaptive response surface grid. The dose contribution value and its spatial coordinates at multiple sampling points within each macroscopic response surface unit are recorded. The particle flux count within each voxel in the high-resolution voxel grid is also recorded. The specific operations are as follows:

[0086] Before the formal simulation calculation process begins, the total planned number of particles required for the overall dose simulation is pre-set. This value represents the full number of particles required to meet the accuracy requirements of clinical dose statistics and is the core particle configuration parameter for the main simulation phase. A small sample particle set is formed by extracting a portion of the total planned particles according to a pre-set fixed proportion. This pre-set proportion is a fixed simulation configuration parameter that can be fixed in advance according to the needs of accuracy and efficiency adaptation, without requiring real-time dynamic adjustment. The small sample particle set is emitted normally from the accelerator beam source plane, following standard medical particle transport physics rules. It completes physical processes such as transport, scattering, and energy attenuation within a complete adaptive response surface grid space. The particle transport range simultaneously covers the high-resolution voxel grid in the high-sensitivity region and the macroscopic response surface unit in the low-sensitivity region.

[0087] During the entire transport process of small-sample particles, a partitioned synchronous data recording mechanism is implemented. For all macroscopic response surface units, multiple spatial sampling points are evenly distributed within the geometric space of each unit. The dose contribution value generated when particles are transported through the sampling point is captured in real time. At the same time, the three-dimensional spatial coordinates corresponding to each effective sampling point are accurately recorded. Multiple sets of spatial position and dose contribution data within a single macroscopic unit are continuously accumulated to ensure that the number of samples in each macroscopic unit meets the basic requirements for subsequent polynomial fitting. For high-resolution voxel grids in high-sensitivity regions, the number of particles passing through each voxel is accumulated in real time as an independent statistical unit to complete the particle flux count statistics for each fine voxel. After all small-sample particles complete the complete transport process and trigger the termination condition, the multi-point sampling dataset of all macroscopic response surface units and the flux count dataset of all high-resolution voxels are retained to complete the preliminary data acquisition work.

[0088] Step 32: Based on the sampling point data within each macroscopic response surface unit, perform polynomial least squares fitting on the unit to generate response surface model coefficients, and calculate the initial flux value of each voxel based on the particle flux count within each high-resolution voxel. Output the response surface model of each macroscopic response surface unit and the initial flux distribution in the high-sensitivity region. The specific operations are as follows:

[0089] For each independent macroscopic response surface unit, the three-dimensional spatial coordinates and corresponding dose contribution values ​​of all sampling points within the unit are retrieved. A polynomial least squares fitting algorithm is used to construct a continuous mapping relationship between spatial coordinates and dose contribution values. The logic of this fitting method is to adjust the coefficients of each polynomial so that the sum of squared residuals between the model-fitted dose calculation value and the actual sampled dose value of all sampling points reaches the minimum value. This ensures that the fitted response surface model can fit the real dose distribution law of the low-sensitivity region to the greatest extent and avoid the model distortion problem caused by single-point sampling error.

[0090] This section introduces the calculation expression for least squares fitting, which is derived from the mathematical logic of optimal error matching. The optimal model coefficients are solved using the global residual minimization criterion, serving as the criterion for judging polynomial response surface fitting. The specific formula is as follows:

[0091] This formula constructs a global error function by accumulating the squared errors of all sampling points. It then solves for the optimal polynomial coefficients with the goal of minimizing the global error, thus avoiding the problem of positive and negative errors canceling each other out and ensuring the overall fitting accuracy of the model. In the formula, S represents the sum of squared residuals from all sampling points, and n represents the total number of sampling points within a single macroscopic response surface unit. This represents the actual dose contribution value collected at the i-th sampling point. This represents the fitted dose value at the i-th sampling point calculated using a polynomial model. By iteratively solving this optimal solution problem, all polynomial coefficients suitable for the macroscopic unit are finally determined, forming an independent unit response surface model that can be directly used for rapid analytical solution of doses in subsequent low-sensitivity regions.

[0092] While fitting the macroscopic unit model, the system processes flux data in the high-sensitivity region, retrieves the cumulative particle flux counts of each high-resolution voxel, and combines this with particle transport conditions simulated in a small sample. It then performs statistical conversion on the cumulative flux data of individual voxels to obtain the initial flux value for each high-resolution voxel. The initial flux values ​​of all high-sensitivity voxels are integrated across the entire region to form a continuous and complete initial flux distribution result for the high-sensitivity region, accurately characterizing the initial particle transport state of the high-precision grid region. After completing all calculations, the system outputs a unified set of polynomial response surface model coefficients for all macroscopic response surface units and the initial flux distribution data for the high-sensitivity region.

[0093] In a preferred embodiment of the present invention, step 4 is further included: based on the response surface model and the initial flux distribution, a master simulation with dynamic weight allocation is performed, wherein the weight of each particle is determined in real time by a dynamic weight function. The dynamic weight function takes the response surface model and the local statistical uncertainty of the currently simulated particles as input, terminates the particles in the low-sensitivity region, and analyzes the contribution dose through the response surface model. The specific operation is as follows:

[0094] Based on the polynomial response surface models of each macroscopic response surface unit output in step 3 and the initial flux distribution data of the high-sensitivity region, a piecewise dynamic weight function is used to complete the real-time adaptive adjustment of the weights for all particles, thereby conducting differentiated main simulation particle transport calculations. The dynamic weight function can output the effective weight value of particles in real time according to the changes in the statistical convergence state of local regions during particle transport. During the simulation, the two types of particle handling logic are divided based on the grid attribute discrimination results. The complete particle transport process is retained in the high-sensitivity region, and the flux and statistical error are updated synchronously. In the low-sensitivity macroscopic unit region, the particle selection and rejection is completed by using the dynamic weight value in conjunction with the Russian roulette algorithm. Particles that are determined to be terminated are no longer transported, and instead, the pre-built response surface model is called to quickly calculate the dose contribution through analytical calculation. By relying on the combination of dynamic weight and response surface analytical calculation, the iteration overhead of invalid particles in the low-sensitivity region is compressed, and the overall simulation time is shortened while maintaining the dose statistical accuracy of the core target area. The specific expression of the dynamic weight function is as follows:

[0095] This formula, starting from the physical requirement of variance reduction, introduces only local uncertainty parameters for error-driven weight fine-tuning in the high-sensitivity partition, while the macroscopic partition additionally embeds response surface coefficients to correlate the intrinsic dose response level of the unit. Both types of partitions share an initial weight benchmark to ensure dose expectation conservation. In the formula, The real-time computational weights of the particles after dynamic correction represent the weights. This represents the initial intrinsic weight assigned to a particle when it is emitted from the beam source. Represents the instantaneous three-dimensional spatial coordinates of the particle. Representing coordinates The local statistical uncertainty is updated in real time at the corresponding location, where k is a weight correction constant specific to the high-sensitivity grid. This represents the arithmetic mean of the polynomial coefficients of all response surfaces of the m-th macroscopic response surface unit. This is the correction constant corresponding to the response surface coefficients. k is the correction constant corresponding to the uncertainty within the macroscopic unit. , Parameters were assigned in advance during the simulation initialization phase.

[0096] Step 4 specifically includes the following steps:

[0097] Step 41: Based on the polynomial coefficients and geometric boundaries of each macroscopic response surface unit, and the initial flux distribution of each high-sensitivity voxel, construct the initial state dataset for the main simulation. This dataset contains the polynomial coefficients and geometric boundaries of each macroscopic response surface unit, as well as the initial statistical uncertainty estimate and the zeroed cumulative flux counter for each high-sensitivity voxel. The specific operations are as follows:

[0098] First, we collect the two types of intrinsic parameters of all macroscopic response surface units in batches. One type is the complete set of polynomial coefficients obtained by least squares fitting in step 32, which is the dynamic weighting function. One type of data source is the original data source, and the other is the geometric boundary coordinate data of each macroscopic unit retained in step 24. The geometric boundary is used to quickly determine whether a particle falls into the corresponding macroscopic unit range during particle transport. Simultaneously, the initial flux distribution data of each high-resolution voxel in the high-sensitivity region, obtained through small-sample simulation calculations, is retrieved. Based on the initial flux values, the initial statistical uncertainty of each fine voxel is estimated. This initial uncertainty is the initial uncertainty of the dynamic weighting function during its initial calculation. The initial values ​​are set. To isolate the residual flux statistics generated in the small sample preliminary simulation stage from interfering with the main simulation count, the cumulative flux counters associated with all high-resolution voxels are reset to zero. The above three types of parameters are integrated into a unified initial state dataset after structured aggregation. The baseline coefficients and initial error parameters required for the subsequent operation of the dynamic weight function are all retrieved from this dataset, ensuring that the source of the weight calculation parameters throughout the entire process is unified and traceable.

[0099] Step 42: Launch the main simulated particle from the beam source plane and determine the type of the mesh element where the particle is currently located in real time. The specific operation is as follows:

[0100] According to the pre-set total number of particles, the main simulated particles are generated and emitted sequentially from the source plane corresponding to the accelerator beam exit. The initial energy, initial direction of motion, and other physical parameters of the particles are matched with the conventional output parameters of the radiotherapy beam of a medical linear accelerator, and a fixed initial weight is uniformly assigned at the moment of emission. At any given moment, as the particle continues to move and transport within the equivalent medium space of the human body, its current three-dimensional coordinates are continuously read. The coordinate values ​​are compared one by one with the high-resolution voxel boundaries and the geometric boundaries of each macroscopic response surface unit stored in the adaptive response surface mesh. This quickly distinguishes whether the particle is currently within the coverage area of ​​the high-resolution voxel or the coverage area of ​​the macroscopic response surface unit. The spatial discrimination result will serve as the trigger condition for the selection of the segmented calculation formula of the dynamic weight function. When it is determined to be a high-sensitivity region, the upper segment formula of the function is activated, and when it is determined to be a macroscopic unit region, the lower segment formula of the function is activated. The automatic switching of the weight calculation rules is achieved by relying on real-time spatial discrimination.

[0101] Step 43: When the particle is located within a high-sensitivity voxel, update the flux counter and statistical uncertainty of that voxel, and use the updated local statistical uncertainty as feedback input. The specific operation is as follows:

[0102] When the spatial determination result confirms that the particle falls within the fine voxel corresponding to the high sensitivity, the particle termination or weight scaling operation is not triggered. The cumulative flux counter corresponding to the target voxel is updated according to the particle's current real-time weight, completing the single-step flux increment statistics. After the flux value is refreshed, the new local statistical uncertainty is iteratively solved based on the effective flux sample size accumulated for that voxel. The updated uncertainty value is immediately stored in the parameter storage location corresponding to the initial state dataset, serving as the input for the dynamic weight function calculation of subsequent particles passing near that voxel. The dynamic weight function directly reads the updated value when calculating the weights of subsequent particles. If local uncertainty increases, the function appropriately increases the computational weight of particles passing through the region by relying on the internal coefficient k, thereby increasing the effective sampling ratio of particles in the high-error region and gradually reducing the local statistical bias.

[0103] Step 44: When the particle is located within a macroscopic response surface cell, the Russian Roulette algorithm is used to determine whether the particle should terminate, based on the sensitivity coefficient of that cell and the current local statistical uncertainty. The specific operation is as follows:

[0104] When a particle is determined to have entered any macroscopic response surface element, the mean value of the polynomial coefficients of that element is first retrieved. and the real-time local statistical uncertainty of the particle's location. The two parameters are substituted into the macroscopic partitioning expression of the dynamic weight function to obtain the real-time weight w of the particle within the macroscopic unit. Using a pre-set fixed reference weight as a baseline value, the survival probability of the particle corresponding to this round of the Russian Roulette algorithm is calculated by combining it with the obtained dynamic weight. The survival probability fluctuates synchronously with the dynamic weight value; the larger the dynamic weight, the lower the corresponding survival probability. Probability comparison is completed using uniform random number sampling, and the final output is either a decision to terminate or continue transport of the particle, achieving dynamic screening of redundant particles in low-sensitivity regions.

[0105] Step 45: If the particle is determined to have terminated, the polynomial response surface model of the macroscopic response surface unit where the particle is located is called to calculate the dose contribution value to each voxel in the target area based on the particle's current spatial coordinates, and the value is accumulated to the dose accumulator of the corresponding voxel. Then, the particle is terminated. If the particle is determined not to have terminated, the particle weight is scaled according to the survival probability of Russian roulette and the transport continues. The specific operation is as follows:

[0106] If the judgment command requires the current particle to be terminated, extract the particle's instantaneous three-dimensional coordinates. The polynomial response surface model of the macroscopic unit to which the coordinates belong is matched. The spatial coordinates are substituted into the fitted polynomial analytical expression to directly solve for the dose contribution value of the particle to each target voxel in the entire space. The calculated doses are then added and stored in the global dose accumulator. After the dose collection is completed, the entire transport process of the particle is terminated directly, and no further computing power is consumed for subsequent particle step tracking. If the judgment command allows the particle to continue transporting, the original dynamic weight of the particle is inversely scaled using the survival probability obtained from the aforementioned conversion. The scaled new weight can ensure that the mathematical expectation dose corresponding to the particle remains unchanged, avoiding the additional statistical error introduced by the probability screening operation. After the weight is corrected, the particle continues to travel in space according to the physical transport law, carrying the new weight. When passing through different grid regions, the dynamic weight calculation and partitioning process is repeated.

[0107] In a preferred embodiment of the present invention, step 5 is further included: during the main simulation, the statistical uncertainty within the high-sensitivity voxel is monitored; when the statistical error in a local region exceeds a predetermined threshold, local source particles are emitted from that region, and steps 3 to 4 are repeated for the local region to perform a sub-loop resimulation. The specific operation is as follows:

[0108] After the routine global master simulation completes the batch particle transport, some high-sensitivity fine voxel regions may experience uneven local statistical convergence speeds due to the random sampling characteristics of particles. The statistical uncertainty in some regions may not fall back to the clinically preset accuracy range. If the calculation results are directly output, there will be local dose calculation deviations. This step accurately captures the regions with substandard accuracy through periodic error inspection. Only for the local spatial regions with excessive errors, a dedicated resimulation sub-loop is started. The particle scale is adaptively configured according to the degree of error exceeding the standard in the local region. The local dose sample is completed through local reverse source particle transport and lightweight sub-loop calculation. Without affecting the global calculation efficiency, the consistent convergence of statistical accuracy in all high-sensitivity key regions is achieved, ensuring the uniformity of global accuracy of the overall dose calculation results.

[0109] Step 5 specifically includes the following steps:

[0110] Step 51: During the main simulation, at predetermined time intervals or after completing a predetermined number of particles, traverse all high-sensitivity voxels, read the current statistical uncertainty of each voxel, compare it with a preset local statistical error threshold, and mark high-sensitivity voxels whose statistical uncertainty exceeds the threshold. The specific operation is as follows:

[0111] This step enables periodic, global screening and anomaly marking of high-sensitivity voxel statistical errors during the main simulation process, serving as a pre-sensing step for local accuracy correction. The system incorporates a dual-dimensional inspection trigger mechanism, which can use either a preset fixed program execution time interval or a preset total number of completed simulation particles as the inspection start condition. These two trigger conditions are independent of each other; any condition met during the main simulation will immediately initiate a global error detection process, achieving dynamic accuracy monitoring throughout the entire simulation process. After each inspection process is initiated, all high-sensitivity voxels contained in the high-resolution voxel grid are traversed. The real-time statistical uncertainty value of each voxel, which is continuously updated during the main simulation, is retrieved one by one. This value is obtained by statistically solving the cumulative flux sample and effective particle count of the voxel, and can intuitively represent the statistical fluctuation level of the current voxel dose calculation result. The real-time statistical uncertainty of each voxel is compared with the pre-fixed local statistical error threshold. This threshold is a fixed critical parameter that limits the dose calculation accuracy of the high-sensitivity area. Any high-sensitivity voxel whose real-time statistical uncertainty value exceeds this critical parameter will be marked by the system as an abnormal voxel to be corrected, thus completing the accuracy screening and anomaly collection of high-sensitivity voxels in the entire domain.

[0112] Step 52: Perform spatial connectivity labeling on the marked high-sensitivity voxels to form one or more local regions to be corrected. The specific operations are as follows:

[0113] After the global screening and labeling in step 51, anomalous voxels are distributed in a discrete and scattered manner within the high-sensitivity region. Directly resimulating a single voxel would generate a large amount of repetitive boundary calculation overhead, reducing correction efficiency. A three-dimensional spatial connected region labeling algorithm, consistent with the aforementioned high-sensitivity region construction logic, is adopted to determine the neighborhood association of all labeled anomalous voxels. Anomalous voxels that are directly adjacent vertically, horizontally, vertically, or vertically in three-dimensional space are grouped into a single continuous whole. Spatially isolated clusters of anomalous voxels without adjacency are independently divided into separate correction units, ultimately forming one or more local regions to be corrected with complete boundaries and clearly defined spatial ranges. Each local region to be corrected retains its own set of voxels and spatial boundary range, ensuring that subsequent local resimulation can be carried out in regular spatial units, avoiding computational redundancy caused by discrete single-point correction.

[0114] Step 53: For each local region to be corrected, determine the number of additional local source particles to be emitted based on the difference between the average statistical error of the region and the threshold. Generate the initial state of local source particles by back-projecting from the region onto the beam source plane. Repeat steps 31 to 43 for the local region and its spatial neighborhood. Accumulate the dose contribution obtained in this sub-cycle into the global dose accumulator and update the statistical uncertainty of the high-sensitivity voxels in the region. The specific operation is as follows:

[0115] For each independent local region to be corrected, the average statistical uncertainty of all anomalous voxels within the region is first calculated. The deviation of this average value from a preset local statistical error threshold is used to adaptively determine the number of local source particles required for this local resimulation. The configuration logic for this particle number conforms to the inherent statistical laws of Monte Carlo simulation; the greater the error exceeding the threshold, the more additional sampling particles are required, thus ensuring that the local error can quickly converge to the threshold range. Here, a formula for calculating the number of local source particles is introduced. This formula relies on Monte Carlo particle transport statistics theory, where the convergence speed of statistical uncertainty is inversely proportional to the square root of the number of simulated particles. To ensure that the exceeding error converges quickly to the acceptable range, the number of additional particles needs to be increased by a factor of the square of the error ratio to match the correction requirements of different degrees of local error. By adaptively matching the particle deployment scale through the error deviation ratio, a dynamic balance between accuracy and computing power is achieved. The specific formula is as follows:

[0116] In the formula, This represents the total number of local source particles that need to be deployed to the local area to be corrected. This represents the preset local correction reference particle number. This represents the average statistical uncertainty of the local region to be corrected. This represents the preset local statistical error threshold.

[0117] After determining the number of local particles, the particle trajectory is back-projected from the spatial range of the local region to be corrected to the accelerator beam source plane based on the reverse transport mapping logic. Combined with the physical parameters of the beam source, complete initial state parameters such as the initial position, initial energy, initial direction, and initial weight of the local source particles are generated, constructing a dedicated local correction particle source. Subsequently, for the local region to be corrected and its adjacent spatial neighborhood, the sub-loop simulation process of steps 31 to 43 is repeated. Within the local area, the entire process of small-sample pre-sampling, response surface model parameter updating, dynamic weighted particle transport, flux counting, and iterative updating of local uncertainty is completed. Only the target local region is subjected to refined re-simulation, eliminating the need for repeated global calculations. All dose contribution data generated during the sub-loop simulation are accumulated in real time into the global dose accumulator, supplementing and improving the dose statistics samples of the local region. After the sub-loop operation is completed, the statistical uncertainty of all high-sensitivity voxels in the region is recalculated, completing the update and iteration of local accuracy parameters. This brings the previously excessive local statistical error back to the preset acceptable range, completing the single-loop local accuracy correction.

[0118] In a preferred embodiment of the present invention, step 6 is further included: reading the cumulative dose total value and effective particle count of each voxel in the global dose accumulator, calculating the expected dose value of each voxel, and calculating the upper and lower boundaries of the confidence interval based on the statistical uncertainty of each voxel and the preset confidence level, and outputting the expected dose value and its confidence interval of each voxel. The specific operation is as follows:

[0119] After the main simulation continues running and all local error correction sub-loops have finished, the particle sampling and model analysis calculation process of the global dose field is completely terminated. The system no longer adds particle transport and dose accumulation operations and enters the steady-state data settlement stage. First, all structured data collected by the global dose accumulator during the entire simulation process is read. Two types of core statistical parameters are extracted for each voxel: the cumulative total dose value and the effective particle count for a single voxel. The cumulative total dose value integrates the total contribution of the direct particle transport deposition dose in the high-sensitivity region and the dose calculated by the response surface model in the low-sensitivity region. The effective particle count is the total number of particle samples that participate in the dose statistics of that voxel and have an effective weighted contribution. These two types of parameters completely record the global dose accumulation statistics process of a single voxel.

[0120] The expected dose value of each grid voxel is calculated based on voxel cumulative dose data and effective particle sample data. This value is the true steady-state dose value after eliminating random sampling fluctuations and is the core benchmark value for evaluating radiotherapy dose distribution. This calculation logic is based on the mean convergence principle of statistical simulation. By eliminating statistical fluctuations caused by random particle transport through the ratio of cumulative total dose to effective statistical sample size, the converged stable dose result is obtained. The corresponding mathematical expression is: This formula is derived from the principle of sample mean in mathematical statistics. The dose of a single voxel in Monte Carlo simulation is formed by the superposition of random contributions from a large number of independent particles. By dividing the total cumulative dose by the effective particle count, random errors can be offset to obtain the statistically significant expected dose value. In the formula, The expected dose value representing a single voxel. The total dose contribution value representing voxel accumulation. The effective particle count represents the number of particles corresponding to a voxel.

[0121] After solving for the expected dose, the final local statistical uncertainty obtained from the full iterative convergence of each voxel, along with a preset fixed confidence level parameter, is used to calculate the confidence interval of the dose result for a single voxel. This is used to quantify the statistical fluctuation range and data reliability of the expected dose. The confidence level is a pre-defined simulation configuration parameter used to define the reliable coverage of the statistical results; fixed standard values ​​are typically used in clinical radiotherapy dose calculation scenarios. Based on large-sample statistical theory, the voxel dose statistical results simulated by this method approximately follow a normal distribution. Based on this, the upper and lower boundaries of the two-sided confidence interval can be solved, and the corresponding mathematical expressions are: This formula is derived from the principle of normal distribution interval estimation. Centered on the converged dose expectation value, it scales the relative statistical error using a confidence coefficient to ultimately define the reliable fluctuation range of the dose result, thus adapting to the error quantification requirements of Monte Carlo simulation's stochastic statistical characteristics. In the formula, This represents the upper boundary value of the confidence interval. denoted by , z represents the lower boundary value of the confidence interval, z represents the critical coefficient of the normal distribution corresponding to the preset confidence level, and U represents the local statistical uncertainty of the voxel's final convergence.

[0122] By traversing all voxels within the 3D dose calculation grid, the expected dose value and the upper and lower boundaries of the confidence interval are calculated and collected for each voxel. The dose value and corresponding statistical error range of each spatial location are fully preserved. Finally, the expected dose value and matching confidence interval results of all voxels in the global domain are uniformly output. This output mode differs from the traditional Monte Carlo simulation, which only outputs a single dose value. It simultaneously carries dose amplitude and statistical convergence accuracy information, and fully presents the calculation accuracy distribution of the global dose field. This not only meets the numerical output requirements of radiotherapy dose calculation, but also provides complete data support for clinical reliability assessment, error tracing, and accuracy verification of dose results.

[0123] Example 2:

[0124] Please see Figure 6 Based on Example 1, this embodiment provides a Monte Carlo simulation system for accelerator target dose distribution, including:

[0125] The reverse accompaniment module is used to emit pseudo-particles in the reverse direction from the target area before the start of the Monte Carlo simulation to perform accompaniment Monte Carlo transport and generate a dose sensitivity field in the whole space.

[0126] The adaptive mesh module generates an adaptive response surface mesh based on the dose sensitivity field. This mesh contains a high-resolution voxel mesh in the high-sensitivity region and macroscopic response surface elements in the low-sensitivity region.

[0127] The sample fitting module is used to perform small-sample advance simulations on the adaptive response surface grid, fit the response surface model within the macroscopic response surface cells, and obtain the initial flux distribution in the high-sensitivity region.

[0128] The dynamic weighting module performs the main simulation of dynamic weight allocation based on the response surface model and the initial flux distribution. The weight of each particle is determined in real time by the dynamic weighting function, which takes the response surface model and the local statistical uncertainty of the currently simulated particles as inputs. It terminates particles in the low-sensitivity region and analyzes the contribution dose through the response surface model.

[0129] The local correction module is used to monitor the statistical uncertainty within the high-sensitivity voxel during the main simulation. When the statistical error in a local region exceeds a predetermined threshold, local source particles are emitted from that region, and the sample fitting module to the dynamic weighting module is repeatedly executed for the local region to perform a sub-loop resimulation.

[0130] The dose output module is used to output the expected dose value and its confidence interval for each voxel.

[0131] The above description is merely a preferred embodiment of the present invention; however, the scope of protection of the present invention is not limited thereto. Any equivalent substitutions or modifications made by those skilled in the art within the scope of the technology disclosed in the present invention, based on the technical solution and its improved concept, should be covered within the scope of protection of the present invention.

Claims

1. A Monte Carlo simulation method for dose distribution in an accelerator target area, characterized in that, include: Step 1: Before the Monte Carlo simulation begins, pseudo-particles are emitted from the target region in the reverse direction to perform accompanying Monte Carlo transport, generating a dose sensitivity field across the entire space. Step 2: Generate an adaptive response surface mesh based on the dose sensitivity field. This mesh contains a high-resolution voxel mesh in the high-sensitivity region and macroscopic response surface elements in the low-sensitivity region. Step 3: Perform small-sample preliminary simulation on the adaptive response surface grid, fit the response surface model within the macroscopic response surface cells, and obtain the initial flux distribution in the high-sensitivity region; Step 4: Based on the response surface model and the initial flux distribution, perform a master simulation with dynamic weight allocation. The weight of each particle is determined in real-time by a dynamic weighting function, which takes the response surface model and the local statistical uncertainty of the currently simulated particles as input. Particles are terminated in low-sensitivity regions, and the contribution dose is analyzed using the response surface model. Step 4 includes: Step 41: Based on the polynomial coefficients and geometric boundaries of each macroscopic response surface unit, and the initial flux distribution of each high-sensitivity voxel, construct the initial state dataset of the main simulation. This dataset contains the polynomial coefficients and geometric boundaries of each macroscopic response surface unit, as well as the initial statistical uncertainty estimate and the zeroed cumulative flux counter of each high-sensitivity voxel. Step 42: Emit the main simulated particle from the beam source plane and determine the type of the grid element where the particle is currently located in real time; Step 43: When the particle is located within a high-sensitivity voxel, update the flux counter and statistical uncertainty of that voxel, and use the updated local statistical uncertainty as feedback input. Step 5: During the main simulation, monitor the statistical uncertainty within the high-sensitivity voxel. When the statistical error in a local region exceeds a predetermined threshold, emit local source particles from that region and repeat steps 3 to 4 for a sub-loop resimulation of that local region; Step 5 includes: Step 51: During the main simulation, at predetermined time intervals or after completing a predetermined number of particles, all high-sensitivity voxels are traversed, the current statistical uncertainty of each voxel is read, and it is compared with a preset local statistical error threshold. High-sensitivity voxels whose statistical uncertainty exceeds the threshold are marked. Step 52: Perform spatial connectivity region labeling on the marked high-sensitivity voxels to form one or more local regions to be corrected; Step 53: For each local region to be corrected, determine the number of additional local source particles to be emitted based on the difference between the average statistical error of the region and the threshold. Project the initial state of the local source particles from the region back onto the beam source plane. Repeat steps 3 to 4 for the local region and its spatial neighborhood. Add the dose contribution obtained from this repeated execution of steps 3 to 4 to the global dose accumulator and update the statistical uncertainty of the high-sensitivity voxels in the region. Step 6: Output the expected dose value and its confidence interval for each voxel.

2. The Monte Carlo simulation method for accelerator target dose distribution according to claim 1, characterized in that, Step 1 includes: Step 11: Convert the three-dimensional target area voxel matrix into a reverse-emission accompanying pseudo-particle source term. The number of pseudo-particles emitted is allocated according to the prescription dose weight of the corresponding voxel. Randomly generate an initial position, initial direction, initial energy and initial weight for each pseudo-particle. Perform accompanying Monte Carlo transport from the target area to the beam exit plane. During the transport process, record the three-dimensional spatial track point sequence of each pseudo-particle and the remaining energy and direction weight corresponding to each track point until the pseudo-particle reaches the beam exit plane or the remaining energy is lower than the cutoff threshold, and generate a reverse track database. Step 12: Construct a three-dimensional empty counting matrix with the same resolution as the dose calculation grid. Traverse each track point in the reverse track database, determine the grid voxel to which the track point belongs based on its spatial coordinates, and accumulate the current weight corresponding to the track point into the counting matrix element of the voxel. After accumulating all track points, normalize the accumulated count value of each voxel based on the global maximum value of the counting matrix, and output the normalized three-dimensional matrix as the full-space dose sensitivity field.

3. The Monte Carlo simulation method for accelerator target dose distribution according to claim 2, characterized in that, Step 2 includes: Step 21: Compare the standardized sensitivity coefficient of each voxel in the full-space dose sensitivity field with the preset high sensitivity threshold, and mark voxels with sensitivity coefficients greater than or equal to the threshold as high sensitivity voxels, and mark the rest of the voxels as low sensitivity voxels. Step 22: Perform connected component labeling on all high-sensitivity voxels to form a spatially continuous set of high-sensitivity regions, each of which contains its voxel list and spatial bounding box.

4. The Monte Carlo simulation method for accelerator target dose distribution according to claim 3, characterized in that, Step 2 also includes: Step 23: Based on the set of high-sensitivity regions, generate a high-resolution voxel mesh for the voxels within each high-sensitivity region according to a preset high-resolution mesh step size. Step 24: Based on the labeling results of low-sensitivity voxels, spatially adjacent low-sensitivity voxels are merged into macroscopic response surface units. Each macroscopic response surface unit stores the average sensitivity coefficient of its geometric boundary and internal voxels, and outputs an adaptive response surface mesh composed of high-resolution voxel mesh and macroscopic response surface units.

5. The Monte Carlo simulation method for accelerator target dose distribution according to claim 4, characterized in that, Step 3 includes: Step 31: A small sample set of particles, accounting for a predetermined proportion of the total planned number of particles, is emitted on the adaptive response surface grid. The dose contribution value and its spatial coordinates at multiple sampling points within each macroscopic response surface unit are recorded, and the particle flux count within each voxel in the high-resolution voxel grid is recorded. Step 32: Based on the sampling point data inside each macroscopic response surface unit, perform polynomial least squares fitting on the unit to generate response surface model coefficients, and calculate the initial flux value of each voxel based on the particle flux count in each high-resolution voxel, and output the response surface model of each macroscopic response surface unit and the initial flux distribution of the high-sensitivity region.

6. The Monte Carlo simulation method for accelerator target dose distribution according to claim 5, characterized in that, Step 4 also includes: Step 44: When the particle is located within a macroscopic response surface unit, determine whether the particle should terminate based on the sensitivity coefficient of that unit and the current local statistical uncertainty using the Russian Roulette algorithm. Step 45: If the particle is determined to have terminated, the polynomial response surface model of the macroscopic response surface unit where the particle is located is called according to the particle's current spatial coordinates to calculate the dose contribution value to each voxel in the target area and accumulate it to the dose accumulator of the corresponding voxel, and then the particle is terminated; if the particle is determined not to have terminated, the particle weight is scaled according to the survival probability of Russian roulette and then transport continues.

7. The Monte Carlo simulation method for accelerator target dose distribution according to claim 6, characterized in that, Step 6 includes: Step 61: Read the total cumulative dose and effective particle count of each voxel in the global dose accumulator, calculate the expected dose value of each voxel, and calculate the upper and lower boundaries of the confidence interval based on the statistical uncertainty of each voxel and the preset confidence level, and output the expected dose value and its confidence interval of each voxel.

8. A Monte Carlo simulation system for accelerator target area dose distribution, applied to the Monte Carlo simulation method for accelerator target area dose distribution according to any one of claims 1-7, characterized in that, include: The reverse accompaniment module is used to emit pseudo-particles in the reverse direction from the target area before the start of the Monte Carlo simulation to perform accompaniment Monte Carlo transport and generate a dose sensitivity field in the whole space. The adaptive mesh module generates an adaptive response surface mesh based on the dose sensitivity field. This mesh contains a high-resolution voxel mesh in the high-sensitivity region and macroscopic response surface elements in the low-sensitivity region. The sample fitting module is used to perform small-sample advance simulations on the adaptive response surface grid, fit the response surface model within the macroscopic response surface cells, and obtain the initial flux distribution in the high-sensitivity region. The dynamic weighting module performs the main simulation of dynamic weight allocation based on the response surface model and the initial flux distribution. The weight of each particle is determined in real time by the dynamic weighting function, which takes the response surface model and the local statistical uncertainty of the currently simulated particles as inputs. It terminates particles in the low-sensitivity region and analyzes the contribution dose through the response surface model. The local correction module is used to monitor the statistical uncertainty within the high-sensitivity voxel during the main simulation. When the statistical error in a local region exceeds a predetermined threshold, local source particles are emitted from that region, and the sample fitting module to the dynamic weighting module is repeatedly executed for the local region to perform a sub-loop resimulation. The dose output module is used to output the expected dose value and its confidence interval for each voxel.

Images