基于阻隔因素的耕地质量分等因素插值方法研究

杨永侠 郭雅萍 张 函 张丽红 桑 婧

(1.中国农业大学土地科学与技术学院， 北京 100083; 2.自然资源部农用地质量与监控重点实验室， 北京 100035)

0 引言

耕地质量调查评价是新时期下对耕地数量、质量和生态“三位一体”保护的重要基础，同时也是有效保护耕地资源、保障国家粮食安全、维持农业可持续发展的重要基石，可为解决国家耕地的占补平衡、基本农田划定、土地利用规划等问题提供科学依据[1-3]。

耕地质量评价是以县为单位划分耕地质量分等单元，县域内动辄有过万分等单元，其数量大、分布广，所以分等因素一般通过样点(样地)数据采用空间插值的方法获得。对于平原地区或者分等因素变化不大的区域，一般插值方法可以较好地拟合实际情况，但当区域内地形复杂、海拔落差较大时，利用克里金或者反距离权重插值的方法不能很好地拟合实际情况。目前，对于此类情况主要是通过增加调查样点提高插值的准确性，但是位于山地丘陵区耕地的地形变化过于复杂，增加控制样点的方法同样也会消耗过多的人力物力，且插值精度的提高并不明显。

在复杂地区土壤属性空间预测方面，国内学者对空间预测模型精度的提高[4-6]及预测性制图[7-12]等进行了研究，并已取得多方面的研究进展。经研究证实，利用辅助变量可以提高土壤属性的空间预测精度，但是不能很好再现原始区域变量的空间结构及变量之间的相互关系，而且不能进行定量的不确定性分析，这使与辅助变量结合的预测结果应用受到一定的限制。

基于以上问题，本文在耕地质量评价体系的基础上，提出山区分等因素基于阻隔因素的插值估算方法，以期提高耕地质量等级评价结果的准确性，为第三次全国土地调查耕地等别调查评价工作提供借鉴。

1 研究区和数据

1.1 研究区概况

研究区选择青海省海东市平安区湟水流域中游南侧，海东市中心腹地，西距西宁市35 km处，地理坐标为东经101°49′～102°10′、北纬36°15′～36°34′之间，境内南北长约33.3 km，东西宽约23 km，总面积769 km2。

平安区属大陆性气候，春季干旱多风，夏季凉爽，冬季寒冷。年平均气温7.6℃，年内降雨分布不均，无霜期218 d。平安区境内多为山区，且有东沟、六道沟、老虎沟等12条沟岔，成为山峦起伏、山川相间的主要河谷地带。平安区境内的大部分地区海拔在2 038～4 127 m之间(图1)，起伏连绵的地势以及复杂多样的气候，构成了农业生产在空间上比较明显的地区差异，最终形成了川水、浅山、脑山3种地貌类型区。

图1 平安区海拔分布Fig.1 Elevation map of Ping’an District

1.2 数据来源与处理

采用数据来源于青海省平安区2007—2011年5年的测土配方施肥项目的调查数据，在综合分析本研究区内耕地分布特点的基础上，均匀选取了280个调查样点数据，利用ArcGIS中子集要素工具(Subset Features)随机确定240个点作为插值点，40个点作为验证点进行精度评估(图2)。

图2 插值样点分布Fig.2 Interpolation sample distribution map

2 研究方法

2.1 可插值分等因素分析

《农用地质量分等规程》[13]中指出，农用地质量分等因素分为推荐因素和自选因素两类。推荐因素由国家统一确定，并且分区、分地貌类型给出，自选因素由省级土地行政主管部门确定[14]。总的来说推荐因素和自选因素可以分为4种类型，即土壤条件、地貌条件、水文类以及农田基本建设类，每种类型所涉及的分等因素如表1所示。

