基于空间插值的区域尺度土壤全盐量时空变异规律研究

徐存东，谷丰佑，朱兴林，王 燕,3，张 鹏，高德坤

(1.华北水利水电大学水利学院，郑州 450046；2.水资源高效利用与保障工程河南省协同创新中心，郑州 450046；3.河南省水工结构安全工程技术研究中心，郑州 450046)

0 引 言

我国西北地区土地广阔、光热充足，制约大片宜耕土地利用的主要因素是可用水资源短缺，为了改善耕地环境，提高土地利用效率，发展区域农业经济，必须大力发展引水灌溉工程[1]。同时，长期以来人工提灌工程的大面积实施，使得区域水文循环规律随着灌溉方式和引水规模变化而不断改变，也主导着区域尺度的土壤水盐分异过程，改变并控制着区域的水盐时空演化和分布特征[2]。由于该地区高蒸发低降雨的气候条件，促使土壤深层的水分携盐分向土壤表层迁移，随着水分的不断蒸发，盐分滞留在土壤表层，最终导致大面积的土地盐渍化；为了灌区的可持续发展，有必要对区域内土壤盐分的空间分布和演化态势进行科学的评估和预测。当前，土壤全盐量野外监测点数量较为稀少，无法形成对整个区域土壤盐分的大规模监测，从中获取的有限数据不能满足对其时空变化的研究，且监测点在空间上分布不均匀，通过空间插值模型可以有效地对土壤盐分空间分布特征进行预测[3]。如姚荣江、杨劲松等学者利用普通克里金法对苏北海涂围垦区土壤盐分空间分布进行了模拟、插值和预测，为滨海地区盐渍土壤的精确评估与高效改良提供了参考依据[4]；刘鑫等运用反距离权重法对河套灌区的土壤进行空间插值，实现了盐化分布的空间可视化[5]；吴亚坤等运用反距离权重法对南疆绿洲区土壤盐分进行了分布解析得到了较好的插值效果[6]；张飞以新疆精河绿洲为研究区，运用反距离权重法对表层土壤盐分进行空间插值，得到了盐渍化程度与分布情况[7]；姜涵等运用GIS 空间分析法对研究区土壤盐分特征进行研究发现可以通过插值得到未检测地区的数据，能够有效的提高对土壤全盐量的研究[8]。上述研究仅是使用了单一插值方法对研究区土壤全盐量进行分析，没有涉及不同插值方法下所得模拟结果的精度比较，不能验证插值结果的准确性；同时由于研究区的土壤盐分监测点分布较为稀疏，无法全面监测研究区，构建出的土壤盐分变化空间插值模型，在进行土壤盐分时空变化特征研究时有一定的局限性。本文以甘肃省景电一期灌区(以下称“景电灌区”)为研究区，利用反距离权重法、趋势面法、样条函数法、普通克里金法四种不同的空间插值方法对研究区2015年土壤全盐量数据进行插值，将插值结果与监测点的实际观测数据进行对比分析，优选出区域土壤全盐量最佳空间插值方法。并分别对1994年、2001年、2008年和2015年监测数据进行空间插值，对研究区土壤全盐量的空间分布及变化规律进行分析，以期为灌区的土壤盐渍化研究提供数据和理论支持。

1 材料与方法

1.1 研究区概况

景电灌区是地处西北干旱荒漠区的大型扬水灌区，东经103°20′～104°04′，北纬37°26′～38°41′，位于甘、宁、蒙的三省的交界地带。北倚腾格里沙漠，东临黄河, 南靠昌岭山, 西连石羊河流域。地层区划属河西走廊六盘山分区, 武威中宁小区。总灌溉面积6.13 万hm2，分布高程在1 596～1 906 m之间。属于典型的内陆寒旱区气候，少雨干旱、蒸发量大、植被稀疏、风积风蚀严重是其主要气候特征。研究区地表水资源十分匮乏，多年平均年降雨量仅为186 mm，且集中在5-9 月份，蒸发量却高达2 516 mm。研究区地表土层结构松散，为土壤盐分溶质运移提供了良好的条件。加之低降雨高蒸发的气候特征，盐分因子极易留存于土层表面形成积盐，日积月累造成土壤盐渍化[9]。

图1 研究区土壤全盐量监测点位置图

1.2 数据来源

根据《景泰川电力提灌工程第一期工程技施设计报告》《景泰川灌区历年土地调查报告》和《景泰县荒地资源及其开发利用报告》，结合相关资料可以得到1994、2001、2008和2015年的0～100 cm土层的土壤全盐量和50个非插值点土壤全盐量数据，表1为各监测点点位坐标及各个点位地下100 cm内土壤全盐量数据。

1.3 研究原理与方法

1.3.1 插值方法

