土壤有机质(SOM)是泛指以多种形式和状态存在于土壤中的各种含碳有机化合物，是土壤的重要组成部分，与土壤养分、 土壤理化性质及土壤侵蚀有着重要关系，对保持土壤肥力具有积极作用，是反映土壤质量的一项重要指标^{[1, 2, 3, 4]}. 进行SOM空间变异及分布格局研究，并对其影响因素做进一步地探讨，对掌握区域SOM空间分布现状，进行SOM精准培育及土壤肥力的有效调控与提升意义重大. 地统计学被证明是土壤属性空间分布特征及其变异规律研究中最为有效的方法之一，利用地统计学的方法，结合现代“3S”技术强大的空间信息获取与分析功能，对SOM的空间变异和分布格局及其影响因素展开研究已成为当今土壤学领域研究的热点问题^{[3, 4, 5, 6, 7]}. 当前研究多见于南方紫土丘陵、 西北黄土区等区域^{[8, 9, 10]}，且多停留在以行政区划或单一地类地块为单元的层面^{[11, 12, 13]}，容易忽略区域SOM形成和演变过程中自然要素的影响，一定程度上也制约了SOM研究的客观真实性. 流域作为一种集社会经济和自然过程于一体的复合型系统，是基于自然因素形成的典型区域，以流域为单元进行SOM空间特征研究可以很好地揭示流域内社会经济及自然因素与SOM形成及其空间特征的内在关系，具有一定的科学研究价值，而目前基于小流域尺度进行黑土区SOM空间异质性和分布格局及其影响因素分析的相关研究还较为少见^{[5, 14, 15, 16, 17, 18, 19]}.

我国东北黑土区以土壤肥沃、 作物产量高著称，是我国重要的商品粮、 畜牧基地. 大量研究表明，近几十年来东北黑土区在大量产出的同时，黑土区SOM正以不同的速度不断流失，土壤肥力逐年下降，黑土质量状况不容乐观，保护与拯救黑土资源的形势十分紧迫^{[6, 7, 20]}. 如何科学准确地掌握黑土区SOM的空间变异和现状格局及其影响因素是进行黑土区SOM保护的基础工作. 本研究综合地统计学和“3S”技术方法，对东北黑土丘陵区海沟河小流域土壤表层(0～20 cm)中有机质的空间变异和分布格局进行研究，并对流域内影响SOM空间分布的主控因素进行探讨，进一步把握黑土区SOM空间分布特征及主控因子的识别，以期为区域内SOM的精准修复与培育、 农业生产增产等提供基础资料.

海沟河是阿什河的支流，经阿什河在哈尔滨东郊汇入松花江. 海沟河小流域属于松花江流域的一部分，地理位置为126°55′45″～127°10′05″ E，北纬45°34′18″～45°40′50″ N，面积为119.76 km²，行政区划以阿城区的料甸满族乡为主(图 1). 流域海拔高度为54～366 m，东高西低呈阶梯状分布，东部地势起伏较大. 土壤类型主要是黑土，属于典型的黑土丘陵区. 气候类型属于中温带大陆性季风气候，冬季干旱漫长，夏季温热湿润，多年平均气温3.9℃，年降水量600～800 mm，多集中在每年的4～10月，占全年降雨量的90%. 流域内土地资源开发历史较早，农业生产活动频繁，西部以农业生产为主，耕地面积占总面积的70.97%，东部低山丘陵区多林业分布，森林覆盖面积占总面积的24%. 交通便利，受周边城市发展的辐射影响，流域内社会经济状况发展迅速，但城市发展的同时也给流域内生态环境带来一定的负面影响，土壤退化、 水源污染、 生态恶化等环境问题凸显.

图 1

Fig. 1

图 1 研究区地理位置及样点分布示意 Fig. 1 Sketch map of geography location for the study area and distribution of the sampling sites

