海底滑坡是近海大陆架常见的海洋地质灾害，通常指在重力和外部触发力作用下海底斜坡沉积物向下滑动的过程^[1-2]. 海底滑坡极具破坏力，可以破坏海洋工程结构、扰乱海洋生态系统以及诱发海啸，威胁沿岸居民的生命和财产安全^[3-5]. 1969年Camille飓风引发的密西西比河三角洲海底滑坡，造成高达1亿美元的经济损失^[6]；1998年巴布亚新几内亚由于地震和海底滑坡诱发的海啸，夺去了2 000人的生命^[7]；2006年吕宋海峡的海底光缆因海底滑坡而被破坏，使我国与中南亚地区的通信中断了12 h ^[8].

目前国内外海底斜坡稳定性评价的方法有：极限平衡法（limit equilibrium method，LEM）、强度折减法（shear strength reduction method，SSR）和剪切带扩展法（shear band propagation method，SBP）. LEM是边坡稳定性分析常用的方法，根据作用于岩土体中潜在破裂面上抗滑力与下滑力之比，求该块体的安全系数（factor of safety，FS），该方法理论简单、计算方便^[9-10]. SSR将土体的抗剪强度指标黏聚力c和内摩擦角φ′比上折减系数F，不断增加F的值使土体的抗剪强度减小，直到土体发生破坏，此时的F就是边坡的FS. SSR通常结合有限元单元法或有限差分法使用，该法考虑了土体的本构关系，可以显示边坡的变形破坏特征^[11-12].

LEM假设当土体所受的切应力τ大于土体峰值抗剪强度τ_p时，破坏沿整个滑动面同时发生. SSR假设土体所有部位的τ_p同时折减至残余抗剪强度τ_r直到边坡失稳. 然而，自然界中常见的大型海底滑坡具有厚度浅（10~100 m）、坡体长（10~100 km）、体积大（可达10⁴ km³）的特征^[13]，目前已知大型的海底滑坡（Storegga滑坡）涉及的沉积物破坏体积约10³~10⁴ km³，比最大体积的陆地滑坡大2个数量级^[14-15]. 对已发生的海底滑坡的反分析表明，在滑动范围如此之大的滑坡中，滑坡并不是瞬时整体发生的，潜在破裂面上土体的抗剪强度参数也不是同时劣化的. 因此，采用LEM和SSR分析规模巨大的海底滑坡有其内在的局限性.

相比于以上2种方法，SBP提供了将海底斜坡不稳定带（τ>τ_p）扩展到斜坡上其他区域（τ_r<τ<τ_p）的扩展准则^[16-17]. SBP扩展准则使海底斜坡剪切带上某个最危险的不稳定带先发生破坏，并通过剪切带的增长扩展达到整体破坏的情况具有可行性，为自然界中常见的大型海底滑坡提供了简单合理的解释.

目前，SBP在国内还未见报道，本文首先对SBP的原理及其在海底滑坡稳定性分析的应用做简明的阐述；再基于高斯函数的SBP计算公式，编写海底斜坡稳定性分析软件；以舟山六横岛典型海底斜坡为分析案例，开展基于SBP和LEM的斜坡稳定性对比分析研究，并开展SBP输入参数对剪切带扩展系数R的敏感度分析，研究成果可为SBP法应用于海底斜坡稳定性计算提供理论参考.

SBP中的剪切带扩展思想起源于1973年Palmer等^[16]发表的经典文献，认为在已发生初始破坏的边坡中，初始破坏区域末端土体的抗剪强度会随着相对位移的增加而降低，导致边坡的进一步破坏. 如图1所示为土体强度的应变软化示意图，在剪切带中的土体发生相对位移而出现强度降低. 图中，δ_r为从τ_p降低到τ_r时土体的相对位移，土体特征位移

(1) $\bar \delta=\frac{{\displaystyle\int {\left( {\tau - {\tau _{\rm{r}}}} \right){\rm{d}}\delta } }}{{{\tau _{\rm{p}}} - {\tau _{\rm{r}}}}}.$

