放射治疗作为癌症治疗中的主要手段，其目标是在保证杀死癌细胞的同时，最大限度保护周围正常组织。相对光子放疗，带电粒子具有更加明显的优势，在穿过介质过程中，带电粒子最大能量损失于有限射程末端，即布拉格峰处，峰位深度后端几乎没有剂量。1955年美国加州大学Lawrence Berkeley实验室开展了世界上第一例质子治疗乳腺癌的病例^[1]。随后的几十年，德国、俄罗斯、瑞士、瑞典、日本等国家相继开展粒子治疗的研究。最新的研究表明，粒子治疗具有替代外科手术的潜力^[2]，只是现阶段还缺少大量的随机临床数据对比。我国的山东万杰医院于2004年治疗首例病人，上海的质子重离子医院、台湾长庚医院也分别于2015年5月和11月正式收治患者；上海交通大学附属瑞金医院的国产首台质子治疗装置正于中国科学院上海应用物理研究所研发，预计2018年开展临床试验。据PTCOG数据，截止2016年3月，全世界正在运行的质子中心有53家，正在建设31家，计划建设16家^[3]，已有超过十万人接受质子治疗。近期国内涿州、天津、广州、武汉、合肥、长沙、郑州等地区内数家单位均有建设质子中心意向。

为充分利用质子治疗的物理特性优势和实现精准放疗，研究人员开发了三维的铅笔束(Pencil Beam)扫描照射方法，这一技术已经成功运用在瑞士PSI (Paul Scherrer Institute) ^[4]、德国GSI (Gesellschaft für Schwerionenforschung)^[5]、日本NIRS (National Institute of Radiological Sciences)^[6]、美国MD安德森(MD Anderson Cancer Center)^[7]。相比被动散射体照射方法^[8]，扫描照射在不需要准直器和补偿器情况下可以实现更加适形的靶区剂量分布。

本文采用自主开发的仿真程序对点扫描照射进行仿真模拟，分别考虑靶区静止和受呼吸运动影响两种情况，为正在研制的装置和将来的临床应用提供必要的参考。

质子治疗装置主要包括：加速器、束流输运系统、治疗头、治疗控制系统、摆位系统、治疗计划系统、剂量验证系统、质子放疗信息系统等组成。用于质子治疗的加速器目前主要有同步加速器和回旋加速器两种。基于同步加速器的质子治疗系统的基本构成如图 1所示。质子在离子源中引出后，经过射频加速器RFQ (Radio Frequency Quadrupole)，DTL (Drift Tube Linac)加速至一定能量，进入主加速器，然后经输运线被引出到不同治疗室，通过照射系统对患者实施治疗。

图 1

图 1 质子治疗装置系统示意图 Figure 1 Schematic of proton treatment system.

通常情况下，从加速器直接引出的质子束不能直接用于患者治疗，因为其能量单一，束斑尺寸小。此时需要利用治疗头把引出的束流分别实施横向扩展和纵向扩展，形成所需的照射野。在使用扫描照射方法的质子治疗系统中，纵向扩展通过加速器主动改变能量使束流停在不同深度，最终使靶区体积内形成均匀的剂量分布。

治疗控制系统是基于计算机网络的分布式实时控制系统，采用大量的软硬件实施治疗控制。控制系统的主要目标^[9]是：治疗过程中，确保患者和工作人员的安全；将满足治疗性能的质子束按治疗计划精确传输到肿瘤区。

摆位系统是用来保证和验证照射野位置。摆位的准确性依赖多模式图像采集，患者摆位、固定，摆位验证等步骤。利用摆位系统将照射系统的等中心点调整到治疗计划指出的等中心点位置，并且保证病人的姿态与治疗计划中一致。

治疗计划系统(Treatment Planning System, TPS)是协助医生设计和制定患者治疗计划与方案的应用软件，包括完成特定功能的应用程序。目前的TPS主要包括图像导入、图像处理（多种影像融合、靶区勾画等）、摆位用的数字重建影像(Digitally Reconstructed Radiography, DRR)的生成，计划设计、剂量计算和计划优化、评估以及精度验证用计划的生成QA-plan (Quality Assurance)等功能。

剂量验证系统是照射质量控制的重要工具，利用其实施计划和照射系统的QA。质子放疗信息系统(Radiation Information System, RIS)主要包括病人排程、诊疗和日常管理信息系统，放疗记录验证系统，患者影像存档与通讯系统，还能够传输治疗计划和其他治疗信息。此系统对实现安全高效的质子治疗具有重要作用。

