CN102221708B - 基于分数阶傅里叶变换的随机噪声压制方法 - Google Patents

Abstract

提供一种基于分数阶傅里叶变换的随机噪声压制方法，所述方法包括：针对采集的地震数据，利用正交多项式逼近法，采用最大相似系数准则，计算描述地震数据中同相轴走势的相关系数；利用计算的相关系数对地震数据在空间上进行拟线性化变换；根据分数阶傅里叶变换的最优阶次求取法，求取地震数据中信号最集中的分数阶域最优阶次；利用最优阶次下的分数阶傅里叶变换将地震数据从时空域变换到分数阶域；求解地震道的预测算子，利用所述预测算子预测出地震道的地震数据，由此消除无规律分布的随机噪声；利用分数阶傅里叶逆变换将去噪后的地震数据从分数阶域变换到时空域；利用计算的相关系数对地震数据进行反拟线性变换，完成地震数据的随机噪声压制。

Description

基于分数阶傅里叶变换的随机噪声压制方法

技术领域

本发明涉及石油地震勘探，属于地震勘探资料处理领域。更具体地讲，本发明涉及一种对低信噪比资料中的随机噪声进行压制来提高地震资料的信噪比的方法，该方法能够使得去噪效果与信号保真达到最佳的折中。

背景技术

利用地下介质弹性和密度的差异，通过观测和分析大地对人工激发地震波的响应，推断地下岩层的性质和形态的地球物理勘探方法称为地震勘探。地震勘探是钻探前勘测石油与天然气资源的重要手段，在煤田和工程地质勘查、区域地质研究和地壳研究等方面，也得到广泛应用。在地震勘探中，最为基础和重要的工作是如何采集到野外地震数据并对其进行有效的处理和解释。为了能够获得有效的地震数据，需要在野外设置大量的监测仪器，其中，每一个监测仪器采集到的数据称作一道地震数据。

复杂山地地震数据噪声来源十分复杂，噪声对有效波有较强的干扰，尤其是无处不在的随机干扰，导致地震资料的信噪比低，严重影响了地震数据的分析和处理。而信噪比是基础，提高信噪比是地震数据处理的首要任务，所以必须对地震数据的随机干扰进行衰减和压制。

目前，现有的与频率域相关的去噪方法都是基于经典傅里叶变换或在其基础上延伸出来的去噪方法，如频率波数(F-K)域滤波、频率空间(F-X)域去噪等。然而，经典傅里叶变换存在明显的缺陷，它是一种全局性变换，得到的是地震信号的整体频谱，因而无法表述信号的时频局部特性，而这种特性正是非平稳信号最根本和最关键的性质。有效的去噪方法是寻求一种在时域(空域)或频域上的信噪自适应分离。一般而言，地震信号的噪声覆盖了整个频谱段，因此，不论在时域或频域都难以将二者区分开来。实际数据处理中，通常是采用了折中的办法，即，在去掉噪声的同时，也会牺牲部分有效信号。

发明内容

本发明的目的在于提供一种能够针对复杂山地低信噪比数据中的随机噪声进行基于分数阶下拟线性空间预测去噪方法。

根据本发明的一方面，提供一种基于分数阶傅里叶变换的随机噪声压制方法，所述方法包括：(1)针对采集的地震数据，利用正交多项式逼近法，采用最大相似系数准则，计算描述地震数据中同相轴走势的相关系数；(2)利用计算的相关系数对地震数据在空间上进行拟线性化变换；(3)根据分数阶傅里叶变换的最优阶次求取法，求取地震数据中信号最集中的分数阶域最优阶次；(4)利用最优阶次下的分数阶傅里叶变换将地震数据从时空域变换到分数阶域；(5)求解地震道的预测算子，利用所述预测算子预测出地震道的地震数据，由此来消除无规律分布的随机噪声；(6)利用分数阶傅里叶逆变换将去噪后的地震数据从分数阶域变换到时空域；(7)利用计算的相关系数对地震数据进行反拟线性变换，完成地震数据的随机噪声压制，得到消除随机噪声后的有效信号的地震剖面。

