Proton radiotherapy dose and biological effect calculation method and system based on dual-energy CT

Abstract

The application discloses a kind of based on dual-energy CT's proton radiotherapy dose and biological effect calculation method and system, belong to medical physics, radiotherapy and artificial intelligence cross technical field.The application obtains relative electron density and effective atomic number by high, low X-ray energy spectrum CT image, obtains effective RSP by deterministic differentiable physical layer;Treatment plan is split into beam unit, and obtains each beam unit gradually beam physical dose, beam range characteristic and residual energy representation, and according to the beam propagation direction, the interface upstream and downstream medium is ordered sampling, and constitutes direction related interface variable;From the variable, the gate strength that attenuates with the distance increase from voxel to interface is determined, and the modulation amplitude is bounded and can be positive or negative gradually beam LET d Residual error;Each beam physical dose is added to obtain total dose, and each beam LET d Total LET is obtained according to physical dose weighting d Mask is executed to low-dose voxel, and RBE weighted dose is calculated according to each fraction physical dose.

Description

Technical Field

[0001] This invention belongs to the interdisciplinary field of medical physics, radiotherapy and artificial intelligence, and particularly relates to a method and system for calculating proton radiotherapy dose and biological effects based on dual-energy CT. Background Technology

[0002] The energy deposition and range of a proton beam are influenced by both the electron density and elemental composition of the medium. The conversion from monoenergetic CT value to relative stopping power (RSP) is not unique; abrupt changes in the medium surrounding the lung-chest wall, bone-soft tissue, and implant materials simultaneously alter the energy loss rate, multiple Coulomb scattering, and secondary particle contributions from nuclear reactions, thus affecting the physical dose and dose-average linear energy transfer (LET) at the end of the range and in the interface neighborhood. d An accurate description is even more difficult.

[0003] In existing technologies, dual-energy CT learning mapping to RSP, physical constraints based on Bethe-Bloch relations, beam-by-beam prediction of physical dose in the beam coordinate system, and prediction of LET from physical dose distribution are all employed. d And the LET uses a variable RBE (Relative Biological Effectiveness) parameterized model. d Both physical dosimetry and biological dosimetry have corresponding protocols in place. Furthermore, LET... d According to its definition of dose averaging, multiple beams should be combined using dose-weighted summation, and LET should not be reported in the low-dose region. d To avoid numerical instability, this is a standard practice in the field. Therefore, the aforementioned individual techniques, or their simple combination, are insufficient to define continuous transport relationships in non-uniform media.

[0004] In situations involving multiple illumination fields, multiple energy layers, or overlapping contributions from multiple pencil beams, the same physical interface exhibits different upstream and downstream relationships and transport states for different beam units. When a proton enters a medium with high relative electron density or high effective atomic number from a medium with low relative electron density or low effective atomic number, the scattering conditions, spectral distortion, and secondary particle equilibrium states in the interface neighborhood differ compared to when it traverses the same interface in the reverse order. d The deviation from the extrapolated value relative to a homogeneous medium may differ in both sign and magnitude. Existing LET... d Predictive models typically use CT images and dose distribution as input, but their inputs do not include the beam propagation direction or distinguish the order of materials on both sides of the interface, thus failing to express the aforementioned directional asymmetry; further methods can only enable LET. dA unidirectional increase in multiplicative correction will result in a systematic overestimation in areas where a negative bias should occur, and this bias cannot be canceled out in the interface region of multi-field cross-irradiation.

[0005] Furthermore, if the interface correction amount is not constrained in space, its effect may extend to a uniform medium region far from the interface, which contradicts the physical intent of interface modulation; if the correction amount is not constrained in amplitude, it is difficult to limit obviously non-physical outputs. Summary of the Invention

[0006] The technical problem to be solved by this invention is: without using conventional network backbone, Monte Carlo training data, beam-by-beam dose prediction itself, or LET. d Based on the dose-weighted definition as a technical contribution, a physical-biological co-computation method is established that can transmit information on medium state, residual energy, range position, and interface orientation on a beam-by-beam scale, enabling LET in the interface neighborhood. d The correction possesses directional discrimination capability and bounded locality, and correctly generates total physical dose and total LET under multi-beam conditions. d Weighted dose with RBE.

[0007] The first aspect discloses a method for calculating proton radiotherapy dose and biological effects based on dual-energy CT, the method comprising:

[0008] The relative electron density and effective atomic number of each voxel are calculated based on the high-energy X-ray spectral CT images and low-energy X-ray spectral CT images of the object to be calculated, and the effective RSP is obtained through a deterministic differentiable physical layer.

[0009] The treatment plan is divided into multiple beam units, and the physical dose and residual energy characterization of each beam unit are obtained. The upstream and downstream media of the interface are sampled in an orderly manner according to the beam propagation direction to form direction-related interface variables. The beam unit is a set of beams with a common incident direction and known initial energy distribution. The residual energy characterization is determined by the effective RSP path integral and the pre-calibrated range-energy relationship.

[0010] The gating strength, which decreases with increasing voxel-to-interface distance, is determined by the direction-dependent interface variables, and the LET beam-by-beam modulation is performed using the gating strength. d The residuals are used to obtain LET bundles. d,b ;

[0011] The total physical dose is obtained by adding up the physical doses of each bundle, and the LET is calculated by adding up the physical doses of each bundle. d,b Total LET was obtained by weighting based on physical dose. d Within a voxel range where the total physical dose is not lower than a preset dose threshold, the RBE-weighted dose is calculated based on the physical dose per fraction, and the effective RSP, total physical dose, and total LET are output. dAnd RBE weighted dose.

[0012] The second aspect discloses a proton radiotherapy dose and biological effect calculation system based on dual-energy CT, the system comprising:

[0013] The dual-energy spectral physics parameter calculation module is used to calculate the relative electron density and effective atomic number of each voxel based on the high-energy spectral CT images and low-energy spectral CT images of the object to be calculated, and obtain the effective RSP through a deterministic differentiable physical layer.

[0014] The module for calculating beam-by-beam physical dose and direction-dependent interface variables is used to break down the treatment plan into multiple beam units, obtain beam-by-beam physical dose and residual energy characterization for each beam unit, and sample the upstream and downstream media of the interface in an orderly manner according to the beam propagation direction to form direction-dependent interface variables. The beam unit is a set of beams with a common incident direction and a known initial energy distribution. The residual energy characterization is determined by the effective RSP path integral and the pre-calibrated range-energy relationship.

[0015] LET by bundle d,b The calculation module is used to determine the gating intensity, which decreases with increasing voxel-to-interface distance, from the direction-related interface variables, and to modulate the beam-by-beam LET using the gating intensity. d The residuals are used to obtain LET bundles. d,b ;

