三维波动方程混合网格有限差分数值模拟方法及装置

摘要

本申请提供一种三维波动方程混合网格有限差分数值模拟方法及装置。所述方法包括：构建三维混合网格有限差分格式，所述三维混合网格有限差分格式包括三维直角坐标系中坐标轴上的网格点和非坐标轴上的网格点；根据所述三维混合网格有限差分格式，对三维波动方程进行有限差分离散，获得所述三维波动方程的有限差分离散方程；根据所述有限差分离散方程和平面波理论，计算所述有限差分离散方程的差分系数；根据所述差分系数求解所述有限差分离散方程，实现所述三维波动方程的数值模拟。利用本申请所述的混合网格有限差分数值模拟方法有效降低了三维波动方程数值模拟的数值频散，提高了三维波动方程的数值模拟精度。

权利要求

1.一种地震波三维波动方程混合网格有限差分数值模拟方法，其特征在于，构建三维混合网格有限差分格式，所述三维混合网格有限差分格式包括三维直角坐标系中坐标轴上的网格点和非坐标轴上的网格点；

根据所述三维混合网格有限差分格式，对三维波动方程进行有限差分离散，获得所述三维波动方程的有限差分离散方程；

根据所述有限差分离散方程和平面波理论，计算所述有限差分离散方程的差分系数；

根据所述差分系数求解所述有限差分离散方程，实现所述三维波动方程的数值模拟；

所述根据所述三维混合网格有限差分格式，对三维波动方程进行有限差分离散，获得三维波动方程的有限差分离散方程，包括：若所述非坐标轴上的网格点位于所述三维直角坐标系的坐标平面内，则将所述非坐标轴上的网格点与所述三维混合网格有限差分格式中的差分中心点进行差分离散，获得所述非坐标轴上的网格点对应的二维拉普拉斯算子；

根据所述二维拉普拉斯算子获得所述非坐标轴上的网格点对应的三维拉普拉斯算子；

利用所述非坐标轴上的网格点对应的三维拉普拉斯算子和所述坐标轴上的网格点对应的三维拉普拉斯算子，计算获得所述三维波动方程的拉普拉斯算子；

根据所述三维波动方程的拉普拉斯算子对所述三维波动方程进行有限差分离散，获得所述三维波动方程的有限差分离散方程。

2.如权利要求1所述的方法，其特征在于，所述将所述非坐标轴上的网格点与所述三维混合网格有限差分格式中的差分中心点进行差分离散，获得所述非坐标轴上的网格点对应的二维拉普拉斯算子，包括：将所述三维直角坐标系的位于同一个坐标平面内的所述非坐标轴上的网格点与所述差分中心点进行差分离散，获得三个坐标平面分别对应的二维拉普拉斯算子；

相应地，根据所述二维拉普拉斯算子获得所述非坐标轴上的网格点对应的三维拉普拉斯算子，包括：将所述三个二维拉普拉斯算子相加，获得所述非坐标轴上的网格点对应的三维拉普拉斯算子。

3.如权利要求1所述的方法，其特征在于，所述根据所述三维混合网格有限差分格式，对三维波动方程进行有限差分离散，获得三维波动方程的有限差分离散方程，包括：若所述非坐标轴上的网格点位于所述三维直角坐标系的坐标平面外，则将所述非坐标轴上的网格点应用三元函数的泰勒级数展开；

将所述非坐标轴上的网格点的泰勒级数展开结果相加后与所述三维混合网格有限差分格式中的差分中心点进行差分离散，获得所述非坐标轴上的网格点对应的三维拉普拉斯算子；

利用所述非坐标轴上的网格点对应的三维拉普拉斯算子和所述坐标轴上的网格点对应的三维拉普拉斯算子，计算获得所述三维波动方程的拉普拉斯算子；

根据所述三维波动方程的拉普拉斯算子对所述三维波动方程进行有限差分离散，获得所述三维波动方程的有限差分离散方程。

4.如权利要求2-3任一项所述的方法，其特征在于，所述利用所述非坐标轴上的网格点对应的三维拉普拉斯算子和所述坐标轴上的网格点对应的三维拉普拉斯算子，计算获得所述三维波动方程的拉普拉斯算子，包括：将所述非坐标轴上的网格点对应的三维拉普拉斯算子和所述坐标轴上的网格点对应的三维拉普拉斯算子进行加权平均，获得所述三维波动方程的拉普拉斯算子。

5.如权利要求1所述的方法，其特征在于，所述方法还包括：

获得所述差分系数后，对所述有限差分离散方程进行频散分析，获得所述有限差分离散方程的数值频散；

判断所述数值频散是否大于预设频散阈值，若是，则调整所述三维混合网格有限差分格式；

根据调整后的三维混合网格差分格式，对所述三维波动方程进行有限差分离散，获得调整后的三维混合网格差分格式对应的调整有限差分离散方程、调整差分系数；

根据所述调整差分系数、调整有限差分离散方程，重新进行频散分析，并判断所述数值频散是否大于预设频散阈值，若是，则继续调整所述三维混合网格有限差分格式，直至所述频散数值小于等于所述预设频散阈值；

将所述频散数值小于等于所述预设频散阈值时对应的调整差分系数、调整有限差分离散方程作为数值模拟的差分系数和有限差分离散方程。

6.如权利要求1所述的方法，其特征在于，所述方法还包括：

利用所述三维波动方程的数值模拟结果至少用于：优化野外地震观测系统、检验处理方法的合理性、验证解释结果的正确性、直接应用于逆时偏移和全波形反演。

7.一种地震波三维波动方程混合网格有限差分数值模拟装置，其特征在于，包括：差分格式构建模块，用于构建三维混合网格有限差分格式，所述三维混合网格有限差分格式包括三维直角坐标系中坐标轴上的网格点和非坐标轴上的网格点；

有限差分方程建立模块，用于根据所述三维混合网格有限差分格式，对三维波动方程进行有限差分离散，获得所述三维波动方程的有限差分离散方程；

有限差分系数计算模块，用于根据所述有限差分离散方程和平面波理论，计算所述有限差分离散方程的差分系数；

数值模拟模块，用于根据所述差分系数求解所述有限差分离散方程，实现所述三维波动方程的数值模拟；

所述根据所述三维混合网格有限差分格式，对三维波动方程进行有限差分离散，获得三维波动方程的有限差分离散方程的方法包括：若所述非坐标轴上的网格点位于所述三维直角坐标系的坐标平面内，则将所述非坐标轴上的网格点与所述三维混合网格有限差分格式中的差分中心点进行差分离散，获得所述非坐标轴上的网格点对应的二维拉普拉斯算子；

根据所述二维拉普拉斯算子获得所述非坐标轴上的网格点对应的三维拉普拉斯算子；

利用所述非坐标轴上的网格点对应的三维拉普拉斯算子和所述坐标轴上的网格点对应的三维拉普拉斯算子，计算获得所述三维波动方程的拉普拉斯算子；

根据所述三维波动方程的拉普拉斯算子对所述三维波动方程进行有限差分离散，获得所述三维波动方程的有限差分离散方程。

8.如权利要求7所述的装置，其特征在于，所述三维波动方程混合网格有限差分数值模拟装置还包括：频散分析模块，用于在有限差分系数计算模块获得所述差分系数后，对所述有限差分离散方程进行频散分析，获得所述有限差分离散方程的数值频散；

判断所述数值频散是否大于预设频散阈值，若是，则差分格式构建模块调整所述三维混合网格有限差分格式；

相应地，所述有限差分方程建立模块根据调整后的三维混合网格差分格式，对所述三维波动方程进行有限差分离散，获得调整后的三维混合网格差分格式对应的调整有限差分离散方程、调整差分系数；

所述频散分析模块用于根据所述调整差分系数、调整有限差分离散方程，重新进行频散分析，并判断所述数值频散是否大于预设频散阈值，若是，则所述差分格式构建模块继续调整所述三维混合网格有限差分格式，直至所述频散数值小于等于所述预设频散阈值；

将所述频散数值小于等于所述预设频散阈值时对应的调整差分系数、调整有限差分离散方程作为数值模拟的差分系数和有限差分离散方程。

9.一种地震波三维波动方程混合网格有限差分数值模拟装置，其特征在于，包括处理器以及用于存储处理器可执行指令的存储器，所述处理器执行所述指令时实现如权利要求

1至6中任意一项所述方法的步骤。

说明书

三维波动方程混合网格有限差分数值模拟方法及装置

技术领域

[0001] 本申请属于地震波数值模拟技术领域，尤其涉及一种三维波动方程混合网格有限差分数值模拟方法及装置。

背景技术

