基于多源遥感数据的黄河数字流域模型改进

薛 源，覃 超，吴保生，张 玍，李 丹，傅旭东

(1.清华大学 水圈科学与水利工程全国重点实验室，北京 100084；2.清华大学 水利部水圈科学重点实验室，北京 100084；3.清华大学 水利水电工程系，北京 100084；4.中国煤炭科工集团有限公司 煤炭科学研究总院有限公司 应急科学研究院，北京 100013)

1 研究背景

我国地形复杂，山区约占国土面积的89%，发育着众多的山区河流。但山区河流往往因地势险要，无法通过现场量测获得流量等基础资料，是典型的水文缺资料地区，制约了水资源调查、河流水沙通量等重要研究的开展。为缓解这一现状，学者们多通过卫星遥感观测、经验估算或反演山区河流的径流过程[1]，但所得数据往往针对单一断面或河段，丰度有限、精度较差，方法也难以推广[2-3]。因此，利用水文模型能够模拟流域径流过程的优势，结合多源遥感数据大范围提取的河流关键参数，模拟山区等缺资料地区河流的径流过程是当今的研究热点[4-5]。

分布式水文模型，能较充分考虑水沙运动的物理机制及下垫面的空间差异，反映了接近真实的水文过程[6]。其中黄河数字流域模型(Digital Yellow River Integrated Model，DYRIM)，以多源遥感及地理信息系统等技术模拟全流域的降雨-径流过程[7]。模型依托DEM数据，以流域分级理论为依据，将全流域分为4个尺度：坡面-沟道单元、单一子流域、多个子流域的组合和全流域，通过基于坡面-沟道单元的产流、汇流模型，在子流域、流区和全流域水系模拟水沙过程[7]。近年来，DYRIM被诸多学者改进并应用，如李铁键、刘家宏等[8-9]利用DYRIM模拟了黄河流域的径流过程，验证了DYRIM模拟不同下垫面条件、多种类型河流水沙过程可靠性；张昂等[10]在考虑降雨、河宽等参数空间分布的基础上改进了DYRIM，模拟了黄河源区1956—2014年的径流过程，进一步提高了缺资料地区径流模拟的精度，文献[11]也通过考虑植被、淤地坝的影响改进DYRIM。由于推广性较好，DYRIM还被用于模拟和预测珠江流域、闽江流域、南方山区河流、长江上游等区域的径流过程[12-14]。

尽管DYRIM能模拟不同尺度、类型流域较真实的径流过程，但仍有较大提升空间。如DYRIM将河道断面简化为“V”形，虽然能较好反映低级别河流的断面特征[15]，但不能用以模拟高级别河流[16]，即使通过模型率定、调参等方式降低因断面形态失真带来的误差[12，17]，仍会使模拟结果不能反映真实的下垫面条件、水流情况，因此对土壤特性、水力要素的模拟偏离实际。可见，虽然DYRIM的模拟效果整体较好，但仍需基于河流几何信息，获取更真实的断面形貌沿程变化，改进模型、提升模拟效果。

本文以位于黄土高原的黄河一级支流皇甫川流域为研究区域，以多源遥感数据为主要数据源，提取了河流表面信息及各级河流的断面形态，与DEM河网提取工具得到的线状河流水系组合，得到了流域河道的边界条件。基于各级河流的断面形态，建立了DYRIM断面判别模块，改进了其汇流模型，提升了径流过程的模拟效果，提高了洪水过程与水力要素的模拟精度。本研究可为探究基于多源遥感数据获得较完整的河道边界条件、提升水文模拟效果及河流地貌演化等研究提供参考。

2 基于断面判别模块的模型改进

DYRIM是以DEM河网为核心，基于物理机理、分布式参数和连续计算时段，能较准确描述黄河流域下垫面环境，较精确模拟黄河径流过程的分布式水文模型[7]。DYRIM模拟流域径流模拟时，涉及其产流、汇流模型，本文主要针对汇流模型进行改进(图1)。

图1 DYRIM产流及改进的汇流模型框架

