基于GPU的岩石碎屑流与拦砂坝交互场景的三维建模与可视化

叶健，陶和平，陈锦雄，陈晓清

(中国科学院 成都山地灾害与环境研究所，四川 成都，610041)

我国幅员辽阔，地貌差异显著，地质环境复杂，灾害性气象及气候条件出现频繁，再加上不断扩大的人类活动因素，使得山地灾害频繁发生，严重威胁山区建设和人民生命财产安全。其中，高速远程滑坡-碎屑流是一种特殊的地质灾害，它是指高速远程滑坡或崩塌在运动过程中转化而成的碎屑流体[1]。目前，很多矿山、采石场等将废弃物放于附近的山坡上，当遇到地震或附近的爆破作业等情况时，或者当滑坡体与山体碰撞粉碎后，散体堆积物就会沿坡滑下[2-4]，形成碎屑流。形成后的碎屑流以高速度和大的位移引发灾难性事故，往往会造成严重的生命财产损失[1]。2008年5月12日汶川地震后，文家沟顶部距离汶川地震发震断层仅3.6 km处，山体在长时间强烈地震震动作用下失稳而形成滑坡。滑坡在高速下滑过程中转化为碎屑流[5]。2009年5月发生在我国重庆武隆的鸡尾山滑坡-碎屑流，运行了1.5 km，等值摩擦因数仅为0.27，高速运行的碎屑流掩埋了下游的铁矿，造成 10人死亡，64人失踪[6]。由于碎屑流灾害近年来的频繁发生，相关部门已着手研究碎屑流防灾减灾对策。在碎屑流防灾减灾的措施中，工程治理是防治岩石碎屑流灾害的主要措施之一，应用极其广泛，对消除和减轻地质灾害对人民生命财产的危害发挥了重要的作用。拦砂坝是工程防治中较常采用的手段，利用拦砂坝可抑制碎屑流的发生，消减碎屑流流量，延缓碎屑流到达时间，减少碎屑流土砂总流出量，阻滞碎屑流颗粒输移、减势、稳定沟床及减轻两岸侵蚀[7]。在进行碎屑流工程防治的过程中，研究碎屑流运动规律，逐步建立和完善碎屑流模型，并根据碎屑流灾害防治工程设计规范建立拦沙坝模型，使其能够以三维可视化的方式再现碎屑流运动和堆积的过程，将有助于人们加深对碎屑流灾害的认识，并有助于进行灾害的防灾减灾措施。尽管目前在地质学领域里，已经有不少学者利用离散元法对碎屑流进行数值模拟，但把离散元法运用到三维仿真领域并进行碎屑流与拦砂坝、排导槽交互的可视化工作却较少。因此，本研究旨在将三维可视化理论与技术引入到碎屑流现象及其防灾减灾的三维建模与可视化中，建立碎屑流模型，不仅可以利用传统的数值模拟技术，同时可以借助动态、三维交互的三维仿真技术来模拟碎屑流现象及其碎屑流的防灾减灾过程。在利用计算机模拟离散物体现象的方法中，基于散体理论建立的离散元(discrete element method, DEM)法具有不依赖网格和描述破坏问题便利的特点，能够更好地反映一些自然现象发生过程的本质，它是散体力学分析的一种有效的工具[8-9]，该算法非常适合于进行碎屑流现象的研究。而另一方面，随着图形技术的飞速发展,近年来传统的图形硬件已经演化为可编程图形处理单元(GPU)，其强大的计算性能、灵活的可编程性以及面向个人消费市场的低廉价格，吸引了越来越多的研究者将 GPU 应用于解决图形渲染之外的通用计算(GPGPU)[10-11]，并利用GPU高效的并行计算来进行复杂的物理计算。本研究基于散体理论的离散元法建立的碎屑流模型，呈现出高度的数据并行性，因此，可以充分利用GPU的高度并行性来进行碎屑流的建模与三维可视化的实时渲染。可视化场景中拦砂坝和排导槽的建模及两者与碎屑流之间边界问题的正确处理在本研究中较重要，直接影响到岩石碎屑流与拦砂坝和排导槽交互场景的可视化表达。为此，本文作者首先分析碎屑流在场景中运动的整个过程；然后，构建拦砂坝和排导槽模型，并确立它们的边界；最后，根据碎屑流颗粒与拦砂坝和排导槽边界之间的受力分析，建立有效的碎屑流颗粒与拦砂坝和排导槽边界问题的处理方法，并充分利用GPU的高度并行性和可编程性，实现岩石碎屑流与拦砂坝和排导槽交互场景的三维建模与可视化。

