CN110749929B - 一种基于机器学习的复杂地区地震反射层识别与追踪方法 - Google Patents

Abstract

本发明公开了一种基于机器学习的复杂地区地震反射层识别与追踪方法，其包括步骤：I.确定需要追踪解释的某一地震反射层位；II.提取目标层位处的局部地震波形数据作为标准地震道；利用动态时间规整算法对标准地震道和各个地震道滑动时窗内的波形数据进行匹配，获得最短弯曲路径，将最短弯曲路径值赋值给滑动时窗中心位置；原始地震剖面中所有滑动时窗计算完成后，计算获得相似性地震属性剖面；III.确定地震剖面层位追踪的纵向深度范围；IV.以相似性地震属性剖面作为输入，以追踪范围作为约束条件，利用Viterbi算法进行约束优化处理，对目标层位进行识别与追踪。本发明能够实现复杂地区地震反射层位识别与追踪，有利于提高地震反射层位解释追踪的精度与效率。

Description

一种基于机器学习的复杂地区地震反射层识别与追踪方法

技术领域

本发明涉及一种基于机器学习的复杂地区地震反射层识别与追踪方法。

背景技术

地震层位的追踪与对比是地震资料解释的一项基本工作，是地震构造解释及相关反演的重要基础。近年来，随着复杂地区地震勘探程度加大，地震勘探采集到的数据日益增多。传统的人工层位追踪方法不仅存在效率低和受人为因素影响显著等问题，而且追踪精度低，地层对比分析难度大，影响了地震勘探效率。因此，有必要研发一种适用于复杂地区稳定高效的地震反射层识别、追踪与对比分析方法，这对地震资料的解释、反演具有重要的意义。

对此，国内外学者曾进行过大量研究，目前常用的地震层位自动追踪方法主要有模糊聚类、AR自动追踪技术、非监督学习的自组织神经网络、遗传算法、边缘检测、地震DNA算法等，这些方法在实际应用中均已取得了一定的预测效果，但是这些方法通常存在稳定性弱、适应性差等问题，对于复杂地区的地震资料，当层位形态较为复杂时，会影响到层位追踪的精度与效率。机器学习作为人工智能技术的重要组成部分，因其具有特征提取明显、鲁棒性强、泛化能力强、效率高等特点，在语音识别、图像处理、行为检测等领域得到了广泛的应用。机器学习方法在其他领域的快速发展极大促进了在地球物理勘探领域的应用。

发明内容

本发明的目的在于提出一种基于机器学习的复杂地区地震反射层识别与追踪方法，以提高复杂地区地震发射层位追踪的精度与效率，便于地层的对比与分析。

本发明为了实现上述目的，采用如下技术方案：

一种基于机器学习的复杂地区地震反射层识别与追踪方法，包括如下步骤：

I.输入原始地震剖面，确定需要追踪解释的某一地震反射层位，即目标层位；

II.提取目标层位处的局部地震波形数据作为标准地震道；

标准地震道的生成过程如下：

通过钻井层位标定，提取以目标层位中点O为中心，长度为n的地震道以及与该地震道邻近的多个地震道的地震波形数据作为标准地震道X，如公式(1)所示：

X＝(x₁,x₂,…,x_n) (1)

其中，i＝1,2,…,n，n为标准地震道的长度；标准地震道X中第i个数据点的值为：

其中，l＝1,2,…,F，F为邻近的地震道数；x_l(i)表示第l个邻近地震道第i个数据点的值；

利用动态时间规整算法对标准地震道和各个地震道滑动时窗内的波形数据进行匹配，获得最短弯曲路径，最短弯曲路径的计算过程如下：

设定输入的原始地震剖面为M×N的矩阵I，则：

其中，p＝1,2,…,M，q＝1,2,…,N，M为采样点数，N为地震道数；

设定滑动时窗内的波形数据的长度与标准地震道的长度相同且均为n，每次移动的步长为1，则第q道某个滑动时窗内的波形数据Y为：

