混合域高分辨率Radon变换及其在绕射波分离与成像中的应用

2020-11-24罗腾腾徐基祥孙夕平

石油物探 2020年6期

罗腾腾,徐基祥,秦臻,孙夕平

(1.中国石油勘探开发研究院,北京100083;2.三峡大学,湖北宜昌443002)

常规的偏移算法主要记录由连续反射层或主要不连续点(大断层)定义的“高能事件”,它是地震解释人员通常使用的地震信息,常被称作“镜像”能量。随着油气勘探开发研究的深入,勘探目标逐步转向断层、裂缝、河道、粗糙岩丘边缘等地下小尺度不连续地质体[1],这些小尺度地质体的地震响应通常表现为能量较弱的绕射波。

20世纪50年代,KREY[2]首先利用绕射波信息来研究地下地质异常体如断层、断点等。KLEM-MUSATOV等[3]指出,绕射波振幅通常比反射波振幅弱一到两个数量级,且在常规处理流程中多以反射波作为有效信号,而那些高分辨率、低能量的绕射波信号常被破坏或压制。为了充分利用地震记录中的绕射波信息来刻画地下非均质地质体,需要在全波场记录中将绕射波场分离并单独成像[4]。基于反射波和绕射波在运动学和动力学方面的特征差异[5-6],国内外学者提出了许多不同的绕射波分离成像方法。按照绕射波分离成像的方法原理以及处理手段的差异可将现有方法分为:叠加、滤波、平面波解构滤波、聚焦-反聚焦、Radon变换-反Radon变换等。

1) 叠加方法:LANDA等[7]、KANASEWICH等[8]、TSINGAS等[9]提出在共偏移距道集和共断层点道集上使用绕射波走时曲线来描述绕射波相干的总和。LANDA等[10]利用在共绕射点剖面上叠加来自绕射点处的地震信号识别局部非均质地质体。

2) 滤波方法:KOZLOV等[11]进一步修改Kirchhoff偏移算子中的锥形权函数使其更好地压制反射能量,对地震响应中的弱绕射波信息进行准确成像;BANSAL等[12]讨论了多种分离反射和绕射波场的方法,主要包括共偏移距倾角滤波、正常时差校正倾角滤波、正常时差-倾角时差校正倾角滤波、特征向量滤波、Radon滤波等,并由合成记录阐述了不同分离方法的优缺点;MOSER等[13]在Kirchhoff偏移积分公式中通过修改加权函数来实现反稳相绕射波提取;KOREN等[14]提出一种绕射波成像新方法,即利用方向角分解设计加权叠加滤波器实现各向同性/各向异性介质中绕射波和反射波的分离,同时提高连续构造和小尺度地质体的成像分辨率;李晓峰等[15]在传统的Kirchhoff偏移算子中引入反稳相滤波器来压制反射波能量,凸显绕射波能量,并在成像过程中引入校正极性算子实现绕射波准确成像;刘培君等[16]借助相关分析以及高斯束偏移的射线追踪构建反稳相滤波算子,采用反稳相偏移得到绕射波成像剖面。

3) 平面波解构滤波方法:CLAERBOUT[17]最先提出平面波解构滤波器的概念,并于1994年[18]进一步完善了该滤波器方法;FOMEL[19]改进了平面波解构滤波器,使所需唯一参数为局部平面波场的斜率;TANER等[20]模拟了平面波剖面和平面波分解滤波器,并依据反射波和绕射波在平面波地震记录上的同相轴连续性差异在叠前道集中进行绕射波分离成像;黄建平等[21]对平面波分解(PWD)解构滤波器使用方法及其在叠前、叠后道集上分离成像绕射波进行了系统总结;孔雪等[22]和刘玉金等[23]证实平面波剖面本质上具有区分绕射波和反射波特性,并利用平面波解构滤波技术压制反射波来获得绕射波;朱生旺等[24]利用局部倾角滤波技术和预测反演技术修正绕射波场分离时低倾角信息估计不足对平面波分解滤波器的影响,从而提高绕射波成像结果的横向分辨率;DECKER等[25]在叠前偏移倾角道集上对每个固定倾角的成像剖面应用PWD解构滤波技术并叠加绕射波场获得单独的绕射波成像剖面。

