一種抗脈沖干擾的數學形態濾波方法

0數學形態濾波地磁觀測容易受到外部環境的干擾。地磁站的日變觀測數據通常會受到車輛、閃電和未知鐵磁性材料的極大干擾。在某些情況下，噪聲干擾可能大于真實的地磁信號本身。地磁觀測儀器產生的儀器噪聲也會干擾地磁數據。對地磁數據的噪聲處理可以減少干擾對地磁質數據質量的影響，成為地磁數據預處理研究的重要內容之一。數字濾波器作為最常用的信號處理手段可對地磁信號進行預處理,常用的數字濾波器有IIR濾波器、FIR濾波器等.但在數字濾波器設計中存在時滯、相移等缺陷,并且當信號與噪聲差異較小且難以區分,或者干擾信號的頻帶分布較寬等情況下,數字濾波方法的應用受到限制.因此,我們嘗試從干擾形態上尋找一種新的方法—數學形態濾波,來對受干擾地磁信號進行預處理.數學形態學是從圖像處理理論處理演變而來的新方法.作為一種重要的非線性濾波器,與傳統的數字濾波器相比,數字形態濾波通過數學形態變換將復雜信號分解為具有物理意義的各個部分,在處理伴隨有較強噪聲干擾,甚至發生嚴重畸變的情況下,其基本形態依舊能別識別、重構及增強.Maragos等人在數學形態框架下系統總結形態濾波器集合理論,并闡明形態濾波器與線性平移不變、排序統計和層疊濾波器的關系,為數學形態學的發展奠定了堅實基礎.目前數學形態學主要的應用領域包括圖像處理、模式識別、電力系統監測、心(腦)電信號干擾抑制、機械振動信號特征提取等方面[9,10,11,12,13,14,15,16].但是,數學形態濾波應用在地磁信號預處理方面的工作開展較少.地磁臺站每天所產出的地磁數據包含平靜變化信號和疊加在平靜變化上的不規則擾動變化信號.平靜期變化的頻率周期一般以24、12、8、6小時的周期為主;而不規則擾動主要有磁擾和脈動兩大類,其頻率分布較廣,成分較為復雜,頻率分布從10-3～102Hz的數量級.通過儀器采集的地磁數據中除了地磁信號本身以外,還會記錄下外界環境干擾所產生的電磁場,以及儀器本身噪聲.外界環境干擾主要包括車輛等鐵磁性物質在距離儀器探頭較近處運動產生的脈沖狀干擾;雷暴雨天氣所引起的連續的脈沖狀干擾;部分大城市軌道交通距離地磁觀測臺站較近時,由于軌道交通漏電電流所產生的軌道交通干擾;以其一些未知的人工干擾等.其中由車輛等鐵磁性物質運動所產生的脈沖狀干擾最為常見,其干擾幅度一般較大,干擾持續時間從幾秒到幾十秒不等,而且干擾發生較為隨機.由于是脈沖狀干擾,其頻譜分布較寬,采用通常的數字濾波器進行抑制效果并不理想.本文根據數學形態濾波原理,提出一種地磁日變數據干擾抑制的新方法.討論了不同結構元素對儀器干擾和暫態脈沖干擾的抑制效果.通過數值模擬,對地磁日變數據中添加隨機噪聲和不同形態的脈沖干擾噪聲進行模擬抑制.最后對實際地磁觀測儀器的儀器噪聲和由車輛開行所產生的脈沖干擾進行處理,獲得較好的應用效果.鑒于數學形態學用于地磁數據處理是一個新的應用領域,這方面的研究工作更需要進一步深入.1數學濾波的基本原理1.1腐蝕運算定義數學形態學是一種以集合理論與邏輯運算為基礎的非線性圖像處理理論.其基本思想是用具有一定形態的結構元素作為“探針”元素對目標信號進行變換和匹配,保持信號基本形態,抑制不相干信號,以達到提取信號、保持細節和抑制噪聲的目的.假設地磁信號f(n)與結構元素g(n)分別定義在實數上的一維離散函數.f(n)的定義域D[f]∈{l,2,3,…,N},g(n)的定義域D[g]∈{0,1,2,3,…,P},且N?P.g(n)關于原點對稱函數gs(n)=g(-n):則腐蝕運算定義:膨脹運算定義式中,符號inf和sup分別表示函數的極小和極大值運算.開運算定義:閉運算定義這里符號?和·分別表示開運算和閉運算.開運算采用先膨脹后腐蝕運算,其作用是分離信號中比結構元素小的孤立部分,用于抑制信號中的正脈沖噪聲;閉合運算采用先腐蝕后膨脹運算,其作用是補缺與內部連通,用于抑制信號的負脈沖噪聲.1.2結構元素的選取數學形態學方法是基于結構元素填充探測的思想,結構元素的選取是影響濾波器性能的關鍵因素,常見的結構元素為扁平型、矩形、三角形、半圓形、正(余)弦型等.大量仿真經驗表明,較長的結構元素能夠獲得較好的低通效果,濾波的效果也較為明顯,但結構元素的結構越長,復雜性越高,計算量也會相應增大.地磁信號易受電磁瞬態脈沖干擾,因此在結構元素的選取上,結構元素應大于瞬態脈沖長度.而對于儀器背景噪聲等隨機噪聲干擾,結構元素長度為地磁信號采樣周期的5至10倍左右可獲得較好的效果.在下一章節中將會具體討論不同結構元素對于這些干擾的抑制效果..1.3混合濾波器與交替混合濾波器級聯開、閉運算可構建交替濾波器foc(n)和fco(n).同時,還能通過開、閉基本運算構建混合濾波器:通過交替濾波器還能構建較為復雜的交替混合濾波器:通常經過交替濾波器的計算,信號會向上或者向下偏移,這跟開閉運算輸出信號的特性有關,為解決這個問題,一般采用混合濾波器及交替混合濾波器.對于相同長度和寬度的結構元素而言,計算最為快速的是交替濾波器,而交替混合濾波器計算量最大,當結構元素更為復雜時,這樣的運算性能上的差異則更為明顯.本文在對干擾進行抑制時,考慮到對噪聲抑制能取得最優的抑制效果,采用的是交替混合濾波器.2模擬計算2.1仿真模擬日變形態曲線地磁場的正常日變化主要由1～6次諧波組成,其周期分別為24,12,8,6,4.8,4h,振幅依次減小.因此,通過傅立葉級數進行前6次諧波可較好地擬合地磁場日變形態.數值模擬后的日變形態曲線作為對干擾抑制效果評判的基礎.地磁日變曲線的擬合函數f(t)通過傅立葉級數,可表示成如下形式:式中a為直流分量,ck表示k次諧波振幅大小,φk表示k次諧波的初相.以中國地磁臺網中心獅泉河臺2008年3月25日GM-4磁通門磁力儀產出的地磁日變數據為例,1Hz采樣的數據通過傅立葉級數擬合,能夠非常好的模擬出當天的日變形態特征.如圖1所示,紅色實線為6次諧波擬合數據,虛線為原始觀測數據.通過對地磁日變信號仿真模擬,可以定量的對數學形態濾波器的濾波效果進行評判.2.2結構元素的數學形態濾波地磁儀器在數據采集過程中常常伴有具有白噪聲特性的儀器噪聲.在通過6階擬合的地磁日變信號中加入白噪聲隨機干擾信號.定義信噪比公式為:式中f(t)為信號,為干擾抑制后的信號.下面探討扁平型、余弦型、三角型這三種結構元素對隨機干擾的抑制效果的影響比較.在仿真地磁日變信號中加入均值為0.018,方差為0.1的高斯分布序列.首先,以扁平結構元素為例,在不同長度的扁平元素結構下,含有隨機干擾的地磁日變信號的抑制效果如圖所示.圖2(a)中為混有隨機干擾的地磁日變信號,信噪比為28.17dB.圖2(b)～圖2(f)分別為施加了寬度為20,50,100,150,200的扁平元素進行數學形態濾波.其中以圖2(e)的濾波效果最好.可見,選擇寬度為100～150左右扁平結構元素能夠抑制隨機干擾,保存有效地磁日變信號.其次,以寬度為150,高度不同的余弦結構元素對混有隨機噪聲地磁日變信號進行抑制.如圖所示,圖3(a)～(d)分別為施加了寬度為150,高度為0.1,0.07,0.05,0.02余弦結構元素進行數學形態濾波.其中以圖3(c)的濾波效果最好.可見,當噪聲方差δ與余弦結構元素的高度h比在2:1左右,即δ/h=2時,能夠較好的抑制隨機干擾,保存有效地磁日變信號.通過選取三種不同高度、寬度結構元素進行數學形態濾波,對包含隨機干擾的地磁日變數據的抑制,我們發現,使用交替混合濾波器的情況下,結構元素為余弦序列或三角序列,在寬度為150左右,高度分別為0.05和0.1左右的情況下,對隨機噪聲的抑制效果最好,其抑制后地磁信號的信噪比能提高將近一倍左右,并且能夠保留有效的地磁日變信號.2.3脈沖狀干擾下的地磁日變曲線地磁觀測臺站易受鐵磁性物質的運動所產生的脈沖狀干擾的影響,這種脈沖狀干擾尤其以離觀測臺站較近的車輛開行過最為常見.其他種類的脈沖狀干擾包括高速鐵路的影響,雷暴天氣的原因等.一般的脈沖狀干擾4～20s左右完成.脈沖干擾的幅度與鐵磁性物體的磁化率,體積以及運行速度以及與觀測儀器探頭之間距離等參數有關.干擾幅度從幾個nT到幾十nT都有可能.因此我們在擬合出的地磁日變曲線中添加脈沖噪聲.如圖5所示,在擬合的地磁日變形態數據的添加脈沖狀干擾,其中①號脈,沖持續時間4s,干擾幅度為3nT;②號脈,沖持續時間為19s,干擾幅度為15nT;③號脈沖,干擾幅度為-3nT,干擾時間持續8s根據對①號脈沖的處理結果,選擇了濾波誤差較小的三角結構元素對②號脈沖進行處理.根據②號脈沖的持續持劍較長的特性.選取三角元素的寬度20,圖7為濾波效果圖,通過選取合適的三角元素,最大濾波誤差只有0.024nT.最后對③號脈沖選取三角元素的寬度為10的結構元素進行形態濾波,濾波效果如圖8所示,其最大濾波誤差只有0.015nT.可見對于濾除尖脈沖干擾,所選的形態濾波結構元素在形態上應接近干擾的形態特征,并且結構元素寬度T與干擾持續的時間Ts較為接近時,即T/Ts≈1時,濾波效果最為理想.3地震信號的實際處理3.1基于構造模擬的儀器噪聲抑制算法目前十五地磁數字前兆采集儀器中使用較為廣泛的地磁秒數據相對觀測儀為GM-4磁通門磁力儀.圖9為儀器檢測中的GM-4磁通門磁力儀在零磁空間中7200s的儀器噪聲曲線.圖中儀器噪聲曲線的方差為0.05,均值接近0.假設GM-4采集的信號由地磁日變信號和儀器的噪聲所組成,則y(t)=f(t)+s(t),(t∈{1,2,3…N})(11)式中,f(t)為地磁日變信號,s(f)為儀器噪聲,我們使用圖9中的儀器噪聲數據添加到f(t)中構成y(t)為含有儀器噪聲的地磁日變信號.從前小結的模擬抑制方法中可以得出,數學形態濾波可以很好的抑制隨機噪聲的干擾.抑制算法可按如下流程進行:1.估計儀器噪聲的均值和方差2.選取合適的結構元素.采用100～150寬度的余弦結構,結構高度為儀器噪聲方差的1/2.3.采用交替混合濾波器對含有儀器噪聲的地磁日變信號進行形態濾波.得到形態濾波日變曲線4.用濾波前的數據與濾波的數據做差,得到差值曲線5.對差值曲線再次進行形態濾波,并將濾波結果回填至步驟3中所產生的形態濾波日變曲線中.在通過諧波擬合的地磁日變信號中加入儀器噪聲,如圖10所示,通過抑制算法,對含有儀器噪聲的地磁日變曲線進行數學形態濾波后,信號信噪比從42.56dB提高到59.75dB.誤差曲線的方差提高到由0.05提高到0.01,考慮到儀器噪聲除去類高斯分布的高頻干擾外,由于環境溫度等外部條件的影響,會有一定量的直流成分的干擾.這也會對通過干擾抑制結果產生一定的影響.3.2形態濾波控制脈沖范圍我們通過實驗,獲得一輛重量約為1噸的汽車以時速20km在離GM-4探頭正前方10m處通過的Z分量干擾曲線,如圖11,從紅色虛線區域可以看出,這個由車輛開行所產生的脈沖干擾的幅度大約在33.2nT,持續時間從第300s到第335s,約35s左右.通過選取長度為35左右,高度0.01的三角元素進行數字形態濾波運算,對脈沖信號進行濾波,圖12為濾波抑制該車輛干擾的效果圖.圖中黑色虛線為受車輛干擾的地磁信號,紅色實線為通過形態濾波器抑制后的地磁信號.4數學形態學在地磁干擾抑制中的應用本文首次將數學形態濾波算法應用在地磁干擾抑制應用中.針對地磁觀測信號中容易受到儀器噪聲和脈沖干擾等環境干擾的影響,通過數值仿真,研究了扁平、余弦、三角等三種不同結構元素的幅值和寬度對于構造形態濾波器抑制此類噪聲的影響,可以得到如下結論:(1)采用數字形態濾波器,通過選取合適的結構元素,可以有效的去除儀器噪聲、脈沖干擾等地磁日變觀測數據中的干擾.(2)針對不同的干擾,需要選取與之相適應的最優結構元素,一般需要根據地磁信號與干擾的具體特征,通過數值仿真來確定.(3)對于隨機背景噪聲等儀器干擾噪聲,應優先選取余弦元素的交替混合濾波器;對于脈沖干擾噪聲,應優先選取三角元素的交替混合濾波器.(4)選取的結構元素的寬度應根據信號采樣頻率的不同而區別對待,對于1Hz采樣的地磁信號,結構元素的長度應控制在100～150之間能夠獲得較優的濾波特性.數學形態學作為一種非線性濾波方法,可以利用靈活的結構元素對信號進行匹配,從而達到抑制相關噪聲,提取信號特征的目的.它是一種完全自適應的,無需人工干預的濾波算法.運算實現簡單,計算負載小,非常適合在儀器內部實現.通過選擇合適的形態結構元素對地磁秒數據處理,其在進行干擾抑制過程中僅對噪聲所在的頻譜段敏感,對于隨機的瞬時脈沖干擾,由于脈沖干擾頻帶分布較廣,應該首先通過算法識別出脈沖干擾的起止時間長度,然后通過形態濾波器對特定長度的脈沖進行形態匹配后再進行抑制處理,最后對處理后的干擾時間段數據進行回填.從而能最大限度避免影響其他未受干擾時間段的地磁觀測數據受到抑制,削弱其頻譜特征.本文對數學形態學在地磁干擾抑制中的應用提出了一些仿真模