[0016] The multi-beam aggregation and bio-dosage calculation module is used to add up the physical doses of each beam to obtain the total physical dose, and to calculate the LET of each beam. d,b Total LET was obtained by weighting based on physical dose. d Within a voxel range where the total physical dose is not lower than a preset dose threshold, the RBE-weighted dose is calculated based on the physical dose per fraction, and the effective RSP, total physical dose, and total LET are output. d And RBE weighted dose.

[0017] As can be seen from the above technical solutions, the present invention has the following beneficial effects:

[0018] This invention obtains relative electron density and effective atomic number from high- and low-energy X-ray spectral CT images, and obtains effective RSP through deterministic differentiable physical layers. The treatment plan is broken down into beam units; for each beam unit, beam-by-beam physical dose, beam-by-beam range characteristics, and residual energy characterization are obtained. The upstream and downstream media of the interface are sampled in an orderly manner according to the beam propagation direction, forming a direction-dependent interface variable. This variable determines the gating intensity that decreases with increasing voxel-to-interface distance, resulting in a bounded and potentially positive / negative beam-by-beam LET. d Residual. The total dose is obtained by adding the physical doses of each bundle, and the LET of each bundle. d Total LET was obtained by weighting based on physical dose.d Masking was performed on low-dose voxels, and the RBE-weighted dose was calculated based on the physical dose per fraction.

[0019] This invention characterizes the residual energy determined by the effective RSP and beam parameters, the beam-by-beam physical dose and range characteristics, the direction-dependent and ordered interface variables, and the LET with bounded amplitude and a sign that can be positive or negative, and is limited to acting on the interface neighborhood. d Residual corrections are written into the same continuous calculation chain indexed by the beam element, and the correction results of each beam are correctly aggregated under their own dose weights. Attached Figure Description

[0020] Figure 1 The present invention provides an overall flowchart of a method for calculating proton radiotherapy dose and biological effects based on dual-energy CT.

[0021] Figure 2 This invention provides a schematic diagram of a proton radiotherapy dose and biological effect calculation system based on dual-energy CT. Detailed Implementation

[0022] To make the objectives, features, and advantages of the present invention more apparent and understandable, specific embodiments of the present invention will be described in detail below with reference to the accompanying drawings. Several embodiments of the present invention are shown in the drawings. However, the present invention can be implemented in many different forms and is not limited to the embodiments described herein. Rather, these embodiments are provided to make the disclosure of the present invention more thorough and complete. Obviously, the described embodiments are only some, not all, of the embodiments of the present invention. All other embodiments obtained by those skilled in the art based on the embodiments of the present invention without inventive effort are within the scope of protection of the present invention.

[0023] The terms "first," "second," etc. (if present) in the specification, claims, and accompanying drawings of this invention are used to distinguish similar objects and are not necessarily used to describe a particular order or sequence. It should be understood that such data used in this way can be interchanged where appropriate so that the embodiments described herein can be implemented in orders other than those illustrated or described herein. Furthermore, the terms "comprising" and "having," and any variations thereof, are intended to cover a non-exclusive inclusion; for example, a process, method, system, product, or apparatus that comprises a series of steps or units is not necessarily limited to those steps or units explicitly listed, but may include other steps or units not explicitly listed or inherent to such processes, methods, products, or apparatus.

[0024] Unless otherwise defined, all technical and scientific terms used in this invention have the same meaning as commonly understood by one of ordinary skill in the art to which this application belongs; the terminology used in the specification of this application is for the purpose of describing particular embodiments only and is not intended to be limiting of this application. The reference to "embodiment" in this invention means that a particular feature, structure, or characteristic described in connection with an embodiment may be included in at least one embodiment of this application. The appearance of this phrase in various places in the specification does not necessarily refer to the same embodiment, nor is it a separate or alternative embodiment mutually exclusive with other embodiments. It will be explicitly and implicitly understood by those skilled in the art that the embodiments described herein can be combined with other embodiments.

[0025] The following embodiments are used to illustrate the minimum operable computational chain of the present invention and do not constitute a limitation on specific scanning devices, network backbones, Monte Carlo platforms, or variable RBE parameterization models. During implementation, the physical quantity definitions, spatial grids, and units of the training data, reference labels, and application data should be kept consistent. Parameters and function forms that have not been experimentally calibrated are given as optional embodiments.

[0026] To facilitate understanding of the embodiments of the present invention, the basic terms and physical quantity definitions will be explained below.

[0027] Effective atomic number Z eff The calculation, based on the power-law mixed rule, is shown in the following formula:

[0028] ;

[0029] in, For the first The electron share of an element For the first The atomic number of the element, the exponent m is selected before training, and can be 3.1 or 3.3. Once selected, this exponent must be within (Z). eff I) Calibration data construction, network Z eff The generation of supervisory labels, interpolation lookup tables, and application stages should be kept consistent; the same material under different indices should be monitored. Using different values ​​interchangeably will lead to systematic biases in the table lookup results.

[0030] Dose-mean linear energy transfer (LET) d Scoring is performed within voxels using the following formula:

[0031] ;

[0032] Where j represents the number of charged particle steps involved in the counting within the voxel. For linear energy transfer in this step, Let this be the step size. This step involves energy deposition within the voxel, and satisfies the formula... The scoring definition includes the following elements, which must be consistent throughout the scheme and recorded with the model: (a) the types of charged particles involved in the counting, such as only primary and secondary protons, or all charged particles including alpha particles and recoil nuclei; (b) the medium reference on which LET is based, i.e., water or a local medium; (c) whether track truncation is used, i.e., using a restricted LET with a cutoff energy of Δ (correspondingly excluding the energy carried by secondary electrons with energies higher than Δ), or using a non-restricted LET without truncation; (d) whether dose averaging is performed after averaging the complete track of each particle, or whether dose averaging is performed directly on all steps within the voxel. This invention's LET d In keV / μm units, training labels, network output, and upper limit of residual amplitude A max This unit is used for both inputs to the RBE model and inputs.

[0033] The effective RSP used in the treatment plan refers to the energy E within the calibrated energy range, based on a preset representative energy. ref Calculated RSP(E) ref ), or the RSP diagram obtained by integrating or discretizing the RSP(E) according to preset energy weights. E ref The preset energy weights are kept consistent throughout the training, label generation, and application phases of the same model. Specifically, if the beam unit is a single energy layer, the representative energy is the nominal initial energy of that energy layer, or the energy calculated from the pre-calibrated range-energy relationship according to the preset path depth. If the beam unit contains multiple energy layers, the particle flux, monitoring unit, or planned relative weight of each energy layer is normalized and used as the energy weight, and the RSP at the corresponding energy is weighted layer by layer. When using a single representative energy to approximate multiple energy layers, the energy that minimizes the deviation between the RSP on the calibration material set and the layer-by-layer weighting result is selected from the candidate energies of the calibration energy region before training. The representative energy, energy weights, and calculation order are saved together with the model and remain consistent throughout the label generation, training, and application phases.

