CN107590354A - 基于稳定水流场的流域地表径流汇流模拟方法 - Google Patents

Abstract

本发明提供了一种基于稳定水流场的流域地表径流汇流模拟方法及实时计算装置，该方法包括：S1、建立栅格模式下的N‑S方程；S2、在径流汇流模型中进行水流混合模拟；S3、利用所述S2得到的径流汇流模型对流域内水流场进行迭代，从而获得稳定水流场。该方法能为流域暴雨洪涝灾害的快速模拟和评估提供辅助决策支持，具有较高的应用价值。

Description

基于稳定水流场的流域地表径流汇流模拟方法

技术领域

本发明涉及水动力学、计算机仿真和GPU通用计算领域。具体涉及流域径流汇流模拟系列关键技术，包括稳定水流场的建立和地表径流汇流的GPU并行计算两个方面。

背景技术

我国中西部地区以山地地貌居多，地处亚热带湿润气候，夏季湿润多雨。在这种地貌以及气候条件下，每年因强降雨诱发的山洪地质灾害频繁发生，对当地的生命和财产安全造成极大的威胁，因而通过对流域地表径流汇流实时模拟，以完成暴雨洪涝灾害的快速模拟和评估具有十分重要的意义。

径流汇流模型是径流汇流模拟的基础。现有技术可大致分为集总式模型和分布式与半分布式模型。

(1)集总式模型

集总式水文模型(Lumped Hydrologic Model)，即不考虑水文现象或要素空间分布，将整个流域作为一个整体进行研究的水文模型。20世纪60年代至80年代，集总式水文模型得到了蓬勃的发展，其中包括新安江模型、水箱模型和SWMM模型等。

新安江模型是由河海大学的赵人俊教授提出的，其在国内有很高的知名度，同时也有较多的应用案例。该模型在产流计算时采用蓄满产流方法，坡面汇流采用单位线法，河道汇流采用马斯京根分河段演算法。朱炬明将新安江与SWAT和BTOPMC模型的应用比较；赵丽平等为了解决一般概念性水文模型参数率定结果不稳定的问题,以新安江模型为例，提出了新安江模型日模参数的线性化率定方法。

水箱模型又称坦克模型，是由日本管原正巳博士于1961年提出，其把降雨转换为径流的复杂过程，描述为蓄水容量与流出的关系进行来模拟，其简化了复杂的流域汇流过程，在全球范围内都有较大的推广。Dressman等人利用混合式的水箱模型对物质的溶解和吸收速率进行了研究；Chen等人利用分布式的水箱模型研究了稻田内降雨与径流汇流间的相互关系。

SWMM模型是由美国环境保护署开发的暴雨洪水管理模型。模型中对各个子流域进行概化，把每个子流域作为单独的水力学单元，每个子流域最终汇流到唯一的出水口，各个出水口通过假定的连接点进行连接。SWMM模型在地表汇流计算时，将各个部分近似看作非线性水库来处理，最后通过联立曼宁方程和连续方程来求解。Ren等人对SWMM模型在霞凝港区进行了应用和研究；余嵘针对城镇内涝问题，以陕西省汉中市某镇为例，利用SWMM模型合成了该地区的典型降雨过程线。

(2)分布式与半分布式模型

分布式水文模型是通过水循环的动力学机制来描述和模拟流域水文过程的数学模型。分布式水文模型由于具有明确物理意义的参数结构和对空间分异性的全面反映，可以更加准确详尽的描述和反映流域内真实的水文过程。分布式水文模型主要包括TOPMODEL模型、SHE模型、FloodArea模型和SWAT模型等。

TOPMODEL模型是一个以地形为基础的半分布式水文模型，由Beven和Kirkby由1979年提出，模型的理论基础是变动产流面积的概念，描述的是结构化水流运动。模型把全流域按DEM网格进行分块，大的流域分又被分为若干个子流域，通过等流时线法进行汇流演算。TOPMODEL模型应用广泛，李抗彬等对TOPMODEL模型的蒸发产流模块以及汇流模块进行改进，以适应于半湿润地区或半干旱地区的径流过程模拟；Niu等人将TOPMODE模型集成到了全球气候模拟模型中。

SHE模型是第一个具有代表性的分布式水文模型，由英国、法国和丹麦的科学家联合研制而成，该模型主要的径流汇流物理过程都使用偏微分方程的差分形式来表示，其中包括质量守恒、能量守恒和动量守恒等方程。其相关应用研究也有许多，Mcmichael等人便将SHE模型应用于半干旱灌丛盆地径流的非确定性预测当中。

SWAT模型是由美国农业部农业中心于1998年开发出的流域尺度模型，模型功能全面，其在流域水量平衡、地表径流、农药输移以及污染物的输移等方面都有较为广泛的应用。Spruill等利用SWAT模型模拟了小流域的日月径流变化，Jayakrishnan等将SWAT模型应用到了水资源的管理中。

FloodArea模型是由德国Geomer公司研制的基于GIS的水动力模型。该模型可以于界定洪水淹没范围和淹没深度，其可以以插件的方式直接与ArcGIS集成，方便业务人员进行洪水淹没模拟。FloodArea在我国应用广泛，Xue等人便将FloodArea模型应用到了城市内涝的模拟当中。

上述现有技术主要存在以下缺点：

集总式水文模型往往先在实验室中进行模拟，再推广到实际流域中，故对于不同的流域特征需调整不同的模型参数。由于集总式水文模型未立足于径流汇流的物理机理，故其无法模拟连续时间段内的流域水文特征；同时，统计型模型的模拟精度不高，已无法满足现今对降雨径流快速实时模拟的要求。

分布式与半分布式水文模型将水动力学引入到实际汇流模拟中，能研究流域出水口或流域各地点的径流汇流。因此在水文模拟的精度上，相较于集总式水文模型有了进一步的提升。但多数分布式水文模型计算量较大，并运行在CPU平台上，其计算效率成为进一步应用的瓶颈。以FloodArea模型为例，对于一个中等大小的流域，24小时径流汇流计算一般需要计算几十分钟甚至一到两个小时，严重影响了淹没风险的预报和评估工作。

另一方面，大部分物理模型未充分考虑流域径流在模拟时间段内的持续水量收入，如上游来水、地下水和地表水的水量交换等，很难反应流域在一定时期内、一定条件下的水位保持相对稳定的现实规律。

发明内容

针对现有主流基于物理机制的汇流模型未充分考虑模拟时间段内的流域持续水量收入，本方案提出了稳定水流场概念，并在径流模拟中增加水流混合模拟环节，实现了栅格模式下的流域稳定水流场建立；针对汇流模型计算量大，将径流汇流模型进行并行化设计，以提高径流汇流计算的时间效率，实现流域径流汇流模型的实时性。

具体而言，本发明提供了如下技术方案：

一方面，本发明提供了一种基于稳定水流场的流域地表径流汇流模拟方法，该方法包括：

S1、建立栅格模式下的N-S方程，其中，每一栅格单元仅能与周围8个相邻栅格进行水流交换，并将N-S方程分解为平流项、压力项和外力项；

S2、在径流汇流模型中进行水流混合模拟，其中，在每一时间片上，每一栅格单元均计算与周围8邻域栅格的总动量和总水量，以得到所述每一栅格单元新的流速；

S3、利用所述S2得到的径流汇流模型对流域内水流场进行迭代，直至整个流域在所有位置上水深和流速不再继续增加，并保持相对稳定，从而获得稳定水流场。

优选地，所述S1中，所述每一栅格单元的水流来源包括水流流入、降雨补给等；水流指出包括水流流出、植被截留、下渗、蒸发等。

优选地，所述下渗的下渗率为：

式中，S、A为反映下渗影响因素的参数；

进一步优选地，同一位置上的植被截留、蒸发对应的蒸发率在模拟时间段内保持不变。

优选地，所述S1中的平流项，可采用如下方式进行模拟：

使用偏移矢量分别标识中心栅格像元的8个相邻像元，则邻域栅格流入到中心栅格的水量比例p为：

式中，为偏移矢量，为邻域栅格的水流速度矢量，且该水流速度归一化到[-1,1]之间；

然后，计算中心栅格时间片上的最终速度v′_c：

υ′_C＝(p_c+Δp)/m'

(8)其中，p_c为初始动量，Δ_p为动量增量，m′为平流计算更新后的水量。

优选地，所述S1中的压力项，采用如下方式进行模拟：

邻域栅格对中心栅格的压力梯度产生的速度矢量增量ΔV表示为：

中心栅格流入或流出到特定邻域栅格的水量Q为：

式中，D为水流交换率，表示偏移矢量为的邻域栅格与中心栅格的水体高度差，dr为邻域栅格单元与中心栅格在水平面上的水流距离相关的常量，为偏移矢量，α为正数常量。

优选地，所述S1中的外力项，通过以下方式模拟：

水速变化表示为：

式中，k为大于0的比例系数；

将水深与水速结合，当水深低于给定阈值时，将水速乘以衰减系数以减小水流速度，衰减系数ε为：

