岷江河谷深层倾倒形成演化特征及堵江灾害链效应

2022-11-04蒙明辉崔圣华张文刘康林杨学之覃亮蒋涛

科学技术与工程 2022年27期

蒙明辉，崔圣华，张文，刘康林，杨学之,2，覃亮,2，蒋涛

(1.四川省华地建设工程有限责任公司，成都 610081； 2.四川省地质灾害防治工程技术研究中心，成都 610081；3.成都理工大学地质灾害防治与地质环境保护国家重点实验室，成都 610059； 4.中国地质科学院探矿工艺研究所，成都 611734)

龙门山强震山区在经历多次地震后，诱发了大量地质灾害，其数量、规模、类型及造成的损失举世罕见[1-2]。岷江河谷地处龙门山强震区，特殊的深切峡谷地貌使得灾害往往发育高位且具有链式效应[3]，如茂县新磨村滑坡。同时，深切的峡谷也为反倾岩层的倾倒变形提供了良好的发育条件。

反倾斜坡因层状岩体倾倒引起斜坡变形的研究已较为系统，胡振辉等[4]基于不同地貌类型的边坡倾倒变形体数量、面积、体积与厚度分析对边坡倾倒变形的贡献率；Liu等[5]以二古溪斜坡为例，分析了在变质岩中深层倾倒变形的形成机制；Adhikary等[6]、Leandro等[7]对岩石边坡的弯曲倾倒破坏机理进行了研究；邢一飞等[8]结合典型实例，模拟了开挖后边坡倾倒变形机制及演化过程，并优化了现有防治措施；宁奕冰等[9]以澜沧江上游某深层倾倒体为研究对象，采用物理模拟和数值模拟，揭示了该类深层倾倒体在河流下切及水库蓄水作用下的失稳机制及演化过程；黄达等[10]基于离心模型试验，研究了坡角、倾角对边坡稳定及破坏模式的影响，进一步分析了反倾边坡变形机理；郑达等[11]基于离心模型试验分析了边坡倾倒变形弯折带的演化发育过程与特征。

上述学者从倾倒变形的机理、倾倒边坡的演化和稳定性方面均有大量研究。但倾倒变形体的折断带完全风化成“土”且转化成滑带的案例相对较少。

现在对茂县梯子槽倾倒采用地质测绘、InSAR、钻探、物探、井探等手段进行现场勘查的基础上，对其特征、成因机制做初步探讨。同时，针对倾倒可能造成的灾害链效应，采用水利上的溃坝洪水模型，分析预测溃坝洪水对下游造成的风险，为短时间内制定合理可行的应急处置方案提供有力的技术支撑。

1 梯子槽深层倾倒概况

区域构造上属于龙门山中央推覆构造带，地壳抬升强烈，地震频发，研究区周边发生Ms4.7级以上地震达60余次[11]，其中1933年叠溪7.5级地震、1976年松潘7.2级地震和2008年汶川8.0级地震影响较大(图1)。

图1 区域位置和历史地震Fig.1 Regional location and historical earthquakes

图2 梯子槽倾倒全貌Fig.2 Overall view of Tizicao landslide

倾倒区属于高山峡谷地貌(图2)，两岸谷坡陡峻，分水岭高程3 700～3 900 m，谷底高程约1 730 m，相对高程约2 000 m。梯子槽倾倒位于岷江右岸谷坡下部，倾倒前缘底部为基岩陡崖，坡度约70°，临空条件极好，倾倒区整体平均坡度约30°，倾倒顶部为缓坡平台，坡度20°～25°。坡表冲沟发育，据钻孔揭露，地下水活动微弱。

梯子槽倾倒平面形态呈不规则的四边形，后缘至滑动形成的深大拉张裂缝后侧陡壁；前缘剪出口位于海拔2 000 m左右，坡度由陡变缓的地形转折部位；北侧至大槽沟附近，南侧至老熊洞沟。

根据调查确定的倾倒周界，倾倒区平均纵向长度540 m，平均宽度约560 m，倾倒区平面面积30.78×104m2，据钻孔和探井揭露情况统计，倾倒平均厚度(含滑带)45.1 m，据此计算该倾倒总方量为1 388.2×104m3，其规模分级为特大型。

2 研究方法

为查清梯子槽倾倒的地质环境条件、结构特征、变形历史、成因机制及灾害链，采用了多种研究方法(图3)。

图3 研究方法示意图Fig.3 Schematic diagram of research methods

2.1 地质调查

