Monte carlo particle transport simulation and constrained inverse planning method for food irradiation dose field

Abstract

The present application relates to the technical field of food irradiation simulation and planning, and discloses a Monte Carlo particle transport simulation and constrained inverse planning method for food irradiation dose field. The method obtains the net mass of food, the candidate region of food and the limited-angle transmission scanning data, and establishes a mass-constrained voxel model through ray-by-ray updating and projection-maintained mass closure correction; the mean and standard error of dose response are calculated for the candidate irradiation unit, the confidence dose-constrained inverse planning is performed in the irradiation planning template under the condition of fixed equipment discretization, and the particle batches are added in a targeted manner according to the dose variance contribution of active constraint voxels; the continuous planning parameters are converted into equipment-implementation candidate combinations, and the equipment-executable irradiation planning parameters are output or the verification feedback constraint is formed after independent verification of the complete irradiation process. The present application takes into account transmission projection, food mass, Monte Carlo statistical error and equipment control resolution, and can reduce the dose cross-border risk caused by invalid particle calculation and equipment quantization.

Description

Technical Field

[0001] This invention relates to big data processing, and more particularly to the field of food irradiation simulation and irradiation process planning technology, specifically to Monte Carlo particle transport simulation and constrained inverse planning method for food irradiation dose field. Background Technology

[0003] Food irradiation typically utilizes electron beams, X-rays, or gamma rays to inhibit microorganisms, pests, and germination processes. To ensure the absorbed dose within the food is within the process window, existing technologies can establish a density model based on transmission imaging, then calculate the dose field using Monte Carlo particle transport, and improve dose uniformity by adjusting beam current, irradiation time, transport speed, or flipping method. However, industrial transport equipment is limited by space and cycle time, often resulting in transmission data acquisition only at finite angles, leading to underdetermined reconstructions prone to quality deviations. Using global scaling to compensate for weighing mass alters the measured projection, and the mass of the transport material in low-density porous foods may also lose consistency with the dose denominator. Furthermore, conventional Monte Carlo programming usually distributes particle history evenly across all irradiation conditions, failing to directly incorporate sampling errors of unit dose responses into cold and hot spot constraints, resulting in low computational resource utilization. Continuous optimization results also need to be converted into discrete dose levels achievable by the equipment; simple rounding may cause previously feasible dose boundaries to cross, while verifying all discrete combinations one by one leads to a rapid increase in the number of combinations. Existing technologies still need to address issues such as the consistency of quality in angle-limited reconstruction, the directional allocation of Monte Carlo particle batches, and the coordination between the discrete execution of continuous planning results. Summary of the Invention

[0004] This invention aims to address the challenges of simultaneously maintaining consistency between the equivalent mass density distribution of food and the measured net mass of food under limited-angle transmission conditions, the failure to utilize Monte Carlo sampling errors in resource allocation planning, and the tendency for dose crossovers to occur after continuous irradiation planning parameters are converted into discrete control values ​​for the equipment. It provides a method that couples mass-constrained voxel reconstruction, confidence-dose-constrained inverse planning, and independent verification of the complete irradiation process of the equipment.

[0005] To achieve the above objectives, the technical solution adopted by the present invention is as follows: A Monte Carlo particle transport simulation and constrained inverse programming method for food irradiation dose fields includes: acquiring the net mass of the food, candidate regions of the food, and angle-limited transmission scan data; discretizing the candidate regions of the food into candidate voxels, updating the equivalent mass density of each candidate voxel ray-by-ray based on the angle-limited transmission scan data, performing projection-preserving mass closure correction after each ray-by-ray update to obtain the equivalent mass density distribution of the food; constructing reduced-density equivalent voxels of candidate voxels with equivalent mass densities greater than 0 and less than the solid density of the food, using the product of the equivalent mass density and the voxel volume as the voxel food mass to obtain a mass-constrained voxel model; and considering factors such as the radiation source side, food carrier attitude, transport path segment, and source... Candidate irradiation units, composed of energy spectrum, source output level, and baseline delivery quantity, are used to simulate particle emission, motion, collision, scattering, and energy deposition to obtain the mean dose response and standard error of the unit dose response. Irradiation plan templates with discrete conditions for fixed equipment are enumerated, and confidence-constrained inverse programming with continuous delivery weights is performed within each irradiation plan template. Particle batches are added to the corresponding candidate irradiation units based on the dose variance contribution of active constraint voxels. The continuous irradiation planning parameters are converted into a set of candidate combinations for equipment implementation. The complete irradiation process of each candidate combination for equipment implementation is simulated using independent verification particle batches. Based on the verification results, the executable irradiation planning parameters for the equipment are output, or the confidence-constrained inverse programming is re-executed after forming verification feedback constraints.

[0006] Furthermore, acquiring angle-limited transmission scanning data includes: obtaining the net weight of the food based on the total cargo mass, the tare weight of the transport vehicle, and the registered weight of the packaging; using a depth camera to collect point clouds of the outer surface of the food loading body and determine the packaging area, the food candidate area, and the external air area; acquiring transmission intensity images, dark field intensity images, and air field intensity images at -50°, -25°, 0°, 25°, and 50°, and obtaining the initial attenuation value after dark field correction and air field correction; measuring food equivalent material calibration blocks of multiple thickness levels using a beam blocking array, and establishing a scattering distribution correction table according to the scanning angle, detection channel, and the attenuation interval where the initial attenuation value is located; placing the food equivalent material calibration blocks in the same packaging and transport vehicle background as the food loading body, establishing a background matching calibration curve between the scattering-corrected transmission measurement value and the newly added food equivalent mass thickness according to the scanning angle and detection channel, and obtaining the measured food equivalent mass thickness using the scattering distribution correction table and the background matching calibration curve.

[0007] Furthermore, the ray-by-ray update includes: determining the mass capacity range based on the volume, lower density limit, and upper density limit of the food candidate voxels; when the net mass of the food is within the mass capacity range, the result obtained by dividing the net mass of the food by the total volume of the food candidate voxels is set as the initial equivalent mass density of the food; when the net mass of the food is outside the mass capacity range, modeling condition conflict information is output; establishing a projection matrix based on the ray's travel length within each food candidate voxel; for each ray, calculating the projection residual between the measured equivalent mass thickness of the food and the sum of the products of each travel length and the corresponding current equivalent mass density of the food; for each food candidate voxel traversed by the current ray, dividing the product of the projection residual, the corresponding travel length, and a relaxation coefficient greater than 0 and less than 2 by the sum of the squares of the travel lengths corresponding to the current ray and the dimensionless positive stability quantity, and superimposing the resulting density increment onto the current equivalent mass density of the food, and limiting the superposition result to between the corresponding lower density limit and upper density limit.

[0008] Furthermore, the projection-preserving quality closure correction includes: weighting the projection changes of each ray according to the reciprocal of the standard deviation of the corresponding ray's repeated calibration; using the weighted sum of the squares of the weighted projection changes caused by the voxel density correction and the squares of the differences in the voxel 6-neighborhood of the corrected food equivalent mass density as the objective function; and using the constraints that the total corrected food candidate voxel mass equals the net food mass and that the corrected food equivalent mass density is between the corresponding lower and upper density limits. The primal-dual interior-point method is used to solve the convex quadratic programming problem. Ray-by-ray updates and projection-preserving quality closure corrections are alternately executed until the weighted projection... The ratio of the residual norm to the weighted measured food equivalent mass thickness norm is less than or equal to the preset projection residual threshold, and the relative difference between the total mass of food candidate voxels and the net mass of the food is less than or equal to the preset mass residual threshold; when the preset maximum number of reconstruction rounds is reached and any ratio exceeds the corresponding threshold, modeling condition conflict information is output; food candidate voxels with an equivalent mass density of 0 are assigned an air material label, and food candidate voxels with an equivalent mass density greater than 0 are established as particle transport materials according to the corresponding food equivalent mass density and food element mass fraction, and all deposition energy within the voxels is included in the food deposition energy.

[0009] Furthermore, the electron beam equipment uses a reference beam current charge as the reference delivery quantity, while the continuous X-ray and gamma-ray equipment use a reference irradiation time as the reference delivery quantity. Candidate voxels whose voxel food quality reaches a preset dose measurement quality threshold are grouped into a food dose-constrained voxel set. Each candidate irradiation unit generates at least 30 batches of particles with equal historical numbers using mutually separate random number subsequences. Based on the rigid body transformation relationship of the food carrier attitude, the energy deposition position is inversely transformed to the food's body coordinate system. The food deposition energy is divided by the voxel food quality and the reference delivery quantity to obtain the reference delivery quantity dose response, and the unit dose response mean and unit dose response standard error are calculated. During the continuous delivery weighting solution, the material state, geometric pose, and source term energy spectrum of the candidate irradiation unit are kept fixed. When adding particle batches, subsequent random number subsequences are called and the unit dose response mean and unit dose response standard error of the corresponding candidate irradiation unit on all food dose-constrained voxels are updated.

