专利名称：一种鲁棒的地震信号瞬时频率的估计方法
技术领域：
本发明涉及信号处理领域，具体地说涉及到地震信号处理领域。
背景技术：
地震勘探是一种利用人工地震技术激发地震波，然后通过地面布置的检波器得到由地下地层分界面反射回来的声波，通过对数据进行处理，提取有用的地震属性以查明地下的地质构造或岩性特征，为寻找油气藏或其他勘探目的提供重要的工具。
地震信号的瞬时频率估计是地震信号的属性提取中的关键问题之一，是其他一些属性提取的基础。地震信号的瞬时频率估计在识别岩性和储层分析方面具有重要的作用，因此准确估计地震信号的瞬时频率具有重要的应用价值。
目前，常用的瞬时频率估计方法是基于希尔伯特变换的解析信号方法(B.Boashash，“Estimating and interpreting the instantaneous frequency of asignal-PartlFundamentals，”Proc.IEEE，vol.80，no.4，pp.520-538，1992.4)。该方法介绍如下解析信号(又称复信号)由实部和虚部组成。它的实部是原信号，虚部是原信号的希尔伯特变换。对于一个普通的信号x(t)＝A(t)cos(θ(t))，其中A(t)是振幅，θ(t)是相位，则其希尔伯特变换为
根据解析信号的定义，解析信号f(t)可以表示为
利用欧拉公式可以得到f(t)＝A(t)exp(iθ(t))，其中exp(·)表示自然指数函数(下同)，i为虚数单位。瞬时频率S(t)定义为相位的时间变化率
(J.Ville，“Theorie etapplication de la notion de signal analytic”，Cables et Transmissions，vol.2A，pp.61-74，1948)。至此就是利用普通的希尔伯特变换的方法估计信号的瞬时频率的全部过程。
该方法在信噪比较高的环境下可以得到可信度较高的结果，但是在信噪比较低的环境下它的结果受噪声的影响很严重。在努力降低噪声对结果带来的影响方面，Luo Yi曾经用一种称为广义希尔伯特变换的方法来估计信号的瞬时相位，从而更好的显示河道和断裂带(Luo Yi.Seismic data processing method toenhance fault and channel displayUS006963804[P].2005-11-08)，这种技术对噪声环境下瞬时相位的估计有一定的效果。
针对上述现有技术存在的问题，本发明的目的是提出一种鲁棒的地震信号瞬时频率的估计方法。使用该方法能够信噪比低的情况下能够准确的估计出信号的瞬时频率。
发明内容
本发明采用的技术方案是一种鲁棒的地震信号瞬时频率的估计方法，该方法包括以下步骤 步骤(1)利用地震勘探设备，采集地震数据 在反射波地震勘探中采用一点放炮多点接收的常规方法，在预勘探区域设置多个接收点，所述多个接收点散布在包括炮点在内的地面上一个勘探目标区域范围的矩形网格点上，使用地震检波设备获取爆炸波及地震波的波形； 步骤(2)对步骤(1)中获取的波形数据进行地震资料的处理，经过反褶积、抽取共中心点道集、动校正、叠加和偏移过程之后获得叠后的每道数据x(t)，t＝1，2，…，H，采样时间间隔T，T为1ms、2ms或4ms，设定 窗长L和加权阶数N，其中L取20～30，N取2～5； 步骤(3)计算机按以下步骤计算新的解析信号
步骤(3.1)沿时间点i＝1，2，…，H-L+1逐点取出一段窗长为L的观测信号zi＝{x(i+j-1)，j＝1，2，…，L}，进行加窗处理yi＝zi·hi，其中hi为长度L的高斯窗函数，求出该信号的频率域响应
其中FFT为快速傅立叶变换； 步骤(3.2)按下述公式计算该信号对应的解析信号的频率域响应
步骤(3.3)变换
为
利用加权振幅几何级数的方法计算新解析信号的频率域函数
其中N为加权阶数，max(·)为取最大值函数； 步骤(3.4)利用快速傅里叶逆变换计算xi(t)对应的新的解析信号
IFFT为快速傅立叶逆变换； 步骤(3.5)对于时间点i，得到多个解析信号的估计
k＝i-L+1，i-L+2，…，i，将该点对应多个估计的均值作为该点对应的新的解析信号
其中n为时间点i所估计出的解析信号的个数，
表示求和； 步骤(4)利用新的解析信号
该解析信号表示为
按下面的公式求出地震信号的瞬时相位θ(t)，
tan-1(·)表示反正切函数； 步骤(5)相位展开，如果|θ(i+1)-θ(i)|＞σ，θ(i+1)＝θ(i+1)+2πk，使得相邻两个相位之间的差值不超过σ，其中σ是阈值，取为π，k为任意整数； 步骤(6)利用估计出来的瞬时相位，利用相位差分方法来计算地震信号的瞬时频率S(t)，采用下面的这个公式
其中T为采样 时间间隔，单位是秒。
由上述本发明采用的技术方案可以看出，本发明是在地震观测信号中逐点选取一定长度的信号进行处理，然后再估计整个观测信号对应的解析信号；在解析信号的频率域对振幅进行了几何级数的加权，这种方法很好的压制了噪声；对于在每个时间点上得到的多个估计，利用多个估计的均值作为该点的解析信号的值，使得估计出来的新的解析信号更加的准确。这样通过解析信号估计地震信号的瞬时相位进而估计地震信号的瞬时频率将更加的准确。
图1是本发明方法的流程示意图。
在我们的仿真实验中，分别产生两种类型的信号，第一种是普通的正弦信号，分别在有噪声和没有噪声的情况下；第二种是普通的chirp(调频)信号，分别在有噪声和没有噪声的情况下。在比较中，分别给出了原信号，原信号真实的瞬时频率，用普通的希尔伯特变换得到的瞬时频率，以及用本发明方法得到的瞬时频率。
图2、图3、图4、图5分别给出了不同情况下两种方法的结果。可以看出在没有噪声的情况下，本方法和普通的方法一样能准确的估计出信号的瞬时频率；在有噪声的情况下，普通的方法无法有效的估计出信号的瞬时频率，但是本发明方法却能够准确的估计出信号的瞬时频率。
图1、本发明所属方法的计算机程序流程图。
图2、正弦信号在没噪声情况下的瞬时频率估计 a，瞬时频率为100Hz的正弦信号； b，正弦信号的真实瞬时频率； c，利用普通希尔伯特变换得到的瞬时频率； d，利用本发明方法得到的瞬时频率。
图3、正弦信号在有噪声情况下的瞬时频率估计 a，添加了15dB高斯白噪声的瞬时频率为100Hz的正弦信号； b，正弦信号的真实瞬时频率； c，利用普通希尔伯特变换得到的瞬时频率； d，利用本发明方法得到的瞬时频率。
图4、调频信号在无噪声情况下的瞬时频率估计 a，瞬时频率从100Hz到500Hz的二次调频信号； b，调频信号的真实瞬时频率； c，利用普通希尔伯特变换得到的瞬时频率； d，利用本发明方法得到的瞬时频率。
图5、调频信号在有噪声情况下的瞬时频率估计 a，添加了15dB高斯白噪声的瞬时频率从100Hz到500Hz的二次调频信号； b，正弦信号的真实瞬时频率； c，利用普通希尔伯特变换得到的瞬时频率； d，利用本发明方法得到的瞬时频率。
具体实施例方式 普通的希尔伯特变换是将一个信号进行-90°的相移。可以将其理解为将一个信号通过一个全通滤波器
得到该信号的希尔伯特变换。该滤波器具有如下特性
综合起来即为
对于一个观测信号对应的解析信号的频率域响应，可以通过如下方法来实现给定一个信号x(t)，其频率域响应为
它的解析信号表示为
其中
为它的希尔伯特变换，i为虚数单位；它的频域表达式为
而
综合起来就是
即我们对该观测信号先求解它的快速傅里叶变换，然后将正频率的幅值乘以2，负频率的幅值置0即可得到该观测信号对应的解析信号的频率域响应。
具体来说，本发明的方法是在一个范围内逐点依次对整个信号中的每一小段地震信号进行解析信号的估计，然后再估计整段信号的解析信号。对于获得的地震观测信号x(t)，t＝1，2，…，H，沿时间点i＝1，2，…，H-L+1依次取窗长L＝20的观测信号zi(t)＝{x(i)，x(i+1)，…，x(i+L-1)}，进行加窗处理yi＝zi·hi，其中hi为长度L的高斯窗函数，
计算yi(t)的频率域响应
利用(1)式求出该信号对应的解析信号的频率域响应
变换该频域函数
其中
是振幅谱，
是相位谱。对获得的频域函数进行振幅加权处理
其中加权阶数N＝3。新获得的频率域函数
为信号xi(t)对应的新的解析信号的频率域响应，对其进行傅里叶逆变换即可得到xi(t)对应的新的解析信号
其中IFFT为快速傅里叶逆变换。由于是逐点选取长度为L的观测信号的，因此对于每个时间点i上会出现多个解析信号的估计
其中k一般为i-L+1，i-L+2，…，i，对于多个估计，利用它们的均值作为该点上的解析信号的值 n为每个点上获得的估计信号的个数。获得了新的解析信号
之后，可以将其写成
其中
在获得信号的瞬时相位之后，由于求出的相位的范围在(-π，π]之间，相邻的两个相位会出现跳变，使得相邻的相位不连续，因此需要对相位进行展开，使得相邻的两个相位之间的差值不超过一个固定的阈值σ(J.M.Tribolet，a new phase unwrapping algorithm，IEEE Transactionson acoustics，speech，and signal processing，vol.ASSP-23，No.2，1977.4)，这里阈值σ＝π。
目前在利用解析信号计算其瞬时频率方面，出现了多种计算方法(E.Barnes，The calculation of instantaneous frequency and instantaneous bandwidthGeophysics，57，1992.11)，有基于两点的有限脉冲响应的微分器(Scheuer.T.E，Oldenburg.D.W，Local phase velocity from complex seismic data，Geophysics，1988，53pp.1503)，有基于三点的有限脉冲响应的微分器(Boashash.B，O’Shea.P，Ristic.B，Statisticalcomputation comparison of some estimators for instantaneous frequency，Proc IEEEICASSP-91，1991，3193～3196)，相位差分法包括前向有限差分法、后向有限差分法、中心有限差分法(Boashash.B，Estimating and interpreting the instantaneousfrequency of a signal-Part2algorithms and applications.Proc.IEEE，vol.80，no.4，pp.540-568，1992.4)。在本发明里利用相位后向差分法计算信号的瞬时频率S(t)，即利用求出的相位然后进行后向差分
普通的希尔伯特变换是对整个信号进行变换的，在频域没有做任何处理；而本发明的方法是先逐点选择一小段信号，然后对每一段信号进行如下处理先进行加窗处理，然后进行频率域的振幅几何级数加权，接着利用傅里叶逆变换获得新的解析信号。再进行求均值获得鲁棒的估计结果。该方法的好处是由于噪声在频率域表现为一些很小的值，这样利用振幅几何级数加权时可以大大的压制噪声；同时反变换后对同一个点上多个值进行求平均，因此这样得到的结果将会更加的鲁棒，受噪声的干扰小。
总之，本发明适用于地震信号的瞬时频率的估计，不仅可以在高信噪比的环境中准确的估计出信号的瞬时频率，在低信噪比的环境下仍然能准确的估计出信号的瞬时频率。
权利要求
1.一种鲁棒的地震信号瞬时频率的估计方法，其特征在于，该方法依次含有以下步骤
步骤(1)利用地震勘探设备，采集地震数据
在反射波地震勘探中采用一点放炮多点接收的方法，在预勘探区域设置多个接收点，所述多个接收点散布在包括炮点在内的地面上一个勘探目标区域范围的矩形网格点上，使用地震检波设备获取爆炸波及地震波的波形；
步骤(2)对步骤(1)中获取的波形数据进行地震资料的常规处理，经过反褶积、抽取共中心点道集、动校正、叠加和偏移过程之后获得叠后的每道数据x(t)，t＝1，2，…，H，采样时间间隔T，T为1ms、2ms或4ms，设定窗长L和加权阶数N；
步骤(3)计算机按以下步骤计算新的解析信号
步骤(3.1)沿时间点i＝1，2，…，H-L+1逐点取出一段窗长为L的观测信号zi＝{x(i+j-1)，j＝1，2，…，L}，进行加窗处理yi＝zi·hi，其中hi为长度L的高斯窗函数，求出该信号的频率域响应
其中FFT为快速傅立叶变换；
步骤(3.2)按下述公式计算该信号对应的解析信号的频率域响应
步骤(3.3)变换
为
利用加权振幅几何级数的方法计算新解析信号的频率域函数
其中N为加权阶数，max(·)为取最大值函数；
步骤(3.4)利用快速傅里叶逆变换计算xi(t)对应的新的解析信号
IFFT为快速傅立叶逆变换；
步骤(3.5)对于时间点i，得到多个解析信号的估计
k＝i-L+1，i-L+2，…，i，将该点对应多个估计的均值作为该点对应的新的解析信号
其中n为时间点i所估计出的解析信号的个数，
表示求和；
步骤(4)利用新的解析信号
该解析信号表示为
按下面的公式求出地震信号的瞬时相位θ(t)，
tan-1(·)表示反正切函数；
步骤(5)相位展开，如|θ(i+1)-θ(i)|＞σ，θ(i+1)＝θ(i+1)+2πk，使得相邻两个相位之间的差值不超过σ，其中σ是阈值，取为π，k为任意整数；
步骤(6)利用估计出来的瞬时相位，利用相位差分方法来计算地震信号的瞬时频率S(t)，采用下面的这个公式
其中T为采样时间间隔，单位是秒。
2.根据权利要求1所述的一种鲁棒的瞬时频率的估计方法，其特征在于，加权阶数N为2～5之间的整数。
3.根据权利要求1所述的一种鲁棒的瞬时频率的估计方法，其特征在于，滤波窗长L根据采样时间间隔来选取，L为20～30。
4.根据权利要求1所述的一种鲁棒的瞬时频率的估计方法，其特征在于，所选择的窗函数是矩形窗、高斯窗、汉宁窗、汉明窗、三角窗之一。
全文摘要
一种鲁棒的地震信号瞬时频率的估计方法属于地震信号处理技术领域。其特征在于该方法包括下列步骤利用地震勘探设备，采集地震数据；获得叠后地震数据；根据每道地震数据的长度，设定窗长和加权阶数；逐点取一段窗长的观测信号，进行加窗处理，计算该信号对应的解析信号的频率域响应；再对其进行振幅的几何级数加权，然后利用快速傅里叶逆变换得到新的解析信号；每个时间点会得到多段解析信号的估计值，利用它们的均值作为该点对应的新的解析信号；利用新的解析信号计算信号的瞬时相位，再计算瞬时频率。本发明不仅可以准确地计算地震信号的瞬时频率，而且信噪比低的情况下依然能够准确的估计瞬时频率。
文档编号G01V1/28GK101825722SQ20101013335
公开日2010年9月8日 申请日期2010年3月25日 优先权日2010年3月25日
发明者陆文凯, 张长开 申请人:清华大学