0 引言

利用被动源记录对地下构造进行成像是被动源领域的重要研究点之一.无论是微震事件还是环境噪声，目前常规的对构造进行成像的方法是：先进行主动源记录的重构，再利用重构出的记录进行成像(陶毅等, 2010; 朱恒等, 2012).利用被动源记录重构主动源记录起源于Claerbout(1968)提出的“声波日光成像”，该方法证明了：对由地下传到地表的透射波记录进行自相关，其结果等价于地表自激自收的模拟记录.Schuster(2001)将该方法命名为地震干涉，多位学者对其原理和应用做出了系统的研究(Draganov et al., 2007; Snieder et al., 2009; Wapenaar, 2003; Wapenaar et al., 2004).互相关方法(Schuster and Rickett, 2000; Wapenaar, 2004)是地震干涉中最常用的方法，它可以用于处理单个微震记录，也可以用于处理噪声记录以获得有效信息(Bakulin and Calvert, 2006; Schuster et al., 2004).随后发展起来的反褶积方法(Vasconcelos and Snieder, 2008a, b)以及多维反褶积方法(Ravasi et al., 2015; van Dalen et al., 2014, 2015; Wapenaar et al., 2008, 2011)分别解决了互相关方法的某些缺陷.

目前所有基于波动方程的偏移方法都是以逆散射为理论基础发展而来的(Bleistein et al., 2001)，这其中包括了单程波偏移(Claerbout and Johnson, 2010; 张关泉, 1986; Zhang et al., 2003)、逆时偏移等(Deng and McMechan, 2007; Fletcher et al., 2009; Sun et al., 2006).地震数据的散射波传播包括入射波传播、散射、散射波传播三个过程.这一类偏移成像方法需要先重构出入射波与散射波的传播过程，再根据“时间一致性”成像原理(Claerbout, 1985)，来得到散射位置处的构造成像.广义上讲，反射波也属于散射波.入射场与散射场的重构过程是一个数值求解微分方程(即波动方程)的过程，其初始值一般取波场值为零，边界条件就是在某些位置的入射波和散射波的记录.在常规的主动源勘探当中，这两者都容易得到：入射波的记录可用估计出的震源子波，位置在震源激发处；散射波的记录可用去除直达波的共炮点道集，位置在检波点处.但在被动源勘探中，寻找满足散射关系的入射波与散射波记录并不容易，因此这是利用被动源记录对构造进行成像的最关键的问题.

另外，这里还需要注意区分：在“重构主动源记录”和“重构入射场与散射场”中，两个重构所指的是完全不同的过程.

在先重构主动源记录再成像的方法中，被动源成像被转化成为了常规的主动源成像问题，由此解决了重构入射场和散射场的边界条件问题.然而在此情况下：一方面，入射场与散射场重构的边界条件的位置都是地表；另一方面，地震干涉重构实际上是提取了被动源记录中的一部分自由表面多次波(Wapenaar et al., 2010)，而忽略了透射波所携带的来自震源下方的反射波信息.因此，深部构造的成像效果不如浅部好.

主动源记录的重构要求震源包围观测系统，这通常是个不易达到的条件.对被动源记录进行直接偏移成像则无此限制，但这方面的相关研究目前还很少.Artman (2006)提出了将微震记录既当作入射场重构的边界条件，也当作散射场重构的边界条件，从而实现了直接偏移成像.该方法可以不用分离出单个的微震事件，但其利用的有效信息实际上仍然是自由表面多次波.这一点与常规地震勘探中的自由表面多次波成像(Guitton, 2002; Muijs et al., 2007; Liu et al., 2011; Lu et al., 2011; Tu and Herrmann, 2018)是一致的，即：在进行波场正传时，作为震源而进行传播的数据并非震源子波，而是包含反射波和自由表面多次波的观测数据.

从微震的角度出发，考虑到震源位置可以利用微震记录和速度模型确定(Artman et al., 2010; 葛奇鑫等, 2017)，震源子波也可以通过反演等方法求得(胡勇等, 2017).于是，在本文中，我们假设已知微震事件的震源位置及其子波，并将偏移成像过程中的入射场的激发位置移动到了地下真实震源处.这样，在对震源下方的构造进行成像时，有效信息并非由自由表面多次波带来的散射波，而是由震源波场直接与深层介质作用而产生的散射波.因此该方法有能力得到深部构造较为清晰的成像.这样所得的成像结果中，除了低频噪声外还存在一种主要干扰：由震源发出的上行波在与上方的介质作用后产生的多级散射波(包括层间多次波等).它与有效信息耦合严重，无法在记录中直接分离.我们分析了其在传播过程中的特点，利用左右行波分离对该干扰进行了压制.对于子波未知的情况，我们采用了镜像法与激发时间成像条件(丁仁伟等, 2008)相结合，进行成像值提取.最后，利用数值实验验证了本文方法的有效性.

1 方法原理 1.1 常规的微震数据直接成像方法

常规的微震记录直接成像的原理可用图 1来说明.以一个二维非均匀介质为例，地下存在一个微震震源S，在符号U_N=i^-/+中：上角标的减号表示波场上行，加号表示波场下行；下角标中的N表示经过自由表面反射的次数.z=0表示地表的位置.地表处的波场U_N=0^-即是检波点接收到的所有未经自由表面反射的波场.经过一次自由表面反射的下行波场U_N=1⁺在与地下介质相互作用后，可以在地表接收到其对应的散射波场U_N=1^-.这样，U_N=1⁺和U_N=1^-就是一对可以用来重构入射场和散射场的边界条件，散射发生的位置，即成像目标，为地表以下所有的散射体.

U_N=1⁺与U_N=0^-相比，除了方向成镜像关系外，其余均相同，故U_N=1⁺可以用U_N=0^-来代替.以此类推，所有的U_N=n^-和U_N=n+1^-都可以构成一对用于重构入射场和散射场的边界条件，且成像目标也完全相同，理论上都是地表以下所有的散射体.由此可知：除了U_N=0^-只能用于重构入射场外，其余U_N=n^-均可同时用于重构入射场和散射场.将所有对结合在一起，可知：用于重构入射场的边界条件就是全体U_N=n^-的集合，即时间轴正向的、完整的微震记录；而用于重构散射场的边界条件则是时间轴倒转的、去除U_N=0^-之后的微震记录.它们的位置都是在地表检波点处.但是，在微震记录中，U_N=0^-和其余U_N=n^-是耦合在一起的，在尾波较长时二者难以分离，故在该方法中，U_N=0^-并未被分离或压制.

综上所述，在常规的微震数据直接成像方法中，入射场和散射场的重构使用的是相同的数据，即微震记录，只是二者时间轴的方向是相反的.接下来选择逆时偏移，并使用互相关成像条件，则地下某点x_i的成像值I(x_i)可以表示为

(1)

式中，U_s为微震记录正向传播的波场值，U_r为微震记录逆时反传的波场值；x_r和x′_r表示波场激发的位置(即检波点处)；t为波场传播时间，T表示微震记录总时间.

该方法在选择入射场与散射场重构的边界条件时，着手点是自由表面的反射，它提取了自由表面多次波作为有效信息以实现偏移成像.其基本思想可以简单表述为：地表接收到的所有波场都会以镜像方向再次向地下传播.由于该方法不需要有关震源的任何信息，故无需对微震记录中的单个微震事件进行分离.即使多个微震事件在记录上耦合在一起，该方法仍然有效.然而，由于其波场重构的边界条件均位于地表，故对于深层介质的成像效果欠佳；另外，相对于常规主动源勘探中的震源子波，该方法的入射场过于复杂，成像结果中可能会因此出现高频噪声.

1.2 震源下方构造成像

微震事件发生在地下，它能够携带很多深层信息而被地表检波点接收.微震成像理应对深层信息有所体现.关于微震定位的研究由来已久，目前已比较成熟.关于微震子波估计的研究相对较少，但通过反演的方法也可以较为准确地确定，或者也可以通过某些反褶积方法去除记录中的子波影响.因此，我们假设微震的位置和子波都是已知的，并将其当作入射场重构的边界条件.自然地，我们也希望将地表的微震记录当作散射场重构的边界条件.然而，在此情况下，微震记录中的成分与上一节会有所不同，成像的目标区域也因此需要有所改变.