收集了研究区相关的历史地震研究报告，基于ArcGIS软件利用1∶10 000地形图生成倾倒滑动前的数字高程模型(digital elevation model，DEM)，充分收集利用了研究区区域地质资料，获取了研究区地质构造、地层岩性等基本地质资料。共完成钻探23个孔1 256.7 m(图4)，井探4口132.3 m，槽探4条85.6 m3，同时结合电镜扫描等手段，充分查明了倾倒地质结构特征。

2.2 变形调查

利用无人机航拍1.4 km2并进行数据处理获得了高分辨率的倾倒滑后DEM模型，结合地面调查，查明了倾倒变形迹象分布位置；利用2013年9月—2016年7月四期期差分干涉探测变形(interferometric synthetic aperture radar，InSAR)，得到了倾倒历史变形过程。

图4 倾倒区各钻孔倾倒结构统计柱状图Fig.4 Statistical histogram of landslide structure of each borehole in the landslide area

2.3 灾害链分析

基于数值模拟分析，分析了倾倒稳定性状态并结合数值计算，获取了倾倒堵江的体积。利用成熟的数学模型谢任之公式研究计算了溃口流量，基于一维圣维南方程预测了溃决洪水下游演进过程和风险。

3 结果

3.1 倾倒地质结构特征

3.1.1 变形体特征

梯子槽倾倒由表层松散土体和下部碎裂岩块组成。松散土体主要为粉质黏土夹碎块石，厚度3～10 m，碎块石以千枚岩为主，碎块石(质量分数)占40%～50%，块度大小不均匀，一般2～15 cm，局部见直径30～50 cm的转石，以棱角状或次棱角状为主，无明显分选性和磨圆度特征。碎裂岩块分布于整个倾倒区，厚度10～50 m不等(图5)。

据勘查，梯子槽倾倒平均厚度为45.1 m，厚度的空间变化较大(图6)。总体上倾倒中部厚度较大，后缘和前缘相对较薄，南侧C区中前部较厚，后缘较薄的特点，B区中后部较厚，前缘较薄，可见倾倒厚度的变化较大。

图5 ZK12变形体碎裂特征Fig.5 Fragmentation characteristics of ZK12 deformed body

图6 倾倒厚度等值线图Fig.6 Contour map of landslide thickness

图7 ZK05滑动带特征Fig.7 Characteristics of ZK05 sliding zone

3.1.2 滑动带特征

根据钻孔和探井，滑带土主要为灰色或黄褐色粉质黏土夹角砾、碎块石(图7)。粉质黏土多呈硬塑状，稍湿-湿，密实；角砾、碎块石成分以千枚岩为主，直径多为2～20 mm，以棱角状和次棱角状为主。滑带土厚度2～5 m，局部厚度较大，ZK05滑带达到10.3 m。滑带土中常见镜面与擦痕等现象(图8)，电镜扫描(SEM)同样显示滑带土具有明显的擦痕(图9)。

图8 倾倒后缘探槽揭露滑带镜面与擦痕Fig.8 The sliding belt mirror surface and scratches are exposed by dumping the trailing edge

图9 电镜扫描揭示滑带擦痕Fig.9 Scanning electron microscope reveals slip band scratches detection groove

由滑面等高线图(图10)可知，滑动面倾向与坡向总体一致，但滑面略呈宽阔的槽状负地形，因此各区滑动方向有所不同，虽总体都是向东滑动，但A、B区滑动方向略向南偏，C区则略向北偏，这也是B区在前缘存在挤压，造成B1分区鼓胀解体变形区的重要原因。

3.1.3 基床特征

图10 滑动带等高线图Fig.10 Contour map of sliding zone

梯子槽倾倒滑床为泥盆系危关群上组(Dwg2)炭质千枚岩(图11)，滑床基岩风化程度总体上较倾倒碎裂岩体低，为中风化-弱风化岩体，岩芯以长柱状为主，探井揭露滑床基岩总体完整(图12)。

图11 钻孔揭露滑床千枚岩Fig.11 Drilling exposed phyllite of sliding bed

3.2 倾倒变形破坏特征

3.2.1 地表变形破坏

总体上倾倒区北侧变形量大于南侧，主要受倾倒区地形和滑面形态的控制。根据倾倒变形破坏强度、变形现状及趋势、成因的初步判断和分区边界特征，将倾倒区划分为北侧滑动变形区(A区)、中部变形区(B区)和南侧影响区(C区)等三个分区(图13)，各区特征如表1。

图13 梯子槽倾倒变形迹象及变形分区图Fig.13 Deformation signs and deformation zoning map of Ti Zicao landslide

