引言

礼乐盆地有着巨大的油气资源潜力（孙龙涛等，2010；李颂等，2012），但其位于南海南部，处于环太平洋地震带上，地震活动频繁。为保障我国对海洋资源的开发，研究海底盆地的地震动规律具有重要意义。目前用来研究地震波在复杂介质中传播与散射的主要方法有：边界元法（Lee，2013）、有限元法（苏波等，2019）、有限差分法（Graves等，1998；姚铭等，2017）、伪谱法（谭文卓等，2020）以及谱元法（刘少林等，2021）。

谱元法由 Patera（1984）提出，它结合有限元法和伪谱法的思想，先在每个单元上进行谱展开，然后按照有限元的思想进行集成运算，因此该方法不但可以模拟任意复杂介质，而且具有较高的精度。相较边界元法，谱元法的使用条件更为宽松，无需求出问题的基本解；与有限元法相比（Antonietti等，2018），谱元法具有伪谱法的优势，即采用快速傅里叶变换，把波场中的空间域转化到波数域，避免了对空间求导的差分运算，可以有效减少数值频散和伪波；与有限差分法不同的是（彭浩天，2018），谱元法可以在模拟任意复杂形状的模型时，避免产生较多的散射波；比伪谱法更方便的是，由于具有有限元法的特点，所以谱元法可以自动满足自由边界条件。因此本文采用谱元法进行海底盆地地震响应的研究。

在流固耦合方面，李伟华等（2003）和赵成刚等（2008）给出了不同场地模型的波动解析解，但一般只能应用于简单规则地形。Komatitsch等（2000）将谱元法引入场地地震动的流固耦合分析中，开发了相应的显式并行计算软件Specfem 2D 与Specfem 3D。孔曦骏等（2021）使用软件Specfem 2D研究了几何地形、水深、地震波卓越频率以及土体性质对江水地震动水压力的影响规律。以上学者对流固耦合问题进行了大量研究，但并未对考虑流固耦合作用的海底盆地进行深入研究。在边界条件的处理上，采用开源软件Specfem 2D的内置完美匹配层PML（Perfectly Matched Layer）边界（董阳，2016），PML的思想是在模型区域外部加1层吸收层，将模型区域内部和吸收层内的波动方程耦合起来，使得在接触面处没有反射，并在吸收层上使用指数形式的衰减层，当波场传播至吸收层时会被衰减，不会发生反射。

大量研究表明，震源入射角对场地地震响应具有显著影响，丁海平等（2017）采用有限元法，对凹陷地形的地震动峰值放大系数进行了研究，发现放大系数受震源入射角的影响较大，但其并未考虑到盆地内软土覆盖层的影响。朱重洋等（2016）采用有限元法，在考虑覆盖层影响的情况下，研究了震源入射角度对放大系数的影响，但其仅考虑了1层覆盖层，对于覆盖层的研究依然不够深入。魏成前（2021）在考虑盆地内多层覆盖层的情况下，采用二维有限元法研究了震源入射角度对盆地地震动放大特征的影响，发现这些因素的影响较为显著，但其盆地模型为规则的梯形，并不能完全符合真实情况。以上学者在研究震源入射角对地震响应的影响时，建立的模型都比较理想化，不能足够准确地反映真实情况，因此本文使用Fortran语言编写了针对开源软件Specfem 2D的数据处理程序，以南海礼乐盆地为背景，建立了考虑海水作用和盆地内不规则地形以及多层介质影响的盆地模型，研究了震源入射角对盆地地震响应的影响，本文结论可为实际工程分析与设计提供有效的参考。

1.   计算方法以及边界条件

1.1   谱元法理论

谱元法的基本思想是先将计算区域划分为有限个单元，然后在每个单元上使用伪谱法，通过在每个单元中配置不均匀的分布节点，把单元的近似解表示为截断的正交多项式展开式，最后采用Galerkin方法求解得到全局的近似解，从而完成对区域内波传播过程的数值模拟（万子轩，2020）。

