CN114460640A - 有限差分模拟弹性波全波形反演方法和装置 - Google Patents

有限差分模拟弹性波全波形反演方法和装置 Download PDF

Info

Publication number
CN114460640A
CN114460640A CN202011239881.4A CN202011239881A CN114460640A CN 114460640 A CN114460640 A CN 114460640A CN 202011239881 A CN202011239881 A CN 202011239881A CN 114460640 A CN114460640 A CN 114460640A
Authority
CN
China
Prior art keywords
longitudinal
wave
difference
velocity
elastic wave
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
CN202011239881.4A
Other languages
English (en)
Other versions
CN114460640B (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.)
China National Petroleum Corp
BGP Inc
Original Assignee
China National Petroleum Corp
BGP Inc
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 China National Petroleum Corp, BGP Inc filed Critical China National Petroleum Corp
Priority to CN202011239881.4A priority Critical patent/CN114460640B/zh
Publication of CN114460640A publication Critical patent/CN114460640A/zh
Application granted granted Critical
Publication of CN114460640B publication Critical patent/CN114460640B/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
    • G01V1/30—Analysis
    • G01V1/306—Analysis for determining physical properties of the subsurface, e.g. impedance, porosity or attenuation profiles
    • 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/62—Physical property of subsurface
    • 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/62—Physical property of subsurface
    • G01V2210/624—Reservoir parameters
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
    • Y02T90/00—Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation

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

本发明提供一种有限差分模拟弹性波全波形反演方法和装置，该方法包括：根据纵横波速度模型，确定速度‑拟应力弹性波动方程的纵横波差分参数；所述速度‑拟应力弹性波动方程的纵横波差分参数，包括：纵横波有限差分离散格式，纵横波差分系数，纵横波稳定性条件，纵横波频散条件，初始条件，边界条件；采用GPU共享内存优化策略，结合速度‑拟应力弹性波动方程的纵横波差分参数，建立速度‑拟应力弹性波动方程，进行时间四阶精度有限差分模拟，获得时间四阶精度有限差分模拟结果；根据时间四阶精度有限差分模拟结果和速度‑拟应力弹性波动方程，进行时间域弹性波全波形反演，确定纵横波速度模型。实现了高效的弹性波场模拟和纵横波速度模型反演。

Description