本研究采用Google Earth高分辨率遥感影像，通过选取地面控制点进行几何校正、 配准，在ArcGIS 10.0平台中采取人工目视判读法进行室内解译，获取海沟河流域2014年土地利用数据，通过野外定点采样调研进行验证，室内解译准确率达90%以上，满足研究需要. 然后以研究区2014年土地利用数据为基础，采用美国陆地资源卫星Landsat5影像数据，结合研究区国土资源资料，进一步提取研究区1992年土地利用数据.

以研究区数字高程模型(30 m×30 m)、 土地利用数据为辅助信息，根据地形、 土地利用方式及距居民点、 交通主干道、 河流水系的距离等要素，在室内进行网格布点. 为避免作物生长、 耕作施肥等对SOM含量的影响，于2014年10月农作物收获后，通过手持GPS进行空间定位，与室内预布点位尽可能一致，进行野外采样，样点主要集中于流域西部耕地区(图 1). 采样过程中，每个样点按照不同方位等距离采样5次，采集土壤表层0-20 cm耕层土壤，用“四分法”取混合均匀的样品约1 kg，留足实验用量，并分别记录采样点坡位及周边环境等相关信息，共采集120个土壤样品. 采样完成后将土壤样品带回实验室，经自然风干，用木棒压碎挑出砾石草根等杂物，研磨过筛，采用重铬酸钾容量法测定每个土壤的有机质含量.

本研究在考虑采样及指标测试误差的基础上，采用格拉布斯(Grubbs)法对所得数据进行最严格的异常值识别和处理，剔除异常数值，后续相关计算分析均采用异常值处理后的数据进行^[21]. 检验数据是否呈正态分布是进行克里格空间分析的基础^[22]，本研究数据正态分布性检验在SPSS 19.0中完成，用K-S检验法进行非参数检验(α=0.05)； 空间自相关分析、 半方差函数的计算及理论模型的拟合等相关参数的计算均采用地统计学软件GS+9.0进行. 其中，半方差函数是SOM空间异质性研究的最有效方法，是SOM进行克里格空间插值的基础，其表达式如下^[21]：

克里格(Kriging)插值法在土壤属性的空间分布中得到广泛应用，并取得了很好的效果. 本研究根据半方差函数的拟合模型，在ArcGIS 10.0平台中进行研究区SOM的普通克里格(Ordinary Kriging)空间插值，然后在普通克里格插值的基础上，选取与土壤属性空间分布显著相关的地形因子作为协变量，对海沟河小流域SOM进行协同克里格(Co-Kriging)插值，并与普通克里格插值进行精度比较. 在精度验证方面，本研究综合克里格模型及其参数的适宜程度，以均方根预测误差最小，平均标准误差最接近于均方根预测误差、 标准均方根预测误差最接近于1为最优插值标准^{[21, 23]}，进行插值结果精度的验证.

通过对研究区SOM测试数据进行经典统计学分析，得到海沟河流域SOM描述性统计特征(表 1). 由表 1可以看出：流域内SOM含量平均水平为19.35 g ·kg^-1，与黑土区相关研究比较处于较低水平，且有机质含量的极大值、 极小值也均处于较低水平^{[23, 24]}. 研究区SOM变异系数为51.11%，结合一般的CV评估标准(当CV<0.1时，为弱变异程度，当0.11时为强变异程度)^{[18, 21]}，研究区SOM空间变异程度为中等变异程度，说明在海沟河小流域SOM的空间变异受人类活动、 耕作措施、 土地利用等随机因素的影响较为明显. 同时通过K-S验证法进行非参数检验(α=0.05)，表明SOM不符合正态分布，故对其进行对数转化后以符合正态分布.

表 1(Table 1)

表 1 SOM描述性统计特征值 Table 1 Descriptive statistical characteristic values of SOM 项目 均值/g ·kg-1 极小值/g ·kg-1 极大值/g ·kg-1 标准差 偏度 峰度 变异系数/% K-S检验 转换前 19.35 6.29 39.94 9.89 0.92 -0.31 51.11 0.011 转换后 1.23 0.80 1.60 0.21 0.07 -0.51 17.07 0.357

表 1 SOM描述性统计特征值 Table 1 Descriptive statistical characteristic values of SOM

