研究区位于地震活动相对较弱的长江中游地区。 该区域内的地质构造主要表现为NE向隔挡式褶皱及其伴生断裂, 断裂多为早、 中更新世活动断裂或前第四纪断裂, 主要包括齐曜山-金佛山断裂、 建始断裂和黔江断裂等。 其中, 黔江断裂是一条规模较大、 活动性较强的区域性断裂, 一般被认为是1856年大路坝地震的发震构造。 该地震的极震区与黔江断裂最近, 并位于该断裂的上盘。 黔江断裂的总体走向为NE, 倾向NW, 兼有右旋走滑性质, 主要发育于古生代地层中, 属切割较深的背斜扩展断裂(图1)。

图1

Fig. 1

	Figure Option ViewDownloadNew Window
	图1 崩滑体的位置与地层分布图(周鑫等, 2018)Fig. 1 Map showing locations of slump masses and strata(after ZHOU Xin et al., 2018).

崩塌滑坡区内主要出露的地层有志留系泥岩与砂岩、 二叠系灰岩, 第四系不发育。 崩滑体的原始岩体为志留系中统罗惹坪组紫红色砂岩与青灰色泥岩, 不等厚互层。 砂岩层厚6~10cm, 泥岩层理发育, 呈片状, 厚2~3cm, 泥岩夹层硬度低, 可污手。 岩层层面平直延伸, 产状为130° ~140° ∠5° ~8° , 向SE缓倾。 岩层中主要发育2组陡立的构造节理, 呈X形(图2), 2组节理的产状分别为290° ~330° ∠75° ~85° 和50° ~70° ∠80° ~85° (周鑫等, 2018)。

图2

Fig. 2

	Figure Option ViewDownloadNew Window
	图2 节理走向玫瑰花图(周鑫等, 2018)Fig. 2 Rose diagram of joint strikes(after ZHOU Xin et al., 2018).

1856年咸丰、 黔江地震诱发了一系列崩塌滑坡, 其中有4处崩塌规模较大, 分别位于重庆黔江区的小南海, 湖北咸丰大路坝乡的掌上界、 汪大海与活龙坪(图3)。 结合遥感影像可以清晰地看出4处崩塌滑坡整体沿NNE向分布, 与已知的黔江断裂或其分支断裂的延伸方向大体一致(图1)。 4处崩塌滑坡体的崩滑方向均为SE向, 与地层倾向一致, 为顺层滑坡。 其中, 小南海滑坡的崩滑规模最大, 崩滑体宽约1 600m, 滑移距离约为1 000m, 并堵塞河道, 形成堰塞湖; 掌上界次之, 崩滑体宽约400m, 长约600m; 汪大海位于两河交汇处, 两侧山体被沟谷切割较为陡峻, 受地形控制, 崩滑体由NW向SE逐渐变窄, 呈舌状展布, 滑移距离约为500m; 活龙坪滑坡为一块由奥陶系灰岩组成的大石板, 长约140m, 宽约100m, 倾向为125° , 倾角为15° , 为“ 山剥皮” 构造。

图3

Fig. 3

	Figure Option ViewDownloadNew Window
	图3 大路坝地震崩滑体的平面分布图Fig. 3 Distribution of landslides(yellow polygons)triggered by Daluba earthquake.

2.1.1 无人机摄影测量技术的应用

摄影测量技术是指运用摄影机和胶片组合测量目标物的形状、 大小和空间位置的技术, 目前已经进入数字摄影测量发展阶段。 数字摄影测量是传统摄影测量与计算机视觉相结合的产物, 可从数字影像中自动提取所摄对象的空间信息。 按平台距离可将数字摄影测量分为航天摄影测量、 航空摄影测量、 地面摄影测量、 近景摄影测量与显微摄影测量(高帅坡, 2017)。 无人机航空摄影测量技术常应用于滑坡研究中, 为了获取高精度的定性地貌数据, 其相对飞行高度一般在数十m至几百m, 属于摄影测量技术中近地航空地形测量的范畴。

本次工作中采用的DJI phantom3型无人机为多旋翼无人机, 自带GPS定位系统与4K镜头。 高清晰度的4K镜头极大提高了摄影测量数据的精度。 由高精度差分GPS完成数据转换是二维光学图像提取并分析三维地形数据的关键。