[0034] A beam unit refers to a collection of beams that share a common incident direction and have a single initial energy value or a known initial energy distribution. A preferred embodiment uses energy layers as beam units. When using an irradiation field as a beam unit, because the residual energy of different energy layers within the same irradiation field varies significantly, it is necessary to calculate the residual energy characterization for each energy layer separately and then weight it according to its relative weight, or to calculate the LET (Lead Time) per beam for each energy layer. d Then, the energy is weighted by physical dose; otherwise, the remaining energy characterization is not unique. The highest accuracy is achieved when using pencil beams as beam units, but the number of beam units and the storage and computation costs of each beam tag increase accordingly. In implementation, the energy layer and pencil beams can be selected based on available resources, and beam units with very low weights can be merged.

[0035] For obtaining training labels, it is often difficult to directly obtain the voxel-by-voxel relative electron density ρ of real objects. e Effective atomic number Z eff Given the media's average ionization energy I and the true value of RSP, one feasible method for generating training data is as follows: Starting with a digital material phantom with known elemental composition and density, on the one hand, orthogonal projection is performed according to the X-ray energy spectrum corresponding to the selected high and low tube voltages and the corresponding elemental mass attenuation coefficients. After adding noise that conforms to the detector's statistical characteristics, reconstruction is performed to obtain the HU corresponding to each voxel of the material map. H With HU L Image; on the other hand, ρ is calculated from the same material diagram according to its elemental composition. e Calculate Z according to the power law rule eff The relative electron density (REP) label is calculated according to the Bragg additive rule, and the RSP label is calculated according to the RSP formula in this paper. Based on this, beam hardening, scattering, and metal artifact models can be introduced to expand the training distribution. The generated synthetic images and measured images can be mixed for training at a preset ratio and verified using measured tissue equivalent phantom data. The specific calculation process for the REP label is as follows: For each material voxel, the material mass density is multiplied by the weighted sum of the ratios of the mass fraction, atomic number, and relative atomic mass of each element to obtain a quantity proportional to the number of electrons per unit volume, which is then divided by the corresponding amount of water obtained using the same elemental data and calculation rules; Avogadro's constant is canceled out in this ratio. For tissue equivalent phantom materials with material certificates, calibration reports, or product technical data, the relative electron density parameters given therein can be used; for materials without provided relative electron density parameters, the relative electron density is calculated according to the measured mass density and known elemental composition, following the rules described above. The parameter sources, measurement conditions, and calculation rules used are kept consistent throughout the training label generation and model verification process.

[0036] Material calibration tables use known elemental composition or known ρ e Z eff Establishment of Standard Organizational Equivalent Materials (Z) eff I) Calibration Table. Piecewise linear interpolation or conformal piecewise cubic interpolation is used for the calibration data, with a constraint that I is greater than zero; this interpolation function is fixed during the main network training. Z outside the calibration range eff Instead of direct extrapolation, it generates out-of-range hints. For metal implants, independent material branches or extended calibration tables can be established, but calibration relationships of ordinary tissues must not be unconditionally applied to uncalibrated metals.

[0037] In one embodiment, the present invention discloses a method for calculating proton radiotherapy dose and biological effects based on dual-energy CT, such as... Figure 1 As shown, the specific steps include:

[0038] S101. Calculate the relative electron density and effective atomic number of each voxel based on the high X-ray energy spectrum CT image and low X-ray energy spectrum CT image of the object to be calculated, and obtain the effective RSP through a deterministic differentiable physical layer.

[0039] One implementation uses spatially registered 80 kVp low-spectral X-ray CT images and 140 kVp high-spectral X-ray CT images. Coordinate unification, voxel resampling, spatial registration, noise suppression, and CT value consistency calibration are performed on the high- and low-spectral X-ray CT images. Both sets of images are resampled to the same object coordinate system, for example, using a 2 mm × 2 mm × 2 mm voxel grid. The training labels and beam feature maps use the same voxel grid. The input tensor contains at least the high-spectral CT value HU. H Low-energy spectral CT value HU L The difference or ratio between the two can be supplemented with neighborhood statistics, metal artifact masking, or anatomical region coding. Tube voltage and voxel size are executable examples; other scanning and resampling methods capable of forming two sets of effective energy spectra are not limited. Once selected, tube voltage combinations, filtering conditions, and reconstruction kernels should be kept consistent during training and application, or normalized through CT value consistency calibration. Artifact masking should be preserved for areas affected by metallic implants.

[0040] S102. The treatment plan is divided into multiple beam units, and the physical dose and residual energy characterization of each beam unit are obtained. The upstream and downstream media of the interface are sampled in an orderly manner according to the beam propagation direction to form direction-related interface variables. The beam unit is a set of beams with a common incident direction and a known initial energy distribution. The residual energy characterization is determined by the effective RSP path integral and the pre-calibrated range-energy relationship.

[0041] Specifically, within the calibrated proton energy range, the deterministic differentiable physical layer, based on the electron stopping process generated by the interaction between protons and dielectric electrons, calculates the proton relative stopping power ratio related to the effective energy using the following formula. , i.e., effective RSP:

[0042] ;

[0043] In the formula, Represents voxels, Represents the relative electron density relative to water; with I w Representing voxels The average ionization energy of the corresponding medium and water, where Equivalent materials from standard organizations (Z) eff I) The positive value constraint piecewise interpolation relationship constructed from the calibration data is determined; γ is the ratio of the proton's velocity to the speed of light. p This represents the Lorentz factor corresponding to the proton velocity. With γ p All are determined by the proton energy E; m e c 2 This represents the electron rest energy. This formula is a first-order approximation of the electron stopping power ratio applicable within the proton energy range. The density effect correction term δ for the medium and water has a relatively small impact on the ratio within the therapeutic energy range and is therefore cancelled. For higher precision, density effect corrections, shell corrections, and higher-order z-corrections (Barkas-Andersen and Bloch terms) can be introduced at a consistent level for the medium and water. Note that the shell correction does not completely cancel out between high atomic number materials and water, and its effect increases as the proton energy decreases. For high atomic number materials such as metals, independent material calibration relationships should be established, rather than relying solely on the extrapolation of the above formula. The material calibration relationships and physical parameters of the deterministic differentiable physics layer remain fixed during training, but their output depends on ρ. e With Z eff Differentiable; therefore, the gradient of the RSP physical consistency loss can be backpropagated to the relative electron density prediction branch, the effective atomic number prediction branch, and the shared feature encoder according to the chain rule, and used to update the trainable parameters of the above network based on the deviation between the predicted RSP and the reference RSP, so that the outputs of the two physical parameter branches are jointly constrained under the same deterministic RSP relationship, without updating the physical layer itself.

