Collapsed Cone光子束剂量计算方法研究

周正东 宋 威

(南京航空航天大学核科学与工程系，南京 210016)

引言

放射治疗是目前治疗恶性肿瘤的主要手段之一，随着计算机、医学影像、放射物理以及放射生物等技术的发展，以X射线为放射源的适形和调强放射治疗技术迅速发展。放射治疗计划系统作为实现放射治疗精确定位、精确计划、精确治疗目标的必要手段具有重要的意义。放射治疗计划的目的是尽可能选择最好的照射方案，以达到对靶区施行高剂量照射，在产生不可恢复性摧毁的同时，使周围正常组织及关键或敏感组织所受影响最小。剂量计算是放射治疗计划系统的核心，collapsed cone卷积/叠加剂量计算方法[1]是一种基于点核模型的剂量计算方法，该方法较常用的笔形束剂量计算方法计算精度高，但计算速度相对较慢，自从该方法提出后，有关此算法的研究一直受到众多学者的关注[2-8]。我们对collapsed cone卷积/叠加剂量计算方法进行了深入研究，在此基础上开发了一套光子束放射治疗计划系统，并运用蒙特卡罗方法对剂量计算的准确性进行了验证。

1 Collapsed Cone剂量计算方法

卷积/叠加算法是目前放射治疗计划系统常用的光子束剂量计算方法，它是一种基于核模型的算法。常用的卷积核包括点核(point kernel)，笔形束核(pencil beam kernel)、面核(planar kernel)等。与笔形束核模型相比较，点核模型可以处理电子失衡区域的剂量计算问题[9]，计算精度更高，不足之处是计算速度较慢，但是通过算法的改进，在保证一定计算精度的前提下，能够满足临床对计算速度的要求。以点核为模型的卷积算法分为两步[10]，首先由质量衰减系数计算原射线穿过人体组织时单位质量释放的总能量 (total energy released per unit mass，TERMA)，然后根据事先计算得到的点核，经过组织密度不均匀修正，与TERMA卷积，最终得到人体组织中的剂量分布。

对于具有n3个体素的三维数据，以点核为模型的标准卷积剂量计算算法的时间复杂度为O(n6)，如果再考虑组织密度的修正，计算复杂度将达到O(n7)，难以满足快速计算的要求。Ahnesjö提出了collapsed cone卷积算法[1]，此算法的时间复杂度为O(mn3)，其中 m为能量沉积核的沉积方向数。collapsed cone卷积算法将散射剂量的积分划分为具有一定立体角的同轴锥形筒串，位于筒串轴的散射体积单元释放的能量沿筒串轴线性传递、衰减、沉积。

式(1)为点核 hm(r)的表达式，其中 am、bm、Am，Bm为某一能量沉积方向的沉积系数，可以通过对蒙特卡罗程序计算得到的能量沉积核进行拟合得到。

对各个能量沉积方向，按照式(2)分别计算原射线剂量和散射线剂量并求和

式中，m、n分别为相对于射束入射方向的天顶角和方位角的离散方向索引，天顶角和方位角的角增量分别记为Δσθ和Δφ，r为剂量计算区域中的任意一点的位置向量，分别表示原射线和散射线剂量。

对于每条能量沉积射线上体素的原射线剂量和散射线剂量，可分别通过式(3)和式(4)迭代求出[1]，此两式等号右侧第一项为某一方向能量沉积线上前方体素ri-1对当前计算体素ri的剂量贡献，第二项为当前体素对自身的剂量贡献。

式中，Ti为当前体素单位质量释放的总能量，Ωmn为某一能量沉积方向筒串的立体角，ηi为当前体素相对水的电子密度，用于对非均匀介质的密度修正。

在collapsed cone算法中用平行于剂量沉积方向的一系列射线，来替代每一个体素的三维卷积方向进行加速计算。为了快速计算平行的能量沉积射束与病人体数据网格的交点以及相交线段的长度，Ahnesjö提出了种子平面的概念[3]，穿过种子平面的平行射束与病人网格的交点，可以简化为求其中一条射线与网格的交点坐标，然后与射线在种子平面上的入射点坐标相加获得。整个剂量计算有以下步骤。

步骤1：计算单位质量释放的总能量(TERMA)

步骤1-1：根据机床角、机架角、机头角、射野的形状计算待计算区域的矩形包围盒;

步骤1-2：计算包围盒内各体素的 CT值、质量衰减系数、质量密度、电子密度;

