李于波，劉祜，段書新，汪碩

（核工業(yè)北京地質(zhì)研究院 中核集團鈾資源勘查與評價技術(shù)重點實驗室，北京 100029）

鹿井地區(qū)位于湖南省郴州市汝城縣境內(nèi)，處于湖南、江西、廣東三省交界處，具有鐵、鎢、銅、稀土、鈾等多種礦產(chǎn)資源。為查明鹿井鈾礦田深部成礦要素，開展了音頻大地電磁（Audio Magnetotelluric Method，簡 稱AMT）測量。但由于工作區(qū)內(nèi)電磁干擾嚴重，部分采集到的時間序列上疊加了與地質(zhì)體無關(guān)的噪音信號，需進行去噪處理以便反演處理得到地下的真實電性結(jié)構(gòu)。

針對音頻大地電磁去噪，前人從時間域或頻率域的角度先后開展了互功率譜法、小波分析法、遠參考技術(shù)、小波自適應閾值、改進閾值的TI 小波、時間序列依賴關(guān)系等去噪技術(shù)［1-7］。上述系列去噪方法取得了一定的應用效果，但亦存在若干問題，如：互功率譜法只能去除不相關(guān)噪聲、小波變換參數(shù)選取不當會導致去噪效果適得其反等。為實現(xiàn)對大量實測數(shù)據(jù)的快速去噪，筆者以鹿井地區(qū)多個測點的時間序列為研究對象，探討了組合廣義形態(tài)濾波方法在該區(qū)AMT 去噪中的應用。

1 鹿井地區(qū)AMT 工程部署及電磁干擾特征

1.1 AMT 工程部署

為查明鹿井鈾礦田深部地質(zhì)情況，在該區(qū)布設(shè)了4 條AMT 探測剖面（圖1）。野外利用V8（主機）結(jié)合AMTC30（磁探頭）的方式開展數(shù)據(jù)采集，獲取到TS2（高頻）、TS3（中頻）、TS4（低頻）3 個時間序列文件，對應的采樣率分別為24 000、2 400、150 Hz。

圖1 鹿井地區(qū)AMT 測線布置圖Fig.1 Layout of AMT survey line in Lujing area

1.2 AMT 電磁干擾特征

據(jù)現(xiàn)場實際踏勘，測區(qū)內(nèi)有高速、省道、鄉(xiāng)道等公路貫穿，同時分布了較多的廣播電臺、高壓線、發(fā)射塔等，預計會對AMT 實測信號產(chǎn)生較多干擾。對實測的時間序列進行分析，證實了上述干擾的存在。

圖2a 中，Ex電場信號由于受到工頻噪聲干擾，其背景值明顯高于其他測點的背景水平。圖2b 中，Ey電場信 號在180 和400 采 樣點處出現(xiàn)明顯的尖波，為三角波噪音在時間序列上的反映。圖2c 中，Hx磁場信號在30 和330 采樣點處亦出現(xiàn)2 個格值超過150 000 的脈沖。圖2d中，Hx磁場信號出現(xiàn)規(guī)律性的似充放電三角波噪聲干擾。為獲取地下真實的電性結(jié)構(gòu)，對上述含噪音數(shù)據(jù)進行去噪處理就顯得尤為重要。

圖2 鹿井地區(qū)采集電場信號受多種噪聲干擾圖Fig.2 The acquired signal interfered by various noises in Lujing area

2 組合廣義形態(tài)濾波原理及其在AMT去噪中的具體實現(xiàn)

形態(tài)濾波是由法國數(shù)學家Matheron G和Serra J于1964年共同創(chuàng)立的信號分析方法，其核心思想是建立適當?shù)慕Y(jié)構(gòu)元素，讓待處理的序列（或集合）依次穿過結(jié)構(gòu)元素，從而提取到該序列（或集合）的數(shù)學形態(tài)。在數(shù)據(jù)分析過程中，若同時使用不同的結(jié)構(gòu)元素進行組合處理，即為組合廣義形態(tài)濾波。

