CN110471104A - 基于智能特征学习的叠后地震反射模式识别方法 - Google Patents

Abstract

本发明提出了一种基于智能特征学习的叠后地震反射模式识别方法，属于地震反射模式识别技术领域。本发明通过同步压缩小波变换对叠后地震信号进行时频域变换，得到高精度的频谱图，构建叠后地震数据矩阵，并对其进行非负矩阵分解从而得到叠后地震信号的特征，对得到的叠后特征进行分类并生成叠后地震相图。本发明对于叠后地震信号有非常好的特征表征能力，能从大量无标签数据中学习数据分布特性，从而达到了较高的叠后地震相反射模式识别能力。

Description

基于智能特征学习的叠后地震反射模式识别方法

技术领域

本发明属于地震反射模式识别技术领域，特别涉及一种基于智能特征学习的叠后地震反射模式识别方法。

背景技术

随着地震地层学在勘探领域的普及，利用地震相进行地质研究的做法受到了广泛的应用。地震相反映了同一区域地震剖面上沉积相表现的总和。地震相划分是依据剖面上地震反射模式识别所确定的地震相类型开展的，它是进行沉积相研究的一种强有力的方法。基于地震信号开展地震反射模式识别，是划分地震相的一种重要方法。如今业界采集地震信号的技术水平不断提升，地震数据中包含的地震信息也愈加丰富，所以通过人工方式对地震反射模式进行识别来描述地震相，具有很大的主观性和不确定性。为了实现自动定量分析地震相的目的，需要依靠地震数据处理技术、信息处理技术和计算机技术，并通过数学表示方法对地震反射模式加以定量表征和精准识别。

准确的特征是完整表征地震反射模式的基础，通过智能提取地震反射信号的特征来增强地震反射模式识别能力的方法，能够进一步提升地震相划分的精度。其结果可以展示地下地质结构的分布，为探井方位的判定提供参考依据，降低钻井失败的可能性，节约勘探成本。

地震反射模式识别的流程，首先要对地震信号进行预处理，在层位上选择合适的时窗大小，获得目标层位段的地震信号。然后结合地震信号特征提取方法和机器学习分类算法，对目标层段地震信号反射特征进行分类，并通过类标签的方式区分对应的地震相，进而就可以对各种目标地质结构的分布情况进行研究。在地震反射模式识别的流程中，提取特征和特征分类是最为关键的两个步骤，其中对地震反射模式进行识别的前提是可以提取到完整表征目标信号的特征。

由于地震波形具有直观的解释意义，因此在地震相划分过程中通常会采用波形分类技术，使其成为当前地震反射模式识别领域内中最关键的方法。最简洁时域波形特征就是直接的波形时间序列，直接输入波形时间序列作为特征，然后使用竞争神经网络进行波形模式识别。进一步地，利用邻域内的多道波形时间序列作为输入进行分类，通过在时域中利用多道特性，可以减小噪声对波形特征提取过程的影响。提取波形的峰值、过零点、峰值时间等统计信息作为时域波形，相比直接时间序列更精确，也更稳健。时域波形特征缺点是稳定差，不确定性强，尤其是对层位解释误差敏感。在功率谱中提取了功率谱总能量，以达到10％，20％，30％，...功率谱总能量所需时间作为频域波形特征。国内研究院认为地震波和语音信号在介质中传播具有相同的物理本质，所以提出了用语音识别中应用非常成功的Mel系数来表示地震波形特征。无论时域还是频域特征对波形的特征表征都是粗糙的、非稳健的，所以现在波形分类方法的研究重点是在时频域波形特征提取上。利用人工免疫神经网络优越的降维特性，对小波变换特征进行降维去除冗余特征，再利用神经网络进行聚类，取得了十分稳健效果。基于小波变换的时频域波形特征并非最好的，利用经验模态分解方法来提取波形的时频特征，其比小波分解模型能够更多地保留波形特征，所以具有更高的波形分类分辨率。地震波形和语音信号物理本质是类似的，不同的是采用了音乐领域里面的特征更精细地表征地震波形，具体思路是首先采用时频变换(比如小波变化，S变换等)获得地震波形的谱图，然后从谱图中提取出地震波形音乐属性。在国内，在地震反射波形分类研究中提出加入分数阶导数，构造出由子波分数阶导数组成的波形集合，其能够精确匹配地震信号中每一道地震波形，从而对地震反射模式进行识别。

