在被动源地震勘探中，由水力压裂、钻井或者断层位移等因素引起的微震，对了解地下构造具有重要的意义，其时空坐标是最需要关注的信息.微震定位主要包括两大类方法：基于波动方程的逆时反传定位成像方法(time-reverse imaging, TRI)以及基于程函方程的走时反演方法.

走时反演根本上属于最优化问题，它将震源坐标设定为未知参量，以某种时差为目标函数.最早的经典走时反演方法(Geiger, 1912) 利用模拟初至和实际初至走时之差，建立了最小二乘目标函数.双差定位方法(Waldhauser and Ellsworth, 2000) 建立目标函数的基础是两个微震事件的模拟初至与实际初至走时之差的差值.这些方法在刚提出时均采用了局部最优化算法，随着计算机技术的发展，全局最优化算法，比如遗传算法(宋维琪和杨晓东, 2011)、模拟退火算法(Pei et al., 2009)、粒子群算法(Deep et al., 2012)等，逐渐在走时反演中扮演了越来越重要的角色.

然而，初至拾取一直是此类方法存在的最大问题，同时也是最值得研究的方面(宋维琪和冯超, 2013;谭玉阳等, 2016).事实上，微震数据通常都具有体积大、信噪比低的特点(王晨龙等, 2013)，这限制了传统走时反演方法的应用.在参考了逆时偏移的原理和思路后，TRI方法由此得到发展(McMechan, 2007; Artman et al, 2010).该方法无需对初至时间或相位进行拾取，并且对低信噪比有着良好的稳定性.

TRI的基本思想可以简单描述为：波场的传播在数学意义上是可逆的(即逆时不变性).据此，将微震记录时间轴倒转，以检波点为震源将记录反传入速度模型，波场会在震源的空间位置和激发时刻聚焦(即相长干涉)(Claerbout, 1985).这种思想与逆时偏移是一致的，故逆时偏移中的成像条件也被引入了此领域.但由于记录中只有检波点波场，且在成像过程中只存在反传波场，因此，TRI中的成像条件描述的并非正传波场与反传波场之间的计算关系，而是来自相同或不同检波点的反传波场之间的计算关系.在逆时偏移中作为一个整体的反传波场，被拆成分散的部分来处理，这使得成像条件的应用更加灵活，比如算术平均成像条件(Gajewski and Tessmer, 2005)、几何平均成像条件(Nakata et al., 2016)、混合成像条件(Sun et al., 2015)等.

走时反演通过迭代更新来求取的目标参量就是震源坐标本身，所以在反演结束后无需后续的处理；而TRI方法在获得成像结果之前都不需要人工干预，但最终需要通过观察成像值的大小和分布来判断震源的存在以及位置，因此，成像结果中的假象，或者不同震源间成像值的差异，都会直接影响到震源识别的精度和正确性.能够给TRI结果带来假象或错误的因素有很多，比如低信噪比(Zhu, 2014; Li et al., 2015)、震源低主频、震源低能量、震源大深度、高速体屏蔽、小尺度速度扰动(李振春等, 2014) 以及速度模型误差(Zhang et al., 2015)等.针对假象问题，Sava和Poliannikov (2008) 提出一种干涉成像条件，能够压制由稀疏观测带来的串扰噪声和假象；Witten和Artman (2011) 提出在数据域中估计噪声模型，并结合阈值压制假象.针对不同震源间成像值的差异问题，考虑大量出现于全局最优化算法中的随机策略，其突出优点之一为：能够获得在常规方法中易被掩盖的信息，比如在众多局部极小值中寻找到全局极小值(马昌凤, 2010).类似地，在TRI中，需要将不易定位的震源(弱震源)从易定位的震源(强震源)的掩盖中分辨出来，但这方面目前仍鲜有研究.

本文从不同成像条件的结合与对比出发，并讨论了高通滤波在加强成像效果方面的作用；然后将随机策略引入了反传检波点的选择方法中，使得水平方向不同位置的聚焦强度出现差异，能够将不易定位的震源(弱震源)从易定位的震源(强震源)中分辨出来；其次，提出了根据波场反传时间序列来度量震源存在的可能性的方法，将震源判断空间由简单的成像域扩展到成像-时间域，以此提高震源识别的精度与正确性；最后，针对成像结果的加强以及假象的消除，利用声波模拟数据验证了该方法在多种因素干扰下的可行性与稳定性.

在逆时偏移中，每一个成像点都会经历一个波场正传序列和一个波场反传序列，而在微震定位中，由于震源的未知性，成像过程中的正传波场是无法准确模拟的.因此，成像条件在由逆时偏移引入TRI时需要做一定的调整，可以表述为：在某一成像点，来自相同或不同检波点的反传波场序列之间的计算处理方法.对于两个序列的对应(相同)时刻的波场值，三种基本的处理方法包括：加、乘以及平方，扩展到整个序列(沿时间轴叠加)则分别对应于：直接叠加、零延迟互相关以及零延迟自相关.