图1（b）中深色区域为不稳定带，浅色区域为不稳定带破坏区域的2个末端，当上部土体与剪切带发生相对位移时，该部分的土体会发生强度弱化. y=f（x）为边坡几何轮廓拟合的函数表达式；L_i即不稳定带长度，L_q为剪切带发生扩展后的剪切带长度，即准稳定带长度；h为剪切带上部土层的厚度，目前所勘测的海底滑坡滑动距离长，滑动面积大，滑动深度浅，因此假设h $\ll $L_i. τ 随着δ增加，逐渐从τ_p降低到τ_r. 当滑体开始变形滑移时，不稳定带末端区域土体的抗剪强度会随着相对位移的增加而降低，剪切带进一步扩展增长成为可能. Palmer等^[16]通过断裂力学推导出剪切带在无限线性边坡中的扩展准则，给出边坡临界剪切带长度L_cr的表达式，并得出结论当L_i大于L_cr时，剪切带会发生扩展直至斜坡整体破坏.

Palmer等^[16]提出的SBP适用于超固结粘土，2005年Puzrin等^[17]基于能量平衡准则将该法推广到砂土和正常固结土中. SBP发展至今已有46年，剪切带扩展思想已被很多学者广泛研究使用. Kvalstad等^[18]使用该法分析了挪威Storegga滑坡的触发因素和滑动机制. Locat等^[19]分析了加拿大东部斯堪纳维亚敏感黏土的渐进性破坏. Zhang等^[20]研究了边坡破坏的动态扩展准则. Bernander等^[21]研究了长自然边坡的渐进下滑，并用有限差分法分析了斜坡破坏的影响因素. SBP将海底滑坡看作是渐进破坏的过程，是分析海底滑坡稳定性的又一重要方法，并逐渐成为近几年国外相关领域的研究热点.

由于海底勘探设备的限制，海底滑坡不稳定带长度不易测定，在一定程度上影响了SBP应用于海底斜坡稳定性的研究. Puzrin等^[22]将SBP推广到S型曲线海底斜坡体，推导了正常固结土的非线性海底斜坡的剪切带扩展准则，根据地形条件、土体物理力学和外力荷载等参数计算出海底斜坡的不稳定带长度，解决了不稳定带长度难以确定的问题. 但是Puzrin等^[22]提出的基于指数函数和双曲线的SBP，只能处理对称图形. Adams等^[23]对大陆架海底斜坡轮廓进行统计，超过80%的海底斜坡轮廓为曲线，在非线性海底斜坡中，曲线S型斜坡由高斯函数控制的超过50%. 因此，本文将高斯函数与SBP相结合，介绍SBP在海底斜坡稳定性分析中的基本原理及步骤. 高斯函数表达式为

(2) $y=f(x)=a\exp \left( { - {{\left( {\frac{{x - b}}{c}} \right)}^2}} \right).$

式中：a为函数曲线的最大y值，b为函数曲线水平顶部所对应的x值，c为函数拐点处的x值.

SBP分析海底斜坡稳定性的步骤如下. 1）建立海底斜坡计算模型：根据海底斜坡地质勘测资料，确定海底斜坡几何参数、土体物理力学参数和外力荷载参数，选择合适的函数y=f（x）拟合海底坡体轮廓线. 2）确定切应力比r：根据勘测得的几何参数、土体物理力学参数和外力荷载参数计算斜坡土体r，确定海底斜坡上不稳定带（r>1）、准稳定带（0<r≤1）和稳定带（r≤0）的分布. 3）确定剪切带长度：分别令r=0、1，求出斜坡上L_i、L_q. 4）计算L_cr：根据临界剪切带长度表达式确定L_cr，并判断L_i是否大于L_cr. 若L_i>L_cr则不稳定带发生扩展，最终海底斜坡破坏长度为L_q；若L_i≤L_cr，不稳定带仅在靠近图1（b）的浅色区域发生微小的扩展，不会引发斜坡的整体破坏，该区域在整个剪切带上可以忽略不计，在这种情况下可以视为不发生剪切带扩展，最终海底斜坡破坏长度为L_i. 5）计算斜坡FS：用下滑力与不同剪切区域的抗滑力之比计算不稳定带、准稳定带和稳定带的FS，得到海底斜坡的FS分布图.

Palmer等^[16]提出无限斜坡在自重条件下的切应力比概念，随后Puzrin等^{[22, 24-25]}给出在自重和外力荷载作用下切应力比的表达式. r是SBP判断剪切带分布的重要判据，与斜坡几何参数、外力荷载和土体物理力学参数有关.

