年廷凯, 王国栋, 郑德凤, 王大伟

1.大连理工大学海岸和近海工程国家重点实验室, 辽宁 大连 116024

2.辽宁师范大学自然地理与空间信息科学辽宁省重点实验室, 辽宁 大连 116029

3.中国科学院深海科学与工程研究所, 海南 三亚 572000

0 引言

随着能源需求的快速增加,海洋开发已进入深水区,特别是油气资源储量丰富的南海,已成为我国海底能源大规模开发的重要战略基地[1]。但南海地处欧亚板块、太平洋板块和印度洋板块的交会处,其地质构造复杂、地震活动性强,滑坡地质灾害活跃;滑坡一旦发生,将会危害水下基础设施(如海底管缆、深水基础、油气开采井等)的安全,阻碍深海能源的开发进程。触发海底滑坡的因素很多,且机制复杂,总体来说地震是导致海底斜坡失稳的主要因素[2-4],故受到了广泛关注。特别是在我国南海北部陆坡,更是地震易发多发区域,其西起海南岛东南海域,向东延伸至台湾岛西南海域,这一地带地形起伏较大,形成了局域分布的多个峡谷陡坡群。在较小地震作用下可能发生局部地段的斜坡失稳,当遭受较大的地震荷载时,则可能诱发大面积的海底滑坡,因此科学评估海底斜坡地震稳定性十分重要[1,4]。

目前海底斜坡地震稳定性分析方法可归纳为数值分析方法、概率分析方法、极限分析方法和极限平衡方法4类。这些方法均有各自特点及适用性,如:数值分析方法可视化程度高,无需假定滑动面位置,在处理复杂海底斜坡稳定性问题上具有优势[5],但其计算精度很大程度上依赖于岩土本构模型的合理性;概率分析方法虽然考虑了土体参数的不确定性,但其需要庞大的海底地质调查资料,通常难以获得足够的样本[6];极限分析方法往往采用上限定理,考虑外力做功的总功率等于速度间断面上的内能耗散率作为临界条件,采用对数螺线破坏模式更接近实际[7-8],但实践中应用较少;极限平衡方法是通过临界滑动面上的抗剪强度和剪应力之间的关系评价边坡的整体或局部稳定性,具有理论简单、计算高效的优点,易于在工程中使用[9-10]。综合考虑峡谷区陡坡群稳定性分析的需求,聚焦区域尺度下海底斜坡群稳定性,极限平衡方法仍是一种适合且高效的分析方法。

受海底原位探测技术、探测成本等多种因素制约,精确获取海底地形及土体参数等相关数据仍较为困难[11]。海底斜坡稳定性分析,通常选取某一典型剖面做二维简化分析,忽略其侧向边界引起的三维效应。然而通过对剪切应变、超孔隙水压力比和预测位移等数据的深入研究发现,海底斜坡地震稳定性的三维效应是显著的,与二维分析结果的差值可达15%[12]。此外,陆上边坡稳定性分析也体现了显著的三维效应[13-14]。因此开展三维稳定性分析变得更加重要。

以南海北部陆坡神狐海域一典型海底峡谷陡坡群为例,根据Scoops3D建模原理,结合GIS技术对研究区数字高程模型进行适配化处理;在此基础上,利用获得的地形数据和土层力学参数,结合Scoops3D构建海底陡坡群三维计算模型,进而开展基于极限平衡原理的海底陡坡群地震稳定性分析,得出各工况下斜坡稳定性安全系数及其分布图,以期为相关海域水下基础设施的选址、设计、施工及防灾减灾提供技术支撑。

1 南海北部陆坡典型峡谷陡坡群概况

1.1 南海北部陆坡概况

南海北部陆坡是一个呈条带状分布的大规模海底斜坡(图1),西起海南岛东南部海域,向东延伸至台湾岛西南部海域。这一区域地质构造复杂、地震易发多发,滑坡地质灾害活跃;特别是区域地形起伏较大,形成了局域分布的多个峡谷陡坡群,为海底滑坡的发生提供了孕灾环境[15-16]。

数据来源见脚注①。

1.2 区域地震地质构造特征

