技术领域
[0001] 本
发明涉及地球物理勘探技术领域，特别是关于
地震波场数值模拟技术领域，具体的讲是一种实现地震波场数值模拟方法。
背景技术
[0002] 地震
波数值模拟或者地震波
波动方程正演是地球物理勘探领域的重要课题，在研究地下复杂介质的地震响应、复杂介质的高
精度地震成像（例如逆时偏移）、高精度的地下介质速度建模（例如全
波形反演）等领域有着广泛的应用。
现有技术中进行地震波数值模拟时主要采用如下几种方法：
[0003] 1、有限元法
[0004] 基于变分原理和网格插值的有限元法具有灵活的不规则网格剖分方案，能从几何形状上适应不规则地下介质或起伏的地表，但是
算法复杂、计算存储量大、计算精度较低和计算效率低下，不适合于大规模的地震波场数值模拟。
[0005] 2、伪谱法
[0006] 利用傅里叶正变换将空间偏导数的计算转换到波数域进行，再通过傅里叶反变换得到空间域的空间偏导数值。理论上，伪谱法计算空间偏导数具有最高的数值精度。一般来说，伪谱法的计算效率要低于有限差分法。
[0007] 3、谱元法
[0008] 将有限元法和伪谱法结合就产生了谱元法，它充分利用了有限元法的灵活网格剖分特点和伪谱法数值模拟精度高的优点，因此在求解波动方程时，谱元法所用的网格
节点少、数值精度高。但是，谱元法涉及到矩阵求逆，计算效率仍然非常低下，比较适合于小规模和中等规模的波动方程数值模拟。
[0009] 4、有限差分法
[0010] 有限差分法是波动方程数值模拟中实际应用中使用最多的方法，算法相对比较简单、与前面三种方法相比计算效率较高，适合于大规模复杂介质的地震波场数值模拟。常见的有限差分法一般分为显式有限差分法和隐式有限差分法，隐式有限差分法一般在
频率域进行，常常用于波动方程反演（如全波形反演），由于只需要少数几个频率
采样点，因此计算效率相对比较高，但是对于三维问题往往带来巨大的存储量和计算量；显式差分法一般在时间-空间域进行，针对所有频率采样点同时计算，常常用于波动方程正演、波动方程成像和波动方程反演。
[0011] 对于高阶显式有限差分法，随着差分阶数的增加，计算精度也随之增加。理论的空间二阶偏导数对应的差分算子响应曲线为一抛物线，高阶显式差分算子响应曲线对抛物线的逼近程度越高，实现数值模拟时的计算精度就越高。
附图3为2M阶差分精度的高阶显式差分算子在波数域的响应，随着M值的增大，对应的曲线越来越逼近理论的抛物线。
[0012] 但是，随着差分阶数的增加，高阶显式有限差分的计算效率也随着降低。高阶显式有限差分法的计算效率与差分阶数呈反比例关系。以x坐标轴方向的空间偏导数为例，2M阶差分精度的高阶显式有限差分法其公式为：
[0013]
[0014] 从公式可以看出，对于N个采样点的空间二阶偏导数计算，2M阶差分精度的高阶显式有限差分法涉及的乘除法运算量约为(M+1)N。乘除法运算量越大，计算效率就越低，实际应用中不能将差分精度设置得过高。兼顾计算精度和计算效率，M值一般取1到6之间。
发明内容
[0015] 为了解决现有技术的地震波场数值模拟中数值模拟时间和数值模拟精度不能兼顾的问题，提出了一种实现地震波场数值模拟方法，可以使用阶数较低的有理分式有限差分算子去逼近阶数很高的超高阶显式有限差分算子，使其以较短的数值模拟时间而具备较高的数值模拟精度。
[0016] 本发明
实施例提供了一种实现地震波场数值模拟方法，包括，
[0017] 步骤101，在波数域使用低阶有理分式有限差分算子逼近超高阶显式有限差分算子；
[0018] 步骤102，计算所述低阶有理分式有限差分算子的系数；
[0019] 步骤103，将有理分式有限差分算子进行分解，利用傅里叶变换将低阶有理分式有限差分算子从波数域变换到空间域；
[0020] 步骤104，在空间域使用有理分式有限差分法分别计算不同方向上的空间二阶偏导数项；
[0021] 步骤105，在时间方向上进行
迭代实现地震波场的数值模拟。
[0022] 根据本发明实施例所述的一种实现地震波场数值模拟方法的一个进一步的方面，使用如下公式逼近超高阶显式有限差分算子：
[0023]
[0024] 其中，方程左端项为波数域的具备2M阶差分精度的超高阶显式差分算子，右端项的分子项为波数域的2L阶差分精度显式差分算子的求和，在空间域通过显式差分法计算；方程右端项的分母项为波数域的二阶隐式差分算子的连乘，在空间域通过隐式差分法求解L个三对
角线性方程组，cl和dl为待求的低阶有理分式有限差分的系数，kh为波数域的自变量，Δh为空间采样间隔，bm为超高阶显式差分算子的系数，其中L、M为正整数。
[0025] 根据本发明实施例所述的一种实现地震波场数值模拟方法的再一个进一步的方4 8
面，M的取值为2 到2 之间，L的取值为1到5之间。
[0026] 根据本发明实施例所述的一种实现地震波场数值模拟方法的另一个进一步的方面，所述bm通过如下公式计算得到：
[0027]
[0028]
[0029] 根据本发明实施例所述的一种实现地震波场数值模拟方法的另一个进一步的方面，通过最小二乘法、模拟
退火，神经网络或
遗传算法计算来计算所述低阶有理分式有限差分算子的系数cl和dl。
[0030] 根据本发明实施例所述的一种实现地震波场数值模拟方法的另一个进一步的方面，使用所述最小二乘法计算所述低阶有理分式有限差分算子的系数cl和dl，具体包括，[0031] 步骤201，将有理分式有限差分算子的分母项转化为等价的L阶多项式求和，所得到的计算公式为：
[0032]
[0033] 其中fl和el为引入的临时变量，满足关系：
[0034]
[0035] 步骤202，
整理有理分式、将非线性的逼近拟合问题转化为线性逼近拟合问题，通过如下公式进行；
[0036]
[0037] 步骤203，通过最小二乘法计算线性方程组的解fl和el，所使用的线性方程组为：
[0038]
[0039] 其中为ψi,l和ψi,L+l为基函数，GM,i为自由项，Nh为空间采样点数；
[0040] 自由项GM,i定义为
[0041]
[0042] 基函数ψi,l和ψi,L+l定义为
[0043] ψi,l＝1-cos(lθi)
[0044] ψi,L+l＝-GM,icosl(θi),(i＝1,2,…,Nh;l＝1,2,…,L)
[0045] 其中θi为离散化后的角度采样，它与波数kh的关系为：
[0046]
[0047] 步骤204，利用上述的fl和el，结合求根公式法由el计算dl，然后由dl和fl计算cl；
[0048] 求根公式法是计算关于变量z的L阶多项式的根，对应着因式分解，由el计算dl所使用的公式为：
[0049]
[0050] 由dl和fl计算cl的公式为：
[0051] 根据本发明实施例所述的一种实现地震波场数值模拟方法的另一个进一步的方面，在将有理分式有限差分算子进行分解，利用傅里叶变换将低阶有理分式有限差分算子从波数域变换到空间域中进一步包括：
[0052] 步骤301，使用有理分式有限差分计算空间二阶偏导数，在波数域通过下面公式来实现：
[0053]
[0054] 其中 表示波数域的波动方程波场变量；
[0055] 步骤302，引入临时变量 在波数域将有理分式有限差分算子进行分解，得到一个2L阶显式差分方程、L个二阶隐式差分方程和二阶偏导数计算方程；
[0056] 其中，2L阶显式差分方程为：
[0057]
[0058] 第l个二阶隐式差分方程为：
[0059]
[0060] 二阶偏导数计算方程为：
[0061]
[0062] 和 为引入的L+1个临时变量；
[0063] 步骤303，通过傅里叶变换将上述步骤302中的所有方程从波数域变换到空间域，所得到的2L阶显式差分方程为：
[0064]
[0065] 所使用的第l个二阶隐式差分方程为：
[0066] ql(h-Δh)+dlql(h)+ql(h+Δh)＝ql-1(h),(l＝1,2,…,L)
[0067] 以及二阶偏导数计算方程：
[0068]
[0069] 其中，其中p(h)表示空间域的波动方程波场变量，q0(h)和ql(h)(l＝1,2,…,L)为引入的L+1个临时变量，h表示空间域的空间变量x、y或z；Δh为空间域采样间隔Δx、Δy或Δz；系数c0通过公式 计算。
[0070] 根据本发明实施例所述的一种实现地震波场数值模拟方法的另一个进一步的方面，所述在空间域使用有理分式有限差分法分别计算不同方向上的空间二阶偏导数项进一步包括：
[0071] （说明：所述的计算空间二阶偏导数，是指直接地或间接地由p(h)计算[0072] 步骤401，结合已知的空间域的波动方程波场变量p(h)，使用显式有限差分法计算步骤103的步骤303中的第一个方程得到q0(h)；
[0073] 步骤402，结合上述步骤401计算的q0(h)，依次由ql-1(h),(l＝1,2,…,L)计算步骤103的步骤303中的L个二阶隐式方程得到ql(h),(l＝1,2,…,L)；
[0074] 步骤403，结合上述步骤402最终计算的qL(h)，通过步骤103的步骤303中的第三个公式，得到二阶空间偏导数结果
[0075] 步骤404，通过上述步骤401～403，计算x方向的空间偏导数
[0076] 步骤405，通过步骤401～步骤403，计算y方向的空间偏导数
[0077] 步骤406，通过步骤401～步骤403，计算z方向的空间偏导数
[0078] 根据本发明实施例所述的一种实现地震波场数值模拟方法的另一个进一步的方面，在所述步骤402中通过追赶法、高斯消元法或矩阵求逆法来求解L个二阶隐式方程。
[0079] 根据本发明实施例所述的一种实现地震波场数值模拟方法的另一个进一步的方面，在所述步骤402中通过改进追赶法求解L个二阶隐式方程：
[0080] 步骤601，预先计算求解L个三对角线性方程组所需的L个常向量αl，通过以下递推公式来实现：
[0081]
[0082]
[0083] 其中dl为前面所求的隐式差分系数；αl表示采样长度为Nh的向量，Nh取空间采样网格数Nx、Ny和Nz中的最大值，Nx、Ny和Nz分别为x、y和z坐标轴方向的空间网格数；
[0084] 步骤602，实现改进追赶法的追过程，通过以下递推公式进行：
[0085] w(1)＝αl(1)ql-1(1)
[0086] w(i)＝αl(i)[ql-1(i)-w(i-1)],(i=1,2,…,Nh)
[0087] 其中w为实现追赶法引入的中间变量，ql-1为前面所述第l-1阶隐式差分方程中的变量，αl为预先计算的常向量，上述公式由ql-1计算w；
[0088] 步骤603，实现改进追赶法的赶过程，通过以下递推公式进行：
[0089] ql(N)＝w(N)
[0090] ql(i)＝w(i)-αl(i)ql(i+1),(i＝Nh-1,Nh-2,…,1)
[0091] 其中w为实现追赶法引入的中间变量，ql为前面所述第l-1阶隐式差分方程中的变量，αl为预先计算的常向量，上述公式由w计算ql。
[0092] 根据本发明实施例所述的一种实现地震波场数值模拟方法的另一个进一步的方面，在时间方向上进行迭代实现地震波场的数值模拟进一步包括：
[0093]n n-1 n+1
[0094] 其中p 表示t＝nΔt时刻传播的波场、p 表示前一时刻传播的波场，p 表示后一时刻传播的波场，Δt为时间采样间隔，Vp为弹性介质的纵
波速度，时间方向上的每一次迭代计算需通过所述步骤104中的所述步骤404～406来完成，共需要迭代Nt次。
[0095] 通过本发明实施例的方法，使用低阶有理分式有限差分法来实现地震波数值模拟，提高了地震波场数值模拟的精度和计算效率。
附图说明
[0096] 结合以下附图阅读对实施例的详细描述，本发明的上述特征和优点，以及额外的特征和优点，将会更加清楚。
[0097] 图1为本发明实施例提供的一种实现地震波场数值模拟方法的
流程图；
[0098] 图2为本发明实施例中使用有理分式有限差分法计算二阶空间偏导数的流程图；
[0099] 图3为现有技术中2M阶差分精度的高阶显式有限差分算子在波数域的响应曲线比较示意图；
[0100] 图4a为本发明实施例中四阶有理分式有限差分算子去拟合64阶差分精度的显式高阶差分算子中两个算子对应的响应曲线示意图；
[0101] 图4b为本发明实施例中四阶有理分式有限差分算子去拟合64阶差分精度的显式高阶差分算子中拟合相对误差曲线示意图；
[0102] 图5a-图5d为本发明实施例中不同有限差分方法所模拟的波场快照的比较示意图；
[0103] 图6a-图6d为本发明实施例中不同有限差分方法所模拟的单道记录的比较示意图；
[0104] 图7为本发明实施例中一个复杂的速度介质模型示意图；
[0105] 图8为现有技术中使用12阶高阶显式有限差分法对复杂速度介质模型所模拟的单炮记录示意图；
[0106] 图9为本发明实施例中使用四阶有理分式有限差分差分法对复杂速度介质模型所模拟的单炮记录示意图。
具体实施方式
[0107] 下面的描述可以使任何本领域技术人员利用本发明。具体实施例和应用中所提供的描述信息仅为示例。这里所描述的实施例的各种延伸和组合对于本领域的技术人员是显而易见的，在不脱离本发明的实质和范围的情况下，本发明定义的一般原则可以应用到其他实施例和应用中。因此，本发明不只限于所示的实施例，本发明涵盖与本文所示原理和特征相一致的最大范围。
[0108] 如图1所示为本发明实施例一种实现地震波场数值模拟方法的流程图。
[0109] 包括步骤101，在波数域使用低阶有理分式有限差分算子逼近超高阶显式有限差分算子。
[0110]
[0111] 其中，右端项的分子项为波数域的2L阶差分精度显式差分算子的求和，在空间域通过显式差分法计算；方程右端项的分母项为波数域的二阶隐式差分算子的连乘，在空间域通过隐式差分法求解L个三对角线性方程组，cl和dl为待求的低阶有理分式有限差分的系数，kh为波数域的自变量，Δh为空间采样间隔，L的取值一般在1到5之间，M的取值一般在8到64之间；bm为超高阶显式差分算子的系数。
[0112] 作为本发明的一个实施例，方程左端项为波数域的具备2M阶差分精度的超高阶显式差分算子，M的取值一般在24到28之间。方程的右端项为2L阶有理分式有限差分算子，L的取值一般在1到5之间。
[0113] 作为本发明的一个实施例，所述bm通过如下公式计算得到：
[0114]
[0115]
[0116] 步骤102，计算所述低阶有理分式有限差分算子的系数cl和dl。
[0117] 作为本发明的一个实施例，可以通过最小二乘法来计算所述地接有理分式有限差分算子的系数cl和dl，还可以使用模拟退火，神经网络，遗传算法计算所述低阶有理分式有限差分算子的系数cl和dl。
[0118] 作为本发明的一个实施例，使用所述最小二乘法计算所述低阶有理分式有限差分算子的系数cl和dl，具体包括，
[0119] 201.将有理分式有限差分算子的分母项转化为等价的L阶多项式求和，所得到的计算公式为：
[0120]
[0121] 其中fl和el为引入的临时变量，满足关系：
[0122]
[0123] 202.整理有理分式、将非线性的逼近拟合问题转化为线性逼近拟合问题，通过如下公式进行；
[0124]
[0125] 203.通过最小二乘法计算线性方程组的解fl和el，所使用的线性方程组为：
[0126]
[0127] 其中为ψi,l和ψi,L+l为基函数，GM,i为自由项，Nh为空间采样点数；
[0128] 自由项GM,i定义为
[0129]
[0130] 基函数ψi,l和ψi,L+l定义为
[0131] ψi,l＝1-cos(lθi)
[0132] ψi,L+l＝-GM,icosl(θi),(i＝1,2,…,Nh;l＝1,2,…,L)
[0133] 其中θi为离散化后的角度采样，它与波数kh的关系为：
[0134]
[0135] 204.利用步骤203中计算的fl和el，结合求根公式法由el计算dl，然后由dl和fl计算cl。
[0136] 求根公式法是计算关于变量z的L阶多项式的根，对应着因式分解。由el计算dl所使用的公式为：
[0137]
[0138] 由dl和fl计算cl的公式为：
[0139] 步骤103，将有理分式有限差分算子进行分解，利用傅里叶变换将低阶有理分式有限差分算子从波数域变换到空间域。
[0140] 作为本发明的一个实施例，在本步骤103中具体包括，
[0141] 301.使用有理分式有限差分来计算空间二阶偏导数，在波数域通过下面公式来实现：
[0142]
[0143] 其中 表示波数域的波动方程波场变量；
[0144] 302.引入临时变量 在波数域将有理分式有限差分算子进行分解，得到一个2L阶显式差分方程、L个二阶隐式差分方程和二阶偏导数计算方程；
[0145] 其中，2L阶显式差分方程为：
[0146]
[0147] 第l个二阶隐式差分方程为：
[0148]
[0149] 二阶偏导数计算方程为：
[0150]
[0151] 和 为引入的L+1个临时变量；
[0152] 303.通过傅里叶变换技术，将步骤302中的所有方程从波数域变换到空间域，所得到的2L阶显式差分方程为：
[0153]
[0154] 所使用的第l个二阶隐式差分方程为：
[0155] ql(h-Δh)+dlql(h)+ql(h+Δh)＝ql-1(h),(l＝1,2,…,L)
[0156] 以及二阶偏导数计算方程：
[0157]
[0158] 其中，其中p(h)表示空间域的波动方程波场变量，q0(h)和ql(h)(l＝1,2,…,L)为引入的L+1个临时变量。h表示空间域的空间变量x、y或z；Δh为空间域采样间隔Δx、Δy或Δz；系数c0通过公式 计算。
[0159] 步骤104，在空间域使用有理分式有限差分法来实现步骤103的步骤303的所有公式，由p(h)来计算不同方向上地震波波动方程中的空间二阶偏导数项 和[0160] 如图2所示为本发明实施例使用有理分式有限差分法计算二阶空间偏导数的流程图，具体包括，
[0161] 401.结合已知的空间域的波动方程波场变量p(h)，使用常规的显式有限差分法计算步骤103中步骤303中的第一个方程，得到q0(h)；
[0162] 402.结合上述步骤401中q0(h)，依次求解L个三对角线性方程组，对应于依次计算步骤103步骤303中的L个二阶隐式方程得到ql(h),(l＝1,2,…,L)；所述的二阶隐式方程本质上表现为三对角线性方程组。
[0163] 其中可以通过追赶法、高斯消元法或矩阵求逆法来求解L个三对角线性方程组。
[0164] 其中，步骤104步骤401中的计算输出结果q0(h)作为第一个二阶隐式差分方程的右端输入项，输出结果为q1(h)；第l-1个二阶隐式方程计算输出结果ql-1(h)作为第l个二阶隐式差分方程的右端输入项，输出结果为ql(h)；第L-1个二阶隐式方程计算输出结果qL-1(h)作为第L个二阶隐式差分方程的右端输入项，输出结果为qL(h)。
[0165] 403.结合上述步骤402中qL(h)，通过步骤103的步骤303中的第三个公式，最终得到二阶空间偏导数结果
[0166] 404.通过步骤401～步骤403，计算x方向的空间偏导数
[0167] 405.通过步骤401～步骤403，计算y方向的空间偏导数
[0168] 406.通过步骤401～步骤403，计算z方向的空间偏导数
[0169] 作为本发明的一个实施例，还可以通过改进追赶法求解L个三对角线性方程组，例如图2所示为所述改进追赶法的流程图。
[0170] 所述改进追赶法包括：
[0171] 601.预先计算求解L个三对角线性方程组所需的L个常向量αl，通过以下递推公式来实现：
[0172]
[0173]
[0174] 其中dl为前面所求的隐式差分系数；αl表示采样长度为Nh的向量，Nh取空间采样网格数Nx、Ny和Nz中的最大值，Nx、Ny和Nz分别为x、y和z坐标轴方向的空间网格数。
[0175] 在本发明中，如果使用常规的追赶法，需要反复大量地计算向量αl，降低了计算效率；由于αl为常数向量，只需要预先计算一次即可，其计算量在整个有理分式有限差分计算过程中可以忽略不计。通过这一步骤，三对角线性方程组的浮点数乘除运算量可以由5Nh降低到3Nh。
[0176] 602.实现改进追赶法的追过程，通过以下递推公式进行：
[0177] w(1)＝αl(1)ql-1(1)
[0178] w(i)＝αl(i)[ql-1(i)-w(i-1)],(i=1,2,…,Nh)
[0179] 其中w为实现追赶法引入的中间变量，ql-1为前面所述第l-1阶隐式差分方程中的变量，αl为预先计算的常向量。上述公式由ql-1计算w。使用常规的追赶法，浮点数乘除法的运算量的追过程需要3Nh；结合所述步骤601，改进的追赶法这一过程只需要Nh次。
[0180] 603.实现改进追赶法的赶过程，通过以下递推公式进行：
[0181] ql(N)＝w(N)
[0182] ql(i)＝w(i)-αl(i)ql(i+1),(i＝Nh-1,Nh-2,…,1)
[0183] 其中w为实现追赶法引入的中间变量，ql为前面所述第l-1阶隐式差分方程中的变量，αl为预先计算的常向量。上述公式由w计算ql。该赶过程的实现所需浮点数乘除法运算量为Nh-1。
[0184] 通过所述改进的追赶法,的浮点数乘除法的运算量由大约5Nh次降低到大约2Nh次，计算效率比常规的追赶法提高了2.5倍。
[0185] 步骤105，在时间方向上进行迭代实现地震波场的数值模拟，通过如下公式进行：
[0186]
[0187] 其中pn表示t＝nΔt时刻传播的波场、pn-1表示前一时刻传播的波场，pn+1表示后一时刻传播的波场，Δt为时间采样间隔，Vp为弹性介质的纵波速度。时间方向上的每一次迭代计算需通过所述步骤104中的所述步骤404～406来完成，共需要迭代Nt次。
[0188] 实施例1
[0189] 理论模型为一个1500m/s的二维均匀介质模型。网格的空间采样步长为20mx20m，网格大小为150x150；时间采样步长为1ms，
震源使用雷克子波激发，峰值频率为18Hz。采用四阶有理分式有限差分法去实现二维的地震波数值模拟。具体实现步骤为：
[0190] 1.在波数域使用四阶的有理分式有限差分算子去逼近64阶差分精度的显式高阶有限差分算子；
[0191] 2.计算64阶差分精度的显式高阶有限差分算子的系数bm(m＝0,1,…,32)；
[0192] 3.通过曲线拟合法计算四阶有理分式有限差分算子的系数c1,c2,d1,d2，其结果为：
[0193]
[0194] 其中
[0195] c1＝-12.4603758,c2＝-12.6803855,
[0196] d1＝11.8369989,d2＝2.5880464.
[0197] 4.在空间域使用四阶有理分式有限差分法去实现地震波数值模拟。首先计算x坐标轴方向的空间偏导数，然后计算z坐标轴方向的空间偏导数。然后在时间方向上进行迭代完成地震波的数值模拟。
[0198] 其中计算x坐标轴方向的空间偏导数的步骤为：
[0199] ①结合系数c1和c2，计算x坐标轴方向的四阶显式有限差分；
[0200] q0,x(i,k)＝2(c1+c2)p(i,k)+c1[p(i-1,k)+p(i+1,k)]+c2[p(i-2,k)+p(i+2,k)][0201] ②结合系数d1，使用改进的追赶法求解第一个二阶隐式差分方程；
[0202] q1,x(i-1,k)+d1q1,x(i,k)+q1,x(i+1,k)＝q0,x(i,k),(i＝1,2,…,Nx)[0203] ③结合系数d2，使用改进的追赶法求解第二个二阶隐式差分方程；
[0204] q2,x(i-1,k)+d2q2,x(i,k)+q2,x(i+1,k)＝q1,x(i,k),(i＝1,2,…,Nx)[0205] ④获取x坐标轴方向的空间二阶偏导数：
[0206]
[0207] 仿照步骤①～④，可以计算z坐标轴方向的空间偏导数。
[0208] 时间方向上的迭代计算公式为：
[0209]
[0210] 图4a和图4b为四阶有理分式差分算子拟合64阶差分精度的高阶显式差分算子示意图。图4a为两个算子对应的响应曲线，在该图中两条曲线几乎重合难以区分，图4b为拟合相对误差曲线，在该图中拟合误差曲线表明最大相对拟合误差不超过0.5%，说明两种算子非常接近、吻合程度非常高，即四阶有理分式有限差分算子具有64阶差分精度。
[0211] 图5a-图5d为本发明实施例中不同有限差分方法所模拟的波场快照的比较示意图，图6a-图6d为本发明实施例中不同有限差分方法所模拟的单道记录的比较示意图。其中如图5a和图6a所示为采用伪谱法模拟的波场快照，数值理论解采用传统的伪谱法实施，其计算精度最高、几乎无数值频散，但计算效率非常低下；如图5b和图6b所示为4阶显式差分法模拟的波场快照，四阶显式有限差分法计算效率很高，但数值频散非常严重、地震子波拖尾很长；如图5c和图6c所示为12阶显式差分法模拟的波场快照，12阶差分精度的高阶显式差分法在一定程度上降低了数值频散，但地震波形的保幅效果不好；如图5d和图6d所示为四阶有理分式有限差分法模拟的波场快照，本发明实施例所采用的方法在波场快照几乎看不见数值频散，在单道记录结果上与理论的数值解也非常接近。
[0212] 从阶数来讲，四阶有理分式有限差分与四阶显式有限差分相同；从计算效率来讲，四阶有理分式有限差分与12阶显式有限差分相同；从计算精度来讲，四阶有理分式有限差分与64阶显式有限差分相同，非常接近于数值理论解。
[0213] 实施例2
[0214] 理论模型为一个二维复杂介质模型。网格的空间采样步长为20mx20m，网格大小为340x250；时间采样步长为1ms，采样长度为4s，震源使用雷克子波激发，峰值频率为18Hz。
采用四阶有理分式有限差分法去实现二维的地震波数值模拟。具体实现步骤为与实施例1相同。
[0215] 图7为本发明实施例的复杂的速度介质模型示意图；
[0216] 图8为现有技术中使用12阶高阶显式有限差分法对复杂速度介质模型所模拟的单炮记录，图中的正上方和右方
能量较强处皆出现了高频振荡拖尾现象，数值频散比较严重；
[0217] 图9为本发明实施例中使用四阶有理分式有限差分差分法对复杂速度介质模型所模拟的单炮记录示，从剖面中几乎看不见数值频散现象，与图8中的方法相比，数值模拟精度得到了很大的改善。
[0218] 通过本发明实例中的方法，通过使用低阶有理分式有限差分法来实现地震波数值模拟，提高了地震波场数值模拟的精度和计算效率。
[0219] 从数值模拟的计算精度角度来讲，阶数很低的有理分式有限差分具备阶数很高的超高阶显式有限差分的数值模拟精度。从附图2可以看出，传统的显式高阶差分算子只有当阶数非常高的时候，才能很好地逼近理论的差分算子。附图3是使用四阶有理分式有限差分算子去逼近64阶差分精度的超高阶显式差分算子，可以看出，两者吻合得非常好，观察相对误差曲线，其最大值小于0.5%，说明低阶有理分式有限差分算子对超高阶显式有限差分算子的逼近程度非常高。
[0220] 从数值模拟的计算效率角度来将，由于使用了改进的追赶法，使得求解一个Nh阶三对角线性方程组的乘法运算量从5Nh次降低到2Nh次，在保持相同计算精度的前提下，有理分式有限差分法的计算效率要大大高于超高阶显式差分法的计算效率。以64阶差分精度为例，四阶有理分式有限差分需要计算一个4阶的显式差分和求解两个二阶隐式差分：计算Nh个采样点四阶显式差分涉及的乘除运算量为3Nh；采用改进的追赶法计算Nh个采样点的两个隐式差分涉及的乘除运算量为4Nh；计算Nh个采样点的64阶显式差分涉及的乘除运算量为33Nh。因此，四阶有理分式有限差分的计算效率是64阶超高阶显式差分计算效率的4.7倍。
[0221] 如果保持相同的计算效率，四阶有理分式有限差分的计算量相当于12阶显式有限差分法的计算量，但计算精度却提高明显、数值频散大大降低。
[0222] 本发明可以以任何适当的形式实现，包括
硬件、
软件、
固件或它们的任意组合。本发明可以根据情况有选择的部分实现，比如计算机软件执行于一个或多个
数据处理器以及数字
信号处理器。本文的每个实施例的元素和组件可以在物理上、功能上、逻辑上以任何适当的方式实现。事实上，一个功能可以在独立单元中、在一组单元中、或作为其他功能单元的一部分来实现。因此，该系统和方法既可以在独立单元中实现，也可以在物理上和功能上分布于不同的单元和处理器之间。
[0223] 在相关领域中的技术人员将会认识到，本发明的实施例有许多可能的
修改和组合，虽然形式略有不同，仍采用相同的基本机制和方法。为了解释的目的，前述描述参考了几个特定的实施例。然而，上述的说明性讨论不旨在穷举或限制本文所发明的精确形式。前文所示，许多修改和变化是可能的。所选和所描述的实施例，用以解释本发明的原理及其实际应用，用以使本领域技术人员能够最好地利用本发明和各个实施例的针对特定应用的修改、
变形。
