基于地形分类的黄土丘陵区土壤有机质插值方法

张婧 毕如田 丁皓希 许月凤 文伟杰

摘要：黄土丘陵山区地形复杂，耕地细碎，采样点密度不足，因此研究适宜该区域的空间插值方法具有重要意义。以山西省偏关县为研究区，基于数字高程模型（digital elevation model，简称DEM）数据生成地形因子，依据地形因子将研究区划分为不同地形单元，利用Kriging插值法对各地形单元进行分层插值，通过叠加不同地形单元的插值结果，得到全区域的土壤有机质含量空间分布，并利用均方根误差和决定系数评价空间插值方法的预测精度。结果表明，采用划分地形单元分层克里金（Kriging）插值法得到的预测值与实测值的决定系数为0.375 3，明显高于全局插值法;均方根误差精度高于全局插值法，预测值无偏性好。因此，基于地形单元的空间插值方法可以更精确有效地获取土壤有机质空间分布特征，为复杂山区低密度采样下的耕地质量调查与评价提供技术参考。

关键词：丘陵山区;地形单元;土壤有机质;Kriging插值;DEM数据

中图分类号： S159文献标志码： A

文章编号：1002-1302（2021）08-0182-07

收稿日期：2020-11-25

基金项目：国土资源部公益性行业项目（编号：201411007） 。

作者简介：张 婧（1995—），女，山西吕梁人，硕士研究生，研究方向为土地利用与信息技术。E-mail：zj972775@163.com。

通信作者：毕如田，博士，教授，博士生导师，主要从事资源环境信息技术、土地复垦与规划研究。E-mail：brt@sxau.edu.cn。

有机质是土壤的重要成分和肥力指标，精确掌握耕地土壤有机质的空间变异规律，有利于耕地精准施肥和耕地养分的科学管理[1-3]。土壤在时间和空间是一个连续的变异体，在空间上具有相关度很高的异质性，同时具有尺度依赖性[4]，而作为土壤的重要养分之一，有机质也具有同样的时空特征[5-6]。

大量研究表明，土壤有机质的空间变异受多种因素的影响，例如地形因子、成土母质、耕作措施、土地利用方式等，地统计学的研究中，常忽略这些因素对有机质质量比空间再分配的协同作用。付晓婷提出了基于地貌类型分类，研究土壤养分空间变异特征方法，但其划分标准不是很明确，对于一个县域来说尺度比较大，且指标单一，这是传统地形分类现存的问题[7]。对于地貌类型划分，韦金丽等以数字高程模型（digital elevation model，简称DEM）为基础数据，提取坡度，设计分级标准然后分类组合进行大致的地形因子归类，此方法划分结果比较粗糙，但为划分地形单元提供了一定的参考[8]。耕地土壤有机质空间变异规律主要受地形显著影响，在田块尺度上，Yanai等选择哈萨克斯坦北部半干旱区探讨影响耕地有机质变异的因子，结果表明地形是主要因素[9];在空间多尺度上，朱洪芬等选择太原盆地为研究区，分析不同地形部位土壤有机质变异特征与影响因子的空间多尺度关系，并且得到了较好的预测结果[10]。赵业婷等运用描述性统计结合地统计理论，研究了耕地养分空间特征规律以及相互之间的影响关系[11];耿广坡等利用克里金（Kriging）插值结合遥感的方法探讨了各种地形因子与耕层土壤养分含量空间分布规律的关系[12]。对土壤属性空间预测最常用的方法是Kriging插值法[13-14]，但是普通Kriging插值法比较适用于平坦、地形单一的区域，对于地形复杂地区插值效果一般。根据地理学第一定律（土壤屬性信息比较类似是因为具有相同空间位置属性的因子特征越相近造成的，土壤属性差异比较大是因为空间位置不同）可知，相同的地形单元里，在其地形条件相似的情况下，空间距离上呈现出较细微的差异，因此土壤有机质空间变异也服从这一规律;而对于不同的地形单元，由于地形条件不同，造成空间距离很大差异，因此各种地形因子对土壤有机质空间分布规律有很大的影响。田瑞云等基于多种地形因素，对DEM数据进行分析，分类组合划分出典型地形单元[15]。黄魏等选取地形因子研究土壤有机质空间变异规律特征，提出基于地形单元的空间插值方法，获取较高精度的有机质空间分布以及在地统计中地形的协同影响，为精细研究土壤属性局部空间规律特征提供了具有实践性的参考价值[16]。

