遺傳算法優化加速器硼中子俘獲治療束流整形設計

程凡杰，李開健，張 巍,*，楊歷軍，李 巖

(1.中國原子能科學研究院，北京 102413；2.北京中科天馬信息技術有限公司，北京 100089)

硼中子俘獲治療(BNCT)是一種二元治療癌癥的方法，經過幾十年的研究，被認為是目前治療癌癥的有效方法之一。其原理是將含硼藥物注入人體后，含硼藥物會在腫瘤區域富集，然后利用中子束照射腫瘤部位，腫瘤中的硼原子與中子反應后產生高傳能線密度粒子α和7Li，這兩種粒子在細胞中的射程分別為9 μm和5 μm，粒子產生的電離輻射損傷作用僅限于含10B的細胞，并且可以引起腫瘤細胞DNA不可修復的損傷，從而達到選擇性地殺死癌細胞的目的。

BNCT治療一般采用超熱中子進行治療，超熱中子束穿透能力強，可以治療深部腫瘤，是各國科學家研究的重點[1]。IAEA-TECDOC-1223報告[2]給出了BNCT裝置出口束流參數，常被作為束流整形裝置(BSA)優化的參考標準。

早期BNCT使用反應堆做中子源，反應堆中子源具有中子注量高、穩定性好的特點，但由于造價高和核安全要求嚴等問題，難以在醫院普及使用。相比于反應堆中子源，加速器中子源具有安裝簡單、維護費用低、相對安全的優點，更容易在醫院中普及。近年來，世界各國都在致力于研制基于加速器的硼中子俘獲治療(AB-BNCT)設備。

加速器打靶產生的中子能量較高，而且含有大量γ無法直接用于治療，必須經過特定的系統進行處理，才能滿足BNCT要求，這一系統即為BSA。BSA可能的設計方案多，設計變量主要有材料種類、材料順序以及幾何尺寸。

傳統的BSA設計方法需要在設計之前比較不同種類的材料在慢化、準直、反射以及屏蔽γ射線等方面的能力[3-9]，以找到不同區域合適的材料，選出材料后還需要進行材料的組合和厚度的調整。但是由于材料的散射、吸收等各種反應截面隨能量是變化的，當中子穿過材料時散射、吸收同時都在發生，不同尺寸材料效果也不同。這導致優化BSA時各參數會相互影響，如在降低快中子成分時，往往需要快中子截面大的材料，結果也會導致超熱中子注量率下降。因此，傳統方法在不同區域推薦的材料也只能是定性的結論，而且這些結論往往只對某一能量區間的中子有效，對于優化設計作用較小。由于材料種類較多、尺寸可調整的范圍巨大，同時調整尺寸和材料得到的方案幾乎是無限的，這個過程將花費大量的時間。得到初步方案后，再通過理論計算驗證設計方案的可行性，如果不滿足要求，則針對性地提出優化方案。

傳統優化方法本質上是從少數個體迭代尋求最優解，這容易陷入局部最優解，難以兼顧多個參數同時優化。利用傳統方法進行BSA設計工作量大、設計周期長，優化效果并不好。

為解決AB-BNCT BSA優化設計問題，將BSA設計轉化為確定各柵元的材料種類、順序以及幾何尺寸的問題，優化目標是達到IAEA推薦的4個參數的要求。BSA設計轉化為多個變量、多個目標參數的優化問題。像這種問題可通過智能算法來解決，遺傳算法就是其中的一種。

遺傳算法是一種模仿生物進化機制發展起來的隨機全局搜索和優化方法，它借鑒了遺傳、突變、自然選擇以及雜交等進化生物學中的一些思想，是解決搜索問題的一種通用算法。該算法對于被求解問題要求低，只需要合適的適應度函數。該算法從初始解的1個集合開始搜索，同時處理群體中的多個個體，按照多種路線進行平行搜索，并且采用輪盤賭的方法進行選擇。遺傳算法的優勢是搜索范圍大，有利于全局擇優。遺傳算法還具有魯棒性、可并行、可移植性的特點，遺傳算法在核電廠燃料管理[10-12]、屏蔽設計[13-15]以及反應堆BNCT BSA設計[16]等方面均有應用。

本文基于中國原子能科學研究院14 MeV質子回旋加速器生成中子的特性，研制遺傳算法與蒙特卡羅程序相耦合的設計程序，對BSA的材料種類、材料順序、幾何尺寸進行優化設計，以得到滿足AB-BNCT治療要求的超熱中子束。

1 優化方法

