抚顺土壤重金属空间分布预测制图

辽宁宏源测绘规划建设有限公司　田艳萍

一、前言

1.研究背景及意义

土壤是人类赖以生存的最基本的自然资源之一。土壤重金属的污染能够影响土壤的正常功能,造成食物链的污染,并通过在生物体内的富集对人类和环境产生威胁。因此,对土壤重金属的空间分布研究,对减少重金属污染,提高人们的生活质量具有重要的意义。鉴于重金属Pb是一种对人体许多系统都有毒性的重金属元素,本文以ArcGIS软件为平台,应用地理统计学方法对抚顺市抚顺县的土壤含Pb量进行了空间分析,旨在通过对多种空间插值方法的研究，找出最优的空间插值模型，为了解该地区土壤污染空间分布情况提供科学依据。研究中应用较多的插值方法有反距离加权法、克里格插值法、局部多项式插值法等。

土壤重金属污染空间分布研究中，由于成本及采样方式的限制很难做到高密度采样，因此，在操作中要根据有限的离散点绘制区域土壤重金属空间分布图。插值模型的选择是土壤重金属空间分布预测的基础，但是研究发现，不同的插值方法在预测精度上存在着较大的误差，因此通过研究选出最优的插值模型在土壤重金属污染研究中至关重要。本文通过不同插值方法对抚顺县土壤重金属Pb空间分布情况进行预测，选出最适合该地区的插值方法，做成土壤重金属Pb空间分布预测图，有助于揭示区域内重金属Pb空间变异的趋势及其影响因素，对分析重金属污染物的来源与扩散、进行土壤环境区规划、制定土壤环境保护策略及进行重金属污染风险评估具有重要意义。

2.国内外研究进展

近年来国内外科学家在利用地统计学对区域土壤重金属空间分布进行研究上，做了大量的工作。Ersoy等运用地统计学方法对英国carsington牧场的土壤重金属Pb、Zn、Cu等重金属进行了空间分析，指出历史上的该地区矿业开采活动造成了4种重金属在该牧场土壤长期污染，并对控制该地区重金属污染风险和污染土壤修复提出了建议。saby等对法国巴黎附近的土壤Pb含量进行了大尺度研究，并运用协同克里格、普通克里格和对数正态克里格法分别进行了空间插值，比较发现对数正态克里格法插值结果最为理想。

邢怀学等人以合肥地区典型中心城镇—大兴镇为例,选取了土壤中Cd、Hg、Pb、As4种重金属元素，统计分析了4种重金属元素含量的统计学特征在地理信息系统（GlS）的支持下，分别采用反距离加权法（IDW）、径向基函数法（RBF）、普通克里格法（OK）三种常用的插值方法进行空间插值，为验证这三种插值方法对插值结果的影响，研究采用交叉验证法，即分别假设每一采样点的要素值未知，值与估算值的误差大小评判插值方法的优劣。研究采用交叉验证法来验证其插值用周围采样点的值来估算，然后根据所有采样点实际观测值与估测值的误差来评判优劣，结果表明三种方法的结果都基本上反映了土壤重金属元素的空间变异规律，而普通克里格法的插值效果最好，精度较高且能反映空间变异方向，反距离加权法和径向基函数法次之。

二、技术路线

三、研究方法与理论

1.研究区概况

抚顺县地处东经123°40′至124°27′，北纬41°27′至42°04′之间。东接新宾满族自治县，南连本溪市溪湖区和本溪满族自治县，西靠沈阳市东陵区、苏家屯区、抚顺市区、顺城区、北邻铁岭县，东北与清原满族自治县接壤。全县下辖4个镇，8个乡。面积约2200平方公里，乡村人口约20万人。抚顺县地区处于低山丘陵向平原过渡地带。东部、东南、东北地势高峻，而西部、西南、西北稍平缓。整个地势由东向西缓倾，中部为浑河谷地。平均海拔100—300米，境内山峦起伏，森林茂密，河流纵横。地貌特征为“七山一水半分田，半分道路和庄园”。抚顺县处在南温带亚湿润区内，属大陆性季风气候。雨热同季，四季分明。年平均气温7.8℃。一月平均气温-14℃，最低气温-35.2℃；七月平均气温24.5℃，最高气温35.8℃。年平均降水量800毫米，多集中在七、八、九月份，无霜期145天左右。境内有沈抚线、沈吉线两条铁路贯穿北部；省级干线公路沈环线、沈通线越县境南部；国级干线公路黑大线横跨北部。县乡公路发达，半数以上乡镇通上了柏油路。

