0 引言

地震定位技术广泛应用于目前的天然地震工作中，包括工程地震、地震预警、微震监测在内的许多领域都需要探知震源的大地坐标及深度信息.近年来，国内外学者不断地对定位方法进行改进和创新，出现了很多前沿定位技术(田玥和成晓非，2002；杨文东和金星，2005).微震定位需要从地震数据中提取出部分已知条件，通常是初至拾取时间以及地震台网布设数据等等(谭玉阳等，2017).

最初在1912年，德国物理学家Geiger采用最小二乘方法求解非线性方程组，第一次采用系统的方法实现了定位这一概念，现代生产技术使用的定位方法也大都起源于此.随着20世纪计算机技术的快速发展，定位方法得到了多次改进.如Powell法、双差定位法以及现在主要使用的基于计算机技术的智能化自动定位方法都形成了各自的算法体系(Lienert et al., 1986).进入21世纪后，生产行业的蓬勃发展促进了相应理论的深入研究，地震干涉成像技术逐步开始用于震源定位(Schuster et al., 1997; Schuster et al., 2003; Schuster et al., 2004).

本文在三维模型下对一类定位方法进行了相关研究(Waldhauser and Ellsworth, 2000; Yu and Schuster, 2006).主要思想是根据光学干涉的原理，将光波替换为数值模拟采用的声波，进而完成一系列的干涉成像步骤(王璐琛等，2016).首先对地震数据进行互相关处理，计算出成像点到检波器组合的时间差值.根据成像条件，随时间差累加后的结果在震源处会明显大于其他区域.这一计算过程略去了常规方法所需要的初至时间，适用于压裂范围内若干震源延时激发且信噪比偏低的微震定位.在数值模拟过程中逐步验证了这一方法在三维模型中的可行性(Nelson and Vidale, 1990)，并对误差进行了分析(Aki K and Lee, 1976; Aki K et al., 1977).最后根据一组三维的微震数据定位出了多个震源位置，证实了这一定位方法在微震定位工作中的有效性(Prugger and Gendzwill, 1988; Daneshvar et al., 1995; Gharti et al., 2010).

1 方法原理

图 1为二维模型下提取时间差值的示意图(Claerbout, 1968).已知反演后地下介质速度模型以及地面检波器排列位置，对地下成像范围进行网格化，计算每一网格点到不同检波器的直达波传播时间, 根据传播时间差值取互相关结果的成像值(Campillo and Paul, 2003; Gouédard et al., 2008).网格点靠近真实震源时，计算得到的时间差值趋近实际时间, 即取各互相关曲线峰值, 累加后得到最大的成像值.基于这种成像理论提出了干涉定位方法.

1.1 二维定位方法

在二维模型中，地面处横向设置检波器排列，接收地震数据.把A、B两道记录进行互相关运算，得到结果φ(A, B)，公式为

(1)

其中F[]和F^-1[]分别是傅里叶变换和傅里叶逆变换.由公式(1)进行所有数据的互相关运算.

计算地下介质中一点到各道的直达波传播时间τ，得到任意两个检波器组合的时间差值τ_i(Spence，1980)，公式为

(2)

根据时间差值在互相关表格中进行取值，得到成像值.所有成像值累加后得到地下一点处的总成像值M(x, z)(Vasconcelos and Snieder, 2008a，b)，公式为

(3)

震源附近计算得到的M(x, z)会明显高于外侧区域，即前文所述的成像值.进行归一化后，可提取成像值最大网格点的坐标为定位结果.

1.2 三维定位方法

在地面设置多条测线，进行所有检波器数据的互相关运算，同时计算地下各点到各检波器的直达波传播时间，整理出所有组合的时间差.根据时间差在所有组合的互相关运算结果中进行累加，得到的M(x, y, z)即为干涉成像值.与上述原理相同，震源处的累加取值点趋近各列数据的峰值位置，累加后的结果是最大值，公式为

(4)

不同于上一节二维成像理论中公式(3)计算出的纵向剖面结果，公式(4)可以得到三维模型中水平方向的定位结果.已知水平位置后，再结合二维定位理论，进行纵向剖面的定位就能够得到三维空间中的水平位置及深度信息.

2 数值算例 2.1 模型试算

图 2为一个二维三层速度模型.布设范围为2000 m×2000 m的计算区域.测线的道间距为50 m，共31道检波器，第一道检波器位于x方向250 m处.采样时间为0~2000 ms, 采样率为1 ms.设置三层P波速度分别为2000 m/s、2200 m/s、2500 m/s.震源设置在(1000 m, 1800 m)处，震源子波选取Ricker子波，主频为50 Hz.

根据干涉成像的方法, 首先对31道数据进行互相关运算，得到以时间差值为纵坐标的数据.再计算不同检波器的时间差值.因相邻两层之间存在速度差异和密度差异，波的传播受到波阻抗的影响发生透射，传播方向有所改变，所以在第二层利用二分法，在第三层利用射线追踪法确定直达波的传播路径(赵爱华等，2015)，由此得到地下每一网格点与不同检波器组合的时间差值.根据公式(3)对所有网格点进行成像，结果如图 3所示.

