CN106646598A - 一种fast‑aic法微地震信号拾取方法 - Google Patents

Abstract

本发明公开了一种基于FAST‑AIC法的微地震信号初至自动拾取方法：先根据水力压裂射孔信号与微地震信号传播路径相似，走时信息相近，提出了一种微地震事件识别新方法；利用水力压裂检波器布阵方式及单道信号局部极值在微地震触发时刻附近设计精确时窗；用Curvelet变换对精确时窗内数据滤波；通过FAST‑AIC算法在精确时窗内进行初至拾取。FAST‑AIC算法是将传统AIC算法优化、变换，求取离散实数序列的算术和、平方和的线性组合形式,减少重复计算,极大地提高计算效率。采用新方法对模拟数据和实际水力压裂微地震监测数据进行初至拾取，验证新方法的可行性和正确性。同时与常规微地震方法初至拾取结果对比，证明提出的新方法较常规方法具有明显的优越性。

Description

一种FAST-AIC法微地震信号拾取方法

技术领域：

本发明涉及一种微地震信号自动拾取方法，尤其是水力压裂微地震勘探中基于FAST-AIC法的微地震信号初至自动拾取方法。

背景技术：

微地震震源定位技术是水力压裂技术的核心，精确获取微地震信号初至是提高微地震震源定位精度的关键。通常情况下，地面观测的微地震资料数据量大、信噪比低、甚至完全淹没在噪声之中，这在很大程度上制约了微地震信号的初至拾取精度。准确地识别微震事件为高精度初至拾取提供了基础，但由于水力压裂震源发生时刻、震源数量的不确定性，致使准确地识别微地震事件也十分困难。

CN104459789A公开了《一种地震波初至拾取的方法》，根据地震剖面选取直达波，作为参考地震子波；参考地震子波与地震数据进行互信息量的计算；在整个时间轴上搜索互信息量极大值，此极大值所对应的时间即为所求的此道地震波初至。其方法包括根据地震剖面选取直达波，作为参考地震子波；地震子波与地震数据进行互信息量的计算；搜索得到互信息量的极大值，此极大值所对应的时间即为所求的地震波初至。本发明的地震波初至拾取的方法所得结果精度高，在低信噪比地区也能取得很好的结果

CN103837891A公开了《一种微地震初至的高精度拾取方法》，利用高阶矩方法增强微地震资料的信噪比，根据当前微地震波是否存在满足振幅大小、空间相关性、同相轴视视速度大小和偏振特征等约束条件来判断该微地震波是否为微地震波初至；根据射孔资料实际事件资料的变化规律，采用统计分析方法，确定约束条件的门槛值；在振幅和视速度约束下，围绕微地震事件走时初至的时空变化特征，利用固定大时窗和滑动小时窗内走时初至点能量、视速度、相关信息的变化规律，采用多次迭代的方法，进行微地震时间初至的自动拾取。

2011年宋维琪，吕世超在《基于小波分析与Akaike信息准则的微地震初至拾取方法》中根据微地震信号与噪声在小波多尺度分解后相邻尺度上连续性不同的特性，提出了识别微地震有效事件的小波-AIC方法。2013年魏路路在《水力压裂微地震事件自动识别与拾取》中联合AIC法和相关法用于水力压裂射孔信号初至波自动拾取，2013年张唤兰在《微地震数值模拟及震源定位方法研究》中提出基于时窗能量比和AIC算法的两步法拾取微震初至到时。

上述现有技术虽在微地震信号初至拾取中取得一定有益效果，但拾取过程复杂，计算效率较低，同时对于微地震事件识别都是基于单道信号处理之上的，忽略了各道之间的相关性，容易造成识别中误检或漏检的问题。

发明内容：

本发明的目的就在于针对以上微地震事件识别与初至拾取方法中存在的问题,首先根据水力压裂射孔信号与微地震信号传播路径相似，走时信息相近，提出了一种微地震事件识别新方法；其次，利用水力压裂检波器布阵方式及单道信号局部极值在微地震触发时刻附近设计精确时窗；本发明对传统AIC算法进行优化,将公式变换为求取离散实数序列的算术和、平方和的线性组合形式,减少重复计算,极大地提高计算效率。最终形成一套系统的基于FAST-AIC法的微地震信号初至拾取方法。

