基于DEM与河网密度的水系提取与应用——以澜湄流域为例

曹 爽，王 蒙，包红军，张雨凤

(1.国家气象中心，北京 100081；2.中国气象局-河海大学水文气象研究联合实验室，北京 100081；3.江苏省水文水资源勘测局苏州分局，江苏 苏州 215011)

0 引 言

流域水系和流域边界能有效反映流域地形地貌特征，是流域面雨量监测预报和分布式水文模型构建的重要基础数据，是数字孪生流域建设的主要依据[1]。它们也是完成面雨量实况监测、预报预测及预警服务的关键环节和前提条件，数据精度对相关计算和模拟结果有着直接影响[2]。随着地理信息系统(GIS)和数字高程模型(DEM)的迅猛发展，国内外不断有学者对提高提取流域河网精度和子流域划分进行深入研究和改进，以期得到最接近真实河网的数字化水系。一些学者[3-6]以DEM数据源为基础通过优化改进算法提出了移动窗口法、坡面径流模拟法、谷线搜索法、从DEM直接提取河网与划分子流域的方法等；还有一些学者[7-12]引入矢量河网作为DEM控制条件，通过Agree算法、Burn-in算法、Stream Burning算法等对主要河流进行高程与流向修正，从而解决平原区和受人类活动干扰的城市化的水系提取问题；另外，设置合理的网格数和阈值范围也是控制河网精度的决定性因素之一。有大量研究表明[13-19]，可以通过河网密度法、水系分形维数法、均值变点分析法、河道平均坡降法、流域宽度分布法等方法确定最佳集水面积阈值；在子流域划分方面，有些学者提出了利用河流分汊和拓扑关系[20]、多阈值虚拟河网融合技术[21]等划分子流域，还提出了考虑高山流域[22]、湖泊水库范围[23]、山区平原地貌差异[24]、反映地表下垫面类型特征[25]的子流域划分和编码方法。

在前人研究成果的基础上，本文以澜沧江-湄公河流域(以下简称“澜湄流域”)为研究对象，从遥感影像获取矢量河网数据，采用Agree算法、应用Arcgis软件叠加矢量河网从而修正主要河道DEM高程数据；分析河网密度与集水面积阈值的变化关系，以曲线割线斜率法定量确定最佳集水面积阈值，提取流域水系、边界并划分子流域；基于河网“套合差”法进行精度评价，验证提取河网的可靠性；再以一次降雨过程为例分析流域面雨量的空间分布，将结果应用于流域面雨量监测、预报服务业务，为水文气象预报预警工作提供技术支撑。

1 研究区概况与数据

本文以澜湄流域作为研究区，流域范围介于9°26′～34°2′ N，93°46′～108°52′ E，河长4 900 km左右，流域面积约81万km2，是东南亚地区最重要也是最大的国际河流。该流域发源于中国青海省唐古拉山脉，流经中国、缅甸、老挝、泰国、柬埔寨、越南等六国，以云南省南腊河口为界，中国境内河段称为澜沧江，境外河段称为湄公河，终由越南胡志明市流入南海。澜湄流域地形地势呈北高南低态势，高度落差大，表现出很强的垂直地带性，高程随纬度减小逐渐下降，地貌由北向南依次是高山峡谷区、中低山宽谷区和冲积平原，海拔跨度从3 500～5 000 m降至1 000～3 000 m，再降至1 000 m以下。受地形地貌影响，澜沧江-湄公河表现为典型的南北向狭长型河流，流域形状似“帚状”；上游支流短小且少，呈 “树枝状”、“羽状”水系；中游东岸支流水系发育较好且西岸少有大支流，呈“梳状”水系；下游多为平原和三角洲，河网特别发育，呈“辫状”、“格状”水系特征[26-27]，见图1。此外，在气候变化和人类活动的共同作用下，澜湄流域频发旱涝灾害[28]，坡地山区易发山洪常伴有泥石流，河谷坝区、平原和三角洲的耕地城镇被洪水侵袭，澜沧江中游常有“焚风”效应导致干旱河谷等灾害事件。

图1 澜湄流域概况

本文采用NASA SRTM的空间分辨率为90 m的DEM数据，以及国家基础地理信息中心提供的全球1∶1 000 000包含行政区划、湖泊、水系等信息的矢量数据，而降水数据采用中央气象台5 km×5 km分辨率的智能网格降水。