平安区在农用地质量分等过程中选用的分等因素有：表层土壤质地、地形坡度、灌溉保证率、土壤pH值、土壤有机质含量、有效土层厚度。由于地形坡度这个分等因素的空间相关性较差且不满足地理学第一定律；而表层土壤质地和灌溉保证率的属性值在空间上也不连续，故均不适合利用插值对未知样点值进行估算；因此本文考虑到分等因素的含义、属性值类型、插值方法，选取了平安区有机质含量、土壤pH值和有效土层厚度3个分等因素进行插值方法研究。虽仅根据3个因子不能对耕地质量进行分等定级，但这3个因子在山地丘陵区耕地质量评价时所占的权重均比较大，且其属性值会很大程度影响最终的评价结果，故需保证这三者的准确性。

表1 分等因素统计Tab.1 Statistics of factors

2.2 传统的空间插值方法

2.2.1反距离权重插值

反距离权重法是空间插值的常用而又简单的方法。它以样本点间距离为权重进行加权平均，距离插值点越近的样本所占权重越大，这种方法简单方便，效果直观而且效率较高[15]。但如果样本中有极端值则会影响较大范围的插值结果。其计算公式为

(1)

反距离权重法和克里金法的计算公式相同，只是在反距离权重法中权重λi仅取决于预测位置的距离。

2.2.2样条函数插值

样条函数法是指通过使用最小化整体表面曲率的函数来估计值，生成恰好经过样本点的平滑表面，如同将一个软膜插入并经过各个已知样点，并且表面的总曲率最小[17]。此方法最适合生成平缓变化的表面，比如海拔、地下水位等。其计算公式为

(2)

式中S——所求点函数值

T——原始点函数值

R——距离函数N——样点数

λj——通过求解线性方程组而获得的系数

rj——点(x,y)与点j之间的距离

2.3 基于阻隔因素的插值方法

在山地丘陵区利用插值对未知样点进行估值时，区域内极有可能会有较大的山体或者河流分布，对插值效果影响较大，所以在进行插值分析时应考虑这一因素。已有的研究表明，土壤属性具有高程地带性的分布规律，并且可以用高程作为辅助变量对土壤的属性进行模拟预测，但是此类研究侧重于对原数据的优化，不能反映数据本身的结构与特征；因此本文在不改变原数据的状态下，从模拟预测过程中引入阻隔因素，便于更好地贴合土壤属性的实际分布状况。

阻隔因素是指在研究区中会中断表面连续性的线状要素的位置。山地丘陵区实际工作中阻隔因素可能为山脊、山谷、河流，这样的地形因素存在时间长，对周围土地影响大，会使阻隔因素两侧的数据不连续，拟合的数据曲面容易发生突变。尤其山体阻隔影响，迎风坡与背风坡风化差异，阳坡与阴坡光照辐射差异，气流和能量对土壤变化影响尤为显著。而且在山地丘陵区，耕地多分布于山谷和河谷两侧，所以本研究选取海拔落差较大的山体与大型河流作为阻隔因素进行研究。

基于阻隔因素的插值方法具体步骤为：首先提取阻隔因素，本文中阻隔因素提取涉及河流和山体。大型河流分布数据可以通过实地调查或者河流分布图获得，方法较为简单，本文不做赘述；山体阻隔因素提取，山脊线的获取采用水文分析的方法[16]，算法思路为：山脊线同时也是分水线，对于在分水线上的栅格，是水流的起源，通过地表径流模拟计算，应当只有流出方向的表面径流，没有流入方向，也就是此处栅格的汇流量为零。提取汇流值为零的栅格，就可以得到山脊线。然后提取不穿越耕地，而且其两侧距离耕地较远的山脊线为山体阻隔线。通过矢量化、连接、平滑等操作确定分水线的位置。利用ArcGIS软件对DEM计算地表径流，提取汇流量为零的栅格，经过筛选、矢量化、平滑得到平安区阻隔因素分布图，如图3所示，将提取的阻隔因素加入插值的过程中。ArcGIS软件中的反距离权重插值工具提供“输入折线障碍要素”选项，此处的折线要素是指在搜索输入采样点时用作中断或限制的折线要素，因此可以将提取的阻隔因素作为其折线要素，其他参数采用默认设置进行插值分析；ArcGIS软件也提供了“含障碍的样条函数”插值工具，同样将阻隔因素作为障碍要素进行输入，其他参数保持不变。

