基于神經網絡的粒子輸運問題高效計算方法*

馬銳垚 王鑫2) 李樹 勇珩 上官丹驊?

1) (北京應用物理與計算數學研究所,北京 100094)

2) (中物院高性能數值模擬軟件中心,北京 100088)

蒙特卡羅方法是求解粒子輸運問題的有力工具之一,其局限性在于為達到精度要求需模擬大量粒子,計算耗時長,這阻礙了該方法的進一步應用,尤其在需快速響應的情形.本文結合神經網絡和若干蒙特卡羅方法基本原理發展了一種計算方法,能夠實現源分布可變,幾何、材料和目標計數不變的中子輸運問題的快速準確求解.首先,為高效生成用于神經網絡訓練的數據,利用重要性原理實現在同樣模擬次數基礎上有效擴充訓練數據集容量,在一定程度上克服了使用蒙特卡羅計算獲取訓練數據耗時長的缺點.進而,基于目標計數是源分布與重要性函數乘積積分的事實,設計了利用神經網絡實現快速輸運計算的策略.該網絡的輸入是中子源項,輸出是目標計數,在幾何、材料和目標計數固定的情況下,該神經網絡可重復使用,根據新的源項快速準確得到目標計數.本文所提出方法的原理和框架同樣適用于其他種類粒子的同類型輸運問題.基于若干基準模型的驗證表明,訓練得到的神經網絡能在不到1 s 的時間內得到目標計數,且與蒙特卡羅大樣本模擬得到基準結果的平均相對偏差均低于5%.

1 引言

蒙特卡羅(Monte Carlo,MC)方法和確定論方法是求解粒子輸運問題的兩類主要方法[1].隨著核技術的發展,各類新型核設施具有更復雜的幾何結構和反應機制,各向異性更強,同時對計算精度的要求也在不斷提高.為解決這類復雜粒子輸運問題,在計算機硬件性能不斷提高的基礎上,MC 方法由于能夠對物理過程和幾何構型進行精確建模且并行效率高,因此得到廣泛的發展和應用[2].同時,pin-by-pin 反應堆模型輸運[3]、復雜多物理場耦合計算[4]等大規模精細模型快速計算要求也在不斷推動MC 方法相關研究的開展[5].由于MC方法求解目標計數的相對誤差與樣本數平方根成反比,而計算時間與樣本數成正比,為使誤差達到要求,其模擬過程非常耗時[3,6-8].國內外研究團隊已發展多種減方差方法[9-19],目標是以較少的樣本數達到更低的統計誤差.權窗算法是其中重要的一類,一般認為比幾何重要性等算法效率更高,其算法參數可由重要性原理產生.重要性原理最早由Booth 和Hendricks[14]提出,網格的重要性定義為進入該網格的單位權重的粒子及其后代對計數器的總貢獻.雖然存在大量研究,但MC 輸運的絕對計算時間仍然相當可觀,且各種減方差方法的適應條件不同,在應用上有一定的局限性.

深度學習中的神經網絡是一種高度并行的信息處理系統,具有很強的自適應學習能力,能處理復雜的多輸入、多輸出非線性系統[20-23],訓練完成的網絡可根據輸入快速輸出目標值[24-27],在諸多領域有廣泛應用[28-32].近年來,在核相關領域利用數據驅動的神經網絡進行中子輸運計算的趨勢日益明顯.Berry 等[33]研究利用人工神經網絡和隨機森林分類器評估給定能群結構的適用性,以準確預測輕水堆中的中子倍增因子.Cao 等[34]提出基于神經網絡的重構中子場方法,用于重建反應堆壓力容器和堆芯區域內的中子場.Zhang 等[35]提出一種基于神經網絡的源項重建方法,實現復雜的三維源項快速重建.張海明等[36]建立了基于卷積神經網絡的深度學習模型,可用于直接預測反應堆堆芯有效增殖因子.針對多組分中子屏蔽材料優化設計中MC 模擬計算時間長的制約問題,林海鵬等[37]討論了利用神經網絡算法快速預測不同材料中子屏蔽效果的方法.Osborne 等[38]通過上采樣方法實現從少量粒子數MC 模擬結果到精確結果的映射.Kim 等[39]構建了一個以堆芯結構為輸入的神經網絡,可快速優化堆芯結構設計.與其他類型粒子輸運結合的研究,主要是在醫學圖像處理以及人體劑量學領域[40-42].

