CN109239776B - 一种地震波传播正演模拟方法和装置 - Google Patents

Abstract

本发明实施例提供一种地震波传播正演模拟方法和装置，该方法包括基于纵向坐标变换算法将用于地震波模拟的预设地质模型上的不规则区域从第一坐标系转换为可用于差分计算的第二坐标系下的规则区域；基于交错网格高阶差分算法求解在第二坐标系下的规则区域中的每个网格点上的物理量；根据所述第一坐标系与第二坐标系之间的映射关系，将各所述网格点上的物理量转换为第一坐标系下的物理量并作为波动方程正演模拟的输出。本发明能够实现对起伏地表等地表模型的地震波模拟，尤其是对于正演的初至波，利用本发明可以进行静校正效果的分析研究，且波形保持良好，模拟精度高，数值频散很小。

Description

一种地震波传播正演模拟方法和装置

技术领域

本发明涉及勘探地球物理学领域，具体而言，涉及一种地震波传播正演模拟方法和装置。

背景技术

地形起伏对地震波传播影响非常大，使得地震波在近地表的传播变得异常复杂，造成地震波振动强度变化以及地震波型之间的相互转化，引起面波、散射体波、多重散射面波之间的共振耦合。此外，天然地震检测中区域性地震波能量异常，以及山地油气勘探中地震数据低信噪比和“静校正”问题均是由地形起伏以及近地表速度异常造成的。现有的在地形起伏情况下的地震波传播数值模拟方法主要有以下几种：

有限元(FE)模拟方法：该模拟方法计算效率低，但处理不规则地表比较方便，成为模拟起伏地表情况下地震波传播的有效方法。为此，发展了一些有限元和其它方法相结合的混合方法，例如，Moczo等用离散波数方法模拟震源激发和下部介质中地震波的传播时，通过有限元方法来模拟沿起伏地表波的传播(DW-FE)；黄自萍等用有限元和有限差分(FE-FD)相结合的方法来模拟起伏地表地震波传播，但前述的各方法中存在的问题是：FE模拟方法和其它方法结合计算区域交界处易产生人为反射，吸收边界也不易处理。

边界积分或边界元(BE)方法：该方法将散射波场通过地表的一个半解析的积分来表示，其中积分项中Green函数一般在频率波数域中计算。这种方法在研究起伏地表地震波传播时使用比较多，但BE方法的半解析性质决定了该方法不能适用于地表速度变化较大情况，而实际应用中，由于后期地质作用造成浅部地层速度变化更为剧烈，因此限制了BE方法的实际应用。

有限差分(FD)方法：该方法计算效率高，在模拟复杂模型中地震波传播时应用最为广泛，但该方法的一个重要缺陷是处理复杂地形比较困难。

因此，研究起伏地表情况下地震波传播规律，研究起伏地表及复杂近地表结构对地震波造成的延迟，并消除这些延迟对地下深层成像造成的影响，使地震波场能够更好更准确的成像，恢复真实的地下构造特征，对于本领域技术人员而言十分迫切。

发明内容

有鉴于此，本发明提供一种地震波传播正演模拟方法和装置，以解决上述问题。

一方面，本发明较佳实施例提供一种地震波传播正演模拟方法，所述方法包括：

基于纵向坐标变换算法将具有起伏地表且用于地震波模拟的预设地质模型上的不规则区域从第一坐标系转换为可用于差分计算的第二坐标系下的规则区域；

基于交错网格高阶差分算法求解在第二坐标系下的规则区域中的每个网格点上的物理量；

根据所述第一坐标系与第二坐标系之间的映射关系，将各所述网格点上的物理量转换为第一坐标系下的物理量，并作为波动方程正演模拟的输出。

在本方面较佳实施例的选择中，所述第一坐标系为xoz坐标系，所述第二坐标系为ξoη坐标系，基于纵向坐标变换算法将用于地震波模拟的预设地质模型上的不规则区域从第一坐标系转换为可用于差分计算的第二坐标系下的规则区域的步骤包括：