2.采样点分析

本文研究所选样点共40个，样点随机采样，基本均匀分布在抚顺县各乡镇。

3.插值方法

（1）反距离加权插值

反距离加权法是以插值点与样本点之间的距离为权重的插值方法，插值点越近的样本点赋予的权重越大，其权重贡献与距离成反比。可表示为

其中Z是插值点估计值；Zi（i=1^n）是实测样本值；n为参与计算的实测样本数；Di为插值点与第i个站点间的距离；p是距离的幂，它显著影响插值的结果，它可以通过均方根误差的最小值确定其最佳值。均方根预测误差是一种通过交叉验证计算得到的统计量。权重与预测点和已知样点间距离的p次幂成反比，因此，随着距离的增加，权重迅速减小；权重减小的速度取决于p值大小。通常p在1—3之间取值，大多数情况下p=2。当取p=2时，即称作反距离平方加权法。

（2）局部多项式插值

局部多项式是将全局多项式方法和移动平均过程结合起来的一种内插方法。局部多项式需要应用最小二乘多项式拟合数据，通常选择一次、二次或者三次多项式。局部多项式只用于拟合一个窗口定义的局部样点区域，通过多个多项式作为局部方程式来拟合研究区域，每个多项式都处在特定重叠的邻近区域内，根据已知样点拟合的多项式表达式和插值点的坐标，就可以计算插值点的变量值。

（3）克里格插值

克里格插值（Kriging）又称空间局部插值法，是以变异函数理论和结构分析为基础，在有限区域内对区域化变量进行无偏最优估计的一种方法，是地统计学的主要内容之一。南非矿产工程师D.R.Krige（1951年）在寻找金矿时首次运用这种方法，法国著名统计学家G.Matheron随后将该方法理论化、系统化，并命名为Kriging，即克里格方法。

克里格方法的适用范围为区域化变量存在空间相关性，即如果变异函数和结构分析的结果表明区域化变量存在空间相关性，则可以利用克里格方法进行内插或外推；否则反之。其实就是利用区域化变量的原始数据和变异函数的结构特点，对未知样点进行线性无偏、最优估计。无偏是指偏差的数学期望为0，最优是指估计值与实际值之差的平方和最小。也就是说，克里格方法是根据未知样点有限邻域内的若干已知样本点数据，在考虑了样本点的形状、大小和空间方位，与未知样点的相互空间位置关系，以及变异函数提供的结构信息之后，对未知样点进行的一种线性无偏最优估计。

本文选用了普通克里格插值方法对土壤重金属Pb空间分布情况进行预测。

普通克里格法（OK）是地统计学中的一种重要插值方法，在满足固有假设的条件下，其估计公式：

其中，Z*（x0）为估计值，Z（xi）为位于区域xi位置的观测值，λi为权重，权重之和等于1。权重的确定使用了半变异函数。变异函数的计算和拟合是空间结构分析的基础，它反映了区域化变量的空间自相关性。良好的半变异函数是普通克里格法精确插值的保障。当区域化变量满足二阶平稳性假设和本征假设时，半变异函数可用以下公式计算：

其中，γ（h）为半变异函数；h为滞后距离或步长；N（h）为距离等于h的样点对数；Z（xi）和Z（xi+h）分别为区域变化量Z（x）在位置xi和（xi+h）处的实测值。若以h为横坐标，以γ（h）为纵坐标绘制函数曲线图，称为半方差函数曲线图，即：

半变异函数图

此图可直接展示Z（xi）的空间变异特点。

在半变异曲线图中有两个非常重要的点：间隔为0时的点和半变异函数趋近平稳时的拐点，由这两个点产生四个相应的参数：块金值（Nugget）、变程（Range）、基台值（Sill）和偏基台值（Partial Sill）它们的含义表示如下：

块金值（Nugget）：理论上，当采样点间的距离为0时，半变异函数值应为0，但由于存在测量误差和空间变异，使得两采样点非常接近时，它们的半变异函数值不为0，即存在块金值。测量误差是仪器内在误差引起的，空间变异是自然现象在一定空间范围内的变化。它们任意一方或两者共同作用产生了块金值。

