重庆市梁滩河汇流阈值分析与地貌发育研究

胡煜佳，郑莉萍，张森林，邵景安，2

(1.重庆师范大学地理与旅游学院，重庆 401331；2.三峡库区地表过程与环境遥感重庆市重点实验室，重庆 401331)

0 引 言

随着GIS技术的不断发展以及在水文模拟中的广泛应用，对流域河网的提取大体分为两种思路：一种是基于已有的遥感影像进行自动提取；另一种则是利用DEM数据进行提取[1]。以DEM和ArcGIS为载体进行河网数字化、自动化提取是众多研究者采取的方法之一；而河网提取以及水文特征信息获取的关键在于每个格网单元的汇流累积量。汇流累积量取决于每个格网单元的流向以及流量分配[2]，采取不同计算流向的方法以及汇流阈值的拟定其自动提取的河网和流域水文特征也不一。单流向算法和多流向算法是研究中常见的两种流向算法[3]，单流向算法如D8算法[4-5]、Lea算法[6]、Rho8算法[7]等；而多流向算法则认为水流具有分散性，分散流入其他下游格网单元[8]，例如FD8算法[9]、Dinf算法[5]等。针对汇流阈值的拟定，学者也采取了不同方法进行研究。左颖等[10]采取多阈值约束法、单一t值确定法、河网密度法综合确定合理阈值；唐杰等[11]采用河网和河源密度法确定集水面积阈值，函数曲线的拐点则为最佳集水面积阈值；张晓娇等[12]采用河网密度法和水系分形法综合确定阈值；许斌等[13]结合分形理论分析了不同汇流累积阈值下河网分形维数和密度变化率特征。综上可看出，大部分学者确定汇流阈值采用频率较高的是河网密度法，并且都用多种方法综合确定。由于基于DEM数字高程模型提取河网以及流域水文信息可行性高、易于操作。因此，本文选取了重庆市污染情况最复杂的次级河流之一——梁滩河流域，基于研究区12.5 m分辨率的DEM，以ArcGIS软件为平台提取流域水文信息，运用数理方法确定最佳阈值，最终提取的河网与TM影像的归一化水体指数提取的河网与三调数据成果提取的河网进行对比分析，并运用Horton算法与起伏比法对流域地貌发育阶段进行识别判定。

1 研究区域概况

梁滩河流域位于重庆主城区缙云山与中梁山之间的北碚丘陵谷地及中梁山的狭长槽谷地带，地貌分区为川东平行岭谷区，地势北低南高，东西高中间凹。该河流发源于九龙坡区走马镇廖家沟水库，流经九龙坡、沙坪坝区、北碚3区，最后汇入嘉陵江，是嘉陵江下游右岸的一条主要支流、主城区流域面积第三大的次级支流，整条流域共有55条支流，干流长约88.7 km，流域面积516.18 km2。流域内常住人口超过20万，并且由于梁滩河流域内曾经发展的养殖业和小型加工业兴盛，其水质受农业活动、工业、城镇污水的影响较大，是全市污染情况最复杂、治理难度最大的次级河流之一。因此，选取该流域为研究对象，提取河网等河流水文信息以期为研究该流域范围内的农业面源污染或水质评价奠定基础。其DEM影像见图1。

2 数据来源与预处理

本文所用的DEM数据来源于重庆市生态环境大数据应用中心，其分辨率为12.5 m，根据梁滩河流域边界对其进行裁剪处理。由于原始DEM存在凹陷单元格以及在平坦地区的水流方向存在不确定性；因此需对原始DEM进行预处理。即填洼操作，削峰填低，创建无凹陷DEM。这是进行流域水文分析的基础。在填洼过程中，限制值Z为指定凹陷点深度和倾泻点间的最大允许差值并确定要填充的凹陷点和保持不变的凹陷点，由于Z值无法精确限定，在本次操作过程中未对Z限制作出任何指定，因此将移除所有峰值。本文所用第三次国土调查(以下简称“三调”)数据来源于重庆市国土局。TM影像来源于地理空间数据云，其分辨率为30 m，利用ENVI5.3对影像进行拼接、大气校正、辐射定标、裁剪等预处理。

3 流域水文特征信息提取

3.1 流域水文特征信息提取

3.1.1 水流方向的确定

在对原始的DEM进行填洼操作后，需要确定每个栅格单元的水流方向，而流向判断主要分为单流向算法和多流向算法[14]。单流向算法目前应用最广泛的则是D8算法。该算法将水流方向简化为8个方向，即东北、东、东南、南、西南、西、西北和北被定义为有效的水流方向，分别用128、1、2、4、8、16、32、64这8个有效特征码表示；再根据高程值从而决定其每个栅格的水流方向[4]。多流向算法如FD8算法[9]、Dinf算法[5]等。本文利用ArcGIS10.2中D8算法计算水流方向，但此时要检验填洼是否完全；若不完全，则需反复填洼。

