近年来，由于人类活动强度的提高，农田土壤重金属污染已成为全球关注的环境问题之一^[1-2]。在污染土壤的重金属元素中，Cd、Pb、Hg、As、Cr、Zn、Cu、Ni等8种元素的研究报道较多^[3-6]。值得关注的是，即使在较低含量水平条件下，土壤Cd、Pb、Hg和As也不仅会影响农产品的产量和品质，还会通过食物链传递，对人体健康造成潜在威胁^[5]。作为我国粮食、油料等作物主要生产区之一的四川盆地丘陵地区，其农田土壤重金属的研究已有少量报道^[7-9]。例如，在尚二萍等^[7]的研究结果中，四川盆地耕地点位超标率为43.55%，而严明书等^[8]和刘洁等^[9]对重庆都市经济圈土壤及沱江中游土壤质量的研究表明，其点位重金属超标率分别为13.59%和0.6%。尽管样点超标率可以描述土壤重金属的污染状况，但不同研究结果较大的空间异质性也增大了准确评估区域尺度重金属污染分布特征的难度。因此，进一步开展该区域土壤重金属空间分布及污染评价研究，为区域土壤环境质量管理提供参考依据是十分必要的。

关于区域土壤重金属污染评价方法，传统上多采用基于样点的单因子污染指数法（P_i）、内梅罗综合指数法（P_N）、地累积指数法（I_geo）和潜在生态风险指数法（RI）等方法进行评价^[10-12]。然而，重金属在土壤中的移动性较差，基于样点的指数法难于提供其污染空间分布的真实信息^[13]。特别是在丘陵地区，重金属在同一农田地块内及不同地块间的迁移易受地形的限制。随着丘陵地区土壤重金属污染防治进入规划与工程实施阶段，提供更准确的污染空间分布信息成为必然要求。因此，近年来一些研究者基于GIS、GPS和地统计学等空间分析方法研究了区域土壤重金属的空间分布特征^{[9, 14-16]}。同时，也有研究者将这些空间分析法与指数法结合对其进行评价^[17]。然而，迄今基于点位的污染超标率与基于空间分布的超标率之间的关系研究仍少见报道。此外，在相同的研究尺度（丘陵区）、采样方案及软件平台下，系统性地对比不同空间插值方法对农田土壤重金属超标率的影响仍鲜见报道。

凯江流域为四川盆地丘陵区涪江流域的一级支流，是四川盆地典型的粮油及蔬菜生产基地之一。流域内种植业和畜禽养殖业发达，人口密度高达500~ 700人·km^-2。近年来，随着农业集约化程度的提高和乡镇企业的发展，其农田土壤重金属的含量水平及污染状况尚不清楚。因此，本文根据我国《土壤环境质量农用地土壤污染风险管控标准（试行）》（GB 15618—2018）^[18]，采用单因子指数法，在GIS和GPS软件平台上，系统对比研究反比距离权重法（Inverse distance weighted，IDW）、径向基函数法（Radial basis function，RBF）、普通克里格法（Ordinary kriging，OK）、析取克里格法（Disjunctive kriging，DK）和概率克里格法（Probability kriging，PK）5种空间插值方法对该区域农田土壤重金属污染评价结果的影响，旨在为区域土壤污染治理提供科学参考。

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

研究区域地处东经104°40′ ~105°6′和北纬30° 56′ ~31°4′，面积约471 km²，属亚热带湿润性季风性气候，年平均气温15 ℃，年均降水量800 mm，具有气候温和、四季分明、雨水充沛等特点。研究河段处于四川盆地丘陵低山地区，其土地利用方式主要为耕地，以种植水稻、小麦、棉花和油菜等为主。成土母质以侏罗系蓬莱组紫色砂泥岩、白垩系城墙岩群黄色砂岩、厚砂薄页岩风化坡残积物为主，凯江主流及支流沿岸有部分地区为第四系新冲积物。

1.2 土壤样品采集与分析方法

根据该区地形地貌、成土母质、土壤类型、土地利用等条件，按随机网格法和分层抽样法进行样点的布设，共187个采样点及20个检查点（图 1）。在每个样点按五点采样法采集同一地块0~20 cm表层土样，混合均匀后用四分法保留土样约1 kg；样点位置采用GPS进行定位，并调查记录地形、母质、土壤类型等自然环境条件，以及作物、灌溉、施肥等社会经济信息。

土壤样品于实验室自然风干后，去掉砾石、动植物残骸等，进行研磨，再分别过20目、100目筛，待测。土壤pH值在水土比为2.5:1条件下采用复合电极法进行测定；土壤Cd、Pb的总量用HNO₃-HF-HClO₄体系（V:V:V=1:2:2）消解，消解液冷却后，加入5 mL 2 mol·L^-1 HCl溶解，并转移至25 mL容量瓶，定容、过滤。待测液采用石墨炉原子吸收分光光度法测定；土壤Hg和As总量采用HCl与HNO₃混合液体系（V:V=1: 1）消解，原子荧光法测定^[19]。在分析过程中所用试剂均为分析纯，所有土壤样品中每个重金属全量测定时均设置3个重复，并设置空白样，检测过程中采用国家标准土壤样品（GSS-3）进行质量控制，Cd、Hg、As和Pb的平均回收率分别为106.8%、92.5%、96.8%和102.3%。

