一种考虑道路影响的流域泥沙连通性评价方法

1.本发明属于生态水文技术领域，尤其涉及一种考虑道路影响的流域泥沙连通性评价方法。


背景技术：

2.泥沙连通性是指流域不同地貌或景观单元间的泥沙级联关系，用于表征流域内泥沙克服区域阻力通过分离和输移从源到汇区(从坡面到沟道)的运输程度，从而衡量泥沙在地貌单元运动的潜能。这既取决于地貌或景观单元的空间格局，也受控于气候、地貌、生态、水文、土壤侵蚀等地表过程及其时空变异。
3.泥沙连通性是近十年来国际上的研究热点，与地形地貌、生态水文、土壤侵蚀及水质保护等多个研究领域密切相关，它强调不同地貌或景观单元间的泥沙输移、分离和沉积等过程，是控制流域景观演变的最重要属性之一，可用于理解土壤侵蚀与泥沙输移过程及其时空分布。泥沙连通性重点关注流域内不同地貌或景观单元间的泥沙级联关系，是决定流域泥沙输移方式、效率和规模的重要因素，是探究泥沙来源、识别水土保持重点区域的基础，在径流产生、泥沙运移相关的水文过程中起着重要作用。因此，开展流域泥沙连通性研究将有助于揭示流域水沙运移过程变化的影响机制，进而为流域生态修复提供科学依据。
4.近年来，由于人类活动强度加剧，原生地带性稳定群落退化严重、森林资源锐减、水土流失问题日趋严峻。特别是，受森林采伐、旅游开发、矿产开采等因素的影响，道路网络持续扩张，如对原有林区采伐废旧道路的恢复、现有道路的翻修和新规划道路的建设等，进一步加剧了水土流失、改变了流域原有的水沙输移过程。泥沙连通性作为识别与量化泥沙输移路径及其时空变化的基础，在探究泥沙来源热点、识别水土流失治理重点区域等方面发挥着重要作用，且与地貌演变、水文过程、景观生态、土壤侵蚀、水质演变及水资源管理等密切相关。因此，研究泥沙连通性具有重要的理论意义与应用价值。
5.野外调查法作为定性的研究方法，需要在发生明显的次降水或侵蚀事件后在流域内调查土壤侵蚀及泥沙沉积等情况，且需要详细的遥感影像、航片和土地利用等资料，但定性的判断只能用来评价或解释定量化指标的计算结果。
6.近年来，随着科技手段的不断进步，能表征连通性的数据的分辨率得到大大的提升。基于高精度数据，量化泥沙连通性的方法亦逐渐被提出，大致分为图论法、指标法和模型模拟法。如heckmann等人应用图论法分析了阿尔卑斯山脉中部流域泥沙级联及其路径的空间分布。高海东等人基于图论的原理研究了淤地坝如何对径流造成影响。这种方法给泥沙通道的网络结构以及上下坡泥沙连通性的分析提供了极大的帮助。然而，它只能适用于泥沙连通性的空间格局分析，对输沙能力的分析较为缺乏。指标法中broselli等人开发的连通性指数(ic)应用最为广泛，如foerster等将其用于评估西班牙土地利用变化对连通性的影响。模型模拟法中，lesschen等人应用lapsus模型发现植被斑块和农业梯田的空间分布是影响流域泥沙连通性的主要影响因素。但是这些方法在针对道路影响下泥沙连通性的空间分布及其变化的影响尚不明确。尤其是，道路网络在改变泥沙传输方向的同时捕获部
分泥沙和径流，使其沿道路延伸方向向下游传输，并在道路与河流交汇点向流域出水口传输，进而改变原有的泥沙传输路径。
7.综上，图论法针对不同的研究对象需要选择不同的定量指标进行表达，虽能直观展示不同景观单元间泥沙级联的空间分布与动态变化，但适用的空间尺度有限，在较大尺度的流域上对小区、坡面、切沟和细沟等较小景观单元的连通性表达能力有限。指标法和模型模拟法都可用于模拟流域侵蚀泥沙源和汇，定量化表征泥沙连通性，尤其是基于连通性指数的指标法是目前应用最为广泛的研究方法。其中，连通性指数能表征流域泥沙源-汇间的潜在联系。但上述技术在评估流域泥沙连通性时均较少考虑道路作为传输通道对泥沙传输方向、汇流路径长度、有效汇流面积等的影响。


技术实现要素：

