基于EGSnrc程序的锥形束CT模拟仿真研究

付 娟，孟慧鹏，*，艾克文，顾 欣*，王克强，张艳龙，丁红军，段敬豪，孙 倩

（1.武警特色医学中心，天津 300162；2.天津大学精密仪器与光电子工程学院，天津 300072）

0 引言

机载式锥形束CT已经广泛用于调强放射治疗（intensity modulated radiotherapy，IMRT）的图像引导，成为位置验证及跟踪肿瘤变化等方面的有力工具[1-3]。随着放疗技术向自适应放射治疗（adaptive radiotherapy，ART）发展，越来越多的学者开始关注基于锥形束CT进行形变配准、剂量计算等新的需求[4-6]。然而为方便此类研究的顺利开展，首先需要解决的问题是如何对现有锥形束CT进行有效的模拟仿真。因此本研究拟通过Monte Carlo模拟程序EGSnrc对锥形束CT进行模块化建模，并对所建立的系统进行验证，进而探讨所建立的仿真系统是否能够准确模拟锥形束CT的主要技术指标，及其是否具有进一步研究其他锥形束CT相关问题的潜能。

1 材料和方法

1.1 仪器设备

实验设备为Varian IX直线加速器（美国，瓦里安公司）；三方能谱文件由开源程序Spektr（Version 3.0）生成；测量时采用型号为MT-VW的红固体水，密度为1.03 g/cm3；百分深度剂量测量使用德国PTW公司的MP3-M大水箱，离轴剂量测量使用该公司的OCTAVIUS Detector 729二维电离室矩阵；模拟仿真程序的版本为EGSnrc v2016（2016年1月发布），在Ubuntu 16.04的64位系统下运行[7]。

1.2 方法

1.2.1 主要部件参数

拟模拟仿真的锥形束CT系统集成于Varian IX直线加速器上，主要组成构件为X射线源和非晶硅平板探测器。X射线源球管为Varian公司生产，型号为G242，峰值电压和阳极最大热容量分别为125 kVp和445 kJ，阳极材质采用钨合金，靶光束角为14°，初级滤波为2 mm厚的铝合金；平板探测器同样为Varian公司生产，型号为4030CB，有效探测面积为39.7 cm×29.8 cm，矩阵尺寸为 1 024×768，像素尺寸为0.388 mm×0.388 mm，探测器可沿探测平面横向向外偏移 16 cm，最大扫描视野（field of view，FOV）可达到46.6 cm[8]。

1.2.2 锥形束CT的模拟仿真

本研究通过EGSnrc的BEAMnrc子程序将锥形束CT系统分为X射线源、射线源组件、探测器、Monte Carlo模拟用模体共计4个模块建模，所建立的系统示意图如图1所示。根据扫描时采用的领结滤波不同，分为全扇束模式和半扇束模式。实验使用的锥形束CT系统扫描协议共有3种扫描电压，而125 kVp对应2种领结滤波模式，因此共建立4种源模型：100 kVp全领结滤波、110 kVp半领结滤波、125 kVp全领结滤波和半领结滤波。

图1 EGSnrc模拟的CBCT系统示意图

EGSnrc的BEAMnrc子程序可提供多种建模组件模块，通过不同组件的组合可实现锥形束CT的模拟仿真。可用于模拟锥形束CT系统的组件有XTUBE、CONSTAK、SLABS、BLOCK、JAWS和PYRAMIDS等，模拟仿真时，由SLABS组件构建初级滤波、射束硬化滤波、机载影像系统窗和探测器4种结构；领结滤波既可以用PYRAMIDS组件构建，也可以用JAWS组件构建，本研究中采用PYRAMIDS组件。每个组件的最终输入参数如下：

（1）X射线源。

①XTUBE组件构建靶：由Varian公司提供的球管参数文档可知，阳极靶物质为钨合金，其中主要成分为金属钨（占比95%）、金属铑（占比5%），等效物理密度为18.68 g/cm3。此外，靶光束角为14°，入射电子的焦距为0.4~0.8 mm，本实验模拟定义的焦斑直径为0.5mm。②CONSTAK组件构建球管出窗：玻璃窗的密度设置为2.23 g/cm3，厚度设置为0.9 mm；窗口与靶之间的距离设置为2.8cm，出窗与靶之间设置为真空。

（2）射线源组件。

