地球物理学报  2011, Vol. 54 Issue (6): 1592-1599   PDF    
单震相微地震事件识别与反演
宋维琪 , 杨晓东     
中国石油大学(华东)地球资源与信息学院, 山东 东营 257061
摘要: 为了对单震相微地震事件进行识别, 同时将识别出来的微地震事件进行定位.根据单一震相任意两道到时差与微地震事件、检波器空间位置及震相速度关系的特征规律, 研究了单震相微地震事件识别方法.首先分析到时差与以上各变量的内在变化规律, 建立起到时差与各变量之间的定量计算关系, 然后就相邻道到时差和检波器排列的首尾两道到时差, 研究了具体的定量计算关系表达式.利用以上关系式能够计算出上述两种到时差的分布区间.把计算的实际微地震事件的两种到时差与设定的事件的到时差的分布区间对比, 以落入区间判断识别单一震相.综合搜索法和遗传算法的特点, 提出基于正演模型迭代的解域约束下微地震事件联合反演方法.从迭代误差解的期望概率密度分布、不同解方向的变化特征分析, 研究了解域约束的界定条件.在确定的解域范围内, 利用遗传算法进行微地震事件的反演.将研究的单震相识别方法, 通过对实际资料的应用验证分析, 极大地减小了震相的误识率.根据识别资料, 利用研究的反演方法对实际资料进行反演.经实际情况对比分析, 反演结果和实际的裂缝分布带是一致的, 并且反演结果聚敛性效果较好.
关键词: 微地震      单震相识别      震相时差      解域约束      联合反演     
Micro-seismic events recognition and inversion in the case of single seismic phase
SONG Wei-Qi, YANG Xiao-Dong     
College of Resources and Information, China University of Petroleum (East China), Dongying 257061, China
Abstract: In order to identify and locate the micro-seismic events of single seismic phase, we have developed a method on the basis of the arrival time difference between any two traces in the case of single seismic phase and the characteristic relationship of micro-seismic events, spatial location of detector and the seismic phase velocity. First, by analyzing the variation of arrival time difference with the change of above variables, we establish the quantitative calculation relationship between the arrival time difference and every variable. Then, with regard to the arrival time differences of the adjacent trace and the first and the last trace of detector array, we research the specific expression of the quantitative calculation relationship. Given the spatial location of the event, we can calculate the distribution range of the above two kinds of arrival time difference using the relation expression above. Comparing the distribution range of two arrival time differences of the actual micro-seismic events with that calculated for the given event, we can identify the single seismic phase with the range which they fall to. Integrating the advantages of the search method and the genetic algorithm, we put forward a micro-seismic joint inversion method based on the solution domain constrained iterative forward model. From the analysis of the expectation probability distribution of the iterative error solution and the variations of different direction of the solutions, we establish the definition conditions of the domain of solution constraints. In the determinate scope of the solution domain, we carry on the inversion of micro-seismic events using the genetic algorithm. Through the analysis of application of practical data, the method of identifying the micro-seismic events reduces the error rate greatly. According to the identified information, we do the inversion of actual data using the inversion method. After comparing with the actual situation, the result is consistent with the actual crack distribution, and is more convergent..
Key words: Micro-seismic      Single seismic phase recognition      Seismic phase time difference      Constrained solution domain      Integrating inversion     
1 引言

微地震监测技术,是一种通过观测、分析生产活动中所产生的微小地震事件来监测生产活动的影响、效果及地下状态,在油气田开发阶段,能够实现储层的压裂增产,注水注气后的油藏驱动监测,对油田开发稳产、高产具有重要的意义.无论P 波、S 波联合反演还是P波或S波单震相事件反演,震相识别及初至拾取都是微地震事件反演的基础环节.天然地震领域提出并讨论了震相识别问题,但是微地震事件的震相识别问题,既和天然地震有相似之处,又存在着很大的区别[1].水力压裂裂缝微地震监测过程中,裂缝开启或闭合产生的微地震事件震相类型是随机的,既可能同时产生P 波S 波,也可能只产生单一的P 波或S 波.由于地层因素、井周边环境影响及各类噪音的干扰,使实际记录到的微地震资料更为复杂,如何从复杂的微地震记录中识别分离出有效微地震事件,同时确定这些有效微地震事件的震相类型,对今后的反演工作至关重要[2-5].

