A sparse proton arc radiotherapy plan planning method based on irradiation time homogenization

Abstract

The present application belongs to the technical field of proton radiotherapy, and particularly relates to a sparse proton arc radiotherapy plan planning method based on irradiation time homogenization, which comprises the following steps: obtaining an initial radiotherapy plan, including a beam spot weight vector, beam information, and a dose influence matrix corresponding to the beam information; constructing a beam spot 01 grouping matrix B according to the beam information, with rows corresponding to gantry angles and columns corresponding to beam spots; taking the discrete degree of total beam spot weights between different gantry angles as a, wherein x represents the beam spot weight vector; P is a centralization matrix; taking the product of the dose influence matrix and the beam spot weight vector as an actual dose distribution, and constructing beam spot distribution sparsity based on the beam spot weight vector; and taking the weighted sum of the deviation between the actual dose distribution and a prescribed dose distribution, the beam spot distribution sparsity, and the discrete degree as a minimum as an objective, and optimizing and solving the beam spot weight vector. The present application can simultaneously improve overall treatment delivery efficiency and gantry operation stability.

Description

Technical Field

[0001] This invention belongs to the technical field of proton radiotherapy, and more specifically, relates to a sparse proton arc radiotherapy planning method based on irradiation time homogenization. Background Technology

[0002] Proton arc therapy is an emerging proton radiotherapy technique that delivers proton beams by continuously rotating the gantry (up to 360°) and controlling the beam at multiple angles. This increases the freedom of dose modulation, resulting in better target dose coverage and protection of normal tissues. Due to the Bragg peak characteristics of the proton beam and the selectable control angles, proton arc therapy offers significant advantages over traditional fixed-field proton therapy in terms of dose conformity and protection of normal tissues.

[0003] However, proton arc therapy typically involves numerous control angles, energy layers, and beam spots. The number of beam spots and their weight distribution vary significantly across different control angles, leading to uneven delivery times at each angle. This, in turn, causes fluctuations in gantry rotation speed and reduces overall treatment delivery efficiency. Existing optimization methods primarily focus on dose distribution quality and beam spot weight sparsity, lacking synergistic optimization for treatment delivery efficiency and gantry motion stability. While optimized treatment plans may meet dosimetric requirements, they may still exhibit significant gantry speed fluctuations and marked differences in control point delivery times, resulting in increased total treatment time, frequent gantry acceleration and deceleration, and increased mechanical load on the equipment. Furthermore, longer treatment times increase the impact of patient positional drift and physiological movement on dose delivery accuracy, reducing the stability and repeatability of the treatment process. Therefore, while ensuring dosimetric quality and beam spot sparsity, further improving treatment delivery efficiency and maintaining stable gantry motion are crucial for enhancing treatment accuracy, improving patient experience, increasing equipment utilization, and enhancing clinical application value. Summary of the Invention

[0004] In view of the above-mentioned defects or improvement needs of the existing technology, the present invention provides a sparse proton arc radiotherapy planning method based on irradiation time homogenization, which aims to improve the overall treatment delivery efficiency and gantry operation stability at the same time.

[0005] To achieve the above objectives, according to one aspect of the present invention, a method for planning sparse proton arc radiotherapy based on irradiation time homogenization is provided, comprising: Obtain an initial radiotherapy plan, including a beam weight vector, beam information, and a dose influence matrix corresponding to the beam information. The beam information includes the energy layer corresponding to each gantry angle and the beam spot corresponding to each energy layer. The dimension of the beam weight vector is the same as the total number of known beam spots. The number of rows in the dose influence matrix is ​​the same as the total number of voxels divided in the target area, and the number of columns is the same as the total number of known beam spots. Based on the beam information, a beam spot 0-1 grouping matrix B is constructed, where the rows correspond to the frame angles and the vertical columns correspond to the beam spots. A 1 indicates that the beam spot corresponding to the frame angle participates in illumination, and a 0 indicates that the beam spot corresponding to the frame angle does not participate in illumination. Multiplying P, B, and x together yields... , representing the degree of dispersion of the total beam spot weights across different frame angles, where x represents the beam spot weight vector; P is the centering matrix, , This represents a diagonal matrix with a diagonal value of 1. Represents a vector consisting entirely of 1s; The product of the dose influence matrix and the beam spot weight vector is used as the actual dose distribution, and the beam spot distribution sparsity is constructed based on the beam spot weight vector. The beam spot weight vector is optimized by minimizing the weighted sum of the deviation between the actual dose distribution and the prescription dose distribution, the beam spot distribution sparsity, and the aforementioned dispersion.

[0006] Furthermore, the method for determining the beam information in the initial radiotherapy plan and the corresponding dose-effect matrix is ​​as follows: Generate an initial energy layer set corresponding to each preset rack angle and the beam spot information corresponding to each energy layer in the set to form the initial beam information; The energy layers are arranged according to the rack angle order. Only one energy layer is selected for each rack angle, and the energy layers selected from the rack angles from smallest to largest are selected from highest to lowest. If there is no energy layer in the initial energy layer set corresponding to a certain rack angle that is lower than the energy layer already selected in the adjacent previous rack angle, then the highest energy layer is selected from the initial energy layer set, and the selection operation of selecting energy layers from the rack angles from smallest to largest in descending order is repeated. Based on the energy layer corresponding to each rack angle and the beam spot information corresponding to that energy layer, the dose contribution relationship of each beam spot to each voxel in the target area is calculated, and the dose influence matrix is ​​obtained.

