Monte Carlo方法的最优源项偏倚抽样密度函数
苏健, 曾志, 程建平, 李君利
清华大学 工程物理系, 粒子技术与辐射成像教育部重点实验室, 北京 100084
曾志, 副研究员, E-mail:zengzhi@tsinghua.edu.cn

作者简介: 苏健(1987—), 男(汉), 福建, 博士研究生。

摘要

Monte Carlo粒子输运中的源项偏倚抽样方法可以减小方差、提高计算效率。该文通过建立一个多区域分权重数学投篮模型,模拟了输运过程中的源项信息,得到了源项偏倚抽样方差最小时的最佳抽样密度函数解析式。采用随机数值方法对模型进行了计算,验证了函数的正确性,并举一例实际的粒子输运模拟问题,表明最优偏倚抽样方法对减小方差的效果显著。该方法可作为一种普适的减方差技巧应用于Monte Carlo粒子输运中,可用于构造粒子源参数(如位置、发射方向等)的最佳偏倚密度函数,尤其在分层抽样时能给出方差最小的最优各层比例系数。

关键词: Monte Carlo方法; 源项偏倚抽样; 减方差; 最优偏倚密度函数; 分层抽样
中图分类号:TL99 文献标志码:A 文章编号:1000-0054(2014)02-0164-08
Optimal source biased sampling density function for the Monte Carlo method
Jian SU, Zhi ZENG, Jianping CHENG, Junli LI
Key Laboratory of Particle and Radiation Imaging of Ministry of Education, Department of Engineering Physics, Tsinghua University, Beijing 100084, China
Abstract

The source biased sampling method reduces the variance and improves the efficiency of Monte Carlo particle transport calculations. This paper gives a multi-region, multi-weight shooting model to simulate the source in transport processes. The density function gives the minimum variance for source biased sampling. The model is solved using a random numerical method to verify the correctness of the function. A particle transportation problem is then simulated to show the significant effect of the variance reduction. This method can be used as a general variance reduction technique in Monte Carlo particle transport analyses to construct the density function for biased source sampling for various particle source parameters, such as the transmission location and direction. It gives the best partition coefficient with the minimum variance in stratified sampling.

Keyword: Monte Carlo method; source biasing sampling; variance reduction; optimal bias density function; stratified sampling

Monte Carlo (MC)方法是一种随机实验方法,其基本原理是利用大量随机实验的平均结果来近似所要求解的问题的解[1]。由于MC方法与所求解的问题维数无关,收敛速度可以预测,因此在粒子输运相关的物理研究中得到广泛应用[2,3]。在应用MC的过程中,为了得到较精确的模拟结果,统计方差要尽可能地小; 同时,为了提高计算效率,计算时间要尽可能地少,但通常研究以前者为研究重点,所研究的方法也称为“减方差技巧”[4,5,6]。当前,已经研究出了多种减方差技巧,如区域粒子分裂-轮盘赌[7]、 Dxtran球[7,8]、 相关抽样[8,11]、自动重要抽样[12]、指数变换[13,14,15]、权重窗[14]、方向偏移[15]等,这些技巧已经在MCNP、 FLUKA、 GEANT4等大型MC程序中得到实现[16,17,18,19,20,21]。在实际应用中,经常使用源偏倚抽样方法,通过改变源项的分布函数,对贡献大的区域多抽样,对贡献小的区域少抽样,再以纠偏权重使统计结果无偏,从而在计算时间基本相同的情况下实现减小方差的目的。本质上,源偏倚抽样属于重要抽样方法,已有多篇文献给出了重要抽样的最佳重要性函数[4,5],但是这些函数都是理论上的“零方差”表达式,在实际的模拟中很难直接应用。本文以投篮数学模型为基础,得到了源偏倚抽样条件下使方差最小的密度函数,它对于源抽样的各类参数(如位置、发射方向等)具有普适性,是可以应用于MC粒子输运的减方差技巧。

1 问题的提出