[0002] 波动方程是一种重要的偏微分方程，它通常描述波在介质中的传播过程。波动方程数值模拟是地震勘探和地震学的重要基础，可以模拟研究地震波在地下各种介质中的传播规律。波动方程数值模拟可以用于野外观测系统的设计和评估、可以检验各种反演方法的正确性、可以对地震解释结果的正确性进行检验，是逆时偏移和全波形反演的基础和关键环节。

[0003] 现有技术中，有限差分法是目前应用最普遍的一种波动方程数值模拟方法，针对三维波动方程有限差分数值模拟，现有技术中的传统2M(M可以表示大于0的任意整数)阶和时空域2M阶有限差分格式，只利用了坐标轴上的网格点差分近似波动方程中的Laplace(拉普拉斯)算子。在进行波动方程的数值模拟时，主要通过增大M的取值来提高数值模拟的精度，然而随着M取值的增大，新增加的在坐标轴上的网格点距离中心的距离越来越远，对提高模拟精度的贡献越来越小，并且存在较严重的数值频散，模拟精度较低。

发明内容

[0004] 本申请目的在于提供一种三维波动方程混合网格有限差分数值模拟方法及装置，综合利用了三维直角坐标系中坐标轴上的网格点和非坐标轴上的网格点进行差分近似出波动方程的Laplace算子，实现对波动方程的有限差分数值模拟，提高了三维波动方程有限差分数值模拟精度。

[0005] 一方面本申请提供了一种三维波动方程混合网格有限差分数值模拟方法，包括：

[0006] 构建三维混合网格有限差分格式，所述三维混合网格有限差分格式包括三维直角坐标系中坐标轴上的网格点和非坐标轴上的网格点；

[0007] 根据所述三维混合网格有限差分格式，对三维波动方程进行有限差分离散，获得所述三维波动方程的有限差分离散方程；

[0008] 根据所述有限差分离散方程和平面波理论，计算所述有限差分离散方程的差分系数；

[0009] 根据所述差分系数求解所述有限差分离散方程，实现所述三维波动方程的数值模拟。

[0010] 进一步地，所述方法的另一个实施例中，所述根据所述三维混合网格有限差分格式，对三维波动方程进行有限差分离散，获得三维波动方程的有限差分离散方程，包括：

[0011] 若所述非坐标轴上的网格点位于所述三维坐直角坐标系的坐标平面内，则将所述非坐标轴上的网格点与所述三维混合网格有限差分格式中的差分中心点进行差分离散，获得所述非坐标轴上的网格点对应的二维拉普拉斯算子；

[0012] 根据所述二维拉普拉斯算子获得所述非坐标轴上的网格点对应的三维拉普拉斯算子；

[0013] 利用所述非坐标轴上的网格点对应的三维拉普拉斯算子和所述坐标轴上的网格点对应的三维拉普拉斯算子，计算获得所述三维波动方程的拉普拉斯算子；

[0014] 根据所述三维波动方程的拉普拉斯算子对所述三维波动方程进行有限差分离散离散，获得所述三维波动方程的有限差分离散方程。

[0015] 进一步地，所述方法的另一个实施例中，所述将所述非坐标轴上的网格点与所述三维混合网格有限差分格式中的差分中心点进行差分离散，获得所述非坐标轴上的网格点对应的二维拉普拉斯算子，包括：

[0016] 将所述三维直角坐标系的位于同一个坐标平面内的所述非坐标轴上的网格点与所述差分中心点进行差分离散，获得三个坐标平面分别对应的二维拉普拉斯算子；

[0017] 相应地，根据所述二维拉普拉斯算子获得所述非坐标轴上的网格点对应的三维拉普拉斯算子，包括：

[0018] 将所述三个二维拉普拉斯算子相加，获得所述非坐标轴上的网格点对应的三维拉普拉斯算子。

[0019] 进一步地，所述方法的另一个实施例中，所述根据所述三维混合网格有限差分格式，对三维波动方程进行有限差分离散，获得三维波动方程的有限差分离散方程，包括：

[0020] 若所述非坐标轴上的网格点位于所述三维直角坐标系的坐标平面外，则将所述非坐标轴上的网格点应用三元函数的泰勒级数展开；

[0021] 将所述非坐标轴上的网格点的泰勒级数展开结果相加后与所述三维混合网格有限差分格式中的差分中心点进行差分离散离散，获得所述非坐标轴上的网格点对应的三维拉普拉斯算子；

[0022] 利用所述非坐标轴上的网格点对应的三维拉普拉斯算子和所述坐标轴上的网格点对应的三维拉普拉斯算子，计算获得所述三维波动方程的拉普拉斯算子；

[0023] 根据所述三维波动方程的拉普拉斯算子对所述三维波动方程进行有限差分离散离散，获得所述三维波动方程的有限差分离散方程。

[0024] 进一步地，所述方法的另一个实施例中，所述利用所述非坐标轴上的网格点对应的三维拉普拉斯算子和所述坐标轴上的网格点对应的三维拉普拉斯算子，计算获得所述三维波动方程的拉普拉斯算子，包括：

[0025] 将所述非坐标轴上的网格点对应的三维拉普拉斯算子和所述坐标轴上的网格点对应的三维拉普拉斯算子进行加权平均，获得所述三维波动方程的拉普拉斯算子。

[0026] 进一步地，所述方法的另一个实施例中，所述方法还包括：

[0027] 获得所述差分系数后，对所述有限差分离散方程进行频散分析，获得所述有限差分离散方程的数值频散；

[0028] 判断所述数值频散是否大于预设频散阈值，若是，则调整所述三维混合网格有限差分格式；

[0029] 根据调整后的三维混合网格差分格式，对所述三维波动方程进行有限差分离散，获得调整后的三维混合网格差分格式对应的调整有限差分离散方程、调整差分系数；

[0030] 根据所述调整差分系数、调整有限差分离散方程，重新进行频散分析，并判断所述数值频散是否大于预设频散阈值，若是，则继续调整所述三维混合网格有限差分格式，直至所述频散数值小于等于所述预设频散阈值；

[0031] 将所述频散数值小于等于所述预设频散阈值时对应的调整差分系数、调整有限差分离散方程作为数值模拟的差分系数和有限差分离散方程。

[0032] 进一步地，所述方法的另一个实施例中，所述方法还包括：

[0033] 利用所述三维波动方程的数值模拟结果至少用于：优化野外地震观测系统、检验处理方法的合理性、验证解释结果的正确性、直接应用于逆时偏移和全波形反演。

[0034] 另一方面，本申请提供了一种三维波动方程混合网格有限差分数值模拟装置，包括：

[0035] 差分格式构建模块，用于构建三维混合网格有限差分格式，所述三维混合网格有限差分格式包括三维直角坐标系中坐标轴上的网格点和非坐标轴上的网格点；

[0036] 有限差分方程建立模块，用于根据所述三维混合网格有限差分格式，对三维波动方程进行有限差分离散，获得所述三维波动方程的有限差分离散方程；

[0037] 有限差分系数计算模块，用于根据所述有限差分离散方程和平面波理论，计算所述有限差分离散方程的差分系数；

[0038] 数值模拟模块，用于根据所述差分系数求解所述有限差分离散方程，实现所述三维波动方程的数值模拟。

[0039] 进一步地，所述装置的另一个实施例中，所述三维波动方程混合网格有限差分数值模拟装置还包括：

[0040] 频散分析模块，用于在有限差分系数计算模块获得所述差分系数后，对所述有限差分离散方程进行频散分析，获得所述有限差分离散方程的数值频散；

[0041] 判断所述数值频散是否大于预设频散阈值，若是，则差分格式构建模块调整所述三维混合网格有限差分格式；

[0042] 相应地，所述有限差分方程建立模块根据调整后的三维混合网格差分格式，对所述三维波动方程进行有限差分离散，获得调整后的三维混合网格差分格式对应的调整有限差分离散方程、调整差分系数；

[0043] 所述频散分析模块用于根据所述调整差分系数、调整有限差分离散方程，重新进行频散分析，并判断所述数值频散是否大于预设频散阈值，若是，则所述差分格式构建模块继续调整所述三维混合网格有限差分格式，直至所述频散数值小于等于所述预设频散阈值；

[0044] 将所述频散数值小于等于所述预设频散阈值时对应的调整差分系数、调整有限差分离散方程作为数值模拟的差分系数和有限差分离散方程。

[0045] 再一方面，本申请还提供了一种三维波动方程混合网格有限差分数值模拟装置，包括：处理器以及用于存储处理器可执行指令的存储器，所述处理器执行所述指令时实现上述三维波动方程混合网格有限差分数值模拟方法。

