中子及次級γ在高空長距離蒙特卡羅輸運模擬中的減方差方法

劉 利,左應(yīng)紅,牛勝利,朱金輝,商 鵬,李夏至

(西北核技術(shù)研究所,西安 710024)

高空核爆炸核輻射環(huán)境研究[1-3]對低軌道衛(wèi)星、臨近空間(20～100 km)飛行器及通信傳播[4]等都有著重要意義。模擬中子及次級γ在高空大氣中的大空間長距離輸運時，輻射源與記錄點的距離可達幾十公里，甚至上百公里，使到達記錄點附近的粒子概率極小。以50 km為例，真空中各向同性點源向外出射的粒子落在50 km外1 m2區(qū)域內(nèi)的概率約為5×10-10。同時中子由高空向低層大氣輸運時，大氣密度迅速增加，中子平均自由程迅速減小，中子與大氣的多次散射作用使粒子到達記錄點的概率進一步減小。高度為20 km時，1 MeV中子和γ在水平方向的自由程約為1.1 km和1.4 km，遠小于輻射源與記錄點的距離。當(dāng)記錄點附近的粒子數(shù)太少時，蒙特卡羅方法模擬結(jié)果的相對偏差較大，結(jié)果不可信。對瞬發(fā)中子、γ和X射線在高空中輸運模擬已有研究[4-7]，但針對的是空氣質(zhì)量厚度較小及大氣散射作用較弱的情況，針對中子及次級γ在高空長距離輸運至低層(如20 km處)這種空氣質(zhì)量厚度較大的大氣環(huán)境的研究還未見報道。中子及次級γ在空氣質(zhì)量厚度較大的大氣環(huán)境中的輸運模擬計算對研究核電磁脈沖環(huán)境和通信傳播等有重要意義。本文開展了中子在高空大氣中的長距離輸運模擬的減方差方法研究，建立了中子在高空非均勻大氣輸運的幾何模型，采用MCNP程序，比較了源偏倚、幾何分裂、截斷、權(quán)窗、強迫碰撞和DXTRAN球等減方差方法[8-11]計算結(jié)果的相對偏差與計算效率，并提出源偏倚、幾何分裂和DXTRAN球組合方法為最佳減方差方法。

1 中子輸運模擬幾何模型

圖1給出了中子在高空大氣中長距離輸運模擬的幾何模型。高空大氣是典型的非均勻連續(xù)介質(zhì)，密度隨高度變化非常劇烈。模型將大氣結(jié)構(gòu)進行分層處理，用于描述大氣密度隨高度的變化。高度越高大氣密度越小，中子平均自由程越大，分層處理時大氣厚度也越大。根據(jù)大氣密度的特點與計算精度需要，將高度為0～50 km的大氣分為50層，每層厚度為1 km；將高度為50～110 km的大氣分為6層，每層厚度為10 km。每層大氣的密度根據(jù)美國標(biāo)準(zhǔn)大氣模式取值[12]。地平面以下5 km為土壤層，密度為2.31 g·cm-3。

設(shè)置中子輻射源為位于高度為80 km爆心處的各向同性點源，中子權(quán)重設(shè)置為1，中子能譜取典型氫彈出殼中子能譜[3]。記錄點位于爆心正下方。由于中子飛行到20～40 km高度時才會與大氣發(fā)生明顯的相互作用，記錄點的高度取為20～40 km。爆心與最低記錄點的距離為60 km，且中子在大氣中輸運時受到多次散射作用，使直接模擬的粒子能到達最低記錄點的概率非常小。

設(shè)置模擬源粒子數(shù)目N為1×107。在20～40 km高度范圍內(nèi)每2 km設(shè)置一個記錄點，采用點計數(shù)方法統(tǒng)計記錄點的中子及次級γ注量。圖2為直接模擬給出的中子與γ的注量、相對偏差σ和FOM因子隨高度的變化關(guān)系。FOM因子可表示為

(1)

其中，t為計算時間，min。

