承德市土壤重金属空间结构与分布特征

安永龙，万利勤，李 霞，殷志强，卫晓峰，何泽新，贾凤超

(1.中国地质环境监测院，北京 100081；2.北京矿产地质研究院，北京 100012)

土壤是成土母岩经过风化、剥蚀、搬运、沉积等地质作用后叠加成土作用形成的自然物质体。在成土的每个阶段，都可能会受到外界环境(气候、海拔、地貌、人类活动等[1～2])的作用，因此在区域化土壤的物理化学性质和空间序列中形成了特有的空间变异性。传统的统计学方法往往忽略了土壤样本的固有属性数据和空间属性数据之间的交互式影响[3]，导致元素的异质性和空间依赖性无法很好地表达[4]。目前已有不少国内外学者针对不同区域时空尺度，利用地统计学结合GIS的方法解决土壤元素空间变异性问题。如崔萌等[5]运用地统计学和主成分分析法对北京市大桃主产区土壤5种重金属含量和空间结构特征进行了研究，结果显示Pb、As、Cr、Hg含量受结构性因素影响，而Cd含量受结构性因素和随机性因素共同影响；汪璇等[6]运用地统计学方法对三峡库区表层土壤部分微量元素进行了空间变异特征分析，结果表明7种微量元素为中等自相关性，2种为弱自相关性，且4个方向上变异程度不明显；Roger等[7]通过环境因子作为数据对土壤养分指数完成克里格插值，表明克里格插值与线性回归共同分析可以一定程度提高土壤养分的预测精度。

承德作为北京的“菜篮子”基地，土壤安全直接关系着首都人民的健康，在这样严峻的背景下摸清该地区土壤重金属本底及空间分布状况具有重要意义[8-9]。已有部分学者在区内按不同流域开展了土壤重金属调查研究工作，如孙厚云等[10]对承德市滦河流域12种重金属地球化学基线值进行厘定，发现表层土壤Cd、V、Ti等基线值高于河北省背景值，且Cd、Pb、Cu的累积程度较大；王梦雨等[11]对承德市柳河流域农田土壤重金属进行调查，发现Cd、Zn、Ni存在不同程度积累现象。本次研究基于地统计学和GIS技术，对承德市全域范围内表层土壤重金属As、Cd、Cr、Cu、Hg、Ni、Pb、Zn等的空间结构及分布特征进行解析，并对土壤重金属与Se之间的相关性进行了探索研究，旨在掌握承德地区表层土壤重金属空间分布状况及其特征，并进一步判断重金属与Se空间上是否存在一定关联，为承德地区着力打造健康安全的生态特色农果产业链，支撑“承德样板”[12]落地，开展国土空间规划和用途管控工作提供基础支撑。

1 研究区概况

承德市处于华北和东北两个地区的过渡地带，南邻北京，东靠辽宁，北接内蒙古高原，地理位置十分重要。地势北高南低，北部为坝上高原地貌，中南部为山地、丘陵地貌，地貌展布具有独特的地理优势，为特色农果产品和道地药材种植提供了有利的地域条件，杏仁、山楂、板栗、金莲花、黄芩等都是承德的名优特产。承德市河流分属三个流域四大水系，即滦河流域、三河流域、辽河流域，年产水量37.6×108m3，是京津唐的重要供水源地。区内属温带大陆性季风型半湿润山地气候，四季分明，雨热同期。年平均温度5.6 ℃，无霜期60～180 d，年降水量451～850 mm。土壤类型以褐土、棕壤、潮土、灰色森林土、粗骨土、栗钙土等为主。

承德矿产资源分布较多，以铁矿资源最为丰富，是我国仅次于攀枝花的第二大钒钛磁铁矿资源基地[13]。铁矿和金矿集中分布于双滦区和滦平县北部以及宽城满族自治县周边，银矿、钼矿、铜矿等分布较为分散。

图1 研究区土壤样点分布Fig.1 Distribution of soil sampling sites in the study area

2 材料与方法

2.1 样品采集与分析方法

