获得稳定的近最佳匹配层吸收边界条件的方法与流程

本发明是应用在山前带等复杂地区的高精度地震传播数值模拟领域，具体可以为间断galerkin有限元等地震波数值模拟算法提供良好且稳定的吸收边界条件。

背景技术：

间断galerkin有限元法(dgfem)是一种高阶高精度有限元法，特别适用于复杂地质条件(尤其是剧烈地形起伏条件下)的地震波传播数值模拟。地震波数值模拟一般模拟有限区域内地震波的传播，因此需要引入人工吸收边界条件来降低来至边界的反射。在应用间断galerkin有限元地震数值模拟算法时，单元边界处的波场流动可以分为两部分：流出单元部分和进入单元部分，进入单元部分由相邻单元控制。因此，可以通过将计算区域边界处数值流的进入单元部分设置为零来控制边界反射，该吸收边界条件被称为开放边界条件。开放边界条件可以吸收一定的边界反射波，为了进一步提高吸收效果，可以在边界处外加一定厚度的阻尼吸收层，该组合吸收边界条件对反射波有很好的吸收效果，但对于面波、倏逝波等并不理想。对吸收边界条件作进一步改善需要采用最佳匹配层(pml)吸收边界条件。

pml目前主要有分裂和不分裂两种完成方式。berenger最早将pml运用于吸收电磁波。该pml技术将磁场分裂为垂直和平行于边界的两个分量，然后对每个分量运用吸收边界条件。从编程角度看，由于要考虑众多pml区域，编程较为复杂。从数学角度看，分裂式pml是弱适定的，容易引起虚假反射等不稳定性现象。对pml技术，理论上入射波可以无反射的进入该层，并被完全吸收，但经过空间离散后，该方法依然存在来自起飞角处的边界反射。为了提高空间离散后pml对来自起飞角处边界反射的吸收效果，kozougluandmittra(1996)引入了复频移(cfs)拉伸坐标变换，并采用了不分裂式pml。基于复频移拉伸坐标变换的pml，对于面波、倏逝波等有良好的吸收效果。不分裂pml不需要分裂波场，但在时间域需要做大量褶积运算，计算量较大。komatitschandmartin(2007)发展了一种无需显式褶积计算、基于迭代格式的不分裂pml技术——cpml(convolutionalperfectmatchedlayer)，大大降低了不分裂pml的时间域计算代价。v.etienne等(2010)最先考虑将cpml与dg-fem相结合模拟弹性波传播，并取得了良好的吸收效果。

cummer(2003)介绍了一种新的不分裂pml技术——近最佳匹配层(npml/nearlyperfectlymatchedlayer)。npml的吸收效果与标准的pml相当，但相比标准pml，npml在完成上更加简洁，且更省内存。薛昭等(2013)基于近最佳匹配层(npml)的思想和复频移拉伸坐标系变换(cfs)，通过定义新的辅助变量，推导了一种新的最佳匹配层技术(cfs-npml)，cfs-npml可以方便地融合到弹性波传播的间断galerkin有限元地震波传播数值模拟中。对比cpml和cfs-npml在dg-fem模拟弹性波传播中的表现也发现，二者都能有效吸收面波等近平行入射界面的地震波。在同样吸收参数条件下，两者吸收效果相当，但cfs-npml明显具有更高的计算效率。但结合间断galerkin有限元地震数值模拟算法使用cfs-npml时，在长时间模拟的情况下，算法存在明显的来至边界的不稳定性现象，这阻碍了cfs-npml的应用。

技术实现要素：

本发明的目的在于克服在间断galerkin有限元地震波数值模拟算法中使用基于复拉伸坐标变换的近最佳匹配层吸收边界条件(cfs-npml)时存在的长时间模拟应用不稳定现象，提出了获得稳定的近最佳匹配层吸收边界条件的方法，从而为复杂地区地震波高精度有限元数值模拟提供提供稳定良好的吸收边界条件。

本发明的目的及解决其技术问题是采用以下技术方案实现的：

一种获得稳定的近最佳匹配层吸收边界条件的方法，包括以下步骤：

一、将时间域的弹性波一阶速度-应力方程变换到频率域；

二、对频率域弹性波一阶速度-应力方程进行坐标变换；

三、设置复频移-近最佳匹配层中的复频移拉伸函数参数；

四、将复频移动拉伸函数移入到被求导波场函数上；

五、对复频移拉伸函数进行变量代换重新表示为：并带入到拉伸后的频率域弹性波一阶速度-应力方程中；

六、引入辅助波场变量对拉伸后的频率域弹性波一阶速度-应力方程进行变量代换；

七、将拉伸后的频率域弹性波一阶速度-应力方程及辅助方程反变换到时间域，获得时间域最佳匹配层吸收边界条件的弹性波一阶速度-应力方程及其辅助常微分方程；

步骤一中时间域的弹性波一阶速度-应力方程为二维的。

步骤二中，坐标变换是利用复频移拉伸函数实现。

步骤三是通过以下步骤完成：

