砒砂岩区集中分布于黄河流域晋陕蒙接壤的中心地带，砒砂岩是一种成岩程度较低的岩石互层，矿物颗粒胶结强度差，在风水复合侵蚀作用下极易风化溃散^[1]。砒砂岩区土壤贫瘠，气候干旱，天然植被稀疏，侵蚀模数高达3×10⁴ ~ 4×10⁴ t/(km²·a)，是黄河流域中游段主要的粗砂来源区^[2]。开展该区域土壤有机质空间自相关性的研究，对于科学指导植被建设、改善区域生态环境质量、减少入黄泥沙有重要现实意义。从现有的研究来看，国内外关于土壤养分空间异质性的研究成果较为丰富^[3-5]，但针对土壤养分空间自相关性的研究相对较少。2000年前后主要包括对土壤微量元素、土壤重金属和土壤有机碳的空间自相关性的研究^[6-8]，以及近年高凤杰等^[9-10]对黑土区土壤含水量和pH空间自相关性的研究，而砒砂岩区地形和土地利用方式等因素对表土有机质空间自相关性的影响机理尚不明确。基于以上背景，本研究选取砒砂岩区典型小流域为研究对象，应用地统计和灰色关联度的理论方法，研究砒砂岩区小流域表土有机质空间自相关格局及影响因素，以期为砒砂岩区土壤养分的空间预测提供数据支撑，并为该区生态修复工作提供科学依据。

空间自相关性是指在特定空间范围内的观测变量之间存在的潜在相互依赖性。任何变量间都存在相互联系，且这种联系随着距离的缩短而愈加紧密，这个理论便是具有划时代意义的地理学第一定律^[11]。土壤是地表上结构和功能最为复杂的生态系统之一，同时具备生态要素和地理要素的特征，地理学第一定律同样适用于土壤。地统计学与土壤学的结合为土壤有机质的空间预测和不确定性分析提供了有效工具^[12-13]，然而大多数专家学者的研究多集中在土壤空间异质性(即现象描述)的研究上^[14-15]，却忽视了产生空间异质现象的本质正是土壤的空间自相关性。土壤有机质作为土壤系统中功能和性质最活跃的组成物质之一，对土壤中养分循环和生态功能的稳定发挥有重要影响。表层土壤有机质更容易受到气候、地形和人类活动等外环境的干扰，从而表现出局部空间的聚集或离散特征^[16-17]。研究表层土有机质空间自相关性及其对区域环境的响应机制，是土地资源可持续利用和区域生态环境修复等重大策略决策的基础。

1 材料与方法 1.1 研究区概况

研究区鲍家沟流域行政上隶属于内蒙古准格尔旗暖水乡，位于110′31′  110°35′E，39°46′  39°48′N，海拔1 110  1 300 m，流域面积12.64 km²，地形北高南低；属于温带大陆性气候，冬季漫长干燥，夏季短暂温热，年均气温7.2 ℃，年降水量400 mm；流域面积较小，土壤类型为栗钙土。研究区自20世纪80年代末陆续开展生态移民搬迁，完全隔离人为干扰，封育禁牧并大面积营造人工植被，土地利用类型以林地、天然草地和裸地为主。主要人工植被有：油松(Pinus tabuliformis)、侧柏(Platycladus orientalis)、柠条(Caragana korshinskii)以及沙棘(Hippophae rhamnoides)等；草本植物主要有羊草(Leymus chinensis)、猪毛菜(Salsola nitraria)和阿尔泰狗娃花(Heteropappus altaicus)等。

1.2 样品采集及测定

以研究区1∶10 000地形图和2017年土地利用现状图(1∶10 000)为基础，根据研究流域地形特征和各类用地面积，布设取样点150个，样点布设情况如图 1所示。野外土壤样品采集于2018年7月，用手持GPS采集样点坐标，并详细记录周围地形条件、基岩裸露情况、植被盖度和土层厚度信息。研究流域基岩出露面积较大，部分地区土层厚度小于10 cm，下伏砒砂岩。用环刀和铝盒采集0  10 cm土样，每一取样点随机采集3份土样，土壤样品带回室内风干、处理以供分析测定。土壤有机质的测定采用K₂CrO₇外加热法，土壤容重和含水量的测定采用烘干法。

