CN105607121B - 一种煤炭陷落柱识别方法和装置 - Google Patents

一种煤炭陷落柱识别方法和装置 Download PDF

Info

Publication number
CN105607121B
CN105607121B CN201610073638.7A CN201610073638A CN105607121B CN 105607121 B CN105607121 B CN 105607121B CN 201610073638 A CN201610073638 A CN 201610073638A CN 105607121 B CN105607121 B CN 105607121B
Authority
CN
China
Prior art keywords
amplitude value
seismic
data
imaging
shot data
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.)
Active
Application number
CN201610073638.7A
Other languages
English (en)
Other versions
CN105607121A (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.)
China University of Mining and Technology Beijing CUMTB
Original Assignee
China University of Mining and Technology Beijing CUMTB
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 China University of Mining and Technology Beijing CUMTB filed Critical China University of Mining and Technology Beijing CUMTB
Priority to CN201610073638.7A priority Critical patent/CN105607121B/zh
Publication of CN105607121A publication Critical patent/CN105607121A/zh
Application granted granted Critical
Publication of CN105607121B publication Critical patent/CN105607121B/zh
Priority to US15/758,559 priority patent/US10302788B2/en
Priority to PCT/CN2016/113064 priority patent/WO2017133361A1/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
    • G01V1/30—Analysis
    • G01V1/307—Analysis for determining seismic attributes, e.g. amplitude, instantaneous phase or frequency, reflection strength or polarity
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
    • G01V1/30—Analysis
    • G01V1/301—Analysis for determining seismic cross-sections or geostructures
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
    • G01V1/36—Effecting static or dynamic corrections on records, e.g. correcting spread; Correlating seismic signals; Eliminating effects of unwanted energy
    • G01V1/362—Effecting static or dynamic corrections; Stacking
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
    • G01V1/30—Analysis
    • G01V1/303—Analysis for determining velocity profiles or travel times
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00—Details of seismic processing or analysis
    • G01V2210/50—Corrections or adjustments related to wave propagation
    • G01V2210/51—Migration
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00—Details of seismic processing or analysis
    • G01V2210/60—Analysis
    • G01V2210/62—Physical property of subsurface
    • G01V2210/622—Velocity, density or impedance
    • G01V2210/6222—Velocity; travel time
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00—Details of seismic processing or analysis
    • G01V2210/60—Analysis
    • G01V2210/64—Geostructures, e.g. in 3D data cubes
    • G01V2210/641—Continuity of geobodies
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00—Details of seismic processing or analysis
    • G01V2210/60—Analysis
    • G01V2210/64—Geostructures, e.g. in 3D data cubes
    • G01V2210/642—Faults
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00—Details of seismic processing or analysis
    • G01V2210/60—Analysis
    • G01V2210/64—Geostructures, e.g. in 3D data cubes
    • G01V2210/646—Fractures

Landscapes

  • Engineering & Computer Science (AREA)
  • Remote Sensing (AREA)
  • Physics & Mathematics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Acoustics & Sound (AREA)
  • Environmental & Geological Engineering (AREA)
  • Geology (AREA)
  • General Life Sciences & Earth Sciences (AREA)
  • General Physics & Mathematics (AREA)
  • Geophysics (AREA)
  • Geophysics And Detection Of Objects (AREA)

Abstract

本发明提供了一种煤炭陷落柱识别方法和装置，包括：获取目标区域的地震炮集数据和地震波偏移速度文件；根据上述数据计算每个单炮数据在不同成像点的绕射波走时；对每个单炮数据及其绕射波走时进行马氏距离计算处理，获取每个单炮数据的绕射波；分别对每个单炮数据的绕射波进行成像处理；将地震炮集数据对应的所有单炮数据的成像处理结果进行叠加处理，得到地震炮集数据的绕射波成像结果；根据该绕射波成像结果进行煤炭陷落柱识别；其通过马氏距离提取地震炮集数据对应的绕射波，由于绕射波携带的是更小尺度的地质信息，故通过成像后的绕射波能够准确的识别陷落柱地质体，从而能够降低煤炭开采中由陷落柱而引发的突水、瓦斯泄露等事故的发生风险。

Description

一种煤炭陷落柱识别方法和装置
技术领域
本发明涉及地震数据成像领域，具体而言，涉及一种基于地震绕射波成像的煤炭陷落柱识别方法和装置。
背景技术
在煤炭开采过程中，断层、断裂、陷落柱等不连续地质体会破坏煤层的连续性，易于诱发透水、瓦斯突出等事故，严重威胁着煤矿作业的安全，因此，有效识别断层、断裂和陷落柱等不连续地质体具有重要的意义。
目前，用于识别断层、断裂的方法有很多种方法，包括地震相干体技术(Bahorichand Farmer,1995；Marfurt,et al.,1998,1999；Gersztenkorn and Marfurt,1999),匹配追踪算法(Mallat and Zhang,1993；Castagna et al.,2003；Liu and Marfurt,2005)和谱分解算法(Partyka et al.,1999；Puryear et al.,2012；Gao et al.,2013)等。但是，由于陷落柱本身的大小不等、分布规律差等特点，而使得上述的几种方法均不能有效识别陷落柱，因此，陷落柱识别一直都是地震勘探的难题。
为了解决上述难题，相关技术采取地震波成像技术对陷落柱进行识别，具体方法包括，采集地震炮集数据，通过现有地震波成像技术针对地震炮集数据中的反射波进行成像处理，通过对反射波的成像处理结果进行煤炭陷落柱识别，但是，反射波为宏观尺度地质单元响应，通常只能解决空间展布大于一个地震波长的地质体勘探问题，其对于煤炭陷落柱等小尺度的地质信息的识别效果较差。因此针对反射波研发的地震波成像技术针对陷落柱等小尺度的地质信息成像上效果不佳，不能很好的对煤炭陷落柱等小尺度的地质信息进行识别。
发明内容
本发明的目的在于提供一种煤炭陷落柱识别方法和装置，能够提取获取的地震炮集数据中携带的是更小尺度的地质信息的绕射波，通过成像后的绕射波准确的评价陷落柱地质体，从而准确的进行目标区域中煤炭陷落柱的识别。
第一方面，本发明实施例提供了一种煤炭陷落柱识别方法，包括：
获取目标区域的地震炮集数据和地震炮集数据对应的地震波偏移速度文件；其中，地震炮集数据包括多个单炮数据；每个单炮数据均包括地下岩层界面反射或者折射的多个地震波；地震波携带有地质信息；
根据每个单炮数据对应的检波器空间位置和地震波偏移速度文件，计算由每个单炮数据的地震波炮点位置经地下成像空间任意一个成像点位置到每个单炮数据的地震波检波点位置的每个单炮数据的绕射波走时；
对每个单炮数据及单炮数据的绕射波走时进行马氏距离计算处理，获取每个单炮数据的绕射波振幅值样点；其中，绕射波振幅值样点携带有小尺度地质信息；小尺度地质信息至少包括：地层层位信息、断层信息和陷落柱信息；
分别对每个单炮数据的绕射波振幅值样点进行成像处理，得到每个单炮数据的绕射波成像结果；
将地震炮集数据对应的所有单炮数据的绕射波成像结果进行叠加处理，得到地震炮集数据的绕射波成像结果，以便根据地震炮集数据的绕射波成像结果进行煤炭陷落柱识别。
结合第一方面，本发明实施例提供了第一方面的第一种可能的实施方式，其中，根据每个单炮数据对应的检波器空间位置和地震波偏移速度文件，计算由每个单炮数据的地震波炮点位置经地下成像空间任意一个成像点位置到每个单炮数据的地震波检波点位置的每个单炮数据的绕射波走时包括：
根据每个单炮数据对应的检波器空间位置和地震波偏移速度文件，计算由每个单炮数据的地震波炮点位置到每一个地下成像空间成像点位置的走时tS；
根据每个单炮数据对应的检波器空间位置和地震波偏移速度文件，计算由每一个地下成像空间成像点位置到每个单炮数据的地震波检波点位置的走时tR；其中，每个单炮数据在一个成像点的走时tR与走时tS相对应；
分别对每个单炮数据的走时tS以及与走时tS对应的走时tR进行求和处理，得到对应于每个单炮数据的多个绕射波走时。
结合第一方面的第一种可能的实施方式，本发明实施例提供了第一方面的第二种可能的实施方式，其中，对每个单炮数据及单炮数据的绕射波走时进行马氏距离计算处理，获取每个单炮数据的绕射波振幅值样点包括：
对于任意一个单炮数据，根据单炮数据和单炮数据的多个绕射波走时，计算单炮数据对应的地下成像空间中每个成像点的振幅值样点序列，得到对应于单炮数据的多个振幅值样点序列；
根据振幅值样点的能量大小，分别对得到的每一个振幅值样点序列中的振幅值样点进行排序处理；
分别对排序处理后的每一个振幅值样点序列进行马氏距离计算处理，得到每一个振幅值样点序列中按预设顺序排列的所有振幅值样点的马氏距离；
根据每个振幅值样点序列中按预设顺序排列的所有振幅值样点的马氏距离，提取对应的每个振幅值样点序列中符合预设条件的振幅值样点，得到对应于每个单炮数据的多组绕射波振幅值样点。
结合第一方面的第二种可能的实施方式，本发明实施例提供了第一方面的第三种可能的实施方式，其中，分别对排序处理后的每一个振幅值样点序列进行马氏距离计算处理，得到每一个振幅值样点序列中按预设顺序排列的所有振幅值样点的马氏距离包括：
通过以下公式计算每一个按预设顺序排列的振幅值样点序列的中值其中，μx为振幅值样点序列的中值，median为中值运算，为按预设顺序排序后的振幅样点序列，N为振幅值样点序列中样点个数；
通过以下公式计算每一个按预设顺序排列的振幅值样点序列的中值绝对值偏差平方其中，为中值绝对值偏差；
根据中值和中值绝对值偏差平方，计算每一个振幅值样点序列中按预设顺序排列的所有振幅值样点的马氏距离其中，为马氏距离。
结合第一方面的第三种可能的实施方式，本发明实施例提供了第一方面的第四种可能的实施方式，其中，根据每个振幅值样点序列中按预设顺序排列的所有振幅值样点的马氏距离，提取对应的每个振幅值样点序列中符合预设条件的振幅值样点包括：
对于任意一个振幅值样点序列，根据振幅值样点序列中按预设顺序排列的每个振幅值样点的马氏距离以及反射波振幅值样点的预设马氏距离阈值，计算振幅值样点序列中的样点去除数量；
根据预设马氏距离阈值和样点去除数量，对振幅值样点序列中按预设顺序排列的振幅值样点进行样点去除处理，得到振幅值样点序列中剩余的振幅值样点；
提取剩余的振幅值样点作为振幅值样点序列中符合预设条件的振幅值样点。
结合第一方面的第四种可能的实施方式，本发明实施例提供了第一方面的第五种可能的实施方式，其中，分别对每个单炮数据的绕射波振幅值样点进行成像处理，得到每个单炮数据的绕射波成像结果包括：
将在每个振幅值样点序列中提取的符合预设条件的振幅值样点分别进行求和处理，得到每个振幅值样点序列对应的绕射波成像结果；
将单炮数据对应的所有振幅值样点序列的绕射波成像结果进行叠加处理，得到对应于单炮数据的绕射波成像结果。
结合第一方面的第五种可能的实施方式，本发明实施例提供了第一方面的第六种可能的实施方式，其中，将在每个振幅值样点序列中提取的符合预设条件的振幅值样点分别进行求和处理，得到每个振幅值样点序列对应的绕射波成像结果包括：
通过以下公式将在每个振幅值样点序列中提取的符合预设条件的振幅值样点分别进行求和处理得到每个振幅值样点序列对应的绕射波成像结果；其中，V(m)为绕射波成像结果，m＝m(x,y,z)为地下成像空间每个成像点位置，u(r,tS+tR)为预处理后的地震炮集数据，r(x,y,z)为检波点位置，tS,tR分别为由炮点位置到成像点位置和由成像点位置到检波器位置的走时，A为偏移成像孔径，w(m,r)为几何扩散因子，TriM为基于马氏距离和反射波振幅去除的绕射波成像算子。
结合第一方面，本发明实施例提供了第一方面的第七种可能的实施方式，其中，获取目标区域的地震炮集数据和地震炮集数据对应的地震波偏移速度文件包括：
获取目标区域的地震炮集数据；
对地震炮集数据进行地震预处理，得到可用于偏移成像的地震炮集数据；其中，地震预处理至少包括：去除噪声处理和静校正处理；
对地震炮集数据进行偏移速度分析处理，得到地震炮集数据对应的地震波偏移速度文件。
第二方面，本发明实施例还提供了一种煤炭陷落柱识别装置，装置包括：
获取模块，用于获取目标区域的地震炮集数据和地震炮集数据对应的地震波偏移速度文件；其中，地震炮集数据包括多个单炮数据；每个单炮数据均包括地下岩层界面反射或者折射的多个地震波；地震波携带有地质信息；
绕射波走时计算模块，用于根据每个单炮数据对应的检波器空间位置和地震波偏移速度文件，计算由每个单炮数据的地震波炮点位置经地下成像空间任意一个成像点位置到每个单炮数据的地震波检波点位置的每个单炮数据的绕射波走时；
马氏距离计算处理模块，用于对每个单炮数据及单炮数据的绕射波走时进行马氏距离计算处理，获取每个单炮数据的绕射波振幅值样点；其中，绕射波振幅值样点携带有小尺度地质信息；小尺度地质信息至少包括：地层层位信息、断层信息和陷落柱信息；
成像处理模块，用于分别对每个单炮数据的绕射波振幅值样点进行成像处理，得到每个单炮数据的绕射波成像结果；
叠加处理模块，用于将地震炮集数据对应的所有单炮数据的绕射波成像结果进行叠加处理，得到地震炮集数据的绕射波成像结果，以便根据地震炮集数据的绕射波成像结果进行煤炭陷落柱识别。
结合第二方面，本发明实施例提供了第二方面的第一种可能的实施方式，其中，绕射波走时计算模块包括：
第一走时计算单元，用于根据每个单炮数据对应的检波器空间位置和地震波偏移速度文件，计算由每个单炮数据的地震波炮点位置到每一个地下成像空间成像点位置的走时tS；
第二走时计算单元，用于根据每个单炮数据对应的检波器空间位置和地震波偏移速度文件，计算由每一个地下成像空间成像点位置到每个单炮数据的地震波检波点位置的走时tR；其中，每个单炮数据在一个成像点的走时tR与走时tS相对应；
第一求和处理单元，用于分别对每个单炮数据的走时tS以及与走时tS对应的走时tR进行求和处理，得到对应于每个单炮数据的多个绕射波走时。
本发明实施例提供的一种煤炭陷落柱识别方法和装置，包括：获取目标区域的地震炮集数据和地震波偏移速度文件；根据每个单炮数据对应的检波器空间位置和地震波偏移速度文件，计算每个单炮数据对应在不同成像点的绕射波走时；对每个单炮数据及单炮数据的绕射波走时进行马氏距离计算处理，获取每个单炮数据的绕射波；分别对每个单炮数据的绕射波进行成像处理，得到每个单炮数据的绕射波成像结果；将地震炮集数据对应的所有单炮数据的绕射波成像结果进行叠加处理，得到地震炮集数据的绕射波成像结果，并根据地震炮集数据的绕射波成像结果进行煤炭陷落柱识别，
与现有技术中的针对反射波研发的地震波成像技术进行陷落柱识别的效果不佳相比，其首先计算地震炮集数据中每个单炮数据对应的绕射波走时，并通过马氏距离对单炮数据及单炮数据的绕射波走时进行计算处理，获取所有单炮数据对应的绕射波，并对绕射波进行成像处理；由于绕射波携带的是更小尺度的地质信息，故通过成像后的绕射波能够准确的评价陷落柱地质体，即能够准确的进行目标区域中煤炭陷落柱的识别，从而能够降低煤炭开采中由陷落柱而引发的突水、瓦斯泄露等事故发生风险，减少不必要的人员伤亡及经济损失。
为使本发明的上述目的、特征和优点能更明显易懂，下文特举较佳实施例，并配合所附附图，作详细说明如下。
附图说明
为了更清楚地说明本发明实施例的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，应当理解，以下附图仅示出了本发明的某些实施例，因此不应被看作是对范围的限定，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他相关的附图。
图1示出了本发明实施例所提供的一种煤炭陷落柱识别方法的流程图；
图2示出了本发明实施例所提供的另一种煤炭陷落柱识别方法的流程图；
图3示出了本发明实施例所提供的另一种煤炭陷落柱识别方法的流程图；
图4示出了本发明实施例所提供的另一种煤炭陷落柱识别方法的流程图；
图5示出了本发明实施例所提供的另一种煤炭陷落柱识别方法的流程图；
图6示出了本发明实施例所提供的常规反射波成像结果的示意图；
图7示出了应用本发明实施例所提供的一种煤炭陷落柱识别方法得到的绕射波成像结果的示意图；其中，包括了相应位置的剖面和切片图；
图8示出了本发明实施例所提供的一种煤炭陷落柱识别装置的结构示意图；
图9示出了本发明实施例所提供的一种煤炭陷落柱识别装置中绕射波走时计算模块和马氏距离计算处理模块的结构示意图。
具体实施方式
下面将结合本发明实施例中附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。通常在此处附图中描述和示出的本发明实施例的组件可以以各种不同的配置来布置和设计。因此，以下对在附图中提供的本发明的实施例的详细描述并非旨在限制要求保护的本发明的范围，而是仅仅表示本发明的选定实施例。基于本发明的实施例，本领域技术人员在没有做出创造性劳动的前提下所获得的所有其他实施例，都属于本发明保护的范围。
目前采取地震波成像技术识别陷落柱的方法主要依赖于反射波的成像结果，从而不能很好的对煤炭陷落柱等小尺度的地质信息进行识别。实际中，绕射波携带的是小尺度的地质信息，因此可以探测到空间展布小于一个地震波长的地质信息；随着地震波成像技术的完善和计算机处理能力的发展，目前已发展到利用绕射波直接成像，以探测小尺度地质体，如岩溶洞穴、裂缝等。在绕射波成像技术研究上，主要分为三大类，包括基于信号的绕射波分离及成像方法(Harlan,et al.,1984；Taner et al.,2006；Fomel et al.,2006,2007)；Bansal and Inhof,2005)、基于聚焦思想的绕射波成像方法(Berkovitch et al.,2009；Dell and Gajewski,2011；Asgedom et al.,2011)和修改Kirchhoff成像函数的绕射波成像方法(Zhang,2004；Moser and Howard,2008；Figueiredo et al.,2013；Zhao etal.,2015)。
本发明实施例提供的一种煤炭陷落柱识别方法和装置(也可以称为基于地震绕射波成像的煤炭陷落柱识别方法和装置)，属于上述第三类绕射波成像技术，原理上依据绕射波与反射波振幅特征差异，通过统计学马氏距离准则实现绕射波成像，由于绕射波携带的是更小尺度的地质信息，故直接应用绕射波成像评价煤炭陷落柱，能够更加真实反映地下地质情况，从而准确的进行目标区域中煤炭陷落柱的识别，且无需事先去除反射波，具有高效计算能力。
参考图1，本发明实施例提供了一种煤炭陷落柱识别方法，所述方法包括：
S101、获取目标区域的地震炮集数据和所述地震炮集数据对应的地震波偏移速度文件；其中，地震炮集数据包括多个单炮数据；每个单炮数据均包括地下岩层界面反射或者折射的多个地震波；地震波携带有地质信息。
具体的，目标区域为预先选定的待进行煤炭陷落柱识别的区域，该目标区域包括多个震源和多个检波器，一个震源对应多个检波器；其中，多个震源对应多个地震波炮点位置，多个检波器对应多个地震波检波点位置；其中，多个震源可以发射多个单炮数据，多个单炮数据组成地震炮集数据；对应的，多个检波器能够采集多个单炮数据，多个单炮数据组成地震炮集数据。
具体的，上述地震波偏移速度文件是对上述地震炮集数据进行偏移速度分析处理得到的；其中，上述地震波偏移速度文件，对应于地层速度参数，用于计算地震波传播旅行时(即时间)，其是由对地震炮集数据进行偏移速度分析获得；而上述偏移速度分析是一种速度建模技术，该速度建模技术通过分析地震采集数据中的地震波聚焦性确定速度参数，即获得地震波偏移速度文件。
其中，上述地震波中携带的地质信息包括小尺度地质信息，还包括大尺度地质信息；其中，上述小尺度地址信息包括：地层层位信息、断层信息和陷落柱信息。
S102、根据每个所述单炮数据对应的检波器空间位置和所述地震波偏移速度文件，计算由每个单炮数据的地震波炮点位置经地下成像空间任意一个成像点位置到每个单炮数据的地震波检波点位置的每个单炮数据的绕射波走时；其中，每一个所述单炮数据对应的检波器为多个，且每个所述检波器均设置在对应的所述检波点位置上。
具体的，根据每个单炮数据对应的检波器空间位置和该单炮数据对应的地震波偏移速度文件进行绕射波走时计算，即根据射线追踪程函方程进行绕射波走时计算，而射线追踪程函方程是速度(即单炮数据对应的地震波偏移速度文件)，空间位置(即单炮数据对应的检波器空间位置)和地震波走时的微分关系。
由于获取的目标区域的地震炮集数据包括多个单炮数据，那么对于每一个单炮数据，根据该单炮数据和对应的地震波偏移速度文件，首先计算该单炮数据的地震波炮点位置到地下成像空间任意一个成像点位置对应的走时(为了区分下面的走时，此处成为第一走时)，由于地下成像空间中包括多个成像点，故计算的该单炮数据的地震波炮点位置到地下成像空间任意一个成像点位置对应的第一走时也为多个；
然后再计算地下成像空间任意一个成像点位置到该单炮数据的地震波检波点位置的多个绕射波走时(此处称为第二走时)，其中，每一个成像点处的一个第二走时均与一个第一走时相对应。
最后将上述第一走时和第一走时对应的第二走时相加即可得到一个单炮数据对应的绕射波走时。
对于地震炮集数据中的每个单炮数据均可以按照上述单炮数据的绕射波走时方式进行计算，得到地震炮集数据中的每个单炮数据的绕射波走时。
S103、对每个所述单炮数据及所述单炮数据的绕射波走时进行马氏距离计算处理，获取每个所述单炮数据的绕射波振幅值样点；其中，所述绕射波振幅值样点携带有小尺度地质信息；所述小尺度地质信息至少包括：地层层位信息、断层信息和陷落柱信息。
具体的，对于每一个单炮数据，根据该单炮数据及该单炮数据的绕射波走时可以计算得出该单炮数据对应的地下成像空间每个成像点的多个振幅值样点序列(即多个成像振幅值样点序列)；
对于每一个振幅值样点序列，计算该振幅值样点序列中每个振幅值样点的马氏距离，最后根据计算的马氏距离剔除该振幅值样点序列中的反射波对应振幅值样点，得到的剩余的振幅值样点即为该振幅值样点序列中对应的绕射波振幅值样点。
其中，一个单炮数据包括多个振幅值样点序列，故根据上述一个振幅值样点序列中绕射波振幅值样点的计算方法即可以计算出一个单炮数据对应的多组振幅值样点序列。
按照上述单炮数据对应的多组振幅值样点序列的计算方法，可以计算地震炮集数据中所有单炮数据对应的振幅值样点序列。
S104、分别对每个所述单炮数据的绕射波振幅值样点进行成像处理，得到每个所述单炮数据的绕射波成像结果。
具体的，对于每个单炮数据，每个单炮数据包括多个振幅值样点序列；将每一个振幅值样点序列中的所有的绕射波振幅值样点进行求和，即可得到每一个振幅值样点序列对应的绕射波成像结果；
对于每个单炮数据，将其包括的每一个振幅值样点序列对应的绕射波成像结果进行叠加处理，即可得到该单炮数据对应的绕射波成像结果。
S105、将所述地震炮集数据对应的所有单炮数据的绕射波成像结果进行叠加处理，得到所述地震炮集数据的绕射波成像结果，以便根据所述地震炮集数据的绕射波成像结果进行煤炭陷落柱识别。
由于获取的地震炮集数据才是携带反应目标区域的地质信息的数据，而地震炮集数据中的每个单炮数据只是携带反应目标区域的地质信息的一部分数据；在实际计算过程中是针对每个单炮数据单独进行成像处理，最后将每一个单炮数据的成像处理结果进行叠加处理，即可得到地震炮集数据的绕射波成像结果。
根据最终得到的地震炮集数据的绕射波成像结果进行煤炭陷落柱识别，即可准确的评价陷落柱地质体，从而准确的进行目标区域中煤炭陷落柱的识别。
本发明实施例提供的一种煤炭陷落柱识别方法，与现有技术中的针对反射波研发的地震波成像技术进行陷落柱识别的效果不佳相比，其首先计算地震炮集数据中每个单炮数据对应的绕射波走时，并通过马氏距离对单炮数据及单炮数据的绕射波走时进行计算处理，获取所有单炮数据对应的绕射波，并对绕射波进行成像处理；由于绕射波携带的是更小尺度的地质信息，故通过成像后的绕射波能够准确的评价陷落柱地质体，即能够准确的进行目标区域中煤炭陷落柱的识别，从而能够降低煤炭开采中由陷落柱而引发的突水、瓦斯泄露等事故发生风险，减少不必要的人员伤亡及经济损失。
具体的，由于获取的目标区域的地震炮集数据包括多个单炮数据，本发明实施例中技术每一个单炮数据对应的绕射波走时，其中，每一个单炮数据对应的地下成像空间包括多个成像点，故每一个单炮数据对应的绕射波走时也为多个，参考图2，计算每个单炮数据对应的多个绕射波走时的具体方法包括如下步骤：
S201、根据每个所述单炮数据对应的检波器空间位置和所述地震波偏移速度文件，计算由每个单炮数据的地震波炮点位置到每一个地下成像空间成像点位置的走时tS。
具体的，由于每个单炮数据对应的地下成像空间成像点有多个，故对应的地下成像空间成像点位置也有多个，故计算该单炮数据到每一个地下成像空间成像点位置的走时tS也对应有多个；即不同的地下成像空间成像点位置对应不同的走时tS。
S202、根据每个所述单炮数据对应的检波器空间位置和所述地震波偏移速度文件，计算由每一个地下成像空间成像点位置到每个单炮数据的地震波检波点位置的走时tR；其中，每个所述单炮数据在一个成像点的走时tR与走时tS相对应。
同样的，由于每个单炮数据对应的地下成像空间成像点有多个，故对应的地下成像空间成像点位置也有多个，故每一个地下成像空间成像点位置到每个单炮数据的地震波检波点位置的走时tR有对应有多个，即不同的地下成像空间成像点位置对应不同的走时tR。
需要说明的是，每个所述单炮数据在一个成像点的走时tR与走时tS相对应。
S203、分别对每个单炮数据的走时tS以及与所述走时tS对应的走时tR进行求和处理，得到对应于每个所述单炮数据的多个绕射波走时。
参考图3，本发明实施例中，上述步骤103中对每个所述单炮数据及所述单炮数据的绕射波走时进行马氏距离计算处理，获取每个所述单炮数据的绕射波振幅值样点的具体步骤包括：
S301、对于任意一个单炮数据，根据所述单炮数据和所述单炮数据的多个绕射波走时，计算所述单炮数据对应的地下成像空间中每个成像点的振幅值样点序列，得到对应于所述单炮数据的多个振幅值样点序列。
具体的，每一个单炮数据对应于不同的地下成像空间中的成像点有不同的绕射波走时；对于任意一个单炮数据，根据该单炮数据和该单炮数据在一个成像点处的一个绕射波走时，可以计算处该单炮数据在该成像点的一个振幅值样点序列，而对于该单炮数据的每一个成像点，均能计算出对应的振幅值样点序列，而一个单炮数据又对应多个成像点，故可以计算得到对应于所述单炮数据的多个振幅值样点序列。
S302、根据振幅值样点的能量大小，分别对得到的每一个振幅值样点序列中的振幅值样点进行排序处理。
其中，对应于每个单炮数据，每个单炮数据包括多个振幅值样点序列，每个振幅值样点序列中包括多个振幅值样点；由于每个振幅值样点对应的能量大小不同，故根据振幅值样点的能量大小，可以分别对得到的每一个振幅值样点序列中的振幅值样点进行排序处理，使得每一个振幅值样点序列中的振幅值样点均按预设顺序进行排列。
S303、分别对排序处理后的每一个振幅值样点序列进行马氏距离计算处理，得到每一个振幅值样点序列中按预设顺序排列的所有振幅值样点的马氏距离。
具体的，每个单炮数据均包括多个振幅值样点序列，对于每一个振幅值样点序列，计算该振幅值样点序列中的按照预设顺序进行排列的每一个振幅值样点的马氏距离。
对应的，将每个单炮数据包括多个振幅值样点序列均按照上述马氏距离计算方法进行计算处理，即可得到每一个振幅值样点序列中按预设顺序排列的每个振幅值样点的马氏距离。
S304、根据每个振幅值样点序列中按预设顺序排列的所有振幅值样点的马氏距离，提取对应的每个振幅值样点序列中符合预设条件的振幅值样点，得到对应于每个所述单炮数据的多组绕射波振幅值样点。
具体的，绕射波的振幅值样点的马氏距离和反射波的振幅值样点的马氏距离均是有预设阈值的，该预设阈值是根据多次实验计算得到，因此，可以根据绕射波的振幅值样点的第一预设马氏距离，提取对应的每个振幅值样点序列中符合该第一预设马氏距离的振幅值样点，将提取的振幅值样点作为绕射波振幅值样点；也可以根据反射波的振幅值样点的第二预设马氏距离，剔除对应的每个振幅值样点序列中符合该第二预设马氏距离的振幅值样点，得到的剩余的振幅值样点即为符合预设条件的振幅值样点，同样将提取的振幅值样点作为绕射波振幅值样点。
由于每个单炮数据对应多个振幅值样点序列，故每个单炮数据也对应多组绕射波振幅值样点。即每个单炮数据的绕射波振幅值样点的组数与该单炮数据的振幅值样点序列的个数相同。
本发明实施例中，上述步骤303中计算振幅值样点的马氏距离的方法具体如下：
通过以下公式计算每一个按预设顺序排列的振幅值样点序列的中值其中，μx为振幅值样点序列的中值，median为中值运算，为按预设顺序排序后的振幅样点序列，N为振幅值样点序列中样点个数；
通过以下公式计算每一个按预设顺序排列的所述振幅值样点序列的中值绝对值偏差平方其中，为中值绝对值偏差；
根据所述中值和所述中值绝对值偏差平方，计算每一个振幅值样点序列中按预设顺序排列的所有振幅值样点的马氏距离 d ( x ( i ) o , μ x ) = [ x ( i ) o - μ x ] 2 σ x 2 ; 其中，为马氏距离。
参考图4，本发明实施例中的上述步骤304中，计算每个单炮数据的多组绕射波振幅值样点的具体方法包括：
S401、对于任意一个振幅值样点序列，根据所述振幅值样点序列中按预设顺序排列的每个振幅值样点的马氏距离以及反射波振幅值样点的预设马氏距离阈值，计算所述振幅值样点序列中的样点去除数量。
具体的，对于任意一个振幅值样点序列，将计算的马氏距离与反射波振幅值样点的预设马氏距离阈值(即上述第二预设马氏距离)进行对比，将符合反射波振幅值样点的预设马氏距离阈值的计算结果作为剔除处理对象，按照该种方式即可以计算出振幅值样点序列中的样点去除数量。
S402、根据所述预设马氏距离阈值和所述样点去除数量，对所述振幅值样点序列中按预设顺序排列的振幅值样点进行样点去除处理，得到所述振幅值样点序列中剩余的振幅值样点。
具体的，根据样点去除数量在排序后的振幅值样点序列中提取绕射波振幅值样点，即对排序后的振幅值样点序列，去除序列首段和末端对应的不符合条件的振幅值样点，去除个数则根据上述计算的样点去除数量决定。
S403、提取所述剩余的振幅值样点作为所述振幅值样点序列中符合预设条件的振幅值样点。
具体的，剩余的振幅值样点都是满足绕射波的振幅值样点的第一预设马氏距离，故提取这些剩余的振幅值样点作为符合预设条件的振幅值样点。
本发明实施例中，步骤104中分别对每个所述单炮数据的绕射波振幅值样点进行成像处理，得到每个所述单炮数据的绕射波成像结果的具体方法包括：
将在每个振幅值样点序列中提取的符合预设条件的振幅值样点分别进行求和处理，得到每个振幅值样点序列对应的绕射波成像结果；
将所述单炮数据对应的所有振幅值样点序列的绕射波成像结果进行叠加处理，得到对应于所述单炮数据的绕射波成像结果。
具体的，由于每个单炮数据包括多个振幅值样点序列，针对每个振幅值样点序列中提取出了预设条件的振幅值样点(即绕射波的振幅值样点)，然后将每个振幅值样点序列的绕射波的振幅值样点进行求和处理，即可得到该振幅值样点序列对应的绕射波成像结果；
在求出每个单炮数据中的所有的该振幅值样点序列对应的绕射波成像结果之后，将这些绕射波成像结果进行叠加处理，即可得到对应于该单炮数据的绕射波成像结果。
对于地震炮集数据，则将其包括的每个单炮数据的绕射波成像结果同样进行叠加处理，即可得到地震炮集数据的绕射波成像结果。
具体的，上述计算每个振幅值样点序列对应的绕射波成像结果包括：
通过以下公式将在每个振幅值样点序列中提取的符合预设条件的振幅值样点分别进行求和处理得到每个振幅值样点序列对应的绕射波成像结果；其中，V(m)为绕射波成像结果，m＝m(x,y,z)为地下成像空间每个成像点位置，u(r,tS+tR)为预处理后的地震炮集数据，r(x,y,z)为检波点位置，tS,tR分别为由炮点位置到成像点位置和由成像点位置到检波器位置的走时，A为偏移成像孔径，w(m,r)为几何扩散因子，TriM为基于马氏距离和反射波振幅去除的绕射波成像算子。
本发明实施例中，为了保证更好在获取的地震炮集数据中提取绕射波，还需要对获取的地震炮集数据进行处理，参考图5，具体方法包括：
S501、获取目标区域的地震炮集数据；
S502、对所述地震炮集数据进行地震预处理，得到可用于偏移成像的地震炮集数据；其中，所述地震预处理至少包括：去除噪声处理和静校正处理。
S503、对所述地震炮集数据进行偏移速度分析处理，得到所述地震炮集数据对应的地震波偏移速度文件。
结合上述步骤501-步骤503，由检波器接收炮集数据(即地震波)，加载观测系统，将获取的炮集数据发送至观测系统，由观测系统对地震炮集数据进行地震预处理，即观测系统加载接收的地震炮集数据，并对加载的地震波进行噪声去除、速度分析和偏移后获得地震炮集数据和地震炮集数据对应的地震波偏移速度文件。
下面结合具体实施例对本发明实施例提供的一种煤炭陷落柱识别方法进行说明
通过煤炭三维地震资料，说明一种基于地震绕射波成像的煤炭陷落柱评价技术及装置在煤炭陷落柱评价中的应用效果。
(1)读入地震炮集数据和偏移速度文件；
(2)由输入的偏移速度模型，计算地震波走时表；
(3)根据走时表和地震炮集数据，得出绕射波成像结果；
(4)为对比绕射波成像结果和常规反射波成像结果在陷落柱评价上的应用效果，本实施例提供了两者对比图，如图6为常规反射波成像结果，包括了主测线方向和联络线方向剖面以及0.24秒切片图，如图7为绕射波成像结果，包括了相应位置的剖面和切片图。
如图6和图7所示，绕射波成像结果在切片主测线和联络线交叉位置很好的显示了已验证陷落柱的形态，而该陷落柱在反射波成像切片上并没有出现。
本发明实施例提供的一种煤炭陷落柱识别方法和装置，与现有技术中的针对反射波研发的地震波成像技术进行陷落柱识别的效果不佳相比，其首先计算地震炮集数据中每个单炮数据对应的绕射波走时，并通过马氏距离对单炮数据及单炮数据的绕射波走时进行计算处理，获取所有单炮数据对应的绕射波，并对绕射波进行成像处理；由于绕射波携带的是更小尺度的地质信息，故通过成像后的绕射波能够准确的评价陷落柱地质体，即能够准确的进行目标区域中煤炭陷落柱的识别，从而能够降低煤炭开采中由陷落柱而引发的突水、瓦斯泄露等事故发生风险，减少不必要的人员伤亡及经济损失。
本发明实施例还提供了一种煤炭陷落柱识别装置，参考图8，所述装置包括：
获取模块11，用于获取目标区域的地震炮集数据和所述地震炮集数据对应的地震波偏移速度文件；其中，所述地震炮集数据包括多个单炮数据；每个所述单炮数据均包括地下岩层界面反射或者折射的多个地震波；所述地震波携带有地质信息；
绕射波走时计算模块12，用于根据每个单炮数据对应的检波器空间位置和地震波偏移速度文件，计算由每个单炮数据的地震波炮点位置经地下成像空间任意一个成像点位置到每个单炮数据的地震波检波点位置的每个单炮数据的绕射波走时；
马氏距离计算处理模块13，用于对每个单炮数据及单炮数据的绕射波走时进行马氏距离计算处理，获取每个单炮数据的绕射波振幅值样点；其中，绕射波振幅值样点携带有小尺度地质信息；小尺度地质信息至少包括：地层层位信息、断层信息和陷落柱信息；
成像处理模块14，用于分别对每个单炮数据的绕射波振幅值样点进行成像处理，得到每个单炮数据的绕射波成像结果；
叠加处理模块15，用于将地震炮集数据对应的所有单炮数据的绕射波成像结果进行叠加处理，得到地震炮集数据的绕射波成像结果，以便根据地震炮集数据的绕射波成像结果进行煤炭陷落柱识别。
进一步的，参考图9，所述基于地震绕射波成像的煤炭陷落柱识别装置中，绕射波走时计算模块12包括：
第一走时计算单元121，用于根据每个单炮数据对应的检波器空间位置和地震波偏移速度文件，计算由每个单炮数据的地震波炮点位置到每一个地下成像空间成像点位置的走时tS；
第二走时计算单元122，用于根据每个单炮数据对应的检波器空间位置和地震波偏移速度文件，计算由每一个地下成像空间成像点位置到每个单炮数据的地震波检波点位置的走时tR；其中，每个单炮数据在一个成像点的走时tR与走时tS相对应；
第一求和处理单元123，用于分别对每个单炮数据的走时tS以及与走时tS对应的走时tR进行求和处理，得到对应于每个单炮数据的多个绕射波走时。
进一步的，参考图9，所述基于地震绕射波成像的煤炭陷落柱识别装置中，马氏距离计算处理模块13包括：
振幅值样点序列计算单元131，用于对于任意一个单炮数据，根据单炮数据和单炮数据的多个绕射波走时，计算单炮数据对应的地下成像空间中每个成像点的振幅值样点序列，得到对应于单炮数据的多个振幅值样点序列；
排序处理单元132，用于根据振幅值样点的能量大小，分别对得到的每一个振幅值样点序列中的振幅值样点进行排序处理；
马氏距离计算单元133，用于分别对排序处理后的每一个振幅值样点序列进行马氏距离计算处理，得到每一个振幅值样点序列中按预设顺序排列的所有振幅值样点的马氏距离；
提取单元134，用于根据每个振幅值样点序列中按预设顺序排列的所有振幅值样点的马氏距离，提取对应的每个振幅值样点序列中符合预设条件的振幅值样点，得到对应于每个所述单炮数据的多组绕射波振幅值样点。
进一步的，所述基于地震绕射波成像的煤炭陷落柱识别装置中，马氏距离计算单元133包括：
第一计算子单元，用于通过以下公式计算每一个按预设顺序排列的振幅值样点序列的中值其中，μx为振幅值样点序列的中值，median为中值运算，为按预设顺序排序后的振幅样点序列，N为振幅值样点序列中样点个数；
第二计算字单元，用于通过以下公式计算每一个按预设顺序排列的振幅值样点序列的中值绝对值偏差平方其中，为中值绝对值偏差；
第三计算子单元，用于根据中值和中值绝对值偏差平方，计算每一个振幅值样点序列中按预设顺序排列的所有振幅值样点的马氏距离 d ( x ( i ) o , μ x ) = [ x ( i ) o - μ x ] 2 σ x 2 ; 其中，为马氏距离。
进一步的，所述基于地震绕射波成像的煤炭陷落柱识别装置中，提取单元334包括：
第四计算子单元，用于对于任意一个振幅值样点序列，根据振幅值样点序列中按预设顺序排列的每个振幅值样点的马氏距离以及反射波振幅值样点的预设马氏距离阈值，计算振幅值样点序列中的样点去除数量；
样点去除处理子单元，用于根据所述预设马氏距离阈值和所述样点去除数量，对所述振幅值样点序列中按预设顺序排列的振幅值样点进行样点去除处理，得到所述振幅值样点序列中剩余的振幅值样点；
提取子单元，用于提取所述剩余的振幅值样点作为所述振幅值样点序列中符合预设条件的振幅值样点。
进一步的，所述基于地震绕射波成像的煤炭陷落柱识别装置中，成像处理模块14包括：
第二求和处理单元，用于将在每个振幅值样点序列中提取的符合预设条件的振幅值样点分别进行求和处理，得到每个振幅值样点序列对应的绕射波成像结果；
叠加处理单元，用于将所述单炮数据对应的所有振幅值样点序列的绕射波成像结果进行叠加处理，得到对应于所述单炮数据的绕射波成像结果。
进一步的，所述基于地震绕射波成像的煤炭陷落柱识别装置中，所述第二求和处理单元还用于，通过以下公式分别对每个单炮数据的绕射波进行成像处理 V ( m ) = - 1 2 π ∫ ∫ A Tri M { w ( m , r ) u ( r , t S + t R ) } d r ; 其中，V(m)为绕射波成像结果，m＝m(x,y,z)为地下成像空间每个成像点位置，u(r,tS+tR)为预处理后的地震炮集数据，r(x,y,z)为检波点位置，tS,tR分别为由炮点位置到成像点位置和由成像点位置到检波器位置的走时，A为偏移成像孔径，w(m,r)为几何扩散因子，TriM为基于马氏距离和反射波振幅去除的绕射波成像算子。
进一步的，所述基于地震绕射波成像的煤炭陷落柱识别装置中，获取模块11包括：
获取单元，用于获取目标区域的地震炮集数据；
地震预处理单元，用于对地震炮集数据进行地震预处理，得到可用于偏移成像的地震炮集数据；其中，地震预处理至少包括：去除噪声处理和静校正处理；
偏移速度分析处理单元，用于对地震炮集数据进行偏移速度分析处理，得到地震炮集数据对应的地震波偏移速度文件。
本发明实施例提供的一种煤炭陷落柱识别装置，与现有技术中的针对反射波研发的地震波成像技术进行陷落柱识别的效果不佳相比，其首先计算地震炮集数据中每个单炮数据对应的绕射波走时，并通过马氏距离对单炮数据及单炮数据的绕射波走时进行计算处理，获取所有单炮数据对应的绕射波，并对绕射波进行成像处理；由于绕射波携带的是更小尺度的地质信息，故通过成像后的绕射波能够准确的评价陷落柱地质体，即能够准确的进行目标区域中煤炭陷落柱的识别，从而能够降低煤炭开采中由陷落柱而引发的突水、瓦斯泄露等事故发生风险，减少不必要的人员伤亡及经济损失。
在本申请所提供的几个实施例中，应该理解到，所揭露的系统、装置和方法，可以通过其它的方式实现。以上所描述的装置实施例仅仅是示意性的，例如，所述单元的划分，仅仅为一种逻辑功能划分，实际实现时可以有另外的划分方式，又例如，多个单元或组件可以结合或者可以集成到另一个系统，或一些特征可以忽略，或不执行。另一点，所显示或讨论的相互之间的耦合或直接耦合或通信连接可以是通过一些通信接口，装置或单元的间接耦合或通信连接，可以是电性，机械或其它的形式。
所述作为分离部件说明的单元可以是或者也可以不是物理上分开的，作为单元显示的部件可以是或者也可以不是物理单元，即可以位于一个地方，或者也可以分布到多个网络单元上。可以根据实际的需要选择其中的部分或者全部单元来实现本实施例方案的目的。
另外，在本发明各个实施例中的各功能单元可以集成在一个处理单元中，也可以是各个单元单独物理存在，也可以两个或两个以上单元集成在一个单元中。
所述功能如果以软件功能单元的形式实现并作为独立的产品销售或使用时，可以存储在一个计算机可读取存储介质中。基于这样的理解，本发明的技术方案本质上或者说对现有技术做出贡献的部分或者该技术方案的部分可以以软件产品的形式体现出来，该计算机软件产品存储在一个存储介质中，包括若干指令用以使得一台计算机设备(可以是个人计算机，服务器，或者网络设备等)执行本发明各个实施例所述方法的全部或部分步骤。而前述的存储介质包括：U盘、移动硬盘、只读存储器(ROM，Read-Only Memory)、随机存取存储器(RAM，Random Access Memory)、磁碟或者光盘等各种可以存储程序代码的介质。
以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应以所述权利要求的保护范围为准。

Claims (10)

1.一种煤炭陷落柱识别方法，其特征在于，包括：
获取目标区域的地震炮集数据和所述地震炮集数据对应的地震波偏移速度文件；其中，所述地震炮集数据包括多个单炮数据；每个所述单炮数据均包括地下岩层界面反射或者折射的多个地震波；所述地震波携带有地质信息；
根据每个所述单炮数据对应的检波器空间位置和所述地震波偏移速度文件，计算由每个单炮数据的地震波炮点位置经地下成像空间任意一个成像点位置到每个单炮数据的地震波检波点位置的每个单炮数据的绕射波走时；
对每个所述单炮数据及所述单炮数据的绕射波走时进行马氏距离计算处理，获取每个所述单炮数据的绕射波振幅值样点；其中，所述绕射波振幅值样点携带有小尺度地质信息；所述小尺度地质信息至少包括：地层层位信息、断层信息和陷落柱信息；
分别对每个所述单炮数据的绕射波振幅值样点进行成像处理，得到每个所述单炮数据的绕射波成像结果；
将所述地震炮集数据对应的所有单炮数据的绕射波成像结果进行叠加处理，得到所述地震炮集数据的绕射波成像结果，以便根据所述地震炮集数据的绕射波成像结果进行煤炭陷落柱识别。
2.根据权利要求1所述的方法，其特征在于，所述根据每个所述单炮数据对应的检波器空间位置和所述地震波偏移速度文件，计算由每个单炮数据的地震波炮点位置经地下成像空间任意一个成像点位置到每个单炮数据的地震波检波点位置的每个单炮数据的绕射波走时包括：
根据每个所述单炮数据对应的检波器空间位置和所述地震波偏移速度文件，计算由每个单炮数据的地震波炮点位置到每一个地下成像空间成像点位置的走时tS；
根据每个所述单炮数据对应的检波器空间位置和所述地震波偏移速度文件，计算由每一个地下成像空间成像点位置到每个单炮数据的地震波检波点位置的走时tR；其中，每个所述单炮数据在一个成像点的走时tR与走时tS相对应；
分别对每个单炮数据的走时tS以及与所述走时tS对应的走时tR进行求和处理，得到对应于每个所述单炮数据的多个绕射波走时。
3.根据权利要求2所述的方法，其特征在于，所述对每个所述单炮数据及所述单炮数据的绕射波走时进行马氏距离计算处理，获取每个所述单炮数据的绕射波振幅值样点包括：
对于任意一个单炮数据，根据所述单炮数据和所述单炮数据的多个绕射波走时，计算所述单炮数据对应的地下成像空间中每个成像点的振幅值样点序列，得到对应于所述单炮数据的多个振幅值样点序列；
根据振幅值样点的能量大小，分别对得到的每一个振幅值样点序列中的振幅值样点进行排序处理；
分别对排序处理后的每一个振幅值样点序列进行马氏距离计算处理，得到每一个振幅值样点序列中按预设顺序排列的所有振幅值样点的马氏距离；
根据每个振幅值样点序列中按预设顺序排列的所有振幅值样点的马氏距离，提取对应的每个振幅值样点序列中符合预设条件的振幅值样点，得到对应于每个所述单炮数据的多组绕射波振幅值样点。
4.根据权利要求3所述的方法，其特征在于，所述分别对排序处理后的每一个振幅值样点序列进行马氏距离计算处理，得到每一个振幅值样点序列中按预设顺序排列的所有振幅值样点的马氏距离包括：
通过以下公式计算每一个按预设顺序排列的振幅值样点序列的中值其中，μx为振幅值样点序列的中值，median为中值运算，i＝1,2,…,N为按预设顺序排序后的振幅样点序列，N为振幅值样点序列中样点个数；
通过以下公式计算每一个按预设顺序排列的所述振幅值样点序列的中值绝对值偏差平方其中，为中值绝对值偏差平方；
根据所述中值和所述中值绝对值偏差平方，计算每一个振幅值样点序列中按预设顺序排列的所有振幅值样点的马氏距离其中，为马氏距离。
5.根据权利要求4所述的方法，其特征在于，所述根据每个振幅值样点序列中按预设顺序排列的所有振幅值样点的马氏距离，提取对应的每个振幅值样点序列中符合预设条件的振幅值样点包括：
对于任意一个振幅值样点序列，根据所述振幅值样点序列中按预设顺序排列的每个振幅值样点的马氏距离以及反射波振幅值样点的预设马氏距离阈值，计算所述振幅值样点序列中的样点去除数量；
根据所述预设马氏距离阈值和所述样点去除数量，对所述振幅值样点序列中按预设顺序排列的振幅值样点进行样点去除处理，得到所述振幅值样点序列中剩余的振幅值样点；
提取所述剩余的振幅值样点作为所述振幅值样点序列中符合预设条件的振幅值样点。
6.根据权利要求5所述的方法，其特征在于，所述分别对每个所述单炮数据的绕射波振幅值样点进行成像处理，得到每个所述单炮数据的绕射波成像结果包括：
将在每个振幅值样点序列中提取的符合预设条件的振幅值样点分别进行求和处理，得到每个振幅值样点序列对应的绕射波成像结果；
将所述单炮数据对应的所有振幅值样点序列的绕射波成像结果进行叠加处理，得到对应于所述单炮数据的绕射波成像结果。
7.根据权利要求6所述的方法，其特征在于，所述将在每个振幅值样点序列中提取的符合预设条件的振幅值样点分别进行求和处理，得到每个振幅值样点序列对应的绕射波成像结果包括：
通过以下公式将在每个振幅值样点序列中提取的符合预设条件的振幅值样点分别进行求和处理得到每个振幅值样点序列对应的绕射波成像结果；其中，V(m)为绕射波成像结果，m＝m(x,y,z)为地下成像空间每个成像点位置，u(r,tS+tR)为预处理后的地震炮集数据，r(x,y,z)为检波点位置，tS,tR分别为由炮点位置到成像点位置和由成像点位置到检波器位置的走时，A为偏移成像孔径，w(m,r)为几何扩散因子，TriM为基于马氏距离和反射波振幅去除的绕射波成像算子。
8.根据权利要求1所述的方法，其特征在于，所述获取目标区域的地震炮集数据和所述地震炮集数据对应的地震波偏移速度文件包括：
获取目标区域的地震炮集数据；
对所述地震炮集数据进行地震预处理，得到可用于偏移成像的地震炮集数据；其中，所述地震预处理至少包括：去除噪声处理和静校正处理；
对所述地震炮集数据进行偏移速度分析处理，得到所述地震炮集数据对应的地震波偏移速度文件。
9.一种煤炭陷落柱识别装置，其特征在于，所述装置包括：
获取模块，用于获取目标区域的地震炮集数据和所述地震炮集数据对应的地震波偏移速度文件；其中，所述地震炮集数据包括多个单炮数据；每个所述单炮数据均包括地下岩层界面反射或者折射的多个地震波；所述地震波携带有地质信息；
绕射波走时计算模块，用于根据每个所述单炮数据对应的检波器空间位置和所述地震波偏移速度文件，计算由每个单炮数据的地震波炮点位置经地下成像空间任意一个成像点位置到每个单炮数据的地震波检波点位置的每个单炮数据的绕射波走时；
马氏距离计算处理模块，用于对每个所述单炮数据及所述单炮数据的绕射波走时进行马氏距离计算处理，获取每个所述单炮数据的绕射波振幅值样点；其中，所述绕射波振幅值样点携带有小尺度地质信息；所述小尺度地质信息至少包括：地层层位信息、断层信息和陷落柱信息；
成像处理模块，用于分别对每个所述单炮数据的绕射波振幅值样点进行成像处理，得到每个所述单炮数据的绕射波成像结果；
叠加处理模块，用于将所述地震炮集数据对应的所有单炮数据的绕射波成像结果进行叠加处理，得到所述地震炮集数据的绕射波成像结果，以便根据所述地震炮集数据的绕射波成像结果进行煤炭陷落柱识别。
10.根据权利要求9所述的装置，其特征在于，所述绕射波走时计算模块包括：
第一走时计算单元，用于根据每个所述单炮数据对应的检波器空间位置和所述地震波偏移速度文件，计算由每个单炮数据的地震波炮点位置到每一个地下成像空间成像点位置的走时tS；
第二走时计算单元，用于根据每个所述单炮数据对应的检波器空间位置和所述地震波偏移速度文件，计算由每一个地下成像空间成像点位置到每个单炮数据的地震波检波点位置的走时tR；其中，每个所述单炮数据在一个成像点的走时tR与走时tS相对应；
第一求和处理单元，用于分别对每个单炮数据的走时tS以及与所述走时tS对应的走时tR进行求和处理，得到对应于每个所述单炮数据的多个绕射波走时。
CN201610073638.7A 2016-02-02 2016-02-02 一种煤炭陷落柱识别方法和装置 Active CN105607121B (zh)

Priority Applications (3)

Application Number Priority Date Filing Date Title
CN201610073638.7A CN105607121B (zh) 2016-02-02 2016-02-02 一种煤炭陷落柱识别方法和装置
US15/758,559 US10302788B2 (en) 2016-02-02 2016-12-29 Method and apparatus for identifying collapsed coal column
PCT/CN2016/113064 WO2017133361A1 (zh) 2016-02-02 2016-12-29 一种煤炭陷落柱识别方法和装置

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201610073638.7A CN105607121B (zh) 2016-02-02 2016-02-02 一种煤炭陷落柱识别方法和装置

Publications (2)

Publication Number Publication Date
CN105607121A CN105607121A (zh) 2016-05-25
CN105607121B true CN105607121B (zh) 2016-12-21

Family

ID=55987197

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201610073638.7A Active CN105607121B (zh) 2016-02-02 2016-02-02 一种煤炭陷落柱识别方法和装置

Country Status (3)

Country Link
US (1) US10302788B2 (zh)
CN (1) CN105607121B (zh)
WO (1) WO2017133361A1 (zh)

Families Citing this family (28)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105607121B (zh) 2016-02-02 2016-12-21 中国矿业大学(北京) 一种煤炭陷落柱识别方法和装置
CN106646662B (zh) * 2016-11-10 2018-09-18 中国矿业大学(北京) 瓦斯突出区域的预测方法和装置
CN107015275B (zh) * 2017-04-14 2019-04-19 中国矿业大学(北京) 陷落柱检测方法和装置
CN106990436B (zh) * 2017-04-14 2019-03-29 中国矿业大学(北京) 陷落柱的识别方法和装置
CN107450098B (zh) * 2017-08-03 2018-11-23 中煤科工集团西安研究院有限公司 一种煤层底板隐伏突水陷落柱动态定位方法
CN107817523B (zh) * 2017-10-30 2018-10-09 中国矿业大学(北京) 绕射波偏移速度的分析方法及装置
CN107918147B (zh) * 2017-11-20 2018-12-14 中国矿业大学(北京) 绕射波成像方法和装置
CN111381274A (zh) * 2018-12-29 2020-07-07 中国石油天然气股份有限公司 一种传输故障道识别方法及装置
CN109934409B (zh) * 2019-03-18 2022-12-27 西安科技大学 一种大采高煤壁工作面片帮预测网络及其预测方法
CN109839663B (zh) * 2019-03-20 2020-04-10 山西山地物探技术有限公司 一种隐伏陷落柱的地震识别方法和装置
CN110362948B (zh) * 2019-07-23 2023-07-18 河北省交通规划设计研究院有限公司 一种基于云模型的岩溶地面塌陷的治理方法
CN112305603B (zh) * 2019-08-02 2024-03-26 中国石油天然气股份有限公司 一种节点仪器质量监控方法及系统
CN110531414B (zh) * 2019-08-21 2020-10-30 中国矿业大学 一种高倾角多层界面的反射地震断层精细探测方法
CN112444869B (zh) * 2019-08-30 2024-06-18 中国石油化工股份有限公司 一种用于压制外源干扰波的地震数据处理方法及存储介质
CN111812714B (zh) * 2020-06-08 2022-11-15 中煤科工集团西安研究院有限公司 基于折射纵波与高频槽波的煤层纵横波速度求取方法
CN111929730B (zh) * 2020-08-20 2021-04-16 中国矿业大学(北京) 小尺度地质异常体检测方法及装置
CN111929729B (zh) * 2020-08-20 2021-04-06 中国矿业大学(北京) 绕射波成像方法、装置和电子设备
CN112230281A (zh) * 2020-09-17 2021-01-15 陕西省煤田地质集团有限公司 一种快速识别陷落柱的地震方法
CN112255679B (zh) * 2020-10-26 2023-09-26 中国石油天然气集团有限公司 地震资料绕射深度偏移处理方法及装置
CN112285777B (zh) * 2020-10-26 2024-07-26 中国石油天然气集团有限公司 地震道串感应识别量化方法及装置
CN112415591B (zh) * 2020-10-30 2024-08-27 中国石油天然气集团有限公司 绕射波的成像方法、装置、电子设备及存储介质
CN112379420B (zh) * 2020-10-30 2024-05-28 中国石油天然气集团有限公司 高精度弯曲测线叠前时间域成像方法及装置
CN112379428B (zh) * 2020-11-02 2024-08-27 中国石油天然气集团有限公司 地震数据一致性处理方法及装置
CN113052208B (zh) * 2021-03-10 2023-08-25 神华神东煤炭集团有限责任公司 基于视觉的煤岩识别方法、存储介质及电子设备
CN113156504B (zh) * 2021-04-09 2021-10-15 中国科学院地理科学与资源研究所 一种地震波速度的确定方法及装置
CN113608261B (zh) * 2021-07-30 2021-12-28 中国矿业大学(北京) 绕射波成像方法、装置和电子设备
CN114185082B (zh) * 2021-12-02 2023-04-21 中国矿业大学 一种基于工作面透射地震观测的煤层下伏陷落柱探测方法
CN117908122A (zh) * 2024-01-18 2024-04-19 广西有色勘察设计研究院 基于绕射波探测隧道围岩工程钻探孔旁地质情况的方法