4) 聚焦-反聚焦方法:KHAIDUKOV等[26]在叠前共炮点道集将反射波聚焦于镜像虚震源点,而绕射波仍保持发散的状态,利用聚焦切除-反聚焦方法成功提取断层、断点等小尺度地质体,并对其成像;MOSER等[13]提出两种绕射波分离成像的方法,一种方法是在叠前深度域利用聚焦滤除反射波实现绕射波分离与成像,另一种方法是修改偏移算法的核函数,采用反稳相偏移压制反射能量;BERKOVITCH等[27]利用绕射多次聚焦叠加的方法将选取的绕射波最优同相轴叠加从而得到主要包含绕射波的叠加剖面;RESHEF等[28]利用角道集上的绕射能量进行高分辨率速度分析。

5) Radon变换-反Radon变换方法:KLOKOV等[29]根据反射波同相轴在倾角域CIGs中表现为开口向上的“笑脸”状,利用混合域Radon变换识别反射波顶点并消除反射能量,增强绕射信号;KLOKOV等[30]推导了偏移倾角域反射波和绕射波的解析表达式,并利用混合域Radon变换实现绕射波场和反射波场分离。该方法计算成本低,但是利用常规分辨率Radon变换分离的绕射波场中会存在大量残余反射波,国内外学者致力于提高Radon变换的分辨率来改善波场分离效果。

Radon变换(Radon Transform,RT)首先由奥地利著名数学家J.RADON在1917年提出。20世纪70年代,美国斯坦福大学以CLAERBOUT为代表的地球物理小组首次将Radon变换应用于地震勘探中,并对线性Radon变换(τ-p变换)进行了重点研究。THORSON等[31]提出的随机反演方法提高了时间域Radon变换的分辨率;HAMPSON[32]提出频率域最小平方抛物线Radon变换,该方法因其高效的计算效率迅速成为当时工业界的标准;SCALES等[33]将迭代重加权最小二乘算法与预条件共轭梯度法相结合,有效求解了线性方程组中的大型稀疏矩阵;SACCHI等[34]提出高精度Radon变换,即以有效信号的稀疏性作为约束条件利用迭代高分辨率算法求解;CARY[35]指出时间域Radon变换结果具有更强的稀疏性而频率域Radon变换能更准确地刻画有效信号的频率特征;HERRMANN[36]提出去假频高分辨率抛物线Radon变换(DHR),该方法利用低频数据的无假频模型解构造加权算子来压制高频数据Radon空间的假频能量,获得了较高的信噪比和分辨率。

在国内,同样有许多学者致力于Radon变换研究。吴律[37]全面分析讨论了τ-p变换的原理和应用;牛滨华等[38]首次提出多项式Radon变换。许多国内学者[39-44]在提高Radon变换分辨率方面做了大量系统深入的研究工作。

为了解决Radon变换能量团泄漏,绕射波场分离不彻底的问题,本文在时间-频率混合域,引入预条件算子并构建了新的时变稀疏模型权,在倾角域CIGs上发展了一种基于预条件共轭梯度(PCG)算法的混合域高分辨率Radon变换绕射波分离成像方法。文中对比分析了本文方法与常规分辨率Radon变换、频率域高分辨率Radon变换在绕射波分离中的应用效果。

1 方法原理

1.1 倾角域共成像点道集波场特征分析

常速介质叠后偏移反射面的深度z(x)可以表示为(图1)[30]:

z(x)=z0+xtanα0

(1)

式中:z0表示倾斜反射界面与z轴交点的深度;α0为反射界面的倾角。在地表y处的自激自收地震响应可表示为:

(2)

式中:v表示介质速度,t为旅行时间。

图1 零偏移距反射示意

由模型坐标{x,z}和数据坐标{y,t}之间的映射关系可以得到偏移公式:

(3)

(4)

式中:vm是偏移速度;α为偏移倾角。将公式(2)代入公式(3)和公式(4),并消去y得到倾斜反射界面在倾角域的成像表达式:

(5)

对(5)式求偏导可以得到:

(6)

从公式(6)可以看出,当偏移速度准确时,即vm=v且α=α0时,(6)式的偏导数值为0,即α0为(5)式的一个极小值点。因此,对于准确的偏移速度,在倾角域CIGs上,反射响应是一个具有稳相顶点的“笑脸”状曲线,且稳相顶点的位置与地层倾角相对应。

类似地,可以得到地下绕射点(x0,z0)处的倾角域响应表达式:

(7)

zα(α)=z0

(8)

由(8)式可以看出绕射响应为一水平直线。