基台值（Sill）：当采样点间的距离h增大时，半变异函数r（h）从初始的块金值达到一个相对稳定的常数时，该常数值称为基台值。当半方差函数值超过基台值时，即函数值不随采样点间隔距离而改变时，空间相关性不存在。

偏基台值（Partial Sill）：基台值与块金值的差值

变程（Range）：当半方差函数的取值由初始的块金值达到基台值时，采样点的间隔距离称为变程。变程表示了在某种观测尺度下，空间相关性的作用范围，其大小受观测尺度的限定。在变程范围内，样点间的距离越小，其相似性，即空间相关性越大。当h>R时，区域化变量Z（x）的空间相关性不存在，即当某点与已知点的距离大于变程时，该点数据不能用于内插或外推。

当限定的样本点间隔过小时，可能出现曲线图上所有r（h）≈Nugget，即曲线为一近似平行于横坐标的直线，此时半方差函数表现为纯块金效应。这是由于所限定的样本间隔内，点与点的变化很大，即各个样点是随机的，不具备空间相关性，区域内样点的平均值即是最佳估计值。此时只有增大样本间隔，才能反映出样本间的空间相关性。

空间相关性的强弱可由Partial_Sill/Sill来反映，该值越大，空间相关性越强。相应地，Nugget/Sill称为基底效应，表示样本间的变异特征，该值越大，表示样本间的变异更多得是由随机因素引起的。

半变异函数的模型包括球状、高斯、指数和线性等模型等。

4.插值精度检验

对于不同空间插值方法估计效果的检验，一般采用交叉验证法来检验插值效果。即先假定每一个采样点的含量值未知，利用周围样点的值来估算，然后计算估计值与实际测定值的误差，根据误差统计结果评估插值方法的优劣。常用的误差统计指标有平均误差（ME），均方差（MSE），均方根误差（RMSE）。ME越接近于0，插值误差越小；MSE的值越小，精度越高；RMSE的值越小，精度越高。ME总体反映估计误差的大小，MAE可以反映估计值可能的误差范围，RMSE可以反映利用样点数据的估值灵敏度和极值效应。用交叉验证法进行精度的验证时可以应用所有数据交叉验证也可以应用独立数据集验证。本文采用的是将全部的40个样点的数据既作为插值数据有作为检验数据的交叉验证法。

四、结果与分析

1.土壤Pb含量统计和空间结构特征

（1）统计特征分析。

对所有采样的数据进行统计处理结果如表1，从中得出研究区土壤重金属Pb含量最大值为47.52毫克/千克，最小值为2.43毫克/千克，平均数为17.20毫克/千克，中位数13.59毫克/千克。经过对数转换后，中位数和平均数接近，峰度接近3，可以近似看成正态分布。从数据直方图中也可以看出采样点土壤中Pb含量数据经过对数变换基本呈正态分布，如图1。由于地统计的半方差函数的计算要求数据符合正态或近似正态分布，否则可能产生比例效应，从而可能使实验方差函数产生畸变，抬高基台值和块金值，增加估计误差，甚至会掩盖其固有的结构，因此在本文中为了可以应用地统计插值对土壤重金属Pb空间分布情况进行预测，应先将数据进行对数转换。

表1　数据统计特征表

图1　数据直方图

（2）数据空间结构特征分析。

对土壤重金属含量的常规统计分析能够概括土壤重金属含量的整体特征，但不能够反映其局部的变化特征，即只在一定程度上反映样本全体，而不能定量地刻画土壤重金属含量的随机性和结构性、独立性和相关性。为了更好地了解土壤重金属含量的空间变异特征，必须进一步采用地统计方法对土壤重金属含量的空间变异结构进行分析和探讨。

在计算半方差时，若依据的是不规则格网数据，则应设置Lag步长和最大步长这两个参数。此最大步长亦为分离距离h的最大值。当点对之间的距离落在所设定的步长范围之内时，就以此范围之内的所有点对的距离的平均值作为分离距离，并以这些点作为一组来计算半方差。故最大步长一旦确定，增加步长就必然要减少分组数。本文中的样点数据并不是用规则网格法选取，因此需要设置合理的步长及分组数，步长过大短程的自相关性将被掩盖；步长过小，就会产生许多空的步长组，步长确定的基本原则是步长值和步长组数乘积为最大步长值，即样点间最远距离的一半。基于上述分析，本文中分离距离为34000米，步长设置为3000米。