1 基于离散元法的岩石碎屑流建模

1.1 岩石碎屑流运动特点

岩石碎屑流常常是由高速滑坡的发生而引起的，但岩石碎屑流与滑坡不同，不是沿滑动面的整体，而是呈颗粒介质的运动，液体对运动不起重要作用[12]。从启动到静止，高速远程滑坡或崩滑在运动过程中逐渐转化为岩石碎屑流，最终以高速远程滑坡-碎屑流的复合形式出现[13]。

岩石碎屑流发生的整个过程通常可归纳为启程阶段、加速阶段和堆积阶段。

第一阶段(启程阶段)：当遇到地震或附近的爆破作业等情况时，或者当滑坡体与山体碰撞粉碎后，散体堆积物就会沿坡滑下。

第二阶段(加速阶段)：高速岩石崩滑体在经历了启程阶段后，与运动路径两侧或前方突出的阻挡山体发生强烈的冲击碰撞崩解形成岩石碎屑流，并以这种形态进入加速阶段，此时，岩石碎屑流势能不断释放，一部分转化为岩石碎屑流动能，一部分转化为摩擦热能。

第三阶段(堆积阶段)：岩石碎屑流的运动随着山体的坡角减小，阻力变大，这时随着动能慢慢变小，岩石碎屑流运动速度也越来越小，并开始堆积，直至所有碎屑物质堆积完毕为止。

岩石碎屑流在运动过程中常常呈现“流态化”的特征，这是它能高速运动的主要原因之一[14-18]。在岩石碎屑流的启动、加速和堆积的整个过程中，岩石碎屑流不会产生大体积的“飞行”运动，且岩石碎屑流空隙度较大，具有离散性、流动性和高运动性的特点。

1.2 基于离散元法的岩石碎屑流模型

离散元法最早是由Cundall于1971年提出的一种不连续的数值方法。其基本思想是把不连续体分离为刚性元素的集合，使各个刚性元素满足运动方程，用时步迭代的方法求解各刚性元素的运动方程，继而求得不连续体的整体运动形态[19-20]。离散元法自问世以来，在岩土工程和颗粒散体工程这两大传统的应用领域中发挥了其他数值算法不可替代的作用，特别是在计算和模拟材料具体的破坏形式和破坏的整个过程方面，离散元法显示出明显的优势[21-22]。离散元法能模拟岩石颗粒之间存在的滑动、平移和转动等复杂过程，具有宏观上的不连续性，可以较真实地模拟碎屑流颗粒的运动，也可以模拟高速大位移运动。由于高速远程滑坡后所形成的岩石碎屑流具有离散度大、流动显著且高速大位移运动等特点，因此，适合于利用离散元法进行模拟。

离散元法是基于牛顿第二运动定理和颗粒在接触点上的力与位移之间的本构关系所建立的数值模拟方法[23-24]。在颗粒运动的过程中，首先由外边界条件决定颗粒所受外力，然后根据颗粒之间的位置关系来确定颗粒间的相互作用，最后求出颗粒所受的合力，根据这个合力求出颗粒的加速度和速度，进而求出各个颗粒的位移。

当2个颗粒之间发生碰撞，颗粒在前进方向受到阻力，此时颗粒之间会产生接触力Fc，并且颗粒在接触点处发生变形，当颗粒达到最大位移变形时，颗粒处于静止状态，然后沿运动的方向反弹。若此时颗粒间的碰撞为偏心碰撞，则相撞点处的接触力可以分解为垂直于接触面的法向接触力Fcn和平行于接触面的切向接触力Fcs。颗粒之间的接触模型是离散元法的核心，本模型在颗粒之间施加了变形元件弹簧和阻尼器，采用的是弹簧-阻尼器模型[25]，模型中颗粒j和颗粒i之间的法向接触力Fcn可以表示为：