下面用3个理论模型来直观阐明倾角域道集曲线的特点,并进一步分析反射波和绕射波的差异。图2a 至图2c显示了3个理论模型,均包含一个反射界面和一个绕射体,绕射体位于地下1km处。图2a中水平反射界面位于地下0.5km处,图2b和图2c中反射界面分别向左倾斜和向右倾斜。我们从理论模型0.6km位置处(绕射点正上方)提取倾角域CIGs(图2d至图2f),图2d至图2f中虚线表示反射界面(水平和30°倾斜反射界面);实线表示绕射点处的倾角域响应曲线。

从图2d至图2f可以看出,对于反射界面来说,不论界面是否倾斜,在倾角域道集中响应曲线总是保持“微笑”的形状(如图中虚线所示),且曲线的稳相顶点刚好指示界面的倾角;而在绕射点位置处的绕射响应为水平直线(如图中实线所示)。

综上可知,在倾角域CIGs中,反射与绕射响应曲线存在明显的几何形态差异,其中反射响应为具有稳相顶点的凹形曲线,绕射响应为水平直线。我们可以利用两者之间形状差异将绕射波同相轴单独分离。下面以此原理我们采用抛物线Radon变换来分离倾角域CIGs中反射波和绕射波。

图2 理论模型及偏移速度正确时的倾角域共成像点道集示意a,b,c 包含一个反射界面和一个绕射点的理论模型; d a 在0.6km处的倾角域CIGs; e b 在0.6km处的倾角域CIGs; f c 在0.6km处的倾角域CIGs

1.2 Radon变换原理

二维Radon变换的定义式为:

(9)

式中:d(t,x)为时间域地震数据;m(τ,q)为Radon域数据;t,τ分别为时空域的时间和Radon域的截距时间;q表示曲率参数;x为炮检距。将(9)式进行离散傅里叶变换,转换到频率域,可得:

(10)

式中:D(ω,x)和M(ω,q)分别为d(t,x)和m(τ,q)的傅里叶变换结果;ω为频率;n为Radon域的曲率个数。由(10)式可知,对于某一频率ω,Radon变换可以表示为一个线性方程组的求解问题:

D=LM

(11)

式中:L是Radon变换算子,由地震数据道的相对空间位置x0,x1,…,xm-1和Radon变换参数q0,q1,…,qn-1来确定,可表示为:

(12)

(13)

(14)

常规Radon变换采用最小二乘反演计算,定义目标函数:

(15)

对上述目标函数求极小值,可以得到M的常规最小二乘解,即:

M=(LTL+βI)-1LTD

(16)

1.3 基于预条件共轭梯度(PCG)迭代的混合域高分辨率Radon变换

为了提高Radon变换的分辨率,在模型中引入稀疏约束条件,将目标函数改写为:

(17)

式中:R(M)表示稀疏约束项。考虑到模型空间的噪声对Radon变换分辨率的影响比实际地震数据中的非线性噪声影响更为严重,因此仅引入模型加权矩阵Wm来提高模型的拟合度,如分辨率和平滑度。求解方程

(18)

(19)

(20)

令A=F-1LTLF,b=F-1LTFD,则方程(20)可以写为:

(21)

可得残差量表达式为:

(22)

常规Radon变换采用共轭梯度(CG)算法求解,高分辨率Radon变换采用PCG反演方法,本文基于PCG算法的高分辨率Radon变换处理步骤如下。

1) 由方程(11)计算得到频率域常规分辨率解,经傅里叶变换由频率域变换到时间域,作为PCG高分辨率迭代的初始解m0。

3)内循环迭代。用PCG法[45]迭代求解方程(20),其中关键求解步骤如下。

② 求解迭代步长。αj:=(rj,zj)/(Apj,pj),其中,j=0,1,…,nmax为内部迭代次数,nmax表示设置的最大迭代次数。

③ 计算新的迭代解。mj+1:=mj+αjpj。

④ 计算更新残差。rj+1:=rj-αjApj。

⑥ 求解共轭方向的搜索步长。βj:=(rj+1,zj+1)/(rj,zj)。

⑦ 计算下一次迭代的共轭向量。pj+1:=zj+1+βjpj。

当多次迭代后达到事先设置好的误差极限时,即可得到(τ,q)的解。

4) 继续进行外部迭代,求取满意的高分辨率解。利用步骤3)中的(τ,q)解,重新计算稀疏约束矩阵,重复步骤2)和步骤3),继续迭代。外部迭代的主要作用是更新稀疏权,增加(τ,q)解的稀疏性。

