多种空间插值法在土地质量地球化学评价中对比分析——以四川省叙永县麻城镇和摩尼镇为例

徐克全 ， 金立新， 郑勇军， 张 华， 包雨函 ， 时章亮

(1.四川省地质调查院，成都 610081；2.稀有稀土战略资源评价与利用四川省重点实验室，成都 610081)

0 引言

土壤元素含量在空间上具有自相关性，也存在多元素间交叉相关[1-2]，在空间分布上总体受到地质背景的控制，局部可能会受到人为因素的影响。为了获取特定空间范围数据，并对特定空间对象进行评价，需要对有限的空间数据进行插值处理。常用到反距离加权(IDW)、半变异函数的克里金(Kriging)等方法进行插值。但不同方法选取，或半变异函数中的块金、变程、偏基台值的设置，都将对插值结果造成影响，最终影响土壤质量地球化学综合评价面积统计的准确性，不利于土地风险管控和科学利用。笔者通过R语言(运用到sp、gstat、geoR、ggplot2、raster、intensity.analysis等R语言包)，使用麻城镇和摩尼镇1:10 000土地质量地球化学调查评价结果的表层土壤数据[3]，考虑元素空间自相关性可能受到地质背景、土地利用类型、元素的地球化学亲和性、空间距离的影响，评价单元对应赋值方法分别采用地层属性分层克里金法插值、土地利用类型分层克里金法插值、多变量协同克里金法插值、反距离加权平均插值法；考虑“有样”评价单元赋原值以及“无样”评价单元插值方法的更合理性，也是主流方法，增加对“有样”单元均采用局部加权平均(权重为公式3)赋值，而对“无样”单元插值分别采用“地层属性分层克里金插值”、“土地利用类型分层克里金插值”、“多变量协同克里金插值”和“反距离加权平均法(权重公式(2))插值”，共8种赋值方法进行养分单指标和土壤质量地球化学综合评价。借鉴土地利用强度分析方法，对以上8种赋值结果进行对比分析，探讨各种方法对评价结果的差异。

1 研究区概况

研究区位于四川省叙永县东南部麻城镇与摩尼镇境内，面积为54.59 km2，其中麻城镇于2015年撤乡设镇(图1)。共采集1 791(含重复样40件)件样品，基本密度为32.07件/km2；重复样数为40件，占2.29%；按网格化方法进行采样，采样单元以0.25 km2为采样大格，250 m×250 m采样小格，以采样小格作为评价单元。

图1 研究区地理位置图[3]

出露地层复杂多样，主要有寒武系(ε)、奥陶系(O)、志留系(S)、二叠系(P)和第四系(Q)(图2)。

图2 研究区地质图[3]

1)寒武系。娄山关组(∈2-O1l)，主要岩性为浅灰、深灰色白云岩、泥云岩，顶部夹页岩、含磷介壳灰岩、磷块岩。

2)奥陶系。红花园组(Oh)：灰、深灰色生物碎屑灰岩、鲕粒灰岩；湄潭组(Om)：上部灰黄色长石石英砂岩、砂岩，中部灰色介壳灰岩、石英砂岩，下部灰色页岩夹生物屑灰岩透镜体；艾家山组(Oa)：上部黄灰色页状泥灰岩、泥质灰岩；中部鲕粒灰岩，下部泥质生物屑灰岩；宝塔组(Ob)：灰、深灰色泥质网纹状灰岩，顶部为瘤状灰岩、泥灰岩；龙马溪组(O3-S1l)：深灰、黑色笔石页岩夹介壳泥灰岩。

3)志留系。桥沟组(Sq)：深灰色薄层状泥灰岩与泥、页岩不等厚互层；石牛栏组(Ss)：上部生物碎屑灰岩，下部瘤状泥灰岩夹珊瑚层礁；韩家店组(Sh)：灰、黄绿色页岩、粉砂岩夹生物碎屑灰岩及长石石英砂岩透镜体。

4)二叠系。梁山组(Pl)：摩尼一带上部深灰色、灰色泥质粉砂岩，下部炭质页岩、灰白色粘土岩，东部仅发育泥质粉砂岩；栖霞组(Pq)：灰、深灰色燧石结构生物屑灰岩，下部眼球状灰岩；茅口组(Pm)：浅灰、深灰色燧石条带生物碎屑灰岩，中上部夹灰黑色硅质岩，下部眼球。

5)第四系(Qesl)。残坡积层：棕红、黄褐色岩块、亚砂土、亚粘土。

