喜马拉雅北麓郭如错冰湖泥石流沟发育特征及对拟建中国至尼泊尔铁路的影响研究

苗晓岐

(轨道交通工程信息化国家重点实验室(中铁一院)，西安 710043)

引言

中国至尼泊尔铁路起于西藏自治区日喀则市，终于尼泊尔首都加德满都，全长580 km。首次穿越喜马拉雅山，是全球冰湖分布最为集中的区域。喜马拉雅北麓对中尼铁路影响较大的冰湖12处，冰湖面积大小不等，其中最为典型的是郭如错冰湖，中国至尼泊尔铁路桥梁工程位于其溃决通道上(图1)。

图1 喜马拉雅北麓冰湖与中国至尼泊尔铁路位置关系(红色为郭如错)

据统计，仅喜马拉雅山在近50年来发生过20余次较大的冰湖溃决灾害，其中2/3发生在我国境内[1]，而目前冰湖溃决型泥石流特征值计算还处在研究阶段，指标选取不够合理，数值模拟亦不够精确。国内外标准、规范未见针对线性工程的溃决特征值计算，加之冰湖溃决后形成的洪水或泥石流具有爆发性强、规模大，危害范围广等特点，对溃决影响范围内的线性工程危害极大，成为制约中尼铁路选线和工程设置的主要工程地质问题。以喜马拉雅北麓郭如错冰湖为例，通过孕灾环境的地质调查，结合1984年至今多期历史卫星遥感影像、SBAS-InSAR形变分析、雨洪法、RAMMS数值模拟等，计算冰湖溃决后形成泥石流的运动学参数，并将参数进行对比分析，最后选取合理的参数，给出桥梁建议高度。为铁路、公路等穿越冰湖发育区，如何开展冰湖溃决型泥石流特征值计算工作进行了初步探索，在中尼、中巴、新藏等交通廊道选线上具有较高的应用和推广价值。

1 郭如错冰湖泥石流沟沟域基本特征

郭如错冰湖泥石流沟位于聂拉木县波绒乡，沟口地理坐标：N28°31′12.87″，E85°39′12.97″。中尼铁路线路方案从该泥石流沟口堆积区下游18.82 km以桥梁形式通过。郭如错冰湖泥石流沟流域面积约85.16 km2，主沟纵长19.70 km，流域范围内最高海拔约6 640 m，最低海拔约4 980 m，相对高差1 660 m，主沟平均纵坡降约为84.26‰，主沟总体顺直，主沟内发育郭如错1处大型冰湖。

2 孕灾地质环境条件

2.1 地形地貌

郭如错冰湖泥石流沟流域内呈南西高北东低，沟道呈宽浅“U”形，谷坡20°～60°，谷宽450～900 m。冰湖两岸岸坡30°～45°，由下而上为侧碛垄、坡积物，侧碛垄平缓处高出湖面70 m左右。沟道后缘为2条母冰川，右侧冰川坡度20°～25°，坡向10°～15°，左侧冰川坡度15°～30°，坡向90°～100°。冰湖下游主沟呈“V”形谷，主沟地势较平缓，早期形成的侧碛垄分布在主沟两侧，坡度35°～45°，高70 m左右[2]。

2.2 地层岩性

郭如错冰湖泥石流沟沟内被第四系滑坡堆积、泥石流堆积、冲洪积、湖积、冰碛堆积等成因的粉土、粉质黏土、砂类土和碎石类土覆盖，沟域后缘及侧缘出露二叠系片麻岩、石英片岩。

2.3 地震动参数区划

根据GB 18306-2015《中国地震动参数区划图》，泥石流沟流域抗震设防烈度为Ⅶ度，地震动峰值加速度0.15g，反应谱特征周期0.45 s。

2.4 郭如错冰湖基本特征

