一种无径流资料地区汇流计算方法与流程

本发明涉及一种无径流资料地区汇流计算方法，属于水文预报技术领域，可应用于中小河流无径流资料地区暴雨洪水预报。

背景技术：

水文预报是一项重要的防洪减灾非工程措施。近年来，随着大江大河防灾减灾治理措施的完善和加强，中小河流因降雨引发洪水灾害问题日益突出，全国约有2/3的中小河流达不到规定的防洪标准，常遇洪灾以及诱发的山洪、泥石流等山地灾害造成的人员伤亡约占水灾伤亡总人数的2/3以上。为完善中小河流水文监测站网和预测预报预警体系，提高中小河流防灾减灾能力，2010、2011年国家相继出台了《关于切实加强全国中小河流治理和山洪地质灾害防治若干意见》、《全国中小河流治理和中小水库除险加固、山洪地质灾害防治、易灾地区生态环境综合治理总体规划》指出中小河流水文监测预报预警体系是中小河流治理的重要非工程措施。

无资料流域水文预报方法主要思想一般是选取特定有资料流域作为参证流域，将参证流域预报方法移用于无资料流域。在我国广泛采用以降雨径流相关图为代表的产流方法，结合经验单位线为代表的汇流方法进行洪水预报。一般对于一个特定的区域，在气候、下垫面条件基本一致的条件下，流域产流特性相对稳定，可直接移用降雨径流相关图。但对于汇流计算来说，流域汇流过程具有高度的时空敏感性，要进行汇流方法的移用必须突破时间尺度和空间尺度的变换难题。

现有技术存在如下缺点：

相似流域选取缺乏客观针对性：现有研究通常选用一种或几种地形、地貌特征参数的综合指数(如地形指数ln(α/tanβ))作为属性判别寻找相似流域，而未针对特定的汇流方法所采用的相似准则进行判断。因不同方法在描述研究对象时，均存在一定的假设前提，即抓住影响该过程的事物矛盾主要方面而对次要方面进行一定的简化处理。不对特定方法选用的相似准则进行分析，仅依据某种或几种流域特征参数进行判别不具有客观针对性。

现有技术仅实现了经验单位线时间尺度上的变换，而对空间尺度变换未能突破，导致传统经验单位线方法无法移用于缺资料地区。

技术实现要素：

为弥补现有技术缺陷，本发明提供了一种无径流资料地区汇流计算方法，通过对单位线方法相似准则的分析及空间尺度变换技术的突破，为无径流资料流域洪水过程预报提供一套单位线汇流预报方法。本发明是采用如下技术手段实现发明目的的：一种无径流资料地区汇流计算方法，包括如下步骤：

1)、收集精度在乙等以上的流域汇流单位线，无径流资料流域和有径流资料流域的DEM资料和土地利用资料，以及流域内或流域所在水系区域的断面平均水深资料；通过所述DEM资料获取地表坡度数据和流向数据，通过土地利用资料获取地表糙率数据；

2)、根据曼宁公式计算各流域的空间流速场；

3)、针对各流域，根据流域空间流速场和对应DEM资料的网格长度得到网格单元汇流时间场，再结合流向数据，得到单元格到流域出口的汇流累积时间场；

4)、针对各流域，假定流域内各网格均有一独立运动净雨水滴，并按照步骤4)计算得到的汇流累积时间场到达流域出口，统计流域内不同时刻到达流域出口净雨水滴数目的数量-时间分布曲线，记录为N-T曲线；

5)、将N-T曲线的起始点P₀、峰值点P_m和结束点P_T进行连接，得到各流域的N-T曲线三角形概化过程线；

6)、将待确定汇流单位线的无径流资料流域三角形概化过程线与流域汇流单位线精度在乙等以上的流域的三角形概化过程线进行相似比较，找出相似度较高的数个相似流域作为相似流域集；

7)、依据相似流域间流域面积比的0.5次方作为比例因子，按照提出的转换公式，将相似流域的汇流单位线转换至无径流资料流域，得到若干汇流单位线集，按照水量平衡和峰值大、峰时短的最不利原则进行综合，确定无径流资料地区的汇流单位线，根据汇流单位线进行无径流资料地区的汇流计算。

