Kriging、反距离权重和泰森多边形空间插值法在土壤苯并（a）芘超标区域模拟中的应用比较

王 君

（上海兴东环保科技有限公司，上海200000）

0 引言

苯并（a）芘作为多环芳烃中一种高环的污染物，具有较高的致畸、致癌作用[1~3]。伴随着我国经济的高速发展和工业化程度的提高，化石燃料燃烧和生物质燃烧的排放量逐渐增加，进而导致大量的苯并（a）芘排放到大气中。这些污染物会通过干湿沉降的方式进入地表，对土壤-植物生态系统造成严重的危害，甚至会通过食物摄取的方式进入人体，对人体健康产生严重影响。部分地区由于受到历史或现存工业活动的影响，土壤中存在高浓度的苯并（a）芘，这种区域被称为苯并（a）芘的污染场地[4~5]。近年来，污染场地引发的危害事件已经引起了国内外学者的高度关注。因此，准确了解污染场地土壤苯并（a）芘的空间分布特征，可为区域土壤污染治理提供理论依据，进而更精准地降低其潜在的健康风险。

为了准确揭示苯并（a）芘的空间分异特征，一些学者通过土壤样品采集的方式来评估各个场地土壤苯并（a）芘浓度的分布格局[6~7]。但是，这种方法存在一定的局限性。首先，每个采样点仅有有限的空间代表性[8]（Shi et al.，2018），该方法不能准确评估离采样点距离较远区域的苯并（a）芘的污染状况。更为重要的是污染场地通常具有污染历史的长期性和污染类型的多样性，这些因素都可能导致土壤苯并（a）芘存在较大的空间异质性[9]，仅仅使用采样点的苯并（a）芘浓度来反映整个区域的污染格局会显著增加评估结果的不确定性。因此，有必要采用一些空间建模的方法来推算模拟采样点以外区域的苯并（a）芘的浓度。

近年来，一些学者开始使用空间插值的方法来捕获高空间分辨率的污染物浓度[10~11]。孙等使用反距离权重、样条插值等方法来预测土壤Cd的浓度分布，发现这些方法模拟效果存在很大差异[12]。因此，有必要对比分析多种空间插值方法，确定最优化的建模方案，才能准确评估土壤污染物的空间特征。谢等使用局部多项式和普通克里格插值等方法来分析土壤中Cd的时空分异特征，发现克里格插值能准确捕获土壤中Cd的污染状况[13]。但是，普通克里格方法也存在一定的缺陷，它过分依赖样点数据数量和质量，忽略了与被预测属性相关的环境要素[14]。而协同克里格考虑了环境要素的影响，能明显改进估计精度及采样效率，有着广泛的应用前景和价值。

因此，本项目基于上海某苯并（a）芘污染场地的检测数据，使用协同Kriging、普通Kriging、反距离权重和泰森多边形四种方法，通过比较四种方法对污染场地土壤苯并（a）芘浓度的模拟效果，确定一种最优的建模方案来准确模拟污染场地土壤中苯并（a）芘的超标污染范围，以此为基准划定修复范围，为后续场地修复施工提供科学依据。

1 项目概况

1.1 研究区概况

本研究选择的苯并（a）芘污染场地总面积为27 000m2，为工业用地，位于上海市北部。钻探深度6 m范围内地质划分为3层：0~2 m为杂填土层；2~3 m为粉质黏土层；3~6 m为淤泥质粉质黏土层。地下水水位埋深约为1 m。

1.2 样品采集与分析

项目地块水土监测工作分为初步调查和详细调查两个阶段，两次调查共布设了57个土壤监测点，采集深度为地下0~6 m，共采集171个土壤样品。所有样品送至具有CNAS检测资质的实验室进行分析。查看样品检测结果发现，地块内苯并（a）芘含量超标较为严重，污染深度主要集中在地面以下0~3 m土壤中苯并（a）芘含量检测如表1所示。

表1 土壤中苯并（a）芘含量检测情况统计

2 建模方法

本项目采用协同Kriging、普通Kriging、反距离权重和泰森多边形4种方法对土壤样品中的苯并（a）芘浓度进行空间插值，使用十折交叉验证法[平均预测误差（Mean）和均方根误差（RMSE）]来评价4种方法的优劣。平均预测误差（Mean）越小，插值预测的无偏性越好；均方根误差（RMSE）越小，插值结果与实际结果越接近，绝对误差越小。最后选择最优方法来获取网格化的苯并（a）芘浓度并根据污染标准，划定苯并（a）芘超标（0.55 mg/kg）的范围，进而开展修复工作。四种空间插值方法的详细原理和公式如下：