8.本发明针对泥沙连通性评价过程中，现有方法未考虑适应山区森林区域的土地覆盖和地表特征，以及道路作为传输通道对流域流向及径流与泥沙传输路径的影响，从而导致泥沙连通性评价过程中道路影响的准确性差的问题，提出一种考虑道路影响的流域泥沙连通性评价方法，基于指标法中的连通性指数(ic)，在计算过程中考虑适用于山区森林流域的相对平滑度指数(rs)作为权重因子；考虑道路作为传输路径对流向的影响，并在影响流域河网的基础上影响汇流路径长度、上坡有效集水区面积等，使得计算结果既能反映山区森林流域的土壤覆盖与地形特征，又能充分反映道路网络对流域泥沙连通性的影响。
9.为了实现上述目的，本发明采用以下技术方案：
10.一种考虑道路影响的流域泥沙连通性评价方法，包括：
11.步骤1、基于arcgis平台，进行数据的获取和处理；所述数据包括数字高程模型dem数据、谷歌卫星影像数据、根据谷歌卫星影像数据目视解译得到的流域内道路网络与河流网络；
12.步骤2、利用arcgis平台，基于流域内dem数据，运用空间分析中的坡度工具计算得到研究区坡度因子；
13.步骤3、计算相对平滑度指数作为权重因子，以反映山地森林地区土地覆盖和地表特征对泥沙连通性的影响；
14.步骤4、以谷歌卫星影像数据矢量化得到的流域内道路网络为基础，基于arcgis平台，运用要素转栅格工具和重分类工具区分流域内的道路和其它土地利用类型，并将道路栅格赋值为0，其它土地利用类型赋值为1，进而生成道路掩膜图层；
15.步骤5、以流域内dem数据为基础，基于arcgis平台，获取水流流向与汇流累积量图层；
16.步骤6、通过空间分析中的水流长度工具计算沿流路径从每个像元到栅格边上的汇点或出水口的下坡距离，并以流域泥沙传输的最终流向作为水流流向的栅格数据、作为权重栅格数据，从而得到下坡组分；
17.步骤7、借助汇流累积量工具分别以相对平滑度指数和坡度为权重栅格数据，计算不同权重栅格数据下的汇流累积量，根据不同权重栅格数据下的汇流累积量，通过栅格计算器工具计算上坡平均权重因子和上坡平均坡度，从而得到上坡组分；
18.步骤8、通过栅格计算器工具中的log
10
函数，计算得到泥沙连通性指数。
19.进一步地，所述步骤1中按照如下方式得到流域内道路网络与河流网络：
20.利用arcgis平台，基于2m分辨率的谷歌卫星影像数据，在catalog下新建要素类型为折线的矢量图层，调用编辑器在编辑命令下描绘谷歌卫星影像数据中呈现的所有道路，描绘完成后点击编辑器下的“保存编辑内容”与“停止编辑”，即生成道路网络；重复上述步骤，生成河流网络。
21.进一步地，所述步骤2还包括：
22.通过栅格计算器工具为坡度＝0的栅格字段赋值为0.005。
23.进一步地，所述步骤3包括：
24.以流域内dem为基础，通过spatial analyst tools中的表面分析-坡度，计算流域内的坡度图层；随后，借助栅格计算器，通过地面粗糙度的计算公式，即内的坡度图层；随后，借助栅格计算器，通过地面粗糙度的计算公式，即得到粗糙度图层；然后，根据粗糙度指数与曼宁系数的关系式，计算得到曼宁系数：
25.y＝2.4234
×e0.3005ln(n)
26.式中：y为地表粗糙度，n为曼宁系数；
27.随后，按照下式得到相对平滑度指数：
[0028][0029]
式中：n
min
为曼宁系数的最小值；rs为相对平滑度指数，是一个无量纲数字，从0.0～1.0不等。
[0030]
进一步地，所述步骤5包括：
[0031]
步骤5.1、以流域dem数据为基础，利用arcgis平台的填洼工具先对dem进行填洼处理，随后运用空间分析工具中的水流流向工具提取dem填洼后的水流方向，即flowdir1；
[0032]
步骤5.2、通过栅格计算器工具，将道路掩膜图层叠加到起始flowdir1上，即道路掩膜*flowdir1，从而得到路网影响下的水流流向，即flowdir2；
[0033]
步骤5.3、基于路网影响下的flowdir2，运用空间分析中的汇流累积量工具计算栅格单元所流经的水量大小，即汇流累积量flowacc1；
[0034]
步骤5.4、根据flowacc1，通过栅格计算器工具，结合矢量化河网确定阈值，从而提取流域河网掩膜图层；再次通过栅格计算器，将河网掩膜图层叠加在flowdir2上，得到流域泥沙传输的最终流向，即flowdir3；
[0035]
步骤5.5、利用汇流累积量工具得到水流流向fowdir3的汇流累积量，即flowacc2；利用栅格计算器，使flowacc2每个栅格的汇流累积量值 1，得到最终的汇流累积量图层，即flowacc3。
[0036]
进一步地，所述步骤6包括：
[0037]
通过空间分析中的水流长度工具计算沿流路径从每个像元到栅格边上的汇点或出水口的下坡距离，并以flowdir3作为水流流向的栅格数据、作为权重栅格数据，得到下坡组分，具体公式如下：
[0038][0039]
式中：d
dn
表示径流和泥沙沿路线到达指定流域出水口或沟道的可能性；di为研究单元到泥沙汇的距离；wi是第i个单元的权重因子，此处用相对平滑度指数代替；si是第i个单元的坡度因子；
[0040]
借助栅格计算器工具为上述计算结果为0的栅格单元赋值。
[0041]
进一步地，所述步骤7包括：
[0042]
借助汇流累积量工具分别以相对平滑度指数和坡度为权重栅格数据，计算不同权重栅格数据下的汇流累积量，即flow
rs
和flows；随后根据flow
rs
和flows，通过栅格计算器工具计算上坡平均权重因子和上坡平均坡度，并根据下式、借助栅格计算器，计算得到上坡组分：
[0043][0044]
式中：d
up
表示上坡的连通分量，决定了上坡产生的径流和泥沙向下游流动的潜力；是上坡贡献面积的平均权重因子，代表山区森林流域径流和泥沙在输移过程中可能遇到的阻抗；为上坡贡献区域面积的平均坡度；为上坡的贡献面积。
[0045]
与现有技术相比，本发明具有的有益效果：
[0046]
本发明通过改进连通性指数(ic)的计算，即在计算过程中以能更好地反映山区森林区域的泥沙连通性的相对平滑度指数(rs)作为权重因子，同时，将道路作为“掩膜”图层参与ic的计算过程，从而考虑路网影响下“河网掩膜”图层的生成，进而反映道路网络对地表径流方向的影响、对汇流路径长度的改变等，以弥补道路作为排水通道对泥沙连通性影响的不足，使其适用于道路影响的山区森林区域。本发明可以弥补道路作为排水通道对泥沙连通性影响的不足，使得计算结果既能反映山区森林流域的土壤覆盖与地形特征，又能充分反映道路网络对流域泥沙连通性的影响。
附图说明
[0047]
图1为本发明实施例一种考虑道路影响的流域泥沙连通性评价方法的整体流程图；
[0048]
图2为本发明实施例计算生成的坡度因子图；
[0049]
图3为本发明实施例计算生成的相对平滑度指数图；
[0050]
图4为本发明实施例计算生成的“道路掩膜”图；
[0051]
图5为本发明实施例计算生成的道路影响下的“水流流向”图；
[0052]
图6为本发明实施例计算生成的道路影响下的“汇流累积量”图；
[0053]
图7为本发明实施例计算生成的下坡组分图；
[0054]
图8为本发明实施例计算生成的平均权重因子图；
[0055]
图9为本发明实施例计算生成的平均坡度因子图；
[0056]
图10为本发明实施例计算生成的上坡组分图；
[0057]
图11为本发明实施例计算生成的泥沙连通性指数图。
具体实施方式
[0058]
下面结合附图和具体的实施例对本发明做进一步的解释说明：
[0059]
如图1所示，一种考虑道路影响的流域泥沙连通性评价方法，包括：
[0060]
步骤1、基于arcgis平台，进行数据的获取和处理；所述数据包括数字高程模型dem数据、谷歌卫星影像数据、根据谷歌卫星影像数据目视解译得到的流域内道路网络与河流网络；
[0061]
步骤2、利用arcgis平台，基于流域内dem数据，运用空间分析中的坡度工具计算得到研究区坡度因子；
[0062]
步骤3、计算相对平滑度指数作为权重因子，以反映山地森林地区土地覆盖和地表特征对泥沙连通性的影响；
[0063]
步骤4、以谷歌卫星影像数据矢量化得到的流域内道路网络为基础，基于arcgis平台，运用要素转栅格工具和重分类工具区分流域内的道路和其它土地利用类型，并将道路栅格赋值为0，其它土地利用类型赋值为1，进而生成道路掩膜图层；
[0064]
步骤5、以流域内dem数据为基础，基于arcgis平台，获取水流流向与汇流累积量图层；
[0065]
步骤6、通过空间分析中的水流长度工具计算沿流路径从每个像元到栅格边上的汇点或出水口的下坡距离，并以流域泥沙传输的最终流向作为水流流向的栅格数据、作为权重栅格数据，从而得到下坡组分；
[0066]
步骤7、借助汇流累积量工具分别以相对平滑度指数和坡度为权重栅格数据，计算不同权重栅格数据下的汇流累积量，根据不同权重栅格数据下的汇流累积量，通过栅格计算器工具计算上坡平均权重因子和上坡平均坡度，从而得到上坡组分；
[0067]
步骤8、通过栅格计算器工具中的log
10
函数，计算得到泥沙连通性指数。
[0068]
进一步地，步骤1中数据的获取与处理包括：
[0069]
本发明所使用的数据包括：15m分辨率的数字高程模型(dem)数据、2m谷歌卫星影像数据、根据谷歌卫星影像数据目视解译得到的流域内道路网络与河流网络。
[0070]
流域内道路网络与河流网络的提取步骤包括：利用arcgis平台，基于2m分辨率的谷歌卫星影像数据，在catalog下新建要素类型为折线的矢量图层，调用“编辑器”在编辑命令下手工描绘谷歌卫星影像数据中呈现的所有道路，描绘完成后点击“编辑器”下的“保存编辑内容”与“停止编辑”，即生成“道路”网络。重复上述步骤，生成“河流”网络。
[0071]
进一步地，步骤2中坡度因子(s)的计算过程包括：
[0072]
以流域内高程数据dem为基础，基于arcgis平台，运用空间分析中的坡度工具计算得到研究区“坡度”因子。为避免最小值为0时引起的连通性指数在计算过程中出现不收敛的情况，通过栅格计算器工具为“坡度”＝“0”的栅格字段赋值为0.005，即[(“坡度”＝＝0)*0.05 “坡度”](如图2)。
[0073]
进一步地，步骤3中权重因子(w)的计算过程包括：
[0074]
权重因子是反映流域地表状况和土地利用的加权系数，直接决定泥沙输送的过程和效率。本发明选择相对平滑度指数(rs)作为权重因子，以反映山地森林地区土地覆盖和地表特征对泥沙连通性的影响。该指数由曼宁系数推导出来的，同时保留泥沙连通性指数
(ic)中权重因子(w)在计算过程的无量纲特性。
[0075]
计算过程以流域dem为基础，通过“spatial analyst tools”中的“表面分析
”‑“
坡度”，计算流域内的“坡度图层”。随后，借助“栅格计算器”，通过地面粗糙度的计算公式，即得到“粗糙度图层”。然后，根据粗糙度指数与曼宁系数的关系式，计算得到曼宁系数。具体公式如下：
[0076]
y＝2.4234
×e0.3005ln(n)
[0077]
式中：y为地表粗糙度，n为曼宁系数。
[0078]
随后，借助相对平滑度指数与曼宁系数之间的关系式，计算得到相对平滑度指数rs(如图3)。具体公式如下：
[0079][0080]
式中：n为曼宁系数；n
min
为曼宁系数的最小值；rs是一个无量纲数字，从0.0～1.0不等，保留了ic计算的无量纲性。这里使用rs作为权重因子，能够反映森林土壤表面积累的森林凋落物对径流与泥沙传输的影响。
[0081]
综上，本发明对权重因子进行了改进：原有泥沙连通性指数在计算过程中侧重于考虑植被覆盖管理因子(c因子)对农业和森林环境中径流和泥沙通量的量化，本发明基于对山区森林环境的考虑，选择相对平滑度指数(rs)对权重因子加以改进，改进后的权重因子能更好地反映土地覆盖和地表特征对山区森林流域泥沙连通性的影响。
[0082]
进一步地，步骤4中“道路掩膜”图层的获取过程包括：
[0083]
以谷歌卫星影像数据矢量化得到的道路网络为基础，基于arcgis平台，运用要素转栅格工具和重分类工具区分流域内的道路和其它土地利用类型，并给道路栅格赋值为0，其它土地利用类型赋值为1，进而生成“道路掩膜”图层(图4)。
[0084]
进一步地，步骤5中“水流流向”与“汇流累积量”图层的获取过程包括：
[0085]
步骤5.1、以流域数字高程dem为基础，为避免计算过程中伪地形对径流路径造成的影响，在数据处理时，利用arcgis平台的填洼工具先对dem进行填洼处理。随后运用空间分析工具中的水流流向工具提取dem填洼后的水流方向，即“flowdir1”。
[0086]
步骤5.2、通过栅格计算器工具，将“道路掩膜”图层叠加到起始“flowdir1”上，即“道路掩膜*flowdir1”，从而得到路网影响下的水流流向，即“flowdir2”。
[0087]
步骤5.3、基于路网影响下的“flowdir2”，运用空间分析中的汇流累积量工具计算栅格单元所流经的水量大小，即汇流累积量“flowacc1”。
[0088]
步骤5.4、根据“flowacc1”，通过栅格计算器工具，结合矢量化河网确定阈值，从而提取流域“河网掩膜”图层。再次通过栅格计算器，将“河网掩膜”图层叠加在“flowdir2”上，得到流域泥沙传输的最终流向，即“flowdir3”(图5)。
[0089]
步骤5.5、利用汇流累积量工具得到水流流向“fowdir3”的“汇流累积量”，即“flowacc2”。为避免在后续计算过程中出现零和无穷大的情况，利用栅格计算器，使“flowacc2”每个栅格的“汇流累积量”值做“ 1”处理，得到最终的“汇流累积量”图层，即“flowacc3”(图6)。
[0090]
综上，本发明对径流与泥沙传输通道进行了改进：通过“道路掩膜”图层的叠加，改变径流与泥沙的传输方向，在此基础之上生成“河网掩膜”图层，影响径流与泥沙从源到汇
的传输路径长度，使得到达各个汇区的汇流累积量有所不同。在此过程中，道路的添加如何改变流向、汇流路径长度、汇流累积量，进而影响河网是泥沙连通性计算过程的关键点。这主要包括道路网络的栅格化，并在计算过程中将其考虑到道路影响下径流和泥沙沿着道路方向传输，以及在道路-河流交汇处向流域出水口的输移。
[0091]
进一步地，步骤6中下坡组分(d
dn
)的计算过程包括：
[0092]
通过空间分析中的水流长度工具计算沿流路径从每个像元到栅格边上的汇点或出水口的下坡距离，并以“flowdir3”作为水流流向的栅格数据、作为栅格数据的权重栅格数据，得到“flowlength”，此时即得到下坡组分，具体公式如下：
[0093][0094]
式中：d
dn
表示径流和泥沙沿路线到达指定流域出水口或沟道的可能性；di为研究单元到泥沙汇的距离；wi是第i个单元的权重因子，此处用相对平滑度指数(rs因子)代替；si是第i个单元的坡度因子。
[0095]
为避免后续计算过程中出现无穷大等不收敛的情况以及更准确的表达道路所经处的实际情况，借助栅格计算器工具为上述计算结果为0的栅格单元赋值，最终计算得到下坡组分(d
dn
)(图7)。
[0096]
进一步地，步骤7中上坡组分(d
up
)的计算过程包括：
[0097]
为计算上坡平均权重因子和平均坡度借助汇流累积量工具分别以“相对平滑度指数”和“坡度”为权重栅格数据，计算不同权重栅格数据下的“汇流累积量”，即flow
rs
和flows，此时的“汇流累积量”分别受“相对平滑度指数”以及“坡度”的影响。随后根据flow
rs
和flows，通过栅格计算器工具计算(图8)和(图9)，并根据上坡组分的计算公式(具体公式如下)借助“栅格计算器”，计算得到d
up
(图10)。
[0098][0099]
式中：d
up
表示上坡的连通分量，决定了上坡产生的径流和泥沙向下游流动的潜力；是上坡贡献面积的平均权重因子，代表山区森林流域径流和泥沙在输移过程中可能遇到的阻抗；为上坡贡献区域面积的平均坡度；为上坡的贡献面积。
[0100]
进一步地，步骤8中泥沙连通性指数(ic)的计算过程包括：
[0101]
通过栅格计算器工具中的“log
10”函数，计算得到泥沙连通性指数(图11)，具体公式如下：
[0102][0103]
综上，本发明通过改进连通性指数(ic)的计算，即在计算过程中以能更好地反映山区森林区域的泥沙连通性的相对平滑度指数(rs)作为权重因子，同时，将道路作为“掩膜”图层参与ic的计算过程，从而考虑路网影响下“河网掩膜”图层的生成，进而反映道路网络对地表径流方向的影响、对汇流路径长度的改变等，以弥补道路作为排水通道对泥沙连通性影响的不足，使其适用于道路影响的山区森林区域。本发明可以弥补道路作为排水通道对泥沙连通性影响的不足，使得计算结果既能反映山区森林流域的土壤覆盖与地形特
征，又能充分反映道路网络对流域泥沙连通性的影响。
[0104]
以上所示仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。