本次研究利用了2017年承德地区土地质量地球化学调查数据，采样点分布情况见图1。土壤样品采集过程中采用GPS定位结合地形图定点方式，充分考虑土壤类型、土地利用方式。以土地利用类型为采样的基本单元，土壤采样深度为0～20 cm，并以“S”或“N”形采样，去除杂草、砾石、虫壳、粪便等杂物，一个组合样由3～5个子样点等组分混合而成，子样点与采样单元内采集的土壤类型一致。采集样品的原始重量不低于1 kg，过20目筛后样品重量不少于500 g。尽量选择在平稳耕地、农用菜地、林果地、山坡下侧土层较厚等地带采样，避开明显点状污染或表土已被破坏的干扰地段，保证了样品的典型性和代表性。

测试的指标为类重金属As，重金属元素Cd、Cr、Cu、Hg、Ni、Pb、Zn，微量元素Se。采用氢化物发生原子荧光仪测定As和Se，采用ICP-OES (PE，USA)测定Cd、Cr、Cu、Hg、Ni、Pb、Zn。分析测试中准确度和精密度采用国家一级土壤标准物质(GBW07349)控制，加10%空白样与平行样控制，且合格率符合规范要求，指标的加标回收率均符合国家标准。

2.2 数据处理

应用Excel 2016进行土壤数据基本参数统计分析，应用SPSS 19.0完成土壤数据正态分布性检验、相关系数分析、聚类分析、因子分析等。为了避免在利用统计学和地统计学方法时可能出现的比例效应[14]，所分析的数据必须全部服从正态分布或近似正态分布，否则需要将数据进行对数变换或者Box-Cox变换。聚类分析方法采用最近邻元素，选择区间Pearson相关性作为度量标准。因子分析采用最大四次方值法，抽取基于特征值大于1，最大收敛性迭代次数为25次。

目前国内外判断数据是否服从正态分布的方法较多[15]，主要与所使用的统计学软件类型(SAS、STATA、SPSS等)和样本量的多少有关。使用SPSS软件时，规则为探索性检验中基于理论假定的概率图、Q-Q图[16]、P-P图等不限定样本量；偏度峰度检验法适用于样本量N>200时；样本量72 000时，以Lilliefors检验(K-S检验[17])为准。基于处理软件和样本数量的限定，本文选用偏度峰度检验法和W检验法，结果表明，偏度值和峰度值远大于1，W检验的显著性均小于0.01；对9项元素进行对数变换后偏度值和峰度值基本在1附近，W检验的显著性均大于0.05，符合正态分布的要求，满足地统计学分析的假设条件，可以进行半变异函数的计算。

应用ArcGIS 10.2中Geostatistics analysis模块的普通克里格法结合半方差模型，完成最优内插过程。

应用GS+9.0软件完成9项元素半变异函数的计算和多种理论模型的拟合。半变异函数作为描述随机场和随机过程空间相关性的重要统计量，是分析土壤指标空间变异结构性和随机性的有效性方法[18]。该函数如下：

(1)

式中：γ(h)——半变异函数值；

N(h)——样点对的数目；

h——步长，即两分隔样点的距离；

Z(xi)——Z(xi)在空间位置xi上的土壤观测值；

Z(xi+h)——Z(xi)在空间位置xi+h上的土壤观测值。

在土壤地球化学工作中，通常区域化变量往往在不同的方向上会表现出不同的变异性，即便是在同一方向上也存在不同尺度的多维度的变异性。为了精准反映区域化变量的主要空间变异性，需要根据半变异函数的决定系数(R2)和残差(RSS)按照土壤元素所对应的最优理论模型进行半变异函数的拟合。本次所涉及的半方差函数模型主要有以下两种：

指数模型(Exponential)：

(2)

线性模型(Linear)：

γ(h)=C0+C·h/α，h>0

(3)

式中：C0——块金常数(间距为0时的半方差，通常由随机因素引起的变异)；

C——拱高(通常由系统因素引起的变异)；

α——变程(半方差达到基台值的样本间距，体现随机变量在空间上的自相关程度)。

3 结果分析

3.1 土壤重金属元素的统计分析

