CN103645504A - 基于广义瞬时相位及p范数负模的地震弱信号处理方法 - Google Patents

Abstract

本发明涉及一种基于广义瞬时相位及P范数负模的地震弱信号处理方法，其包括以下内容：1）对地震信号进行识别，具体过程为：对输入的地震信号进行P阶广义希尔伯特变换并得到广义瞬时相位；求取广义瞬时相位P范数负模；根据广义瞬时相位P范数负模提取地震弱信号的相关属性参数，通过提取的属性参数对地震有效信号进行识别；2）对地震信号的噪声进行去除，得到信噪比比较高的地震信号，具体过程为：对输入的地震信号进行P阶广义希尔伯特变换并得到广义瞬时相位；对计算得到的广义瞬时相位进行平滑滤波处理，得到滤波后的广义瞬时相位信号；采用逐步迭代的方法求取地震信号的瞬时振幅；根据求取的广义瞬时相位信号和瞬时振幅重构地震信号；并对重构地震信号的合理性进行判断。

Description

基于广义瞬时相位及P范数负模的地震弱信号处理方法

技术领域

本发明涉及一种油气地震勘探领域中地震资料的弱振幅识别和去噪方法，特别是关于一种适用于为后续地震资料岩性反演和高精度解释提供高信噪比数据的基于广义瞬时相位及P范数负模的地震弱信号处理方法。

背景技术

在油气地球物理勘探领域，地震资料的“三高”（高分辨率、高信噪比、高保真度）处理一直是研究的难点和热点，特别是在岩性反演广泛应用于隐蔽型油气藏预测的今天。其中，地震资料的信噪比水平在很大程度上制约着分辨率的提高（根据分辨率定量估算公式，分辨率高低主要与资料带宽和信噪比有关，信噪比对分辨率的贡献权重为2/5）。另一方面，地震振幅包含了大量的储层信息，噪音的存在干扰了地震振幅的有效识别，由此确定的储层圈闭和井位的可信度也受到很大影响，进而影响到钻井的成功率。因此，强噪背景下的弱信号识别和提取在资料处理中起着重要的作用，也是应用地球物理学和自然科学的研究热点。

目前去噪方法有很多种，由于适用条件不一样，不同类型的噪声需要选取合适的去噪方法进行滤除：其中一种去噪方法是滤波去噪，从利用有效信号和噪声不同频带分布的经典滤波，到利用二者统计特性区别的维纳滤波理论，都能够在某种程度上对信号进行去噪处理，但是对信号本身的要求比较高；另一种信号去噪方法是时频分析方法，能够对非平稳地震信号进行去噪处理，如小波变换，引入尺度因子和平移因子，实现了信号特征的局域化分析，是信号去噪处理的有力工具之一。希尔伯特变换是信号分析的另一个有力工具，在地震信号处理和解释中应用了很多年，但是存在对噪声太敏感（即抗噪能力较弱）的缺点，Yi Luo等人在2003年提出了广义希尔伯特变换（GHT），这种变换的抗噪性较好，对噪声不像其他时频分析方法那样有明显的敏感性，能够有效的探测地震剖面中的突变信息，为弱信号的识别和提取研究提供了一个较好的工具，根据GHT提取的信号瞬时属性（包括瞬时振幅、瞬时频率和瞬时相位），一般称为广义瞬时属性，广义瞬时属性高鲁棒性这种特征目前还没有被引入到弱信号识别和去噪处理中。

发明内容

针对上述问题，本发明的目的是提供一种能够对地震弱信号进行高精度识别和去噪处理，且能够为后续地震资料的处理和解释特别是岩性反演提供高品质资料的基于广义瞬时相位及P范数负模的弱信号处理方法。

为实现上述目的，本发明采取以下技术方案：一种基于广义瞬时相位及P范数负模的地震弱信号处理方法，其包括以下内容：1）对地震信号进行识别，具体过程为：1.1）对输入的地震信号x(n)进行P阶广义希尔伯特变换并得到广义瞬时相位φ(t)；1.2）求取广义瞬时相位P范数负模γ(t)；1.3）根据广义瞬时相位P范数负模γ(t)提取地震弱信号的相关属性参数，通过提取的属性参数对地震有效信号进行识别；2）对地震信号的噪声进行去除，得到信噪比比较高的地震信号，具体过程为：2.1）对输入的地震信号x(n)进行P阶广义希尔伯特变换并得到广义瞬时相位φ(t)；2.2）对计算得到的广义瞬时相位φ(t)进行平滑滤波处理，得到滤波后的广义瞬时相位信号