1.3 单因子污染指数法

单因子污染指数法作为最简单的评价方法之一，通过实际含量和评价标准的比值直接表现出重金属的污染级别及危害程度。其计算公式^[7]为：

(1)

式中：P_i为第i种土壤重金属的环境质量指数；C_i为土样中第i种重金属的实测含量，mg·kg^-1；S_i为第i种重金属的农用地土壤污染风险筛选值（试行）（GB 15618—2018）^[18]，mg·kg^-1。本文评价分为两个等级：当P_i < 1时，土样重金属含量低于筛选值，样点基本未受到污染；当P_i≥1时，土壤样点受到污染，指数值越高，污染的程度越高。

1.4 空间插值方法

本研究基于样点GPS坐标和ArcGIS 10.2地理信息系统软件平台，采用地统计学的普通克里格法进行空间插值获得区域土壤4种重金属含量的空间分布图。同时，为系统对比研究不同插值方法对空间分布的影响，在相同的样点GPS坐标和ArcGIS 10.2平台上，分别采用反比距离权重法和径向基函数法两种确定性方法，以及普通克里格法和析取克里格法两种地统计学方法对单因子指数P_i进行空间插值^[20]，以研究不同空间插值方法对该区域农田土壤重金属污染评价结果的影响。其中，基于地统计学的普通克里格法和析取克里格法均要求数据为正态分布。此外，本研究还采用概率克里格法^[20]对单因子指数P_i进行了不同污染水平的风险分析。最后，根据以上5种空间插值方法获得的信息，综合讨论该区域农田土壤重金属污染超标状况。

1.5 数据处理方法

本研究所有数据统计均在SPSS 19.0平台上进行；样点分布图和各类空间插值图利用ArcGIS 10.2软件绘制。

2 结果与讨论 2.1 土壤中重金属含量统计

研究区土壤样点pH范围为4.90~8.25。土壤中Cd、Hg、As和Pb 4种重金属均呈正态分布（表 1），平均含量分别为0.38、0.10、7.77 mg·kg^-1和22.54 mg· kg^-1。其中，土壤Cd、Hg、As和Pb的平均含量分别为当地土壤环境背景值的1.58、1.43、1.27倍和1.11倍，表明该区域存在不同程度的重金属累积情况。根据GB 15618—2018^[18]，4种重金属元素的平均值均低于其风险筛选值，但有部分土壤样点的Cd和As含量超过其污染风险筛选值。变异系数（CV）相较于标准差更能反映出数据的离散性^[21]，土壤Cd、Hg、As、Pb的变异系数分别为50%、40%、29%和13%，均处于中度变异程度（10% < CV < 100%）^{[13, 22]}，表明这几种元素数据的离散性较小，分布相对较为均匀。

2.2 土壤重金属含量空间特征分析

由于该区域土壤Cd、Hg、As和Pb含量均符合正态分布，因此本研究在ArcGIS 10.2平台上选择地统计学模块中的普通克里格法进行空间插值，获得它们含量的空间分布图（图 2）。研究区土壤Cd（图 2a）含量在西北角和凯河一带有两个高值中心（>0.8 mg· kg^-1），同时呈现向四周递减的趋势。研究区西北角靠近县城，工业活动频繁；同时，在凯河附近曾有一养猪场，其养殖废弃物可能是造成污染的主要原因。土壤Hg的空间分布（图 2b）呈现以南华北部为高值区（>0.17 mg·kg^-1），向东、向北逐渐减少到0.06~0.10 mg·kg^-1的低值区，并在东南部最低（< 0.06 mg·kg^-1）。由于研究区土壤Hg平均含量较低，因此该分布可能与土壤母质因素有较大关系。土壤As的空间分布（图 2c）呈现以南华和万安的两个高值中心（>11 mg· kg^-1），向四周递减的趋势。研究区内土壤Pb（图 2d）的高值中心（>25 mg·kg^-1）位于南华、凯河处，并向周围呈现递减趋势，其东北及西北部属于低值区（< 20 mg·kg^-1）。有研究表明，土壤中As和Pb主要源于交通污染和大气沉降^[23-24]。研究区土壤As和Pb的高值区均处于交通较好的区域，同时在南华附近存在一个垃圾焚烧厂，这也可能是造成高值的原因。

2.3 土壤重金属污染评价 2.3.1 基于样点的单因子污染指数评价

研究区内土壤Cd、Hg、As和Pb的单项污染指数的平均值分别为0.487、0.111、0.381和0.102（表 2），其中仅土壤Cd和As存在P_i≥1的超过风险筛选值的样点，其占比分别为7.5%和0.5%。因此，Cd是研究区农田土壤中主要的重金属污染物。

2.3.2 基于空间的单因子污染指数评价

为进一步探究研究区重金属污染指数的空间分布特征，分别采用反比距离权重法、径向基函数法、普通克里格法和析取克里格法对4种重金属的P_i进行空间插值，并分级统计其面积占总面积的百分比（表 3）。

