李 贵 郑华庆 吴宜灿,3 FDS团队

1（中国科学院等离子体物理研究所 合肥 230031）

2（安徽省精确放疗工程技术研究中心 合肥 230031）

3（中国科学技术大学核科学技术学院 合肥 230027）

基于多算法的电子剂量场反演方法

李 贵1,2郑华庆1,2吴宜灿1,2,3FDS团队

1（中国科学院等离子体物理研究所 合肥 230031）

2（安徽省精确放疗工程技术研究中心 合肥 230031）

3（中国科学技术大学核科学技术学院 合肥 230027）

为了研究精确放射治疗中多算法的电子剂量场反演方法及反演效果，建立非线性规划电子剂量场反演模型，利用体外透射剂量重建全空间电子剂量场，并通过多种算法实现。采用 7个例题测试结果表明：除NMinimize算法不能实现该模型外，其他算法都能实现，其中Levenberg-Marquardt算法最精确快速。进一步测试表明，利用接近电子射程的透射剂量会带来较大的反演误差（尾差），原因可能是能谱平均带来的，因此改进方法需考虑改进对平均能量的近似。但由于反演时间较短，精度高，可实现实时剂量场反演；并且每一种算法都有特点，为实现提供多方法的选择性，满足临床上对实时剂量验证的需要。

电子剂量场，剂量验证，能谱，精确放射治疗

精确放射治疗中，蒙特卡罗模拟等精确剂量计算可达到很高的正向计算精度[1−3]，但实际治疗操作中的不可控误差[4,5]，会使照射区域的剂量分布偏移计划剂量分布，因此需实施实时的剂量场的重建即剂量场反演，以确保患者体内肿瘤区域的剂量符合计划要求。目前，精确放疗剂量验证的方法基本考虑光子照射，但也须验证电子束照射的剂量，方能确保“精确放疗”的实现。此外，电子剂量场反演还可应用于工业探测，如用高能电子束照射工件，测量透射剂量来反推工件内部剂量场分布，从而获得工件的损伤与材料分布等情况。

电子束照射剂量验证的难点，一是电子射程短，射程之外当然无电子剂量分布；二是光子污染问题，光子剂量场可以很好地控制电子污染，而电子剂量场中的光子不能很好地消除。第一个问题的解决方案，只有提高电子束能量，如20 MeV电子束能穿透10 cm水模，可穿透人体一些部位（如乳腺癌的治疗），可实施剂量验证，因此本文把电子剂量场限定在高能区。对于第二个问题，本文将把光子场也考虑在内，建立一个混合场模型。

本文发展一种新型的电子剂量场的反演方法，在混合笔形束模型基础上，建立非线性规划反演模型并利用多种非线性反演方法实现快速精确的电子束剂量场实时反演，最后与蒙特卡罗模拟对比[6−10]进行验证。

1 基于混合笔形束的非线性规划反演模型

1.1 剂量场误差分析与反演模型的建立

Luo等[11−13]根据 Fermi-Eyges多级散射理论，发展了混合笔形束模型，用于精确放射治疗的精确解析电子剂量计算，其模型决定电子光子混合剂量场的参数有：电子能谱、矩形野尺寸、散射参数、光子污染函数与“距离平方反比修正”函数中的若干系数。本文将在此模型基础上建立电子剂量场反演模型。

电子剂量场的反演如图1的 C→B过程，根据体外C处的透射剂量反推B的剂量分布。其中要求电子能穿透B穿过空气到达C处，这样才能根据C处测得到的透射剂量分布反推出全空间的剂量分布。此法把图中 ABC三部分剂量场看作一体，把问题归结为求解非均匀模体的反向剂量计算问题。

基于以上分析，对如下几种情况的变化造成空间的剂量场重新分布进行讨论。

(1) 矩形野变化。设变形后的矩形野为(x∈[D2,D1],y∈[C2,C1])，未形变的矩形野为(x∈[d2,d1],y∈[c2,c1])，由此引起的矩形变化(负号表示变短)为：

野长变化：

野宽变化：

图1 电子剂量场反演Fig.1 Inversion of electron dose map.