需要注意：在常见的波的类型中，直达波不属于散射波，不可用于散射场重构；单散射是基于波动方程的偏移成像方法的假设前提(Oristaglio, 1989)，单散射场是成像所用到的主要有效信息，也是散射场重构的必要信息；多级散射波，在单散射假设前提下，同样会被当作单散射波来进行处理，因而会造成假象，故不适合单独用来进行散射场重构.

基于此，在本方法中，微震记录的成分可大致分为四类.如图 2所示，震源到地表之间存在散射点A和B，以及平缓的反射界面C和D，震源下方也是非均匀介质，但其构造形态没有具体给出.U₁为直达波，U₂表示震源下行波与更深层介质作用而产生的散射波，U₃表示震源上行波与散射点作用而产生的散射波，U₄和U₅用来代表各类多级散射波，也包括自由表面多次波.

结合单散射假设，可知，当以真实微震事件为入射场重构的边界条件时：对于震源与地面之间的介质，微震记录只包含来自散射点(尺度远小于波长)的单散射波，因而完全无法对平缓的界面进行成像；但是，对于震源下方的介质，微震记录则包含了来自所有类型的散射体的单散射波，因此可以用于散射场的重构.同样选择逆时偏移，并使用互相关成像条件，则震源下方某点x_i的成像值I(x_i)可以表示为

(2)

式中，U_s为微震事件正向模拟的波场值，U_r为微震记录逆时反传的波场值；x_s表示震源位置，x_r表示检波点位置；t为波场传播时间，T表示微震记录总时间.

至此，针对震源下方的介质，我们构造出了相应的用于入射场和散射场重构的边界条件，实现了利用微震数据对目标散射体进行直接偏移成像的目的.

1.3 尾波干扰压制

图 2所示的这些波场在被地表检波点接收时并无绝对的到达先后顺序，除直达波一般最先到达且较容易识别外，若速度结构较复杂，其余波场基本都是耦合在一起的.而在本方法中我们认为只有U₂是有效信息，其他波场(由震源上行波产生的各类波场)均为干扰.在常规方法中，微震记录中的U_N=0^-在重构散射场时，是残留的入射直达波(相对其入射波而言)，属于干扰场.但由于其对成像结果影响并不大，而且压制困难，故并未被处理.同样的，在本文方法中，我们首先需要分析这些干扰场对成像过程的影响方式及程度，以确定是否有必要对这些干扰进行处理，再对影响较大的干扰进行压制.

在本节中，我们将这些干扰波简单分为两类：直达波以及由震源上行波产生的直达波后的尾波.如图 3a所示，U_s1为有效信息，即来自震源下方的散射波，U_s2代表各类尾波干扰，而直达波并未绘出.需要注意到：在对散射场的传播过程进行重构时，微震记录逆时反传的最终时刻相当于微震事件的激发时刻，此时，逆时反传的直达波恰好聚焦于震源点.由此可知：直达波在逆时反传过程中只存在于震源及其上方(震源和检波点之间)，不会对成像目标区域造成影响；但是，其尾波则会完全穿过震源而到达成像目标区域，因此会造成较为严重的干扰.如图 3b所示，U_s表示入射波，即震源下行波，U_r1和U_r2表示逆时反传的检波点波场.其中U_r1是我们所需要的散射波，U_r2为尾波干扰，它穿过震源到达了成像目标区域.直达波仍然没有绘出.

在原始的微震记录中，由于直达波一般最先被接收到且较容易识别，故其分离相对容易.但若要在记录中对U_s1与U_s2 (或者说U_r1与U_r2)直接进行分离则是非常困难的，因为其耦合程度很高，且在几何形态、振幅、相位等属性上并无明显差异.对此，观察它们的传播方向，发现：在震源上方，U_r2的传播方向大致与直达波相同，是向震源点汇聚的；穿过震源到达震源下方后，U_r1和U_r2的传播方向是不同的.在震源左下方，入射波U_s左行，检波点波场中的有效信息U_r1右行而干扰波U_r2左行；而在震源右下方则变为U_s右行，U_r1左行而U_r2右行.将(2)式中的U_s和U_r分离成左行和右行部分，可得：

(3)

式中U的下角标中的第二个字母，L表示左行波，R表示右行波.在分解所得的四项之中，根据图 3b以及相关的分析可知，第一项和第四项(由不同向波场所得的成像值)是对成像结果有贡献的成分，其余两项即为尾波干扰.基于此，对重构出的入射场和散射场进行左右行波分离，提取出有效信息，则能够有效压制尾波干扰.

