全局減方差方法在大空間γ輻射場計(jì)算中的應(yīng)用

劉利 左應(yīng)紅 牛勝利 朱金輝 商鵬 王學(xué)棟

（西北核技術(shù)研究所 西安 710024）

γ射線在近地面大空間長距離輸運(yùn)中會(huì)與空氣和地面土壤發(fā)生多次散射，使得輸運(yùn)至距輻射源公里及以上范圍內(nèi)的γ粒子數(shù)較少，致使該處的蒙特卡羅模擬統(tǒng)計(jì)漲落較大。左應(yīng)紅等［1］建立了柵元權(quán)窗結(jié)合密度逼近迭代的局域減方差方法（Local Variance Reduction，LVR）［2-3］，能夠求解離源5 km的觀測點(diǎn)處的γ射線注量。而在核設(shè)施廠房輻射環(huán)境［4］、早期核輻射環(huán)境［5］、核電磁脈沖環(huán)境［5］、人員輻射劑量學(xué)等方面的研究中，常常需要計(jì)算大空間中γ輻射場的整體分布。傳統(tǒng)的LVR方法通常只針對特定目標(biāo)探測器，難以給出理想的全空間γ輻射場分布。

相比之下，全局減方差方法（Global Variance Reduction，GVR）［6-9］更適于求解輻射場的整體分布。GVR方法通過全局權(quán)窗來控制模擬粒子的空間分布，使得模擬結(jié)果的相對偏差隨空間的分布盡可能減小，從而整體提高全局空間內(nèi)的計(jì)算效率。Davis和Andrew等［6-7］提出的基于蒙特卡羅正向計(jì)算的GVR方法，通過蒙特卡羅模擬計(jì)算得到粒子通量、數(shù)量和權(quán)重等粒子分布信息，用于產(chǎn)生基于柵元/網(wǎng)格的全局權(quán)窗參數(shù)，利用全局權(quán)窗參數(shù)開展蒙特卡羅模擬，循環(huán)迭代直至結(jié)果收斂。該方法思路清晰，適用性強(qiáng)，已廣泛用于反應(yīng)堆裝置［10-12］、托卡馬克裝置［13-14］和厚屏蔽體［15-16］等復(fù)雜模型的全局輻射環(huán)境計(jì)算。

大空間γ輻射場模擬具有空間尺度跨度大和計(jì)數(shù)區(qū)外散射較強(qiáng)等特點(diǎn)。直接應(yīng)用已有的GVR方法面臨計(jì)數(shù)柵元/網(wǎng)格體積差異和非計(jì)數(shù)區(qū)引起的“過度分裂”問題。一個(gè)是輻射場空間尺度跨度大，計(jì)數(shù)柵元/網(wǎng)格體積相差較大。而現(xiàn)有全局減方差方法更適用于等體積的柵元/網(wǎng)格，其引導(dǎo)粒子輸運(yùn)使得模擬粒子在全空間范圍內(nèi)分布較為均勻，相對誤差在全空間范圍內(nèi)相對均勻。當(dāng)輻射場內(nèi)計(jì)數(shù)柵元/網(wǎng)格體積差異較大時(shí)，體積大的柵元/網(wǎng)格內(nèi)模擬粒子數(shù)較多，體積小的柵元/網(wǎng)格內(nèi)模擬粒子數(shù)較少，使得模擬相對誤差在柵元/網(wǎng)格內(nèi)分布不均，從而降低全局計(jì)算效率。另一個(gè)是計(jì)算邊界附近的非計(jì)數(shù)區(qū)占用的計(jì)算資源較多。由于計(jì)數(shù)區(qū)外土壤和大氣對γ射線的散射作用較強(qiáng)，大空間γ輻射場求解時(shí)在計(jì)數(shù)區(qū)外還需要非計(jì)數(shù)區(qū)。直接將權(quán)窗參數(shù)統(tǒng)一應(yīng)用于計(jì)數(shù)區(qū)與非計(jì)數(shù)區(qū)，非計(jì)數(shù)區(qū)內(nèi)（尤其截面較大的吸收介質(zhì)如土壤）模擬粒子過度分裂，占用較多計(jì)算資源，從而降低計(jì)算效率。