图3 阻隔线分布图Fig.3 Block line distribution map

2.4 精度检验方法

采用交叉验证[18-20]对反距离权重插值、样条函数插值、基于阻隔因素的反距离权重插值、基于阻隔因素样条函数插值这4种插值方法的精度进行对比分析，采用平均绝对误差(MAE)、平均相对误差(MRE)、均方根误差(RMSE)、相对均方根误差(RMSEr)4个评判标准来评价插值结果。其中，平均绝对误差能够估算出获取的估计值的误差范围；平均相对误差可以判断出估计值和实测值之间的误差；均方根误差反映利用样点数据的估值灵敏度和极值效应，其值越小越好；相对均方根误差越小则误差越小，精度越高[21-23]，计算公式如下

(3)

(4)

(5)

(6)

式中xi1——验证样点的测量真值

xi2——验证样点预测值

n——验证样本总数

3 结果与分析

3.1 原始数据描述性统计分析

对原始样本的有效土层厚度、有机质含量、土壤pH值进行基本的描述性统计，从表2可知，有效土层厚度的范围为27～153 cm，且极差最大；有机质质量比范围是9～62 g/kg，土壤pH值极差最小；偏度和峰度是描述数据是否符合正态分布的2个重要指标，当二者都为0时，表明数据符合标准的正态分布。当偏度大于0时，表明数据集中分布在左侧，向右延伸，小于0时，集中分布在右侧，向左延伸。当峰度大于0时，整个数据分布形态比标准的正态分布高耸，数据多集中分布在平均值附近，小于0时，整个数据分布形态比标准正态分布平坦，数据分布较为分散。由表2可知，3个分等因素均不符合标准的正态分布；变异系数是反映数据分布状态的一项重要指标，主要反映数据的离散程度。由表2可知，有效土层厚度和有机质含量属于中等变异性，pH值属于弱变异性，所以有效土层厚度和有机质含量较pH值更易受环境变量的影响而变化。

表2 原始样本描述性统计特征值Tab.2 Descriptive statistical eigenvalues of original sample

3.2 插值结果分析

3.2.1土壤有机质含量结果分析

调查样点有机质含量通过反距离权重、基于阻隔因素的反距离权重、样条函数、基于阻隔因素的样条函数4种插值方法获得的插值结果如图4所示。

图4 有机质含量插值结果对比Fig.4 Comparisons of organic matter content interpolation results

土壤有机质不仅是植物养分供给的源泉之一，而且是保持土壤良好的物理性质的物质。平安区全区耕地土壤的有机质质量比的范围在9～62 g/kg之间，就全区而言，平安区的有机质水平均较低。由图4可以看出，这4种不同的插值方式均能较好地模拟平安区土壤有机质的分布状况，呈现出由北向南逐渐递增的趋势。

图5 土壤pH值插值结果对比Fig.5 Comparison of soil pH value interpolation results

4种插值方法得到的有机质含量的空间分布相一致，但其插值结果却存在一定的差异，通过反距离插值、样条函数插值、基于阻隔因素的反距离权重、基于阻隔因素的样条函数插值得到的有机质含量范围分别是9.065 94～61.966 1 g/kg、8.989 31～73.701 6 g/kg、9.064 9～61.966 1 g/kg、8.896 7～72.529 7 g/kg，通过对比可知，反距离权重插值的结果与有机质真实值的范围最为接近；从平滑度和局部变异性来看，反距离权重插值出现了“牛眼”现象，样条函数插值结果图较平滑、连续；且对比基于阻隔因素与未基于阻隔因素的反距离插值和样条函数插值结果，可以明显看出，基于阻隔因素后对其两侧数据的插值产生的影响，能更真实反映土壤有机质分布情况的细节特征，因此基于阻隔因素的插值方法适用性更强。