2.1 产流、汇流模型及数据输入坡面产流是流域径流过程的起点，为了适应黄河中游超渗产流的特点，DYRIM产流模型采用土水势概念模拟地表入渗过程的变化，获得较精确的产流过程[8]。其后，以每一个坡面及与其相邻沟道组成DYRIM的最小汇流单元——坡面-沟道单元，并假定坡面-沟道单元产流直接作用于单元出口，按自上游至下游的拓扑关系流入下一级河网中完成水流的汇集与演进[16]。支流入汇直接由各级径流叠加，水流演进过程通过断面判别模块调用扩散波、运动波法计算。基于DYRIM产流、汇流模型的径流模拟需输入流域降雨、地形几何参数、下垫面参数，参数率定是兼顾了模型实际意义的部分下垫面参数[10]。

2.2 断面判别模块山区河流多不具备实测断面资料，仅能基于中、低分辨率的DEM提取有限几何参数，使得改进前的汇流模型对水流运动的模拟主要基于简化后对称的“V”型断面，忽略了断面随河流级别由窄“V”型至宽“U”型的沿程变化。虽然原模型对尺度本就较小的低级别河流仍有较好的模拟效果[8-9，16]，但为最大限度保留原模型对低级别河流的模拟效果，并使其能考虑高级别河流的河道特征，建立了DYRIM断面判别模块，自动考虑各级河流的断面形态。

断面判别模块通过综合判别河流级别及断面形态，调用扩散波法[8]、运动波法[18-19]模拟不同级别河流径流(图1)。与以往DYRIM基于“V”型断面的扩散波法模拟水流演进不同，本研究建立的断面判别模块则主要基于概化断面的边界条件和适用于一般梯形断面的优化的运动波法[19]，自动判别断面形态并调用对应的演进方程。

3 研究区域及数据处理

3.1 研究区域概况皇甫川流域位于黄土高原北麓(110°20′E—111°15′E，39°12′N—39°59′N)，面积为3246 km2，于陕西省榆林市府谷县汇入黄河[20]。流域地处半干旱高原地区，砒砂岩大面积裸露，沙地面积大，土壤类型以栗钙土、黄绵土和风沙土为主[17]。流域内多年平均降雨量为365 mm且多以短时强降雨为主，年均蒸发量达1500 mm，2010—2018年流域内主要分布有草地、农作物、灌木、沙地植被等[20]。流域内设皇甫、沙圪堵水文站和12个雨量站，分布见图2(a)，2010—2018年皇甫川流域10 m和30 m分辨率的土地利用叠加数据见图2(b)。

图2 皇甫川流域基本情况

3.2 研究数据及预处理为满足改进后模型需要(表1)，本研究采用了多种类型的数据，包括国产高分1号卫星(GF-1)、国产资源3号卫星(ZY-3)影像，Jason1/2号卫星影像、Landsat系列卫星影像、MODIS产品、SRTM 90 m DEM数据、高分辨率UAV DEM产品、CMPA国产卫星融合降雨数据等多源遥感数据。水文年鉴数据涵盖2010—2016年皇甫川流域皇甫站和沙圪堵站实测资料，以及用于概化河流断面的2010—2015年黄河中游紧邻于皇甫川、分布于9个流域的32个水文站的实测断面(图2(a))。

表1 水文站实测及高分辨率DEM概化断面几何参数

3.2.1 平滩流量及枯水流量确定 河流断面几何信息提取及形态概化的关键是断面边坡坡度，该坡度可由最低、最高/平滩水位对应的河宽确定。这两类河宽资料往往需通过实测获取[20]，十分匮乏。为了能通过遥感影像获取缺资料地区这两类河宽数据，需确定平滩及枯水流量日期。采用皇甫、沙圪堵两站1990—2016年的汛期径流及场次洪水资料，并验证资料一致性[21]。由于皇甫站断面不具备滩槽结构，难以直接确定平滩流量及日期，因此依据缺资料流域平滩流量估算方法及流域干支流关系[22]，以沙圪堵站断面的平滩流量、日期代表皇甫川流域。确定水位为100 m左右的洪水流量为沙圪堵站平滩流量，流量为113 m3/s，重现期约2年。最枯流量的确定，为保证一致性仍以沙圪堵站为准，为0.1 m3/s。