Y＝(y₁,y₂,…,y_n) (4)

其中，j＝1,2,…,n；下面对标准地震道X＝(x₁,x₂,…,x_n)与第q道某个滑动时窗内的波形数据Y＝(y₁,y₂,…,y_n)进行波形动态时间规整，即：

定义距离函数d_ij，该距离函数d_ij采用欧拉距离，表达式如下：

根据X、Y内不同数据之间点与点的欧拉距离，构成n×n的动态时间弯曲距离矩阵d：

在矩阵d中，相邻矩阵元素的集合构成一条弯曲路径，如公式(7)所示：

W＝{w₁,…,w_k,…,w_K} (7)

式中，n≤K≤2n-1；w_k＝(d_ij)_k，表示W中第k个元素x_i与y_j两点的距离；

从求解得到的弯曲路径中选取最短弯曲路径W_best，表达式如公式(8)所示：

为求解公式(8)，由递归思想，构造累积代价矩阵γ，累积距离γ(i,j)为当前格点距离d(i,j)与能够到达该点的最小的邻近元素的累积距离之和：即：

其中，i＝1,2,…,n，j＝1,2,…,n,n为标准地震道与滑动时窗内的波形数据的长度；

因此，X和Y的动态时间规整距离为：DTW(X,Y)＝γ(n,n) (10)

经过上述计算，获得标准地震道X与滑动时窗内波形数据Y的最佳匹配关系，计算出最短弯曲路径W_best，将最短弯曲路径值赋值给滑动时窗中心位置；

原始地震剖面中所有滑动时窗计算完成后，得到最短路径距离矩阵P，公式如下：

其中，最短路径距离矩阵P为M×N的矩阵；

W_pq表示最短路径距离，W_pq值越大表示相似性越差，W_pq值越小表示相似性越高；

对最短路径距离矩阵P进行如下处理：

式中，W_max表示矩阵P中的最大值，W_min表示矩阵P中的最小值；

Q_pq值越大表示相似度越高；经过公式(12)处理后获得相似性地震属性剖面矩阵Q：

III.根据所要追踪解释的目标层位位置，确定地震剖面层位追踪的纵向深度范围；

设定目标层位的顶部界面深度为h₁，底界面深度为h₂，则定义如下的追踪范围：

h+h₁＜σ＜h+h₂ (14)

其中，σ为层位追踪搜寻的纵向深度范围，h为调节深度；

IV.对将步骤II得到的相似性地震属性剖面矩阵作为输入，以步骤III获得的追踪范围作为约束条件，利用Viterbi算法进行约束优化处理，对目标层位进行识别与追踪；

具体过程如下：

IV.1.地震反射特征最大相似度做向前积分，并记录累加路径，公式如下：

其中，ω(p,q)为步骤II获得的相似性地震属性剖面矩阵Q；

表示以(ik,q-1)为坐标，ik为取值范围，取q-1之前累积的结果中最大值

该值与下一个相似性地震属性剖面矩阵ω(p,q)求和即得ω′(p,q)；

L(p,q)＝ik|_q-1表示求和累积的过程路径；

s为Viterbi算法追踪的时窗大小；

IV.2.搜寻地震反射特征累加路径最大相似度值，向后递归搜寻最短路径，追踪地震反射层位；

当公式(15)的累积达到q＝N-1时，搜寻到最大值ω′_max(p,N)，并记下所对应的路径L(p,N-1)＝ik|_N-2为x_ipick|_N-1＝p；

其中，L(p,N-1)表示公式(15)的累积达到q＝N-1，搜寻到最大值ω′_max(p,N)时所对应的路径位置，x_ipick|_N-1表示拾取路径的纵坐标；

根据向后递减跟踪原理搜寻最短路径并计算所需要的值v(q′)，得如式：

其中，q′＝N-1,N-2,…,2,1；x_ipick|_q′＝L(x_ipick|_q′+1,q′)表示所选择累积积分最大相似度值所对应路径中的每个点坐标，即为层位追踪的路径；v(q′)为最短路径中每个点的坐标对应的相似性属性值，ω(x_ipick|_q′)表示最短路径对应的相似性属性值。