其中，σ为根据流域特征待定的摩擦比例系数，水流摩擦力的上界深度dmax，下界深度为dmin，中心栅格水深为d。

优选地，所述S2中，所述时间片的时长Δt为：

Δt＝C/μ

式中，Δt为计算时间步长，单位为s；C表示流域栅格的分辨率，单位为m；μ为最大流速，单位为m/s。

另一方面，本发明还提供了一种基于稳定水流场的流域地表径流汇流模拟实时计算装置，所述装置利用GPU进行并行计算，并将建模过程设置为栅格单元可独立运行的计算单元；所述装置包括：

稳定水流场生成模块，用于设定流域稳定水源供应，利用径流汇流模型进行迭代计算，并获得稳定水流场；

流域面雨量序列模块，利用时间序列的气象观测数据，插值获得流域面雨量序列；

流域地表径流汇流并行计算模块，基于稳定水流场，将径流汇流模拟的时间轴分解为时间片元，基于所述GPU进行径流汇流迭代计算。

优选地，对于每一时间片元，径流汇流计算依据以下方式进行：

第一步、计算水量和水速的增量；第二步、平流项模拟，依据每一个栅格单元的水速大小和方向，在时间片内计算水流的流入和流出量，并计算水流混合后新的速度矢量；第三步、水流混合模拟，对流域各位置流速大小和方向进行修正。

优选地，所述水源供应的供应量由流域空间特征和时间特征确定，并利用实际的水文数据进行率定。

与现有技术相比，本发明技术方案能够更加精确地模拟水域情况，为流域暴雨洪涝灾害的快速模拟和评估提供辅助决策支持，具有较高的应用价值。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其它的附图。

图1为本发明实施例的中心栅格及8邻域像元的偏移矢量；

图2为本发明实施例的稳定水流场生成过程中流速大小和水深变化；

图3为本发明实施例的稳定水流场效果图；

图4为本发明实施例的典型样例的降水时间分配直方图；

图5为本发明实施例的水文站模拟水位变化曲线与水文站观测水位变化曲线典型样例,；

图6为本发明实施例的方法流程图。

具体实施方式

下面结合附图对本发明实施例进行详细描述。应当明确，所描述的实施例仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其它实施例，都属于本发明保护的范围。本领域技术人员应当知晓，下述具体实施例或具体实施方式，是本发明为进一步解释具体的发明内容而列举的一系列优化的设置方式，而该些设置方式之间均是可以相互结合或者相互关联使用的，同时，下述的具体实施例或实施方式仅作为最优化的设置方式，而不作为限定本发明的保护范围的理解。

实施例1：

本发明提供了一种基于稳定水流场的流域地表径流汇流模拟方法，结合图6，该方法包括：

S1、建立栅格模式下的N-S方程，其中，每一栅格单元仅能与周围8个相邻栅格进行水流交换，并将N-S方程分解为平流项、压力项和外力项；

S2、在径流汇流模型中进行水流混合模拟，其中，在每一时间片上，每一栅格单元均计算与周围8邻域栅格的总动量和总水量，以得到所述每一栅格单元新的流速；

S3、利用所述S2得到的径流汇流模型对流域内水流场进行迭代，直至整个流域在所有位置上水深和流速不再继续增加，并保持相对稳定，从而获得稳定水流场。

优选地，所述S1中，所述每一栅格单元的水流来源包括水流流入、降雨补给等；水流指出包括水流流出、植被截留、下渗、蒸发等。

优选地，所述下渗的下渗率为：

式中，S、A为反映下渗影响因素的参数；

进一步优选地，同一位置上的植被截留、蒸发对应的蒸发率在模拟时间段内保持不变。

优选地，所述S1中的平流项，可采用如下方式进行模拟：

使用偏移矢量分别标识中心栅格像元的8个相邻像元，则邻域栅格流入到中心栅格的水量比例p为：

式中，为偏移矢量，为邻域栅格的水流速度矢量，且该水流速度归一化到[-1,1]之间；

然后，计算中心栅格时间片上的最终速度v′_c：

υ′_C＝(p_c+Δ_p)/m′

(8)其中，p_c为初始动量，Δ_p为动量增量，m′为平流计算更新后的水量。

优选地，所述S1中的压力项，采用如下方式进行模拟：

邻域栅格对中心栅格的压力梯度产生的速度矢量增量ΔV表示为：

中心栅格流入或流出到特定邻域栅格的水量Q为：

式中，D为水流交换率，表示偏移矢量为的邻域栅格与中心栅格的水体高度差，dr为邻域栅格单元与中心栅格在水平面上的水流距离相关的常量，为偏移矢量，a为正数常量。