对于多个序列在同一时刻的波场值，这三种基本方法可表示为

(1)

(2)

(3)

式中，I_l(x, t)(l=1, 2, 3)代表在空间位置 x 与反传时刻t处的临时成像值(temporary imaging value, TIV)，R_i(x, t)代表在(x, t)处来自第i个检波点的波场值.

这三种基本方法的结合也存在多种形式(Sun et al., 2015).若将检波点分组，则组内和组间可以采用不同的处理方法，现给出三种混合形式如下：

(4)

(5)

(6)

式中，G为分组的数量，N(j)为第j组中检波点的数量.在(4) 式中使用连乘是为了拉大震源点与非震源点之间TIV的量级差距，从图像角度讲，即是为了提高对比度；此外，这也是为了拉大震源激发时刻与其他时刻之间TIV的量级差距.(5) 式中的平方能够继续加强这两个效果，而且移除了TIV中的负值，使得结果更易处理.但(5) 式只是(4) 式的一种简单扩展，其指数并不是固定的，而且该指数不能过大，否则计算机的计算过程会很不稳定.(6) 式将(5) 式中的连乘改换为连加，在这三者中拥有最小的TIV量级差与对比度.

图 1展示了同一份模拟数据在同一个成像点分别利用(4) 式和(6) 式所得到的TIV序列.二者在大致形态上(图 1(a、c))基本相同，但从细节上(图 1(b、d))看明显是前者更优秀，因为TIV的峰值越清晰，震源存在的判断标准就会越简单.

图 1

Fig. 1

图 1 (a)、(c)相同数据、相同空间位置分别使用(4) 式和(6) 式作为成像条件得到的TIV序列；(b)、(d)图(a)和图(c)中峰值附近的局部放大 Fig. 1 (a)、(c) The time sequences at a source point with the same gather and using equation (4) and (6) as the imaging conditions respectively. (b)、(d) Zoom in the peak of (a) and (c) respectively

(1) 式到(6) 式中的I(x, t)只是某一时刻的成像值(TIV)，要求得成像域完整的成像值(whole imaging value, WIV)，还需要将其沿时间轴叠加：

(7)

即使使用了效果较好的成像条件，对于由单次成像过程得到的结果而言，如果没有其他信息作为对照或参考，其可靠性也可能会受到怀疑；另外，若同时对许多震源进行定位，难定位的震源容易被易定位的震源所掩盖.为此提出了一种反传检波点随机选择的策略：首先，在总体中选择随机数量、随机位置的检波点(称其为选取过程)；之后，将所选择的检波点进行随机分组(称其为分配过程).对这两步处理进行重复(假设重复D次)，每一次新的分组都需要一次额外的传统定位成像过程.在(4) 到(6) 式中，组内波场值的连加与实际物理过程中的波场叠加是等价的，故组内所有检波点同时反传即可得到连加结果；但组间的连乘却需要进行多次反传(组数减一，即(G－1)).这样，反传模拟的总次数将会达到(G－1)·D.考虑到效率，G和D的取值应该谨慎.

图 2给出了一次分组的简单例子.总共100个检波点，选择了其中的96个(选取过程)，并将其随机分为4组(分配过程)，直线上的黑点代表检波点的位置.显然，每一个检波点在全部组中至多出现一次.

图 2

Fig. 2

图 2 一次分组的简单实例共100个检波点，选择了其中的96个，并将其随机分为4组(四条线)，直线上的黑点代表检波点的位置. Fig. 2 A simple example of random receiver grouping at a time 96 of total 100 receivers are selected and divided into 4 groups randomly; black dots on the lines denote receivers′ positions.

不同的分组造成检波点分布出现疏密差距，使得反传波场在水平方向上对地下的聚焦出现强弱之分.当某个震源受到强聚焦，则该震源更易被分辨，该方法也因此有能力突出弱震源.又因为随机选择时的均匀概率分布，故在聚焦弱的位置也会有些许信息冗余.它们能够补偿不同次分组之间的信息缺失并提高成像的可靠性.

图 3展示了使用相同数据、不同随机分组方法所得到的四份定位结果.震源分布如图 3e所示，使用均匀模型，网格点距10 m，检波点以20 m的点距均匀分布在地面0~7000 m的范围.成像条件为(4) 式，为了显示效果，对最终成像值WIV取对数.明显可以看出，每一份结果都有不同的强聚焦突出区域：(a)两侧；(b)中间；(c)右侧；(d)左侧.

图 3

Fig. 3

