CN105477789A - 基于二次规划模型抑制总出束时间的动态调强放疗方法 - Google Patents

Abstract

本发明公开了一种基于二次规划模型抑制总出束时间的动态调强放疗优化方法。该方法首先依据临床剂量目标和约束进行建模；然后通过前驱、后驱叶片移动矩阵对射野强度图进行数学表述，并在此基础上建立了带有上下叶片碰撞约束和典型凸凹槽效应约束的凸规划条件；通过分析叶片运动速度及总出束时间与叶片运动矩阵的关系，建立了相应的凸规划条件；最后针对叶片单向运动调制剂量的局限性和叶片随机往复运动对叶片产生磨损的缺点，建立了叶片双向运动模式。本发明成功地得到了一种能够满足绝大多数临床实施软硬件方面约束和需求的动态调强放疗优化方法。如果将照射野角度提升至一定数量后，可以直接用于容积调强放疗的数学建模和求解中。

Description

基于二次规划模型抑制总出束时间的动态调强放疗方法

技术领域

本发明属于放射治疗数学规划技术领域，涉及调强放疗逆向计划求取方法，特别是涉及一种用于多叶准直器实现动态调强放疗逆向计划的优化方法。

背景技术

目前，恶性肿瘤的最主要的治疗手段包括外科手术治疗、化学类治疗和放射类治疗，尤其对于中晚期癌症病例，放射类治疗所占比例较大；而在放射类治疗手段中，近二十年国内外应用最为广泛的一种治疗方法为调强放疗，即在适形放疗中加入强度可调这个维度。

全球范围中主要的调强放疗设备生产厂商包括瓦里安(Varian)、医科达(Elekta)、菲利普(Philips)、西门子(Siemens)等公司。这些厂商推出的调强放疗软硬件系统基本上垄断了全球在该领域的市场份额。尤其是瓦里安和医科达，由于采用了先进的逆向计划系统，它们总是走在调强放疗技术的最前沿，如最新的容积调强治疗计划系统(VMAT，IMAT)。

调强放疗逆向计划系统中最关键的一个步骤就是逆向计划的求取。逆向计划的输入包括临床医生对靶区的剂量需求及对正常组织和危险区域待保护组织的剂量约束，以及所采用射线输出设备的硬件约束等。而逆向技术的输出为可供放疗硬件设备执行的子野位置、大小和照射机器跳数(出束时间)。由于逆向计划求取所建立数学模型的超大规模和复杂程度，目前国内外知名厂商大多使用启发式算法来实现整个逆向优化过程。如瓦里安的治疗计划系统采用直接子野优化算法(DirectApertureOptimization，DAO)以实现多叶准直器叶片位置的优化计算和射野强度的计算；在每一次叶片位置改变后，算法会自动确认计量覆盖约束是否符合要求，一旦违背约束条件算法将自动重新优化计算；医科达则采用模拟退火优化算法(SimulatedAnnealingOptimization，SAO)直接优化多叶准直器位置，并得到序列的强度图，在产生预定目标剂量失败时需要进一步采样并计算；菲利普则首先设置照射弧参数，然后进行调强放疗逆向优化得到强度图，接着考虑滑窗参数进行子野分割；以得到的子野分割结果作为初始参数输入优化模型中，并结合机器参数约束(包括叶片移动、计量率和机架速度等)进行二次优化得到优化权重系统，从而得到优化后的逆向计划。启发式算法虽然能够很快收敛到某个结果，但是却存在所得结果不确定不唯一、可能存在更优结果等缺点。与启发式算法不同，确定型算法所得到的优化结果唯一，且在所涉及数学模型中为全局最优解；然而确定型算法所适用的数学模型非常有限，目前仅凸规划类数学模型中的几种数学模型有较为成熟的求解算法。而对于上述的逆向计划求取过程中所需要考虑的各种需求、约束和变量，想建立一个比较理想的凸规划数学模型几乎是一件不可能完成的事情。