优选地，所述S1中的外力项，通过以下方式模拟：

水速变化表示为：

式中，k为大于0的比例系数；

将水深与水速结合，当水深低于给定阈值时，将水速乘以衰减系数以减小水流速度，衰减系数ε为：

其中，σ为根据流域特征待定的摩擦比例系数，水流摩擦力的上界深度dmax，下界深度为dmin，中心栅格水深为d。

优选地，所述S2中，所述时间片的时长Δt为：

Δt＝C/μ

式中，Δt为计算时间步长，单位为s；C表示流域栅格的分辨率，单位为m；μ为最大流速，单位为m/s。

实施例2：

在本实施例中，以更为详细的方式阐述本发明的技术方案。

本发明采用的N-S方程算法如下：

纳维-斯托克斯方程(简称N-S方程)，为汇流模型的控制方程。N-S方程是水动力学中描述流体的运动的控制方程，是根据牛顿第二定理推理得到。因为径流汇流中的水流运动是同性不可压的流体运动。因此N-S公式为：

公式(1)中u为速度场，t为时间，ρ为流体密度常数，P为压力，ν为运动粘性系数，F为外力通常是作用于流体上的质量力。并且在公式中运用到了(微分算数符)，其可运用在梯度、散度和旋度中。如表示梯度，表示散度，表示旋度。

根据梯度和散度的定义和原理，公式(1)中的有关梯度和旋度部分改写为：

二维水动力模型模拟汇流过程是在二维平面之上且降雨水流的粘性较小,因此汇流模拟过程中的水流可看作是无粘性的流体，即v＝0。则公式(1)简写为:

公式(4)中每项都有其实际的物理意义：

平流项:方程右边第一项。此项表示流体本身及属性随速度而迁移，是流体各属性往四周的分配过程。

压力项:方程右边第二项。流体分子是***,它们会挤压和碰撞。压强受到外力作用改变不均匀时，产生压强梯度，进而产生压力。因为压力是单位力,所以在流体受到的任何压力都会存在加速度。方程第二项就代表压力梯度产生的加速度。

外力项:方程右边第三项。由外力施加而产生的加速度。这些力可以是只对局部有作用的局部力，也可以是对全局都有作用的外部力。

在一个具体的实施例中，本发明提出的流域稳定水流场，通过如下方式建立：

流域稳定水流场的建立关键在于合理的径流汇流模型。本发明以N-S方程为指导方程，首先对栅格模式下的N-S方程进行了具体实现；然后在径流汇流模拟中加入水流混合模拟，对水流场的速度大小和方向进行修正。通过为流域每一个网格单元提供稳定均匀的水量增量，以模拟流域的上游来水、地下水和地表水交换等，最终建立流域稳定水流场。

(1)栅格模式下的N-S方程具体实现

纳维-斯托克斯方程，描述粘性不可压缩流体动量守恒的运动方程。该方程只有在极少数简单流动的情况下才有解析解，绝大部分流动都不能直接对N-S方程解析。另一方面，本发明的径流汇流模拟是以栅格模式为基础进行汇流演算。在栅格模式中，每一个栅格单元仅能与周围8个相邻栅格进行水流交换。为使模型更真实地反演实际的径流汇流过程，将N-S方程分解为平流项、压力项和外力项，针对栅格的特点探讨了实现方法。在具体实现中，将其设计为栅格单元可独立运算的形式，为后续GPU实时计算奠定基础。

A、栅格单元水流的水量收支

在地表径流汇流过程中，栅格单元的水流主要来源为水流流入、降雨补给；水流支出包括水流流出、植被截留、下渗和蒸发等。依据水量平衡原理，在同一时间段内，水量增量＝降水补给量+流入量-流出量–植被截留–下渗量–蒸发量。

GPU径流模拟以时间片为单位，需要将各项水量收支分摊到各时间片上。其中，降水补给数据来源于气象观测站、雷达降水产品等。作为目前主要的降水数据源，气象观测数据一般以小时累计雨量为计量单位。假定在两个观测时间点之间降水时间分配均匀，则小时降水量除以1个小时的时间片数为单个时间片的降水补给量。由于气象观测数据在整个流域中离散分布，需要将其空间插值得到整个流域的降水栅格，并逐小时更新。

植被截留、下渗和蒸发为地表径流损耗的重要途径。在模型中，同一位置上的植被截留、蒸发率在模拟时间段内保持不变，下渗率采用如下方式进行模拟：

式中，S、A为反映下渗影响因素的参数，与土壤性质有关。