前述步骤1)中通过坡度分析和水文分析，从DEM资料获取地表坡度数据和流向数据。坡度分析和水文分析过程可采用ArcGis、MapGis等商业软件相应模块实现或利用QGIS、Taudem等开源软件或工具中实现。不同土地利用情况对应的糙率系数如表1所示。

表1不同土地利用条件糙率系数统计表

前述步骤2)中的曼宁公式为：

式①中，i为地表坡度，水深h依据流域内或流域所在水系区域已知点位率定得到的水深同汇水面积之间的函数关系进行计算获得，所述函数关系为：

h＝αA^β ②

式②中A为单元格上游流域汇水面积，α、β为反映流域形态和水系结构的参数，依据流域已有水深资料的点率定得到。

前述步骤3)中网格汇流时间计算公式为：

△τ＝L/v ③

其中L为网格边长，式③用于水流方向与网格边线平行的情况，式④用于水流方向同网格对角线平行的情况。水流流向通过D8法计算，水流方向分为两类，一类与网格边线平行，另一类近同网格对角线平行。D8法参见文献Computer Vision,Graphics,and Image Processing，vol.28,No.3,PAGES:323-344,December 1984,“The extraction of drainage networks from digital elevation data”。

根据径流路径，网格累积汇流时间计算公式为：

式⑤中m为j网格径流路径上的网格数。

前述步骤5)中采用三角形概化N-T曲线原因如下：

N-T曲线反映了单位线的物理本质，在一定程度上代表了单位线形状。控制单位线形状的指标有单位线洪峰流量q_p、洪峰滞时T_p及单位线总历时T。其中(T_p,q_p)即为单位线峰值点P_m，总历时T可以从单位线终点P_T(T,0)得到，起始点P₀通常为原点，属于已知点。因此采用P₀、P_m、P_T三点概化的三角形抓住了单位线过程的主要矛盾，且三角形间相似性比较方法成熟，便于执行。

前述步骤6)中采用的三角形相似定理为：两边成比例且夹角相等的两个三角形相似。将待确定汇流单位线的无径流资料流域三角形概化过程线两边比例及夹角与已有乙等预报方案流域的三角形概化过程线相应两边比例及夹角进行比较，采用距离系数法计算相似距离。

步骤6)中，与待确定汇流单位线的无径流资料流域三角形概化过程线相似距离小于相似度阈值的流域，判定为相似度较高，作为相似流域加入到相似流域集中，两流域间三角形概化过程线的相似距离计算公式如下：

式中，λ₁、λ₂分别为无资料流域和有资料流域两相应邻边的比例，λ₃为无资料流域和有资料流域两邻边夹角的比例。所述相似度阈值优选的取值为0.1。

步骤7)中，单位线变换公式分析推导如下：

单位线是指在特定的流域上，单位时段内时空分布均匀的一次降雨产生的单位净雨量，在流域出口断面形成的地面直接径流过程线。分析上述定义可知，在特定流域上，单位线反映了净雨与流域出口径流的响应关系，而径流总量与净雨量一致，径流过程是净雨经流域调配后在流域出口重新汇集的过程。天热条件下，流域调配的方式即是流域下垫面在重力作用下运移净雨的速度。水文学家从系统论、物理学和概率论所作的揭示表明，流域瞬时单位线实质上就是均匀分布降落在流域上的无穷个净雨水滴汇集至流域出口断面所形成的时间-面积曲线或汇流时间概率密度。由此可得出，单位线是流域汇流速度场的反映，单位线方法的相似因子为流域汇流速度。

单位线空间转换方法分析如下：

在流域尺度上，流域的汇流速度可表示为

式⑦中S为流域面积(量纲为长度平方)，T为汇流时间(量纲为时间)，两者相比与速度量纲一致。

依据前述步骤1)分析，则若流域A与流域B具有相似性，有：

令流域A、B相似因子比例系数为常数C，即：