[0007] Furthermore, the method for determining the beam information in the initial radiotherapy plan and the corresponding dose-effect matrix is ​​as follows: Generate the energy layer and its beam spot information corresponding to each preset rack angle to form the initial beam information; Based on the initial beam information and treatment equipment parameters, the water equivalent path length vector corresponding to each gantry angle is determined, where each element represents the water equivalent path length under the corresponding energy layer; the water equivalent path length vectors corresponding to different angles are completed to form vectors of the same dimension, and the gantry angles are clustered by the deviation between vectors, and the gantry angles whose deviation between vectors meets the preset threshold are divided into the same angle group. Based on the clustering results, the energy layers were optimized and organized according to the rack angle order; Based on the energy layer and its corresponding beam spot information corresponding to each rack angle after processing, the dose contribution relationship of each beam spot to each voxel in the target area is calculated to obtain the dose influence matrix.

[0008] Furthermore, the optimization problem of the beam spot weight vector is as follows:

[0009]

[0010] In the formula, This represents the dose effect matrix; This represents the prescription dose vector composed of the doses at each voxel in the target region; Represents the beam spot weight vector; Represents the 0-1 grouping matrix of the beam spot; Represents a centered matrix; Indicates the total number of beam spots; Indicates the weight of the sparse constraint on the beam pattern; This represents the weight of the angle delivery time balance term, i.e., the weight of the dispersion term; This represents the minimum number of hops threshold.

[0011] Furthermore, the alternating direction multiplier method is used to optimize the solution of the beam spot weight vector, as follows: Optimize the following issues:

[0012]

[0013] In the formula, This represents the dose effect matrix; This represents the prescription dose vector composed of the doses at each voxel in the target region; Represents the beam spot weight vector; This indicates the sparse distribution of the fascicles; Represents the 0-1 grouping matrix of the beam spot; Represents a centered matrix; Indicates the total number of beam spots; Indicates the weight of the sparse constraint on the beam pattern; This represents the weight of the angle delivery time balance term, i.e., the weight of the dispersion term; Indicates taking The L2 norm; This represents the minimum number of hops threshold; Rewritten as:

[0014]

[0015] Among them, auxiliary variables Used for separation Sparse terms, auxiliary variables Used to separate minimum hop count constraints;

[0016] The optimal beam pattern distribution is obtained by solving the alternating direction multiplier method: initialization parameters. Lagrange multipliers The superscript indicates initialization, with a preset maximum iteration parameter. Convergence tolerance and penalty parameters Then iterative updates are performed, at the [number]th [time]. In the next iteration, the following solution is obtained: Sub-problems Subproblems and Sub-problems updated :

[0017]

[0018]

[0019] Determine if the termination condition has been met; if so, complete the solution; otherwise, continue to the next iteration.

[0020] According to another aspect of the present invention, an electronic device is provided, including a memory and a processor, wherein the memory stores a computer program, and the processor executes the computer program to implement the steps of the method described above.

[0021] According to another aspect of the invention, a computer-readable storage medium is provided, the computer-readable storage medium including a stored computer program, wherein, when the computer program is executed by a processor, it controls the device where the storage medium is located to perform the steps of the method described above.

[0022] According to another aspect of the invention, a computer program product is provided, comprising a computer program or instructions that, when executed by a processor, implement the steps of the method described above.

[0023] In summary, compared with the prior art, the technical solutions conceived by this invention have the following main advantages: 1. This invention provides a treatment planning method for proton arc radiotherapy. A beam spot 0-1 grouping matrix B is constructed based on beam information, where rows correspond to gantry angles and vertical columns correspond to beam spots. A 1 indicates that the beam spot corresponding to the gantry angle participates in irradiation, and a 0 indicates that the beam spot corresponding to the gantry angle does not participate in irradiation. Multiplying P, B, and x yields... , representing the degree of dispersion of the total beam spot weights across different frame angles, where x represents the beam spot weight vector; P is the centering matrix, , This represents a diagonal matrix with a diagonal value of 1. This represents a vector of all ones. By using the degree of dispersion as a constraint on beam delivery time consistency, and co-optimizing it with the sparse distribution of beam spot weights and dose distribution, this method improves delivery efficiency and gantry rotation speed stability while ensuring the dosimetric quality of the radiotherapy plan. Specifically, gantry rotation speed stability is achieved by uniformly distributing beam spot irradiation time between control points. This invention's method, through mathematical problem-solving, simultaneously improves overall treatment delivery efficiency and gantry operational stability. Compared to existing optimization methods, this method effectively reduces frequent changes in gantry speed, improves equipment mechanical stability and control precision, and reduces overall treatment time. It provides a technical foundation for achieving uniform rotation speed dynamic delivery in proton arc radiotherapy and has promising clinical application prospects. Attached Figure Description

