一种基于时空分数阶滤波的地震资料噪声压制方法与流程

本发明属于地球物理勘探领域，涉及一种地震资料噪声压制方法，特别涉及一种基于时空分数阶各向异性滤波的地震资料噪声方法。

背景技术：

在地震勘探过程中，由于地面的微震、仪器在接收或处理过程中的噪声、复杂地表和地形条件，如山地、沙漠和戈壁等因素，均会对接收地震记录产生干扰。这些干扰会降低地震资料的信噪比及质量，使得后续地震资料解释和参数反演难度变大。因此，提高地震资料的信噪比，改善地震资料品质，是地震资料处理中需要解决的关键问题之一。由于复杂的地质沉积作用，地震信号的显著特点是层状特征明显，且含有丰富的不连续结构特征，如断层、河道等。这些结构信息对后续地震资料解释尤为重要。因此，地震资料在进行噪声压制的过程中，应充分保持这些结构特征。

按照不同特征，地震资料中的噪声可分为不同类型。在实际地震资料处理中，应针对具体问题区分有效信号与噪声，并采取相应的噪声压制方法。常规的地震资料噪声压制方法主要包括f-x域预测滤波、kl变换、拉东变换、小波变换去噪、contourlet变换去噪、curvelet变换去噪、svd分解等。为了在压制地震资料噪声的同时，充分保护地震资料的结构特征，luo等基于次序统计思想提出一种地震资料保边滤波方法。该方法利用kuwahara窗分析技术，将待分析点周围方差最小的窗口内的均值代替该点的值。albinhassan等将该方法推广到三维地震资料保边滤波。hoeber等综合分析了基于次序统计思想的滤波方法，提出利用级联滤波器和尺度自适应高斯滤波器来压制三维地震资料噪声。fehmers和hocker于2003年首次将各向异性扩散滤波方法应用于地震资料保边缘噪声压制，利用两种不同尺度的结构张量来度量地层结构的不连续性，从而控制扩散滤波在地震结构体边缘的滤波方向和程度。孙夕平等研究了有限差分法求解扩散方程的特点，提出具有最优旋转不变性的各向异性扩散滤波方法。王旭松和杨长春在非线性各向异性扩散方程中引入二阶导数项来提高对地震结构的保护能力。lavielle等改进了三维地震资料噪声压制中扩散张量特征值的构造方法，将线型扩散和面型扩散相结合，使得扩散过程不仅能对面型反射结构进行滤波，同时能保持和增强地震断层结构。张尔华等提出利用结构张量计算地震数据的不连续性参数，以控制各向异性扩散滤波器的边缘保持能力。

基于偏微分方程的地震资料噪声压制方法最早来源于图像处理领域。自20世纪90年代以来，扩散滤波方法得到迅速发展，并在图像去噪、分割和增强等领域得到广泛应用。perona和malik在1990提出一种非线性扩散滤波方法

式中，x为空间坐标，ω为图像所在空间区域，u0为原始含噪图像，u(t,x)为t时刻扩散滤波结果，▽u为图像的梯度，▽·为散度算子，g(·)为扩散率函数，满足非负性和单调递减性，同时满足g(0)＝1，g(∞)＝0。该方法根据图像的局部特征确定扩散系数，引导扩散过程在图像较为平坦的区域具有较快的扩散速度，在图像边缘附近具有较慢的扩散速度，从而兼顾边缘结构信息和噪声压制。由于扩散过程有标量的扩散率函数控制，其本质上依然是各向同性的。

weickert在1996年提出了基于结构张量分析的各向异性扩散滤波方法，

式中，d为扩散张量。该方法根据扩散张量控制沿不同方向上的扩散速率和程度，考虑了图像的局部一致性结构，使得平滑过程主要沿着各向异性较弱的方向进行，实现了各向异性扩散滤波。

以上现有地震资料噪声压制方法存在以下缺点：

(1)常规扩散滤波方程本质上是一个热传导方程，其解呈指数衰减，在压制噪声的同时，易损害有效地震信号的振幅；