本发明具有如下优点：

如上所述，本发明提供了一种基于机器学习的复杂地区地震反射层识别与追踪方法，该方法是基于动态时间弯曲算法以及最短路径积分法(Viterbi算法)等机器学习法，并结合地质构造与地震数据特点而提出的，该方法自动追踪的层位连续性好、精度高、稳定性强，在一定程度上提高了地震层位解释的精度与效率，适用于复杂地区地震反射层位高效率追踪和地震资料的高精度构造解释。

附图说明

图1为本发明实施例中基于机器学习的复杂地区地震反射层识别与追踪方法的流程图。

图2为本发明实施例中两序列所对应的某条弯曲路径示意图。

图3为本发明实施例中标准地震道提取的示意图。

图4为本发明实施例中某探区常规二维地震的剖面图。

图5为本发明实施例中时窗为3时的地震层位追踪结果示意图。

图6为本发明实施例中时窗不加限制时的地震层位追踪结果示意图。

图7为本发明实施例中时窗为12时的地震层位追踪结果示意图。

图8为本发明实施例中层位识别与追踪结果示意图。

具体实施方式

下面结合附图以及具体实施方式对本发明作进一步详细说明：

结合图1所示，基于机器学习的复杂地区地震反射层识别与追踪方法，包括如下步骤：

I.输入原始地震剖面，确定需要追踪解释的某一地震反射层位，即目标层位。

II.利用动态时间规整算法计算相似性地震属性剖面(矩阵)，以压制非目标区地震反射层特征，同时增强目标地震反射层特征，为Viterbi算法进行地震反射层位的识别与追踪提供高精度数据源，以增加层位识别追踪的稳定性与准确性。

基于互相关的传统层位自动追踪技术，是通过直接计算地震道与地震道之间的相关性进行同相轴的自动追踪，该方法计算简单，对于构造简单地区层位追踪效果较好。

但是，对于构造复杂，横向变化快的区域(如断层的阻隔)，同一层位在复杂地区的起伏变化，就会造成波形的拉伸及压缩(同一层位的地震波形特征仍具有一定的相似性)，波形变化大，往往无法进行地震反射层位的连续追踪，自动识别追踪的结果连续性差。

同时，传统的层位追踪方法容易出现串层现象，不能按照正确的目标层位进行追踪，层位追踪的精度较低。因此，传统的相关计算方法在构造复杂地区的层位追踪难度大，精度低。

针对传统层位追踪存在的问题(复杂地区横向变化快、串层现象)，本发明提出如下的解决方案：首先引入动态时间弯曲算法，对标准地震道波形与各道滑动时窗内的局部波形数据进行拉伸和压缩，获得一个最短的弯曲路径。然后，将最短弯曲路径值赋值给滑动时窗中心位置，原始地震剖面中所有滑动时窗计算完成后，计算获得相似性地震属性剖面。通过该过程以增强目标层的地震反射特征，增加目标层位地震反射特征的连续性，压制非目标层地震反射特征，有效解决复杂地区波形横向变化快，导致层位追踪不连续、串层现象等，为后续步骤利用Viterbi算法进行地震反射层位的识别与追踪提供高精度数据源，增加层位识别追踪的稳定性与准确性。本实施例计算相似性地震属性的基本原理如下：

首先利用波形的相似性原理，提取目标层位处的局部地震波形数据作为标准地震道。

再使用动态时间规整算法(Dynamic Time Warping，简称DTW)对各地震道滑动时窗内的波形数据进行匹配，以获得最短弯曲路径，将最短弯曲路径值赋值给滑动时窗中心位置；

原始地震剖面中所有滑动时窗计算完成后，经处理获得相似性地震属性剖面，以增强目标地震反射层特征，增加目标层位地震反射特征的连续性，压制非目标层地震反射特征。

