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

摘要

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

权利要求

1.一种基于模糊域的地震信号分数域最优阶计算方法，其特征在于，所述方法包括以下步骤：步骤1：输入一单道地震信号s(t)，对信号s(t)做模糊变换，得到信号的模糊函数As；

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

步骤3：利用窗函数的模糊函数Ag对信号的模糊函数As做低通滤波，得到信号自项的模糊函数A0；

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

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

具体理论如下：

T

G＝USV  (1)

式中，U是N×2的矩阵，V是2×2的矩阵，S是2×2对角的奇异值矩阵，矩阵V的第一行向量为V1＝[V1，1，V1，2]信号自项的模糊函数A0的主方向θ的计算公式为：步骤6：利用信号自项的模糊函数A0的主方向θ计算信号s(t)的分数域最优阶popt，

2.根据权利要求1所述的一种基于模糊域的地震信号分数域最优阶计算方法，其特征在于，所述步骤1中对信号s(t)做的模糊变换所采用的变换公式如下：式中，τ是时延，θ是频移，j是虚数单位，e是自然常数，t是时间。

3.根据权利要求1所述的一种基于模糊域的地震信号分数域最优阶计算方法，其特征在于，所述步骤2中对窗函数hTBP(t)做模糊变换，得到窗函数hTBP(t)的模糊函数Ag，变换公式如下：

4.根据权利要求1所述的一种基于模糊域的地震信号分数域最优阶计算方法，其特征在于，所述步骤3中的滤波方法，计算方法如下：Ao＝AsAg    (7)。

5.根据权利要求1所述的一种基于模糊域的地震信号分数域最优阶计算方法，其特征在于，所述步骤4中信号自项的模糊函数A0是N×N的矩阵信号，其方向梯度为：其中gx(k)表示矩阵A0在点(xk，yk)水平方向的梯度，gy(k)表示矩阵A0在点(xk，yk)竖直方向的梯度。

说明书

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

技术领域

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

背景技术

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

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

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

[0005] 对信号进行FrFT变换，本质上是对信号以线性调频信号进行分解。地震信号是多分量信号，信号所包含的线性调频信号的调频率是未知的，因此不能由信号本身的线性调频率直接计算分数域最优阶，因此上述的分数域时频分析方法都采用了遍历搜索的方法，运算量较大。对地震信号做FrFT变换，总能找到某一个线性调频信号是占优的，这个占优的线性调频信号的调频率即为地震信号的最优阶。单个的线性调频信号的模糊函数在模糊平面上为一过原点的直线，且直线的斜率就是线性调频信号的调频率。本发明把地震信号转化到模糊域，利用地震信号模糊函数的主方向直接计算地震信号的最优阶。

发明内容

[0006] 本发明提供了一种基于模糊域的地震信号分数域最优阶计算方法,旨在把地震信号转换到模糊域，可以直接计算地震信号最优阶。

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

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

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

[0010] 步骤2：计算地震信号s(t)的时宽Ts和带宽Bs，求取地震信号s(t)基于时频带宽积最小的窗函数hTBP(t)，并对窗函数hTBP(t)做模糊变换，得到窗函数hTBP(t)的模糊函数Ag；

[0011] 步骤3：利用窗函数的模糊函数Ag对信号的模糊函数As做低通滤波，得到信号自项的模糊函数A0；

[0012] 步骤4：计算信号自项的模糊函数A0的方向梯度矩阵G；

[0013] 步骤5：对方向梯度矩阵G进行奇异值分解，利用奇异值分解结果计算信号自项的模糊函数A0的主方向θ；

[0014] 步骤6：利用信号自项的模糊函数A0的主方向θ计算信号s(t)的分数域最优阶popt。

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

[0016]

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

[0018] 进一步的，所述步骤2中

[0019]

[0020] 对窗函数hTBP(t)做模糊变换，得到窗函数hTBP(t)的模糊函数Ag，变换公式如下：

[0021]

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

[0023] A0＝AsAg   (4)。

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

[0025]

[0026] 其中gx(k)表示矩阵A0在点(xk，yk)水平方向的梯度，gy(k)表示矩阵A0在点(xk，yk)竖直方向的梯度。

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

[0028] G＝USVT   (6)

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

[0030] 信号自项的模糊函数A0的主方向θ的计算公式为：

[0031]

[0032] 进一步的，所述步骤6中

[0033]

[0034] 综上所述，由于采用了上述技术方案，本实用新型的有益效果是：

[0035] 利用信号时频带宽积最小的高斯窗函数滤除地震信号的交叉项，比采用一般的窗函数滤波效果更好，可以更好的保留信号内在的有效特征，提高信号的时频分辨率。

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

附图说明

[0037] 图1为方法流程图；

[0038] 图2为单道地震信号；

[0039] 图3为信号的模糊函数；

[0040] 图4为窗函数的模糊函数；

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

具体实施方式

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

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

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

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

[0046]

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

[0048] 步骤2：计算信号s(t)的时宽Ts和带宽Bs，求取信号s(t)基于时频带宽积最小的窗函数hTBP(t)

[0049]

[0050] 对窗函数hTBP(t)做模糊变换，得到窗函数hTBP(t)的模糊函数Ag(如图4所示)[0051]

[0052] 步骤3：利用模糊函数Ag对模糊函数As做低通滤波，得到信号自项的模糊函数A0(如图5所示)

[0053] A0＝AsAg   (11)

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

[0055]

[0056] 其中gx(k)表示矩阵A0在点(xk，yk)水平方向的梯度，gy(k)表示矩阵A0在点(xk，yk)竖直方向的梯度。

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

[0058] G＝USVT   (13)

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

[0060] 信号自项的模糊函数A0的主方向θ的计算公式为：

[0061]

[0062] 步骤6：利用信号自项的模糊函数A0的主方向θ计算信号s(t)的分数域最优阶popt[0063]

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