陽林峰，楊宏偉，駱志平，吳建平，程衛(wèi)亞

（1．中國原子能科學研究院輻射安全研究所，北京 102413；2．成都理工大學核技術與自動化學院，四川成都 610059）

利用核函數提高隨機游走模式中污染物濃度計算的效率

陽林峰1，2，楊宏偉1，駱志平1，吳建平2，程衛(wèi)亞1

（1．中國原子能科學研究院輻射安全研究所，北京 102413；2．成都理工大學核技術與自動化學院，四川成都 610059）

本文在拉格朗日隨機游走模式中引入核函數法代替質點法進行濃度計算，并提出了對核函數的頻寬進行變頻寬處理，將其與對流擴散時間相關聯(lián)，形成變頻寬核函數。通過對固定頻寬高斯核函數、固定頻寬伊番科尼可夫核函數、變頻寬高斯核函數和變頻寬伊番科尼可夫核函數的濃度計算結果與質點法以及解析解的濃度計算結果進行比較，發(fā)現變頻寬伊番科尼可夫核函數在不增加模擬粒子數的情況下能有效提高濃度計算結果的精確性。

核事故后果評價；隨機游走；濃度計算；核函數法；變頻寬

隨著福島第一核電站泄漏事故的發(fā)生，污染物在海洋中的對流擴散特性愈發(fā)受到外界關注［13］。在事故后，通過數值模擬模式及時、準確地給出污染物的濃度分布成為當前研究的熱點。目前，國內外用于研究污染物擴散的數值模擬模式主要有歐拉模式和拉格朗日模式兩大類［4－8］。因拉格朗日模式在數值模擬計算過程中避免了處理歐拉模式的閉合問題而被廣泛采用。在拉格朗日隨機游走模式濃度計算部分中常采用質點法來計算污染物的濃度分布［9－10］，即在統(tǒng)計粒子最終位置分布的過程中將每個數值粒子看成帶有質量的粒子，因此污染物的最終濃度分布可通過相應網格內的粒子數乘以每個粒子所攜帶的質量除以該網格的體積得到，但這種方法在粒子數較少時常因粒子對流擴散距離增加，導致計算網格內粒子數稀疏，使得最終的濃度計算結果出現大幅漲落，增加了數值計算結果的不確定度。為提高該方法計算濃度結果的精確性，可通過采取增加粒子數的方式，但增加粒子數的同時也相應地增加了數值模擬計算時間。因此，為在不增加數值模擬計算時間的基礎上，提高濃度計算結果的精確性，本工作考慮引入核函數法用于拉格朗日隨機游走模式中的濃度計算。

1 核函數計算方法

核函數法計算污染物濃度，即將對流擴散后每個粒子對最終污染物濃度分布的貢獻看成是某一函數的分布［11］。因此，采用核函數后拉格朗日隨機游走模式中粒子對空間各位置總的濃度分布貢獻可表示為：

高斯核函數：

伊番科尼可夫核函數：

式中，νd為維度d的單位球體積，νd＝2πd／2／dΓ（d／2），Γ（x）為伽馬函數。

在核函數中頻寬主要是用來度量單個粒子在進行濃度分布貢獻計算過程中的有效距離［12］。污染物輸運過程中，兩相鄰粒子間距隨對流擴散時間的增加不斷增大，若繼續(xù)采用固定頻寬的核函數，則當兩相鄰粒子間距超過兩倍核函數的頻寬時，其每個粒子類似于質點法被看成單個帶有質量的粒子。為避免這種情況，在固定頻寬的基礎上，提出對核函數中的頻寬進行變頻寬處理，即將核函數中的頻寬與對流擴散時間相關聯(lián)，其處理方式如下式：

其中，ci（i為x、y、z）為常數。

由此，通過上述處理，在固定頻寬高斯核函數和固定頻寬伊番科尼可夫核函數的基礎上，分別得到變頻寬高斯核函數和變頻寬伊番科尼可夫核函數。在拉格朗日隨機游走濃度計算中，通過這4種核函數可將單個粒子對空間各位置處的濃度分布進行函數化。

2 基于核函數的隨機游走濃度計算

2．1 濃度計算模型構建

為此，當t0＝dt時，隨機游走模型無量綱化處理后為：

式中，dW′x、dW′y、dW′z分別為平均值為0、方差為1的高斯隨機數。

對僅考慮二階全反射的相應解析解無量綱處理后為（下式中各變量省略了上標）：

式中：h為污染源排放點垂直高度；（xs，ys，zs）為污染源的釋放點位置。

