农田土壤重金属空间变异多尺度分析<br/>——以北京顺义土壤Cd为例

农田土壤重金属空间变异多尺度分析
——以北京顺义土壤Cd为例

2019-01-23郜允兵周艳兵潘瑜春戴华阳高秉博阎跃观

农业环境科学学报 2019年1期

刘伟，郜允兵，周艳兵，潘瑜春，戴华阳，高秉博，阎跃观

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

土壤重金属是表征土壤环境质量的重要因素，受错综复杂的自然条件、工矿企业生产、长期小规模/细碎化的农业生产管理方式和交通运输等因素的影响，我国土壤重金属空间变异大，空间结构复杂，当前的采样调查结果通常不足以准确地反映实际情况，更难用于对某个具体地块的污染范围划定和防治指导[1]。因此，掌握土壤重金属的空间分异性及其空间结构特征，对区域农田土壤环境的全方位防治与修复具有重要的现实意义。

影响农田土壤重金属空间变化的环境因素存在不同的空间作用范围[2]。工矿业生产排放的废气、废水通常具有相对较小的空间影响范围，而土壤母质、土壤类型以及地形地貌等因素通常在中观宏观尺度等较大范围内体现出趋势性[3]。上述特征使得土壤重金属表现出多尺度效应，存在复杂的空间多尺度变异性特征[4]，局限在某单一尺度上的研究无法准确有效地揭示土壤重金属的空间结构特征[5]。近年来，国内外研究学者在多尺度效应方面做了较多研究，王幼奇等[6]研究了引黄灌区不同采样尺度下农田土壤重金属的空间分布特征，分析了两种尺度的空间变异差异。霍霄妮等[7]和黄银华等[8]分析了北京、广州农业耕作层土壤重金属多尺度效应，使用多尺度嵌套模型刻画重金属的空间结构特征。王景云等[9]、胡孙等[10]和Lv等[11]分析了城郊农业土壤重金属不同尺度上的污染来源及尺度差异。上述研究重在分析论证重金属存在多尺度空间异质性，但在土壤重金属多尺度定量刻画分析方面略显不足。由于土壤重金属尺度定量化研究的不充分，给采样点布设带来一定困难，同时也增加了基于监测样点的空间插值制图的不确定性，不利于土壤重金属污染分级分类监管、防治和修复等工作的开展。因此，有必要开展土壤重金属含量空间分异的多尺度定量分析研究，确定土壤重金属含量空间异质性的特征尺度。采样尺度设计是精准刻画土壤重金属空间异质性的关键点[12]，特征尺度是一定区域内能够表征土壤重金属空间异质性的最优空间单元[13]，一旦确定了特征尺度，就可以此确定适宜的采样间距[14]。本文以北京市顺义区农田土壤重金属Cd采样点为数据源，采用尺度方差法研究该区域土壤Cd含量的空间结构特征，确定土壤Cd含量空间变异的特征尺度，以期较全面地揭示农田土壤Cd含量的空间多尺度特征，并为特定的采样需求确定适宜的采样尺度提供参考依据。

1 材料和方法

1.1 研究区概况

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

1.2 数据来源与测定

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

1.3 研究方法

1.3.1 尺度方差分析

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

式中：Xijk...z为等级系统最低层次上组成单元的值，μ表示当前层次上等级系统基本组成单元的总体平均值，αi、βij、γijk、ωijk...z为系统中各等级水平上的影响。选择与变量空间梯度变化相似的尺度划分方式能比较合理地刻画空间变量的尺度特征[4]。对实测样点进行趋势分析，结果表明土壤Cd在南北方向上均有梯度变化，因此采用N-S尺度划分方式。顺义区菜地、果园等零碎农用地斑块面积平均值为0.36 km2，且细碎图斑在研究区占比不大，因此本文确定以0.36 km2作为研究的最小尺度单元，以确保不遗漏小尺度的空间变异。根据最小尺度单元及顺义区的总面积进行计算，采用12个尺度水平（表1）可覆盖顺义全区，各尺度水平的计算公式如表2。

