基于多属性引导下的三维全层位追踪方法

1.本发明属于地震图像解释领域，特别涉及一种层位追踪技术。


背景技术：

2.随着社会的进步与发展，各行各业对石油，天然气等资源的需求不断增大，石油行业得到了高度的关注和迅速的发展，因此对石油和天然气相关地质勘探的要求也越来越高。虽然新能源的开发和利用越来越纯熟，人们对传统能源的依赖性可能会越来越小，但是在现阶段石油和天然气依然是人类赖以生存的最重要的能源。然而要进行能源开采，首先必须清楚地下的地质构造，因此地震解释在地质勘探过程中发挥了重要作用。往往地震解释工作都是基于勘探前期采集到的地震图像进行的，因此在地震解释工作中地震图像的处理将是影响地震解释工作的重要环节。油气藏的开采离不开地震解释等工作，地质勘探工作结果越精准，资源开采工作就越顺利。层位解释工作是地址勘探的重要组成部分，也是地震图像解释的基础工作。
3.传统层位追踪手段多数是在二维地震数据上进行的，为了更好的了解地震数据层位信息，现今多数要求三维的层位数据，基于此目标，本发明研究的即为三维层位追踪方法。dbscan(density-based spatial clustering of applications with noise)是一种基于密度的聚类方法，之前也有使用dbscan进行层位追踪的方法，但多数使用的是二维数据，使用三维数据的方法也是先进行二维空间上的层位追踪，再扩展到三维。本发明提出的方法是一种全新的基于dbscan的三维层位追踪算法，首先阐述了各向同性的三维层位追踪方法，其次在原算法的基础上考虑了实际地震数据的属性特征，加入了地震的倾角及断层属性，形成了各向异性的三维层位追踪算法，其中，地震倾角等属性的计算使用了gst算法。本发明在实际地震数据进行了具体实验，并与商业软件效果进行了对比，体现了本发明算法的优越性。
4.在国内研究层位自动追踪起步晚，并且没有形成较为系统的方法，多是参照国外现有的技术方法再加以改进，与国外相比具有较大差距。地震层位追踪所使用的方法主要有：线性插值法、相关算法、边缘检测、人工神经网络和模式识别等。
5.线性插值方法较为简单方便，1998年赵改善等在《基于旅行时线性插值的地震射线追踪算法》一文中将线性插值方法应用于地震射线追踪，并提出了反射波lti算法。2007年赵成喜等在《线性插值法地震剖面反射层位追踪》提出了通过基于目的层种子点选取一定大小的时窗和道窗为临区，采用线性插值的方法对相邻地震道进行识别，追踪地震层位。线性插值方法程序设计简单，操作方便，成功实现了地震层位在二维地质剖面上的追踪。
6.杨峰等人于2003年提出了基于相关性的层位追踪方法。此方法首先提取了数据实域段，即真实数据，进而提取其他域的特征向量，将这些特征向量与确定层位的特征向量进行匹配分析，求其相关系数并进行模糊均值处理，得到综合相关系数后，提取出综合低速最大的实域段，即层位数据。2008年温庆庆采用了相似的算法完成了在地震数据中的层位追踪。虽然局部效果良好，但整体层位面是有残缺的。
7.p.alberts等(2002)提出了一种基于人工神经网络层位追踪算法，主要将跨断层等不连续地质结构的层位追踪问题视为模式识别的问题，然后采用神经网络进行模式识别。
8.深度学习方法需要算法收敛并且计算复杂，2020年，yang等人使用了cnn卷积神经网络进行同相轴的识别与追踪，但也是基于二维图像的实验。类似的，gao等人采用图像检测方法进行二维地震数据同相轴的追踪。
9.现有的层位追踪方法主要存在以下问题：
10.第一是这些方法并未考虑到实际地震数据的物理特征，如倾角、方位角等，不同的地震数据含有的地震特征也是千差万别，利用好这些属性能更好的推进层位追踪工作；
11.第二是这些方法多数是对三维地震数据进行了二维切片，在二维矩阵上进行基于深度学习或图像检测的聚类分析任务，有一定的局限性。


