土壤有机质是土壤的重要组成部分，对改善土壤理化性质、调节土壤养分循环具有重要作用，是评价土地生产力的关键指标^[1]。掌握土壤有机质的空间分布特征及其关联因素，可为制定合理的农业种植规划及农田管理措施提供依据，是实现农田土壤可持续利用的前提^[2]。全国第二次土壤普查以来，我国大部分地区农田土壤有机质含量均呈现上升的趋势^[3]。了解现阶段农田土壤有机质含量的现状及空间分布特征，对下一步开展土壤有机质管理、土壤培肥等具有重要意义。为此，学者们开展了大量研究工作^[2-4]。马利芳等^[5]以新疆阜康市为研究区，综合运用经典统计学和地统计学方法，研究不同深度土壤盐分和有机质的空间变异特征。刘忆莹等^[6]采用地理空间插值方法，探讨了东北典型黑土区耕地有机质的空间分布规律。已有的研究多是基于地统计学方法评价土壤有机质的空间变异规律，并进一步采用不同的插值方法进行区域可视化。

黄淮海平原作为我国重要的农业生产基地之一，针对该地区土壤有机质时空变化特征的研究已有诸多报道。陶雯等^[7]对华北平原引黄灌区农业土壤资源质量的研究结果表明，灌区农田土壤有机质含量存在明显时空分布差异性。胡景辉等^[8]探讨了河北山前平原高产类型区土壤养分状况，结果表明土壤有机质含量由第二次土壤普查的10.45 g·kg^-1提升至2014—2015年的19.50 g·kg^-1。贾良良等^[9]通过分析河北省2009—2014年测土配方施肥数据，发现河北全省土壤有机质含量均值为17 g·kg^-1，且自西北向东南方向逐渐降低。然而，前人的研究多是基于田块尺度或者小区域尺度展开，并且单纯注重描述土壤有机质含量的时间变化或空间分布^[10-12]。此外，调研取样点范围小、取样时期不统一等是目前研究的普遍现状。在黄淮海平原全区域尺度内开展农田土壤有机质含量现状调查，并深入分析其空间分布特征，对该地区农田土壤碳固持和储量评估、在宏观尺度上把握土壤有机质含量的动态分布特征、因地制宜制定相应的培肥管理措施具有重要意义。

鉴于此，本研究依托2017年黄淮海旱作区土壤肥力调查项目，从空间自相关的角度探讨了该区域旱作农田土壤有机质含量的空间分布特征，并进一步分析了土壤有机质与黏粒含量和容重的空间自相关关系，基于划分的空间自相关类型，浅析了每种空间结构可能出现的原因，以期为确定合理的宏观土壤养分管理单元及土地利用规划等提供一定的参考。

本研究中黄淮海平原旱作农田的确定：基于2015年黄淮海地区土地利用现状图，根据坡度分级标准，将地形坡度<5°、1 km²网格内旱地占耕地40%以上区域、并考虑已有野外布点情况的区域作为旱作农田（图 1）。研究区位于32°13′~41°06′N，112°24′~119°85′E之间，涵盖北京、天津、河北、山东、河南以及安徽的274个县市，总面积约270 000 km²。该区域地处暖温带大陆性季风气候带，雨热同期特点明显，春冬寒冷干燥，夏季高温多雨，年降水量在500~900 mm之间，年际降水差异较大，年内降水季节分配不均，多集中于夏季，约占全年降水的70%，区域内土壤类型以潮土、褐土、砂姜黑土为主，农业种植制度主要为冬小麦-夏玉米轮作，是我国重要的农业生产区域之一^[13-14]。

图 1 研究区高程和样点分布图 Fig. 1 Map of elevation and sample sites distribution in the study area