3.1.2 汇流累积量的计算

每个栅格都会产生一定的水流方向，某一栅格的汇流累积量就意味着从上游流至该栅格的流量累计；因此，通过计算汇流累积量从而判断每个栅格的汇流能力。汇流累积量为栅格单元数目与栅格单元面积乘积之和[1]。

3.1.3 河网提取与分级

河网提取是在计算汇流累积量的基础上进行的。即，将汇流累积量大于某一设定阈值的栅格提取出来的便是河网；因此，阈值拟定对最终提取河网尤为关键。随后再利用水文分析工具集下的Strahler河网分级法进行分级。

3.1.4 栅格河网矢量化与水文特征信息提取

利用水文分析工具集下的栅格河网矢量化工具，将表示线状网络的栅格转换为表示线状网络的要素；再使用河流链接、分水岭工具生成集水区；最后，利用转换工具集下的栅格转面工具将栅格数据集转换为面要素，从而计算不同阈值下的流域面积、河流长度、河道个数等水文信息。

3.2 TM影像的流域水文信息提取

利用遥感影像提取流域水体信息的方法主要有目视解译和计算机自动解译2种[15]。本文利用常用的水体信息提取方法——归一化水体指数(NDWI)进行自动提取。该方法由Mcfeeters[16]在NDWI的基础上提出，利用不同波段对地物反射的特征——近红外波段对水体的强吸收而植被强反射的特点，采用绿波段和近红外波段的比值可以较大程度上抑制植被信息，突出水体，从而较好地提取流域水体(见图2)。这里归一化水体指数

NDWI=(Green-NIR)/(Green+NIR)

(1)

式中，Green为绿光波段反射率；NIR为近红外波段反射率。

从图2可看出，利用ENVI对TM影像的波段识别最终得出水体归一化处理效果不是很理想，对河流以及小面积水域提取效果不佳或存在漏提或断流现象，导致诸类现象的原因可能是该遥感影像分辨率不高，对于地物的识别较模糊，有待采用更高分辨率的遥感影像进行提取与研究。

4 结果与分析

4.1 汇流阈值对流域水文特征的影响

阈值直接影响河网的生成以及流域水文信息，不同汇流阈值条件下其提取出来的河流长度、流域面积、河道个数等有所差异。图3显示，随着汇流阈值的增大河网越稀疏，支流数量越少。

通过属性表计算出不同汇流阈值下的河流长度、流域面积、河网密度以及河道个数(见表1)。

由表1可知，当汇流阈值由10 000增加至30 000时，河道个数由200个减至60个，减少了70%，河流长度、流域面积、河网密度分别降低了46.72%、3.42%、44.83%；当汇流阈值由30 000增加至50 000时，河道个数、河流长度、流域面积、河网密度分别降低了48.33%、25.87%、1.53%、24.72%；当汇流阈值由50 000增加至80 000时，河道个数由31个减至20个，减少了35.48%，其余分别降低了18.29%、1.85%、16.75%；汇流阈值由80 000增加至100 000时，河道个数由20个减至13个，减少了35%，其余分别降低了9.76%、2.7%、7.25%。由以上分析可知，随着汇流阈值的增加流域各水文特征值都在不断降低，但是变幅最小的为流域面积，其他3类特征值波动明显，并且可发现汇流阈值由10 000增加至30 000 时流域各水文特征值变化最为明显，各特征值急剧减小；阈值大于30 000后，各水文特征值变化幅度明显减小。

由于流域面积受阈值波动影响远小于其他特征值，因此，本文分别用对数、幂函数、指数函数、多项式函数等函数对汇流累积量与河流长度、河网密度、河道个数进行趋势分析，寻找最佳函数类型以及最佳阈值，经过多次实验发现幂函数拟合性最优。汇流阈值与河流长度、河网密度、河道个数的幂函数关系如下：

Y1=55 755.948X-0.562(R2=0.997 4)

(2)

Y2=74.479X-0.527(R2=0.995 1)

(3)

Y3=6.909X-1.134(R2=0.998 8)

(4)

式中，Y1为河流长度，km；Y2为河网密度，km/km2；Y3为河道个数；X为汇流阈值；R2为拟合度。

对汇流阈值与河流长度、河网密度、河道个数的幂函数关系求二阶导数，并做非线性拟合与趋势分析(见图4)。

由图4可看出，河流长度、河网密度、河道个数随阈值变化的二阶导数趋势一致，曲线的变化趋势由急剧下降至缓慢减少再至趋于平缓，都是从10 000 增至20 000时，流域各水文特征值急剧下降，20 000增至30 000时，各特征值降幅减小但降幅也较其他值段明显。因此，阈值30 000左右是该流域明显的一个特殊点。汇流阈值与河流长度、河网密度、河道个数的二阶导数关系曲线趋于平缓时所对应的点为合理的汇流阈值。从图4还可看出，当阈值在35 000时为曲线趋于平缓时的变点，其二阶导数值趋向0，各特征值趋于稳定，后与实际河网进行对比验证。即