因此，本文在原有GVR方法上引入了體積修正和非計(jì)數(shù)區(qū)修正，解決柵元/網(wǎng)格體積差異和非計(jì)數(shù)區(qū)造成的粒子過度分裂問題，提高GVR方法在大空間γ輻射場的適用性。同時(shí)開展基于粒子誤差、權(quán)重、徑跡、數(shù)目、能量、碰撞和通量的共7種GVR方法在大空間γ輻射場的應(yīng)用研究，分析比對各個(gè)GVR的計(jì)算效率與特征參數(shù)，得到減方差效果最好的GVR方法。最后在基于通量的GVR方法上引入光滑因子SI，研究了SI對于GVR方法計(jì)算的全局權(quán)窗參數(shù)和計(jì)算效率的影響。

1 大空間γ輻射場計(jì)算幾何模型

圖1為本文建立的大空間γ輻射場全局計(jì)算幾何模型。模型主要有空氣和土壤兩種介質(zhì)組成，其中空氣密度取標(biāo)準(zhǔn)大氣密度1.205×10-3g·cm-3，土壤密度取2.35 g·cm-3。設(shè)置γ注量的計(jì)數(shù)區(qū)域?yàn)楦叨?～1 km、距源心水平距離5 km范圍內(nèi)（圖1中藍(lán)色區(qū)域）。由于計(jì)數(shù)區(qū)外的土壤和空氣對于γ射線有較強(qiáng)的散射作用，計(jì)數(shù)區(qū)外應(yīng)建立非計(jì)數(shù)區(qū)。故模型最大高度取2 km，距源心水平最大距離取6 km，土壤厚度取1 m。計(jì)算表明，繼續(xù)增加邊界尺寸對計(jì)數(shù)區(qū)內(nèi)γ注量影響可以忽略。利用MCNP（Monte Carlo N Particle Transport Code，MCNP）［2］開展計(jì)算，設(shè)置γ輻射源各向同性發(fā)射，權(quán)重為1，能譜為watt平衡裂變譜，源心距地面高度設(shè)置為20 m。為防止過度分裂造成的機(jī)時(shí)太長，文中所有模擬設(shè)置CPU總運(yùn)行時(shí)間截?cái)酁?00 min。

圖1 大空間γ輻射場計(jì)算幾何模型（彩圖見網(wǎng)絡(luò)版）Fig.1 Geometric model for the computation of large-space γ radiation field (color online)

由于GVR方法需要用到柵元/網(wǎng)格內(nèi)粒子通量、數(shù)量、權(quán)重等粒子分布信息作為輸入?yún)?shù)，本文模型將空氣與大地進(jìn)行分層處理。對于空氣，在水平方向上，每100 m分一層，共分為60層，計(jì)算所關(guān)心的0～5 km范圍內(nèi)共分50層；在垂直方向上，每100 m分一層，共分為20層，計(jì)算所關(guān)心的0～1 km高度范圍內(nèi)共分為10層。對于大地，在垂直方向上，每0.125 m分一層，共分8層；在水平方向上與大氣保持一致，每100 m分一層，共分60層。模型整體被劃分為28×60=1680個(gè)柵元，藍(lán)色計(jì)數(shù)區(qū)域共被劃分為10×50=500個(gè)柵元。

能量為1 MeV的γ射線在標(biāo)準(zhǔn)大氣中的自由程約100 m。本文γ射線平均能量略大于1 MeV，保守按照平均自由程125 m來估算，γ射線在模型中最大輸運(yùn)距離超過40個(gè)平均自由程，不采用減方差方法的情況下實(shí)際能輸運(yùn)至遠(yuǎn)處的γ射線數(shù)目極少。由于γ射線各向同性向外發(fā)射，同時(shí)射線與空氣和大地會(huì)發(fā)生多次散射作用，距離源心越遠(yuǎn)的柵元內(nèi)粒子數(shù)越少，模擬統(tǒng)計(jì)漲落越大，使得整體計(jì)算誤差呈中心小邊緣大的特點(diǎn)。中心柵元內(nèi)γ注量比邊緣柵元內(nèi)γ注量高10個(gè)數(shù)量級(jí)以上，粒子分布呈中心多邊緣少的特點(diǎn)。直接增加模擬的粒子數(shù)使得大量計(jì)算時(shí)間用于模擬誤差小的中心處的粒子輸運(yùn)，而誤差較大的邊緣網(wǎng)格處模擬粒子數(shù)不能有效增加，計(jì)算效率低。因此有必要研究GVR方法，通過設(shè)置全局權(quán)窗來增加邊緣柵元處的模擬粒子數(shù)，減小中心網(wǎng)格處的模擬粒子數(shù)，從而在保證整體無偏的情況下降低全局計(jì)算的相對偏差。

2 GVR方法與評估參數(shù)