有限差分模拟弹性波全波形反演方法和装置
技术领域
本发明涉及地震数据处理技术领域，尤其涉及一种有限差分模拟弹性波全波形反演方法和装置。
背景技术
本部分旨在为权利要求书中陈述的本发明的实施方式提供背景或上下文。此处的描述不因为包括在本部分中就承认是现有技术。
随着勘探开发环境的复杂化，地震数据处理技术正面临着新的严峻考验。地震资料处理的关键任务就是对地下介质构造进行清晰的高精度成像，实现这一目标的前提条件是要有准确的速度参数。为了获得准确的地下模型信息，近年来许多方法得以发展，例如旅行时层析成像、立体层析、偏移速度分析、波形反演等。层析成像基于高频近似假设，其分辨率往往较低，而全波形反演可以利用地震走时、振幅、相位等多种波形信息，因此能够得到精确的高分辨率模型，而且理论上全波形反演适合所有类型的波，这使得全波形反演成为近年来人们研究的热点。
由于当前地震成像需要解决的问题变得越来越复杂化和精细化，单一纵波成像已经无法满足工业生产的需求，这主要原因是：1)单一纵波成像无法准确刻画地下构造信息；2)单一纵波成像结果不能充分反映地下岩石参数。为了充分估计地下岩石参数的性质，需提供弹性多波的成像信息。因此，为了弥补单独依靠纵波进行成像的不足，弹性波全波形反演显得十分重要，而且具有广阔的应用前景。
弹性波全波形反演的核心部分是波动方程数值模拟技术。波动方程正演求解的方法有很多，如有限差分法、有限元法、边界元法、伪谱法和谱元法等。相对于其他方法而言，有限差分法具有计算速度快、占用内存小，且可以方便地实现并行计算等优点，因此该方法被广泛地应用于波动方程数值模拟。目前，有限差分方法在空间差分精度优化方面已经趋于成熟，基于交错网格技术的有限差分方法可以实现空间任意偶数阶精度。然而，有限差分模拟中的时间频散误差对波场模拟结果有比较大的影响。该现象主要表现为模拟的波场相位超前和同相轴增多，为了提高模拟结果的精度，很多学者研究时间高阶有限差分压制波场模拟中的时间频散的方法，从而增大模拟波场的时间步长，在提高计算精度的同时又保证了计算效率。目前大多数时间高阶优化方法是基于声波方程的，对弹性波时间高阶优化方法的研究较少。相比于声波方程而言，弹性波方程可以模拟纵波(P)和横波(S)信息，并且纵横波的频散关系不同。有学者提出了一种基于S波频散关系的弹性波有限差分系数，但模拟P波的结果存在误差。后来发展了同时使用P和S波的频散关系来优化差分系数的方法，能达到减少时间频散的目的。此外，还有学者提出了一种基于速度-拟应力的弹性波方程，并给出了其优化后的时间四阶和六阶精度的差分系数。该方程既可以自然地实现P波和S波分离，又可以实现纵横波的高精度模拟。从计算效率和计算精度上看，该模拟方法优于传统方法。
然而，当前的弹性波全波形反演都是基于时间二阶有限差分算子的，当采用较大的时间步长并进行长时间的波场模拟时，会引入比较明显的时间频散，从而影响弹性波正演模拟和全波形反演的精度和效率。综上所述，现有的基于有限差分模拟的弹性波全波形反演方法，存在时间模拟精度低，计算效率低等问题。
因此，如何提供一种新的方案，其能够解决弹性波全波形反演存在的上述技术问题是本领域亟待解决的技术难题。
发明内容
本发明实施例提供一种有限差分模拟弹性波全波形反演方法，实现了高效的弹性波场模拟和纵横波速度模型反演，该方法包括：
根据纵横波速度模型，确定速度-拟应力弹性波动方程的纵横波差分参数；所述速度-拟应力弹性波动方程的纵横波差分参数，包括：纵横波有限差分离散格式，纵横波差分系数，纵横波稳定性条件，纵横波频散条件，初始条件，边界条件；
采用GPU共享内存优化策略，结合速度-拟应力弹性波动方程的纵横波差分参数，建立速度-拟应力弹性波动方程，进行时间四阶精度有限差分模拟，获得时间四阶精度有限差分模拟结果；
根据时间四阶精度有限差分模拟结果和速度-拟应力弹性波动方程，进行时间域弹性波全波形反演，确定纵横波速度模型。
本发明实施例还提供一种有限差分模拟弹性波全波形反演装置，包括：
速度-拟应力弹性波动方程的纵横波差分参数确定模块，用于根据纵横波速度模型，确定速度-拟应力弹性波动方程的纵横波差分参数；所述速度-拟应力弹性波动方程的纵横波差分参数，包括：纵横波有限差分离散格式，纵横波差分系数，纵横波稳定性条件，纵横波频散条件，初始条件，边界条件；
时间四阶精度有限差分模拟结果确定模块，用于采用GPU共享内存优化策略，结合速度-拟应力弹性波动方程的纵横波差分参数，建立速度-拟应力弹性波动方程，进行时间四阶精度有限差分模拟，获得时间四阶精度有限差分模拟结果；
纵横波速度模型确定模块，用于根据时间四阶精度有限差分模拟结果和速度-拟应力弹性波动方程，进行时间域弹性波全波形反演，确定纵横波速度模型。
本发明实施例还提供一种计算机设备，包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序，所述处理器执行所述计算机程序时实现上述有限差分模拟弹性波全波形反演方法。
本发明实施例还提供一种计算机可读存储介质，所述计算机可读存储介质存储有执行实现上述有限差分模拟弹性波全波形反演方法的计算机程序。
本发明实施例提供的一种有限差分模拟弹性波全波形反演方法和装置，采用速度-拟应力弹性波动方程，可以实现时间四阶精度有限差分模拟，并且可以对纵波和横波采用不同的空间差分阶数进行模拟，进而避免在一个波长内采样点不足或者过采样的问题。同时，该模拟方法可以实现纵横波的解耦，对基于波场分离类的反演算法具有优势。时间高阶差分允许使用较大的时间离散步长进行波场模拟，同时基于图形处理器(GPU)的共享内存的优化方式能进一步提高数据访问速度，从而加速弹性波动方程有限差分计算过程，实现了高效的弹性波场模拟和纵横波速度模型反演。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。在附图中：
图1为本发明实施例一种有限差分模拟弹性波全波形反演方法示意图。
图2为本发明实施例一种有限差分模拟弹性波全波形反演方法的时间四阶、空间任意偶数阶精度的差分离散格式示意图。
图3a-图3c为本发明实施例一种有限差分模拟弹性波全波形反演方法的常规时间二阶和本发明实施例的时间四阶有限差分模拟的波场切片对比图。
图4为本发明实施例一种有限差分模拟弹性波全波形反演方法的波剖面对比示意图。
图5a-图5b为本发明实施例一种有限差分模拟弹性波全波形反演方法的GPU计算及其共享内存的优化示意图。
图6a-图6d为本发明实施例一种有限差分模拟弹性波全波形反演方法的基于GPU共享内存的加速比测试结果的示意图。
图7a-图7b为本发明实施例一种有限差分模拟弹性波全波形反演方法的模型大小对基于GPU共享内存的正演模拟计算效率的影响图。
图8a-图8h为本发明实施例一种有限差分模拟弹性波全波形反演方法的时间四阶弹性波全波形反演运用到Overthrust模型反演结果与常规方法反演结果的对比示意图。
图9a-图9b为本发明实施例一种有限差分模拟弹性波全波形反演方法的利用图8a-8图h的结果绘制的局部放大显示示意图。
图10为运行本发明实施的一种有限差分模拟弹性波全波形反演方法的计算机装置示意图。
图11为本发明实施例一种有限差分模拟弹性波全波形反演装置示意图。
具体实施方式
为使本发明实施例的目的、技术方案和优点更加清楚明白，下面结合附图对本发明实施例做进一步详细说明。在此，本发明的示意性实施例及其说明用于解释本发明，但并不作为对本发明的限定。
图1为本发明实施例一种有限差分模拟弹性波全波形反演方法示意图，如图1所示，本发明实施例提供一种有限差分模拟弹性波全波形反演方法，实现了纵横波速度模型的高效模拟和反演，该方法包括：
步骤101：根据纵横波速度模型，确定速度-拟应力弹性波动方程的纵横波差分参数；所述速度-拟应力弹性波动方程的纵横波差分参数，包括：纵横波有限差分离散格式，纵横波差分系数，纵横波稳定性条件，纵横波频散条件，初始条件，边界条件；
步骤102：采用GPU共享内存优化策略，结合速度-拟应力弹性波动方程的纵横波差分参数，建立速度-拟应力弹性波动方程，进行时间四阶精度有限差分模拟，获得时间四阶精度有限差分模拟结果；
步骤103：根据时间四阶精度有限差分模拟结果和速度-拟应力弹性波动方程，进行时间域弹性波全波形反演，确定纵横波速度模型。
本发明实施例提供的一种有限差分模拟弹性波全波形反演方法，采用速度-拟应力弹性波动方程，可以实现时间四阶精度有限差分模拟，并且可以对纵波和横波采用不同的空间差分阶数进行模拟，进而避免在一个波长内采样点不足或者过采样的问题。同时，该模拟方法可以实现纵横波的解耦，对基于波场分离类的反演算法具有优势。时间高阶差分允许使用较大的时间离散步长进行波场模拟，同时基于图形处理器(GPU)的共享内存的优化方式能进一步提高数据访问速度，从而加速弹性波动方程有限差分计算过程，实现了高效的弹性波场模拟和纵横波速度模型反演。
本发明实施例提供的一种有限差分模拟弹性波全波形反演方法，可以包括：
根据纵横波速度模型，确定速度-拟应力弹性波动方程的纵横波差分参数；所述速度-拟应力弹性波动方程的纵横波差分参数，包括：纵横波有限差分离散格式，纵横波差分系数，纵横波稳定性条件，纵横波频散条件，初始条件，边界条件；采用GPU共享内存优化策略，结合速度-拟应力弹性波动方程的纵横波差分参数，建立速度-拟应力弹性波动方程，进行时间四阶精度有限差分模拟，获得时间四阶精度有限差分模拟结果；根据时间四阶精度有限差分模拟结果和速度-拟应力弹性波动方程，进行时间域弹性波全波形反演，确定纵横波速度模型。
在本发明实施例的一个实例中，为了实现时间四阶有限差分模拟，本发明实施例建立了一种速度-拟应力弹性波动方程，通过有限差分模拟方法进行波场外推。该方程的表达式为：
Figure BDA0002768009270000061
式(1)中，x和z表示二维(2D)空间坐标系下的坐标，t表示时间坐标；λ和μ表示拉梅常数，ρ表示密度；λ+2μ＝ρα2，μ＝ρβ2，α和β分别表示纵横波速度；vx和vz分别表示质点的水平和垂直振动速度，τ111，τ11，τ12，τ22表示质点的应力分量；
Figure BDA0002768009270000062
表示时间一阶偏导数算子，
Figure BDA0002768009270000063
和
Figure BDA0002768009270000064
表示沿着x和z的一阶空间偏导数。式(1)中的空间偏导数的波数域响应式可表示为：
Figure BDA0002768009270000065
式(2)中，
Figure BDA0002768009270000066
为虚数单位，kr(r＝x,z)表示r方向的波数，
Figure BDA0002768009270000067
表示波数域的波场，F-1表示反傅里叶变换；Ll的表达式为：
Lα＝sinc(αkΔt/2),Lβ＝sinc(βkΔt/2) (3)
式(3)中，Δt是采用二阶中心差分近似式(1)中时间偏导数所采用的时间步长，并且波数
Figure BDA0002768009270000068
式(1)给出了全波形反演中所使用的弹性波方程，式(2)和(3)给出了空间偏导数的波数域响应。式(1)中的时间导数可以通过交错网格有限差分来求解。
在具体实施本发明实施例提供的一种有限差分模拟弹性波全波形反演方法时，一个实施例中，按照如下方式，确定纵横波有限差分离散格式：
Figure BDA0002768009270000069
Figure BDA00027680092700000610
其中，
Figure BDA00027680092700000611
为纵波有限差分离散格式；
Figure BDA00027680092700000612
为横波有限差分离散格式；α为纵波速度；β为横波速度；u为待差分波场变量，r,φ∈{x,z}，并且r≠φ；hr表示沿着空间r方向的空间网格步长，Nr表示沿空间r方向的单边差分格式长度；u中的下标r表示沿着差分方向的网格点，u中的上标φ表示非差分轴向的网格点；c表示差分系数，其上标之一r表示差分沿着r方向，p表示差分系数由纵波速度计算，其另一上标s表示差分系数由横波速度计算。
前述提到的确定纵横波有限差分离散格式的表达式为举例说明，本领域技术人员可以理解，在实施时还可以根据需要对上述公式进行一定形式的变形和添加其它的参数或数据，或者提供其它的具体公式，这些变化例均应落入本发明的保护范围。
在本发明实施例的一个实例中，由于
Figure BDA0002768009270000071
和
Figure BDA0002768009270000072
分别与纵波、横波的速度相关，所以需要两套差分系数来求解该方程。以
Figure BDA0002768009270000073
为例，其基于时间四阶交错网格有限差分的离散格式如上述公式(4)所示。
式(4)所描述的差分模板如图2所示，待计算点的坐标为(x,z)，这里以
Figure BDA0002768009270000074
详细说明，
Figure BDA0002768009270000075
表示波场变量位于坐标点(x+0.5hx,z+hz)。在图2中：
Figure BDA0002768009270000076
表示关于坐标轴r方向的偏导数符号r，在二维笛卡尔xoz坐标系下，r＝x或z，u表示波场变量，Nr表示空间差分阶数，hr表示空间沿着r方向的网格间隔，空心的圆圈表示参与差分计算并分布在r轴上网格节点上的波场值，四个实心的小圆圈表示参与差分计算但不位于r轴上网格节点的波场值。
在具体实施本发明实施例提供的一种有限差分模拟弹性波全波形反演方法时，一个实施例中，按照如下方式，确定纵横波差分系数：
Figure BDA0002768009270000077
其中，c表示差分系数；在网格比γr,p＝αΔt/hr时，通过上式获得纵波差分系数，其中Δt表示二阶中心差分近似时间偏导数所使用的时间步长；相应地，将上式中γr,p相应替换为网格比γr,S＝βΔt/hr即可获得横波差分系数；α为纵波速度；β为横波速度；m为轴向差分系数次序索引，n为上式求和引入的变量。
前述提到的确定纵横波差分系数的表达式为举例说明，本领域技术人员可以理解，在实施时还可以根据需要对上述公式进行一定形式的变形和添加其它的参数或数据，或者提供其它的具体公式，这些变化例均应落入本发明的保护范围。
差分系数可以通过泰勒展开或者最小二乘优化的方法求取。需要注意的是，c中的上标p表示了该系数是用来求解与纵波相关的导数的。此处给出了使用泰勒展开方法求取的与纵波相关的导数的差分系数，其表达式如上述公式(5)所示。
在公式(5)中，γr,p＝αΔt/hr。在式(1)中，与横波相关的空间导数
Figure BDA0002768009270000081
和
Figure BDA0002768009270000082
可以用相同的差分模板求得相应的差分系数。通过对比纵横波的差分系数可知，横波相关的差分系数使用与横波相关的γr,S＝βΔt/hr，其余的表达式与式(5)一致。由于对于式(1)采用了与纵横波相关的两套独立的差分系数，可以采用不同的差分精度来分别求取纵波和横波相关的空间导数，避免常规弹性波方程数值模拟中采用相同差分系数而造成的采样不足或者过采样问题。
在具体实施本发明实施例提供的一种有限差分模拟弹性波全波形反演方法时，一个实施例中，按照如下方式，确定纵波稳定性条件：
Figure BDA0002768009270000083
其中，
Figure BDA0002768009270000084
r＝x,z；α为纵波速度；hx为x轴向网格间距，hz为z轴向网格间距；sx和sz的表达式见sr；
按照如下方式，确定横波稳定性条件：
Figure BDA0002768009270000085
其中，
Figure BDA0002768009270000086
r＝x,z；β为横波速度；hx为x轴向网格间距，hz为z轴向网格间距；sx和sz的表达式见sr。
前述提到的确定纵波稳定性条件和横波稳定性条件的表达式为举例说明，本领域技术人员可以理解，在实施时还可以根据需要对上述公式进行一定形式的变形和添加其它的参数或数据，或者提供其它的具体公式，这些变化例均应落入本发明的保护范围。
为了保证差分计算的稳定性，纵波的差分稳定性条件应该满足如上述公式(6)-1；
在上述公式(6)-1中，
Figure BDA0002768009270000091
r＝x,z。横波的差分格式稳定性条件是在式(6)-2中使用
Figure BDA0002768009270000092
代替
Figure BDA0002768009270000093
即可求取。
在具体实施本发明实施例提供的一种有限差分模拟弹性波全波形反演方法时，一个实施例中，按照如下方式，确定纵波频散条件：
Figure BDA0002768009270000094
其中，α为纵波速度；max{·}表示选取最大值，f0是子波的主频，N是某一差分精度下所允许的最小波长最少采样点数；
按照如下方式，确定横波频散条件：
Figure BDA0002768009270000095
如下方式，β为横波速度；max{·}表示选取最大值，f0是子波的主频，N是某一差分精度下所允许的最小波长的最少采样点数。
前述提到的确定纵波频散条件和横波频散条件的表达式为举例说明，本领域技术人员可以理解，在实施时还可以根据需要对上述公式进行一定形式的变形和添加其它的参数或数据，或者提供其它的具体公式，这些变化例均应落入本发明的保护范围。
为了避免明显的空间数值频散，需要保证每个波长有足够的采样点。通常，空间二阶精度需要保证最小波长至少有10个采样点，空间四阶精度需要保证最小波长至少有5个采样点。
在具体实施本发明实施例提供的一种有限差分模拟弹性波全波形反演方法时，一个实施例中，按照如下方式，确定初始条件：
Figure BDA0002768009270000096
其中，v＝[vx vz]T；τ＝[τ111 τ11 τ12 τ22]Tτ111；τ11，τ12，τ22表示质点的应力分量；s(r0)表示震源函项，f(t)表示震源函数。
前述提到的确定初始条件的表达式为举例说明，本领域技术人员可以理解，在实施时还可以根据需要对上述公式进行一定形式的变形和添加其它的参数或数据，或者提供其它的具体公式，这些变化例均应落入本发明的保护范围。
为了模拟弹性波动方程，需要初始条件和边界条件。其中初始条件表示上述公式(8)。
在具体实施本发明实施例提供的一种有限差分模拟弹性波全波形反演方法时，一个实施例中，按照如下方式，确定边界条件：
Figure BDA0002768009270000101
Figure BDA0002768009270000102
其中，r∈{x,z}，δ是边界条件的层数，R是PML边界反射系数，其取值范围为10-6～10-8；α为纵波速度；β为横波速度；在上述边界条件中将vx分裂成x方向和z方向两部分，且将x方向的部分进一步分裂成纵波分量
Figure BDA0002768009270000103
和横波分量
Figure BDA0002768009270000104
即
Figure BDA0002768009270000105
前述提到的确定边界条件的表达式为举例说明，本领域技术人员可以理解，在实施时还可以根据需要对上述公式进行一定形式的变形和添加其它的参数或数据，或者提供其它的具体公式，这些变化例均应落入本发明的保护范围。
吸收边界设置在模型的四周，用来吸收截断边界反射，消除边界反射对反演的影响。在本发明中，采用完美匹配层(PML)来消除边界反射。以式(1)中的第一项为例，使用PML边界条件的方程可以表示上述公式(9)和公式(10)。
基于上述时间四阶有限差分模拟算法的描述，使用模拟结果来说明该算法的模拟精度。采用一个两层的速度模型，密度为常密度2000千克每立方米。网格大小为800×400，网格间距为10米。模型的第一层对应的纵波和横波速度分别为3000米每秒和1800米每秒，第二层分别为4000米每秒和2400米每秒。速度界面位于地下2千米处。采用横波的空间差分算子长度为10阶，纵波的差分算子长度为6阶。图3a-图3c为本发明实施例一种有限差分模拟弹性波全波形反演方法的常规时间二阶和本发明实施例的时间四阶有限差分模拟的波场切片对比图，表示时刻为1.68秒模拟的波场切片示意图。图3a是常规时间二阶精度并采用时间步长为0.6毫秒的模拟结果，图3b是常规时间二阶精度并采用时间步长为1.2毫秒的模拟结果，图3c是时间四阶精度并采用时间步长为1.4毫秒的模拟结果。从图3a-图3c中可知常规时间二阶精度模拟方法在采用较大的时间步长时，会存在明显的时间频散现象，而本发明的时间四阶模拟方法在使用较大的时间步长的条件下，仍然可以得到准确的模拟结果。图4为本发明实施例一种有限差分模拟弹性波全波形反演方法的波剖面对比示意图。T2表示时间二阶差分精度，T4表示时间四阶差分精度。
图4是来自图3a-图3c集合的一组曲线，由图可知时间四阶算法在使用时间离散步长为1.4毫秒时，得到的结果与时间二阶精度方法使用0.6毫米时间步长的结果一致。采用大的时间步长模拟，可以在保证计算精度的前提下，提高正演模拟效率。
在具体实施本发明实施例提供的一种有限差分模拟弹性波全波形反演方法时，一个实施例中，所述GPU共享内存优化策略，包括：
从全局内存中将时间四阶精度有限差分模拟所需的数据拷贝至每个GPU块的共享内存中；
为每个GPU块沿着差分方向扩出差分长度的区域，使差分方向的每个点满足时间四阶精度差分格式；
根据时间四阶精度有限差分格式，从共享内存中读取时间四阶精度有限差分模拟数据进行计算。
上述步骤介绍了速度-拟应力弹性波动方程的离散格式、差分系数、稳定性及频散条件、初始和边界条件。此处描述其GPU的实现及其共享内存的优化策略。GPU的计算能力优越，基于GPU的编程实现包括：CPU控制数据的输入输出以及简单的逻辑判断，GPU承担核心的差分计算。图5a-图5b为本发明实施例一种有限差分模拟弹性波全波形反演方法的GPU计算及其共享内存的优化示意图。
图5a展示了常规GPU全局内存的使用，采用图2所示的差分格式计算弹性波动方程(1)中的空间偏导数。计算每一个空间偏导数
Figure BDA0002768009270000111
需要的全局内存为：
C1＝{nxz×(2Nr+4)}g (11)
式(11)中，nxz＝nx×nz，nx和nz表示沿着x和z方向的网格点数；式(11)中下标g表示读取的数据来自于GPU全局内存。尽管GPU计算能力很强，但是频繁地从全局内存读取数据会影响计算效率。为此，采用GPU共享内存的方式来减少频繁的全局内存数据读取。由于GPU的共享内存小但其数据访问更高效，所以合理使用该内存可以提升效率。GPU的共享内存是GPU块(Block)私有的，在同一个GPU数据块中，所有的线程共享该内存。图5b展示了GPU共享内存数据访问策略的实现方法。第一步，从全局内存中拷贝数据到每个GPU块的共享内存中。在拷贝数据过程中，为了使得导数方向的每个点满足图1的差分格式，需要为每个GPU块沿着导数方向扩出差分长度的区域。第二步，根据差分格式计算波场值，此时多次数据读取是通过访问共享内存完成的。在GPU共享内存的数据访问方式下，有限差分计算一个空间偏导数的数据访问次数为：
Figure BDA0002768009270000121
式(12)中，Bxz＝Bx×Bz，下标s表示读取的数据来自于GPU共享内存。
图6a-图6d为本发明实施例一种有限差分模拟弹性波全波形反演方法的基于GPU共享内存的加速比测试结果的示意图；图6a是使用GPU的块大小为32×32的不同差分阶数的弹性波正演模拟的耗时结果；图6b是使用GPU的块大小为16×16的不同差分阶数的弹性波正演模拟的耗时结果；图6c是使用GPU的块大小为8×8的不同差分阶数的弹性波正演模拟的耗时结果；图6d是三种GPU块大小情况下，不同差分阶数对应的弹性波正演模拟的加速度比结果。图7a-图7b为本发明实施例一种有限差分模拟弹性波全波形反演方法的模型大小对基于GPU共享内存的正演模拟计算效率的影响图；图7a是弹性波正演模拟的耗时结果；图7b是其加速比曲线。在图中，所用的GPU块的大小是16×16，纵横波的空间差分阶数分别是6和10；图6a-图6d和图7a-图7b测试了基于GPU共享内存数据访问策略下，不同差分阶数，不同的GPU块的大小(块所包含的核心数)，以及不同大小的正演模型所对应的模拟时长及其加速比情况。其中图6a-图6c为不同GPU块大小下数值模拟计算耗时，图6d是由这三种耗时计算出的加速比。图6a-图6d中结果证实，此类加速比与差分阶数有着明显的关系，随着空间阶数的增加，未优化的基于GPU全局内存的算法需要多次从GPU全局内存中读取数据并计算差分近似，而优化后的GPU共享内存算法则是从GPU的共享内存中读取数据并进行差分计算，因此差分算子越长，所带来的计算效率提升越明显。图6a-图6d说明了基于GPU共享内存的有限差分算法与模型的大小没有明显的依赖性关系，这种优化算法仅与空间差分算子长度有关，当空间差分长度一定时，该加速比趋于一定值。
在具体实施本发明实施例提供的一种有限差分模拟弹性波全波形反演方法时，一个实施例中，按照如下方式，建立速度-拟应力弹性波动方程，也即是前述的式(1)：
Figure BDA0002768009270000131
其中：x和z表示二维(2D)空间坐标系下的坐标，t表示时间坐标；λ和μ表示拉梅常数，ρ表示密度；λ+2μ＝ρα2，μ＝ρβ2，α和β分别表示纵横波速度；vx和vz分别表示质点的水平和垂直振动速度，τ111，τ11，τ12，τ22表示质点的应力分量；
Figure BDA0002768009270000132
表示时间一阶偏导数算子，
Figure BDA0002768009270000133
和
Figure BDA0002768009270000134
表示沿着x和z的一阶空间偏导数。
前述提到的建立速度-拟应力弹性波动方程的表达式为举例说明，本领域技术人员可以理解，在实施时还可以根据需要对上述公式进行一定形式的变形和添加其它的参数或数据，或者提供其它的具体公式，这些变化例均应落入本发明的保护范围。
在具体实施本发明实施例提供的一种有限差分模拟弹性波全波形反演方法时，一个实施例中，根据时间四阶精度有限差分模拟结果和速度-拟应力弹性波动方程，进行时间域弹性波全波形反演，确定纵横波速度模型，包括：
匹配观测数据和合成数据，构建目标函数；
根据速度-拟应力弹性波动方程和目标函数，采用拉格朗日伴随方法获得目标函数对模型参数的梯度方程；
根据时间四阶精度有限差分模拟结果进行正演和反传计算，获得震源和伴随波场变量；
利用目标函数对模型参数的梯度方程和震源、伴随波场变量获取模型参数梯度；
采用共轭梯度方法求取更新方向，对模型参数进行迭代更新，直至残差达到设定值，输出纵横波速度模型反演参数，确定纵横波速度模型。
在实现了基于速度-拟应力弹性波动方程的时间四阶有限差分模拟的基础上，进一步发展全波形反演算法。其原理是构建最小二乘目标函数来匹配观测数据和合成数据，通过迭代的方法，修改模型并逐步减小数据残差，最终求取反演的模型参数。为了增强算法的鲁棒性，采用一种基于卷积类型的目标函数。
在具体实施本发明实施例提供的一种有限差分模拟弹性波全波形反演方法时，一个实施例中，按照如下方式，构建目标函数：
Figure BDA0002768009270000141
其中，d是观测数据，u是合成数据，xr是检波点的位置，*是时间卷积算子，xref表示要提取参考道的位置；
Figure BDA0002768009270000142
表示欧氏范数，m＝(α，β)是模型参数。
前述提到的构建目标函数的表达式为举例说明，本领域技术人员可以理解，在实施时还可以根据需要对上述公式进行一定形式的变形和添加其它的参数或数据，或者提供其它的具体公式，这些变化例均应落入本发明的保护范围。
由于卷积运算的滤波作用，使用由低频到高频的子波模拟，从而实现多尺度全波形反演。通常，在式(13)中的速度分量作为观测数据，即u＝(vx,vz)。
在具体实施本发明实施例提供的一种有限差分模拟弹性波全波形反演方法时，一个实施例中，按照如下方式，获得目标函数对模型参数的梯度方程：
Figure BDA0002768009270000143
Figure BDA0002768009270000144
其中，α为纵波速度；β为横波速度；<,>t表示关于时间的内积，
Figure BDA0002768009270000145
表示伴随波场变量；ρ表示密度。
前述提到的获得目标函数对模型参数的梯度方程的表达式为举例说明，本领域技术人员可以理解，在实施时还可以根据需要对上述公式进行一定形式的变形和添加其它的参数或数据，或者提供其它的具体公式，这些变化例均应落入本发明的保护范围。
为了构建速度-拟应力弹性波动方程的数据对模型的梯度，使用拉格朗日伴随方法推导目标函数对模型参数的梯度。考虑到编程实现的方便性，推导正传波动方程和反传波动方程具有相同形式时目标函数对模型参数的梯度表达式如公式(14)所示。
在公式(14)和式(15)中，<,>t表示关于时间的内积，
Figure BDA0002768009270000151
表示伴随波场变量。通过反传检波点波场可以求取伴随波场变量。基于反传波场的波动方程仍然采用式(1)，但是使用新的震源，其定义如下：
Figure BDA0002768009270000152
式(16)中的
Figure BDA0002768009270000153
表示互相关运算。
在具体实施本发明实施例提供的一种有限差分模拟弹性波全波形反演方法时，一个实施例中，按照如下方式，求取更新方向：
Figure BDA0002768009270000154
其中，下标k表示迭代次数；上标T表示矩阵转置，并且
Figure BDA0002768009270000155
sk为共轭因子；m为模型参数向量，在此式中为纵横波速度组成的向量。
前述提到的求取更新方向的表达式为举例说明，本领域技术人员可以理解，在实施时还可以根据需要对上述公式进行一定形式的变形和添加其它的参数或数据，或者提供其它的具体公式，这些变化例均应落入本发明的保护范围。
在具体实施本发明实施例提供的一种有限差分模拟弹性波全波形反演方法时，一个实施例中，按照如下方式，进行迭代更新：
mk+1＝mk+tkΔmk+1 (18)
其中，tk表示选取的更新步长；下标k表示迭代次数；m为模型参数向量，在此式中为纵横波速度组成的向量。前述提到的进行迭代更新的表达式为举例说明，本领域技术人员可以理解，在实施时还可以根据需要对上述公式进行一定形式的变形和添加其它的参数或数据，或者提供其它的具体公式，这些变化例均应落入本发明的保护范围。
通过正演和反传计算出相应的波场变量，利用式(14)求取梯度，采用公式(17)的共轭梯度方法求取更新方向。
图8a-图8h为本发明实施例一种有限差分模拟弹性波全波形反演方法的时间四阶弹性波全波形反演运用到Overthrust模型反演结果与常规方法反演结果的对比示意图。
图8a和图8b是真实的纵波和横波的Overthrust速度模型；图8c和图8d是初始的纵横波速度模型；该模型对应的密度是常密度模型。该模型的网格数为800×170，网格间距是12.5米。采用式(1)合成30炮观测数据进行反演，其中炮间距324米，炮点深度12.5米，每一炮有800个检波点接收数据。在时间四阶有限差分数值模拟过程中，分别采用8阶和12阶空间差分阶数来近似纵横波对应的偏导数项。
图8e和图8f是时间四阶精度反演的纵横波速度模型；图8g和图8h是传统的时间二阶精度反演的纵横波速度模型。
图8e-图8h是反演结果示意图。在该反演过程中，采用了主频为4赫兹、7赫兹、10赫兹、14赫兹和19赫兹的雷克子波进行多尺度反演，其中每个频带迭代15次。为了减少波场存储，采用有效边界在反传过程中重构波场。图8e和图8f是基于时间四阶差分法反演得到的结果，其中采用的时间离散步长是1.3毫秒，为稳定性条件所允许时间步长的99％；图8e和图8f是基于时间二阶精度差分法反演得到的结果，其中采用的时间离散步长是1.1毫秒，为稳定性条件所允许时间步长的98％。其中图8e和图8g是纵波速度反演结果，图8f和8h是横波速度反演结果。从反演结果上看，尽管基于时间二阶差分算子的反演结果和基于时间四阶差分算子的反演结果比较接近，但是基于时间四阶差分算子的反演结果在层内比较光滑，层内速度分布比较均匀。
图9a-图9b为本发明实施例一种有限差分模拟弹性波全波形反演方法的利用图8a-8图h的结果绘制的局部放大显示示意图。图9a和图9b分别是位于z＝0.12到z＝0.8公里、x＝2到x＝6公里范围内的纵横波反演结果的局部显示对比。在图9a和图9b中，从左往右依次为真实的速度结果，时间二阶差分算子的反演结果和时间四阶差分算子的反演结果。
图9a-图9b是图8a-图8h中位于z＝0.12到z＝0.8公里、x＝2到x＝6公里范围内的反演结果的局部显示示意图。图9a是纵波速度反演结果，图9b是横波速度反演结果。从局部结果看到，两种方法可以取得比较好的结果，但是基于时间四阶差分算子的反演方法耗时10.83小时，基于时间二阶差分算子的反演方法耗时12.42小时。从运行时间上可以看出，基于时间四阶差分算子的反演方法允许使用较大的时间步长，因此计算效率比常规时间二阶差分算子反演方法的效率要高。
在具体实施时，本发明实施例包括如下三个步骤：1)时间四阶精度有限差分弹性波模拟实现方法；2)基于GPU加速及其共享内存的优化策略；3)基于速度-拟应力弹性波动方程的全波形反演技术。
所述步骤1)中，采用速度-拟应力弹性波动方程，使用有限差分方法进行波场模拟。给定差分离散格式，利用泰勒展开求取时间四阶精度、空间任意偶数阶精度的有限差分系数，在满足纵横波稳定性条件开展波场模拟。在该步骤中，纵横波使用独立的差分系数，并针对纵横波设置不同的空间精度进行波场模拟，避免常规弹性波方程数值模拟中采用相同的差分系数而造成的采样不足或者过采样的问题。
所述步骤2)中，主要对全波形反演算法的效率进行优化。相比于传统的中央处理器(CPU)，GPU具有更高的计算效率。因此采用CPU和GPU协同完成高效的计算。CPU主要负责逻辑控制，GPU负责核心计算。通过将计算区域映射到GPU的不同块中，在每一个GPU块中可以同时启用一定数量的线程，多线程以及多GPU的同时运行，完成了关键的差分计算部分。此外，传统的GPU计算中没有使用GPU的共享内存，为了充分利用GPU的资源，采用GPU的共享内存进行优化。共享内存具有小内存，高速读取等特点，对于一些差分计算中反复读取的数值，可以使用共享内存的变量类型，从而达到高效计算的目的。
所述步骤3)中，主要阐述了基于卷积型目标函数的反演算法。通过使用该目标函数，可以自然地实现多尺度反演，同时可以削弱子波对反演结果的影响。该目标函数是分别从观测数据和合成数据各选一个参考道，然后用数据和不同类型的参考道的卷积(观测数据与合成数据的参考道的卷积以及合成数据与观测数据的卷积)的差来衡量数据的拟合程度。其中参考道的选取规则是：尽量选取炮点附近的具有较高信噪比的道作为参考道，对于信噪比较差的资料，可以选取炮点附近相邻的数道，较正到同一时间并叠加形成参考道。由于采用了步骤1)中描述的新的速度-拟应力弹性波动方程，推导了基于该弹性波动方程的模型参数的梯度表达式。通过求取模型参数的梯度，采用共轭梯度法来更新模型。
本发明由于采取以上技术方案，其具有以下优点：1)本发明采用速度-拟应力弹性波动方程，可以实现时间四阶精度的有限差分模拟，并且可以对纵横波采用不同的空间差分阶数进行波场模拟，进而避免在一个波长内采样点不足或者过采样的问题。同时，该模拟方法可以实现纵横波的解耦，对基于波场分离类的反演算法具有独特的优势。2)本发明采用GPU共享内存对有限差分计算进行优化，进一步提高弹性波动方程有限差分数值模拟的计算效率。3)本发明在弹性波全波形反演中首次引入速度-拟应力方程作为控制方程，并推导出该方程对应的梯度表达式，形成了一套基于该波动方程的全波形反演算法和流程。4)在3)中推导的梯度表达，允许正传算子和反传算子使用相同的控制方程和不同的震源项，从而使得编程简单并易于实现。基于以上优点，本发明可以广泛应用于弹性波参数建模中。
图10为运行本发明实施的一种有限差分模拟弹性波全波形反演方法的计算机装置示意图。如图10所示，本发明实施例还提供一种计算机设备，包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序，所述处理器执行所述计算机程序时实现上述有限差分模拟弹性波全波形反演方法。
本发明实施例还提供一种计算机可读存储介质，所述计算机可读存储介质存储有执行实现上述有限差分模拟弹性波全波形反演方法的计算机程序。
本发明实施例中还提供了一种有限差分模拟弹性波全波形反演装置，如下面的实施例所述。由于该装置解决问题的原理与一种有限差分模拟弹性波全波形反演方法相似，因此该装置的实施可以参见一种有限差分模拟弹性波全波形反演方法的实施，重复之处不再赘述。
图11为本发明实施例一种有限差分模拟弹性波全波形反演装置示意图。如图11所示，本发明实施例还提供一种有限差分模拟弹性波全波形反演装置，包括：
速度-拟应力弹性波动方程的纵横波差分参数确定模块1101，用于根据纵横波速度模型，确定速度-拟应力弹性波动方程的纵横波差分参数；所述速度-拟应力弹性波动方程的纵横波差分参数，包括：纵横波有限差分离散格式，纵横波差分系数，纵横波稳定性条件，纵横波频散条件，初始条件，边界条件；
时间四阶精度有限差分模拟结果确定模块1102，用于采用GPU共享内存优化策略，结合速度-拟应力弹性波动方程的纵横波差分参数，建立速度-拟应力弹性波动方程，进行时间四阶精度有限差分模拟，获得时间四阶精度有限差分模拟结果；
纵横波速度模型确定模块1103，用于根据时间四阶精度有限差分模拟结果和速度-拟应力弹性波动方程，进行时间域弹性波全波形反演，确定纵横波速度模型。
在具体实施本发明实施例提供的一种有限差分模拟弹性波全波形反演装置时，一个实施例中，速度-拟应力弹性波动方程的纵横波差分参数确定模块，用于按照如下方式，确定纵横波有限差分离散格式：
Figure BDA0002768009270000191
Figure BDA0002768009270000192
其中，
Figure BDA0002768009270000193
为纵波有限差分离散格式；
Figure BDA0002768009270000194
为横波有限差分离散格式；α为纵波速度；β为横波速度；u为待差分波场变量，r,φ∈{x,z}，并且r≠φ；hr表示沿着空间r方向的空间网格步长，Nr表示沿空间r方向的单边差分格式长度；u中的下标r表示沿着差分方向的网格点，u中的上标φ表示非差分轴向的网格点；c表示差分系数，其上标之一r表示差分沿着r方向，p表示差分系数由纵波速度计算，其另一上标s表示差分系数由横波速度计算。
在具体实施本发明实施例提供的一种有限差分模拟弹性波全波形反演装置时，一个实施例中，速度-拟应力弹性波动方程的纵横波差分参数确定模块，用于按照如下方式，确定纵横波差分系数：
Figure BDA0002768009270000195
其中，c表示差分系数；在网格比γr,p＝αΔt/hr时，通过上式获得纵波差分系数，其中Δt表示二阶中心差分近似时间偏导数所使用的时间步长；相应地，将上式中γr,p相应替换为网格比γr,S＝βΔt/hr即可获得横波差分系数；α为纵波速度；β为横波速度；m为轴向差分系数次序索引，n为上式求和引入的变量。
在具体实施本发明实施例提供的一种有限差分模拟弹性波全波形反演装置时，一个实施例中，速度-拟应力弹性波动方程的纵横波差分参数确定模块，用于按照如下方式，确定纵波稳定性条件：
Figure BDA0002768009270000201
其中，
Figure BDA0002768009270000202
r＝x,z；α为纵波速度；hx为x轴向网格间距，hz为z轴向网格间距；
速度-拟应力弹性波动方程的纵横波差分参数确定模块，还用于按照如下方式，确定横波稳定性条件：
Figure BDA0002768009270000203
其中，
Figure BDA0002768009270000204
r＝x,z；β为横波速度；hx为x轴向网格间距，hz为z轴向网格间距。
在具体实施本发明实施例提供的一种有限差分模拟弹性波全波形反演装置时，一个实施例中，速度-拟应力弹性波动方程的纵横波差分参数确定模块，用于按照如下方式，确定纵波频散条件：
Figure BDA0002768009270000205
其中，α为纵波速度；max{·}表示选取最大值，f0是子波的主频，N是某一差分精度下所允许的最小波长的最少采样点数；
速度-拟应力弹性波动方程的纵横波差分参数确定模块，用于按照如下方式，确定横波频散条件：
Figure BDA0002768009270000206
如下方式，β为横波速度；max{·}表示选取最大值，f0是子波的主频，N是某一差分精度下所允许的最小波长的最少采样点数。
在具体实施本发明实施例提供的一种有限差分模拟弹性波全波形反演装置时，一个实施例中，速度-拟应力弹性波动方程的纵横波差分参数确定模块，用于按照如下方式，确定初始条件：
Figure BDA0002768009270000211
其中，v＝[vx vz]T；τ＝[τ111 τ11 τ12 τ22]Tτ111；τ11，τ12，τ22表示质点的应力分量；s(r0)表示震源函项，f(t)表示震源函数。
在具体实施本发明实施例提供的一种有限差分模拟弹性波全波形反演装置时，一个实施例中，速度-拟应力弹性波动方程的纵横波差分参数确定模块，用于按照如下方式，确定边界条件：
Figure BDA0002768009270000212
Figure BDA0002768009270000213
其中，r∈{x,z}，δ是边界条件的层数，R是边界反射系数，其取值范围为10-6～10-8；α为纵波速度；β为横波速度；在上述边界条件中将vx分裂成x方向和z方向两部分，且将x方向的部分进一步分裂成纵波分量
Figure BDA0002768009270000214
和横波分量
Figure BDA0002768009270000215
即
Figure BDA0002768009270000216
在具体实施本发明实施例提供的一种有限差分模拟弹性波全波形反演装置时，一个实施例中，时间四阶精度有限差分模拟结果确定模块，具体用于：
从全局内存中将时间四阶精度有限差分模拟所需的数据拷贝至每个GPU块的共享内存中；
为每个GPU块沿着差分方向扩出差分长度的区域，使差分方向的每个点满足时间四阶精度差分格式；
根据时间四阶精度有限差分格式，从共享内存中读取时间四阶精度有限差分模拟数据进行计算。
在具体实施本发明实施例提供的一种有限差分模拟弹性波全波形反演装置时，一个实施例中，时间四阶精度有限差分模拟结果确定模块，用于按照如下方式，建立速度-拟应力弹性波动方程：
Figure BDA0002768009270000221
其中：x和z表示二维(2D)空间坐标系下的坐标，t表示时间坐标；λ和μ表示拉梅常数，ρ表示密度；λ+2μ＝ρα2，μ＝ρβ2，α和β分别表示纵横波速度；vx和vz分别表示质点的水平和垂直振动速度，τ111，τ11，τ12，τ22表示质点的应力分量；
Figure BDA0002768009270000222
表示时间一阶偏导数算子，
Figure BDA0002768009270000223
和
Figure BDA0002768009270000224
表示沿着x和z的一阶空间偏导数。
在具体实施本发明实施例提供的一种有限差分模拟弹性波全波形反演装置时，一个实施例中，纵横波速度模型确定模块，用于：
匹配观测数据和合成数据，构建目标函数；
根据速度-拟应力弹性波动方程和目标函数，采用拉格朗日伴随方法获得目标函数对模型参数的梯度方程；
根据时间四阶精度有限差分模拟结果进行正演和反传计算，获得震源和伴随波场变量；
利用目标函数对模型参数的梯度方程和震源、伴随波场变量获取模型参数梯度；
采用共轭梯度方法求取更新方向，对模型参数进行迭代更新，直至残差达到设定值，输出纵横波速度模型反演参数，确定纵横波速度模型。
在具体实施本发明实施例提供的一种有限差分模拟弹性波全波形反演装置时，一个实施例中，纵横波速度模型确定模块，用于按照如下方式，构建目标函数：
Figure BDA0002768009270000225
其中，d是观测数据，u是合成数据，xr是检波点的位置，*是时间卷积算子，xref表示要提取参考道的位置；
Figure BDA0002768009270000226
表示欧氏范数，m＝(α，β)是模型参数。
在具体实施本发明实施例提供的一种有限差分模拟弹性波全波形反演装置时，一个实施例中，纵横波速度模型确定模块，用于按照如下方式，获得目标函数对模型参数的梯度方程：
Figure BDA0002768009270000231
Figure BDA0002768009270000232
其中，α为纵波速度；β为横波速度；<,>t表示关于时间的内积，
Figure BDA0002768009270000233
表示伴随波场变量；ρ表示密度。
在具体实施本发明实施例提供的一种有限差分模拟弹性波全波形反演装置时，一个实施例中，纵横波速度模型确定模块，用于按照如下方式，求取更新方向：
Figure BDA0002768009270000234
其中，下标k表示迭代次数；上标T表示矩阵转置，并且
Figure BDA0002768009270000235
sk为共轭因子；m为模型参数向量，在此式中为纵横波速度组成的向量。
在具体实施本发明实施例提供的一种有限差分模拟弹性波全波形反演装置时，一个实施例中，纵横波速度模型确定模块，用于按照如下方式，进行迭代更新：
mk+1＝mk+tkΔmk+1
其中，tk表示选取的更新步长；下标k表示迭代次数；m为模型参数向量，在此式中为纵横波速度组成的向量。
综上，本发明实施例提供的一种有限差分模拟弹性波全波形反演方法和装置，采用速度-拟应力弹性波动方程，可以实现时间四阶精度的有限差分模拟，并且在该方程的差分模拟中，可以对纵横波采用不同的空间差分阶数进行波场模拟，避免常规弹性波方程数值模拟中采用相同的差分系数而造成的采样不足或者过采样的问题。同时，该模拟方法可以实现纵横波的解耦，对基于波场分离类的反演算法具有独特的优势。由于高阶时间差分允许使用较大的时间离散步长进行波场模拟，同时基于图形处理器(GPU)的共享内存的优化方式能进一步提高数据访问速度，从而加速弹性波动方程有限差分计算过程。
本领域内的技术人员应明白，本发明的实施例可提供为方法、系统、或计算机程序产品。因此，本发明可采用完全硬件实施例、完全软件实施例、或结合软件和硬件方面的实施例的形式。而且，本发明可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、CD-ROM、光学存储器等)上实施的计算机程序产品的形式。
本发明是参照根据本发明实施例的方法、设备(系统)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。
这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。
这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。
以上所述的具体实施例，对本发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施例而已，并不用于限定本发明的保护范围，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (17)