在步骤(1)中，可利用正交多项式逼近法对表示叠加地震数据的反射同相轴的时间多项式进行变换，逐项扫描变换后的时间多项式的系数，将扫描得到的系数作为所述相关系数。

在逐项扫描变换后的时间多项式的系数的过程中，每次扫描新的系数时使先前确定的系数保持不变。

可在叠加地震数据的剖面上设置多个时窗，在每个时窗内通过多地震道互相关扫描来确定该时窗内的时间多项式系数。

在步骤(3)中，可根据下面的等式来计算所述最优阶次：

$\tan \left(2\mathrm{p\π}\right)=\frac{{4\mathrm{\μ}}_0(w_0-w_1)}{4{\mathrm{\μ}}_0^2-{(w_0-w_1)}^2}$

其中，p表示分数阶，μ₀＝(w₀+w₁)/2-ww_0.5，w₀是分数阶p＝0阶时地震信号在分数阶域的二阶矩，w₁是分数阶p＝1阶时地震信号在分数阶域的二阶矩，w_0.5是分数阶p＝0.5阶时地震信号在分数阶域的二阶矩。

步骤(4)可包括：从连续傅里叶变换的特征函数出发，通过对特征函数进行离散化近似和正交投影，得到一组与特征函数形状相似的离散傅里叶变换矩阵的正交化离散特征函数的特征向量；根据连续分数阶傅里叶变换的核函数谱分解表达式，构造离散分数阶傅里叶变换核矩阵；将地震数据与最优阶次下的分数傅里叶变换核矩阵相乘，得到变换后的分数阶域下的地震数据。

在步骤(5)中，可根据地震数据中相干信号可预测，随机噪声不可预测，并且同相轴为线性，来求解所述预测算子。

可对预测算子的求取进行线性变换，通过复数维纳滤波求得预测算子。

在利用所述预测算子预测出地震道的地震数据的步骤中，在特定频率上，对于一个预测算子长度内的地震道，可利用所述预测算子长度内的所有地震道的地震数据分别与对应的预测算子相乘，，再将所述预测算子长度内的相乘结果相加作为下一地震道的地震数据的预测结果。

附图说明

通过结合附图，从下面的实施例的描述中，本发明这些和/或其它方面及优点将会变得清楚，并且更易于理解，其中：

图1是根据本发明示例性实施例的基于分数阶傅里叶变换的随机噪声压制方法的流程图；

图2A至图2F是根据本发明的理论模型的测试效果图；

图3A至图3F分别是图2A至图2F所示的理论模型数据中第10地震道的频谱；

图4A至图4C示出了根据本发明的地质情况更为复杂的实际地震数据的应用效果图。

具体实施方式

以下，参照附图来详细描述本发明的实施例。在附图中，相同的标号表示相同的部件。

本发明将分数阶傅里叶变换(fractional Fourier transform，FRFT)的方法引入到拟线性的F-X域预测滤波技术中，把二维地震数据变换到任一分数阶空间域中，进行预测处理，可实现预测滤波效果与信号特征保持的最佳折中。

图1是根据本发明示例性实施例的基于分数阶傅里叶变换的随机噪声压制方法的流程图。

参照图1，在步骤101，针对采集的地震数据，利用正交多项式逼近法，采用最大相似系数准则，计算描述地震数据中同相轴走势的相关系数。

在进行预测去噪之前，需要计算出描述地震数据中同相轴走势的相关系数，以便于对地震数据在空间上做拟线性化变换。计算相关系数的具体步骤如下：

叠加地震数据的反射同相轴可用时间多项式来描述：

t(x)＝c₀+c₁x+c₂x²+c₃x³+…(1)

如果已知等式(1)中的系数就可容易将其进行线性化变换。对地震数据做空间方向拟线性化变换，可在地震数据上进行扫描，通过使相似系数最大来得到多项式系数。这里，相似系数定义如下：