(1)北侧滑动变形区(A区)。倾倒区北侧地形陡峻，前缘呈陡岩地貌，临空条件好，滑面形态总体平直，仅前缘稍缓，非常有利于倾倒的复活剪出，因而变形强烈。当滑动岩土体在重力推动下滑动，在前缘剪出后，直接以塌落方式掉落至陡岩下方，并牵引后缘出现LF01、LF02等弧形拉张裂缝，呈现目前A1分区的变形破坏特征；由于倾倒本身的滑移破坏，加之A1分区的塌落变形后，倾倒前缘失去支撑，牵引后缘强烈的拉裂沉陷变形加剧，从而发育了LF07和LF04等裂缝，构成A2变形分区；同时，由于A区的强烈滑动，与北侧壁岩体的摩擦，因而在北侧倾倒边界附近发育LF11的一系列羽状剪切裂缝，构成A3变形分区(图14)。

图12 探井揭露滑床千枚岩Fig.12 Phyllite of sliding bed exposed by exploration well

(2)中部变形区(B区)。倾倒区中部地形总体上也较陡峻，但前缘发育相对较长的缓坡台地或宽缓山脊，滑动面呈弧形，前缘抗滑体分布相对较长，倾倒稳定性好于北侧，变形相对较弱，形成了倾倒B区。由于倾倒的整体滑移，其中部和后缘出现LF03、LF05、LF06、LF07等一系列拉张裂缝及沿裂缝的下坐变形，形成B2分区；随着坡体中部和后缘的拉裂沉陷，前缘抗力体段受阻发生隆起，加之各区的差异滑动及滑动方向的不同造成B区前缘挤压和应力集中，纵向和横向两个方向的挤压隆起，导致该区出现表面的张裂解体现象，形成B1分区(图15)。

图14 A区配套变形特征Fig.14 Matching deformation characteristics of zone A

(3)南侧影响区(C区)。倾倒区南侧前缘抗滑体分布最长，倾倒稳定性最好，仅在A、B区滑动影响下出现一些变形迹象。由于倾倒的整体滑移，其中部和后缘出现LF07、LF08拉张裂缝及沿裂缝的下坐变形，形成C2分区；随坡体中部和后缘的拉裂沉陷，前缘抗力段受阻发生隆起，加之差异滑动和九道拐沟局部地形临空面控制，发育LF09和LF10两条纵张裂缝，形成C1分区(图16)。

3.2.2 内部变形破坏

1)地层倾倒迹象

通过调查发现岩体的倾倒变形迹象有两处，分别位于倾倒后缘及滑塌槽顶部出露岩体中。从图17、图18中可以看到，倾倒岩体发生了较为明显的倾倒弯曲迹象，表现为下部岩层陡立，上部岩层倾角有所变缓。

岩体倾倒变形主要受岩体结构控制，由于变形体整体为反倾上游斜向坡，岩性为千枚岩，性质软弱，且为薄层状结构，刚度小，故易导致岩体在重力作用下向上游方向发生倾倒弯曲。

图15 B区配套变形特征Fig.15 Matching deformation characteristics of zone B

图16 C区配套变形特征Fig.16 Matching deformation characteristics of Zone C

图17 前缘岩体倾倒变形迹象Fig.17 Evidence of toppling deformation of front rock mass

表1 梯子槽倾倒区分区特征统计表Table 1 Statistics of zoning characteristics of Tizicao landslide area

2)岩层方位与角度变化

在对倾倒周边20个岩性出露点调查发现(表2、图19)，梯子槽倾倒范围内基岩为泥盆系危关群上组(Dwg2)灰色千枚岩与灰色薄层石英岩互层，坡型为反倾上游斜向坡，岩层产状由低高程向高高程逐渐变缓，由前缘沿河侧的195°∠71°转变为后缘的213°∠37°(图20)。

图19 梯子槽倾倒周边岩层产状分布图Fig.19 Occurrence distribution of rock strata around Ti Zicao

根据对梯子槽倾倒各处产状调查结果可知，由低高程至高高程，倾向逐渐由190°转变为215°，倾角逐渐由75°转变为40°，甚至在坡体后缘倾角转变为35°，斜坡上游侧大槽沟地形较陡，坡向约40°，说明斜坡岩体存在倾倒变形。

图18 后缘岩体倾倒变形迹象Fig.18 Evidence of toppling deformation of trailing edge rock mass

在对梯子槽倾倒进行勘探时共布设了23个钻孔，钻孔深度35～68 m不等，均揭露了倾倒(上部粉质黏土夹碎块石，下部碎裂岩块)、滑带和滑床，完整显示了梯子槽坡体结构中的岩层倾角变化规律：以ZK01为例，随着钻孔深度增加，岩层倾角从23°逐渐整大到71°(图21)。