1.3.2 影响因素的获取与处理

图1 研究区土壤采样点分布Figure 1 Location of soils sample sites in studied area

表1 尺度方差分析的划区方案（a=0.6 km）Table 1 Zonal systems for scale variance analysis（a=0.6 km）

表2 尺度方差构成Table 2 Scale variance components

为了研究影响因素对土壤中Cd含量空间变异的作用范围，本文结合相关土壤重金属来源成因分析论文[17-19]，筛选了成土母质、土壤类型、土壤质地、土地利用类型、农业管理措施等作为候选因素，其中北京地区耕作层的成土母质均为第四世纪黄土，从土壤发生学看并无明显岩石岩性上的差异，经分析土壤质地对重金属含量无显著差异性，具体可见表5。土地利用类型和农业管理措施可以通过土地利用强度进行表征，因此本文选取2因素5变量进行尺度方差分析（表3）。其中，土壤类型数据来自北京市1∶50 000的土壤类型图。土壤类型分为普通褐土、潮褐土、普通潮土、褐潮土4类，本文以每平方公里内各土壤类型所占面积表示各土壤类型的密度。土地利用数据来源于北京市1∶50 000的土地利用现状图（2006年），根据不同土地的农业化学品、有机肥料的使用量来表示土地利用强度，对土地利用数据分类得到4种不同利用强度的土地：高强度（设施农业、菜地、果园）、中强度（耕地、苗圃）、低强度（林地、草地）、闲置地，然后给不同利用强度的土地赋予分级指数进行计算每平方公里土地的利用强度（强度分级指数见文献[20]），计算公式见表3。影响因素空间化及因子变量的量化均由ArcGIS 10.1的空间分析模块计算得到，因素因子数据均采用栅格数据，像元空间分辨率为100 m×100 m。

1.3.3 数据处理

本文数据处理方法及步骤如下：（1）运用SPSS 18.0对土壤重金属Cd的采样数据进行基础统计分析和正态分布检验。（2）采用SPSS进行单因素方差分析，研究土壤类型、土壤质地、土地利用类型等因子对研究区土壤Cd含量的影响。（3）使用GS+9.0软件分别对Cd采样点进行半变异函数分析和高斯序贯模拟，其中序贯高斯模拟设置不同的种子（Seed）进行8组单独模拟，8组模拟计算依次产生了100、200、500、800、1000、2000、3000、5000模拟数据，取8组的平均值作为最终的结果。（4）将土壤中Cd的模拟数据导入ArcGIS生成点位图后，与农田空间分布图叠置分析，筛选落入农田区域内的有效数据点，然后采用尺度方差法对农田区域内有效数据点进行多尺度分析。尺度方差分析通过R语言编程实现。（5）使用Geoda1.12进行空间自相关性分析，研究土壤类型和土地利用强度对研究区土壤Cd含量空间变异的作用范围。

表3 土壤重金属含量影响因素Table 3 Influencing factors of soil heavy metals

2 结果与讨论

2.1 土壤Cd实测数据的统计分析

经分析（表4），顺义区土壤Cd含量变化范围为0.015～0.469 mg·kg-1，变异系数为44.9%，属于中等强度变异，与王纪华等研究成果相近[21]。Cd的平均值为0.136 mg·kg-1，高于北京市土壤背景值0.119 mg·kg-1，说明顺义区Cd的含量受到人类活动的影响，有一定的累积。自然背景下土壤重金属含量通常符合正态分布，由于受外源的影响，Cd的偏态系数大于0.5，呈现出正偏态分布，峰度系数亦达到6.36。参照K-S检验结果，K-S为3.03，双侧显著性小于0.05，不符合正态分布，这可能是Cd受人类活动影响而改变了自然状况下分布的又一佐证。

2.2 不同类型农业土壤样品间Cd含量分析

土壤类型、土壤质地和农业土地利用对土壤中Cd的空间分布都有一定的影响[22]。其中土壤类型包含普通潮土、褐潮土、普通褐土、潮褐土4种分类；土壤质地分为轻壤质、中壤质、砂壤质；农业土地利用按照菜地、果园、林地、设施农业用地和粮田分类。

