李存兴，莫建新，王 成

(中交第二航务工程勘察设计院有限公司，湖北 武汉 430060)

近几年来，伴随着国家“一带一路”倡议的拓展，我国工程建设企业在国际工程领域的影响力逐步扩大，进而参与建设的项目也逐渐增多。然而，不少工程位于环太平洋地震带或欧亚地震带，地震烈度高，水工结构抗震设计成为工程设计的重中之重。

斜坡式护岸作为围海造地项目常见的结构形式，其合理的抗震设计将有利于节约工程总投资。国内规范对斜坡式护岸地震作用下稳定性分析采用的方法是拟静力法，但拟静力法只提供一个安全系数，不能反映边坡最终变形情况。当安全系数小于1.0时，国际通用规范均允许采用位移法评价护岸工程在地震作用下的安全性。美国土木工程师协会和海岸、海洋、港口及内河研究所联合发布的码头抗震设计规范ASCE/COPRI61-14[1]，国际航运协会的抗震设计规范PIANC WG34[2]，英国建筑工业研究与情报协会、美国陆军工程兵团和法国生态可持续发展和能源部等联合发布的国际堤防工程手册USACE&CIRIA C731[3]等都明确规定：当安全系数小于规定时，可以通过计算堤体位移来评估结构的安全性，比传统拟静力分析更能反映结构的破坏程度。

本文结合工程实例，介绍了几种边坡地震位移简化计算方法，分别选取临界地震系数、峰值地震系数、峰值速度、地震震级、结构周期、地震反应谱或阿里亚斯烈度等指标获取最终的边坡永久位移，并基于美国地质地震局SLAMMER程序的SLIDE2018软件位移分析进行对比。

1 地震位移简化计算

1.1 理论简介

Newmark滑块位移法[4]最初由Newmark于1965年提出，用于岸坡和土石坝的永久变形计算。半个世纪以来，经过众多研究者不断探索和实践，这种计算方法和模型在理论和工程应用等方面都得到极大的发展。如图1所示的刚塑性滑块模型，滑块代表具有任意形状滑动面的潜在滑动量，a(t)为地震波加速度时程。最初的Newmark滑块位移法假定，坡体的地震荷载可以用一个拟静力荷载代替，拟静力荷载为拟静力地震系数k与潜在滑动体重力的乘积，地震中岸坡破坏会形成明显的滑动面，当地震荷载超过滑动体的极限承载力时，滑动体会沿着潜在的滑动面产生塑性位移。

假定屈服加速度(临界加速度)在整个地震历时中保持恒定(地震过程中土体强度不会发生明显退化)。当施加在滑块上的地震加速度超过屈服加速度(岸坡拟静力稳定安全系数为1时的加速度)时，即t1时刻滑块位移发生，滑块的速度可以通过对阴影部分积分得到，滑块的速度会一直增加到t2时刻，此时加速度再次降到屈服加速度以下，并且随着加速度反向，速度最终在t3时刻减为零，滑块位移可以通过对速度时间的关系进行积分得到。由图1可知，滑块位移的大小取决于施加加速度的幅值(地震动振幅)、单次历时时间(频率)和在整个地震历时中加速度超过屈服加速度的次数(持时)。

图1 Newmark算法原理

1.2 地震位移简化计算公式

近年来，国外诸多学者和行业协会通过对地震数据进行收集、分析和研究，基于Newmark滑块位移概念，总结得出了多种典型的地震下永久位移的简化经验计算公式。经过归纳分析，根据考虑要素的不同，可以分为6大类。

1.2.1考虑临界加速度和峰值加速度

Ambraseys & Menu[5]分析11次地震的50个强震记录，提出了各种再回归方程估算Newmark位移作为临界加速度的函数，并根据研究结果提出了位移计算公式：

(1)

式中:DN为地震下永久位移；kmax为峰值地震系数，即地震峰值加速度amax与重力加速度g的比值；kc为临界地震系数，即临界加速度ac与重力加速度g的比值，可以通过试算确定或软件工具求得。