(2)当地震噪声干扰严重时，扩散张量难以得到有效估计，无法正确引导扩散滤波的方向和程度。

技术实现要素：

本发明的目的在于克服上述现有技术的缺点，提供一种基于时空分数阶各向异性扩散滤波方程的地震资料噪声压制方法，该算法基于时间和空间分数阶微积分理论，实现对扩散滤波方向和强度的精细控制，引导扩散过程沿着地震资料各向异性较弱区域进行，能有效压制叠前和叠后地震资料中的随机噪声，增强地震同相轴的空间一致性，并且能够保护有效信号和断层、地层边缘等细节结构，该技术方案易于实现，可操作性强。

本发明的目的是通过以下技术方案来解决的：

本发明是通过求解一个时间-空间分数阶各向异性扩散滤波方程实现二维叠前或叠后地震资料噪声压制。连续情况下，分数阶扩散滤波方程的具体形式为

式中，t为时间，u0(x)为二维含噪地震记录，u(t；x)是经过扩散时间t后的地震记录，x＝(x,y)，g(·)为扩散率函数，α(t；x)是时间扩散分数阶阶次函数，α∈[1,2)，|·|表示矩阵frobenius范数，和分别为x方向和y方向分数阶微分算子，和分别为和的伴随算子。

特别地，空间分数阶微分算子和的计算公式为

式中，β和γ分别为空间分数阶阶次；ξ1和ξ2分别为水平方向和垂直方向波数；i为虚数单位，即是u的二维傅里叶变换。

riemann-liouville意义下的时间分数阶微分算子的计算公式为

式中，γ(·)为gamma函数，n为整数，且满足n-1＜α≤n。

进一步，本发明用于地震资料噪声衰减的时空间分数阶各向异性扩散滤波方法，包括以下步骤：

1)记二维叠前或叠后地震数据为u(x)，其中，x＝(x,y)∈[0,x]×[0,y]，x为列指标，代表偏移距方向，y为行指标，代表时间方向。对u沿水平方向和垂直方向进行离散，记为u(p,l)，u∈r^n×m，离散间隔δx＝x/m,δy＝y/n，其中p∈{1,2,…,m}，l∈{1,2,…,n}；

2)建立时间-空间分数阶各向异性扩散滤波方程；

3)确定时间迭代步长δt和最大迭代时间t；

4)确定水平方向滤波分数阶阶次β和垂直方向滤波分数阶阶次γ；

5)对地震数据做高斯滤波得到uσ(p,l)，然后确定自适应变时间分数阶阶次函数α(t；p,l)；

6)利用快速离散傅里叶变换计算水平方向和垂直方向分数阶微分和然后计算和的伴随算子和

7)分别计算横向和纵向扩散率函数和

8)利用预测-校正算法求解时间分数阶微分方程；

9)进行迭代k＝k+1；判断是否满足终止条件，若不满足，则返回步骤4)；否则，执行步骤9)；

10)输出最终的地震记录u^clean(p,l)。

进一步，上述步骤2)中，建立时间-空间分数阶各向异性扩散滤波方程如下

式中，t为时间，u0(x)为二维含噪地震记录，u(t；x)是经过扩散时间t后的地震记录，x＝(x,y)，g(·)为扩散率函数，α(t；x)是时间扩散分数阶阶次函数，α∈[1,2)，|·|表示矩阵frobenius范数，和分别为x方向和y方向分数阶微分算子，和分别为和的伴随算子。

特别地，空间分数阶微分算子和的计算公式为

式中，β和γ分别为空间分数阶阶次；ξ1和ξ2分别为水平方向和垂直方向波数；i为虚数单位，即是u的二维傅里叶变换。

riemann-liouville意义下的时间分数阶微分算子的计算公式为

式中，γ(·)为gamma函数，n为整数，且满足n-1＜α≤n。

上述步骤5)中，对当前时刻的地震数据u(t；p,l)做高斯滤波得到uσ(t；p,l)，即

式中，为卷积符号，kσ(p,l)为标准差为σ的二维空间高斯核函数

