基于微分格式的微地震走时反演方法研究

胡 飞， 苏 奥

(1.中国科学技术大学 地球与空间科学学院，合肥 230026；2.东方地球物理公司 研究院，涿州 072750)



基于微分格式的微地震走时反演方法研究

胡 飞1，2， 苏 奥2

(1.中国科学技术大学 地球与空间科学学院，合肥 230026；2.东方地球物理公司 研究院，涿州 072750)

在页岩气开发涉及的微地震监测中，只利用声波测井资料建立的速度模型并不足够准确，通常的做法是利用测井资料与射孔数据、井下爆炸索数据、井中下落球震数据结合的方法构建模型，但射孔、爆炸索等数据的震源初始时间很难准确测量。初始时间的准确性将会影响速度模型的准确性，进一步影响定位结果的准确性。这里提出一种不需要震源初始时间的反演方法，与传统的速度反演方法相比，该方法基于不同震相的到时差信息。对同一震源，提取P-P，S-S和P-S震相到时差，这三部分基于微分格式的到时差信息，都可完全消除初始时间的影响。对合成数据测试、分析和讨论的结果，证实了该方法实现速度模型校正的可行性，并给出了该方法在实际数据应用中的示例。

微地震； 微分格式； 走时反演； 速度模型校正

0 引言

微地震监测技术可以追溯到上世纪70年代，在过去几十年里已经被证明是一种了解地下结构等信息的不可或缺的技术。在微地震成像中，速度模型的获取最为重要，因为更加准确的目标区域速度结构，可以提高微地震事件定位的精确度。另外在微地震定位结果的不确定性分析中，速度模型误差带来的定位错误要比走时拾取误差带来的定位错误大得多[1]。考虑到页岩气所在区域的沉积地质特性，工业中通用的做法是利用声波测井资料建立层状的速度模型。但是这也带来一个重要问题——该速度模型并不是足够准确的，这将会影响到后续处理，比如地震定位。速度模型不够准确主要有四方面原因，①频率不同，声波测井的频率可以高达几千赫兹，而微地震的频率一般只有几十赫兹，更高的频率带来更大的衰减系数，导致测量不准；②利用测井资料获取的速度信息只是纵向的速度，实际需要的是穿过目标区域的横向速度信息，目标区域的横向速度值相对于测井资料结果，有10%～20%的速度误差[2]；③射孔数据，井下爆炸索数据的初始时间很难准确测量[2-3]，Warpinski等[4]研究发现是由于安全射孔工艺技术导致的，该技术中的电容器放电和火花间隙会导致导火索点火被延迟；④水力压裂过程会改变目标区域的速度结构，比如该区域已经进行过一段时间的生产，会导致压力下降，或者地下流体改变，这些都会改变该区域的速度，测井得到的资料相对于生产监测时，速度信息已经不再正确[5]。因此在微地震监测工作开始前，必须要进行速度模型的校正。

一些评估水力压裂速度模型的方法已经被提出，Block等[6]提出一种参数分离方法来联合反演破裂区速度信息、震源位置、发震时刻；Warpinsk等[4]基于射孔数据和初始时间的评估，提出一种最小二乘的线性反演方法。这些方法都提高了微地震成像的效果，但是他们的方法依然存在初始时间无法精确测量的问题。Pei等[5]利用快速模拟退火法和线性奥卡姆反演方法去校正速度模型。他们把P、S震相利用联合目标函数结合起来，不过初始时间无法准确测量的问题依旧存在。

这里提出一种基于微分格式的非线性最小二乘反演方法，通过拟合微分格式的旅行时而不是直接拟合旅行时，来校正测井资料构建的速度模型。利用微分格式的走时信息，其优势在于能很好解决震源初始时刻难以准确测量的问题。将该方法应用到合成数据和实际数据的实例中，可以实现：只反演P波速度、只反演S波速度、同时反演P波和S波速度。

1 反演原理和方法

在页岩气的工业生产中，由于研究尺度较小，地层多为沉积特性，故通常建立层状结构模型。该模型多为结合测井资料与射孔，爆炸索等数据获得，但是存在震源起爆初始时间很难精确测量的问题。基于微分格式的微地震走时反演方法，试图解决这一难题。该方法通过建立最小二乘目标函数，既实现对微分格式走时信息的拟合，也实现对模型平滑的约束。以测井资料提供的信息为初始速度结构，迭代更新速度模型，至目标函数前、后两次残差达到最小且保持不变时停止迭代，从而达到校正测井速度模型的目的。