1.2.2考虑临界加速度、峰值加速度和地震震级

Jibson[6]对1994年美国加州Northridge市地震诱发的滑坡事故建模计算进行危险性评估，考虑到地震震级的影响，得出了一个适用于震级在[5.3, 7.6]范围的位移计算公式：

0.424Mw±0.454

(2)

式中：Mw为地震震级。

1.2.3考虑峰值速度、峰值加速度和临界加速度

美国公路合作研究组织的《挡土墙，埋入式结构，斜坡和堤防的抗震分析与设计手册》NCHRP 611[7]对地震数据进行相关性分析，给出了永久位移与峰值速度、临界地震系数、峰值地震系数的关系，常用的计算公式为：

lgDN=-1.51-0.74lg(kC/kmax)+3.27lg(1-kC/kmax)-

0.80lgkmax+1.59lg(PGV)

(3)

式中：PGV为地震峰值速度(以英寸计)。

美国联邦公路局国家公路研究所出版的《交通工程岩土特性与结构地基抗震设计参考手册》(FHWA-NHI-11-032)[8]也采用上述公式。

1.2.4考虑滑动体基本周期、地震反应谱、临界加速度和地震震级

Bray & Travasarou[9]通过收集大量的地震资料，运用非线性耦合黏滑块模型，采用概率方法分析地质稳定性，得出了不同条件下非零位移的计算公式：

lnDN=-1.10-2.83lnkc-0.333(lnkc)2+0.566lnkc·

ln[Sa(1.5Ts)]+3.04ln[Sa(1.5Ts)]-0.244·

ln[Sa(1.5Ts)]2+1.5Ts+0.278(Mw-7)±ε

(4)

式中：Ts为滑动体基本周期，按照图2模式计算；Sa(1.5Ts)为1.5倍基本周期Ts对应的地震加速度反应谱；ε一般取0，正负偏差0.66。

注：对于滑坡体高度大于边坡高度的情况，Ts=4H/vs，其中H为滑动体高度，vs为滑坡区域土体的剪切波速。

图2潜在滑动体的基本周期计算模式

1.2.5考虑滑坡体高度、坡顶峰值加速度、临界加速度和地震震级

Makdisi & Seed[10]在Newmark模型的基础上进行了改进，采用等效线性模型计算潜在的滑动体平均地震反应。利用静力分析得到相关计算参数，从而近似评估岸坡在地震作用下的永久位移。基本步骤归纳如下：

1)根据地震峰值基底加速度PBA(peak base acceleration)查图3a)求得坡顶峰值加速度PCA(peak crest acceleration),为方便对比计算，根据曲线拟合了两者关系:

(5)

2)根据滑弧高度与堤身高度的比值，查图3b)求得滑块的平均加速度。由于滑弧高度大于堤身高度，按照y/h=1.0对应的上限求得滑块的平均加速度u0max= 0.47PCA。滑动体的平均最大地震系数k0max=u0max/g。

3)求得滑块滑动临界地震系数kc。

4)计算临界地震系数与平均最大地震系数的比值：

λ=kc/k0max

(6)

5)计算滑坡体结构周期Ts。

6)查图4求得滑体位移，图4a)为位移范围，4b)为位移中值。为方便对比计算，对图4b)数据进行了曲线拟合分析。

图4 不同地震震级下永久位移与临界地震系数的对应关系

3条主要曲线结果如下：

①地震震级为7.5级时，位移采用公式(7)计算：

(7)

式中：D75为7.5级地震下的位移。

②地震震级为6.5级时，位移采用公式(8)计算：

(8)

式中：D65为6.5级地震下的位移。

③地震震级为8.25级时，位移采用公式(9)计算：

(9)

式中：D825为8.25级地震下的位移。

1.2.6考虑阿里亚斯烈度、临界加速度、峰值加速度

1970年，智利工程师Arturo Arias提出了一种地震动强度的参数[11]，定义为地面加速度的平方在时间上的积分，包含了地震动振幅、频率以及持时的地震动参数，在表征地震诱发滑坡等地震灾害危险性方面具有重要的意义，其表达式为：

(10)

式中：Ia为阿里亚斯烈度；g为重力加速度；a(t)为加速度(m/s2)；t为时间；Td为有效震动持续时间。一般来讲，地震的加速度越大、频率越低、持续时间越长，阿里亚斯烈度越大。

Jibson[12]通过进一步研究，提出了考虑阿里亚斯烈度、临界加速度及峰值加速度的位移计算公式：

lgDN=0.561 lgIa-3.833 lg(kc/kmax)-1.474±0.616

(11)

Shang-Yu Hsieh和Chyi-Tyi Lee[13]通过收集、分析和研究世界范围内的地震加速度数据，给出了岩基(式(12))和土基(式(13))的最大永久位移计算公式：

lgDn=0.788 lgIa-10.166kc+5.95kclgIa+1.779±0.294

(12)

lgDn=0.802 lgIa-10.981kc+7.377kclgIa+1.914±0.274

(13)

式中：Ia为阿里亚斯烈度。

1.3 地震位移数值模拟

随着计算机科学的快速发展，利用软件进行数值模拟计算地震下的永久位移成为分析地震位移最有效的工具。SLIDE 2018是由加拿大RocScience公司研发的一款功能强大的二维边坡稳定性的分析软件，其在条分法的基础上根据极限平衡进行边坡的稳定性分析，适用于各种类型的土质和岩质；其地震位移分析的核心计算基于美国地质地震局的SLAMMER程序，是目前较为权威的算法；利用该软件来分析地震作用下边坡的Newmark位移分析，可以同时考虑刚性、耦合和解耦算法，为地震下边坡风险分析和工程设计提供客观依据。

2 工程案例

2.1 工程概况

菲律宾某滨海新城项目，拟吹填造陆形成318万m2地块。工程水工建筑物主要为围堰及护岸结构，永久工程护岸总长度约为10.2 km；护岸结构形式根据使用功能要求和水文地质等因素，不同区段采用斜坡结构或半斜坡半直立结构。

拟建人工岛护岸区域水底泥面高程为-13.35～-4.86 m；地貌类型为海岛海岸地貌；工程区地质条件复杂，表层软土具有强度低且压缩性高的特性，分布不均匀，工程场地土自上而下划分为多个土层，其中②-1层含砂高液限黏土、②-2层含砂黏土、③-1层高液限黏土层物理力学指标较差，须进行地基处理；深层土的多个亚土层为含砂高液限黏土、高液限黏土、含砾黏土质砂、胶结砂、砂质弹性粉土等[14]，埋深较深且力学性能指标相对较高，不进行地基处理。

堤身主体采用常规的抛石堤结构，堤顶设挡浪墙和景观平台，外侧采用人工块体护面，内侧设倒滤设施，斜坡段结构的典型断面见图5。地基处理主要采用挤密砂桩，根据水深和地质条件等不同，处理宽度为40～75 m，处理深度北岛平均为14 m、最大为27 m，南岛平均为8 m、最大为13 m；软土层厚度较薄的区域采用开挖换填中粗砂振冲处理、局部采用抛石挤淤。

根据本项目地震灾害专项评估报告[15]，对于项目场地，level 1(72 a一遇)地震加速度为0.28g(地表加速度，已经考虑场地放大,下同)，level 2(475 a一遇)地震加速度为0.44g。

注：N为标贯击数。

2.2 地震波选取

边坡的永久位移(地震的破坏性)不仅与峰值加速度PGA有关，同时与地震持续的有效时间、频率以及地震的阿里亚斯烈度密切等相关。地震波波形不仅仅与地震震级相关，而且与震源深度、震中距离、断裂带的距离、场地特征等密切相关。整体稳定分析时，为了尽可能接近场地特征，从数据库中选取实测地震波，并根据场地特征进行谱匹配。根据地震灾害评估报告的建议，选取8组地震32条地震波(其中16条为72 a重现期，16条为475 a重现期)经过目标谱匹配后生成的人工地震波可用于设计。其中475 a一遇地震波的基本参数见表1，不同波对应的加速度反应谱都接近目标谱，见图6，不同波的时程曲线见图7。