在精度要求较低或者无土地覆盖数据的场合下，初始下渗率、植被截留率和蒸发率一般设置为常数；在精度要求高的场合下，需要叠加流域的土地覆盖数据并设置不同土地覆盖类型的初始下渗率、植被截留率和蒸发率。

单位时间内的水流流入量和流出量可依据中心栅格和邻域栅格的水流速度矢量计算得到，具体在平流项环节中实现。

B、栅格模式下的平流项模拟

在N-S方程中，平流项表示流体本身及属性随速度而迁移。以地表径流为例，平流项表示在当前流场中，每一个栅格单元在单位时间内流入流出的水量及其水速的变化。在GIS的水文分析中，常采用单流向法进行水流方向模拟，这与实际的径流汇流过程存在一定的误差。本发明采用多流向法对栅格水流进行分配。为便于模型计算的并行化，以中心栅格为研究对象计算其流入和流出的水量，并依据动量守恒定律计算水流交换后新的流速大小。

为定量研究水流的流动，使用偏移矢量分别标识中心栅格像元的8个相邻像元。偏移矢量采用笛卡尔坐标系，其原点为中心栅格，x方向向右为正，y方向向上为正。设偏移矢量为则中心栅格及8邻域像元的偏移矢量如图1所示。

设邻域栅格的水流速度矢量为并将水流速度归一化到[-1,1]之间，则该栅格流入到中心栅格的水量比例p为：

可以证明，式6满足质量守恒定律。因此，利用该式进行平流项计算时，只需要计算每一个栅格单元的流入量和流出量，而不考虑邻域栅格的水量增加或扣减，从而允许以栅格单元为基本运算单元进行并行计算。

另外，式6只适用于归一化后的水流速度，并将1m/s作为水流流入流出的最大值。如果邻域栅格指向中心栅格的水流速度达到-1或1m/s，则水流全部流入到中心栅格；如果中心栅格水流速度达到-1或1m/s，则中心栅格的流水全部流出。因此，本发明的径流模拟模型每个时间片的时长Δt为：

Δt＝C/μ

(7)式中，Δt为计算时间步长，单位为s；C表示流域栅格的分辨率，单位为m；μ为最大流速，单位为m/s。理论上a取1m/s。但模型中使用了归一化流水速度，比真实水速偏小，其μ值一般设定为小于1的常数，具体依据流域状况通过率定确定。

计算出中心栅格水量的同时，可依据动量守恒原理计算中心栅格时间片上的最终速度v′_C，计算方法为：

υ′_C＝(p_c+Δ_p)/m′

(8)其中，p_c为初始动量，Δ_p为动量增量，m′为平流计算更新后的水量。

C、栅格模式下的压力项模拟

水流随速度的迁移过程中，由于水压不同产生水流内部的压力梯度力。压力梯度力并不是真正意义上了“力”，它其实是由于水压不同而产生的加速度。在GPU径流模拟中，计算单位为栅格像元，压力梯度力采用8邻域近似计算方法。设相邻栅格水体高度差产生的压力梯度ΔG，则：

其中，ΔP为邻域栅格单元与中心栅格的水压差；Δd为邻域栅格单元与中心栅格在水平面上的水流距离。

由液体压强公式可知，水体压强与水体高度成正比，水体局部压力由相邻栅格的水体高度差决定。设中心栅格高程和水深为H₀、W₀，邻域栅格高程为H_n、W_n,则水体高度差ΔH为：

依据式9和式10，邻域栅格对中心栅格的压力梯度产生的速度矢量增量ΔV表示为：

式中，表示偏移矢量为的邻域栅格与中心栅格的水体高度差；a为与水体密度、重力加速度、栅格分辨率等相关的正数常量；dr为邻域栅格单元与中心栅格在水平面上的水流距离相关的常量,定义如下：

压力梯度力产生的水流加速度，使中心栅格的水团与邻域栅格发生水流交换。设水流交换率为D,则中心栅格流入或流出到特定邻域栅格的水量Q定义如下：

依据公式11和13进行压力项的水量增减和流速变化计算。为实现压力项的栅格计算并行化，单个栅格仅计算该栅格的速度和水量变化。由于输入速度和水深数据在时间片的计算中保持不变，相邻栅格的水量交换和为0，满足径流模拟过程中的质量守恒。

D栅格模式下的外力项模拟

由于地表水流在坡面上或河床上流动，地表径流除受到重力作用外，还存在地表摩擦力、相邻水团的作用力等。其中，重力为径流的根本动力。在重力作用下，水流从高处向低处流动。水流从一个栅格流向另一个栅格，在水面高度降低后重力势能部分转换为动能，表现为水流速度的增加；在水面高度升高后部分动能转换为重力势能，表现为水流速度的减小，甚至反向流动。

