中子及次级γ在高空长距离蒙特卡罗输运模拟中的减方差方法

刘 利,左应红,牛胜利,朱金辉,商 鹏,李夏至

(西北核技术研究所,西安 710024)

高空核爆炸核辐射环境研究[1-3]对低轨道卫星、临近空间(20～100 km)飞行器及通信传播[4]等都有着重要意义。模拟中子及次级γ在高空大气中的大空间长距离输运时，辐射源与记录点的距离可达几十公里，甚至上百公里，使到达记录点附近的粒子概率极小。以50 km为例，真空中各向同性点源向外出射的粒子落在50 km外1 m2区域内的概率约为5×10-10。同时中子由高空向低层大气输运时，大气密度迅速增加，中子平均自由程迅速减小，中子与大气的多次散射作用使粒子到达记录点的概率进一步减小。高度为20 km时，1 MeV中子和γ在水平方向的自由程约为1.1 km和1.4 km，远小于辐射源与记录点的距离。当记录点附近的粒子数太少时，蒙特卡罗方法模拟结果的相对偏差较大，结果不可信。对瞬发中子、γ和X射线在高空中输运模拟已有研究[4-7]，但针对的是空气质量厚度较小及大气散射作用较弱的情况，针对中子及次级γ在高空长距离输运至低层(如20 km处)这种空气质量厚度较大的大气环境的研究还未见报道。中子及次级γ在空气质量厚度较大的大气环境中的输运模拟计算对研究核电磁脉冲环境和通信传播等有重要意义。本文开展了中子在高空大气中的长距离输运模拟的减方差方法研究，建立了中子在高空非均匀大气输运的几何模型，采用MCNP程序，比较了源偏倚、几何分裂、截断、权窗、强迫碰撞和DXTRAN球等减方差方法[8-11]计算结果的相对偏差与计算效率，并提出源偏倚、几何分裂和DXTRAN球组合方法为最佳减方差方法。

1 中子输运模拟几何模型

图1给出了中子在高空大气中长距离输运模拟的几何模型。高空大气是典型的非均匀连续介质，密度随高度变化非常剧烈。模型将大气结构进行分层处理，用于描述大气密度随高度的变化。高度越高大气密度越小，中子平均自由程越大，分层处理时大气厚度也越大。根据大气密度的特点与计算精度需要，将高度为0～50 km的大气分为50层，每层厚度为1 km；将高度为50～110 km的大气分为6层，每层厚度为10 km。每层大气的密度根据美国标准大气模式取值[12]。地平面以下5 km为土壤层，密度为2.31 g·cm-3。

设置中子辐射源为位于高度为80 km爆心处的各向同性点源，中子权重设置为1，中子能谱取典型氢弹出壳中子能谱[3]。记录点位于爆心正下方。由于中子飞行到20～40 km高度时才会与大气发生明显的相互作用，记录点的高度取为20～40 km。爆心与最低记录点的距离为60 km，且中子在大气中输运时受到多次散射作用，使直接模拟的粒子能到达最低记录点的概率非常小。

设置模拟源粒子数目N为1×107。在20～40 km高度范围内每2 km设置一个记录点，采用点计数方法统计记录点的中子及次级γ注量。图2为直接模拟给出的中子与γ的注量、相对偏差σ和FOM因子随高度的变化关系。FOM因子可表示为

(1)

其中，t为计算时间，min。

FOM因子可反映出模拟结果的计算效率，FOM因子越大，减方差效果越好。由图2(a)可见，中子注量随记录点高度降低而降低，γ注量随记录点高度变化幅度相对较小。由图2(b)可见，中子注量和γ注量的相对偏差随记录点高度降低而增加，FOM因子随高度降低而降低。当记录点高度小于30 km时，中子注量和γ注量的相对标准偏差均大于5%，结果不可信。即使模拟粒子数逐渐增加至1×108，直接模拟得到的中子和γ注量依然没有收敛。图3为直接模拟和采用减方差(variance reduction，VR)方法计算给出的20 km高度处中子与γ的注量、相对偏差和FOM因子随粒子数的变化关系。由图3可见，直接模拟的中子与γ注量、相对偏差和FOM因子随粒子数的增加而大幅度波动；采用减方差方法(图3中为源偏倚+几何分裂+DXTRAN球)可有效降低统计涨落，减小相对偏差，提高FOM因子。

2 减方差方法

