全局减方差方法在大空间γ辐射场计算中的应用

刘利 左应红 牛胜利 朱金辉 商鹏 王学栋

（西北核技术研究所 西安 710024）

γ射线在近地面大空间长距离输运中会与空气和地面土壤发生多次散射，使得输运至距辐射源公里及以上范围内的γ粒子数较少，致使该处的蒙特卡罗模拟统计涨落较大。左应红等［1］建立了栅元权窗结合密度逼近迭代的局域减方差方法（Local Variance Reduction，LVR）［2-3］，能够求解离源5 km的观测点处的γ射线注量。而在核设施厂房辐射环境［4］、早期核辐射环境［5］、核电磁脉冲环境［5］、人员辐射剂量学等方面的研究中，常常需要计算大空间中γ辐射场的整体分布。传统的LVR方法通常只针对特定目标探测器，难以给出理想的全空间γ辐射场分布。

相比之下，全局减方差方法（Global Variance Reduction，GVR）［6-9］更适于求解辐射场的整体分布。GVR方法通过全局权窗来控制模拟粒子的空间分布，使得模拟结果的相对偏差随空间的分布尽可能减小，从而整体提高全局空间内的计算效率。Davis和Andrew等［6-7］提出的基于蒙特卡罗正向计算的GVR方法，通过蒙特卡罗模拟计算得到粒子通量、数量和权重等粒子分布信息，用于产生基于栅元/网格的全局权窗参数，利用全局权窗参数开展蒙特卡罗模拟，循环迭代直至结果收敛。该方法思路清晰，适用性强，已广泛用于反应堆装置［10-12］、托卡马克装置［13-14］和厚屏蔽体［15-16］等复杂模型的全局辐射环境计算。

大空间γ辐射场模拟具有空间尺度跨度大和计数区外散射较强等特点。直接应用已有的GVR方法面临计数栅元/网格体积差异和非计数区引起的“过度分裂”问题。一个是辐射场空间尺度跨度大，计数栅元/网格体积相差较大。而现有全局减方差方法更适用于等体积的栅元/网格，其引导粒子输运使得模拟粒子在全空间范围内分布较为均匀，相对误差在全空间范围内相对均匀。当辐射场内计数栅元/网格体积差异较大时，体积大的栅元/网格内模拟粒子数较多，体积小的栅元/网格内模拟粒子数较少，使得模拟相对误差在栅元/网格内分布不均，从而降低全局计算效率。另一个是计算边界附近的非计数区占用的计算资源较多。由于计数区外土壤和大气对γ射线的散射作用较强，大空间γ辐射场求解时在计数区外还需要非计数区。直接将权窗参数统一应用于计数区与非计数区，非计数区内（尤其截面较大的吸收介质如土壤）模拟粒子过度分裂，占用较多计算资源，从而降低计算效率。

因此，本文在原有GVR方法上引入了体积修正和非计数区修正，解决栅元/网格体积差异和非计数区造成的粒子过度分裂问题，提高GVR方法在大空间γ辐射场的适用性。同时开展基于粒子误差、权重、径迹、数目、能量、碰撞和通量的共7种GVR方法在大空间γ辐射场的应用研究，分析比对各个GVR的计算效率与特征参数，得到减方差效果最好的GVR方法。最后在基于通量的GVR方法上引入光滑因子SI，研究了SI对于GVR方法计算的全局权窗参数和计算效率的影响。

1 大空间γ辐射场计算几何模型

图1为本文建立的大空间γ辐射场全局计算几何模型。模型主要有空气和土壤两种介质组成，其中空气密度取标准大气密度1.205×10-3g·cm-3，土壤密度取2.35 g·cm-3。设置γ注量的计数区域为高度0～1 km、距源心水平距离5 km范围内（图1中蓝色区域）。由于计数区外的土壤和空气对于γ射线有较强的散射作用，计数区外应建立非计数区。故模型最大高度取2 km，距源心水平最大距离取6 km，土壤厚度取1 m。计算表明，继续增加边界尺寸对计数区内γ注量影响可以忽略。利用MCNP（Monte Carlo N Particle Transport Code，MCNP）［2］开展计算，设置γ辐射源各向同性发射，权重为1，能谱为watt平衡裂变谱，源心距地面高度设置为20 m。为防止过度分裂造成的机时太长，文中所有模拟设置CPU总运行时间截断为400 min。