1.1 目标函数

通过拟合检波器记录的地震波走时和速度模型正演计算的走时，我们可以获取地下结构的层速度，作者提出的反演方法是通过拟合微分格式的走时而不是直接拟合走时来校正速度模型。具体来讲，对同一震源，考虑不同检波器P-P震相到时差，S-S震相到时差，P-S震相到时差，这三部分基于微分格式的到时差信息，都可完全消除初始时间的影响，所以目标函数中要加入这三个震相时差的约束项。考虑到测井资料模型虽然不够精确，但已是很好的初始速度模型，也是较为准确的真实速度结构，所以目标函数要加入初始速度模型(测井资料)的约束项。同时考虑到问题的稳定性和模型的平滑性，我们应用吉洪诺夫正则化[7]到目标函数中，并选择L2范数作为最小二乘评估算子。综上所述，可采取如下具有物理意义的目标函数Φ(m)：

Φ(m)=a0‖(Δdp-ΔG(mp))‖2+a1‖(Δds-ΔG(ms))‖2+a2‖[(ds-dp)-(G(ms)-G(mp))]‖2+τp‖L(mp-mp0)‖2+τs‖L(ms-ms0)‖2

(1)

式中：Δdp和ΔG(mp)分别是观测和计算的P-P震相走时微分；Δds和ΔG(ms)分别是观测和计算的S-S震相走时微分；mp和ms分别是P波和S波慢度；dp和ds分别是观测的P波到时、S波到时，都包括了震源初始时间；G(mp)和G(ms)分别是利用P波慢度和S波慢度计算的合成到时；a0,、a1和a2分别是P-P、S-S、P-S震相权重因子；τp和τs是模型平滑控制参数；mp0和ms0分别是从声波测井数据获取的初始P波慢度和初始S波慢度；L是正则化算子。

目标函数Φ中前三项分别给出了P-P、S-S和P-S震相的走时残差信息(观测到的地震波震相走时差与计算得到的震相走时差的残差)，后两项分别给出了对P、S波的新速度模型相对测井提供的初始模型变化量的约束。算子L采用二阶Laplace算子，考虑到模型空间或者实际地球介质的速度结构存在非线性变化，所以采用具有非线性插值约束的二阶梯度Laplace算子是合理的。通常来讲，平滑控制因子τp和τs的选取具有一定的主观性，初始值可以给1.0。权重因子a0、a1和a2可以根据地震资料各震相的信噪比选择，或者根据用户拟合某些震相的要求选取权重大小。在此作者根据P-P、S-S和P-S三震相的走时残差均方根不同，提供一种求解方案，具体公式如下：

(2)

即a0=1，a1取P-P震相与S-S震相旅行时差的均方根之比，a2取P-P震相与P-S震相旅行时差的均方根之比，其中RMS代表求解均方根值。这样选择权重的出发点是三组微分格式走时中，P-S震相到时差的量级要比P-P震相到时差，S-S震相到时差大很多。

1.2 反演方法

求解使得目标函数取值最小的模型值即新的速度结构，我们使用最小二乘的分析方法来求解目标函数的极值问题，这需要进行两部分的工作：①正演计算炮点到检波点的走时；②求解目标函数的极值。

在第一部分正演问题中需要计算炮点到检波器的走时，传统上有两种方案计算旅行时:射线追踪法和波前追踪法。射线追踪包括试射法和弯曲法，波前追踪包括有限差分法和图法。试射法和弯曲法是追踪射线从源点到检波器点路径的方法，有限差分是用程函方程去求解走时的方法，图法是基于惠更斯原理的方法。在微地震合成数据的测试中，模拟微地震生产中产生射孔数据的观测系统，首先建立层状模型，在井下进行震源激发，在井中布置检波器接收，并利用图法中的最短射线路径法(SPR)[8]计算P波、S波旅行时。考虑到实际资料处理是从检波器记录的波形上拾取到时信息作为观测走时，所以为了结论的严谨性和正确性，在合成数据测试中，采用广义反透射系数方法(GRTM)[9]先计算波形，然后手动拾取初至走时作为观测数据，在迭代过程中使用SPR方法得到计算走时数据。