研究区土地利用类型大部分为旱地、其次水田和果园(图3)。研究区土壤环境和养分指标背景值普遍高于全国土壤(A层)和成都经济区背景值(表1)。

图3 研究区土地利用图

表1 研究区土壤相关元素地球化学特征参数统计表[3]

2 空间插值方法与强度分析原理

2.1 反距离加权法(IDW)

在统计中，通常用一组数据的均值代表研究对象的真实值，但由于数据组中各数据对真实值的贡献大小不同，均值不能较好地反应真实值，因而增加了权重的概念，即加权平均法。该种方法增加了各数据对真实值的贡献大小(权重大小)后求平均值。

对于空间数据，离插值点越近的空间样本点数据越接近插值点，表现为离插值点越近的样本点所占权重越高。权重受到距离的影响，且距离越近权重越高，因而采用加权平均法对未知点插值，被称为反距离加权平均法。实现方法为以插值点与样本点间的距离为权重计算插值点S0处的加权平均值，见式(1)。

(1)

ω(si)=‖si-s0‖-p

(2)

(3)

在对研究区域内进行插值时，各样本点权重ω(si)选择式(2)进行计算，‖si-s0‖表示欧几里得距离，p为反距离加权幂，本文使用p为2。与克里金插值相比，反距离插值只考虑了样本和插值点的距离，而忽略了样本点的空间结构。

2.2 克里金插值法(Kriging)

空间数据往往存在自相关关系，且越相邻的越相似，含量会随着空间距离的变化呈现规律的变化。如成矿元素和伴生元素在原生晕或次生晕中存在规律的变化，靠近矿体元素含量高，远离矿体元素含量低；同时受到大地构造、地质背景等影响，区域上可能呈带状分布，而非距离越近含量越高，因而只考虑距离的反距离加权法不能较好地反应真实情况。克里金插值法不仅考虑了空间数据受距离的影响，也考虑了空间数据间的相关性结构，更真实地反应空间数据的分布情况。克里金插值法通过构建半变异函数模型来实现，方法如下：

(4)

图4 半变异函数拟合模型

2.3 强度分析原理

在GIS领域，用土地利用强度分析法(LULCC)分析各年度土地利用覆盖遥感数据的年度、土地类型和过度水平变化强度情况，是基于通过各年度土地利用栅格数据间交叉表，计算各年度间交叉表数据的年度水平、类型水平和过渡水平的土地转变强度来实现[4-6]。笔者将强度分析做了一些改动(修改了R语言包“intensity.analysis”中的部分函数算法)，基于各插值方法评价后的评价单元栅格数据间的交叉表(交叉表反应了两种插值方法由第一种插值方法转变为第二种插值方法后，前一种插值评价结果等级转变为后一种插值评价结果等级单元个数)，计算交叉表数据的插值方法间、评价等级间和等级增量来源情况的转变强度。式(5)为计算t类插值方法间的评价结果的差异强度，用St表示，其中t表示插值方法组tn-m(前者n类方法转变到后者m类方法)，Ctij代表t类插值方法间(例如，t1分层克里金插值法转变到t2反距离加权平均插值法)交叉表中，由评价等级i转变为等级j的评价单元(栅格)个数。式(6)为计算所有插值方法间的评价结果差异的平均强度，用U表示；式(7)为计算t类插值方法间中后者插值方法(tm)j类评价等级增加的平均比例，用Gtj表示；式(8)为计算t类插值方法间前者插值方法(tn)i类评价等级减少的平均比例，用Lti表示；式(9)为计算t类插值方法间i类评价等级增加情况，是来自于非i类评价等级转变的结果，用Rtin表示 ；式(10)为计算t类插值方法间n类评价等级，由其他等级转变而来的平均增加情况，用Wtn表示。

(5)

(6)

(7)

(8)

(9)

(10)

3 空间插值与土地质量地球化学评价

为了更好地反应评价单元的真实情况，应尽量保证评价单元内有样品数据，本研究调查的最小采样单元为250 m×250 m格子，因而将其作为本次的评价单元。评价单元共1 908个，其中“有样”单元1 794个，占比94.03%；“无样”单元114个，占比5.97%(图5)。

图5 研究区评价单元分布

图6 半变异函数值空间分布图

图7 各插值方法土壤N元素分级评价图

图8 各插值方法土壤P元素分级评价图

图9 各插值方法土壤K元素分级评价图

土壤质量地球化学综合评价涉及指标包括：养分元素指标N、P、K2O，环境元素指标As、 Cd、 Cr、 Cu、 Hg、 Ni、 Pb、 Zn以及土壤pH。现主流评价软件均采用原值赋原评价单元的方式，对“无样单元”采取插值或赋值的方式。为了对比不同插值方法，增加了单一插值的方式，此外使用主流方法时，对2个及以上样品数据的“有样”单元，采用了“局部反距离加权”的方式进行赋值。共涉及8种赋值方式，见method1-8。

