基于临床质子加速器的DeepPlan笔形束模型构建

张增鹏，陶寅，刘红冬，陈志，徐榭，4，裴曦，4

1.中国科学技术大学核科学技术学院，安徽合肥，230025；2.中国科学技术大学核医学物理研究所，安徽合肥230025；3.中山大学附属肿瘤医院放疗科，广东广州510060；4.安徽慧软科技有限公司，安徽合肥230088

前言

相比光子，质子放射治疗具有独特的Bragg 峰，在靶区给予更高剂量的同时，降低对肿瘤后端器官的辐射损伤［1-2］，因此在全球范围高速发展。根据国际粒子治疗协作组（Particle Therapy Co-Operative Group,PTCOG）的统计，截止2020年3月份，全球正在运营的质子治疗中心已达91 个，我国正在运营的质子放疗中心仅有2家，仍有7家正在建设中［3-4］。可以看到，随着越来越多质子中心的投入使用，质子治疗在放射治疗中的比重会逐渐增大［5］。与此同时，与质子治疗相配套的治疗计划系统也会越来越受到重视。目前，质子治疗计划系统几乎被国外各大厂商垄断，如瑞典RaySearch 的RayStation，美国Varian 的Eclipse，瑞典Elekta 的Xio 和荷兰Philips 的Pinnacle等［6-7］。我国质子放疗计划研究起步较晚，国内拥有自主知识产权的放疗系统主要有中国科学技术大学放射医学物理中心研发的DeepPlan 和中国科学院核能安全技术研究所研发的KylinRay［8-9］。两者质子模块均采用笔形束算法，通过质子束在体模中的纵向积分深度剂量和横线束斑展宽计算体模中各体素的剂量沉积。目前由于缺乏质子放疗临床数据，国产质子放疗计划系统开展的临床数据对比实验较少。

本研究旨在研究DeepPlan 质子模块中建立佛罗里达大学质子放疗中心（University of Florida Health Proton Therapy Institute，UFHPTI）质子加速器笔形束模型，通过对比UFHPTI质子加速器在笔形束扫描模式下的临床实验数据，验证模型构建的准确性并与UFHPTI配套的商用放疗计划软件RayStation就前列腺癌临床数据进行剂量对比。

1 材料与方法

1.1 DeepPlan质子模块

1.1.1 模块简介DeepPlan 是由中国科学技术大学放射医学物理中心自主研发的放疗计划系统［9-12］。其质子模块采用笔形束算法，公式如下：

其中，D（x,y,z,x0,y0,E）表示能量为E，中心坐标为（x0,y0）的质子束在模体位置（x,y,z）的剂量沉积，IDD是质子束在纵向的积分深度剂量，利用蒙卡软件TOPAS模拟得到，w表示核散射权重，Gmcs表示多库伦散射的高斯分布，Gnuc表示核散射的高斯分布，σmcs是能量为E的质子束由于多库伦散射在x轴或y轴方向上高斯分布的标准差。σnuc是能量为E的质子束由于核散射在x轴或y轴方向上高斯分布的标准差，这些参数根据Soukup等［13］、Gottschalk等［14-15］的工作计算得到。

DeepPlan 质子模块的工作流程为：首先通过数据读取得到模体信息和质子束信息，包括模体位置、大小、体素材料、等中心点位置、质子能量、入射方向、初始束斑大小和质子束权重等，并导入计算所需要的数据库，如上述计算所需要的IDD 数据库、多库伦散射数据库、核散射数据库和核散射权重数据库等。接着光线追踪通过设定两个位置点P1、P2计算质子束穿过体素的几何距离和等效水深度。最后剂量计算根据等效水深度得到计算所需的数据参数，对每个质子束利用公式（1）～（3）得到体模的总剂量。这3个部分实现笔形束扫描下质子剂量的快速计算。

1.1.2 模型构建

