CN111142157A - 一种三维非均匀介质弹性波的处理方法、装置和设备 - Google Patents

一种三维非均匀介质弹性波的处理方法、装置和设备 Download PDF

Info

Publication number
CN111142157A
CN111142157A CN202010021263.6A CN202010021263A CN111142157A CN 111142157 A CN111142157 A CN 111142157A CN 202010021263 A CN202010021263 A CN 202010021263A CN 111142157 A CN111142157 A CN 111142157A
Authority
CN
China
Prior art keywords
format
determining
spatial
sprk
time
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN202010021263.6A
Other languages
English (en)
Other versions
CN111142157B (zh
Inventor
王健
杨顶辉
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Tsinghua University
Original Assignee
Tsinghua University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Tsinghua University filed Critical Tsinghua University
Priority to CN202010021263.6A priority Critical patent/CN111142157B/zh
Publication of CN111142157A publication Critical patent/CN111142157A/zh
Application granted granted Critical
Publication of CN111142157B publication Critical patent/CN111142157B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00—Details of seismic processing or analysis
    • G01V2210/60—Analysis
    • G01V2210/67—Wave propagation modeling
    • G01V2210/673—Finite-element; Finite-difference
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00—Details of seismic processing or analysis
    • G01V2210/60—Analysis
    • G01V2210/67—Wave propagation modeling
    • G01V2210/675—Wave equation; Green's functions

Landscapes

  • Engineering & Computer Science (AREA)
  • Remote Sensing (AREA)
  • Physics & Mathematics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Acoustics & Sound (AREA)
  • Environmental & Geological Engineering (AREA)
  • Geology (AREA)
  • General Life Sciences & Earth Sciences (AREA)
  • General Physics & Mathematics (AREA)
  • Geophysics (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本说明书实施例公开了一种三维非均匀介质弹性波的处理方法、装置和设备，该方法包括：根据待勘探的目标区域、当前具备的计算资源和最大库朗数，确定所述目标区域对应的时间步长、空间步长和空间精度；根据所述空间精度确定有限差分格式优化系数，并根据所述目标区域对应的空间步长和所述有限差分格式优化系数，确定高阶空间偏导数的离散格式；基于所述目标区域对应的时间步长、预定系数向量、预定空间算子，确定修正保辛分部龙格库塔SPRK时间格式；基于所述修正SPRK时间格式和所述高阶空间偏导数的离散格式，确定三维非均匀介质弹性波方程的数值解。

Description

一种三维非均匀介质弹性波的处理方法、装置和设备
技术领域
本说明书涉及油气地震勘探的技术领域，尤其涉及一种三维非均匀介质弹性波的处理方法、装置和设备。
背景技术
弹性波在地球介质中的传播实际上是一个波动方程的正演过程，为了获取地球内部信息，地震波正演模拟技术一直是地球物理学中非常重要的技术。地震波正演模拟技术是利用数值方法模拟地震波在已知地下介质中传播的技术，也是油气勘探反演成像的关键技术和反演问题的基础技术。同时，地震波正演模拟对于解释地震信息、评估观测系统等也具有重要意义。
有限差分方法是目前应用最为广泛的地震波正演模拟算法。有限差分方法会将计算区域划分为有限个规则的或不规则的网格点，利用泰勒展开的方式将波动方程中的微分算子转化为网格点上函数值的差分形式，通过求解网格点上的函数值近似得到原波动方程的解。
然而，有限差分方法也面临着可能产生数值频散的问题。数值频散的产生是由于利用差分形式近似微分算子时，不同频率波的相速度产生了差别，特别是当网格步长很大或波场梯度很大时，这种现象尤为明显。因此当波传播时间较长时，不同频率的波由于相速度不同而会发生分离。这显然与实际物理定律相违背，故而产生了虚假信息，影响了正演模拟的精度。
发明内容
本说明书实施例的目的是提供一种三维非均匀介质弹性波的处理方法、装置和设备，以解决求解三维非均匀介质弹性波的数值解中存在的数值频散问题。
为解决上述技术问题，本说明书实施例是这样实现的：
本说明书实施例提供的一种三维非均匀介质弹性波的处理方法，所述方法包括：
根据待勘探的目标区域、当前具备的计算资源和最大库朗数，确定所述目标区域对应的时间步长、空间步长和空间精度；
根据所述空间精度确定有限差分格式优化系数，并根据所述目标区域对应的空间步长和所述有限差分格式优化系数，确定高阶空间偏导数的离散格式；
基于所述目标区域对应的时间步长、预定系数向量、预定空间算子，确定修正保辛分部龙格库塔SPRK时间格式；
基于所述修正SPRK时间格式和所述高阶空间偏导数的离散格式，确定三维非均匀介质弹性波方程的数值解。
可选地，所述根据所述目标区域对应的空间步长和所述有限差分格式优化系数，确定高阶空间偏导数的离散格式，包括：
根据所述目标区域对应的空间步长和所述有限差分格式优化系数，分别通过公式
Figure BDA0002360884700000021
Figure BDA0002360884700000022
Figure BDA0002360884700000023
Figure BDA0002360884700000024
Figure BDA0002360884700000025
Figure BDA0002360884700000026
Figure BDA0002360884700000027
Figure BDA0002360884700000028
Figure BDA0002360884700000029
Figure BDA0002360884700000031
Figure BDA0002360884700000032
Figure BDA0002360884700000033
确定高阶空间偏导数的离散格式，其中，c、λ、μ表示预定系数，Δx、Δy和Δz分别表示所述目标区域中x方向的空间步长、y方向的空间步长和z方向的空间步长，
Figure BDA0002360884700000034
和
Figure BDA0002360884700000035
为所述有限差分格式化系数。
可选地，所述基于所述目标区域对应的时间步长、预定系数向量、预定空间算子，确定修正保辛分部龙格库塔SPRK时间格式，包括：
基于所述目标区域对应的时间步长、预定系数向量、预定空间算子，确定以下修正保辛分部龙格库塔SPRK时间格式：
Figure BDA0002360884700000036
其中，Δt表示所述目标区域对应的时间步长，c＝(c1,c2,c3)和d＝(d1,d2)表示所述预定系数向量，L表示所述预定空间算子，un和vn表示在第n个时间步长的数值解，un+1和vn +1表示在第n+1个时间步长的数值解，u1和v1表示中间变量。
可选地，所述方法还包括：
根据所述目标区域中声波或弹性波的初始波速、所述预定空间算子对应的特征值，确定所述修正SPRK时间格式的最大库朗数，以确定所述修正SPRK时间格式的稳定性条件。
可选地，所述根据所述目标区域中声波或弹性波的初始波速、所述预定空间算子对应的特征值，确定所述修正SPRK时间格式的最大库朗数，以确定所述修正SPRK时间格式的稳定性条件，包括：
根据所述目标区域中声波或弹性波的初始波速、所述预定空间算子对应的特征值，通过下述公式
Figure BDA0002360884700000037
确定所述修正SPRK时间格式的最大库朗数，以确定所述修正SPRK时间格式的稳定性条件，其中，c0表示所述目标区域中P波的最大波速，αmax表示最大库朗数，χmin表示所述预定空间算子对应的最小特征值。
本说明书实施例提供的一种三维非均匀介质弹性波的处理装置，所述装置包括：
区域信息确定模块，用于根据待勘探的目标区域、当前具备的计算资源和最大库朗数，确定所述目标区域对应的时间步长、空间步长和空间精度；
离散格式确定模块，用于根据所述空间精度确定有限差分格式优化系数，并根据所述目标区域对应的空间步长和所述有限差分格式优化系数，确定高阶空间偏导数的离散格式；
修正格式确定模块，用于基于所述目标区域对应的时间步长、预定系数向量、预定空间算子，确定修正保辛分部龙格库塔SPRK时间格式；
求解模块，用于基于所述修正SPRK时间格式和所述高阶空间偏导数的离散格式，确定三维非均匀介质弹性波方程的数值解。
可选地，所述离散格式确定模块，用于
根据所述目标区域对应的空间步长和所述有限差分格式优化系数，分别通过公式
Figure BDA0002360884700000041
Figure BDA0002360884700000042
Figure BDA0002360884700000043
Figure BDA0002360884700000044
Figure BDA0002360884700000045
Figure BDA0002360884700000046
Figure BDA0002360884700000051
Figure BDA0002360884700000052
Figure BDA0002360884700000053
Figure BDA0002360884700000054
Figure BDA0002360884700000055
Figure BDA0002360884700000056
确定高阶空间偏导数的离散格式，其中，c、λ、μ表示预定系数，Δx、Δy和Δz分别表示所述目标区域中x方向的空间步长、y方向的空间步长和z方向的空间步长，
Figure BDA0002360884700000057
和
Figure BDA0002360884700000058
为所述有限差分格式化系数。
可选地，所述修正格式确定模块，用于基于所述目标区域对应的时间步长、预定系数向量、预定空间算子，确定以下修正保辛分部龙格库塔SPRK时间格式：
Figure BDA0002360884700000059
其中，Δt表示所述目标区域对应的时间步长，c＝(c1,c2,c3)和d＝(d1,d2)表示所述预定系数向量，L表示所述预定空间算子，un和vn表示在第n个时间步长的数值解，un+1和vn +1表示在第n+1个时间步长的数值解，u1和v1表示中间变量。
可选地，所述装置还包括：
稳定性条件确定模块，用于根据所述目标区域中声波或弹性波的初始波速、所述预定空间算子对应的特征值，确定所述修正SPRK时间格式的最大库朗数，以确定所述修正SPRK时间格式的稳定性条件。
本说明书实施例提供的一种三维非均匀介质弹性波的处理设备，所述三维非均匀介质弹性波的处理设备包括：
处理器；以及
被安排成存储计算机可执行指令的存储器，所述可执行指令在被执行时使所述处理器：
根据待勘探的目标区域、当前具备的计算资源和最大库朗数，确定所述目标区域对应的时间步长、空间步长和空间精度；
根据所述空间精度确定有限差分格式优化系数，并根据所述目标区域对应的空间步长和所述有限差分格式优化系数，确定高阶空间偏导数的离散格式；
基于所述目标区域对应的时间步长、预定系数向量、预定空间算子，确定修正保辛分部龙格库塔SPRK时间格式；
基于所述修正SPRK时间格式和所述高阶空间偏导数的离散格式，确定三维非均匀介质弹性波方程的数值解。
由以上本说明书实施例提供的技术方案可见，本说明书实施例通过根据待勘探的目标区域、当前具备的计算资源和预定的最大库朗数，确定目标区域对应的时间步长、空间步长和空间精度，根据空间精度确定有限差分格式优化系数，并根据目标区域对应的空间步长和有限差分格式优化系数，确定高阶空间偏导数的离散格式，基于目标区域对应的时间步长、预定系数向量、预定空间算子，确定修正保辛分部龙格库塔SPRK时间格式，基于修正SPRK时间格式和高阶空间偏导数的离散格式，确定三维非均匀介质弹性波方程的数值解，这样，通过优化函数，最终得到最优的有限差分格式优化系数，从而可以大大节省内存需求和计算成本，而且，修正SPRK格式通过引入修正项，可以使一个二级格式达到了三阶时间精度。
附图说明
为了更清楚地说明本说明书实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本说明书中记载的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。
图1为本说明书一种三维非均匀介质弹性波的处理方法实施例；
图2为本说明书另一种三维非均匀介质弹性波的处理方法实施例；
图3为本说明书一种在α＝0.1,Sp＝0.4时不同空间精度的MTSOS方法在不同方向上的频散误差的示意图；
图4为本说明书一种在α＝0.1,Sp＝0.4时不同空间精度的SPRK方法在不同方向上的频散误差的示意图；
图5为本说明书一种t＝3.0s时刻位移场不同分量分别在震源所在XY、XZ和YZ平面的波场快照的示意图；
图6为本说明书一种在0-5s内位移三分量波形与解析解的比较的示意图；
图7为本说明书一种t＝1.75s时刻的位移场的不同分量分别在震源所在的XY、XZ和YZ平面的波场快照示意图；
图8为本说明书一种t＝1.75s时位移场x分量在XZ平面的波场快照示意图；
图9为本说明书一种以位移场的x分量在XZ平面为例给出了从t＝0.5s至t＝2.5s时刻的波场快照示意图；
图10为本说明书一种在t＝0.25s时刻位移场不同分量于震源所在的XY、XZ和YZ平面的波场快照示意图；
图11为本说明书一种波形示意图；
图12为本说明书一种波形示意图；
图13为本说明书一种三维非均匀介质弹性波的处理装置实施例；
图14为本说明书一种三维非均匀介质弹性波的处理设备实施例。
具体实施方式
本说明书实施例提供一种三维非均匀介质弹性波的处理方法、装置和设备。
为了使本技术领域的人员更好地理解本说明书中的技术方案，下面将结合本说明书实施例中的附图，对本说明书实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本说明书一部分实施例，而不是全部的实施例。基于本说明书中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都应当属于本文件保护的范围。
如图1所示，本说明书实施例提供一种三维非均匀介质弹性波的处理方法，该方法的执行主体可以为终端设备或服务器，其中，该终端设备可以如个人计算机等终端设备，还可以是如手机或平板电脑等移动终端设备。服务器可以是一个独立的服务器，还可以是由多个服务器构成的服务器集群。为了提高处理效率，本说明书实施例中的执行主体以服务器为例进行说明，对于执行主体为终端设备的情况，可以参见下述服务器的相关内容进行处理，在此不再赘述。该方法具体可以包括以下步骤：
在步骤S102中，根据待勘探的目标区域、当前具备的计算资源和最大库朗数，确定目标区域对应的时间步长、空间步长和空间精度。
其中，目标区域可以是需要进行地球物理勘探(如油气勘探或地震波探测等)的任意区域。计算资源可以包括服务器的CPU资源、GPU资源、存储资源等。库朗数可以是指时间步长和空间步长的相对关系，通常，随着库朗数从小到大的变化，数值模拟效率逐渐加快，但是稳定性逐渐降低。
在实施中，弹性波在地球介质中的传播实际上是一个波动方程的正演过程，为了获取地球内部信息，地震波正演模拟技术一直是地球物理学中非常重要的技术。地震波正演模拟技术是利用数值方法模拟地震波在已知地下介质中传播的技术，也是油气勘探反演成像的关键技术和反演问题的基础技术。同时，地震波正演模拟对于解释地震信息、评估观测系统等也具有重要意义。有限差分方法是目前应用最为广泛的地震波正演模拟算法。有限差分方法会将计算区域划分为有限个规则的或不规则的网格点，利用泰勒展开的方式将波动方程中的微分算子转化为网格点上函数值的差分形式，通过求解网格点上的函数值近似得到原波动方程的解。此外，不同的差分形式对应不同的有限差分方法。有限差分方法具有易于实现、存储量小、运算速度快、易于并行等优点。
然而，有限差分方法也面临着可能产生数值频散的问题。数值频散的产生是由于利用差分形式近似微分算子时，不同频率波的相速度产生了差别，特别是当网格步长很大或波场梯度很大时，这种现象尤为明显。因此当波传播时间较长时，不同频率的波由于相速度不同而会发生分离。这显然与现实的物理定律相违背，故而产生了虚假信息，影响了正演模拟的精度。
为解决上述问题，可以构造优化有限差分格式，将数值波数和真实波数之差作为优化函数，从而得到最优的差分离散系数。其中，可以将数值波数和真实波数的无穷范数误差作为目标函数，通过模拟退火法优化目标函数，最终得到最优的差分离散系数。数值结果表明，使用优化的高阶有限差分方法，可以大大节省内存需求和计算成本。
波动方程可以被看作是一个无限维线性可分的哈密尔顿系统(HamiltonSystem)，其相空间存在着辛几何结构，可以表示物体的运动规律。在哈密尔顿体系下构造的辛算法(Symplectic Method)可以保持哈密尔顿系统的辛几何结构不变，因此具有重要意义。修正保辛分部龙格库塔(Symplectic Partitioned Runge-Kutta,SPRK)格式通过引入修正项，可以使一个二级格式达到了三阶时间精度。
基于上述内容，当需要对某一区域(即目标区域)进行地球物理勘探时，可以确定目标区域，可以包括获取目标区域的坐标位置、目标区域周围的环境信息等，可以将获取的上述信息输入到服务器中，服务器可以根据输入的数据构建目标区域所在的地理模型。可以统计当前所具备的计算资源，具体地，服务器可以获取自身和其它需要使用的服务器的资源信息，可以将获取的资源信息对应的资源作为当前所具备的计算资源。此外，还可以根据实际情况设定适当的库朗数(本说明书实施例中可以使用最大库朗数)。服务器得到上述信息后，可以针对目标区域的范围、计算资源的量和计算资源对应的数据计算能力，选择合适的库朗数(即最大库朗数)，进而可以确定目标区域对应的时间步长、空间步长和空间精度等。
在步骤S104中，根据上述空间精度确定有限差分格式优化系数，并根据目标区域对应的空间步长和有限差分格式优化系数，确定高阶空间偏导数的离散格式。
其中，有限差分格式优化系数可以是一种求偏微分(或常微分)方程和方程组定解问题的数值解中所使用的优化系数。
在实施中，为了对目标区域内的三维非均匀介质弹性波方程进行求解，服务器通过上述方式确定目标区域对应的空间精度后，可以基于需要使用的有限差分离散方法和空间精度等确定有限差分格式优化系数。然后，服务器可以建立针对目标区域的弹性波方程，并可以将目标区域对应的空间步长和得到的有限差分格式优化系数分别代入到上述弹性波方程中，从而得到相应的高阶空间偏导数的离散格式。
需要说明的是，可以利用不同的有限差分格式优化系数构建不同的有限差分离散方法，例如可以采用如下的有限差分离散方法，即利用u(x0+jΔx)在x＝x0处泰勒级数展开，以保证相应的有限差分格式具有预定的空间精度，并可以通过求解线性方程组以确定有限差分格式系数。
在步骤S106中，基于目标区域对应的时间步长、预定系数向量、预定空间算子，确定修正保辛分部龙格库塔SPRK时间格式。
其中，预定系数向量可以通过泰勒级数展开等方式确定。预定空间算子可以是为了使得二级格式具有三阶时间精度而引入的修正信息。
在实施中，对于2n维哈密尔顿系统，如果存在两个关于广义动量和广义坐标的线性函数，则上述哈密尔顿系统即为线性可分的哈密尔顿系统，其中，线性可分的哈密尔顿系统中包括目标区域对应的时间步长和预定空间算子，为此，服务器可以将弹性波方程表示为线性可分的哈密尔顿系统，进而可以得到显式SPRK时间格式，其中，该显式SPRK时间格式中还包括预定系数向量等。为了保证二级格式具有三阶时间精度，可以引入相应的修正项，进而可以得到修正保辛分部龙格库塔SPRK时间格式。
在步骤S108中，基于修正SPRK时间格式和高阶空间偏导数的离散格式，确定三维非均匀介质弹性波方程的数值解。
在实施中，服务器可以将修正SPRK时间格式与高阶空间偏导数的离散格式相结合，将空间格式拓展到三维情况，同时利用等效介质参数，边界修正法和完全匹配层吸收边界条件，即可获得了三维非均匀介质弹性波方程正演模拟的修正时空优化保辛(MTSOS)方法，并可以通过上述修正时空优化保辛(MTSOS)方法计算三维非均匀介质弹性波方程的数值解，从而得到三维非均匀介质弹性波方程的数值解。
本说明书实施例提供一种三维非均匀介质弹性波的处理方法，通过根据待勘探的目标区域、当前具备的计算资源和预定的库朗数，确定目标区域对应的时间步长、空间步长和空间精度，根据空间精度确定有限差分格式优化系数，并根据目标区域对应的空间步长和有限差分格式优化系数，确定高阶空间偏导数的离散格式，基于目标区域对应的时间步长、预定系数向量、预定空间算子，确定修正保辛分部龙格库塔SPRK时间格式，基于修正SPRK时间格式和高阶空间偏导数的离散格式，确定三维非均匀介质弹性波方程的数值解，这样，通过优化函数，最终得到最优的有限差分格式优化系数，从而可以大大节省内存需求和计算成本，而且，修正SPRK格式通过引入修正项，可以使一个二级格式达到了三阶时间精度。
实施例二
如图2所示，本说明书实施例提供一种三维非均匀介质弹性波的处理方法，该方法的执行主体可以为终端设备或服务器，其中，该终端设备可以如个人计算机等终端设备，还可以是如手机或平板电脑等移动终端设备。服务器可以是一个独立的服务器，还可以是由多个服务器构成的服务器集群。为了提高处理效率，本说明书实施例中的执行主体以服务器为了进行说明，对于执行主体为终端设备的情况，可以参见下述服务器的相关内容进行处理，在此不再赘述。该方法具体可以包括以下步骤：
在步骤S202中，根据待勘探的目标区域、当前具备的计算资源和最大库朗数，确定目标区域对应的时间步长、空间步长和空间精度。
在步骤S204中，根据上述空间精度确定有限差分格式优化系数，并根据目标区域对应的空间步长和有限差分格式优化系数，分别通过公式
Figure BDA0002360884700000101
Figure BDA0002360884700000102
Figure BDA0002360884700000103
Figure BDA0002360884700000111
Figure BDA0002360884700000112
Figure BDA0002360884700000113
Figure BDA0002360884700000114
Figure BDA0002360884700000115
Figure BDA0002360884700000116
Figure BDA0002360884700000117
Figure BDA0002360884700000118
确定高阶空间偏导数的离散格式，其中，c、λ、μ表示预定系数，Δx、Δy和Δz分别表示所述目标区域中x方向的空间步长、y方向的空间步长和z方向的空间步长，
Figure BDA0002360884700000119
和
Figure BDA00023608847000001110
为所述有限差分格式化系数
在实施中，为了得到弹性波方程的数值解，需要对空间算子L进行离散。有限差分方法是一种易于编程及并行的数值离散方法。有限差分格式对函数的空间高阶偏导数，利用该点及周围点的函数值进行离散，即
Figure BDA00023608847000001111
其中Δx表示离散网格的空间步长。可利用不同的有限差分格式优化系数
Figure BDA00023608847000001112
发展不同的有限差分离散方法。通常，有限差分离散方法可以利用u(x0+jΔx)在x＝x0处的泰勒级数展开，以保证相应的格式具有预先给定的空间精度，通过求解线性方程组确定有限差分格式优化系数
Figure BDA0002360884700000121
但上述方法在粗网格下可能会产生强烈的数值频散，如果利用减小空间步长的方法来压制数值频散，又会使计算费用剧增。为了解决上述问题，可以将数值波数与真实波数的无穷范数误差作为目标函数。对于目标函数的一阶空间偏导数，定义误差函数为：
Figure BDA0002360884700000122
其中，kc是给定误差限后的最大真实波数。通过模拟退火法(SimulatedAnnealing Algorithm)，可以得到最优有限差分格式优化系数。对于二阶空间偏导数，类似地定义目标函数为：
Figure BDA0002360884700000123
通过模拟退火法，可以得到最优有限差分格式优化系数。这样，可以得到公式(4)所需的有限差分格式优化系数。使用如上定义的有限差分格式优化系数所形成的有限差分格式，可以减少经典高阶显式有限差分格式的数值频散。
其中，不同空间精度的一阶空间偏导数的有限差分格式优化系数
Figure BDA0002360884700000124
可以如下表1所示
表1
Figure BDA0002360884700000125
不同空间精度的二阶空间偏导数的有限差分格式优化系数
Figure BDA0002360884700000126
可以如下表2所示
表2
Figure BDA0002360884700000131
在非均匀介质中，需要考虑介质的速度或弹性参数，对于弹性波方程，还要考虑混合偏导数。为此可以利用二维非均匀介质中弹性波方程空间算子离散的有限差分优化格式。对于声波方程，以
Figure BDA0002360884700000132
为例，其离散格式可定义为：
Figure BDA0002360884700000133
其中，
Figure BDA0002360884700000134
和
Figure BDA0002360884700000135
分别表示在
Figure BDA0002360884700000136
和
Figure BDA0002360884700000137
位置处的波速，ui,j,k表示在(i,j,k)处的位移。
对于弹性波方程
Figure BDA0002360884700000141
其中，λ和μ分别表示拉梅系数，u＝(u1,u2,u3)T表示位移场。此时，需要处理单向二阶空间偏导数和混合二阶空间偏导数。对于单向二阶空间偏导数，以
Figure BDA0002360884700000142
为例，弹性波方程对应的离散格式可定义为：
Figure BDA0002360884700000143
其中，
Figure BDA0002360884700000144
和
Figure BDA0002360884700000145
分别表示在
Figure BDA0002360884700000146
和
Figure BDA0002360884700000147
位置处的拉梅系数。
对于混合二阶空间偏导数，以
Figure BDA0002360884700000148
为例，弹性波方程对应的离散格式可定义为：
Figure BDA0002360884700000149
对于其它的二阶空间偏导数，可以通过上述类似的方式得到相应的离散格式，进而可以得到高阶空间偏导数的离散格式。
在步骤S206中，基于目标区域对应的时间步长、预定系数向量、预定空间算子，确定以下修正保辛分部龙格库塔SPRK时间格式：
Figure BDA0002360884700000151
其中，Δt表示目标区域对应的时间步长，c＝(c1,c2,c3)和d＝(d1,d2)表示预定系数向量，L表示预定空间算子，un和vn表示在第n个时间步长的数值解，un+1和vn+1表示在第n+1个时间步长的数值解，u1和v1表示中间变量。
在实施中，可以将2n维哈密尔顿系统表示为如下公式
Figure BDA0002360884700000152
其中，
Figure BDA0002360884700000153
和
Figure BDA0002360884700000154
分别表示广义动量和广义坐标，H为系统的哈密尔顿函数。如果存在
Figure BDA0002360884700000155
且f和g分别是关于u和v的线性函数，则称上述公式(9)对应的哈密尔顿系统为线性可分的哈密尔顿系统。对于弹性波方程(7)，可将上述弹性波方程表示为线性可分的哈密尔顿系统，即
Figure BDA0002360884700000156
其中，
Figure BDA0002360884700000157
Figure BDA0002360884700000158
Figure BDA0002360884700000159
空间算子L的其它项也可以通过上述方式可得，s阶显式SPRK格式可表示为：
Figure BDA0002360884700000161
其中，un和vn表示该哈密尔顿系统在第n个时间步长的数值解，c＝(c1,c2,...,cs)和d＝(d1,d2,...,ds)是该s阶显式SPRK格式的系数向量，Δt表示该s阶显式SPRK格式的时间步长，ui和vi表示中间变量。
为保证该s阶显式SPRK格式具有s阶精度，可以通过泰勒级数展开的方式确定系数向量c和d。为了使得s阶显式SPRK格式达到s阶精度，至少需要s级格式才能满足要求。修正SPRK格式如下：
Figure BDA0002360884700000162
其中，引入了修正项Lv1，可以通过泰勒级数展开的方式确定系数向量c和d，以使得一个二级格式具有三阶时间精度。
通过计算，修正SPRK格式的系数向量可以为：
Figure BDA0002360884700000163
在步骤S208中，基于修正SPRK时间格式和高阶空间偏导数的离散格式，确定三维非均匀介质弹性波方程的数值解。
通过上述处理过程，可以得到三维非均匀介质弹性波方程正演模拟的修正时空优化保辛(MTSOS)方法。利用傅里叶分析的方法，可以得到时间离散算子的谱半径，进而得到MTSOS方法的稳定性条件，具体可以参见下述步骤S210的处理。
在步骤S210中，根据目标区域中声波或弹性波的初始波速、预定空间算子对应的特征值，确定修正SPRK时间格式的最大库朗数，以确定修正SPRK时间格式的稳定性条件。
在实施中，对于三维均匀介质中的声波方程，将谐波解
Figure BDA0002360884700000164
代入声波方程，其中，ωnum表示数值角频率，k＝(kx,ky,kz)表示波矢量，Δt表示时间步长，h＝Δx＝Δy＝Δz表示空间步长。可以得到
Figure BDA0002360884700000171
其中，增长矩阵G为
Figure BDA0002360884700000172
I表示单位算子，L表示空间算子。设矩阵G的特征值为ψ，算子L的特征值为χ。由于弹性波方程是双曲型偏微分方程，故其空间离散算子的特征值χ＜0。由公式(20)，ψ的特征多项式f(ψ)和χ满足关系：
Figure BDA0002360884700000173
为保证格式稳定，需要满足|ψ|≤1，即
Figure BDA0002360884700000174
对上述不等式进行求解，可以得到
-12≤Δt2χ≤0 (23)
设库朗数
Figure BDA0002360884700000175
则有：
Figure BDA0002360884700000176
上式对空间算子L的任意特征值χ均成立，因此，修正SPRK时间格式的稳定性条件为
Figure BDA0002360884700000177
其中，c0表示所述目标区域中声波的初始波速，αmax表示最大库朗数，χmin表示所述预定空间算子对应的最小特征值。
此外，上述步骤S210的处理可以多种多样，以下提供一种可选的处理方式，具体可以包括以下内容：根据所述目标区域中声波的初始波速、所述预定空间算子对应的特征值和库朗数，通过下述公式
Figure BDA0002360884700000178
确定所述修正SPRK时间格式的稳定性条件，其中，c0表示所述目标区域中声波的初始波速，αmax表示最大库朗数，χmin表示所述预定空间算子对应的最小特征值。
可以根据不同精度的空间算子L的最小特征值，得到与之对应的最大库朗数。在表3中分别给出了三维情形下不同精度的MTSOS及TSOS方法的最大库朗数。
表3
Figure BDA0002360884700000181
由表3可知，MTSOS方法的稳定性相比TSOS方法有了明显的提升。
对于非均匀情形，可采用冻结系数法的思路，将库朗数定义为
Figure BDA0002360884700000182
其中cmax表示目标区域中的最大波速。而对于弹性波情形，可以将库朗数定义为
Figure BDA0002360884700000183
其中vp表示P波的最大波速。
为了对修正SPRK时间格式进行数值频散分析，首先定义数值频散比率为如下形式：
Figure BDA0002360884700000184
其中，cnum为数值波速，λgrid为数值波长，ωnum为数值角频率，
Figure BDA0002360884700000185
为空间采样率，且由Nyquist-Shannon采样定理可知有，Sp≤0.5时格式才有意义。
将谐波解公式(18)代入公式(19)，可以得到
Figure BDA0002360884700000186
由于方程必有非零解，因此有：
det(exp(iωnumΔt)I-G)＝0 (29)
即
ωnumΔt＝arccos(Re(ψ)) (30)
又由公式(21)可知：
Figure BDA0002360884700000191
因此，最终可得到
Figure BDA0002360884700000192
上述公式(32)即为数值频散关系式，它表示了数值频散比率R与库朗数、采样率及空间算子的特征值之间的关系。R越接近于1，说明数值频散越小。
为了说明MTSOS方法的优势，以下将比较三维情形下MTSOS方法与SPRK方法的数值频散比率。
在图3和图4中分别给出了在α＝0.1,Sp＝0.4时不同空间精度的MTSOS方法和SPRK方法在不同方向上的频散误差|1-R|。从图中可以得到，MTSOS方法与SPRK方法的数值频散误差都表现出了各向异性的特点，且数值频散误差的最大值都出现在沿坐标轴的方向，而在与坐标轴呈45度角的方向，数值频散误差最小。同时，随着空间精度的提升，数值频散误差也呈现出下降的趋势。
从图中还可以得到，从总体上看，MTSOS方法的数值频散误差要小于SPRK方法。为了量化比较两种方法的数值频散误差，在表4中分别给出了不同空间精度的MTSOS方法与SPRK方法沿各个方向的最大数值频散误差。可以得到，不同空间精度的MTSOS方法的最大频散误差均小于相同空间精度的SPRK方法的最大频散误差，而且随着空间精度的提升，MTSOS方法对于最大频散误差的改进更加明显。平均而言，在α＝0.1,Sp＝0.4的情形中，MTSOS方法的最大数值频散误差为SPRK方法的57.5％。
表4
Figure BDA0002360884700000193
下面将MTSOS方法与边界修正法、吸收边界条件等结合，并选取两个模型，以考察MTSOS方法和吸收边界条件的有效性，以及边界修正法和吸收边界条件对计算区域边界的处理效果。
首先，选取三维均匀介质模型，利用MTSOS方法来进行波场模拟。目标区域为40km×40km×40km，介质密度为2660kg/m3，P波速度为5.8km/s，S波速度为3.199km/s。震源位于(20km,20km,18km)处。震源取为主频2Hz的点震源，震源函数为雷克子波
Figure BDA0002360884700000201
震源位于区域中心。时间步长为2.5ms，空间步长为0.2km。计算区域的上方为自由地表边界条件，其它五个边界采用吸收边界条件，吸收层取为15层。采用六阶空间精度的MTSOS方法进行波场模拟。
图5给出了t＝3.0s时刻位移场不同分量分别在震源所在XY、XZ和YZ平面的波场快照。从图可以清晰地发现直达P波及直达S波，且没有可见的数值频散。
为了验证MTSOS方法的准确性，在(23km,24km,18km)处设置一个接收器，并将接收器处的波形记录与解析解进行对比。图6是在0-5s内位移三分量波形与解析解的比较，其中实线表示解析解，虚线表示数值解。从图中可以看出，接收器处的波形记录与解析解能够很好地吻合，说明数值模拟结果是可信的。
下面选取一个三维双层介质模型，来检验MTSOS方法在利用等效介质参数法处理模型间断以及吸收边界条件的有效性。目标区域为20km×20km×20km，上层介质密度为2600kg/m3，P波速度为5.8km/s，S波速度为3.2km/s；下层介质密度为3723kg/m3，P波速度为9.134km/s，S波速度为4.932km/s，间断面位于z＝20km处。震源取为主频8Hz的点震源，震源函数为雷克子波。震源位于(10km,10km,8km)处。时间步长为1s，空间步长为0.08km。计算区域的上方为自由地表边界条件，其他五个边界采用吸收边界条件，吸收层取为15层。采用六阶空间精度的MTSOS方法进行波场模拟。
图7给出了t＝1.75s时刻的位移场的不同分量分别在震源所在的XY、XZ和YZ平面的波场快照。从图可以清晰地发现直达波、反射波及透射波的各个震相以及产生的转换波。在XY平面，由于不存在间断面，波场呈同心圆结构，不同的圆表示了以不同速度传播的波。其中，最外层为直达P波(P)，向内依次为反射P波(PP)、S波反射后的转换P波(SP)、直达S波(S)、P波反射后的转换S波(PS)、反射S波(SS)等。在XZ及YZ平面，可以看到经间断面的反射、透射和转换波，而且没有可见的数值频散。图8选取了t＝1.75s时位移场x分量在XZ平面的波场快照，并标注出所有产生的震相。
在图9中以位移场的x分量在XZ平面为例给出了从t＝0.5s至t＝2.5s时刻的波场快照。从图中可以看出波场随时间稳定的传播：最初时刻仅有直达P波及直达S波。当波传播到间断面，产生了各种反射、透射及转换波。各种震相均稳定且无可见的数值频散。约t＝1.5s时，透射P波(Pp)传播到边界；约t＝2.0s时，其它波的震相也陆续传播到边界。其中，自由地表处的直达P波经过自由地表的反射形成反射P波和转换S波，而在其他的边界处，波场均被吸收边界条件很好地吸收，未产生可见的反射波。该数值实验很好地说明等效介质参数可以准确地处理模型的间断，边界修正法与吸收边界条件可以很好地与数值算法结合，给出准确的边界条件。
下面选取一个三维含水裂隙介质模型来考察MTSOS方法对于细微的非均匀结构的分辨能力。含水裂隙介质是在石油勘探领域经常会接触到的地质结构。可以向井中注水，用以驱油并减少开采时间。因此，模拟弹性波在含水裂隙介质模型中的传播对于石油勘探具有重要的意义。
在这个例子中，将目标区域取为2km×2km×2km，背景介质密度为2170kg/m3，P波速度为2.618km/s，S波速度为1.421km/s。震源取为主频40Hz的点震源，震源函数为雷克子波。震源位于(1km,1km,1km)处。时间步长为0.2ms，空间步长为5m。目标区域的上方为自由地表边界条件，其它五个平面采用吸收边界条件，吸收层取为15层。在(0.9km,0.9km,1km)至(1km,1km,1.1km)有一个宽度为5m的含水裂隙，裂隙的两端恰好分别位于震源所在的XY和YZ、XZ平面内。在含水裂隙内部，介质密度为1000kg/m3，P波速度为1.5km/s，S波速度为0km/s。采用六阶空间精度的MTSOS方法进行波场模拟。在t＝0.25s时刻位移场不同分量于震源所在的XY、XZ和YZ平面的波场快照如图10所示。在图10中，可以清晰地发现由于含水裂隙产生的散射波。为了更准确地显示散射波的形成过程，分别在(0.925km,0.92km,1km)和(0.885km,0.88km,1km)各放置一个接收器。注意到第一个接收器和震源位于含水裂隙的同侧，而第二个接收器和震源分别位于含水裂隙的两侧。图11为六阶MTSOS方法计算得到的位于接收器处的波形图，其中，(a)为x分量；(b)为x分量局部放大图；(c)为y分量；(d)为y分量局部放大图；(e)为z分量；(f)为z分量局部放大图，实线为不包含含水裂隙时的波形图，虚线为包含含水裂隙时的波形图。图12为六阶MTSOS方法计算得到的位于(0.885km,0.88km,1km)的接收器的波形图。(a)为x分量；(b)为x分量局部放大图；(c)为y分量；(d)为y分量局部放大图；(e)为z分量；(f)为z分量局部放大图，其中实线为不包含含水裂隙时的波形图，虚线为包含含水裂隙时的波形图。右侧为波形图的局部放大图。可以看出，对于位于震源同侧的接收器，由于含水裂隙的存在，在直达S波后可清晰地发现散射波所形成的一系列尾波。而对于位于震源异侧的接收器，散射P波与直达S波几乎同时到达，使得直达S波的振幅与没有含水裂隙时的情形有了明显的区别。该数值实验说明MTSOS方法可以准确模拟勘探尺度下模型内的细微结构变化对于波形的影响。
本说明书实施例应用于三维非均匀介质弹性波方程正演模拟的新方法。该方法在时间上应用修正的SPRK格式，利用一个二级格式达到了三阶时间精度；在空间上应用优化有限差分格式，并将其拓展到三维，同时将等效介质参数、吸收边界条件等结合，进而获得了新的三维MTSOS方法。由于优化有限差分格式及等效介质参数的引入，三维MTSOS方法更适用于非均匀介质模型的正演模拟。同时，本说明书实施例还给出了MTSOS方法的稳定性条件及数值频散关系。理论及数值结果均表明，MTSOS方法的数值频散误差要低于相同空间精度的SPRK方法，而且随着空间精度的提升，这种优势也越来越明显。总体上，新方法的数值频散误差约为SPRK方法的57.5％。波场模拟实验进一步验证了利用MTSOS方法可以精确给出波场模拟结果，能够清晰显示波传播过程中产生的各种震相，且没有可见的数值频散，可以模拟非均匀结构对于波场的影响，说明了MTSOS方法的有效性和准确性。
本说明书实施例提供一种三维非均匀介质弹性波的处理方法，通过根据待勘探的目标区域、当前具备的计算资源和预定的库朗数，确定目标区域对应的时间步长、空间步长和空间精度，根据空间精度确定有限差分格式优化系数，并根据目标区域对应的空间步长和有限差分格式优化系数，确定高阶空间偏导数的离散格式，基于目标区域对应的时间步长、预定系数向量、预定空间算子，确定修正保辛分部龙格库塔SPRK时间格式，基于修正SPRK时间格式和高阶空间偏导数的离散格式，确定三维非均匀介质弹性波方程的数值解，这样，通过优化函数，最终得到最优的有限差分格式优化系数，从而可以大大节省内存需求和计算成本，而且，修正SPRK格式通过引入修正项，可以使一个二级格式达到了三阶时间精度。
实施例三
以上为本说明书实施例提供的三维非均匀介质弹性波的处理方法，基于同样的思路，本说明书实施例还提供一种三维非均匀介质弹性波的处理装置，如图13所示。
所述三维非均匀介质弹性波的处理装置包括：区域信息确定模块1301、离散格式确定模块1302、修正格式确定模块1303和求解模块1304，其中：
区域信息确定模块1301，用于根据待勘探的目标区域、当前具备的计算资源和最大库朗数，确定所述目标区域对应的时间步长、空间步长和空间精度；
离散格式确定模块1302，用于根据所述空间精度确定有限差分格式优化系数，并根据所述目标区域对应的空间步长和所述有限差分格式优化系数，确定高阶空间偏导数的离散格式；
修正格式确定模块1303，用于基于所述目标区域对应的时间步长、预定系数向量、预定空间算子，确定修正保辛分部龙格库塔SPRK时间格式；
求解模块1304，用于基于所述修正SPRK时间格式和所述高阶空间偏导数的离散格式，确定三维非均匀介质弹性波方程的数值解。
本说明书实施例中，所述离散格式确定模块1302，用于根据所述目标区域对应的空间步长和所述有限差分格式优化系数，分别通过公式
Figure BDA0002360884700000231
Figure BDA0002360884700000232
Figure BDA0002360884700000233
Figure BDA0002360884700000234
Figure BDA0002360884700000235
Figure BDA0002360884700000236
Figure BDA0002360884700000237
Figure BDA0002360884700000238
Figure BDA0002360884700000241
Figure BDA0002360884700000242
Figure BDA0002360884700000243
Figure BDA0002360884700000244
确定高阶空间偏导数的离散格式，其中，c、λ、μ表示预定系数，Δx、Δy和Δz分别表示所述目标区域中x方向的空间步长、y方向的空间步长和z方向的空间步长，
Figure BDA0002360884700000245
和
Figure BDA0002360884700000246
为所述有限差分格式化系数。
本说明书实施例中，所述修正格式确定模块1303，用于基于所述目标区域对应的时间步长、预定系数向量、预定空间算子，确定以下修正保辛分部龙格库塔SPRK时间格式：
Figure BDA0002360884700000247
其中，Δt表示所述目标区域对应的时间步长，c＝(c1,c2,c3)和d＝(d1,d2)表示所述预定系数向量，L表示所述预定空间算子，un和vn表示在第n个时间步长的数值解，un+1和vn +1表示在第n+1个时间步长的数值解，u1和v1表示中间变量。
本说明书实施例中，所述装置还包括：
稳定性条件确定模块，用于根据所述目标区域中声波或弹性波的初始波速、所述预定空间算子对应的特征值，确定所述修正SPRK时间格式的最大库朗数，以确定所述修正SPRK时间格式的稳定性条件。
本说明书实施例中，所述稳定性条件确定模块，用于根据所述目标区域中声波或弹性波的初始波速、所述预定空间算子对应的特征值，通过下述公式
Figure BDA0002360884700000248
确定所述修正SPRK时间格式的最大库朗数，以确定所述修正SPRK时间格式的稳定性条件，其中，c0表示所述目标区域中P波的的最大波速，αmax表示最大库郎数，χmin表示所述预定空间算子对应的最小特征值。
本说明书实施例提供一种三维非均匀介质弹性波的处理装置，通过根据待勘探的目标区域、当前具备的计算资源和最大库朗数，确定目标区域对应的时间步长、空间步长和空间精度，根据空间精度确定有限差分格式优化系数，并根据目标区域对应的空间步长和有限差分格式优化系数，确定高阶空间偏导数的离散格式，基于目标区域对应的时间步长、预定系数向量、预定空间算子，确定修正保辛分部龙格库塔SPRK时间格式，基于修正SPRK时间格式和高阶空间偏导数的离散格式，确定三维非均匀介质弹性波方程的数值解，这样，通过优化函数，最终得到最优的有限差分格式优化系数，从而可以大大节省内存需求和计算成本，而且，修正SPRK格式通过引入修正项，可以使一个二级格式达到了三阶时间精度。
实施例四
基于同样的思路，本说明书实施例还提供一种三维非均匀介质弹性波的处理设备，如图14所示。
该三维非均匀介质弹性波的处理设备可以为上述实施例提供的终端设备或服务器等。
三维非均匀介质弹性波的处理设备可因配置或性能不同而产生比较大的差异，可以包括一个或一个以上的处理器1401和存储器1402，存储器1402中可以存储有一个或一个以上存储应用程序或数据。其中，存储器1402可以是短暂存储或持久存储。存储在存储器1402的应用程序可以包括一个或一个以上模块(图示未示出)，每个模块可以包括对三维非均匀介质弹性波的处理设备中的一系列计算机可执行指令。更进一步地，处理器1401可以设置为与存储器1402通信，在三维非均匀介质弹性波的处理设备上执行存储器1402中的一系列计算机可执行指令。三维非均匀介质弹性波的处理设备还可以包括一个或一个以上电源1403，一个或一个以上有线或无线网络接口1404，一个或一个以上输入输出接口1405，一个或一个以上键盘1406。
具体在本实施例中，三维非均匀介质弹性波的处理设备包括有存储器，以及一个或一个以上的程序，其中一个或者一个以上程序存储于存储器中，且一个或者一个以上程序可以包括一个或一个以上模块，且每个模块可以包括对三维非均匀介质弹性波的处理设备中的一系列计算机可执行指令，且经配置以由一个或者一个以上处理器执行该一个或者一个以上程序包含用于进行以下计算机可执行指令：
根据待勘探的目标区域、当前具备的计算资源和最大库朗数，确定所述目标区域对应的时间步长、空间步长和空间精度；
根据所述空间精度确定有限差分格式优化系数，并根据所述目标区域对应的空间步长和所述有限差分格式优化系数，确定高阶空间偏导数的离散格式；
基于所述目标区域对应的时间步长、预定系数向量、预定空间算子，确定修正保辛分部龙格库塔SPRK时间格式；
基于所述修正SPRK时间格式和所述高阶空间偏导数的离散格式，确定三维非均匀介质弹性波方程的数值解。
本说明书实施例中，所述根据所述目标区域对应的空间步长和所述有限差分格式优化系数，确定高阶空间偏导数的离散格式，包括：
根据所述目标区域对应的空间步长和所述有限差分格式优化系数，分别通过公式
Figure BDA0002360884700000261
Figure BDA0002360884700000262
Figure BDA0002360884700000263
Figure BDA0002360884700000264
Figure BDA0002360884700000265
Figure BDA0002360884700000266
Figure BDA0002360884700000267
Figure BDA0002360884700000268
Figure BDA0002360884700000271
Figure BDA0002360884700000272
Figure BDA0002360884700000273
Figure BDA0002360884700000274
确定高阶空间偏导数的离散格式，其中，c、λ、μ表示预定系数，Δx、Δy和Δz分别表示所述目标区域中x方向的空间步长、y方向的空间步长和z方向的空间步长，
Figure BDA0002360884700000275
和
Figure BDA0002360884700000276
为所述有限差分格式化系数。
本说明书实施例中，所述基于所述目标区域对应的时间步长、预定系数向量、预定空间算子，确定修正保辛分部龙格库塔SPRK时间格式，包括：
基于所述目标区域对应的时间步长、预定系数向量、预定空间算子，确定以下修正保辛分部龙格库塔SPRK时间格式：
Figure BDA0002360884700000277
其中，Δt表示所述目标区域对应的时间步长，c＝(c1,c2,c3)和d＝(d1,d2)表示所述预定系数向量，L表示所述预定空间算子，un和vn表示在第n个时间步长的数值解，un+1和vn +1表示在第n+1个时间步长的数值解，u1和v1表示中间变量。
本说明书实施例中，还包括：
根据所述目标区域中声波或弹性波的初始波速、所述预定空间算子对应的特征值和，确定所述修正SPRK时间格式的最大库朗数，以确定所述修正SPRK时间格式的稳定性条件。
本说明书实施例中，所述根据所述目标区域中声波或弹性波的初始波速、所述预定空间算子对应的特征值，确定所述修正SPRK时间格式的最大库朗数，以确定所述修正SPRK时间格式的稳定性条件，包括：
根据所述目标区域中声波或弹性波的初始波速、所述预定空间算子对应的特征值，通过下述公式
Figure BDA0002360884700000278
确定所述修正SPRK时间格式的最大库朗数，以确定所述修正SPRK时间格式的稳定性条件，其中，c0表示所述目标区域中P波的最大波速，αmax表示最大库朗数，χmin表示所述预定空间算子对应的最小特征值。
本说明书实施例提供一种三维非均匀介质弹性波的处理设备，通过根据待勘探的目标区域、当前具备的计算资源和最大库朗数，确定目标区域对应的时间步长、空间步长和空间精度，根据空间精度确定有限差分格式优化系数，并根据目标区域对应的空间步长和有限差分格式优化系数，确定高阶空间偏导数的离散格式，基于目标区域对应的时间步长、预定系数向量、预定空间算子，确定修正保辛分部龙格库塔SPRK时间格式，基于修正SPRK时间格式和高阶空间偏导数的离散格式，确定三维非均匀介质弹性波方程的数值解，这样，通过优化函数，最终得到最优的有限差分格式优化系数，从而可以大大节省内存需求和计算成本，而且，修正SPRK格式通过引入修正项，可以使一个二级格式达到了三阶时间精度。
上述对本说明书特定实施例进行了描述。其它实施例在所附权利要求书的范围内。在一些情况下，在权利要求书中记载的动作或步骤可以按照不同于实施例中的顺序来执行并且仍然可以实现期望的结果。另外，在附图中描绘的过程不一定要求示出的特定顺序或者连续顺序才能实现期望的结果。在某些实施方式中，多任务处理和并行处理也是可以的或者可能是有利的。
在20世纪90年代，对于一个技术的改进可以很明显地区分是硬件上的改进(例如，对二极管、晶体管、开关等电路结构的改进)还是软件上的改进(对于方法流程的改进)。然而，随着技术的发展，当今的很多方法流程的改进已经可以视为硬件电路结构的直接改进。设计人员几乎都通过将改进的方法流程编程到硬件电路中来得到相应的硬件电路结构。因此，不能说一个方法流程的改进就不能用硬件实体模块来实现。例如，可编程逻辑器件(Programmable Logic Device，PLD)(例如现场可编程门阵列(Field Programmable GateArray，FPGA))就是这样一种集成电路，其逻辑功能由用户对器件编程来确定。由设计人员自行编程来把一个数字系统“集成”在一片PLD上，而不需要请芯片制造厂商来设计和制作专用的集成电路芯片。而且，如今，取代手工地制作集成电路芯片，这种编程也多半改用“逻辑编译器(logic compiler)”软件来实现，它与程序开发撰写时所用的软件编译器相类似，而要编译之前的原始代码也得用特定的编程语言来撰写，此称之为硬件描述语言(Hardware Description Language，HDL)，而HDL也并非仅有一种，而是有许多种，如ABEL(Advanced Boolean Expression Language)、AHDL(Altera Hardware DescriptionLanguage)、Confluence、CUPL(Cornell University Programming Language)、HDCal、JHDL(Java Hardware Description Language)、Lava、Lola、MyHDL、PALASM、RHDL(RubyHardware Description Language)等，目前最普遍使用的是VHDL(Very-High-SpeedIntegrated Circuit Hardware Description Language)与Verilog。本领域技术人员也应该清楚，只需要将方法流程用上述几种硬件描述语言稍作逻辑编程并编程到集成电路中，就可以很容易得到实现该逻辑方法流程的硬件电路。
控制器可以按任何适当的方式实现，例如，控制器可以采取例如微处理器或处理器以及存储可由该(微)处理器执行的计算机可读程序代码(例如软件或固件)的计算机可读介质、逻辑门、开关、专用集成电路(Application Specific Integrated Circuit，ASIC)、可编程逻辑控制器和嵌入微控制器的形式，控制器的例子包括但不限于以下微控制器：ARC 625D、Atmel AT91SAM、Microchip PIC18F26K20以及Silicone Labs C8051F320，存储器控制器还可以被实现为存储器的控制逻辑的一部分。本领域技术人员也知道，除了以纯计算机可读程序代码方式实现控制器以外，完全可以通过将方法步骤进行逻辑编程来使得控制器以逻辑门、开关、专用集成电路、可编程逻辑控制器和嵌入微控制器等的形式来实现相同功能。因此这种控制器可以被认为是一种硬件部件，而对其内包括的用于实现各种功能的装置也可以视为硬件部件内的结构。或者甚至，可以将用于实现各种功能的装置视为既可以是实现方法的软件模块又可以是硬件部件内的结构。
上述实施例阐明的系统、装置、模块或单元，具体可以由计算机芯片或实体实现，或者由具有某种功能的产品来实现。一种典型的实现设备为计算机。具体的，计算机例如可以为个人计算机、膝上型计算机、蜂窝电话、相机电话、智能电话、个人数字助理、媒体播放器、导航设备、电子邮件设备、游戏控制台、平板计算机、可穿戴设备或者这些设备中的任何设备的组合。
为了描述的方便，描述以上装置时以功能分为各种单元分别描述。当然，在实施本说明书一个或多个实施例时可以把各单元的功能在同一个或多个软件和/或硬件中实现。
本领域内的技术人员应明白，本说明书的实施例可提供为方法、系统、或计算机程序产品。因此，本说明书一个或多个实施例可采用完全硬件实施例、完全软件实施例、或结合软件和硬件方面的实施例的形式。而且，本说明书一个或多个实施例可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、CD-ROM、光学存储器等)上实施的计算机程序产品的形式。
本说明书的实施例是参照根据本说明书实施例的方法、设备(系统)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程基于区块链的用户目标处理设备的处理器以产生一个机器，使得通过计算机或其他可编程基于区块链的用户目标处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。
这些计算机程序指令也可存储在能引导计算机或其他可编程基于区块链的用户目标处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。
这些计算机程序指令也可装载到计算机或其他可编程基于区块链的用户目标处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。
在一个典型的配置中，计算设备包括一个或多个处理器(CPU)、输入/输出接口、网络接口和内存。
内存可能包括计算机可读介质中的非永久性存储器，随机存取存储器(RAM)和/或非易失性内存等形式，如只读存储器(ROM)或闪存(flash RAM)。内存是计算机可读介质的示例。
计算机可读介质包括永久性和非永久性、可移动和非可移动媒体可以由任何方法或技术来实现信息存储。信息可以是计算机可读指令、数据结构、程序的模块或其他数据。计算机的存储介质的例子包括，但不限于相变内存(PRAM)、静态随机存取存储器(SRAM)、动态随机存取存储器(DRAM)、其他类型的随机存取存储器(RAM)、只读存储器(ROM)、电可擦除可编程只读存储器(EEPROM)、快闪记忆体或其他内存技术、只读光盘只读存储器(CD-ROM)、数字多功能光盘(DVD)或其他光学存储、磁盒式磁带，磁带磁磁盘存储或其他磁性存储设备或任何其他非传输介质，可用于存储可以被计算设备访问的信息。按照本文中的界定，计算机可读介质不包括暂存电脑可读媒体(transitory media)，如调制的数据信号和载波。
还需要说明的是，术语“包括”、“包含”或者其任何其他变体意在涵盖非排他性的包含，从而使得包括一系列要素的过程、方法、商品或者设备不仅包括那些要素，而且还包括没有明确列出的其他要素，或者是还包括为这种过程、方法、商品或者设备所固有的要素。在没有更多限制的情况下，由语句“包括一个……”限定的要素，并不排除在包括所述要素的过程、方法、商品或者设备中还存在另外的相同要素。
本领域技术人员应明白，本说明书的实施例可提供为方法、系统或计算机程序产品。因此，本说明书一个或多个实施例可采用完全硬件实施例、完全软件实施例或结合软件和硬件方面的实施例的形式。而且，本说明书一个或多个实施例可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、CD-ROM、光学存储器等)上实施的计算机程序产品的形式。
本说明书一个或多个实施例可以在由计算机执行的计算机可执行指令的一般上下文中描述，例如程序模块。一般地，程序模块包括执行特定任务或实现特定抽象数据类型的例程、程序、对象、组件、数据结构等等。也可以在分布式计算环境中实践本说明书一个或多个实施例，在这些分布式计算环境中，由通过通信网络而被连接的远程处理设备来执行任务。在分布式计算环境中，程序模块可以位于包括存储设备在内的本地和远程计算机存储介质中。
本说明书中的各个实施例均采用递进的方式描述，各个实施例之间相同相似的部分互相参见即可，每个实施例重点说明的都是与其他实施例的不同之处。尤其，对于系统实施例而言，由于其基本相似于方法实施例，所以描述的比较简单，相关之处参见方法实施例的部分说明即可。
以上所述仅为本说明书的实施例而已，并不用于限制本说明书。对于本领域技术人员来说，本说明书可以有各种更改和变化。凡在本说明书的精神和原理之内所作的任何修改、等同替换、改进等，均应包含在本说明书的权利要求范围之内。