选择泰森多边形算法又称为最近邻点插值法，是由荷兰气象学家A.H.Thiessen提出的一种插值分析方法，最初用于气象学中进行平均降水量的估算。其原理是将各个采样点连成三角形，作这些三角形各边的垂直平分线，将每个三角形的三条边的垂直平分线的交点（也就是外接圆的圆心）连接起来得到一个多边形。用这个多边形内所包含的采样点的苯并（a）芘浓度作为多边形区域内的苯并（a）芘浓度。

IDW法是一种局部插值法，其假设前提是未知值的点受较近控制点的影响比远控制点的影响更大。它以插值点与样本点间的距离为权重进行加权平均，离插值点越近的样品点赋予的权重越大。适用于呈均匀分布且密集程度足以反映局部差异的样点数据集。

普通Kriging法在数值具有相同的变异性的基础上，以变异函数理论和结构分析为基础，对区域化变量进行无偏最优估计，对数据属性的正态分布要求较高。

协同Kriging法协要求有一个已知的相关函数，这就需要在很多地点同时采样，测定两个函数间的相互关系，它比普通克里格法能明显改进估计精度及采样效率。本项目中使用球形模型，3个不同深度的土壤苯并芘的基台值分别为0.16、0.12和0.10，块金值为0.85、0.74和0.69，最大变程为167、168和172 m。

Z^为模拟的苯并（a）芘浓度；λ1i和 λ2i为权重；Z1和Z2为实测值；p为样本个数。

3 结果与分析

本项目采用协同Kriging、普通Kriging、反距离权重和泰森多边形四种方法对土壤样品中的苯并（a）芘浓度进行空间插值，使用十折交叉验证法[平均预测误差（Mean）和均方根误差（RMSE）]来评价4种方法的优劣。结果显示如表2所示，从0~0.5 m深度来看，苯并（a）芘浓度的平均偏差的绝对值呈现协同Kriging（0.01 mg/kg）<普通Kriging（0.02 mg/kg）<反距离权重（0.05 mg/kg）<泰森多边形（0.07 mg/kg）的趋势。RMSE也呈现协同Kriging（0.25 mg/kg）<普通Kriging（0.27 mg/kg）<反距离权重 （0.30 mg/kg）<泰森多边形（0.35 mg/kg）的趋势。从0.5~2 m和2~3 m的空间插值结果来看，其平均偏差和RMSE呈现的趋势与0~0.5 m的结果相似。结果表明，这4种空间插值方法中，协同Kriging方法具有最好的模拟效果，泰森多边形法效果最差。这主要是因为泰森多边形法将三角形覆盖区域都假设为恒定值，这样会显著增加其估算结果的不确定性。协同Kriging方法相比普通的Kriging方法考虑了模拟函数之间的相关性，能更加准确地反映苯并（a）芘浓度真实的分布状况。

表2 苯并a芘各种插值法评价指标对比表

此外，本项目还根据模拟的苯并（a）芘浓度划分了超标区域和非超标区域的范围（见图1）。结果显示，0~0.5 m、0.5~2 m和2-3 m模拟的苯并（a）芘浓度超标 （>0.55 mg/kg） 比例分别为24.65%、17.71%和7.99%。0~0.5 m、0.5~2 m和2~3 m实测的苯并（a）芘浓度超标（>0.55 mg/kg）比例分别为28.07%、15.79%和5.26%。模拟值的超标比例与实测值的超标比例非常接近，说明协同Kriging插值方法的可靠性较高，能准确反映污染场地的苯并（a）芘超标状况。

图1 基于Co-Kriging、Kriging、IDW和Voronoi方法模拟的不同深度苯并a芘的浓度

4 结论

本项目采用协同Kriging、普通Kriging、反距离权重和泰森多边形4种方法对土壤样品中的苯并（a）芘浓度进行空间插值，使用十折交叉验证法来确定最优的建模方法。结果显示，对于不同深度的土壤苯并（a）芘浓度的模拟，协同Kriging方法均取得了最好的模拟效果。基于协同Kriging插值方法的估算，0~0.5 m、0.5~2 m和2~3 m模拟的苯并（a）芘浓度超标（>0.55 mg/kg）比例分别为24.65%、17.71%和7.99%。表层土壤苯并（a）芘的超标比例最高并随深度逐渐递减。