与上下行波分离(徐兴荣等, 2012)类似，以x(水平轴)-z(深度轴)-t(时间轴)数据体为例，左右行波的分离可以在对t轴和x轴进行二维傅里叶变换后的ω(频率)-k_x(水平波数)域实现.分离的方程可表示为

(4)

(5)

(6)

(7)

简单来讲，(4)-(7)这四个方程说明：在ω-k_x域，一三象限对应的是左行波，而二四象限对应的是右行波.再将其反变换回时空(t-x)域，即可完成左右行波的分离.

在去除了由同向波场所得的成像值之后，最终得到的成像结果为

(8)

图 4给出了一个数值算例来展示压制尾波干扰的效果.一个微震震源位于Marmousi模型的1100 m深度的位置(图 4a)，地下微震事件采用30 Hz的雷克子波进行模拟.可以明显看出，相比原始成像结果(图 4b)，在利用左右行波分离压制尾波干扰后(图 4c)，高频噪声得到了有效的压制，同相轴连续性也得到了加强，总体上来讲，整个成像结果更加清晰.但是，也可以观察到：零偏移距附近的成像效果受到了削弱.这是因为零偏移距附近的波接近竖直传播，在左右行波分离中会有所丢失.另外，震源及其附近的成像值过大且十分复杂，该处的构造基本被掩盖而不可识别.这是因为：在微震刚刚激发的时刻，图 3中的U_s1集中在震源周围，而U_r1也聚焦在震源处.因此，在取最终成像结果时应从震源下方一定距离开始(向下).

另外，注意到：在图 4b、c中，震源上方的成像结果也体现出了构造的形态特征.这一部分成像值的来源大致有两种：

第一，与传统方法的来源相似.在散射场重构的流程上，本文方法与传统方法完全相同，都是将微震记录逆时反传.而对于入射场重构，本文方法模拟了完整的微震事件，自然也包括在地表(自由表面)发生的反射.在偏移速度模型足够精确的情况下，入射场重构过程中的地表处的记录应该与真实的微震记录基本相同，这样本文方法的成像结果就完全包含了传统方法的结果.即使在一般情况下(指偏移速度模型不够精确)，其地表处的记录也包含了真实微震记录的部分有效信息(主要是图 1的U_N=0^-中的震源直达波)，因此可以得到震源上方的构造信息.

第二，由直达波这种不满足散射关系的波场“成像”得来.震源上行直达波在重构的入射场和散射场当中都是存在的，它们之间不满足散射关系，不能用来成像，但却会依赖于偏移速度模型而产生相应的类似构造成像结果的假象.这种假象完全取决于偏移速度模型，不可当作有效信息.因此，后文中，在重构散射场时，首先对直达波进行了切除.

2 数值模拟算例

本节利用二维声波数值模拟实验对本文方法的效果进行了验证，并主要针对微震事件本身的两个要素将其与传统的方法进行了对比：一是微震事件的数量(或者说疏密程度)；二是微震记录的耦合程度(单个微震的记录是否可分离).

若无特殊指定，本节所用的模型、观测系统等各类参数如下：正演速度模型以及偏移速度模型如图 5a所示，为Marmousi模型中界面较为平缓的一部分；模型大小为230×500，网格间距为10 m，250个检波点以20 m为间距均匀分布在地表，时间采样间隔为1 ms，地下微震事件的子波均采用30 Hz的雷克子波，并且将微震的位置、激发时间与子波均设置为已知量.

由于本文方法仅将震源下方的成像结果视为有效信息，故本章所有算例的震源位置均设定在地下浅部，具体为深度不超过700 m，以增加成像目标区域的范围；相应地，在对比成像结果时，只展示深度在800 m以上的部分，这也是为了避开如图 4所示的震源附近的强干扰.其对应的速度模型如图 5b所示.

前文所述的尾波干扰是本文方法所特有的一种噪声，而除此之外，逆时偏移中最常见也最严重的噪声还属低频噪声.对此，我们使用了上下行波分离与拉普拉斯滤波相结合的策略来对该噪声进行压制.

2.1 震源稀疏，记录可分离

