李于波,劉祜,段書新,汪碩
(核工業(yè)北京地質(zhì)研究院 中核集團(tuán)鈾資源勘查與評(píng)價(jià)技術(shù)重點(diǎn)實(shí)驗(yàn)室,北京 100029)
鹿井地區(qū)位于湖南省郴州市汝城縣境內(nèi),處于湖南、江西、廣東三省交界處,具有鐵、鎢、銅、稀土、鈾等多種礦產(chǎn)資源。為查明鹿井鈾礦田深部成礦要素,開展了音頻大地電磁(Audio Magnetotelluric Method,簡(jiǎn) 稱AMT)測(cè)量。但由于工作區(qū)內(nèi)電磁干擾嚴(yán)重,部分采集到的時(shí)間序列上疊加了與地質(zhì)體無(wú)關(guān)的噪音信號(hào),需進(jìn)行去噪處理以便反演處理得到地下的真實(shí)電性結(jié)構(gòu)。
針對(duì)音頻大地電磁去噪,前人從時(shí)間域或頻率域的角度先后開展了互功率譜法、小波分析法、遠(yuǎn)參考技術(shù)、小波自適應(yīng)閾值、改進(jìn)閾值的TI 小波、時(shí)間序列依賴關(guān)系等去噪技術(shù)[1-7]。上述系列去噪方法取得了一定的應(yīng)用效果,但亦存在若干問(wèn)題,如:互功率譜法只能去除不相關(guān)噪聲、小波變換參數(shù)選取不當(dāng)會(huì)導(dǎo)致去噪效果適得其反等。為實(shí)現(xiàn)對(duì)大量實(shí)測(cè)數(shù)據(jù)的快速去噪,筆者以鹿井地區(qū)多個(gè)測(cè)點(diǎn)的時(shí)間序列為研究對(duì)象,探討了組合廣義形態(tài)濾波方法在該區(qū)AMT 去噪中的應(yīng)用。
為查明鹿井鈾礦田深部地質(zhì)情況,在該區(qū)布設(shè)了4 條AMT 探測(cè)剖面(圖1)。野外利用V8(主機(jī))結(jié)合AMTC30(磁探頭)的方式開展數(shù)據(jù)采集,獲取到TS2(高頻)、TS3(中頻)、TS4(低頻)3 個(gè)時(shí)間序列文件,對(duì)應(yīng)的采樣率分別為24 000、2 400、150 Hz。
圖1 鹿井地區(qū)AMT 測(cè)線布置圖Fig.1 Layout of AMT survey line in Lujing area
據(jù)現(xiàn)場(chǎng)實(shí)際踏勘,測(cè)區(qū)內(nèi)有高速、省道、鄉(xiāng)道等公路貫穿,同時(shí)分布了較多的廣播電臺(tái)、高壓線、發(fā)射塔等,預(yù)計(jì)會(huì)對(duì)AMT 實(shí)測(cè)信號(hào)產(chǎn)生較多干擾。對(duì)實(shí)測(cè)的時(shí)間序列進(jìn)行分析,證實(shí)了上述干擾的存在。
圖2a 中,Ex電場(chǎng)信號(hào)由于受到工頻噪聲干擾,其背景值明顯高于其他測(cè)點(diǎn)的背景水平。圖2b 中,Ey電場(chǎng)信 號(hào)在180 和400 采 樣點(diǎn)處出現(xiàn)明顯的尖波,為三角波噪音在時(shí)間序列上的反映。圖2c 中,Hx磁場(chǎng)信號(hào)在30 和330 采樣點(diǎn)處亦出現(xiàn)2 個(gè)格值超過(guò)150 000 的脈沖。圖2d中,Hx磁場(chǎng)信號(hào)出現(xiàn)規(guī)律性的似充放電三角波噪聲干擾。為獲取地下真實(shí)的電性結(jié)構(gòu),對(duì)上述含噪音數(shù)據(jù)進(jìn)行去噪處理就顯得尤為重要。
圖2 鹿井地區(qū)采集電場(chǎng)信號(hào)受多種噪聲干擾圖Fig.2 The acquired signal interfered by various noises in Lujing area
形態(tài)濾波是由法國(guó)數(shù)學(xué)家Matheron G和Serra J于1964年共同創(chuàng)立的信號(hào)分析方法,其核心思想是建立適當(dāng)?shù)慕Y(jié)構(gòu)元素,讓待處理的序列(或集合)依次穿過(guò)結(jié)構(gòu)元素,從而提取到該序列(或集合)的數(shù)學(xué)形態(tài)。在數(shù)據(jù)分析過(guò)程中,若同時(shí)使用不同的結(jié)構(gòu)元素進(jìn)行組合處理,即為組合廣義形態(tài)濾波。
在AMT 時(shí)間序列處理過(guò)程中,為克服時(shí)間序列的基線漂移,采用正負(fù)結(jié)構(gòu)元素串行的方法構(gòu)建組合廣義形態(tài)濾波器[8],其流程如圖3 所示。以Ex 電場(chǎng)信號(hào)為例,將信號(hào)Input_Ex 依次通過(guò)ΨGOC(GCO)(s1,s2)和ΨGOC(GCO)(-s1,-s2)兩個(gè)基本濾波單元(其中S1,S2 為結(jié)構(gòu)元素),即得到Ex 電場(chǎng)中所包含的噪聲序列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)開運(yùn)算,“?”為形態(tài)閉運(yùn)算,其定義如式(2)、(3)所示:
式(2)和式(3)中的Θ 為腐蝕運(yùn)算,其目的是減少待處理序列的峰值,即剔除掉不平滑的峰值數(shù)據(jù),縮小峰域和拓寬谷域,使目標(biāo)序列的值收縮;⊕為膨脹運(yùn)算,其目的是拓寬峰域和縮小谷域,使目標(biāo)序列的值膨脹,兩運(yùn)算的公式分別如式(4)、(5)所示:
在本次AMT 去噪處理中,將S1 和S2 分別設(shè)定為5 點(diǎn)圓盤型結(jié)構(gòu)元素(圖4a)和5 點(diǎn)拋物線型結(jié)構(gòu)元素(圖4b),兩者的函數(shù)表達(dá)式分別如式(6)、(7)所示:
圖4 五點(diǎn)圓盤型(a)和五點(diǎn)拋物線型(b)結(jié)構(gòu)元素示意圖Fig.4 Diagram of structuring elements of disc type(a)and parabolic type(b)
據(jù)圖3 獲取到噪聲序列Exlb(n)后,將其從Ex電場(chǎng)信號(hào)中剔除,即得到重構(gòu)后的時(shí)間序列信號(hào)Exdn(n)。
為實(shí)現(xiàn)對(duì)大量實(shí)測(cè)數(shù)據(jù)的快速去噪,筆者以Python 編程語(yǔ)言為平臺(tái),按照?qǐng)D5 所示流程實(shí)現(xiàn)了鹿井地區(qū)AMT 實(shí)測(cè)數(shù)據(jù)的組合廣義形態(tài)濾波去噪。整個(gè)流程的核心是對(duì)4個(gè)電磁場(chǎng)分量做正、負(fù)結(jié)構(gòu)元素的廣義形態(tài)濾波處理,并重構(gòu)TS 文件。其中TS2 時(shí)間序列的結(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 測(cè)點(diǎn)為例,重構(gòu)后的時(shí)間序列(圖6、7、8)分別去除了階躍噪聲、三角波、似充放電三角波的干擾,表明組合廣義形態(tài)濾波對(duì)該類電磁干擾具有較好的去噪效果。
圖6 鹿井地區(qū)采集的含噪時(shí)間序列去噪前、去噪后信號(hào)對(duì)比圖(line2-05)Fig.6 The comparison of original series and denoised series in Lujing area(line2-05)
圖7 鹿井地區(qū)采集的含噪時(shí)間序列去噪前、去噪后信號(hào)對(duì)比圖(line2-07)Fig.7 The comparison of original series and denoised series in Lujing area(line2-07)
圖8 鹿井地區(qū)采集的含噪時(shí)間序列去噪前、去噪后信號(hào)對(duì)比圖(line2-27)Fig.8 The comparison of original series and denoised series in Lujing area(line2-27)
上述去噪效果在各個(gè)測(cè)點(diǎn)的功率譜密度上亦有所表現(xiàn)。以Line2-07 測(cè)點(diǎn)為例,圖9a 展示了TS2 序列中各分量的功率譜密度,曲線在50 Hz 及其倍頻處出現(xiàn)了明顯的峰值情況,表明該點(diǎn)電磁場(chǎng)信號(hào)受到了工業(yè)諧波干擾;同時(shí),低頻的功率譜密度小、高頻的功率譜密度大,推測(cè)是三角波等大格值噪聲淹沒(méi)了低頻有用信號(hào)所致。對(duì)其進(jìn)行組合廣義形態(tài)濾波處理,得到重構(gòu)后的時(shí)間序列功率譜密度曲線如圖9b 所示。可以看出,去噪后的時(shí)間序列在150、250、350、450 Hz 等頻率處的工業(yè)噪聲得到了壓制,且低頻1~300 Hz 范圍內(nèi)的功率譜密度得到了加強(qiáng),Ex、Ey、HxHy 四道信號(hào)的功率譜密度均較為穩(wěn)定,分別保持在1 000、100、10 000、1 000 w/Hz 左右。
圖9 Line2-07 測(cè)點(diǎn)TS2、TS3 序列去噪前、后功率譜密度對(duì)比圖Fig.9 The comparison of original and denoised TS2 and TS3 PSD of Line 2-07 point
相比較而言,該測(cè)點(diǎn)處TS3 時(shí)間序列受噪聲影響較小,經(jīng)組合廣義形態(tài)濾波處理前后的功率譜密度曲線變化不大。
圖10 為測(cè)點(diǎn)Line2-05 去噪前后的視電阻率曲線對(duì)比。由圖可知:在小于1 Hz 的頻率范圍,rhoxy 和rhoyx 經(jīng)廣義形態(tài)濾波去噪后變得更加光滑、穩(wěn)定。針對(duì)去噪后rhoxy 曲線在10 Hz 左右出現(xiàn)的跳變,推測(cè)是由于形態(tài)濾波時(shí)將部分有用信號(hào)進(jìn)行了剔除所致,需要在后續(xù)數(shù)據(jù)處理過(guò)程中修改濾波參數(shù),以減少有效信號(hào)的損失,進(jìn)一步提高廣義形態(tài)濾波效果。
圖10 鹿井Line2-05 測(cè)點(diǎn)去噪前、去噪后視電阻率對(duì)比圖Fig.10 The comparison of original and denoised apparent resistivity of Line 2-05 point
圖11為測(cè)點(diǎn)Line2-07去噪前后的視電阻率曲線對(duì)比。由圖可知:去噪后的視電阻率曲線更平滑穩(wěn)定,該現(xiàn)象在rhoxy曲線10 Hz以低頻率表現(xiàn)的更為明顯。整體來(lái)看,去噪后50 Hz處的干擾依舊存在,后續(xù)可采用FIR數(shù)字濾波、傅里葉逆變換[9]、希爾伯特-黃變換[10]等方法進(jìn)行進(jìn)一步的去除。
圖11 Line2-07 測(cè)點(diǎn)去噪前、去噪后視電阻率對(duì)比Fig.11 The comparison of original and denoised apparent resistivity of point line 2-07
圖12 為測(cè)點(diǎn)Line2-27 去噪前后的視電阻率曲線對(duì)比??梢钥吹剑家曤娮杪是€在3 Hz 以后出現(xiàn)了上下波動(dòng)的不穩(wěn)定現(xiàn)象,去噪后的視電阻率曲線則變得連續(xù)、穩(wěn)定。
圖12 鹿井Line2-27 測(cè)點(diǎn)去噪前、去噪后視電阻率對(duì)比Fig.12 The comparison of original and denoised apparent resistivity of point Line 2-27
通過(guò)開展組合廣義形態(tài)濾波在鹿井地區(qū)音頻大地電磁去噪中的應(yīng)用研究,取得以下成果:
1)采用圓盤型結(jié)構(gòu)元素和拋物線型結(jié)構(gòu)元素構(gòu)建了組合廣義形態(tài)濾波器,在處理鹿井地區(qū)AMT實(shí)測(cè)數(shù)據(jù)時(shí)分別賦以不同的結(jié)構(gòu)元素參數(shù)K、L,對(duì)不同頻率采集的時(shí)間序列取得了較好的去噪效果。
2)數(shù)學(xué)形態(tài)濾波方法在做濾波處理后能提取到噪聲的波形形態(tài),對(duì)于三角波、似充放電三角波、方波、階躍噪聲等波形形態(tài)明顯的電磁干擾,具有較好的去噪效果。
3)將組合廣義形態(tài)濾波法應(yīng)用于鹿井地區(qū)AMT實(shí)測(cè)數(shù)據(jù),能便捷高效地完成時(shí)間域TSn時(shí)間序列的去噪處理,為AMT數(shù)據(jù)去噪提供了一種新的思路。