图 3 使用相同数据、不同随机分组方法所得到的四份定位结果 每一份结果都有不同的强聚焦突出区域：(a)两侧；(b)中间；(c)右侧；(d)左侧; (e)震源真实位置. Fig. 3 Four location results with different grouping of a same synthetic gather using equation (4) as the imaging condition The velocity model is homogenous and position is the only distinction between these sources. Each result has its own emphasis: (a) Two sides; (b) Middle; (c) Right; (d) Left; (e) The true sources distribution.

震源处的成像值WIV通常会呈现出峰值形态，且比非震源位置的成像值要大，这些峰值即为判断震源存在和位置的依据.相对于全模型的成像结果图，将这些峰值看做其中的高频成分，可以使用简单的高通滤波来提取、加强这些高频成分.若成像结果中出现假象，则震源的判断会受到严重的干扰，滤波也可能会加强这些干扰.注意到WIV实际上是由TIV序列沿时间轴叠加得到的，本节将对于震源的判断由简单的WIV扩展到TIV序列，以此提高识别的精度与正确性.

图 4a和4d为同一个模拟实验中两个网格点的完整TIV序列，其中(a)对应某一震源点，(d)对应某一非震源点，使用(4) 式作为成像条件.在震源点的TIV序列(图 4a)中，存在唯一一个数量级远大于其他位置的峰值；而在非震源点的TIV序列(图 4d)中，存在众多数量级相同的峰值.宏观来看，二者最大值所在的峰并无明显差异；而在此峰以外，前者基本在零值线附近徘徊，明显要比后者“干净”.基于此现象，希望度量TIV序列相对于零值的离散程度，即序列的“干净”程度，由此判断该点是否可能存在震源.因此，若震源点与非震源点的离散程度差异越大，同时震源点的TIV序列越“干净”，则判断标准越容易确定.

图 4

Fig. 4

图 4 (a)某一震源点的完整TIV序列；(b)对(a)进行局部放大；(c)去除(a)序列中最大值所在的峰；(d)、(e)、(f)与(a)、(b)、(c)意义相同但对应于某一非震源点 Fig. 4 (a) A whole time sequence of a source. (b) Zoom in the peak of (a). (c) Remove the peak of (a). (d)、(e)、(f) have the same meanings but refer to a non-source point

为此，将最大值所在峰去除(以零值代替)(图 4c和4f)，则震源点TIV序列几乎不再存在偏离零值的点，而非震源点TIV序列仍然十分复杂.经此步骤，二者的离散程度下降了大致相同的数值，但在比值上的差距却极大地增加了，有利于震源存在性的判断.另外，如图 4b和4e所示，此峰的宽度并不一定相等，多出一个偏离零值的点对非震源点TIV影响较小，但对震源点TIV影响较大，此为需要去除该峰的另一个原因.

基于上述思想，建立了两个函数，来大致度量TIV序列相对于零值以及峰值所在时刻的离散程度：

(8)

(9)

式中I(x, t)表示移除峰之后的TIV序列，K为TIV峰值对应时刻.(8) 式和(9) 式可以在波场反传的过程中进行计算，无需将TIV序列(大体积的三维数据体)存储在计算机内存中.

在计算出模型各网格点的离散程度之后，结合经验阈值(震源处与非震源处的离散值基本不在一个数量级上，因此阈值比较好确定)来构造筛选模型(例如图 9d).此筛选模型是由1和0组成的，离散程度小于阈值则赋值为1，代表可能存在震源，保留该点WIV；离散程度大于阈值则赋值为0，代表不可能存在震源，去除该点WIV.筛选模型也可以理解为特殊的权值分配，权值仅由1和0组成.

对于一次定位成像，由(8) 式和(9) 式可以构建出两个筛选模型，将其进行结合，结合的方式为：两个筛选模型同时显示该点可能存在震源，则保留，否则移除；再考虑到2.2节随机方法所述的多次定位成像，各次的结果也需要进行结合，结合方式同上，以此来消除偶然性的影响.

为了验证本文方法的效果及其影响因素，进行了多组声波数值模拟实验.若无特殊指定，实验参数设定如下(有特殊指定则以该实验具体设定为准)：模型采用1043×302的Marmousi模型或相同大小、速度为3000 m·s^-1的均匀模型，网格间距10 m，时间采样间隔1 ms，使用一种改进有限差分法(Liu and Sen, 2009) 来进行正演以及反传模拟，记录总时间8000 ms，地表均匀分布522个检波点，点距20 m，定位成像过程共重复3次(即2.2节中的参数D)，选取过程中选择的检波点为522的80%到90%，分配过程中的分组数量期望为5，使用(4) 式作为成像条件.