在上述數據驅動的神經網絡應用中,網絡通過分析所提供數據集的模式和輸入輸出之間的映射關系來學習和進行快速預測,因此訓練數據集的數量和質量對網絡的魯棒性和有效性有很大影響,而訓練數據的獲取通常采用MC 模擬或其他數值計算方法得到,難以避免大量計算時間和計算資源的消耗.對于數據驅動型的神經網絡,如何利用神經網絡代替特定問題的傳統計算,以及如何利用傳統方法高效獲取訓練數據都具有重要的研究意義.

本文的研究對象是源分布可變,幾何構型、材料和目標計數不變的問題,目標是在結果合理的前提下提高計算速度.這類問題在國防科技和國民經濟領域存在多種應用場景.例如,中子輸運是戰略國防科技研究所必須的多物理耦合計算中最消耗計算資源的一環,計算時間可以占總體時間的90%以上,在裝置變化緩慢階段,可認為每一步中子輸運的源項都有變化而裝置的幾何、材料變化可忽略不計,因此呈現為變源問題;在民用經濟方面,反應堆高置信度數值模擬需要快速計算壓力容器的中子注量率,而反應堆構型復雜、尺度大的特點導致其求解非常耗時,且工況復雜、燃料布置方案多樣,以堆芯整體為源項進行模擬計算時是典型的變源問題,面臨著快速求解難題;另外,變源問題的快速實時求解也是為“臟彈”恐怖襲擊的科學預防提供備選方案的重要技術手段.

本文的主要貢獻和創新之處在于: 1)以高效產生訓練數據為目標,從重要性的定義出發,首次將其應用在等效MC 計算的訓練樣本生成中,發展出神經網絡訓練數據高效生成方法,達到在同樣模擬次數基礎上有效擴充訓練數據集容量的效果.2)基于目標計數是源分布與重要性函數乘積積分的事實,構建可重復使用的高效神經網絡,以快速準確獲得各種源分布對應的目標計數,一定程度上突破了傳統MC 輸運計算的效率瓶頸.雖然本文的研究對象是中子,但該方法可適用于其他種類粒子的同類型輸運問題中.

2 方法介紹

本文建立了一種利用神經網絡加速MC 模擬中子輸運問題的方法(圖1 為技術路線圖),該神經網絡所必需的訓練數據集是通過對相應的輸運問題進行MC 模擬而獲得的.對于新的中子源分布,將其輸入給訓練好的神經網絡就能快速輸出相應的目標計數.數據驅動的神經網絡需要足夠的高質量訓練數據,但通過傳統的MC 模擬來獲取這些訓練數據非常耗時.針對這一問題,將重要性原理作為一種高效生成訓練數據的手段,由此在一定程度上提高了神經網絡的訓練效率.

圖1 本文技術路線圖Fig.1.Study framework of this paper.

2.1 理論基礎

對于中子固定源問題,統計體通量、能量沉積、劑量等物理量都與源項信息直接相關,其中源項信息包括位置分布、能量分布、粒子發射方向等多維信息.中子的屬性用p=(r,E,w) 表示,其中r表示粒子空間位置,E表示能量,w表示發射方向.中子源分布滿足歸一條件:

MC 模擬中目標計數可以表示為

其中,S(p) 為源分布,g(p) 為重要性函數.將相空間Ω進行網格劃分,即Ω=(M足夠大,表示劃分足夠精細),當源分布所在相空間劃分足夠精細的假設成立時,相空間內單個網格的重要性函數可以用該網格內一個點的重要性函數表示,則(2)式可被表示為

對于固定源問題,同一個模型(包括幾何及材料)和目標量對應同一套,對于不同的源分布,將S(p) 離散化為代入(4)式,即可得到計數I.

本文以體平均通量作為目標計數(相同方法可適用于任何計數).對于不考慮時間變量的固定源問題,定義

