基于地统计学和空间自相关的土壤有机质空间异质性分析及分区*——土地整治视角

左昕弘，赖佳鑫，刘 峰，石孝均

（西南大学资源环境学院，重庆 400715）

0 引言

耕地是土地资源的精华，是保障粮食安全的重要资源[1-2]。土地整治是保障我国粮食安全和生态安全的人类活动，也是耕地质量变化的重要驱动因素[3]。土地整治是指对低效利用、不合理利用、未利用及生产建设活动和自然灾害损毁的土地进行整治、提高土地利用效率的活动[4-5]，农村土地整治更加强调通过田、水、路、林、村综合整治，以增加有效耕地面积并提高耕地质量[6]。土地整治措施的实施对整治区域的土壤理化特性[7]、水资源空间分布[8]、景观环境[9]、生物多样性[10]等有重要影响[3]。如土地平整工程中田块翻耕和表土剥离与回填等工程措施，通过增加土壤透气性和保持土壤耕作层质量来改善土壤理化性质，从而提高耕地质量[11]。2021年4月通过的《中华人民共和国乡村振兴促进法》指出，县级以上地方人民政府应当推进农村土地整治和农用地安全利用，加强农田水利基础设施建设，改善农业生产条件。然而，未结合耕地土壤指标空间异质性的土地整治措施不仅不利于耕地质量改善，而且会增加整治成本[7]。因此，依据耕地质量指标空间异质性进行分区，并采取针对性的土地整治措施，是土地整治项目设计的重要基础，对提高土地整治水平具有重要的意义。

土地整治分区是将实施土地整治工作的整个区域划分为若干整治类型与方向相对一致的子区域。近年来，土地整治分区研究主要集中于研究尺度、评价指标和分区方法等方面。在研究尺度上，以宏观区域（省[12]、市/县[13]、自然区域[14]）居多，微观（如项目区尺度[15]）的研究不足。在评价指标上，从自然质量、利用质量和经济质量[16]逐步拓展到乡村地域功能[17]、生态位[18]、生态系统服务[19]、生命共同体[5]、生态敏感性[20]等方面。但耕地质量指标，特别是土壤指标，仍是土地整治分区需要考虑的基础和核心内容[21]。研究方法方面，在GIS空间分析的基础上，多依据不同土地整治分区的侧重点采用多因素综合评价法[12]和属性特征层次组合法[22]获得分区，当研究对象事前分类方案不明确时，多采用聚类分析法[23]、空间自相关分析法[24]等进行分区。上述研究对提高土地整治效益，提升耕地质量提供了理论依据和技术支持，但因研究尺度较大，对土壤质量指标在土壤剖面上的异质性及水平方向上的空间自相关性少有考虑，从而忽略了评价指标的空间集聚程度及其显著性，导致具有显著空间集聚性的区域可能在分区时被割裂，不利于土地平整工程中表土剥离、移土培肥和土壤改良等关键技术的实施。

土壤有机质（soil organic matter，SOM）作为土壤质量和生产力最重要的指标之一，对改善土壤质地、提高土壤肥力、促进植物生长发育有重要作用[25]，也是基于耕地质量的土地整治分区中多选用的评价指标[26]。文章以土壤有机质为例，引入全局和局部空间自相关分析法，在渝西丘陵区典型土地整治项目区尺度上，探讨基于土壤有机质垂直和水平空间上的异质性和空间自相关性的分区方法及其在土地整治分区中的应用。

1 研究区概况与研究方法

1.1 研究区概况

潼南国家农业科技园区位于渝西紫色土丘陵区。为建设耕地分布集中、基础设施健全、农业发展条件较好的粮食主产区，从2017年底开始进行大规模的基本农田建设。该文研究区为2017年该园区基本农田集中建设区（105°44′39″E～105°45′26″E，30°04′03″N～30°04′53″N），在行政区划上属于重庆市潼南区柏梓镇小岭村，面积约120hm2（图1）。研究区为浅丘带坝地貌类型，海拔250～271 m；属亚热带湿润季风气候，年均气温17.9℃，年均无霜期335d，年均降雨量975.5mm。土壤类型以侏罗纪中统遂宁组厚层泥岩母质发育的红棕紫泥土和红棕紫色水稻土为主，该土壤富含钙质，呈微碱性。主要粮食作物为水稻和玉米。

1.2 土样采集与测定