其中：fi,s为颗粒i和颗粒j碰撞的排斥力，颗粒在碰撞过程中，颗粒阻尼力fi,d使用阻尼器模拟；κ，η，d，ri,j和vi,j分别为弹簧弹性系数、阻尼系数、颗粒直径，相对位置和颗粒j对颗粒i的相对速度。

碰撞后，颗粒的动能会产生一定的损失，损失的大小与颗粒的弹性阻尼系数及颗粒间的相对速度有关。而切向接触力Fcs与相对切向速度vij,t可描述为：

相对切向速度vij,t的计算公式为：

作用于刚体的力Fc和力矩Tc为施加于所用刚体颗粒之上的力和力矩的和：

其中：ri为当前颗粒i质心的相对位置。

1.3 离散元法在GPU通用计算平台上的实现

为了更好地模拟岩石碎屑流的细节，碎屑流颗粒的数目往往达到数十万个，若此时将碎屑流颗粒间相互作用的计算完全放在CPU上，碎屑流运动的模拟几乎不能满足实时性和交互性的要求，而利用离散元法模拟碎屑流可以实现碎屑流颗粒之间的交互计算，非常适合在GPU上实现。在GPU通用计算平台上，由于将每个颗粒的计算都分配到每个GPU线程上执行，这样就可以充分利用GPU的并行性，可极大地加快算法的计算速度，从而满足数十万个碎屑流颗粒相互作用的模拟场景实时性和交互性的要求。

在本研究所构建的碎屑流模型中，碎屑流颗粒间相互作用力的计算需要建立高效的颗粒查找和排序机制以及颗粒间碰撞检测机制，因此，首先需要查询邻近颗粒并计算它们之间的碰撞。为了避免对邻近颗粒进行n2次搜索，采用邻近颗粒搜索算法，该算法与文献[26]中的相似，可以概括为：

(1) 划分仿真区域为统一网格；

(2) 利用每一颗粒的空间位置去查找它所属的单元；

(3) 把颗粒单元位置作为哈希函数的输入；

(4) 根据空间哈希排序颗粒；

(5) 根据它们的哈希值在线性缓冲区中重新排序颗粒。

因此，把三维空间划分为统一均匀的网格单元，每一颗粒被赋予一个基于中心点的网格单元，颗粒的索引被存储在这个网格单元中，并以该网格单元索引作为该颗粒的哈希值。一旦建立了网格结构，就能使用它来加速颗粒和颗粒之间的交互计算。

当计算颗粒所在的网格单元时，可以循环遍历相邻的26个网格单元(3×3×3个单元)并检测是否与这26个网格单元中的颗粒之间发生碰撞：若发生碰撞，则通过计算便可得到该颗粒所在的网格单元，继而写入该网格单元索引到正确位置的网格数组中，而此时颗粒和网格的对应关系已经发生了变化。由于更新颗粒位置后的颗粒-网格对应关系是乱序的，必须通过排序算法以网格为序排列起来才能在计算中快速地取得邻近颗粒的信息。颗粒的排序是基于哈希值来排序颗粒的，采用的排序算法为并行基数排序[27]。利用并行基数算法可以得到颗粒和网格单元的正确对应关系，从而保证在计算中可以快速地取得邻近颗粒的信息。由于本研究将离散元法的全部模拟计算都分配到GPU中处理，因此，充分利用了GPU的高并行性和可编程性，可利用离散元法实现碎屑流的实时计算和模拟。

2 拦砂坝和排导槽的建模

工程治理是防治岩石碎屑流灾害的主要措施之一，应用极其广泛，在消除和减轻地质灾害对人民生命财产的危害中发挥了重要的作用。修建于沟道中的拦砂坝是工程治理的重要手段，通过在岩石碎屑流运动的通道中构建层层拦截的拦砂坝，以拦砂坝的规模效应形成了较强的拦沙、削峰、滞洪能力和上拦下保的作用。不仅如此，拦砂坝还具有固定河床，抬高河底高程，稳定边坡，遏制沟岸扩张、沟底下切和沟头前进，减轻沟道侵蚀，大大减少泥沙下泄量的作用[28]。而排导槽由于结构简单、就地取材、施工方便和效果显著等优点得到广泛应用[29]。