3.3 倾倒形成演化机制

3.3.1 地表变形演化历史

据现场调查访问，倾倒变形始于2013年，主要表现为中前部坡体和房屋的开裂下错；2014年9月倾倒前缘中部始出现垮塌，逐渐形成宽约60 m，高约200 m的塌落槽，堆积体进入岷江成为雍塞体，与此同时，倾倒后缘开始出现裂缝；2015年8—9月，倾倒变形明显加剧，裂缝数量随之增多、裂缝加宽、增长；2016年和2017年雨季，倾倒变形继续发展，至2017年7月，倾倒后缘、前缘局部下错已达10余米，前缘剪出口附近塌落体部位不断发生垮塌和落石。

表2 岩层产状调查点统计表Table 2 Statistical table of survey points of rock stratum occurrence

根据2013年9月—2016年7月差分干涉探测变形(InSAR)成果图(图22)，倾倒后缘和前缘干涉条纹呈现黄-紫色调，与周边区域的浅蓝色调区别明显，表明倾倒后缘及前缘在该时间内存在较为明显的变形，与现场调查访问的结果一致。

图20 各调查点岩层产状特征Fig.20 Occurrence characteristics of rock strata at each survey point

图21 钻孔岩心岩层倾角变化Fig.21 Variation of rock dip angle of borehole core

3.3.2 倾倒演化过程机制

梯子槽倾倒历经了从岷江的河谷下切、边坡岩体发生卸荷现象到岩体倾倒弯折变形折断的漫长的地质动力发展过程，这种河谷的地质改造作用与发育的薄层状软质千枚岩、重力作用和良好的临空条件有关，其变形发展阶段可概况为五个阶段(图23)。

(1)河谷下切、卸荷回弹，岩体初始变形阶段：岷江河谷下切，反倾坡外的层状千枚岩受到卸荷回弹及重力双重作用下产生弯矩，开始向临空面发生由表及里得初始倾倒变形迹象。

图22 2013年9月—2016年7月差分干涉结果Fig.22 Results of differential interference from September 2013 to July 2016

图23 梯子槽倾倒演化过程Fig.23 Evolution process of Tizicao landslide

(2)河谷继续下切，岩体倾倒-层内拉张发展阶段：当河谷继续下切，变形累积，倾倒变形的变形位移进一步加大，岩体底部的软弱部开始产生向坡外的剪切位移并最终起主导作用，岩体将会发生拉张破坏。

(3)岩体倾倒变形出现折断面阶段：岩体内部出现宏观倾倒弯曲变形，当弯曲段弯矩大于岩体的极限抗弯强度时，岩体将会被折断，逐步形成断续的折断面。同时，受软弱带剪切力的影响，岩体向坡外剪切位移增大，导致折断面逐渐连续贯通。

(4)倾倒变形继续发展出现折断带阶段：岩体倾倒和剪切变形进一步发展，沿着最大的弯折面逐步形成的折断带，折断带受剪切力作用下，岩体异常破碎。

(5)折断带风化成“土”-底部滑移阶段：折断带上部岩体在倾倒变形过程中出现众多裂隙，降水及地下水沿裂隙和岩层面可直达折断带，造成折断带风化加剧形成目前粉质黏土夹角砾的梯子槽倾倒的滑带。且倾倒地处强震山区，多次受诸如叠溪地震、汶川地震、九寨沟地震的损伤叠加，以“滑移-拉裂”的形式在变形破坏，即目前状态下的梯子槽倾倒。

3.4 倾倒体失稳-堵江灾害链分析

3.4.1 堵江体积

叠溪—茂县一带多为高山峡谷地貌，岷江在该段河谷狭窄，一般仅150～300 m，两岸的倾倒易堵塞河道，形成堰塞湖，历史上岷江多次发生倾倒、泥石流等地质灾害堵塞河道情况，造成重大损失[12]。

河谷宽度、河水流量，对倾倒堵江的成生有控制作用，它取决于河床条件、河水流量、倾倒入江体积[13]。采用堵江最小土石方方法评价梯子槽倾倒堵江可能性，估算梯子槽变形体堵江的最小方量的公式为

(1)

式(1)中：Vmin为完全堵江所需最小土石方量，m3；Hr为河水深度,m，取25 m；Br为河床宽度,m，取235 m；φs为堵江岩土体饱水状态下的内摩擦角(取30°)。