郭如错冰湖位于佩枯错东南侧上游，喜马拉雅山北坡，下游湖口地理坐标：N28°30′59.1″，E85°39′23.8″，湖口距离郭如错冰湖泥石流沟沟口约8.5 km。冰湖形似长条状，沿轴向长约6.50 km，平均宽约750 m，轴向40°，湖面未见浮冰，湖水较清澈。湖面靠近出水口有一平面似椭圆状的小岛，长50 m左右，宽约30 m，将湖水隔开分成两股水流，在小岛前段(近跌水处)3 m处汇成一股，进入跌水段。两岸靠近冰川处有大量基岩出露，裂隙发育，其下有风化碎片。该湖母冰川冰舌与湖相连，冰舌宽700～900 m，厚10～15 m[2](图2～图5)。

图2 郭如错冰湖卫星俯视图

图3 郭如错冰湖侧碛垄及后缘冰川

图4 郭如错冰湖出水口

图5 郭如错冰湖终碛垄堆积物

2.5 终碛垄特征

对郭如错冰湖终碛垄调查发现，该终碛垄以冰碛孤、块石为主，底部为二叠系片麻岩和石英片岩。终碛垄长700～900 m，宽350 m，高50～60 m，物质组成为块石70%，碎石25%，砂3%，粉土和粉质黏土 2%。出水口底宽30 m，顶宽100 m，高约30 m，水流方向40°，溃水口沟槽纵比降200‰，背水坡度10°[2]。根据图2及图6，目前终碛垄估算总体积约1.54×107m3。

图6 郭如错冰湖终碛垄纵剖面5-5′

3 物源分布特征

根据野外调查和遥感解译，郭如错冰湖泥石流主沟流域受冰川运动和常年积雪影响，泥石流沟域内物源储量丰富，松散物源类型为冰碛物、冻融物源、沟道物源、坡面物源、崩滑物源(图7)，主要集中在沟道内及冰川处，并以冰碛物及沟道物源为主，其中沟域内可启动物源共21处，具体分布及规模见表1和图8。

图7 郭如错冰湖沟道、冰碛、崩滑、冻融物源

图8 线路方案及郭如错冰湖物源分布示意

表1 郭如错冰湖泥石流沟物源统计

4 郭如错冰湖溃决危险性评价

4.1 库容和面积变化趋势分析

通过对郭如错冰湖1984年～2021年影像(来源GF-2、Landsat-8、OLI陆地成像仪、TIRS热红外传感器)的分析解译，得到郭如错冰湖历年湖水面积，由于受云量影响，可用影像较少，最终按照年份顺序筛选出8组影像[3-6]，见图9。

图9 郭如错冰湖1984年～2017年范围解译(红线)

统计结果显示，从1984年～2011年，郭如错冰湖面积变化较大(表2)，湖水面积4.531～4.833 km2，面积扩大约6.7%，从2011年～2017年，郭如错冰湖面积变化较小，湖水面积为4.833～4.841 km2。在统计年份内，湖水面积最大的年份为2017年。

冰湖库容估算采用冰湖面积乘以冰湖平均深度，冰湖平均深度的估算依据郭如错冰湖泥石流沟主沟纵断面图进行估算。假设冰湖上游边界到终碛垄下游堆积边界的沟道为顺直沟道，在纵剖面上将其连接为一条直线，同时假设终碛垄上游及下游的纵坡基本一致，结合Google高程所切现状终碛垄及冰湖纵剖面，可以推测出冰湖在剖面上的平均水深以及终碛垄高度。根据纵剖面推测，郭如错冰湖在主沟纵断面(图10)上沿纵剖面最长线统计的面积为0.338 km2，长度为6.60 km，因此其平均深度约为51.22 m，由此估算的郭如错冰湖库容见表2。

表2 郭如错冰湖1984年～2017年面积、库容

图10 郭如错冰湖泥石流主沟纵断面

统计结果显示，从1984年～2017年，库容为2.321～2.480亿m3，在统计年份内，库容最大的年份为2017年。

4.2 SBAS-InSAR郭如错冰湖形变分析

SBAS-InSAR技术最早由Berardino P等提出，其主要通过设置一定的时空基线阈值进行多组影像的组合，从而生成一定数量的差分干涉图，之后针对干涉结果进行滤波、解缠操作，解缠时主要以位于研究区的某一稳定点进行解缠操作，最后基于像元逐个计算分析并最终获取时间序列上的位移量[7-12]。本文通过SBAS-InSAR技术进行郭如错冰湖形变分析。如图11分别展示了郭如错冰湖2个时间段的典型差分干涉图(图11(a)、图11(b))，InSAR年均形变速率结果(图11(c))，由差分干涉图以及年均形变速率结果可看出，冰湖并未发生明显形变，并绘制冰湖典型形变点的历史形变曲线(图11(d))。

