崔亞彤，王勝侯，蔡忠賢

(1.天津市勘察設(shè)計(jì)院集團(tuán)有限公司，天津 300191；2.中國地質(zhì)大學(xué)(武漢) 資源學(xué)院，湖北 武漢 430074)

0 引言

噪聲壓制在地震數(shù)據(jù)處理中起著至關(guān)重要的作用，去噪效果將直接影響地震資料解釋工作。為了有效壓制地震資料隨機(jī)噪聲，提高地震資料的信噪比，國內(nèi)外許多學(xué)者提出了不同的方法，本文將這些方法大致分為5類。首先是基于數(shù)學(xué)變換的方法，如Fourier變換[1]、Radon變換[2-3]、Curvelet變換[4-6]、Dreamlet變換[7]和Wavelet變換[8-9]；第二種是基于預(yù)測(cè)濾波的方法，如F-X預(yù)測(cè)濾波[10-11]和自適應(yīng)預(yù)測(cè)方法[12]；第三種是基于非局部均值(NLM)濾波的方法，如傳統(tǒng)NLM方法[13-14]和自適應(yīng)NLM濾波[15-16]；第四種是基于降秩的方法，如矩陣降秩法[17-19]和張量降秩法[20-21]；第五種是基于人工智能的方法，如數(shù)據(jù)驅(qū)動(dòng)緊密框架[22]和機(jī)器學(xué)習(xí)[23-27]。

與預(yù)測(cè)濾波方法和降秩濾波方法不同的是，NLM方法并不基于線性假設(shè)，因此在處理彎曲同相軸時(shí)，該方法可以有效地保護(hù)有效信號(hào)，壓制隨機(jī)噪聲。傳統(tǒng)NLM方法起源于圖像的隨機(jī)噪聲壓制處理[28]，隨后Bonar成功地將該方法引入到地震數(shù)據(jù)噪聲壓制處理中[13]。然而，傳統(tǒng)NLM方法在應(yīng)用上也存在一定局限性。相比于矩陣降秩或預(yù)測(cè)濾波等方法，該方法計(jì)算時(shí)間較長(zhǎng)，在處理大型地震數(shù)據(jù)時(shí)效率低下。為了解決這一問題，前人提出了分塊NLM法[28]、并行分塊NLM法[29]、基于隨機(jī)投影算法的NLM法[30]、下采樣的快速NLM法[31]、變窗口的快速NLM法[32]等。但由于前述方法為了提高計(jì)算效率，在計(jì)算過程中并沒有完全遍歷每個(gè)數(shù)據(jù)點(diǎn)，因此可能會(huì)犧牲計(jì)算精度。而基于數(shù)據(jù)積分算法的快速NLM法[33]，原理上等價(jià)于在計(jì)算過程中遍歷所有數(shù)據(jù)點(diǎn)，因此在提高計(jì)算速度的同時(shí)避免了犧牲精度的可能。傳統(tǒng)NLM方法在處理實(shí)際地震數(shù)據(jù)時(shí)，所選擇的濾波參數(shù)值通常為一常數(shù)，為了進(jìn)一步提高去噪效果，前人利用結(jié)構(gòu)張量算法[34-35]、矩陣本征特性算法[36]、灰色關(guān)聯(lián)分析算法[37]等方法，針對(duì)不同區(qū)域可選擇不同的濾波參數(shù)來提高去噪效果，但會(huì)明顯增加計(jì)算量。NLM方法的去噪效果很大程度上受濾波器參數(shù)的選擇影響，如果參數(shù)太大，就會(huì)丟失細(xì)節(jié)信息，使圖像模糊，如果參數(shù)過小，則噪聲不能被完全抑制。因此，一些學(xué)者通過局部數(shù)據(jù)的隨機(jī)噪聲估計(jì)來自適應(yīng)的計(jì)算濾波參數(shù)，例如高通濾波法[38]和最小方差估計(jì)法[15-16]。上述這些方法在提高傳統(tǒng)NLM方法去噪效果的同時(shí)，可以減少濾波參數(shù)對(duì)計(jì)算速度的影響，但該方法在計(jì)算效率上仍有提升空間。