步骤1-3：根据射野的形状和计算分辨率，计算射束内射线的离散数目n;

步骤1-4：计算每条射线在包围盒内沉积的TERMA。

步骤2：计算剂量

步骤2-1：根据设定的天顶角增量Δθ和方位角增量 Δφ，将方位角和天顶角离散为(θi，φj);

步骤2-2：对于每一个能量沉积方向(θi，φj)计算种子平面的大小、能量沉积射线的数量、过原点的能量沉积射线与包围盒内体数据网格的交点和相交线段长度，然后根据坐标平移关系，计算每一条能量沉积射线与体数据网格的交点坐标，进而计算每一条能量沉积射线在包围盒内沉积的剂量，并累加。

2 软件系统结构及功能

为了验证上述剂量计算方法的准确性，我们在VS 2005平台上开发了一套光子束放射治疗计划系统，其软件结构如图1所示。主要包括DICOM数据的读取、图像分割、三维组织重建、计划设计、剂量计算、剂量评估、数据库管理七个模块，每个模块又可以细分为若干子模块。系统首先读入DICOM序列，以位图格式显示在软件界面上;通过手动分割或者智能剪刀(live wire)分割方法[11-13]提取靶区和关键组织的轮廓线后，使用轮廓线重建算法、marching cubes表面重建算法或体绘制算法进行三维重建;根据勾画出的靶区，设定加速器及射野参数等，进行剂量计算;将计算得到剂量数据以等剂量面、等剂量线或者以伪彩色的方式显示在病人断层数据上，并计算靶区及敏感组织的剂量体积直方图，以便进行剂量评估。最后把轮廓线数据、计划设计的参数、剂量计算的参数及结果等存放到指定的路径下，并把文件路径保存到数据库中，软件系统界面示例图如图2所示。

图1 放射治疗计划系统软件结构Fig.1 Software structure of radiation treatment planning system

3 实验结果与分析

为了验证上述剂量计算方法的准确性，采用蒙特卡罗程序DOSXYZnrc进行比较，在相同能谱、射野大小、水体模等条件下进行了模拟计算，比较射线能量分别为6 MV和15 MV条件下的百分深度剂量(PDD)，及20 cm×20 cm方野6 MV射束条件下不同深度的离轴比。源皮距SSD=100 cm，水体模体积为30 cm×30 cm×30 cm，射野大小为3 cm×3 cm、5 cm ×5 cm、10 cm ×10 cm、20 cm ×20 cm，DOSXYZnrc程序的历史次数为108，collapsed cone算法中天顶角增量和方位角增量(Δθ、Δφ)取(7.5°、22.5°)、(3.75°、22.5°)、(3.75°、15°)等 3 种情况，用来分析离散角对计算精度的影响。剂量计算平均误差见表1，可以看出：射束能量增加误差减少;射野增大误差增加;计算分辨率提高误差降低，当 Δθ、Δφ 取 3.75°、15°时剂量计算误差最小，图 3给出了该情况下，本算法与蒙特卡罗程序DOSXYZnrc计算结果比较的示例图，图中“DOSXYZnrc”表示蒙特卡罗程序 DOSXYZnrc计算的结果，“CC”表示本 collapsed cone算法计算得到的结果。图3(a)是射线能量为6 MV，照射野分别为 3 cm ×3 cm、5 cm ×5 cm、10 cm ×10 cm、20 cm ×20 cm条件下的百分深度剂量。图3(b)是射线能量为15 MV，照射野分别为 3 cm×3 cm、5 cm×5 cm、10 cm×10 cm、20 cm×20 cm条件下的百分深度剂量。图3(c)是射线能量为6 MV，照射野大小为 20 cm ×20 cm，深度为 1.5、10、30 cm 条件下的离轴比。从图3中的(a)和(b)可以看出：1)随着射野或能量增大，百分深度剂量增加;2)当深度较小时，collapsed cone算法的计算结果与蒙特卡罗模拟的结果一致，从图3(c)亦可看出并可得出此结论，但随着深度的增加，由于射束随深度增加而硬化，导致能量沉积核发生变化，因而误差略有增加;3)相比于高能射线，低能射线硬化效应相对明显，因而深度较大时的误差相对较大;4)在射野较大的情况下，由于每一能量沉积方向需要计算的距离增大，因而相对于小射野，射野大时误差略大。

表1 深度剂量计算平均误差(%)表Tab.1 Mean calculated dose error

