CN112444849A - 一种基于交错网格低秩有限差分的弹性逆时偏移成像方法 - Google Patents

一种基于交错网格低秩有限差分的弹性逆时偏移成像方法 Download PDF

Info

Publication number
CN112444849A
CN112444849A CN201910807274.4A CN201910807274A CN112444849A CN 112444849 A CN112444849 A CN 112444849A CN 201910807274 A CN201910807274 A CN 201910807274A CN 112444849 A CN112444849 A CN 112444849A
Authority
CN
China
Prior art keywords
wave
representing
elastic
longitudinal
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.)
Pending
Application number
CN201910807274.4A
Other languages
English (en)
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 Petroleum and Chemical Corp
Sinopec Geophysical Research Institute
Original Assignee
China Petroleum and Chemical Corp
Sinopec Geophysical Research Institute
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 Petroleum and Chemical Corp, Sinopec Geophysical Research Institute filed Critical China Petroleum and Chemical Corp
Priority to CN201910807274.4A priority Critical patent/CN112444849A/zh
Publication of CN112444849A publication Critical patent/CN112444849A/zh
Pending legal-status Critical Current

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/282—Application of seismic models, synthetic seismograms
    • 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/301—Analysis for determining seismic cross-sections or geostructures
    • 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/303—Analysis for determining velocity profiles or travel times

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)
  • Geophysics And Detection Of Objects (AREA)

Abstract

本发明提供了一种基于交错网格低秩有限差分的弹性逆时偏移成像方法和计算机存储介质以及计算机系统。该方法借助平面波场分解、低秩分解方法及傅里叶变换等方法，在时间‑空间域中构建了时空兼顾的高精度弹性波场延拓算子，基于此算子可实现弹性波场延拓并得到纵横波矢量波场，延拓结果可以直接应用于弹性逆时偏移成像，从而使得弹性逆时偏移无需额外进行波场分离和极性校正，极大地提升弹性逆时偏移成像的质量。

Description

