近年来，随着城市化进程的加速，城市关停、搬迁工业企业后原址场地的污染问题凸显，对生态环境、食品安全和人体健康构成了严重威胁^[1-3]。基于这一现状，国家和地方越来越重视污染场地的环境管理，相继出台场地相关法律法规及政策标准。按照《中华人民共和国土壤污染防治法》第三十五条的规定，土壤污染风险管控和修复包括土壤污染状况调查和土壤污染风险评估、风险管控、修复、风险管控效果评估、修复效果评估、后期管理等活动。依据场地调查和风险评估的结果，对于超过可接受风险水平的点位需要建立其风险管控范围或修复范围(下统称为“修复范围”)，为场地后续修复方案的建立及修复工程的实施提供依据。由于场地污染物及污染状况的差异，从管理层面上国家尚未出台修复范围建立的技术规定，因此有必要开展这一方面的研究，以期为污染场地风险管控及修复工作提供参考。

目前常用的土壤属性制图方法有地统计学插值模型、确定性插值模型、离散点的定性分析方法等^[4-10]。地统计学插值模型是利用样本点的统计规律，使样本点之间的空间自相关定量化，从而在待预测点周围构建样本点的空间结构模型，比如克里格插值法^[11]。克里格是从变量相关性和变异性出发，在有限区域内对区域化变量的取值进行无偏、最优估计的一种方法^[12]。按照空间场是否存在漂移可将克里格插值分为普通克里格和泛克里格，其中普通克里格常称作局部最优线性无偏估计，当数据服从二阶平稳假设时，可以采用普通克里格进行插值；泛克里格是在区域化变量非平稳条件下的一种最优化无偏线性估值方法，它是从区域化变量的角度考虑数据的空间结构性，内插值更接近客观实际^[13]。确定性插值方法是基于信息点之间的相似程度或者整个曲面的光滑性来创建一个拟合曲面，如反距离加权平均插值法、趋势面法、样条函数法等。泰森多边形是常用的空间分析方法之一^[14-17]，它是由一组连接两相邻点线段的垂直平分线组成的连续多边形组成。由于环境事件发生的空间集聚性，采用地统计方法进行空间分析多数情况下会导致结果有偏，因此泰森多边形也经常用于反映区域环境因子的集聚效应^{[8, 16-17]}。

由于受历史生产活动和污染源分布等因素的影响，工业污染场地土壤中的关注污染物往往呈现局部污染严重的特点，这些区域的污染物含量极高，具有很强的空间变异性^{[2, 18]}。如何根据场地污染物空间分布选择适用的空间插值方法是影响土壤修复范围预测质量的最主要因素之一。克里格插值在数据服从二阶平稳假设时，预测精度更高；反距离插值在克里格插值无法取得可信赖的半变异函数时更占优势；泰森多边形依据污染点位的空间位置进行空间划分，可以更好地反映区域的集聚效应和空间结构。研究表明，对于某种场地土壤污染物浓度呈偏态情况下，分区预测(克里格+TIN)的预测精度更高^[19]。

锑是土壤中广泛分布的一种有毒重金属，其对人体及生物具有毒性及致癌性。中国的锑矿储量和锑生产量为世界之最，也是场地上经常出现的重金属之一^[20]，而环境中锑污染的报道却较少。锑是子弹制作时常用的硬化剂，子弹大约含锑2% ~ 5%，研究表明，在打靶场废弃的子弹被腐蚀后，分解释放的锑会进入土壤或者经淋溶进入土壤^[21]。本研究选择该场地土壤特征污染物锑作为研究对象，分析其含量的空间变异性，分别采用普通克里格插值(OK)、泛克里格插值(UK)、反距离插值(IDW)、泰森多边形(TIN)、分区预测(OK+TIN)建立场地土壤锑修复范围，从预测精度和修复范围的角度评价几种预测模型的适用性，以期为土壤锑污染状况及修复范围划定方法的选择提供参考。

1 材料与方法 1.1 研究区域概况

研究场地为我国某大型军工场地，场地占地面积约27 hm²，按照当地规划的要求，企业目前已停产搬迁。企业有近70 a的生产历史，生产初期工艺设备相对落后、经营管理较为粗放、环保设施缺少或不完善，在多年生产经营过程中，导致一些对人体有害物质的排放。根据场地生产工艺及现场调查的需要，本次调查将场地划分为A(污水排放区)、B(底火、火工区)、C(曳光弹生产区)、D(射钉弹生产区)4个区域(图 1、表 1)。该企业主要生产原料包括雷汞、铅、铝粉、铜、钾锑和石油烃等。场地地势在总体上是东北、东南两侧向河谷逐渐降低，河谷部位则呈由西向东逐渐降低的趋势。

