污染场地的环境调查、风险评估以及污染分布确定等工作都是基于场地采样样点数据来进行^{[1, 2, 3]}，场地采样样点钻孔数据由于其带有空间坐标信息，是一种典型的地学变量，可以进行地学上的相关模拟和预测^[4]。随着计算机技术和空间信息技术的发展，各种插值模型和地理信息系统等行业软件被广泛应用于污染场地的空间分布预测和污染分布制图^{[5, 6, 7, 8]}。污染场地污染物空间分布的界定和污染修复土方量计算的准确与否对修复方法的选择和修复成本都有着重要的影响^[9]。目前，污染场地污染物的空间分布范围界定大多基于2维平面，利用场地表层或某单一地层土壤钻孔样点数据，采用地统计学模型或确定性插值模型进行空间分布预测并进行污染制图^{[10, 11]}，对于单一地层污染空间分布界定能够取得较好的预测效果^[12]，但在整个场地污染范围界定过程中没有考虑到不同地层之间污染物迁移的特性，受土壤异质性等因素的影响，预测结果存在着一定的不确定性^{[13, 14]}，不同地层之间误差的累积，导致最后污染范围的确定和污染土方量的计算存在着误差。

自20世纪90年代以来，随着对3维空间插值技术的深入研究以及3维可视化系统的出现和功能的不断完善，使得地学空间变量的真实3维模拟和3维插值成为了可能^{[15, 16]}。3维空间插值和可视化技术根据地层内的采样点空间分布信息，可以将区域化变量插值成空间体，能够获取变量在真实3维环境下的空间分布，使地层建模和3维插值技术在相关领域取得了较好的应用效果^{[17, 18, 19, 20]}。通过地层建模和3维插值，可以直观获取污染物的3维空间分布信息，降低2维插值模型在高程变化较大地区和污染分布离散性强的地层中产生的预测误差。受不同3维插值模型算法及适用范围等因素影响，不同插值模型对污染范围的界定存在一定的不确定性，因此选取最优插值模型，有效规避对预测结果的不确定性影响，对准确界定污染场地污染范围和受污染土壤土方量具有重要的理论和实际意义。本研究以某典型焦化企业污染场地为目标研究区域，选择场地中特征污染物Bap(苯并(a)芘，Benzo(a)pyrene)为研究对象，对比研究不同3维插值模型对预测土壤中Bap分布范围界定不确定性的影响，为场地修复治理提供科学依据。

本研究目标污染场地所在地区曾是北京市以化工行业为主的工业区，随着污染企业的转型和改造搬迁，该焦化企业也已关停，周边土地利用类型逐渐改为公建与居住用地。该焦化企业曾是我国规模最大的独立炼焦化学工业企业之一，是国内最大的商品焦炭供应和出口基地，主要为大型冶金、化工企业提供各种规格的优质焦炭，在建厂初期生产工艺技术条件和污染治理水平所限，环境污染相对较严重，其中不乏对人体有害和致癌物质排放，对焦化厂厂区和周围环境造成一定影响。该焦化厂总面积1.35 km²，主要生产车间包括备煤分厂、炼焦一分厂、二分厂、三分厂、筛焦分厂、回收一分厂、二分厂、煤气精制分厂、精苯分厂、焦油分厂等。产品包括焦炭、焦炉煤气，同时还生产焦油、苯、硫铵、沥青、萘等40多种化工产品。通过场地的初步环境调查，确定多环芳烃为场地内特征污染物。

根据对目标场地的实地勘察，结合原有水文地质和岩土工程勘察资料，将本场地内最大勘探深度19.30 m范围的土层划分为人工堆积层和第四纪沉积层两大类，并按地层岩性及其物理力学性质指标，进一步划分为6个大层，第1大层为人工堆积的粘质粉土填土、粉质粘土填土；第2大层为粘质粉土、砂质粉土；第3大层为粉砂、细砂；第4大层为粉质粘土、粘质粉土；第5大层为细砂、粉砂；第6大层为粉质粘土、粘质粉土。

第1次现场采样采用判断布点的原则，其目的是在场地污染识别的基础上，选择潜在污染区域进行土壤布点采样，共布设64个土壤采样点；第2次现场采样是在第1次采样点布设的基础上，采用近似网格布点和判断布点相结合的方法，目的是全面了解场地污染分布情况，判断布点是结合第1次采样分析结果，对污染严重区进行加密。鉴于场地中有部分硬化地面和建筑物，两次采样没有完全按照网格布点的方法进行。两次共采集有效土壤样点114个。为了判断土壤中污染物含量随深度的变化情况进行不同深度的取样，取样深度为0.3—17.0 m，每个钻孔的深度根据污染源的情况有所不同。本次采样将取样深度分为6层：0.3—1.5 m，1.5—4.0 m，4.0—8.0 m，8.0—10.0 m，10.0—13.0 m，和13—17.0 m。钻探采用本地区常用的SH-30型钻机，采样使用原状土取土器按照方案设计深度取土，取土后采样。多环芳烃的测定方法以及质量控制参照USEPA8270D标准^[21]。场地污染土壤钻孔采样分布如图 1所示。

