2017, Vol. 60
2. 北京大学地球与空间科学学院石油与天然气研究中心, 北京 100871;
3. 中国石化石油勘探开发研究院, 北京 100083
2. Institute of Oil & Gas, School of Earth and Space Sciences, Peking University, Beijing 100871, China;
3. Petroleum Exploration & Production Research Institute, SINOPEC, Beijing 100083, China
水力压裂微地震监测技术是通过观测注水压裂过程所产生的微弱地震事件来描述地下裂缝的形态特征及分布规律的一项技术.由于该技术直接利用微地震的震源位置描述裂缝，因此，震源定位是该技术要解决的最主要问题.目前，在微地震震源定位方面常用的方法包括两大类，即走时拟合法和扫描叠加法(Bardainne et al.，2009;Pesicek et al.，2014).这两类方法分别利用微地震事件的走时和波形信息来反演震源位置，并适用于不同的观测条件.对于井下监测来说，其优点是观测系统更加靠近震源，接受到的微地震信号较强且环境噪声水平较低，微地震事件P波和S波初至能够较准确地识别出来，故一般采用走时拟合法进行震源定位(Maxwell et al.，2010;Jones et al.，2014).而对于地面观测系统，由于距离震源区域较远，微地震信号能量在传播过程中衰减严重，加之环境噪声的影响，造成地面记录中难从识别出有效的微地震事件.由于扫描叠加法无需提取有效事件，因此，常用于地面微地震资料的处理(Chambers et al.，2010;Duncan and Eisner，2010).
本文研究对象为走时拟合法.该方法包括两个主要步骤：一是根据P波偏振特征确定震源方位角，二是利用P波和S波到时计算震源的径向距离和纵向深度；之后，通过坐标转换即可得到震源的空间位置(Jones et al.，2010;Jones et al.，2014).震源方位角可利用矢端图法求得，而震源的距离和深度则通常利用网格搜索算法得到.此外，一些学者还将各类非线性反演方法用于震源定位.例如宋维琪和杨晓东(2011)结合网格搜索算法与遗传算法的特点，提出了一种解域约束下微地震事件联合反演方法；宋维琪等(2013)研究了微地震贝叶斯差分进化反演方法；盛冠群等(2014)提出了基于粒子群差分进化算法的定位方法.实际应用表明，上述方法的定位精度优于传统的网格搜索算法，然而它们也存在某些不足之处，例如反演结果容易陷入局部极小值，每定位一个事件需要重新计算理论到时等.网格搜索算法则不存在上述问题.理论上讲，只要网格划分的足够小，可得到任意精度的定位结果，但这样做会极大增加处理时间.总的来说，网格搜索算法由于其原理简单、容易实现等优点，目前仍然是微地震资料处理中最常用的定位方法.
本文针对基于网格搜索的微地震震源定位方法做出了几点改进，并采用模型数据和实际资料对改进算法的应用效果进行了验证，同时讨论了影响震源定位结果精度的主要因素.
2 方法原理 2.1 震源方位角的计算矢端图法采用三分量检波器两个水平分量的记录建立P波振动轨迹在水平面内的投影，对该投影轨迹进行线性拟合得到P波的偏振方向.在实际应用中，不同位置的检波器接收到同一微地震信号的信噪比不同，造成各道记录的P波偏振特征存在一定差别.一般来说，对于信噪比较高的道，其P波振动轨迹接近线性偏振，采用线性拟合求得的方位角较准确；反之，对于信噪比较低的道，其P波振动轨迹趋于椭圆偏振，此时通过线性拟合求得的方位角则可能存在较大误差.为了解决该问题，本文在采用偏振分析方法的基础上，通过构建一个概率密度函数来求取震源的方位角.
对于微地震事件中任意一道记录，本文方法首先选用一个短时窗截取包含其P波初至的一小段记录.假设xi和yi(i=1，2，…，N)分别表示该段记录的两个水平分量，可构造如下协方差矩阵：
|
(1) |
利用该协方差矩阵的特征值和特征向量可计算该段记录的线性度η和偏振角α(Flinn，1965; Montalbetti and Kanasewich，1970)，即
|
(2) |
|
(3) |
其中，λ1，λ2(λ1>λ2)为协方差矩阵M的特征值，u1=[ux，uy]为λ1所对应的特征向量.η的范围为0~1，其中η=1表示线性偏振，η=0表示圆偏振.
利用线性度η和偏振角α可定义如下概率密度函数：
|
(4) |
|
(5) |
其中，γ为转换系数，其值约为0.025；θ为方位角，其取值范围为0~180°.由(4)式可知，当η的值接近1时，P(θ)近似为脉冲信号，且仅在θ=α时存在非0值，表明方位角即为α；而当η的值接近0时，P(θ)近似为白谱，此时方位角可取任意值.在实际应用中，由于环境噪声的存在，微地震事件P波振动轨迹通常为椭圆，则η的值为一个介于0~1之间的实数.此时P(θ)满足高斯分布，其最大值位于α，表明α为实际震源方位角的概率最大.图 1显示了当η取不同值时的概率密度函数P(θ).
|
图 1 η取不同值时的概率密度函数P(θ) Fig. 1 Probability density function P(θ) for different η |
采用上述步骤对微地震事件中各道记录进行处理后，可得到其对应的概率密度函数P(θ)j(j=1，2，…，M).定义总的概率密度函数为
|
(6) |
Γ(θ)最大值所对应的角度即为求得的震源方位角.
2.2 震源距离及深度的计算基于网格搜索的微地震震源定位方法通过计算走时残差作为目标函数来确定震源位置.走时残差即为观测走时与理论走时之差，其中观测走时可由检波器记录到的微地震事件P波和S波的初至到时减去发震时刻得到.由于微地震事件的发震时刻和震源位置一样，都是未知(待求)参数，某些定位方法选择将最早记录到的P波初至到时作为发震时刻，但该方法过于依赖所选取的参考发震时刻的精度(Prugger and Gendzwill，1988).本文采用观测到时与理论走时之差的平均值作为发震时刻的估值，得到单一震相的走时残差计算公式为(Nelson and Vidale，1990)
|
(7) |
|
(8) |
其中，Tcj为实际观测到时(c表示P或S)，tcj为理论计算走时.
一般来说，微地震事件既包含P波，又包含S波，因此，可利用二者的到时差来计算走时残差，即
|
(9) |
本文将F1与F2结合，提出一种新的走时残差计算方法，即
|
(10) |
其中ρ为均衡系数，其取值范围为0~1，在实际应用时需要通过模型试验给出ρ的估计值.首先采用模型算例比较上述三类走时残差计算方法.该算例所采用的地层模型为一个均匀介质模型，其P波和S波速度分别为2000 m·s-1和1200 m·s-1.震源位置为(600，1500)m，观测系统为一个20级的井下检波器串，其级间距为10 m，首级位置为(100，1300)m.假设震源的发震时刻为100 s，利用射线追踪算法正演计算地震波的理论P波和S波到时，之后，根据(7)式、(9)式和(10)式计算得到震源附近区域内的目标函数等值线图如图 2所示.由该图可知，F1在水平方向上的收敛性较差，而F2在垂向上的收敛性较差，由此可推断当观测到时存在误差时，二者的定位结果在距离和深度方面可能存在较大误差.相比较而言，F3的收敛性最好，采用该式作为目标函数可以得到更加准确的定位结果.
|
图 2 三类目标函数对比 (a)F1;(b)F2;(c)F3. Fig. 2 Comparison of the three objective functions |
传统的网格搜索算法需要对搜索目标区域进行网格剖分，并计算各个网格点的目标函数值.本文在借鉴Font等(2004)提出的最大交汇算法的基础上，提出了一种改进的搜索算法.假设TPj和TSj分别表示实际观测的P波和S波初至到时，tPj和tSj为理论计算走时.对于空间中的任意一点，可定义
|
(11) |
|
(12) |
|
(13) |
在观测到时不含任何误差的情况下，空间中所有满足S=0的点可构成一系列曲面，称为EDT面(Zhou，1994)，震源位置即位于它们共同的交点处(如图 3a所示).
|
图 3 EDT面示意图 (a)EDT面厚度为0时的情况;(b)EDT面厚度不为0时的情况. Fig. 3 Illustration of the EDT surface (a)The case that the thickness of EDT surface is zero;(b)The case that the thickness of EDT surface is not zero. |
由于实际观测到时不可避免地会存在拾取误差，因此，所有EDT面通常并不交于一点.为了解决该问题，一种方法是通过对这些曲面进行‘增厚’使其产生共同的交点，即满足
|
(14) |
其中δ为误差阈值，其大小决定了EDT面的厚度.一般来说，(14)式定义的EDT面相交于一个较小的区域(如图 3b所示)，震源位置即位于该区域内.在实际应用中，若某个观测点的初至到时存在较大误差(即存在异常值)，则可能造成与该观测点相关的EDT面不会与其他EDT面存在共同的交汇区域.根据这一点可以判断实际观测到时中是否存在异常值.
本文算法的具体步骤为(图 4)：
|
图 4 本文方法流程图 Fig. 4 Flowchart of the proposed method |
(1) 根据(6)式计算震源方位角，确定震源所在垂直平面.
(2) 设定搜索目标区域，对该区域采用均匀网格进行剖分；采用射线追踪算法正演计算每个网格点到达各个检波器的理论P波、S波走时.
(3) 对每个网格点，根据(14)式判断各个EDT面是否经过该点，并统计经过该点的EDT面的数目及标号.
(4) 筛选出经过最多EDT面的所有网格点；若仅存在一个网格点，则该点即为震源位置；若存在多个网格点，则进行下一步.
(5) 对步骤(4)筛选出的各个网格点，判断其是否与所有EDT面相交；若是，进行下一步；否则，根据经过该点的EDT面的标号判断观测到时中异常值的位置，并将其剔除.
(6) 采用(10)式及剔除异常值后的观测到时计算步骤(4)筛选出的各个网格点的目标函数值，选取目标函数值最小点为震源位置.
3 模型数据处理结果首先采用模型数据检验本文方法的应用效果.图 5为模拟的微地震事件P波初至波形.该波形记录由一个线性偏振的地震子波与随机噪声叠加得到，子波的偏振角为45°，噪声随道号逐渐增大(其均值为0，方差分别为0.01、0.015、0.02、0.025、0.03、0.035、0.04、0.045).图 6显示了各道记录的振动轨迹，采用矢端图法求得的各道记录的方位角见表 1.通过对比分析图 6及表 1中的数据可知，随着噪声水平的增大，信号的振动轨迹逐渐接近椭圆偏振，根据矢端图法求得的方位角的误差也逐渐增大，表明该方法受噪声影响较大，在低信噪比条件下难以取得较理想的结果.
|
图 5 模拟微地震事件P波初至波形 (a)X分量;(b)Y分量. Fig. 5 Simulated P-wave waveforms of microseismic event (a)X-component;(b)Y-component. |
|
图 6 图 5中各道记录的振动轨迹及概率密度函数 Fig. 6 Trajectories and probability density functions of the waveform data shown in Fig. 5 |
|
表 1 利用矢端图法求得的图 5中各道记录的震源方位角 Table 1 Source azimuths of the waveform data shown in Fig. 5 estimated by the hodogram method |
利用(2)式、(3)式计算图 5中各道波形记录的偏振特征参数，结果见表 2.根据各道偏振特征参数计算的概率密度函数如图 6所示.图 7为根据(6)式计算得到的总的概率密度函数，其最大值(用“*”表示)所对应的方位角为44.8°，与真实值仅相差0.2°.为了得到更一般的结论，将上述试验重复进行8次，每次试验采用100组具有相同噪声水平的合成记录，得到的结果如图 8所示.通过该图可以看出，随着噪声水平的增大，采用矢端图法求得的方位角逐渐偏离真实值(用虚线表示)，而本文方法的结果则在真实值附件摆动，表明利用本文方法能够更加准确地求出震源方位角.
|
图 7 总的概率密度函数 Fig. 7 Overall probability density function |
|
图 8 不同噪声水平时矢端图法与本文方法计算得到的震源方位角对比 Fig. 8 Comparison of the source azimuths estimated by the hodogram method and the proposed method for different noise levels |
研究表明，影响微地震震源定位精度的主要因素包括观测系统布设位置、初至拾取误差和速度模型误差等(Eisner，2009;Maxwell，2009).具有较大覆盖范围和良好观测角度的观测系统对于提高震源定位精度至关重要，这是由于不同方位的射线交叉接收对于解是“天然”的约束.然而在实际应用时，尤其针对水力压裂微地震监测，受到施工条件、花费成本以及处理方法等因素制约，往往很难采用最优的观测系统.单井监测是目前水力压裂微地震监测普遍采用的观测方式，因此，本文将主要采用该观测系统进行分析讨论.本文在对比分析不同目标函数定位结果的同时，讨论初至拾取误差和速度模型误差对定位结果的影响.图 9为一个一维水平层状速度模型，其模型参数见表 3.本次试验采用的观测系统为一个25级的井下检波器串.假设震源位置为(50，2250)m，发震时刻为100 s，利用射线追踪算法正演得到由该震源位置传播到各个检波器的理论P波和S波到时如图 10所示.
|
图 9 地层模型及震源、观测系统的布设位置 Fig. 9 Layer model and the locations of the source and receivers |
|
|
表 3 地层模型参数 Table 3 Parameters of the layer model |
|
图 10 理论观测到时 (a)P波到时;(b)S波到时. Fig. 10 Theoretical arrival times (a)P-wave arrival time;(b)S-wave arrival time. |
本次试验采用一组随机数(其平均值为0，方差为0.002)表示初至拾取误差，并将该组随机数加到正演到时中模拟实际观测数据.为了避免随机性，重复上述步骤得到50组观测到时数据.对这50组数据分别采用(7)—(10)式作为目标函数进行震源定位的结果如图 11所示.通过对比图 2与图 11可知，反演震源位置与目标函数等值线的分布规律基本一致：(7)式在水平方向上的收敛性较差，其定位结果在距离方面存在较大误差；(9)式在垂向上的收敛性较差，其定位结果在深度方面存在较大误差；而(10)式的收敛性最好，其定位结果的误差最小.
|
图 11 三类目标函数定位结果对比(含初至拾取误差) (a)F1;(b)F2;(c)F3. Fig. 11 Comparison of the location results of the three objective functions(with arrival picking errors) |
本次试验通过对真实地层速度模型进行随机扰动模拟速度模型误差.同样地，为了避免结果的随机性，生成了50组模型数据.根据这50组模型数据及理论观测到时得到的定位结果如图 12所示.该结果与图 11所示结果具有相似的规律，进一步证明了本文提出的目标函数具有较好的收敛性.
|
图 12 三类目标函数定位结果对比(含速度模型误差) (a)F1;(b)F2;(c)F3. Fig. 12 Comparison of the location results of the three objective functions(with velocity model errors) |
为了检验本文提出的改进搜索算法的有效性，本次试验在搜索目标区域内随机选取了50个点作为震源位置，采用上述方法得到50组观测到时数据.此外，为了模拟由于某几道数据信噪比过低而出现误拾的情况，对每组观测到时添加了若干个异常值.对这50组数据分别采用传统的网格搜索算法及本文提出的改进算法进行定位后得到的结果如图 13所示.图 14为两种方法定位误差的比较.通过分析图 13与14可知，本文算法的定位结果更加接近真实震源位置，其平均定位距离误差仅为2.95 m；而传统网格搜索算法的定位结果则存在较大误差，其平均定位误差高达18.5 m.
|
图 13 传统网格搜索算法与本文提出改进算法定位结果对比 Fig. 13 Comparison of the location results of the conventional grid search algorithm and the improved algorithm |
|
图 14 传统网格搜索算法与本文提出改进算法定位误差对比 Fig. 14 Comparison of the location errors of the conventional grid search algorithm and the improved algorithm |
本文处理的实际资料为某油田对一口水平井进行的11段水力压裂施工的微地震监测数据.在采集该资料时所使用的观测系统为一个15级的井下检波器串，其级间距为10 m，时间采样间隔为0.0005 s.实际监测时，该观测系统被布设在压裂井附近的一口直井中进行数据采集.在进行震源定位前，首先采用工区内一口井的声波测井曲线建立初始速度模型(如图 15所示)，其模型参数见表 3；之后，利用各个压裂段的射孔记录对该初始速度模型进行校正；最后，采用校正后的速度模型对微地震震源进行定位.
|
图 15 测井曲线及地层模型 Fig. 15 Sonic logging data and layer model |
|
|
表 4 地层模型参数 Table 4 Parameters of the layer model |
图 16为采用本文方法得到的微地震震源分布图.通过该图可以看出裂缝的延伸方向为西北-东南方向，且基本穿越压裂目的层段.
|
图 16 本文方法震源定位结果 (a)俯视图;(b)侧视图. Fig. 16 Source location results using the proposed method (a)Plan view;(b)Side view. |
为了验证本文方法的应用效果，我们对实际资料进行了多次定位处理.在每次定位过程中，对某个步骤采用不同的方法进行处理，并与本文方法的定位结果进行对比分析.例如，图 17为采用矢端图法求取震源方位角得到的定位结果，图中震源位置基本符合裂缝的延伸方向.与图 16所示结果对比可以看出，图 16中震源位置主要集中在井筒附近，而图 17中震源则分布较均匀.该结果可能造成解释的裂缝尺寸大于实际情况.
|
图 17 矢端图法震源定位结果 (a)俯视图;(b)侧视图. Fig. 17 Source location results using the hodogram method (a)Plan view;(b)Side view. |
图 18和19分别为采用目标函数F1和F2得到的定位结果.根据图 18无法确定裂缝的延伸方向.图中大部分震源位置集中在监测井附近，其余震源则分布较分散.该结果进一步证明了目标函数F1在确定震源距离方面存在较大误差.采用目标函数F2得到的定位结果基本符合实际情况.然而，通过对比图 16b与图 19b可以看出，图 19b中部分震源位置位于压裂目的层段之外，表明其深度可能存在一定误差.
|
图 18 采用目标函数F1得到的震源定位结果 (a)俯视图;(b)侧视图. Fig. 18 Source location results using the objective function F1 (a)Plan view;(b)Side view. |
|
图 19 采用目标函数F2得到的震源定位结果 (a)俯视图;(b)侧视图. Fig. 19 Source location results using the objective function F2 (a)Plan view;(b)Side view. |
图 20为传统网格搜索算法的定位结果.图中除个别事件外，大部分震源位置与图 16所示结果相近.这是由于本文处理的实际微地震事件均为初至较清晰的事件，其初至拾取误差较小，故而造成两种方法的定位结果并无明显差异.
|
图 20 传统网格搜索算法震源定位结果 (a)俯视图;(b)侧视图. Fig. 20 Source location results using the conventional grid-search algorithm (a)Plan view;(b)Side view. |
本文提出了一种改进的基于网格搜索的微地震震源定位方法.在求取震源方位角时，本文方法考虑了环境噪声对P波振动轨迹的影响，根据P波的偏振特征参数计算概率密度函数求取震源方位角.在确定震源方位角后，本文方法采用改进的目标函数和搜索算法求取震源的径向距离和深度.模型数据和实际资料的处理结果表明，本文方法具有较强的抗噪性，计算得到的震源方位角更加接近真实值；本文提出的目标函数具有较好的收敛性，其定位结果受初至拾取误差和速度模型误差影响较小；本文提出的搜索算法能够消除由于错误拾取造成的观测到时中的异常值对定位结果的影响.
| Bardainne T, Gaucher E, Cerda F, et al. 2009. Comparison of picking-based and waveform-based location methods of microseismic events:Application to a fracturing job.//SEG Technical Program Expanded Abstract 2009. Houston, Texas, 1547-1551, doi:10.1190/1.3255144. | |
| Chambers K, Kendall J M, Brandsberg-Dahl S, et al. 2010. Testing the ability of surface arrays to monitor microseismic activity. Geophysical Prospecting, 58(5): 817-826. DOI:10.1111/j.1365-2478.2010.00893.x | |
| Duncan P M, Eisner L. 2010. Reservoir characterization using surface microseismic monitoring. Geophysics, 75(5): 75A139-75A146. DOI:10.1190/1.3467760 | |
| Eisner L, Duncan P M, Heigl W M, et al. 2009. Uncertainties in passive seismic monitoring. The Leading Edge, 28(6): 648-655. DOI:10.1190/1.3148403 | |
| Flinn E A. 1965. Signal analysis using rectilinearity and direction of particle motion. Proceedings of the IEEE, 53(12): 1874-1876. DOI:10.1109/PROC.1965.4462 | |
| Font Y, Kao H, Lallemand S, et al. 2004. Hypocentre determination offshore of eastern Taiwan using the maximum intersection method. Geophysical Journal International, 158(2): 655-675. DOI:10.1111/j.1365-246X.2004.02317.x | |
| Jones G A, Raymer D, Chambers K, et al. 2010. Improved microseismic event location by inclusion of a priori dip particle motion:a case study from Ekofisk. Geophysical Prospecting, 58(5): 727-737. DOI:10.1111/j.1365-2478.2010.00873.x | |
| Jones G A, Kendall J M, Bastow I D, et al. 2014. Locating microseismic events using borehole data. Geophysical Prospecting, 62(1): 34-49. DOI:10.1111/1365-2478.12076 | |
| Maxwell S C. 2009. Microseismic location uncertainty. CSEG Recorder, 34(4): 41-46. | |
| Maxwell S C, Rutledge J, Jones R, et al. 2010. Petroleum reservoir characterization using downhole microseismic monitoring. Geophysics, 75(5): 75A129-75A137. DOI:10.1190/1.3477966 | |
| Montalbetti J F, Kanasewich E R. 1970. Enhancement of teleseismic body phases with a polarization filter. Geophysical Journal of the Royal Astronomical Society, 21(2): 119-129. DOI:10.1111/j.1365-246X.1970.tb01771.x | |
| Nelson G D, Vidale J E. 1990. Earthquake locations by 3-D finite-difference travel times. Bulletin of the Seismological Society of America, 80(2): 395-410. | |
| Pesicek J D, Child D, Artman B, et al. 2014. Picking versus stacking in a modern microearthquake location:Comparison of results from a surface passive seismic monitoring array in Oklahoma. Geophysics, 79(6): KS61-KS68. DOI:10.1190/GEO2013-0404.1 | |
| Prugger A F, Gendzwill D J. 1988. Microearthquake location:a nonlinear approach that makes use of a simplex stepping procedure. Bulletin of the Seismological Society of America, 78(2): 799-815. | |
| Sheng G Q, Li Z C, Wang W B, et al. 2014. A source location method for microseismic monitoring based on particle swarm optimization combined with differential evolution algorithm. Acta Petrolei Sinica (in Chinese), 35(6): 1172-1181. DOI:10.7623/syxb201406015 | |
| Song W Q, Yang X D. 2011. Micro-seismic events recognition and inversion in the case of single seismic phase. Chinese J. Geophys. (in Chinese), 54(6): 1592-1599. DOI:10.3969/j.issn.0001-5733.2011.06.019 | |
| Song W Q, Gao Y K, Zhu H W. 2013. The differential evolution inversion method based on Bayesian theory for micro-seismic data. Chinese J. Geophys. (in Chinese), 56(4): 1331-1339. DOI:10.6038/cjg20130427 | |
| Zhou H W. 1994. Rapid three-dimensional hypocentral determination using a master station method. Journal of Geophysical Research, 99(B8): 15439-15455. DOI:10.1029/94JB00934 | |
| 盛冠群, 李振春, 王维波, 等. 2014. 水力压裂微地震粒子群差分进化定位算法. 石油学报, 35(6): 1172–1181. DOI:10.7623/syxb201406015 | |
| 宋维琪, 杨东晓. 2011. 单震相微地震事件识别与反演. 地球物理学报, 54(6): 1592–1599. DOI:10.3969/j.issn.0001-5733.2011.06.019 | |
| 宋维琪, 高艳珂, 朱海伟. 2013. 微地震资料贝叶斯理论差分进化反演方法. 地球物理学报, 56(4): 1331–1339. DOI:10.6038/cjg20130427 | |