$R(c_0,c_1,...)=\frac{\underset{i\≠j}{\mathrm{\Σ}}\underset{l}{\mathrm{\Σ}}s(t+l,x_i)s(t+l,x_j)}{\underset{i\≠j}{\mathrm{\Σ}}\sqrt{\underset{l}{\mathrm{\Σ}}{s}^2(t+l,x_i)\underset{l}{\mathrm{\Σ}}{s}^2(t+l,x_j)}}---\left(2\right)$

在等式(2)中，s(t，x)为地震道数据，i，j＝1、2、...、N，N为地震道总道数，l为用于计算相关的时窗样点数。

在等式(1)中各个系数是相关联的，如果简单地通过扫描来获得各个系数，则计算量相当大。为此，利用正交多项式逼近法对表示叠加地震数据的反射同相轴的时间多项式进行变换。具体地，这里采用Gram-Schmidt正交多项式逼近法，令

g₀(x)＝1                                    (3)

$g_i\left(x\right)=x^i+\underset{l=0}{\overset{i-1}{\mathrm{\Σ}}}{k}_{\mathrm{il}}{g}_l\left(x\right)---\left(4\right)$

$[g_i\left(x\right),g_j\left(x\right)]=\underset{n=0}{\overset{N}{\mathrm{\Σ}}}{g}_i\left(x_n\right)g_j\left(x_n\right)---\left(5\right)$

在等式(5)中，当i≠j时，g_i(x_n)g_j(x_n)＝0；当i＝j时，g_i(x_n)g_j(x_n)＞0；n为地震道序号，N为地震道总道数，x_n为窗口内第n地震道的相对横坐标。

等式(4)中的k_ij由下式求得：

$k_{\mathrm{il}}=-\frac{\mathrm{\Σ}{x}^i{g}_l\left(x\right)}{\mathrm{\Σ}{g}_l^2\left(x\right)}---\left(6\right)$

因此，时间多项式变为：

t(x)＝a₀g₀(x)+a₁g₁(x)+a₂g₂(x)+a₃g₃(x)+… (7)

逐项扫描a₀，a₁，....，，在每次扫描新的系数时使先前确定的系数保持不变，这样就完成了多项式系数的快速求取过程。将求取的多项式系数作为相关系数。多项式的次数由构造复杂程度决定，即，构造越复杂，次数越高。

同相轴相位时间的拟合是在叠加地震数据的剖面上设置多个小的时窗，在每个时窗中进行同相轴时间拟合。设时窗内在x方向上有N＝2M+1道(N为横向半时窗地震道数)，在时间t方向上有2L+1个样点(L为纵向半时窗样点数)，窗内各个地震道中点的时间不一样，通常可以用N-1次多项式表示。由于可以认为N-1次多项式的前几项反映信号，后几项反映噪声，所以一般取三次多项式进行拟合，即，等式(1)，而等式(1)式中的相关系数可通过等式(2)求取，但是求取比较麻烦，效率不高。通过等式(3)至等式(6)的变换，则等式(1)可由等式(7)表示，因此就只需求取等式(7)中一系列系数a₀、a₁、....。对于一系列正交多项式g(x)，可由等式(3)至等式(6)求取。

在每个时窗内使用上述等式(3)至等式(6)的多项式，通过多地震道互相关扫描来确定该时窗内的时间多项式系数a₀、a₁、....。时窗滑动特定的步长，重复下去，就可以得到每个时窗内的多项式系数a₀、a₁、....。然后，可根据各个地震道的空间坐标拟合出各个地震道的各个时窗内信号到达的时间(即，得到时间多项式t(x))。

因此，根据时间多项式系数a₀、a₁、....，可获得时间多项式t(x)，从而能够拟合地震数据。也就是说，先获得每个时窗内的多项式系数a₀、a₁、....，根据空间上各个地震道的空间坐标x，能够得到各个地震道的时间多项式t(x)。

