一种基于下垫面特征的重配置产流模拟方法与流程

本发明属于水文技术领域，具体涉及一种基于下垫面特征的重配置产流模拟方法。

背景技术：

近年来，受气候变化影响，由局地强降水造成的中小河流突发性洪水频繁发生，已成为造成人员伤亡的主要灾种。我国河流众多，流域面积200至3000km²的中小流域近9000个。中小流域由于通常处于地形复杂、坡度陡峻的偏远山区，急促上涨的洪水极易形成危害当地的居民人身安全和社会经济的洪涝灾害，因此对中小流域突发性洪水进行准确可靠地预警预报成为亟待解决的重要问题。

水文模型是进行水文预报的重要工具，而产流计算是构建水文模型时的关键一环。目前常用的产流计算方法主要包括蓄满产流法与超渗产流法两种。其中蓄满产流主要适用于植被茂盛，坡度平缓的山间谷地，而超渗产流主要适用于山坡陡峭，下渗能力常常小于降雨强度的坡地。对于山区性流域，下垫面特征空间变异性较大，山坡与谷地相间分布。常常具有河流上游山岭陡峭，下游地势平缓；远离河流坡度陡峭，靠近河流地势平坦的特征。因此采用单一的产流计算方法往往不能很好地考虑到流域下垫面特征的空间变异性对于产流过程的影响，从而对水文模型的产流计算带来一定的偏差，进而影响到最终的预警预报效果，为山洪灾害防治带来一定的隐患。基于遥感、地理信息以及数字流域等技术的发展，采用数值矩阵描述地表高程变化的栅格数字高程模型(dem，digitalelevationmodel)逐步成熟，并得到了广泛的应用。尤其是在地形复杂的山区性中小流域中，dem数据因其能够较为准确地考虑流域内地形变化而具有重要的应用价值。同时全球植被土壤遥感资料也在不断地丰富完善，进而为考虑流域内土壤植被对于产流过程的影响提供了重要的支撑。如何利用dem数据和相应的植被土壤数据来对山区性流域中的产流计算方法进行完善，从而在一定程度上考虑下垫面特征空间变异性对于产流计算的影响，也是水文模型走向分布化，精细化过程中的重点和难点之一。

为了进一步促进水文模型中产流计算的发展，需要更深入理解下垫面特征空间变异性对于产流计算的影响，研究基于下垫面特征的重配置产流模拟方法。

超渗产流是指地面径流产生的原因是同期的降水量大于同期植物截留量、填洼量、雨期蒸发量及下渗量等的总和，多余出来的水量产生了地面径流。一般来说，植物截留量、雨期蒸发量、填洼量一般较小；而下渗量一般较大、且变化幅度也很大，它从初渗到稳渗、在时程上具有急变特性，空间上也具有多变的特性。下渗量的空间变化一般表现为：土壤质地疏松，植被茂盛，坡度平缓时下渗能力强；土壤质地紧密，植被稀疏，坡度陡峭时，下渗能力弱；由此可见，地形，植被与土壤分布对下渗能力产生的影响很大，进而影响到超渗产流发生的主要区域范围。蓄满产流是因降水使土壤包气带和饱水带基本饱和而产生径流的方式，是降雨径流的产流方式之一。在降雨量较充沛的湿润、半湿润地区，地下潜水位较高，土壤前期含水量大，由于一次降雨量大，历时长，降水满足植物截留、入渗、填洼损失后，损失不再随降雨延续而显著增加，土壤基本饱和，从而广泛产生地表径流。此时的地表径流不仅包括地面径流，也包括壤中流和其它形式的浅层地下水产流。蓄满产流方式往往不能在坡度陡峭的区域上普遍实现，在平缓地带则容易发生。因此地形，植被覆盖和土壤质地是影响蓄满和超渗发生区域的主要因素，若是忽略这些因素对于产流过程的影响，则会导致产流过程模拟的不准确，从而为山区性中小流域的山洪防治留下隐患。

因而对于下垫面特征空间变异性的忽视不利于水文模型产流方法的发展。

针对以上不足，如何考虑下垫面特征空间变异性对于产流过程的影响，判别发生蓄满产流和超渗产流的区域，在不同的区域中采用合适的产流方法进行产流计算，从而构建一种基于下垫面特征的重配置产流模拟方法，正是发明人需要解决的问题。

