一种基于模糊域的地震信号分数域最优阶计算方法

摘要

一种基于模糊域的地震信号分数域最优阶计算方法，属于地震信号处理领域，旨在利用模糊域快速确定地震信号分数域的最优阶，简化最优阶计算过程。首先读入单道地震信号，并对输入信号做模糊变换；然后计算信号的时宽和带宽，求取信号基于时频带宽积最小的窗函数，并对窗函数做模糊变换；利用窗函数的模糊函数对信号的模糊变换做低通滤波；计算出地震信号自项的模糊函数的水平和竖直两个方向的方向梯度；对方向梯度矩阵进行奇异值分解，计算地震信号自项模糊函数的主方向；最后通过主方向计算出地震信号的分数域最优阶。

说明书

技术领域

本发明属于地震信号处理领域，具体涉及一种基于模糊域的地震信号分数域最优阶计算方法。

背景技术

分数阶傅里叶变换(fractional Fourier transform,FrFT)也称为角度傅里叶变换，是傅里叶变换的广义形式。对信号做FrFT，也就是把信号转化到分数域。FrFT具有多个变换阶次，不同的变换阶次有不同的变换结果，信号的傅里叶变换只有唯一的一种变换结果。Almeida(1994)将FrFT引入信号处理领域，给出了FrFT的物理意义，其意义是信号的FrFT运算为信号在时频平面上的旋转，并初步提出了FrFT和传统时频分析方法结合得到分数域时频分析方法的思想。把FrFT引入到传统时频分析，由于FrFT的变换多样性，为提供高分辨率的时频分析结果提供了可能性。

一般来说，在众多的变换阶次中，有一个阶次对应的分数域时频分析结果的时频分辨率最高，这个阶次被称为分数域时频分析的最优阶。分数域时频分析方法的最优阶有多种计算方法。Capus和Brown(2003)提出了全局最优和局部最优的短时分数阶傅里叶变换，并通过计算信号的线性调频率来计算最优阶；Durak和Arikan(2003)利用分数域最小时频带宽积计算信号的最优阶；Akan等(2007)提出了利用时频分析结果的谱范数比值最小计算最优阶的方法。史军等(2010)运用分数域小波变换去除高斯白噪声，并把去噪均方误差最小的阶次作为最优阶。

由于分数域时频分析方法具有较高的分辨率，陈红等(2011)首先把分数阶Gabor变换运用到地震信号频谱成像技术，并利用广义时频带宽积最小的方法搜索分数阶Gabor变换的最优阶。随着分数域时频分析方法在地震信号领域的研究不断深入，陈颖频等(2011)将分数域时频分析方法运用于地震信号处理，利用分数域幅值最大化方法搜索地震信号的最优阶；田琳等利用分数域峰度系数最大化方法搜索地震信号的最优阶；王雨青等(2015)搜索地震信号分数域二阶距最小的阶次定义为最优阶。

对信号进行FrFT变换，本质上是对信号以线性调频信号进行分解。求解地震信号的最优阶，不同于其他信号，地震信号是多分量信号，而且不知道信号所包含的线性调频率，因此不能由信号本身直接计算线性调频率。因此上述的分数域时频分析方法都采用了遍历搜索的方法，不能直接计算地震信号的最优阶，目前有研究从数学角度出发将遍历搜索转化为计算过程，但是运算量大，不简洁。

发明内容

本发明提供了一种基于模糊域的地震信号分数域最优阶计算方法,旨在利用模糊域快速确定地震信号分数域的最优阶，简化最优阶计算过程。

为了解决上述技术问题，达到上述目的，本发明采用如下技术方案：

一种基于模糊域的地震信号分数域最优阶计算方法，所述方法包括以下步骤：

步骤1：输入一单道地震信号s(t)，对信号s(t)做模糊变换,得到信号的模糊函数A_s；

步骤2：计算地震信号s(t)的时宽T_s和带宽B_s，求取地震信号s(t)基于时频带宽积最小的窗函数h_TBP(t)，并对窗函数h_TBP(t)做模糊变换，得到窗函数h_TBP(t)的模糊函数A_g；

步骤3：利用模糊函数A_g对模糊函数A_s做低通滤波，得到信号自项的模糊函数A₀；

