本发明属于水利工程技术领域，尤其涉及山洪防治工程技术领域，具体为一种山洪数值模拟方法。
背景技术：
山洪是指山丘区小流域强降雨引起的溪沟洪水，具有突发、暴涨暴落的特性，常同时诱发山体滑坡和泥石流，往往给局部地区带来严重灾害，造成人员伤亡。山区小流域暴雨洪水计算技术是进行山洪灾害防治工作中的一项核心技术，目前，对山区小流域的暴雨洪水分析多是采用水文模型的方法，如HBV模型、分布式的新安江模型、KINEROS(KINEROS2)模型、GBHM模型、CASC2D模型和PIHM模型等。这些模型都能用于山洪的分析计算，但是在应用到山洪灾害防治工作中时又都存在着各自的不足，例如它们基本都不能灵活处理流域内存在的堰、闸、涵洞等特殊构筑物的影响；有的模型对水文资料依赖高，水文资料匮乏地区应用困难；有的模型仅能提供控制断面处的水力要素信息，无法提供面上的水力要素信息，在进行山洪危险等级划分时受到限制；有的模型采用简化的运动波或扩散波方程模拟山洪运动，水流复杂时会带来较大的计算误差。
技术实现要素：
本发明的目的在于提供一种基于Godunov格式一维、二维耦合技术的山洪过程数值模拟方法，以解决山区洪水在缺乏水文资料情况下，水文模型适用性不强的问题，同时可以充分发挥水动力模型的优势来解决山洪的一些复杂技术问题，该方法能够用于山区小流域的洪水风险分析和实时预报预警业务中。本发明的目的是通过以下技术方案实现的：一种基于Godunov格式一、二维耦合技术的山洪数值模拟方法，该方法采用完全的水动力学方法来模拟山洪过程，包括以下步骤：1)获取流域的数字高程模型数据、沟道断面数据、流域土地利用类型遥感影像数据和降雨输入数据；2)采用四边形非结构网格离散山区流域坡面单元，采用一维有限控制体单元离散沟道部分；3)采用显格式进行坡面二维计算和沟道一维计算，首先确定一、二维计算采用的时间步长Δt；4)采用完整二维浅水方程组描述山洪坡面流动，采用Roe格式显式计算网格界面处的数值通量，变量定义在网格单元型心；由t时刻的净雨强度值和t时刻流域坡面水力要素值，计算t时刻二维网格单元普通单元边界面处的数值通量值；5)采用堰流公式计算二维坡面交互单元边处在t时刻的流量值；6)采用一维浅水方程组描述沟道山洪运动，采用HLL格式显式计算控制体界面处的数值通量，变量定义在控制体单元型心；由t时刻的坡面入流边界条件和t时刻沟道的水力要素值，计算t时刻沟道一维控制体单元界面处的数值通量；7)根据t时刻二维网格单元普通边界面处的数值通量和t时刻交互单元边处的流量值，计算二维坡面网格单元在t+1时刻的水力要素值；8)根据t时刻一维控制体单元界面处的数值通量和t时刻交互单元边处的流量值，计算一维控制体单元在t+1时刻的水力要素值；9)获取t+1时刻的净雨强度值，重复步骤4)～8)，直到计算结束。进一步的，步骤4)二维浅水方程的连续方程为如式(1)：式中：h为水深，u，v分别为x，y方向的流速，q2r为净雨源项，q2c为坡面与沟道水流交互源项。进一步的，步骤5)中采用如式(2)的堰流公式计算二维坡面和一维沟道间的水流交互：式中，h上＝max(Z坡,Z沟)-Z连；h下＝min(Z坡,Z沟)-Z连；q为通过二维单元连接边处的单宽流量；Z坡为二维连接单元水位，Z沟为对应的一维沟道单元水位，Z连为连接边的高程值；g为重力加速度。进一步的，步骤6)一维浅水方程的连续方程为如式(3)：式中：B为断面宽度，Z为水位，Q为流量，q1r为净雨源项，q1c为沟道与坡面水流交互源项。有益效果：整个模拟采用基于完整一、二维浅水方程的水动力学方法来描述山洪的地表运动过程，其中，坡面流采用二维方法来模拟，沟道洪水过程采用一维方法来模拟，坡面和沟道间的水流交互通过堰流公式来描述，在数值求解方法上坡面和沟道部分均采用了Godunov格式来求解，该方法能够自动适应流态变化，捕捉激波，非常适合模拟山洪这种水面梯度大的水流运动。相比于水文学方法，水动力学方法具有一些独特的优势，1)能够更真实的反应小流域情况，在山区小流域内，常存在一些水工建筑物(如堰、闸、涵洞等)，另外，铁路、公路、旅游景区等基础设施的建设也比较常见，采用水动力学方法可以很好地考虑这些因素的影响。2)提供更详细的水力信息数据。山洪灾害不仅会发生在溪沟两侧的村庄，在流域内较低洼的地区也可能发生山洪灾害。采用水动力学模型能够方便地模拟出流域内所有的风险点，提供整个流域面上的水深和流速信息，更利于完成山洪风险分析工作。3)能够适用于水文资料缺乏地区。水动力学模型需要的参数相对较少，且具有明确物理意义，没有实测水文资料验证模型时，可以根据下垫面情况选择较合理的参数。Godunov格式不仅适用于光滑的古典解，同时适用于大梯度的水面流动模拟，能够自动捕捉激波和水面间断，基于Godunov格式的暴雨洪水水动力模型逐渐成为小流域降雨径流计算的一种新的技术手段。为了适应现有基础数据的资料精度，提高水动力模型模拟山洪过程的实用性，将坡面汇流过程和沟道汇流过程分开考虑是非常必要的，研究成果具有广泛的推广应用价值。附图说明图1是基于Godunov格式一维、二维耦合技术的山洪数值模拟方法计算流程图。具体实施方式下面结合附图1对本发明作进一步的描述。本发明提供的方法是一种基于Godunov格式来模拟山洪地表径流过程的水动力方法。用二维非结构网格离散流域坡面，采用二维水动力模型来计算坡面汇流过程；沟道部分采用一维控制体单元进行离散，采用一维水动力模型来计算沟道洪水的汇流过程；坡面和沟道部分采用堰流公式来进行二维坡面流和一维沟道流之间的水量交互过程。该方法的具体过程如下：1)获取流域的数字高程模型数据、沟道断面数据、流域土地利用类型遥感影像数据和降雨输入数据；；2)采用四边形非结构网格离散山区流域坡面单元，采用一维有限控制体单元离散沟道部分。在坡面与沟道交界面处，单个二维网格不能跨越两个一维控制体单元。四边形非结构单元边按功能分为普通单元边与交互单元边，与一维控制体单元相连接，有水量交互的二维单元边称为交互单元边，其余为普通单元边，交互单元边所在的单元称为连接单元，其余称为普通单元。3)坡面二维计算和沟道一维计算均采用显格式进行，因此，计算步长受到CFL(Courant-Friedrichs-Lewy)条件限制，在计算中一、二维计算采用相同的时间步长Δt。二维坡面流计算时CFL限制条件如下：其中，其中u，v为x，y方向的流速分量，h为水深，g为重力加速度，Ncfl为CFL数，推荐在计算中取为0.8，Δt2D为二维计算时间步长，dL,LR为二维控制单元中心到对应边中点之间的距离。一维沟道流计算时受到CFL限制条件如下：其中，V为沟道流速，c为波速，Δx为一维控制体单元长度，Δt1D为一维计算时间步长。整个计算过程中二维计算和一维计算采用相同的时间步长，采用的时间步长Δt的取值如下：Δt＝min(Δt2D,Δt1D)4)采用完整二维浅水方程组描述山洪坡面流动，采用Roe格式显式计算网格界面处的数值通量，变量定义在网格单元型心(CC式)。由t时刻的净雨强度值和t时刻流域坡面水力要素值，计算t时刻二维网格单元普通单元边界面处的数值通量值。该步骤中二维浅水方程的连续方程为如下形式：其中h为水深，u，v分别为x，y方向的流速，q2r为净雨源项。与普通的二维浅水方程的连续方程相比，在连续方程中引入了q2c，q2c为坡面与沟道水流交互源项，在二维连接网格单元(与沟道相连接的单元)中作源项处理。采用的完整的二维浅水方程来描述山洪坡面运动，采用的二维浅水方程的向量守恒形式如下：其中：h为水深，u，v分别为x，y方向的流速；分别为x，y方向的底坡项，其中Zb为底高程。分别为x，y方向的摩阻坡降，其中n为曼宁糙率系数，q2r为净雨源项，q2c为坡面与沟道水流交互源项。5)采用堰流公式计算二维坡面交互单元边处在t时刻的流量值。采用的堰流公式如下：其中，h上＝max(Z坡,Z沟)-Z连；h下＝min(Z坡,Z沟)-Z连；q为通过二维单元连接边处的单宽流量；Z坡为二维连接单元水位，Z沟为对应的一维沟道单元水位，Z连为连接边的高程值；g为重力加速度。6)采用一维浅水方程组描述沟道山洪运动，采用HLL格式显式计算控制体界面处的数值通量，变量定义在控制体单元型心(CC式)。由t时刻的坡面入流边界条件和t时刻沟道的水力要素值，计算t时刻沟道一维控制体单元界面处的数值通量。该步骤中一维浅水方程的连续方程为如下形式：其中B为断面宽度，Z为水位，Q为流量，q1r为净雨源项。为考虑沟道一维与坡面二维水流的交互，在一维浅水方程的连续方程中引入了q1c项，该项在一维模型中作集中入流处理。此处的q1c对应步骤4)中提到的二维浅水连续方程中的q2c项，从物理意义上讲，q1c与q2c描述的都是二维坡面与一维沟道间的水流交互项，只是量纲表述有所不同。该项具体的计算方法见步骤5)。采用的一维浅水方程的向量形式如下：其中式中B为水面宽度，Z为水位，Q为断面流量，A为过水断面面积，为便于随后的表述，f1和f2分别代表向量F(U)的两个分量，g为重力加速度，t为时间变量，J为沿程阻力损失，其表达式为J＝(n2Q|Q|)/(A2R4/3)，R为水力半径，n为Manning糙率系数，q1r为净雨源项，q1c为沟道一维与坡面二维水流的交互项。7)根据t时刻二维网格单元普通边界面处的数值通量和t时刻交互单元边处的流量值，计算二维坡面网格单元在t+1(t+Δt)时刻的水力要素值。8)根据t时刻一维控制体单元界面处的数值通量和t时刻交互单元边处的流量值，计算一维控制体单元在t+1(t+Δt)时刻的水力要素值。9)获取t+1时刻的净雨强度值，重复步骤4)～8)，直到计算结束。上述的实施例仅是本发明的部分体现，并不能涵盖本发明的全部，在上述实施例以及附图的基础上，本领域技术人员在不付出创造性劳动的前提下可获得更多的实施方式，因此这些不付出创造性劳动的前提下获得的实施方式均应包含在本发明的保护范围内。