裸露砒砂岩区小流域土壤侵蚀空间自相关特征及影响因素

杨振奇, 郭建英, 秦富仓, 刘铁军, 张铁钢

(1.水利部 牧区水利科学研究所, 内蒙古 呼和浩特 010020; 2.内蒙古阴山北麓荒漠草原生态水文国家野外科学观测研究站, 内蒙古 呼和浩特 010020; 3.内蒙古农业大学 沙漠治理学院, 内蒙古 呼和浩特 010019)

空间自相关性是指在特定空间范围内的观测变量之间存在的潜在相互依赖性。任何变量间都存在相互联系，且这种联系随着距离的缩短而愈加紧密，这个理论便是地理学第一定律[1]。土壤是结构和功能最为复杂的生态系统之一，同时具备生态要素和地理要素的特征，地理学第一定律同样适用于土壤学。对于一个封闭的小流域而言，裸地、耕地及沟道等充当土壤侵蚀的“源”，林地和草地等则是土壤侵蚀的“汇”，但“源”和“汇”对象的空间分布存在一定的不确定性，从而导致流域各点侵蚀的发生机理和调控机制的差异。地统计学与景观生态学的结合为土壤侵蚀的空间预测和不确定性分析提供了有效工具[2-3]，干旱和半干旱区的土壤侵蚀过程更容易受到地形特征和植被景观格局的干扰，从而表现出局部空间的聚集或离散特征[4-6]，然而目前应用空间自相关理论和方法描述土壤侵蚀空间特征，并基于此特征分析区域环境因素对其影响机制的研究鲜见报道。

黄河流域的生态保护和高质量发展，是我国新时代生态文明建设的重要内容。砒砂岩区是黄河流域上中游地带的粗泥沙集中来源区，该区基岩成岩程度低，裸露基岩极易风化溃散，引发严重水土流失作用[7]，砒砂岩区按照基岩出露程度和覆土类型被分为盖土区、盖沙区和裸露区，以裸露砒砂岩区的基岩出露面积最大，土壤侵蚀作用也最为剧烈[8]。而裸露砒砂岩区小流域土壤侵蚀在空间上是否具有自相关特征？小流域土壤侵蚀空间自相关性受何种因素影响等问题尚不明确。为此，本研究选取裸露砒砂岩区鲍家沟小流域为研究对象，应用地统计和灰色关联的理论方法，研究小流域土壤侵蚀空间自相关格局及影响因素，以期为砒砂岩区的土壤侵蚀的预报提供数据支撑以及该区生态修复工作提供科学依据。

1 材料与方法

1.1 研究区概况

研究区鲍家沟流域位于鄂尔多斯市准格尔旗，地理坐标为110°31′—110°35′E，39°46′—39°48′N，海拔在1 110～1 300 m之间，流域面积12.67 km2；流域地形北高南低；气候属典型大陆性季风气候，冬季漫长干燥，夏季短暂温热，年均气温7.2 ℃，年降雨量388.3 mm；研究区自20世纪末实施生态移民政策，流域内目前无人居住，土地利用类型以林地和草地为主。植被以人工植被为主，常见青杨(Populuscathayana)、旱柳(Salixmatsudana)、油松(Pinustabuliformis)、侧柏(Platycladusorientalis)、柠条(Caraganakorshinskii)、沙棘(Hippophaerhamnoides)、山杏(Prunusarmeniaca)等；草本植物主要有羊草(Leymuschinensis)、猪毛菜(Salsolanitraria)、阿尔泰狗娃花(Heteropappusaltaicus)等。

1.2 数据来源与预处理

本研究基于研究区1∶1万地形图矢量化后生成DEM数据，采取航拍照片与实地调查相结合的方式解译研究区土地利用现状，利用Landsat遥感影像计算归一化植被指数(NDVI)和植被覆盖度数据。研究流域内设典型样地150个，用于植被调查和土壤样品采集，此外还布设径流小区6个，分别为裸露基岩、油松、山杏、沙棘、柠条和草地，坡度均为20°，小区规格均为2 m×5 m，观测内容为次降雨的降雨过程及产流产沙量，降雨过程采用径流小区配备气象站(HOBO-U30)进行观测，雨后收集径流收集桶中的径流泥沙测量体积，烘干后进行称重。

由于研究面积相对较小且无泥沙资料，因此采用参数相对简明的修正后的通用流失方程(RUSLE)对流域的侵蚀现状进行估算，其中R因子计算方法，基于研究区2014—2019年降雨的自记观测资料，通过降雨的I30可以计算出流域的年平均降雨侵蚀力为1 033.98 MJ·mm/(hm2·h·a)；土壤可蚀性K因子基于采集土壤样品团粒体的几何平均直径，采取网格法与随机法相结合的方法，于2018年7—8月，在流域内布设典型样地150个，采集0—20 cm土壤样品，每点采集3份，供450份样品，应用Le Bissonnais法测定土壤团粒体含量，采用Shiraz模型[9]计算，经克里格插值生成：LS因子基于DEM数据，根据黄土丘陵区地形特征[10]，采用累积流量法[11]提取；采用史培军等[12]在砒砂岩区地区的研究成果计算C；P因子按土地利用类型进行赋值[13]；计算结果依据《土壤侵蚀分类分级标准》[14]进行分级，图1—2为流域各因子及土壤侵蚀的空间分布特征。

图1 鲍家沟流域RUSLE模型因子的空间分布特征

图2 鲍家沟流域土壤侵蚀强度分级

1.3 数据分析方法

地统计学中常用Moran指数来判断变量是否在空间中存在聚集特征，用p值和Z得分来半段是否拒绝零假设，p值表示概率，Z值为标准差的倍数，综合p值和Z得分值来判断要素是否在统计意义上的显著性聚类或离散模式。Moran散点图表现某个变量的观测值向量与它的空间滞后向量之间的相关关系，通过散点图的形式表现出来。其中横轴对应观测值向量，纵轴对应空间滞后向量，空间滞后向量是空间权重矩阵与因变量向量相乘的结果，即为空间滞后向量即该观测值邻域的加权平均。全局空间相关指数Moran’sI就是空间滞后向量对观测值向量线性回归的斜率系数。空间自相关性常用全局和局部两类指标度量，全局指标可以反映研究区整体的空间模式，局部指标则体现了每一个空间单元与邻近单元就某一属性的相关程度[15]，全局性IG和局部性IL指数的计算过程[16]为：

(1)

(2)

参考前人的研究[17-18]，从面积指标、形状指标和聚集指标3方面选取了斑块面积(class area)、面积加权的平均形状指数(area weight shape index distribution)、面积加权的平均斑块分形维数(area weight fractal index distribution)、平均邻近指数(mean proximity index)4项指标。

各地形因子的提取运用ArcGIS的空间分析功能，提取了包含海拔、坡度、坡位、坡向和曲率，同时考虑到研究区地形破碎且侵蚀剧烈的特性，加入能反映地形和径流侵蚀的地表切割深度、汇流动力指数(stream power index， SPI)和沉积物运移指数(sediment transport index， STI)。SPI和STI的计算公式[19]为：

SPI=ln(Ac×tanβ×100)

(3)

(4)

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

基于计算得出的小流域土壤侵蚀图，应用ArcGIS和GeoDa软件的空间统计功能分别采用全局空间自相关和局部空间自相关，对土壤侵蚀的空间自相关性进行分析，提取聚集区和异常区的地形因素和植被格局因素，运用灰色关联方法分析三者间的关系。

2 结果与分析

2.1 小流域土壤侵蚀的空间格局和空间自相关性