本发明是一种基于FAST-AIC法的微地震信号初至自动拾取方法，包括以下步骤：

(1)输入微地震监测数据以及射孔信号；

(2)人工拾取每条测线上检波器记录的射孔信号初至时刻，设置每条测线靠近射孔点的检波器记录信号为参考道，计算射孔信号初至时差文件Δt_i，即使每条测线上检波器记录的射孔信号初至时刻与该条侧线上参考道的初至时刻的差值；

(3)利用射孔信号初至时差文件Δt_i校正微地震数据；

(4)将校正后的多道微地震数据叠加获得一个模型道M(i)；

(5)对模型道M(i)做希尔伯特(Hilbert)变换求取能量瞬时包络E(i)；

其中为模型道M(i)的Hilbert变换

(6)采用STA/LTA法基于模型道能量包络E(i)识别微地震事件，当

STA/LTA比值超过阈值Thr时，认为检测到微地震事件，并且此时对应的时刻T_i记为微地震触发时刻，STA/LTA法定义如下：

其中，E_i代表微地震信号，m表示短时窗长度，n表示长时窗长度,N表示模型道长度,Thr为根据实际数据调节的经验值。

(7)采用检波器布阵方式在触发时刻T_i处设计初步时窗。

首先，分析震源S相对测线L1最近及最远距离位置(震源S到检波器A₁与An走时差的最大值与最小值)，即当公式(3)Δe取最大值和最小值时分别对应的震源坐标。

其中

其中(x₁,y₁，z₁)为检波器A₁坐标,(x_n,y_n,z_n)为A_n坐标,S为震源分布球体中一点，坐标为(x,y,z),分别代表震源点S到检波器A_n、A₁的距离。

由于震源点S在以射孔点(0,0,h)为圆心，半径为R的球形内，可知震源点S满足以下约束条件：

0≤x²+y²+(z-h)²≤R² (6)

根据约束条件公式(6)，联合公式(3)-(5)可求得Δe取最大值对应震源S₂坐标和最小值时分别对应的震源S₁坐标。

其次，根据两个极限震源坐标定义初步时窗为[T_i-Δ_i1,T_i+Δ_i2]。

分别代表他们两点间的距离，v代表地层速度。

(8)根据初步时窗内各道信号最大值确定精确时窗。

精确时窗获取包括以下步骤：

a、寻找时窗[T_i-Δ_i1,T_i+Δ_i2]内各道信号最大值Max_i；

b、更改时窗下限为各道信号幅值最大值对应时刻之后C1s处；

c、在各道时窗内最大值对应时刻向前C2s作为时窗起始点，即精确时窗范围为[Max_i-C2,Max_i+C1],其中C1、C2为经验值。

(9)利用离散Curvelet变化对时窗有效区域进行多尺度分解，将不同频率的信号剥离开；

(10)采用FAST-AIC法算法在精确时窗内进行准确初至拾取。

为提高计算效率，在传统AIC方法基础上进行改进,提出FAST-AIC法，将传统AIC法如公式(9)所示：

AIC(k)＝k×lg(Var(x[1,k]))+(N-k-1)lg(Var(x[k+1,N])) (9)

将公式(1)变换为求取离散实数序列的算术和、平方和的线性组合形式得到公式(10)和公式(11)。

其中并且设当地震数据长度N固定时,sm即为定值因此:

设因此：

式中x_i为微地震数据，k的取值对应x_i中采样点，N为微地震记录长度。采用(13)式计算的优点是,当地震波长度N固定时，sm、为定值，只需要计算一次，不随时窗分界点k的移动而改变。当k移动到新的位置时，只需计算新的sk1,sk2,并且在计算sk1，sk2的过程中不必重复计算时窗内所有点的值,只须在上一个时窗值的基础上,加上新进入时窗采样点的值或值的平方即可，大大减少计算量。