图1 大空间γ辐射场计算几何模型（彩图见网络版）Fig.1 Geometric model for the computation of large-space γ radiation field (color online)

由于GVR方法需要用到栅元/网格内粒子通量、数量、权重等粒子分布信息作为输入参数，本文模型将空气与大地进行分层处理。对于空气，在水平方向上，每100 m分一层，共分为60层，计算所关心的0～5 km范围内共分50层；在垂直方向上，每100 m分一层，共分为20层，计算所关心的0～1 km高度范围内共分为10层。对于大地，在垂直方向上，每0.125 m分一层，共分8层；在水平方向上与大气保持一致，每100 m分一层，共分60层。模型整体被划分为28×60=1680个栅元，蓝色计数区域共被划分为10×50=500个栅元。

能量为1 MeV的γ射线在标准大气中的自由程约100 m。本文γ射线平均能量略大于1 MeV，保守按照平均自由程125 m来估算，γ射线在模型中最大输运距离超过40个平均自由程，不采用减方差方法的情况下实际能输运至远处的γ射线数目极少。由于γ射线各向同性向外发射，同时射线与空气和大地会发生多次散射作用，距离源心越远的栅元内粒子数越少，模拟统计涨落越大，使得整体计算误差呈中心小边缘大的特点。中心栅元内γ注量比边缘栅元内γ注量高10个数量级以上，粒子分布呈中心多边缘少的特点。直接增加模拟的粒子数使得大量计算时间用于模拟误差小的中心处的粒子输运，而误差较大的边缘网格处模拟粒子数不能有效增加，计算效率低。因此有必要研究GVR方法，通过设置全局权窗来增加边缘栅元处的模拟粒子数，减小中心网格处的模拟粒子数，从而在保证整体无偏的情况下降低全局计算的相对偏差。

2 GVR方法与评估参数

在蒙特卡罗模拟中，问题通常被分为三类：1）源-点探测器响应；2）源-区域通量或者响应；3）整个问题空间的全局通量或响应。前两类问题通常采用局域减方差方法，而后一种应采用GVR方法。GVR方法中最重要的减方差技巧就是全局权窗。通过设置权窗来控制粒子的整体分布使得整个问题空间内的通量平均计算效率提高。通常利用确定论离散纵标法（discrete ordinates Shifted N2Box，SN）或者蒙特卡罗（Monte Carlo，MC）方法进行正向输运计算或者伴随输运计算，将计算的粒子空间分布信息或粒子对目标探测器的响应贡献转化为权窗参数，进而采用权窗进行输运计算。根据不同的计算方法与计算内容的两两组合，GVR方法通常可以分为4类：1）SN方法正向输运［17］；2）MC方法正向输运；3）SN方法伴随输运［18］；4）MC方法伴随输运［19］。其中正向输运计算可以得到粒子通量的空间分布，从而设置输运偏倚参数，适用于源项小、屏蔽体厚的计算模型。大空间γ辐射场计算是典型的第三类全局通量求解问题，具有点源、深穿透问题等特征，适合基于正向输运的计算方法。

本文采用的基于MC方法正向输运的减方差方法通过蒙特卡罗模拟正向计算得到粒子通量、权重和误差等粒子分布信息，进而设置全局权窗参数来指导粒子输运。设置全局权窗的原则是在高重要性区域内设置较小的权窗下限使得更多的粒子分裂，降低统计方差；在低重要性区域内设置较大的权窗下限（Lower limit of the weight window）wth，对不太重要的粒子进行轮盘赌，提高计算效率。当前栅元/网格的重要性为［20-21］：

式中：P为所有曾进入该栅元/网格的粒子对于计算结果的贡献之和；wall为所有曾进入该栅元/网格的粒子权重之和。

权窗下限与重要性成反比：

对于GVR方法来说，每一个栅元/网格对整体计算空间贡献相等，即所有曾进入每一个栅元/网格的粒子对于计算结果的贡献之和相等。则有权窗下限：