技术实现要素：

发明目的：本发明的目的在于提供一种基于下垫面特征的重配置产流模拟方法，具有数据来源稳定可靠、计算效率高、结果客观合理等优点，有利于基于下垫面特征的产流计算方法的快速计算，值得推广。

技术方案：为实现上述目的，本发明提供如下技术方案：

一种基于下垫面特征的重配置产流模拟方法，包括以下步骤:

s1，封装产流计算模式，并分别对其进行编号；

s2，基于流域dem数据，计算得到地形指数栅格；

s3，提取得到流域中的植被栅格与土壤栅格；

s4，基于地形指数栅格，土壤栅格与植被栅格将流域划分为gx，dx，gc与fc四个区域，在不同区域中分别匹配不同的产流模式进行产流计算。

进一步地，所述的s1包括如下步骤：

1)基于概念性方法，构建第一蓄满产流模式，gx；

2)基于达西下渗原理，构建第二蓄满产流模式，dx；

3)基于格林安普特下渗原理，构建第一超渗产流模式，gc；

4)基于菲利普下渗原理，构建第二超渗产流模式，fc。

进一步地，所述的步骤1)中，时段产流量通过下式确定：

当pe≤0或pe+w0≤wm时：

r＝0；

当pe+w0＞wm时：

r＝pe+w0-wm；

式中：r为时段产流量；w0为栅格单元实际的张力水含量；wm为栅格单元张力水容量；pe为净雨量，pe＝p-et，其中p为降雨，et为蒸发，由三层蒸散发原理计算得到。

进一步地，所述的步骤2)中，所述的第二蓄满产流模式为：

r＝ra+rb；

式中：ra为基流，rb为饱和坡面流，在栅格单元的土壤缺水量d＜0的时候产生：

式中：r0为等于零时的流量，为全流域所有栅格单元的土壤缺水量d的平均值，m为参数；

rb＝|d|；

和d的关系如下：

式中：tanβ为该栅格单元的地形坡度，acc为该栅格单元以上的累积汇流面积，∈为全流域所有栅格单元的值的平均。

进一步地，所述的步骤3)中，第一超渗产流模式为：

式中：fm为栅格单元格林安普特下渗速率；k为栅格单元饱和水力传导度；z为栅格单元饱和层厚度；h为栅格单元地面滞水深；为栅格单元湿润锋面的毛管水压力；

r＝pe-fm。

进一步地，所述的步骤4)中，构建第二超渗产流模式：

式中：fp为栅格单元菲利普下渗速率；t为时间；fc为栅格单元稳定下渗率；a为参数；

r＝pe-fp。

进一步地，所述的s2包括以下步骤：

1)利用流域dem数据计算出流域中每一个栅格单元的坡度、流向和汇流累计值，得到坡度栅格raster_slope；

以栅格单元cell为中心，通过周围栅格单元的高程值与该栅格单元的高程值的对比，找出与其相比最低的栅格单元celld，并计算cell和celld之间的高程差dhmax和水平投影距离dis，结合dhmax和dis计算栅格单元cell的坡度β；

β＝dhmax/；

2)将cell作为出流栅格单元，celld作为入流栅格单元，入流栅格单元汇流累计值加1，逐栅格循环，计算出每一个栅格单元中的汇流累计值acc；同时依据cell和celld之间的相对位置关系，确定cell中的流向；

若celld在cell的右方，则cell中的流向设置为1；若celld在cell的右下方，则cell中的流向设置为2；若celld在cell的下方，则cell中的流向设置为3；若celld在cell的左下方，则cell中的流向设置为4；若celld在cell的左方，则cell中的流向设置为5；若celld在cell的左上方，则cell中的流向设置为6；若celld在cell的上方，则cell中的流向设置为7；若celld在cell的右上方，则cell中的流向设置为0；按照以上方法遍历流域中每一个栅格单元，从而得到流向栅格raster_dir和汇流累计栅格raster_acc。

3)结合流域资料和卫星影像确定面积汇流阈值，提取流域中的河道栅格，并将流域中的栅格单元划分为坡面栅格单元与河道栅格单元；