针对黄土高原丘陵山区，由于气候条件与耕作方式基本相同，因此，地形是引起土壤有机质空间变异的重要原因。本研究尝试依据地形因子特征属性将研究区划分为不同地形单元，单独分析各地形单元下的土壤有机质属性的空间分布特征，以期为耕地土壤质量调查、土地资源管理以及精准农业的实施提供参考依据。

1 材料与方法

1.1 研究区概况

选择山西省忻州市偏关县为研究区，地处晋西北的晋蒙交界处（39°12′～39°40′N，111°22′～112°01′E），北依长城，西邻黄河，与内蒙古准格尔旗隔河相望，属北温带大陆性气候，年均气温为3～8 ℃，年均降水量为425.3 mm。研究区处于黄河中游黄土丘陵区，县域内丘陵起伏，沟壑纵横，地势东高西低，地形复杂多样，海拔高度为877～1 844 m。研究区耕地面积为37 131.03 hm2，耕地细碎，分布于全县各种地形类型，如图1-a所示。偏关县成土母质主要为黄土母质、人工堆垫、沙质黄土，土壤类型主要为风沙土、黄绵土和栗褐土。收集3 007个采样点，约每12.3 hm2 有1 个采样点，不同地形采样间距不均匀，如图1-b所示。

1.2 数据来源及处理

本研究所需数据包括分辨率为12.5 m的DEM数据[来源于中国科学院计算机网络信息中心地理空间数据云网站（www.gscloud.cn）]、土壤有机质数据（来源于山西省偏关县2010年测土配方项目采样）、耕地数据（来源于2015年1 ∶10 000偏关县土地利用变更调查数据库以及2017年偏关县耕地质量等级更新成果）。以上数据均经空间数据格式转换、 校正并投影转化为2000国家大地坐标系（CGCS2000）[17]。利用ArcGIS软件对数字高程模型进行地形因子提取;采用SPSS软件对样点数据进行描述性统计分析。

2 研究方法

2.1 地形单元划分及设计

地形因子代表的空间特征与环境效应相比有很大的差异，所以在土壤属性的空间变异规律中，地形可以作为衡量整个空间特征的主要影响因素[18-19]。本研究根据研究区实际情况选取坡度（slope）、坡向（aspect）、地形位置指数（topographic position index，简称TPI）等3个地形指标，利用 12.5 m 精度高程数据，在ArcGIS软件中分析得到3个地形指标，其中运用反距离加权修正法对地形位置指数进行修正[15]。地形因子在各类地形条件下分布特征具有一定的规律性，依据这一属性的特征，可将偏关县耕地划分为不同的地形单元。本研究主要从DEM中提取出各因子，依据各因子属性特征重分类组合，运用叠加分析划分出地形单元。由于受自然条件的影响，研究区地形复杂，地形类型多元，依据研究区的实际情况，本着真实性、典型性、全面性的原则选取11种地形类型，运用栅格计算器按照图2中地形因子分类组合形成目标单元类型。其中，依据耕地坡度分级标准将坡度分为0°～6°、6°～15°、15°～25° 3级;根据光照方向将坡向划分为阳坡、半阳坡、阴坡和半阴坡4个坡向;根据地形位置指数（TPI）结合坡度、坡向和土壤有机质空间属性，将耕地划分为坡脚、平地、平顶峰、沟谷、阴坡、半阴坡、阳坡、半阳坡、陡坡、山脊、坡肩等共11个地形单元。根据表征出来的实际情况，提出地形位置指数值的分级标准对TPI进行分级，其中SD表示它的标准偏差[20]，划分指标见图2。最后依据地形单元结果分区划分采样点，获取各种地形下的土壤采样点，为分地形单元插值作基础。

2.2 全局空间自相关

全局空间自相关通常用来反映地理实体的空间关系，而耕地土壤有机质在空间分布上变异规律也具有差异性。本研究选择全局莫兰指数（I）来衡量全局空间自相关，公式如下：

