深貫穿問題脈沖高度計數(shù)減方差技巧

李長楷 湯曉斌 岳愛忠

1（南京航空航天大學 核科學與工程系 南京 210016）

2（中石油測井有限公司 西安 710077）

深貫穿問題脈沖高度計數(shù)減方差技巧

李長楷1湯曉斌1岳愛忠2

1（南京航空航天大學 核科學與工程系 南京 210016）

2（中石油測井有限公司 西安 710077）

在處理核測井、輻射防護以及核醫(yī)學等深貫穿問題時簡單使用無偏方法會導致計算時間長且誤差大。較新版本的MCNP程序在進行脈沖高度計數(shù)時可使用多種減方差技巧。本文使用蒙特卡羅程序代碼MCNP5 1.51對包括Mesh-based權(quán)窗在內(nèi)的多種減方差技巧的計算效率進行測試研究，旨在尋求較好的脈沖高度減方差技巧(Pulse height variance reduction techniques, PHVRT)來解決深貫穿問題中計數(shù)效率低的問題。計算表明使用DX球配合強迫碰撞生成的權(quán)窗可大幅提高計算效率，對各參數(shù)進行優(yōu)化后又可進一步大幅度提高計算效率。最終得出使用Mesh-based權(quán)窗配合其他減方差技巧可以很好解決深貫穿問題計算效率低問題。

深貫穿，脈沖高度計數(shù)，減方差，Mesh-based權(quán)窗

蒙特卡羅方法被公認為處理放射性輸運問題最有效的方法[1]，其本質(zhì)是通過模擬大量粒子輸運過程來推斷粒子平均輸運過程，這使其在計算時很耗費時間[2?6]，MCNP記數(shù)結(jié)果可以表示為：

式中，σ是計算值和平均值的偏差；T是計算時間，抽樣較好時和抽樣粒子數(shù)成反比，而T和抽樣粒子數(shù)成正比，因此對于較好抽樣的問題FOM應該趨于穩(wěn)定。

MCNP脈沖高度計數(shù)是粒子輸運過程中進入柵元和離開柵元時的能量變化，包括源粒子以及所有次級粒子在柵元的能量沉積，該能量沉積就是粒子輸運過程中在柵元的能量脈沖，脈沖高度計數(shù)是MCNP唯一需要在粒子輸運過程完成才能實現(xiàn)的計數(shù)[7?8]。常規(guī)的玻爾茲曼通量計數(shù)可以在粒子離開柵元或者在柵元里發(fā)生碰撞時使用減方差技巧，而在對脈沖高度計數(shù)使用減方差技巧時需要額外儲存整個粒子輸運過程，這使得脈沖高度計數(shù)減方差更為復雜[9?12]。

前人已在關(guān)于深貫穿問題減方差技巧方面做了大量優(yōu)秀的工作[13?17]，但很少涉及關(guān)于這類問題具體解決辦法，以及各脈沖高度減方差技巧(Pulse height variance reduction techniques, PHVRT)實際效果。本文對各減方差技巧的計算效率進行了詳細的研究。

1 研究方法

1.1 可用于脈沖高度計數(shù)的減方差技巧

1.1.1 源偏倚

源偏倚可以通過降低權(quán)重來抽樣較多感興趣的源粒子，包括能量、時間、角度偏倚等。使用該偏倚會按照偏倚的分布進行抽樣，而不是按照源的真實分布來抽樣，而抽樣粒子的權(quán)重將被修改為真實分布概率除以偏倚概率分布的值，因此能保證總權(quán)重守恒和結(jié)果無偏。

1.1.2 幾何分裂與輪盤賭

通過給每個柵元賦一個重要性值來控制粒子在柵元間輸運的行為，以增加感興趣柵元的抽樣，同時減少不感興趣柵元的抽樣。如前一柵元重要性為i1，第二柵元的重要性為i2，如i2＞i1，則分裂為i2/i1個粒子數(shù)，每個粒子權(quán)重為Wgt×(i1/i2)，當i2＜i1，則進行輪盤賭，粒子存活的幾率為i1/i2，存活粒子權(quán)重為Wgt×(i2/i1)，通過設置合適的權(quán)重該方法也可增加感興趣區(qū)域的抽樣。

1.1.3 暗含俘獲

暗含俘獲也稱為存活偏倚，即碰撞不吸收（殺死）粒子而是降低其權(quán)重，如碰撞總截面為X，吸收截面為X1，則碰撞后粒子權(quán)重降為Wgt(1?X1/X)，而不是按X1/X的概率被吸收，從而保證了較多的有效抽樣粒子，但這樣做的缺點是保留了大量權(quán)重極低的粒子，因此該方法必須與權(quán)重截斷或權(quán)窗一起使用來控制低權(quán)重粒子的數(shù)量。