在本节的三个算例中，均有十个震源分布在均匀模型中.前两个算例中各自的十个震源深度分别相同(图 5a，图 6a)，前十个主频从左到右逐渐由25 Hz提高到35 Hz，后十个最大振幅从左到右逐渐由1提高到10；第三个算例中的十个震源沿着一条斜线分布(图 7a)，子波完全相同.不同震源的激发时间是分散的.

图 5

Fig. 5

图 5 (a)十个震源分布在相同深度，主频从左到右递增；(b)初步定位成像结果沿震源深度(1510 m)的截面；(c)经过滤波的截面；(d)应用筛选模型后的截面 Fig. 5 (a) Ten sources are located at a same depth whose DFs are progressively increasing from left to right; (b) The cross section of the preliminary result at the common depth (1510 m); (c) The cross section of the filtered result; (d) The cross section after sifting

图 6

Fig. 6

图 6 (a)十个震源分布在相同深度，最大振幅从左到右递增；(b)初步定位成像结果沿震源深度(1510 m)的截面；(c)经过滤波的截面；(d)应用筛选模型后的截面 Fig. 6 (a) Ten sources are located at a same depth whose maximum amplitudes are progressively increasing from left to right; (b) The cross section of the preliminary result at the common depth (1510 m); (c) The cross section of the filtered result; (d) The cross section after sifting

图 5b为第一个算例初步成像结果(未经过滤波和筛选)在震源所在深度的截面.可以看出，WIV与震源主频之间成正相关.这表明，当对拥有不同主频的多个震源同时进行定位时，低主频的震源更难得到好的结果.图 5b中最左侧的WIV峰值几乎等于最右侧的背景值，这致使其很容易被识别为假象.这种现象主要源于低频信号的低分辨率特性.为了平衡不同WIV峰值间的差异，使用高通滤波来压制低频成分.由图 5c可见，弱震源的成像值得到了很好的补偿.图 5d为应用筛选模型后的截面.

图 6为第二个算例的模型与结果.从结果中可观察到其规律与前一个算例基本相同：在初步成像结果中，弱震源处的成像值基本等于强震源处的背景值(图 6a)；采用高通滤波使弱震源的成像值得到了很好的补偿(图 6b)；应用筛选模型基本将非震源处的成像值归零了(图 6c).

震源深度的影响需要从三个方面来考虑.一方面，若模型较复杂，则震源深度越大，波场传播的路径会越长且越复杂.考虑到波场反传时同样存在反射过程致使波场能量减弱，因此，深度越大，定位应当越困难.第二方面，震源深度越小，远偏移距接收到的波场的主频就会越低.根据前文主频的影响方式可知，震源深度越大，定位应当越容易.第三方面，在检波点分布一定的前提下，震源深度的增加会减小其相对检波点分布范围的张角，即减小了观测孔径，这会降低定位成像的分辨率.

此处将震源深度的第一种影响归类于模型因素，只对后两种影响在均匀模型中做了测试.图 7b为第二个算例的初步成像结果沿着震源所在斜线的截面，不同峰值之间的差距并不明显，使用滤波能够将这种差距进一步抹平.图 7d为应用筛选模型后的截面.由此可知，与模型复杂程度相比，纯粹的深度对定位成像的影响并不大.

图 7

Fig. 7

图 7 (a)十个震源沿一条斜线分布，主频相同；(b)初步定位成像结果沿震源所在斜线的截面；(c)经过滤波的截面；(d)应用筛选模型后的截面 Fig. 7 (a) Ten sources are located along an oblique line whose DFs are the same; (b) The cross section of the preliminary result along the line containing all the ten sources; (c) The cross section of the filtered result; (d) The cross section after sifting.

筛选模型根据微震能量在时间轴上的集中性来去除假象，首先测试其对于镜像效应这种常见假象的消除效果.图 8a所示为300×300、网格间距10 m的二层速度模型，上层速度2000 m·s^-1，下层速度4000 m·s^-1，单一震源位于(150, 200) 处，地表均匀分布150个检波点，点距20 m.为了成像显示效果，对WIV取对数，再将过小(小于最大值的50%~60%)的WIV对数值充零(下同).未加筛选模型时(图 8b)，定位成像结果中除了真实震源以外，在模型上部还出现了一处非震源聚焦现象.这是因为反传波场在界面产生的反射波同样会聚焦到一点，若该界面反射系数较大，则可以得到与真实震源强度相近甚至更强的聚焦.在这种简单模型情况下，筛选模型能够较为完美地去除这种假象(图 8c).

图 8

Fig. 8

图 8 (a)二层介质速度模型，上层速度2000 m·s-1，下层速度4000 m·s-1，单一震源位于(150, 200) 处，以黑点指示；(b)未经筛选的定位成像结果；(c)应用筛选模型后的定位成像结果 Fig. 8 (a) A two-layer velocity model, with its upper velocity 2000 m·s-1 and lower velocity 4000 m·s-1. A single source is placed at (150, 200). (b) Location result without sifting. (c) Location result after sifting

