本发明属于水文技术领域，具体涉及一种基于动态沟道的水文响应单元出流过程的计算方法。
背景技术
我国河流众多，流域面积200至3000km2的中小流域近9000个。近年来，受气候变化影响，由局地强降水造成的中小河流突发性洪水频繁发生，已成为造成人员伤亡的主要灾种。中小流域由于通常处于地形复杂、坡度陡峻的偏远山区，一方面山洪具有陡涨陡落的特点，另一方面缺少长期的水文气象观测资料，因此急促上涨的洪水极易形成危害当地的居民人身安全和社会经济的洪涝灾害，对中小流域突发性洪水进行快速预报成为亟待解决的重要问题。
随着遥感、地理信息以及数字流域等技术的发展，基于栅格数字高程模型(dem，digitalelevationmodel)的分布式水文模型以其充分考虑降雨和下垫面条件空间变化的特点，现已成为流域水文模型的发展趋势，尤其是在地形复杂的山区性中小流域中，分布式水文模型比传统的不能考虑流域内地形变化的集总式模型更有优势。在构建分布式的水文模型时，通常将流域划分为若干个正交的网格，以河源点以上的集水区域中若干网格构成的集合作为水文响应单元。在每一个水文响应单元中通过产流模型计算产生的净雨深，然后使用一定的坡面汇流方法得到计算单元出口处的出流过程，并将其作为河源点处的入流，再在河道中进行河道汇流计算得到整个流域出口处的径流过程。在汇流计算中十分重要的一环就是水文响应单元中的坡面汇流演算，这也是分布式水文模型建模的重点和难点之一。
为了进一步促进流域水文模型的发展，需要更深入研究水文响应单元中的坡面汇流演算的计算方法。
坡面汇流是坡面净雨沿着地表向水文响应单元出口断面汇集的过程，是水文径流模型的一项重要内容，坡面汇流计算的合理与否对计算结果有着重要影响。对于汇流计算的研究主要以整个流域或者是河道为研究对象，而对于坡面的汇流规律研究较少。目前对于坡面汇流的计算主要有两种方法，一种是简单地将坡面看作是一个大的水库，降雨落在地面上产生的径流作为这个水库的输入，经过这个水库的调蓄作用后形成水文响应单元出口断面处的出流过程；另一种方法则将流域中每一个网格作为计算单元，使用水力学的方法在栅格之间进行扩散波或者是运动波演算，从而得到河源处的径流过程。这两种方法均事先将流域划分为特征具有鲜明差别的坡地与河道，但是在实际的降雨过程中坡面会随着降雨强度的变化而出现可以让水流快速通过的沟道。这些沟道并不像河道一样常年过水且有稳定的河岸边界，它们只会在降雨过程中出现，同时随着降雨强度的变化而迅速发展和消退。沟道的出现，会使得坡面上原本运动缓慢的水分汇集成束状的水流，并会因为在一般较为陡峻的河源区山坡上受到重力作用而快速运动，其速度将因此提高数个数量级。因此沟道的出现对于坡面中水流的运动具有极其重要的影响。而事先将流域划分为固定的坡地与河道这种做法无法充分地考虑在降雨过程中发展变化的沟道对于汇流的重要影响。
因而对沟道在降雨过程中发展变化的忽视不利于分布式水文模型的发展。
针对以上不足，如何考虑降雨过程中不断动态发展变化的沟道，量化沟道对于水文响应单元中汇流的影响，正是发明人需要解决的问题。
技术实现要素：
为了解决现有技术中存在的不足，本发明提供了一种基于动态沟道的水文响应单元出流过程的计算方法，具有数据来源稳定可靠、计算效率高、结果客观合理等优点，有利于水文响应单元出流过程的快速计算，值得推广。
为解决上述问题，本发明具体采用以下技术方案：
一种基于动态沟道的水文响应单元出流过程的计算方法，其特征在于，包括以下步骤：
步骤1，基于dem数据提取水文响应单元；
步骤2，基于降雨过程中的时段降雨量计算水文响应单元中各个栅格单元的水流动能指标，并由此判别该时段中的坡面栅格和沟道栅格；
步骤3，估算水文响应单元中每一个栅格单元中的汇流流速；
步骤4，在水文响应单元中计算从每一个栅格单元出发，沿着一定的路径运动到水文响应单元出口所需要的时间，即该栅格单元的汇流时间，从而得到汇流时间栅格；
步骤5，基于汇流时间栅格，提取水文响应单元中的汇流单位线；
步骤6，对于降雨过程中的每一个时段重复步骤2到步骤5之间的操作，得到每一个时段中水文响应单元的汇流单位线，并结合降雨过程进行卷积计算，得到水文响应单元的出流过程。
前述的一种基于动态沟道的水文响应单元出流过程的计算方法，其特征在于，所述步骤1中基于dem数据提取水文响应单元，具体包括以下步骤：
步骤1.1，利用流域dem数据计算出流域中每一个栅格单元的坡度、流向和汇流累计值，得到坡度栅格、流向栅格和汇流累计栅格；
以栅格单元cell为中心，通过周围栅格单元的高程值与该栅格单元的高程值的对比，找出与其相比最低的栅格单元celld，并计算cell和celld之间的高程差dhmax和水平投影距离dis，结合dhmax和dis计算栅格单元cell的坡度j：
j＝dhmax/dis
将cell作为出流栅格单元，celld作为入流栅格单元，入流栅格单元汇流累计值加1，逐栅格循环，计算出每一个栅格单元中的汇流累计值acc，同时依据cell和celld之间的相对位置关系，采用d8流向法确定cell中的流向，按照以上方法遍历流域中每一个栅格单元，从而得到坡度栅格raster_slope、流向栅格raster_dir和汇流累计栅格raster_acc；
步骤1.2，基于面积汇流阈值提取流域中的河道栅格，并结合流域资料和卫星影像确定固定河源点；
首先初步假定值t，将raster_acc中acc高于t的栅格单元判定为河道栅格单元，并将提取得到的河道栅格单元的位置与googleearth上的卫星影像中的实际河道相对比，反复调整t，使得提取出来的河道栅格单元的位置与卫星影像中的实际河道的位置重合，从而去确定汇流阈值t，同时提取出流域中的河道栅格，并进一步得到各个河道的固定河源点的集合，即固定河源点的集合points；
步骤1.3，以points中的一个固定河源点point为出口点，结合raster_dir，提取出所有汇流到point的栅格单元，并将这些栅格单元所在的区域作为point所对应的水文响应单元subaisn，遍历points中的所有固定河源点，从而提取出流域中所有的水文响应单元。为简化说明，本文中步骤2至步骤6中均以一个水文响应单元subaisn为例。
前述的一种基于动态沟道的水文响应单元出流过程的计算方法，其特征在于，所述步骤2中基于降雨过程中的时段降雨量计算水文响应单元中各个栅格单元的水流动能指标，并由此判别该时段中的坡面栅格和沟道栅格，具体包括以下步骤：
步骤2.1，计算水文响应单元中该时段内每一个栅格单元内的水流动能指标；
etp＝a+ptp×a×j
a＝csize2×acc
式中：tp为降雨过程时中的时段编号，从1到tp，tp为降雨过程时段总数；etp为tp时段内栅格单元中的水流动能指标；ptp为tp时段内栅格单元中的降雨量；a为栅格单元中的累计汇流面积；j为栅格单元内的坡度；csize为栅格单元的边长；acc为栅格单元中的汇流累计值；
步骤2.2，在tp时段内的水文响应单元中，基于水流动能指标并结合流域实际的自然地理情况确定形成沟道的水流动能指标阈值te；
式中：bum为汇流到固定河源点points的区域中的栅格单元的数目；s为汇流到固定河源点points的区域中的栅格单元的编号，从1到bum；js为编号为s的栅格单元中的坡度；aum为流域中的所有栅格单元的数目；g为流域中的所有栅格单元的编号，从1到aum；jg为编号为g的栅格单元中的坡度；t为步骤1.2中确定的汇流阈值。
步骤2.3，对于tp时段内水文响应单元中的水流动能指标进行汇流累计计算，得到水流动能指标汇流累计栅格；
在水文响应单元中以栅格单元wcell为中心，通过周围栅格单元的高程值与该栅格单元的高程值的对比，找出与其相比最低的栅格单元wcelld；以wcell作为出流栅格单元，wcelld作为入流栅格单元，将wcelld中的水流动能指标加上wcell中的水流动能指标作为wcelld中的水流动能指标汇流累计值，逐栅格循环，计算水文响应单元中的每一个栅格单元的水流动能指标汇流累计值we，得到tp时段内的水流动能指标汇流累计栅格raster_we；
步骤2.4，利用te对步骤2.3中计算得到的raster_we中的栅格单元进行重分类，raster_we中的栅格单元的we高于te的判定为沟道栅格单元，we低于te的判定为坡面栅格单元。
前述的一种基于动态沟道的水文响应单元出流过程的计算方法，其特征在于，所述步骤3中估算水文响应单元中每一个栅格单元中的汇流流速，具体包括以下步骤：
步骤3.1，估算水文响应单元内每一个坡面栅格单元中的汇流速度；
vp＝k×q0.5×jp0.125
式中：vp为坡面栅格中的汇流速度；q为坡面栅格中的单宽流量，在降雨过程中以tp时段内的降雨量为输入，通过green—ampt方法计算得到；jp为坡面栅格中的坡度；k为经验系数，一般结合流域实际情况取值；
步骤3.2，估算水文响应单元内每一个沟道栅格单元中的汇流速度；
vr＝a×aε×jrβ
式中：vr为沟道栅格单元中的汇流速度；a为沟道栅格单元中的累计汇流面积，jr为沟道栅格单元中的坡度；a、ε以及β为经验系数，结合流域实际情况取值。
前述的一种基于动态沟道的水文响应单元出流过程的计算方法，其特征在于，所述步骤4中在水文响应单元中计算从每一个栅格单元出发，沿着一定的路径运动到水文响应单元出口所需要的时间，即该栅格单元的汇流时间，从而得到汇流时间栅格，具体包括以下步骤：
步骤4.1，结合水文响应单元中每一个栅格单元的边长以及该栅格单元的流向计算水流在该栅格单元中运动时所经过的路径的长度clenth；
若栅格单元中的流向为1、3、5、7，则：
clenth＝csize
若栅格单元中的流向为2、4、6、8，则：
步骤4.2，结合每一个栅格单元中的clenth以及汇流速度计算水流经过该栅格单元所需要的时间ctime；
在坡面栅格中：
在沟道栅格中：
步骤4.3，依据raster_dir中每一个栅格单元内的流向，将水流从栅格单元pcell出发，运动到水文响应单元的出口时需要经过的栅格单元依次连接起来得到该栅格单元的汇流运动路径path；并将水流经过path中各个栅格单元所需要的时间累加得到水流从pcell出发运动到水文响应单元出口所需要的时间，即汇流时间htime：
式中：n为汇流运动路径path中栅格单元的总数；i为汇流运动路径path中栅格单元的编号，从1到n；ctimei为水流经过path中编号为i的栅格单元所需要的时间；
步骤4.4，按照步骤4.3中的方法遍历水文响应单元中的每一个栅格单元，计算出每一个栅格单元的汇流时间，得到汇流时间栅格raster_htime。
前述的一种基于动态沟道的水文响应单元出流过程的计算方法，其特征在于，所述步骤5中基于汇流时间栅格，提取水文响应单元中的汇流单位线，具体包括以下步骤：
步骤5.1，基于步骤4.4中计算得到的汇流时间栅格统计得到水文响应单元中的栅格数—时间关系；
在raster_htime的基础上，以δtime为间隔进行分组统计，可以得到汇流时间落在每一组的时间范围内的栅格单元的数量，从而得到水文响应单元中的栅格数—时间关系；
式中：tcountz表示汇流时间大于等于z×δtime、小于(z+1)×δtime的栅格单元的数量；z为汇流时间分组的编号，从0到sz，其中sz＝max(τc)/δtime；δtime为时间间隔；τc为编号为c的栅格单元的汇流时间；max(τc)为流域水文响应单元中各个栅格单元的汇流时间中的最大值；m为水文响应单元中栅格单元的总数；c为水文响应单元中栅格单元的编号，从1到m；l为指示函数，当τc∈[z×δtime,(z+1)×δtime)时l为1，否则l为0；
步骤5.2，将步骤5.1中计算得到的水文响应单元中的栅格数—时间关系转换为面积—时间关系；
将汇流时间在[z×δtime,(z+1)×δtime)范围的栅格单元的数量乘以栅格单元的面积得到空间范围tareaz，即表示落在tareaz范围内的降雨会在[z×δtime,(z+1)×δtime)时间范围运动至水文响应单元的出口处；
tareaz＝tcountz×csize2
即得到水文响应单元中的面积—时间：
步骤5.3，将步骤5.2中得到的面积—时间关系转换成流量—时间关系；
式中：unitr为单位降雨深，一般取1mm；tflowz为在tp时段内tareaz范围内的单位降雨于[z×δtime,(z+1)×δtime)时间范围在水文响应单元出口处形成的径流，
即得到水文响应单元内的流量—时间关系：
步骤5.4，在流量—时间关系的基础上结合线性水库调蓄演算公式得到水文响应单元的汇流单位线；
uhtu＝y×tflowtu+(1-y)×uhtu-1
式中：tu为单位线时间序列编号，由1到uhtu小于0.01时的时刻tu；uhtu为tu时的单位线纵坐标值；uhtu-1为tu-1时单位线纵坐标值；tflowtu的值由tu与0、sz以及tu的相对大小决定，
y为调蓄系数，由以下公式计算：
式中：f为线性水库的演算系数，一般取sz×δtime的0.5倍。
前述的一种基于动态沟道的水文响应单元出流过程的计算方法，其特征在于，所述步骤6中对于降雨过程中的每一个时段重复步骤2到步骤5之间的操作，得到每一个时段中水文响应单元的汇流单位线，并结合降雨过程进行卷积计算，得到水文响应单元的出流过程，具体包括以下步骤：
步骤6.1，基于每一个时段的降雨，计算水文响应单元中逐时段的平均径流深；
式中：c为水文响应单元中栅格单元的编号，从1到m；m为水文响应单元中栅格单元的总数；rc,tp为以tp时段中编号为c的栅格单元中的降雨为输入，通过green—ampt方法计算得到的净雨深；rave,tp为tp时段中水文响应单元的平均净雨深；
步骤6.2，对水文响应单元中逐时段的汇流单位线和平均净雨深进行卷积计算，得到出口处的径流过程：
式中：t为直接径流量时序，从1到tp+tu-1；utp(δtime,t-tp)表示与tp时段的雨量相对应的，以t-tp为时间起点、以δtime为时间间隔的汇流单位线，该单位线是将在tp时段的降雨条件下形成的步骤5.4中的单位线将其开始的时间向后推移t-tp得到的。
本发明的有益效果：本发明提供的一种基于动态沟道的水文响应单元出流过程的计算方法，以影响沟道形成的物理因子为基础，刻画了降雨过程中沟道在水文响应单元中的动态发展变化，量化了沟道对于汇流的影响作用，进而提取了水文响应单元中的汇流单位线，并结合时段降雨量与该时段的汇流单位线进行卷积计算得到了水文响应单元出口处的出口过程。这样既保证了计算结果的精度与可靠性，同时解决了缺少观测资料的中小流域中水文响应单元的坡面汇流的计算问题。且本方法主要应用流域数字高程模型，数据来源稳定可靠，方法中变量之间的函数关系明确，有利于水文响应单元中沟道的自动判别与汇流单位线的自动化生成，通过数字流域技术以简化提取步骤，同时，保证了结果的客观合理性，有利于分布式水文模型的直接调用，可以进一步促进数字水文学以及分布式模型的深入发展。
附图说明
图1计算流程示意图；
图2水文响应单元示意图；
图3沟道栅格单元与坡面栅格单元示意图；
图4水文响应单元内汇流流速分布示意图；
图5水文响应单元内水流经过栅格单元时的路径长度分布示意图；
图6水文响应单元内水流经过栅格单元所需要的时间分布示意图；
图7水文响应单元内计算出栅格单元的汇流时间分布示意图；
图8水文响应单元中的栅格数—时间关系；
图9水文响应单元中的面积—时间关系；
图10水文响应单元中的流量—时间关系；
图11水文响应单元的汇流单位线示意图；
图12水文响应单元的出流过程示意图。
具体施方式
下面结合附图和具体实施例对本发明作进一步描述。
如图1所示，本发明提供的一种基于动态沟道的水文响应单元出流过程的计算方法，包括以下步骤：
步骤1，基于dem数据提取水文响应单元，具体包括以下步骤：
步骤1.1，利用流域dem数据计算出流域中每一个栅格单元的坡度、流向和汇流累计值，得到坡度栅格、流向栅格和汇流累计栅格；
以栅格单元cell为中心，通过周围栅格单元的高程值与该栅格单元的高程值的对比，找出与其相比最低的栅格单元celld，并计算cell和celld之间的高程差dhmax和水平投影距离dis，结合dhmax和dis计算栅格单元cell的坡度j：
j＝dhmax/dis
将cell作为出流栅格单元，celld作为入流栅格单元，入流栅格单元汇流累计值加1，逐栅格循环，计算出每一个栅格单元中的汇流累计值acc，同时依据cell和celld之间的相对位置关系，采用d8流向法确定cell中的流向，按照以上方法遍历流域中每一个栅格单元，从而得到坡度栅格raster_slope、流向栅格raster_dir和汇流累计栅格raster_acc；
步骤1.2，基于面积汇流阈值提取流域中的河道栅格，并结合流域资料和卫星影像确定固定河源点；
首先初步假定值t，将raster_acc中acc高于t的栅格单元判定为河道栅格单元，并将提取得到的河道栅格单元的位置与googleearth上的卫星影像中的实际河道相对比，反复调整t，使得提取出来的河道栅格单元的位置与卫星影像中的实际河道的位置重合，从而去确定汇流阈值t，同时提取出流域中的河道栅格，并进一步得到各个河道的固定河源点的集合，即固定河源点的集合points；
步骤1.3，以points中的一个固定河源点point为出口点，结合raster_dir，提取出所有汇流到point的栅格单元，并将这些栅格单元所在的区域作为point所对应的水文响应单元subaisn，遍历points中的所有固定河源点，从而提取出流域中所有的水文响应单元。为简化说明，本文中步骤2至步骤6中均以一个水文响应单元subaisn为例。
步骤2，基于降雨过程中的时段降雨量计算水文响应单元中各个栅格单元的水流动能指标，并由此判别该时段中的坡面栅格和沟道栅格，具体包括以下步骤：
步骤2.1，计算水文响应单元中该时段内每一个栅格单元内的水流动能指标；
etp＝a+ptp×a×j
a＝csize2×acc
式中：tp为降雨过程时中的时段编号，从1到tp，tp为降雨过程时段总数；etp为tp时段内栅格单元中的水流动能指标；ptp为tp时段内栅格单元中的降雨量；a为栅格单元中的累计汇流面积；j为栅格单元内的坡度；csize为栅格单元的边长；acc为栅格单元中的汇流累计值；
步骤2.2，在tp时段内的水文响应单元中，基于水流动能指标并结合流域实际的自然地理情况确定形成沟道的水流动能指标阈值te；
式中：bum为汇流到固定河源点points的区域中的栅格单元的数目；s为汇流到固定河源点points的区域中的栅格单元的编号，从1到bum；js为编号为s的栅格单元中的坡度；aum为流域中的所有栅格单元的数目；g为流域中的所有栅格单元的编号，从1到aum；jg为编号为g的栅格单元中的坡度；t为步骤1.2中确定的汇流阈值。
步骤2.3，对于tp时段内水文响应单元中的水流动能指标进行汇流累计计算，得到水流动能指标汇流累计栅格；
在水文响应单元中以栅格单元wcell为中心，通过周围栅格单元的高程值与该栅格单元的高程值的对比，找出与其相比最低的栅格单元wcelld；以wcell作为出流栅格单元，wcelld作为入流栅格单元，将wcelld中的水流动能指标加上wcell中的水流动能指标作为wcelld中的水流动能指标汇流累计值，逐栅格循环，计算水文响应单元中的每一个栅格单元的水流动能指标汇流累计值we，得到tp时段内的水流动能指标汇流累计栅格raster_we；
步骤2.4，利用te对步骤2.3中计算得到的raster_we中的栅格单元进行重分类，raster_we中的栅格单元的we高于te的判定为沟道栅格单元，we低于te的判定为坡面栅格单元。
步骤3，估算水文响应单元中每一个栅格单元中的汇流流速，具体包括以下步骤：
步骤3.1，估算水文响应单元内每一个坡面栅格单元中的汇流速度；
vp＝k×q0.5×jp0.125
式中：vp为坡面栅格中的汇流速度；q为坡面栅格中的单宽流量，在降雨过程中以tp时段内的降雨量为输入，通过green—ampt方法计算得到；jp为坡面栅格中的坡度；k为经验系数，一般结合流域实际情况取值；
步骤3.2，估算水文响应单元内每一个沟道栅格单元中的汇流速度；
vr＝a×aε×jrβ
式中：vr为沟道栅格单元中的汇流速度；a为沟道栅格单元中的累计汇流面积，jr为沟道栅格单元中的坡度；a、ε以及β为经验系数，结合流域实际情况取值。
步骤4，在水文响应单元中计算从每一个栅格单元出发，沿着一定的路径运动到水文响应单元出口所需要的时间，即该栅格单元的汇流时间，从而得到汇流时间栅格，具体包括以下步骤：
步骤4.1，结合水文响应单元中每一个栅格单元的边长以及该栅格单元的流向计算水流在该栅格单元中运动时所经过的路径的长度clenth；
若栅格单元中的流向为1、3、5、7，则：
clenth＝csize
若栅格单元中的流向为2、4、6、8，则：
步骤4.2，结合每一个栅格单元中的clenth以及汇流速度计算水流经过该栅格单元所需要的时间ctime；
在坡面栅格中：
在沟道栅格中：
步骤4.3，依据raster_dir中每一个栅格单元内的流向，将水流从栅格单元pcell出发，运动到水文响应单元的出口时需要经过的栅格单元依次连接起来得到该栅格单元的汇流运动路径path；并将水流经过path中各个栅格单元所需要的时间累加得到水流从pcell出发运动到水文响应单元出口所需要的时间，即汇流时间htime：
式中：n为汇流运动路径path中栅格单元的总数；i为汇流运动路径path中栅格单元的编号，从1到n；ctimei为水流经过path中编号为i的栅格单元所需要的时间；
步骤4.4，按照步骤4.3中的方法遍历水文响应单元中的每一个栅格单元，计算出每一个栅格单元的汇流时间，得到汇流时间栅格raster_htime。
步骤5，基于汇流时间栅格，提取水文响应单元中的汇流单位线，具体包括以下步骤：
步骤5.1，基于步骤4.4中计算得到的汇流时间栅格统计得到水文响应单元中的栅格数—时间关系；
在raster_htime的基础上，以δtime为间隔进行分组统计，可以得到汇流时间落在每一组的时间范围内的栅格单元的数量，从而得到水文响应单元中的栅格数—时间关系；
式中：tcountz表示汇流时间大于等于z×δtime、小于(z+1)×δtime的栅格单元的数量；z为汇流时间分组的编号，从0到sz，其中sz＝max(τc)/δtime；δtime为时间间隔；τc为编号为c的栅格单元的汇流时间；max(τc)为流域水文响应单元中各个栅格单元的汇流时间中的最大值；m为水文响应单元中栅格单元的总数；c为水文响应单元中栅格单元的编号，从1到m；l为指示函数，当τc∈[z×δtime,(z+1)×δtime)时l为1，否则l为0；
步骤5.2，将步骤5.1中计算得到的水文响应单元中的栅格数—时间关系转换为面积—时间关系；
将汇流时间在[z×δtime,(z+1)×δtime)范围的栅格单元的数量乘以栅格单元的面积得到空间范围tareaz，即表示落在tareaz范围内的降雨会在[z×δtime,(z+1)×δtime)时间范围运动至水文响应单元的出口处；
tareaz＝tcountz×csize2
即得到水文响应单元中的面积—时间：
步骤5.3，将步骤5.2中得到的面积—时间关系转换成流量—时间关系；
式中：unitr为单位降雨深，一般取1mm；tflowz为在tp时段内tareaz范围内的单位降雨于[z×δtime,(z+1)×δtime)时间范围在水文响应单元出口处形成的径流，
即得到水文响应单元内的流量—时间关系：
步骤5.4，在流量—时间关系的基础上结合线性水库调蓄演算公式得到水文响应单元的汇流单位线；
uhtu＝y×tflowtu+(1-y)×uhtu-1
式中：tu为单位线时间序列编号，由1到uhtu小于0.01时的时刻tu；uhtu为tu时的单位线纵坐标值；uhtu-1为tu-1时单位线纵坐标值；tflowtu的值由tu与0、sz以及tu的相对大小决定，
y为调蓄系数，由以下公式计算：
式中：f为线性水库的演算系数，一般取sz×δtime的0.5倍。
步骤6，对于降雨过程中的每一个时段重复步骤2到步骤5之间的操作，得到每一个时段中水文响应单元的汇流单位线，并结合降雨过程进行卷积计算，得到水文响应单元的出流过程，具体包括以下步骤：
步骤6.1，基于每一个时段的降雨，计算水文响应单元中逐时段的平均径流深；
式中：c为水文响应单元中栅格单元的编号，从1到m；m为水文响应单元中栅格单元的总数；rc,tp为以tp时段中编号为c的栅格单元中的降雨为输入，通过green—ampt方法计算得到的净雨深；rave,tp为tp时段中水文响应单元的平均净雨深；
步骤6.2，对水文响应单元中逐时段的汇流单位线和平均净雨深进行卷积计算，得到出口处的径流过程：
式中：t为直接径流量时序，从1到tp+tu-1；utp(δtime,t-tp)表示与tp时段的雨量相对应的，以t-tp为时间起点、以δtime为时间间隔的汇流单位线，该单位线是将在tp时段的降雨条件下形成的步骤5.4中的单位线将其开始的时间向后推移t-tp得到的。
以陕西省大河坝流域为例，研究区dem原始数据采用美国太空总署(nasa)与国防部国家测绘局(nima)联合提供的90m分辨率srtm(shuttleradartopographymission)数据。
步骤1、基于dem数据提取水文响应单元：
1)利用流域dem数据计算出流域中每一个栅格单元的坡度、流向和汇流累计值，得到坡度栅格、流向栅格和汇流累计栅格；
以栅格单元cell为中心，通过周围栅格单元的高程值与该栅格单元的高程值的对比，找出与其相比最低的栅格单元celld，并计算cell和celld之间的高程差dhmax和水平投影距离dis，结合dhmax和dis计算栅格单元cell的坡度j：
j＝dhmax/dis
将cell作为出流栅格单元，celld作为入流栅格单元，入流栅格单元汇流累计值加1，逐栅格循环，计算出每一个栅格单元中的汇流累计值acc。同时依据cell和celld之间的相对位置关系，采用d8流向法确定cell中的流向。按照以上方法遍历流域中每一个栅格单元，从而得到坡度栅格raster_slope、流向栅格raster_dir和汇流累计栅格raster_acc；
2)基于面积汇流阈值提取流域中的河道栅格，并结合流域资料和卫星影像确定固定河源点；
首先初步假定汇流阈值t，将raster_acc中acc高于t的栅格单元判定为河道栅格单元，并将提取得到的河道栅格单元的位置与googleearth上的卫星影像中的实际河道相对比，反复调整t，使得提取出来的河道栅格单元的位置与卫星影像中的实际河道的位置重合，从而去确定汇流阈值t，在陕西省大河坝流域中，当采用90m分辨率的dem数据时t的值为500，即认为当有500个栅格单元的水流汇流到一起时则认为形成河道。同时提取出流域中的河道栅格，并进一步得到各个河道的固定河源点的集合，即固定河源点的集合points；
3)以points中的一个固定河源点point为出口点，结合raster_dir，提取出所有汇流到point的栅格单元，并将这些栅格单元所在的区域作为point所对应的水文响应单元subaisn，遍历points中的所有固定河源点，从而可以提取出流域中所有的水文响应单元。为简化说明，本文中步骤2至步骤6中均以一个水文响应单元subaisn为例，如图2所示。
步骤2、基于降雨过程中的时段降雨量计算水文响应单元中各个栅格单元的水流动能指标，并由此判别该时段中的坡面栅格和沟道栅格，具体包括以下步骤：
1)结合时段降雨量，汇流面积以及坡度计算水文响应单元中该时段内每一个栅格单元内的水流动能指标；
etp＝a+ptp×a×j
a＝csize2×acc
式中：tp为降雨过程时中的时段编号，从1到tp，tp为降雨过程时段总数；etp为tp时段内栅格单元中的水流动能指标；ptp为tp时段内栅格单元中的降雨量；a为栅格单元中的累计汇流面积；j为栅格单元内的坡度；csize为栅格单元的边长；acc为栅格单元中的汇流累计值；
2)在tp时段内的水文响应单元中，基于水流动能指标并结合流域实际的自然地理情况确定形成沟道的水流动能指标阈值te；
式中：bum为汇流到固定河源点points的区域中的栅格单元的数目；s为汇流到固定河源点points的区域中的栅格单元的编号，从1到bum；js为编号为s的栅格单元中的坡度；aum为流域中的所有栅格单元的数目；g为流域中的所有栅格单元的编号，从1到aum；jg为编号为g的栅格单元中的坡度；t为步骤1中确定的汇流阈值。由此计算得到的te的值为550。
3)对于tp时段内水文响应单元中的水流动能指标进行汇流累计计算，得到水流动能指标汇流累计栅格；
在水文响应单元中以栅格单元wcell为中心，通过周围栅格单元的高程值与该栅格单元的高程值的对比，找出与其相比最低的栅格单元wcelld。以wcell作为出流栅格单元，wcelld作为入流栅格单元，将wcelld中的水流动能指标加上wcell中的水流动能指标作为wcelld中的水流动能指标汇流累计值，逐栅格循环，计算水文响应单元中的每一个栅格单元的水流动能指标汇流累计值we，得到tp时段内的水流动能指标汇流累计栅格raster_we；
4)利用te对步骤2.3中计算得到的raster_we中的栅格单元进行重分类，raster_we中的栅格单元的we高于te的判定为沟道栅格单元，we低于te的判定为坡面栅格单元，如图3所示。
步骤3、估算水文响应单元中每一个栅格单元中的汇流流速，如图4所示，具体包括以下步骤：
1)估算水文响应单元内每一个坡面栅格单元中的汇流速度；
vp＝k×q0.5×jp0.125
式中：vp为坡面栅格中的汇流速度；q为坡面栅格中的单宽流量，在降雨过程中以tp时段内的降雨量为输入，通过green—ampt方法计算得到；jp为坡面栅格中的坡度；k为经验系数，一般结合流域实际情况取值，在本例中k取1.25；
2)估算水文响应单元内每一个沟道栅格单元中的汇流速度；
vr＝a×aε×jrβ
式中：vr为沟道栅格单元中的汇流速度；a为沟道栅格单元中的累计汇流面积，jr为沟道栅格单元中的坡度；a，ε以及β为经验系数，一般结合流域实际情况取值，在本例中a取2.25，ε取0.67，β取0.33。
步骤4、在水文响应单元中计算从每一个栅格单元出发，沿着一定的路径运动到水文响应单元出口所需要的时间，即该栅格单元的汇流时间，从而得到汇流时间栅格，具体包括以下步骤：
1)结合水文响应单元中每一个栅格单元的边长以及该栅格单元的流向计算水流在该栅格单元中运动时所经过的路径的长度clenth，如图5所示；
若栅格单元中的流向为1、3、5、7，则：
clenth＝csize
在本例中，此时的clenth为90m。
若栅格单元中的流向为2、4、6、8，则：
在本例中，此时的clenth为127m。
2)结合每一个栅格单元中的clenth以及汇流速度计算水流经过该栅格单元所需要的时间ctime，如图6所示；
在坡面栅格中：
在沟道栅格中：
3)依据raster_dir中每一个栅格单元内的流向，将水流从栅格单元pcell出发，运动到水文响应单元的出口时需要经过的栅格单元依次连接起来得到该栅格单元的汇流运动路径path。并将水流经过path中各个栅格单元所需要的时间累加得到水流从pcell出发运动到水文响应单元出口所需要的时间，即汇流时间htime：
式中：n为汇流运动路径path中栅格单元的总数；i为汇流运动路径path中栅格单元的编号，从1到n；ctimei为水流经过path中编号为i的栅格单元所需要的时间；
4)按照4.3中的方法遍历水文响应单元中的每一个栅格单元，计算出每一个栅格单元的汇流时间，得到汇流时间栅格raster_htime，如图7所示。
步骤5、基于汇流时间栅格，提取水文响应单元中的汇流单位线，具体包括以下步骤：
1)基于步骤4.4中计算得到的汇流时间栅格统计得到水文响应单元中的栅格数—时间关系；
在raster_htime的基础上，以δtime为间隔进行分组统计，可以得到汇流时间落在每一组的时间范围内的栅格单元的数量，从而得到水文响应单元中的栅格数—时间关系，如图8所示；
式中：tcountz表示汇流时间大于等于z×δtime，小于(z+1)×δtime的栅格单元的数量；z为汇流时间分组的编号，从0到sz，其中sz＝max(τc)/δtime，本例中sz＝5，从tcount0到tcount5分别为【9，30，35，52，32，7】；δtime为时间间隔，本例中δtime＝0.01h；τc为编号为c的栅格单元的汇流时间；max(τc)为流域水文响应单元中各个栅格单元的汇流时间中的最大值；m为水文响应单元中栅格单元的总数，本例中为165；c为水文响应单元中栅格单元的编号，从1到m；l为指示函数，当τc∈[z×δtime,(z+1)×δtime)时l为1，否则l为0；
2)将步骤5.1中计算得到的水文响应单元中的栅格数—时间关系转换为面积—时间关系，如图9所示；
将汇流时间在[z×δtime,(z+1)×δtime)范围的栅格单元的数量乘以栅格单元的面积得到空间范围tareaz，即表示落在tareaz范围内的降雨会在[z×δtime,(z+1)×δtime)时间范围运动至水文响应单元的出口处；
tareaz＝tcountz×csize2
即得到水文响应单元中的面积—时间：
本例中计算得到的tarea0到tarea5分别为【0.0729，0.243，0.2835，0.4212，0.2592，0.0567】，单位为km2。
3)将步骤5.2中得到的面积—时间关系转换成流量—时间关系，图10所示；
式中：unitr为单位降雨深，一般取1mm；tflowz为在tp时段内tareaz范围内的单位降雨于[z×δtime,(z+1)×δtime)时间范围在水文响应单元出口处形成的径流；
即得到水文响应单元内的流量—时间关系：
在本例中计算得到的tflow0到tflow5的值分别为【0.121，0.405，0.472，0.702，0.432，0.094】单位为m3/s。
4)在流量—时间关系的基础上结合线性水库调蓄演算公式得到水文响应单元的汇流单位线；
uhtu＝y×tflowtu+(1-y)×uhtu-1
式中：tu为单位线时间序列编号，由1到uhtu小于0.01时的时刻tu；uhtu为tu时的单位线纵坐标值；uhtu-1为tu-1时单位线纵坐标值；tflowtu的值由tu与0，sz以及tu的相对大小决定；
y为调蓄系数，由以下公式计算：
式中：f为线性水库的演算系数，一般取sz×δtime的0.5倍，在本例中f为0.025，y为0.33。
本例中计算得到的单位线为【0，0.040，0.162，0.265，0.411，0.418，0.310，0.206，0.137，0.091，0.061，0.040，0.027，0.018，0.012，0.008，0.005，0.003，0.002，0.001，0】，如图11所示。
步骤6、对于降雨过程中的每一个时段重复步骤2到步骤5之间的操作，得到每一个时段中水文响应单元的汇流单位线，并结合降雨过程进行卷积计算，得到水文响应单元的出流过程，具体包括以下步骤：
1)基于每一个时段的降雨，计算水文响应单元中逐时段的平均径流深；
式中：c为水文响应单元中栅格单元的编号，从1到m；m为水文响应单元中栅格单元的总数；rc,tp为以tp时段中编号为c的栅格单元中的降雨为输入，通过green—ampt方法计算得到的净雨深；rave,tp为tp时段中水文响应单元的平均净雨深；
本例中计算得到的水文响应单元的平均净雨深过程为【1,2,1】。
2)对水文响应单元中逐时段的汇流单位线和平均净雨深进行卷积计算，得到出口处的径流过程：
式中：t为直接径流量时序，从1到tp+tu-1；utp(δtime,t-tp)表示与tp时段的雨量相对应的，以t-tp为时间起点，以δtime为时间间隔的汇流单位线。该单位线是将在tp时段的降雨条件下形成的步骤5.4中的单位线将其开始的时间向后推移t-tp得到的。
本例中计算得到的水文响应单元出口处的径流过程为【0，0.023，0.159，0.487，0.837，1.318，1.437，1.269，0.996，0.719，0.500，0.348，0.243，0.169，0.118，0.083，0.058，0.040，0.028，0.020，0.014，0.009，0.003，0.002，0】，如图12所示。
另外需要说明的是，在上述具体实施方式中所描述的各个具体技术特征，在不矛盾的情况下，可以通过任何合适的方式进行组合。为了避免不必要的重复，本发明对各种可能的组合方式不再另行说明。
以上显示和描述了本发明的基本原理、主要特征及优点。本行业的技术人员应该了解，本发明不受上述实施例的限制，上述实施例和说明书中描述的只是说明本发明的原理，在不脱离本发明精神和范围的前提下，本发明还会有各种变化和改进，这些变化和改进都落入要求保护的本发明范围内。本发明要求保护范围由所附的权利要求书及其等效物界定。