基于簇檢測的地震數(shù)據(jù)脈沖噪聲壓制方法

姜弢， 胡秋月

(吉林大學儀器科學與電氣工程學院， 長春 130012)

在地震勘探[1]、地震監(jiān)測[2-3]和使用地震數(shù)據(jù)對火車、飛機和船舶等進行檢測[4-5]時，地震記錄中常常會伴有強烈的脈沖噪聲干擾，大大降低了資料信噪比和分辨率，并嚴重影響地震數(shù)據(jù)分析的后續(xù)工作。因此，脈沖噪聲抑制是地震數(shù)據(jù)處理的重要步驟之一[6]。脈沖噪聲主要來自地震儀附近的噪聲源，例如近處的人為活動、物體的突然放電、風吹樹動、地面交通事件或其他物理過程等[7]。此外，這些脈沖噪聲常表現(xiàn)為聚集狀態(tài)，持續(xù)時間可能為幾十到幾百毫秒。

史維等[8]提出一種自適應多尺度形態(tài)濾波算法壓制脈沖噪聲，但該方法僅適合于處理幅值遠高于有效信號且持續(xù)時間極短的脈沖噪聲。劉軍等[9]提出基于混沌蟻群優(yōu)化小波閾值法壓制地磁信號噪聲，但其僅適應于數(shù)據(jù)中隨機噪聲為主導噪聲的情況。江友華等[10]使用改進的最小均方算法壓制局部放電脈沖信號中的窄帶干擾，換言之，該算法可以抑制脈沖信號，但算法有效的前提是信號與噪聲之間存在明顯的相關性差異，否則算法失效。中值濾波器是被廣泛應用的脈沖噪聲抑制方法[11-15]。但是，中值濾波器對數(shù)據(jù)中所有采樣點都采用無差別濾波策略，當應用長濾波窗口去除長持續(xù)時間的脈沖噪聲時，會嚴重損壞有效信號。近些年來，許多改進算法被提出，主要有兩類。第一類是使用某種手段將有效信號的波形數(shù)據(jù)展平，再對展平后的數(shù)據(jù)應用中值濾波器則可將損壞最小化。如Gan等[16]根據(jù)噪聲數(shù)據(jù)估算出局部斜率，提出面向結構的中值濾波器。Chen 等[17]基于相鄰跡線之間的可預測性將波形數(shù)據(jù)展平后進行濾波。另一類是使用某種方法對有效信號和噪聲進行判別后應用中值濾波方法。如Liu等[6]采用經(jīng)參考中值濾波器處理后的數(shù)據(jù)的絕對平均值作為噪聲檢測閾值，然后對數(shù)據(jù)進行平滑。Xu等[18]將當前像素與參考圖像的像素進行比較，區(qū)分出損壞的像素后進行濾波。Chen等[19]提出了一種利用正態(tài)分布原理和局部強度統(tǒng)計的噪聲檢測器。Arastehfar等[20]提出判決中值濾波器，它根據(jù)相鄰采樣點的幅值變化對脈沖噪聲和有效信號進行判別，然后應用中值濾波器平滑噪聲部分。Hashemi等[21]應用該判決中值濾波器壓制地震數(shù)據(jù)脈沖噪聲。但當有效信號和脈沖噪聲的相鄰采樣點幅值變化相近時，則判決中值濾波器決策失敗。

顯然，噪聲的判決和檢測是壓制聚集脈沖噪聲的關鍵?；诿}沖噪聲和有效信號在短時能量變化上的差異，現(xiàn)應用長短時平均值比法(short-term average/long-term average, STA/LTA)[22-24]對噪聲進行檢測，并分析該方法的參數(shù)選取原則和存在的問題，然后針對噪聲漏檢問題，提出同一聚集的脈沖噪聲判別方法，最后，結合中值濾波器自適應的伸縮濾波窗口對脈沖噪聲進行濾波。通過將該脈沖噪聲壓制方法應用于合成數(shù)據(jù)和真實數(shù)據(jù)，以有效抑制地震數(shù)據(jù)脈沖噪聲簇，并保留有效信號能量。

1 脈沖噪聲簇特征