本发明从建立凸规划模型着手，同时考虑各种临床需求和主要的多叶准直器硬件需求和约束，以动态多叶准直器输出调强射线为治疗手段，最终得到全局最优的动态调强放疗逆向计划。其中，临床需求包括医生对靶区的最低剂量需求，以及对各种不同种类非靶区组织的最高剂量约束、平均剂量约束。多叶准直器硬件需求和约束包括对准直器射线总输出时间的约束，叶片凸凹槽效应所导致的欠剂量约束，上下相邻叶片运动碰撞约束，叶片最大运动速度约束等。

基于该发明所实现的动态调强放疗逆向计划系统很容易在现有的瓦里安和医科达多叶准直器硬件平台上得到应用。

发明内容

考虑到传统非凸数学规划模型在动态调强放疗中应用的诸多问题，如求解结果不唯一、所得解非全局最优。本发明提供了一种基于二次规划模型抑制总出束时间的动态调强放疗方法。

为实现上述目的，本发明采用的技术方案：该基于二次规划模型抑制总出束时间的动态调强放疗方法，其特征在于按下述步骤构建：

(1)首先依据传统的二次规划数学模型对临床剂量目标和约束进行建模，二次规划模型表达式为：

min{x^TQx+c·x+e}

$S.t.\left\{\begin{array}{c}{d}^{TP}\≤F^T\·x\\ F^N\·x\≤d^{NU}\\ Average(F^N\·x)\≤d^A\\ 0\≤x\end{array}\right.$

其中，Q＝Σp^T(F^T)'·F^T+Σp^N(F^N)'·F^N，c＝-2Σp^T(d^TP)'·F^T，e＝Σp^T(d^TP)'·d^TP，d^TP为临床对目标靶的理想杀灭剂量，d^TL为最低剂量约束，p^T为相应的重要性权重系数；d^NU为对非靶区器官的最大剂量约束，d^A为均值剂量约束，p^N为相应的重要性权重系数。为定照射野离散化照射强度图i的变量矩阵，向量x为所有离散化强度图变量，F为给定照射的放射线对体模照射的效应剂量传递矩阵；

(2)然后在通过前驱、后驱叶片移动矩阵对射野强度图进行数学表述，其中，基于叶片移动矩阵带有相邻叶片碰撞约束条件的凸规划约束条件如下：

$\stackrel{\‾}{X}\·W=A^{+}-A^{-}$

0≤Α⁺

0≤Α^-

I^T＝A⁺·W^-1

I^L＝A^-·W^-1

I^L(j,k)≤I^T(j+1,k)

I^L(j+1,k)≤I^T(j,k)

其中，叶片移动矩阵为两个整数矩阵I^L和I^T，来分别描述对于不同位置的射束块的前驱叶片的叶片打开时间和后驱叶片的叶片关闭时间，矩阵I^L中的第i行第j列的元素代表了第i行前驱叶片运动至第j列的位置时的打开时刻，而矩阵I^T中的第i行第j列的元素代表了第i行后驱叶片运动至第j列的位置时的关闭时刻；

(3)在此基础上建立了带有上下叶片碰撞约束和典型凸凹槽效应约束的凸规划条件满足凸凹槽效应约束的约束条件如下：

I^T(j,k)+I^L(j,k)＝I^T(j+1,k)+I^L(j+1,k)；

(4)接着通过分析叶片运动速度与叶片运动矩阵的关系以及总出束时间与叶片运动矩阵的关系，建立带有最大叶片运动速度约束、总出束时间约束的凸规划条件；

叶片运动速度约束条件可以描述如下：

I^T(j,k+1)-I^T(j,k)≥τ₀，I^L(j,k+1)-I^L(j,k)≥τ₀

总出束时间约束条件，描述如下：

$\underset{i=1}{\overset{l}{\Σ}}\left\{\underset{j}{\max }\right\{I_i^T(j,n)\left\}\right\}\≤{\mathrm{TNMU}}_0$

(5)最后根据叶片单向运动的局限性和叶片随机往复运动对叶片的磨损，建立了叶片双向运动模式，在每个照射野上，前驱和后驱叶片首先从最左端出发，全部运动至最右端形成第一次调制，然后叶片从最右端出发，全部运动至最左端形成第二次调制，双向运动调制完毕后即完成该射野剂量调制，完毕后，照射头运动至下一个照射野方向继续进行下一轮的调制；叶片双向运动模式约束可以将强图分布图分为两部分组成，第一部分为第一次调制结果，第二部分为第二次调制结果，由此得到相应的约束条件：

$\stackrel{\‾}{X}=\stackrel{\‾}{Y}+\stackrel{\‾}{Z}$

$\stackrel{\‾}{Y}\≥0$

$\stackrel{\‾}{Z}\≥0$

(6)建立能够满足绝大多数动态调强放疗在软硬件方面约束和需求的数学模型，得到了总的凸规划数学模型，描述如下：

$min\left\{{\mathrm{\Σp}}^T\right||d^T-d^{TP}|{|}_2^2+{\mathrm{\Σp}}^N\left|\right|d^N-0\left|{|}_2^2\right\}$

其中，x为展开的向量表达形式， $W={\left[\begin{array}{ccccc}1& -1& 0& & 0\\ 0& 1& -1& L& 0\\ 0& 0& 1& & M\\ & M& & O& -1\\ 0& 0& L& 0& 1\end{array}\right]}_{n\×n},$ $W^{-1}={\left[\begin{array}{ccccc}1& 1& 1& & 1\\ 0& 1& 1& L& 1\\ 0& 0& 1& & M\\ & M& & O& 1\\ 0& 0& L& 0& 1\end{array}\right]}_{n\×n},V={\left[\begin{array}{ccccc}1& 0& 0& & 0\\ -1& 1& 0& L& 0\\ 0& -1& 0& & M\\ & M& & O& 0\\ 0& 0& L& -1& 1\end{array}\right]}_{n\×n},V^{-1}={\left[\begin{array}{ccccc}1& 0& 0& & 0\\ 1& 1& 0& L& 0\\ 1& 1& 1& & M\\ & M& & O& 0\\ 1& 1& L& 1& 1\end{array}\right]}_{n\×n}.$

上述数学模式展开即为线性约束二次规划模型，通过内点法或共轭梯度法凸规划求解方法可以快速得到全局最优解，求解得到的I^T，I^L以及J^T，J^L即为动态调强逆向计划中的多叶准直器运动矩阵，可以直接写入驱动文件对叶片进行驱动控制，完成相应的剂量调制过程。

采用上述技术方案的有益效果：该基于二次规划模型，通过线性约束二次规划的凸规划数学模型得到全局最优的叶片驱动计划，该凸规划数学模型能够以给定的总出束时间为约束条件获取临床剂量需求的最优分布结果，能够对多叶准直器中的上下叶片运动碰撞进行严格约束，能够消除多叶准直器在剂量调制中所产生的典型的凸凹槽效应，能够按照多叶准直器的叶片最大运动速度对叶片驱动计划进行约束，能够让叶片在单照射野中进行双向运动，从而让叶片回至起始位置，克服随机往复运动所导致的叶片磨损和定位不准确的问题；该方法可以直接升级至容积调强放疗优化逆向计划获取中。既能满足多叶准直器凸凹槽效应约束、上下叶片碰撞约束及叶片最大运动速度要求约束，又能够最大程度抑制总出束时间的具备双向运动模式的动态调强放疗优化方法，该方法可以借助于凸规划内点法进行快速求解，从而获取全局最优解。

附图说明

图1为基于多叶准直器的动态调强放疗整体流程图。

具体实施方式

如图1所示，基于多叶准直器的动态调强放疗逆向计划的整体工作流程主要包括待治疗区域的三维医学影像的获取及靶区勾勒、待治疗组织区域放射剂量的模拟计算、基于多叶准直器的动态调强放疗逆向优化、最优化结果的仿真评估、叶片驱动文件输出并执行剂量调制等。

由于采用了凸规划模型，优化过程的算法是确定型的，因此本发明中所涉及的动态调强放疗逆向优化的最重要步骤为数学模型的建立，即确定模型的输入。输入条件包括靶区组织和非靶区组织的勾勒确定。在获取了所需区域的三维医学影像后(包括CT成像、PET成像等)，临床医生就可以进行GTV勾勒和修改，并在此基础上对组织进行离散化。临床医生会讨论并确定目标靶区的处方剂量以及待保护组织的剂量约束。在确定待治疗靶区特性后，临床医生可以依据经验设定照射野的数量及相应的方向。与此同时，还必须通过蒙特卡洛等模拟计算的方式得到照射—剂量输出矩阵。另一方面，还需要确定所采用的多叶准直器的诸多物理特性，如叶片数量、相应的宽度、运动控制精度、最大运动速度等等。临床医生会依据经验给出总出束时间限定条件。在获取上述输入后，依据本发明所开发的逆向治疗系统将生成相应的输入数据矩阵及向量；在此基础上，采用内点法等确定型求解算法能够很快求得全局最优的叶片驱动文件(若存在最优解)。若输入参数要求过于苛刻，则可能导致输出无解，此时系统将提示无法满足的输入条件；一旦获取最优解，系统将依据该输出文件对所获取的放疗计划的仿真结果进行三维可视化显示，并对该计划优劣进行客观评估，包括描绘计算剂量体积曲线、等剂量曲线图等。评估结果将直接供临床医生进行参考，医生可能对相应的输入参数进行修改并重新输入系统中进行优化输出。优化—评估的循环过程将依据医生对评估结果的分析重复多次直至满意结果生成为止。一旦该循环结束，所获得的最优叶片驱动文件将直接通过上位机对多叶准直器叶片进行控制，并完成相应的剂量调制过程。

基于二次规划模型抑制总出束时间的动态调强放疗方法，按下述步骤构建：

(1)首先依据传统的二次规划数学模型对临床剂量目标和约束进行建模，二次规划模型表达式为：

min{x^TQx+c·x+e}

$S.t.\left\{\begin{array}{c}{d}^{TP}\≤F^T\·x\\ F^N\·x\≤d^{NU}\\ Average(F^N\·x)\≤d^A\\ 0\≤x\end{array}\right.$

其中，Q＝Σp^T(F^T)'·F^T+Σp^N(F^N)'·F^N，c＝-2Σp^T(d^TP)'·F^T，e＝Σp^T(d^TP)'·d^TP，d^TP为临床对目标靶的理想杀灭剂量，d^TL为最低剂量约束，p^T为相应的重要性权重系数；d^NU为对非靶区器官的最大剂量约束，d^A为均值剂量约束，p^N为相应的重要性权重系数。为定照射野离散化照射强度图i的变量矩阵，向量x为所有离散化强度图变量，F为给定照射的放射线对体模照射的效应剂量传递矩阵；

(2)然后在通过前驱、后驱叶片移动矩阵对射野强度图进行数学表述，其中，基于叶片移动矩阵带有相邻叶片碰撞约束条件的凸规划约束条件如下：

$\stackrel{\‾}{X}\·W=A^{+}-A^{-}$

0≤Α⁺

0≤Α^-

I^T＝A⁺·W^-1

I^L＝A^-·W^-1

I^L(j,k)≤I^T(j+1,k)

I^L(j+1,k)≤I^T(j,k)

其中，叶片移动矩阵为两个整数矩阵I^L和I^T，来分别描述对于不同位置的射束块的前驱叶片的叶片打开时间和后驱叶片的叶片关闭时间，矩阵I^L中的第i行第j列的元素代表了第i行前驱叶片运动至第j列的位置时的打开时刻，而矩阵I^T中的第i行第j列的元素代表了第i行后驱叶片运动至第j列的位置时的关闭时刻；

(3)在此基础上建立了带有上下叶片碰撞约束和典型凸凹槽效应约束的凸规划条件满足凸凹槽效应约束的约束条件如下：

I^T(j,k)+I^L(j,k)＝I^T(j+1,k)+I^L(j+1,k)；

(4)接着通过分析叶片运动速度与叶片运动矩阵的关系以及总出束时间与叶片运动矩阵的关系，建立带有最大叶片运动速度约束、总出束时间约束的凸规划条件；

叶片运动速度约束条件可以描述如下：

I^T(j,k+1)-I^T(j,k)≥τ₀，I^L(j,k+1)-I^L(j,k)≥τ₀

总出束时间约束条件，描述如下：

$\underset{i=1}{\overset{l}{\Σ}}\left\{\underset{j}{\max }\right\{I_i^T(j,n)\left\}\right\}\≤{\mathrm{TNMU}}_0$

最后根据叶片单向运动的局限性和叶片随机往复运动对叶片的磨损，建立了叶片双向运动模式，在每个照射野上，前驱和后驱叶片首先从最左端出发，全部运动至最右端形成第一次调制，然后叶片从最右端出发，全部运动至最左端形成第二次调制，双向运动调制完毕后即完成该射野剂量调制，完毕后，照射头运动至下一个照射野方向继续进行下一轮的调制；叶片双向运动模式约束可以将强图分布图分为两部分组成，第一部分为第一次调制结果，第二部分为第二次调制结果，由此得到相应的约束条件：

$\stackrel{\‾}{X}=\stackrel{\‾}{Y}+\stackrel{\‾}{Z}$

$\stackrel{\‾}{Y}\≥0$

$\stackrel{\‾}{Z}\≥0$

(6)建立能够满足绝大多数动态调强放疗在软硬件方面约束和需求的数学模型，得到了总的凸规划数学模型，描述如下：

$min\left\{{\mathrm{\Σp}}^T\right||d^T-d^{TP}|{|}_2^2+{\mathrm{\Σp}}^N\left|\right|d^N-0\left|{|}_2^2\right\}$

其中，x为展开的向量表达形式， $W={\left[\begin{array}{ccccc}1& -1& 0& & 0\\ 0& 1& -1& L& 0\\ 0& 0& 1& & M\\ & M& & O& -1\\ 0& 0& L& 0& 1\end{array}\right]}_{n\×n},$ $W^{-1}={\left[\begin{array}{ccccc}1& 1& 1& & 1\\ 0& 1& 1& L& 1\\ 0& 0& 1& & M\\ & M& & O& 1\\ 0& 0& L& 0& 1\end{array}\right]}_{n\×n},V={\left[\begin{array}{ccccc}1& 0& 0& & 0\\ -1& 1& 0& L& 0\\ 0& -1& 0& & M\\ & M& & O& 0\\ 0& 0& L& -1& 1\end{array}\right]}_{n\×n},V^{-1}={\left[\begin{array}{ccccc}1& 0& 0& & 0\\ 1& 1& 0& L& 0\\ 1& 1& 1& & M\\ & M& & O& 0\\ 1& 1& L& 1& 1\end{array}\right]}_{n\×n}.$

上述数学模式展开即为线性约束二次规划模型，通过内点法或共轭梯度法凸规划求解方法可以快速得到全局最优解，求解得到的I^T，I^L以及J^T，J^L即为动态调强逆向计划中的多叶准直器运动矩阵，可以直接写入驱动文件对叶片进行驱动控制，完成相应的剂量调制过程。

Claims (1)

1.一种基于二次规划模型抑制总出束时间的动态调强放疗方法，其特征在于按下述步骤构建：

(1)首先依据传统的二次规划数学模型对临床剂量目标和约束进行建模，二次规划模型表达式为：

min{x^TQx+c·x+e}

$S.t.\left\{\begin{array}{c}{d}^{TP}\≤F^T\·x\\ F^N\·x\≤d^{NU}\\ Average(F^N\·x)\≤d^A\\ 0\≤x\end{array}\right.$

其中，Q＝Σp^T(F^T)'·F^T+Σp^N(F^N)'·F^N，c＝-2Σp^T(d^TP)'·F^T，e＝Σp^T(d^TP)'·d^TP，d^TP为临床对目标靶的理想杀灭剂量，d^TL为最低剂量约束，p^T为相应的重要性权重系数；d^NU为对非靶区器官的最大剂量约束，d^A为均值剂量约束，p^N为相应的重要性权重系数。为定照射野离散化照射强度图i的变量矩阵，向量x为所有离散化强度图变量，F为给定照射的放射线对体模照射的效应剂量传递矩阵；

(2)然后在通过前驱、后驱叶片移动矩阵对射野强度图进行数学表述，其中，基于叶片移动矩阵带有相邻叶片碰撞约束条件的凸规划约束条件如下：

$\stackrel{\‾}{X}\·W=A^{+}-A^{-}$

0≤Α⁺

0≤Α^-

I^T＝A⁺·W^-1

I^L＝A^-·W^-1

I^L(j,k)≤I^T(j+1,k)

I^L(j+1,k)≤I^T(j,k)

其中，叶片移动矩阵为两个整数矩阵I^L和I^T，来分别描述对于不同位置的射束块的前驱叶片的叶片打开时间和后驱叶片的叶片关闭时间，矩阵I^L中的第i行第j列的元素代表了第i行前驱叶片运动至第j列的位置时的打开时刻，而矩阵I^T中的第i行第j列的元素代表了第i行后驱叶片运动至第j列的位置时的关闭时刻；

(3)在此基础上建立了带有上下叶片碰撞约束和典型凸凹槽效应约束的凸规划条件满足凸凹槽效应约束的约束条件如下：

I^T(j,k)+I^L(j,k)＝I^T(j+1,k)+I^L(j+1,k)；

(4)接着通过分析叶片运动速度与叶片运动矩阵的关系以及总出束时间与叶片运动矩阵的关系，建立带有最大叶片运动速度约束、总出束时间约束的凸规划条件；

叶片运动速度约束条件可以描述如下：

I^T(j,k+1)-I^T(j,k)≥τ₀，I^L(j,k+1)-I^L(j,k)≥τ₀

总出束时间约束条件，描述如下：

$\underset{i=1}{\overset{l}{\Σ}}\left\{\underset{j}{\max }\left\{I_i^T\left(j,n\right)\right\}\right\}\≤{\mathrm{TNMU}}_0$

(5)最后根据叶片单向运动的局限性和叶片随机往复运动对叶片的磨损，建立了叶片双向运动模式，在每个照射野上，前驱和后驱叶片首先从最左端出发，全部运动至最右端形成第一次调制，然后叶片从最右端出发，全部运动至最左端形成第二次调制，双向运动调制完毕后即完成该射野剂量调制，完毕后，照射头运动至下一个照射野方向继续进行下一轮的调制；叶片双向运动模式约束可以将强图分布图分为两部分组成，第一部分为第一次调制结果，第二部分为第二次调制结果，由此得到相应的约束条件：

$\stackrel{\‾}{X}=\stackrel{\‾}{Y}+\stackrel{\‾}{Z}$

$\stackrel{\‾}{Y}\≥0$

$\stackrel{\‾}{Z}\≥0$

(6)建立能够满足绝大多数动态调强放疗在软硬件方面约束和需求的数学模型，得到了总的凸规划数学模型，描述如下：

$min\left\{{\mathrm{\Σp}}^T\right||d^T-d^{TP}|{|}_2^2+{\mathrm{\Σp}}^N\left|\right|d^N-0\left|{|}_2^2\right\}$

其中，x为展开的向量表达形式，

上述数学模式展开即为线性约束二次规划模型，通过内点法或共轭梯度法凸规划求解方法可以快速得到全局最优解，求解得到的I^T，I^L以及J^T，J^L即为动态调强逆向计划中的多叶准直器运动矩阵，可以直接写入驱动文件对叶片进行驱动控制，完成相应的剂量调制过程。