将用于表征所述预设地质模型上的不规则区域的地表高程变化函数z₀＝f(x)由xoz坐标系线性变换为ξoη坐标系下的

其中，z为高程值；

将xoz坐标系下的变密度声波方程

转换为ξoη坐标系下的声波方程

其中，

K为介质体积模量，ρ为介质密度，U'＝(u',p',q')为ξoη坐标系下的列向量，

和

为坐标拉伸系数，

为x＝ξ处的地形坡度。

在本方面较佳实施例的选择中，基于交错网格高阶差分算法求解在第二坐标系下的规则区域中的每个网格点上的物理量的步骤包括：

对所述第二坐标系下的规则区域中的二维声波方程进行降维处理，以将二维声波方程中的位移对时间的任意奇数阶高阶导数转换为中间变量对空间的导数，以及将中间变量对时间的任意奇数阶高阶导数转换为位移对空间的导数；

基于交错网格高阶差分算法求解降维处理后的声波方程中的空间导数，进而确定差分系数；

基于所述差分系数计算第二坐标系下的规则区域中的每个网格点上的物理量。

在本方面较佳实施例的选择中，所述地震波包括直达波、折射波中的一种或多种。

在本方面较佳实施例的选择中，所述方法还包括:

确定预设地质模型上用于地震波传播模拟的吸收边界条件,并基于该吸收边界条件划定的区域执行上述基于纵向坐标变换算法，将第一坐标系下的预设地质模型上的不规则区域转换为可用于差分计算的第二坐标系下的规则区域的步骤。

在本方面较佳实施例的选择中，确定预设地质模型上用于地震波传播模拟的吸收边界条件的步骤包括：

在第一坐标系下确定用于地震波传播模拟的所述预设地质模型的边界条件；

基于所述纵向坐标变换算法，将所述第一坐标系下的边界条件转换为第二坐标系下的边界条件。

在本方面较佳实施例的选择中，所述第一坐标系为xoz坐标系，所述第二坐标系为ξoη坐标系，所述第二坐标系下的所述吸收边界条件包括：

其中，U'＝(u',p',q')为ξoη坐标系下的列向量，

和

为坐标拉伸系数，

为x＝ξ处的地形坡度。

另一方面，本发明较佳实施例还提供一种基于波动方程的地震波传播正演模拟装置，所述装置包括：

区域拉伸模块，用于基于纵向坐标变换算法将具有起伏地表且用于地震波模拟的预设地质模型上的不规则区域从第一坐标系转换为可用于差分计算的第二坐标系下的规则区域；

物理量计算模块，用于基于交错网格高阶差分算法求解在第二坐标系下的规则区域中的每个网格点上的物理量；

正演数据输出模块，用于根据所述第一坐标系与第二坐标系之间的映射关系，将各所述网格点上的物理量转换为第一坐标系下的物理量并作为波动方程正演模拟的输出。

在本方面较佳实施例的选择中，所述第一坐标系为xoz坐标系，所述第二坐标系为ξoη坐标系，所述区域拉伸模块包括：

线性变换单元，用于将用于表征所述预设地质模型上的不规则区域的地表高程变化函数z₀＝f(x)由xoz坐标系线性变换为ξoη坐标系下的

其中，z为高程值；

方程变换单元，用于将xoz坐标系下的变密度声波方程

转换为ξoη坐标系下的声波方程

其中，

K为介质体积模量，ρ为介质密度，U'＝(u',p',q')为ξoη坐标系下的列向量，

和

为坐标拉伸系数，

为x＝ξ处的地形坡度。

在本方面较佳实施例的选择中，所述物理量计算模块包括：

降维处理单元，用于对所述第二坐标系下的规则区域中的二维声波方程进行降维处理，以将二维声波方程中的位移对时间的任意奇数阶高阶导数转换为中间变量对空间的导数，以及将中间变量对时间的任意奇数阶高阶导数转换为位移对空间的导数；

系数确定单元，用于基于交错网格高阶差分算法求解降维处理后的声波方程中的空间导数，进而确定差分系数；