在通量计算中,一个源粒子不论经过多么复杂的输运过程,最终要么被计数器接受,要么未被计数器接受(对于多道计数器来说,实际上每一道都是一个计数器,只接受某个能量范围内的粒子)。这与投篮问题有相似之处。如果在含有不同位置权重的区域投篮,可以与粒子输运问题中含有面积权重、体积权重的源项区域相类比。下面以一个多区域的投篮问题说明这个问题。

1.1 多区域投篮问题的概率论描述

某运动员进行投篮水平测试。篮球场有A、 B、 C 3个区域,面积比为0.6∶0.3∶0.1。下面用MC模拟的方法求在整个篮球场上平均命中率(即数学期望), 并给出模拟结果的均方差。设模拟的投篮总次数为 N=30 000, 且已知投篮者在3个区域的命中率大约分别为0.1、 0.3、 0.8。

这是一个含有2个随机变量的问题。记投篮区域为 X, 命中与否为 Y, 则 X的可能取值有{A,B,C}, Y的可能取值有{1,0}。并且,已知在3个区域的命中率,这是条件概率,即

PY|X(y=1|x=A)=0.1,PY|X(y=1|x=B)=0.3,PY|X(y=1|x=C)=0.8.

3个区域的面积之比是 X的边际概率,即

PX(x=A)=0.6,PX(x=B)=0.3,PX(x=C)=0.1.

根据条件概率公式可以求出所有的联合概率,如

PXY(Y=1,X=A)=PY|X(Y=1|X=A)·PX(X=A)=0.1×0.6=0.06.

表1中列出了 X Y联合概率。

表1 X Y的联合概率表

最终要求的是 Y的数学期望:

Ey=Y=01YPYY

Y=01YX=ACPX,Y=

0×(0.54+0.27+0.02)+

1×(0.06+0.09+0.08)=

0×0.77+1×0.23=0.23.

在未采用偏倚抽样时,各区域投篮次数实际上是按照3个区域的自然比例0.6∶0.3∶0.1进行分配的。而进一步要做的是,通过偏倚抽样,改变投篮次数在3个区域的分配比例,使方差达到最小。为此,需要先对一般化的情形进行推导和证明。

1.2 偏倚抽样期望、二阶矩和方差的数学推导

考虑一般性的数学描述,二元联合分布的数学期望可以用二重积分表达,

Eg(x,y)=x1x2y1y2g(x,y)f(x,y)dydx.

引入新的密度函数 f*( x, y), 根据条件概率公式,有

f(x,y)=fY|X(y|x)fX(x),f*(x,y)=fY|X(y|x)fX*(x).

其中: fX*( x)就是要寻找的最优偏倚密度函数。则期望(一阶矩)可以表达为

Eg(x,y)=x1x2y1y2g(x,y)f(x,y)dydx=x1x2y1y2g(x,y)f(x,y)f*(x,y)f*(x,y)dydx=x1x2y1y2g(x,y)fY|X(y|x)fX(x)fY|X(y|x)fX*(x)f*(x,y)dydx=x1x2y1y2g(x,y)fX(x)fX*(x)f*(x,y)dydx.

因此,二阶矩为

Eg2(x,y)=x1x2y1y2g(x,y)fX(x)fX*(x)2·f*(x,y)dydx=x1x2y1y2g(x,y)fX(x)fX*(x)2·fY|X(y|x)fX*(x)dydx=x1x2y1y2g2(x,y)fX2(x)fY|X(y|x)fX*(x)dydx.

方差为

Dg(x,y)=Eg2(x,y)-(Eg(x,y))2.

由于 Eg( x, y)不随偏倚抽样而改变,因此减小方差就是要减小二阶矩。二阶矩可化为

Eg2(x,y)=x1x2y1y2g2(x,y)fX2(x)fY|X(y|x)dyfX*(x)dx.

其中: 分子的被积式 g2( x, y) fX2( x) fY|X( y|x)都是固定的函数,而分母 fX*( x)正是要寻找的最优偏倚密度函数。

