黄土高原区域样带小流域泥沙连通性的空间变化

2.1 研究区域

黄土高原(33.72?-41.27?N, 100.90?-114.55?E)位于中国黄河中上游，是中国第二大高原，也是世界上最大的黄土矿床。流域面积64万平方公里，海拔85~5210米(图1)。这个地区属于大陆性季风气候。年平均降水量具有明显的时空异质性，从东南部的800 mm减少到西北部的150 mm, 60-80%在夏季以暴雨形式下降。年平均气温范围为0~13℃，由东南向西北递减。黄土高原的地形分为黄土台地和丘陵沟壑区。本研究在丘陵沟壑区进行，该地区为典型的黄土地貌，受水蚀和风蚀影响较大，逐渐形成了破碎的景观，沟壑发育。沟壑密度范围为4.5~6.0kg km ^?2 。

主要的土壤类型是黄土，由第四纪泥沙发展而来，其特点是物质组成均匀，松散，易受侵蚀。其深度从30米到180米不等，易受侵蚀、疏松、贫瘠。土壤质地东南部较细，西北部逐渐变粗。从东南到西北，植被范围从森林、森林草原、草原和沙漠草原带。然而，大部分原始植被已被不可持续的农业实践破坏，现在主要是由人工植被覆盖。为了控制严重的水土流失，1999年中央政府实施了著名的“退耕还林”工程。该项目的主要措施是将陡峭的农田改造为草地、灌丛或林地。因此，黄土高原植被覆盖率从36%增加到63%。植被的快速恢复肯定会改变植物群落的近土壤表面特征，即凋落物、生物结皮、根系和土壤理化性质，这可能会影响小流域的泥沙连通性及其空间格局。因此，选取6个典型小流域作为南北样带，研究黄土高原小流域泥沙连通性的空间变化规律，并进一步确定其影响因素。所选的小流域的排水面积从36.5到72.8平方公里，平均坡度从7.06°增加到19.64°。年平均降水量由南(620.23 mm)向北(365.73 mm)递减。所选6个小流域的详细信息见表1。

图1. 黄土高原试验小流域空间分布特征。LYG、NXHG、YG、JYG、XHG和WLG分别是罗玉沟、南小河沟、洋沟、九源沟、小河沟和乌兰沟。

2.2 数据收集