流域A与流域B具有相似性，则C接近于1，则有：

式中S_A、S_B分别为流域A、B面积，T_A、T_B分别为流域A、B汇流时间。由式⑩可知在流域间汇流速度相似的情景下，单位线汇流时段长与流域面积0.5次方成正比。此外，式⑩创新性的建立了单位线时段因子与空间尺度的转换关系。

在实际计算中，单位线常取等时段间隔进行演算，如Δt＝1h,2h…6h等。令流域A、B划分的时段间隔数相等，均为n，即

式中，Δt_A、Δt_B分别为流域A、B单位线的时段间隔，其余同上述。

另外，由水量平衡关系有

式中q_A,i、q_B,i分别为流域A、B单位线流量，其余同上述。由式⑧可知，在流域间汇流速度相似的情景下，流域间单位线累积流量与流域面积0.5次方也成正比，但比例系数是汇流时长的倒数。式⑧创新性的建立了单位线流量因子与空间尺度的转换关系。式④和式⑧实现了任意相似流域单位线空间尺度转换，在非线性的汇流单位线移植领域取得了突破性的进展。

在实际应用中，认为各时段的汇流速度也较为一致，此时各时段对应的单位线流量比与总流量比例相等，即：

至此，联立上述⑩、两式，即可完成将有资料流域A的单位线转换至无资料的B流域。转换后，流域B的时段间隔Δt_B通常为小数，为满足实际要求，可对转换结果进行插值得到整时段单位线结果。

本发明的有益效果为：

1)、建立了将历史已完成的大量具有较高精度的单位线汇流方案移用于无径流资料地区技术方案，实现了无径流资料地区也可采用传统经验方案进行洪水过程预报的设想；

2)、相似因子选取具有方法针对性：从单位线汇流方法的物理本质出发，提出了将流域汇流速度作为相似因子，且相似因子计算均采用所有可以获取的信息；

3)、相似流域判断科学、简便：依据单位线三要素构建概化三角形，抓住了反映单位线过程的主要矛盾，采用三角形进行相似判断，方法成熟、简便。

附图说明

图1为本发明无径流资料地区汇流计算方法技术框架图。

图2为本发明三角形概化N-T曲线图。

图3为本发明相似度计算图。

图4为本发明单位线综合图。

具体实施方式

下面结合附图对本发明的较佳实施例进行详细阐述，以使本发明的优点和特征能更易于被本领域技术人员理解，从而对本发明的保护范围做出更为清楚明确的界定。

如图1所示，展现了本发明方法的理论分析、技术实现和具体的方法流程。

本发明理论分析如下：

从单位线方法的基本概念出发，结合系统论、物理学和概率论理论对单位线方法的本质进行揭示，分析得出单位线方法相似因子为流域汇流速度。

从系统论角度来看，流域汇流系统可以用线性动力系统微分方程进行描述，既

式中Q(t)为出口断面流量过程，h(t)为相应于出口断面流量过程的净雨过程，a₀、a₁、…、a_m,b₀、b₁、…、b_n均为常数或者仅是时间t的函数。引入微分算子D＝d/dt，有

根据单位线定义若净雨输入h(t)是单位净雨量，则Q(t)就变成单位线过程了，既单位线反映了流域系统的作用，与流域地形地貌具有相关性。

从物理学角度来看，假定净雨是在τ时段内呈空间分布均匀状态注入流域的，且由n个水滴组成，那么均匀分布在流域上的n个水滴都是从τ时刻开始从流域上不同位置向流域出口断面汇集的。由于各水滴到达出口的汇流时间不一样，τ时刻注入的n个水滴必将在不同时刻到达流域出口断面。显然，只有那些流域汇流时间正好等于t-τ的水滴能在t时刻到达出口断面，从而组成了t时刻的流域出口断面流量。如果把满足这个条件的水滴所占据的面积计算出来，并用A(t-τ)表示，那么这个面积与τ时刻降雨强度i(τ)的乘积，就是τ时刻净雨对t时刻流域出口断面流量的贡献，即

