一种治疗计划中质子重离子能量快速反演的方法

摘要

本发明公开了一种治疗计划中质子重离子能量快速反演的方法，涉及治疗计划系统领域。质子重离子能量快速反演算法基于HU值(又称CT值)把CT断层图像划分成多种介质；利用蒙特卡洛模拟建立质子重离子在这些介质中的能量衰减曲线函数；编写能量快速反演程序，提取给定束流路径上的介质，采用循环迭代的方法从粒子末端布拉格峰沉积位置处开始，基于路径上的介质种类进行能量反演积分，从而得到最终的入射能量；把反演得到的入射能量作为参数进行蒙特卡洛模拟，给出反演结果的相对误差。克服了组织使用CT值取平均后从宏观层面计算入射能量带来较大误差的影响，从CT片分辨率层面的微观尺度计算入射能量，提高了利用CT值计算质子碳离子入射能量的准确性。

说明书

技术领域

本发明涉及治疗计划系统领域，具体涉及一种治疗计划中质子重离子能量快速反演的方法。

背景技术

治疗计划系统(TPS)是一种帮助物理师制定放疗计划的应用软件，物理师首先将病人的CT断层扫描图像导入TPS中，再在三维CT视图中勾勒靶区和器官，在此基础上进行照射野的布置，并设置剂量参数信息。目前临床上使用的TPS主要由国外设计，比如RaySearch公司的RaySearch公司的RayStation TPS、Varian公司的Eclipse TPS、Elekta公司的XIO和Philips公司的Pinnacle TPS等，而我国TPS研发起步较晚，目前还没有应用于临床治疗阶段。

质子束的射程由入射能量及其所经过的人体组织相对阻止本领(relativestopping power，RSP)决定。质子成像技术可直接获得人体组织的三维相对阻止本领，但该技术仍处于实验研究阶段。因此，临床上只能借助CT图像，通过建立人体组织CT值(Hounsfield units，HU)与RSP的转换关系(HU-RSP)，间接获得相对阻止本领值。HU-RSP转换关系是导致射程精度降低的最主要因素。现在主流质子模拟定位采用水等效深度算法。水等效深度(Water-equivalent Thickness,WET)通常被用来作为质子束的射程单位。当一束相同能量的质子束穿过一定厚度的水和介质，并在最后的水箱中有相同的剩余射程，我们称这个厚度的水为该介质的水等效深度。通过建立HU-RSP的关系，再通过积分在束流路经上这些组织的水等效水深度于是可以得出入射质子能量。当前此转换关系通常使用组织替代物、动物组织和人体组织进行CT扫描和基于Bethe-Bloch的理论值建立的。在水等效深度算法中，质子模拟定位系统并不单独考虑每个像素的CT值，而是直接用一个窄束的穿透长度总和的均值来代替。这样大量减少了数据输入量与计算时间。但它会导致定位误差较大，且末端区域的剂量分布失真。

曹吉生,张伟军,彭红.治疗计划系统中利用CT值计算碳离子入射能量的可行性研究.核电子学与探测技术,2011(12):79-83。

崔相利,Jee K,王宏志,et al.质子放疗中CT值到相对阻止本领转换关系依赖因素研究及其在脂肪区的改进[J].中华放射医学与防护杂志,2019,39(10):772-777。

上述文献介绍了现有的利用HU-RSP转换关系和水等效深度进行质子定位的方法，重离子定位与质子相似。其中崔相利、王宏志等人指出了现存利用HU-RSP转换系数进行射程测量方法的不足，即同一组织中的HU值也会相差很大，所以如果只用一套使用标准膜体制作出的标准，可能会导致误差交大情况。

发明内容

本发明的目的是针对利用HU-RSP转换关系和水等效深度进行质子定位方法存在的不足，提出了一种通过蒙特卡洛模拟算法的结论，建立多个HU区间介质的层数-能量曲线，利用循环迭代，对某一固定布拉格峰沉积点的治疗计划中质子重离子入射能量快速反演的方法。

本发明具体采用如下技术方案：

一种治疗计划中质子重离子能量快速反演的方法，包括以下步骤：

a.基于HU值的有限份数的介质划分，与器官或组织无关，只要HU值属于划分的同一区间，就认为这两种器官或组织属于同一介质。

b.对划分出的介质中的能量衰减函数进行标准函数的制定。

