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

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

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

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

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

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

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

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

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

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

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

1 優(yōu)化方法

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

1.1 程序架構(gòu)

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

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

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

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

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

輸出文件處理模塊主要功能是處理蒙特卡羅程序結(jié)果文件，得到對應(yīng)計數(shù)卡的結(jié)果，根據(jù)適應(yīng)度公式計算適應(yīng)度并返回給適應(yīng)度計算模塊，同時產(chǎn)生記錄文件。在輸出文件處理模塊中，通過調(diào)整每個參數(shù)的權(quán)重可以針對某個參數(shù)進行優(yōu)化。記錄文件中的信息包括蒙特卡羅程序輸出文件中的原始數(shù)據(jù)、權(quán)重因子和適應(yīng)度等參數(shù)。

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

1.2 適應(yīng)度函數(shù)

遺傳算法中適應(yīng)度函數(shù)是決定結(jié)果的最關(guān)鍵因素，是描述個體優(yōu)劣的主要指標(biāo)，遺傳算法根據(jù)個體適應(yīng)度大小進行優(yōu)勝劣汰[16]。根據(jù)優(yōu)化的問題，定義適應(yīng)度函數(shù)為：

(1)

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

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

表1 IAEA推薦的參數(shù)Table 1 IAEA recommended parameter

2 BSA優(yōu)化設(shè)計

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

圖2 優(yōu)化前相關(guān)結(jié)構(gòu)Fig.2 Related structure before optimization

2.1 材料種類優(yōu)化

如前文所述，必須使用適當(dāng)?shù)穆牧蠈a(chǎn)生中子慢化到超熱范圍。這些材料必須在超熱能段具有低散射截面，但在較高的能量段具有高散射截面。某些具有較高吸收截面的慢化材料在中子到達患者病灶之前會吸收較多中子，發(fā)生(n，γ)反應(yīng)，導(dǎo)致γ射線污染。但是，傳統(tǒng)方法很難給出滿足BSA全部要求的材料組合。

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

利用程序進行優(yōu)化，得到每代個體適應(yīng)度的變化，如圖3所示。

圖3 適應(yīng)度隨每一代的變化Fig.3 Fitness change with each generation

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

表2 材料種類優(yōu)化后結(jié)果Table 2 Result after material type optimization

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

2.2 材料順序優(yōu)化

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

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

表3 材料順序優(yōu)化后的結(jié)果Table 3 Result after material sequence optimization

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

2.3 幾何尺寸優(yōu)化

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

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

圖5 優(yōu)化后相關(guān)結(jié)構(gòu)Fig.5 Related structure after optimization

表4 幾何尺寸優(yōu)化后的結(jié)果Table 4 Result after geometry optimization

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

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

圖7 BSA優(yōu)化后的中子能譜Fig.7 Neutron spectrum after BSA optimization

3 總結(jié)和展望

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