高精度微地震震源反演通常是基于模型计算的走时约束反演,决定反演精度的主要因素有:拾取初至的精度、速度模型和选择的反演算法.在实际生产中,大多是纵波声波时差资料,横波声波时差资料较少,在这种情况下,实际S波微地震反演横波速度初始模型的建立,一般采用纵波资料乘以某个系数的办法,研究发现这种横波速度建模近似程度大,不能满足高精度反演定位的需要[6-10].

基于模型正演的微地震反演参数变量主要包括:径向距离、纵向深度、层位厚度及速度,属于多变量反演问题.针对这种反演问题,人们采用了各式各样的线性、非线性反演方法,基本概括为三大类:(1)迭代修正法是以最小迭代误差作为解寻优判别方法,但易陷入局部极小.(2)贝叶斯反演方法,虽然考虑了全局最优问题,对解的评价采用了迭代误差和迭代解与初始解之误差相加权的统计方法,但只是对总误差进行了评价,没有将评价后的结果再返回到迭代计算过程中,使问题的解不能完全最优.(3)搜索法是微地震反演中应用较多的方法,一般是采用粗细网格相结合的解搜索法,虽然也采用了迭代误差的概率密度的统计,找出误差分布期望较大的区域再进行细网格搜索,但是这又容易限于局部极小问题中去[11-14].

针对微地震震相识别、速度模型建立与校正、反演解的评价等方面存在的问题,结合微地震资料特点,通过理论分析与实际生产相结合,提出了单震相微地震事件识别、解域约束下微地震联合反演技术,并对相应的微地震有效事件进行处理验证,结果与实际资料相吻合.

2 微地震事件中单震相识别 2.1 单震相同相轴的判别

对于已识别出的有效微地震事件,并且只有一种震相情况即只有P波或S波,如何判断它是何种震相波?常规的能量比法、视速度法等都无法从中判断并识别出是S波还是P 波.因此结合微地震记录的特点,提出通过到时时差分析来识别P、S波震相的方法.

2.2 相邻道时差的变化规律

对于一个有效微地震事件记录,如果多余两个震相,即在所研究的记录范围内,出现两个或两个以上的同相轴,这种情况下微地震震相的识别,利用两个震相到时时差分析来识别到达波波场的性质,已经为加拿大的学者研究应用.该方法的基本思想是:假设2C为相邻两检波器之间的距离,通过2C/VP,SΔt比较分析可判断震相类型,其中2C为相邻道检波器的距离,VPVS 分别为P 波速度和S 波速度.例如,如果检测到的到达时差超过了2C/VS,且没有别的证据表明是外来杂波,则初至波一定为纵波,次至波是横波.

但是对于实际资料中大量存在的单一震相的有效微地震事件来说,上述方法并不适用,因此考虑到时差分析的优点,结合裂缝监测观测系统的特点,提出综合相邻道时差和首尾道时差对比的单一震相识别方法.

假定震源的位置为(x0z0),相邻道检波器的位置为确定的(xizi),(xi+1zi+1),其中zi+1 =zi+dii∈[1,n],则相邻道时差可以表示为Δt=l/v,时差的大小与v成反比,一般情况下,vp/vs ≈1.73,也就是说,其他条件确定的情况下,P 波的相邻道时差Δt要小于S波时差,因此对同一震源发出的且由同一检波器序列接收到的PS 波通过时差的对比可以识别,但是,值得注意的是,对任意单震相的微地震有效事件而言,(x0z0)并不已知,检波器的位置(xizi)对于不同的观测系统而言,也是不同的,且P、S波的时差对比,必须要已知P波或者S波的时差变化规律,显然上述方法不适用于所讨论的问题.

如式(1)所示,将Δt表示为xzziv的函数,Δt的取值与xzziv有关,单纯依靠速度差异对时差大小的影响难以对单一震相的微地震事件进行有效地识别.

(1)