区内地貌依成因类型可划分为堆积地形及河漫滩，河漫滩可分为高漫滩和低漫滩。场地土壤分为4层：第一层(①层)：人工堆积土，埋深为0 ~ 3.5 m，主要为低液限黏土夹炉灰、碎砖及建筑垃圾，将土壤岩性概化为低液限黏土进行风险评估；第二层(② ~ ③层)：老土层，埋深为3.5 ~ 5.5 m，主要为低液限黏土，将岩性统一概化为黏土；第三层(④ ~ ⑦层)：含水层，埋深为5.5 ~ 22.0 m，主要为级配不良中砂、粗砂组成，将岩性统一概化为砂土；第四层(⑧层)：含水层底板，埋深在22 m以下，为低液限黏土，将岩性统一概化为黏土。

该地块土壤锑污染主要集中在埋深3.8 m以上的人工堆积层，随着深度的增加土壤污染程度减轻，同时考虑到土壤岩性对污染物迁移风险的影响，结合地质分层(图 2)，本场地土壤锑风险评估及修复范围划定将地层划为0 ~ 1.0、1.0 3.8、> 3.8 m。

1.2 样品采集与分析

根据我国目前污染场地的管理要求，本次调查按照初步调查和详细调查两个步骤进行。本次初步调查主要采用判断布点，在可能造成潜在污染的热点区域进行采样点布设。结合初步调查结果，进行详细调查的布点。整个调查区域通过开展初步和详细调查共布置土壤采样点159个，最大钻探深度为31.5 m。

为了判断不同深度土壤污染情况，按照地层分布进行了不同深度的采样，取样深度范围为0.2 ~ 31 m。本研究实际的钻孔数目为159个，现场钻探的同时采用RTK对各采样点地面高程和地理位置进行测量，对采样点的坐标和高程进行记录。场地土壤钻探采用传统的钢索冲击钻法，获取的土壤剖面样品全部放置于木制岩心箱，表层土壤采样时，先清除地表石块、植被等，再进行样品采集。本次土壤调查共采集样品数量为1 264个(含100个平行样)，采集的样品装于200 ml的广口玻璃瓶内进行保存。所有样品采集完成后，置于硬质塑料样品转运箱内，加入蓝冰以确保运输过程中维持低温。样品采集完毕后送至实验室进行检测。

土壤锑的检测方法采用USEPA 6020B-2014电感耦合等离子体质谱法进行测定。样品分析取0.1 ~ 0.2 g风干土，过100目筛后，加入5 ml硝酸、2 ml高氯酸、2 ml氢氟酸，加盖并放入160 ℃石墨炉消解，开盖于约160 ℃加热至透明液滴，加入0.5 ml硼酸溶液冷却定容至50 ml待测。在样品分析时进行加标回收率和平行样的测定，标样(编号GSB04-1767-2004，批号174001)由国家有色金属及电子材料分析测试中心提供，加标样品回收率为95% ~ 97%，现场平行样的相对偏差为0.0% ~ 20.69%。

1.3 风险控制值

根据健康风险评估结果，确定锑为该场地土壤的目标污染物，需基于健康风险评估结果采用剂量-效应模型反算其土壤风险控制值。基于设定为1的可接受非致癌风险水平，通过风险评估模型反推获得其风险控制值为20 mg/kg。我国现行《土壤环境质量建设用地土壤污染风险管控标准(试行)》(GB36600—2018)^[22]规定的第一类用地锑筛选值也为20 mg/kg，因此确定场地土壤锑的风险控制值为20 mg/kg。

1.4 研究方法 1.4.1 土壤锑修复范围预测方法

以土壤锑超标地层为研究对象，结合土壤每层锑含量的分布特点，本研究选用以下5种方法分析土壤超标锑修复范围及精度：①普通克里格插值(OK)：基于普通克里格插值的“二阶平稳假设”，将原始数据进行对数转换后进行插值；②泛克里格插值(UK)：通过数据分析，去除数据来自于背景值的主导趋势，采用泛克里格插值对剔除趋势后数据进行插值；③反距离插值：采用原始数据进行反距离插值；④泰森多边形方法：作为污染场地常用修复范围的确定方法，基于采样点建立泰森多边形，以超过修复目标值的泰森多边形作为修复范围；⑤分区预测方法：基于正态分布的“3σ”原则，以锑含量超过平均值加4倍标准差的数值作为异常值，并以原始数据中值代替其高值部分，进行拆除高值后数据的普通克里格插值(OK)，最后将两部分插值数据叠加，将超过修复目标值的区域划定为修复范围。