图11 郭如错冰湖2个时间段典型差分干涉

4.3 终碛垄结构及稳定性评价

郭如错冰湖终碛垄为冰碛物形成的块碎石堆积体，终碛垄坝基为片麻岩和石英片岩，岩性坚硬，透水性差，整个终碛垄面坡坡度约15°，背坡坡度约10°，块石之间较密实、摩擦系数较高。郭如错冰湖常年处于稳定流水状态，在坝顶处终碛垄右侧溃口有水流出，溃口过水面底宽30 m，说明该坝体具有稳定的出水状态(图2)，且根据近37年对冰湖及其终碛垄影像分析，堰塞坝体无变形迹象，稳定性良好。

4.4 冰湖溃决危险性评价

冰湖溃决的机制一般分成两类:(1)终碛垄内部埋藏的冰核融化引起终碛垄管涌破坏；(2)冰川崩塌或者滑坡导致湖水位上涨并引起涌浪，造成终碛垄溢流型溃决。这两种机制中，后者发生的次数多、灾害大[13]。郭如错后缘发育现代冰川，顺山坡呈舌状分布，面积2.3 km2，冰舌与湖相连，冰舌宽700～900 m，厚10～15 m，前缘冰川坡度约23.6°(图12(a))。后缘山体受冰川强烈侵蚀，山壁陡立，高达1 000 m，终年积雪。在冰川上发育大量冰裂缝(图12(b)～图12(d))，裂缝宽度最大可达20 m，深入冰川内部，将冰川切割成不连续的块体。在气温急剧上升的情况下，冰川融水沿裂缝通道渗入冰川内部，对冰川滑动起到促进作用，大量冰体可能坠落湖中，使湖水急剧上升，产生冲击力，导致冰湖终碛垄瞬间溃决，产生冰湖溃决洪水，在冰碛物的参与下，最终形成冰湖溃决泥石流。

图12 冰舌、冰滑坡、冰崩、冰裂缝

本文取临界冰川厚度最大值h作为冰川滑动的判别依据，当冰湖后缘冰川平均厚度>h时，冰川发生滑动[13]。取冰川与基岩接触面内聚力c=100 kPa，内摩擦角φ=0°[14]。

(1)

式中，ρ为冰川密度，取900 kg/m3；g为重力加速度，取10 m/s2；θ为冰川坡角。

研究区冰川类型属于亚海洋型冰川[15]，因此危险冰川的平均厚度H计算公式如下[17]

H=-11.32+53.21F0.3

(2)

式中，H为危险冰川平均厚度，m；F为危险冰川面积，km2。计算参数选取和结果如表3所示。

表3 郭如错冰湖危险冰川平均厚度计算结果[14]

因此，郭如错冰湖危险冰川平均厚度57 m大于临界冰川厚度27.8 m，冰川可能发生滑动，稳定性较低，发生冰川崩塌或者滑坡导致湖水位上涨并引起涌浪[16]，造成终碛垄溢流型溃决可能性较大。

5 郭如错冰湖泥石流沟特征值计算[18]

郭如错泥石流沟汇水区域面积约为193.14 km2，主沟长38.00 km，平均纵坡降为53.42‰。取郭如错冰湖泥石流沟沟口和中尼铁路工程所在位置为计算断面，结合铁路工程全寿命周期，本文设计频率考虑100年一遇对泥石流特征值进行计算和分析。

5.1 郭如错冰湖泥石流流速

根据流体形态现场调查，郭如错冰湖泥石流沟为黏性泥石流，密度为1.621 g/cm3，考虑堵溃效应，沟床糙率系数1/n取10，黏性泥石流沟采用西藏古乡沟、东川蒋家沟、武都火烧沟的通用公式计算流速

(3)

式中，VC为泥石流流速，m/s；HC为泥石流平均泥深，m；IC为泥位纵坡降；1/n为沟床糙率系数。