①SLABS组件构建初级滤波：初级滤波为2 mm厚的铝合金，与靶的距离为5.72 cm。其中铝合金的等效物理密度为3.23 g/cm3，包含96.5%的铝、1.3%的硅、1.2%的镁和1%的锰。②BLOCK组件构建初级准直器：初级准直器的材质为2 cm厚的铅，其上端开口（靠近源方向）和下端开口（远离源方向）分别为3.2cm×3.2cm和4.2cm×4.2cm，与靶的距离为6.52cm。③JAWS组件构建挡片：挡片的材质为0.2 mm厚的铅，同准直器的设置类似，上下端开口分别设置为4.72 cm×4.72 cm和4.86 cm×4.86 cm，与靶的距离为9.59 cm。④SLABS组件构建射束硬化滤波：射束硬化滤波为0.5 mm厚的钛，物理密度为4.54 g/cm3，与靶的距离设置为11.19 cm。⑤SLABS组件构建机载影像系统窗：机载影像系统窗为1 mm厚的透光有机物聚碳酸酯，与靶的距离设置为14.9 cm。⑥PYRAMIDS组件构建领结滤波：构建全扇束领结滤波和半扇束领结滤波，领结滤波的材质为铝合金，与初级滤波的铝合金成分相同。建模时由6层铝层组成（每层中心开孔不同），从上到下的铝层厚度分别为 1.5、1.5、7、13、3 和 2 mm，与靶的距离为 15 cm。

（3）探测器。

实验采用的是Varian 4030CB平板探测器，BEAMnrc子程序建模时主要考虑两个部分：一个是装配于探测器阵列上方的10∶1的抗散射栅，另一个是平板探测器铜层和由无定形硅构成的平板薄膜晶体管阵列（闪烁体层）。摆位采集到的锥形束CT图像原始投影信号实际上就是探测器阵列对到达探测器粒子的响应，建模时X射线的响应可通过提取该层的能量沉积实现。对于千伏级别的X射线，抗散射栅自身和平板探测器铜层与粒子的相互作用不可忽略，因此将探测器建模成1 mm厚的铜层，采用SLABS组件构建，物理密度设置为8.9 g/cm3，与靶的距离为150cm，设置平板探测器的大小为39.7cm×29.8cm，矩阵尺寸为1024×768，像素尺寸为0.388mm×0.388mm。

（4）Monte Carlo模拟用模体。

为了验证BEAMnrc子程序建立系统的准确性及后续实验，建立的固体水模尺寸为30 cm×30 cm×30 cm，同时建立相同尺寸的标准水模体。

1.2.3 仿真系统的验证

为了验证建立系统的几何和光束特性是否正确，通过X射线能谱、百分深度剂量、离轴剂量等进行分析比较。Monte Carlo模拟上述参数的结果分两步生成：第一，在BEAMnrc子程序中相应平面生成相空间文件；第二，由相应子程序分析相空间文件得到对应结果，其中BEAMDP子程序可获得X射线能谱，DOSXYZnrc子程序可获得百分深度剂量和离轴剂量。比对数据的生成如下：由开源程序Spektr生成X射线能谱，由实测数据获取百分深度剂量和离轴剂量。

（1）X射线能谱。

用于比对的能谱文件由开源程序Spektr（Version 3.0）生成，分别生成100、110和125 kVp的能谱文件。具体设置时，在Spektr程序的UI界面中选择Tube_09，并在tubeSettings函数中将Tube_09的值自定义为[13 2；74 0]，在“ADDED FILTERATION”模块中，增加Ti（钛，Z=22），“Thickness”设置为0.5 mm，以确保生成的能谱与本实验设备的参数尽可能一致。

Monte Carlo模拟X射线能谱时，对应锥形束CT系统扫描协议的3种X射线管电压，实验验证进行了3次主要的Monte Carlo模拟。每次运行均得到对应的相空间文件，该文件生成于模拟仿真中的一个虚拟平面，可存储所有通过该平面的粒子的能量、位置、方向等相关数据，由BEAMDP子程序分析后得到X射线能谱。由于Spektr程序只能生成纯钨靶X射线能谱，且只能生成距离源100 cm处的能谱文件，为了便于比较，在生成分析X射线能谱的相空间文件模拟中，本实验将所建立系统中的靶材料替换为纯钨（其余模拟均为钨合金），相空间文件记录平面设置为距离靶100 cm处。X射线能谱的Monte Carlo模拟采用的粒子数为1×109，在距离靶100 cm处模拟的射野大小为50 cm×50 cm，模拟时未加载领结滤波组件。

（2）百分深度剂量和离轴剂量。