I=∑ni=1∑nj=1Wij（Xi-X）（Xj-X）S2∑ni=1∑nj=1Wij。（1）

式中，S2=∑ni=1（Xi-X）2n为样本方差;Wij是要素i和j之间的空间权重;∑ni=1∑nj=1Wij为所有空间权重之和;n为采样点数量;Xi、Xj为空间中i和j要素的属性值;X表示土壤有机质实测值的平均值。

将各类地形单元下的土壤有机质作为观测变量，通过此方法分析偏关县不同地形单元内的土壤有机质空间变异特征。由于各地形单元面积大小不等，各地形单元采样点数量不一样，土壤有机质蕴含内在方差不稳定可能不满足空间自相关分析中隐含的数据稳定性假设，推断出错误的结果，所以进行全局自相关分析时，首先对土壤有机质数据进行标准化[21]。

2.3 聚类分析

利用统计软件SPSS采用欧氏距离和组平均法对11个地形单元进行系统聚类，将地形单元聚类结果在ArcGIS软件中显示。

2.4 基于地形单元的Kriging空间插值方法

首先采用地统计学中普通克里金插值法对整个研究区有机质进行空间插值，得到预测结果;然后单独对不同地形单元下的土壤采样点进行分层插值，利用耕地图斑提取与之相对范围内采样点的预测结果，得到土壤有机质空间分布图，再将各地形单元下的空间预测图拼接，获得整个县域的空间分布预测图。最后设置验证点，各地形单元驗证点个数不能低于采样点个数的5%，运用验证点位上的土壤有机质的实测值和预测值结果验证分地形插值和直接插值拼接获得的预测结果精度。

南非矿山工程师Krige最早提出了克里金插值方法，然后经由法国地理数学家Matheron进行了插值方法的优化[22]。克里金插值是在有限区域内对区域化变量进行无偏最优估计的一种方法[23]，这主要是由变异函数和结构分析为理论支撑。

普通克里金法由于假设条件相对较少和需求参数较简单，所以在克里金方法中是最为普遍的，公式为

Z*（X0）=∑ni=1λiZ（Xi）。（2）

式中：Z*（X0）为估算点土壤有机质含量的预测值;λi为参与插值的位置对预测点的权重;Xi表示土壤采样点的实际位置;Z（Xi）表示采样点的实测值。

2.5 预测结果准确性检验

本研究选取均方根误差（RMSE）、平均预测误差（ME）作为指标，对于研究结果合理性的验证是通过比较不同地形单元验证点位的预测值与实测值之间的差异，利用2个指标分析2种插值方法的差异和好坏，展开精度评价，从而验证本研究提出方法的可行性。其中RMSE值越小，说明所预测精度越高，ME接近0，满足预测的无偏性最优，反之则说明预测精度低，不满足预测的无偏性。

3 结果与分析

3.1 地形单元结果聚类分析

首先把11个地形单元的数据标准化，然后采用欧氏距离和组平均法进行系统聚类，将11个地形单元分成4类。沟谷单独为一类;坡脚、阴坡、半阴坡、平地合并为坡脚类;陡坡、阳坡、半阳坡、坡肩合并为坡肩类;山脊和平顶峰合并为山脊类，如图3所示。根据聚类结果结合实地数据把规律具有一致性的合并为一个单元，最终将研究区耕地划分为4个地形单元（图4）。根据调查与实地数据相符合。

3.2 土壤有机质含量统计分析

依据地形单元划分采样点，土壤有机质含量描述性统计特征见表1。区域内土壤有机质含量呈正态分布特征，由平均值和标准差可知，不同地形单元内土壤有机质含量离散度较相近;由偏度和峰度

可知不同地形单元的有机质含量空间分布特征规律比较明显，符合克里金插值的要求，也为下一步对土壤有机质含量进行半方差函数结构研究以及相应地形单元插值提供了依据。

3.3 土壤有机质空间自相关分析

在ArcGIS中运用空间统计模块分别计算各地形单元内和全局土壤有机质含量的空间自相关程度，结果见表2。