[0024] Figure 1 A flowchart of a sparse proton arc radiotherapy planning method based on irradiation time homogenization is provided for an embodiment of the present invention; Figure 2 This is a schematic diagram of the water equivalent path angle distribution corresponding to a specific brain tumor, provided as an embodiment of the present invention. Figure 3 Provided for embodiments of the present invention Figure 2 Actual dose-volume histogram for corresponding brain tumor cases; Figure 4 Provided for embodiments of the present invention Figure 2 Dose-slice diagrams corresponding to brain tumor cases; Figure 5 This is a schematic diagram showing the comparison results of the time smoothness of the treatment plan before and after optimization provided in the embodiments of the present invention. Among them, (a) is a time distribution comparison diagram, (b) is a time absolute deviation comparison diagram, (c) is a time distribution comparison diagram after sorting, (d) is a radiation time distribution comparison diagram, (e) is a beam extraction (i.e. radiation) time absolute deviation comparison diagram, (f) is a beam extraction (i.e. radiation) time distribution comparison diagram after sorting, (g) is a total time variation coefficient comparison diagram, (h) is a beam extraction (i.e. radiation) time variation coefficient comparison diagram, and (i) is a sparsity comparison diagram. Detailed Implementation

[0025] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in detail below with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are merely illustrative and not intended to limit the invention. Furthermore, the technical features involved in the various embodiments of this invention described below can be combined with each other as long as they do not conflict with each other.

[0026] Example 1 A sparse proton arc radiotherapy planning method based on irradiation time homogenization includes: Obtain an initial radiotherapy plan, including a beam weight vector, beam information, and a dose influence matrix corresponding to the beam information. The beam information includes the energy layer corresponding to each gantry angle and the beam spot corresponding to each energy layer. The dimension of the beam weight vector is the same as the total number of known beam spots. The number of rows in the dose influence matrix is ​​the same as the total number of voxels divided in the target area, and the number of columns is the same as the total number of known beam spots. Based on the beam information, a beam spot 0-1 grouping matrix B is constructed, where the rows correspond to the frame angles and the vertical columns correspond to the beam spots. A 1 indicates that the beam spot corresponding to the frame angle participates in illumination, and a 0 indicates that the beam spot corresponding to the frame angle does not participate in illumination. Multiplying P, B, and x together yields... , representing the degree of dispersion of the total beam spot weights across different frame angles, where x represents the beam spot weight vector; P is the centering matrix, , This represents a diagonal matrix with a diagonal value of 1. Represents a vector consisting entirely of 1s; The product of the dose influence matrix and the beam spot weight vector is used as the actual dose distribution, and the beam spot distribution sparsity is constructed based on the beam spot weight vector. The beam spot weight vector is optimized by minimizing the weighted sum of the deviation between the actual dose distribution and the prescription dose distribution, the beam spot distribution sparsity, and the aforementioned dispersion.

[0027] This embodiment proposes a method based on the beam angle... The beam spots are grouped to construct a beam spot 0-1 grouping matrix. This is used to calculate the total beam spot weights corresponding to each angle. Additionally, a centering matrix P is defined. This is used to describe the dispersion of the total beamwidth weight across different rack angles. Since the delivery time at each angle is approximately proportional to the corresponding total beamwidth weight, therefore... It can also characterize the degree of fluctuation in delivery time between different angles. Under the premise of satisfying dosimetric constraints, the degree of dispersion of the total beam spot weight between different gantry angles is constructed by establishing the coupling relationship between beam spot irradiation time and gantry rotation speed. Therefore, by introducing a beam spot distribution sparsity term and a dispersion term, the sparsity constraint and the beam spot irradiation time consistency constraint are introduced to achieve the synergistic optimization of beam spot weight and irradiation time to obtain a uniform beam delivery time distribution.

[0028] The technical solution described in this embodiment optimizes beam spot weight sparse distribution and beam delivery time consistency, improving delivery efficiency and gantry rotation speed stability while ensuring the dosimetric quality of the radiotherapy plan. Specifically, gantry rotation speed stability is achieved by uniformly distributing beam spot irradiation time between control points. Compared to existing optimization methods, this method effectively reduces frequent gantry speed changes, improves equipment mechanical stability and control precision, and reduces overall treatment time. It provides a technical foundation for achieving uniform rotation speed dynamic delivery in proton arc radiotherapy and has promising clinical application prospects.

[0029] In fact, to simultaneously improve overall treatment delivery efficiency and gantry operational stability, the technical approach proposed in this embodiment is primarily based on the following considerations: From an engineering control perspective, the delivery time can be uniformized by directly optimizing the gantry rotation speed or the time allocation of each control point. However, such methods are typically adjusted after the treatment plan is generated, failing to fully consider the impact of speed adjustments on patient dose distribution and the quality of the original plan. On the other hand, while traditional beam sparsity optimization can reduce the number of effective beams and improve treatment delivery efficiency, it may result in uneven distribution of beam number and beam intensity among different control points, leading to significant differences in the required delivery time for each control point. This makes it difficult to meet the requirements of the treatment equipment for gantry speed stability and actual motion control.

[0030] To address the aforementioned issues, this embodiment takes a mathematical optimization approach, unifying dosimetric targets, treatment delivery efficiency, and delivery time uniformity into a single objective function. While ensuring target dose coverage, it also considers spot sparsity and delivery time uniformity between different control points, thereby achieving joint optimization of treatment plan quality, delivery efficiency, and device executability.

[0031] To achieve optimization from a mathematical perspective, this embodiment creatively proposes... This describes the dispersion of the total beam spot weight across different gantry angles. By establishing the coupling relationship between beam spot irradiation time and gantry rotation speed, the dispersion of the total beam spot weight across different gantry angles is constructed, as theoretically explained below: The total delivery time of a proton beam mainly consists of spill time, beam spot switching time, energy layer switching time, and other equipment response times. Spill time accounts for a significant portion of the total delivery time and is the primary factor causing differences in delivery time at different control points. Within the operating range where the beam intensity and operating status remain stable, the spill time of a single beam spot has a linear relationship with its beam spot weight, which can be expressed as: ,in, The overflow time coefficient is the unit beam spot weight.

