编码与解码框架下的局部平面波域浅水多次波压制方法

徐 鹏,刘亚辉,王华忠

(1.波现象与智能反演成像研究组(WPI),同济大学海洋与地球科学学院,上海200092;2.中国石油大学(华东)地球科学与技术学院,山东青岛266580)

随着地震采集技术的不断发展,“两宽一高”地震勘探日渐常态化,但是地震数据的品质与成像结果的分辨率仍然受到多种因素的制约,多次波是其中重要因素之一。当前勘探地震学的核心问题是地震波反演成像,其目标在于利用低波数的光滑速度场估计高波数散射体。现有成像方法多数假设观测数据中的反射波仅在地下经历一次反射,而地震波在真实介质中的传播通常包含多个波阻抗界面的反射,即包含多次波,常规针对一次波的成像结果会含有假象,严重影响后续的数据分析与解释。

对于多次波的处理可以分为“压制”和“利用”两类。“压制”类是将数据中的多次波视为噪声去除后,利用常规的一次波成像方法处理;“利用”类则将多次波也视为有效信号,使用更复杂的成像方法对一次波、多次波进行联合成像。后者能更多地挖掘出数据中的地下介质信息,是更理想的成像方式。例如刘伊克等[1]对南海深水区域含多次波的数据进行偏移成像,成功弥补了目标区域的照明不足;王永强[2]利用最小二乘逆时偏移(LS-RTM)对一次波-多次波进行联合成像,验证了一次波-多次波联合成像方法的可行性。

但是,就目前研究现状而言,多次波成像方法适用范围有限,因为它对偏移速度场的要求过于严格,同时为了消除不同阶次的反射波之间的错误相关值,引入最小二乘(LS)迭代流程,使得总体计算量增加。因此,先去除多次波,再利用成熟的常规算法进行成像是目前更为合理的选择。

勘探地震技术发展至今,已有多种多次波压制方法在实际生产中得到应用。根据其基本原理,可将这些方法分为两类:第一类是滤波方法,基于多次波与一次波在道集中的形态特性差异,将多次波剔除;第二类是褶积方法,基于多次波的物理成因,由波动理论建立一次波与多次波的高维褶积模型,预测出多次波的位置、形态,实现对多次波的精确压制。第一类方法主要利用一次波与多次波(在变换域)的可分离性压制多次波,具有成本低、计算量小等特点,实际应用广泛,代表方法有最优叠加、特征值谱滤波、f-k域滤波等。然而,由于多次波与一次波并非完全可分,因此压制效果差强人意。第二类方法采用“预测-减去”的策略,预测阶段主要考虑多次波的运动学信息,减去阶段则主要用动力学信息进行匹配。代表方法有经典预测反褶积、纯数据驱动的自由表面多次波消除(SRME)[3-6]、基于模型的浅水多次波压制(简称SWD)[7]等。该方法将多次波压制分解为多次波模型预测与自适应相减两步,降低了实现难度,但是引入了额外的中间数据,降低了自动化程度[8-9]。

浅海环境为多次波压制问题引入了额外的难点。浅海环境主要的多次波是水体相关多次波,由于海底-海面路径短,且海面、海底都是强反射界面,因此该类多次波的特点是周期短、能量强,对有效信号干扰极强[10-11]。然而,由于现有的观测系统大多是针对中深层目标体设计的,用于浅海底的反射波观测时采样间隔偏大,同时更易造成近偏移距数据的缺失,故常规的多次波压制方法不能直接用于浅水数据的多次波压制。例如,经典的SRME作为一种数据驱动的高维褶积类方法,对数据采样间距、完备性有较高要求,因而无法在浅水数据中取得令人满意的效果[12-15]。SWD方法引入了海底形态模型,对海底反射的刻画更准确,但其本质仍是Rayleigh积分的数值实现,故对积分孔径、积分网格都有较高要求,在近偏移距数据缺失时仍有局限性。因此,需要发展针对浅海环境的多次波压制技术。