分别选取了球状模型、高斯模型、指数模型进行半变异函数拟合，并根据数据采用不同类型的半方差模型进行拟合，根据残差值（RSS）越小越好（主要指标），拟合优度（R2）越接近于1越优的标准选取最适当的拟合模型，得到的结果及有关参数如下，见表2。

表2　 半变异函数拟合模型

表2中的C0表示块金值；C0+C表示基台值；a是变程

从表2可以看出它们的块金值C0与基台值（C+C0）之比［C0/（C+C0）］的顺序为指数模型（72%）>球状模型（56.1%）>高斯模型（50.8%），均小于75%。按照区域化变量空间相关性程度分级标准，当块金系数分别为≤25%、25%～75%、>75%时，分别提示变量空间自相关程度为强烈、中等及微弱，本文中重金属Pb空间自相关性为中等强度自相关。通过拟合精度检验结果来看，球状模型的残差值（RSS）最小，拟合优度（R2）最接近1，拟合效果最好。因此，在用普通克里格插值对研究区土壤重金属空间分布制图是应选用球状模型。

2.土壤重金属Pb的趋势效应分析

图2中，X轴表示正东方向，Y表示正北方向，Z轴表示个点测定值的大小；左后投影面上的线表示东西向全局性的趋势变化情况，有后投影面上的线表示南北向全局性的趋势效应变化情况。从图2中可以看出，抚顺县土壤重金属Pb的含量分布不具有很明显的趋势效应，东西方向基本上呈平直性或零趋势，南北方向轻度的呈二阶多项式变化。因此，直接从趋势图上不能准确地判断出趋势效应，需要进行各种趋势效应的比较。

图2　趋势分析图

3.土壤重金属Pb空间插值方法分析

（1）普通克里格插值（OK）。

由于不能判准确地判断出趋势效应，为了比较不同趋势效应的插值情况，在普通克里格方法下，假定半方差函数拟合模型分别为球状模型；趋势效应分别选0（无趋势）、一阶（线性）、二阶（二阶多项式）。通过不同模型插值结果精度的对比，选出最适合的插值模型。根据克里格输出预测表面所生成的交叉验证参数评定标准：ME、MSE越接近与0，精度越高，RMSE和ASE尽可能的接近，RMSSE越接近1的插值方法精度越高。从表3可以看出，平均误差ME的绝对值最小的是二阶趋势效应的内插，其次是0阶趋势；MSE绝对值最小的是0阶趋势效应的内插方法，其次是一阶趋势；RMSSE最接近于1以及ASE与RMSE最接近的均是采用二阶趋势效应的内插方法；RMSE最小的是0阶趋势效应的内插方法；综合比较各种误差的大小并结合空间预测图来看，二阶趋势效应的内插方法比较好。

表3　土壤重金属Pb采用不同趋势指数的OK法预测误差

因此，以OK法进行土壤重金属Pb含量空间插值时，二阶趋势的球状模型效果最好。

（2）反距离加权插值（IDW）。

反距离加权法是最常用的插值方法之一，在本文中IDW加权系数P取值1-3。利用交叉验证中不同参数P产生的均方根预测误差（RMSE）的大小，其值越小，参数P值越优的原则，通过对不同参数P计算出来的插值预测表面进行交叉验证的结果见表4。

表4　IDW交叉验证误差统计结果

从表4中可以看出研究区内土壤重金属Pb以一次反距离权重插值效果最好。

（3）局部多项式插值（LPI）。

与反距离加权法中的P值的选择标准相同，局部多项式插值也以均方根预测误差（RMSE）较小的模型较优，插值效果较好。通过对不同参数P计算出来的插值预测表面进行交叉验证的计算结果见表5。

从表5中可以看出研究区土壤重金属Pb以1次局部多项式插值效果相对较好。而回归系数为3时插值产生的误差较大。

表5　LPI交叉验证误差统计结果

4.不同插值模型的土壤重金属Pb含量空间插值精度比较

本文分别应用了OK、IDW及LPI三种插值方法对土壤重金属Pb含量空间分布进行分析。综合三种方法来看，IDW插值交叉检验精度误差最小，插值的精度在三种方法中最高，OK法次之，LPI方法误差最大精度最低。三种方法预测的土壤重金属Pb空间分布如图3，虽然IDW方法交叉验证误差小，但是在插值图上可以看出IDW方法虽然能总体反映，然而在局部现象的反映上没有OK方法好。