图3 与蒙特卡罗程序DOSXYZnrc的剂量比较示例图。(a)6 MV射线下百分深度剂量;(b)15 MV射线下百分深度剂量;(c)6 MV射线20 cm×20 cm射野下离轴比Fig.3 Comparision of dose with DOSXYZnrc.(a)PDDs for 6 MV beam;(b)PDDs for 15 MV beam;(c)OARs for 6 MV beam with field size of 20 cm×20 cm

4 结论与讨论

对collapsed cone剂量计算的步骤以及加速算法进行了深入研究，在VS 2005平台上开发了一套功能相对完善的光子束放射治疗计划系统，运用蒙特卡罗方法和水体模对该剂量计算方法的准确性进行验证。与蒙特卡罗方法比较结果表明，本文研究开发的collapsed cone剂量计算方法准确可靠，具有重要的临床应用价值。百分深度剂量与射线的能量有关，临床应根据病灶的具体情况选择合适能量的射线进行治疗。与常用的笔形束剂量计算方法相比，collapsed cone剂量计算方法具有更高的计算精度，其剂量计算精度与剂量计算分辨率(即collapsed cone算法的离散参数)及射束硬化等因素有关。随着计算分辨率的提高，计算时间相应增加，计算复杂度增加、速度降低;射束硬化会导致剂量计算误差增加。对于射束硬化引起的计算误差，通过对能谱进行硬化校正，可以减小计算误差，提高剂量计算的精度;由于剂量需要在三维空间进行计算，且在放射治疗计划优化的过程中需要反复进行剂量计算，因此在确保剂量计算精度的前提下，剂量计算的速度也至关重要，GPU具有高效的并行计算能力且成本低，近几年来在诸多领域得到众多研究人员的关注，对于剂量计算同样具有较高的实用价值。今后我们将研究基于GPU的剂量计算方法，以提高剂量计算的速度，同时研究剂量计算中的射束硬化校正方法，进一步提高剂量计算的精度，并完善软件系统功能。

[1]Ahnesjö A.Collapsed cone convolution of radiant energy for photon dose calculation in heterogeneous media[J].Med Phys，1989，16(4)：577-592.

[2]Hobant PW，Murray DC，Round WH.Photon beam convolutionusing polyenergetic energy deposition kernels[J].Phys Med Biol，1994，39(4)：669-685.

[3]Ahnesjö A.Cone discretization for the collapsed cone algorithm[C]//Proceedings of International Conference on the Use of Computers in Radiation Therapy.XII ICCR.Madison：Leavitt DD and Starkschall G Medical Physics Publishing，1997：114-116.

[4]Maspradakis M， Morrison RH， Richmond N， et al.Experimental verification of convolution/superposition photon dose calculations for radio therapy treatment planning[J].Phys Med Biol，2003，48：2873-2893.

[5]Lu Weiguo， Olivera1 GH， Chen Mingli， et al. Accurate convolution/superposition for multi-resolution dose calculation using cumulative tabulated kernels[J].Phys Med Biol，2005，50(4)：655-680.

[6]Zhou B，Hu XS，Chen DZ，et al.Hardware acceleration for 3-D radiation dose calculation[C]//Proceedings of 18 th IEEE International Conference on Application-specific Systems，Architectures and Processors.Montreal：IEEE，2007：290-296.

[7]Carlsson TA，Ahnesjö A.Optimization of the computational efficiency of a 3 D，collapsed cone dose calculation algorithm for brachytherapy[J].Med Phys，2008，35(4)：1611-1618.

[8]Hissoiny S， Ozell B. Fast convolution-superposition dose calculation on graphics hardware[J].Med Phys，2009，36(6)：1998-2005.

[9]陈朝斌，黄群英，吴宜灿，等.蒙特卡罗方法在放疗计划中的应用[J].核技术，2006，29(1)：22-28.

[10]Ahnesjö A，Aspradakis MM.Dose calculations for external photon beams in radiotherapy[J].Phys Med Biol，1999，44(11)：99-155.

[11]Mortensen EN， Barrett WA. Interactive segmentation with intelligent scissors[J].Graphical Models and Image Processing，1998，60(5)：349-384.

[12]高新波，雷云，姬红兵.一种改进的 Live-Wire交互式图像分割算法[J].系统工程与电子技术，2003，25(8)：915-918.

[13]杨金柱，赵姝颖，胡英，等.序列医学图像三维分割的一种方法[J].系统仿真学报，2005，17(12)：2896-2900.