无限元边界条件下的浅海海底地震波场建模

邓 峰，程广利，刘 宝

(1. 海军工程大学电子工程学院，湖北武汉 430033；2. 91001部队，北京 100089)

0 引 言

随着舰船减振降噪技术的改进和提高，舰船辐射噪声呈现低频化、甚低频化的特点，浅海水声信道的低频截止效应使得波长较长的 (甚) 低频噪声向海底耦合，因此利用舰船航行激发的海底地震波成为获取水中目标信息的一个重要途径[1]。舰船激发的海底地震波场主要包括水中简正波和地声[2]，其中 Scholte 波的能量主要聚集在海底界面附近，其在水平方向上衰减较慢，且不易受到海水水文条件的影响，在水雷引信设计、水中目标探测等军事领域具有较大的应用潜力[3]。

舰船航行激发的海底地震波波场的数值求解本质上是采用不同方法求解波动方程，主要包括三类方法。第一类是采用声场计算方法，文献[4]采用抛物方程结合反转算子法，对浅海弹性海底界面位移场进行求解，分析了不同海水环境和不规则海底对声传播和甚低频海底位移场的影响；文献[5]基于波数积分方法并通过快速场计算程序，分析了不同浅海环境下海底声压、位移和加速度的频率特性；文献[6]将传统的简正波模型进行改进后用于求解海底界面波的传播损失，与Kraken的计算结果一致性较好。文献[4-6]都是从传统的水声学模型出发研究地震波场，并没有对实际海底的弹性力学属性进行分析。第二类采用高阶交错网格有限差分法[7]，该方法的主要原理是将原始声场波动方程在时间和空间域上直接离散成差分方程组，配以相应的物理连续条件约束，通过数值计算可获得水中点声源激发的全波场解；海底边界的处理是该方法的一个研究难点，通常采用完美匹配层(Perfectly Matched Layer, PML)来实现，以模拟海底半无限空间海底，PML边界层将层中的介质阻抗与弹性介质阻抗进行匹配，波在吸收层中的传播按照添加衰减因子的波动方程来计算，通过合理设置衰减因子和PML层厚度实现入射波场的吸收[8]，PML边界层的数量决定了吸收效果，但层数的增加会大大降低整体计算效率。第三类是有限元法，该方法的基础是将物理域分成一些子域，在有限个自由度内求得子域中的精确解或近似解，其广泛应用于地质学、机械工程、水力工程等领域[9]，也可应用于舰船地震波场的求解。研究表明，能否精确地模拟实际中趋近无限空间的海底，对反射波的仿真结果影响较大，这使得对有限空间边界的处理成为海底地震波场建模研究的一个热点和难点；文献[10]在讨论不同激励源的海底地震波特性时未对边界进行处理，只研究了地震波未到达边界时的空间波场，从而回避了边界反射对地震波的影响；文献[11]也采用PML对边界进行处理，考虑其对波动成分的吸收能力有限，采取增大波场计算的水平距离以及Scholte波在水平方向的衰减系数的方法，以此实现Scholte波反射最小化的目的，同时对纵波和横波都施加衰减，以减小其对Scholte波的干扰；文献[12]使用黏弹性人工边界对海底边界进行处理，与波数积分法结果进行比对，总体变化趋势基本吻合，但频率较高时海底表面上的声传播损失存在较大误差。综上所述，从现有文献来看，关于海底地震波场边界的处理和研究并不多，且当前的处理结果并不能令人满意。

无限元法最早由Zienkiewicz提出[13]，是一种几何上趋近于无限远的单元类型，在使用时多与有限元相结合，因此无限元具有方向性。本质上，它是有限元的类型之一，作为有限单元的有效补充，能够更好地解决复杂的无限域问题[14]。经过多年的发展，该方法广泛应用于辐射声学、岩土力学、和电磁学等领域[15-17]，但尚未见到将其应用于浅海海底地震波场建模的研究中。

本文采用声学有限元和结构有限元构建海水-海底耦合模型，应用无限元方法对模型进行截断处理，模拟点声源在浅海中激发的海底地震波场，计算仿真海水中声压场的水平距离传播损失、海底地震波的应力场波场快照，与已有结论对比验证建模的正确性，最后分析了无限元边界的吸收效果。

1 基础理论

1.1 声压波动方程的有限元解

简谐声压p(x,t)的波动方程为[9]

其中： ∇2为拉普拉斯算子；P为声压振幅；k为波数。

对于无穷远问题还需声压满足索末菲(Sommerfeld)辐射条件：

以式(1)～(3)为基础，参照有限元理论[18]，常规等参单元上任一点的声压P和坐标x通过等参变换，可写成：

其中：Pe和Xe分别为单元节点的声压向量和坐标向量；N(ξ,η)为形函数。

由变分原理[18]可将方程(1)～(2)转化为动刚度矩阵方程：

其中：K和M分别为刚度矩阵和质量矩阵；P为整体的声压向量；F为载荷向量。

离散后的有限单元的各矩阵可表示为