1.一种有限差分模拟弹性波全波形反演方法，其特征在于，包括：
根据纵横波速度模型，确定速度-拟应力弹性波动方程的纵横波差分参数；所述速度-拟应力弹性波动方程的纵横波差分参数，包括：纵横波有限差分离散格式，纵横波差分系数，纵横波稳定性条件，纵横波频散条件，初始条件，边界条件；
采用GPU共享内存优化策略，结合速度-拟应力弹性波动方程的纵横波差分参数，建立速度-拟应力弹性波动方程，进行时间四阶精度有限差分模拟，获得时间四阶精度有限差分模拟结果；
根据时间四阶精度有限差分模拟结果和速度-拟应力弹性波动方程，进行时间域弹性波全波形反演，确定纵横波速度模型。
2.如权利要求1所述的方法，其特征在于，按照如下方式，确定纵横波有限差分离散格式：
Figure FDA0002768009260000011
Figure FDA0002768009260000012
其中，
Figure FDA0002768009260000013
为纵波有限差分离散格式；
Figure FDA0002768009260000014
为横波有限差分离散格式；α为纵波速度；β为横波速度；u为待差分波场变量，r,φ∈{x,z}，并且r≠φ；hr表示沿着空间r方向的空间网格步长，Nr表示沿空间r方向的单边差分格式长度；u中的下标r表示沿着差分方向的网格点，u中的上标φ表示非差分轴向的网格点；c表示差分系数，其上标之一r表示差分沿着r方向，p表示差分系数由纵波速度计算，其另一上标s表示差分系数由横波速度计算。
3.如权利要求1所述的方法，其特征在于，按照如下方式，确定纵横波差分系数：
Figure FDA0002768009260000021
其中，c表示差分系数；在网格比γr,p＝αΔt/hr时，通过上式获得纵波差分系数，其中Δt表示二阶中心差分近似时间偏导数所使用的时间步长；相应地，将上式中γr,p相应替换为网格比γr,S＝βΔt/hr即可获得横波差分系数；α为纵波速度；β为横波速度；m为轴向差分系数次序索引，n为上式求和引入的变量。
4.如权利要求1所述的方法，其特征在于，按照如下方式，确定纵波稳定性条件：
Figure FDA0002768009260000022
其中，
Figure FDA0002768009260000023
r＝x,z；α为纵波速度；hx为x轴向网格间距，hz为z轴向网格间距；
按照如下方式，确定横波稳定性条件：
Figure FDA0002768009260000024
其中，
Figure FDA0002768009260000025
r＝x,z；β为横波速度；hx为x轴向网格间距，hz为z轴向网格间距。
5.如权利要求1所述的方法，其特征在于，按照如下方式，确定纵波频散条件：
Figure FDA0002768009260000026
其中，α为纵波速度；max{·}表示选取最大值，f0是子波的主频，N是某一差分精度下所允许的最小波长的最少采样点数；
按照如下方式，确定横波频散条件：
Figure FDA0002768009260000031
如下方式，β为横波速度；max{·}表示选取最大值，f0是子波的主频，N是某一差分精度下所允许的最小波长的最少采样点数。
6.如权利要求1所述的方法，其特征在于，按照如下方式，确定初始条件：
Figure FDA0002768009260000032
其中，v＝[vx vz]T；τ＝[τ111 τ11 τ12 τ22]T，τ111；τ11，τ12，τ22表示质点的应力分量；s(r0)表示震源函项，f(t)表示震源函数。
7.如权利要求1所述的方法，其特征在于，按照如下方式，确定边界条件：
Figure FDA0002768009260000033
Figure FDA0002768009260000034
其中，r∈{x,z}，δ是边界条件的层数，R是边界反射系数，其取值范围为10-6～10-8；α为纵波速度；β为横波速度；在上述边界条件中将vx分裂成x方向和z方向两部分，且将x方向的部分进一步分裂成纵波分量
Figure FDA0002768009260000035
和横波分量
Figure FDA0002768009260000036
即
Figure FDA0002768009260000037
8.如权利要求1所述的方法，其特征在于，所述GPU共享内存优化策略，包括：
从全局内存中将时间四阶精度有限差分模拟所需的数据拷贝至每个GPU块的共享内存中；
为每个GPU块沿着差分方向扩出差分长度的区域，使差分方向的每个点满足时间四阶精度差分格式；
根据时间四阶精度有限差分格式，从共享内存中读取时间四阶精度有限差分模拟数据进行计算。
9.如权利要求1所述的方法，其特征在于，按照如下方式，建立速度-拟应力弹性波动方程：
Figure FDA0002768009260000041
其中：x和z表示二维空间坐标系下的坐标，t表示时间坐标；λ和μ表示拉梅常数，ρ表示密度；λ+2μ＝ρα2，μ＝ρβ2，α和β分别表示纵横波速度；vx和vz分别表示质点的水平和垂直振动速度，τ111，τ11，τ12，τ22表示质点的应力分量；
Figure FDA0002768009260000042
表示时间一阶偏导数算子，
Figure FDA0002768009260000043
和
Figure FDA0002768009260000044
表示沿着x和z的一阶空间偏导数。
10.如权利要求1所述的方法，其特征在于，根据时间四阶精度有限差分模拟结果和速度-拟应力弹性波动方程，进行时间域弹性波全波形反演，确定纵横波速度模型，包括：
匹配观测数据和合成数据，构建目标函数；
根据速度-拟应力弹性波动方程和目标函数，采用拉格朗日伴随方法获得目标函数对模型参数的梯度方程；
根据时间四阶精度有限差分模拟结果进行正演和反传计算，获得震源和伴随波场变量；
利用目标函数对模型参数的梯度方程和震源、伴随波场变量获取模型参数梯度；
采用共轭梯度方法求取更新方向，对模型参数进行迭代更新，直至残差达到设定值，输出纵横波速度模型反演参数，确定纵横波速度模型。
11.如权利要求10所述的方法，其特征在于，按照如下方式，构建目标函数：
Figure FDA0002768009260000045
其中，d是观测数据，u是合成数据，xr是检波点的位置，*是时间卷积算子，xref表示要提取参考道的位置；
Figure FDA0002768009260000046
表示欧氏范数，m＝(α，β)是模型参数。
12.如权利要求10所述的方法，其特征在于，按照如下方式，获得目标函数对模型参数的梯度方程：
Figure FDA0002768009260000051
Figure FDA0002768009260000052
其中，α为纵波速度；β为横波速度；<,>t表示关于时间的内积，
Figure FDA0002768009260000053
表示伴随波场变量；ρ表示密度。
13.如权利要求10所述的方法，其特征在于，按照如下方式，求取更新方向：
Figure FDA0002768009260000054
其中，下标k表示迭代次数；上标T表示矩阵转置，并且
Figure FDA0002768009260000055
sk为共轭因子；m为模型参数向量，在此式中为纵横波速度组成的向量。
14.如权利要求10所述的方法，其特征在于，按照如下方式，进行迭代更新：
mk+1＝mk+tkΔmk+1
其中，tk表示选取的更新步长；下标k表示迭代次数；m为模型参数向量，在此式中为纵横波速度组成的向量。
15.一种有限差分模拟弹性波全波形反演装置，其特征在于，包括：
速度-拟应力弹性波动方程的纵横波差分参数确定模块，用于根据纵横波速度模型，确定速度-拟应力弹性波动方程的纵横波差分参数；所述速度-拟应力弹性波动方程的纵横波差分参数，包括：纵横波有限差分离散格式，纵横波差分系数，纵横波稳定性条件，纵横波频散条件，初始条件，边界条件；
时间四阶精度有限差分模拟结果确定模块，用于采用GPU共享内存优化策略，结合速度-拟应力弹性波动方程的纵横波差分参数，建立速度-拟应力弹性波动方程，进行时间四阶精度有限差分模拟，获得时间四阶精度有限差分模拟结果；
纵横波速度模型确定模块，用于根据时间四阶精度有限差分模拟结果和速度-拟应力弹性波动方程，进行时间域弹性波全波形反演，确定纵横波速度模型。
16.一种计算机设备，包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序，所述处理器执行所述计算机程序时实现权利要求1至14任一项所述有限差分模拟弹性波全波形反演方法。
17.一种计算机可读存储介质，所述计算机可读存储介质存储有执行实现权利要求1至14任一项所述有限差分模拟弹性波全波形反演方法的计算机程序。
CN202011239881.4A 2020-11-09 2020-11-09 有限差分模拟弹性波全波形反演方法和装置 Active CN114460640B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202011239881.4A CN114460640B (zh) 2020-11-09 2020-11-09 有限差分模拟弹性波全波形反演方法和装置

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202011239881.4A CN114460640B (zh) 2020-11-09 2020-11-09 有限差分模拟弹性波全波形反演方法和装置