本文提出了一種基于快速自適應(yīng)非局部均值濾波的地震隨機(jī)噪聲壓制方法，可快速且有效地壓制地震隨機(jī)噪聲。首先，為了提高NLM方法的計(jì)算效率，本文給出了一種中心對(duì)稱數(shù)據(jù)積分算法，有效地降低了計(jì)算成本。其次，通過計(jì)算兩個(gè)鄰域窗口的相似度，本文給出了利用均勻度來自適應(yīng)調(diào)整濾波參數(shù)分布的方法。最后，通過模型數(shù)據(jù)和實(shí)際數(shù)據(jù)驗(yàn)證了該方法的有效性、實(shí)用性。

1 方法

1.1 非局部均值濾波方法

噪聲數(shù)據(jù)Dnoise(t,x)可由下式表示：

Dnoise(t,x)=Dtrue(t,x)+n

(1)

式中：Dtrue(t,x)表示無噪聲數(shù)據(jù)；n表示隨機(jī)噪聲。假設(shè)(Ds,Ds)表示搜索窗口半徑，(ds,ds)表示鄰域窗口半徑，Ddenoise(t,x)表示去噪后的數(shù)據(jù)，通過對(duì)Dtrue(t,x)進(jìn)行加權(quán)平均計(jì)算即可得到每個(gè)點(diǎn)處的Ddenoise(t,x)數(shù)據(jù)[13, 28]。因此，Ddenoise(t,x)可以表示為：

(3)

Ω={i,j∈Ω:|i-t|≤Ds,|j-x|≤Ds}。

(4)

NLM方法可以有效地處理彎曲同相軸，但計(jì)算成本較高。假設(shè)Dnoise(t,x)總點(diǎn)數(shù)為N=NtNx，兩個(gè)正方形鄰域窗口之間的相似度計(jì)算量為d2=(2ds+1)2。對(duì)于Dnoise(t,x)中的每個(gè)數(shù)據(jù)點(diǎn)，必須在搜索窗內(nèi)計(jì)算D2=(2Ds+1)2次相似度，NLM濾波的計(jì)算復(fù)雜度為O(ND2d2)，計(jì)算量巨大。因此，傳統(tǒng)NLM方法在處理大型實(shí)際地震數(shù)據(jù)時(shí)，會(huì)明顯受到計(jì)算效率的制約。

1.2 中心對(duì)稱數(shù)據(jù)積分算法

={r1,r2∈；|r1|≤Ds,|r2|≤Ds} ，

(5)

={t,x∈;1≤t≤t1,1≤x≤x1} ，

(6)

[St(t+ds,x+ds)+St(t-ds-1,x-ds-1)-

St(t+ds,x-ds-1)-St(t-ds-1,x+ds)]，

(7)

通過觀察式(5)，發(fā)現(xiàn)建立的差分矩陣s(t,x)為中心對(duì)稱的，因此本文給出了一種中心對(duì)稱數(shù)據(jù)積分算法，利用數(shù)據(jù)積分的數(shù)學(xué)對(duì)稱性來提高NLM方法的計(jì)算效率。假設(shè)

進(jìn)一步的，上式可推導(dǎo)為:

(9)

據(jù)此s(t,x)[-r]的值可被表示為s(t,x)[+r]的中心對(duì)稱值，在計(jì)算過程中，即可避免構(gòu)建s(t,x)[-r]矩陣，進(jìn)而減少了一半的計(jì)算量，將式(9)代入式(6)中，

St(t1-r1,x1-r2)[+r]=St(t1,x1)[-r]，

(10)

再將式(10)代入式(7)中，得到：

St(t+ds-r1,x-ds-r2-1)[+r]-St(t-ds-r1-1,x+ds-r2)[+r]]。

(11)