目前, 无人机航空摄影测量技术一般基于Structure from Motion(SfM)算法, 可在无序的数字图像中自动检测需匹配的特征点集, 并恢复、 优化摄像机之间的相互位置以及目标物的坐标, 获取目标物三维空间信息(Snavely et al., 2008), 最终利用Photoscan 软件生成DEM高程数据与DOM正射影像。

2.1.2 基于无人机的小南海崩滑体数据采集

在本次数据采集工作中, 对小南海崩塌滑坡体以及相邻地区的地貌形态数据进行采集。 手动操作DJI phantom3型无人机以蛇形前进的形式飞行和拍照, 确保每张航空图片正对地面且图片的重叠率> 50% 。 随后将所有图片导入Photoscan软件, 生成DEM高程数据与DOM正射影像。

无人机数据采集的理想区域为植被稀疏、 地貌面较为平坦的地区。 而本研究涉及的区域存在一些高差较大的地区, 可能会产生模型畸变或数据空洞, 且局部地区植被茂盛, 植物的枝叶在相邻航片里可能会表现出不同形态, 无法配对识别。 这些因素可能会使部分数据产生偏差, 对滑坡前的地貌形态恢复产生影响。 具体的数据校正方法将在下文中进行介绍。

由于利用无人机无法测量水域覆盖部分的湖底滑坡体的准确高程值, 只能采用传统的人工测量方法采集堰塞湖湖底的地形数据。

为得到水域覆盖部分的湖底滑坡体的准确高程值, 先利用GPS确定湖面的高程值为671m, 然后租用船只在湖面上进行人工测量, 具体的测量步骤为: 在湖中每隔5~10m进行GPS定点, 利用带有刻度的线绳悬挂重物扔至湖底测量水深, 并做好记录; 用湖面高程671m减去水深, 即得到各GPS点的湖底高程值(图4, 表1); 对所有带有高程值的GPS点进行插值, 可得到湖底的高程等值线图(图5)。

图4

Fig. 4

	Figure Option ViewDownloadNew Window
	图4 水域深度测点的GPS点位Fig. 4 Distribution of GPS sites to measure water depth.

表1

Table 1

表1(Table 1)

表1 湖底GPS点的高程值 Table 1 Elevation values of GPS sites on lakebed 序号 东经/(° ) 北纬/(° ) 水深/m 海拔/m 序号 东经/(° ) 北纬/(° ) 水深/m 海拔/m 1 108.746 76 29.654 35 11 660 27 108.741 76 29.649 80 42 629 2 108.747 49 29.654 33 6 665 28 108.742 61 29.649 75 36 635 3 108.746 31 29.653 86 10 661 29 108.743 23 29.649 70 28 643 4 108.747 22 29.653 82 25 646 30 108.743 98 29.649 71 23 648 5 108.747 93 29.653 85 15 656 31 108.744 82 29.649 76 14 657 6 108.745 40 29.653 10 18 653 32 108.745 69 29.650 45 7 664 7 108.746 26 29.653 08 22 649 33 108.746 64 29.650 44 19 652 8 108.746 96 29.653 09 26 645 34 108.747 46 29.650 47 13 658 9 108.747 80 29.653 13 21 650 35 108.748 23 29.650 45 10 661 10 108.748 45 29.653 17 14 657 36 108.748 88 29.650 47 8 663 11 108.744 93 29.652 59 21 650 37 108.745 69 29.649 55 18 653 12 108.745 85 29.652 57 30 641 38 108.746 56 29.649 58 26 645 13 108.746 61 29.652 56 32 639 39 108.747 42 29.649 57 22 649 14 108.747 36 29.652 55 25 646 40 108.748 17 29.649 57 20 651 15 108.748 09 29.652 52 2 669 41 108.742 50 29.648 83 43 628 16 108.745 10 29.651 88 10 661 42 108.743 66 29.648 77 39 632 17 108.746 02 29.651 91 25 646 43 108.744 47 29.648 80 32 639 18 108.746 84 29.651 90 12 659 44 108.745 41 29.648 76 35 636 19 108.747 92 29.651 89 4 667 45 108.746 31 29.648 79 31 640 20 108.745 82 29.651 14 7 664 46 108.747 52 29.648 74 18 653 21 108.746 76 29.651 15 26 645 47 108.745 39 29.648 12 42 629 22 108.747 68 29.651 21 22 649 48 108.746 44 29.648 11 43 628 23 108.742 55 29.651 08 18 653 49 108.747 67 29.648 14 17 654 24 108.742 02 29.650 54 19 652 50 108.746 60 29.647 57 41 630 25 108.743 24 29.650 92 37 634 51 108.747 78 29.647 51 12 659 26 108.743 30 29.650 49 40 631

