星-机一体的水力几何形态流量估算方法

王鹏飞，杨胜天，王 娟，潘天力，章四龙，巩同梁

（1. 北京师范大学 水科学研究院，北京 100875；2. 城市水循环与海绵城市技术北京市重点实验室，北京 100875；3. 北京师范大学 地理科学学部，北京 100875；4. 西藏自治区水利厅，西藏 拉萨 850000）

1 研究背景

水文资料是现代社会发展所需的基础数据［1］，但很多地区因自然环境恶劣、经济发展滞后，难以建立传统的实测水文站，造成水文资料匮乏。为了在这些地区获取有效的水文资料，在理论和技术方面应用遥感数据监测河流流量已开展了大量的研究［2-4］，并取得了一系列的成果［5-7］。这些成果可归纳为两类：第一类是利用遥感观测河流中的流量指示物，建立指示物与流量的关系，通过观测这些指示物的变化来估算流量。如Ling等［8］选取长江中的江心洲，利用长序列遥感数据获取江心洲不同时期的出露面积，建立“岛屿面积-径流”的关系曲线监测长江流量变化。此类方法的主要限制是在实际应用中难以找到理想的流量指示物。另一类是利用遥感观测河流水面宽度、水位、河道坡度、波纹等表征河流水力几何形态的对象，并利用水力学方程估算流量。如Huang等［7］在雅鲁藏布江上游段，应用遥感观测河宽与水深，通过改进的曼宁公式计算断面流量。相比第一类方法，因水力几何形态参数广泛存在于河流中且更容易获取，因此第二类方法具有更大的推广应用价值［9-12］。

围绕如何利用水力几何形态估算流量问题，Leopold 和Maddock［13］提出了水力几何形态与流量的关系式，为利用河流水力特性估算流量奠定了算法基础。Leopold发现虽然不同河流有不同的水力几何形态特征，但仍可用河流的水力几何形态参数来表征流量的变化［13］；随后Knighton 深入研究了不同河流水力几何形态差异的原因，认为断面形态变化、边滩类型、江心洲发育等是造成这种差异的主要因素［14-15］；Park统计分析了139条河流的断面资料，给出了水力几何形态参数的取值范围，促进了该方法的推广应用［16］；Ferguson根据河流断面类型，用矩形、三角形、抛物线形断面来表征河道水力几何形态参数，并应用于流量估算，使水力几何形态估算流量的方法更加容易［17］。Rhodes对全球587组河流进行流量估算的结果，也证明用水面宽度来估算流量具有可靠性和实用性［18］。在此基础上，应用遥感技术获取河流水深和水面宽度，基于水力几何特征计算流量的方法得以探索［7，9，15，18］，如Gleason和Smith通过Landsat卫星遥感数据获取河流水面宽，模拟洪水的峰值流量及变化过程［19］。

但当应用水力几何形态理论开展遥感估算流量时，也遇到两方面的困难。首先是通过遥感提取河流水深信息仍然存在较大误差，特别是中小河流，遥感数据的水深误差往往超过了河流正常水深变化的范围；其次，尽管水面宽度比水深的遥感信息提取更容易准确，遥感数据仅需要区分水体与堤岸的边界，即可提取到水面宽度信息，然而由于卫星遥感数据分辨率低的原因，现有研究仅适用在大型河流上［16，20-21］。因此，充分发挥低空遥感数据高精度的优势［22］，以及卫星遥感数据时间序列数据丰富的特点，在水文资料匮乏区探索基于水力几何形态理论的流量估算方法具有一定的实际应用价值［23-25］。为此本文将卫星遥感数据与低空遥感数据相结合，主要开展以下研究：（1）依据水力几何形态理论，通过低空遥感数据建立表征河流水力几何形态的数字模型；（2）利用卫星遥感长序列数据提取不同时期水面宽度，在建立的水力几何形态数字模型上估算流量；（3）用实测数据评价计算误差，分析星-机一体的水力几何形体流量估算方法的可靠性。

2 研究区