土壤样品采集于土地整治前（2017 年5 月），依据研究区农用地地块大小及形状特征随机布点采样。研究区以河流（张家堰河）为界呈不同的地貌类型，河流以北为平坝区，农用地地块面积较大，适当减少采样密度，河流以南为丘陵区，受地形限制，农用地地块较小，适当加大采样密度。各采样点按0～20cm、20～40cm、40～60 cm 进行分层采样，分别表示为表层（Upper layer，U）、中层（Middle layer，M）和下层（Lower layer，L）。采样时先去除地表凋落物，每个样点向四周辐射共采集至少5 个分样点，各土层土样分别混合后四分法留取1 kg 土样装入样品袋，共采集样点53 个，分层土样144 个（部分采样点因土层薄，未采集到中、下层土样）。同时用GPS记录样点经纬度和海拔（图1）。

图1 研究区地理位置、土地利用类型及采样点分布

将采集的土壤样品挑出砾石、草根等杂物，自然风干，研磨后在室内进行测定。SOM 含量采用重铬酸钾容量法—外加热法测定[27]。

1.3 研究方法

1.3.1 地统计分析

普通Kriging 法（ordinary kriging）是地统计学中最常用的插值法，是建立在变异函数理论及结构分析基础上，在有限区域内对区域化变量的取值进行无偏最优估计的一种方法。在GS+9.0中进行SOM 的半方差分析，以决定系数R2最大、残差RSS最小确定最优理论插值模型。根据半方差拟合模型和参数，在Arc-GIS 10.2.2中采用普通Kriging法进行SOM的空间插值。

半方差函数的理论模型为[28]：

估算未采样点的理论模型为[29-30]：

式（1）（2）中，ϒ(h)为半方差函数，h为样点空间距离，N(h)为空间距离为h的样点数，r(xi)和r(xi+h)分别为变量r(xi)在空间位置xi和xi+h处的实测值，r(x0)为x0处的指标估计值，λi为第i个样点的权重。

1.3.2 空间自相关分析

空间自相关分析（spatial autocorrelation analysis）用于检验具有空间位置的某要素的观测值是否显著地与其相邻空间点上的观测值相关联[30-31]。包括全局空间自相关和局部空间自相关，通常采用Moran′s I指数表示空间要素的自相关程度。I 的取值范围为-1～1，I＞0 表示变量在空间上呈正相关，表现出空间集聚特征，反之为负相关，表现出空间孤立特征，I=0 表示无空间自相关性。通常用Z值检验自相关显著性，当Z＞1.96或Z＜-1.96时，变量的空间自相关显著，否则不显著，变量呈随机分布。采用GeoDa 1.14进行Moran′s I分析，并在ArcGIS 10.2.2中结合Moran′s I分析结果制作LISA聚类图。

全局Moran′s I指数的表达式为[30]：

局部Moran′s I指数的表达式为[30]：

式（4）中，N为空间变量个数，zi、zj分别为变量在i、j处的值，-z为变量z的均值Z值的计算公式为[32]：

其中：

式（6）中，Wij为zi和zj之间的空间权重函数，i≠j。

1.3.3 基础数据处理

通过目视解译及实地调查获取2017 年研究区土地利用数据。运用Excel 2016 和SPSS 25.0对SOM 进行描述性统计分析和正态分布检验。采用ArcGIS10.2.2的Zonal Statistics工具生成各农用地斑块的SOM含量。

2 结果与分析

2.1 SOM含量统计特征分析

研究区各土层SOM 含量的描述性统计特征如表1 所示。3 层SOM 的含量均值范围为12.86～19.06 g/kg，表现为表层最高，且随土层深度的增加而递减。依据全国第二次土壤普查养分分级标准，各土层SOM 含量的均值均位于土壤养分分级标准的4 级。表层SOM 含量介于1～5 级，中、下层介于2～6 级（表2）。3 层SOM 的空间变异系数范围为40.6%～60.8%，以表层最小，且随土层深度的增加而逐渐增加，均属于中等变异。

表1 不同土层SOM的描述性统计特征

表2 中国土壤养分分级标准（全国第二次土壤普查）

2.2 SOM的空间变异性分析