表1 湖底GPS点的高程值 Table 1 Elevation values of GPS sites on lakebed

图5

Fig. 5

	Figure Option ViewDownloadNew Window
	图5 湖底地形推测图Fig. 5 Inferred terrain of lakebed.

利用无人机采集数据生成高分辨率DOM影像(图6)与25m间距的地形图(图7), 通过观察、 分析DOM影像、 地形图以及人工水域测深推测的湖底地形图, 结合野外观察的滑坡边界的GPS点位, 可清晰地划定滑坡范围(图6)。 大垮岩崩滑体主要堆积于堰塞湖的北岸, 宽约700m, 长约1 000m, 平面特征为舌形。 堰塞湖东岸的堆积体由小垮岩和另一组崩滑体组成, 这与以往研究认为的崩滑体由大垮岩与小垮岩2处崩滑体组成不同。 在野外调查中发现紧邻小垮岩的东侧还有2处规模相对较小的陡壁(图8), 与大垮岩、 小垮岩后缘的陡壁相似; 陡壁面走向相同, 崩滑堆积体的风化程度相当, 可以认为是同一次地震形成的崩滑体。 小垮岩是堰塞坝的主要构成部分, 宽约500m, 长约1 000m, 平面特征为扇形。 新发现的崩滑体位于小垮岩东侧, 宽约500m, 长约600m, 平面特征为舌形。 大垮岩、 小垮岩及新发现崩滑体的主滑线均沿NW-SE向展布。

图6

Fig. 6

	Figure Option ViewDownloadNew Window
	图6 圈定的滑坡范围Fig. 6 Delineated scope of the Xiaonanhai landslide.

图7

Fig. 7

	Figure Option ViewDownloadNew Window
	图7 小南海主要的崩滑体分布图Fig. 7 Distribution of landslides at Xiaonanhai.

图8

Fig. 8

	Figure Option ViewDownloadNew Window
	图8 滑坡后缘的陡崖(周鑫等, 2018)Fig. 8 Rear scarps of landslides(after ZHOU Xin et al., 2018).

通过观察已获取的DOM影像并结合现场实地调查, 认为研究区周边有2种山体形态可能与小南海滑坡前的山体原始形态较为相似。 第一种如图9a所示, 山体斜坡较为平缓, 坡度为40° ~50° ; 第二种如图9b所示, 山体高耸陡立, 斜坡坡度为70° ~80° 。

图9

Fig. 9

	Figure Option ViewDownloadNew Window
	图9 小南海崩滑体周围的山体形态Fig. 9 Photo showing shape of mountains around the Xiaonanhai landslide.

由于2种山体形态的恢复过程相同, 故以图9a坡度为45° 的山体为例进行介绍。 观察DOM影像(图10)可以看出, 山体两翼形态不一, 产状有所变化。 在GIS软件中根据山体两翼的形态特征, 将两翼划分为几个与山脊线呈一定交角的矩形区域。 在每个区域中, 以5m为间隔提取1个高程数据, 作为后期山脊曲面拟合的已知控制点。

图10

Fig. 10

	Figure Option ViewDownloadNew Window
	图10 山体两翼的分区图Fig. 10 Zoning map of the two mountain flanks.

同样, 运用GIS软件在没有滑坡崩塌体覆盖的堰塞湖上、 下游的山谷地区以5m为间隔提取各个点位的高程数据, 作为后期山谷曲面拟合的已知控制点。

至此, 获得包含经度x、 纬度y以及高程ξ 等信息的全部曲面拟合的已知控制点。

