结合边缘检测和非局部均值的地震数据降噪

李晓璐，周亚同，何静飞，翁丽源，李书华

（1.河北工业大学 电子信息工程学院，天津 300401；2.天津市第三烟草专卖局 天津 300131）

0 引言

油气开采过程中，地震数据通常含有大量噪声，会造成有效数据丢失[1-2]。因此，要对地震数据进行降噪处理。去除地震数据中的高斯随机噪声，对改善数据信噪比和分辨率具有现实意义[3-4]。常用的空间域降噪算法有均值滤波[5]、中值滤波[6]、高斯滤波[7]、维纳滤波[8]等，变换域降噪算法有小波变换[9]、曲波变换[10]、高斯混合模型[11]等。而地震数据细节丰富且有冗余信息，普通空间域降噪算法会平滑边缘细节，变换域降噪运算量大。黄英等[12]提出将非局部均值算法用于地震数据降噪，利用数据中特有的冗余信息，实现了有效降噪。常德宽等[13]利用高斯模型代替非局部均值算法中使用的全部相似数据块加权平均，实现了对高斯随机噪声的抑制。但由于地震数据结构信息丰富[14]、降噪要求高，上述改进算法仍然不能有效处理地震数据中的高斯随机噪声。

针对以上问题，本文将八方向Sobel算子与非局部均值算法[15]结合，并应用于地震数据降噪。该算法可准确检测同相轴边缘，并利用地震数据的相似性，改进非局部均值算法的权值函数，扩大结构相差较小邻域的权值，缩小结构相差较大邻域的权值。与NLM[16]、Sobel2-NLM、Canny-NLM[17]等算法相比，Sobel8-NLM用于地震数据降噪时，其峰值信噪比和结构相似度更高。

1 八方向Sobel算子的同相轴边缘检测

设(m,n)是地震数据中的任意一点，f(m,n)表示点(m,n)的地震数据，以(m,n)为中心的地震数据邻域矩阵为

Sobel算子可表示为

式中：hx检测水平边缘；hy检测垂直边缘。

Gx和Gy分别是像素点(m,n)为中心的邻域矩阵v[N(m,n)]与hx、hy进行卷积运算得到的在水平方向和垂直方向的梯度值，⊗表示卷积运算，中心像素点的梯度幅值为：

中心点的梯度方向为

Sobel算子利用像素点上下和左右邻点的灰度加权平均，并用边缘点取极值的方法提取边缘信息，其算法步骤如下：

1）分别将地震数据中的数据点与方向模板做卷积运算，模板中心点与矩阵中心点重合，得到2个方向的梯度值；

2）取上述2个梯度值的均方根距离，并将其作为对应点的新灰度值；

3）选择阈值Th，若新灰度值大于Th，则该点为地震数据的边缘点，将边缘点集合即得到地震数据的同相轴边缘。

Sobel算子边缘检测用两方向梯度模板，对复杂同相轴边缘检测效果并不理想。本文使用八方向Sobel边缘检测算子对地震数据进行同相轴边缘检测，从2个方向增加到8个方向，分别为0°、45°、90°、135°、180°、225°、270°、315°，算子模板如图1所示。

图1 八方向边缘检测算子模板Fig.1 Eight edges detection operator template

将地震数据中每一个点与八方向模板做卷积运算，得到八方向上的梯度矩阵，并对每一点取运算结果的最大值为梯度值，取最大值对应的模板方向作为该点的边缘方向，最后聚合得到整个梯度幅值图像G(x,y)。

方向模板的增加使得Sobel算子可以更加有效地提取到同相轴边缘，为规避噪声对其造成的影响，需要选取合适的阈值对G(x,y)进行二值化处理，选取步骤如下：

1）选取G(x,y)中所有点的灰度均值作为阈值初值T0；

2）将G(x,y)中所有点的灰度值与T0比较，将G(x,y)分为2个部分，G1是灰度值大于T0的点集合，G2是灰度值小于T0的点集合；

3）分别计算这个两个部分的平均灰度值，记为u1和u2；

4）定义一个新的阈值T1=(u1+u2)∕2；

6）利用Th对G(x,y)进行二值化处理：当G(x,y)＞Th时，令Gh(x,y)为1；当G(x,y)＜Th时，令Gh(x,y)为0；

用Sobel算子和八方向Sobel算子对地震数据做边缘信息检测，得到的同相轴边缘结果如图2所示。

图2 地震数据的边缘信息检测Fig.2 Edge information detection of seismic data

对比图2b）和图2c），可以很明显看出对于具有复杂的地震数据，后者可以得到较前者更清晰、更接近原始地震数据的同相轴边缘。

2 非局部均值算法算法降噪

对于含噪地震数据v(i)，其观测模型为

式中：u(i)是不含噪声的原始地震数据；n(i)是均值为0，方差为δ2的高斯白噪声。降噪后得到的原始地震数据的估计值u^(i)、u^(i)越接近原始地震数据u(i)，降噪效果越好。

设Ni是以i点为中心的邻域，Nj是以j点为中心的邻域，滤波前i点的灰度值为y(i)，滤波后i点的灰度估计值y^(i)可表示为

式中：w(i,j)表示在搜索邻域中Ni和Nj之间的权值，0≤w(i,j)≤1，∑j∈Iw(i,j)=1权对比示例如图3所示。

图3 地震数据中的自相似性块Fig.3 Selfsimilarity blocks in seismic data

以点q1和q2为中心的邻域和以点p为中心的邻域具有相似性，而以点p为中心的邻域和以点q3为中心的邻域与他们相差较大，在非局部均值算法中，q1和q2被赋予较大的权值，但q3被赋予较小的权值，这样就可以滤除与目标点不相似的点，有利于目标点数值的估计。

高斯加权欧式距离定义为

相似点之间的权重计算表达式为

式中：h为滤波参数，通过控制指数函数的衰减速度，从而控制平滑程度；D(i,j)为高斯加权欧式距离；Ga是标准差为a的高斯核函数，用于确定Ni与Nj之间欧氏距离，距离较远的两点权值较小，距离较近的两点权值较大。

3 结合边缘检测和非局部均值的地震数据降噪

在数字化后的图像中，对于图像中的每个点，都有与之对应的8个邻域，本文选用八方向Sobel边缘检测得到的结果改进权值，定义如下：

式中：ed(Ni)表示边缘图像中以i为中心的邻域；ed(Nj)表示边缘图像中以j为中心的邻域；‖ed(Ni)-ed(Nj)‖表示利用标准差为b的高斯核计算得到的i和j之间的边缘距离，计算方法与计算高斯加权欧式距离相似，并且a=b。

由权值函数可得，在地震数据的平滑区域，各点的边缘距离近似为0，对权值的影响不大，而在地震数据的同相轴边缘，距离变大，此时权值变小。在地震数据的估计中，相关度高的邻域权值变化不大，相关度低的邻域权值减小，这样有利于在降噪的同时较完整恢复数据。

综上，利用结合边缘检测和非局部均值的地震数据降噪步骤如下：

1）对大小为(M,N)的含噪地震数据做高斯平滑预处理；

2）分别计算每一个点与八方向模板的卷积，得到8个结果，取最大的卷积结果作为地震数据的梯度值；

3）设置精度δ为150，通过迭代得到最佳阈值Th，对于地震数据的梯度幅值矩阵的每一个点进行阈值分割，得到同相轴边缘；

4）取邻域窗口边长为( 2f+1)，搜索窗口边长为( 2t+1)，对地震数据进行对称扩充，扩充后数据大小为(M+f+t)×(N+f+t)，计算每个点的邻域与搜索窗口中所有点的高斯加权欧式距离；

5）根据欧式距离计算其归一化常数，用改进后的权值函数计算目标点与搜索窗口中所有像素点之间的权值；

6）通过加权平均得到目标点的新灰度值；

7）重复步骤3～5，遍历含噪地震数据的每个采样点。

4 地震数据降噪实验结果与分析

实验对象为合成地震数据、实际地震数据，实验所用算法为NLM、Sobel2-NLM、Canny-NLM和Sobel8-NLM，所有算法参数相同。实验在一台CPU主频为2.40 GHz、内存为8 GB、安装有Microsoft Windows 10家庭中文版、64位操作系统的个人笔记本电脑上进行，运行环境为MATLAB（R2014b）。