从同步加速器引出的质子束的尺寸一般只有几个毫米，为能够满足不同大小肿瘤的治疗，必须使束流横向扩展到整个照射野，并且使能量覆盖肿瘤的纵向深度范围。束流横向扩展有多种不同方法，主要分为两类：被动散射体法、动态扫描法^[10]。被动散射体法是束流经散射体和射程调制器的作用扩展为均匀覆盖整个治疗体积的扩展布拉格峰(Spread-out Bragg Peak, SOBP)的方法；动态扫描中束流每次只照射靶体积的一小部分，然后改变束流的位置直至照射所有的靶体积。相比动态扫描，被动散射技术简单易行，仍在多数质子中心得到应用。但是动态扫描无需病人专用的射野挡块和射程调制器，粒子利用率高，同时产生次级碎片和中子也少，被新的质子治疗中心广泛采用。动态扫描可以形成三维适形的剂量分布，还能够实现调强质子束治疗(Intensity Modulated Proton Therapy, IMPT)，对复杂形状的肿瘤更加有效。现在及未来的趋势是更多的利用动态扫描束流输运法。

动态扫描法分为两类：点扫描(Spot Scanning)^[11]和光栅扫描(Raster Scanning)^[12]。在使用这两种技术时，沿深度方向把肿瘤分成不同的等能量层，每一层用相同能量的质子照射，每层有一定数量的照射点，一般从最深的那一层开始照射，待该层所有照射点达到预设剂量，再改变质子束流能量，照射下一层，依次逐层扫描照射，直至全部的等能量层及照射点得到照射。

点扫描时如图 2所示，当一个照射点照完之后，关断束流，改变扫描铁电流，移动束流至下一个照射点的位置，继续照射直至照完计划给出的所有照射点；在光栅扫描中，在点与点间移动时，不关断束流，而是快速移动束流位置，以缩短扫描照射所需的时间。

图 2

图 2 点扫描照射示意图 Figure 2 Schematic of spot scanning irradiation.

为实现目标的剂量分布，点扫描中每个照射点的粒子数必须预先计算，束流能量和质子数决定了每个照射点的三维剂量分布，束流能量确定后，需要计算每个照射点的粒子数，即权重w。靶体积被分为众多的计算网格(voxels)，每个voxel处剂量D_i由所有照射点对它剂量贡献之和得到，如式(1) ：

式中：d_j、w_j代表第j个照射点的剂量分布、粒子数。

考虑所有计算网格后，通过逆向优化目标函数式(2) ，计算每个照射点的粒子数。

式中：D_{p, i}表示靶区体积V中第i个网格点处的计划剂量（考虑物理剂量优化或生物剂量优化时对应的剂量），计算时一般使用拟牛顿算法迭代，当目标函数迭代至预设条件时，得到了各个照射点的权重因子w_j。

本文中计算权重因子时采取如下简化处理：1) 不考虑目标函数式(2) 中靶区周围危机器官(Organs at risk, OAR)的贡献；2) 靶内各等能量层中的照射点均匀排布，且同一层中所有照射点采用相同的权重，因模拟所用的靶为长方体水靶，密度处处相等，虽然靶区边缘的半影非最优，但不影响本文结论；3) 利用特定能量的质子在水中的深度剂量分布，通过移动此曲线在靶内深度方向上的位置来替代加速器主动换能，这样三维的权重优化问题就简化为只在深度方向上优化该曲线的权重的一维优化问题。

扫描束缺点是对有位置变动的器官很敏感^[13]，因此照射运动靶区时需要谨慎应对。器官运动的种类有多种^[14]，包括位置依赖的运动、压力导致的运动（如膀胱充盈度）、血液搏动引起位置变化（主要是心脏周围）、呼吸运动（肺癌、肝癌等）。呼吸运动被认为是胸腔和腹腔内器官运动的重要因素，常运用三角函数来描述器官受呼吸影响的运动状态：

式中：A、T、f₀、r₀、r(t)分别表示呼吸运动的幅度、周期、初始相位、初始位置、t时刻位置。

器官运动时，照射点位置误差会带来较大的剂量误差，本文采用门控技术^[15]解决呼吸运动对剂量分布的影响。门控技术基本思想如图 3所示：1) 为呼吸运动设置一个幅度阈值，小于此值所对应时间段才有可能照射，找到满足要求的呼吸门时间；2) 加速器束流的周期内可供束流的时间和呼吸门时间的交集作为真正照射时间。

图 3

图 3 门控技术示意图 (a) 患者呼吸信号，(b) 呼吸阈值，(c) 满足阈值的呼吸信号，(d) 加速器束流时间结构，(e) 实际照射时间 Figure 3 Schematic of gated technique. (a) Patients’ respiration waveform, (b) Threshold level for gated irradiation, (c) Gated respiration signal, (d) Periodic spill structure extracted from synchrotron, (e) Real irradiation time