有益效果：本发明采用微地震事件识别和波至到时拾取两步实现微地震信号初至自动拾取，实现了多道微地震记录准确识别微地震事件，有效地解决了单道微地震事件识别中出现漏检或误检的问题。同时，通过对局部数据滤波和快速算法提高了初至拾取精度和计算效率，且具有良好的抗噪性能。

附图说明：

图1监测系统检波器布阵方式示意图。

图2一种FAST-AIC法微地震信号初至到时自动拾取方法流程图。

图3利用射孔信号校正微地震数据示意图。

图4极限震源S1、S2模拟数据校正图。

图5传统AIC方法与改进快速AIC方法计算效率对比图。

图6模拟数据DATA1道集示意图。

图7数据DATA1校正结果图。

图8微地震事件识别过程示意图。

图9精确时窗内数据放大图。

图10精确时窗内数据滤波后结果图以及拾取结果图。

图11现有STA/LTA法识别微地震事件结果图。

图12 FAST-AIC法初至拾取结果图。

图13实际野外检波器布阵图。

图14山西某地实际微地震信号记录图。

图15根据射孔信号初至时差文件校正微地震数据结果图。

图16微地震事件识别过程示意图。

图17不同方法对实际为地震记录初至到时拾取拾取结果图。

图18实际数据信噪比与初至拾取误差关系散点图。

具体实施方式：

为了使本技术领域的人员更好地理解本发明中的技术方案,本发明的检波器布阵方式采用“米”字形，如图1所示，且考虑到实际地层变化不一等因素，每条测线单独拾取，一次拾取一条测线上多道微地震信号初至。

下面结合说明书附图对本发明提出的一种FAST-AIC法微地震信号初至自动拾取方法进行详细说明，如图2所示，包括以下步骤：

(1)微地震数据及射孔信号输入；

(2)利用已知高信噪比射孔信号的初至时差文件Δt_i(各道信号初至时刻与参考道初至时刻的差值)对整体微地震数据静校正，使单条测线上各检波器记录的微地震信号初至到时基本一致。图3为利用射孔信号校正微地震数据示意图。图4为极限震源S1、S2模拟数据校正图。

(3)将校正后的单条侧线上多道微地震数据叠加获得一个高信噪比模型道M(i)；由于随机噪声没有规律，叠加后随机噪声减弱，而有效信号加强，极大地增加了模型道M(i)信噪比。

(4)为进一步克服水力压裂微地震信号低信噪比的特点，保证微地震事件识别的准确性，对模型道M(i)做希尔伯特变换以求取能量瞬时包络E(i)，使其幅值突出更加明显，有利于微地震事件识别。

其中为模型道M(i)的Hilbert变换。

(5)采用STA/LTA法基于模型道能量包络E(i)识别微地震事件。当STA/LTA比值超过阈值Thr时，认为检测到微地震事件，并且此时对应的时刻T_i为微地震触发时刻。STA/LTA法定义如下：

其中，E_i代表微地震信号，m表示短时窗长度，n表示长时窗长度,N表示模型道长度。Thr为经验值，根据实际微地震信号调整。

(6)采用检波器布阵方式先设计粗时窗，再在粗时窗内进一步获取精确时窗。

①粗时窗起止范围为[T_i-Δ_i1,T_i+Δ_i2]，其中：

t₁为粗时窗起始时刻，T_i为STA/LTA法识别微地震事件触发时刻，Δ1为震源点S₂到检波器A₁与检波器A_n的走时差和震源点S₀到检波器A₁与检波器A_n的走时差的差值。分别代表他们两点间的距离，v代表地层速度，可根据测井资料获得。

分别代表他们两点间的距离，v代表地层速度，可根据测井资料获得。

②精确时窗范围设计如下，首先寻找时窗[t₁,T_i+Δ2]内各道信号最大值Max_i，更改时窗下限为各道信号幅值最大值对应时刻之后C1s处,保证幅值最值处波形完全包含在精确时窗内，有利于计算精度的提高，在各道时窗内最大值对应时刻向前C2s作为时窗起始点，即精确时窗范围为[Max_i-C2,Max_i+C1],,这样有效避免以上问题出现,C1、C2大小根据实际数据调整，经试验表明，C1取值0.01s-0.03s，C2取值0.04s-0.14s效果较好。