由于拦砂坝和排导槽在山地灾害工程治理中发挥了重要的作用，本研究将以岩石碎屑流作为研究对象，通过利用离散元法来建立岩石碎屑流模型，然后构建拦砂坝和排导槽两大岩石碎屑流防护措施模型，重点模拟岩石碎屑流形成以后,岩石碎屑流蓄满在第一级拦砂坝后被第二级拦砂坝拦截的防灾减灾过程。

为了构建碎屑流与拦沙坝、排导槽交互的虚拟场景，需要构建拦砂坝和排导槽的防护措施模型。要构建该防护措施模型，其关键在于如何确定和处理拦砂坝及其排导槽的边界。为了更好地再现虚拟场景中拦砂坝的防灾减灾过程，在虚拟场景中构建了二级岩石碎屑流拦砂坝，当岩石碎屑流发生时，拦砂坝能将岩石碎屑流拦截并堆积在第一级拦砂坝库内，直到该级拦砂坝蓄满。若大量岩石碎屑流从坝的上方溢出，则溢出的岩石碎屑流开始在第二级拦砂坝内蓄积，由于溢出的岩石碎屑流经过上一级拦砂坝的拦挡，其流量和规模将大幅度减少，流速也减小。经过二级拦砂坝拦蓄与消能后，其危害可以大大减小。

2.1 边界的确定

岩石碎屑流模型建立后，需要确定拦砂坝和排导槽的边界。其边界的俯视图如图1所示。

图1 拦沙坝和排导槽的边界平面图Fig.1 Boundary planform of debris dam and drainage canal

图2 拦沙坝与排导槽的剖面示意图Fig.2 The diagrammatic crosssection of debris dam and drainage canal

图2中排导槽的斜坡段L1的坡度α为30°，L2段为排导槽的水平段，岩石碎屑流自排导槽顶端在重力的作用下冲向排导槽底端，岩石碎屑流首先要触碰排导槽L1段的底端和两侧，然后被拦砂坝拦截并堆积在第一级拦砂坝A的库内，直到该级拦砂坝蓄满为止；然后，岩石碎屑流开始从拦砂坝A溢出并蓄积在拦砂坝B内，经过两级拦砂坝的拦截，岩石碎屑流基本被拦截在拦砂坝内，只有少量岩石碎屑流颗粒会落入排导槽L2段。

结合图1和图2可知：首先需要确定L1和L2段排导槽底端及其两侧的边界，然后，必须确定拦砂坝A和拦砂坝B的边界。根据以上分析，在场景中构建了拦砂坝和排导槽的模型，其效果图如图3所示。

图3 拦砂坝和排导槽的三维建模效果Fig.3 Effect of three-dimensional modeling of debris dam and drainage canal

2.2 边界的处理

在确定拦砂坝和排导槽的边界后，需要处理碎屑流颗粒与边界的碰撞。碎屑流颗粒与边界的碰撞过程如图4所示，颗粒首先从边界面的某一方向撞击边界面，在碰撞的过程中，当颗粒碰撞边界面的那一时刻，颗粒发生形变，在回弹的那一瞬间，由于能量损失，反弹后的速度减小。

假设颗粒与边界在碰撞前一瞬间的速度为v，这时，可以把颗粒沿着边界面的法向和切向分解为vn和vt。其颗粒与边界碰撞前和碰撞后速度分解示意图如图5所示，碰撞后的法向和切向速度将变为：

图4 岩石碎屑流颗粒与边界的碰撞过程示意图Fig.4 Sketch map of collision process of rock avalanche granule and boundary

图5 岩石碎屑流颗粒与边界碰撞前后Fig.5 Rock avalanche granule before and after colliding with boundary