[0046] 本申请提供的三维波动方程混合网格有限差分数值模拟方法及装置，综合利用了三维直角坐标系中坐标轴上的网格点和非坐标轴上的网格点，构建出三维混合网格有限差分格式。基于三维混合有限差分格式，可以将坐标轴上的网格点以及非坐标轴上的网格点分别与差分中心点差分近似出波动方程的Laplace算子，并将坐标轴上的网格点对应的Laplace算子与非坐标轴上的网格点对应的Laplace算子进行加权平均，获得三维波动方程的Laplace算子，进一步获得三维波动方程的有限差分离散方程。根据获得的有限差分离散方程和平面波解，求解出有限差分离散方程的差分系数，进一步可以对有限差分离散方程进行迭代求解，实现三维波动方程的数值模拟。本说明书实施例提供的三维混合网格有限差分，充分利用了距离差分中心点更近的位于非坐标轴上的网格点，差分近似波动方程中的三维Laplace算子，理论上更为合理，进一步降低了三维波动方程数值模拟的数值频散，提高了三维波动方程的数值模拟精度。

附图说明

[0047] 为了更清楚地说明本申请实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本申请中记载的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。

[0048] 图1是本申请提供的一种三维波动方程混合网格有限差分数值模拟方法一个实施例的方法流程示意图；

[0049] 图2是本申请实施例中常规的2M阶有限差分格式的结构示意图；

[0050] 图3(a)是本申请实施例中非坐标轴上的网格点(12个非坐标轴上的网格点)构建的有限差分格式的结构示意图；

[0051] 图3(b)是本申请实施例中非坐标轴上的网格点(8个非坐标轴上的网格点)构建的有限差分格式的结构示意图；

[0052] 图3(c)是本申请实施例中非坐标轴上的网格点(24个非坐标轴上的网格点)构建的有限差分格式的结构示意图；

[0053] 图4(a)是本申请实施例中2M+N(N＝1)型三维混合网格有限差分格式的结构示意图；

[0054] 图4(b)是本申请实施例中2M+N(N＝2)型三维混合网格有限差分格式的结构示意图；

[0055] 图4(c)是本申请实施例中2M+N(N＝3)型三维混合网格有限差分格式的结构示意图；

[0056] 图5是本申请又一个实施例中三维波动方程混合网格有限差分数值模拟方法的流程示意图；

[0057] 图6(a)是本申请实施例中2M+N(M＝6，N＝1)型三维混合网格有限差分格式对应的频散关系曲线示意图；

[0058] 图6(b)是本申请实施例中2M+N(M＝6，N＝2)型三维混合网格有限差分格式对应的频散关系曲线示意图；

[0059] 图6(c)是本申请实施例中2M+N(M＝6，N＝3)型三维混合网格有限差分格式对应的频散关系曲线示意图；

[0060] 图6(d)是本申请实施例中2M+N(M＝12，N＝1)型三维混合网格有限差分格式对应的频散关系曲线示意图；

[0061] 图6(e)是本申请实施例中2M+N(M＝12，N＝2)型三维混合网格有限差分格式对应的频散关系曲线示意图；

[0062] 图6(f)是本申请实施例中2M+N(M＝12，N＝3)型三维混合网格有限差分格式对应的频散关系曲线示意图；

[0063] 图7(a)是本申请实施例中常规的2M(M＝8)阶有限差分格式频散关系曲线示意图；

[0064] 图7(b)是本申请实施例中时空域的2M(M＝8)阶有限差分格式频散关系曲线示意图；

[0065] 图7(c)是本申请实施例中2M+N(M＝6；N＝1)型三维混合网格有限差分格式频散关系曲线示意图；

[0066] 图8(a)-图8(f)是本申请实施例中三种不同有限差分格式对应的介质模型的波场快照示意图；

[0067] 图9是本申请提供的三维波动方程混合网格有限差分数值模拟装置一个实施例的模块结构示意图；

[0068] 图10是本申请又一个实施例中三维波动方程混合网格有限差分数值模拟装置的结构示意图；

[0069] 图11是本申请提供的另一种三维波动方程混合网格有限差分数值模拟装置实施例的模块结构示意图。

具体实施方式

[0070] 为了使本技术领域的人员更好地理解本申请中的技术方案，下面将结合本申请实施例中的附图，对本申请实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本申请一部分实施例，而不是全部的实施例。基于本申请中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都应当属于本申请保护的范围。

[0071] 波动方程的数值模拟在地震勘探、油气开发过程中有着重要的作用，通过波动方程的数值模拟可以获得地震波在介质中的传播规律。本发明实施例采用的是有限差分法对三维波动方程进行数值模拟，有限差分法是一种用差分代替微分求解微分方程的方法。有限差分法的基本思想包括把连续的定解区域用有限个离散点构成的网格来代替，这些离散点可以称作网格的节点或可以称为网格点，把连续定解区域上的连续变量的函数用在网格上定义的离散变量函数来近似。

[0072] 本申请实施例采用三维混合网格有限差分格式对波动方程进行有限差分离散，构建出三维混合网格有限差分格式，在三维混合网格有限差分格式的基础上，结合三维混合网格中位于坐标轴上的网格点和非坐标轴上的网格点对波动方程进行差分离散离散，获得声波方程数值模拟的有限差分离散方程，提高了波动方程数值模拟的精度。

[0073] 图1是本申请提供的一种三维波动方程混合网格有限差分数值模拟方法一个实施例的方法流程示意图，本申请提供的三维波动方程混合网格有限差分数值模拟方法包括：

[0074] S1、构建三维混合网格有限差分格式，所述三维混合网格有限差分格式包括三维直角坐标系中坐标轴上的网格点和非坐标轴上的网格点。

[0075] 本申请一个实施例中采用三维混合网格有限差分格式进行波动方程的数值模拟，三维混合网格有限差分格式可以理解为一种有限差分的网格模型，可以使用三维混合网格有限差分格式中距离差分中心点距离相等的一组网格点差分近似波动方程中的空间偏导数，Laplace算子。三维混合网格有限差分格式中包括坐标轴上的网格点和非坐标轴上的网格点。

[0076] 图2是本申请实施例中常规的2M阶有限差分格式的结构示意图，M可以取任意大于0的整数，如图2所示，现有技术中，常规的有限差分格式中的网格点均位于坐标轴上，随着M取值的增大，新增加的网格点距离差分中心点的距离越来越远，对提高模拟精度的贡献越来越小。图3(a)-图3(c)是本申请实施例中非坐标轴上的网格点构建的有限差分格式的结构示意图，如图3(a)所示，图中包括12个非坐标轴上的网格点，分布在三个坐标轴平面(xoy，xoz，yoz)内，每个坐标轴平面包括4个网格点，且距离差分中心点的距离均为 如图3(b)所示，图中包括8个非坐标轴上的网格点，分布在以坐标原点即差分中心点为中心的立方体的8个顶点上，8个非坐标轴上的网格点距离差分中心点的距离均为 如图3(c)所示，图中包括24个非坐标轴上的网格点，分布在三个坐标轴平面(xoy，xoz，yoz)内，每个坐标轴平面包括8个网格点，且距离差分中心点的距离均为 图3(a)-图3(c)中的非坐标轴上的网格点按照距离差分中心点的距离依次增加，根据需要还可以利用距离差分中心点距离更远的非坐标轴上的网格点构造更为丰富的差分格式。

[0077] 将图3(a)-图3(c)中的非坐标轴上的网格点构建的有限差分格式与图2中的坐标轴上的网格点构建的有限差分格式进行组合，可以获得三维混合网格有限差分格式。图4(a)-图4(c)是本申请实施例中三维混合网格有限差分格式的结构示意图，图4(a)中包括图2中的坐标轴上的网格点和图3(a)中的非坐标轴上的网格点，图4(b)中包括图2中的坐标轴上的网格点和图3(a)、图3(b)中的非坐标轴上的网格点，图4(c)中包括图2中的坐标轴上的网格点和图3(a)、图3(b)、图3(c)中的非坐标轴上的网格点。可以看出，本申请实施例中的三维混合网格有限差分格式可以根据需要将坐标轴上的网格点和非坐标轴上的网格点进行组合，可以利用三维混合网格有限差分格式中的坐标轴上的网格点和非坐标轴上的网格点差分近似波动方程中的Laplace算子。图2至图4(a)-图4(c)中带有数字的圆圈可以表示网格点。

[0078] 三维混合网格有限差分格式中位于非坐标轴上的网格点的数量通常是4的倍数，每次增加的非坐标轴上的网格点都是与差分中心点距离相等的一组网格点，并且按照距离从近到远，依次增加的。需要说明的是，三维混合网格有限差分格式中的三维坐标轴可以根据需要建立，本申请实施例不作具体限定。

[0079] S2、根据所述三维混合网格有限差分格式，对三维波动方程进行有限差分离散，获得所述三维波动方程的有限差分离散方程。

[0080] 三维直角坐标系中的时间-空间域常密度波动方程可以表示为公式(1)：

