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

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

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

1 中子輸運(yùn)模擬幾何模型

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

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

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

(1)

其中，t為計(jì)算時(shí)間，min。

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

2 減方差方法

2.1 單一減方差方法

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

2.1.1 源偏倚

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

2.1.2 幾何分裂與輪盤(pán)賭

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

2.1.3 強(qiáng)迫碰撞

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

2.1.4 指數(shù)變換

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

(2)

2.1.5 DXTRAN球

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

2.1.6 權(quán)窗

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

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

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

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

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

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

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

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

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

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

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

3 時(shí)間譜的計(jì)算

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

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

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

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

4 結(jié)論

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