1.2 弹性波动方程的有限元解

已知二维弹性介质中的波动方程为

其中：u和w分别为每个单元结点的水平和垂直位移；λ和μ为拉梅(Lame)常数。

单元内任意点的位移可以表示为

其中：Ue为单元结点上的位移矢量。

应变和位移的关系写成向量形式为

将式(10)、(11)代入拉格朗日泛函L后[19]，可得：

其中：B和D分别为无限单元的应变矩阵和材料本构矩阵。

由哈密顿原理在时间段[t1,t2]上对L积分，并使其变分为0，最终整理为

其中：Me和Ke分别为单元质量矩阵和刚度矩阵，Re为单元结点力矢量。

1.3 二维映射动力无限单元

为模拟无限大的海水-海底域，本文采用无限元边界截断无限大空间。声学无限元和弹性介质无限元原理[20]相似，这里以弹性介质的二维映射动力无限元为例简单介绍无限元理论。映射动力无限元的建立主要分为两步，一是建立单元在整体坐标系和局部坐标系的映射关系，二是构造位移形函数来确定单元在局部坐标系中的位移模式。

六结点映射动力无限元示意图如图1所示。构造六结点二维动力无限元的坐标系映射关系为

图1 六结点映射动力无限元示意图Fig.1 Schematic diagram of the mapped dynamic infinite element with 6 nodes

将无限元与有限元进行连接，根据两单元共同边上的位移连续条件，可以得到六结点二维动力无限元的位移形函数为

其中：Pj(ξ) (j=1 ,2,… , 6)是六结点二维动力无限元的波传播函数；Nj(ξ)(j=1 ,2,… , 6)是六结点二维动力无限元的位移形函数。

采用与一般有限单元相同的方法可以得到六结点二维动力无限元的单元质量矩阵和单元刚度矩阵为

2 基于无限元边界的浅海海底地震波场建模方法

为了构建浅海海底地震波场，建立点声源在海水-海底半无限空间中激发的地震波场模型，对海水和海底分别采用声学无限元和弹性介质无限元方法。海水-海底半无限空间模型示意图如图2所示。

图2 海水-海底半无限空间模型示意图Fig.2 Schematic diagram of seawater-seabed semi-infinite space model

图2中，声学无限元边界为虚线部分，其设置在海水层两侧，首先设置一个参考点作为声场分析的原点，然后对结构进行网格划分，将单元类型更改为声学无限单元，即完成了海水无限边界的设置。海底环境的无限元边界实现过程与声学无限元略有差异，实际中的海底为弹性基底，因此在建立海底地震波场模型时，将海底单元类型设置为平面应力单元，这种单元能较好地模拟一般弹性介质中的应力、应变、位移等物理量，以及弹性波在介质中的传播，在修改为无限元时不需要设置原点，只需将修正无限单元方向使其统一指向外域。

另外，在分析以Scholte界面波为主的海底地震波场时，能否准确地实现流体部分和固体部分的耦合，直接影响最终的仿真结果。在建模中，过“绑定”约束条件修正总体刚度矩阵，两界面间应力和位移连续，确保绑定区域不发生形变和相对位移。

由于平面应力单元的单元属性主要包括密度和弹性力学等参数，而非直接输入各类波速的数值。因此，需根据弹性力学参数和弹性波速之间的关系，计算获得各类波速。

假设海底底质是理想的弹性介质，根据虎克定律，无限虎克介质中的横波波速cs和纵波波速cp分别为

其中：E为弹性模量；υ为泊松比。

根据以上关系对常见的几类海底介质力学参数[21]和地声参数进行换算，结果如表1所示。

表1 不同海底介质类型的部分力学参数和地声参数Table 1 Mechanical and geoacoustic parameters of different seabed media

3 模型验证与特性分析

3.1 海底地震波场

与图2模型相似，在海水层中，海面为压力释放边界，海水的深度为400 m，海水中的声速c1为1500m·s-1，海水的密度ρ1为1000kg·m-3，海水左右两侧设置声学无限元来模拟无限元海域。海底地震波场的建模以砾岩和石灰岩两类基底为例，相关参数详见表1，整体仿真区域大小为800 m×800 m，网格大小为1 m×1 m，时间间隔为1 ms。声源位置设置在(400 m, 395 m)，声源信号采用频率为 15 Hz的雷克子波[7]。另外，为更好地反映海底地震波场的时空特性，本算例在时域上进行[22]。