空间自相关性是检验某一类要素是否显著地与其相邻要素相关联的指标，是进行地统计插值的基础. 大量研究也已表明地理现象具有普遍的空间自相关特性，一般通过Moran's I指数的大小表示^{[21, 25]}.

本研究中通过GS+9.0对海沟河小流域SOM整体及其在各方向上的空间相关性进行分析，结果见图 2. 从中可以看出：海沟河小流域SOM空间自相关性整体上随着间隔距离的增加先呈正相关性，后转变为负相关性(间隔距离为3 200 m处)，在约7 000 m处又变化为正相关，变化比较复杂[图 2(a)]，Moran's I主要分布在-0.167-0.126之间，其变化随间隔距离的变化而出现较大变化，整体自相关性较弱. 通过各向自相关性分析[图 2(b)]发现，在0°(南-北)方向上，Moran's I值随间隔距离的变化在1 000 m处大于0，在横坐标轴上方呈现出一系列的波动，在间隔距离为6 400 m处有最大值出现，表现出较强的空间正自相关性； 在45°(东北-西南)方向则以负相关为主，在小范围内出现正相关性(间隔距离为1 200 m处)，整体规律性较差； 在90°(东-西)方向，Moran's I值的变化趋势与图 2(a)类似，呈现先正相关后负相关之势，可以反映整体自相关性特征； 在135°(东南-西北)方向上，当间隔距离小于3 000 m时，表现较强的正相关性，随着距离的增加Moran's I值以负数为主，表现出较强的空间负自相关性. 综观以上分析，海沟河小流域SOM含 量存在空间自相关性，且在不同方向上具有差异性，可以进行地统计插值.

图 2

Fig. 2

图 2 SOM空间自相关性 Fig. 2 Spatial autocorrelation charts of SOM

对海沟河小流域SOM进行半方差函数计算与模型拟合，得到研究区SOM半方差函数模型[图 3(a)]及相关拟合参数(表 2). 从表 2、 图 3可以看出海沟河小流域SOM半方差函数拟合模型符合高斯(Gaussian)模型，决定系数R²为0.751，能够有效反映研究区SOM的空间结构特征^{[25, 26]}.

图 3

Fig. 3

图 3 SOM各向半方差函数 Fig. 3 Sei-variogram charts for SOM

表 2(Table 2)

表 2 SOM半方差函数模型及参数 Table 2 Parameters of the semi-variogram models for SOM 理论模型 块金值(C0) 基台值(C0+C) 块金效应(C/C0+C) 变程/m 决定系数(R2) 残差(RSS) Gaussian 0.036 0.253 0.875 1 800 0.751 5.723E-03

表 2 SOM半方差函数模型及参数 Table 2 Parameters of the semi-variogram models for SOM

一般认为块金值(C₀)能够揭示区域变量在小于抽样尺度上发生突然变异的程度. 研究区SOM的块金值较小，为0.036，表明在小于抽样尺度上几乎不存在某种不可忽略的生态过程对有机质的空间分布产生影响. 基台值(C₀+C)通常用来表征系统内总的变异程度，块金效应(C/C₀+C)表征系统内非随机部分引起的变异占系统总变异的比例. 研究区块金效应值为0.875，该值明显高于黄土高原等区域的相关研究值^{[8, 12]}，表明由于受人类开发利用时间较短，海沟河小流域SOM空间变异受人类活动主要是对资源的开发利用、 耕作措施(施肥、 秸秆处理等)等随机因素影响很小，受地形、 成土母质等结构性要素影响较大，与黄土区相关研究有较大的差异性. 变程是衡量土壤属性特征空间自相关性范围的定量指标，海沟河小流域SOM的变程为1 800 m，说明半方差函数在间隔距离为1 800 m时达到平稳阶段，SOM在该尺度内存在较强的连续性.

