农田土壤重金属空间变异多尺度研究

刘 伟， 郜允兵， 潘瑜春

[1.中国矿业大学(北京)地球科学与测绘工程学院，北京 100083； 2.国家农业信息化工程技术研究中心，北京 100097；3.农业部农业信息技术重点实验室，北京 100097； 4.北京市农业物联网工程技术研究中心，北京 100097]

土壤重金属是表征土壤环境质量的重要因素，如何掌握重金属含量的空间异质性及其空间结构，对区域农田土壤环境的监控具有重要的生产意义。受错综复杂的自然因素和长期小规模、破碎化的土地利用生产方式的影响，我国土壤重金属空间变异大，空间结构复杂，当前的采样调查结果通常不足以准确地反映实际情况，更难用于对某个具体地块的污染范围划定和防治指导[1]。这些特征意味着：采样尺度过大时，小尺度的空间变异信息常常被减弱或丢失[2]；采样尺度过小时，大尺度上的空间变异由于比较微弱而被作为“随机成分”处理[3]，同时也会导致采样工作成本高、效率低下。因此，采样尺度设计的合理性是能否精准刻画土壤重金属空间异质性的关键[4]。与此同时，影响农田土壤重金属空间分布的环境因子的尺度作用范围存在差异，土壤重金属亦表现出多尺度效应[5]。而采样问题的复杂性根本上源于农田土壤重金属含量空间变异的多尺度特征[2]，单纯某一尺度上的研究并无法真实地揭示土壤重金属的空间结构特征[5-6]，有必要将土壤重金属含量空间变异的多尺度分析作为重点研究方向，其要点是确定土壤重金属含量空间异质性的特征尺度。国内外研究学者通过多尺度的采样和对比分析，对土壤重金属含量空间变异的多尺度结构作诸多的探讨[7-10]。然而，这些研究多集中在土壤重金属在不同采样尺度上的空间变异规律，将重金属的空间异质性定量地刻画到相应的尺度水平上的研究则涉及较少，尺度定量方面研究略显不足，这些不足对重金属污染区域范围的确定与风险评估显然是不利的[11]。

特征尺度是一定区域内能够表征土壤重金属空间异质性的最优空间单元[12]，一旦确定了特征尺度，就可以此确定适宜的采样间距[13]。本研究以北京市顺义区农田土壤重金属砷(As)采样点为数据源，采用尺度方差法研究该区域土壤As含量的空间结构特征，确定土壤As含量异质性的特征尺度，以期较全面地揭示农田土壤As含量的空间多尺度特征，并为特定的采样需求确定适宜的采样尺度提供参考依据。

1 材料与方法

1.1 研究区概况

顺义区地处北京市城区东北方向，地理位置为 40°00′～40°18′N、116°28′～116°58′E，总面积 1 019.89 km2。该地区属于华北平原北端，北接燕山南麓，全境地势北高南低，海拔24～637 m。境内平原为潮白河冲积扇下段，为河流洪水携带沉积物质造成，表面堆积物主要是沙、亚沙土；土壤类型主要为潮土和褐土，且潮土分布于全区大部分地区、褐土主要分布在东北部。顺义区是京郊重要的农业生产区，据统计该区2004年农业用地面积6.36万hm2，主要的种植作物为小麦、玉米、蔬菜//水果，冬小麦—夏玉米为主要的粮食种植模式，粮食生产用地占据大部分农业用地，其中南部是主要的菜地集中区域，东北部则是果园的聚集地。

1.2 数据来源与测定

土壤样点数据于2007年春季采集，采集耕层(0～20 cm)土壤样品294个，采用GPS定位记录样点中心位置，采样点主要分布于粮田、菜地、果园、苗圃等农业用地，采用ArcMap 10.1将原始采样点数据编辑生成样点分布图(图1)。所有土样在室内自然风干，碾压磨碎后，过100目尼龙网筛，采用原子吸收法测定土壤重金属As的含量。

1.3 研究方法

1.3.1 尺度方差分析 尺度方差是一种空间等级分析方法，对空间变量的多尺度结构比较敏感[14]。该方法是将研究对象的方差按尺度等级或尺度嵌套的水平逐步分解， 观察空间变量的尺度方差随尺度增加是否会发生突变。一般地说，尺度方差发生突变的尺度也是空间变异性突出的尺度，表征了该等级水平上的特征尺度，与此同时尺度方差突变的相对大小可以反映不同尺度水平空间变异对系统总体变异的贡献程度[15]。尺度方差的计算公式为：