3.2.2 卫星遥感数据及预处理

(1)光学卫星影像获取及预处理：为通过遥感影像提取流域平滩、枯水河宽，选用能在山区更完整提取细小河流的高分辨率光学卫星影像，故联合使用GF-1和ZY-3影像[23]。根据平滩日期，选择了7景1A级GF-1影像，成像时间为分别为2013年7月30日、8月23日，2016年7月26日、9月1日。选择了3景1A级ZY-3影像，成像时间分别为2013年7月30日和2014年9月3日。热红外数据在夏季能有效区分水体及其他地物[24]，因此采集了同期的6景Landsat-8 TIRS产品。在枯水日期下，GF-1影像足够覆盖，直接使用14景1A级影像，成像时间为2013年11月20日、11月28日。

(2)测高卫星数据：基于星载雷达观测的水位将用于提取缺资料地区河流断面形态。通过观测平滩、枯水河宽及其对应的水位信息，即可获取概化河道断面形态[25]。以Jason-1/2卫星产品为数据源，依据两类流量日期提取对应水位。经查询2009年2月—2013年6月，Jason-1卫星的第103号轨道位于流域上空，但受限于时空分辨率，选取与平滩流量相近的2012年8月6日的Jason-1卫星数据，与枯水流量相近的2017年4月17日的Jason-2卫星数据提取相应水位。

3.2.3 土地利用数据及采样 通过本团队开发的河宽自动提取算法ARWE(Automated river widths extraction method)准确提取河流中心线及河宽，该算法因基于随机森林及神经网络算法，故需依赖准确的机器学习样本及合适的分类指数[26]。为确定研究区河流周围地物类型，基于土地利用类型数据[27-28]分析了2010—2018年皇甫川流域地物类型(图2(b))，结果显示2010—2018年流域内土地利用类型变化不大，确定了紧邻河流且易混淆的4类地物，遵循随机原则[29]完成了250余万学习样本的采集，样本像元总面积约占流域的3%，其中河流、农田耕地、山体林地、草场、城市的样本比例为4∶2∶2∶2∶1。2017年研究团队在研究区域内进行了2次野外考察(图2(b))，考察中记录下80个样本点的位置，测量了平滩河宽。因研究区域地物类型变化不大，故这80个实测样本可作为验证样本。

3.2.4 DEM数据及河网提取 DEM河网作为DRYIM的结构核心以及河道边界条件的重要组成，其精度至关重要。低级别河流多因河宽极细，难以通过卫星产品提取[23]，因此其断面概化、河宽提取均依赖高精度DEM数据辅助。数据分为以下几类：

(1)高分辨率DEM数据：团队于2016年使用点位精度6 mm、距离精度4 mm、扫描密度1 cm的徕卡三维激光扫描仪(Leica scan station 2)，经过配准、Cyclone点云处理等流程[30]生成了岔巴沟流域0.1 m DEM数据。于2006年、2017年联合使用机载激光扫仪与实时动态差分定位技术(Real time kinematic)生成了皇甫川流域上游1 m分辨率的DEM和桥沟流域、油房渠流域0.3 m的DEM数据。

(2)基础DEM数据：用于提取黄河中游、皇甫川流域河网的DEM基础数据为2015年公布的4.1版SRTM 90 m DEM，该数据经过NASA、USGS人工处理，高程精度优于15 m，局部可达3 m[31]。在纸坊沟流域使用由陕西省测绘局提供的5 m DEM数据。

基于SRTM 90 m DEM v4.1，以清华大学开发的河网提取工具DNET(Drainage Network Extraction Tool)提取黄河中游线状河流水系。经验证，DNET能在提取流域水系时最大程度的消除平行错误，配合SRTM 90 m DEM v4.1使用提取结果更准[23，32]。使用DNET时流域最小汇流面积经预实验设定为0.1 km2。