结合流域实际自然地理情况，设置面积汇流阈值t，将raster_acc中acc高于t的栅格单元判定为河道栅格单元，否则认为该栅格单元为坡面栅格单元；将提取得到的河道栅格单元的位置与googleearth上的卫星影像中的实际河道相对比，一般而言，t的值越小，提取得到的河道栅格单元的密度越大，t的值越大，提取得到的河道栅格单元的密度越小。反复调整t，使得提取出来的河道栅格单元的位置与卫星影像中的实际河道的位置基本重合，提取出流域中的河道栅格raster_river；同时依据strahler原理：初始发源的河道级别为1，相同级别的河道交汇后形成的新河道级别提升1级，不同级别的河道交汇后形成的新河道的级别与交汇前的高级别的河道一致。

4)基于坡度栅格raster_slope与汇流累计栅格raster_acc，在每一个坡面栅格单元中计算地形指数，得到地形指数栅格raster_ti；

在每一个坡面栅格单元中：

式中：acc为编号坡面栅格单元中的汇流累积值；β为编号为坡面栅格单元中的坡度值；按照以上方法，逐个计算所有坡面栅格单元中的地形指数，得到地形指数栅格raster_ti。

进一步地，所述的s3包括以下步骤：

1)基于流域中的植被栅格，将流域中植被类型划分为林地，灌木，草地，耕地，水面，并分别进行赋值100，70，50,30,0；构建得到植被下渗能力栅格raster_plantinfli；

2)基于流域中的土壤栅格，将流域中土壤类型划分为砂土，壤土，黏壤土，黏土；并分别进行赋值100，80，60,40,20；构建得到土壤下渗能力栅格raster_soilinfli。

进一步地，所述的s4包括以下步骤：

1)将每一个栅格单元中地形指数，植被下渗能力和土壤下渗能力相乘，得到流域下垫面下渗特性值infiltration，遍历流域中所有的栅格单元，得到流域下垫面下渗特性栅格raster_infiltration；

2)基于流域自然地理特征设置下渗特性阈值t_infi，infiltration大于t_infi且地形指数大于10的栅格单元中匹配gx产流模式；infiltration大于t_infi且地形指数小于10的栅格单元中匹配dx产流模式；infiltration小于t_infi且地形指数小于10的栅格单元中匹配gc产流模式；infiltration小于t_infi且地形指数大于10的栅格单元中匹配fc产流模式。

有益效果：与现有技术相比，本发明的一种基于下垫面特征的重配置产流模拟方法，以影响产流的物理因子为基础，量化了地形，植被覆盖和土壤质地对于产流过程的影响作用，在不同区域采用合适的产流方法进行计算，进而提出一种基于下垫面特征的产流计算方法重配置技术；这样既保证了计算结果的精度与可靠性，同时解决了下垫面特征空间变异明显的山区性中小流域中的产流计算问题；且本方法主要应用流域数字高程模型，数据来源稳定可靠，方法中变量之间的函数关系明确，有利于流域中适用于不同产流计算方法区域的快速自动判别，通过数字流域技术以简化提取步骤，同时，保证了结果的客观合理性，有利于基于下垫面特征的产流计算方法的直接调用，可以进一步促进数字水文学以及山区性中小流域山洪防治研究的深入发展。

附图说明

图1计算流程示意；

图2马渡王流域高程分布栅格；

图3马渡王流域坡度栅格；

图4马渡王流域河道栅格和坡面栅格；

图5马渡王流域地形指数栅格；

图6马渡王流域植被分布；

图7马渡王流域土壤分布；

图8马渡王流域产流模式重配置分布示意图。

具体实施方式

下面结合附图和具体实施例对本发明作更进一步的说明。

如图1-8所示，一种基于下垫面特征的重配置产流模拟方法，包括以下步骤：

s1、封装产流计算模式，并分别对其进行编号，具体包括以下步骤：

1)基于概念性方法，构建第一蓄满产流模式，gx；

当pe≤0或pe+w0≤wm时：

r＝0

当pe+w0＞wm时：

r＝pe+w0-wm