砾岩海底、石灰岩海底(具体声学参数见表1)条件下本文方法仿真结果见图3(a)、图4(a)，中间黑色实线代表海水-海底分界面，界面以上部分为海水声压场，界面以下部分为海底地震波场，海底区域外阴影设置无限单元来吸收海底中各类弹性波。为了对比仿真结果，两种海底条件下，采用阶数为6、PML层数为10时的高阶交错网格有限差分法在0.2 s时的波场快照如图3(b)、图4(b)所示。图3、4中，“0”为声源位置，“1”为纵波，“2”为横波，“3”为泄露瑞利波，“4”为直达声波，“5”为反射声波，“6”为表面波，“7”为横波相关的侧面波，“8”为纵波相关的侧面波，图中S表示应力，S11表示11方向的应力。值得说明的是图3(a)、图4(a)中，海水-海底界面上以声压场为表征量，其下以应力场为表征量，两个物理量单位不同，故色棒不能统一，而图3(b)、图4(b)中，由有限差分法得到的仿真结果均以应力场表示。水中声源在弹性介质中能够激发出压缩波和剪切波，且这两种波彼此独立传播。由于界面反射产生的非均匀压缩波和非均匀剪切波相互作用形成表面波[4]，此类波在真空-弹性半空间界面中称为瑞利波(Rayleigh wave)，在液体-弹性半空间中称为Scholte波，而在液体-弹性介质环境中瑞利波依然存在，但由于上介质的反作用使其稍有改变，将其归为“泄露”波类。侧面波则由入射角大于全反射角的平面波在两种介质的界面处产生[1]。

从图3(a)、3(b)中可以清楚地看到，海底横波在分界面上的传播距离略小于水中声波的距离，我们称这种横波波速小于水中声波波速的介质为“软介质”，而类似于图4(a)、4(b)中横波波速大于水中声波波速的介质为“硬介质”。还可以看出，无限元边界法能清晰地反映出地震波的各组成成分，得到的波场成分波阵面与交错网格有限差分相同，且在同一时刻各波场成分到达位置相同，这也说明两种方法得到的波场成分传播速度一致，从而证明了本文方法的正确性。

图3 本文方法和高阶交错网格有限差分法计算得到砾岩海底地震波波场快照对比Fig.3 Comparisons between conglomerate-based seismic wave fields obtained by the proposed method and high- order staggered-grid finite difference method

图4 本文方法和高阶交错网格有限差分计算得到石灰岩海底地震波波场快照Fig.4 Comparisons between limestone-based seismic wave fields obtained by the proposed method and high- order staggered-grid finite difference

3.2 边界层吸收效果

以砾岩海底为例，说明无限元边界的吸收效果，建模方式、声源参数与环境参数与3.1节相同，计算的空间间隔为1 m，时间间隔1 ms，得到0.3 s和0.45 s时的波场快照如图5所示。在海水层中，与施加自由边界的海面相比，可以看到左右边界施加声学无限元边界后对水中声波吸收效果较好，几乎没有反射波存在。在海底层中，图5(a)箭头所示位置的坐标为(184 m, -737 m)，可以看到幅度较大的纵波到达边界后并没有出现明显的反射波，该点处应力幅值为0.01 Pa，占该时刻Scholte波幅值的0.89%，说明边界对纵波吸收效果较好。在图5(b)中能观察箭头处出现了较弱的横波反射波。读取两点的应力幅值，“1”点(214 m, -561 m)和“2”点(195 m, -760 m)的应力值分别为0.038、0.053 Pa，分别占该时刻Scholte波幅值的3.38%和4.71%，可见Scholte波反射的影响较小，可以忽略不计。

图5 在0.3 s和0.45 s时刻的无限元边界下砾岩海底地震波波场快照Fig.5 Conglomerate-based seismic wave fields obtained by the proposed method at 0.3 s and 0.45 s

此外，前文已提到文献[7]中有限差分法边界吸收效果与PML层数有关，且PML层数的增加导致计算效率降低，与高阶交错网格有限差分法相比，本文方法在不降低计算速度的基础上，能保证较好的边界吸收效果。为比较两者计算速度差异，在相同硬件环境下，对图3中仿真条件的三次运行时间进行比对，结果如表2所示。从表2中的结果可以看出，本文方法可大幅度缩短计算时间，提高运算效率。

表2 两种方法运行时间比对结果Table 2 Comparison of the running time of two methods

4 结 论

本文将无限元边界应用到浅海海底的地震波场建模中，对一般的浅海-海底半无限空间环境，将声学无限元施加于海水两侧模拟无限水平距离海域，将映射动力无限元施加于海底部分外围模拟无限大的海底区域，模拟了海底垂直正应力的波场快照，结果表明：对于浅海海底地震波场模型，无限元边界法仿真得到的垂直正应力波场快照能清楚地反映出波场各成分的存在，同时与高阶交错网格有限差分法仿真结果比对后，两者波场成分的波速一致性好。在海水层中，声学无限元边界可以较好地吸收声波；在海底层中，当海底纵波、横波抵达无限元边界时，对这两种波的吸收效果好，反射波幅值相较于Scholte波可忽略不计，这说明弹性介质无限元可以应用到浅海海底地震波场仿真中。另外，与高阶交错网格有限差分法相比，本文方法在保证结果准确性的前提下可大幅度提高计算速度。

真实水下目标与点声源存在较大差异，下一步将重点研究无限元边界下多点源、线源和体积源产生的海底地震波及其特性，分析不同声源类型对海底地震波场的影响。