2 研究方法

2.1 河网提取方法

2.1.1 河网修正Agree算法

DEM数据是高程信息的反演产品，其本身就存在一定误差。在提取水系研究河网分布时无法综合考虑地形地貌、雨水冲刷、人工河道等其他因素的共同作用，会直接影响河网提取精度。尤其是在地势平坦、受人类活动影响大的地区表现明显[9,11]。

澜湄流域下游地区即为地势落差小的冲积平原，且耕地、人工开凿水渠、城市化等人类活动影响大。以DEM作为单一要素提取的河网水系与实际分布会存在较大差异，遂引入Agree算法修正DEM。该算法原理是通过叠加矢量河网，降低实际河网所在栅格的高程，控制水流方向，增加提取河网栅格的“汇流”能力。

算法主要步骤包括：①通过卫星遥感影像解译、Google earth清绘河道中线或国家基础地理信息中心提供的地理信息数据等途径获取所需的实际主要河网数字化水系的矢量文件；②叠加矢量河网和填洼后的DEM，以矢量线要素为中心做缓冲区，缓冲区邻域半径的设置应不小于1/2个DEM分辨率；③运用转换工具将矢量缓冲区范围转为栅格文件，进而运算降低重叠部分DEM栅格高程，获取修正后DEM，命名为DEM-A。

2.1.2 流域水系提取

基于DEM的水系和流域边界提取技术路线见图2。工作原理是基于地表径流漫流模型，关键步骤包括：迭代计算洼地填平；D8单流向算法确定水流流向；沿水流流向计算每个栅格单元的上游汇流能力即汇流累积量；设置汇流阈值，汇流量不小于阈值的栅格是潜在河网；根据Strahler水系分级法对河流进行分级；捕捉倾斜点和计算分水岭，获取自然全流域和自然子流域出口和流域边界；水系河网和流域边界矢量化[14]。

图2 水系和流域边界提取技术路线

2.2 确定最佳集水阈值方法

在流域提取过程中有一个重要的变量参数是集水面积阈值。阈值的变化控制着水系河网的疏密程度，目视判读的方式存在一定的主观性，计算结果难以达成统一。于是，本文采用河网密度法，运用拟合指数函数曲线的割线所对应的斜率与拟合曲线相切的数学方法来确定切点，该切点的物理意义即为河网密度变化趋于平稳时的最佳集水面积阈值。

拟合集水面积阈值与河网密度关系曲线，构建集水面积阈值与河网密度关系方程f(x)，定义域为[x1，x2]，值域为[f(x2)，f(x1)]，令f(x)的一阶导数为f′(x)。在定义域内，曲线割线的斜率k=[f(x1)-f(x2)]/(x1-x2)，以该斜率作为拟合曲线的切线，即k=f′(x0)，求解x0，切线与拟合曲线的切点(x0，y0)代表了河网密度随集水面积阈值由剧烈变化变为平缓变化的转折点。即，x0为目标阈值[16]。

2.3 河网精度评价

本文对河网精度评价分为两个方面：一方面，检验提取的水系和数字化水系主要河道河长的相对误差，有效反映偏离真值的实际大小；另一方面，引入“河网套合差”的概念，指的是叠加提取水系和数字化水系两部分，两水系由于位置偏移会形成细碎的多边形，计算这些细碎多边形面积占流域总面积的比值即为套合差，该比值越小代表两条水系的位置偏差越小[14,29](见图3)。

图3 河网套合差示意

2.4 流域面雨量计算

流域面雨量是指某一流域或区域整个面上的平均降雨量。它能客观地反映流域的降水情况，在水文模型和水情预报等工作中应用广泛。本文的降水资料数据采用的是中央气象台5 km×5 km网格的智能网格预报降水产品，以算术平均法估算流域面雨量预报结果适用可行，计算简单。具体计算原理及公式不再详述[30]。

3 结果与分析

3.1 集水面积阈值与水系河网的变化关系

本文选择1 000、10 000、15 000、20 000、25 000、30 000、40 000、50 000、60 000、100 000共10个河网栅格数为阈值样本，在提取水系河网的过程中以此判断汇流累积量栅格数量的可变参数影响河网疏密程度。不同集水面积阈值对水系提取的影响特征见表1。