[0010] Furthermore, the confidence dose-constrained inverse programming includes: enumerating a finite number of irradiation plan templates based on the radiation source side separation order, food carrier flipping order, number of irradiation passes, and source output level; setting the range of values ​​for the continuous delivery weight according to the range of delivery speed or residence time, where the continuous delivery weight represents the ratio of the total delivery amount of each irradiation pass on the same candidate irradiation unit to the baseline delivery amount; before the start of each planning round, determining the conservative confidence coefficient of the Student's t-distribution based on the minimum value of the degrees of freedom obtained by subtracting 1 from the number of particle batches in each candidate irradiation unit in the current irradiation plan template, the planning tail parameters, the number of food dose-constrained voxels, and the maximum number of planning rounds; and multiplying the unit dose response mean by the continuous delivery weight... The mean of the planned dose is determined, the square root of the sum of the squares of the products of the unit dose response standard error and the continuous delivery weights is determined as the standard error of the planned dose, the sum of the products of the mean of the planned dose and the conservative confidence coefficient of the student t-distribution and the standard error of the planned dose is determined as the upper bound of the planned Monte Carlo statistical confidence dose, and the subtraction result is determined as the lower bound of the planned Monte Carlo statistical confidence dose. The food dose constraint is that the lower bound of the planned Monte Carlo statistical confidence dose reaches the lower limit of the target dose and the upper bound of the planned Monte Carlo statistical confidence dose is within the upper limit of the target dose. The constraints of equipment power, conveying speed and delivery weight changes of adjacent conveying path segments are added, and the primal-dual interior point method is used to solve the continuous second-order cone programming.

[0011] Furthermore, the planning objective of the continuous second-order cone programming is the sum of the products of continuous delivery weights and equipment time conversion factors or equipment energy consumption conversion factors, and the weighted sum of the absolute differences in continuous delivery weights between adjacent transport path segments. When food dose constraints are met, food dose constraint voxels whose non-negative difference between the planned Monte Carlo statistical confidence dose lower bound and the target dose lower bound, or whose non-negative difference between the target dose upper bound and the planned Monte Carlo statistical confidence dose upper bound, is within a preset judgment margin are determined as active constraint voxels. The square of the product of the unit dose response standard error and the continuous delivery weight is determined as the dose variance contribution. Select the fewest candidate irradiation units whose cumulative contribution reaches the preset variance ratio according to the descending order of dose variance contribution. Add particle batches to the union of the candidate irradiation units selected for each active constraint voxel. Update the mean dose response of the corresponding whole unit, the standard error of the whole unit dose response, and the conservative confidence coefficient of the Student t-distribution, and then resolve the continuous second-order cone programming. When the continuous second-order cone programming forms a conflict state, add a unified dose relaxation amount to obtain the minimum relaxation solution, and add particle batches according to the cross-boundary voxels of the minimum relaxation solution. When the maximum number of planning rounds or the maximum number of particle batches is reached, record the current irradiation plan template as the planning conflict template.

[0012] Furthermore, generating a set of candidate combinations for device implementation includes: obtaining the lower and upper adjacent realizable values ​​of each continuous delivery weight based on the device control resolution, and establishing a search tree for selecting adjacent realizable values ​​layer by layer; determining the lower bound of a node by the sum of the cumulative weighted absolute deviation between the selected adjacent realizable values ​​and the corresponding continuous delivery weights, and the weighted absolute deviations formed by the remaining continuous delivery weights using adjacent realizable values ​​with smaller deviations; establishing a table of verified combinations for storing completed verification combinations, expanding the search tree nodes in ascending order according to the lower bound of the nodes using a priority queue, and retaining candidate combinations for device implementation that satisfy the constraints of device power, delivery speed, changes in delivery weights of adjacent delivery path segments, and device synchronization at the leaf nodes, and that are outside the table of verified combinations, until the number of candidate combinations for device implementation reaches the preset maximum number of combinations, the priority queue is cleared, or the number of expanded nodes reaches the preset maximum number of expanded nodes.

[0013] Furthermore, before verification, a maximum total number of verifications is set that covers all irradiation plan templates, all candidate combinations of equipment implementation, and all feedback loops, allowing for statistical verification decisions. For the first verification of each candidate combination of equipment implementation, a complete process verification particle batch of at least 30 identical particles is used. The complete irradiation process is simulated according to the corresponding radiation source output, delivery speed, residence time, food carrier attitude, radiation source side order, food carrier flipping order, and number of irradiation passes. Each candidate combination of equipment implementation and each verification uses a mutually separate random number subsequence. The overall verification tail probability is divided by the number of food dose constraint voxels and the maximum verification total. The product of the number of times and 2 is used to obtain the one-sided tail probability. The student t-distribution quantile is determined by subtracting the one-sided tail probability from 1 to obtain the cumulative probability and subtracting 1 from the number of batches of particles verified in the complete process. The upper bound of the verification Monte Carlo statistical confidence dose is obtained by adding the voxel dose mean to the student t-distribution quantile and the voxel dose standard error. The lower bound of the verification Monte Carlo statistical confidence dose is obtained by subtracting the corresponding product from the voxel dose mean. Each calculation of the verification Monte Carlo statistical confidence dose bound increases the verification count by 1. When all food dose constraints are met, the output device can execute irradiation planning parameters, predict cold spot locations, predict hot spot locations, and recommend dosimeter placement locations.

[0014] Furthermore, the candidate combinations of equipment implementations that have completed verification are written into the verified combination table. When all candidate combinations of equipment implementations in this round form cross-boundary results, the ratio of the positive difference between the lower limit of the target dose and the lower limit of the verification Monte Carlo statistical confidence dose to the lower limit of the target dose, and the ratio of the positive difference between the upper limit of the verification Monte Carlo statistical confidence dose and the upper limit of the target dose to the upper limit of the target dose are calculated. The maximum value among all ratios is determined as the maximum normalized dose cross-boundary, and the candidate combination of equipment implementations with the smallest maximum normalized dose cross-boundary is selected as the feedback target combination. Food dose constraint voxels in the feedback target combination that simultaneously exhibit cold point cross-boundary and hot point cross-boundary are determined as statistical width cross-boundary voxels, and a subsequent preset number of complete process verification particle batches are added. The cumulative complete process cross-boundary results are then used to determine the statistical width cross-boundary voxels. The entire irradiation process verification is re-executed for the particle batch verification. When the number of cross-boundary voxels is 0, verification feedback constraints are formed for cold cross-boundary voxels, where the average planned dose of subsequent rounds reaches the sum of the average planned dose of the feedback target combination and the positive feedback dose step size. For hot cross-boundary voxels, verification feedback constraints are formed, where the average planned dose of subsequent rounds is within the difference between the average planned dose of the feedback target combination and the positive feedback dose step size. The verification feedback constraints of this round are replaced with the verification feedback constraints of the previous round, and the confidence dose constraint reverse planning is re-executed. When the verification count reaches the maximum total number of verifications or all irradiation plan templates are planning conflict templates, all equipment implementation of candidate combination verification is stopped. Irradiation plan templates that have not yet been verified are recorded as planning conflict templates, and the planning conflict results corresponding to the cross-boundary voxels and equipment constraints are output.

[0015] By employing the above-described technical content, the present invention achieves the following beneficial effects: This invention uses the weighing results of food carriers, angle-limited transmission scanning data, and the control capabilities of irradiation equipment as a unified data foundation. After updating each ray, it implements projection-preserving mass closure correction, ensuring that mass compensation is directed towards density directions that have minimal impact on the measured projection and suppressing isolated density spikes. It simultaneously determines the particle transport material density and voxel food mass using the equivalent mass density of the food, thus unifying the weighing mass, transport material mass, and absorbed dose denominator. After obtaining the mean and standard error of the unit dose response, Monte Carlo sampling errors are incorporated into the two-sided confidence dose constraints of cold and hot spots. The particle batches of corresponding candidate irradiation units are increased according to the dose variance contribution of active constraint voxels, concentrating computational resources on statistical error sources that determine the feasibility of the planning. By pre-fixing the equipment discrete conditions and solving the continuous programming parameters within the template, the solvability of the optimization problem and the equipment operating conditions are considered. The continuous programming parameters are converted into candidate combinations for equipment realization through minimum deviation combination search, and then the complete irradiation process is simulated using non-overlapping random number subsequences. When verification fails, the feedback target combination with the smallest cross-limit of the maximum normalized dose is selected. When the statistical interval is too wide, additional batches of verification particles are added. This can reduce the risk of erroneous feasibility or human conflict caused by equipment quantization, repeated statistical judgment and opposite feedback, and output cold spots, hot spots and recommended dosimeter placement locations. Attached Figure Description