南海北部陆坡地质构造复杂,其东北部陆坡更是地震易发、多发地带[15-16],如图2所示。由图2可知,截至2018年,南海北部陆坡区已发生4.0～5.0级地震558次,5.0～6.0级地震155次,6.0～7.0级地震32次,7.0级以上地震6次[17]。研究[16]表明,该区发生的地震多为浅源地震,占比达94.5%。总体来说,南海北部陆坡地震活动性分布不均匀,表现出东部强于西部、北部强于南部的分布特征[17]。其中,震源在南海北部陆坡区东北部的分布较为密集,且密集区附近普遍发育海底滑坡,可见这两者之间存在密切关系[18-19]。

据文献[17]修编。

1.3 典型海底峡谷陡坡群

研究区位于南海北部陆坡神狐海域(图1红色框内),区域内海底滑塌体和断层-褶皱体系发育。选取一典型海底峡谷陡坡群,利用GIS处理DEM(数字高程模型)高程数据,建立可视化模型,其三维地形地貌如图3所示。

C1—C8为海底峡谷。

通过ArcGIS对地形图分析可知,研究区面积约为1 300 km2,水深为500～1 500 m,在此区域发育有一典型海底峡谷地貌,沟-脊相间排列,北西部地势较高,南东部地势较低,总体上地形起伏较大,极易发生滑坡地质灾害。

海底峡谷作为特殊的海底地貌单元,有着特殊的地质动力背景[20],是未来深海油气资源开发的重点,受到了国内外众多学者的广泛关注。根据Jobe等[21]对海洋探测资料和谷中岩心样本的研究,海底峡谷从形态学上可分为两类:Ⅰ型和Ⅱ型海底峡谷(图4)。Ⅰ型峡谷的谷道剖面多为“V”字形,支脉多、谷道蜿蜒曲折;Ⅱ型峡谷多为“U”字形,谷道顺直,谷壁平滑[21]。从图3可见,本区峡谷剖面大多为“U”字形剖面,支流较少、谷道较为顺直光滑,可推断研究区内峡谷为Ⅱ型海底峡谷系统。进一步分析图3可见,研究区内发育有C1—C8共8条海底峡谷,总体呈SW—NE向展布,地形地貌较为复杂,地形坡度主要集中在0°～10°之间,但海底峡谷两侧谷壁坡度基本在10°以上,局部区域大于15°,最大坡度甚至可以达到73°。将坡度数据通过GIS技术处理,将DEM高程数据中的坡度以云图的形式建立为可视化模型,详见图5的研究区坡度分布图。

图5 研究区坡度分布图

为了深入研究各海底峡谷的地形地貌与形态特征,在8条峡谷的上段、下段分别截取横断面(图6)开展分析,结果列于表1,其中以最陡的C1峡谷为例说明其特征。C1峡谷呈NNW—SSE走向,全长约24 km,上段水深约1 050 m,宽约3.0 km,下切深度约180 m,下段水深约1 350 m,宽约5.5 km,下切深度约350 m,谷底宽度由上段的800 m加宽至下段的1 600 m。上段两侧谷壁均较陡,最大坡度可达73°,为整个研究区坡度最大处,推测此处为海底基岩,沉积物层较薄,才形成如此陡峭的峡谷谷壁;下段谷壁坡度较缓,但最大处仍大于30°,集中在东侧谷壁,呈带状分布。峡谷上段西侧有2处相邻的大型冲沟(图3),冲沟头部较宽,北侧冲沟长为2.8 km,宽约650 m,南侧冲沟长2.4 km,宽约400 m。整体上,上段谷底地形剖面近似“V”形,而下段则转为典型“U”形峡谷,谷底坡度小于3°(表1),其他各峡谷有类似或局部不同特征,详见文献[18]。

表1 8条峡谷上、下段剖面基本信息

图6 8条峡谷上、下段剖面位置图

2 基于Scoops3D的三维斜坡稳定性分析方法

2.1 基于三维极限平衡法的Scoops3D简介