3.2.5 CMPA融合降雨数据 DYRIM需降雨数据驱动[7]，以缺资料地区有限的雨量站插值估算流域降水进而模拟径流过程往往使模拟失真[10]。为解决这一问题，本研究采用了分布式卫星降雨融合数据CMPA(China Merged Precipitation Analysis)，通过与皇甫川流域2011—2015年的实测资料比对，验证了CMPA数据精度及可用性[33]。

4 河流断面形态的概化

以一般梯形等概化河流断面形态已有一定的研究基础[25]，因此通过融合有限实测资料和多源遥感数据，建立了研究区各级河流断面概化模型：(1)绘制2010—2015年皇甫川流域及与其紧邻流域共34个水文站的河流实测断面，按梯形断面概化并验证；(2)将实测断面对应的最低河流级别作为断面概化模型中高、低级别河流的分界；(3)基于研究区域及紧邻流域的高分辨率DEM，提取了65个河流断面，按三角形(底宽为0的梯形)概化低级别河流断面并验证；(4)基于多源遥感资料，提取了皇甫川流域的3个梯形断面。基于上述流程，完成了对皇甫川流域各级河流断面的概化。

4.1 基于实测数据概化断面形态在黄河中游基于SRTM 90 m DEM数据共获得8级DNET河流，依此将断面分级。由高分辨率DEM所提断面知1—3级低级别河流断面近似“V”型；水文站实测断面对应河流级别为4—8级，近似呈梯形。

统计实测断面，综合比对断面所处级别后，按梯形(三角形为底宽为0的梯形)概化[25，34]。改进后的DYRIM模拟径流过程时，将按河流级别调用概化断面。概化后的断面几何参数见表1。

4.2 基于遥感观测提取断面形态现有研究对缺资料河流断面的提取多基于河相关系[36]或由水槽实验推导出的经验公式[36]，研究区域多为大江大河的下游，对山区河流、低级别河流的研究相对较少。皇甫川流域暂无相关研究且仅有2个实测断面，为更准确地概化各级断面，除参考实测资料外，还尝试基于遥感观测的河宽、水位提取河流断面。

4.2.1 多源遥感观测的河道断面提取原理 遥感观测断面形态，需提取断面关键几何参数，再依据图3所示方法获取。由于研究区域内的低级别河流河道极窄，即便通过高分辨率卫星也难以提取[37]，因此本方法适用于级别相对较高河流。

图3 多源遥感观测的河流断面示意

基于多源遥感观测提取河流断面形态的步骤如下：

(1)基于PRF-ANN(Parallel random forest and artificial neural network algorithm)河流表面信息提取方法[26]及ARWE河宽自动提取算法，提取流域高分辨率河流表面信息、河流中心线及平滩、最枯河宽。

(2)基于Birkett算法[38]提取经河流表面信息去噪的研究区域Jason-1/2卫星的平滩及枯水水位。

(3)确定Jason卫星轨道与步骤(1)提取的河流表面信息的交线，交线的端点坐标为A、B且认为A、B处水位相等。读取该处水位h(图3)。

(4)按水位递减顺序，从左岸至右岸依次连接端点得图3中最低水位以上断面。

4.2.2 河宽自动提取 基于PRF-ANN及ARWE法提取河宽的流程为：(1)基于PRF-ANN法提取平滩、最枯流量下的河流表面信息。(2)基于ARWE算法，提取对应河流表面信息的河流中心线及水面边界并校正。(3)ARWE基于校正后的河流中心线及水面边界，自动提取河宽。

4.2.3 平滩河宽及枯水河宽提取结果 基于PRF-ANN法，通过预实验确定的流域丰、枯水期能较好区分水体及其他地物的指数，依据这些指数提取流域内水体表面信息(包含河流、湖泊、水库等水体的初步产品)。如丰水期(6—9月)选用了三价铁离子指数(Fe3+)、归一化水体指数(NDWI)、增强水体指数(EVI)、色彩渲染指数550(CRI550)4种指数、2种Landsat-8 TIRS热红外数据和4种灰度共生矩阵(GLCM)纹理指数(绿光、近红外波段均值，红光、近红外波段熵)提取平滩流量下的水体表面信息。枯水期使用形状指数(IF)、氧化铁指数(IO)、显色指数(CI)、Fe3+及蓝光、绿光、红光波段的4种纹理指数提取最枯流量下的水体表面信息。