根据公式(3)，计算结果见表4。

表4 郭如错冰湖泥石流沟口平均流速计算结果

5.2 郭如错冰湖泥石流冲击力计算

泥石流冲击力包括整体冲击力和单块最大冲击力，整体冲击力计算公式为

(4)

式中，F为泥石流整体冲积力，kN/m2；g为重力加速度，9.8 m/s2；α为受力面与泥石流冲压力方向的夹角，(°)；λ为受力体形状系数，方形为1.47，矩形为1.33；圆形、尖端、圆端形为1.00。

单块块石的最大撞击力按下式计算

(5)

式中，FS为单块巨石的撞击力，kN/m2；γ为动能折减系数，正面撞击时取0.3；α为受力面与泥石流撞击面撞击角，(°)；C1、C2为巨石与建筑物的弹性变形系数，若采用船筏与桥墩台的撞击系数，C1+C2=0.005；W为巨石重力，kN。冲击力计算均考虑最危险情况，即有关参数选取为冲击力最大的数值，并只考虑正面撞击。泥石流整体冲击力和单块最大冲击力，根据公式(4)和公式(5)计算郭如错冰湖沟口断面和工程所在位置断面的整体和单块最大冲击力，结果见表5。

表5 泥石流冲击力计算 kN/m2

5.3 郭如错冰湖泥石流最大冲刷深度

根据公式(6)计算郭如错冰湖泥石流沟在沟口和工程所在位置断面处局部最大冲刷深度，相关参数取值和计算结果见表6。

表6 沟口泥石流局部最大冲刷深度计算

(6)

式中，HB为局部最大冲刷深度，m；P为冲刷系数；HC为泥石流泥深，m；VC为泥石流流速，m/s；VH为土壤不冲刷流速，m/s；n为与堤岸平面形状有关的系数，一般取值1/4～1/2；K为泥石流平均流速增大系数，根据内插法确定。

5.4 郭如错冰湖泥石流沟最大冲起高度和爬高

根据公式(7)和公式(8)计算郭如错冰湖泥石流沟在沟口和工程所在位置断面处最大冲起高度和爬高，相关参数取值和计算结果见表7。泥石流最大冲起高度

表7 沟口泥石流最大冲起高度计算

(7)

泥石流在爬高过程中由于受到沟床阻力的影响，其爬高

(8)

式中，b为迎面坡度的函数，取1.6。

6 郭如错冰湖溃决全过程数值模拟

为对雨洪法进行验证，同时进一步摸清郭如错冰湖溃决过程的运动特征，采用瑞士RAMMS(Rapid Mass Movements Simulation)三维数值模型[19-21]，模拟冰湖溃决全过程，计算终碛垄溃决后沟域不同断面的特征值。此次模拟，选择佩枯错至郭如错冰湖作为RAMMS模拟计算范围，参数选取同雨洪法，对冰湖面积4.841 km2(2017年遥感解译数据)的工况进行RAMMS模拟，获取其影响范围和冲击力等参数。

根据模拟结果，将郭如错冰湖溃决后洪水运动过程划分为4个阶段(图13)，具体包括：(1)开始运动阶段：终碛垄在地震或降雨条件下发生破坏，强度降低，郭如错冰湖开始发生溃决；(2)加速运动阶段：随着洪水流量的增加，洪水携带沟道内泥沙、碎石向沟道下游运动；(3)减速运动阶段：冰湖下游纵坡降仅有20.3‰，坡度很缓，沟道宽阔，动能变小，向下游运动，流速慢慢减小；(4)结束阶段：模拟结果显示，溃决洪水总量约2.48×108m3。

图13 郭如错冰湖溃决运动过程模拟

根据数值模拟结果，可知溃决洪水影响范围(图14)，拟建中尼铁路路线通过冰湖溃决洪峰流经区域，线路附近区域流速最大值为3.45 m/s，平均0.84 m/s，流深最大值为1.28 m，平均值为0.31 m，线路处冲击力最大值为19.00 kN/m2，平均值为1.53 kN/m2，可见溃决洪水对拟建铁路有影响。

图14 冰湖溃决流速、泥深、最大冲击力空间分布

