技术领域
[0001] 本
发明属于水利工程技术领域，特别涉及一种利用完整二维浅水方程组获得流域单位线的方法。
背景技术
[0002] 所谓单位线是指单位时段内在流域上均匀分布的单位净雨量在流域出口断面形成的流量过程线，单位线法是流域汇流计算最常用的方法之一。传统业务应用上，常采用分析法和试错法来推求一个流域的单位线，这些方法的应用需要一个前提条件，那就是流域有充分的实测降雨和径流数据。但是，在我国大部分中小流域中，并无充分的降雨径流数据。
[0003] 当前，我国目前正在开展山洪灾害治理项目和中小河流治理项目，在项目建设过程中，有很多科学技术问题亟待解决，如山洪危险区的划分问题、山洪预警指标
阈值的确定问题、山洪
风险图成图技术、中小河流的预报预警问题等，这些问题的解决都离不开流域暴雨洪水分析计算技术。很多传统的
水文学模型在进行暴雨洪水分析计算时往往会用到单位线的方法，虽然在山洪项目和中小河流项目中建设了一大批雨量站和水位站，但这些
站点的主要作用是用于预报预警，项目中并没有监测流量过程数据，因此采用传统方法推求单位线依然困难，传统单位线法的应用受到限制。
发明内容
[0004] 本发明的目的是提供一种利用完整二维浅水方程组推求单位线的分析方法，该方法区别于传统的单位线获取方法，不需要大量的降雨径流数据，而是采用完整的二维浅水方程组构建数学模型直接模拟连续净雨过程，通过获取流域出口S曲线来得到流域的单位线过程。该方法在无资料地区依然有很好的适用性，克服了现在中小流域由于水文资料匮乏而导致的单位线汇流计
算法难以使用的问题。
[0005] 本发明的目的是通过以下方案实现的：
[0006] 本发明为一种利用完整二维浅水方程组推求单位线的方法，该方法充分利用水动
力学的优势解决一个水文学的问题，采用能够捕捉激波的高
分辨率的数值格式离散完整的二维浅水方程组来构建二维浅水数学模型，利用该模型模拟连续时段净雨在流域坡面上的运动过程，从而得到流域出口处的S曲线，根据S曲线进一步得到流域单位线，其具体方法步骤如下：
[0007] 1)采集流域高
精度地形数据和流域
地貌类型数据，流域地形分辨率不低于10m。
[0008] 2)完整二维浅水方程组的离散。采用能够捕捉激波的高分辨率的Roe格式离散完整的二维浅水方程组，建立二维浅水数学模型，确保该数学模型能够准确模拟流域坡面水流的运动过程。
[0009] 3)流域计算范围离散及参数赋值。采用三
角形非结构单元离散流域计算范围，
对流域内的挡水构筑物(如道路)和沟渠等部分进行适当加密。采用地形数据对网格型心和网格
节点进行高程插值，获取网格型心和单元节点处的高程数据。根据流域地貌类型数据对网格单元的糙率进行赋值，区分坡面和
沟道糙率差异。
[0010] 4)设置净雨输入条件。根据S曲线定义，设置一定时段(标记为DT)内的净雨量10mm，假定该净雨过程持续不断，以此作为二维浅水数学模型的输入条件。
[0011] 5)设置流域出口边界条件。在进行单位线推求时，流域出口处并无水力约束条件，因此设定流域出口边界为自由跌水出流。
[0012] 6)启动二维浅水数学模型计算。所有网格单元的初始流速和水深均设置为0；流域‐6干湿阈值的判别指标设置为10 m；模型为显格式开发，其数值
稳定性受到CFL(Courant‐Friedrichs‐Lewy)条件限制，为保证计算数值稳定，设置CFL<1.0；运行二维浅水数学模型。
[0013] 7)获取该流域时段为DT的单位线。当流域出口流量值达到稳态的F×i后(其中F为流域面积，i为降雨强度)，二维浅水数学模型计算结束，随后输出流域出口边界处的流量过程即为该流域时段为DT的S曲线，根据公式(2)q(DT,t)＝S(DT,t)-S(DT,t-DT)得到时段为DT的单位线，q即为单位线过程上t时刻对应的流量值。
[0014] 进一步的，步骤(2)中，首先对完整二维浅水方程组进行离散。完整二维浅水方程组的守恒形式如下：
[0015]
[0016] 其中：
[0017]
[0018] h为水深，u，v分别为x，y方向的流速，t为时间。
[0019] 分别为x，y方向的坡度，Zb为地面高程，g为重力
加速度。
[0020] 分别为x，y方向的摩阻项，其中n为Manning糙率系数，r为净雨源项。
[0021] 为适应任意复杂边界，选用三角形非结构网格作为基本离散单元，为解决坡面流模拟时遇到的小水深问题，在每个单元格内假定坡面流为水深均匀的薄
层流动。采用三角形网格三点的坐标值计算该单元x，y方向的坡度，采用高分辨率的Roe格式计算单元界面间数值通量，为保证格式和谐性，需对数值通量进行修正。为避免小水深下出现不符合物理意义的大流速值引起的数值格式不稳定问题，摩阻项采用半隐式格式处理。
[0022] 进一步的，步骤(4)中设置净雨输入条件，该处设置的净雨过程为一持续的净雨过程，流域出口对应的出流过程为S曲线，之所以不直接用单一时段净雨作为输入条件，主要是为避免流域低洼地带蓄水造成流域出口水量不守恒问题。另外，根据洪水量级的不同，可以灵活设置不同量级的连续时段净雨作为输入条件，得到不同雨强条件下对应的单位线组，尽可能减小流域采用单一雨强单位线带来的误差。
[0023] 进一步的，步骤(5)中设置流域出口边界条件，在进行单位线计算时，流域出口处并不受水力条件约束，因此采用自由跌水出流状态，即假定与流域出口单元相连接的虚单元水深为0，单元高程相等，单元出口界面处的数值通量利用Roe格式求解。
[0024] 进一步的，步骤(7)中获取该流域时段为DT的单位线，该单位线过程与水文方法得到的单位线过程在细节上有所不同，用水文方法获取的时段为DT的单位线，其过程线上的数值也是以DT为时段的，本方法获得的时段为DT的单位线，其过程线上的数值理论上可以是任意时段dt(一般远小于DT)，这样可以把整个单位线过程描述的更加细致。对流域汇流时间小于DT的流域，这一点尤为重要。
[0025] 相比于水文学方法，
水动力学方法的优势如下：1)更真实的反应流域情况。在中小流域内，常存在一些水工
建筑物(如堰、闸、涵洞等)，另外，
铁路、公路、旅游景区等
基础设施的建设也比较常见，采用水动力学方法可以很好地考虑这些因素的影响。2)能够适用于水文资料缺乏地区。水动力学模型需要的参数相对较少，且具有明确物理意义，没有实测水文资料验证模型时，可以根据下垫面情况选择较合理的参数。
[0026] 由于扩散波和运动
波形式简单，数值稳定性相对较好，因此在工程上得到了较广泛的应用。从理论角度讲，控制方程中的压力项和惯性项能够引起洪水波的衰减，因此，忽略这些项，在模型中采用扩散波和运动波方程计算时是受到一定限制的，尤其是对于地形复杂的情况下，当水流出现激波运动现象时，采用简化的运动波或扩散波计算方法会带来较大的计算误差，采用完整的二维浅水方程对坡面流进行模拟可避免以上问题，保证数值计算的精度，但是计算速度又限制了其广泛应用。基于此，本发明提出一种新的方法，采用完整二维浅水方程组来获取流域单位线，既利用了水动力模型的优势，又规避了其计算效率的问题。
[0027] 本发明的有益效果：
[0028] 本发明充分发挥二维浅水模型的优势，利用二维浅水模型直接模拟连续净雨过程，得到流域出口的S曲线，进而得到流域的单位线过程。该方法克服了单位线过程的推求一般需要大量降雨径流资料的问题，使用本发明提供的方法，在无水文资料地区依然能够获得可靠的单位线过程，整个过程仅有糙率参数需要率定，而糙率值的选取在水利工程领域已相对比较成熟。该方法为单位线的推求提供了一种新的思路，能够扩大单位线汇流计算方法的适用范围，有效弥补现有水文方法的不足。
附图说明
[0029] 图1为本发明的一种利用完整二维浅水方程组推求时段单位线的方法
流程图；
[0030] 图2为某室内降雨径流试验装置三维效果图；
[0031] 图3为某室内试验降雨径流过程的模拟结果；
[0032] 图4为解家湾流域地形图；
[0033] 图5为解家湾流域S曲线(DT＝10min)；
[0034] 图6为解家湾流域时段单位线(DT＝10min)。
具体实施方式
[0035] 下面结合附图1和
实施例对本发明作进一步说明。
[0036] 本发明提供的是一种利用完整二维浅水方程组推求单位线的方法。完整的二维浅水方程组能够很准确的描述坡面径流的运动过程，这为该方法的应用打下了理论基础。二维浅水数学模型推求得到流域单位线，克服了在无水文资料地区难以获取单位线的问题，扩大了单位线汇流计算法的适用范围，为无水文资料地区的洪水汇流计算提供了一种新思路。该方法具体过程如下：
[0037] 1)该方法使用之前需要采集流域地形数据和土地地貌类型数据。流域地形数据可以从测绘部
门获取或者采用LiDAR等测绘技术去现场实测获取，土地地貌类型数据可以从专业
网站免费下载，为获得更高分辨率的土地地貌类型数据也可以到测绘部门获取或者无人机航拍获取。
[0038] 2)完整二维浅水方程组的离散。首先采用高分辨率的Roe格式离散完整的二维浅水方程组，建立二维浅水数学模型，确保该数学模型能够准确模拟流域坡面水流的运动过程。完整二维浅水方程组的守恒形式如下：
[0039]
[0040] 其中：
[0041]
[0042] h为水深，u，v分别为x，y方向的流速，t为时间。
[0043] 分别为x，y方向的坡度，Zb为地面高程，g为重力加速度。
[0044] 分别为x，y方向的摩阻项，其中n为Manning糙率系数，r为净雨源项。
[0045] 为适应任意复杂边界，选用三角形非结构网格作为基本离散单元。为解决坡面流模拟时遇到的小水深问题，在每个单元格内假定坡面流为水深均匀的薄层流动。采用三角形网格三点的坐标值计算该单元x，y方向的坡度，采用高分辨率的Roe格式计算单元界面间数值通量，为保证格式和谐性，需对数值通量进行修正，其修正方法可参见该文献描述(Kim J,Warnock A,Ivanov V Y et al.Coupled modeling of hydrologic and hydrodynamic processes including overland and channel flow.Advances in water resources,2012,37:104-126.)。为避免小水深下出现不符合物理意义的大流速值引起的数值格式不稳定问题，摩阻项采用半隐式格式处理，具体形式如下：
[0046]
[0047] 其中h为水深，u，v分别为x，y方向的流速，g为
重力加速度，n指n时刻，n+1指n+1时刻。
[0048] 3)流域计算范围离散及参数赋值。采用三角形非结构单元离散流域计算范围，对流域内的挡水构筑物(如道路)和沟渠等部分进行适当加密。采用地形数据对网格型心和网格节点进行高程插值，获取网格型心和单元节点处的高程数据。根据流域地貌类型数据对网格单元的糙率进行赋值，区分坡面和沟道糙率差异。
[0049] 4)设置净雨输入条件。根据S曲线定义，设置一定时段DT内的净雨量10mm，假定该净雨过程持续不断，以此作为二维浅水数学模型的输入条件。该处设置的净雨过程为一持续的净雨过程，流域出口对应的出流过程为S曲线，之所以不直接用单一时段净雨作为输入条件，主要是为避免流域低洼地带蓄水造成流域出口水量不守恒问题。另外，根据洪水量级的不同，可以灵活设置不同量级的连续时段净雨作为输入条件，得到不同雨强条件下对应的单位线组，尽可能减小流域采用单一雨强单位线带来的误差。
[0050] 5)设置流域出口边界条件。在进行单位线推求时，流域出口处并无水力约束条件，因此设定流域出口边界为自由跌水出流，即假定与流域出口单元相连接的虚单元水深为0，单元高程相等，单元出口界面处的数值通量利用Roe格式求解。该自由跌水出流边界与水动力学模型常用的水面梯度为0的自由出流边界是有明显不同的。
[0051] 6)启动二维浅水数学模型计算。所有网格单元的初始流速和水深均设置为0；流域干湿阈值的判别指标设置为10‐6m；模型为显格式开发，其数值稳定性受到CFL(Courant‐Friedrichs‐Lewy)条件限制，为保证计算数值稳定，设置CFL<1.0；运行二维浅水数学模型。
[0052] 7)获取该流域时段为DT的单位线。当流域出口流量值达到稳态的F×i后(其中F为流域面积，i为雨强)，二维浅水数学模型计算结束，随后输出流域出口边界处的流量过程即为该流域时段为DT的S曲线，根据公式(2)q(DT,t)＝S(DT,t)-S(DT,t-DT)得到时段为DT的单位线，q即为单位线过程上t时刻对应的流量值。该单位线过程与水文方法得到的单位线过程在细节上有所不同，用水文方法获取的时段为DT的单位线，其过程线上的数值也是以DT为时段的，本方法获得的时段为DT的单位线，其过程线上的数值理论上可以是任意时段dt(一般远小于DT)，这样可以把整个单位线过程描述的更加细致。对流域汇流时间小于DT的流域，这一点尤为重要。
[0053] 图2为某室内降雨径流试验装置三维效果图。流域装置由3片不锈
钢片组成，三
块钢片的坡度均为0.05，流域几何尺寸为2m×2.5m。在流域内部布置2堵墙来增加流域的汇流时间。采用喷淋系统模拟降雨过程，在流域出口测量出流过程用来检验数学模型的精度。降雨强度是317mm/h，降雨历时45s。图3为该场降雨径流的模拟结果，通过该结果可以看出，二维浅水数学模型的计算值与实测值吻合的非常好，说明本发明采用的二维浅水数学模型是可靠的。
[0054] 图4为解家湾流域地形图，该流域在四川省遂宁市安居镇境内。流域面积为65827m2，最高高程330m，最低高程275m，属于典型的丘陵地貌，流域平均坡度为0.27。假定时段DT＝10min，每个时段的净雨为10mm，净雨持续不断，则在流域出口形成的流量过程线即S曲线如图5所示。S曲线理论上的稳定值为1.097m3/s，本发明采用的二维浅水数学模型计算的数值解也为1.097m3/s，说明S曲线的计算是合理的。图6为采用该S曲线计算得到的单位线，计算的流域出口总水量为658.2m3，流域出口总水量的理论值为658.3m3，水量相对误差为0.02％，这说明该单位线是合理的，该发明提供的方法是可行的。从另一个角度可以看出本发明提供的单位线的优点，由于该流域面积较小，所给降雨时段为10min，如果采用传统单位线计算，则只能给出10min、20min和30min处有限的几个点，为一条折线，不能代表该流域真实的出流过程，而采用本发明提供的方法，同样降雨时段为10min，给出的是一条光滑的单位线，代表了流域真正的出流过程。
[0055] 上述的实施例仅是本发明的部分体现，并不能涵盖本发明的全部，在上述实施例以及附图的基础上，本领域技术人员在不付出创造性劳动的前提下可获得更多的实施方式，因此这些不付出创造性劳动的前提下获得的实施方式均应包含在本发明的保护范围内。
