基于FLO-2D的西藏若如村泥石流危险性分析

张奋翔， 张路青， 周 剑， 王 颂

(1.中国科学院 地质与地球物理研究所 页岩气与地质工程重点实验室，北京 100029； 2.中国科学院 地球科学研究院， 北京 100029； 3.中国科学院大学， 北京 100049)

1 研究背景

泥石流是山区沟谷或坡面在暴雨、洪水等自然或人为因素作用下发生的一种挟带大量泥砂、石块或巨砾等固体物质的特殊洪流，会对人类生命财产或生存环境造成较大危害[1]。青藏高原地区地貌多变，峡谷、河谷、宽谷等交替出现，在喜马拉雅块体、拉萨块体和羌塘块体的相互作用下，新构造运动活跃[2]。其中，青藏铁路、拉日铁路、G318国道、G109国道等交通工程沿线风化剧烈，岩土体完整性较差，发育着大量泥石流，对交通安全及线路运维造成了极大威胁。因此，对此区域泥石流动力过程的研究十分必要。

随着泥石流动力模型研究的深入以及计算机技术的发展和完善，数值模拟已逐渐成为泥石流动力过程研究的重要手段[3]。利用已有的或新建动力模型，选择合适的数值方法进行计算分析是泥石流动力过程数值模拟的基础。在数值算法和本构模型的研究相对成熟后，出现了描述泥石流运动的数值模拟软件[4]。Debris-2D，FLO-2D和Geoflow等是国际上应用比较成熟的数值模拟软件[5]，能够开展较大尺度的工程问题分析。其中，FLO-2D软件通过建立微分形式的质量守恒和动量守恒方程，采用显式中心差分法进行求解，可以计算流深、流速和影响范围随时间的变化。作为二维洪灾模型，FLO-2D的计算时间可以由用户控制，操作简单且实用，近年来应用广泛。Bertolo等[6]运用FLO-2D模拟了小流域泥石流运动淤积过程；阮德修等[7]将FLO-2D与3DMine耦合模拟了尾矿库溃坝灾害下的泥石流流动过程；龚柯等[8]利用FLO-2D对四川省汶川县绵虒镇簇头沟进行了泥石流危险性评价。

西藏地区由于地理位置偏远、人口居住密度小、人类活动工程少等原因，泥石流灾害的研究程度较低，且仅限于拉萨、日喀则、林芝等地区。戚国庆等[9]通过建立降雨型泥石流预测模型，对贡觉县的3条泥石流沟进行了灾害分析；韦方强等[10]进行了西藏古乡沟泥石流的数值模拟；刘伟朋等[11]研究了拉日铁路沿线(年木乡-日喀则段)泥石流的分布特征。

若如村泥石流位于拉萨市堆龙德庆区德庆乡，青藏公路、青藏铁路及在建的青藏高速公路位于该泥石流堆积扇的前缘，绕避或下穿泥石流扇。2016年7月2日，德庆乡暴发降雨频率为5%的持续性暴雨，当日20时至次日20时累计降雨量为26.6 mm，计算最大雨强为14.09 mm/h，此次降雨使得原若如村泥石流的堆积范围进一步扩大，对上述工程区域产生极大威胁。在详细野外调查的基础上，采用FLO-2D软件对青藏铁路、青藏高速公路及若如村的泥石流危险性进行了研究。

2 若如村泥石流概况

若如村泥石流距离拉萨市中心60 km，地处当雄至拉萨峡谷出口段。此地段属高原寒温带湿润季风气候，年平均降水量近460 mm；堆龙曲从区内经过，流域面积约5 100 km2；处于雅鲁藏布江缝合线与班公错-怒江缝合线之间，岩性主要是花岗闪长岩和以砂岩为代表的坚硬陆源碎屑岩。G109公路沿线南侧边坡较高，平均坡度大于30°，在风化侵蚀和降雨的综合作用下，坡面岩体破碎，有大量松散体堆积，因此，泥石流的发育具备良好的动力和物源条件。除了若如村泥石流，区段内还发育有近15条规模较大的泥石流。