返回参照图1，在步骤102，利用计算的相关系数对地震数据在空间上进行拟线性化变换。

通过等式(1)至等式(7)的操作，采用了正交多项式逼近法，对输入的整个地震数据进行拟线性变换，使得地震数据中同相轴的线性特征更明显(即，是拟线性的)。

然后，在步骤103，根据分数阶傅里叶变换的最优阶次求取法，求取地震数据中信号最集中的分数阶域最优阶次。

在求得地震数据的相关系数并进行拟线性变换后，在做分数阶傅里叶变换之前，需要求取地震数据中信号最集中的分数阶域最优阶次，具体步骤如下：

信号在时域或频域的宽度是由其二阶中心矩确定，与此相似，信号在分数阶域的宽度应与其二阶中心分数阶傅氏变换矩相关。

分数阶域S_p(u)的二阶矩W_p定义为：

$w_p={\left[{\∫}_{-\∞}^{\∞}{(u-{\mathrm{\η}}_p)}^2{\left|S_p\left(u\right)\right|}^2\mathrm{du}\right]}^{1/2}---\left(8\right)$

在等式(8)中：

μ为时间域中时间t变换到分数阶域中的变量。

由于在进行分数阶域变换计算前进行了尺度归一化，所以有η_p＝0。如果分数阶域的信号能量也进行了归一化，即，S_p(u)＝S_p(u)/|S_p(u)|，则分数阶域S_p(u)的二阶矩就是分数阶域的宽度W_p，即：

$w_p={\left[{\∫}_{-\∞}^{\∞}{u}^2{\left|S_p\left(u\right)\right|}^2\mathrm{du}\right]}^{1/2}---\left(9\right)$

应用分数阶傅里叶变换的旋转特性，可得到：

w_p＝w₀cos²(pπ/2)+w₁sin²(pπ/2)-μ₀sin(pπ)  (10)

其中，p表示分数阶，μ₀是中间运算算子，μ₀＝(w₀+w₁)/2-w_0.5，w₀是分数阶p＝0阶时地震信号在分数阶域的二阶矩，w₁是分数阶p＝1阶时地震信号在分数阶域的二阶矩，w_0.5是分数阶p＝0.5阶时地震信号在分数阶域的二阶矩。由此可知，只需要计算出p＝0、p＝0.5和p＝1三个分数阶域的信号宽度就可计算各个分数阶域的信号宽度，极大地减小了计算量。由基于分数傅里叶变换的最优高斯窗函数：g_GTP(t)＝exp(-πt²B_fp/T_fp)，可变为：

$\mathrm{GTBP}\left\{x\left(t\right)\right\}=\underset{0\mathrm{\πp\π}2}{\min }\mathrm{TBP}\left\{x\left(t\right)\right\}=\underset{0\mathrm{\πp\π}2}{\min }\{w_p\·w_{p+1}\}---\left(11\right)$

其中：

$w_p\·w_{p+1}=w_0{w}_1+\frac{1}{4}[{(w_0-w_1)}^2-4{\mathrm{\μ}}_0^2]{\sin }^2\left(\mathrm{p\π}\right)+\frac{1}{2}{\mathrm{\μ}}_0(w_0-w_1)\sin \left(2\mathrm{p\π}\right)---\left(12\right)$

对等式(12)求解关于分数阶次p的导数，并令其等于零，可得：

$\tan \left(2\mathrm{p\π}\right)=\frac{{4\mathrm{\μ}}_0(w_0-w_1)}{4{\mathrm{\μ}}_0^2-{(w_0-w_1)}^2}---\left(13\right)$

求解等式(13)即可得到适合信号最集中的分数阶域最优阶次p。这种求分数阶域最优阶次的方式过程简单，计算量小。

接着，在步骤104，利用最优阶次下的分数阶傅里叶变换将地震数据从时空域变换到分数阶域中。

将地震数据从时空域变换到分数阶域中的处理涉及到分数阶傅里叶变换的特征分解型离散算法，具体步骤如下：