1.3 环境因子的选取

本研究小流域面积较小，流域内气候和土壤类型基本一致，且流域内无村民居住且无生产开发建设项目，人类活动较少，因此本研究不考虑气候、土壤类型和人类活动因素。参考前人的研究成果^[18-19]，地形和土地利用因素是影响干旱半干旱区小流域土壤养分变化的主要因素，因此本研究将从地形和土地利用2个方面选取环境因子。砒砂岩区小流域具有丘陵沟壑地貌特征，局部区域地形的起伏和褶皱会引起水热资源的再分配，可以通过坡面所处的空间位置和方位角即坡位、坡向来反映。因此地形方面选取了海拔、坡位和坡向3个因子，坡向按照方位角可以划分为阴坡、阳坡和平地，坡位按照地形地貌特征可以分为沟谷、坡面和山脊^[20]。同时，考虑到研究区地形破碎且侵蚀剧烈的特性，加入能综合反映坡度和径流侵蚀特征的TWI(地形湿度指数)共5个因子。TWI的计算公式如下^[21]：

$ {\rm{TWI}} = \ln \frac{{Ac}}{{{{\tan }_\alpha }}} $

(1)

式中：Ac为垂直与特定水流方向的汇流面积，α为坡度。

1.4 分析方法

空间自相关性常用全局和局部两类指标度量，全局指标可以反映研究区整体的空间模式，局部指标则体现了每一个空间单元与邻近单元就某一属性的相关程度。本研究应用Moran指数的全局性I_G和局部性I_L对表土有机质的全局空间自相关性和局部空间自相关性进行描述。

全局性I_G、局部性I_L指数和Z值得分的计算过程如下：

$ {I_G} = \frac{{n\sum\limits_{i = 1}^n {\sum\limits_{j = 1}^n {{w_{ij}}({x_i} - \bar x)({x_j} - \bar x)} } }}{{(\sum\limits_{i = 1}^n {\sum\limits_{j \ne i}^n {{w_{ij}}} } )\sum\limits_{i = 1}^n {{{({x_i} - \bar x)}^2}} }} $

(2)

$ {I_L} = \frac{{n({x_i} - x)\sum\limits_{j \ne i}^n {{w_{ij}}({x_j} - x)} }}{{\sum\limits_{i = 1}^n {{{({x_i} - x)}^2}} }} $

(3)

$ E(I) = \frac{{ - 1}}{{n - 1}} $

(4)

$ V(I) = E({I^2}) - E{(I)^2} $

(5)

$ Z = \frac{{1 - E(I)}}{{\sqrt {V(I)} }} $

(6)

式中：n为空间数据的个数，x_i和x_j分别为i区和j区的要素属性值，$\bar x$为所有数据的平均值，w_ij为空间权重矩阵元素。Moran’s I > 0表示空间正相关性，Moran’s I < 0表示空间负相关性，其值越大，空间相关性越明显，相反，相关性越差；而当Moran’s I = 0，空间呈随机性。

采用地统计学和灰色关联分析法相结合，运用GeoDa1.12和GS+9.0进行半方差函数和空间自相关性分析，利用ArcGIS10.2进行克里格插值成图，并进行空间交叉验证，应用DPS15.10软件子进行灰色关联度分析，数据可视化在Orgin9.0下完成。

2 结果与分析 2.1 表土有机质空间分布特征