(3) $r=\frac{{{\tau _{\rm{g}}} + {\tau _{\rm{h}}} - {\tau _{\rm{r}}}}}{{{\tau _{\rm{p}}} - {\tau _{\rm{r}}}}}{\rm=}\frac{{s\left( {{\tau _{\rm{g}}} + {\tau _{\rm{h}}}} \right)/{\tau _{\rm{p}}} - 1}}{{s - 1}},$

其中

(4) $s=\frac{{{\tau _{\rm{p}}}}}{{{\tau _{\rm{r}}}}},$

(5) ${\tau _{\rm{g}}}=\gamma 'h\sin \alpha .$

式中：τ_g为重力产生的下滑切应力，τ_h为外力荷载产生的下滑切应力，s为土的峰值抗剪强度与残余抗剪强度之比，γ′为土体的有效重度，α为斜坡的倾角. 根据Hance^[26]对366例海底滑坡触发因素的统计，地震是触发海底滑坡的主要因素，约占45%，选取地震荷载作为外部荷载，τ_h的表达式为

(6) ${\tau _{\rm{h}}}={k_{\rm{h}} }\gamma 'h\cos \alpha .$

式中：k_h为地震影响系数，与地震烈度和水平最大加速度a_max有关^{[22, 24, 27]}. τ_p的计算式为

(7) ${\tau _{\rm{p}}}={\delta _{\rm{d}} }\tan\; (\varphi '){\sigma '_{\rm{n}} }={\delta _{\rm{d}}}k\gamma '(1 - {r_{\rm{u}} })h\cos \alpha .$

式中：σ′_n为土体所受有效正应力； $\varphi ' $为土的内摩擦角，令k=tan $\varphi ' $，称k为不排水抗剪强度系数，对海底正常固结土一般取0.20~0.30；δ_d为地震折减系数，当地震烈度为6度和7度时可取0.6，8、9度时可取0.7^[27-28]；r_u为距海底斜坡表面h处的归一化超孔隙水压力，r_u=u_e/(γ′h)，其中u_e为超孔隙水压力，假设u_e随着h线性增加，此时r_u为常数.

将式（5）~（7）代入式（3）可得

(8) $r=\frac{{s(\tan \alpha + {k_{\rm{h}}}) - J}}{{J(s - 1)}}.$

其中J=δ_dk（1−r_u）是中间变量. 对于非线性斜坡地形，r随着斜坡角度α改变，tan α=f ′（x），将高斯函数代入式（8）得

(9) $r=\dfrac{{s\left( {a\exp \left( {\dfrac{{ - {{\left( {b - x} \right)}^2}}}{{{c^2}}}} \right)\left( {2b - 2x} \right) + {c^2}{k_{\rm{h}} }} \right) - {c^2}J}}{{{c^2}J(s - 1)}}.$

式（9）可计算海底斜坡任意位置处的r.

对于式（3），在海底斜坡上r>1的区域，τ_g＋τ_h>τ_p，此时下滑切应力大于τ_p，为不稳定带，图2红色区域；在0<r≤1的区域，τ_r<τ_g＋τ_h≤τ_p，下滑切应力小于τ_p而大于τ_r，随着土体相对位移的增加，剪切带末端区域土体的抗剪强度τ弱化，因此该区域弱化后的土体抗剪强度有可能小于下滑切应力τ_g＋τ_h，使得土层进一步下滑，剪切带扩展增长，也有可能大于下滑切应力τ_g＋τ_h，剪切带不发生扩展. 如图2所示为斜坡剪切带和切应力比的分布示意图. 图中α_i是不稳定带2个末端处的坡度，α_q是准稳定带2个末端处的坡度. 本文将剪切带扩展后的区域称为准稳定带，对于r≤0的区域，τ_g＋τ_h≤τ_r，下滑切应力小于τ_r，该区域为稳定带. 对于式（9），在土体物理力学参数和地震参数确定的情况下，r只与斜坡轮廓函数的斜率f′（x）有关，r是关于x的一元函数，并随着斜率f′（x）的增大而增大，海底斜坡切应力比r分布如图2所示. 分别令r=1，r=0，求出不稳定带和准稳定带两端的坐标式，确定L_i和L_q.