Publications (2)

Publication Number Publication Date
CN114460640A true CN114460640A (zh) 2022-05-10
CN114460640B CN114460640B (zh) 2024-06-25

Family

ID=81404558

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202011239881.4A Active CN114460640B (zh) 2020-11-09 2020-11-09 有限差分模拟弹性波全波形反演方法和装置

Country Status (1)

Country Link
CN (1) CN114460640B (zh)

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115421190A (zh) * 2022-07-29 2022-12-02 南方科技大学 一种地下结构弹性波正演声波反演成像方法及设备
CN115453621A (zh) * 2022-09-14 2022-12-09 成都理工大学 基于一阶速度-应力方程的纵横波解耦分离假象去除方法

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103091711A (zh) * 2013-01-24 2013-05-08 中国石油天然气集团公司 全波形反演方法及装置
CN104965223A (zh) * 2015-05-29 2015-10-07 中国石油天然气股份有限公司 粘声波全波形反演方法及装置
CN105158797A (zh) * 2015-10-16 2015-12-16 成都理工大学 一种基于实际地震资料的交错网格波动方程正演的方法
CN105319581A (zh) * 2014-07-31 2016-02-10 中国石油化工股份有限公司 一种高效的时间域全波形反演方法
CN105467444A (zh) * 2015-12-10 2016-04-06 中国石油天然气集团公司 一种弹性波全波形反演方法及装置

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103091711A (zh) * 2013-01-24 2013-05-08 中国石油天然气集团公司 全波形反演方法及装置
CN105319581A (zh) * 2014-07-31 2016-02-10 中国石油化工股份有限公司 一种高效的时间域全波形反演方法
CN104965223A (zh) * 2015-05-29 2015-10-07 中国石油天然气股份有限公司 粘声波全波形反演方法及装置
CN105158797A (zh) * 2015-10-16 2015-12-16 成都理工大学 一种基于实际地震资料的交错网格波动方程正演的方法
CN105467444A (zh) * 2015-12-10 2016-04-06 中国石油天然气集团公司 一种弹性波全波形反演方法及装置