本研究土壤样点数据来源于2017年旱作区土壤肥力调查项目。土壤采样点布设原则和方法如下：通过采用网格布点与分层抽样相结合，在研究区以15 km × 15 km网格布点并抽取旱地斑块（含水浇地），结合土壤亚类进行抽样，同时考虑种植体系、分布面积和集中程度等因素，主要在潮土、褐土、砂姜黑土等土类布设，且保证每种土壤亚类均有布设，根据上述规则，黄淮海旱作农田共设置267个土壤采样点（图 1）。土壤样品采集时间为2017年5月，样地均选择作物长势良好的小麦地块，采样深度为0~40 cm，每10 cm一个层次，每点取原状土和扰动土各一种，每种设3个重复。采集的土样利用重铬酸钾外加热法测定土壤有机质含量，激光粒度仪法测定土壤机械组成，环刀法测定土壤容重。

空间自相关是指地理对象的某一属性值在相邻空间位置处取值之间的相互关系，是对该属性空间聚集或分散程度的度量。空间自相关现象广泛存在于地理空间数据中^[15]，因此，进行空间自相关分析对深入理解地理对象的空间分布特征具有重要意义。Moran's I指数在地理对象的空间自相关性研究中具有广泛应用，其不仅能在总体水平上对空间自相关现象做出评价，而且能够对地理属性值的聚集和离散程度作进一步的区域可视化。目前，Moran's I指数在耕地质量保护^[16]、土壤重金属污染^[17]、土壤养分与肥力^[18]、土壤有机质^[19]等研究中得到普遍应用。空间自相关分为全局空间自相关和局部空间自相关。全局空间自相关对属性值的整体自相关程度进行描述，反映空间邻域单元属性值的总体相似程度^[16]。本研究中全局空间自相关采用全局Moran's I指数描述^[20]：

式中，N为采样点总数目，W_ij为空间权重，X_i和X_j分别为i处和j处的有机质含量，X为有机质含量均值；Z为标准化统计量值，E（I）为观测变量自相关性的期望值，var（I）为方差。Moran's I指数的取值范围为[-1, 1]，I> 0代表属性值呈现空间正相关，趋于空间聚集；I<0代表属性值呈现空间负相关，趋于空间分散；I= 0代表属性值趋于空间随机分布。研究采用Z值进行显著性检验，当Z≥1.96或Z≤-1.96时，表示空间相关性显著（P<0.05）。局部空间自相关用于进一步度量每个地域单元与其临近空间单元的属性值之间的相似性和相关性。本研究中局部空间自相关采用局部Moran's I指数进行描述^[20]：

式中，等号右边各项含义同式（1）。将局部Moran's I指数进行区域可视化，通过构建LISA（Local Indicators of Spatial Association）集聚图，揭示属性值的不同空间分布类型。本研究采用LISA图将土壤有机质的空间自相关类型分为五种结构^[21]：HH（高-高型）、LL（低-低型）、HL（高-低型）、LH（低-高型）以及NN（非显著型）。其中，HH和LL为正相关类型，表示土壤有机质空间聚合程度较高，即土壤有机质含量高/低的地区周围含量也高/低，HL和LH为负相关类型，表示土壤有机质呈空间离散特性，即土壤有机质含量高/低的地区周围含量低/高，NN为非显著型，表示土壤有机质无明显聚集或分散特征。此外，研究进一步运用双变量局部空间自相关分析方法^[22]，探究土壤有机质含量与土壤黏粒含量/容重的空间自相关关系，类似单变量空间分析，双变量局部空间自相关采用双变量局部Moran's I指数进行描述，计算公式如下：

式中，Y_i为i处的土壤黏粒含量/土壤容重，Y为土壤黏粒含量/土壤容重均值，其余各项含义同式（1）、式（4）。研究运用ArcGIS10.2软件进行地理数据空间结果展示，运用GeoDa1.12软件进行空间自相关分析。

黄淮海平原旱作农田不同土层有机质含量的描述统计特征如表 1所示。由表 1可知，不同土层有机质含量的区域均值分别为20.11±6.46（0~10 cm），14.76±5.11（10~20 cm），9.96±4.14（20~30 cm），8.03±3.45（30~40 cm）g·kg^-1，变异系数分别为32%、35%、42%、43%。有机质含量随土层深度的增加而减少，变异系数则不断增大，但均属中等变异范围^[23]。表层土壤（0~10 cm）有机质含量的极差最大，为33.99 g·kg^-1，底层土壤（30~40 cm）最小，为21.81 g·kg^-1。