2.1 单一减方差方法

常用的蒙特卡罗模拟减方差方法可分为3类：(1)改进抽样方法，包括源偏倚、强迫碰撞、指数变换、重要性抽样和统计估计抽样等；(2)总体分布控制法，包括权窗、分裂与轮盘赌技巧、截断、分层处理和两步法等；(3)半解析法，包括点探测器、环探测器和DXTRAN球等。本文以20 km高度处中子和γ注量计算优化为目标，开展中子及次级γ在高空中长距离输运的减方差方法研究。本文采用的减方差方法主要有：

2.1.1 源偏倚

对粒子源的方向和能量等参数做偏倚，增加重要区域的抽样粒子数，并相应减小粒子的权重。在图1所示的输运几何模型中，爆心高度处中子平均自由程很大，由爆心向上发射的粒子对爆心正下方的记录点几乎不产生计数贡献。为提高抽样效率，采用源方向偏倚提高粒子到达记录点的概率。根据模型几何结构，设置源偏倚使源出射方向与垂直地面向下方向的夹角余弦在0.948 7～1范围内和范围外的概率各为50%，相应修改粒子的权重以保证模拟的无偏性。夹角余弦为0.948 7时，对应的夹角为18.4°，可覆盖几何模型中以20 km高度处记录点为中心，水平方向距离为0～20 km的范围。

2.1.2 几何分裂与轮盘赌

将模型划分为多个栅元，设置重要栅元的重要性较高，不重要栅元的重要性较低。当粒子进入重要性高的栅元内，有一定的概率分裂为多个粒子；当粒子进入重要性低的栅元，则只有一定的概率存活，从而增加重要性较高的栅元内的模拟粒子数。对图1所示几何模型，设置高度为40 km以上的栅元重要性为1，40 km以下的栅元重要性每2层翻一倍，直到栅元重要性达到28；高度为20 km以下的栅元重要性每2层降低至一半，直到栅元重要性为1。

2.1.3 强迫碰撞

在指定的栅元内人为增加粒子碰撞抽样概率，并将进入栅元的粒子分为碰撞与不碰撞2部分。不碰撞部分直接离开栅元，粒子权重相应减少。碰撞部分具有粒子剩余权重，并使用轮盘赌来控制粒子数量，防止出现大量低权重粒子而增加计算时间。针对图1所示几何模型，设置粒子在高度为18 km以上的栅元内发生强迫碰撞，18 km以下时，不发生强迫碰撞，从而引导粒子进入20 km高度处的记录点附近。

2.1.4 指数变换

在特定方向上减小微观截面并在相反方向上增加微观截面，从变换后的非模拟密度函数抽样碰撞距离，相应修改粒子权重，从而增加指向特定方向的模拟粒子数。修改后的截面为

(2)

2.1.5 DXTRAN球

使用DXTRAN球时，设置一包围感兴趣区域的小球。在球外的每次碰撞都能产生特殊的DXTRAN粒子，按照确定的概率进入DXTRAN球内，权重利用确定论方法计算。碰撞产生的粒子正常输运，只是在进入DXTRAN球范围时被终止模拟，保证结果的无偏性。对图1所示几何模型，设置DXTRAN球的球心位于20 km高度处，半径取1.5 km。

2.1.6 权窗

权窗是与位置和能量相关的分裂和轮盘赌技巧，是最常用的减方差方法之一。设置圆柱形栅格权窗，在高度方向上按图1所示模型划分，半径方向上网格宽度逐渐增加。对中子与γ联合输运，同时使用中子与γ的权窗，并取20 km处的γ注量来产生权窗。设置模拟源粒子数目为1×108个，表1为单独采用某一减方差方法的计算效率。表1第3行还给出了采用模拟粒子数为1×1010时的直接模拟结果，用于验证其他减方差方法的无偏性。当模拟粒子数增加至1×1010时，计算相对偏差有所降低，但表征计算效率的FOM因子也大幅度降低。

表1 单独采用某一减方差方法的计算效率Tab.1 Computational efficiency by using a single VR method

由表1可知，与指数源偏倚方法相比，分布式源偏倚方法将FOM因子分别从4.9和0.2提高至23和13，计算偏差分别从6.4%和32%降低至3.0%和4.0%，所以本文采用分布式源偏倚方法；与直接模拟相比，采用源偏倚、几何分裂和DXTRAN球方法计算的相对偏差都有所降低，同时FOM因子有所提高，这3种方法都引导粒子向记录点位置偏倚，大大提高了计数区域内的粒子抽样概率，从而大幅度提高了计算效率，同时降低了计算偏差。