经典统计学方法仅在一定程度上反映研究区SOM 的样本总体特征及变异状况，不能定量地刻画其随机性和结构性、独立性和相关性，需进一步采用地统计学方法对其空间变异结构进行分析和探讨[33]。经K-S 检验，3 层SOM 除表层呈正态分布外，其余土层经对数转换后呈正态分布，满足地统计学分析要求。该文采用GS+9.0对不同土层SOM 的空间变异性进行分析，用球状、指数和高斯等模型进行拟合，得到模型的相关参数值，当决定系数（R2）接近于1且残差（RSS）较小时，模型拟合度最优。如表3所示，3层土壤中，除中层SOM 的最优拟合模型为球面模型外，其余均为高斯模型。最优拟合模型结构参数中块金值（C0）反映的是最小抽样尺度下变量的测量误差和随机因素引起的变异。基台值（C0+C）表示系统内的总变异量，基台值较大表明系统变量具有较强的空间变异。块基比C0/（C0+C）表示系统变量的空间相关性程度，若比值＞25%，表明系统具有强烈的空间相关性，若比值介于25%～75%，表明系统具有中等的空间相关性，若比值＞75%，表明系统空间相关性较弱，若比值接近于1，表明在整个尺度上具有恒定的变异[34-35]。研究区3 层SOM 的C0/（C0+C）介于28.6%～33.3%，具有中等的空间相关性，表明各土层SOM 的空间变异受结构性因素（成土母质、土壤类型、地形地貌等）和随机性因素（耕作、施肥管理、种植制度等）的共同影响。变程指空间相关性的作用范围，反映区域化变量影响范围的大小，同时也反映变量自相关范围的大小。研究区SOM的变程介于88.33～102.19 m，表明SOM含量在该范围内存在空间相关性。

表3 不同土层SOM的半方差函数理论模型和相关参数

2.3 SOM的空间分布特征分析

为更直观地揭示各土层SOM 在研究区内的空间分布格局，该文基于最优半方差拟合模型所得参数，采用普通Kriging 法在ArcGIS 中进行插值，得到3 层SOM 的空间分布图。如图2 所示，各土层SOM 呈不规则的斑状与块状分布。SOM 含量在表层表现为西北部和东南部的部分区域较高，中部和北部较低；在中、下层表现为由北向南递减。此外，表层和中层SOM 含量介于10～20 g/kg之间的区域面积最大，下层则以＞10 g/kg的区域为主。

2.4 SOM的空间自相关分析

2.4.1 全局空间自相关分析

为计算研究区农用地斑块的全局与局部Moran′s I指数并绘制散点图及LISA聚类图，该文基于SOM空间分布图（图2），采用ArcGIS的Zonal Statistics工具，生成各农用地斑块的SOM平均含量作为基础图层（图3）。

图2 不同土层SOM的空间分布

图3 不同土层农用地斑块SOM平均含量的空间分布

为描述SOM 的空间自相关性随滞后距的变化规律，该文采用GeoDa 软件和Excel 计算并绘制了各土层SOM 在不同滞后距下的空间自相关程度。如图4 所示，各土层SOM 的全局Moran′s I 指数在滞后距位于68.03～1 000.00 m 的范围内大于0，表现为较强的正相关性，表明SOM 在上述范围内呈正空间自相关，存在明显的空间集聚性。当滞后距为68.03 m 时，在保证每个对象都有至少1 个邻居的情况下，各土层SOM 的Moran′s I 指数均达到峰值，此时SOM 的正空间自相关集聚性最强。随着滞后距的增加，Moran′s I 指数下降，且下降趋势逐渐减缓，表明以68.03 m 为阈值构建的空间权重和计算得到的空间自相关结果具有更高的精确度。因此该文选择68.03 m 为阈值对SOM进行空间自相关分析。

图4 不同土层SOM的空间自相关性随滞后距的变化

为探讨各土层SOM 的空间关联程度，该文以农用地斑块为基本单元，对各土层SOM 进行全局空间自相关分析，计算全局Moran′s I 指数并绘制Moran′s I 散点图。如表4 所示，各土层SOM 的Moran′s I指数均大于0.8，呈较强的正空间自相关集聚性，表明各土层SOM 的空间分布表现为相似值的空间集聚，即高值趋向于与高值相邻，低值趋向于与低值相邻，且相邻观测值间相似性较强，空间变异性则相对较弱。各土层SOM的Z值均大于1.96，呈显著空间自相关，表明观测值与周围观测值相似。显著性检验P值均小于0.01，表明各土层SOM含量均表现为较强的正空间自相关集聚态势。

表4 不同土层SOM的空间自相关参数

Moran′s I散点图中第一、三象限为高—高（HH）和低—低（LL）集聚区域，表明指标存在正的空间关联，HH 和LL 分别表示高值斑块被其他高值斑块所包围和低值斑块被其他低值斑块所包围；第二、四象限为低—高（LH）和高—低（HL）集聚区域，表明指标存在负的空间关联，LH 和HL 分别表示低值斑块被其他高值斑块所包围和高值斑块被其他低值斑块所包围。如图5 所示，各土层的SOM 含量主要分布于一、三象限，表现为明显的HH或LL集聚特征，表明研究区SOM的空间自相关性整体呈现为两大趋势，即SOM 含量高值区较集中，形成高值与高值自相关趋势，同时，低值区也较为集中，形成低值与低值自相关趋势。

