元素空间分布插值方法的对比研究_以铜陵地区土壤中的重金属元素为例

1、第15卷第5期2008年9月地学前缘(中国地质大学(北京;北京大学Eart h Science Frontiers (China University of Geosciences ,Beijing ;Peking University 收稿日期:2008203230;修回日期:2008206230基金项目:安徽省科技攻关计划项目(08010302200;安徽省优秀青年科技基金项目(08040106907,04045063;中国地质调查局项目(20042012作者简介:袁峰(1971,男,博士,教授,矿床学专业。Email :yf_hfut 元素空间分布插值方法的对比研究:以铜陵地区土壤中的重

2、金属元素为例袁峰1,白晓宇1,周涛发1,李湘凌1,张鑫1,陈兴仁2,陈永宁2,陈富荣2,贾十军211合肥工业大学资源与环境工程学院,安徽合肥23000921安徽省地质调查院,安徽合肥230009Yuan Feng 1,Bai Xiaoyu 1,Zhou Taofa 1,Li Xiangling 1,Zhang Xin 1,Chen Xingren 2,Chen Yongning 2,Chen Furong 2,Jia Shijun 211S chool of Resources and Envi ronmental Engineering ,Hef ei Universit y of Tech

3、nolog y ,Hef ei 230009,China 21A nhui I nstit ute of Geological S urvey ,Hef ei 230009,ChinaYu an Feng ,B ai Xiaoyu,Zhou T aofa ,et al.Comparison betw een methods for interpolation of studying sp atial distribution of el 2ements :a case study of soil heavy metals in T ongling area ,south China.Ea rt

4、 h Science F rontiers ,2008,15(5:1032109Abstract :Taking six kinds of elements of soil contamination including As ,Cd ,Cu ,Pb ,Tl and Zn in Tongling area as an example ,we have made spatial interpolation for soil elements by using four types of representative interpolation methods ,including inverse

5、 distance weighted ,radial basis f unction ,ordinary Kriging and multi 2f ractal interpolation ,and have validated and appraised the results 1Comparison is made between methods with optimized parameters 1337soil samples in a total of 372samples were used for interpolation ,and the rest 35samples not

6、 having been interpolated were used for validation 1The results of comparison indicate that ordinary Kriging exhibits best effect in characterizing spatial distribution trend of soil elements ,yet the optimized mod 2els and parameters of semi 2variogram are still pending f urther study ;multif racta

7、l method is the best in charac 2terizing local abnormality and contamination ,whereas its f unction of reflecting the general texture of soil ele 2ment distribution needs in 2depth study ;inverse distance weighted and radial basis f unction interpolate less ac 2curately for the spatial distribution

8、of soil elements ,but are easy to use and to select optimized parameters.K ey w ords :spatial interpolation ;cross validation ;heavy metals ;soil ;Tongling area摘要:文中以铜陵地区As 、Cd 、Cu 、Pb 、Tl 、Zn 等6种土壤污染元素为例,选取常用且具有代表性的反距离加权法、径向基函数法、普通克里格法、多维分形法4种空间插值方法,进行土壤元素空间插值,并对其结果进行验证分析和评价。各方法均选取最优参数进行插值对比,土壤样本数共

9、372个,其中337个用于插值计算,35个不参与插值计算而用于验证插值结果。对比研究显示,普通克里格法对刻画区域土壤元素的空间分布趋势效果最佳,但其半变异函数模型及参数的优选仍有待进一步研究;多维分形法对刻画土壤元素局部异常和污染效果最佳,但其对土壤元素分布普遍特征的反映仍需深入研究;反距离加权法和径向基函数法对土壤元素分布的空间插值精度一般,但其简单易用、插值最优参数易于选择。关键词:空间插值;交叉验证;重金属元素;土壤;铜陵地区中图分类号:P595;X142文献标识码:A文章编号:100522321(2008050103207土壤中元素空间分布的研究常利用空间插值方法,将离散的点数据转化为

10、连续的面数据,以达到研究区域内元素的空间分布趋势和圈出局部区域内元素异常的目的。由于估值的不确定性加之不同目标对估值结果的要求不同,使得任何插值方法都无法满足所有要求1。因此,空间插值方法的合理选择是土壤元素空间分析的重要环节。铜陵地区是我国重要的矿业经济区,因矿产资源开发利用和高地球化学背景导致的土壤重金属污染问题日益严重223。近年来,对铜陵地区的区域土壤重金属元素空间结构、局部异常以及土壤元素环境污染评价与治理等方面的研究逐渐开展,积累了较丰富的数据和资料3211。铜陵地区土壤元素空间插值方法的对比、评价与筛选研究,对确定该地区土壤污染元素的空间分布,进行土壤元素污染的评价、预测、治理具

11、有重要意义,同时也将为其他地区的土壤污染元素研究提供有价值的参考。本文选取常用且具有代表性的反距离加权、径向基函数、普通克里格法及多维分形4种空间插值方法,以铜陵地区土壤中As、Cd、Cu、Pb、Tl、Zn等6种元素为例进行空间插值,并对其结果进行验证分析和评价。1数据来源及插值方法111数据来源铜陵地区土壤元素含量数据来源于安徽省地质调查院“安徽省江淮流域多目标区域地球化学调查”表层土壤样品分析数据,数据网度为1个数据/4 km2(图1。样点共372个,其中337个用于插值计算,35个用于插值结果的验证(不参与插值计算。其中As数据的变化区间为(61261016×10-6,均值为1

12、811×10-6;Cd为(0111151866×10-6,均值为01552×10-6;Cu为(1417169211×10-6,均值为6017×10-6;Pb为(201751919×10-6,均值为5011×10-6;Tl为(01472178×10-6,均值为0176×10-6;Zn为(451880813×10-6,均值为11317×10-6。图1研究区域及土壤采样点分布图Fig.1Study area and soil sample sites112插值方法11211反距离加权法反距离

13、加权法(inverse distance weighted,以下简称IDW是最常用的加权平均插值法之一。它考虑了空间数据值与距离成反相关的空间特性,其公式为12Z3(x0=Ni=1i Z(x i(1式(1中,Z3(x0为x0点处的估计值;N为用于插值的邻近样点的个数;Z(x i为样点x i处的实测值;i为第i个样点对估值点贡献的权重,它由下式计算:i=d i0-pNi=1d i0-p,Ni=1i=1(2式(2中,d i0为待估点与样点间的距离;p为距离的幂,它控制随着距离增大权系数的变化,其选择标准是最小平均误差,同时p也决定着内插结果的平滑效果。11212径向基函数法径向基函数法(radia

14、l basis f unction,以下简称RB F通过基函数计算待估节点的一组权系数,从而实现平滑插值。常见的5种基函数为:反转多重二次曲面(inverse multiquadric,多重二次曲面(mul2tiquadratic ,张力样条(spline wit h tension ,完全规则样条(completely regularized spline ,薄板样条(t hin 2plate spline ;它们类似于克里格插值法中的变异函数模型,通过调整基函数中的平滑因子可以控制插值面的光滑程度及估值精度13217。11213普通克里格法普通克里格法(ordinary Kriging ,

15、以下简称Kriging 是地统计学中的一种重要插值方法,在满足固有假设的条件下,其估计公式与IDW 法的公式相同,所不同的是其权系数i 的确定使用了半变异函数18。变异函数的计算和拟合是空间结构分析的基础,它反映了区域化变量的空间自相关性。当区域化变量满足二阶平稳性假设和本征假设时,半变异函数可用下式计算:(h =12N (h N (h i =1Z (x i-Z (x i +h 2(3式(3中,(h 为半变异函数;h 为滞后距离或步长;N (h 为距离等于h 的样点对数;Z (x i 和Z (x i +h 分别为区域化变量Z (x 在位置x i 和x i +h处的实测值。半变异函数可以用理论模

16、型进行拟合,如:Nugget 模型、线性模型、球形模型、指数模型、高斯模型等,所拟合变异函数模型反映了区域化变量的结构性(如方向性、空间自相关距离等和随机性(如随机误差等,同时也为普通克里格法的插值提供了有用信息19。良好的半变异函数是普通克里格法精确插值的保障20。为提高插值精度,在变异函数拟合过程中,应优先满足近距离的点的拟合效果,经过反复对比筛选得出最佳模型21。11214多维分形插值法Cheng 提出一种多维分形插值方法22(multi 2f ractal interpolation ,以下简称M FI ,将变量的空间自相关性与奇异性相结合,可以保持和增强数据的局部结构信息。该方法的实

17、质是在滑动平均加权的基础上加入了奇异性因子,其公式为Z (x 0=(x 0-2(x 0,(x -x 0Z (x (4式(4中,表示分形尺度,(x 0表示x 0处的奇异性指数,为滑动平均加权的权系数。它具有这样的特点:如果是背景场和非奇异场时,(x 0=2,那么通过该方法所计算的加权平均值与通常的加权平均相同;而当处于含量富集地段且局部场具有奇异性时,(x 0<2,该方法所得的结果将高于通常的加权平均结果;相反地,当处于贫化地段时,(x 0>2,该方法所得结果将低于通常的加权平均结果。该方法在矿产勘察、土壤、遥感图像处理等领域已有相关应用23226。本文所使用的多维分形法采用IDW

18、法的滑动平均加权公式。113插值计算实现方法本文使用Arc GIS 912的Geostatistical Analy 2sis Extensio n 模块进行IDW 、RBF 、Kriging 插值及交叉验证;使用Matlab 710编写MFI 法的奇异性因子算法程序,并结合surfer 810实现MFI 法插值。所有元素插值后的等值线填充图在Arc GIS 912下完成。2结果与讨论211最优参数的确定插值方法在参数选取不同的情况下会产生不同的结果,故插值方法的比较应建立在方法参数最优的前提下。IDW 法和RB F 法的参数较少,易于挑选最优参数。本文针对铜陵地区土壤元素含量数据,采用交叉验

19、证方法进行检验和优选,确定了IDW 法和RB F 法对于不同元素的最优参数(表1。Kriging 法的参数较为复杂,主要体现在变异函数的拟合上,通常需根据经验反复尝试、验证,不断优化;针对铜陵地区土壤元素含量数据,得到As 、Cd 、Zn 元素采用球形模型拟合最好,Cu 、Pb 、Tl 元素以指数模型拟合最佳(表1。M FI 法中滑动加权平均系数与IDW 法的权系数相同。212插值结果的验证与分析对比35个验证数据,插值方法的验证结果见表2,各方法得到的土壤元素空间分布见图2。考察6种元素的平均误差和均方差可以看出,IDW 法、RB F 法对6种元素估值的绝对平均误差较小(IDW 平均误差范围

20、为-7109801210,RB F 为-3151101059。这可能主要得益于Arc GIS 中对这两种方法参数的估计算法中考虑了交叉验证的平均误差。但这两种方法对As 、Cd 、Zn 元素的估值均方差相对较大(IDW 法的分别为10121278、01157、21621112,RB F 法的分别为11081767、01162、22181502,这表明平均误差小所反映出的只是正负误差的抵消结果,其多数点的估值误差较大。其中IDW 法所作的Pb 、Zn 、Tl 等值线图的平滑性相表1各插值方法的最佳参数Table1Optimum parameters in interpolation method

21、s块金值9136E-043195E-017165E-042189E-055191E-023193E-05基台值1168E-035130E-011188E-035173E-053176E-024161E-05a:CRS代表完全规则样条(completely regularized spline;b:Box2Cox变换幂。表2铜陵地区土壤元素四种插值方法的验证结果Table2Validation of four methods for interpolation元素插值方法平均误差均方差平均相对误差25275a10b90cAsCdCuPbTlZna:上下四分位点间的平均误差;b:小于第10百分位点

22、的平均误差;c:大于第90百分位点的平均误差。对较差,主要由于所算得的3种元素权系数偏高,而导致插值结果受局部临近的样本影响较大,而欠缺整体的平滑趋势。Kriging法的平均误差并不理想(-19132801002,但对As、Cd、Zn元素插值的均方差尚可(分别为9581977、01114、20601352。而M FI法虽在As、Cd、Cu元素的插值结果中表现出较好的平均误差和均方差(平均误差分别为-11686、01058、-11636,均方差分别为211769、01052、26631122,但对Pb、Tl、Zn元素估值误差较大(平均误差分别为21943、-01015、71512,均方差分别为12231954、01075、30921528。这些特征反映出,4种方法的插值结果没有绝对的最优者,对不同的元素会表现出不同的插值效果。这可能与方法对数据特性(如统计特征和空间分布特征的适应性有关。在对插值的精度要求不高的情况下,IDW 法和RB F法由于其简单性和易操作性是一种较好的选择。相对误差反映误差相对于实测值的大小,对于图2