标准地震道的生成过程如下：

通过钻井层位标定，假设所要追踪的层位如图3中的黑色线条所示。

提取以目标层位中点O为中心，长度为n的地震道以及与该地震道邻近的多个地震道的地震波形数据作为标准地震道X，如公式(1)所示：

X＝(x₁,x₂,…,x_n) (1)

其中，i＝1,2,…,n，n为标准地震道的长度，一般为奇数。

由于低信噪比或者噪音的干扰，获得的地震数据难免会存在局部异常值，为增加层位追踪的精度，在标准地震道提取过程中，通常要融合井孔地震道以及相邻地震道的信息，以降低局部异常值的影响。标准地震道X中第i个数据点的值为：

其中，l＝1,2,…,F，F为邻近的地震道数；x_l(i)表示第l个邻近地震道第i个数据点的值。

本实施例利用动态时间规整算法(Dynamic Time Warping，简称DTW)对标准地震道和各个地震道滑动时窗内的波形数据进行匹配，以获得最短弯曲路径。

动态时间规整算法由于具有很好的鲁棒性，较高的匹配精度，通常被应用在语音识别等模式识别领域，它能够将时间规整和距离测度两者结合起来，是一种非线性规整技术。

动态时间规整DTW是一个典型的优化问题，它用满足一定条件的时间规整函数描述测试模板和参考模板的时间对应关系，通过对时间序列进行拉伸和压缩，求解两模板匹配时累计距离最小所对应的规整函数，来计算两个时间序列性之间的相似性。

动态时间规整DTW算法是利用动态规划原理，将问题中的全局最优化问题分解为多个局部最优化问题，并进行逐步的解决。最短弯曲路径的计算过程如下：

设定输入的原始地震剖面为M×N的矩阵I，则：

其中，p＝1,2,…,M，q＝1,2,…,N，M为采样点数，N为地震道数。

设定滑动时窗内的波形数据的长度与标准地震道的长度相同且均为n，每次移动的步长为1，则第q道某个滑动时窗内的波形数据Y为：

Y＝(y₁,y₂,…,y_n) (4)

其中，j＝1,2,…,n；下面对标准地震道X＝(x₁,x₂,…,x_n)与第q道某个滑动时窗内的波形数据Y＝(y₁,y₂,…,y_n)进行波形动态时间规整。即：

定义距离函数d_ij，该距离函数d_ij采用欧拉距离，表达式如下：

根据X、Y内不同数据之间点与点的欧拉距离，构成n×n的动态时间弯曲距离矩阵d：

在矩阵d中，相邻矩阵元素的集合构成一条弯曲路径，如公式(7)所示：

W＝{w₁,…,w_k,…,w_K} (7)

式中，n≤K≤2n-1。

W中第k个元素的表达式为：w_k＝(d_ij)_k，表示W中第k个元素x_i与y_j两点的距离。两序列X和Y所计算出的某条弯曲路径图2所示，在计算和求解弯曲路径时，需要考虑弯曲路径的具体要求。弯曲路径需满足以下条件：

①有界性：弯曲路径W中的元素需满足：w₁＝d₁₁,w_k＝d_nn，这说明弯曲路径的起始点和终止点是确定的。

②连续性：所求解的弯曲路径W，需要满足连续性的要求。

③单调性：弯曲路径W需要满足单调性的要求。

④边界条件：假设弯曲路径W的表达式为W＝{w₁,…,w_k,…,w_K}，则弯曲路径下的元素总数K需满足n≤K≤2n-1。

由弯曲路径的必要条件可知，满足上述条件的弯曲路径不止一条。

从求解得到的弯曲路径中选取最短弯曲路径，表达式如公式(8)所示：

为求解公式(8)，由递归思想，构造累积代价矩阵γ，即：

其中，i＝1,2,…,n，j＝1,2,…,n,n为标准地震道与滑动时窗内的波形数据的长度。

因此，X和Y的动态时间规整距离为：DTW(X,Y)＝γ(n,n) (10)