5) 进行Radon反变换,完成反射波切除,得到绕射波分离结果。

2 数据测试

应用的技术流程主要包括:对倾角域CIGs进行Radon变换,在Radon域将开口向上的反射波聚焦于曲率小于0的位置,而绕射波的曲率在零值附近。在进行反射波切除时,选择的Radon域切除的q值在0值附近,以减少反射波能量的干扰并尽可能完整地保留绕射波信号。最后采用反Radon变换将分离出的绕射波同相轴恢复为时空域数据,即从全波场中分离出绕射波。

2.1 模拟数据测试

在模拟的倾角域CIGs中,包含两个反射波和一个绕射波响应曲线,反射波顶点分别位于1000m、1500m,绕射波则处于地下2000m,如图3所示。我们对比分析了最小二乘Radon变换方法、频率域高分辨率Radon变换方法和本文方法在Radon域中对绕射波的分离效果。采用最小二乘Radon变换对倾角域CIGs进行绕射波分离,其Radon变换结果以及绕射波分离后的结果如图4所示,基于该方法所得到的Radon域结果具有明显的剪刀状发散,这些有限孔径产生的端点假象严重影响了Radon变换的信噪比,进而导致在绕射波分离结果中存在大量残余的反射能量(图4b中红色箭头所示)。

图3 模拟数据的CIGs记录

采用频率域高分辨率Radon变换对模拟CIGs记录进行绕射波分离,其结果如图5所示。从图5可以看出,频率域高分辨率Radon变换结果相较于最小二乘Radon变换结果,剪刀状发散能量减弱,分辨率有了明显提高,但在绕射波上方的假象仍然存在。我们将内循环迭代次数分别选取为10次、20次、40次,谱白化因子分别选取为0.01、0.05、0.10,该假象依然存在。该方法得到的绕射波分离后的结果中同样存在反射波的残余。

图4 最小二乘Radon变换对绕射波的分离结果a Radon域结果; b 绕射波分离后的结果

图5 频率域高分辨率Radon变换对绕射波的分离结果a Radon域结果; b 绕射波分离后的结果

采用本文方法对该模拟数据进行绕射波分离,结果如图6所示。从图6可以看出,基于混合域高分辨率Radon变换的子波幅度无损失且分辨率很高,绕射波分离后的结果也优于其它两种方法。

绕射波场分离效果的好坏直接影响绕射波成像效果,对比3种方法的分离结果可见,基于混合域高分辨率Radon变换的绕射目标成像方法相较于其它两种方法得到的结果分辨率更高,子波保持得更好,这将有利于后续小尺度地质异常体的准确识别和定位。

2.2 实际资料测试

选取中国西部M地区实际地震资料进行测试,该地区溶蚀孔洞和裂缝储层十分发育,利用传统的成像处理流程难以对储层内的多种绕射目标体如速度异常体、河道、断层、生物礁等进行精确成像。首先利用本文方法对某一CDP点下Kirchhoff叠前深度偏移(PSDM)中产生的倾角域CIGs进行反射波和绕射波分离。

图6 混合域高分辨率Radon变换对绕射波的分离结果a Radon域结果; b 绕射波分离后的结果

图7a和图7b分别为实际地震数据的倾角域CIGs和采用本文方法分离后的绕射波数据。从图7可以看出,原始道集中的强反射能量被明显压制,由于大部分反射能量被去除,淹没在全波场中的弱绕射能量凸显出来,位于地下约5.5km处的绕射波同相轴能量明显加强,如图中红色箭头所指部位。

进一步将本文方法应用于该地区某条测线的倾角域CIGs,其波场分离前的PSDM剖面如图8a所示,由于强反射能量的干扰,绕射点、断层点的位置难以识别。采用本文方法分离后的绕射波成像剖面如图8b所示,岩丘边界刻画更加清晰,并且强反射同相轴被滤除,淹没的绕射点、断层点位置凸显出来,如图中红色箭头所指之处。图8a和图8b中蓝色方框部分的局部放大结果如图8c和图8d所示。从图8d中能明显区分出绕射点,如红色椭圆标示,地下小型绕射体也被凸显出来。因此,经过混合域高分辨率Radon变换滤波后得到的绕射波成像剖面可以联合反射波成像剖面进行准确、高确定性的地震解释。