CN106650766A - 基于固有特征分析的三维体波形分类方法 - Google Patents

Abstract

本发明公开了一种基于固有特征分析的三维体波形分类方法，包括以下步骤：S1、沿解释层位截取数据体，获得原始数据；S2、对原始数据采用结构导向滤波方法进行去噪处理；S3、对去噪处理后的数据进行IFA特征降维处理；S4、对降维后的数据进行SOM聚类；S5、利用SOM聚类数据生成波形分类相图。本发明通过对地震数据进行IFA降维处理，降维后得到的特征是每道波形的固有特征，能在降低特征之间的冗余度、减小数据量的同时，保持数据点之间的原始近邻关系，并且挖掘出隐藏在高维空间中的低维流形结构，保留最具有区分性的特征，能更有效地区分出不同的波形，提高了波形的区分度，有效地提高了波形分类结果的准确率和分辨率。

Description

基于固有特征分析的三维体波形分类方法

技术领域

本发明属于三维体波形分类技术领域，特别涉及一种基于固有特征分析的三维体波形分类方法。

背景技术

三维体波形分类技术是以地震解释层位为基准，沿给定的上下时窗从三维地震数据中截取数据体，对该数据体对应的波形进行预处理及特征选择来建立样本集，并对该样本集进行训练以建立合适有效的分类器，将目的层段的地震道波形进行分类，达到划分地震相的目的。由于相同的波形对应的是相同的地震相，不同的波形对应的是不同的地震相，所以对波形进行分类的过程也是进行地震相划分的过程。划分的结果是属于同一个簇中的波形相似度较高，而不同簇中的波形相似度较小，不同的簇对应不同的波形类别，最终达到区分不同波形的目的。根据划分结果可以识别地震相单元，分析和解释得到的地震相图，识别区域沉积相和沉积体系，并进行储层的预测。波形分类技术广泛应用于各个领域，对科技的发展和人们的日常生活都产生了很多积极的影响。

发明内容

本发明的目的在于克服现有技术的不足，提供一种通过对地震数据进行IFA降维处理，能在降低特征之间的冗余度、减小数据量的同时，保留最具有区分性的特征，能更有效地区分出不同的波形，有效地提高波形分类结果的准确率的基于固有特征分析的三维体波形分类方法。

本发明的目的是通过以下技术方案来实现的：基于固有特征分析的三维体波形分类方法，包括以下步骤：

S1、沿解释层位截取数据体，获得原始数据；

S2、对原始数据采用结构导向滤波方法进行去噪处理；

S3、对去噪处理后的数据进行IFA特征降维处理；

S4、对降维后的数据进行SOM聚类；

S5、利用SOM聚类数据生成波形分类相图。

进一步地，所述步骤S2包括以下子步骤：

S21、采用了梯度结构张量来对原始三维体数据进行特征提取，定义梯度结构张量为：

其中，u(n)为原始三维体数据，x、y、z分别表示三维体数据三个维度上的坐标；

S22、根据区域结构性信息设计扩散矩阵：

v₁、v₂、v₃分别表示当下迭代时的梯度结构张量的特征值对应的特征向量；

S23、设立连续性因子：

ε的值域为[0,1]，在平滑区域接近1，在图像边缘区域接近0；S⁰代表初始的(也就是进行迭代滤波前的)梯度结构张量矩阵；S_ρ代表当下迭代时的梯度结构张量；Tr(·)表求括号内的矩阵的迹，即求主对角元素的和；

S24、根据扩散方程进行迭代滤波，得到滤波后的三维体数据u(n+1)：

G_σ为一个方差为σ的高斯滤波函数；表示梯度；Δt表示迭代步长；表示求解的散度。

进一步地，所述步骤S3具体实现方法为：将步骤S2滤波的到的三维体数据u(n+1)作为IFA特征降维的待分类样本的原始高维特征集其中，i表示样本序号，n表示样本总数，d表示样本的原始特征维度；设表示降维之后的特征集，其中，d'表示降维之后的特征维度，d＜＜d'；则IFA的问题表述为：求解一组非线性函数将X映射到低维度空间X'；