物理量计算单元，用于基于所述差分系数计算第二坐标系下的规则区域中的每个网格点上的物理量。

与现有技术相比，本发明实施例提供一种地震波传播正演模拟方法和装置，其中，通过纵向坐标变换算法对具有起伏地表的地质模型进行规则区域变换，并基于变换后的规则区域进行地震波模拟，从而实现对起伏地表等地表模型的地震波模拟，且波形保持良好，模拟精度高，数值频散很小，能够为检验静校正方法提供理论模型依据。

另外，本发明还通过吸收边界条件的设置，有效阻止了正演模拟中存在的复杂地表产生的多次波及其各种复杂散射，进一步提高了地震波模拟过程中的直达波以及折射波等初至波的模拟精度。

附图说明

为了更清楚地说明本发明实施例的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，应当理解，以下附图仅示出了本发明的某些实施例，因此不应被看作是对范围的限定，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他相关的附图。

图1为本发明实施例提供的地震波传播正演模拟方法和装置的应用场景示意图。

图2为本发明实施例提供的地震波传播正演模拟方法的流程示意图。

图3(a)为本发明实施例给出的预设地质模型的剖面结构示意图.

图3(b)为本发明实施例给出的预设地质模型的地震传播速度的设计结果示意图。

图4为图2中所示的步骤S11的子流程示意图。

图5为纵向拉伸变换结果示意图。

图6为图2中所示的步骤S12的子流程示意图。

图7为精度为ο(Δt⁴+Δx¹⁰)的交错差分网格示意图。

图8为本发明实施例提供的地震波传播正演模拟方法的另一流程示意图。

图9为基于本发明实施例提供的地震波传播正演模拟方法和预设地质模型实现的有限差分数值模拟原始炮记录结果示意图。

图10为本发明实施例提供的地震波传播正演模拟装置的方框结构示意图。

图标：10-电子终端；100-地震波传播正演模拟装置；110-区域拉伸模块；111-线性变换单元；112-方程变换单元；120-物理量计算模块；121-降维处理单元；122-系数确定单元；123-物理量计算单元；130-正演数据输出模块；200-存储器；300-存储控制器；400-处理器。

具体实施方式

为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。通常在此处附图中描述和示出的本发明实施例的组件可以以各种不同的配置来布置和设计。

因此，以下对在附图中提供的本发明的实施例的详细描述并非旨在限制要求保护的本发明的范围，而是仅仅表示本发明的选定实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

应注意到：相似的标号和字母在下面的附图中表示类似项，因此，一旦某一项在一个附图中被定义，则在随后的附图中不需要对其进行进一步定义和解释。

如图1所示，为本发明实施例提供的地震波传播正演模拟装置100的应用场景示意图。其中，电子终端10包括地震波传播正演模拟装置100、存储器200、存储控制器300以及处理器400。其中，所述电子终端10可以是，但不限于，电脑、移动上网设备(mobile Internetdevice，MID)等具有处理功能的电子设备，还可以是服务器等。

可选地，所述存储器200、存储控制器300、处理器400各元件相互之间直接或间接地电性连接，以实现数据的传输或交互。例如，这些元件之间通过一条或多条通讯总线或信号线实现电性连接。所述地震波传播正演模拟装置100包括至少一个可以软件或固件的形式存储于所述存储器200中或固化在所述电子终端10的操作系统中的软件功能模块。所述处理器400在所述存储控制器300的控制下访问所述存储器200，以用于执行所述存储器200中存储的可执行模块，例如所述地震波传播正演模拟装置100所包括的软件功能模块及计算机程序等。

可以理解，图1所示的结构仅为示意，所述电子终端10还可包括比图1中所示更多或者更少的组件，或者具有与图1所示不同的配置。图1中所示的各组件可以采用硬件、软件或其组合实现。