[0016] Figure 1 A schematic diagram illustrating the principle of angle-limited transmission scanning and food candidate region determination.

[0017] Figure 2 This is a schematic diagram illustrating the characteristics of scattering distribution correction and background matching calibration, where, Figure 2 (a) in the figure shows the characteristic graph of the scattering correction ratio as a function of the initial attenuation value at different scanning angles. Figure 2 (b) in the figure is a background matching calibration characteristic diagram between the attenuation integral after scattering correction and the equivalent mass thickness of the newly added food.

[0018] Figure 3 This is a comparison curve of the equivalent mass density of food after ray-by-ray updating and after projection-preserving mass closure correction.

[0019] Figure 4 A schematic diagram illustrating the particle transport and food deposition energy measurement principle of equivalent voxels for dedensed foods.

[0020] Figure 5 This is a schematic diagram illustrating the principle of inverse transformation of candidate irradiation units and energy deposition locations to the food body coordinate system.

[0021] Figure 6A schematic diagram illustrating the Monte Carlo statistical confidence dose bound and active constraint voxel determination characteristics.

[0022] Figure 7 A schematic diagram illustrating the dose variance contribution of active constraint voxels and the selection characteristics of particle batch addition targets.

[0023] Figure 8 This diagram illustrates the relationship between continuous irradiation planning parameters, adjacent feasible combinations, and verification feedback constraints. Detailed Implementation

[0024] A Monte Carlo particle transport simulation and constrained inverse programming method for food irradiation dose fields includes: acquiring the net mass of the food, candidate regions of the food, and angle-limited transmission scan data; discretizing the candidate regions of the food into candidate voxels, updating the equivalent mass density of each candidate voxel ray-by-ray based on the angle-limited transmission scan data, performing projection-preserving mass closure correction after each ray-by-ray update to obtain the equivalent mass density distribution of the food; constructing reduced-density equivalent voxels of candidate voxels with equivalent mass densities greater than 0 and less than the solid density of the food, using the product of the equivalent mass density and the voxel volume as the voxel food mass to obtain a mass-constrained voxel model; and considering factors such as the radiation source side, food carrier attitude, transport path segment, and source... Candidate irradiation units, composed of energy spectrum, source output level, and baseline delivery quantity, are used to simulate particle emission, motion, collision, scattering, and energy deposition to obtain the mean dose response and standard error of the unit dose response. Irradiation plan templates with discrete conditions for fixed equipment are enumerated, and confidence-constrained inverse programming with continuous delivery weights is performed within each irradiation plan template. Particle batches are added to the corresponding candidate irradiation units based on the dose variance contribution of active constraint voxels. The continuous irradiation planning parameters are converted into a set of candidate combinations for equipment implementation. The complete irradiation process of each candidate combination for equipment implementation is simulated using independent verification particle batches. Based on the verification results, the executable irradiation planning parameters for the equipment are output, or the confidence-constrained inverse programming is re-executed after forming verification feedback constraints.

[0025] The present invention will be further described below with reference to the accompanying drawings and specific embodiments. The numerical values ​​used herein are intended to indicate repeatable implementation conditions and do not constitute a limitation on the type of equipment, food category, or processing scale.

[0026] This embodiment focuses on granular dehydrated food loaded in a rectangular turnover box, using a continuous X-ray device as the irradiation equipment. The food carrier is an integral unit formed by the packaging and the food inside, and is fixed to a conveyor with four positioning marks. A conveyor coordinate system, an equipment coordinate system, and a food-based coordinate system are established: the conveyor coordinate system is used to unify the angle-limited transmission scanning and surface point cloud; the equipment coordinate system is fixed to the irradiation equipment; and the food-based coordinate system translates and rotates along with the food carrier. All three coordinate systems are right-handed rectangular coordinate systems, with the origin of the food-based coordinate system set at the geometric center of the food candidate region.

[0027] The total weight of the cargo is 8.88 kg, the tare weight of the conveyor is 0.52 kg, and the registered weight of the packaging is 0.36 kg. Subtracting the tare weight of the conveyor and the registered weight of the packaging from the total weight yields a net weight of 8.00 kg for the food. The effective internal dimensions of the packaging area are 300 mm × 220 mm × 150 mm. A depth camera collects point clouds above the food container and from two oblique sides, retaining at least 150,000 effective depth points at each viewpoint. The rigid body transformation relationship is first solved using any three non-collinear positioning markers from the four positioning markers, and then the registration error is checked using the fourth positioning marker. The point cloud is accepted if the root mean square error of the registered positioning markers is no greater than 0.8 mm; otherwise, it is re-acquired.

[0028] The packaging area is determined based on the wall thickness, internal cavity, and partition structure data of the packaging material. The space enclosed by the inner surface of the packaging area, after deducting the fixed partitions and supports, is defined as the food candidate area, and the remaining space is defined as the external air area. In this embodiment, the turnover box does not have fixed partitions or supports, and the food candidate area covers an internal cavity of 300 mm × 220 mm × 150 mm. The food candidate area is used to characterize the planned space that may contain food entities or food pores, and is not used to characterize the actual microstructure uniquely recovered from five scanning angles. Constraining the food candidate area with external surface point clouds can exclude ray paths outside the food carrier; retaining the low-density space inside the food candidate area allows subsequent mass-constrained voxel reconstruction to represent the porosity effect as a reduced-density food equivalent voxel.

[0029] like Figure 1As shown, the imaging X-ray source is located sequentially at -50°, -25°, 0°, 25°, and 50°. A linear array detector receives the transmitted signal on the other side of the food carrier, with five scanning angles covering a 100° angular limit. The imaging X-ray source uses a 320 kV peak voltage and a 5 mA tube current. The distance from the source to the center of the food carrier is 750 mm, and the distance from the source to the linear array detector is 1250 mm. The linear array detector has 1024 detection channels with a channel spacing of 0.8 mm. Twenty frames of transmission intensity images are acquired for each scanning angle, and 50 frames each of dark field intensity and open field intensity images are acquired. The arithmetic mean of all frames from the same detection channel is taken to reduce the impact of instantaneous electronic noise and source output fluctuations on the initial attenuation value.

[0030] Dark field correction, empty field correction, and initial attenuation value are calculated according to the following formula:

[0031] In the formula, Indicates the X-ray exposure number; Indicates the detection channel number; This indicates the transmission intensity of the corresponding X-ray exposure and detection channels; This indicates the dark field intensity of the corresponding detection channel; This indicates the air field intensity of the corresponding detection channel; Indicates the initial transmittance; Indicates the initial decay value; This indicates a dimensionless positive number to prevent the logarithmic independent variable from being zero; in this implementation, it is taken as 0.000001. When the denominator is less than 0.5% of the average intensity of the air field, the corresponding detection channel is marked as an invalid detection channel, and interpolation is performed using adjacent valid detection channels to avoid amplifying noise with a very small denominator.

[0032] Scattering correction was performed using a pre-established scattering distribution correction table based on food equivalent material calibration blocks and a beam blocking array. The food equivalent material calibration blocks used elemental mass fractions close to those of the food being treated, with thickness levels set to 0 mm, 20 mm, 40 mm, 60 mm, 80 mm, 100 mm, 120 mm, and 150 mm. The beam blocking array consisted of tungsten discs with a diameter of 3 mm and a center-to-center distance of 24 mm. The residual signal after subtracting the dark field within the shadow of the tungsten discs was used as the local scattering signal. Piecewise cubic interpolation along the probe channel direction, maintaining non-negativity, was used to reconstruct the scattering distribution at the unblocked location. Each thickness level was measured 10 times, and the standard deviation was recorded.

[0033] like Figure 2As shown in (a), the scattering correction ratio changes with the initial attenuation value and the scanning angle. The scattering distribution correction table uses the scanning angle, the probe channel, and the attenuation interval containing the initial attenuation value as a three-level index, with each attenuation interval divided into 0.25 intervals based on the initial attenuation value. After calling the corresponding table entry, the scattering estimate is subtracted from the dark field correction intensity, and then open field normalization and negative logarithmic transformation are performed. Using the initial attenuation values ​​obtained from dark field correction and open field correction to select table entries, instead of selecting entries again based on the scattering correction result, avoids changing the index and causing circular lookups. The table entries are interpolated using a linear interpolation method with respect to the scanning angle and probe channel, respectively, and a monotonic interpolation method with respect to the attenuation interval. The interpolation result is limited to non-negative values ​​and less than the dark field correction intensity.