(2) 射线平移。只考虑水平面的平移，射线上下平移的实质是源皮距变化，这里用临床常用的平行源，因此源皮距的变化对剂量场的影响很不敏感，可忽略。由于射线位置与射野是相对的，射线平移导致的后果就等同于射野偏离了原射野中心：[(d1+d2)/2, (c1+c2)/2]。新的射野中心可通过变形后的矩形野求得。新的射野中心为：[(D1+D2)/2, (C1+C2)/2]，则新射野中心对原射野中心的偏移量为：

(3) 射线转动。射线发生转动等同于坐标系转动。根据坐标转动公式，设θ和φ角分别是法线方向n的极偏移角和方位偏移角，则旧坐标(x,y,z)和新坐标(X,Y,Z) (实验坐标系)的变换关系为

(4) 光子污染变化。我们曾在先前的电子能谱反演研究中提出[14,15]：光子污染由两个指数函数描述。假设高能散射光平均衰减系数µp、轫致辐射光子和次级光子衰减系数µe、散射光子与轫致辐射光子的比例系数ν等决定光子物理性质的参数不变，并设电子剂量场的重新分布使光子污染改变为原来的β倍，则只要在电子剂量场中考虑光子污染项，且该光子污染项为一个平均值，就可基本消除光子污染对电子剂量场的误差影响。β实质就是实际的光子污染对原计划中光子污染的调整值，用于获得更精确的电子剂量场。光子污染P'(z)为[14]，

经过以上的误差修正，这时所得的剂量场可认为是接近实际的剂量场：

其中，x-y-z符合式(4)。反演的具体目的是求得该模型中的修正因子，从而求得实际的全空间剂量分布。

反演法就是利用体外透射剂量曲线反演全空间的剂量分布。设深度为ZN处测得的透射剂量为D(X,Y,ZN)，反演的透射剂量为D'(X,Y,ZN)，这里采用均方根误差来衡量反演的计算精度：其中，m×n是用于反演的数据数目，N是反演的参数数目(二维N=5；三维N=8)。这样，该数学模型是一个非线性规划模型：目标函数是反演的透射剂量与测量透射剂量的最小均方根误差，约束条件是非负的光子污染倍数β，引入绝对值，该模型就是一个无约束的最优化问题(Minσ)：

这时反演的具体实施过程是用方程(9)左边的测量透射剂量拟合方程右边的若干参数，将算得参数代入式(9)，得到实际的剂量场分布。用非线性拟合法获得式(9)的参数求解。

最后，假设加速器相对稳定，在实施放疗中电子能谱的变化忽略不计，但在放疗前，通过反演法获得能谱与光子污染函数。

1.2 反演模型的实现

为进行对比验证，这里选用来自蒙特卡罗EGS模拟(即EGSnrc，其中野内统计误差小于0.5%)的测试例题对比数据。EGS模拟计算时，输入能谱与文献[14]电子能谱一致。同时考虑电子剂量场存在严重的光子污染问题，绝大部分光子污染来自入射电子源(约5%)，EGS程序无此计算功能。因此在反演模型中添加光子污染，数据来自 Deng J采用BEAMnrc模拟20 MeV Varian Clinac 1800电子加速器的光子污染[16]。

计算模型如图2所示，这里采用Mathematica6.0编程，实现求解如上的非线性规划模型；调用该软件自带软件包 Nonlinear Regression进行非线性反演，有7种反演算法可用，包括：C1：Levenberg-Marquardt；C2：Quasi-Newton；C3：Gradient；C4：Conjugate-Gradient；C5：Newton；C6：Principal-Axis；C7：NMinimize算法。计算机配置为Win XP，P-4 CPU：2.0 GHz，RAM：256 Mb。测试内容采用7个计算例题：可行性测试、射野大小影响、纵向材料不同、大偏差的影响及复杂几何模型下存在大测量误差的影响。其中，大偏差是指实际剂量场严重偏移原计划的剂量场；小偏差则是指这两者差别不大。采用的测试设计如下（坐标单位均为cm）：

图2 几何模型Fig.2 Geometric model.

(1) 可行性测试

例题 1水-空气-水(小偏差，10 cm ×10 cm野)

