琼州海峡海床地震反应特性的一维非线性分析

陈国兴，夏高旭，王彦臻，金丹丹

(1. 南京工业大学岩土工程研究所，江苏，南京 210009；2. 江苏省土木工程防震技术研究中心，江苏，南京 210009；3. 江苏大学土木工程与力学学院，江苏，镇江 212013)

海南岛地处南海之滨，扼守太平洋与印度洋连接之要冲，与东南亚诸多国家地缘相近，是中国唯一的国家级全省经济特区。海南全省建国家级经济特区30 周年之际，国家赋予海南“全面深化改革开放试验区、国家重大战略服务保障区”的战略定位，是着眼于当前国际国内发展作出的重大决策。这意味着海南将在国家整体内外战略的实施中承载重大而特殊的使命。拟建的包头-海口高铁是中国“十纵十横”规划中的综合运输大通道，是一条实现南北融合、海陆互连、东西贯通的国家战略通道。因琼州海峡天堑，海南岛与内陆隔海相望。修建琼州海峡跨海通道，对海南自贸港建设、海南岛经济发展、海洋国防建设等具有重大的战略意义，对中国增强南海战略投送能力，履行南海维权、维稳、资源开发的重要使命具有不可替代的作用。

琼州海峡拟建跨海通道场址区域内发育了近东西向、北东向和北西向3 组断裂。北部陆域主要发育北东向断裂，南部海域主要发育北东东向、近东西向断裂。晚第四纪活动断裂主要位于琼州海峡附近，包括马袅—铺前断裂、海秀断裂、长流—仙沟断裂与琼华—莲塘村断裂。6 级以上地震主要沿区域中部北西向断裂与其他方向断裂交汇的地带发生，其发震构造为晚第四纪活动断裂或控制第四纪盆地的区域性早中更新世断裂。琼州海峡所在中南沿海地震带在整个活动周期(公元1400 年至今)内，共发生M≥ 4.7 级地震32 次，目前处于第二活动周期(1711 年至今)后期，未来遭遇概率较高的破坏性地震为震中距较近或适中的中强地震。拟建海底隧道沿线附近历史上曾数次发生中-强、强地震，如：1605 年琼山M7½级地震(震中距39 km)，1605 年湛江东南M6½级地震(震中距72 km)，1618 年澄迈M5½级地震(震中距17 km)，1871 年琼州海峡M5½级地震(震中距8 km)。其中，1605 年琼山地震主震造成的隧道沿线场地最大影响烈度达到Ⅷ度。因此，近场、中远场的中-强、强地震对拟建海底隧道安全存在严重的现实威胁。

经多年不懈研究，已积累了大量相关资料，如：海峡表层沉积物多为砂砾土[1]；海床软土及砂层互层沉积物深厚，易触发液化、震陷等地震地质灾害[2]；海峡中部四塘-天尾角线位是相对合适的盾构隧道线位[3]。孙田等[4-5]试验研究了琼州海峡海床钻孔84 个原状土样的动力特性，给出了海床100 m 以浅的砂类土、粉质黏土等土类的动剪切模量比和阻尼比，及海床100 m 以深的粉质黏土动剪切模量比和阻尼比。Chen 等[6]研究了福州河口盆地非线性地震效应特征；陈国兴等[7]研究了某海湾滨海平原的巨厚第四系地层非线性地震效应特征；Chen 等[8]、Ruan 等[9]研究了汕头海湾跨海剖面的非线性地震效应特征。这些成果对研究琼州海峡海床地震效应特征具有借鉴意义。

全长近30 km 的琼州海峡海底隧道是一项巨型工程。保证强地震时穿越琼州海峡海床复杂岩土介质的海底隧道的工程安全是一项重大的挑战，至今还未见琼州海峡跨海通道场址非线性地震效应研究的文献报道。本文以琼州海峡拟建跨海通道沿线海床为研究背景，选取4 个典型钻孔，基于通用软件ABAQUS，采用新近提出的non-Masing 本构模型[8,10]描述海床土的动力非线性特性，进行海床非线性地震反应分析，并与基于扩展Masing 类土体动力本构模型的专业软件DEEPSOIL V6.0 的计算结果进行对比，探讨土体动力本构模型的选取对琼州海峡海床非线性地震反应的影响及海床的非线性地震反应特征，期望能为琼州海峡跨海通道抗震设计提供科学依据。