1.4.2 插值精度验证

土壤污染物空间插值精度的验证受验证方法和验证指标的影响。本研究通过调整半方差模型的参数进行交叉验证，确定特定插值方法下的适用参数。交叉验证经常用来比较不同方法或同一方法不同参数选择的插值结果，如可以通过半方差模型中选择最佳搜索半径和步长距离或邻域点个数使克里格方差最小。

插值精度验证指标中平均误差(ME)、均方根误差(RMSE)为常用的衡量已知数据和预测数据差别的指标，根据每个验证点位置的观测值和预测值计算，公式如下：

$ {\text{ME}} = \frac{1}{n}\sum\limits_{i = 1}^n {\left[ {Z\left( {{x_i}} \right) - {Z^*}({x_i})} \right]} $

$ {\text{RMSE}} = \sqrt {\frac{1}{n}\sum\limits_{i = 1}^n {{{\left[ {z({x_i}) - {z^*}({x_i})} \right]}^2}} } $

式中：$Z\left( {{x_i}} \right)$为预测值，${Z^*}({x_i})$为样点实测值。ME越接近于0、RMSE越小表明插值精度就越高。

1.5 数据分析与处理

本研究土壤锑含量的描述性统计分析采用软件SPSS 19.0，空间插值方法及精度验证采用ESRI ArcGIS10.5进行处理。

2 结果与分析 2.1 场地土壤锑含量的统计与分析

从表 2中可以看出，场地0 ~ 1 m土壤锑含量范围为0.40 ~ 507.39 mg/kg，平均值为23.64 mg/kg，中位数为5.50 mg/kg，锑平均含量超过我国《土壤环境质量建设用地土壤污染风险管控标准(试行)》(GB36600—2018)第一类用地筛选值标准(20 mg/kg)，204个土壤样品中有37个土壤锑超标，超标率为18.14%；1 ~ 3.8 m土壤锑含量范围为0.10 ~ 260.83 mg/kg，平均值为3.78 mg/kg，418个样品中有10个土壤锑超标，超标率为2.39%。与0 ~ 1 m土层相比，1 ~ 3.8 m层土壤锑含量最大值和平均值均呈现下降趋势。

0 ~ 1 m和1 ~ 3.8 m土壤锑含量均呈偏态分布(表 2，图 3)，其偏度分别为5.15、14.4，峰度分别为30.04、246.38，极差分别为506.99、260.73 mg/kg，说明土壤锑含量具有明显的集聚性。基于OK+TIN方法识别的0 ~ 1 m土壤锑含量的异常值有4个，1 ~ 3.8 m土壤锑含量的异常值有2个，分别用其中位数替代。0 ~ 1 m、1 ~ 3.8 m替换前后的土壤锑含量对数正态累积频率如图 4所示。从图中可以看出，中值替换对每层的数据分布基本没有影响，0 ~ 1 m土壤锑含量基本均位于95% 置信区间内，数据呈对数正态分布；而1 ~ 3.8 m土壤锑含量均不服从对数正态分布，原因在于1 ~ 3.8 m土壤锑浓度整体较低，其超标率为2.4%，而最大超标倍数为12，说明1 ~ 3.8 m土壤锑含量的极端效应十分明显。

2.2 场地土壤锑含量空间分布

将土壤各采样点锑浓度与原厂区平面布置图进行空间叠加后可以发现(图 5)，0 ~ 1 m、1 ~ 3.8 m土壤锑含量具有明显的空间差异性，以0 ~ 1 m土壤尤为明显。0 ~ 1 m土壤锑含量超过筛选值特别是超标严重的采样点主要分布在以锑为原料的生产与加工区域，即场区西南部的曳光弹生产区、射钉弹生产区和东中部的底火火工区。从图 5中可以看出，土壤锑超标严重的采样点分布在B、C两个分区，锑最高浓度507.39 mg/kg出现在曳光弹生产区域0 ~ 1 m地层，1 ~ 3.8 m土壤锑含量最大值出现在C区底火加工区。原因在于这些区域都属于场地的子弹生产区，其生产过程中采用的起爆药含有硫化锑成分，生产过程中炸药的遗撒和污水排放都会导致锑在表层土壤累积。

3 土壤锑修复范围预测 3.1 模型拟合