对研究区SOM进行更加深入的各向异性分析，绘制在0°(北-南)、 45°(东北-西南)、 90°(东-西)和135°(东南-西北)方向上的半方差函数图[图 3(b)]，结果表明：研究区SOM在不同方向上的变异程度呈现较大差异性，其中，0°方向变异值以0.2为中心做上下变化，变异程度小，当间隔距离持续增加时，变化幅度较大，这与研究区中西部以平原为主，东部小范围内为山地为主有关，南北方向则呈同一梯度趋势； 在东西45°方向上的变异程度随着间隔距离的增加变化幅度较为平缓，变异值在0.3左右变动，这主要是因为研究区东北-西南方向以平原耕作为主，地势及地表景观变化小； 在90°方向上半方差函数变化最为剧烈，随着间隔距离的增加而变大，在约4 000 m处有机质变异程度达到最大值，之后出现交替减小趋势，与研究区东西方向地形的剧烈变化相关； 在135°方向上有机质的半方差函数变化较大，且分别在约2 200 m和6 400 m处出现不同拐点，整体上随间隔距离的增加而变大，在间隔距离接近8 000 m处达到最大值. 总言之，研究区各向半方差函数模型表现出不同的特征，各向异性显著. 这对指导该区域采样具有十分重要的意义，在后续研究布设采样点时，在0°(南-北)方向上，特别是同一海拔梯度和土地利用方式下，采样点可以稀疏些，其他方向要适当地密集布设采样点.

在对研究区SOM进行半方差函数分析的基础上，选择高斯(Gaussian)模型，对海沟河小流域SOM进行普通克里格(Ordinary Kriging)插值，根据高程与SOM相关性分析结果(表 4)，SOM分布与高程相关系数为0.427，为显著水平，因此选取高程作为协变量进行协同克里格(Co-Kriging)插值(图 4)，两种插值均通过ArcGIS 10.0自身的交叉验证进行插值效果的检验. 如表 3所示，综合各项验证指标参数及参考文献^[25]，两种插值结果的均方根预测误差均比较小，且均方根预测误差值与平均标准误差接近，插值效果良好. 普通克里格和协同克里格插值结果的标准平均值误差分别为0.021 3和0.009 8，均接近于0，效果很好，标准均方根误差分别为1.055 7和1.034 9，均较接近于1，表明对海沟河小流域SOM进行普通克里格和以高程为协变量的协同克里格插值效果均较好，选择显著相关的高程作为协变量进行协同克里格插值能够提高插值精度.

表 3(Table 3)

表 3 克里格空间插值交叉验证结果 Table 3 Cross validation results of Kriging spatial interpolation 插值类型 均方根误差 平均标准误差 标准平均值误差 标准均方根误差 普通克里格 11.16 10.55 0.021 3 1.055 7 协同克里格 10.74 10.36 0.009 8 1.034 9

表 3 克里格空间插值交叉验证结果 Table 3 Cross validation results of Kriging spatial interpolation

表 4(Table 4)

表 4 研究区SOM含量与外部因子相关性分析结果 1) Table 4 Results of correlation analysis between SOM and external indicators in the study area 项目 高程 坡度 平面曲率 剖面曲率 距主干道距离 距居民点距离 距水系距离 转换前 0.427** -0.091 0.128 0.158 -0.087 -0.113 -0.346** 转换后 0.303** -0.130 0.028 0.111 -0.107 0.028 -0.220* 1)**表示在 0.01 水平(双侧)上显著相关； *表示在 0.05 水平(双侧)上显著相关

表 4 研究区SOM含量与外部因子相关性分析结果 1) Table 4 Results of correlation analysis between SOM and external indicators in the study area

图 4

Fig. 4

图 4 研究区SOM空间分布 Fig. 4 Spatial distribution of SOM in the study area

由图 4可以看出：两种插值结果在8.99-31.53 g ·kg^-1之间，且均呈现出东高西低呈阶梯状分布趋势，空间特征明显，与研究区地形变化具有较强的一致性. 南北方向SOM含量变化微弱，与该流域南北方向地形变化平缓关系密切. SOM的高值区均出现在东部海拔较高的山地林区，随海拔降低而逐渐减少，同时，可以看到SOM在研究区中西部地区含量变化小，空间分布较均一，这可能与研究区中西部多平原耕作区、 土地利用方式相对简单有很大的关系，说明海沟河小流域SOM受流域内地形因素以及由其产生的地表环境差异的影响较大. 综合比较普通克里格和以高程为协变量的协同克里格插值效果发现，两种插值效果整体趋势相对一致，在东部海拔较高的山区差异较为显著，且协同克里格插值结果较普通克里格法在局部区域表现更加精细，有机质含量范围变大，能够准确反映有机质空间分布的局部特征.

