青藏高原泥石流的形成及运动过程研究

2019-03-02王彬张全贵

广西水利水电 2019年1期

王彬，张全贵

（青海省有色地质矿产勘查局八队，西宁810001）

0 引言

泥石流主要是由砂、土、石等固态物质与水以及气体组成的混合物，在重力的驱动下沿着山坡或者沟谷进行运动。泥石流具有来势凶猛、爆发突然且破坏力极强的特点。泥石流属于地质灾害，其分布范围很广并且运动时间很长[1]。每年都会有泥石流发生，并且程度规模大不相同，给人们的生活带来了重大的影响。根据调查发现，青藏高原是我国泥石流高发的地区，青藏高原的构造运动与气候变化是引起泥石流高发的根本原因。青藏高原泥石流的发生，造成了巨大的经济损失[2]。因此，本文主要对青藏高原泥石流的形成过程进行研究。

泥石流的破坏性较大，受到了社会以及学术界的重视，关于对泥石流的研究也越来越多[3]。传统的泥石流模拟模型具有与实际情况吻合度较差的缺陷。因此根据收集到的资料与现场调查为数据，本文构建泥石流形成及运动模拟模型，对青藏高原泥石流的运动过程进行模拟与分析，为泥石流的防御与治理提供数据支撑。

1 泥石流形成运动模拟模型的构建

泥石流发生比较突然且破坏力极强，并且发生的地点往往较为偏远，对于同一地点来说，其重复发生的间隔时间会很长，因此，泥石流发生的全过程很难被观察到，这对泥石流的防御与治理及其不利。为了解决上述难题，构建泥石流形成、运动模拟模型[4]。首先建立泥石流形成、运动控制方程，然后采用数值分析法构建泥石流的模拟模型，通过计算机来模拟计算泥石流的形成、运动过程，将模拟结果呈现出来，为泥石流的防御与治理提供参考。

1.1 泥石流形成运动控制方程的建立

泥石流形成、运动模拟模型的前提是建设泥石流在运动的过程中质量不变，也就是说其遵守质量守恒原则，同时将泥石流中的流体与固体之间的相对运动忽略，将其看作单相流体[5]。将泥石流整个过程放在二维模型中，因此需要考虑x方向和y方向的质量守恒。另外，在泥石流的形成运动过程中，会有堆积和侵蚀的现象发生，并且堆积、侵蚀以及降雨都会对泥石流产生较大的影响，不容忽视。

根据上述影响因素，设定泥石流连续方程为：

其中：ξ=，ρ为石流密度，单位为kg/m3；T为泥石流的厚度，单位为m；u,v分别为x方向与y方向的泥石流速度，单位为m/s；t为时间；ξ为修正系数；E为侵蚀率；Tb为河床的高度。

若将侵蚀与降雨的影响忽略，则公式为：

而在实际情况中，降雨量对泥石流的影响非常大，因此，本文在建立泥石流形成、运动控制方程时，增加了降雨量与流失量两个项，则上述公式表示为：

其中：P为降雨的增加量；Q为降雨的流失量，单位为m/s。

假设泥石流在形成、运动过程中没有间隙，也就是泥石流的运动是连续的，因此，泥石流形成运动过程中的动量守恒方程为：

其中：β为动量分配系数；gx,gy,gz分别为在x,y,z轴上的重力分量；K为泥石流土的压力系数；ε为泥石流的摩擦阻力。经过动量守恒方程的建立，可以减少模型计算的时间，使模型计算简化，加快模型计算的效率。

采用Voellmy模型来对泥石流运动过程中遇到的抗剪应力进行计算，抗剪应力主要是摩擦力与阻力之和，其表达式为：

其中：σ为正应力；μ为摩擦系数；ε为泥石流的剪应力。

根据青藏高原泥石流的运动特征与运动模式，为了尽可能使计算变得简洁高效，本文建立的泥石流形成运动过程控制方程将不考虑泥石流密度的变化、侵蚀的情况与动量分配系数，也就是将β设定为1，将控制方程进行简化为

1.2 泥石流厚度与流速的计算

本文主要采用有限差分数值方法对控制方程进行计算，将上述公式转化为矢量格式为：

其中：X=

采用算子分裂方法将上述方程进行分解，将二维方程降为一维方程进行计算，其公式为：

该方程具有稳定高效的特点，保证了计算的稳健性与有效性。

2 青藏高原泥石流形成运动过程模拟的实现

首先，对模型参数进行选取，主要包含两部分内容，一是泥石流运动的初始条件，二是泥石流的运动条件。

初始条件主要包括泥石流的位置、厚度、模拟的计算范围以及泥石流的初始位置与初始量。

通过DEM数字高程模型来对青藏高原的地形条件进行测量与输入[6]。根据所得的地形情况，将易发生泥石流的位置进行切割，设定为研究区域，只针对这部分进行计算，其他部分不参与计算，这样可以提升计算的效率。

对曾经发生的泥石流进行研究，发现青藏高原的物源条件主要来源于青藏高原坡面上风化的基岩，岩石经过风化会产生大量的泥土与石块，还会带有枯木等杂物，泥石流发生后，这些物源会产生堆积，若是再次发生泥石流，物源就会发生再次转移[7]。因此，根据上述描述对初始条件进行设置，将平均物源厚度定义为0.01 m，在主泥石流运动方向内的物源厚度定义为0.1 m，在沟口处的物源厚度定义为1 m。