实测数据均需锥形束CT射线源在0°处，因此选择设备的服务模式，该模式可自定义相关扫描参数，如管电压、管电流、扫描时间、有无领结滤波、挡片开口大小等，所有数据均在透视模式下源皮距（source skin distance，SSD）100 cm处测量得到。为验证所建立的仿真系统的准确性，分别实测了100kVp全领结滤波、110 kVp半领结滤波、125 kVp全领结滤波和半领结滤波的百分深度剂量和离轴剂量。百分深度剂量测量时设置挡片野大小为26 cm×26 cm，在射束稳定状态下，采集射束中心轴0~20 cm范围内的剂量；离轴剂量测量时将0.5 cm厚的固体水放置于电离室矩阵之上，采集固体水下0.75 cm平面内等中心位置横向和纵向的剂量，射野大小为20 cm×20 cm。

Monte Carlo模拟百分深度剂量和离轴剂量时，对应锥形束CT系统实测数据的4种电压-滤波模式进行4次主要的Monte Carlo模拟，共记录4组相空间文件，记录位置为源到等中心距离（source isocentre distance，SID）100 cm 处，Monte Carlo 模拟使用的粒子数为1×109。记录的相空间文件使用DOSXYZnrc子程序分析计算，使用的模体大小为30 cm×30 cm×30 cm，计算体素大小为 0.5 cm×0.5 cm×0.5 cm，粒子数为5×108，光子分裂数设置为100，对所有感兴趣区域的体素计算中统计不确定性控制在2%以内。DOSXYZnrc子程序的输出文件格式为*.3ddose，该文件可由STATDOSE子程序分析。提取出射束中心轴0~20 cm范围内的剂量和模体深度为1.25 cm平面内等中心位置横向和纵向的剂量，绘制生成百分深度剂量和横向、纵向离轴剂量曲线与实测结果进行比较。

2 结果

在SID 100 cm平面处，3种管电压100、110和125 kVp的Monte Carlo模拟结果与Spektr程序生成的对应比较结果如图2所示，所有能谱信息均与峰值信息做了归一化处理。结果表明，2种能谱具有良好的一致性，差异为0.6%~3.0%，平均差异小于2.0%。由图2还可知，3种管电压的最大差异区域均出现在30~55 keV范围内，其中125 kVp管电压的36 keV能量处差异达到了3.0%。

图2 仿真系统生成的能谱信息与Spektr程序生成的能谱信息的比较

Monte Carlo模拟的百分深度剂量和离轴剂量结果与实测的结果一致性良好，如图3~5所示，所有的数据均对最大值做了归一化处理。图3为百分深度剂量分布的比较，4种电压-滤波模式的射束中心轴剂量差异均较小，剂量差异平均值在±2.0%以内。其中100 kVp全领结滤波差异最大，且表现为全领结滤波模式[如 3 图（a）和（c）所示）]的差异大于半领结滤波模式[如图 3（b）和（d）所示]，且水下 5 cm厚度以内的Monte Carlo模拟值与实测值的差异大于5 cm以外的。差异最大处在100 kVp全领结滤波的水下5 cm范围内，达到了2.7%；4种模式水下5 cm以外差异值均逐渐减小，均小于1.0%。图4为横向离轴剂量分布，用于验证系统加载全领结滤波和半领结滤波过滤器的几何设计，4种电压-滤波模式横向离轴剂量的差异平均值在±1.1%以内。由图4可知，110 kVp半领结滤波模式Monte Carlo模拟值与实测值在4种模式中差异最大，最大处达到了4.0%，该模式中的差异平均值为1.3%。图5为纵向离轴剂量分布比较，4种电压-滤波模式纵向离轴剂量的差异平均值在±0.8%以内，最大差异为2.7%。由图5可知，4种模式的曲线有类似的形状，其中125 kVp全领结滤波和半领结滤波模式的曲线肉眼无法区分其差异。

图3 仿真系统中心轴百分深度剂量曲线的Monte Carlo模拟数据与实测数据比较

图4 仿真系统等中心位置横向离轴剂量曲线的Monte Carlo模拟数据与二维电离室矩阵实测数据比较

3 讨论

ART技术最终向剂量引导放射治疗（dose-guided radiotherapy，DGRT）[9]发展，实现实际患者体内剂量沉积的实时追踪[4]。通过治疗时采集的患者的摆位锥形束CT图像可间接实现患者体内累积剂量的追踪，是ART技术当前较易实现的一种方法。然而，基于锥形束CT进行形变配准、剂量计算的研究[5-6]必须基于一套准确高效的模拟仿真系统。本研究在充分了解Varian OBI系统（锥形束CT系统）的构造及特点的基础上，基于Monte Carlo模拟软件的EGSnrc程序实现了锥形束CT系统的模拟仿真。所建立系统计算得到的百分深度剂量、离轴剂量等物理参数与实测值误差平均值均在±2.0%以内，与文献[8，10-11]中其他学者的建模方法相比，各种几何结构及物质组成更接近真实情况，如采用钨合金和铝合金代替其他学者采用纯钨和纯铝的模拟、考虑平板探测器的1 mm铜层的模拟等。此外，通过改变各组件有限的参数，所建立的系统可对其他主流锥形束CT系统进行模拟，具有良好的通用性，可用于模拟研究其他锥形束CT系统及锥形束CT系统的其他问题。

