王振宇，黄伟奇,*，孙 健，彭广宇

(1.陆军防化学院，北京 102205；2.湖南交通职业技术学院，湖南 长沙 410132)

掌握辐射剂量场的分布对实施辐射防护行动具有重要意义，掌握辐射剂量场分布的方法主要可分为计算和监测两大类：计算类有点核积分、蒙特卡罗、离散纵标等方法；监测一般指部署固定或移动辐射测量设备来获取辐射场数据。当放射源的活度、位置等信息未知时，辐射剂量场分布的调查主要依靠辐射探测器监测，监测方法的缺点是监测数据分布离散、成本高。而插值算法能在少量监测数据的基础上，通过构建数学模型对未知辐射剂量场分布进行重构，研究插值算法有助于快速、低成本地解决辐射剂量场分布问题。

插值算法是一种常见的数据处理方法，在各领域均有不同程度的应用。在进行辐射监测时，合理运用插值算法可提高计算效率。Ródenas等[1]和Wang等[2-3]分别采用线性插值和网函数插值完成了监测数据呈网格均匀分布时的辐射剂量场的插值重建。但实际上监测数据呈离散分布的情况更普遍，此时相对复杂的反距离权重法[4-6]、克里金(Kriging)法[7]和径向基函数法[8-9]能进行辐射剂量场的插值重建。此外还有一些考虑到辐射场特性的改进插值方法，如张敏[10]在反距离平方加权法的基础上提出的辐射剂量值占优构造还原法。

泊松克里金法属于改进的插值算法，最早由Monestiez等[11]在调查长须鲸数量分布时提出，Zhao等[12]将其应用于辐射场的插值重构。实际上，泊松克里金法并没有改变克里金法的求解方法，全局插值优化时仍存在最优解无法收敛的问题，此外泊松克里金法在辐射监测领域的应用较少，参数影响效果、使用范围与效果有待进一步研究讨论。本文尝试在泊松克里金优化计算中采用代理模型的方法，通过实验分析参数影响的效果，并用不同条件下的辐射场进行验证，为相关研究提供参考。

1 泊松克里金算法

1.1 泊松克里金算法原理

(1)

(2)

泊松克里金算法是克里金插值算法的一种特殊应用，主要区别是泊松克里金算法用观测数据取代真实数据进行无偏估计。泊松克里金算法主要应用于真实样本数据难以获取的情况，其假设是观测值与真实值之间符合泊松分布的关系。在辐射场测量中，由于放射性衰变规律和辐射探测器原理的限制，真实的剂量率无法测得，只能获得其观测值x。根据辐射场计数泊松分布理论，辐射探测器获得的剂量率观测值也服从泊松分布关系[12]，即：

x|d～Poisson(d)

(3)

故本文采用泊松克里金算法，即用样本的剂量率观测值xi对辐射场真实剂量率进行无偏估计，估计值的表达式示于式(4)。

(4)

求解出满足条件的一组最优权重系数ci，即可求得辐射场真实剂量率的估计值。

1.2 基于代理模型的泊松克里金算法求解

克里金插值算法的权重系数方程组的求解通常采用拉格朗日乘数法，该方法只能得到局部最优解。前期计算中发现，当辐射场剂量率数据多时，会出现计算时间长、最优解无法收敛的情况。本文将克里金代理模型用于泊松克里金算法的求解。代理模型是指可替代复杂真实模型的一类简化近似数学模型，克里金代理模型是代理模型的一种。本文在泊松克里金理论的基础上对克里金代理模型进行适当修改，同时将模式搜索算法用于最优解的收敛优化。

在构建的泊松克里金代理模型中，某点的辐射场真实剂量率d可表示为回归模型mTβ与随机误差e的和[13-14]：

d=mTβ+e

(5)

式中：mT为回归模型基函数，为行向量；β为回归系数，为列向量；e为误差函数。同时在泊松克里金算法中，辐射场剂量率观测值x服从期望为真实值d的泊松分布，根据全期望公式E(x)=E[E(x│d)]与总方差公式var(x)=E[var(x│d)]+var[E(x│d)]，可求得剂量率观测值x的期望与方差：E(x)=mTβ、var(x)=mTβ+e。因此剂量率观测值x也可表示为回归模型与整体误差之和(式(6))，其中的mTβ和e的取值与式(5)不同。

x=mTβ+e

(6)

因此，在泊松克里金代理模型中，真实值的估计值可用矩阵表示为式(7)。

(7)

式中：回归模型的基函数mT和误差函数e为可选择的已知参数，基函数mT为可变化阶次的多项式函数，后续计算中用矩阵M表示一组基函数mT；误差函数e之间的相关性用相关系数表示，其中矩阵R表示样本数据的相关系数；向量cT为权重系数。回归系数向量β与权重系数向量cT为未知参数，需单独求解。

