基于模擬退火算法的加速器能譜重建

郭金森GUO Jin-sen；李傳龍LI Chuan-long；王薇WANG Wei

（中國原子能科學研究院，北京102413）

0 引言

在利用加速器打靶的韌致輻射光子束進行治療時，劑量計算軟件為加速器TPS 治療計劃系統的基礎。任何劑量計算程序都需要用到光子能譜數據。然而在實際的治療系統中，有些情況下很難直接測量能譜數據。因此針對醫用加速器，不同的研究者一直在探求獲取準確能譜的方法[1]。主流方法有：一是利用蒙特拉羅方法模擬加速器治療機頭得到能譜，如Deng J 等人[2]在2000 年利用EGS4 程序模擬了來自瓦里安Clinac2100c 和2300c/d 加速器的4、6和15MV 光子束能譜。2015 年Fátima Padilla-Cabal 等人[3]利用MCNPX 和EGS 程序建立了醫科達加速器的精確模型并做了劑量計算。這種方法能夠保證所建能譜的精度，但前提是需要精確地知道準直器、均整器、電離室等部件的尺寸、材料等數據。二是利用射束穿過不同介質的透射數據得到X 射線能譜，典型如[4，5]等人的工作，該方法并不適用于高能光子射束。三是在射束中軸測量PDD 數據得到精確測量的前提下[6]，利用數學方法進行能譜重建。如2001 年，Deng J[7]等人在研究了從電子束中心軸百分比深度劑量（PDD）曲線，采取“隨機蠕變”（random creep）算法推導了臨床電子和光子能譜，但該方法有可能收斂至局部最優解。2010 年北京大學姚杏紅[8]等人采用Cimmino 方法重建瓦里安600CD 直線加速器X 射線能譜，該方法基于雙源模型，結果比較精確，但大射野情況下重建能譜出現了病態解。2013 年陳元華[9]等人針對瓦里安15MV 光子束利用遺傳算法優化重建了光子能譜，得到的能譜與蒙特卡洛模擬得到的能譜具有較好的一致性。2013 年劉娟[10]等人根據利用模擬退火算法重建了西門子醫用加速器6MV X射線能譜。

模擬退火算法為化工、冶金工業術語，現在已經作為最優化算法用于各個行業。模擬退火算法的核心是在爬山算法的基礎上增加Metropolis 準則，即從當前狀態i→新狀態j 時，即使新狀態j 在某個規則下優于狀態i，也以一定的概率接受狀態j，否則仍保留狀態i。模擬退火算法有一定的概率跳出局部最優解從而找到全局最優解。本文擬利用模擬退火算法，結合實驗PDD 數據，對國內某型醫用加速器6MV、18MV 能量下X 射線能譜進行重建，重建結果和蒙特卡羅模擬得到的能譜進行對比以驗證可靠性。并為后續劑量計算程序開發奠定基礎。

1 材料與方法

1.1 問題的數學描述

由測量PDD 和單能光子PDD 求解光子能譜的問題，可以描述為解線性方程組問題。假設加速器出射能譜可以離散化為E1，E2…En，其中E1~En為每個離散化能量區間的平均能量。n 組能量的按照注量加權的權重為ω1，ω2，…ωn，則應有如下如下線性方程組成立：

其中dij代表深度為i 處，Ei能量下的單能PDD 的劑量數據。Dj為深度為Hj處不同權重下的單能PDD 加權相加后的合成PDD 數據，在ω1，ω2，…ωn為實際能譜情況下，D1…Dm即為測量PDD 數據。該方程組一般為超定方程組，其解不止一個。因此需要找到最符合物理意義的一組解，即為能譜數據。可見，最優化的能譜可以使合成PDD 數據最接近測量PDD，也即相似度最高。為描述合成PDD（用量F 表示）和測量PDD（用M 表示）x 相似度，引入相關系數：

進行能譜重建。其中f（x）為目標函數，Fi為深度為Hi處的合成PDD，Mi為深度為Hi處的測量PDD。

1.2 單能光子PDD 以及能譜模擬