一种基于交错网格低秩有限差分的弹性逆时偏移成像方法
技术领域
本发明属于地震资料叠前深度偏移成像技术领域，尤其涉及一种基于交错网格低秩有限差分的弹性逆时偏移成像方法和计算机存储介质以及计算机系统。
背景技术
弹性逆时偏移隶属于逆时偏移的范畴，是一种面向多分量地震资料的处理技术。它不受倾角限制及横向变速的影响，能够充分保持地下波场的矢量特性。目前，作为前沿研究课题的弹性逆时偏移理论与方法体系有待深入研究。为满足高精度地震勘探的要求，需要提高成像精度，高精度的波场延拓方法对提高弹性逆时偏移的成像精度具有重要的意义。常规交错网格有限差分方法和伪谱法是两种常用的弹性波场延拓方法。其中，就常规交错网格有限差分方法而言，受限于Courant-Friedrch-Lewy稳定性条件，时间和空间方向需采用合适的延拓步长才能够保证波场延拓的精度和稳定性；虽然伪谱法可视作有限差分算子的空间极值，但是它难以适应不均匀介质波场延拓，并且不满足大时间步长的波场延拓要求。
发明内容
针对上述问题，本发明提出了一种基于交错网格低秩有限差分的弹性逆时偏移方法。
根据本发明的第一个方面，本发明的一种基于交错网格低秩有限差分的弹性逆时偏移成像方法，包括以下步骤：
S100，利用低秩分解法构建适应于一阶弹性波方程的交错网格有限差分系数；
S200，根据所述交错网格有限差分系数确定交错网格低秩有限差分算子，其中，所述交错网格低秩有限差分算子包括沿坐标轴负方向交错并与P波相关及S波相关的空间偏导数以及沿坐标轴正方向交错并与P波相关及与S波相关的空间偏导数；
S300，将所述交错网格低秩有限差分算子代入递归时间波场延拓方程，获得基于交错网格低秩有限差分离散化的弹性波动方程；
S400，根据速度模型和所述弹性波动方程进行正向弹性波场延拓，获得正向延拓的纵横波矢量波场；
S500，根据地震资料和所述弹性波动方程进行反向弹性波场延拓，获得反向延拓的纵横波矢量波场；
S600，根据基于正向延拓和反向延拓的纵横波矢量波场完成弹性逆时偏移成像。
根据本发明的实施例，上述步骤S100包括以下步骤：
S110，利用傅里叶变换将一阶弹性波方程从时间域转换到波数域，以获得相应的波数域方程；
S120，根据表征纵波和横波的质点振动的速度矢量、整个波场质点振动的速度矢量同纵横波质点振动的速度矢量之间的关系以及弹性递归时间延拓算子，对所述波数域方程进行整理，以获得一阶弹性波递归时间延拓方程；
S130，根据速度和应力在交错网格的定义关系，对一阶弹性波递归时间延拓方程进行整理，以获得交错网格一阶弹性波递归时间积分方程；
S140，利用经由傅里叶变换得到的适用于一阶弹性波递归时间积分方程的交错网格弹性递归时间积分算子，并通过低秩分解法对交错网格弹性递归时间积分算子进行分解近似；
S150，基于分解近似后的交错网格弹性递归时间积分算子确定适用于交错网格一阶弹性递归时间积分方程的交错网格有限差分系数。
根据本发明的实施例，上述步骤S110中，
所述各向同性介质中的一阶弹性波方程为：
Figure BDA0002182862640000021
其中，λ和μ分别表示拉梅系数，ρ为密度，τjl表示应力张量，vl、vj及vm表示速度分量，δjl为克罗内克算符，t表示时间变量，x表示空间变量，下标j，l，m表示笛卡尔坐标下的坐标分量；
所述相应的波数域方程为：
Figure BDA0002182862640000031
Figure BDA0002182862640000032
及
Figure BDA0002182862640000033
表示波数域中速度分量，kj、km及kl表示波数变量，
Figure BDA0002182862640000034
表示波数域中的应力张量。
根据本发明的实施例，上述步骤S120中，整个波场质点振动的速度矢量同纵横波质点振动的速度矢量之间的关系为：
Figure BDA0002182862640000035
vj为整个波场质点振动的速度矢量，
Figure BDA0002182862640000036
为纵波质点振动的速度矢量，
Figure BDA0002182862640000037
为横波质点振动的速度矢量。
根据本发明的实施例，上述基于低秩分解方法的交错网格有限差分系数为：
Figure BDA0002182862640000038
式中，Gpj(x，l)及Gsj(x，l)分别表示基于低秩分解法得到的纵波或横波相关交错网格有限差分系数，wpj(x，km)及wsj(x，km)分别表示纵波或横波空间相关矩阵，amn表示中间矩阵，Cpj(xn，l)及Csj(xn，l)分别表示纵波或横波相关三角函数系数矩阵。
根据本发明的实施例，上述数值离散化弹性波动方程具体为：
Figure BDA0002182862640000039
Figure BDA00021828626400000310
Figure BDA0002182862640000041
上式中，Δt表示时间步长，vx表示x方向整个波场质点振动速度矢量，
Figure BDA0002182862640000042
表示x方向纵波相关质点振动速度矢量，
Figure BDA0002182862640000043
表示x方向横波相关质点振动速度矢量，vz表示z方向整个波场质点振动速度矢量，
Figure BDA0002182862640000044
表示z方向纵波相关质点振动速度矢量，
Figure BDA0002182862640000045
表示z方向横波相关质点振动速度矢量，
Figure BDA0002182862640000046
表示t+Δt/2时刻处x方向纵波相关质点振动速度矢量，
Figure BDA0002182862640000047
表示t-Δt/2时刻处x方向纵波相关质点振动速度矢量，
Figure BDA0002182862640000048
表示t+Δt/2时刻处x方向横波相关质点振动速度矢量，
Figure BDA0002182862640000049
表示t-Δt/2时刻处x方向横波相关质点振动速度矢量，
Figure BDA00021828626400000410
表示t+Δt/2时刻处z方向纵波相关质点振动速度矢量，
Figure BDA00021828626400000411
表示t-Δt/2时刻处z方向纵波相关质点振动速度矢量，
Figure BDA00021828626400000412
表示t+Δt/2时刻处z方向横波相关质点振动速度矢量，
Figure BDA00021828626400000413
表示t-Δt/2时刻处z方向横波相关质点振动速度矢量，ρ1表示与参考网格点x方向相差一个空间步长处的密度值，ρ2表示表示与参考网格点z方向相差一个空间步长处的密度值，τp表示纵波应力，
Figure BDA00021828626400000414
表示横波相关x方向正应力，
Figure BDA00021828626400000415
表示横波相关切应力，
Figure BDA00021828626400000416
表示横波相关z方向正应力，
Figure BDA00021828626400000417
表示t+Δt时刻处横波相关z方向正应力，
Figure BDA00021828626400000418
表示t时刻处横波相关z方向正应力，
Figure BDA00021828626400000419
表示t+Δt时刻处横波相关x方向正应力，
Figure BDA00021828626400000420
表示t时刻处横波相关x方向正应力，
Figure BDA00021828626400000421
表示t+Δt时刻横波相关切应力，
Figure BDA00021828626400000422
表示t时刻处横波相关切应力，τp(t+Δt)表示t+Δt时刻处纵波应力，τp(t)表示t时刻处纵波应力，μ0及λ0表示参考网格点处的拉梅系数，μ1，2表示与参考网格点x方向和z方向均相差一个空间步长的拉梅系数，
Figure BDA00021828626400000423
及
Figure BDA00021828626400000424
分别表示与纵波相关且沿x坐标轴正负方向交错的空间偏导数，
Figure BDA00021828626400000425
及
Figure BDA00021828626400000426
分别表示与纵波相关且沿z坐标轴正负方向交错的空间偏导数，
Figure BDA00021828626400000427
及
Figure BDA00021828626400000428
分别表示与横波相关且沿x坐标轴正负方向交错的空间偏导数，
Figure BDA00021828626400000429
及
Figure BDA0002182862640000051
分别表示与横波相关且沿z坐标轴正负方向交错的空间偏导数。
根据本发明的实施例，上述步骤S600具体为：
针对源检波场，采用矢量内积归一化方法，根据正向延拓和反向延拓的纵横波矢量波场获得PS及PP成像结果。
根据本发明的实施例，上述步骤S600具体为：针对源检波场，采用矢量内积归一化方法，按照下式根据正向延拓和反向延拓的纵横波矢量波场获得PS及PP成像结果：
Figure BDA0002182862640000052
式中，Ipp表示PP共成像点，Ips表示PS共成像点，Sp(x，t)为P波震源波场，Rp(x，t)为P波检波波场，Rs(x，t)为S波检波波场。
根据本发明的实施例，本发明还提供一种计算机存储介质，其中存储有用于实现上述方法的计算机程序。
根据本发明的实施例，本发明还提供一种计算机系统，包括存储器和处理器，所述处理器用于执行所述存储器中存储的计算机程序，所述计算机程序用于实现上述方法。
与现有技术相比，本发明提供的基于交错网格低秩有限差分的弹性逆时偏移成像方法和系统等具有如下优点或有益效果：
该方法借鉴交错网格的思想及平面波分解原理，构建能够兼顾时间和空间高精度的交错网格空间-波数域波场延拓算子。由于常规的空间-波数域波场延拓算子计算量巨大，借助低秩分解算法及傅里叶变换特性可将其作分解近似，得到基于交错网格低秩有限差分的弹性波场延拓算子。该算子兼具时间和空间上的高精度，同时与空间-波数域波场延拓算子相比，其计算复杂度明显较低。此外，在利用基于交错网格地质有限差分弹性波场延拓算子进行弹性波场延拓时，可以直接得到解耦的纵横波，该解耦的纵横波并可直接应用于弹性逆时偏移成像，使得在弹性逆时偏移成像过程中无需额外进行波场分离(分解)和极性校正，并且获得的偏移成像剖面质量更高。本发明的方法对于推动高精度成像具有重要意义。
本发明的其它特征和优点将在随后的说明书中阐述，并且，部分地从说明书中变得显而易见，或者通过实施本发明而了解。本发明的目的和其他优点可通过在说明书、权利要求书以及附图中所特别指出的结构来实现和获得。
附图说明
从下面描述的实施例并参考附图，本发明的其它优点和细节将变得显而易见。
以下是示意图并示出：
图1是本发明实施例一的基于交错网格低秩有限差分的弹性逆时偏移成像方法的步骤流程图；
图2是本发明实施例二的单层介质模型及模型参数；
图3是本发明实施例二的利用单层介质模型进行弹性波场延拓得到的1.5s处的弹性波场快照：(a)为利用4阶常规交错网格有限差分法得到的P波水平分量波场快照，(b)为利用4阶交错网格低秩有限差分得到的P波水平分量波场快照，(c)为利用4阶常规交错网格有限差分得到的S波水平分量波场快照，(d)为利用4阶交错网格低秩有限差分得到的S波水平分量波场快照；
图4是从图3所示的波场快照中抽取的单道数据：(a)为利用常规交错网格有限差分法得到的与P波相关的单道数据，(b)为利用交错网格低秩有限差分法得到的与P波相关的单道数据，(c)为利用常规交错网格有限差分法得到的与S波相关的单道数据，(d)为利用交错网格低秩有限差分法得到的与S波相关的单道数据；
图5是本发明实施例二的Marmousi2模型；
图6是利用图5的Marmousi2模型得到的波场快照：(a)为P波水平分量波场快照，(b)为P波垂直分量波场快照，(c)为S波水平分量波场快照，(d)为S波垂直分量波场快照；
图7是本发明实施例二的共炮点道集：(a)为水平分量，(b)为垂直分量；
图8是本发明实施例二的Marmousi2多炮叠加成像剖面：(a)为PP成像剖面，(b)为PS成像剖面。
具体实施方式
针对常规交错网格有限差分方法难以同时满足空间和时间方向高精度弹性逆时偏移的要求，本发明介绍了一种基于交错网格低秩有限差分的弹性逆时偏移方法，其借助平面波场分解、低秩分解方法及傅里叶变换等方法，在时间-空间域中构建了时空兼顾的高精度弹性波场延拓算子，基于此算子可实现弹性波场延拓并得到纵横波矢量波场，延拓结果可直接应用于弹性逆时偏移过程，从而使得弹性逆时偏移无需额外进行波场分离(分解)及极性校正过程，可极大提升弹性逆时偏移的成像质量。
实施例一
下面详细地介绍本发明的基于交错网格低秩有限差分的弹性逆时偏移方法的原理。
在各向同性介质中，一阶弹性波方程的表达式为：
Figure BDA0002182862640000071
其中，λ和μ分别表示拉梅系数，ρ为密度，τjl表示应力张量，vjl表示速度分量，δjl为克罗内克算符，t表示时间变量，x表示空间变量，下标j，l，m表示笛卡尔坐标下的坐标分量。波数域中，可将式(1)表示为：
Figure BDA0002182862640000072
其中，i为复数单位，k是波数变量，
Figure BDA0002182862640000073
和
Figure BDA0002182862640000074
分别是波数域中质点振动的速度分量与应力张量。可使用伪谱法对式(2)进行数值计算，伪谱法是一种比较间接的求解空间导数的方法，其空间精度高于有限差分方法，然而在时间方向上，伪谱算子的精度并不高。将耦合速度、波数及时间步长的弹性递归时间延拓算子
Figure BDA0002182862640000075
其中，cp，s表示纵波或横波速度)替换式(2)中的波数变量可提升波场延拓精度。因式(2)中质点振动的速度矢量以及应力张量既与纵波相关又与横波也相关，难以直接利用上述算子。为此，定义纵波应力张量
Figure BDA0002182862640000076
和横波的应力张量
Figure BDA0002182862640000077
并引入一组表征纵波和横波的质点振动的速度变量
Figure BDA0002182862640000078
和
Figure BDA0002182862640000079
在线性假设下，整个波场质点振动的速度矢量同纵横波质点振动的速度矢量之间的关系可以表示为：
Figure BDA0002182862640000081
利用表征纵波和横波的质点振动的速度变量、整个波场质点振动的速度矢量同纵横波质点振动的速度矢量之间的关系及弹性递归时间延拓算子对式(2)进行整理，可得一阶弹性波递归时间延拓方程：
Figure BDA0002182862640000082
Figure BDA0002182862640000083
根据各物理参量(速度和应力)在交错网格定义的关系，对式(4)和式(5)作进一步整理，可得交错网格一阶弹性递归时间积分方程。其中，τx、τxx和τzz位于网格参考点位置处，切应力τxz同上述应力在x轴方向和z轴方向分别相差半个网格点；垂直方向上质点振动的速度分量
Figure BDA0002182862640000084
同参考点在x轴方向相差半个网格点；水平方向上质点振动速度分量
Figure BDA0002182862640000085
同参考点在z轴方向相差半个网格点。
在波数域中，可利用傅里叶变换平移性质对上述交错网格中各物理参量的对应关系进行表述。为表述简洁，定义了如下的辅助变量：
Figure BDA0002182862640000091
利用局部傅里叶变换可对式(6)进行数值求解。此方法需要在混合域中进行，每个时间步长内均需数次正反傅里叶变换，此举将极大增加计算量。令Nx表示总的离散网格点数，其计算复杂度为
Figure BDA0002182862640000092
从数学意义上讲，上述弹性递归时间积分算子是光滑、震荡的，一般情况下非满秩，可做低秩分解近似。以与P波相关且沿坐标轴正方向交错的空间偏导数，
Figure BDA0002182862640000093
为例，对交错网格弹性递归时间积分算子的分解近似过程做详细说可明。
Figure BDA0002182862640000094
的具体表达式为
Figure BDA0002182862640000095
令wpj(x，k)＝kj sinc(cp(x)|k|Δt/2)，其为Nx×Nx大小的矩阵；可将其分解为
Figure BDA0002182862640000096
式(5)中m和n均代表矩阵wpj(x，k)的秩，wpj(x，km)为空间相关矩阵，wpj(xn，k)为波数相关矩阵，amn为中间矩阵，矩阵wpj(x，km)只有波数相关，
j为指标集，j＝1，3，可做进一步分解：
Figure BDA0002182862640000097
其中，
Figure BDA0002182862640000098
是L×Nx的子矩阵，L是交错网格差分算子长度，
Figure BDA0002182862640000099
是三角函数系数矩阵，
Figure BDA00021828626400000910
表示伪逆算子，定义Nx×L的矩阵：
Figure BDA00021828626400000911
经过低秩分解后，式(7)可表示为
Figure BDA0002182862640000101
根据欧拉式及傅里叶变换的性质，可做进一步算术整理，得
Figure BDA0002182862640000102
沿坐标轴负方向交错并与P波相关及S波相关的空间偏导数
Figure BDA0002182862640000103
以及沿坐标轴正方向交错并与P波相关及与S波相关的空间偏导数
Figure BDA0002182862640000104
的求取过程同上述过程相类似，在此不再赘述。
将上述空间偏导数(即交错网格低秩有限差分算子)带入式(4)及式(5)中，可得基于交错网格低秩有限差分离散后的一阶弹性波动方程：
Figure BDA0002182862640000105
Figure BDA0002182862640000106
Figure BDA0002182862640000111
其中，μ0≡μ(x，z)，μ1，2≡μ(x+Δx，z+Δz)，ρ1≡ρ(x+Δx，z)，ρ2≡ρ(x，z+Δz)及λ0≡λ(x，z)。
在现有的弹性逆时偏移成像过程中，主要是利用常规有限差分方法对弹性波场进行正向和反向波场延拓。这种常规有限差分法的波场延拓精度不高，不利于复杂模型波场延拓，同时，为避免偏移成像过程中不同波型之间的串扰现象，需在应用成像条件前进行波场分离，且针对PS成像剖面需要进行极性校正。
与上述方法不同，本发明的方法利用基于交错网格低秩有限差分算子进行正向和反向弹性波场延拓，可以直接获得高精度的正向延拓和反向延拓的纵横波矢量波场(P/S波矢量波场)，这种矢量波场可以直接应用于弹性逆时偏移成像过程中，所得成像结果无需进行极性校正。
因此，在具体实施时，可以通过式(11)分别构建沿坐标轴负方向交错并与P波相关及S波相关的空间偏导数以及沿坐标轴正方向交错并与P波相关及与S波相关的空间偏导数，利用式(12)、式(13)、式(14)以及速度模型、地震资料分别构建正向延拓和反向延拓的高精度纵横波矢量波场。然后针对源检波场，采用矢量内积归一化方法，基于正向延拓和反向延拓的高精度纵横波矢量波场完成弹性逆时偏移成像过程。该方法既可以保证纵横波矢量波场的精度，又可以避免了成像过程中纵横波矢量波场串扰现象以及转换波剖面极性反转现象，可以有效提高成像精度。
具体的实现步骤为：
步骤一：借助低秩分解方法构建交错网格有限差分系数
Figure BDA0002182862640000121
步骤二：采用基于交错网格低秩有限差分算子得到的离散化波动方程，分别实现正向及反向解耦的弹性波场延拓
Figure BDA0002182862640000122
Figure BDA0002182862640000123
Figure BDA0002182862640000124
步骤三：采用矢量内积归一化方法得到PP及PS成像剖面
Figure BDA0002182862640000125
式中，Ipp表示PP共成像点，Ips表示PS共成像点，Sp(x，t)为P波震源波场，Rp(x，t)为P波检波波场，Rs(x，t)为S波检波波场。
实施例二
下面结合实际应用来说明，本发明提出的方法的有效性和进步性。
图2为选取模型大小为600×600的单层介质模型，空间采样间隔为15m，时间采样间隔为1.5ms，界面位于3000m处。上层密度、纵波速度及横波速度分别为2.0g/cm3，1700m/s及1000m/s，下层密度、纵波速度及横波速度分别为2.2g/cm3、2550m/s、1500m/s。以主频为35Hz的雷克子波作为纵波震源，并分别采用常规交错网格有限差分方法和交错网格低秩有限差分方法进行弹性波场数值模拟，所得波场快照如图3所示。
图3是利用单层介质模型进行弹性波场延拓得到的1.5s处的弹性波场快照：(a)为利用4阶常规交错网格有限差分法得到的P波水平分量波场快照，(b)为利用4阶交错网格低秩有限差分得到的P波水平分量波场快照，(c)为利用4阶常规交错网格有限差分得到的S波水平分量波场快照，(d)为利用4阶交错网格低秩有限差分得到的S波水平分量波场快照。
图4从图3所示的波场快照中抽取的单道数据：(a)为利用常规交错网格有限差分法得到的与P波相关的单道数据，(b)为利用交错网格低秩有限差分法得到的与P波相关的单道数据，(c)为利用常规交错网格有限差分法得到的与S波相关的单道数据，(d)为利用交错网格低秩有限差分法得到的与S波相关的单道数据。
从图3中抽取对应的单道地震记录与常规交错网格有限差分方法得到的波场快照及单道地震记录相比，交错网格低秩有限差分方法得到波场快照几乎不受频散影响，子波波形特征保持良好，由此可说明当空间采样间隔和时间采样间隔较大时交错网格有限差分法得到的波场快照受频散影响程度较小，满足弹性波场数值模拟的需求。
图5表示的是Marmousi2模型，模型大小为1700×350，水平方向和垂直方向的空间采样间隔为10m。利用图中白色选项方框中的模型(白色方框四个顶角所在的具体位置为(6.0，0)km、(13.0，0)km、(6.0，3.5)km、(13.0，3.5)km。模型网格大小为700×350，利用交错网格低秩有限差分方法进行弹性波场数值模拟得到图6及图7结果。图6是利用Marmousi2模型得到的波场快照：(a)为P波水平分量波场快照，(b)为P波垂直分量波场快照，(c)为S波水平分量波场快照，(d)为S波垂直分量波场快照。图7显示的是共炮点道集：(a)水平分量，(b)垂直分量。
由图可知，P波水平分量波场快照、P波垂直分量波场快照、S波水平分量波场快照、S波垂直分量波场快照及共炮点道集均未受数值频散影响，在波场快照中P波及S波形态清晰可见。此例可以充分说明本发明的交错网格低秩有限差分方法能够适用复杂介质模型的弹性波场数值模拟，同时也说明S波对地下构造的刻画更为细致，是描绘地下构造的重要的工具。
图8表示的是利用交错网格低秩有限差分方法进行正向及反向弹性波场延拓并利用矢量内积归一化方法得到的PS及PP成像结果，在图中可见薄层、三大断层、不整合面及背斜等构造均清晰展现。
上述模型测试可表明，本发明的基于交错网格低秩有限差分的弹性逆时偏移成像方法能够满足高精度的弹性波偏移成像要求。该方法利用平面波理论及低秩分解方法构建了同时考虑时间和空间精度的交错网格低秩有限差分算子，能够满足高精度弹性波场延拓的需求。同时，该方法能够充分利用多分量地震资料，在延拓过程中可得高精度纵横波矢量，充分保持了地下波场的矢量特性，所得的纵横波波场矢量可以直接应用于成像过程，避免了常规弹性逆时偏移需额外进行波场分离及极性校正等工作，成像结果清晰可靠，十分利于后续处理解释工作的开展。
实施例三
在本实施例中，本发明还提供一种计算机存储介质，其中存储有用于实现上述方法的计算机程序。
实施例四
在本实施例中，本发明还提供一种计算机设备，包括存储器和处理器，所述处理器用于执行所述存储器中存储的计算机程序，其中所述计算机程序用于实现上述方法。
应该理解的是，本发明所公开的实施例不限于这里所公开的特定处理步骤或材料，而应当延伸到相关领域的普通技术人员所理解的这些特征的等同替代。还应当理解的是，在此使用的术语仅用于描述特定实施例的目的，而并不意味着限制。
说明书中提到的“实施例”意指结合实施例描述的特定特征、或特性包括在本发明的至少一个实施例中。因此，说明书通篇各个地方出现的短语“实施例”并不一定均指同一个实施例。
本领域的技术人员应该明白，结合本文中所公开的实施例描述的各示例的单元及算法步骤，能够以电子硬件、计算机软件或者二者的结合来实现，为了清楚地说明硬件和软件的可互换性，在上述说明中已经按照功能一般性地描述了各示例的组成及步骤。这些功能究竟以硬件还是软件方式来执行，取决于技术方案的特定应用和设计约束条件。本领域的技术人员可以对每个特定的应用来使用不同方法来实现所描述的功能，但是这种实现不应认为超出本发明的范围。
结合本文中所公开的实施例描述的方法或算法的步骤可以直接用硬件、处理器执行的软件模块，或者二者的结合来实施。软件模块可以置于随机存储介质(RAM)、内存、只读存储介质(ROM)、电可编程ROM、电可擦除可编程ROM、寄存器、硬盘、可移动磁盘、CD-ROM或技术领域内所公知的任意其它形式的存储介质中。
虽然本发明所公开的实施方式如上，但所述的内容只是为了便于理解本发明而采用的实施方式，并非用以限定本发明。任何本发明所属技术领域内的技术人员，在不脱离本发明所公开的精神和范围的前提下，可以在实施的形式上及细节上作任何的修改与变化，但本发明的保护范围，仍须以所附的权利要求书所界定的范围为准。

Claims (10)

1.一种基于交错网格低秩有限差分的弹性逆时偏移成像方法，其特征在于，包括以下步骤：
S100，利用低秩分解法构建适应于一阶弹性波方程的交错网格有限差分系数；
S200，根据所述交错网格有限差分系数确定交错网格低秩有限差分算子，其中，所述交错网格低秩有限差分算子包括沿坐标轴负方向交错并与P波相关及S波相关的空间偏导数以及沿坐标轴正方向交错并与P波相关及与S波相关的空间偏导数；
S300，将所述交错网格低秩有限差分算子代入递归时间波场延拓方程，获得基于交错网格低秩有限差分离散化的弹性波动方程；
S400，根据速度模型和所述弹性波动方程进行正向弹性波场延拓，获得正向延拓的纵横波矢量波场；
S500，根据地震资料和所述弹性波动方程进行反向弹性波场延拓，获得反向延拓的纵横波矢量波场；
S600，根据基于正向延拓和反向延拓的纵横波矢量波场完成弹性逆时偏移成像。
2.根据权利要求1所述的方法，其特征在于，所述步骤S100包括以下步骤：
S110，利用傅里叶变换将一阶弹性波方程从时间域转换到波数域，以获得相应的波数域方程；
S120，根据表征纵波和横波的质点振动的速度矢量、整个波场质点振动的速度矢量同纵横波质点振动的速度矢量之间的关系以及弹性递归时间延拓算子，对所述波数域方程进行整理，以获得一阶弹性波递归时间延拓方程；
S130，根据速度和应力在交错网格的定义关系，对一阶弹性波递归时间延拓方程进行整理，以获得交错网格一阶弹性波递归时间积分方程；
S140，利用经由傅里叶变换得到的适用于一阶弹性波递归时间积分方程的交错网格弹性递归时间积分算子，并通过低秩分解法对交错网格弹性递归时间积分算子进行分解近似；
S150，基于分解近似后的交错网格弹性递归时间积分算子确定适用于交错网格一阶弹性递归时间积分方程的交错网格有限差分系数。
3.根据权利要求2所述的方法，其特征在于，所述步骤S110中，
所述各向同性介质中的一阶弹性波方程为：
Figure FDA0002182862630000021
其中，λ和μ分别表示拉梅系数，ρ为密度，τjl表示应力张量，vl、vj及vm表示速度分量，δjl为克罗内克算符，t表示时间变量，x表示空间变量，下标j，l，m表示笛卡尔坐标下的坐标分量；
所述相应的波数域方程为：
Figure FDA0002182862630000022
Figure FDA0002182862630000023
及
Figure FDA0002182862630000024
表示波数域中速度分量，kj、km及kl表示波数变量，
Figure FDA0002182862630000025
表示波数域中的应力张量。
4.根据权利要求3所述的方法，其特征在于，所述步骤S120中，整个波场质点振动的速度矢量同纵横波质点振动的速度矢量之间的关系为：
Figure FDA0002182862630000026
vj为整个波场质点振动的速度矢量，
Figure FDA0002182862630000027
为纵波质点振动的速度矢量，
Figure FDA0002182862630000028
为横波质点振动的速度矢量。
5.根据权利要求4所述的方法，其特征在于，所述基于低秩分解方法的交错网格有限差分系数为：
Figure FDA0002182862630000029
式中，Gpj(x，l)及Gsj(x，l)分别表示基于低秩分解法得到的纵波或横波相关交错网格有限差分系数，wpj(x，km)及wsj(x，km)分别表示纵波或横波空间相关矩阵，amn表示中间矩阵，Cpj(xn，l)及Csj(xn，l)分别表示纵波或横波相关三角函数系数矩阵。
6.根据权利要求3所述的方法，其特征在于，所述数值离散化弹性波动方程具体为：
Figure FDA0002182862630000031
Figure FDA0002182862630000032
Figure FDA0002182862630000033
上式中，Δt表示时间步长，vx表示x方向整个波场质点振动速度矢量，
Figure FDA0002182862630000034
表示x方向纵波相关质点振动速度矢量，
Figure FDA0002182862630000035
表示x方向横波相关质点振动速度矢量，vz表示z方向整个波场质点振动速度矢量，
Figure FDA0002182862630000036
表示z方向纵波相关质点振动速度矢量，
Figure FDA0002182862630000037
表示z方向横波相关质点振动速度矢量，
Figure FDA0002182862630000038
表示t+Δt/2时刻处x方向纵波相关质点振动速度矢量，
Figure FDA0002182862630000039
表示t-Δt/2时刻处x方向纵波相关质点振动速度矢量，
Figure FDA00021828626300000310
表示t+Δt/2时刻处x方向横波相关质点振动速度矢量，
Figure FDA00021828626300000311
表示t-Δt/2时刻处x方向横波相关质点振动速度矢量，
Figure FDA00021828626300000312
表示t+Δt/2时刻处z方向纵波相关质点振动速度矢量，
Figure FDA00021828626300000313
表示t-Δt/2时刻处z方向纵波相关质点振动速度矢量，
Figure FDA00021828626300000314
表示t+Δt/2时刻处z方向横波相关质点振动速度矢量，
Figure FDA0002182862630000041
表示t-Δt/2时刻处z方向横波相关质点振动速度矢量，ρ1表示与参考网格点x方向相差一个空间步长处的密度值，ρ2表示表示与参考网格点z方向相差一个空间步长处的密度值，τp表示纵波应力，
Figure FDA0002182862630000042
表示横波相关x方向正应力，
Figure FDA0002182862630000043
表示横波相关切应力，
Figure FDA0002182862630000044
表示横波相关z方向正应力，
Figure FDA0002182862630000045
表示t+Δt时刻处横波相关z方向正应力，
Figure FDA0002182862630000046
表示t时刻处横波相关z方向正应力，
Figure FDA0002182862630000047
表示t+Δt时刻处横波相关x方向正应力，
Figure FDA0002182862630000048
表示t时刻处横波相关x方向正应力，
Figure FDA0002182862630000049
表示t+Δt时刻横波相关切应力，
Figure FDA00021828626300000410
表示t时刻处横波相关切应力，τp(t+Δt)表示t+Δt时刻处纵波应力，τp(t)表示t时刻处纵波应力，μ0及λ0表示参考网格点处的拉梅系数，μ1，2表示与参考网格点x方向和z方向均相差一个空间步长的拉梅系数，
Figure FDA00021828626300000411
及
Figure FDA00021828626300000412
分别表示与纵波相关且沿x坐标轴正负方向交错的空间偏导数，
Figure FDA00021828626300000413
及
Figure FDA00021828626300000414
分别表示与纵波相关且沿z坐标轴正负方向交错的空间偏导数，
Figure FDA00021828626300000415
及
Figure FDA00021828626300000416
分别表示与横波相关且沿x坐标轴正负方向交错的空间偏导数，
Figure FDA00021828626300000417
及
Figure FDA00021828626300000418
分别表示与横波相关且沿z坐标轴正负方向交错的空间偏导数。
7.根据权利要求1所述的方法，其特征在于，所述步骤S600具体为：
针对源检波场，采用矢量内积归一化方法，根据正向延拓和反向延拓的纵横波矢量波场获得PS及PP成像结果。
8.根据权利要求7所述的方法，其特征在于，针对源检波场，采用矢量内积归一化方法，按照下式根据正向延拓和反向延拓的纵横波矢量波场获得PS及PP成像结果：
Figure FDA00021828626300000419
式中，Ipp表示PP共成像点，Ips表示PS共成像点，Sp(x，t)为P波震源波场，Rp(x，t)为P波检波波场，Rs(x，t)为S波检波波场。
9.一种计算机存储介质，其特征在于，其中存储有用于实现上述权利要求1至8中任意一项所述方法的计算机程序。
10.一种计算机系统，其特征在于，包括存储器和处理器，所述处理器用于执行所述存储器中存储的计算机程序，所述计算机程序用于实现上述权利要求1至8中任意一项所述方法。
CN201910807274.4A 2019-08-28 2019-08-28 一种基于交错网格低秩有限差分的弹性逆时偏移成像方法 Pending CN112444849A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201910807274.4A CN112444849A (zh) 2019-08-28 2019-08-28 一种基于交错网格低秩有限差分的弹性逆时偏移成像方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201910807274.4A CN112444849A (zh) 2019-08-28 2019-08-28 一种基于交错网格低秩有限差分的弹性逆时偏移成像方法

Publications (1)

Publication Number Publication Date
CN112444849A true CN112444849A (zh) 2021-03-05

Family

ID=74742183

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201910807274.4A Pending CN112444849A (zh) 2019-08-28 2019-08-28 一种基于交错网格低秩有限差分的弹性逆时偏移成像方法

Country Status (1)

Country Link
CN (1) CN112444849A (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114660659A (zh) * 2022-03-30 2022-06-24 中国矿业大学 一种纵横波解耦的高精度微震定位方法

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104122585A (zh) * 2014-08-08 2014-10-29 中国石油大学(华东) 基于弹性波场矢量分解与低秩分解的地震正演模拟方法
CN104459773A (zh) * 2014-08-08 2015-03-25 中国石油天然气集团公司 基于交错网格Lowrank分解的无条件稳定地震波场延拓方法
CN108828659A (zh) * 2018-07-12 2018-11-16 中国石油天然气集团有限公司 地震波场延拓方法及装置
US20190018155A1 (en) * 2017-07-13 2019-01-17 Praveen Nakshatrala Visco-Pseudo-Elastic TTI FWI/RTM Formulation and Implementation

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104122585A (zh) * 2014-08-08 2014-10-29 中国石油大学(华东) 基于弹性波场矢量分解与低秩分解的地震正演模拟方法
CN104459773A (zh) * 2014-08-08 2015-03-25 中国石油天然气集团公司 基于交错网格Lowrank分解的无条件稳定地震波场延拓方法
US20190018155A1 (en) * 2017-07-13 2019-01-17 Praveen Nakshatrala Visco-Pseudo-Elastic TTI FWI/RTM Formulation and Implementation
CN108828659A (zh) * 2018-07-12 2018-11-16 中国石油天然气集团有限公司 地震波场延拓方法及装置

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
DONG HAN ET AL.: "Elastic wave propagation on a staggered-grid using lowrank finite-difference method", 《SEG INTERNATIONAL EXPOSITION AND 86TH ANNUAL MEETING》 *
韩冬: "基于Lowrank有限差分方法的弹性逆时偏移", 《中国优秀博硕士学位论文全文数据库（硕士) 基础科学辑》 *

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114660659A (zh) * 2022-03-30 2022-06-24 中国矿业大学 一种纵横波解耦的高精度微震定位方法

Similar Documents

Publication Publication Date Title
CA2726453C (en) Removal of surface-wave noise in seismic data
CA2726462A1 (en) Estimation of soil properties using waveforms of seismic surface waves
CN111025387B (zh) 一种页岩储层的叠前地震多参数反演方法
US20120016592A1 (en) Image domain signal to noise estimate with borehole data
CN112327358B (zh) 一种粘滞性介质中声波地震数据正演模拟方法
CN109946742B (zh) 一种TTI介质中纯qP波地震数据模拟方法
EP2513675A1 (en) Energy density and stress imaging conditions for source localization and characterization
Qu et al. Topographic elastic least‐squares reverse time migration based on vector P‐and S‐wave equations in the curvilinear coordinates
Shragge Acoustic wave propagation in tilted transversely isotropic media: Incorporating topography
Yao et al. Accurate seabed modeling using finite difference methods
Feng et al. Multiscale phase inversion for vertical transverse isotropic media
Li et al. Prestack waveform inversion based on analytical solution of the viscoelastic wave equation
CN112444849A (zh) 一种基于交错网格低秩有限差分的弹性逆时偏移成像方法
Ouyang et al. Multi-parameter true-amplitude generalized radon transform inversion for acoustic transversely isotropic media with a vertical symmetry axis
Liang et al. Born scattering integral, scattering radiation pattern, and generalized radon transform inversion in acoustic tilted transversely isotropic media
Ha et al. Seismic random noise attenuation in the laplace domain using singular value decomposition
CN115373022B (zh) 一种基于振幅相位校正的弹性波场Helmholtz分解方法
CN113687417B (zh) 一种三维叠前地震数据层间多次波预测和压制方法
Luo et al. A Born–WKBJ pre-stack seismic inversion based on a 3-D structural-geology model building
MX2011003852A (es) Atributos de procesamiento de imagen con inversion de tiempo.
Przebindowska Acoustic full waveform inversion of marine reflection seismic data
Tschache et al. On the accuracy and spatial sampling of finite-difference modelling in discontinuous models
Nita et al. Pseudo-depth/intercept-time monotonicity requirements in the inverse scattering algorithm for predicting internal multiple reflections
Shin et al. Laplace-domain waveform modeling and inversion for the 3D acoustic–elastic coupled media
Fang et al. Seismic wavefield modeling based on time-domain symplectic and Fourier finite-difference method

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
RJ01 Rejection of invention patent application after publication
RJ01 Rejection of invention patent application after publication

Application publication date: 20210305