技术领域
[0001] 本
发明涉及地球物理勘探，尤其是涉及一种基于高阶空间差分方法的弹性波方程逆时偏移成像方法，适用于复杂地下结构的高
精度成像，并能够高效快速求解大尺度问题。
背景技术
[0002] 逆时偏移(Reverse time migration,RTM)方法是最近10多年来油气勘探领域的研究热点之一，也是目前最高精度的
地震波成像方法之一。与其他成像方法相比，逆时偏移成像方法具有成像精度高、适应性广等特点，尤其是在高倾
角、
水平速度变化大的复杂地质结构成像中具有明显的优势，能够对常规地震成像方法中无法成像的复杂地质结构进行高精度成像。逆时偏移的基本思想最早在1978年由Hemon提出[1]。随后，Baysal[2]，Whitmore[3]及Mcmechan[4]别提出了对逆时偏移思想的具体应用。近年来，逆时偏移成像的研究热点不断发展，从叠后偏移逐步发展到叠前偏移，从2D成像推广到3D成像，从
声波方程偏移扩展到弹性波方程偏移，从各项同性介质扩展到倾斜横向各向同性介质(tilted transversely isotropic media,TTI)和其它
各向异性介质，都取得了一系列的重要成果。
[0003] 在逆时偏移成像中所需的庞大计算资源仍然是制约逆时偏移成像技术的广泛应用，迫切需要提高逆时偏移的效率，降低计算成本。要缓解计算
瓶颈，可以考虑使用高效快速
算法来更快更好的求解
波动方程，从而减少所需的计算资源。有限差分方法(finite difference method，FDM)是用于求解偏微分方程最常用的数值方法[5]。而时域有限差分方法(Finite difference time domain method，FDTD)的基本思想是用差分方法来近似时间域和空间域的偏微分，并通过时间域的递推来模拟波的传播过程。FDTD的计算需要涉及到网格剖分、差分格式、差分阶数等问题。对于逆时偏移成像，时域有限差分方法是最常用的技术方法。
[0004] IO带宽瓶颈也是逆时偏移成像中需要考虑的重要问题。传统的逆时偏移成像需要把其中一个波场在所有时刻的快照保存到磁盘，在成像过程中再读取出来，这很大程度上成为整个成像过程的瓶颈。为了解决这个问题，采用Clapp[7]和Dussaud[8]提出波场保存方法，该方法先将源波场进行外推，每个时间步只保存源波场中靠近计算区域边界的若干层，有效地缓解了IO瓶颈。
发明内容
[0005] 本发明的目的在于提供基于高阶差分方法的弹性波方程逆时偏移成像方法。
[0006] 本发明包括以下步骤：
[0007] 1)使用初始模型，在
指定炮点放置
震源信号，并应用一阶速度-应
力方程进行数值仿真，对震源波场进行正向波场延拓直到最大时间步，即从t0到tmax，同时保存每个时间步的边界波场(根据模型大小，每个方向各保存3～5层不等)；
[0008] 在步骤1)中，所述一阶速度-
应力方程，是指质点速度矢量v和应力张量τ的一阶弹性介质偏微分方程，在非均匀各向同性弹性介质中，介质差异主要通过介质
密度ρ(r)、
拉梅常数λ(r)和μ(r)来进行描述；采用了高阶时域有限差分方法作为数值仿真方法，实现4阶、6阶、8阶、10阶和12阶空间差分格式，对于高阶差分方法，将求导算子Dη作用在函数u(jη)(η＝x,y,z)上，定义如下：
[0009]
[0010] 其中，order为FDTD的空间阶数。
[0011] 2)使用初始模型，并利用保存的边界波场作为边界条件，使用数值方法仿真，对震源波场进行逆时外推，即从tmax到t0；
[0012] 仿真使用的公式和仿真方法同步骤1)；
[0013] 3)使用初始模型，使用
检波器波场作为边界条件注入，使用数值方法仿真，对检波器波场进行逆时外推，即从tmax到t0。
[0014] 仿真使用的公式和仿真方法同步骤1)；检波器信号在现场实验中获得或使用真实模型仿真产生；
[0015] 4)对步骤2)、3)中的震源波场和检波器波场应用震源补偿归一化互相关成像条件；所述震源补偿归一化互相关成像条件如下：
[0016]
[0017] 其中S(x,y,z,t)是震源波场，R(x,y,z,t)是检波器波场，time从t0到tmax是时间步数。
[0018] 与现有的逆时偏移成像方法相比，本发明具有以下技术效果：
[0019] 1、本发明使用高阶时域有限差分方法进行数值仿真。数值仿真是逆时偏移成像的最重要组成部分。逆时偏移成像需要对双程波动方程进
行波场延拓，求解过程中需要同时实现震源波场和检波器波场的数值
迭代和互相关，这都需要使用数值方法来实现。因此，数值方法的实现很大程度上决定了逆时偏移成像的效率和
质量。高阶方法与传统的二阶差分方法相比，不仅能保证计算所需的精度，而且能有效的降低空间和时间
采样率，从而节省存储空间、提高计算效率。本发明中还采用了Holberg提出的对FDTD差分系数进行优化的方法，能进一步程度上减少对存储的需求。
[0020] 2、针对成像过程中所需的存储资源和输入输出带宽。先将源波场进行正向外推，每个时间步只保存源波场中靠近计算区域边界的若干层，然后再把源波场逆推，并把之前保存的边界波场注入以恢复源波场。这种方法有效的减少了所需的存储空间并缓解了IO瓶颈。
附图说明
[0022] 图2为真实发层模型(P波)。
[0023] 图3为真实分层模型(S波)。
[0024] 图4为初始分层模型(P波)。
[0025] 图5为初始分层模型(S波)。
[0026] 图6为分层模型成像结果。
[0027] 图7为Marmousi II真实模型(P波)。
[0028] 图8为Marmousi II真实模型(S波)。
[0029] 图9为Marmousi II初始模型(P波)。
[0030] 图10为Marmousi II真实模型(S波)。
[0031] 图11为Marmousi II模型成像结果。
具体实施方式
[0032] 以下将结合附图对本发明的技术方案及其突出效果作进一步说明。
[0033] 图说明本发明中逆时偏移成像方法的主要流程。
[0034] 1)使用初始模型，在指定炮点放置震源信号，并应用一阶速度-应力方程进行数值仿真，对震源波场进行正向波场延拓直到最大时间步，即从t0到tmax，同时保存每个时间步的边界波场(根据模型大小，每个方向各保存3～5层不等)。
[0035] 2)使用初始模型，并利用保存的边界波场作为边界条件，使用数值方法仿真，对震源波场进行逆时外推，即从tmax到t0；
[0036] 3)使用初始模型，使用检波器波场作为边界条件注入，使用数值方法仿真，对检波器波场进行逆时外推，即从tmax到t0。
[0037] 检波器信号在现场实验中获得或使用真实模型仿真产生。
[0038] 4)对步骤2)、3)中的震源波场和检波器波场应用震源补偿归一化互相关成像条件。
[0039] 本发明可应用于非均匀各向同性介质中的弹性波逆时偏移成像。经实践证明，采用高阶差分方法实现的逆时偏移成像能大幅度降低所需的存储资源和计算时间，并能得到高精度的成像结果。
[0040] 图2～5为一个简单分层模型的P波和S
波速度分布。图2、图3是真实模型，包含4个水平分层和3个不同形状的反射体；图4、图5为仿真时使用的初始模型，仅包含4个水平分层无反射体。
[0041] 图6是采用本发明中的逆时偏移成像方法，采用十二阶差分(Taylor展开)在最低空间采样下得到的成像结果。仿真时使用图4、图5所示的速度模型作为初始模型，并使用中心
频率为20Hz的Ricker子波作为震源信号。
[0042] 从结果来看，各主要的反射体都得到了较清晰的图像。成像效果良好。这说明高阶差分方法的逆时偏移可以在低采样下保持成像精度，从而大幅度提高成像效率，降低成像的成本。
[0043] 使用12阶差分(Taylor展开)与使用2阶差分相比，空间采样率降低到了三分之一。对于这样一个二维问题，每个维度降低到三分之一，两个维度
叠加就下降到了九分之一，节省了九分之八，即88.9％的存储空间；十二阶差分方法所需的计算时间步数降低到传统二阶中心差分方法的44.6％。使用CPU单核计算时，十二阶差分方法(Taylor展开)与二阶差分方法相比，总体计算时间在二维空间问题中可节省约89％。
[0044] 图至图为Marmousi II模型的P波和S波速度分布。图、图是真实模型；图、图为仿真时使用的初始模型，由真实模型各点在周围200米矩形区域内取平均值得到。
[0045] 图是采用本发明中的逆时偏移成像方法，采用十二阶差分(Holberg优化系数)并在最低空间采样下得到的成像结果。仿真时使用图、图所示的速度模型作为初始模型，并使用中心频率为5Hz的Ricker子波作为震源信号。
[0046] 从结果来看，各主要的反射层都得到了较清晰的图像，包括位于速度模型深处的弱反射层。总体成像效果良好。经实践证明，采用高阶差分方法实现的逆时偏移成像能大幅度降低所需的存储资源和计算时间，并能得到高精度的成像结果。
[0047] 从效率上来看，使用12阶差分(Holberg优化系数)与使用传统的2阶中心差分相比，空间采样率降低到了24.3％(约四分之一)。对于这样一个二维问题，每个维度降低到四分之一，两个维度叠加就下降到了十六分之一，节省了93.75％的存储空间；十二阶差分方法所需的计算时间步数降低到传统二阶中心差分方法的35.4％。使用CPU单核计算时，十二阶差分方法与二阶差分方法相比，总体计算时间在二维空间问题中可节省约93％。