本文针对浅海环境中多次波问题的特点,引入信号领域的“编码-解码”框架,设计一种针对浅水多次波的“预测-减去”类压制方法,借助稀疏τ-p变换克服不规则空间采样问题,再利用射线理论高效生成多次波模型,最终在局部平面波域完成观测数据与多次波模型的匹配相减。

1 水体多次波问题描述

水体中存在海水-空气、海水-海底两个强反射界面,因此海上地震数据中多次波往往更为发育。海上勘探中多次波种类主要有:鬼波、水体相关多次波、其它自由表面相关多次波以及层间多次波。前三者成因与海上观测系统紧密相关,层间多次波基本与陆上勘探相同。

鬼波可分为源鬼波、检鬼波和源检鬼波。震源激发后,地震波上行至水面发生反射继而下行,稍迟于震源直接激发的下行波传播,即源鬼波。地震波与地下介质相互作用后上行至海面反射一次后,再传至检波点,即检鬼波。若同时符合源鬼波、检鬼波的路径特点,则称为源检鬼波。

水体相关多次波是自由表面相关多次波的一种特例,以海底向上反射开始或/和结束,可分为源端、检端和源检两端三种。源端水体相关多次波是指地震波自震源激发后,先在海底-海面间震荡一次,再自海面向深层介质传播的波现象。检端水体相关多次波是指地震波穿过地层介质后,最终被检波器接收之前,在海面-海底间反射一次所形成的波现象。同时满足源端水体相关多次波、检端水体相关多次波路径特点的波现象被称为源检两端水体相关多次波。另外,仅在水层中震荡传播的波现象被称为纯水体鸣震[16-17]。

除鬼波、水体相关多次波之外的自由表面多次波统称为其它自由表面相关多次波。此类多次波含义较广,例如在海底之下波阻抗界面与水面之间的多次反射、多次折射等。

在海上勘探中(尤其是浅海探区),首先必须面对、同时也是干扰最强的多次波是水体相关多次波。水体相关多次波的存在,掩盖了深层反射信号,影响了微断层、盐丘下部构造的成像[12,18]。

2 正问题:用编码观点看待多次波形成正过程

从各种多次波成因来看,多次波实质上是一次波的延拓结果。因此可从多个角度对多次波进行表述,如将多次波表述为波动方程的积分解、镜像虚震源激发的波场等。本文从信号分析角度出发,将观测数据视为输入信号,将地下介质整体视为线性滤波器(即“地层滤波系统”),多次波视为输出信号。即将多次波视为一次波的编码结果,而编码过程是一次波在地下介质中震荡传播的过程,或者说取“地层滤波系统”的脉冲响应作为编码算子。将多次波形成的正过程表述为编码过程,则多次波的消除过程可理解为解码过程。引入信号分析领域的“编码-解码”框架使得输入信号、编码算子、输出信号等概念更加直观,简化整个多次波形成过程的描述,有利于多次波压制算法的设计与实现。

图1a为真实地下介质(含海水和空气)示意图,为方便讨论,引入参考介质,如图1b所示,即去除海面的自由边界并在原海面及以上空间填充海水的假想介质。基于图1a介质的地震波传播情况如图1c所示,可见由于海水自由表面的存在,一次反射波又被向下反射,继续与地下介质发生作用,进而产生一阶多次波;一阶多次波继续产生二阶多次波,直至无穷阶多次波。由于反射系数绝对值小于1,且地震波传播会发生衰减,高阶多次波可以忽略不计。

在参考介质中,震源直接激发不会产生自由表面多次波。真实观测数据中的自由表面多次波可看作将真实观测数据作为边值条件加载在深度z=0平面上所引起的波场(图1d)。可用经典的波动方程积分解解释:以所有上行波在海面的反射R0U(x′;xs,ω)为边值条件,以海面+半径无穷大的下半球为封闭面,最终在检波点引起的波场。即:

(1)

式中:M(xr;xs,ω)表示自由表面相关多次波;U(x′;xs,ω)表示观测数据(含一次波、多次波);Gref(x′;xr,ω)表示参考介质中的格林函数;ω表示频率;xs,xr,x′分别表示震源点、检波点以及积分微元的坐标;R0表示海面反射系数;积分曲面S取为海平面;n表示积分微元处的单位法向量。公式(1)即水体相关多次波在“编码-解码”框架下的表达式。