Xijk…z=μ+α1+βij+γijk+…+ωijk…z。

(1)

式中：Xijk…z为等级系统最低层次上组成单元的值；μ表示当前层次上等级系统基本组成单元的总体平均值；αi、βij、γijk、ωijk…z为不同尺度层次的方差。选择与变量空间梯度变化相似的尺度划分方式能比较合理地刻画空间变量的尺度特征[5]。对实测样点进行趋势分析，结果表明，土壤As含量在东西和南北方向上均有梯度变化，因此采用W-E和N-S等2种尺度划分方式，设置最小尺度单元为0.36 km2，共12个尺度水平(表1)。

表1 尺度方差分析的划区方案

注：a=0.6 km。

1.3.2 多环缓冲区分析 分别以河流、畜禽养殖基地为中心按照不同缓冲半径进行多环分析，发现以河流为中心时，0.5 km 的缓冲半径可以使各缓冲区内有适宜数量的样点，以畜禽养殖基地为中心时0.3 km比较合适。因此，分别以 0.5、0.3 km的缓冲半径对河流、畜禽养殖基地进行多环缓冲区分析，将落入同一缓冲区的土壤样点归为一类，并进行统计分析。

1.3.3 数据处理 本研究数据处理方法及步骤如下：(1)运用SPSS 18.0对土壤重金属As含量样点集进行基础统计分析和正态分布检验。(2)使用GS+7.0软件分别对As含量样点进行半变异函数分析和高斯序贯模拟，其中序贯高斯模拟设置不同的种子(seed)进行8次单独模拟(G1～G8表示第1次至第8次模拟)。(3)将土壤As含量的模拟数据导入ArcGIS生成点位图后，与农田空间分布图叠置分析，筛选落入农田区域内的有效数据点，然后采用尺度方差法对农田区域内有效数据点进行多尺度分析。尺度方差分析通过R语言实现。(4)采用SPSS进行单因素方差分析，研究土壤类型、土壤质地、土地利用方式等因子对研究区土壤As含量的影响；应用ArcGIS 10.1进行多环缓冲区分析，研究河流和畜禽养殖基地对研究区土壤As含量的影响。

2 结果与分析

2.1 土壤砷实测数据及条件模拟的统计特征

由表2可知，顺义区土壤As含量变化范围为3.85～17.34 mg/kg，变异系数为27.25%，属于中等强度变异，与现有研究成果[16]相近。As含量的平均值为7.79 mg/kg，高于北京市土壤背景值7.09 mg/kg，说明顺义区As的含量受到人类活动的影响，有一定的累积。自然背景下土壤重金属含量通常符合正态分布[17]，由于受外源的影响，As含量的偏态系数大于0.5，呈现出正偏态分布，峰度系数亦达到1.97。参照Kolmogorov-Smirnov(K-S)检验结果，K-S值为1.22，双侧显著性大于0.05，符合正态分布要求。

由图2可知，模拟结果和实测数据的分布相近，变化范围相同，较充分地表征了土壤As含量的极值数据。模拟数据的均值在实测均值附近波动，这体现了条件模拟对空间数据的概率模拟，同时也表征出模拟数据和实测数据的分布总体是相同的。

表2 土壤As含量的统计结果

注：显著性水平为0.05，K-S项括号中为双侧显著性。

2.2 土壤砷空间结构分析

在土壤属性的多尺度分析中，要充分考虑变量的空间结构特征，即采样点的空间自相关性和空间变异等级结构[2]。土壤重金属As含量半方差函数参数如表3所示。采用GS+软件多次拟合半方差模型，得到Spherical为最佳模型，块金值C0为 1.853，基台值(C0+C)为4.705，基底比C0/(C0+C)为 39.4%，变程为7 km。块金值表示空间变量受随机因素引起的变异，基台值是空间变量总的变异。基底比C0/(C0+C) 表示变量随机部分引起的变异与变量总变异的比例。从结果可知，采样点基底比为0.394，具有中等空间自相关性，可进行空间变异分析。

表3 As含量半方差函数参数