其中,Ψ(r,w,E) 是中子角通量密度,則體平均通量可以被表示為

ΦV可以通過徑跡長度估計法計算,粒子一旦通過計數區域就會對目標計數造成貢獻,一般而言較其他估計法效率更高,因此被多數MC 輸運程序采用.

2.2 基于神經網絡的粒子輸運問題高效計算方法

對于源分布多變,幾何、材料構型及目標計數穩定的一系列中子輸運問題,MC 模擬可看作一個黑箱,其輸入是源分布,輸出是固定的目標物理量,這個黑箱通過精確度很高,但是速度較慢的模擬將輸入轉換為輸出.通過設計合適的神經網絡并以合理的形式表示輸入端的源分布函數,在獲得充足的訓練數據基礎上,訓練得到的神經網絡將以很快的響應速度得到新的中子源分布對應的目標物理量.

在實際應用中,難以實現將連續形式的、不同分布的源項處理為神經網絡的統一輸入形式.為了解決此問題,本文首先對源分布進行N次采樣,從而產生N個源粒子.然后將這些粒子的六維屬性(r,E,w)作為點分布輸入網絡.最終網絡產生N個輸出的平均值被視為目標統計量的估計值.

為實現上述目標,將源分布在相空間離散化處理,進而以(4)式表達目標量后,可用訓練得到合適的神經網絡模型代替(4)式中權值,方法和網絡結構示意圖如圖2 所示.本文采用多層感知神經網絡(MLP)[21]的模型結構.神經網絡中單層感知器的擬合能力有限,而多層感知器的擬合能力更強,可用于解決復雜問題,因而MLP 的優勢在于它可以學習將任何輸入映射到輸出的權值.最典型的MLP 由三個網絡結構組成: 輸入層、隱藏層和輸出層.本文使用的網絡包括一個輸入層、多個全連接層和一個輸出層.在全連接層之間加入了整流線性單元 (ReLU)激活函數,激活函數為網絡引入了非線性特征,有助于網絡學習輸入和輸出之間的復雜關系,并能將神經元的輸出限制在一定范圍內.常用的激活函數包括Tanh,Sigmoid和ReLU,ReLU 的優點是能有效緩解梯度消失問題,從而提高訓練效率.各層之間還引入了一個dropout 層,以減少過擬合,提高神經網絡的泛化能力,使其更加穩健和精確.

圖2 用于代理加速MC 模擬的神經網絡結構Fig.2.Data-driven neural network for Monte Carlo simulation.

神經網絡需要大量高質量的訓練數據,以減少過擬合和網絡泛化性不足等問題.然而通過MC模擬生成如此大量的訓練樣本必定成本高昂.因此,本文提出了一種基于重要性原理[14]的方法來高效獲取足夠的訓練數據,此方法可在單次MC模擬中獲取盡可能多的訓練數據.重要性原理最早由Booth 和Hendricks[14]提出,用于輸運計算中的權重窗參數生成算法.在本文中,重要性不再是作為生成權重窗參數的基礎,而是作為相應網格中的源對目標計數的最終貢獻,從而擴大了從單次算例模擬中獲得的網絡訓練樣本的數量.具體說明如下: 將相空間劃分為網格時,網格i的重要性定義為通過該網格的粒子在單位權重下對目標計數的貢獻.

當Ω空間內網格劃分足夠細時,對Ωi內的點pi與其所在網格內定點距離足夠近,則有

上述方法在開源軟件OpenMC[43]中得以實現.對于每個網格,定義了兩個新的向量,分別記錄網格的總權重和總計數.對于每個粒子,都定義了一個新的向量來記錄粒子在整個軌跡中每一步所穿過的網格以及當前的權重.當前粒子被殺死時,記錄的網格和相應的權重會被累加到對應網格的總權重向量中.此外,如果粒子被判定已到達最終統計區域,當前粒子的權重將被累加到它(及其子代粒子)穿過的所有網格的計數向量中.當所有源粒子及其次級粒子的模擬結束后,根據(7)式計算所有網格的重要性.

3 結果與討論