经过上述计算，获得标准地震道X与滑动时窗内波形数据Y的最佳匹配关系，计算出最短弯曲路径W_best，将最短弯曲路径值赋值给滑动时窗中心位置。

原始地震剖面中所有滑动时窗计算完成后，获得最短弯曲路径矩阵P，公式如下：

其中，最短弯曲路径矩阵P为M×N的矩阵。

W_pq表示最短路径距离，W_pq值越大表示相似性越差，W_pq值越小表示相似性越高。

为了便于步骤IV利用Viterbi算法对层位进行识别与追踪，本发明实施例还对由公式(11)得到的最短弯曲路径矩阵P中的元素进行如下处理：

式中，W_max表示矩阵P中的最大值，W_min表示矩阵P中的最小值。

Q_pq值越大表示相似度越高；经过公式(12)处理后获得相似性地震属性剖面矩阵Q：

III.根据所要追踪解释的目标层位位置，确定地震剖面层位追踪的纵向深度范围。

为提高Viterbi算法进行层位识别追踪的准确性与效率，减少计算量，首先要确定层位追踪范围，根据所要解释地震反射层位的位置确定图3所示解释地震剖面层位的顶底界面。

设定目标层位的顶部界面深度为h₁，底界面深度为h₂，则定义如下的追踪范围：

h+h₁＜σ＜h+h₂ (14)

其中，σ为层位追踪搜寻的纵向深度范围，h为调节深度，h的取值不宜过大。

IV.对将经过步骤II得到的相似性地震属性剖面作为输入，以步骤III获得的追踪范围作为约束条件，利用Viterbi算法进行约束优化处理，对目标层位进行识别与追踪。

具体过程如下：

IV.1.地震反射特征最大相似度做向前积分，并记录累加路径，公式如下：

其中，ω(p,q)为步骤II获得的相似性地震属性剖面矩阵Q；

表示以(ik,q-1)为坐标，ik为取值范围，取q-1之前累积的结果中最大值

该值与下一个相似性地震属性剖面矩阵ω(p,q)求和即得ω′(p,q)；

L(p,q)＝ik|_q-1表示求和累积的过程路径；

s为Viterbi算法追踪的时窗大小，一般为奇数。

对公式(15)的物理意义阐述如下：对经过步骤II处理后获得的目标层位地震反射特征得到明显增强的相似性地震属性剖面作为Viterbi算法的输入，以步骤III获得的纵向深度范围作为约束条件，在全区范围内做最大相似度向前积分，并记录累加路径。

IV.2.搜寻地震反射特征累加路径最大相似度值，向后递归搜寻最短路径，追踪地震反射层位。

当公式(15)的累积达到q＝N-1时，搜寻到最大值ω′_max(p,N)，并记下所对应的路径L(p,N-1)＝ik|_N-2为x_ipick|_N-1＝p；

其中，L(p,N-1)表示公式(15)的累积达到q＝N-1，搜寻到最大值ω′_max(p,N)时所对应的路径位置，x_ipick|_N-1表示拾取路径的纵坐标。

根据向后递减跟踪原理搜寻最短路径并计算所需要的值v(q′)，得如式：

其中，q′＝N-1,N-2,…,2,1；x_ipick|_q′＝L(x_ipick|_q′+1,q′)表示所选择累积积分最大相似度值所对应路径中的每个点坐标，即为层位追踪的路径；v(q′)为最短路径中每个点的坐标对应的相似性属性值，ω(x_ipick|_q′)表示最短路径对应的相似性属性值。

对公式(16)的物理意义阐述如下：经过公式(15)计算获得全区范围内多条最大相似度积分路径后，搜寻地震反射特征累加路径最大的相似度值，并向后递归搜寻最短路径，并追踪地震反射层位；由于在步骤II对目标层位的反射特征进行了增强处理，在该过程中，搜寻的地震反射特征最大相似度累加路径，即为目标层位，以此完成层位的识别与追踪。