BSA設計的本質是一個帶有約束條件的多目標尋優過程。一種優化方法是將慢化區進行精細劃分為若干個區，然后用多種材料進填充。這種方法從原理上可行，但是會導致可能的解過多，計算量巨大，難以得到最優解。另外一種方法是先將慢化區初步劃分為較少數量的區域，通過優化算法選出最可能的材料；然后使用可能的材料進行材料順序的優化；最后調整每個柵元的尺寸。這種BSA設計方法轉化為材料種類、材料排列順序、材料尺寸的優化，利用遺傳算法可較快得到最優解，提高設計效率。

1.1 程序架構

遺傳算法主要通過編碼和解碼、交叉、變異、適應度評估來進行尋優，根據遺傳算法原理編制了遺傳算法與蒙特卡羅方法耦合的設計程序，程序流程圖如圖1所示。該程序主要包括4個模塊：遺傳算法模塊、適應度計算模塊、輸入文件生成模塊和輸出文件處理模塊。

圖1 程序流程圖Fig.1 Program flow chart

遺傳算法模塊主要功能是設置參數，參數包括變量數目，變量上限、下限，種群規模，交叉概率，變異概率，初始種群，保留到下一代的精英數量，終止條件等。返回的是最優個體的基因及其適應度。

適應度計算模塊主要功能是調用輸入文件生成模塊，執行蒙特卡羅程序以及調用輸出文件處理模塊，得到適應度。該模塊在產生輸入文件之前判斷是否存在該文件輸入文件、輸出文件，如果存在就跳過輸入文件生成和蒙特卡羅程序計算兩個步驟，直接處理輸出文件得到適應度。由于蒙特卡羅程序計算耗費了絕大多數時間，這樣做的優勢是可以避免重復計算、節約時間。

輸入文件生成模塊主要功能是將個體基因帶入到模板文件中產生蒙特卡羅程序計算的輸入文件。參數個數不受限制，可以通過合適的編碼實現材料尺寸和材料種類同時優化，但這樣會導致可能的解更多，優化時間更長。實際優化時按照材料種類、順序和尺寸分步進行優化，每一步變量數目較少，這樣可以提高效率。

輸出文件處理模塊主要功能是處理蒙特卡羅程序結果文件，得到對應計數卡的結果，根據適應度公式計算適應度并返回給適應度計算模塊，同時產生記錄文件。在輸出文件處理模塊中，通過調整每個參數的權重可以針對某個參數進行優化。記錄文件中的信息包括蒙特卡羅程序輸出文件中的原始數據、權重因子和適應度等參數。

該程序中斷運行后可以繼續運行，運行結束后的蒙特卡羅程序輸出文件和記錄文件保留在文件夾中。初始種群可以根據記錄文件個體信息選擇加快收斂速度。

1.2 適應度函數

遺傳算法中適應度函數是決定結果的最關鍵因素，是描述個體優劣的主要指標，遺傳算法根據個體適應度大小進行優勝劣汰[16]。根據優化的問題，定義適應度函數為：

(1)

式中：fitness為目標函數；wi為權重因子；S為解空間；k為優化的參數個數；Ti為每個目標函數的目標值，計算時采用IAEA的推薦值；x為向量，進行材料優化時，x為代表材料順序的向量，進行尺寸優化時，x為代表材料尺寸的向量；fi(x)代表方案一定情況下利用蒙特卡羅程序計算得到的各參數值。

IAEA的推薦值如表1所列，其中：En為超熱中子能量范圍；df為快中子產生的劑量；dr為γ產生的劑量；J為中子流量；φepi為超熱中子注量率，可以看出各目標參數值有數量級的差異，該適應度函數可消除參數之間數量級的差異，實現多個目標參數同時優化。

表1 IAEA推薦的參數Table 1 IAEA recommended parameter

2 BSA優化設計

結合14 MeV回旋加速器生成中子的特性，按照裝置輻射防護要求，確定整型裝置內慢化材料為含硼聚乙烯，屏蔽材料為鉛。然后將慢化區域分成9個區域，準直區分為2個區域，慢化區中每個區厚度約為8 cm，優化前相關部分結構如圖2所示。然后，利用程序依次進行材料種類、材料順序、幾何尺寸的優化。

圖2 優化前相關結構Fig.2 Related structure before optimization

2.1 材料種類優化

如前文所述，必須使用適當的慢化材料將產生中子慢化到超熱范圍。這些材料必須在超熱能段具有低散射截面，但在較高的能量段具有高散射截面。某些具有較高吸收截面的慢化材料在中子到達患者病灶之前會吸收較多中子，發生(n，γ)反應，導致γ射線污染。但是，傳統方法很難給出滿足BSA全部要求的材料組合。