本研究的回归方程拟合流程设计如图11 所示。 首先, 设某已知点的高程值ξ 与其平面坐标x、 y之间的关系式为ξ =f(x, y)+ε , 式中, ξ 为高程值, f(x, y)为ξ 的趋势值, ε 为误差值。 对每个已知点列出上述方程, 当已知点的个数大于参数个数时, 由最小二乘法拟合出参数, 得到回归方程(胡菊菊等, 2012); 然后, 再将已知点的二维坐标数据(经度值x, 纬度值y)代入回归方程, 得到回归高程值ξ ₁。 若实际高程值ξ 与回归高程值ξ ₁呈现良好的线性关系(图12), 则认为回归方程可信; 若实际高程值ξ 与回归高程值ξ ₁的偏差较大, 则可能的原因是测量区的植被密集, 对无人机采集数据产生了影响。 此时, 则需设置置信区间, 对数据进行去噪处理, 再次进行多项式拟合, 重新得出回归方程, 直至回归方程可信。 通过这种回归方程拟合程序设计流程, 求出每个区域的回归方程。

图11

Fig. 11

	Figure Option ViewDownloadNew Window
	图11 拟合程序设计流程图Fig. 11 Flow chart of fitting program design.

图12

Fig. 12

	Figure Option ViewDownloadNew Window
	图12 原始高程值与回归高程值的对比图Fig. 12 Comparison diagram of original elevation value and regression elevation value.

在大垮岩、 小垮岩及新发现的陡崖前划定需要拟合的区域, 如图13 所示。 结合野外观测的滑坡分界线与规模, 将图10 中划定的山脊线与矩形区域按比例放大, 使之分别与大垮岩、 小垮岩及新发现的崩滑体相当。 在山脊各区域与山脊外的山谷区域中, 以5m为间隔提取各个需要拟合高程的控制点, 获取该点的二维坐标信息(经度值x、 纬度值y)。 将需要拟合区域的所有控制点的二维坐标代入回归方程拟合出高程值ξ , 即可得到恢复原始地形的所有控制点的数据。

图13

Fig. 13

	Figure Option ViewDownloadNew Window
	图13 拟合恢复区域Fig. 13 Fitting reconstruction area.

将上述拟合生成的所有控制点的三维数据在GIS软件中利用插值法生成小南海发生滑坡前的原始地貌面, 如图14 所示。 图14a为模型A: 山体斜坡较为平缓、 坡度为40° ~50° 的缓坡; 图14b为模型B: 山体高耸陡立、 坡度为70° ~80° 的陡坡。

图14

Fig. 14

	Figure Option ViewDownloadNew Window
	图14 小南海滑坡前原始形态的DEM数据Fig. 14 DEM data of original shape before Xiaonanhai landslide.

通过上述恢复过程, 得到小南海滑坡前的2种地貌原始形态。 那么, 哪种原始地形与现实更为接近? 2008年汶川地震、 2013年芦山地震等的震后调查表明, 有些局部场地的地形对地震作用下的次生灾害发育具有明显的影响或控制作用(张建毅等, 2012)。 因此, 本节将通过Geostudio软件进行数值模拟, 从山体原始形态边坡稳定性的角度定性分析哪种原始地形与现实更为近似, 并通过在堰塞坝钻孔资料中获取的原始谷底高程数据, 从另一个角度验证所恢复的原始形态高程数据的可信度。

3.5.1 运用Geostudio分析边坡的稳定性

引起滑坡的根本原因是土体内部某个滑动面上的剪应力达到了土的抗剪强度, 使稳定平衡遭到破坏(张克恭等, 2010), 而边坡的稳定性安全系数是衡量边坡稳定性的一项重要指标。 边坡稳定安全系数(Fs)定义为滑动面上土的抗剪强度(Mf)与实际产生的剪应力M之比, 即Fs=Mf/M。 一般情况下, 当边坡的安全系数> 1时, 边坡处于稳定状态; 当边坡的安全系数< 1时, 就会发生失稳破坏, 其数值越小, 发生滑坡的危险越大。

Geostudio是一套专业、 高效而且功能强大的适用于地质工程和地质环境模拟计算的有限元仿真软件, 被广泛地应用于滑坡稳定性分析中(Gasmo et al., 2000; Rahardjo et al., 2001)。 运用软件中的SLOPE/W、 QUAKE/W等模块, 可对边坡的静态稳定性、 地震动响应等方面进行分析。