步骤4：计算信号自项的模糊函数A₀的方向梯度矩阵G；

步骤5：对方向梯度矩阵G进行奇异值分解，计算信号自项的模糊函数A₀的主方向θ；

步骤6：利用信号自项的模糊函数A₀的主方向θ计算信号s(t)的分数域最优阶p_opt。

进一步的，所述步骤1中对信号s(t)做的模糊变换所采用的变换公式如下：

式中，τ是时延，θ是频移，j是虚数单位，e是自然常数，t是时间。

进一步的，所述步骤2中

对窗函数h_TBP(t)做模糊变换，得到窗函数h_TBP(t)的模糊函数A_g，变换公式如下：

进一步的，所述步骤3中的滤波方法，计算方法如下：

A₀＝A_sA_g(4)。

进一步的，所述步骤4中信号自项的模糊函数A₀是N×N的矩阵信号,其方向梯度为：

其中g_x(k)表示矩阵A₀在点(x_k，y_k)水平方向的梯度，g_y(k)表示矩阵A₀在点(x_k，y_k)竖直方向的梯度。

进一步的，所述步骤5中对方向梯度矩阵G进行奇异值分解，具体理论如下：

G＝USV^T(6)

式中，U是N×2的矩阵，V是2×2的矩阵，S是2×2对角的奇异值矩阵，矩阵V的第一行向量为V₁＝[V_1,1,V_1,2]

信号自项的模糊函数A₀的主方向θ的计算公式为：

进一步的，所述步骤6中

综上所述，由于采用了上述技术方案，本发明的有益效果是：

将地震多分量信号转换到模糊域的时延-频偏平面上，地震信号的自项分布在二维平面的原点，交叉项分布于远离原点的位置，利用信号时频带宽积最小的高斯窗函数滤除地震信号的交叉项，比采用一般的窗函数滤波效果更好，可以更好的保留信号内在的有效特征，提高信号的时频分辨率。

该方法将信号转换到模糊域，在模糊域可直接计算出地震信号的最优阶，避免了现有技术的遍历搜索的弊端，简化计算过程，可以有效地运用于地震信号分数域时频分析。

附图说明

图1为方法流程图；

图2为单道地震信号；

图3为信号的模糊函数；

图4为窗函数的模糊函数；

图5信号自项的模糊函数。

具体实施方式

本说明书中公开的所有特征，除了互相排斥的特征和/或步骤以外，均可以以任何方式组合。

下面结合附图对本发明作详细说明。

步骤1：输入一单道地震信号s(t)(如图2所示)；

对信号s(t)做模糊变换,得到信号的模糊函数A_s(如图3所示)，所采用的变换公式如下：

式中，τ是时延，θ是频移，j是虚数单位，e是自然常数

步骤2：计算信号s(t)的时宽T_s和带宽B_s，求取信号s(t)基于时频带宽积最小的窗函数h_TBP(t)

对窗函数h_TBP(t)做模糊变换，得到窗函数h_TBP(t)的模糊函数A_g(如图4所示)

步骤3：利用模糊函数A_g对模糊函数A_s做低通滤波，得到信号自项的模糊函数A₀(如图5所示)

A₀＝A_sA_g(11)

步骤4：计算信号自项的模糊函数A₀的方向梯度矩阵G，信号自项的模糊函数A₀是N×N的矩阵信号,其方向梯度为：

其中g_x(k)表示矩阵A₀在点(x_k，y_k)水平方向的梯度，g_y(k)表示矩阵A₀在点(x_k，y_k)竖直方向的梯度。

步骤5：对方向梯度矩阵G进行奇异值分解，分解公式如下：

G＝USV^T(13)

式中，U是N×2的矩阵，V是2×2的矩阵，S是2×2对角的奇异值矩阵，矩阵V的第一行向量为V₁＝[V_1,1,V_1,2]

信号自项的模糊函数A₀的主方向θ的计算公式为：

步骤6：利用信号自项的模糊函数A₀的主方向θ计算信号s(t)的分数域最优阶p_opt

如上所述即为本发明的实施例。本发明不局限于上述实施方式，任何人应该得知在本发明的启示下做出的结构变化，凡是与本发明具有相同或相近的技术方案，均落入本发明的保护范围之内。