其中：α为法向速度修改系数，α=1时，颗粒的碰撞无能量损失；α=0时，颗粒的法向分量速度为 0。本实验α取为0.6，碰撞后法向速度反向。β为表面摩擦因数，控制障碍物表面的滑移特性。当β=1 时，颗粒在切向碰撞时，无摩擦；当β=0 时，颗粒沿边界面的切向速度将变为0。本实验β取为0.75。为了进行颗粒与边界的碰撞检测及颗粒碰撞后速度的计算，需在程序初始化之前保存拦砂坝和排导槽边界距离场的信息，当颗粒与边界的距离d小于边界碰撞阈值T时，可以确定颗粒与边界发生碰撞，然后进行受力和速度的计算。由于本研究仅处理静止的边界，因此，当检测到碰撞且颗粒已进入边界的内部时，边界不会发生任何的移动或变化，仅仅需要处理颗粒的速度。

3 实验结果

为了验证拦砂坝拦截岩石碎屑流的效果，把实验分为有拦砂坝拦截和无拦砂坝拦截2组。

本实验采用Intel(R) Core (TM)2 Duo 3.00 GHz，2.0 GB内存，NVIDIA GeForce 8800GT 显卡的PC机。实验软件平台为 Windows XP和 VS2005，并使用NVIDIA CUDA作为GPU计算架构。

第1组实验如图6所示，在重力作用下，岩石碎屑流从排导槽的顶部滑下，岩石碎屑流的势能逐渐转换为动能，在没有阻挡物阻挡的情况下，岩石碎屑流的速度越来越快。从实验结果可以看出：岩石碎屑流的运动将对排导槽的水平段造成巨大的冲击。

图6 岩石碎屑流在没有拦砂坝拦截下的全过程Fig.6 Overall process of rock avalanche granule without debris dam

在第2组实验中，添加了拦砂坝，实验如图7所示。尽管岩石碎屑流在滑落的过程中速度越来越快，但由于受到第一级拦砂坝的消能和拦截，部分岩石碎屑流将蓄积在第一级拦砂坝内，而另一部分岩石碎屑流将继续冲击排导槽的水平段，在经过第二级拦砂坝的消能和拦截后，岩石碎屑流将基本蓄积在拦砂坝内，只有部分岩石碎屑流会飞溅到拦砂坝的水平段。从图7(b)可以看出：岩石碎屑流与拦沙坝的碰撞伴随着撞击和爬高，其可视化效果能够较真实地再现自然界中以下特点：当岩石碎屑流势能基本转化为动能后，岩石碎屑流运动速度快，惯性力大；当前进遇障碍阻挡时，可冲击爬高、翻越障碍等。并再现了岩石碎屑流堆积的过程。另外，为了模拟出较丰富的细节效果，当本模型初始化的粒子数达到122.2万个时，实时模拟速度仍然能达到10帧/s以上。

图7 岩石碎屑流被拦砂坝拦截的全过程Fig.6 Overall process of rock avalanche granule intercepted by debris flow

4 结论

(1) 利用离散元法构建了岩石碎屑流模型，并充分利用GPU的高度并行性特点，将岩石碎屑流模拟的大部分模拟放在GPU中计算，使得岩石碎屑流颗粒数为数百万时仍能满足实时的要求。

(2) 在绘制岩石碎屑流与拦砂坝和排导槽交互的场景中，通过有效地确定和处理岩石碎屑流与拦砂坝和排导槽的边界条件，模拟出岩石碎屑流蓄满在第一级拦砂坝后被第二级拦砂坝拦截的过程，同时也展示了岩石碎屑流与拦沙坝和排导槽交互过程中撞击和爬高的细节，并实现了岩石碎屑流堆积的过程。

(3) 为岩石碎屑流防灾工程的处理设计与施工提供了直观的可视化分析平台，此平台可结合岩石碎屑流灾害防治工程设计规范和岩石碎屑流运动参数的计算结果设计出的岩石碎屑流防灾减灾工程，再现岩石碎屑流防灾减灾的过程。

在岩石碎屑流的动态绘制中，尚未考虑岩石碎屑流颗粒级配的问题，所绘制的岩石碎屑流颗粒粒径是均匀的，这与自然现象中的岩石碎屑流有所不同。为了提高岩石碎屑流模拟的效果，应该考虑模拟不同颗粒半径的岩石碎屑流。另外，应充分发掘GPU更加广阔的通用计算性能，利用能够表达连续流体的物理模型绘制流体的动态场景，并将模型应用于工程设计和灾害评估。