其中,(x,z)表示为有效微地震事件的位置,zizi+1依次为任意相邻两道检波器的纵向深度,xi为检波器的横坐标,d为任意相邻道检波器间的纵向距离.

因此,通过对影响时差Δtxzxiziv等因素的讨论,计算出P 波、S 波的相邻道Δt时差的分布区间,同时利用得到的相邻道时差变化规律可以计算出P波、S波首、尾两级检波器的时差Δt′ 的分布区间,然后结合ΔtΔt′的分布规律识别单一震相的微地震事件.

现在,讨论Δt相对xzxiziv的变化规律,设定检波器为n级三分量检波器,垂直排列,检波器的坐标设为(xizi),任意微地震事件的坐标为(x,z),其中纵坐标z>>zixi=0,i=1,…,nx∈ (Xmin, Xmax),z∈ (Zmin, Zmax),则

(2)

(3)

(4)

对于式(2),当x∈ (Xmin, Xmax),恒小于0,说明时差Δt大小在径向上随着x值增大逐渐减小,即有Δtmin =f(Xmax),反之亦然.

对于式(3),当即在深度方向上,震源在某相邻两级检波器的中心位置上方时,上式(3)恒小于0,反之,当式(3)恒大于0.

假定事件深度在检波器排列下方,即z>zii=1,…,n.对于z∈ (Zmin, Zmax),Zmin 恒大于式(3)恒大于0,即Δtmin =f(Zmin),反之亦然;

对于式(4),同理我们可以得出,Δt随着zi的增大而减小,结合上述各种因素,相邻道时差Δt=f(XZzizi+1).

相邻道时差极大和极小值为:

(5)

(6)

依据式(5)、(6)可分别求出P 波对应的Δt区间,

(7)

(8)

同理可得S波时差区间:

(9)

(10)

同理可以得到P 波对应的首尾两级检波器对应的时差Δt′ 区间:

(11)

(12)

S波对应的首尾两级检波器对应的时差如式(13)、(14)所示:

(13)

(14)

因此,综合考虑震源的位置、检波器的位置以及P、S波速度差异条件下及已知压裂范围条件下,即[Xmin, Xmax],[Zmin, Zmax]区间限定后,可以计算出P、S波的相邻道时差ΔtPΔtS,首尾道时差ΔtPΔtS的分布区间.对于单一震相的微地震有效事件可以得到相应的ΔtΔt′,与计算的出ΔtPΔtSΔtPΔtS分布区间范围相比较,判断该事件为P波还是S波.

ΔtmaxPΔt′ ≤ ΔtmaxS, 或m-1个相邻道时差中存在ΔtmaxPΔtΔtmaxS, 如果满足以上条件为S波,否则为P波.

3 P、S波速度反演

P、S波反演,关键在于速度结构标定,尤其在只得到P波测井资料和射孔资料的条件下,如何准确地标定S波的速度结构,将直接影响定位精度.常规处理中,通常认为同一地层的P、S波速度间存在常数比例关系,约为1.73,但是对微地震反演而言,由于水力压裂产生事件的能量较弱,而且监测区域分布在几百米的范围内,因此,对速度结构的精度要求很高,而常规的近似处理根本无法保证其精度.本文提出利用射孔资料反演S波速度.首先,对射孔资料进行P、S波震相识别并分别拾取初至.然后利用测井声波时差资料建立P、S波的初始速度模型.针对拾取的P、S波的初至,利用射孔P、S波走时差来约束反演P、S波速度.

4 单震相微地震事件反演 4.1 正演算法

射线路径的计算方法大致分为三类[15-22]:试射法、弯曲法和迭代法.结合计算模型的特点,考虑算法各自的优缺点,试射法射线追踪方法适用于水平层状模型,其具有很快的收敛速度,且计算精度高,能够满足水平层状模型正演的计算精度要求,但是,当界面较为复杂时,试射法射线计算效率非常低,且路径搜索的收敛性比较差;弯曲法适用于比较复杂的地质情况,将复杂地质构造网格化,通过结点的选择来计算走时,但是其计算的精度与效率依赖于选择网格的大小,对基于模型正演微地震反演来说,网格选的过大计算就会受到影响;网格尺度选择得太小,计算效率就会大大降低.为了能够满足以下反演问题的要求,研究应用斜层模型下迭代法射线路径追踪算法.