在AMT 時間序列處理過程中，為克服時間序列的基線漂移，采用正負結(jié)構(gòu)元素串行的方法構(gòu)建組合廣義形態(tài)濾波器［8］，其流程如圖3 所示。以Ex 電場信號為例，將信號Input_Ex 依次通過ΨGOC(GCO)(s1，s2)和ΨGOC(GCO)(-s1，-s2)兩個基本濾波單元（其中S1，S2 為結(jié)構(gòu)元素），即得到Ex 電場中所包含的噪聲序列Exlb(n)。

圖3 組合廣義形態(tài)濾波流程圖Fig.3 Flow chart of combined generalized morphology filtering

圖（3）中，ΨGOC(GCO)(s1，s2) 濾波單元的定義為

其中GOC(f(n))、GCO(f(n))分別為廣義形態(tài)開-閉濾波器和廣義形態(tài)閉-開濾波器，“°”為形態(tài)開運算，“?”為形態(tài)閉運算，其定義如式（2）、（3）所示：

式（2）和式（3）中的Θ 為腐蝕運算，其目的是減少待處理序列的峰值，即剔除掉不平滑的峰值數(shù)據(jù)，縮小峰域和拓寬谷域，使目標序列的值收縮；⊕為膨脹運算，其目的是拓寬峰域和縮小谷域，使目標序列的值膨脹，兩運算的公式分別如式（4）、（5）所示：

在本次AMT 去噪處理中，將S1 和S2 分別設(shè)定為5 點圓盤型結(jié)構(gòu)元素（圖4a）和5 點拋物線型結(jié)構(gòu)元素（圖4b），兩者的函數(shù)表達式分別如式（6）、（7）所示：

圖4 五點圓盤型（a）和五點拋物線型（b）結(jié)構(gòu)元素示意圖Fig.4 Diagram of structuring elements of disc type（a）and parabolic type（b）

據(jù)圖3 獲取到噪聲序列Exlb(n)后，將其從Ex電場信號中剔除，即得到重構(gòu)后的時間序列信號Exdn(n)。

3 實測數(shù)據(jù)的去噪應用

為實現(xiàn)對大量實測數(shù)據(jù)的快速去噪，筆者以Python 編程語言為平臺，按照圖5 所示流程實現(xiàn)了鹿井地區(qū)AMT 實測數(shù)據(jù)的組合廣義形態(tài)濾波去噪。整個流程的核心是對4個電磁場分量做正、負結(jié)構(gòu)元素的廣義形態(tài)濾波處理，并重構(gòu)TS 文件。其中TS2 時間序列的結(jié)構(gòu)元素參數(shù)L2=3，k2=0.5；TS3 的結(jié)構(gòu)元素參數(shù)L3=5，k3=0.4；TS4 的結(jié)構(gòu)元素參數(shù)L4=2，k4=0.3。

圖5 數(shù)據(jù)處理及組合廣義形態(tài)濾波去噪流程圖Fig.5 Flow chart of data processing and filtering denoising by generalized morphology combination

以line2-05、line2-07、line2-27 測點為例，重構(gòu)后的時間序列（圖6、7、8）分別去除了階躍噪聲、三角波、似充放電三角波的干擾，表明組合廣義形態(tài)濾波對該類電磁干擾具有較好的去噪效果。

圖6 鹿井地區(qū)采集的含噪時間序列去噪前、去噪后信號對比圖（line2-05）Fig.6 The comparison of original series and denoised series in Lujing area（line2-05）

圖7 鹿井地區(qū)采集的含噪時間序列去噪前、去噪后信號對比圖（line2-07）Fig.7 The comparison of original series and denoised series in Lujing area（line2-07）

圖8 鹿井地區(qū)采集的含噪時間序列去噪前、去噪后信號對比圖（line2-27）Fig.8 The comparison of original series and denoised series in Lujing area（line2-27）

