一种实现地震波场数值模拟方法技术

本发明专利技术涉及地球物理勘探技术领域，特别是关于地震波场数值模拟技术领域，具体的讲是一种地震波场数值模拟的方法，包括在波数域使用低阶有理分式有限差分算子逼近超高阶显式有限差分算子；计算所述低阶有理分式有限差分算子的系数；将有理分式有限差分算子进行分解，利用傅里叶变换将低阶有理分式有限差分算子从波数域变换到空间域；在空间域使用有理分式有限差分法分别计算不同方向上的空间二阶偏导数项；在时间方向上进行迭代实现地震波场的数值模拟。通过本发明专利技术实施例的方法，使用低阶有理分式有限差分法来实现地震波数值模拟，提高了地震波场数值模拟的精度和计算效率。

【技术实现步骤摘要】
一种实现地震波场数值模拟方法
本专利技术涉及地球物理勘探
，特别是关于地震波场数值模拟
，具体的讲是一种实现地震波场数值模拟方法。
技术介绍
地震波数值模拟或者地震波波动方程正演是地球物理勘探领域的重要课题，在研究地下复杂介质的地震响应、复杂介质的高精度地震成像（例如逆时偏移）、高精度的地下介质速度建模（例如全波形反演）等领域有着广泛的应用。现有技术中进行地震波数值模拟时主要采用如下几种方法：1、有限元法基于变分原理和网格插值的有限元法具有灵活的不规则网格剖分方案，能从几何形状上适应不规则地下介质或起伏的地表，但是算法复杂、计算存储量大、计算精度较低和计算效率低下，不适合于大规模的地震波场数值模拟。2、伪谱法利用傅里叶正变换将空间偏导数的计算转换到波数域进行，再通过傅里叶反变换得到空间域的空间偏导数值。理论上，伪谱法计算空间偏导数具有最高的数值精度。一般来说，伪谱法的计算效率要低于有限差分法。3、谱元法将有限元法和伪谱法结合就产生了谱元法，它充分利用了有限元法的灵活网格剖分特点和伪谱法数值模拟精度高的优点，因此在求解波动方程时，谱元法所用的网格节点少、数值精度高。但是，谱元法涉及到矩阵求逆，计算效率仍然非常低下，比较适合于小规模和中等规模的波动方程数值模拟。4、有限差分法有限差分法是波动方程数值模拟中实际应用中使用最多的方法，算法相对比较简单、与前面三种方法相比计算效率较高，适合于大规模复杂介质的地震波场数值模拟。常见的有限差分法一般分为显式有限差分法和隐式有限差分法，隐式有限差分法一般在频率域进行，常常用于波动方程反演（如全波形反演），由于只需要少数几个频率采样点，因此计算效率相对比较高，但是对于三维问题往往带来巨大的存储量和计算量；显式差分法一般在时间-空间域进行，针对所有频率采样点同时计算，常常用于波动方程正演、波动方程成像和波动方程反演。对于高阶显式有限差分法，随着差分阶数的增加，计算精度也随之增加。理论的空间二阶偏导数对应的差分算子响应曲线为一抛物线，高阶显式差分算子响应曲线对抛物线的逼近程度越高，实现数值模拟时的计算精度就越高。附图3为2M阶差分精度的高阶显式差分算子在波数域的响应，随着M值的增大，对应的曲线越来越逼近理论的抛物线。但是，随着差分阶数的增加，高阶显式有限差分的计算效率也随着降低。高阶显式有限差分法的计算效率与差分阶数呈反比例关系。以x坐标轴方向的空间偏导数为例，2M阶差分精度的高阶显式有限差分法其公式为：从公式可以看出，对于N个采样点的空间二阶偏导数计算，2M阶差分精度的高阶显式有限差分法涉及的乘除法运算量约为(M+1)N。乘除法运算量越大，计算效率就越低，实际应用中不能将差分精度设置得过高。兼顾计算精度和计算效率，M值一般取1到6之间。
技术实现思路
为了解决现有技术的地震波场数值模拟中数值模拟时间和数值模拟精度不能兼顾的问题，提出了一种实现地震波场数值模拟方法，可以使用阶数较低的有理分式有限差分算子去逼近阶数很高的超高阶显式有限差分算子，使其以较短的数值模拟时间而具备较高的数值模拟精度。本专利技术实施例提供了一种实现地震波场数值模拟方法，包括，步骤101，在波数域使用低阶有理分式有限差分算子逼近超高阶显式有限差分算子；步骤102，计算所述低阶有理分式有限差分算子的系数；步骤103，将有理分式有限差分算子进行分解，利用傅里叶变换将低阶有理分式有限差分算子从波数域变换到空间域；步骤104，在空间域使用有理分式有限差分法分别计算不同方向上的空间二阶偏导数项；步骤105，在时间方向上进行迭代实现地震波场的数值模拟。根据本专利技术实施例所述的一种实现地震波场数值模拟方法的一个进一步的方面，使用如下公式逼近超高阶显式有限差分算子：其中，方程左端项为波数域的具备2M阶差分精度的超高阶显式差分算子，右端项的分子项为波数域的2L阶差分精度显式差分算子的求和，在空间域通过显式差分法计算；方程右端项的分母项为波数域的二阶隐式差分算子的连乘，在空间域通过隐式差分法求解L个三对角线性方程组，cl和dl为待求的低阶有理分式有限差分的系数，kh为波数域的自变量，Δh为空间采样间隔，bm为超高阶显式差分算子的系数，其中L、M为正整数。根据本专利技术实施例所述的一种实现地震波场数值模拟方法的再一个进一步的方面，M的取值为24到28之间，L的取值为1到5之间。根据本专利技术实施例所述的一种实现地震波场数值模拟方法的另一个进一步的方面，所述bm通过如下公式计算得到：根据本专利技术实施例所述的一种实现地震波场数值模拟方法的另一个进一步的方面，通过最小二乘法、模拟退火，神经网络或遗传算法计算来计算所述低阶有理分式有限差分算子的系数cl和dl。根据本专利技术实施例所述的一种实现地震波场数值模拟方法的另一个进一步的方面，使用所述最小二乘法计算所述低阶有理分式有限差分算子的系数cl和dl，具体包括，步骤201，将有理分式有限差分算子的分母项转化为等价的L阶多项式求和，所得到的计算公式为：其中fl和el为引入的临时变量，满足关系：步骤202，整理有理分式、将非线性的逼近拟合问题转化为线性逼近拟合问题，通过如下公式进行；步骤203，通过最小二乘法计算线性方程组的解fl和el，所使用的线性方程组为：其中为ψi,l和ψi,L+l为基函数，GM,i为自由项，Nh为空间采样点数；自由项GM,i定义为...