由于栅元/网格内粒子通量通常与粒子权重成正比，wth∝ϕi，ϕi为第i栅元/网格内粒子平均通量。粒子通量是蒙特卡罗模拟中的常见物理量，通常采用粒子通量来计算全局权窗参数。相应基于粒子通量的GVR方法（Flux-based GVR，GVR-F）的权窗下限为：

式中：Max(ϕi)为所有栅元/网格中最大粒子平均通量。每一个栅元/网格权窗上限下限之比β一般固定，比如本文所有计算取5。

另一种常用于GVR方法的粒子信息就是相对误差。误差大说明该处粒子数少统计涨落较大，应降低权窗下限增加粒子分裂。则有基于误差的GVR方法（Error-based GVR，GVR-R）的权窗下限为：

其中：Rei为第i栅元/网格内粒子通量的误差。在直接模拟中相对误差与粒子通量平方根成反比，。所以基于误差的权窗下限又可以表示为：

利用计算得到的全局权窗参数开展新的蒙特卡罗模拟，得到新的粒子分布信息。重复迭代直到结果收敛，从而得到全局精确结果。

为评价GVR方法的优劣引入全局品质因子FOMG和基于误差的标准差σ。全局品质因子为：

式中：N为栅元/网格的总个数；T为模拟的CPU计算时间。

基于相对误差的标准差为：

一般来说，FOMG值越大计算效率越高，σ值越小平均统计误差越小。GVR方法的加速比定义为GVR方法的全局品质因子FOMG与直接模拟的FOMG之比。

3 GVR方法在大空间γ辐射场的应用

3.1 体积修正

图1中大空间γ辐射场计算几何模型在径向和垂直方向分别被划分为60层和28层，生成相应的非等体积的圆环形栅元。中心处栅元体积最小，径向边缘处栅元体积较大。采用圆环形的栅元相比于长方体栅元可以有效增加远离中心栅元的体积，从而增加栅元内模拟粒子数，降低计算误差提高计算效率。图2给出了采用圆环形栅元和长方形栅元计算的γ注量和误差随距离的变化。其中长方形栅元边长为100 m，与圆环形栅元径向间隔一致。采用等体积的长方形栅元统计计算的栅元内γ注量平均误差远大于圆环形栅元。

图2 采用圆环形栅元和长方形栅元计算的γ注量和误差随距离的变化Fig.2 Calculated γ fluence and error in the circular and rectangular grid cells with respect to distance

式（4）和（6）的全局权窗下限适用于等体积的栅元和权窗，其引导粒子输运使得蒙卡模拟粒子在全空间范围内分布较为均匀，相对误差在全空间范围内相对均匀。对于体积非均匀的栅元/网格，采用式（4）和（6）的全局权窗下限，体积大的栅元/网格内模拟粒子数较多，体积小的栅元/网格内模拟粒子数较少，使得模拟相对误差在全局空间范围内分布不均。为使每一个非等体积的栅元/网格内的模拟粒子数趋近相等，本文引入体积修正因子fvol，对权窗参数进行体积修正。基于通量的权窗下限式（4）修正后为：

式中：Vi为第i个栅元/网格的体积；Vs为点源所在的栅元/网格的体积。

表1给出了直接模拟，采用LVR方法、GVR-F方法和引入体积修正的GVR-F方法迭代一次模拟得到的主要参数。其中采用LVR方法时参考文献［1］设置水平5 km、高度500 m处的环探测器计数作为权窗生成的目标探测器。采用LVR方法并合理设置权窗生成器可以有效降低整体误差水平，并提高FOMG因子，降低标准差σ。由于LVR方法只针对局域进行优化，水平距离5 km、高度1 km附近计算误差较大，且水平距离小于5 km、高度1 km处某些点误差比直接模拟还大（图3（b））。即使迭代三次后，LVR方法模拟得到的计算区域内最大误差依然有38%。

表1 采用不同减方差方法模拟所得的主要参数Table 1 Main parameters simulated using different variance reduction methods

图3 采用不同减方差方法模拟所得γ注量误差分布Fig.3 Error distribution of γ flux obtained using different variance reduction methods