[0044] Furthermore, regarding the treatment of nuclear reaction contributions: The aforementioned deterministic differentiable physics layer only provides the electron stopping power ratio and does not include the primary proton flux decay and secondary particle energy deposition caused by nuclear reactions. This contribution is implicitly characterized by a beam-by-beam physical dose branch trained with Monte Carlo physical dose labels; the range-energy relationship upon which the remaining energy characterization is based is also labeled with Monte Carlo data including the nuclear reaction process. Therefore, the approximation of the deterministic differentiable physics layer does not lead to the omission of nuclear reaction contributions at the total physical dose level.

[0045] Furthermore, regarding the processing of mixed voxels, when there is partial volume mixing of interface voxels, ρ e It can be linearly superimposed according to the number of electrons, while the average ionization energy should follow the Bragg additive relationship of ln I weighted by electron density, that is, it is determined by the "average Z" of the mixed voxels. eff"The I obtained by directly looking up the table deviates from the I obtained by mixing according to the additive rule, and this deviation is larger at high atomic number contrast interfaces. One possible mitigation method is to determine I by weighted interpolation of ln I within the interface mask; another method is to add a mixing fraction to the effective atomic number prediction branch at the interface voxels and reconstruct I according to the Bragg additive rule. Without the above processing, the RSP of the interface voxels has a deviation related to the volume degree, which can be limited by increasing the resampling resolution or adding labels to the interface voxels."

[0046] Furthermore, the treatment plan is divided into B beam units. The feature tensor of each beam unit b encodes at least the unit propagation direction u. b Initial energy E 0,b The model includes the pencil beam center position, monitoring units or relative weights, and flux distribution; when using a range shifter, the equivalent water thickness is encoded; and the air gap is further encoded when changes in the air gap affect the beam spot. To enable the model to be migrated across accelerators and beam transport systems, the beam spot size parameters and initial energy broadening parameters in the air at the isocenter can be further encoded to characterize the initial phase space of the device.

[0047] For each beam element b, the effective RSP field is transformed to the beam coordinate system and input into the physical dose branch along with the beam characteristic tensor to obtain the beam-by-beam physical dose D. b (x). One specific implementation involves constructing a beam feature map using the initial energy of the beam unit, the pencil beam position, flux or weight, beam spot size, and energy broadening parameters. An analytical initial dose is generated based on a pre-calibrated integral depth-dose curve, the water equivalent depth obtained from the effective RSP path integration, and the transverse beam spot broadening model. The analytical initial dose, the effective RSP field, and the beam feature map are then input into a 3D encoder-decoder channel by channel, and the output is used for additive or multiplicative correction to obtain the beam-by-beam physical dose. Another implementation allows direct regression of the beam-by-beam physical dose from the effective RSP field and the beam feature map. Both implementations use the Monte Carlo physical dose under the same beam unit, the same spatial grid, and the same flux normalization conditions as the supervision label, and the output results are returned to the object coordinate system via inverse spatial transformation. Beam-by-beam range feature P b Depth dose can be determined from the central axis, dose-weighted rays, or a preset ray set: in the distal descent phase, the point where the dose first drops to 80% of the peak value of the curve is taken as the distal 80% dose position R. 80 The location resolution is improved by linear interpolation of adjacent sampling points; the location of the minimum first-order difference between dose and depth in the distal descent segment is taken as the location of the distal dose gradient extremum. After calculation for multiple rays, the integrated dose of each ray is normalized and weighted, and the obtained location, the signed distance from the voxel to that location, and the neighborhood dose and gradient features are encoded as P. b .

[0048] The network backbone can employ 3D convolution, encoder-decoder, or attention structures, but must retain the beam unit b and its aggregation relationships within the same object space. To control the data scale of each beam label, the magnitude of B can be limited in the phantom training set. A sub-network with shared weights can be used for each beam unit, and sparse beam units can be merged according to energy levels. Specifically, spatially registered high- and low-energy X-ray spectral CT images and their differences, ratios, or neighborhood statistics serve as dual-energy voxel features input to the shared feature encoder. The shared feature encoder forms multi-scale shared features through 3D convolution or equivalent spatial feature extraction structures. The relative electron density prediction branch and the effective atomic number prediction branch each use a voxel-by-voxel regression output head to decode the multi-scale shared features into ρ. e With Z eff Voxel plot. During training, the two output heads use supervision labels provided by digital material phantoms or material calibration data, and simultaneously receive RSP physical consistency gradients propagated back through a deterministic differentiable physics layer.

[0049] Furthermore, along the path from beam element b to voxel x The effective RSP of the integral yields the equivalent path length of water. The remaining energy characterization is determined by the range-energy function F. The specific formula is as follows:

[0050] ;

[0051] ;

[0052] in, Representing a path The unit of length integration, The effective RSP is indicated; the range-energy function F can be derived using Monte Carlo methods or a single survey table based on measurement data, and interpolated between beam energy layers. This represents the inverse of the range-energy function F. If the range shifter is located upstream of the object, its equivalent water thickness is taken into account. ;when minus When the energy is less than zero, the remaining energy is set to zero and the voxel is marked as outside the beam element's range. It should be noted that... For characterizing the remaining energy corresponding to beam unit b, the path integral is insufficient to fully recover the voxel-by-voxel energy spectrum and secondary particle composition; in implementations requiring more complete energy spectrum information, an additional [function / mechanism] can be added. Additional input channels are formed by the energy spectrum broadening characterization quantity determined by the interface gating intensity, or by the spectral shape parameters calibrated by Monte Carlo.

[0053] Furthermore, the localized non-uniformity of the medium referred to in this invention refers to the characteristic of ρ... e With Zeff Spatial gradient magnitude, gradient direction and ρ in the neighborhood of voxel e or Z eff The feature quantity is composed of at least one of the variances. The candidate set of interfaces can be formed by ρ. e or Z eff The spatial gradient magnitude is determined by voxels, material partition boundaries, or a combination of both, where the voxel size exceeds a preset threshold. Before calculating the gradient, local smoothing matching the voxel size can be used to suppress noise. For each voxel to be calculated, the nearest interface point is determined from the candidate interface set, and the normalized direction of the smoothed gradient at that interface point is used as the initial interface normal. If the inner product of the initial interface normal and the beam propagation direction is negative, the opposite direction of the normal is taken, so that the inner product of the unified interface normal and the beam propagation direction is non-negative. This process only unifies the sign of the normal and does not change the position of the physical interface.

[0054] For beam element b, the upstream medium ρ is obtained in the opposite direction of the beam propagation direction within a preset sampling distance or sampling window. e With Z eff The ρ of the downstream medium is obtained in the positive direction of the beam propagation within the corresponding sampling distance or sampling window. e With Z eff When using a sampling window, the mean or median within the window can be taken. The sampling order between upstream and downstream is uniquely determined by the beam propagation direction and cannot be interchanged. The signed distance d from the voxel to the nearest interface. s The downstream side is positive and the upstream side is negative, and the incident angle is determined by the unified interface normal and the beam propagation direction.