在第二部分求解目标函数极值的过程中，为避免计算繁杂的海赛矩阵，通过Gauss-Newton (GN) 法对目标函数进行线性化处理，以期求得迭代形式的更新公式。首先对目标函数(公式)进行形式变换，利用L2范数的定义，目标函数可改写为：

(3)

其中：M是检波器记录的总道数；N是模型层数。公式中的各矩阵定义如下：

公式(3)相对公式(1)形式上得到大幅简化，很好地将观测数据d与计算数据G进行了结合，这将使得反演过程得到简化。根据GN方法，当目标函数Φ(m)取得极值时，有递推公式如下：

(4)

对公式(3)求导，可求得如下结果：

(5)

将公式(5)代入公式(4)可得到递推公式(6)。

(6)

在公式(6)中，各矩阵定义如下：

在矩阵Ak中,第i行第j列元素代表了第i条射线在模型第j层中的射线长度。矩阵Wk是对不同数量级的P-P、S-S和P-S震相数据进行权重控制的因子。考虑到问题的非线性情况和求解的稳定性，我们在GN方法得到的公式(6)左边加入可调节的阻尼项εkI。对于阻尼最小二乘法，最大的困难是如何选取阻尼εk值。麻省理工学院的Morgan教授[10]提出了一种方法:εk=α×σrms， 其中α是一个经验性的参数(大约0.1)，σrms是公式(6)右边部分的均方根误差。这个方法的优势在于随着每次迭代的进行，阻尼因子εk会一直变化并且被自动赋值。具体来讲，如果目标函数的值不是足够小或者拟合效果较差时，σrms将会变大，这进一步导致了一个足够大的εk被自动应用到阻尼系数中，最终会导致相对小的模型改变。另一方面，随着反演过程进行，结果越来越收敛，σrms就会减小，即一个足够小的εk将会被自动应用到阻尼系数中。该方法在一些已经发表的文章中已经被成功应用[11-13]。综上所述我们得到最终的迭代公式：

(7)

我们可以应用共轭梯度法等方法通过多次迭代来更新速度模型，迭代过程中求解观测数据和计算数据残差的均方根误差，当该均方根误差降低到稳定不变时迭代结束，迭代所得结果即新速度模型。

基于微分格式的走时反演可以分为三种情况，①只反演P波速度；②只反演S波速度；③同时反演P波和S波速度。在上面的公式推导中我们以同时反演P波和S波为例给出了公式(1)，对于另外两种情况，我们只需将对应的系数和权重因子设为“0”。

2 合成数据测试

通过上面的理论分析，我们将用合成数据进一步阐释该方法的可行性。图1是在文献[14]中使用的德克萨斯州微地震工区模型基础上，建立的一个微地震层状模型。观测系统由2个射孔震源和24个检波器组成(图1(a))，其中检波器序列位于垂直井中，间隔15 m，深度从2 070 m至2 415 m。图1(b)展示了震源到检波器的射线路径,图中震源位置为(400 m,2 600 m)和(800 m,2 600 m)。

表1给出了该速度模型的具体层数，厚度，P波和S波真实速度, 其中3 000 m以下为无限空间。从表1中可以看到，该模型有高速层，有低速层，在目标区域(可能存在页岩气的区域)附近的层位中(深度2 171 m～2 457 m)层厚最小34 m，最大119 m，目标区域厚度大约300 m，符合常见页岩气地层的主要速度特征。

测井数据提供的模型很接近真实结构，但通常与目标区域横向速度存在10%～20%[2]的误差。因此在合成数据测试中，将真实模型每层的速度进行±[10%, 20%]随机扰动后，可当作测井资料提供的速度信息，同时也可作为反演的初始模型。然后利用微分格式的走时反演方法，进行只反演P，只反演S，同时反演P、S速度。图2显示利用GRTM方法计算的两个震源A和B的波形，红色竖线和蓝色竖线分别是手动拾取的P波、S波到时。图2中横轴代表走时，单位s,纵轴是检波器序号，每个检波器记录中的三种颜色，代表了三分量的波形记录。

图1 速度速度信息及射线路径

图3分别给出了只反演P、只反演S和同时反演P、S速度的结果。图3中横轴是速度，纵轴是深度，虚线是初始速度模型，红线是真实速度模型，蓝线是反演结果。从图3中可以看到，三种反演方案在目标区域层都得到了很好的速度恢复，而顶、底两层的恢复效果较差，从图1(b)看到射线路径未经过顶、底两层，故无法进行这两层的速度校正，这与图3中结果一致。而图3显示顶、底两层的速度仍发生了变化，这是因为目标函数中二阶Laplace算子对模型慢度参数的运算所致。在实际处理中，考虑到该运算对慢度改变量的未知性，建议这类层位在校正后仍沿用测井资料提供的速度信息。

