基于Jacobi算法求解结构张量的置信扩散滤波方法

饶 溯， 李录明， 刘力辉， 胡 滨， 冯 鑫

(1.中海石油国际能源服务(北京)有限公司，北京 100028；2.油气藏地质及开发工程国家重点实验室(成都理工大学)，成都 610059；3.北京诺克斯达石油科技有限公司，北京 100192)

0 引言

地震资料去噪在叠前和叠后资料应用中均非常重要。压制叠前角道集的噪声一直是地震数据叠前处理中的重要问题，其对叠前振幅属性分析及应用有重要影响(如噪声压制的好坏将直接影响叠前地震资料振幅的应用(如叠前AVA反演)效果)。而压制叠后地震资料的噪声也将更有利于断层断点准确位置的解释、层位的自动追踪及构造圈闭的准确落实。

地震勘探中噪声一般包括随机噪声和相干噪声。通常相干噪声属低频信号，而随机噪声属高频信号，并具有统计特性。噪声的压制通常采用滤波的方法，针对噪声的不同特征及叠前、叠后处理，采用多种滤波方法，其中具有代表性方法包括奇异值分解(SVD)[1]、F-xy域图像特征值滤波[2]、F-k域滤波[3]、τ-p域滤波[4]、神经网络等一些模糊滤波去噪[5]和扩散滤波[6-9]等。

Koenderink[10]发现扩散滤波等价于一个热传导的物理过程。而非线性各向异性扩散方程(P-M方程)是由Perona等[11]把它作为图像处理中有选择性地保留边界信息的工具而首先提出。P-M方程有可能是病态的， You[12]通过正则化的方法消除了这种病态性。

Weickert[13-14]对P-M方程进行了改进，引进了结构张量D，并提出最优旋转不变性的新型算法;孙夕平等[6]对Weickert的新型算法公式采用有限差分优化求解并设计近似旋转不变性模板代替常规的平滑算子；王益艳等[15]在非线性各向异性扩散滤波基础上构建梯度和切线方向上的扩散滤波系数，并加入了数据保真项。他们所提出的方法，笔者都将其归为常规基于结构张量的扩散滤波方法，而其中结构张量的设计和求解方面还存在一定的改进空间，笔者是在常规方法的基础上进行改进。

在非线性各向异性扩散滤波基础上，采用Jacobi算法来求解结构张量D中的特征向量，并引入线状置信度Cline来共同来构建扩散滤波系数λ，最后来求解扩散方程，实现置信扩散虑波。本文方法相比常规基于结构张量的扩散滤波方法，在构建和求解结构张量D方面及构建扩散滤波系数有所差异，最终效果较常规基于结构张量的扩散滤波方法也所提高：随机噪声压制效果较好、地震同相轴连续性更强，断层成像更加清晰准确，明显改善资料品质。

1 常规方法

1.1 非线性各向异性扩散滤波原理

非线性各向异性扩散可模拟一个热传导的物理过程，它在没有产生和损失热量的前提下平衡温度。将扩散应用于地震道集，用地震振幅代替热量，在扩散过程中时间是一个固定参数。扩散进行的时间越长，地震道集滤波程度则越强，所以滤波过程相当于是提高信噪比的过程。地震振幅的梯度变化相当于热传导的通量，其扩散方程可以表示为[13-14]：

(1)

其中:D(|▽uσ|)为结构张量;|▽uσ|为正则化后的▽uσ，正则化过程相当于高斯滤波;t为扩散时间(也可计为迭代次数);u0(x,y)为原始地震道集；u(x,y,t)为迭代次数为t次后的滤波地震道集。

▽uσ=▽(Kσ*u(x,y,t))

(2)

(3)

式中:σ为尺度函数，可以根据噪声的均方误差自适应调节。

1.2 最优旋转不变算法

公式(1)可在正交坐标系下展开成如下形式：

(4)

其中结构张量可表示为式(5)。

(5)

公式(4)可表示为沿着x、y两个方向(即地震道集道号和时间方向)上进行扩散，而最优旋转不变算法是需要对公式(5)的结构张量D进行一定程度上的旋转，使其沿着地震道集局部振幅值变化最快和最慢方向上扩散，如王益艳等[15]再用扩散系数来控制在不同变化方向上平滑程度，即常规基于结构张量的扩散滤波方法。

公式(4)旋转及引入扩散系数之后可转换成如下形式：

(6)

式中:λ为扩散系数，扩散系数小则表示扩散较慢，反之扩散系数大则表示扩散较快；η表示地震道集局部振幅值变化较快和变化剧烈处，其扩散程度较慢；ξ表示地震道集局部振幅值变化较慢和变化均匀处，其扩散程度较快。

2 本文方法

2.1 结构张量的改进