采用不考虑体积修正的GVR-F方法计算得到的FOMG因子比直接模拟结果还要小。其原因为该方法引导模拟粒子在全局空间均匀分布，但同一高度处径向最外层栅元体积为最内层栅元体积的119倍，若不考虑空气的衰减，最外层栅元内模拟粒子数远大于内层栅元模拟粒子数，使得模拟粒子在外层栅元内过度分裂。图3（c）中GVR-F方法模拟得到径向4～5 km范围内的γ注量误差明显小于直接模拟。但由于“过度分裂”问题，在相同的CPU运行时间内模拟的总粒子数比直接模拟小3～4个数量级，使得径向4 km以内γ注量误差比直接模拟还大。增加模拟粒子数或者提高CPU运行时间也不能进一步提高该方法的FOMG因子。引入体积修正后，GVR方法引导粒子在全栅元/网格内均匀分布，防止出现“过度分裂”现象。体积修正的GVR-F方法计算的FOMG因子比直接模拟提高39倍，基于误差的标准差σ降低约两个数量级。

3.2 非计数区修正

大空间γ辐射场计算模型在计数区外还存在非计数区域。非计数区域内粒子可以输运至计数区，影响计数区粒子通量计算。大尺度空间辐射场计算必须要合理设置计数区与非计数区。一般来说，非计数区的厚度应达到粒子的饱和反射厚度。将GVR方法计算的权窗下限统一应用于计数区与非计数区，虽然计数区计算精度得以保证，但非计数区内（尤其本文算例中的土壤等吸收介质内）模拟粒子数较多，占用较多计算资源，从而降低计算效率。

根据§2，计数区内栅元/网格的权窗下限与重要性成反比，与粒子平均通量成正比。非计数区内模拟粒子只有输运至计数区才能对计数区粒子通量有贡献。非计数区栅元/网格内粒子输运至计数区栅元/网格的概率为e-s/λ，其中，s为非计数区栅元/网格到计数区栅元/网格的距离；λ为粒子的平均吸收自由程。以离非计数区最近的计数区栅元/网格为基准，可以设置非计数区的权限下限为计数区权窗下限与粒子从非计数区输运至计数区的概率：

式中：wth，none为非计数区栅元/网格的权窗下限；wth，n为离非计数区最近的计数区栅元/网格的权窗下限。本文中空气中γ的平均吸收自由程λ取250 m。该修正改变了权窗下限的数值，实际输运模拟时根据权窗来分裂或者杀死粒子，同时改变粒子权重来保证结果的无偏性。

图4 给出了采用考虑体积修正的不同GVR方法计算的最靠近地面一层栅元内的权窗下限wth分布。权窗下限wth随水平距离增加而迅速下降，不考虑非计数区修正，wth一直下降至6 km的计算边界。以式（11）引入非计数区修正后，wth在计数5 km处达到最小值，而后在非计数区随水平距离而增加。引入非计数区修正后的wth变化规律与文献［1］中以5 km处注量为目标的LVR方法计算的wth变化规律类似。

图4 采用体积修正的不同GVR方法计算的近地面栅元内权窗下限wth（彩图见网络版）Fig.4 Lower limit of the weight window wth in the nearground cell computed using different GVR methods with volume modification (color online)

表2给出了采用体积修正和不同非计数区修正的GVR方法迭代一次所得的主要参数。直接设置非计数区内权窗下限wth，none=0关闭权窗，模拟得到的FOMG因子不增反降，基于误差的标准差σ有所增加，最大误差大幅度增加。采用非计数区修正的GVR-F方法相比于无计数区的GVR-F方法计算的大空间γ辐射场的FOMG因子提高约39%，基于误差的标准差σ降低约27%；采用非计数区修正的GVRR方法计算的FOMG因子相应提高约18%，基于误差的标准差σ相应降低约18%。由于体积修正与非计数区修正使得本文的大空间γ辐射场计算效率显著提高，下文对GVR方法的对比研究中默认都采用了体积修正与非计数区修正。

表2 采用体积修正和不同非计数区修正的GVR方法计算的主要参数Table 2 Main parameters computed using volume modification and different non-count area modifications

4 不同GVR方法的对比

4.1 粒子信息