使用水体模大小为30 cm×30 cm×30 cm。图2为电子照射的几何坐标示意图，各电子剂量计算程序计算的结果的坐标要求与图2相一致。第一层为30 cm×30 cm×6 cm的水(ρ= 1 g/cm3，下同)；第二层为 30 cm×30 cm×2 cm 的空气(ρ= 0.00121 g/cm3，下同)；第三层为30 cm×30 cm×22 cm 的水。

(2) 射野大小影响

例题 2 水-空气-水(小偏差，20 cm ×20 cm 野)

例题1基础上，射野扩大为20 cm ×20 cm，其它不变。对比数据来源：EGS计算，参数与例题1一致。

(3) 纵向材料不同

例题 3 水-肺-空气-水(小偏差，10cm ×10 cm野)

例题1基础上，空气层上方30 cm×30 cm×6 cm的水模替换为30 cm×30 cm×2 cm水层与30 cm×30 cm×4 cm肺层(ρ= 0.258 g/cm3)，其它不变。对比数据来源：EGS计算，参数与例题1一致。

(4) 大偏差的影响

例题 4 水-空气-水(大偏差，10 cm ×10 cm 野)

例题1基础上，将入射电子束沿着坐标原点，与X-Z平面水平，逆时针转动5°，射野由射野坐标(–5,5)，(–5,5)，(5, –5)，(5,5)变化为(–4,5)，(–4,5)，(4.2, –5)，(4.2,5)；光子污染比例题1放大2倍，其他不变。

例题 5 水-空气-水(大偏差，20 cm ×20 cm 野)

几何模型在例题 2基础上，射野扩大为 20 cm×20 cm，然后将入射电子束沿坐标原点，与X-Z平面水平，顺时针转动 10°，射野由射野坐标(–10,10)，(–10,10)，(10, –10)，(10,10)变化为(–12,10)，(–12,10)，(11, –10)，(11,10)；光子污染比例题 3 放大5倍，其他不变。

例题6 水-肺-空气-水(大偏差，10 cm×10 cm野)

在例题水-肺-空气-水(小偏差，10 cm×10 cm 野)的基础上，将入射电子束沿着坐标原点与X-Z平面水平，顺时针转动 10°，射野由坐标(–5,5)，(–5,5)，(5,–5)，(5,5)变化为(–7,5)，(–7,5)，(6,–5)，(6,5)；光子污染比例题3放大5倍，其他不变。

(7) 复杂模型(横向与纵向材料不同，透射剂量在电子射程附近)

例题 7：水-空腔-水-空气-水(10 cm×10 cm 野)

例题1基础上，沿着中心轴自上而下，分别放置30 cm×30 cm×3 cm的水层，与2 cm×2 cm×2 cm的空气腔(周围是水)，30 cm×30 cm×3 cm的水，30 cm×30 cm×2 cm空气，30 cm×30 cm×20 cm水。

2 结果与讨论

反演结果见表1、图3与图4，这里只给出了典型结果的剂量场，图3与图4是例题4与例题7的离轴剂量曲线。其中，表1的“–”表示对应的算法不能实现计算。

除例题 3(但误差大)，其他例题的 NMinimize算法都不能实现，原因可能是该算法对约束条件不敏感，或该算法计算内存消耗太大，超过了Mathematica6.0程序包Nonlinear Regression所允许的计算要求。反演中能计算出射野大小、射线中心与偏斜角这些小偏差。在例题7中，重建的剂量场与EGS的对比表明，除了空腔处有一定的偏差外，其他部位重合得很好。空腔的偏差可能主要来自两方面，其一，EGS的计算网格稀疏带来的测量平面误差及对比 EGS计算曲线本身就有误差(0.5%–0.8%)；其二，反演模型对能量取平均的结果使得计算的尾部剂量场范围变短——即电子在尾部还没有完全衰减，但平均后就忽略了电子对尾部的剂量贡献，由此产生尾部误差（尾差）。因此改进的反演模型需要考虑改进对平均能量的近似。对于PrincipalAxis的失效可能与该算法的搜索方法有关，由于PrincipalAxis沿着均匀主轴有良好的收敛性，对于复杂模型就有可能失效。