震源包围观测系统是利用被动源记录重构主动源记录的假设前提之一，稀疏的震源分布不满足此条件，是该方法无法处理的一种情况，也正是直接成像方法的优势之一.因此首先对震源稀疏的情况做了测试.在本节中，单个微震记录被设定为可分离的，即对每一个微震事件的记录进行单独成像，再进行叠加，由此得到最终的成像结果.单个微震记录可分离也是主动源记录重构的另一条假设前提.

为了保证在最终叠加时，水平方向不同位置的覆盖次数是基本一致的，在本算例中，假设在地下500 m深度的位置等间隔分布有9个微震震源(图 6a).由该模型得到9份单事件微震记录，对每份记录分别成像再进行叠加，传统方法的结果如图 6b所示，本文方法的结果如图 6c所示.单从结果来看，两图对于构造界面的刻画都是比较连续而精细的，且没有出现明显的错误；但对比来看，尤其是图像左下方这一片区域，图 6c明显比图 6b要更加完整和清晰，这在一定程度上体现出了本文方法对深层界面进行成像的能力.

图 6d展示的是未对尾波干扰进行压制的成像结果.与该结果相比，在对该噪声进行压制之后的结果(图 6c)中，杂乱的高频噪声得到了衰减，同相轴的连续性与清晰度均得到了提高.观察到图 6d本身所受的干扰并不强，故由尾波干扰压制所带来的成像质量的提高也并不十分明显.

2.2 震源稀疏，记录不可分离

微震属于瞬态源，虽然不同微震事件的震源子波在形态上有些许差别，但它们之间始终具有很强的相干性.在主动源记录重构中，之所以要假设单个微震记录可分离，就是为了避免由这种相干性带来的串扰噪声.这种噪声具体表现为重构记录中的大量假同相轴.

上一节的算例是满足这一假设的，而在本节中，为了测试这条假设对于直接成像方法的影响，仍然使用图 6a所示的震源分布模型，但将微震事件设定为同时激发，即所有微震事件的记录都耦合在一起，无法分离，仅使用该记录进行一次成像.

在此情况下，所得的成像结果如图 7所示.图 7a-c的意义分别对应于图 6b-d.从不同假设条件的对比来看，相比于图 6、图 7中的各结果在成像质量上均有较大幅度的下降，具体表现在：出现了大量的高频噪声和假同相轴，这些都可以理解为串扰噪声；而且，同相轴的连续性和清晰度也都有所下降.从不同方法的对比来看，相比于传统方法的结果(图 7a)，在本文方法的结果(图 7b)中：高频噪声较少，图像左下方区域同相轴的连续性和清晰度仍然略优；然而，其右侧陡倾角界面的成像效果则较差.从尾波干扰压制的效果来看，相比于未压制该噪声之前的结果(图 7c)，图 7b的成像质量有了明显的提升.

串扰噪声的成因可做如下简单解释：从地震干涉的角度来讲，对于某一散射界面，其入射波和散射波具有因果关系，是一一对应的；在成像过程中，这样的一对入射-散射波在散射界面处发生相长干涉，从而得到较大的成像值，此为真实同相轴；但若成像过程中同时存在较多的震源，由于微震信号的相干性，许多不具有因果关系的入射波和散射波在非散射界面处也能够发生相长干涉，从而造成假同相轴，即串扰噪声.

2.3 震源稠密，记录不可分离

本节的实验在上一节的基础上，将震源数量由10增加到300左右，模型如图 8a所示，震源随机分布在其中的白色虚线方框中，其余设定均保持不变.在此情况下，所得的成像结果如图 8b-d所示，其意义分别对应于图 7a-c.可以看出：在震源稠密时，两种成像方法的对比以及尾波干扰压制的效果与上一节基本一致；但相比于图 7a-c，在本节的成像结果中，假同相轴的连续性和延伸性变差了，即这些假象更加接近于随机噪声，这既在一定程度上突出了真实界面信息，也有利于后续的噪声压制.

这种现象可以通过以下几点分析来解释.成像结果中的真实同相轴的位置与观测系统，或者说震源和检波点的位置，是无关的，它总是出现在地下真实的散射界面处；但假同相轴的位置是与观测系统有关联的，它在成像结果中并没有固定的分布规律.尤其是对于微震信号，其震源不仅横向位置不同，纵向位置也是不同的，更增强了这种无序性.也就是说：震源越多，真实同相轴处的叠加效应越强，有效信息也越突出；而假同相轴并不会出现叠加效应.因此随震源数量的增加，假同相轴分布的随机性会增强，从而逐渐向白噪声趋近.由于子波的影响，这种噪声的频谱应当是带限的.