对40个实际调查样点数据的插值精度采用交叉验证的方法进行分析，从表3可看出，基于阻隔因素的样条函数插值方法精度最高，比样条函数插值的平均绝对误差、平均相对误差、均方根误差、相对均方根误差分别提高了11.23%、10.98%、7.54%和9.20%；比反距离权重插值方法的平均绝对误差、平均相对误差、均方根误差、相对均方根误差分别提高了31.27%、25.60%、31.74%和21.39%。反距离权重、样条函数、基于阻隔因素的反距离权重和基于阻隔因素的样条函数4种插值方式的实测值与预测值的决定系数分别为：0.876、0.920、0.878和0.927；综上可知基于阻隔因素的样条函数插值方式最优。

表3 有机质含量插值结果精度Tab.3 Accuracy analysis of organic matter content interpolation results

3.2.2土壤pH值结果分析

调查样点土壤pH值通过反距离权重、基于阻隔因素的反距离权重、样条函数、基于阻隔因素的样条函数4种插值方法获得的插值结果，如图5所示。

平安区的全区耕地的土壤pH值范围在7.83～8.56之间，呈现出从中部分别向南北递减的分布趋势。由图5可看出，4种插值方法均能较好地反映平安区全区的土壤pH值的空间分布特征。

4种插值方法得到的土壤pH值的空间分布相一致，但其插值结果范围却存在一定的差异，通过反距离插值、样条函数插值、基于阻隔因素的反距离权重、基于阻隔因素的样条函数插值得到的土壤pH值范围分别是：7.831 34～8.539 96、7.256 99～8.645 68、7.830 00～8.539 96、7.758 4～8.553 67，通过对比可知，反距离权重法对pH值估算的结果与真实值的范围最为接近；但是采用反距离权重进行插值的结果在插值点高值附近会出现“牛眼”，这是由于反距离权重插值是一种精确插值方式，不改变插值点原始值；耕地质量调查评价工作要求，在插值过程中应该保证调查样点的值不变，才能保证耕地质量评价的准确性；分别对比基于阻隔因素与未添加阻隔因素的反距离插值和样条函数插值结果，可以得知基于阻隔因素的插值结果反映细节比较清晰，因此基于阻隔因素的反距离插值方法对土壤pH值进行插值的适用性最优。

从表4可知，基于阻隔因素的反距离权重插值方法精度最高，比反距离权重插值的平均绝对误差、平均相对误差、均方根误差、相对均方根误差分别提高了2.94%、2.56%、1.85%和1.90%；比样条函数插值方法的平均绝对误差、平均相对误差、均方根误差、相对均方根误差分别提高了11.19%、11.63%、14.39%和14.88%；反距离权重、样条函数、基于阻隔因素的反距离权重和基于阻隔因素的样条函数4种插值方式的实测值与预测值的决定系数分别为0.527、0.560、0.554和0.476；虽然样条函数和基于阻隔因素的反距离插值方法两者的决定系数均比较高，且相差不大，但是后者插值精度更高，故研究区内对土壤pH值最适合的插值方法是基于阻隔因素的反距离权重插值方法。

表4 土壤pH值插值精度Tab.4 Soil pH value interpolation accuracy

3.2.3有效土层厚度结果分析

调查样点有效土层厚度通过反距离权重、基于阻隔因素的反距离权重、样条函数、基于阻隔因素的样条函数4种插值方法获得的插值结果如图6所示。

图6 有效土层厚度插值结果对比Fig.6 Comparison of effective soil thickness interpolation results

