楊 磊，劉 燕，羅志平，李文革，周毅吉，張美琴

(1.煙臺大學，山東 煙臺 264005；2.湖南省農(nóng)業(yè)科學院，湖南 長沙 410300)

60Co γ輻照加工技術(shù)在農(nóng)業(yè)、工業(yè)、醫(yī)療衛(wèi)生、食品、環(huán)境保護等多個領(lǐng)域應用十分廣泛[1],目前總體發(fā)展平穩(wěn)，新建裝置呈現(xiàn)大型化、自動化及泛用化等特點[2]。因60Co源半衰期較短，為維持和擴大加工能力，需定期補充和更換60Co源棒，而源價格昂貴，為此提高源的整體利用率十分重要。生產(chǎn)實踐表明，優(yōu)化60Co輻照裝置排源設(shè)計是提高源利用率的核心技術(shù)因素，對優(yōu)化裝置劑量學性能、改善貨物加工劑量、降低單位成本、提高裝置加工量和經(jīng)濟效益等均有顯著影響。

排源方法主要包括人工經(jīng)驗排源和計算機自動排源兩大類，行業(yè)公司仍以使用前者為多，此類方法對工作人員專業(yè)素質(zhì)要求高，人工和時間花費多[3]，研究快速、可靠、結(jié)果方案足夠優(yōu)化的自動排源方法成為必然需要。目前國內(nèi)外有基于貪心算法、神經(jīng)網(wǎng)絡算法(深度學習算法)等的排源方法，貪心算法快速有效但全局搜索能力不夠[4]。李磊等[5]結(jié)合經(jīng)驗方法提出了一種混合求解方法，提高了全局搜索水平且排源效果較好。神經(jīng)網(wǎng)絡算法潛力巨大，但其難點在于需搜集足夠多的優(yōu)秀排源方案對算法進行預訓練[6]。文獻[7-8]提出了一種全局隨機搜索型方法，針對輻照裝置劑量不均勻度(DUR)進行優(yōu)化排源設(shè)計并證實其可行有效，文獻[9-10]在此基礎(chǔ)上對相關(guān)算法進行并行化改進并利用GPU加速，可滿足任意裝機規(guī)模輻照裝置的排源工作。近年來通過與多家單位的交流合作發(fā)現(xiàn)，企業(yè)不僅對DUR，而且對射線能量利用率、源棒排布工作量等更多目標均有需求，因此本文對原方法的評價數(shù)學模型、求解算法及代碼均進行改進，形成一種新型多目標優(yōu)化排源方法。

1 排源評價數(shù)學模型

影響60Co γ輻照裝置性能的主要因素及相互關(guān)系如圖1所示[11-13]，優(yōu)化選源、排源是保證和提高60Co γ輻照裝置加工質(zhì)量和加工效率的關(guān)鍵，因目前相關(guān)單位多直接購買成品源棒，源的活度分布選擇余地小，為此排源成為工作核心。相關(guān)單位在評判排源設(shè)計方案時會綜合考慮DUR、射線能量利用率及排布工作量等，而主流評價方法——參考平面法[14-15]難以直接評判貨箱中劑量分布DUR和射線能量利用率。本文從真實輻照情景出發(fā)，通過計算貨箱在輻照室中完成動態(tài)輻照后的累積劑量分布，可直接得到上述參量值。

1.1 輻射場分布和貨箱累積劑量分布的計算方法

(1)

圖1 影響60Co γ輻照裝置性能的主要因素及相互關(guān)系Fig.1 Main factor affecting 60Co gamma irradiation facility and interrelation

常數(shù)，對60Co，Γ=2.503×10-18C·m2·kg-1·Bq-1·s-1；A為源棒活度，Bq；L為源棒長度，m；a為待計算空間點到線源垂直距離的絕對值，m；b為垂足到線源最近點的距離，m；C為校準因子(測量值/理論值)，以修訂源棒包殼、劑量計、輻照裝置附屬機構(gòu)(如導源管、護源板(罩))帶來的射線散射和吸收影響，對已定裝置而言C近似為常數(shù)且歸一化處理后實際對當前研究無影響。

圖2 源架、貨箱及劑量計的空間布置Fig.2 Spatial arrangement of source rack, container and dosimeter

(2)

其中：oi為流水線層編號；si為流水線上貨箱工位編號。

進一步考慮貨箱換面照射，則貨箱中參考點的最終累積劑量Di,j,k可表示為：

(3)