研究区内共采集土壤表层样品401件，表层土壤重金属元素是以统计学为基本依据，通过系统性计算各指标的最小值、最大值、中位数、平均值、标准差等统计特征值，来描述表层土壤重金属元素的含量区间和分布规律，见表1。表层土壤重金属元素的中位数与平均值相差较大，表明其中心趋向分布可能被异常值影响而使其呈非标准正态分布。将各指标平均值与河北省背景值[19]对比，发现土壤中Cd、Cu和Pb的平均含量均高于河北省背景值，其中以Cd最高，是河北省背景值的2倍；As、Cr、Hg、Ni和Zn的平均含量均与河北省背景值基本一致。

表1 土壤重金属的描述性统计

大量研究表明，变异系数(Cv)不仅可以反映总体土壤样本中各样点元素含量的平均变异程度，也可反映土壤元素的离散程度[20]。变异系数值越大，元素的空间分布越不均匀，离散程度越高。参照本次测试数据特点，将变异程度划分为三种类型：Cv<70%为均匀分布、70%≤Cv<100%为中强分异(中高起伏)、Cv≥100%为强分异(很大起伏)。依据该分类标准，Cd、Cu、Hg和Pb的变异系数分别为385%、143%、350%、118%，表明表层土壤中Cd、Cu、Hg和Pb属于强分异，其中以Cd和Hg的分异程度最强，表明局部地区可能受到人为活动影响，存在点源污染；As、Cr、Ni的变异系数分别为96%、93%、86%，属于中强分异；Zn的变异系数为64%，分异程度最弱。

在置信区间为95%内，不同土壤类型中单因素方差分析结果为F=2.913、p=0.006，判断土壤类型对土壤Zn含量均值具有显著性影响(表2)；不同土地利用类型中单因素方差分析结果表明，Cr的F=3.408、p=0.003，Cu的F=3.097、p=0.006，Ni的F=2.355、p=0.030，判断土地利用类型对土壤Cr、Cu、Ni含量均值具有显著性影响(表3)。

因此按照不同土壤类型，对Zn含量均值进行统计分析，发现潮土>栗钙土>粗骨土>褐土>棕壤>草甸土>灰色森林土>山地草甸土。按照不同土地利用类型，对Cr、Cu、Ni含量均值进行统计分析，发现Cu在内陆滩涂中含量高，Cr和Ni在建设用地中含量高，可能与人为活动影响有关。

表2 9种元素在不同土壤类型中的均值特征及单方差分析结果

表3 9种元素在不同土地利用类型的均值特征及单方差分析结果

3.2 土壤元素之间的相关性分析

对土壤重金属元素和Se进行相关性分析(spearman分析)，相关系数R越大，元素之间的关联程度越高；相关系数R越小，元素之间的关联程度越低或同源性越差。由各元素间的相关系数表可知(表4)，具有显著性相关关系的约占总关系数的25%，极显著相关关系约占总关系数的47.2%，其中Cr和Ni、Cd和Pb、Cd和Zn、Pb和Zn、Cu和Hg呈极显著的正相关，相关系数分别为0.832，0.698，0.739，0.731，0.635。Se与As、Hg、Pb、Zn呈极显著的正相关，与Cu和Ni呈显著的正相关，而与Cd和Cr之间的相关性较弱。Cr与As和Pb之间呈一定的负相关关系。由此推断Cr和Ni、Cd和Pb和Zn、Cu和Hg、Se和As之间可能具有相同的来源。

表4 9个土壤元素的相关性分析

3.3 土壤元素聚类及主成分分析

通过对研究区表层土壤中9种元素数据进行聚类分析，树状图结果表明，Cr和Ni最先聚合为一类，随后Cd、Pb、Zn聚合为第一类，Cu和Hg聚合为第二类，而As和Se聚合为第三类，如图2所示。

在一定程度上通常采用特征值来表征主成分影响力度大小，本次提取特征值大于1的主成分。由主成分分析结果可知(表5)，表层土壤中9种元素可被划分为四个主成分，第一个主成分(F1)的方差贡献率为34.50%，占四种主成分中较大比例，可以作为衡量该区重金属状况的一个综合性指标。第二个主成分(F2)的方差贡献率为20.55%，第三个主成分(F3)的方差贡献率为14.25%，第四个主成分(F4)的方差贡献率为12.58%，这四个成分的累积方差贡献率为81.87%，具备原始数据的大部分信息。