[0055] The direction-dependent interface variable q(x,b) includes at least ρ e,up ρ e,down Z eff,up Z eff,down d s Angle with incident light, ρ e,up ρ represents the relative electron density of the upstream medium. e,down Z represents the relative electron density of the downstream medium. eff,up Z represents the effective atomic number of the upstream medium. eff,down Indicates the effective atomic number of the downstream medium, and may include the remaining energy characterization E. rem,b Each channel is normalized according to the statistics of the training set and then spliced ​​in a fixed order. Since the interface normal orientation and the upstream and downstream sampling order change when the beam is reversed, the same physical interface will have different q(x,b) under opposite incident directions.

[0056] S103. Determine the gating strength that decreases with increasing voxel-to-interface distance based on the direction-related interface variables, and modulate the LET beam by beam using the gating strength. d The residuals are used to obtain LET bundles.d,b ;

[0057] Specifically, interface gating strength Determine using the following formula:

[0058] ;

[0059] In the formula, satisfy σ is the sigmoid function or other bounded function taking values ​​from 0 to 1; h is a differentiable mapping, trained or calibrated from the interface Monte Carlo reference data. Indicates to Perform differentiable mapping; Φ is the spatial localization factor. Follow Monotonically non-increasing, and in Exceeding the preset radius of action R int When the time approaches zero, Represents the absolute value symbol. Indicates the voxel in beam unit b The effective symbolic distance, one possible form is:

[0060] ;

[0061] or, ;

[0062] in, R is the characteristic decay length of the Gaussian spatial localization factor. int To truncate the preset radius of action of the linear spatial localization factor. Considering the LET within the interface neighborhood. d The spatial scale of the deviation from the extrapolated value of the homogeneous medium is mainly determined by the transverse and longitudinal perturbation scales caused by multiple Coulomb scattering and the renormalization distance of the energy spectrum behind the interface, R int The interface Monte Carlo reference data is calibrated according to material pairs, incident direction, and energy range, and remains consistent throughout the training, validation, and application phases. The gating strength represents the degree and spatial range of the interface correction, without a preset correction direction. The purpose of setting Φ is to ensure that the interface correction decays to zero with distance from the interface, avoiding interface-independent biases in deep homogeneous media.

[0063] Characterize the remaining energy E rem,b Localized non-uniformity of the medium, P b The direction-dependent interface variable q(x,b) is then input into LET after being stitched together by channel on the same beam coordinate grid. d Proxy branch. LET d The proxy branch can employ a 3D encoder-decoder that shares parameters among all beam units, and set up a base value output head and a residual output head: the base value output head uses a non-negative output to obtain the base value LET.d,b,base (x); The residual output head uses linear output to obtain f. b (x), then through amplitude-bounded mapping, forms a residual whose sign can be positive or negative. The basic value output head can first be supervised by a uniform medium phantom or a voxel located beyond a preset radius of action from the interface; the residual output head is generated by voxel-by-voxel LET of the interface phantom relative to the matched uniform medium phantom. d The difference is monitored, and a constraint is applied to homogeneous media and voxels far from the interface, making the interface correction approach zero, to reduce mutual compensation between the fundamental and residual terms. The residuals, after amplitude-bounded mapping, work in conjunction with the gating strength, as described in the following formula:

[0064] ;

[0065] ;

[0066] in, f represents the residual formed by the amplitude-bounded mapping. b (x) is LET d The residual output of the proxy branch, A max The upper limit of the residual magnitude is determined by the Monte Carlo reference data of the interface. It is the hyperbolic tangent function. This indicates the LET corresponding to beam element b. d , from 0≤ ≤1 and | |≤A max The interface modulation amount can be obtained to satisfy | |≤A max Furthermore, it is limited to the interface neighborhood, meaning the interface correction is bounded in magnitude and spatially confined to the interface neighborhood, while its sign is unrestricted. A max Determined by interface Monte Carlo reference data (unit: keV / μm), and can be calibrated separately for each material pair.

[0067] It can take positive or negative values, and can be determined by the interface Monte Carlo LET under different material sequences, different incident angles, and positive and negative incident conditions. d Label training or calibration. Therefore, interfaces where the effective atomic number undergoes a positive or negative step transition, such as those between low and high relative electron density media, and between high and low relative electron density media, as well as voxels upstream and downstream of these interfaces, are not forced to adopt corrections in the same direction. In another embodiment, the gating intensity can be applied to the LET channel by channel. d The intermediate characteristics of the proxy branch, but the final output still corresponds to the beam-by-beam physical dose D. b One-to-one corresponding LET d,b .

[0068] S104. Summing the individual physical doses to obtain the total physical dose, LET for each bundle. d,b Total LET was obtained by weighting based on physical dose. d Within a voxel range where the total physical dose is not lower than a preset dose threshold, the RBE-weighted dose is calculated based on the physical dose per fraction, and the effective RSP, total physical dose, and total LET are output. d And RBE weighted dose.

[0069] Specifically, after all beam elements are transformed back to the object coordinate system, the total physical dose... With the total LET d Determine using the following formula:

[0070] ;

[0071] ;

[0072] in, The physical dose of beam unit b, this weighted form is given by LET. d The dose-average value is directly defined, and its function is to ensure that the direction-related corrections of each beam element are correctly synthesized under their own dose weights. This is a necessary step for this invention to work under multi-field conditions. In the formula... A positive number is set to prevent the denominator of the total physical dose weighting formula from being zero.

[0073] Simultaneously set the preset dose threshold D thr The proton beam exhibits low dose and LET levels beyond the end of its range in the lateral penumbra. d In regions with high and significant Monte Carlo statistical fluctuations, a low threshold will retain a large number of high-noise voxels. One implementation method uses D... thr The specific ratio is determined according to a preset proportion of the total maximum physical dose or prescribed dose, and remains consistent throughout the training, validation, and application phases. Compared to D thr Smaller by at least two orders of magnitude, making It only serves to prevent division by zero and avoid non-numerical propagation in voxels where the total physical dose approaches zero, and in reporting LET... d The influence of voxels on the weighted result is negligible; for D(x) <D thr Voxels do not report LET d The RBE is used, but is marked with a mask. This preset dose threshold is used for numerical stability and display range limitation, and does not represent a clinical dose criterion; the threshold value remains consistent during training, validation, and application phases.

[0074] Within the voxel not masked by low dose, the total physical dose According to the treatment plan, the number of sessions n fxThe physical dose per fraction, d(x), is calculated using the following formula:

[0075] ;

[0076] Again with the general LET d The reference photon α / β parameters are also input into the variable RBE parameterization model. Taking the phenomenological model based on a linear quadratic framework as an example, the specific formula is as follows:

[0077] ;

[0078] In the formula, α / β is a reference parameter corresponding to the selected tissue type and biological endpoint, and RBE max With RBE min The largest and minimum The parameterization relationships preset by each model depend on LET. d And α / β. It should be noted that the RBE expression based on the linear quadratic framework above is derived from the cell survival relationship of a single irradiation, and its dose variable is the physical dose per fraction; if the total dose of the treatment course is mistakenly substituted, the RBE will approximate the given RBE. min This leads to a significant underestimation of the biological dose increase at the end of the range and at the interface.

[0079] This invention does not present incomplete accuracy, speed, or clinical benefits as achieved results. Differences between candidate variable RBE parameterized models only characterize model selection sensitivity; voxel generation cue masks that exceed material calibration ranges or training distributions are not automatically considered as directly acceptable calculation results.

[0080] RBE weighted dose Give the formula as follows, and indicate the corresponding fractional scheme:

[0081] ;

[0082] Furthermore, for the matching of the LET definition, the LET definition upon which the selected variable RBE parameterization model is based during its calibration must be consistent with the LET definition of this invention. d The scoring definition is consistent in three aspects: the type of charged particles involved in the counting, the medium reference (water or local medium), and whether track truncation is used. If they are inconsistent, the LET should be adjusted according to the pre-defined conversion relationship. d The values ​​are then converted to the definition used in this model and substituted back into the model. A mismatch between the two will introduce a systematic bias in the high LET region, which is precisely the range terminus and interface neighborhood of interest in this invention. One feasible conversion and calibration method is as follows: In the same batch of Monte Carlo simulations, for the same irradiation conditions, two sets of LET values ​​are output in parallel according to the scoring definition of this invention and the scoring definition used in the target RBE model. dBased on this, a conversion factor table is established according to material category and remaining energy range. When applying it, the conversion factor is determined by the material of the voxel and the E... rem,b Interpolation is used; the conversion factor table is determined before training and saved with the model.

[0083] Furthermore, α / β is a reference parameter corresponding to the selected tissue type or biological endpoint, rather than a definitive true value for the individual subject. One implementation involves assigning α / β parameter maps structurally according to the structural profile of the treatment plan; when structural information is unavailable, a uniform preset value is used, and this hypothesis is output along with the results; when α / β is unavailable, the interval result of the RBE-weighted dose is output within a preset value range.

[0084] When multiple candidate variable RBE models are used in parallel, each model is calculated independently. The range, standard deviation or quantile interval of the results only characterize the sensitivity of model selection and are used to indicate that the calculation results in this region are more sensitive to model selection. Typically, this is located near the end of the range and the strong non-uniform interface.

[0085] It should also be noted that before implementing the proton radiotherapy dose and biological effect calculation method based on dual-energy CT of this invention, corresponding modeling should be performed. The modeling includes two stages: model establishment and model application. In the model establishment stage, the relative electron density and effective atomic number are decoupled, deterministic relative stopping power is calculated, beam-by-beam physical dose and residual energy are characterized, direction-dependent interface dose is determined, average linear energy transfer is corrected, and the weighted aggregation relationship of multi-beam dose is determined. In the model application stage, the physical dose, dose-average linear energy transfer, and relative biological effect weighted dose are output based on the established model.

[0086] In addition, the minimum training data for this invention includes: spatially consistent dual-energy CT images or generated synthetic images; ρ e Z eff Or material calibration data; effective RSP tags for treatment planning; beam-by-beam physical dose D of at least one homogeneous medium phantom and one two-material interface phantom in multiple beam directions. b With LET d,b Tags. Monte Carlo ratings use a unified LET system. d By definition, the training and validation sets use the same particle selection, cutoff energy, and scoring rules. The interface phantom must contain illumination data with opposite incident directions and different incident angles to calibrate the relationship between direction-dependent interface variables and potentially positive or negative residuals, and to determine A accordingly. max With R int Interface residual labels can be generated from the LET bundles of the interface modulus. d Subtract the beam-by-beam LET of the matched homogeneous medium modulus at the same initial beam energy, the same flux normalization conditions, and the same water equivalent depth. dWhen it is not possible to construct strictly paired samples, zero interface correction constraints are set for uniform media and voxels that are more than the radius of action away from the interface, and the basic value loss and residual regularization terms are recorded respectively.

[0087] The overall network model of this invention includes at least a shared feature encoder, a relative electron density prediction branch, an effective atomic number prediction branch, a deterministic differentiable physical layer, a physical dose branch, and a LET. d Modules such as proxy branches optimize and train the parameters of each module by constructing objective functions.

[0088] The input to the network model can be represented as an object voxel tensor. With the beam characteristic tensor X of group B beam,b ,in, , , The dimensions of the voxel mesh are defined in three directions, and C is the number of input channels; the output includes ρ e Z eff Effective RSP, D b P b LET d,b and the aggregated D and LET d .

[0089] The training objective function includes at least the first three terms of the following formula. In one embodiment, the training objective function is shown in the following formula:

[0090] ;

[0091] Among them, L RSP The RSP physical consistency loss is used to compare the deterministic differentiable physical layer output with the valid RSP label, L. D For physical dose monitoring loss, used to compare beam-by-beam or total physical dose, L LET For LET d Supervised loss, used to compare bundle-by-bundle LET d or dose-weighted total LET d L ρ With L Z For optional physical parameter monitoring. , , , , L respectively RSP、 L D、 L LET、 L ρ、 L ZThe weights corresponding to each loss can be set according to actual needs. Edge-aware RBE smoothing regularization is not a necessary training term to avoid imposing unfounded smoothing constraints on real material boundaries, α / β variation regions, or low-dose regions.

[0092] One training sequence is: first train ρ e Z eff With the RSP chain; then fix all or part of its parameters, train the beam-by-beam physical dose branch and LET. d Agent branch; final joint fine-tuning. Course learning, attention module, hard sample weighting, and dynamic loss weights can all be used. Training LET d When branching, the basic value output head can be trained with uniform medium data and its parameters fixed first. Then, the residual output head and the interface gating mapping can be trained with paired interface residual labels. Finally, joint fine-tuning is performed while retaining the zero interface correction constraint.

[0093] Furthermore, the application phase only calls the trained model and the material calibration relationships, effective RSP definitions, and range-energy relationships fixed during training. The dual-energy CT images of the object to be calculated and the planned beam parameters are input, and ρ is output sequentially. e Z eff Effective RSP, beam-by-beam and total physical dose, beam-by-beam and total LET d RBE and RBE-weighted dose; model parameters are not refitted during the application phase.

