GIS与地统计学的土壤水分空间插值方法

张 优, 王 娟, 张 杰, 彭文甫*, 朱亚兰

(1. 四川师范大学 地理与资源科学学院, 四川 成都 610066；2. 四川师范大学 西南土地资源评价与监测教育部重点实验室, 四川 成都 610066)

土壤是由多种因素相互作用而形成的,其属性在空间分布上有一定的差异[1].实现精确农业养分管理和解决全球变化等环境问题的关键在于准确掌握土壤理化性状的空间变异规律[2].土壤干湿状况对监测土壤退化、土地干旱有着巨大作用,是监测土地状况的重要指标[3].土壤水分的时空分布研究对农业干旱监测和生态环境监测具有重大意义.目前,揭示土壤属性空间分布特征的主要手段仍是通过野外采样和室内测定.但其问题在于,无论采样密度多大,均无法得到所有区域的土壤属性值.而空间插值技术可以将离散的数据点转化为连续的数据曲面从而实现研究区土壤属性的全覆盖.作为土壤属性变化的时空定量监测方法,土壤属性空间插值方法及其精度是数字土壤研究的重点.姚雪玲等[4]以黄土丘陵沟壑区洋圈沟小流域为例,介绍了基于ArcGIS与地统计学的回归克里格模型插值方法,详细介绍了基于虚拟变量技术将分类变量引入模型的过程,同时将回归克里格模型、回归线性模型和普通克里格模型的插值精度进行对比分析.史文娇等[5]介绍了土壤属性空间插值的常用方法与精度验证的方法和指标,并且总结了能够提高土壤属性插值精度的6种途径.陈思萱等[6]针对广东省某地区土壤中采样样品点的重金属污染物砷(As),选取统计学中具有代表性的3种空间插值方法,对比分析不同空间插值方法下土壤重金属浓度估算精度和污染格局空间区划差异.龙军等[1]根据不同地貌类型特征在福建省各地级市选取了9个典型县,利用2008年采集的样点数据评价不同空间插值方法对土壤有机质含量插值精度的影响.

本文以山地平原过渡地区四川省绵竹市为研究区,选择了几种常用的空间插值方法,包括反距离权重法(IDW)、径向基函数空间插值方法、普通克里格插值方法以及结合线性回归模型和普通克里格的回归克里格插值方法,探究不同插值方法对山地平原过渡地带土壤水分的影响.

1 研究区概况

绵竹市位于四川盆地西北部,位于山地平原过渡地区,地处30°09′—31°42′N,103°54′—104°20′E之间,面积1 245.3 km2,辖19镇2乡,总人口52万.本文研究区为绵竹市部分区域(图1).该区域属四川盆地中亚热带湿润气候区,气候温和,年均温15.7℃,年降水量为1 053.2 mm,雨量充沛,四季分明,大陆季风性气候显著.地貌形态区域分异十分明显,西北部为山地,东南部为平原,地势西北高,东南低,由西北向东南倾斜.森林占幅员面积49.8%.土地利用类型主要以林地、耕地、草地、建设用地、水体和未利用地6类为主.境内平原区土壤包括13个土属,8个亚类,土地肥沃,物产丰富.

图 1 研究区位置图

2 数据准备及预处理