可以看出M值较大的位置集中在(1000 m, 1800 m)附近，呈椭圆形范围.观察累加值较大的区域，在x方向较为集中，但是在深度方向的拉伸变形较为明显.横向定位的结果范围较小，相比之下深度定位的范围较大，导致成像区域形状是椭圆形，而不是圆形.

将震源附近区域网格化，以网格大小为1 m进行成像，得到定位坐标为(995 m, 1806 m).定位结果存在误差，但是在可接受范围之内.

进行三维模型中的数值模拟.布设大小为2000 m×2000 m×1000 m的三维计算区域，在地面设置20条检波器测线.采样时间为0~2000 ms, 采样率为0.5 ms.为减少时间差值引起的误差，取2000 m/s的单一速度模型.设置震源在(1000 m, 1000 m, 500 m)处，选取主频为20 Hz的Ricker子波，并加入随机噪声.

根据所有测线的地震数据进行互相关运算.图 4a为其中一条测线的地震数据.图 4b为第10道数据与其他道的互相关运算结果，其中各列数据的峰值成一圆弧状.由上述定位原理可知，若根据震源位置累加后的结果数值达到最大值，则得到的时间差值取值点应该贴近各道的峰值位置.

将三维模型网格化，得到全区域成像结果.如图 5所示，抽取其中500 m深度以及x=1000 m、y=1000 m处的三张切片，交点附近为定位结果的范围.

由于测线均布设在地面上，所以水平方向的成像结果受到两个方向数据的约束，收敛性较好.而纵向成像时相当于只有一个方向数据的约束，只在横向方向上收敛，所以会产生图中椭球形的成像结果.若加入钻井数据会改善这一现象(He et al., 2007; Hornby and Yu, 2007).

设置网格间距为1 m，分析定位结果的误差.图 6为震源附近成像结果的俯视图和侧视图.其中定位坐标为(997 m, 997 m, 506 m).可知水平方向绝对误差为3 m，纵向绝对误差为6 m.由于纵向集中度较为不好，在表 1中列出了震源在不同深度时的纵向误差(蔡明军等，2004).

2.2 微震定位

本文采用了三维空间中的一组微震数据来测试这种干涉成像方法的定位效果(Pei et al., 2009).井下发生水力压裂事件，裂缝随时间不断扩张，激发了压裂范围内的若干震源.随着压裂的进行，若干震源在不同时间激发，且随机分布于压裂中心点附近的一定范围内.已知三维速度模型，设置分布于地下某一范围内的14个震源，位置随机.地面设置网状单分量检波器81个，接收微震数据.

根据该定位原理，同时对14处震源的数据进行运算定位，得到网格化后的定位结果.如图 7所示，红色球体为微地震震源真实位置，蓝色球体为震源定位结果，所有微震均在不同的随机时刻激发.图 8为定位区域的俯视成像图，可以看出这种成像方法对较近的震源具有一定分辨率，当震源间距离小于20 m时，二者附近的较大成像点会相互贴近，且在震源集中范围内存在高于远区成像值的背景值.由于根据速度模型计算得到的时间存在误差，导致定位结果也有一定的误差，但可以观察到定位结果收敛于真实微震的区域，证实了该成像方法可以对小范围多震源同时进行定位，适用于微震定位工作.

表 2是14处微震的定位结果，表中给出了三个方向的误差.结合图 8可以发现在保证传播时间的精度后，由于多个震源之间距离较近导致定位坐标有靠近的趋势.说明在多震源的微震定位时，较近震源之间也会产生相互影响，这一现象是除时间精度外导致误差存在的另一主要因素.观察到引起的误差在10 m量级上，在可接受范围之内.

3 结论 3.1

本文基于光学干涉原理提出了三维震源成像方法，通过数值模拟和微震定位验证了这种成像方法的可行性，算法的已知条件需要地下介质反演后的速度模型和测线数据，根据互相关运算结果以及时间差值表格可得到成像结果.文中通过数值模拟逐步将成像理论在三维模型中进行实现，经过成像方法的几个步骤，最后得到了可靠的震源定位坐标，并通过误差分析验证了成像方法的稳定性.对一组三维模型下的微震数据进行应用，也得到了多个震源的定位结果.

3.2

从两例三维数据的验证可以看出，该成像方法适用于有随机噪声的地震数据，也适用于压裂范围内若干震源随机时刻激发的情况.但也存在一定的问题，对大区域直接进行高精度的震源成像会产生较大的计算成本.可以考虑先进行分辨率较低的成像，再根据结果对单震源或多震源再进行小范围的高精度成像，提高计算过程的有效性.解决纵向收敛性的问题需要结合钻井数据.而定位误差主要由时间差值的计算精度不够引起，需要反演后更加精确的速度模型.综上所述可知此三维成像方法可用于微震数据的震源定位，计算稳定且适用性较强.

致谢 感谢审稿专家提出的宝贵修改意见和编辑部的大力支持！