进一步地，请结合参阅图2，本发明实施例还提供一种可应用于所述地震波传播正演模拟装置100的地震波传播正演模拟方法。所应说明的是，本发明所述的地震波传播正演模拟方法并不以图2以及以下所述的具体顺序为限制。应当理解，本发明所述的地震波传播正演模拟方法其中部分步骤的顺序可以根据实际需要相互交换，或者其中的部分步骤也可以省略或删除。

步骤S11，基于纵向坐标变换算法将具有起伏地表且用于地震波模拟的预设地质模型上的不规则区域从第一坐标系转换为可用于差分计算的第二坐标系下的规则区域；

步骤S12，基于交错网格高阶差分算法求解在第二坐标系下的规则区域中的每个网格点上的物理量；

步骤S13，根据所述第一坐标系与第二坐标系之间的映射关系，将各所述网格点上的物理量转换为第一坐标系下的物理量并作为波动方程正演模拟的输出。

本发明实施例中，通过上述步骤S11-步骤S13中给出的地震波传播正演模拟方法，能够对具有起伏地表的地质模型进行规则区域变换，并基于变换后的规则区域进行地震波模拟，从而实现对起伏地表等地表模型的地震波模拟，且波形保持良好，模拟精度高，数值频散很小，为检验静校正方法提供了理论模型依据。

详细地，在步骤S11中，所述预设地质模型可以是但不限于准噶尔盆地的地质模型等，如图3所示为按照典型的近地表结构特征设计的地质剖面，该地质剖面中包含典型近地表结构特征的低降速带结构。

可选地，所述预设地质模型可以是如图3所示的剖面总长度为40Km、剖面深度6Km的剖面模型，其中，纵坐标为模型深度，深度1000m设置一水平界面，之上为低降速带介质。低降速带地质模型模拟各种典型地表条件，其中有冲沟、地表高程快速变化带、低降速带底界速度横向变化带、高速夹层带、山前高速带。另外，在所述预设地质模型中，可在深度1000m之下设置两个水平界面，使得在基准面静校正后这两个界面的成像将也为水平状态，以用于检验不同静校正方法中长波长的误差。可以理解的是，所述预设地质模型的实际设计结构可根据实际需求进行灵活选取，如随着深度的增加又设立了复杂构造带、水平地层带，目的是不同基准面选取方法的研究、叠前偏移方法的基准面研究。

实际实施时，可基于所述预设地质模型，按照地震数据采集的观测方式模拟野外地震采集过程，并应用声波波动理论得到模型正演的地震数据。

作为一种实施方式，如图4所示，所述步骤S11可通过下述步骤S111-步骤S112实现，其中，所述第一坐标系为xoz坐标系，所述第二坐标系为ξoη坐标系。

步骤S111，将用于表征所述预设地质模型上的不规则区域的地表高程变化函数z₀＝f(x)由xoz坐标系线性变换为ξoη坐标系下的

其中，z为高程值；

步骤S112，将xoz坐标系下的变密度声波方程

转换为ξoη坐标系下的声波方程

其中，

K为介质体积模量，ρ为介质密度，U'＝(u',p',q')为ξoη坐标系下的列向量，

和

为坐标拉伸系数，

为x＝ξ处的地形坡度。

另外，由于起伏地表条件下复杂的地形变化，在近地表区域难以构造合适的差分方程，因此，本实施例中利用纵向坐标拉伸的途径将不规则计算区域变为适合于差分计算的矩形规则区域，即采用如上述步骤S111-步骤S112中给出的纵向拉伸变换，可将原来顶界面起伏的不规则计算区域(图5a)就变换为顶界面水平的规则计算区域(图5b)，曲网格也相应地变换为矩形网格，为进行有限差分数值计算提供了方便。

在步骤S12和步骤S13中，由于本实施例中的xoz坐标系和ξoη坐标系这两个坐标系之间的坐标变换是一个线性变换，从而决定了ξoη坐标系下的物理量U'和xoz坐标系下物理量U的是一个一一对应的映射关系。因此，只要在规则ξoη坐标系下求得每个网格点上的U'，就可以确定不规则(弯曲)坐标系xoz下每个网格点上的物理量U。换言之，可通过复合求导来数值计算xoz坐标系下的不规则(弯曲)网格点处的波场的空间导数。