平安区耕地的有效土层厚度在27～153 cm之间，其分布特征是：东北部和西南部最低，西北部的有效土层厚度最高；且呈现出从中部向东西两侧递减的趋势。由图6可以看出，4种插值方法均能较好地反映土壤有效土层厚度的空间分布特征。

4种插值方法得到的有效土层厚度的空间分布相一致，但其插值结果却存在一定的差异，通过反距离权重、样条函数插值、基于阻隔因素的反距离权重、基于阻隔因素的样条函数插值得到的有效土层厚度范围分别是：27.567～152.824 cm、26.002 7～204.752 0 cm、27.567～152.872 cm、26.147 3～166.074 0 cm，通过对比可知，反距离权重结果与真实值的范围最为接近，样条函数插值次之；同样采用反距离权重进行插值的结果与土壤pH值类似，在插值点高值附近会出现“牛眼”，且其空间分布图过渡平滑性欠佳，可观性较差；样条函数插值后的结果较好，图形表面平滑渐进；且分别对比基于阻隔因素与未基于阻隔因素的反距离权重和样条函数插值结果图，基于阻隔因素的插值结果反映了阻隔因素对于其两侧影响的细节，形成明显的分界线，其插值总体保留了丰富的细节，还反映了分等因素的空间渐变特征，更符合其实际的状况，因此基于阻隔因素的样条函数插值方法对土壤有效土层厚度的适用性最优。

从表5可知，基于阻隔因素的样条函数插值方法精度最高，比样条函数插值的平均绝对误差、平均相对误差、均方根误差、相对均方根误差分别提高了15.08%、11.74%、17.41%和9.40%；比反距离权重插值方法的平均绝对误差、平均相对误差、均方根误差、相对均方根误差分别提高了26.33%、23.03%、25.05%和24.63%；反距离权重、样条函数、基于阻隔因素的反距离和基于阻隔因素的样条函数4种插值方式的实测值与预测值的决定系数分别为0.869、 0.852 、0.865和0.901；故基于阻隔因素的样条函数的插值方法最适宜有效土层厚度的区域空间预测。

表5 有效土层厚度插值精度Tab.5 Interpolation accuracy of effective soil thickness values

4 讨论

4.1 基于阻隔因素插值结果与耕地质量年度更新评价结果对比

将土壤有机质含量、土壤pH值和有效土层厚度3个分等因素插值结果与耕地质量评价成果中因素指标值的值域和均值对比可得，分等因素的值域与均值均发生明显变化(表6)。基于阻隔因素的样条函数插值得到的耕地有机质含量为11.84～35.61 g/kg，均值为19.63 g/kg，比年度评价成果均值26.45 g/kg低6.82 g/kg，其平均绝对误差、平均相对误差、均方根误差及相对均方根误差分别提高了 87.35%、87.69%、87.78%和88.52%；基于阻隔因素的样条函数插值得到的耕地有效土层厚度为43.74～149.15 cm，均值为73.59 cm，比年度评价成果均值76.25 cm低2.66 cm，其平均绝对误差、平均相对误差、均方根误差及相对均方根误差分别提高了51.88%、52.80%、63.18%和61.89%；基于阻隔因素的反距离权重插值得到土壤pH值的值域为7.98～8.37，均值为8.26，比年度评价成果均值8.28低0.02，其平均绝对误差、平均相对误差、均方根误差及相对均方根误差分别提高了72.32%、72.36%、69.19%和69.25%。3种分等因素利用基于阻隔因素得到的结果，耕地图斑地块的分等因素值均有减少，但总体来看，值域和均值变化不大，且对土壤有机质、有效土层厚度和土壤pH值3个分等因素的插值精度均有较大提高。在本文研究中发现，基于阻隔因素的确定性插值方法对不同的土壤性质表现不一致，土壤有机质和有效土层厚度最适宜的插值方式是基于阻隔因素的样条函数插值，而在土壤pH值插值中，基于阻隔因素的反距离权重插值精度最高，这是由于反距离权重和样条函数插值原理不同和3个分等因素原始数据的特点各异造成的。由表2可知，土壤pH值属于弱变异性，方差最小，相较于另外两个分等因素其数据分布较为集中；而反距离权重插值由于其根据插值点与样本点之间的距离为权重进行加权平均，在样本分布较为集中、均匀时效果最佳。

