2. 中国科学院地质与地球物理研究所, 北京 100029;
3. 洛恩能源科技(北京)有限公司, 北京 100012
2. Institute of Geology and Geophysics, Chinese Academy of Sciences, Beijing 100029, China;
3. Kerogen Energy Services, Beijing 100012, China
在海洋地震勘探中,由于海表面是一个强阻抗界面,海上地震资料中往往存在强烈的多次波.常规地震成像技术只对一次波进行成像,而多次波的存在产生严重假象.近年来发展的多次波成像技术虽可利用多次波进行成像(Liu et al., 2011, 2015, 2016),但多次波成像技术采用了不同的成像条件,因此需在成像前分离一次波和多次波.综上所述,无论对常规一次波成像技术还是先进的多次波成像方法,多次波去除或一次波和多次波分离的效果严重影响成像的精度,因此对该技术的研究具有重要的应用价值.
多次波去除方法分为滤波类方法和基于波动方程的反演类方法.滤波类方法直接通过信号分析,利用多次波和一次波的差异特性直接滤除多次波(Robinson,1957;Taner,1980),对于复杂实际地震资料,这类方法的利用效果有限.基于波动方程的反演方法(Berkhout,1982;Verschuur,1991;Weglein et al., 1997;van Groenestijn and Verschuur, 2009;Ypma and Verschuur, 2013)利用地震波场的动力学特征预测多次波.由于其理论上的先进性,该类方法是目前多次波去除方法研究的重点和主流.
基于波动方程的反演类方法通常的实现方式是先预测再匹配相减.即先基于波动方程预测多次波,该预测多次波的过程主要考虑多次波的走时信息,因此其预测结果在波形、振幅和相位上与实际多次波存在较大的差异,需通过最小二乘意义下的自适应匹配滤波后才还能获得只含一次波的地震记录.这类方法包括模型驱动的波场延拓多次波预测去除方法、数据驱动的自由表面相关多次波消除法(Surface-Related Multiple Elimination, SRME)(Verschuur,1991;Berkhout and Verschuur, 1997;李鹏,2007)和逆散射级数法(Inverse Scattering Serious method, ISS)(Weglein et al., 1997;金德刚等,2008)等.其中SRME方法利用反馈迭代模型(Berkhout,1982;Verschuur,1991)预测自由表面相关多次波,再进行自适应的匹配相减.由于其理论上具有先进性且先预测再相减的算法容易实现,SRME是目前学术界和工业界研究和应用的主流方法.但该方法仍存在局限性,包括:1)要求地震数据必须完整、规则,对于缺失数据必须进行数据重建、规则化等预处理(Herrmann and Hennenfent, 2008;白兰淑等,2014)等.2)当近道缺失时,SRME无法预测出浅水多次波.3)匹配滤波会导致与多次波发生重叠的一次波产生畸变或损失.
针对经典SRME方法存在的问题,直接分离一次波和多次波的EPSI方法得到了研究.该类方法同样利用地震波场的动力学特征建立目标函数并通过求解其最小化问题获得一次波和多次波,包括稀疏反演一次波估计方法(Estimation of Primaries by Sparsity Inversion,EPSI)(van Groenestijn and Verschuur, 2009;van Groenestijn,2010;宋家文等,2014)和广义EPSI(Generalized Estimation of Primaries by Sparsity Inversion)方法(Ypma and Verschuur, 2013).EPSI同样以反馈迭代模型作为理论基础,但是与SRME不同,该方法建立使一次波和多次波估计总残差最小、稀疏约束一次波脉冲响应的目标函数,通过最速下降算法获得一次波脉冲响应和地震子波.之所以称其为稀疏反演是因为在反演过程中需对一次波脉冲响应进行稀疏约束,使EPSI求解的欠定问题在此约束下得到可靠的结果.一次波可通过一次波脉冲响应和地震子波的褶积直接获得,或从原始数据中将多次波去除获得.多次波通过一次波脉冲响应、自由界面反射系数、总地震波场多维褶积获得.因为EPSI在反演过程中无匹配相减过程,因此相比SRME而言可以更好地保留一次波,减少一次波的损失.另外,EPSI方法是基于全波形数据的反演,因此可以利用多次波信息重建出缺失道上的一次波.但EPSI方法的计算量相比SRME要大几倍甚至更多.Lin等2013年提出了基于稀疏反演的稳健一次波估计方法(Robust Estimation of Primaries by Sparsity Inversion,REPSI)(Lin and Herrmann, 2013).该方法采用L1范数最小化策略对一次波脉冲响应进行稀疏约束,数据测试结果中一次波脉冲响应清晰、干净,反演效果好,但仍然未能克服计算量大的问题.
文章首先分析SRME和EPSI的理论和算法实现的关键技术,通过一个水平层状模型和Pluto复杂盐丘模型对比分析两种方法在不同条件下的应用效果.最后针对一个海上实际资料进行了试处理,深入探讨了两种方法的优缺点.
2 理论基础 2.1 反馈迭代正演模型反馈迭代模型通过基于数据自身的计算获得地震波场正演结果,其频率域正演模型表达式如下:
(1) |
其中,P-是接收点记录的全波场数据,P0是表面相关一次反射波,M0是表面相关多次波,A(ω)I是匹配算子,ω是单一频率,下标“0”代表自由表面,“-”代表地震波的上行方向.
2.2 表面相关多次波迭代去除方法(Iterative Surface-Related Multiple Elimination)表面相关多次波迭代去除方法(Iterative Surface-Related Multiple Elimination,以下简称SRME)(Verschuur,1991;Berkhout and Verschuur, 1997)从反馈迭代正演模型(Berkhout, 1982; Verschuur, 1991)出发预测多次波.(1)式中,P0P为多次波预测结果,与真实多次波波形在振幅、相位上存在不一致,因此需要A(ω)算子对其进行波形匹配;M0=A(ω)P0P是波形匹配后的多次波预测结果.
SRME方法基于一次波能量最小假设,其反演目标函数可表示为如下形式:
(2) |
其中,is、ir分别是激发点、接收点序号,NS、NR分别是总炮数、总道数.由于P0为待求未知数,(2)式无法直接求解.SRME方法实际通过迭代求解P0,其迭代公式为
(3) |
其中P0(i)代表第i次迭代后的一次波结果,Niter为总迭代次数.(3)式中有两个未知量,P0和A(ω)I.每次迭代中,首先在频率域求解P0(i)P,再通过在时间域求解(2)式的最小二乘问题求出地表算子A(i+1)(ω)I,再从P中减去A(i + 1)(ω)P0(i)P获得一次波P0(i+1).匹配波形的滤波算子有多种方式,如单道匹配、多道匹配、均衡拟多道匹配等.
SRME方法在工业界已得到广泛应用,但实际应用仍存在诸多局限.首先,SRME无法直接处理含缺失道的或不规则的地震数据,必须先对数据进行缺失道差值或数据规则化之后才能预测多次波,其预测效果依赖于插值的准确度.其次,匹配算子的解存在不确定性导致匹配相减的结果不佳,尤其是当多次波与一次波发生重叠时,相减结果的一次波存在能量损失和波形畸变.
2.3 稀疏反演一次波估计方法(Estimation of Primaries by Sparsity Inversion,EPSI)Groenextijn和Verscchuur在2009年提出了基于稀疏反演的一次波估计方法(Estimation of Primaries by Sparsity Inversion,EPSI)(van Groenestijn and Verschuur, 2009).与SRME一样,EPSI也基于反馈迭代正演模型,因此EPSI同样无需提供地下介质信息.与SRME不同之处在于EPSI方法实现过程不是采用先预测再相减的两步法,而是直接分离一次波和多次波.
EPSI方法对地震总波场的频率域表达形式为
(4) |
其中,X0为一次波脉冲响应,S+为地震子波,R-为自由表面反射系数,“+”代表地震波的下行方向(以下省略S+的上标“+”和R-中的“-”).一次波P0用X0与S的多维褶积表示,其中S=S(ω)I;第二项X0RP-为多次波,即X0、R、P-的多维褶积.
基于(4)式,EPSI的目标函数J写成如下形式:
(5) |
其中,‖·‖0为L0“范数”.(5)式中等号右侧第一项为总波场减去一次波、多次波的波场残差,第二项为一次波脉冲响应的L0 “范数”稀疏约束.式中含有两个未知数X0和S,EPSI通过对(5)中的目标函数关于X0求一阶偏导数得到最速下降法反演的校正公式,如下:
(6) |
其中ΔX0为X0的梯度,i代表迭代次数,iis、is、ir代表不同的炮、检点下标,(·)H符号代表矩阵的伴随矩阵,即共轭转置,共轭意味着时间相位发生反转,转置意味着炮检位置相互对调.(P--X0, iSi-X0, iRP-)为上一次迭代后的波场残差.该式的物理意义可从图 1中得到更直观的理解.图 1a中黑色实线代表需要求解的一次波脉冲响应(ΔX0)ω, ir, is;黑色虚线代表(P--X0S-X0RP-)ω, iis, is,由Sris出发在Reciis到达;红色实线代表(S+RP-)ω, ir, iis,从Reciis处出发Recir处到达.图 1b中,(S+RP-)ω, ir, iisH,从Srir出发在Reciis处到达,且走时相位发生了反转.图 1(a—b)表达的意思相同,即黑色实线通过将红色实线路径上的走时、波形从黑色虚线中去除获得.
可以注意到(6)式中的(·)H里包含两项,即Si和(RP-)H,两项具有不同的作用.(P--X0, iSi-X0, iRP-)ω, iis, is与(Si)ω, iis, irH的褶积代表残差项的去子波效应,该部分对应的一次波脉冲响应ΔX0的更新量从一次波自身获得;(P--X0, iSi-X0, iRP-)ω, iis, is与(RP-)ω, ir, iisH的褶积是沿波传播的反方向对残差项中地震数据进行反向延拓.从而,两者中的同阶地震波场反向延拓聚焦至炮点附近,而相差一阶次的地震反射信息则反向延拓至一次波位置,也就是该项对一次波脉冲响应更新量的贡献来源于多次波信息.EPSI方法同时从一次波和多次波信息中获得一次波信息,因此即使记录中含有缺失道,只要数据中含有较强能量的多次波,EPSI可从多次波信息中重建出一次波信息.
EPSI假设一次波脉冲响应ΔX0在时空域由一个个尖脉冲组成,具有稀疏性,且一次波能量比多次波或其他噪声强.基于此假设,在时间域对每个ΔX0共炮记录选取一个时空窗,在窗口内拾取每一道中振幅最大的一个或多个点作为稀疏约束下的一次波脉冲响应值ΔX0(包括走时和脉冲振幅值),初始选择时间窗要将聚焦波场和多次波排除在外.随着迭代进行,窗口逐渐扩大至整个数据.加窗拾取一次波脉冲响应是为了提高算法的稳定性.首先,相关波场中除了真正的一次波信息,还含有相关计算过程产生的一些假象,这些假象包括同阶地震波场相关所得到的炮点附近的聚焦波场以及不同地层对应的同相轴之间相关所得到的串扰波场.聚焦波场能量集中在炮点附近,能量往往又比一次波强,选取稀疏脉冲前需对其进行切除,或将其排除在时空窗范围之外.串扰波场的错误拾取则在后续的反演过程中可以得到很大改善.其二,初期迭代时,ΔX0中除了一次波脉冲响应还包含多次波信息.具体原因是因为(P--X0, iSi-X0, iRP-)(Si)H会使残差中的多次波仍然残留在其原始位置,而(P--X0, iSi-X0, iRP-)(RP-)H则会使残差中的高阶多次波通过相关计算,反向延拓至低阶的多次波位置.此时如果地下有强反射层,那么该反射层对应的多次波能量非常大,甚至比大多数一次波能量更强.一旦多次波脉冲响应被当作一次波脉冲响应拾取,这些错误拾取在后续迭代中无法消除.为了避免这种错误拾取,迭代初期选择较小时窗将强能量多次波排除在外.一至两次迭代后大部分强能量的多次波已经从波场中去除或者能量很弱,此后时间窗便可扩大至整个数据范围.
获得稀疏约束的一次波脉冲响应ΔX0之后,将其加到上一次迭代结果中更新一次波脉冲响应:
(7) |
其中,α是迭代步长,本文按照下式求取:
(8) |
求解得到X0, i+1之后,EPSI求解地震子波S.地震子波求取一般在时间域获得,即将X0, i+1与(P--X0, i+1RP-)在最小二乘意义下匹配获得.
通过迭代重复上述过程可获得最终的
多次波特点是其传播路径比一次波更长,多次波中包含更多的地下介质相关信息.在数据缺失情况下,SRME方法本身无法获得缺失道处的一次波,因此缺失道附近的多次波去除效果也不理想.这时只能先通过插值处理重构缺失数据,再用SRME算法预测多次波,然而这一预测结果依赖于插值方法的准确程度,对于连续缺失多道的情况,预测效果定会受较大影响.前面提到EPSI获取一次波的来源有两部分组成,即一部分来自数据中的一次波本身,还有一部分来自多次波中包含的一次波信息(因为一次波是多次波传播路径上的一部分),因此EPSI可以重建出缺失道上的一次波.图 2中,灰色区域表示近道缺失数据,近道缺失区域的一次波脉冲响应(ΔX0)ω, ir, is可以由虚线路径所代表的(P--X0S-X0R-P-)ω, iis, is和红色路径代表的(S+R-P-)ω, ir, iisH多维相关获得,其中Reciis处坐标iis=1, 2, …, NS.图 2描述的仅是近道缺失情况的一次波重建原理,EPSI重建原理也适用于地震道不规则缺失的地震数据.
SRME和EPSI建立的目标函数不同,获得一次波的算法也不同.两者的算法框架分别如下.
SRME算法框架
1) 求解P0(i-1)P-,其中P0(0)=P-;
2) 在时间域最小二乘匹配P0(i-1)P-与P-获得匹配算子A(ω)(i);
3) 从P中减去A(i)(ω)P0(i-1)P获得一次波P0(i);
4) 返回1)迭代多次.
EPSI算法框架
1) 在频率域按(6)式获得一次波脉冲响应更新量ΔX0,第一次迭代时X0=0,因此直接计算ΔX0=P-(RP-)H,其结果如图 3b所示(图 3a为原始地震记录);
2) 将1)中的结果变换至时间域,在一个时间窗内拾取每道最大振幅值作为稀疏约束后的一次波脉冲响应更新量ΔX0;
3) 根据(8)式确定迭代步长α;
4) 根据(7)式获得更新后的一次波脉冲响应X0, i+1.第1次和第10次迭代后的一次波脉冲响应分别见图 3c和图 3d;
5) 获得地震子波S.在时间域最小二乘匹配X0, i+1与P--X0, i+1RP-获得.根据已求取的X0, i+1和S,通过
6) 判断是否满足迭代结束条件(最大迭代次数或者残差小于给定阈值).若不满足,继续重复第(1)步-第(6)步;若满足,则迭代结束.
图 4(a—c)、(d—f)分别展示了第40、70、110炮的一次波脉冲响应和一次波估计结果,共30次迭代,可以发现表面相关多次波去除效果很好.图 4(d—f)的一次波结果中含有层间多次波,这部分多次波常规EPSI无法处理.
Pluto模型地震数据是用来检验多次波去除方法的典型数值模拟数据,其速度模型如图 5所示.该速度模型顶部有海水层,海底包括很多沉积层和多个盐丘.盐丘是高速异常体,因此该模型中除了海底强反射界面之外还含有盐丘顶底对应的强反射界面.图 6a、6d展示了Pluto数据的典型单炮记录(已做了直达波切除处理).图中两个能量最强的地震同相轴分别是海底界面和海底盐丘对应的一次反射波.由于海底面以及盐丘顶底界面的反射率大,因此在这两个界面发生反射的表面相关多次波能量也非常强(图 6a、6d中红色箭头所指处),其能量盖过了与其发生重叠的一次波.图 6b、6e展示了通过EPSI方法从图 6a、6d所示原始记录中得到的一次波估计结果.可以看到表面相关多次波得到了很好的压制,原本被强能量多次波覆盖的一次波同相轴也变得更加清晰,没有产生波形畸变.多次波结果展示在图 6c—6f中.图 7a展示了第360炮原始单炮记录,图中黑色箭头处多次波和一次波产生重叠.图 7b—7c分别展示了EPSI方法和SRME方法的处理结果,可以发现EPSI方法处理结果中一次波保留得更好.SRME方法先预测多次波再对其进行匹配相减造成一次波损失.EPSI方法基于(5)式直接获得一次波,不涉及多次波的匹配相减过程,因此大大减少了一次波能量损失.EPSI方法准确性依赖于对一次波脉冲响应的准确估计,如果随着迭代计算无法更正遗漏或者误拾的一次波脉冲响应将会导致一次波预测结果产生一定程度的误差.另外,EPSI在计算成本上相比SRME要大.采用8个并行进程处理540个Pluto共炮点数据的具体计算参数如表 1所示.EPSI方法采用了40次迭代,整个计算消耗了235.85 min,近SRME的4.5倍.EPSI的内存消耗也比SRME大,约为1.25倍.进一步发展快速稳定的一次波脉冲响应反演算法是提高EPSI的计算效率的有效途径.
前文中提到EPSI对一次波的估计由两部分来源获得,一部分来自一次波本身,另一部分来自多次波.因此,当数据不完整时EPSI可从多次波重建出缺失道上的一次波.为了展示EPSI方法对近道缺失数据的处理效果,我们用Pluto部分数据(179个共炮集记录的前580个时间采样点的数据)进行了测试.图 8a展示了原始完整单炮记录,图 8b展示了近道缺失的单炮记录,图 8c展示了EPSI得到的一次波脉冲响应估计结果,可以看到一次波脉冲响应连续、位置准确.图 8d展示了EPSI的一次波估计结果,通过将图 8c中所示的一次波脉冲响应与地震子波褶积获得.可以看到原本近道缺失的一次波得到了很好的重建.图 8e为SRME方法得到的结果,可以看到SRME无法重建出近道缺失的一次波.SRME方法常常要借助插值技术先对近道进行数据插值,而插值之后的多次波估计结果依赖于插值方法的准确度.
SRME通过一次波和波场的多维褶积运算预测多次波,因此原始数据中包含完整的一次波数据完整才能预测出完整的多次波.EPSI方法则通过一次波脉冲响应和波场的多维褶积运算预测多次波.该方法可以从全波场信息中获得一次波脉冲响应,因此即使数据缺失一定数量的地震道EPSI也可以获得完整的、可靠的一次波和多次波结果.
5 实际资料测试大量的模拟测试验证了EPSI方法可以准确有效地预测一次波和多次波,且在数据不完整情况下的也可以很好地重建缺失的一次波.我们进行的实际资料测验采用中国某海域海上地震数据.该数据炮间距50 m、道间距25 m、时间采样间隔8 ms,即每隔一个检波点激发一炮.图 9a展示了该海上资料的某典型单炮记录,图中多次波主要分布在24 s至40 s的时间区间.由于炮间距是道间距的两倍,因此采用常规SRME方法压制多次波时,首先需要通过插值技术将共炮道集加密一倍.
我们在未对炮记录进行任何插值处理的情况下分别用EPSI和SRME对该海上资料直接进行了处理,结果分别展示在图 9b和图 9c中.可以看到无需进行任何插值,EPSI就可以得到好的去除效果.而SRME由于未进行炮插值处理,多次波能量未能全部预测出来,因此一次波结果中仍残留较多的多次波能量.图 10a—10c分别展示了图 9a—9c中多次波分布区域的局部放大,从中可以看到更清晰的对比结果.
通过对比分析SRME和EPSI两种多次波去除方法在标准模型合成数据和海上实际资料的应用效果,我们得到如下的结论:SRME只能单向地从一次波预测出多次波,而EPSI则可以双向地互相预测,即从一次波可以预测多次波,也同时从多次波信息中预测出一次波.SRME处理效果依赖于数据完整性以及一次波多次波的分布情况,而EPSI无论数据完整或缺失、一次波多次波是否重叠,始终可以获得稳定可靠的结果.SRME方法因其算法实现简单,对计算硬件要求不高,在工业界已获得广泛应用,是目前多次波去除的主流技术.但随计算机技术的发展和高性能计算的普及,EPSI计算瓶颈终将被突破,EPSI在多次波去除处理中具有广阔的应用前景.
Bai L S, Liu Y K, Lu H Y, et al.
2014. Curvelet-domain joint iterative seismic data reconstruction based on compressed sensing. Chinese J. Geophys., 57(9): 2937-2945.
DOI:10.6038/cjg20140919 |
|
Berkhout A J. 1982. Seismic Migration:Imaging of Acoustic Energy by Wave Field Extrapolation. Volume A. Theoretical Aspects. 2nd ed. New York:Elsevier Science Publ. Co., Inc., 211-218.
|
|
Berkhout A J, Verschuur D J.
1997. Estimation of multiple scattering by iterative inversion, Part Ⅰ:Theoretical considerations. Geophysics, 62(5): 1586-1595.
DOI:10.1190/1.1444261 |
|
Herrmann F J, Hennenfent G.
2008. Non-parametric seismic data recovery with curvelet frames. Geophysical Journal International, 173(1): 233-248.
DOI:10.1111/gji.2008.173.issue-1 |
|
Jin D G, Chang X, Liu Y K.
2008. Algorithm improvement and strategy of internal multiples prediction based on inverse scattering series method. Chinese J. Geophys., 51(4): 1209-1217.
DOI:10.3321/j.issn:0001-5733.2008.04.032 |
|
Lin T T Y, Herrmann F J.
2013. Robust estimation of primaries by sparse inversion via one-norm minimization. Geophysics, 78(3): R133-R150.
DOI:10.1190/geo2012-0097.1 |
|
Liu Y K, Chang X, Jin D G, et al.
2011. Reverse time migration of multiples for subsalt imaging. Geophysics, 76(5): WB209-WB216.
DOI:10.1190/geo2010-0312.1 |
|
Liu Y K, Hu H, Xie X B, et al.
2015. Reverse time migration of internal multiples for subsalt imaging. Geophysics, 80(5): S175-S185.
DOI:10.1190/geo2014-0429.1 |
|
Liu Y K, Liu X J, Osen A, et al.
2016. Least-squares reverse time migration using controlled-order multiple reflections. Geophysics, 81(5): S347-S357.
DOI:10.1190/geo2015-0479.1 |
|
Robinson E A.
1957. Predictive decomposition of seismic traces. Geophysics, 22(4): 767-778.
DOI:10.1190/1.1438415 |
|
Song J W, Verschuur D J, Chen X H.
2014. Research status and progress in multiple elimination. Progress in Geophysics, 29(1): 240-247.
DOI:10.6038/pg20140134 |
|
Taner M T.
1980. Long period sea-floor multiples and their suppression. Geophysical Prospecting, 28(1): 30-48.
DOI:10.1111/gpr.1980.28.issue-1 |
|
van Groenestijn G J A, Verschuur D J.
2009. Estimating primaries by sparse inversion and application to near-offset data reconstruction. Geophysics, 74(3): A23-A28.
DOI:10.1190/1.3111115 |
|
van Groenestijn G J A. 2010. Estimation of primaries and multiples by sparse inversion[Ph. D. thesis]. Delft:Delft University of Technology.
|
|
Verschuur D J. 1991. Surface-related multiple elimination, an inversion approach[Ph.D. thesis]. Netherlands:Technische Universiteit Delft.
|
|
Weglein A B, Gasparotto F A, Carvalho P M, et al.
1997. An inverse-scattering series method for attenuating multiples in seismic reflection data. Geophysics, 62(6): 1975-1989.
DOI:10.1190/1.1444298 |
|
Ypma F H C, Verschuur D J.
2013. Estimating primaries by sparse inversion, a generalized approach. Geophysical Prospecting, 61(1): 94-108.
|
|
白兰淑, 刘伊克, 卢回忆, 等.
2014. 基于压缩感知的Curvelet域联合迭代地震数据重建. 地球物理学报, 57(9): 2937–2945.
DOI:10.6038/cjg20140919 |
|
金德刚, 常旭, 刘伊克.
2008. 逆散射级数法预测层间多次波的算法改进及其策略. 地球物理学报, 51(4): 1209–1217.
DOI:10.3321/j.issn:0001-5733.2008.04.032 |
|
李鹏. 2007. 复杂介质多次波处理方法研究[博士论文]. 北京: 中国科学院地质与地球物理研究所.
|
|
宋家文, VerschuurD J, 陈小宏.
2014. 多次波压制的研究现状与进展. 地球物理学进展, 29(1): 240–247.
DOI:10.6038/pg20140134 |
|