1 海床地震反应的计算模型

1.1 场地模型

琼州海峡海底地貌为深水盆地，东西长约80 km，宽约30 km，海水等深线大于50 m。海床地形起伏大，土层分布不均匀，多有隆起与凹陷。

海床上部覆盖第四系土层：淤泥、粉砂、粉质黏土、中砂及粉砂互层等，厚约几十米；下部沉积第三系土层，多为厚层黏土，部分区域为粉砂质互层，厚达上百米。选取的典型钻孔ZK-06、ZK-08、ZK-11、ZK-13 的剖面信息如图1 所示。

图1 琼州海峡地层的典型钻孔和剪切波速剖面Fig. 1 Shear-wave velocity profiles and the soil lithology with depth in the typical boreholes, Qiongzhou Strait

ZK-06 孔深150.65 m，顶面标高-26.94 m，上部为粉质黏土、粉砂、中砂及下部粉质黏土土层。ZK-08 孔深120.20 m，顶面标高-55.60 m，上部含有较厚粉质黏土层，下部为粉质黏土和粉砂质互层；ZK-11 孔深200.20 m，顶面标高-84.80 m，上部为粉砂互层、粉质黏土、中砂及下部为厚的层状黏土分布，为海床土层条件的控制性钻孔。ZK-13 孔深120.50 m，顶面标高-81.60 m，钻孔剖面上部含有较厚细砂层，下部为厚层状粉质黏土。

一维场地地震反应分析时，将场地视为水平成层介质，输入地震动视为垂直入射的剪切波。

《中国地震动参数区划图》(GB 18306-2015)[11]规定剪切波速vs≥ 500 m/s 且其下不存在更低波速的岩土层界面可作为基岩。由图1 可知：4 个钻孔的终孔深度处均可作为地震动输入界面的基岩。

基于ABAQUS/Explicit 软件平台的土柱地震反应视为平面应变问题，采用二维水平成层模型计算一维土柱地震反应；土单元采用4 节点双线性插值平面应变缩减积分单元(CPE4R)；土柱底部和侧边界设置等效黏弹性人工边界单元[12-13]，土柱模型如图2 所示。基于ABAQUS 平台显式算法和并行计算技术，采用中心差分算法求解动力方程，利用每一时步起始时刻的平衡条件计算时步结束时刻的位移平衡解，运动方程解的显式表达式为：

图2 场地反应计算模型Fig. 2 Simulation model of site response

式中：u¨为加速度矢量；M为集中质量矩阵；F为施加的外荷载矢量；I为单元内力矢量；上标i为显式动力分析中第i步增量。为确保强非线性场地地震反应分求解的收敛性，显式动力分析的时步不超过10-5s。二维和三维非线性地震反应分析结果佐证了上述方法的可用性和通用性[6,9,14-16]。

采用DEEPSOIL 软件计算水平成层土层的地震反应时，土层简化为由质点、非线性弹簧和粘滞阻尼器组成的集中质量模型[17]，底部为弹性支座，如图2 所示。

土层的计算分层厚度h按下式确定[9]：

式中：vs为土层剪切波速；fmax为最大截止频率，本文取为25 Hz。

海底强震观测与理论研究表明，对表面水平的海床，海水对海床水平向地震反应的影响很小[18-19]。故本文未考虑海水与海床的动力耦合效应。

1.2 土体非线性本构模型

ABAQUS/Explicit 中缺乏土体动力非线性本构模型，故采用Chen 等[8,10]提出的基于Davidenkov骨架曲线、不规则加卸载准则构造的非曼辛(non-Mashing)滞回本构模型(简称DCZ 模型)描述土动应力-应变关系(图3)。Davidenkov 骨架曲线表达式：

图3 DCZ 模型描述的应力-应变关系示意图Fig. 3 Stress-strain curves of the DCZ Model under irregular loading-unloading-reloading