2.4 震源子波未知

从地震记录中获取有效的震源子波信息通常不是一个容易的过程，尤其是在实际资料中，估计出的子波精度不高.震源定位不需要对子波进行估计，但子波误差会在逆时偏移的成像结果中造成相位误差.对此，在无法获得准确子波信息时，我们提出在成像时可以使用激发时间成像条件，而非互相关成像条件.激发时间成像条件只需要重构散射场，不需要重构入射场，即无需其波形信息，而是利用入射场的走时来提取散射场中的有效信息.

计算入射场的走时要求已知微震的位置和激发时间，因此也是以微震定位为前提的.以图 9a所示的模型为简单示例，地下浅层存在一个震源(白色圆点).按照本文的成像方法，入射场重构的边界条件是近真实的微震事件，于是，此时目标区域(震源下方)的入射场走时就是以真实震源为震源点计算得到的初至走时(图 9b).

另外，由图 4的相关解释可知，在使用互相关成像条件时，本文方法实际上近似包含了传统方法的成像结果.这一点在激发时间成像条件中也可以实现.传统方法的入射场(图 1)十分复杂，并非单一同相轴的波场，无法计算出其中所有同相轴的走时；然而，在这所有的同相轴中，由震源上行直达波经一次自由表面反射得到的反射波是其入射场中能量最强、也是对成像结果贡献最多的成分.这条同相轴的走时可以利用简单的镜像法来求得(图 9c).

使用2.1节的模型和参数设定，假设震源子波未知，采用激发时间成像条件.对每一份单事件微震记录，从重构的散射场中将图 9b、c所代表的两类有效信息提取出来并叠加，再对9份成像结果进行叠加，所得结果如图 10所示.图 10与图 6c的成像质量没有明显差别，在同相轴连续性和界面清晰度上基本一致.

3 讨论与结论

利用微震记录对地下构造进行成像的传统方法是先重构主动源记录，再利用主动源成像方法进行成像.与该方法相比，常规的微震记录直接成像方法主要解决了重构法中的以下几点缺陷：

(1) 主动源记录的重构要求震源密集且包围观测系统，而直接成像方法则无此限制，零散微震事件的记录也可以用来成像.

(2) 所谓重构并非创造信息，而是提取有效信息.这两种方法所利用的有效信息都来源于自由表面多次波.重构法提取出的有效信息一般只是自由表面多次波中的一小部分，而直接成像方法则能够对所有的自由表面多次波进行有效利用.

(3) 重构法需要对单个微震事件的记录进行分离，以避免由微震信号相干性所导致的串扰噪声；而直接成像法则不需要这一步(Artman, 2006).然而，本文的数值实验显示：在直接成像法中，若不对单个微震事件的记录进行分离，则串扰噪声也会出现.随着微震数量的增加，串扰噪声在成像结果中的分布会向带限的随机噪声趋近.

本文在已知微震定位结果的情况下，对常规的微震记录直接成像方法进行了修改：将入射场的重构方式修改为近真实的微震事件，将成像的目标区域修改为震源下方.这种修改所带来的优缺点如下：

(1) 入射场重构的边界条件的位置更接近深部介质，重构出的震源下方的入射场因而更加精确.不同条件下的数值实验均显示该方法对深部界面的成像效果相比常规方法有所提高.

(2) 震源上行直达波后的尾波在重构出的散射场中成为干扰成分，在成像结果中表现为高频噪声.通过左右行波分离可以压制该噪声，但数值实验显示，这种压制方法也会在一定程度上降低陡倾角界面的成像效果.

(3) 重构出的入射场近似包含了常规方法的入射场的主要信息，所以，由此得到的成像结果实际上接近于两种方法的结合.

(4) 除了入射场和散射场的重构，微震定位还需要消耗额外的计算资源，这限制了本文方法的应用.然而，除了互相关成像条件，本文方法还可以使用激发时间成像条件来进行成像，这一点是常规方法所不能做到的.激发时间成像条件不需要重构入射场，节省了大量的计算资源.