4.2 解域约束下的遗传算法反演

实际微地震资料经过去噪、滤波、极化等一系列的处理后,拾取的初至与理论值已经产生了局部的差异,且常规的网格搜索法对初至变化的敏感度很高,直接导致真解的寻找变得十分困难,给实际反演实际微地震事件增加了不确定性.因此,考虑到网格搜索法简单、快速的特点,能够快速确定真解的分布区间,但是真解的确定又易受初至、速度模型等因素的影响,很难确定,而遗传算法对初至的适应程度很高,而且全局搜索能力强[23-25],但是任意设定搜索范围,则搜索过程会繁琐、费时,且解的可靠性大大降低,结合二者的优缺点,提出了解域约束下的遗传算法联合反演方法.

4.2.1 解域约束及其在遗传算法中的界定条件分析

以实际记录到的微地震事件并且震相只有S波为例.设定震源位置为(207.4 m, 2140.0 m),分析网格搜索法反演结果发现,图 1a所示,径向R方向上的解域误差分布呈现单调性,根本无法找到一个确切的真解,而在Z方向上的解域误差分布呈现多个收敛子区间,可见单纯依赖误差值的大小来评价真解会忽略掉许多关于解的其他信息,同时也无法准确判断真解的位置.但是通过分析反演解的分布规律后,发现解分布的概率密度函数区域分布特征中同样隐含着真解信息.如何通过分析反演解的误差分布规律,确定真解的定义域是研究问题的重点.具体思路是:

图 1 网格搜索法实际得到的解域(a)R方向解域分布特征;(b)Z方向解域分布特征. Fig. 1 The solution domain of using the method grid search (a) The solution domain in R direction;(b) The solution domain in Z direction.

首先,将所有解对应的误差值进行统计,考虑到大于误差期望值所对应的(R,Z)不可能为真解,即将小于误差期望所对应的(R,Z)筛选出来.然后,将筛选出的误差划分为多个子区间,将最大出现概率的误差子区间估计为真解的误差区间,实际反演过程中,真解及其附近的解对应的误差肯定落在真解误差区间内,而且占很大比例(图 2).最后,在选定真解误差区间的基础之上,对(R,Z)的分布进行概率统计,选择R、Z分布概率最大值对应的解分布区间,作为设定遗传算法反演参数设定的依据.

图 2 解域约束过程(a)R方向解域约束过程;(b)Z方向解域约束过程. Fig. 2 The progress of constrained solution domain (a) R direction (b) Z direction.

因此,可以将初始种群数目选为20,遗传代数为100代,将R、Z搜索范围分别设为(180.0 m, 230.0m),(2120.0m, 2170.0m)作为搜索范围,定义染色体的数目为2,进行反演计算,搜索结果R、Z的分布如图 3(a, b)所示.图 3中横坐标表示R的数值分布,纵坐标表示遗传运算过程中迭代计算对应的适应值,适应值越大,说明搜索到的解越可信,反演得到真解在(207.0m, 2143.0m),与模型设定的震源(207.4m, 2140.0m)基本一致.反演过程中真解寻找过程如图 4 所示,可见解域约束处理后,遗传反演搜索范围一直在真解的附近波动,有效地避免了奇异解的产生.

图 3 最优解分布特征(a)R方向最优解;(b)Z方向最优解. Fig. 3 The distribution of optimal solution (a) The optimal solution distribution of R;(b) The optimal solution distribution of Z.
图 4 搜索解域分布 Fig. 4 The distribution of searching solution domain
4.2.2 实际微地震资料反演

实际反演资料来自国内某地区水力压裂微地震资料.在其压裂后期,开启的裂缝一般都远离压裂井,离监测井较近,裂缝震源对应的微地震事件能量较弱,记录的信噪比也较低,且多为单震相微地震记录,在这种情况下,事件震相的正确识别非常关键.当S波事件误识为P 波情况下,裂缝发育趋势在靠近监测井一侧发生畸变(如图 5所示),利用本文研究的单震相识别技术对相应的微地震事件重新识别后,反演得到如图 6所示的结果,对比两次反演结果,可以发现震相识别后图中右侧裂缝分布明显收敛.通过分析该压裂区域构造特征,可以判断出主裂缝发育趋势应该为明显的狭长裂隙,而图 5中监测井侧的事件分布杂乱无序,且裂缝发育规律有明显地拓宽趋势,与实际情况明显不符.其次,从反演结果的聚散特征也能判断图 5 结果的右半部分是不合理的.

图 5 错误识别S波反演结果 Fig. 5 The inversion results of wrong identification of wave S
图 6 正确识别S波后反演结果 Fig. 6 The inversion results of right identification of wave S
5 结论

通过对影响时差微地震事件的位置检波器位置及速度等因素讨论,总结出上述诸因素对震相到时差的变化规律,根据变化规律,经数学公式的严谨推导,建立了震相时差与微地震事件、检波器空间位置及速度变化的定量关系.以此计算出震相的相邻道的时差和首尾两级检波器时差的分布区间,结合其分布规律识别单一震相的微地震事件.对实际微地震事件进行识别,不但识别出一些新的微地震事件,而且改进了原来的误识结果.

网格搜索法简单、快速,能够迅速确定真解的大概位置,但是真解的确定又易受初至、速度模型等因素的影响,很难确定,而遗传算法对初至的适应程度很高,而且全局搜索能力强,结合二者的优缺点提出的解域约束下的遗传算法联合反演方法,把寻解过程中的约束条件分析、解的搜索及解的评价问题综合为一种定量确定真解的分析方法,使得反演求取的解更加可靠.在微地震事件反演过程中,利用搜索法寻优时,纵向和径向方向上解的迭代误差分布具有不同的变化特征,径向上呈现单调变化趋势,而纵向上具有多个收敛子区间的变化特征.这些变化特征与接收排列与微地震事件源的相对空间位置密切相关,即蕴含着真解信息.通过研究发现,最可信解域是迭代误差期望概率密度大的区域.在迭代误差期望概率密度相同的情况下,真解区域应选择解为平坦或光滑解对应的区域.

综合微单震相识别和联合反演技术,对实际的微地震事件反演,反演结果表明:实际裂缝的分布和反演结果是一致的,同时反演结果聚敛效果较好.

致谢

论文研究得到了中石化科技部、中石化胜利地球物理勘探开发公司、中石化胜利油田分公司及中国石油大学物探实验室的大力支持,在此对其表示衷心地感谢.