(7)利用离散Curvelet变化对时窗有效区域进行多尺度分解，将不同频率的信号剥离开来，进一步提高有效数据信噪比。

(8)采用FAST-AIC法算法在精确时窗内进行精确初至拾取。对一个长度为N的地震记录x(i),(i＝1,2,...,N),AIC检测器定义为:

AIC(k)＝k×lg(Var(x[1,k]))+(N-k-1)lg(Var(x[k+1,N]))

将上式变换为求取离散实数序列的算术和、平方和的线性组合形式,减少重复计算,极大地提高了计算效率。

将上式展开并变换得到：

其中令当微地震数据长度一定时，sm为定值，则有：

综上可得快速AIC算法如下：

其中式中x_i为微地震波形序列，k的取值对应x_i中采样点，N为微地震记录长度。图5所示为改进的FAST-AIC算法的计算效率与传统的AIC算法的对比。

本发明利用震源S₁、S₀、S₂依次产生一组正演模拟数据DATA，取测线L1的12道数据进行可行性和准确性的验证。DATA1是模拟数据DATA加入随机噪声(SNR＝4)后的合成数据。图6为模拟数据DATA1道集示意图。图7为利用射孔信号初至时差文件Δt_i校正模拟数据DATA1，图8(a)为将校正后的模拟数据各道信号横向叠加获得模型道，图8(b)为利用hilbert变换求取模型道能量瞬时包络E(i)，图8(c)为采用STA/LTA法基于模型道能量包络E(i)识别微地震事件，图9为精确时窗内数据放大图，图10精确时窗内数据滤波后结果图以及拾取结果。

本发明采用常规方法和本发明方法分别对四组加入不同强度噪声的数据进行初至拾取，由表1可知，在相同信噪比情况下，新方法在精度上具有明显优势。在计算效率上，信噪比的改变对两种方法耗时影响较小，但本文方法相较于常规方法耗时更短，由于本文采用数据较短，表现不明显，当处理海量微地震监测数据时，计算效率会有显著提高。

表1常规方法与本文方法对不同信噪比数据拾取对比图

为了验证新方法的优越性，利用常规方法和本发明方法对模拟数据进行初至拾取对比。图11为现有STA/LTA法识别微地震事件结果图，由图可知微地震事件识别漏检、误检情况的发生对后续初至拾取结果精度会产生直接影响，拾取精度较低。图12为FAST-AIC法拾取初至结果图，具有更好的拾取精度。

为验证新的微地震信号初至拾取方法的明显优势，本发明对山西西北部某盆地进行实际压裂资料初至拾取，图13为地面实际阵列式观测系统排布图。图14是山西某地微地震信号的完整记录，图15为根据射孔信号初至时差文件Δt_i校正微地震数据结果图，图16(a)是校正后的模拟数据各道信号横向叠加获得模型道，图16(b)所示是利用Hilbert变换求取模型道M(i)的能量包络E(i)，图16(c)为STA/LTA法基于模型道能量包络E(i)识别微地震事件结果，图17为采用本发明提出的方法与STA/LTA法、AR-AIC法、分形维法、以及人工拾取方法进行对比，17(a)为AR-AIC方法拾取结果图，17(b)为STA/LTA方法拾取结果图，17(c)为分形维法拾取结果图，17(d)为本文方法拾取结果图，17(e)为人工拾取结果图。图18绘制了四种方法拾取结果与人工拾取结果的绝对误差与实际数据信噪比(SNR)关系散点图。

Claims (4)

1.一种基于FAST-AIC法的微地震信号初至自动拾取方法，其特征在于，包括以下步骤：

A、输入微地震监测数据以及射孔信号；

B、人工拾取每条测线上检波器记录的射孔信号初至时刻，设置每条测线靠近射孔点的检波器记录信号为参考道，计算射孔信号初至时差文件Δt_i，即使每条测线上检波器记录的射孔信号初至时刻与该条侧线上参考道的初至时刻的差值；

C、利用射孔信号初至时差文件Δt_i校正微地震数据；