表1 模型参数

图2 利用GRTM方法计算得到两个震源的波形

图3 不同波形反演结果

图4 震源初始时间T0的影响

3 分析与讨论

基于微分格式的微地震走时反演方法主要特点，解决震源起爆初始时间难以准确测量的困难，以测井资料提供的初始速度结构为约束，来进行地震波层速度的校正。考虑到初始时间和初始模型在反演过程中的重要性，将对初始时间T0和初始速度模型对反演结果的影响分别讨论。

3.1 震源初始时间T0对反演结果的影响

微地震研究区域通常较小，同时在井下激发和接收，使得微地震观测走时量级上偏小，较小的初始时间误差也会对速度模型的正确性产生较大影响。图4给出了震源初始时间T0对反演结果的影响。测试中，对正确T0加上[-2,2] ms随机扰动作为错误T0。图4的结果显示，较小的震源初始时间偏差，就会导致利用绝对走时反演得到的速度模型出现错误，这进一步说明了初始时间的重要性，同时可以看到作者提出的方法在震源初始时间存在错误时，仍然可以得到准确的速度结构，表明了该方法校正速度模型的可行性。

3.2 初始模型对反演结果的影响

在合成数据的测试中，用来模拟测井资料的初始速度结构是通过真实模型随机扰动得到的。为了进一步说明反演方法的可行性，进行了不同常速度初始模型的讨论。图5给出了不同常速度初始模型所对应的反演结果。从图5的结果可以看到，基于微分格式的反演方法对初始模型并没有太强的依赖，对初始模型有一定的抵抗性。图5中初始速度模型分别为500 m/s、 1 000 m/s、4 000 m/s、8 000 m/s。

4 应用实例

将微分格式走时反演方法应用到一个微地震的下落球震(dropball data)实际资料中，并给出相应的反演结果。图6是从声波测井资料得到的初始速度模型(速度-深度图)，图6中左边为S波速度曲线，右边为P波速度曲线。图7(a)给出了此次工程中四个球震震源和12个检波器组成的观测系统图，

图5 初始速度模型对反演结果影响

图7中星号(*)为四个下落小球震源，倒三角(▽)是12个检波器。由于是层状模型，将3维震源坐标 旋转到x-z平面内后，四个震源坐标从左至右依次为(118.56 m，1 423.44 m)，(230.29 m,1 419.96m)，(534.41 m,1 422.93 m)，(863.86 m,1 425.73 m)。12个检波器坐标从1 120 m至1 285 m，间隔为15 m。

图6 由声波测井得到的P，S波初始速度模型

此次研究的目标区域位于深度1 100 m至1 450 m，将图7(a)中目标区域(红色方框所示)放大得到图7(b)，利用SPR计算的射线路径也展示在图7(b)中。目标区域以外并无射线经过，所以建议这些无射线经过的层位最终速度仍沿用测井资料提供的速度信息。

图7 实际资料模型和射线路径

图8显示了12个三分量检波器记录到的波形信息，这里只展示检波器记录到的最左边震源的波形。图8中横轴代表走时，单位s，纵轴是检波器序号，每个检波器记录中的三种颜色，代表了三分量的波形记录，红竖线和蓝竖线分别为手动拾取的P、S波初至时间。

图8 球震数据实际波形和手动拾取的P、S到时

图9给出了实际数据只反演P、只反演S、同时反演P、S的结果。由于不知道真实速度模型，故图9中只有初始模型和反演结果。