本节第二个算例测试了一个较为复杂的综合情况，共有十个震源随机分布在Marmousi模型中(图 9a)，其参数由表 1给出.由于最大振幅与主频的影响相似(3.1节)，故只设定了主频变化.为各震源设置不同的主频和深度是为了提高定位结果对于定位方法的灵敏度.如果拥有最低主频和最大深度(这里代表最复杂的传播路径)的震源能够被正确定位和识别，那其他震源也可以.如果某几个震源无法正确定位，也能够较为容易地分析出震源丢失的原因.

表 1(Table 1)

表 1 本算例中的震源参数 Table 1 Arguments of the sources in the experiment

表 1 本算例中的震源参数 Table 1 Arguments of the sources in the experiment

图 9b和9c分别展示了未经筛选以及经过筛选的结果.两份结果均能对全部十个震源进行正确且清晰的定位成像，可见复杂模型本身对镜像效应有一定的抵抗力，但只有在使用了筛选模型之后，图 9b中由黑色箭头所指示的假象才会消失(图 9c).这体现了该方法在消除假象上的能力.图 9d给出了此算例中的筛选模型示意图，黑色表示1，白色表示0.

图 9

Fig. 9

图 9 (a)十个震源(黑点)分布在Marmousi模型中；(b)未经筛选的结果，黑色箭头指示假象；(c)应用筛选模型后的结果；(d)本算例中求得的筛选模型，黑色表示1，白色表示0 Fig. 9 (a) Ten sources (black dots) distributed in Marmousi model; (b) Location result without sifting. The artifacts are pointed with black an arrowhead; (c) Location result after sifting. Artifacts disappeared; (d) Sifting model of this experiment. Black denotes 1; white denotes 0

利用定位出的震源位置以及对应点TIV序列中的峰值时刻((9) 式中的K值)，可以确定微震的激发时间(没有展示).另外，在下文中的结果图展示中，将会更多地直接使用经过筛选的结果.

仅从表达式来看，(4)、(5)、(6) 三式的计算稳定性排序为(6) > (4) > (5)，灵敏度或对比度的排序为(5) > (4) > (6).图 8a和8b分别展示了采用(6) 式和(5) 式作为成像条件的定位成像结果.本算例使用了与上一节相同的数据.这里所谓的计算稳定性是计算过程的一种反映，它与收敛性无关，只是表示数值的指数应当被限制在计算机字长的取值范围内.

当震源分布较为简单时，这种计算稳定性通常不会出现问题.(6) 式只包含一次乘法运算，即平方，拥有最小的对比度，这导致了图 10a中左下角两个震源的丢失.而(5) 式拥有最多的乘积运算以及最大的对比度，与图 9c相比，此结果的细节甚至要更好一些，比如右下角震源的聚焦范围更小、更精确了.

图 10

Fig. 10

图 10 (a)、(b)分别采用(6) 式和(5) 式作为成像条件的定位成像结果 Fig. 10 (a)、(b) The imaging results with equation (6) and (5) respectively

利用3.2节的模拟数据，在其中加入高斯白噪声以测试本方法的抗噪能力.首先，使噪声基本掩盖住弱同相轴(图 11b).此算例的结果(图 11c)与图 9c相比基本没有任何质量降低.接下来引入更强的高斯白噪声以基本掩盖住所有可见同相轴(图 11d).由图 11e可以观察到深部震源出现的轻微焦散，而假象控制依然出色.这体现了本方法在白噪声条件下的较好稳定性.理论上，对于纯粹的白噪声，由于其互不相干的性质，经成像条件计算所得结果(WIV)基本为零值附近的随机噪声，对于有效信息的干扰非常小.但大地滤波会使白噪声也具有子波特性，不相干性质受到破坏，由此影响成像精度.

图 11

Fig. 11

图 11 (a)截取的部分原始地震记录；(b)、(d)原始记录添加高斯白噪声；(c)、(e)分别使用数据(b)和(d)进行定位成像的结果 Fig. 11 (a) The noise-free record in the first experiment; (b)、(d) Records with strong Gaussian white noise. Gain is applied for display; (c)、(e) Location result using the record (b) and (d) respectively

实际上地下介质并非平滑的，其中会存在许多小尺度的扰动.将这种扰动引入Marmousi模型中(图 10a)来获得更接近现实的正演数据.最大扰动量设为10%.由图 12c可以识别出全部的十个震源，但也出现了少量的假象.然而，与图 12b相比，筛选模型对于假象的消除效果依然很明显.

图 12

Fig. 12