[0081]

[0082] 上式中，P＝P(x,y,z,t)可以表示标量声波波场，v＝v(x,y,z)可以表示声波在介质中的传播速度。

[0083] 则， 表示波动方程的三维Laplace算子。

[0084] 对公式(1)中的波动方程中的二阶时间偏导数进行二阶有限差分离散近似，可以得到：

[0085]

[0086] 上式中， m、l、n分别表示离散的三维空间坐标，j表示离散的时间采样坐标。h表示空间采样间隔，τ表示时间采样间隔，则
表示在任意参考时刻t和参考位置(x,y,z)处的波场值。

[0087] 可以看出，波动方程是一个偏微分方程，本申请实施例利用差分代替微分近似求解波动方程，即可以基于三维混合网格有限差分格式，对三维波动方程进行有限差分离散。可以利用三维混合网格有限差分格式中的坐标轴上的网格点和非坐标轴上的网格点差分近似出波动方程的拉普拉斯算子即Laplace算子。将三维混合网格有限差分格式中的坐标轴上的网格点和非坐标轴上的网格点差分近似获得的Laplace算子代入波动方程中，结合上述公式(2)，可以获得三维波动方程的有限差分离散方程。

[0088] S3、根据所述有限差分离散方程和平面波理论，计算所述有限差分离散方程的差分系数。

[0089] 三维介质模型可以理解为波动方程数值模拟的实施介质，具体可以根据需要进行设置，可以预先设置声波(如：地震波)在三维介质模型中的传播速度。平面波理论可以表示在计算差分系数的过程中利用了波动方程的平面波解。基于三维混合网格有限差分格式，对三维波动方程进行有限差分离散，获得波动方程的有限差分离散方程后，可以根据有限差分离散方程和平面波理论(平面波解)求解出有限差分离散方程的差分系数。可以将波动方程的平面波解的离散形式代入有限差分离散方程，并将代入平面波解的有限差分离散方程进行变形处理，可以获得有限差分离散方程的差分系数。也可以参考二维波动方程有限差分数值模拟中，差分系数的计算方法。

[0090] S4、根据所述差分系数求解所述有限差分离散方程，实现所述三维波动方程的数值模拟。

[0091] 获得波动方程的有限差分离散方程和对应的差分系数后，给定初始条件，如：空间采样间隔、时间采样间隔、声波在三维介质模型中的传播速度等，对有限差分离散方程进行迭代求解，可以将有限差分离散方程的解，近似代替波动方程的解，获得地震波(声波)在三维介质模型中的传播规律等，实现三维波动方程的数值模拟。

[0092] 本申请一个实施例中，获得三维波波动方程的数值模拟结果后，可以用于优化野外地震观测系统、检验处理方法的合理性、验证解释结果的正确性、直接应用于逆时偏移和全波形反演。当然，根据实际需要还可以用于其他地震勘探或地震学中，本申请实施例不作具体限定。

[0093] 本申请实施例提供的三维波动方程混合网格有限差分数值模拟方法，提出了三维混合网格有限差分格式，综合利用了坐标轴上的网格点以及非坐标轴上的网格点进行差分近似出波动方程的Laplace算子，实现对波动方程的有限差分数值模拟，提高了三维波动方程有限差分数值模拟精度。

[0094] 在上述实施例的基础上，本申请一个实施例中，所述根据所述三维混合网格有限差分格式，对三维波动方程进行有限差分离散，获得三维波动方程的有限差分离散方程，可以包括：

[0095] 若所述非坐标轴上的网格点位于所述三维直角坐标系的坐标平面内，则将所述非坐标轴上的网格点与所述三维混合网格有限差分格式中的差分中心点进行差分离散，获得所述非坐标轴上的网格点对应的二维拉普拉斯算子；

[0096] 根据所述二维拉普拉斯算子获得所述非坐标轴上的网格点对应的三维拉普拉斯算子；

[0097] 利用所述非坐标轴上的网格点对应的三维拉普拉斯算子和所述坐标轴上的网格点对应的三维拉普拉斯算子，计算获得所述三维波动方程的拉普拉斯算子；

[0098] 根据所述三维波动方程的拉普拉斯算子对所述三维波动方程进行有限差分离散，获得所述三维波动方程的有限差分离散方程。

[0099] 如图3(a)-图3(c)，位于非坐标轴上的网格点可以位于坐标平面内(如图3(a)、图3(c))或坐标平面外(如图3(b))，若位于非坐标轴上的网格点分别设置在三维直角坐标系中三个坐标平面内，则在将非坐标轴上的网格点进行差分近似波动方程的Laplace算子时，可以将非坐标轴上的网格点的Laplace算子，分解为二维的Laplace算子。可以将非坐标轴上的网格点与三维混合网格有限差分格式中的差分中心点即三维坐标系的坐标原点进行差分离散，差分近似出非坐标轴上的网格点对应的波动方程的二维拉普拉斯算子。将非坐标轴上的网格点对应的二维拉普拉斯算子进行组合，可以获得非坐标轴上的网格点对应波动方程的三维拉普拉斯算子。

[0100] 本申请一个实施例中，当非坐标轴上的网格点位于坐标平面内时，将计算获得所述非坐标轴上的网格点对应的二维拉普拉斯算子，可以包括：

[0101] 将所述三维坐标轴的位于同一个坐标平面内的所述非坐标轴上的网格点与所述差分中心点进行差分离散，获得三个坐标平面分别对应的二维拉普拉斯算子。

[0102] 将非坐标轴上的网格点根据分布的坐标平面，分别进行差分离散。可以将位于同一个坐标平面内的非坐标轴上的网格点看作二维网格点，将位于同一个坐标平面内的非坐标轴上的网格点与差分中心点进行差分近似，获得同一坐标平面内的非坐标轴上的网格点差分近似出的二维Laplace算子表达式。

[0103] 例如：在图3(a)中有12个非坐标轴上的网格点和1个差分中心点，12个非坐标轴上的网格点分布在三个坐标平面xOy，yOz和zOx内，且距离差分中心点的距离均为 可以参考二维旋转直角坐标系的方式，利用坐标平面xOy内的4个非坐标轴上的网格点和差分中心点差分近似二维Laplace算子可以得到：

[0104]

[0105] 同样地，在坐标平面yOz内有：

[0106]

[0107] 在坐标平面zOx内有：

[0108]

[0109] 可以将三个坐标平面对应的三个二维拉普拉斯算子相加，即将上述公式(3)(4)(5)相加，可以获得利用图3(a)中非坐标轴上的网格点对应的三维Laplace算子，可以表示为：

[0110]

[0111] 同样的方法，可以获得图3(c)中24个位于坐标平面内的非坐标轴上的网格点与差分中心点差分近似的三维拉普拉斯算子，此处不再赘述。

[0112] 对于位于坐标轴上的网格点，可以采用图2的有限差分格式，利用现有技术可以获得图2中坐标轴上的网格点差分近似的Laplace算子，可以表示为：

[0113]

[0114] 上述各式中，h可以表示空间采样间隔，cm(m＝1,2,…,M)可以表示权系数，本质上可以理解为将Laplace算子表示为常规直角坐标系中M个Laplace算子的加权平均。

[0115] 结合非坐标轴上的网格点对应的波动方程的三维拉普拉斯算子和坐标轴上的网格点对应的波动方程的三维拉普拉斯算子，可以获得三维波动方程的拉普拉斯算子。将三维波动方程的拉普拉斯算子代入三维波动方程，可以获得三维波动方程的有限差分离散方程。

[0116] 例如：可以将非坐标轴上的网格点对应的波动方程的三维拉普拉斯算子和坐标轴上的网格点对应的波动方程的三维拉普拉斯算子相加，获得三维波动方程的拉普拉斯算子。本申请一个实施例中，也可以将非坐标轴上的网格点对应的波动方程的三维拉普拉斯算子和坐标轴上的网格点对应的波动方程的三维拉普拉斯算子进行加权平均，获得三维波动方程的拉普拉斯算子。

[0117] 例如：可以将图2所示的常规2M阶有限差分格式与图3(a)所示的非坐标轴上的网格点的差分格式进行组合，可以得到图4(a)所示的混合网格有限差分格式。图4(a)中三维混合网格有限差分格式中的波动方程的Laplace算子，可以将上述公式(6)和公式(7)进行加权平均，对应的Laplace的差分表达式可以表示为：

[0118]

[0119] 上式中，c1,c2,…,cM；c1,1,0可以表示权系数。

[0120] 将三维混合网格有限差分格式中网格点对应的三维Laplace算子即公式(8)，结合上述公式(2)，代入上述公式(1)中，可以获得图4(a)所示的三维混合网格有限差分格式对应的波动方程的有限差分离散方程，可以表示为如下：