图5 不同土层SOM的全局Moran's I散点

2.4.2 局部空间自相关分析

全局空间自相关的Moran′s I 指数可判断各指标是否存在空间集聚区及空间离散区，但无法获得其在空间上的分布状况，而局部空间自相关指标结合LISA 聚类图能够揭示各指标在空间上具体的分布位置。因此，为进一步明确各土层SOM 含量空间集聚区和空间离散区的区域分布特征，采用局部空间自相关指标及LISA 聚类图将农用地斑块划分为HH、HL、LH、LL和“不显著”（NN）5种类型并绘制空间分布图。各类型斑块数量与面积分别占总斑块数量与面积的比例如表5所示，SOMM表现为HH型占显著性斑块的比例最大，表明中层SOM主要表现为高—高集聚特征；其余均为LL型所占比例最大，即表、下层SOM主要表现为低—低集聚特征，也表明在土地整治中保护表层养分含量较高地块的重要性。此外，各土层SOM均无LH型，表明研究区无SOM低值被高值所包围的农用地斑块。

表5 各类型斑块数量与面积占比 %

LISA 聚类图表明（图6），表层SOM 的HH 集聚区主要分布于平坝区西北部和东南部丘间谷地，LL 集聚区主要分布于南部丘陵区，部分位于平坝区西部；中、下层SOM 的HH 集聚区则主要集中于北部平坝区，LL集聚区主要分布于南部丘陵区。

图6 不同土层SOM的LISA聚类

2.5 SOM分区

基于SOM 空间自相关分析结果，依据以下原则对耕地进行分区：（1）以表层SOM 空间自相关分布特征为主，中、下层为辅，将3 层SOM 的高值集聚分布、低值集聚分布、高低离散分布、低高离散分布和空间随机分布叠加组合进行分区；（2）3层SOM空间自相关类型叠加后，将斑块面积较小的类别根据表层SOM 空间自相关类型合并到其他斑块面积较大的分区。结果表明，研究区可划分为高SOM 集聚区、表层高SOM集聚区、低SOM集聚区和SOM随机分布区4个分区（表6），空间分布如图7所示。

高SOM集聚区主要分布在水稻种植区，共有31个耕地斑块，面积为4.98hm2，占研究区耕地总面积的5.2%，该区域宜作为优先整治区域（表6 和图7）。由图8 可知，该区域表层SOM 含量显著高于中、下层（P＞0.05），3 层SOM 含量均显著高于其他分区（P＞0.05），是各土层SOM 含量最高的耕地集聚区。该区域自然条件较好，耕地集中连片分布，为耕地利用提供了良好的基础。在土地整治中应重视表土剥离与回填，同时应保护和维持现有自然条件和利用管理水平，避免经济建设对耕地质量的影响。

表层高SOM 集聚区主要分布在东南部丘间谷地，共有31 个耕地斑块，面积为5.43 hm2，占研究区耕地总面积的5.6%，该区域宜作为适度整治区（表6 和图7）。由图8 可知，该区域表层SOM 含量显著高于中、下层（P＞0.05），且表层SOM 含量显著高于除高SOM 集聚区外的其他分区（P＞0.05），中层SOM 含量显著高于低SOM 集聚区（P＞0.05），下层SOM 含量略高于低SOM 集聚区，是表层SOM 含量相对较高的耕地集聚区。在土地整治中应同样重视对表层的保护，同时应加强对中、下层SOM 含量的提升（如有机物料回填、机械改土等）。此外，应保护和维持现有自然条件和利用管理水平，尽量避免经济建设对耕地质量的影响。

图7 SOM分区

低SOM 集聚区主要分布在平坝西部和南部丘陵区坡耕地区域，共有111 个斑块，面积为14.56hm2，占研究区耕地总面积的15.1%，该区域宜作为重点整治区（表6 和图7）。由图8 可知，该区域表层SOM 含量显著高于中、下层（P＞0.05），3层SOM 含量均低于其他分区，且除下层外均表现出显著性（P＞0.05），是各层SOM 含量最低的耕地集聚区。在土地整治中应结合自然条件和利用管理特征进行中长期改良措施（如坡改梯、改良培肥等），同时改善农田生产条件，以提升区域整体耕地质量水平。