本研究以塔里木河的重要支流玉龙喀什河为研究对象，选取该河中上游的同古孜洛克水文断面进行流量估算研究。玉龙喀什河发源于昆仑山南区，河流主要补给水源为冰川融水。河流自高山区至中山区为主要产流区，低山丘陵地区水源补充甚少，出山口至各绿洲灌区为散流段，最终与喀拉喀什河、叶尔羌河汇于塔里木河。同古孜洛克水文断面至河源之间河长340 km，控制上游流域面积15 098.5 km2。玉龙喀什河水系如图1所示。

图1 玉龙喀什河水系分布

3 方法

以水力几何形态为理论基础，利用低空遥感数据和卫星遥感数据按以下3 个步骤进行流量估算（图2）。通过无人机获取的低空遥感数据构建河流的3D模型，并将河段分割，概化出河段水力几何形态关系；在河段水力几何形态关系基础上，通过卫星遥感数据获取流量估算所需参数值；进行流量估算并进行精度验证。

图2 水力几何形态流量估算过程

图3 水力几何形态关系建立过程（单位：m）

3.1 水力几何形态关系的建立

（1）天然河道数字化。使用无人机对研究河段及周围地形进行扫描，获取低空遥感数据，使用Pix4D、ArcGIS等软件完成天然河道的数字化（图3（a））。（2）天然河道的分割。天然河道在发育过程中由于水流的侵蚀、搬运、堆积作用，难以使用统一的关系表述不同位置的断面形状特征。目标河段的分割产生多个较短的河段，在这些较短的河段内可以归纳出表征断面形状关系的水面宽-水深关系曲线（图3（b））。（3）分割河段断面的水面宽-水深关系曲线。提取分割河段的高程数据，获取不同水深对应的水面宽度数据，绘制分割河段断面的水面宽-水深关系曲线（图3（c））。（4）概化河段。在获取分割河段断面水面宽-水深关系曲线的基础上，对目标河流进行整体的概化。计算同一水深不同分割河段水面宽的平均值，最终使用单一的水面宽-水深关系曲线表示研究河段的断面形状特征。对这一关系曲线使用确定的方程进行拟合，实现输入任意水深（水面宽）计算对应水面宽（水深）的结果（图3（d））［26-27］。（5）流量估算。依据水力几何形态理论，建立流量与水面宽度、水深的对应关系，通过监测水面宽度或水深的变化即可进行流量估算的研究。

3.2 平均水面宽计算 现有卫星数据对以线为特征的水面宽监测存在较大的误差，而以面为特征的水面面积监测则较为成熟，因此通过测量面积间接计算河宽是提升数据精度的有效方法［28-29］。使用研究河段水面面积与河段长的比值作为该河段的平均宽度：

式中：w 为概化河段宽度；A为研究河段的水面面积；L为研究河段的长度。

3.3 水体表面提取 水体表面提取的方法是以水体对入射能量的吸收特征为基础［30-31］，其中McFeeters［32］提出的归一化水体指数（Normalized Difference Water Index，NDWI）利用水体对绿光与近红外吸收的差异计算水体范围，是一种使用率较高的方法［33-35］。NDWI对水体的计算方法如下式（2）所示：

式中：Gb为绿光波段的表面反射率；NIR为近红外波段的表面反射率。

当研究河段上空覆盖薄雾时，NDWI的水体表面提取不准确，参考已有成果使用了手动勾绘水面范围的方法［35］。

3.4 流量估算方法 传统的水文水力计算中，糙率系数是一个重要的动态参数，而以往的研究表明，在糙率系数缺失时仍然可以实现流量的精确估算［5，9，36］；同时作为附加比降存在的水面梯度是和水位、河底坡度密切相关的参数，Bjerklie 等［5］提出使用易于测量的水面坡度代替地形坡度的方法。因此在固定的研究河段上，将水力梯度及糙率系数视为常量，避免动态变量的增加，在此基础上对曼宁公式进行整理，合并相关参数以适应不同的数据源和断面形状，实现以水深或水面宽为输入参数进行流量估算。Sichangi 等［9］在水力几何形态理论的基础上，提出一种当河道断面为三角形、梯形、圆弧形时的流量估算方法，如式（3）—式（5）所示：

式中：q为流量；w为水面宽度；h为水深；a，b，c为待标定的综合参数。