[0032] The total beam spot weight of each control point can be uniformly expressed as: Correspondingly, the overflow time vector for each control point is: therefore, Established beam spot weight vector The linear correspondence between the overflow time and the overflow time at each control point. Since overflow time accounts for a major proportion of the total delivery time, It can be used as a proxy quantity to characterize the main delivery time of each control point and its relative differences.

[0033] To measure the global uniformity of beam delivery time among different control points, this embodiment uses the coefficient of variation (CV) for modeling. To simplify notation, the main delivery time vector of the control points in the model is defined as: time vector The average value is: .

[0034] According to the above formula, we have:

[0035] Using the definition of population standard deviation, delivery time vector The standard deviation is:

[0036] Therefore, the coefficient of variation of delivery time is defined as:

[0037] The square of the coefficient of variation is:

[0038] Will Substituting into the above equation, we get:

[0039] in,

[0040] The average value of the total beam spot weight at each control point is defined as:

[0041] Then the expression can be written as , which represents the degree of dispersion of the total beam spot weight between different frame angles.

[0042] further, The i The components are: ,therefore, .Right now This represents the sum of squared deviations of the total beamwidth weights at each control point from their global average. When the optimization process primarily focuses on the temporal dispersion between control points, this can be represented as... proportionality coefficient in the formula Merge into the corresponding regularization parameters Therefore, in this embodiment, a further preferred method is to use... This serves as a secondary penalty term for beam delivery time non-uniformity.

[0043] This penalty term is equivalent to the sum of squared deviations of the total beam spot weights at each control point from their average values. Minimizing this term can make the total beam spot weights at different control points and their corresponding linear overflow times more consistent, thereby reducing the main delivery time differences between control points and improving the smoothness of beam delivery and gantry movement during proton arc therapy.

[0044] As a preferred embodiment, the method for determining the beam information in the initial radiotherapy plan and the dose influence matrix corresponding to that beam information is as follows: Generate an initial energy layer set corresponding to each preset rack angle and the beam spot information corresponding to each energy layer in the set to form the initial beam information; The energy layers are arranged according to the rack angle order. Only one energy layer is selected for each rack angle, and the energy layers selected from the rack angles from smallest to largest are selected from highest to lowest. If there is no energy layer in the initial energy layer set corresponding to a certain rack angle that is lower than the energy layer already selected in the adjacent previous rack angle, then the highest energy layer is selected from the initial energy layer set, and the selection operation of selecting energy layers from the rack angles from smallest to largest in descending order is repeated. Based on the energy layer corresponding to each rack angle and the beam spot information corresponding to that energy layer, the dose contribution relationship of each beam spot to each voxel in the target area is calculated, and the dose influence matrix is ​​obtained.

[0045] As a preferred embodiment, the method for determining the beam information in the initial radiotherapy plan and the dose influence matrix corresponding to that beam information is as follows: Generate the energy layer and its beam spot information corresponding to each preset rack angle to form the initial beam information; Based on the initial beam information and treatment equipment parameters, determine the water equivalent path length vector corresponding to each gantry angle, such as... Figure 2 As shown, each element represents the water equivalent path length under the corresponding energy layer. The water equivalent path length vectors corresponding to different angles are completed to form vectors of the same dimension. The gantry angles are clustered by the deviation between vectors, and the gantry angles whose deviation between vectors meets the preset threshold are divided into the same angle group. Based on the clustering results, the particle swarm optimization algorithm is used to optimize and organize the energy layers according to the order of the gantry angles. According to the energy layer corresponding to each gantry angle and its corresponding beam spot information, the dose contribution relationship of each beam spot to each voxel in the target area is calculated to obtain the dose influence matrix.

[0046] In the specific implementation of the above preferred method, treatment plan parameters, such as dosage prescriptions, are established based on the patient's CT images, target area, and information on the structures of organs at risk, for the subsequent construction of the objective function; additionally, a set of gantry angles is set.

[0047] in, This represents the total number of control points.

[0048] For each rack angle The corresponding energy layer and beam spot information are generated, denoted as:

[0049] in, For frame angle The corresponding number of energy layers.

[0050] Based on the initial beam information (including gantry angles and their corresponding energy layers and beam spot information) and the treatment equipment parameters (known), calculate the water equivalent path length corresponding to each gantry angle:

[0051] in, Indicates frame angle The corresponding water equivalent path length is a vector; The sequence is composed of the water equivalent path lengths corresponding to all rack angles:

[0052] Based on the water equivalent path length characteristics corresponding to each rack angle, a clustering method is used to group the rack angles, and control points with similar water equivalent path lengths (deviations meeting a preset threshold) are grouped into the same angle group:

[0053] in, This refers to the number of frame angle groups.

[0054] Angular clustering can characterize the similarity of proton beam penetration paths at different gantry angles, providing a basis for subsequent beam information processing and dose calculation.

[0055] For the control points in each angle group, based on the clustering results, the energy layers and beam spots are optimized and organized according to the rack angle order to form an energy layer sequence under continuous rack angles:

[0056] in, Indicates frame angle The corresponding set of energy layers.

[0057] Simultaneously organize the corresponding beam pattern set:

[0058] in, Indicates frame angle The set of clusters of spots below.

[0059] Based on the processed energy layer and beam spot information, the beam information is regenerated (forming new beam information), and the dose contribution relationship of each beam spot to the patient voxel is calculated to obtain the dose influence matrix: .in, Indicates the first The first spot size to the first Dosage contribution of individual factors.

[0060] Further establish dose mapping relationships: ;in, Let be the dose distribution vector. This is the beam spot weight vector.

[0061] An initial radiotherapy plan is formed based on the new beam information, dose influence matrix, and beam spot weight vector (which can be an initialized all-1 vector), providing a basis for subsequent beam spot weight distribution optimization and control point delivery time distribution calculation.

[0062] This preferred method optimizes the energy layer and beam organization based on the water equivalent path length characteristics, thereby improving the utilization efficiency of beam information.

[0063] The two preferred methods described above can comprehensively consider energy layer distribution, beam spot distribution, and treatment delivery efficiency. Through energy layer optimization, beam spot distribution optimization, and time characteristic evaluation, a more balanced beam spot weight distribution and control angle time distribution can be achieved, thereby improving treatment transmission efficiency, reducing transmission time fluctuations between gantry angles, and providing technical support for the efficient implementation of proton arc radiotherapy.

[0064] During proton arc therapy, the number of beam spots and the distribution of beam spot weights differ significantly at different gantry angles, leading to uneven delivery times at each angle, fluctuations in gantry rotation speed, and a decrease in overall treatment efficiency. To achieve uniform or near-uniform gantry delivery, while ensuring the dose distribution meets clinical prescription requirements (the actual dose distribution Ax, which is made to approximate b), an angle delivery time balance constraint is further introduced. Preferably, the following optimization model is constructed:

[0065]

[0066] In the formula, This represents the dose effect matrix; This represents the prescription dose vector composed of the doses at each voxel in the target region; Represents the beam spot weight vector; Represents the 0-1 grouping matrix of the beam spot; Represents a centered matrix; Indicates the total number of beam spots; Indicates the weight of the sparse constraint on the beam pattern; This represents the weight of the angle delivery time balance term, i.e., the weight of the dispersion term; This represents the minimum number of hops threshold.

[0067] Dose effect matrix Each element of the matrix represents the dose contribution of a specific beam spot to a specific voxel in the patient's body under unit weight. The cumulative dose of each voxel in the patient's body under the combined effect of all beam spots is obtained by multiplying the dose influence matrix with the beam spot weight vector. The prescribed dose vector is denoted as... (A vector consisting of the doses at each voxel in the target region, known); the beam spot weight vector is denoted as... The total number of spots is recorded as (Known quantity).

[0068] The first term in the optimization model is the dose fidelity term, used to ensure that the optimized dose distribution meets clinical prescription requirements. The second term is the beam sparsity constraint term, used to compress beam spot weights and reduce the number of beam spots involved in irradiation; using the L1 norm can improve efficiency. The third term is the angle delivery time equalization term, which uses the L2 norm for convenient calculation. Since the delivery time corresponding to different angles is approximately proportional to the total beam spot weight at that angle, minimizing the dispersion between the total beam spot weights at different angles can reduce the delivery time fluctuations between angles, achieving balanced control of the gantry rotation speed.

[0069] Simultaneously, a minimum hop count constraint is introduced: That is, for any speckle, when its weight is zero, it means that the speckle has been deleted; when its weight is greater than zero, it must satisfy the following condition: .

[0070] As a preferred implementation method, the beam spot weight vector can be optimized using the alternating direction multiplier method, and the implementation method is as follows: Optimize the following issues:

[0071]

[0072] In the formula, This represents the dose effect matrix; This represents the prescription dose vector composed of the doses at each voxel in the target region; Represents the beam spot weight vector; Represents the 0-1 grouping matrix of the beam spot; Represents a centered matrix; Indicates the total number of beam spots; Indicates the weight of the sparse constraint on the beam pattern; This represents the weight of the angle delivery time balance term, i.e., the weight of the dispersion term; This represents the minimum number of hops threshold; Rewritten as:

[0073]

[0074] Among them, auxiliary variables Used for separation Sparse terms, auxiliary variables Used to separate minimum hop count constraints;

[0075] The optimal beam pattern distribution is obtained by solving the alternating direction multiplier method: initialization parameters. Lagrange multipliers The superscript indicates initialization, with a preset maximum iteration parameter. Convergence tolerance and penalty parameters Then iterative updates are performed, at the [number]th [time]. In the next iteration, by solving Sub-problems Subproblems and Sub-problems updated ,like Let x be the L1 norm, then:

[0076]

[0077]

[0078] in, The subproblem is used to update the beam spot weight variable. The subproblem is used to handle the beam speckle sparsity constraint term. The subproblem is used to handle the minimum hop count constraint. The solution to the subproblem is a quadratic convex optimization problem, and its closed-form solution is:

[0079] The subproblem can be solved using the soft thresholding operator, i.e.:

[0080] make ,but The subproblem can be solved using the projection operator, and its update form is:

[0081] Then update the Lagrange multipliers:

[0082]