在蒙特卡羅模擬中，問題通常被分為三類：1）源-點(diǎn)探測器響應(yīng)；2）源-區(qū)域通量或者響應(yīng)；3）整個(gè)問題空間的全局通量或響應(yīng)。前兩類問題通常采用局域減方差方法，而后一種應(yīng)采用GVR方法。GVR方法中最重要的減方差技巧就是全局權(quán)窗。通過設(shè)置權(quán)窗來控制粒子的整體分布使得整個(gè)問題空間內(nèi)的通量平均計(jì)算效率提高。通常利用確定論離散縱標(biāo)法（discrete ordinates Shifted N2Box，SN）或者蒙特卡羅（Monte Carlo，MC）方法進(jìn)行正向輸運(yùn)計(jì)算或者伴隨輸運(yùn)計(jì)算，將計(jì)算的粒子空間分布信息或粒子對目標(biāo)探測器的響應(yīng)貢獻(xiàn)轉(zhuǎn)化為權(quán)窗參數(shù)，進(jìn)而采用權(quán)窗進(jìn)行輸運(yùn)計(jì)算。根據(jù)不同的計(jì)算方法與計(jì)算內(nèi)容的兩兩組合，GVR方法通?？梢苑譃?類：1）SN方法正向輸運(yùn)［17］；2）MC方法正向輸運(yùn)；3）SN方法伴隨輸運(yùn)［18］；4）MC方法伴隨輸運(yùn)［19］。其中正向輸運(yùn)計(jì)算可以得到粒子通量的空間分布，從而設(shè)置輸運(yùn)偏倚參數(shù)，適用于源項(xiàng)小、屏蔽體厚的計(jì)算模型。大空間γ輻射場計(jì)算是典型的第三類全局通量求解問題，具有點(diǎn)源、深穿透問題等特征，適合基于正向輸運(yùn)的計(jì)算方法。

本文采用的基于MC方法正向輸運(yùn)的減方差方法通過蒙特卡羅模擬正向計(jì)算得到粒子通量、權(quán)重和誤差等粒子分布信息，進(jìn)而設(shè)置全局權(quán)窗參數(shù)來指導(dǎo)粒子輸運(yùn)。設(shè)置全局權(quán)窗的原則是在高重要性區(qū)域內(nèi)設(shè)置較小的權(quán)窗下限使得更多的粒子分裂，降低統(tǒng)計(jì)方差；在低重要性區(qū)域內(nèi)設(shè)置較大的權(quán)窗下限（Lower limit of the weight window）wth，對不太重要的粒子進(jìn)行輪盤賭，提高計(jì)算效率。當(dāng)前柵元/網(wǎng)格的重要性為［20-21］：

式中：P為所有曾進(jìn)入該柵元/網(wǎng)格的粒子對于計(jì)算結(jié)果的貢獻(xiàn)之和；wall為所有曾進(jìn)入該柵元/網(wǎng)格的粒子權(quán)重之和。

權(quán)窗下限與重要性成反比：

對于GVR方法來說，每一個(gè)柵元/網(wǎng)格對整體計(jì)算空間貢獻(xiàn)相等，即所有曾進(jìn)入每一個(gè)柵元/網(wǎng)格的粒子對于計(jì)算結(jié)果的貢獻(xiàn)之和相等。則有權(quán)窗下限：

由于柵元/網(wǎng)格內(nèi)粒子通量通常與粒子權(quán)重成正比，wth∝?i，?i為第i柵元/網(wǎng)格內(nèi)粒子平均通量。粒子通量是蒙特卡羅模擬中的常見物理量，通常采用粒子通量來計(jì)算全局權(quán)窗參數(shù)。相應(yīng)基于粒子通量的GVR方法（Flux-based GVR，GVR-F）的權(quán)窗下限為：

式中：Max(?i)為所有柵元/網(wǎng)格中最大粒子平均通量。每一個(gè)柵元/網(wǎng)格權(quán)窗上限下限之比β一般固定，比如本文所有計(jì)算取5。

另一種常用于GVR方法的粒子信息就是相對誤差。誤差大說明該處粒子數(shù)少統(tǒng)計(jì)漲落較大，應(yīng)降低權(quán)窗下限增加粒子分裂。則有基于誤差的GVR方法（Error-based GVR，GVR-R）的權(quán)窗下限為：

其中：Rei為第i柵元/網(wǎng)格內(nèi)粒子通量的誤差。在直接模擬中相對誤差與粒子通量平方根成反比，。所以基于誤差的權(quán)窗下限又可以表示為：