2 计算方法
2.1 最优偏倚抽样密度函数的求解及数学证明

为了解出使二阶矩 Eg2(x,y)达到最小时 fX*( x)的表达式,需要利用泛函这一数学工具。记

h(x)=y1y2g2(x,y)fX2(x)fY|X(y|x)dy,f(x)=fX*(x).

则二阶矩 Eg2(x,y)表示为

Eg2(x,y)=x1x2h(x)f(x)dx.

令泛函 J[ f( x)]等于二阶矩 Eg2(x,y), 记为

J[f(x)]=x1x2h(x)f(x)dx.

再令

Fx,f(x)=h(x)f(x),ϕ(x)=f(x).

积分约束条件为

x1x2ϕ(x)dx=x1x2f(x)dx=1.

则二阶矩 Eg2(x,y)的极小化问题可化为如下的“含积分约束条件的泛函极值问题”,

Jf(x)=x1x2Fx,f(x)dx,x1x2φ(x)dx=1.

解此问题的方法为“Lagrange乘数法”,通过引入Lagrange乘数 λ, 把“含积分约束条件的泛函极值问题”转化成新泛函 J*的“无条件极值问题”,

J*[f(x),λ]=x1x2F*[x,f(x),λ]dx.

其中,

F*[x,f(x),λ]=F[x,f(x)]+λϕ(x).

根据“泛函的Euler方程”,当泛函取极值时满足

F*f-ddxF*f'=0.

由于 F*并不含 f'( x)项,因此 ddxF*f' =0, 因此有

F*f=0,F*f=h(x)f(x)+λf(x)f=-h(x)f2(x)+λ=0.

于是,

再根据积分约束条件得到

解得

因此使泛函取极值时,

这就是 fX*( x)的表达式。随后,将 h(x)=y1y2g2(x,y)fX2(x)fY|X(y|x)dy代回,最优偏倚密度函数即为

2.2 将最优偏倚抽样密度函数应用于多区域投篮问题

第2.1节求出的 fX*( x)的表达式是一般性的结果,现将其应用到投篮问题中。投篮问题是一个单抽样变量( X)、 单计数变量( Y)的求期望问题。对于投篮问题, g( x, y)的表达式很简单,

g(x,y)=y.

而且 y的取值只有0或1, 将 fX*( x)的表达式中分子的积分变为求和,

这里的 PY|X( y=1 |x)正是表示“在 X区域投中的概率”,是已知条件。于是,

两边同乘以d x即得

再将分母的积分变为求和

至此即可直接求出:

这就是区域离散情况下的最优抽样分配比例。此时,在统计投中次数时,每次投中对计数的贡献就不再是1, 而是相应的纠偏计数权重 w。纠偏权重的表达式是

w(x)=PX(x)PX*(x).

表2和3分别列出了无偏倚抽样和最优偏倚抽样时,投篮问题的均方差计算结果。

表2 未采用偏倚抽样时投篮问题的均方差
表3 采用最优偏倚抽样时投篮问题的均方差

可见采用源项最优偏倚抽样后,方差确实显著地减小了。并且,第2.1节的推导过程已经证明了,该投篮问题采用上述的源项偏倚抽样是全局最优的(即均方差最小只能达到0.379 2)。

利用Matlab画出图1所示的函数图。 x1 x2表示区域A和区域B的分配比例,区域C的比例则必为1 -x1 -x2, 以均方差为竖坐标,画出均方差随 x1 x2变化的3维曲面图。可以看到,当 x1 =0 .427 8, x2 =0 .370 5时,均方差确实是全局最小的。

图1 方差σ随分配比例x<sub>1</sub>、x<sub>2</sub>的变化趋势

3 源项最优偏倚抽样在实际粒子输运模拟中的实例分析

下面举一个实际的粒子输运模拟问题实例说明源项最优偏倚抽样的应用。