Claims (10)

1.一种三维非均匀介质弹性波的处理方法，其特征在于，所述方法包括：
根据待勘探的目标区域、当前具备的计算资源和最大库朗数，确定所述目标区域对应的时间步长、空间步长和空间精度；
根据所述空间精度确定有限差分格式优化系数，并根据所述目标区域对应的空间步长和所述有限差分格式优化系数，确定高阶空间偏导数的离散格式；
基于所述目标区域对应的时间步长、预定系数向量、预定空间算子，确定修正保辛分部龙格库塔SPRK时间格式；
基于所述修正SPRK时间格式和所述高阶空间偏导数的离散格式，确定三维非均匀介质弹性波方程的数值解。
2.根据权利要求1所述的方法，其特征在于，所述根据所述目标区域对应的空间步长和所述有限差分格式优化系数，确定高阶空间偏导数的离散格式，包括：
根据所述目标区域对应的空间步长和所述有限差分格式优化系数，分别通过公式
Figure FDA0002360884690000011
Figure FDA0002360884690000012
Figure FDA0002360884690000013
Figure FDA0002360884690000014
Figure FDA0002360884690000015
Figure FDA0002360884690000021
Figure FDA0002360884690000022
Figure FDA0002360884690000023
Figure FDA0002360884690000024
Figure FDA0002360884690000025
Figure FDA0002360884690000026
Figure FDA0002360884690000027
确定高阶空间偏导数的离散格式，其中，c、λ、μ表示预定系数，Δx、Δy和Δz分别表示所述目标区域中x方向的空间步长、y方向的空间步长和z方向的空间步长，
Figure FDA0002360884690000028
和
Figure FDA0002360884690000029
为所述有限差分格式化系数。
3.根据权利要求2所述的方法，其特征在于，所述基于所述目标区域对应的时间步长、预定系数向量、预定空间算子，确定修正保辛分部龙格库塔SPRK时间格式，包括：
基于所述目标区域对应的时间步长、预定系数向量、预定空间算子，确定以下修正保辛分部龙格库塔SPRK时间格式：
Figure FDA00023608846900000210
其中，Δt表示所述目标区域对应的时间步长，c＝(c1,c2,c3)和d＝(d1,d2)表示所述预定系数向量，L表示所述预定空间算子，un和vn表示在第n个时间步长的数值解，un+1和vn+1表示在第n+1个时间步长的数值解，u1和v1表示中间变量。
4.根据权利要求3所述的方法，其特征在于，所述方法还包括：
根据所述目标区域中声波或弹性波的初始波速、所述预定空间算子对应的特征值，确定所述修正SPRK时间格式的最大库朗数，以确定所述修正SPRK时间格式的稳定性条件。
5.根据权利要求4所述的方法，其特征在于，所述根据所述目标区域中声波或弹性波的初始波速、所述预定空间算子对应的特征值，确定所述修正SPRK时间格式的最大库朗数，以确定所述修正SPRK时间格式的稳定性条件，包括：
根据所述目标区域中声波或弹性波的初始波速、所述预定空间算子对应的特征值，通过下述公式
Figure FDA0002360884690000031
确定所述修正SPRK时间格式的最大库朗数，以确定所述修正SPRK时间格式的稳定性条件，其中，c0表示所述目标区域中P波的最大波速，αmax表示最大库朗数，χmin表示所述预定空间算子对应的最小特征值。
6.一种三维非均匀介质弹性波的处理装置，其特征在于，所述装置包括：
区域信息确定模块，用于根据待勘探的目标区域、当前具备的计算资源和最大库朗数，确定所述目标区域对应的时间步长、空间步长和空间精度；
离散格式确定模块，用于根据所述空间精度确定有限差分格式优化系数，并根据所述目标区域对应的空间步长和所述有限差分格式优化系数，确定高阶空间偏导数的离散格式；
修正格式确定模块，用于基于所述目标区域对应的时间步长、预定系数向量、预定空间算子，确定修正保辛分部龙格库塔SPRK时间格式；
求解模块，用于基于所述修正SPRK时间格式和所述高阶空间偏导数的离散格式，确定三维非均匀介质弹性波方程的数值解。
7.根据权利要求6所述的装置，其特征在于，所述离散格式确定模块，用于根据所述目标区域对应的空间步长和所述有限差分格式优化系数，分别通过公式
Figure FDA0002360884690000032
Figure FDA0002360884690000041
Figure FDA0002360884690000042
Figure FDA0002360884690000043
Figure FDA0002360884690000044
Figure FDA0002360884690000045
Figure FDA0002360884690000046
Figure FDA0002360884690000047
Figure FDA0002360884690000048
Figure FDA0002360884690000049
Figure FDA00023608846900000410
Figure FDA00023608846900000411
确定高阶空间偏导数的离散格式，其中，c、λ、μ表示预定系数，Δx、Δy和Δz分别表示所述目标区域中x方向的空间步长、y方向的空间步长和z方向的空间步长，
Figure FDA00023608846900000412
和
Figure FDA00023608846900000413
为所述有限差分格式化系数。
8.根据权利要求7所述的装置，其特征在于，所述修正格式确定模块，用于基于所述目标区域对应的时间步长、预定系数向量、预定空间算子，确定以下修正保辛分部龙格库塔SPRK时间格式：
Figure FDA0002360884690000051
其中，Δt表示所述目标区域对应的时间步长，c＝(c1,c2,c3)和d＝(d1,d2)表示所述预定系数向量，L表示所述预定空间算子，un和vn表示在第n个时间步长的数值解，un+1和vn+1表示在第n+1个时间步长的数值解，u1和v1表示中间变量。
9.根据权利要求8所述的装置，其特征在于，所述装置还包括：
稳定性条件确定模块，用于根据所述目标区域中声波或弹性波的初始波速、所述预定空间算子对应的特征值，确定所述修正SPRK时间格式的最大库朗数，以确定所述修正SPRK时间格式的稳定性条件。
10.一种三维非均匀介质弹性波的处理设备，其特征在于，所述三维非均匀介质弹性波的处理设备包括：
处理器；以及
被安排成存储计算机可执行指令的存储器，所述可执行指令在被执行时使所述处理器：
根据待勘探的目标区域、当前具备的计算资源和最大库朗数，确定所述目标区域对应的时间步长、空间步长和空间精度；
根据所述空间精度确定有限差分格式优化系数，并根据所述目标区域对应的空间步长和所述有限差分格式优化系数，确定高阶空间偏导数的离散格式；
基于所述目标区域对应的时间步长、预定系数向量、预定空间算子，确定修正保辛分部龙格库塔SPRK时间格式；
基于所述修正SPRK时间格式和所述高阶空间偏导数的离散格式，确定三维非均匀介质弹性波方程的数值解。
CN202010021263.6A 2020-01-09 2020-01-09 一种三维非均匀介质弹性波的处理方法、装置和设备 Active CN111142157B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010021263.6A CN111142157B (zh) 2020-01-09 2020-01-09 一种三维非均匀介质弹性波的处理方法、装置和设备

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010021263.6A CN111142157B (zh) 2020-01-09 2020-01-09 一种三维非均匀介质弹性波的处理方法、装置和设备