现有的监测点数量相较于整个研究区面积略显稀少，且分布不匀，有限的数据无法准确反映研究区域土壤全盐量空间分布状况及时空变化规律。因此根据现有的研究区监测点数据，利用ArcGIS中的空间插值模块得出整个研究区域的土壤全盐量插值图，进而分析区域尺度土壤盐分发展特征[10]。ArcGIS空间分析工具中有四种常用的空间插值方法：反距离权重法(IDW)、趋势面法(TL)、样条插值法(Splines)和克里金(Kriging)插值法[11]。反距离权重法(IDW)通过对已知点与预测点之间的距离进行赋权得到插值结果，赋权大小与距离远近呈负相关[12]；趋势面法(TL)以多项式回归分析为基础，遵循最小二乘法原理，对输入的监测点数据进行拟合，反映出土壤盐分的空间发展趋势[13]；样条插值法分为规则样条插值(RS)和张力样条插值(TS)，是由相邻的两个数据点通过数学函数确定一个多项式，多项式两两相连组成一条经过所有数据点的光滑样条曲线[14]；克里金法是在区域化变量和半变异函数基础上，对区域内的变量取值进行线性无偏最小二乘估计，以此反映区域化变量在空间的真实变化特性，其中又以普通克里金法(OK)的应用最为广泛[15]。

表1 研究区监测点点位及其土壤全盐量

1.3.2 精度评价方法

空间插值技术的发展为土壤全盐量空间分析提供了新的方法，然而如何从多种插值方法中选取该研究区域的最优插值模型是空间插值在土壤全盐量分析中的一个关键问题[16]。本文运用交叉验证法对几种插值方法进行精度评价。首先从16个监测点中选出一个点，设其数据为未知值，使用其他15个监测点实测数据构建插值模型，读取该点的插值结果，计算插值结果与所选监测点实际数据之间的误差值。依次从16个监测点中选取下一个点作为未知站点，重复进行操作，获取这16个监测点误差值。即将现有监测点数据分为已知值和验证值，通过不同插值方法对已知值进行插值，读取验证值点位对应的插值结果，与验证值进行对比，插值结果与实际数据之间的误差越小则插值方法越精确[17]。

采用平均绝对误差MAE，平均相对误差MRE，均方根误差RMSE来评价插值结果的精度[18]。

(1)

(2)

(3)

式中：Pei、Pai分别为第i个监测点的预测值和实测值；n为监测点数量。

2 结果与分析

2.1 不同插值方法比较与优选

为了更为准确地显示研究区100 cm内土壤全盐量在水平空间的分布情况，得出土壤全盐量最优插值方法。运用ArcGIS中的Spatial Analyst模块对2015年土壤全盐量数据进行四种空间插值分析，绘制出各插值方法100 cm内土壤盐分的空间分布图。

样条函数插值结果见图2、图3，从图中可以看出运用样条函数法进行插值时，规则样条函数插值中权重为1的时候正负向外延较大，而权重为0时正负向外延较小，出现此现象主要是由于在插值时缺乏数据基础，监测点数据分布不均匀且各点之间变化不够平缓；而张力样条曲线中权重越大越接近真实区间。但是两种函数类型的插值结果负向数值始终为负数，与土壤全盐量必须大于等于零的实际不符。

趋势面法插值结果见图4，趋势面法通过数学函数(各阶多项式)进行插值，可以看出在运用趋势面插值法时，采用二阶多项式插值时与真实数据拟合最好，一阶多项式插值结果区间甚至小于实测数据真实区间。但是由于趋势面法过于注重整体发展趋势，导致插值结果负向数值为负数，不符合实际。

普通克里金插值结果见图5，可以看出普通克里金插值结果综合考虑了距离、方位的影响，对极值点的处理较好，所有插值结果都没有出现负值，且区间与实验数据真实区间基本吻合，仅有OK(指数函数)的区间略小于真实区间。

图2 张力样条函数插值图

图3 规则样条函数插值图

反距离权重法见图6，综合所有插值结果可得运用反距离权重法进行插值时，结果直观，与真实区间符合。

运用误差矩阵(Error Matrix)对3种反距离权重法和其他4种普通克里金法插值结果进行精度评价，以土壤全盐量0.5%为一阶段进行划分，将误差矩阵划分为7个阶段区。对50个非插值点的数据进行检验。用基于误差矩阵的精度评价方法对土壤盐分插值结果与非插值点数据进行评价，kappa系数高于0.60为高度一致性，低于0.60为中度一致性。由表2可知，插值总精度反距离权重法(幂3)最高(为76.00%)，其次是普通克里金法三角函数(为66.00%)，反距离权重法(幂1)最低(为56.00%)。以Kappa系数进行判定，仅有反距离权重法(幂3)属于高度一致性，其余插值结果均属于中等一致性，故反距离权重法(幂3)对土壤盐分插值精度最高。

