本发明涉及地震预测技术领域，特别涉及一种基于hht-mfcc的地震声纹流体预测方法。
背景技术：
希尔伯特-黄变换(hilbert-huangtransform，hht)是n.e.huang于1998年提出的一种全新的时频分析方法，该方法主要包括经验模态分解(empiricalmodedecomposition，emd)和希尔伯特谱分析(hilbertspectralanalysi)两大部分。其中emd为hht的核心，emd将将信号内相互叠加的不同尺度的波动或变化趋势逐一分解开，得到反映信号内部不同震动模式的一组固有本征模态函数(itrinsicmodefunction,imf)。在imf的基础之上利用希尔伯特谱分析方法得到各imf的希尔伯特边际谱，希尔伯特边际谱是具有统计意义的频谱，它体现了每个频率成分对于信号整体总能量的贡献情况，可在一定程度上较为精细地体现出高频、中频、低频子带上能量的分布和变换规律，是分析不同频率成分的能量在频域上分布特征的有效手段。
基于听觉特性的mel频率倒谱系数mfcc在语音识别领域被广泛使用，并且在识别率上体现出了很好的优势，mfcc特征能够很好的表征频谱信息的特点，以一种小数据的方式尽可能的保留了最多的原始频谱信息，并且对于噪声具有较高的鲁棒性，mfcc特征是一种非常有效的表达频谱信息的特征参量。mel频率mel(f)与真实频率f的关系如下：
mel(f)＝1125×ln(1+f/700)；
mel-1(f)＝700×(exp(f/1125)-1)；
传统的倒谱分析算法都采用傅里叶变换将信号从时间域转换到频率域，傅里叶变换是基于平稳线性信号发展了来，实质上是一种全局性的变换，语音信号以及地震信号通常都是复杂的非线性非平稳信号，傅里叶变换的短板日益显现。近年来迅猛发展地希尔伯特黄变换作为一种全新的时频分析方法，可以完全自适应的分解时变的非平稳信号，并且可以同时具备较高的时间分辨率和频率分辨率。通过希尔伯特变换得到信号的希尔伯特边际谱相较于傅里叶谱具有更明确的物理意义，更高的时频分辨率。本发明用具有更高时频分辨率的希尔伯特变换替代mel频率倒谱计算中的傅里叶变换，提取出一种长时特征。
目前储层流体预测已经发展了很多的新技术方法，如亮点、平点、暗点技术，avo、ava技术，弹性阻抗反演等，也取得了一些实质性的效果。另外，近年来低频阴影现象作为储层含流体特征也越来越被人们所重视，地震信号的时频谱分析是进行储层含流体的低频阴影现象识别的基础，常用的时频谱分析技术包括短时傅里叶变换、wigner-ville分布、小波变换、s变换、广义s变换等等。
例如亮点、平点、暗点技术，avo、ava技术，弹性阻抗反演等方法主要对成层性好、均质性强的碎屑岩储层流体有较好地指示。并且受一些严苛条件的限制，要求高保真的地震数据以及储层具有一定的规模且储层界面相对平缓，然而随着勘探目标的日益复杂，这些条件不可能一一满足，因此限制了这些方法在实际工作中的应用。针对低频阴影技术，短时傅里叶变换、小波变换、s变换以及广义s变换均受海森堡测不准原理约束，即信号的持续时间和频带宽度不能同时任意窄，不能同时具备高的时间分辨率和频率分辨率。wigner-ville分布的受海森堡测不准现象相对较小，但其引入的交叉项会影响时频局部精度。
技术实现要素：
本发明针对现有技术的缺陷，提供了一种基于hht-mfcc的地震声纹流体预测方法，能有效的解决上述现有技术存在的问题。
为了实现以上发明目的，本发明采取的技术方案如下：
一种基于hht-mfcc的地震声纹流体预测方法，包括以下步骤：
步骤一：输入叠后地震单道cdp数据，对输入数据进行emd经验模态分解，得到不同频带范围的imf分量:
上式中s(t)为输入的叠后地震单道cdp数据，imfi(t)为经验模态分解得到的各imf分量，rn(t)为分解最后的残余趋势项；
步骤二：对各imf分量分别进行分帧加窗处理；采用有限长度的时窗来截取各imf分量，将其划分成一个一个的短时段信号即一帧。为了保证各帧之间的连续性，相邻帧之间保留一段重叠区域，重叠区域长度取窗函数的1/3-1/2窗长；采用汉明窗，其定义如下：
其中w(n)为窗函数，n为窗长；
步骤三：对分帧加窗后的各imf分量逐帧进行希尔伯特变换得到各帧的希尔伯特边际谱，将各imf分量从时域转换到频域；
步骤四：将步骤三得到的频率域信号通过m个mel三角带通滤波器组获得mel频谱，再将得到的mel频谱取对数能量，得到相应频带的对数能量e(m)，mel三角带通滤波器hm(k)的中心频率f(m)由下式求取：
其中，fs为地震信号的采样频率，fh为带通滤波器的最高频率，一般取fs的一半。fl为带通滤波器的最低频率，一般取值趋于零值。m表示mel三角带通滤波器组中的第m个带通滤波器，n为步骤二中的窗长。
每个滤波器的输出的对数能量e(m)为：
e(m)＝log10(h2(f)hm(f))；
其中h(f)为信号希尔伯特边际谱，hm(f)为三角带通滤波器的频率响应；
步骤五：对步骤四得到的对数能量e(m)进行离散余弦变换dct，去除e(m)间的相关性，得到m阶mel频率倒谱系数mfcc:
步骤六：对输入地震的每个cdp信号循环步骤一到步骤五，即可得到原始二维地震信号的m阶地震声纹参数；
步骤七：对步骤六的结果进行k均值聚类分析：首先从包含n个元素的对象集随机选择k个元素作为初始的聚类中心；然后根据其他不同元素与聚类中心的距离，将其划分进最相似的簇中；接着重新更新各聚类中心；循环上述过程直至准则函数收敛。
进一步地，步骤二中窗函数要求主瓣无限窄且没有旁瓣，从而避免频谱泄露。
进一步地，所述步骤三中各帧的希尔伯特谱和边际谱为：
h(ω)＝∫h(ω,t)dω；
其中，h(ω,t)为希尔伯特谱，h(ω)为希尔伯特边际谱，h为希尔伯特变换运算符，re表示取实部运算；
进一步地，所述步骤四中mel三角带通滤波器组包含m个带通滤波器hm(k)，并且均为三角形滤波器。
进一步地，所述步骤七的聚类中心由以下公式求取：聚类样本为x＝{xi∈rd,i＝1,2,3,…,n},划分成不同类c＝{c1,c2,…,ck}，xi为d为向量，k为聚类的个数；
停止准则函数为：
其中d(xi,cj)表示各元素xi与聚类中心cj的距离；j表示簇内距的和。
与现有技术相比本发明的优点在于：充分考虑了地震信号的非平稳性，采用hht分解地震信号，求取的地震信号分量的希尔伯特边际谱，时频分辨率更高且物理意义更加明确。基于希尔伯特边际谱采用mfcc方法提取地震信号的地震声纹特征参数对储层含流体具有较强的标识性以及更好地容错性，能够更加有效且准确的反映储层流体性质，有利于提高储层流体预测的精度。
具体实施方式
为使本发明的目的、技术方案及优点更加清楚明白，以下举实施例，对本发明做进一步详细说明。
实施例1
一种基于hht-mfcc的地震声纹流体预测方法，包括以下步骤：
步骤一：输入叠后地震单道cdp数据，对输入数据进行emd经验模态分解，得到不同频带范围的imf分量:
上式中s(t)为输入的叠后地震单道cdp数据，imfi(t)为经验模态分解得到的各imf分量，rn(t)为分解最后的残余趋势项。
步骤二：对各imf分量分别进行分帧加窗处理；采用有限长度的时窗来截取各imf分量，将其划分成一个一个的短时段信号即一帧。为了保证各帧之间的连续性，相邻帧之间保留一段重叠区域，重叠区域长度取窗函数的1/3-1/2窗长。窗函数要求主瓣无限窄且没有旁瓣，从而避免频谱泄露。本发明采用的是汉明窗，其定义如下：
其中w(n)为窗函数，n为窗长。
步骤三：对分帧加窗后的各imf分量逐帧进行希尔伯特变换得到各帧的希尔伯特边际谱，将各imf分量从时域转换到频域。
各帧的希尔伯特谱和边际谱为：
h(ω)＝∫h(ω,t)dω；
其中，h(ω,t)为希尔伯特谱，h(ω)为希尔伯特边际谱，h为希尔伯特变换运算符，re表示取实部运算；
步骤四：将步骤三得到的频率域信号通过m个mel三角带通滤波器组获得mel频谱，再将得到的mel频谱取对数能量，得到相应频带的对数能量e(m)。其中mel三角带通滤波器组包含m个带通滤波器hm(k)，并且均为三角形滤波器。
mel三角带通滤波器hm(k)的中心频率f(m)可由下式求取：
其中，fs为地震信号的采样频率，fh为带通滤波器的最高频率，一般取fs的一半。fl为带通滤波器的最低频率，一般取值趋于零值。m表示mel三角带通滤波器组中的第m个带通滤波器，n为步骤二中的窗长。
每个滤波器的输出的对数能量e(m)为：
e(m)＝log10(h2(f)hm(f))；
其中h(f)为信号希尔伯特边际谱，hm(f)为三角带通滤波器的频率响应；
步骤五：对步骤四得到的对数能量e(m)进行离散余弦变换(dct)，去除e(m)间的相关性，得到m阶mel频率倒谱系数(mfcc):
步骤六：对输入地震的每个cdp信号循环步骤一到步骤五，即可得到原始二维地震信号的m阶地震声纹参数。
步骤七：对步骤六的结果进行k均值聚类分析：k均值聚类的基本思想是将包含n个元素的对象集划分成k个不同的类(簇)，期望簇内元素相似度最高，簇间的相似度最低。算法的实现过程如下，首先从包含n个元素的对象集随机选择k个元素作为初始的聚类中心；然后根据其他不同元素与聚类中心的距离，将其划分进最相似的簇中；接着重新更新各聚类中心；循环上述过程直至准则函数收敛。
聚类中心由以下公式求取：聚类样本为x＝{xi∈rd,i＝1,2,3,…,n,划分成不同类c＝c1,c2,…,ck，xi为d为向量，k为聚类的个数。
停止准则函数为：
其中d(xi,cj)表示各元素xi与聚类中心cj的距离；j表示簇内距的和。
实施例2
本实施例只针对实施例1的不同之处进行描述，相同之处不再阐述；
本方法步骤三中可以通过快速傅里叶变换(fft)对分帧加窗后的信号进行时间域到频率域的转换，得到各帧信号的傅里叶谱。然而傅里叶谱表示频率全局性的分布，对于非线性非平稳信号，傅里叶谱缺乏实际的物理意义。本方法采取的希尔伯特边际谱则反映了不同频率点上幅值的分布，表示不同频率在信号整体上的幅度贡献。相较于傅里叶谱，本方法采用的希尔伯特谱具有更高的时频分辨率。
本领域的普通技术人员将会意识到，这里所述的实施例是为了帮助读者理解本发明的实施方法，应被理解为本发明的保护范围并不局限于这样的特别陈述和实施例。本领域的普通技术人员可以根据本发明公开的这些技术启示做出各种不脱离本发明实质的其它各种具体变形和组合，这些变形和组合仍然在本发明的保护范围内。