概率累积曲线图不仅可以判别变量是否符合正态分布特征(正态分布情况下概率累积曲线表现为一条直线)，而且通过曲线的突变拐点，某种程度上能够揭示土壤重金属含量是否来自不同的样本总体，从而定性地判断土壤重金属含量是否具有等级结构[2，17]。由土壤重金属As含量的对数概率累计分布图(图3)可知，As含量累计概率曲线存在明显的拐点，初步判断土壤中As含量来自2个特征明显不同的总体：A总体含量水平较低，初步定位于自然背景来源；B总体含量水平较高，可能来源于人类活动的排放。其次，半方差函数图也可以识别出土壤As含量的空间等级结构特征，Robertson等通过半方差分析识别了农田土壤中pH值的空间变异等级结构[18]。由于空间变量的变异可能存在巢式等级结构，因此该变量的半方差值随着距离的增加表现为台阶式上升的趋势，而半方差值突变的拐点则刻画了不同水平上的特征尺度[14]。由图4可知，2.5 km处为半方差图突变转折点，说明土壤As含量的空间变异在2.5 km处发生了突变，表明土壤As含量空间变异存在等级结构。

2.3 多尺度空间变异识别

如上所述，研究区存在空间多尺度结构，因此，本研究采用尺度方差法计算土壤As含量在不同等级水平上的方差，并识别土壤As含量的特征尺度。在Excel导入8次条件模拟在不同尺度水平上的方差(散点)，并计算出8次条件模拟的方差在同一水平上的平均值(折线)，结果如图5所示。

由图5可知，尺度方差随着尺度的增大而表征出不同的特征，即不同尺度上具有不同的空间异质性。整体上尺度方差曲线呈现倒U形，在小尺度上随尺度增加空间变异性增大，大尺度上则相反。先看W-E划区(图5-A)的尺度方差图，随着尺度的增大，尺度方差在等级5(4a×4a)和等级9(16a×16a)处有明显的峰值，说明土壤As含量在等级5和等级9的空间异质性发生了突变，因此识别出特征尺度为 2.4、9.6 km；再看N-S划区(图5-B)的尺度方差图，尺度方差在等级3(2a×2a)、等级5(4a×4a)和等级9(16a×16a)处表现为峰值，但是等级3突变程度相对较弱，特征尺度不明显，因此只识别出特征尺度2.4、9.6 km。综合得到小尺度上特征尺度为2.4 km，与半方差函数2.5 km拐点相吻合；大尺度上特征尺度为9.6 km，与As含量半方差函数的变程 7 km 相接近。

2.4 土壤As含量多尺度效应成因分析

土壤重金属空间异质性大小的尺度效应受控于不同尺度下控制土壤重金属含量变异的各种生态过程的重要程度，即区域内不同影响因素对重金属As含量变异的作用范围不同，实际采样工作中可以根据取样的目的和关注的影响因素，选择接近影响因素作用范围的尺度作为采样尺度。例如，交通运输污染源对土壤重金属Pb含量的影响区域一般在交通线两侧几十米的范围内，影响尺度比较小[19-20]。也就是说，如果某一地区的主要污染来源于交通运输，那么道路两侧土壤重金属含量的空间变异范围(变程)是在污染源的影响距之内的，为了通过采样来全面地刻画重金属的空间变异特征，采样间距须小于影响距，这一点在柳云龙等的研究[21]中已得到证实。

本研究前述的统计分析及半方差分析都显示As含量存在中等强度空间变异，且As含量存在等级结构，空间变异受到了土壤自然背景等结构因子和人类活动等随机因子的共同影响。因此，本研究以土壤类型、土壤质地、土地利用、离河边距、离畜禽养殖基地距离等探讨各影响因素对土壤As含量空间变异的影响范围。其中，土壤类型、土壤质地、土地利用采用单因素方差分析(ANOVA)，结果显示土壤质地、土地利用对土壤As含量的空间变异存在显著影响(P<0.05)；离河边距、离畜禽养殖基地距离则采用多环缓冲区分析。