式中：r为时段产流量；w0为栅格单元实际的张力水含量；wm为栅格单元张力水容量；pe为净雨量，pe＝p-et，其中p为降雨，et为蒸发，由三层蒸散发原理计算得到。

2)基于达西下渗原理，构建第二蓄满产流模式，dx；

r＝ra+rb

式中：ra为基流，rb为饱和坡面流，在栅格单元的土壤缺水量d＜0的时候产生。

式中：r0为等于零时的流量，为全流域所有栅格单元的土壤缺水量d的平均值，m为参数。

rb＝|d|

和d的关系如下：

式中：tanβ为该栅格单元的地形坡度，acc为该栅格单元以上的累积汇流面积，∈为全流域所有栅格单元的值的平均。

3)基于格林安普特下渗原理，构建第一超渗产流模式，gc；

式中：fm为栅格单元格林安普特下渗速率；k为栅格单元饱和水力传导度；z为栅格单元饱和层厚度；h为栅格单元地面滞水深；为栅格单元湿润锋面的毛管水压力。

r＝pe-fm

4)基于菲利普下渗原理，构建第二超渗产流模式，fc；

式中：fp为栅格单元菲利普下渗速率；t为时间；fc为栅格单元稳定下渗率；a为参数。

r＝pe-fp

s2、基于流域dem数据，计算得到地形指数栅格，具体包括以下步骤：

1)利用流域dem数据(图2)计算出流域中每一个栅格单元的坡度、流向和汇流累计值，得到坡度栅格raster_slope(图3)；

以栅格单元cell为中心，通过周围栅格单元的高程值与该栅格单元的高程值的对比，找出与其相比最低的栅格单元celld，并计算cell和celld之间的高程差dhmax和水平投影距离dis，结合dhmax和dis计算栅格单元cell的坡度β；

β＝dhmax/

2)将cell作为出流栅格单元，celld作为入流栅格单元，入流栅格单元汇流累计值加1，逐栅格循环，计算出每一个栅格单元中的汇流累计值acc；同时依据cell和celld之间的相对位置关系，确定cell中的流向；

若celld在cell的右方，则cell中的流向设置为1；若celld在cell的右下方，则cell中的流向设置为2；若celld在cell的下方，则cell中的流向设置为3；若celld在cell的左下方，则cell中的流向设置为4；若celld在cell的左方，则cell中的流向设置为5；若celld在cell的左上方，则cell中的流向设置为6；若celld在cell的上方，则cell中的流向设置为7；若celld在cell的右上方，则cell中的流向设置为0。

按照以上方法遍历流域中每一个栅格单元，从而得到流向栅格raster_dir和汇流累计栅格raster_acc。

3)结合流域资料和卫星影像确定面积汇流阈值，提取流域中的河道栅格，并将流域中的栅格单元划分为坡面栅格单元与河道栅格单元；

结合流域实际自然地理情况，设置面积汇流阈值t，将raster_acc中acc高于t的栅格单元判定为河道栅格单元，否则认为该栅格单元为坡面栅格单元；将提取得到的河道栅格单元的位置与googleearth上的卫星影像中的实际河道相对比，一般而言，t的值越小，提取得到的河道栅格单元的密度越大，t的值越大，提取得到的河道栅格单元的密度越小。反复调整t，使得提取出来的河道栅格单元的位置与卫星影像中的实际河道的位置基本重合，提取出流域中的河道栅格raster_river(图4)；同时依据strahler原理：

初始发源的河道级别为1，相同级别的河道交汇后形成的新河道级别提升1级，不同级别的河道交汇后形成的新河道的级别与交汇前的高级别的河道一致。

4)基于坡度栅格raster_slope与汇流累计栅格raster_acc，在每一个坡面栅格单元中计算地形指数，得到地形指数栅格raster_ti(图5)；

在每一个坡面栅格单元中：

式中：acc为编号坡面栅格单元中的汇流累积值；β为编号为坡面栅格单元中的坡度值；按照以上方法，逐个计算所有坡面栅格单元中的地形指数，得到地形指数栅格raster_ti。

s3、基于流域中的植被栅格(图6)与土壤栅格(图7)，计算得到植被下渗能力栅格和土壤下渗能力栅格，具体包括以下步骤：