表1 抗震分析计算选用地震波

注：持续加速度SMA为地震波第3大峰值加速度。

图6 场地475 a一遇地震反应谱

图7 相同反应谱下不同地震波对应加速度时程曲线

不同地震波虽然差异巨大，峰值加速度、持续加速度、峰值速度、持续时间、阿里亚斯烈度都不一样，例如1b波峰值加速度为0.584g、8a为0.364g，但对应的反应谱是一致的(反应谱的基本加速度都是0.44g)，都是可以适用于场地特征的，可以用来做工程分析。

2.3 数值模拟分析

位移分析采用SLIDE2018软件，其核心程序基于美国地质地震局的SLAMMER程序。首先利用SLIDE软件分析该护岸结构在地震作用下(PGA=0.44g)的整体稳定性，计算方法选用简化Bishop法。以某钻孔为例，地震工况下护岸结构的拟静力分析(地震系数kh=0.5 PGA)整体稳定系数FOS只有0.586，不满足FOS> 1的规定，因此，需要对地震下结构位移进行分析评价。对护岸临界地震系数进行模拟分析，计算得出kc=0.098。利用SLIDE软件自带的基于Newmark法的地震位移分析模块，输入上述拟合的16条地震波，土体模型采用等效线性，阻尼比取0.05，参考应变取值0.05，分别选取刚体算法、耦合算法以及解耦算法3种模型进行计算，得到不同地震波下每个滑动面的永久位移。SLIDE地震位移计算见图8。

图8 SLIDE地震位移计算

2.4 软件与简化公式计算结果对比

根据SLIDE计算得到的临界加速度以及地震波主要特征参数，按照上述7种位移简化公式进行计算，得到相应的最大永久位移值。简化公式位移计算参数为：峰值加速度平均0.44g，临界加速度为0.098g(最危险断面)，阿里亚斯烈度平均值 Ia=5.2，平均峰值速度PGV=69 cm/s，地震震级 Mw=7.6，滑动体高度为25 m，剪切波速vs=163 m/s，滑动体基本周期Ts=0.613 s,1.5 倍周期0.92 s对应反应谱Sa=0.66g。

将SLIDE软件计算和上述几种简化位移计算方法得到的结果进行分析拟合，对比关系见图9。

对16条地震波进行分析，尽管16条地震对应谱相同，但是波形相差巨大，不同方法的计算差异较大。根据分析可知，目前主流简化算法考虑的因素主要为统计值，Ambraseys & Menu 公式、Jibson 公式、NCHRP 611算法、Bray and Travasarou法、Makdisi & seed法，均无法反映具体地震波的破坏性。Jibson公式考虑了阿里亚斯烈度，计算结果可以明显反映不同地震波的区别，Shang-Yu Hsieh、Chyi-Tyi Lee对Jibson公式进行了修正，更好地反映了破坏趋势，与SLIDE软件计算结果更加接近。

图9 不同地震波下不同计算方法的位移结果对比

根据国内外相关文献，地震位移分析仍然以刚体分析为主。也有学者开展了耦合算法和解耦算法的研究，本次抗震设计时，除了采用刚体算法，同时考虑了耦合算法和解耦算法的对比计算，土体模型采用目前国内外学者主流采用的等效线性模型，阻尼比取值5%，参考应变取值0.05。分析可知，采用耦合算法或解耦算法，虽然部分结果更接近于目前常用的NCHRP 611方法和Makdisi & seed法，但是参数设定对结果影响巨大。采用刚体计算方法，与目前国际主流计算结果的大值基本一致，总体上偏于保守的。因此抗震位移分析，采用刚体算法作为基准，计算结论可信。

考虑到计算结果的离散性，本文采用位移计算结果的平均值作为工程设计的依据，即取每目标反应谱下16条人工地震波在同一临界加速度下滑块位移的平均值作为控制标准。根据计算结果，设计方案的稳定性满足本次抗震设计的目标。