3.5 精度评估方法 为评价该方法在流量估算上的合理性与可靠性，选用均方根误差（RMSE）、纳什效率系数（NSE）、相对精度（RA）和平均百分比误差（MPE）作为精度评价方法，计算公式如下：

式中：Qe为估算流量；Qm为水文站记录的实测数据；Qam为实测值的平均值；n为模拟次数。

RMSE被用来描述估算流量与实测值之间的偏差，NSE用以评价模拟结果的优劣，MPE是一种误差整体分析方法，判断估测结果是否在合理误差范围内。上述3 种指数是对流量估算结果的整体评价，作为补充，添加RA对单次估算结果进行精度评价，选用实测值的20%作为许可误差范围［37］。

4 数据

4.1 低空遥感数据 以大疆精灵4无人机（表1）为低空遥感数据获取平台，配合Pix4D软件规划飞行任务，获取研究区的低空遥感数据。低空遥感数据用于获取高精度的河段地形数据，包括了厘米级的高程数据（DSM）和正射影像（DOM）［22］。

表1 Phantom-4-pro无人机基本参数

使用无人机获取目标河段低空遥感数据分为3步：（1）规划航线。根据目标河段所处位置及周围环境设置无人机飞行的高度和飞行范围。（2）实施飞行任务。无人机飞行任务在软件的控制下可以实现完全的自动化。将上一步中设置的航线导入飞控系统中，无人机即可自主完成起飞、巡航、拍摄、降落过程，工作期间获取的数据将储存在机身携带的储存卡或传输到与控制器相连的移动设备中。（3）数据处理与影像拼接。无人机获取的影片经Pix4D mapper进行拼接、点云加密处理，得到同古孜洛克河段的DOM、DSM和3D模型（图4）。

图4 同古孜洛克河段DOM影像、DSM影像与3D数字影像

4.2 卫星遥感数据 哨兵系列卫星中，Sentinel-2提供的数据及处理方法较为成熟，在地物识别及水域分析方面应用广泛［35］，研究中结合NDWI 计算水体表面的需要，选用Sentinel-2 10 m 分辨率的波段3（Gb）和波段8（NIR）数据进行水面提取［38］。

4.3 水文站数据 为验证基于水力几何形态的流量估算结果的可靠性，计算式（3）—式（5）所需的综合参数，整理同古孜洛克水文站2006—2018 年的实测流量和水深数据。将该项资料分割为率定期（2006—2015年）和验证期（2016—2018年）两部分，前者用以率定流量估算所需的综合参数a、b、c，后者用以验证估算流量的可靠性。

5 结果

5.1 分割断面水面宽-水深关系曲线 同古孜洛克水文站所处的440 m 河段均匀的划分出8 个断面，通过ArcGIS的3D Analyst模块在DSM影像上提取分割断面的高程，绘制各个断面的水面宽-水深关系曲线（图5）。

图5 分割断面形状及水面宽-水深关系曲线

5.2 概化河段水力几何形态特征 在分割断面形状及其水面宽-水深关系曲线基础上，合并计算同一水深不同分割河段水面宽的平均值，最终使用水面宽-水深关系平均曲线作为研究河段的断面形状特征（图6）。

图6 概化河段形状特征及拟合曲线

研究中使用了数值分析软件1stopt计算水面宽-水深关系曲线的拟合方程。该软件提供的通用全局优化算法（Universal Global Optimization，UGO），依据其数据库中丰富的数学模型进行迭代计算，求出最优的拟合结果。结果显示，依据麦考特迭代法（Levenberg-Marquardt）的全局优化算法对数据点的拟合效果良好，得出拟合方程及残差分析（表2）。

由于目前卫星遥感数据对水深的测量仍然存在较大误差，这种误差往往超过中小河流的正常水深变化范围［39］，而水面宽度的测量较为简单，利用h=f（w）方程即可在卫星遥感数据获取的水面宽的基础上计算水深，弥补卫星遥感数据对水深测量的不足；获取的水文资料中仅记录了水深与流量数据，没有记录水面宽度的变化。而基于水力几何形态的流量估算方法（式（3）和式（5））在率定期需要水面宽参与综合参数a、b、c的率定，因此利用w=g（h）方程，由实测的水深数据计算水面宽，完成参数率定的需求。