一般认为，SOM受区域成土母质、 地形、 气候等自然要素的主导作用较大^{[1, 2, 8, 9, 10]}，而区域环境因素及其变化被认为在短期内对SOM的空间分布也产生着强烈的影响. 考虑到小流域内土壤成土母质及气候要素的相对一致性，本研究在对海沟河小流域SOM进行空间异质性分析的基础上，结合研究区土地利用现状、 人类活动、 交通干扰、 地貌形态及外部环境等因素对海沟河小流域SOM的空间变异和分布格局进行主控因子分析.

根据SOM的形成机制过程，地球表面植被分布对其形成与累积有着强烈的影响，以海沟河流域1992-2014年间土地利用变化情况(图 5)为基础可以看出，研究区主要地类未发生大的变化，东部山区长期以来分布着大面积的森林植被，在较长的时间里该区域的人类活动较少，各类植被的凋零物、 死亡植物残体及根系经过长期稳定的演变过程，最终促进了该区域SOM的大量累积过程，这也是东部山地林区有机质富集的主要原因. 而在西部平原耕作区，自1992年以来水田与旱地之间发生了明显变化，旱改水现象较多，这种耕作方式的变化以及由此产生的育苗、 施肥、 秸秆处理等频繁的人类活动，一定程度上减缓甚至阻止了有机质的形成过程，加速了其分解及释放过程，导致了SOM含量下降，可见长久稳定的土地利用方式也有助于SOM及相关养分的累积^{[27, 28]}. 从1992-2014年土地利用变化中可以看出在东西方向上土地利用方式存在长期的复杂交替现象(图 5). 促进形成了原始的东西方向强的变异性，近些年来，虽然土地利用向着规模化方向发展，土地利用呈现集中连片现象，但东部林地、 中部旱地、 西部水田的格局未发生变化，SOM在东西方向的强变异性也未发生变化. 综合以上分析，除受研究区地形等结构性要素影响外，从土地资源的开发伊始，土地利用的方式及强度也在很大程度上影响了研究区SOM的空间变异及其分布.

图 5

Fig. 5

图 5 研究区1992-2014年土地利用变化 Fig. 5 Land use change from 1992 to 2014 in the study area

其次，对SOM空间分布(图 4)与研究区地形(图 5)进行对比，可以发现研究区东高西低的地形分布与SOM的分布呈现强烈的相关性，随着海拔高度的增加土壤中有机质含量也逐步增加. 本研究参考已有的土壤养分分布及影响因素的相关分析^{[8, 22, 29, 30]}，结合研究区地形地貌、 城镇居民点、 交通及水系分布现状，对SOM含量及对数转换后的值 与相关地形环境因素进行相关性分析(表 4)，并进一步对影响区域内SOM空间变异及分布的主要因素进行分析.

从表 4可以看出，SOM含量与高程呈显著性正相关水平(r=0.427)，与水系距离呈显著性负相关水平(r=-0.346)，是影响研究区SOM含量空间分布的主要因子. 这一方面是因为受地形等结构性要素的长期影响形成了地表较稳定或较易形成地带性森林土壤的气候条件，有利于SOM的积累； 另一方面在海拔较高的山地人类活动干扰小，对土壤中自然状态下有机质的破坏十分有限，因此保留了较高的有机质含量. 总之，地形要素通过直接地或者间接地形式在很大程度上主导了海沟河小流域SOM的空间分布，高程是其显著相关的外部因素，但不能直接决定其形成过程. 水系的分布在土壤养分的长久运移过程中影响了土壤养分的运移方向及累积过程，距离水系近的地方地势一般较低，一定范围内是土壤水分及养分的“汇”区域，极大地促进相关累积过程，同时，河流水系附近土壤厌氧呼吸也是SOM含量较高的一个原因，水系分布是大时间尺度研究土壤养分运移规律不可忽视的指标因素.