波动方程可以表示为：

式中，$ \rho $ 是密度；$ \boldsymbol{u} $ 是待求解的位移向量；${\boldsymbol{\nabla}}$是梯度算子；$ \boldsymbol{\sigma } $是应力张量，$ \boldsymbol{\sigma } $与位移的梯度呈线性关系；$ \boldsymbol{f} $ 为震源项。

在数值求解时，首先需要对式（1）中的微分形式进行积分。积分的具体方式为将式（1）乘以权函数$ \boldsymbol{w} $，并在整个求解区域上进行分步积分，从而得到相应的积分方程：

式中，$ \varOmega $是计算区域；$ \varGamma $是其边界；$ \boldsymbol{t} $为边界上的力。

与有限元法一样，谱元法在求解时，需要将计算区域离散为一个个单元，对于单元内待求的未知量u，使用一系列在基点上的$ {u}_{i}=u\left({\xi }_{i},{\eta }_{i}\right) $，$ i $=1, …, n和一系列的形函数$ {N}_{i}\left(\xi ,\eta \right)，i $=1, …, n来表示。

映射单元有2套坐标系。1套是局部坐标系$ \left(\xi ,\eta \right) $，另1套是全局的笛卡尔坐标系$ \left(x,y\right) $，坐标系$ \left(x,y\right) $和坐标系$ \left(\xi ,\eta \right) $通过坐标映射关系进行转换。

局部坐标系和整体坐标系的转换关系为：

式中，$ \boldsymbol{J} $为雅可比矩阵。

面积元的计算：

由于谱元法基于对微分方程的积分形式来处理，因此刚度矩阵、阻尼矩阵、质量矩阵形式均为：

式中，S为计算区域；$ \boldsymbol{G} $表示刚度矩阵、阻尼矩阵或质量矩阵。

积分区间选取[−1，1]内的点$ {\xi }_{1} $, $ {\xi }_{2} $, …, $ {\xi }_{n} $和权重值$ {H}_{1} $, $ {H}_{2} $, …, $ {H}_{n} $，使得：

谱元法的取值点$ {\xi }_{i} $（i=0, …, n）是Lagrange插值多项式中所选取的节点，同时满足n+1阶Legendre多项式的根：

该取值点被称为Gassus-Lobatto-Legendre点（简称“GLL点”）。$ {P}'_{n} $是n阶Legendre多项式的导数。GLL点具有正交的特点，所以质量矩阵为：

式中，${\boldsymbol{\varOmega }}_{\bf{e}}$为计算区域。

最终得到：

式中，$ \boldsymbol{M} $代表质量矩阵；$ \boldsymbol{K} $代表刚度矩阵；$ \boldsymbol{F} $代表外力。

1.2   考虑流-固耦合影响的运动方程

由于Legendre谱元法形成的质量矩阵为对角阵形式，便于建立显式动力算法，故在大型波动分析中应用广泛。采用处理流体动力方程的速度势格式，基于谱元法的流-固整体质量矩阵仍然保持对角阵形式，运动方程可表示为（孔曦骏等，2021）：

式中，$ \boldsymbol{U} $为固体位移矢量；$ \varPhi $为流体速度势；速度矢量$ \boldsymbol{V}=\nabla \varPhi $；$ {\boldsymbol{M}}_{\mathbf{s}} $为固体质量矩阵；$ {\boldsymbol{M}}_{\mathbf{f}} $为流体质量矩阵；$ {\boldsymbol{D}}_{\mathbf{s}} $为固体吸收边界相关的矩阵；$ {\boldsymbol{D}}_{\mathbf{f}} $为流体吸收边界相关的矩阵；$ \boldsymbol{A} $为流固2种介质的耦合矩阵，这是由交界面的位移和应力的连续条件产生；$ {\boldsymbol{K}}_{\mathbf{s}} $为固体刚度矩阵；$ {\boldsymbol{K}}_{\mathbf{f}} $为流体刚度矩阵；$ {\boldsymbol{F}}_{\mathbf{w}\mathbf{s}} $为固体或流体中波源的等效荷载矢量。