1.2 基于多目標的評價函數(shù)定義

(4)

2 排源模型的局部復合求解算法

從評價模型可知，排源問題為典型的非線性規(guī)劃問題，相關(guān)研究表明模擬植物生長算法(PGSA)求解排源問題可行有效，該算法通用性強，原則上只需替換評價函數(shù)即可[16]。為減少現(xiàn)場布源工作量，本文改進了PGSA而形成局部復合求解算法。

2.1 局部復合求解的基本原理

PGSA潛在需至少1個在所有變量維數(shù)上盡量均勻分布的初始方案，文獻[16]表明強魯棒性雖使得初始方案對最終結(jié)果的優(yōu)化影響不大，但計算量變化劇烈。為保證初始方案的合理性及降低計算量，在求解過程中將源架分割為若干區(qū)域，源棒按總活度進行比例分配，分別排源后組成1個次優(yōu)化方案，再進行整體優(yōu)化計算。

2.2 算法的求解流程

算法的求解流程及功能如下。

1) 讀入計算所需的各種初始數(shù)據(jù)，包括源架、棒位、源棒、貨箱及填充材料、相應空間位置及控制算法本身等的系列參數(shù)。

2) 源架可按子源板自然分為N0個子區(qū)域，或按源架對稱軸劃分；將源棒按總活度對應分為N0組，接近退役的源棒優(yōu)先分配到一起。

3) 依次將源棒分組排布到各源板子區(qū)域，若全部排布完，則進行步驟11；若沒有，則進行步驟4。

5) 以xB為起點，以λ(≥1的正整數(shù))為步長，在各維變量上沿坐標軸的正負方向?qū)ふ覞M足約束條件[a,b]的生長點xTemp，若f(xTemp)

6) 求出X內(nèi)目標函數(shù)值最小點xmin，若f(xmin)≤f(Xmin)且不重復，則把該點放入最小生長點集合Xmin。

7) 若生長點集合X容量太大，則以一定的修剪比例對X內(nèi)的劣勢生長點進行裁減，經(jīng)驗值取0.4～0.6為佳；本輪產(chǎn)生的xTemp應全部保留從而保持生長點集合X的生物多樣性以防算法早熟而導致優(yōu)化停滯。

8) 若當前生長周期內(nèi)的f(xmin)相比f(Xmin)沒有進步即視為生長停滯，且停滯次數(shù)超過設(shè)定值，或總的生長周期數(shù)大于設(shè)定值，那么本分組計算結(jié)束，保存Xmin、f(Xmin)等結(jié)果待合并，進入步驟3開始下一分組排源，若不滿足則進入步驟9。

9) 計算X內(nèi)各生長點的形態(tài)素濃度，即各自的狀態(tài)空間概率[16]。

10) 產(chǎn)生隨機數(shù)β并計算落入的狀態(tài)空間，將所落區(qū)間對應的生長點記為新生長基點xB，轉(zhuǎn)回步驟5開始新的生長周期。

11) 將所有分組的Xmin合并，完成所有源棒在整體源架上的排放。

12) 若不需進一步整體排源，則輸出Xmin形成系列最終排源方案，進而精細計算得到這些方案在貨箱中的DUR、能量利用率、平均累積劑量等評價用結(jié)果。若進一步需要整體計算，則將Xmin作為整體排源用初始生長點X0，依次調(diào)用步驟4～10完成整體優(yōu)化計算，最終得到系列優(yōu)化排源方案供客戶選用。

2.3 基于GPU的并行加速計算

本文提出的新評價模型相比原模型，其理論復雜度更高而計算復雜度偏低，原模型所耗計算時間占總計算時間的90%以上[7]，新模型經(jīng)初步測試計算占45%左右，需進一步提升其計算效率。文獻[17]中基于GPU通用計算技術(shù)實現(xiàn)了PGSA的并行化，在費用增加不明顯的前提下使得原方法的計算效率顯著提高?；谕瑯铀悸?，改進后的PGSA也進行了基于GPU的并行化計算改進。

限于篇幅，本文僅給出當前算法中最重要的改進部分，即目標函數(shù)f的GPU并行化實現(xiàn)原理，如圖3所示，其中步驟1、2、4、5計算壓力小，以GPU最大Block和Thread數(shù)執(zhí)行即可，步驟3計算壓力大，執(zhí)行時設(shè)計為每個Block負責計算生長點集合1個生長點的目標函數(shù)值，Block中每個Thread計算1根源棒造成的劑量。