发明内容

本发明的目的在于实现一种利用矩阵分解的方法对叠后地震信号进行特征学习，从而提升叠后地震反射模式的识别效果。本发明提出了一种基于智能特征学习的叠后地震反射模式识别方法，引入同步压缩小波变换对叠后地震信号进行时频变换，在时频域中构造数据矩阵并对其进行非负分解，从而将计算得到的稀疏表达矩阵作为学习到的特征，然后对学习到的特征进行分类，并生成最终的叠后地震相图。

一种基于智能特征学习的叠后地震反射模式识别方法，包括以下步骤：

S1、获取叠后地震信号，采用同步压缩小波变换(Synchrosqueezed WaveletTransform,SST)对所述叠后地震信号进行时频变换；

S2、在时频域中构建数据矩阵；

S3、对所述数据矩阵进行非负分解，得到叠后地震信号特征。

进一步地，所述步骤S1包括：

S11、获取叠后地震信号，对所述叠后地震信号中的一道数据s(t)进行连续小波变换(Continuous Wavelet Transform,CWT)，结果为

其中，ψ^*表示母小波ψ(t)复数共轭，b表示母小波的时间平移因子，a表示母小波的尺度因子；W_s(a,b)表示小波系数，将所述小波系数映射到二维平面得到时间-尺度聚合的图像；

S12、根据时间域和频率域的能量相等，重写所述小波系数为

其中，表示s(t)的傅里叶变换，ξ表示s(t)的角频率；

所述母小波的傅里叶变换的能量集中在ψ＝ω₀处，s(t)的连续小波变换W_s(a,b)在时间-尺度平面内聚集在上尺度的直线，其中，ω₀表示小波中心频率，ω表示信号中心频率；在任意点(a,b)处通过计算小波变换的偏导得到其瞬时频率ω_a(a,b)，对于满足W_s(a,b)≠0的小波系数，得到

S13、计算所述叠后地震信号的SST量值；

通过ω_a(a,b)公式将信号参数从时间-尺度平面映射到时间-频率平面上，把每个点(b,a)转换为(b,ω_a(a,b))，完成同步压缩，得到基于连续小波变换的时频分布；通过尺度离散化计算W_s(a,b)，尺度间隔为Δa_k＝a_k-1-a_k；当信号从时间-尺度平面映射到时间-频率平面的过程中，SST的量值T_s(ω,b)以ω_l为中心的范围决定

其中，a_k表示第k个离散尺度，满足区间|ω_s(a_k,b)-ω_l|≤Δω/2，ω_l表示第l个离散角频率。

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

将每一道时频图像S按照时间轴的方向分成M个长度为Q＝T/M的切片，每个切片互相不重叠，其中，T表示时间轴方向长度；S_m表示第m个谱图切片，信号经过同步压缩小波变换得到时频图像S可以表示为一组连续的时频切片S＝[S₀,...S_M-1]；分别对每个切片进行时间平均，每个时频图像转变为M个向量的组合；所述叠后地震信号有L个道即L个时频图像，第l个时频图像由一组向量表示，其中表示通过对时频切片求时间平均而获得的大小为P向量；对所有L个时频图像提取L组向量D^(l)，并逐列堆叠构建数据矩阵其中，D＝[D⁽¹⁾,...D^(L)]，N＝ML。

进一步地，所述步骤S3包括：

基于非负矩阵分解(Nonnegative Matrix Factor,NMF)的叠后地震信号特征学习；

计算所述数据矩阵的模并进行归一化处理，令其中，表示进行归一化处理后的数据矩阵，P表示每一个样本的维度，N表示样本数目；

对X进行NMF运算，获得基矩阵U和系数矩阵V，其中，V中的列向量为特征向量，即用于表征叠后地震反射模式的叠后地震信号特征；

NMF的损失函数为

对所述损失函数求解获得最优解，目标函数为