由于样点数据不符合正态分布，因此本文采用非参数检验方法（Kruskal-Wallis）比较不同类别土壤样品间Cd含量是否存在差异，检验结果如表5。结果显示Cd在不同的农业土地利用中存在显著差异，通过图2直观地表示出这种差异，整体上看不同农业土地利用类型中土壤样品间Cd含量大致为：菜地＞设施农业用地＞果园＞大田＞林地，且菜地含量显著高于林地，这主要是因为菜地利用率极高，常年摄入较多的化肥、农药；其次通过图3可知，Cd在不同的土壤分类中亦存在显著差异，普通褐土＞普通潮土＞褐潮土＞潮褐土，表明土壤类型可能是引起Cd空间变异的因素之一；Cd在不同的土壤质地中差异不显著。Zheng等[23]的研究结果提到Cd空间分布特征受人为因素和自然成因控制，与本研究结果相一致。

表4 土壤Cd含量的统计特征Table 4 Statistical characteristics of soil Cd content

2.3 土壤Cd空间结构分析

在土壤属性的多尺度分析中，要充分考虑变量的空间结构特征，即采样点的空间自相关性和空间变异等级结构[24]。土壤重金属Cd半方差函数参数见表6。结合采样点数据采用GS+软件进行半方差函数模型拟合，多次拟合统计结果比较，采用椭球模型误差较小。块金值表示空间变量受随机因素引起的变异，基台值是空间变量总的变异。基底比C0/（C0+C）表示变量随机部分引起的变异与变量总变异的比例。从结果可知，采样点基底比为49.8%，具有中等空间自相关性，可进行空间变异分析，同时也说明随机因素在总的系统变异中占较大的比重。

表5 不同分组土壤样品间差异检验结果Table 5 Results of ANOVA and Kruskal-Wallis test for the soil samples among different defined groups

图2 不同土地利用类型土壤样品箱线图Figure 2 The boxplot of soil sample by different land use

图3 不同土壤类型土壤样品箱线图Figure 3 The boxplot of soil sample by different soil types

概率累积曲线图可以判别变量分布特征，定性判断土壤重金属含量是否具有等级结构[24-25]。由土壤重金属Cd的对数概率累计分布图（图4）可知，Cd累计概率曲线存在明显的拐点，初步判断土壤中Cd来自两个特征明显不同的总体：A总体，含量水平较低，初步定位于自然背景来源；B总体，含量水平较高，可能来源于人类活动的排放。其次，半方差函数图也可以识别出土壤Cd的空间等级结构特征，Robertson等[26]通过半方差分析识别了农田土壤中pH值的空间变异等级结构。由于空间变量的变异可能存在巢式等级结构，因此该变量的半方差值随着距离的增加表现为台阶式上升的趋势，而半方差值突变的拐点则刻画了不同水平上的特征尺度[15]。观察半方差图（图5）可知，3.5 km处为半方差图突变转折点，说明土壤Cd的空间变异在3.5 km处发生了突变，表明土壤Cd空间变异存在等级结构。

2.4 多尺度空间变异识别

采用尺度方差法计算土壤Cd在不同等级水平上的方差，并识别土壤Cd的特征尺度。8次单独模拟结果（散点）及其平均值（折线），如图6所示。尺度方差随着尺度的增大而表征出不同的特征，即不同尺度上具有不同的空间异质性。随着尺度的增大，尺度方差分别在等级6（4a×2a）、等级9（16a×16a）处表现为波峰，暗示着土壤Cd空间变异的特征尺度位于特征尺度2.4～4.8 km、9.6 km左右。综合尺度方差的波峰和半方差函数在3.5 km处的拐点确定小尺度上特征尺度为3.5 km；大尺度上特征尺度为9.6 km，与Cd含量半方差函数的变程10.6 km相接近。

表6 半方差函数参数Table 6 Semivariance function parameters

图4 Cd含量的概率累计曲线Figure 4 Probability cumulative curve of Cd content