对0-t时刻所用输入净雨进行积分有

式中若净雨输入i是单位净雨量，则Q(t)就转变成单位线过程了，既单位线反映了净雨水滴汇集至出口断面的时间-面积变化过程。

从概率论角度来看，流域汇流包括两个事件过程，一是雨滴降落到流域表面，二是雨滴在流域表面向出口汇集。雨滴降落到流域表明可视为一随机事件，雨滴降落的时间τ是随机变量，对净雨强度h(t)过程做归一化处理即为雨滴降落时间τ的分布密度函数f_h(τ)。

不同雨滴到达流域出口的时间T也作为随机变量，对流域出口断面流量过程Q(t)作归一化处理，得到雨滴汇集时间T的分布密度函数f_Q(t)。

另一方面，雨滴到达出口的时间T由雨滴注入流域表面时间τ和水滴汇流时间ω组成，既

若用f_ω(t)表述雨滴流域汇流时间ω的分布密度函数，依据概率论随机函数分布理论有

带入以上两式，并由水量平关系得到

上式中若净雨输入h是单位净雨量，则Q(t)又转变成单位线过程了，既单位线反映了净雨水滴汇集至出口断面的汇流时间分布密度函数。

上述揭示表明，流域瞬时单位线实质上就是均匀分布降落在流域上的无穷个净雨水滴汇集至流域出口断面所形成的时间-面积曲线或汇流时间概率密度，而汇流时间在流域形状即汇流路径确定后，直接受控于汇流速度场。由此可得出，单位线方法的相似因子为流域汇流速度。

技术实现如下：

采用水文相似理论分析方法——量纲分析法，即对物理量的量纲进行考察。通过步骤1)的分析，单位线方法的物理本质是对净雨水滴汇流速度的响应，而速度V的量纲为长度L与时间T的比，即

从流域角度来看，单位线过程是流域面上所有净雨水滴按照一定的汇流速度汇集至出口断面产生的出流过程线。对单个雨滴i来说，其汇流速度V_i是路径长度L_i与时间T_i的比值，但对所有水滴来说，无数个水滴路径的累积叠加，使得研究对象发生了量变到质变的转换，即水滴路径不再是简单的线而是面——流域表面。从量纲角度来看，流域面S的根号与时间T的比与速度V量纲一致，即

依据前述理论分析，则若流域A与流域B具有相似性，有：

令流域A、B相似因子比例系数为常数C，即：

流域A与流域B具有相似性，则C接近于1，则有：

式中v_A、v_B分别为流域A、B整体汇流速度，S_A、S_B分别为流域A、B面积，T_A、T_B分别为流域A、B汇流时间。由式可知在流域间汇流速度相似的情景下，单位线汇流时段长与流域面积0.5次方成正比。此外，式创新性的建立了单位线时段因子与空间尺度的转换关系。

在实际计算中，单位线常取等时段间隔进行演算，如Δt＝1h,2h…6h等。令流域A、B划分的时段间隔数相等，均为n，即

式中，Δt_A、Δt_B分别为流域A、B单位线的时段间隔，其余同上述。

另外，由水量平衡关系有

式中q_A,i、q_B,i分别为流域A、B单位线流量，其余同上述。由式可知，在流域间汇流速度相似的情景下，流域间单位线累积流量与流域面积0.5次方也成正比，但比例系数是汇流时长的倒数。式创新性的建立了单位线流量因子与空间尺度的转换关系。式和式实现了任意相似流域单位线空间尺度转换，在非线性的汇流单位线移植领域取得了突破性的进展。

在实际应用中，认为各时段的汇流速度也较为一致，此时各时段对应的单位线流量比与总流量比例相等，即：

至此，联立上述两式，即可完成将有资料流域A的单位线转换至无资料的B流域。转换后，流域B的时段间隔Δt_B通常为小数，为满足实际要求，可对转换结果进行插值得到整时段单位线结果。

本实施例无径流资料地区汇流计算方法，具体步骤如下：