[0034] like Figure 2 As shown in (b), the background matching calibration curve is established after scattering correction. Food equivalent material calibration blocks of various thickness grades are placed in a background of packaging and transport vehicles identical to the food packaging. The newly added known areal density of the calibration blocks is used as the new food equivalent mass thickness. For each scanning angle and detection channel, a piecewise cubic interpolation curve maintaining monotonicity between the attenuation integral after scattering correction and the new food equivalent mass thickness is established. During the actual scan, the corresponding attenuation integral after scattering correction is substituted into the curve to obtain the measured food equivalent mass thickness. The transmission relationship of multi-energy X-rays undergoes beam hardening with material thickness. The attenuation of packaging and transport vehicles cannot be directly subtracted from the negative logarithmic projection without linearization. Background matching calibration allows the combined effects of the fixed structure and beam hardening to enter the same calibration mapping, thus the obtained quantity can be used as the equivalent mass thickness along the food path.

[0035] The food candidate region was discretized into cubes with sides of 10 mm, forming 9900 food candidate voxels, each with a volume of 1 cubic centimeter. The density of the food entity was set to 1050 kg / m³, with a lower limit of 0 kg / m³ and an upper limit of 1200 kg / m³. Before reconstruction, the mass capacity was checked using the following formula:

[0036] In the formula, Represents a set of candidate voxels for food; Indicates the voxel number of the food candidate; This indicates the volume of the corresponding food candidate voxel; and These represent the lower and upper limits of density for the corresponding food candidate voxels, respectively. This represents the net weight of the food. In this embodiment, the weight capacity range is 0 kg to 11.88 kg. Since the net weight of the food falls within this range, the initial equivalent mass density of the food is set to 808.08 kg / m³, obtained by dividing the net weight of the food by the total volume of the candidate voxels. If the net weight of the food is outside the weight capacity range, the packaging area, the candidate food area, the weighing data, and the density boundary are checked, and modeling condition conflict information is output; ray-by-ray updating is not performed.

[0037] The traversal length of each ray within each candidate food voxel is calculated using a face-by-face intersection method between the ray and the voxel boundary. For any ray, the incident and exit parameters of the ray and the three sets of opposite surfaces of the voxel are calculated. The positive difference between the maximum and minimum incident parameters is taken as the traversal segment; if no positive difference exists, the traversal length is 0. The projection matrix is ​​stored in a sparse row format, with each row recording only the food candidate voxel number and traversal length traversed by the current ray, thus reducing the storage requirements of the angle-limited voxel model. The measured equivalent mass thickness of the food is used to form the projection vector according to the same ray sequence.

[0038] Calculate the projection residuals sequentially according to the valid ray numbers, and only update the food candidate voxels that the current ray passes through:

[0039]

[0040] In the formula, This represents the projection residual of the current ray; This indicates the measured equivalent mass thickness of the food sample under the current radiation. This indicates that the current ray is at the number... The travel length of the food candidate body element is converted to meters; This represents the current food equivalent mass density of the corresponding food candidate voxel; This represents the relaxation coefficient, which is set to 0.8 in this implementation. This represents a positively stable quantity with the same dimensions as the square of the travel length; in this implementation, it is taken as 0.0000000001 square meters. This indicates a projection operator that limits the updated result to the density range between the corresponding lower and upper density limits. The projection residual is expressed in kilograms per square meter, and the travel length in meters, and the calculated density increment is in kilograms per cubic meter. A relaxation coefficient between 0 and 2 controls the ray-by-ray correction amplitude, while a positive stability factor prevents extremely short paths from causing the denominator to approach zero.

[0041] One ray-by-ray update is completed after all effective rays have been processed, followed by a projection-preserving mass closure correction for the equivalent mass density of the intermediate food. The correction problem is:

[0042]

[0043] In the formula, This represents the equivalent mass density vector of intermediate food after one round of ray-by-ray updates; Represents the equivalent mass density correction vector for food; This represents the projection matrix consisting of the travel lengths; This represents a weighted matrix, where the diagonal elements are the reciprocals of the standard deviation of the repeated calibration of the equivalent mass thickness of the food corresponding to the ray. This represents the normalized voxel 6-neighborhood difference matrix. Each row outputs the difference in the equivalent mass density of two adjacent food candidate voxels divided by 1000 kg per cubic meter. The space regularization weight is 0.015 in this implementation. Normalization ensures that both objective terms are dimensionless. The mass equation ensures that the sum of the masses of the food candidate voxels is precisely closed, and the density interval constraint excludes non-physical densities.

[0044] The original-dual interior-point method is used to solve the above convex quadratic programming problem. The weighted projection variation term suppresses the change in measured rays caused by mass compensation, distributing the mass difference along the food equivalent mass density correction direction, which has a smaller impact on the limited-angle transmission projection; the normalized voxel 6-neighborhood difference term suppresses isolated spikes in the corrected food equivalent mass density. For example... Figure 3 As shown, the density curves after ray-by-ray updating and the projection-preserving mass closure correction only exhibit constrained changes at local weak observation locations, with the total mass of food candidate voxels closing from 7.84 kg to 8.00 kg. The projection preservation effect is constrained by the convergence condition of the weighted projection residuals. Compared to multiplying all food candidate voxels by the same proportion, this method reduces projection mismatch while satisfying the weighing mass requirement.

[0045] Ray-by-ray updates and projection-preserving quality closure corrections are performed alternately, with the convergence condition being:

[0046] In the formula, This represents the current equivalent mass density vector of the food. This represents the equivalent mass thickness vector of the measured food product. This represents the preset projection residual threshold, which is set to 0.03 in this implementation. The preset quality residual threshold is 0.00001 in this implementation. The preset maximum number of reconstruction rounds is 60. After the first round of ray-by-ray update, the total mass of the food candidate voxels is 7.84 kg, and the weighted projection residual norm ratio is 0.137. After 18 rounds of alternating calculations, the total mass of the food candidate voxels is 8.0000 kg, the weighted projection residual norm ratio is 0.019, and the relative mass difference is 0.0000026, thus obtaining the equivalent mass density distribution of food for dose-oriented calculation. When 60 rounds are reached and any residual exceeds the corresponding threshold, modeling condition conflict information is output.

[0047] A mass-constrained voxel model is established based on the equivalent mass density distribution of food. Candidate voxels of food with an equivalent mass density of 0 are assigned the label of air material; candidate voxels of food with an equivalent mass density greater than 0 but less than the solid density of the food are constructed as reduced-density food equivalent voxels; candidate voxels of food with an equivalent mass density equal to or exceeding the solid density of the food are assigned the label of food material. For candidate voxels of food with an equivalent mass density greater than 0, particle transport material is established using the mass fraction of food elements corresponding to the food category and the current equivalent mass density of the food; the corresponding candidate voxels are not changed to air material.

[0048] The quality of voxel-based food and the voxel absorption dose are calculated using the following formula:

[0049] In the formula, This indicates the voxel quality of the corresponding food candidate voxel; Indicates the absorption dose of the corresponding food candidate voxel; This represents the food deposition energy accumulated to the corresponding food candidate voxel during particle transport. When calculating the absorbed dose, the food deposition energy is converted to joules, and the voxel food mass is converted to kilograms. The de-density food equivalent voxel can be considered as a homogeneous equivalent material of porous food; the porosity effect is reflected in the overall density, and all deposition energy within the voxel is included in the food deposition energy. Therefore, the reconstruction mass, transport material mass, and dose denominator remain consistent. Figure 4 The process of incident particles passing through food material voxels and de-density food equivalent voxels is shown, as well as how food deposition energy is attributed to metric voxels.

[0050] Food candidate voxels with a mass of 0.05 grams or less were grouped into a food dose-constrained voxel set. Food candidate voxels with a mass less than 0.05 grams were retained in the particle transport geometry model and participated in scattering and energy loss, but were not included as food dose-constrained voxels to avoid numerical amplification of voxel absorbed dose due to extremely small food masses near the air boundary. Particle transport materials were established using registered material compositions and densities for the packaging area, transport vehicle, and external air area, respectively.

[0051] This implementation uses a 7.5 MeV continuous X-ray device, with source output levels set to 70% and 100% of rated output. The baseline delivery time is 0.20 seconds. The delivery path is divided into 12 25mm segments along a 300mm direction. The food container attitude includes a 0° attitude and an attitude rotated 180° around the vertical axis. The radiation source sides include radiation source side 1 and radiation source side 2. Ninety-six candidate irradiation units are formed by two radiation source sides, two food container attitudes, 12 delivery path segments, one set of source term energy spectra, two source output levels, and corresponding baseline delivery times. Each fixed source output level's irradiation plan template calls upon 48 of these candidate irradiation units.