一个空心铅球,外半径 R2 =70 cm, 内半径 R1 =20 cm, 铅材料中均匀含有40K。设定中心半径 R0 =10 cm范围内为计数区域,用MC模拟求40K的1 .461 MeV光子未发生作用而直接进入计数区域的概率。

像这样的粒子输运问题,根据先验知识,可以对计数概率有一个大致的预测。根据光子与物质作用规律,铅中某个位置的源光子,进入计数区域主要受到两个因素的影响: 1) 它与计数区域之间隔着多厚的铅, 2) 发射方向是否指向计数区域。换句话说,假如对于给定源位置、源发射方向的情形,可以估计它的计数概率,那么就可以认为计数条件概率 PY|X( y=1 |x)是“已知”的,可以把它作为源项最优偏倚密度函数的计算依据。

下面只考虑对源位置进行偏倚抽样,因此发射方向始终为各向同性。

设铅的自吸收因素、立体角因素导致的计数贡献概率减弱因子分别为 p1( r)、 p2( r)。

p1(r)=e-μ(r-R1)=eμR1·e-μr=C1e-μr.

其中, μ为1 .461 MeV光子在铅中的线衰减系数,查表得知 μ=0 .435 3 cm -1

p2(r)=Ω计数区域4π=S球冠4πr2=2πrh4πr2=h2r=r02/r2r=C21r2,p(r)=p1(r)p2(r)=C3e-μrr2.

即计数变量 Y的条件概率为

PY|R(y=1|r)=p(r)=C3e-μrr2.

对于球壳内均匀体积源,抽样半径 r本应服从的概率密度函数为(如图2所示)

fR(r)=3R23-R13r2, (R1<r<R2).

图2 未偏倚抽样时源位置半径的概率密度函数

进行最优偏倚抽样要服从的概率密度函数为

其中,

对于这个 fR*( r), 有分层抽样和连续函数抽样两种操作方法。

3.1 用分层抽样的方法进行源项最优偏倚抽样

对于源项抽样不方便自主编写的MC程序(如MCNP), 可以不算出 f( r)归一化系数的表达式,只需根据

fR*(r)r·e-μ2r

来分配抽样粒子数,然后归一化即可。不妨将铅球分为10层球壳,层号作为变量记为 X, 每层厚 5 cm, 在每层内部按体积均匀抽样,取每层的半径中间值作为名义半径 ri来计算各层的源粒子分配比例 PX*( x)。

设第 i层的名义分配比例(未归一化)为

fi=ri·e-μ2ri.

fi归一化即得到模拟粒子数份额,

各层的体积份额 vi为该层体积与总体积的比值,再算出计数权重,

wi=viPX*(x=i).

计算结果如表2所示,概率密度函数及纠偏权重函数如图3和4所示。

表2 分层偏倚抽样参数

图3 分层偏倚抽样时源位置半径的概率密度函数

图4 分层偏倚抽样时的纠偏权重函数

应当指出的是,分层抽样只是对最优的偏倚抽样的近似,层数越多,分得越细,减方差效果就越好。

3.2 用连续的密度函数进行源项最优偏倚抽样

对于源项抽样可以自主编写的MC程序(如GEANT4), 可以直接采用 fR*( r)对半径 r进行偏倚抽样,需先求出归一化系数 的表达式。

, 根据 R1R2fR*(r)dr=1可求出

C4=1R1R2r·e-μ2rdr=12070r·e-0.43532rdr=0.6875.

即偏倚抽样使用的密度函数为(函数图像如图5所示)

fR*(r)=C4·r·e-μ2r.

图5 连续偏倚抽样时源位置半径的概率密度函数

纠偏权重为(函数图像如图6所示)

图6 连续偏倚抽样时的纠偏权重函数

3.3 计算结果比较

利用GEANT4建模并计算,分别采用无偏倚抽样、分10层偏倚抽样、连续型偏倚抽样进行模拟,结果如表3所示,其中的不确定度为95%置信水平。

表3 3种偏倚抽样类型模拟结果