若如村泥石流流域面积约3.2 km2，主沟道长2.5 km，最高海拔5 440 m，相对高差1 330 m，两侧山坡的平均坡度为30°，主沟平均坡降30.6%。坡体上部有大量松散堆积体，冲沟内也有碎屑物堆积，为泥石流的发育提供了充足的物源。先前泥石流活动在若如村口形成了一个较大的泥石流堆积扇，堆积扇顶角约40°，扇长350 m，扇宽410 m，面积约0.07 km2。若如村泥石流概况见图1。

图1 若如村泥石流概况

在泥石流沟口处，青藏铁路设置涵洞穿过堆积体，远离堆积体后铁路建立高架桥；堆积体前缘和青藏铁路之间建有三级重力式挡土墙来保障行车安全，墙体高15 m左右，基底宽度为6～7 m；青藏高速公路目前正在修建中，堆积体前方为明洞开挖区(以下简称明挖区)，开挖深度近5 m；沿着泥石流沟口至高速路施工区设置主排导槽，槽宽1.5 m，深约0.5 m。该泥石流的主要威胁对象为青藏铁路、青藏高速公路和若如村，此外，主排导槽在泥石流流动过程中可能面临淤积失效的风险。

3 数值模型构建

3.1 高精度DEM网格划分

根据FLO-2D泥石流运动数值模拟的方法和理论[8]建立模型。为了进行本次研究，在若如村泥石流堆积区前缘飞行无人机获得了泥石流流域的航拍照片。将航拍照片的高程和GPS信息提取后，获得流域的DEM数据进而转化为FLO-2D软件能够识别的ASCII文件。主要步骤可分为：

(1)获取DEM。将航拍照片利用Agisoft PhotoScan Professional软件进行处理，得到流域的DEM数据。

(2)ASCII文件转化。运用Uedit64中工具箱的转换功能进行文件格式的转换，得到FLO-2D建模所需的DTM数据。

(3)在FLO-2D中新建项目，将带有高程数据的DTM文件加载进去作为基础地形并进行网格划分。选取的网格尺寸为2 m，网格数量146 972个。计算网格如图2所示，划分出计算区域，最后对选定的计算区域进行高程插值。模型的入水口设置在泥石流沟上游的物源区，出水边界设置为堆龙曲河道。

图2 FLO-2D计算网格

3.2 模型参数取值

FLO-2D数值模拟选取的主要参数有[12]：降雨强度i，等效曼宁系数n，泥石流土石材料比重Gs和体积浓度CV，层流阻滞系数K，屈服应力τ及黏滞系数η，其中τ、η用相关系数α1、β1、α2、β2表征。

η=α1eβ1CV

(1)

τ=α2eβ2CV

(2)

根据现场情况，将泥石流物源分为坡体上部和冲沟内两部分。实地调查得到坡体上部松散体平均厚度约0.25 m，冲沟内碎屑物的平均厚度约0.15 m，结合卫星图进一步得出：坡体上部松散体面积约1.7×105m2，冲沟内碎屑物面积约3×105m2。日降雨量按照实际工况26.6 mm计算，该泥石流流域内的集水量为8.5×104m3。按照以上数据得到泥石流土石材料的体积浓度为50.6%，取0.51。根据 FLO-2D 手册结合现场情况，若如村泥石流土石材料的比重取值为2。

由土石材料重度和体积浓度，划分为中阻泥石流，则泥砂比Rns=0.75，现场调查若如村泥石流的堆积扇厚度约5 m。根据王裕宜等[13]提出的不同泥砂比条件下的统一公式：

(3)

τ=0.021Rns-0.3exp(19.64Rns0.25·CV)

(4)

η=1.39×10-4Rns-0.47exp(18.07Rns0.13·CV)

(5)

式中：h为泥深，m。将h=5，Rns=0.75代入公式(3)，得到流动区域曼宁系数n=0.09；将公式(1)、(2)、(4)、(5)联立，得：

α1=1.39×10-4Rns-0.47,β1=18.07Rns0.13

(6)

α2=0.021Rns-0.3,β2=19.64Rns0.25

(7)

将Rns=0.75代入公式(6)、(7)，得到α1=0.000159，β1=17.41，α2=0.0229，β2=18.28。

层流阻滞系数K的取值范围为 24～50 000，鉴于土石材料的体积浓度值处于0.48～0.55，且具有黏性泥石流的特点，参考FLO-2D手册中对于层流阻滞系数K的建议，取K=2 285。模拟参数取值见表1。