Scoops3D是基于FORTRAN语言编制的开源计算程序[22],可基于真实地形要素实现三维边坡稳定性分析。图7显示了DEM和一个潜在滑坡体滑动面的三维透视图,还显示了2个留作参照的三维柱体,每个柱从柱基到旋转轴的距离Ri,j不同,其中i,j表示DEM单元位置,δ表示滑动方位角。先将每个用于计算的三维柱体单元的水平截面尺寸设定为10 m×10 m,在每个栅格上方形成虚拟球心;然后对每个球心按顺序搜索不同半径的球面截断斜坡表面,并针对涉及的三维栅格开展极限平衡分析,得到最危险的潜在滑动面;再计算得到稳定性安全系数,并认定为滑动面包含所有三维栅格的稳定性安全系数;之后依次计算其他虚拟球心及各自最危险的潜在滑动面,每个栅格都因球心位置的不同而有若干个最小稳定性安全系数,取最小值作为最终的计算结果;接着对于边界外的区域,建立虚拟DEM网格,利用边界内数据插值;最后计算边界区域的稳定性安全系数。

ei,j. 柱体底端到旋转轴的水平投影距离;x轴也称为i轴,y轴也称为j轴,z轴也称为k轴。据文献[22]修编。

图8a中,滑动方位角方向δ不同于滑动面的真实倾斜方向(ε),这种差异在三维滑动表面中是常见的。这里要注意一点,滑动面近似为每个柱底部的平面。图8b中W、N和s都作用在滑动面上;若施加地震荷载keqW,其作用点位于柱的中心(图8b),作用方向沿滑移方位角方向(δ),并在极限平衡计算中使用柱体底端到旋转轴的水平投影距离еi, j来计算。

基于刚体极限平衡理论,Scoops3D模型可采用Bishop简化法等三维扩展方法来计算边坡的整体或局部稳定性。边坡稳定性安全系数Fs可通过滑动面上的抗剪强度τ与抗滑剪切力s的比值获得,其计算公式见式(1),其中抗剪强度τ可根据式(2)的摩尔-库伦准则得出。

Fs=τ/s;

(1)

τ=c′+(σn-u)tanφ′。

(2)

式中:c′为土的有效黏聚力;φ′为有效内摩擦角;σn为正应力;u为孔隙水压力。

考虑地震荷载的三维边坡稳定性安全系数Fs,其计算公式见式(3)。采用基于有效应力原理的Bishop法(有效应力及抗剪强度参数),认为土体骨架受力,无需考虑水深及坡面受到的海水压力[2],采用饱和重度计算地震力,利用三维极限平衡式(式(3))直接进行海底斜坡稳定性安全系数的计算:

(3)

式中:Ahi,j=Ai,jcosεi,j,其中Ai,j为潜在滑动面面积,εi,j为滑动面倾斜角度;ci,j和φi,j分别为三维柱体有效黏聚力和摩擦角;Wi,j为三维柱体失效部分的重力;ui,j为三维柱体在滑动面处的孔隙水压力;αi,j为滑坡体滑移方向在滑动面内的视倾角;keq为拟静力地震系数;为方便起见,mαi,j=cosεi,j+(sinαi,jtanφi,j)/Fs,所有带(i,j)角标的参数都随着三维柱体位置的改变而不断发生变化。至此,考虑了三维地形与地震荷载对区域斜坡稳定性的共同影响,完成了峡谷陡坡群地震稳定性评价。

2.2 基于Scoops3D的海底陡坡群三维计算模型

利用Scoops3D构建三维域模型步骤如下。首先,对于三维域地形高程数据,由于本文研究区土层分为3层,依据分层情况处理10 m×10 m分辨率的原始栅格地形数据,并转换为ASCII格式的DEM数据文件(这也就意味着每个用于计算的三维柱体单元的截面为10 m×10 m),依次输入3层DEM数据构建基本地形和土层分布(图9)。3层DEM数据的区别在于高程值的不同,由于该区域存在高低起伏众多陡坡,计算量较大,且通过以往工程地质调查认为地层力学参数相对均匀单一,故设置3层地质模型,包括Layer_1_DEM、Layer_2_DEM和Layer_3_DEM,其中,Layer_2_DEM相较于地表Layer_1_DEM总体减少2 m,而Layer_3_DEM相较于地表Layer_1_DEM总体减少10 m。其次,利用已收集到的土体强度参数定义地下3层材料,按照上文中所提出的水平拟静力地震系数Kh分别取0(自重工况)、0.05(地震烈度Ⅵ度)、0.10(Ⅶ度)、0.20(Ⅷ度)、0.40(Ⅸ度)输入这5种工况。最后,基于Bishop简化法,设定体积和球体半径等搜索方式,分别计算这5种工况下海底峡谷陡坡群的地震稳定性。