利用計(jì)算得到的全局權(quán)窗參數(shù)開展新的蒙特卡羅模擬，得到新的粒子分布信息。重復(fù)迭代直到結(jié)果收斂，從而得到全局精確結(jié)果。

為評價(jià)GVR方法的優(yōu)劣引入全局品質(zhì)因子FOMG和基于誤差的標(biāo)準(zhǔn)差σ。全局品質(zhì)因子為：

式中：N為柵元/網(wǎng)格的總個(gè)數(shù)；T為模擬的CPU計(jì)算時(shí)間。

基于相對誤差的標(biāo)準(zhǔn)差為：

一般來說，F(xiàn)OMG值越大計(jì)算效率越高，σ值越小平均統(tǒng)計(jì)誤差越小。GVR方法的加速比定義為GVR方法的全局品質(zhì)因子FOMG與直接模擬的FOMG之比。

3 GVR方法在大空間γ輻射場的應(yīng)用

3.1 體積修正

圖1中大空間γ輻射場計(jì)算幾何模型在徑向和垂直方向分別被劃分為60層和28層，生成相應(yīng)的非等體積的圓環(huán)形柵元。中心處柵元體積最小，徑向邊緣處柵元體積較大。采用圓環(huán)形的柵元相比于長方體柵元可以有效增加遠(yuǎn)離中心柵元的體積，從而增加?xùn)旁獌?nèi)模擬粒子數(shù)，降低計(jì)算誤差提高計(jì)算效率。圖2給出了采用圓環(huán)形柵元和長方形柵元計(jì)算的γ注量和誤差隨距離的變化。其中長方形柵元邊長為100 m，與圓環(huán)形柵元徑向間隔一致。采用等體積的長方形柵元統(tǒng)計(jì)計(jì)算的柵元內(nèi)γ注量平均誤差遠(yuǎn)大于圓環(huán)形柵元。

圖2 采用圓環(huán)形柵元和長方形柵元計(jì)算的γ注量和誤差隨距離的變化Fig.2 Calculated γ fluence and error in the circular and rectangular grid cells with respect to distance

式（4）和（6）的全局權(quán)窗下限適用于等體積的柵元和權(quán)窗，其引導(dǎo)粒子輸運(yùn)使得蒙卡模擬粒子在全空間范圍內(nèi)分布較為均勻，相對誤差在全空間范圍內(nèi)相對均勻。對于體積非均勻的柵元/網(wǎng)格，采用式（4）和（6）的全局權(quán)窗下限，體積大的柵元/網(wǎng)格內(nèi)模擬粒子數(shù)較多，體積小的柵元/網(wǎng)格內(nèi)模擬粒子數(shù)較少，使得模擬相對誤差在全局空間范圍內(nèi)分布不均。為使每一個(gè)非等體積的柵元/網(wǎng)格內(nèi)的模擬粒子數(shù)趨近相等，本文引入體積修正因子fvol，對權(quán)窗參數(shù)進(jìn)行體積修正?；谕康臋?quán)窗下限式（4）修正后為：

式中：Vi為第i個(gè)柵元/網(wǎng)格的體積；Vs為點(diǎn)源所在的柵元/網(wǎng)格的體積。

表1給出了直接模擬，采用LVR方法、GVR-F方法和引入體積修正的GVR-F方法迭代一次模擬得到的主要參數(shù)。其中采用LVR方法時(shí)參考文獻(xiàn)［1］設(shè)置水平5 km、高度500 m處的環(huán)探測器計(jì)數(shù)作為權(quán)窗生成的目標(biāo)探測器。采用LVR方法并合理設(shè)置權(quán)窗生成器可以有效降低整體誤差水平，并提高FOMG因子，降低標(biāo)準(zhǔn)差σ。由于LVR方法只針對局域進(jìn)行優(yōu)化，水平距離5 km、高度1 km附近計(jì)算誤差較大，且水平距離小于5 km、高度1 km處某些點(diǎn)誤差比直接模擬還大（圖3（b））。即使迭代三次后，LVR方法模擬得到的計(jì)算區(qū)域內(nèi)最大誤差依然有38%。