表2 概化河段形状特征方程

5.3 水面面积、水面宽及水深 在卫星遥感数据的基础上，通过NDWI计算研究区水面范围（图7）。其中2017/4/10和2018/6/4两日由于同古孜洛克河道上空被薄雾覆盖，NDWI计算结果与实际水面范围相差较大，使用了手动勾绘水面范围的方法。

在提取同古孜洛克河段水面范围基础上，计算了水面面积、水面宽度和水深（表3）。在提取水面范围时，由于低空遥感数据边界与卫星遥感数据边界并未重合，以卫星遥感数据为准，上下游边界各扩展一个像素宽（10 m），此处河长实际为460 m。

5.4 流量估算 已取得验证期流量估算所需的水面宽及水深数据，综合性参数a、b、c 由率定期实测水文数据计算。验证期实测数据及式（3）、式（4）、式（5）的流量估算结果如图8 所示，其中图8（a）、图8（b）和图8（c）分别展示了式（3）、式（4）和式（5）的估算值与实测值的差异；图8（d）则着重于式（3）、式（4）和式（5）之间的差异分析。

图8 流量估算结果与同时期水文站实测值

6 讨论

6.1 精度评价 考虑到验证期实测流量的变化区间过大（最小值为7.32 m3/s，最大值为900 m3/s），精度评价中采用分级的方法对估算结果进行更为详细的模拟效果分析（表4）。按照验证期实测流量的大小将估算结果分为两个等级：≥200 m3/s 和＜200 m3/s。前者包括2016/6/24、2016/8/3、2017/7/14、2017/7/19和2017/8/28的估算案例，代表了较高的流量等级；后者是剩余的18次估算案例，代表了低水平的流量等级。

表4 分级精度评价结果

由分级的精度评价结果可知：（1）当流量等级＜200 m3/s，式（3）的各项评价结果均优于式（4）和式（5）。注重于总体评价效果的RMSE、NSE、MPE 均表明式（3）在低水平流量等级下具有较好的模拟结果。特别是RMSE的结果为4.65 m3/s，在水文资料匮乏区中小型河流普遍较低水平的流量等级的条件下具有很强的应用价值；关注于单次估算精度的RA，其平均合格率为94.44%，18个估算案例中有17个在允许误差范围内，同样表明了基于水面宽的水力几何形态流量估算方法在该地区具有较好的适用性。（2）当流量等级≥200 m3/s，耦合了水深与水面宽的式（5）展示了更好的估算结果。RMSE为均值的17.48%；同样关注于总体评价的NSE、MPE也表明式（5）的模拟结果较为可靠；对于关注单次估算精度的RA而言，使用式（5）进行流量估算的结果中，80%的流量估算结果在许可误差的范围内。总体而言，在引入水深和水面宽作为水力几何形态估算流量的参数时，式（5）在同古孜洛克河段高流量等级时具有较好的计算结果。

6.2 断面几何形状特征及水下地形获取 Leopold和Maddock［13］在水力几何形态理论上提出水面宽与流量的关系至今，众多研究在建立“水面宽-流量”关系时都将研究断面选在具有三角形、梯形、抛物线型等具有明显宽深关系的河段［5，7，9］，本研究所提出的方法也同样受制于河段断面类型的限制，难以应用在矩形断面或宽深关系不明显的梯形断面。然而水文资料匮乏区河流的中、上游由于水动力及两岸岩体的性质，河道普遍是三角形或梯形断面，有利于水力几何形态理论流量估算方法在该地区的推广。