由表2可知，各地形单元下土壤有机质含量全局Morans I 均大于0，P均小于0.001，因此可以推断县域内不同地形单元下的土壤有机质含量的空间特征并不是随机分布的，而是空间聚集性表现出非常显著的特点，但是聚集效应不同。土壤有机质含量在各类地形单元下有明显的空间正相关关系，其中在坡脚的空间聚集性最高，坡肩的聚集性最低，这也可能是由于坡肩坡度较大，土壤侵蚀严重，造成土壤有机质流失。由以上分析可知，不同地形单元下的土壤有机质含量空间分布规律不同，但是各单元特点都与实际情况相符，这也说明了基于划分地形单元研究空间插值方法的可行性。

3.4 地统计空间分析及Kriging空间插值

利用地统计GS+软件对各地形单元采样点土壤有机质含量进行半方差函数及拟合参数计算[22-24]，结果见表3。

表3是根据最优拟合模型以决定系数大、残差小为依据，比较各模型参数，选择土壤有机质含量半方差函数最优拟合指数参数的[25]。根据块金值的相关定义，原理上采样点之间距离等于0时，半变异函数值也是0，由于存在测量误差以及空间变异特征的影响，土壤样点十分相近时，半变异函数值会发生变化，就会产生块金值效应，因此最小距离尺度下变量的变异特征及误差是由块金值反映的。由表3可知，采样点整体、划分地形单元的有机质块金值不一樣，说明划分地形单元空间插值与全局插值具有差异性。

C是结构方差，表示系统属性区域变量最大空间变异，是由土壤属性本身、海拔和气候等结构性因素引起的变异。由于该区域气候差异小，土壤有机质含量的空间变异特征基本是由于地形差异等结构因素引起的。不同地形单元下结构方差不同，所以地形因素是影响该区域土壤有机质空间变异规律的主要原因。C/（Co+C）即表示块金值与基台值的比，被称为空间相关度，该值的大小是系统变量的空间相关性的大小。根据Cambardella提出的划分依据[26]，块金基台比小于25%，表明空间相关性强;大于75%则变量之间相关性较弱;处于二者之中，说明具有中等强度相关性。由表3可知，沟谷空间相关度大于75%，变异差异小，可能与结构因素密切相关性较小，其他地形单元的空间相关度处于25%～75%之间，呈中等变异。

由表3可知，与全局有机质含量相比，坡肩块金值、基台值都降低，变程也降低，这说明不同地形单元土壤有机质含量的变化相对比较均匀，这也证明对该区域土壤有机质含量进行预测时考虑地形因子的合理性，从另一方面说明了以地形单元对偏关县土壤有机质含量进行Kriging插值结果的有效性和可行性。

3.5 预测精度对比分析

基于划分地形单元验证点Kriging插值土壤有机质含量实测值、预测值精度与全局插值精度分析结果进行比较。根据2种方法所得到的线性回归方程系数，散点分布图见图5，划分地形单元进行Kriging分层插值得到的预测值与实测值的决定系数（0.375 3）明显高于全局Kriging插值法的决定系数（0.132 3）。

表4是对土壤有机质含量空间特征分地形单元和全局预测合理结果的验证。根据均方根误差值越小精度越高，平均预测误差接近于0，满足预测值的无偏性[27]，得出划分地形单元比全局预测精度高，而且插值结果的无偏性好。

3.6 土壤有机质含量的空间分布特征

采用Kriging插值对不同地形单元土壤有机质含量进行空间插值，根据地形单元耕地图斑裁剪各地形单元内的有机质含量空间预测结果分布，并将其拼接得到整个研究区的空间分布见图6-A，而不划分地形单元采用全局采样点直接插值得到的土壤有机质含量空间分布结果见图6-B。

由图6可知，划分地形单元Kriging插值得到耕地土壤有机质含量范围为4.15～14.93 g/kg，全局Kriging得到的含量为4.40～11.72 g/kg。总体来看，研究区2种方法预测结果范围大体一致，空间分布格局基本相同，都呈现出偏东部和南部区域有机质含量较高，东西方向偏低。可见土壤有机质在空间上的特征规律与实际地形相符合[28]，东西方向是沟谷，受人为影响因素比较大，所以有机质含量预测结果偏低，这也进一步验证了地形是引起空间变异存在的必然因素。