从连续傅里叶变换的特征函数(Hermite函数)出发，通过对特征函数进行离散化近似和正交投影，得到一组与特征函数形状相似的离散傅里叶变换矩阵的正交化离散特征函数的特征向量；然后，根据连续分数阶傅里叶变换的核函数谱分解表达式，构造离散分数阶傅里叶变换核矩阵；最后，将地震数据与最优阶次下的分数傅里叶变换核矩阵相乘，得到变换后的分数阶域下的地震数据。

接着，在步骤105，求解地震道的预测算子，利用所述预测算子预测出地震道的地震数据，由此来消除无规律分布的随机噪声。

具体地，通过等式(1)至等式(7)的操作，采用了正交多项式逼近法，对输入的整个地震数据进行拟线性变换，使得地震数据中同相轴的线性特征更明显。因此，可根据地震数据中相干信号可预测，随机噪声不可预测，并且同相轴为线性，来求解预测算子。

本发明在常规的F-X域滤波离散算法中进行了改进。常规F-X域预测滤波去噪技术假定在F-X域相干信号是可预测的，而随机噪声是不可预测的，同时假定同相轴是线性的。这些假定条件在处理复杂地区地震资料时难以满足。为此提出一种F-X域适合拟线性变换的离散方法，具体步骤如下：

设地震数据S(t，x)由L组一共N个地震道组成(L＜N)，各个地震道的倾角可以不同，S(t，x)在时空域可表示为：

$S(t,x)=\stackrel{L}{\mathrm{\Σ}}{W}_k\left(t\right)\mathrm{\σ}(t-b_{k}x)---\left(14\right)$

L组地震道对应于L个时窗，即，每个时窗覆盖一组地震道。

对式(14)做时间方向的傅氏变换，有

$S(f,x)=\stackrel{L}{\mathrm{\Σ}}{W}_k\left(f\right)e^{i2\mathrm{\π}{\mathrm{fb}}_{k}x}---\left(15\right)$

等式(14)和等式(15)中，为地震子波(K为线性同相轴序号)，

为地震子波的傅氏变换，b_k为同相轴的斜率。

如果对式(15)在x方向进行离散，并且用b_k表示第K个同相轴的同一斜率间相邻地震道之间的时间差，那么数据可表示为矩阵形式：

$\left[\begin{array}{c}{S}_1\left(f\right)\\ S_2\left(f\right)\\ ...\\ S_N\left(f\right)\end{array}\right]=\left[\begin{array}{cccc}1& 1& ...& 1\\ Z_1& Z_2& ...& Z_L\\ ...& ...& ...& ...\\ Z_1^{N-1}& Z_2^{N-1}& ...& Z_L^{N-1}\end{array}\right]\left[\begin{array}{c}{W}_1\left(f\right)\\ W_2\left(f\right)\\ ...\\ W_L\left(f\right)\end{array}\right]---\left(16\right)$

或S(f)＝Z(f)W(f)                                    (17)

等式(16)中，S_i(f)为数据S_i(f，x)在特定频率的频谱；

等式(17)等号右边第一个矩阵中的每一行可用前面的行的线性组合求得，即

$(Z_1^{K-1},Z_2^{K-1},Z_3^{K-1},\·\·\·,Z_L^{K-1})=[P_L\left(f\right),P_{L-1}\left(f\right),\·\·\·P_1\left(f\right)]\×\left[\begin{array}{ccc}{Z}_1^{K-L-1}& ...& Z_L^{K-L-1}\\ ...& ...& ...\\ Z_1^{K-2}& ...& Z_L^{K-2}\end{array}\right]---\left(18\right)$

这里的P₁(f)，P₂(f)，...，P_L(f)与行号无关。如果将式(18)两边同乘以W(f)可得：

$S_i\left(f\right)=\left[\begin{array}{cccc}{P}_L\left(f\right),& P_{L-1}\left(f\right),& \·\·\·,& P_1\left(f\right)\end{array}\right]\×\left[\begin{array}{c}{S}_{K-L}\left(f\right)\\ S_{K-L-1}\left(f\right)\\ ...\\ S_{K-1}\left(f\right)\end{array}\right]$