值的一提的是，实验中模拟生成的能谱未能与实际测量的能谱数据进行比对，主要原因是能谱的测量复杂，实际测量较困难，且测量精度受多种因素制约，因此本实验采用很多研究使用的成熟的第三方能谱生成软件来生成比对能谱信息。能谱生成程序采用的是Spektr最新版本，其基于三次样条插值算法的钨阳极光谱模型（tungsten anode spectral model using interpolating cubic splines，TASMICS）的算法，相对于之前采用多项式插值的钨阳极光谱模型（tungsten anode spectral model using interpolating polynomials，TASMIP）算法的老版本进行了改进，3.0版本允许用户自定义球管参数、不同厚度材料的过滤器及对生成的能谱文件优化，可生成20~150 keV能量范围内的X射线能谱，其模拟出的能谱信息与实际数据的差异在4.0%以内，平均差异在2.7%以内[10]。本中心Ming等[11-12]的前期工作表明，10%的能谱差异可导致约4%的射束中心轴百分深度剂量差异和小于1%的器官剂量差异，因此仿真系统生成的能谱文件中射束中心轴百分深度剂量差异均较小，可满足锥形束CT系统的相关研究。对于百分深度剂量曲线，在5 cm以内深度处差异较大的原因与模拟时未考虑实际照射的多向散射及电离室的精度有关，其均会带来浅表区域的差异增加；对于能谱信息30~55 keV能量范围内的差异大于其他区域的可能原因是初级滤波和射束硬化滤波对低能射线的过滤放大了该区域的差异，此外本实验建模的1 mm铜层探测器也对低能部分影响较大。

图5 仿真系统等中心位置纵向离轴剂量曲线的Monte Carlo模拟数据与二维电离室矩阵实测数据比较

横向离轴剂量曲线全领结滤波模式在离轴距离[-10 cm，10 cm]区间内差异略大于[-13 cm，-10 cm）和（10 cm，13 cm]区间，可能的原因是模块化建模时忽略支撑领结滤波的钢支架（厚度小于1 mm）导致的。全领结滤波结构为中间开孔，两边厚，钢支架对中心位置的剂量影响被放大，从而导致领结滤波开孔区域内的剂量差异变大。半领结滤波模式在离轴距离[-13 cm，6 cm]区间内的差异略大于（6 cm，13 cm]区间，原因与全领结滤波类似。尽管存在差异，但是横向离轴剂量的差异整体在较小的范围内，对使用该仿真系统进行锥形束CT的其他研究影响甚微。

纵向离轴剂量曲线中，125 kVp全领结滤波和125 kVp半领结滤波模式的曲线肉眼无法区分其差异，究其原因与领结滤波的结构有关。领结滤波的中心开孔方向与纵向平行，导致125 kVp全领结滤波采集数据时领结滤波对其造成的影响与半领结滤波相同（仅与钢支架发生作用）。领结滤波在[-10 cm，10 cm]区间内纵向离轴剂量模拟值与实测值吻合性更好，其余区域差异略大，究其原因是实际测量使用的OCTAVIUS Detector 729二维电离室矩阵相邻2个电离室的中心间距为1 cm，而模拟时的体素计算间距为0.5 cm，分辨力的差异造成剂量变化较快的区域出现更大误差，此外差异可能还与二维电离室矩阵的测量精度有关。分析该组数据还发现，Monte Carlo模拟值与实测值差异最大的模式为125 kVp全领结滤波，出现在离轴距离约-12 cm处，差异值为2.7%。同样，纵向离轴剂量的差异整体上也在较小的范围内，4种电压-滤波模式纵向离轴剂量的差异平均值在±0.8%以内，对使用该仿真系统进行锥形束CT的其他研究时的影响比横向的影响更小。

综上所述，本研究建立的模拟仿真系统可实现锥形束CT的准确模拟，可为研究锥形束CT的相关问题奠定基础。下一步将基于该系统对如何提高锥形束CT图像质量及基于该系统的形变配准算法等展开研究。