[0083] Then check the residuals and the number of iterations to determine if convergence has occurred. If the conditions are met... Or reach the maximum number of iterations If the iteration fails, the iteration stops; otherwise, the next iteration continues.

[0084] By using the alternating direction multiplier method described above, the optimal beam distribution can be optimized while ensuring the dosimetric target. It also satisfies the beam sparsity constraint, the minimum number of jumps constraint, and the equalization constraint of delivery time at each angle, thereby reducing the delivery time fluctuation between different gantry angles and improving the delivery efficiency and gantry rotation speed stability of proton arc therapy planning.

[0085] Specifically, in When x is the L1 norm, the beam spot distribution and delivery time are efficiently optimized by the joint modeling of the least squares data fitting term (i.e., the L2 norm term) and the global L1 regularization term (i.e., the L1 norm term), thereby obtaining an arc radiotherapy plan with uniform beam delivery time, i.e., the beam spot weight vector.

[0086] Regarding the sparseness of the beam pattern distribution, ideally, it can be achieved by using... The norm is directly used as a sparse regularization term in the objective function. However, Norms are non-convex and discontinuous, and direct solutions often involve combinatorial optimization. When proton arc therapy plans contain a large number of candidate beam spots, the computational complexity is high, making it difficult to obtain efficient and stable solutions. Although Norms can be introduced into the ADMM solution framework, but the corresponding subproblems typically employ hard thresholding operations, and the speckle weights may jump between zero and non-zero states, making the algorithm quite sensitive to initial values ​​and penalty parameters, and only guaranteeing convergence to a local solution. Therefore, as a preferred embodiment, this method adopts... Norm as The convex relaxation form of the norm, i.e. Norms can compress some low-weight beam patterns to zero while preserving the convexity and solvability of the optimization problem, thereby promoting sparsity in the beam pattern weight distribution. In the ADMM solution process, The subproblems corresponding to the norm have soft-threshold closed-form solutions, which have the advantages of simple computation, stable iteration, and clear convergence properties.

[0087] Furthermore, due to limitations in the beam output accuracy and control capabilities of proton therapy equipment, actual machines typically have a minimum number of deliverable hops. To reflect these machine characteristics, constraints are imposed on the beam spot weights:

[0088] in, This represents the minimum number of hops allowed for delivery by the actual treatment device. This constraint stipulates that the plaque can only be in a non-delivery state, i.e. Or it is in a deliverable state and meets the following conditions. Therefore, the beam spot weight cannot be located in the interval. Internally, this constraint stems from the actual beam output characteristics and delivery accuracy requirements of the machine. If the optimization results contain values ​​lower than... If a non-zero spot is present, the corresponding spot may not be accurately delivered by the machine, resulting in a deviation between the planned dose and the actual dose delivered. Introducing this constraint ensures that all retained non-zero spots meet the minimum delivery requirements of the equipment, thereby improving the consistency between the treatment plan and the actual machine performance.

[0089] Therefore, the role of the minimum hop count constraint is to directly incorporate the physical delivery constraints of the actual machine into the treatment plan optimization process, avoid generating low-weight beams that the machine cannot execute stably, improve the machine executability and actual delivery accuracy of the optimization results, and thus achieve the invention objective of generating a treatment plan that can be accurately and stably executed by the proton therapy device.

[0090] This preferred method establishes a joint optimization model including dose fidelity, beam spot sparsity, and angle delivery time equalization terms; the joint optimization model is solved using the alternating direction multiplier method to obtain the optimal beam spot distribution that satisfies dosimetric constraints, beam spot sparsity constraints, and angle delivery time equalization constraints. Figure 3 As shown, a dose-volume histogram (DVH) for a brain tumor case is presented. The DVH indicates that the target area receives adequate and uniform dose coverage, while the dose received by organs at risk, such as the brainstem, is significantly reduced, demonstrating that the treatment plan achieves good protection of normal tissues while ensuring tumor control. Furthermore, from... Figure 4 The dose slice diagrams and dose distribution maps of the same brain tumor case shown indicate that the high-dose area is mainly concentrated in the clinical target area, achieving good dose conformity; the dose decays rapidly outside the target area, with most normal tissue receiving only a low dose, indicating that the treatment plan achieves good protection of normal tissue while ensuring target coverage. It should be noted that... Figure 3 and Figure 4 The results were obtained under the following conditions: Initial beam information in the initial radiotherapy plan and the corresponding dose influence matrix acquisition method were used: An initial energy layer set corresponding to each preset gantry angle and beam spot information corresponding to each energy layer in the set were generated to form the initial beam information; the energy layers were organized according to the gantry angle order, where only one energy layer was selected for each gantry angle, and the selected energy layers were arranged from high to low order from the smallest to the largest gantry angle. If there was no energy layer in the initial energy layer set corresponding to a certain gantry angle that was lower than the energy layer already selected for its adjacent previous gantry angle, the highest energy layer was selected from that initial energy layer set, and the selection operation was repeated from high to low order from the smallest to the largest gantry angle; based on the organized energy layer corresponding to each gantry angle and the beam spot information corresponding to that energy layer, the dose contribution relationship of each beam spot to each voxel in the target area was calculated to obtain the dose influence matrix.

[0091] To further illustrate the effectiveness of the present invention, it can be further performed based on optimizing the beam spot weight (i.e., the beam spot weight vector). Calculation and evaluation of treatment transmission time.