表1 受集水面积阈值影响提取水系河网特征变化

由表1可看出，随着阈值面积由8.1 km2逐渐增加至810 km2的变化过程中，河道总数从44 429条减少为433条，河源数从22 261个减少为217个，总河长从198 854.57 km减至23 363.51 km，河道级别从八级河流减少至五级河流；在集水面积阈值扩大了100倍的情况下，河道总数、河源数相应减少为1/100，总河长则减少成了近原来的1/9，流域面积虽受到一定影响有一些减少但变化不大。另外，流域河网分级结果也存在一定变化规律。在同一集水面积阈值下，河道数和河长随河网级别增加逐渐减少且减小幅度逐渐变缓，一级河道数占河道总数50%左右，二级河道数和河长是一级河道数和河长的1/2，一级和二级的河长和河道数在总数的占比达到75%左右。

3.2 流域最佳阈值的确定与精度评价

依据表1中河网阈值栅格数对应的河网密度绘制图4关系曲线。由图4可知，随阈值增大，河网密度逐渐减小，阈值范围为103～104，最大减小幅度呈断崖式；阈值范围为104～4×104，仍呈减小趋势但下降幅度逐渐变缓；阈值范围为4×104～105，变化逐渐趋于稳定。

图4 河网栅格数阈值—河网密度关系曲线

本文引入多种函数关系拟合河网栅格数阈值和河网密度，最终确定采用拟合效果最好的幂函数关系，得到拟合方程

y=6.165 7x-0.465(R2=0999 3)

(1)

式中，y为河网密度，km/km2；x为河网栅格数阈值；R2为相关系数，代表拟合程度，R2=0.999 3表明拟合程度高，两者相关性强。

如图5所示，连接阈值样本首尾两端点做割线，以割线斜率为拟合曲线切线斜率做切线，求得切点即为拟合曲线由剧烈到平缓变化的转折点。对拟合方程(1)求导，得到拟合曲线一阶导数

图5 最佳阈值的确定

k=y′=-2.867x-1.465

(2)

进一步求解方程(2)与斜率k的关系，获得切点坐标为(15 933，0.068)，此时河网栅格数阈值15 933为提取水系的最佳阈值，即最佳集水面积阈值约为129.1 km2，以该阈值提取水系结果如图6所示。切点纵坐标0.068其物理意义为，拟合曲线上河网栅格数阈值15 933对应的河网密度是0.068 km/km2；利用Arcgis水文分析工具以15 933为阈值提取水系，计算得到河网密度为0.069 km/km2，相对误差为1.5%，拟合效果好。

图6 90 m分辨率下澜湄流域最佳集水面积阈值提取水系

为验证提取水系精度，反映提取水系与矢量化数字水系偏移程度，对河流总长的相对误差和河网“套合差”两个指标进行计算分析。提取水系河流总长52 061.6 km，数字水系河流总长48 395.5 km，相对误差为7.5%。基于ArcGIS图层加载1∶1 000 000矢量化数字水系图作为标准做检验，对比数字化水系和提取水系之间的位置偏移，统计由此偏移产生的两水系之间的细碎多边形面积(见图7)。统计结果显示，整个澜湄流域内的细碎多边形面积约为20 000 km2，河网套合差为2.5%，小于3%，吻合程度较好。由图7a可见，上游地区海拔较高、地势变化明显且无较大湖泊，两水系叠加形成的多边形相对较少，多边形面积约为2 100 km2，对应的流域面积为166 000 km2，河网套河差为1.3%；图7b显示为流域中下游地区，该范围内地势逐渐变缓，水系丰富，在柬埔寨境内有较大湖泊洞里萨湖，两水系之间的多边形较上游明显增多且面积增大，计算得中下游的河网套河差为2.9%，上游和中下游的河网套河差结果均在有效范围内，但中下游地区河网偏移程度较上游大，提取水系精度略差于上游地区。总体而言，本文提取水系结果与实际水系较为吻合，能有效、客观地反映澜湄流域水系整体特征，可为后续研究提供数据支撑。

图7 提取水系与数字化水系河网套河差

3.3 河网分级和子流域划分