2.4.1 土壤质地 研究区土壤质地类型主要包括轻壤质、中壤质、沙壤质、沙质，294个样点中有2个样点位于其他几种类型中，为了更好地进行方差分析，本研究剔除了这2个样点。由表4可知，沙质土壤显著低于其他几类，且经过多重比较分析得出轻壤质、中壤质差异不显著，沙壤质、沙质差异也无显著差异。因此，将轻壤质、中壤质归为一类，沙壤质、沙质归一类，并将图斑合并(图6)。采用ArcGIS计算矢量斑块边长，结果显示土壤质地主要斑块边长范围4～19 km，斑块平均边长为9.1 km，与As含量空间变异在9.6 km的特征尺度相近，说明在大尺度上土壤质地是导致土壤As含量空间变异的主要因素。

表4 不同土壤质地中As的含量

注：同列数据后不同字母表示差异显著(P<0.05)。表5同。

2.4.2 土地利用类型 由表5可知，不同土地利用类型的土壤As含量依次为粮田>果园>苗圃>菜地，果园和粮田的As含量显著高于菜地。因此，将果园、粮田归为一类，并将相邻图斑合并。采用ArcGIS计算矢量斑块边长，结果显示，土地利用类型主要斑块边长范围为300～4 860 m，而土壤As含量小尺度上的特征尺度为2.5 km，说明在小尺度上土地利用类型是导致土壤As含量空间变异的重要因素之一。

表5 不同土地利用土壤中As的含量

2.4.3 河流 河流沿岸土地由于开发历史较长，土壤环境受人为的干扰较剧烈，河岸带表层土壤受重金属的污染[22]。本研究采用多环缓冲区分析法对样点进行分组统计以描述采样点与河流距离的变化关系(图7)，整体上随着样点与河流距离的增加，As的含量逐渐递减，拟合曲线在0～2.6 km范围内迅速下降，2.6 km之后逐渐趋于平稳。基于空间分析的结果显示，河流对土壤As含量的空间分布产生了一定的影响，且控制范围约为离河流2.6 km的区域内。

2.4.4 畜禽养殖 规模化、集约化的畜禽养殖会产生大量的排泄物(主要包括畜禽粪便)，由于运输距离有限，大量未经处理的畜禽粪便主要以有机肥的形式施用在养殖基地附近的农田[23]。有研究表明，华北地区畜禽粪便中As含量存在超标现象[24]。图8为采样点与畜禽养殖不同距离范围内土壤As含量变化曲线，总体上随着距离的增加，重金属含量逐渐下降，拟合曲线在0～2.2 km范围内下降比较快，2.2 km之后趋于平缓。由此推断，畜禽粪便主要投放在以畜禽养殖场为中心，方圆2.2 km的范围内，对土壤As含量空间变异的影响范围在 2.2 km 左右。

3 结论与讨论

本研究经过多尺度分析得到以下结论：(1)顺义区土壤As含量变化范围为3.85～17.34 mg/kg，平均值为 7.79 mg/kg，变异系数为27.25%，属于中等强度变异。(2)顺义区农田土壤重金属As含量具有较强的空间自相关性，空间异质性存在多尺度结构。(3)顺义区农田土壤重金属As含量空间异质性存在2个特征尺度，分别是2.4、9.6 km，并且不同的特征尺度受控于不同的影响因素，在小尺度上受土地利用、畜禽养殖和河流的影响，大尺度上主要受土壤质地的影响。

虽然本研究采用尺度方差研究了顺义区农业土壤中As含量空间变异的尺度效应，但是也存在以下不足：(1)本研究采用的基本划区尺度为0.6 km，每1级尺度上推的间距划分都是下1级的2倍，这就造成了尺度等级的不连续性，这可能会使部分特征尺度没有被挖掘出来或者特征尺度所在的尺度等级无法精确地表达，这种情况在大尺度上时会对特征尺度识别十分不利。因此，在接下来的研究可以通过改变基本划区尺度，进一步多尺度分析，以提高特征尺度识别的精度。(2)尺度涵盖粒度和幅度2个方面，本研究所述尺度仅考虑了空间粒度。为了更加全面地掌握采样尺度与土壤重金属空间异质性间的相互关系，后续的研究将就不同采样幅度下的情况进行探讨，当幅度充足时，有利于全面揭示农田土壤重金属含量的空间结构信息。(3)由于相关数据的获取难度比较大，如大气沉降对土壤As含量的影响，本研究对顺义区As含量影响因素的分析不够充分，接下来的研究中将就此展开深入的探讨。