除了常用的基于通量的GVR-F方法和基于误差的GVR-R方法，聂星辰、李佳和Andrew Davis等［6，22-23］还采用了基于栅元/网格内粒子能量（Energy-based GVR，GVR-E）、数目（Populationbased GVR，GVR-P）、径迹（Track-based GVR，GVRT）和权重（Weight-based GVR，GVR-W）的GVR方法，并在国际热核聚变实验堆（International Thermonuclear Experimental Reactor，ITER）和中国聚变工程实验堆（China Fusion Engineering Test Reactor，CFETR）托卡马克装置中子环境计算中取得良好的减方差效果。以上方法认为全局权窗下限与栅元/网格内粒子能量、数目、径迹、权重和碰撞次数成正比，其权窗下限计算方法与基于粒子通量的权窗下限基本一致，只需将式（4）中ϕi分别替换成第i栅元/网格内粒子能量Ei、数目Ni、径迹Ti和权重wi等。基于类似原理，本文还提出了基于栅元/网格内粒子碰撞的GVR（Collision-based GVR，GVR-C）方法，同样将式（4）中ϕi替换成第i栅元/网格内粒子碰撞数Ci即可。

表3给出了采用不同的GVR方法模拟迭代3次后得到的主要参数。通过不断迭代，GVR方法增加统计涨落较大区域内的模拟粒子数，使得模拟粒子数在整个计数区空间中偏向于均匀分布，尤其是迭代2～3次后GVR-P方法模拟的最靠近地面栅元内模拟粒子数基本相等（图5）。迭代2～3次后，不同GVR方法给出的主要参数已经收敛。GVR-W方法模拟得到的FOMG因子提高约2300倍，基于误差的标准差σ降低约5750倍，模拟得到计数区栅元内γ注量平均误差为0.6%，最大误差为1.8%，在所有GVR方法中减方差效果最好。GVR-T、GVR-P、GVR-E、GVR-C和GVR-F方法模拟得到的FOMG因子提高至2～3，加速比达到500～700倍，模拟得到的基于误差的标准差σ降至10-5量级。GVR-R方法在所有GVR方法中减方差效果最差，但FOMG因子依然提高至0.87，加速比能达到190倍。对比图6中迭代3次的GVR-R、GVR-W和GVR-F三种方法计算γ注量误差分布，GVR-W模拟得到的计数区栅元内γ注量误差整体比GVR-R和GVR-F的误差小。三种方法计算的栅元误差最大值都位于水平距离5 km靠近地面处。

表3 采用7种GVR方法模拟所得的主要参数（迭代3次）Table 3 Main parameters simulated using seven GVR methods (3rd iteration)

图5 GVR-P方法模拟的栅元内模拟粒子数（彩图见网络版）Fig.5 Simulated particle number Nps in the near-ground cells simulated using GVR-P method (color online)

图6 三种GVR方法计算的γ注量误差分布Fig.6 Error distribution of γ flux computed using three GVR methods

图7和图8给出了GVR-W方法计算的γ注量分布和权窗下限wth分布。大空间γ辐射场计数区栅元内γ注量呈中心高边缘低的特征，计数区γ注量总衰减超过10个数量级。对于注量变化如此大的深穿透厚屏蔽体问题，GVR-W方法模拟的所有栅元内γ注量误差最大仅为1.8%，验证了合理的GVR方法在大空间辐射场计算中的实用性。GVR-W方法计算的计数区内权窗下限wth同样呈现中心高边缘低的特征，计数区外wth逐渐增加，wth最低衰减至10-7量级。

图7 GVR-W方法计算的γ注量分布Fig.7 γ flux distribution computed using three GVR-W methods

图8 GVR-W方法计算的权窗下限wth分布Fig.8 Distribution of wth computed using GVR-W methods

不同GVR方法模拟的FOMG因子及相应加速比不同的主要原因在于不同GVR方法计算的权窗下限wth不同。如图9所示，GVR-R方法计算的wth最低衰减至10-4量级，GVR-W方法计算的wth最低衰减至10-7量级，GVR-T、GVR-P、GVR-E、GVR-C和GVR-F方法计算的wth最低衰减至10-9量级。由于计算的wth分属于三个水平，不同GVR方法计算的FOMG因子也分属于三个水平，其中GVR-W方法计算的wth衰减至10-7量级，对应FOMG因子提高最多，全局减方差效果最好。

图9 基于不同粒子信息的GVR方法计算的近地面栅元内权窗下限wth（彩图见网络版）Fig.9 wth in the near-ground cell computed using GVR methods based on information of different particles(color online)