流体压力$ p $与速度势$ \varPhi $的关系为：

式中，$ \rho $为流体密度。

自由液面处的压力为0，该处的边界条件为：

流固交界面的法向速度和法向力连续：

式中，$ \boldsymbol{n} $为流固交界面处法向单位向量；$ {\boldsymbol{\sigma }}_{\mathbf{s}} $为固体介质中的应力矢量。

1.3   PML边界条件

为了避免内部入射波在界面发生反射，本文引入了PML（完全匹配层）吸收边界条件。如图1所示，根据PML原理沿水平方向x引入衰减因子$ {d}^{x} $；沿竖直方向y引入衰减因子$ {d}^{y} $。$ {d}^{x} $和$ {d}^{y} $一般选择为（彭浩天，2018）：

图 1  PML吸收边界示意图

Figure 1.  Schematic diagram of PML absorption boundary

下载: 全尺寸图片 幻灯片

式中，$ {x}_{i} $、$ {y}_{i} $分别为匹配层区域与内部区域界面的横向距离、纵向距离。$ {v}_{\mathrm{m}\mathrm{a}\mathrm{x}} $为最大纵波速度值；L为匹配层厚度，在不考虑流固耦合的情况下，通常PML厚度为半个波长时，则可达到满意的吸收效果。一般$ \alpha $=1e-6，$ a= $0.25，$ b= $0.75（彭浩天，2018）。

1.4   方法验证

本文选取Cristini等（2012）中基于谱元法的水声传播数值模型进行验证，并研究了网格大小对计算精度的影响。计算简图如图2所示，模型长8000 m，海底深2400 m，海水深400 m，声源作用在点M（400，395）处，接收器位于点N（800，380）处。假设海水为无黏、无旋的理想流体，海底介质为均匀和各向同性材料。海水材料参数为：密度$ {\rho }_{1}= $1000 $ \mathrm{k}\mathrm{g}/{\mathrm{m}}^{3} $，压缩波声速$ {V}_{{\rm{p}}1}=1\,500\;\mathrm{ }\mathrm{m}/\mathrm{s} $。海底材料参数为：密度$ {\rho }_{2}= $2000 $ \mathrm{k}\mathrm{g}/{\mathrm{m}}^{3} $，压缩波声速$ {V}_{{\rm{p}}2}=2\,400\;\mathrm{m}/\mathrm{s} $，剪切波声速$ {V}_{{\rm{s}}2}=1\,200\;\mathrm{m}/\mathrm{s} $。模型上侧采用自由边界，左侧、右侧和下侧边界采用PML吸收单元。声源选择主频率为30 Hz的Ricker波，波长λ为40 m，波源时程曲线如图3所示。模型1、2、3网格大小分别为5 m、10 m和20 m，边界均为PML边界。计算得到N点的压强曲线，将其峰值取为1，并与Cristini等（2012）结果进行对比。由图4可以看出，当网格尺寸为λ/2（20 m）时（模型3），计算结果存在误差，误差主要为曲线的不规则波动；当网格尺寸为λ/4波长（10 m）时（模型2），计算结果与Cristini等（2012）结果比较接近，也与网格尺寸为λ/8（5 m）时（模型1）的结果一致，说明本文所用研究方法是正确的，且考虑流固耦合问题时，网格尺寸为λ/4时即可取得较高精度。

图 2  验证算例计算简图（单位：米）

Figure 2.  Sketch of calculation example for verification（Unit: m）

下载: 全尺寸图片 幻灯片

图 3  波源时程曲线

Figure 3.  Source time-history curve

下载: 全尺寸图片 幻灯片

图 4  压强对比

Figure 4.  Pressure contrast

下载: 全尺寸图片 幻灯片

1.5   PML边界参数分析