[0121]

[0122] 上式中，am＝cm/m2；a1,1,0＝c1,1,0/4，a1,a2,…,aM；a1,1,0可以表示差分系数。

[0123] 对于图4(c)所示的三维混合网格有限差分格式中网格点对应的波动方程的近似Laplace算子以及对应的有限差分离散方程的计算方法，可以参考图4(a)的处理方法，此处不再赘述。

[0124] 对于非坐标轴上的网格点位于三维坐标平面外，如图3(b)所示，此时对于非坐标轴上的网格点的差分离散可以包括：

[0125] 将非坐标轴上的网格点应用三元函数的泰勒级数展开，例如：可以对图3(b)中位于以差分中心点为对称中心的8个位于坐标平面外的非坐标轴上的网格点中的 采用三元函数的泰勒展开，可以得到如下公式：

[0126]

[0127] 同样的方法，可以对图3(b)中其他7个顶点进行三元函数的泰勒展开，此处不再赘述。

[0128] 获得位于坐标平面外的非坐标轴的网格点的三元函数的泰勒展开后，可以对8个泰勒展开结果进行相加和整理。将8个非坐标轴上的网格点的泰勒展开结果相加后，可以推导出利用图3(b)中8个非坐标轴上的网格点和差分中心点一起差分近似三维Laplace算子的表达式为：

[0129]

[0130] 同样的，也可以将图2所示的常规2M阶有限差分格式与图3(b)所示的位于坐标平面外的非坐标轴上的网格点的差分格式进行组合，具体可以将上述公式(7)和公式(11)进行加权平均，获得对应的Laplace算子。具体过程可以参考上述实施例的介绍，此处不再赘述。

[0131] 图4(b)是将图2中的坐标轴上的网格点和图3(a)、图3(b)中的非坐标轴上的网格点组合，构建的三维混合网格有限差分格式。则相应的可以将图2中的常规直接坐标系中坐标轴上的网格点表示的M个Laplace算子即公式(7)，与图3(a)中位于非坐标轴上的网格点对应的Laplace算子即公式(6)、以及图3(b)中位于非坐标轴上的网格点对应的Laplace算子即公式(11)进行加权平均，可以获得图4(b)的三维混合网格有限差分格式对应的波动方程的Laplace算子。

[0132] 同样的，图4(c)是将图2中的坐标轴上的网格点和图3(a)、图3(b)、图3(c)中的非坐标轴上的网格点组合，构建的三维混合网格有限差分格式。则相应的可以将图2中的常规直角坐标系中坐标轴上的网格点表示的M个Laplace算子即公式(7)，与图3(a)、图3(b)、图3(c)中位于非坐标轴上的网格点对应的Laplace算子进行加权平均，可以获得图4(c)的三维混合网格有限差分格式对应的波动方程的Laplace算子。

[0133] 本申请实施例中，混合网格2M+N型有限差分格式的基本思想就是将Laplace算子表示为常规直角坐标系中坐标轴上的网格点表示的M个Laplace算子和非坐标轴上的网格点表示的N个Laplace算子的加权平均，这是本申请实施例中构建混合网格有限差分法的根本出发点。

[0134] 获得三维混合网格有限差分格式中网格点对应的波动方程的近似Laplace算子后，结合上述公式(2)，代入上述公式(1)中，可以获得三维混合网格有限差分格式对应的波动方程的有限差分离散方程，具体可以参考上述图4(a)对应的三维混合网格有限差分格式的波动方程的有限差分离散方程的建立方法。

[0135] 现有技术中的2M阶和时空域2M有限差分格式，只是将Laplace算子表示为常规直角坐标系中坐标轴上的网格点表示的M个Laplace算子的加权平均，本申请实施例中的三维混合网格2M+N型有限差分格式的基本思想是将Laplace算子表示为常规直角坐标系中坐标轴上的网格点表示的M个Laplace算子和非坐标轴上的网格点表示的N个Laplace算子的加权平均。图4(a)-图4(c)依次对应三维混合网格2M+N(N＝1，2，3)型有限差分格式，继续增大N的取值，可以得到更为丰富的混合网格有限差分格式。

[0136] 现有技术中的波动方程的有限差分格式，只利用了坐标轴上的网格点差分近似Laplace算子，主要通过增大M的取值来减小数值频散，提高模拟精度，但随着M取值的增大，新增的网格点距离有限差分点越来越远，对提高模拟精度的贡献越来越小。本申请实施例，综合了非坐标轴上的网格点对波动方程进行有限差分离散，位于非坐标轴上的网格点与有限差分点之间的距离增加的比较缓慢，可以通过增加非坐标轴上的网格点，提高波动方程的数值模拟精度。

[0137] 获得三维混合网格有限差分格式中网格点对应的Laplace算子后，结合上述公式(2)，可以获得基于三维混合网格有限差分格式的波动方程的有限差分离散方程，基于有限差分离散方程计算出有限差分离散方程的差分系数即上述实施例中的a1,a2,…,aM；a1,1,0。

[0138] 计算差分系数是混合网格有限差分数值模拟的关键环节和重要研究内容，声波方程在均匀介质中存在平面波解，本申请一个实施例中，三维混合网格有限差分格式对应的差分系数的求解方法可以基于平面波解和时空域频散关系的差分系数计算，具体可以参考如下：

[0139] 步骤1：可以将波动方程的平面波解的离散形式代入基于三维混合网格有限差分格式获得的三维波动方程的有限差分离散方程，获得三维波动方程的余弦有限差分离散方程。

[0140] 具体过程可以参考如下：

[0141] 波动方程即上述公式(1)的平面波解的离散形式可以表示为：

[0142]

[0143] kx＝k sinφcosθ,ky＝k sinφsinθ，ky＝k cosφ   (13)

[0144] 上式中，ω可以表示圆频率，k可以表示波数，φ可以表示平面波传播方向与三维混合网格有限差分格式中三维坐标轴中的z轴正方向的夹角，θ可以平面波传播方向在与三维混合网格有限差分格式中的三维坐标平面xOy平面的投影与x轴正方向的夹角，则上式中，m、l、n分别表示离散的三维空间坐标，j表示离散的时间采样坐标。h表示空间采样间隔，τ表示时间采样间隔，i可以表示网格点。

[0145] 将公式(12)代入上述公式(9)可以得到余弦有限差分离散方程，如下式：

[0146]

[0147] 上式中，r＝vτ/h，可以表示Courant(人名)稳定性条件数，v＝v(x,y,z)可以表示声波在介质中的传播速度，τ可以表示时间采样间隔。

[0148] 步骤2：对余弦有限差分离散方程中的余弦函数进行泰勒展开，获得三维波动方程的泰勒展开有限差分离散方程，如将上述公式(14)中的余弦函数进行泰勒展开，可以获得：

[0149]

[0150] 步骤3：将泰勒展开有限差分离散方程中任意两个坐标轴方向的波数的平方的乘积对应的系数设置为相等，简化泰勒展开有限差分离散方程，获得差分系数方程。例如：将上述公式(15)中左右两边的 (或 或 )的系数对应相等，可以得到：

[0151]

[0152] 步骤4：将泰勒展开有限差分离散方程中任意一个坐标轴方向的波数的偶数次幂对应的系数设置为相等，简化泰勒展开有限差分离散方程，获得差分系数方程二。例如：将上述公式(15)中左右两边的 (或 或 )的系数对应相等，可以得到：

[0153]

[0154] 步骤5：求解差分系数方程一和差分系数方程二，获得有限差分离散方程的差分系数。例如：方程(16)和(17)共M+1个方程，刚好可以求解出混合网格2M+N(N＝1)型有限差分格式的M+1个差分系数a1,a2,…,aM；a1,1,0。

[0155] 可以看出，差分系数与模型速度(即声波在介质中的传播速度)有关，对于变速模型，每个模型网格点可以对应一个速度，也将对应一组差分系数。

[0156] 利用相同的方法可以推导出传统2M阶，时空域2M阶，和其他混合网格有限差分格式的差分系数计算方法。

[0157] 表1给出了传统的2M(M＝7)阶，时空域2M(M＝7)阶，和混合网格2M+N(M＝5；N＝1，2，3)型共五种有限差分格式的差分系数。计算差分系数时采用的速度(v)，时间采样间隔(τ)，和空间采样间隔(h)参数分别为v＝3000m/s,τ＝0.001s,h＝10m。

[0158] 表1五种有限差分格式的差分系数表

[0159]

[0160] 获得三维波动方程的有限差分离散方程以及有限差分离散方程的差分系数后，可以对有限差分离散方程进行迭代求解，获得波动方程的解，即获得声波在介质中到的传播规律，实现波动方程的数值模拟。