具体包括以下子步骤：

S31、构建一个k最近邻图，k最近邻图是一个无向有权图，并计算k最近邻图的拉普拉斯矩阵L；

S32、求解下式中的广义特征值问题，得到最小的d'个特征值对应的特征向量集

(I+εLK)α＝λKα；

其中，I∈R^n×n，是一个单位矩阵；K是定义在原始高维特征集X上的格拉姆矩阵，其元素K_i,j＝K(x_i,x_j)；

S33、将求得的特征向量集代入下式，求解非线性映射函数：

v表示出训练样本外的新样本；

S34、通过步骤S33求得的非线性映射函数将原始高维数据映射到低维空间R^d'中。

进一步地，所述步骤S4实现方法为：实现SOM聚类的目标函数如下所示：

其中，m是任何比1大的实数，x_i是第d维数据，U_ij是x_i在簇j中的度，C_j是簇的d维中心；

簇的中心点的计算公式如下：

将以上过程结合DBI指数，作为聚类的评估指标，对数据进行更精确的预测和制图，当DBI值最小时，预示着分类结果最好；

DBI的计算公式如下：

s_i代表第i个簇内各样本的均值；s_j代表第j个簇内各样本的均值；d_ij代表第i个簇和第j个簇中心的距离；表示求解当i≠j时，的最大值；N是簇的个数。

本发明的有益效果是：

1、由于采集地震资料的范围很广，数据量巨大，地震道的特征较多，直接使用这些数据进行波形分类，计算复杂度和时间复杂度都很高，通过IFA算法可以降低特征集维度，提高算法的适用性和计算效率；

2、地震资料采集过程中，环境复杂，存在诸多不确定因素，地震波形中会掺杂一定的噪声，在进行波形分类时，通过对地震数据进行IFA降维处理，降维后得到的特征是每道波形的固有特征，能在降低特征之间的冗余度、减小数据量的同时，保持数据点之间的原始近邻关系，并且挖掘出隐藏在高维空间中的低维流形结构，保留最具有区分性的特征，能更有效地区分出不同的波形，提高了波形的区分度，有效地提高了波形分类结果的准确率和分辨率。

附图说明

图1为本发明的三维体波形分类方法流程图；

图2为一维神经网络结构示意图；

图3为二维神经网络结构示意图；

图4为三维神经网络结构示意图；

图5为SOM训练流程图。

具体实施方式

下面结合附图进一步说明本发明的技术方案。

如图1所示，本发明的一种基于固有特征分析的三维体波形分类方法，包括以下步骤：

S1、沿解释层位截取数据体，获得原始数据；

S2、对原始数据采用结构导向滤波方法进行去噪处理；包括以下子步骤：

S21、采用了梯度结构张量来对原始三维体数据进行特征提取，定义梯度结构张量为：

其中，u(n)为原始三维体数据，x、y、z分别表示三维体数据三个维度上的坐标；

S22、根据区域结构性信息设计扩散矩阵：

v₁、v₂、v₃分别表示当下迭代时的梯度结构张量的特征值对应的特征向量；

S23、设立连续性因子：

ε的值域为[0,1]，在平滑区域接近1，在图像边缘区域接近0；S⁰代表初始的(也就是进行迭代滤波前的)梯度结构张量矩阵；S_ρ代表当下迭代时的梯度结构张量；Tr(·)表求括号内的矩阵的迹，即求主对角元素的和；

S24、根据扩散方程进行迭代滤波，得到滤波后的三维体数据u(n+1)：

G_σ为一个方差为σ的高斯滤波函数；表示梯度；Δt表示迭代步长；表示求解的散度。

结构导向滤波算法通过构建扩散张量矩阵来提取三维图像中的局部区域结构信息，在计算该矩阵的过程中，能够找出相应的地震信号最大相干性的走向，从而确定后续的扩散滤波过程。同时，为了保护图像的边缘结构信息，该算法利用结构张量矩阵来确定图像的边缘，以在适当的位置停止扩散过程，从而更好地保留有用的信息。

