[发明专利]基于广义Morse标架的地震瞬时属性提取方法

摘要

本发明公开了一种基于广义Morse标架的地震瞬时属性提取方法，采用广义Morse标架，将有效信号能量分布空间的确定转化成一个优化问题，利用迭代萎缩阈值算法和快速迭代萎缩阈值算法求解。在确定了有效信号能量分布空间的基础上，利用小波变换与Hilbert变换的关系，提出了含噪信号瞬时属性分析的方法。本发明的方法具有良好的抗噪性能和准确性，得到的瞬时频率剖面能够更加清楚地反映瞬时频率的变化，指示异常区域。

主权项

基于广义Morse标架的地震瞬时属性提取方法，其特征在于包括以下步骤：步骤1：采集地震数据体；步骤2：计算各道含噪地震数据对应的广义Morse标架系数；实地震道s(t)的小波变换定义为：$W_{\mathrm{\&psi;}}(t,a)\&equiv;\frac{1}{a}{\&Integral;}_{-\mathrm{\&infin;}}^{+\mathrm{\&infin;}}s\left(\mathrm{\&tau;}\right){\mathrm{\&psi;}}^{*}\left(\frac{\mathrm{\&tau;}-t}{a}\right)d\mathrm{\&tau;},---\left(1\right)$式中基本小波ψ(t)选用广义Morse小波；广义Morse小波频域表达式为：${\mathrm{\&Psi;}}_{\mathrm{\&beta;},\mathrm{\&gamma;}}\left(\mathrm{\&omega;}\right)=U\left(\mathrm{\&omega;}\right){\mathrm{\&alpha;}}_{\mathrm{\&beta;},\mathrm{\&gamma;}}{\mathrm{\&omega;}}^{\mathrm{\&beta;}}{e}^{-{\mathrm{\&omega;}}^{\mathrm{\&gamma;}}},---\left(2\right)$其中U(ω)为单位阶跃函数，β和γ为小波的参数，且β＞0，γ＞0，α_β,γ为归一化常量，α_β,γ≡2(eγ/β)^β/γ；离散化后的小波族可以表示为：${\mathrm{\&psi;}}_{m,n}\left(x\right)=a_0^{-m}\mathrm{\&psi;}\&lsqb;a_0^{-m}(x-{\mathrm{nt}}_0)\&rsqb;---\left(3\right)$其中a₀尺度因子a的离散化步长且a₀＞1，t₀为平移因子t的离散化步长；和公式(1)相对应的离散化后的小波变换表示为：C_m,n＝<s,ψ_m,n>, (4)其中C_m,n为标架系数；步骤3：对于各道的广义Morse标架系数迭代得到有效信号对应的系数；定义K为从映射到的算子，将标架系数映射为上的信号，即算子K的伴随算子K*为从映射到的算子，将上的信号投影到标架系数上：K^*s＝<s,ψ′_m,n>, (6)则s＝KK^*s, (7)算子K为合成算子，K^*为分析算子和小波标架{ψ′_m,n}相对应；将含噪信号表示为y＝s+n＝Kx+n, (8)其中y表示含噪信号，s为不含噪的有效信号，n为高斯白噪声，K为公式(7)中的合成算子，x表示变换域的系数，通过求解下列优化问题得到有效信号对应的系数$\begin{array}{ccc}\tilde{x}=\arg \underset{x}{min}\left|\right|x|{|}_1& s.t.& \left|\right|y-Kx|{|}_2\&le;\mathrm{\&epsiv;}\end{array},---\left(9\right)$其中ε和待分析信号的噪声水平有关，用Lagrange乘子λ将该问题转化为以下的无约束问题：${\tilde{x}}_{\mathrm{\&lambda;}}=\arg \underset{x}{min}\frac{1}{2}\left|\right|y-Kx|{|}_2^2+\mathrm{\&lambda;}|\left|x\right|{|}_1,---\left(10\right)$其中λ被称为正则化参数，采用下列公式所示的迭代方法对系数x进行更新x^(k+1)＝T_λ[x^(k)+K^*(y‑Kx^(k))],k＝1,2,…,N, (11)其中T_λ为阈值函数，定义为$T_{\mathrm{\&lambda;}}\left(x\right)=\left\{\begin{array}{cc}x-\mathrm{\&lambda;}\frac{x}{\left|x\right|}& \left|x\right|\&GreaterEqual;\mathrm{\&lambda;},\\ 0& \left|x\right|\&lt;\mathrm{\&lambda;}\end{array}\right.---\left(12\right);$步骤4：由公式(13)计算解析信号；$c\left(t\right)\&equiv;\frac{1}{K_{\mathrm{\&psi;}}}{\&Integral;}_0^{\mathrm{\&infin;}}{W}_{\mathrm{\&psi;}}(t,a)\frac{da}{a}=s\left(t\right)+ih\left(t\right),---\left(13\right)$其中W_ψ(t,a)是有效信号对应的系数，h(t)是s(t)的Hilbert变换，c(t)是s(t)对应的解析信号；步骤5：利用c(t)计算瞬时振幅、瞬时相位和瞬时频率：$\hat{e}\left(t\right)=\sqrt{{\mathrm{Re}}^2\&lsqb;c\left(t\right)\&rsqb;+{\mathrm{Im}}^2\&lsqb;c\left(t\right)\&rsqb;},---\left(14\right)$$\widehat{\mathrm{\&theta;}}\left(t\right)=\arctan \&lsqb;\frac{\mathrm{Im}\&lsqb;c\left(t\right)\&rsqb;}{\mathrm{Re}\&lsqb;c\left(t\right)\&rsqb;}\&rsqb;,---\left(15\right)$$\hat{f}\left(t\right)=\frac{1}{2\mathrm{\&pi;}}\frac{d}{dt}\&lsqb;\arctan \&lsqb;\frac{\mathrm{Im}\&lsqb;c\left(t\right)\&rsqb;}{\mathrm{Re}\&lsqb;c\left(t\right)\&rsqb;}\&rsqb;\&rsqb;,---\left(16\right)$其中Re[c(t)]和Im[c(t)]分别表示c(t)的实部和虚部。