水体表面信息提取完成后，需基于由流域DNET河流水系生产的缓冲区(半径300～1000 m)去除湖泊、水库、水塘等非河流水体和其他地物，再通过形态学算法、适当人工编辑后，得到精确的河流表面信息[20]。基于精确的河流表面信息，ARWE算法结合DYRIM的子流域划分结果共提取了1360个平滩、枯水河宽数据对(图4(a))。提取的平滩河宽通过80组野外实测数据及15组经高分辨率DEM断面算得的平滩河宽检验，枯水期最枯水位河宽则通过水文站及Google Earth历史影像检验。平滩河宽提取的R2、RMSE、平均偏差MBE分别为0.99、1.49、0.06，最小河宽对应DNET的2级河流；枯水河宽提取的R2、RMSE、MBE分别为0.91、1.61、0.08，最小河宽对应DNET的3级河流。提取所得流域枯、平滩河宽之比为0.31，实测值为0.25，整体结果较准，提取结果的误差分析示意见图4(b)。平滩、枯水河宽数据对将用以提取相应水位。

图4 河宽提取结果及检验

4.2.4 平滩及枯水水位的自动提取 提取Jason系列卫星103号轨道与皇甫川流域河流交线，据交线与河流边界交点确定水位提取点，假定每时期、每断面左右交点水位相等，概化河流断面。受限于Jason卫星的测量幅宽和过境时间等，在2012年8月6日、2017年4月17日仅从皇甫川流域自动提取到3个断面，如图4(a)，断面关键参数见表2。

表2 遥感观测的河道断面

4.3 皇甫川流域断面概化结果验证综合实测断面概化结果(表1)及多源遥感断面概化结果(表2)结果，将皇甫川流域5级河道的左、右岸坡度调整为0.44、0.35，6级河道调整为0.45、0.26。概化后的各级断面与皇甫川流域实测断面对比见图5，验证断面区域见图4。

图5 皇甫川流域断面概化结果示意

以两岸坡度平均偏差定量评价，验证得1—3级低级别河道概化的坡度平均偏差+4.6%，4—7级河道概化边坡坡度平均偏差-3.8%，故概化断面在一定程度上能较好反映断面真实形态。图5(h)(i)为全流域概化模型示意。将沙圪堵、皇甫站按本模型概化后(图5(i))可知除皇甫站因设在乡镇，左岸受人为活动影响使坡角接近90°外，概化断面与实测断面边坡较吻合，可供模型使用。概化模型由皇甫川及与其紧邻流域建立，概化断面的相关参数可供黄河中游上段其他无资料流域参考。

因河流表面信息由DNET水系辅助提取，包含了DNET拓扑信息在内的所有参数，概化断面与河流表面信息的组合可提供流域河道连续的边界条件。

5 改进的数字流域模型模拟结果

5.1 研究区域参数率定基于皇甫川流域河道边界条件，建立DYRIM运行所需数据库。模型率定基于CMPA降雨数据、皇甫及沙圪堵站2010—2015年汛期逐日流量，主要评价指标为径流模拟纳什效率系数(Nash-sutcliffe efficiency coefficient)[39]及《水文情报预报规范》[40]中的洪峰流量相对误差、峰现时间误差。根据双层参数率定法[41]，模型率定结果见表3。

表3 DYRIM参数率定结果

以流域出口的皇甫站为例，率定期径流过程模拟结果如图6。改进后模型模拟2010—2015年汛期日径流的平均NSE为0.83，从Ⅰ—Ⅲ区可得改进后模型模拟日径流过程趋势更平缓、结果更准。

图6 率定期夏汛期径流过程结果及对比

