一种高精度的逆时偏移成像方法和装置

1.本发明涉及一种逆时偏移成像方法和装置，尤其涉及一种高精度的逆时偏移成像方法和装置，属于石油勘探技术领域。


背景技术：

2.随着我国石油勘探向深层及超深层转移，深部目标层高陡复杂构造的高精度成像问题已成为当前制约我国石油勘探和开发迫在眉睫的瓶颈。成像技术虽在过去的数十年中取得了巨大的成功，但已无法满足高分辨率地震成像的需求。因此，亟需发展高精度的成像方法以满足复杂构造成像的需求。
3.波动方程逆时偏移成像(rtm)方法是一种高精度的成像方法。由于弹性波方法能够综合考虑纵横波对成像的贡献，弹性波逆时偏移成像方法被广泛研究和应用。
4.数值模拟是逆时偏移的关键，很大程度上决定着偏移成像的精度。目前在弹性波逆时偏移中普遍使用的交错网格有限差分模拟，而常规的交错网格有限差分求解方法尽管在空间离散上可以实现任意偶数阶精度，但是在时间离散方向上仅有二阶精度。低的时间模拟精度会导致数值模拟中产生严重的时间频散，诸如在波前面出现同相轴超前，同相轴增多的现象。同时，常规方法逆时偏移成像方法的稳定性条件比较苛刻，也会影响弹性波逆时偏移的计算效率。
5.因此，如何提供一种能够同时压制时间和空间频散的高精度逆时偏移成像方法是本领域技术人员亟需解决的技术难题。


技术实现要素：

6.本发明要解决的技术问题是提供一种高精度的逆时偏移成像方法和装置。
7.为解决上述技术问题，本发明所采取的技术方案一是：
8.一种高精度的逆时偏移成像方法，包括以下步骤：
9.步骤1：建立纵横波解耦形式的准速度应力方程：
[0010][0011]
其中：ρ为密度；α和β分别为p波速和s波速；v
x
和vz分别为x方向的质点振动的速度分量和z方向的质点振动的速度分量；τ
111
、τ
11
、τ
22
和τ
12
分别为质点沿着四个不同方向振动的应力分量；为一阶时间微分算子,r＝x,z表示空间r方向求导算子；和为与纵波相关的质点沿着x和z方向振动的速度分量，和为与横波相关的质点沿着x和z
方向振动的速度分量，j＝x,z分别为与纵波和横波相关的空间方向求导算子；
[0012]
步骤2：求解纵横波解耦形式的准速度应力方程：
[0013]
利用采集的震源数据作为已知数据，求解步骤1建立的纵横波解耦形式的准速度应力方程，得到正传波场和
[0014]
再以监测点监测的地震数据作为已知数据，求解步骤1建立的纵横波解耦形式的准速度应力方程，得到反传波场和
[0015]
采用有限差分算法求解所述纵横波解耦形式的准速度应力方程，求解过程中使用的差分模板为时间四阶空间任意偶数阶差分模板；
[0016]
步骤3：计算成像结果：
[0017]
对rtm的点积成像条件进行归一化成像：
[0018][0019][0020][0021][0022]
进一步，步骤2中所述的时间四阶空间任意偶数阶差分模板为：
[0023][0024][0025][0026][0027][0028]
其中，u表示波场变量，可以为τ
111
、τ
11
、τ
22
、τ
12
、和nr表示差分模板单边长度，hr和h
φ
表示r和φ方向的网格步长，r表示与导数相关的坐标轴，φ表示平行于r轴的坐标轴；γ
r,χ
＝χδt/hr，γ
φ,χ
＝χδt/h
φ
；δt表示时间采样间隔；χ分别取α和β。
[0029]iαα
(x)和i
αβ
(x)分别表示pp波图像和ps波图像，其中p表示纵波，s表示横波。
[0030]
本发明所采取的技术方案二是：
[0031]
一种高精度的逆时偏移成像设备，包括纵横波解耦形式的准速度应力方程构建求解模块、正传波场模块、反传波场模块和成像计算模块；纵横波解耦形式的准速度应力方程构建求解模块用于构建求解所述纵横波解耦形式的准速度应力方程；
[0032]
震源数据输入正向波场模块，由正向波场模块调用纵横波解耦形式的准速度应力方程构建求解模块获得正向波场；
[0033]
监测点观测数据输入反传波场模块，由反传波场模块调用纵横波解耦形式的准速度应力方程构建求解模块获得反向波场；
[0034]
成像计算模块根据解耦的正传波场、反传波场和矢量成像条件计算纵横波成像结果。
[0035]
采用上述技术方案所产生的有益效果在于：
[0036]
(1)本发明构建了纵横波解耦形式的准速度应力方程，在应力波场不解耦的前提下，实现仅解耦速度波场的分解，显示的解耦方程编程易于实现且不会显著增加计算量。
[0037]
(2)采用时间四阶和空间任意偶数阶有限差分方法，能够有效地抑制传统时间二阶有限差分方法中出现的时间色散和波形畸变问题，从而从数值模拟的角度保证了高精度的成像结果。
[0038]
(3)本发明采用时间四阶和空间任意偶数阶有限差分方法求纵横波解耦形式的准速度应力方程，有着更宽松的稳定性条件，并且解耦的速度波场，可以简便地实现偏移成像中的点积成像条件的归一化成像。
附图说明
[0039]
图1是本发明的流程图；
[0040]
图2是本发明中时间四阶空间任意偶数阶差分模板示意图；
[0041]
图3是本发明实施例1中计算波场结果图；
[0042]
图4是本发明实施例1中参考波场图；
[0043]
图5是本发明实施例1中计算波场结果、常规时间二阶精度差分模拟结果与参考波场的剖面对比图；
[0044]
图6是本发明实施例2中pp波成像结果；
[0045]
图7是本发明实施例2中ps波成像结果；
[0046]
图8是本发明实施例2中计算结果与常规方法在x＝4km处pp波成像结果对比图；
[0047]
图9是本发明实施例2中计算结果与常规方法在x＝10km处ps波成像结果对比图。
具体实施方式
[0048]
实施例1：
[0049]
一种高精度的逆时偏移成像方法，包括以下步骤：
[0050]
步骤1：建立纵横波解耦形式的准速度应力方程：
[0051][0052]
其中：ρ为密度；α和β分别为p波速和s波速；v
x
和vz分别为x方向的质点振动的速度分量和z方向的质点振动的速度分量；τ
111
、τ
11
、τ
22
和τ
12
分别为质点沿着四个不同方向振动的应力分量；为一阶时间微分算子,r＝x,z表示空间r方向求导算子；和为与纵波相关的质点沿着x和z方向振动的速度分量，和为与横波相关的质点沿着x和z方向振动的速度分量，j＝x,z分别为与纵波和横波相关的空间方向求导算子；
[0053]
步骤2：求解纵横波解耦形式的准速度应力方程：
[0054]
利用采集的震源数据作为已知数据，求解步骤1建立的纵横波解耦形式的准速度应力方程，得到正传波场和
[0055]
再以监测点监测的地震数据作为已知数据，求解步骤1建立的纵横波解耦形式的准速度应力方程，得到反传波场和
[0056]
采用有限差分算法求解所述纵横波解耦形式的准速度应力方程，求解过程中使用的差分模板为时间四阶空间任意偶数阶差分模板；
[0057]
步骤3：计算成像结果：
[0058]
对rtm的点积成像条件进行归一化成像：
[0059][0060][0061][0062][0063]
进一步，步骤2中所述的时间四阶空间任意偶数阶差分模板为：
[0064]
[0065][0066][0067][0068][0069]
其中，u表示波场变量，可以为τ
111
、τ
11
、τ
22
、τ
12
、和nr表示差分模板单边长度，hr和h
φ
表示r和φ方向的网格步长，r表示与导数相关的坐标轴，φ表示平行于r轴的坐标轴；γ
r,χ
＝χδt/hr，γ
φ,χ
＝χδt/h
φ
；δt表示时间采样间隔；χ分别取α和β可得到。不同的模板，分别用于计算纵波和横波的空间倒导数。i
αα
(x)和i
αβ
(x)分别表示pp波图像和ps波图像，其中p表示纵波，s表示横波。
[0070]
本实施例中采用了双层介质模型进行试算。模型划分为1500
×
500的网格，网格间隔为20m。第一层的p波和s波的速度分别为3.0km/s和1.8km/s,第二层p波和s波的速度分别为4.0km/s和2.4km/s。差分计算的空间阶数2nr＝8，时间间隔δt＝2ms，时间阶数为四阶；其在3.35秒处的波场如图3所示。参考波场如图4所示，其结果通过采用二阶时间导数计算，δt＝0.5ms，2nr＝12。图5对比了本实施例模拟结果、常规时间二阶精度差分模拟结果与参考波场的剖面，本实施例相对参考波场的均方根误差为0.0025，而常规时间二阶精度差分模拟结果相对参考波场的均方根误差为0.042。本实例证明了本发明采用的高精度模拟方法能显著地提升数值模拟精度。
[0071]
实施例2：
[0072]
本实施例应用于偏移成像。选取的公开模型为加拿大山丘模型。本实施例中将该模型划分为了500
×
834的网格，网格间隔为20m。空间阶数2nr＝10，时间间隔δt＝1.8ms，时间阶数为四阶。参考波场采用二阶时间导数计算，δt＝0.45ms，2nr＝12。获得的pp波和ps波成像结果分别见图6和图7。为了说明本发明的成像精度，将本实施例结果、常规时间二阶精度的成像结果与参考解做比较。位于x＝4km的pp波成像结果和位于x＝10km的ps波成像结果对比分别见图8和图9。对比不难发现本发明能够得到与参考解高度一致的成像结果，然而使用二阶时间精度常规方法在成像位置会引入了可见成像误差，并产生了虚假。
[0073]
实施例3：
[0074]
一种采用实施例1或2所述高精度的逆时偏移成像方法的成像设备，包括纵横波解耦形式的准速度应力方程构建求解模块、正传波场模块、反传波场模块和成像计算模块；纵横波解耦形式的准速度应力方程构建求解模块用于构建求解所述纵横波解耦形式的准速度应力方程；
[0075]
震源数据输入正向波场模块，由正向波场模块调用纵横波解耦形式的准速度应力方程构建求解模块获得正向波场；
[0076]
监测点观测数据输入反传波场模块，由反传波场模块调用纵横波解耦形式的准速
度应力方程构建求解模块获得反向波场；
[0077]
成像计算模块根据解耦的正传波场、反传波场和矢量成像条件计算纵横波成像结果。