1.1.4 指數(shù)變換

通過修改粒子在物質(zhì)中的微觀截面來改變指向特定方向的自由程，同時粒子權(quán)重也會被相應修改。

式中，φ*是修改后的總截面；φ是修改前的總截面；μ是粒子方向和拉伸方向夾角的余弦。P＞0時，會使其在物質(zhì)中的穿透性增強，該方法與強迫碰撞不可發(fā)生在同一柵元里，本文下面所用指數(shù)變換和強迫碰撞結(jié)合使用是指二者在不同柵元里使用。

1.1.5 強迫碰撞

強迫碰撞與暗含俘獲相似，將進入柵元的粒子分為碰撞和不碰撞兩類，并由不碰撞的粒子攜帶剩下的權(quán)重離開柵元。該方法不僅適用于增加光疏介質(zhì)柵元碰撞，也同時適用于增加光密介質(zhì)的通過性。由于在使用強迫碰撞時權(quán)重截斷是關(guān)閉的，因此在相鄰幾個柵元同時使用強迫碰撞會導致大量低權(quán)重粒子產(chǎn)生并極大增加計算時間。

1.1.6 DX球(DXTRAN)

DX球一般用于不能很好抽樣的小區(qū)域，通過設置DX球，在球外每次碰撞中都會產(chǎn)生一個DXTRAN粒子(DXTRAN particle)，該粒子按確定幾率散射到DX球表面，期間不再參與碰撞，權(quán)重按確定性輸運計算。由于強迫碰撞會增加球外碰撞次數(shù)，DX球粒子數(shù)也會相應增加，因此在DX球外柵元使用強迫碰撞可以增強DX球效果。

1.1.7 權(quán)窗

權(quán)窗原理是在源與計數(shù)區(qū)之間設立一系列權(quán)重區(qū)間，權(quán)重窗上限值取下限值的固定倍數(shù)，粒子權(quán)重高于能窗上限時會分裂，低于能窗下限時粒子會進行輪盤賭。該方法與幾何分裂相似，權(quán)窗與幾何分裂的區(qū)別在于幾何分裂不考慮權(quán)重值的大小只考慮相鄰柵元間的重要性比值，而權(quán)窗中需要將粒子的權(quán)重值與權(quán)窗上、下限進行比較來控制其分裂和輪盤賭，其特有的抑制過度偏倚功能使其和其他減方差技巧有很好的兼容性。可以將權(quán)窗形象比喻為“由源通向計數(shù)區(qū)域的走廊”。

MCNP提供的權(quán)窗生成器能自動計算權(quán)窗重要性函數(shù)，通過模擬一定數(shù)量的粒子，重要性的計算可表示為：重要性=（進入該柵元粒子對計數(shù)的貢獻/進入柵元的總權(quán)重），然后按重要性倒數(shù)的分布指定權(quán)窗下限值。權(quán)窗的生成過程也涉及到粒子統(tǒng)計過程，通過使用偏倚技巧可較快生成有效的權(quán)窗。權(quán)窗生成器有效避免了人為因素，且其生成的權(quán)窗效果往往好于人為設定的重要性值。

1.2 計算模型和計算方法

本文計算模型是γ密度測井(圖1)。直徑10 cm的鎢鎳合金測井工具放置于直徑20 cm的井眼內(nèi)，井內(nèi)充滿水，各向同性的137Cs點源位于一個45°斜向上的準直裝置中，遠探測器為?50×48 mm NaI(TI)探測器，位于距源40 cm處，晶體表面鍍0.03 cm的鎘膜。由于鎢鎳合金對γ射線屏蔽系數(shù)很高，對探測器計數(shù)貢獻最大的是經(jīng)斜準值孔射入地層的粒子，這部分粒子占總粒子份額很少，這會引起計數(shù)區(qū)域抽樣不足，該問題屬于深貫穿問題，有必要使用減方差技巧。

圖1 γ密度測井的蒙特卡羅模型Fig.1 Monte Carlo model of gamma-ray density well.