BSA常用的材料有鋁、Fluental、氟化鋁、鎳、鉍、氟化鋰、氟化鎂、硫、鉛、含硼聚乙烯、鈹、聚乙烯、不銹鋼、石墨、氟化鈣等[3-9]，共計15種。這些材料多數是反應堆或加速器BNCT慢化材料，也有屏蔽材料，這里不進行人為區分，全部作為備選材料。優化時對超熱中子注量率、快中子成分、γ成分、流量注量率比值4個參數的權重wi分別進行設定。

利用程序進行優化，得到每代個體適應度的變化，如圖3所示。

圖3 適應度隨每一代的變化Fig.3 Fitness change with each generation

根據優化結果，得到最優的10種方案均含有鋁、氟化鋁、石墨、鉛、Fluental，其中最優方案的出口束流的參數列于表2。由表2可看出，各項參數較靶后的參數均有很大提升，但此方案的超熱中子注量率較低，γ成分偏高，不能滿足要求，需要繼續優化。

表2 材料種類優化后結果Table 2 Result after material type optimization

此次優化得到了合適的材料，與最初材料相比可看出含氫和硼的材料均被舍棄。通過分析記錄文件可看出，方案中含有這兩種元素會大幅度降低超熱中子注量率，在優化過程中含硼聚乙烯、聚乙烯等材料很容易被淘汰，這與實事相符。

2.2 材料順序優化

根據初步優化選出的這些材料進行材料順序的優化。由于超熱中子注量率低，快中子和γ成分高，需要提高兩參數的絕對權重。進行優化時，每一代適應度的變化如圖4所示。由圖4可看出，在200代左右優化停止，得到最優的材料排列順序，BSA第2次優化結果列于表3。

圖4 材料順序優化適應度隨每一代的變化Fig.4 Fitness change with each generation for material sequence optimization

表3 材料順序優化后的結果Table 3 Result after material sequence optimization

Fluental的主要成分是27Al和19F，第2次優化后BSA主體元素排序基本上是208Pb、27Al和19F，在出口處用的主要是208Pb和12C。經過分析，208Pb、27Al、19F 3種核素共振峰所在的能量區間逐漸降低，其彈性散射截面也逐漸下降，這種排序可以使能量高于10 keV的中子都能得到慢化，大幅度降低了快中子成分，而在出口處采用208Pb可以減少γ成分。

2.3 幾何尺寸優化

進行材料種類和材料順序優化時，每個區域幾何尺寸是固定的，兩次優化后快中子成分仍然偏高，其他參數滿足IAEA要求。然后優化各層之間的尺寸，直到滿足參數要求。

優化后相關結構如圖5所示。調整后的方案與之前方案相比減少了氟化鋁厚度，增加了鉛和Fluental的厚度，減少了出口鉛的厚度。這是由于鉛和Fluental在相應的能量范圍慢化性能優于氟化鋁，這樣調整可進一步減少快中子份額。由于屏蔽γ所需要的鉛比較少，減少鉛厚度可以減少超熱中子的損失。幾何尺寸優化后的出口束流參數列于表4。

圖5 優化后相關結構Fig.5 Related structure after optimization

表4 幾何尺寸優化后的結果Table 4 Result after geometry optimization

另外模擬了加速器打靶產生的中子能譜和BSA出口處的中子能譜，如圖6所示。優化后在0.5 eV～10 keV能量范圍內，超熱中子份額約為84.2%，如圖7所示。對比圖6可看出，利用該方法優化BSA設計可以大幅度提高超熱中子份額，減小了快中子份額。

圖6 加速器打靶產生的中子能譜Fig.6 Neutron spectrum behind cyclotron target

圖7 BSA優化后的中子能譜Fig.7 Neutron spectrum after BSA optimization

3 總結和展望

本文將遺傳算法與蒙特卡羅程序進行耦合，并利用該程序開展了AB-BNCT BSA優化設計，包括材料種類、材料順序和幾何尺寸。優化前只需要修改生成蒙特卡羅計算輸入文件的模板和各變量的上下限，優化過程自動完成。計算結果表明，利用該程序優化得到的BSA模型滿足IAEA的推薦參數要求，該方法為AB-BNCT BSA設計提供了一種有效的技術方案。該程序通用性強，也可應用于反應堆設計、核臨界安全和屏蔽設計中。