一种应用于低频舰船地震波数值模拟的联合吸收边界方法

本发明属于舰船地震波数值模拟领域，具体涉及一种应用于低频舰船地震波数值模拟的联合吸收边界方法。

背景技术：

1、舰船地震波数值模拟是一种通过将真实海洋环境简化为数学模型，并采用数值计算来模拟舰船地震波在模型中传播的方法。通过对各种模型的数值模拟，可分析舰船地震波在不同海底介质中的传播规律，从而为实际基于舰船地震波的水下远距离目标探测提供理论基础。在基于计算机实现的舰船地震波有限差分数值模拟中，受限于计算机的有限内存，需要人为界定计算区域，而简单地截断边界将会导致边界处产生强烈的反射波干扰，显著降低模拟精度。因此，通常需要对人工边界进行适当处理以消除计算边界处的虚假反射，从而使数值模拟结果更加符合实际真实情况。

2、单程波吸收边界、完全匹配层(pml)吸收边界是当前人工边界处理的两类主要方法。单程波吸收边界是在边界区域内使用单程波方程吸收入射波能量的边界条件方法，主要包括clayton-enquist(ce)吸收边界和higdon吸收边界两类。相比于ce吸收边界，higdon吸收边界通常具有更优的吸收效果，且其差分计算更为便捷。pml吸收边界则是通过在边界区域内添加完全匹配层，并在层内采用含有衰减函数的方程从而实现对虚假反射的吸收。该方法理论上能够吸收各个角度的入射波，因此近年来受到业界的广泛关注。孟路稳、朱旭芳等人先后采用pml吸收边界实现了舰船地震波有限差分数值模拟计算，在压制人工边界反射的基础上获得了较好的仿真结果。

3、然而，现有的pml吸收边界在实际应用过程中仍存在以下不足：其一，常规pml吸收边界选取的指数型衰减函数在内边界处增幅过快，这会导致中心波场与内边界的波场耦合效果不佳，从而产生内边界虚假反射；其二，入射波经过pml衰减层后通常不能被完全吸收，pml外边界处仍会产生部分残余边界反射。综上，当前常规pml吸收边界尚未能兼顾解决边界层中内、外边界处的虚假反射问题，其综合吸收效果仍有待于进一步提高。

技术实现思路

1、本发明要解决的技术问题在于提供一种应用于低频舰船地震波数值模拟的联合吸收边界方法。本发明在基于有限差分法进行舰船地震波数值模拟时，首先在pml边界区域内采用高斯型衰减函数，以优化中心波场与内边界的耦合效果，有效减少内边界反射；同时在人工边界的外边界处联合使用二阶higdon吸收边界，进一步吸收残余的入射波。相比于常规pml边界，本发明采用的联合吸收边界能够在不增加存储量的情况下综合压制内、外边界的虚假反射，从而显著提高舰船地震波的数值模拟精度。

2、本发明采取以下技术方案：

3、一种应用于低频舰船地震波数值模拟的联合吸收边界方法，包括以下步骤：

4、(1)在数值模拟过程中，首先将模拟区域划分为中心波场区域和人工边界区域。在中心波场区域内，将雷克子波作为舰船噪声声源，并采用交错网格有限差分法进行舰船地震波数值模拟；舰船地震波波动方程由二维一阶速度-应力弹性波方程来描述；

5、作为具体的实施方式，所述的二维一阶速度-应力弹性波方程其具体表达式如下：

6、

7、式(1)中，t为时刻，ρ为模型介质密度，w为激发声源，x、z分别表示空间坐标中的水平、垂直方向，vx、vz分别表示速度的水平、垂直分量，τxx、τzz、τxz则分别表示应力的三个分量，λ和μ为拉梅系数，其计算公式如下：

8、

9、式(2)中，vp、vs分别为纵波速度、横波速度。

10、(2)在中心波场区域内，对弹性波方程进行交错网格有限差分离散，采用时间2阶、空间2n阶精度的展开形式；

11、作为具体的实施方式，对弹性波方程进行交错网格有限差分离散表示为：

12、

13、式(3)中，δt为时间采样间隔，bm为空间差分系数，k为时间离散点序号，i、j分别表示空间坐标中x、z方向上的采样点序号，δx、δz分别表示空间网格步长。

14、(3)在人工边界区域内，添加l层完全匹配层作为吸收边界。为了优化中心波场与内边界处的耦合效果，进一步压制内边界反射，在pml层内采用高斯型衰减函数。作为具体的实施方式，所述高斯型衰减函数其表达式如下：

15、

16、式(4)中，d(h)为随距离变化的衰减函数的数值，vp为纵波速度，h为pml内计算点到pml最外层边界的距离，r为理论反射系数；将二维一阶速度-应力弹性波方程中各物理量在pml边界层内分解为水平分量和垂直分量两部分，并通过有限差分离散得到以下形式：

17、