总体上4种插值方法的平均误差（ME）、标准化平均误差（MSE）和20个校验点的检验误差（MEP）均较小（表 3），表明它们用于插值是可行的，但普通克里格法的误差相对较小。从土壤重金属种类分析，Hg和Pb的4种空间插值分析方法的统计结果完全相同，全区域均为P_i < 0.4等级的面积分布。因此，这两种土壤重金属在研究区内无污染风险。然而，不同空间插值方法所得的土壤Cd和As的面积统计结果却呈现一定差异（表 3）。

土壤As仅存在P_i < 0.4、0.4≤P_i < 0.6和0.6≤P_i < 0.8这3个污染指数等级的空间分布。因此，尽管这4种空间插值方法所得的分布等级相同，但在前3个等级各自分布的面积却呈现一定差异。例如，在0.6≤P_i < 0.8等级的分布中，占区域面积百分数由小到大依次为析取克里格法面积（仅0.07%） < 径向基函数法（0.10%） < 普通克里格法（0.67%） < 反比距离权重法（0.89%），与基于样点污染指数的评价结果（2.7%）不同。此外，4种空间插值方法皆因仅有1个超标样点而无法形成超标等级的面积分布。土壤Cd污染指数的空间分布则因插值方法不同呈现明显差异。但无论哪一种空间插值方法，其超标区域占研究区总面积的百分比（0.14%~1.06%）均明显低于基于样点统计的百分比（7.5%）。

为进一步观察4种空间插值方法对其污染指数，特别是超标指数空间分布的影响，本文选取在土壤质量环境评价研究中采用最频繁的反比距离权重法^[25]和普通克里格法^[9]作图，进行辅助分析（图 3和图 4）。对比图 3、图 4和表 3可知，反比距离权重法和径向基函数法中土壤Cd的P_i≥1超标等级面积（1.06%和0.67%）明显高于普通克里格法（0.42%）和析取克里格法（0.14%），这是由于算法公式上的差异所致^[26-27]。前两种方法属确定性插值方法，均突出高值点对其邻近空间区域影响，其权重更大，所以插值土壤Cd的超标面积偏大（图 3a）；反之，基于地统计学的普通克里格法和析取克里格法并不强调高值点对其邻近空间区域的影响权重，且不同含量样点取等权重，因而插值结果具有一定的平滑效应^{[20, 28]}。特别是析取克里格法，其平滑效应明显，导致其超标面积仅为反比距离权重法和普通克里格法的13.2%和33.3%。

2.4 单因子指数的风险概率评价

与前面4种研究单因子污染指数值空间分布特征的插值方法不同，概率克里格法^[20]能基于污染指数的某一阈值水平来观察研究区超过该阈值水平概率的空间分布特征^[29]。毛竹等^[30]利用该方法分析了汉源富泉铅锌矿区土壤重金属的概率空间分布，预测出研究区各区域重金属的超标概率。本研究基于土壤Cd、Hg、As、Pb的含量水平及单因子污染指数统计结果（表 2），通过预测概率分析，将Cd和As的阈值设定为1.0，而Hg和Pb的最大阈值分别为0.3和0.2，以超过这些阈值为条件，用概率克里格法依次进行空间插值获得它们的风险概率分布图（图 5）。其中，Hg在P_i>0.3及Pb在P_i>0.2条件下均无高风险概率，因此整个研究区农田土壤受Hg和Pb污染的概率很低，该高值区域评价结果与前面基于样点或空间的单因子指数评价结果均一致。

土壤Cd与As污染超标（P_i>1）的风险概率区域占整个研究区面积的比例则与前面几种评价不完全一致。其中，在万安附近，有一小块农田土壤As的超标概率达到0.8~1.0的水平（图 5c，占研究区面积的1.05%），即该区域土壤有80%的概率受到As污染。其评价结果与基于样点的单因子污染指数评价相似，而未像基于空间的前4种插值方法一样漏掉有1个样点土壤As超标的信息。同时，土壤Cd存在两个概率为1.0的高值区域，其空间分布与P_i值的超标区域空间分布特征相似，但面积约为区域面积的0.16%，更接近基于空间的析取克里格法的评价结果。

3 结论

（1）研究区土壤Cd、Hg、As和Pb均呈正态分布，平均含量依次是当地土壤环境背景值的1.58、1.43、1.27倍和1.11倍，存在不同程度的重金属累积情况，其中有部分土壤点位Cd和As含量超过GB 15618—2018中的风险筛选值。

（2）采用普通克里格法获得区域土壤中4种重金属含量的空间分布图。研究区内土壤Cd、As和Pb的含量总体上呈现以两个高值中心向四周递减的趋势；土壤Hg则总体呈现西部高、东部低的趋势。

（3）基于样点的单因子指数评价得出土壤污染程度为Cd>As>Hg>Pb。基于空间的面积统计结果表明，研究区仅土壤Cd存在超标面积。4种方法中，普通克里格插值方法的误差最小，更接近于真实情况。

（4）概率克里格空间插值表明研究区仅有土壤Cd和As存在高污染风险概率区域，其空间分布特征与P_i的高值区域空间分布特征基本一致。