根据野外调查结果与本研究推测的山体原始形态, 建立图15 所示的岩石力学概念模型A与B。 模型长约1 100m, 高约350m, 所需的具体岩体物理力学参数见表2(黄青松, 2014)。 下文将分别介绍具体的分析结果。

图15

Fig. 15

	Figure Option ViewDownloadNew Window
	图15 岩石力学概念模型及静态分析结果 a 模型A; b 模型BFig. 15 Conceptual model of rock mechanics and static analysis results.

表2

Table 2

表2(Table 2)

表2 岩石力学参数 Table 2 Rock mechanical parameters 物理力学参数 泥岩 砂岩 重度/kN· m-3 24.35 25.6 内摩擦角/(° ) 34 36 内聚力/kPa 300 500 泊松比 0.23 0.183 阻尼比 0.02 0.02 杨氏模量/kPa 350 000 450 000 抗拉强度/kPa 1 320 2 020

表2 岩石力学参数 Table 2 Rock mechanical parameters

(1)静态稳定性分析

Geostudio的SLOPE/W模块使用极限平衡理论, 以摩尔-库仑准则作为基础, 通过假想的简单破坏面, 将滑动面以上的滑体划分成若干垂直的土条, 建立静力平衡方程, 计算每个垂直土条的安全系数, 再从所有土条中找到最小安全系数(刘华强等, 2010)。 如图15 所示, 在无外力作用的静态下, 通过模拟计算, 模型A的边坡安全系数> 1, 边坡稳定; 模型B的安全系数< 1, 可能发生滑坡。

(2)地震动力响应分析

图16

Fig. 16

	Figure Option ViewDownloadNew Window
	图16 水平向加速度时程曲线Fig. 16 Horizontal acceleration time-history curve.

模拟中采用QUAKE/W与SLOPE/W模块耦合建模。 QUAKE/W模块采用等效线性模型进行地震响应分析, 计算地震动作用下边坡的动应力状态, 并将所得动应力代入SLOPE/W模块, 结合极限平衡理论求解安全系数(王立纬等, 2018)。 将上述地震动水平向加速度数据输入模型进行模拟分析。 同时, 边坡对地震波具有垂直和水平放大作用, 坡体加速度在垂直方向随着高程增加而增大, 到坡顶时急剧增大(徐光兴等, 2008)。 不同地形凸起坡度所对应的地震动峰值加速度(PGA)的放大倍数不同, 坡角为45° 的PGA放大倍数为1.39倍, 坡角> 60° 的PGA放大倍数为1.547倍(郝明辉等, 2019)。 分别将地震动的水平向峰值加速度按倍数放大为0.57g与0.64g输入模型A与模型B进行模拟分析。

分析结果如图17 所示, 图17a与17b代表在原始地震动力响应下的分析结果。 模型A的安全系数为0.929; 模型B的安全系数为0.789, 与模型A相比下降了15%; 且通过图中绿色区域代表的最优滑动面可以得出, 模型B的滑坡规模约为模型A的3倍。 图17c与17d代表在地形放大效应下地震动力响应的分析结果。 模型A的安全系数为0.828; 模型B的安全系数为0.783, 与模型A相比下降了5%; 且模型B的滑坡规模扩大为模型A的5倍。

图17

Fig. 17

	Figure Option ViewDownloadNew Window
	图17 地震动力响应分析Fig. 17 Seismic dynamic response analysis.

因此, 通过分析边坡静态稳定状态、 地震动力响应以及考虑地形放大效应地震动力响应等多情景状态下边坡的稳定性安全系数以及滑坡的规模, 可以认为模型B发生大规模滑坡的可能性更大, 可信度更高。

3.5.2 利用堰塞坝钻孔资料验证恢复地形高程数据的可信度

在本次调查中作者曾走访黔江区政府, 查阅了《小南海水库病害整治配套工程可行性研究报告》工程地质部分的钻孔资料^①(① 水利部四川水利水电勘测设计研究院, 1999, 小南海水库病害整治配套工程可行性研究报告。)。 如图18 所示, 工勘部门沿堰塞坝布设了11处钻孔, 沿B1、 B2 2条测线绘制了坝体地质剖面图2 幅(图19)。 从图19 看, 钻孔K1、 K2、 K3、 K5、 K6、 K7、 K9和K10均探测至谷底, 滑坡前谷底高程基本约为610m。 沿B1测线对比滑坡前山体地质剖面与勘探剖面(图20), 可以看出本段谷底的拟合高程与钻孔资料所揭示的谷底高程基本一致, 说明拟合数据可信。