研究中所选的同古孜洛克水文断面具有明显的干枯两季，在枯季时河道存水极少，河床暴露，便于获取河段完整的低空遥感数据。但是对于有水淹没的其他河道，提出如下可能的解决方案：（1）实测数据补充缺失的水下地形数据。水深较浅时，可以通过人工测量的方法将实测的水下断面、低空遥感数据提取的水上断面相结合，补充出完整的断面形状（图9（a））。（2）拟合水下地形。在高山峡谷、水流湍急等不具实测条件的地区，对水下地形进行拟合是优先推荐的方法。根据河流所处的地理位置及水上地形信息，使用几何形状拟合水下断面（图9（b））。除抛物线拟合之外，三角形、梯形、矩形、圆弧形亦为常用拟合形状，根据河流所处位置及实际情况选择［40］。（3）估算流量与实测流量的结合。在河流水量较少的时期，只对水面以上的断面进行数字化并记录当时流量的实测值Q0，Q0将作为一个基础值参与到之后流量估算中（图9（c））。该方法较好地解决了水量超出Q0时的情况，但当水量小于Q0时，重新建立断面关系并校正Q0即可。

图9 水下断面获取方法

6.3 分割河段的数量 当目标河流被分割的越短时，单一断面所表示的河槽特征越为准确。研究中所使用的低空遥感数据单像素分辨率为3.4 cm，可以每像素宽提取一次断面形状，刻画一次河槽的形状特征。但过于精细的划分任务限制了该研究的推广和应用价值，同时作为一种探索性尝试，研究中先通过较为稀疏的划分方法进行流量估算研究，论证该方法的可行性。

天然河流作为自然界的一套复杂系统，在岸边约束条件不一致、水动力复杂多变的情况下，使用单一断面堆叠出河段的整体形状特征是不准确的。因此本研究在保留分割河段物理形态的基础上，简化非必要信息，在低空遥感数据获取的440 m河段中划分8个断面进行流量估算研究。流量估算及精度评价结果说明这样的划分方式是有价值的，但其他形式的划分是否合理仍然需要更多的研究证明。

6.4 复杂河流背景下的适用性 同古孜洛克水文断面是三角形或宽深比较大的梯形断面的代表，不同时期水面宽度具有指示流量变化的能力。对于更为复杂的河流，当水面宽具有指示流量变化特征的功能时，本研究所提供的方法仍然具有价值。典型河道过流断面如图10所示，其中图10（a）为常见的山区河流断面形状，具有明显的三角形特征，可以使用水面宽表征流量的变化过程；图10（b）为常见的山区或平原河流断面形状，表现为两侧边坡的不同，水面宽度与流量之间呈现正相关关系；图10（c）和图10（d）为平原区河流断面形状代表，其中图10（c）具有使用水面宽度表征流量的能力，图10（d）近似于矩形的断面则不具有此能力；图10（e）为辫状河流或河道有不可忽略突起时的形状，水面宽同样具有表征流量变化的能力［41］。

图10 典型河道过流断面

本研究所提出的方法是在低空遥感数据的基础上，建立水面宽-水深关系曲线，在卫星数据的支持下进行流量估算，实现时间序列的扩展。在可以使用水面宽表征流量变化的复杂河段，该方法仍然具有适用性；但对于矩形河道或宽深比不明显的梯形河道，水面宽度随着流量的变化难以被监测，该方法的使用受到限制。

7 结论

为解决水文资料匮乏区中小河流流量监测难题，本文建立了星-机一体的水力几何形态流量估算方法，具体为：以大疆精灵4无人机作为低空遥感高精度数据采集平台，建立河段的水力几何形态数字模型；在数字模型基础上，分割河段绘制断面形状，建立水面宽-水深关系曲线；应用卫星遥感数据计算河段水面宽度，并由水面宽-水深关系曲线计算不同时期的水深；以水面宽和水深作为水力几何形态理论的输入参数进行流量估算。验证结果表明，当流量小于200 m3/s时，以水面宽作为关键参数的水力几何形态流量估算效果较好，平均合格率为94.44%；当流量大于200 m3/s 时，以水深和水面宽共同作为流量估算参数时具有更好的模拟效果，平均合格率为80%。

研究结果表明该方法立足于水力几何形态理论基础，发挥低空遥感数据高精度优势和卫星遥感数据时间序列长的特点，为水文资料匮乏区提供一种可以连续进行流量监测的方法。同时，研究结果也表明该方法在矩形断面、宽深比较小的梯型断面以及其他更为复杂的河流断面时仍然有一定的局限性，需要进一步的评估与验证。