图 1 场地土壤钻孔样点分布 Fig. 1 The distribution of soil drilling samples

图选项

本研究选择插值精度评价中常用的交叉验证法来评价各种插值模型的精度，采用平均误差(ME)和均方根误差(RMSE)作为误差统计指标，ME越接近于0，RMSE值越小，插值精度就越高。

点击浏览公式

式中，u(x_i)为预测值，u^*(x_i)为原始样点含量值。

本研究对钻孔数据的地层建模、3维空间插值和3维可视化表征均在软件EVS．Pro中实现。对钻孔样点数据进行前期预处理，将钻孔数据转换成EVS．Pro可读取的地层文件，用于构建研究区的地层建模，主要架构过程为，插值地表面及各个地层的底面，每个地层的底面即为下层的顶面，再将每两个相邻的地层间填充为某一地层，结合构建的3维地层模型，采用不同的3维空间插值方法对污染物钻孔样点含量数据进行3维插值和分布表征。

不同地层土壤钻孔样点含量数据的描述性统计特征结果如表 1所示。从表 1可以看出，在每层中样点含量的最小值与最大值差异较大。前3层样点含量数据的最大值超标严重，均大于后3层最大值的数百倍，表明污染物在场地分布状况受一定迁移转化因素的影响。将采样样点数据叠加到原厂区平面图上可以看出，样点含量的高值点主要分布在焦油分厂、炼焦分厂等重点污染厂区，在局部区域存在严重污染现象。6层样点含量数据的偏度和峰度都较大，变异系数均超过了200%，表明样点含量数据具有很高的偏倚性和空间变异特征。

在污染场地土壤及地下水调查中，仅依靠柱状图或2维平面图提供的信息，不能完全满足获取整个场地地层分布状况的需求，通过钻孔数据可以获取钻孔样点的准确地层分布信息，但受钻孔样点个数的限制，并不能掌握整个场地的地层状况，通过有限的钻孔样点进行建模，是模拟整个场地地层分布的主要手段，因此，3维地层建模技术在污染场地的环境调查中具有很强的技术优势，通过构建3维地层模型，可以快速获取不同地层的分布状况和分布规律，加强对场地地质特征的认识，同时也可以辅助分析污染物在不同地层的分布和迁移规律。为说明钻孔样点数据对模拟地层几何体的过程，将钻孔高程数据与地层几何体显示在同一坐标系当中。根据本焦化场地的钻孔样点数据，首先提取出不同层的高程点，在EVS．Pro中通过Post sample模块生成钻孔样点分布图，3D Geology模块可以读取记录土壤类型的钻孔数据文件，自动进行空间网格grid的划分，创建不同地层，地层属性包括土壤类型、高程等，Geologic surfaces模块用来接收3D geology 模块生成的网格，节点模型，展现地层的曲面，通过上述方法构建的底层模型即第6大层粉质粘土、粘质粉土层如图 2a所示，同样方法，建立不同层的层面，根据初步地质调查结果，本场地内地层岩性较好，相邻土层之间没有空隙，对相邻土层之间进行填充，即可获得上层的地层模型，同理可构建不同土层的层面模型，最终形成了场地不同土层的3维地层建模，如图 2b所示，图 2b中从上到下蓝色、浅蓝、浅绿、绿色、橙色、红色分别代表粘质粉土填土、粉质粘土填土层；粘质粉土、砂质粉土层；粉砂、细砂层；粉质粘土、粘质粉土层；细砂、粉砂层和粉质粘土、粘质粉土层。通过地层建模，可以实现地层的3维立体显示和查看，也可以在任意方向对地质体截取切片，了解地质结构在不同平面上的分布特征，结合特征污染物的3维插值结果，通过切片的截取，能够显示污染物在不同剖面上的污染情况。

图 2 场地不同地层的三维地层模型 Fig. 2 3D stratum model of the contaminated site

图选项

在对场地污染土壤钻孔数据预处理以及3维体模型构建基础上，选用该建模系统中Krig-3D、IDW-Shepard、IDW-(Franke/Nielson)以及Nearest Neighbor 4种3维插值方法，对Bap的钻孔样点含量数据进行3维插值。应用克里格模型插值之前，需要对待插值数据进行半变异函数的拟合，根据理论半变异函数拟合的结果，对特征污染物进行3维空间插值计算。苯并(a)芘Bap污染含量钻孔数据最终3维插值结果如图 3所示。从污染物的3维空间插值分布图可以看出，4种方法的插值结果差异较大，图 3d中显示表层污染严重区域最大，图 3a、3c、3b次之，表明污染插值结果受插值方法的影响较大。通过切片显示可以看出，第2层的插值结果具有和表层相同的规律。污染严重区域的钻孔样点为本场地高风险的集中点，4种插值结果的污染超标严重区域主要集中在场地的中下部，结合该焦化厂原有厂区平面布局图以及历史生产活动可知，该厂区中下部主要有炼焦分厂和焦油分厂等生产车间，是多环芳烃污染产生的主要原因，个别车间在生产、存储过程中的泄露、遗洒等原因使得局部地区污染严重。另外，污染物在不同地层中也出现不均匀性分布特征，如个别钻孔样点的下层含量超过上层土壤的含量，说明污染物在场地中的分布受累积释放因素以及不同地层迁移特征的影响。