(10) $\tan \left( {{\alpha _{\rm{i}}}} \right)=f'({x_{\rm{i}}} - {L_{\rm{i}}})=f'({x_{\rm{i}}})=J - {k_{\rm{h}}}.$

(11) $\tan \left( {{\alpha _{\rm{q}}}} \right)=f'({x_{\rm{q}}} - {L_{\rm{q}}})=f'({x_{\rm{q}}})=\frac{J}{s} - {k_{\rm{h}}}.$

Puzrin等^[22]在Palmer等^[16]的基础上，基于能量平衡准则，假设土体为弹塑性，进一步推导出在海底正常固结土非线性斜坡上剪切带的扩展准则，临界剪切带长度^{[22, 24-25]}

(12) ${L_{{\rm{cr}} }}=\frac{1}{{\bar r}}\left( {1 + \sqrt {\frac{{{E_1}}}{{{E_{\rm{u}}}}}} } \right)\sqrt {\frac{{s\cos {\alpha _{\rm{i}}}}}{{(s - 1)J}}\frac{{2{E_{\rm{u}} }\bar \delta }}{{\gamma '}}} .$

式中：E₁和E_u分别为土体在加荷和卸荷下的弹性模量， $\bar r$为不稳定带的平均切应力比，

(13) $\bar r=\frac{{s\left( {(f({x_{\rm{2}}}) - f({x_{\rm{1}}}))/({x_2} - {x_1}) + {k_{\rm{h}}}} \right) - J}}{{(s - 1)J}}.$

L_i与L_cr的大小关系是SBP判断剪切带是否扩展的关键，为了表达方便，将两者之比定义为剪切带扩展系数R=L_i/L_cr，当R>1时，剪切带发生扩展.

1）不稳定带，斜坡发生破坏，τ_p只能短暂地保持不变，随后会被外部荷载力弱化为τ_r. 因此不稳定带以τ_r计算安全系数

(14) ${{\rm{FS}} _{\rm{i}}}(x)=\dfrac{{{c^2}J}}{{s\left( {a\exp \left( {\dfrac{{ - {{\left( {b - x} \right)}^2}}}{{{c^2}}}} \right)\left( {2b - 2x} \right) + {c^2}{k_{\rm{h}}}} \right)}}.$

2）稳定带，由重力和地震荷载产生的切应力已经小于土体的残余抗剪强度（τ_g＋τ_h≤τ_r），该区域土体稳定，并不受到扰动而发生强度弱化，稳定带安全系数应以τ_p计算：

(15) ${{\rm{FS}} _{\rm{s}}}(x)=\dfrac{{{c^2}J}}{{a\exp \left( {\dfrac{{ - {{\left( {b - x} \right)}^2}}}{{{c^2}}}} \right)\left( {2b - 2x} \right) + {c^2}{k_{\rm{h}}}}}.$

3）准稳定带，确定安全系数的表达式分2种情况：当不稳定带发生扩展时（R>1），准稳定带土体抗剪强度降低到残余抗剪强度，因此准稳定带安全系数FS_q1（x）=FS_i（x）；当不稳定带不发生扩展时（R≤1），剪切带扩展区域的土体抗剪强度不发生弱化，FS_q2（x）=FS_s（x）.

在极限平衡法中，一旦τ_g+τ_h>τ_p，海底斜坡将沿剪切带同时发生破坏. 因此以τ_p计算海底斜坡整体的安全系数FS_LE（x）=FS_s（x）.当使用SBP计算安全系数时，首先利用r对斜坡进行分区，得到L_i再与L_cr对比以判断剪切带是否发生扩展，若发生扩展则在准稳定带采用τ_r计算安全系数，若不发生扩展则采用τ_p计算安全系数. 极限平衡法不对斜坡进行分区，认为斜坡发生破坏时沿整个剪切带滑移，以τ_p计算安全系数.

Puzrin制作了一个试验设备来验证海底斜坡的剪切带扩展^{[24, 29]}. 对于剪切带的扩展长度，实验分析结果显示，SBP分析结果和数值分析结果的匹配误差在10%以内，表明它们都合理地拟合了测量的剪切带长度. 该试验提供了斜坡SBP的直接试验证据，也验证了SBP能够对剪切带扩展的长度进行量化.

以舟山六横岛一处典型的海底斜坡为案例^[30]，分别采用基于高斯函数的SBP和LEM开展海底斜坡稳定性分析研究. 六横岛是舟山群岛的第三大岛，位于浙江省舟山群岛南部，面积98.0 km². 六横岛深水岸线资源丰富，全岛10 m水深以上且有一定腹地配套的岸线36.3 km，占舟山市可用深水岸线的21.9%. 依托优越的深水岸线资源，六横岛已兴建较多深水码头泊位. 选取一个该岛典型的海底斜坡，如图3所示，土体物理力学参数及地震参数取值，如表1所示.

舟山市的抗震设防烈度为7度，本文设计地震烈度为7度，此时k_h取0.08，δ_d可取0.6^[27-28]；六横岛土质天然含水率为35.50%~62.80%^[30]，内摩擦角4.50°~15.90°，k取0.25；舟山海域土体灵敏度一般为4~8^[31]，假设s取5.0；E₁=300 τ_p，E₁ / E_u=0.5. 由于地质资料较少，且缺乏实验数据，并未考虑弹性模量在土体强度衰减过程中土体弹性模量的变化. 基于地震参数、土体物理力学参数、斜坡几何参数以及拟合的高斯曲线函数y=f（x），结合前文的SBP分析流程，计算该海底斜坡的剪切带分布和安全系数分布情况.

1）建立海底斜坡计算模型：采用式（2）拟合海底斜坡轮廓，可得高斯表达式系数a=77.95，b=1 011，c=707.5，海底斜坡地形及高斯函数拟合曲线如图4所示. 2）将土体物理力学参数和地震荷载参数代入式（9），计算斜坡任意处的r，r随着x的分布曲线如图5所示. 3）计算L_i、L_q：根据计算出的r，令r=1，由式（10）可以求出不稳定带两端的坐标x_i=754.74 m，x_i − L_i=215.64 m，得出L_i=539.10 m. 对于准稳定带，根据r的分布情况可知，整个斜坡上r>0，因此当不稳定带发生扩展时，案例中整个海底斜坡都会发生破坏，L_q=1 000 m. 随后需要计算R以判断剪切带是否发生扩展. 4）计算临界剪切带长度L_cr：由式（12）、（13）求出L_cr=130.33 m，R = 4.14，不稳定带发生扩展，故案例中整个海底斜坡最终都会发生破坏，破坏长度L_q=1 000 m. 5）计算FS. 根据式（14）、（15）分别计算该斜坡的安全系数，计算结果如图6所示. 同时，采用LEM计算该斜坡的安全系数，令FS_LE（x）=1可以求得不稳定区域两端的坐标x_le、x_le−L_le，计算可得x_le=754.74 m，x_le−L_le=215.64 m. 因此，LEM计算下得到破坏长度L_le=L_i=539.10 m.

基于以上分析，在斜坡剪切带破坏长度方面，当R≤1时，不稳定带不发生扩展，LEM和SBP得出的结果都表明斜坡只在不稳定带区域发生破坏，最终破坏长度一致为L_i. 当R>1时，SBP计算结果表明斜坡除了在不稳定带区域发生破坏之外，剪切带会继续扩展，最终破坏长度为L_q，SBP计算的最终破坏区域比LEM计算的大得多. 本例LEM计算的破坏带长度为539.10 m，采用SBP计算时斜坡整体都会发生破坏，破坏长度为1 000 m. 在安全系数方面，不稳定带和准稳定带中SBP比LEM计算得到的安全系数低，这是由于SBP在破坏的区域以残余抗剪强度计算安全系数，而LEM采用峰值剪切强度计算安全系数.

用Matlab APP Designer开发基于SBP的海底斜坡稳定性分析程序，如图7所示. 将该程序与Puzrin等^[22]案例进行对比分析，结果吻合.

该程序具备如下功能：1）输入海底斜坡地形数据，显示海底斜坡几何地形，选择相应的拟合函数，计算拟合参数；2）输入海底斜坡的土体物理力学参数和地震参数，可计算海底斜坡的不稳定带长度、准稳定带长度、临界剪切带长度以及剪切带扩展系数；3）得到海底斜坡切应力比分布曲线，安全系数分布曲线以及海底斜坡剪切带分布图；4）可开展敏感度分析，研究土体物理力学参数和地震参数对剪切带扩展系数的影响，选取需要分析的敏感性参数，输入最小值、最大值和变化量可得到敏感度分析图.

SBP分析海底斜坡稳定性的过程中涉及很多输入参数，有必要开展参数敏感度分析研究，分析参数对剪切带扩展的影响程度. 以表1的参数取值作为敏感度分析的基准值，敏感度分析的参数取值范围如表2所示. 由前文可知，判断剪切带是否发生扩展是SBP最重要的步骤，因此本文开展SBP输入参数对R的敏感度研究. 使用海底斜坡稳定性分析程序分析土体物理力学参数k、s、γ′、 $\overline \delta $，地震参数k_h、δ_d对R的影响程度，结果如图8所示。由图可以看出，随着输入参数取值的增大，R有2种趋势：1）R随着参数取值的增大而增大。图8（a）、（b）、（c）对应的s、k_h、γ′. s越大，准稳定带土体的抗剪强度受到扰动后降低得越严重，不稳定带更易扩展到准稳定带；k_h越大代表地震荷载越大，地震荷载产生的下滑力τ_h增大，引起L_i增大，使R增大；γ′增大引起土体自重产生的下滑力τ_g增大，增大L_i且减小海底斜坡的L_cr，使R增大.2）R随着参数取值的增大而减小。图8（d）、（e）、（f）对应的k、δ_d、 $ \bar{\delta }$. 随着k增大，土体抗剪强度增大，使得L_i减小，R减小；δ_d与k作用相同，随着δ_d增大，土体抗剪强度增大，不稳定带不易发生扩展；随着 $ \bar{\delta }$的增大，土体抗剪强度弱化到τ_r所需要的位移量增大，使得土体强度不易发生弱化，提高了准稳定带的抗剪强度，使R减小.

采用数理统计的方法计算敏感度系数S，量化各个输入参数对R的影响程度. 敏感度系数的定义如下^[32].

(16) $\begin{split} \\ S=\frac{{{{\left| {\Delta F} \right|} / {F{_0}}}}}{{{{\left| {{\Delta}X} \right|} / {\left( {{X_{\max }} - {X_{\min }}} \right)}}}}. \end{split}$

式中： ${ {{\Delta}X} } $所示为参数X的变化量，X_max、X_min分别为参数X的最大值和最小值，F₀为各参数的基准值所计算出的R基准值， ${\Delta F}$为某一参数变化所引起的R变化量. 将每个步长过程中的敏感度系数相加求出平均敏感度系数

(17) $\bar S=\frac{{{S_1} + \cdots + {S_{\rm{N}} }}}{N}.$

由（17）分别计算SBP输入参数的平均敏感度系数，结果如表3所示. 各参数敏感度系数由大到小排序如下：为k_h>δ_d>k> $\bar \delta $>γ′>s. k_h、δ_d、 $\bar \delta $对R的影响最大（ $\overline S \geqslant 0.85$），远高于其他输入参数（ $\overline S \leqslant 0.48$），输入参数s对R的影响最小（ $\overline S=0.09$）. 因此，在应用SBP分析海底斜坡稳定性时，应尤其注意k_h、δ_d、 $\bar \delta $的选取.

（1）介绍SBP的原理及其在海底斜坡稳定性分析中的应用. 相比于LEM，SBP将海底滑坡看作是渐进破坏的过程，允许剪切带发生扩展，能更加合理地解释自然界中大型海底斜坡失稳机理.

（2）编写基于高斯函数的SBP海底斜坡稳定性分析程序，以舟山六横岛典型海底斜坡为例，开展基于SBP和LEM的边坡稳定性对比分析. 研究发现，在斜坡剪切带破坏长度方面，当R≤1时，LEM和SBP得到的最终破坏长度一致. 当R>1时，SBP计算得到的最终破坏区域比LEM计算得到的破坏区域大. 在安全系数方面，SBP比LEM计算得到的海底斜坡破坏区域安全系数低.

（3）开展SBP输入参数对R的敏感研究，计算敏感度系数，量化输入参数对海底滑坡剪切带扩展的影响程度，研究发现k_h、δ_d、 $\bar \delta $对R的影响最大，s对R的影响最小.