[0092] In practice, a corresponding proton arc therapy beam sequence is constructed based on the obtained optimal beam spot weight distribution. The beam delivery sequence during actual treatment is generated based on the beam spot information, beam spot weight, and energy layer information at each gantry angle, and the energy layer configuration, beam spot scanning path, and particle fluence information corresponding to each gantry angle are extracted.

[0093] Calculate the transmission time for each control angle. For the first... i The total transmission time for each control angle consists of the spot switch time and the spot spray time, i.e.: ,in, , indicating the first i The beam switching time caused by the deflection of the scanning magnet between all beams at each control angle; , indicating the first Particle transport time corresponding to all beam spots under each control angle; Indicates the first The set of beams contained in each control angle.

[0094] Further, the time distribution of all control angles was obtained: ,in, To control the total number of angles.

[0095] Calculate the treatment time evaluation index. Based on the transmission time corresponding to each control angle, calculate the total transmission time of the entire proton arc therapy plan: Calculate the transmission time for the maximum control angle: The transmission time required to characterize the control angle with the highest load during treatment is used. Simultaneously, to evaluate the uniformity of the time distribution among different control angles, the coefficient of time variation is calculated. ,in, Indicates the average control angle transmission time; This represents the standard deviation of the control angle transmission time.

[0096] Based on the above calculation results, the time distribution characteristics of the proton arc therapy plan are obtained, and evaluation indicators such as the transmission time of each control angle, total transmission time, maximum control angle transmission time, and time variation coefficient are output. Among them, the time variation coefficient is used to measure the consistency of transmission time among control angles during treatment; the maximum control angle transmission time is used to reflect the time requirement of the slowest control angle during gantry rotation; and the total transmission time is used to evaluate the overall treatment efficiency.

[0097] The above indicators can be used to quantitatively assess the time characteristics of the optimized treatment plan and provide a basis for subsequent gantry speed adjustment and uniform operation analysis.

[0098] Figure 5 Showing Figure 3 The results show the comparison of the temporal smoothness of the treatment plan before and after optimization for the corresponding brain tumor case. Among them, (a) is the comparison of time distribution, (b) is the comparison of absolute time deviation, (c) is the comparison of time distribution after sorting, (d) is the comparison of beam extraction (radiation) time distribution, (e) is the comparison of absolute deviation of beam extraction (radiation) time, (f) is the comparison of beam extraction (radiation) time distribution after sorting, (g) is the comparison of total time coefficient of variation, (h) is the comparison of beam extraction (radiation) time coefficient of variation, and (i) is the sparsity comparison.

[0099] Example 2 This application also relates to an electronic device, including a memory and a processor, wherein the memory stores a computer program, and the processor executes the computer program to implement the steps of the method described above.

[0100] The electronic device can be a desktop computer, laptop, handheld computer, or cloud server, etc. The processor can be a Central Processing Unit (CPU), or other general-purpose processors, digital signal processors (DSPs), application-specific integrated circuits (ASICs), field-programmable gate arrays (FPGAs), or other programmable logic devices, discrete gate or transistor logic devices, discrete hardware components, etc. The memory can be used to store computer programs and / or modules. The processor performs various functions of the electronic device by running or executing the computer programs and / or modules stored in the memory, and by accessing data stored in the memory.

[0101] The relevant technical solutions are the same as above, and will not be repeated here.

[0102] Example 3 This application also relates to a computer-readable storage medium having a computer program stored thereon, which, when executed by a processor, implements the steps of the method described above.

[0103] Specifically, the memory may include high-speed random access memory, as well as non-volatile memory, such as hard disks, RAM, plug-in hard disks, smart media cards (SMC), secure digital (SD) cards, flash cards, at least one disk storage device, flash memory device, or other volatile solid-state storage devices.

[0104] The relevant technical solutions are the same as above, and will not be repeated here.

[0105] Example 4 This application provides a computer program product or computer program that includes computer instructions stored in a computer-readable storage medium. A processor of a computer device reads the computer instructions from the computer-readable storage medium and executes the computer instructions, causing the computer device to perform the steps of the method described in the above embodiments of this application.

[0106] The relevant technical solutions are the same as above, and will not be repeated here.

[0107] In summary, the method presented in this embodiment achieves synergistic optimization of dosimetric targets, angle-time balance, and sparsity through joint optimization, reducing transmission time fluctuations between different gantry angles. By establishing a treatment time evaluation system, it enables quantitative assessment of the transmission efficiency of proton arc radiotherapy plans. The proposed method can improve treatment efficiency and gantry operational stability while ensuring dose distribution quality, providing an effective planning method for the clinical application of proton arc radiotherapy technology and contributing to its development.

[0108] Those skilled in the art will readily understand that the above description is merely a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, and improvements made within the spirit and principles of the present invention should be included within the scope of protection of the present invention.

Claims

