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

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

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

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

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

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

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

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

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

1 核函數(shù)計算方法

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

高斯核函數(shù)：

伊番科尼可夫核函數(shù)：

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

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

其中，ci（i為x、y、z）為常數(shù)。

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

2 基于核函數(shù)的隨機游走濃度計算

2．1 濃度計算模型構(gòu)建

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

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

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

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

2．2 不同核函數(shù)法的濃度計算結(jié)果比較

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

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

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

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

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

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

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

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

3 結(jié)論

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

［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］ 喬方利，王關(guān)鎖，趙偉，等．2011年3月日本福島核泄漏物質(zhì)運輸擴散路徑的情景模擬和預(yù)測［J］．科學(xué)通報，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è)