Chen 等[10]采用Fortran 语言编写用户自定义的DCZ 模型材料子程序VUMAT，实现了定义单元内高斯点材料状态的应力更新算法。DCZ 模型具有构造应力-应变滞回圈的加、卸载法则简单且精度高、计算效率高的优点，数值模拟的日本阪神地震液化场地井下台阵不同深度的水平向、竖向地震动与强震记录高度吻合，验证了该模型的有效性[10]。VUMAT 子程序已在海床场地、地下结构、海底隧道、软岩上核电厂结构等一系列二维和三维非线性地震反应分析中进行了有效的测试，佐证了该DCZ 模型的适用性[9,14-16]。

DEEPSOIL 采用Matasovic 骨架曲线及扩展Masing 法则构造不规则加、卸载的滞回本构模型(简称MKZ 模型)[17]，以描述土的动力非线性和滞回特性。Matasovic 骨架曲线表达式：

式中：A、B、β、s为无量纲常数。

Davidenkov 和Matasovic 骨架曲线参数可由拟合动剪切模量比G/Gmax与阻尼比λ 的试验曲线给出。孙田等[4-5]试验给出的琼州海峡各类土深度相关的G/Gmax与λ 的均值曲线，如图4 所示。

图4 琼州海峡海洋土的动剪切模量比与阻尼比曲线Fig. 4 Variations of shear modulus reduction and damping ratio curves of marine soils in the Qiongzhou Strait

1.3 输入地震动

由于缺乏工程场址及邻近地区历史地震记录，参考历史地震震级、震中距等信息，在日本Kik-net 网站选取3 次近场、中远场的中-强、强震记录NS 分量作为基岩输入地震动，详细信息见表1。依据海底隧道场址地震安全性评价结果，50 年超越概率63%(小震)、10%(中震)和2%(大震)的基岩峰值加速度(PBA)约为0.075g、0.20g和0.40g，g为重力加速度。据此调整强震记录的峰值加速度，以探讨PBA 对场址地震效应的影响。图5 为加速度时程及其傅氏谱、谱加速度(阻尼比5%，下同)、频域累积能量曲线。TCG012、IWT013 和HRS005 记录的5%～95%累积频域能量频段分别为3.4 Hz～9.0 Hz、2.4 Hz～17.8 Hz、4.3 Hz～18.0 Hz。TCG012 和HRS005 记录的频谱特征呈单峰特征，TCG012 记录呈窄频带特征、HRS005 记录呈宽频带特征；而IWT013 记录的频谱特征总体上呈双峰特征、宽频带特征。

表1 基岩输入地震动的原始地震记录信息Table 1 Information of original earthquake recordings for bedrock input motions

图5 输入基岩地震动的加速度时程、傅氏谱值谱、5%阻尼比谱加速度及频域累积能量比曲线(PGA = 0.1 g)Fig. 5 Bedrock motion (NS) time histories, Fourier spectra, 5% damping spectral accelerations, and cumulative distributions for frequency domain energy (0.1 g level)

2 海床地震反应的模拟结果与分析

2.1 场地基本周期及类别

海床30 m 以浅平均走时等效剪切波速vs30：

式中：vsi和di为海床30 m 以浅第i层土的剪切波速和厚度。钻孔ZK-06、ZK-08、ZK-11、ZK-13剖面的vs30分别为152.3 m/s、168.7 m/s、270.8 m/s、204.5 m/s。

采用HVSR 法[20]计算场地基本频率fs，取不同输入地震动计算的均值。以钻孔ZK-11 为例，计算结果如图6 所示。钻孔ZK-06、ZK-08、ZK-11 和ZK-13 剖面的fs分别为0.93 Hz、1.02 Hz、1.17 Hz 和1.16 Hz，相应的周期Ts分别为1.08 s、0.98 s、0.85 s 和0.86 s。由图1 可知：钻孔ZK-06、ZK-08、ZK-11 和ZK-13 表面至基岩vs= 500 m/s的覆盖层深度H分别为134 m、110 m、154 m 和87.6 m。据此，基于vs30、H和Ts的场地分类法[21]：钻孔ZK-06、ZK-08 剖面为Ⅳ类场地；钻孔ZK-11、ZK-13 剖面为Ⅲ类场地，但邻近Ⅲ类、Ⅳ类场地的界限，尤其钻孔ZK-13 剖面。