由90 m分辨率高程数据提取最优集水面积阈值的澜湄流域水系情况如图6所示，河网分为1～6级，一级为河道最低级别的河源水系，共计1 342条，占总河道数的50%，一级河网发育系数约为2.9，河系不均匀系数约为1.3。

基于Arcgis自动化提取流域出口，共划分小流域3 259个，流域面积普遍在50～900 km2范围。为满足业务应用需要，依据水系上下游汇流关系，建立河道水系与小流域以及各小流域间的拓扑关系；参考流域内地形特征以及湖泊范围，满足水系连续性、分水线完整性、出水口准确性等关键条件；对流域面积在5 km2以下的微小流域或中间汇流区，以汇流关系为主要依据合并到相邻小流域内[31-32]；另外，又以流域内干流上的8个主要水文站点作为流域出口对其进行流域划分；进而完成人工合并小流域的工作。将澜湄流域的子流域划分出3个等级，根据子流域面积由小到大依次定义为一级、二级、三级流域。一级流域共计9个子流域，子流域面积最小的是昌都—旧州区间，约为1.6万km2，最大的是穆达汉—上丁区间，约为24万km2。二级流域共计61个子流域，子流域面积小于5 000 km2的占比25%，在5 000～10 000 km2的占23%，在10 000～20 000 km2的占33%，大于20 000 km2的占19%。三级流域共计328个子流域，子流域面积普遍控制在500～5 000 km2，占比达到90%。其中，小于500 km2的为不适合与相邻合并的中间汇流区，另存在一个大于5 000 km2的子流域，受洞里萨湖水域限制，包含湖域区域，子流域面积约6 700 km2。

3.4 不同子流域尺度流域面雨量预报分析

选取2021年7月24日～25日的一次降水过程预报为例，即为24日8时起报的24 h面雨量，分析不划分子流域的澜湄流域以及划分出一级流域、二级流域的面雨量结果对此次过程在空间分布上的响应。未划分子流域的全流域经算术平均计算面雨量为25 mm，全流域面积大仅以一个数值表达面雨量，在降水分布均匀且降雨量较小时可起到一定参考作用，但在降雨空间分布不均时无法体现降雨中心。图8a中一级流域面雨量计算结果显示，昌都—旧州和穆达汉—上丁面雨量量级为小雨，旧州—允景洪量级为中雨，允景洪—清盛、清盛—琅勃拉邦和琅勃拉邦—万象段量级为大雨，旧州—允景洪和万象—穆达汉量级达到暴雨。图8b显示，二级流域进一步缩小子流域面积，面雨量结果与智能网格预报降雨结果在暴雨中心的空间体现上更为一致，降雨较大范围集中于中下游地区，泰国中东部和老挝中南部交界的位置，有6个子流域面雨量超过60 mm，量级为大暴雨，最大在泰国境内达到143 mm，与一级流域万象—穆达汉段暴雨级面雨量相呼应。

图8 面雨量空间分布

4 结论与讨论

本文分析了澜湄流域集水面积阈值与水系特征变化的数量关系，建立了集水面积阈值与河网密度之间的幂函数，通过曲线割线斜率法客服主观性，确定了90 m分辨率高程下的最佳集水面积阈值为129.1km2。并以1∶1 000 000的矢量化数字水系图做验证依据，总河长相对误差为7.5%，河网套河差为2.5%，均说明提取水系与实际情况吻合程度较高，提取水系有效，进一步划分了澜湄流域1～6级的水系和3个等级的子流域。基于智能网格预报降水，以一次降水过程展现面雨量在不同等级子流域上的空间分布，子流域等级越大、划分越细致对暴雨中心的体现越明显。本文的研究成果，可在水文气象业务上为澜湄流域的面雨量监测、预报和水文模型的构建提供参考，为数字孪生流域建设与流域防洪减灾奠定基础。

虽然本文对DEM进行了修正工作，但是在处理较大湖泊、入海口三角洲等特殊地形上与真实情况方面仍存在差异；研究区域面积大，在高精度高程数据情况下运行速度慢，选取的集水阈值数据量有限。下一步工作可以考虑按照地形、水系类型等因素先将整个流域划分为几个不同的小区域，以求在保证高程精度情况下提高运算效率，进而综合模拟流域水系。