圖1(a)為某震源產(chǎn)生的持續(xù)時間為2～15 s 的chirp信號，圖1(b)為某道地震檢波器記錄的實測地震信號，可以看出chirp信號能量大幅衰減，并存在強烈的脈沖噪聲干擾。分析圖1(b)發(fā)現(xiàn)，脈沖噪聲可能表現(xiàn)為孤立脈沖或脈沖聚集的形式。孤立脈沖由一個或幾個峰值組成，通常持續(xù)幾毫秒，如圖1(c)所示。聚集的脈沖由多個幅值和持續(xù)時間各異的孤立脈沖在時間軸上無規(guī)則緊密排列形成，持續(xù)時間可能為幾十到幾百毫秒，將其稱之為脈沖噪聲簇，如圖1(d)所示。圖1(b)中存在多個持續(xù)時間各異的脈沖噪聲簇。顯然，應用中值濾波器濾除長持續(xù)時間的脈沖噪聲簇會損壞有效信號。此外，脈沖噪聲簇的波形隨機起伏，噪聲簇中不同數(shù)據(jù)段的能量變化速率不同。

圖1 脈沖噪聲簇的特征Fig.1 Characteristics of impulse noise clusters

雖然脈沖噪聲簇相較于孤立脈沖持續(xù)時間長、能量強，但相較于有效信號仍具有持續(xù)時間短、幅值高的特點，即脈沖噪聲簇的短時能量變化相較于有效信號的變化的快，故而應用STA/LTA方法進行脈沖噪聲簇檢測。

2 脈沖噪聲簇壓制方法

地震數(shù)據(jù)脈沖噪聲簇壓制方法分為以下三步：①應用STA/LTA方法對脈沖噪聲簇進行檢測；②進行同一脈沖噪聲簇判別；③結合中值濾波器，取長度為2L+1的中值濾波窗口對檢測到的脈沖噪聲簇逐一進行濾波，其中L為拾取的脈沖噪聲簇的持續(xù)時間。

步驟①和步驟②的目的是完成脈沖噪聲簇檢測，步驟③是完成平滑處理。中值濾波器的應用與常規(guī)無異，本文不再贅述。下面對該方法的前兩個步驟分別進行分析和說明。

2.1 基于STA/LTA方法的脈沖噪聲簇檢測

STA/LTA方法的基本原理是：用STA和LTA之比來反映信號幅度、頻率等特征的變化[25]。被用于脈沖噪聲簇檢測時，噪聲部分的STA/LTA值R會有一個大的突變，當R大于某一設定閾值(Thr)時，則被判定為脈沖噪聲。

為了避免監(jiān)測數(shù)據(jù)長久占用內(nèi)存，建議使用遞歸STA/LTA方法[26]。但其中Thr設置與大多數(shù)文獻中地震事件檢測和震相拾取所應用的經(jīng)驗閾值不同，為增強算法的自主性，采用統(tǒng)計學思想，將Thr用R的平均值μ和標準差表示，即

(1)

式(1)中：n為地震數(shù)據(jù)采樣點序號，n=1,2,…,N；N為總采樣點數(shù)；a為常數(shù)。當脈沖噪聲簇能量強，信噪比低時，a取小值。反之，當脈沖噪聲簇能量較弱，信噪比高時，a取大值。

應用不同類型的特征函數(shù)對同一脈沖噪聲簇進行檢測，含噪數(shù)據(jù)如圖2(a)所示，結果如表1所示。其中，脈沖噪聲漏檢率(MDR)計算方法為

(2)

式(2)中：IN為脈沖噪聲數(shù)據(jù)；IND為檢測到的脈沖噪聲數(shù)據(jù)；m為噪聲采樣點序號，m=1,2,…,Q,…,M；M為噪聲總采樣點數(shù)，且1≤Q≤M。

特征函數(shù)(CF)如式(3)～式(5)所示

CF1=x2(n)-x(n-1)x(n+1)

(3)

CF2=x2(n)+[x(n)-x(n-1)]2

(4)

CF3=x2(n)+K[x(n)-x(n-1)]2,

(5)

式中：xn為地震數(shù)據(jù)。

它們包含數(shù)據(jù)的幅值和頻率變化特征[27]，CF1是一種常用的增強數(shù)據(jù)幅值和頻率變化的特征函數(shù)，CF2是利用數(shù)據(jù)的幅值和一階導數(shù)計算得到的特征函數(shù)，CF3是對數(shù)據(jù)的幅值和一階導數(shù)加權后計算得到的特征函數(shù)，加權因子為K。注意，選取短的短時窗口長度(Nsta)，以增強該算法對數(shù)據(jù)能量變化的敏感度。這有利于精準檢測噪聲和避免有效信號的誤檢，適應于數(shù)據(jù)中存在多個持續(xù)時間各不相同的脈沖噪聲簇的情況。因為Nsta長，會減小噪聲和有效信號的R的差距，不利于區(qū)分噪聲和有效信號。此外，對數(shù)據(jù)中短持續(xù)時間的脈沖噪聲簇來說，Nsta長，使得噪聲首部數(shù)據(jù)對應的R緩慢增大，噪聲尾部數(shù)據(jù)對應的R緩慢減小，從而導致噪聲首部漏檢和噪聲尾部相鄰有效信號誤檢，如圖2(b)中實線框和虛線框內(nèi)數(shù)據(jù)所示。當Nsta過長時，甚至會造成噪聲完全漏檢。經(jīng)驗證，Nsta約取短持續(xù)時間脈沖噪聲簇的長度的1/20～1/100。

不同特征函數(shù)的脈沖噪聲漏檢率如表1所示。由表1得到，各特征函數(shù)檢測結果相近，脈沖噪聲簇能量漏檢20%左右。這是由于脈沖噪聲簇具有不同數(shù)據(jù)段能量變化速率不同的特征，從而導致使用STA/LTA方法時，高能量變化的噪聲部分被檢測，而漏檢能量變化相對較低的噪聲，如圖2(c)所示。此外，若此時應用中值濾波器對檢測到的噪聲進行平滑，則由于漏檢噪聲的存在，會使得噪聲數(shù)據(jù)的替代值表現(xiàn)為噪聲的能量特征。

表1 不同特征函數(shù)的脈沖噪聲漏檢率

2.2 同一脈沖噪聲簇判別

分析被檢測到的噪聲和漏檢的噪聲的時間長度，得到漏檢的噪聲的持續(xù)時間基本小于其相鄰的已被檢測到的噪聲的最長持續(xù)時間，如圖3所示。因此，根據(jù)時間特性對檢測到的噪聲是否為同一脈沖噪聲簇進行判斷，以解決漏檢及由漏檢引起的中值替代值的問題。方法如下：當兩個相鄰的已檢測到的噪聲的時間間隔小于這兩個噪聲的最長持續(xù)時間，即若

(6)

則將這兩個噪聲以及其中間間隔部分的數(shù)據(jù)共同視為一個脈沖噪聲簇。

加入判斷后的脈沖噪聲簇檢測結果如表1和圖2(d)所示。此時脈沖噪聲簇中能量變化相對較低的噪聲部分被成功檢測，漏檢率與未進行同一脈沖噪聲簇判別前相比大大降低，僅在脈沖噪聲簇首部和尾部有少量漏檢，滿足信號處理需求。CF2的漏檢率相對更低，故而在遞歸STA/LTA方法中建議使用該特征函數(shù)。此外，應用中值濾波器對噪聲進行平滑，此時中值替代值主要表現(xiàn)為有效信號的能量特征。

圖2 脈沖噪聲簇檢測結果(以CF3為例)Fig.2 Detection results of impulse noise clusters (take CF3 as an example)

圖3 已檢到噪聲和漏檢噪聲的時間長度對比Fig.3 Comparison of the length of time between detected noise and missed noise

3 實驗與分析

為驗證該算法的去噪效果，進行合成數(shù)據(jù)和真實數(shù)據(jù)實驗，并用式(7)和式(8)所示的信噪比(SNR)和信號相似度(XCOR)衡量算法性能，XCOR在0～1，且越接近1表示結果越好。

(7)

(8)

式中：s(n)為不含噪聲的地震信號；y(n)為去噪后信號。

3.1 合成地震數(shù)據(jù)降噪

使用合成數(shù)據(jù)來評估該方法的降噪效果。給定不含噪聲的震地信號為chirp信號，掃描時間為2～10 s、頻率為0～200 Hz，如圖4(a)所示。在其中加入噪聲干擾，分別為脈沖噪聲簇、孤立脈沖以及隨機噪聲。考慮到脈沖噪聲簇的多樣性，則加入持續(xù)時間分別為0.338 s和0.12 s的長持續(xù)時間脈沖噪聲簇和短持續(xù)時間脈沖噪聲簇，加入噪聲后數(shù)據(jù)如圖4(b)所示。

圖4 各方法對合成地震數(shù)據(jù)降噪效果對比Fig.4 Comparison of noise reduction effects of various methods on synthetic seismic data

對此合成地震數(shù)據(jù)記錄分別使用中值濾波器、決策中值濾波器和本文方法進行降噪處理，3種方法處理后的定性和定量結果分別如圖4(c)～圖4(e)和表2所示。

由圖4可以明顯看出，使用本文方法進行噪聲壓制后，脈沖噪聲簇稍有殘留，且非脈沖噪聲簇部分的有效信號沒有被損壞，但使用決策中值濾波器和中值濾波器后，脈沖噪聲不但有很大的殘留，而且非脈沖噪聲部分的有效信號也存在不同程度的損壞。對比表2結果可以得到，使用本文方法去噪后數(shù)據(jù)的SNR和XCOR遠優(yōu)于其他兩種常用方法。顯然，本文方法對脈沖噪聲簇的壓制效果更好。此外，由圖4(e)得到孤立脈沖噪聲也被成功壓制，故而本文方法也適用于孤立脈沖噪聲的抑制。

表2 不同方法壓制合成數(shù)據(jù)脈沖噪聲簇定量結果

3.2 實測地震數(shù)據(jù)降噪

應用真實地震記錄以進一步研究本文方法的降噪性能，實驗數(shù)據(jù)來源為北京時間2019年4月4日在吉林省長春市龍嘉機場應用地震檢波器實際監(jiān)測的飛機信號。圖5(a)為某道地震檢波器記錄的低空飛機引起的真實振動信號。此段數(shù)據(jù)中共有6次飛機事件，如圖中點線框內(nèi)所示。圖5(a)中實線框內(nèi)為某一脈沖噪聲簇，其放大圖如圖5(b)所示。

對此真實低空飛機數(shù)據(jù)分別使用中值濾波器、決策中值濾波器和本文方法進行脈沖噪聲簇衰減，定性和定量結果分別如圖5(c)～圖5(e)和表3所示。其中，取圖5(a)中6個已知的飛機信號作為不含噪信號，從而進行表3數(shù)據(jù)計算。

圖5 各方法對實測地震數(shù)據(jù)降噪效果對比Fig.5 Comparison of noise reduction effects of various methods on measured seismic data

表3 不同方法壓制實測數(shù)據(jù)脈沖噪聲簇定量結果

由圖5和表3可以得到，中值濾波器和決策中值濾波器雖然對脈沖噪聲簇有抑制作用，但會嚴重損壞飛機信號，導致濾波后飛機信號淹沒于噪聲之中。而本文方法基本完全壓制脈沖噪聲簇，且僅對飛機信號稍有損壞。本文方法相比于其他方法具有明顯優(yōu)勢，可在確保有效信號損壞程度最小的同時，大大衰減脈沖噪聲簇能量。

4 結論

地震數(shù)據(jù)中脈沖噪聲簇干擾嚴重?？紤]到使用中值濾波器進行噪聲抑制會大大損壞有效信號，改進的決策中值濾波器無法區(qū)分相鄰噪聲點和相鄰信號點幅度值變化相近的數(shù)據(jù)，提出一種基于簇檢測的脈沖噪聲壓制方法，經(jīng)理論分析及實驗，得出如下結論。

(1)提出的基于STA/LTA和同一噪聲簇判別方法可以有效判別和檢測脈沖噪聲簇，保證噪聲的低漏檢和有效信號的低誤檢。

(2)合成和真實數(shù)據(jù)實例表明，本文方法相比于其他方法可以顯著提高數(shù)據(jù)質量，在較大程度地去除脈沖噪聲簇的同時，確保有效信號損壞程度最小。

(3)本文方法能夠自主壓制地震數(shù)據(jù)中多樣的脈沖噪聲簇，并同樣適用于孤立脈沖噪聲的抑制。此外，對于其他領域，若其數(shù)據(jù)存在有效信號與脈沖噪聲能量變化差異大的特點，則也可使用本文方法進行脈沖噪聲壓制。