等价地,可以用向量形式描述上述过程:

式中,输入信号为所有上行波,编码算子为海水层的脉冲响应,输出信号为检波点观测到的波场。应该指出,引入编码算子A是为了简洁表达(1)式的线性关系,相当于离散化表示的积分核。后续计算中不需要将该算子显式表达为矩阵,仅需要实现A的矩阵向量乘。

地震数据中的一次波及各阶多次波表达式可以拆分,得到更直观的编码预测关系:以一次波在海面的反射作为边值条件,在参考介质中传播得到的波场即为一阶检端多次波;以一阶检端多次波在海面的反射作为边值条件,在参考介质中传播得到的波场即为二阶检端多次波;以此类推,以k阶检端多次波在海面的反射作为边值条件,在参考介质中传播得到的波场即为k+1阶检端多次波。

图1 地下介质及地震波传播a 真实地下介质; b 参考介质; c 地震波在真实地下介质中的传播; d 在参考介质中借助面源描述地震波传播

(4)

式中,Mi表示第i阶自由表面相关多次波,P表示一次反射波。总之,算子A的作用是海水层对海面处上行波的编码“升阶”。

分别累加各式等号两端,得:

(5)

式中,Mtotal为观测到的总体多次波(各阶多次波之和),右侧为观测数据。

同样地,若将各阶多次波都写为关于一次波P的函数:

(6)

可得总体多次波与一次波的关系:

(7)

其中,总体多次波不是一个直接观测的量。为了建立一次波P与直接观测结果Uobs之间的关系,在(7)式左右两端都加上P,因此得到描述一次波与含多次波的实际观测数据之间的线性编码关系:

(8)

其中,I表示单位矩阵。

3 反问题:用解码的观点看待多次波压制

(9)

其中,L为迭代次数。

令:

(10)

将公式(10)代入公式(9),得:

(11)

(12)

公式(12)即为水体相关多次波的求解公式,它与SRME的“预测+减去”公式有类似之处。相对SRME使用观测数据来代替(1)式的Rayleigh积分,本文使用了特征波域的射线束传播算子更为精确地实现编码算子A,完成了(12)式中的A(ω)Uobs(ω)计算,即总体多次波的预测。

公式(12)描述了理想情况下的多次波压制或一次波估计的数学模型。但事实上,由于对水体介质各属性(如海面和海底的反射系数分布、水体速度模型、海底起伏形态)不能进行精确描述,实际数值实现的A算子都是近似的。

多次波压制过程具体可描述为:

1) 利用公式(13)计算或预测多次波模型。

(13)

2) 利用公式(14)估计一个整形滤波器f,使得输出的一次波能量最小。

(14)

3) 采用公式(15)输出一次波。

(15)

借鉴SRME方法,公式(12)的减法是自适应减法,且在局部平面波域进行。相对于空间域减法,局部平面波域自适应减法能更好地区分多次波与一次波的斜率差别,也更能容忍两者的交叉,使得输出的去噪结果质量更好。

4 在平面波域实现多次波的编码预测与一次波的解码

公式(1)是将海面每个微元作为参考介质中激发的偶极震源,积分计算所有微元在检波点处所产生的波场。作为波动方程的积分解,理想情况下公式(1)在数学上是精确成立的,并有明确的物理意义。然而在实际应用中,有限的积分孔径、不规则的面元分布、近偏移距数据的缺失等问题,使得公式(1)的精确性下降。编码算子A是公式(1)中积分所描述的线性关系的紧凑表达。在局部平面波域实现编码算子A可以缓解上述问题。在每个局部区域,作为输入信号的入射波上行波场都可分解为不同方向的平面波[19-20]。局部区域激发的二次波场的传播过程可以看作一组不同方向的平面波传播过程。该方法在地震成像领域应用广泛,如冯波[21]通过局部平面波分解解决波形反演中的“周期跳跃”问题;刘少勇[22]通过局部平面波分解解决常规偏移中的画弧问题;孙维蔷[23]将局部平面波分解用于鬼波与自由表面多次波的统一压制。