同一插值方法中，不同参数的设定对插值精度有显著影响。影响IDW插值的因素有指数p值的选择、估计的邻域样点数目及搜索半径；影响克里格插值的因素包括数据的变异性和空间结构、半方差函数模型和搜索半径，这些方法大部分参数的确定没有太多可遵循的规律，需通过交叉验证等方法来确定。

3.1.1 IDW拟合参数

本研究通过选取不同的指数p值进行IDW插值，邻域样点数据和搜索半径根据场地生产工艺分布确定合理的数值，对插值结果进行交叉验证并计算均方根误差(RMSE)。由图 6可以看出，RMSE随着指数p值增大呈现先减小后增大的趋势，0 ~ 1 m土壤锑含量RMSE值在p=1.88处达到最小值，1 ~ 3.8 m土壤锑含量则在p=9.59处达到最小值，此时其IDW插值效果最优，因此本研究分别选取p=1.88及p=9.59进行IDW插值。

3.1.2 土壤锑含量半变异函数

普通克里格和泛克里格插值拟合的关键是变量的变异函数分析，因此对0 ~ 1 m和1 ~ 3.8 m土壤锑含量整体数据经对数转换后的数据、去掉主导趋势的对数转换数据和拆分后待插值数据分别进行半变异函数分析。块金值和基台值是半变异函数的两个重要的参数，块金值反映的是最小抽样尺度以下变量的变异性及测量误差，块金值与基台值的比值表示随机部分引起的空间异质性占系统总变异的比例。如果该比值高，说明样本间的变异更多地是由随机因素引起的，经交叉验证确定的最优模型及参数如表 3。从表中可以看出，稳态模型可较好地模拟3种情形下不同深度土壤锑含量的空间变异；0 ~ 1 m土壤锑数据去掉主导趋势后半变异模型参数与整体数据保持一致，说明0 ~ 1 m土壤锑含量主要受场地活动的影响，土壤背景含量影响有限；拆分后数据模型基台值降低，块金效应增加，说明场地土壤锑含量空间变异减少，空间相关性增加。1 ~ 3.8 m土壤锑数据去掉主导趋势后块金效应增加明显，说明1 ~ 3.8 m土壤锑主要受场地土壤随机因素的影响，来自于场地生产活动的影响减少。

3.2 不同方法对土壤锑修复范围的预测精度分析

表 4列出了各个插值方法预测统计结果与预测误差。从表中可以看出，0 ~ 1 m土壤深度OK、UK、IDW、OK+TIN预测锑含量最大值为160.70、160.69、238.48、41.72 mg/kg，最小值为2.47、2.47、1.15、2.84 mg/kg，平均值为23.36、27.84、28.30、18.23 mg/kg。克里格插值对于高值和低值的平滑作用特别明显，由于去除了极端高值的影响，数据明显呈现左侧拖尾效应，从而导致OK+TIN的平滑作用更为明显，估计的平均值也明显低于其他模型。从预测的平均误差和均方根误差来看，OK+TIN的分区预测方法要优于其他3种插值方法，说明OK+TIN对于污染物空间分布趋势有更好的精度，但由于正态分布转换过程中引起的平滑效应，导致其高值的估计误差太大，因此修复范围确定的误差较大。IDW方法对样点分布位置和样点数量要求较高，样点越多且分布越均匀其预测精度要高，而场地由于受建筑物、污染识别环节及成本的影响，很难做到高密度和均匀网格布点采样，因此插值方法对于污染范围的预测精度较低。

从表 4中可以看出1 ~ 3.8 m土壤深度OK、UK、IDW、OK+TIN预测锑含量最大值为26.72、29.16、57.42、18.89 mg/kg，预测的最小值1.68、1.69、0.26、1.50 mg/kg，平均值为6.38、7.70、6.97、5.04 mg/kg。4种插值方法对于高值的平滑作用均超过10倍，这可能与土壤锑含量不服从对数正态分布有关。OK、UK、IDW 3种插值方法的均方根误差差别不大。尽管OK+TIN具有相对低的均方根误差，但是由于本地层土壤大部分采样点锑含量较低，加上OK+TIN对于高值的平滑作用，导致插值后土壤锑含量均低于风险控制值20 mg/kg，这进一步表明OK+TIN插值方法跟OK方法一样受限于数据分布的要求。

3.3 不同方法对于土壤锑修复范围预测结果比较

