井间地震波场数值模拟和弹性波逆时偏移*

2021-07-27李梦真

应用声学 2021年3期

关键词：井间检波器波场

李梦真何晓陈浩李超

(1 中国科学院声学研究所声场声信息国家重点实验室北京 100190)

(2 北京市海洋深部钻探测量工程技术研究中心北京 100190)

(3 中国科学院大学北京 100049)

0 引言

井间地震，顾名思义，就是在两口井中分别放置震源(激发信号)和检波器(接收信号)，通过对检波器获得的信号进行分析处理，来得到两井之间地质结构的一种地震勘探方法。井间地震的资料处理主要分为两类：一类是层析成像，主要利用接收信号中直达波的到时来得到井间的速度模型；另一类是反射波成像，它是在层析成像所得的速度场基础上，将接收信号中的反射波归位到地层的正确位置。与层析成像相比，反射波成像的垂直向分辨率更高、勘探区域更广，是井间地震资料处理中不可缺少的环节[1-2]。

井间地震反射波成像的算法可以分为两类。一类是利用波的几何射线理论，选用不同的道集并根据其特点进行成像的方法。其代表性的算法有XSP-CDP转化法[3]、CMD叠加法[4]以及共测向点叠加法[5]等。这类方法在实际中较为常用，对比较简单的地层结构能够获得好的效果。另一类是基于波动理论的偏移方法。相比于几何理论，波动方程理论更贴近于波的实际传播情况，能够对更为复杂的地层结构进行成像。许多的学者都对其进行了研究，如左建军等[6]将Gazdag 相移法应用到井间地震中；严建文等[7]针对井间地震资料的特点提出了一种计算效率高的炮域波动方程偏移法；张辉[8]实现了井间地震反射波成像技术的叠前深度偏移等。在常规的波动方程偏移法中，相移法适合于纵向变速而横向变速不大的介质[6]。而基于单程波波动方程偏移法能够很好地描述近垂直方向传播的波并适应速度变化剧烈的介质，但由于其在描述大角度传播的波时存在相位改变和振幅削弱的问题，导致其难以对陡倾角界面成像[9]，在本文算例中也验证说明了这一点。与这些方法相比，基于双程波的弹性逆时偏移方法虽然对计算性能和存储量要求较高[10]，但因其原理简单、成像精确、无反射倾角限制、适应任意复杂速度模型等优点[9]，近年来已受到了广泛关注。本文将针对井间地震弹性波逆时偏移的若干问题进行探讨，包括多层介质模型中井间地震信号的成分识别、波场分离、成像条件的对比以及对不同倾角反射体的成像结果对比等。

1 井间地震波场模拟与逆时偏移

井间地震弹性波的波场模拟是使用计算效率高的时域有限差分法(Finite-difference timedomain，FDTD)[11]来完成的。所用到的波场计算公式是弹性波的一阶速度-应力方程，区域边界产生的强反射采用完全匹配层(Perfect matched layer，PML)[12]进行吸收。以图1的多层地质模型为例，进行井间地震的数值模拟。

在图1模型中，将模型的左上角设为坐标原点，向右为x轴的正方向，表示水平的距离；向下为z轴的正方向，表示所处的深度。模型大小为300 m×600 m，网格大小为Δx=Δz=0.5 m，时间间隔约为57 μs，最大的记录时间是0.4 s。震源坐标在(0 m，300 m)处，如图1中左侧星号所对应的位置。震源使用200 Hz 的雷克子波，加载到Vx分量上。在模型的右侧均匀布满检波器，如图1右侧的点线。两个相邻检波器的间隔为2 m，共301道数据，并每隔3 道进行显示。模型的每一层地层参数在表1中都有说明。弹性波在震源处激发，经介质传播到达检波器。被接收的垂直速度分量Vz的示意图如图2所示。在图2中标出了主要的波场成分，其中横轴表示接收的时间，纵坐标是检波器被放置的深度。

图1 地质模型Fig.1 The geologic model

图2 接收的Vz 分量Fig.2 The Vz component received

表1 地层参数Table 1 The formation parameters

本文对数值模拟获得的人工合成信号进行处理和反射波偏移。首先是对接收到的信号进行波场分离。由于直达波的幅度较大，对界面成像起干扰作用，需要被去除。常用的方法是中值滤波，先将需要去除掉的直达波信号进行排齐，再对其进行中值滤波。用排齐后的信号减去中值滤波的信号后再解排，就可得到去除直达波后的波列，如图3。

图3 消除直达波后的Vz 分量Fig.3 The Vz component after eliminating direct wave

获得处理过的井间地震信号后，即可进行逆时偏移。逆时偏移是利用了波的时间一致性原理。当波传播到地层界面时会发生反射，在反射界面同时同点处产生入射波(上行波)和反射波(下行波)。逆时偏移就是根据这样的原理来对地层界面进行成像的。它的实现通常分为三步：一是震源波场的正向计算，来模拟地震当中的下行波场；二是检波器波场的逆向延拓，来模拟地震当中的上行波场；三是选用合适的成像条件，将上行和下行波场两者结合起来，来成像地层界面。其具体的实现方法可以参考文献[13-14]。

