关志远 绪梅
陆军防化学院 北京 102205
摘要：本文首先阐述了蒙特卡罗方法的相关概述，并详细分析了相关抽样方法在中子、光子输运中的应用。
关键词：蒙特卡罗方法；相关抽样；中子、光子输运
在计算结果本身精度不高的情况下，相关抽样法可使计算结果差的精度足够高，充分反映探测器响应的变化。因此，采用相关抽样法研究样品组成变化时探测器响应谱的变化，能节省大量的计算时间，同时也可在一次计算中得到多个样品结果，从而提高了计算效率。
一、蒙特卡罗方法概述
蒙特卡罗方法又称统计模拟法、随机抽样技术，是一种随机模拟方法，以概率和统计理论方法为基础的一种计算方法，是使用随机数(或更常见的伪随机数)来解决很多计算问题的方法。将所求解的问题同一定的概率模型相联系，用电子计算机实现统计模拟或抽样，以获得问题的近似解。为象征性地表明这一方法的概率统计特征，故借用赌城蒙特卡罗命名。
该方法于20世纪40年代美国在第二次世界大战中研制原子弹的“曼哈顿计划”计划的成员S.M.乌拉姆和J.冯·诺伊曼首先提出。数学家冯·诺伊曼用驰名世界的赌城-摩纳哥的Monte Carlo-来命名这种方法，为它蒙上了一层神秘色彩。在这之前，蒙特卡罗方法就已存在。1777年，法国Buffon提出用投针实验的方法求圆周率∏。这被认为是蒙特卡罗方法的起源。
由概率定义知，某事件的概率可用大量试验中该事件发生的频率来估算，当样本容量足够大时，可认为该事件的发生频率即为其概率。因此，可先对影响其可靠度的随机变量进行大量的随机抽样，然后把这些抽样值一组一组地代入功能函数式，确定结构是否失效，最后从中求得结构的失效概率。蒙特卡罗法正是基于此思路进行分析的。此外，该方法回避了结构可靠度分析中的数学困难，不管状态函数是否非线性、随机变量是否非正态，只要模拟的次数足够多，就可得到一个比较精确的失效概率和可靠度指标。
二、相关抽样
1、相关抽样方法原理。在具体计算中，按通常方式计算一个样品(称为参考样品)，同时，计算另外一组样品，这些样品与参考样品相比在密度或元素成分上稍有不同(称为相关样品)。在粒子输运过程中，强迫各相关样品中输运的粒子与参考样品中输运的粒子经历同样的物理过程(径迹、发生碰撞的元素)，并通过调整相关样品中粒子的权重来修正相关样品中的结果。这样，一方面节省了大量的时间(各元素微观截面的处理，几何位置的处理，反应元素的处理，反应类型的处理等对多个样品只需一次)，与所有样品单独计算相比
效率大幅提高；另一方面，这些样品的计算过程呈正相关，对样品发生微小改变造成的探测器响应差异的计算结果可获得较高的精度。
2、对相关样品中粒子权重的调整。由于相关样品的抽样不是采用它真正的概率密度函数，而是根据参考样品的概率密度函数来进行的，因此必须对相关样品中粒子的权重做相应调整，以保证结果无偏。
1)抽样自由程时的调整()。设D为沿中子运动方向到达均匀介质边缘的距离，为参考样品中抽样得到的中子的自由飞行距离，权重调整因子()为：
A、＜D时
.png)
式中：-为参考样品的宏观总截面；-为相关样品的宏观总截面；-为参考样品中第种元素的宏观总截面；-为相关样品中第种元素的宏观总截面。
3)抽样反应类型时的调整。确定碰撞元素后，由于中子与此元素发生何种反应只与该元素本身有关，和样品的成分等其他参数无关，因此权重无需调整。
以上三个权重调整方式适应于中子和光子两种情况。
4)中子反应放出光子时的调整()。为高效率的获得库谱，在产生光子时采用了期望值技术，中子与任何元素反应时所需元素的光子都放出。比如，要计算元素H的库谱，则在中子历史中，不管中子实际上与何种元素反应，都假定中子与H元素反应且放出光子并通过调整放出光子的权重来保证结果的无偏。这样，权重的调整包括了两部分：一是针对期望值技术的调整，二是针对相关样品的调整。权重因子(m)的调整与上述类似。
5)采用方差减小技术时的调整()。这种情况下权重因子的调整与采用何种方差减小技术有关。若使用的方差减小技术与具体的样品成分无关，则无需对相关样品的权重进行调整，否则视具体情况调整。
通过对相关样品中中子、光子的权重进行调整，最终可获得无偏的结果。
3、计算结果的记录。假定对某个中子历史的结果进行记录时，参考样品中子的权重为W，则相关样品中中子的权重为
.png)
前一公式可获得相关样品的结果，后一公式可获得参考样品和相关样品的结果的差值。由于这两种样品的计算过程呈正相关，因此二者差值的方差也可达到较好的结果。此外，光子历史结果的记录与中子类似。
三、相关抽样的计算实例
一个为计算参数发生微小变化时系统结果的变化量；另一个为采用相关抽样法，在一次计算中同时获得两种样品(两种样品的成分相差不大时)的C元素的库谱。使用程序为NPTrams程序。计算模型：
①取半径为30cm的煤球，假定煤中元素均匀分布.
②假定中子源为241Am-Be点源，且中子源位于煤球的球心。
③对光子的计算结果未经探测器[NaI(Tl)]响应展宽(并非实际的库谱)。
计算所用的参考煤样成分为：H元素质量百分数为5.00；C元素质量百分数为69.3；N元素质量百分数为1.00；O元素质量百分数为17.00；Na元素质量百分数为0.20；Mg元素质量百分数为0.02；Al元素质量百分数为3.00；Si元素质量百分数为2.30；S元素质量百分数为1.00；Cl元素质量百分数为0.03；K和Ti元素质量百分数为0.05；Ga元素质量百分数为0.50；Fe元素质量百分数为0.55。
实例一：计算煤密度由0.67552g/cm3变为0.68228g/cm3(煤密度增加1%)时模型中球面出射的中子流量变化量。计算结果表明，在煤样密度增加1%时，煤球面出射的中子流量变化小。采用相关抽样及两次独立计算结果相减，两种方法求出的中子流量的变化量如图1所示。由图1可看出，采用相关抽样法得到的曲线平滑，结果明显比采用直接相减方法好。
图1 两种方法得到的中子流量的变化量
.png)
把相关抽样得到的结果作为准确值，采用直接相减法的结果完全不可信。实际上，在计算3264000个源中子的情况下，相关抽样得到的结果的相对统计涨落约为2～4%(由程序给出)，而根据相关公式估计的直接相减法结果的相对统计涨落多数大于40%，甚至达到100%以上，结果根本就没有意义。采用相关抽样法得到的最终结果为：煤样密度增加1%后，球面出射的中子流量在热中子区大约减小1%，在快中子区大约减小2%。
由此可知，采用相关抽样法后，计算效率大幅提高。
实例二：计算煤密度为1.3650g/cm3,C元素含量不变(69.3%)与相对减少5%(65.83%)两种情况下模型中C元素的库谱。
采用两种方法来计算。一是采用相关抽样法，在一次计算中同时给出两个结果(C元素含量为69.3%的样品作为参考样品)；二是单独计算两次。为验证相关抽样的正确性，单独计算采用的程序是MCNP。简单起见，表1仅给出了C元素的4种特征光子的结果，并与MCNP的结果进行了比较。
表1 采用相关抽样与单独计算结果的比较
.png)
从表中可知，采用相关抽样方法的结果与单独计算的结果相比，无差异。但采用相关抽样法在一次计算中给出了两个结果，而所用时刻与单独计算一次相比仅增加了约5%，计算效率提高了仅一倍。
综上所述，相关抽样法是一种有效的方差减小技术，其适应于计算系统参数发生微小改变时输出结果的变化量及在一次计算中同时给出多个相似系统(系统参数差异不大)的结果，尤其在第一种情况下，由于蒙特卡罗计算结果的概率性质，若不使用相关抽样法，计算结果往往被统计涨落淹没，难以达到所需精度。
参考文献：
[1]陈伯显.中子感生瞬发射线煤多元素分析研究[J].核电子学与探测技术，2015(01).
[2]李铁柱.相关抽样方法在中子、光子输运中的应用[J].核电子学与探测技术，2015(04).