综合考虑各种外力因素，水速变化近似使用下式表示：

式中，k为大于0的比例系数，是一个综合考虑重力、地表摩擦力等外力的常数系数。

地表水流由于受到地表摩擦力的影响，使水流速度发生垂直变化。在靠近河床的水层水流速度小，上层水流速度大。从整体上看，河道的中心水深,平均流速大,靠岸的地方水浅,平均流速小。在地表径流汇流模拟中，将水速与水深相结合。当水深低于给定阈值时，将水速乘以衰减系数以减小水流速度。设水流摩擦力的上界深度dmax，下界深度为dmin，中心栅格水深为d，衰减系数ε近似定义为：

其中，σ为根据流域特征待定的摩擦比例系数。

(2)水流混合模拟

实践发现，将N-S方程按上述平流项、压力项和外力项运用于流域径流汇流时，随着迭代次数的增加，流域内水流场将变得越来越不稳定，其水深和水速变化剧烈，甚至出现河道上相邻上下游栅格流速方向相反并相互交替的现象。流场不稳定现象首先发生在支流汇入的主河道位置，并逐渐向河道其他位置扩展。不稳定流场不符合真实河道水流运动的客观规律，影响水流在流域中正常迁移，模拟结果与实际状况存在较大误差。

不稳定流场产生的根本原因在于对流场中流速的不正确模拟。在径流模拟过程中，如果不能保持同一河段纵向上的流量稳定，将导致水流在局部位置拥堵或流空，水面高度变化剧烈，流水在河道内部做无意义的快速往返流动。

本发明认为不稳定流场产生的原因在于模型中未充分考虑水流汇聚时水团融合产生的速度变化，并在模型中添加水流混合模拟环节，以实现稳定水流场的建立。具体实现办法为：在每一个时间片上，平流项、压力项和外力项模拟完成后，每一个栅格单元均计算与周围8邻域栅格的总动量和总水量，利用总动量除以总水量得到该栅格单元新的流速。新的流速综合考虑了与相邻水团混合后的速度变化，使流速分布接近更真实的状态。

(3)稳定水流场的生成

现实流域中，降水并不是流域地表径流的唯一水源。除此之外，还存在上游河道、冰川融水、地下水补给等，这些水源为流域持续不断供水，从而保证了流域即使在短期内不降水，其河道水位也能在一定程度上保持稳定。

由于受到数据获取能力的限制，这些稳定水源尤其是地下水源往往在模型中很难精确描述。本发明假定流域内所有位置上均有一个固定的水量供应τ(单位为mm)，利用径流汇流模型对流域内水流场进行不断迭代计算，直至整个流域在所有位置上水深和流速大小不再继续增加并保持相对稳定，则生成了流域的稳定水流场。稳定水流场生成过程中，同一位置上的水深和水速变化曲线如图2所示。其中，图2中的左侧图为流速大小的变化曲线，流速的两次增加分别表示了河道从坡面流变为河道径流的过程；图2中的右侧图，为水深的变化曲线，随着迭代次数的增加，流速大小和水深最终趋于稳定数值，此时表明流域整体的水流收支保持大致平衡。

通过研究发现，稳定水流场的水深和水速分布接近自然流域的真实状态。具体表现在以下几个方面：

A、流域水深分布接近真实状况，流域水深产生的水力坡度与水速大小呈正相关；

B、流域水速方向从上游指向下游；

C、流速在河道横向上变化与真实状况一致，靠近河岸水速较小，河道中心水速较大；

D、流速在河道纵向上变化平稳；在同一河段上(无支流汇入)，狭窄河道水速大，宽阔河道水速小，流量大小保持恒定。

E、将稳定水流场应用于实际的径流汇流模拟，发现在同一河段上(无支流汇入)水深变化曲线基本一致，表明模型算法的稳定性。

图3为本发明的稳定水流场的地图展示效果。

实施例3：

在有一个实施例中，本发明提供了一种流域地表径流汇流模拟实时计算装置。

将GPU通用计算引入到流域地表径流汇流模拟中，实现径流汇流模拟的实时性。其核心内容是将径流汇流演算各步骤，包括水量收支、平流项、压力项、外力项和水流混合，均设计为栅格单元可独立运行的计算单元，并编写GPU并行计算代码。

在一个具体的实施方式中，该装置设置如下：

该装置利用GPU进行并行计算，并将建模过程设置为栅格单元可独立运行的计算单元；所述装置包括：

稳定水流场生成模块，用于设定流域稳定水源供应，利用径流汇流模型进行迭代计算，并获得稳定水流场；