5.2 模拟结果及验证以率定后参数模拟2016年7月12日—20日、7月25日—27日、8月13日—9月5日的径流过程及2016年洪水过程，仍以流域出口的皇甫站为例，模拟结果如图7(a)(b)所示，水力要素模拟较改进前有较大提升。径流模拟NSE达0.84，且模型在Ⅰ—Ⅲ区径流过程更接近实测过程(图7(a))，涨、落水趋势更平稳。进一步分析2010—2016年模拟所得断面水深h及平均流速v，与实测数据验证率定期和模拟期所得水深h的平均误差为9%，平均测小0.04 m、R2为0.86、RMSE为0.10；流速v模拟的平均误差为7%，平均测小0.08 m/s、R2为0.80、RMSE为0.31，水力要素模拟结果均较准确。

图7 2016年夏汛期径流过程及水力要素模拟

对洪水过程的模拟，以模拟难度较大的多峰洪水过程展示效果。模拟时间步长为匹配水文站洪水要素表的6 min。模拟洪峰流量相对误差分别为15.1%、14.9%，均处在《水文情报预报规范》规定的20%范围内，模拟效果见图8，改进后的模型能更真实模拟主洪峰过程中的子过程，且峰现时间与实际接近。

图8 典型洪水过程模拟效果及对比

模拟2010—2016年发生的共计21场洪水时，率定期18场洪水峰现时间较真实值平均提前约1.49 h；模拟期则为1.17 h，均处于±3 h范围内。模拟的峰现时间整体提前可能源于两方面：(1)卫星降雨融合产品在皇甫川流域仍存在一定的偏差；(2)建立数字流域模型网格时密度仍然稀疏，导致汇流时间缩短。

6 结论与展望

本文基于PRF-ANN河流表面信息提取方法及ARWE河宽提取方法，提取了皇甫川流域1360对平滩、最枯河宽数据对。结合有限的水文站实测资料及多源遥感数据，建立了皇甫川流域各级河流的断面概化模型，得到了流域各级较完整的河道边界条件。根据其包含的断面几何信息，建立了DYRIM断面判别模块，改进了DYRIM汇流模型。以2010—2015年为率定期率定DYRIM参数，模拟了2016年皇甫川流域的汛期日径流过程、水力要素及年内的3场洪水过程。改进后的模型模拟效果提升明显，获得了更准确的水力要素。研究成果为水文模型改进提供了重要基础数据和参考，为模拟山区等缺资料地区河流的径流过程及洪水预报提供技术支持。取得的主要结论如下：

(1)建立了基于多源遥感数据的河道边界条件提取方法。该方法以传统线状DEM河流水系为基础，通过PRF-ANN、ARWE法自动提取流域内的河宽数据，结合Jason系列卫星提取的水位数据获得了河流断面。融合3个多源遥感提取断面、34个水文站实测断面及65个高分辨率DEM断面，以一般梯形为基础建立了皇甫川流域1—7级河流的断面形态概化模型。经实测资料验证，最枯河宽提取的R2和RMSE分别为0.91、1.61；断面概化模型的左右岸坡度与实测坡度的相对误差分别为7.2%、4.5%。

(2)建立数字流域模型断面判别模块，该模块可根据河流级别自动识别河道断面形态，并调用对应的水流演进方程综合模拟流域内各级河流的径流过程。对2010—2016年日径流过程模拟及验证的NSE均优于0.8；流速模拟的相对误差、R2、RMSE分别为7%、0.8、0.31；水深模拟的相对误差、R2、RMSE分别为9%、0.86、0.1。

(3)对2016年3场典型多峰洪水的模拟结果表明，改进后的模型模拟洪峰流量的相对误差为14.9%，峰现时间误差约为1.1 h。整个率定、模拟期，21场洪水的洪峰流量相对误差约15%，峰现时间提前约1.4 h。

虽然本研究改进了DYRIM的汇流模型，提升了径流、洪水过程及水力要素的模拟效果，但未能联合不易受云雨干扰的雷达卫星，暂未实现对河流的全天候观测、提取。对断面的概化尚基于较单一的一般梯形，暂未考虑更复杂的复式断面及最枯水位以下形态，较适用于河流的中上游段，对最枯水位下仍存在较大水深的河道断面适用性有待验证。对DYRIM的改进暂引入了较简单的运动波法，还应继续尝试建立更适合体现山区河流洪水及其输沙特征的计算方法。