method1、method2分别为按地层属性和按土地利用类型分层克里金插值法，原理为：除了考虑了各元素的半变异函数模型差异与各向异性外，还兼顾了不同分类区域存在不同的半变异函数(协方差结构)；method3为多变量协同克里金法，原理为：元素不仅自相关，也与其他变量存在协相关，因而该方法构建了自相关和协相关半变异函数模型，但需要设定这些变量具有相同的半变异函数模型和正定的偏基台值矩阵，所以采用各元素最优半变异函数模型参数的中位数，对评价单元网格进行插值[10]。

Method4为反距离加权平均法(权重为式(2))。method5至method8为主流赋值方法，其对“有样”单元均采用局部加权平均(权重为式(3))赋值，而对“无样”单元插值分别采用“地层属性分层克里金插值”、“土地利用类型分层克里金插值”、“多变量协同克里金插值”和“反距离加权平均法(权重式(2))插值”。

克里金法插值前，需构建半变异函数模型，其半变异函数值空间分布图(图6)可知： As半变异函数值在NW-SE方向上含量低于0.4，呈NW-SE向条带分布，各方向上存在明显的差异性；其余元素指标半变异函数值在各方向变化基本一致。因此，为了获得更准确地评价元素指标的克里金半变异函数模型，As较其余元素指标需要多考虑方向差异性。最终获得As在135°方向变程为3 310，在45°方向的变程为3 310×0.7=2 317，其余元素指标半变异函数模型参数见表2。method1、method2、method5、method6采用了不同的半变异函数模型参数并考虑了As元素方向差异性(模型参数见表2)；method3、method7使用的半变异函数模型参数为各元素最优偏基台和变程的中位数，分别为0.05和1 584，计算模型为球状模型“Sph”。

表2 半变异函数模型最优参数表

根据《土壤环境质量 农用地土壤污染风险管控标准(试行)》(GB15618-2018)和《土地质量地球化学评价规范》(DZ/T 0295-2016)规范，将最优模型插值后的数据，进行单指标养分元素N、P、K和土壤质量地球化学综合评价。评价结果见图7～图10，养分单指标和综合评价结果的空间分布主要受地质背景控制，各插值方法的评价结果，具有高度一致性。K元素分级结果很单一，显示研究区大面积富K元素。N和P元素分级结果差异最为明显，综合评价结果次之，表现出N元素丰富、P元素较缺乏和综合评价的优质土壤区域与志留系高桥组分布相一致；局部少数综合评价劣等土壤受人类活动影响，分布在研究区西南处摩尼镇附近。同时，表现出各指标不同插值法结果，局部存在明显差异，集中体现在研究区的中北部。虽然“无样”单元比率不足6%，但不同插值方法造成的评价结果面积的误差可能很大，最终影响评价结果的可信度以及政策制定。为了更好的进一步探讨不同插值方法间的差异，使用差异度分析方法。

图10 各插值方法土壤质量地球化学综合评价图

4 差异度分析

采用R语言“intensity.analysis”代码包，对各插值法对照组依次采用差异度分析、等级转变强度分析和增量强度来源分析。method1至method8在图中分别用Ⅰ至Ⅷ的罗马数值表示。

4.1 差异度分析

根据式(5)、式(6)得图11～图14。

图11 氮元素各插值方法差异度图

图12 磷元素各插值方法差异度图

图13 钾元素各插值方法差异度图

图14 土壤综合各插值方法差异度图

总体上钾元素不同插值法间差异率最低，普遍低于2%和图8吻合；土壤质量地球化学综合评价结果各插值法间的差异率较低，普遍低于5%；而单指标氮和磷元素除了method5～method8间，其余方法间差异率普遍高于10%，同图6和图7直观感受一致。

4种指标method1、method2、method3两两间差异值都在平均值附近，而method3与method4间差异度远高于平均值。这说明：克里金法和反距离加权平均法插值的评价结果差异很大，而克里金插值法的不同类型方法差异相对较小，但也存在明显差异。

method1与method5、method2与method6、method3与method7相比较的差异度远高于平均差异强度，存在强烈的差异。method4与method8相比较的差异度与平均差异度相当，但也存在明显的差异。这说明：单一用地层分层克里金法、土地利用分层克里金法或多变量协同克里金法，分别与不同克里金插值法+反距离加权平均法的差异度远高于平均差异强度。