回归系数β：为确定最优参数β，可根据样本点辐射剂量率观测值X进行最大似然估计(X为x的矩阵形式)。在模型中，辐射场剂量率是一个n维二元高斯分布问题，样本观测剂量的似然函数可构建为式(8)。

(8)

式中：L为似然函数；σ2为高斯分布的方差；n为坐标向量的维度(本文取n=2)，求解式(8)后可得：

(9)

(10)

将式(9)、(10)代入似然函数(式(8))中，忽略常数项后可将似然函数化简为：

(11)

在最大似然估计中，似然函数(式(11))为目标函数，相关系数矩阵的行列式|R|为自变量。通过模式搜寻算法，可逐步逼近最优相关系数矩阵R，进而也可求出回归系数β。通过相关系数矩阵R也能确定相关系数模型的具体参数，利用该相关系数模型可求解待估计点的相关系数rT。

(12)

同样通过拉格朗日乘数法求解后，辐射场真实剂量率的估计值表达式为：

(MTR-1M)-1MTR-1X

(13)

(14)

2 泊松克里金法参数影响分析

为研究泊松克里金算法中参数的影响，基于SuperMC软件构建一个室内多放射源虚拟辐射场。该虚拟辐射场尺寸为40 m(长)×30 m(宽)，含两个活度为3.7×108Bq 的137Cs和活度为1.85×108Bq的60Co放射源。通过计算得到距地面1 m高度处的当量剂量率，如图1所示。采用datasample函数从图1中随机均匀抽取100个数据作为样本数据。该样本数据用于分析泊松克里金代理模型中各参数变量对辐射场插值还原结果的影响。

2.1 多项式基函数的影响

多项式基函数的阶次影响回归模型拟合的整体趋势，拟合结果通过判定系数(R2)来评价，R2是线性回归中对估计的回归方程拟合优度的度量，其计算公式如式(15)所示。不同阶次多项式的拟合结果列于表1。

(15)

式中：SSE为残差平方和；SST为总离差平方和。

表1 回归模型拟合结果Table 1 Regression model fitting result

由表1可知，多项式阶次越高，判定系数R2越大，即拟合效果越好，但系数个数增加，计算更加复杂。其中二次多项式拟合与三次多项式拟合结果中R2相差较小，系数个数却相差较大。综合上述因素，选择二次多项式为回归模型的基函数。

2.2 误差相关模型的影响

不同点位处的误差具有相关性，这也是不同位置辐射剂量率的相互关联在数学模型上的具体体现。常用的相关模型列于表2。

图1 距地面1 m高度处的当量剂量率Fig.1 Equivalent dose rate at height of 1 m above ground

表2 相关模型函数Table 2 Correlation model function

图2 各相关系数模型函数曲线Fig.2 Correlation coefficient model function curve

由图2可看出，总体上相关系数R随距离d的增大而减小，参数θ主要影响距离d与相关系数R之间的函数关系。当θ<0.1时，无论距离d大小，相关系数始终接近于1；当θ增加时，呈现距离近处的相关系数大，距离远处的相关系数小，当θ增加至一定值时，相关系数随距离的增大迅速减小，且距离稍远处的相关系数为0，此时继续增加θ对模型函数的影响较小。如在EXP、LIN、SPHERICAL、CUBIC、SPLINE模型中，θ=20的函数曲线与θ=50的函数曲线十分接近，GAUSS模型中θ=50的函数曲线与θ=100的函数曲线接近。因此在后续泊松克里金代理模型的最优相关系数求解中，可将模式搜索算法中θ的取值限定在一定范围内，在EXP、LIN、SPHERICAL、CUBIC、SPLINE模型中θ∈(0.1,20)，GAUSS模型中θ∈(0.1,50)。在此基础上，采用不同的相关系数模型对原辐射场进行插值还原。主要参数如下：采用二阶多项式回归，GAUSS模型θ∈(0.1,50)，初始θ=25，其余模型θ∈(0.1,20)，初始θ=10。各模型插值还原的效果示于图3，图中黑点为样本数据。

图3 不同相关模型对辐射场插值还原效果Fig.3 Radiation field reconstruction effect of different related models

在CUBIC和SPLINE两种模型结果中，相关系数随距离的增加迅速减小，与其他模型相比，在同一距离处有更小的相关系数，因此同一峰上的不同数据被识别成多个单峰。EXP、LIN和SPHERICAL三种模型形成的数据不够平滑，不符合实际辐射场连续场的特征，其中EXP模型形成的辐射峰值底部较大，不符合辐射峰较尖锐的特征，会错误地增加辐射场区域面积。而GAUSS模型形成的峰与原始数据较一致，且平滑效果较好。此外，对比样本数据点可发现，GAUSS模型还原的辐射场最高剂量率为2.5 mSv/h，明显高于样本数据2 mSv/h，且与原始数据趋势一致。在实际应用中，这种寻找样本数据范围外的热点区域的功能十分重要。综上所述，GAUSS相关系数模型较适用于γ射线在空气中形成的辐射场插值还原。