c.能量快速反演时的积分单位为CT图像的分辨率。

d.能量的初始值为质子重离子发生布拉格峰沉积后的能量值。

e.拟合函数的能量取值为每个断层质子重离子最大通量所在能量区间的中心值。

f.采用循环迭代的方法，从质子重离子发生布拉格峰沉积时的能量反演到入射能量。

g.反演时所用的每一层介质取入射路径上经过的唯一体素的介质，(入射束流无截面积的情况下)，即当入射粒子经过一层时，认为只经过一种介质。如果束流有有限的截面积，那么取几何中心轴线路径上经过的介质序列，这是因为根据蒙特卡洛的模拟可以看出在束流路径上的粒子源分部是以轴线为中心的高斯分布，即使束流在行进的过程中发散，当质子碳离子穿过一层时绝大部分的通量依然会落在单一的体素中，附录图中可以看出，当然前提是入射坐标合理的(x0,y0)取值为离散的像素几何中心点。

上述只讨论束流垂直射入断层扫描三维几何区域的场景。

优选地，基于HU值的有限份数的介质划分时，与器官或组织无关，只要HU值属于划分的同一区间，就认为这两种器官或组织属于同一介质。

优选地，能量反演时以CT图像像素为积分单元。

优选地，基于蒙特卡洛模拟的数据拟合出的质子重离子在各种介质中的衰减函数作为反演函数。

优选地，使用粒子发生布拉格峰沉积后的能量作为初值。

优选地，先利用剖面为点状的束流进行验证，然后推广到束流剖面为有有限面积的束流。

本发明具有如下有益效果：

通过本发明公开的算法，克服了组织使用CT值取平均后从宏观层面计算入射能量带来较大误差的影响，从CT片分辨率层面的微观尺度计算入射能量，提高了利用CT值计算质子碳离子入射能量的准确性，且使用个人电脑计算一束能量的时间也在ms量级，并没有显著增加时间耗费。

本发明公开的算法能够快速确定质子重离子在确定位置处进行布拉格峰沉积时所需入射单能质子重离子束的能量。解决了现有TPS中基于人体组织进行等效水深度换算导致误差过大的劣势，利用了蒙特卡洛模拟这个治疗评估的黄金标准的结论，既拥有蒙特卡洛模拟的高精度优点，又拥有比蒙特卡洛模拟快得多的优势。

附图说明

图1为治疗计划中质子重离子能量快速反演的方法的流程示意图；

图2为某能量碳离子入射到某介质中中发生布拉格峰沉积所在断层的碳离子通量径向分布；

图3为碳离子在介质15-介质38中的能量衰减曲线和传能线密度(LET)曲线；

图4为6MeV的碳离子入射到某介质中时能谱随能量的变化；

图5为碳离子经过肺组织时的能量和LET反演曲线；

图6为利用图5反演出的参数进行正向蒙特卡罗模拟后碳离子射程末端的剂量分布。

具体实施方式

下面结合附图和具体实施例对本发明的具体实施方式做进一步说明：

结合图1，治疗计划中质子重离子能量快速反演的方法，包括以下步骤：

a.基于HU值的有限份数的介质划分，与器官或组织无关，只要HU值属于划分的同一区间，就认为这两种器官或组织属于同一介质。例如把-1000～3000的HU值区间划分成40份，每一份认为是一种介质，具有固定的元素组成，在这里使用了FLUKA的DICOM插件中自带的分配方案，如表1所示。如果想要更高的精度可以划分的更细致一些。

表1

b.对划分出的介质中的能量衰减函数进行标准函数的制定；利用蒙特卡洛软件FLUKA对250MeV的质子(用于质子治疗的能量为70MeV～250MeV)和500MeV/u(用于碳离子治疗的能量范围为85MeV/u～430MeV/u)的碳离子在所有介质中进行模拟，得到能量随穿透断层数变化的数据。需要说明的是，由于粒子在介质中的发散效应，导致能谱(能量-粒子数直方图)分布为一高斯分布，如图2所示，这时取高斯分布中心的峰值能量，即粒子数最多的能量区间的中心值为此时的能量值，可以证明此能量值的粒子处于束流轴线与当前断层相交的唯一像素点内部，如图3所示。

c.能量快速反演时的积分单位为CT图像的分辨率；把上述数据进行拟合，得到40条能量-断层序数曲线，作为反演时调用的源函数，如图4所示。