图3　不同插值方法插值图

因此，OK法在本次研究中用来进行土壤重金属插值更为合适，取得的效果更好。

5.土壤重金属Pb空间分布制图中分级色彩确定

所谓分级色彩表示方法就是将要素属性数值按照一定的分级方法分成若干级别之后，用不同的颜色来表示不同级别。每个级别用来表示数值的一个范围，从而可以明确反映制图要素得定量差异。色彩选择和分级方案是分级色彩表示方法中的重要环节，因为颜色的选择和分级的设置要取决于制图要素的特征，只有合理的配色方案和科学的分级方法才能将地图中要素的宏观分布规律体现得清晰明确。本文在研究土壤重金属Pb的空间分布特征时选用几何间隔法分级，这种方法将相似性最大的数据分在同一级，差异性最大的数据分在不同级，这种方法可以较好保持数据的统计特性。

五、土壤重金属空间分布及污染评价

1.土壤重金属Pb的空间分布特征

图4　研究区土壤重金属空间分布图

通过最优的插值模型得到研究区土壤重金属Pb空间分布从分布图（图4），可以看出，研究区域土壤重金属Pb分布以条带状为主斑块状为辅，含量高的地区集中在中部及东北部地区，含量33.57—47.52毫克/千克主要分布在这两个地区，并以此两个高值中心向四周递减，西部地区及南部地区土壤Pb含量较低。

2.土壤Pb自然背景值及污染评价

（1）土壤重金属Pb环境背景值

土壤环境背景值是指在不受或者很少受现代工业污染与破坏的情况下，土壤原来固有的化学组成和结构特征。但是现代工业和人类的活动的影响已遍布全球，很难找到绝对不受人类影响的土壤，所以自然背景值是在一定时空条件下的相对含量。本文引用的是“中国土壤背景值研究”中的辽宁省A层土壤背景值。

（2）单因子指数法

通过单因子评价，可以确定主要的重金属污染物及其危害程度。一般以污染指数来表示，以重金属含量实测值和环境背景值之比来计算污染指数。

公式中：Pi为重金属的污染指数；Ci为重金属实测值；Si土壤环境背景值。污染评价标准如表6。

表6　单项污染程度分级标准

3.土壤重金属Pb污染评价制图

图5　土壤重金属污染评价图

通过Arcgis以污染指数作为分级标准对得到的土壤重金属Pb空间分布图重新进行分级，得到污染评级图（图5），从而能直观的了解研究区内重金属Pb的污染情况。

从土壤重金属污染评价图中可以看出，研究区土壤重金属铅污染不严重，大部分面积为非污染土壤，没有重污染地区。轻度污染地区较多，轻污染地区主要集中在救兵乡东北部地区、上马乡西南部地区及汤图乡大部分地区。经调查发现研究区域矿产丰富，污染相对严重的几个乡镇分布有大大小小的洗煤厂，采矿场等污染企业。同时，救兵乡又是抚顺工业发展比较发达的乡镇，除了采矿场及洗煤厂外还有大量的木制品生成厂。

六、结论

土壤重金属污染已成为威胁环境安全及人类健康的要素，土壤的环境质量也影响地区的经济发展，采用GIS空间分析技术来研究土壤问题已成为发展趋势。本文应用地统计学与GIS相结合，以ArcGIS为平台研究抚顺县土壤重金属Pb空间分布问题。得到以下结论：

1.本文比较了三种插值方法的精度，得出就研究区域而言OK法插值效果优于IDW方法，LPI方法插值效果最差。

2.插值模型选择对空间分布插值结果的影响较大，就同一种插值方法而言，不同的模型参数设置，得到的插值精度也存在着较大的差异。本文应用多种插值模型对比，从中选出最优的插值模型——OK方法球二阶趋势球状模型，得到的重金属空间分布图与实际情况更为符合，对于了解当地的重金属污染分布有一定的参考意义。

3.半变异函数拟合的精度对OK方法的精度影响较高，在样本量较小的情况下，很难获得理想的半变异拟合函数。同时，半变异函数的拟合存在较大的主观性，这些因素对本文的空间分布预测精度降低，要想提高精度要适当增加采样点的数量，而这样就增加了操作成本。因此合适的采样很重要。这也是本文有所欠缺需改进的部分。