表1 FLO-2D输入参数值

为全面分析流域内的危险性，以德庆乡1998-2018年的降雨数据为基础，分别对降雨频率为5%、2%和1%下的若如村泥石流进行模拟计算，并以上述频率设计泥石流流量过程线。

根据水文地质手册，对降雨数据分析计算得到2%和1%频率下的日降雨强度分别为48.1mm和57.8mm，采用雨洪法计算出该泥石流在不同降雨频率下集水点的流量，结果见表 2。

由于 FLO-2D数值模拟时无法考虑泥石流对沟道的侵蚀效应，因此在清水流量的基础上乘以流量放大系数，得到泥石流流量再生成流量过程线输入到FLO-2D中：

(8)

式中：BF为流量放大系数。

表2 若如村泥石流不同降雨频率下的集水点流量

根据表2数据，采用典型洪水过程放大法求取泥石流流量过程线[14]如图3。

4 模拟结果验证

由于实际24 h降雨量达到20年一遇的级别，将5%降雨频率下的模拟结果与实际测量结果进行对比，对模拟结果的准确性进行评价。采用堆积范

围、沟口最大冲出距离与横向最大堆积宽度3个指标，将现场调查与遥感影像结合得到的实际堆积数据与模拟情况进行对比，见表3。根据表3数据，计算泥石流数值模拟精度因子A，公式如下：

(9)

式中：S1为实际堆积面积，m2；S2为模拟堆积面积，m2；S0为模拟与实际堆积范围重叠的面积，m2。

图3 不同降雨频率下泥石流流量概化过程线

表3 模拟堆积与实际堆积数据对比

将表3数据代入公式(9)进行计算，得到此次模拟的精度因子为68.7%。通过野外调查发现，高速公路明挖坑周边及隧道口、青藏铁路路基下方都有残留的泥石流松散体，若如村西北角房屋前方仍可见0.5 m厚的土石堆积体，与模拟结果基本相符。由于沟口最大冲出距离的模拟误差率为19.4%，横向最大堆积宽度的模拟误差率为8.4%，与精度因子结合来看，此次模拟结果比较符合实际。

5 模拟结果分析

该泥石流的危险性分析主要包括两方面：(1)主排导槽在不同降雨强度下的导流情况分析；(2)堆积区下穿交通工程和若如村在不同降雨频率下的危险性分析。主排导槽从泥石流沟口向下延伸，经过铁路涵洞上方，将冲沟上游的洪水引导至高速公路隧道明挖区。主排导槽深度约0.5 m，当其边界的泥石流流深大于0.5 m时，流体会从导槽内溢出，此时主排导槽失效。下面将结合不同日降雨强度下泥石流流深、流动速度以及冲击力，对主要威胁区域的危险性进行分析。

分别选取各威胁区域内最易产生堆积的网格记录不同降雨工况下泥石流的流动情况，各威胁区域代表性点的具体位置如表4所示，取点位置如图4。

表4 各威胁区域危险性分析的选点网格单元及位置

5.1 泥石流的流动深度分析

图5给出了3种降雨频率下模拟得到的泥石流最大流深。从堆积范围来看，5%降雨频率下的堆积范围约15.8×104m2，2%频率下约17×104m2，1%频率下为18×104m2。在5%、2%和1% 3种降雨频率下，青藏铁路涵洞西出口及其以西铁路段、铁路涵洞、高速公路隧道明挖工程区以及若如村的大部分区域都被泥石流所覆盖。铁路涵洞东出口、高速路东侧工程区以及河漫滩在5%降雨频率下未见泥石流覆盖，在2%降雨频率下开始有流体堆积；1%频率下的泥石流堆积范围更广，铁路涵洞东出口完全被覆盖。

图6给出了上述各代表性点在3种不同降雨频率下泥石流流深随时间的变化曲线。结合图5、6分析如下：

(1)在5%、2%和1%降雨频率下，铁路涵洞西出口的流深分别为1.4、1.9、2.1 m；与涵洞西出口相比，涵洞东出口受泥石流影响较小，在5%和2%降雨频率下几乎无流体堆积，在1%降雨频率下的流深小于0.1 m。根据铁路行车规范，当流体覆盖道床时铁路将无法通行，由此判断5%降雨频率下的泥石流流深就会影响铁路正常通车。

