2018, Vol. 61
2. 中国石油塔里木油田公司, 新疆库尔勒 841000
2. Tarim Oilfield Company, PetroChina, Korla Xinjiang 841000, China
多次波指在地下经过一次以上反射的地震波.多次波一直被人们视为影响地震偏移、层析成像、速度分析的重要因素.多次波会使一次波发生畸变，使地震成像条件更加复杂，从而使地震资料的解释结果出现偏差.多次波分为表面多次波和层间多次波两种.层间多次波是指向下散射只发生在地下界面间的多次波，多见于从表层风化的陆地采集的地震数据.在中东等地区，由于地表附近存在高、低速层垂直交错分布的情况，且相邻层之间的速度差异较大，层间多次波发育良好，甚至存在多次波振幅高于一次反射波的情况(Verschuur and Berkhout, 2015).因此层间多次波压制对于陆地数据处理和解释工作有着不可替代的重要意义.但在通常情况下，层间多次波的动校正量、叠加速度和频率与一次波并无明显差异，这给识别、预测和压制多次波工作带来了极大挑战.
学界对于表面多次波的压制有了很多显著成果(Berkhout and Verschuur, 1997；Weglein et al., 1997；van Groenestijn and Verschuur, 2009；Lopez and Verschuur, 2015)，然而对于层间多次波的讨论深度和广度却明显不如表面多次波.主要原因在于层间多次波的产生机理不清楚，层间多次波的预测难度大、压制难度更大.Weglein等(1997)提出逆散射级数(ISS)的概念和方法，从理论上给出了层间多次波的预测模型.Innanen(2017)针对原始ISS方法静态参数无法处理横向延展较长的散射波产生界面的问题，提出了基于1D和1.5D模型的ISS非静态参数，但并未在完整二维地震道上成功应用.ISS及其后续改进算法普遍受限于庞大的计算量，直到目前仍不能大范围推广到二维和三维叠前数据中.
Jakubowicz(1998)基于Keydar等(1997)的工作提出了利用一次反射波构建一阶层间多次波的设想，并在合成数据中取得了一定效果.然而该方法对于同相轴分离操作的要求太过苛刻，无法在实际数据中得到满足.Hung和Wang(2012)对此做出了有益的改进，按照类似ISS中的“低-高-低”分层方法将单个同相轴的分离操作简化为三个区域的分割，再通过叠加实现全部层间多次波的预测.然而这种方法的问题也是显而易见的，简单叠加多次波模型破坏了它们与真实振幅和相位之间的关系，从而使得多次波模型过分依赖匹配算法，增加了一次反射信号被压制的风险.
Berkhout和Verschuur(2005)将经典的表面多次波去除方法(SRME)推广到了地下界面当中，提出了扩展SRME方法，但该方法高度依赖于共聚焦点道集的构建、宏观速度模型的准确性以及波场延拓.Ypma和Verschuur(2013)基于稀疏反演预测一次波方法(EPSI, van Groenestijn and Verschuur, 2009)在估计子波的同时完成对层间多次波的预测，该方法的优点是基于闭环的反演算法，不需要自适应相减就可以用一次脉冲、多次脉冲和地震子波估计三个分量解释原始数据，但缺点是把EPSI方法推广到层间多次波后，原本就不稳定的子波估计变得更加不稳定，导致方法在实际数据中没有取得显著效果，同时也具有反演结果收敛较慢的问题.近年来广受关注的Marchenko方法和地震干涉法(Wapenaar, 2004；Broggini et al., 2012；Wapenaar and Slob, 2014)在层间多次波方面也有最新的进展，da Costa Filho等(2017)将Marchenko方法在理论上推广到了弹性波场情况，但由于该方法所需的物理量在实际地震数据中并不能测得，因此只在合成数据中应用了简化为声弹性介质的Marchenko方法，且未取得实际数据效果.
此外，近年来国内学者对于层间多次波等各类多次波的压制技术研究也有显著进展(金德刚等，2008；刘伊克等，2008；李翔和胡天跃，2009；吴静等，2013；宋家文等，2014；叶月明等，2015).然而目前学界对于层间多次波去除技术的探索和实际数据应用远未达到成熟水平.虚同相轴概念及其层间多次波压制方法的提出始于Ikelle(2006)，虚同相轴是指在实际地震资料当中并不存在的同相轴.借助虚同相轴可以将散射点从地表移动至产生层间多次波的向下散射界面，从而顺利预测出层间多次波.Ikelle(2006)和Ikelle等(2009)认为通过将虚同相轴与原始数据褶积，就可以得到多次波模型，再通过匹配并逐层重复进行即可压制所有层间多次波.这个思路虽然没有得到很多学者的关注，但其方法具有很高的应用前景(吴静等，2013；Liu and Hu, 2017a, 2017b).然而，该方法并没有建立在严格的数学物理推导基础之上，其理论本身难免存在缺陷.虚同相轴本身包含的多次波信息与原始数据存在大量重复，导致其褶积运算结果中的多次波模型振幅和相位存在扭曲.对于阶次越高的多次波，其在预测模型中的组合种类越多，因此叠加后的振幅理论上会高于那些组合种类更少的低阶多次波.预测模型在振幅和相位上的畸变导致最终的多次波压制匹配效果并不理想；同时，原始虚同相轴方法的同相轴选取要求十分苛刻，这也导致目前尚未有人将其运用到实际数据中.
对于虚同相轴方法在理论和应用中的缺陷，本文通过严格的理论推导，得出基于自适应虚同相轴的层间多次波预测方法.该方法在一维和二维合成数据上的应用效果证明，相比于传统虚同相轴方法，自适应虚同相轴方法可以得到精度更高的多次波预测模型，减少对匹配算法的依赖，多次波压制比例至少提高30%.在PLUTO复杂模型和中国西部实际陆地资料应用的实例证明了自适应虚同相轴方法对去除层间多次波、恢复并突出目标储层同相轴、提高地震成像分辨率的显著作用.
1 自适应虚同相轴理论 1.1 传统虚同相轴理论借助虚同相轴，可以将层间多次波的向下散射点地下移至地表，从而顺利预测出层间多次波.Ikelle(2006)认为只需两步即可构建出层间多次波.假设d0为与某一层地下界面相关的一次反射波，d′0为不含d0的原始数据，那么虚同相轴即可定义为逆时一次反射波与不含该一次波的原始数据的褶积
|
(1) |
其中，上标*代表复共轭.图 1直观地展示了虚同相轴的物理意义，它是一种无法在地震记录中被观察到的波，但可以用来构建与该地下反射界面相关的层间多次波.本文将虚同相轴与切除了该界面一次反射波后的原始数据做褶积，就可以得到层间多次波的预测模型：
|
图 1 传统虚同相轴原理示意图 (a)构造传统虚同相轴；(b)用传统虚同相轴构造层间多次波.叉号表示褶积，全文同. Fig. 1 Sketch of traditional virtual events method (a) Construction of traditional virtual events; (b) Estimation of internal multiples by traditional virtual events. |
|
(2) |
值得注意的是，图 1展现的只是较为简单的一类虚同相轴，实际上传统虚同相轴本身不仅包含一次波，还包含多次波信息，能够预测更高阶的多次波.此外，传统虚同相轴不仅包含一次波的信息，还包含与当前层位有关的多次波信息，本文稍后将讨论这一特性对预测精度的影响并对此进行理论修正.
1.2 自适应虚同相轴理论Berkhout(2012)给出的波场传播矩阵模型(通常简称WRW模型)以传播算子而非传播参数的形式，描述了叠前地震数据离散空间.本文将利用此模型的描述方法给出自适应虚同相轴的理论推导和层间多次波预测模型公式.
在WRW模型中，一组二维地震数据用一个三维矩阵代表，矩阵x分量代表检波点的水平方向坐标xR，y分量代表炮点的水平方向坐标xS，z坐标代表时间t，因此该三维矩阵的每一个t分量下的二维矩阵代表了各个时间点的地震数据切片.沿着t时间轴方向做傅里叶变换，可以将z坐标变换为圆频率ω，从而得到矩阵P.新的三维矩阵P的每一个ω分量代表着一个单频率的地震数据切片.这种表示方式的优势在于其既可以简化描述也可以用于实际编程计算.基于上述表示方法，SRME方法可以表示为(Berkhout and Verschuur, 1997, 2005)
|
(3) |
其中，D表示检波点排列和鬼波效应对上行反射波的改变作用，S表示炮点排列和震源特征，检波点和炮点矩阵括号内的变量z代表检波线和炮线的所在深度.R表示反射算子，^表示反射算子的作用方向朝下，也就是将上行波转换为下行波，括号内的变量z表示反射算子所在深度.X表示脉冲响应，它与地震子波褶积后得到P，P和X矩阵括号内的两个z变量分别代表检波点和炮点所在深度.假设检波点一侧的鬼波已被去除，即D(z0)=I，那么表面多次波生成模型、去除表面多次波后的地震数据体和表面多次波可以分别表示为
|
(4) |
|
(5) |
|
(6) |
值得注意的是，本文在此用了新的表示方法给出SRME迭代模型，以便于与下文表述形式统一.矩阵P和X的上划线和下标0代表在第0层的所有地震波已经被切除(即切除直达波)；矩阵P和X的上标0表示代表与第0层界面相关的多次波——即“表面多次波”——已经去除；对于原始数据体P0来说，上标既可以记为“-1”以与下文保持统一，也可以空缺，本文选择简洁表示方式.P0(z0, z0)虽然不包括表面多次波，但它不仅含有一次波，还包括层间多次波.此外，由于空气和海水的波阻抗差异足够大，(6)式中的向下反射算子一般可以近似地取负单位矩阵；为了不失一般性，在此不作近似.
保持炮点和接收点不变，如果将向下反射算子从z0层下移至zn层，那么(5)和(6)式可改写为针对特定层位的层间多次波扩展SRME模型：
|
(7) |
|
(8) |
对该式的理解可以从两个方面入手.首先暂不考虑括号内标识的激发和接收点，Pnn－1的横线和上标n-1代表去除了表面多次以及与第1~n-1层相关的层间多次波的数据体，这里的“相关”指的是在该层发生过至少一次向下散射，下标n表示第n层及第n层以上的所有反射波被切除；假设Pnn－1已知，此数据体减去与第n层相关的层间多次波Mn即可得到求解目标Pnn，它与子波做反褶积运算后得到不含层间多次波的脉冲响应Xnn.根据向下反射算子的定义，不难得到(8)式中的
因此，本文希望提出一种不依赖速度结构、波场延拓等任何地下信息的完全数据驱动的层间多次波预测方法.该方法只需要地表激发地表接收的地震数据.为了实现这个目标，本文引入Γ(zn, z0)和Γ(z0, zn)分别作为检波点端和炮点端的波场延拓算子，即
|
(9) |
|
(10) |
|
(11) |
为了将预测模型所需数据全部替换为地表激发地表接收的数据，将(9)—(11)式代入(7)和(8)式，从而得到基于自适应虚同相轴的层间多次波模型：
|
(12) |
|
(13) |
其中Vnn为自适应虚同相轴.自适应虚同相轴的定义是逆时一次反射脉冲ΔXn*和不含与该层相关多次波的脉冲响应Xnn的褶积，它表示的是实际地震记录中不存在的波场，它与地震数据做褶积运算即可得到层间多次波.而传统虚同相轴则为逆时一次波ΔPn*(z0, z0)和去掉了该一次波的整体波场Pnn－1(z0, z0)的褶积，后者包含与该层相关的层间多次波.此外，虽然本文利用扩展SRME模型推导自适应虚同相轴方法，但正如之前提到的，两者具有本质的区别.自适应虚同相轴只需要利用地表激发地表接收的数据、不需要宏观速度估计、不需要波场延拓.而扩展SRME则需要波场延拓、宏观地下速度信息，并且需要每次更新CFP算子.
1.3 传统虚同相轴与自适应虚同相轴的本质区别自适应虚同相轴与传统虚同相轴本质主要有以下两个方面不同.首先，自适应虚同相轴只包含一次波信息，而传统虚同相轴包含所有阶次的多次波信息.自适应虚同相轴与地震数据的褶积时，将一次波转换为层间多次波，将n阶层间多次波转换为n+1阶层间多次波.虚同相轴则是由一次反射波和各阶次多次波构建，包含所有阶次波信息，虽然它与地震数据褶积的结果也产生更高阶的多次波，但由于参与褶积的双方都包含所有阶次的多次波，因此并没有上述自适应虚同相轴的严格对应转化关系，这会导致严重的振幅相位畸变(图 2).
|
图 2 自适应虚同相轴模型与传统虚同相轴模型比较 (a)自适应虚同相轴模型；(b)传统虚同相轴模型. Fig. 2 Illustrations of adaptive virtual events model and traditional virtual events model (a) Adaptive virtual events model; (b) Traditional virtual events model. |
其次，自适应虚同相轴的定义基于脉冲信号褶积，而传统虚同相轴则基于常规信号褶积.如果一次反射脉冲响应已知，那么构建的自适应虚同相轴可以直接与原始数据褶积得到实际层间多次波，大大减少对匹配相减算法的依赖.传统虚同相轴则无法做到这一点.关于子波估计，EPSI方法及扩展EPSI方法(van Groenestijn and Verschuur, 2009；Ypma and Verschuur, 2013)提出的反演方法有一定效果，但对于复杂波场数据，其子波估计并不完全可靠.对于实际地震数据，如果无法得到较为可靠的地震子波估计，那么自适应虚同相轴方法仍可以通过含有子波成分的一次反射波褶积构建得到.为了消除子波影响，需要将模型与实际数据做匹配相减(Guitton and Verschuur, 2004；Abma et al., 2005；Fomel, 2008).本文第2节证明，在使用相同匹配相减算法下，新方法的模型匹配结果比传统方法的结果具有更少的多次波残留.
下面直观地解释传统虚同相轴和自适应虚同相轴的本质区别，同时展示新方法如何消除传统方法不可避免的振幅相位畸变.图 2展示了对于任意一种层间多次波，无论其阶次多高，自适应虚同相轴方法都能够给出唯一的褶积构建方法；然而传统虚同相轴则会将所有构建方式简单相加.例如，对于图 2中的三阶多次波，假设第一层分界面为当前向下散射界面.此时，可以利用1个自适应虚同相轴和1个二阶多次波构建这个多次波(图 2a)，也可以利用1个一阶虚同相轴(即含有一阶多次波信息)和1个一阶多次波构建(图 2b)，还可以利用1个二阶多次波虚同相轴和1个一次波构建(图 2b).传统虚同相轴模型的做法是将上述所有组合方式加总(图 2b)，把求和结果作为多次波预测模型，而新的自适应虚同相轴方法则将第一种组合视作多次波预测模型(图 2a)，不需要加总操作.
实际上，对于二阶及二阶以上的层间多次波，阶数越高，构建组合方式就越多.这意味着通过传统虚同相轴方法得到的层间多次波模型是所有组合方式的简单相加，阶次越高的多次波，其相加次数也就越多.在相加的过程中，多次波模型的振幅、相位等信息发生严重畸变，而且它与其他阶次的层间多次波的能量相对关系也偏离了实际数据和真实物理规律.对于振幅来说，由于高阶多次波的组合方式多，参与加总的预测模型多，因此在传统虚同相轴方法中，高阶多次波的振幅会得到加强，而具有较少组合方式的低阶多次波的振幅就会被相对削弱，这与高阶多次波弱于低阶多次波的真实情况不符.不论层间多次波阶次高低，自适应虚同相轴方法却总能给出唯一的一个褶积预测模型(例如图 2)，不存在次数不等的模型组合加总，因而该方法预测的多次波具有更准确的相对振幅大小关系.因此，从物理图像的角度看，对任意阶次多次波都存在一一对应褶积构建关系的自适应虚同相轴模型具有更明确的物理意义和更准确的预测精度，有效避免了传统方法对振幅相位的扭曲作用.本文第2节定量验证了这一重要结论.
1.4 自适应虚同相轴算法构建自适应虚同相轴的困难之处在于它所依赖的一次波信息并不能通过地震数据直接得到，一个可行的办法是将传统多次波压制方法得到的结果作为一次波的近似数据输入自适应虚同相轴模型，另一个方法是以原始数据为起点，逐步迭代逼近求解出一次波，本文选择第二种思路，推导其具体算法.
对于一般实际地震数据，假设子波信息未知，采用匹配算子A(z0)作为反子波算子，则(13)式等价于
|
(14) |
其中暗含了A(z0)=[S(z0)]－2的假设.根据诺依曼级数展开可得
|
(15) |
(15)式与(7)式、(12)式比较可得
|
(16) |
该式表明自适应虚同相轴等价于诺伊曼展开无穷级数之和.这个结论十分重要，因为它表明了只要从初始值开始，通过自适应迭代的方式就可以逐步逼近自适应虚同相轴的真实值.综合考虑(14)—(16)式，得到了基于迭代构建虚同相轴的层间多次波预测算法：
|
(17) |
其中，m表示迭代次数，m=0表示初始值.(17)式的证明并不困难，按照迭代过程运算求和即可得到(15)式所示的无穷级数.如果令n=0、自适应虚同相轴取负单位矩阵，公式退化为表面多次波模型.值得注意的是，经过(17)式迭代过程，构建的自适应虚同相轴和层间多次波均与第n层界面(即当前向下散射界面)有关，因此并不包含其他层位的多次波.通过下移层位重新定义向下散射界面并重复上述过程即可预测并压制所有层间多次波.
2 数值验证 2.1 一维数值验证为检验自适应虚同相轴模型的准确性，本节通过一个已知一维模型产生的合成数据展示新方法与传统虚同相轴方法结果的优劣，同时定量给出效果提升幅度.图 3a给出了速度模型，包含4个强波阻抗差异界面，2个速度反转区域，有利于产生较强层间多次波.模型密度为常数2000 kg·m-3.为了排除表面多次波的影响，正演计算采用吸收边界.图 3(b—c)给出了原始道集，并用箭头标出了一次反射波，总道长3500 ms，采样率1 ms.从中可以看出，层间多次波的振幅较强，个别走时与一次波接近，传统方法去除难度大.图 3d给出了传统虚同相轴方法预测的一次波，可以看到大部分多次波能量已被压制，但残余能量有可能干扰后续对一次波的判断.图 3e给出了在相同匹配算法条件下的自适应虚同相轴方法的压制结果，可以看到绝大部分多次波能量均已被压制，相比图 3f的真实一次波，新方法的结果中只有极个别的多次波残余.两组结果使用的全局匹配算子长度均为21 ms，且均对200 ms局部窗口采用3 ms局部算子作局部振幅微调.较短的匹配算子从一个侧面说明了新方法相较于传统方法，减少了对匹配相减的依赖.
|
图 3 自适应虚同相轴应用于一维层间多次波压制结果 (a)速度模型；(b)原始记录(箭头表示一次波)；(c)放大后的原始记录；(d)传统虚同相轴处理结果；(e)自适应虚同相轴处理结果；(f)真实一次波. Fig. 3 Results of adaptive virtual events on a 1D synthetic data set (a) Velocity model; (b) Raw data(arrows indicate primaries); (c) Enlargements of raw data; (d) Result by traditional virtual events method; (e) Results by adaptive virtual events method; (f) True primaries. |
定量统计原始道集和一次波预测结果中的多次波能量，能量以地震道绝对值之和计.将真实层间多次波能量作为分母，将经过层间多次波压制、且减去了真实一次波后的地震记录的能量作为分子，通过除法得到残余能量相对于真实多次波能量的比例.用100%减去这个数值即可得出多次波去除比例.比例越高，意味着残余能量越少，也就意味着方法对一次波的损伤小、引入的额外噪声少、预测多次波的精度高.相较于把预测多次波能量与真实多次波能量简单相除的做法，这个指标更能客观全面地反映真实压制情况.经过上述计算，传统虚同相轴方法去除61.09%的多次波，而新的自适应虚同相轴算法可以去除高达91.06%的多次波能量.相比传统方法，新方法的多次波压制比例提高了30%，效果显著.
2.2 二维数值验证本节将一维模型扩展到二维情况，检验新的自适应虚同相轴方法的多次波压制效果，并与传统方法对比分析.本节继续沿用上一节的模型，并将其扩展为横向长4.8 km的二维平层模型.炮点间隔和检波点间隔均为20 m，采用240×240个台站的激发和观测系统，四周采用吸收边界消除表面多次波影响.图 4a给出了其中一炮的原始记录，可以看到层间多次波振幅依然很强，由于加入了偏移距的因素，多次波与一次波的走时交叉更多，去除难度更大.两种方法均使用了相同的全局匹配算子，长度为31 ms，未使用局部算子.图 4b给出了用传统方法去除第1层相关的多次波的结果，图 4c是用新方法去除第1层相关多次波的结果.比较图 4b和4c可以看出，新方法对于多次波的压制更加彻底，效果更加显著.图 4d是用新方法去除所有多次波的最终结果，不难看出除了少许残留，大部分多次波已被压制.自适应虚同相轴方法在二维情况的效果得到成功验证.
|
图 4 自适应虚同相轴应用于二维层间多次波压制结果 (a)原始炮集；(b)用传统虚同相轴方法去除第1层相关多次波结果；(c)用自适应虚同相轴方法去除第1层相关多次波结果；(d)用自适应虚同相轴方法去除所有多次波后的最终结果. Fig. 4 Results of adaptive virtual events on a 2D synthetic data set (a) Raw data; (b) Results of 1st-layer-related internal multiples attenuation by traditional virtual events method; (c) Results of 1st-layer-related internal multiples attenuation by adaptive virtual events method; (d) Final result by adaptive virtual events method. |
图 5a和5b展示的是与第1层相关的自适应虚同相轴构建过程.图 5c和5d展示的是与第2层相关的自适应虚同相轴构建过程自适应虚同相轴构建过程，可以看出迭代计算的过程就是将与该层相关的多次波信息不断消除的过程.自适应虚同相轴顶部的叉型噪声来自有限的积分孔径，该噪声对消除层间多次波的处理过程没有影响，且可以用常规滤波方法消除.
|
图 5 自适应虚同相轴构建过程图 (a)与第1层相关的虚同相轴；(b)与第1层相关的自适应虚同相轴；(c)与第2层相关的虚同相轴；(d)与第2层相关的自适应虚同相轴. Fig. 5 Illustration of construction of adaptive virtual events (a) Virtual events related to the first layer; (b) Adaptive virtual events related to the first layer; (c) Virtual events related to the second layer; (d) Adaptive virtual events related to the second layer. |
本节将新的自适应虚同相轴方法运用于复杂合成数据并进行验证.本节将该方法运用于SMAART团队公开发布的多次波测试模型PLUTO.在运用SRME方法去除了表面多次波之后，采用自适应虚同相轴方法进行层间多次波压制.图 6以第637炮数据为例，给出了用自适应虚同相轴方法去除层间多次波前后的效果对比，箭头指出了层间多次波的主要分布位置.不同于表面多次波，层间多次波往往难以在炮域上直接通过周期性识别出来，因此传统依靠走时或变换切除的方法很难奏效.运用自适应虚同相轴方法，可以在不识别多次波走时的情况下有效压制多次波；同时，对比放大图 6d—e可以看出在压制多次波的同时，一次波得到了完整保留，一次波没有出现任何损失，没有引入新的噪声.在这里只选择部分炮集而非全部炮集参与运算，这一点也有别于表面多次波的预测.
|
图 6 自适应虚同相轴应用于PLUTO模型压制层间多次波结果 (a)经过表面多次波压制的第637炮数据；(b)经过表面多次波和层间多次波压制的第637炮数据; (c)自适应虚同相轴方法预测的第637炮多次波；(d)(e)(f)分别为(a)(b)(c)的局部放大图.箭头表示层间多次波. Fig. 6 Results of adaptive virtual events on PLUTO data set (a) Shot record 637 after surface multiples attenuation; (b) Shot record 637 after surface and internal multiples attenuation; (c) Estimated internal multiples in shot 637 by adaptive virtual events method; (d) (e) and (f) are enlargements of (a) (b) and (c) respectively.Arrows indicate internal multiples. |
由于不知道一次波和多次波的理论走时，且层间多次波与一次波的部分特征十分相近，因此从炮集上无法直接判断去除的成分是否是多次波，需要借助图 7给出的叠加剖面和多次波叠加剖面结果，结合从原始剖面读取到的地层基本特征，综合判定多次波压制效果.图 7用箭头标出了被压制的层间多次波.从剖面结果看，尤其是结合图 7c给出的噪声剖面看，被去除的噪声具有明显的多次波特征，它们具有与“眼眉状”盐丘上方边界和下方边界相似的空间延展特征，横向分布涉及的CMP范围也与盐丘类似.同时，自适应虚同相轴方法对于层内短路径多次波也具有良好的预测压制效果.例如在CMP 1627、2.8 s附近的与真实地层同相轴具有相近的走时的层间多次波同样能够被成功地预测和压制.此外从结果中还可以看出该方法对于高倾角多次波的预测压制能力，例如在CMP 2107、5.1s附近具有高倾角的深部盐丘层间多次波被明显压制.
|
图 7 自适应虚同相轴应用于PLUTO模型压制层间多次波结果 (a)去除表面多次波后的叠加剖面；(b)去除表面和层间多次波后的叠加剖面；(c)自适应虚同相轴方法预测的多次波剖面.箭头表示层间多次波. Fig. 7 Results of adaptive virtual events on PLUTO data set (a) Stacked section after surface multiples attenuation; (b) Stacked section after surface and internal multiples attenuation; (c) Stacked internal multiples estimated by adaptive virtual events method. Arrows indicate internal multiples. |
最后，将自适应虚同相轴方法运用于陆地实际地震资料的多次波压制，展该方法在去除层间多次波、恢复并突出目标储层同相轴、提高地震成像分辨率方面的显著效果.数据来自中国西北部地区，最大偏移距近6 km，目的储层深度最深处可达8 km，采集区域总长度达74 km，层间多次波发育、振幅强，且层间多次波的叠加速度和走时与一次波十分相近，这类多次波往往在剖面上“附着”于一次波附近，传统方法对这类多次波的去除效果并不理想，因此尝试用自适应虚同相轴方法进行多次波压制.
首先选取CMP范围，定义参与多次波预测的炮集数据范围.CMP范围的选定依据是，被选取部分包含较为连续的、波阻抗差异较大的能够产生向下散射的关键界面，这样给出的多次波预测模型才能有坚实的物理产生机制.选取的应用范围是垂直长度8 km、横向延展20 km的广大地下空间，涉及676炮数据.其次根据炮集和原始叠加剖面观察的结果，选择1.10、1.41、2.41 s和3.45 s附近的同相轴作为四个关键散射体，依次迭代计算虚同相轴和与之相关的层间多次波.这一点与之前的合成数据差异很大.一般来说，实际地层分层十分精细，同相轴数目很多，尤其是浅层数据尤为如此.浅层同相轴虽然对于该算法极为重要，但浅层数据质量普遍不高，同相轴连续性差，因此需要选取有代表性的连续同相轴作为其他层位的近似.如果考虑时间成本，也必须选择那些近似可替代其他散射面的强分界面作为散射算子.此外，在应用自适应虚同相轴方法之前，静校正、振幅恢复、去噪等一些预处理技术将有助于提升预测模型的准确度和多次波压制的最终效果.
图 8a给出了去噪后的原始炮集记录，从记录中可以看出，虽然该地区数据的整体信噪比较高，但具有多次波特征的同相轴发育明显，且出现了精细分层现象.相邻同相轴的视速度和走时非常接近，因此难以通过肉眼区分一次波和多次波.同时，同相轴数量众多，走时差异较小，同相轴之间的边缘分界不明显，很难用时窗将单条同相轴与原始数据完美分离，这进一步增加了去除多次波的难度，然而自适应虚同相轴方法却在不伤害有效波的情况下给出了很好的压制效果(图 8b).从炮集上能够观察到许多令人感兴趣的特征，例如图 8c所示的被压制的层间多次波具有明显的周期性特征和精细分层特征，多次波分布区域遍布整个时间域的地震道，无论周围是否有一次波，都有振幅较强的多次波发育.从图 8b压制结果中可以看出，多次波泄露的问题依然存在.出现这个现象的首要原因是经过新方法改进后的多次波模型的振幅和相位与真实值仍有一定差距；其次，陆地资料噪声和深部数据较低的信噪比也是影响多次波压制的重要因素；第三，部分多次波的走时与目的层位一次波较为相近，去除这部分多次波且不伤害有效波的难度较大.从图 8c的多次波炮集上看，一次波附近区域并无明显强于邻近区域的能量集中区，说明一次波信号得到了很好地保留，并未被算法误判为多次波.
|
图 8 自适应虚同相轴压制实际陆地层间多次波结果 (a)预处理后的第338炮数据；(b)经过层间多次波压制的第338炮数据; (c)自适应虚同相轴方法预测的第338炮多次波.黑色箭头表示目的层位一次波，浅色箭头表示主要层间多次波. Fig. 8 Results of adaptive virtual events on land field data (a) Shot record 338; (b) Shot record 338 after internal multiples attenuation; (c) Estimated internal multiples in shot 338 by adaptive virtual events method.Dark arrows indicate the primaries of objective. Light arrows indicate the internal multiples. |
实际数据的处理难度虽然比合成数据高出很多，但它带来的结果和启示却远比合成数据丰富.炮集上的结果分析并不足以证明方法的有效性，因此在图 9中分别给出了原始叠加剖面、自适应虚同相轴方法压制后的剖面、多次波剖面以及传统方法的压制结果，在图 10中给出了图 9a和9b的深部放大图.(1)虽然参与计算的数据在时空域分布很广，但自适应虚同相轴算法在整个叠加剖面中都取得了十分显著的多次波压制效果.(2)大量与一次波振幅相当、走时相近的多次波同相轴被去除，例如在1.5~3.5 s之间很多层间多次波振幅与精细分层一次反射波相当，且走时也重合(图 9a—b).(3)一次波同相轴的连续性得到大幅改善(见图 10圆圈标记处).(4)与层间多次波产生机理类似的散射波和绕射波也被充分压制，这些散射波主要分布在3.5 s和5.0 s等主要一次波反射轴附近，这些散射的形态特点不同于平层同相轴，具有十分鲜明的弧形形态和绕射波的特点(图 9c).(5)剖面分辨率和质量也随着多次波的去除而明显提高了，例如在2.5 s、3.5 s和5.0 s附近的一次波信号在去除多次波后得到了加强，分辨率得到了提高.这是由于自适应虚同相轴方法能够在有效压制多次波的同时，尽可能保留和恢复一次波的振幅相位信息，最终能够使去除多次波畸变影响后的一次波实现同相叠加.(6)在多次波被压制的一些区域，例如1.5~2.5 s和5~6 s的全CMP范围内，原始数据中被多次波掩盖的真实层位同相轴得到显著恢复，这些同相轴的视倾角小于层间多次波.
|
图 9 自适应虚同相轴压制实际陆地层间多次波结果 (a)预处理后的叠加剖面；(b)经过自适应虚同相轴方法层间多次波压制后的叠加剖面；(c)自适应虚同相轴方法预测的层间多次波叠加剖面；(d)传统虚同相轴方法压制多次波后的剖面. Fig. 9 Results of adaptive virtual events on land field data (a) Stacked section before multiples attenuation; (b) Stacked section after internal multiples attenuation by adaptive virtual events method; (c) Stacked internal multiples estimated by adaptive virtual events method; (d) Stacked section after internal multiples attenuation by traditional virtual events method. |
|
图 10 自适应虚同相轴对实际陆地资料一次反射同相轴的改善效果 (a)预处理后的叠加剖面；(b)经过层间多次波压制的叠加剖面.圆圈表示一次波同相轴连续性改善的区域，浅蓝色箭头表示被压制的层间多次波，深蓝色斜箭头表示被恢复的真实层位，恢复层位的视倾角小于层间多次波. Fig. 10 Improvement results of primaries in land field data by adaptive virtual events method (a) Stacked section before multiples attenuation; (b) Stacked section after internal multiples attenuation.Circles indicate the improved continuity of primary events. Light arrows indicate the suppressed internal multiples. Dark arrows indicate the recovered subsurface layers which have smaller apparent dip than internal multiples. |
图 9d给出了传统虚同相轴方法在同一组陆地实际资料上的多次波压制剖面结果，可以看出虽然5 s附近的部分层间多次波得到压制，但剖面整体并未得到明显改观，被多次波覆盖的层位并未恢复，5 s附近存在横向间断的一次波同相轴连续性也未得到提升.
实际陆地地震资料实例证明了自适应虚同相轴方法在去除层间多次波、恢复并突出目标储层同相轴、提高地震成像分辨率方面的显著作用.
4 讨论自适应虚同相轴方法利用简单明了的物理图像修正了传统虚同相轴模型，依靠坚实的理论推导建立起了新的模型.从理论上看，由于新模型公式中的一次波信息未知，为了避免非线性反演诺依曼展开系数，提出了基于迭代计算的求解思路.虽然迭代公式是严格的，但由于实际总会面临迭代次数有限和迭代收敛速度不够的问题，总会存在一些近似.上述数据测试结果表明，自适应虚同相轴方法只需3至4次迭代即可收敛.
从实现方式上来看，虽然已经通过上述工作将人工介入减少到最小，该算法实质上未摆脱人工操作的依赖.目前算法只需通过观察炮集和原始叠加剖面，给出向下散射界面的大致时窗即可使算法自行计算得到本文的所有结果.界面选取影响结果.
最后讨论数据规则化的问题.实际上，对于本文采用的PLUTO数据和陆地实际地震资料，其原始状态都不满足数据规则化要求，即要求每一个接收点都必须有一组以它为炮点的共炮道集.实际二维陆地地震资料的检波点密度一般是炮点密度的倍数，而三维实际数据的炮点更加稀疏，这就要求在共检波点域进行插值，使得炮检点密度相同.另一种普遍存在的情况是检波点密度也不满足基尔霍夫积分精度要求，此时需要在共炮点域进行插值.此外，海上数据还可能面临浅层近偏移距数据缺失的问题.本文的后两组结果都是在简单插值的基础上得到的，说明自适应虚同相轴方法对于数据规则化的精度要求并不高.
5 结论本文通过严格的理论推导，对传统虚同相轴方法进行了理论修正，提出了新的自适应虚同相轴理论，提高了层间多次波预测模型的精度，减少了对匹配相减的依赖.将该方法用于一维和二维合成数据的层间多次波压制，定量统计得出传统虚同相轴方法去除61.09%的多次波，而新的自适应虚同相轴算法可以去除高达91.06%的多次波能量.相比传统虚同相轴方法，本文研究的自适应虚同相轴方法可以显著提升多次波压制效果.在PLUTO复杂模型和中国西部实际陆地资料应用的实例证明了自适应虚同相轴方法对去除层间多次波、恢复并突出目标储层同相轴、提高地震成像分辨率的显著作用.
Abma R, Kabir N, Matson K H, et al.
2005. Comparisons of adaptive subtraction methods for multiple attenuation. The Leading Edge, 24(3): 277-280.
DOI:10.1190/1.1895312 |
|
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 |
|
Berkhout A J, Verschuur D J.
2005. Removal of internal multiples with the common-focus-point (CFP) approach:Part 1-Explanation of the theory. Geophysics, 70(3): V45-V60.
DOI:10.1190/1.1925753 |
|
Berkhout A J. 2012.
Seismic Migration:Imaging of Acoustic Energy by Wave Field Extrapolation. New York: Elsevier.
|
|
Broggini F, Snieder R, Wapenaar K.
2012. Focusing the wavefield inside an unknown 1D medium:Beyond seismic interferometry. Geophysics, 77(5): A25-A28.
DOI:10.1190/geo2012-0060.1 |
|
da Costa Filho C A, Meles G A, Curtis A.
2017. Elastic internal multiple analysis and attenuation using Marchenko and interferometric methods. Geophysics, 82(2): Q1-Q12.
DOI:10.1190/geo2016-0162.1 |
|
Fomel S.
2008. Adaptive multiple subtraction using regularized nonstationary regression. Geophysics, 74(1): V25-V33.
|
|
Guitton A, Verschuur D J.
2004. Adaptive subtraction of multiples using the L1-norm. Geophysical Prospecting, 52(1): 27-38.
DOI:10.1046/j.1365-2478.2004.00401.x |
|
Hung B, Wang M. 2012. Internal demultiple methodology without identifying the multiple generators. //82nd Ann. Internat Mtg., Soc. Expi. Geophys. . Expanded Abstracts, 1-5.
|
|
Ikelle L T.
2006. A construct of internal multiples from surface data only:the concept of virtual seismic events. Geophysical Journal International, 164(2): 383-393.
DOI:10.1111/gji.2006.164.issue-2 |
|
Ikelle L T, Erez I, Yang X J.
2009. Scattering diagrams in seismic imaging:More insights into the construction of virtual events and internal multiples. Journal of Applied Geophysics, 67(2): 150-170.
DOI:10.1016/j.jappgeo.2008.10.009 |
|
Innanen K A.
2017. Time-and offset-domain internal multiple prediction with nonstationary parameters. Geophysics, 82(2): V105-V116.
DOI:10.1190/geo2016-0220.1 |
|
Jakubowicz H. 1998. Wave equation prediction and removal of interbed multiples. //60th EAGE Conference and Exhibition. EAGE, 1527-1530.
|
|
Jin D G, Chang X, Liu Y K.
2008. Algorithm improvement and strategy of internal multiples prediction based on inverse scattering series method. Chinese Journal of Geophysics, 51(4): 1209-1217.
|
|
Keydar S, Lande E, Gurevich B, et al. 1997. Multiple prediction using wavefront characteristics of primary reflections. //59th EAGE Conference & Exhibition, Extended Abstracts. EAGE.
|
|
Li X, Hu T Y.
2009. Surface-related multiple removal with inverse scattering series method. Chinese Journal of Geophysics, 52(6): 1633-1640.
|
|
Liu J, Hu T. 2017a. Inversion based internal multiple attenuation study. //International Geophysical Conference, Soc. Expi. Geophys. and Chinese Petroleum Society, Expanded Abstracts, 561-564.
|
|
Liu J, Hu T. 2017b. Internal multiple attenuation by iterative construction of virtual events. //79th EAGE Conference & Exhibition, Extended Abstracts. EAGE.
|
|
Liu Y K, Chang X, Wang H, et al.
2008. Internal multiple removal and its application by wavepath migration. Chinese Journal of Geophysics, 51(2): 589-595.
|
|
Lopez G A, Verschuur D J.
2015. Closed-loop surface-related multiple elimination and its application to simultaneous data reconstruction. Geophysics, 80(6): V189-V199.
DOI:10.1190/geo2015-0287.1 |
|
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 |
|
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 |
|
VerschuurD J, Berkhout A J.
2015. From removing to using multiples in closed-loop imaging. The Leading Edge, 34(7): 744-759.
DOI:10.1190/tle34070744.1 |
|
Wapenaar K.
2004. Retrieving the elastodynamic Green's function of an arbitrary inhomogeneous medium by cross correlation. Physical Review Letters, 93(25): 254301.
DOI:10.1103/PhysRevLett.93.254301 |
|
Wapenaar K, Slob E.
2014. On the Marchenko equation for multicomponent single-sided reflection data. Geophysical Journal International, 199(3): 1367-1371.
DOI:10.1093/gji/ggu313 |
|
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 |
|
Wu J, Wu Z Q, Hu T Y, et al.
2013. Seismic internal multiple attenuation based on constructing virtual events. Chinese Journal of Geophysics, 56(3): 985-994.
DOI:10.6038/cjg20130326 |
|
Ye Y M, Yao G S, Zhao C L, et al.
2015. Sea bottom peg-leg multiple suppression based on seismic interferometry. Oil Geophysical Prospecting, 50(2): 225-231.
|
|
Ypma F H C, Verschuur D J.
2013. Estimating primaries by sparse inversion, a generalized approach. Geophysical Prospecting, 61(s1): 94-108.
|
|
金德刚, 常旭, 刘伊克.
2008. 逆散射级数法预测层间多次波的算法改进及其策略. 地球物理学报, 51(4): 1209–1217.
|
|
李翔, 胡天跃.
2009. 逆散射级数法去除自由表面多次波. 地球物理学报, 52(6): 1633–1640.
|
|
刘伊克, 常旭, 王辉, 等.
2008. 波路径偏移压制层间多次波的理论与应用. 地球物理学报, 51(2): 589–595.
|
|
宋家文, VerschuurD J, 陈小宏.
2014. 多次波压制的研究现状与进展. 地球物理学进展, 29(1): 240–247.
DOI:10.6038/pg20140134 |
|
吴静, 吴志强, 胡天跃, 等.
2013. 基于构建虚同相轴压制地震层间多次波. 地球物理学报, 56(3): 985–994.
DOI:10.6038/cjg20130326 |
|
叶月明, 姚根顺, 赵昌垒, 等.
2015. 利用地震干涉法衰减海底相关层间多次波. 石油地球物理勘探, 50(2): 225–231.
|
|