本文基于点核叠加的剂量计算模型,提出了一种由点核叠加构建的笔形束核进行剂量计算的方法,这种方法可大大提高调强放射治疗(IMRT)优化迭代过程的剂量计算速度,使得基于点核叠加技术的计划系统得以集成直接孔隙的逆向调强技术。通过对模体及临床实际病例进行试验,结果表明此方法在提高剂量计算速度的同时,也保证了迭代过程剂量计算的精度,从而保证了与使用精确剂量计算模型得到的优化结果相一致,完全可用于调强优化过程中的剂量计算。在基于点核叠加剂量计算模型的计划系统中使用该方法,可以避免重新研发笔形束剂量计算模型以及产品维护造成的成本大量增加。
引用本文： 邹豪, 王玉. 基于点核构建的笔形束模型在逆向调强中的应用. 生物医学工程学杂志, 2014, 31(1): 166-171. doi: 10.7507/1001-5515.20140033 复制
引言
调强放射治疗(intensity modified radiation therapy,IMRT)[1]是通过调整多个照射野内射束的强度分布，在靶区得到高度适形剂量的同时使周围正常组织尽可能受到较少的照射，从而达到提高治疗增益比的目的。IMRT对提高肿瘤的局部控制率以及降低正常组织并发症发生率起到了非常大的作用，因此临床应用也越来越多。
调强优化的算法根据其优化求解方式，基本可以分为确定算法和随机算法两类。确定算法主要包括梯度法[2-3] 、二次规划、线性规划、迭代法等。随机算法包括模拟退火[4]和遗传算法[5]等。每种方法都有各自的优点和缺点：梯度法优化速度快，但容易陷入局部极小点；模拟退火和遗传算法是全局优化算法，可以避免陷入局部极小点，但是迭代次数多，优化速度慢。
IMRT从实现方式上可分为两大类：一类是连续或螺旋式的断层放疗，此类又可分为将射束准直为窄条形或用机架旋转进行射束调强两种，此类调强也叫扇束调强；第二大类是基于多叶准直器(multileaf collimator,MLC)的射束调强，此类调强也叫锥束调强。锥束调强可分为静态和动态两种，其区别是出束时叶片是否在伴随运动。在静态锥束调强中，调强束是用很多MLC 的均匀束分段迭加实现的，所以也叫“停-射”(Step and Shoot)模式；
传统的调强计划使用的是两步法，即首先对每个射野进行强度分布优化，得到射野强度分布图，然后将强度分布图转换为子野序列。两步法优点是可以支持动态锥束调强，缺点是靶区适形度不好。而且对于形状不规则的靶区，如鼻咽癌等，两步法优化得到的子野数及机器跳数(monitor units,MU)非常多，所需治疗时间很长。而直接机器参数优化法(direct machine parameter optimization,DMPO)通过直接计算子野的形状及权重来得到一个可实施的调强治疗计划。它的特点是物理师可以预先设置计划的子野数，从而有效地减少子野数量，减少治疗时间，提高加速器的工作效率。同时由于MLC叶片位置可直接优化得到，避免了子野分割造成的适形度误差。因此，近年来基于随机算法的直接子野优化方法越来越受到重视，其中尤以基于模拟退火优化算法的直接孔隙优化(direct aperture optimization,DAO)方法为代表。但由于随机算法迭代次数很多，而每次迭代都需要更新剂量分布，因此提高迭代过程的剂量计算速度对优化时间影响最大。在迭代过程中，为了提高剂量计算速度，一般会采用笔形束剂量计算模型，但由于基于点核叠加的剂量计算模型具有更高的精度，因此研究如何基于点核叠加构建笔形束进行优化迭代过程的剂量计算有非常大的价值，在保证系统兼容性的前提下，有望提高迭代过程中剂量计算的精度和计算速度。
本文基于已有的点核叠加模型的计划系统，在调强优化迭代过程，提出了一种由点核叠加构建的笔形束进行剂量计算的方法，在提高了剂量计算速度的同时，也保证了迭代剂量计算的算法精度及其与精确计算结果的一致性，具体见2.3节，使得DMPO在基于点核叠加剂量计算模型系统中得以实现，避免了重新研发笔形束剂量计算模型的大量工作。
1 直接孔径优化
直接孔径调强是一步调强法，利用模拟退火算法直接优化子野形状和对应照射强度，从而得到计划所需的子野信息，这种方法比传统的两步法减少了子野数和MU，提高了实施照射的效率。
直接孔径调强主要有三方面工作：① 模拟退火优化方法；② 剂量计算；③ 目标函数与约束。
1.1 模拟退火优化方法
模拟退火算法是20世纪80年代发展起来的一种用于求解大规模优化问题的随机搜索算法，此问题的自变量组合状态空间是MLC叶片的位置和此子野的强度。优化的目标函数是处方剂量与模体剂量之间的差值函数。设定初始温度、子野叶片初始状态及权重，自变量状态改变时计算模体的剂量分布，然后得到目标函数值，与上个状态的目标函数值进行比较，如果目标函数值减小，接受这个状态，否则以概率P接受。进行一定次数的状态变换后，降低温度，重复上述步骤，直至温度或目标函数值达到停止条件。
影响模拟退火算法速度和质量有三个关键因素：① 状态跳跃方式；② 接受方式；③ 退火机制。
根据以上三个因素的不同可以构造不同的模拟退火算法。
① 跳跃方式
叶片位置和大小的改变服从高斯分布，高斯宽度的改变如下：
| $\sigma =1+\left( A-1 \right){{e}^{-\log \left( {{n}_{succ}}+1 \right)/T_{0}^{Step}}},$ |
式中A为初始宽度，nsucc为成功迭代次数，T0step为冷却比率[4, 6]。
② 接受方式
如果叶片位置的改变与MLC的限制条件相抵触，则重新调整叶片位置，直到可以接受为止。如果叶片位置可以接受，则计算剂量分布，确定目标函数的值，与上次迭代的目标函数值进行比价，如果减小，则接受这次迭代，如果增加，则接受概率P以标准玻尔兹曼模拟退火公式给出：
| $P=2B\frac{1}{1+{{e}^{-\log \left( {{n}_{succ}}+1 \right)/T_{0}^{prob}}}},$ |
式中B为初始概率，nsucc为成功迭代数，T0prob为降温比率。
③ 退火机制
| ${{T}_{n}}={{\alpha }^{n}}{{T}_{0}},$ |
式中T0为初始温度。
1.2 剂量计算
由于模拟退火算法需要进行至少几千次的迭代，需要快速的剂量计算算法保证整个优化的速度，因此对于迭代过程中的剂量计算一般都采用笔形束模型。
笔形束算法的基本思想是利用单一方向的笔形束核做二维卷积。基于这种思想，发展了一系列的基于笔形束核的二维卷积算法[7-9]。由于笔形束核是二维卷积，与点核叠加相比，速度更快。关于迭代过程中基于笔形束模型的剂量计算将在第二节详细叙述。
1.3 目标函数与约束
目标函数是将治疗计划输入的约束简化为某一函数数学表达式，优化过程即是寻找这一目标函数的极小值，其基本形式由公式(4) 给出：
| $\underset{m\ge 0}{\mathop{\min }}\,{{\sum\limits_{m}{{{\theta }_{m}}\left( d_{l}^{c}-d_{l}^{p} \right)}}^{2}},$ |
其中θm为感兴趣区(region of interest,ROI)的权重，dlc为ROI体元计算得到的剂量，dlp为ROI体元处方剂量。
对于MLC物理限制，如叶片最大行程、相邻叶片位置差等约束，如在新的状态下超出物理限制即不接受。同样地，对一些特殊器官，如脊髓，即使是一小部分的缺失也会造成整个器官的损伤，我们就要设置一个最大剂量确保剂量不能超出。图 1是直接孔径调强算法流程。
图1
直接孔径调强算法流程图
Figure1.
Process flowchart of DAO IMRT algorithm
2 基于点核叠加构建笔形束模型
2.1 点核叠加模型
迭加的原理可以被用于外部射束剂量分布的计算[10-11]。利用这个方法，单能光子入射到一个均匀吸收介质中得到的剂量，可以由射束的密度分布与一个空间不变核卷积得到：
| $D\left( r \right)=\int\limits_{v}{T\left( {{r}'} \right)A\left( r-{r}' \right)}{{d}^{3}}{r}'$ |
式中D(r)表示的是吸收介质中点r处的剂量，T(r′)表示Terma，是主光子在点r′处的体元d3r′中作用释放的所有能量，A(r-r′)为剂量伸展核。几何位置关系如图 2所示。
图2
点核卷积/叠加原理
Figure2.
Point kernel convolution/superposition theory
单能光子在点r′处的Terma值T(r′)可以按照如下公式解析地计算：
| $T\left( {{r}'} \right)={{\left( \frac{{{r}_{0}}}{{{r}'}} \right)}^{2}}\left( \frac{\mu }{\rho } \right)E{{\Phi }_{0}}\exp \left( -\mu {{\left| {{r}'} \right|}_{true}} \right),$ |
式中μ为介质的线性衰减因子，E为光子能量，Φ0为沿r′方向上模体表面的光子能量注量，ρ为介质的密度，｜r′｜true为在r′方向上在模体中的实际长度。
2.2 笔形束核叠加模型
在不同的算法中，笔形束核的取法也不同，可以是解析核，或者直接由蒙特卡洛程序生成的线核。下面简要介绍一下有限笔形束算法的思想[12]：
模型假设：
(1) 射束平行照射；
(2) 射束能被从几何意义上分割成有限个有限尺寸大小的射束，叫做有限大小笔形束(finite size pecil beam,FSPB)；
(3 )FSPB由蒙特卡洛方法产生；
(4 )FSPB剂量分布的迭加产生整个射束的剂量分布；
(5) 每个FSPB 都有一个与它在整个射束中的位置有关的权重注量；
(6) 每个放射束都有一个能谱，可以由有限个离散的能量段加权来模拟。
平行射束被分成多个固定大小的窄射束(笔形束)，每个窄射束可以作为一个独立的射束，根据它的能谱，就可以利用蒙特卡洛方法计算出笔形束在均匀模体中单位注量下的剂量分布，这个剂量分布就是窄射束对应的笔形束核。如图 3所示，分割后的窄射束(I,J)对应不同的能量通量强度，但具有相同的笔形束核。因此在知道各个笔形束位置的注量后，通过对各个核的迭加就可以得到整体射束下的剂量分布。
图3
笔形束核卷积/叠加原理
Figure3.
Pencil beam convolution/superposition theory
模体内每个体元的剂量Dp等于每个笔形束核对这个体元的剂量贡献的和，这个求和的过程就是剂量迭加过程。假设射束被分成n×m个笔形束，则某点P的剂量为：
| ${{D}^{p}}=\sum\limits_{IJ=0}^{nm}{D_{IJ}^{p}{{W}_{IJ}}}$ |
这个权重WIJ可以理解为笔形束(I,J)的能量注量强度。
2.3 基于点核构建笔形束模型
利用蒙特卡洛方法产生的单能笔形束核或解析核进行剂量叠加运算模型，需要重新处理多能核、核倾斜、密度修正等问题。对于那些已经建立的基于点核叠加剂量计算模型的计划系统，重新建立一套笔形束模型的研发投入是非常大的。笔形束模型假设条件与点核叠加模型有区别，尤其是对于非均匀修正区别很大。两种模型的自身差别会造成基于笔形束模型优化得到的剂量分布结果与系统原有的基于点核模型得到的结果有一定差别，这种差别会因为迭代过程的剂量计算模型简化被放大，也就造成优化结果不精确，从而无法稳定达到优化目标。
图 4(a)是射野针对整个靶区适形后得到的MLC成野，图 4(b)是在射野MLC成野范围内基于MLC叶片宽度对能量通量进行离散处理。单位能量通量下，每一个离散单元野Bij在受照射体内都对应一个剂量分布Kij。二级准直器成野到射野MLC适形野的跟随野，如图 4(a)二级准直器开野位置，MLC成封闭野状态，在单位能量通量下，在受照射体内对应剂量分布S。
图4
基于MLC适形野离散能量通量
(a)MLC适形野;(b)能量通量离散
Figure4. Discrete energy fluence based on MLC conform field(a) MLC conform field; (b) discrete energy fluence
考虑到本射野内的调强子野是本射野适形野的子集，对于不同子野下受照射体内点的剂量分布计算如下：
| ${{D}^{p}}={{S}^{p}}+\sum\limits_{ij=0}^{nm}{\left( {{K}_{ij}}-{{S}_{ij}} \right)}{{M}_{ij}},$ |
式中Dp为模体内P点处体元的剂量值，Sp为MLC封闭状态下模体内P点处体元的剂量值。Kij为离散单元野Bij对模体内P点处体元的剂量贡献。Sij为MLC封闭状态下离散单元野Bij对模体内P点处体元的剂量贡献。Mij对应离散单元野Bij状态，且
| ${{M}_{ij}}=\left\{ \begin{matrix} 0\left( {{B}_{ij}}被MLC叶片遮挡 \right) \\ 1\left( {{B}_{ij}}未被MLC叶片遮挡 \right) \\ \end{matrix} \right.$ |
考虑到调强优化过程迭代次数较多，为了减少迭代过程的剂量计算时间，一般Kij大小限制在Bij沿射束投影范围内的体元或小范围扩展。此外考虑到迭代过程基于简化的笔形束模型剂量计算是有误差的，为了减小与精确模型剂量计算的误差，在迭代过程中基于当前的子野形状会进行精确模型的剂量计算，得到在受照射体内对应的剂量分布Cij。基于当前子野下精确模型的剂量分布，下一次迭代的剂量计算为
| ${{D}^{n+1}}=C_{ij}^{n}+\sum\limits_{ij=0}^{nm}{\left( {{K}_{ij}}-{{S}_{ij}} \right)}{{L}_{ij}},$ |
式中Cijn为第n次精确剂量计算得到的模体内P点处体元的剂量值，Dn+1为第n次精确剂量计算后的下一次迭代模体内P点处体元的剂量值，Lij对应第n+1次迭代离散单元野Bij状态，且
| ${{L}_{ij}}\left\{ \begin{matrix} -1,\left( {{B}_{ij}}由打开变为遮挡 \right) \\ 0,\left( {{B}_{ij}}状态未改变 \right) \\ 1,\left( {{B}_{ij}}由遮挡变为打开 \right) \\ \end{matrix} \right.$ |
3 试验数据与结果
基于实际的患者计划数据，我们分别用点核卷积叠加和本文提出的方法做了相应的剂量计算，并对两种剂量计算得到的结果做了对比，对比结果如图 5所示。
图5
剂量比较
(a)中心轴深度剂量比较;(b)X轴离轴剂量比较;(c)Y轴离轴剂量比较
Figure5. Results of dose comparison(a) results of central axis depth dose comparison; (b) results of X off axis dose comparison; (c) results of Y off axis dose comparison
图 5(a)表示经过等中心中心轴深度剂量。图 5(b)表示经过等中心X轴离轴剂量分布。图 5(c)表示经过等中心Y轴离轴剂量分布。其中X、Y轴为与中心轴形成笛卡儿坐标系的垂直轴，纵坐标表示本文2.3节提供的方法计算得到的结果与精确点核叠加模型计算结果的比值乘以100%。绿色曲线为点核卷积叠加得到的结果，红色曲线为本文提出的方法得到的结果。
为了进一步验证本文提到的剂量计算方法在逆向优化设计中的应用，我们对模体数据进行了相应的实际优化计算。模体CT图像数据如图 6所示，其中，黄色区域为靶区计划靶体积(plan target volume,PTV)，绿色区域为危及器官Cord。
图6
模体CT图像数据
Figure6.
Phantom CT image
在优化过程中我们选择了8个野进行照射，每个野设定为6个子野。优化结果的剂量体积直方图(dose volume histogram,DVH)如图 7所示，其中红色曲线为优化前DVH曲线，蓝色曲线为优化后DVH曲线。从图中可以看出，危及器官和靶区的DVH曲线都有了相应的改进，尤其是危及器官的改进非常明显。
图7
优化前后DVH结果
Figure7.
DVH results before and after optimization
为了进行定量评估，对本文提及的剂量计算方法和精确点核叠加剂量计算的时间进行了对比试验。患者CT图像数据图像层数为76，层间距0.5 cm，剂量网格0.4 cm，计划按Gantry角度分布设计了8个射野，每个射野为10×10方野照射，分别用精确点核叠加及本文2.3方法进行剂量计算。软件运行环境为台式机Windows XP系统，硬件配置为Intel® CoreTM2 Duo CPU E8400@3.00GHz，内存2 GB，表 1是相关对比数据。
表1
剂量计算时间比较
Table1.
Comparison of dose calculation time
从表 1中数据可以看出，本文提及的方法基本上可以把迭代过程中单野的剂量计算控制在1 s以内。
4 结论
逆向调强的优化过程是基于目标函数确定子野形状及权重的最优解，在迭代过程中对于剂量计算并不要求完全精确，只要基于剂量计算结果的目标函数值与要求的收敛预期一致，就可接近或达到优化目标的设定。
考虑到上述原因及在优化结束后可基于优化结果进行精确剂量计算后评估优化结果是否满足要求，我们认为在优化迭代过程中的剂量计算更重要的是保证一定精度的前提下让剂量计算速度更快，这样才能提高逆向调强的优化速度。从试验数据及结果看，本文提出的方法可以保证优化迭代过程中的剂量计算精度满足上述要求。使得直接机器参数优化法(DMPO)在基于点核叠加剂量计算模型系统中得以实现，避免了重新研发笔形束剂量计算模型的大量工作。
引言
调强放射治疗(intensity modified radiation therapy,IMRT)[1]是通过调整多个照射野内射束的强度分布，在靶区得到高度适形剂量的同时使周围正常组织尽可能受到较少的照射，从而达到提高治疗增益比的目的。IMRT对提高肿瘤的局部控制率以及降低正常组织并发症发生率起到了非常大的作用，因此临床应用也越来越多。
调强优化的算法根据其优化求解方式，基本可以分为确定算法和随机算法两类。确定算法主要包括梯度法[2-3] 、二次规划、线性规划、迭代法等。随机算法包括模拟退火[4]和遗传算法[5]等。每种方法都有各自的优点和缺点：梯度法优化速度快，但容易陷入局部极小点；模拟退火和遗传算法是全局优化算法，可以避免陷入局部极小点，但是迭代次数多，优化速度慢。
IMRT从实现方式上可分为两大类：一类是连续或螺旋式的断层放疗，此类又可分为将射束准直为窄条形或用机架旋转进行射束调强两种，此类调强也叫扇束调强；第二大类是基于多叶准直器(multileaf collimator,MLC)的射束调强，此类调强也叫锥束调强。锥束调强可分为静态和动态两种，其区别是出束时叶片是否在伴随运动。在静态锥束调强中，调强束是用很多MLC 的均匀束分段迭加实现的，所以也叫“停-射”(Step and Shoot)模式；
传统的调强计划使用的是两步法，即首先对每个射野进行强度分布优化，得到射野强度分布图，然后将强度分布图转换为子野序列。两步法优点是可以支持动态锥束调强，缺点是靶区适形度不好。而且对于形状不规则的靶区，如鼻咽癌等，两步法优化得到的子野数及机器跳数(monitor units,MU)非常多，所需治疗时间很长。而直接机器参数优化法(direct machine parameter optimization,DMPO)通过直接计算子野的形状及权重来得到一个可实施的调强治疗计划。它的特点是物理师可以预先设置计划的子野数，从而有效地减少子野数量，减少治疗时间，提高加速器的工作效率。同时由于MLC叶片位置可直接优化得到，避免了子野分割造成的适形度误差。因此，近年来基于随机算法的直接子野优化方法越来越受到重视，其中尤以基于模拟退火优化算法的直接孔隙优化(direct aperture optimization,DAO)方法为代表。但由于随机算法迭代次数很多，而每次迭代都需要更新剂量分布，因此提高迭代过程的剂量计算速度对优化时间影响最大。在迭代过程中，为了提高剂量计算速度，一般会采用笔形束剂量计算模型，但由于基于点核叠加的剂量计算模型具有更高的精度，因此研究如何基于点核叠加构建笔形束进行优化迭代过程的剂量计算有非常大的价值，在保证系统兼容性的前提下，有望提高迭代过程中剂量计算的精度和计算速度。
本文基于已有的点核叠加模型的计划系统，在调强优化迭代过程，提出了一种由点核叠加构建的笔形束进行剂量计算的方法，在提高了剂量计算速度的同时，也保证了迭代剂量计算的算法精度及其与精确计算结果的一致性，具体见2.3节，使得DMPO在基于点核叠加剂量计算模型系统中得以实现，避免了重新研发笔形束剂量计算模型的大量工作。
1 直接孔径优化
直接孔径调强是一步调强法，利用模拟退火算法直接优化子野形状和对应照射强度，从而得到计划所需的子野信息，这种方法比传统的两步法减少了子野数和MU，提高了实施照射的效率。
直接孔径调强主要有三方面工作：① 模拟退火优化方法；② 剂量计算；③ 目标函数与约束。
1.1 模拟退火优化方法
模拟退火算法是20世纪80年代发展起来的一种用于求解大规模优化问题的随机搜索算法，此问题的自变量组合状态空间是MLC叶片的位置和此子野的强度。优化的目标函数是处方剂量与模体剂量之间的差值函数。设定初始温度、子野叶片初始状态及权重，自变量状态改变时计算模体的剂量分布，然后得到目标函数值，与上个状态的目标函数值进行比较，如果目标函数值减小，接受这个状态，否则以概率P接受。进行一定次数的状态变换后，降低温度，重复上述步骤，直至温度或目标函数值达到停止条件。
影响模拟退火算法速度和质量有三个关键因素：① 状态跳跃方式；② 接受方式；③ 退火机制。
根据以上三个因素的不同可以构造不同的模拟退火算法。
① 跳跃方式
叶片位置和大小的改变服从高斯分布，高斯宽度的改变如下：
| $\sigma =1+\left( A-1 \right){{e}^{-\log \left( {{n}_{succ}}+1 \right)/T_{0}^{Step}}},$ |
式中A为初始宽度，nsucc为成功迭代次数，T0step为冷却比率[4, 6]。
② 接受方式
如果叶片位置的改变与MLC的限制条件相抵触，则重新调整叶片位置，直到可以接受为止。如果叶片位置可以接受，则计算剂量分布，确定目标函数的值，与上次迭代的目标函数值进行比价，如果减小，则接受这次迭代，如果增加，则接受概率P以标准玻尔兹曼模拟退火公式给出：
| $P=2B\frac{1}{1+{{e}^{-\log \left( {{n}_{succ}}+1 \right)/T_{0}^{prob}}}},$ |
式中B为初始概率，nsucc为成功迭代数，T0prob为降温比率。
③ 退火机制
| ${{T}_{n}}={{\alpha }^{n}}{{T}_{0}},$ |
式中T0为初始温度。
1.2 剂量计算
由于模拟退火算法需要进行至少几千次的迭代，需要快速的剂量计算算法保证整个优化的速度，因此对于迭代过程中的剂量计算一般都采用笔形束模型。
笔形束算法的基本思想是利用单一方向的笔形束核做二维卷积。基于这种思想，发展了一系列的基于笔形束核的二维卷积算法[7-9]。由于笔形束核是二维卷积，与点核叠加相比，速度更快。关于迭代过程中基于笔形束模型的剂量计算将在第二节详细叙述。
1.3 目标函数与约束
目标函数是将治疗计划输入的约束简化为某一函数数学表达式，优化过程即是寻找这一目标函数的极小值，其基本形式由公式(4) 给出：
| $\underset{m\ge 0}{\mathop{\min }}\,{{\sum\limits_{m}{{{\theta }_{m}}\left( d_{l}^{c}-d_{l}^{p} \right)}}^{2}},$ |
其中θm为感兴趣区(region of interest,ROI)的权重，dlc为ROI体元计算得到的剂量，dlp为ROI体元处方剂量。
对于MLC物理限制，如叶片最大行程、相邻叶片位置差等约束，如在新的状态下超出物理限制即不接受。同样地，对一些特殊器官，如脊髓，即使是一小部分的缺失也会造成整个器官的损伤，我们就要设置一个最大剂量确保剂量不能超出。图 1是直接孔径调强算法流程。
图1
直接孔径调强算法流程图
Figure1.
Process flowchart of DAO IMRT algorithm
2 基于点核叠加构建笔形束模型
2.1 点核叠加模型
迭加的原理可以被用于外部射束剂量分布的计算[10-11]。利用这个方法，单能光子入射到一个均匀吸收介质中得到的剂量，可以由射束的密度分布与一个空间不变核卷积得到：
| $D\left( r \right)=\int\limits_{v}{T\left( {{r}'} \right)A\left( r-{r}' \right)}{{d}^{3}}{r}'$ |
式中D(r)表示的是吸收介质中点r处的剂量，T(r′)表示Terma，是主光子在点r′处的体元d3r′中作用释放的所有能量，A(r-r′)为剂量伸展核。几何位置关系如图 2所示。
图2
点核卷积/叠加原理
Figure2.
Point kernel convolution/superposition theory
单能光子在点r′处的Terma值T(r′)可以按照如下公式解析地计算：
| $T\left( {{r}'} \right)={{\left( \frac{{{r}_{0}}}{{{r}'}} \right)}^{2}}\left( \frac{\mu }{\rho } \right)E{{\Phi }_{0}}\exp \left( -\mu {{\left| {{r}'} \right|}_{true}} \right),$ |
式中μ为介质的线性衰减因子，E为光子能量，Φ0为沿r′方向上模体表面的光子能量注量，ρ为介质的密度，｜r′｜true为在r′方向上在模体中的实际长度。
2.2 笔形束核叠加模型
在不同的算法中，笔形束核的取法也不同，可以是解析核，或者直接由蒙特卡洛程序生成的线核。下面简要介绍一下有限笔形束算法的思想[12]：
模型假设：
(1) 射束平行照射；
(2) 射束能被从几何意义上分割成有限个有限尺寸大小的射束，叫做有限大小笔形束(finite size pecil beam,FSPB)；
(3 )FSPB由蒙特卡洛方法产生；
(4 )FSPB剂量分布的迭加产生整个射束的剂量分布；
(5) 每个FSPB 都有一个与它在整个射束中的位置有关的权重注量；
(6) 每个放射束都有一个能谱，可以由有限个离散的能量段加权来模拟。
平行射束被分成多个固定大小的窄射束(笔形束)，每个窄射束可以作为一个独立的射束，根据它的能谱，就可以利用蒙特卡洛方法计算出笔形束在均匀模体中单位注量下的剂量分布，这个剂量分布就是窄射束对应的笔形束核。如图 3所示，分割后的窄射束(I,J)对应不同的能量通量强度，但具有相同的笔形束核。因此在知道各个笔形束位置的注量后，通过对各个核的迭加就可以得到整体射束下的剂量分布。
图3
笔形束核卷积/叠加原理
Figure3.
Pencil beam convolution/superposition theory
模体内每个体元的剂量Dp等于每个笔形束核对这个体元的剂量贡献的和，这个求和的过程就是剂量迭加过程。假设射束被分成n×m个笔形束，则某点P的剂量为：
| ${{D}^{p}}=\sum\limits_{IJ=0}^{nm}{D_{IJ}^{p}{{W}_{IJ}}}$ |
这个权重WIJ可以理解为笔形束(I,J)的能量注量强度。
2.3 基于点核构建笔形束模型
利用蒙特卡洛方法产生的单能笔形束核或解析核进行剂量叠加运算模型，需要重新处理多能核、核倾斜、密度修正等问题。对于那些已经建立的基于点核叠加剂量计算模型的计划系统，重新建立一套笔形束模型的研发投入是非常大的。笔形束模型假设条件与点核叠加模型有区别，尤其是对于非均匀修正区别很大。两种模型的自身差别会造成基于笔形束模型优化得到的剂量分布结果与系统原有的基于点核模型得到的结果有一定差别，这种差别会因为迭代过程的剂量计算模型简化被放大，也就造成优化结果不精确，从而无法稳定达到优化目标。
图 4(a)是射野针对整个靶区适形后得到的MLC成野，图 4(b)是在射野MLC成野范围内基于MLC叶片宽度对能量通量进行离散处理。单位能量通量下，每一个离散单元野Bij在受照射体内都对应一个剂量分布Kij。二级准直器成野到射野MLC适形野的跟随野，如图 4(a)二级准直器开野位置，MLC成封闭野状态，在单位能量通量下，在受照射体内对应剂量分布S。
图4
基于MLC适形野离散能量通量
(a)MLC适形野;(b)能量通量离散
Figure4. Discrete energy fluence based on MLC conform field(a) MLC conform field; (b) discrete energy fluence
考虑到本射野内的调强子野是本射野适形野的子集，对于不同子野下受照射体内点的剂量分布计算如下：
| ${{D}^{p}}={{S}^{p}}+\sum\limits_{ij=0}^{nm}{\left( {{K}_{ij}}-{{S}_{ij}} \right)}{{M}_{ij}},$ |
式中Dp为模体内P点处体元的剂量值，Sp为MLC封闭状态下模体内P点处体元的剂量值。Kij为离散单元野Bij对模体内P点处体元的剂量贡献。Sij为MLC封闭状态下离散单元野Bij对模体内P点处体元的剂量贡献。Mij对应离散单元野Bij状态，且
| ${{M}_{ij}}=\left\{ \begin{matrix} 0\left( {{B}_{ij}}被MLC叶片遮挡 \right) \\ 1\left( {{B}_{ij}}未被MLC叶片遮挡 \right) \\ \end{matrix} \right.$ |
考虑到调强优化过程迭代次数较多，为了减少迭代过程的剂量计算时间，一般Kij大小限制在Bij沿射束投影范围内的体元或小范围扩展。此外考虑到迭代过程基于简化的笔形束模型剂量计算是有误差的，为了减小与精确模型剂量计算的误差，在迭代过程中基于当前的子野形状会进行精确模型的剂量计算，得到在受照射体内对应的剂量分布Cij。基于当前子野下精确模型的剂量分布，下一次迭代的剂量计算为
| ${{D}^{n+1}}=C_{ij}^{n}+\sum\limits_{ij=0}^{nm}{\left( {{K}_{ij}}-{{S}_{ij}} \right)}{{L}_{ij}},$ |
式中Cijn为第n次精确剂量计算得到的模体内P点处体元的剂量值，Dn+1为第n次精确剂量计算后的下一次迭代模体内P点处体元的剂量值，Lij对应第n+1次迭代离散单元野Bij状态，且
| ${{L}_{ij}}\left\{ \begin{matrix} -1,\left( {{B}_{ij}}由打开变为遮挡 \right) \\ 0,\left( {{B}_{ij}}状态未改变 \right) \\ 1,\left( {{B}_{ij}}由遮挡变为打开 \right) \\ \end{matrix} \right.$ |
3 试验数据与结果
基于实际的患者计划数据，我们分别用点核卷积叠加和本文提出的方法做了相应的剂量计算，并对两种剂量计算得到的结果做了对比，对比结果如图 5所示。
图5
剂量比较
(a)中心轴深度剂量比较;(b)X轴离轴剂量比较;(c)Y轴离轴剂量比较
Figure5. Results of dose comparison(a) results of central axis depth dose comparison; (b) results of X off axis dose comparison; (c) results of Y off axis dose comparison
图 5(a)表示经过等中心中心轴深度剂量。图 5(b)表示经过等中心X轴离轴剂量分布。图 5(c)表示经过等中心Y轴离轴剂量分布。其中X、Y轴为与中心轴形成笛卡儿坐标系的垂直轴，纵坐标表示本文2.3节提供的方法计算得到的结果与精确点核叠加模型计算结果的比值乘以100%。绿色曲线为点核卷积叠加得到的结果，红色曲线为本文提出的方法得到的结果。
为了进一步验证本文提到的剂量计算方法在逆向优化设计中的应用，我们对模体数据进行了相应的实际优化计算。模体CT图像数据如图 6所示，其中，黄色区域为靶区计划靶体积(plan target volume,PTV)，绿色区域为危及器官Cord。
图6
模体CT图像数据
Figure6.
Phantom CT image
在优化过程中我们选择了8个野进行照射，每个野设定为6个子野。优化结果的剂量体积直方图(dose volume histogram,DVH)如图 7所示，其中红色曲线为优化前DVH曲线，蓝色曲线为优化后DVH曲线。从图中可以看出，危及器官和靶区的DVH曲线都有了相应的改进，尤其是危及器官的改进非常明显。
图7
优化前后DVH结果
Figure7.
DVH results before and after optimization
为了进行定量评估，对本文提及的剂量计算方法和精确点核叠加剂量计算的时间进行了对比试验。患者CT图像数据图像层数为76，层间距0.5 cm，剂量网格0.4 cm，计划按Gantry角度分布设计了8个射野，每个射野为10×10方野照射，分别用精确点核叠加及本文2.3方法进行剂量计算。软件运行环境为台式机Windows XP系统，硬件配置为Intel® CoreTM2 Duo CPU E8400@3.00GHz，内存2 GB，表 1是相关对比数据。
表1
剂量计算时间比较
Table1.
Comparison of dose calculation time
从表 1中数据可以看出，本文提及的方法基本上可以把迭代过程中单野的剂量计算控制在1 s以内。
4 结论
逆向调强的优化过程是基于目标函数确定子野形状及权重的最优解，在迭代过程中对于剂量计算并不要求完全精确，只要基于剂量计算结果的目标函数值与要求的收敛预期一致，就可接近或达到优化目标的设定。
考虑到上述原因及在优化结束后可基于优化结果进行精确剂量计算后评估优化结果是否满足要求，我们认为在优化迭代过程中的剂量计算更重要的是保证一定精度的前提下让剂量计算速度更快，这样才能提高逆向调强的优化速度。从试验数据及结果看，本文提出的方法可以保证优化迭代过程中的剂量计算精度满足上述要求。使得直接机器参数优化法(DMPO)在基于点核叠加剂量计算模型系统中得以实现，避免了重新研发笔形束剂量计算模型的大量工作。