FOM因子可反映出模擬結(jié)果的計算效率，F(xiàn)OM因子越大，減方差效果越好。由圖2(a)可見，中子注量隨記錄點高度降低而降低，γ注量隨記錄點高度變化幅度相對較小。由圖2(b)可見，中子注量和γ注量的相對偏差隨記錄點高度降低而增加，F(xiàn)OM因子隨高度降低而降低。當(dāng)記錄點高度小于30 km時，中子注量和γ注量的相對標(biāo)準(zhǔn)偏差均大于5%，結(jié)果不可信。即使模擬粒子數(shù)逐漸增加至1×108，直接模擬得到的中子和γ注量依然沒有收斂。圖3為直接模擬和采用減方差(variance reduction，VR)方法計算給出的20 km高度處中子與γ的注量、相對偏差和FOM因子隨粒子數(shù)的變化關(guān)系。由圖3可見，直接模擬的中子與γ注量、相對偏差和FOM因子隨粒子數(shù)的增加而大幅度波動；采用減方差方法(圖3中為源偏倚+幾何分裂+DXTRAN球)可有效降低統(tǒng)計漲落，減小相對偏差，提高FOM因子。

2 減方差方法

2.1 單一減方差方法

常用的蒙特卡羅模擬減方差方法可分為3類：(1)改進抽樣方法，包括源偏倚、強迫碰撞、指數(shù)變換、重要性抽樣和統(tǒng)計估計抽樣等；(2)總體分布控制法，包括權(quán)窗、分裂與輪盤賭技巧、截斷、分層處理和兩步法等；(3)半解析法，包括點探測器、環(huán)探測器和DXTRAN球等。本文以20 km高度處中子和γ注量計算優(yōu)化為目標(biāo)，開展中子及次級γ在高空中長距離輸運的減方差方法研究。本文采用的減方差方法主要有：

2.1.1 源偏倚

對粒子源的方向和能量等參數(shù)做偏倚，增加重要區(qū)域的抽樣粒子數(shù)，并相應(yīng)減小粒子的權(quán)重。在圖1所示的輸運幾何模型中，爆心高度處中子平均自由程很大，由爆心向上發(fā)射的粒子對爆心正下方的記錄點幾乎不產(chǎn)生計數(shù)貢獻。為提高抽樣效率，采用源方向偏倚提高粒子到達記錄點的概率。根據(jù)模型幾何結(jié)構(gòu)，設(shè)置源偏倚使源出射方向與垂直地面向下方向的夾角余弦在0.948 7～1范圍內(nèi)和范圍外的概率各為50%，相應(yīng)修改粒子的權(quán)重以保證模擬的無偏性。夾角余弦為0.948 7時，對應(yīng)的夾角為18.4°，可覆蓋幾何模型中以20 km高度處記錄點為中心，水平方向距離為0～20 km的范圍。

2.1.2 幾何分裂與輪盤賭

將模型劃分為多個柵元，設(shè)置重要柵元的重要性較高，不重要柵元的重要性較低。當(dāng)粒子進入重要性高的柵元內(nèi)，有一定的概率分裂為多個粒子；當(dāng)粒子進入重要性低的柵元，則只有一定的概率存活，從而增加重要性較高的柵元內(nèi)的模擬粒子數(shù)。對圖1所示幾何模型，設(shè)置高度為40 km以上的柵元重要性為1，40 km以下的柵元重要性每2層翻一倍，直到柵元重要性達到28；高度為20 km以下的柵元重要性每2層降低至一半，直到柵元重要性為1。

2.1.3 強迫碰撞

在指定的柵元內(nèi)人為增加粒子碰撞抽樣概率，并將進入柵元的粒子分為碰撞與不碰撞2部分。不碰撞部分直接離開柵元，粒子權(quán)重相應(yīng)減少。碰撞部分具有粒子剩余權(quán)重，并使用輪盤賭來控制粒子數(shù)量，防止出現(xiàn)大量低權(quán)重粒子而增加計算時間。針對圖1所示幾何模型，設(shè)置粒子在高度為18 km以上的柵元內(nèi)發(fā)生強迫碰撞，18 km以下時，不發(fā)生強迫碰撞，從而引導(dǎo)粒子進入20 km高度處的記錄點附近。

2.1.4 指數(shù)變換