(5)

(6)

(7)

4.2 提取河网结果的分析

通过最佳汇流阈值提取的河网需要检验与实际河网的差距。本文以三调数据成果提取的河流作验证分析(见图5)，发现通过数理方法确定的最佳汇流阈值35 000提取的河网(见图6)与实际河网较为接近，但仍有一定出入；从DEM中提取的河网很好地保留了河流主干道信息，但存在有一小部分支流不能很好地体现或存在“伪支流”现象，尤其是流域中西部存在部分伪河道。通过其他阈值提取的河网发现，当阈值低于35 000时，河网越密集，提取的河网发育较好，设定的阈值较小时生成的河网与实际河网对比发现很多都是“伪支流”，在实际中并非真正的河网。当阈值高于35 000时，河流长度变短，虽然都较好地保留了主干道信息，但绝大部分支流无法体现，支流发育程度较差，与实际情况也不符合。综上所述，该研究区最佳汇流阈值为35 000，提取的河网在保留河流主干道信息之外其支流也能较好地与实际河网吻合。

4.3 地貌发育阶段识别

4.3.1 Horton算法

Horton算法是依据地形图或数字高程数据利用ArcGIS提取水系，再用Strahler河网分级法对水系进行分级[17]。即

RB=NW-1/NW

(8)

RL=LW/LW-1

(9)

(10)

式中，RB为水系分叉比；RL为河长比；w为水系级别号；NW为第w级河道数目；LW为第w级河道的平均长度；Db为Horton-Strahler水系分维值。

由上述表达式得出梁滩河水系分维值为1.072 4，根据王玉成等[18]和任娟等[19]的相关研究，提出当水系分维值Db≤1.6、1.6

4.3.2 起伏比法

判断地貌发育阶段的方法除了Horton法，还有经典的起伏比法，该方法由Pike和Wilson通过数学方法推导出的估算面积高程积分值的典型方法，其计算过程方便快捷、精度也较高[17]。计算式为

HI=(Hmean-Hmin)/(Hmax-Hmin)

(11)

式中，HI为流域的起伏比；Hmean为流域的平均高程，m；Hmax为流域的最大高程，m；Hmin为流域的最小高程，m。

依据地貌发育理论与Strahler对地貌发育阶段的定量划分(见表2)。梁滩河流域HI值为0.514 9，处于地貌发育阶段的壮年(偏幼)期。由以上2种方法可推断出，梁滩河流域的地貌发育阶段可能为幼龄期向壮年(偏幼)期发展。

表2 地貌发育阶段划分标准

5 结论与讨论

目前，利用DEM和ArcGIS快速提取流域河网的方法广泛适用于各个流域的水文分析，是众多研究者采取的主要方法。本文以重庆市梁滩河为例，利用ArcGIS10.2水文分析模块的相关工具提取河网等水文信息，并通过数理方法确定最佳汇流阈值与流域地貌发育阶段判定，得出了以下结论。

(1)随着汇流阈值的增大，河流长度、流域面积、河网密度、河道个数等流域水文特征值不断减小，提取的河网越稀疏，支流数量越少。

(2)根据汇流阈值与河流长度、河网密度、河道个数的二阶导数关系曲线确定其最佳汇流阈值为35 000。

(3)通过汇流阈值35 000提取的河网与TM影像的水体归一化处理结果、三调数据成果提取的河网进行对比发现该结果具有可靠性，但仍有一定出入。笔者认为原因可能如下：①在使用软件计算与提取的过程中，由于算法的局限性以及某些参数的设定不同，会对运算结果产生很大的影响，如在填洼中Z值的限定、计算流量中汇流阈值的拟定都会对最后提取的结果产生很大的影响，众多学者尤其对于Z值的限定无法做出较好的回答，这便成为日后研究中所需关注的一大重点；②利用TM影像进行水体归一化处理效果并不太理想，存在漏提或断流现象。这可能是由于TM遥感影像分辨率不高，因此有待采用更高分辨率的遥感影像进行提取与研究。

(4)针对梁滩河流域地貌发育阶段的识别，本文通过Horton算法与起伏比法，最终推断出该流域地貌发育阶段可能为幼龄期向壮年(偏幼)期发展。

就总体而言，在流域提取或水文分析中越来越广泛运用GIS技术，相关理论研究和提取方法也日益成熟。但由于DEM数据精度的差异、处理方法不同和人为因素影响等，流域水文特征提取的精度会受到很大影响，如何不断进行技术创新、提高精度是日后研究的重点。