或写成 $S_K\left(f\right)=\underset{j=1}{\overset{L}{\mathrm{\Σ}}}{P}_j\left(f\right)S_{K-j}\left(f\right)---\left(19\right)$

其中，S_K(f)表示第K地震道在频率f处的地震数据，K＝L+1，.....，N，预测算子p_j(f)可用复数维纳滤波求得，等式(19)的含义是S_K(f)在x方向可通过向前(或向后)一步预测滤波求得。

通过等式(14)至等式(18)的操作，在求取预测算子的时候，对预测算子的求取进行线性变换，通过复数维纳滤波求得预测算子，这样的求取方式对复杂地区较差的地震资料的预测效果会更好。

返回参照图1，在步骤106，利用分数阶傅里叶逆变换将去噪后的地震数据从分数阶域变换到时空域。

在步骤107，利用计算的相关系数对地震数据进行反拟线性变换，完成地震数据的随机噪声压制，得到消除随机噪声后的有效信号的地震剖面。

下面描述将本发明的基于分数阶傅里叶变换的随机噪声压制方法应用于理论数据的示例。

图2A至图2F是根据本发明的理论模型的测试效果图，其中，图2A是具有三个不同倾角地层的理论模型。在该理论模型中，没有任何噪声干扰，存在三个连续追踪的同相轴，包括一个平层、一个斜层和一个半椭圆层。

图2B是在图2A中加入随机噪声后的理论模型。可以看出，加入随机噪声后的理论模型的信噪比较低，有效信号被淹没在随机噪声中。

将采集到的地震数据运用正交多项式逼近法，并利用最大相似系数准则，扫描出描述地震数据中同相轴走势的相关系数(a₀，a₁，....)，然后利用这些相关系数对地震数据进行拟线性化变换，使得地震数据的同相轴是拟线性的。

然后，利用地震数据在分数阶域的二阶中心矩，计算地震数据中信号最集中的最优阶次。根据前面的等式(13)来计算所述最优阶次p：

$\tan \left(2\mathrm{p\π}\right)=\frac{{4\mathrm{\μ}}_0(w_0-w_1)}{4{\mathrm{\μ}}_0^2-{(w_0-w_1)}^2}$

然后利用最优阶次p，构造地震数据分数阶傅里叶变换的核矩阵，对最优阶次下的分数阶傅里叶算法进行离散，将地震数据从时空域变换到分数阶域中。

接着，在分数阶域中求取维纳滤波算子(即，预测算子P_j(f))，然后通过前面的等式(19)预测地震数据：

$S_K\left(f\right)=\underset{j=1}{\overset{L}{\mathrm{\Σ}}}{P}_j\left(f\right)S_{K-j}\left(f\right)$

具体地，在x方向做一个向前(后)预测滤波，即，在特定频率上，对于一个预测算子长度(即，时窗个数L的值)内的所有地震道，利用该预测算子长度内的所有地震道的地震数据分别与对应的预测算子相乘，再将该预测算子长度内的相乘结果(一共L个)相加作为下一地震道(即，第K地震道)的地震数据的预测结果。采用同样的方式可预测出其它地震道的地震数据。可采用前向预测或后向预测方式。在前向预测方式中，从预测算子长度内的最后一个地震道开始向前预测，直到预测出该预测算子长度内的第一地震道的地震数据；在后向预测方式中，从预测算子长度内的第一地震道开始向后预测，直到预测出该预测算子长度内的最后一个地震道的地震数据。

然后，利用分数阶傅里叶逆变换将分数阶域滤波(去噪)后的地震数据反变换到时空域中，并再次利用相关系数对地震数据进行拟线性反变换，恢复地震数据同相轴的本来面貌，以此来进行分数阶傅里叶变换的随机噪声压制处理。