[0052] The mass-constrained voxel model, packaging structure, transport vehicle structure, source focus, collimator, shielding structure, and source term energy spectrum are converted into a particle transport geometry model. For each transport path segment, the initial particle position is continuously and uniformly sampled along the direction of movement of the food packaging carrier, so that the unit dose response includes the average influence of the moving position throughout the entire transport path segment. Photon transport includes photoelectric absorption, Compton scattering, and electron-electron pair effects; secondary electron and positron transport includes ionization, bremsstrahlung, and multiple scattering. The photon cutoff energy is set to 5 kEV, the electron cutoff energy to 10 kEV, and the maximum electron transport step size to 2 mm. Residual energy below the cutoff energy is locally deposited within the current material voxel. Optionally, the above physical process settings can be replaced with validated particle transport procedures and cross-section databases.

[0053] like Figure 5 As shown, different radiation source sides, food carrier orientations, and transport path segments constitute different candidate irradiation units. The energy deposition locations in the equipment coordinate system are incorporated into a unified food body coordinate system through inverse rigid body transformation.

[0054] In the formula, This represents the energy deposition location vector in the device coordinate system; This represents the energy deposition location vector in the food's body coordinate system; This represents the food carrier attitude rotation matrix corresponding to the candidate irradiation unit; Represents the translation vector of the corresponding candidate irradiation unit; superscript This indicates the matrix transpose. The rotation matrix is ​​an orthogonal matrix, and its transpose is equal to its inverse. After the inverse transformation, the food deposition energy generated by the posture of each food carrier falls into the same food candidate voxel number, avoiding direct superposition in the device coordinate system and causing food position mismatch.

[0055] Each candidate irradiation unit initially generates 30 particle batches with equal history numbers, each batch containing 200,000 primary particle histories. Equal history number particle batches refer to batches within the same candidate irradiation unit having the same number of primary particle histories. The random number subsequence is mapped from the candidate irradiation unit number, particle batch number, and particle history number to a non-overlapping subsequence; different candidate irradiation units do not share subsequences, thus the planned dose standard error can be synthesized from independent variance components. The computing workstation is configured with two 64-core processors, 256 gigabytes of memory, and two 48-gigabyte graphics processors. Parallel processing is implemented between candidate irradiation units, and within each unit, parallel processing is performed on a batch-by-batch basis.

[0056] The baseline delivery dose response for a single particle batch is calculated using the following formula:

[0057] In the formula, Indicates the number is The candidate irradiation unit Each batch of particles is numbered as follows The baseline delivery dose response of food candidate voxels; This indicates the number of physical source particles corresponding to the baseline delivery quantity obtained from the device output calibration. This represents the cumulative food deposition energy of the corresponding particle batch within the corresponding food candidate voxel; This indicates the primary particle history number for the corresponding particle batch. After converting energy to joules and voxel food mass to kilograms, the unit for baseline delivery dose response is gray per baseline delivery.

[0058] The mean of unit dose response and the standard error of unit dose response are calculated according to the following formulas:

[0059]

[0060] In the formula, This represents the mean dose response of a single unit. Indicates the standard error of the unit dose response; This represents the number of particle batches for the corresponding candidate irradiation unit. During the continuous delivery weighting process, the material state, geometric pose, and source term energy spectrum of each candidate irradiation unit are kept constant, and the baseline delivery quantity is scaled only according to the continuous delivery weighting. Each time a particle batch is added, the subsequent random number subsequence is called, and the mean of the array dose response and the standard error of the array dose response of the corresponding candidate irradiation unit on all food dose-constrained voxels are updated.

[0061] Two radiation source identification sequences are defined: radiation source identification 1 followed by radiation source identification 2, and radiation source identification 2 followed by radiation source identification 1. Two food carrier flipping sequences are defined: 0° attitude followed by 180° attitude, and 180° attitude followed by 0° attitude. These two sequences, along with the food carrier flipping sequences, one or two irradiations, and 70% or 100% source output settings, are combined to form 16 irradiation plan templates. Each template pre-fixes discrete equipment conditions, and only the continuous delivery weight is calculated within the template. The number of irradiations and the source output setting are not considered variables in the continuous second-order cone programming. When using segmented dwell time delivery, the dwell time per irradiation is set to 0.10 seconds to 0.32 seconds. The number of irradiations represents the number of complete executions of the corresponding radiation source identification sequence and food carrier flipping sequence; the continuous delivery weight represents the ratio of the total delivery volume of each irradiation cycle on the same candidate irradiation unit to the baseline delivery volume. The continuous delivery weight ranges from 0.50 to 1.60 within a single irradiation template, and from 1.00 to 3.20 within a two-irradiation template; the difference in continuous delivery weight between adjacent delivery path segments shall not exceed 0.20. The lower limit of the target dose is 3.0 kGy, and the upper limit of the target dose is 7.0 kGy.

[0062] Before each planning round, the minimum particle batch number minus 1 within the current irradiation plan template is used as the degree of freedom. The planning tail parameter 0.05 of the current irradiation plan template is divided by the product of the number of food dose constraint voxels, the maximum planning round number 8, and 2 to obtain the planning one-sided tail probability. Then, the conservative confidence coefficient of the Student's t-distribution is determined based on the cumulative probability obtained by subtracting the planning one-sided tail probability from 1 and the corresponding degree of freedom. The conservative confidence coefficient of the Student's t-distribution is fixed during the continuous second-order cone programming solution of this round and is only updated after adding particle batches and before entering the next round of solution, to avoid the confidence coefficient from destroying the second-order cone structure due to changes in the continuous delivery weights. The planned dose and the planning Monte Carlo statistical confidence dose bounds are:

[0063]

[0064] In the formula, This indicates the sequential delivery weight of the corresponding candidate irradiation unit; This represents the average planned dose of the corresponding food dose-restrained voxel; Indicates the standard error of the planned dose; Indicates the first Conservative confidence coefficients for the student t-distribution used in round planning; and These represent the lower and upper bounds of the planned Monte Carlo statistical confidence dose, respectively. During the planning phase, the lower bound of the confidence dose is no lower than the target dose lower limit, and the upper bound of the confidence dose is no higher than the target dose upper limit, serving as the food dose constraint. The planned Monte Carlo statistical confidence dose boundaries are used to plan, screen, and determine the subjects for additional particle batches; the final statistical determination is based on the independent verification results of the complete irradiation process.

[0065] The planning objective is a weighted sum of the sum of equipment times corresponding to continuous delivery weights and the absolute difference of continuous delivery weights in adjacent transport path segments. The absolute difference is represented by two linear inequalities and non-negative auxiliary variables. The equipment time conversion factor is 0.20 seconds, and the smoothing weight is 0.04 seconds. When the planning objective uses equipment energy consumption, the equipment energy consumption conversion factor is the sum of the radiation source input power, transport drive power, and auxiliary equipment power of the corresponding candidate irradiation unit under the baseline delivery quantity. In addition to food dose constraints, the following constraints are added: radiation source power limit, dwell time range, transport speed range, changes in delivery weights in adjacent transport path segments, and equipment synchronization constraints. The equipment synchronization constraint stipulates that only one radiation source side performs irradiation within the same transport path segment, and the dwell time begins after the food carrier positioning completion signal and the radiation source stable output signal are simultaneously valid. The primal-dual interior point method is used to solve the continuous second-order cone programming. For each irradiation plan template, if the food dose constraint holds and the product of the standard error of all planned doses and the conservative confidence coefficient of the Student's t-distribution does not exceed 0.05 kGy, the template feasible solution is recorded; if multiple templates are feasible, the template feasible solution with the smallest planning objective value is selected as the continuous irradiation planning parameter.

[0066] Figure 6 The vertical position of the food component is used as the horizontal axis to show the planned mean dose, the upper bound of the planned Monte Carlo statistical confidence dose, and the lower bound of the planned Monte Carlo statistical confidence dose. When the food dose constraint is valid, food dose constraint voxels with a non-negative difference between the lower bound of the planned Monte Carlo statistical confidence dose and the lower limit of the target dose, or a non-negative difference between the upper limit of the target dose and the upper bound of the planned Monte Carlo statistical confidence dose, are defined as active constraint voxels. In the figure, the voxels closer to the lower limit of the target dose are the lower limit active constraint voxels, and the voxels closer to the upper limit of the target dose are the upper limit active constraint voxels. For each active constraint voxel, the dose variance contribution of each candidate irradiation unit is calculated:

[0067] In the formula, This represents the dose variance contribution of the corresponding candidate irradiation unit to the corresponding active constrained voxel. Candidate irradiation units with the fewest cumulative contributions reaching 80% of the preset variance ratio are selected in descending order of dose variance contribution. The union of all selected active constrained voxels is used as an additional object, and each additional object adds 5 batches of particles with equal historical data, with a maximum batch size of 70. Figure 7 The cumulative contribution of the top three candidate irradiation units exceeds 80%, therefore, their particle batches are prioritized for addition. When using the above selection method, the newly added particle history focuses on the sources of statistical error that determine the confidence dose boundary, and low-contribution candidate irradiation units are not averaged up.