2.5 多尺度效应成因分析

土壤重金属空间异质性大小的尺度效应受控于不同尺度下控制土壤重金属变异的各种生态过程的重要程度，即影响Cd含量空间变异的环境因子具有不同的作用范围，实际采样工作中，可以根据取样的目的和关注的影响因素，选择接近影响因素作用范围的尺度作为采样尺度。比如，工业生产造成的点源污染通常影响的范围较小，交通运输污染源对土壤重金属Pb的影响区域一般在交通线两侧几十米的范围内，影响尺度比较小[27-28]，而土壤母质通常在大尺度上影响Cd的空间变异。也就是说，如果某一地区的主要污染来源于交通运输，那么道路两侧土壤重金属的空间变异范围（变程）是在污染源的影响距之内的，为了通过采样来全面地刻画重金属的空间变异特征，采样间距须小于影响距，这一点在柳云龙等[29]基于上海市3个区的采样研究中得到了证实。

图5 土壤Cd半方差图Figure 5 Semi-variance maps of Cd in soils

图6 Cd的尺度方差Figure 6 Scale variance of Cd

本文前述统计分析及半方差分析显示Cd存在中等强度空间变异且Cd空间变异存在等级结构，空间变异受到了土壤自然背景（土壤类型）等结构因子和人类活动（农业土地利用）等随机因子的共同影响。因此本文分别探讨土壤类型、土地利用强度对土壤Cd空间变异的影响范围。图7表明土壤类型、土地利用强度表现出一定的空间正相关性，说明2种影响因素在空间上的分布存在结构性。首先，观察土地利用强度的空间自相关系数由显著的空间正相关转为显著的空间负相关发生在3～4 km之间，在特征尺度2.4～4.8 km之间，因此推测土地利用在小尺度上影响了Cd含量的空间变异。图8为菜地、果园、设施农业用地空间分布与Cd插值图的叠加图，图中标识区域（画圈区）为Cd高值区，同时也是菜地、果园的集聚地。由图8可知：高值区范围较小，分布较破碎；菜地、果园集中分布区与Cd含量高值区吻合程度较好，这进一步佐证了土壤中Cd含量在小尺度上与土地利用之间的内在关系。Lv等[11]以日照市为研究区，采用因子克里金识别出土壤中Cd在小尺度上的空间变异受控于土地利用，与本文的研究结果一致。其一致源于两点：（1）我国市售的农药、化肥、有机肥料中Cd普遍存在超标现象[30]；（2）我国农业分散经营地块较小的实际情况造成土壤重金属空间异质性较强。其次，主要的3种土壤类型密度由显著的空间正相关转为显著的空间负相关发生在8～13 km之间，与特征尺度9.6 km以及土壤中Cd的相关距10.6 km相近，因此土壤类型可能是Cd大尺度上的主导因素。Nikos等[31]在杜埃布罗河流域土壤重金属污染源识别研究结果表明：在更大的空间尺度（流域尺度）上，自然因素最大限度地影响重金属的分布。但是，该研究结果显示自然因素的影响尺度达到130 km，研究区域涉及到的范围达到97.29 km2。因此，尺度分析不仅包含粒度层面，与幅度亦存在紧密的联系[13]，不同幅度的多尺度分析亦值得深究。土壤环境中Cd的来源比较广泛，比如交通运输污染、污水灌溉以及工业污染等，这些污染源均对土壤Cd的空间变异有影响。但是限于相关环境数据的获取难度比较大，本文选取的自变量数目并不能覆盖所有可能影响因素，在土壤中Cd空间变异的影响因素分析上还存在一定的局限性。因此，后续仍可进一步在此方面深入研究。其次，本文基于顺义区菜地、果园等零碎斑块大小的平均值确定0.36 km2，并以0.6 km为基本划区尺度，其目的是以较小的尺度划分防止小尺度空间变异的遗漏。但是，实际尺度划分操作缺乏严格数学推导，处理相对简单，因此，关于基本尺度划分已有进一步探讨的空间。