另外，在ξoη坐标系下进行数值计算时，还需要进行两方面的工作：(1)根据映射关系和xoz坐标系下的预设地质模型，通过插值确定各网格点上的介质物性参数；(2)如果采用纵波或横波震源，需要根据不同震源模拟方式和映射关系对其作相应变换。

进一步地，如图6所示，下面结合步骤S121-步骤S123对步骤S12的实现过程进行介绍。

步骤S121，对所述第二坐标系下的规则区域中的二维声波方程进行降维处理，以将二维声波方程中的位移对时间的任意奇数阶高阶导数转换为中间变量对空间的导数，以及将中间变量对时间的任意奇数阶高阶导数转换为位移对空间的导数；

步骤S122，基于交错网格高阶差分算法求解降维处理后的声波方程中的空间导数，进而确定差分系数；

步骤S123，基于所述差分系数计算第二坐标系下的规则区域中的每个网格点上的物理量。

其中，在密度变化较大情况下，至少需要使用标量波动方程

来描述地震波的传播。如果直接离散计算，必然遇到物性参数(密度)对空间导数的计算问题，那么在物性变化强烈的区域，一方面可能引起计算不稳定，另一方面吸收边界的处理也不方便，对此，本发明实施例中采用交错网格高阶差分解法以解决上述问题，同时可以提高模拟的精度，减小数值频散。具体地，可首先将上述二维声波方程化为下列一阶方程组：

其中，

且矩阵A、B由介质的物性参数决定K为介质体积模量，ρ为介质密度，U＝(u,p,q)^T为位移u和引入的中间变量p、q组成的列向量。

上述式(1)中的一阶方程组可以很方便地利用交错网格差分法求解，这样就避免了对密度求导这个问题，而是将其转化为求辅助变量p、q的空间导数(值得注意的是，p、q为连续函数)。同时，为了提高差分精度，作如下高阶差分近似。

利用Taylor展开，可以得到2M阶时间精度的差分近似，如：

其中，Δt为时间步长。当m＝1时(2)式即为传统的二阶精度差分近似。

直接计算(2)式中的

要涉及过多的时间层，需要较大的内存量。为此，利用方程(1)可以完全准确地将位移对时间的任意奇数阶高阶导数转嫁到中间变量对空间的导数上去，将中间变量对时间的任意奇数阶高阶导数转嫁到位移对空间的导数上去。

如2M＝2时，式(2)变为：

同理可得方程(3.1)中其它两个方程的二阶时间精度近似：

关于更高阶(2M>3)时间差分近似方程同理可得，这里不再赘述。

为了提高空间差分精度，方程(2)-(5)中的一阶空间导数利用下面的交错网格技术来计算：

关于(6)式中差分系数C_n的确定方法本实施例在此不做赘述。下面给出声波方程2阶时间差分、2N阶空间差分精度的差分格式，相应的o(Δt⁴+Δx¹⁰)精度差分网格见图7。

上述声波方程交错网格(2M,2N)阶差分解法的稳定性条件：

设

根据(7)式给出几种差分精度差分方程的稳定性条件，见表1。

进一步地，为了有效阻止复杂地表产生的多次波及其各种复杂散射，提高地震波模拟过程中的直达波以及折射波等初至波的模拟精度。在本实施例中，所述地震波传播正演模拟方法还包括：确定预设地质模型上用于地震波传播模拟的吸收边界条件,并基于该吸收边界条件划定的区域执行上述基于纵向坐标变换算法，将第一坐标系下的预设地质模型上的不规则区域转换为可用于差分计算的第二坐标系下的规则区域的步骤。

表1

详细地，如图8所示，可通过步骤S140-步骤S141实现确定预设地质模型上用于地震波传播模拟的吸收边界条件这一过程，具体如下。

步骤S140，在第一坐标系下确定用于地震波传播模拟的所述预设地质模型的边界条件；