为了增强同向轴连续性和提高压制噪声的效果，在本文方法中对结构张量的设计进行了改进，并引入线状置信度来改进和自适应调节扩散系数。

该方法以公式(4)为基础，并对公式(5)中的结构张量改进为如下形式：

(7)

指示地震剖面局部振幅值变化剧烈位置处扩散系数较小，记为λ1，指示地震剖面局部振幅值变化均匀段扩散系数较大，记为λ2，有如下形式:

λ1=α

(8)

(9)

(10)

最后根据特征值u、扩散系数λ和线状置信度Cline可以将公式(4)转换成形式如下：

(11)

当η值和ξ值比较接近时，从公式(11)可以看出两个方向上扩散速率较为一致。公式(11)中振幅剧烈变化处位置值η值必然大于振幅平稳变化处位置值ξ，而当ξ值远小于η值时即代表振幅平稳位置有很快的扩散速率。

2.2 Jacobi算法原理

考虑到结构张量D矩阵是一个2×2阶的对称矩阵，采用2×2的对称Schur分解的Jacobi算法来求解效率较快。Jacobi算法的思想是逐步地减小特征分量矩阵v中非对角位置上值的大小[17]，也可称为Givens旋转变换：

(12)

2.3 方法的实现步骤

公式(11)可化为离散差分形式如下：

(13)

本文方法即是对公式(13)进行求解，完成对输入原始地震道集u0的t次迭代运算。主要可分为以下几步：

1)先对原始地震道集u0进行高斯滤波，这个高斯滤波相当于是正则化，然后分别求道号和时间方向上的二阶导数和混合偏导数v。

2)根据公式(7)中构建特征分量v，并用Jacobi算法来求解其特征值u，再根据公式(10)计算线状置信度Cline，将Cline代入扩散系数公式(9)中计算出λ1和λ2。

3)将原始地震道集u0、及其在道号和时间方向上的二阶导数λ1和λ2代入到公式(13)，即完成一次扩散滤波计算过程，有扩散滤波计算结果ut。

3 模型试验及效果分析

3.1 理论模型正演及去噪试验

表1给出了单层满足四类AVA曲线的各向同性介质理论模型，图1(a)～图1(c)为由表1中给出的模型组合在一个道集上的AVA正演角道集，入射角为1°～45°,AVA正演角道集采用合成记录正演。从图1(c)扩散滤波方法的结果中仍然可见四类AVA曲线地震振幅特征，即该方法有一定程度的保幅效果，且地震同相轴也较为连续。

图1 含四类AVA曲线的正演模型角道集Fig.1 Forward modeling of angle gather with four kinds of AVA curves(a)原始正演模型；(b)加20%随机噪声模型；(c)20%随机噪声模型去噪后

表1 四类AVA曲线的各向同性介质理论模型

图2 SNR和MSE随迭代次数变化的曲线Fig.2 Curves of SNR and MSE with iteration times

3.2 去噪试验的信噪比及误差分析

图像滤波和信号降噪通常采用SNR(信噪比)和MSE(均方误差)两项指标来进行衡量。本文分别以含10%和20%随机噪声的模型去噪试验为例，对比分析了SNR和MSE曲线的变化。如图2所示，随着迭代次数的增加，SNR值合理的变大，MSE值有效的减小，在满足一定迭代次数的时候，10%和20%噪声下的SNR值和MSE值分别达到极值，即扩散到一定程度，扩散过程停止，地震振幅滤波效果将不再发生变化。

3.3 去噪试验的时频分析

去噪试验压制噪声效果的好坏也可通过时频域的响应特性来观察。如图3(a)～图3(c)所示，对比了加噪正演模型带通滤波和本文方法去噪后的时频图，其中时频图采用广义S变换方法计算[18]。通过对比图3(b)和图3(c)的压制噪声效果好坏，可观察到图3(c)中红色方框内时频能量团更为聚焦，其对应的时间刻度更为准确。图3(b)为带通滤波方法的结果，它仅仅压制了频带范围外的噪声，对于有效信号频带范围内噪音没有起到压制作用。而本文方法可以消除部分有效信号频带范围内的噪声，但仍有部分噪声信号未被消除。对比图3(b)和图3(c)的浅层黑色方框位置，可见带通滤波方法仍可见时频能量团的信息，而本文方法中弱的时频能量已经掩盖在噪声中。这原因可能归结于扩散滤波方法仅是在时间域上进行相干增强滤波，有效的强地震同向轴可明显增强，但要保证弱地震同向轴信息可尝试变换域内去除噪声，可参考小波域扩散滤波方法，但其计算速度较慢[19-20]。

图3 含四类AVA曲线的正演模型的时频图Fig.3 Forward modeling of the time-frequency diagrams with four kinds of AVA curves(a)原始模型；(b)加20%随机噪声模型带通滤波结果；(c)加20%随机噪声模型扩散滤波结果