Non-Patent Citations (6)

* Cited by examiner, † Cited by third party
Title
PETROV, PV 等: "3D finite-difference modeling of elastic wave propagation in the Laplace-Fourier domain", GEOPHYSICS, vol. 77, no. 4, 5 July 2012 (2012-07-05), pages 137 - 155, XP001576794, DOI: 10.1190/geo2011-0238.1 *
任志明: "声波和弹性波波动方程有限差分正反演方法研究", pages 011 - 80 *
岳晓鹏: "全波形反演方法技术研究", pages 011 - 16 *
张光超: "基于GPU并行计算的时间域弹性波全波形反演方法研究", pages 011 - 496 *
王守进: "三维弹性波高阶有限差分模拟及GPU加速" *
陈汉明: "波动方程数值模拟与粘滞波形反演方法研究", pages 011 - 58 *

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115421190A (zh) * 2022-07-29 2022-12-02 南方科技大学 一种地下结构弹性波正演声波反演成像方法及设备
CN115453621A (zh) * 2022-09-14 2022-12-09 成都理工大学 基于一阶速度-应力方程的纵横波解耦分离假象去除方法
CN115453621B (zh) * 2022-09-14 2024-03-22 成都理工大学 基于一阶速度-应力方程的纵横波解耦分离假象去除方法