4种指标的method5与method6、method6与method7、method7与method8的相比较的差异率不足1%，差异度远低于平均差异度，这说明“有样”单元赋反距离加权平均值与“无样”单元赋不同方法的插值结果，都能有效降低方法插值间造成的误差。

4.2 土壤质量地球化学综合等级转变强度分析

4种指标评价结果中，最受关注的为综合评价指标，为了探讨其等级转变方向，因而进一步讨论土壤质量地球化学综合等级转变强度。根据式(7)和式(8)所得各方法组土壤质量地球化学综合等级改变强度图(图15)。图15中红色条形表示前一种插值方法土地质量地球化学综合等级减少强度，绿色条形表示后一种插值方法土地质量地球化学综合等级增加强度；蓝色虚线表示前一种插值方法转变为后一种插值方法后所有等级增加的平均强度，虚线左侧dormant，表示转变强度不明显的土地质量地球化学综合等级，虚线右侧active表示转变强度明显的土地质量地球化学综合等级。

图15 各方法组土壤质量地球化学综合等级改变强度图

图16 各方法组各土壤质量地球化学综合等级增加来源强度图

图17 各方法组各土壤质量地球化学综合等级增加来源强度图

由图14可知,method5与method6、method6与method7、method7与method8相比较的土壤质量地球化学综合等级转变强度很低，百分比不到15%；method1与method2、method2与method3相比较的等级转变强度最高达到40%；method3～method4、method1～method5、method2～ method6、method3～ method7、method4～ method8的等级转变强度高达80%。

同时，反应出土壤质量地球化学综合评价的优质和良好等级在两方法对比的转变强度很大，都表现出了后者方法在优质和良好的等级较前者方法得到较强的增加。说明：插值方法的选择会较强地影响到土壤质量地球化学综合评价中的优质和良好土壤评价。

4.3 增量强度来源分析

根据式(9)和式(10)所得各方法组各土壤质量地球化学综合等级增量来源强度图(图16)。图16中每列表示增加的土壤质量地球化学综合等级，从左往右依次为优质、良好、中等和劣等，每行为相同方法组，红线表示当前等级增加的平均强度，红线上方target表示当前等级明显增加的来源等级，红线下方avoid表示当前等级不明显增加的来源等级。

method1与method5、method2与method6、method3与method7、method4与method8相比较，后者方法的优质增加明显来自于前者方法的良好，优质土壤增加强度较高达到10%～50%间，呈现了后者方法优质等级评价单元数量明显增加。说明前者插值方法部分评价单元数据被平均化，经重新加权平均赋值后，后者方法表现为优质土地增加。

Method1与method5相比较，后者方法劣等主要来源于前者的优级，也同样说明前者插值方法部分评价单元数据被平均化，经重新加权平均赋值后，后者方法表现为劣等土壤增加。

Method5与method6、method6与method7、method7与method8相比较的等级间转变强度很弱，增加强度普遍在2.5%以内，且基本与增加强度平均值相当，说明先插值再对有样品数据单元重新加权平均赋值后，各方法间土壤质量地球化学综合等级接近一致。

5 结论

1)几种插值后的评价结果都能反映其空间分布受到地质背景控制的影响。

2)通过借鉴土地利用强度分析法(LULCC)得知：①差异度分析,“有样”单元赋反距离加权平均值与“无样”单元赋不同方法的插值结果，都能有效降低方法插值间造成的误差；②等级转变强度分析,方法间差异度较大的原因在于优质、良好和劣等间差异度较大，其等级内差异度普遍在20%～80%间，插值法+反距离加权平均结合法其优良和劣等的增加强度明显，单一插值法其优良减少强度明显；③等级增量来源强度分析,插值法+反距离加权平均结合法其优质、劣等的增加量，主要分别源于单一插值法的良好和优质；总体呈现单一插值法靠中间的等级减少向插值法+反距离加权平均法两端等级增加的方向转变的趋势。

3)“有样”单元赋反距离加权平均值与“无样”单元赋不同方法的插值结果相结合，能避免因空间插值计算范围过大，促使高值或低值数据被平均后，导致较小的评价单元评价等级降低或升高的情况发生。

4)笔者只评价了250 m×250 m的栅格单元，有必要进一步在不同尺度栅格(Raster)和多边形(Polgon)图斑的评价单元上探讨，以便使不同尺度调查评价结果更趋于真实地反映客观情况。

致谢

感谢单位同事提供数据支持。感谢北京师范大学，韩春萌博士，对本人ARCGIS的使用提供帮助。