在特定方向上減小微觀截面并在相反方向上增加微觀截面，從變換后的非模擬密度函數(shù)抽樣碰撞距離，相應(yīng)修改粒子權(quán)重，從而增加指向特定方向的模擬粒子數(shù)。修改后的截面為

(2)

2.1.5 DXTRAN球

使用DXTRAN球時，設(shè)置一包圍感興趣區(qū)域的小球。在球外的每次碰撞都能產(chǎn)生特殊的DXTRAN粒子，按照確定的概率進入DXTRAN球內(nèi)，權(quán)重利用確定論方法計算。碰撞產(chǎn)生的粒子正常輸運，只是在進入DXTRAN球范圍時被終止模擬，保證結(jié)果的無偏性。對圖1所示幾何模型，設(shè)置DXTRAN球的球心位于20 km高度處，半徑取1.5 km。

2.1.6 權(quán)窗

權(quán)窗是與位置和能量相關(guān)的分裂和輪盤賭技巧，是最常用的減方差方法之一。設(shè)置圓柱形柵格權(quán)窗，在高度方向上按圖1所示模型劃分，半徑方向上網(wǎng)格寬度逐漸增加。對中子與γ聯(lián)合輸運，同時使用中子與γ的權(quán)窗，并取20 km處的γ注量來產(chǎn)生權(quán)窗。設(shè)置模擬源粒子數(shù)目為1×108個，表1為單獨采用某一減方差方法的計算效率。表1第3行還給出了采用模擬粒子數(shù)為1×1010時的直接模擬結(jié)果，用于驗證其他減方差方法的無偏性。當(dāng)模擬粒子數(shù)增加至1×1010時，計算相對偏差有所降低，但表征計算效率的FOM因子也大幅度降低。

表1 單獨采用某一減方差方法的計算效率Tab.1 Computational efficiency by using a single VR method

由表1可知，與指數(shù)源偏倚方法相比，分布式源偏倚方法將FOM因子分別從4.9和0.2提高至23和13，計算偏差分別從6.4%和32%降低至3.0%和4.0%，所以本文采用分布式源偏倚方法；與直接模擬相比，采用源偏倚、幾何分裂和DXTRAN球方法計算的相對偏差都有所降低，同時FOM因子有所提高，這3種方法都引導(dǎo)粒子向記錄點位置偏倚，大大提高了計數(shù)區(qū)域內(nèi)的粒子抽樣概率，從而大幅度提高了計算效率，同時降低了計算偏差。

由表1還可知，與直接模擬相比，DXTRAN球方法的FOM因子提高最多，由4.6和3.2提高至52和39；幾何分裂方法的相對標(biāo)準(zhǔn)偏差降低最多，由9.1%和10.8%降低至0.9%和1.3%；強迫碰撞方法的相對偏差有所降低，但FOM因子同樣降低；指數(shù)變換方法的相對偏差較大，且FOM因子偏低，減方差效果較差，是因為收斂性太差，修改拉伸系數(shù)沒有使減方差效果明顯好轉(zhuǎn)，與其他減方差方法結(jié)合使用可能會提高減方差效果；采用權(quán)窗進行一次迭代的計算結(jié)果將中子注量計算的FOM因子提高了2個量級，相對偏差降低至1%，但計算結(jié)果明顯偏離正常計算值，與采用模擬粒子數(shù)為1×1010時直接模擬的結(jié)果相差較大，且γ注量計算相對偏差較大；采用權(quán)窗進行二次迭代的中子與γ的計算相對偏差都較大，且FOM因子較低，這說明一次迭代計算中子注量的相對偏差為偽偏差，計算結(jié)果已完全不可信。

2.2 減方差方法的組合優(yōu)化

組合使用多種減方差方法可綜合利用各種減方差方法的優(yōu)點，通常比單獨采用一種減方差方法具有更高的計算效率[10]。對于圖1所示幾何模型，各向同性發(fā)射的源粒子輸運至記錄點位置的概率極低，采用方向源偏倚方法使源粒子更多地朝指定方向發(fā)射。表2為采用源偏倚方法和其他一種減方差方法組合的計算效率。由表2可知：源偏倚與幾何分裂或指數(shù)變換組合的方法大幅度提高了計算效率；源偏倚與DXTRAN球組合的方法比單獨采用DXTRAN球方法的減方差效果略有提高；源偏倚與DXTRAN球方法通過調(diào)整粒子方向或產(chǎn)生DXTRAN粒子都使粒子朝記錄點運動，與源偏倚與幾何分裂等組合方法相比，源偏倚與DXTRAN球組合的方法計算效果較差一點；源偏倚與強迫碰撞組合的方法有效減小了相對偏差，但也降低了FOM因子。