1)基于流域中的植被栅格，将流域中植被类型划分为林地，灌木，草地，耕地，水面，并分别进行赋值100，70，50,30,0；构建得到植被下渗能力栅格raster_plantinfli；

2)基于流域中的土壤栅格，将流域中土壤类型划分为砂土，壤土，黏壤土，黏土；并分别进行赋值100，80，60,40,20；构建得到土壤下渗能力栅格raster_soilinfli。

s4、基于地形指数栅格，土壤栅格与植被栅格将流域划分为gx，dx，gc与fc四个区域(图8)，在不同区域中分别匹配不同的产流模式进行产流计算，具体包括以下步骤：

1)将每一个栅格单元中地形指数，植被下渗能力和土壤下渗能力相乘，得到流域下垫面下渗特性值infiltration，遍历流域中所有的栅格单元，得到流域下垫面下渗特性栅格raster_infiltration；

2)基于流域自然地理特征设置下渗特性阈值t_infi，infiltration大于t_infi且地形指数大于10的栅格单元中匹配gx产流模式；infiltration大于t_infi且地形指数小于10的栅格单元中匹配dx产流模式；infiltration小于t_infi且地形指数小于10的栅格单元中匹配gc产流模式；infiltration小于t_infi且地形指数大于10的栅格单元中匹配fc产流模式。

实施例1

以陕西省马渡王流域为例，该流域面积1601km²，位于陕西省西部，属暖温带半湿润大陆性季风气候，四季冷暖干湿分明。暴雨中心多集中在流域的中上游。研究区dem原始数据采用美国太空总署(nasa)与国防部国家测绘局(nima)联合提供的90m分辨率srtm(shuttleradartopographymission)数据。植被和土壤数据采用美国马里兰大学发布的全球1km分辨率的土壤类型与土壤利用类型数据。

步骤一、封装产流计算模式，并分别对其进行编号，具体包括以下步骤：

1)基于基于概念性方法，构建第一蓄满产流模式，gx；

当pe≤0或pe+w0≤wm时：

r＝0

当pe+w0＞wm时：

r＝pe+w0-wm

式中：r为时段产流量；w0为栅格单元实际的张力水含量；wm为栅格单元张力水容量；pe为净雨量，pe＝p-et，其中p为降雨，et为蒸发，由三层蒸散发原理计算得到。

2)基于达西下渗原理，构建第二蓄满产流模式，dx；

r＝ra+rb

式中：ra为基流，rb为饱和坡面流，在栅格单元的土壤缺水量d＜0的时候产生。

式中：r0为等于零时的流量，为全流域所有栅格单元的土壤缺水量d的平均值，m为参数。

rb＝|d|

和d的关系如下：

式中：tanβ为该栅格单元的地形坡度，acc为该栅格单元以上的累积汇流面积，∈为全流域所有栅格单元的值的平均。

3)基于格林安普特下渗原理，构建第一超渗产流模式，gc；

式中：fm为栅格单元格林安普特下渗速率；k为栅格单元饱和水力传导度；z为栅格单元饱和层厚度；h为栅格单元地面滞水深；为栅格单元湿润锋面的毛管水压力。

r＝pe-fm

4)基于菲利普下渗原理，构建第二超渗产流模式，fc；

式中：fp为栅格单元菲利普下渗速率；t为时间；fc为栅格单元稳定下渗率；a为参数。

r＝pe-fp

步骤二、基于流域dem数据(图2)，计算得到地形指数栅格，具体包括以下步骤：

1)利用流域dem数据计算出流域中每一个栅格单元的坡度、流向和汇流累计值，得到坡度栅格raster_slope(图3)；

以栅格单元cell为中心，通过周围栅格单元的高程值与该栅格单元的高程值的对比，找出与其相比最低的栅格单元celld，并计算cell和celld之间的高程差dhmax和水平投影距离dis，结合dhmax和dis计算栅格单元cell的坡度β；

β＝dhmax/

2)将cell作为出流栅格单元，celld作为入流栅格单元，入流栅格单元汇流累计值加1，逐栅格循环，计算出每一个栅格单元中的汇流累计值acc；同时依据cell和celld之间的相对位置关系，确定cell中的流向；