D、将校正后的多道微地震数据叠加获得一个模型道M(i)；

E、对模型道M(i)做希尔伯特(Hilbert)变换求取能量瞬时包络E(i)；

$E\left(i\right)=\sqrt{M{\left(i\right)}^2+\stackrel{\‾}{M}{\left(i\right)}^2}---\left(1\right)$

其中为模型道M(i)的Hilbert变换

F、采用STA/LTA法基于模型道能量包络E(i)识别微地震事件，当

STA/LTA比值超过阈值Thr时，认为检测到微地震事件，并且此时对应的时刻T_i记为微地震触发时刻，STA/LTA法定义如下：

$R=\frac{STA}{LTA}=\frac{\underset{i=T_i}{\overset{T_i+m-1}{\Σ}}{E}_i^2/m}{\underset{i=T_i-n-1}{\overset{T_i}{\Σ}}{E}_i^2/n}---\left(2\right)$

其中，E_i代表微地震信号，m表示短时窗长度，n表示长时窗长度,N表示模型道长度,Thr为根据实际数据调节的经验值。

G、采用检波器布阵方式在触发时刻T_i处设计初步时窗。

H、根据初步时窗内各道信号最大值确定精确时窗。

I、利用离散Curvelet变化对时窗有效区域进行多尺度分解，将不同频率的信号剥离开；

J、采用FAST-AIC法算法在精确时窗内进行准确初至拾取。

2.按照权利要求1所述的一种FAST-AIC法微地震信号初至自动拾取方法，其特征在于，G步骤所述的初步时窗设计如下：

首先，分析震源S相对测线L1最近及最远距离位置(震源S到检波器A₁与A_n走时差的最大值与最小值)，即当公式(3)Δe取最大值和最小值时分别对应的震源坐标。

$\mathrm{\Δ}e=\stackrel{\‾}{{\mathrm{SA}}_n}-\stackrel{\‾}{{\mathrm{SA}}_1}---\left(3\right)$

其中

$\stackrel{\‾}{{\mathrm{SA}}_1}=\sqrt{{(x-x_1)}^2+{(y-y_1)}^2+{(z-z_1)}^2}---\left(5\right)$

其中(x₁,y₁,z₁)为检波器A₁坐标,(x_n,y_n,z_n)为A_n坐标,S为震源分布球体中一点，坐标为(x,y,z),分别代表震源点S到检波器A_n、A₁的距离。

由于震源点S在以射孔点(0,0,h)为圆心，半径为R的球形内，可知震源点S满足以下约束条件：

0≤x²+y²+(z-h)²≤R² (6)

根据约束条件公式(6)，联合公式(3)-(5)可求得Δe取最大值对应震源S₂坐标和最小值时分别对应的震源S₁坐标。

其次，根据两个极限震源坐标定义初步时窗为[T_i-Δ_i1,T_i+Δ_i2]。

${\mathrm{\Δ}}_{i1}=\frac{|\stackrel{\‾}{S_2{A}_n}-\stackrel{\‾}{S_2{A}_1}|-|\stackrel{\‾}{S_0{A}_n}-\stackrel{\‾}{S_0{A}_1}|}{v}---\left(7\right)$

${\mathrm{\Δ}}_{i2}=\frac{|\stackrel{\‾}{S_2{A}_n}-\stackrel{\‾}{S_2{A}_1}|}{v}---\left(8\right)$

分别代表他们两点间的距离，v代表地层速度。

3.按照权利要求1所述的一种FAST-AIC法微地震信号初至自动拾取方法，其特征在于，H步骤所述的精确时窗获取；

精确时窗获取包括以下步骤：

a、寻找时窗[T_i-Δ_i1,T_i+Δ_i2]内各道信号最大值Max_i；

b、更改时窗下限为各道信号幅值最大值对应时刻之后C1s处；

c、在各道时窗内最大值对应时刻向前C2s作为时窗起始点，即精确时窗范围为[Max_i-C2,Max_i+C1],其中C1、C2为经验值。