7 郭如错冰湖泥石流特征值对比分析

通过雨洪法和数值模拟，计算郭如错冰湖全溃后铁路桥梁位置的流速、泥深、对桥墩整体冲击力如表8所示。

表8 铁路桥梁位置流速、泥深、对桥墩整体冲击力计算结果汇总

通过将雨洪法与RAMMS数值模拟方法进行对比可知，在铁路桥梁位置断面处，雨洪法计算的100年一遇降雨条件下的断面流速及泥深与数值模拟结果接近，说明将冰湖堵塞溃决看成沟道堵塞溃决的一部分，采用雨洪法计算的结果具备合理性。雨洪法计算的泥石流冲击力相对于数值模拟结果偏大，但是为了安全起见，以更为保守的雨洪法计算结果为准。

8 郭如错冰湖溃决对拟建线路影响评价

雨洪法和RAMMS数值模拟计算结果均表明，郭如错冰湖溃决后对铁路桥梁有较大影响，根据DZ/T0220—2006《泥石流灾害治理工程勘查规范》及在拉日铁路、青藏铁路泥石流方面的研究经验，流域面积≥50 km2的泥石流沟的建议安全高度为5 m。因此，对于铁路工程所在位置断面桥跨净空的建议高度采用泥深、冲起高度、爬高、安全高度之和7.89 m，可满足溃决后的净空需求。此外，应加强墩台的防护，减轻洪水裹挟块石等的冲击、掏蚀。

9 结论

(1)喜马拉雅北麓郭如错冰湖泥石流沟位于聂拉木县波绒乡，中尼铁路线路方案从该泥石流沟口堆积区下游18.82 km以桥形式通过，郭如错冰湖终碛垄估算总体积约1.54×107m3，沟内以冰碛物及沟道物源为主，可启动物源共21处，总静储量约28 049.82万m3，总动储量约3640.34万m3。

(2)通过1984年～2017年高精度遥感解译，郭如错冰湖面积变化较大，湖水面积为4.531～4.833 km2，面积扩大约6.7%，从2011年～2017年，郭如错冰湖面积变化较小，湖水面积为4.833～4.841 km2，平均深度51.22 m，库容为2.321～2.480亿m3，在统计年份内，湖水面积和库容最大的年份为2017年。

(3)对郭如错冰湖SBAS-InSAR典型时间段的差分干涉图、InSAR年均形变速率结果图、冰湖典型形变点的历史形变曲线进行分析，可以看出冰湖并未发生明显形变。

(4)郭如错冰湖有稳定的出水口，常年处于稳定流水状态，根据近37年对冰湖及其终碛垄影像分析，终碛垄无变形迹象，稳定性良好。

(5)郭如错冰湖危险冰川平均厚度57 m大于临界冰川厚度27.8 m，冰川可能发生滑动，稳定性较低，发生冰川崩塌或者滑坡导致湖水位上涨并引起涌浪，造成终碛垄溢流型溃决可能性较大。

(6)首次针对线形工程，计算郭如错冰湖泥石流沟铁路桥梁位置特征值计算结果，雨洪法：流速Vc=3.51 m/s，泥深Hc=1.26 m，整体冲击力F=29.30 kN/m2，单块最大冲击力FS=33.91 kN/m2，局部最大冲刷深度HB=3.92 m，泥石流最大冲起高度ΔH=0.63 m，爬高ΔHc=1.0 m；RAMMS数值模拟：流速Vc=3.54 m/s，泥深Hc=1.28 m，整体冲击力F=19.00 kN/m2。首次将雨洪法与RAMMS数值模拟方法进行计算结果进行对比可知，在铁路桥梁位置断面处，雨洪法计算的100年一遇降雨条件下的断面流速及泥深与数值模拟结果接近，证明采用雨洪法计算的结果具备合理性。雨洪法计算的泥石流冲击力相对于数值模拟结果偏大，但是为了工程安全，以更为保守的雨洪法计算结果为准。

(7)郭如错冰湖溃决对铁路桥梁的净空和墩台防护影响较大，对于铁路工程所在位置断面桥跨净空的建议高度采用7.89 m，并加强墩台防护，可满足工程安全需求。