S3、对去噪处理后的数据进行IFA特征降维处理；IFA属于流形学习的范畴，是一种数据驱动的无监督学习方法，属于信号处理技术。IFA利用图的拉普拉斯谱聚类机制，去近似一个在固有坐标系中的非线性投影图。在该固有坐标系统中，对于一个d维流形，可以提取出其在R^d中的基本嵌入，并能将隐藏在原始高维数据中的低维流形结构进行明确的参数化表示，恢复数据本身潜在的位形空间，同时测地距离(流形上的最短距离的长度)可以用简单的欧几里得距离来代替。IFA可以保持数据中的非线性信息，从原始的高维数据中提取出隐藏的低维流形结构，该结构可以从本质上反映样本之间的内在关系，表征原始数据的固有特征，并且特征维度较低。进行IFA特征降维处理具体实现方法为：将步骤S2滤波的到的三维体数据u(n+1)作为IFA特征降维的待分类样本的原始高维特征集其中，i表示样本序号，n表示样本总数，d表示样本的原始特征维度；设表示降维之后的特征集，其中，d'表示降维之后的特征维度，d＜＜d'；则IFA的问题表述为：训练一组非线性函数将X映射到低维度空间X'；并且对于新的待降维样本，可以利用该非线性函数直接求得其低维度表示。

具体包括以下子步骤：

S31、构建一个k最近邻图，k最近邻图是一个无向有权图，并计算k最近邻图的拉普拉斯矩阵L；

S32、求解下式中的广义特征值问题，得到最小的d'个特征值对应的特征向量集

(I+εLK)α＝λKα；

其中，I∈R^n×n，是一个单位矩阵；K是定义在原始高维特征集X上的格拉姆矩阵，其元素K_i,j＝K(x_i,x_j)；

S33、将求得的特征向量集代入下式，求解非线性映射函数：

v表示出训练样本外的新样本；

S34、通过步骤S33求得的非线性映射函数将原始高维数据映射到低维空间R^d'中。

在IFA的四个步骤中，前三步可以利用待分类的部分数据集X在线下计算，得到非线性基函数中的特征向量集对于流形上的其他样本数据，可以直接利用这些求得的结果进行后续映射基函数以及映射结果的计算。在这一点上，该方法明显优于拉普拉斯特征映射以及其它流形学习方法。这样不仅仅减少了训练集数量和训练过程中的计算量，同时，对于训练集外的新样本，不需要再次重新进行整个流程的计算，体现了很好的泛化能力，对于海量数据的处理具有很好的适用性。

S4、对降维后的数据进行SOM聚类；SOM聚类的目的在于把相似的数据聚成一类，并与其他簇分隔开来。

一般而言，一个最优的聚类算法的目标是最小化每一类中各自的元素的总距离，同时，最大化不同簇之间的距离。自组织映射神经网络由输入层和输出层构成，是一个两层前馈神经网络。该算法把高维输入数据之间的统计关系简化成单纯的输出层中输出神经节点之间的集合关系。根据排列方式的不同，输出层神经节点可以分为多种结构，不管是何种结构，在输出网络上都保持了高维输入数据的相似性度量关系和拓扑关系。当输出层分别为一维、二维、三维时，SOM神经网络的基本结构分别如图2、图3、图4所示。

对于输入节点即地震属性样本，相当于n维空间中的一个点，SOM旨在重新组织这些高维数据点，将这些数据描画到一个几何结构上，该几何结构称为SOM，本发明中该几何结构即为二维平面输出层。每一个SOM单元被定义为一个输出神经元，其与周围相邻的点连接，通常形成一个六边形或者矩形结构图。

假设输出层平面中有P个元素，即有P个输出神经元mi,i＝1,2,…,P，每个输出神经元可以表示为一个一维向量：m_i＝[m_i1,m_i2,...,m_in]，其中n是训练样本的维度，即输入的地震属性的数量。在进行了SOM训练之后，输出神经元可以用来代表输入的地震属性数据集。