在Viterbi算法进行层位追踪的过程中，时窗大小s将会影响层位追踪的精度，时窗大小s过小，追踪路径易偏离正确的方向，时窗大小s过大，追踪精度将会降低。

因此，必须确定层位追踪最佳的时窗大小s，具体过程如下：

首先定义一初始窗口s后，进行层位追踪，若满足精度要求，则停止层位追踪；

若不满足精度要求，将时窗大小s修改为：

s＝s+λ (17)

在此基础上重新进行层位追踪，直至满足精度要求为止，以获得最佳的追踪窗口。

根据地质地震任务，对步骤IV获得的层位追踪结果进行对比分析，评估层位追踪效果。

本发明基于动态时间弯曲算法及最短路径积分法(Viterbi算法)等机器学习法，并结合地质构造及地震数据特点，提出一种适用于复杂地区地震反射层的识别与追踪方法，以提高复杂地区地震发射层位追踪的精度与效率，该方法中各个步骤的组合以及它们之间的特定的组合顺序共同构成了本发明对现有技术做出贡献的技术特征，缺一不可，具体分析如下：

步骤I首先输入原始地震剖面，确定需要追踪解释的某一地震反射层位。步骤II在步骤I确定目标层位的基础上，为解决传统层位拾取过程中，构造复杂、横向变化快的区域存在层位追踪连续性差、精度低以及串层等现象，本发明利用动态时间规整算法对标准地震道和各个地震道滑动时窗内的波形数据进行匹配，获得最短弯曲路径，将最短弯曲路径值赋值给滑动时窗中心位置；原始地震剖面中所有滑动时窗计算完成后，计算获得相似性地震属性剖面；通过该过程以增强目标层的地震反射特征，增加目标层位地震反射特征的连续性，压制非目标层地震反射特征，有效解决复杂地区波形横向变化快导致层位追踪不连续、串层现象等，为后续步骤利用Viterbi算法进行地震反射层位的识别与追踪提供高精度数据源，增加了层位识别追踪的稳定性与准确性。步骤III根据所要追踪解释的目标层位位置，确定地震剖面层位追踪的纵向深度范围，该过程能够提高Viterbi算法进行层位识别追踪的准确性与效率，减少计算量。步骤IV对经过步骤II处理后获得的目标层位地震反射特征得到明显增强的相似性地震属性剖面作为Viterbi算法的输入，以步骤III获得的纵向深度范围作为约束条件，在全区范围内做最大相似度向前积分，并记录累加路径，搜寻地震反射特征累加路径最大的相似度值，并向后递归搜寻最短路径，并追踪地震层位；由于在步骤II对目标层位的反射特征进行了增强处理，在该过程中，搜寻的地震反射特征最大相似度累加路径，即为目标层位，以此完成层位的识别与追踪。最后，根据地质地震任务，对步骤IV获得的层位追踪结果进行对比分析，评估层位追踪效果。

将本发明应用于某地区，以此来验证本发明方法的有效性与可行性。

根据区域资料及钻井数据分析，该探区探测目的层为海相地层，且存在一套煤系地层，即该区上二叠统主要为含煤地层，但由于该区域勘探程度低，钻井少，煤层的密度较围岩低，导致地震波频率降低，地震反射表现为强振幅，低频率。

由于断层原因，断层两侧(以cdp960道为界，如图4所示)地层是否同属一地层不易判断。为此，对该研究区采用本发明所提出的方法进行目标层位的识别、追踪与对比分析，对断层两侧地层的亲缘关系进行判断。

将本发明方法应用于该地区进行地震反射层位识别与追踪的具体步骤如下：

1.如图4所示是某地区经过一系列处理之后得到的原始地震剖面，目标层位由于断层阻隔(以cdp960道为界，如图4所示)，两侧的地层是否同属一地层不易判断。

因此，首先对目标层位的地震反射特征进行增强处理，选取cdp580井旁地震道二叠系(矩形框)的一段数据作为标准地震道，以整个地震剖面作为输入进行目标层位地震反射特征增强，得到目标层位地震反射特征增强的二维相似地震属性剖面图。