[0161] 图5是本申请又一个实施例中三维波动方程混合网格有限差分数值模拟方法的流程示意图，如图5所示，本申请一个实施例中三维波动方程混合网格有限差分数值模拟方法还可以包括：

[0162] S20、构建三维混合网格有限差分格式。具体构建三维混合网格有限差分格式的方法可以从参考上述实施例的介绍，此处不再赘述。

[0163] S21、基于三维混合网格有限差分格式，对三维声波方程进行有限差分离散得到相应的有限差分离散方程。具体构建有限差分离散方程的方法可以从参考上述实施例的介绍，此处不再赘述。

[0164] S22、根据有限差分离散方法和平面波理论求解差分系数。具体的求解方法可以参考上述实施例的介绍，此处不再赘述。

[0165] S23、对有限差分离散方程进行频散分析，获得所有限差分离散方程的数值频散。数值频散可以分为时间数值频散和空间数值频散，时间数值频散使得相速度变大，模拟波场中会出现“相位超前”的现象；空间数值频散使得相速度变小，模拟波场中会出现“相位滞后”的现象。

[0166] 数值频散是有限差分法求解波动方程的固有特性，无法完全消除，只能减小，通常可以用数值频散的大小来衡量有限差分法的数值模拟精度。本申请一个实施例中，可以定义归一化相速度δ来描述有限差分格式的数值频散，参考下式：

[0167]

[0168] 上式中，vph可以表示相速度，v可以表示声波(如地震波)在介质中传播的真实速度。

[0169] 结合相速度的定义vph＝ω/k，公式(14)和(18)可以给出混合网格2M+N(N＝1)型有限差分格式的归一化相速度关系式：

[0170]

[0171]

[0172] 上式中，G＝λ/h，λ可以表示波长，则G可以表示单位波长内网格点的个数。

[0173] δ的值与1越接近，可以表明数值频散误差越小；δ＞1可以表示存在时间数值频散，相速度偏大，模拟波场中会出现“相位超前”的现象；δ＜1可以表示存在空间数值频散，相速度偏小，模拟波场中会出现“相位滞后”的现象。利用相同的方法可以推导出常规2M阶，时空域2M阶，和其他混合网格有限差分格式的频散关系。

[0174] S24、判断获得的数值频散是否大于预设频散阈值，若是则执行步骤S25，否则执行步骤S26。可以根据要求的波动方程的数值模拟精度设置预设频散阈值的大小，本申请实施例不作具体限定。

[0175] S25、调整三维混合网格有限差分格式，并返回步骤S21。即调整三维混合网格有限差分格式中的坐标轴上的网格点对应的Laplace算子的数量M和非坐标轴上的网格点对应的Laplace算子的数量N。

[0176] 调整三维混合网格有限差分格式后，返回步骤S21。即根据调整后的三维混合网格差分格式，对三维波动方程进行有限差分离散，获得调整后的三维混合网格差分格式对应的调整有限差分离散方程、调整差分系数，并进行新的频散分析，直至数值频散小于等于预设频散阈值。将频散数值小于等于预设频散阈值时对应的调整差分系数、调整有限差分离散方程作为数值模拟的差分系数和有限差分离散方程，执行步骤S26。

[0177] S26、根据差分系数求解有限差分离散方程，实现所述三维波动方程的数值模拟。具体方法可以参考上述实施例，此处不再赘述。

[0178] 通过频散分析，可以获得数值模拟精度符合要求的三维混合网格有限差分格式，并可以基于此时的三维混合网格有限差分格式对波动方程进行有限差分离散，实现三维波动方程的数值模拟。

[0179] 图6(a)-图6(f)是本申请实施例中不同的三维混合网格有限差分格式对应的频散关系曲线示意图，图中部分曲线重合度比较大。图6(a)、图6(b)、图6(c)分别可以表示三维混合网格2M+N(M＝6；N＝1，2，3)型有限差分格式对应的频散关系曲线；图6(d)、图6(e)、图6(f)分别可以表示三维混合网格2M+N(M＝12；N＝1，2，3)型有限差分格式对应的频散关系曲线。从公式(19)和(20)可以看出，数值频散与地震波的传播方向φ和θ有关，同时还与参数1/G有关。图6(a)-图6(f)可以表示不同有限差分格式在θ＝π/8；φ＝0，π/8，2π/8，3π/8，4π/
8时的归一化相速度δ随参数1/G变化的五条频散曲线，图中部分曲线重合度比较大。绘制频散曲线时采用的速度(v)，时间采样间隔(τ)，和空间采样间隔(h)参数可以分别为v＝
3000m/s,τ＝0.001s,h＝10m。

[0180] 图6(a)-图6(c)中的三维混合网格2M+N(M＝6；N＝1，2，3)型有限差分格式的频散曲线特征基本相同，数值模拟精度基本相当，但增大N的取值，会增加数值模拟的计算量，降低计算效率。因此，当M取值不是很大时，例如M＝6，可以采用三维混合网格2M+N(N＝1)型有限差分格式，这样既可以保证数值模拟精度，同时可以兼顾计算效率。

[0181] 对比图6(d)-图6(f)中的三维混合网格2M+N(M＝12；N＝1，2，3)型有限差分格式的频散曲线，可以看出，与三维混合网格2M+N(M＝12；N＝1，2)型有限差分格式相比，三维混合网格2M+N(M＝12；N＝3)型有限差分格式的数值频散误差明显减小。因此，当M取值较大时，例如M＝12，可以采用混合网格2M+N(N＝3)型有限差分格式来进一步提高模拟精度。

[0182] 图7(a)-图7(c)是本申请实施例中有限差分格式频散关系曲线对比示意图，图7(a)可以表示常规的2M(M＝8)阶有限差分格式对应的频散曲线关系示意图，图7(b)可以表示时空域2M(M＝8)阶有限差分格式对应的频散曲线关系示意图，图7(c)可以表示三维混合网格2M+N(M＝6；N＝1)型有限差分格式对应的频散曲线关系示意图。

[0183] 如图7(a)所示，图中部分曲线重合度比较大，传统的2M(M＝8)阶有限差分格式的归一化相速度δ约等于1对应1/G的取值区间很小，大致为(0，0.075)，并且随着1/G取值的增大，δ大于1的现象越来越明显，表明传统2M阶有限差分格式存在严重的时间数值频散。

[0184] 如图7(b)所示，图中部分曲线重合度比较大，时空域2M(M＝8)阶有限差分格式的归一化相速度δ约等于1对应1/G的取值区间大致为(0，0.125)，与传统2M(M＝8)阶有限差分格式相比，这一区间范围明显增大。并且时空域2M(M＝8)阶比传统2M(M＝8)阶有限差分格式的数值频散幅值明显减小，这两点均表明时空域2M(M＝8)阶比传统2M(M＝8)阶有限差分格式能更有效地压制数值频散。但是，时空域2M(M＝8)阶有限差分格式的频散曲线仍然较发散，表明时空域2M阶有限差分格式存在一定的时间和空间数值频散。

[0185] 如图7(c)所示，图中部分曲线重合度比较大，混合网格2M+N(M＝6；N＝1)型有限差分格式归一化相速度δ约等于1对应1/G的取值区间大致为(0，0.225)，它是时空域2M(M＝8)阶有限差分格式的1.8倍，是传统2M(M＝8)阶有限差分格式的3倍。因此混合网格2M+N型有限差分格式能更好地压制数值频散，具有最高的模拟精度，同时，它的频散曲线收敛性更好。

[0186] 因此，计算量基本相等时，传统2M阶有限差分格式模拟精度最低，具有严重的时间数值频散；时空域2M阶有限差分格式模拟精度中等，具有一定的时间和空间数值频散；三维混合网格2M+N型有限差分格式模拟精度最高，数值频散最小。

[0187] 图8(a)-图8(f)是本申请实施例中不同三维混合网格有限差分格式对应的介质模型的波场快照示意图，图8(a)-图8(f)中下面的图例可以表示色标，数字可以表示色标刻度。本申请一个实施例中图8(a)-图8(f)中的介质模型的规模可以为：nx×ny×nz＝1215×1215×1215，nx、ny和nz可以分别表示介质模型在x方向、y方向和z方向离散网格点的个数，声波在介质模型中的速度可以为v＝3000m/s，时间采样间隔τ＝0.001s，空间采样间隔h＝
10m，震源子波为30Hz雷克子波。介质模型是数值模拟实施的介质，利用传统2M(M＝8)阶、时空域2M(M＝8)阶和三维混合网格2M+N(M＝6；N＝1)型有限差分格式三种方法，在介质模型上进行数值模拟实验。