[0068] When consecutive second-order cone planning results in a conflict, a unified dose relaxation is introduced: the sum of the confidence dose lower bound and the unified dose relaxation is constrained to be no less than the target dose lower bound, and the difference between the confidence dose upper bound and the unified dose relaxation is constrained to be no more than the target dose upper bound, with the minimum unified dose relaxation being the priority objective. Additional targets are selected based on the cross-boundary voxels of the minimum relaxation solution and their corresponding dose variance contributions. If a conflict still occurs after reaching the maximum number of planning rounds or the maximum number of particle batches, the current irradiation plan template is recorded as the planning conflict template, and the process moves to the next irradiation plan template; when all irradiation plan templates result in a conflict, the planning conflict results corresponding to the cross-boundary voxels and equipment constraints are output, and the planning process is then terminated.

[0069] The single-pass dwell time resolution of the device controller is 0.02 seconds, corresponding to a device realization resolution of 0.10 for continuous delivery weights. For a two-pass irradiation template, each device realization delivery weight is decomposed into two single-pass dwell times within the range of 0.10 seconds to 0.32 seconds and meeting the 0.02-second resolution. The sum of the two single-pass dwell times equals the product of 0.20 seconds and the device realization delivery weight. Device realization delivery weights that cannot complete the above decomposition are discarded. For each continuous delivery weight, the next adjacent realized value and the next adjacent realized value are obtained respectively; if two adjacent realized values ​​are the same, only one value is retained. A search tree is built according to the continuous delivery weight order to select adjacent realized values ​​layer by layer. The sum of the cumulative weighted absolute deviation of the selected values ​​and the minimum possible weighted absolute deviation of the unselected values ​​is used as the lower bound of the node. The node lower bound is decomposed as follows:

[0070]

[0071]

[0072] In the formula, Represents a search tree node The lower bound of the node; This represents the cumulative weighted absolute deviation of the deliverable weights for the selected equipment; This represents the minimum possible weighted absolute deviation of the delivery weight that can be achieved by the unselected equipment; and These represent the sets of selected and unselected candidate irradiation units, respectively; Indicates the equipment time conversion factor; This indicates the delivery weight that the selected equipment can achieve; and These represent the lower adjacent realizable value and the upper adjacent realizable value, respectively.

[0073] The search tree is expanded in ascending order of node lower bounds using a priority queue. For nodes with only partial assignments, the minimum and maximum reachable values ​​of the left-hand side of each device constraint are calculated using the lower and upper adjacent reachable values ​​of the unselected device's reachable delivery weight. Nodes are pruned when the allowable range of device power, conveying speed, or device synchronization constraints does not intersect with the corresponding reachable interval. Nodes are also pruned when both adjacent conveying path segments have been assigned values ​​and the delivery weight difference violates the constraint of changing delivery weights in adjacent conveying path segments. Leaf nodes that simultaneously satisfy all device constraints and are not in the verified combination table are retained as candidate combinations for device realization. The verified combination table uses the integer codes of each device control level to build a hash index. The preset maximum number of combinations is 32, and the preset maximum number of expanded nodes is 10,000. The search stops when there are 32 candidate combinations for device realization, the priority queue is cleared, or the number of expanded nodes reaches 10,000. This generates combinations to be verified in ascending order of deviation, without explicitly generating the Cartesian product of all adjacent values.

[0074] Figure 8 Taking the continuous delivery weights of two candidate irradiation units as an example, this illustrates the continuous irradiation planning parameters, adjacent feasible combinations, equipment constraint feasible region, and verification feedback constraints. The continuous irradiation planning parameters lie within the continuous feasible region, and equipment feasible candidate combinations can only be taken from discrete grid points; the minimum deviation combination search prioritizes neighboring grid points that satisfy the equipment constraints. When a dose crosses the limit during the complete process verification of the combination, the verification feedback constraint limits the continuous feasible region for the next round, ensuring that subsequent continuous irradiation planning parameters meet the planned mean dose limit corresponding to the dose crossover direction. This figure is only used to illustrate the two-dimensional relationship; the actual search is performed layer by layer according to all continuous delivery weights within the current irradiation plan template.

[0075] Before validation, a maximum total number of validations of 64 is set. This maximum number of validations covers all irradiation plan templates, all candidate combinations of equipment implementations, and all statistical validation decisions allowed in all feedback loops. Each calculation of the validation Monte Carlo statistical confidence dose bound for a candidate combination of equipment implementations increments the validation count by one. Each candidate combination of equipment implementations initially employs 40 complete process validation particle batches, each containing 2,000,000 primary particle histories. Each complete process validation particle batch uses a random number subsequence that does not overlap with the unit dose response calculations and does not overlap between different candidate combinations of equipment implementations or different validation counts, thus constituting an independent validation particle batch.

[0076] The complete process validation does not employ a re-linear combination of unit dose responses. Within each complete process validation particle batch, the complete irradiation process is directly simulated according to the radiation source output, delivery speed, single-pass dwell time, food carrier attitude, radiation source side order, food carrier flipping order, and number of irradiation passes corresponding to the candidate combinations realized by the equipment. For a 2-pass irradiation template, two irradiations are simulated sequentially according to the two single-pass dwell times corresponding to the deliverable weights of the equipment. Dose samples of food dose-constrained voxels are directly obtained from each complete process validation particle batch, and the voxel dose mean and voxel dose standard error are calculated. The validation population tail probability is set to 0.05, and its one-sided tail probability is equal to the validation population tail probability divided by the product of the number of food dose-constrained voxels, the maximum total number of validations, and 2. The validation Monte Carlo statistical confidence dose bound is:

[0077] In the formula, and These represent the lower bound and upper bound of the statistical confidence dose for validating Monte Carlo methods, respectively. This represents the average voxel dose of the particle batch during the complete process verification. Indicates the standard error of voxel dosage; This indicates the current cumulative number of particle batches verified throughout the complete process; This represents the cumulative probability obtained by subtracting the one-sided tail probability from 1; The degrees of freedom are the number of batches of particles verified in the complete process minus 1, and the cumulative probability is... The quantiles of the student's t-distribution. When the lower bound of the fully validated Monte Carlo statistical confidence dose reaches the lower limit of the target dose and the upper bound of the fully validated Monte Carlo statistical confidence dose is within the upper limit of the target dose, the output device can execute the irradiation planning parameters.

[0078] Verified device implementation candidate combinations are added to the verified combination table. When all candidate combinations in this round of device implementation result in cross-boundary outcomes, the maximum normalized dose cross-boundary amount for each combination is calculated:

[0079] In the formula, Indicates the device achieves candidate combinations Maximum normalized dose crossing; and These represent the lower limit of the target dose and the upper limit of the target dose, respectively. and These represent the lower and upper bounds of the validation Monte Carlo statistical confidence dose for the corresponding device implementation candidate combinations on the corresponding food dose-constrained voxels, respectively. The two inner-layer maximum value operations set the negative values ​​of the cold-point and hot-point crossover values ​​to 0, respectively. The device implementation candidate combination with the smallest maximum normalized dose crossover value is selected as the feedback target combination; subsequent feedback only uses the crossover results of the feedback target combination.

[0080] When the same food dose constraint voxel simultaneously exhibits a lower bound of the validation Monte Carlo statistical confidence dose below the target dose lower limit and an upper bound of the validation Monte Carlo statistical confidence dose above the target dose upper limit, it is identified as a statistical width cross-boundary voxel. A statistical width cross-boundary voxel indicates that the statistical interval width is greater than the dose window, and it cannot simultaneously form directional constraints for increasing and decreasing the dose. In this case, 20 complete process validation particle batches are added to the feedback target combination, with the cumulative number of complete process validation particle batches being 40, 60, 80, and 100 respectively. The validation Monte Carlo statistical confidence dose boundary is recalculated, and the validation count is increased. When the cumulative number of complete process validation particle batches reaches 100 and statistical width cross-boundary voxels still exist, the current irradiation plan template is recorded as a planning conflict template and transferred to the next irradiation plan template. When the number of statistical width cross-boundary voxels is 0, cold spot cross-boundary voxels and hot spot cross-boundary voxels respectively form the following validation feedback constraints:

[0081]