1.3 均勻性估計(jì)

式(2)中的濾波參數(shù)h的選擇，對(duì)于NLM方法的去噪效果是至關(guān)重要的。大量地震資料表明，不同地層的有效信號(hào)能量差異較大，甚至隨機(jī)噪聲的分布也不是完全隨機(jī)的，因此濾波參數(shù)h的取值通常難以確定。相比于整個(gè)區(qū)域的地震數(shù)據(jù)取相同的濾波參數(shù)h， 利用一個(gè)參數(shù)矩陣h來控制不同區(qū)域數(shù)據(jù)的去噪水平可以提高噪聲壓制效果。基于最小方差估計(jì)的自適應(yīng)算法可以自適應(yīng)選擇濾波參數(shù)h[15-16]，該方法認(rèn)為參數(shù)h的估計(jì)是噪聲σ的標(biāo)準(zhǔn)差，滿足下式

2σ2，

(12)

其中，V0代表無噪聲數(shù)據(jù)。自適應(yīng)濾波參數(shù)矩陣h2的估計(jì)即可表示為

(13)

通過利用最小方差估計(jì)的方法即可調(diào)整整個(gè)數(shù)據(jù)區(qū)域內(nèi)部的濾波參數(shù)值。圖1給出了一個(gè)含噪聲數(shù)據(jù)估計(jì)的濾波參數(shù)h2的分布情況，如圖1c所示，利用最小方差估計(jì)方法得到的濾波參數(shù)h2在有效信號(hào)結(jié)構(gòu)邊緣處相對(duì)較大，而背景均勻區(qū)域的濾波參數(shù)值相對(duì)較小。

a—無噪聲的數(shù)據(jù)；b—加噪數(shù)據(jù)；c—最小方差估計(jì)的參數(shù)分布情況；d—相似度標(biāo)準(zhǔn)差的參數(shù)分布情況a—noise-free data;b—noisy data;c—parameter distribution of minimum variance estimation;d—parameter distribution of standard deviation of similarity圖1 濾波參數(shù)h2分布Fig.1 Distribution of the filtering parameter h2

(14)

式中，STDs(t,x)代表相似度的均方差，

Ω={i,j∈Ω；|i-t|≤Ds,|j-x|≤Ds} 。

(15)

本文給出的基于快速自適應(yīng)非局部均值濾波方法，通過引入中心對(duì)稱數(shù)據(jù)積分算法進(jìn)一步提高了傳統(tǒng)NLM方法的計(jì)算效率，同時(shí)利用相似度標(biāo)準(zhǔn)差來估計(jì)均勻性，實(shí)現(xiàn)了自適應(yīng)濾波參數(shù)調(diào)整，有效地提高了去噪效果和計(jì)算效率。

2 模型數(shù)據(jù)試驗(yàn)

在地震數(shù)據(jù)噪聲壓制處理中，通常采用信噪比(SNR)[40]、峰值信噪比(PSNR)和均方誤差(MSE)[41]來定量分析方法的地震數(shù)據(jù)去噪效果，其定義如下:

(16)

(17)

(18)

式中：Dtrue和Ddenoise分別代表原始無噪聲數(shù)據(jù)和噪聲壓制后數(shù)據(jù)；N=NtNx表示噪聲數(shù)據(jù)Dnoise(t,x)的總點(diǎn)數(shù)。結(jié)合傳統(tǒng)NLM方法[13]和基于最小方差估計(jì)的NLM方法[15]與本文快速自適應(yīng)NLM方法進(jìn)行對(duì)比，利用2個(gè)模型地震數(shù)據(jù)，從計(jì)算效率和精度兩方面分別驗(yàn)證本文方法的有效性。

模型數(shù)據(jù)試驗(yàn)1為簡(jiǎn)單合成地震剖面，包括兩個(gè)線性同相軸和兩個(gè)彎曲同相軸。圖2a為60道無噪聲模型數(shù)據(jù)，圖2b為加噪數(shù)據(jù)，信噪比為-4.55 dB，圖2c～h分別為傳統(tǒng)NLM法、基于最小方差估計(jì)的NLM法和本文方法的去噪結(jié)果及相應(yīng)的殘差剖面。由于F-X預(yù)測(cè)濾波方法[10]和降秩方法[18]算法本身基于線性假設(shè)，因此對(duì)于彎曲同相軸的處理存在誤差，而NLM方法對(duì)于彎曲同相軸的計(jì)算精度就要明顯高于上述兩種方法(圖2c～h)。由于從圖中很難看出3種NLM方法的不同之處，因此分別計(jì)算上述方法去噪結(jié)果的SNR、PSNR和MSE，如表1?？梢钥闯觯帽疚姆椒ǖ玫降腟NR和PSNR都高于其他方法，且MSE明顯低于其他方法，該方法對(duì)于隨機(jī)噪聲的壓制效果要優(yōu)于上述其他方法。

a—無噪聲的數(shù)據(jù)；b—噪聲數(shù)據(jù)(SNR=-4.55 dB)；c—傳統(tǒng)NLM方法去噪結(jié)果；d—a與c的差剖面；e—基于最小方差估計(jì)的NLM方法去噪結(jié)果；f—a與e的差剖面；g—本文方法去噪的結(jié)果；h—a與g的差剖面a—noise-free data;b—noisy data(SNR=-4.55 dB);c—denoised result by using conventional NLM method;d—difference between a and c;e—denoised result by using NLM based on minimum variance estimation;f—difference between a and e;g—denoised result by using proposed method;h—difference between a and g圖2 模型試驗(yàn)1的合成模型數(shù)據(jù)噪聲壓制結(jié)果Fig.2 Noise attenuation results of synthetic model data from model data test 1

表1 不同方法的SNR、PSNR和MSE對(duì)比Table 1 Comparison of SNR,PSNR and MSE using different methods

模型數(shù)據(jù)試驗(yàn)2為正演偏移剖面。圖3a為100道無噪聲模型數(shù)據(jù)，圖3b為噪聲數(shù)據(jù)，其信噪比為-3.01 dB，圖3c～h分別展示了傳統(tǒng)NLM方法、基于最小方差估計(jì)的NLM方法和本文方法去噪后的結(jié)果及相應(yīng)的殘差剖面。上述方法去噪后的SNR分別為17.254 5 dB、20.727 4 dB、22.557 9 dB，利用本文方法去噪結(jié)果的SNR最高。結(jié)果表明，本文方法除在合成數(shù)據(jù)的殘差剖面中泄漏少量能量外，在壓制隨機(jī)噪聲和保留有效信號(hào)方面具有較好的效果。

a—無噪聲的數(shù)據(jù)；b—噪聲數(shù)據(jù)(SNR=-3.01 dB)；c—傳統(tǒng)NLM方法去噪結(jié)果；d—a與c的差剖面；e—基于最小方差估計(jì)的NLM方法去噪結(jié)果；f—a與e的差剖面；g—本文方法去噪的結(jié)果；h—a與g的差剖面a—noise-free data;b—noisy data (SNR=-3.01 dB);c—denoised result by using conventional NLM method;d—difference between a and c;e—denoised result by using NLM based on minimum variance estimation;f—difference between a and e;g—denoised result by using proposed method;h—difference between a and g圖3 模型試驗(yàn)2的合成模型數(shù)據(jù)噪聲壓制結(jié)果Fig.3 Noise attenuation results of synthetic model data from model data test 2

最后，利用一系列不同大小的地震數(shù)據(jù)對(duì)比不同方法的計(jì)算時(shí)間。圖4a顯示了本文快速自適應(yīng)NLM方法、基于數(shù)據(jù)積分方法的NLM方法、基于最小方差估計(jì)的分塊NLM方法和傳統(tǒng)NLM方法的計(jì)算時(shí)間對(duì)比。圖4b顯示了不同信噪比的噪聲地震模型數(shù)據(jù)下，這些方法的去噪質(zhì)量。與其他方法相比，本文方法計(jì)算效率明顯高于其他方法，且當(dāng)利用不同信噪比的模型進(jìn)行去噪處理后，本文方法也具有更好的去噪質(zhì)量。

a—計(jì)算時(shí)間對(duì)比；b—去噪質(zhì)量對(duì)比a—computational time comparison;b—denoising quality comparison圖4 基于模型數(shù)據(jù)2不同方法去噪計(jì)算時(shí)間及效果對(duì)比Fig.4 Comparison of calculation time and effect of different denoising methods based on model data test 2

3 實(shí)際數(shù)據(jù)試驗(yàn)

為了驗(yàn)證本文快速自適應(yīng)非局部均值濾波方法的實(shí)用性，將該方法應(yīng)用于實(shí)際地震數(shù)據(jù)的隨機(jī)噪聲壓制處理中，結(jié)果如圖5所示。從圖5a中方框區(qū)域可以看出，隨機(jī)噪聲污染了實(shí)際地震資料，使得地層模糊，斷層構(gòu)造不清晰，難以進(jìn)行后續(xù)地震資料解釋工作。圖5b～d分別展示了傳統(tǒng)NLM方法、基于最小方差估計(jì)的NLM方法以及本文方法的去噪結(jié)果。圖5e～g顯示了這些方法的差剖面。從圖5b～g可以看出，傳統(tǒng)的NLM方法、基于最小方差估計(jì)的NLM方法以及本文提出的方法都可以有效壓制隨機(jī)噪聲。但從圖5b、5c和5d中的方框區(qū)域可以看出，本文方法處理結(jié)果中的斷層結(jié)構(gòu)相比于另外兩種方法更加清晰，有效信號(hào)的信息得到了很好的保存，有效地壓制了隨機(jī)噪聲。在實(shí)際地震數(shù)據(jù)的噪聲壓制處理中，本文快速自適應(yīng)NLM方法、基于數(shù)據(jù)積分方法的NLM方法、基于最小方差估計(jì)的分塊NLM方法和傳統(tǒng)NLM方法的計(jì)算時(shí)間分別為：0.63 s、1.14 s、58.42 s和94.22 s。本文提出方法的計(jì)算時(shí)間也有明顯的提高。因此，在處理大型地震數(shù)據(jù)資料時(shí)，本文快速自適應(yīng)非局部均值濾波方法具有更好的實(shí)用性。

a—實(shí)際地震數(shù)據(jù);b—傳統(tǒng)NLM方法去噪結(jié)果;c—基于最小方差估計(jì)的NLM去噪結(jié)果;d—本文方法去噪結(jié)果;e—a與b的差剖面;f—a與c的差剖面;g—a與e的差剖面a—field data;b—denoised result by using conventional NLM method;c—denoised result by using NLM based on minimum variance estimation;d—denoised result by using proposed method;e—difference between a and b;f—difference between a and c;g—difference between a and e圖5 實(shí)際地震數(shù)據(jù)去噪結(jié)果Fig.5 Noise attenuation results of field data

4 結(jié)論

本文給出了一種用于地震資料隨機(jī)噪聲壓制的快速自適應(yīng)非局部均值濾波方法。首先，基于數(shù)據(jù)積分算法的數(shù)學(xué)對(duì)稱性，利用中心對(duì)稱數(shù)據(jù)積分算法來加速傳統(tǒng)NLM濾波方法。其次，利用相似度標(biāo)準(zhǔn)差來估計(jì)均勻性，自適應(yīng)地計(jì)算NLM濾波參數(shù)，進(jìn)一步提高噪聲壓制效果。因此，本文方法在有效提高計(jì)算效率的同時(shí)，又提高了噪聲壓制的效果。最后，通過對(duì)模型數(shù)據(jù)和實(shí)際數(shù)據(jù)的處理，驗(yàn)證了該方法的有效性、實(shí)用性。