经过该步骤处理后，与标准地震道(二叠系)同为一时期的地层突出显示，且地层横向连续性好；与标准地震道波形不相似区域，波形被压制，压制效果明显。

2.经过上述步骤处理后，目标层位的反射特征得到了明显的突出，非目标层的地震反射特征得到了明显的压制，为利用Viterbi算法进行层位的识别与追踪提供高精度数据源，增加了层位追踪的稳定性与准确性。在Viterbi算法进行层位追踪的过程中，必须确定层位追踪最佳的时窗大小。本发明实施例在层位识别与追踪过程中，采用不同的Viterbi算子计算时窗，对目标层位进行识别与追踪，分析时窗大小对于地震层位搜寻精度的影响。

首先将Viterbi算子计算时窗限定为3，层位追踪的结果如图5所示。

从图5中能够看出Viterbi算子搜寻的路径并不适合实际层位的分布，搜寻的路径存在偏差，搜寻的结果并未达到理想的结果。

然后将Viterbi算子搜寻层位的时窗不加限制，那么在搜寻过程中，路径将不受约束地在整张剖面搜寻，进行运算之后的结果如图6所示。

从图6中能够看出，搜寻数值点遍布整张剖面，并未达到理想状态。

经历过一系列修改时窗大小的尝试之后，将Viterbi算子计算时窗修改为12，程序运行之后的结果如图7所示，准确地搜寻出了目标层的地震层位，达到了理想状态。

图8为采用本发明方法进行目标层位识别与追踪的结果示意图。

从图8中能够看出，层位追踪具有较好的稳定性，地层横向连续性好，断层两侧的地层，具有较强的相似性，因此，能够判断断层两侧地层是同属一地层。

与传统的地震反射层位追踪方法相比，本发明方法不仅能够实现复杂地区地震反射层位的识别与追踪，而且所追踪的层位连续性好，精度高，即本发明方法适应于复杂地区地震反射层位高效率识别与追踪和地震资料的高精度构造解释。

当然，以上说明仅仅为本发明的较佳实施例，本发明并不限于列举上述实施例，应当说明的是，任何熟悉本领域的技术人员在本说明书的教导下，所做出的所有等同替代、明显变形形式，均落在本说明书的实质范围之内，理应受到本发明的保护。

Claims (1)

1.一种基于机器学习的复杂地区地震反射层识别与追踪方法，其特征在于，

包括如下步骤：

I.输入原始地震剖面，确定需要追踪解释的某一地震反射层位，即目标层位；

II.提取目标层位处的局部地震波形数据作为标准地震道；

标准地震道的生成过程如下：

通过钻井层位标定，提取以目标层位中点O为中心，长度为n的地震道以及与该地震道邻近的多个地震道的地震波形数据作为标准地震道X，如公式(1)所示：

X＝(x₁,x₂,…,x_n) (1)

其中，标准地震道X中第i个数据点的值为：

其中，i＝1,2,…,n，n为标准地震道的长度，l＝1,2,…,F，F为邻近的地震道数；x_l(i)表示第l个邻近地震道第i个数据点的值；

利用动态时间规整算法对标准地震道和各个地震道滑动时窗内的波形数据进行匹配，获得最短弯曲路径，最短弯曲路径的计算过程如下：

设定输入的原始地震剖面为M×N的矩阵I，则：

其中，M为采样点数，N为地震道数；

设定滑动时窗内的波形数据的长度与标准地震道的长度相同且均为n，每次移动的步长为1，则第q道某个滑动时窗内的波形数据Y为：

Y＝(y₁,y₂,…,y_n) (4)

其中，q＝1,2,…,N，N为地震道数；下面对标准地震道X＝(x₁,x₂,…,x_n)与第q道某个滑动时窗内的波形数据Y＝(y₁,y₂,…,y_n)进行波形动态时间规整，即：

定义距离函数d_ij，该距离函数d_ij采用欧拉距离，表达式如下：

其中，j＝1,2,…,n，n为标准地震道的长度；

根据X、Y内不同数据之间点与点的欧拉距离，构成n×n的动态时间弯曲距离矩阵d：

在矩阵d中，相邻矩阵元素的集合构成一条弯曲路径，如公式(7)所示：

W＝{w₁,…,w_k,…,w_K} (7)

式中，n≤K≤2n-1；w_k＝(d_ij)_k，表示W中第k个元素x_i与y_j两点的距离；

从求解得到的弯曲路径中选取最短弯曲路径W_best，表达式如公式(8)所示：

为求解公式(8)，由递归思想，构造累积代价矩阵γ，累积距离γ(i,j)为当前格点距离d(i,j)与能够到达该点的最小的邻近元素的累积距离之和：即：

其中，i＝1,2,…,n，j＝1,2,…,n,n为标准地震道与滑动时窗内的波形数据的长度；

因此，X和Y的动态时间规整距离为：DTW(X,Y)＝γ(n,n) (10)

经过上述计算，获得标准地震道X与滑动时窗内波形数据Y的最佳匹配关系，计算出最短弯曲路径W_best，将最短弯曲路径值赋值给滑动时窗中心位置；

原始地震剖面中所有滑动时窗计算完成后，得到最短路径距离矩阵P，公式如下：

其中，最短路径距离矩阵P为M×N的矩阵；

W_pq表示最短路径距离，W_pq值越大表示相似性越差，W_pq值越小表示相似性越高；

其中，p＝1,2,…,M，M为采样点数；

对最短路径距离矩阵P进行如下处理：

式中，W_max表示矩阵P中的最大值，W_min表示矩阵P中的最小值；

Q_pq值越大表示相似度越高；经过公式(12)处理后获得相似性地震属性剖面矩阵Q：

III.根据所要追踪解释的目标层位位置，确定地震剖面层位追踪的纵向深度范围；

设定目标层位的顶部界面深度为h₁，底界面深度为h₂，则定义如下的追踪范围：

h+h₁＜σ＜h+h₂ (14)

其中，σ为层位追踪搜寻的纵向深度范围，h为调节深度；

IV.对将步骤II得到的相似性地震属性剖面矩阵作为输入，以步骤III获得的追踪范围作为约束条件，利用Viterbi算法进行约束优化处理，对目标层位进行识别与追踪；

具体过程如下：

IV.1.地震反射特征最大相似度做向前积分，并记录累加路径，公式如下：

其中，ω(p,q)为步骤II获得的相似性地震属性剖面矩阵Q；

表示以(ik,q-1)为坐标，ik为取值范围，取q-1之前累积的结果中最大值

该值与下一个相似性地震属性剖面矩阵ω(p,q)求和即得ω′(p,q)；

L(p,q)＝ik|_q-1表示求和累积的过程路径；

s为Viterbi算法追踪的时窗大小；

IV.2.搜寻地震反射特征累加路径最大相似度值，向后递归搜寻最短路径，追踪地震反射层；

当公式(15)的累积达到q＝N-1时，搜寻到最大值ω′_max(p,N)，并记下所对应的路径L(p,N-1)＝ik|_N-2为x_ipick|_N-1＝p；

其中，L(p,N-1)表示公式(15)的累积达到q＝N-1，搜寻到最大值ω′_max(p,N)时所对应的路径位置，x_ipick|_N-1表示拾取路径的纵坐标；

根据向后递减跟踪原理搜寻最短路径并计算所需要的值v(q′)，得如式：

其中，q′＝N-1,N-2,…,2,1；x_ipick|_q′＝L(x_ipick|_q′+1,q′)表示所选择累积积分最大相似度值所对应路径中的每个点坐标，即为层位追踪的路径；v(q′)为最短路径中每个点的坐标对应的相似性属性值，ω(x_ipick|_q′)表示最短路径对应的相似性属性值。

Images