由于減方差技巧的偏倚性，不能同時使遠近探測器獲得較好的抽樣，本文只對遠探測器計數(shù)進行討論。本文各減方差技巧的參數(shù)選擇均是使用MCNP_pstudy功能，在參數(shù)可選區(qū)間做均勻線性插值分別計算并對比結(jié)果得到的最優(yōu)值，MCNP插值命令是c @@@ p=uniform 10 min max。對0?0.665MeV能量區(qū)間設100個能量箱，對探測器進行脈沖高度計數(shù)并對計數(shù)結(jié)果進行高斯展寬。由于脈沖高度計數(shù)不可以和權(quán)窗生成器WWG同時使用，本文使用F4（徑跡長度計數(shù)）生成權(quán)窗，并將該權(quán)窗用于脈沖高度計數(shù)。

圖2為加Mesh權(quán)窗效果圖，本文選用直角坐標Mesh，虛擬網(wǎng)格數(shù)量為20×20×30，采用均勻分布。

圖2 權(quán)窗虛擬網(wǎng)格軸向和徑向分布(a) xz方向，(b) xy方向Fig.2 Radial and axial cross-section views of weight windowvirtual grid distribution.(a) Axial cross-section views, (b) Radial cross-section views

2 計算結(jié)果與分析

2.1 不使用任何減方差技巧

圖3 不使用減方差技巧的探測器響應曲線Fig.3 Corresponding curve got with no variance reduction techniques.

由圖3可見，計算結(jié)果各能量箱誤差太大，導致計算結(jié)果不可信。

2.2 使用常規(guī)減方差技巧

運行2×109粒子周期各減方差技巧的計算效率值如表1所示。本文的相對誤差是指總計數(shù)的相對誤差。

不使用任何減方差技巧運行2×109個粒子周期，相應探測器響應曲線如圖3所示。

表1 單個減方差技巧計算效率Table 1 Computational efficiency with one PHVRT.

由表1，使用DX球計算效率最高，因為該方差技巧可以極大程度提高指定區(qū)域的抽樣，指數(shù)變換、源偏倚和幾何分裂也可以通過提高在指定方向的抽樣從而較大幅度提高計數(shù)效率，暗含俘獲由于其全局性特點，在提高對計數(shù)貢獻大的區(qū)域粒子生存幾率的同時也提高了其他區(qū)域的粒子生存幾率，單獨使用時其對計數(shù)效率影響不大。強迫碰撞也可以提高粒子生存幾率但其可以發(fā)生在指定區(qū)域，因此在對計數(shù)貢獻大的區(qū)域使用強迫碰撞可以提高計數(shù)效率。

將DX球和其他減方差技巧結(jié)合使用運行2×109粒子周期，計算效率如表2所示。

表2 DX球與其他減方差技巧結(jié)合計算效率Table 2 Computational efficiency of DXTRAN used in conjunction with the other variance reduction techniques.

由表2可見，DX球和其他所有減方差技巧結(jié)合使用都會較大幅度提高其計算效率，其中DX球和源偏倚同時使用有較高的計算效率，因為該問題中幾乎只有通過準直孔進入地層的粒子才可能對計數(shù)區(qū)有貢獻，其對應響應曲線如圖4所示。

圖4 DX球和源偏倚配合使用計算結(jié)果Fig.4 Corresponding response curve achieved by using DXTRAN with source bias.

圖4 中其各能量箱計算誤差約為20%，一些區(qū)間誤差甚至更大，結(jié)果依然不可信?？梢姵Ｒ?guī)減方差技巧不能很好地解決深貫穿問題。下面的計算會采用一種自動減方差技巧-權(quán)窗減方差。

2.3 使用Mesh權(quán)窗

權(quán)窗值只與計數(shù)相關(guān)而和柵元密度無關(guān)，為較快生成有效權(quán)窗本文對主要材料做降低密度處理，將地層和鎢鎳合金的密度分別降為實際值的1/10。為驗證這樣做的可行性，分別用降低后的密度和實際密度條件下生成的權(quán)窗用于實際密度的計算，并對計算結(jié)果進行對比。如圖5所示，兩條響應曲線吻合度很高，其中兩條曲線各能量箱誤差都在0.5%以下，可見密度對權(quán)窗有效性沒有影響。

圖5 減密度和真實密度所得權(quán)窗結(jié)果對比Fig.5 Corresponding response curve achieved by using weight window generated by reduced density and real density respectively.

表3 Mesh權(quán)窗與其他一種減方差技巧結(jié)合使用計算效率Table 3 Computational efficiency achieved by using mesh-based weight window in conjunction with one of the other PHVRT.

使用Mesh-based權(quán)窗，并使用一種其他減方差技巧，在1/10主要物質(zhì)密度條件下運行2×107粒子周期生成權(quán)窗，將該權(quán)窗用于接下來的真實密度條件計算，相應權(quán)窗計算效率值見表3。