表 1 Table 1

表 1 不同土层有机质含量的统计特征 Table 1 Statistical characteristics of soil organic matter content in different soil layers 土层 Soil layers /cm 最小值 Minimum /（g·kg-1） 25%分位数 25% quantile /（g·kg-1） 平均值 Mean /（g·kg-1） 中位数 Median /（g·kg-1） 75%分位数 75% quantile /（g·kg-1） 最大值 Maximum /（g·kg-1） 标准差 SD /（g·kg-1） 变异系数 CV/% 0~10 4.96 15.65 20.11 19.87 24.08 38.95 6.46 32 10~20 3.91 11.44 14.76 14.18 17.87 35.14 5.11 35 20~30 2.02 6.98 9.96 9.39 12.11 26.79 4.14 42 30~40 1.04 5.61 8.03 7.58 9.93 22.85 3.45 43

表 1 不同土层有机质含量的统计特征 Table 1 Statistical characteristics of soil organic matter content in different soil layers

不同土层有机质含量的空间分布如图 2所示。0~10 cm土层有机质含量多数在15 g·kg^-1以上，高值区域主要分布在河北西部、山东、河南南部和安徽北部地区，有机质含量在20~30 g·kg^-1之间；10~20 cm土层有机质含量多数在10~20 g·kg^-1之间，高值区域分布类似0~10 cm土层；20~30 cm、30~40 cm土层有机质含量多数在15 g·kg^-1以下，且高值区域主要分布在研究区南部，多数样点在10~15 g·kg^-1之间，低值区域主要分布在研究区北部，多数样点在10 g·kg^-1以下。20~30 cm、30~40 cm土层有机质含量总体呈现出“南高北低”的特点。

图 2 不同土层有机质含量的空间分布 Fig. 2 Spatial distribution of soil organic matter content in different soil layers

根据全国第二次土壤普查养分分级标准可知，研究区0~10 cm土层有机质含量总体处于三级水平（20~30 g·kg^-1），10~20 cm土层有机质含量总体处于四级水平（10~20 g·kg^-1），20~30 cm、30~40 cm土层则处于五级水平（6~10 g·kg^-1）。从空间分布来看，不同土层之间存在较大差异。0~10 cm土层中，河北西部、山东、河南南部和安徽北部地区有机质含量均达到了三级水平，上述区域在10~20 cm土层的有机质含量多数在15~20 g·kg^-1之间，接近三级水平。从地理位置看，河北西部地处太行山山前平原地带，水资源条件优越，是我国高集约农业生产区域之一，为获得作物高产，大量投入水肥和有机物料^[24]，同时实施秸秆还田^[8]，是导致土壤有机质含量较高的重要原因。研究区山东县市地处引黄灌区，水资源条件较好，主要土壤类型为潮土，农业资源条件优越，完善的农田管理措施使土壤有机质含量处在一个较高的水平。位于引黄灌区的中国科学院禹城试验站长期定位试验结果表明，长期平衡施肥显著增加了该地区土壤有机质含量^[25]，2016年德州市土壤有机质含量为16.3 g·kg^-1，与第二次土壤普查相比，增加了5.47 g·kg^{-1 [7]}。0~20 cm土层中，河南南部和安徽北部的部分区域土壤有机质含量同样维持在一个较高的水平，该地区土壤类型主要为砂姜黑土，是一种典型的中低产土壤，该区域土壤有机质含量相对较高，与近年来针对砂姜黑土改良的增施有机肥、秸秆还田等农田管理措施有一定的关系^[26]。