【技术保护点】
一种实现地震波场数值模拟方法，其特征在于包括，步骤101，在波数域使用低阶有理分式有限差分算子逼近超高阶显式有限差分算子；步骤102，计算所述低阶有理分式有限差分算子的系数；步骤103，将有理分式有限差分算子进行分解，利用傅里叶变换将低阶有理分式有限差分算子从波数域变换到空间域；步骤104，在空间域使用有理分式有限差分法分别计算不同方向上的空间二阶偏导数项；步骤105，在时间方向上进行迭代实现地震波场的数值模拟。

【技术特征摘要】
1.一种实现地震波场数值模拟方法，其特征在于包括，步骤101，在波数域使用低阶有理分式有限差分算子逼近超高阶显式有限差分算子；其中，使用如下公式逼近超高阶显式有限差分算子：其中，方程左端项为波数域的具备2M阶差分精度的超高阶显式差分算子，右端项的分子项为波数域的2L阶差分精度显式差分算子的求和，在空间域通过显式差分法计算；方程右端项的分母项为波数域的二阶隐式差分算子的连乘，在空间域通过隐式差分法求解L个三对角线性方程组，cl和dl为待求的低阶有理分式有限差分的系数，kh为波数域的自变量，Δh为空间采样间隔，bm为超高阶显式差分算子的系数，其中L、M为正整数；步骤102，计算所述低阶有理分式有限差分算子的系数；步骤103，将有理分式有限差分算子进行分解，利用傅里叶变换将低阶有理分式有限差分算子从波数域变换到空间域；步骤104，在空间域使用有理分式有限差分法分别计算不同方向上的空间二阶偏导数项；步骤105，在时间方向上进行迭代实现地震波场的数值模拟。2.根据权利要求1所述的一种实现地震波场数值模拟方法，其特征在于，M的取值为24到28之间，L的取值为1到5之间。3.根据权利要求1所述的一种实现地震波场数值模拟方法，其特征在于，所述bm通过如下公式计算得到：4.根据权利要求1所述的一种实现地震波场数值模拟方法，其特征在于，通过最小二乘法、模拟退火、神经网络或遗传算法来计算所述低阶有理分式有限差分算子的系数cl和dl。5.根据权利要求4所述的一种实现地震波场数值模拟方法，其特征在于，使用所述最小二乘法计算所述低阶有理分式有限差分算子的系数cl和dl，具体包括，步骤201，将有理分式有限差分算子的分母项转化为等价的L阶多项式求和，所得到的计算公式为：其中fl和el为引入的临时变量，满足关系：步骤202，整理有理分式，将非线性的逼近拟合问题转化为线性逼近拟合问题，通过如下公式进行；步骤203，通过最小二乘法计算线性方程组的解fl和el，所使用的线性方程组为：其中，ψi,l和ψi,L+l为基函数，GM,i为自由项，Nh为空间采样点数；自由项GM,i定义为基函数ψi,l和ψi,L+l定义为ψi,l＝1-cos(lθi)ψi,L+l＝-GM,icosl(θi),(i＝1,2,…,Nh；l＝1,2,…,L)其中θi为离散化后的角度采样，它与波数kh的关系为：步骤204，利用上述的fl和el，结合求根公式法由el计算dl，然后由dl和fl计算cl；求根公式法是计算关于变量z的L阶多项式的根，对应着因式分解，由el计算dl所使用的公式为：由...

【专利技术属性】
技术研发人员：秦臻，宋建勇，马晓宇，宋雪娟，
申请(专利权)人：中国石油天然气股份有限公司，
类型：发明
国别省市：北京;11