本文使用峰值信噪比（PSNR）、均方误差（MSE）和平均结构相似性（MSSIM）3个指标评价地震数据的降噪效果。PSNR和MSSIM越高，MSE越低，表示算法说明残留噪声在降噪数据中的比重越小，算法的降噪能力越好。

4.1 合成地震数据降噪

合成地震数据为已经过预处理的地震数据，共计90道，每道160个采样点，给其加入均值为0噪声水平分别为信号幅值15%的高斯随机噪声。用NLM、Sobel2-NLM、Canny-NLM和Sobel8-NLM对其降噪，结果如图4所示。

图4 合成地震数据降噪结果对比Fig.4 Comparison of noise reduction results of synthetic seismic data

从图4可以看到，Canny-NLM和Sobel8-NLM的降噪结果优于NLM和Sobel2-NLM，但Canny-NLM对噪声过于敏感，同相轴边缘附近在降噪后出现较多不规则斑块，背景区域不够平滑，Sobel8-NLM算法降噪效果优于以上3种算法，地震剖面同相轴较为清晰，背景区域较平滑。

为定量分析各算法的降噪效果，现给该信号添加均值为0、噪声水平分别为原信号幅值5%、10%、15%、20%的高斯白噪声，并将各降噪算法得到的PSNR、MSE、MSSIM指标列于表1中。

表1 合成地震数据的不同降噪算法性能指标Tab.1 Performance indicators of different noise reduction algorithms for synthetic seismic data

随着噪声水平增加，4种算法降噪性能均有所下降，在噪声水平达到20%时，Sobel2-NLM、Canny-NLM、Sobel8-NLM的PSNR比NLM算法高，且MSSIM一直高于0.999 0，说明在地震数据降噪中，结合边缘检测的非局部均值算法降噪效果好，和原始数据相似度高。其中，Sobel8-NLM的PSNR明显高于其它算法且MSE最低，验证了本文算法的有效性。

4.2 实际地震数据降噪实验

对实际地震数据进行降噪实验，共计600道，单道含240个采样点，给其加入均值为0、噪声水平分别为信号幅值10%的高斯随机噪声，分别用Sobel2-NLM、Canny-NLM和Sobel8-NLM对其降噪，抽取实际地震数据第90道和其他3种算法降噪后第90道数据，结果对比如图5所示。

图5 实际地震数据降噪（第90道）Fig.5 Noise reduction of actual seismic data（the 90th channel）

从图5可以看出，对于单道数据而言，Sobel2-NLM降噪后振幅被平滑，Canny-NLM也出现了小部分振幅压缩或消失，Sobel8-NLM降噪后结果与原数据更加接近，能够更加完整保留地震数据的信息。设置均值为0，噪声水平分别为原信号幅值5%、10%、15%、20%的加性高斯白噪声，并将各降噪算法得到的PSNR、MSE、MSSIM等性能指标展示于表2中。

表2 实际地震数据的不同降噪算法性能指标Tab.2 Performance indicators of different noise reduction algorithms for field seismic data

在实际地震的降噪实验中，Sobel8-NLM算法的PSNR最高，MSE较其他算法较低，噪声残留更少。噪声水平大于25%时，Canny-NLM、Sobel2-NLM的性能逐渐接近NLM，但Sobel8-NLM依然可以保持较高的PSNR和MSSIM值，验证了本文算法对实际地震数据降噪效果更好。

5 结论

本文利用八方向Sobel算子检测同相轴边缘，并改进非局部均值算法的权值函数，使得相关度低的邻域之间权值减小，实现有效降噪。在合成地震数据降噪实验中，验证了结合边缘检测的非局部均值算法具有更高的峰值信噪比及更低的均方误差，与原始地震数据相似度高。在实际地震数据降噪实验中，验证了结合八方向Sobel算子的非局部均值算法，对地震数据具有更高的降噪性能。对于地震数据中其他类型噪声的降噪算法仍需继续研究。