1.1.2.1 数据读取数据库构建的准确性是剂量计算的前提条件，比如IDD 数据库构建时不同初始能散和角度发散质子束的IDD 是不一样的。图1是用蒙卡软件TOPAS模拟的200 MeV质子在不同能散下的IDD数据曲线。由图可知，相同能量不同能散质子束的IDD 曲线在Bragg 峰前端有一定差距，能散越小，IDD 曲线越低，这样根据式（1）计算的剂量会比实际情况偏低。本研究结合UFHPTI 质子加速器信息并结合刘红冬等［16］的研究，插值计算得到加速器输出各种能量的参数信息，构建包括IDD在内的数据库。

图1 200 MeV质子束在不同能散的积分深度剂量Fig.1 Integrated depth dose of 200 MeV proton beam in different energy spreads

1.1.2.2 光线追踪光线追踪是剂量计算的重要一环，算法程序根据Siddon［17］、Han 等［18］的工作实现。该算法认为构成体模的每一个体素不是独立的，而是一系列等距平行平面的正交集焦点体积，加上质子束穿过体模时的两个位置点信息，可实现位置点连线上各体素几何距离和等效水深度的快速计算。需要注意的是在实际笔形束扫描模型中，质子加速器通过X 和Y 方向扫描磁铁实现对质子束位置的精确控制，因而有X 和Y 两个焦点，如图2所示。以UFHPTI质子加速器为例，其X、Y 两个焦点与等中心平面的距离分别为2 329 和1 934 mm。因此，在设置两个位置点进行光线追踪时必须考虑X、Y 方向之间的差异。

1.1.2.3 剂量计算剂量计算过程需要考虑的因素较多，包括空气处理和散射修正等，公式如下：

图2 笔形束扫描模型虚拟源示意图Fig.2 Diagram of virtual source in pencil beam scanning model

其中，式（4）用于在空气中进行剂量计算，主要是因为在空气中等效水深度为0，且在空气中只有多库伦散射没有核散射。式（5）用于在非空气体模中进行剂量计算，与式（1）相比，该公式增加了多库伦散射修正项δmcs和核散射修正项δnuc，表达式如上所示。在（8）～（9）式中多库伦散射修正因子cmcs和核散射修正因子cnuc分别是与多库伦散射标准差σmcs和核散射标准差σnuc有关的常数项。

图3是剂量计算示意图，在参数导入之后首先应该判断各体素是否为空气，判别方法为计算各体素的质量组织本领之和T，若T＜10-6，则认为在空气中，采用式（4）进行剂量计算；否则，采用式（5）进行剂量计算。

图3 剂量计算示意图Fig.3 Diagram of dose calculation

1.1.3 测试环境本文完成DeepPlan质子模块临床测试的计算机配置如下：操作系统为Windows10 专业版（64 位），CPU 为Intel(R)CoreTMi7-87100K、六核十二线程、16 G内存，GPU 为NVIDIA GeForce GTX10606GB。DeepPlan 质子模块有c++和c#两个版本，本研究基于c#版本在软件Visual Studio 2017上运行。

1.2 UFHPTI质子加速器临床实验数据测量

1.2.1 IDD 数据测量UFHPTI 使用IBA 公司的BPC探测器（Bragg Peak Chamber,model 34070,PTW,Freiburg,Germany）在IBA 水箱中进行积分深度的测量，共有30 组质子能量，分别为：226.35、223.91、219.06、214.12、209.16、204.54、199.63、195.1、190.18、185.74、180.9、176.21、171.89、167.79、163.02、158.09、153.04、149.07、144.81、140.63、136.07、131.51、127.38、122.46、118.69、115.09、110.04、106.57、102.36、99.59 MeV 等。针对测试数据，在DeepPlan质子模块中计算上述30组能量的质子束在水模中的剂量沉积，体模大小为100 mm×100 mm×400 mm，每个体素的大小为1 mm×1 mm×1 mm，对垂直于质子束方向的每个平面进行剂量求和，利用最大值进行归一即可得到在纵向的积分深度剂量曲线。

