一种鲁棒的地震信号瞬时频率的估计方法技术

一种鲁棒的地震信号瞬时频率的估计方法属于地震信号处理技术领域。其特征在于该方法包括下列步骤：利用地震勘探设备，采集地震数据；获得叠后地震数据；根据每道地震数据的长度，设定窗长和加权阶数；逐点取一段窗长的观测信号，进行加窗处理，计算该信号对应的解析信号的频率域响应；再对其进行振幅的几何级数加权，然后利用快速傅里叶逆变换得到新的解析信号；每个时间点会得到多段解析信号的估计值，利用它们的均值作为该点对应的新的解析信号；利用新的解析信号计算信号的瞬时相位，再计算瞬时频率。本发明专利技术不仅可以准确地计算地震信号的瞬时频率，而且信噪比低的情况下依然能够准确的估计瞬时频率。

【技术实现步骤摘要】

本专利技术涉及信号处理领域，具体地说涉及到地震信号处理领域。
技术介绍
地震勘探是一种利用人工地震技术激发地震波，然后通过地面布置的检波器得到由地下地层分界面反射回来的声波，通过对数据进行处理，提取有用的地震属性以查明地下的地质构造或岩性特征，为寻找油气藏或其他勘探目的提供重要的工具。 地震信号的瞬时频率估计是地震信号的属性提取中的关键问题之一，是其他一些属性提取的基础。地震信号的瞬时频率估计在识别岩性和储层分析方面具有重要的作用，因此准确估计地震信号的瞬时频率具有重要的应用价值。 目前，常用的瞬时频率估计方法是基于希尔伯特变换的解析信号方法(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.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的二次调频信本文档来自技高网...

【技术保护点】
一种鲁棒的地震信号瞬时频率的估计方法，其特征在于，该方法依次含有以下步骤步骤（１）利用地震勘探设备，采集地震数据：在反射波地震勘探中采用一点放炮多点接收的方法，在预勘探区域设置多个接收点，所述多个接收点散布在包括炮点在内的地面上一个勘探目标区域范围的矩形网格点上，使用地震检波设备获取爆炸波及地震波的波形；步骤（２）对步骤（１）中获取的波形数据进行地震资料的常规处理，经过反褶积、抽取共中心点道集、动校正、叠加和偏移过程之后获得叠后的每道数据ｘ（ｔ），ｔ＝１，２，…，Ｈ，采样时间间隔Ｔ，Ｔ为１ｍｓ、２ｍｓ或４ｍｓ，设定窗长Ｌ和加权阶数Ｎ；步骤（３）计算机按以下步骤计算新的解析信号＊（ｔ）：步骤（３．１）沿时间点ｉ＝１，２，…，Ｈ－Ｌ＋１逐点取出一段窗长为Ｌ的观测信号ｚ↓［ｉ］＝｛ｘ（ｉ＋ｊ－１），ｊ＝１，２，…，Ｌ｝，进行加窗处理ｙ↓［ｉ］＝ｚ↓［ｉ］.ｈ↓［ｉ］，其中ｈ↓［ｉ］为长度Ｌ的高斯窗函数，求出该信号的频率域响应Ｘ（＊），Ｘ↓［ｉ］（＊）＝ＦＦＴ（ｘ↓［ｉ］（ｔ）），其中ＦＦＴ为快速傅立叶变换；步骤（３．２）按下述公式计算该信号对应的解析信号的频率域响应＊＊＊步骤（３．３）变换Ｆ↓［ｉ］（＊）为Ｆ↓［ｉ］（＊）＝Ａ（＊）ｅｘｐ（ｉθ（＊）），利用加权振幅几何级数的方法计算新解析信号的频率域函数＊↓［ｉ］（＊），＊↓［ｉ］（＊）＝Ｆ↓［ｉ］（＊）Ａ↑［Ｎ－１］（＊）／ｍａｘ（Ａ↑［Ｎ－１］（＊）），其中Ｎ为加权阶数，ｍａｘ（.）为取最大值函数；步骤（３．４）利用快速傅里叶逆变换计算ｘ↓［ｉ］（ｔ）对应的新的解析信号＊↓［ｉ］（ｔ），＊↓［ｉ］（ｔ）＝ＩＦＦＴ（＊↓［ｉ］（＊）），ＩＦＦＴ为快速傅立叶逆变换；步骤（３．５）对于时间点ｉ，得到多个解析信号的估计＊↓［ｋ］（ｉ），ｋ＝ｉ－Ｌ＋１，ｉ－Ｌ＋２，…，ｉ，将该点对应多个估计的均值作为该点对应的新的解析信号＊（ｉ）＝１／ｎ＊＊↓［ｋ］（ｉ），其中ｎ为时间点ｉ所估计出的解析信号的个数，＊（.）表示求和；步骤（４）利用新的解析信号＊（ｔ），该解析信号表示为＊（ｔ）＝＊（ｔ）＋ｉ＊（ｔ），按下面的公式求出地震信号的瞬时相位θ（ｔ），θ（ｔ）＝ｔａｎ↑［－１］（＊（ｔ）／＊（ｔ）），ｔａｎ↑［－１］（.）表示反正切函数；步骤（５）相位展开，如｜θ（ｉ＋１）－θ（ｉ）｜＞σ，θ（ｉ＋１）＝θ（ｉ＋１）＋２πｋ，使得相邻两个相位之间的差值不超过σ，其中σ是阈值，取为π，ｋ...

【技术特征摘要】

【专利技术属性】
技术研发人员：陆文凯，张长开，
申请(专利权)人：清华大学，
类型：发明
国别省市：11[中国|北京]