為獲取（1）中的單能PDD，采用蒙特卡羅模擬程序BEAMnrc 的劑量計算程序xyznrc 進行模擬。模體為30cm×30cm×30cm 厚的水模，密度為1.0g/cm3。源皮距為100cm，射野為5cm×5cm。沿射野中心軸取1.5cm×1.5cm×0.1cm 體素共299 個。針對6MV 計算了17 個能量點的單能PDD 數據，針對18MV 則計算了23 個能量點，即公式（1）中n=23。所有模擬結果誤差最大均不超過0.4%。

為驗證重建能譜的準確性，利用蒙特卡羅模擬程序EGSnrc 建立了加速器治療頭模型，獲取了治療頭6MV 和18MV 的出射能譜。加速器治療頭結構包括靶、初級準直器、均整器、電離室，次級準直器（上、下兩個）等結構。通過控制上下次級準直器控制射野大小為5cm×5cm，分別模擬了6MeV 與18MeV 窄電子束入射情況下的韌致輻射譜。見圖2，圖3。

1.3 利用模擬退火算法重建光子能譜

在熱力學上，一塊被加熱至高溫的物體的降溫過程被稱之為“退火”。退火過程滿足Metropolis 準則，即溫度為T時，出現能量差為dE 的降溫概率P（E）為：

即溫度越高，出現一次能量差為dE 的降溫的概率就越大；溫度越低，相應概率就越小。

模擬退火算法的基本流程圖見圖1。

圖1 模擬退火算法流程圖

模擬退火算法基本流程為：

①對光子各能量箱賦予一定的權重，經驗表明該權重不可過分偏離實際權重，否則可能導致收斂過慢甚至收斂至病態解。

②開始迭代過程，每一次迭代在每個能量箱上加一個小的隨機量，即

其中，η 為攝動系數，取0.00001，εi取[-1，1]之間的隨機數。

③如果新的目標函數f（x）（見公式（3））小與當前目標函數，則接受新解和新的目標函數。如果新目標函數大于當前目標函數，則以概率P（C）=e-dC/T接受新解和新的目標函數。容易看出隨著T 逐漸降低，越難接受“壞”的新解。迭代終止條件為目標函數低于某一個截止值，或模擬退火溫度低于某個截止值。

模擬退火降溫過程由初始溫度T 及溫度控制參數a表示，如下

這里初始溫度取1×10-9，a 取0.999，每個溫度出迭代100 次，取目標函數截止值Scut為4.6×10-5，即合成PDD 和測量PDD 之間的相關系數C=1-Scut即0.99954。

2 結果

圖2、圖3 分別為優化前后6MV、18MV 能譜和MC 模擬能譜對比圖。二者均作了歸一化處理。可以看出優化能譜和MC 模擬能譜峰位相同，譜形基本一致。計算得到相關系數均為0.99 以上。另外可以看出，圖2 優化后能譜尾端有微小抬升，懷疑此處算法陷入局部最優解。但低能光子由于權重小對總劑量貢獻較低。而高能部分兩條能譜基本一致。圖3 優化后能譜和MC 模擬能譜除個別點外也基本一致。

圖2 優化后6MV 能譜和MC 模擬能譜對比圖

圖3 優化后18MV 能譜和MC 模擬能譜對比圖

圖4-圖5 為優化前后6MV、18MV PDD 數據相對于測量數據的偏差。可以看出優化后絕大部分點相對誤差在0.5%以下。

圖4 優化前后后6MV PDD 相對測量數據的偏差

圖5 優化前后18MV PDD 相對測量數據的偏差

3 討論

本研究利用模擬退火算法，基于蒙特卡羅模擬的單能光子PDD 數據和測量得到的PDD，對醫用電子加速器6MV、18MV X 射線能譜進行重建。計算得到的能譜與蒙卡程序直接模擬治療頭得到的能譜形狀基本一致。

模擬退火算法具有全局搜索性，相比簡單的爬山算法，有更大的概率得到全局最優解。需要指出的是，選擇合適的初值仍然是有必要的，首先可以減小搜索時間，其次進一步降低收斂至局部最優解的風險。今后的工作將分為兩步進行：①調整算法，嘗試在迭代計算初期以更大的步長進行搜索，在目標函數f（x）滿足一定條件時進行小步長精細化搜索，提高算法魯棒性；②利用建成能譜開發卷積劑量計算程序。