第四章 蒙特卡罗方法解粒子输运问题
1,屏蔽问题模型
2,直接模拟方法
3,简单加权法
4,统计估计法
5,指数变换法
6,蒙特卡罗方法的效率
? 作 业
第四章 蒙特卡罗方法解辐射屏蔽问题
辐射（光子和中子）屏蔽问题是蒙特卡罗
方法最早广泛应用的领域之一。本章主要从物
理直观出发，说明蒙特卡罗方法解决这类粒子
输运问题的基本方法和技巧。而这些方法和技
巧对于诸如辐射传播、多次散射和通量计算等
一般粒子输运问题都是适用的。
1,屏蔽问题模型
在反应堆工程和辐射的测量与应用中, 常
常要用一些吸收材料做成屏蔽物挡住光子或中
子 。 我们所关心的是经过屏蔽后射线的强度及
其能量分布, 这就是屏蔽问题 。
当屏蔽物的形状复杂, 散射各向异性, 材
料介质不均匀,核反应截面与能量, 位置有关
时, 难以用数值方法求解, 用蒙特卡罗方法能
够得到满意的结果 。
粒子的输运问题带有明显的随机性质，粒
子的输运过程是一个随机过程。粒子的运动规
律是根据大量粒子的运动状况总结出来的，是
一种统计规律。蒙特卡罗模拟，实际上就是模
拟相当数量的粒子在介质中运动的状况，使粒
子运动的统计规律得以重现。不过，这种模拟
不是用实验方法，而是利用数值方法和技巧，
即利用随机数来实现的。
为方便起见，选用平板屏蔽模型，在厚度为 a，
长、宽无限的平板左侧放置一个强度已知，具有已知
能量、方向分布的辐射源 S 。求粒子穿透屏蔽概率
（穿透率）及其能量、方向分布。穿透率就是由源发
出的平均一个粒子穿透屏蔽的数目。
同时，假定粒子在两次碰撞之间按直线运动,且粒
子之间的相互作用可以忽略。
2,直接模拟方法
直接模拟方法就是直接从物理问题出发,
模拟粒子的真实物理过程 。
1) 状态参数与状态序列
2) 模拟运动过程
3) 记录结果
粒子在介质中的运动的状态，可用一组参数来描
述，称之为 状态参数 。它通常包括：粒子的空间位置 r,
能量 E 和运动方向 Ω，以 S＝ ( r,E,Ω ) 表示。
有时还需要其他的参数，如粒子的 时间 t 和附带
的权重 W,这时 状态参数 为 S'＝ ( r,E,Ω,t,W ) 。
状态参数 通常要根据所求问题的类型和所用的方
法来确定。
对于无限平板几何，取 S＝ ( z,E,cosα)
其中 z 为粒子的位置坐标,α为粒子的运动方向与
Z 轴的夹角。
对于球对称几何,取 S＝ ( r,E,cosθ)
其中 r 表示粒子所在位置到球心的距离,θ为粒子
的运动方向与其所在位置的径向夹角。
1) 状态参数与状态序列
粒子第 m 次碰撞后的状态参数为
或
它表示一个由源发出的粒子, 在介质中经过 m 次碰撞
后的状态, 其中
rm,粒子在第 m 次碰撞点的位置
Em,粒子第 m 次碰撞后的能量
Ωm：粒子第 m 次碰撞后的运动方向
tm,粒子到第 m 次碰撞时所经历的时间
Wm,粒子第 m 次碰撞后的权重
有时, 也可选为粒子进入第 m 次碰撞时的状态参数 。
),,( mmmm E ΩrS ?
),,,,( mmmmmm WtE ΩrS ??
一个由源发出的粒子在介质中运动, 经过若干次
碰撞后, 直到其运动历史结束 （ 如逃出系统或被吸收
等 ） 。 假定粒子在两次碰撞之间按直线运动, 其运动
方向与能量均不改变, 则粒子在介质中的运动过程可
用以下碰撞点的 状态序列 描述,
S0, S1, …, SM-1, SM
或者更详细些,用
来描述 。 这里 S0 为粒子由源出发的状态, 称为初态,
SM 为粒子的终止状态 。 M 称为粒子运动的链长 。
这样的序列称为粒子随机运动的历史, 模拟一个
粒子的运动过程, 就变成确定状态序列的问题 。
?
?
?
?
?
?
?
?
?
?
?
?
?
MM
MM
MM
EEEE
ΩΩΩΩ
rrrr
,,,,
,,,,
,,,,
110
110
110
?
?
?
为简单起见，这里以中子穿透均匀平板的模型来
说明，这时 状态参数 取 S＝ ( z,E,cosα)。
模拟的步骤如下,
(1) 确定初始状态 S0,
确定粒子的初始状态，实际上就是要从中子源的
空间位置、能量和方向分布中抽样。设源分布为
则分别从各自的分布中抽样确定初始状态。
对于平板情况,
抽样得到 z0＝ 0。
2) 模拟运动过程
)( c o s)()()c o s,,( 030201000 ?? fEfzfEzf ?
)()( 001 zzf ??
(2) 确定下一个碰撞点,
已知状态 Sm-1,要确定状态 Sm,首先要确定下一个
碰撞点的位置 zm。 在相邻两次碰撞之间, 中子的输运
长度 l 服从如下分布,
对于平板模型, l 服从分布,
其中, Σ t 为介质的中子宏观总截面,
积分 称为粒子输运的自由程数,
系统的大小通常就是用系统的自由程数表示的 。
? ?? ?????????? ?????? l mmmtmmmt ldElEllf 0 111111 ),(e x p),()( ΩrΩr
? ?? ?????????? ?????? l mmmtmmmt ldElzElzlf 0 111111 ),c o s(e x p),c o s()( ??
? ????? ???l mmmt ldEl0 111 ),( Ωr
显然, 粒子输运的自由程数服从指数分布,
因此从 f ( l ) 中抽样确定 l,就是要从积分方程
中解出 l。
对于单一介质
则 下一个碰撞点的位置
如果 zm≥ a,则中子穿透屏蔽, 若 zm≤ 0,则中子被反
射出屏蔽 。 这两种情况, 均视为中子历史终止 。
?ln),(0 111 ???????? ???l mmmt ldEl Ωr
)(
ln
1??
??
mt E
l ?
1
1
1
11
c o s)(ln
c o s
?
?
?
??
???
???
m
mt
m
mmm
Ez
lzz
??
?
(3) 确定被碰撞的原子核,
通常介质由几种原子核组成, 中子与核碰撞时,
要确定与哪一种核碰撞 。 设介质由 A,B,C 三种原子
核组成, 其核密度分别为 NA,NB,NC,则介质的宏观
总截面为,
其中 分别为核 A,B,C 的宏观总截面 。 其定
义如下,
分别表示 (·)核的宏观总截
面, 核密度和微观总截面 。
)()()()( 1111 ???? ??????? mCtmBtmAtmt EEEE
CtBtAt ???,,
)()( 1)()(1)( ????? ?? mtmt ENE ?
)()( 1)()(1)( ?????? mtmt ENE ?、、
由于中子截面表示中子与核碰撞可能性的大小,
因此, 很自然地, 中子与 A,B,C 核发生碰撞的几率
分别为,
利用离散型随机变量的抽样方法, 确定碰撞核种类,
)(
)(,
)(
)(,
)(
)(
1
1
1
1
1
1
?
?
?
?
?
?
?
??
?
??
?
??
mt
m
C
t
C
mt
m
B
t
B
mt
m
A
t
A E
EP
E
EP
E
EP
核碰撞与
核碰撞与
核碰撞与
C
BPP
AP
BA
A
??
?
?
?
＞
≤
＞
≤
(4) 确定碰撞类型,
确定了碰撞的核 (比如 B核 )后, 就要进一步确定碰
撞类型 。 中子与核的反应类型有弹性散射, 非弹性散
射, (n,2n)反应, 裂变和俘获等, 它们的微观截面分
别为
则有
各种反应发生的几率分别为
)()()()()( 111)2,(11 ????? mBcmBfmB nnmBinmBel EEEEE ????? 和、、、
)()()()()()( 111)2,(111 ?????? ????? mBcmBfmB nnmBinmBelmBt EEEEEE ??????
)()(
)()(
)()(
)()(
)()(
11
11
11)2,()2,(
11
11
??
??
??
??
??
?
?
?
?
?
m
B
tm
B
cc
m
B
tm
B
ff
m
B
tm
B
nnnn
m
B
tm
B
inin
m
B
tm
B
elel
EEP
EEP
EEP
EEP
EEP
??
??
??
??
??
利用离散型随机变量的抽样方法, 确定反应类型 。
在屏蔽问题中, 中子与核反应常只有弹性散射和
吸收两种类型, 吸收截面为,
这时, 总截面为,
发生弹性散射的几率为,
若, 则为弹性散射；否则为吸收, 发生吸收反
应意味着中子的历史终止 。
)(
)(
1
1
?
??
m
B
t
m
B
el
el E
EP
?
?
elP??
)()()( 111 ??? ?? mBcmBfmBa EEE ???
)()()( 111 ??? ?? mBamBelmBt EEE ???
(5) 确定碰撞后的能量与运动方向,
如果中子被碰撞核吸收, 则其输运历史结束 。 如
果发生弹性散射, 需要确定散射后中子的能量和运动
方向 。 中子能量 Em 为,
A是碰撞核的质量与中子质量之比, 一般就取元素的原
子量； θC 为质心系中 中子散射前后方向间的夹角, 即
偏转角 。 可从 质心系中 弹性散射角分布
fC(μC) 中抽样产生 。 实验室系散射角 θL的余弦 μL为,
? ? ? ?? ?
2
1
)
1
1(
c o s11
2
?
??
???? ?
A
Ar
rr
E
E Cmm ?
CC ?? co s?
C
C
L AA
A
?
??
21
1
2 ??
??
如果给出实验室系散射角余弦分布 fL(μL),可直接
从 fL(μL)中抽取 μL,此时能量 Em与 μL的关系式为,
确定了实验室系散射角 θL后,
再使用球面三角公式
确定 cosαm,
其中 χ 为在 ［ 0,2π ］ 上均匀分布的方位角 。
?????? c o ss i ns i nc o sc o sc o s 11 ????? ?? LmLmm
? ?22221 1)1( LLmm AAEE ?? ????? ?
至此, 由 Sm-1完全可以确定 Sm。 因此, 当中子由
源出发后, 即 S0确定后, 重复步骤 (2)～ (5),直到中子
游动历史终止 。 于是得到了一个中子的随机游动历史
S0, S1, …, SM-1, SM,即
也就是模拟了一个由源发出的中子的运动过程 。
?
?
?
?
?
?
?
?
?
?
?
?
?
MM
MM
MM
EEEE
zzzz
???? c o s,c o s,,c o s,c o s
,,,,
,,,,
110
110
110
?
?
?
以上模拟过程可分为两大步：第一步确定粒子的
初始状态 S0,第二步由状态 Sm-1来确定 状态 Sm。 这第二
步又分为两个过程：第一个过程是确定碰撞点位置 zm,
称为输运过程；第二个过程是确定碰撞后粒子的能量
及运动方向, 称为碰撞过程 。 对于中子而言, 碰撞过
程是先确定散射角, 进而确定能量和运动方向；而对
于光子, 碰撞过程是先确定能量, 再确定散射角以及
运动方向 。 重复这两个过程, 直至粒子的历史终止 。
这种模拟过程, 是解任何类型的粒子输运问题所
共有的, 它是蒙特卡罗方法解题的基本手段 。
在获得中子的随机游动历史后，我们要对所要计
算的物理量进行估计。对于屏蔽问题，我们要计算中
子的穿透率。考察每个中子的随机游动历史，它可能
穿透屏蔽（ zM≥ a），可能被屏蔽发射回来（ zM≤ 0），
或者被吸收。设第 n 个中子对穿透的贡献为 ηn,则
如果我们共跟踪了 N 个中子，则穿透屏蔽的中子数为,
3) 记录结果
?
?
?
?
??
，或者被吸收当
当
0,0
,1
M
M
n z
az?
?
?
?
N
n
nN
1
1 ?
则穿透屏蔽概率的近似值为,
它是穿透率的一个无偏估计 。
我们称这种直观地模拟过程和估计方法为直接模
拟方法 。 在置信水平 1－ α ＝ 0.95 时, 的误差为,
其中 为 ηn的均方差, 由于 ηn是一个服从二项分布的随
机变量, 所以
或
?
?
??
N
n
nN NN
NP
1
1)1( 1? ?
)1(?NP
NNPP N
??? ??? 96.1? )1( ???
??
)?1(??
)1(
)1()1(2
2
NN PP
PP
??
??
?
?
?
?
为得到中子穿透屏蔽的能量, 角分布, 将能量,
角度范围分成若干个间隔,
其中 Emax,Emin分别表示能量的上, 下限, 对于穿透屏
蔽的中子按其能量, 方向分间隔记录 。 设一穿透屏蔽
的中子能量为 EM,其运动方向与 Z轴夹角为 αM,若能
量 EM属于第 i 个能量间隔 Δ Ei,角度 αM属于第 j 个角
度间隔 Δ αj,则分别在第 i 个能量计数器及第 j 个角
度计数器中加 "1"。
20 10
m a x01m i n
???? ?????
?????
J
I EEEEE
?
?
跟踪 N 个中子后, 则
分别为穿透中子的能量分布和角分布 。 其中 N1,i 和
N2,i 分别为第 i 个能量和第 j 个角度间隔的穿透中
子数 。 归一后分别为,
Jj
N
N
P
Ii
EN
N
P
j
j
N
i
i
N
j
i
,,2,1?
,,2,1?
,2)1(
,1)1(
,2
,1
?
?
?
??
?
?
??
?
?
Jj
N
N
P
P
P
Ii
EN
N
P
P
P
j
j
N
N
N
i
i
N
N
N
j
j
i
i
,,2,1
?
?
?
,,2,1
?
?
?
1
,2
)*1(
)1(
)*1(
1
,1
)*1(
)1(
)*1(
,2
,2
,1
,1
?
?
?
??
??
?
??
??
?
3,简单加权法
从模拟物理过程来说, 直接模拟法是最简单, 也
是最基本的方法 。 但是, 在直接模拟法中, 不管中子
在屏蔽中经过多少次碰撞, 只要在介质中被吸收, 对
穿透的贡献就为零；因此在所跟踪的粒子中绝大部分
都对穿透没有贡献 。 而在许多屏蔽问题中, 穿透率的
数量级在 10-6到 10-8。 进一步, 如果我们要求穿透率达
相对误差小于 1％, 即
那么, N 要大到惊人的数量级 1010到 1012。 显然, 这时
用直接模拟法计算不是很有效 。
%1)1(|?| ???????? PNPN PPPNP PP ??? ?????
屏蔽物一般是由吸收强的介质组成，因此在每次
碰撞时，粒子很有可能被吸收而停止跟踪。现在改变
模拟方法，在判断碰撞类型 时，可以认为粒子
的 部分是弹性散射，而其余部
分被吸收，即人为地把中子分成两部分，一部分弹性
散射，一部分吸收。弹性散射这部分继续跟踪；吸收
部分则停止跟踪。也就是说，我们利用中子权重的变
化来反应继续弹性散射的部分。这就是简单加权法的
基本思想。
1) 简单加权法
elP??
)()( 11 ?? ??? mBtmBelel EEP
显然, 在加权法中中子的权重 W 已成为中子状态
参数的组成部分 。 这时, 中子历史成为,
对源中子, 取 W0=1。 经过碰撞中子权重的变化为,
因子 称为尚存因子 。
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
,,,,,
c o s,c o s,,c o s,c o s
,,,,
,,,,
110
110
110
110
MM
MM
MM
MM
WWWW
EEEE
zzzz
?
?
?
?
????
)(
)(
1
1
1
?
?
? ?
??
m
B
t
m
B
el
mm E
EWW
)()( 11 ?? ?? mBtmBel EE
这时, 第 n 个中子对穿透的贡献为,
如果我们共跟踪了 N个中子, 则穿透率 P的无偏估计为,
类似地, 可以得到 穿透中子的能量分布和角分布 。 只
不过在对各计数器进行的加 "1" 操作改为加 WM。
??
?
?
??
0,0
,
M
MM
n z
azW
当
当?
?
?
?
N
n
nN NP
1
)2( 1? ?
简单加权法的方差估计为,
与直接模拟法相比，有
注意到 ζn≤1，有
这表明简单加权法的方差小于直接模拟法的方差。这
是因为加权法比直接模拟法减少了一次随机抽样。
2) 简单加权法的方差
2)2(
1
22 )?(1?
N
N
n
n PN ?? ?
?
?? ?
?
?
???
N
n
nnN
1
222 )(1?? ????
??
22 ?? ?? ?? ?
加权法的思想在蒙特卡罗方法中用途很广泛。例
如，对于具有中子增殖反应，如裂变,(n,2n),(n,3n)
反应的中子输运问题，一个中子与核发生碰撞后，根
据反应的类型会产生不同数量的次级中子，每个次级
中子又会产生新的次级中子，这样链锁反应 下去，使
得用直接模拟法模拟每一个中子是非常困难的。这种
情况可以利用加权法来处理。
3) 权重方法的其它应用
中子与核发生碰撞 后，产生的次级中子平均数为,
这里 νf 为裂变次级中子数 。 于是, 碰撞后的权重为,
而决定碰撞类型的几率分别为,
其中
加权法的思想, 还可以应用到连续分布情况和偏
倚抽样的问题
t
ffnnnninel
?
??????? ?????????? )3,()2,( 3211
??? ?1mm WW
t
ff
t
nn
t
nn
t
in
t
el
?
??
?
?
?
?
?
?
?
? ???,3,2,,)3,()2,(
ffnnnninelt ??????? ?????????? )3,()2,( 3211
4,统计估计法
加权法虽然改进了直接模拟法, 但它同样只关心
中子是否穿透屏蔽这一信息, 因此对每一个中子历史
的信息利用得很不充分 。 统计估计法能够较多地利用
中子的历史信息, 因而能得到更好的结果 。
一个中子, 可能在介质内不发生碰撞而直接穿透
屏蔽, 也可能在介质内发生一次碰撞后再穿透屏蔽,
或经过二次碰撞穿
透屏蔽, 等等,
这些事件是互不相
容的, 因此穿透概
率 P 可表示为,
其中 Pm 是中子恰好经过 m 次碰撞而穿透屏蔽的概率 。
这表明, 可以用求 Pm (m=0,1,… ) 的方法得到 P。 这样,
中子对穿透概率的贡献就不只限于末次碰撞了 。
??
?
?
0m
mPP
α0 α1 S0
S1
Sm αm
P1
P0
Pm
Z a 0
设中子的历史为,
根据该中子的历史, 我们可以估计出中子恰好经
过 m 次碰撞后, 穿透屏蔽的部分
显然, 具有初态 S0＝ ( 0,E0,cosα0,W0 ) 的中子, 未
经碰撞直接穿透的部分是,
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
,,,,,
c o s,c o s,,c o s,c o s
,,,,
,,,,
110
110
110
110
MM
MM
MM
MM
WWWW
EEEE
zzzz
?
?
?
?
????
1,,1,0,? ?? MmP m ?
?
?
?
?
?
? ???
0
000 c o s)(e x p
?
?
aEWP
t
类似地, 在经过了第 m 次碰撞后的中子具有状态
Sm＝ ( zm,Em,cosαm,Wm ), 其可能穿透的部分, 正好是
一个中子恰好经过 m 次碰撞穿透的部分,
这里的这种估计技巧, 由于是对每次碰撞后的状
态, 求其后未经碰撞直接穿透的贡献, 因此该方法也
称为最后自由飞行估计 。
1,,1,0
0
0c o s
c o s
)(e x p?
??
?
?
?
?
?
?
?
?
?
?
?
? ?
??
?
Mm
za
EW
P m
m
m
mtm
m
?
其它，
,?
?
于是得到该中子对穿透的贡献,
如果我们共跟踪了 N个中子, 则穿透率 P的估计为,
其方差估计为,
??
?
?
1
0
?? M
m
mPP
?
?
?
N
n
N nPNP
1
)3( )(?1?
? ? ? ?2)3(
1
22 ?)(?1)?(
N
N
n
PnPNP ?? ?
?
?
在直接模拟方法中，相对误差为
其中 为与置信水平 1－ α相应的量。
如果构造一个新的概率模型，使得该模型的穿透
率 P*与原模型的穿透率 P之间存在关系,
使用直接模拟方法,相对误差为
5,指数变换法
PNPN
PP
PN ??
??? ??? ????? )1(
??
PKP ??*
*
*
PN ?
? ???
如果令 ε*＝ ε,即
这意味着, 达到同样的相对误差, 跟踪粒子的数
目缩小 K 倍, 从而减少 K 倍的计算量 。 指数变换法就
是构造一个新的概率模型的一个有效方法 。
K
N
N
PKNPNPN
?
??
?
?
?
?
*
***
??? ???
构造如下伪过程：宏观总截面为
散射截面仍为 Σel(E)。 其中 Emin,Emax 分别为能量的下
限和上限, α为粒子的运动方向与 Z 轴的夹角 。 可以证
明这个伪过程的穿透概率 P* 与原过程的穿透概率 P之
间有如下关系,
显然, 。 因伪过程与原过程的结果相差 e 指
数, 所以该方法称为指数变换法 。
)(m i n
c o s)()(
m a xm i n
*
EC
CEE
tEEE
tt
??
?????
??
?
PeP aC ?? ?*
1?? ?aCeK
分析一下伪过程的定义,可以明显看出 P*增大的原
因 。 当 cosα=1 时, 粒子运动方向与 Z 轴方向一致, 其
截面最小, 粒子沿 Z 轴方向输运的距离较远；而当
cosα=－ 1 时, 粒子运动方向背向 Z 轴方向, 这时其截
面最大, 粒子向后输运的距离较短 。 因此, 截面变换
的结果是加强了粒子向前运动的能力, 因而使穿透概
率增大 。
伪过程的构造与几何形状及所考虑的问题有关 。
比如, 对球形几何, 使用指数变换法求穿透概率时,所
构造的宏观总截面与平板屏蔽的情况不同, 粒子的模
拟方法也较复杂 。
6,蒙特卡罗方法的效率
衡量一种蒙特卡罗技巧的好坏, 除了看其方差大
小外, 还要看其所需费用 （ 计算时间 ） 多少, 即从该
技巧的效率 Ef（ 方差与费用乘积的倒数 ） 全面考虑,
其中 σ2 为方差, T 为所需费用 。 Ef 大时, 所用方
法的效率高；否则, 效率低 。
在一般情况下, 有些方法虽然减小了方差, 却增
加了费用 。 例如, 加权法, 统计估计法虽然较直接模
拟方法减小了方差, 却使每个粒子的运动链长增加,
或记录贡献的计算时间增加 。 因 此, 不能认为方差小
的方法一定好, 要从方法的效率全面考虑 。 在有些情
况下, 直接模拟方法仍然是一个被广泛使用的方法 。
TE f ?? 2
1
?
? 作 业
1) 光子散射后能量分布的抽样
把光子散射能量分布改写成如下形式进行抽样,
??
?
?
?
??
?
?
?
?
?
??
?
? ??
?
?
?
?
?
?
?
?
??
?
??
?
? ???
22
2 11111
)(
1)(
xxxx
x
K
xf
?
?
?
?
在 ［ 1,1+2α］ 上定义如下函数,
3
2
2
2
1
2
21
)1(
)(
2
)(
1
1
)21)((
2
)(
2
1
)(
1
2
21
)(
x
x
K
xH
x
K
xH
xf
x
xf
?
?
?
?
?
?
?
?
?
?
??
?
?
?
?
? ??
?
?
?
?
?
?
?
?
?
?
?
??
?
?
?
?
?
?
?
? 作 业
1) 给出分布密度函数
的抽样方法。
＞
≤
1,屏蔽问题模型
2,直接模拟方法
3,简单加权法
4,统计估计法
5,指数变换法
6,蒙特卡罗方法的效率
? 作 业
第四章 蒙特卡罗方法解辐射屏蔽问题
辐射（光子和中子）屏蔽问题是蒙特卡罗
方法最早广泛应用的领域之一。本章主要从物
理直观出发，说明蒙特卡罗方法解决这类粒子
输运问题的基本方法和技巧。而这些方法和技
巧对于诸如辐射传播、多次散射和通量计算等
一般粒子输运问题都是适用的。
1,屏蔽问题模型
在反应堆工程和辐射的测量与应用中, 常
常要用一些吸收材料做成屏蔽物挡住光子或中
子 。 我们所关心的是经过屏蔽后射线的强度及
其能量分布, 这就是屏蔽问题 。
当屏蔽物的形状复杂, 散射各向异性, 材
料介质不均匀,核反应截面与能量, 位置有关
时, 难以用数值方法求解, 用蒙特卡罗方法能
够得到满意的结果 。
粒子的输运问题带有明显的随机性质，粒
子的输运过程是一个随机过程。粒子的运动规
律是根据大量粒子的运动状况总结出来的，是
一种统计规律。蒙特卡罗模拟，实际上就是模
拟相当数量的粒子在介质中运动的状况，使粒
子运动的统计规律得以重现。不过，这种模拟
不是用实验方法，而是利用数值方法和技巧，
即利用随机数来实现的。
为方便起见，选用平板屏蔽模型，在厚度为 a，
长、宽无限的平板左侧放置一个强度已知，具有已知
能量、方向分布的辐射源 S 。求粒子穿透屏蔽概率
（穿透率）及其能量、方向分布。穿透率就是由源发
出的平均一个粒子穿透屏蔽的数目。
同时，假定粒子在两次碰撞之间按直线运动,且粒
子之间的相互作用可以忽略。
2,直接模拟方法
直接模拟方法就是直接从物理问题出发,
模拟粒子的真实物理过程 。
1) 状态参数与状态序列
2) 模拟运动过程
3) 记录结果
粒子在介质中的运动的状态，可用一组参数来描
述，称之为 状态参数 。它通常包括：粒子的空间位置 r,
能量 E 和运动方向 Ω，以 S＝ ( r,E,Ω ) 表示。
有时还需要其他的参数，如粒子的 时间 t 和附带
的权重 W,这时 状态参数 为 S'＝ ( r,E,Ω,t,W ) 。
状态参数 通常要根据所求问题的类型和所用的方
法来确定。
对于无限平板几何，取 S＝ ( z,E,cosα)
其中 z 为粒子的位置坐标,α为粒子的运动方向与
Z 轴的夹角。
对于球对称几何,取 S＝ ( r,E,cosθ)
其中 r 表示粒子所在位置到球心的距离,θ为粒子
的运动方向与其所在位置的径向夹角。
1) 状态参数与状态序列
粒子第 m 次碰撞后的状态参数为
或
它表示一个由源发出的粒子, 在介质中经过 m 次碰撞
后的状态, 其中
rm,粒子在第 m 次碰撞点的位置
Em,粒子第 m 次碰撞后的能量
Ωm：粒子第 m 次碰撞后的运动方向
tm,粒子到第 m 次碰撞时所经历的时间
Wm,粒子第 m 次碰撞后的权重
有时, 也可选为粒子进入第 m 次碰撞时的状态参数 。
),,( mmmm E ΩrS ?
),,,,( mmmmmm WtE ΩrS ??
一个由源发出的粒子在介质中运动, 经过若干次
碰撞后, 直到其运动历史结束 （ 如逃出系统或被吸收
等 ） 。 假定粒子在两次碰撞之间按直线运动, 其运动
方向与能量均不改变, 则粒子在介质中的运动过程可
用以下碰撞点的 状态序列 描述,
S0, S1, …, SM-1, SM
或者更详细些,用
来描述 。 这里 S0 为粒子由源出发的状态, 称为初态,
SM 为粒子的终止状态 。 M 称为粒子运动的链长 。
这样的序列称为粒子随机运动的历史, 模拟一个
粒子的运动过程, 就变成确定状态序列的问题 。
?
?
?
?
?
?
?
?
?
?
?
?
?
MM
MM
MM
EEEE
ΩΩΩΩ
rrrr
,,,,
,,,,
,,,,
110
110
110
?
?
?
为简单起见，这里以中子穿透均匀平板的模型来
说明，这时 状态参数 取 S＝ ( z,E,cosα)。
模拟的步骤如下,
(1) 确定初始状态 S0,
确定粒子的初始状态，实际上就是要从中子源的
空间位置、能量和方向分布中抽样。设源分布为
则分别从各自的分布中抽样确定初始状态。
对于平板情况,
抽样得到 z0＝ 0。
2) 模拟运动过程
)( c o s)()()c o s,,( 030201000 ?? fEfzfEzf ?
)()( 001 zzf ??
(2) 确定下一个碰撞点,
已知状态 Sm-1,要确定状态 Sm,首先要确定下一个
碰撞点的位置 zm。 在相邻两次碰撞之间, 中子的输运
长度 l 服从如下分布,
对于平板模型, l 服从分布,
其中, Σ t 为介质的中子宏观总截面,
积分 称为粒子输运的自由程数,
系统的大小通常就是用系统的自由程数表示的 。
? ?? ?????????? ?????? l mmmtmmmt ldElEllf 0 111111 ),(e x p),()( ΩrΩr
? ?? ?????????? ?????? l mmmtmmmt ldElzElzlf 0 111111 ),c o s(e x p),c o s()( ??
? ????? ???l mmmt ldEl0 111 ),( Ωr
显然, 粒子输运的自由程数服从指数分布,
因此从 f ( l ) 中抽样确定 l,就是要从积分方程
中解出 l。
对于单一介质
则 下一个碰撞点的位置
如果 zm≥ a,则中子穿透屏蔽, 若 zm≤ 0,则中子被反
射出屏蔽 。 这两种情况, 均视为中子历史终止 。
?ln),(0 111 ???????? ???l mmmt ldEl Ωr
)(
ln
1??
??
mt E
l ?
1
1
1
11
c o s)(ln
c o s
?
?
?
??
???
???
m
mt
m
mmm
Ez
lzz
??
?
(3) 确定被碰撞的原子核,
通常介质由几种原子核组成, 中子与核碰撞时,
要确定与哪一种核碰撞 。 设介质由 A,B,C 三种原子
核组成, 其核密度分别为 NA,NB,NC,则介质的宏观
总截面为,
其中 分别为核 A,B,C 的宏观总截面 。 其定
义如下,
分别表示 (·)核的宏观总截
面, 核密度和微观总截面 。
)()()()( 1111 ???? ??????? mCtmBtmAtmt EEEE
CtBtAt ???,,
)()( 1)()(1)( ????? ?? mtmt ENE ?
)()( 1)()(1)( ?????? mtmt ENE ?、、
由于中子截面表示中子与核碰撞可能性的大小,
因此, 很自然地, 中子与 A,B,C 核发生碰撞的几率
分别为,
利用离散型随机变量的抽样方法, 确定碰撞核种类,
)(
)(,
)(
)(,
)(
)(
1
1
1
1
1
1
?
?
?
?
?
?
?
??
?
??
?
??
mt
m
C
t
C
mt
m
B
t
B
mt
m
A
t
A E
EP
E
EP
E
EP
核碰撞与
核碰撞与
核碰撞与
C
BPP
AP
BA
A
??
?
?
?
＞
≤
＞
≤
(4) 确定碰撞类型,
确定了碰撞的核 (比如 B核 )后, 就要进一步确定碰
撞类型 。 中子与核的反应类型有弹性散射, 非弹性散
射, (n,2n)反应, 裂变和俘获等, 它们的微观截面分
别为
则有
各种反应发生的几率分别为
)()()()()( 111)2,(11 ????? mBcmBfmB nnmBinmBel EEEEE ????? 和、、、
)()()()()()( 111)2,(111 ?????? ????? mBcmBfmB nnmBinmBelmBt EEEEEE ??????
)()(
)()(
)()(
)()(
)()(
11
11
11)2,()2,(
11
11
??
??
??
??
??
?
?
?
?
?
m
B
tm
B
cc
m
B
tm
B
ff
m
B
tm
B
nnnn
m
B
tm
B
inin
m
B
tm
B
elel
EEP
EEP
EEP
EEP
EEP
??
??
??
??
??
利用离散型随机变量的抽样方法, 确定反应类型 。
在屏蔽问题中, 中子与核反应常只有弹性散射和
吸收两种类型, 吸收截面为,
这时, 总截面为,
发生弹性散射的几率为,
若, 则为弹性散射；否则为吸收, 发生吸收反
应意味着中子的历史终止 。
)(
)(
1
1
?
??
m
B
t
m
B
el
el E
EP
?
?
elP??
)()()( 111 ??? ?? mBcmBfmBa EEE ???
)()()( 111 ??? ?? mBamBelmBt EEE ???
(5) 确定碰撞后的能量与运动方向,
如果中子被碰撞核吸收, 则其输运历史结束 。 如
果发生弹性散射, 需要确定散射后中子的能量和运动
方向 。 中子能量 Em 为,
A是碰撞核的质量与中子质量之比, 一般就取元素的原
子量； θC 为质心系中 中子散射前后方向间的夹角, 即
偏转角 。 可从 质心系中 弹性散射角分布
fC(μC) 中抽样产生 。 实验室系散射角 θL的余弦 μL为,
? ? ? ?? ?
2
1
)
1
1(
c o s11
2
?
??
???? ?
A
Ar
rr
E
E Cmm ?
CC ?? co s?
C
C
L AA
A
?
??
21
1
2 ??
??
如果给出实验室系散射角余弦分布 fL(μL),可直接
从 fL(μL)中抽取 μL,此时能量 Em与 μL的关系式为,
确定了实验室系散射角 θL后,
再使用球面三角公式
确定 cosαm,
其中 χ 为在 ［ 0,2π ］ 上均匀分布的方位角 。
?????? c o ss i ns i nc o sc o sc o s 11 ????? ?? LmLmm
? ?22221 1)1( LLmm AAEE ?? ????? ?
至此, 由 Sm-1完全可以确定 Sm。 因此, 当中子由
源出发后, 即 S0确定后, 重复步骤 (2)～ (5),直到中子
游动历史终止 。 于是得到了一个中子的随机游动历史
S0, S1, …, SM-1, SM,即
也就是模拟了一个由源发出的中子的运动过程 。
?
?
?
?
?
?
?
?
?
?
?
?
?
MM
MM
MM
EEEE
zzzz
???? c o s,c o s,,c o s,c o s
,,,,
,,,,
110
110
110
?
?
?
以上模拟过程可分为两大步：第一步确定粒子的
初始状态 S0,第二步由状态 Sm-1来确定 状态 Sm。 这第二
步又分为两个过程：第一个过程是确定碰撞点位置 zm,
称为输运过程；第二个过程是确定碰撞后粒子的能量
及运动方向, 称为碰撞过程 。 对于中子而言, 碰撞过
程是先确定散射角, 进而确定能量和运动方向；而对
于光子, 碰撞过程是先确定能量, 再确定散射角以及
运动方向 。 重复这两个过程, 直至粒子的历史终止 。
这种模拟过程, 是解任何类型的粒子输运问题所
共有的, 它是蒙特卡罗方法解题的基本手段 。
在获得中子的随机游动历史后，我们要对所要计
算的物理量进行估计。对于屏蔽问题，我们要计算中
子的穿透率。考察每个中子的随机游动历史，它可能
穿透屏蔽（ zM≥ a），可能被屏蔽发射回来（ zM≤ 0），
或者被吸收。设第 n 个中子对穿透的贡献为 ηn,则
如果我们共跟踪了 N 个中子，则穿透屏蔽的中子数为,
3) 记录结果
?
?
?
?
??
，或者被吸收当
当
0,0
,1
M
M
n z
az?
?
?
?
N
n
nN
1
1 ?
则穿透屏蔽概率的近似值为,
它是穿透率的一个无偏估计 。
我们称这种直观地模拟过程和估计方法为直接模
拟方法 。 在置信水平 1－ α ＝ 0.95 时, 的误差为,
其中 为 ηn的均方差, 由于 ηn是一个服从二项分布的随
机变量, 所以
或
?
?
??
N
n
nN NN
NP
1
1)1( 1? ?
)1(?NP
NNPP N
??? ??? 96.1? )1( ???
??
)?1(??
)1(
)1()1(2
2
NN PP
PP
??
??
?
?
?
?
为得到中子穿透屏蔽的能量, 角分布, 将能量,
角度范围分成若干个间隔,
其中 Emax,Emin分别表示能量的上, 下限, 对于穿透屏
蔽的中子按其能量, 方向分间隔记录 。 设一穿透屏蔽
的中子能量为 EM,其运动方向与 Z轴夹角为 αM,若能
量 EM属于第 i 个能量间隔 Δ Ei,角度 αM属于第 j 个角
度间隔 Δ αj,则分别在第 i 个能量计数器及第 j 个角
度计数器中加 "1"。
20 10
m a x01m i n
???? ?????
?????
J
I EEEEE
?
?
跟踪 N 个中子后, 则
分别为穿透中子的能量分布和角分布 。 其中 N1,i 和
N2,i 分别为第 i 个能量和第 j 个角度间隔的穿透中
子数 。 归一后分别为,
Jj
N
N
P
Ii
EN
N
P
j
j
N
i
i
N
j
i
,,2,1?
,,2,1?
,2)1(
,1)1(
,2
,1
?
?
?
??
?
?
??
?
?
Jj
N
N
P
P
P
Ii
EN
N
P
P
P
j
j
N
N
N
i
i
N
N
N
j
j
i
i
,,2,1
?
?
?
,,2,1
?
?
?
1
,2
)*1(
)1(
)*1(
1
,1
)*1(
)1(
)*1(
,2
,2
,1
,1
?
?
?
??
??
?
??
??
?
3,简单加权法
从模拟物理过程来说, 直接模拟法是最简单, 也
是最基本的方法 。 但是, 在直接模拟法中, 不管中子
在屏蔽中经过多少次碰撞, 只要在介质中被吸收, 对
穿透的贡献就为零；因此在所跟踪的粒子中绝大部分
都对穿透没有贡献 。 而在许多屏蔽问题中, 穿透率的
数量级在 10-6到 10-8。 进一步, 如果我们要求穿透率达
相对误差小于 1％, 即
那么, N 要大到惊人的数量级 1010到 1012。 显然, 这时
用直接模拟法计算不是很有效 。
%1)1(|?| ???????? PNPN PPPNP PP ??? ?????
屏蔽物一般是由吸收强的介质组成，因此在每次
碰撞时，粒子很有可能被吸收而停止跟踪。现在改变
模拟方法，在判断碰撞类型 时，可以认为粒子
的 部分是弹性散射，而其余部
分被吸收，即人为地把中子分成两部分，一部分弹性
散射，一部分吸收。弹性散射这部分继续跟踪；吸收
部分则停止跟踪。也就是说，我们利用中子权重的变
化来反应继续弹性散射的部分。这就是简单加权法的
基本思想。
1) 简单加权法
elP??
)()( 11 ?? ??? mBtmBelel EEP
显然, 在加权法中中子的权重 W 已成为中子状态
参数的组成部分 。 这时, 中子历史成为,
对源中子, 取 W0=1。 经过碰撞中子权重的变化为,
因子 称为尚存因子 。
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
,,,,,
c o s,c o s,,c o s,c o s
,,,,
,,,,
110
110
110
110
MM
MM
MM
MM
WWWW
EEEE
zzzz
?
?
?
?
????
)(
)(
1
1
1
?
?
? ?
??
m
B
t
m
B
el
mm E
EWW
)()( 11 ?? ?? mBtmBel EE
这时, 第 n 个中子对穿透的贡献为,
如果我们共跟踪了 N个中子, 则穿透率 P的无偏估计为,
类似地, 可以得到 穿透中子的能量分布和角分布 。 只
不过在对各计数器进行的加 "1" 操作改为加 WM。
??
?
?
??
0,0
,
M
MM
n z
azW
当
当?
?
?
?
N
n
nN NP
1
)2( 1? ?
简单加权法的方差估计为,
与直接模拟法相比，有
注意到 ζn≤1，有
这表明简单加权法的方差小于直接模拟法的方差。这
是因为加权法比直接模拟法减少了一次随机抽样。
2) 简单加权法的方差
2)2(
1
22 )?(1?
N
N
n
n PN ?? ?
?
?? ?
?
?
???
N
n
nnN
1
222 )(1?? ????
??
22 ?? ?? ?? ?
加权法的思想在蒙特卡罗方法中用途很广泛。例
如，对于具有中子增殖反应，如裂变,(n,2n),(n,3n)
反应的中子输运问题，一个中子与核发生碰撞后，根
据反应的类型会产生不同数量的次级中子，每个次级
中子又会产生新的次级中子，这样链锁反应 下去，使
得用直接模拟法模拟每一个中子是非常困难的。这种
情况可以利用加权法来处理。
3) 权重方法的其它应用
中子与核发生碰撞 后，产生的次级中子平均数为,
这里 νf 为裂变次级中子数 。 于是, 碰撞后的权重为,
而决定碰撞类型的几率分别为,
其中
加权法的思想, 还可以应用到连续分布情况和偏
倚抽样的问题
t
ffnnnninel
?
??????? ?????????? )3,()2,( 3211
??? ?1mm WW
t
ff
t
nn
t
nn
t
in
t
el
?
??
?
?
?
?
?
?
?
? ???,3,2,,)3,()2,(
ffnnnninelt ??????? ?????????? )3,()2,( 3211
4,统计估计法
加权法虽然改进了直接模拟法, 但它同样只关心
中子是否穿透屏蔽这一信息, 因此对每一个中子历史
的信息利用得很不充分 。 统计估计法能够较多地利用
中子的历史信息, 因而能得到更好的结果 。
一个中子, 可能在介质内不发生碰撞而直接穿透
屏蔽, 也可能在介质内发生一次碰撞后再穿透屏蔽,
或经过二次碰撞穿
透屏蔽, 等等,
这些事件是互不相
容的, 因此穿透概
率 P 可表示为,
其中 Pm 是中子恰好经过 m 次碰撞而穿透屏蔽的概率 。
这表明, 可以用求 Pm (m=0,1,… ) 的方法得到 P。 这样,
中子对穿透概率的贡献就不只限于末次碰撞了 。
??
?
?
0m
mPP
α0 α1 S0
S1
Sm αm
P1
P0
Pm
Z a 0
设中子的历史为,
根据该中子的历史, 我们可以估计出中子恰好经
过 m 次碰撞后, 穿透屏蔽的部分
显然, 具有初态 S0＝ ( 0,E0,cosα0,W0 ) 的中子, 未
经碰撞直接穿透的部分是,
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
,,,,,
c o s,c o s,,c o s,c o s
,,,,
,,,,
110
110
110
110
MM
MM
MM
MM
WWWW
EEEE
zzzz
?
?
?
?
????
1,,1,0,? ?? MmP m ?
?
?
?
?
?
? ???
0
000 c o s)(e x p
?
?
aEWP
t
类似地, 在经过了第 m 次碰撞后的中子具有状态
Sm＝ ( zm,Em,cosαm,Wm ), 其可能穿透的部分, 正好是
一个中子恰好经过 m 次碰撞穿透的部分,
这里的这种估计技巧, 由于是对每次碰撞后的状
态, 求其后未经碰撞直接穿透的贡献, 因此该方法也
称为最后自由飞行估计 。
1,,1,0
0
0c o s
c o s
)(e x p?
??
?
?
?
?
?
?
?
?
?
?
?
? ?
??
?
Mm
za
EW
P m
m
m
mtm
m
?
其它，
,?
?
于是得到该中子对穿透的贡献,
如果我们共跟踪了 N个中子, 则穿透率 P的估计为,
其方差估计为,
??
?
?
1
0
?? M
m
mPP
?
?
?
N
n
N nPNP
1
)3( )(?1?
? ? ? ?2)3(
1
22 ?)(?1)?(
N
N
n
PnPNP ?? ?
?
?
在直接模拟方法中，相对误差为
其中 为与置信水平 1－ α相应的量。
如果构造一个新的概率模型，使得该模型的穿透
率 P*与原模型的穿透率 P之间存在关系,
使用直接模拟方法,相对误差为
5,指数变换法
PNPN
PP
PN ??
??? ??? ????? )1(
??
PKP ??*
*
*
PN ?
? ???
如果令 ε*＝ ε,即
这意味着, 达到同样的相对误差, 跟踪粒子的数
目缩小 K 倍, 从而减少 K 倍的计算量 。 指数变换法就
是构造一个新的概率模型的一个有效方法 。
K
N
N
PKNPNPN
?
??
?
?
?
?
*
***
??? ???
构造如下伪过程：宏观总截面为
散射截面仍为 Σel(E)。 其中 Emin,Emax 分别为能量的下
限和上限, α为粒子的运动方向与 Z 轴的夹角 。 可以证
明这个伪过程的穿透概率 P* 与原过程的穿透概率 P之
间有如下关系,
显然, 。 因伪过程与原过程的结果相差 e 指
数, 所以该方法称为指数变换法 。
)(m i n
c o s)()(
m a xm i n
*
EC
CEE
tEEE
tt
??
?????
??
?
PeP aC ?? ?*
1?? ?aCeK
分析一下伪过程的定义,可以明显看出 P*增大的原
因 。 当 cosα=1 时, 粒子运动方向与 Z 轴方向一致, 其
截面最小, 粒子沿 Z 轴方向输运的距离较远；而当
cosα=－ 1 时, 粒子运动方向背向 Z 轴方向, 这时其截
面最大, 粒子向后输运的距离较短 。 因此, 截面变换
的结果是加强了粒子向前运动的能力, 因而使穿透概
率增大 。
伪过程的构造与几何形状及所考虑的问题有关 。
比如, 对球形几何, 使用指数变换法求穿透概率时,所
构造的宏观总截面与平板屏蔽的情况不同, 粒子的模
拟方法也较复杂 。
6,蒙特卡罗方法的效率
衡量一种蒙特卡罗技巧的好坏, 除了看其方差大
小外, 还要看其所需费用 （ 计算时间 ） 多少, 即从该
技巧的效率 Ef（ 方差与费用乘积的倒数 ） 全面考虑,
其中 σ2 为方差, T 为所需费用 。 Ef 大时, 所用方
法的效率高；否则, 效率低 。
在一般情况下, 有些方法虽然减小了方差, 却增
加了费用 。 例如, 加权法, 统计估计法虽然较直接模
拟方法减小了方差, 却使每个粒子的运动链长增加,
或记录贡献的计算时间增加 。 因 此, 不能认为方差小
的方法一定好, 要从方法的效率全面考虑 。 在有些情
况下, 直接模拟方法仍然是一个被广泛使用的方法 。
TE f ?? 2
1
?
? 作 业
1) 光子散射后能量分布的抽样
把光子散射能量分布改写成如下形式进行抽样,
??
?
?
?
??
?
?
?
?
?
??
?
? ??
?
?
?
?
?
?
?
?
??
?
??
?
? ???
22
2 11111
)(
1)(
xxxx
x
K
xf
?
?
?
?
在 ［ 1,1+2α］ 上定义如下函数,
3
2
2
2
1
2
21
)1(
)(
2
)(
1
1
)21)((
2
)(
2
1
)(
1
2
21
)(
x
x
K
xH
x
K
xH
xf
x
xf
?
?
?
?
?
?
?
?
?
?
??
?
?
?
?
? ??
?
?
?
?
?
?
?
?
?
?
?
??
?
?
?
?
?
?
?
? 作 业
1) 给出分布密度函数
的抽样方法。
＞
≤