2.3）采用逐步迭代的方法求取地震信号的瞬时振幅a(n)；2.4）根据求取的广义瞬时相位信号

和瞬时振幅a(n)重构地震信号2.5）将重构地震信号

与测井资料得到的数据进行比对，对重构地震信号

的合理性进行判断，如果合理则输出重构信号

如果不合理，则返回步骤2.2）。

广义瞬时相位

式中，hr(t)和hi(t)是复数道h(t)的实部和虚部。

广义瞬时相位P范数负模γ(t)：

$\mathrm{\γ}\left(t\right)=-\left|\arctan \frac{\mathrm{hi}\left(t\right)}{\mathrm{hr}\left(t\right)}\right|$

式中，hr(t)和hi(t)是复数道h(t)的实部和虚部。

所述步骤2.3）采用逐步迭代的方法求取地震信号的瞬时振幅a(n)，包括以下步骤：

①将地震信号表示为广义瞬时相位信号

的形式，即：

②采用一假设的正数c代替①中振幅的初始值迭代值|a₀(n)|，令

③计算

的N点DFT得到A₁(k)，由于解析信号的傅立叶变换是的单边函数，在

处取值为0，所以有 $A_1^{\′}\left(k\right)=\left\{\begin{array}{cc}{A}_1\left(k\right)& 0\≤\frac{N}{2}-1\\ 0& \frac{N}{2}\≤k\≤N-1\end{array}\right.,$ 式中N都取2^M形式；

④求A′₁(k)的傅立叶逆变换IDFT得到

其中|a′₁(n)|是a(n)的新估算值；

⑤判断|a′₁(n)|是否满足精度要求，如果不满足精度要求，则返回步骤③继续迭代；假设迭代次数为p，则迭代时采用

代替

迭代结果如果满足精度要求，则输出得到结果a(n)=|a'_p(n)|。

本发明由于采取以上技术方案，其具有以下优点：1、本发明的地震弱信号识别方法是基于广义瞬时相位负模，根据广义瞬时相位范数负模提取地震弱信号的相关属性参数，通过得到的属性参数对地震有效信号进行识别和应用，因此能够对强噪背景下的弱信号进行较好的识别，具有较好的鲁棒性。2、本发明的地震去噪方法首先通过P阶广义希尔伯特变换得到广义瞬时相位并对其进行滤波处理，然后采用逐步迭代的方法求取地震信号的瞬时振幅，最后重构地震信号，因此本发明能够较好地滤除随机噪音，并能够对两个过程对噪音进行压制：1）广义瞬时相位的提取过程具有自适应压噪；2）广义瞬时相位曲线的平滑滤波，去除噪音引起的曲线振荡。3、由于地震振幅包含了大量的储层信息，噪音的存在干扰了地震振幅的有效识别，因此本发明通过迭代的方法得到比较接近真实的振幅，因此通过此振幅能够确定储层圈闭和井位的可信度，进一步提高钻井的成功率。本发明可以广泛应用于地震资料的分析处理中。

附图说明

以下结合附图来对本发明进行详细的描绘。然而应当理解，附图的提供仅为了更好地理解本发明，它们不应该理解成对本发明的限制。

图1为本发明的瞬时相位提取示意图，图1（A）为无噪信号，图1（B）为含噪信号；

图2为本发明的P范数负模提取示意图，图2（A）为无噪信号，图2（B）为含噪信号；

图3（A）为地震合成记录（含单层、双层和三层模型）在含噪情况下广义瞬时相位P范数负模的提取图例，图3（B）是无噪合成记录和相位负模的对比示意图；

图4为本发明基于广义瞬时相位重构信号的迭代流程示意图；

图5（A）为无噪信号基于广义瞬时相位的重构图，图5（B）为含噪信号的去噪效果图；

图6（A）为合成的无噪CMP道集示意图，图6（B）为加噪的CMP道集示意图，图6（C）为去噪后道集示意图，图6（D）为滤除的噪音示意图；

图7（A）为去噪前的时频谱示意图，图7（B）为去噪后的时频谱示意图；

图8（A）为含噪地震剖面示意图，图8（B）为去噪后剖面示意图，图8C为滤除的噪音示意图。

具体实施方式

下面结合附图和实施例对本发明进行详细的描述。

本发明的地震信号处理方法是基于Li You提出的广义希尔伯特变换（GHT），下面从传统的希尔伯特变换（HT）出发，引出GHT的定义。

对于给定的连续时间信号x(t)，希尔伯特变换H[x(t)]定义为：

$H\left[x\left(t\right)\right]=x\left(t\right)*\frac{1}{\mathrm{\πt}}=\frac{1}{\mathrm{\π}}{\∫}_{-\∞}^{\∞}\frac{g\left(\mathrm{\τ}\right)}{t-\mathrm{\τ}}\mathrm{d\τ}={\∫}_{-\∞}^{\∞}\frac{g(t-\mathrm{\τ})}{\mathrm{\τ}}\mathrm{d\τ}---\left(1\right)$

式中，x(t)的希尔伯特变换是x(t)和信号1/πt的褶积，它是信号x(t)对线性时不变滤波器（称为希尔伯特变换器）的响应，该滤波器的单位脉冲响应为h(t)＝1/πt。希尔伯特变换H[x(t)]通常记为

由傅立叶变换理论可以得到，jh(t)＝j/πt的傅立叶变换是符号函数sgn(Ω)，因此希尔伯特变换器的频率响应表示为：

$H\left(\mathrm{j\&Omega;}\right)=-\mathrm{jsgn}\left(\mathrm{\&Omega;}\right)=\left\{\begin{array}{cc}-j& \mathrm{\&Omega;}>0\\ 0& \mathrm{\&Omega;}=0\\ j& \mathrm{\&Omega;}<0\end{array}\right.---\left(2\right)$

传统的希尔伯特变换在频率域可以表示为：

hi(ω)＝X(ω)sgn(ω)i    （3）

式中，X(ω)是信号x(t)的傅立叶变换，hi(ω)是x(t)的在频率域的希尔伯特变换，sgn(ω)是符号函数，其中

对hi(ω)应用傅立叶逆变换并令t＝0，得到：

hr(0)＝2Σ_ωRe[X(ω)]+Re[X(0)]

hi(0)＝2Σ_ωIm[X(ω)]    （4）

h(0)＝hr(0)+i·hi(0)

式中，Σ_ω是正频率的和，hr(0)和hi(0)是复数道h(0)的实部和虚部。

为了使上式适用于所有时间值t，广义希尔伯特变换（GHT）定义为：

$\begin{array}{c}\mathrm{hr}\left(t\right)={\{2{\mathrm{\&Sigma;}}_{\mathrm{\&omega;}}{\left\{\mathrm{Re}\right[X(t,\mathrm{\&omega;})\left]\right\}}^n+\mathrm{Re}{\left[X(t,0)\right]}^n\}}^{\frac{1}{n}}\\ \mathrm{hi}\left(t\right)={\left\{2{\mathrm{\&Sigma;}}_{\mathrm{\&omega;}}{\left\{\mathrm{Im}\right[X(t,\mathrm{\&omega;})\left]\right\}}^n\right\}}^{\frac{1}{n}}\\ h\left(t\right)=\mathrm{hr}\left(t\right)+i\&CenterDot;\mathrm{hi}\left(t\right)\end{array}---\left(5\right)$

式中，hi(t)是输入信号的Lⁿ阶广义希尔伯特变换，h(t)是复数道，hr(t)和hi(t)是复数道h(t)的实部和虚部。X(t,ω)是输入信号在以t为中心的时窗内的加窗傅立叶变换。当采用高斯窗时，X(t,ω)就指的是信号的STFT。所以从某种意义上也可以说，广义希尔伯特变换是短时傅立叶变换的推广。

由GHT的定义可以看出，当阶次n＝1并且窗函数是无限长时，GHT和传统HT的结果是相同的；当阶次n＝2时，GHT的包络就是常规的均方根振幅，与均方根振幅相关的相位也能够从2阶GHT中衍生出来。广义希尔伯特变换是对传统希尔伯特变换的扩展，它不仅具有HT的功能，还具有HT所不具备的自适应压噪功能，这和GHT的基本原理是相对应的。

广义希尔伯特变换的基本思想是：在频域对实信号谱和希尔伯特变换后的谱与傅立叶核函数乘积后做提反差处理，即给每个采样值进行n次方计算，使能量大小的差距进一步拉大，再通过傅立叶逆变换的积分（求和）运算的统计效应压噪，这种压噪方法由于统计效应的应用而具有自适应压噪的功能。

这里给出与本发明有关的两个瞬时相位定义：广义瞬时相位及其P范数负模，其中广义瞬时相位为：

$\mathrm{\&phi;}\left(t\right)=\arctan \frac{\mathrm{hi}\left(t\right)}{\mathrm{hr}\left(t\right)}---\left(6\right)$

广义瞬时相位P范数负模γ(t)：

$\mathrm{\&gamma;}\left(t\right)=-\left|\arctan \frac{\mathrm{hi}\left(t\right)}{\mathrm{hr}\left(t\right)}\right|---\left(7\right)$

式中，P是广义希尔伯特变换取的阶次n，P范数的概率来自于GHT实部和虚部的提取与数列P范数的求取过程类似，γ(t)称为信号的广义瞬时相位P范数负模。

信号在时域表示为f(t)，在频域表示为F(ω)，其中F(ω)是f(t)的傅立叶变换，二者是完全等价的，即

$f\left(t\right)\&LeftRightArrow;F\left(\mathrm{\&omega;}\right)---\left(8\right)$

对于离散信号，同样有时域和变换域表现的等价关系，即离散信号与其离散傅立叶变换（DFT）的等价关系：

$f\left(n\right)\&LeftRightArrow;F\left(k\right)---\left(9\right)$

公式（8）和式（9）可以这样理解：等价的双方可以互相转换而不改变或者丢失应有的信息，并且是互相对应的，即时域信息对应唯一的频域信息，由一方可以重建另一方。

由瞬时相位的定义可以证明φ(t)和信号x(n)是等价的，或者说φ(t)包含了原始信号x(n)的全部信息，下面通过由φ(t)重构x(n)的迭代算法对这个结论进行再次进行证明：

1）通过 $\tilde{a}\left(n\right)=x\left(n\right)+\mathrm{jH}\left[x\left(n\right)\right]=\left|a\left(n\right)\right|e^{\mathrm{j\&phi;}\left(n\right)}$ 计算φ(t)；

2）用一个正数c代替|a₀(n)|，令

3）计算

的N点DFT得到A₁(k)，由于解析信号的傅立叶变换是φ的单边函数，在φ<0处取值为0，所以有 $A_1^{\&prime;}\left(k\right)=\left\{\begin{array}{cc}{A}_1\left(k\right)& 0\&le;\frac{N}{2}-1\\ 0& \frac{N}{2}\&le;k\&le;N-1\end{array}\right.,$ 式中N都取2^M形式；

4）求A′₁(k)的傅立叶逆变换IDFT得到

其中|a′₁(n)|是a(n)的新估算值，p为迭代次数，每次迭代用

代替

只要解析信号z变换不满足

即z变换的零点中没有互为共轭倒数的零点，此时算法就是收敛的。为了确定迭代次数，引入相对平方误差概念，即：

$E_p=\frac{\underset{n=0}{\overset{N-1}{\mathrm{\&Sigma;}}}{[a_{p+1}\left(n\right)-a_p\left(n\right)]}^2}{\underset{n=0}{\overset{N-1}{\mathrm{\&Sigma;}}}{\left[a_p\left(n\right)\right]}^2}---\left(10\right)$

迭代次数越大，运算精度越高，但是迭代次数的增加会加大计算量，因此在实际运算过程中，为了提高运算效率，可以取相对平方误差小于10^-3。

基于以上广义希尔伯特变换的原理，本发明的基于广义瞬时相位及P范数负模的地震弱信号处理方法，包括以下内容：

1、对地震信号x(n)的弱信号进行识别，具体过程为：

1）对输入的地震信号x(n)进行P阶广义希尔伯特变换并根据公式（6）得到广义瞬时相位φ(t)；如图1所示，图1（A）为无噪信号，图1（B）为含噪信号，其中（A）为仿真信号，（B）为传统HT得到的瞬时相位，(c)为2阶GHT得到的瞬时相位，由图1（A）可以看到，在无噪情况下，两种变换提取的瞬时相位具有明显的分辨特征和相同的变化规律，图1（B）是含噪情况下（信噪比SNR为0.12）提取的瞬时相位，在此信噪比水平下，单频信号几乎完全淹没在随机噪声中，基本看不出单频特征；根据HT提取的瞬时相位不能反应信号原有特征，根据GHT得到的广义瞬时相位能很好地反应唔噪信号的瞬时相位变换趋势，具有较高的抗噪性，通常在不同信噪比水平下，GHT提取的广义瞬时相位抗噪性都比HT高，图1（B）中信号的信噪比已经远远小于实际地震资料的信噪比水平，因此GHT的抗噪性能在多数情况下满足实际地震弱信号的有效识别。

2）根据公式（7）求取广义瞬时相位P范数负模γ(t)；如图2所示，图2（A）显示了单频信号的广义瞬时相位及其P范数负模曲线，（A）为单频信号，（B）为广义瞬时相位，（c）为广义瞬时相位P范数负模。由于广义瞬时相位无法直观地看出瞬时相位和信号本身之间的关系，因此广义瞬时相位负模则较为直观地表现出和信号的关系，两个曲线峰值的相关性较高。图2（B）是含噪信号提取的广义瞬时相位负模（信噪比为0.02）曲线，即使噪音的加入干扰了信号的准确识别，广义瞬时相位P范数负模仍能较好地对信号峰值进行指示。

3）根据广义瞬时相位P范数负模γ(t)提取地震弱信号的有关属性参数，通过得到的有关属性参数可以对地震有效信号进行识别和应用，属性参数可以使地震弱信号的振幅、周期和峰值。如图3所示，其中（A）为无噪记录，（B）为含噪记录（SNR=2.53），（c）为广义瞬时相位，（d）为广义瞬时相位P范数负模。图3（B）是无噪合成记录和相位负模的对比，从中可以看到，相位负模对信号峰值的指示性较好。

2、如图4所示，对输入的地震信号x(n)的噪声进行去除，得到信噪比比较高的地震信号，具体过程为：

1）对输入的地震信号x(n)进行P阶广义希尔伯特变换并根据公式（6）得到广义瞬时相位φ(t)；

2）选择合适的滤波方法对计算得到的广义瞬时相位φ(t)进行平滑滤波处理，得到更加光滑连续的广义瞬时相位信号

滤波方法可以根据实际需要采用现有技术的各种滤波方法，在此不作限定；

3）采用逐步迭代的方法求取地震信号的瞬时振幅a(n)，具体过程为：

①将地震信号表示为广义瞬时相位信号

的形式，即：

②采用一假设的正数c代替①中振幅的初始值迭代值|a₀(n)|，令

③计算

的N点DFT得到A₁(k)，由于解析信号的傅立叶变换是

的单边函数，在

处取值为0，所以有 $A_1^{\&prime;}\left(k\right)=\left\{\begin{array}{cc}{A}_1\left(k\right)& 0\&le;\frac{N}{2}-1\\ 0& \frac{N}{2}\&le;k\&le;N-1\end{array}\right.,$ 式中N都取2^M形式；

④求A′₁(k)的傅立叶逆变换IDFT得到

其中|a′₁(n)|是a(n)的新估算值；

⑤判断|a′₁(n)|是否满足精度要求（即前后两次迭代得到的瞬时振幅差小于预先给定的数值），如果不满足精度要求，则返回步骤③继续迭代；假设迭代次数为p，则

每次迭代时采用

代替

迭代结果如果满足精度要求，（即前后两次迭代得到的瞬时振幅差小于预先给定的数值），则输出得到结果a(n)=|a'_p(n)|；

4)根据求取的广义瞬时相位信号

和瞬时振幅a(n)重构地震信号

5）将重构地震信号与现有测井资料得到的数据进行比对，判断重构地震信号的合理性，如果判断合理则输出重构信号

如果不合理，那么可以对广义瞬时相位提取参数n、瞬时相位曲线滤波参数、瞬时振幅迭代精度等进行修改，使效果达到最佳。

如图5所示，图5（A）是无噪信号的重构效果，从上到下分别为原始合成地震记录、重构的地震记录和重构的误差。从中可以看到，基于广义瞬时相位能较好的重构信号，重构误差较小。图5（B）为含噪信号的重构去噪效果（原始信号的SNR为4.72），从上往下分别为原始合成记录、加噪合成记录、重构的地震记录和重构的误差。重构记录的信噪比水平得到明显提高，且对原始信号的损伤较小，达到噪音滤除的效果。表1给出了图5（A）所示的模型信号在不同信噪比下重构的相对误差统计结果，从中可以看出，该种方法对原始信号的损伤随着信噪比的降低而增加。但是实际在地震资料处理中，常见的信噪比往往在15以上，即该重构方法从理论上来讲可以适用于普通地震资料信噪比范围的重构。

表1相对平方误差统计

含噪情况	无噪信号	SNR=31.46	SNR=16.41	SNR=9.12	SNR=4.72
						相对平方误差	0.627%	0.935%	1.023%	1.287%	1.373%

如图6所示，图6（A）是一个模型剖面CMP道集，分别有四个反射界面，其中界面一、三、四是正向的，界面二是反向的，除了界面一是水平的，其他界面都有一定角度。图6（B）是含噪CMP道集（SNR为0.653），图6(C)为去噪后的结果，图6(D)为滤除的噪音。可以看到，基于广义瞬时相位的去噪方法较好地去除了随机噪音，对有效信号的损伤较小。图7（A）是图7（B）分别为去噪前后资料的时频谱，高频噪音得到了有效的滤除，有效信号的时频特征得到较好的保持。图8（A）、图8（B）和图8(C)分别为实际资料去噪前、后及滤除噪音的地震剖面。去噪后的剖面同相轴的连续性增强，地震信号能量得到加强，整个剖面显得干净了很多，随机噪声的滤除效果较好。

上述各实施例仅用于说明本发明，其中方法的各实施步骤等都是可以有所变化的，凡是在本发明技术方案的基础上进行的等同变换和改进，均不应排除在本发明的保护范围之外。

Claims (5)

1.一种基于广义瞬时相位及P范数负模的地震弱信号处理方法，其包括以下内容：

1）对地震信号进行识别，具体过程为：

1.1）对输入的地震信号x(n)进行P阶广义希尔伯特变换并得到广义瞬时相位φ(t)；

1.2）求取广义瞬时相位P范数负模γ(t)；

1.3）根据广义瞬时相位P范数负模γ(t)提取地震弱信号的相关属性参数，通过提取的属性参数对地震有效信号进行识别；

2）对地震信号的噪声进行去除，得到信噪比比较高的地震信号，具体过程为：

2.1）对输入的地震信号x(n)进行P阶广义希尔伯特变换并得到广义瞬时相位φ(t)；

2.2）对计算得到的广义瞬时相位φ(t)进行平滑滤波处理，得到滤波后的广义瞬时相位信号

2.3）采用逐步迭代的方法求取地震信号的瞬时振幅a(n)；

2.4）根据求取的广义瞬时相位信号

和瞬时振幅a(n)重构地震信号

2.5）将重构地震信号与测井资料得到的数据进行比对，对重构地震信号的合理性进行判断，如果合理则输出重构信号

如果不合理，则返回步骤2.2）。

2.如权利要求1所述的基于广义瞬时相位及P范数负模的地震弱信号处理方法，其特征在于：广义瞬时相位

式中，hr(t)和hi(t)是复数道h(t)的实部和虚部。

3.如权利要求1所述的基于广义瞬时相位及P范数负模的地震弱信号处理方法，其特征在于：广义瞬时相位P范数负模γ(t)：

$\mathrm{\&gamma;}\left(t\right)=-\left|\arctan \frac{\mathrm{hi}\left(t\right)}{\mathrm{hr}\left(t\right)}\right|$

式中，hr(t)和hi(t)是复数道h(t)的实部和虚部。

4.如权利要求2所述的基于广义瞬时相位及P范数负模的地震弱信号处理方法，其特征在于：广义瞬时相位P范数负模γ(t)：

$\mathrm{\&gamma;}\left(t\right)=-\left|\arctan \frac{\mathrm{hi}\left(t\right)}{\mathrm{hr}\left(t\right)}\right|$

式中，hr(t)和hi(t)是复数道h(t)的实部和虚部。

5.如权利要求1或2或3或4所述的基于广义瞬时相位及P范数负模的地震弱信号处理方法，其特征在于：所述步骤2.3）采用逐步迭代的方法求取地震信号的瞬时振幅a(n)，包括以下步骤：

①将地震信号表示为广义瞬时相位信号

的形式，即：

②采用一假设的正数c代替①中振幅的初始值迭代值|a₀(n)|，令

③计算

的N点DFT得到A₁(k)，由于解析信号的傅立叶变换是

的单边函数，在处取值为0，所以有 $A_1^{\&prime;}\left(k\right)=\left\{\begin{array}{cc}{A}_1\left(k\right)& 0\&le;\frac{N}{2}-1\\ 0& \frac{N}{2}\&le;k\&le;N-1\end{array}\right.,$ 式中N都取2^M形式；

④求A′₁(k)的傅立叶逆变换IDFT得到

其中|a′₁(n)|是a(n)的新估算值；

⑤判断|a′₁(n)|是否满足精度要求，如果不满足精度要求，则返回步骤③继续迭代；假设迭代次数为p，则

迭代时采用

代替

迭代结果如果满足精度要求，则输出得到结果a(n)=|a'_p(n)|。

Images