点扫描模拟计算采用自主开发的仿真程序，设计思想如下：首先设置靶区（确定计算网格尺寸，分配照射点位置）；设定照射点的权重，权重对应各点将要照射的绝对粒子数；绝对粒子数和束流信息确定照射点的照射时间；当在对应的照射时间内完成对束流流强的积分，移动至下一个点照射；每个计算网格内的剂量是通过所有照射点的粒子数和其剂量分布在此处的贡献的卷积运算得到，最终按要求输出剂量分布和相关结果。

1) 靶区和计算范围：目标靶区选择长方体，材料为水靶，尺寸为100 mm×100 mm×60 mm，靶中心距水模表面70 mm，相对位置以及照射方向如图 4所示。

图 4

图 4 靶区示意图 Figure 4 Location and size of target.

2) 扫描参数：X、Y、Z方向上照射点扫描间隔均为3 mm，照射点共有33×33×20=21780个。束斑形状为X、Y方向标准偏差都是为5 mm的高斯分布。X、Y方向的扫描速度分别为20 mm·ms^-1和5mm·ms^-1。在靶区内总计给予2 Gy的处方剂量。

3) 质子束深度剂量分布：选用动能为190 MeV的质子在水中的深度剂量分布的数据，将对应权重的该剂量分布移动至所有等能量层，如图 5所示。以此方式近似替代加速器主动改变能量，形成宽度为6 cm的SOBP，各等能量层的厚度为3 mm，每层对应的该深度剂量分布曲线的相对权重见表 1。

表 1

表 1 剂量深度分布在20个等能量层中的相对权重 Table 1 Relative weights of the depth dose curve in 20 iso-energy layers. i ii iii iv v vi vii viii ix x xi xii xiii xiv xv xvi xvii xviii xix xx 120 50 35 28 25 20 19 17 16 15 13 12 11 11 10 9 10 9 8 9

表 1 剂量深度分布在20个等能量层中的相对权重 Table 1 Relative weights of the depth dose curve in 20 iso-energy layers.

图 5

图 5 不同权重的深度剂量分布平移后形成的6 cm的SOBP Figure 5 The SOBP of 6-cm overlapped by a series of depth dose curves with the optimized weight.

4) 加速器束流参数：参照国产首台质子治疗装

置的主要设计指标，单次引出质子数目(4-8) ×10¹⁰，重复频率0.5-0.1 Hz；模拟时选取加速器运行周期为3.5 s，可供引出的flattop时间为1.5 s，假设1.5 s内引出4×10¹⁰个质子，理想情况下，引出束流没有波动，流强为(4×10¹⁰/1.5) /6.241×10⁹≈4.3 nA（1 nA流强时，每秒6.241×10⁹个质子通过）。做保守的处理，模拟时采用3.0 nA的流强。

5) 加速器关断响应延迟：关断延迟为200 μs。

靶区静止时，不考虑§3.1中加速器关断响应延迟的剂量分布，见图 6。本文所显示剂量均以相对靶区中心点(X=75 mm、Y=75 mm、Z=70 mm)的剂量百分比给出，图 6(a)显示深度Z方向上靶区前端路径上剂量逐渐变大，在靶区内形成60 mm的SOBP，靶区尾部剂量迅速跌落；图 6(b)为垂直于靶中心Z=70 mm的横向XY平面内剂量分布，可以看出剂量分布平坦；图 6(c, d)分别为上述平面内X、Y方向的剂量分布，图 6(c)中W为照射野的宽度（以射束中心轴剂量为100%，相对两边50%等剂量线间距离），0.8 W内剂量均匀度均小于2%，与照射野内均匀的剂量相比，边缘剂量迅速跌落，形成约6 mm半影区（80%-20%等剂量线之间的宽度）。

图 6

图 6 靶区剂量分布 (a) 靶中心X=75 mm、Y=75 mm处纵向深度剂量分布，(b) Z=70 mm处横截面内剂量分布，(c) Z=70 mm、 Y=75 mm处X方向剂量分布，(d) Z=70 mm、X=75 mm处Y方向剂量分布 Figure 6 Dose distribution in target. (a) Depth dose at the target center of X=75 mm, Y=75 mm, (b) Dose distribution at the cross section of Z=70 mm, (c) Dose distribution in the X direction at Z=70 mm, Y=75 mm, (d) Dose distribution in the Y direction at Z=70 mm, X=75 mm

