CN110187386A - 一种自动快速识别地质构造的dtw地震体属性分析方法 - Google Patents
一种自动快速识别地质构造的dtw地震体属性分析方法 Download PDFInfo
- Publication number
- CN110187386A CN110187386A CN201910441148.1A CN201910441148A CN110187386A CN 110187386 A CN110187386 A CN 110187386A CN 201910441148 A CN201910441148 A CN 201910441148A CN 110187386 A CN110187386 A CN 110187386A
- Authority
- CN
- China
- Prior art keywords
- dtw
- data
- xline
- inline
- seismic
- 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
Links
- 238000004458 analytical method Methods 0.000 title claims abstract description 22
- 238000000034 method Methods 0.000 claims abstract description 23
- 238000013517 stratification Methods 0.000 claims abstract description 10
- 238000005070 sampling Methods 0.000 claims description 44
- 230000005284 excitation Effects 0.000 claims description 20
- BULVZWIRKLYCBC-UHFFFAOYSA-N phorate Chemical compound CCOP(=S)(OCC)SCSCC BULVZWIRKLYCBC-UHFFFAOYSA-N 0.000 claims description 10
- 230000015572 biosynthetic process Effects 0.000 claims description 9
- 238000010276 construction Methods 0.000 claims description 6
- 238000013508 migration Methods 0.000 claims description 6
- 230000005012 migration Effects 0.000 claims description 6
- 238000003786 synthesis reaction Methods 0.000 claims description 6
- 239000000284 extract Substances 0.000 claims description 5
- 230000011218 segmentation Effects 0.000 claims description 5
- 241001269238 Data Species 0.000 claims description 3
- 238000005452 bending Methods 0.000 claims description 3
- 230000002123 temporal effect Effects 0.000 claims description 3
- 238000007794 visualization technique Methods 0.000 claims description 3
- 230000005856 abnormality Effects 0.000 claims description 2
- 239000010410 layer Substances 0.000 abstract description 33
- 230000008901 benefit Effects 0.000 abstract description 3
- 239000011229 interlayer Substances 0.000 abstract description 2
- 238000003325 tomography Methods 0.000 description 5
- 238000010586 diagram Methods 0.000 description 4
- 238000011160 research Methods 0.000 description 3
- 230000002547 anomalous effect Effects 0.000 description 2
- 230000008859 change Effects 0.000 description 2
- 230000008569 process Effects 0.000 description 2
- 208000035126 Facies Diseases 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000004364 calculation method Methods 0.000 description 1
- 239000003245 coal Substances 0.000 description 1
- 230000001427 coherent effect Effects 0.000 description 1
- 238000000354 decomposition reaction Methods 0.000 description 1
- 230000007812 deficiency Effects 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 230000008676 import Effects 0.000 description 1
- 239000000463 material Substances 0.000 description 1
- 238000005065 mining Methods 0.000 description 1
- 238000012545 processing Methods 0.000 description 1
- 238000003672 processing method Methods 0.000 description 1
- 230000000630 rising effect Effects 0.000 description 1
- 238000012163 sequencing technique Methods 0.000 description 1
- 230000003595 spectral effect Effects 0.000 description 1
- 238000012800 visualization Methods 0.000 description 1
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/306—Analysis for determining physical properties of the subsurface, e.g. impedance, porosity or attenuation profiles
-
- 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
- G01V1/34—Displaying seismic recordings or visualisation of seismic data or attributes
- G01V1/345—Visualisation of seismic data or attributes, e.g. in 3D cubes
-
- 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
- G01V2210/512—Pre-stack
-
- 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
- G01V2210/514—Post-stack
-
- 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/63—Seismic attributes, e.g. amplitude, polarity, instant phase
-
- 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
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V2210/00—Details of seismic processing or analysis
- G01V2210/70—Other details related to processing
- G01V2210/74—Visualisation of seismic data
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
本发明公开了一种自动快速识别地质构造的DTW地震体属性分析方法，采用分布式并行算法，划分工作小区，在粗网格密度下，通过计算各地震道与模型参考道特征点序列DTW偏移时间量，能够快速解释层位和构建地震等时地层格架，在井分层的基础上进行层间等比例法剖分形成地震分层数据体，减少DTW计算的样点数，保证数据之间独立性，可有效发挥分布式运算的优势，大幅缩短DTW属性的计算时间。采用Inline和Xline两个方向进行相邻道的DTW计算，突出目标层段地震数据之间的弱连续性、低相似性的异常特征，从而实现不同尺度构造的精细识别。
Description
技术领域
本发明一种自动快速识别地质构造的DTW地震体属性分析方法，属于地震勘探技术领域。
背景技术
地震属性是指从人工地震时间序列中提取的运动学、动力学及统计学的度量值，与初始地震剖面数据相比，属性数据可突出断层、陷落柱等地质构造引起的地震波旅行时、振幅、频率和相位等异常的横向、垂向的不连续性和低相似性。根据属性维度，可分为层面二维和三维体两类。二维属性通常是沿某个层面计算的瞬时或固定小时窗属性，体属性则是三维空间的数据体。基于三瞬属性、频谱分解、方差或相干体属性、蚂蚁体、断层概率体等三维体属性及其组合分析，受人工干预较少，成果更客观。但是，目前上述属性均采用比较两个等长的时间序列之间相似性的计算思路，而DTW （动态时间弯曲）算法由于能够比较两个不等长的时间序列之间的相似性，算法更先进，从而受到了地球物理业界的广泛关注，称之为“动态波形匹配”或“动态时间规整”算法。
2003年，Keysers等提出了DTW算法在密集采样数据中求取精确解的高耗时这一核心问题，针对这一问题，国内外学者进行了大量的改进算法的研究，并应用于地震资料处理、解释。在资料处理方面，Hale（2013,2014）研究了平滑动态图像规整算法，认为在信号变化较慢的情况下更容易获得准确解；国内的研究主要集中在叠前道集拉平技术（钱峰等，2017）、纵波或横波多波资料匹配处理方法（陈茂山，2017；姚兴苗等，2018）、噪声压制（逯宇佳等,2018）、高精度速度分析（曹俊兴等,2018）等领域。在资料解释方面，则主要用于层位自动匹配追踪（曹俊兴等，2018）、地震波形分类（昌艳，2016）及地震相解释（钱峰等，2018）等方面，涉及地质构造精细识别及相关属性分析的研究极少。
目前利用反射地震勘探资料，通过地震多属性分析技术进行构造解释,可以查明落差5m及以上的断层、30m及以上直径的陷落柱，然而距离现代煤矿化综合机械化开采的要求仍相差甚远，小于上述尺度的断层、陷落柱等构造解释的多解性强，受经验因素影响大。
动态时间弯曲算法（DTW）能够更精细分析两个时间序列的相似性，但是，其精细解随着时间序列长度变大呈几何级数增长，计算效率极低，对于三维地震海量数据还未见可行的解决方案。
发明内容
本发明为克服现有技术的不足，目的是提供一种自动快速识别地质构造的DTW地震体属性分析方法，采用分布式并行算法，快速完成海量三维地震数据的DTW属性体计算，提高DTW地震属性体分析的可行性和达到商业化标准，并降低了小尺度地质构造和异常体解释的多解性和人为因素影响。
本发明通过以下技术方案实现：
一种自动快速识别地质构造的DTW地震体属性分析方法，包括以下步骤：
步骤1，结合井网密度进行工作区分割，利用基准井的井震层位标定结果，通过粗网格数据的DTW计算，完成层位自动追踪，快速构建整个工作区的地震等时地层格架；
步骤2，根据步骤1得到的地震等时地层格架，采用分布式并行算法，在每个工作小区内，采用等比例法分离出相邻两层位之间地震道数据，构建不等样点数量的地震分层数据体，每个分层体在Inline和Xline两个方向分别按照相邻道对比的方法，计算DTW时间偏移量；并将其按分层体的起止时间顺序存放；
步骤3，将不同分层体数据进行垂向叠加，再选取Inline和Xline方向权重系数，计算获得各工作小区综合Inline和Xline两个方向的DWT属性数据体，最后将不同工作小区DTW属性数据体进行横向叠加，获得全工作区的DTW属性体；
步骤4，将整个工作区的DWT属性体数据导入地震资料解释软件，与常规地震属性分析方法相仿，采用地层切片或三维可视化技术，种子点追踪方法，实现地质构造的三维自动追踪和解释。
优选地，所述步骤1包括如下步骤：
步骤1.1，利用地震叠后或叠前偏移数据，进行工作区井震标定，根据标定结果，确定每口井中多个目标反射层的时间信息；
步骤1.2，根据整个工作区井位分布情况，建立不小于500m×500m网度的井网；基于井间距的1/3-1/2，按Inline和Xline线道号将整个工作区分割为不同工作小区，记作B 1 ， B 2 ，……Bm，标记其起止inline号和Xline号；每个区域中优选一口反射层位完整的井，作为基准井，记作BW 1 ，BW 2 ，……，BW m ；
步骤1.3，在基准井处选取井旁道，提取极值和过零点等特征值，包括极大值与极小值，正过零点和负过零点，并将所述特征值按时间依次排序，建立各井区的模型道，记作GM w1 ， GM w2 ，…… ，GM wm ，存储每个特征点的时间、振幅参数，样点数控制在200个以内；存储每个模型道中各反射层T 1 ,T 2 , ……Tn对应的自激自收时间T 0 ，作为基准井中每个地质分层体的起止时间；
步骤1.4，在全区建立不小于100m×100m的层位解释网格，解释线尽可能过基准井，以保证分层结果与井震联合标定结果相匹配；对解释网格上的每个地震道提取极值和过零点等特征值，包括极大值与极小值，正过零点和负过零点，并将所述特征值按时间依次排序，分别形成各井区的目标对比道集；存储每个地震道中所有特征点的Inline和Xline号以及自激自收时间T 0 和振幅等4个参数；
步骤1.5，在每个工作小区，基于各自的模型道，对每个对比道求DTW距离（式1），求取每个地震道相对于模型道上各个特征点的最佳动态弯曲匹配时间，使特征点与模型参考道各样点完全对齐，存储每个对比道的Inline和Xline号，与模型道中各反射层对应的特征点自激自收时间；
所述DTW距离的递归算法如下：
（1）
其中：D base (x a ，y b )表示向量点x a 和y b 之间的欧氏距离，D dtw (x a-1，y b )表示向量点x a-1 和y b 之间的DTW距离，D dtw (x a-1，y b-1)表示向量点x a-1 和y b-1 之间的DTW距离，D dtw (x a ，y b-1)表示向量点x a 和y b-1 之间的DTW距离；
步骤1.6，根据各分区基准井模型参考道上标定的反射波层时间，按照振幅值相同的要求，搜索获取每个对比道上与模型道各反射层相对应的样点值及匹配时间，按式2求取各对比道中与各反射层的对应的自激自收时间T 0 ；
步骤1.7，将每个对比道中各反射层的T0时间，输入地震解释软件，进行异常值剔除，平滑，完成粗网格层位三维闭合解释；进一步在整个工作区做平面插值，建立全区的地震等时地层格架。
所述各对比道中与各反射层的对应的自激自收时间T 0 通过如下公式求得：
（2）
式中，i为Inline号，j为Xline号，n为层位号，B m 为工作小区号，BW m 为基准井号；T 0 (i，j，n)为i,j道n反射层的自激自收时间；T 0(B m ，BW m ，n)为B m工作小区内BW m 模型道上n反射层对应的T 0 时间；T dtw (i,j,n)为Bm工作小区内i，j道n反射层对应的动态弯曲匹配时间。
优选地，所述步骤2包括如下步骤：
步骤2.1，根据工作小区总数，建立m个计算机集群，在主机上发起一个job，选择某一个计算机集群负责相对应的Bm工作小区的DTW计算；设置要开启的总线程数（workers）；并保证集群中所有计算机具有一致的k个worker数量或CPU核数进行计算；
步骤2.2，采用等比例法剖分方法，在每个工作小区对应的计算机上，分割成n个相互独立的地震分层体，作为相应的独立的数据；各地震分层体每个地震道的样点数一般小于或等于20个样点，采样间隔与原地震数据相同，记作Bm(i,j,n),存储记录每个样点的Inline号（i）和Xline号(j)，自激自收时间T 0 和对应的振幅；
步骤2.3，在m台计算机上，k个worker同时对n个分层体的每道数据按Inline方向做相邻道间DTW计算；
步骤2.4，每个分层体的DTW属性计算结果，按分层体起止自激自收时间T 0 进行排序，起止时间外的数据全部置零，使每道的采样长度与原地震数据的采样长度相同；
步骤2.5，重复步骤2.3和2.4，完成Xline方向的每个工作小区的DTW属性体的计算。
所述2.2-2.4各步骤计算数据和存储相互独立，可以有效实现分布式并行运算。
优选地，所述步骤3包括如下步骤：
步骤3.1，以工作小区为单位，将各个分层体内部地震数据沿Inline方向进行DTW运算，并将不同分层体计算得到的有效数据在垂向上叠加，获得各工作小区Inline方向的全采样长度的DTW属性体，并对各样点的计算结果进行绝对值运算；
步骤3.2，以工作小区为单位，将各个分层体内部地震数据沿Xline方向进行DTW运算，并将不同分层体计算得到的有效数据在垂向上叠加，获得各工作小区Xline方向的全采样长度的DTW属性体，并对各样点的计算结果进行绝对值运算；
步骤3.3，以工作小区为单位，根据区域构造的走向与地震数据方向，选取相同或不同的权重系数，计算获得Inline和Xline两个方向的DWT属性综合体；
所述Inline和Xline两个方向的DWT属性综合体通过下式计算获得：
（3）
式中，i为Inline号，j为Xline号；DTW(i，j)为i，j道上每个样点对应的DTW偏移时间；DTW(i)为i，j道上每个样点与Inline方向相邻道上相同时间样点所对应的DTW偏移时间；DTW(j)为i，j道上每个样点与Xline方向相邻道相同时间样点所对应的DTW偏移时间；W(i)为Inline方向的构造权系数；W(j)为Xline方向的构造权系数；e为自然常数。
步骤3.4，将所有数据传回主机，按照各工作小区的Inline和Xline号进行横向数据拼合，获得全工作区的DWT属性体数据。
所述步骤3.1和3.2的数据相互独立，各个工作小区的计算过程可同时分布式实现。
DTW属性体可以是DTW偏移时间，也可以是DTW距离。
与现有技术相比，本发明具有如下有益效果：
本发明对海量地震数据进行区块划分，在粗网格上对各地震道与井模型参考道的特征点时间序列进行DTW运算，分布式并行算法，在与目前常见属性体计算耗时相当的时间内，能够快速完成数十万道乃至上百万道海量三维地震数据的DTW属性体计算，提高DTW地震属性体分析的可行性和达到商业化标准，能快速建立整个工作区的地震等时地层格架。经过垂向的地震分层数据体剥离，保证数据之间独立性，减少DTW计算样点数，充发发挥分布式并行运算的高效优势，大幅缩短DTW计算的时间。
本发明进一步沿Inline和Xline两个方向对相邻道进行相同采样时间序列的二次DTW计算，突出地震数据体中弱连续性、低相似性的小尺度构造异常特征。采用并通过多尺度、多类型地质构造的自动解释和三维可视化，降低小尺度地质构造和异常体解释的多解性和人为因素影响。
附图说明
图1为本发明的DTW地震属性体计算方法流程示意图；
图2为本发明的DTW地震属性体计算装置示意图；
图3为本发明工作小区及层位解释网格划分示意图；
图4为本发明地震分层数据体等比例剖分示意图；
图5为本发明实施例中分层数据体及DTW属性体异常三维示例图。
具体实施方式
下面结合具体实施例对本发明做进一步的详细说明，但是本发明的保护范围并不限于这些实施例，凡是不背离本发明构思的改变或等同替代均包括在本发明的保护范围之内。
一种自动快速识别地质构造的DTW地震体属性分析方法，包括以下步骤：
步骤1，结合井网密度进行工作区分割，利用基准井的井震层位标定结果，通过粗网格数据的DTW计算，完成层位自动追踪，快速构建整个工作区的地震等时地层格架；
所述步骤1包括如下步骤：
步骤1.1，利用地震叠后或叠前偏移数据，进行工作区井震标定，根据标定结果，确定每口井中多个目标反射层的时间信息；
步骤1.2，根据整个工作区井位分布情况，建立不小于500m×500m网度的井网；基于井间距的1/3-1/2，按Inline和Xline线道号将整个工作区分割为不同工作小区；如图3所示，记作B 1 ，B 2 ，……Bm，标记其起止inline号和Xline号；每个区域中优选一口反射层位完整的井，作为基准井，记作BW 1 ，BW 2 ，……，BW m ；
步骤1.3，在基准井处选取井旁道，提取极值和过零点等特征值，包括极大值与极小值，正过零点和负过零点，并将所述特征值按时间依次排序，建立各井区的模型道，记作GM w1 ， GM w2 ， …… ，GM wm ，存储每个特征点的时间、振幅参数，样点数控制在200个以内；存储每个模型道中各反射层T 1 ,T 2 , ……Tn对应的自激自收时间T 0 ，作为基准井中每个地质分层体的起止时间；
本发明是面向海量地震数据设计的，如果整个工作区的总地震道数如果小于100万道，在步骤1.2和步骤1.3中可不进行工作小区的划分，仅进行分层体的分割。
DTW属性体可以是DTW偏移时间，也可以是DTW距离。
步骤1.4，在全区建立不小于100m×100m的层位解释网格，解释线尽可能过基准井，以保证分层结果与井震联合标定结果相匹配；如图3所示；对解释网格上的每个地震道提取极值和过零点等特征值，包括极大值与极小值，正过零点和负过零点，并将所述特征值按时间依次排序，分别形成各井区的目标对比道集；存储每个地震道中所有特征点的Inline和Xline号以及自激自收时间T 0 和振幅等4个参数；
步骤1.5，在每个工作小区，基于各自的模型道，对每个对比道求DTW距离，式（1），求取每个地震道相对于模型道上各个特征点的最佳动态弯曲匹配时间，使特征点与模型参考道各样点完全对齐，存储每个对比道的Inline和Xline号，与模型道中各反射层对应的特征点自激自收时间；
所述DTW距离的算法如下：
（1）
其中：D base (x a ，y b )表示向量点x a 和y b 之间的欧氏距离，D dtw (x a-1，y b )表示x a-1 和y b 之间的DTW距离，D dtw (x a-1，y b-1)表示x a-1 和y b-1 之间的DTW距离，D dtw (x a ，y b-1)表示x a 和y b-1 之间的DTW距离；
步骤1.6，根据各分区基准井模型参考道上标定的反射波层时间，按照振幅值相同的要求，搜索获取每个对比道上与模型道各反射层相对应的样点值及匹配时间，按式（2）求取各对比道中与各反射层的对应的自激自收时间T 0 ；所述各对比道中与各反射层的对应的自激自收时间T 0 通过如下公式求得：
（2）
式中，i为Inline号，j为Xline号，n为层位号，B m 为工作小区号，BW m 为基准井号；T 0 (i，j，n)为i,j道n反射层的自激自收时间；T 0(B m ，BW m ，n)为B m工作小区内BW m 模型道上n反射层对应的T 0 时间；T dtw (i,j,n)为Bm工作小区内i，j道n反射层对应的动态弯曲匹配时间。
步骤1.7，将每个对比道中各反射层的T0时间，输入地震解释软件，进行异常值剔除，平滑，完成粗网格层位三维闭合解释；进一步在整个工作区做平面插值，建立全区的地震等时地层格架。
步骤2，根据步骤1得到的地震等时地层格架，采用分布式并行算法，在每个工作小区内，采用等比例法分离出相邻两层位之间地震道数据，构建不等样点数量的地震分层数据体，每个分层体在Inline和Xline两个方向分别按照相邻道对比的方法，计算DTW时间偏移量；并将其按分层体起止时间将各样点DTW时间顺序存放；
优选地，所述步骤2包括如下步骤：
步骤2.1，根据工作小区总数，建立m个计算机集群，在主机上发起一个job，选择某一个计算机集群负责相对应的Bm工作小区的DTW计算；设置要开启的总线程数（workers）；并保证集群中所有计算机具有一致的k个worker数量或CPU核数进行计算；
步骤2.2，采用如图4所示的等比例法剖分方法，在每个工作小区对应的计算机上，分割成n个相互独立的地震分层体，作为相应的独立的数据；各地震分层体每个地震道的样点数一般小于或等于20个样点，采样间隔与原地震数据相同，记作Bm(i,j,n),存储记录每个样点的Inline号（i）和Xline号(j)，自激自收时间T 0 和对应的振幅；
步骤2.3，在m台计算机上，k个worker同时对n个分层体的每道数据按Inline方向做相邻道间DTW计算；
步骤2.4，每个分层体的DTW属性计算结果，按分层体起止自激自收时间T 0 进行排序，起止时间外的数据全部置零，使每道的采样长度与原地震数据的采样长度相同；
步骤2.5，重复步骤2.3和2.4，完成Xline方向的每个工作小区的DTW属性体的计算。
所述2.2-2.4各步骤计算数据和存储相互独立，可以有效实现分布式并行运算。
步骤3，将不同分层体数据进行垂向叠加，再选取Inline和Xline方向权重系数，计算获得各工作小区综合Inline和Xline两个方向的DWT属性数据体，最后将不同工作小区DTW属性数据体进行横向叠加，获得全工作区的DTW属性体；
优选地，所述步骤3包括如下步骤：
步骤3.1，以工作小区为单位，将各个分层体内部地震数据沿Inline方向进行DTW运算，并将不同分层体计算得到的有效数据在垂向上叠加，获得各工作小区Inline方向的全采样长度的DTW属性体，并对各样点的计算结果进行绝对值运算；
步骤3.2，以工作小区为单位，将各个分层体内部地震数据沿Xline方向进行DTW运算，并将不同分层体计算得到的有效数据在垂向上叠加，获得各工作小区Xline方向的全采样长度的DTW属性体，并对各样点的计算结果进行绝对值运算；
步骤3.3，以工作小区为单位，根据区域构造的走向与地震数据方向，选取相同或不同的权重系数，计算获得Inline和Xline两个方向的DWT属性综合体；
所述Inline和Xline两个方向的DWT属性综合体通过下式计算获得：
（3）
式中，i为Inline号，j为Xline号；DTW(i，j)为i，j道上每个样点对应的DTW偏移时间；DTW(i)为i，j道上每个样点与Inline方向相邻道上相同时间样点所对应的DTW偏移时间；DTW(j)为i，j道上每个样点与Xline方向相邻道相同时间样点所对应的DTW偏移时间；W(i)为Inline方向的构造权系数；W(j)为Xline方向的构造权系数；e为自然常数。
步骤3.4，将所有数据传回主机，按照各工作小区的Inline和Xline号进行横向数据拼合，获得全工作区的DWT属性体数据。
所述步骤3.1和3.2的数据相互独立，各个工作小区的计算过程可同时分布式实现。
步骤4，将整个工作区的DWT属性体数据（标准sgy格式文件），导入地震资料解释软件，与常规地震属性分析方法相仿，采用地层切片或三维可视化技术，种子点追踪方法，实现地质构造的三维自动追踪和解释。
本发明划分工作小区，在粗网格密度下，通过计算各地震道与模型参考道特征点序列DTW偏移时间量，能够快速解释层位和构建地震等时地层格架。在井分层的基础上进行层间等比例法剖分形成地震分层数据体，减少DTW计算的样点数，保证数据之间独立性，可有效发挥分布式运算的优势，大幅缩短DTW属性的计算时间。采用Inline和Xline两个方向进行相邻道的DTW计算，突出目标层段地震数据之间的弱连续性、低相似性的异常特征，从而实现不同尺度构造的精细识别。
实施例为煤田三维地震工作小区，面积约7.9km2；叠前时间偏移数据体中共包含678,935道，每道1000个样点数；在不做分层体的情况下，采用8节点并行计算，总计算耗时约32小时；将上述数据依等比例剖分成10个分层体，每道样点数约92-135个；同样采用8节点并行计算，总计算耗时20.3分钟；将数据体依等比例剖分成20个分层体，每道样点数约6-67个，同样采用8节点计算，计算总耗时仅4.2分钟，这一结果与前二者相比分别减少了约455倍和3.8倍；表明本发明对于计算效率的提高非常明显。
如图5所示，图中多条不同落差断层（5-20m）的三维空间属性特征清晰，可用种子点法快速实现断层面的三维自动追踪。
本发明不会限制于本文所示的实施例，而是要符合与本文所公开的原理和新颖性特点相一致的最宽范围。
Claims (9)
1.一种自动快速识别地质构造的DTW地震体属性分析方法，其特征在于，包括以下步骤：
步骤1，结合井网密度进行工作区分割，利用基准井的井震层位标定结果，通过粗网格数据的DTW计算，完成层位自动追踪，快速构建整个工作区的地震等时地层格架；
步骤2，根据步骤1得到的地震等时地层格架，采用分布式并行算法，在每个工作小区内，采用等比例法分离出相邻两层位之间地震道数据，构建不等样点数量的地震分层数据体，每个分层体在Inline和Xline两个方向分别按照相邻道对比的方法，计算DTW时间偏移量；并将其按分层体的起止时间顺序存放；
步骤3，将不同分层体数据进行垂向叠加，再选取Inline和Xline方向权重系数，计算获得各工作小区综合Inline和Xline两个方向的DWT属性数据体，最后将不同工作小区DTW属性数据体进行横向叠加，获得全工作区的DTW属性体；
步骤4，将整个工作区的DWT属性体数据导入地震资料解释软件，采用地层切片或三维可视化技术，种子点追踪方法，实现地质构造的三维自动追踪和解释。
2.根据权利要求1所述的一种自动快速识别地质构造的DTW地震体属性分析方法，其特征在于，所述步骤1包括如下步骤：
步骤1.1，利用地震叠后或叠前偏移数据，进行工作区井震标定，根据标定结果，确定每口井中多个目标反射层的时间信息；
步骤1.2，根据整个工作区井位分布情况，建立不小于500m×500m网度的井网；基于井间距的1/3-1/2，按Inline和Xline线道号将整个工作区分割为不同工作小区，标记其起止inline号和Xline号；每个区域中优选一口反射层位完整的井，作为基准井；
步骤1.3，在基准井处选取井旁道，提取特征值，并将所述特征值按时间依次排序，建立各井区的模型道，存储每个特征点的时间、振幅参数，样点数控制在200个以内；存储每个模型道中各反射层对应的自激自收时间T 0 ，作为基准井中每个地质分层体的起止时间；
步骤1.4，在全区建立不小于100m×100m的层位解释网格，解释线过基准井，对解释网格上的每个地震道提取特征值，并将所述特征值按时间依次排序，分别形成各井区的目标对比道集；存储每个地震道中所有特征点的Inline和Xline号以及自激自收时间T 0 和振幅；
步骤1.5，在每个工作小区，基于各自的模型道，对每个对比道求DTW距离，求取每个地震道相对于模型道上各个特征点的最佳动态弯曲匹配时间，使特征点与模型参考道各样点完全对齐，存储每个对比道的Inline和Xline号，与模型道中各反射层对应的特征点自激自收时间；
步骤1.6，根据各分区基准井模型参考道上标定的反射波层时间，按照振幅值相同的要求，搜索获取每个对比道上与模型道各反射层相对应的样点值及匹配时间，求取各对比道中与各反射层的对应的自激自收时间T 0 ；
步骤1.7，将每个对比道中各反射层的T0时间，输入地震解释软件，进行异常值剔除，平滑，完成粗网格层位三维闭合解释；进一步在整个工作区做平面插值，建立全区的地震等时地层格架。
3.根据权利要求2所述的一种自动快速识别地质构造的DTW地震体属性分析方法，其特征在于，所述特征值包括但不限于极大值与极小值，正过零点和负过零点。
4.根据权利要求2所述的一种自动快速识别地质构造的DTW地震体属性分析方法，其特征在于，所述各对比道中与各反射层的对应的自激自收时间T 0 通过如下公式求得：
式中，i为Inline号，j为Xline号，n为层位号，B m 为工作小区号，BW m 为基准井号；T 0 (i，j，n)为i,j道n反射层的自激自收时间；T 0(B m ，BW m ，n)为B m工作小区内BW m 模型道上n反射层对应的T 0 时间；T dtw (i,j,n)为Bm工作小区内i，j道n反射层对应的动态弯曲匹配时间。
5.根据权利要求1所述的一种自动快速识别地质构造的DTW地震体属性分析方法，其特征在于，所述步骤2包括如下步骤：
步骤2.1，根据工作小区总数，建立m个计算机集群，在主机上发起一个job，选择某一个计算机集群负责相对应的Bm工作小区的DTW计算；设置要开启的总线程数workers；并保证集群中所有计算机具有一致的k个worker数量或CPU核数进行计算；
步骤2.2，采用等比例法剖分方法，在每个工作小区对应的计算机上，分割成n个相互独立的地震分层体，作为相应的独立的数据；各地震分层体每个地震道的样点数小于或等于20个样点，采样间隔与原地震数据相同，存储记录每个样点的Inline号和Xline号，自激自收时间T 0 和对应的振幅；
步骤2.3，在m台计算机上，k个worker同时对n个分层体的每道数据按Inline方向做相邻道间DTW计算；
步骤2.4，每个分层体的DTW属性计算结果，按分层体起止自激自收时间T 0 进行排序，起止时间外的数据全部置零，使每道的采样长度与原地震数据的采样长度相同；
步骤2.5，重复步骤2.3和2.4，完成Xline方向的每个工作小区的DTW属性体的计算。
6.根据权利要求5所述的一种自动快速识别地质构造的DTW地震体属性分析方法，其特征在于，所述2.2-2.4各步骤计算数据和存储相互独立。
7.根据权利要求1所述的一种自动快速识别地质构造的DTW地震体属性分析方法，其特征在于，所述步骤3包括如下步骤：
步骤3.1，以工作小区为单位，将各个分层体内部地震数据沿Inline方向进行DTW运算，并将不同分层体计算得到的有效数据在垂向上叠加，获得各工作小区Inline方向的全采样长度的DTW属性体，并对各样点的计算结果进行绝对值运算；
步骤3.2，以工作小区为单位，将各个分层体内部地震数据沿Xline方向进行DTW运算，并将不同分层体计算得到的有效数据在垂向上叠加，获得各工作小区Xline方向的全采样长度的DTW属性体，并对各样点的计算结果进行绝对值运算；
步骤3.3，以工作小区为单位，根据区域构造的走向与地震数据方向，选取相同或不同的权重系数，计算获得Inline和Xline两个方向的DWT属性综合体；
步骤3.4，将所有数据传回主机，按照各工作小区的Inline和Xline号进行横向数据拼合，获得全工作区的DWT属性体数据。
8.根据权利要求6所述的一种自动快速识别地质构造的DTW地震体属性分析方法，其特征在于，所述步骤3.1和3.2的数据相互独立。
9.根据权利要求6所述的一种自动快速识别地质构造的DTW地震体属性分析方法，其特征在于，所述Inline和Xline两个方向的DWT属性综合体通过下式计算获得：
式中，i为Inline号，j为Xline号；DTW(i，j)为i，j道上每个样点对应的DTW偏移时间；DTW(i)为i，j道上每个样点与Inline方向相邻道上相同时间样点所对应的DTW偏移时间；DTW(j)为i，j道上每个样点与Xline方向相邻道相同时间样点所对应的DTW偏移时间；W(i)为Inline方向的构造权系数；W(j)为Xline方向的构造权系数；e为自然常数。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910441148.1A CN110187386B (zh) | 2019-05-24 | 2019-05-24 | 一种自动快速识别地质构造的dtw地震体属性分析方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910441148.1A CN110187386B (zh) | 2019-05-24 | 2019-05-24 | 一种自动快速识别地质构造的dtw地震体属性分析方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN110187386A true CN110187386A (zh) | 2019-08-30 |
| CN110187386B CN110187386B (zh) | 2020-08-18 |
Family
ID=67717732
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201910441148.1A Active CN110187386B (zh) | 2019-05-24 | 2019-05-24 | 一种自动快速识别地质构造的dtw地震体属性分析方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN110187386B (zh) |
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110749929A (zh) * | 2019-10-28 | 2020-02-04 | 山东科技大学 | 一种基于机器学习的复杂地区地震反射层识别与追踪方法 |
| CN112511162A (zh) * | 2020-11-11 | 2021-03-16 | 许继集团有限公司 | 一种模拟量采集动态补偿方法及系统 |
Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102012524A (zh) * | 2010-09-29 | 2011-04-13 | 中国海洋石油总公司 | 一种海上三维地震观测系统羽状漂移定量评估方法 |
| CN102053270A (zh) * | 2009-10-30 | 2011-05-11 | 中国石油化工股份有限公司 | 一种基于沉积地层单元的地震相分析方法 |
| CN105243388A (zh) * | 2015-09-09 | 2016-01-13 | 电子科技大学 | 基于动态时间规整和划分算法的波形分类方法 |
| US20160377754A1 (en) * | 2015-06-24 | 2016-12-29 | Chevron U.S.A. Inc. | Well Placement Using Closure Stress Based Landing Map |
| CN109633746A (zh) * | 2017-10-09 | 2019-04-16 | 中国石油化工股份有限公司 | 一种反射倾角及其相似系数的自动拾取方法 |
| US10295688B2 (en) * | 2010-08-10 | 2019-05-21 | Westerngeco L.L.C. | Attenuating internal multiples from seismic data |
-
2019
- 2019-05-24 CN CN201910441148.1A patent/CN110187386B/zh active Active
Patent Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102053270A (zh) * | 2009-10-30 | 2011-05-11 | 中国石油化工股份有限公司 | 一种基于沉积地层单元的地震相分析方法 |
| US10295688B2 (en) * | 2010-08-10 | 2019-05-21 | Westerngeco L.L.C. | Attenuating internal multiples from seismic data |
| CN102012524A (zh) * | 2010-09-29 | 2011-04-13 | 中国海洋石油总公司 | 一种海上三维地震观测系统羽状漂移定量评估方法 |
| US20160377754A1 (en) * | 2015-06-24 | 2016-12-29 | Chevron U.S.A. Inc. | Well Placement Using Closure Stress Based Landing Map |
| CN105243388A (zh) * | 2015-09-09 | 2016-01-13 | 电子科技大学 | 基于动态时间规整和划分算法的波形分类方法 |
| CN109633746A (zh) * | 2017-10-09 | 2019-04-16 | 中国石油化工股份有限公司 | 一种反射倾角及其相似系数的自动拾取方法 |
Non-Patent Citations (1)
| Title |
|---|
| 于会臻: "反射系数反演方法研究及其在车排子薄储层描述中的应用", 《中国地球科学联合学术年会 2017》 * |
Cited By (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110749929A (zh) * | 2019-10-28 | 2020-02-04 | 山东科技大学 | 一种基于机器学习的复杂地区地震反射层识别与追踪方法 |
| CN110749929B (zh) * | 2019-10-28 | 2020-06-05 | 山东科技大学 | 一种基于机器学习的复杂地区地震反射层识别与追踪方法 |
| WO2021082080A1 (zh) * | 2019-10-28 | 2021-05-06 | 山东科技大学 | 一种基于机器学习的复杂地区地震反射层识别与追踪方法 |
| CN112511162A (zh) * | 2020-11-11 | 2021-03-16 | 许继集团有限公司 | 一种模拟量采集动态补偿方法及系统 |
| CN112511162B (zh) * | 2020-11-11 | 2023-05-02 | 许继集团有限公司 | 一种模拟量采集动态补偿方法及系统 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN110187386B (zh) | 2020-08-18 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| US6018498A (en) | Automated seismic fault detection and picking | |
| CN104297789B (zh) | 一种三维倾角域稳相叠前时间偏移方法及系统 | |
| US6574566B2 (en) | Automated feature identification in data displays | |
| US10073190B2 (en) | Method and system for geophysical modeling of subsurface volumes based on computed vectors | |
| CN105974479B (zh) | Gpu空间网格体的层析2d/3d各向异性深度域速度建模方法 | |
| CN103969683B (zh) | 一种三维地震解释中基于约束的批量拾取层位面的方法 | |
| US9529115B2 (en) | Geophysical modeling of subsurface volumes based on horizon extraction | |
| CN102176052B (zh) | 一种面向三维层面网格生成的层序分析方法 | |
| CN108897066A (zh) | 碳酸盐岩裂缝密度定量预测方法及装置 | |
| CN109100795A (zh) | 一种面元的炮检点布设方法、装置及系统 | |
| CN110187386A (zh) | 一种自动快速识别地质构造的dtw地震体属性分析方法 | |
| CN106415321A (zh) | 用于储层建模的基于瞬时等时属性的地质体识别 | |
| CN108828668A (zh) | 一种叠前时间偏移数据处理方法及装置 | |
| CN109765615A (zh) | 一种地层品质因子反演方法及装置 | |
| CN110187383A (zh) | 一种海上宽方位地震数据cov道集快速分选方法 | |
| CN110031898A (zh) | 数据优化方法及积分法叠前深度偏移方法 | |
| CN111505713B (zh) | 基于多点地质统计的叠前地震反演方法 | |
| CN108663714A (zh) | 一种沉积微相刻画方法 | |
| CN104267432B (zh) | 一种基于规则化的转换波共转换点道集高精度抽取方法 | |
| CN108226997A (zh) | 一种基于叠前地震数据的地震相划分方法 | |
| US7020558B2 (en) | Method of measuring local similarities between several seismic trace cubes | |
| CN105301638B (zh) | 一种提取风化层底界面的方法和装置 | |
| CN114861515A (zh) | 层速度数据体的计算方法、装置、设备及介质 | |
| CN102914790B (zh) | 二维观测系统和三维观测系统一次采集的观测系统方法 | |
| US5265068A (en) | Efficient generation of traveltime tables for complex velocity models for depth migration |
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 |