图18

Fig. 18

	Figure Option ViewDownloadNew Window
	图18 钻孔分布图Fig. 18 Map showing locations of boreholes.

图19

Fig. 19

	Figure Option ViewDownloadNew Window
	图19 崩滑体钻孔剖面图Fig. 19 Geologic cross section of landslide and collapse accumulations.

图20

Fig. 20

	Figure Option ViewDownloadNew Window
	图20 根据地貌形态恢复推测的滑坡前地质剖面图Fig. 20 The geological section map of pre-landslide based on the recovery of landform.

综上所述, 小南海滑坡发生前山体的原始形态应如图21 所示: 山体高耸陡立, 属坡度为70° ~80° 的陡坡, 山谷向NEE缓倾。

图21

Fig. 21

	Figure Option ViewDownloadNew Window
	图21 小南海滑坡前山体原始形态示意图(a)和小南海附近陡坡地形(大路坝乡)(b)Fig. 21 Conceptual model of rock mechanics and analysis of results(a), steep terrain near the Xiaonanhai landslide(b).

如图22 所示, 山体滑坡后会形成滑坡凹地与滑坡堆积体。 假设滑坡后不发生物质流失(如遇河水冲刷将引起物质减少)的前提下, 滑坡凹地的物质总量等于滑坡堆积体的物质总量, 理论上这种等量关系反映了物质平衡原理。

图22

Fig. 22

	Figure Option ViewDownloadNew Window
	图22 滑坡纵向剖面简化示意图Fig. 22 Simplified sketch of landslide profile.

依据物质平衡原理, 可在ArcGIS中进行填方与挖方分析。 将滑坡发生后的高程数据与滑坡前的高程数据相减, 数据值为负的区域是挖方, 为滑坡凹地; 数据值为正的区域是填方, 为滑坡堆积物。 填方与挖方的计算公式均为

V=Δ S₁× Δ h₁+…+Δ S_n× Δ h_n=Δ S ∑Δ h_i (1)

式中, V(m³)为填方或挖方的体积, Δ S(m²)为离散的计算单元水平投影面积, Δ h(m)为离散的计算单元高程变化值, n为离散的计算单元个数。 挖方总量与填方总量基本相等, 即为滑坡体的总体积。

本工作中, 在2.3节圈定的滑坡范围内, 利用ArcGIS对滑坡前、 后的高程数据进行1m× 1m精度的采样, 每个离散的计算单元的水平投影面积为1m², 即Δ S=1m²。 然后将重采样的2种高程数据进行对比, 开展填方与挖方分析。 分析结果为: 挖方离散的计算单元高程变化值总和为4.3× 10⁷, 即 ∑i=1nΔh=4.3×107m, V=ΔS∑i=1nΔhi=4.3×107m3, 即挖方为 4.3×107m3; 填方离散的计算单元高程变化值总和为 4.6×107, 即 ∑i=1nΔh=4.6×107m, V=ΔS∑i=1nΔhi=4.6×107m3, 即填方为4.6× 10⁷m³。

从上文中可以看出, 挖方(滑坡凹地)与填方(滑坡堆积物)的计算结果存在一定的差距, 通过分析认为是数据采集、 原始地形准确度以及体积计算范围的圈定等方面带来的误差影响所致。 除此之外, 分析认为挖方量的数值作为体积计算结果较为可信, 具体原因包括:

(1)在现场考察中发现的崩塌岩块体积巨大, 且堆积杂乱无章, 岩块之间有很多空隙, 孔隙度的变化必然导致滑坡堆积体的体积大于实际崩滑体的体积;

(2)小南海崩塌滑坡属于历史滑坡, 距今已超过150a, 滑坡遗址已被改造为地质公园, 且在滑坡堆积体上已经开始种植农作物, 考虑到人为搬运等后期改造因素, 堆积体体积必然存在一定的误差;

(3)小南海为地震堰塞湖, 堆积体可能存在局部溃决流失与泥沙淤积的情况。

综上所述, 认为挖方(滑坡凹地)的计算结果更加可信, 即崩塌滑坡体的体积为4.3× 10⁷m³。