图6 场地基本频率的确定Fig. 6 Determination of site fundamental frequency

2.2 加速度传递函数

图7 给出钻孔ZK-11 海床地震反应的加速度传递函数[土层地震反应与基岩地震动的傅里叶幅值谱的比值(Fourier spectrum amplitude ratio, FSR)]随深度的变化。MKZ 模型计算的FSR 色谱图显示：距海床地表约80 m～100 m 以浅土层对0.3 Hz～0.8 Hz 频率地震波的放大效应显著。DCZ 模型计算的FSR 色谱图显示：约0.6 Hz～1.1 Hz 和1.5 Hz～2.0 Hz 为海床地震波传播放大效应显著的优势频带，距海床地表约35 m 以浅和深度80 m～165 m的土层对近1 Hz 频率地震波的放大效应尤为显著。DCZ 模型计算的FSR 色谱图在更宽的频带存在较多分散的、颜色较深的波纹，即FSR 值更大。

图7 不同基岩地震动激励时钻孔ZK-11 剖面的地震动加速度传递函数随深度的变化Fig. 7 Variation of the Fourier spectrum amplitude ratios with depth at the ZK-11 borehole profile subjected to various bedrock motions

这意味着：与MKZ 模型相比，DCZ 模型在深软土层中能传播更宽频带、较高频率的地震波。对相同频谱特性的基岩地震动，随着PBA 增大，较高频率的FSR 色谱颜色变浅，即地震波在土层传播过程中的放大效应减弱。究其原因，应为强震动激励下土的非线性和滞回效应更显著，导致土体刚度降低、阻尼增大，更多的地震波能量在土层中的传播过程中被吸收。

2.3 峰值加速度反应

图8 展示了钻孔ZK-11 处海床的峰值加速度(peak acceleration, PA)反应。由图可以看出：海床浅层的PA 放大显著，且海床的PA 反应随PBA 的增大而增大。总体上，DCZ 模型计算的PA 反应大于MKZ 模型计计算的PA 反应，两者计算的PA 反应沿深度的变化趋势，小震时基本相同、中震时相近、大震时差异显著。DCZ 模型计算的海床地表峰值加速度PGA 显著大于MKZ 模型计算的PGA，且两者的差异随PBA 的增大而增大。定义PGA 放大系数等于PGA/PBA，则两者计算的PGA 放大系数均随PBA 的增大而减小；当PBA 相同时，呈宽频带、双峰特征的IWT013记录激励时海床地表PGA 放大最小，能量集中在低频带的TCG012 记录激励下海床地表PGA 放大最大。在不同PBA、不同频谱特征的基岩地震动激励下，DCZ 模型计算的海床PGA 放大系数均大于MKZ 模型计算值，且在IWT013 记录激励下DCZ 模型和MKZ 模型计算的PGA 放大系数的差异最大。相应于0.075g、0.2g和0.4g的基岩地震动水平，DCZ 模型计算的PGA 分别为0.139g、0.297g和0.488g。由此可见，较之MKZ 模型，DCZ 模型能将更多的地震波能量传播到海床浅层土。

图8 不同基岩地震动激励时钻孔ZK-11 剖面的PA 随土层深度的变化Fig. 8 Variation of peak accelerations with depth at the ZK-11 borehole profile subjected to various bedrock motions

2.4 谱加速度反应