表1 计算误差与时间Table 1 Calculated error and time

图3 水-空气-水(大偏差，10 cm×10 cm野)反演结果(a) 重建离轴剂量曲线(C1–C6算法)与EGSnrc计算结果(点线)对比，(b) 重建全空间等剂量线(C1–C6)，其中C1–C6算法的结果互相重合，C7算法无法实现Fig.3 Inversion results of water-air-water model (10 cm×10 cm field).(a) comparison between reconstructed off-axis curve (by C1–C6 algorithms) and simulation results of EGSnrc;(b) reconstructed results of the whole dose map (by C1–C6 algorithms), where the reconstructed results of C1–C6 but C7 algorithms agreed well with the EGSnrc's.

在这 7个例题中(表 1)，从计算精度来看，Levenberg-Marquardt、Quasi-Newton与 Newton算法具有共性。这是因为Levenberg-Marquardt算法是一种改进的高斯牛顿(Gauss-Newton)算法，也就是搜索最小均方的方法；Quasi-Newton算法又称为变尺度法，是牛顿法与梯度法的一种混合方法：如果取尺度为单位阵，此时搜索的方向就是梯度法方向；而进一步，当多次搜索完成时，其方向就是牛顿法方向。因此 Levenberg-Marquardt、Quasi-Newto算法都具有Newton算法的精度特性；但计算速度上，Levenberg-Marquardt算法>Quasi-Newton算法>Newton算法。其他算法在精度与速度上差别较大，精度上比不上前面三种方法，计算时间上各有差异。表1还表明，在所有的测试中Levenberg- Marquardt算法表现出很好的特性，精度快速，可作为电子剂量场反演的首选算法。

图4 水-空腔-水-空气-水(10 cm×10 cm野)模型反演结果(a) 重建离轴剂量曲线(C1–C5)与EGSnrc计算结果(点线)的对比，(b) 重建的全空间等剂量线(C1–C5)；C1–C5算法的结果互相重合，C6与C7算法无法实现Fig.4 Inversion results of water-cavity-water-air-water model (10 cm×10 cm field).(a) comparison between reconstructed off-axis curve (by C1–C5 algorithms) and simulation results of EGSnrc;(b) reconstructed results of the whole dose map (by C1–C5 algorithms), where the reconstructed results of C1–C5 but C6 and C7 algorithms agreed well with the EGSnrc's.

3 结论

本文建立了利用体外透射剂量重建全空间剂量分布的非线性规划模型。该模型考虑了电子放疗实施过程中的主要误差来源，包括射线平移、射线偏转、射野变形与光子污染。在原放疗计划电子剂量计算的基础上添加这些误差因子，利用体外透射剂量曲线反演出这些误差因子，从而获得经过修正的原计划，获得实际的全空间剂量分布。其中采用了7种算法包括Levenberg-Marquardt、Quasi-Newton、Gradient、Conjugate-Gradient、Newton、Principal-Axis与NMinimize算法来实现该模型。多个例题(非均匀介质的、与原计划有小偏差与大偏差等情况下的剂量场)测试结果表明，除了 NMinimize算法(只有一算例，误差大)不能实现该模型外，其他算法都能实现，反演误差小于0.9%；Levenberg-Marquardt算法表现出很好的特性，精度快速，可以作为电子剂量场反演的首选算法。进一步测试表明，利用接近电子射程的透射剂量会带来较大的反演误差（尾差），这可能是反演模型对能量取平均的结果使得计算的尾部剂量场范围变窄——即电子在尾部还没有完全衰减，但平均后就忽略了电子对尾部的剂量贡献。因此改进的反演模型需要考虑改进对平均能量的近似。

总之，每一种算法实现非线性规划反演模型都有特点，这为在实际实现中提供多方法的选择可能性；由于反演时间较短，精度高，可以实现实时剂量验证与实时剂量场可视化，可满足临床上实时的需要。

1 Wu Y, Song G, Cao R,et al. Chinese Physics C (HEP &NP), 2008, 32(Suppl. II): 177–182