20~30 cm、30~40 cm土层中，土壤有机质含量多数降到了15 g·kg^-1以下，处在四级水平或五级水平的标准。全国第二次土壤普查以来，我国大部分地区农田土壤有机质含量均呈上升趋势^[3]，耕作措施、施肥方式、灌溉、秸秆还田等农田管理措施对土壤有机质含量的变化起决定性的作用。有研究指出，长期均衡施肥，可显著增加土壤有机质储量，但这种增加主要集中在0~20 cm的土层^[25]。化肥通过使作物生长繁茂，根茬、枝叶等残留量增多，间接增加土壤有机质含量^[27]。然而，黄淮海平原种植制度比较单一，尤其是河北和山东地区，连年的冬小麦-夏玉米轮作，使得根系分布和下扎深度年际变化小^[28]，连年的玉米季免耕、小麦季旋耕导致土壤耕层变浅、犁底层变厚变硬^[29]，同时，农业机械化程度高，土壤亚表层长期受机械耕作压实^[30]，减弱了根系对下层土壤的穿透，从而使根茬在深层土壤的残留减少。此外，现行的旋耕耕作制度深度一般在15 cm左右^[31]，秸秆旋耕还田、免耕还田^[26]也使得表层土壤有机质含量高于深层土壤。

不同土层有机质含量的空间自相关分析结果表明：0~10 cm土层有机质含量的全局空间自相关特征并不显著（-1.96<Z=-1.9099<1.96）；10~20 cm、20~30 cm、30~40 cm土层有机质含量的全局Moran's I值分别为0.0388、0.0389、0.0591，且相应的标准化Z值均大于检验临界值1.96，表明各层次均呈现出显著的空间正相关特性（P<0.05），土壤有机质含量有空间聚集的趋势。

图 3揭示了在95%的置信水平下，0~10 cm土层大部分土壤样点表现出随机分布特点，并无显著的局部空间自相关特征，这与全局空间自相关结果一致。10~20 cm、20~30 cm、30~40 cm土层中，HL、LL型样点主要分布在河北平原和鲁西北地区，且以LL型为主，表明这些区域除少数样点（HL型）外，有机质含量呈现低值聚集的空间分布结构。HH和LH型在10~20 cm土层中集中分布在河南中南部区域，在20~30 cm土层中主要分布在河南南部，在30~40 cm土层中则主要分布在河南南部和安徽北部区域，不同土层中的具体分布差异较大。

注：HH（高-高型）、LH（低-高型）、HL（高-低型）、LL（低-低型）、NN（非显著型）。   Note: HH(High-High), LH(Low-High), HL(High-Low), LL(Low-Low), NN(Nonsignificant) 图 3 不同土层有机质含量局部空间自相关LISA（Local Indicators of Spatial Association）集聚图 Fig. 3 LISA(Local Indicators of Spatial Association)maps of local spatial autocorrelation of soil organic matter content in different soil layers

从空间自相关的分析结果来看，0~10 cm土层有机质含量在区域分布上并无地理位置上的相关性。由于表层土壤（0~10 cm）更容易受到人类活动和自然因素的影响，其有机质含量并无显著空间自相关性，反映了表层土壤有机质含量在空间结构特征上的随机性。10~20 cm、20~30 cm、30~40 cm土层中，有机质含量均呈现出显著的空间正相关特征，并且随着土层深度的增加，正相关性增强，空间结构趋于明显。陈思明等^[19]的研究结果表明，随着土层深度的增加，结构性因素增强，随机性因素减少，土壤有机质含量的空间分布趋向规则。此外，10~20 cm、20~30 cm、30~40 cm土层中，LL/HL型样点集中分布在河北平原和鲁西北地区，且在不同土层的分布具有一定的相似性，表明这些地区土壤有机质含量在空间水平方向以及土壤垂直分层上均具有一定的结构稳定性。农业种植制度单一（冬小麦-夏玉米轮作），以及所带来的连年固定的农业管理措施，是造成上述分布特征的重要原因。相反，HH/LH型样点在不同土层的空间分布差异较大。研究区南部的河南、安徽地区土壤类型包括潮土、砂姜黑土、黄褐土、褐土等，丰富的土壤类型是导致有机质含量空间差异较大的重要原因之一。此外，种植结构多元化、农业管理措施的地域差异性也会对土壤有机质含量分布产生重要影响。