2．2 不同核函數法的濃度計算結果比較

在濃度結果比較過程中，模型的輸入參數均相同，其中水平流速u＝1m／s，垂直高度h＝30m，污染源釋放點位置為（0，0，－10）。4種核函數法的濃度計算結果與質點法以及解析解所得到的濃度計算結果在中心線上時間t＝10、100、500、800個單位時的比較情況如圖1～4所示（為便于觀察，對縱坐標污染物濃度均放大1 000倍）。

圖1 t為10個單位時中心線上的濃度分布Fig．1 Concentration distribution in center line at t＝10units

從圖1可見，當t為10個單位時，即在污染物釋放源附近，由于粒子剛開始對流擴散，粒子相對集中。因此采用質點法、固定頻寬高斯核函數法、固定頻寬伊番科尼可夫核函數法、變頻寬高斯核函數法以及變頻寬伊番科尼可夫核函數法，其濃度計算結果曲線均未出現震蕩；變頻寬高斯核函數法因高斯函數本身考慮全空間區(qū)域在計算過程中有截斷，因此其濃度結果計算值較解析解的計算值低。

圖2 t為100個單位時中心線上的濃度分布Fig．2 Concentration distribution in center line at t＝100units

圖3 t為500個單位時中心線上的濃度分布Fig．3 Concentration distribution in center line at t＝500units

由圖2b可知，隨著對流擴散的進行，當對流擴散時間t增加到100個單位時，分別采用質點法和固定頻寬伊番科尼可夫核函數法計算濃度時，其濃度結果曲線已開始出現較大的震蕩，前者主要是因為隨粒子對流擴散距離的增加在部分網格內粒子數稀疏，因此采用質點法計算濃度時會出現相鄰網格之間的濃度值大幅漲落；后者盡管在計算過程中將每個粒子對濃度的貢獻函數化，但隨著相鄰兩粒子間距超過兩倍固定頻寬伊番科尼可夫核函數的頻寬時，在計算過程中每個粒子被當作單個帶有質量的粒子來計算，因此此時采用固定頻寬伊番科尼可夫核函數，其濃度結果曲線會出現震蕩。此外，從圖2a、d可見，此時采用固定頻寬高斯核函數和變頻寬伊番科尼可夫核函數法計算濃度，其濃度結果曲線與解析解的濃度結果曲線符合較好。

由圖3b可知，隨著對流擴散距離的繼續(xù)增加，當對流擴散時間增加到500個單位時，采用質點法和固定頻寬伊番科尼可夫核函數法計算的濃度值曲線震蕩加劇。此時，如圖3a、d所示，若采用固定頻寬高斯核函數和變頻寬伊番科尼可夫核函數法計算濃度，其濃度結果曲線與解析解的濃度結果曲線依然符合較好。

當對流擴散時間t增加到800個單位時，由圖4a可見，采用固定頻寬高斯核函數法計算的濃度值曲線出現震蕩，這是因為隨著相鄰兩粒子間距超過兩倍固定頻寬高斯核函數的頻寬時，粒子在計算過程中被當作單個攜帶質量的粒子，因此其相鄰濃度值之間會出現大幅漲落。相比之下，從圖4c可知，當時間增加到800個單位時采用變頻寬高斯核函數法計算濃度能有效地克服對流擴散距離增加帶來的震蕩，但高斯函數在計算過程中因質量虧損其濃度計算值會較解析解的值偏低；與此同時如圖4d所示，若采用變頻寬伊番科尼可夫核函數法計算的濃度值曲線與解析解的濃度值曲線依然符合較好。

綜上所述，在污染物釋放源附近或短距離處，因粒子對流擴散距離短、相對集中，因此無論采用質點法還是核函數法其濃度計算結果均較理想，但隨著對流擴散距離增加，質點法所計算的濃度值曲線開始出現震蕩并逐步加劇。若采用固定頻寬核函數，其抗震能力主要由核函數的頻寬決定；若采用變頻寬核函數則能有效地克服震蕩，但在使用過程中因高斯函數的質量虧損，所以變頻寬伊番科尼可夫核函數是較理想地選擇。

圖4 t為800個單位時中心線上的濃度分布Fig．4 Concentration distribution in center line at t＝800units

3 結論

在隨機游走模式濃度計算中用核函數法代替?zhèn)鹘y(tǒng)的質點法計算濃度時，通過對相應的核函數本身進行變頻寬處理，能有效克服濃度計算過程中結果隨粒子對流擴散距離增加帶來的漲落，這樣在不增加模擬粒子數的前提下，能有效提高拉格朗日隨機游走模式中濃度計算結果的精確性。