4.2 光滑因子

§4.1中7种GVR方法减方差效果差异的原因为计算的wth差异明显。对于本文中的大空间γ辐射场计算模型，GVR-W方法全局减方差效果最好。但如果模型变化，比如辐射源由γ换成中子、计数区增大等，GVR-W方法计算的wth并不一定最优。

为进一步研究wth对全局减方差效果的影响，在式（4）中引入光滑因子SI，相应权窗下限为：

采用光滑因子SI可以方便快捷地修改wth。当SI=1对应基于通量的GVR-F方法，SI=0.5对应基于误差的GVR-R方法。

在本文计算模型中，光滑因子SI由0.4增加至1，对应的权窗下限wth最小值由10-3量级降低至10-9量级（图10）。不同光滑因子的GVR方法计算的FOMG因子和加速比随SI的增加先增加后减小，计数区栅元内γ注量的平均误差Rave、最大误差Rmax和标准差σ随SI的增加先减小后增加（表4）。

表4 不同光滑因子的GVR方法计算的主要参数（迭代3次）Table 4 Main parameters computed using GVR methods with different SI values (3rd iteration)

图10 不同光滑因子的GVR方法计算的近地面栅元内权窗下限wth（彩图见网络版）Fig.10 wth in the near-ground cell computed using GVR methods with different SI values (color online)

当SI过小时，wth较大，模拟粒子在计数区内分裂数目较少，统计涨落较大，模拟误差较大。当SI过大时，wth较小，模拟粒子在计数区内分裂数目较多，分裂过度，计算效率下降。因此，存在一个合适的光滑因子SI，使得模拟计算效率最高。对于本文大空间γ辐射场，使得计算效率最高的光滑因子SI为0.8。当SI=0.8时，GVR方法计算的FOMG因子最大为14.78，平均误差最小为0.5%，标准差σ最小为1.4×10-5，加速比最大为3246，减方差效果最好。SI=0.8的GVR方法计算FOMG因子比GVR-W方法还提高40%。合理设置光滑因子SI，采用式（12）计算权窗下限wth，使用权窗开展蒙特卡罗模拟，迭代至收敛，就可以实现全局辐射环境的精确计算。

5 结语

为解决大空间γ辐射场的全局精确计算问题，本文开展了多种GVR方法在大空间γ辐射场的应用研究。针对计数栅元/网格体积差异较大的情况，提出了体积修正方法，适用于非等体积栅元/网格的GVR方法的权窗参数的计算。在本文大空间γ辐射场计算中，采用体积修正的GVR-F方法计算的FOMG因子比直接模拟提高39倍，基于误差的标准差σ降低约两个数量级，优于LVR方法。针对非计数区内粒子过度分裂的情况，提出了非计数区修正方法。采用非计数区修正后，GVR-F方法计算的FOMG因子在体积修正的基础上进一步提高40%。

在体积修正和非计数区修正的基础上，本文将基于粒子误差、权重、径迹、数目、能量、碰撞和通量共7种GVR方法应用于大空间γ辐射场环境计算。基于权重的GVR-W方法模拟得到的FOMG因子比直接模拟提高约2304倍，在所有GVR方法中减方差效果最好。其他GVR方法模拟得到FOMG因子也能提高两个多量级。在基于通量的GVR-F方法上，引入光滑因子SI，模拟计算的FOMG因子随SI的增加先增加后减小。当SI=0.8时，该方法模拟得到的FOMG因子最大，达到14.78，比直接模拟提高3246倍。对于常见的全局计算算例，存在一个合适的光滑因子SI，使得模拟计算效率最高。

本文提出的对GVR方法的体积修正和非计数区修正方法具有通用性，不仅仅适用于大空间γ辐射场计算场景，在其他非等体积计数栅元/网格和含非计数区的计算场景中同样适用。

作者贡献声明刘利负责计算模型开发，程序运行，结果分析，文章撰写；左应红负责参与计算模型开发，结果分析，文章修改；牛胜利负责参与计算模型开发，指导结果分析，文章修改；朱金辉负责提供总体技术指导，文章修改；商鹏负责参与结果分析；王学栋负责参与计算模型开发。