图2C至图2F示出了随机噪声的压制效果，其中，图2C是针对图2B的数据进行p＝0阶时分数阶域去噪后的数据，图2D是针对图2B的数据进行p＝0.3阶时分数阶域去噪后的数据，图2E是针对图2B的数据进行p＝0.94阶时分数阶域去噪后的数据，图2F是针对图2B的数据进行p＝1阶时分数阶域去噪后的数据。

如图2C所示，当p＝0阶时，滤波后对随机噪声的压制效果最差。如图2D所示，当p＝0.3阶时，虽然压制了大部分随机噪声，但是损害了很多有效信号的波形。如图2F所示，当p＝1阶时，虽然随机噪声压制效果较好，但不能很好地保留有效信号特征。如图2E所示，当p＝0.94阶时，压制了大部分随机噪声，被淹没的有效信号得到很好的恢复。

从图2C至图2F可以看出：进行分数阶域预测去噪时，对于不同的阶数，取得不同的滤波结果；只有取得最优阶次，才能兼顾噪声压制和有效信号特征保持两个方面，由此取得较好的滤波效果。

图3A至图3F分别是图2A至图2F所示的理论模型数据中第10地震道的频谱。可从图2C至图2F的地震数据中分别抽取第10地震道进行频谱分析，得到图3C至图3F所示的频谱。

如图3C所示，当p＝0阶时，滤波后频谱效果最差。如图3D所示，当p＝0.3阶时，虽然压制了大部分随机噪声，但是损害了很多有效信号的频谱。如图3F所示，当p＝1阶时，虽然随机噪音压制效果好，但不能很好地保留有效信号频谱特征。如图3E所示，当p＝0.94阶进行分数阶域预测滤波时，压制了大部分的随机噪声，被淹没的有效信号频谱特征得到很好的恢复。从图3C至图3F再次可以看出：进行分数阶域预测去噪时，对于不同的阶数，取得不同的滤波结果；只有取得最优阶次，才能兼顾噪声压制和有效信号特征保持两个方面，由此取得较好的滤波效果。

图4A至图4C示出了地质情况更为复杂的实际地震数据的应用效果图，其中，图4A是实际叠加后的地震数据的剖面图，图4B是针对图4A的地震数据进行p＝0.2阶时分数阶域去噪后的地震数据的剖面图，图4C是针对图4A的地震数据进行p＝0.92阶时分数阶域去噪后的地震数据的剖面图。

对于地质情况更为复杂的实际地震数据，如图4A所示，虽然实际地震数据本身的复杂性导致可预测性不高，地震数据的规律性差，但是可以通过一个合适的旋转角度，找到信号最集中的最优阶次进行分数阶域预测滤波处理。

旋转角度与阶次的关系如下：a＝pπ/2，a为旋转角度，p为分数阶次，这是分数阶傅里叶变换的旋转特性。

找到一个合适的旋转角度就等同于找到一个合适的分数阶次，而寻找合适的分数阶次就是最优阶次的求取方法，可按照前面的等式(8)至等式(13)来求取最优阶次。

从噪声压制效果图4B和图4C可以看出，FRFT预测滤波去噪取得了预期的滤波效果，并且在分数阶p＝0.92(图4C)的滤波效果比p＝0.2(图4B)时更好，噪声得到更好的压制，同相轴的连续性也更好，有效信号也突出，也能更好地保持实际地震记录所反映的地层信息和地质特征。

本发明根据地震信号在不同阶次下的分数阶域具有不同的广义时间-带宽积(Generelized Time-Bandwidth Product，GTBP)的原理，将地震信号变换到不同的分数阶域下，然后求取每个分数阶域下的广义时间-带宽积，最后根据带宽积越小信号的时频聚集度越高的理论搜索出最小的GTBP，该GTBP对应的阶次就是分数阶域最优阶次。进而在该最优阶次下进行随机噪声消除，有利于更好地提高地震信号的分辨率，并且损失更少的有效信号。

