基于多算法的電子劑量場反演方法

李 貴 鄭華慶 吳宜燦,3 FDS團隊

1（中國科學院等離子體物理研究所 合肥 230031）

2（安徽省精確放療工程技術研究中心 合肥 230031）

3（中國科學技術大學核科學技術學院 合肥 230027）

基于多算法的電子劑量場反演方法

李 貴1,2鄭華慶1,2吳宜燦1,2,3FDS團隊

1（中國科學院等離子體物理研究所 合肥 230031）

2（安徽省精確放療工程技術研究中心 合肥 230031）

3（中國科學技術大學核科學技術學院 合肥 230027）

為了研究精確放射治療中多算法的電子劑量場反演方法及反演效果，建立非線性規(guī)劃電子劑量場反演模型，利用體外透射劑量重建全空間電子劑量場，并通過多種算法實現。采用 7個例題測試結果表明：除NMinimize算法不能實現該模型外，其他算法都能實現，其中Levenberg-Marquardt算法最精確快速。進一步測試表明，利用接近電子射程的透射劑量會帶來較大的反演誤差（尾差），原因可能是能譜平均帶來的，因此改進方法需考慮改進對平均能量的近似。但由于反演時間較短，精度高，可實現實時劑量場反演；并且每一種算法都有特點，為實現提供多方法的選擇性，滿足臨床上對實時劑量驗證的需要。

電子劑量場，劑量驗證，能譜，精確放射治療

精確放射治療中，蒙特卡羅模擬等精確劑量計算可達到很高的正向計算精度[1?3]，但實際治療操作中的不可控誤差[4,5]，會使照射區(qū)域的劑量分布偏移計劃劑量分布，因此需實施實時的劑量場的重建即劑量場反演，以確保患者體內腫瘤區(qū)域的劑量符合計劃要求。目前，精確放療劑量驗證的方法基本考慮光子照射，但也須驗證電子束照射的劑量，方能確保“精確放療”的實現。此外，電子劑量場反演還可應用于工業(yè)探測，如用高能電子束照射工件，測量透射劑量來反推工件內部劑量場分布，從而獲得工件的損傷與材料分布等情況。

電子束照射劑量驗證的難點，一是電子射程短，射程之外當然無電子劑量分布；二是光子污染問題，光子劑量場可以很好地控制電子污染，而電子劑量場中的光子不能很好地消除。第一個問題的解決方案，只有提高電子束能量，如20 MeV電子束能穿透10 cm水模，可穿透人體一些部位（如乳腺癌的治療），可實施劑量驗證，因此本文把電子劑量場限定在高能區(qū)。對于第二個問題，本文將把光子場也考慮在內，建立一個混合場模型。

本文發(fā)展一種新型的電子劑量場的反演方法，在混合筆形束模型基礎上，建立非線性規(guī)劃反演模型并利用多種非線性反演方法實現快速精確的電子束劑量場實時反演，最后與蒙特卡羅模擬對比[6?10]進行驗證。

1 基于混合筆形束的非線性規(guī)劃反演模型

1.1 劑量場誤差分析與反演模型的建立

Luo等[11?13]根據 Fermi-Eyges多級散射理論，發(fā)展了混合筆形束模型，用于精確放射治療的精確解析電子劑量計算，其模型決定電子光子混合劑量場的參數有：電子能譜、矩形野尺寸、散射參數、光子污染函數與“距離平方反比修正”函數中的若干系數。本文將在此模型基礎上建立電子劑量場反演模型。

電子劑量場的反演如圖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) 射線轉動。射線發(fā)生轉動等同于坐標系轉動。根據坐標轉動公式，設θ和φ角分別是法線方向n的極偏移角和方位偏移角，則舊坐標(x,y,z)和新坐標(X,Y,Z) (實驗坐標系)的變換關系為