若celld在cell的右方，则cell中的流向设置为1；若celld在cell的右下方，则cell中的流向设置为2；若celld在cell的下方，则cell中的流向设置为3；若celld在cell的左下方，则cell中的流向设置为4；若celld在cell的左方，则cell中的流向设置为5；若celld在cell的左上方，则cell中的流向设置为6；若celld在cell的上方，则cell中的流向设置为7；若celld在cell的右上方，则cell中的流向设置为0。

按照以上方法遍历流域中每一个栅格单元，从而得到流向栅格raster_dir和汇流累计栅格raster_acc。

3)结合流域资料和卫星影像确定面积汇流阈值，提取流域中的河道栅格，并将流域中的栅格单元划分为坡面栅格单元与河道栅格单元；

结合流域实际自然地理情况，设置面积汇流阈值t，将raster_acc中acc高于t的栅格单元判定为河道栅格单元，否则认为该栅格单元为坡面栅格单元；将提取得到的河道栅格单元的位置与googleearth上的卫星影像中的实际河道相对比，一般而言，t的值越小，提取得到的河道栅格单元的密度越大，t的值越大，提取得到的河道栅格单元的密度越小。反复调整t，使得提取出来的河道栅格单元的位置与卫星影像中的实际河道的位置基本重合，提取出流域中的河道栅格raster_river(图4)。在马渡王流域，t的值最终设置为20km²；同时依据strahler原理：

初始发源的河道级别为1，相同级别的河道交汇后形成的新河道级别提升1级，不同级别的河道交汇后形成的新河道的级别与交汇前的高级别的河道一致。

4)基于坡度栅格raster_slope与汇流累计栅格raster_acc，在每一个坡面栅格单元中计算地形指数，得到地形指数栅格raster_ti；

在每一个坡面栅格单元中：

式中：acc为编号坡面栅格单元中的汇流累积值；β为编号为坡面栅格单元中的坡度值；按照以上方法，逐个计算所有坡面栅格单元中的地形指数，得到地形指数栅格raster_ti。

按照以上方法，逐个计算所有坡面栅格单元中的地形指数，得到地形指数栅格raster_ti(图5)。

步骤三、基于流域中的植被栅格与土壤栅格，计算得到植被下渗能力栅格和土壤下渗能力栅格，具体包括以下步骤：

1)基于流域中的植被栅格(图6)，将流域中植被类型划分为林地，灌木，草地，耕地，水面，并分别进行赋值100，70，50,30,0。构建得到植被下渗能力栅格raster_plantinfli；

2)基于流域中的土壤栅格(图7)，将流域中土壤类型划分为砂土，壤土，黏壤土，黏土。并分别进行赋值100，80，60,40,20。构建得到土壤下渗能力栅格raster_soilinfli。

步骤四、基于地形指数栅格，土壤栅格与植被栅格将流域划分为gx，dx，gc与fc四个区域(图8)，在不同区域中分别匹配不同的产流模式进行产流计算，具体包括以下步骤：

1)将每一个栅格单元中地形指数，植被下渗能力和土壤下渗能力相乘，得到流域下垫面下渗特性值infiltration，遍历流域中所有的栅格单元，得到流域下垫面下渗特性栅格raster_infiltration；

2)基于流域自然地理特征设置下渗特性阈值t_infi，infiltration大于t_infi且地形指数大于10的栅格单元中匹配gx产流模式；infiltration大于t_infi且地形指数小于10的栅格单元中匹配dx产流模式；infiltration小于t_infi且地形指数小于10的栅格单元中匹配gc产流模式；infiltration小于t_infi且地形指数大于10的栅格单元中匹配fc产流模式。

另外需要说明的是，在上述具体实施方式中所描述的各个具体技术特征，在不矛盾的情况下，可以通过任何合适的方式进行组合。为了避免不必要的重复，本发明对各种可能的组合方式不再另行说明。

以上显示和描述了本发明的基本原理、主要特征及优点。本行业的技术人员应该了解，本发明不受上述实施例的限制，上述实施例和说明书中描述的只是说明本发明的原理，在不脱离本发明精神和范围的前提下，本发明还会有各种变化和改进，这些变化和改进都落入要求保护的本发明范围内。本发明要求保护范围由所附的权利要求书及其等效物界定。