步骤S141，基于所述纵向坐标变换算法，将所述第一坐标系下的边界条件转换为第二坐标系下的边界条件。

其中，所述第一坐标系为xoz坐标系，所述第二坐标系为ξoη坐标系，在利用上述纵向坐标变换算法进行预设地质模型中的地震波模拟时，对于吸收边界条件，同样可利用式

将xoz坐标系下的吸收边界条件变换至ξoη坐标系下并进行求解，然后利用xoz坐标系与ξoη坐标系之间的映射关系将计算得到的波场值变换至xoz坐标系下。

例如，在各向同性介质中，右边界吸收边界条件为：

为了在ξoη坐标系下使用上述吸收边界条件，同样需要据

作坐标拉伸变换，得到右边界的吸收边界条件：

那么，其它计算区域吸收边界条件在ξoη坐标系下分别变换为：

左边界：

底边界：

右下角：

左上角：

右上角：

左下角：

进一步地，基于上述地震波传播正演模拟方法的描述，下面以预设地质模型为新疆起伏地表模型为例对地震波传播数值进行正演模拟，如图9所示为模拟结果示意图，从模拟单炮记录上可以发现，波形保持良好，模拟精度高，数值频散很小，直达波和折射波比较清楚。另外，从150炮开始，由于排列所在地表起伏较大，对地震波走时的影响在模拟记录上表现明显。在750-800炮范围内，由于高速层出露，浅部无折射层存在，模拟记录上几乎看不到折射波，初至拾取比较困难。

由上模拟结果可以得出：在本发明给出的地震传播正演模拟方法中，直达波和折射波比较清晰，能够为静校正方法研究提供了理论数据。同时不同的近地表速度结构对直达波和折射波影响较大，在高速地层出露地区，没有折射波存在。

进一步地，如图10所示，为本发明实施例给出的地震波传播正演模拟装置100的方框结构示意图，所述地震波传播正演模拟装置100包括区域拉伸模块110、物理量计算模块120和正演数据输出模块130。

所述区域拉伸模块110，用于基于纵向坐标变换算法将具有起伏地表且用于地震波模拟的预设地质模型上的不规则区域从第一坐标系转换为可用于差分计算的第二坐标系下的规则区域；本实施例中，步骤S11可由所述区域拉伸模块110执行，具体过程请参考步骤S11，在此不再赘述。可选地，在本实施例中，所述区域拉伸模块110包括线性变换单元111和方程变换单元112。

所述线性变换单元111，用于将用于表征所述预设地质模型上的不规则区域的地表高程变化函数z₀＝f(x)由xoz坐标系线性变换为ξoη坐标系下的

其中，z为高程值；本实施例中，步骤S111可由所述线性变换单元111执行，具体过程请参考步骤S111，在此不再赘述。

所述方程变换单元112，用于将xoz坐标系下的变密度声波方程

转换为ξoη坐标系下的声波方程

本实施例中，步骤S112可由所述方程变换单元112执行，具体过程请参考步骤S112，在此不再赘述。

所述物理量计算模块120，用于基于交错网格高阶差分算法求解在第二坐标系下的规则区域中的每个网格点上的物理量；本实施例中，步骤S12可由所述物理量计算模块120执行，具体过程请参考步骤S12，在此不再赘述。可选地，在本实施例中，所述物理量计算模块120包括降维处理单元121、系数确定单元122和物理量计算单元123。

所述降维处理单元121，用于对所述第二坐标系下的规则区域中的二维声波方程进行降维处理，以将二维声波方程中的位移对时间的任意奇数阶高阶导数转换为中间变量对空间的导数，以及将中间变量对时间的任意奇数阶高阶导数转换为位移对空间的导数；本实施例中，步骤S121可由所述降维处理单元121执行，具体过程请参考步骤S121，在此不再赘述。

所述系数确定单元122，用于基于交错网格高阶差分算法求解降维处理后的声波方程中的空间导数，进而确定差分系数；本实施例中，步骤S122可由所述系数确定单元122执行，具体过程请参考步骤S122，在此不再赘述。