(4) 光子污染變化。我們曾在先前的電子能譜反演研究中提出[14,15]：光子污染由兩個指數函數描述。假設高能散射光平均衰減系數μp、軔致輻射光子和次級光子衰減系數μe、散射光子與軔致輻射光子的比例系數ν等決定光子物理性質的參數不變，并設電子劑量場的重新分布使光子污染改變?yōu)樵瓉淼摩卤?，則只要在電子劑量場中考慮光子污染項，且該光子污染項為一個平均值，就可基本消除光子污染對電子劑量場的誤差影響。β實質就是實際的光子污染對原計劃中光子污染的調整值，用于獲得更精確的電子劑量場。光子污染P'(z)為[14]，

經過以上的誤差修正，這時所得的劑量場可認為是接近實際的劑量場：

其中，x-y-z符合式(4)。反演的具體目的是求得該模型中的修正因子，從而求得實際的全空間劑量分布。

反演法就是利用體外透射劑量曲線反演全空間的劑量分布。設深度為ZN處測得的透射劑量為D(X,Y,ZN)，反演的透射劑量為D'(X,Y,ZN)，這里采用均方根誤差來衡量反演的計算精度：其中，m×n是用于反演的數據數目，N是反演的參數數目(二維N=5；三維N=8)。這樣，該數學模型是一個非線性規(guī)劃模型：目標函數是反演的透射劑量與測量透射劑量的最小均方根誤差，約束條件是非負的光子污染倍數β，引入絕對值，該模型就是一個無約束的最優(yōu)化問題(Minσ)：

這時反演的具體實施過程是用方程(9)左邊的測量透射劑量擬合方程右邊的若干參數，將算得參數代入式(9)，得到實際的劑量場分布。用非線性擬合法獲得式(9)的參數求解。

最后，假設加速器相對穩(wěn)定，在實施放療中電子能譜的變化忽略不計，但在放療前，通過反演法獲得能譜與光子污染函數。

1.2 反演模型的實現

為進行對比驗證，這里選用來自蒙特卡羅EGS模擬(即EGSnrc，其中野內統(tǒng)計誤差小于0.5%)的測試例題對比數據。EGS模擬計算時，輸入能譜與文獻[14]電子能譜一致。同時考慮電子劑量場存在嚴重的光子污染問題，絕大部分光子污染來自入射電子源(約5%)，EGS程序無此計算功能。因此在反演模型中添加光子污染，數據來自 Deng J采用BEAMnrc模擬20 MeV Varian Clinac 1800電子加速器的光子污染[16]。

計算模型如圖2所示，這里采用Mathematica6.0編程，實現求解如上的非線性規(guī)劃模型；調用該軟件自帶軟件包 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 結論

本文建立了利用體外透射劑量重建全空間劑量分布的非線性規(guī)劃模型。該模型考慮了電子放療實施過程中的主要誤差來源，包括射線平移、射線偏轉、射野變形與光子污染。在原放療計劃電子劑量計算的基礎上添加這些誤差因子，利用體外透射劑量曲線反演出這些誤差因子，從而獲得經過修正的原計劃，獲得實際的全空間劑量分布。其中采用了7種算法包括Levenberg-Marquardt、Quasi-Newton、Gradient、Conjugate-Gradient、Newton、Principal-Axis與NMinimize算法來實現該模型。多個例題(非均勻介質的、與原計劃有小偏差與大偏差等情況下的劑量場)測試結果表明，除了 NMinimize算法(只有一算例，誤差大)不能實現該模型外，其他算法都能實現，反演誤差小于0.9%；Levenberg-Marquardt算法表現出很好的特性，精度快速，可以作為電子劑量場反演的首選算法。進一步測試表明，利用接近電子射程的透射劑量會帶來較大的反演誤差（尾差），這可能是反演模型對能量取平均的結果使得計算的尾部劑量場范圍變窄——即電子在尾部還沒有完全衰減，但平均后就忽略了電子對尾部的劑量貢獻。因此改進的反演模型需要考慮改進對平均能量的近似。

總之，每一種算法實現非線性規(guī)劃反演模型都有特點，這為在實際實現中提供多方法的選擇可能性；由于反演時間較短，精度高，可以實現實時劑量驗證與實時劑量場可視化，可滿足臨床上實時的需要。

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

2 吳宜燦, 李國麗, 陶聲祥, 等. 中國醫(yī)學物理學雜志,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年出生，在職博士，從事醫(yī)學物理相關研究

2009-12-16，

2009-12-16