3.1 重要性原理產生數據的正確性驗證

為驗證重要性結果計算的正確性,建立簡單幾何算例: 幾何為10 cm×10 cm×10 cm、材料為56Fe的正方體.源項為位于坐標(0,0,0)的中子點源,發射方向為沿X軸正方向的單一方向,能量為10 MeV 的單一能量.在幾何相空間劃分0.1 cm×0.1 cm×0.1 cm 的虛網格,統計目標區域探測器(1 cm×1 cm×1 cm)的通量計數.分別建立10 個探測器,其中心位于X軸的不同位置.本數值實驗基于相空間內一個點源的在進行MC 模擬的同時統計重要性結果,重要性結果取源項所在網格的重要性值,即以該網格為源對目標區域產生的目標值貢獻,驗證MC 計算結和重要性結果如圖3 所示.由重要性值的定義可知,本算例中重要性的統計過程為由該相空間網格內的此點源產生的樣本對目標區域的貢獻除以總粒子權重,而針對這一點源MC 模擬同樣為統計由此點源產生的目標結果除以總粒子數的平均值,因此對于此驗證算例重要性與MC 模擬結果完全一致,說明本研究在OpenMC中重要性數據產生方法實現的正確性.

圖3 重要性結果和MC 模擬結果對比(MC-MC 計算結果,IMP-重要性原理計算結果)Fig.3.Comparison of results of MC simulation and importance data (MC-results of Monte Carlo simulation,IMPresults of importance data).

3.2 神經網絡結果驗證

3.2.1 簡單模型

基于3.1 節所設置簡單算例對網絡進行訓練和測試.相空間網格的劃分參數如表1 所列.在相空間中隨機生成1000 個MC 算例(每個算例都是不同的點源),每個算例模擬106個粒子,以確保結果的最大統計誤差低于5%.同時,用重要性生成方法計算生成1000 套重要性數據及統計誤差,以重要性結果的統計誤差小于2%為標準選取重要性數據,使用傳統MC 模擬計算1000 個算例僅能產生1000 個源項及其目標統計量作為神經網絡的樣本,即訓練集容量為1000,使用本文提出的基于重要性原理的高效訓練數據產生方法后,在1000 個算例獲得MC 模擬結果的同時產生了24285 個重要性數據,均可以視為該網格中點源產生的目標統計量,即將1000 次模擬得到的訓練集容量提升了24 倍,顯著提升了單次模擬獲得的訓練樣本數.分別采用以下兩組數據作為數據集對神經網絡進行訓練和測試: 1)僅采用1000 個算例MC 模擬的結果;2)結合重要性數據和MC 模擬結果.預測通量與真值的相對偏差被用作損失,并按(9)式計算:

表1 簡單模型相空間劃分參數Table 1.Phase space meshing parameters of the Fe model.

其中,x是作為真值的MC 計算結果;y是網絡預測的通量結果.

神經網絡的構建和計算使用TensorFlow[44]實現.采用自適應學習率優化算法Adam[45]作為優化器,學習率為10-4.Adam 優化器的優勢在于無需手動調整,能根據每個參數的梯度和動量變化來調整其學習率,可以糾正其學習率消失、收斂過慢或是高方差的參數更新導致損失函數波動較大等問題.

圖4 顯示了網絡訓練過程中的損失函數曲線.在每個網絡的訓練測試中,80%的數據集用于訓練,10%用于驗證,10%用于測試.當基于驗證數據的損失開始增加時,即終止訓練.兩個網絡分別在5000 和10000 個算例時終止訓練.最終網絡輸出結果與MC 計算給出的真值的相對偏差分別收斂到19.75%和0.97%.這表明,采用重要性原理補充數據集可以使網絡更好地收斂.

圖4 網絡訓練損失曲線 (a) 使用 1000 個MC 模擬結果樣本進行訓練;(b) 結合重要性數據和MC 結果進行訓練Fig.4.Network training loss curves: (a) Trained with 1000 results samples of MC simulation;(b) trained with a combination of importance data and MC results.