求取uσ(t；p,l)的归一化梯度▽uσ(t；p,l)，然后构造自适应变时间分数阶阶次函数α(t；p,l)如下

式中，λ为阈值参数，λ＞1以保证α∈(1,2)，这里选取λ＝1.5。

上述步骤6)中，离散情况下，利用分数阶空间差分和逼近和具体的计算公式为

式中，β和γ分别为空间分数阶阶次；dft和idft分别表示离散傅里叶变换和离散逆傅里叶变换；ωx∈{0,1,…,m-1}为水平方向离散化空间频率；ωy∈{0,1,…,n-1}为垂直方向离散化空间频率；i为虚数单位，即

和的伴随算子和的计算公式为

式中，“.”为矩阵乘法符号，kx和ky为对角矩阵

上述步骤7)中，扩散率函数为

上述步骤8)中，为了便于叙述，记右端的空间分数阶差分项为

则时间分数阶微分方程可转化为如下积分方程形式

式中，γ(·)为gamma函数。

本发明具有以下有益效果：

本发明提出一种新的有效压制叠前和叠后地震资料中随机噪声干扰和采集“脚印”噪声方法。本算法利用时间和空间分数阶各向异性扩散滤波方程，实现了对扩散方向和强度的精细控制，克服了传统扩散方程易损害有效信号振幅的缺点。本发明算法可有效引导扩散滤波沿着各向异性较弱区域进行，从而达到保持地震资料纹理结构、断层、裂缝等地层边缘结构，有效地改善地震同相轴的空间连续性，增强地震同相轴的空间一致性，特别是对噪声干扰严重的深层弱反射信号亦取得较好的去噪效果。

附图说明

图1是本发明流程示意图；

图2是某陆上实际地震资料去噪前后对比显示图；

图3是去噪前后地震记录局部放大显示图；

图4是去噪前后地震记录归一化振幅谱对比图。

具体实施方式

下面结合附图对本发明做进一步详细描述：

提高地震资料的信噪比，改善地震资料品质，是地震资料处理中需要解决的关键问题之一。在地震资料在进行噪声压制的过程中，应充分保持地震资料的纹理特征与不连续结构特征。本发明提供一种基于时间-空间分数阶各向异性扩散滤波方程的地震记录噪声压制方法。

本发明的物质基础是地震采集设备采集到的叠前地震资料，或者结果预处理后得到的叠后地震资料。本发明是通过求解一个时间-空间分数阶各向异性扩散滤波方程实现二维地震资料噪声压制。连续情况下的扩散方程的具体形式为

式中，t为时间，u0(x)为原始含噪地震记录，u(t；x)是经过扩散时间t后的地震记录，x＝(x,y)，g(·)为扩散率函数，α(t；x)是时间扩散分数阶阶次函数，α∈[1,2)，|·|表示矩阵frobenius范数，和分别为x方向和y方向分数阶微分算子，和分别为和的伴随算子。

特别地，空间分数阶微分算子和的计算公式为

式中，β和γ分别为空间分数阶阶次；ξ1和ξ2分别为水平方向和垂直方向波数；i为虚数单位，即是u的二维傅里叶变换。

riemann-liouville意义下的时间分数阶微分算子的计算公式为

式中，γ(·)为gamma函数，n为整数，且n-1＜α≤n。

本发明的基于时空分数阶扩散方程的地震资料噪声压制算法的框架如图1所示，具体步骤分别为：

1)记二维叠前或叠后地震数据为u(x)，其中，x为二维连续变量，x＝(x,y)∈[0,x]×[0,y]，x为列指标，代表偏移距方向，y为行指标，代表时间方向。对u沿水平方向和垂直方向分别进行离散，记为u(p,l)，u∈r^n×m，水平和垂直方向离散间隔分别为δx＝x/m,δy＝y/n，其中p∈{1,2,…,m}，l∈{1,2,…,n}。

2)确定时间迭代步长δt和最大迭代时间t。为了保证偏微分方程数值求解的稳定性，时间迭代步长δt应满足稳定性条件，这里选取δt＝0.05，其大小满足l2范数意义下有限差分稳定条件。最大迭代时间t应由地震资料的噪声水平决定，即当地震记录噪声干扰严重时，取较大的t；反之，取较小的t值。