新評價模型對每個Thread的任務壓力明顯下降，求解代碼不僅計算效率顯著提升，且在GPU單精度計算模式下即可獲得滿意結(jié)果，可避免購買昂貴的專業(yè)雙精度計算卡。

圖3 PGSA中目標函數(shù)f基于GPU的并行化原理Fig.3 GPU parallelization principle of target function f by PGSA

3 實際應用及結(jié)果分析

3.1 相關(guān)裝置參數(shù)

以不同大型輻照裝置的實踐排源設(shè)計結(jié)果為例，驗證算法的可靠性和有效性，表1列出了輻照裝置的有關(guān)參數(shù)。

表1 輻照裝置參數(shù)Table 1 Parameter of irradiation facility

算法代碼的測試平臺主要參數(shù)為：CPU，Xeon E3 1230 V2(數(shù)量1)；GPU，GTX285(數(shù)量1，硬件規(guī)格1.3)；內(nèi)存，DDR3 1333(8 GB)；軟件平臺，Windows 7 64bit、VS2008、CUDA4.2、C++。

算法的魯棒性較好，參數(shù)限制少。計算時源板按每層板塊數(shù)劃分，源棒按總活度等比例劃分并由程序隨機抽取后均勻布放到棒位上，貨箱中參考點的空間間隔為5 cm，算法修剪比例為50%，編寫參數(shù)輸入文件，利用不同版本的代碼程序計算，其耗費時間情況列于表2。從表2可知，計算同等規(guī)模裝置，本文提出的新方法對應代碼的計算時間明顯縮短，可顯著縮短排源方案設(shè)計周期。

表2 不同版本排源代碼的計算時間Table 2 Computation time for source pencil deployment of various codes

3.2 結(jié)果及分析

現(xiàn)場部署排源方案后，按國標采用重鉻酸銀劑量計布置測試用貨箱，表3為測試貨箱參數(shù)、裝置劑量場的計算結(jié)果及實測結(jié)果。

從表3可知，計算采取均衡策略，在測試貨箱充滿模擬貨物前提下，對于裝置1和裝置2，新模型的理論DUR相比原模型分別增大了7.0%和16.5%，實測DUR相比計算結(jié)果分別增加了4.195%和0.70%，均顯著低于國標要求限值而實際影響不大，同時能量利用率顯著提高5.72倍和3.02倍。實測DUR相比理論值升高，其主要原因是計算模型相比實際輻照裝置做了合理簡化，未考慮護源板、貨箱壁，另外計算模型以空氣為受照介質(zhì)，而實測采用了密度更大的材料，射線的衰減和散射差別較大，IAEA證實這類問題可采用ρ-DUR關(guān)系曲線修訂而不影響排源[18]，且排源工作核心點是尋找優(yōu)化方案，DUR和能量利用率的值僅用于評判方案的優(yōu)化程度。

優(yōu)化的排源方案在裝置1、2參考面(圖2)上的等劑量線分布如圖4所示。

表3 測試貨箱參數(shù)、裝置劑量場的計算結(jié)果及實測結(jié)果Table 3 Test container parameter and result of facility dose field calculation and measurement

注：( )內(nèi)為模型改進后的數(shù)據(jù)；[ ]內(nèi)為相對偏差

圖4 裝置1(a)、2(b)參考面上的劑量分布Fig.4 Dose distribution along reference plane in facility 1 (a) and facility 2 (b)

從圖4可知，排源方案在裝置1、2參考面上的劑量分布很均勻，基本以裝置中心為圓心呈同心圓分布且梯度均勻，這對降低DUR有利，證實了排源方案設(shè)計的可靠性。

4 結(jié)論

本文提出了一種新型多目標優(yōu)化排源方法，建立了綜合考慮貨物DUR、射線能量利用率、現(xiàn)場布源工作量等多目標參數(shù)的新數(shù)學模型，設(shè)計了新型的GPU并行PGSA算法用于模型求解，并對兩座大型γ輻照裝置進行了實踐排源。計算和測量結(jié)果表明：新方法計算精確度得到保證，對昂貴專業(yè)計算卡的依賴度下降；計算效率提升顯著，計算時間從原有方法的h量級下降到min量級；新方法的理論和實測DUR比原方法增大但仍顯著低于國標控制標準，且能量利用率顯著提升，同時分區(qū)域排源降低了現(xiàn)場布源工作量。檢測結(jié)果驗證了新方法的可行性和有效性。