[0082] In the formula, Indicates the weight of the next round of continuous delivery; This represents the average planned dose of the feedback target combination on the corresponding food dose-constrained voxel; This represents the positive feedback dose step size, which is 0.05 kGray in this implementation. The current round of validation feedback constraints replaces the previous round's validation feedback constraints in the current irradiation plan template, without accumulating constraints from different failed combinations. This avoids infeasibility caused by the simultaneous existence of opposite-direction validation feedback constraints for candidate combinations implemented by different devices. When the validation count reaches 64, validation of all candidate combinations implemented by all devices is stopped. The current irradiation plan template and the irradiation plan templates that have not yet been validated are recorded as planning conflict templates, and the planning conflict results corresponding to cross-boundary voxels and device constraints are output.

[0083] In one calculation of this embodiment, the minimum lower bound of the planned Monte Carlo statistical confidence dose corresponding to the continuous irradiation planning parameters is 3.12 kGy, and the maximum upper bound of the planned Monte Carlo statistical confidence dose is 6.82 kGy. The first complete process verification of the minimum deviation device implementation candidate combination yielded a verification Monte Carlo statistical confidence dose lower bound of 2.96 kGy at a food dose constraint voxel near the center of the food carrier. After selecting the device implementation candidate combination with the smallest maximum normalized dose crossover as the feedback target combination and forming the verification feedback constraint, the minimum lower bound of the verification Monte Carlo statistical confidence dose for the next round of device implementation candidate combinations is 3.08 kGy, and the maximum upper bound of the verification Monte Carlo statistical confidence dose is 6.91 kGy, satisfying the target dose window. The output includes the source output level, the dwell time per pass for each transport path segment within each irradiation pass, the radiation source side order, the food carrier flipping order, and the number of irradiation passes. The positions of the minimum and maximum voxels of the verification dose are converted to the food body coordinate system and used as the predicted cold spot and predicted hot spot positions, respectively. The predicted cold spot and predicted hot spot positions are listed together as the recommended placement locations for the dosimeter.

[0084] To verify the effect of projection-preserving quality closure correction, global scaling was used as a control. Under the condition that the total mass of food candidate voxels reached 8.00 kg in both methods, the weighted projection residual norm ratio after global scaling was 0.054, while the ratio after projection-preserving quality closure correction was 0.019. An average number of additional particle batches was used as a control: adding 5 particle batches to each of the 48 candidate irradiation units in the current template required 240 additional particle batches; after selecting based on the dose variance contribution of active constraint voxels, adding a total of 55 particle batches to 11 candidate irradiation units was sufficient to reach the statistical error threshold. Compared with global scaling, projection-preserving quality closure correction reduced the weighted projection residual norm ratio; compared with average additional particle batches, selecting the additional targets based on dose variance contribution reduced the number of new particle batches required to reach the statistical error threshold.

[0085] As an optional implementation, the irradiation equipment can be a 10 MeV electron beam device, with the food carrier thickness set to 50 mm, and the reference delivery amount using a 0.10 microcoulomb reference beam charge. The number of primary electrons corresponding to the reference beam charge is obtained by dividing the reference beam charge by the elementary charge; the candidate irradiation unit still consists of the radiation source side, food carrier attitude, delivery path segment, source term energy spectrum, source output level, and reference delivery amount. When the food carrier is irradiated from the top and bottom sides and subjected to 180° rotation, the depth dose drop caused by the electron range can be compensated. Particle batch division, food body coordinate transformation, unit dose response statistics, confidence dose constraint inverse programming, and equipment implementation candidate combination verification remain unchanged.

[0086] As another optional implementation, the irradiation equipment can be a cobalt-60 gamma-ray device, and the baseline delivery quantity uses a 1-second baseline irradiation time. When there is only one source output position, the fixed source output state is used as the only source output position. The irradiation plan template still consists of the radiation source side order, the food carrier flipping order, and the number of irradiation passes. The continuous delivery weight corresponds to the dwell time or conveying speed. For food carriers with large internal cavities or mixed food packaging, the size of the food candidate voxels can be reduced, and the food element mass fraction and density boundary corresponding to the food category can be set separately. When all food candidate voxels with an equivalent mass density greater than 0 are determined with the same density to determine the transport material mass and voxel food mass, the net food mass, particle transport material mass, and voxel food mass maintain mass consistency.

[0087] The predicted cold spot locations, predicted hot spot locations, and recommended dosimeter placement locations output by this method are used for irradiation process design and dosimeter placement, but do not replace irradiation equipment performance verification, physical dose distribution measurement, or routine process monitoring. When changing the food category, packaging structure, source term energy spectrum, or food carrier orientation, the corresponding mass constraint voxel model and unit dose response matrix should be re-established to maintain the linear superposition conditions of the planning stage.

Claims

1. A Monte Carlo particle transport simulation and constrained inverse programming method for food irradiation dose fields, characterized in that, include: The process involves acquiring the net mass of the food, candidate regions of the food, and angle-limited transmission scanning data. The candidate regions are discretized into candidate voxels. Based on the angle-limited transmission scanning data, the equivalent mass density of each candidate voxel is updated ray-by-ray. After each ray-by-ray update, a projection-preserving mass closure correction is performed to obtain the equivalent mass density distribution. Candidate voxels with equivalent mass densities greater than 0 and less than the solid density of the food are constructed as de-density equivalent voxels. The product of the equivalent mass density and the voxel volume is used as the voxel's food mass, resulting in a mass-constrained voxel model. For candidate irradiation units composed of radiation source side, food carrier attitude, transport path segment, source term energy spectrum, source output level, and baseline delivery quantity, particle emission, motion, collision, scattering, and energy deposition are simulated to obtain the unit dose response mean and unit dose response standard error. Enumerate the irradiation plan templates for discrete conditions of fixed equipment, perform confidence dose constraint backward programming with continuous delivery weight as variable in each irradiation plan template, and add particle batches to the corresponding candidate irradiation units according to the dose variance contribution of active constraint voxels. The continuous irradiation planning parameters are converted into a set of candidate combinations for equipment implementation. The complete irradiation process of each candidate combination for equipment implementation is simulated using independent verification particle batches. Based on the verification results, the executable irradiation planning parameters of the equipment are output, or the confidence dose constraint inverse planning is re-executed after forming verification feedback constraints.

2. The method according to claim 1, characterized in that, Acquiring angle-limited transmission scanning data includes: obtaining the net weight of the food based on the total cargo mass, the tare weight of the transport vehicle, and the registered weight of the packaging; using a depth camera to collect point clouds of the outer surface of the food loading body and determine the packaging area, the food candidate area, and the external air area; acquiring transmission intensity images, dark field intensity images, and air field intensity images at -50°, -25°, 0°, 25°, and 50°, and obtaining the initial attenuation value after dark field correction and air field correction; measuring food equivalent material calibration blocks of multiple thickness levels using a beam blocking array, and establishing a scattering distribution correction table according to the scanning angle, detection channel, and the attenuation interval where the initial attenuation value is located; placing the food equivalent material calibration blocks in the same packaging and transport vehicle background as the food loading body, establishing a background matching calibration curve between the scattering-corrected transmission measurement value and the newly added food equivalent mass thickness according to the scanning angle and detection channel, and obtaining the measured food equivalent mass thickness using the scattering distribution correction table and the background matching calibration curve.

3. The method according to claim 1, characterized in that, Ray-by-ray updates include: determining the mass capacity range based on the volume, lower density limit, and upper density limit of food candidate voxels; when the net mass of the food is within the mass capacity range, dividing the net mass of the food by the total volume of the food candidate voxels sets the initial equivalent mass density of the food; when the net mass of the food is outside the mass capacity range, outputting modeling condition conflict information; establishing a projection matrix based on the ray's travel length within each food candidate voxel; for each ray, calculating the projection residual between the measured equivalent mass thickness of the food and the sum of the products of each travel length and the corresponding current equivalent mass density of the food; for each food candidate voxel traversed by the current ray, dividing the product of the projection residual, the corresponding travel length, and a relaxation coefficient greater than 0 and less than 2 by the sum of the squares of the current ray's corresponding travel length and the dimensionless positive stability quantity, and superimposing the resulting density increment onto the current equivalent mass density of the food, limiting the superposition result to between the corresponding lower and upper density limits.

4. The method according to claim 1, characterized in that, The projection-preserving quality closure correction includes: weighting the projection changes of each ray according to the reciprocal of the standard deviation of the corresponding ray's repeated calibration; using the weighted sum of the squares of the weighted projection changes caused by the voxel density correction and the squares of the voxel 6-neighborhood differences of the corrected food equivalent mass density as the objective function; and using the constraints that the total corrected food candidate voxel mass equals the net food mass and that the corrected food equivalent mass density is between the corresponding lower and upper density limits. The primal-dual interior-point method is used to solve the convex quadratic programming problem. Ray-by-ray updates and projection-preserving quality closure corrections are performed alternately until the weighted projection residuals are within the range... The ratio of the number to the weighted measured food equivalent mass thickness norm is less than or equal to the preset projection residual threshold, and the relative difference between the total mass of food candidate voxels and the net mass of food is less than or equal to the preset mass residual threshold; when the preset maximum number of reconstruction rounds is reached and any ratio exceeds the corresponding threshold, modeling condition conflict information is output; food candidate voxels with an equivalent mass density of 0 are assigned an air material label, and food candidate voxels with an equivalent mass density greater than 0 are established as particle transport materials according to the corresponding food equivalent mass density and food element mass fraction, and all deposition energy within the voxels is included in the food deposition energy.