3 基于SLIDE结果分析的简化计算公式

前文的7种简化算法虽然简单易用，但相对于SLIDE软件的精确积分计算而言，部分地震波工程不能很好地反映最终位移与地震参数之间的关系。Ambraseys & Menu 公式、Jibson公式、NCHRP 611算法、Bray and Travasarou法、Makdisi & seed法，可以计算某一地震反应谱基本参数下的地震位移，但是均无法反映具体地震波下的位移，对于相同反应谱不同地震波适应性较差。Jibso公式和Shang-Yu Hsieh公式虽然体现了阿里斯亚烈度，体现了不同地震波的差异，但是由于选取地震波参数的不同，有时误差较大，未能直接体现阿里亚斯烈度中超出临界加速度造成破坏的部分。

本文根据项目计算结果的归纳分析，综合以上简化算法的优点，提出一种基于有效阿里亚斯烈度和有效破坏时长，基本参数采用阿里亚斯烈度、临界地震系数、峰值地震系数的永久位移表达式：

(14)

式中：Dm为简化算法预估平均位移；A为经验系数；Ia为地震波阿里亚斯烈度；kc为临界地震系数；kmax为峰值地震系数；γ为经验系数，反映最终位移与临界地震系数和峰值地震系数的比值的因素；β为经验系数，反映同一地震波加速度超出临界加速度累积频率相关的因素。

对菲律宾项目不同水深(4.5～12.5 m)、滑弧深度在原泥面以下4～13 m、滑动体高度15～30 m的20钻孔断面近720组计算工况数据的计算结果进行拟合和分析总结。由于不同地震波离散性较大，拟合时，kmax取值为0.6 PGA/g+0.4 SMA/g，即0.6倍最大峰值地震系数与0.4倍第3大峰值地震系数之和。根据试算的分析情况，提出两组近似解：第1组A=5.484，β=3.0，γ=1.0;第2组A=6.221，β=3.0，γ=0.964，式(14)可以进一步写为：

(15)

(16)

对于不同钻孔、不同水深和不同计算断面的16条波位移平均值，大部分断面的SLIDE计算平均值与简化算法式(15)平均值误差在5%以内，个别钻孔相差10%左右；对于单一地震波，计算位移值的误差一般在25%以内，个别地震波误差1倍左右，如图10所示，主要原因为地震波比较特殊，PGA与SMA相差过大，地震波峰值悬殊。

图10 简化计算与SLIDE软件模拟对比

对于同一断面采用多条地震波平均值统计时，本文推荐的简化算法适应性较好，类似条件可以参照使用。

4 结论

1)国内规范对于地震工况下的稳定分析常采用拟静力法，而对于海外项目尤其是位于高震区的工程，从节约工程投资的角度出发，当地震情况下整体稳定系数小于规范要求时，可采用位移分析法对结构的安全性进行评价。采用位移控制的抗震设计方法，更加真实清晰地反映破坏程度，为节省高震区的护岸结构工程造价奠定了基础。

2)通过对国际上常用的7个地震位移简化公式和SLIDE软件的计算结果进行拟合对比分析，验证了利用SLIDE软件和其中部分简化公式进行地震位移计算的合理性和可靠性。

3)地震分析时，SLIDE软件推荐采用刚体算法，计算结果虽然相对偏保守，但与目前国际主流公式的计算结果基本一致。采用耦合算法或解耦算法时，土体模型的选择和参数选取至关重要，对计算结果比较敏感，在无试验支撑时，建议慎重选用。

4)在无地震波的情况下，可以依据地震反应谱的基本参数，使用Makdisi & seed法、NCHRP 611法、Bray and Travasarou法其中之一进行快速评估，但是具体到特定地震波可能会有较大偏差；有具体地震波参数时，可以采用Jibson、Shang-Yu Hsieh、Chyi-Tyi Lee公式对位移进行快速评估，有条件时也可采用Slide软件进行精确计算。