Publications (2)

Publication Number Publication Date
CN111142157A true CN111142157A (zh) 2020-05-12
CN111142157B CN111142157B (zh) 2021-01-26

Family

ID=70524189

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010021263.6A Active CN111142157B (zh) 2020-01-09 2020-01-09 一种三维非均匀介质弹性波的处理方法、装置和设备

Country Status (1)

Country Link
CN (1) CN111142157B (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116893447A (zh) * 2023-09-11 2023-10-17 中国科学院武汉岩土力学研究所 一种基于瞬变时域地质雷达频散差分高效分析方法

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102053258A (zh) * 2010-12-15 2011-05-11 中国石油集团川庆钻探工程有限公司 基于复杂地质构造的自适应三维射线追踪方法
CN102162859A (zh) * 2011-01-10 2011-08-24 中国海洋石油总公司 一种斜井井间地震波场的成像方法
WO2014172011A2 (en) * 2013-04-16 2014-10-23 Exxonmobil Upstream Research Company Seismic velocity model updating and imaging with elastic wave imaging
CN106772589A (zh) * 2017-02-27 2017-05-31 中国石油大学(北京) 一种叠前地震反演方法及装置
CN110703322A (zh) * 2019-10-10 2020-01-17 清华大学 一种波传播的处理方法、装置和设备

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102053258A (zh) * 2010-12-15 2011-05-11 中国石油集团川庆钻探工程有限公司 基于复杂地质构造的自适应三维射线追踪方法
CN102162859A (zh) * 2011-01-10 2011-08-24 中国海洋石油总公司 一种斜井井间地震波场的成像方法
WO2014172011A2 (en) * 2013-04-16 2014-10-23 Exxonmobil Upstream Research Company Seismic velocity model updating and imaging with elastic wave imaging
CN106772589A (zh) * 2017-02-27 2017-05-31 中国石油大学(北京) 一种叠前地震反演方法及装置
CN110703322A (zh) * 2019-10-10 2020-01-17 清华大学 一种波传播的处理方法、装置和设备

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
马啸: "波动方程保辛近似解析离散化算法研究", 《中国博士学位论文全文数据库 基础科学辑》 *

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116893447A (zh) * 2023-09-11 2023-10-17 中国科学院武汉岩土力学研究所 一种基于瞬变时域地质雷达频散差分高效分析方法

Also Published As

Publication number Publication date
CN111142157B (zh) 2021-01-26

Similar Documents

Publication Publication Date Title
Florinski et al. Galactic cosmic ray transport in the global heliosphere
CN108983285B (zh) 一种基于矩张量的多种震源波场模拟方法及装置
Padovani et al. Low and high order finite element method: experience in seismic modeling
Marburg Boundary element method for time-harmonic acoustic problems
CN108387926B (zh) 一种确定气枪阵列远场子波的方法及装置
Ely et al. A support-operator method for viscoelastic wave modelling in 3-D heterogeneous media
van Driel et al. Accelerating numerical wave propagation using wavefield adapted meshes. Part I: forward and adjoint modelling
Zhang et al. Accelerating 3D Fourier migration with graphics processing units
Lambrecht et al. A nodal discontinuous Galerkin approach to 3-D viscoelastic wave propagation in complex geological media
Le Bouteiller et al. A discontinuous Galerkin fast-sweeping eikonal solver for fast and accurate traveltime computation in 3D tilted anisotropic media
Aseev et al. Numerical applications of the advective‐diffusive codes for the inner magnetosphere
Yang et al. An optimal nearly analytic discrete-weighted Runge-Kutta discontinuous Galerkin hybrid method for acoustic wavefield modeling
Baldassari et al. Numerical performances of a hybrid local‐time stepping strategy applied to the reverse time migration
CN111142157B (zh) 一种三维非均匀介质弹性波的处理方法、装置和设备
Akbudak et al. Asynchronous computations for solving the acoustic wave propagation equation
CN108828659B (zh) 基于傅里叶有限差分低秩分解的地震波场延拓方法及装置
van Driel et al. On the modelling of self-gravitation for full 3-D global seismic wave propagation
CN112505775B (zh) 地震波正演模拟方法、装置、存储介质及处理器
Zhang et al. A new spectral finite volume method for elastic wave modelling on unstructured meshes
CN109725346B (zh) 一种时间-空间高阶vti矩形网格有限差分方法及装置
CN111007565B (zh) 三维频率域全声波成像方法及装置
He et al. A numerical dispersion-dissipation analysis of discontinuous Galerkin methods based on quadrilateral and triangular elements
Gorodnitskiy et al. Depth migration with Gaussian wave packets based on Poincaré wavelets
Chen The LMARS based shallow‐water dynamical core on generic gnomonic cubed‐sphere geometry
Li et al. A 3D reflection ray-tracing method based on linear traveltime perturbation interpolation

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant