王 偉，高靜懷，陳文超＊，朱振宇

1 西安交通大學(xué)電子與信息工程學(xué)院波動(dòng)與信息研究所，西安 710049

2 中海石油研究總院，北京 100027

基于結(jié)構(gòu)自適應(yīng)中值濾波器的隨機(jī)噪聲衰減方法

王 偉1，高靜懷1，陳文超1＊，朱振宇2

1 西安交通大學(xué)電子與信息工程學(xué)院波動(dòng)與信息研究所，西安 710049

2 中海石油研究總院，北京 100027

本文提出一種保護(hù)斷層、裂縫等地層邊緣特征的結(jié)構(gòu)自適應(yīng)中值濾波器，用于衰減地震資料中的隨機(jī)噪聲.基于地震反射同相軸局部呈線(xiàn)型結(jié)構(gòu)的假設(shè)，采用梯度結(jié)構(gòu)張量估計(jì)地層傾向，分析地層結(jié)構(gòu)的規(guī)則程度，在此基礎(chǔ)上引入地震剖面中線(xiàn)型和橫向不連續(xù)性?xún)煞N結(jié)構(gòu)特征的置信度量.結(jié)構(gòu)自適應(yīng)中值濾波器根據(jù)這兩種置信度量調(diào)整濾波器窗函數(shù)的尺度和形狀，根據(jù)地層傾角調(diào)整濾波器窗函數(shù)的方向，從而使得濾波操作窗能夠最佳匹配信號(hào)的局部結(jié)構(gòu)特征.將本文方法用于合成和實(shí)際數(shù)據(jù)的處理，并與兩種常用中值濾波方法進(jìn)行對(duì)比，結(jié)果表明，該方法能夠更好地解決地震剖面的隨機(jī)噪聲衰減和有效信號(hào)保真的問(wèn)題，在增強(qiáng)反射同相軸的橫向一致性的同時(shí)有效保持了剖面內(nèi)的地層邊緣和細(xì)節(jié)特征，顯著改善了地震資料的品質(zhì).

中值濾波，自適應(yīng)濾波，噪聲衰減，構(gòu)造保護(hù)，梯度結(jié)構(gòu)張量

1 引 言

在石油天然氣勘探開(kāi)發(fā)中，斷層、通道和裂縫等結(jié)構(gòu)特征以及河道砂體等沉積特征是發(fā)現(xiàn)和描述油藏的基礎(chǔ)，分辨和分析地震資料中這些結(jié)構(gòu)特征具有重要的意義.然而，由于這些區(qū)域地質(zhì)形態(tài)的復(fù)雜性，使得地震數(shù)據(jù)中這些區(qū)域易受噪聲的干擾，無(wú)論采用人工解釋還是相干和邊緣檢測(cè)等自動(dòng)檢測(cè)算法來(lái)描述和刻畫(huà)這些區(qū)域都變得較為困難.傳統(tǒng)的地震資料去噪方法，如F－K域傾向?yàn)V波、f－x域預(yù)測(cè)濾波和多道相干濾波等，多是利用地震反射同相軸的橫向一致性，易造成對(duì)傾斜同相軸振幅的壓制，也會(huì)模糊小的斷層和裂縫等不連續(xù)結(jié)構(gòu)，甚至?xí)鸫蟮臄鄬觾蓚?cè)同相軸的錯(cuò)誤連接.如何能夠較好地平衡地震噪聲衰減和有效信號(hào)保護(hù)的關(guān)系受到了廣泛的關(guān)注［1－6］，而能夠保護(hù)和增強(qiáng)斷裂等邊緣結(jié)構(gòu)的疊后噪聲衰減技術(shù)尤其受到地震解釋工作者的重視［7－9］.

Luo等［10］基于次序統(tǒng)計(jì)的思想給出一種地震數(shù)據(jù)保邊濾波（EPS，Edge－Preserving Smoothing）算法，該方法采用Kuwahara多窗分析技術(shù)，計(jì)算當(dāng)前分析點(diǎn)周?chē)鱾€(gè)子窗內(nèi)數(shù)據(jù)的均值和方差，將最小方差對(duì)應(yīng)的子窗的均值作為當(dāng)前點(diǎn)的濾波輸出.AlBinHassan等［11］發(fā)展了這種保邊濾波思想，給出了三維情況下的子窗剖分方法，將EPS算法用于三維地震數(shù)據(jù)的保邊濾波處理.Lu等［12］基于多項(xiàng)式擬合的思想推廣了一維的EPS算法，該方法通過(guò)計(jì)算包含當(dāng)前分析點(diǎn)的一系列平移子窗內(nèi)信號(hào)對(duì)給定階次的多項(xiàng)式擬合誤差，選擇出最小擬合誤差對(duì)應(yīng)的子窗的信號(hào)估計(jì)值作為當(dāng)前點(diǎn)的濾波輸出.楊培杰等［13］也基于這種多窗分析思想給出一種方向性邊界保持?jǐn)鄬釉鰪?qiáng)技術(shù)，使得相干體中的斷層圖像得到更加清晰顯示.Liu等［14－15］應(yīng)用圖像處理中的二維多級(jí)中值濾波器消除地震隨機(jī)噪聲，并分析了濾波器尺度參數(shù)選擇對(duì)噪聲衰減和信號(hào)細(xì)節(jié)保護(hù)的影響.此外，F(xiàn)ehmers等［16］引入微分方程濾波法用于地震資料的保邊濾波處理，該方法采用梯度結(jié)構(gòu)張量估計(jì)地層方向，度量地層的連續(xù)性，以此約束擴(kuò)散濾波方程的濾波方向和濾波程度.Lavialle等［17］結(jié)合線(xiàn)型和面型擴(kuò)散濾波模型給出專(zhuān)門(mén)針對(duì)保護(hù)和增強(qiáng)斷層構(gòu)造的斷層保持?jǐn)U散濾波器.孫夕平等［18］和王旭松等［19］通過(guò)在擴(kuò)散方程中引入二階導(dǎo)數(shù)項(xiàng)、張爾華等［20］通過(guò)修正擴(kuò)散方程中的保邊約束項(xiàng)進(jìn)一步改善非線(xiàn)性各向異性擴(kuò)散濾波器的邊緣保護(hù)性能.然而，基于各向異性擴(kuò)散方程的濾波方法以原始地震數(shù)據(jù)作為初始條件，通過(guò)擴(kuò)散時(shí)間的迭代推移形成濾波作用，由于不好估計(jì)出獲得最佳濾波效果的擴(kuò)散時(shí)間，導(dǎo)致不能準(zhǔn)確預(yù)測(cè)其濾波性態(tài)［21］.

經(jīng)典的二維中值濾波法具有優(yōu)異的噪聲衰減性能，尤其適用于消除非平穩(wěn)信號(hào)中的峰值噪聲.但是，由于在濾波過(guò)程中其濾波窗口保持恒定，且濾波操作沒(méi)有方向特性，使得采用中值濾波法處理具有層狀結(jié)構(gòu)的地震信號(hào)時(shí)，必然會(huì)引起有用信號(hào)的損傷和不連續(xù)結(jié)構(gòu)特征的模糊.一些改進(jìn)的中值濾波方法如二維多級(jí)中值濾波器［14－15］，雖然有較強(qiáng)的信號(hào)細(xì)節(jié)保護(hù)能力，但是其對(duì)噪聲的衰減能力不夠.本文結(jié)合微分方程濾波法中利用梯度結(jié)構(gòu)張量刻畫(huà)地震反射層局部結(jié)構(gòu)特征的方法，定義了帶有能量變化強(qiáng)度信息的橫向不連續(xù)性置信度量，結(jié)合線(xiàn)型置信度量和橫向不連續(xù)性置信度量，給出一種二維結(jié)構(gòu)自適應(yīng)中值濾波方法，使中值濾波操作窗函數(shù)根據(jù)信號(hào)局部結(jié)構(gòu)特征自適應(yīng)地作方向延伸和尺度伸縮，從而既能最大限度地消除隨機(jī)噪聲，又能有效保護(hù)有用信號(hào)和地層邊緣等細(xì)節(jié)特征.

2 結(jié)構(gòu)自適應(yīng)中值濾波原理

2.1 地層傾向估計(jì)及橫向不連續(xù)性度量

有多種方法可用于地層結(jié)構(gòu)的規(guī)則性分析［16，22－25］，梯度結(jié)構(gòu)張量法是其中一種簡(jiǎn)單而又有效的途徑.記二維地震剖面為u（x，t），則該二維剖面對(duì)應(yīng)的梯度結(jié)構(gòu)張量由下式描述

其中，和T分別為梯度算子和轉(zhuǎn)置算子，?表示卷積運(yùn)算，Gρ表示尺度為ρ的二維高斯函數(shù)Gρ（x，t）＝exp（－（x2＋t2）／2ρ2）.在梯度結(jié)構(gòu)張量的定義中，由梯度向量Δu及其轉(zhuǎn)置向量形成的張量矩陣與尺度為ρ的高斯低通濾波器作用，確定了梯度結(jié)構(gòu)張量可度量的信號(hào)特征的最大尺度.將該對(duì)稱(chēng)半正定矩陣Sρ（Δu）作矩陣特征分解

其中，特征值μ1≥μ2≥0反應(yīng)了信號(hào)在位置（x，t）的高斯鄰域Gρ內(nèi)沿特征方向的振幅變化強(qiáng)度；特征向量v1和v2給出局部正交方位，v1對(duì)應(yīng)于局部振幅變化最大的方向，即信號(hào)梯度方向，而v2對(duì)應(yīng)于局部振幅變化最小的方向，即反射同相軸的取向.根據(jù)梯度結(jié)構(gòu)張量?jī)蓚€(gè)特征值的取值情況，可分辨不同的信號(hào)結(jié)構(gòu)特征，文獻(xiàn)［26］中引入線(xiàn)型信號(hào)結(jié)構(gòu)的置信度量：

該置信度量在區(qū)間［0，1］之間取值，衡量局部信號(hào)特征與線(xiàn)型結(jié)構(gòu)的相似程度，同時(shí)亦表征局部方位估計(jì)的準(zhǔn)確程度.若CL→1，則表示局部信號(hào)特征趨于線(xiàn)型結(jié)構(gòu)，估計(jì)的地層傾角的準(zhǔn)確度越高；反之，若CL→0，則局部信號(hào)結(jié)構(gòu)背離了線(xiàn)型假設(shè)，這可能由干擾噪聲和斷層、裂縫等橫向不連續(xù)性造成.根據(jù)線(xiàn)型結(jié)構(gòu)置信度量可較好地識(shí)別局部呈線(xiàn)型結(jié)構(gòu)的反射同相軸；但是，由于該置信度量值判斷的是信號(hào)能量的相對(duì)變化程度，容易受到噪聲的干擾，并且在無(wú)明顯信號(hào)結(jié)構(gòu)的區(qū)域內(nèi)取值沒(méi)有意義，從而很難從干擾背景中明確地識(shí)別斷層、裂縫等橫向不連續(xù)結(jié)構(gòu).

基于上述分析，本文在線(xiàn)型結(jié)構(gòu)置信度量的基礎(chǔ)上，給出一種含有信號(hào)橫向能量變化強(qiáng)度信息的置信度量值，用于識(shí)別斷層、裂縫等橫向不連續(xù)結(jié)構(gòu).

（4）式中（1－CL）項(xiàng)表示信號(hào)局部結(jié)構(gòu)特征相對(duì)線(xiàn)型結(jié)構(gòu)的背離程度，μ2表征在給定度量鄰域內(nèi)在均方誤差最小意義下信號(hào)沿局部一致方向的能量變化強(qiáng)度［24］.在反射同相軸呈線(xiàn)型結(jié)構(gòu)的區(qū)域，由于CL→1，即（1－CL）→0，此時(shí)沿地層傾向的信號(hào)變化相對(duì)較小，即μ2→0，從而使得CI→0.在含有斷層、裂縫等橫向不連續(xù)性的區(qū)域，由于CL→0，即（1－CL）→1；且在橫向不連續(xù)區(qū)域，由于信號(hào)結(jié)構(gòu)背離了局部線(xiàn)型結(jié)構(gòu)假設(shè)，此處信號(hào)的一致性取向不明確，因而沿局部一致性方向的信號(hào)變化強(qiáng)度也相對(duì)較大，即μ2?0，從而使得CI?0.因此，根據(jù)橫向不連續(xù)性置信度量CI可以很好地識(shí)別斷層和裂縫等橫向不連續(xù)結(jié)構(gòu).由于信號(hào)橫向能量變化項(xiàng)μ2的存在，對(duì)于噪聲能量弱于信號(hào)能量的情況，CI具有良好的抗噪能力；此外，CI也能夠判別無(wú)明顯信號(hào)結(jié)構(gòu)的區(qū)域（此時(shí)μ2→0，即CI→0），這對(duì)采用CI指導(dǎo)濾波過(guò)程尤為重要，可以避免一些濾波假象的產(chǎn)生.

2.2 結(jié)構(gòu)自適應(yīng)中值濾波器

對(duì)于二維地震剖面u（x），記x＝（x，t）為采樣點(diǎn)的空間、時(shí)間位置，定義位置x處的結(jié)構(gòu)自適應(yīng)中值濾波操作為：

式中，M（x，y）表示位置x處的結(jié)構(gòu)自適應(yīng)濾波窗口（見(jiàn)圖1a），y為當(dāng)前濾波窗口包含的采樣點(diǎn)的空間、時(shí)間位置，（x）表示位置x處的結(jié)構(gòu)自適應(yīng)中值濾波的輸出結(jié)果.結(jié)構(gòu)自適應(yīng)濾波窗M（x，y）是濾波位置x的函數(shù)，其大小和形狀取決于x處分析鄰域內(nèi)的信號(hào)結(jié)構(gòu)特征：

其中，·表示內(nèi)積運(yùn)算；n（x）和n⊥（x）分別對(duì)應(yīng)于x處的地層傾向及梯度方向，取為梯度結(jié)構(gòu)張量的兩個(gè)特征向量，即n（x）＝v2（x），n⊥（x）＝v1（x）；σ1（x）和σ2（x）分別表示橢圓形濾波窗口的長(zhǎng)軸和短軸，它們共同決定了結(jié)構(gòu)自適應(yīng)中值濾波器的濾波尺度以及濾波操作的方向選擇性.濾波器尺度參數(shù)的設(shè)計(jì)是影響結(jié)構(gòu)自適應(yīng)中值濾波器性態(tài)的關(guān)鍵，根據(jù)結(jié)構(gòu)自適應(yīng)濾波操作的要求，本文結(jié)合線(xiàn)型結(jié)構(gòu)置信度量CL和橫向不連續(xù)性置信度量CI，給出如下的濾波器尺度參數(shù)自適應(yīng)選擇策略：

式中Rmax規(guī)定了橢圓形濾波窗口的最大尺寸，g（·）為關(guān)于CI（x）的單調(diào)減函數(shù)，其取值范圍限定為（0，1］，文中取g（·）為指數(shù)函數(shù)

其中，β為針對(duì)CI（x）的閾值參數(shù)，控制指數(shù)函數(shù)的衰減速度，β取值越小，指數(shù)函數(shù)的下降越快；反之，指數(shù)函數(shù)的衰減變緩.

結(jié)構(gòu)自適應(yīng)中值濾波器根據(jù)信號(hào)的局部結(jié)構(gòu)特征具有靈活的濾波尺度調(diào)節(jié)功能和方向選擇特性，控制最大濾波尺度的σ1（x）由橫向不連續(xù)性置信度量CI（x）決定，而調(diào)整方向選擇特性的σ2（x）由線(xiàn)型結(jié)構(gòu)置信度量CL（x）決定.在斷層、裂縫等橫向不連續(xù)區(qū)域由于CI（x）?0，使得橢圓型窗函數(shù)的長(zhǎng)軸σ1（x）趨短；而在其他區(qū)域隨CI（x）→0，有σ1（x）→Rmax.在反射同相軸呈線(xiàn)型結(jié)構(gòu)區(qū)域，由于CL（x）→1，使得橢圓型窗函數(shù)的短軸σ2（x）萎縮；而在其他區(qū)域隨CL（x）→0，有σ2（x）→σ1（x）.因此，如圖1b所示，在反射同相軸局部呈線(xiàn)型延伸的區(qū)域，濾波窗口沿同相軸取向拉伸為狹長(zhǎng)的各向異性窗口（位置x1區(qū)域）；在斷層等橫向不連續(xù)區(qū)域，濾波窗口自適應(yīng)地萎縮為小窗以匹配此處的不連續(xù)結(jié)構(gòu)特征（位置x2區(qū)域）；而在無(wú)明顯信號(hào)結(jié)構(gòu)的噪聲干擾區(qū)域，濾波窗口延展為大的近似各向同性窗口（位置x3區(qū)域）.

圖1 結(jié)構(gòu)自適應(yīng)中值濾波器窗口的構(gòu)造（a）濾波器窗口函數(shù)；（b）信號(hào)結(jié)構(gòu)自適應(yīng)濾波策略.Fig.1 Filter window construction of structure－adaptive median filter（a）Filter window function；（b）Signal structure－consistent filtering mechanism.

聯(lián)合式（7）和（8）可見(jiàn)，結(jié)構(gòu)自適應(yīng)中值濾波器的濾波性能受到參數(shù)Rmax和β共同約束，最大尺寸參數(shù)Rmax控制濾波器的平滑性能，閾值參數(shù)β則調(diào)整濾波器對(duì)斷層等邊緣結(jié)構(gòu)的保持程度.考慮到文中采用梯度結(jié)構(gòu)張量估計(jì)地層的傾角和規(guī)則程度，梯度結(jié)構(gòu)張量的尺度參數(shù)ρ決定了其對(duì)局部信號(hào)結(jié)構(gòu)的度量孔徑，因此為了確?；趦煞N置信度量約束的濾波過(guò)程可信，要求濾波窗口的最大尺寸Rmax不超過(guò)2ρ，即取Rmax≤2ρ.對(duì)于閾值參數(shù)β，由圖2可見(jiàn)，當(dāng)β取值偏小時(shí)，稍大的橫向不連續(xù)性置信度量CI（x）即引起濾波窗口的尺寸減小，因而會(huì)使一些弱的斷層、裂縫等邊緣結(jié)構(gòu)得到有效保護(hù)，但同時(shí)也會(huì)減弱濾波器在橫向不連續(xù)性區(qū)域的濾波效果；當(dāng)β取值偏大時(shí)，只有很大的CI（x）才會(huì)導(dǎo)致濾波窗口萎縮到足夠小，因而在這種情況下濾波器的噪聲衰減能力得到增強(qiáng)，但在一定程度上會(huì)造成弱邊緣結(jié)構(gòu)的模糊和壓制.

圖2 結(jié)構(gòu)自適應(yīng)中值濾波窗的長(zhǎng)軸控制函數(shù)Fig.2 Controlling functions of the main axis of the structure－adaptive median filter kernel

2.3 保邊性態(tài)自適應(yīng)控制方法

結(jié)構(gòu)自適應(yīng)中值濾波器的邊緣結(jié)構(gòu)保持性能受閾值參數(shù)β控制，選擇合適的β是保證濾波過(guò)程中能夠有效保持?jǐn)鄬印⒘芽p等地層邊緣的關(guān)鍵.然而，地震資料通常從淺層到深層，即使在同一層段內(nèi)信號(hào)橫向能量變化的動(dòng)態(tài)范圍波動(dòng)也很大，很難統(tǒng)一地選擇一個(gè)全局閾值參數(shù)β.本文給出一種區(qū)域分割的閾值參數(shù)自適應(yīng)選擇策略，即根據(jù)待處理剖面內(nèi)橫向能量變化在時(shí)間和空間上的分布情況，將剖面Ω分割為若干具有一定邊界重疊的子區(qū)域Ω＝∪Ωi，在每個(gè)子區(qū)域Ωi內(nèi)取橫向不連續(xù)性置信度量CI（x）的極大值，即

從而通過(guò)下式得到該區(qū)域內(nèi)的自適應(yīng)閾值參數(shù)β為

式中α為根據(jù)保邊緣濾波要求而全局確定的百分比因子，thr是由剖面的整體噪聲干擾水平所確定的基底噪聲閾值.

通過(guò)將待處理剖面分割為若干子區(qū)域，由（9）和（10）式將（8）式中對(duì)閾值參數(shù)β的絕對(duì)取值轉(zhuǎn)化為對(duì)百分比因子α的相對(duì)取值，從而能夠統(tǒng)一調(diào)整結(jié)構(gòu)自適應(yīng)中值濾波器的邊緣保持性態(tài).在實(shí)際資料處理中，并不需要精確地確定子區(qū)域的大小，一種典型的選擇是每個(gè)子區(qū)域內(nèi)包含100道，每道包含150個(gè)采樣點(diǎn)，即可滿(mǎn)足濾波過(guò)程保持?jǐn)鄬?、裂縫等邊緣特征的要求.

3 合成模型分析

為了驗(yàn)證結(jié)構(gòu)自適應(yīng)中值濾波器的保護(hù)邊緣去噪的有效性，本文選取主頻為20Hz的Ricker子波，采用射線(xiàn)追蹤方法得到一個(gè)含有斷層的復(fù)雜地層結(jié)構(gòu)模型的合成地震數(shù)據(jù)（見(jiàn)圖3a），在合成數(shù)據(jù)中加入5dB的帶限高斯隨機(jī)噪聲，得到圖3b的含噪數(shù)據(jù).

圖3 含有斷層的復(fù)雜地層結(jié)構(gòu)的合成地震數(shù)據(jù)（a）原始數(shù)據(jù)；（b）加噪數(shù)據(jù).Fig.3 Synthetic seismic data of complex geology model with a fault（a）Clean data；（b）Noisy data.

采用結(jié)構(gòu)自適應(yīng)中值濾波器處理加噪的合成數(shù)據(jù)，考慮到合成數(shù)據(jù)的信噪比很低，在梯度結(jié)構(gòu)張量的計(jì)算中，二維高斯低通濾波函數(shù)取較大的尺度參數(shù)ρ＝4.0，結(jié)構(gòu)自適應(yīng)中值濾波器的最大濾波尺寸取為Rmax＝ρ＝4.0，即對(duì)應(yīng)最大尺度為9點(diǎn)的濾波窗口.針對(duì)圖3b中黑色長(zhǎng)方形框所選定的典型區(qū)域，分析閾值百分比因子α對(duì)濾波結(jié)果的影響.調(diào)整α從10%開(kāi)始以10%的間隔變化到200%，記錄每次濾波結(jié)果的信噪比，得到該區(qū)域的濾波結(jié)果的信噪比隨閾值百分比因子α的變化曲線(xiàn)，如圖4所示.從圖中可以看出，信噪比曲線(xiàn)隨α的增大而升高，且曲線(xiàn)在α＝90%后升高逐漸變緩.這是由于隨著α的增加，濾波器的保邊緣性能減弱，當(dāng)α＞90%后，濾波過(guò)程已造成斷面處波形的損傷，在一定程度上抵消了濾波性能增強(qiáng)帶來(lái)的信噪比增加.因而在處理該加噪的合成數(shù)據(jù)時(shí)，選擇閾值百分比因子α＝90%，此時(shí)對(duì)應(yīng)的結(jié)構(gòu)自適應(yīng)中值濾波器的各向異性窗函數(shù)的長(zhǎng)軸σ1和短軸σ2的取值分布如圖5所示.對(duì)比圖5a和圖5b可見(jiàn)，結(jié)構(gòu)自適應(yīng)中值濾波器的尺度參數(shù)能夠很好地匹配合成數(shù)據(jù)的信號(hào)結(jié)構(gòu)特征，在局部呈線(xiàn)型結(jié)構(gòu)的同相軸區(qū)域橢圓形濾波窗口的長(zhǎng)軸拉伸、短軸萎縮，形成有強(qiáng)方向選擇特性的各向異性窗口；在斷層區(qū)域長(zhǎng)軸和短軸同時(shí)萎縮，形成小尺度分析窗口；而在無(wú)信號(hào)特征的區(qū)域，長(zhǎng)軸和短軸均被拉伸，形成大尺度的各向同性窗口.

圖4 不同α對(duì)應(yīng)的結(jié)構(gòu)自適應(yīng)中值濾波結(jié)果的信噪比Fig.4 SNR of filtered results by structure－adaptive median filter using differentα

圖5 針對(duì)加噪合成數(shù)據(jù)的結(jié)構(gòu)自適應(yīng)中值濾波器的窗口參數(shù)（a）橢圓形濾波窗口的長(zhǎng)軸σ1；（b）橢圓形濾波窗口的短軸σ2.Fig.5 Window parameters of structure－adaptive median filter corresponding to the noisy synthetic data（a）Main axesσ1of the elliptic filter window；（b）The second axesσ2of the elliptic filter window.

加噪合成數(shù)據(jù)的結(jié)構(gòu)自適應(yīng)中值濾波結(jié)果及其得到的噪聲剖面分別如圖6a和圖6d所示.為了對(duì)比本文方法的處理效果，使用經(jīng)典的二維中值濾波器和具有信號(hào)細(xì)節(jié)保護(hù)能力的二維多級(jí)中值濾波器處理加噪合成數(shù)據(jù).取兩種中值濾波器的窗口尺寸與結(jié)構(gòu)自適應(yīng)中值濾波器的最大窗口尺寸相同，即二維中值濾波器選擇9×9的濾波窗口，二維多級(jí)中值濾波器取9點(diǎn)的濾波長(zhǎng)度，兩者對(duì)應(yīng)的濾波結(jié)果以及得到的噪聲剖面分別如圖6b－c和圖6e－f所示.

從圖6a可見(jiàn)，經(jīng)過(guò)結(jié)構(gòu)自適應(yīng)中值濾波后，隨機(jī)噪聲得到很大程度地衰減，反射同相軸的一致性得到增強(qiáng)，剖面整體變得非常干凈；從圖6d可見(jiàn)，結(jié)構(gòu)自適應(yīng)中值濾波器得到的噪聲剖面中只有隨機(jī)分布的干擾噪聲，無(wú)明顯的有效信號(hào)結(jié)構(gòu)，斷層結(jié)構(gòu)在濾波過(guò)程中得到了有效保持.對(duì)比圖6b－c中兩種中值濾波方法的濾波結(jié)果，9×9的二維中值濾波器消除了大部分的隨機(jī)噪聲，但也造成有效信號(hào)的損傷，從其得到的噪聲剖面圖6e可見(jiàn)，濾波過(guò)程引起斷層信息的模糊，也造成反射同相軸的損傷，尤其對(duì)傾斜反射結(jié)構(gòu)的損傷更為明顯；9點(diǎn)濾波長(zhǎng)度的二維多級(jí)中值濾波器雖然有較強(qiáng)的信號(hào)細(xì)節(jié)保護(hù)能力，但是從其得到的噪聲剖面圖6f可見(jiàn)，二維多級(jí)中值濾波器對(duì)隨機(jī)噪聲衰減程度有限.

從上述分析、比較可以看出：在取相同濾波窗口尺度的條件下，相比二維中值濾波器和二維多級(jí)中值濾波器，結(jié)構(gòu)自適應(yīng)中值濾波器能夠根據(jù)地震剖面的信號(hào)結(jié)構(gòu)特征，自適應(yīng)地構(gòu)造與信號(hào)結(jié)構(gòu)相吻合的濾波窗口，在消除隨機(jī)噪聲和保護(hù)斷層等邊緣結(jié)構(gòu)之間找到了很好折中；結(jié)構(gòu)自適應(yīng)中值濾波器的噪聲衰減性能接近9×9的二維中值濾波器，而其對(duì)斷層及有效信號(hào)保護(hù)能力又優(yōu)于9點(diǎn)濾波長(zhǎng)度的二維多級(jí)中值濾波器.

4 實(shí)際資料算例

圖6 加噪合成數(shù)據(jù)的濾波結(jié)果比較（a）結(jié)構(gòu)自適應(yīng)中值濾波結(jié)果；（b）二維中值濾波結(jié)果；（c）二維多級(jí)中值濾波結(jié)果；（d）結(jié)構(gòu)自適應(yīng)中值濾波器得到的噪聲剖面；（e）二維中值濾波器得到的噪聲剖面；（f）二維多級(jí)中值濾波器得到的噪聲剖面.Fig.6 Comparison of filtering results of the noisy synthetic data（a）Result obtained by structure－adaptive median filter；（b）Result obtained by 2Dmedian filter；（c）Result obtained by 2Dmultistage median filter；（d）Noise removed by structure－adaptive median filter；（e）Noise removed by 2Dmedian filter；（f）Noise removed by 2D multistage median filter.

圖7 XX油田實(shí)際地震剖面Fig.7 Real seismic profile from XX oil field

基于理論分析結(jié)果，進(jìn)一步將結(jié)構(gòu)自適應(yīng)中值濾波器用于疊后地震剖面的濾波處理.圖7中給出某油田的疊后純波地震資料，剖面中含有豐富的不連續(xù)結(jié)構(gòu)特征，除了兩條主要的斷層外，在中央部位有一處雜亂反射區(qū)域，在右側(cè)偏下部有一處振幅異常區(qū)域.由于受到隨機(jī)噪聲和采集腳印的干擾，同相軸的連續(xù)性受到一定影響，不利于層位的橫向追蹤、不連續(xù)結(jié)構(gòu)的自動(dòng)檢測(cè)和拾取等后繼的解釋和處理.運(yùn)用本文的結(jié)構(gòu)自適應(yīng)中值濾波器處理該剖面，選擇100道、每道150個(gè)采樣點(diǎn)的子區(qū)域分割；在梯度結(jié)構(gòu)張量的計(jì)算中取高斯濾波器的尺度參數(shù)ρ＝3.0；結(jié)構(gòu)自適應(yīng)中值濾波器的最大濾波尺寸取為Rmax＝4.0，即對(duì)應(yīng)最大為9點(diǎn)的濾波窗口；為了達(dá)到較好的邊緣保持性能，取閾值百分比因子α＝50%，對(duì)應(yīng)的濾波結(jié)果和得到的噪聲剖面分別如圖8a和圖8d所示.作為對(duì)比，同樣采用兩種中值濾波方法處理該剖面.由于二維中值濾波器的濾波窗口選擇得過(guò)大，對(duì)有效信號(hào)損傷嚴(yán)重，而選擇得過(guò)小，會(huì)導(dǎo)致噪聲衰減不夠，故此處取二維中值濾波器的窗口尺寸為5×5，而二維多級(jí)中值濾波器采用9點(diǎn)濾波長(zhǎng)度，即與結(jié)構(gòu)自適應(yīng)中值濾波器的最大濾波尺寸相當(dāng)，對(duì)應(yīng)的濾波結(jié)果和它們得到的噪聲剖面分別如圖8b－c、e－f所示.

圖8 實(shí)際地震資料濾波結(jié)果比較（a）結(jié)構(gòu)自適應(yīng)中值濾波結(jié)果；（b）二維中值濾波結(jié)果；（c）二維多級(jí)中值濾波結(jié)果；（d）結(jié)構(gòu)自適應(yīng)中值濾波器得到的噪聲剖面；（e）二維中值濾波器得到的噪聲剖面；（f）二維多級(jí)中值濾波器得到的噪聲剖面.Fig.8 Comparison of filtering results of real seismic profile（a）Result obtained by structure－adaptive median filter；（b）Result obtained by 2Dmedian filter；（c）Result obtained by 2Dmultistage median filter；（d）Noise removed by structure－adaptive median filter；（e）Noise removed by 2Dmedian filter；（f）Noise removed by 2Dmultistage median filter.

對(duì)比結(jié)構(gòu)自適應(yīng)中值濾波前后的剖面圖7和圖8a可見(jiàn)，結(jié)構(gòu)自適應(yīng)中值濾波過(guò)程消除了大部分的干擾噪聲，改善了反射同相軸的一致性，斷層構(gòu)造得到增強(qiáng)，振幅異常區(qū)域變得清晰.分析結(jié)構(gòu)自適應(yīng)中值濾波得到的噪聲剖面（圖8d）可見(jiàn)，淺層的干擾噪聲非常強(qiáng)，中深層的干擾噪聲相對(duì)較弱，剖面中除隨機(jī)噪聲外還有一些大角度的相干噪聲；對(duì)比原始地震剖面圖7，噪聲剖面中看不到明顯的有效信號(hào)結(jié)構(gòu).分析圖8b中二維中值濾波器的濾波結(jié)果，從視覺(jué)上看噪聲基本被壓制，信噪比得到很大的提高，但是信號(hào)的邊緣結(jié)構(gòu)尤其斷層、雜亂反射等被嚴(yán)重模糊，這一點(diǎn)可以從其得到的噪聲剖面（圖8e）得到印證.從圖8c和圖8f可見(jiàn)，二維多級(jí)中值濾波器僅僅沿同相軸走向有微弱的噪聲衰減作用，濾波能力非常有限.以上結(jié)果表明，結(jié)構(gòu)自適應(yīng)中值濾波器對(duì)淺層的強(qiáng)噪聲和深層的弱噪聲均能夠有效壓制，不僅能夠衰減隨機(jī)噪聲，也能夠消除部分相干噪聲，并且很好地保持了有效信號(hào)和斷層、振幅異常等地層邊緣和細(xì)節(jié)特征；結(jié)構(gòu)自適應(yīng)中值濾波器不僅具有二維中值濾波器優(yōu)異的噪聲衰減能力，而且又不會(huì)因保護(hù)地層邊緣和細(xì)節(jié)特征的需求而像二維多級(jí)中值濾波器那樣引起濾波能力的明顯降低.

對(duì)原始地震剖面和三種濾波方法的處理結(jié)果進(jìn)行頻譜分析，圖9給出相應(yīng)剖面的平均振幅譜.對(duì)比處理前后的頻譜發(fā)現(xiàn)，三種方法處理后的高頻段隨機(jī)噪聲幅度均要低于處理前，且頻譜的主要結(jié)構(gòu)特點(diǎn)沒(méi)有被破壞；然而結(jié)構(gòu)自適應(yīng)中值濾波器和二維多級(jí)中值濾波器均較好地保持了主頻和有用信息頻段的振幅譜，而二維中值濾波器對(duì)主頻和有用信息頻段的振幅譜的損傷較為嚴(yán)重.以上結(jié)果說(shuō)明結(jié)構(gòu)自適應(yīng)中值濾波器能夠有效壓制地震隨機(jī)噪聲，并且?guī)缀醪粨p傷有用信息頻段的頻譜.

圖9 實(shí)際地震資料濾波結(jié)果平均振幅譜對(duì)比（a）圖7的平均振幅譜；（b）圖8（a）的平均振幅譜；（c）圖8（b）的平均振幅譜；（d）圖8（c）的平均振幅譜.Fig.9 Comparison of average spectrum of filtering results of real seismic profile（a）Average spectrum of Fig.7；（b）Average spectrum of Fig.8（a）；（c）Average spectrum of Fig.8（b）；（d）Average spectrum of Fig.8（c）.

5 結(jié) 論

在實(shí)際地震資料的處理中，中值濾波器的應(yīng)用已經(jīng)比較成熟.本文針對(duì)現(xiàn)有的一些中值濾波方法在疊后地震資料隨機(jī)噪聲衰減方面存在的缺陷，提出了一種結(jié)構(gòu)自適應(yīng)中值濾波器.該方法依據(jù)梯度結(jié)構(gòu)張量對(duì)分析鄰域內(nèi)地層結(jié)構(gòu)規(guī)則程度的估計(jì)情況，引入了地震信號(hào)的兩種結(jié)構(gòu)特征的置信度量，利用它們調(diào)控結(jié)構(gòu)自適應(yīng)中值濾波器的濾波窗口的尺度和形狀，使得濾波窗口能夠最佳地匹配處理鄰域內(nèi)的地層結(jié)構(gòu)特征.通過(guò)合成模型和實(shí)際地震資料的處理結(jié)果表明，相比經(jīng)典的二維中值濾波器及二維多級(jí)中值濾波器，本文方法不僅能夠有效地衰減隨機(jī)噪聲、提高剖面的信噪比，而且能夠較好地保持有效信號(hào)和斷層、裂縫等地層邊緣和細(xì)節(jié)特征.本文方法可作為一個(gè)預(yù)處理過(guò)程，濾波結(jié)果進(jìn)一步用于層位的自動(dòng)追蹤、對(duì)噪聲敏感屬性的提取等后繼的解釋和處理流程.

（References）

［1］ 李月，楊寶俊，趙雪平等.檢測(cè)地震勘探微弱同相軸的混沌振子算法.地球物理學(xué)報(bào)，2005，48（6）：1428－1433.Li Y，Yang B J，Zhao X P，et al.An algorithm of chaotic vibrator to detect weak events in seismic prospecting records.Chinese J.Geophys.（in Chinese），2005，48（6）：1428－1433.

［2］ 高靜懷，毛劍，滿(mǎn)蔚仕等.疊前地震資料噪聲衰減的小波域方法研究.地球物理學(xué)報(bào)，2006，49（4）：1155－1163.Gao J H，Mao J，Man W S，et al.On the denoising method of prestack seismic data in wavelet domain.Chinese J.Geophys.（in Chinese），2006，49（4）：1155－1163.

［3］ Neelamani R，Baumstein A I，Gillard D G，et al.Coherent and random noise attenuation using the curvelet transform.The Leading Edge，2008，27（2）：240－248.

［4］ 劉洋，F(xiàn)omel S，劉財(cái)?shù)?高階seislet變換及其在隨機(jī)噪聲消除中的應(yīng)用.地球物理學(xué)報(bào)，2009，52（8）：2142－2151.Liu Y，F(xiàn)omel S，Liu C，et al.High－order seislet transform and its application of random noise attenuation.Chinese J.Geophys.（in Chinese），2009，52（8）：2142－2151.

［5］ 陸文凱.基于離散余弦變換的地震隨機(jī)噪聲壓制技術(shù).石油地球物理勘探，2011，46（2）：202－206.Lu W K.Seismic random noise suppression based on the discrete cosine transform.Oil Geophysical Prospecting（in Chinese），2011，46（2）：202－206.

［6］ 林紅波，李月，徐學(xué)純.壓制地震勘探隨機(jī)噪聲的分段時(shí)頻峰值濾波方法.地球物理學(xué)報(bào)，2011，54（5）：1358－1366.Lin H B，Li Y，Xu X C.Segmenting time－frequency peak filtering method to attenuation of seismic random noise.Chinese J.Geophys.（in Chinese），2011，54（5）：1358－1366.

［7］ Chopra S，Marfurt K J.Seismic attributes–a historical perspective.Geophysics，2005，70（5）：3SO－28SO.

［8］ Chopra S，Marfurt K J.Emerging and future trends in seismic attributes.The Leading Edge，2008，27（3）：298－318.

［9］ Chopra S，Misra S，Marfurt K J.Coherence and curvature attributes on preconditioned seismic data.The Leading Edge，2011，30（4）386－393.

［10］ Luo Y，Marhoon M，Al－Dossary S，et al.Acquistion Processing—Edge－preserving smoothing and applications.The Leading Edge，2002，21（2）：136－158.

［11］ AlBinHassan N M，Luo Y，Al－Faraj M N.3Dedgepreserving smoothing and applications.Geophysics，2006，71（4）：P5－P11.

［12］ Lu Y H，Lu W K.Edge－preserving polynomial fitting method to suppress random seismic noise.Geophysics，2009，74（4）：V69－V73.

［13］ 楊培杰，穆星，張景濤.方向性邊界保持?jǐn)鄬釉鰪?qiáng)技術(shù).地球物理學(xué)報(bào)，2010，53（12）：2992－2997.Yang P J，Mu X，Zhang J T.Orientational edge preserving fault enhance.Chinese J.Gephys.（in Chinese），2010，53（12）：2992－2997.

［14］ 劉財(cái)，王典，劉洋等.二維多級(jí)中值濾波技術(shù)在隨機(jī)噪聲消除中的應(yīng)用初探.石油地球物理勘探，2005，40（2）：163－167.Liu C，Wang D，Liu Y，et al.Preliminary study of using 2D multi－level median filtering technique to eliminate random noises.Oil Geophysical Prospecting（in Chinese），2005，40（2）：163－167.

［15］ Liu C，Liu Y，Yang B J，et al.A 2Dmultistage median filter to reduce random seismic noise.Geophysics，2006，71（5）：V105－V110.

［16］ Fehmers G C，H?cker C F W.Fast structural interpretation with structure－oriented filtering.Geophysics，2003，68（4）：1286－1293.

［17］ Lavialle O，Pop S，Germain C，et al.Seismic fault preserving diffusion.Journal of Applied Geophysics，2007，61（2）：132－141.

［18］ 孫夕平，杜世通，湯磊.相干增強(qiáng)各向異性擴(kuò)散濾波技術(shù).石油地球物理勘探，2004，39（6）：651－655.Sun X P，Du S T，Tang L.Coherent－enhancing anisotropic diffusion filtering technique.Oil Geophysical Prospecting（in Chinese），2004，39（6）：651－655.

［19］ 王旭松，楊長(zhǎng)春.對(duì)地震圖像進(jìn)行保邊濾波的非線(xiàn)性各向異性擴(kuò)散算法.地球物理學(xué)進(jìn)展，2006，21（2）：452－457.Wang X S，Yang C C.An edge－preserving smoothing algorithm of seismic image using nonlinear anisotropic diffusion equation.Progress in Geophys.（in Chinese），2006，21（2）：452－457.

［20］ 張爾華，王偉，高靜懷等.非線(xiàn)性各向異性擴(kuò)散濾波器用于三維地震資料噪聲衰減與結(jié)構(gòu)特征增強(qiáng).地球物理學(xué)進(jìn)展，2010，25（3）：866－870.Zhang E H，Wang W，Gao J H，et al.Non－linear anisotropic diffusion filtering for 3Dseismic noise removal and structure enhancement.Progress in Geophys.（in Chinese），2010，25（3）：866－870.

［21］ Gómez G.Local smoothness in terms of variance.Proceedings of the BMVC，2000：815－824.

［22］ Marfurt K J，Sudhaker V，Gersztenkorn A，et al.Coherency calculations in the presence of structural dip.Geophysics，1999，64（1）：104－111.

［23］ Randen T，Monsen E，Signer C，et al.Three－dimensional texture attributes for seismic data analysis.70th Annual International Meeting，SEG Expanded Abstracts，2000：668－671.

［24］ Wang Y E，Luo Y，Al－faraj M N.Inverse－vector approach for smoothing dips and azimuths.Geophysics，2008，73（6）：P9－P14.

［25］ Gao D L.Latest developments in seismic texture analysis for subsurface structure，facies，and reservoir characterization：A review.Geophysics，2011，76（2）：W1－W13.

［26］ Bakker P.Image structure analysis for seismic interpretation［Ph.D′s thesis］.Netherlands：Technische Universiteit Delft，2002.

Random seismic noise suppression via structure－adaptive median filter

WANG Wei1，GAO Jing－Huai1，CHEN Wen－Chao1＊，ZHU Zhen－Yu2
1 School of Electronic and Information Engineering，Xi′an Jiaotong University，Xi′an 710049，China
2 CNOOC Research Institute，Beijing100027，China

This paper presents a structure－adaptive median filter for reducing random seismic noise which preserves seismic edges and details such as faults and fractures.By considering seismic horizons as line－like structures，this filtering framework relies on the computation of the Gradient Structure Tensors，which provides dips of seismic events and the regularity of local seismic structures，based on which two confidence measures are defined for different seismic structures.The structure－adaptive median filter adjusts the shape of the filter kernel according to the two confidence measures and aligns the filtering kernel along local dips to make the filtering kernel optimally matched with different local geological features.The proposed filter has been applied to both the synthetic and real data，and compared with two widely used median filtering schemes.The processing and comparison results show that the proposed structure－adaptive median filter is more suitable to balance suppressing random seismic noise and preserving signals which preserves seismic edges and details while enhancing the coherence of seismic horizons，and improves the quality of seismic images significantly.

Median filter，Adaptive filtering，Noise suppression，Structure preserving，Gradient structure tensor

10.6038／j.issn.0001－5733.2012.05.030

P631

2011－04－06，2011－06－28收修定稿

國(guó)家自然科學(xué)基金重點(diǎn)項(xiàng)目（40730424）、國(guó)家自然科學(xué)基金面上項(xiàng)目（40674064）、國(guó)家油氣重大專(zhuān)項(xiàng)（2011ZX05023－005－009）資助.

王偉，男，1983年生，現(xiàn)為西安交通大學(xué)在讀博士生，主要從事地震信號(hào)處理和地震屬性分析方面的研究.E－mail：feitianliuhuo＠gmail.com

＊通訊作者陳文超，男，1970年生，副教授，2000年于西安交通大學(xué)獲博士學(xué)位，主要從事小波分析及盲信號(hào)分離在地震信號(hào)處理、解釋中的應(yīng)用研究.E－mail：wencchen＠m(xù)ail.xjtu.edu.cn

王偉，高靜懷，陳文超等.基于結(jié)構(gòu)自適應(yīng)中值濾波器的隨機(jī)噪聲衰減方法.地球物理學(xué)報(bào)，2012，55（5）：1732－1741，

10.6038／j.issn.0001－5733.2012.05.030.

Wang W，Gao J H，Chen W C，et al.Random seismic noise suppression via structure－adaptive median filter.Chinese J.Geophys.（in Chinese），2012，55（5）：1732－1741，doi：10.6038／j.issn.0001－5733.2012.05.030.

（本文編輯 劉少華）