图 12 (a)向Marmousi模型引入至多10%的速度扰动；(b)未经筛选的定位成像结果；(c)应用筛选模型后的定位成像结果 Fig. 12 (a) Introducing at most 10% velocity perturbations into a Marmousi model; (b) Location result without sifting; (c) Location result after sifting

对于TRI类方法，波场的逆时不变性需要以精确的速度模型为前提，但即使通过FWI等方法获得了速度模型，该模型也很难与真实信息完全吻合.此算例中，对Marmousi模型应用了5个网格点尺度的平滑(图 13a)，用以进行波场反传.所得结果(图 13c)的假象控制相比于图 13b来讲仍然较为出色.然而，表 1中的第三个震源(坐标(337, 261))在图 13b和13c中均丢失了.显然，这并非筛选模型的问题.注意到该震源的最低的主频(25 Hz)以及较大的深度(2610 m)，可以认为是该算例给出了本方法的一个临界情况.若将准确的模型看做背景模型与扰动模型两部分的叠加，则平滑后的模型会在一定程度上失去其扰动部分.在波场反传时，这对微震直达波的影响较小，但对层间多次波的正确归位影响较大，故模型越复杂(即层间多次波越丰富)，模型平滑对于定位精度的影响越大.

图 13

Fig. 13

图 13 (a)对Marmousi模型进行5个网格点尺度的平滑；(b)未经筛选的定位成像结果；(c)应用筛选模型后的定位成像结果 Fig. 13 (a) Applying a 5-grid-scale smoothing to Marmousi model; (b) Location result without sifting; (c) Location result after sifting

过多地考虑计算效率问题限制了分配过程中分组的数量，而且如此数量的检波点(522) 确实过于理想化，本节放宽了这些限制.仍然采用3.2节的数据，本节前两个算例对参数分别做了以下改动：将分组数量期望增加到10(检波点数不变)；减少约一半的选取检波点数(分组数量不变).

图 14a和14b分别展示了将分组数量期望增加到10以及减少近一半检波点数的结果，这两份结果的质量与图 9c基本相同，尤其是右下方震源的聚焦状况.由图 14a可见，在检波点数量充分的情况下，相比于增加分组数量，使用更合适的成像条件(图 10b)对于定位质量的提高更有效.

图 14

Fig. 14

图 14 定位结果：(a)分组数量期望为10；(b)检波点选取数量为总体的40%到50%；(c)检波点选取数量为总体的10%到20% Fig. 14 Location results with (a) the group expectation of 10; (b) 40%~50% selected receivers; (c) 10%~20% selected receivers

通常来说，检波点越多，定位结果的质量就会越高.但此处存在一个事实，即不同震源的定位存在难易之分，这一点前文已有所讨论.诚然，成像所需的信息都来自检波点反传，更多的检波点可以带来更多的成像信息.但选取的检波点越多，强震源与弱震源之间用于成像的信息的差距也就会越大，因此弱震源就越有可能被掩盖.这里所谓的高质量，其对象只是强震源.由图 14b对比图 9c可知，在检波点数量充分的情况下，定位质量与检波点数量之间并没有绝对的正相关关系.

理论上，在此条件(检波点数量充分)下，从检波点的角度来讲，观测孔径(或者说震源点相对检波点分布范围的张角)的大小才是影响成像分辨率的主要因素.一般来讲，张角越大，分辨率越高，且沿张角平分线方向分辨率最低，垂直于张角平分线方向分辨率最高.因此，由地面观测数据所得定位成像结果的纵向分辨率较低，这一点在本文各结果中普遍存在.

然而，如果选取的检波点数量持续下降，至总体(522) 的10%到20%，由图 14c的右下方可以观察到明显的焦散现象，此时检波点数量充分的条件在此算例中不再满足.

本节将讨论的重点转向震源群，分别对震源群的时间和空间分布的疏密程度进行了对比实验.选用均匀模型以及相同的震源主频(30 Hz)，以去除它们的干扰.

第一个与第二个算例分别测试了较少数量(图 15b)以及较多数量(图 16b)的、激发时间与空间位置均稀疏的微震事件.两份结果均未出现任何的假象与焦散.即使将所有震源设定为同时激发(图 16c)，即激发时间稠密，定位质量也完全没有降低.然而，注意到在图 16c的黑色圆圈中，成像稍显模糊，图 16b相同位置的成像值也非常弱，这与震源时空分布同时较为稠密有关.

图 15

Fig. 15

图 15 (a)均匀模型中分布有14个震源，激发时间分散；(b)由(a)的模拟数据得到的定位结果 Fig. 15 (a) 14 sources distributed in homogeneous media with discrete origin times; (b) Location result of (a)

图 16

Fig. 16