其中，Y₁和Y₂为拉格朗日乘子，tr(·)表示矩阵的迹，上标^T表示矩阵的转置；针对所述目标函数，对U和V求偏导数，为

根据所述目标函数的KKT条件，有(Y₁ ^TU)_ik＝0，其中i,j,k表示矩阵元素下标，得到U和V的迭代更新公式，对U和V进行迭代更新；所述迭代更新公式为

进一步地，其特征在于，所述方法还包括：

S4、对所述叠后地震信号特征进行分类，并生成叠后地震相图。

本发明的有益效果：本发明提供了一种基于智能特征学习的叠后地震反射模式识别方法，本发明通过同步压缩小波变换对叠后地震信号进行时频域变换，得到高精度的频谱图，构建叠后地震数据矩阵，并对其进行非负矩阵分解从而得到叠后地震信号的特征，对得到的叠后特征进行分类并生成叠后地震相图。本发明对于叠后地震信号有非常好的特征表征能力，能从大量无标签数据中学习数据分布特性，从而达到了较高的叠后地震相反射模式识别能力。

附图说明

图1为本发明实施例的流程图。

图2为本发明实施例的步骤S2的示意图。

图3为本发明实施例的步骤S3的示意图。

图4为本发明实施例的工区振幅属性图。

图5为本发明方法的结果图。

图6为PCA方法的结果图。

图7为本发明实施例的实际工区振幅属性图。

图8为本发明方法的另一结果图。

图9为PCA+SOM方法的结果图。

具体实施方式

下面结合附图对本发明的实施例做进一步的说明。

请参与图1，本发明提出了一种基于智能特征学习的叠后地震反射模式识别方法，本发明首先将对叠后地震信号进行时频变换，本方法引入同步压缩小波变换(Synchrosqueezed Wavelet Transform,SST)，SST能够有效地提取出信号的频率分类，对复杂多成分谐波信号的分解效果很好，为后续步骤奠定了良好的基础。地震信号在经过时频变化之后通常会使数据的维度增加，比如一道地震信号会从向量变为矩阵，由于地震信号具有稀疏的特性，所以这些矩阵中有用信息只集中在很小的区域内，大部分矩阵元素为0。本发明提出数据矩阵的构建方法，对数据的时间维度进行适应性的预处理，通过切片和池化操作减少数据的维度，从而降低后续提取特征的计算量。最后引入非负矩阵分解算法对叠后地震数据矩阵进行分解和优化，该算法从处理好的数据矩阵中学习字典，然后使用该字典上的数据的投影作为学习到的叠后地震信号特征。非负矩阵分解是一种非线性的分解方式，能够细致学习叠后地震信号的复杂特征。具体通过以下步骤实现：

S1、获取叠后地震信号，采用同步压缩小波变换对叠后地震信号进行时频变换。

本实施例中，由于真实地下地质构造复杂多变，采集到的叠后地震信号数据往往十分杂乱，数据中所包含的地质特征信息很容易被掩盖或产生畸变。如果直接对采集到数据按照典型的信号处理方法进行特征提取，可能存在特征空间维数过高、冗余信息过以及过度依赖专家知识的问题。地球物理学家经常用时频表示来处理地震信号，通过不同的时频表示，地震数据的频谱中包含的深度信息可以被挖掘出来。传统的时频表示如连续小波变换(Continuous Wavelet Transform,CWT)能够将一个信号分解成时间-频率成分，但无法同时在时间方向和频率方向得到最佳的清晰度。SST是一种基于小波变换的时频表示方法，在小波尺度方向上，SST通过对时间-尺度平面上的能量进行重新分配，并将尺度维度变换为频率维度，最终使其频率曲线更加集中。采用SST处理叠后地震信号包括以下子步骤：

S11、采用连续小波变换处理叠后地震信号

s(t)为叠后地震信号中的一道数据，其CWT为

其中，ψ^*表示母小波ψ(t)复数共轭，b表示母小波的时间平移因子，a表示母小波的尺度因子。小波变换可以理解为信号s(t)和原始母小波通过波形压缩或扩展后在时间轴上平移后的新小波的互相关。W_s(a,b)表示小波系数，将计算得到的小波系数映射到二维平面，即可得到一幅时间-尺度(频率)聚合的图像。

S12、计算叠后地震信号的瞬时频率根据Parseval理论，根据时间域和频率域的能量相等，重写小波系数为

其中，表示s(t)的傅里叶变换，ξ表示s(t)的角频率。经过变换，CWT从公式(1)中的卷积形式变成了公式(2)中的乘积形式。为了便于理解，以信号s(t)＝Ac(oωs)为例进行解释，此时s(t)经过傅里叶变换为：将该式代入公式(2)中有：

母小波的傅里叶变换的能量集中在ψ＝ω₀处，理想情况下s(t)的连续小波变换W_s(a,b)在时间-尺度平面内聚集在上尺度的直线，其中，ω₀表示小波中心频率，ω表示信号中心频率。但通常在实际情况下，CWT的能量是以这条直线为中心向周边范围分散，在频率(尺度)方向出现了一些的不清晰的纹路，使准确观察信号时频表示时的难度增加。如果将时间轴周边不清晰的纹路忽略，就可以在任意点(a,b)处通过计算小波变换的偏导得到其瞬时频率ω_a(a,b)，对于满足W_s(a,b)≠0的小波系数，得到

S13、计算叠后地震信号的SST量值

通过公式(4)将信号参数从时间-尺度平面映射到新的时间-频率平面上，把每个点(b,a)转换为(b,ω_a(a,b))，这个步骤为同步压缩。这样就得到基于连续小波变换的时频分布，通过尺度离散化计算W_s(a,b)，其尺度间隔为Δa_k＝a_k-1-a_k。当信号从时间-尺度平面映射到时间-频率平面的过程中，SST的量值T_s(ω,b)以ω_l为中心的范围[ω_l-Δω/2,ω_l+Δω/2]决定

其中，a_k表示第k个离散尺度，满足区间|ω_s(a_k,b)-ω_l|≤Δω/2，ω_l表示第l个离散角频率。

通过以上公式和约束条件，可以在尺度(频率)方向将叠后地震信号中每一道数据CWT值的模糊区域挤压到十分逼近真实的频率值的范围，这很大程度提升了信号在时频方向上的集中能力，可以清晰、精准地描述地震信号，从而提高叠后地震信号特征提取效率和准确率。

S2、在时频域中构建数据矩阵。

本实施例中，为了最终得到叠后地震信号的特征，在得到每一道信号数据对应的时频图像后，需要将这些图像拼接在一起来进行分解。由于叠后地震信号通常含有上万道数据，直接拼接后得到的矩阵通常维度巨大，分解这种不合理矩阵大小的时频图像计算量巨大，所以对这些时频图像进行降维是十分有必要的。现有的矩阵分解技术一般在其原始公式中不包括任何形式的时间建模，即对数据的时间维度进行适应性的预处理。因此，我们应用两个简单的切片和池化步骤，旨在减少数据的维度，从而降低计算量，同时为特征学习步骤提供合适的时频表示。

如图2所示，首先将每一道时频图像S按照时间轴的方向分成M个长度为Q＝T/M的切片，每个切片互相不重叠，其中，T表示时间轴方向长度。S_m表示第m个谱图切片，信号经过同步压缩小波变换得到时频图像S可以表示为一组连续的时频切片S＝[S₀,...S_M-1]。然后分别对每个切片进行时间平均，这时每个时频图像转变为M个向量的组合。令某叠后地震信号有L个道即L个时频图像，则第l个时频图像由一组向量表示，其中表示通过对时频切片求时间平均而获得的大小为P向量。最后，对所有L个时频图像提取L组向量D^(l)，并逐列堆叠构建数据矩阵其中，D＝[D⁽¹⁾,...D^(L)]，N＝ML。

以这种方式改变时频图像的尺寸有助于表示不同时间段内场景的频率特征。一方面，每个时频切片包含场景中发生的各种事件的时频信息。矩阵分解步骤将在字典中收集最频繁事件的表示，即平均后的时频切片。另一方面，发生频率较低的事件不会对字典的构造产生很大影响。因此，仅对表征场景的最相关的事件进行建模，并且对于足够大小的字典，对这些基础事件的投影将能够识别绝大多数叠后地震反射模式。