对原始数据进行了反距离插值，同时对经对数转换的数据进行普通克里格和泛克里格插值，将不同方法插值后超过修复目标值20 mg/kg的区域划定为修复范围，如图 7所示。从图中可以看出，针对0 ~ 1 m土壤锑修复范围，普通克里格和泛克里格方法确定的修复范围一致，进一步验证了土壤锑主要来自于场地的生产活动。OK、UK、IDW对于低值的平滑效应也特别明显，未超标点位被纳入修复范围，从而带来修复规模的增加。与其他方法相比，OK+TIN方法可以更好地刻画污染严重区域的修复范围，但对于周围大部分清洁点存在的情况下，超标点极易被平滑为未超标点而导致修复范围的缩小。以土壤钻孔建立泰森多边形，并将超过20 mg/kg的点位所在的泰森多边形划为修复范围，表 5统计了各种方法确定的修复面积。从表中可以看出，各方法确定的修复面积大小顺序为IDW > OK=UK > TIN > OK+TIN，修复面积依次为74 154、57 427、37 338、32 707 m²，表明由于受采样方案的限制，IDW方法平滑效益最明显，将会带来修复成本的增加；OK+TIN与TIN确定的修复范围接近，但OK+TIN仍存在一定的平滑效应。从1 ~ 3.8 m各插值方法确定的场地土壤锑的修复范围如图 8从图中可以看出，OK、UK、IDW方法对于高值的平滑作用较为明显，以OK、OK+TIN尤为突出，OK+TIN仅保留了TIN的部分，OK插值导致了8个超标点均平滑掉，从而导致修复范围的减少，OK确定的修复面积也最小(表 5)；IDW方法相比较其他两种方法较为合理，但其对于未超标点的平滑作用带来修复面积的增加。

从不同深度土壤锑的模型预测精度和修复范围可以看出，克里格插值由于其对于高值和低值具有较强的平滑效应，预测结果精度较低；反距离插值由于受场地数据布点方法的影响，插值结果存在一定的局限性；OK+TIN在一定程度上弥补了克里格插值的平滑效应，但其对于高值的平滑作用仍存在，仍可能存在遗漏高风险区域范围的情况。从插值效果来看，考虑空间相关性的克里格插值对于数据质量要求较高，在数据不符合假设条件的情况下，容易带来修复面积的增加或者遗漏，这表明了地统计学插值方法对于样本数量和质量的依赖性很难满足污染物空间分布预测对于精度的要求^[23-24]。泰森多边形为反距离插值的一种极端情况，通过以土壤钻孔代表特定修复地块的方法，从几何原理上给出了较为合理的修复范围，其优势在于确定样点较少的深层土壤修复范围时，可以以钻孔代表的地块为修复范围，从而保持了不同深度超标点位修复范围的一致性。

4 结论与建议

受场地生产活动的影响，场地土壤锑含量呈明显的空间变异性，在污染源区域锑含量呈现明显的集聚性。0 ~ 1 m深度土壤锑含量主要受场地活动影响，基本接近对数正态分布，1 ~ 3.8 m土壤锑含量仅个别点位受场地活动影响，呈偏态分布。

采用普通克里格插值(OK)、泛克里格插值(UK)、分区预测(OK+TIN)、反距离插值(IDW)、泰森多边形(TIN)5种模型确定场地土壤需修复的范围。OK、UK方法可以结合数据的空间分布确定预测精度和预测误差，但其对高值和低值的平滑作用特别明显，容易导致超标点位被平滑掉或未超标点纳入修复范围，从而导致修复范围的减少或增加。OK+TIN本质上是克里格插值和TIN的结合，其通过调整数据分布进行克里格插值，因此估计精度较高，但其在大部分点位都是清洁点(1 ~ 3.8 m)的情况下很难克服OK的局限，容易导致修复范围的减少甚至遗漏。IDW的估计精度受限于样点分布位置和样本数量，在数据不均匀的情况下容易将未超标点纳入修复范围。在受成本和污染物空间分布制约的情况下，建议采用TIN方法确定污染物超标范围，其从几何原理的角度可以较为准确地划定不同深度土壤的修复范围，其突出优势在于可以实现不同深度超标点位修复范围的一致性且可操作性强。尽管TIN方法无法提供估计精度，但从业者可以结合地块的污染识别情况和增加采样点的方法进一步优化修复边界，进而较为准确地确定污染物的修复范围。

污染场地调查中，准确地反映场地污染修复范围十分重要。由于场地的复杂性和成本的约束，调查数据往往很难满足克里格插值或其变身(OK+TIN)及IDW的假设条件，因此其确定的修复范围往往与实际情况不相符。采用几何(TIN)和布点优化相结合的方法不仅可以准确地划定场地修复范围，而且对污染地块的修复有着重要的意义。