(2)主排导槽左边界在5%、2%和1%降雨频率下的流深分别为2.39、2.74和2.91 m，导槽深度仅有0.5 m，由此判断在5%降雨频率下主排导槽就会失效。从流深曲线来看，主排导槽在最初10 min里的流深小于0.5 m，因此，当泥石流流量较小时，主排导槽可以较好地发挥排导作用。

(3)由于地势较低，高速公路明挖点成为各个方向的汇水区，其流深变化曲线的规律性不强。5%降雨频率下的流深约5.4 m，而明挖坑深5 m，处于溢满状态，对施工人员安全和工程稳定有显著影响；2%和1%降雨频率下的流深超过6 m，对明挖区域的安全和稳定影响更大。

(4)在5%降雨频率下，若如村西北角流深的最大值为0.67 m，村民的居住和通行将受到严重影响；2%和1%降雨频率下的流深超过1 m，居民安全将面临威胁。

(5)随着时间的变化，各代表性点的流深都经历了无堆积-突然增大-逐渐减小的过程，符合泥石流突然暴发和逐渐消退的特点。

5.2 泥石流流动速度与冲击力分析

图7给出了不同降雨频率下各代表性点模拟过程中的泥石流流速变化，由图7可以得出：

(1)5%降雨频率下，泥石流在铁路涵洞西出口的最大流速为0.52 m/s；在铁路涵洞东出口的流速为0，这与流深的变化相对应；主排导槽左边界的最大流速达到1.35 m/s，为工程区域内最大流速点，再次验证了主排导槽的导流功能；与流深变化曲线类似，高速路明挖区的流速变化呈无规律曲线，可以得出其最大流速为0.65 m/s；若如村西北角房屋边界的最大流速为0.1 m/s，威胁相对较小。

(2)与流深变化相同，各代表性点的模拟流速基本呈现不流动-突然增大-逐渐减小的过程；与流深变化不同的是，各点最大流速并非随着日降雨强度的提高而增大，如主排导槽左边界的最大流速在1%降雨频率下为1.44 m/s，小于2%降雨频率下的1.59 m/s，而高速路明挖区在5%降雨频率下的最大流速也大于2%降雨频率，说明流深增大产生的淤积效应可能引起泥石流流速的降低。

图8为不同降雨频率下各代表性点模拟过程中泥石流冲击力变化。根据图8可以得出，在5%降雨频率下，铁路涵洞西出口的最大泥石流冲击力约1.1 kPa，东出口冲击力为0；主排导槽左边界在泥石流运移过程中的最大冲击力为16.29 kPa；高速路明挖点的最大冲击力为7.51 kPa；若如村西北角房屋承受的最大冲击力仅25.06 Pa。

5.3 其他极端降雨条件下的泥石流危险性分析

为了更加全面地分析各威胁区域在极端降雨天气下的泥石流危险性，以研究区域近20年的年最大24 h降雨量数据为基础，选取其中的最小值19.5 mm作为极端降雨区间的下限，将1%降雨频率下的24 h降雨量57.8 mm作为上限，对降雨量每隔5 mm进行插值，进行了19.5、24.5、29.5、34.5、39.5、44.5、49.5、54.5、57.8 mm共9种日降雨强度的模拟。将模拟结果进行处理，得到各威胁区域的最大流深和最终流深随降雨量的变化曲线，如图9所示。

由图9可得，当日降雨强度在19.5～57.8 mm变化时：

(1)泥石流在铁路涵洞西出口的最大瞬时流深大于1 m，在日降雨强度为54.5 mm时可达2.2 m，根据铁路行车规范，列车在某个时间段内无法安全通行；然而，从最终流深变化曲线分析，涵洞西出口的泥石流流深会逐渐消退至0.4 m以下，对列车通行安全影响有所减小，19.5～34.5 mm日降雨强度下的最终流深在0.2m左右，对列车通行几乎无影响。铁路涵洞东出口的泥石流最大流深几乎为零，对列车通行无影响。

(2)主排导槽左边界的最大瞬时流深大于2m，在日降雨强度为54.5 mm时接近3 m；不同日降雨强度下主排导槽左边界的最终流深在0.7～1.1 m之间变化，流深始终大于主排导槽深度，故判断为失效。

图6 3种降雨频率下各代表性点的泥石流流深变化曲线

图7 3种降雨频率下各代表性点的泥石流流动速度变化曲线

(3)当日降雨强度大于19.5 mm时，明挖坑的瞬时流深大于5 m，此时坑内已被泥石流淤满，施工安全及工程稳定将会受到严重影响；从最终流深来看，坑内泥石流流深会逐渐消退至2～4 m，可能延误工期。

(4)日降雨强度在19.5～29.5 mm变化时，若如村内的最大瞬时流深小于0.8 m，可能影响居民出行；当日降雨强度大于29.5 mm时，村内的最大瞬时流深大于1m，可能危及居民的生命财产安全。从最终流深变化曲线分析，若如村的泥石流流深会逐渐消退至0.4 m以下，此时仍可能影响居民交通，对生命财产安全几乎无影响。

(5)随着日降雨强度逐渐增大，各威胁区域的最大瞬时流深均呈现先增大后减小的趋势，在日降雨强度为54.5 mm时达到最大值，说明此降雨强度下的若如村泥石流危险性最强。

图8 3种降雨频率下各代表性点的泥石流冲击力变化曲线

图9 各威胁区域流深随日降雨量的变化曲线

6 人类活动区的泥石流危险性评价

基于上述模拟获得的泥石流流深与流速，对人类活动区整体的危险性进行评价。为了使评价结果更具代表性，将1%、2%和5% 3种降雨频率进行综合分析。根据FLO-2D手册，泥石流的危险性与灾害程度及暴发频率相关，评价标准如图10。泥石流的灾害程度受流体的流深和流动速度控制，采用汶川泥石流强度划分标准[8]，如表5所示。

图10 泥石流危险性评价标准

各颜色表示意义如下：透明色(无标注色)为无危险区，表征安全区域；黄色为低危险区，表示该区域内，泥石流可导致工程建筑物轻微损伤；橙色为中危险区，表示该区域内工程建筑物部分摧毁，可能出现人员伤亡；红色为高危险区，表示泥石流会导致该区域内工程建筑摧毁、人员伤亡惨重。从透明色到红色依次分级为0～3。根据上述标准，得到若如村泥石流危险性评价图见图11。

表5 泥石流灾害程度评价标准

图11 若如村泥石流危险性评价图(人类活动区)

从图11来看，高危险区主要分布在主排导槽、涵洞西出口以及高速路明挖区，此区域内包括居住、铁路通行以及高速公路施工在内的一系列人类活动将面临极大的安全威胁，这与前文分析一致；铁路涵洞、高速公路隧道以北的工程区以及若如村的小部分区域处在中危险区，可能出现工程设施摧毁和人员伤亡；涵洞东出口及其以东的青藏铁路段、高速公路隧道以东的工程区基本无危险。

7 结 论

(1)在实际降雨频率5%的情况下，泥石流堆积面积约1.58×105m2，最大流深为5.7 m，最大流速13.1 m/s，人类活动区域会出现铁路无法正常通车、主排导槽失效和高速路明挖区溢满等问题，若如村村民的住行安全也将受到影响。

(2)当日降雨强度在19.5～57.8 mm变化时：青藏铁路若如村段在降雨期间无法通车，雨后列车通行不受影响；主排导槽始终失效；高速公路隧道明挖坑在降雨期间达到淤满状态，雨后淤积深度2～4 m，可能延误工期；日降雨强度小于29.5 mm时，若如村民出行不便，大于29.5 mm时，泥石流会危及居民的生命财产安全；54.5 mm日降雨强度下的泥石流危险性最高，在防治工作中应高度重视。

(3)根据危险性分区评价图，高危险区主要分布在主排导槽、涵洞西出口以及高速路明挖区，铁路涵洞、高速路隧道以北工程区以及若如村的小部分区域处在中危险区，涵洞东出口及其以东铁路段、高速公路隧道以东工程区基本无危险。

(4)在中、高危险区，需要加强泥石流防治工程建设，建议在铁路涵洞西出口或若如村西北角布置降雨观察点，根据雨情采取相应的预防措施。此外，建议增加主排导槽坡降以增大流体流速，考虑采用速流结构。