1. A sparse proton arc radiotherapy planning method based on irradiation time homogenization, characterized in that, include: Obtain an initial radiotherapy plan, including a beam weight vector, beam information, and a dose influence matrix corresponding to the beam information. The beam information includes the energy layer corresponding to each gantry angle and the beam spot corresponding to each energy layer. The dimension of the beam weight vector is the same as the total number of known beam spots. The number of rows in the dose influence matrix is ​​the same as the total number of voxels divided in the target area, and the number of columns is the same as the total number of known beam spots. Based on the beam information, a beam spot 0-1 grouping matrix B is constructed, where the rows correspond to the frame angles and the vertical columns correspond to the beam spots. A 1 indicates that the beam spot corresponding to the frame angle participates in illumination, and a 0 indicates that the beam spot corresponding to the frame angle does not participate in illumination. Multiplying P, B, and x together yields... , representing the degree of dispersion of the total beam spot weights across different frame angles, where x represents the beam spot weight vector; P is the centering matrix, , This represents a diagonal matrix with a diagonal value of 1. Represents a vector consisting entirely of 1s; The product of the dose influence matrix and the beam spot weight vector is used as the actual dose distribution, and the beam spot distribution sparsity is constructed based on the beam spot weight vector. The beam spot weight vector is optimized by minimizing the weighted sum of the deviation between the actual dose distribution and the prescription dose distribution, the beam spot distribution sparsity, and the aforementioned degree of dispersion.

2. The sparse proton arc radiotherapy planning method as described in claim 1, characterized in that, The method for determining the beam information in the initial radiotherapy plan and the corresponding dose-effect matrix is ​​as follows: Generate an initial energy layer set corresponding to each preset rack angle and the beam spot information corresponding to each energy layer in the set to form the initial beam information; The energy layers are arranged according to the rack angle order. Only one energy layer is selected for each rack angle, and the energy layers selected from the rack angles from smallest to largest are selected from highest to lowest. If there is no energy layer in the initial energy layer set corresponding to a certain rack angle that is lower than the energy layer already selected in the adjacent previous rack angle, then the highest energy layer is selected from the initial energy layer set, and the selection operation of selecting energy layers from the rack angles from smallest to largest in descending order is repeated. Based on the energy layer corresponding to each rack angle and the beam spot information corresponding to that energy layer, the dose contribution relationship of each beam spot to each voxel in the target area is calculated, and the dose influence matrix is ​​obtained.

3. The sparse proton arc radiotherapy planning method as described in claim 1, characterized in that, The method for determining the beam information in the initial radiotherapy plan and the corresponding dose-effect matrix is ​​as follows: Generate the energy layer and its beam spot information corresponding to each preset rack angle to form the initial beam information; Based on the initial beam information and treatment equipment parameters, the water equivalent path length vector corresponding to each gantry angle is determined, where each element represents the water equivalent path length under the corresponding energy layer; the water equivalent path length vectors corresponding to different angles are completed to form vectors of the same dimension, and the gantry angles are clustered by the deviation between vectors, and the gantry angles whose deviation between vectors meets the preset threshold are divided into the same angle group. Based on the clustering results, the energy layers are optimized and organized according to the rack angle order using the particle swarm optimization algorithm; Based on the energy layer and its corresponding beam spot information corresponding to each rack angle after processing, the dose contribution relationship of each beam spot to each voxel in the target area is calculated to obtain the dose influence matrix.

4. The sparse proton arc radiotherapy planning method as described in claim 1, characterized in that, The optimization problem of the beam spot weight vector is: In the formula, This represents the dose effect matrix; This represents the prescription dose vector composed of the doses at each voxel in the target region; Represents the beam spot weight vector; Represents the 0-1 grouping matrix of the beam spot; Represents a centered matrix; Indicates the total number of beam spots; Indicates the weight of the sparse constraint on the beam pattern; Indicates taking The L1 norm characterizes the sparsity of the beam pattern distribution; This represents the weight of the angle delivery time balance term, i.e., the weight of the dispersion term; Indicates taking The L2 norm; This represents the minimum number of hops threshold.

5. The sparse proton arc radiotherapy planning method as described in claim 1, characterized in that, The beam spot weight vector is optimized using the alternating direction multiplier method, and the implementation is as follows: Optimize the problem: In the formula, This represents the dose effect matrix; This represents the prescription dose vector composed of the doses at each voxel in the target region; Represents the beam spot weight vector; This indicates the sparseness of the beam pattern distribution; Represents the 0-1 grouping matrix of the beam spot; Represents a centered matrix; Indicates the total number of beam spots; Indicates the weight of the sparse constraint on the beam pattern; This represents the weight of the angle delivery time balance term, i.e., the weight of the dispersion term; Indicates taking The L2 norm; This represents the minimum number of hops threshold; Rewritten as: Among them, auxiliary variables Used for separation Sparse terms, auxiliary variables Used to separate minimum hop count constraints; The optimal beam pattern distribution is obtained by solving the alternating direction multiplier method: initialization parameters. Lagrange multipliers The superscript indicates initialization, with a preset maximum iteration parameter. Convergence tolerance and penalty parameters Then iterative updates are performed, at the [number]th [time]. In the next iteration, the following solution is obtained: Sub-problems Subproblems and Sub-problems updated : Determine if the termination condition has been met; if so, complete the solution; otherwise, continue to the next iteration.

6. An electronic device comprising a memory and a processor, wherein the memory stores a computer program, characterized in that, When the processor executes the computer program, it implements the steps of the method as described in any one of claims 1 to 5.

7. A computer-readable storage medium, characterized in that, The computer-readable storage medium includes a stored computer program, wherein the computer program, when executed by a processor, controls the device on which the storage medium is located to perform the steps of the method as described in any one of claims 1 to 5.

8. A computer program product, comprising a computer program or instructions, characterized in that, When the computer program or instructions are executed by a processor, they implement the steps of the method as described in any one of claims 1 to 5.

Images