由表1还可知，与直接模拟相比，DXTRAN球方法的FOM因子提高最多，由4.6和3.2提高至52和39；几何分裂方法的相对标准偏差降低最多，由9.1%和10.8%降低至0.9%和1.3%；强迫碰撞方法的相对偏差有所降低，但FOM因子同样降低；指数变换方法的相对偏差较大，且FOM因子偏低，减方差效果较差，是因为收敛性太差，修改拉伸系数没有使减方差效果明显好转，与其他减方差方法结合使用可能会提高减方差效果；采用权窗进行一次迭代的计算结果将中子注量计算的FOM因子提高了2个量级，相对偏差降低至1%，但计算结果明显偏离正常计算值，与采用模拟粒子数为1×1010时直接模拟的结果相差较大，且γ注量计算相对偏差较大；采用权窗进行二次迭代的中子与γ的计算相对偏差都较大，且FOM因子较低，这说明一次迭代计算中子注量的相对偏差为伪偏差，计算结果已完全不可信。

2.2 减方差方法的组合优化

组合使用多种减方差方法可综合利用各种减方差方法的优点，通常比单独采用一种减方差方法具有更高的计算效率[10]。对于图1所示几何模型，各向同性发射的源粒子输运至记录点位置的概率极低，采用方向源偏倚方法使源粒子更多地朝指定方向发射。表2为采用源偏倚方法和其他一种减方差方法组合的计算效率。由表2可知：源偏倚与几何分裂或指数变换组合的方法大幅度提高了计算效率；源偏倚与DXTRAN球组合的方法比单独采用DXTRAN球方法的减方差效果略有提高；源偏倚与DXTRAN球方法通过调整粒子方向或产生DXTRAN粒子都使粒子朝记录点运动，与源偏倚与几何分裂等组合方法相比，源偏倚与DXTRAN球组合的方法计算效果较差一点；源偏倚与强迫碰撞组合的方法有效减小了相对偏差，但也降低了FOM因子。

表2 源偏倚与其他一种减方差方法组合的计算效率Tab.2 Computational efficiency by using source bias method and another VR methods

由表2还可知，源偏倚与权窗组合方法给出的相对偏差和FOM因子随着迭代次数的增加并未收敛，反而出现强烈的波动，即使某一步迭代的中子注量相对偏差较小，中子注量计算值也与其他方法相差较大，同时对应的γ注量相对偏差较大，相应的权窗下限也出现了明显的波动，且有时会出现较多的权窗下限为“0”的权窗；源偏倚+几何分裂+权窗、源偏倚+DXTRAN球+权窗和SGD(source bias+geometry splitting+DXTRAN)+权窗等组合方法的计算结果同样无法收敛且相对偏差较大。在高空中子及次级γ长距离输运模拟中，权窗下限中出现大范围的“0”权窗下限，使模拟粒子过度朝向非“0”权窗下限栅元偏倚，导致权窗无法准确描述粒子对探测器计数的贡献。权窗不适用的原因可能有2个：一是空间尺度太大，较远栅元内的模拟粒子对探测器计数的贡献被忽略，从而使权窗下限为“0”，额外计算表明，在采用权窗方法模拟高度为10 km范围内的高空中子输运可有效减小相对偏差；二是高空大气密度变化剧烈，高度为20～80 km时，空气密度变化3～4个量级，高密度区域内权窗网格尺寸应该比低密度权窗网格尺寸要小得多，但MCNP自带的网格权窗参数难以同时满足高密度与低密度区域权窗网格尺寸的需求。此类“0”权窗下限问题在不少深穿透问题的研究中都有出现[13]，有关文献还提出采用源偏倚、指数变换、密度迭代和猜测初始权窗等方法来获得合理的权窗参数，但对于百公里量级非均匀大气中中子及次级γ长距离输运问题，采用以上方法和本文所述方法均未能得到合理可靠的结果，还需采用几何迭代、自适应权窗和全局权窗等新的权窗生成方法进一步研究。