图9 为钻孔ZK-11 剖面谱加速度(spectrum acceleration, SA)反应沿海床深度的色谱图。可以发现，地震波从基岩经海床传播至地表，不同频率分量的地震波在土层的传播过程中存在复杂的滤波或放大效应。总体而言，海床浅层、中部以浅的SA 反应更强烈，呈现谱加速度放大现象；不同强度的IWT013 记录激励时，在周期0.1 s～0.2 s和0.3 s～0.5 s 区间的海床的SA 反应较大，这与IWT013 记录反应谱的双峰周期基本对应；PBA =0.4g时，频谱特性总体呈单峰特征的TCG012 和HRS005 记录激励时在周期0.1 s～0.3 s 内的海床SA反应较大。这表明：钻孔ZK-11 处海床SA 反应放大显著的优势周期段(频带)与基岩输入地震动的卓越周期和频域累积能量分布存在极高的相关性。相同频谱特性的基岩地震动激励下，海床的SA 反应随着PBA 的增大而增大。与MKZ 模型计算的SA 色谱相比，DCZ 模型计算的SA 色谱颜色较深部分的频带更宽而且颜色更深，尤其是近海床地表及深约35 m～85 m 土层。呈宽频带、双峰特征的IWT013 记录激励下，DCZ 模型计算的周期0.1 s 附近的SA 色谱明显比MKZ 模型计算的SA 色谱颜色更深，也即SA 反应要大。

图9 不同基岩地震动激励时钻孔ZK-11 剖面的谱加速度随深度的变化Fig. 9 Variation of spectrum accelerations with depth at the ZK-11 borehole profile subjected to various bedrock motions

图10 给出了PBA = 0.4g时海床地表的谱加速度曲线。对不同频谱特性的基岩地震动，DCZ模型均能很好地模拟周期小于0.1 s 的高频地震动在海床中的传播，呈现高频放大现象，但MKZ 模型模拟周期小于0.1 s 的高频地震动在海床中的传播时被显著过滤或无放大、过滤效应。中低频含量丰富、总体单峰频谱特征的TCG012 记录激励时，DCZ 模型计算的中-长周期0.1 s～1.0 s 的地表SA 值大于MKZ 模型的计算值。宽频带、双峰频谱特征的IWT013 记录激励时，MKZ 模型计算的周期0.3 s ～0.4 s 的地表SA 值显著大于DCZ 模型的计算值。宽频带、单峰特征的HRS005 记录激励时，除在邻近HRS005 记录的谱加速度卓越周期0.15 s 处，总体上，DCZ 模型计算的中-长周期0.15 s～1.0 s 的地表SA 值大于MKZ 模型的计算值，尤其周期0.15 s ～0.3 s 段的海床地表SA 值。

图10 不同地震动激励时钻孔ZK-11 剖面的地表谱加速度(PBA = 0.4 g)Fig. 10 Spectral accelerations at the ZK-11 borehole surface subjected to various bedrock motions with PBA = 0.4 g

文献[22]表明：DEEPSOIL 软件对软弱场地存在地震波高频率成分计算失真的现象，计算的地表峰值加速度明显低于实际记录。文献[10]基于DCZ 模型模拟了强震台阵液化场地的地震反应，在强震仪深度处的水平向和竖向地震动模拟和记录之间显示出很好的一致性。文献[23]将DCZ 模型用于拟强震台阵非液化和液化场地的地震反应分析，在强震仪深度处的水平地震动模拟和记录之间显示出很好的一致性，计算值与强震记录的误差明显小于DEEPSOIL 计算的结果。究其原因：DCZ 模型的不规则加卸载准则与时步很小的显式算法对高频地震动具有更强的传播能力。海床地表SA 值的大小与地震波经土层传递至浅表层土的高频分量有很大关系，如高频地震动传递充分，则海床地表高频段的SA 较大。综上，可合理推测，基于ABAQUS 平台的DCZ 模型计算的海床地表地震动更接近场地的真实地震反应。

2.5 海床地表设计地震动参数

习惯上采用峰值加速度(PGA)表征地震动强度，但PGA 未反映地震动的频率及持续时间的影响。近期研究表明[24]：累积绝对速度CAV (cumulative absolute velocity)与场地地震灾害严重性程度密切相关。CAV 定义为：