研究流域土壤侵蚀模数全局空间自相关性分析结果显示，全局Moran指数Z值得分为136.87，土壤侵蚀模数在空间上随机聚类的可能性小于1%(p<0.01)，具有显著的空间自相关性。对其空间自相关性与空间距离的关系进行分析(图3)，各向同性分析结果显示，表现为随着滞后距离的增加，空间自相关性由正相关转向负相关，表明土壤侵蚀模数的空间关系由聚集转向离散，这是由于随着空间距离的延伸，区域下垫面条件复杂多变，侵蚀程度也随之变化。当滞后距离为146.17时，土壤侵蚀模数Moran指数最高分别为0.443，其空间相关性最强。当滞后距离处于1 595.18～1 812.25范围时，土壤侵蚀模数的空间自相关性最弱，Moran指数接近0，土壤侵蚀模数在空间上呈随机分布，随后Moran指数均转为负值，土壤侵蚀模数由相对聚集转向离散分布。各向异性分析(图3)结果显示，整体来看，各方向上土壤侵蚀模数的空间自相关性的最高值分布在107.27～161.54处，在741.11～1 600.42处的空间自相关性最弱。从各向异性的角度结合流域的地形特征来看，土壤侵蚀模数的空间聚集特性可以归纳为两类，一类是沿西北至东南(135°)即主沟的径流汇集方向随滞后距离逐渐减弱，另一方向是东西、东北至西南(45°)即各支沟的径流汇集方向上，随着空间滞后距离的增加，各支沟间的侵蚀程度的差异也随之增大，因而空间自相关性由聚集转向离散。而在0°即正南正北方向上，Moran指数围绕坐标轴周期性波动，这可能与流域的地形走势及南北方向不是流域径流的汇流方向有关，导致了土壤侵蚀模数的空间自相关性在南北方向上没有表现出明显的变化趋势。

图3 鲍家沟流域土壤侵蚀全局空间自相关性分析

进一步对裸露砒砂岩区土壤侵蚀模数在的局部自相关性进行分析，图4a为以贝叶斯标准化土壤侵蚀模数为X轴，空间滞后向量为Y轴，生成局部Moran散点图，该图的4个象限分别代表了高值聚集(HH)、低高值异常(LH)、低值聚集(LL)和高低值异常(HL)4种局部关联形式，散点图的斜率即为局部空间自相关的Moran指数(0.425)，可以看出图4a中大部分散点都分布于HH区和LL区，高值聚集区分布相对零散，而低值聚集区分布较为集中。分布在二四象限的散点相对较少，说明出现高值和低值聚集的异常现象的概率较低。土壤侵蚀模数的局部聚集特征可以通过空间关联局部指标的表示，计算Moran指数的Z值，将在显著性水平(a=0.05)的条件下通过检验的为正的Z值区域标注在图上。图4b直观的反映了土壤侵蚀模数高低值聚集区域和异常区域的具体位置，可以看出土壤侵蚀模数高值聚集区主要位于基岩大幅出露的区域，空间上呈小聚集大分散的分布状态，而低值聚集区的土壤侵蚀模数的空间关联性相对较强，呈连片的聚集形式，此外还存在少量的高低值聚集异常区，散落分布在沟道区域。

图4 研究区土壤侵蚀局部空间自相关性

2.2 小流域植被景观格局分析

结合图5和表1可以对研究流域的各斑块单元和斑块类型的景观指数进行完整的解读。斑块面积是景观斑块最基本的属性也是计算其他指标的基础，由研究流域各景观的斑块面积分布图可以看出(图5)，流域中的大面积斑块主要有两类，一类是分布在流域的边界地带的平缓坡面上的林地斑块，另一类则是成片裸露基岩斑块。由表1可知，沙棘林和裸露基岩景观呈集中连片的分布格局，其斑块面积分别为441.15和365.92 hm2，而其他景观则相对破碎，斑块面积较小。

表1 鲍家沟流域斑块类型景观指数

图5 研究区斑块景观格局指数空间特征

面积加权的斑块形状指数反映了斑块形状与正方形的偏离程度，其值越大表明斑块形状越趋于扁长或复杂化，研究流域各类景观斑块形状指数在2.51～5.59之间，按照数值由小到大排列为公路、草地、柠条林、油松林、裸露基岩、沙棘林、河道；面积加权的斑块分形指数则从分形理论区描述斑块的复杂性和规则性，其值越趋近于1，表明斑块形状越规则，各斑块类型的分形维数在1.19～1.29之间，由大到小依次为河道、沙棘林、裸露基岩、油松林、山杏林、柠条林、公路、草地，这也证实了各斑块类型的规则程度。结合图4来看，形状指数和分形维数较小的斑块如草地多分散的镶嵌在大面积的景观中，其多具备面积小而形状规则的特点。形状指数和分形维数较大的斑块多集中在裸露基岩区，由于基岩出露地带，多是形态曲折的沟道，因而其形状指数较高。