[0094] Consistency checks include: ρ e The nonnegativity of RSP and physical dose; Z eff Whether I is within the material calibration range; low-dose mask; out-of-distribution detection during training. A clear out-of-distribution detection method is as follows: calculate the mean and covariance of the training sample features in the preset feature layer of the shared encoder, calculate the Mahalanobis distance of voxel or local block features during application, and generate an out-of-distribution cue mask when it exceeds the validation set quantile threshold.

[0095] For tissue equivalent material or phantom calibration scenarios, the effective RSP can be further compared with independently measured or validated calculation results. Three-dimensional Gamma pass rate, R... 80 Deviation, dosage difference, LET d Differences and multi-beam aggregation errors are both usable verification metrics; in the absence of actual measurements or reference calculation results, accuracy or clinical efficacy values ​​should not be entered.

[0096] Compared to simply combining multiple existing modules side-by-side, this invention, through beam unit b, integrates the remaining energy state determined by the effective RSP and beam parameters, the range and far-end gradient characteristics extracted by the physical dose branch, the ordered material information on both sides of the interface, and the LET. dThe contribution connection is a causal chain with a clear data source and computational destination; and the positive or negative residual correction driven by the direction-related interface variable, with bounded amplitude and limited to the interface neighborhood, is transferred across tasks from bundle transport state to biological effector.

[0097] This invention also provides a proton radiotherapy dose and biological effect calculation system based on dual-energy CT, corresponding to the above-described method embodiments. Since the system embodiments are basically similar to the method embodiments, the description is relatively simple. For details of the relevant technical features and the effects achieved, please refer to the corresponding descriptions of the above-described method embodiments. This invention discloses a proton radiotherapy dose and biological effect calculation system based on dual-energy CT, such as... Figure 2 As shown, the system includes:

[0098] The dual-energy spectral physics parameter calculation module is used to calculate the relative electron density and effective atomic number of each voxel based on the high-energy spectral CT images and low-energy spectral CT images of the object to be calculated, and obtain the effective RSP through a deterministic differentiable physical layer.

[0099] The module for calculating beam-by-beam physical dose and direction-dependent interface variables is used to break down the treatment plan into multiple beam units, obtain beam-by-beam physical dose and residual energy characterization for each beam unit, and orderly sample the upstream and downstream media of the interface according to the beam propagation direction to form direction-dependent interface variables. The beam unit is a set of beams with a common incident direction and known initial energy distribution. The residual energy characterization is determined by the effective RSP path integral and the pre-calibrated range-energy relationship.

[0100] LET by bundle d,b The calculation module is used to determine the gating intensity, which decreases with increasing voxel-to-interface distance, from the direction-related interface variables, and to modulate the beam-by-beam LET using the gating intensity. d The residuals are used to obtain LET bundles. d,b ;

[0101] The multi-beam aggregation and bio-dosage calculation module is used to add up the physical doses of each beam to obtain the total physical dose, and to calculate the LET of each beam. d,b Total LET was obtained by weighting based on physical dose. d Within a voxel range where the total physical dose is not lower than a preset dose threshold, the RBE-weighted dose is calculated based on the physical dose per fraction, and the effective RSP, total physical dose, and total LET are output. d And RBE weighted dose.

[0102] The system includes a dual-spectrum physical parameter calculation module, a beam-by-beam physical dose and direction-related interface variable calculation module, and a beam-by-beam LET module. d,bThe system comprises a calculation module, a multi-beam convergence and bio-dosage calculation module, and a result output and consistency check module. The beam-by-beam physical dose and direction-related interface variable calculation module performs calculations for each beam element. After each beam element is transformed into the same object coordinate system, the multi-beam convergence and bio-dosage calculation module performs physical dose summation and LET calculations. d Dose-weighted aggregation. Modules can exchange images, structures, beam currents, doses, and cue masks using DICOM, DICOM-RT, NIfTI, or preset array formats. The system section describes the calculation modules and their data relationships, excluding claims involving electronic devices or computer-readable storage media.

[0103] Without changing the dual-energy spectrum ρ e / Z eff Decoupling, deterministic RSP definition, beam-by-beam physical state and beam-by-beam range characteristics, direction-dependent interface variables sampled in an ordered manner according to beam direction, and amplitude-bounded LETs that can be positive or negative and limited to the interface neighborhood. d Residual correction, and multi-beam LET d Under the premise of dose-weighted relationships, a specific network backbone, Monte Carlo platform, range-energy lookup table, material interpolation method, and input-compatible RBE parameterized model can be used to replace the existing model. The scope of protection of this invention is defined by the claims.

[0104] It should be noted that the order of the above embodiments of the present invention is merely for descriptive purposes and does not represent the superiority or inferiority of the embodiments. Furthermore, specific embodiments have been described above. Other embodiments are within the scope of the appended claims. In some cases, the actions or steps described in the claims can be performed in a different order than that shown in the embodiments and still achieve the desired result. Additionally, the processes depicted in the drawings do not necessarily require a specific or sequential order to achieve the desired result. In some embodiments, multitasking and parallel processing are also possible or may be advantageous.

[0105] It should be understood that the above description of the preferred embodiments is quite detailed and should not be construed as a limitation on the scope of protection of the present invention. Those skilled in the art can make substitutions or modifications under the guidance of the present invention without departing from the scope of protection of the claims of the present invention, and all such substitutions or modifications fall within the scope of protection of the present invention. The scope of protection of the present invention should be determined by the appended claims.

Claims

1. A method for calculating proton radiotherapy dose and biological effects based on dual-energy CT, characterized in that, The method includes: The relative electron density and effective atomic number of each voxel are calculated based on the high-energy X-ray spectral CT images and low-energy X-ray spectral CT images of the object to be calculated, and the effective RSP is obtained through a deterministic differentiable physical layer. The treatment plan is divided into multiple beam units, and the physical dose and residual energy characterization of each beam unit are obtained. The upstream and downstream media of the interface are sampled in an orderly manner according to the beam propagation direction to form direction-related interface variables. The beam unit is a set of beams with a common incident direction and known initial energy distribution. The residual energy characterization is determined by the effective RSP path integral and the pre-calibrated range-energy relationship. determining a gating strength that decays with increasing voxel-to-interface distance from the direction-dependent interface variable, and modulating the per-beam LET by the gating strength d residual, resulting in a per-beam LET d,b ; The total physical dose is obtained by adding up the physical doses of each bundle, and the LET is calculated by adding up the physical doses of each bundle. d,b Total LET was obtained by weighting based on physical dose. d Within a voxel range where the total physical dose is not lower than a preset dose threshold, the RBE-weighted dose is calculated based on the physical dose per fraction, and the effective RSP, total physical dose, and total LET are output. d And RBE weighted dose.