图9 三维模型示意图

2.3 海底土性参数及拟静力地震系数

研究区分为3层均质土层,第一层为0～2 m,第二层为2～10 m,第三层为10 m以上,每层土体强度参数随深度依次增大,对于土体分层和强度参数分布考虑如下。

1)目前各海域新近发生的海底滑坡,其深度普遍较浅,如:台湾岛西南部小琉球海底滑坡的平均深度约为11 m[23],黄河三角洲海底滑坡平均深度约为6 m[24]。由此可见,近代出现的海底滑坡与古、老滑坡在滑动深度上有着显著的区别,近代海底滑坡多以浅层滑坡为主,深度在10 m左右。故本文侧重于研究海底陡坡群发生浅层滑坡的可能性。

2)海洋工程地质调查[25]表明,研究区内地层比较单一,一定深度内土层的不排水抗剪强度可认为随深度近似呈线性增大,即海底土体强度参数随深度增大。

3)各土层有效抗剪强度指标的取值,参见图9内部标注,主要基于南海北部陆坡神狐海域原状土样的室内土工试验结果和相关文献[26-27],并兼顾土体强度随深度增大而线性增加的变化模式[25],未考虑地震力引发的土体强度弱化效应。

4)拟静力地震系数的选取,综合考虑了研究区可能遭受到的Ⅵ—Ⅸ度地震烈度,具体地震系数取值见2.2节。

3 计算结果与分析

3.1 基于稳定性系数的海底陡坡群安全性分区原则

Scoops3D计算结果为海底陡坡的稳定性安全系数Fs。根据水工规范[28],对于破坏后给社会、经济和环境带来重大影响的1级边坡,在正常运行条件下其稳定安全系数可取1.30～1.50;同时根据水电规范[29],特殊工况下边坡抗滑稳定安全系数可取2.50。参考上述安全系数标准,并考虑海底峡谷区陡坡群完全处于水下,其水动力环境和地形地质条件复杂,还具有区别于陆上滑坡的不可见性等诸多特点,特别是为了保障海洋基础设施的顺利建设和长期安全运营,故提出海底峡谷区陡坡群稳定性安全系数标准:Fs<1.50时,视为不稳定状态,将其划分为高风险等级;Fs为1.50～2.50时,视为基本稳定状态,将其划分为低风险等级;Fs≥2.50时,则视为完全稳定状态,将其划分为安全等级。

3.2 峡谷区陡坡群地震稳定性分析

考虑到研究区范围很大,一次性计算整个区域会导致计算效率过低,故将研究区的8个海底峡谷按各峡谷两侧海底山脊线进行划分(图10),分别进行稳定性评估,以此解决算力不足问题。限于篇幅,以下仅介绍C1峡谷区陡坡群的稳定性分析。

图10 研究区C1峡谷分区图

图11显示了5种拟静力地震工况下C1峡谷区陡坡群的稳定性评价结果,以安全系数分布图的形式展示,其背景为峡谷地形的高程分布图或水深分布图。其中:Fs<1.50的栅格单元以红色作为标识,所组成的区域为高风险区;Fs为1.50～2.50时以黄色栅格单元作为标识,所组成的区域为低风险区;Fs≥2.50的栅格单元未标记颜色,显示为背景色(即高程分布图底色),所组成的区域为安全区。总体来看,随着水平拟静力地震系数(Kh)的增加,风险区范围在不断扩大,当Kh=0.20(工况四)风险区面积较前3种工况略有扩张,当Kh=0.40(工况五)低风险区面积扩张最为显著。对比分析可见,低风险区总体分布在峡谷两侧的谷壁上,而高风险区主要集中在峡谷头部西侧谷壁波折处、南侧冲沟谷壁和上述两处中间相隔的西侧谷壁上,东侧谷壁高风险区也有出现但面积较小,不过在Kh=0.40时峡谷下段东侧谷壁出现了面积较大高风险区,呈带状南北向分布。表2列出了各工况下各级风险区面积占比及稳定性安全系数分布范围。

表2 研究区C1峡谷各工况下风险区面积占比及安全系数分布范围

a. Kh=0;b. Kh=0.05;c. Kh=0.10;d. Kh=0.20;e. Kh=0.40。