表1 采用不同減方差方法模擬所得的主要參數(shù)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方法計(jì)算得到的FOMG因子比直接模擬結(jié)果還要小。其原因?yàn)樵摲椒ㄒ龑?dǎo)模擬粒子在全局空間均勻分布，但同一高度處徑向最外層?xùn)旁w積為最內(nèi)層?xùn)旁w積的119倍，若不考慮空氣的衰減，最外層?xùn)旁獌?nèi)模擬粒子數(shù)遠(yuǎn)大于內(nèi)層?xùn)旁M粒子數(shù)，使得模擬粒子在外層?xùn)旁獌?nèi)過度分裂。圖3（c）中GVR-F方法模擬得到徑向4～5 km范圍內(nèi)的γ注量誤差明顯小于直接模擬。但由于“過度分裂”問題，在相同的CPU運(yùn)行時(shí)間內(nèi)模擬的總粒子數(shù)比直接模擬小3～4個(gè)數(shù)量級(jí)，使得徑向4 km以內(nèi)γ注量誤差比直接模擬還大。增加模擬粒子數(shù)或者提高CPU運(yùn)行時(shí)間也不能進(jìn)一步提高該方法的FOMG因子。引入體積修正后，GVR方法引導(dǎo)粒子在全柵元/網(wǎng)格內(nèi)均勻分布，防止出現(xiàn)“過度分裂”現(xiàn)象。體積修正的GVR-F方法計(jì)算的FOMG因子比直接模擬提高39倍，基于誤差的標(biāo)準(zhǔn)差σ降低約兩個(gè)數(shù)量級(jí)。

3.2 非計(jì)數(shù)區(qū)修正

大空間γ輻射場計(jì)算模型在計(jì)數(shù)區(qū)外還存在非計(jì)數(shù)區(qū)域。非計(jì)數(shù)區(qū)域內(nèi)粒子可以輸運(yùn)至計(jì)數(shù)區(qū)，影響計(jì)數(shù)區(qū)粒子通量計(jì)算。大尺度空間輻射場計(jì)算必須要合理設(shè)置計(jì)數(shù)區(qū)與非計(jì)數(shù)區(qū)。一般來說，非計(jì)數(shù)區(qū)的厚度應(yīng)達(dá)到粒子的飽和反射厚度。將GVR方法計(jì)算的權(quán)窗下限統(tǒng)一應(yīng)用于計(jì)數(shù)區(qū)與非計(jì)數(shù)區(qū)，雖然計(jì)數(shù)區(qū)計(jì)算精度得以保證，但非計(jì)數(shù)區(qū)內(nèi)（尤其本文算例中的土壤等吸收介質(zhì)內(nèi)）模擬粒子數(shù)較多，占用較多計(jì)算資源，從而降低計(jì)算效率。

根據(jù)§2，計(jì)數(shù)區(qū)內(nèi)柵元/網(wǎng)格的權(quán)窗下限與重要性成反比，與粒子平均通量成正比。非計(jì)數(shù)區(qū)內(nèi)模擬粒子只有輸運(yùn)至計(jì)數(shù)區(qū)才能對計(jì)數(shù)區(qū)粒子通量有貢獻(xiàn)。非計(jì)數(shù)區(qū)柵元/網(wǎng)格內(nèi)粒子輸運(yùn)至計(jì)數(shù)區(qū)柵元/網(wǎng)格的概率為e-s/λ，其中，s為非計(jì)數(shù)區(qū)柵元/網(wǎng)格到計(jì)數(shù)區(qū)柵元/網(wǎng)格的距離；λ為粒子的平均吸收自由程。以離非計(jì)數(shù)區(qū)最近的計(jì)數(shù)區(qū)柵元/網(wǎng)格為基準(zhǔn)，可以設(shè)置非計(jì)數(shù)區(qū)的權(quán)限下限為計(jì)數(shù)區(qū)權(quán)窗下限與粒子從非計(jì)數(shù)區(qū)輸運(yùn)至計(jì)數(shù)區(qū)的概率：

式中：wth，none為非計(jì)數(shù)區(qū)柵元/網(wǎng)格的權(quán)窗下限；wth，n為離非計(jì)數(shù)區(qū)最近的計(jì)數(shù)區(qū)柵元/網(wǎng)格的權(quán)窗下限。本文中空氣中γ的平均吸收自由程λ取250 m。該修正改變了權(quán)窗下限的數(shù)值，實(shí)際輸運(yùn)模擬時(shí)根據(jù)權(quán)窗來分裂或者殺死粒子，同時(shí)改變粒子權(quán)重來保證結(jié)果的無偏性。

圖4 給出了采用考慮體積修正的不同GVR方法計(jì)算的最靠近地面一層?xùn)旁獌?nèi)的權(quán)窗下限wth分布。權(quán)窗下限wth隨水平距離增加而迅速下降，不考慮非計(jì)數(shù)區(qū)修正，wth一直下降至6 km的計(jì)算邊界。以式（11）引入非計(jì)數(shù)區(qū)修正后，wth在計(jì)數(shù)5 km處達(dá)到最小值，而后在非計(jì)數(shù)區(qū)隨水平距離而增加。引入非計(jì)數(shù)區(qū)修正后的wth變化規(guī)律與文獻(xiàn)［1］中以5 km處注量為目標(biāo)的LVR方法計(jì)算的wth變化規(guī)律類似。