为研究PML边界吸收效率和厚度的关系，采用图5所示计算模型，该模型除了震源、观测点位置以及材料参数外，其余参数均与图2一样。震源作用在点M（4000，1400）处，以集中力的方式施加，观测点为点N（5000，1400），海水材料参数为：密度$ {\rho }_{{\rm{w}}}= $1000 $ \mathrm{k}\mathrm{g}/{\mathrm{m}}^{3} $，压缩波速$ {V}_{{\rm{pw}}}=1\,500\;\mathrm{ }\mathrm{m}/\mathrm{s} $。海底材料参数为：密度$ {\rho }_{{\rm{s}}}= $2000 $ \mathrm{k}\mathrm{g}/{\mathrm{m}}^{3} $，压缩波速$ {V}_{{\rm{ps}}}=1\,000\;\mathrm{m}/\mathrm{s} $，剪切波速$ {V}_{{\rm{ss}}}=200\;\mathrm{m}/\mathrm{s} $。输入主频率为2 、2.5 、3 Hz的Ricker波，波长分别为100、80、67 m，网格大小为10 m，时程曲线如图6所示。PML边界厚度h分别取60、90、120、150 m，得到观测点的加速度时程曲线如图7所示。从图中可以看到，当h为60 m时，会存在较多反射波，导致加速度不断增大，随着h增大，反射波逐渐减少，直至h为120 m时反射波消失。这说明在考虑流固耦合的情况下，PML边界厚度为2个波长时方可达到理想吸能效果。

图 5  计算简图（单位：米）

Figure 5.  Calculation sketch（Unit: m）

下载: 全尺寸图片 幻灯片

图 6  波源时程曲线

Figure 6.  Source time-history curve

下载: 全尺寸图片 幻灯片

图 7  PML吸收效率对比

Figure 7.  Comparison of PML absorption efficiency

下载: 全尺寸图片 幻灯片

2.   礼乐盆地模型的建立

礼乐盆地模型如图8所示，地形采用张田升等（2019）对礼乐盆地的勘探数据，并将模型等比例缩小10倍。介质采用Chen等（2021）的数据，如表1所示。假定海水为无黏、无旋的理想流体，海底介质为均匀性和各向同性材料。模型长5000 m，海底介质分别为黏土、砂土和岩石，厚度均为60 m，海水深度为0 ～143.75 m，震源坐标分别为I（2500，500）和G（0，500），以集中力的方式施加，震源从盆地中部入射时，震源为I点，震源从盆地左侧入射时，震源为G点，震源入射角为α，分别取15°、30°、45°和60°。模型最大网格尺寸为10 m，时间步长为0. 2 ms。模型上侧采用自由边界，左侧、右侧、下侧均为PML吸收边界。流固耦合面上布置122个观测点，用以观测流固耦合面上不同位置的地震动峰值放大系数，取A、B、C、D和E 5个特征点，用以观测加速度时程以及频谱。震源分别采用Friuli波、Imperial_Valley波和Ricker子波，输入P_SV型地震波，时程以及频谱曲线如图9所示。

图 8  礼乐盆地计算简图（单位：米）

Figure 8.  Calculation sketch of the Lile basin（Unit: m）

下载: 全尺寸图片 幻灯片

表 1  模型介质参数

Table 1.  Model media parameters

介质	$ \rho $/（kg·m−3）	Vp/（m·s−1）	Vs/（m·s−1）
海水	1000	1500	—
黏土	1650	1650	218
砂土	1800	1697	264
岩石	2100	2135	485

下载: 导出CSV 

| 显示表格

图 9  震源时程以及频谱

Figure 9.  Source time histories and spectrum

下载: 全尺寸图片 幻灯片

3.   震源入射角对海底盆地地震响应的影响