从局部细节分布特征来看，2种预测结果之间存在着一定的差异。经过对比分析发现，划分地形单元的预测结果更能展示细节信息，而且细节信息与地形单元的分布规律一致，边界平滑且为渐变过渡;而全局Kriging插值得到的土壤有机质含量空间分布特征比较明晰，但是局部本身属性细节展示不充分，这显然与研究区的实际情况有很大的差异，结果只能反映出整个研究区的大致情况，不能精确反映研究区局部实际情况，说明平滑效应并没有使不同地形单元下的土壤有机质含量差异性变小。比如个别地形单元下耕地由于个别采样点的土壤有机质含量偏高，造成该区域的整体偏高，全局Kriging插值并不会详细体现出高的地方。全局Kriging插值是基于采样点的空间分布特征以及样点之间的空间位置关系来进行插值，没有根据地形单元对采样点进行分类，所以空间插值结果会受距离较近的个别预测样点的影响[29]，而划分地形单元Kriging插值因为按照地形因素与土壤采样点之间的关系进行分区，就可以很好地消除个别样点对插值结果的误差影响。总之，划分地形单元Kriging插值的预测结果精度较高，更能反映土壤有机质含量的空间信息以及局部区域分布情况，而直接Kriging插值忽略了土壤有机质含量空间分布规律平稳过渡性以及误差影响。

4 讨论与结论

对于地形条件复杂的丘陵山区，进行高密度的土壤采样需要大量的人力物力，提前设计好有效土壤样点采样的方案能够充分获取土壤属性的空间信息，还能大大降低成本[30]。所以综合更多个地形因子分类组合精确划分地形单元，研究各个地形单元变异规律，为局部样点更新、最佳采样、以及测土配方精准施肥方案提供地理科学依据。本研究取得的主要研究结果如下：利用DEM数据结合坡度、坡向和地形位置指数综合来划分地形单元，基于地形单元，利用Kriging分层插值对复杂丘陵山区土壤有机质含量进行空间预测，可更详尽地表征出复杂山区的特定地形条件分异特征规律，将各个地形单元预测结果图拼接为全域的空间分布图。设置验证点，并将该方法与全局Kriging插值法得到的预测结果进行比较。划分地形单元与全局预测结果范围和规律一致，但是预测精度高，而且展示了更详细的局部区域空间细节信息。

由于研究区成土因素和母质的影响，土壤有机质含量的变化较小，所以与全局插值相比沟谷和坡肩的变程影响比较大。基于地形单元的土壤有机质含量空间变异规律，是一种综合地形因素的研究方法，更是为复杂的丘陵山区精确插值提供了一种参考。

参考文献：

[1]Massey P A，Creamer R E，Schulte R，et al. The effects of earthworms，botanical diversity and fertiliser type on the vertical distribution of soil nutrients and plant nutrient acquisition[J]. Biology and Fertility of Soils，2013，49（8）：1189-1201.

[2]韩 丹，程先富，谢金红，等. 大别山区江子河流域土壤有机质的空间变异及其影响因素[J]. 土壤学报，2012，49（2）：403-408.

[3]张学雷，冯婉婉，钟国敏.豫中褐土耕地土壤性质空间分异及质量评价[J]. 应用生态学报，2011，22（1）：121-128.

[4]張 泽，吕 新，吕 宁，等. 基于GIS、RS的滴灌棉田土壤养分精确管理分区研究[J]. 农业机械学报，2014，45（7）：125-132.

[5]Zhu H F，Ying Z，Feng N，et al. Relative influence of soil chemistry and topography on soil available micronutrients by structural equation modeling[J]. Journal of Soil Science and Plant Nutrition，2016，16（4）：1038-1051.

[6]Mabit L，Bernard C. Spatial distribution and content of soil organic matter in an agricultural field in eastern Canada，as estimated from geostatistical tools[J]. Earth Surface Processes and Landforms，2010，35（3）：278-283.

[7]付晓婷.基于地貌分区的土壤养分空间变异性分析及采样设计研究——以陕西省蓝田县为例[D]. 西安：西北大学，2012.

[8]韦金丽，王国波，凌子燕.基于高分辨率DEM的地形特征提取与分析[J]. 测绘与空间地理信息，2012，35（1）：33-36.