本文将局部平面波分解用于描述多次波的形成过程。单个平面波分量的传播过程如图2a所示,相比传统预测方式(图2b),该方法考虑了各个平面波分量的方向,同时简化了地震波在水体中的传播过程。局部平面波域水体相关多次波预测方法的实现步骤如下:

1) 将计算区域的海面分解为若干局部窗。

2) 通过稀疏求解局部窗τ-p系数,将局部窗内的入射波上行波场分解为不同方向的平面波。

3) 从局部窗中心向检波面进行反射射线追踪(从窗中心出发,在海底反射,终点在检波面,见图2a),并记录各条射线出射方向、走时(T)。

4) 将局部平面波当作有方向、有宽度的射线来传播,各个τ-p系数按照与其方向对应的射线的走时计算延迟,然后“放置”到射线路径终点对应小邻域内检波点的数据上(需考虑相对于各检波点的时差)。

5) 对所有局部窗重复步骤2)～4),完成所有计算。

预测过程主要利用运动学信息构建了多次波模型,其波形、振幅等动力学信息与真实数据仍有差异,因此精细的“匹配相减”是不可或缺的措施。得益于本文方法中的“平面波分解”步骤,可以将“匹配相减”过程放在局部平面波域,从而能额外考虑方向信息,使得匹配过程更为精细。

图2 观测数据作为二次源对于多次波的不同贡献方式a 本文方法中的局部平面波源方式; b 积分解形式中的偶极点源方式

5 实际应用

为了验证本文方法的有效性及实用性,设计、实施了浅海环境的OBC数据实验和拖缆数据实验。

5.1 浅海OBC数据

浅海OBC数据共含有8条检波线,每条检波线有240个检波器。炮点、检波点间距为25m,炮线间距250m,检波线间距400m,观测系统如图3所示。采集环境平均水深约37m,水体相关多次波周期约为50ms(该数据中单个震相的波形时长约40ms),严重干扰了有效信号。

图3 OBC数据观测系统

采用以下步骤对多次波进行了压制:

1) 通过“近道自相关”提取工区网格内的水深,并插值为光滑曲面作为海底模型。该海底深度模型是深度关于坐标的二阶光滑函数。

2) 以单炮道集为独立单位进行多次波预测。将单炮数据按检波点坐标划分为若干局部空间窗口,并利用稀疏τ-p变换进行平面波分解。

3) 在每个空间窗内,将所有平面波分量“传播”到检波阵列中。即将每个平面波分量当作一束有宽度的射线,按照其起飞角(由参数p计算),沿“海底-海面”反射路径,贡献到检波阵列(图2a),最终得到该单炮道集的多次波预测结果。

4) 进行“局部平面波域匹配相减”,减去多次波,输出结果即为一次波。

考虑各个单炮以及各个空间窗之间的独立性,本文在各个单炮道集间采用MPI多进程并行计算,每炮内各空间窗之间采用OpenMP多线程并行的计算策略,从而提高计算效率。

本文方法中的稀疏估计平面波系数方法可以利用有限样点还原出局部波前面,从而缓解采样不足问题,因此比常规方法具有明显优势。其中局部平面波传播在三维坐标系下进行,能够自动适应实际三维OBC观测系统。OBC数据实验结果表明,本文方法能够准确地预测并消除水体相关多次波。图4和图5 对比了多次波压制处理前后OBC炮集和OBC数据体,图6为原始OBC数据及其平面波系数,图7为多次波压制处理后OBC数据及其平面波系数,可以看出多次波已被成功分离。图8对比了OBC炮集多次波压制处理前后自相关结果,可见数据自相关中的周期性成分已被成功消除。

图5 原始OBC数据体(a)与多次波压制处理后OBC数据体(b)对比

图6 原始OBC数据(a)及其局部平面波系数(b)

图7 多次波压制处理后OBC数据(a)及其局部平面波系数(b)