所述物理量计算单元123，用于基于所述差分系数计算第二坐标系下的规则区域中的每个网格点上的物理量。本实施例中，步骤S123可由所述物理量计算单元123执行，具体过程请参考步骤S123，在此不再赘述。

所述正演数据输出模块130，用于根据所述第一坐标系与第二坐标系之间的映射关系，将各所述网格点上的物理量转换为第一坐标系下的物理量并作为波动方程正演模拟的输出。本实施例中，步骤S13可由所述正演数据输出模块130执行，具体过程请参考步骤S13，在此不再赘述。

综上所述，本发明实施例提供一种地震波传播正演模拟方法和装置，其中，通过纵向坐标变换算法对具有起伏地表的地质模型进行规则区域变换，并基于变换后的规则区域进行地震波模拟，从而实现对起伏地表等地表模型的地震波模拟，且波形保持良好，模拟精度高，数值频散很小，为检验静校正方法提供了理论模型依据。

另外，本发明还通过吸收边界调节的设置，有效阻止了复杂地表产生的多次波及其各种复杂散射，进一步提高了地震波模拟过程中的直达波以及折射波等初至波的模拟精度。

在本发明的描述中，术语“设置”、“相连”、“连接”应做广义理解，例如，可以是固定连接，也可以是可拆卸连接，或一体地连接；可以是机械连接，也可以是电连接；可以是直接相连，也可以通过中间媒介间接相连，可以是两个元件内部的连通。对于本领域的普通技术人员而言，可以具体情况理解上述术语在本发明中的具体含义。

在本发明实施例所提供的几个实施例中，应该理解到，所揭露的装置和方法，也可以通过其他方式实现。以上所描述的装置和方法实施例仅仅是示意性的，例如，附图中的流程图和框图显示了根据本发明的预设数量个实施例的装置、方法和计算机程序产品可能实现的体系架构、功能和操作。在这点上，流程图或框图中的每个方框可以代表一个模块、程序段或代码的一部分。所述模块、程序段或代码的一部分包含一个或预设数量个用于实现规定的逻辑功能。

也应当注意，在有些作为替换的实现方式中，方框中所标注的功能也可以以不同于附图中所标注的顺序发生。例如，两个连续的方框实际上可以基本并行地执行，它们有时也可以按相反的顺序执行，这依所涉及的功能而定。也要注意的是，框图和/或流程图中的每个方框、以及框图和/或流程图中的方框的组合，可以用执行规定的功能或动作的专用的基于硬件的系统来实现，或者可以用专用硬件与计算机指令的组合来实现。

以上所述仅为本发明的优选实施例而已，并不用于限制本发明，对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (8)

1.一种地震波传播正演模拟方法，其特征在于，所述方法包括：

基于纵向坐标变换算法将具有起伏地表且用于地震波模拟的预设地质模型上的不规则区域从第一坐标系转换为可用于差分计算的第二坐标系下的规则区域，具体的步骤包括：

所述第一坐标系为xoz坐标系，所述第二坐标系为ξoη坐标系；

将用于表征所述预设地质模型上的不规则区域的地表高程变化函数z₀＝f(x)由xoz坐标系线性变换为ξoη坐标系下的

其中，z为高程值；

将xoz坐标系下的变密度声波方程

转换为ξoη坐标系下的声波方程

其中，

K为介质体积模量，ρ为介质密度，U'＝(u',p',q')为ξoη坐标系下的列向量，

为坐标拉伸系数，

为x＝ξ处的地形坡度；

基于交错网格高阶差分算法求解在第二坐标系下的规则区域中的每个网格点上的物理量；

根据所述第一坐标系与第二坐标系之间的映射关系，将各所述网格点上的物理量转换为第一坐标系下的物理量，并作为波动方程正演模拟的输出。

2.根据权利要求1所述的地震波传播正演模拟方法，其特征在于，基于交错网格高阶差分算法求解在第二坐标系下的规则区域中的每个网格点上的物理量的步骤包括：