研究涉及数据主要包括：绵竹市行政区划矢量图,来源于四川省测绘地理信息局；绵竹市2016年土地利用类型数据为Landsat 8 OLI解译数据,解译为耕地、林地、草地、水域、建设用地和未利用地共6类[7],经实地样点验证,一级土地利用解译精度满足后续分析需要(总体精度达94.3%)；绵竹市数字地面模型(DEM)来自USGS(https://www.usgs.gov/),空间分辨率为30 m×30 m.

绵竹市土壤含水量样点数据通过野外采集获取,考虑到样点空间分布对插值结果的影响,在样点选取时按以下思路进行.首先在研究区行政区划图上创建2 km×2 km的网格,以各个格网中心作为预采样点,然后叠合土地利用数据,对预采样点进行适当调节以规避诸如水体、城市不透水面等区域.此外,为适应研究区西北部复杂地形条件下土壤水分空间分布的潜在异质性,采样时对样点进行了适当加密,最终确定在空间上尽可能均匀分布的样点共184个.结合手持GPS,土壤采样工作在2016年4～5月间进行,采样时选取有代表性的新鲜土壤,采样深度为0～20 cm.取样时,先铲出一个耕层断面,再平行于断面下铲取土,混合土样太多时用四分法将多余土壤弃去[8],采样的同时记录样点的土地利用类型.采样当天即采用烘干法进行土样含水量测定,具体过程：将鲜土捏碎放入已知准确质量的铝盒中,对其进行称量,精确至0.001 g；然后,将其置于已预热至(105±2) ℃的烘箱中烘烤12 h,取出冷却至室温后立即称重,每个样点做3个平行测定.若第一次测定结果中3个平行样本之间的差距过大(超过0.05%),则进行二次烘干再次称重[9].

利用ArcGIS对土壤含水量样点数据进行数据子集提取,为满足不同插值方法对数据的要求,本文按照7∶3进行插值模型输入数据和插值精度验证数据子集提取.插值模型输入数据集(记为Ⅰ)共129个样点,精度验证数据集(记为Ⅱ)共55个样点.表1为土壤含水量数据的主要统计量信息.从表1可看出,山地平原过渡地带的土壤水分质量分数的变异系数为19%,属于中等变异.

表 1 土壤含水量数据统计结果

3 研究方法

3.1 空间插值方法选择空间插值方法的选择是影响插值精度的主要因素之一[10].不同空间插值方法在数据要求、参数设置等方面均有较大差异[11].本文选择了几种常用的空间插值方法,包括：反距离权重法(IDW)、径向基函数函数空间插值方法、普通克里格插值方法以及结合线性回归模型和普通克里格的回归克里格插值方法.

3.1.1反距离权重(IDW)和克里格插值法(Kriging) 反距离权重法是一种常用又简单的空间确定性插值方法,待插值点的预测值为搜索半径中已知点的值以与待插值点的距离为权重的加权平均,距离待插值点越近的已知样点的权重越大[12].与此类似,克里格空间插值(Kriging)方法亦通过一定范围内的实测样点数据的线性加权预测待插值点,所不同的是,各点的权重由待插值变量的半方差函数得到.IDW和Kriging空间插值方法的一般模型可表示为

(1)

Wi=d-p

(2)

其中,dSi-S0为搜索半径中的点Si到待插值点S0的距离,p为参数.因此,当通过搜索半径或空间插值样点数目确定参数N时,参数p成为IDW插值方法的主要影响因素.

Kriging方法中权重Wi的一般形式为

W

Z(xi+h)]2,

(3)

其中,N(h)为空间距离为h时点对数目,Z(xi)和Z(xi+h)分别为位置xi和(xi+h)处的属性采样值.实际上,Wi(h)可以看作是属性空间自相关的度量函数,称为半变异函数.

3.1.2径向基函数插值方法(RBF) RBF径向基函数插值方法是一种确定性空间插值方法,其基本思路是用一个通过所有属性样点的曲面来预测待插值点,这个拟合曲面拥有最小的空间曲率.径向基函数包括5种常用的基本函数,包括规则样条函数、张力样条函数、高次曲面函数、反高次曲面样条函数和平面样条函数.本文以交叉验证精度为标准选取最优的径向基函数插值结果作为对比.

3.1.3回归克里格插值方法 地物属性受多自然因素甚至社会人文因素的影响,且这些影响因影响因素的空间异质、时间变率不定而对不同区域具有不同的影响深度.一些影响因素空间分布广、影响时间尺度大、空间影响较为均一,使得受其影响的地物属性具有明显的空间趋势.而一些区域性或短时间尺度的因子则多体现为局部影响,这就在由大尺度因素构建的“趋势”上叠加了短程变异.回归克里格插值方法将这种多尺度因素影响下的地物属性进行分离,其概念模型为

(4)

研究表明,太阳辐射、地形、土地利用类型、坡度、坡向、海拔高程等是影响土壤水分的重要因子[13].结合数据的可获取性和研究区特点,本文将土地利用类型、太阳辐射和坡度因子作为土壤含水量“趋势”拟合的考虑因子.需要注意的是,土地利用类型数据是一种“类别集合”形式的空间“软数据”,利用其进行土壤含水量多因素线性回归时需要对其进行“硬化”[14].本文引入虚拟变量对土地利用类型软数据进行“硬化“.在研究区中,结合采样时记录的土地利用现状,最终确定6种土地利用类型：耕地、河滩、未利用地、草地、林地、湿地.根据显著性差异,分为4类:林地归为第一类(C1),耕地归为第二类(C2),荒地归为第三类(C3),河滩地和草地归为第四类(C4).对各个土地利用类型进行土壤含水量方差分析,以此为参考将土地利用类型做聚合归并,并进行虚拟变量处理(表2)和回归因子显著性检验(表3).

因此多元回归预测结果为

μ(fi)=17.780+[landusetype].weight+

0.006*[solar]+0.05*[slope].

(5)

多元回归预测结果和样点实测结果的差值即为局部、短时间尺度因素引起的变异部分(ε(gj)),对其进行空间插值即得到待差值点的变异值,本文采用普通克里格插值方法作为土壤含水量变异值的空间分布.最后根据(5)式得到最终的空间插值结果.

