逆时偏移方法是基于双程波波动方程, 通过对波场正向和反向延拓来实现偏移。逆时偏移是较Kirchhoff偏移方法而言能有效利用全波场信息, 较单程波偏移方法而言, 无传播方向限制, 不需要对波场进行上、下行波的分离, 而且不受倾角限制, 并能较好的成像回转波、多次波, 理论上在给定的速度场足够准确的情况下能对任意复杂介质精准成像, 是一种具有极高成像精度的偏移方法^[2]。

但逆时偏移有着与生俱来的三大缺陷:计算量大、存储量大和存在虚假成像^[3]。如何有效压制虚假的噪声是获得正确逆时偏移成像结果的关键所在。根据互相关成像条件公式^[4]

I(x, z)=∫0tmaxS(x, z, t1)R(x, z, t-t2)dt, (1)t=t1+t2。

只要满足:震源正向延拓波场的旅行时间t₁与检波点反向延拓波场的旅行时间t₂之和, 等于在从震源出发至检波点的总旅行时间t, 就满足互相关成像条件, 因此都可以成像。如图1所示的两条射线的路径上, 任何一点的走时之和都是等于t的, 因此在这两条射线的每个点上, 都满足互相关成像条件, 在偏移时都会成像, 而真实的成像点只在反射界面处, 这就是低频噪声产生的原因。

图1

	Figure Option ViewDownloadNew Window
	图1 逆时偏移低频噪声产生的原因

图1中方框处是震源正向延拓波场的与检波点反向延拓波场夹角等于180° 的成像噪声产生的位置, 圆圈处是震源正向延拓波场的与检波点反向延拓波场夹角小于180° 的成像噪声产生的位置^[5]。因此低频噪声分为两只种:第一种是波场夹角等于180° 的成像噪声, 为图2白色箭头所指的位置, 这种噪声振幅最强, 干扰严重。第二种是波场夹角小于180° 的成像噪声, 为图2黑色箭头所指的位置, 这种噪声振幅较弱, 一般有两条。

图2

	Figure Option ViewDownloadNew Window
	图2 水平模型逆时偏移结果

去除低频噪声的方法主要可以分为两类^[6]。第一类是滤波, 针对噪声低频的特点直接成像后进行滤波处理, 如高通滤波^[7], 拉普拉斯滤波^[8], 角度滤波^{[9, 10, 11]}等, 杨仁虎等^[12]提出了一种振幅补偿滤波方法, 在滤波前后, 尽可能保持有效振幅。Guitton^[13]比较了几种滤波方法, 提出最小平方滤波方法。第二类是控制方向成像。只在正确的成像点处有选择性的进行成像, 如Baysal^[14]和Etgen^[15]假设波阻抗为常数推导出了无反射波动方程; Fletcher^[16]指出无反射波动方程只有当入射波垂直入射时, 其反射能量才为零, 在波动方程中加入方向衰减因子, 压制反向的反射; Yoon^[17]利用坡印廷矢量判断波场方向, 只对相反方向的波场进行成像, 压制了大角度的虚假成像。文中重点分析了波场分离滤波方法并比较了坡印廷矢量滤波方法和拉普拉斯滤波方法三种滤波方法的去噪效果。

Liu^[18]提出通过一种波场分离的方法, 其主要思路是将波场拆分成不同方向的单成波分量, 在成像时对波场分量进行人工选择, 提取所需要的波场分量, 只在有效路径上进行互相关成像, 从而达到去除噪声的目的。

在逆时偏移中使用的是双程波动方程, 它包含了上行波和下行波两个组成部分, 如果在z方向进行上、下行波的分解, 则震源波场和检波点波场可分解为

S(x, t)=Sz+(x, t)+Sz-(x, t), (2)R(x, t)=Rz+(x, t)+Rz-(x, t)。(3)

其中:“ +” 号代表下行波场, “ -” 号代表上行波场。将式(2)、式(3)代入互相关成像条件式(1)中, 得到新的成像条件

I(x, t)= ∫0tmaxSz+(x, t)Rz-(x, t)dt+∫0tmaxSz-(x, t)Rz+(x, t)dt+∫0tmaxSz+(x, t)Rz+(x, t)dt+∫0tmaxSz-(x, t)Rz-(x, t)dt=Iz1+Iz2+Iz3+Iz4。(4)

经过真实的反射界面时, 入射波场和反射波场的传播方向在界面的法线方向会发生变化。因此震源波场和检波点波场在反射界面上都会产生相对于界面向上传播的波和相对于界面向下传播的波。按照震源波场和检波点波场的传播方向, 将其波场分离成其各自的传播方向不同的波场, 只有两者传播方向相反的波场才会在成像时产生真实的成像点。也就是说, 只需要I_z1和I_z2这两项就可以实现正确的逆时偏移成像。在相同传播方向所成的像是虚假的成像点, 因此I_z3和I_z4这两项其成像结果是低频噪声, 所以应该将其删去。这样就从成因上消除了低频噪声。