對MC 模擬和神經網絡預測的時間耗費進行了比較.在MC 模擬中,模擬106個粒子以確保統計誤差保持在5%以下,使用1000 個實例的平均時間作為MC 模擬用.將1000 個算例輸入已訓練好的神經網絡,并計算平均預測時間作為神經網絡用時.結果表明,神經網絡將計算時間大約從70 s減少到0.3 s,這意味著計算用時減少了約200 倍,顯著提高了計算速度.

3.2.2 Kobayashi 模型

使用Kobayashi 模型[46]中的問題1 測試了本文所提出的方法.該模型包括三個區域: 源區域、真空區域和屏蔽區域(圖5).X=0 和Y=0 平面是該模型的兩個反射面.該例題給定單一截面,則僅在幾何和角度維度進行相空間劃分,劃分參數如表2 所列.

表2 Kobayashi 模型相空間劃分參數Table 2.Phase space division parameters of Kobayashi benchmark.

圖5 Kobayashi-1 模型幾何示意圖Fig.5.Diagram of Kobayashi-1 benchmark.

該 Kobayash 問題一的基準包括兩套不同截面設置: 在 Kobayash-1-i 中,源區和屏蔽區的材料是純吸收體;在 Kobayash-1-ii 中,源區和屏蔽區的材料設定為50%的散射截面.這種材料的總截面為0.1 cm-1,而空隙區域的總截面為10-4cm-1.對于每種截面,分別設置四個探測器位于坐標(15 cm,15 cm,15 cm),(25 cm,25 cm,25 cm),(35 cm,35 cm,35 cm)和(45 cm,45 cm,45 cm)處.MC 模擬使用107個源粒子,以確保所有探測器結果的最大統計誤差保持在5%以下.

針對不同的探測器訓練了8 個網絡,每個網絡的訓練各基于1000 個MC 算例.隨機抽取MC 和重要性數據的80%樣本用作訓練網絡的訓練集,10%樣本作為用于優化網絡參數的驗證集,10%的樣本作為最終測試訓練好的網絡預測效果的測試集合,按照(9)式分別計算訓練樣本、驗證樣本、測試樣本與真值的偏差作為訓練偏差、驗證偏差、測試偏差如表3 所列,測試偏差的分布如表4 和圖6 所示,由表3 和表4 可知所有情況的測試偏差均值均低于5%,小于5%的樣本占總測試樣本量的比例均超過90%,本數值實驗中的神經網絡輸入是相空間內給定能量、位置、方向六維坐標的單個離散點,對應到傳統MC 模擬中的單個樣本,而本方法利用神經網絡進行目標物理量的統計為抽樣所得源粒子對目標量貢獻的均值,因此預測均值偏差低于5%可以說明網絡給出的結果精度符合要求.偏差分布中最大偏差結果出現在源項距離探測器最遠的位置,分析原因是該位置的訓練數據本身統計誤差較大,使得訓練得到的網絡對該位置單個點源給出的偏差較大,這一特點與傳統MC 單個樣本模擬結果的特點相同.此外,與原始MC 模擬相比,獲得這些結果所需的時間幾乎可以忽略不計,用時從大約5 min 縮短到不到1 s,說明神經網絡對輸運計算效率的提升效果顯著.作為參考,Igor 等[47]使用離散縱標法在不同硬件平臺下計算完整Kobayashi 基準題的時間花費在數小時量級.

表3 Kobayashi 模型不同探測器的網絡預測結果Table 3.Prediction results for different detectors by networks of Kobayashi-1.

表4 Kobayashi-1 模型不同探測器的網絡預測偏差分布Table 4.Deviation distributions of predicted-results for different detectors by networks of Kobayashi-1.

圖6 Kobayashi-1 模型不同探測器的網絡預測偏差分布Fig.6.Deviation distributions of predicted-results for different detectors by networks of Kobayashi-1.

3.2.3 HBR2 反應堆基準題

H.B.Robinson-2 (HBR-2)核電站是Westinghouse 公司設計的商用壓水式輕水反應堆,于1971 年開始運行,歸Carolina Power and Light Company 所有.H.B.Robinson-2 壓力容器測定基準題[48]主要測量 HBR2 核電站壓力容器內部和外部的中子通量.