本研究使用的数据集包括蒸发量(ET)、湿度指数(HI)、潜在蒸发量(PET)、降水(PRE)、降雨侵蚀力(R)、增强植被指数(EVI)、叶面积指数(LAI)、归一化植被指数(NDVI)、高程(ELE)、坡度(S)、地表粗糙度(SR)、地形湿润指数(TWI)、土壤结构稳定指数(SSI)、土壤可蚀性因子(K)，各指标的详细描述及计算如表2所示。利用中分辨率成像光谱仪(MODIS)数据集获取的ET(MODIS16A2)、PET(MODIS16A2)、EVI (MOD13A1)、LAI (MOD15A2H)和NDVI (MOD13A1)数据集，检索国家航空航天局(nasa)地球观测系统(https://modis.gsfc.nasa.gov) 2000-2020年的数据集。利用MODIS重投影工具(MRT)对下载的MODIS数据集进行格式转换、投影转换和掩码裁剪预处理。MODIS (500 m × 500 m)的空间分辨率足以在小流域尺度上确定上述参数，已有研究证实了这一点(Yang et al.， 2010;Hammond等人，2018)。ET、PET和LAI的时间分辨率为8天，EVI和NDVI的时间分辨率为16天。为了消除云和气溶胶污染的影响，对每个像素采用最大值复合法确定EVI、LAI和NDVI。

年降水量从中国国家气象局(http://data.cma.cn)获得，并进一步插值以覆盖研究流域。降雨R因子和土壤可蚀性K因子来源于国家地球系统科学数据中心和国家科技基础设施(http://www.geodata.cn)。空间分辨率分别为250 m × 250 m和30 m × 30 m。先进陆地观测卫星(ALOS) PALSAR图像，从阿拉斯加卫星设施分布式活动档案中心(ASF DAAC, https://searchasf.alaska.edu/#/)下载，用于测定ELE、S、SR和TWI，空间分辨率为12.5 m。该产品的局限性是不能完全代表裸地表面，因为它是由ASTER GDEM V3复制的，原始分辨率为30 m。土壤颗粒分布和有机质(SOM)数据来源于中国土壤数据库(http://soil.geodata.cn/)，尺度为1:1M。

表1. 6个被测小流域基本情况

表2. 所选流域环境因子的缩写和描述

2.3  方法

2.3.1 泥沙连通性

连通性指数(IC)，首先由Borselli等人提出(2008)，利用ArcGIS 10.6软件根据上坡产径流概率和下坡输沙概率计算。对应的方程为(Borselli et al.，2008):

其中，Dup和Ddn分别为上坡贡献区和下坡贡献区(m2)， W为上坡贡献区平均权重(无因次)，S为上坡贡献区平均坡度(m m?1)，A为上坡贡献区(m2)， di为第i个单元的下坡流道长度(m)， Wi和Si (m m?1)分别为权重因子和第i个单元的坡度。在集成电路计算中，采用D8算法计算流动方向。

本研究采用源汇景观指数(SSLI)计算W。SSLI指标是从源汇理论发展而来的。为确定W，将通用水土流失方程的C因子和K因子分为Csource和Csink, Ksour和Ksink，分别表示水文响应单元(HRU)中的源汇效应。具体算法为:

其中，Csour和Csink是C因子的源效应和汇效应，由NDVI计算，使用Van der Knijff等人(2000)提出的方程。Ksour为K因子的源效应(t hm2 h hm?2 MJ mm?1)，Ksink为K因子的汇效应，K为平均值K (t hm2 h hm?2 MJ mm?1)，S为HRU的斜率梯度(m m?1)，A为HRU的面积(m2)。

将归一化的SSLI作为W，计算如下:

其中，SSLInor是标准化的SSLI, SSLImax是被测小流域内的最大SSLI, SSLImin是最小SSLI。详细的程序可以在我们之前的论文中找到。

2.3.2. 空间自相关

通过空间自相关分析确定Moran’s Index (Moran’s I)，该指数用于检验IC的整体分布，并检测是否存在空间上的聚集特征。Moran’s I取值范围为?1 ~ 1。Moran’s I越大，正自相关越强，属于典型的聚类空间分布格局。p值用来表示统计显著性，显著的Moran’s I表示向聚类的空间模式(Li et al.， 2010)。Z评分也用于确定显著性。Z值为正反映了泥沙连通性在空间上的正自相关。利用ArcGIS 10.6 (https://www.esri.com)对泥沙连通性进行空间自相关分析。

2.4. 冗余分析的方差划分

方差划分是一种解决不同解释变量与响应变量之间的解释能力的方法。考虑到的14个环境因素被分为三组。第一组为气候变量，包括ET、HI、PET、PRE和R；第二组为植被变量，包括EVI、LAI和NDVI。第三组是地形和土壤特征，如ELE、S、SR、TWI、K和SSI。采用方差分割冗余分析(RDA)计算三个环境因子组对每个流域IC的相对贡献。通过不同因素组的相对贡献来量化环境变量对IC的影响。从环境因子组中确定了影响泥沙连通性空间变化的方差比例。从环境因子组中确定了影响泥沙连通性空间变化的方差比例。泥沙连通性的变化可以用每组环境变量及其相互作用来解释。基于500个随机重采样数据，采用蒙特卡罗置换法评价RDA的统计学显著性，p < 0.5为有统计学意义。

3.1 泥沙连通性的空间变化

2000-2020年计算的平均IC在每个测试小流域内表现出显著的空间异质性(图2)。虽然这些小流域位于不同的地方，但小流域内IC的空间分布相似。在流域出水线附近IC较高，而在流域边界附近IC较低(图2)。所选6个小流域的计算Moran’s I均大于0，且由南向北波动，且有较弱的上升趋势。这些结果表明，在黄土高原，IC具有空间依赖性，相关性随纬度的增加而增加(表3)。IC在区域横断面上有明显的空间变化。黄土高原自南向北的平均IC显著线性增加(R2 = 0.85, p < 0.01)(图3)。在6个测试流域中，LYG的平均IC最低(?2.02)，XHG的平均IC最高(2.30)。此外，在黄土高原地区，IC的空间分布密度在南部地区比在北部地区更为分散(图3)，这表明小流域内IC的空间异质性随纬度的增加而减弱。

3.2  环境因子的空间变化

可能影响泥沙连通性的区域样带小流域环境因子的空间变化如图4所示。黄土高原的ET、HI、PRE和R由南向北分别降低了59.17%、52.38%、40.86%和57.58%，而PET由南的1152 mm a - 1增加到北的1423 mm a - 1(图4a、b、c、d和e)。同样，植被相关因子(EVI、LAI和NDVI)由南向北分别降低了40.35%、59.28%和28%(图4f、g和h)。相比而言，地形和土壤因子的空间变化更为多样。在黄土高原，ELE、TWI和SSI由南向北先减小后增大，分别在970.55 ~ 1464.40 m、4.86 ~ 6.4和1.20 ~ 4.29之间(图4i、l和n)。S和SR在黄土高原由南向北先增大后减小，在6.41 ~ 19.73?和1.04 ~ 1.08之间(图4j和k)。黄土高原土壤可蚀性k因子在0.013 ~ 0.018之间波动(图4m)。这些结果表明，这些环境因子的空间格局存在显著差异。

3.3.   泥沙连通性与环境因子的关系

Pearson相关分析结果揭示了影响泥沙连通性的潜在环境因子(表4)。IC与ET (r = - 0.66)、HI (r = 0.65)和ELE (r = 0.63)的相关性最高，其次是PET (r = - 0.59)、LAI (r = - 0.57)、EVI (r = - 0.56)、PRE (r = 0.53)、NDVI (r = - 0.52)和S (r = 0.50)。总体而言，除降雨侵蚀力(R)和土壤可蚀性(K)外，不同小流域间环境变量与土壤侵蚀力(IC)的显著相关(p < 0.05)基本相似，南部流域IC与土壤侵蚀力(R)的回归系数最大(0.69)。北部流域逐渐下降至0.16。也就是说，降雨侵蚀力对黄土高原泥沙连通性的影响由南向北递减。在一些位于所选样带北部的测试小流域中，降雨侵蚀力不是影响泥沙连通性的显著因素。相反，IC与K的相关性从南部的0.15上升到北部的0.29。在靠近南部的一些小流域，相关性不显著。IC与HI、PRE、R、ELE、S、SR、TWI、K呈显著正相关关系，与ET、PET、EVI、LAI、NDVI、SSI呈显著负相关关系(表4)。总的来说，不同小流域泥沙连通性与环境因子的相关性相似。但相关程度不同(表4)。这一结果表明，不同位置的小流域影响泥沙连通性的主导环境因子不同，因此可能需要单独分析。

一些环境因子之间也存在显著的相关性(图5)。一般来说，属于同一类的环境因子，如气候因子(ET、HI、PET、PRE和R)和植被因子(EVI、LAI和NDVI)之间的相关性较高，达到统计学显著水平。不同环境因子之间的相关系数随纬度的增加而降低，尤其是陆面变量。不同环境因子之间的相关性可能影响了它们对泥沙连通性的影响，从而影响了选定区域样带上IC的空间变化。

图2所示 LYG (a)、NXHG (b)、YG (c)、JYG (d)、XHG (e)、WLG (f)流域IC空间格局。流域简称见表1。

表3 Z值和Moran’s I反映了被测小流域IC的空间自相关性(**在0.01水平上显著)。

图3. 黄土高原受试小流域IC的空间变异.

不同的小写字母表示小流域间IC的差异(Duncan, p≤0.05)。流域的缩写见表1。

3.4 影响IC空间变化的主要因素及其贡献

RDA的变化分区分析结果表明，不同位置的小流域影响泥沙连通性的主导因子不同(图6)，总体上，环境因子对泥沙连通性变化的解释在南部小流域大于北部小流域。小流域LYG、NXHG、YG、JYG、XHG和WLG的IC变化量分别为90.8%、88.0%、80.5%、75.6%、66.3%和70.1%。这些结果表明，这些环境因子及其相互作用对泥沙连通性的影响随着纬度的增加而减小。在本研究中，由于6个被测小流域的地质地貌条件相同，因此可以忽略地质地貌因素的潜在影响。此外，所选小流域IC的变化主要是由植被和气候因子组解释的。但随着纬度的增加，其解释作用减弱(图6)。相反，地形和土壤因素对IC的影响由南向北增大。

4.1 泥沙连通性的空间变化及其影响因素

泥沙连通性与流域径流输沙直接相关，因此其空间变化必然与影响径流和泥沙过程的环境因素密切相关。本研究结果表明，沿区域样带所测小流域的IC从南部的- 2.01增加到北部的1.76(图2和图3)。黄土高原河流总悬浮含沙量随纬度的增加而减少的观测结果证实了这一结果。这一结果与之前几项研究的发现一致，其中IC的空间变化直接由高度分化的气候条件、不同的植被性质和复杂的地表特征(包括地形和土壤)造成的径流和泥沙过程的空间异质性引发。

气候特征(如降水、湿度和蒸发量)被认为是促进土壤侵蚀和产沙的主要外部因素，其空间格局是影响流域泥沙连通性空间变化的主要因素。结果表明，土壤IC的空间变化受降水、蒸发量等气候因子的影响。IC与ET、PET呈负相关(p < 0.01)，与HI、PRE、R呈正相关(p < 0.05);这一结果与Wang and Alimohammadi(2012)和Haghighi et al.(2018)的结果一致，他们的研究表明气候特征(即降水、蒸发)的变化直接影响径流和泥沙过程，从而触发了流域泥沙连通性的变化。但气候特征与泥沙导流率的相关系数随纬度的增加而降低(表4)。黄土高原南部降水和降雨侵蚀力较高，产生的径流较多，侵蚀严重，输沙量大，因此区域样带南部小流域泥沙连通性大，与PRE和R关系密切。同样，沿样带辐射差异引起的大气湿度和蒸发量的变化也在一定程度上影响了IC的空间变化。特别是，当所选小流域的地理位置从半干旱区向干旱区移动时，气候条件的变化进一步降低了泥沙连通性的空间异质性(表3)。

IC的空间变异还受到植被的影响。植被在减少径流和侵蚀方面的作用与植被类型、物种和覆盖度的空间分布密切相关。本研究结果表明，在测试的小流域中，植被因子与泥沙连通性均呈负相关(表4)。这一结果与Wei等（2010)的结果一致。Li et al.(2016)和Zhao et al.(2022)，他们发现植被通过对小流域径流和泥沙的影响降低了泥沙连通性。结果还表明，在黄土高原，植被覆盖度与植被因子的相关系数随纬度的增加而减小(表4)。植被覆盖度的空间变化受水热条件的控制。黄土高原植被NDVI和绿化率均由南向北递减。众所周知，植被在防止土壤侵蚀方面的作用随着植被覆盖的增加而减少。因此，植被特征的具体空间格局影响着黄土高原IC的空间变异(图2)。此外，由于“退耕还林”工程的刺激，土地利用类型从陡峭的农田向草地和林地的显著变化，黄土高原小流域泥沙连通性的空间变化可能增强。

沿着测试样带的不同地形和土壤通过改变径流和泥沙之间的关系促进了泥沙连通性的空间异质性。研究结果表明，沿区域样带的小流域IC随ELE、S、SR、K和TWI的增加而增加，随SSI的增加而减小(表4)。地形和土壤对IC的影响随纬度的增加而增强。充分发育的沟谷极大地增强了小流域部分地区的泥沙连通性，促进了泥沙连通性的空间变异。高程和地表粗糙度的空间变化也会影响土壤侵蚀和泥沙沉积，从而影响泥沙连通性的空间异质性。此外，黄土高原南部土壤相对细粒，北部土壤相对粗粒(Wang et al.， 2021)，导致黄土高原南部土壤可蚀性低，北部土壤可蚀性高，且IC与K因子呈正相关。除上述环境因素外，土壤侵蚀类型的空间分布、土地利用变化和水土保持措施(如拦河坝和梯田)也对黄土高原泥沙连通性的空间变化有较大影响。

图4 影响泥沙连通性的环境因子的空间变化。(a) ET为蒸散量，(b) HI为湿度指数，(c) PET为潜在蒸散量，(d) PRE为降水，(e) R为降雨侵蚀力，(f) EVI为增强植被指数，(g) LAI为叶面积指数，(h) NDVI为归一化植被指数，(i) ELE为高程，(j) S为坡度，(k) SR为地表粗糙度，(l)TWI为地形湿润指数，(m) k为土壤可蚀性，(n) SSI为土壤结构稳定性。误差条为1.5标准差。误差条上不同的小写字母表示小流域间IC的差异(Duncan, p≤0.05)。流域的缩写见表1。

表4 测试小流域IC与环境因子的相关性。

图5所示。(a) LYG、(b) NXHG、(c) YG、(d) JYG、(e)   XHG、(f) WLG流域环境因子之间的相关系数。**在0.01水平上显著相关，*在0.05水平上显著相关。流域缩写可以在表2中找到。

4.2 影响泥沙连通性空间变化的主要因素

确定泥沙连通性的主要影响因素对于理解下垫面条件的一系列变化的复杂影响以及制定土壤侵蚀缓解策略至关重要。RDA的结果表明，影响泥沙连通性空间变化的主要因素是植被和气候(图6)。这一结果与植被及其形态在减少径流和侵蚀方面作用显著的共识一致，并被“退耕还林”以来黄土高原土壤侵蚀强度的大幅下降直接证实。应该指出的是，这里讨论的结果特别集中在小流域。在这个尺度上，有许多强有力的证据表明，植被在减少泥沙通量方面具有重要作用。此外，结果还表明，黄土高原植被对泥沙连通性的影响由南向北递减(图6)，这是由于植被盖度、生物量随纬度的增加而减少所致。由于黄土高原植被盖度的差异，土壤侵蚀北部比南部严重。植被对泥沙连通性的强烈影响可以用来解释为什么黄土高原南部高降雨侵蚀力小流域的IC小于北部低降雨侵蚀力小流域(图2和图3)。Zheng等(2021)关于黄土高原河流含沙量南低北高的结果间接证实了我们的结果。

此外，RDA的结果还表明，在黄土高原，气候和植被因子对IC空间变化的贡献逐渐减弱，而陆面因子对IC空间变化的影响随着纬度的增加而增强(图6)。这些结果是由受测环境因子的空间变化引起的。降水和植被覆盖度南高北低，土壤质地西北粗东南细。RDA结果表明，环境因子对泥沙连通性的解释由南部的90.8%下降到北部的70.1%。这种差异可能是由于黄土高原不同地区土壤侵蚀类型的差异造成的，随着纬度的增加，侵蚀类型由水蚀过渡到水-风混合侵蚀。显然，水蚀的影响因素与水蚀风蚀的影响因素有很大的不同。然而，本研究所考虑的环境因素仅与水侵蚀有关。因此，对黄土高原北部小流域泥沙连通性空间变化的解释相对较差。在水和水-风固定侵蚀控制的区域，需要进一步的研究来充分确定影响IC空间变化的因素。

图6所示。LYG、NXHG、YG、JYG、XHG、WLG流域的变化分区(a) ET、HI、PET、PRE和R气候因子;(b) EVI、LAI和NDVI植被因子;(c) DEM、ELE、S、SR、TWI、K、SSI的地形土壤因子。(d)、(e)和(f)是(a)和(b)、(a)和(c)、(b)和(c)之间的相互作用;(g)是(a)， (b)和(c)之间的相互作用。

本文以黄土高原小流域为研究对象，分析了小流域泥沙连通性的空间变化规律，确定了影响小流域泥沙连通性的环境因子。IC量化的泥沙连通性由南向北增加，范围为?2.02 ~ 1.76。小流域泥沙连通性空间变化存在显著的空间正自相关(p < 0.01)，且随纬度增加而增大。大多数环境因子在区域样带的小流域中也表现出显著的空间异质性。IC与环境因子ET、HI、PET、PRE、R、EVI、LAI、NDVI、ELE、S、SR、TWI、SSI、K在0.01或0.05水平上具有显著相关性。受环境因子空间异质性的影响，黄土高原环境因子对泥沙连通性的影响表现出明显的空间差异性。此外，影响泥沙连通性的主要因素是植被。气候因子、植被因子和陆面因子对泥沙连通性的相对贡献率分别为40.9 ~ 71.2%、45.5 ~ 80.3%和7.9 ~ 49.3%，平均值分别为57.8%、61.6%和34.8%。在黄土高原，气候和植被因子对泥沙连通性的影响逐渐减弱，而陆面因子对泥沙连通性的影响从南向北逐渐增强。