4.按照权利要求1所述的一种FAST-AIC法微地震信号初至自动拾取方法，其特征在于，J步骤所述的FAST-AIC法是将AIC公式(9)变换为求取离散实数序列的算术和、平方和的线性组合形式的FAST-AIC公式(13)所示:

AIC(k)＝k×lg(Var(x[1,k]))+(N-k-1)lg(Var(x[k+1,N])) (9)

将公式(9)变换为求取离散实数序列的算术和、平方和的线性组合形式得到公式(10)和公式(11)。

$\begin{array}{c}\lg \left(Var\right(x\[1,k\]\left)\right)=\lg \left(\frac{1}{k}\underset{i=1}{\overset{k}{\Σ}}{\left(x_i-\stackrel{\‾}{x}\right)}^2\right)\\ =\lg (\frac{1}{k}\underset{i=1}{\overset{k}{\Σ}}{x}_i^2+\frac{1}{k}\underset{i=1}{\overset{k}{\Σ}}{\stackrel{\‾}{x}}^2-\frac{2\stackrel{\‾}{x}}{k}\underset{i=1}{\overset{k}{\Σ}}{x}_i)\\ =\lg (\frac{1}{k}\underset{i=1}{\overset{k}{\Σ}}{x}_i^2-2\stackrel{\‾}{x}\frac{1}{k}\underset{i=1}{\overset{k}{\Σ}}{x}_i+{\stackrel{\‾}{x}}^2)\end{array}---\left(10\right)$

$\begin{array}{c}\lg \left(Var\right(x\[k+1,N\]\left)\right)=\lg \left(\frac{1}{N-k}\underset{i=k+1}{\overset{N}{\Σ}}{\left(x_i-\stackrel{\‾}{x}\right)}^2\right)\\ =\lg (\frac{1}{N-k}\underset{i=k+1}{\overset{N}{\Σ}}{x}_i^2+\frac{1}{N-k}\underset{i=k+1}{\overset{N}{\Σ}}{\stackrel{\‾}{x}}^2-\frac{2\stackrel{\‾}{x}}{N-k}\underset{i=k+1}{\overset{N}{\Σ}}{x}_i)\\ =\lg (\frac{1}{N-k}\underset{i=k+1}{\overset{N}{\Σ}}{x}_i^2-2\stackrel{\‾}{x}\frac{1}{N-k}\underset{i=k+1}{\overset{N}{\Σ}}{x}_i+{\stackrel{\‾}{x}}^2)\end{array}---\left(11\right)$

其中并且设当地震数据长度N固定时,sm即为定值因此:

$\begin{array}{c}\lg \left(Var\right(x\[k+1,N\]\left)\right)=\lg (\frac{1}{N-k}\×(sm-\underset{i=1}{\overset{k}{\Σ}}{x}_i^2)\\ -\frac{2\stackrel{\‾}{x}}{N-k}\×(N\stackrel{\‾}{x}-\underset{i=1}{\overset{k}{\Σ}}{x}_i)+{\stackrel{\‾}{x}}^2)\end{array}---\left(12\right)$

设因此：

$\begin{array}{c}AIC\left(k\right)=k\×\lg (\frac{1}{k}{\mathrm{sk}}_2-2\stackrel{\‾}{x}\frac{1}{k}{\mathrm{sk}}_1+{\stackrel{\‾}{x}}^2)\\ +(N-k-1)\×\lg (\frac{1}{N-k}\×(sm-{\mathrm{sk}}_2)-\frac{2\stackrel{\‾}{x}}{N-k}\×(N\stackrel{\‾}{x}-{\mathrm{sk}}_1)+{\stackrel{\‾}{x}}^2)\end{array}---\left(13\right)$

式中x_i为微地震数据，k的取值对应x_i中采样点，N为微地震记录长度。采用(13)式计算的优点是,当地震波长度N固定时，sm、为定值，只需要计算一次，不随时窗分界点k的移动而改变。当k移动到新的位置时，只需计算新的sk1,sk2,并且在计算sk1,sk2的过程中不必重复计算时窗内所有点的值,只须在上一个时窗值的基础上,加上新进入时窗采样点的值或值的平方即可，大大减少计算量。