[0188] 图8(a)-图8(f)可以表示不同有限差分格式方法，在均匀介质模型(nx×ny×nz＝1215×1215×1215)上进行数值模拟，获得的三维体1/8部分(608≤nx≤1215；1≤ny≤608；
1≤nz≤608)和截面nx＝608(x轴中心截面)左上角1/4部分的波场快照。图8(a)、图8(b)可以分别表示传统的2M(M＝8)阶有限差分格式的三维体1/8部分和截面nx＝608左上角1/4部分的波场快照；图8(c)、图8(d)可以分别表示时空域2M(M＝8)阶有限差分格式三维体1/8部分和截面nx＝608左上角1/4部分的波场快照；图8(e)、图8(f)可以分别表示三维混合网格
2M+N(M＝6；N＝1)型有限差分格式三维体1/8部分和截面nx＝608左上角1/4部分的波场快照。如图8(a)-图8(f)所示，对比三种有限差分格式数值模拟产生的波场产生的波场快照中的数值频散，对比可以看出，三维混合网格2M+N(M＝6；N＝1)无明显的数值频散，模拟精度最高。

[0189] 本申请实施例提供的三维波动方程混合网格有限差分数值模拟方法，综合利用了位于坐标轴上的网格点以及非坐标轴上的网格点，构建出三维混合网格有限差分格式。基于三维混合有限差分格式，可以将坐标轴上的网格点以及非坐标轴上的网格点分别与差分中心点差分近似出波动方程的Laplace算子，并将坐标轴上的网格点对应的Laplace算子于非坐标轴上的网格点对应的Laplace算子进行加权平均，获得三为波动方程的Laplace算子，进一步获得三为波动方程的有限差分离散方程。根据获得的有限差分离散方程，求解出有限差分离散方程的差分系数，进一步可以对有限差分离散方程进行迭代求解，近似获得波动方程的解，完成三维波动方程的数值模拟。综合利用非坐标轴上的网格点近似差分出三维波动方程的Laplace算子，降低了三为波动方程数值模拟的数值频散，提高了三维波动方程的数值模拟精度。

[0190] 基于上述所述的三维波动方程混合网格有限差分数值模拟方法，本说明书一个或多个实施例还提供一种三维波动方程混合网格有限差分数值模拟装置。所述的装置可以包括使用了本说明书实施例所述方法的系统(包括分布式系统)、软件(应用)、模块、组件、服务器、客户端等并结合必要的实施硬件的装置。基于同一创新构思，本说明书实施例提供的一个或多个实施例中的装置如下面的实施例所述。由于装置解决问题的实现方案与方法相似，因此本说明书实施例具体的装置的实施可以参见前述方法的实施，重复之处不再赘述。以下所使用的，术语“单元”或者“模块”可以实现预定功能的软件和/或硬件的组合。尽管以下实施例所描述的装置较佳地以软件来实现，但是硬件，或者软件和硬件的组合的实现也是可能并被构想的。

[0191] 具体地，图9是本申请提供的三维波动方程混合网格有限差分数值模拟装置一个实施例的模块结构示意图，如图9所示，本申请中提供的三维波动方程混合网格有限差分数值模拟装置包括：差分格式构建模块91，有限差分方程建立模块92，有限差分系数计算模块93，数值模拟模块94。

[0192] 差分格式构建模块91，可以用于构建三维混合网格有限差分格式，所述三维混合网格有限差分格式包括坐标轴上的网格点和非坐标轴上的网格点；

[0193] 有限差分方程建立模块92，可以用于根据所述三维混合网格有限差分格式，对三维波动方程进行有限差分离散，获得所述三维波动方程的有限差分离散方程；

[0194] 有限差分系数计算模块93，可以用于根据所述有限差分离散方程和平面波理论，计算所述有限差分离散方程的差分系数；

[0195] 数值模拟模块94，可以用于根据所述差分系数求解所述有限差分离散方程，实现所述三维波动方程的数值模拟。

[0196] 本申请提供的三维波动方程混合网格有限差分数值模拟装置，提出了三维混合网格有限差分格式，综合利用了坐标轴上的网格点以及非坐标轴上的网格点进行差分近似出波动方程的Laplace算子，实现对波动方程的有限差分数值模拟，提高了三维波动方程有限差分数值模拟精度。

[0197] 图10是本申请又一个实施例中三维波动方程混合网格有限差分数值模拟装置的结构示意图，如图10所示，所述三维波动方程混合网格有限差分数值模拟装置还包括：

[0198] 频散分析模块101，用于在有限差分系数计算模块获得所述差分系数后，对所述有限差分离散方程进行频散分析，获得所述有限差分离散方程的数值频散；

[0199] 判断所述数值频散是否大于预设频散阈值，若是，则差分格式构建模块91调整所述三维混合网格有限差分格式；

[0200] 相应地，所述有限差分方程建立模块92用于根据调整后的三维混合网格差分格式，对所述三维波动方程进行有限差分离散，获得调整后的三维混合网格差分格式对应的调整有限差分离散方程、调整差分系数；

[0201] 相应地，所述频散分析模块101用于根据所述调整差分系数、调整有限差分离散方程，重新进行频散分析，并判断所述数值频散是否大于预设频散阈值，若是，则所述差分格式构建模块91继续调整所述三维混合网格有限差分格式，直至所述频散数值小于等于所述预设频散阈值；

[0202] 将所述频散数值小于等于所述预设频散阈值时对应的调整差分系数、调整有限差分离散方程作为数值模拟的差分系数和有限差分离散方程。本申请提供的三维波动方程混合网格有限差分数值模拟装置，给出了位于三维坐标轴的坐标平面内的非坐标轴上的网格点近似出波动方程的Laplace算子的方法，为后续获得三维波动方程的有限差分离散方程提供了准确的数据基础。

[0203] 本申请提供的三维波动方程混合网格有限差分数值模拟装置，通过频散分析调整三维混合网格有限差分格式，获得满足数值模拟精度的三维混合网格有限差分格式，进一步确保三维波动方程的数值模拟精度。

[0204] 需要说明的，上述所述的装置根据方法实施例的描述还可以包括其他的实施方式，具体的实现方式可以参照相关方法实施例的描述，在此不作一一赘述。

[0205] 上述对本说明书特定实施例进行了描述。其它实施例在所附权利要求书的范围内。在一些情况下，在权利要求书中记载的动作或步骤可以按照不同于实施例中的顺序来执行并且仍然可以实现期望的结果。另外，在附图中描绘的过程不一定要求示出的特定顺序或者连续顺序才能实现期望的结果。在某些实施方式中，多任务处理和并行处理也是可以的或者可能是有利的。

[0206] 本说明书实施例提供的上述三维波动方程混合网格有限差分数值模拟方法或装置可以在计算机中由处理器执行相应的程序指令来实现，如使用windows操作系统的c++语言在PC端实现、linux系统实现，或其他例如使用android、iOS系统程序设计语言在智能终端实现，以及基于量子计算机的处理逻辑实现等。本说明书提供的一种三维波动方程混合网格有限差分数值模拟装置的另一个实施例中，图11是本申请提供的另一种三维波动方程混合网格有限差分数值模拟装置实施例的模块结构示意图，如图11所示，本申请另一实施例提供的三维波动方程混合网格有限差分数值模拟装置可以包括处理器111以及用于存储处理器可执行指令的存储器112，

[0207] 处理器111和存储器112通过总线113完成相互间的通信；

[0208] 所述处理器111用于调用所述存储器112中的程序指令，以执行上述各三维波动方程混合网格有限差分数值模拟方法实施例所提供的方法，例如包括：构建三维混合网格有限差分格式，所述三维混合网格有限差分格式包括坐标轴上的网格点和非坐标轴上的网格点；根据所述三维混合网格有限差分格式，对三维波动方程进行有限差分离散，获得所述三维波动方程的有限差分离散方程；根据所述有限差分离散方程和三维介质模型中声波的传播速度，计算所述有限差分离散方程的差分系数；根据所述差分系数求解所述有限差分离散方程，实现所述三维波动方程的数值模拟。

[0209] 需要说明的是说明书上述所述的装置根据相关方法实施例的描述还可以包括其他的实施方式，具体的实现方式可以参照方法实施例的描述，在此不作一一赘述。本申请中的各个实施例均采用递进的方式描述，各个实施例之间相同相似的部分互相参见即可，每个实施例重点说明的都是与其他实施例的不同之处。尤其，对于硬件+程序类实施例而言，由于其基本相似于方法实施例，所以描述的比较简单，相关之处参见方法实施例的部分说明即可。

[0210] 本说明书实施例并不局限于必须是符合行业通信标准、标准计算机数据处理和数据存储规则或本说明书一个或多个实施例所描述的情况。某些行业标准或者使用自定义方式或实施例描述的实施基础上略加修改后的实施方案也可以实现上述实施例相同、等同或相近、或变形后可预料的实施效果。应用这些修改或变形后的数据获取、存储、判断、处理方式等获取的实施例，仍然可以属于本说明书实施例的可选实施方案范围之内。

[0211] 在20世纪90年代，对于一个技术的改进可以很明显地区分是硬件上的改进(例如，对二极管、晶体管、开关等电路结构的改进)还是软件上的改进(对于方法流程的改进)。然而，随着技术的发展，当今的很多方法流程的改进已经可以视为硬件电路结构的直接改进。设计人员几乎都通过将改进的方法流程编程到硬件电路中来得到相应的硬件电路结构。因此，不能说一个方法流程的改进就不能用硬件实体模块来实现。例如，可编程逻辑器件(Programmable Logic Device，PLD)(例如现场可编程门阵列(Field Programmable Gate Array，FPGA))就是这样一种集成电路，其逻辑功能由用户对器件编程来确定。由设计人员自行编程来把一个数字系统“集成”在一片PLD上，而不需要请芯片制造厂商来设计和制作专用的集成电路芯片。而且，如今，取代手工地制作集成电路芯片，这种编程也多半改用“逻辑编译器(logic compiler)”软件来实现，它与程序开发撰写时所用的软件编译器相类似，而要编译之前的原始代码也得用特定的编程语言来撰写，此称之为硬件描述语言
(Hardware Description Language，HDL)，而HDL也并非仅有一种，而是有许多种，如ABEL(Advanced Boolean Expression Language)、AHDL(Altera Hardware Description Language)、Confluence、CUPL(Cornell University Programming Language)、HDCal、JHDL(Java Hardware Description Language)、Lava、Lola、MyHDL、PALASM、RHDL(Ruby Hardware Description Language)等，目前最普遍使用的是VHDL(Very-High-Speed Integrated Circuit Hardware Description Language)与Verilog。本领域技术人员也应该清楚，只需要将方法流程用上述几种硬件描述语言稍作逻辑编程并编程到集成电路中，就可以很容易得到实现该逻辑方法流程的硬件电路。

[0212] 控制器可以按任何适当的方式实现，例如，控制器可以采取例如微处理器或处理器以及存储可由该(微)处理器执行的计算机可读程序代码(例如软件或固件)的计算机可读介质、逻辑门、开关、专用集成电路(Application Specific Integrated Circuit，ASIC)、可编程逻辑控制器和嵌入微控制器的形式，控制器的例子包括但不限于以下微控制器：ARC 625D、Atmel AT91SAM、Microchip PIC18F26K20以及Silicone Labs C8051F320，存储器控制器还可以被实现为存储器的控制逻辑的一部分。本领域技术人员也知道，除了以纯计算机可读程序代码方式实现控制器以外，完全可以通过将方法步骤进行逻辑编程来使得控制器以逻辑门、开关、专用集成电路、可编程逻辑控制器和嵌入微控制器等的形式来实现相同功能。因此这种控制器可以被认为是一种硬件部件，而对其内包括的用于实现各种功能的装置也可以视为硬件部件内的结构。或者甚至，可以将用于实现各种功能的装置视为既可以是实现方法的软件模块又可以是硬件部件内的结构。

[0213] 上述实施例阐明的系统、装置、模块或单元，具体可以由计算机芯片或实体实现，或者由具有某种功能的产品来实现。一种典型的实现设备为计算机。具体的，计算机例如可以为个人计算机、膝上型计算机、车载人机交互设备、蜂窝电话、相机电话、智能电话、个人数字助理、媒体播放器、导航设备、电子邮件设备、游戏控制台、平板计算机、可穿戴设备或者这些设备中的任何设备的组合。

[0214] 虽然本说明书一个或多个实施例提供了如实施例或流程图所述的方法操作步骤，但基于常规或者无创造性的手段可以包括更多或者更少的操作步骤。实施例中列举的步骤顺序仅仅为众多步骤执行顺序中的一种方式，不代表唯一的执行顺序。在实际中的装置或终端产品执行时，可以按照实施例或者附图所示的方法顺序执行或者并行执行(例如并行处理器或者多线程处理的环境，甚至为分布式数据处理环境)。术语“包括”、“包含”或者其任何其他变体意在涵盖非排他性的包含，从而使得包括一系列要素的过程、方法、产品或者设备不仅包括那些要素，而且还包括没有明确列出的其他要素，或者是还包括为这种过程、方法、产品或者设备所固有的要素。在没有更多限制的情况下，并不排除在包括所述要素的过程、方法、产品或者设备中还存在另外的相同或等同要素。第一，第二等词语用来表示名称，而并不表示任何特定的顺序。

[0215] 为了描述的方便，描述以上装置时以功能分为各种模块分别描述。当然，在实施本说明书一个或多个时可以把各模块的功能在同一个或多个软件和/或硬件中实现，也可以将实现同一功能的模块由多个子模块或子单元的组合实现等。以上所描述的装置实施例仅仅是示意性的，例如，所述单元的划分，仅仅为一种逻辑功能划分，实际实现时可以有另外的划分方式，例如多个单元或组件可以结合或者可以集成到另一个系统，或一些特征可以忽略，或不执行。另一点，所显示或讨论的相互之间的耦合或直接耦合或通信连接可以是通过一些接口，装置或单元的间接耦合或通信连接，可以是电性，机械或其它的形式。

[0216] 本发明是参照根据本发明实施例的方法、装置(系统)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。

[0217] 这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。

[0218] 这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。

[0219] 在一个典型的配置中，计算设备包括一个或多个处理器(CPU)、输入/输出接口、网格接口和内存。

[0220] 内存可能包括计算机可读介质中的非永久性存储器，随机存取存储器(RAM)和/或非易失性内存等形式，如只读存储器(ROM)或闪存(flash RAM)。内存是计算机可读介质的示例。

[0221] 计算机可读介质包括永久性和非永久性、可移动和非可移动媒体可以由任何方法或技术来实现信息存储。信息可以是计算机可读指令、数据结构、程序的模块或其他数据。计算机的存储介质的例子包括，但不限于相变内存(PRAM)、静态随机存取存储器(SRAM)、动态随机存取存储器(DRAM)、其他类型的随机存取存储器(RAM)、只读存储器(ROM)、电可擦除可编程只读存储器(EEPROM)、快闪记忆体或其他内存技术、只读光盘只读存储器(CD-ROM)、数字多功能光盘(DVD)或其他光学存储、磁盒式磁带，磁带磁磁盘存储、石墨烯存储或其他磁性存储设备或任何其他非传输介质，可用于存储可以被计算设备访问的信息。按照本文中的界定，计算机可读介质不包括暂存电脑可读媒体(transitory media)，如调制的数据信号和载波。

[0222] 本领域技术人员应明白，本说明书一个或多个实施例可提供为方法、系统或计算机程序产品。因此，本说明书一个或多个实施例可采用完全硬件实施例、完全软件实施例或结合软件和硬件方面的实施例的形式。而且，本说明书一个或多个实施例可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、CD-ROM、光学存储器等)上实施的计算机程序产品的形式。

[0223] 本说明书一个或多个实施例可以在由计算机执行的计算机可执行指令的一般上下文中描述，例如程序模块。一般地，程序模块包括执行特定任务或实现特定抽象数据类型的例程、程序、对象、组件、数据结构等等。也可以在分布式计算环境中实践本本说明书一个或多个实施例，在这些分布式计算环境中，由通过通信网格而被连接的远程处理设备来执行任务。在分布式计算环境中，程序模块可以位于包括存储设备在内的本地和远程计算机存储介质中。

[0224] 本说明书中的各个实施例均采用递进的方式描述，各个实施例之间相同相似的部分互相参见即可，每个实施例重点说明的都是与其他实施例的不同之处。尤其，对于系统实施例而言，由于其基本相似于方法实施例，所以描述的比较简单，相关之处参见方法实施例的部分说明即可。在本说明书的描述中，参考术语“一个实施例”、“一些实施例”、“示例”、“具体示例”、或“一些示例”等的描述意指结合该实施例或示例描述的具体特征、结构、材料或者特点包含于本说明书的至少一个实施例或示例中。在本说明书中，对上述术语的示意性表述不必须针对的是相同的实施例或示例。而且，描述的具体特征、结构、材料或者特点可以在任一个或多个实施例或示例中以合适的方式结合。此外，在不相互矛盾的情况下，本领域的技术人员可以将本说明书中描述的不同实施例或示例以及不同实施例或示例的特征进行结合和组合。

[0225] 以上所述仅为本说明书一个或多个实施例的实施例而已，并不用于限制本本说明书一个或多个实施例。对于本领域技术人员来说，本说明书一个或多个实施例可以有各种更改和变化。凡在本申请的精神和原理之内所作的任何修改、等同替换、改进等，均应包含在权利要求范围之内。