图 3 场地土壤多环芳烃三维空间分布预测 Fig. 3 3-Dimensional spatial distribution prediction of soil PAHsa：采用Krig-3D; b: IDW-Shepard; c: IDW-(Franke/Nielson); d: Nearest Neighbor模型计算的3维空间插值结果; Bap: 苯并(a)芘Benzo(a)pyrene

图选项

对不同插值模型插值结果的内部特征进行统计分析，可以确定污染物在3维空间不同地层上的污染分布范围以及污染土方量的估算。EVS．Pro中的Plum volume模块接受插值结果，按照临界值的标准界定出需要表达的污染范围。EVS．Pro将3维应用分析与ESRI ArcGIS进行了无缝集成，污染体可通过网格表面的形式在ArcGIS中呈现，体现了对数据编辑管理的优势，也可以对污染土壤的边界点进行标示。参照北京市污染场地土壤筛选值所规定的标准，定义苯并(a)芘Bap污染超标的临界值为0.4 mg/kg，4种插值模型界定的污染物在不同地层污染超标的范围如图 4所示。从图 4可以看出，4种插值模型界定的污染超标范围主要集中在第1和第2层上。局部污染严重的区域主要分布在表层土壤中，位于场地中下部，叠加土壤钻孔样点数据，通过不同方向上的切片显示可以看出，含量超标严重的样点均在插值结果中局部污染严重的区域，说明污染物在场地中的分布状况和污染超标范围受场地历史生产、管理布局和污染成因的影响。大气排放是本焦化场地的主要污染途径，排放产生的污染物通过扩散对整个场地表层土壤产生污染，场地的中下部为回收一分厂、焦油分厂、炼焦分厂等车间，是产生污染的主要原因，因此界定的污染超标范围也主要分布在该区域的表层土壤中，另外，由于储存了含多环芳烃的煤化工产品，储罐渗漏、遗洒也会对表层土壤造成局部严重污染，并可能渗透到下层土壤中。由于多环芳烃在土壤和地下水中很难进行迁移转化，在本场地中其污染主要分布在第1和第2地层当中。4种插值模型界定的污染范围和受污染土壤的土方量差异较大，参照北京市污染场地土壤筛选值所规定的标准，计算4种插值模型在场地中污染超标土壤的土方量分别为8.51×10⁵、5.62×10⁵、7.12×10⁵、1.09×10⁶m³，Krig-3D和IDW-(Franke/Nielson)模型计算受污染土壤土方量较为接近，但IDW-Shepard和Nearest Neighbor模型计算结果与上述两种模型有较大的偏差。

图 4 场地土壤多环芳烃污染污染空间分布范围 Fig. 4 Scope of soil PAHs pollutiona: 采用Krig-3D; b:IDW-Shepard; c: IDW-(Franke/Nielson); d: Nearest Neighbor模型界定的污染范围; Bap: 苯并(a)芘Benzo(a)pyrene

图选项

不同的插值模型虽然可以应用于该焦化场地污染物的3维空间分布预测，但预测的污染范围以及计算出受污染土壤的土方量差异较大，对预测结果产生一定的不确定性，通过对预测结果的精度验证，可以选取最优插值模型，提高预测结果的精度。本文采用交叉验证方法中的平均值误差(ME)和均方根误差(RMSE)两个指标来比较不同插值模型的精度。4种插值模型插值精度的交叉验证结果如表 2所示。通过交叉验证结果可以看出，Krig-3D的ME和RMSE最小，取得了最高的预测精度，而Nearest Neighbor计算出的插值精度最低。插值结果的精度验证表明Krig-3D模型界定的污染范围与场地污染的实际情况最为接近，取得了最优的预测效果，其界定的污染范围能够有效指导场地修复范围和开挖边界的确定。相对于其它3种3维空间插值模型，Krig-3D模型最适合于本场地土壤多环芳烃的3维空间分布预测。

利用场地污染土壤的钻孔样点数据进行场地的地层建模并对污染样点含量数据3维空间插值，可以实现污染物在真3维环境下的污染分布表征，界定污染物在不同地层上的污染范围和受污染土壤的土方量。3维空间插值及污染分布的可视化可以直观体现污染物在不同地层上的污染特征，以及污染特性在地层间的相互关系和变化规律。

不同插值模型受其算法和场地钻孔数据特征的影响，插值结果存在较大差异。Krig-3D模型取得了最高的预测精度，能有效降低对预测结果不确定性的影响。界定的污染范围与实际情况较为接近，为场地修复治理和开挖边界的确定提供了重要参考。

地层建模和3维插值可视化技术为污染场地的风险评估、修复治理等工作提供了一个新的思路，改变了传统的2维平面工作模式，在真3维环境下揭示污染物的分布特征和分布规律，能够更为真实掌握污染物在场地土壤中的分布状况，为后续的污染场地修复治理等相关工作提供科学指导。