LISA集聚图中，HL和LH型样点的存在反映了局部地区土壤有机质含量的空间异向性。除了土壤母质、土壤类型的影响，黄淮海地区耕地地块破碎、小农户种植管理比重大，导致农田管理措施具有较大的地域差异性，从而使得土壤有机质含量出现“低值环绕高值”（HL型）以及“高值环绕低值”（LH型）的空间分布结构。在实际农业生产中，HL型样点土壤有机质含量高于周边地区的原因可做进一步分析，以期为周边低有机质含量区培肥土壤、提升地力提供参考；LH型则可以借鉴周边高有机质含量农田的管理措施等，促使向HH型转化。

基于空间自相关理论，进一步探讨有机质含量与黏粒含量/土壤容重的双变量空间自相关特征，相应结果解释如下：全局Moran's I值表示土壤有机质与黏粒含量/土壤容重的全局空间自相关性；局部空间自相关中，HH表示土壤有机质含量高的区域，其周围土壤样点的黏粒含量/土壤容重也高，其他4种自相关类型的含义与此类似^[22]。与传统的数理相关分析相比，双变量空间自相关分析对识别两个地理要素在空间上的同向/异向分布特征，揭示其空间分布结构具有独特优势。

各个土层有机质含量和黏粒含量之间均呈现出显著的数理相关关系（P<0.01），相关系数分别为0.39、0.47、0.43、0.60。土壤有机质和黏粒含量的全局空间自相关结果表明：除0~10 cm土层两者无显著空间自相关关系外，其余层次有机质含量与黏粒含量均呈现显著的全局空间正相关性，表明有机质与黏粒含量的空间分布具有一定的同向性。有机质和黏粒含量的LISA集聚图（图 4）进一步揭示了两者的空间分布特征：不同土层中，LL型和HL型主要分布在河北平原和山东西北部地区，尤其在0~20 cm土层中，HL型集中分布在太行山山前平原以及山东引黄灌区一带；HH型和LH型则主要分布在河南东南部和安徽北部地区。

图 4 不同土层有机质含量和黏粒含量的局部空间自相关LISA（Local Indicators of Spatial Association）集聚图。注：HH（高-高型）、LH（低-高型）、HL（高-低型）、LL（低-低型）、NN（非显著型） Fig. 4 LISA(Local Indicators of Spatial Association)maps of local spatial autocorrelation of soil organic matter content and clay content. Note: HH(High-High), LH(Low-High), HL(High-Low), LL(Low-Low), NN(Nonsignificant)

一般而言，土壤中黏粒含量越多，有机质越容易被吸附和积累，两者之间具有较大的关联性^[32]。本研究中，土壤有机质含量与黏粒含量具有较好的数理相关关系；LISA集聚图中，HH/LL型土壤样点在空间上印证了上述理论。然而也有研究表明，土壤黏粒、粉粒含量与土壤有机质呈负相关关系，而土壤砂粒与土壤有机质呈正相关关系^[33]，造成这种结果差异的原因可能是研究区域的不同，以及所带来的农田管理措施的差异性。本研究中，研究区北部的河北平原和鲁西北地区土壤有机质含量较低与土壤黏粒含量偏低有较大的关系（LL型主导），HL型样点的存在则表明局部地区土壤有机质管理的优越性，使之维持在一个较高的水平。研究区南部的LH型地区土壤黏粒含量较高，相应的有机质含量却处在一个较低的水平，这些地区应被视为农业生产中的“有机质含量异常点”，不适当的人类活动，如频繁耕作、秸秆移除或焚烧等^[25]，一方面会加速土壤有机质的消耗，另一方面使其得不到补充，从而使有机质含量相对偏低。