2 吴宜灿, 李国丽, 陶声祥, 等. 中国医学物理学杂志,2005, 22(6): 683–690

WU Yican, LI Guoli, TAO Shengxiang,et al. Chin J Med Phys, 2005, 22(6): 683–690

3 Lin H, Wu Y, Chen Y. Phys Med Biol, 2006, 51: L13–L15

4 Tao S, Wu A, Wu Y. Clinical Oncology, 2006, 18(4):363–366

5 Tao S, Wu Y, Chen Y,et al. Computerized Medical Imaging and Graphics, 2006, 30(5): 273–278

6 International Commission on Radiation Units and Measurements. ICRU report 35. Radiation dosimetry:electron beams with energies between 1 and 50 MeV.Bethesda: ICRU, 1984

7 WU Yican, CU George. The Need for Further Development of CAD/MCNP Interface Codes. Presented at the American Nuclear Society Annual Meeting, Boston,Massachusetts, June 24–28, 2007; Transactions of American Nuclear Society, Vol.96, TRANSAO 961-882(2007), ISSN: 0003-018X. (invited presentation)

8 吴宜灿, 李莹, 卢磊, 等. 核科学与工程, 2006, 26(2):12–20

WU Yican, LI Ying, LU Rei,et al. Nucl Sci Eng, 2006,26(2): 12–20

9 HU Haimin, WU Yican, CHEN Mingliang,et al.Automatic Conversion of CAD Model into Neutronics Model. Presented at the 15th International Conference on Nuclear Engineering (ICONE-15), April 22–26, 2007,Nagoya, Japan, ISBN: 978-4-88898-159-0

10 Intensity Modulated Radiation Therapy Collaborative Working Group. Int J Radiat Oncol Biol Phys, 2001, 51(4):880–914

11 Luo Z. Radiat Phys Chem, 1998, 53: 305–327

12 Luo Jette, Walker. Med Phys, 1998, 25(10): 1958

13 GOU Chengjun, WU Zhangwen, LUO Zhengming,et al.Med Phys, 2003, 30: 415–423

14 LI Gui, LIN Hui, WU Aidong,et al. Chin Phys Lett, 2008,25(7): 2710–2713

15 LI Gui, WU Aidong, LIN Hui,et al. Electron energy spectrum reconstruction as nonlinear programming model using micro-adjusting algorithm, oral and full paper in the Seventh Asian-Pacific Conference on Medical and Biological Engineering (APCMBE), Beijing, 2008

16 George X Ding. Phys Med Biol, 2002, 47: 1025–1046

CLCR730.55, R811

Electron dose map inversion based on several algorithms

LI Gui1,2ZHENG Huaqing1,2WU Yican1,2,3FDS Team
1(Institute of Plasma Physics, Chinese Academy of Sciences, Hefei 230031, China)
2(Engineering Technology Research Center of Accurate Radiotherapy of AnHui Province, Hefei 230031, China)
3(School of Nuclear Science and Technology, University of Science and Technology of China, Hefei 230027, China)

The reconstruction to the electron dose map in radiation therapy was investigated by constructing the inversion model of electron dose map with different algorithms. The inversion model of electron dose map based on nonlinear programming was used, and this model was applied the penetration dose map to invert the total space one.The realization of this inversion model was by several inversion algorithms. The test results with seven samples show that except the NMinimize algorithm, which worked for just one sample, with great error, though, all the inversion algorithms could be realized to our inversion model rapidly and accurately. The Levenberg-Marquardt algorithm,having the greatest accuracy and speed, could be considered as the first choice in electron dose map inversion. Further tests show that more error would be created when the data close to the electron range was used (tail error). The tail error might be caused by the approximation of mean energy spectra, and this should be considered to improve the method. The time-saving and accurate algorithms could be used to achieve real-time dose map inversion. By selecting the best inversion algorithm, the clinical need in real-time dose verification can be satisfied.

Electron dose map, Inversion algorithm, Monte Carlo, Accurate radiation therapy

R730.55，R811

安徽省自然科学基金(090413095)资助

李 贵，男，1981年出生，在职博士，从事医学物理相关研究

2009-12-16，

2009-12-16