上述去噪效果在各個測點的功率譜密度上亦有所表現(xiàn)。以Line2-07 測點為例，圖9a 展示了TS2 序列中各分量的功率譜密度，曲線在50 Hz 及其倍頻處出現(xiàn)了明顯的峰值情況，表明該點電磁場信號受到了工業(yè)諧波干擾；同時，低頻的功率譜密度小、高頻的功率譜密度大，推測是三角波等大格值噪聲淹沒了低頻有用信號所致。對其進行組合廣義形態(tài)濾波處理，得到重構(gòu)后的時間序列功率譜密度曲線如圖9b 所示。可以看出，去噪后的時間序列在150、250、350、450 Hz 等頻率處的工業(yè)噪聲得到了壓制，且低頻1～300 Hz 范圍內(nèi)的功率譜密度得到了加強，Ex、Ey、HxHy 四道信號的功率譜密度均較為穩(wěn)定，分別保持在1 000、100、10 000、1 000 w/Hz 左右。

圖9 Line2-07 測點TS2、TS3 序列去噪前、后功率譜密度對比圖Fig.9 The comparison of original and denoised TS2 and TS3 PSD of Line 2-07 point

相比較而言，該測點處TS3 時間序列受噪聲影響較小，經(jīng)組合廣義形態(tài)濾波處理前后的功率譜密度曲線變化不大。

圖10 為測點Line2-05 去噪前后的視電阻率曲線對比。由圖可知：在小于1 Hz 的頻率范圍，rhoxy 和rhoyx 經(jīng)廣義形態(tài)濾波去噪后變得更加光滑、穩(wěn)定。針對去噪后rhoxy 曲線在10 Hz 左右出現(xiàn)的跳變，推測是由于形態(tài)濾波時將部分有用信號進行了剔除所致，需要在后續(xù)數(shù)據(jù)處理過程中修改濾波參數(shù)，以減少有效信號的損失，進一步提高廣義形態(tài)濾波效果。

圖10 鹿井Line2-05 測點去噪前、去噪后視電阻率對比圖Fig.10 The comparison of original and denoised apparent resistivity of Line 2-05 point

圖11為測點Line2-07去噪前后的視電阻率曲線對比。由圖可知：去噪后的視電阻率曲線更平滑穩(wěn)定，該現(xiàn)象在rhoxy曲線10 Hz以低頻率表現(xiàn)的更為明顯。整體來看，去噪后50 Hz處的干擾依舊存在，后續(xù)可采用FIR數(shù)字濾波、傅里葉逆變換［9］、希爾伯特-黃變換［10］等方法進行進一步的去除。

圖11 Line2-07 測點去噪前、去噪后視電阻率對比Fig.11 The comparison of original and denoised apparent resistivity of point line 2-07

圖12 為測點Line2-27 去噪前后的視電阻率曲線對比?？梢钥吹?，原始視電阻率曲線在3 Hz 以后出現(xiàn)了上下波動的不穩(wěn)定現(xiàn)象，去噪后的視電阻率曲線則變得連續(xù)、穩(wěn)定。

圖12 鹿井Line2-27 測點去噪前、去噪后視電阻率對比Fig.12 The comparison of original and denoised apparent resistivity of point Line 2-27

4 結(jié)論

通過開展組合廣義形態(tài)濾波在鹿井地區(qū)音頻大地電磁去噪中的應用研究，取得以下成果：

1）采用圓盤型結(jié)構(gòu)元素和拋物線型結(jié)構(gòu)元素構(gòu)建了組合廣義形態(tài)濾波器，在處理鹿井地區(qū)AMT實測數(shù)據(jù)時分別賦以不同的結(jié)構(gòu)元素參數(shù)K、L，對不同頻率采集的時間序列取得了較好的去噪效果。

2）數(shù)學形態(tài)濾波方法在做濾波處理后能提取到噪聲的波形形態(tài)，對于三角波、似充放電三角波、方波、階躍噪聲等波形形態(tài)明顯的電磁干擾，具有較好的去噪效果。

3）將組合廣義形態(tài)濾波法應用于鹿井地區(qū)AMT實測數(shù)據(jù)，能便捷高效地完成時間域TSn時間序列的去噪處理，為AMT數(shù)據(jù)去噪提供了一種新的思路。