本发明采用了分数域预测滤波技术，将地震数据变换到分数阶域中，并结合预测滤波方法，在分数空间域进行随机噪声的衰减。这有利于更好地提高地震信号的分辨率且损失更少的有效信号，有利于预测滤波效果与信号特征保持的最佳折中。

虽然本发明是参照其示例性的实施例被具体描述和显示的，但是本领域的普通技术人员应该理解，在不脱离由权利要求限定的本发明的精神和范围的情况下，可以对其进行形式和细节的各种改变。

Claims (7)

1.一种基于分数阶傅里叶变换的随机噪声压制方法，包括：

(1)针对采集的地震数据，利用正交多项式逼近法，采用最大相似系数准则，计算描述地震数据中同相轴走势的相关系数；

(2)利用计算的相关系数对地震数据在空间上进行拟线性化变换；

(3)根据分数阶傅里叶变换的最优阶次求取法，求取地震数据中信号最集中的分数阶域最优阶次；

(4)利用最优阶次下的分数阶傅里叶变换将地震数据从时空域变换到分数阶域；

(5)求解地震道的预测算子，利用所述预测算子预测出地震道的地震数据，由此来消除无规律分布的随机噪声；

(6)利用分数阶傅里叶逆变换将去噪后的地震数据从分数阶域变换到时空域；

(7)利用计算的相关系数对地震数据进行反拟线性变换，完成地震数据的随机噪声压制，得到消除随机噪声后的有效信号的地震剖面，

其中，在步骤(3)中，根据下面的等式来计算所述最优阶次：

$\tan \left(2\mathrm{p\π}\right)=\frac{4{\mathrm{\μ}}_0(w_0-w_1)}{4{\mathrm{\μ}}_0^2-{(w_0-w_1)}^2}$

其中，p表示分数阶，μ₀＝(w₀+w₁)/2-w_0.5，w₀是分数阶p＝0阶时地震信号在分数阶域的二阶矩，w₁是分数阶p＝1阶时地震信号在分数阶域的二阶矩，w_0.5是分数阶p＝0.5阶时地震信号在分数阶域的二阶矩，

其中，步骤(4)包括：从连续傅里叶变换的特征函数出发，通过对特征函数进行离散化近似和正交投影，得到一组与特征函数形状相似的离散傅里叶变换矩阵的正交化离散特征函数的特征向量；根据连续分数阶傅里叶变换的核函数谱分解表达式，构造离散分数阶傅里叶变换核矩阵；将地震数据与最优阶次下的分数傅里叶变换核矩阵相乘，得到变换后的分数阶域下的地震数据。

2.根据权利要求1所述的随机噪声压制方法，其中，在步骤(1)中，利用正交多项式逼近法对表示叠加地震数据的反射同相轴的时间多项式进行变换，逐项扫描变换后的时间多项式的系数，将扫描得到的系数作为所述相关系数。

3.根据权利要求2所述的随机噪声压制方法，其中，在逐项扫描变换后的时间多项式的系数的过程中，每次扫描新的系数时使先前确定的系数保持不变。

4.根据权利要求3所述的随机噪声压制方法，其中，在叠加地震数据的剖面上设置多个时窗，在每个时窗内通过多地震道互相关扫描来确定该时窗内的时间多项式系数。

5.根据权利要求1所述的随机噪声压制方法，其中，在步骤(5)中，根据地震数据中相干信号可预测，随机噪声不可预测，并且同相轴为线性，来求解所述预测算子。

6.根据权利要求5所述的随机噪声压制方法，其中，对预测算子的求取进行线性变换，通过复数维纳滤波求得预测算子。

7.根据权利要求6所述的随机噪声压制方法，其中，在利用所述预测算子预测出地震道的地震数据的步骤中，在特定频率上，对于一个预测算子长度内的地震道，利用所述预测算子长度内的所有地震道的地震数据分别与对应的预测算子相乘，再将所述预测算子长度内的相乘结果相加作为下一地震道的地震数据的预测结果。

Images