［1］ CHINO M，NAKAYAMA H，NAGAI H，etal．Preliminary estimation of release amounts of131I and137Cs accidentally discharged from the Fukushima Daiichi Nuclear Power Plant into the atmosphere［J］．J Nucl Sci Technol，2011，48：1 129－1 134．

［2］ BUTLER D．Radioactivity spreads in Japan［J］．Nature，2001，471：555－556．

［3］ 喬方利，王關鎖，趙偉，等．2011年3月日本福島核泄漏物質運輸擴散路徑的情景模擬和預測［J］．科學通報，2011，56（12）：964－971．

QIAO Fangli，WANG Guansuo，ZHAO Wei，et al．Predicting the spread of nuclear radiation from the damaged Fukushima Nuclear Power Plant［J］．Chinese Science Bulletin，2011，56（12）：964－971（in Chinese）．

［4］ KOSHIZUKA S，TAMAKO H，OKA Y．A particle method for incompressible viscous flow with fluid fragmentation［J］．Journal of Computational Fluid Dynamics，1995，4（1）：29－46．

［5］ VALE L M，DIAS J M．Coupling of a Lagrangian particle tracking module to a numerical hydrodynamic model：Simulation of pollution events inside an estuarine port area［J］．Journal of Coastal Research，2011，64：1 609－1 613．

［6］ THOMSON D J．Criteria for the selection of stochastic models of particle trajectories in turbulent flows［J］．J Fluid Mech，1987，180：529－556．

［7］ GRIFFA A，KIRWAN A D，Jr，MARIANO A J．Lagrangian analysis and prediction of coastal and ocean dynamics［M］．Cambridge：Cambridge University Press，2007．

［8］ WANG Cui，SUN Yinglan，ZHANG Xueqing．Study of the three－dimensional Lagrangian model and its application in Jiaozhou Bay［C］∥International Workshop on Education Technology and Training ＆International Workshop on Geoscience and Remote Sensing．［S．l．］：［s．n．］，2008．

［9］ RIDDLE A M．The specification of mixing in random walk models for dispersion in the sea［J］．Cont Shelf Res，1998，18：441－456．

［10］LIU Aihua，KUAI Linping．A review on radionuclides atmospheric dispersion modes［J］．Journal of Meteorology and Environment，2011，27（4）：59－65．

［11］SPIVAKOVSKAYA D，HEEMINK A W，DELEERSNILDER E．Lagrangianmodelling of multidimensional advection－diffusion with space－varying diffusivities：Theory and idealized test cases［J］．Ocean Dynamics，2007，57：189－203．

［12］AIZERMAN M，BRAVERMAN R L．Theoretical foundations of the potentialfunction method in pattern recognitionlearning［J］．Automation and Remote Control，1964，25：821－837．

［13］WADE E H．Simplified approach to particle tracking methods for contaminant transport［J］．Hydraulic Engineering，1997，123（12）：1 157－1 160．

［14］CRISTIANINI N，SHAWE－TAYLOR J．Kernel methods for pattern recognition［M］．Cambridge：Cambridge University Press，2004．

Using Kernel Function to Improve Efficiency of Contaminated Concentration Calculation in Radom Walk Model

YANG Lin－feng1，2，YANG Hong－wei1，LUO Zhi－ping1，WU Jian－ping2，CHENG Wei－ya1
（1．China Institute of Atomic Energy，P．O．Box275－15，Beijing102413，China；
2．College of Nuclear Technology and Automation Engineering，
Chengdu University of Technology，Chengdu610059，China）

In this paper，the kernel function method was introduced into the concentration calculation of the random walk model to replace the box method．And a change to the bandwidth of the kernel function was proposed to form a frequency bandwidth kernel function by relating it to convection diffusion time．By comparing the calculated concentration of analytical solution with the fixed bandwidth Gaussian kernel function，the fixed bandwidth Epanechnikov kernel function，the frequency bandwidth Gaussian kernel function，the frequency bandwidth Epanevhnikov kernel function and the box method，the conclusion was obtained that the frequency bandwidth Epanevhnikov kernel function can efficiently improve the accuracy of the contaminated concentration calculation without increasing simulated particles．

consequence assessment of nuclear accident；random walk；concentration calculation；kernel function method；frequency bandwidth

TL32

：A

1000－6931（2015）03－0485－06

10．7538／yzk．2015．49．03．0485

2013－11－28；

2014－01－21

陽林峰（1986—），男，四川自貢人，碩士研究生，輻射防護與環(huán)境保護專業(yè)