S3、对数据矩阵进行非负分解，得到叠后地震信号特征。

本实施例中，叠后地震信号的特征提取是识别叠后地震反射模式的关键，由于地下地质结构的复杂性，地震信号信号通常为非线性非平稳性信号。为了更贴合叠后地震信号的特点，本发明引用非负矩阵分解(NMF)作为矩阵分解的方法来学习叠后信号的特征。

NMF是一种矩阵分解算法，以矩阵当中各元素非负性为约束条件，并进行非线性降维，所以提取到的特征更具实际物理意义。相比于其他传统的降维方法，利用NMF学习叠后地震信号特征主要有三个优点：非线性分解，能够细致学习数据的复杂特征，降维后所抽取的特征更容易解释；收敛速度快，能够快速处理以矩阵形式存储的大规模数据；分解结果具有稀疏性，更贴近地震信号的特性。

NMF原理如图3所示，令原始数据矩阵为X每列均表示数据样本。NMF通过用两个非负数据矩阵U和V的乘积形式来表示原始数据矩阵X。其中，V表示X映射在基空间的组合系数，V中的每一列为X对应的特征向量。U则是基空间，其中每一列都是基向量。一般来说r＜＜min(m,n)，并且X和U、V之间须满足X≈UV。通过U空间，将高维数据X映射到维数较低的V是NMF的主要作用。NMF本质是一种非线性矩阵投影，其中非线性是由于NMF的非负约束使得在分解中只能进行基的非减操作。虽然在计算中NMF看起来像线性操作，但是由于它的维度会在投影过程发生改变，因此NMF是非线性的，这一方面是为了使计算具有一定的稀疏度；另一方面，也可以使分解后的数据更容易解释。很明显，如果没有有效性的限制，这种分解方式能够用任意低维空间中的表示代替高维数据。

本实施例中，采用NMF方法对前两步骤得到的叠后地震信号的数据矩阵提取特征。首先，计算数据矩阵的模并进行归一化处理以满足NMF的初始条件，令新的数据矩阵其中，表示进行归一化处理后的数据矩阵，P表示每一个样本的维度，N表示样本数目。然后对这个数据矩阵进行NMF运算，获得基矩阵U和稀疏矩阵V，其中，V中的列向量为所需的特征向量，由此便通过NMF算法学习到用于表征叠后地震反射模式的叠后地震信号特征。

NMF的损失函数为

对损失函数采用增广拉格朗日乘子法(ALM)求得最优解，目标函数为

其中，Y₁和Y₂为拉格朗日乘子，tr(·)表示矩阵的迹，上标^T表示矩阵的转置。针对(7)中的目标函数，分别对U和V求偏导数，有

根据目标函数(7)的KKT条件，有(Y₁ ^TU)_ik＝0，其中i,j,k表示矩阵元素下标，得到U和V的迭代更新公式(10)和(11)，对U和V进行迭代更新。迭代更新公式为

S4、对叠后地震信号特征进行分类，并生成叠后地震相图。

本实施例中，对得到的特征采用FSOM算法进行分类，并生成最终的叠后地震相图。

本实施例中，以工区数据对本发明方法进行验证。

已知的工区数据振幅属性如图4所示，图5和图6为本发明方法和传统PCA方法对比结果，可以发现，本发明生成的叠后地震相图(图5)更能清晰地识别地下地质结构分布，表明本发明提出的叠后地震反射模式识别效果更好，更说明了本发明提出的智能特征学习方法更能在叠后地震信号中提取揭示其反射模式的特征。

本实施例中，将本发明方法应用于实际工区，图7为实际工区的振幅属性，图8为本发明方法结果，图9为为PCA+SOM方法结果，通过对比可以看出，本发明所提方法对裂隙的走势进行的精确的刻画，并且清晰、平滑地反映了右边圈内结构的变化。这表明本发明提出的方法对叠后地震反射模式的识别效果很好，与真实地质状况比较吻合。

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

Claims (5)

1.一种基于智能特征学习的叠后地震反射模式识别方法，其特征在于，包括以下步骤：