Family Cites Families (11)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2004135868A (ja) * 2002-10-17 2004-05-13 Fuji Photo Film Co Ltd 異常陰影候補検出処理システム
WO2008111037A2 (en) * 2007-03-12 2008-09-18 Geomage 2003 Ltd A method for identifying and analyzing faults/fractures using reflected and diffracted waves
US7788070B2 (en) * 2007-07-30 2010-08-31 Caterpillar Inc. Product design optimization method and system
WO2009077440A2 (en) * 2007-12-14 2009-06-25 Shell Internationale Research Maatschappij B.V. Method of processing data obtained from seismic prospecting
CN103076631B (zh) 2011-10-26 2016-03-23 中国石油化工股份有限公司 一种基于零脉冲反褶积提频技术的煤层气田陷落柱预测方法
CN102520444B (zh) * 2011-12-13 2014-10-01 中国科学院地质与地球物理研究所 一种叠后地震波中绕射波信息提取方法
CN102841063B (zh) * 2012-08-30 2014-09-03 浙江工业大学 一种基于光谱技术的生物炭溯源鉴别方法
CN103984012B (zh) * 2014-04-16 2016-06-29 刘豫宝 基于叠前高斯束深度偏移的绕射波场分离方法
CN104237940B (zh) * 2014-09-29 2017-03-08 中国石油天然气股份有限公司 一种基于动力学特征的绕射波成像方法及装置
CN104730571A (zh) * 2015-03-11 2015-06-24 中国科学院地质与地球物理研究所 一种利用绕射再聚焦识别小尺度地质体的方法与装置
CN105607121B (zh) * 2016-02-02 2016-12-21 中国矿业大学(北京) 一种煤炭陷落柱识别方法和装置

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
地震勘探中的边缘绕射波及其动力学识别方法;赵惊涛 等;《地球物理学进展》;20110215;第26卷(第1期);194-206 *
地震绕射波振幅特征研究;崔伟雄 等;《能源技术与管理》;20140828;第39卷(第4期);9-11 *
陷落柱地震数值模拟响应特征分析;李之奋等;《能源技术与管理》;20071228(第06期);36-39 *

Also Published As

Publication number Publication date
US10302788B2 (en) 2019-05-28
CN105607121A (zh) 2016-05-25
US20180246241A1 (en) 2018-08-30
WO2017133361A1 (zh) 2017-08-10

Similar Documents

Publication Publication Date Title
CN105607121B (zh) 一种煤炭陷落柱识别方法和装置
Ma et al. Discrimination of seismic sources in an underground mine using full waveform inversion
CN104237940B (zh) 一种基于动力学特征的绕射波成像方法及装置
Wollin et al. A unified earthquake catalogue for the Sea of Marmara Region, Turkey, based on automatized phase picking and travel-time inversion: Seismotectonic implications
CN103926617B (zh) 一种缝洞储集体的检测方法和检测装置
WO2016141598A1 (zh) 共姿态道集地表一致性振幅补偿方法及装置
McBrearty et al. Earthquake arrival association with backprojection and graph theory
Gvishiani et al. FCAZm intelligent recognition system for locating areas prone to strong earthquakes in the Andean and Caucasian mountain belts
CN107783185B (zh) 一种层析静校正的处理方法及装置
CN103792576A (zh) 一种基于梯度结构张量的储层非均质检测方法及设备
CN104122587B (zh) 一种基于炮检距域的异常初至识别方法以及系统
CN105785439A (zh) 小尺度非均匀地质体空间分布位置的预测方法和装置
WO2012145113A1 (en) Methods for deblending of seismic shot gathers
Chamarczuk et al. Unsupervised learning used in automatic detection and classification of ambient‐noise recordings from a large‐N array
CN113189644B (zh) 一种微震震源定位方法及系统
Vinard et al. Localizing microseismic events on field data using a U-Net-based convolutional neural network trained on synthetic data
CN112305591A (zh) 隧道超前地质预报方法、计算机可读存储介质
CN102914799A (zh) 非等效体波场正演模拟方法及装置
Jiang et al. Microseismic event classification with time-, frequency-, and wavelet-domain convolutional neural networks
Nguyen et al. Spatial cluster and outlier identification of geochemical association of elements: A case study in Juirui copper mining area
Tobely et al. Position detection of unexploded ordnance from airborne magnetic anomaly data using 3-D self organized feature map
CN108693560B (zh) 一种基于互相关道的散射波成像方法及系统
CN115690772A (zh) 基于深度学习的重力断层自动识别方法
CN115331096B (zh) 一种采矿沉陷区域识别方法、系统、存储介质和电子设备
Hofman et al. A new seismicity catalogue of the eastern Alps using the temporary Swath-D network

Legal Events

Date Code Title Description
C06 Publication
PB01 Publication
C10 Entry into substantive examination
SE01 Entry into force of request for substantive examination
C14 Grant of patent or utility model
GR01 Patent grant