对于实际数据，由于不知道真实的地下结构，虽然不能通过反演后走时曲线的拟合图来判断速度结构正确性，但是走时拟合图可以作为检验结果正确性的一种辅助方式。图10 (分别给出了只反演P、只反演S、同时反演P和S时走时曲线拟合图，在图中上下左右四块分别对应了沿x方向的四个震源。考虑到实际数据处理中存在噪音，拾取误差等，可认为走时曲线拟合效果较好。

5 结论

这里介绍了页岩气开采中，以测井资料为基础，校正微地震速度模型的一种可行方案。利用微分格式的走时信息，很好地解决了射孔数据等震源初始时间不能测量或者很难准确测量的困难。作者尽可能组合更多走时信息：P-P震相到时差、S-S震相到时差、P-S震相到时差，这三部分信息都可以排除初始时间的影响。合成数据的测试结果，以及初始时间、初始模型的讨论结果，都证明了该方法的可行性。还通过一个球震实际资料给出了实际应用的示范。该方法主要有两方面的贡献：①解决了震源初始发震时间很难准确测量的问题，②将旅行时以微分格式应用到速度反演中。

图9 实际球震数据反演结果图

图10 四个震源的走时曲线拟合图

[1] JONES G A, KENDALL J M, BASTOW I D, et al. Locating microseismic events using borehole data [J]. Geophysical Prospecting, 2014, 62(1): 34-49.

[2] WARPINSKI N. Microseismic monitoring: Inside and out [J]. Journal of Petroleum Technology, 2009, 61(11): 80-85.

[3] PEI D, QUIREIN J A, COMISH B E, et al. Velocity calibration for microseismic monitoring: Applying smooth layered models with and without perforation timing measurements [C]. SPE Annual Technical Conference and Exhibition. Society of Petroleum Engineers, 2008.

[4] WARPINSKI N R, SULLIVAN R B, UHL J, et al. Improved microseismic fracture mapping using perforation timing measurements for velocity calibration [J]. SPE Journal, 2005, 10(01): 14-23.

[5] PEI D, QUIREIN J A, COMISH B E, et al. Velocity calibration for microseismic monitoring: A very fast simulated annealing (VFSA) approach for joint-objective optimization [J]. Geophysics, 2009, 74(6): WCB47-WCB55.

[6] BLOCK L V, CHENG C H, FEHLER M C, et al. Seismic imaging using microearthquakes induced by hydraulic fracturing [J]. Geophysics, 1994, 59(1): 102-112.

[7] TIKHONOV A N, ARSENIN V I A, JOHN F. Solutions of ill-posed problems [M]. Washington, DC: Winston, 1977.

[8] MOSER T J. Shortest path calculation of seismic rays [J]. Geophysics, 1991, 56(1): 59-67.

[9] CHEN X. Seismogram synthesis for multi- layered media with irregular interfaces by global generalized reflection/transmission matrices method. I. Theory of two-dimensional SH case [J]. Bulletin of the Seismological Society of America, 1990, 80(6A): 1696-1724.

[10] MORGAN F D O. Electronics of sulfide minerals: implications for induced polarization [D].Massachusetts Institute of Technology, 1981.

[11] MACKIE R L, MADDEN T R. Three- dimensional magnetotelluric inversion using conjugate gradients [J]. Geophysical Journal International, 1993, 115(1): 215-229.

[12] ZHANG J, MACKIE R L, MADDEN T R. 3-D resistivity forward modeling and inversion using conjugate gradients [J]. Geophysics, 1995, 60(5): 1313-1325.

[13] ZHANG J, TOKSOZ M N. Nonlinear refraction traveltime tomography [J]. Geophysics, 1998, 63(5): 1726-1737.

[14] WONG J, MANNING P M, HAN L, et al. Synthetic microseismic datasets . Focus, 2011.

Microseismic travel-time inversion based on time difference

HU Fei1,2, SU Ao2

(1.School of Earth and Space Science, University of Science and Technology of China, Hefei 230026, China;Geophysical Research Institute, BGP, CNPC, Zhuozhou 072750, China)

During microseismic monitoring for unconventional shale gas exploitation, initial velocity model is constantly built from sonic data or a combination of sonic data and perforation shot (or string shot, drop ball). Unfortunately, it's quite difficult to obtain the origin time. The origin time of perforation shot affects the accuracy of velocity structure, consequently the accuracy of subsequent location results. Therefore, we demonstrate an inversion procedure to calibrate the velocity model without origin time but using differential arrival times. Comparing with traditional methods, we combine travel-time data, such as differential P-P arrival times, differential S-S arrival times, and differential S-P arrival times. Both synthetic examples and real data results show the improvement and effectiveness of the inversion method in velocity model calibration.

microseimic; differential; travel-time inversion; velocity calibration

2014-07-25 改回日期：2014-10-11

胡飞(1988-)，男，硕士，从事微地震地球物理研究，E-mail:hufei@mail.ustc.edu.cn 。

1001-1749(2015)04-0478-10

P 631.4

A

10.3969/j.issn.1001-1749.2015.04.11