[1]张明, 殷跃平, 吴树仁, 等.高速远程滑坡-碎屑流运动机理研究发展现状与展望[J].工程地质学报, 2010, 18(6)∶ 805-812.ZHANG Ming, YIN Yueping, WU Shuren, et al.Development status and prospects of studies on kinematics of long runout rock avalanches[J].Journal of Engineering Geology, 2010, 18(6)∶805-812.

[2]贺可强, 安振远.崩滑碎屑流的形成条件与形成类型[J].河北地质学院学报, 1996, 19(3/4)∶ 344-351.HE Keqiang, AN Zhenyuan.A preliminary analysis on the forming conditions and the Forming types of debris flow caused by landslides and falls[J].Journal of Hebei College of Geology,1996, 19(3/4)∶ 344-351.

[3]陈自生, 张晓刚.1994-04-30 四川省武隆鸡冠岭滑坡→崩塌→碎屑流→堵江灾害链[J].山地研究, 1994, 12(4)∶ 225-229.CHEN Zisheng, ZHANG Xiaogang.A hazard-chain of landslide→collapse→debris flow→river stoppage in Wulong county, Sichuan Province on April 30, 1994[J].Mountain Research, 1994, 12(4)∶ 225-229.

[4]朱平一, 王成华, 唐邦兴.西藏特大规模碎屑流堆积特征[J].山地学报, 2000, 18(5)∶ 453-456.ZHU Pingyi, WANG Chenghua, TANG Bangxing.The deposition characteristic of supper debris flow in Tibet[J].Journal of Mountain Science, 2000, 18(5)∶ 453-456.

[5]黄河清, 赵其华.汶川地震诱发文家沟巨型滑坡-碎屑流基本特征及成因机制初步分析[J].工程地质学报，2010, 18(2)∶168-177.HUANG Heqing, ZHAO Qihua.Basic characteristics and preliminary mechanlism analysis of large scale rockslide-sturzstrom at WenJiagou triggered by Wenchuan Earthquake[J].Journal of Engineering Geology, 2010, 18(2)∶168-177.

[6]殷跃平.斜倾厚层山体滑坡视向滑动机理研究∶ 以重庆武隆鸡尾山滑坡为例[J].岩石力学与工程学报, 2010, 29(2)∶217-226.YIN Yueping.Mechanism of apparent dip slide of inclined bedding rockslide∶ A case study of Jiweishan rockslide in Wulong, Chongqing[J].Chinese Journal of Rock Mechanics and Engineering, 2010, 29(2)∶ 217-226.

[7]熊道锟, 徐世民.泥石流拦挡坝之虞[J].中国地质灾害与防治学报, 2010, 21(4)∶ 136-138.XIONG Daokun, XU Shiming.The fear for debris dam[J].The Chinese Journal of Geological Hazard and Control, 2010, 21(4)∶136-138.

[8]成名, 刘维甫, 刘凯欣.高速冲击问题的离散元法数值模拟[J].计算力学学报, 2009, 26(4)∶ 591-594.CHENG Ming, LIU Weipu, LIU Kaixin.Numerical simulation for high-speed impact problems by discrete element method[J].Chinese Journal of Computational Mechanics, 2009, 26(4)∶591-594.

[9]杨洋, 唐寿高.颗粒流的离散元法模拟及其进展[J].中国粉体技术, 2006, 38(5)∶ 38-43.YANG Yang, TANG Gaoshou.Discrete element method simulation of granula flow and Its advances[J].China Powder Science and Technology, 2006, 38(5)∶ 38-43.

[10]Owens J D, Luebke D, Govindaraju N, et al.A survey of general-purpose computation on graphics hardware[J].Computer Graphics Forum, 2007, 26(1)∶ 802113.

[11]吴恩华, 柳有权.基于图形处理器(GPU)的通用计算[J].计算机辅助设计与图形学学报, 2004, 16(5)∶ 6012612.WU Enhua, LIU Youquan.General purpose computation on GPU[J].Journal of Computer-Aided Design & Computer Graphics, 2004, 16(5)∶ 6012612.