土壤有机质在空间上的变异规律可以用半方差函数的理论模型来描述。表 1为研究区表土有机质含量的基本描述性统计特征，研究区土壤有机质含量的平均值为9.56 g/kg，K-S检验的P值为0.442 > 0.05，数据符合正态分布，可以进行下一步分析。图 2A为半方差函数的拟合图，研究小流域土壤有机质空间变异的半方差函数拟合模型为指数模型，拟合的决定系数R²为0.792，能够较为准确地反映研究区土壤有机质的空间变异特征。块金值可以反映空间距离小于抽样尺度时区域变量的突异程度，基台值则表征了总体的变异程度，研究流域不同土层深度半方差函数的块金值为5.85，基台值为14.87；块金效应可以说明变量的空间相关程度，研究区块金效应值为60.66%，属中等程度的空间相关性。变程衡量了相似斑块的空间连续范围即土壤有机质空间相关性范围，土壤有机质空间变异的变程为1 579 m。综合来看，研究区土壤有机质的空间变异是结构因素和随机因素综合作用的结果。

克里格插值生成表土有机质的空间分布图(图 2B)，可直观地反映出土壤有机质的空间分布特征。交叉验证结果显示，标准平均值误差为0.028，标准均方根误差为1.046，插值效果良好。研究流域地形西高东低，有机质的高值区集中在地势平缓的流域中部，有机质的低值区集中在分水岭地带和流域出口区域，总体呈沟道 > 坡面 > 山脊的规律。表层土壤有机质含量空间差异明显，呈明显的斑块化和岛状聚集趋势，这可能与研究区的土地利用方式和基岩裸露程度有关。

2.2 表土有机质整体空间格局分析

地统计学中常用Moran指数来判断变量是否在空间中存在聚集特征。研究区表土有机质全局空间自相关性分析结果显示，全局Moran指数Z值得分为2.80，大于临界值1.96(P < 0.01)，具有显著的空间自相关性。进一步对其空间自相关性与空间距离的关系进行分析，各向同性分析结果(图 3A)显示，随着滞后距离的增加，空间自相关性由正相关转向负相关；当滞后距离为0.2 km时，表土有机质Moran指数最高为0.606，其空间相关性最强；当滞后距离处于1.1 ~ 1.9 km范围时，表土有机质的空间自相关性最弱，Moran指数接近0，有机质在空间上呈随机分布；滞后距离大于2.2 km时，Moran指数均转为负值，土壤有机质含量由相对聚集转向离散分布。

各向异性分析结果(图 3B)显示，在90°和135° (西北、东南和正东正西)方向上，在0.15 km处表土有机质具有强烈的正向空间自相关性，在1.3 km处的空间自相关性最弱，而后随着滞后距离增加，空间自相关性由正相关转向负相关，结合流域的地形特征来看，主沟的径流流向为由西北至东南，此方向表土有机质空间自相关性与主沟径流流向一致。在45°即东北、西南方向上，表土有机质整体呈正向自相关关系，当距离达到3 km时，空间自相关性达到最大，Moran指数为0.689，而支沟径流流向为由东北至西南以及西南至东北方向，此方向表土有机质空间自相关性与支沟径流流向一致。而在0°即正南正北方向上，Moran指数围绕坐标轴周期性波动，表土有机质的空间自相关性没有表现出明显的变化趋势。综上，全局空间自相关性的分析结果，可以确定砒砂岩区表土有机质在空间上存在聚集趋势，且这种聚集趋势与流域的径流流动方向有关，但并不能说明产生聚集的具体位置和聚集机制，需要通过进一步的局部自相关性的分析。

2.3 表土有机质局部空间格局分析

表土有机质的局部空间自相关特征可以用Moran散点图来表示。图 4A为以贝叶斯标准化土壤有机质为X轴，空间滞后向量为Y轴，制作的Moran散点图，该图的4个象限分别代表了高值聚集(HH)、低高值异常(LH)、低值聚集(LL)和高低值异常(HL)4种局部关联形式，散点图的斜率即为局部空间自相关的Moran指数(0.186 1)，可以看出图中大部分散点都分布于HH区和LL区，表明相邻点位的土壤有机质差异较小，空间自相关性较强，也表明土壤有机质呈明显的高值聚集和低值聚集趋势。相反，分布在二四象限的散点相对较少，说明出现高值和低值聚集的异常现象的概率较低。