1.2.2 空气中束斑尺寸测量UFHPTI 使用IBA 公司的Lynx PT 闪烁探测器进行束斑尺寸的测量。30 组能量在上面已经提过，5 个测量位置分别为距离等中心平面20、10、0、-10、-20 cm。为了将二维数据转成一维，需要将测试数据导入IBA的MyQA软件提取X或Y方向上的横向剂量分布，并做归一处理。

针对测试数据，在DeepPlan 质子模块中计算30组能量的质子束在体模为空气中的剂量沉积，平面大小为40 mm×40 mm，平面精度和Lynx PT 的精度保持一致，为0.25 mm×0.25 mm，将计算结果在X 或Y方向求和得到一维数据，并做归一处理。

1.2.3 SOBP 和横向剂量的测量为保证剂量分布在纵向和横向上与靶区保持适形，需要使用不同能量、不同权重、不同位置的质子束。在SOBP 测量时UFHPTI 使用IBA 公司的多层电离室Zebra 测量中心轴上的绝对剂量；在横向剂量测量时UFHPTI 使用Lynx PT 闪烁探测器进行测量。针对测量数据，在DeepPlan 质子剂量模块根据提供的质子束能量信息、权重信息和位置信息等，在水模中进行剂量计算，水模大小为200 mm×200 mm×400 mm，体素大小为2 mm×2 mm×2 mm。

1.2.4 前列腺癌的剂量计算RayStation 是由RaySearch 公司于2008年开始研发，2011年市场发布的治疗计划系统［19］。RayStation 的技术主要有轮廓的勾画、形变配准及定量评估、自适应治疗、剂量追踪、多标准优化和后备计划等。作为与UFHPTI质子加速器相配套的治疗计划系统，RayStation 在笔形束模型下的质子剂量计算有两种算法，一种是笔形束解析算法，另一种是蒙卡算法。笔形束算法比蒙卡算法计算速度快，但准确性低于蒙卡算法：对于大量非均匀介质体模的计算和所有质子束使用射程移位器的情况等，主要以蒙卡算法为主。本文以RayStation 蒙卡算法的计算结果进行对比。在UFHPTI 的两个前列腺癌临床放疗计划中，相关信息如下：病例1的体模大小为512 mm×512 mm×361 mm，体素大小为1.171 9 mm×1.171 9 mm×1.000 0 mm，270°和90°方向的质子束总数分别为1 039 和977；病例2的体模大小为512 mm×512 mm×293 mm，体素大小为1.171 9 mm×1.171 9 mm×1.000 0 mm，270°和90°方向的质子束总数分别为1 307 和1 303。由于RayStation 中给出的三维剂量信息体素大小为2 mm×2 mm×2 mm，为了进行gamma 通过率对比，需要对病例1 和2 的体模按相同体素大小进行重采样。最后使用PTW 公司的VeriSoft 软件对二者的计算结果进行gamma分析。

2 结果与讨论

2.1 IDD结果比较

为展示方便，图4所示将30 组IDD 数据从高到低分为6 组（a～f 组），每组5 个能量进行比较，每组能量的IDD 曲线都利用该能量下最大值进行归一，深度方向进行适当截取，DeepPlan 质子模块计算结果用DeepPlan 表示，UFHPTI质子加速器临床实验数据用Measure 表示。从图中可以看到，DeepPlan 质子模块计算结果与UFHPTI 质子加速器临床实验数据在Bragg 峰前后曲线相吻合。两者平均相对误差为0.01%，少数能量如图4f 组中的99.59、102.36 MeV等，DeepPlan 计算结果与实验数据在最大值位置附近出现一个0.23%的相对误差。出现这个情况的原因可能是测试数据最小深度间隔为3 mm，精度较低，导致在测量峰值最大位置时出现一定误差。

图4 30组不同能量下IDD结果对比Fig.4 Comparison of integrated depth dose under different energies

2.2 空气中束斑尺寸结果对比