HBR-2 反應堆堆芯包含157 個燃料組件、堆芯外包殼、堆芯吊籃、隔熱罩、壓力容器和生物屏蔽等部件.在本研究中,在OpenMC 中構建了HBR-2 反應堆結構.統計以下兩個探測器(①和②)的通量為目標統計量,圖7 展示了 OpenMC 中HBR2 基準題的幾何結構及探測器示意圖,其中左半部分為水平面示意圖,右半部分為矢狀面示意圖.在OpenMC 中使用六層立方體來表示源項分布,每層由157 個立方體組成.

圖7 HBR-2 模型幾何結構示意圖Fig.7.Detector geometry diagram of HBR-2 benchmark.

幾何網格尺寸設定為 10.752 mm×10.752 mm×60.960 mm.源粒子發射方向為各向同性,能譜為瓦特譜.MC 模擬中模擬的粒子數量因源與探測器之間的距離而異,在源與探測器距離較遠的算例中,有必要模擬更多的粒子數量以確保獲得的通量統計誤差達到要求,因此每個算例都模擬了 4.0×106—6.0×108個粒子,以保證最大統計誤差保持在5%以下.最終利用MC 模擬共建立了942 個算例來生成真值數據集,其中754 個算例被用作為訓練樣本,92 個算例被用作為驗證樣本,92 個算例被用作為測試樣本,按照(9)式分別計算測試樣本、驗證樣本、預測樣本與真值的偏差作為訓練偏差、驗證偏差、測試偏差(表5).

表5 HBR-2 模型不同探測器的網絡預測結果Table 5.Prediction results for different location detector networks of HBR-2 benchmark.

表5 列出了網絡對這兩個探測器的預測結果,測試偏差的分布見表6 和圖8.網絡預測結果的平均相對偏差分別為3.98%和4.42%,小于5%的樣本占總測試樣本量的比例均超過90%,同前算例一樣,最大偏差結果出現在靠近堆芯區域即源項距離探測器最遠的位置,由于該位置的訓練數據本身統計誤差相對較大而使得單個預測結果與真值偏差超過5%.考慮到目標物理量的統計過程為從源分布中多次抽樣后輸入神經網絡并計算對目標量貢獻的均值,因此預測均值偏差低于5%表明網絡輸出結果精度符合要求.在時間成本方面,相對于MC模擬將計算時間從數小時縮短至約0.3 s,實現了顯著提速.作為參考,Roberto[49]在不同硬件平臺下使用TORT 程序基于離散縱標法對HBR-2 基準題進行計算,單CPU 時間花費在數十小時量級.可見本文提出方法可以給出滿足工程精度標準的預測結果,而計算時間與原始MC 模擬相比可以忽略不計,從而獲得顯著的效率提升.

表6 HBR-2 模型不同探測器的網絡預測偏差分布Table 6.Deviation distributions of predicted-results for different detectors by networks of HBR-2.

圖8 HBR-2 模型不同探測器的網絡預測偏差分布Fig.8.Deviation distributions of predicted-results for different detectors by networks of HBR-2.

4 結論

本文設計了用于加速中子輸運問題MC 模擬的神經網絡,只需一次性訓練即可針對不同源項反復使用.在訓練該神經網絡時,使用重要性原理來高效獲取足夠多的高質量訓練數據,從而一定程度上克服了數據生成成本高的難題.上述方法通過Kobayashi 和HBR-2 等基準問題進行了驗證.網絡預測結果的相對偏差始終低于5%.而且與原始MC 模擬相比,獲得這些結果所需的計算時間可忽略不計.

本文研究的一個局限性是只適用于源分布有變化,而幾何、材料及目標計數保持不變的情形.如果幾何、材料或目標計數發生變化,則需要重新訓練網絡,后續研究計劃將幾何、材料以及目標計數的變化也作為輸入信息引入網絡訓練過程,以增強網絡的泛化能力,并以遷移學習的訓練方法提高訓練效率.

感謝北京應用物理與計算數學研究所勇珩研究員牽頭的“ AI++”團隊成員的討論和幫助.