表土有机质的局部聚集特征可以通过空间关联局部指标的表示，计算Moran指数的Z值，将在显著性水平(α=0.05)的条件下通过检验的为正的Z值区域标注在图上。图 4B直观地反映了土壤有机质高低值的聚集区域和异常区域的具体位置，研究流域大部分区域土壤有机质的局部Moran指数未达到显著性水平(α=0.05)，空间分布呈随机分布状态。而达到显著性水平的区域其土壤有机质含量与邻近区域的关联性较强，以高值聚集和低值聚集类型为主。此外，还存在少量的高低值聚集异常区。

2.4 表土有机质空间自相关性与环境因子的灰色关联分析

基于表土有机质局部空间自相关特征(图 4)，对聚集区域和异常区域相应的地形因子和土地利用因子进行统计，见表 2。由表 2可以看出：聚集区和异常区主要分布在海拔1 146.56 ~ 1 289.85 m范围区域，其余海拔高度区域为差异不显著区域。其中，高值聚集区其海拔在1 146.56 ~ 1 220.64 m范围，地形湿度指数为3.96 ~ 13.84，坡向特征包括阳坡和平面2类，坡位特征包含沟谷和坡面2类，土地利用方式包含林地类型1类；低值聚集区其海拔在1 190.00 ~ 1 260.00 m范围，地形湿度指数为4.43 ~ 6.49，坡向特征包括阴坡和平面2类，坡位特征包含沟谷、坡面和山脊3类，土地利用方式包含林地、草地和裸地类型3类；高低值异常区，其海拔在1 162.50 ~ 1 289.85 m范围，坡向在8.85 ~ 247.30°，地形湿度指数为3.50 ~ 14.89，坡向特征包括阳坡、阴坡和平面3类，坡位特征包含沟谷、坡面和山脊3类，土地利用方式包含林地和草地类型2类。综合来看，高值聚集区主要出现在流域中部水分条件较好的沟谷且分布着乔木林地的地带；低值聚集区主要集中在流域的下部以及出口处水分条件相对较差的坡面的裸露区域；而高低值聚集的异常区在流域的各坡位特征上均有分布，主要集中在坡面的草地类型区。说明高低值聚集异常区虽然在空间上分布较为随机，但其产生与局部区域环境的转变有关，如土壤侵蚀引起的地形变化和养分流失。

土壤有机质空间自相关性与环境因子具有各自的变化态势，对二者间的关系进行定量描述可以应用灰色关联分析法来进行分析。表 3中给出了表土有机质聚集和异常区的局部Moran指数(I_L)与地形因子和土地利用因子的灰色关联系数。高值聚集区地形和土地利用因子与土壤有机质空间自相关性的关联系数均为正值，按关联程度大小依次为地形湿度指数(0.498) > 海拔(0.323) > 坡向(0.267) > 土地利用方式(0.202) > 坡位(0.137)；低值聚集区地形和土地利用因子与局部Moran指数呈正相关关系，关联程度呈坡位(0.489) > 土地利用方式(0.463) > 海拔(0.417) > 地形湿度指数(0.327) > 坡向(0.205)的排列关系；而高低值异常区与I_L关联程度最大的因子是土地利用方式(0.329)，其余依次为坡位、地形湿度指数、海拔和坡向。综上可知，地形湿度指数是土壤有机质高值聚集区形成的主导因素，坡位因子是土壤有机质低值聚集区形成的主要因素，而土地利用方式则是干扰土壤有机质高/低值空间聚集的关键因素。

3 讨论