(1)设定参数kθ的值；

(2)按照多轴最佳匹配层同时在多个方向上引入了衰减因子dθ＝dθ(x，z)，表示为：

(3)求得复频移拉伸函数参数αθ＝παmax(1-δθ/lpml)

其中，步骤(1)中参数kθ为1。

其中步骤五中，变量

其中，步骤四中，复频移拉伸函数为常数。

步骤六后将辅助变量化为一阶常微分方程的形式：

还包括步骤八：利用间断galerkin有限元地震数值模拟技术求解步骤七中的时间域最佳匹配层吸收边界条件弹性波一阶速度-应力方程及其辅助常微分方程，获得长时间模拟保持稳定的模拟结果。

步骤四是利用近最佳匹配层的原则实现的。

本发明具有以下优点：针对在间断galerkin有限元地震波数值模拟算法中使用基于复拉伸坐标变换的近最佳匹配层吸收边界条件(cfs-npml)时存在的长时间模拟应用不稳定的现象，通过借鉴多轴最佳匹配层(mpml)的思想，在cfs-npml中多个方向上同时引入了衰减因子，使在间断galerkin有限元地震波数值模拟算法中使用cfs-npml能保持长时间模拟稳定，且能保持原cfs-npml良好的边界吸收效果，特别是对在复杂地表地震波传播模拟中面波等震相的吸收，为复杂地区地震波高精度有限元数值模拟提供提供稳定良好的吸收边界条件。

附图说明

在下文中将基于实施例并参考附图来对本发明进行更详细的描述。其中：

图1是pml区域参数分布示意图。

图2是采用基于复拉伸坐标变换的近最佳匹配层吸收边界条件(cfs-npml)的波场快照结果图。

图3是采用稳定的基于复拉伸坐标变换的近最佳匹配层吸收边界条件(cfs-npml)的波场快照结果图。

图4自由表面均匀模型示意图。

具体实施方式

下面将结合附图对本发明作进一步说明。

本发明针对在间断galerkin有限元地震波数值模拟算法中使用基于复拉伸坐标变换的近最佳匹配层吸收边界条件(cfs-npml)时存在的长时间模拟应用不稳定现，通过借鉴多轴最佳匹配层(mpml)的思想，在cfs-npml中多个方向上同时引入了衰减因子，达到了cfs-npml的长时间模拟的稳定化，为复杂地区地震波高精度有限元数值模拟提供提供稳定良好的吸收边界条件。该发明的具体步骤包括：

将(二维)时间域的弹性波一阶速度-应力方程变换到频率域，表示为：其中a、b为物性参数矩阵，为频率域波场，ω为角频率，

利用复频移拉伸函数对频率域弹性波一阶速度-应力方程进行坐标变换：

按照多轴最佳匹配层(mpml)设置pml参数的方法设置cfs-npml中的复频移拉伸函数参数：①参数kθ的值一般取为1；②按照mpml同时在多个方向上引入了衰减因子dθ＝dθ(x，z)，表示为

其中，p一般为常系数(见附图1)，为常规pml方法中x和z方向上的衰减系数其中，δθ为单元中心到pml边界的深度，lpml为pml的厚度，③最后，参数αθ＝παmax(1-δθ/lpml)，其中，αmax＝f0为子波的主频，r为理论反射系数，通常取10^-4～10^-6。

利用近最佳匹配层的思想，假定复频移拉伸函数为常数，将拉伸函数移入到被求导波场函数上：

对复频移拉伸函数进行变量代换，重新表示为其中变量并带入到拉伸后的频率域弹性波一阶速度-应力方程中：

引入辅助波场变量其中对拉伸后的频率域弹性波一阶速度-应力方程进行变量代换：

将辅助变量化为一阶常微分方程的形式：

将拉伸后的频率域弹性波一阶速度-应力方程及辅助方程反变换到时间域，获得时间域cfs-npml弹性波一阶速度-应力方程及其辅助常微分方程：

利用间断galerkin有限元地震数值模拟技术求解步骤(7)中的时间域cfs-npml弹性波一阶速度-应力方程及其辅助常微分方程，获得长时间模拟保持稳定的模拟结果。

实施例：

附图1为自由表面均匀模型，其中模型计算区域大小为4000m*4000m，pml区域厚度为400m，采用下面的吸收边界条件结合间断galerkin有限元法(dgfem)进行地震波数值模拟，采用空间3阶模拟，在自由表面激发爆炸震源。；附图2为采用基于复拉伸坐标变换的近最佳匹配层吸收边界条件(cfs-npml)的波场快照结果，左为0.4s时的结果，右为2s时的结果。从图中可以看，模拟计算到2s时已出现明显不稳定的现象，这种不稳定主要来至吸收边界层内部；附图3采用稳定的基于复拉伸坐标变换的近最佳匹配层吸收边界条件(cfs-npml)的波场快照结果，左图为0.4s时的结果，右图为2s时的结果。从图中可以看到2s时外传波场已经基本被吸收干净，达到了稳定的模拟。