CN106646610A - 一种利用偏振约束aic算法自动拾取微震初至的算法 - Google Patents

Abstract

本发明公开了一种利用偏振约束AIC算法自动拾取微震初至的算法，解决了现有技术拾取微震初至精确不够、动化程度低、以及时窗选取多且受时窗选取长短影响大的问题。本发明首先判断微震数据是单分量或三分量微震数据，三分量对应三个特征值序列，单分量对应一个特征值序列；其次利用偏振分析方法计算出偏振协方差矩阵的特征值序列，然后计算出三分量最大特征值序列的最大值点或单分量特征值序列的最大值点，微震事件初至点在上述计算出的最大值点之前；最后计算特征值序列最大值点之前数据的AIC函数，此段AIC函数的最小值点即为微震事件初至点。本发明拾取微震初至精确、动化程度高、以及时窗选取少且几乎不受时窗选取长短的影响。

Description

一种利用偏振约束AIC算法自动拾取微震初至的算法

技术领域

本发明涉及一种利用偏振约束AIC算法自动拾取微震初至的算法。

背景技术

微震监测技术已经被广泛的应用于大坝矿山的安全监测、页岩气开采中水力压裂监测、深埋隧洞稳定性监测等领域，并且取得了许多显著的研究成果，取得这些成果的基础是对微震事件的精确定位，而微震的初至准确拾取是进行微震事件进行定位的必要先决条件之一，因此对微震初至的精确拾取是微震信号处理中的一项必要工作。

微震事件初至的拾取主要是在环境噪声中识别出微震信号的初到时，目前常用的初至拾取算法主要是借鉴天然地震数据处理中的地震波初至拾取算法，如STA/LTA(ShortTime Average/Long Time Average，长短时平均)算法、AIC算法、基于高阶统计量偏斜度和峰度的PAI-S/K算法、分形维数算法、散射积分算法、偏振分析算法等。上述方法虽然在一些工程中成功应用，但相比于传统的地震数据，微震数据中经常会包含机械振动，环境噪声，电磁噪声等多种干扰信号，信噪比普遍偏低，采用单一算法对微震数据进行初至拾取均难以达到理想的效果。

针对微震数据的初至拾取，国内外学者也做了较多的相关研究，如谭玉阳等人针对低信噪比的微震数据的初至拾取，提出了SLPEA算法，其相当于是STA/LTA算法和偏振分析方法的结合，利用STA/LTA算法计算出的能量比函数Ri和偏振法计算出的偏振度函数Pi归一化后相乘得到Ki，通过边缘检测函数Di检测Ki的变化快慢，Ki变化最快的点存在一个局部极大值点，微震的初至点被认为是在Di的极值点附近，再通过一定的约束条件，精准的拾取到微震事件的初至，根据其试验可以看出，该方法具有优良的拾取效果，但是通过对算法的分析可以发现，算法在计算过程中，需要设置6个不同时窗和3个不同阈值，同时计算量偏大，耗时较长，难以满足实时处理的要求。吕世超等人利用偏振约束的STA/LTA方法实现对有效微震事件的自动拾取，并将该方法应用于整套水力压裂微震监测资料中，取得了较好的效果，其方法是利用偏振度函数和STA/LTA算法的能量比函数相乘得到一个整体的判断因子，通过其来识别微震事件，由于其基础算法中包含STA/LTA算法和偏振度函数，算法在实现中仍然至少需要设置3个时窗，同时其对信噪比极低的数据可靠性依然不高。贾瑞生等人通过Huang变换将微震信号分解成从高频到低频的多个基本时间序列的组合，对含噪的微震信号进行降噪处理，对降噪后的微震信号进行Hilbert变换和包络分析，并进行归一化处理，沿时间轴顺序搜索，当微震信号的均一化包络值大于设定阈值时，即可确定微震信号初至的大致时刻，然后计算该时刻前后N个采样点时窗范围内的AIC函数，即可求出微震事件的初至，该方法利用归一化处理后的包络值首次超过阈值的方法对AIC函数进行约束，尽管牺牲了部分计算时间，但是对于低信噪比的微震信号初至拾取具有良好的效果，不过该方法仍需要人工设置阈值和时窗的大小。可以看出以上方法，均能不同程度的提高微震初至拾取的准确度，但是在计算时需要人为的设置时窗或阈值，降低了微震数据处理的自动性，同时精度的提高也牺牲了算法的实时性。

发明内容

本发明要解决的技术问题是：提供利用偏振约束AIC算法自动拾取微震初至的算法，该方法拾取微震初至精确、自动化程度高、能满足实时性要求、以及时窗选取少且不受时窗选取长短的影响。

为实现上述目的，本发明采用的技术方案如下：

一种利用偏振约束AIC算法自动拾取微震初至的算法，包括以下步骤：

(1)判断微震信号数据是单分量微震信号数据或三分量微震信号数据，如果微震信号数据是三分量微震信号数据，则采用步骤(2)估算出微震初至的范围，如果微震信号数据是单分量微震信号数据，则采用步骤(3)估算出微震初至的范围；

(2)利用偏振分析方法计算出微震信号数据的偏振协方差矩阵的特征值序列，三分量微震信号数据中的三个分量分别对应三个特征值序列，运用微震波拾取的特征函数计算三个特征值序列中最大特征值序列的最大值点，所述微震初至在该最大值点之前，计算协方差矩阵的特征值λ₁，λ₂，λ₃，存在λ₁≥λ₂≥λ₃的关系，取其最大特征值λ₁，求λ₁序列最大值对应的时间点p，依此点作为基准点，利用步骤(4)精确计算微震初至；

(3)微震信号数据为单分量微震信号数据时，其可视作为三分量数据特征矩阵仅存在一个分量的特殊情况，以该分量所对应的特征值序列作为其特征值，计算该序列最大值对应的时间点p，利用步骤(4)精确计算出微震初至；

(4)以特征值序列最大值点对应的时间点p为基准时间，采用AIC函数计算出微震数据在时间点p前的AIC函数值，该AIC函数计算出的最小值点即为微震事件的初至点。

进一步地，步骤(2)中所述协方差矩阵的计算公式为：

其中，x_i、y_i和z_i为i时刻对应的三分量微地震数据，cov(x,y)表示数据x和y的协方差，var(x)表示数据x的方差，M为三分量数据的协方差矩阵。

进一步地，步骤(4)中所述AIC函数为：