表2 源偏倚與其他一種減方差方法組合的計算效率Tab.2 Computational efficiency by using source bias method and another VR methods

由表2還可知，源偏倚與權(quán)窗組合方法給出的相對偏差和FOM因子隨著迭代次數(shù)的增加并未收斂，反而出現(xiàn)強烈的波動，即使某一步迭代的中子注量相對偏差較小，中子注量計算值也與其他方法相差較大，同時對應(yīng)的γ注量相對偏差較大，相應(yīng)的權(quán)窗下限也出現(xiàn)了明顯的波動，且有時會出現(xiàn)較多的權(quán)窗下限為“0”的權(quán)窗；源偏倚+幾何分裂+權(quán)窗、源偏倚+DXTRAN球+權(quán)窗和SGD(source bias+geometry splitting+DXTRAN)+權(quán)窗等組合方法的計算結(jié)果同樣無法收斂且相對偏差較大。在高空中子及次級γ長距離輸運模擬中，權(quán)窗下限中出現(xiàn)大范圍的“0”權(quán)窗下限，使模擬粒子過度朝向非“0”權(quán)窗下限柵元偏倚，導(dǎo)致權(quán)窗無法準(zhǔn)確描述粒子對探測器計數(shù)的貢獻。權(quán)窗不適用的原因可能有2個：一是空間尺度太大，較遠柵元內(nèi)的模擬粒子對探測器計數(shù)的貢獻被忽略，從而使權(quán)窗下限為“0”，額外計算表明，在采用權(quán)窗方法模擬高度為10 km范圍內(nèi)的高空中子輸運可有效減小相對偏差；二是高空大氣密度變化劇烈，高度為20～80 km時，空氣密度變化3～4個量級，高密度區(qū)域內(nèi)權(quán)窗網(wǎng)格尺寸應(yīng)該比低密度權(quán)窗網(wǎng)格尺寸要小得多，但MCNP自帶的網(wǎng)格權(quán)窗參數(shù)難以同時滿足高密度與低密度區(qū)域權(quán)窗網(wǎng)格尺寸的需求。此類“0”權(quán)窗下限問題在不少深穿透問題的研究中都有出現(xiàn)[13]，有關(guān)文獻還提出采用源偏倚、指數(shù)變換、密度迭代和猜測初始權(quán)窗等方法來獲得合理的權(quán)窗參數(shù)，但對于百公里量級非均勻大氣中中子及次級γ長距離輸運問題，采用以上方法和本文所述方法均未能得到合理可靠的結(jié)果，還需采用幾何迭代、自適應(yīng)權(quán)窗和全局權(quán)窗等新的權(quán)窗生成方法進一步研究。

表3為采用源偏倚與其他多種減方差方法組合的計算效率。結(jié)合表1、表2和表3可知，采用強迫碰撞或在其他方法上加入強迫碰撞可有效降低計算結(jié)果的相對偏差，但同時FOM因子也明顯降低。這是因為對于圖1所示幾何模型，為引導(dǎo)粒子輸運至記錄點，需在相鄰柵元內(nèi)設(shè)置強迫碰撞，導(dǎo)致產(chǎn)生大量低權(quán)重粒子，大大增加了計算時間；在源偏倚+幾何分裂、源偏倚+DXTRAN球和SGD組合方法上再加上指數(shù)變換的減方差效果并不理想。結(jié)合表1、表2和表3還可知，SGD組合方法的減方差效果最好，中子與γ的FOM因子分別由直接模擬法的4.6和3.2提高至165和41，分別提高了39倍和13倍，相對偏差由9.1%和10.8%降至0.2%和0.5%。SGD組合方法與模擬粒子數(shù)為1×1010時直接模擬的結(jié)果較為接近，相對差異在模擬結(jié)果的相對偏差范圍之內(nèi)，驗證了該組合方法的無偏性。圖3中的減方差方法就是SGD組合方法。由圖3可見，SGD組合方法的模擬結(jié)果收斂性較好，隨著粒子數(shù)的增加，注量的計算值基本保持不變，且相對偏差和FOM因子變化幅度較小。