可见,在本计算实例中,采用源项位置最优偏倚抽样的方法,不论是分层偏倚还是连续型偏倚,均方差估计值都能降低到无偏倚时的约20%, 而且连续型偏倚比分层偏倚的均方差更小。

4 结 论

大型几何MC粒子输运问题中,源项区域往往是非常大的面源或者体源,而探测器相对于源项区域分布而言,往往很小。如果采用常规MC粒子输运方法对源项区域分布进行抽样,效率很低,而且不能保证位置抽样均匀性; 若要保证位置抽样均匀性,则计算量会非常大,同时对最终模拟结果的改善也不显著。本文研究了一种源项位置最优偏倚抽样方法,通过建立数学模型,进行了理论推导和模拟验证。结果表明,当源项各区域的计数贡献概率可以估计时,如果采用本文提出的源项位置最优偏倚抽样,方差可以达到源项偏倚意义上的最小值。该源项位置最优偏倚抽样方法能够直接应用在大型几何MC粒子输运问题中。

The authors have declared that no competing interests exist.

参考文献
[1] 裴鹿成. 蒙特卡罗方法中的若干问题[J]. 原子能科学技术, 1963(6): 422-431.
PEI Lucheng. Several issues in Monte Carlo method[J]. Atomic Energy Science and Technology, 1963(6): 422-431. (in Chinese) [本文引用:1] [CJCR: 0.426]
[2] 董秀芳. 蒙特卡罗方法及其基本特点[J]. 原子能科学技术, 1978(3): 277-289.
DONG Xiufang. Monte Carlo method and its basic characteristics[J]. Atomic Energy Science and Technology, 1978(3): 277-289. (in Chinese) [本文引用:1] [CJCR: 0.426]
[3] 朱辉, 刘义保, 游运. 蒙特卡罗方法与拟蒙特卡罗方法的历史、现状及展望[J]. 东华理工大学学报: 自然科学版, 2010, 33(4): 357-362.
ZHU Hui, LIU Yibao, YOU Yun. Monte Carlo method and quasi-Monte Carlo method[J]. Journal of East China Institute of Technology: Natural Science, 2010, 33(4): 357-362. (in Chinese) [本文引用:1] [CJCR: 0.607]
[4] 裴鹿成, 张孝泽. 蒙特卡罗方法及其在粒子输运问题中的应用 [M]. 北京: 科学出版社, 1986.
PEI Lucheng, ZHANG Xiaoze. Monte Carlo Methods and Application in Particle Transportation Problem [M]. Beijing: Science Press, 1986. (in Chinese) [本文引用:2]
[5] 裴鹿成. 计算机随机模拟 [M]. 长沙: 湖南科学出版社, 1989.
PEI Lucheng. Computer Random Simulation [M]. Changsha: Hunan Science & Technology Press, 1989. (in Chinese) [本文引用:1]
[6] 朱永生. 实验物理中的概率和统计 [M]. 北京: 原子能出版社, 2006.
ZHU Yongsheng. Probability and Statistics in Experimental Physics [M]. Beijing: Atomic Energy Press, 2006. (in Chinese) [本文引用:1]
[7] 张崭, 李君利, 武祯, . 内照射小器官剂量计算中减方差技巧的比较和应用[J]. 清华大学学报: 自然科学版, 2007, 47(S1): 1051-1056.
ZHANG Zhan, LI Junli, WU Zhen, et al. Comparison and application of variance-reduction techniques used in internal radiation dose calculations for small organs[J]. Journal of Tsinghua University: Science and Technology, 2007, 47(S1): 1051-1056. (in Chinese) [本文引用:2] [CJCR: 0.609]
[8] 武祯, 李君利. 用于探测器校正因子计算的Monte Carlo方法[J]. 清华大学学报: 自然科学版, 2006, 46(9): 1585-1588.
WU Zhen, LI Junli. Monte Carlo method for calculating particle radiation detector correction factors[J]. Journal of Tsinghua University: Science and Technology, 2006, 46(9): 1585-1588. (in Chinese) [本文引用:2] [CJCR: 0.609]
[9] Bielajew A F. Correction factors for thick-walled ionization chambers in point-source photon beams[J]. Phys Med Biol, 1990, 35(4): 501-516. [本文引用:1] [JCR: 2.922]
[10] Hedtjarn H, Carlsson G, Williamson J F. Accelerated Monte Carlo-based dose calculations for brachytherapy planning using correlated sampling[J]. Phys Med Biol, 2002, 47(3): 351-376. [本文引用:1] [JCR: 2.922]
[11] Buckley L A, Kawrakow I, Rogers D W O. CSnrc: Correlated sampling Monte Carlo calculations using EGSnrc[J]. Med Phys, 2004, 31(12): 3425-3435. [本文引用:1] [JCR: 3.012]
[12] 武祯, 李君利, 程建平. 一种改进的计算探测器校正因子的相关抽样方法[J]. 高能物理与核物理, 2006, 30(8): 771-775.
WU Zhen, LI Junli, CHENG Jianping. An improved correlated sampling method for calculating correction factor of detector[J]. High Energy Physics and Nuclear Physics, 2006, 30(8): 771-775. (in Chinese) [本文引用:1] [JCR: 0.233] [CJCR: 0.281]
[13] 许淑艳. 蒙特卡罗方法在实验核物理中的应用 [M]. 北京: 原子能出版社, 1996.
XU Shuyan. Monte Carlo Method in Experimental Nuclear Physics [M]. Beijing: Atomic Energy Press, 1996. (in Chinese) [本文引用:1]
[14] 邱睿, 李君利, 曾志. 邮件辐照系统屏蔽的Monte Carlo计算[J]. 清华大学学报: 自然科学版, 2004, 44(3): 297-300.
QIU Rui, LI Junli, ZENG Zhi. Monte Carlo calculation of the shielding of mail irradiation systems[J]. Journal of Tsinghua University: Science and Technology, 2004, 44(3): 297-300. (in Chinese) [本文引用:1] [CJCR: 0.609]
[15] 王汝赡, 姜宏宇. 方向偏移法及其在光子输运模拟中的应用[J]. 核电子学与探测技术, 1998, 18(3): 177-181.
WANG Rushan, JIANG Hongyu. Direction biasing method and its application in γ ray transportation simulation[J]. Nuctear Eleetronics & Detection Technology, 1998, 18(3): 177-181. (in Chinese) [本文引用:2]
[16] Olsher R H. A practical look at Monte Carlo variance reduction methods in radiation shielding[J]. Nuclear Engineering and Technology, 2006, 38(4): 225-230. [本文引用:1] [JCR: 0.758]
[17] MCNP Monte Carlo Team, X-5. MCNP5_RSICC_1. 30, LA-UR-04. 5921 [R]. Los Alamos, NM: Los Alamos National Laboratory, 2004. [本文引用:1]
[18] Hendricks J S, McKinney G W, Waters L S, et al. MCNPX, Version 2. 5. 3, LA-UR-04-0569 [R]. Los Alamos, NM: Los Alamos National Laboratory, 2004. [本文引用:1]
[19] Booth T E. A Sample Problem in Variance Reduction in MCNP, LA-10363-MS [R]. Los Alamos, NM: Los Alamos National Laboratory, 1985. [本文引用:1]
[20] Fassò A, Ferrari A, Ranft J, et al. FLUKA: Performances and applications in the intermediate energy range [C]//Proc 1st AEN/NEA Specialists' Meeting on Shielding Aspects of Accelerators, Targets and Irradiation Facilities (SATIF 1). Arlington, TX: OECD Documents, 1995: 287-304. [本文引用:1]
[21] Allison J, Amako K, Apostolakis J, et al. GEANT4 developments and applications[J]. IEEE Transactions on Nuclear Science, 2006, 53(1): 270-277. [本文引用:1] [JCR: 1.455]