逆时偏移最常用到互相关成像条件，它的计算公式如下：

式(1)中，Sn代表第n个炮点的震源波场，Rn是第n个炮点的检波器波场，T为最大的记录时间，N代表最大的炮点数。image是成像的结果，(x，z)是坐标位置，t表示时间。应用互相关成像条件进行逆时偏移后得到单个炮点(N=1)Txx分量的成像结果，如图4所示。在图4中，可以看到单炮的逆时偏移是能够在界面位置(倾斜和水平界面)处正确成像的。图4的成像幅值整体较小且有附加噪声，针对这种情况可以通过多个炮点叠加成像的方法来进行解决，图1模型多炮叠加的结果对应图5(c)。

图4 单炮的逆时偏移结果Fig.4 The reverse time migration imaging result of single shot

2 不同成像条件的逆时偏移结果

逆时偏移的成像条件有多种，如激发时间成像条件、反褶积成像条件等。不同的成像条件，它的实现方式和成像效果是不同的[15-16]。在众多的成像条件中，互相关成像条件也就是式(1)，以其成像稳定和实现方便而获得了广泛的应用。在互相关成像条件的应用过程中，又衍生出了检波器波场归一化和震源波场归一化两种成像条件，分别对应式(2)和式(3)：

可以看到，归一化的成像条件就是在常规互相关成像的基础上用震源波场或检波器波场来对其进行归一化。为了防止分母为零的情况，也可以在分母处添加一个小数。由于互相关成像条件需要对波场进行积分，就使得其得到的成像值物理意义变得模糊，没有合适的量纲与之匹配。而使用不同的波场对互相关条件成像值进行归一化后，得到的结果就具有了量纲。归一化成像条件的量纲接近反射系数，从而使其结果具有了保幅性。这里对这3 种成像条件在井间地震中的成像效果进行探讨。

仍旧使用图1的地层模型和参数进行试验。所不同的是图4是震源在300 m 深度处的单炮偏移结果，而下面的偏移结果都是多炮叠加后的成像。在模型的左侧均匀放置震源。模型的垂直长度是600 m，每炮的间隔是10 m，共激发61 炮。检波器相邻两个的间距为2 m。每激发一炮共接收301 道数据并对所得的信号进行偏移成像，将61炮偏移成像的结果进行叠加，就得到了多炮逆时偏移的成像。将图1模型分别对不同的成像条件进行偏移，就得到了不同成像条件下的偏移结果，显示在图5。

图5分别是检波器归一化成像条件、震源归一化成像条件和常规互相关成像条件的结果。从图5(a)中可以看到靠近震源一侧的界面可以被成像出来，同时在震源一测产生了噪声。图5(b)和图5(c)相比，图5(b)中靠近检波器一侧的界面成像被保留和增强，而靠近震源附近的界面成像被减弱。这里对产生这种现象的原因进行解释：归一化的成像条件其实就是对互相关成像条件后的结果进行加权，只是加权的系数由震源或检波器波场进行提供。例如式(2)的分母，由检波器波场计算得到，可以将式(2)的分母看作是一幅在检波器一侧的像素值较大、而在震源附近幅度值较小的图像。用这样的一幅图与互相关成像的结果进行相除，可以使互相关成像中震源一侧的像被增强，而检波器一端的像被减弱，这就是造成在图5(a)中只成像出了靠近震源这侧的界面，而靠近检波器端界面的像几乎消失的原因。由于这种加权作用不仅对界面成像有效，同时也对互相关结果中单个炮点震源附近的波场进行增强，导致在图5(a)震源位置附近产生了噪声，这种噪声可以通过使用更多的炮点叠加来进行降低。同理，震源归一化成像条件会导致震源附近的像被减弱，检波器一侧的像被增强，如图5(b)所示。

图5 不同成像条件的对比Fig.5 Comparison of different imaging conditions

通过对图1模型的试算可以看到，无论是使用检波器波场还是震源波场来进行归一化，成像后都使对应波场这一侧的界面成像被削弱，而使另一侧的成像被增强；对于检波器波场归一化成像条件，成像结果相对较差，在成像出震源一侧的界面同时，也在震源附近产生了噪声。

3 井间地震狭缝模型的偏移成像

井间地震的发射和观测方式有别于地面地震，其震源和检波器的位置需分别置于目标地层的两侧。对于倾角很大的界面，如垂直方向的地层界面，在震源一侧激发的波场经过竖直地层反射后，反射波又回到了震源的一测，而没有传播到检波器端。这样就使得检波器获得的反射波信号大大减少，从而导致陡倾角的地层界面无法在井间地震中成像的结果。针对这样的情况，本文将设计狭缝的地层模型，来探讨其在井间地震中不同倾角下的逆时偏移成像问题。