邻近指数反映了各类型间的相邻程度，常用来描述景观的破碎程度，其值越小，表明斑块间越离散，相应的景观的破碎程度也就越高。各景观类型邻近指数的平均值在5.25～157.06之间，表现为：沙棘林>裸岩基岩>油松林>山杏林>柠条林>河道>草地>公路的趋势。由图5可以看出，流域中各斑块间平均邻近指数较低的区域主要分布在流域中下部，该区域是狭长流域中相对宽阔的地带，区域地形也相对复杂，平缓的坡面和陡峭的沟道交错分布，平缓的坡面上分布着不同种类的人工林，沟道内植被覆盖较高，相比流域中上部，裸露基岩的面积较小，分布也较为零散。

2.3 土壤侵蚀空间自相关性与植被和地形的关系

基于土壤侵蚀模数的局部空间自相关特征图，对聚集区域和异常区域相应的植被格局因子进行统计。由图6可以看出，有56.44%的高值聚集区位于流域内的基岩裸露区，同时也有44.22%的高低值异常区分别位于该区内，表明流域内基岩裸露地带总体上的侵蚀作用较为剧烈，但也零散分布着水肥条件相对较好，植被覆盖相对较高的局部地区，从而造成了本该是剧烈侵蚀的大面积景观中出现轻度侵蚀的细碎斑块。有96.93%的低值聚集区位于林地和草地类型中，其中30.39%集中在沙棘林地区，同时也有23.96%和31.02%的高值聚集区和高低值异常区位于该区内，表明沙棘林虽然是流域内面积较大的土地利用类型之一，但是并非所有沙棘林都长势良好，部分地区也存在沙棘林的退化和死亡现象，在大面积沙棘林中出现裸露空地，造成局部地区侵蚀高低值聚集异常。从地形上看，有92.81%的低值聚集区位于坡面地带，其中56.37%位于坡面中部，表明坡面的侵蚀作用相对较弱，这可能与坡面的地形条件和植被生长状况有关。有24.06%的高值聚集区位于沟谷地带，同时也有48.51%的高低值异常区位于该地带，表明沟谷的侵蚀作用相对严重，但沟谷地带的地形条件更为复杂多变，可能存在部分植被覆盖相对较高，侵蚀作用相对较弱的地区。

图6 聚集(异常区)与土地利用类型和地形的关系

土壤侵蚀模数、地形因子与景观格局因子都具有各自的变化趋势，且不一定具备简单的线性关系，因此对三者间的关系的定量描述可以应用灰色关联分析法来进行分析。表2中给出了土壤侵蚀模数聚集和异常区的局部Moran指数与地形和景观格局因子的灰色关联程度。可见高值聚集区的Moran指数与地形和景观格局因子间的关联程度由大到小依次为坡度、高程变异系数、汇流动力指数、斑块间最近距离、斑块面积分形维数、地表切割深度、斑块邻近指数、斑块形状指数；低值聚集区的Moran指数与地形和景观格局因子间的关联程度由大到小依次为斑块面积分形维数、斑块间最近距离、地表切割深度、斑块形状指数、高程变异系数、坡度、斑块邻近指数、斑块形状指数；高低异常聚集区的Moran指数与地形和景观格局因子间的关联程度由大到小依次为斑块形状指数、斑块面积分形维数、地表切割深度、高程变异系数、斑块间最近距离、汇流动力指数、坡度、斑块邻近指数。

表2 聚集和异常区局部Moran指数与环境因子灰色关联度

3 讨 论