照射所用质子总数约为1.5×10¹¹个，束流周期数21个，总照射时间约74 s。最深的等能量层用时两个束流周期，每个点照射时间1970 μs；最浅的等能量层用时一个束流周期，每个点照射时间148μs；剂量率为1.6 Gy∙min^-1（靶0.6 L时）。

考虑200 μs加速器关断延迟后，剂量分布明显变差。靶区中心沿深度Z（X=75 mm、Y=75 mm时）方向的剂量分布如图 7所示，在深度40-100 mm范围内，越靠近水模表面，模拟所得到的剂量与计划剂量相差越大。造成剂量分布变差原因是：照射点达到计划剂量后，关断束流移至下一个照射点，在加速器关断响应延迟时间内，仍有额外的质子照射在该点上。而权重随靶等能量层深度Z变浅呈下降趋势，所以额外剂量与计划剂量的比值随深度变浅逐渐变大。当流强降低后，此比值变小。

图 7

图 7 靶中心纵向深度方向剂量比较 Figure 7 Comparison of depth doses at target center.

实际照射过程中，束流流强是有波动的，且在关束延迟这段时间内，流强呈逐渐下降趋势，额外剂量可能比上述模拟值小。当流强为1 nA时，以更短的100 μs的关断延迟来近似等效关断束流时的流强下降，可以看到中心深度剂量与计划值的差均小于2%。消除额外剂量影响的办法是：选择合适的引出流强大小，或者把额外剂量的影响反馈到照射点的权重计算中。前者可能会加长照射时间，降低剂量率，后者实施过程会比较复杂，是值得进一步研究的课题。

对于平躺的患者，放松安静的呼吸时，肺等器官受呼吸运动影响最大的是从头到脚S-I (superior-inferior)方向，参照图 5中X方向，该方向运动幅度最大可达30 mm (peak-to-peak)，而左右 (left-right, L-R)方向和平行于束流深度 (anterior- posterior, A-P)方向峰值幅度最大仅为8-10 mm^[15]，分别对应图 5中Y和Z方向。X方向对应扫描速度较快的方向，Y方向对应扫描速度较慢的方向。为验证门控技术的效果，本文模拟X方向运动对剂量分布的影响，以式(2) 随时间变化的周期函数近似替代，其中运动幅度、周期、初始相位、初始位置分别为20 mm、4 s、0 rad、0 mm，呼吸门控阈值为5mm。小于该阈值的适合照射时间为1s。呼吸运动时，是否采用呼吸门控的剂量分布结果对比如图 8所示，在靶区中心(Z=70 mm)横截面内，剂量分布在X方向随呼吸运动平移，靶区内靠近边界一侧的照射点剂量低于计划值，差值最大可达近80%，而另一侧额外部分体积却受到照射。在所选的运动参数下，照射该层时靶区正向图 8(a)中剂量不足一侧运动。模拟结果显示，当呼吸运动初始相位、运动幅度、周期变化时，剂量分布结果随之变化。直接照射运动靶时，靶区内没有得到足量照射，而附近正常组织和危机器官却得到超出计划的照射剂量。当采取呼吸门控后，靶内X方向剂量分布明显改善，边缘很快下降，满足计划要求。

图 8

图 8 呼吸和采取门控时剂量分布对比 (a) 有呼吸运动靶中心Z=70 mm处横截面XY内剂量分布，(b) 采取呼吸门控后Z=70 mm处横截面XY内剂量分布， (c) Z=70 mm、Y=75 mm处X方向剂量分布，(d) X=75 mm、Y=75 mm处沿深度剂量分布 Figure 8 Results of dose distribution on the moving target. (a) Dose distribution of moving target at the cross section of Z=70 mm, (b) Dose distribution utilizing gated irradiation at the cross section of Z=70 mm, (c) Dose distribution in the X direction at Z=70 mm, Y=75 mm, (d) Depth dose at the target center of X=75 mm, Y=75 mm

模拟结果显示总照射时间为90 s，剂量率为1.3Gy∙min^-1（靶0.6 L时）。采用呼吸门控，靶区只有满足呼吸阈值时才可能得到照射，因此相比照射静止靶区，耗时更长，剂量率变低。

本文利用仿真程序模拟质子治疗中点扫描照射技术。不考虑加速器关断响应延迟情况下，照射静止靶得到均匀的剂量分布；照射受呼吸影响的运动靶区时，采用门控技术显著改善剂量分布，达到计划要求。考虑关断束流响应延迟条件下，纵向剂量分布随深度变浅逐渐变差。照射剂量率信息对正在研制的装置提供重要参考。目前仿真程序网格剂量叠加耗时较长，下一步考虑使用GPU并行加速这一过程。