2017, Vol. 60
在地震成像(Claerbout，1985)与反演过程中，传统的地震成像方法对含有高阻抗上覆地层的区域往往由于强散射、高速屏蔽作用导致照明不足、成像效果较差.低频信息对于地下介质的构造成像尤为重要，低频信号具有较强的穿透高速屏蔽层的能力，可以利用低频信号来降低高速屏蔽层对地震信号的散射和屏蔽作用，提高深层信号的信噪比.实际应用中，由于被动源数据含有丰富的低频信息，常常通过地震干涉法对被动源数据进行处理，得到有用的低频信息.Claerbout提出地震干涉理论，证明了水平层状介质条件下，对在地表接收到的从底部传来的透射波记录进行自相关运算与其自激自收的模拟记录等价.Schuster (2001)正式将该方法命名为“地震干涉技术”，在其发表的一系列文章中对该方法的原理和应用价值进行了系统阐述，并使用该方法对多次波进行偏移成像，得到了反射界面的构造特征.Calvert (2004)提出“虚震源法”，可以将震源重置到检波点位置，该方法在勘探地震数据上得到了广泛应用.Nakata等(2011)用地震干涉技术对实际数据处理和分析，对地震干涉法的实用性进行了大量研究.Wapenaar等(2012)用积分理论和互易定理证明了地震干涉技术对不同的非衰减介质和不同的震源条件下均成立，并用理论和实际数据进行了验证.国内的研究人员也对地震干涉法开展了研究工作.朱恒等(2012)对地下介质随机分布的被动源发出的透射波记录合成虚拟炮集，探索了被动源成像的可行性.许卓(2012)认为被动源和主动源地震干涉法能够对重构前的数据予以补充和改善.张盼等(2015)对被动源地震数据进行了地震干涉处理，并研究了被动源重构记录与主动源地震记录的联合插值方法.叶月明等(2015)研究了利用地震干涉法进行了衰减海底层间多次波的研究.
常规的地震干涉法(Seismic Interferometry，SI)构建出虚拟震源必须置于检波器的位置上，这极大地限制了地震干涉法的应用范围，而地震超越干涉法(Beyond Seismic Interferometry，BSI)不同于常规的地震干涉法，其实现了在无需知道地下介质信息的情况下在非检波器位置构建虚拟震源的新方法，即虚拟震源的位置可以分布在地下介质的任意位置.BSI最早由Filippo Broggini和Kees Wapenaar于2012年提出，它利用地表接收到的主动源反射响应和在不精确的平滑模型上估计的直达波来重构包含一次波和多次波的格林函数场(Broggini et al., 2012; Wapenaar et al., 2013).利用重构得到的上传和下传格林函数进行Marchenko成像，能够在无需了解上覆介质构造的情况准确对目标区域的介质进行精确成像.基于单次散射理论的传统偏移成像方法的构造假象来自于其不能很好地解释层间多次波反射，而Marchenko成像很好地解释了多次波，Broggini等(2014)和Slob等(2014)通过计算重构出的上传和下传格林函数场的互相关函数或者多维反褶积(Wapenaar et al., 2008)构造出没有层间多次波假象的成像结果.Marchenko成像的成像结果比传统的互相关逆时偏移成像(Baysal et al., 1983)结果更加准确.Meles (2015)结合了SI和Marchenko成像进行层间多次波预测和消去.本文的方法也是二者的结合，不过重点放在SI低频信息对主动源BSI低频信息的重构上.
基于低频缺失的地震数据重构地下介质处的格林函数会导致在最终的Marchenko成像结果中出现类似于虚假反射界面的地下界面不够聚焦的问题，这会对后续地震勘探的解释造成影响.现有的地震数据低频重构研究中，一类是从数学计算的方式得到低频信息：韩立国等(2012)基于压缩感知和稀疏反演方法重构数据低频信息，并在FWI中进行了应用，取得了一定的效果；Shen (2013)等研究了基于高灵敏度检波器数据的对高频端进行分频匹配滤波的拓频方法；黄超等(2015)用弹性波包络法恢复了地下介质纵横波速度的长波长成分；胡勇等(2017)提出解调包络方法重构地震记录中的低频信息并应用于全波形反演中.此类低频重构方法在实际应用中可能在地震记录上产生虚假信息，最终影响成像结果.另一类是从被动源数据中提取低频信息：Zhang等(2015)、Alali等(2016)从能量角度(振幅)匹配主动源和被动源数据，但即使经过去噪之后的高信噪比低频数据中依然存在噪声，这种局限性也会对最终的成像结果产生影响.另外，此类低频补偿方式是基于权重函数的补偿，通过人为调整能量的增益系数来进行低频补偿，主观性较强而且不能基于数据驱动.由于Marchenko成像中同时利用了一次波和多次波的信息，同时迭代求解Marchenko方程重构格林函数需要宽频带的地震数据作为输入，所以低频信息对Marchenko成像非常重要.这与传统的成像处理中压制低频、增强高频能量进而提高地震资料分辨率有本质的区别.
基于对SI与BSI假设和理论近似的分析，本文提出地震超越干涉法与地震干涉法的结合方法SIBSI，首次实现被动源SI低频对主动源BSI的低频重构，进而基于重构后的宽频带地震数据在平滑的模型上进行格林函数重构，最终对比重构前后的Marchenko成像效果.本文提出基于频率优势的主动源低频重构方法，并给出了SI低频信息重构主动源BSI低频信息的完整流程，该方法在较完整地保留高频信息的基础上，很好地重构了低频信息，同时由于匹配滤波的选择性通过作用，最终的Marchenko成像结果中没有出现由被动源低频噪声导致的虚假反射界面.文章通过数值算例给出对重构后宽频带数据的预处理过程，消去了自由表面多次波和直达波，并将预处理得到的符合Marchenko迭代算法的低频重构后的数据作为输入重构出格林函数，最终的Marchenko成像结果中，目标成像区域的界面更加聚焦，由低频缺失导致的旁瓣效应消失，验证了方法的有效性.由于本文提出的低频重构方法与基于Marchenko方程的格林函数重构都是数据驱动的，所以本文的方法是完全数据驱动的.Marchenko成像涉及巨大的计算量，文章在结论部分对方法的效率进行讨论并给出优化方法.
2 方法原理 2.1 常规地震干涉法理论假设及误差来源分析地震干涉法理论是基于互易定理推导得来的，假设在非均匀介质中的震源x和检波点分别为x0、x″0的格林函数为G (x0, x, t)、G (x″0, x, t)(如图 1)，其频率域形式是G (x0, x, ω)、G (x″0, x, ω)，那么此时互相关法重构在频率域(Wapenaar et al., 2004)可以表示为
|
图 1 互相关法格林函数重构 Fig. 1 Green′s function retrieval by cross-correlation |
|
(1) |
其中，
|
(2) |
式中x0和x″0在∂D内.(2) 式是在任意无损非均匀声波介质的精确表达式，通过上式重构出的格林函数G (x0, x″0)包含了x0和x″0之间的直达波、一次波、非均匀性导致的层间多次波以及自由表面多次波.
在实际的应用当中，为了重构两个检波点之间的格林函数，需要对式(2) 进行近似：1、∂D之外是均匀介质；2、远场近似；3、围绕的介质参数均匀变化(Wapenaar and Follema, 2006).经过近似后所得的公式为
|
(3) |
理论上这些近似会造成振幅错误(Ramírez and Weglein, 2009)，但是近似并不影响式(3) 的相位，所以误差被认为是可接受的(Thorbecke and Draganov, 2011).
地震干涉法基于双边聚焦的假设来重构格林函数.如果中包含自由表面，那么在式(2)、式(3) 中，只需要在剩余不包含自由表面的部分有震源存在就能满足SI的震源包围检波器假设.
2.2 基于频率优势的低频重构方法本文提出了基于频率优势的低频重构方法，即基于被动源的低频优势对主动源低频进行重构.基于被动源在低频端相对于主动源的优势，基本的思路为：匹配被动源和主动源的振幅(能量)，寻找被动源和主动源低频位置的频谱交叉点P，使用低通滤波器，滤出被动源、主动源P点之前的低频部分，用此时的主动源低频部分去对被动源低频部分进行匹配滤波进行主动源低频重建.然后归一化后求和，最后进行振幅能量恢复，得到宽频带地震数据.同时，主动源低频去匹配被动源时，由于匹配滤波的信号选择性通过作用，在匹配滤波的期望输出里不会引入太多的低频噪声.
假设被动源数据经过多域迭代去噪处理(张盼等, 2015)，得到了信噪比较高的被动源重构记录.此时被动源记录记为P (t)，主动源记录记为A (t).在频率域，对P (ω)，A (ω)进行能量匹配，使能量增强后的被动源记录P (ω)与主动源A (ω)记录的频谱曲线交叉于低频端的fcut处.结合低通滤波器对主动源记录A (ω)进行低通滤波得到滤波后的主动源数据Afcut(ω).结合低通滤波器对被动源记录P (ω)进行低通滤波得到滤波后的被动源数据Pfcut(ω).
首先用Afcut(t)对Pfcut(t)进行匹配滤波，设计匹配滤波器m (t)作用于Afcut(t)，使其满足：Lfcut(t)=Afcut(t)*m (t)，*表示褶积.
假设mi(t)是第i个滤波算子，作用于低通滤波后主动源记录的第i道Afcut, i(t)，此时的匹配滤波期望输出为
|
(4) |
定义误差函数
|
(5) |
利用最小二乘原理，令总误差能量对mi(t)的倒数等于0可以得到求解匹配滤波算子的Toeplite矩阵方程(段云卿, 2006):
|
(6) |
其中Rxr为输入道Afcut, i(t)的自相关矩阵，M为匹配滤波算子向量，Rzr为期望输出道Lfcut(t)与输入道Afcut, i(t)的互相关函数向量.求解式(6) 可以得到匹配滤波算子mi(t).定义相关性好、信噪比高的地震道计算出的N个算子进行平均，可以得到所需要的匹配滤波算子
|
(7) |
对Lfcut(t)与A (t)归一化后，并基于主动源在高频段的能量重构宽频带地震数据：
|
(8) |
简化上式可得，
|
(9) |
由此可知，基于频率优势的低频重构完整地保留了高频部分，同时基于归一化，低频能量得到了重构.
2.3 超越地震干涉理论Marchenko方程是一维逆散射问题的基本方程，Wapenaar (2012, 2014)将其扩展到三维形式，三维Marchenko方程(Wapenaar et al., 2014)使得在已知地表接收到的反射响应时可以重构出地下某点虚拟震源的格林函数场.由于传统的地震干涉法需要在物理检波器的位置重构虚拟震源，所以此方程极大地扩展了地震干涉理论，被称为地震超越干涉法(Wapenaar et al., 2014).地震超越干涉法基于无损介质和单边聚焦的假设，仅需要已知地表一侧接收到的反射地震响应以及估计的地下虚拟震源到地表检波器之间格林函数的直达波.SI与BSI的理论对比如表 1.借助三维Marchenko方程重构的格林函数包含了准确的一次波和由于地下介质不均匀性导致的层间多次波.根据互易定理可知，此时的虚拟震源可以看作震源位于地表的虚拟检波器.通过对三维Marchenko方程进行分解，可以得到该虚拟检波器位置的上传和下传格林函数.本文基于互相关函数，求取虚拟检波器位置的反射系数，并对目标地区求取一定间隔的离散虚拟震源位置的上传、下传格林函数进行面向目标的成像，称之为Marchenko成像.由于层间多次波造成的假象在Marchenko成像中得到了压制，所以最终得到不含虚假反射界面的成像结果.
|
表 1 SI与BSI的理论对比 Table 1 Theoretical comparison between SI and BSI |
基于本文的意图假设，推导了求解当前情况下的格林函数重构方法.
2.3.1 基于路径的格林函数分解定义格林函数G (x, x″0, t)为如下实际非均匀介质中标量波动方程的因果解：
|
(10) |
其中，x″0为震源位置，c=c (x)代表传播速度，ρ=ρ(x)代表非均匀介质的密度，t代表时间，此时G (x, x″0, t)代表位于x″0的点源在接收位置x的声波压力响应.
由于介质的非均匀性，依据如下关系(图 2)，G (x, x″0, t)可以被分解为上传格林函数G－(x, x″0, t)和下传格林函数G+(x, x″0, t)
|
图 2 下传和上传格林函数示意 Fig. 2 Decomposed Green′s functions |
|
(11) |
其中，x表示∂D0下的任何位置，上角标“-”代表传播到x位置时波的传播方向为向上，上角标“+”代表传播到x位置时波的传播方向为向下(如图 3).
|
图 3 上传格林函数和下传格林函数(黑线：下传格林函数；黄线：上传格林函数) Fig. 3 Up-going Green′s function and Down-going Green′s function |
在一维散射问题的Marchenko方程推导过程中，f1(x, t)和f2(x, t)称作Schrödinger方程的基本解.推广到三维情形，定义与f1(x, t)中的相同的X，x′i为∂Di面上的聚焦点，此时f1(x, x′i, t)和f2(x, x″0, t)称为聚焦函数(Wapenaar et al., 2014)，“聚焦函数”很好地反映了其聚焦到地下某点的特性.
与格林函数分解方式类似，f1(x, x′i, t)和f2(x, x″0, t)可以被表示为
|
(12) |
|
(13) |
图 4为聚焦函数的定义示意，聚焦函数f1(x, x′i, t)聚焦在x′i上后继续以发散波f1+(x, x′i, t)的形式向下传播到无反射参考半空间.聚焦函数f2(x, x″0, t)聚焦在x″0上后继续以发散波f2－(x, x″0, t)的形式向上传播到均匀半空间.由于没有损耗波和大角度入射，所以聚焦函数是一个有限带宽的Delta函数.聚焦函数f1+(x0, x′i, t)从地表入射到地下介质，聚焦到x′i的位置，此时f1+(x0, x′i, t)可以表示为从x′i点源激发到地表接收响应的时域反转.
|
图 4 聚焦函数示意 Fig. 4 Focusing functions |
聚焦函数的推导可以由地表接收的反射响应和估计的从地下某点到地表检波器之间的直达波推导得到，对于估计的直达波而言，一个平滑的模型即可.而传统的叠前偏移成像方法对速度模型的精度有很高的要求.
2.3.3 地震超越干涉法完全格林函数重构首先推导格林函数与聚焦函数、反射响应之间的关系.
在深度∂D0、∂Di上，单程聚焦函数之间的关系表示如下:
|
(14) |
|
(15) |
在深度∂Di上的单程格林函数与深度上的聚焦函数f1∓和反射响应RU之间的关系表示如下：
|
(16) |
|
(17) |
依据(14)、(15)、(16)、(17) 式，对式(16)、(17) 左右两侧分别相加，可以得到完全格林函数的表达式：
|
(18) |
假设td(x′i, x″0)为完全格林函数的初至波到达时间，那么t＜td(x′i, x″0)时，G (x′i, x″0, t)值为0，即
|
(19) |
假设f2(x′i, x″0, t)可以被写成直达波和散射尾波的和，则
|
(20) |
其中，f2d(x′i, x″0, t)为f2的直达波，M (x′i, x″0, t)为散射尾波.
又由于f2－(x′i, x″0, t)等于∂D0、∂Di之间的透射波的逆，则
|
(21) |
又由射线路径可知(图 4b)，f2－的直达波f2d－就是的直达波，即
|
(22) |
故
|
(23) |
其中，Tdinv(x′i, x″0, t)为透射波的逆的直达波，它的直达波的达到时间为－td(x′i, x″0)，则有
|
(24) |
假设忽略界面上传播损失(Wapenaar and Berkhout, 1989)，则Tdinv(x′i, x″0, t)可以被近似地写成
|
(25) |
将式(23) 代入完全格林函数表达式可得三维Marchenko方程：
|
(26) |
其中，tcε(x′i, x0)=tcε(x′i, x0)－ε，ε是一个很小的正常量，这可以把直达波包含在积分里.
利用三维Marchenko方程求解散射尾波M (x′i, x″0, t)的迭代算法如下：
|
(27) |
|
(28) |
此时，t＜td(x′i, x″0)
当 t≥td(x′i, x″0)时，M (x′i, x″0, －t)=0.
基于以上公式，如果地表反射响应、估计的直达波已知，那么第k阶散射尾波MK(x′i, x″0, t).此时聚焦函数f2－(x′i, x″0, t)的迭代形式可以表示为
|
(29) |
假设经过k次迭代后迭代算法收敛，那么可以把第k次迭代求取的代入完全格林函数表达式(18)，即可求得G (x′i, x″0, t).
值得注意的是，聚焦函数的求取和格林函数重构基于直达波假设，其假设初至波就是直达波，这种情况并不是总是满足，因为有时初至波是折射波.本文所建立的就是满足此条件的、具有有限弯曲的界面、高阻抗地层、向斜构造的模型.更复杂的构造上进行格林函数重构需要更多的研究(Wapenaar et al., 2014).
由此可知，完全格林函数可以从它的直达波和地表记录的反射响应中重构.
2.3.4 地震超越干涉法上传和下传格林函数重构单程格林函数G－(x′i, x″0, t)和G+(x′i, x″0, t)被用来进行虚拟震源基准面重建和Marchenko成像.根据Slob等(2014)提出的方法，利用上文中G－(x′i, x″0, t)和G+(x′i, x″0, t)的公式(16) 和(17)，可以得到：
|
(30) |
|
(31) |
与上文的分解方式类似，把f1+用∂D0与∂Di之间的透射响应中直达波的逆和散射尾波M表示
|
(32) |
对于尾波M+(x″0, x′i, t)有
|
(33) |
又根据式(14)、式(15) 以及f2(x′i, x″0, t)的因果性可得：
|
(34) |
同样，基于上文提到的“直达波假设”，把f1+(x′i, x″0, t)展开成透射响应直达波的逆与散射尾波的和，代入单程格林函数G－(x′i, x″0, t)和G+(x′i, x″0, t)的表达式(16)、式(17)，可得：
|
(35) |
|
(36) |
求解这两个Marchenko方程的迭代算法如下：
当t＜td(x′i, x″0)
|
(37) |
|
(38) |
其中，
|
(39) |
当t≥td(x′i, x″0)
|
(40) |
此时，单程聚焦函数f1+(x″0, x′i, t)的迭代表达式可以写成：
|
(41) |
当求得f1+(x″0, x′i, t)与f1－(x″0, x′i, t)，代入单程格林函数的表达式(16)、(17) 即可求出上传格林函数和下传格林函数.
2.3.5 基于低频补偿的数据驱动面向目标Marchenko成像地震超越干涉法格林函数重构及Marchenko成像要求地震数据不包含自由表面多次波，所以重建低频缺失的BSI数据的障碍就是BSI与SI重构的数据不一致.低频重构后的宽频带地震数据中包含自由表面多次波、层间多次波、直达波以及一次波，Marchenko成像要求的输入是去除了自由表面多次波的地表接收到的反射响应以及估计的地下虚拟震源到地表检波器之间的直达波，在实际地震数据中，可以采用SRME (Verschuur et al., 1992)进行自由表面多次波消除.本文利用数值模拟的方法进行自由表面多次波及直达波的消除，方法如下：
假设F1(t)是地表存在吸收边界条件(如PML)的反射响应，其包含了一次波、直达波和层间多次波.F2(t)是在均匀介质数值模拟得到的直达波响应.那么自由表面多次波F3(t)可以由下式得出：
|
(42) |
此时符合Marchenko成像的反射响应R0可以表示为
|
(43) |
使用上传格林函数和下传格林函数进行Marchenko成像的控制方程如下：
|
(44) |
此时∂Di是地下任意深度，R0是∂Di下介质的反射响应(Claerbout, 1985; Wapenaar et al., 2008)，R0(x′i, xi, t′)是检波器和震源同时位于∂Di上的反射响应.此时地下介质的像是反射响应R0在t=0且偏移距为0，即R0(xi, xi, t=0).
在频率域，上传格林函数和下传格林函数的互相关函数为
|
(45) |
此时，对零偏移距C (xi, xi, ω)在全频率范围的积分是互相关成像条件.相对应地，在时间域的成像条件为零时刻、零偏移距，即C (xi, xi, t=0).通过计算目标区域每一个独立成像点的反射系数，可以求得目标区域的Marchenko成像结果.所以Marchenko成像是面向目标的成像方法.
如图 5流程所示，基于低频信息补偿的Marchenko成像共分为如下几步：(1) 基于互相关法重构与主动源相同观测系统的被动源地震数据；(2) 补偿去噪之后的被动源数据能量，确定优势频率fcut；(3) 结合低通滤波进行主动源低频与被动源低频匹配滤波；(4) 归一化能量匹配合成宽频带地震数据；(5) 进行数据预处理，利用数值模拟方法消去自由表面多次波及直达波后将其以及估计的地下虚拟震源到地表检波器间的直达波作为输入，并基于Marchenko方程迭代求解格林函数；(6) 基于互相关成像条件，循环求取目标区域的所有成像点位置的反射系数完成Marchenko成像.
|
图 5 基于频率优势的低频重构Marchenko成像流程 Fig. 5 Work flow of low frequency compensation based Marchenko imaging |
我们在模型上(如图 6)模拟接收被动源信号并进行地震干涉法重构.模型横向长6000 m，纵向深2000 m，网格间距均为2.5 m.图 7a中星号覆盖的区域为Marchenko目标成像区域，放大后的目标成像区域如图 7b所示.检波器布设在地表(-3000 m, -3000 m)范围内，检波器间距10 m.为了得到高信噪比的被动源重构结果，模拟地下4500个噪声源随机分布在横向(-2900 m, 2900 m)，纵向(1350 m, 1700 m)范围内，检波器接收1200 s的被动源数据.通过常规互相关重构进行被动源重构，地表(0 m, 0 m)位置的虚拟震源的地震记录如图 8a.
|
图 6 本文所用的速度模型与密度模型 (a)速度模型; (b)密度模型. Fig. 6 Models used in SI and BSI (a) Velocity model; (b) Density model. |
|
图 7 Marchenko面向目标成像区域 (a)目标成像区域在地下介质中的位置；(b)放大后的目标成像区域. Fig. 7 Marchenko target oriented imaging area (a) Location of target imaging area; (b) Zoomed in version of target imaging area. |
|
图 8 地震干涉法重构单炮记录和匹配滤波后结果 (a)互相关重构结果; (b)被动源匹配滤波后结果. Fig. 8 Retrieved SI result and result after Matched filtering (a) SI result by cross-correlation; (b) SI result after Matched filtering. |
首先来模拟低频缺失主动源地震数据，模型同被动源相同.为了在后续更好地重构主动源的低频信息，震源和检波器的位置(观测系统)也与被动源重构出的地震记录一致.观测系统一致性的要求在野外实际地震数据采集中可以被满足.得到的低频缺失观测记录如图 9a所示.主动源正演的时候采用的类似于Sinc波的子波作为震源，该子波在频率域有平的振幅谱.低频缺失观测记录经过上文提出的数据预处理方法后作为Marchenko迭代求解格林函数的输入，同时Marchenko迭代还需要从地下虚拟震源位置到地表检波器之间估计的直达波作为输入.Marchenko不需要精确的速度、密度模型作为输入，本文采用平滑的速度、密度模型(如图 10)进行直达波的估计.重构出的完全格林函数如图 11a所示.利用Marchenko迭代求得的上传和下传格林函数进行Marchenko成像结果.
|
图 9 主动源观测记录与低频重构后宽频带记录 (a)低频缺失主动源观测记录; (b)低频重构后宽频带记录. Fig. 9 Shot records comparison before and after low frequency compensation (a) Low frequency loss shot record; (b) Broadband shot record. |
|
图 10 地震超越干涉法格林函数重构所用的速度和密度模型 (a)平滑的速度模型; (b)平滑的密度模型. Fig. 10 Velocity model and density model used in BSI (a) Smoothed velocity model; (b) Smoothed density model. |
|
图 11 (0 m, 1100 m)处的格林函数对比 (a)主动源数据重构的格林函数；(b)主动源低频补偿后重构的格林函数；(c)参考格林函数. Fig. 11 Green′s functions comparison at location (0 m, 1000 m) (a) Low frequency loss shots Green′s function retrieval; (b) Broadband shots Green′s function retrieval; (c) Reference Green′s function. |
下面基于本文提出的频率优势低频重构方法，利用噪声信号中提取的低频信息来重构主动源地震数据中的低频信息.
首先确定对被动源能量补偿的增益系数，使之在频率域与主动源在低频端的频谱有交点.需要注意的是，此时所做的被动源能量补偿不是为了对主动源进行低频补偿，此步骤仅仅是为了确定优势频率的分界点.依此确定的滤波频率为20 Hz，结合低通滤波器对被动源和主动源单炮分别进行低通滤波.然后，对低通滤波之后的主动源和被动源单炮进行匹配滤波，由于主动源数据中的低频能量较弱，所以用低通滤波主动源数据去匹配低通滤波被动源数据，这样能够保留被动源数据中的有效信息.同时，由于匹配滤波有对信号的选择性过滤作用，所以，被动源信号中的噪声信号在匹配滤波之后的期望输出中不会得到大量保留.
由式(9)，首先对重构的低频数据和主动源数据分别进行归一化，然后根据主动源数据中的归一化之前的最大振幅进行振幅恢复，得到宽频带的地震数据如图 9b.通过对频谱分析图 12分析可以看到，原始主动源中的高频成分得到了很好的保留，同时由低频端的频谱趋势可知低频信号得到了有效重构.在数值模拟时，图 12中所示的振幅差异可以通过调整主被动源的权重(Alali et al., 2016)得到更好的结果，但是在实际的勘探过程中，我们并不能直接模拟全频带的地震记录.抽取单炮记录中的单道波形进行对比(如图 13)可知，基于频率优势低频重构方法重构出的相位与全频带正演很好地符合，振幅有一些差异，这对最终的准确成像结果不会产生较大影响.
|
图 12 全频带正演记录、主动源观测记录与低频重构记录频谱分析 Fig. 12 Frequency spectrum analysis |
|
图 13 主动源观测记录与低频重构记录抽道波形对比 Fig. 13 Single trace comparison between BSI data and BSI data after compensation |
宽频带地震数据经过上文提出的数据预处理方法后作为Marchenko迭代求解格林函数的输入，即在本文的数值模拟的情形下，基于完全相同的观测系统设置，通过正演计算得到直达波、自由表面多次波，然后得到在宽频带数据中消去这两类干扰波后的反射响应.同时采用平滑的速度、密度模型(如图 10)进行直达波的估计.图 11b和图 11c中分别给出了重构的完全格林函数和直接正演得到的真实格林函数(参考格林函数)，很直观地看到二者符合地很好.对重构的格林函数抽道对比(如图 14)，可以看到低频重构之后的格林函数与参考格林函数的波形也符合地很好，验证了本文所提出方法的可靠性.上传格林函数和下传格林函数如图 15所示.利用Marchenko迭代求得的上传和下传格林函数进行低频重构后的Marchenko成像结果如图 16d所示.
|
图 14 (0 m, 1100 m)处的格林函数抽道对比 Fig. 14 Sing trace comparison |
|
图 15 (0 m, 1100 m)处分解的格林函数 (a)上传格林函数；(b)下传格林函数. Fig. 15 Decomposed Green′s functions (a) Up-going Green′s function; (b) Down-going Green′s function. |
|
图 16 Marchenko成像结果 (a)面向目标成像区域；(b)逆时偏移成像结果；(c)基于低频缺失炮集的Marchenko成像结果；(d)基于低频重构后宽频带炮集的Marchenko成像结果. Fig. 16 Marchenko imaging result (a) Target-oriented imaging area; (b) RTM imaging result; (c) Low frequency loss Marchenko imaging; (d) Broadband Marchenko imaging. |
通过对比低频缺失Marchenko成像结果(如图 16c)与低频重构Marchenko成像结果(如图 16d)，我们可以看到，低频重构的反射界面更加聚焦，而在低频缺失Marchenko成像中，由于输入的反射响应里低频缺失，导致旁瓣效应影响了Marchenko成像效果.
为了对比Marchenko成像的效果，本文采用传统的逆时偏移成像(RTM)对Marchenko成像的目标区域进行成像，由于传统的偏移算法基于单次散射假设，要求输入的炮记录不能包含层间多次波.当地震记录中含有多次波时，所得到偏移成像结果中含有虚假的地下反射界面.逆时偏移成像的结果如图所示16b，从图中可以看到，模型中高阻抗地层导致的反射能量较弱导致成像结果的反射界面不够清晰，同时虚假的反射界面出现在成像结果中.
值得注意的是，为了验证匹配滤波对噪声信号的过滤效果，将被动源数据匹配主动源数据，得到了更高信噪比的地震记录如图 8b所示，可见匹配滤波在低频重构中起到了抗噪的作用，最终的低频重构Marchenko成像结果也表明，低频补偿没有在成像结果中引入虚假的反射界面.同时，Broggini等(2014)证明了Marchenko方法具有抗噪性，所以我们认为基于优势频率低频重构的Marchenko成像也具有抗噪性.
4 结论与探讨(1) 传统的地震成像方法对含有高阻抗上覆地层的区域往往由于强散射、高速屏蔽作用导致照明不足、成像效果较差.本文提出了SI与BSI的结合方法SIBSI，通过数值算例证明了SI与BSI结合方法的可行性，本方法同时利用了低频信息较强的高速屏蔽层穿透能力和Marchenko成像算法对层间多次波的消去作用，获得了高精度的地下构造成像结果.本方法所用到的主动源和被动源数据在实际生产过程中可以被采集到，所以实际生产能满足使用该方法所需要的条件.
(2) 相对于纯数值低频延拓，本文的低频来源于真实的地震数据，所以补偿得到的宽频带地震数据更加可靠；相对于现有的基于能量的补偿方法，本文提出的方法能够减弱被动源低频信号中的噪声对重构后宽频带数据的影响.基于归一化处理，主动源地震数据中的高频部分得到了完整的保留，低频数据的能量趋势也与全频带地震数据相符.最终的Marchenko成像结果中没有由于低频噪声而造成的虚假反射界面，且真实的反射界面更加聚焦，消除了由低频缺失导致的旁瓣效应，验证了基于频率优势的低频重构方法的优势和有效性.本文提出的方法也可以看成是基于被动源信息的Marchenko成像算法，因为我们重构出的结果也可以进行Marchenko成像，此时主动源信息看成是对被动源信息的补偿.
(3) Marchenko成像算法依赖于目标成像区域的每一个点的格林函数，很小的网格大小会造成巨大的计算量(指数级).本文在计算效率上进行了如下改进：首先，将炮集数据放入内存虚拟硬盘里，减少大规模并行计算时的I/O并发等待时间；其次，BSI格林函数重构是一个高并行化的过程，耗时主要用于对直达波的估计，由于不同的地下虚拟震源对应着不同的正演过程，所以可以基于MPI或者OpenMP对计算过程进行并行化以加快计算速度.
致谢作者感谢国家自然科学基金项目(41674124、41374115) 和国家高技术研究发展计划(863计划)重大项目课题(2014AA06A605) 联合资助.同时特别感谢来自TU Delft的Jan Thorbecke与来自科罗拉多矿业学院CWP研究组的Alex Jia对本文建设性的帮助.作者最后特别感谢两位匿名审稿人对本文提出的宝贵意见.
| Alali A, Van Der Neut J, Draganov D. 2016. Merging active and passive seismic reflection data with interferometry by multidimensional deconvolution.//86th Ann. Internat Mtg., Soc. Expi. Geophys.. Expanded Abstracts, 5182-5186. | |
| Baysal E, Kosloff D, Sherwood J. 1983. Reverse time migration. Geophysics, 48(11): 1514-1524. DOI:10.1190/1.1441434 | |
| 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 | |
| Broggini F, Wapenaar K, Van Der Neut J, et al. 2014. Data-driven Green's function retrieval and application to imaging with multidimensional deconvolution. Journal of Geophysical Research:Solid Earth, 119(1): 425-441. DOI:10.1002/2013JB010544 | |
| Calvert R W, Bakulin A, Jones T C. 2004. Virtual sources, a new way to remove overburden problems.//66th Ann. Internat Mtg., Soc. Expi. Geophys.. Expanded Abstracts 234. | |
| Claerbout J F. 1985. Imaging the Earth's Interior. Oxford:Blackwell Scientific Publications. | |
| Duan Y Q. 2006. Matched filtering and wavelet shaping. Oil Geophysical Prospecting (in Chinese), 41(2): 156-159. | |
| Han L G, Zhang Y, Han L, et al. 2012. Compressed sensing and sparse inversion based low-frequency information compensation of seismic data. Journal of Jilin University (Earth Science Edition) (in Chinese), 42(S3): 259-264. | |
| Hu Y, Han L G, Xu Z, et al. 2017. Demodulation envelope multi-scale full waveform inversion based on precise seismic source function. Chinese Journal of Geophysics (in Chinese), 60(3): 1088-1105. DOI:10.6038/cjg20170321 | |
| Huang C, Dong L G, Chi B X. 2015. Elastic envelope inversion using multicomponent seismic data with filtered-out low frequencies. Applied Geophysics (in Chinese), 12(3): 362-377. DOI:10.1007/s11770-015-0499-8 | |
| Meles G A, L er K, Ravasi M, et al. 2015. Internal multiple prediction and removal using Marchenko autofocusing and seismic interferometry. Geophysics, 80(1): A7-A11. DOI:10.1190/geo2014-0408.1 | |
| Nakata N, Snieder R, Tsuji T, et al. 2011. Shear wave imaging from traffic noise using seismic interferometry by cross-coherence. Geophysics, 76(6): SA97-SA106. DOI:10.1190/geo2010-0188.1 | |
| Ramírez A C, Weglein A B. 2009. Green's theorem as a comprehensive framework for data reconstruction, regularization, wavefield separation, seismic interferometry, and wavelet estimation:A tutorial. Geophysics, 74(6): W35-W62. DOI:10.1190/1.3237118 | |
| Schuster G T. 2001. Theory of Daylight/Interferometric Imaging-Tutorial.//63rd EAGE Conference & Exhibition, Session A32. | |
| Shen H L, Tian G, Shi Z J. 2013. Partial frequency band match filtering based on high-sensitivity data:method and applications. Applied Geophysics, 10(1): 15-24. DOI:10.1007/s11770-013-0363-7 | |
| Slob E, Wapenaar K, Broggini F, et al. 2014. Seismic reflector imaging using internal multiples with Marchenko-type equations. Geophysics, 79(2): S63-S76. DOI:10.1190/geo2013-0095.1 | |
| Thorbecke J W, Draganov D. 2011. Finite-difference modeling experiments for seismic interferometry. Geophysics, 76(6): H1-H18. DOI:10.1190/geo2010-0039.1 | |
| Verschuur D J, Berkhout A J, Wapenaar C P A. 1992. Adaptive surface-related multiple elimination. Geophysics, 57(9): 1166-1177. DOI:10.1190/1.1443330 | |
| Wapenaar C P A, Berkhout A J. 1989. Elastic wave field extrapolation:Elsevier. | |
| 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, Fokkema J. 2006. Green's function representations for seismic interferometry. Geophysics, 71(4): SI33-SI46. DOI:10.1190/1.2213955 | |
| Wapenaar K, Van Der Neut J, Ruigrok E. 2008. Passive seismic interferometry by multidimensional deconvolution. Geophysics, 73(6): A51-A56. DOI:10.1190/1.2976118 | |
| Wapenaar K, Van Der Neut J, Thorbecke J. 2012. On the relation between seismic interferometry and the simultaneous-source method. Geophysical Prospecting, 60(4): 802-823. DOI:10.1111/gpr.2012.60.issue-4 | |
| Wapenaar K, Slob E, Van Der Neut J, et al. 2013. Three-dimensional Marchenko equation for Green's function retrieval beyond seismic interferometry.//83th Ann. Internat Mtg., Soc. Expi. Geophys.. Expanded Abstracts, 4573-4578. | |
| Wapenaar K, Thorbecke J, Van Der Neut J, et al. 2014. Marchenko imaging. Geophysics, 79(3): WA39-WA57. DOI:10.1190/geo2013-0302.1 | |
| Xu Z. 2012. A study of imaging technique and application based on seismic interferometry method[Ph. D. thesis] (in Chinese). Changchun:Jilin University. | |
| 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 (in Chinese), 50(2): 225-231. | |
| Zhang P, Han L G, Zhou Y, et al. 2015. Passive-source multitaper-spectral method based low-frequency data reconstruction for active seismic sources. Applied Geophysics, 12(4): 585-597. DOI:10.1007/s11770-015-0520-2 | |
| Zhang P, Han L G, Liu Q, et al. 2015. Interpolation of seismic data from active and passive sources and their joint migration imaging. Chinese Journal of Geophysics (in Chinese), 58(5): 1754-1766. DOI:10.6038/cjg20150525 | |
| Zhu H, Wang D L, Shi Z A, et al. 2012. Passive seismic imaging of seismic interferometry. Progress in Geophysics (in Chinese), 27(2): 496-502. DOI:10.6038/j.issn.1004-2903.2012.02.012 | |
| 段云卿. 2006. 匹配滤波与子波整形. 石油地球物理勘探, 41(2): 156–159. | |
| 韩立国, 张莹, 韩利, 等. 2012. 基于压缩感知和稀疏反演的地震数据低频补偿. 吉林大学学报(地球科学版), 42(S3): 259–264. | |
| 胡勇, 韩立国, 许卓, 等. 2017. 基于精确震源函数的解调包络多尺度全波形反演. 地球物理学报, 60(3): 1088–1105. DOI:10.6038/cjg20170321 | |
| 黄超, 董良国, 迟本鑫. 2015. 低频缺失情况下的弹性波包络反演. 应用地球物理, 12(3): 362–377. | |
| 许卓. 2012. 基于地震干涉法的波场成像技术及应用研究[博士论文]. 长春: 吉林大学. http://cdmd.cnki.com.cn/Article/CDMD-10183-1012355017.htm | |
| 叶月明, 姚根顺, 赵昌垒, 等. 2015. 利用地震干涉法衰减海底相关层间多次波. 石油地球物理勘探, 50(2): 225–231. | |
| 张盼, 韩立国, 刘强, 等. 2015. 主动源与被动源地震数据插值及联合数据成像. 地球物理学报, 58(5): 1754–1766. DOI:10.6038/cjg20150525 | |
| 朱恒, 王德利, 时志安, 等. 2012. 地震干涉技术被动源地震成像. 地球物理学进展, 27(2): 496–502. DOI:10.6038/j.issn.1004-2903.2012.02.012 | |