据此可估算梯子槽变形体所处河道完全堵江需要的最小土石方Vmin=90.1×104m3，前缘LF02下方A1区变形最大，安全系数(factor of safety，FOS)小于0.9(图24)，其估算体积约130×104m3，即前缘A1区下滑就可能造成岷江堰塞。由于该倾倒物质以碎石土为主，力学强度低，渗透性较强，形成的堰塞坝自身稳定和抗渗能力较差，易溃坝引起次生灾害。

图24 自重状态下倾倒安全系数分区图Fig.24 Zoning diagram of dumping safety factor under self weight

若梯子槽倾倒堵塞河道形成天然坝体，在岷江水流的冲刷及侵蚀下溃坝，将对下游的安全造成极大的威胁。分析堰塞湖的溃决过程，是降低溃坝洪水对下游区域灾害损失的有效途径，可为下游危险区域人员撤离以及制定应急处理措施提供较为科学合理的依据。基于此，取最不利状态，即倾倒A区和B区倾倒全部进入岷江形成堰塞坝瞬时全溃进行灾害链演进及风险分析。根据计算分析，下滑方量为478.7×104m3，形成堰塞坝高度22.1 m，回水长度1 227 m，堰塞体库容271.3×104m3。

堰塞湖溃决洪水灾害链演进及风险分析主要包括：①溃口流量计算；②溃决洪水下游演进过程预测；③下游洪水淹没风险分析[14-16]。

3.4.2 溃决洪水灾害链演进

1)溃口流量计算

溃口流量的大小与溃坝前上下游水深和坝址断面形状及尺寸有关。采用谢任之公式计算溃口流量，即

(2)

式(2)中：qm为溃口最大流量，m3/s；λ为流量参数；B0为溃口顶宽度，m；H0为发生溃坝时坝前水深，m；g为重力加速度,m/s2。

采用四次抛物线法对溃口坝址流量过程线计算。瞬时全溃时，堰塞坝溃口顶宽度162.3 m，起溃水位1 752.7 m，得到溃口流量为8 824.1 m3/s，其过程曲线如图25所示。

图25 瞬时全溃溃口流量过程Fig.25 Instantaneous total burst flow process

2)溃决洪水下游演进过程预测

溃决洪水下游演进过程一般采取对一维圣维南方程进行差分求解来进行模拟，即

(3)

式(3)中：A为过水断面面积，m2；Q为断面流量，m3/s；q为旁侧流量，m2/s；h为断面水深，m；t、x为距离和时间；R为水力半径，m；C为谢才系数。

选取堰塞坝下游重要乡镇及人类聚居区9个特征断面分析洪水演进及沿途淹没情况如下(图26、图27)：溃坝发生后河道下游出现明显的洪水波且随时间向下游演进，受干流和黑水河支流影响，洪水波波峰被削弱但不会随时间将至零点。溃坝后5.4 min洪水传至两河口，最大水位为1 662.88 m，洪峰流量3 567.57 m3/s，水位抬升9.17 m；溃坝后41 min洪水传至宁江堡，最大水位为1 612.39 m，洪峰流量2 086.78 m3/s，水位抬升6.59 m；溃坝后101 min洪水传至茂县县城，最大水位为1 565.09 m，洪峰流量1 791.88 m3/s，水位仅抬升2.09 m(表3)。

3.4.3 下游洪水淹没风险预测

下游洪水淹没风险分析是倾倒灾害链效应研究的重要步骤，对评估可能造成的社会、经济损失，制定人员避险方案，指导下游群众安全撤离等具有指导作用。

图26 各特征断面水位历时过程Fig.26 Water level duration process of each characteristic section

图27 各特征断面流量历时过程Fig.27 Flow duration process of each characteristic section

表3 溃决洪水下游演进过程特征参数预测Table 3 Prediction of characteristic parameters of downstream routing process of outburst flood

基于此，根据溃决洪水下游演进过程预测成果，提取河道各典型断面的最高水位，利用GIS(geographic information science)及遥感技术叠加数字高程模型，分析获取了洪水淹没范围(图28)。为指导地方政府有序撤离群众，按洪水位与国道和居民区的距离h，将下游区域划分为三个区域：危险区(h≤0 m)、警戒区(0

通过洪水淹没风险分析(图29)，溃坝洪水将对下游乡镇及G213线产生较大影响。洪水将会淹没金龙潭水电站、两河口、飞虹乡以及沟口镇，淹没及影响G213国道20.2 km，总危害资产4.5亿元。随着洪水向下游演进，影响逐步减弱，在茂县城区段，虽然存在较大流速带，但溃坝洪水水位淹没范围较小，不会对茂县城区造成影响。基于此划分石大关—沟口镇为危险区、沟口镇—渭门镇为警戒区、渭门镇—茂县县城警示区。