图 16 (a)均匀模型中分布有45个震源；(b)由(a)激发时间分散的模拟数据得到的定位结果；(c)由(a)激发时间相同的模拟数据得到的定位结果 Fig. 16 (a) 45 sources distributed in homogeneous media; (b) Location result of (a) with discrete origin times; (c) Location result of (a) with the same origin time

图 17a展示了大量震源存在于一个小区域内的情况(只展示了存在震源的部分区域)，即空间分布稠密，而激发时间仍是稀疏的.图 17b为其定位结果，总体上看，震源位置基本是正确的，但聚焦情况相互差别较大.可见单纯较大的震源空间密度仅对成像聚焦情况有较大影响，对定位的正确性并无太大影响.造成这一点的主要原因可以归于以下几点：成像过程会使波场在某真实震源A处产生最强聚焦，而在震源附近也会产生一些弱聚焦，其成像值的数量级比真实震源A略小；如果在A附近的弱聚焦带中存在另外的某震源B，则对于同时存在于A与B两个震源的弱聚焦带中的点，其TIV序列就可能出现来自A、B两个震源的峰值；归一化使得其中较大的峰值等于1，较小的峰值略小于1；去除较大的峰后，较小的峰的离散度依然较大，从而在筛选模型计算中被去除；而对于某个相对孤立的震源C，其周围的弱聚焦带在应用筛选模型时同样会被保留，因此从最终结果上看，定位精度差别较大；然而，弱聚焦带TIV序列的峰值一般远小于震源点，故对于震源点，其干扰在很大程度上可以忽略，而定位准确性依然有所保证.

图 17

Fig. 17

图 17 (a)均匀模型小区域内分布有45个震源；(b)由(a)激发时间分散的模拟数据得到的定位结果；(c)由(a)激发时间相同的模拟数据得到的定位结果 Fig. 17 (a) 48 sources distributed in a small-scale district of homogeneous media; (b) Location result of (a) with discrete origin times; (c) Location result of (a) with the same origin time

现在考虑描述微震群时空分布的三个要素：震源数量、空间聚集程度以及激发时间聚集程度.它们根本上影响的是地震记录的复杂程度，进一步讲，即记录中同相轴的重叠程度.震源数量对定位质量的影响是间接的，后两者是直接影响因素.

根据逆时不变性，理论上无论模型多复杂，反传波场都会在真实震源处实现最强聚焦.但是若震源的时空聚集程度处在一个较高水平(图 17c)，那些来自不同微震事件的记录在反传时，多条同相轴的交叉重叠会在某些非震源点也产生相长干涉，致使其成像值甚至可能大于震源点.这些点在计算TIV序列离散度时会出现混乱，从而降低筛选模型消除假象的能力.本节检波点只分布在地表，故此类假象基本只分布在震源群的上下方.震源群内部此类假象出现较少.如果这种干扰相长干涉发生在震源点，则震源点的TIV序列会变复杂，也会导致筛选模型的错误判断.当震源的时空分布较为稀疏时，以上两种效应非常弱，可以忽略.

由本节的算例可知，震源的时间与空间聚集程度的关系可以描述为“与”(“与或非”中的与).单一的时间或空间高聚集程度并不会对定位结果造成太大影响，而当二者均处在高水平时，定位结果的质量则会急剧下降，以至于只能分辨出震源群的大体轮廓(图 17c).

井下微震监测数据信噪比通常较高，并且初至清晰可见.为此，本节测试了使用井中模拟观测数据进行微震定位的效果.考虑到井中检波器分布范围(通常采用10~20个井下检波器，检波器间距为10~20 m)的限制，用其对几千米范围内的大量微震进行同时定位并不现实(观测孔径过小).故本节采用了网格间距为1 m、速度为2500 m·s^-1的均匀模型，以对更小范围内的微震进行探讨.考虑到检波器的位置应尽可能地靠近微震发生区域，并在深度范围上涵盖该区域，设计了如图 18a所示的实验.其中A与B指示了两口监测井，在0 m到500 m的深度范围内各分布有点距为20 m的25个检波点.设置两口监测井是为了降低非震源聚焦的影响(王晨龙等, 2013).在两井间某一区域分布有29个震源，微震最大振幅保持相同且激发时间离散.主频由30 Hz提高到90 Hz，这除了井中信号本身频率较高的原因外，还因为对于该空间尺度，30 Hz子波的分辨率过低，在发生微震较少的情况下，尚能从定位成像结果中分辨各震源的位置.但对于较为密集的微震，即使采用高分辨率的成像条件(比如(5) 式)，本文方法也很难得到震源的精确位置.在以上条件下，可得到如图 18b所示的较为精确的定位结果.

图 18

Fig. 18

图 18 (a)模型以及微震分布，A与B指示了两口监测井，深度均为500 m；(b)由(a)的模拟数据得到的定位结果 Fig. 18 (a) Velocity model and sources distribution. The straight lines A and B denote two monitoring wells, whose depths are both 500 m; (b) Location result of (a)