表5 土壤元素的主成分分析结果

为了更加精准确定各主成分所包含的元素信息，经过方差最大正交旋转后，对9种元素的四种主成分上因子载荷进行统计，如图3。结果发现，Cd、Pb、Zn在第一主成分中为高载荷，而在第二主成分、第三主成分、第四主成分中均为低载荷；Cr和Ni在第二主成分中为高载荷，而在第一主成分、第三主成分、第四主成分中均为低载荷；Cu和Hg在第三主成分中为高载荷，而在第一主成分、第二主成分、第四主成分中均为低载荷；As和Se在第四主成分中为高载荷，而在第一主成分、第二主成分、第三主成分中均为低载荷，见表6。

图3 表层土壤元素主成分载荷Fig.3 Factor loading analysis of elements in surficial soils

综上所述，传统统计学中主成分分析与相关性分析和聚类分析的结果在一定程度上保持一致，即Cr和Ni，Cd、Pb、Zn三者，Cu和Hg，As和Se来源可能相同。

表6 主成分分析旋转后的成分载荷矩阵

3.4 土壤元素的趋势分析

表层土壤一般容易受到外界因素(气候条件、地形地貌、人为活动等[21])的共同影响，因此区域性的表层土壤空间分布具有显著的趋势性和异向性。目前可通过Arcgis 10.2中趋势分析生成的三维趋势面曲线形状以及其变化的陡峭程度来判断空间数据的总体变化趋势。一般分为三种类型，即无(没有趋势效应)、一阶(区域化变量沿一定方向呈直线变化)和二阶(区域化变量沿一定方向呈多项式变化)[22]，通过异向性轴向自动搜索模块，对各元素数据进行处理，从而确定各向异性的趋势效应特征。以研究区As、Cu、Cd的空间分布趋势效应分析图为例，分别代表了无、一阶、二阶趋势效应(图4)，图中东西方向由x轴表示，南北方向由y轴表示，每个样点实测值的大小由z轴表示；正北方向投影面上浅绿色曲线表示东西向的全局性趋势效应变化，正东向投影面上紫色曲线表示南北向的全局性趋势效应变化。从图上可以看出，研究区土壤中Cd表现为在东西方向上先减小后增大的趋势，Cu表现为在东西方向和南北方向上呈直线分布，而As表现为无阶效应。

图4 土壤元素趋势分析图Fig.4 Trend analysis of soil elements

在Kringing插值过程中，选用不同的趋势类型会导致插值过程中产生的误差不同，因此需要综合分析各类参数特点，最大程度将各类误差降为最低，进而使半方差函数模型及其参数最适合。本次研究参照以下标准判断，平均误差(mean error，ME)的绝对值最接近0，均方根误差(root-mean-square error，RMSE)与平均标准误差(average standard error，ASE)应尽量相等且值越小越优，若前者小于后者，表明低估了预测值，反之则说明高估了预测值。标准化均方根误差(root-mean-square standardized error，RMSSE)值最接近1，如果标准化均方根误差小于1，说明高估了预测值，反之为低估预测值[23]，标准化平均误差(mean standardized error，MSE)的绝对值应最接近0。由表7可知，As、Cr、Pb、Zn、Ni、Se的无趋势预测与一阶趋势预测和二阶趋势预测相比，ME的绝对值更加接近0，ASE和RMSE最接近，MSE的绝对值更加接近0，RMSSE的值最接近1，因此Kringing插值时应选择无趋势更加精准，同理Hg和Cu选择一阶，Cd选择二阶。

3.5 土壤元素空间异质性特征

半变异函数(Semivariogram)是在区域化变量满足二阶平稳和本征假设的前提下，按照一定的抽样间隔获得样本方差的数学期望值，是地统计学中分析空间格局变量时最主要的应用工具。通常应用半变异函数及其参数可以对区域化变量的分布进行结构性和随机性的判断，进而结合克里格法完成对无样品空白区的插值和预测[24]。有三个相应的参数可直接反应空间变异性程度，块金值(C0)是在取样距离为0时，半变异函数为一个非0定值，反映的是最小抽样尺度以下由仪器测量误差等随机因素引起的微变异；变程(A0)也称半变异函数达到基台值时空间最大间隔距离，反映了变量空间自相关范围；块金系数C0/(C+C0)是块金值与基台值的比值，反应了空间异质性。而空间变异性程度主要由随机性变异和结构性变异大小所决定，结构性因素会增加变量的空间相关性，随机性因素则会降低变量的空间相关性，当块金系数小于25%时，表明区域变量的空间相关性强烈，此时空间变异主要由地形地貌、土壤类型、土地利用方式等结构性因素控制[25～26]；若块金系数大于75%，则表明区域变量几乎不具有或具有较弱的空间相关性，其空间变异程度受随机性因子控制比重较高；当块金系数介于25%～75%，表明区域变量的空间相关性强度居中，其空间变异由随机性因素与结构性因素共同控制[27]。土壤表层重金属元素空间变异性研究的关键是半变异函数模型的确定，模型拟合的效果可按照R2值最大和RSS值最小的原则检验[28]。

由表8可知，9种土壤元素的R2值除了Pb为0.176，其余均在0.334～0.773之间，说明上述重金属元素和Se的理论半变异与实验半变异模型较为吻合，拟合效果趋优。这9种元素中有8项元素的变程远远超过了8 000 m，最小的变程也在3 000 m以上，因此可对其进行Kriging插值。

As的半方差函数属于指数模型，块金系数为10.8%，变程为3 000 m，因此As具强烈的空间相关性，反映出研究区表层土壤中As主要受到地形地貌、成土母质、地质背景等结构性因素的影响，这与Chen等[29]和何厅厅等[30]结论保持一致，他们认为As的空间分布主控因素为土壤属性和地质背景，受人为活动影响较小；Cr、Cd、Hg、Pb、Zn、Cu、Ni、Se的块金系数分别为83.7%、83.0%、87.4%、100%、81.7%、89.6%、85.7%、86.9%，且半方差函数均属于线性模型(图5)。

表7 不同趋势阶数插值误差比较

表8 土壤元素含量变异函数理论模型及其相关参数

这8项元素具较弱的空间相关性，表明研究区表层土壤中这8项元素主要受到人为活动等随机性因素的影响，结论与Facchinelli等的[31]一致，其中Cu、Zn和Pb的来源与人类活动较为密切。不存在中等强度的空间变异性情况，这可能也与本次研究区边缘地带样品点较为稀疏有一定关系。

图5 土壤元素半方差函数图Fig.5 Semivariogram of the soil elements

3.6 普通克里格(Kriging)空间插值法

GIS在处理海量数据和图形视觉化表达方面具有很突出的特点，选取模型插值成图后可以直观地反映土壤元素的空间变异性。空间插值法能够使用一定量的样本量对区域范围内空白样本进行预测，从而分析整个区域范围的相关特征。空间插值的方法较多，主要分为确定性插值法(反距离加权插值法、径向基函数插值法、局部多项式插值等[32-33])和地统计学插值法(克里格法/协同克里格法、面插值法、经验贝叶斯克里金法等[34-35])。克里格法又称空间局部插值法，是在一定区域内对区域化变量进行无偏最优估计的一种方法，可分为泛克里格插值、简单克里格插值等，由于满足内蕴假设等条件[36]，本文采用地统计学中常用的普通克里格空间插值法[37]。

图6 土壤元素的空间分布Fig.6 Spatial distribution of soil elements