输出层网络中输出神经元的个数决定了算法的有效性和泛化能力。在训练过程中，SOM形成一个弹性网络，并不断进行调节，使其与输入的地震属性数据相适应。在输入空间中相近的数据，在输出的图中依旧是相近的。SOM可以认为是一个将n维输入空间中的点集映射到一个二维网格中，并且能保持原始的拓扑结构。又由于地震数据表达了地质情况的变化，因此SOM保持了地下结构的拓扑关系。

通常，SOM输出神经元是迭代更新的。如图5所示，SOM训练过程详述如下：

1)初始化输出神经元：可以用随机的方式来进行初始化，也可以利用输入数据最大的两个特征值对应的特征向量来进行映射初始化。

2)随机地从输入数据集中选择一个输入向量x，计算x和输出神经元之间的距离；将输出神经元中与输入向量x距离最小的向量，称为获胜神经元，设为m_b，

获胜神经元的计算公式如下：

3)更新获胜神经元和与其邻近的点，第i个神经元的更新规则如下：

m_i(t+1)＝m_i(t)+λ(t)h_bi(t)[x-m_i(t)]

其中，t代表迭代次数，λ(t)是学习率，h_bi(t)是获胜神经元的领域大小。h_bi(t)变量的值随着学习过程中迭代次数的增加而减小，计算公式如下：

其中，r_b和r_i是输出神经元b和i在SOM网格中的位置，δ²(t)定义领域的宽度。

4)重复步骤2)、3)，直到达到迭代上限或者获胜神经元的领域半径达到了预定义的最小值。由此，便得到了输出层的网络结构。

本发明对降维后的数据进行SOM聚类的具体实现方法为：实现SOM聚类的目标函数如下所示：

其中，m是任何比1大的实数，x_i是第d维数据，U_ij是x_i在簇j中的度，C_j是簇的d维中心；

簇的中心点的计算公式如下：

将以上过程结合DBI指数，作为聚类的评估指标，对数据进行更精确的预测和制图，当DBI值最小时，预示着分类结果最好；

DBI的计算公式如下：

s_i代表第i个簇内各样本的均值；s_j代表第j个簇内各样本的均值；d_ij代表第i个簇和第j个簇中心的距离；表示求解当i≠j时，的最大值；N是簇的个数。

S5、利用SOM聚类数据生成波形分类相图。

本领域的普通技术人员将会意识到，这里所述的实施例是为了帮助读者理解本发明的原理，应被理解为本发明的保护范围并不局限于这样的特别陈述和实施例。本领域的普通技术人员可以根据本发明公开的这些技术启示做出各种不脱离本发明实质的其它各种具体变形和组合，这些变形和组合仍然在本发明的保护范围内。

Claims (4)

1.基于固有特征分析的三维体波形分类方法，其特征在于，包括以下步骤：

S1、沿解释层位截取数据体，获得原始数据；

S2、对原始数据采用结构导向滤波方法进行去噪处理；

S3、对去噪处理后的数据进行IFA特征降维处理；

S4、对降维后的数据进行SOM聚类；

S5、利用SOM聚类数据生成波形分类相图。

2.根据权利要求1所述的基于固有特征分析的三维体波形分类方法，其特征在于，所述步骤S2实现方法为：

S21、采用了梯度结构张量来对原始三维体数据进行特征提取，定义梯度结构张量为：

$S=\▿u\·\▿u^T=\left[\begin{array}{ccc}\frac{\∂u}{\∂x}\·\frac{\∂u}{\∂x}& \frac{\∂u}{\∂x}\·\frac{\∂u}{\∂y}& \frac{\∂u}{\∂x}\·\frac{\∂u}{\∂z}\\ \frac{\∂u}{\∂y}\·\frac{\∂u}{\∂x}& \frac{\∂u}{\∂y}\·\frac{\∂u}{\∂y}& \frac{\∂u}{\∂y}\·\frac{\∂u}{\∂z}\\ \frac{\∂u}{\∂z}\·\frac{\∂u}{\∂x}& \frac{\∂u}{\∂z}\·\frac{\∂u}{\∂y}& \frac{\∂u}{\∂z}\·\frac{\∂u}{\∂z}\end{array}\right]$