Also Published As

Publication number Publication date
CN114460640B (zh) 2024-06-25

Similar Documents

Publication Publication Date Title
Trinh et al. Efficient time-domain 3D elastic and viscoelastic full-waveform inversion using a spectral-element method on flexible Cartesian-based mesh
Etgen et al. Computational methods for large-scale 3D acoustic finite-difference modeling: A tutorial
CN103091711B (zh) 基于时间域一阶弹性波动方程的全波形反演方法及装置
US8296069B2 (en) Pseudo-analytical method for the solution of wave equations
KR101948509B1 (ko) 지구 물리학적 데이터의 반복 반전의 아티팩트 감소
CN104122585A (zh) 基于弹性波场矢量分解与低秩分解的地震正演模拟方法
CN105137486A (zh) 各向异性介质中弹性波逆时偏移成像方法及其装置
Fang et al. Elastic full-waveform inversion based on GPU accelerated temporal fourth-order finite-difference approximation
CN113740901B (zh) 基于复杂起伏地表的陆上地震数据全波形反演方法及装置
CN110954945B (zh) 一种基于动态随机震源编码的全波形反演方法
Qu et al. Fluid-solid coupled full-waveform inversion in the curvilinear coordinates for ocean-bottom cable data
MX2011003850A (es) Estimado de señal de dominio de imagen a interferencia.
CN104965222B (zh) 三维纵波阻抗全波形反演方法及装置
CN114460640B (zh) 有限差分模拟弹性波全波形反演方法和装置
Terrana et al. A spectral hybridizable discontinuous Galerkin method for elastic–acoustic wave propagation
Habashy et al. Source-receiver compression scheme for full-waveform seismic inversion
Deschizeaux et al. Imaging earth’s subsurface using CUDA
Wu et al. A new full waveform inversion method based on shifted correlation of the envelope and its implementation based on OPENCL
CN114139335A (zh) 基于单松弛时间格子玻尔兹曼模型的粘滞声波模拟方法
Song et al. An efficient wavefield inversion for transversely isotropic media with a vertical axis of symmetry
US11199641B2 (en) Seismic modeling
CN111665556A (zh) 地层声波传播速度模型构建方法
Fang et al. 3D elastic multisource full waveform inversion on distributed GPU clusters
CN115373022A (zh) 一种基于振幅相位校正的弹性波场Helmholtz分解方法
Wang et al. Rock Fracture Monitoring Based on High‐Precision Microseismic Event Location Using 3D Multiscale Waveform Inversion

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