技术实现要素：

12.为解决上述技术问题，本发明提出一种基于多属性引导下的三维全层位追踪方法，结合实际地震数据具有的倾角属性，则可以在原本dbscan算法中加入倾角属性，使minpts邻域范围的选择更加灵活；而且地震数据的断层属性也是重要的特征，将其融合进算法中可实现更好的层位追踪效果。
13.本发明采用的技术方案为：基于多属性引导下的三维全层位追踪方法，包括：
14.s1、基于采样的三维地震波的振幅值提取种子点；
15.s2、遍历步骤s1提取的所有种子点；
16.s3、对于未标记的种子点，根据距离阈值、密度邻域、倾角属性、断层属性，选取相关种子点；
17.s4、根据步骤s3选取的未标记种子点的相关种子点，提取可信种子点；
18.s5、根据该未标记种子点的可信种子点以及断层属性，判断该未标记种子点为核心种子点或边缘种子点；
19.s6、若为核心种子点，则核心种子点对其可信种子点进行层位标记；
20.s7、若为边缘种子点，且该边缘种子点的可信种子点中包括某个核心种子点，则根据这个核心种子点对该边缘种子点进行层位标记；
21.s8、直至遍历完所有未标记的种子点，完成层位追踪。
22.本发明的有益效果：本发明着力于三维层位追踪算法，本发明如传统方法一样，使用地震数据的极值点为种子点，考虑到dbscan密度聚类方法是一个比较成熟且系统的聚类手段，可以以此为基础进行种子点的聚类，二维的聚类方法需要加以改进形成适合地震数据的三维层位聚类方法，从而形成了t-dhtd(three-dimensional horizon tracking based on dbscan)算法。在t-dhtd算法的基础上，因考虑到地震数据的物理特性，则在聚类过程扩散邻域的范围上做出了一些调整，二维dbscan聚类算法的minpts邻域是无方向扩散的，结合实际地震数据具有的倾角属性，则可以在原本dbscan算法中加入倾角属性，使minpts邻域范围的选择更加灵活；而且地震数据的断层属性也是重要的特征，将其融合进算法中可实现更好的层位追踪效果。
附图说明
23.图1dbscan聚类示意图；
24.图2为三维层位追踪流程图；
25.图3为地震道振幅图；
26.图4为极值点分布图；
27.图5为三维层位聚类算法流程图；
28.图6为地震数据的背景图；
29.图7为震数据slope属性图；
30.图8为层结构图；
31.图9为角属性面示意图；
32.图10为引导种子点扩散图；
33.图11为断层处种子点扩散示意图；
34.图12为断层属性示意图；
35.图13为本发明的属性引导算法流程图。
36.图14为层位追踪inline切面效果对比图；
37.其中，(a)为基于dbscan的层位追踪算法效果图(t-dhtd)；(b)为商业软件paleoscan效果图；(c)为基于dbscan的属性引导层位追踪算法效果图(t-dadht)。
38.图15为层位追踪xline切面效果对比图；
39.其中，(a)为基于dbscan的层位追踪算法效果图(t-dhtd)；(b)为商业软件paleoscan效果图；(c)为基于dbscan的属性引导层位追踪算法效果图(t-dadht)。
40.图16为层位追踪三维效果对比图；
41.其中，(a)为t-dhtd层位1；(b)为商业软件paleoscan层位1；(c)为t-dadht层位1；(c)为t-dhtd层位2；(d)为商业软件paleoscan层位2与层位3；(e)为t-dadht层位2。
具体实施方式
42.为便于本领域技术人员理解本发明的技术内容，下面结合附图对本发明内容进一步阐释。
43.伴随着现代社会科技技术日新月异的发展，人类对能源的需求也日益增长，新能源技术尚未成型，以石油煤矿为首的传统能源仍是主力，这些资源的开采都需要人类对地质情况有详细的了解。其中层位追踪是必不可少的技术，初始采集的叠前地震数据通过有序的排放形成本发明所需的三维地震数据，提取三维地震数据的同相轴可以帮助我们更好的了解目标地区的地质层位情况，普通层位追踪方法多数是基于二维地震数据的，难以对数据的全貌有足够的了解，本发明专注于三维地震数据的层位追踪方法。
44.本发明采用了以下三种技术：
45.1、dbscan算法
46.dbscan(density-based spatial clustering of applications with noise)算法是一种基于密度聚类的算法，顾名思义，本算法与其他聚类方法的区别即其他聚类算法都是基于距离的，基于密度聚类的算法对数据集的密度较为敏感，不被数据集的其它因素，如距离形状等所限制。该算法不受先验信息的影响，最终聚类团的数量是不定的，比较适合
处理新数据集。dbscan将数据集内元素分为三类，核心点，边缘点，噪声点，主要流程为：首先任取一元素p，在p的距离范围minpts内找到密度可达的其他元素点，进而形成簇团，若p是边缘点或噪声点则选取另一个元素点，重复直至聚类完毕，聚类如图1示，其中位于圆圈中心的点点为核心元素，箭头为扩散链路，圆圈内没有位于中心的点为边缘点，没有向外扩散的能力，圆圈外的点是噪声点，不被归为任何一个簇团。
47.本算法将dbscan用于三维地震数据的层位追踪中，形成全新的t-dhtd算法。经处理的地震数据是三维地震波的振幅值，沿z轴方向的一列振幅值数据可看为一个地震道，通过地震波的振幅可计算每个地震道的波峰波谷值，这些点即为本算法所用种子点，通过本发明的算法可形成三维层位。具体算法流程如图2示。
48.本发明算法流程具体为：首先是数据预处理，因为采集而来的地震数据多数含有噪声，造成同相轴不连续或其他情况，一般可采用一些滤波手段进行处理，如结构导向滤波等，可以强化地震数据的关联性。
49.种子点的提取是基于采样的数据振幅值的，前文所述了地震道的物理含义，将地震道的振幅转换为波形图，如图3示，其中向右凸出的点代表极大值，向左凸出的点代表了极小值，地震三维数据体中的数据值就是振幅值，将一个地震道看做二维向量x＝[x0,x1,
…
,xn]，因x内为离散点，无法具体求导，可以对其进行二次差分来近似二阶求导，进而求解极值，根据x计算一阶差分向量y，根据y计算二阶差分向量z，如式(1)、(2)所示：
[0050]
y＝[y0,y1,
…
,y
n-1
]＝[x
1-x0,x
2-x1,
…
,x
n-x
n-1
]
ꢀꢀ
(1)
[0051]
z＝[z0,z1,
…
,z
n-2
]＝[y
1-y0,y
2-y1,
…
,y
n-1-y
n-2
]
ꢀꢀ
(2)
[0052]
其中，x代表一个地震道数据，由多个数据点构成，x0x1xn代表x中的某个数据点。
[0053]
通过二阶导与极值的数学特性可知，z中大于0的元素所对应的x中的元素是极小值，小于零的是极大值。
[0054]
得到所有极值点后，将其视为种子点进行三维层位聚类，其中种子点在三维空间分布结果如图4所示，横纵坐标为crossline，inline，垂直方向为timeline，这是极值点进行抽稀后的显示效果，背景板为三维地震数据的三个横截面。
[0055]
三维层位聚类算法是整个流程的核心，算法类似于二维的dbscan算法，主要有三个参数：距离阈值(distance threshold，简称dthre)，密度邻域(density epsilon，简称denep)，最少点数(minimum points，简称minpts)；围绕着种子点形成了五个核心概念：
[0056]
(1)核心种子点：对应dbscan的核心点；若某个种子点在denep密度邻域与限定的dthre距离范围内包含大于等于minpts个数的可信种子点，则视其为核心种子点，且有向其可信种子点传播层位标记的权利；minpts一般取值为6-8的整数；
[0057]
(2)相关种子点：核心种子点邻域范围内的种子点，被称为核心种子点的相关种子点；这里所述的邻域范围应理解为denep密度邻域与限定的dthre距离范围；
[0058]
(3)可信种子点：若两个种子点互为相关种子点且两点都未进行层位标记，或两点层位标记相同，则两个种子点互为可信种子点。
[0059]
(4)边缘种子点：对应dbscan的边缘点；若某个种子点在denep密度邻域与限定的dthre距离范围内包含可信种子点个数小于minpts，但可信种子点中包含核心种子点，则视其为边缘种子点，边缘种子点有被可信种子点中所包含的这个核心种子点标记层位的权利，但无向denep邻域范围内的其它可信种子点传播自身层位的权利；
[0060]
(5)无关种子点：对应dbscan的噪声点；若某个种子点在denep密度邻域与限定的dthre距离范围内包含可信种子点个数小于minpts，且此范围中无核心种子点，则视其为无关种子点，它没有被核心种子点标记层位的权利，所以无关种子点无层位标记，则也无法向denep邻域范围的可信种子点传播层位标记；
[0061]
本发明算法步骤如下所示：
[0062]
2、基于梯度结构张量的倾角属性
[0063]
首先构造三维数据三个方向梯度值g1，g2，g3，然后构造gst矩阵：
[0064][0065]
u是原始地震数据，xyz分别代表了数据的维度分量，《》符号代表乘，t(x)即为和原地震数据大小相同的三维张量矩阵。
[0066]
再进行矩阵分解，特征值与特征向量提取预处理后，可得到本发明所需的地震slope属性体。公式(3)中，t(x)还可表示为式(4)：
[0067][0068]
t(x)是地震数据体，λu，λv，λw是t(x)的特征值，u，v，w是对应的特征向量。其中u是梯度变化最大的方向，w是梯度变化最小的方向。
[0069]
地震数据的背景图可以看到比较明显的倾向性，如图6示的某工区背景图，此切面图横坐标为crossline，纵坐标为timeline，为了保持图例的物理特性，本发明所有地震数据示意图的零点都在左上角的位置。从图6中可以看到多数同相轴都有一定的倾向性，我们称其为倾角或slope属性。前文所述中，u是地震数据梯度最速下降方向，则可将其沿三个维度分解为uc，ui以及u
t
，分别代表了u沿crossline，inline以及timeline方向的分量，则crossline的slope属性体sc以及timeline的slope属性体si求解如公式(5)、(6)所示：
[0070][0071][0072]
以图6背景图，绘制其种子点的slope属性示意图，如图7示，为了保证图例的清晰度，slope属性进行了抽样处理。图7的slope属性体为三维属性体si的部分切片，可以看出slope属性能够契合背景图中同相轴的走向。合理利用slope属性体可以进一步完善前文的层位追踪算法。
[0073]
3、基于梯度结构张量的断层属性
[0074]
地震数据有多种多样的结构，最为常见的复杂结构就为断层结构，根据gst分解得
到的三个特征值，可以将地震数据结构形象的表示出来。λu，λv，λw三个特征值不同的计算组合可以得到三个关系：线状置信度，面状置信度以及断层置信度：
[0075][0076]
三个特征值与置信度之间的关系如表1所示：
[0077]
表1特征值与置信度及地层结构关系
[0078][0079]
其中面状结构代表同相轴较为连续，地层水平方向的特征结构较为相似。线状结构不在水平方向结构相似，则可能为断层体，各项同性代表地层结构研所有方向都一样，地层分层不明显。本发明使用线状结构来判断断层，提取断层属性体f：
[0080][0081]
若f越接近于0则越可能是断层，越接近于1则代表同相轴越连续。
[0082]
本发明结合上述的三种技术：基于dbscan的层位追踪算法，以及slope倾角与断层属性，进行层位追踪，具体实现过程如图13所示。
[0083]
地震结构是存在倾斜度及倾角属性的，种子点的扩散不能仅使用距离阈值dthre进行可信种子点判断，要在此基础上加入倾角的限制，上述slope属性体有两个，分别为沿crossline方向与沿inline方向的属性体。从空间几何知识可知，三维空间内的两个矢量可以确定一个平面，则可通过种子点的两个方向的slope值确定种子点在地震数据的所有方向的slope倾角走向，如图9示，通过p点的沿crossline方向的slope属性sc以及沿timeline的slope属性si可以在三维空间中确定p点的slope平面。
[0084]
算法实现过程中，加入了slope倾角属性体的种子点扩散过程如图10，三个地震道分别为s1，s0，s5，可以看到p点上下黑色实线所示，大小并没有发生改变，但范围有所不同，整个距离阈值偏离了水平线，水平线与斜线之间的高度差为slope值，即图10中slope所指的虚线所示。若不使用slope倾角值来改变阈值范围，则可信种子点的选取范围只能如p点上下的水平两道虚线所示，采用人工观测的方法可以看出，在水平两道虚线的范围内，p点无法在s1与s5地震道中找到合适的相关种子点，从算法角度要结束验证p点为核心种子点的后续步骤，但从地震数据的物理特性角度出发，若改变距离阈值，p点能在s1与s5上找到波形特征相似的波峰值，且在黑色实线的范围内，则p点经后续算法判断可能为核心中心种
子点。所以，在种子点的扩散算法中，加入种子点的物理特征，可以更好的找到匹配的相关种子点，算法效果也应有进一步的提升。
[0085]
实际地震数据中，倾角特征表现的最为明显，其次为断层属性，图6的背景中，可以看到一些明显的断层痕迹，在这些地方，断层两端的种子点相差较远，无论如何改变距离阈值dthre值都无法将两个种子点匹配为相关种子点，如图11，若按照距离阈值及slope属性引导来选取点p的相关种子点则大概率会选中种子点q2，观测波形可知，p点上下的波形与q2点并不相似，而与q1点较为相似，这是一个比较简单的断层体，则在此处的算法步骤中会出现错误的层位标记传播结果。
[0086]
基于实际数据的考虑，本算法加入了新的约束条件——断层属性体f，算法新加入步骤：检测种子点的断层属性f值的大小，若高于断层阈值f
th
，则加入相关性检测算法。若初始种子点与相关种子点的相关系数低于设定阈值，则否定相关种子点成为可信种子点的资格。如图11中，p点处于断层区域，检测到q2为相关种子点，在此判断p点与q2点的相关系数，低于阈值则不认定q2点为p点的可信种子点。f
th
依据具体数据的断层属性体的值进行取值，一般取值为0.5-0.7。
[0087]
图12的断层属性对应地震数据背景图6，深色所示区域皆为断层属性值较高的区域，颜色越深，f值越小，则越有可能为断层。浅色则代表同相轴较为连续。通过图12可以看出地震数据的断层区域有很大的占比。基于slope倾角属性与f断层属性，可以设计新的基于dbscan的属性引导层位追踪算法。
[0088]
图13即为本发明的基于dbscan的属性引导层位追踪算法。理论分析阶段，本发明的算法在三维空间中实现了层位聚类，且结合了实际地震数据的特性，抽取地震数据的倾角属性与断层属性，更加符合实际数据特征。如图13所示搜索队列为保存选定的未标记核心种子点及其所有可信种子点的一个队列。
[0089]
本发明选取真实地震数据进行实验验证。通过不同的对比算法对工区数据进行三维层位追踪，对效果进行全方位的对比。本节实验效果来自三种算法。算法一：基于dbscan的层位追踪算法；算法二：商业软件的层位追踪效果；算法三：本发明的基于dbscan的属性引导层位追踪算法。
[0090]
实验参数设置如下：距离参数dthre为2.1，可信种子点minpts个数为8，邻域范围denep为1，即扩散范围为种子点周围8个地震道。
[0091]
选取大庆工区的数据集在时间轴上取301个采样点，crossline与inline为251个采样点，数据为一个301*251*251的三维采样点矩阵。本工区将对比基于dbscan的层位追踪算法与基于dbscan的属性引导层位追踪算法效果。实现效果将从三个方面进行对比，分别为inline切面二维图，crossline切面二维图，三维立体层位效果图。
[0092]
如图14所示，此为最终层位效果图的inline切面图层位抽稀结果，选取inline等于200的二维切面进行效果图对比。
[0093]
其中比较明显的差别在椭圆圈标记处，其中，在图14左上方的1号椭圆圈处，图14(a)中1号椭圆圈内部第一条层位不够完整，缺少了一小部分；图14(b)中1号椭圆圈内部第一条层位只有一小段，缺少了大片的层位；只有图14(c)中1号椭圆圈内部第一、二条层位都比较完整。
[0094]
在图14左中部的2号椭圆圈标记处，如图14(b)所示的层位在中间断为两半，说明
在商业软件的追踪过程不够完善，而图14(a)与分图14(c)中2号椭圆圈处的层位是紧密连接在一起的，效果良好。
[0095]
图14中右下方的3号椭圆圈处的部分层位的追踪效果都不一样，如图14(a)的走势较为平缓，图14(b)的走势极为陡峭，图14(c)的走势较为折中。综合分析，图14(c)的效果更符合背景图的同相轴走势。
[0096]
结合图14的三张分图来看，除了圆圈部分的层位效果，其他层位走势大致相同，说明本发明的算法比较符合实际地质层位分布，通过上文分析，t-dadht(3d attribute guided dbscan horizon tracking)算法的层位效果更为完整。
[0097]
图15是追踪结果的xline切面图，其中效果差距较为突出的地方同样用不同标号椭圆圈标记了出来。图15(a)与(c)的1号椭圆圈处的第一条层位在分图(b)中并不存在，结合图14分析，商业软件在该层位的追踪过程中结束的过于突兀，小部分的层位碎片在图14(b)中仍有显示，但在图15(b)中找寻不到。
[0098]
图15中2号椭圆圈处的第一条层位情况与1号椭圆圈处的第一条层位较为相似，本发明的t-dhtd算法与t-dadht算法追踪的2号椭圆圈处的第一条层位出现在了图15(c)所示的左半部分，商业软件只在图15(b)所示的右半部分有少部分的2号椭圆圈处的第一条层位显示。
[0099]
图15中3号椭圆圈处均有层位。由于商业软件在图14(b)2号椭圆圈处的追踪出现了断层，在图15(b)中显示的层位结果即为图14(b)断层左侧的层位，说明商业软件的追踪层位结果有少量的缺失。
[0100]
在图15(a)中4号椭圆圈处的层位出现了窜层，而其他两分图并无此种情况，说明t-dhtd算法确实并不完善。结合图14与15，t-dadht算法效果更为优秀。
[0101]
图14(b)中2号椭圆圈断层左侧的层位标记为层位3，三种方法的三维层位显示效果如图16所示。
[0102]
结合图14与15，对比图16的(a)、(b)、(c)分图，可以看出初始t-dhtd算法的三维层位在crossline位于[0,50]，inline位于[200,250]的范围内有部分的缺失；而商业软件只追出了少部分的层位就截止了，如图16(b)所示。而图16(c)中t-dadht算法的层位效果比较完整。
[0103]
图16(d)(e)(f)是同一层位的对比效果，图16(d)、(f)的效果比价类似，说明地质情况良好，较少受到地震数据属性特征的影响。而图16(e)中有两种层位，标记的1号层位缺失的一角由2号层位补齐，还有少量的层位缺失情况，这与二维切面图14(b)与图15(b)的层位效果符合。
[0104]
本领域的普通技术人员将会意识到，这里所述的实施例是为了帮助读者理解本发明的原理，应被理解为本发明的保护范围并不局限于这样的特别陈述和实施例。对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的权利要求范围之内。