所采用的模型如图6所示。图6(a)、图6(b)和图6(c)中的黑色区域是一条矩形的裂缝，被流体填充。假设流体的声速为1480 m/s，密度为1000 kg/m3。裂缝长度为150 m，宽度选取远小于震源波长的5 m 值(震源频率为200 Hz)。在图6(a)、图6(b)和图6(c)中的背景地层区域介质参数为：纵波速度3000 m/s，横波速度2000 m/s，密度2000 kg/m3。将裂缝放平且与x轴(水平)方向平行，对应的模型如图6(a)所示。将图6(a)的矩形狭缝以计算区域的中心为不动点对其进行顺时针方向的旋转，将狭缝与x轴所夹的锐角作为它的夹角。通过旋转不同的角度就得到了不同倾角下的狭缝模型，如图6(b)为倾角30°的狭缝，图6(c)是倾角90°的狭缝。

将模型计算区域的大小设为300 m×600 m。在模型的最左侧放置震源，每个震源间隔为20 m，从上到下依次共激发31 个炮点。在最右侧放置检波器，相邻两个检波器间隔为2 m，每一炮获取301道数据。在波场分离中，去除全波列中的直达纵波、直达横波以及纵波和横波的透射信号，而保留反射波以及在透射过程中产生的纵波和横波的转换波信号。通过对不同倾角的狭缝模型进行逆时偏移成像，得到图7的多炮偏移结果。数值模拟和成像实验表明，基于弹性波的逆时偏移可以获得0°～60°倾角狭缝的清晰成像。根据射线理论，当井间反射体倾角很大时，在接收井中获得的反射波信号非常少；但对于40°～60°倾角的狭缝，仍可获得与输入模型一致的偏移结果，而且成像效果与低倾角狭缝相当。推测此时在逆时偏移中，用作偏移资料的主要是折射过程中由纵波转换而成的横波以及由横波转换而成的纵波。这些转换波在逆时反向传播过程中，与正向传播的直达波发生互相关，从而在空间域准确获得了发生折射事件的像点。

图7 不同倾角的狭缝模型成像结果Fig.7 The imaging results of slit model with different obliquity

为了进一步说明这一点，本文还特别考察了夹角为90°的狭缝模型，即垂直狭缝的情况，对应图6(c)。在图6(c)模型中将震源放置在z= 200 m的深度处进行数值模拟，得到井间的接收信号，如图8所示，并在其上标注了主要的波场成分。图9是在图8的接收波列基础上去除了直达纵横波和透射P-P 波、S-S 波的结果。将处理后的信号逆时偏移，再对多个炮点叠加成像就得到了图10(b)的偏移结果。而图10(a)是基于单程声波方程F-K 偏移得到的结果。通过观察图10(a)和图10(b)可以看到，基于单程波声波方程的F-K 偏移法无法对狭缝垂直的边进行成像，而基于弹性波逆时偏移的方法就可以成像狭缝中间的部分。这说明对于高倾角的界面，基于单程波声波方程的F-K 偏移法是无法对其成像的，而弹性波逆时偏移仍然可以获得清晰的结果。这进一步体现了弹性波逆时偏移算法的优势。下面通过算例说明，应用于高倾角界面进行成像的波场成分是透射转换波。

图8 接收的Vz 分量Fig.8 The Vz component received

图9 消除直达波和P-P、S-S 折射波后的Vz 分量Fig.9 The Vz component after eliminating direct waves and P-P and S-S waves

图10 垂直狭缝在不同偏移方法下的结果Fig.10 Results of the vertical slit under different migration methods

先对图8中波场成分的名称进行说明。当震源在左侧产生弹性波(纵波和横波)后，经过狭缝中间部分的流体时都变为纵波，再穿透流体层后会产生转换波而被检波器接收。由于狭缝宽度远小于波长，在分析中可忽略弹性波在狭缝中的传播路径。对于穿过狭缝的纵波和横波，可写为P-P波和S-S波。将传播到狭缝之前的纵波和穿过狭缝后又变为横波的转换波，简写为透射P-S波；横波在穿过狭缝后也可产生透射P波，写为S-P波。如果在波场分离中再进一步去掉这些透射转换波，即在图9消去直达波和P-P 波、S-S 波的基础上再消去透射P-S 波和S-P波，其对应的结果显示在图11 上。再对处理过的信号做多炮叠加偏移成像，得到如图12的成像结果。

图11 消去P-S 和S-P 波后的结果Fig.11 The result after eliminating P-S and S-P waves

与图10(b)成像相比，在图12 中可以清晰地看到，狭缝两端点之间的像完全消失了，这证明了狭缝两端点间的像是由透射P-S 和S-P 转换波产生的。通过这样的计算实验，验证了弹性波逆时偏移在井间地震模型中可以获得任意倾角界面的成像结果。其中对于高倾角的界面，在地层界面处产生的透射转换波是成像的主要信号来源，这些透射转换波在波场分离时应该作为有用信号而保留。