土壤容重是土壤物理性状的一项重要指标，反映了土壤的紧实度、孔隙状况和结构特征^[34]。黄淮海旱作农田土壤有机质和容重之间的统计学相关系数表明两者之间并未有显著的相关关系。土壤有机质与容重的全局空间自相关结果表明：总体水平上，各个土层有机质含量与容重均呈现显著的空间正相关性（P<0.05）。有机质与容重的LISA集聚图（图 5）中：0~10 cm、10~20 cm土层HH/LH型主要分布在河南中南部和安徽北部地区；LL/HL型主要分布在河北平原和鲁西北地区，且HL型主要分布在太行山山前平原及引黄灌区一带。20~30 cm土层中，各自相关类型的空间分布较表层变化较大，且NN型的样点增多。30~40 cm土层中，LL/HL型主要分布在河北和山东西北部地区，HH/LH型则主要分布在河南地区。20~40 cm土层各自相关类型的空间分布差异较大。

图 5 不同土层有机质和容重的局部空间自相关LISA（Local Indicators of Spatial Association）集聚图。注：HH（高-高型）、LH（低-高型）、HL（高-低型）、LL（低-低型）、NN（非显著型） Fig. 5 LISA(Local Indicators of Spatial Association)maps of local spatial autocorrelation of soil organic matter and bulk density in different soil layers. Note: HH(High-High), LH(Low-High), HL(High-Low), LL(Low-Low), NN(Nonsignificant)

通常情况下，农田土壤有机质含量越多，越容易促进土壤团聚体的形成，从而使土壤中大小孔隙增多，容重降低。然而本研究中，土壤有机质含量与容重之间并没有显著的数理相关关系。有研究表明，耕作对容重的影响较大，耕作能改善表层土壤孔隙状况，使表层土壤容重减小^[29]，免耕则会使表层土壤容重增加^[35]。冬小麦-夏玉米轮作、夏玉米季免耕冬小麦季旋耕是黄淮海旱作农田普遍的种植制度和耕作制度^[28]。连年固定的种植模式使得作物根系分布和下扎深度年际变化较小，对下层土壤的穿透减少，减弱对底层土壤的疏松作用^[27]。此外，连年的旋耕、农机具的压实等，促使土壤亚表层变厚变硬，容重增大^[28，36]。农业措施对土壤容重、有机质含量带来了垂直分布上的影响，而各地区在具体耕作管理上的不同，则造成了两者的地域差异性，从而使两者之间关系复杂化。本研究中，土壤有机质和容重的LISA集聚图揭示了两者在空间上的同向/异向分布特征，为区域土壤有机质调控提供了一定参考。例如，对HH型地区来说，在维持现有有机质含量的前提下，可考虑深松等耕作方式，疏松深层土壤，增大孔隙度，减小容重，为作物根系生长创造良好的土壤物理环境。LH型地区应在注重提高土壤有机质含量的同时，适当降低土壤容重，比如秸秆深耕还田等。LL型地区则应在维持或降低现有土壤容重的情况下，通过增施有机肥等方法提高土壤有机质含量。HL型地区可认为是相对比较理想的农业生产示范区，例如太行山山前平原一带以及山东引黄灌区一带。需要指出的是，土壤有机质培肥、作物根系生长最适容重等耕层特征的构建，需要在现有实际情况的基础上进行改良和完善，并结合各地区农业生产实际情况，因地制宜实施农业管理。

随着土层深度的增加，有机质含量的空间分布结构趋于明显。有机质含量与黏粒含量之间存在显著的数理相关关系；两者的空间分布结构中，HL/LH（高-低/低-高）型反映了在土壤黏粒含量影响的基础上，人类活动等对土壤有机质含量的增加或减少作用。有机质含量与土壤容重之间无显著的数理相关关系；两者的空间分布结构中，HL（高-低）型样点“有机质含量较高，同时土壤容重较低”，是相对理想的分布结构；LH（低-高）型样点“有机质含量低，同时土壤容重高”，属于相对不理想的空间分布结构。有机质与黏粒含量的空间分布结构更明显；0~10 cm、10~20 cm土层较20~30 cm、30~40 cm土层的空间分布结构更明显；农业管理措施及其地域差异性是造成上述分布结构在空间水平方向及土壤垂直分层方向上差异化的重要原因。