如图5所示，以能量为199.63 MeV 的质子束为例，蓝色曲线表示UFHPTI质子加速器的临床测试结果，红色曲线表示DeepPlan 质子模块的计算结果，图中包含5个不同位置的对比结果，束斑宽度由小到大分别为距离等中心平面20、10、0、-10、-20 cm 处的结果。对于某一位置的计算结果和测量数据，导入MATLAB R2018b 版本中，使用软件自带的高斯拟合函数求出束斑大小。图6是不同能量质子束在等中心平面束斑大小的对比结果，大部分相对误差在0.5%以内，只有163.02 MeV 比实际测试略低一些，相对误差为1.13%。除此之外，表1展示距离等中心平面20、10、-10 和-20 cm 处的计算结果。从表中可知，DeepPlan 质子模块的计算结果与临床实验数据大部分相对误差在0.3%以内，平均相对误差为0.15%。可以看到，由于质子束在空气和体模中的输运不同，在剂量计算时将二者分开可使结果更加准确。

图5 199.63 MeV质子束在不同位置结果对比Fig.5 Comparison of 199.63 MeV proton beams at different positions

图6 不同能量质子束在等中心平面的束斑尺寸Fig.6 Spot size of proton beams with different energies in the isocenter plane

表1 30组不同能量质子束在空气不同位置的束斑比较Tab.1 Comparison of the spot size of 30 groups of proton beams with different energies at different positions in the air

2.3 SOBP和横向剂量结果对比

图7是纵向和横向的结果对比，黑色曲线表示UFHPTI 质子加速器的临床测量结果，红色曲线表示DeepPlan 质子模块的计算结果。由图7a 在纵向扩展布拉格峰的对比中，可以看到实现射程为200 mm，展宽100 mm 目的。DeepPlan 剂量计算结果与临床实验数据平均相对误差1.07%，最大相对误差3.91%。图7b 展示的是等中心平面横向剂量分布的对比，等中心平面在距离水模表面15 cm 的位置，可以看到平坦区域达到100 mm 目的。DeepPlan 剂量计算结果与临床实验数据平均相对误差1.92%，最大相对误差4.09%。

图7 纵向和横向结果对比Fig.7 Comparison of longitudinal and transverse results

2.4 前列腺癌结果对比

DeepPlan与RayStation的结果在该平面二维gamma通过率达到了97.1%，但在末端右下角未通过，出现一些冷点和热点，见图8。可能原因是相比DeepPlan解析算法，RayStation采用蒙卡算法，对非均匀介质部分计算结果更好，这是后续工作需要改进的地方。表2具体展示了两个前列腺癌病例在两个射野方向的三维gamma通过率，测试标准为3 mm/3%。由表2可知每个射野方向在大于10%处方剂量区域里的gamma通过率都达到95%以上，其中病例1在两个方向的gamma通过率分别为96.4%和97.5%，病例2在两个方向的gamma通过率分别为99.3%和98.9%。

图8 DeepPlan与RayStation关于前列腺癌的二维剂量结果比较Fig.8 Comparison of two-dimensional dose calculation results of DeepPlan and RayStation on prostate cancer

表2 DeepPlan与RayStation关于前列腺癌的gamma分析结果Tab.2 Gamma analysis results of DeepPlan and RayStation on prostate cancer

3 结语

本文对DeepPlan 质子模块在数据读取、光线追踪和剂量计算等方面建模，通过与UFHPTI质子加速器在30组IDD对比、30组空气中质子束斑发散对比、SOBP 和横向剂量分布对比，验证DeepPlan 质子模块模型构建的准确性。同时以两个前列腺癌的放疗计划为实例，与商用放疗计划系统RayStation 进行gamma 分析，表明DeepPlan 质子模块在笔形束扫描模式下可以初步应用于临床前列腺癌的剂量计算。后续还将继续对治疗计划系统DeepPlan 质子模块进行改进，包括对非均匀介质的修正和加入射程移位器的剂量计算等，对更多类型的病例进行测试分析，进一步增强DeepPlan 质子模块在临床应用的可行性。

感谢UFHPTI提供的临床实验数据。