最后，通过实地采样并记录采样点周围环境状况，最终与样点SOM含量进行比较分析，发现采样点周边环境如垃圾、 动物粪便堆积、 畜牧养殖等情形会造成采样点附近SOM含量明显变大等异常现象，增加了SOM的空间变异特征，同时也增加了SOM空间插值的不确定性，在未来土壤采样中要做适当的规避.

科学准确地掌握黑土区SOM的空间变异和现状格局及其影响因素是进行黑土区SOM保护的基础工作. 海沟河小流域SOM含量平均水平为19.35 g ·kg^-1，与黑土区其他区域相比，处于较低水平^{[23, 24]}，土壤肥力保护和SOM培育有待进一步加强； 海沟河小流域SOM含量存在空间自相关性，并且不同方向的自相关性表现出显著性差异，其中，90°(东-西)方向上SOM空间自相关性能够反映整个研究区SOM的自相关性特征，这与研究区东西方向上地形地势的高低变化具有密切关系^{[21, 30]}. 在空间异质性方面，SOM空间变异为中等变异程度，且不同方向上的变异程度呈现较大差异性，这一点与黄土高原等区域空间变异性研究相比，其空间变异受人类活动(对资源的开发与利用)、 耕作措施(施肥、 秸秆处理等)等随机因素影响更小，受地形、 成土母质等结构性要素的主导作用更加明显^{[5, 8, 10]}.

以半方差函数分析为基础，可以指导海沟河小流域SOM普通克里格和协同克里格插值及交叉验证，以高程为协变量的协同克里格插值较普通克里格插值能够提高插值的精度，且达到很好的效果，与其他研究结果相近^{[17, 20]}. 在考虑小流域内土壤成土母质及气候要素相对一致的前提下，SOM的空间分布与区域地形以及受地形要素影响较大的水系分布呈极显著的相关性，与交通、 居民点分布等相关性微弱，其原因可能是地形要素通过改变地表植被、 土地利用以及小区域内的湿度温度等环境要素而直接地或者间接地影响了研究区SOM的形成和累积过程. 通过实地采样记录及对应样点SOM含量的比较分析，发现采样点周边环境如垃圾、 动物粪便堆积、 畜牧养殖等情形增加了SOM的空间变异的特征，同时也增加了SOM空间插值的不确定性，是相关研究中的重点，在未来的精细研究中要做适当的规避.

以小流域为研究单元，综合考虑影响SOM空间分布的自然、 社会要素对流域内SOM空间变异、 分布格局及影响因子进行研究具有重要意义. 未来在土壤养分空间相关性及变异特征分析的基础上，如何实现主控因子的进一步量化，确定其与土壤属性间的定量关系，来提高SOM数字制图的精度，还需要进一步研究探讨.

(1)海沟河小流域SOM含量在6.29-39.94 g ·kg^-1之间，均值为19.35 g ·kg^-1，处于较低水平，空间变异程度为中等变异，存在空间自相关性，且在不同方向表现出显著性差异. 半方差函数分析显示SOM的变程为1 800 m，空间变异整体上受人类活动(耕作制度、 施肥措施、 秸秆处理)等随机因素影响小，受地形、 成土母质等结构性要素的影响大，90°(东-西)方向变异程度在空间变异中均占主导作用.

(2)协同克里格插值法可以提高SOM的空间插值精度，插值结果表明，海沟河小流域SOM的空间分布整体上呈现东高西低、 逐步减小的特征，SOM空间分布受地形海拔等自然因素以及受其制约的土地利用的影响较大，受交通道路、 居民点等要素的影响较小，采样点周边环境增加了SOM空间分布的不确定性.

(3)以流域为单元，综合考虑影响SOM空间分布的自然、 社会要素对我国黑土区海沟河小流域SOM的空间变异、 分布格局及影响因子进行研究具有重要意义和科学价值.