Experimental study on seismic wave attenuation compensation based on convolution principle and improved generalized S-transform
-
摘要:
为了恢复震动波能量在传播过程中产生的衰减损耗，提出基于褶积原理求取品质因子Q的方法与改进广义S变换相结合的反Q滤波法。通过震动波衰减补偿模型试验，对试验数据进行改进广义S变换的时频特性分析，得出了信号的能量分布情况以及时间频率对应关系；采用基于褶积原理求取品质因子的方法，得到时变Q值；对试验数据进行反Q滤波处理，使震动波能量得到了补偿。结果表明本文提出的反Q滤波法提高了对震动波能量衰减补偿的效果，拓宽了地震资料的频带，提高了地震资料分辨率，有利于进行高分辨率地震勘探、深部信号增强和油气藏预测工作的开展。
Abstract:In order to recover the attenuation loss caused by the vibration wave energy in the process of propagation, the inverse Q filtering method based on the convolution principle and the improved generalized S-transform is proposed. Through the vibration attenuation compensation model test, the time-frequency characteristics of the test data are analyzed by improving the generalized S-transform, and the energy distribution of the signal and the corresponding relationship between time and frequency are obtained. The method of quality factor based on convolution principle is used to obtain the time-varying Q value. The test data are processed using inverse Q filtering, and the vibration wave energy is compensated. The results show that the inverse Q filtering method proposed in this paper improves the compensation effect of seismic wave energy attenuation, broadens the frequency band of seismic data, improves the resolution of seismic data, and is conducive to the development of high-resolution seismic exploration, deep signal enhancement and oil and gas reservoir prediction.
-
引言
微震监测是通过对微小地震事件的实时监测来实现对地压稳定性或地下施工过程的安全监控与管理，在开采、矿藏勘探、大坝监测、矿山安全等工程领域，尤其是井下水力压裂、救援等方面应用广泛。准确可靠的P波初至到时拾取是微震监测系统的关键技术之一，但井下工作环境复杂，实时监测的微震数据中经常混入较多的噪声干扰，加之微震信号本身能量微弱，频率范围广，且持续时间短，噪声干扰的增加无疑给P波到时拾取带来了巨大挑战，因此，快速、准确、可靠地拾取微震P波到时成为业界研究的重点。
微震P波到时拾取算法，目前应用较广泛的有长短时窗比值（short term averaging/long term averaging，缩写为STA/LTA）（刘晗，张建中，2014；刘晓明等，2017）、赤池信息准则（Akaike Information Criterion，缩写为AIC）（赵大鹏等，2013；朱权洁等，2018）。STA/LTA快速简单，但在信噪比（signal-to-noise ratio，缩写为SNR）较低或微震不明显时准确度低；AIC方法较稳定，但在低SNR时准确度有所下降。针对STA/LTA算法的不足，刘晗和张建中（2014）引入长短时窗内信号幅度的标准差之比的加权系数来放大信号，增加检测灵敏度，但该系数对噪声亦有放大作用，使检测准确度降低；刘晓明等（2017）基于信号一阶导数引入权重因子K构建新特征函数，从而快速反映信号振幅及频率的变化，提高检测灵敏度，但该因子为全局性值，未根据信号幅度变化作灵活调整，存在一定的缺陷。
单一方法进行初至拾取各有适用条件和局限性，近年，利用各种算法优势的综合分析方法备受关注，尤其是将地震信号先进行变换或分解，再用AIC准则或高阶统计量提取P波到时，能大大提高识别准确度。例如：希尔伯特-黄变换与AIC结合（贾瑞生等，2015）；小波分解与AIC结合（Gaci，2014）；小波包分解AIC的P波初至拾取（田优平，赵爱华，2016）；基于离散小波变换（discrete wavelet transform，缩写为DWT）与峰度进行P波到时精确估计（Li et al，2016）；利用多尺度离散平稳小波变换（discrete stationary wavelet transform，缩写为DSWT）和峰度/峭度实现P波自动拾取（Charles，Ge，2018）等。
经验模态分解（empirical mode decomposition，缩写为EMD）是近年发展起来的自适应信号分解方法，非常适合处理非线性非平稳的地震信号，其被运用于地震波初至拾取也处于不断探索和研究中。例如：Liu等（2017）利用EMD方法分解信号，找到本征模函数（intrinsic mode functions，缩写为IMF）分量1（IMF1）的瞬时频率突增位置，再对该段记录利用AIC准则来确认信号初至时间；Li等（2017）对地震信号EMD后，以IMF5的P波到时为原信号的P波到时；Shang等（2018）对地震信号进行EMD，对IMF1和IMF2去噪后重构原信号，再利用余弦函数去除工频噪声，最后基于AIC计算P波到时；Kirbas和Peker （2017）对地震信号进行EMD，选择IMF3计算滑动时间窗口长度的Teager能量算子（Teager-Kalser energy operator，缩写为TKEO），以TKEO大于设定阈值且为最大值的位置为P波到时；李伟等（2018）利用集合经验模态分解（ensemble empirical mode decomposition，缩写为EEMD）进行信号分解后再用Hankel矩阵奇异值分解（singular value decomposition with Hankle matrix，缩写为Hankle SVD）法降噪，最后利用AIC进行微震初至拾取。这些算法通过EMD分解原记录，利用某些IMF的细节特征进行到时拾取，取得了一定的效果，但EMD算法本身存在模态混叠及伪分量的问题，这无疑会影响拾取的准确度。
基于以上背景，为解决低信噪比下初至拾取准确率低的问题，本文拟利用STA/LTA快速和AIC识别稳定的特点，同时结合变分模态分解（variational mode decomposition，缩写为VMD）可有效减少伪分量及抑制模态混叠的优势，提出一种基于改进STA/LTA和自适应VMD-峰度AIC的两步识别P波初至的算法。首先，改进传统STA/LTA算法，引入反映信号幅度实时变化的权重因子，构建能反映信号能量和频率变化的特征函数，提高P波初次识别的准确度；其次，对初次判别的P波到时前后2—3 s信号进行优化VMD，使IMF分量既能准确地描述原信号的细节特征，又可防止出现过分解和模态混叠；最后，对分解后的IMF进行峰度AIC拾取，并根据各IMF的能量比进行综合加权，实现P波到时的精确拾取。
1. 基本算法原理
1.1 STA/LTA算法原理
STA/LTA算法根据噪声与有效地震信号的振幅、能量或频率等特征的变化趋势来拾取P波到时，具体方法为分别计算长、短时窗内定义的特征函数值的STA和LTA，其中STA刻画的是地震信号的变化趋势，而LTA刻画的是背景噪声的变化趋势。当地震事件到来时，引起STA的变化快于LTA，即STA/LTA会出现突增。当STA/LTA大于预先设定的阈值时，则判定为有效地震P波到达。特征函数是基于原始地震记录构建的能灵敏反映信号幅度或频率变换的新序列，目前广泛使用的有
$$\left\{ \begin{split}& CF （ i ） ＝ y{ （ i ） ^2} ， \\& CF （ i ） ＝ y{ （ i ） ^2} － y （ {i － 1} ） y （ {i ＋ 1} ） ， \\& CF （ i ） ＝ y{ （ i ） ^2} ＋ {[y （ i ） － y （ {i － 1} ） ]^2} \;， \end{split}\right. $$ (1) STA/LTA的计算方法为
$$ \left\{\begin{split}& {\rm{ST}}{{\rm{A}}}_{i}＝{\rm{ST}}{{\rm{A}}}_{i－1}＋\frac{CF （ i ） －{\rm{ST}}{{\rm{A}}}_{i－1}}{{L}_{{\rm{STA}}}}，\\& {\rm{LT}}{{\rm{A}}}_{i}＝{\rm{LT}}{{\rm{A}}}_{i－1}＋\frac{CF （ i－{L}_{{\rm{STA}}}－1 ） －{\rm{LT}}{{\rm{A}}}_{i－1}}{{L}_{{\rm{LTA}}}}，\\& \frac{{\rm{STA}}}{{\rm{LTA}}}＝\frac{{\rm{ST}}{{\rm{A}}}_{i}}{{\rm{LT}}{{\rm{A}}}_{i}}，\end{split} \right. $$ (2) 式中，CF（i）为i时刻地震信号的特征函数值，LSTA为短时窗长度，LLTA为长时窗长度。STA/LTA拾取结果的准确性受特征函数、时窗长度及阈值的影响，需要构建灵敏的特征函数，选取合理的阈值和时窗长度。
1.2 AIC算法基本原理
AIC算法假定地震信号可分为若干局部平稳部分，每部分均可模拟为自回归过程。将地震到来前后的时间序列假设为两个不同的平稳序列，则其AIC值的最小值点为两个序列的分界点，以此确定P波的初达时间。对于长度为
$ N $ 的序列，其AIC值可表示为可表示为$$ {\rm{AIC}} ＝ k\lg ［ {{\rm{var}}} （ x ［ 1， k ］ ） ］ ＋ （ N － k － 1 ） \lg ［ {{\rm{var}}} （ x ［ k ＋ 1， {{N}} ］ ） ］ ， $$ (3) 式中，k为数据窗口内的采样点数，var为序列方差，则AIC最小值对应P波初至到时。根据赵大鹏等（2013）的方法，用峰度（Kurtosis）作为特征函数替代式（3）中的方差能更好地识别微震时间，即Kur-AIC为
$$ \qquad\quad{\rm{Kur}}-{\rm{AIC}} ＝ k\lg \Bigg(\frac{1}{k}\mathop \sum \limits_{j ＝ 1}^k CF （ j ） ^2\Bigg) ＋ （ N － k － 1 ） \lg \Bigg(\frac{1}{{N － k － 1}}\mathop \sum \limits_{j ＝ k}^N CF （ j ） ^2\Bigg)\;{\text{.}} $$ (4) 峰度AIC方法不依赖于时窗长度，计算稳定，对微弱P波亦有较好的分辨能力。AIC算法虽稳定，但当信号信噪比较低时其准确度会大大降低，因为即使是噪声记录的AIC值也会出现一个局部最小值，故AIC不能准确判断一段记录中是否存在有效微震信号，只能用于拾取微震P波的初至时间。因此，一般先利用STA/LTA初步确定微震P波到时区间，再利用AIC进行拾取，降低误差。
1.3 变分模态分解的原理
变分模态分解（VMD）是一种自适应的信号分解方法（Dragomiretskiy，Zosso，2014），能将信号分解为若干从高频到低频，带限、准正交的IMF分量之和，且每个IMF都是一个窄带的调幅调频信号。VMD算法包括构造约束变分模型和变分模型的求解，其实质是在“各IMF之和为原信号，且IMF的估计带宽之和最小”的约束条件下搜寻L个IMF分量。不同于EMD分解分量个数不可人为设定，VMD的IMF分量个数可设定。通过迭代搜寻变分模型的最优解，VMD确定L个IMF的中心频率和带宽，从而实现信号从低频到高频的分解。相比EMD，VMD模态稳定性好，收敛速度快，具有更强的抗噪声能力和局部分解能力，且能有效克服EMD模态混叠及存在伪分量的不足，比较适合进行地震资料处理（Xue et al，2016；Li et al，2018）。
VMD能准确地分解出每一个IMF分量，分解后的IMF瞬时频谱和幅值具有较高的分辨率，不仅能刻画原信号的时频特征，也能更精细地分析信号，更好地描述信号的局部突变或渐变特征（郑小霞等，2017）。因此可先基于VMD对地震信号进行分解，再利用AIC分析各IMF的初至时间。
VMD分解有两个重要参数需要预先设定，一个是分解层数L，一个是用于确保准确重构信号的惩罚因子α。随着分解层数L的不同，VMD会出现不同分解结果，L过大会使信号断断续续无规律，过小则易导致模态混叠，直接影响分解的准确性，因此需要选择合理的分解层数。本文拟对VMD进行优化，根据分解后IMF的特点以及IMF与原信号的相关性自适应确定分解层数，使分解后的IMF不混叠、较纯净，能精细刻画原信号的时频特性，从而进行初至识别。
2. 改进STA/LTA及优化VMD的峰度AIC微震初至拾取
2.1 算法基本思想
针对信噪比低导致P波拾取精确度低的问题，本文对传统STA/LTA算法进行改进，基于地震信号幅度的实时变化引入权重因子，可快速、灵敏地从微震监测数据中提取微地震事件，实现P波的粗略提取；然后选取初次拾取时间点前后2—3 s的记录，根据互相关系数和排列熵准则确定VMD分解层数，使分解后的IMF无伪分量和模态混叠，能较好地表征微震信号特征；最后计算各IMF分量的能量比，并进行峰度AIC识别，以各IMF分量的拾取结果及能量比进行综合加权，得到精确的P波到时，算法流程图如图1所示。
![]()
图 1 基于改进STA/LTA及优化VMD的微震初至拾取算法流程图Figure 1. The flow chart of the first break picking of micro-seismic based on imp-roved STA/LTA and adaptive VMD2.2 基于改进STA/LTA的初步拾取
特征函数的选取决定了初至拾取的精度，因此要构建尽可能灵敏反映信号频率或振幅变化的特征函数，使微震到来时，根据特征函数计算所得的STA发生明显变化，STA/LTA出现较大幅度的增大。本文根据信号幅度的实时变化引入加权系数K，其计算方法为
$$ K ＝\sqrt {\left|\frac{{y （ i ） － y （ i － 1 ） }}{{y （ i － 1 ） }}\right|} ， $$ (5) 式中y为地震信号序列。当微震事件未到来时，信号幅度变化较小，K较小（接近于0），而当微震到来，信号幅度发生较大变化时，K将增大，可见K能实时反映信号的幅度变化。利用K构造新的特征函数，即
$$ CF （ i ） ＝ y{ （ i ） ^2} ＋ K{ ［ y （ i ） － y （ {i － 1} ） ］ ^2} {\text{.}} $$ (6) y（i）2反映信号幅度特性，［ y（i）－y（i－1） ］ 2反映信号频率特性。K能根据信号采样频率对信号幅度和频率进行加权分配，K的引入将增加特征函数的敏感度。当微震事件到来时，K较大，对特征函数有放大作用，使特征函数随之发生较大变化，STA/LTA可迅速反映出信号幅度或频率的变化，从而更易识别微震事件。相比刘晓明等（2017）的研究中权重因子不能根据信号幅度变化灵活调整的不足，本文特征函数能更好地反映并放大地震信号的实时变化，从而提高微震信号检测的灵敏度。
改进的STA/LTA可快速可靠地初步识别微震，为了提高识别精度，取初次拾取时间点前后2—3 s的记录进行二次拾取。由于VMD算法能较好地描述信号的局部突变或渐变特征，本文对VMD算法进行优化，对分解后的IMF进行基于峰度-AIC的二次拾取。
2.3 基于互相关系数和排列熵优化的VMD
经VMD后，地震信号s（t）可表示为
$$ s （ t ） ＝ \mathop \sum \limits_{i ＝ 1}^L {s_{{\rm{IM}}{{\rm{F}}_i}}} （ t ） ＋ {r_n} （ t ） ， $$ (7) 式中，sIMF（t）为从高频到低频的IMF分量，rn（t）为残差分量。为确保分解后IMF既不包含伪分量，又能较好地表征原信号的特征，本文根据互相关系数和排列熵准则来确定VMD分解层数。互相关系数的计算方法为
$$ r （ x， y ） ＝ \frac{{{{\rm{cov}}} （ x， y ） }}{{\sqrt {{\rm{var}} （ x ） {\rm{var}} （ y ） } }} ＝ \frac{\displaystyle\sum \limits_{i ＝ 1}^n （ {{x_i} － \bar x} ） （ {{y_i} － \bar y} ） }{\sqrt {\displaystyle\sum \limits_{i ＝ 1}^n {{ （ {{x_i} － \bar x} ） }^2}} \sqrt {\displaystyle\sum \limits_{i ＝ 1}^n {{ （ {{y_i} － \bar y} ） }^2}} } \;\; {\text{.}}$$ (8) 互相关系数衡量的是两个信号的相关程度，根据互相关系数理论，当两个信号间的互相关系数大于0.3时，表明相关性较大（郑近德等，2013），计算分解后IMF与原信号的互相关系数，若互相关系数不小于0.3，则表明不含伪分量，反之则表明出现了过分解。
为了进一步确定分解层数，本文利用排列熵（Bandt，Pompe，2002）进行分解后IMF的随机性和异常检测。排列熵能较好地反映时间序列细节的微小变化，有效进行异常或突变检测。其计算方法为：
① 设原时间序列为x（i）（i＝1，2，···，N），利用相空间重构得到原序列的重构矩阵，即
$$ \left\{\begin{split}& x （ 1 ） ＝ \{ x （ 1 ） ， x （ 1 ＋ \lambda ） ， \cdots， x （ 1 ＋ （ m － 1 ） \lambda ） \} ， \\& \qquad \qquad \qquad \qquad \qquad\vdots \\& x （ i ） ＝ \{ x （ i ） ， x （ i ＋ \lambda ） ， \cdots， x （ i ＋ （ m － 1 ） \lambda ） \} ， \\& x （ N － （ m － 1 ） \lambda ） ＝ \{ x （ N － （ m － 1 ） \lambda ） ， \cdots ， x （ N ） \} ； \end{split} \right.$$ (9) ② 对
$x （ i ） $ 作升序排列，得到序列$ x （ i ） ＝ \{ x （ i ＋ （ {j_1} － 1 ） \lambda \} ， \{ x （ i ＋ （ {j_2} － 1 ） \lambda \} ， \cdots ， \{ x （ i ＋ （ {j_m} － 1 ） \lambda \} $ ，其中，$ \{ x （ i ＋ （ {j_1} － 1 ） \lambda \} {\text{≤}} \{ x （ i ＋ （ {j_2} － 1 ） \lambda \} {\text{≤}} \cdots {\text{≤}} \{ x （ i ＋ （ {j_m} － 1 ） \lambda \} $ ，当序列值相等时，按j的大小排序；③ 对于排序后的序列x（i），得到一个符号序列sg＝{j1，j2，···，jm}，g＝1，2，···，k，k≤m!，即{j1，j2，···，jm}有m!种不同的排列，而sg是这m!种的其中一种，之后计算每一种符号序列出现的概率
$ {P_g} （ {g = 1， 2， \cdots ， k} ） ， \displaystyle\sum \limits_{g ＝ 1}^k {P_g} = 1 $ ，则x（i）的排列熵${H_p} （ m ） ＝ - \displaystyle\sum\limits_{g ＝ 1}^k {{P_g}{\rm{ln}}{P_g}}$ ，当${P_g} ＝ 1/{{m!}}$ 时，${H_p} （ m ） ＝ \ln （ m! ） $ 最大，最后对Hp（m）进行归一化$\;{H_p} ＝ {H_p} （ m ） /\ln （ m! ） $ 。排列熵通过时间序列的相邻数据进行对比而不是考虑数据的具体值来度量序列随机性，其变化反映并放大了原序列x（i）的微弱变化，从而易于检测信号异常。对归一化后的排列熵，其值越趋近于1，表明信号越随机，据此可判断分解后IMF是否为异常信号。根据王志坚等（2018）的研究，分解到某一层后的IMF分量排列熵不大于0.6，则表明分解后的IMF无异常分量，能较好地刻画原信号的特征，否则表明出现过分解。基于互相关系数和排列熵准则优化VMD，确定分解层数L的算法步骤为：① 设定初始分解层数i＝2；② 基于VMD对原信号进行分解，得到一系列IMFj（j＝1，2，···，i），计算IMFj与原信号的互相关系数corj，以及IMFj的排列熵Hp；③ 若corj≤0.3，且Hp≥0.6，则表明出现过分解，存在异常分量，则停止分解，分解层数L＝i－1，否则令i＝i＋1进入第2步继续VMD分解。经优化VMD分解后的IMF既不包含异常分量，又能刻画原信号的起伏特征和细节信息。
2.4 基于优化VMD的峰度AIC精确识别P波初至
为了提高拾取精度，对改进STA/LTA拾取的P波到时tP前后2—3 s记录进行二次拾取，以获得精确的P波到时。根据互相关系数和排列熵确定分解层数，对这段记录进行优化VMD分解，实现信号分解。如前所示，VMD中参数惩罚因子α需提前设定，一般为采样频率的0.25—2倍。α越小，分解后IMF分量带宽越大，反之分解后IMF分量带宽越小 （唐贵基，王晓龙，2015）。鉴于峰度AIC的稳定性和精度高，对分解出的L个IMF根据式（4）计算峰度AIC值，以最小值为IMFi（i≤L）的拾取时间tAi。VMD分解后IMF分量的中心频率由低至高分布，且具有不同的能量，其中低频IMF是包含有效微震信号自身特征信息的主模态分量，占原信号总能量中的主要部分。因此，对L个IMF计算能量比，其定义为
$$ ER_i ＝ \frac{{|{\rm{IM}}{{\rm{F}}_i}{|^2}}}{{\displaystyle\sum\limits_{i ＝ 1}^L {|{\rm{IM}}{{\rm{F}}_i}{|^2}} }} . $$ (10) 根据各IMF的拾取时间及计算所得能量比，按式（11）进行加权，得到最终的P波拾取到时tAi，即
$$ t_{{\rm{A}}i} ＝ \displaystyle\sum\limits_{i ＝ 1}^L {E{R_i}} t_{{\rm{A}}i}{\text{.}} $$ (11) 3. 实验结果及分析
3.1 改进STA/LTA的模拟微震拾取
设雷克（Ricker）子波主频为20 Hz，采样频率为1 000 Hz，基于卷积模型得到的理论地震记录如图2所示，其中加入不同SNR的随机噪声后可得模拟地震记录，当SNR为−5 dB时，由于噪声严重，起震不明显，且有效微震信号被噪声严重干扰，增加拾取难度。
![]()
图 2 基于雷克子波和卷积模型的理论和模拟地震记录Figure 2. Theoretical and simulated seismogram based on Ricker wavelet and convolution model设STA＝0.1 s，LTA＝0.5 s，微震阈值为1.5，分别采用本文改进STA/LTA，传统STA/LTA，刘晓明等（2017）的改进STA/LTA进行初次识别，结果如图3所示，人工拾取时间为3.065 s，本文改进STA/LTA识别时间为3.125 s，传统STA/LTA识别时间为3.218 s，刘晓明等的改进STA/LTA拾取时间为3.203 s。可见，STA/LTA算法拾取时间有延迟，相比传统STA/LTA，刘晓明等的改进STA/LTA和本文的改进STA/LTA算法能更快地反映信号变化，从而快速检测微震事件。相比刘晓明等的改进SAT/LTA，本文改进STA/LTA权重因子K能实时反映信号幅度变化，而非固定值，特征函数能更快速、真实地反映信号变化，可快速增加超过阈值，从而能快速识别微震。以人工拾取结果为参考，本文改进STA/LTA拾取误差为0.06 s，刘晓明等的改进STA/LTA和传统STA/LTA的误差分别为0.138 s和0.153 s，可见本文改进STA/LTA算法能减小初次拾取误差。
为了验证改进STA/LTA的性能及抗噪性，将此理论地震记录随机添加不同SNR噪声进行拾取（SNR为−5—20 dB），并随机仿真100次，计算与人工拾取误差的均值。图4给出了3种STA/LTA算法在不同SNR下的拾取误差曲线，可见，SNR越低，拾取误差越大。但在较低SNR时，本文改进STA/LTA拾取误差更小，在5—−5 dB时拾取结果与人工拾取的误差为0.05—0.28 s，比传统STA/LTA和刘晓明等的改进STA/LTA算法的拾取误小差约0.02—0.1 s；当SNR在5—−20 dB时，本文与人工拾取的误差为0.02—0.05 s，比传统STA/LTA和刘晓明等的改进STA/LTA算法的拾取误差小约0.01—0.02 s，以上说明本文的改进STA/LTA有效，比传统STA/LTA方法降低误差0.01 s以上，适用于强噪声微震初至拾取，且性能较好，尤其是在低信噪比时，能有效降低初次拾取误差。
![]()
图 4 不同SNR下三种STA/LTA算法拾取误差曲线Figure 4. Picking up error curves of three STA/LTA algorithms under different SNRs3.2 优化VMD峰度AIC的模拟微震拾取
为了提高拾取精度，本文利用优化VMD对信号进行分解，并利用峰度AIC再次识别。根据优化VMD确定分解层数，并将初次拾取时间向前后各推2—3 s，精确确定该记录P波到时。图2中的模拟地震记录，经优化VMD确定的分解层数为2，得到IMF1和IMF2。由图5可见，分解后的IMF不包含噪声异常分量，一定程度上消除噪声影响，且不同频率段的各IMF能较准确地反映原信号中不同成分的振荡信息，保持信号细节特征。根据式（4）计算对应的峰度AIC曲线，以曲线最小值点为初至时间，图5给出了峰度AIC曲线拾取时间分别为3.06 s和3.087 s的IMF1和IMF2，按式（10）计算出的IMF1和IMF2的能量比分别为0.55，0.45，最终加权拾取时间为3.07 s，与人工拾取时间3.065 s的误差为0.005 s，可见优化VMD峰度AIC能有效提高拾取精度。
![]()
图 5 优化VMD峰度AIC为3.07 s时的初至拾取结果Figure 5. First arrival picking based on improved VMD when Kurtosis-AIC is 3.07 s为了验证算法性能，将小波包峰度AIC （田优平，赵爱华，2016）、EMD-AIC算法（Liu et al，2017）与本文优化VMD峰度AIC算法进行对比。如图6，7所示，对于图5的原信号地震记录，基于EMD-AIC的初至拾取时间为3.077 s （3个IMF的拾取时间分别为3.098 s，3.052 s，3.081 s）；小波包峰度AIC的初至拾取时间为3.075 s （三个分量的拾取时间分别为3.073 s，3.078 s，3.075 s），与参考结果的误差分别为0.012 s和0.01 s。可见，EMD-AIC算法、小波包峰度AIC和优化VMD峰度AIC算法均能较准确地拾取到时，且本文算法拾取误差更小，精度较高。
![]()
图 7 基于小波包峰度AIC的初至拾取Figure 7. First arrival picking based on wavelet packet Kurtosis-AIC为了分析本文优化VMD峰度AIC算法在降低二次拾取误差方面的性能，取不同SNR （−5—20 dB）的模拟地震信号，同一SNR添加随机噪声仿真100次，并对比初次与二次拾取时间与人工结果的误差均值。如图8所示，当SNR较低（−5—3 dB）时，初次拾取误差偏大，约为0.05—0.28 s；经二次拾取后，误差明显降低，最终结果与人工拾取结果平均相差在0.023 s以内；随着SNR的提高，初次拾取误差逐渐减小且再次拾取误差大大降低，与人工拾取结果非常接近，误差在0.01 s以内。
为了验证VMD峰度AIC算法在降低二次拾取误差方面的性能，添加不同SNR的噪声并随机仿真100次取均值，对比小波包峰度AIC和EMD-AIC算法与人工拾取结果的误差。如图9所示，三种算法都能较好地拾取到初至P波，本文算法在低SNR下表现稳定，平均误差在0.023 s以内。相比小波包分解和EMD，经优化VMD后IMF分量能更好地保留原信号的细节和突变特征，且有效消除噪声影响，经峰度AIC拾取后能有效降低拾取误差。
3.3 实际微震记录拾取
为测试本文算法拾取真实微震记录的准确性，对2013—2019年间广西平果、武鸣、贵港、乐业等地的若干次ML1.3—2.0的矿藏勘探、开采等人工爆破微震记录共计1 137条进行拾取，图10给出了平果地区2013年7月4日15时3分在（23.343°N，107.583°E）发生的ML＝1.3爆破的某条监测记录，其采样率为100 Hz，爆破持续时间较短，且有较多外部环境噪声，SNR较低，该记录人工拾取的参考时间为14.76 s。本文改进STA/LTA初次拾取时间为14.78 s，传统STA/LTA拾取时间为14.81 s，刘晓明等的改进STA/LTA拾取时间为14.79 s （图11）。可见，同样的阈值，本文改进STA/LTA更敏感快速反映信号变化，减少初次拾取误差。
![]()
图 10 广西平果2013年7月4日ML1.3爆破的一条实际微震记录Figure 10. Real micro-seismic record of some coal mine with ML1.3 at Pingguo，Guangxi on July 4th 2013![]()
图 11 本文改进STA/LTA初次拾取结果Figure 11. First picking results of improved STA/LTA in this paper在二次拾取中，本文优化VMD自适应确定的分解层数为3 （图12），各分量能较好地刻画原地震信号的细节特征，各分量拾取时间分别为14.8 s，14.73 s和14.77 s，能量比为0.31，059，0.1，则最终拾取时间为14.76 s，与人工拾取结果一致。
![]()
图 12 本文改进STA/LTA及优化VMD峰度AIC拾取结果Figure 12. First arrival picking by improved STA/LTA and adaptive VMD Kurtosis-AICEMD-AIC 3个IMF分量的拾取时间分别如图13所示，为14.76 s，14.76 s，14.83 s，最终拾取时间为14.78 s，误差为0.02 s。小波包峰度AIC算法分解的3个IMF拾取时间分别为图14所示的14.78 s，14.76 s，14.76 s，最终拾取时间为14.77 s，误差为0.01 s。3种算法都能较准确拾取初至波，误差在0.02 s内，但本算法误差更小。
对其余记录分别用本算法、EMD-AIC，小波包峰度AIC算法进行拾取，以人工拾取结果为参考，3种算法的拾取结果列于表1。
表 1 实际微震记录初至拾取结果Table 1. First arrival picking of real micro-seismics拾取方法 误差 30 ms 误差 20 ms 误差 10 ms 记录条数 占比 记录条数 占比 记录条数 占比 本文 1 117 98.2% 1 090 95.9% 1 031 90.7% EMD-AIC 1 090 95.9% 1 052 92.5% 980 86.2% 小波包峰度AIC 1 075 94.5% 1 037 91.2% 968 85.1% 本文改进STA/LTA及优化VMD峰度AIC算法拾取误差在30，20，10 ms的比例分别为98.2%，95.9%，90.7%，将误差在0.02 s （20 ms）以内的视为准确拾取，结果显示，本文算法有效，拾取准确率为95.9%，分别高于EMD-AIC和小波包峰度AIC约3.4%和4.7%，且本文拾取误差在10 ms的比例比其它两种算法的要高4.5%以上，表明本文算法拾取误差更小，具有更高的精度。
4. 讨论与结论
本文对传统STA/LTA的特征函数进行了改进，提高P波的初次拾取精度；再对VMD进行了优化，能自适应确定分解层数，对分解后IMF基于峰度AIC二次拾取后，基于各IMF的能量比综合加权得到最终二次拾取到时。实验表明：
1） 改进STA/LTA利用信号幅度设计权重因子构建特征函数，能快速反映信号变化，提高初至P波的初次检测灵敏度；
2） 优化多尺度VMD分解，能根据分解IMF分量的互相关系数和排列熵确定VMD分解层数，既能解决模态混叠，又能防止出现过分解，且充分保留微震信号的细节特征；
3） 仿真实验表明，在较低信噪比时,算法仍能较好识别初至P波，比传统STA/LTA和刘晓明等（2017）改进STA/LTA减小初次识别误差约0.01 s以上；基于优化VMD分解后，能再次降低拾取误差，相比小波包峰度AIC和EMD-AIC拾取算法，本文最终误差基本在0.023 s以内，表明本文算法的有效性、准确性和抗噪性；
4） 实际微震拾取表明，本文拾取结果与人工拾取误差在20 ms内的记录占比95.9%，准确率高于EMD-AIC和小波包峰度AIC，表现出较强性能。
-
![]()
图 1 室内试验模型
（a） 模型搭建；（b） 模型A面
Figure 1. The model of indoor experiment
（a） Model building；（b）A-side of the model
![]()
图 3 1号（a）和4号（b）传感器接收的原始信号
Figure 3. Original microseism signals of sensor 1 （a） and 4 （b）
![]()
图 4 1号（a）和4号（b）传感器接收到的原始信号经过改进广义S变换后的时频谱
Figure 4. Time-frequency spectrum of the original signals received by sensor 1 （a） and 4 （b） after the improved generalized S transform
![]()
图 5 1号（a）和4号（b）传感器的重构信号
Figure 5. Reconstructed signal of sensor 1 （a） and sensor 4 （b）
![]()
图 6 1号（a）和4号（b）传感器接收到原始信号补偿前后的幅值谱
Figure 6. Amplitude spectrum of original signal received by sensor 1 （a） and sensor 4 （b） before and after compensation
![]()
图 7 1号（a）和4号（b）传感器接收的原始信号补偿后时频谱
Figure 7. Time-frequency spectrum of original signal received by the sensor 1 （a） and sensor 4 （b） after compensation
表 1 相似材料用量表
Table 1 The list of similar material consumption
岩层名称 模型厚度/cm 各相似材料用量/kg 沙子 石灰 石膏 水 泥岩层 35 963.2 160.6 160.6 183.4 煤层 60 1694.0 101.7 237.3 290.4 细砂岩层 35 1133.3 149.9 133.6 202.4 
下载: 导出CSV
-
白桦,李鲲鹏. 1999. 基于时频分析的地层吸收补偿[J]. 石油地球物理勘探,34(6):642–648. Bai H,Li K P. 1999. Stratigraphic absorption compensation based on time-frequency analysis[J]. Oil Geophysical Prospecting,34(6):642–648 (in Chinese).
丁进杰,戴永寿,张亚南,陈健,魏玉琴,张漫漫. 2013. 基于高频补偿方法提高地震资料分辨率的初步研究[J]. 地球物理学进展,28(6):3214–3221. Ding J J,Dai Y S,Zhang Y N,Chen J,Wei Y Q,Zhang M M. 2013. Improving seismic data resolution based on high frequency compensation[J]. Progress in Geophysics,28(6):3214–3221 (in Chinese).
高军,凌云,周兴元,牟永光. 1996. 时频域球面发散和吸收补偿[J]. 石油地球物理勘探,31(6):856–866. Gao J,Ling Y,Zhou X Y,Mou Y G. 1996. Compensation for spherical divergence and absorption in time-frequency domain[J]. Oil Geophysical Prospecting,31(6):856–866 (in Chinese).
高静怀,汪文秉,朱光明,彭玉华,王玉贵. 1996. 地震资料处理中小波函数的选取研究[J]. 地球物理学报,39(3):392–400. Gao J H,Wang W B,Zhu G M,Peng Y H,Wang Y G. 1996. On the choice of wavelet functions for seismic data processing[J]. Acta Geophysica Sinica,39(3):392–400 (in Chinese).
高静怀,陈文超,李幼铭,田芳. 2003. 广义S变换与薄互层地震响应分析[J]. 地球物理学报,46(4):526–532. Gao J H,Chen W C,Li Y M,Tian F. 2003. Generalized S transform and seismic response analysis of thin interbeds[J]. Chinese Journal of Geophysics,46(4):526–532 (in Chinese).
何书梅,吴吉忠,魏朋朋,杨倩倩,吴吉厚. 2020. 一种井震融合求取反Q滤波Q体的方法[J]. 物探化探计算技术,42(4):442–450. doi: 10.3969/j.issn.1001-1749.2020.04.02 He S M,Wu J Z,Wei P P,Yang Q Q,Wu J H. 2020. A method for the Q estimation of inverse Q filtering based on well seismic fusion[J]. Computing Techniques for Geophysical and Geochemical Exploration,42(4):442–450 (in Chinese).
贾晓东,翟丽娜,陈石. 2021. 辽宁地区时变重力场源变化特征分析[J]. 地震,41(1):180–190. Jia X D,Zhai L N,Chen S. 2021. Analysis on the characteristics of time-varying gravity field source changes in Liaoning area[J]. Earthquake,41(1):180–190 (in Chinese).
刘喜武,年静波,刘洪. 2006. 基于广义S变换的吸收衰减补偿方法[J]. 石油物探,45(1):9–14. doi: 10.3969/j.issn.1000-1441.2006.01.002 Liu X W,Nian J B,Liu H. 2006. Generalized S-transform based compensation for stratigraphic absorption of seismic attenuation[J]. Geophysical Prospecting for Petroleum,45(1):9–14 (in Chinese).
裴江云,何樵登. 1994. 基于Kjartansson模型的反Q滤波[J]. 地球物理学进展,9(1):90–99. Pei J Y,He Q D. 1994. Inverse Q filtering according to Kjartansson model[J]. Progress in Geophysics,9(1):90–99 (in Chinese).
王小杰,栾锡武. 2017. 基于小波分频技术的地层Q值补偿方法研究[J]. 石油物探,56(2):203–209. Wang X J,Luan X W. 2017. The study of formation Q value compensation method based on wavelet frequency division technology[J]. Geophysical Prospecting for Petroleum,56(2):203–209 (in Chinese).
熊晓军,贺振华,黄德济,陈学华. 2006. 广义S变换在地震高分辨处理中的应用[J]. 勘探地球物理进展,29(6):415–418. Xiong X J,He Z H,Huang D J,Chen X H. 2006. Application of generalized S transform in seismic high resolution processing[J]. Progress in Exploration Geophysics,29(6):415–418 (in Chinese).
姚振兴,高星,李维新. 2003. 用于深度域地震剖面衰减与频散补偿的反Q滤波方法[J]. 地球物理学报,46(2):229–233. Yao Z X,Gao X,Li W X. 2003. The forward Q method for compensating attenuation and frequency dispersion used in the seismic profile of depth domain[J]. Chinese Journal of Geophysics,46(2):229–233 (in Chinese).
张固澜,熊晓军,容娇君,张彦斌,蔡志东. 2010. 基于改进的广义S变换的地层吸收衰减补偿[J]. 石油地球物理勘探,45(4):512–515. Zhang G L,Xiong X J,Rong J J,Zhang Y B,Cai Z D. 2010. Stratum absorption and attenuation compensation based on improved generalized S-transform[J]. Oil Geophysical Prospecting,45(4):512–515 (in Chinese).
张固澜,林进,王熙明,贺振华,曹俊兴,张建军,贺锡雷,林凯,薛雅娟. 2015. 一种自适应增益限的反Q滤波[J]. 地球物理学报,58(7):2525–2535. Zhang G L,Lin J,Wang X M,He Z H,Cao J X,Zhang J J,He X L,Lin K,Xue Y J. 2015. A self-adaptive approach for inverse Q-filtering[J]. Chinese Journal of Geophysics,58(7):2525–2535 (in Chinese).
赵岩,毛宁波,陈旭. 2021. 基于时频域信噪比的自适应增益限反Q滤波方法[J]. 岩性油气藏,33(4):85–92. Zhao Y,Mao N B,Chen X. 2021. Self-adaptive gain-limit inverse Q filtering method based on SNR in time-frequency domain[J]. Lithologic Reservoirs,33(4):85–92 (in Chinese).
Braga I L S,Moraes F S. 2013. High-resolution gathers by inverse Q filtering in the wavelet domain[J]. Geophysics,78(2):V53–V61. doi: 10.1190/geo2011-0508.1
Hale D. 1981. An inverse Q-filter[J]. SEP Report,26:231–243.
Wang Y H. 2002. A stable and efficient approach of inverse Q filtering[J]. Geophysics,67(2):657–663. doi: 10.1190/1.1468627
Wang Y H. 2006. Inverse Q-filter for seismic resolution enhancement[J]. Geophysics,71(3):V51–V60. doi: 10.1190/1.2192912
Zhao Y,Mao N B,Xu J. 2019. Generalized stable inverse Q filtering[J]. J Appl Geophys,169:214–225. doi: 10.1016/j.jappgeo.2019.07.007
-
期刊类型引用(1)
1. 马明星，李剑，臧丹枫，曾援，郭陈莉. 基于衰减层析的地下浅层震动能量场重建方法. 电子测量技术. 2024(03): 42-47 . 
百度学术
其他类型引用(1)
计量
- 文章访问数: 293
- HTML全文浏览量: 155
- PDF下载量: 101
- 被引次数: 2
