利用核函数提高随机游走模式中污染物浓度计算的效率

阳林峰，杨宏伟，骆志平，吴建平，程卫亚

（1．中国原子能科学研究院辐射安全研究所，北京 102413；2．成都理工大学核技术与自动化学院，四川成都 610059）

利用核函数提高随机游走模式中污染物浓度计算的效率

阳林峰1，2，杨宏伟1，骆志平1，吴建平2，程卫亚1

（1．中国原子能科学研究院辐射安全研究所，北京 102413；2．成都理工大学核技术与自动化学院，四川成都 610059）

本文在拉格朗日随机游走模式中引入核函数法代替质点法进行浓度计算，并提出了对核函数的频宽进行变频宽处理，将其与对流扩散时间相关联，形成变频宽核函数。通过对固定频宽高斯核函数、固定频宽伊番科尼可夫核函数、变频宽高斯核函数和变频宽伊番科尼可夫核函数的浓度计算结果与质点法以及解析解的浓度计算结果进行比较，发现变频宽伊番科尼可夫核函数在不增加模拟粒子数的情况下能有效提高浓度计算结果的精确性。

核事故后果评价；随机游走；浓度计算；核函数法；变频宽

随着福岛第一核电站泄漏事故的发生，污染物在海洋中的对流扩散特性愈发受到外界关注［13］。在事故后，通过数值模拟模式及时、准确地给出污染物的浓度分布成为当前研究的热点。目前，国内外用于研究污染物扩散的数值模拟模式主要有欧拉模式和拉格朗日模式两大类［4－8］。因拉格朗日模式在数值模拟计算过程中避免了处理欧拉模式的闭合问题而被广泛采用。在拉格朗日随机游走模式浓度计算部分中常采用质点法来计算污染物的浓度分布［9－10］，即在统计粒子最终位置分布的过程中将每个数值粒子看成带有质量的粒子，因此污染物的最终浓度分布可通过相应网格内的粒子数乘以每个粒子所携带的质量除以该网格的体积得到，但这种方法在粒子数较少时常因粒子对流扩散距离增加，导致计算网格内粒子数稀疏，使得最终的浓度计算结果出现大幅涨落，增加了数值计算结果的不确定度。为提高该方法计算浓度结果的精确性，可通过采取增加粒子数的方式，但增加粒子数的同时也相应地增加了数值模拟计算时间。因此，为在不增加数值模拟计算时间的基础上，提高浓度计算结果的精确性，本工作考虑引入核函数法用于拉格朗日随机游走模式中的浓度计算。

1 核函数计算方法

核函数法计算污染物浓度，即将对流扩散后每个粒子对最终污染物浓度分布的贡献看成是某一函数的分布［11］。因此，采用核函数后拉格朗日随机游走模式中粒子对空间各位置总的浓度分布贡献可表示为：

高斯核函数：

伊番科尼可夫核函数：

式中，νd为维度d的单位球体积，νd＝2πd／2／dΓ（d／2），Γ（x）为伽马函数。

在核函数中频宽主要是用来度量单个粒子在进行浓度分布贡献计算过程中的有效距离［12］。污染物输运过程中，两相邻粒子间距随对流扩散时间的增加不断增大，若继续采用固定频宽的核函数，则当两相邻粒子间距超过两倍核函数的频宽时，其每个粒子类似于质点法被看成单个带有质量的粒子。为避免这种情况，在固定频宽的基础上，提出对核函数中的频宽进行变频宽处理，即将核函数中的频宽与对流扩散时间相关联，其处理方式如下式：

其中，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个单位时，即在污染物释放源附近，由于粒子刚开始对流扩散，粒子相对集中。因此采用质点法、固定频宽高斯核函数法、固定频宽伊番科尼可夫核函数法、变频宽高斯核函数法以及变频宽伊番科尼可夫核函数法，其浓度计算结果曲线均未出现震荡；变频宽高斯核函数法因高斯函数本身考虑全空间区域在计算过程中有截断，因此其浓度结果计算值较解析解的计算值低。

图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可知，随着对流扩散距离的继续增加，当对流扩散时间增加到500个单位时，采用质点法和固定频宽伊番科尼可夫核函数法计算的浓度值曲线震荡加剧。此时，如图3a、d所示，若采用固定频宽高斯核函数和变频宽伊番科尼可夫核函数法计算浓度，其浓度结果曲线与解析解的浓度结果曲线依然符合较好。

当对流扩散时间t增加到800个单位时，由图4a可见，采用固定频宽高斯核函数法计算的浓度值曲线出现震荡，这是因为随着相邻两粒子间距超过两倍固定频宽高斯核函数的频宽时，粒子在计算过程中被当作单个携带质量的粒子，因此其相邻浓度值之间会出现大幅涨落。相比之下，从图4c可知，当时间增加到800个单位时采用变频宽高斯核函数法计算浓度能有效地克服对流扩散距离增加带来的震荡，但高斯函数在计算过程中因质量亏损其浓度计算值会较解析解的值偏低；与此同时如图4d所示，若采用变频宽伊番科尼可夫核函数法计算的浓度值曲线与解析解的浓度值曲线依然符合较好。

综上所述，在污染物释放源附近或短距离处，因粒子对流扩散距离短、相对集中，因此无论采用质点法还是核函数法其浓度计算结果均较理想，但随着对流扩散距离增加，质点法所计算的浓度值曲线开始出现震荡并逐步加剧。若采用固定频宽核函数，其抗震能力主要由核函数的频宽决定；若采用变频宽核函数则能有效地克服震荡，但在使用过程中因高斯函数的质量亏损，所以变频宽伊番科尼可夫核函数是较理想地选择。

图4 t为800个单位时中心线上的浓度分布Fig．4 Concentration distribution in center line at t＝800units

3 结论

在随机游走模式浓度计算中用核函数法代替传统的质点法计算浓度时，通过对相应的核函数本身进行变频宽处理，能有效克服浓度计算过程中结果随粒子对流扩散距离增加带来的涨落，这样在不增加模拟粒子数的前提下，能有效提高拉格朗日随机游走模式中浓度计算结果的精确性。

［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—），男，四川自贡人，硕士研究生，辐射防护与环境保护专业