表6 插值结果与年度更新数据对比Tab.6 Comparison of interpolation results and annual update data

4.2 插值结果的进一步修正

在实际情况中对于山体、河流而言，迎风坡与背风坡风化存在着差异，阳坡与阴坡光照的辐射也不相同，气流和能量对土壤变化影响也尤为显著。而且在山地丘陵区，耕地多分布于山谷和河谷两侧。土壤有机质、土壤pH值和有效土层厚度的实际分布特征不仅会受到地形地貌的影响，也会受到其他自然因素与人为因素的影响。土壤有机质的具体分布会受到高程[24]、坡度、地形湿度指数、土壤类型[25]等自然因素和种植模式[26]、施肥结构[27]、灌溉、土地利用方式等人为因素的影响，其中高程在大尺度上与土壤有机质的关系最为明显，土壤类型主要是由于成母土质的差异影响土壤有机质的分布[28]。在河流流经的地区土壤质地对土壤有机质的空间分布也有显著的影响[29]；土壤的pH值影响土壤中某些养分元素的有效性，从而影响植物生产，而施用肥料也会对土壤pH值产生影响，相关研究表明土壤pH值与同样也会受到地形地貌、成母土质[30]、土壤类型[31]等自然因素和农耕、施肥措施、工业发展等人为因素的影响；土层厚度关系到植物的扎根条件，而且与养分和水分的保蓄能力相关，薄层土壤对深根性植物的生长有限制，而有效土层厚度和地形地貌、土壤类型、土地利用类型、岩石裸露率、植被覆盖度有一定的相关性。

对平安区有机质含量、有效土层厚度以及土壤pH值3个分等因素的属性值采用基于阻隔因素的方法对其进行了插值对比。但是分等因素种类繁多且分布类型不一，各个地区采用的分等因素也不尽相同，而本文仅对山地丘陵区耕地质量分等因素属性值的预测提出了新的方法，并非能适用于所有地区的所有分等因素；而且本文研究的范围仅是小尺度(县区级)，若获取更为准确的大尺度范围内(地形复杂区)的分等因素属性值，还需要在本方法的基础上考虑多重环境因子。

5 结论

(1)采用不同的插值方法对土壤有机质含量、土壤pH值、有效土层厚度3种分等因素进行插值后，3种分等因素的插值结果的空间分布趋势的模拟与其实际分布基本一致。

(2)由基于阻隔因素的插值方法与不添加阻隔因素的插值方法进行对比可知，基于阻隔因素的插值结果反映了阻隔因素对于其两侧影响的细节，形成明显的分界线，其插值总体保留了丰富的细节，还反映出分等因素的空间渐变特征，更符合其实际的状况，同时也具有更高的预测精度。

(3)对研究区的土壤pH值进行插值，基于阻隔因素的反距离权重插值方法优于其他3种插值方法。

(4)对研究区的土壤有机质含量和有效土层厚度进行插值，基于阻隔因素的样条函数插值方法的精度均明显优于其他3种插值方法，且二者实测值与预测值之间的决定系数分别为0.927和0.901；对土壤有机质含量的插值，基于阻隔因素的样条函数插值方法比样条函数插值的平均绝对误差、平均相对误差、均方根误差、相对均方根误差分别提高了11.23%、10.98%、7.54%和9.20%；对有效土层厚度的插值，基于阻隔因素的样条函数插值方法比样条函数插值的平均绝对误差、平均相对误差、均方根误差、相对均方根误差分别提高了15.08%、11.74%、17.41%和9.40%。