从表2可以看出:当Kh在0～0.20之间时,高、低风险区面积占比均较小,最大值分别为1.87%和5.88%(工况四);当Kh=0.40时(工况五),高风险区的面积占比增至3.78%,约为工况四的2倍,而低风险区面积占比陡升到35.18%,与工况四相比提高了近5倍。

工况一—工况四条件下,海底陡坡群全域的最高安全系数均大于10.00,这是因为研究区北部区域存在许多平缓斜坡,提升了全域最高安全系数值,依据计算结果可认为局部地区极安全和稳定;工况五条件下最高安全系数降为6.43,这说明强震作用可以大幅削弱海底斜坡稳定性。

各工况下海底陡坡群全域的最低安全系数均小于1.00,表明C1峡谷陡坡群存在局部斜坡失稳的风险。另外,工况一—工况四条件下沿坡向方向高风险区的长度为20～90 m,而工况五条件下沿坡向方向高风险区的长度为30～130 m,局部区域(峡谷中段东侧谷壁)甚至可以达到200 m,足见地震荷载对海底陡坡群稳定性影响之大。特别地,受当前海洋地质调查手段和成本的限制,加之研究区海底土层力学参数获取不足,因此本研究简化了土体强度参数及其分布,认为土体强度参数按图9模式均匀分布,得出了现状(自重)条件下全域最低安全系数<1.00的结果(局部地带处于失稳状态)。而实际上,南海北部陆坡区峡谷陡坡群这一局部地带,存在较大坡度且仍保持自稳,是因为其岩土体强度远高于计算选用的土体强度参数。

3.3 8个峡谷陡坡群稳定性结果对比分析

其他7个峡谷区(C2—C8)陡坡群的地震稳定性分析结果如图12、13所示。对比8个峡谷区可见,地震作用下西部峡谷区(图11、12)陡坡群的失稳范围明显高于东部区域(图13);大部分潜在滑坡集中在两侧谷壁靠近山峰峰顶的位置,此处地势起伏大、坡度较高;另外,也有相当一部分潜在滑坡位于靠近谷底的谷壁上,如C1和C2峡谷下段的东侧谷壁,其上有面积较广的陡坡区域(坡度介于20°～30°之间),呈带状分布。显然,稳定性安全系数<1.50的斜坡(群)构成了区域滑坡易发区。

a. C6峡谷区陡坡区;b. C7峡谷区陡坡区;c. C8峡谷区陡坡区。①Kh=0;②Kh=0.05;③Kh=0.10;④Kh=0.20;⑤Kh=0.40。

4 结论

1)基于GIS技术的海底真实地形调查表明,神狐海域典型海底峡谷剖面大多呈 “U”字形,支流较少、谷道较为顺直光滑,推断研究区峡谷为Ⅱ型海底峡谷系统;谷底地形总体上呈西高东低,且各峡谷西侧谷壁的平均坡度均小于东侧,这种地形的形成可能与洋流流向有关;峡谷陡坡群的坡度集中于0°～10°,局部区域大于15°,平均坡度为5.3°,最大坡度可达73°。

2)基于Scoops3D所提供的三维扩展极限平衡方法,对整个峡谷区三维陡坡群开展了地震稳定性分析。整体上,地震作用下西部峡谷陡坡群的失稳范围明显高于东部区域;大部分滑坡集中在两侧谷壁靠近山峰峰顶的位置,此处地势起伏大、坡度较高;另外,也有相当一部分滑坡位于靠近谷底的谷壁上,如C1和C2峡谷区南段的东侧谷壁,其上分布面积较广的陡坡区域(坡度介于20°～ 30°),呈带状分布。

3) 从三维陡坡群的稳定性安全系数分布图可以直观显示不同地震烈度下峡谷区陡坡群的整体或局域稳定性,其中稳定性安全系数<1.50的斜坡(群)构成了区域滑坡易发区。结合坐标位置,可以为海底基础设施的选址、设计和施工以及防灾减灾等提供参考依据,从而可降低海底滑坡所导致的水下设施破坏、重大经济损失和人员伤亡。

致谢:中国地质科学院地质力学研究所张帅博士、中铁第一勘察设计院集团有限公司付崔伟博士在软件技术方面提供了支持,博士研究生张浩、硕士研究生吴家伟在图文编辑等方面提供了帮助,在此一并致谢。