圖4 采用體積修正的不同GVR方法計(jì)算的近地面柵元內(nèi)權(quán)窗下限wth（彩圖見網(wǎng)絡(luò)版）Fig.4 Lower limit of the weight window wth in the nearground cell computed using different GVR methods with volume modification (color online)

表2給出了采用體積修正和不同非計(jì)數(shù)區(qū)修正的GVR方法迭代一次所得的主要參數(shù)。直接設(shè)置非計(jì)數(shù)區(qū)內(nèi)權(quán)窗下限wth，none=0關(guān)閉權(quán)窗，模擬得到的FOMG因子不增反降，基于誤差的標(biāo)準(zhǔn)差σ有所增加，最大誤差大幅度增加。采用非計(jì)數(shù)區(qū)修正的GVR-F方法相比于無計(jì)數(shù)區(qū)的GVR-F方法計(jì)算的大空間γ輻射場的FOMG因子提高約39%，基于誤差的標(biāo)準(zhǔn)差σ降低約27%；采用非計(jì)數(shù)區(qū)修正的GVRR方法計(jì)算的FOMG因子相應(yīng)提高約18%，基于誤差的標(biāo)準(zhǔn)差σ相應(yīng)降低約18%。由于體積修正與非計(jì)數(shù)區(qū)修正使得本文的大空間γ輻射場計(jì)算效率顯著提高，下文對GVR方法的對比研究中默認(rèn)都采用了體積修正與非計(jì)數(shù)區(qū)修正。

表2 采用體積修正和不同非計(jì)數(shù)區(qū)修正的GVR方法計(jì)算的主要參數(shù)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］還采用了基于柵元/網(wǎng)格內(nèi)粒子能量（Energy-based GVR，GVR-E）、數(shù)目（Populationbased GVR，GVR-P）、徑跡（Track-based GVR，GVRT）和權(quán)重（Weight-based GVR，GVR-W）的GVR方法，并在國際熱核聚變實(shí)驗(yàn)堆（International Thermonuclear Experimental Reactor，ITER）和中國聚變工程實(shí)驗(yàn)堆（China Fusion Engineering Test Reactor，CFETR）托卡馬克裝置中子環(huán)境計(jì)算中取得良好的減方差效果。以上方法認(rèn)為全局權(quán)窗下限與柵元/網(wǎng)格內(nèi)粒子能量、數(shù)目、徑跡、權(quán)重和碰撞次數(shù)成正比，其權(quán)窗下限計(jì)算方法與基于粒子通量的權(quán)窗下限基本一致，只需將式（4）中?i分別替換成第i柵元/網(wǎng)格內(nèi)粒子能量Ei、數(shù)目Ni、徑跡Ti和權(quán)重wi等?；陬愃圃?，本文還提出了基于柵元/網(wǎng)格內(nèi)粒子碰撞的GVR（Collision-based GVR，GVR-C）方法，同樣將式（4）中?i替換成第i柵元/網(wǎng)格內(nèi)粒子碰撞數(shù)Ci即可。

表3給出了采用不同的GVR方法模擬迭代3次后得到的主要參數(shù)。通過不斷迭代，GVR方法增加統(tǒng)計(jì)漲落較大區(qū)域內(nèi)的模擬粒子數(shù)，使得模擬粒子數(shù)在整個(gè)計(jì)數(shù)區(qū)空間中偏向于均勻分布，尤其是迭代2～3次后GVR-P方法模擬的最靠近地面柵元內(nèi)模擬粒子數(shù)基本相等（圖5）。迭代2～3次后，不同GVR方法給出的主要參數(shù)已經(jīng)收斂。GVR-W方法模擬得到的FOMG因子提高約2300倍，基于誤差的標(biāo)準(zhǔn)差σ降低約5750倍，模擬得到計(jì)數(shù)區(qū)柵元內(nèi)γ注量平均誤差為0.6%，最大誤差為1.8%，在所有GVR方法中減方差效果最好。GVR-T、GVR-P、GVR-E、GVR-C和GVR-F方法模擬得到的FOMG因子提高至2～3，加速比達(dá)到500～700倍，模擬得到的基于誤差的標(biāo)準(zhǔn)差σ降至10-5量級(jí)。GVR-R方法在所有GVR方法中減方差效果最差，但FOMG因子依然提高至0.87，加速比能達(dá)到190倍。對比圖6中迭代3次的GVR-R、GVR-W和GVR-F三種方法計(jì)算γ注量誤差分布，GVR-W模擬得到的計(jì)數(shù)區(qū)柵元內(nèi)γ注量誤差整體比GVR-R和GVR-F的誤差小。三種方法計(jì)算的柵元誤差最大值都位于水平距離5 km靠近地面處。

