一种自动快速识别地质构造的DTW地震体属性分析方法

摘要

本发明公开了一种自动快速识别地质构造的DTW地震体属性分析方法，采用分布式并行算法，划分工作小区，在粗网格密度下，通过计算各地震道与模型参考道特征点序列DTW偏移时间量，能够快速解释层位和构建地震等时地层格架，在井分层的基础上进行层间等比例法剖分形成地震分层数据体，减少DTW计算的样点数，保证数据之间独立性，可有效发挥分布式运算的优势，大幅缩短DTW属性的计算时间。采用Inline和Xline两个方向进行相邻道的DTW计算，突出目标层段地震数据之间的弱连续性、低相似性的异常特征，从而实现不同尺度构造的精细识别。

说明书

技术领域

本发明一种自动快速识别地质构造的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₁，B₂，……Bm，标记其起止inline号和Xline号；每个区域中优选一口反射层位完整的井，作为基准井，记作BW₁，BW₂，……，BW_m；

步骤1.3，在基准井处选取井旁道，提取极值和过零点等特征值，包括极大值与极小值，正过零点和负过零点，并将所述特征值按时间依次排序，建立各井区的模型道，记作GM_w1，GM_w2，……>GM_wm，存储每个特征点的时间、振幅参数，样点数控制在200个以内；存储每个模型道中各反射层T₁,T₂,>Tn对应的自激自收时间T₀，作为基准井中每个地质分层体的起止时间；

步骤1.4，在全区建立不小于100m×100m的层位解释网格，解释线尽可能过基准井，以保证分层结果与井震联合标定结果相匹配；对解释网格上的每个地震道提取极值和过零点等特征值，包括极大值与极小值，正过零点和负过零点，并将所述特征值按时间依次排序，分别形成各井区的目标对比道集；存储每个地震道中所有特征点的Inline和Xline号以及自激自收时间T₀和振幅等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₀；

步骤1.7，将每个对比道中各反射层的T₀时间，输入地震解释软件，进行异常值剔除，平滑，完成粗网格层位三维闭合解释；进一步在整个工作区做平面插值，建立全区的地震等时地层格架。

所述各对比道中与各反射层的对应的自激自收时间T₀通过如下公式求得：

（2）

式中，i为Inline号，j为Xline号，n为层位号，B_m为工作小区号，BW_m为基准井号；T₀(i，j，n)为i,j道n反射层的自激自收时间；T₀(B_m，BW_m，n)为B_m工作小区内BW_m模型道上n反射层对应的T₀时间；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₀和对应的振幅；

步骤2.3，在m台计算机上，k个worker同时对n个分层体的每道数据按Inline方向做相邻道间DTW计算；

步骤2.4，每个分层体的DTW属性计算结果，按分层体起止自激自收时间T₀进行排序，起止时间外的数据全部置零，使每道的采样长度与原地震数据的采样长度相同；

步骤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₁，B₂，……Bm，标记其起止inline号和Xline号；每个区域中优选一口反射层位完整的井，作为基准井，记作BW₁，BW₂，……，BW_m；

步骤1.3，在基准井处选取井旁道，提取极值和过零点等特征值，包括极大值与极小值，正过零点和负过零点，并将所述特征值按时间依次排序，建立各井区的模型道，记作GM_w1，GM_w2， ……>GM_wm，存储每个特征点的时间、振幅参数，样点数控制在200个以内；存储每个模型道中各反射层T₁,T₂,>Tn对应的自激自收时间T₀，作为基准井中每个地质分层体的起止时间；

本发明是面向海量地震数据设计的，如果整个工作区的总地震道数如果小于100万道，在步骤1.2和步骤1.3中可不进行工作小区的划分，仅进行分层体的分割。

DTW属性体可以是DTW偏移时间，也可以是DTW距离。

步骤1.4，在全区建立不小于100m×100m的层位解释网格，解释线尽可能过基准井，以保证分层结果与井震联合标定结果相匹配；如图3所示；对解释网格上的每个地震道提取极值和过零点等特征值，包括极大值与极小值，正过零点和负过零点，并将所述特征值按时间依次排序，分别形成各井区的目标对比道集；存储每个地震道中所有特征点的Inline和Xline号以及自激自收时间T₀和振幅等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₀；所述各对比道中与各反射层的对应的自激自收时间T₀通过如下公式求得：

（2）

式中，i为Inline号，j为Xline号，n为层位号，B_m为工作小区号，BW_m为基准井号；T₀(i，j，n)为i,j道n反射层的自激自收时间；T₀(B_m，BW_m，n)为B_m工作小区内BW_m模型道上n反射层对应的T₀时间；T_dtw(i,j,n)为Bm工作小区内i，j道n反射层对应的动态弯曲匹配时间。

步骤1.7，将每个对比道中各反射层的T₀时间，输入地震解释软件，进行异常值剔除，平滑，完成粗网格层位三维闭合解释；进一步在整个工作区做平面插值，建立全区的地震等时地层格架。

步骤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₀和对应的振幅；

步骤2.3，在m台计算机上，k个worker同时对n个分层体的每道数据按Inline方向做相邻道间DTW计算；

步骤2.4，每个分层体的DTW属性计算结果，按分层体起止自激自收时间T₀进行排序，起止时间外的数据全部置零，使每道的采样长度与原地震数据的采样长度相同；

步骤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.9km²；叠前时间偏移数据体中共包含678,935道，每道1000个样点数；在不做分层体的情况下，采用8节点并行计算，总计算耗时约32小时；将上述数据依等比例剖分成10个分层体，每道样点数约92-135个；同样采用8节点并行计算，总计算耗时20.3分钟；将数据体依等比例剖分成20个分层体，每道样点数约6-67个，同样采用8节点计算，计算总耗时仅4.2分钟，这一结果与前二者相比分别减少了约455倍和3.8倍；表明本发明对于计算效率的提高非常明显。

如图5所示，图中多条不同落差断层（5-20m）的三维空间属性特征清晰，可用种子点法快速实现断层面的三维自动追踪。

本发明不会限制于本文所示的实施例，而是要符合与本文所公开的原理和新颖性特点相一致的最宽范围。