参考文献
[1] 宋维琪, 陈泽东, 毛中华. 水力压裂裂缝微地震检测技术. 东营: 中国石油大学出版社, 2008 . Song W Q, Chen Z D, Mao Z H. Hydro-fracturing Break Micro-seismic Monitoring Technology (in Chinese). Dongying: Publishing House of China University of Petroleum, 2008 .
[2] Jean D, Frederique F. Multivariate statistical analyses applied to seismic facies recognition. Geophysics , 1988, 53(9): 1151-1159. DOI:10.1190/1.1442554
[3] Ivan D M, Jean J B, Bruce S H. A visual data-mining methodology for seismic facies analysis: Part 1-Testing and comparison with other unsupervised clustering method. Geophysics , 2009, 74(1): 1-11.
[4] Ivan D M, Jean J B, Bruce S H. A visual data-mining methodology for seismic facies analysis: Part 2-Application to 3D seismic data. Geophysics , 2009, 74(1): 13-23.
[5] Meersman K D, Kendall J M, Baan M. The 1998 Valhall microseismic data set: An integrated study of relocated sources, seismic multiples and S-wave splitting. Geophysics , 2009, 74(6): 183-195.
[6] Lisbeth E S. Inversion of Arrival times of micro-earthquake sources in the north sea using a 3-D velocity structure and prior information PART I. Method Bulletin of the Seismological Society of America , 1991, 81(4): 1183-1194.
[7] Andy J, John C, Rob J. Micro-seismic monitoring: Listen and see the reservoir. World Oil , 2000, 9: 171-173.
[8] Dyer B C, Jones R H, Cowles J F. Micro-seismic Survey of a North Sea reservoir. World Oil , 2002, 3: 74-77.
[9] Willis R B, Fontaine J, Paugh L. Geology and Geometry: A Review of Factors Affecting the Effectiveness of Hydraulic Fracturing. SPE Eastern Regional Meeting, 14-16 September: 97993-MS.
[10] Wolhart S L, Odegard C E. Micro-seismic Fracture Mapping Optimizes Development of Low-Permeability Sands of the Williams Fork Formation in the Piceance Basin. SPE Annual Technical Conference and Exhibition, 9-12 October: 95637-MS.
[11] 徐果明. 反演理论及其应用. 北京: 地震出版社, 2003 . Xu G M. Inversion Theory and Applications (in Chinese). Beijing: Seismological Press, 2003 .
[12] 高尔根, 徐果明, 赵焱. 一种任意界面的逐段迭代射线追踪方法. 石油地球物理勘探 , 1998, 33(1): 54–60. Gao E G, Xu G M, Zhao Y. Segmental iterative ray tracing method for any interface. Oil Geophysical Prospecting (in Chinese) (in Chinese) , 1998, 33(1): 54-60.
[13] 万永革, 李鸿吉. 遗传算法在确定震源位置中的应用. 地震地磁观测与研究 , 1995, 16(6): 1–7. Wan Y G, Li H J. The preliminary study on the seismic hypocenter location using genetic algorithms. Seismological and Geomagnetic Observation and Research (in Chinese) (in Chinese) , 1995, 16(6): 1-7.
[14] 宋维琪, 刘军, 陈伟. 改进射线追踪算法的微震源反演. 物探与化探 , 2008, 32(3): 275–277. Song W Q, Liu J, Chen W. Micro-earthquake source inversion of an improved ray tracing algorithm. Geophysical and Geochemical Exploration (in Chinese) (in Chinese) , 2008, 32(3): 275-277.
[15] 王华忠, 方正茂, 徐兆涛, 等. 地震波旅行时计算. 石油地球物理勘探 , 1999, 34(2): 155–163. Wang H Z, Fang Z M, Xu Z T, et al. Computation of seismic travel time. Oil Geophysical Prospecting (in Chinese) (in Chinese) , 1999, 34(2): 155-163.
[16] Moser T J. Shortest path calculation of seismic rays. Geophysics , 1991, 56(1): 59-67. DOI:10.1190/1.1442958
[17] Leidenfrost A, Ettrich N, Gajewski D, et al. Comparison of six different methods for calculating travel times. Geophysical Prospecting , 1999, 48(4): 269-297.
[18] Chander R. On tracing rays with specified end points. Geophys , 1975, 41: 173-177.
[19] Chander R. On tracing seismic rays with specified end points in layers of constant velocity and plane interface. Geophys Prosp , 1977, 25: 120-124. DOI:10.1111/gpr.1977.25.issue-1
[20] Cerveny V, Molotkov I A, Psencik I. Ray method in seismology. 1977 .
[21] Farra V. Ray tracing in complex media. Applied Geophysics , 1993, 30: 55-73. DOI:10.1016/0926-9851(93)90018-T
[22] Julian B R, Gubbins D. Three dimensional seismic ray tracing. Geophys , 1977, 43: 95-114.
[23] Dong H P, John A Q, Bruce E C, et al. Velocity calibration for micro-seismic monitoring: A very fast simulated annealing (VFSA) approach for joint-objective optimization. Geophysics , 2009, 74(6): 47-55.
[24] Raghu K C, Mrinal K S, Paul L S. Hybrid optimization methods for geophysical inversion. Geophysics , 1997, 62(4): 1196-1207. DOI:10.1190/1.1444220
[25] Paul L S, Mrinal K S. Nonlinear multi-parameter optimization using genetic algorithms: Inversion of plane-wave seismograms. Geophysics , 1991, 56(11): 1794-1810. DOI:10.1190/1.1442992