[12]鲁晓兵, 张旭辉, 崔鹏.碎屑流沿坡面运动的数值模拟[J].岩土力学, 2009, 30(增刊 2)∶ 524-527.LU Xiaobing, ZHANG Xuhui, CUI Peng.Numerical simulation of clastic grain flow along a slope[J].Rock and Soil Mechanics,2009, 30(Suppl 2)∶ 524-527.

[13]孙萍, 张永双, 殷跃平, 等.东河口滑坡-碎屑流高速远程运移机制探讨[J].工程地质学报, 2009, 17(6)∶ 737-744.SUN Ping, ZHANG Yongshuang, YIN Yueping, et al.Discussion on long run-out sliding mechanism of Donghekou landslide-debris flow[J].Journal of Engineering Geology, 2009,17(6)∶ 737-744.

[14]Hsu K J.Catastrophic debris streams(sturzstroms)generated by rockfalls[J].Geological Society of America Bulletin, 1975, 86∶129-140.

[15]Heim A, von Elm der Bergsturz.Deutsch[J].Geol Zeitschr, 1882,34∶ 74-115.

[16]Buss E, Heim A.von Elm der Bergsturz[J].Zurich, Wurster &Cie, 1881∶ 163.

[17]Bagnold R A.Deposition in the process of hydraulic transport[J].Sedimentolgy, 1968, 10∶ 45-56.

[18]Bagnold R A.The nature of saltation and of bed-load transport in water[C]//Proceedings, Royal Society of London, Series A , 1973,332∶ 473-504.

[19]CU NDALL P A.A computer model for simulating progressive large scale movement in block rock system[J].Symposium ISRM, 1971, 2∶ 129-136.

[20]刘凯欣, 高凌天.离散元法研究的评述[J].力学进展, 2003,33(4)∶ 483-490.LIU Kaixin, GAO Lingtian.A review on the discrete element method[J].Advances in Mechanics, 2003, 33(4)∶ 483-490.

[21]Liu K, Zheng W, Gao L, et al.A numerical analysis for stress wave propagation of anisotropic solids by discrete element method[C]//Chiba A, Tanimura s, Hokamoto K,eds.Proceedings of the 4th International Symposium on Impact Engineering.Kumamoto, Japan.UK∶ Elsevier Science Ltd, 2001∶ 589-594.

[22]Sawamoto Y, Tsubota H, Kasai Y, et al.Analytical studies on local damage to reinforced concrete structures under impact loading by discrete element method[J].Nuclear Engineering and Design, 1998, 179∶ 157-177.

[23]Zhu H P, Yu A B.A theoretical analysis of the force models in discrete element method[J].Powder Technol, 2006, 161∶122-129.

[24]Nezami E G, Hashash Y M.A fast contact detection algorithm for 3-D discrete element method[J].Compute Geotechnol, 2004, 31∶575-587.

[25]Mishra B K.A review of computer simulation of tumbling mills by the discrete element method[J].Part 1—Contact Mechanics.International Journal of Mineral Processing, 2003, 71(1)∶ 73-93.

[26]Hubert Nguyen.GPU Gems[M].3ed, Boston∶ Addison Wesley Professional, 2007∶ 531-549.

[27]Satish N, Harris M, Garland M.Designing efficient sorting algorithms for manycore gpus[R].Santa Clara∶ NVIDIA Technical Report NVR-2008-001, NVIDIA Corporation, 2008.

[28]李发斌, 涂建军, 胡凯衡.基于GIS的泥石流拦砂坝优化设计模型[J].中国科学 E 辑∶ 技术科学, 2003, 33(增刊)∶ 146-156.LI Fabing, TU Jianjun, HU Kaiheng.Design model of debris dam Based on GIS[J].Science in China Series E∶ Technological Sciences, 2003, 33(suppl)∶ 146-156.

[29]李德基.泥石流减灾理论与实践[M].北京∶ 科学出版社,1997∶ 178-186.LI Deji.Theory and practice on debris flow hazard reduction[M].Beijing∶ Science Press, 1997∶ 178-186.