表3 源偏倚與其他多種減方差方法組合的計算效率Tab.3 Computational efficiency by using source bias method and other multiple VR methods

采用SGD組合方法計算給出的中子與γ的注量、相對偏差和FOM因子隨高度的變化關(guān)系，如圖5所示。由圖5(b)可見，高度為20 km處記錄點的相對偏差比22 km處還要小，但FOM因子比22 km處數(shù)據(jù)大，是因為在20 km處設(shè)置了DXTRAN球，提高了計算效率與計算精度。采用該組合方法模擬得到由高空輸運至20～40 km高度范圍內(nèi)的中子和γ注量的相對偏差小于1%。若將模擬粒子數(shù)由1×108降低至1×107，計算結(jié)果的相對偏差能保持在5%以內(nèi)。由此可見，采用SGD組合方法也適用于模型中其他高度記錄點的中子和γ注量計算。

3 時間譜的計算

在高空核輻射環(huán)境與效應(yīng)研究中，中子或γ的時間譜是極其重要的參數(shù)[3]。采用不同減方差方法模擬計算得到的高度為20 km處，中子與次級γ的注量率隨時間的變化關(guān)系，如圖6所示。由圖6可見，SGD組合方法計算得到的中子及次級γ的時間譜較平滑，波動較小；直接模擬得到的時間譜波動較大，且收斂性較差。這說明組合方法同樣適用于計算中子及次級γ時間譜。

由圖6還可見，SGD組合方法計算得到的中子次級γ時間譜的相對偏差在γ峰值處較小，在γ峰值后逐漸增加。不同的時間段內(nèi)γ時間譜的相對偏差不一致。針對這種情況，本文采用了時間分裂的減方差方法。與幾何分裂和能量分裂類似，時間分裂通過分裂與輪盤賭來增加關(guān)心時間段內(nèi)的模擬粒子數(shù)，降低計算結(jié)果的相對偏差。設(shè)置模擬時間達到2×10-4s時，粒子被一分為二；時間達到8×10-4s時，粒子以0.5的概率存活；時間達到2×10-3s和5×10-3s時，粒子再次被一分為二。

由圖6還可見，采用SGD+時間分裂組合方法計算的γ時間譜首末兩端更加平滑，相對偏差也有所降低。表4列出了不同減方差方法計算時間譜得到的相對偏差σ。其中σ(>2×10-3s)指的是時間點大于2×10-3s的相對偏差。由表4可知，采用減方差方法后，中子及次級γ時間譜的相對偏差明顯減小，在此基礎(chǔ)上加入時間分裂算法可進一步減小相對偏差，但計算時間略有增加。將模擬粒子數(shù)設(shè)置為1×109時，SGD+時間分裂組合方法計算中子及次級γ時間譜的相對偏差小于5%。

表4 不同減方差方法計算結(jié)果的相對偏差Tab.4 The relative deviation of time spectra of neutron and secondary γ by using different VR methods

4 結(jié)論

針對中子及次級γ在高空長距離蒙特卡羅輸運模擬的效率不高問題，研究了源偏倚、幾何分裂、強迫碰撞、權(quán)窗、指數(shù)變換、DXTRAN球和時間分裂等減方差方法在高空長距離輸運中的應(yīng)用，通過比較各個減方差方法的相對偏差與計算效率，得出SGD組合方法的減方差效果最好，與直接模擬相比，SGD組合方法給出的中子及次級γ的FOM因子分別提高了39倍和13倍，得到高度為20～40 km時的中子及次級γ注量的相對標(biāo)準(zhǔn)偏差小于1%。在該方法基礎(chǔ)上加上合理設(shè)置參數(shù)的時間分裂算法可進一步降低時間譜的相對偏差。