[9]Yanai J，Mishima A，Funakawa S，et al. Spatial variability of organic matter dynamics in the Semi-Arid croplands of northern Kazakhstan[J]. Soil Science & Plant Nutrition，2005，51（2）：261-269.

[10]朱洪芬，南 锋，徐占军，等. 黄土高原盆地土壤有机质与影响因子的空间多尺度关系[J]. 生态学报，2017，37（24）：8348-8360.

[11]赵业婷，常庆瑞，李志鹏，等. 渭北台塬区耕地土壤有机质与全氮空间特征[J]. 农业机械学报，2014，45（8）：140-148.

[12]耿广坡，高 鹏，吕圣桥，等. 鲁中南山区马蹄峪小流域土壤有机质和全氮空间分布特征[J]. 中国水土保持科学，2011，9（6）：99-105.

[13]朱会义，刘述林，贾绍凤.自然地理要素空间插值的几个问题[J]. 地理研究，2004，23（4）：425-432.

[14]宋 轩，李立东，寇长林，等. 黄水河小流域土壤养分分布及其与地形的关系[J]. 应用生态学报，2011，22（12）：3163-3168.

[15]田瑞云，王玉宽，傅 斌，等. 基于DEM的地形单元多样性指数及其算法[J]. 地理科学进展，2013，32（1）：121-129.

[16]黄 魏，韩宗伟，罗 云，等. 基于地形单元的土壤有机质空间变异研究[J]. 农业机械学报，2015，46（4）：161-167.

[17]田惠文，张欣欣，毕如田，等. 非均匀网格环境下基于耕地质量与集中连片度的基本农田划定[J]. 农业工程学报，2020，36（3）：254-262.

[18]李 超，李文峰.高原耕地土壤养分空间分布与影响因子相关性研究[J]. 土壤通报，2014，45（5）：1113-1118.

[19]Tang C Q. Forest vegetation as related to climate and soil conditions at varying altitudes on a humid subtropical mountain，Mount Emei，Sichuan，China[J]. Ecological Research，2006，21（2）：174-180.

[20]Weiss A D. Topographic position and landforms analysis[C]. San Diego：ESRI User Conference，2001.

[21]黄 婷，周冀衡，李 强，等. 不同海拔高度植烟土壤pH值分布情况及其与土壤养分的关系——以云南省曲靖市为例[J]. 土壤通报，2015，46（1）：105-110.

[22]Matheron G. Principles of geostatistics[J]. Economic Geology，1963，58（8）：1246-1266.

[23]楊 东，程军奇.甘肃黄土高原近52年四季极端降水变化特征

[J]. 中国农学通报，2013，29（11）：180-187.

[24]张素梅，王宗明，张 柏，等. 利用地形和遥感数据预测土壤养分空间分布[J]. 农业工程学报，2010，26（5）：188-194.

[25]张志霞，许明祥，师晨迪，等. 黄土丘陵区不同地貌单元土壤有机碳空间变异的尺度效应[J]. 自然资源学报，2014，29（7）：1173-1184.

[26]Cambardella C A，Moorman T B，Novak J M，et al. Field-scale variability of soil properties in central Iowa soils[J]. Soil Science Society of America Journal，1994，58（5）：1501-1511.

[27]杨之江，陈效民，景 峰，等. 基于GIS和地统计学的稻田土壤养分与重金属空间变异[J]. 应用生态学报，2018，29（6）：1893-1901.

[28]刘 芳，塔西甫拉提·特依拜，依力亚斯江·努尔麦麦提，等. 准东煤炭产业区周边土壤重金属污染与健康风险的空间分布特征[J]. 环境科学，2016，37（12）：4815-4829.

[29]Bhunia G S，Shit P K，Chattopadhyay R. Assessment of spatial variability of soil properties using geostatistical approach of lateritic soil （West Bengal，India）[J]. Annals of Agrarian Science，2018，16（4）：436-443.

[30]Taye M，Simane B，Selsssie Y，et al. Analysis of the spatial variability of soil texture in a tropical highland：the case of the Jema watershed，northwestern highlands of Ethiopia[J]. International Journal of Environmental Research and Public Health，2018，15（9）：1903.