基于二阶偏微分波动方程的pml吸收边界条件的方法
【专利摘要】本发明提供了一种基于二阶偏微分波动方程的PML吸收边界条件的方法，包括以下步骤：加载某一采样时刻的震源的纵波波场；计算所述采样时刻的纵波波场在三维空间某一方向的二阶偏导数，并根据所述方向的二阶偏导数对所述方向的完全匹配层区域内的纵波波场的时间一阶导数波场进行吸收衰减；重复上一步骤，分别对所述采样时刻的纵波波场在三维空间另外两个方向的完全匹配层区域内的纵波波场的时间一阶导数波场进行吸收衰减；重复以上步骤，直至完成所述纵波波场所有采样时刻的处理。采用本发明的完全匹配层吸收边界条件的实现复杂度低、计算存储量小。
【专利说明】基于二阶偏微分波动方程的PML吸收边界条件的方法
【技术领域】
[0001]本发明涉及地球物理勘探领域中一种波动方程数值模拟技术，尤其涉及一种基于二阶偏微分波动方程的PML (perfect matched layer,完全匹配层)吸收边界条件的方法。
【背景技术】
[0002]地震波数值模拟技术在地球物理勘探领域得到广泛的应用，常用来设计地震采集观测系统、研究复杂储层的地震波响应、实现地震波逆时偏移成像和全波形反演速度建模技术等。
[0003]在波动方程数值模拟技术中，一项重要的研究内容就是吸收边界条件。在现实世界中，速度介质模型的边界可能为无穷大；采用计算机进行数值模拟实现时，只能对无限大的模型进行截断，这样就会导致地震波传到截断边界时产生不渴望的虚假反射。因此，必须对这种虚假的边界反射进行吸收或消除，一般使用人为地引入吸收边界条件的方法来解决。
[0004]目前应用于波动方程的吸收边界条件主要有三类。第一类是海绵吸收边界条件，在一定层厚的介质中对波场逐渐进行衰减，如阻尼吸收边界条件；第二类是基于波动方程的旁轴近似吸收边界条件，本质上是一种单程波，不足之处是不能很好处理入射角度较大的反射问题；第三类就是基于波动方程的PML吸收边界条件，由于不同完全匹配层之间的复波阻抗相同，波在完全匹配层之间传播时不产生任何反射只产生衰减。到目前为止，完全匹配层吸收边界条件是吸收效果最理想的吸收边界条件，在数值模拟精度较高的场合，一般主要采用完全匹配层吸收边界条件。其中，PML是通过在FDTD (Finite-DifferenceTime-Domain,时域有限差分)区域截断边界设置一种特殊的介质层，该层介质的波阻抗与相邻介质的波阻抗完全匹配，入射波将无反射地穿过分界面而进入PML。
[0005]但是，现有的完全匹配层吸收边界条件只适用于一阶偏微分波动方程数值模拟，主要采用的为分裂式完全匹配层吸收边界条件。在计算过程中，为了节省计算存储，现有技术常常在完全匹配层区域对波场进行分裂，对于三维空间问题，其需要考虑的区域包括6个面部区域、12个棱边区域和8个角点区域，导致实现起来非常复杂。
【发明内容】
[0006]本发明的目的在于提供一种基于二阶偏微分波动方程的PML吸收边界条件的方法，以减小在波动方程数值模拟技术中采用完全匹配层吸收边界条件的实现复杂度。
[0007]为达到上述目的，本发明提供了一种基于二阶偏微分波动方程的PML吸收边界条件的方法，包括以下步骤:
[0008]加载某一米样时刻的震源的纵波波场；
[0009]计算所述采样时刻的纵波波场在三维空间某一方向的二阶偏导数，并根据所述方向的二阶偏导数对所述方向的完全匹配层区域内的纵波波场的时间一阶导数波场进行吸收衰减；[0010]重复上一步骤，分别对所述采样时刻的纵波波场在三维空间另外两个方向的完全匹配层区域内的纵波波场的时间一阶导数波场进行吸收衰减；
[0011]重复以上步骤，直至完成所述纵波波场所有采样时刻的处理。
[0012]本发明的方法，在所述在整个三维数值模拟计算区域对所述纵波波场进行更新之后，重复以上步骤，直至完成所述纵波波场所有采样时刻的处理之前，还包括:
[0013]在整个三维数值模拟计算区域对所述纵波波场进行更新。
[0014]本发明的方法，在所述加载某一采样时刻的震源的纵波波场之前，还包括:
[0015]预先将纵波波场的三维速度模型在三维空间上进行网格离散化，并设置波动方程数值模拟所需参数；
[0016]预先将经过网格离散化处理后的三维速度模型的整个三维空间数值模拟区域划分为主区域和完全匹配层区域，并在所述完全匹配层区域的每个方向各引入两个完全匹配层辅助波场，所述主区域为纵波波场传播时能量不被吸收衰减的区域；
[0017]预先对所述完全匹配层区域进行参数设置。
[0018]本发明的方法，所述波动方程数值模拟所需参数包括采样网格数、三维空间方向的采样步长、三维空间方向的采样点数、时间采样点数、时间采样步长、地震子波类型和主频；所述主区域的波动方程为:
【权利要求】
1.一种基于二阶偏微分波动方程的PML吸收边界条件的方法，其特征在于，包括以下步骤: 加载某一米样时刻的震源的纵波波场； 计算所述采样时刻的纵波波场在三维空间某一方向的二阶偏导数，并根据所述方向的二阶偏导数对所述方向的完全匹配层区域内的纵波波场的时间一阶导数波场进行吸收衰减； 重复上一步骤，分别对所述采样时刻的纵波波场在三维空间另外两个方向的完全匹配层区域内的纵波波场的时间一阶导数波场进行吸收衰减； 重复以上步骤，直至完成所述纵波波场所有采样时刻的处理。
2.根据权利要求1所述的方法，其特征在于，在所述在整个三维数值模拟计算区域对所述纵波波场进行更新之后，重复以上步骤，直至完成所述纵波波场所有采样时刻的处理之前，还包括: 在整个三维数值模拟计算区域对所述纵波波场进行更新。
3.根据权利要求2所述的方法，其特征在于，在所述加载某一采样时刻的震源的纵波波场之前，还包括: 预先将纵波波场的三维速度模型在三维空间上进行网格离散化，并设置波动方程数值模拟所需参数； 预先将经过网格离散化处理后的三维速度模型的整个三维空间数值模拟区域划分为主区域和完全匹配层区域，并在所述完全匹配层区域的每个方向各引入两个完全匹配层辅助波场，所述主区域为纵波波场传播时能量不被吸收衰减的区域； 预先对所述完全匹配层区域进行参数设置。
4.根据权利要求3所述的方法，其特征在于，所述波动方程数值模拟所需参数包括采样网格数、三维空间方向的采样步长、三维空间方向的采样点数、时间采样点数、时间采样步长、地震子波类型和主频；所述主区域的波动方程为:
5.根据权利要求4所述的方法，其特征在于， 所述X方向完全匹配层区域的波动方程为:
6.根据权利要求3所述的方法，其特征在于，所述预先对完全匹配层区域进行参数设置包括: 设置完全匹配层吸收边界条件的层数、最大吸收衰减值和衰减函数。
7.根据权利要求5所述的方法，其特征在于，所述计算所述采样时刻的纵波波场在三维空间某一方向的二阶偏导数，具体为: 采用有限差分法或伪谱法计算X方向、y方向或z方向的二阶偏导数。
8.根据权利要求7所述的方法，其特征在于，所述根据所述方向的二阶偏导数，对所述方向的完全匹配层区域内的纵波波场的时间一阶导数波场进行吸收衰减，具体为: 根据所述方向的二阶偏导数，在所述整个三维空间数值模拟区域更新所述方向的纵波波场的时间一阶导数波场； 根据所述方向的二阶偏导数，在所述方向的完全匹配层区域更新所述方向的完全匹配层辅助波场； 结合更新后的所述方向的完全匹配层辅助波场，对更新后的所述方向的完全匹配层区域内的纵波波场的时间一阶导数波场进行吸收衰减。
9.根据权利要求8所述的方法，其特征在于，所述根据所述方向的二阶偏导数，在所述整个三维空间数值模拟区域更新所述方向的纵波波场的时间一阶导数波场，具体为: 根据X方向的二阶偏导数g以及公式
10.根据权利要求8所述的方法，其特征在于，所述根据所述方向的二阶偏导数，在所述方向的完全匹配层区域更新所述方向的完全匹配层辅助波场，具体为:根据X方向的二阶偏导数g对所述Ex和所述Hx进行更新计算，具体通过公式:
11.根据权利要求8所述的方法，其特征在于，所述结合更新后的所述方向的完全匹配层辅助波场，对更新后的所述方向的完全匹配层区域内的纵波波场的时间一阶导数波场进行吸收衰减，具体为:
12.根据权利要求8所述的方法，其特征在于，在整个三维数值模拟计算区域对所述纵波波场进行更新，具体为:在整个三维数值模拟计算区域根据公式P — P+ Δ tQ对所述纵波波场进行更新。
13.根据权利要求10所述的方法，其特征在于，所述
14.根据权利要求6所述的方法，其特征在于，所述的完全匹配层的层数的取值范围为4~10层。
15.根据权利要求6所述的方法，其特征在于，所述的完全匹配层最大吸收衰减值的取值范围为300~2000。
16.根据权利要求6所述的方法，其特征在于，所述的完全匹配层的衰减函数为
【文档编号】G01V1/28GK103616721SQ201310606742
【公开日】2014年3月5日 申请日期:2013年11月25日 优先权日:2013年11月25日
【发明者】秦臻, 宋建勇, 马晓宇, 宋雪娟 申请人:中国石油天然气股份有限公司