5. The method according to claim 1, characterized in that, The electron beam equipment uses a reference beam current charge as the reference delivery quantity, while the continuous X-ray and gamma-ray equipment use a reference irradiation time as the reference delivery quantity. Candidate voxels whose voxel quality reaches a preset dosimetry quality threshold are grouped into a food dose-constrained voxel set. Each candidate irradiation unit generates at least 30 batches of particles with equal historical data using mutually separate random number subsequences. Based on the rigid body transformation relationship of the food carrier attitude, the energy deposition position is inversely transformed to the food's body coordinate system. The food deposition energy is divided by the voxel food quality and the reference delivery quantity to obtain the reference delivery quantity dose response, and the unit dose response mean and unit dose response standard error are calculated. During the continuous delivery weighting solution, the material state, geometric pose, and source term energy spectrum of the candidate irradiation unit are kept fixed. When adding particle batches, subsequent random number subsequences are called, and the unit dose response mean and unit dose response standard error of the corresponding candidate irradiation unit on all food dose-constrained voxels are updated.

6. The method according to claim 5, characterized in that, The confidence dose-constrained inverse programming includes: enumerating a finite number of irradiation plan templates based on the radiation source side separation order, food carrier flipping order, number of irradiation passes, and source output levels; setting the range of continuous delivery weights based on the delivery speed range or residence time range, where the continuous delivery weight represents the ratio of the total delivery amount of each irradiation pass on the same candidate irradiation unit to the baseline delivery amount; before the start of each planning round, determining the conservative confidence coefficient of the Student's t-distribution based on the minimum degree of freedom obtained by subtracting one from the number of particle batches in each candidate irradiation unit in the current irradiation plan template, the planning tail parameters, the number of food dose-constrained voxels, and the maximum number of planning rounds; and determining the sum of the products of the unit dose response mean and the continuous delivery weights. The planned dose mean is determined by the square root of the sum of squares of the products of the unit dose response standard error and the continuous delivery weights. The upper bound of the planned Monte Carlo statistical confidence dose is determined by the sum of the planned dose mean, the conservative confidence coefficient of the Student's t-distribution, and the planned dose standard error. The subtraction result is determined as the lower bound of the planned Monte Carlo statistical confidence dose. The food dose constraint is defined as the lower bound of the planned Monte Carlo statistical confidence dose reaching the lower limit of the target dose and the upper bound of the planned Monte Carlo statistical confidence dose being within the upper limit of the target dose. Constraints of equipment power, conveying speed, and changes in delivery weights of adjacent conveying path segments are added. The primal-dual interior point method is used to solve the continuous second-order cone programming problem.

7. The method according to claim 6, characterized in that, The planning objective of continuous second-order cone programming is the sum of the products of continuous delivery weight and equipment time conversion factor or equipment energy consumption conversion factor, and the weighted sum of the absolute differences of continuous delivery weights of adjacent transport path segments; when food dose constraints are met, food dose constraint voxels whose non-negative difference between the planned Monte Carlo statistical confidence dose lower bound and the target dose lower bound or the target dose upper bound and the planned Monte Carlo statistical confidence dose upper bound is within a preset judgment margin are determined as active constraint voxels; The square of the product of the unit dose response standard error and the continuous delivery weight is determined as the dose variance contribution. The candidate irradiation units with the fewest cumulative contributions reaching the preset variance ratio are selected in descending order of dose variance contribution. Particle batches are added to the union of the candidate irradiation units selected for each active constraint voxel. The mean of the unit dose response, the standard error of the unit dose response, and the conservative confidence coefficient of the Student t-distribution are updated, and the continuous second-order cone programming is solved again. When the continuous second-order cone programming forms a conflict state, a unified dose relaxation is added to obtain the minimum relaxation solution, and particle batches are added based on the cross-boundary voxels of the minimum relaxation solution. When the maximum number of planning rounds or the maximum number of particle batches are reached, the current irradiation plan template is recorded as the planning conflict template.

8. The method according to claim 1, characterized in that, The generation of the candidate combination set for device implementation includes: obtaining the lower and upper adjacent realizable values ​​of each continuous delivery weight based on the device control resolution, and establishing a search tree for selecting adjacent realizable values ​​layer by layer; determining the lower bound of the node by the sum of the cumulative weighted absolute deviation between the selected adjacent realizable values ​​and the corresponding continuous delivery weights, and the weighted absolute deviations of the remaining continuous delivery weights formed by using adjacent realizable values ​​with smaller deviations; establishing a verified combination table to store the completed verified combinations, expanding the search tree nodes in ascending order according to the lower bound of the nodes using a priority queue, and retaining the device implementation candidate combinations that meet the constraints of device power, delivery speed, changes in delivery weights of adjacent delivery path segments, and device synchronization at the leaf nodes, and that are outside the verified combination table, until the number of device implementation candidate combinations reaches the preset maximum number of combinations, the priority queue is cleared, or the number of expanded nodes reaches the preset maximum number of expanded nodes.

9. The method according to claim 8, characterized in that, Before verification, a maximum total number of verifications is set, covering all irradiation plan templates, all candidate combinations of equipment implementation, and all feedback loops, allowing for statistical verification decisions. For the first verification of each candidate combination of equipment implementation, a complete process verification particle batch of at least 30 particles of the same quantity is used. The complete irradiation process is simulated according to the corresponding radiation source output, delivery speed, residence time, food carrier attitude, radiation source side separation order, food carrier flipping order, and number of irradiation passes. Each candidate combination of equipment implementation and each verification uses a mutually separate random number subsequence. The overall verification is... The tail probability is divided by the product of the number of voxels constrained by the food dose, the maximum total number of validations, and 2 to obtain the one-sided tail probability. The Student's t-distribution quantile is determined by subtracting the one-sided tail probability from the cumulative probability from 1 and subtracting 1 from the number of validation particle batches for the complete process. The upper bound of the validation Monte Carlo statistical confidence dose is obtained by adding the voxel dose mean to the product of the Student's t-distribution quantile and the standard error of the voxel dose. The lower bound of the validation Monte Carlo statistical confidence dose is obtained by subtracting the corresponding product from the voxel dose mean. Each calculation of the validation Monte Carlo statistical confidence dose bound increases the validation count by 1. When all food dose constraints are met, the output device can execute irradiation planning parameters, predict cold spot locations, predict hot spot locations, and recommend dosimeter placement locations.

10. The method according to claim 9, characterized in that, Write the verified candidate combinations of equipment implementations into the verified combination table; when all candidate combinations of equipment implementations in this round form cross-boundary results, calculate the ratio of the positive difference between the lower limit of the target dose and the lower limit of the verified Monte Carlo statistical confidence dose to the lower limit of the target dose, and the ratio of the positive difference between the upper limit of the verified Monte Carlo statistical confidence dose and the upper limit of the target dose to the upper limit of the target dose. Determine the maximum value of all ratios as the maximum normalized dose cross-boundary amount, and select the candidate combination of equipment implementations with the smallest maximum normalized dose cross-boundary amount as the feedback target combination; Food dose constraint voxels in the same food dose constraint combination that simultaneously exhibit both cold-point and hot-point crossovers are identified as statistically wide crossover voxels. A preset number of complete process verification particle batches are added, and the complete irradiation process verification is re-executed using the accumulated complete process verification particle batches. When the number of statistically wide crossover voxels is 0, verification feedback constraints are formed for cold-point crossover voxels, ensuring that the average planned dose of subsequent rounds reaches the sum of the average planned dose of the feedback target combination and the positive feedback dose step size. For hot-point crossover voxels, verification feedback constraints are formed, ensuring that the average planned dose of subsequent rounds is within the difference between the average planned dose of the feedback target combination and the positive feedback dose step size. The verification feedback constraints of this round replace the verification feedback constraints of the previous round, and the confidence dose constraint inverse planning is re-executed. When the verification count reaches the maximum total number of verifications or all irradiation plan templates are planning conflict templates, all equipment candidate combination verification is stopped. Irradiation plan templates that have not yet completed verification are recorded as planning conflict templates, and the planning conflict results corresponding to the crossover voxels and equipment constraints are output.

Images