波导效应源自高速体对波场传播的限制.在图 19a中，震源群左右两侧分别有一块近三角形和近四边形的高速体.这两部分震源在定位结果中均完全丢失了，如图 19b.其根本原因可以追溯到震源激发后的波导效应，在地震记录以及成像过程中很难对此进行补偿.

图 19

Fig. 19

图 19 (a)在Marmousi模型的一片包含高速体的小区域内分布有82个震源，所有震源主频相同，激发时间分散；(b)由(a)的模拟数据得到的定位结果 Fig. 19 (a) 82 sources distributed in a high velocity block contained district of the Marmousi model with the same DF and discrete origin times; (b) Location result of (a)

除了TRI类方法以及前文提及过的走时反演方法之外，基于偏移的定位方法(Gajewski et al., 2007) 以及基于波形反演的定位方法(Zhang et al., 2015) 近年来也均有所发展.

这些方法在实现步骤上互有差异，但在思想上有着紧密的联系.这四类方法所包含的基本要素无非以下几种：波场，几何，成像以及反演.TRI是波场与成像的结合，走时反演包括几何与反演，基于偏移的方法包括波场、几何与成像，基于波形反演的方法包括波场与反演.

反演需要利用最优化算法不断地进行迭代更新，属于耗时较长的一种策略.但是如果每次迭代需要计算的仅为走时信息的话，由于走时计算耗时较小，其整体效率并不会太低.根本上，在这四个要素中，用以获得波场值的数值模拟过程(有限差分等)才是耗时最严重的.因此，包括本文方法在内的多数TRI类方法，在计算机性能不够强的前提下，都不太适合用于微震的实时监测.

另外，这四类方法都存在固有的缺陷，这是由其理论基础决定的.这些缺陷一般表现在定位过程中，那些最易出现问题、最需要人工干预的步骤有，例如TRI中的假象消除、震源判断，走时反演中的走时拾取等.

针对这些缺陷，学者们提出了许多自适应、自动化以及优化方法，但基本都只是在一定程度上解决了某些问题，完全抛弃人脑的判断能力目前仍是不现实的.本文的方法也是基于此，希望能将“一定程度”进一步深化.但该方法与在实际资料处理中取得良好的效果之间还有一定的距离：

(1) 需要先完成向弹性波的扩展.构建筛选模型根本上依据的是：微震属于瞬态源，其能量在时间轴上是集中的.而实际资料中存在的横波若统一按照纵波处理，则无法在反传时间轴上实现唯一聚焦；

(2) TRI类方法普遍对记录白噪声具有较强的抵抗力，这也是此类方法提出的初衷之一.但对于其他影响因素，若无针对性的策略，此类方法并无明显抵抗力.本文方法对某些影响因素进行了针对压制，但仍有若干因素未被纳入讨论范围；

(3) 前文在讨论影响定位质量的因素时，我们大多采用了控制变量法，即每次只对某一个因素进行改变并加以测试，其余条件均保持理想化，而实际资料处理面临的则是所有影响因素的综合.这种综合，尤其是那些未被针对压制的影响因素的综合，会使定位质量大打折扣；

(4) 对于主动源数据的偏移成像，其成像结果的质量除了与成像方法有关外，还在很大程度上依赖于原始地震记录的处理质量.同样地，TRI归根结底只是一种成像方法，其结果也依赖于输入数据的质量.然而，对于被动源数据的处理方法(去噪、补偿、校正等)，除了向主动源勘探借鉴外，由于其极低的信噪比，具有针对性的研究目前仍鲜有报道.

本文提出的微震定位方法属于逆时反传定位成像，包含两个主要关键点，即重复的反传检波点随机选择，以及筛选模型.前者实现了水平方向不同位置的聚焦强度差异化，可以有效突出弱震源的成像值，并带来适当的信息冗余以提高定位的可靠性；后者从地下各点的波场反传序列出发，为震源的判断标准增加了时间维度，以此提高了震源识别的正确性与精度.此外，一些次要的策略，例如成像条件结合、高通滤波等，也在一定程度上加强了定位成像的效果.

理论计算表明，检波点随机选择与筛选模型的结合对于速度扰动、速度平滑、以及记录白噪等干扰因素具有一定的抵抗性；且在检波点数量不充分的条件下，其定位成像质量依然可以保持一个较高的水平.一般来讲，TRI类方法的计算消耗相比走时类方法要更大，但若计算资源较充裕，此类方法对于实时监控的意义仍是不可忽视的.

然而，对于时空聚集度过大的微震事件，其较差的分辨能力依然是令人头疼的问题；而向弹性波以及实际资料处理的扩展也是接下来所需要研究的主要内容之一.