大部分专家学者对流域尺度的研究表明[20-21]，植被景观的空间分布格局是影响流域土壤侵蚀差异的主要因素，然而单纯分析土地利用和植被覆盖的矢量数据，只能描述其景观格局的几何特征，缺乏与地形因子间的耦合联系，因而也就无法准确的反映土壤侵蚀与植被格局的关系[22-23]。相关研究表明，干旱和半干旱地区，在土地利用类型以林草地为主体的小流域中，侵蚀强度大部分被控制在中度侵蚀以下[24]，陈浩[25]的研究证实这一现象。据2019年鄂尔多斯市水土保持公报公布数据，2019年准格尔旗有林地面积3 276.19 km2，草地面积2 835.85 km2，林草面积占总面积的80.94%，中等覆盖以上面积为3 908.41，占总面积的42.12%，中度以下水土流失面积2 334.82 km2，占水土流失总面积的72.95%。本研究区所在的皇甫川流域，中度以下水力侵蚀面积为1 016.78 km2，占水力侵蚀面积的73.22%[26]。据2015年内蒙古自治区水土保持公报公布数据，准格尔旗圪坨店坡面标准径流观测场油松乔木林小区(15°)的侵蚀量为6 t/km2，人工草地小区(15°)的侵蚀量为48 t/km2，沙棘林小区(15°)的侵蚀量为42 t/km2，裸地小区的侵蚀量为3 007 t/km2[27]。本研究也得出类似结果，以鲍家沟小流域为代表的裸露砒砂岩区典型流域，在现有的土地覆被条件下，通过RUSLE计算出的土壤侵蚀强度以微度侵蚀为主，结合研究区内径流小区2019年的侵蚀量数据，油松林167 t/km2，沙棘林148 t/km2，山杏林204 t/km2，柠条林297 t/km2，天然草地964 t/km2，裸地3 185 t/km2。表明砒砂岩区的典型小流域，在退耕还林和生态移民政策的持续推动下，土地利用格局随之改变，充当侵蚀“源”的耕地、裸地被大面积的人工林所取代，原有植被景观被割裂成零散的斑块，特别是沙棘植被，即使在陡峭和基岩裸露的坡面也能展现强大的适应性，促成了连续的人工林景观的出现，从而导致了流域侵蚀强度的转变。

邱扬[28]和沈中原[29]的研究发现，黄土高原地区小流域土壤侵蚀具有明显的空间分布规律，与土地利用格局和流域地形因素关系密切，本研究对现有土壤侵蚀模数的空间自相关性进行研究后发现，土壤侵蚀各向异性的空间自相关程度与流域主沟和支沟径流汇集的方向吻合。在丘陵沟壑地貌类型区，地形因素决定了流域径流汇集和土壤搬运的方向[4]，在坡度陡峭且植被景观破碎的沟道中，径流冲刷动力极强，地表覆盖不均匀，是土壤侵蚀最为剧烈的区域。植被景观集中连片的地带，通常也是地形相对平缓的坡面和川台地带，径流冲刷能力相对较弱，土壤侵蚀强度低。本研究中土壤侵蚀局部空间自相关性的结果也证实了这一点，土壤侵蚀模数高值聚集区的空间自相关性与地形因子的关联程度最高，低值区也多聚集在坡面上连片的植被景观带上，土高低值异常聚集区与景观格局因子的相关程度最高，这是因为在大面积的裸露基岩景观中仍存在着水热条件相对较好的地带，这些地带零散分布着植被的景观斑块，这也造成了本该是连片剧烈侵蚀的地带，仍有个别的微度侵蚀区域存在，这也证实了植被的空间格局对土壤侵蚀调控作用。综上，地形因素是土壤侵蚀模数高值聚集区形成的主导因素，植被格局因素是导致土壤侵蚀模数低值聚集区形成的主要因素。

4 结 论

(1) 以鲍家沟小流域为代表的裸露砒砂岩区典型流域，通过RUSLE模型估算出的研究流域土壤侵蚀模数的最高值为9 127.45 t/(km2·a)，小流域土壤侵蚀强度总体上以微度侵蚀为主，沟坡等基岩大面积裸露地区，侵蚀作用依旧剧烈。流域土壤侵蚀模数的全局Moran指数Z值得分为136.87，土壤侵蚀模数在空间上随机聚类的可能性小于1%(p<0.01)，具有显著的空间自相关性。

(2) 鲍家沟流域景观类型以人工植被景观为主，主要有沙棘林、柠条林、山杏林和油松林和草地5种，优势景观为裸露基岩景观和沙棘林景观。

(3) 土壤侵蚀模数的局部Moran指数与地形和植被景观格局因子关联程度较高，地形因素是土壤侵蚀模数高值聚集区形成的主导因素，植被格局因素是导致土壤侵蚀模数低值聚集区形成的主要因素。