取水平地震动峰值放大系数$ \beta = $|$ {{A}}_{\mathrm{m}\mathrm{a}\mathrm{x}} $/$ {{A}}_{\mathrm{m}\mathrm{a}\mathrm{x},\mathrm{i}\mathrm{n}\mathrm{p}\mathrm{u}\mathrm{t}} $|，即观测点地表水平地震动峰值$ {{A}}_{\mathrm{m}\mathrm{a}\mathrm{x}} $与输入点水平地震动峰值$ {{A}}_{\mathrm{m}\mathrm{a}\mathrm{x},\mathrm{i}\mathrm{n}\mathrm{p}\mathrm{u}\mathrm{t}} $比值的绝对值（朱重洋等，2016）。本研究编写了提取观测点地震响应峰值的Fortran程序，该程序可以批量读取结果文件，并按照观测点顺序提取观测点地震响应峰值，有效减少了工作量。

3.1   地震波由盆地中部入射情形

3.1.1   时域放大特征分析

图10为不同震源入射情况下观测点地震动峰值放大系数，可以发现：（1）震源的不同仅影响放大系数的大小，而不会影响其变化规律。不论采用何种震源，盆地中部的放大系数均为最大，盆地两侧的放大系数则最小，二者最大相差约18倍。盆地内的观测点放大系数分布曲线变化剧烈，这可能是由于模型考虑了流固耦合作用、盆地内大量凹凸地形以及多层介质的影响，导致地震波在传播过程中进行了多次反射和折射，从而使放大系数分布曲线变化剧烈，这是以往研究中没有发现的问题。图10（a）中最大放大系数峰值为0.68，图10（b）中最大放大系数峰值为0.62，图10（c）中最大放大系数峰值则为1.8，图10（c）中的放大系数峰值分别是图10（a）、图10（b）中放大系数峰值的2.65倍、2.9倍，这是因为图10（a）和图10（b）工况下的地震波作用时间较长，达到30 s，因此由盆地表面反射回来的反射波，会与后续震源部位的入射波产生相互作用，导致输入点位置的加速度峰值变大，根据本文对放大系数的定义可知，这种影响会使图10（a）和图10（b）的放大系数较图10（c）小。（2）随着震源入射角的增大，放大系数逐渐减小，该结论与丁海平等（2017）的结论一致。且入射角每增加15°，放大系数减少20%～30%。

图 10  观测点放大系数分布

Figure 10.  Amplification coefficient distribution of observation points

下载: 全尺寸图片 幻灯片

3.1.2   特征点的频谱分析

图11为特征点加速度时程以及频谱的对比图，可以看出：（1）在不同震源入射角下，盆地不同位置的加速度时程差异较大。在盆地左侧（A点），入射角为60°时，加速度峰值具有滞后性，且滞后3 s左右；而在盆地其他部位，不同入射角下的加速度则在同一时间达到峰值，且越靠近盆地中部，加速度到达峰值的时间就越短。（2）盆地不同位置的频谱特征具有明显差异。在盆地边缘部位（A点与E点），谱曲线的波动比较强烈，震源入射角的改变不仅会改变谱峰的大小，也会改变谱曲线的变化规律，且当震源为Ricker波时，这种改变更为明显。这与禹乐等（2020）使用Ricker波进行研究而得出的结论一致。而在盆地中部（C点），谱曲线的波动则更具规律性，震源入射角的改变不会使谱峰位置发生改变，仅会影响谱峰大小，且谱峰随震源入射角的增大而增大。（3）震源对不同震源入射角下的频谱特征具有较大影响。当震源为Friuli 波时，在盆地边缘部位（A点与E点），有效幅值在25 Hz 之内，而在盆地中部（C点），则30 Hz 以上的幅值才可忽略；当震源为Imperial_Valley 波时，仅有A点处25 Hz 以上的幅值可以忽略，其余盆地部位的幅值均在33 Hz 左右才会显著减小；当震源为Ricker 波时，有效幅值的频率区间则无明显变化。（4）结合禹乐等（2020）的相关研究可以发现，由于考虑了海水作用、盆地内不规则地形以及多层介质的影响，特征点的加速度时程存在剧烈变化，这是由地震波的多次反射和折射导致的。

下载: 全尺寸图片 幻灯片