选取最优趋势参数和半方差函数模型进行普通克里格空间插值，得到研究区土壤元素的空间分布图(图6)。参照《土地质量地球化学评价规范》[38]和河北省土壤元素背景值调查[19]， As含量范围为2.9～29.9 mg/kg，整体含量值偏低；Ni含量范围为12.7～61.9 mg/kg，Pb含量范围为16.4～38.6 mg/kg，Se含量范围为0.055～0.401 mg/kg，整体含量值处于中等；Cd含量范围为0.082～0.742 mg/kg，Cr含量范围为28.0～228.4 mg/kg，Cu含量范围为11.3～142.4 mg/kg，Hg含量范围0.013～0.126 mg/kg，Zn含量范围为43.5～199.4 mg/kg，整体含量值偏高。从区域空间分布来看，承德市北部地区9种元素含量整体偏低，而南部地区含量整体偏高，体现出坝上高原地区整体受人为活动影响较低的特点。Cr和Ni空间分布具有很强的一致性，以丰宁县和隆化县为界，低值区主要分布于承德地区的北部，南部鹰手营子矿区和宽城满族自治县附近有低值区的零星分布；Cu和Hg保持了一定的空间分布一致性，高值区主要分布于承德南部地区，而Hg在丰宁县北部地区亦有高值区；Pb在中部形成了一条高值带，这又与Cd分布特征相似，这与区内铁矿、金矿等金属矿床的集中分布状况具有一定相关性。除去Pb中部的高值带，Zn和Pb二者在承德南部以及中部地区的高值区近似分布， Ni的高值区较好地表现出由南向北逐渐递变的过度特性，可能与区内金属矿产分布状况有一定关系；As低值区主要分布于承德中部隆化县—滦平县地区，围场满族蒙古族自治县以北、兴隆县以西地区。虽然上文统计学分析并未提及Se和Cd具有某种同源关系，但是二者高值区都分布在承德南部地区，西北部地区亦呈现局部高值，滦平县—双滦区、宽城满族自治县达到最高值，而西南部边缘地带和北部地区呈现低值区，也表现出弱空间异构性，这可能与Cd和Se都具有亲硫性和亲生物性，土壤中Cd、Se经常共生在一起有关，因此为富硒土地资源的安全开发利用方面带来了严峻挑战。

综上所述，单凭某一种方法去判断元素之间的同源关系是不准确的，需要在传统统计学结论的基础上进一步参考变异函数模型并叠加克里格插值结果进行优化。综合分析认为，研究区Cr和Ni、Cu和Hg、Zn和Pb、Cd和Se之间具有很好的空间分布相关性，来源具有一致性。建议应在现有土地利用现状属性基础上，充分考虑土壤重金属元素空间分布格局要素，在有矿业活动且土壤呈现重金属高值地区(如滦平县—双滦区北部山区)进行周期性土壤监测。同时，进一步开展区内土壤Cd和Se两者定量化控制因素的厘定及模型机理的研究，从而实现农作物富硒降镉的目标。

4 结论

(1)传统统计学分析中Cd、Cu、Hg和Pb属于强分异，变异系数分别为385%、143%、350%、118%，As、Cr和Ni属于中强分异，变异系数分别为96%、93%、86%，Zn属于弱分异，变异系数为64%。Zn含量均值受土壤类型影响显著，Cr、Cu、Ni含量均值则受土地利用类型影响显著。结合相关性分析、聚类分析和主成分结果，可知Cr和Ni、Cd和Pb和Zn、Cu和Hg来源相同。

(2)通过异向性轴向自动搜索功能结合不同趋势阶数插值误差综合对比，确定了As、Cr、Pb、Zn、Ni、Se适宜选择无趋势，Hg和Cu适宜选择一阶，而Cd适宜选择二阶。

(3)地统计学分析结果表明，As采用指数模型可较好地拟合，Cd、Cr、Cu、Hg、Ni、Pb、Zn、Se采用线性模型可较好地拟合。其中As的块金效应小于25%，因此具强烈的空间相关性，受结构性因素影响较大；Cd、Cr、Cu、Hg、Ni、Pb、Zn、Se的块金效应大于75%，具较弱的空间相关性，受随机性因素影响较大。

(4)普通克里金插值图直观地反映了承德地区土壤重金属的空间结构特征。研究区土壤的9种重金属元素都呈现明显的北低南高趋势，Cr和Ni、Cu和Hg、Zn和Pb、Se和Cd之间高值区分布具有很强的空间一致性，Hg在丰宁县北部地区亦有高值区，Zn和Pb高值区都分布于承德南部有向中部地区发展的趋势， Pb在中部地区形成了一条较宽的高值带，而其余分布特征与Cd相似，Ni的高值区较好地表现出由南向北逐渐递变的过度特性，在一定程度上较好地印证了地统计学所得出的规律。