根据相关调查，泥石流的发生频率和规模与降雨频率和强度成正比例关系，因此，在进行泥石流形成运动模拟中对降雨频率进行相应的设置，本文主要采用设计暴雨来对降雨强度进行表示。在研究区域内，设置两处暴雨点，分别为A和B，A的权重为0.8，B的权重为0.2，根据模比系数计算相应的暴雨量[8]。

由于泥石流属突发性灾难，无法测量泥石流的容重。因此，本文根据研究区的地形条件以及泥石流易发程度量化评分标准，对模拟泥石流容重进行计算，得到泥石流的容重为1.9 t/m3。

本文采用的是Voellmy模型，该模型中包括两个参数σ和μ，σ为湍流系数，μ为摩擦系数，通过模拟还原青藏高原泥石流的运动过程，将计算结果与真实情况进行对比分析，获得泥石流的运动参数。

根据记载，青藏高原泥石流持续时间约为5 min，厚度为2～5 m，速度为6～8 m/s，随着降雨量的增加，泥石流持续运动[9]。

通过实际情况设置相应的运动参数进行模拟，将μ范围设置为0.1～0.5，σ范围设置为10～90，其他参数一致。经过多次实验，将得到的结果与实际情况进行比较，发现当μ为0.3，σ为10时，与实际情况最吻合。因此，将μ设定为0.3，σ设定为10。

然后，对青藏高原泥石流进行模拟。根据研究发现，青藏高原泥石流发生的周期大约为20 a。因此，本文对20年一遇的泥石流形成运动过程进行模拟。

该次模拟设定时间为5 min，分别对0.5 min、1 min、3 min与5 min的泥石流厚度分布、流速进行模拟计算。

2.1 泥石流厚度分布情况

当泥石流运动时间为0.5 min时，泥石流最大厚度为0.4 m，主要分布在主沟内，其他支沟的平均厚度为0.2 m，整体来看研究区内的泥石流厚度比较小，由于泥石流形成运动的时间较少，导致泥石流的运动有间隙；当泥石流运动时间为1 min时，泥石流最大厚度为0.54 m，主要分布在主沟内，其他支沟的平均厚度不足0.3 m，整体来看泥石流的运动呈现不连续的状态；当泥石流运动时间为3 min时，泥石流最大厚度为1.6 m，主要分布在上游集水点，其他支沟的平均厚度不足0.7 m，并且处于支沟的下游部分；当泥石流运动时间为5 min时，泥石流最大厚度为2.4 m，主要分布在主沟下游，该部分属于缓坡，容易产生堆积的现象。其他支沟的厚度不足0.6 m，大部分支沟没有泥石流的堆积现象。

2.2 泥石流流速分布情况

整体来看，泥石流的流速是随着时间的增加逐渐减小的。当泥石流运动时间为0.5 min时，流速为24 m/s，主要分布于支沟两侧的坡体上，由于坡度大，高度差较大，产生的重力也会随之增加，泥石流的流速会相对较快，最大值可以达到20～24 m/s，到达支沟下游时，流速会减小到15 m/s，到达主沟后流速会再次减少为10 m/s，在主沟的中心部分流速为5 m/s；当泥石流运动时间为1 min时，泥石流会向主沟聚集，其流速约为5～10 m/s；当泥石流运动时间为3 min时，泥石流几乎全部位于主沟，其流速约为15 m/s；当泥石流运动时间为5 min时，泥石流流速的最大值约为12 m/s[10]。

3 实验结果与分析

为了保证本文构建的泥石流形成、运动模拟模型的有效性，设计实验对其进行验证。在实验过程中，将泥石流作为实验对象，主要是对其形成、运动过程进行模拟与分析。为了保证实验过程与结果的准确性，使用构建的泥石流形成、运动模拟模型与传统模拟模型进行比较，观察实验对比结果。在实验过程中，将传统模拟模型称为对照组，构建的泥石流形成、运动模拟模型称为实验组。

3.1 数据准备

为了尽可能地保障实验结果的准确性，对实验过程中的参数进行相应的设置，本文采用不同的模拟模型对泥石流形成、运动过程进行模拟与分析，由于采用的模型不同，因此，在实验过程中必须保证外部环境参数的一致。本文实验参数设置结果如表1所示。

表1 实验参数设置结果

3.2 实验对比结果

在实验过程中，由于采用的模拟模型不同，因此，本文利用ASO100软件对实验数据进行记录与分析。主要通过模拟结果与实际情况的吻合度来验证模拟模型的有效性。实验对比结果如图1、图2所示。

图1 泥石流最大厚度实验对比图

图2 泥石流最大流速实验对比图

由图1可知，实验组和对照组与实际情况的泥石流最大厚度的变化趋势一致的，但是，实验组与实际情况的吻合度更高。由图2可以看出，实验组与实际情况的泥石流最大流速的变化趋势是一致的，而对照组的变化趋势在运动时间3～4 min之内与实际情况是相反的，实际情况的变化趋势是下降，而对照组的变化趋势是上升，明显与实际情况不吻合。因此，实验组与实际情况吻合度更高。

综上所述，实验组与实际情况的吻合度远远高于对照组，说明本文构建的泥石流形成、运动模拟模型具备较高的有效性。