其中迭代1表示使用降低密度條件下生成的權(quán)窗進行實際密度計算，計算中使用相應的減方差技巧；迭代2表示使用降低密度條件下生成的權(quán)窗進行實際密度計算，計算中不使用任何其他減方差技巧；迭代3表示使用迭代1過程中生成的權(quán)窗進行實際密度計算，計算中和權(quán)窗生成時使用了同樣的減方差技巧；—表示無法得到有效計算結(jié)果。

由表3可見，權(quán)窗和暗含俘獲有較好的相容性，暗含俘獲保證抽樣粒子有較大數(shù)量，其和權(quán)窗配合使用可以去除大量對計數(shù)貢獻小的低權(quán)重粒子，而使計算時間不致太長。定向計數(shù)功能的DX球生成的Mesh-based權(quán)窗也有較好的效果。

DX球結(jié)合其他減方差技巧一起使用生成權(quán)窗，只使用除DX球外的一種減方差技巧。

由表4，通過使用強迫碰撞配合DX球可以在迭代3中大幅提高Mesh權(quán)窗的計算效率，其原因是強迫碰撞可以增加對計數(shù)區(qū)有貢獻柵元的粒子存活幾率，也就增加了在該區(qū)域碰撞產(chǎn)生DXTRAN粒子的數(shù)量，并通過權(quán)窗去除大量低權(quán)重粒子，從而有較大的計算效率。表5為加一種其他減方差技巧和二者配合使用的計算效率。

表4 Mesh權(quán)窗和DX球與其他減方差技巧結(jié)合計算效率Table 4 Computational efficiency of Mesh-based weight window and DXTRAN in conjunction with one of the other PHVRT.

表5 Mesh權(quán)窗和強迫碰撞以及DX球與其他減方差技巧結(jié)合使用計算效率Table 5 Computational efficiency of Mesh-based weight window, DXTRAN and forced collision in conjunction with one of the other PHVRT.

生成權(quán)窗過程中再加上其他減方差技巧反而會降低所生成權(quán)窗的計算效率。其原因是多種減方差計數(shù)互相干擾使生成權(quán)窗下限過低導致大量極低權(quán)重粒子產(chǎn)生而影響計算速度；生成權(quán)窗上限過高使大量粒子因輪盤賭而損失，導致抽樣不足。

目前為止，使用DX球和強迫碰撞配合權(quán)窗在迭代1中生成的權(quán)窗是最好的，為進一步提高計算效率，接下來使用該權(quán)窗配合DX球以及強迫碰撞計算時加入其他減方差技巧，各減方差組合計算效率如表6所示。

表6 Mesh權(quán)窗和強迫碰撞生成的權(quán)窗配合其他減方差技巧一起使用計算效率Table 6 Computational efficiency of improved PHVRT.

如表6所示，加上其他減方差技巧后各組合均較原來有較大提高，說明權(quán)窗和其他減方差技巧有很好的兼容性。

為顯示本文減方差技巧的效果，對比沒使用減方差技巧的二維可視化輸運過程和使用本文減方差技巧的二維可視化輸運過程。

圖6 可視化運行粒子碰撞軌跡，虛線是DX球(a) 不使用減方差技巧運行50 000個粒子周期，(b) 使用減方差技巧運行500個粒子周期Fig.6 Particle collisions in the visualization tool, the dotted box is DXTRAN.(a) Collisions of 50 000 particles histories with no PHVRT, (b) Collisions of 500 particles histories with PHVRT

圖6 (a)中整個輸運過程粒子權(quán)重不變，計數(shù)區(qū)不能得到很好的抽樣。圖6(b)粒子權(quán)重由下到上依次遞減，計數(shù)區(qū)抽樣粒子權(quán)重較低但可以得到很好的抽樣。

3 結(jié)語

本文列出多種可用于脈沖高度計數(shù)的減方差技巧并作簡要介紹說明，通過對比深貫穿問題中脈沖高度計數(shù)多種減方差技巧的計算效率，得出

Mesh-based權(quán)窗與其他減方差技巧結(jié)合使用解決這類問題有很好的效果；提出使用常規(guī)的玻爾茲曼計數(shù)產(chǎn)生的權(quán)窗用于脈沖高度技術(shù)以克服脈沖高度計數(shù)不能和權(quán)窗生成器結(jié)合使用的不足；提出權(quán)窗值和密度無關(guān)，通過降低密度可以較快生成有效權(quán)窗。

本文所有減方差技巧均適用于脈沖高度計數(shù)以外的常規(guī)玻爾茲曼計數(shù)。

