技术领域
[0001] 本
发明属于
数据处理技术领域，尤其涉及一种基于GIS的大范围泥石流数值模拟及数值处理方法。
背景技术
[0002] 目前，最接近的
现有技术：基于
牛顿模型和Herschel-Bulkley模型为代表的复杂黏塑性模型进行求解的连续介质计算方法。其中牛顿模型没有考虑泥石流的黏性特征，无法真实地模拟泥石流的运动特点；复杂黏塑性模型通常有过多的变量，通常只能针对某一实例进行模拟，没有普遍适用性。本发明利用了有限差分的宾汉模型对泥石流进行模拟，计算消耗少，并能充分考虑泥石流的黏性特征。此外对于连续介质常用于小尺度范围的局限，本发明与GIS结合起来，巧妙利用DEM栅格网作为有限差分网格进行计算，可以处理大范围复杂地形上的泥石流运移。
[0003] 泥石流是山区常见的
自然灾害，每年在世界各地都造成大量的人员伤亡和财产损失。泥石流的流动距离可能很长，泛滥的范围可能很大，而且在非常平缓的斜坡上也具有较高的速度和传播能
力。泥石流灾害在防灾减灾工作中已成为一项亟待解决的问题。
[0004] 泥石流运动数值模拟方法可以分为连续介质、离散介质和混合介质三大类：
[0005] (1)连续介质计算方法，又分为基于网格和基于粒子进行离散化，进而求解计算的方法。其中基于网格的计算方法主要有：有限差分法、有限元法和有限体积法。基于网格是指将计算域划分为网格或者体元，控制方程可以转化为含有
节点变量信息的代数方程组。求微分方式的控制方程为有限差分法，求积分方式的控制方程为有限元法和有限体积法。
对于不同类型的泥石流常采用不同的流变本构关系。以黏性
应力为主的泥石流(通常拜格诺数小于40，且
雷诺数小于500)剪胀和湍动特征都不明显，可以基于宾汉模型、Cross模型、Herschel-Bulkley模型以及Voellmy模型等黏塑性本构模型进行求解；而
水石流(通常拜格诺数大于450，且相
对流深小于10)主要以固体颗粒间的相互作用(碰撞、摩擦)为主，剪胀效应明显，可以采用剪胀
体模型进行分析。
[0006] (2)离散介质计算方法，主要是离散元法和格子玻尔兹曼法。采用离散元法模拟三维颗粒流动，颗粒间
接触力考虑弹黏性和阻尼作用，
流体按欧拉法模拟，流固耦合考虑相间力，得到泥石流流动过程中的
应力分布和
能量变化，提出了DEM模拟泥石流的潜在优势。基于格子玻尔兹曼法，结合泥石流流团模型的特点建立特殊格式的格子玻尔兹曼(LB)模型，可以从平衡和非平衡态分布函数两个层次刻画泥石流复杂的流变关系和运移特征，尤其是泥石流的堆积过程。
[0007] (3)混合介质计算方法，是连续介质与离散元耦合的方法。包括混合的
计算流体动力学－离散元法(CFD-DEM)、大涡模拟与离散元法的耦合(LES-DEM)，还包括有限差分法、有限体积法、光滑粒子流体动力学、格子玻尔兹曼法与离散元法的耦合。耦合
算法涉及的变量多且算法复杂，当模型较大时，每一时步的离散元计算结果与其他算法间交换的数据量很大，对计算机资源的消耗很大。
[0008] 对于小尺度的泥石流问题，常采用配合模型试验的方法，常用离散元和连续介质耦合方法进行分析。由于这些模型应用了过多的流变参数，在实际应用中受到了限制。离散元的计算速度跟单元体网格划分有很多关系，过多的接触可使计算量变大，花费大量时间，在大范围的泥石流计算中不存在优势。
[0009] 综上所述，现有技术存在的问题是：(1)现有连续介质方法所采用的本构模型通常无法准确描述泥石流的运动过程或者难以描述所有类型泥石流的运动过程。
[0010] (2)离散元的计算速度跟单元体网格划分有很多关系，过多的接触可使计算量变大，花费大量时间，在大范围的泥石流计算中不存在优势。
[0011] (3)采用配合模型试验的方法，常用离散元和连续介质耦合方法应用了过多的流变参数，这些流变参数通常难以获得，而且耦合算法复杂，计算数据量很大，在实际应用中受到了限制。
[0012] 解决上述技术问题的难度：上述技术问题中，复杂参数的获取通常难以实现，针对不同形态的泥石流不能普遍使用，大范围计算造成的计算成本大幅增加，这些都是制约泥石流数值计算的难点。
[0013] 解决上述技术问题的意义：本试验提出的一种基于GIS的大范围泥石流数值模拟及数值处理方法，基于宾汉模型的有限差分计算可以实现对各种形态泥石流的模拟；所需变量仅为地形数据、泥石流
密度、粘聚力和摩擦
角和速度修正系数，不受过多经验参数准确性的影响；而且作为一种连续介质的方法算法成熟，在大范围泥石流数值模拟计算中可根据需要调整网格大小，计算消耗易于接受。综上，可以完美解决上述技术问题。
发明内容
[0014] 针对现有技术存在的问题，本发明提供了一种基于GIS的大范围泥石流数值模拟及数值处理方法。
[0015] 本发明是这样实现的，一种基于GIS的大范围泥石流数值模拟及数值处理方法，所述基于GIS的大范围泥石流数值模拟及数值处理方法包括以下步骤：
[0016] 第一步，通过现场观测获得地形数据，建立数值标高模型；
[0017] 第二步，采用有限差分法求解泥石流的连续方程和运动方程，利用地球信息系统GIS中的数字高程模型DEM，自动生成地形变量，利用数字高程模型中的栅格网作为有限差分的网格；利用野外试验或实际工程验证模型和进行参数调整；
[0018] 第三步，使用基于DEM栅格网的有限差分公式可以实现数值解。
[0019] 进一步，所述第一步泥石流是均匀连续、不可压缩、非定常的
牛顿流体，流动遵循下列方程，其连续性方程为：
[0020]
[0021] 运动方程为：
[0022]
[0023]
[0024]
[0025] 其中p为水压力，g为重力
加速度，u，v和w分别为x，y，z方向的平均流速，ρd是泥石流的平均密度；ρd＝ρsνs+ρwνw，其中ρs和ρw分别是固体颗粒和水的密度，νs和νw分别是混合物中固体颗粒和水的体积浓度；
[0026] 根据浅水流运动原理，u,v>>w，方程简化为：
[0027] p＝ρdg(η-ηb)；
[0028] 其中η是流体上表面的标高，ηb是流体底面的标高；
[0029] 流体底面的运动边界条件写为：
[0030]
[0031] 流体上表面的边界条件写为：
[0032]
[0033] 得到泥石流的连续方程为：
[0034]
[0035] 其中h为泥石流流深， 和 分别为x和y方向的平均流速；
[0036] 对方程 进行深度上的积分，得到：
[0037]
[0038]
[0039] 最终形式的深度积分运动方程为：
[0040]
[0041]
[0042] 其中α为速度修正系数，提供了有关垂直速度分布与均匀性偏差的信息，β为垂直正应力的比值，H为自由面的高度，ν＝μ/ρ为运动
粘度系数，μ是粘度，θx和θy分别是地形沿x和y方向的倾斜角度， 为动
摩擦系数，t是流动时间。
[0043] 进一步，所述第二步在基于栅格的DEM分析中，每个像元都有左、右、上、下以及四个对角线八个流动方向。
[0044] 进一步，所述第三步使用基于DEM栅格网的有限差分公式可以实现数值解，控制方程的有限差分形式为：
[0045]
[0046] x方向的有限差分形式为：
[0047]
[0048] 宾汉启动条件：
[0049]
宾汉流体在没有达到一定切应力(τ0)之前不会由任何
剪切速率；描述宾汉流体需要两个参数：屈服应力τ0和粘度μ：
[0050]
[0051] 流体当达到：
[0052] τ>τ0；
[0053] 时开始流动，对于
层流，根据粘性流体的牛顿法则，水平方向的切应力：
[0054]
[0055] 对上述方程在深度方向上积分，得到：
[0056]
[0057] 该流动是与水平成θ角的平面上稳定均匀的宾汉流动；屈服应力为：
[0058] τ0＝ρgsinθ；
[0059] 对上述方程在深度方向上积分，得到：
[0060]
[0061] 其中h＝h(x,y,t)＝η-ηb是流体的深度；
[0062] 宾汉流体从静止开始运动，启动条件为：
[0063]
[0064] 对上述启动条件进行有限差分，得到：
[0065]
[0066] 本发明的另一目的在于提供一种实施所述基于GIS的大范围泥石流数值模拟及数值处理方法的基于GIS的大范围泥石流数值模拟及数值处理系统，所述基于GIS的大范围泥石流数值模拟及数值处理系统包括：
[0067] 参数调整模
块，用于通过现场观测获得地形数据，建立数值标高模型；可以利用野外试验或实际工程验证模型和进行参数调整；
[0068] 地形变量生成模块，用于采用有限差分法求解泥石流的连续方程和运动方程，利用地球信息系统GIS中的数字高程模型DEM，自动生成地形变量，利用数字高程模型中的栅格网作为有限差分的网格。
[0069] 数值解获取模块，用于使用基于DEM栅格网的有限差分公式可以实现数值解。
[0070] 本发明的另一目的在于提供一种实现所述基于GIS的大范围泥石流数值模拟及数值处理方法的信息数据处理终端。
[0071] 本发明的另一目的在于提供一种计算机可读存储介质，包括指令，当其在计算机上运行时，使得计算机执行所述的基于GIS的大范围泥石流数值模拟及数值处理方法。
[0072] 综上所述，本发明的优点及积极效果为：本发明在牛顿流体运动模型
基础上，推导宾汉流体的连续方程和运动方程，尤其是宾汉流体启动条件，并将其深度积分以适用于有限差分网格模拟计算。并将宾汉流体模型与GIS耦合，便于利用大范围复杂地形的计算。
[0073] 本发明提出了一个深度积分的数值模型用来模拟大范围泥石流的运动，可以在复杂的地形上得到泥石流的运动路径、运动速度和扩散范围。并采用了与GIS技术相结合的模拟方法，由于GIS中的栅格网络可以作为有限差分的网格，因此可以使用有限差分方法对连续性方程和运动方程进行数值求解。差分网格的大小可以根据实际工程需要任意调整，超大范围内多发生源的泥石流模拟可以使用大网格，以减少计算时间；单发生源的泥石流模拟可以调小网格，以准确判断灾害影响。所有的输入和输出数据都可以在GIS中处理，在可以实现
可视化的前提下，准确模拟大范围泥石流灾害。
[0074] 在大范围泥石流的流动过程和流经范围数值模拟方面的定量研究的突破，将为定量描述泥石流运动特征，提高泥石流预报
精度，精准做好泥石流防护和防治工程提供理论基础。而这些理论和方法的研究成果不仅能提升我国泥石流研究的理论水平，而且能为国家防灾减灾，保障山区人民生命财产安全，以及对交通、
水电等大型基础设施建设的安全提供保障。
附图说明
[0075] 图1是本发明
实施例提供的基于GIS的大范围泥石流数值模拟及数值处理方法
流程图。
[0076] 图2是本发明实施例提供的单个单元格的泥石流流动方向(周边单元格里的数字表示流动方向的角度)示意图。
[0077] 图3是本发明实施例提供的DEM网格示例，数字表示单元格标高，箭头表示泥石流流动方向示意图。
[0078] 图4是本发明实施例提供的泥石流运动的有限差分网格示意图。
[0079] 图5是本发明实施例提供的一个大范围泥石流运动模拟的实例图，在GIS中实现可视化效果示意图。
[0080] 图6是本发明实施例提供的某泥石流流经区域标高变化图。
[0081] 图7是本发明实施例提供的某泥石流模拟流经区域变化图。
具体实施方式
[0082] 为了使本发明的目的、技术方案及优点更加清楚明白，以下结合实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。
[0083] 现有连续介质方法所采用的本构模型通常无法准确描述泥石流的运动过程或者难以描述所有类型泥石流的运动过程。离散元的计算速度跟单元体网格划分有很多关系，过多的接触可使计算量变大，花费大量时间，在大范围的泥石流计算中不存在优势。采用配合模型试验的方法，常用离散元和连续介质耦合方法应用了过多的流变参数，这些流变参数通常难以获得，而且耦合算法复杂，计算数据量很大，在实际应用中受到了限制。
[0084] 针对现有技术存在的问题，本发明提供了一种基于GIS的大范围泥石流数值模拟及数值处理方法，下面结合附图对本发明作详细的描述。
[0085] 如图1所示，本发明实施例提供的基于GIS的大范围泥石流数值模拟及数值处理方法包括以下步骤：
[0086] S101：通过现场观测获得地形数据，建立数值标高模型。
[0087] S102：采用有限差分法求解泥石流的连续方程和运动方程，利用地球信息系统(GIS)中的数字高程模型(DEM)，自动生成地形变量，利用数字高程模型中的栅格网作为有限差分的网格；可以利用野外试验或实际工程验证模型和进行参数调整。
[0088] S103：使用基于DEM栅格网的有限差分公式可以实现数值解。
[0089] 下面结合附图对本发明的技术方案作进一步的描述。
[0090] 本发明实施例提供的基于GIS的大范围泥石流数值模拟及数值处理方法具体包括以下步骤：首先通过现场观测获得地形数据，建立数值标高模型。
[0091] 假定泥石流是均匀连续、不可压缩、非定常的牛顿流体，流动遵循下列方程，其连续性方程为：
[0092]
[0093] 运动方程为：
[0094]
[0095]
[0096]
[0097] 其中p为水压力，g为
重力加速度，u，v和w分别为x，y，z方向的平均流速，ρd是泥石流的平均密度。ρd＝ρsνs+ρwνw，其中ρs和ρw分别是固体颗粒和水的密度。
[0098] 根据浅水流运动原理，u,v>>w，方程(4)可以简化为：
[0099] p＝ρdg(η-ηb) (5)
[0100] 其中η是流体上表面的标高，ηb是流体底面的标高。
[0101] 流体底面的运动边界条件可以写为：
[0102]
[0103] 流体上表面的边界条件可以写为：
[0104]
[0105] 这样可以得到泥石流的连续方程为：
[0106]
[0107] 其中h为泥石流流深， 和 分别为x和y方向的平均流速。
[0108] 对方程(2)进行深度上的积分，可以得到：
[0109]
[0110]
[0111] 因此最终形式的深度积分运动方程为：
[0112]
[0113]
[0114] 其中α为速度修正系数，主要提供了有关垂直速度分布与均匀性偏差的信息，H为自由面的高度，ν＝μ/ρ为运动粘度系数，μ是粘度，θx和θy分别是地形沿x和y方向的倾斜角度， 为
动摩擦系数，t是流动时间。可以利用野外试验或实际工程来验证模型和进行参数调整。
[0115] 本发明采用有限差分法求解泥石流的连续方程和运动方程，利用地球信息系统(GIS)中的数字高程模型(DEM)，可以自动生成地形变量，利用数字高程模型中的栅格网作为有限差分的网格。在基于栅格的DEM分析中，每个像元都有八个可能的流动方向(左、右、上、下以及四个对角线)，如图2所示。单元格的流动方向可表示为：左＝0°，上＝90°，右＝180°，下＝270°，对角线为45°，135°，225°，315°。在单元格中泥石流沿一个方向流动，流动方向为落差最大的相邻单元格，如图3所示。
[0116] 使用基于DEM栅格网的有限差分公式可以实现数值解，如图4所示。控制方程的有限差分形式为：
[0117]
[0118] 运动方程(x方向)的有限差分形式为：
[0119]
[0120] 宾汉启动条件：
[0121] 宾汉流体在没有达到一定切应力(τ0)之前不会由任何剪切速率(没有流动，因此也没有速度)。描述宾汉流体需要两个参数：屈服应力τ0和粘度μ：
[0122]
[0123] 流体当达到：
[0124] τ>τ0 (16)
[0125] 时开始流动。对于层流，根据粘性流体的牛顿法则，水平方向的切应力：
[0126]
[0127] 对上述方程在深度方向上积分，得到：
[0128]
[0129] 该流动是与水平成θ角的平面上稳定均匀的宾汉流动。屈服应力为：
[0130] τ0＝ρgsinθ (19)
[0131] 对上述方程在深度方向上积分，得到：
[0132]
[0133] 其中h＝h(x,y,t)＝η-ηb是流体的深度。
[0134] 宾汉流体从静止开始运动，启动条件为：
[0135]
[0136] 对上述启动条件进行有限差分，得到：
[0137]
[0138] 下面结合附图，对本发明的较佳实施例做进一步的描述。
[0139] 如图5、6、7所示。某泥石流位于日本熊本县阿苏大桥附近，其流域面积为1.5km2，长度达到800m，最大宽度280m，最大深度25m，总崩塌和侵蚀体积大约54万立方米，最终泥石流汇入河谷，并冲毁大桥。
[0140] 通过现场实地调查和现场取样，确定泥石流密度为2200kg/m3，粘度为0.11Pa·s，动摩擦系数为0.6，α确定为1。经过计算得到泥石流流经区域，并得到泥石流到达河谷的最大速度为12.24m/s。计算流经区域跟现场基本吻合；泥石流冲进河谷后为只受重力的抛射体，根据在现场河谷对岸发现的泥石流体，验证计算所得泥石流速度无误。
[0141] 以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何
修改、s等同替换和改进等，均应包含在本发明的保护范围之内。