表3 采用7種GVR方法模擬所得的主要參數(shù)（迭代3次）Table 3 Main parameters simulated using seven GVR methods (3rd iteration)

圖5 GVR-P方法模擬的柵元內(nèi)模擬粒子數(shù)（彩圖見網(wǎng)絡(luò)版）Fig.5 Simulated particle number Nps in the near-ground cells simulated using GVR-P method (color online)

圖6 三種GVR方法計(jì)算的γ注量誤差分布Fig.6 Error distribution of γ flux computed using three GVR methods

圖7和圖8給出了GVR-W方法計(jì)算的γ注量分布和權(quán)窗下限wth分布。大空間γ輻射場計(jì)數(shù)區(qū)柵元內(nèi)γ注量呈中心高邊緣低的特征，計(jì)數(shù)區(qū)γ注量總衰減超過10個(gè)數(shù)量級(jí)。對于注量變化如此大的深穿透厚屏蔽體問題，GVR-W方法模擬的所有柵元內(nèi)γ注量誤差最大僅為1.8%，驗(yàn)證了合理的GVR方法在大空間輻射場計(jì)算中的實(shí)用性。GVR-W方法計(jì)算的計(jì)數(shù)區(qū)內(nèi)權(quán)窗下限wth同樣呈現(xiàn)中心高邊緣低的特征，計(jì)數(shù)區(qū)外wth逐漸增加，wth最低衰減至10-7量級(jí)。

圖7 GVR-W方法計(jì)算的γ注量分布Fig.7 γ flux distribution computed using three GVR-W methods

圖8 GVR-W方法計(jì)算的權(quán)窗下限wth分布Fig.8 Distribution of wth computed using GVR-W methods

不同GVR方法模擬的FOMG因子及相應(yīng)加速比不同的主要原因在于不同GVR方法計(jì)算的權(quán)窗下限wth不同。如圖9所示，GVR-R方法計(jì)算的wth最低衰減至10-4量級(jí)，GVR-W方法計(jì)算的wth最低衰減至10-7量級(jí)，GVR-T、GVR-P、GVR-E、GVR-C和GVR-F方法計(jì)算的wth最低衰減至10-9量級(jí)。由于計(jì)算的wth分屬于三個(gè)水平，不同GVR方法計(jì)算的FOMG因子也分屬于三個(gè)水平，其中GVR-W方法計(jì)算的wth衰減至10-7量級(jí)，對應(yīng)FOMG因子提高最多，全局減方差效果最好。

圖9 基于不同粒子信息的GVR方法計(jì)算的近地面柵元內(nèi)權(quán)窗下限wth（彩圖見網(wǎng)絡(luò)版）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方法減方差效果差異的原因?yàn)橛?jì)算的wth差異明顯。對于本文中的大空間γ輻射場計(jì)算模型，GVR-W方法全局減方差效果最好。但如果模型變化，比如輻射源由γ換成中子、計(jì)數(shù)區(qū)增大等，GVR-W方法計(jì)算的wth并不一定最優(yōu)。

為進(jìn)一步研究wth對全局減方差效果的影響，在式（4）中引入光滑因子SI，相應(yīng)權(quán)窗下限為：

采用光滑因子SI可以方便快捷地修改wth。當(dāng)SI=1對應(yīng)基于通量的GVR-F方法，SI=0.5對應(yīng)基于誤差的GVR-R方法。

在本文計(jì)算模型中，光滑因子SI由0.4增加至1，對應(yīng)的權(quán)窗下限wth最小值由10-3量級(jí)降低至10-9量級(jí)（圖10）。不同光滑因子的GVR方法計(jì)算的FOMG因子和加速比隨SI的增加先增加后減小，計(jì)數(shù)區(qū)柵元內(nèi)γ注量的平均誤差Rave、最大誤差Rmax和標(biāo)準(zhǔn)差σ隨SI的增加先減小后增加（表4）。

表4 不同光滑因子的GVR方法計(jì)算的主要參數(shù)（迭代3次）Table 4 Main parameters computed using GVR methods with different SI values (3rd iteration)

