2. 中国科学院大学, 北京 100049
2. University of Chinese Academy of Sciences, Beijing 100049, China
页岩气和致密油气的开发是当前油气勘探开发的热点.压裂改造是这类油气开发的主要技术手段(Agarwal et al., 1979).在储层压裂改造过程中,储层受压破裂形成微地震信号.一次破裂形成的微地震信号被称为一个微震事件.获取尽可能多的微震事件,计算其震源位置,利用震源间的空间位置关系,可对压裂裂缝的发育方向等信息进行判断.因此,微地震监测是对压裂作业进行直观评价的有效方法(Daniels et al., 2007; Maxwell et al., 2010).
微震事件具有两个主要特点:它是一种低能量高频率信号(Warpinski,1994);相对于持续数小时的压裂监测过程,单个微震事件持续时间很短,仅有几十毫秒.这些特点导致微震监测系统的采样频率较高,通常为2 KHz或4 KHz,记录得到的数据量很大.因此,微震监测的首要工作就是在大数据量的监测数据中识别出尽可能多的微震事件,压制相应记录中的噪音干扰,提取事件初至和极化特征等信息(Tselentis et al., 2012).前人根据微震事件的各种信号特征提出了多种微震事件的判别方法.其中,STA/LTA(Short Time Average/Long Time Average)是常见的微震事件识别方法,它利用地震记录波形的能量变化作为判断标准;Lois等(2012)提出利用微震事件的极化特征进行判别;St-Onge(2010)提出利用AIC(Akaike Information Criteria)来作为判别标准;刘劲松等(2013)利用改进的高阶累计量的偏度和峰度计算方法来识别微震事件.以上方法可提高数据处理的效率,减轻人工处理的负担.然而上述方法仅适合处理信噪比较高的微震事件,当监测数据信噪比过低时,上述方法易出现微震事件的误判.由于微震事件的强度决定于压裂时注入流体的压力、体积以及被压裂地层性质等因素,并极易受到各种噪音的干扰.大多数情况下,监测数据中的微震事件信噪比较低.对于低信噪比的微震事件,即使采用人工方式也较难识别,当采用上述方法时,事件识别的准确度降低,同时能识别出的事件数量相当有限,进而加大了后期微震监测解释的难度.
本文针对在人工方式下较难识别的低信噪比微震记录,利用时频稀疏性方法自动判别微震事件是否存在及存在的时段.本文研究表明,将微震监测记录变换到时频域后,噪音能量被分散,由于微震事件在时频域中的稀疏性,使其能量相对增强.基于以上对微震事件特点的分析,本文提出一种新的微震事件自动识别方法,该方法利用Renyi信息熵计算各个数据道在时频域内的稀疏性,熵值越小表示相应道的微震事件存在的可能性越大,挑选多个低熵值的数据道,对这些数据道在时频域中的最大稀疏系数(假定为微震事件)在时频域平面内做聚类分析,若满足设定的约束条件,即可认为当前记录中有微震事件存在.同时,由聚类分析所得到的中点位置可视为微震事件时频特征的估计值,压制估计区域外的时频数据,对各个数据道做滤波处理,即可恢复出较高信噪比的微震事件. 2 时频稀疏性分析法
时频分析方法最早是Gabor(1946)提出的短时 傅立叶变换STFT(Short Time Fourier Transform),这种方法先利用滑动时窗将信号截成若干小段,再逐段做傅立叶变换,即可确定每个时段的频率特征.小波变换克服STFT分辨率固定的不足,具有可变的分辨率,但小波变换不能很好保留信号的相位信息,另一方面,小波的尺度与信号的频率之间不存在准确的对应关系.Wigner分布具有最为理想的时频分辨率(Cohen,1995),但由于这种方法是非线性分析方法,会产生交叉项干扰,时域信号的相位信息在时频域中完全被消除.S变换(Stockwell et al., 1996)是短时傅里叶变换和小波变换的结合体,其窗函数是一个尺度随频率变化而变化的高斯窗(Pinnegar and Mansinha, 2003).虽然它不具有Wigner分布的高时频分辨率,但它不会引入交叉项干扰,并且保留被分析信号各个频率成份的相位信息,并且S变换完全可逆.高静怀等(2003,2004)对S变换进行了推广,并将其定义的广义S变换应用于薄互层检测和信号识别中.结合微震事件的时频特征,本文采用S变换作为时频分析工具. 2.1 稀疏性测量方法
有多种方法可对时频分布的稀疏性做出定量计算:范数比测量法(Jones and Parks, 1992),该方法利用时频分析结果的不同范数的比值作为判断依据,范数比值越高,时频分析结果的稀疏性越好;Stankovi čc'(2001)测量法计算结果的数值越小,表明时频分析结果越稀疏;Shannon(2001)信息熵测量法是对不确定性的一种度量.如果时频分析结果的稀疏性越好,不确定性越小,相应的熵值越小.但若时频分析中出现了负值时,这种方法就不再适用;Renyi(1961)信息熵测量法(见公式(1))克服了Shannon信息熵测量法的不足,同时Renyi熵法在计算前,对时频分析结果做归一化处理,将原始时频分析结果转化为概率密度函数为
其中R为Renyi熵值,t为时间,f为频率,p(t,f)为时域信号时频表示的归一化结果,即=1,α 为控制系数.
由于各个检波器与井管耦合程度存在差异,微震事件按照不同的传播路径到达各个检波器,这些因素会造成同一微震事件在不同检波器上的地震记录信噪比存在明显差异.因此,对于监测数据做自动处理,单道的特征并不充分可靠,需结合多道数据的信息,以防止出现误判.本文采用Renyi熵值表示单道监测数据的时频稀疏程度,这样可消除各检波器接收信号的能量差异,联合多道信息对微震事件进行自动判别.结合Fland rin等(1994)对Renyi熵参数α的分析结果,将这个控制系数设定为3. 2.2 稀疏约束条件
微震事件到达各个检波器的时差主要受观测系统和地层速度的限制,通常这两点在压裂改造前为已知信息或可做近似估计.同时,虽然各个检波器记录到的微震事件信号主频有所不同,但这种差异有限.因此在时频域平面上,各数据道时频分析的最大稀疏系数(假定为微震事件)的位置应聚集在一起,即同一微震事件在各道间的时频距离总是在特定的范围内:
其中T表示微震事件记录的单道数据,Ti和Tj为任意不同两道,ΔfTi,Tj是微震事件在Ti和Tj两道中的主频频率差,ΔtTi,Tj是微震事件到达Ti和Tj 两道的时差,Cf和Ct为设定的频域和时域约束条件. 2.3 判别微震事件的目标函数计算各个数据道的熵值R,同时为防止在单个数据道中出现的随机强噪音对分析结果的干扰,将背景噪音熵值的均值作为门槛值Rnoise.如果数据道的熵值R小于门槛值Rnoise,则这个数据道被认为可能含有微震事件.低熵值的数据道越多,微震事件存在的可能性越大.
建立微震事件判别的目标函数为
其中T代表单个数据道,NT为小于门槛值Rnoise的低熵值数据道数量,m为事件判别所需的最小道数,m的具体数值大小可由后期微震震源定位所使用方法特点来决定,RT为单个数据道的Renyi熵值,Rnoise为背景噪音的Renyi熵值,Rnoise可在压裂改造前的背景记录中计算获得,公式(3)的第二部分为公式(2)定义的约束条件,s.t.表示约束条件.对于不含微震事件的记录,由于噪音是随机的,在时频域上,最大稀疏系数出现的位置也是随机的.这样,各数据道时频分析结果中的最大稀疏系数在时频域中是分散的,不能满足公式(2)的约束条件.本文利用聚类分析来实现公式(3)的约束判别,用来寻找多数数据道聚集的区域.采用聚类分析方法有两个原因:第一是可降低强噪音的干扰;第二是可减少人工干预.
如果当前监测记录中有足够多的低熵值数据道,并且这些数据道的最大时频稀疏系数点之间满足频差和时差约束条件,即可判定当前监测记录中有微震事件存在.同时,将这些系数点的中心位置作为微震事件主频和初至位置的估计值. 2.4 合成数据分析
合成数据为两端进行指数衰减处理的正弦信号,以此来模拟微震事件.合成信号长度为400个样点,采样频率200 Hz,合成信号的主频为80 Hz.图 1a为合成信号的时域波形,图 1b为相应时频分析结果(横轴为频率,纵轴为时间,色标表示时频点的能量大小),合成信号在时频域中具有明显的稀疏性.在原始合成数据的基础上,加入较强的高斯白噪音干扰(见图 1c),从图 1c中可以看出,在时域波形中,很难对是否存在微震事件做出判断.将图 1c的时域波形变换到时频域中(见图 1d),其时频域最大能量点出现的位置与图 1b保持一致.这一数值计算结果表明,在强噪音背景下,微地震信号完全淹没在噪音中,微地震记录表现出信噪比极低的特征,对于这样的低信噪比微地震记录,几乎无法在时域采用人工方式或计算机自动检测方法判别微震事件是否存在.但是,通过时频稀疏性分析,可在时频域中对其进行判别.
本文选择某一砂岩储层压裂改造的微震监测数据,利用本文提出的方法对微震事件进行自动识别与恢复.观测系统采用井下临近井观测方式,10级 三分量检波器,即30个数据道,检波器的间距为10 m,采样频率为2000 Hz.此砂岩储层的纵波速度均值为4500 m·s-1,横波速度均值为2500 m·s-1.在压裂改造前,对压裂井的相应井段进行射孔处理.在压裂改造过程中,由于微震监测设备和工区环境等多方面原因,实际接收到的微震事件信噪比较低.图 2是该微震监测数据中的一段典型记录,从中可看出监测记录存在很强的背景噪音干扰,微震事件的波形不清晰,并且中间几级检波器的接收信号较弱.简单对比各道数据的波形,只能大致看到第3000个采样点,第5200个采样点和第6500个采样点附近可能有微震事件存在.
图 3中的三组记录分别是此次微震监测中的射孔记录,背景噪音记录和微震监测记录(对应图 2中 的第9233个到第9632个采样点).射孔记录(图 3a)的信噪比较高,射孔事件记录较为清晰;微震监测记录(图 3c)信噪比较低,几乎与背景噪音记录(图 3b)无差异,很难在时域中判断微震事件是否存在.
对于同一检波器的不同分量间往往具有接近的时频稀疏特征,所以在本文中并不区分数据道的分量类型.在图 3a、3b和3c中分别抽取一道数据做时频分析(见图 4).对于射孔记录(图 4a),由于射孔事件能量强,记录的波形信噪比较高,时频域内的能量相对集中(图 4d);对于背景噪音记录(图 4b),其时频域内的能量比较分散(图 4e);对于微震监测记录(图 4c),其时域波形很难分辨出有微震事件存在,而时频域内的能量较为集中,稀疏性较好(图 4f).微震监测记录的时频特征(图 4f)与背景噪音记录的时频特征(图 4e)存在明显的差异,而与射孔记录的时频特征(图 4d)具有一定的相似性.虽然图 4f中的时频稀疏性明显,但还是有多个能量强点出现,仅以这一道的时频特征还不能对微震事件存在是否存在做出判别.
对微震监测记录(图 3c)各道均做时频分析.由 于采集设备等方面原因,原始记录波形特征(见图 2)具有如下特点:两端数据道的资料品质相对较好,而中间数据道的资料品质相对较差.因此抽取第2 级、第3级、第6级和第9级三分量检波器(即第4—6道、 第7—9道、第16—17道和第25—27道)的波形记录和时频分析结果进行显示(见图 5).从图 5中可知,各道时域波形中的噪音较强,很难识别出微震事件.而各道的时频稀疏性存在明显差异:其中图 5b7(微震监测记录中的第16道)的时频稀疏性相对较差,图中有多个能量强点出现,而且出现位置较分散;图 5b1、5b2、5b3、5b6和5b11(分别是微震监测 记录中的第4道,第5道,第6道,第9道和第26道)的时频稀疏性相对较好.综合各个数据道的稀疏性特征,对微震事件是否存在进行判别.判别过程如下:
(1)提取时频稀疏性较好的数据道
利用公式(1),计算出各数据道时频分析的熵值R.从图 5的时频分析结果中可以看出,各道都存在较强的低频噪音,在计算熵值R前,需要对这些低频噪音进行压制.图 6是对图 3中射孔记录、背景噪音记录和微震监测记录的各道时频分析的熵值R计算结果.在图 6中,背景噪音记录的熵值线(紫色)为高值,各道间的熵值保持相对稳定;射孔记录的熵值线(蓝线)明显为低值;微震监测记录的熵值线(绿线)几乎与背景噪音记录的熵值线相重叠,仅在第一级、第二级和第三级检波器上存在明显的差异.如前所述,在原始监测记录中,两端的数据道的资料品质相对较好,中间数据道的资料品质相对较差,数据道的时频稀疏程度测量结果(Renyi熵)也体现了这个数据特点(见图 6).
由于背景噪音记录的熵值相对稳定,因此可以其熵值的均值作为门槛值Rnoise.本实例采用的门槛值Rnoise为15.5.由图 6可知,微震监测记录(图 6绿线)共有14个数据道的熵值小于此门槛值Rnoise,分 别是第2道,第3道,第4道,第5道,第6道,第7道,第8道,第9道,第15道,第17道,第21道,第23道,第26道和第27道.
(2)建立时频约束条件
由本实例的观测系统特点可知,各个检波器间的理论最大距离差为90 m,即当整条测线为直线时,第一级检波器到第十级检波器的距离.低信噪比条件下,由于微震事件纵波部分不易被记录到,仅能记录其横波部分,则在各个检波器间微震事件初至时差理论最大值为36 ms,即72个采样点.这样,时间方向的约束值可设定为72个样点数,对应为公式2中的Ct.由于在压裂改造前,无法提前确定出微震事件的主频范围,所以对频率方向的约束值设定较为宽松.在本实例中,频率约束设定为200 Hz,对应为公式2中的Cf.
(3)多道时频稀疏性的聚类分析
将上述14个低熵值数据道的时频分析最大稀疏系数投影到同一平面中(见图 7).在图 7中,数字表示对应的道号,红点表示相应道的最大稀疏系数的位置.由聚类分析方法可自动得到:第2道,第3道,第6道,第7道,第21道,第23道,第26道和第27道的最大稀疏系数的位置最为聚集.同时,这些数据道的最大稀疏系数之间的位置关系满足设定的时频约束条件要求.聚类分析有效排除了异常点的干扰,如第4道,第9道和第17道.因此可判定图 3c的记录中有微震事件存在,并以聚集区内系数点的坐标均值作为此微震事件的时频特征(事件主频和初至位置)估计值.
沿时间方向,逐段抽取图 2中的监测数据,经过微震事件自动识别处理后,共找到7个微震事件,分 别位于第2200个,第3010个,第4800个,第5270个,第6500个,第8030个和第9400个采样点附近.本文提出的自动识别方法的分析结果不仅与人工识别结果相一致,而且识别出更多在时域内不易判别的微震事件.
(4)微震事件恢复
因为不同的微震事件对应的震源破裂规模和破裂机制存在差别,同时不同震源点到各个检波器的传播距离也有差异,进而震源信号在传播过程中会存在不同程度的衰减,所以不同震源信号的原始主频会有差别,检波器接收到的微震事件的主频也存在差异.利用聚类分析可自动对微震事件时频特征进行估计,减少人工干预.利用时频特征的估计值,可对原始的微震监测数据进行滤波处理,恢复出信噪比较高的微震事件波形.图 8a是原始微震监测记录(见图 2)中的第6400~6620个采样点的波形.虽然人工方式能判断出第6500个采样点附近有微震事件存在,但这段记录的信噪比较低,微震事件波形和初至位置依然不够清晰.采用本文提出的方法进行事件自动识别和恢复后(见图 8b),微震事件的波形特征和初至位置得到明显改善.图 9a是原始微震监测记录中的第9340~9560个采样点的波形,相对图 8a中的波形记录,其信噪比更低,人工方式不能直接对记录中是否存在微震事件做出判断,而使用本文提出的方法,可判断出在这段记录中有微震事件存在,恢复的微震事件波形见图 9a.相对图 9a中的波形,虽然图 9b中恢复的微震事件信噪比相对较低,但恢复的微震事件波形特征相对清晰.恢复微震事件以后,可采用前人提出的各种方法来获取微震事件的初至和极化特征等信息.
本文提出利用Renyi熵值来表示微震监测记录的稀疏程度,并以计算结果为基础,建立了带约束条件的微震事件判别目标函数.它能针对低信噪比的微震监测数据自动判别出微震事件是否存在,并能对微震事件主频率特征和初至位置做出估计.利用这个估计值,通过滤波处理提高微震事件信噪比,为下一步的微震监测资料处理提供准确的数据信息.由于本文提出的方法几乎不需要的人工干预,大大提高数据处理效率,可实现微震监测数据的自动处理.同时,本文提出的事件自动识别方法也可应用于地球物理的其它领域,如天然地震事件的自动识别.
[1] | Agarwal R G, Carter R D, Pollock C B. 1979. Evaluation and performance prediction of low-permeability gas wells stimulated by massive hydraulic fracturing. Journal of Petroleum Technology, 31(3): 362-372, doi:10. 2118/6838-PA. |
[2] | Cohen L. 1995. Time-frequency analysis. New Jersey: Prentice Hall PTR. |
[3] | Daniels J, Waters G, Le Calvez J, et al. 2007. Contacting more of the barnett shale through an integration of real-time microseismic monitoring, petrophysics, and hydraulic fracture design. In SPE Annual Technical Conference and Exhibition, doi:10. 2523/110562-MS. |
[4] | Flandrin P, Baraniuk R G, Michel O. 1994. Time-frequency complexity and information. Proc. IEEE Int. Conf. Acoust., Speech, Signal Processing, 3: 329-332, doi:10.1109/ICASSP.1994. 390031. |
[5] | Gabor D. 1946. Theory of communication. Part 1: The analysis of information. J IEE, 93(26): 429-441, doi: 10.1049/ji-3-2.1946. 0074. |
[6] | Gao J H, Chen W C, Li Y M, et al. 2003. Generalized S transform and seismic response analysis of thin interbeds. Chinese J. Geophys. (in Chinese), 46(4): 526-532. |
[7] | Gao J H, Man W S, Chen S M. 2004. Recognition of signals from colored noise background in generalized S-transformation domain. Chinese J. Geophys. (in Chinese), 47(5): 869-875. |
[8] | Jones D L, Parks T W. 1992. A resolution comparison of several time-frequency representations. IEEE Transactions on Signal Processing, 40(2): 413-420, doi:10.1109/78. 124951. |
[9] | Liu J S, Wang Y, Yao Z X. 2013. On micro-seismic first arrival identification: A case study. Chinese J. Geophys. (in Chinese), 56(5): 1660-1666, doi:10.6038/cjg20130523. |
[10] | Lois A, Sokos, E, Martakis N, et al. 2012. A new automatic S-onset detection technique: Application in local earthquake data. Geophysics, 78(1): P. KS1-KS11, doi:10.1190/GEO2012-0050. 1. |
[11] | 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. |
[12] | Pinnegar C R, Mansinha L. 2003. The S-transform with windows of arbitrary and varying shape. Geophysics, 68(1): 381-385, doi:10.1190/1. 1543223. |
[13] | Renyi A. 1961. On measures of entropy and information. Fourth Berkeley Symposium on Mathematical Statistics and Probability, 547-561. |
[14] | Shannon C E. 2001. A mathematical theory of communication. ACM SIGMOBILE Mobile Computing and Communications Review, 5(1): 3-55, doi:10.1145/584091. 584093. |
[15] | Stankovičc' L. 2001. A measure of some time-frequency distributions concentration. Signal Processing, 81(3): 621-631, doi:10. 1016/S0165-1684(00)00236-X. |
[16] | Stockwell R G, Mansinha L, Lowe R P. 1996. Localization of the complex spectrum: the S transform. IEEE Transactions on Signal Processing, 44(4): 998-1001, doi:10.1109/78. 492555. |
[17] | St-Onge A. 2010. Akaike information criterion applied to detecting first arrival times on microseismic data. GOPH701 project report, Department of Geoscience, University of Calgary, doi:10.1190/1. 3627522. |
[18] | Tselentis G A, Martakis N, Paraskevopoulos P, et al. 2012. Strategy for automated analysis of passive microseismic data based on S-transform, Otsu's thresholding, and higher order statistics. Geophysics, 77(6): KS43-KS54, doi:10.1190/Geo2011-0301. 1. |
[19] | Warpinski N R. 1994. Interpretation of hydraulic fracture mapping experiments. In University of Tulsa Centennial Petroleum Engineering Symposium, 291-300, doi:10. 2118/27985-MS. |
[20] | 高静怀, 陈文超, 李幼铭等. 2003. 广义S变换与薄互层地震响应分析. 地球物理学报, 46(4): 526-532. |
[21] | 高静怀, 满蔚仕, 陈树民. 2004. 广义S变换域有色噪音与信号识别方法. 地球物理学报, 47(5): 869-875. |
[22] | 刘劲松, 王赟, 姚振兴. 2013. 微地震信号到时自动拾取方法. 地球物理学报, 56(5): 1660-1666,doi:10. 6038/cjg20130523. |