其中，u(n)为原始三维体数据，x、y、z分别表示三维体数据三个维度上的坐标；

S22、根据区域结构性信息设计扩散矩阵：

$D=\[\begin{array}{ccc}{v}_1& v_2& v_3\end{array}\]*\left[\begin{array}{ccc}0& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right]*\left[\begin{array}{c}{v}_1^T\\ v_2^T\\ v_3^T\end{array}\right]$

v₁、v₂、v₃分别表示当下迭代时的梯度结构张量的特征值对应的特征向量；

S23、设立连续性因子：

$\mathrm{\ϵ}=\frac{Tr(S^0\·S_{\mathrm{\ρ}})}{Tr\left(S^0\right)\·Tr\left(S_{\mathrm{\ρ}}\right)}$

S⁰代表初始的梯度结构张量矩阵；S_ρ代表当下迭代时的梯度结构张量；Tr(·)表求括号内的矩阵的迹，即求主对角元素的和；

S24、根据扩散方程进行迭代滤波，得到滤波后的三维体数据u(n+1)：

$u(n+1)=u\left(n\right)+\mathrm{\Δ}t\·div\{\mathrm{\ϵ}\·D\[G_{\mathrm{\σ}}*(\▿u\·\▿u^T)\]\▿u\}$

G_σ为一个方差为σ的高斯滤波函数；表示梯度；Δt表示迭代步长；

表示求解的散度。

3.根据权利要求2所述的基于固有特征分析的三维体波形分类方法，其特征在于，所述步骤S3具体实现方法为：将步骤S2滤波的到的三维体数据u(n+1)作为IFA特征降维的待分类样本的原始高维特征集其中，i表示样本序号，n表示样本总数，d表示样本的原始特征维度；设表示降维之后的特征集，其中，d'表示降维之后的特征维度，d＜＜d'；则IFA的问题表述为：求解一组非线性函数将X映射到低维度空间X'；

具体包括以下子步骤：

S31、构建一个k最近邻图，k最近邻图是一个无向有权图，并计算k最近邻图的拉普拉斯矩阵L；

S32、求解下式中的广义特征值问题，得到最小的d'个特征值对应的特征向量集(I+εLK)α＝λKα；

其中，I∈R^n×n，是一个单位矩阵；K是定义在原始高维特征集X上的格拉姆矩阵，其元素K_i,j＝K(x_i,x_j)；

S33、将求得的特征向量集代入下式，求解非线性映射函数：

$f_j^{*}={\mathrm{\Σ}}_{i=1}^n{\mathrm{\α}}_i^{j}K(x_i,v);$

v表示出训练样本外的新样本；

S34、通过步骤S33求得的非线性映射函数将原始高维数据映射到低维空间R^d'中。

4.根据权利要求3所述的基于固有特征分析的三维体波形分类方法，其特征在于，所述步骤S4实现方法为：实现SOM聚类的目标函数如下所示：

$J_m={\mathrm{\Σ}}_{i=1}^N{\mathrm{\Σ}}_{j=1}^c{u}_{ij}^m\left|\right|x_i-C_j|{|}^2$

其中，m是任何比1大的实数，x_i是第d维数据，U_ij是x_i在簇j中的度，C_j是簇的d维中心；

簇的中心点的计算公式如下：

$C_j=\frac{{\mathrm{\Σ}}_{i=1}^N{u}_{ij}^m{X}_i}{{\mathrm{\Σ}}_{i=1}^N{u}_{ij}^m}$

将以上过程结合DBI指数，作为聚类的评估指标，对数据进行更精确的预测和制图，当DBI值最小时，预示着分类结果最好；

DBI的计算公式如下：

$DBI=\frac{1}{N}{\Σ}_{i=1}^N{\max }_{i\≠j}\left\{\frac{s_i+s_j}{d_{ij}}\right\}$

s_i代表第i个簇内各样本的均值；s_j代表第j个簇内各样本的均值；d_ij代表第i个簇和第j个簇中心的距离；表示求解当i≠j时，的最大值；N是簇的个数。