相关研究表明，土壤有机质空间异质性可以通过半变异函数模型进行描述，本研究得出，砒砂岩区小流域表土有机质含量的半变异函数模型为指数模型，这与薛志婧等^[22]的研究结果一致。在此基础上，本研究加强了土壤有机质空间自相关性以及局部空间自相关性的研究，对局部变异的不稳定性进行了分析^[23]，并结合流域径流和土壤侵蚀现象对局部变异事件进行了解释。根据半变异函数和空间插值可以对土壤空间结构和变异规律进行描述和预测，但是无法准确反映出发生变异的具体位置和变异机制，亦不能判别某一环境因素变异过程中的作用效果^[24]。与前人研究有所不同的是，本研究在确定表土有机质半变异函数最优模型为指数模型的基础上，对土壤有机质局部空间自相关性进行了研究，得出砒砂岩区表土有机质空间上呈高值聚集、低值聚集和高/低值异常的状态，并通过计算空间关联局部指标定位了聚集和异常发生的具体区域，最终利用灰色关联度分析法分析了环境因素对表土有机质局部自相关性的影响，实现了对砒砂岩区表土有机质空间变异现象描述的同时，阐述其空间变异的形成机制。

土壤有机质的积累和分解过程对环境条件极为敏感，其空间异质性是土壤在自然界中长期演变的结果。相关研究表明，环境因子对土壤有机质空间变异的影响具有明显的尺度效应，大尺度的研究认为气候和地形因素是决定土壤有机质空间分布特征的主导因素，中小尺度的研究则将影响因素归结为地形和人类活动因素两方面^[25-26]。本研究以砒砂岩区小流域为研究对象，表土有机质半方差函数模型的块基比为60.66%，属中等程度的空间相关，说明表土有机质空间变异来源于结构性因素和随机性因素。参考了前人的研究，在小流域尺度且土壤类型一致的前提下，结构性因素主要来源于地形差异，随机性因素则来源于土地利用方式的演变和土壤侵蚀的干扰^{[18, 27]}。本研究中表土有机质空间自相关性和灰色关联分析的结果表明，土壤有机质整体的空间自相关性与流域径流流向有关，局部性的高值和低值聚集区的形成与地形湿度指数和坡位有关，而土地利用方式则会影响聚集区的形成促成异常区的产生。这与砒砂岩区的地理特征有关，该区沟壑密布，侵蚀剧烈，表土细颗粒物质由山脊和坡面随着径流的冲刷搬运沉积在沟底地带，因此表土有机质空间自相关性与径流流向的关系极为密切，这也是高值聚集区多形成在沟谷平坦地带，而低值区多形成于山脊处和流速较大沟口地带的原因。自20世纪80年代末起，砒砂岩区大面积开展人工植被建设，土地利用类型以林地为主，林地大量蓄积的枯枝落叶为有机质的形成提供了丰富的来源^[28]，所以林地土壤有机质含量高于草地和裸地。因此，在区域地形要素一致时，土地利用方式的改变会造成相邻单元土壤有机质的空间差异，从而干扰空间自相关性的稳定状态，在生产实践中，掌握土壤有机质整体的空间变异规律以及局部聚集和异常区形成机制，可以通过适当的人为干扰，改善造林区域的土壤水肥条件，促使生态修复工作向理想化方向发展。

4 结论

1) 砒砂岩区小流域表土有机质半方差函数为指数模型，有机质空间呈明显的斑块化和岛状聚集趋势，总体呈沟道 > 坡面 > 山脊的规律。

2) 砒砂岩区小流域表土有机质全局空间自相关性与径流流向关系密切，局部空间自相关特征的高值聚集区分布沟谷和林地，低值聚集区分布在水分条件相对较差的裸露坡面，高低值异常区分布在水分和侵蚀特征多变的坡面草地。

3) 地形因素是促使砒砂岩区小流域表土有机质高值聚集区形成的主导因素，坡位对于表土有机质低值聚集区形成起决定作用，土地利用能够影响表土有机质高/低值聚集区的形成。

4) 基于砒砂岩区表土有机质的空间变异规律和局部空间自相关特征，建议在生产实践中，对于土壤养分基础较差的山脊和陡坡地带，应以封育禁牧保护天然草地为主；对于水肥条件相对较好的缓坡坡面，可以营造合理密度的灌木林草带，并配合水平沟和鱼鳞坑等坡面径流调控工程；而对于水肥条件适宜的沟谷地带和乔木林带，可以发展适当的农业生产和林下经济生产。