1)、收集精度在乙等以上的流域汇流单位线，无径流资料流域和有径流资料流域的DEM资料和土地利用资料，以及流域内或流域所在水系区域的断面平均水深资料；通过所述DEM资料获取地表坡度数据、流向数据和汇水面积数据，通过土地利用资料获取地表糙率数据。

通过坡度分析计算得到地表坡度数据，水文分析过程计算得到流向数据，地表糙率数据通过查询土地利用资料获得。地表糙率与土地利用关系表格如下：

2)、根据曼宁公式计算各流域的流域空间流速场；

采用曼宁公式，建立流域单元网格速度计算公式。曼宁公式是明渠道流量和速度的经验公式，水文学上研究表明，可以将其引入流域汇流速度计算。假定流域内每个网格上的水质点相互独立，即每个网格上的水流速度与周边网格上的流速无关，并且不随时间变化，从而建立一个空间上变化而时间上恒定的空间流速场。曼宁公式为：

其中地表坡度i通过DEM获得，地表糙率n可利用土地利用情况获得，水力半径r用水深h替代；所述水深依据流域内或流域所在水系区域已知点位率定得到的水深与汇水面积之间的函数关系进行计算，所述函数关系为：

式中A为单元格上游流域汇水面积，α、β为反映流域形态和水系结构的参数，依据流域已有水深资料的点率定得到。

水深计算方法为现有技术范畴，可参见文献WATER RESOURCES RESEARCH，VIL.34,NO.7,PAGES 1809-1818，JULY 1998,“Energy Dissipation Theories and Optimal Channel Characteristics of River Networks”。

3)依据步骤2)计算得到的流域汇流速度场，计算网格雨滴汇至出口的累积时间场。以网格j为例，其汇流时间计算公式为

其中L_j为网格j边长，式用于水流方向与网格边线平行的情况，式用于水流方向同网格对角线平行的情况。依据流向数据场，以网格j为出发点，依次追踪其达到流域出口所流经的网格单元，如为1、2、…m。则j网格的累积汇流时间τ_j为

4)、针对各流域，假定流域内各网格均有一独立运动净雨水滴，并按照步骤3)计算得到的汇流累积时间场到达流域出口，统计流域内不同时刻到达流域出口净雨水滴数目的数量-时间分布曲线，记录为N-T曲线。

5)、采用三角形对步骤4)计算得到的N-T曲线进行概化。概化选用的控制点分别为起始点P₀、峰值点P_m和结束点P_T，如图2所示。由于N-T曲线反映了单位线的物理本质，在一定程度上代表了单位线形状。采用三角形进行概化一方面可以利用三角形相似定理，便于流域间的相似比较，另一方面抓住了单位线峰值、峰值时间和历时三要素构成的主要矛盾。

6)、依据三角形相似定理——两边成比例且夹角相等的两个三角形相似。将待确定汇流单位线的无径流资料流域三角形概化过程线两边比例及夹角与已有乙等预报方案流域的三角形概化过程线相应两边比例及夹角进行比较，采用距离系数法计算相似距离。如附图3中，λ₁、λ₂分别为无资料流域和有资料流域两相应邻边的比例，λ₃为两邻边夹角的比例，则相似距离为

一般取相似度阈值为小于0.1，原则上寻求的相似流域个数大于等于3为满足个数要求。

7)、依据相似流域间流域面积比的0.5次方作为比例因子，按照公式和将相似流域单位线过程经时空变换转换至无径流资料流域，对得到的若干单位线，采用最不利原则，即洪峰最大、洪峰历时短进行综合，得到无径流资料地区单位线进行汇流预报计算。如附图4中，对上述方法计算得到的4条相似转换单位线集，依据单位线集峰值流量q_m和最早峰现时间t_m进行综合，得到最终的综合单位线。

以上所述，仅是本发明的较佳实施例而已，并非对本发明作任何形式上的限制，在不脱离本发明技术方案范围内,当可利用上述揭示的技术内容做出些许更动或修饰为等同变化的等效实施例,但凡是未脱离本发明技术方案的内容，均仍属于本发明技术方案的范围内。