AIC(n)＝n×lg(var(x[l，n]))+(L-n-1)×lg(var(x[n+1，L])，

其中，x(n)为每个采样点对应的微震数据，n＝1,2,…,L；L＝p+N/10,p为特征值序列最大值对应的点,N为采样点的总长度。

进一步地，步骤(3)中所述单分量所对应的特征值序列为var(x)或var(y)或var(z)。

与现有技术相比，本发明具有以下有益效果：

本发明在微震事件定位分析中对三分量和单分量微震信号数据进行微震初至拾取时均有效，在微震初至拾取时精确度高并且不会受到时窗长短的影响；对于信噪比较低的微震信号数据的初至拾取具有极高的准确性，抗噪性极高，在微震初至拾取过程中仅需设置一个时窗或阀值，极大地提高了微震信号数据处理的自动化程度，很好的满足了微震事件定位分析中的实时性要求。

附图说明

图1为本发明实施流程框图；

图2为本发明实验例16.7dB合成数据采用本发明、STA/LTA算法、Maeda-AIC算法、偏振分析法的拾取效果图；

图3为本发明实验例2.6dB合成数据采用本发明、STA/LTA算法、Maeda-AIC算法、偏振分析法的拾取效果图；

图4为本发明实验例0.2dB合成数据采用本发明、STA/LTA算法、Maeda-AIC算法、偏振分析法的拾取效果图；

图5为本发明实践例微震监测检波器及现场设备布设的示意图；

图6为本发明实践例事件A采用本发明、STA/LTA算法、Maeda-AIC算法、偏振分析法的拾取效果图；

图7为本发明实践例事件B采用本发明、STA/LTA算法、Maeda-AIC算法、偏振分析法的拾取效果图；以及

图8为本发明实践例事件C采用本发明、STA/LTA算法、Maeda-AIC算法、偏振分析法的拾取效果图。

具体实施方式

下面结合附图说明和实施例对本发明作进一步说明，本发明的方式包括但不仅限于以下实施例。

实施例

如图1-8所示，一种利用偏振约束AIC算法自动拾取微震初至的算法，包括以下步骤：

(1)判断微震信号数据是单分量微震信号数据或三分量微震信号数据，如果微震信号数据是三分量微震信号数据则采用步骤(2)估算出微震初至的范围，如果微震信号数据是单分量微震信号数据则采用步骤(3)估算出微震初至的范围；

(2)利用偏振分析方法计算出微震信号数据的偏振协方差矩阵的特征值序列，所述三个分量的微震信号数据分别对应三个特征值序列，运用微震波拾取的特征函数计算最大特征值序列的最大值点，所述微震初至在该最大值点之前，计算协方差矩阵的特征值λ₁，λ₂，λ₃，存在λ₁≥λ₂≥λ₃的关系，取其最大特征值λ₁，求λ₁序列最大值对应的时间点p，依此点作为基准点，利用步骤(4)精确计算微震事件初至，所述协方差矩阵的计算公式为：

其中，x_i、y_i和z_i为i时刻对应的三分量微地震数据，cov(x,y)代表数据x和y的协方差，var(x)代表数据x的方差，M为三分量数据的协方差矩阵；

(3)微震数据为单分量时，其可视作为三分量数据特征矩阵仅存在一个分量var(x)或var(y)或var(z)的特殊情况，以该分量作为其特征值，计算该序列最大值对应的时间点p，利用步骤(4)精确计算出微震事件的初至点；

(4)以特征值序列最大值点对应的时间点p为基准时间，采用下面的公式计算微震数据在时间点p前的AIC函数值，该AIC函数的最小值点即为微震事件的初至点，所述AIC函数为：

AIC(n)＝n×lg(var(x[1，n]))+(L-n-1)×lg(var(x[n+1，L])，

其中，x(n)为每个采样点对应的微震数据，n＝1,2,…,L；L＝p+N/10,p为特征值序列最大值对应的点,N为采样点的总长度。

为了更好地理解本发明，下面采用实验例对本发明做进一步的阐述：

本实例采用6组合成的不同信噪比数据来验证本发明的可靠性，不同信噪比下的稳定性，同时估计本发明计算时间，该6组数据分别用A1、A2、B1、B2、C1、C2表示，其中A1、B1、C1为单分量合成数据，A2、B2、C2为三分量合成数据，6组合成数据均由一个按照指数规律衰减的正弦波叠加不同幅值的随机噪声产生。为了保持和实际微震数据类似，合成数据的采样率也设置为6kHz，数据长度定义为4000个采样点，合成的3对数据信噪比分别为16.7dB、2.6dB和0.2dB。合成地震数据的真实初至在第1000个采样点处。图1为16.7db的单分量数据A1和三分量数据A2的初至拾取情况及其对应的特征值及AIC函数序列，其中，图中竖直虚线为合成数据的实际微震初至；图2为2.6db单分量数据B1和三分量数据B2的初至拾取情况及其对应的特征值及AIC函数序列，其中，图中竖直虚线为合成数据的实际微震初至；图3为0.2dB单分量数据C1和三分量数据C2的初至拾取情况及其对应的特征值及AIC函数序列，其中，图中竖直虚线为合成数据的实际微震初至。表1为计算中每种算法选取的时窗大小，表2为本发明对比STA/LTA算法、Maeda-AIC算法、偏振分析法的拾取结果，表3为本发明对比其他方法的计算时间情况，表4为本发明采用不同长度时窗时对数据B2的拾取精度和计算时间，时窗长度不大于总采样长度的十分之一。

从对合成数据的计算结果中可以看出，由于STA/LTA算法受时窗的大小影响较大，尽管图中选取的为效果较好的时窗结果，但效果依旧不尽人意；由于合成的数据信噪比普遍不高，Maeda-AIC的AIC函数的最小值滞后，导致拾取结果不可靠，偏振分析法具有较高的拾取精度，但是针对低信噪比数据效果不好，本发明相对来说对不同信噪比的数据都有较高的拾取精度，其单分量数据的自动拾取时间比Maeda-AIC方法要短，三分量数据的自动拾取时间也仅比偏振分析法长9％左右。

表1不同算法采用的时窗长度

表2合成数据采用不同方法的拾取结果

表3合成数据采用不同方法的计算时间(单位：秒)

表4不同时窗长度下本发明的拾取精度和计算时间

下面采用实践例对本发明做更进一步地阐述：

本实践例采用的数据为我国西南某水电站尾水隧洞开挖岩爆监测的微震监测数据。仪器采用IMS的20通道微震监测设备，现场布置了10支微震检波器，其中，单分量检波器和三分量检波器各5支，采样频率为6kHz，检波器及现场设备布设的示意图如图4所示。

本实践例选取了A、B、C三次不同信噪比的岩爆事件，分别采用本发明和STA/LTA方法、偏振分析方法以及Maeda-AIC方法进行验证。其中事件A触发了4支检波器，包含2支三分量检波器和2支单分量检波器，事件B触发了4支检波器，包含1支三分量检波器和3支单分量检波器，事件C触发了3支检波器，包含一支三分量传感器和2支单分量检波器；事件A和事件B的微震数据时窗含2499个采样点，时间长度为0.4165s，事件C的数据时窗含2250个采样点，时间长度为0.375s。从波形数据上看，A、B、C三个事件的信噪比依此降低，图5为采用本发明、STA/LTA算法、Maeda-AIC算法、偏振分析法对事件A的拾取结果图，其中，图中a为本发明算法、b为STA/LTA算法、c为Maeda-AIC算法、d为偏振分析法，图中竖直短线为上述相应算法自动拾取微震初至；图6为采用本发明、STA/LTA算法、Maeda-AIC算法、偏振分析法对事件B的拾取效果图，其中，图中a为本发明算法、b为STA/LTA算法、c为Maeda-AIC算法、d为偏振分析法，图中竖直短线为上述相应算法自动拾取微震初至；图7为采用本发明、STA/LTA算法、Maeda-AIC算法、偏振分析法对事件C的拾取结果图，其中，图7中a为本发明算法、b为STA/LTA算法、c为Maeda-AIC算法、d为偏振分析法，图7中竖直短线为上述相应算法自动拾取微震初至。表5为本发明对比STA/LTA算法、Maeda-AIC算法、偏振分析法的计算时间情况。表6(a、b、c)为不同时窗情况下，本发明对A、B、C三个微震事件初至的拾取结果。

表5实际数据采用不同方法的计算时间(单位：秒)

表6(a)不同时窗长度下本发明的初至拾取结果(事件A：数据长度2499)

表6(b)不同时窗长度下本发明的初至拾取结果(事件B：数据长度2499)

表6(c)不同时窗长度下本发明的初至拾取结果(事件C：数据长度2250)

表7 100次实际数据分析中准确率统计

从上述数据的分析中可以看出，微震事件通常会同时触发多支检波器，每个事件中都会存在信噪比较低的数据波形，本发明相比于STA/LTA算法、Maeda-AIC算法、偏振分析法具有较高的抗噪性，同时运算速度也接近Maeda-AIC算法，从表6中可以看出，当时窗长度小于数据长度的十分之一时，本发明的拾取结果几乎不受所选时窗长短的影响，当时窗过长时，对偏振分析的影响会造成拾取结果的错误。表7为100组数据通过本发明拾取的误差统计情况，以人工拾取为参考，从表中可以看出，本发明具有极高的准确性。

本发明通过偏振分析对Maeda-AIC函数的时窗进行约束，缩小了计算Maeda-AIC函数时选取时窗的长度，在一定程度上缩减了程序的计算时间；同时，针对单分量微震数据不存在协方差矩阵的情况，直接求单通道数据的特征值，利用其特征值对Maeda-AIC时窗进行约束。通过与STA/LTA算法、Maeda-AIC算法、偏振分析法的对比，本发明针对低信噪比数据的拾取具有较高的准确率；同时由于本发明不用人工设置阈值，时窗的长短对计算结果不存在影响，只会影响运算的速度，与其他几种改进算法相比，本发明具有更高的自动化程度。

上述实施例仅为本发明的优选实施方式之一，不应当用于限制本发明的保护范围，但凡在本发明的主体设计思想和精神上作出的毫无实质意义的改动或润色，其所解决的技术问题仍然与本发明一致的，均应当包含在本发明的保护范围之内。

Claims (4)

1.一种利用偏振约束AIC算法自动拾取微震初至的算法，其特征在于，包括以下步骤：

(1)判断微震信号数据是单分量微震信号数据或三分量微震信号数据，如果微震信号数据是三分量微震信号数据，则采用步骤(2)估算出微震初至的范围，如果微震信号数据是单分量微震信号数据，则采用步骤(3)估算出微震初至的范围；

(2)利用偏振分析方法计算出微震信号数据的偏振协方差矩阵的特征值序列，三分量微震信号数据中的三个分量分别对应三个特征值序列，运用微震波拾取的特征函数计算三个特征值序列中最大特征值序列的最大值点，所述微震初至在该最大值点之前，计算协方差矩阵的特征值λ₁，λ₂，λ₃，存在λ₁≥λ₂≥λ₃的关系，取其最大特征值λ₁，求λ₁序列最大值对应的时间点p，依此点作为基准点，利用步骤(4)精确计算微震初至；

(3)微震信号数据为单分量微震信号数据时，其可视作为三分量数据特征矩阵仅存在一个分量的特殊情况，以该分量所对应的特征值序列作为其特征值，计算该序列最大值对应的时间点p，利用步骤(4)精确计算出微震初至；

(4)以特征值序列最大值点对应的时间点p为基准时间，采用AIC函数计算出微震数据在时间点p前的AIC函数值，该AIC函数计算出的最小值点即为微震事件的初至点。

2.根据权利要求1所述的一种利用偏振约束AIC算法自动拾取微震初至的算法，其特征在于，步骤(2)中所述协方差矩阵的计算公式为：

$\cos (x,y)=\frac{1}{T}\underset{i=T_1}{\overset{T_2}{\Σ}}(x_i-\stackrel{\‾}{x})(y_i-\stackrel{\‾}{y})$

$var\left(x\right)=\frac{1}{T}\underset{i=T_1}{\overset{T_2}{\Σ}}{(x_i-\stackrel{\‾}{x})}^2$

$M_i=\left[\begin{array}{ccc}\mathrm{var}\left(x\right)& \mathrm{cov}(x,y)& \mathrm{cov}(x,z)\\ \mathrm{cov}(x,y)& \mathrm{var}\left(y\right)& \mathrm{cov}(y,z)\\ \mathrm{cov}(x,z)& \mathrm{cov}(y,z)& \mathrm{var}\left(z\right)\end{array}\right]$

其中，x_i、y_i和z_i为i时刻对应的三分量微地震数据，cov(x，y)表示数据x和y的协方差，var(x)表示数据x的方差，M为三分量数据的协方差矩阵。

3.根据权利要求1或2所述的一种利用偏振约束AIC算法自动拾取微震初至的算法，其特征在于，步骤(4)中所述AIC函数为：

AIC(n)＝n×lg(var(x[1，n]))+(L-n-1)×lg(var(x[n+1，L])，

其中，x(n)为每个采样点对应的微震数据，n＝1，2，…，L；L＝p+N/10，p为特征值序列最大值对应的点，N为采样点的总长度。

4.根据权利要求2所述的一种利用偏振约束AIC算法自动拾取微震初至的算法，其特征在于，步骤(3)中所述单分量所对应的特征值序列为var(x)或var(y)或var(z)。