SOM 随机分布区空间分布较为广泛，共有327 个斑块，面积为71.53 hm2，占研究区耕地总面积的74.1%，该区域宜作为优化调整区（表6 和图7）。由图8 可知，该区域表层SOM 含量显著高于中、下层（P＞0.05），3 层SOM 含量均处于中等水平。在土地整治中需注意对少量优质耕地的保护，以此为中心优化农田布局，同时应考虑适当调整、优化农业结构。

图8 各分区的SOM含量大写字母表示同一分区下SOM在不同土层间的显著性差异，小写字母表示同一土层下SOM在不同分区间的显著性差异；误差线代表±SD

表6 SOM分区及区域特征

3 讨论

半方差函数分析结果表明，研究区各土层SOM 具有中等的空间相关性，受随机性因素与结构性因素的共同影响，空间自相关范围较小。滞后距较小时，SOM 表现为较强的正空间自相关性，存在明显的空间集聚性（图4 和表4）。这与半方差函数分析的结果较为一致（表3）。但二者对空间自相关衡量的严格程度不同，半方差函数的分析结果包括了正负自相关，而全局Moran′s I 指数将正负自相关进行了区分[36]。从空间分布上看，各土层SOM 呈不规则的斑状与块状分布（图2），LISA 聚类图也表现为围绕HH 和LL的多中心分布格局（图6），这一结果可以指导SOM 在分区时的精准划分，进一步为土地整治分区提供依据。

结合土地利用类型来看，研究区水田所处地势平缓，积水条件下SOM 分解较缓慢，易于积累，且水田的秸秆还田率较高，被微生物转化为有机碳的量较多[37]，使得水田表层SOM 最高且显著高于撂荒地和旱地（P＞0.05）（图9），因此SOM 表层的HH 集聚区主要为水田（图1），分布在北部平坝区和东南部丘间谷地。受施肥和表层土壤腐殖质随下渗水在土体中淋溶、迁移、淀积的影响[38]，中、下层SOM 含量也以水田最高且显著高于其他土地利用类型（P＞0.05），其次为旱地，因此中、下层SOM 的HH 集聚区除水田外还有部分旱地，主要分布在北部平坝区。这是土地利用类型影响土壤养分空间变异及剖面异质性的直接体现。但LL 集聚区土地利用类型表现多样化，且主要分布在南部丘陵区，表明除土地利用类型外，SOM 还受到地形因子、成土母质、土壤类型和秸秆还田方式等其他因素的影响[39]。

图9 不同土地利用类型下不同土层的SOM含量大写字母表示同一土地利用类，SOM在不同土层间的显著性差异，小写字母表示同一土层下SOM在不同土地利用类型间的显著性差异；误差线代表±SD

基于SOM 的剖面异质性和空间自相关性划分的分区，可服务于土地整治中的均质单元划分、表土剥离与回填、移土培肥和土壤改良等工程措施的实施。其中，全局空间自相关明确了SOM 的空间相关性与规律性，局部空间自相关进一步描述了SOM 的空间自相关程度，Moran′s I 散点图和LISA 聚类图则展示了上述空间自相关在空间上的具体分布格局[40]。与依据SOM 含量等级划分的分区相比，基于3 层SOM 的LISA 聚类结果的叠加分区，能将SOM 在垂直和水平方向上的分布情况与空间集聚性结合，从而提高评价指标分区及基于此的土地整治分区的合理性，更有利于针对性地采取土地整治措施。

4 结论

在渝西典型紫色土丘陵区，0～60 cm 土层中的SOM 含量总体较低，受结构性和随机性因素的共同影响，各土层SOM 均表现为显著的正空间自相关，且空间自相关范围较小。在水平空间分布上，各土层SOM 呈不规则的斑状与块状分布，LISA 聚类图也表现为围绕HH 和LL 的多中心分布格局。土地利用类型对SOM 的空间分布存在显著影响，HH 集聚区主要为水田，多分布于平坝和东南部丘间谷地，LL 集聚区土地利用类型则表现多样化，主要分布于南部丘陵区。以表层SOM 空间自相关分布特征为主，中、下层为辅，可将研究区划分为高SOM 集聚区、表层高SOM 集聚区、低SOM 集聚区和SOM 随机分布区4 个分区。各分区体现了SOM 含量在剖面和水平上高值与低值的集聚性特征，可作为土地整治分区，尤其是土地平整工程中的表土剥离、移土培肥和土壤改良等关键技术实施的重要参考。