3)确定水平方向滤波分数阶阶次β和垂直方向滤波分数阶阶次γ。空间分数阶β和γ取值大小决定了沿水平方向和垂直方向的滤波强度。β和γ的选取依赖于将要处理的地震资料类型和结构。对于叠后地震数据，若地层结构起伏较小，同相轴倾角小且横向连续性好，则应取较大β值和较小的γ值；反之，β和γ应取较小的值，以保护丰富的地层结构，特别是断层、裂缝等不连续结构。对于叠前地震数据，由于反射同相轴曲率差异较大，应取较小的β和γ。

4)对当前时刻的二维地震数据u(t；p,l)做高斯滤波得到uσ(t；p,l)，即

uσ(t；p,l)＝kσ(p,l)*u(t；p,l)

式中，“*”为二维卷积符号，kσ(p,l)为离散化后标准差为σ的高斯核函数

求取uσ(t；p,l)的归一化梯度▽uσ(t；p,l)，然后构造自适应变时间分数阶阶次函数α(t；p,l)如下

式中，λ为阈值参数，λ＞1以保证α∈[1,2)，这里选取λ＝1.5。

5)利用快速离散傅里叶变换计算水平方向和垂直方向空间分数阶微分和然后计算和的伴随算子和离散情况下，可利用分数阶空间差分和来逼近空间微分算子和具体的计算公式如下

式中，β和γ分别为空间分数阶阶次；dft和idft分别表示离散傅里叶变换和离散逆傅里叶变换；ωx∈{0,1,…,m-1}为水平方向离散化空间频率；ωy∈{0,1,…,n-1}为垂直方向离散化空间频率；i为虚数单位，即

空间分数阶差分算子和的伴随算子和的计算公式为

式中，“·”为矩阵乘法符号，kx和ky的表达式如下

式中，diag{·}表示由大括号内的向量张成的对角矩阵。

6)分别计算横向和纵向扩散率函数和通过引入扩散率函数扩散率函数，实现了非线性各向异性扩散滤波。本发明选用的扩散率函数如下所示

7)利用预测-校正算法求解时间分数阶微分方程。为了便于叙述，记右端的空间分数阶差分项为d(p,l；β,γ)，其表达式为

则时间分数阶微分方程可转化为如下积分方程形式

式中，γ(·)为gamma函数。对于上述积分方程，可采用具有较高数值精度和较好数值稳定性的预测-校正算法进行求解。

8)进行迭代t＝t+δt；判断是否满足终止条件t＝t，若不满足，则返回步骤4)；否则，执行步骤9)。

9)输出最终的地震记录u^clean(p,l)。

实际地震资料测试

本节将本发明提出的基于时空分数阶滤波的地震资料噪声压制方法应用于某陆上油田实际叠后二维地震数据，验证该算法的有效性。图2(a)为原始三维地震数据体抽取的一条纵向测线，共2000道，每道共901个采样点，时间采样间隔为2.0ms，截取时间范围为1.8～3.0s。该勘探区块位于黄土塬，地表复杂，地震剖面采集脚印干扰严重。图2(a)和(b)分别为利用本发明提出算法去噪后地震剖面和相应的差剖面。可以看到，地震剖面中的采集脚印干扰得到有效压制，地震同相轴的空间连续性得到显著改善。截取cdp601～800、时间1.5～2.7s范围内去噪前后地震剖面进行局部放大显示，如图3所示。从结果可以看出，地震剖面的同相轴连续性得到极大改善，边界清晰，且残差剖面中未见明显的同相轴信息，表明该方法在有效压制噪声的同时，有效地震信号的振幅得到较好保持。图4中虚线和实线分别为去噪前后地震剖面的归一化平均振幅谱。可以看出，去噪后地震数据的频谱形态在有效频带内基本没有改变，表明有效地震信号和能量得到较好保持，验证了本发明算法的有效性。