CN112904412A - 一种矿山微震信号p波初至时刻提取方法及系统 - Google Patents

一种矿山微震信号p波初至时刻提取方法及系统 Download PDF

Info

Publication number
CN112904412A
CN112904412A CN201911221797.7A CN201911221797A CN112904412A CN 112904412 A CN112904412 A CN 112904412A CN 201911221797 A CN201911221797 A CN 201911221797A CN 112904412 A CN112904412 A CN 112904412A
Authority
CN
China
Prior art keywords
time
wave
arrival
microseismic signal
signal
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN201911221797.7A
Other languages
English (en)
Other versions
CN112904412B (zh
Inventor
石雅倩
张达
冀虎
戴锐
吴彦博
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
BGRIMM Technology Group Co Ltd
Original Assignee
BGRIMM Technology Group Co Ltd
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by BGRIMM Technology Group Co Ltd filed Critical BGRIMM Technology Group Co Ltd
Priority to CN201911221797.7A priority Critical patent/CN112904412B/zh
Publication of CN112904412A publication Critical patent/CN112904412A/zh
Application granted granted Critical
Publication of CN112904412B publication Critical patent/CN112904412B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • GPHYSICS
    • G01MEASURING; TESTING
    • G01VGEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28Processing seismic data, e.g. for interpretation or for event detection
    • G01V1/288Event detection in seismic signals, e.g. microseismics
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01VGEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28Processing seismic data, e.g. for interpretation or for event detection
    • G01V1/30Analysis
    • G01V1/307Analysis for determining seismic attributes, e.g. amplitude, instantaneous phase or frequency, reflection strength or polarity
    • YGENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02ATECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
    • Y02A90/00Technologies having an indirect contribution to adaptation to climate change
    • Y02A90/30Assessment of water resources

Landscapes

  • Engineering & Computer Science (AREA)
  • Remote Sensing (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Physics & Mathematics (AREA)
  • Environmental & Geological Engineering (AREA)
  • Acoustics & Sound (AREA)
  • Geology (AREA)
  • General Life Sciences & Earth Sciences (AREA)
  • General Physics & Mathematics (AREA)
  • Geophysics (AREA)
  • Business, Economics & Management (AREA)
  • Emergency Management (AREA)
  • Measurement Of Velocity Or Position Using Acoustic Or Ultrasonic Waves (AREA)
  • Measurement Of Mechanical Vibrations Or Ultrasonic Waves (AREA)

Abstract

本发明公开一种矿山微震信号P波初至时刻提取方法及系统,其方法包括:获取微震信号,并对所述微震信号进行升维,得到所述微震信号的时间‑频率‑能量分布图;根据所述时间‑频率‑能量分布图中信号与噪声能量分布,确定所述微震信号的主频带;在所述主频带上确定P波初至大致时刻;在所述微震信号中基于所述P波初至大致时刻内,选取时间窗,进一步确定P波初至精确时刻。本发明提供的方法能够充分利用微震信号的时域和频域的信息,通过能量分析来找到微震发生的真正时刻,从源头上解决P波初次定位产生的严重误差,进而拾取更准确的P波初至,消除突变噪声带来的影响。

Description

一种矿山微震信号P波初至时刻提取方法及系统
技术领域
本公开涉及微震监测技术领域,尤其涉及一种矿山微震信号P波初至时刻提取方法及系统。
背景技术
在地下矿山深部开采中,由于受到“三高一扰动”(高地应力、高温、高渗透压以及强烈的开采扰动)的恶劣环境影响,围岩变形、岩层塌陷、岩爆等动力灾害频发,严重威胁着矿山井下工作人员的安全。地压监测能够对矿井下岩层内部应力变化进行及时监测,在深部开采中占有极其重要的地位。矿山动力灾害发生过程中均伴随着岩石的破裂、裂纹的扩展及摩擦,这些过程均会产生量级较小的震动,称之为微震事件。而微震监测就是通过对物体震动进行监测,得到微震事件发生的时间、位置及震级,从而判断岩体承受压力大小、受到破坏程度等情况,以对其安全程度进行评价,是地压监测中的一种重要的监测手段。
微震监测系统能够为矿山进行安全预警、控制地压灾害提供有效依据,其主要包括微震数据采集系统和微震处理分析系统两大部分。首先微震数据采集系统对现场数据进行实时采集,再将采集到的信号传送到微震处理分析系统中。与地震信号类似,微震信号包括P波和S波。在微震处理分析系统中,必须要先拾取该段采集信号中的微震信号波形的P波初至,再通过对微震信号发生时刻与能量的分析来得到微震信号发生源,最后通过对多个微震信号的统计分析,得到整体矿山压力变化分析,进而指导矿山安全作业,为矿山提供安全预警。因此P波初至是否足够精确的问题在微震监测系统中至关重要,只有更准确地找到P波初至时刻,才能更精准地确定岩石破裂所在的位置与时间,进而判断井下岩体稳定性及矿山风险水平。
现有技术中大多是采用基于时域信号本身变化进行拾取的,其中最为常见也最为基础的P波初至方法是STA(信号短时平均值)/LTA(信号长时平均值)长短时窗法,该方法通过设定一定的阈值来进行震相拾取。但是,该类方法仅从信号在时域上发生的变化来寻找P波初至,忽略了信号在频域上的特征,对信噪比低的信号效果较差,使用场景存在一定的局限性。此外,基于STA/LTA的P波拾取方法会受到长短时间窗长选取、拾取阈值的选取的影响,尤其受突变噪声引起的信号波动的影响较大。因此在井下低信噪比、未知噪声干扰环境下,实现更高精度的P波初至拾取是非常重要的。
发明内容
(一)要解决的技术问题
为了解决现有技术的上述问题,本公开提供一种高精度的矿山微震信号P波初至时刻提取方法及系统,解决现有技术中无法获取P波初至的精确时刻的问题。
(二)技术方案
为了达到上述目的,本公开采用的主要技术方案包括:
一方面,本公开一实施例提供一种矿山微震信号P波初至时刻提取方法,其包括:
步骤S210:获取微震信号,并对所述微震信号进行升维,得到所述微震信号的时间-频率-能量分布图;
步骤S220:根据所述时间-频率-能量分布图中信号与噪声能量分布,确定所述微震信号的主频带;
步骤S230:在所述主频带上确定P波初至大致时刻;
步骤S240:在所述微震信号中基于所述P波初至大致时刻内,选取时间窗,进一步确定P波初至精确时刻。
本公开的一个实施例中,所述步骤S210中获取微震信号,并对所述微震信号进行升维,得到所述微震信号的时间-频率-能量分布图包括:
步骤S301:通过传感器对矿山进行监测,获取所述微震信号,并将其设为f(t);
步骤S302:将所述微震信号在时域上表示为多个谐波分量的叠加形式:
Figure BDA0002301049650000031
其中Ak(t)为第k个谐波分量的瞬时振幅,θk(t)为第k个谐波分量的瞬时相位,e(t)为噪声或误差,K为可分解分量的个数;
步骤S303:利用同步挤压小波变换方法对所述微震信号从时域和频域进行分解,得到所述微震信号的时间-频率-能量分布图。
本公开的一个实施例中,所述步骤S303中利用同步挤压小波变换方法对所述微震信号从时域和频域进行分解,得到所述微震信号的时间-频率-能量分布图包括:
利用同步挤压变换结合所述频域小波变换,将时间-尺度平面的能量重新分配到时间-频率平面,得到所述时间-频率-能量分布图。
本公开的一个实施例中,所述步骤S220中根据所述时间-频率-能量分布图中信号与噪声能量分布,确定所述微震信号的主频带包括:
步骤S601:对所述时间-频率-能量分布图的时域部分进行叠加,得到基于频域的能量曲线图;
步骤S602:针对所述基于频域的能量曲线图中的数据利用预设阈值进行选取,得到所述主频带。
本公开的一个实施例中,所述步骤S602中针对所述基于频域的能量曲线图中的数据利用预设阈值进行选取,得到所述主频带包括:
针对所述基于频域的能量曲线图中的数据按照从小到大重新排列;
对重新排列后的数据按照所述预设阈值的位置所在的数据值进行选取,得到所述主频带,其中所述预设阈值的范围为70%~90%。
本公开的一个实施例中,所述步骤S230中在所述主频带上确定P波初至大致时刻包括:
步骤S801:在一个时间点上对所述主频带间的能量进行求和,得到基于时域的主频带的能量曲线图;
步骤S802:对所述基于时域的主频带的能量曲线图进行滑动平均处理,并对处理后的曲线图采用AIC算法计算AIC值得到AIC曲线;
步骤S803:在所述AIC曲线上多个极小值中选择第一个极小值所在的时刻,作为所述P波初至大致时刻。
本公开的一个实施例中,所述步骤S340中在所述微震信号中基于所述P波初至大致时刻内,选取时间窗,进一步确定P波初至精确时刻包括:
步骤A1:根据所述P波初至大致时刻设置时间窗;
步骤A2:对所述微震信号在所述时间窗内利用AIC算法计算AIC值;
步骤A3:选择所述时间窗内AIC值最小的位置为所述P波初至精确时刻。
本公开的一个实施例中,所述步骤A1中根据所述P波初至大致时刻设置时间窗包括:
根据所述P波初至大致时刻结合预设间隔设置时间窗,所述时间窗为:
[t1-2*int,t1+2*int]
其中t1为所述P波初至大致时刻,int为预设间隔。
另一方面,本公开另一实施例还提供一种矿山微震信号P波初至时刻提取系统,其包括:
分解模块,用于获取微震信号,并对所述微震信号进行升维,得到所述微震信号的时间-频率-能量分布图;
主频带模块,用于根据所述时间-频率-能量分布图中信号与噪声能量分布的,确定所述微震信号的主频带;
粗提取模块,用于在所述主频带上确定P波初至大致时刻;
精确提取模块,用于在所述微震信号中基于所述P波初至大致时刻内,选取时间窗,确定P波初至精确时刻。
本公开的一个实施例中,所述升维模块包括:
数据获取子模块,用于通过传感器对矿山进行监测,获取所述微震信号,并将其设为f(t);
分解子模块,用于将所述微震信号在时域上表示为多个谐波分量的叠加形式:
Figure BDA0002301049650000051
其中Ak(t)为第k个谐波分量的瞬时振幅,θk(t)为第k个谐波分量的瞬时相位,e(t)为噪声或误差,K为可分解分量的个数;
变换子模块,用于利用同步挤压小波变换方法对所述微震信号从时域和频域进行分解,得到所述微震信号的时间-频率-能量分布图。
本公开的一个实施例中,所述主频带模块包括:
叠加子模块,用于对所述时间-频率-能量分布图的时域部分进行叠加,得到基于频域的能量曲线图;
选取子模块,用于针对所述基于频域的能量曲线图中的数据利用预设阈值进行选取,得到所述主频带。
本公开的一个实施例中,所述选取子模块包括:
排序单元,用于针对所述基于频域的能量曲线图中的数据按照从小到大重新排列;
选取单元,用于对重新排列后的数据按照所述预设阈值为80%的位置所在的数据值进行选取,得到所述主频带。
本公开的一个实施例中,所述粗提取模块包括:
求和子模块,用于在一个时间点上对所述主频带间的能量进行求和,得到基于时域的主频带的能量曲线图;
计算子模块,用于对所述基于时域的主频带的能量曲线图进行滑动平均处理,并对处理后的曲线图采用AIC算法计算AIC值得到AIC曲线;
大致时刻确定子模块,用于在所述AIC曲线上多个极小值中选择第一个极小值所在的时刻,作为所述P波初至大致时刻。
本公开的一个实施例中,所述精确提取模块包括:
时间窗设置子模块,用于根据所述P波初至大致时刻设置时间窗;
时间窗内计算子模块,用于对所述微震信号在所述时间窗内利用AIC算法计算AIC值;
精确时刻确定子模块,用于选择所述时间窗内AIC值最小的位置为所述P波初至精确时刻。
(三)有益效果
本公开的有益效果是:本发明提供的一种矿山微震信号P波初至时刻提取方法及系统,同时分析原始的微震信号的时域和频域所对应的信息,通过自动搜索微震信号能量集中的主频带,再从时域上在时频谱中找到微震信号P波初至的大致起始时间,最后回到微震信号在时间域中的变化来进一步确定P波初至的精确时刻。该方法能够充分利用微震信号的时域和频域的信息,通过能量分析来找到微震发生的真正时刻,从源头上解决P波初次定位产生的严重误差,进而拾取更准确的P波初至,消除突变噪声带来的影响。
附图说明
图1为本公开相关实施例中P波拾取的结果示意图;
图2为本公开一实施例提供的一种矿山微震信号P波初至时刻提取方法的流程图;
图3为本发明一实施例图2中步骤S210的流程图;
图4为本发明一实施例中微震信号的波形图;
图5为本发明一实施例中图4的波形变换后得到的时间-频率-能量图;
图6为本发明一实施例图2中步骤S220的流程图;
图7为本发明一实施例中得到的基于频域的能量曲线图;
图8为本发明一实施例图2中步骤S230的流程图;
图9为本发明一实施例中进行压缩的示意图;
图10为本发明一实施例中基于时域的主频带的能量曲线图;
图11为本发明一实施例中得到的AIC曲线示意图;
图12为本发明一实施例中对于微震信号提取P波初至时刻的结果示意图;
图13为本发明一实施例中实现上述方法的信号处理流程示意图;
图14为本发明另一实施例中提供的矿山微震信号P波初至时刻提取系统的示意图;
图15为本发明另一实施例图14中升维模块1410的示意图;
图16为本发明另一实施例图14中主频带模块1420的示意图;
图17为本发明另一实施例图14中粗提取模块1430的示意图;
图18为本发明另一实施例图14中精确提取模块1440的示意图;
图19为本发明中提供的适于用来实现本发明实施例的电子设备的计算机系统的结构示意图。
具体实施方式
为了更好的解释本公开,以便于理解,下面结合附图,通过具体实施方式,对本公开作详细描述。
本文所使用的所有的技术和科学术语与属于本公开的技术领域的技术人员通常理解的含义相同。本文中在本公开的说明书中所使用的术语只是为了描述具体的实施例的目的,不是旨在于限制本公开。本文所使用的术语“和/或”包括一个或多个相关的所列项目的任意的和所有的组合。
针对现有技术中的存在的不足,本公开提供一种矿山微震信号P波初至时刻提取方法,基于时间-频率-能量三维空间分析来提高矿山微震信号P波初至拾取精度的方法,采用自适应方法,自动搜索信号主频带,能够有效适用于不同矿山不同季节下的各类特征波形,确定信号发生的频带,进而能够准确拾取P波初至大致时刻,避免P波初至二次定位产生的严重误差。
图2为本公开一实施例提供的一种矿山微震信号P波初至时刻提取方法的流程图,如图2所示,该方法包括以下步骤:
如图2所示,在步骤S210中,获取微震信号,并对所述微震信号进行升维,得到所述微震信号的时间-频率-能量分布图;
如图2所示,在步骤S220中,根据所述时间-频率-能量分布图中信号与噪声能量分布,确定所述微震信号的主频带;
如图2所示,在步骤S230中,在所述主频带上确定P波初至大致时刻;
如图2所示,在步骤S240中,在所述微震信号中基于所述P波初至大致时刻内,选取时间窗,进一步确定P波初至精确时刻。
基于上述,该方法通过设置一些固定的参数外全程不需要人工参与,比如信号主频带,不需要人工参与查看信号等就可以自动找到主频带的位置,包括后面P波初至拾取,均没有人工参与,因此是一种自动搜索信号、自动提取P波初至的方法,而且该方法比其他自动提取P波初至的精度更高。
以下对图2所示实施例的各个步骤的具体实现进行详细阐述:
在步骤S210中,获取微震信号,并对所述微震信号进行升维,得到所述微震信号的时间-频率-能量分布图。
在本公开的一个实施例中,该步骤中首先获取微震信号,并对所述微震信号进行升维,得到所述微震信号的时间-频率-能量分布图。由于某种情况下在低维空间无法有效分类,当映射到高维时可以进行好的分类,因此本实施例中“升维”为将信号从低维空间映射到高维空间,将微震信号从时频域二维分析升维到时域、频域和能量域的三维分析。由于有效震动信号的特征为能量分布集中,能量值较大,并且该信号仅发生在一定的时间段内;噪声信号的特征多为能量分布分散、能量值普遍较小,部分噪声为连续信号,遍布整个时域,本实施例中通过升维,能够清楚看到信号发生时能量的变化,从而便于后续步骤确定能量集中的部分。
图3为本发明一实施例图2中步骤S210的流程图,具体包括以下步骤:
如图3所示,在步骤S301中,通过传感器对矿山进行监测,获取所述微震信号,并将其设为f(t)。具体为:接收采集到微震信号的波形f(t),其数据序列长度为N,采样率为λ,数据时间长度为N/λ秒。
如图3所示,在步骤S302中,将所述微震信号在时域上表示为多个谐波分量的叠加形式:
Figure BDA0002301049650000091
其中Ak(t)为第k个谐波分量的瞬时振幅,θk(t)为第k个谐波分量的瞬时相位,e(t)为噪声或误差,K为可分解分量的个数。
如图3所示,在步骤S303中,利用同步挤压小波变换方法对所述微震信号从时域和频域进行分解,得到所述微震信号的时间-频率-能量分布图。
在本公开的一个实施例中,该步骤中包括:
首先,对所述微震信号在时域进行连续小波变换,得到小波系数;然后,将所述小波系数变换到频域,得到频域小波系数;最后,利用同步挤压变换结合所述频域小波变换,将时间-尺度平面的能量重新分配到时间-频率平面,得到所述时间-频率-能量分布图。
例如,对矿用微震信号f(t)进行连续小波变换(CWT),得到的小波系数Wf(a,b):
Figure BDA0002301049650000092
式中a为尺度因子,b为平移因子;ψ*为共轭小波函数。其在频率域可等价变换为:
Figure BDA0002301049650000093
式中ξ代表频率;
Figure BDA0002301049650000094
分别表示为f(t)、ψ(t)的傅里叶变换。
利用同步挤压所变换,对其时间-尺度平面的能量重新分配成为时间-频率平面,此压缩变换Tfl,b)的离散公式为:
Figure BDA0002301049650000101
其中,(△a)i=ai-ai-1
Figure BDA0002301049650000102
ωl为中心频率。
这样,经过上述同步挤压小波变换后,微震信号可以被分解为时间-频率-能量图。以新疆某矿采集到的微震信号为例,图4为本发明一实施例中微震信号的波形图,图5为本发明一实施例中图4的波形变换后得到的时间-频率-能量图。
在步骤S220中,根据所述时间-频率-能量分布图中信号与噪声能量分布的特点,自适应确定所述微震信号的主频带。
图6为本发明一实施例图2中步骤S220的流程图,具体包括以下步骤:
如图6所示,在步骤S601中,对所述时间-频率-能量分布图的时域部分进行叠加,得到基于频域的能量曲线图。
图7为本发明一实施例中得到的基于频域的能量曲线图,如图7所示,其中两线W1和W2之间表示该信号频率区间范围为25~520Hz。
如图6所示,在步骤S602中,针对所述基于频域的能量曲线图中的数据利用预设阈值进行选取,得到所述主频带。
在本公开的一个实施例中,该步骤中具体为针对所述基于频域的能量曲线图中的数据按照从小到大重新排列;对重新排列后的数据按照所述预设阈值的位置所在的数据值进行选取,其中所述预设阈值范围为70%~90%。例如预设阈值可以选择80%,仍以图4所示波形图为例,其数据采集系统中设置所采集到的波形都较为完整,且长度较短。阈值选择方法为:对基于频域的能量曲线数据按数据值从小到大重新排列,根据数采预先设定的采集方法及数据长度限制,本实施例选择其80%位置所在的数据值作为阈值来选取主频带。
在步骤S230中,在所述主频带上确定P波初至大致时刻。
图8为本发明一实施例图2中步骤S230的流程图,具体包括以下步骤:
如图8所示,在步骤S801中,在一个时间点上对所述主频带间的能量进行求和,得到基于时域的主频带的能量曲线图。
该步骤中对于所选频带(主频带)内的数据,按照时间尺度,对该频带内能量求和,得到一条基于时域的主频带的能量曲线图,图9为本发明一实施例中进行压缩的示意图,即对图5所示的主频带进行压缩后得到图9所示示意图,其中W3代表一个时间点对应的频带间所有能量,对W3上所有的能量值进行求和,得到该时间点的能量值,每个时间点对应求出一个值后,得到一条基于时间的能量曲线,如图10所示,图10为本发明一实施例中基于时域的主频带的能量曲线图。
如图8所示,在步骤S802中,对所述基于时域的主频带的能量曲线图进行滑动平均处理,并对处理后的曲线图采用AIC算法计算AIC值得到AIC曲线。
该步骤中利用赤池信息量准则(Akaike information criterion,简称AIC)方法进行处理。另外,为了剔除突变噪声和毛刺的影响,设置预设间隔int,对基于时域的能量曲线图进行滑动平均处理,本实施例中预设间隔为:
Figure BDA0002301049650000111
其中N为数据长度。
接着对处理后的能量曲线求其AIC值,AIC方法的基本原理是求取地震信号AIC函数局部最小值的位置。AIC函数定义为:
AIC(k)=klog(var(x[1,k])+(N-k-1)log(var(x[k+1,N]))) 公式(8)
其中x[i],i=1,2,…,N为矿用微震信号f(t)的离散序列,k为阶数。
如图8所示,在步骤S803中,在所述AIC曲线上多个极小值中选择第一个极小值所在的时刻为所述P波初至大致时刻。图11为本发明一实施例中得到的AIC曲线示意图,如图11所示,找出AIC曲线中的第一个极小值点t,即图11中S点所标志的位置,则P波初至大致时刻t1为:
t1=t*int 公式(9)
本实施例中结合时频域进行分析,所能得到的信息要比其他方法直接在时域上观察信号变化得到的信息要多。从时频域上可以通过选取主频带来减少一些突变、毛刺的噪声的干扰,从而找到特征信号准确发生时间。
在步骤S240中,在所述微震信号中基于所述P波初至大致时刻内,选取时间窗,进一步确定P波初至精确时刻。
在本公开的一个实施例中,该步骤具体包括:
步骤A1:根据所述P波初至大致时刻设置时间窗;
步骤A2:对所述微震信号在所述时间窗内利用AIC算法计算AIC值;
步骤A3:选择所述时间窗内AIC值最小的位置为所述P波初至精确时刻。
在本公开的一个实施例中,上述步骤中利用AIC算法确定P波初至精确时刻过程中需要在P波初至时刻附近选择一个合适的时窗进行局部计算。本实施例中根据所述P波初至大致时刻结合预设间隔设置时间窗时,所述时间窗为:
[t1-2*int,t1+2*int] 公式(10)
其中t1为所述P波初至大致时刻,int为预设间隔。
图12为本发明一实施例中对于微震信号提取P波初至时刻的结果示意图,如图12所示,其中W4为采用STA/LTA方法得到的P波初至时刻,W5为采用本发明方法得到的P波初至时刻,根据图12所示,采用本发明实施例提供的方法利用AIC算法对信号趋势的准确辨别的优点,两次采用AIC算法,第一次利用AIC算法确定P波初至大致时刻,进而在该时刻设置时间窗,在时间窗内第二次利用AIC算法确定P波精确时刻,分别从能量变化趋势上和时域信号变化中搜索P波位置,进而得到高精度的P波初至时刻。
图13为本发明一实施例中实现上述方法的信号处理流程示意图,如图13所示:
步骤S1301,输入信号。
步骤S1302,同步挤压小波变化。
步骤S1303,基于频域进行能量压缩。
步骤S1304,对于给定频域段,基于时域进行能量压缩。
步骤S1305,初次AIC计算得到P波初至大致时刻(即粗略P波位置)。
步骤S1306,基于给定的时间窗,基于原始微震信号进行二次AIC计算得到P波初至精确时刻(即精确P波位置)。
需要说明的是,在本实施例中,第一次是对滑动平均后的基于时域的能量曲线图做AIC算法,第二次是对时窗内部分的原始信号做AIC算法,在本发明其他实施例中,还可以利用其他算法代替AIC算法。
综上所述,采用本发明基于时间-频率-能量三维空间分析来提高矿山微震信号P波初至拾取精度的方法。该方法通过对信号同时从时域和频域进行能量分析的方法来提高有效信号的P波初至拾取精度,是一种时域和频域相结合方法,有效避免了因未知噪声干扰而导致的微震信号发生位置定位不准确的问题,进而提高P波初至的拾取精度。针对STA/LTA法中对波形数据长度的限制问题,本发明未采用STA/LTA长短时窗比的方法,避免因选取时窗对采集数据时长的限制,能够自适应于不同长短的波形并有效进行处理,更具有普适性。另外,针对大部分搜索信号的方法均针对固定频带存在的问题,本发明采用自适应方法,自动搜索信号主频带,能够有效适用于不同矿山不同季节下的各类特征波形,确定信号发生的频带,进而能够准确拾取P波初至大致时刻,避免P波初至二次定位产生的严重误差。
另一方面,本发明另一实施例中还提供一种矿山微震信号P波初至时刻提取系统,图14为本发明另一实施例中提供的矿山微震信号P波初至时刻提取系统的示意图,如图14所示,该系统1400包括:升维模块1410、主频带模块1420、粗提取模块1430和精确提取模块1440。
升维模块1410用于获取微震信号,并对所述微震信号进行升维,得到所述微震信号的时间-频率-能量分布图;主频带模块1420用于根据所述时间-频率-能量分布图中信号与噪声能量分布,确定所述微震信号的主频带;粗提取模块1430用于在所述主频带上确定P波初至大致时刻;精确提取模块1440用于在所述微震信号中基于所述P波初至大致时刻内,选取时间窗,进一步确定P波初至精确时刻。
在本公开其中一个实施例中,图15为本发明另一实施例图14中升维模块1410的示意图,升维模块1410中进一步包括:数据获取子模块1411、分解子模块1412和变换子模块1413。
数据获取子模块1411用于通过传感器对矿山进行监测,获取所述微震信号,并将其设为f(t);
分解子模块1412用于将所述微震信号在时域上表示为多个谐波分量的叠加形式:
Figure BDA0002301049650000141
其中Ak(t)为第k个谐波分量的瞬时振幅,θk(t)为第k个谐波分量的瞬时相位,e(t)为噪声或误差,K为可分解分量的个数;
变换子模块1413用于利用同步挤压小波变换方法对所述微震信号从时域和频域进行分解,得到所述微震信号的时间-频率-能量分布图。
其中变换子模块1413具体为:利用同步挤压变换结合所述频域小波变换,将时间-尺度平面的能量重新分配到时间-频率平面,得到所述时间-频率-能量分布图。
在本公开其中一个实施例中,图16为本发明另一实施例图14中主频带模块1420的示意图,主频带模块1420进一步包括:叠加子模块1421和选取子模块1422。其中叠加子模块1421用于对所述时间-频率-能量分布图的时域部分进行叠加,得到基于频域的能量曲线图;选取子模块1422用于针对所述基于频域的能量曲线图中的数据利用预设阈值进行选取,得到所述主频带。
其中选取子模块包括:排序单元,用于针对所述基于频域的能量曲线图中的数据按照从小到大重新排列;选取单元,用于对重新排列后的数据按照所述预设阈值为80%的位置所在的数据值进行选取,得到所述主频带。
在本公开其中一个实施例中,图17为本发明另一实施例图14中粗提取模块1430的示意图,粗提取模块1430进一步包括:求和子模块1431、计算子模块1432和大致时刻确定子模块1433。
其中求和子模块1431用于在一个时间点上对所述主频带间的能量进行求和,得到基于时域的主频带的能量曲线图;计算子模块1432用于对所述基于时域的主频带的能量曲线图进行滑动平均处理,并对处理后的曲线图采用AIC算法计算AIC值得到AIC曲线;大致时刻确定子模块1433用于在所述AIC曲线上多个极小值中选择第一个极小值所在的时刻,作为所述P波初至大致时刻。
在本公开其中一个实施例中,图18为本发明另一实施例图14中精确提取模块1440的示意图,精确提取模块1440进一步包括:时间窗设置子模块1441、时间窗内计算子模块1442和精确时刻确定子模块1443。
时间窗设置子模块1441用于根据所述P波初至大致时刻设置时间窗;时间窗内计算子模块1442用于对所述微震信号在所述时间窗内利用AIC算法计算AIC值;精确时刻确定子模块1443用于选择所述时间窗内AIC值最小的位置为所述P波初至精确时刻。由于本公开的示例实施例的矿山微震信号P波初至时刻提取系统的各个功能模块与上述图2所示的矿山微震信号P波初至时刻提取方法的示例实施例的步骤对应,因此对于本公开装置实施例中未披露的细节,请参照本公开上述的矿山微震信号P波初至时刻提取方法的实施例。
下面参考图19,其示出了适于用来实现本发明实施例的电子设备的计算机系统1900的结构示意图。图19示出的电子设备的计算机系统1900仅是一个示例,不应对本发明实施例的功能和使用范围带来任何限制。
如图19所示,计算机系统1900包括中央处理单元(CPU)1901,其可以根据存储在只读存储器(ROM)1902中的程序或者从存储部分1908加载到随机访问存储器(RAM)1903中的程序而执行各种适当的动作和处理。在RAM 1903中,还存储有系统操作所需的各种程序和数据。CPU 1901、ROM 1902以及RAM 1903通过总线1904彼此相连。输入/输出(I/O)接口1905也连接至总线1904。
以下部件连接至I/O接口1905:包括键盘、鼠标等的输入部分1906;包括诸如阴极射线管(CRT)、液晶显示器(LCD)等以及扬声器等的输出部分1907;包括硬盘等的存储部分1908;以及包括诸如LAN卡、调制解调器等的网络接口卡的通信部分1909。通信部分1909经由诸如因特网的网络执行通信处理。驱动器1910也根据需要连接至I/O接口1905。可拆卸介质1911,诸如磁盘、光盘、磁光盘、半导体存储器等等,根据需要安装在驱动器1910上,以便于从其上读出的计算机程序根据需要被安装入存储部分1908。
特别地,根据本发明的实施例,上文参考流程图描述的过程可以被实现为计算机软件程序。例如,本发明的实施例包括一种计算机程序产品,其包括承载在计算机可读介质上的计算机程序,该计算机程序包含用于执行流程图所示的方法的程序代码。在这样的实施例中,该计算机程序可以通过通信部分1909从网络上被下载和安装,和/或从可拆卸介质1911被安装。在该计算机程序被中央处理单元(CPU)1901执行时,执行本申请的系统中限定的上述功能。
需要说明的是,本发明所示的计算机可读介质可以是计算机可读信号介质或者计算机可读存储介质或者是上述两者的任意组合。计算机可读存储介质例如可以是——但不限于——电、磁、光、电磁、红外线、或半导体的系统、装置或器件,或者任意以上的组合。计算机可读存储介质的更具体的例子可以包括但不限于:具有一个或多个导线的电连接、便携式计算机磁盘、硬盘、随机访问存储器(RAM)、只读存储器(ROM)、可擦式可编程只读存储器(EPROM或闪存)、光纤、便携式紧凑磁盘只读存储器(CD-ROM)、光存储器件、磁存储器件、或者上述的任意合适的组合。在本发明中,计算机可读存储介质可以是任何包含或存储程序的有形介质,该程序可以被指令执行系统、装置或者器件使用或者与其结合使用。而在本发明中,计算机可读的信号介质可以包括在基带中或者作为载波一部分传播的数据信号,其中承载了计算机可读的程序代码。这种传播的数据信号可以采用多种形式,包括但不限于电磁信号、光信号或上述的任意合适的组合。计算机可读的信号介质还可以是计算机可读存储介质以外的任何计算机可读介质,该计算机可读介质可以发送、传播或者传输用于由指令执行系统、装置或者器件使用或者与其结合使用的程序。计算机可读介质上包含的程序代码可以用任何适当的介质传输,包括但不限于:无线、电线、光缆、RF等等,或者上述的任意合适的组合。
附图中的流程图和框图,图示了按照本发明各种实施例的系统、方法和计算机程序产品的可能实现的体系架构、功能和操作。在这点上,流程图或框图中的每个方框可以代表一个模块、程序段、或代码的一部分,上述模块、程序段、或代码的一部分包含一个或多个用于实现规定的逻辑功能的可执行指令。也应当注意,在有些作为替换的实现中,方框中所标注的功能也可以以不同于附图中所标注的顺序发生。例如,两个接连地表示的方框实际上可以基本并行地执行,它们有时也可以按相反的顺序执行,这依所涉及的功能而定。也要注意的是,框图或流程图中的每个方框、以及框图或流程图中的方框的组合,可以用执行规定的功能或操作的专用的基于硬件的系统来实现,或者可以用专用硬件与计算机指令的组合来实现。
描述于本发明实施例中所涉及到的单元可以通过软件的方式实现,也可以通过硬件的方式来实现,所描述的单元也可以设置在处理器中。其中,这些单元的名称在某种情况下并不构成对该单元本身的限定。
作为另一方面,本申请还提供了一种计算机可读介质,该计算机可读介质可以是上述实施例中描述的电子设备中所包含的;也可以是单独存在,而未装配入该电子设备中。上述计算机可读介质承载有一个或者多个程序,当上述一个或者多个程序被一个该电子设备执行时,使得该电子设备实现如上述实施例中所述的矿山微震信号P波初至时刻提取方法。
例如,所述的电子设备可以实现如图2中所示的:在步骤S210中,获取微震信号,并对所述微震信号进行升维,得到所述微震信号的时间-频率-能量分布图;在步骤S220中,根据所述时间-频率-能量分布图中信号与噪声能量分布,确定所述微震信号的主频带;在步骤S230中,在所述主频带上确定P波初至大致时刻;在步骤S240中,在所述微震信号中基于所述P波初至大致时刻内,选取时间窗,进一步确定P波初至精确时刻。
应当注意,尽管在上文详细描述中提及了用于动作执行的设备的若干模块或者单元,但是这种划分并非强制性的。实际上,根据本发明的实施方式,上文描述的两个或更多模块或者单元的特征和功能可以在一个模块或者单元中具体化。反之,上文描述的一个模块或者单元的特征和功能可以进一步划分为由多个模块或者单元来具体化。
通过以上的实施方式的描述,本领域的技术人员易于理解,这里描述的示例实施方式可以通过软件实现,也可以通过软件结合必要的硬件的方式来实现。因此,根据本发明实施方式的技术方案可以以软件产品的形式体现出来,该软件产品可以存储在一个非易失性存储介质(可以是CD-ROM,U盘,移动硬盘等)中或网络上,包括若干指令以使得一台计算设备(可以是个人计算机、服务器、触控终端、或者网络设备等)执行根据本发明实施方式的方法。
应当注意,尽管在上文详细描述中提及了用于动作执行的设备的若干模块或者单元,但是这种划分并非强制性的。实际上,根据本公开的实施方式,上文描述的两个或更多模块或者单元的特征和功能可以在一个模块或者单元中具体化。反之,上文描述的一个模块或者单元的特征和功能可以进一步划分为由多个模块或者单元来具体化。
本领域技术人员在考虑说明书及实践这里公开的发明后,将容易想到本公开的其它实施方案。本申请旨在涵盖本公开的任何变型、用途或者适应性变化,这些变型、用途或者适应性变化遵循本公开的一般性原理并包括本公开未公开的本技术领域中的公知常识或惯用技术手段。说明书和实施例仅被视为示例性的,本公开的真正范围和精神由下面的权利要求指出。
应当理解的是,本公开并不局限于上面已经描述并在附图中示出的精确结构,并且可以在不脱离其范围进行各种修改和改变。本公开的范围仅由所附的权利要求来限制。

Claims (10)

1.一种矿山微震信号P波初至时刻提取方法,其特征在于,其包括:
步骤S210:获取微震信号,并对所述微震信号进行升维,得到所述微震信号的时间-频率-能量分布图;
步骤S220:根据所述时间-频率-能量分布图中信号与噪声能量分布,确定所述微震信号的主频带;
步骤S230:在所述主频带上,确定P波初至大致时刻;
步骤S240:在所述微震信号中,基于所述P波初至大致时刻,选取时间窗,进一步确定P波初至精确时刻。
2.如权利要求1所述的矿山微震信号P波初至时刻提取方法,其特征在于,所述步骤S210中获取微震信号,并对所述微震信号进行升维,得到所述微震信号的时间-频率-能量分布图包括:
步骤S301:通过传感器对矿山进行监测,获取所述微震信号,并将其设为f(t);
步骤S302:将所述微震信号在时域上表示为多个谐波分量的叠加形式:
Figure FDA0002301049640000011
其中Ak(t)为第k个谐波分量的瞬时振幅,θk(t)为第k个谐波分量的瞬时相位,e(t)为噪声或误差,K为可分解分量的个数;
步骤S303:利用同步挤压小波变换方法对所述微震信号从时域和频域进行分解,得到所述微震信号的时间-频率-能量分布图。
3.如权利要求2所述的矿山微震信号P波初至时刻提取方法,其特征在于,所述步骤S303中利用同步挤压小波变换方法对所述微震信号从时域和频域进行分解,得到所述微震信号的时间-频率-能量分布图包括:
利用同步挤压变换结合所述频域小波变换,将时间-尺度平面的能量重新分配到时间-频率平面,得到所述时间-频率-能量分布图。
4.如权利要求1所述的矿山微震信号P波初至时刻提取方法,其特征在于,所述步骤S220中根据所述时间-频率-能量分布图中信号与噪声能量分布,确定所述微震信号的主频带包括:
步骤S601:对所述时间-频率-能量分布图的时域部分进行叠加,得到基于频域的能量曲线图;
步骤S602:针对所述基于频域的能量曲线图中的数据利用预设阈值进行选取,得到所述主频带。
5.如权利要求4所述的矿山微震信号P波初至时刻提取方法,其特征在于,所述步骤S602中针对所述基于频域的能量曲线图中的数据利用预设阈值进行选取,得到所述主频带包括:
针对所述基于频域的能量曲线图中的数据按照从小到大重新排列;
对重新排列后的数据按照所述预设阈值的位置所在的数据值进行选取,得到所述主频带,其中所述预设阈值的范围为70%~90%。
6.如权利要求1所述的矿山微震信号P波初至时刻提取方法,其特征在于,所述步骤S230中在所述主频带上确定P波初至大致时刻包括:
步骤S801:在一个时间点上对所述主频带间的能量进行求和,得到基于时域的主频带的能量曲线图;
步骤S802:对所述基于时域的主频带的能量曲线图进行滑动平均处理,并对处理后的曲线图采用AIC算法计算AIC值得到AIC曲线;
步骤S803:在所述AIC曲线上多个极小值中选择第一个极小值所在的时刻,作为所述P波初至大致时刻。
7.如权利要求1所述的矿山微震信号P波初至时刻提取方法,其特征在于,所述步骤S340中在所述微震信号中,基于所述P波初至大致时刻,选取时间窗,进一步确定P波初至精确时刻包括:
步骤A1:根据所述P波初至大致时刻设置时间窗;
步骤A2:对所述微震信号,在所述时间窗内,利用AIC算法计算AIC值;
步骤A3:选择所述时间窗内AIC值最小的位置为所述P波初至精确时刻。
8.如权利要求7所述的矿山微震信号P波初至时刻提取方法,其特征在于,所述步骤A1中根据所述P波初至大致时刻设置时间窗包括:
根据所述P波初至大致时刻结合预设间隔设置时间窗,所述时间窗为:
[t1-2*int,t1+2*int]
其中t1为所述P波初至大致时刻,int为预设间隔。
9.一种矿山微震信号P波初至时刻提取系统,其特征在于,包括:
升维模块,用于获取微震信号,并对所述微震信号进行升维,得到所述微震信号的时间-频率-能量分布图;
主频带模块,用于根据所述时间-频率-能量分布图中信号与噪声能量分布,确定所述微震信号的主频带;
粗提取模块,用于在所述主频带上确定P波初至大致时刻;
精确提取模块,用于在所述微震信号中基于所述P波初至大致时刻内,选取时间窗,进一步确定P波初至精确时刻。
10.如权利要求9所述的矿山微震信号P波初至时刻提取系统,其特征在于,所述升维模块包括:
数据获取子模块,用于通过传感器对矿山进行监测,获取所述微震信号,并将其设为f(t);
分解子模块,用于将所述微震信号在时域上表示为多个谐波分量的叠加形式:
Figure FDA0002301049640000031
其中Ak(t)为第k个谐波分量的瞬时振幅,θk(t)为第k个谐波分量的瞬时相位,e(t)为噪声或误差,K为可分解分量的个数;
变换子模块,用于利用同步挤压小波变换方法对所述微震信号从时域和频域进行分解,得到所述微震信号的时间-频率-能量分布图。
CN201911221797.7A 2019-12-03 2019-12-03 一种矿山微震信号p波初至时刻提取方法及系统 Active CN112904412B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201911221797.7A CN112904412B (zh) 2019-12-03 2019-12-03 一种矿山微震信号p波初至时刻提取方法及系统

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201911221797.7A CN112904412B (zh) 2019-12-03 2019-12-03 一种矿山微震信号p波初至时刻提取方法及系统

Publications (2)

Publication Number Publication Date
CN112904412A true CN112904412A (zh) 2021-06-04
CN112904412B CN112904412B (zh) 2024-02-23

Family

ID=76104288

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201911221797.7A Active CN112904412B (zh) 2019-12-03 2019-12-03 一种矿山微震信号p波初至时刻提取方法及系统

Country Status (1)

Country Link
CN (1) CN112904412B (zh)

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115016001A (zh) * 2022-08-09 2022-09-06 北京科技大学 一种基于维度变换的微震信号到时自动拾取方法
CN116047604A (zh) * 2023-02-23 2023-05-02 成都理工大学 基于振幅统计和时频分析的深度震相快速拾取方法

Family Cites Families (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104216008B (zh) * 2013-06-05 2017-02-08 中国石油天然气集团公司 一种井中压裂微地震事件识别方法
CA2872289A1 (en) * 2013-11-25 2015-05-25 King Abdullah University Of Science And Technology High repetition rate thermometry system and method
CN104914468B (zh) * 2015-06-09 2017-06-06 中南大学 一种矿山微震信号p波初至时刻联合拾取方法
CN106886044B (zh) * 2017-03-02 2019-02-12 吉林大学 一种基于剪切波与Akaike信息准则的微地震初至拾取方法
CN106896407B (zh) * 2017-03-28 2018-07-13 吉林大学 一种基于近似负熵的微地震信号初至拾取方法
CN107132576B (zh) * 2017-07-05 2018-10-30 西安交通大学 基于二阶同步挤压小波变换的地震资料时频分析方法
CN107884822B (zh) * 2017-11-13 2019-09-27 北京矿冶研究总院 一种提高矿用微震震源定位精度的方法
CN110398775B (zh) * 2019-08-23 2021-04-06 山东大学 隧道突涌水灾害微震事件信号波动初至拾取方法及系统

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115016001A (zh) * 2022-08-09 2022-09-06 北京科技大学 一种基于维度变换的微震信号到时自动拾取方法
CN115016001B (zh) * 2022-08-09 2023-05-16 北京科技大学 一种基于维度变换的微震信号到时自动拾取方法
CN116047604A (zh) * 2023-02-23 2023-05-02 成都理工大学 基于振幅统计和时频分析的深度震相快速拾取方法
CN116047604B (zh) * 2023-02-23 2023-10-17 成都理工大学 基于振幅统计和时频分析的深度震相快速拾取方法

Also Published As

Publication number Publication date
CN112904412B (zh) 2024-02-23

Similar Documents

Publication Publication Date Title
Jian et al. On the denoising method of prestack seismic data in wavelet domain
US20140297188A1 (en) Time-frequency representations of seismic traces using wigner-ville distributions
US11892578B2 (en) Seismic imaging method, system, and device based on pre-stack high-angle fast Fourier transform
CN112904412A (zh) 一种矿山微震信号p波初至时刻提取方法及系统
CN106443787A (zh) 叠前地震道集压噪方法及其装置
CN109061724A (zh) 一种基于自适应变分模态分解的地震数据降噪方法
Wang et al. Relational database for horizontal‐to‐vertical spectral ratios
Zhou et al. An improved joint method for onset picking of acoustic emission signals with noise
Lyubushin Global seismic noise entropy
Chi-Durán et al. Automatic detection of P-and S-wave arrival times: new strategies based on the modified fractal method and basic matching pursuit
CN115270076A (zh) 一种品质因子确定方法、装置、电子设备及存储介质
JPWO2019049263A1 (ja) シンプルプロファイリング(spm)手法の変換プログラムおよびシンプルプロファイリング(spm)手法の変換方法
Bueno et al. Recursive entropy method of segmentation for seismic signals
Li et al. A Reliable Strategy for Improving Automatic First‐Arrival Picking of High‐Noise Three‐Component Microseismic Data
Pardo-Igúzquiza et al. CYSTRATI: a computer program for spectral analysis of stratigraphic successions
Liu et al. A new strategy of finite‐fault inversion using multiscale waveforms and its application to the 2015 Gorkha, Nepal, earthquake
CN109633744B (zh) 地震子波的提取方法、装置、设备及存储介质
CN116804772A (zh) 一种利用mp-vmd算法去除微地震噪声信号的算法
Yang et al. New baseline correction method for near‐fault ground‐motion records based on continuous wavelet transform
CN113568044A (zh) 阵列声波测井首波到时确定方法和装置
Atefi et al. Rapid estimation of earthquake magnitude by a new wavelet‐based proxy
Li et al. Noise reduction method of microseismic signal of water inrush in tunnel based on variational mode method
CN112782763B (zh) 一种地震品质因子估算方法、装置、设备及存储介质
Yin et al. Research on interference signal recognition in p wave pickup and magnitude estimation
Wodecki et al. Combination of Principal Component Analysis and Time‐Frequency Representation for P‐Wave Arrival Detection

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant