結合邊緣檢測和非局部均值的地震數據降噪

李曉璐，周亞同，何靜飛，翁麗源，李書華

（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算子的非局部均值算法，對地震數據具有更高的降噪性能。對于地震數據中其他類型噪聲的降噪算法仍需繼續研究。