图11 给出不同PBA 的基岩地震动激励时海床地表的PGA 放大系数、地震动持时D5～95延长系数(地表/基岩地震动的D5～95之比)和CAV。DCZ模型计算的地表PGA 放大系数、D5～95延长系数和CAV 值均大于MKZ 模型的计算值。地震动的频谱呈单峰、窄带特征和中-低频丰富的近场地震TCG012 记录激励时，DCZ 模型计算的4 个钻孔的地表PGA 放大系数均为最大。地表D5～95延长系数随PBA 增大而增大；DCZ 模型计算的4 个钻孔的地表D5～95延长系数以TCG012 记录激励时最大、HRS005 记录激励时次之、IWT013 记录激励时最小。从地表PGA 和D5～95视角而言，该海底隧道海床地震效应最受近场地震控制。地表CAV随PBA 的增大而增大，且4 个钻孔的地表CAV，IWT013 记录激励时最大、HRS005 记录激励时次之、TCG012 记录激励时最小。究其原因，IWT013记录的有效持时D5～95最长、具宽频带特性；TCG012 记录的有效持时D5～95最短、具窄频频带特性；而HRS005 记录的有效持时D5～95、强震动段持续时间居中，也具有宽频带特性。从地表CAV 视角而言，该海底隧道海床地震效应最受中-远场震地震控制。

图11 不同地震动激励时海床地表峰值加速度PGA、地震动有效持时D5～95 和累积绝对速度CAVFig. 11 Peak ground accelerations (PGAs), 5～95% significant durations (D5～95), and cumulative absolute velocities (CAV) at the seabed surface subjected to various bedrock motions

虽然钻孔ZK-06、ZK-08 剖面为Ⅳ类场地、钻孔ZK-11、ZK-13 剖面为Ⅲ类场地，但这4 个钻孔剖面均处于场地类别分类界限的两侧，场地的总体特征差异不大。从图11 可以发现：输入相同的基岩地震动时，4 个钻孔的PGA 放大系数、D5～95延长系数和CAV 值比较接近。鉴于此，图12 给出3 个输入地震动水平的4 个钻孔ZK-06、ZK-08、ZK-11 和ZK-13 剖面的规准化地表加速度反应谱β 谱及均值线。作为对比，图12 中给出了《中国地震动参数区划图》中Ⅲ类、Ⅳ类场地的罕遇地震设计地震反应谱β 谱。可以发现：T= 0.05 s～0.5 s 的中、短周期段，海床场地相关反应谱β 谱明显大于区划图β 谱的取值。意味着在该周期段，《中国地震动参数区划图》的β 谱偏于不安全。

图12 规准化地表加速度反应谱β 谱Fig. 12 Normalized ground motion acceleration response spectrum β spectra at the seabed surface

3 结论

通过琼州海峡拟建海底隧道沿线4 个典型钻孔剖面的一维非线性地震反应分析，主要结论如下：

(1) 对比DCZ 模型和MKZ 模型计算的深厚海床非线性地震反应特征，结合文献[22]的结果，可以合理推测：与MKZ 模型相比，DCZ 模型能更好地模拟地震动的高频和中-长周期分量通过深厚海床土层的传播，DCZ 模型计算的海床浅层的地震反更接近场地的真实地震反应。近地表35 m以浅和深度80 m～165 m 土层对0.6 Hz～1.1 Hz 频带的地震波传播的放大效应尤为显著。

(2) DCZ 模型计算的海床峰值加速度PA 大于MKZ 模型的计算值，两者计算的PA 沿海床深度的变化，小震时基本相同、中震时相近、大震时差异显著。不论基岩地震动水平，DCZ 模型计算的海床地表PGA 放大系数、有效持时D5～95延长系数和累积绝对速度CAV 值均大于MKZ 模型的计算值。

(3) 相应于0.2g(中震)的基岩地震动激励，海床地表峰值加速0.297g，达到地震烈度8 度水平，这与场地历史最大影响烈度一致；CAV 较之PGA 是一个更合理的地震动强度指标；周期0.05 s～0.5 s 的海床场地相关反应谱β 谱大于《中国地震动参数区划图》反应谱β 谱的取值，依据该反应谱进行跨海通道抗震设计将偏于不安全。

鉴于海床地形地貌及海洋土沉积环境的复杂性，土体动力本构模型对海床场地二维非线性地震反应特性的影响，有待后续深入研究。