本发明涉及地质灾害动力过程数值仿真计算领域，尤其涉及一种大规模泥石流地质灾害快速数值仿真模拟方法及系统。
背景技术：
地质灾害是指在地球的发展演化过程中，由各种地质作用形成的灾害性地质事件。我国是世界上灾害最严重的国家之一，90年代以来，我国每年的灾害造成了成千上万人死亡、造成直接经济损失高达数百亿元至数千亿元，相当国民生产总值(gnp)的3％-6％。随着灾种的增多、灾率的增加、灾度的加重，其危害的严重性和环境问题一样，成为制约我国可持续发展的重要因素之一。而要全面提升我国地质灾害防灾减灾救灾能力，不仅需要知道灾害发生后怎么做，更需要在灾害发生前判识和预测灾害的运动方式、轨迹以及危害范围。地质灾害的动力学机理和运动特性及其复杂，基于动力过程的数值模拟和定量评估在国际上仍然缺乏很好的解决方案，是当前国际前沿性科学问题和亟待突破的重要领域。
现阶段对泥石流地质灾害动力过程数值模拟主要基于连续介质、基于离散和介于二者三种方法。各种数值模拟方法都有其独特的功能和解决诸多问题的能力，并在综合应用中取得了大量卓有成效的研究成果。鉴于泥石流等地质灾害都是在重力驱动下的地表物质运动，它们都具有一个共同特征，即其运动方向的特征尺度远大于其垂直方向上的特征高度。但如何做到对大规模泥石流地质灾害(尤其是灾害链)的动力过程和运动特性的快速仿真模拟，是国内乃至国际上待解决的问题。其主要难点体现在以下几个方面：1.区域大，网格数量多，计算规模较大；2.计算速率慢，计算时间长。而一旦灾害发生后，对灾情的评估刻不容缓，如何快速做到对大规模泥石流灾害的评估甚至对大型灾害发生前快速判识和预测灾害的运动方式、轨迹以及危害范围，为相关防灾减灾部门提供及时的数据支持，从而减少甚至避免灾害给人类及人类工程活动造成的损失仍然是目前亟需解决的问题。
目前，尚无基于区域重划和结构化网格建立对大规模泥石流地质灾害动力过程的快速数值仿真模拟方法及系统的报道。
技术实现要素：
为了克服上述不足，本发明的目的是通过对大规模网格进行区域划分，对物源比例不同的区域采用不同等级的粗细化网格，重点关注区域细化，其它区域网格粗化，提供一种基于区域重划和结构化网格的大规模泥石流地质灾害快速数值仿真模拟方法及系统，克服现有技术中对大规模泥石流地质灾害数值模拟时区域大、网格数量多、计算速率慢的实际问题。
为了实现上述目的，本发明采用如下技术方案：
一种大规模泥石流地质灾害快速数值仿真模拟方法，包括如下步骤：
步骤1：获取泥石流地质灾害地形数据、物源区范围、高程以及物理力学参数；
步骤2：划分大规模泥石流地质灾害计算区域；
步骤3：对粗细化网格条件进行判断，并对计算区域网格粗化、网格细化和网格结构化；
步骤4：利用数值模拟软件对地表灾害动力过程进行数值模拟计算；
步骤5：输出n个计算区域的计算结果，并将n个计算区域的计算结果合并成一个计算区域输出。
优选的，所述步骤1中获取泥石流地质灾害地形数据的方法具体为：利用无人机采集地形数据，并用三维建模软件生成dem数据。
更优选的，所述利用无人机采集地形数据时，根据泥石流地质灾害范围采取定点飞行、s形飞行或回形飞行航线飞行，所述三维建模软件利用smart3d、pix4d或photoscan生成dem数据。
优选的，所述步骤1中获取泥石流地质灾害物源区范围及高程获取的方法具体为：
根据测量的泥石流的流速、泥石流断面平均泥深和泥石流水力坡度，采用泥石流断面处沟床纵坡，由式1泥石流流速公式计算得到泥石流沟床粗糙率系数，并根据式2计算得泥石流物源区范围及高程：
vc＝mc×hc2/3×ic1/2
式1
其中，vc表示泥石流的流速，mc表示泥石流沟床糙率系数；hc表示泥石流断面平均泥深(m)；ic表示泥石流水力坡度(‰)；
h＝k5/3*(q2)0.6
k＝mc*ic1/2
q2＝q/l
式2
其中，h表示物源区泥石流高度；k表示泥石流泥沙修正系数，通过mc、ic进行求解，q2表示泥石流沟道横断面在x方向的速度，通过流量q、l进行求解；q表示通过泥石流沟道横断面的流量；l表示通过泥石流沟道横断面的长度。
优选的，所述步骤2中划分大规模泥石流地质灾害计算区域的方法具体为：利用简单的相除和求余，将得到的大规模泥石流地质灾害dem数据划分为n*n个区域，然后将之前大区域所涉及的计算参数赋值给n*n个区域，保证区域划分后的等效性，其中n根据不同的计算精度和计算范围进行取值。
优选的，所述的步骤3中对粗细化网格条件进行判断的方法具体包括：根据计算对象对不同计算区域进行网格粗细化，利用物源作为判断条件，即没有物源、不需要计算的区域网格粗化，有物源需要计算的区域网格细化，同时选择粗化、细化的等级，比如计算区域物源量超过50％的将网格不断细化两级，计算区域有物源但物源量小于50％的将网格细化一级；或利用计算时步作为判断条件，比如初始时步对网格进行粗化，随着计算进行，物源量增加，后续时步网格粗化；同时根据自己计算需要和计算速率等采用不同的粗细化条件。
优选的，所述步骤3中计算区域网格细化的方法具体包括：
(1)设置细化等级，每一级细化包括对网格节点数、网格大小以及网格虚拟节点数的细化，具体实现如式3：
(nx)'＝2*(nx)-1
(ny)'＝2*(ny)-1
(dx)'＝0.5*(dx)
(dy)'＝0.5*(dy)
(mbc)'＝2*(mbc)
式3
其中，nx、(nx)'表示划分前、后计算网格x方向节点数；ny、(ny)'表示划分前、后计算网(dy)'格y方向节点数；dx、(dx)'表示划分前、后计算网格x方向单元长度；dy表示划分前、后计算网格y方向单元长度；mbc、(mbc)'表示划分前、后计算网格虚拟节点数；
(2)对细化后单元节点赋值：分别包括与细化前单元共节点的点进行赋值，在x和y方向新细化单元节点在两个单元节点之间的节点进行赋值，新细化节点在矩形单元中间的节点进行赋值；
(3)对细划后的区域赋予边界条件，边界条件包括对称边界条件、开边界条件、墙边界条件。
优选的，所述步骤3中计算区域网格粗化的方法具体包括：
(1)设置粗化等级，每一级粗化包括网格节点数、网格大小以及网格虚拟节点数的粗化，每一级粗化分为当前节点数为奇数还是偶数，若当前节点数为偶数时，将最后一个节点作为虚拟节点；
当节点数为奇数时，每一级粗化如式4：
(nx)'＝((nx)+1)/2
(ny)'＝((ny)+1)/2
(dx)'＝2*(dx)
(dy)'＝2*(dy)
(mbc)'＝(mbc)/2
式4
其中，nx、(nx)'表示划分前、后计算网格x方向节点数；ny、(ny)'表示划分前、后计算网格y方向节点数；dx、(dx)'表示划分前、后计算网格x方向单元长度；dy、(dy)'表示划分前、后计算网格y方向单元长度；mbc、(mbc)'表示划分前、后计算网格虚拟节点数；
(2)对粗化后单元节点赋值：主要对新粗化单元与原单元共节点的点进行赋值，考虑到节点数可能为偶数，赋值过程需提前判断粗化后节点是否超过粗化前节点数，然后进行赋值；
(3)对粗化后的区域赋予边界条件，边界条件主要包括对称边界条件、开边界条件、墙边界条件。
优选的，所述步骤3中计算区域网格结构化的方法具体包括：
查询当前区域的虚拟节点单元位于哪个目标区域内；
查询当前区域的虚拟节点单元在目标区域的哪个单元内；
利用式5-式7形函数插值计算当前区域的虚拟节点单元的具体参数：
n1＝β1(x-1)(y-1)
n2＝β2(x+1)(y-1)
n3＝β3(x+1)(y+1)
n4＝β4(x-1)(y+1)
式6
其中，u表示当前区域的虚拟节点单元的相应参数，如高度h，x方向速度的冲量hu，y方向速度的冲量hv等。
优选的，所述的步骤4具体包括：所述数值模拟软件采用massflow软件，根据步骤1、步骤2以及步骤3的结果对地表灾害动力过程进行数值模拟计算。
优选的，所述的步骤5具体包括：输出n个计算区域的计算结果，结果以dem或tecplot格式输出，根据不同格式利用arcgis或tecplot将n个计算区域的计算结果合并成一个区域进行输出。
本发明还提供一种大规模泥石流地质灾害快速数值仿真模拟系统，包括存储器及处理器，其中：
所述存储器用于存储程序指令；
所述处理器用于运行所述程序指令，以执行以下步骤：
获取地质灾害地形数据、物源区范围及高程h和物理力学参数；
划分大规模泥石流地质灾害计算区域；
计算区域网格粗化、网格细化和网格结构化；
利用模拟软件massflow对地表灾害动力过程进行数值模拟计算；
输出n个计算区域的计算结果，并将n个计算区域的计算结果合并成一个计算区域输出。
本发明与现有技术相比，具有如下优点：
(1)大规模网格计算，提高计算速率:本发明通过对大规模网格进行区域划分，对物源比例不同的区域采用不同等级的粗细化网格，重点关注区域细化，其它区域网格粗化，由于许多地质灾害其运动路径往往只占其坡面地形一小部分，上述策略既能避免对大规模区域中的一些非物源区域的计算，将大规模区域缩小到重点关注的小区域，保证计算精度，又能节约计算量，提高计算速率。
(2)动态调整计算区域网格的粗细化：本发明不仅以根据自己计算需要和计算速率等采用不同的粗细化条件，更关键的是，泥石流地质灾害是一个动态演化过程，随着泥石流地质灾害的动态演化，网格根据区域物源量等变化进行动态粗细化，从而网格随着计算的变化而变化，从而更加精确的计算需要重点计算的区域，比如在一个超大流域内，诸如桥墩、拦砂坝和排导槽等细部结构，需要精细网格才能刻画局部特征，因此结构化网格和随计算过程网格重划分优势明显。
(3)结构化网格和网格节点交互：本发明为实现区域的边界计算和边界节点的拟合，需将网格节点结构化，传统的结构化网格生成技术通常利用间接法，通过其他网格进行生成，不急效率低，而且存在一定误差，本发明通过将当前区域的虚拟节点单元用邻近区域的计算节点代替，借助有限元中的形函数理论，通过“查询—交换”实现节点的交换，即查询到虚拟节点所在的目标区域、目标单元后，利用形函数插值来求解当前区域的虚拟节点单元的具体参数。该方法相对于传统方法，计算效率更高，同时尽量减少插值后的误差。
附图说明
为了更清楚地说明本发明实施方式的技术方案，下面将对实施方式中所需要使用的附图作简单地介绍，应当理解，以下附图仅示出了本发明的某些实施例，因此不应被看作是对范围的限定，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他相关的附图。
图1是本发明实施例提供的泥湾沟地形数据处理后的结果图；
图2是本发明实施例中对区域划分和对划分后区域计算区域网格粗细化和网格结构化的结果图；
图3是本发明实施例中通过数值模拟软件计算后的结果图。
具体实施方式
为使本发明实施方式的目的、技术方案和优点更加清楚，下面将结合本发明实施方式中的附图，对本发明实施方式中的技术方案进行清楚、完整地描述，显然，所描述的实施方式是本发明一部分实施方式，而不是全部的实施方式。基于本发明中的实施方式，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施方式，都属于本发明保护的范围。
因此，以下对在附图中提供的本发明的实施方式的详细描述并非旨在限制要求保护的本发明的范围，而是仅仅表示本发明的选定实施方式。基于本发明中的实施方式，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施方式，都属于本发明保护的范围。
实施例
以泥湾沟泥石流为例，一种大规模泥石流地质灾害快速数值仿真模拟方法，包括以下几个步骤：
1、(1)通过对泥湾沟泥石流区域调查，利用无人机采用s形飞行获取影像数据，借助smart3d生成dem泥湾沟影像和dsm影像。如图1所示，其中(a)为泥湾沟地形数字高程(dem)影像，(b)为泥湾沟数字地表模型(dsm)影像。
(2)选取其中一段副沟进行定量模拟计算，同时利用泥石流流速和流量公式获得该泥石流物源区范围及高程h，具体公式如式8：
q2＝(-0.00247*(t-450)2+500.4)/6
h＝0.4885*q20.6
式8
其中，q2代表在x方向速度的冲量，t代表当前时步，h代表泥石流物源区高程。
(3)通过相应的物理力学实验测得材料密度为1825g/cm3，内聚力为10000，基底摩擦系数为0.35。
2、泥湾沟流域面积10.3km2，规模较大，可以将该区域划分成36个小区域，只计算泥石流运动路径部分，其他部分不予计算。
3、对计算区域网格粗细化和结构化网格，对计算区域物源量超过50％的将网格不断细化两级，计算区域有物源但物源量小于50％的将网格细化一级，没有物源且不需要计算的区域粗化两级。同时动态粗细化网格，即保证网格根据区域物源量流动等变化进行动态粗细化，随计算的变化而变化。如图2所示，其中(a)、(b)、(c)、(d)分别输出了时步为1s、3s、5s、7s后的结果。
4、将上述结果采用massflow等数值模拟软件对地表灾害动力过程进行数值模拟计算。
其中，massflow(即地表过程动力学数值模拟软件)是一款地表动力过程数值模拟软件，是中科院研究团队开发完成，并由成都山地环安防灾减灾技术有限公司销售的一款商业化软件。
每隔2s一个计算时步输出一次，将输出输出的36个计算区域的计算结果，利用tecplot合并成一个计算区域输出。通过数值模拟软件计算后的部分计算结果如图3所示。
由于许多地质灾害其运动路径往往只占其坡面地形一小部分，上述策略既能保证精度又能节约计算量，提高计算速率。在一个超大流域内，诸如桥墩、拦砂坝和排导槽等细部结构，需要精细网格才能刻画局部特征，因此结构化网格和随计算过程网格重划分优势明显。
上述实施例，通过对计算区域进行划分，对物源比例不同的区域采用不同等级的粗细化网格，重点关注区域细化，其它区域网格粗化。整个计算时长为30.2s，而利用传统的数值模拟(不用区域划分和网格结构化)计算总时长为40s左右，见明显缩短了计算时长，提高了计算速率。同时网格根据区域物源量流动变化进行动态粗细化调整，重点关注物源量较多的区域，提高了计算精度，迅速模拟大规模地质灾害动力演化过程，从而为防灾减灾部门提供更快更精确的数据支持。
以上所述仅为本发明的优选实施方式而已，并不用于限制本发明，对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。