图 11  特征点加速度时程以及频谱

Figure 11.  Acceleration time histories and spectrum of characteristic points

下载: 全尺寸图片 幻灯片

3.2   地震波由盆地左侧入射情形

3.2.1   时域放大特征分析

图12为不同震源入射情况下观测点地震动峰值放大系数，与地震波从中部入射的区别如下：（1）当地震波由左侧入射时，放大系数最大值均位于盆地左侧靠近盆地边缘部位，且随着入射角的增大向左偏移，而放大系数最小值则位于盆地右侧，二者最大相差14倍；图12（a）、图12（b）、图12（c）中的最大放大系数峰值分别为0.72、0.88、1.4，图12（c）中的放大系数峰值分别是图12（a）和图12（b）放大系数峰值的1.94倍、1.75倍。（2）随着震源入射角的增大，放大系数逐渐增大。与地震波由中部入射情形不同的是，入射角每增加15°，放大系数减少17%～30%。

图 12  观测点放大系数分布

Figure 12.  Amplification coefficient distribution of observation points

下载: 全尺寸图片 幻灯片

3.2.2   特征点的频谱分析

图13为特征点加速度时程以及频谱的对比图，与地震波从中部入射的区别如下：（1）无论在盆地何位置，加速度时程及谱曲线均具有规律性，震源入射角的改变不会使二者峰值位置发生改变，仅会影响峰值大小，且加速度峰值和谱峰均随着震源入射角的增大而减小。越靠近盆地左侧，加速度峰值到来的时间越短。（2）震源对不同震源入射角下的频谱特征具有较大影响。当震源为Friuli 波时，在盆地左侧（A点），有效幅值在32 Hz之内，越靠近右侧的部位，有效幅值区间越小，在盆地右侧（E点），有效幅值区间在20 Hz之内；当震源为Imperial_Valley 波时，在盆地左侧（A点），有效幅值在33 Hz之内，而在盆地右侧（E点），则在20 Hz之内；当震源为Ricker 波时，有效幅值的频率区间无明显变化。

下载: 全尺寸图片 幻灯片

图 13  特征点加速度时程以及频谱

Figure 13.  Acceleration time histories and spectrum of characteristic pointsts

下载: 全尺寸图片 幻灯片

4.   结论

本文运用谱元法理论建立了礼乐盆地地震动分析二维模型，研究了震源入射角对盆地地震响应的影响。得到以下结论：

（1）在考虑流固耦合的情况下，网格尺寸为四分之一波长时即可取得较高精度，PML边界厚度至少为2个波长时方可达到理想吸能效果。

（2）当地震波从盆地中部入射时，随着震源入射角的增大，放大系数逐渐减小，且入射角每增加15°，放大系数就减少20%～30%；盆地中部的放大系数最大，盆地两侧的最小，二者最大相差18倍左右；震源的不同会影响放大系数的大小，而不会影响其变化规律，但是对频谱特征具有较大影响；在不同的震源入射角下，盆地不同位置的加速度时程以及频谱特征差异较大。

（3）当地震波由左侧入射时，随着震源入射角的增大，放大系数逐渐减小，且入射角每增加15°，放大系数减少17%～30%；放大系数最大值均位于盆地左侧，靠近盆地边缘部位，随着入射角的增大而向左偏移，而放大系数最小值则位于盆地右侧，二者最大相差14倍左右；与地震波从中部入射相比，放大系数最大值变化了6%～30%；震源的不同会影响放大系数的大小以及频谱特征；在不同的震源入射角下，盆地不同位置的加速度时程以及频谱特征差异较小。

（4）由于考虑了海水、盆地内大量不规则地形以及多层介质的影响，地震波在传播过程中进行了多次反射和折射，致使盆地内的观测点放大系数分布曲线以及特征点的加速度时程曲线均存在剧烈变化。因此，在研究盆地地震相关问题时，这些因素的作用是不可忽视的。在后续的工作中，会对相关因素进行更加全面深入的研究。