1 Brown F B, Sweezy J E, Bull J S, et al. Verification of MCNP5[R]. Version 1.50, LA-UR-08-3443, USA: Los Alamos National Laboratory, 2008

2 Wagner J C, Peplow D E, Evans T M. Automated variance reduction applied to nuclear well-logging problems[J]. Nuclear Technology, 2009, 168(3): 799?809

3 Wagner J C. An automated deterministic variance reduction generator for Monte Carlo shielding applications[C]. Proceedings of the American Nuclear Society 12thBiennial RPSD Topical Meeting, 2002: 14?18

4 Junli L, Zhi Z, Rui Q, et al. An auto-importance sampling method for deep penetration problems[J]. Progress in Nuclear Science and Technology, 2011, 2: 732–737

5 Booth T E. Sample problem for variance reduction in MCNP[R]. Los Alamos National Laboratory, USA: NM, 1985

6 Cramer S N. Variance reduction methods applied to deep-penetration problems[R]. Oak Ridge National Laboratory, USA: TN, 1984

7 Hughes H G. Variance reduction with pulse-height tallies in MCNP[R]. Rapport Technique LA-UR-07-6660, USA: Los Alamos National Laboratory, 2007

8 Booth T E. Pulse-height tally variance reduction in MCNP[R]. Los Alamos Report LA-13955, USA: Los Alamos National Laboratory, 2002

9 Sood A, Forster R A, Adams B J, et al. Verification of the pulse height tally in MCNP 5[J]. Nuclear Instruments and Methods in Physics Research Section B: Beam Interactions with Materials and Atoms, 2004, 213: 167?171

10 Bull J S, Booth T E, Sood A. Implementation of pulse height tally variance reduction in MCNP5[R]. Technical Report LA-UR-08-2330, USA: Los Alamos National Laboratory, 2008

11 Booth T E. Monte Carlo variance reduction approaches for non-Boltzmann tallies[R]. Los Alamos National Laboratory, USA: NM, Funding Organization: USDOE, Washington DC, 1992

12 Olsher R H. A practical look at Monte Carlo variance reduction methods in radiation shielding[J]. Nuclear Engineering and Technology, 2006, 38: 225?230

13 Haghighat A, Wagner J C. Monte Carlo variance reduction with deterministic importance functions[J]. Progress in Nuclear Energy, 2003, 42(1): 25?53

14 Smith H P, Wagner J C. A case study in manual and automated Monte Carlo variance reduction with a deep penetration reactor shielding problem[J]. Nuclear Science and Engineering, 2005, 149(1): 23?37

15 Hendricks J S, McKinney G W. Pulse-height tallies with variance reduction[C]. Proceedings of the ANS Monte Carlo 2005 Topical Meeting, 2005: 17?21

16 van Riper K A, Urbatsch T J, Soran P D. AVATAR-automatic variance reduction in Monte Carlo calculations[R]. Los Alamos National Laboratory, USA: NM, 1997

17 Wagner J C, Haghighat A. Automated variance reduction of Monte Carlo shielding calculations using the discrete ordinates adjoint function[J]. Nuclear Science and Engineering, 1998, 128(2): 186?208

CLC TL99

Pulse-height tally variance reduction in deep penetration problem

LI Changkai1TANG Xiaobin1YUE Aizhong2

1(Department of Nuclear Science and Engineering, Nanjing University of Aeronautics and Astronautics, Nanjing 210016, China)
2(China Petroleum Logging Co., Ltd, Xi'an 710077, China)

Background: In deep penetration application such as well-logging, medical physics and radiation shielding, analog Monte Carlo method will cause long calculation time and large variance. The newly released Monte Carlo code MCNP5 version 1.51 allows almost all of MCNP's variance reduction techniques to be used with pulse height tally. Purpose: This paper aims to find an optimize pulse height variance reduction technique (PHVRT) to solve low efficiency with pulse-height tallies in deep penetration problems. Methods: A variety of variance reduction techniques including an auto-important sampling method-Mesh-based weight windows were proposed, and all kinds of variance reduction techniques proposed were tested and results with different variance reduction techniques were compared. Results: Mesh-based weight windows used with DXTRAN and forced collision get good calculation efficiency. After optimizing, the efficiency can be improved with large amplitude. Conclusion: Mesh-based weight window used with other variance reduction techniques can resolve deep penetration problem effectively.

Deep penetration, Pulse height tally, Variance reduction, Mesh-based weight windows

TL99

10.11889/j.0253-3219.2015.hjs.38.030501

No.11475087)資助

李長楷，男，1987年出生，2014年于南京航空航天大學獲碩士學位，從事核測井研究

2014-11-04，

2014-12-12