对所述第二坐标系下的规则区域中的二维声波方程进行降维处理，以将二维声波方程中的位移对时间的任意奇数阶高阶导数转换为中间变量对空间的导数，以及将中间变量对时间的任意奇数阶高阶导数转换为位移对空间的导数；

基于交错网格高阶差分算法求解降维处理后的声波方程中的空间导数，进而确定差分系数；

基于所述差分系数计算第二坐标系下的规则区域中的每个网格点上的物理量。

3.根据权利要求1所述的地震波传播正演模拟方法，其特征在于，所述地震波包括直达波、折射波中的一种或多种。

4.根据权利要求1所述的地震波传播正演模拟方法，其特征在于，所述方法还包括:

确定预设地质模型上用于地震波传播模拟的吸收边界条件,并基于该吸收边界条件划定的区域执行上述基于纵向坐标变换算法，将第一坐标系下的预设地质模型上的不规则区域转换为可用于差分计算的第二坐标系下的规则区域的步骤。

5.根据权利要求4所述的地震波传播正演模拟方法，其特征在于，确定预设地质模型上用于地震波传播模拟的吸收边界条件的步骤包括：

在第一坐标系下确定用于地震波传播模拟的所述预设地质模型的边界条件；

基于所述纵向坐标变换算法，将所述第一坐标系下的边界条件转换为第二坐标系下的边界条件。

6.根据权利要求5所述的地震波传播正演模拟方法，其特征在于，所述第一坐标系为xoz坐标系，所述第二坐标系为ξoη坐标系，所述第二坐标系下的所述吸收边界条件包括：

其中，U'＝(u',p',q')为ξoη坐标系下的列向量，

和

为坐标拉伸系数，

为x＝ξ处的地形坡度；A⁽⁺⁾与A^(-)为所述变密度声波方程的参数A经特征分解法所获得的矩阵，分别表示沿ξoη坐标系中x轴两方向的单程波；B⁽⁺⁾与B^(-)为所述变密度声波方程的参数B经特征分解法所获得的矩阵,分别表示沿ξoη坐标系中z轴两方向的单程波。

7.一种地震波传播正演模拟装置，其特征在于，所述装置包括：

区域拉伸模块，用于基于纵向坐标变换算法将具有起伏地表且用于地震波模拟的预设地质模型上的不规则区域从第一坐标系转换为可用于差分计算的第二坐标系下的规则区域,其中，所述第一坐标系为xoz坐标系，所述第二坐标系为ξoη坐标系，所述区域拉伸模块包括：

线性变换单元，用于将用于表征所述预设地质模型上的不规则区域的地表高程变化函数z₀＝f(x)由xoz坐标系线性变换为ξoη坐标系下的

其中，z为高程值；

方程变换单元，用于将xoz坐标系下的变密度声波方程

转换为ξoη坐标系下的声波方程

其中，

K为介质体积模量，ρ为介质密度，U'＝(u',p',q')为ξoη坐标系下的列向量，

和

为坐标拉伸系数，

为x＝ξ处的地形坡度；

物理量计算模块，用于基于交错网格高阶差分算法求解在第二坐标系下的规则区域中的每个网格点上的物理量；

正演数据输出单元，用于根据所述第一坐标系与第二坐标系之间的映射关系，将各所述网格点上的物理量转换为第一坐标系下的物理量并作为波动方程正演模拟的输出。

8.根据权利要求7所述的地震波传播正演模拟装置，其特征在于，所述物理量计算模块包括：

降维处理单元，用于对所述第二坐标系下的规则区域中的二维声波方程进行降维处理，以将二维声波方程中的位移对时间的任意奇数阶高阶导数转换为中间变量对空间的导数，以及将中间变量对时间的任意奇数阶高阶导数转换为位移对空间的导数；

系数确定单元，用于基于交错网格高阶差分算法求解降维处理后的声波方程中的空间导数，进而确定差分系数；

物理量计算单元，用于基于所述差分系数计算第二坐标系下的规则区域中的每个网格点上的物理量。

Images