18、式(5)中，i、j分别为空间水平、垂直方向上的离散点序号，k为时间离散点序号。vx、vz分别为速度的水平分量和垂直分量。σxx、σzz分别为水平、垂直方向上的正应力，σxz为剪应力，m表示空间序号。表示在空间位置(i+1/2,j)处第k个时刻的速度水平分量。表示在空间位置(i+1/2,j)处第k-1个时刻的速度水平分量。表示在空间位置(i+m,j)处第k-1/2个时刻的水平方向上的正应力。表示在空间位置(i-m+1,j)处第k-1/2个时刻的水平方向上的正应力。表示在空间位置(i+1/2,j+m-1/2)处第k-1/2个时刻的剪应力。表示在空间位置(i+1/2,j-m+1/2)处第k-1/2个时刻的剪应力。表示在空间位置(i+1/2,j)处第k个时刻的速度垂直分量。表示在空间位置(i+1/2,j)处第k-1个时刻的速度垂直分量。表示在空间位置(i+m-1/2,j+1/2)处第k-1/2个时刻的剪应力。表示在空间位置(i-m+1/2,j+1/2)处第k-1/2个时刻的剪应力。表示在空间位置(i,j+m)处第k-1/2个时刻的垂直方向上的正应力。表示在空间位置(i,j-m+1)处第k-1/2个时刻的垂直方向上的正应力。表示在空间位置(i,j)处第k+1/2个时刻的水平方向上的正应力。表示在空间位置(i,j)处第k-1/2个时刻的水平方向上的正应力。表示在空间位置(i+m-1/2,j)处第k个时刻的速度水平分量。表示在空间位置(i-m+1/2,j)处第k个时刻的速度水平分量。表示在空间位置(i,j+m-1/2)处第k个时刻的速度垂直分量。表示在空间位置(i,j-m+1/2)处第k个时刻的速度垂直分量。表示在空间位置(i,j)处第k+1/2个时刻的垂直方向上的正应力。表示在空间位置(i,j)处第k-1/2个时刻的垂直方向上的正应力。表示在空间位置(i+1/2,j+1/2)处第k+1/2个时刻的剪应力。表示在空间位置(i+1/2,j+1/2)处第k-1/2个时刻的剪应力。表示在空间位置(i+1/2,j+m)处第k个时刻的水平速度分量。表示在空间位置(i+1/2,j-m+1)处第k个时刻的水平速度分量。表示在空间位置(i+m,j+1/2)处第k个时刻的垂直速度分量。表示在空间位置(i-m+1,j+1/2)处第k个时刻的垂直速度分量。同理，表示在空间位置(i,j)处第k+1个时刻的水平速度分量。表示在空间位置(i,j)处第k+1个时刻的垂直速度分量。表示在空间位置(i,j)处第k+1个时刻的水平方向上的正应力。表示在空间位置(i,j)处第k+1个时刻的垂直方向上的正应力。表示在空间位置(i,j)处第k+1个时刻的剪应力。表示在空间位置(i+1,j)处第k+1个时刻的水平速度分量。表示在空间位置(i+2,j)处第k+1个时刻的水平速度分量。

19、(4)为吸收pml外边界处的残余入射波，在外边界处联合使用二阶higdon吸收边界。作为具体的实施方式，所述的二阶higdon吸收边界条件的左边界方程表达式为：

20、

21、式(6)中，ux、uz分别表示x和z方向上的波场分量，vp、vs分别为纵波速度、横波速度；采用talyor公式对式(6)进行有限差分离散，并将其应用于一阶速度-应力弹性波方程中，可得vx、vz、τxx、τzz、τxz的差分格式分别为：

22、

23、式(7)中，表示在空间位置(i,j)处第k+1个时刻的速度水平分量。表示在空间位置(i,j)处第k-l个时刻的速度水平分量。表示在空间位置(i,j)处第k+1个时刻的速度垂直分量。表示在空间位置(i,j)处第k-l个时刻的速度垂直分量。表示在空间位置(i,j)处第k+1个时刻的水平方向上的正应力。表示在空间位置(i,j)处第k-l个时刻的水平方向上的正应力。表示在空间位置(i,j)处第k+1个时刻的垂直方向上的正应力。表示在空间位置(i,j)处第k-l个时刻的垂直方向上的正应力。表示在空间位置(i,j)处第k+1个时刻的剪应力。表示在空间位置(i,j)处第k-l个时刻的剪应力。

24、式(7)中的系数可以通过以下恒等式得到：

25、

26、式(8)中，b通常为0.3～0.5之间的常数，vp为纵波速度，δt为时间采样间隔，δx为空间网格步长。右边界、上边界和下边界类比以上推导过程可得出。

27、本发明与现有技术相比的有益效果：

28、与常规pml吸收边界中使用的指数型衰减函数相比，高斯型衰减函数在内边界处数值变化更缓慢，这使得中心波场与内边界的耦合效果更佳，显著减少了内边界反射。此外，本发明提出的联合吸收边界是在采用高斯型衰减函数的pml吸收边界的基础上，再联合二阶higdon吸收边界对pml外边界处未衰减完全的入射波进行二次吸收，其可在不增加存储消耗的情况下，有效提高整体边界吸收效果。