3.2 数据探索和插值参数确定为合理设定空间插值模型参数、满足空间插值模型对输入数据的要求(比如克里格插值要求样点数据满足正态分布),对样点数据的统计特征、空间分布特征及异常值进行探索是空间插值的必要过程.正态QQ图是用于反映样点数据与标准正态分布的接近程度,样点数据的分布越接近一条直线,则越接近于服从正态分布.同时,图中远离标准直线的点为潜在的异常值点,因而QQ正态图还可用于异常值剔除.图2为本文土壤含水量样点数据经异常值剔除后的正态QQ图,可以看出,样点数据具有很好的正太分布特点.如前所述,IDW和Kriging插值方法对待插值点的预测是基于空间距离的,对IDW插值模型中距离参数(或点数量)N和Kriging方法中的权重参数Wi(h)的确定均要求对样点的空间分布特征进行探索.图3为本文土壤含水量样点数据的半变异函数拟合结果.从图3可以看出,研究区土壤含水量最大空间自相关距离(变程)约为1 600 m,因此以1 600 m为IDW空间插值方法中的搜索半径阈值.

图 2 土壤含水量QQ正态分布

图 3 土壤含水量半变异函数拟合

3.3 精度评价指标空间插值精度评价指标可分为相对验证指标和绝对验证指标.本文选择绝对验证指标有均方根误差和平均绝对百分比误差,均方根误差(RMSE)表示预测值与实测值偏差的平方和观测次数比值的平方根,均方根误差对预测中的最大值和最小值比较敏感；平均绝对百分比误差(MAPE)表示预测值和实测值的偏差占实测值的百分比的平均根.相对验证指标选取模拟优度指数(G)它是衡量模型模拟效果是否优于仅对实测值取平均值得的指标.指标选取的优劣对于空间插值精度评价有重要的影响.

4 结果及精度对比分析

为了明确不同插值方法对山地平原过渡带土壤水分空间数据的影响,通过几种插值方法得到结果如图4,并将不同的插值精度进行比较,筛选出适宜过渡带的土壤水分空间插值方法.

由表4可知,对于RMSE,预测误差表现为：

表 4 4种插值方法的精度评估

RBF>RK>IDW>Kriging,说明Kriging的预测误差最小,RBF的预测误差最大；对于MAPE,预测误差表现为：RBF>RK>IDW>Kriging,说明Kriging的精度最高,RBF的预测精度最低；对于G,预测误差表现为：RK>Kriging>IDW>RBF,说明回归克里格插值的精度(0.36)高于普通克里格法(0.27)和反距离权重插值法(0.19)与样条插样法(0.08),其预测模型效果最优.不同检验指标精度评价不同,绝对验证指标显示,克里格插值精度最高,其次为反距离权重法,回归克里格插值次之,径向基函数效果最差.而相对验证指标显示,回归克里格插值高于普通克里格与其他2种方法,说明验证指标的选取也会导致评价精度的不同.

5 讨论

山地平原过渡带具有明显地域分异,特殊的地形地貌造成土壤水分含量区域范围内的差异.不同的插值方法插值精度不同,在众多插值方法中,克里格插值仍然是一种比较适用的方法.在目标变量空间自相关性较弱,与环境因子相关性较强的时候,回归克里格是较为优越的插值方法.本文的回归克里格考虑了山地平原过渡带的土壤水分受多种因素的影响,尝试将坡度、年太阳辐射与土地利用类型等环境因子引入模型.反距离权重法是一种空间确定性插值方法,插值表面中的最小值和最大值出现在采样点处,输出表面对异常值十分敏感.目前对此模型使用的合理性证明仍较困难,但反距离权重法假定数据中存在空间自相关,当目标变量出现局部变异性时,此方法较为适宜.径向基函数法可为平缓变化的表面生成很好的结果,但在表面值在短距离内出现剧烈变化或怀疑样本值很可能有测量误差或不确定性时,这种方法不适用.从其插值效果来看,山地平原过渡带的土壤水分空间分布特征为西北山区的土壤水分含量明显高于东南平原,城镇周边土壤水分含量明显低于周围耕地.

土壤水分对农作物的生长发育有着重要的影响,研究区的土壤水分质量分数为13.22%～72.52%,对于中生性作物而言,对土壤含水量的要求各不相同.马铃薯等豆类作物等最适土壤含水量相当于田间持水量的70%～80%,禾谷类作物为60%～70%.土壤含水量低于最适值时,光合作用降低.各种作物光合作用开始降低时的土壤含水量(占田间持水量之百分数)分别为：水稻57%,大豆45%,小麦41%,花生32%.因此,对于绵竹地区而言,马铃薯等豆类植物适宜种植在西北山区,这里土壤水分充足,利于作物的生长；对于东南平原而言,适宜种植水稻、大豆、小麦、花生等粮食作物和经济作物.土壤水分含量的探讨对于农事活动的安排具有重要的作用.

土壤属性的空间连续数据在理论与实践方面有着重要的作用,如何提高其精度是未来发展的一个重要研究方向.影响空间插值精度的原因除了插值方法的选择外,采样点的密度和数目[15]、插值方法中参数的设定以及空间自相关的范围和程度都是影响其精度的因素[16].对于克里格插值方法而言,获得稳定的半方差函数需要50～100个采样点,采样点位100～150个才可达到一个平稳的变异函数[17].插值参数的设定对其精度有一定的影响,但是目前绝大多数参数的确定还未有太多的规律可循.