因此可以将式(4)改写成新的成像条件

I(x, z)=∫0tmaxSz+(x, z, t)Rz-(x, z, t)dt+∫0tmaxSz-(x, z, t)Rz+(x, z, t)dt。(5)

为了成功的应用这个成像条件, 需要将震源波场和检波点波场拆分成上行波场分量和下行波场分量。文中采用的是二维傅里叶变换方法, 先将震源和检波点波场转换到频率波数域, 再对其进行上、下行波分离, 然后反变换到时间域, 这样就可以成功的利用波场分离的方法来压制逆时偏移噪声。

图3是应用波场分离法去噪后的成像效果图。波场分离方法去噪彻底, 绝大部分噪声得到了压制, 去噪之后图像界面清晰。波场分离方法从成因上消除了虚假噪声, 因此完全可以压制波场夹角等于180° 的噪声。但是两条“ 弧形噪声” 只得到了部分压制, 这些对于波场夹角小于180° 的噪声无法完全去除, 因为这些波场的传播路径会和一些有效波场一样, 如果强行去除反而还会破坏真实成像信息。

图3

	Figure Option ViewDownloadNew Window
	图3 波场分离法去噪后成像效果

只进行水平模型的计算, 缺乏实用价值, 无法充分说明逆时偏移的效果, 还需要对标准模型进行测试。但是复杂标准模型的偏移一般需要进行几百炮集以上的计算, 传统的串行计算肯定无法满足需求, 因此需要选用高效率的计算方法。本文采用的是基于消息传递接口(MPI)的并行计算方案^[19]。

文中选用的Marmousi2模型总网格大小为972× 701, 空间网格横向步长为8 m, 纵向网格步长均为5 m, 时间步长为4 ms, 模型最小速度为1 028 m/s, 最大速度为4 700 m/s, 密度都为1 g/cm³, 选取的参数满足有限差分计算所需的稳定性条件。交错网格有限差分的差分精度为时间二阶、空间十阶。采用右端放炮、左端排列接收的观测方式, 每炮的计算区域大小为240× 701, 检波器分布在整个计算区域表面, 最小偏移距为0 m, 每炮接收240道记录, 道间距为8 m, 炮点在模型的(240, 1)处开始右移, 炮间距为24 m。选用主频25 Hz的零相位雷克子波在地表激发, 最大接收时间为3 s, 按照此观测方式可以得到240个地面检波点炮集记录。

如果用串行计算(CPU为AMD FX-8320, 主频3.50 GHz ), 一个炮集的偏移时间为100 min, 如果使用串行计算, 则需要花17 d的时间。为了减少计算时间, 文中应用多核计算机(CPU为AMD FX-8320八核处理器, 主频3.50 GHz)进行计算, 使用其中7核Mpi并行计算(留出1核计算空间处理系统资源), 只需要55 h就可以完成全部240炮的偏移。如果采用更高端的并行机计算, 时间可以进一步缩短。

图4为去噪前Marmousi2模型逆时偏移成像结果, 低频噪声基本覆盖整个模型, 尤其模型上部高振幅的假象最为严重, 遮掩了下部构造, 影响了成像质量。

图4

	Figure Option ViewDownloadNew Window
	图4 去噪前Marmousi2模型成像结果

图5是用了三种不同方法去噪结果。通过细节对比, 用坡印廷方法和拉普拉斯方法去噪后, 噪声去只得到了一定程度上的压制, 反射界面上的噪声还有残留。而使用波场分离方法后成像噪声压制情况良好, 波场夹角等于180° 的噪声完全得到了压制, 低频噪声的影响基本消除。可以看到, 坡印廷方法、拉普拉斯方法、波场分离法的去噪效果越来越好(注意图5中三种方法圆圈处的噪声压制之比较)。经过波场分离方法去噪之后, 地层反射界面显示更准确, 含气透镜体模型、三条大的断层和下部的高速的横向侵入岩更清晰, 地层特征细节更为明显, 成像质量得到了明显的改善。

图5

	Figure Option ViewDownloadNew Window
	图5 三种不同去噪方法后成像结果对比 a— 坡印廷方法去噪后成像结果; b— 拉普拉斯去噪后成像结果; c— 波场分离法去噪后成像结果

基于波场分离的逆时偏移去噪方法从成因上消除了虚假噪声, 在成像时可以完全压制波场夹角等于180° 的噪声。通过波场分离方法在Marmousi2标准模型上的测试, 证明了这种方法的有效性, 对比其他去噪方法表明在噪声的压制效果上该方法更具优势。不过要运用两次二维傅里叶变换一定程度上也增加了计算量和存储量, 但是采用并行计算或CPU/GPU计算平台可以满足计算需求。

在陡倾角复杂的地区如果只分离上、下行波, 可能无法有效的将真实反射点与虚假成像点分离开, 这时需要在x和z方向上进行两次傅里叶变换, 将波场分离成上、下、左、右行波。