表3为采用源偏倚与其他多种减方差方法组合的计算效率。结合表1、表2和表3可知，采用强迫碰撞或在其他方法上加入强迫碰撞可有效降低计算结果的相对偏差，但同时FOM因子也明显降低。这是因为对于图1所示几何模型，为引导粒子输运至记录点，需在相邻栅元内设置强迫碰撞，导致产生大量低权重粒子，大大增加了计算时间；在源偏倚+几何分裂、源偏倚+DXTRAN球和SGD组合方法上再加上指数变换的减方差效果并不理想。结合表1、表2和表3还可知，SGD组合方法的减方差效果最好，中子与γ的FOM因子分别由直接模拟法的4.6和3.2提高至165和41，分别提高了39倍和13倍，相对偏差由9.1%和10.8%降至0.2%和0.5%。SGD组合方法与模拟粒子数为1×1010时直接模拟的结果较为接近，相对差异在模拟结果的相对偏差范围之内，验证了该组合方法的无偏性。图3中的减方差方法就是SGD组合方法。由图3可见，SGD组合方法的模拟结果收敛性较好，随着粒子数的增加，注量的计算值基本保持不变，且相对偏差和FOM因子变化幅度较小。

表3 源偏倚与其他多种减方差方法组合的计算效率Tab.3 Computational efficiency by using source bias method and other multiple VR methods

采用SGD组合方法计算给出的中子与γ的注量、相对偏差和FOM因子随高度的变化关系，如图5所示。由图5(b)可见，高度为20 km处记录点的相对偏差比22 km处还要小，但FOM因子比22 km处数据大，是因为在20 km处设置了DXTRAN球，提高了计算效率与计算精度。采用该组合方法模拟得到由高空输运至20～40 km高度范围内的中子和γ注量的相对偏差小于1%。若将模拟粒子数由1×108降低至1×107，计算结果的相对偏差能保持在5%以内。由此可见，采用SGD组合方法也适用于模型中其他高度记录点的中子和γ注量计算。

3 时间谱的计算

在高空核辐射环境与效应研究中，中子或γ的时间谱是极其重要的参数[3]。采用不同减方差方法模拟计算得到的高度为20 km处，中子与次级γ的注量率随时间的变化关系，如图6所示。由图6可见，SGD组合方法计算得到的中子及次级γ的时间谱较平滑，波动较小；直接模拟得到的时间谱波动较大，且收敛性较差。这说明组合方法同样适用于计算中子及次级γ时间谱。

由图6还可见，SGD组合方法计算得到的中子次级γ时间谱的相对偏差在γ峰值处较小，在γ峰值后逐渐增加。不同的时间段内γ时间谱的相对偏差不一致。针对这种情况，本文采用了时间分裂的减方差方法。与几何分裂和能量分裂类似，时间分裂通过分裂与轮盘赌来增加关心时间段内的模拟粒子数，降低计算结果的相对偏差。设置模拟时间达到2×10-4s时，粒子被一分为二；时间达到8×10-4s时，粒子以0.5的概率存活；时间达到2×10-3s和5×10-3s时，粒子再次被一分为二。

由图6还可见，采用SGD+时间分裂组合方法计算的γ时间谱首末两端更加平滑，相对偏差也有所降低。表4列出了不同减方差方法计算时间谱得到的相对偏差σ。其中σ(>2×10-3s)指的是时间点大于2×10-3s的相对偏差。由表4可知，采用减方差方法后，中子及次级γ时间谱的相对偏差明显减小，在此基础上加入时间分裂算法可进一步减小相对偏差，但计算时间略有增加。将模拟粒子数设置为1×109时，SGD+时间分裂组合方法计算中子及次级γ时间谱的相对偏差小于5%。

表4 不同减方差方法计算结果的相对偏差Tab.4 The relative deviation of time spectra of neutron and secondary γ by using different VR methods

4 结论

针对中子及次级γ在高空长距离蒙特卡罗输运模拟的效率不高问题，研究了源偏倚、几何分裂、强迫碰撞、权窗、指数变换、DXTRAN球和时间分裂等减方差方法在高空长距离输运中的应用，通过比较各个减方差方法的相对偏差与计算效率，得出SGD组合方法的减方差效果最好，与直接模拟相比，SGD组合方法给出的中子及次级γ的FOM因子分别提高了39倍和13倍，得到高度为20～40 km时的中子及次级γ注量的相对标准偏差小于1%。在该方法基础上加上合理设置参数的时间分裂算法可进一步降低时间谱的相对偏差。