2. The method for calculating proton radiotherapy dose and biological effects based on dual-energy CT according to claim 1, characterized in that, The process of calculating the relative electron density and effective atomic number of each voxel based on high-energy X-ray spectral CT images and low-energy X-ray spectral CT images of the object to be calculated, and obtaining the effective RSP through a deterministic differentiable physical layer, includes: The dual-spectral voxel features formed by the high-energy X-ray spectral CT images and low-energy X-ray spectral CT images are input into the shared feature encoder, and the voxel-level relative electron density and effective atomic number are obtained through the relative electron density prediction branch and the effective atomic number prediction branch, respectively. The voxel-level relative electron density and effective atomic number are input into a deterministic differentiable physical layer. The average ionization energy of the medium is determined based on a pre-calibrated material parameter relationship. The effective RSP is then calculated within a calibrated energy range based on the average ionization energy of the medium. The material parameter relationship represents the correspondence between the effective atomic number and the average ionization energy of the medium.

3. The method for calculating proton radiotherapy dose and biological effects based on dual-energy CT according to claim 2, characterized in that, The process involves inputting the voxel-level relative electron density and effective atomic number into a deterministic differentiable physical layer, determining the average ionization energy of the medium based on pre-calibrated material parameter relationships, and calculating the effective RSP within a calibrated energy range based on the average ionization energy of the medium, using the following formula: ； in, The energy of a proton is... At that time, voxels The relative proton stopping power of the corresponding medium relative to water is determined based on the electronic stopping process generated by the interaction between protons and electrons in the medium. Represents voxels, Represents the relative electron density relative to water; with I w Representing voxels The corresponding average ionization energy of the medium and water; γ is the ratio of the proton's velocity to the speed of light. p This represents the Lorentz factor corresponding to the proton velocity. and γ p All are determined by the proton energy E; m e c 2 It represents the rest energy of an electron.

4. The method for calculating proton radiotherapy dose and biological effects based on dual-energy CT according to claim 1, characterized in that, The calculation process for the beam-by-beam physical dose includes: For each beam element, the effective RSP field is transformed to the beam coordinate system and input into the physical dose branch along with the beam feature tensor to obtain the beam-by-beam physical dose. The beam feature tensor is encoded by at least the beam injection direction, initial energy, pencil beam position, weights, and flux distribution.

5. The method for calculating proton radiotherapy dose and biological effects based on dual-energy CT according to claim 1, characterized in that, The remaining energy characterization is determined by the effective RSP path integral and a pre-calibrated range-energy relationship, including: For each beam element, the effective RSP is integrated along the beam path to obtain the water equivalent path length, and the remaining energy characterization is obtained by combining the initial energy of the corresponding beam element and the pre-calibrated range-energy relationship.

6. The method for calculating proton radiotherapy dose and biological effects based on dual-energy CT according to claim 1, characterized in that, The orderly sampling of the upstream and downstream media of the interface according to the beam propagation direction constitutes direction-dependent interface variables, including: The upstream and downstream media of the interface are determined according to the beam propagation direction. The relative electron density, effective atomic number, interface normal, signed distance from the voxel to the interface, and angle between the beam injection direction and the interface normal of the upstream and downstream media constitute the direction-related interface variables. For each voxel to be calculated, the nearest interface point is determined in the interface candidate set, and the normalized direction of the smoothed gradient at the nearest interface point is used as the initial interface normal. If the inner product of the initial interface normal and the beam propagation direction is negative, the opposite direction of the initial interface normal is taken to determine the interface normal.

7. The method for calculating proton radiotherapy dose and biological effects based on dual-energy CT according to claim 6, characterized in that, The gating strength, which decreases with increasing voxel-to-interface distance, as determined by the direction-related interface variable, is obtained by the following formula: ； in, Let represent the gating strength, σ be a bounded function taking values ​​from 0 to 1, and h be a differentiable mapping. Represents the interface variables related to the direction. Perform a differentiable mapping, where Φ is the spatial localization factor. Indicates the voxel in beam unit b The effective symbolic distance, Indicates following Monotonically non-increasing, and in When the value exceeds the preset radius of action, it tends to zero. Represents the absolute value symbol.

8. The method for calculating proton radiotherapy dose and biological effects based on dual-energy CT according to claim 7, characterized in that, LET beam-by-beam modulation via gated intensity modulation d The residuals are used to obtain LET bundles. d,b include: The residual energy characterization, local medium inhomogeneity characteristics, beam-by-beam range characteristics, and direction-dependent interface variables are stitched together on the same beam coordinate voxel grid by channel and then input into LET. d proxy branch, wherein the LET d The proxy branch includes a base value output head and a residual output head. The base value output head outputs the base value under the supervision of a uniform medium phantom or a voxel located at a distance from the interface exceeding a preset radius. The residual output head uses linear output and then performs amplitude-bounded mapping to form a bundle-by-bundle LET with a sign that can be positive or negative. d residual; Based on the aforementioned baseline value and gating strength and bundle-by-bundle LET d The sum of the residual products yields bundle-by-bundle LET. d,b .

9. The method for calculating proton radiotherapy dose and biological effects based on dual-energy CT according to claim 1, characterized in that, The calculation of the RBE weighted dose based on the physical dose per fraction includes: The total physical dose is converted into a per-session physical dose according to the treatment plan, and the per-session physical dose and total LET are calculated. d The RBE is obtained by inputting the reference photon α / β parameters corresponding to the selected tissue into a variable RBE parameterization model; The product of the RBE and the total physical dose is used as the RBE weighted dose.

10. A proton radiotherapy dose and biological effect calculation system based on dual-energy CT, characterized in that, The system includes: The dual-energy spectral physics parameter calculation module is used to calculate the relative electron density and effective atomic number of each voxel based on the high-energy spectral CT images and low-energy spectral CT images of the object to be calculated, and obtain the effective RSP through a deterministic differentiable physical layer. The module for calculating beam-by-beam physical dose and direction-dependent interface variables is used to break down the treatment plan into multiple beam units, obtain beam-by-beam physical dose and residual energy characterization for each beam unit, and sample the upstream and downstream media of the interface in an orderly manner according to the beam propagation direction to form direction-dependent interface variables. The beam unit is a set of beams with a common incident direction and a known initial energy distribution. The residual energy characterization is determined by the effective RSP path integral and the pre-calibrated range-energy relationship. LET by bundle d,b The calculation module is used to determine the gating intensity, which decreases with increasing voxel-to-interface distance, from the direction-related interface variables, and to modulate the beam-by-beam LET using the gating intensity. d The residuals are used to obtain LET bundles. d,b ; The multi-beam aggregation and bio-dosage calculation module is used to add up the physical doses of each beam to obtain the total physical dose, and to calculate the LET of each beam. d,b Total LET was obtained by weighting based on physical dose. d Within a voxel range where the total physical dose is not lower than a preset dose threshold, the RBE-weighted dose is calculated based on the physical dose per fraction, and the effective RSP, total physical dose, and total LET are output. d And RBE weighted dose.

Images