S1、获取叠后地震信号，采用同步压缩小波变换(Synchrosqueezed WaveletTransform,SST)对所述叠后地震信号进行时频变换；

S2、在时频域中构建数据矩阵；

S3、对所述数据矩阵进行非负分解，得到叠后地震信号特征。

2.如权利要求1所述的基于智能特征学习的叠后地震反射模式识别方法，其特征在于，所述步骤S1包括：

S11、获取叠后地震信号，对所述叠后地震信号中的一道数据s(t)进行连续小波变换(Continuous Wavelet Transform,CWT)，结果为

其中，ψ^*表示母小波ψ(t)复数共轭，t表示时间，b表示母小波的时间平移因子，a表示母小波的尺度因子；W_s(a,b)表示小波系数，将所述小波系数映射到二维平面得到时间-尺度聚合的图像；

S12、根据时间域和频率域的能量相等，重写所述小波系数为

其中，表示s(t)的傅里叶变换，ξ表示s(t)的角频率；

所述母小波的傅里叶变换的能量集中在ψ＝ω₀处，s(t)的连续小波变换W_s(a,b)在时间-尺度平面内聚集在上尺度的直线，其中，ω₀表示小波中心频率，ω表示信号中心频率；在任意点(a,b)处通过计算小波变换的偏导得到其瞬时频率ω_a(a,b)，对于满足W_s(a,b)≠0的小波系数，得到

S13、计算所述叠后地震信号的SST量值；

通过ω_a(a,b)公式将信号参数从时间-尺度平面映射到时间-频率平面上，把每个点(b,a)转换为(b,ω_a(a,b))，完成同步压缩，得到基于连续小波变换的时频分布；通过尺度离散化计算W_s(a,b)，尺度间隔为Δa_k＝a_k-1-a_k；当信号从时间-尺度平面映射到时间-频率平面的过程中，SST的量值T_s(ω,b)以ω_l为中心的范围决定

其中，a_k表示第k个离散尺度，满足区间|ω_s(a_k,b)-ω_l|≤Δω/2，ω_l表示第l个离散角频率。

3.如权利要求2所述的基于智能特征学习的叠后地震反射模式识别方法，其特征在于，所述步骤S2包括：

将每一道时频图像S按照时间轴的方向分成M个长度为Q＝T/M的切片，每个切片互相不重叠，其中，T表示时间轴方向长度；S_m表示第m个谱图切片，信号经过同步压缩小波变换得到时频图像S可以表示为一组连续的时频切片S＝[S₀,...S_M-1]；分别对每个切片进行时间平均，每个时频图像转变为M个向量的组合；所述叠后地震信号有L个道即L个时频图像，第l个时频图像由一组向量表示，其中表示通过对时频切片求时间平均而获得的大小为P向量；对所有L个时频图像提取L组向量D^(l)，并逐列堆叠构建数据矩阵其中，D＝[D⁽¹⁾,...D^(L)]，N＝ML。

4.如权利要求3所述的基于智能特征学习的叠后地震反射模式识别方法，其特征在于，所述步骤S3包括：

基于非负矩阵分解(Nonnegative Matrix Factor,NMF)的叠后地震信号特征学习；

计算所述数据矩阵的模并进行归一化处理，令其中，表示进行归一化处理后的数据矩阵，P表示每一个样本的维度，N表示样本数目；

对X进行NMF运算，获得基矩阵U和系数矩阵V，其中，V中的列向量为特征向量，即用于表征叠后地震反射模式的叠后地震信号特征；

NMF的损失函数为

对所述损失函数求解获得最优解，目标函数为

其中，Y₁和Y₂为拉格朗日乘子，tr(·)表示矩阵的迹，上标^T表示矩阵的转置；针对所述目标函数，对U和V求偏导数，为

根据所述目标函数的KKT条件，有(Y₁ ^TU)_ik＝0，其中i,j,k表示矩阵元素下标，得到U和V的迭代更新公式，对U和V进行迭代更新；所述迭代更新公式为

5.如权利要求1-4任一项所述的基于智能特征学习的叠后地震反射模式识别方法，其特征在于，所述方法还包括：

S4、对所述叠后地震信号特征进行分类，并生成叠后地震相图。