圖10 不同光滑因子的GVR方法計(jì)算的近地面柵元內(nèi)權(quán)窗下限wth（彩圖見網(wǎng)絡(luò)版）Fig.10 wth in the near-ground cell computed using GVR methods with different SI values (color online)

當(dāng)SI過小時(shí)，wth較大，模擬粒子在計(jì)數(shù)區(qū)內(nèi)分裂數(shù)目較少，統(tǒng)計(jì)漲落較大，模擬誤差較大。當(dāng)SI過大時(shí)，wth較小，模擬粒子在計(jì)數(shù)區(qū)內(nèi)分裂數(shù)目較多，分裂過度，計(jì)算效率下降。因此，存在一個(gè)合適的光滑因子SI，使得模擬計(jì)算效率最高。對于本文大空間γ輻射場，使得計(jì)算效率最高的光滑因子SI為0.8。當(dāng)SI=0.8時(shí)，GVR方法計(jì)算的FOMG因子最大為14.78，平均誤差最小為0.5%，標(biāo)準(zhǔn)差σ最小為1.4×10-5，加速比最大為3246，減方差效果最好。SI=0.8的GVR方法計(jì)算FOMG因子比GVR-W方法還提高40%。合理設(shè)置光滑因子SI，采用式（12）計(jì)算權(quán)窗下限wth，使用權(quán)窗開展蒙特卡羅模擬，迭代至收斂，就可以實(shí)現(xiàn)全局輻射環(huán)境的精確計(jì)算。

5 結(jié)語

為解決大空間γ輻射場的全局精確計(jì)算問題，本文開展了多種GVR方法在大空間γ輻射場的應(yīng)用研究。針對計(jì)數(shù)柵元/網(wǎng)格體積差異較大的情況，提出了體積修正方法，適用于非等體積柵元/網(wǎng)格的GVR方法的權(quán)窗參數(shù)的計(jì)算。在本文大空間γ輻射場計(jì)算中，采用體積修正的GVR-F方法計(jì)算的FOMG因子比直接模擬提高39倍，基于誤差的標(biāo)準(zhǔn)差σ降低約兩個(gè)數(shù)量級(jí)，優(yōu)于LVR方法。針對非計(jì)數(shù)區(qū)內(nèi)粒子過度分裂的情況，提出了非計(jì)數(shù)區(qū)修正方法。采用非計(jì)數(shù)區(qū)修正后，GVR-F方法計(jì)算的FOMG因子在體積修正的基礎(chǔ)上進(jìn)一步提高40%。

在體積修正和非計(jì)數(shù)區(qū)修正的基礎(chǔ)上，本文將基于粒子誤差、權(quán)重、徑跡、數(shù)目、能量、碰撞和通量共7種GVR方法應(yīng)用于大空間γ輻射場環(huán)境計(jì)算?；跈?quán)重的GVR-W方法模擬得到的FOMG因子比直接模擬提高約2304倍，在所有GVR方法中減方差效果最好。其他GVR方法模擬得到FOMG因子也能提高兩個(gè)多量級(jí)。在基于通量的GVR-F方法上，引入光滑因子SI，模擬計(jì)算的FOMG因子隨SI的增加先增加后減小。當(dāng)SI=0.8時(shí)，該方法模擬得到的FOMG因子最大，達(dá)到14.78，比直接模擬提高3246倍。對于常見的全局計(jì)算算例，存在一個(gè)合適的光滑因子SI，使得模擬計(jì)算效率最高。

本文提出的對GVR方法的體積修正和非計(jì)數(shù)區(qū)修正方法具有通用性，不僅僅適用于大空間γ輻射場計(jì)算場景，在其他非等體積計(jì)數(shù)柵元/網(wǎng)格和含非計(jì)數(shù)區(qū)的計(jì)算場景中同樣適用。

作者貢獻(xiàn)聲明劉利負(fù)責(zé)計(jì)算模型開發(fā)，程序運(yùn)行，結(jié)果分析，文章撰寫；左應(yīng)紅負(fù)責(zé)參與計(jì)算模型開發(fā)，結(jié)果分析，文章修改；牛勝利負(fù)責(zé)參與計(jì)算模型開發(fā)，指導(dǎo)結(jié)果分析，文章修改；朱金輝負(fù)責(zé)提供總體技術(shù)指導(dǎo)，文章修改；商鵬負(fù)責(zé)參與結(jié)果分析；王學(xué)棟負(fù)責(zé)參與計(jì)算模型開發(fā)。