图8 OBC炮集多次波压制处理前(a)、后(b)自相关结果对比

5.2 浅海拖缆数据

浅海拖缆数据共1149炮,炮间距25m,每炮4条拖缆,检波器间距12.5m,拖缆间距100m。采集环境水深90m,多次波周期约120ms。该数据中单个震相的波形时长约40ms,水体相关多次波的存在使得有效数据同相轴下方多出若干“虚影”。相比OBC数据,拖缆数据的观测系统更具挑战性:近偏移距数据缺失、拖缆发生漂移甚至弯曲(图9)。得益于稀疏τ-p反演,本文方法在面对非规则网格时,能根据采样点真实坐标将数据分解为不同方向的平面波,解决了观测系统带来的问题。

拖缆数据多次波压制的步骤与OBC数据实验相同。图10对比了拖缆数据多次波压制处理前后的叠加剖面,可以看到多次波的影响已被成功消除。图11对比了拖缆数据多次波压制处理前后自相关结果,可以看出零延迟之外的旁瓣已被消除,这意味着数据中的周期性成分已被成功压制。

图9 拖缆数据中炮检点分布

图10 拖缆数据多次波压制处理前(a)、后(b)叠加剖面对比

图11 拖缆数据多次波压制处理前(a)、后(b)自相关结果对比

6 讨论与结论

本文分析了浅海环境中多次波问题的特点,引入信号分析领域的“编码-解码”框架,提出了针对性的多次波压制方法。将多次波的产生与压制用一个线性系统及其逆进行描述,即将上行波在水层震荡产生多次波的过程视为编码过程,将多次波压制过程视为解码过程,以精简的方式从线性系统求解角度给出了压制算法。该方法将水体相关多次波的预测转移到局部平面波域,利用稀疏τ-p反演缓解了采样不规则、近偏移距数据缺失造成的影响。利用浅海OBC数据和拖缆数据验证了方法的有效性和高效性,展示了方法在浅海环境地震数据处理中的应用前景。

本文方法相对传统方法具有如下优势:①稀疏求解τ-p系数,能在非均匀网格下得到各个窗口的平面波分解,克服了非均匀采样问题,缓解了近道缺失的影响;②所需的海底模型为光滑曲面,不需要离散建模,因此能应对崎岖海底,且不会引入模型的锯齿和波场的数值频散现象;③传播过程化简为从窗中心到海面的反射射线追踪更高效,由于射线追踪计算次数等于局部窗口数目,远小于检波点数目,故总体效率也高于MWD(model-based water-layer demultiple)算法中从检波点到检波点的射线追踪;④各个单炮道集以及各个局部窗的计算过程相互独立,具有良好的可并行性;⑤将匹配相减放在局部平面波域进行,能够区分一次波、多次波的不同视速度,提高多次波压制精度。

但是,本文方法在多次波预测过程中将每个平面波分量当作一束有宽度的射线(即胖射线)传播,对于射线宽度的选择尚缺理论依据,需要进行手动调试。射线过宽时,所构建的多次波模型的频带将偏低;射线过窄时,构建的多次波模型会出现一些空缺点。可能的解决方法是根据“菲涅尔半径”确定射线宽度,需要开展进一步的探索研究。

展望未来多次波相关技术的发展,“压制”与“利用”应当并重。勘探地震学的重要方法如层析、FWI、成像等,目前仍依赖于一次波假设。在实际应用中,多次波压制的效率与效果,将与最终处理、解释结果紧密相关。因此,发展有效的多次波压制技术具有重要意义。同时,多次波成像应予以重视,一种能够应对多次波的成像方法将大大简化地震数据的预处理,并能更多地挖掘出数据中的地下介质信息,最终呈现出更准确的成像结果。故对待多次波,应当立足“压制”,展望“利用”。

致谢:感谢中海油研究总院和湛江分公司对本文研究的支持。感谢中石油勘探开发研究院及西北分院、中国石油化工股份有限公司石油物探技术研究院和胜利油田分公司对波现象与智能反演成像研究组(WPI)研究工作的资助与支持。