d.能量的初始值为质子重离子发生布拉格峰沉积后的能量值；选取能量反演初值，此初值的选取只与束流粒子种类有关，一般选取发生布拉格峰沉积后源粒子剩余的粒子能量作为初值，当然当发生完布拉格峰沉积后粒子的剩余能量几乎为零，也就是我们平常所说的后照射计量几乎为零，但如果选取0为初值会导致误差会变大，不推荐。

e.拟合函数的能量取值为每个断层质子重离子最大通量所在能量区间的中心值。

f.采用循环迭代的方法，从质子重离子发生布拉格峰沉积时的能量反演到入射能量。选取想要反演的位置，此参数由垂直于束流方向平面的(x0,y0)坐标和深度z0组成，即输入参数为(x0,y0,z0)，此时反演程序会自动把距离单位转换成层数，把距离转换成层数的好处是因为，当进行能量反演的时候，反演积分的最小单位是由体素大小决定的，在这里使用的CT图像序列为512*512*361，361是CT图像序列的张数，而每张CT图像的分辨率为512*512，这种情况下束流方向一定是平行于深度z轴的，故512*512中必定有一个是沿着深度z方向的层数。确定好(x0,y0,z0)后，进行束流路径上的介质的提取，最终我们得到是一个一维数组，数组中每个元素代表着介质的序号。

g.反演时所用的每一层介质取入射路径上经过的唯一体素的介质，(入射束流无截面积的情况下)，即当入射粒子经过一层时，认为只经过一种介质。如果束流有有限的截面积，那么取几何中心轴线路径上经过的介质序列，这是因为根据蒙特卡洛的模拟可以看出在束流路径上的粒子源分部是以轴线为中心的高斯分布，即使束流在行进的过程中发散，当质子碳离子穿过一层时绝大部分的通量依然会落在单一的体素中，附录图中可以看出，当然前提是入射坐标合理的(x0,y0)取值为离散的像素几何中心点。

把初始能量带入数组的第一个元素所代表介质的能量-层数方程中求解得到层数，此层数代表的含义是粒子在第一个元素所代表介质中经过这么多层后能量衰减到此数值，然后对E＝f(x)(其中x代表层数)求一阶倒数，得到dE/dx＝g(x)，把层数带入g(x)中我们得到的dE/dx就是往入射方向的第一层所需要积分的能量，也就是入射束流在经过这一断层时所损耗的能量，初始能量加上dE/dx就是下一层介质的能量-层数函数中需要输入的能量值，然后进行迭代得到入射粒子的能量，碳离子经过肺组织时的反演的曲线如图5所示。

下一步需要对此结果进行验证，把此入射能量E和位置信息作为参数带入FLUKA中进行模拟，得到此能量的粒子在体素中能量沉积的分布。由于肺与骨骼组织的密度是人体中的两个极端，在这里主要考虑了两条路径，一条经过尽可能多的肺组织，一条经过尽可能多的骨骼组织。先进行剖面为点的束流模拟，再把剖面的束流面积从无推广到有，当束流截面积不为零时每层经过的材料认为是此束流剖面几何中心点所在像素的介质。表2是在不同的粒子源形状、不同的路径上进行反演的参数、结果及误差。可以看出误差小于1％。需要说明的是，由于随着深度的增加，质子重离子的发散程度必然增加，精度也会越来越差，所以验证的深度都比较深，在20cm左右，当尺度远小于20cm时，精度也会大大增加。可以看出剖面为半径0.5cm、粒子分布为平均分布的束流，在轴线附近的模拟结果符合的也非常好，如图6所示，这是由于人体组织的在小的尺度内密度等物理量是连续的，一般不会出现突变的情况，除非束流经过器官的交界面位置处，但是这种场景是需要尽量避免在治疗计划中出现的。

表2

上述只讨论束流垂直射入断层扫描三维几何区域的场景。

上述反演算法过程中，基于蒙特卡洛模拟的数据拟合出的质子重离子在各种介质中的衰减函数作为反演函数。使用粒子发生布拉格峰沉积后的能量作为初值。先利用剖面为点状的束流进行验证，然后推广到束流剖面为有有限面积的束流。

当然，上述说明并非是对本发明的限制，本发明也并不仅限于上述举例，本技术领域的技术人员在本发明的实质范围内所做出的变化、改型、添加或替换，也应属于本发明的保护范围。