流域面雨量序列模块，利用时间序列的气象观测数据，插值获得流域面雨量序列；

流域地表径流汇流并行计算模块，基于稳定水流场，将径流汇流模拟的时间轴分解为时间片元，基于所述GPU进行径流汇流迭代计算。

流域地表径流汇流模拟的实时计算装置总体设计及使用方式如下：

(1)稳定水流场生成

设定流域稳定水源供应，该水源供应均匀应用于流域内所有栅格像元，并利用径流模型迭代计算直至稳定水流场生成。水源供应量由流域空间特征和时间特征确定，并利用实际的水文数据进行率定。

(2)流域面雨量序列生成

利用气象站时间序列的观测数据，逐小时插值得到流域面雨量序列。在径流汇流模拟过程中，按照时间先后顺序将面雨量输入到模型中进行汇流演算。

(3)流域地表径流汇流并行计算

以稳定水流场为基础，将流域径流汇流模拟的时间轴按公式(7)所示的步长分解为一定数量的时间片元，通过时间片元的GPU计算完成径流汇流模拟过程的迭代。各时间片元的水量收入包含流域面雨量和稳定水流场的水源供应两部分。

对于每一个时间片元，径流模拟分成三步完成。第一步计算水量和水速的增量，具体包括水量收支的计算、压力梯度力和外力产生的水量和水速的增量计算；第二步为平流项模拟，即依据每一个栅格单元的水速大小和方向，在时间片内计算水流的流入和流出量，并依据动量守恒定律计算水流混合后新的速度矢量；第三步为水流混合模拟，对流域各位置流速大小和方向进行修正。

实施例4：

本实施例对本发明的径流汇流模型计算性能和计算精度评价进行了示例性说明。

由于本发明的径流汇流模型采用了GPU通用计算技术，其性能可提高一至两个数量级。为定量评价算法的运行效率，将同一个流域相同分辨率的地形和降水数据代入到本发明的模型和目前广泛使用的Flood Area模型中，使用同一个计算机进行性能测试，其中测试CPU为Intel CORE i5，GPU为NVIDA GTX 980。模型中仅考虑地形和降水数据，测试结果如表1所示。从表1可见，本发明模型运行效率是Flood Area模型至少6倍以上，能完成流域地表径流汇流的实时模拟。

表1不同流域尺度24小时模拟时间对比

在模型计算精度方面，采用实测水文数据对模型进行评价。选取四川某中小河流域为研究区域，首先依据实际水文数据对模型参数进行率定。经率定，该流域μ＝1/3；σ＝0.995；a＝0.08；D＝0.035；k＝0.002；τ＝0.2。在此基础上，利用流域内不同时间段实测气象站点降水数据对测试流域进行汇流模拟，并采集流域内某水文站的模拟水位和观测水位。测试发现，对于所有测试样例，模型模拟的水位时间变化趋势与降水时间变化趋势均呈现正相关，其水位峰值滞后于降水峰值，与实际情况吻合。图4为其中一个典型样例的降水时间分配直方图；图5的蓝色曲线为模型模拟的水位变化曲线，红色曲线为水文站实际观测水位变化曲线。

本领域普通技术人员可以理解实现上述实施例方法中的全部或部分流程，是可以通过计算机程序来指令相关的硬件来完成，所述的程序可存储于一计算机可读取存储介质中，该程序在执行时，可包括如上述各方法的实施例的流程。其中，所述的存储介质可为磁碟、光盘、只读存储记忆体(Read-Only Memory，ROM)或随机存储记忆体(Random AccessMemory，RAM)等。

以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应以权利要求的保护范围为准。

Claims (10)

1.一种基于稳定水流场的流域地表径流汇流模拟方法，其特征在于，所述方法包括：

S1、建立栅格模式下的N-S方程，其中，每一栅格单元仅能与周围8个相邻栅格进行水流交换，并将N-S方程分解为平流项、压力项和外力项；

S2、在径流汇流模型中进行水流混合模拟，其中，在每一时间片上，每一栅格单元均计算与周围8邻域栅格的总动量和总水量，以得到所述每一栅格单元新的流速；

S3、利用所述S2得到的径流汇流模型对流域内水流场进行迭代，直至整个流域在所有位置上水深和流速不再继续增加，并保持相对稳定，从而获得稳定水流场。

2.根据权利要求1所述的方法，其特征在于，所述S1中，所述每一栅格单元的水流来源包括水流流入、降雨补给；水流指出包括水流流出、植被截留、下渗、蒸发。

3.根据权利要求2所述的方法，其特征在于，所述下渗的下渗率为：