4 实际资料处理

为了验证基于Jacobi算法求解结构张量的置信扩散滤波方法，在实际地震资料处理中的效果，将该方法应用在Geoscope软件平台上进行叠后和叠前地震数据去噪处理[21]。本文方法源自于图像去噪处理，对图像质量没有特定要求。由方法原理可知，该方法只需输入地震数据，设置控制迭代次数的扩散变换率阈值δ0、与扩散系数有关的系数α、高斯窗尺度σ的值，其对其他条件并无特别要求，故对低、高信噪比的二维和三维地震资料均可适用。根据地震资料情况，仅去噪效果有一定的差别；对于信噪比极低的地震资料去噪效果本身就较差，而对于信噪比极高的地震资料去噪效果可能会不太显著。三项关键参数δ0、α、σ值的选取参考方法设计原理和试验测试的结果，δ0和α取值都不宜过大且一般都取接近于零的极小值，略微调节该两参数的大小对实际地震资料处理的影响不大，σ高斯窗尺度选取3或5。

图4(a)～4(b)展示了邻近几道的叠前道集去噪对比图。Geoscope软件中采用作业流形式对逐个道集进行处理。从图4(a)可以看出，这几个道集含有较强的随机噪声，信噪比不高。图4(c)为去噪前后几个道集的差值剖面，与随机噪音比较符合。说明经过本文方法去噪后，该道集中的随机噪声得到了有效地压制，同相轴也变得更加清晰，但是道集上仍有部分有规则形态的斜干扰噪声不能很好去除。图5(a)～图5(c)展示了叠后地震资料的常规方法和本文方法去噪对比剖面图。从中可知，本文方法压制噪声效果明显优于常规方法的效果：更加有效地压制了随机噪声，更加增强了同相轴连续性，更明显地突出了断层位置，削弱了同相轴的蚯蚓状现象。

实际资料处理表明，本文方法不仅在叠前去噪处理中为后续的AVO属性技术和叠加处理提供了更有力地保证，而且在叠后地震数据体处理将更有利于层位自动追踪和断层位置精确解释。

此外，本文方法和常规方法的中δ0、σ值可选同一值，仅α值为本文方法中特有取值。当α和σ值恒定时，本文方法同其常规方法一致，其计算速度和效率仅与控制迭代次数的扩散变换率阈值δ0和待处理的地震数据体大小有关，δ0越小则迭代次数越多、地震数据体越大则计算量增加，导致计算效率有所降低。当δ0一定时，本文方法仅与地震数据体大小有关。在Intel至强E3系列CPU处理器、64G内存和64位操作系统的个人工作电脑上使用Geo-scope软件平台windows版本，计算1.7GB的叠后三维地震数据体的去噪结果，本文方法计算用时约为38 min，常规方法计算用时约为5 min，一定程度上认为本文方法的计算用时在可接收范围，而对较为庞大的叠前地震数据处理情况，可尝试改进算法及采用GPU计算来达到提高计算效率的目的。

图4 扩散滤波在叠前地震资料中的应用效果Fig.4 Application results of diffusion filtering in pre-stack seismic data(a)原始叠前道集；(b)本文方法去噪后叠前道集；(c)去噪前后的差值剖面

图5 扩散滤波在叠后地震资料中的应用效果Fig.5 Application results of diffusion filtering in post-stack seismic data(a)原始叠后剖面图；(b)常规基于结构张量的扩散滤波方法去噪后剖面图；(c)本文方法去噪后剖面图

5 结论

笔者在非线性各向异性扩散滤波基础上，利用Jacobi算法求解结构张量的特征值，并建立线状置信度来计算扩散滤波系数，提出了一种基于Jacobi算法求解结构张量的置信扩散滤波方法：

1)该方法对叠前道集的AVA曲线特征有一定保幅性，保证了后续AVO属性分析和叠加处理过程的顺利进行。

2)该方法不仅压制随机噪声效果好于常规基于结构张量的扩散滤波方法，而且在同相轴的连续性增强和突出断层位置方面都更好的效果。

3)该方法可为层位自动追踪和断层位置精确解释，提供高信噪比资料。

4)该方法随迭代次数的增加，SNR合理的变大，MSE有效的减小，最终达到一定迭代次数，扩散变化率极低，扩散过程停止，该方法将不再有滤波效果。

5)该方法相比帯通滤波方法，部分有效信号频带范围内的噪声能被有效压制，强地震同向轴位置的时频能量团更为聚焦，而弱地震同向轴的时频能量团掩盖在噪声中。

6)该方法的参数选取较为简单，自适应计算扩散系数，在软件应用中不用过度注重去调节和选取参数。