3 实验验证

3.1 真实辐射场的构建与测量

实验场地选在室内，场地尺寸为15 m×20 m×4 m。以1 m为间隔在实验场地面上建立x-y坐标系，用标记贴在网格点处顺次编号标记。本实验共设2组，放射源均为137Cs。单源组放射源活度为1.746 MBq；双源组中一号放射源活度为0.786 MBq，二号放射源活度为0.96 MBq。

采用6150AD-b闪烁体探测器(卡迪诺科技贸易有限公司，国防科技工业电离辐射一级计量站于2013年12月26日校准)进行测量，探测器能量响应范围为23 keV～7 MeV，最大误差为±30%。实验测量高度为0.5 m，与放射源之间的测量角度保持在±80°以内，测量过程中保证仪器与放射源间无障碍物。

经过对原始数据进行判弃筛选，取某点处连续3次所读取数据的平均值作为该点的最终测量剂量率。经过数据处理，单源辐射场、双源辐射场的三维分布和等高线图分别示于图4、5，其中单源辐射场最高剂量率为5.7 μSv/h，双源辐射场2号峰剂量率为2.8 μSv/h、3号峰剂量率为3.3 μSv/h。

图4 单源辐射场分布Fig.4 Radiation field of single source

图5 双源辐射场分布Fig.5 Radiation field of dual sources

3.2 基于小范围测量数据的算法验证

小范围辐射场数据取自上述实验测量数据。为研究样本点数量对辐射场的还原效果，在样本数据量N=20、40、60、80、100、150、200、300时，进行插值计算并绘制等剂量率线图。单源辐射场等剂量率线图示于图6，双源辐射场等剂量率线图示于图7。

平均误差E用于定量对比效果的影响，所有网格点位处误差e的平均值为平均误差E，误差e计算公式如式(16)所示。为消除样本数据随机选取过程对误差带来的影响，对每种情况选取100个样本数据，并进行计算分析，结果示于图8。

(16)

图6 单源辐射场等剂量率线图Fig.6 Contour map of single source radiation field

图7 双源辐射场等剂量率线图Fig.7 Contour map of dual sources radiation field

由图6、7可知，随着样本点数量的增加，辐射场剂量率热点区域逐渐收敛，同时辐射剂量率峰也逐渐尖锐。图8中，当样本点数量为10时，单源和双源辐射场的平均误差分别为27.63%和26.30%；当样本点数量增加到150时，平均误差下降较快，分别降至4.53%和4.24%；样本点数量为150～300时，平均误差下降较慢，最终仍分别存在2.56%和2.46%的平均相对误差。综上，在10%的相对误差允许范围内，泊松克里金算法还原辐射场的最小样本量为网格点数的1/10。

3.3 基于大范围(福岛周边)数据的算法验证

为验证该算法在大范围辐射场的插值重构效果，从日本原子能研究开发机构(JAEA)网站获取了福岛第一核电站80 km范围内空间剂量率的监测结果[15]，监测时间为2016年8—10月，如图9a所示。在插值重构中分别以监测数据坐标经纬度范围(东经140.125 3°～141.037 9°、北纬36.710 6°～38.165 3°)为x、y坐标范围，对该区域重新划分成200×300的网格。然后通过随机函数选取600个剂量率数据点作为样本数据。采用算法插值重构的辐射场分布如图9b所示。

由图9可知，采用算法插值重构的剂量率分布(图9b)中峰值的分布走向与原始剂量率分布(图9a)一致，说明在大范围辐射场应用中，辐射场热点区域能基本识别出来。计算显示，其平均相对误差为141.69%，表明热点定位效果不如小范围辐射场准确。此外，在大范围辐射场计算中，形成可靠的辐射场剂量分布所需的样本点数量较多，分析其主要原因是福岛周边大范围辐射场中放射源的位置分布不规则，放射性核素浓度分布不均匀，形成的真实辐射场较复杂。

图8 平均误差与样本点数量的关系Fig.8 Influence of sample data size on error

图9 原始剂量率分布(a)[14]和采用算法插值重构的剂量率分布(b)Fig.9 Original dose rate distribution (a)[14] and dose rate distribution reconstructed by algorithm (b)

4 结论

本文对代理模型求解的泊松克里金算法进行了参数分析，并通过实测辐射场数据、福岛核电站周边辐射场数据进行了算法应用验证。结果表明，该算法在小范围简单辐射场重构中效果较好，由30个监测点数据构建的辐射场的相对误差在10%左右；在大范围复杂辐射场重构中效果略差，由600个监测点数据构建的辐射场的相对误差为141.69%。本文计算结果初步验证了泊松克里金算法在不同条件下辐射剂量场重构应用中的可行性，证明该方法在快速、低成本解决未知放射源辐射场的重构方面中具有一定的潜力和研究价值。