<mrow> <msub> <mi>f</mi> <mi>t</mi> </msub> <mo>=</mo> <mi>f</mi> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mn>1</mn> <mn>2</mn> </mfrac> <msup> <mi>St</mi> <mrow> <mo>-</mo> <mn>0.5</mn> </mrow> </msup> <mo>+</mo> <mi>A</mi> </mrow>

式中，S、A为反映下渗影响因素的参数；

同一位置上的植被截留、蒸发对应的蒸发率在模拟时间段内保持不变。

4.根据权利要求1所述的方法，其特征在于，所述S1中的平流项，采用如下方式进行模拟：

使用偏移矢量分别标识中心栅格像元的8个相邻像元，则邻域栅格流入到中心栅格的水量比例p为：

式中，为偏移矢量，为邻域栅格的水流速度矢量，且该水流速度归一化到[-1，1]之间；

然后，计算中心栅格时间片上的最终速度v′_c：

v′_c＝(p_c+Δp)/m′

(8)

其中，p_c为初始动量，Δp为动量增量，m′为平流计算更新后的水量。

5.根据权利要求1所述的方法，其特征在于，所述S1中的压力项，采用如下方式进行模拟：

邻域栅格对中心栅格的压力梯度产生的速度矢量增量ΔV表示为：

中心栅格流入或流出到特定邻域栅格的水量Q为：

式中，D为水流交换率，表示偏移矢量为的邻域栅格与中心栅格的水体高度差，dr为邻域栅格单元与中心栅格在水平面上的水流距离相关的常量，为偏移矢量，α为正数常量。

6.根据权利要求1所述的方法，其特征在于，所述S1中的外力项，通过以下方式模拟：

水速变化表示为：

式中，k为大于0的比例系数；

将水深与水速结合，当水深低于给定阈值时，将水速乘以衰减系数以减小水流速度，衰减系数ε为：

<mrow> <mi>&amp;epsiv;</mi> <mo>=</mo> <mfenced open = "{" close = ""> <mtable> <mtr> <mtd> <mrow> <mi>&amp;sigma;</mi> <mo>&amp;times;</mo> <msup> <mrow> <mo>(</mo> <mfrac> <mrow> <mi>d</mi> <mo>-</mo> <mi>d</mi> <mi>min</mi> </mrow> <mrow> <mi>d</mi> <mi>max</mi> <mo>-</mo> <mi>d</mi> <mi>min</mi> </mrow> </mfrac> <mo>)</mo> </mrow> <mn>2</mn> </msup> </mrow> </mtd> <mtd> <mrow> <mi>d</mi> <mo>&amp;le;</mo> <mi>d</mi> <mi>max</mi> </mrow> </mtd> </mtr> <mtr> <mtd> <mi>&amp;sigma;</mi> </mtd> <mtd> <mrow> <mi>d</mi> <mo>&gt;</mo> <mi>d</mi> <mi>max</mi> </mrow> </mtd> </mtr> </mtable> </mfenced> </mrow>

其中，σ为根据流域特征待定的摩擦比例系数，水流摩擦力的上界深度dmax，下界深度为dmin，中心栅格水深为d。

7.根据权利要求1所述的方法，其特征在于，所述S2中，所述时间片的时长Δt为：

Δt＝C/μ

式中，Δt为计算时间步长，单位为s；C表示流域栅格的分辨率，单位为m；μ为最大流速，单位为m/s。

8.一种基于稳定水流场的流域地表径流汇流模拟实时计算装置，其特征在于，所述装置利用GPU进行并行计算，并将建模过程设置为栅格单元可独立运行的计算单元；所述装置包括：

稳定水流场生成模块，用于设定流域稳定水源供应，利用径流汇流模型进行迭代计算，并获得稳定水流场；

流域面雨量序列模块，利用时间序列的气象观测数据，插值获得流域面雨量序列；

流域地表径流汇流并行计算模块，基于稳定水流场，将径流汇流模拟的时间轴分解为时间片元，基于所述GPU进行径流汇流迭代计算。

9.根据权利要求8所述的装置，其特征在于，对于每一时间片元，径流汇流计算依据以下方式进行：

第一步、计算水量和水速的增量；第二步、平流项模拟，依据每一个栅格单元的水速大小和方向，在时间片内计算水流的流入和流出量，并计算水流混合后新的速度矢量；第三步、水流混合模拟，对流域各位置流速大小和方向进行修正。

10.根据权利要求8所述的装置，其特征在于，所述水源供应的供应量由流域空间特征和时间特征确定，并利用实际的水文数据进行率定。