为选取土壤全盐量最优插值方法，本研究采用交叉验证法通过平均相对误差(MAE)、平均绝对误差(MRE)、均方根误差(RMSE)3个指标来评价不同插值方法插值结果精度。从表3中可以看出，从MRE来看，TS>RS>TL>OK>IDW，说明IDW的插值相对误差最小，张力样条函数法最大；从MAE来看，RS>TS>TL>OK>IDW，即规则样条函数法插值误差最大，IDW的绝对误差最小；从RMSE来看，TS>RS>TL>IDW>OK，说明OK的灵敏度最高，张力样条函数法最低，但IDW灵敏度仅低于OK，且在数值上相差不大。综合3个误差评价指标以及误差矩阵分析结果，可得IDW(幂3)相比较于其他方法在插值精度和预测准确度上均为最佳，所以选择IDW(幂3)对土壤全盐量进行空间插值。

图4 趋势面法插值图

表2 土壤盐分插值精度评价表

2.2 区域土壤盐分插值结果验证

运用ArcGIS软件使用反距离权重法(幂3)对2015年土壤盐分进行空间插值，提取插值图中相应监测点的预测值，选取最近的年份2015年，将所得预测值与实测值进行对比，结果如图7所示。可以看出，运用反距离权重法(幂3)对研究区土壤全盐量进行空间插值得到的预测值与所选年份的各个测点的实际值基本一致。说明利用反距离权重法(幂3)得到的结果图可以准确的反映研究区各个年份的土壤全盐量空间分布特征，并可以此分析土壤全盐量的年际动态变化规律。

图5 普通克里金插值图

图6 反距离权重插值图

表3 5种插值方法插值土壤盐分结果的误差统计

注：MRE为平均相对误差；MAE为平均绝对误差；RMSE为均方根误差。

图7 2015年各监测点预测值与实际值对比图

2.3 区域土壤全盐量时空进程分析

运用反距离权重法(幂3)对1994、2001、2008和2015年土壤全盐量进行空间插值，插值结果见图8。可以看出，土壤全盐量呈现西南低东北高分布状况，从1994年到2015年南部低盐分含量地区不断向北发展至边界，西部土壤全盐量在0.5%～1.0%的地区逐年减少为一较小区域，东北部高盐分含量地区面积虽无明显增长，但其土壤全盐量一直在增大，从1994年的2.0%～2.5%增大到2015年的2.5%～3.0%，并且局部区域的土壤全盐量已增大到了3.0%～3.5%。

图8 各年份土壤全盐量空间插值结果图(单位：%)

通过ArcGIS对图8中各年份的土壤全盐量特征值进行提取。由表4可知，研究区土壤全盐量最小值从1994年到2015年前期增长后期下降，截止2015年相比于1994年仅提高了0.000 8%，几乎没有发生变化。最大值由1994至2015年持续增长，到2015年较1994年提高了0.964 6%。土壤全盐量平均值总体上发生了增长，2015年较1994年总体增长了0.116 7%。

表4 研究区各年份土壤盐分特征值 %

3 结 论

(1)文中所述的5种插值方法中，规则样条函数、张力样条函数和趋势面法插值结果负向数值为负数，与实际情况不符；普通克里金法除了以指数函数为插值函数的插值结果外，都与实验数据基本吻合，但是反距离权重法结果更加直观，且与真实区间符合，经交叉验证可知，反距离权重法插值结果与实际的监测结果吻合度较高，结果最为精确，可较准确预测灌区土壤全盐量的变化趋势，为灌区区域尺度盐分空间处理分析提供指导。

(2)使用反距离权重法对研究区1994、2001、2008和2015年土壤盐分含量进行空间插值，插值结果表明，研究区土壤盐分含量平均值由1994年的1.006 1%增大为2015年的1.122 8%，除少数区域由于长期以来的灌溉洗盐，土壤盐分有所降低外，其他区域土壤盐分总体上呈现增大趋势，且土壤盐渍化程度也越来越严重。

(3)根据研究区土壤盐渍化发展趋势，针对轻度盐碱化区域，引进先进的灌溉技术如喷灌、滴灌等，加强灌溉水的利用率，控制地下水位的同时加速土壤表层脱盐；针对中度盐碱化区域，在淋溶洗盐的基础上，大量种植盐生植物扩大地表植被覆盖，降低土体与大气间的水热交换速率，以此来减小土壤表面积盐的动力；针对重度盐碱化区域，前期通过土地平整、农田水利工程等措施有效降低土壤内含盐量，后期通过深耕晒垡结合有机肥、土壤翻耕结合生物菌肥等综合措施改善土壤生态结构，降低土壤盐碱化程度。

□