葛雙超, 李 斌
(1.中北大學 儀器與電子學院,太原 300051;2.太原衛星發射中心技術部,太原 300051)
我國礦產資源種類繁多,為我國經濟發展提供了堅實的能源保障,但這些礦產資源的勘探程序較低,其中大部分資源的存儲情況未得到準確探測,同時淺部資源開發殆盡,礦產資源日益短缺,為國家的安全發展帶來隱患,因此亟需加強深部復雜礦藏勘探能力,提高礦產資源勘查的深度和效率。
在各種地球物理勘探方法中,電磁勘探是目前應用最多的礦產資源勘探手段之一。大地電磁法(MT)通過在地表測量正交電磁場分量的變化情況,來研究地下介質的電性結構,是被動源電磁勘探的一個重要分支。MT法無需大功率供電設備且探測深度大,在礦產資源勘查、地震前兆監測、環境保護及地下水資源勘查等方面得到廣泛應用。但大地電磁測深要求以平面電磁波作為有效場源,場源噪聲、地質噪聲以及人文噪聲等所有不滿足平面電磁波的電磁信號,都會對MT法觀測數據造成影響[1]。尤其是隨著社會經濟不斷發展,高壓輸電網、通訊電纜、地下金屬管網、廣播電臺、信號塔及各種交通工具等帶來的電磁干擾,嚴重污染了大地電磁實測資料,極大地影響了地下物性結構和電性分析的可解釋性及數據本身的可靠性,導致大地電磁阻抗估算偏差過大,測得的視電阻率曲線失真。而實際工作中又無法避開各類干擾源,故開發有效的大地電磁數據信噪分離方法,提高大地電磁測深數據質量,始終是大地電磁數據采集與處理的核心問題[2]。
筆者首先分析了MT觀測資料中的常見人文噪聲及特點,然后對國內、外近幾年發展起來的重要MT法數據處理方法進行了綜述。主要對Hilbert-Huang 變換、小波分析、統計分析和稀疏分解法等方法的原理、特點及其在MT數據處理中的發展和應用進行了討論,最后對MT觀測數據信噪分離方法做出總結與展望。
MT法的觀測對象是瞬息萬變的天然大地電磁信號,頻帶范圍約為104Hz~10-4Hz。MT觀測數據由有效信號和干擾噪聲兩部分構成,其中前者為天然大地電磁信號,由太陽風與地球磁層、電離層之間復雜的相互作用,以及雷電等地球外部不同頻率段的電磁場相疊加混合形成,屬于非平穩信號,具有很強的隨機性;而干擾噪聲主要來自工頻干擾、游散電流、電設備的開關、車輛噪聲、風擾、水流等,通常具有較強的規律性。MT觀測系統和各類強干擾源示意圖見圖1。

圖1 MT觀測系統和干擾噪聲源
工頻干擾來自觀測點附近的高壓輸電線,主要體現在電通道,兩個正交電通道的工頻分量具有較好的相關性,磁道中偶爾也會出現工頻成分。工頻干擾時域波形(圖2)表現為等振幅類正弦波干擾,頻域表現為工頻干擾基頻(我國通常為50 Hz)及奇次諧波干擾。這類干擾通常強度很大,往往會淹沒有用信號,給后續數據處理帶來很大困難。

圖2 實測MT觀測數據中的工頻干擾
電動設備突然開、關或負荷發生突變時,接地電流導入大地會引起游散電流噪聲干擾,是MT法中的一類常見干擾,其波形圖如圖3所示。該類噪聲通常出現在各種采樣率的電道信號和磁道信號中,在時間序列上通常呈正弦阻尼振蕩,見圖3(a)紅色部分局部放大圖,其幅值是正常有用信號的幾個數量級。

圖3 實測MT觀測數據中的游散電流干擾
電子設備開關瞬間會對MT觀測資料帶來強烈干擾,這類干擾類似方波噪聲,是礦集區內影響強度最大的噪聲之一,其時頻波形如圖4所示。該類噪聲通常出現在中低頻段的電場通道中,且兩正交電通道數據時域波形相關性較好,其幅值一般較大,會淹沒正常大地電磁有用信號,導致阻抗估算出現嚴重偏差。

圖4 電子設備開關造成的方波干擾
電動機調速和閥門控制產生的干擾在觀測數據中表現為不規則的三角波形,一般出現在磁通道中,與磁場信號具有較好的相關性(圖5)。電場信號很少受到三角波噪聲的影響。

圖5 電動機馬達干擾
大型礦山機械進行礦產開采作業時,會產生大范圍高強度電磁干擾,噪聲波形類充放電三角波。該類噪聲強度較大,觀測數據時域波形有明顯跳變,各通道具有較好的相關性(圖6)。由圖6可以看出,觀測數據中除車輛干擾外還混有工頻干擾。

圖6 大型礦山機械車輛干擾
除上述幾種常見干擾外,MT觀測資料中還包含周圍人員或動物活動、風擾等一些隨機噪聲,在觀測數據時域波形圖上體現為短時不規則噪聲(圖7)。

圖7 短時隨機干擾
進行大地電磁測深時,通常采集正交電磁場的時域信息,然后進行傅里葉變換并進行阻抗估計,由于傅里葉變換適用于確定性的平穩信號,所以采用這種方法進行MT信號處理時,通常假設MT信號具有高斯分布、地質模型是最小相位系統、電磁場是平面電磁波等。王書明等[3]通過對國內多個不同地區的MT觀測資料進行統計分析,證實MT信號一般具有非高斯、非線性和非最小相位性。因此,直接采用傅里葉變換來分析和處理大地電磁信號具有一定的局限性。為提高阻抗估計準確度,地球物理學家不斷改進和優化MT資料處理和參數計算方法,比如時域處理中的Hilbert-Huang變換、小波分析、統計分析、經驗模態分解等。
Hilbert-Huang變換是20世紀90年代中期發展起來的一種處理非線性、非平穩信號的有效方法[4]。利用Hilbert-Huang變換進行MT觀測資料處理時,首先對MT信號進行經驗模態分解(EMD),獲得各階固有模態函數;然后根據噪聲的不同來源和特征,對各分量進行提取和重構,從而提高大地電磁信號信噪比。
蔡建華等[5]研究了Hilbert-Huang變換在大地電磁信號噪聲識別、功率譜估計和阻抗估算等方面的應用,開發了相應的數據處理算法對大地電磁信號中幾種典型干擾噪聲進行了處理。結果表明,所述方法對基線漂移、高頻噪聲和工頻干擾有較好的抑制效果。

圖8 小波變換去噪流程
小波變換是一種分析瞬時時變信號的有力工具,通過伸縮和平移等運算功能對函數或信號進行多尺度細化分析,能夠成功處理傅立葉變換不能解決的問題[6]。因此小波變換成為大地電磁信號信噪分離的有效方法之一。
小波變換與傅立葉變換的原理類似,但小波變換利用有限長會衰減的小波基代替無限長的三角函數基。連續小波變換公式為式(1)。
(1)
其中:x(t)為原始信號;ψ(t)為小波基函數;α為尺度因子;τ為平移因子;Q控制小波基函數的伸縮;τ控制小波基函數的平移。不同的小波基函數在處理信號時各有不同,為保證最優信噪分離效果,需要事先優選合適的小波基函數。
小波變換結果Wt(α,τ)為分解系數,是信號在所選擇的小波函數空間的投影,與α成反比。
由于規則信號的能量集中在少數小波系數上,而白噪聲在小波變換域上仍然分散在大量小波系數之上。利用非線性閾值對小波變換系數進行處理,將較小的小波系數置零并保留較大的小波系數(硬閾值法),或者對較大的小波系數向零收縮(軟閾值法)。然后通過逆變換對信號進行重構,進而實現信噪分離。
小波變換信噪分離流程見圖8。
小波變換是一個多分辨率分析過程,以三層小波分解為例,其分析樹結構如圖9所示,分解結果如下:

圖9 多分辨分析樹結構圖
x=sss+dss+ds+d
(2)
隨著分解層數的增加細節信號(即高頻信號、隨機信號)會逐漸突出。高層分解的小波系數對應低頻信號分量,通常情況下低頻部分主要由規則信號或干擾構成。因此分解層次越高,去掉的低頻成分越多,去噪效果越明顯,但失真度也增大。所以在實際應用中,應綜合考慮信號的實際變化趨勢及其理論變化規律來合理確定小波分解層數,在保證有效分離不同頻率的信號成分的同時,又要避免過度分解導致信號內在變化規律和趨勢降低。Daniel Trad等[7]采用離散小波變換和Robust阻抗估計相結合的方法處理MT觀測數據,當干擾成分在時域和頻域均表現為脈沖噪聲時,該方法能夠獲得較好的觀測結果。
凌振寶等[8]針對礦集區人文噪聲中的方波、基線漂移的長周期干擾,以及周邊大型機械突然開關產生的高頻干擾對MT觀測資料的影響,基于多分辨率分析算法,先采用db3小波基對觀測信號進行10層分解來去除頻率低于0.2 Hz長周期噪聲;然后利用經驗貝葉斯估計優化小波閾值函數去除頻率高于0.375 Hz的高頻干擾。實測結果證實所提方法可以有效抑制近源效應。Cai Jianhua等[9]為克服小波閾值選擇不合理對MT觀測數據處理的誤差問題,提出了一種基于中點優化算法和多層分解Stein無偏似然估計的自適應小波閾值去噪技術,并用此方法對MT觀測資料中的強脈沖干擾進行了抑制。曹小玲等[10]結合小波分析和盲源分離的相關理論,提出了一種改進的獨立分量分析去噪方法。模擬仿真實驗證明該方法相比傳統小波閾值去噪方法性能更優。
局部均值分解(LMD)是一種新的自適應非平穩信號處理方法。該方法自適應地將一個復雜非平穩信號,分解為若干個瞬時頻率具有物理意義的乘積函數(Product function,PF),其中每一個PF分量由一個包絡信號和一個純調頻信號相乘獲得[11]。前者對應該PF分量的瞬時幅值,后者可用于計算PF分量的瞬時頻率。該算法由兩層循環構成,內循環用于求解某PF分量的包絡和調頻信號,外循環通過迭代處理,將所有PF分量分離出來,最終可獲得原始信號的時頻分布(圖10)。

圖10 LMD算法流程圖
李晉等[12]根據大地電磁信號的特征,結合局域均值分解法(LMD)和小波分析法,提出了一種基于局域均值分解和小波閾值的大地電磁噪聲壓制方法。首先利用局域均值分解將大地電磁信號自適應地分解為一系列PF分量之和;然后保留PF1分量,僅對其余各階PF分量選擇合適的軟閾值進行小波閾值去噪;最后疊加重構獲得大地電磁有用信號。計算機模擬仿真結果表明,所述方法可用于壓制礦集區典型大尺度方波和充放電三角波干擾。Zhang Gang等[13]為補償相關檢測技術在處理MT觀測數據時磁場數據相關性問題,利用磁通道和電通道的相干系數分析法,對龍門山地區的遠參考MT數據進行預處理,然后再用Robust法進行阻抗張力估計。實驗結果證明,在信噪比較低情況下相干分析可以有效改進阻抗估計效果。
傅里葉變換和小波變換都屬于線性濾波,形態濾波屬于非線性濾波器,可以有效提取信號的邊緣輪廓以及形狀特征。
數學形態學由腐蝕和膨脹兩個基本運算單元構成。對原始離散信號序列x(n),其中n=1、2、…、N;定義結構元素g(m),其中m=1、2、…、M,且M≤N。x(n)關于g(m)的腐蝕和膨脹操作為:
腐蝕(xΘg)(n)=min[x(n+m)-g(m)],m∈0,1,…,M-1
膨脹(x⊕g)(n)=max[x(n+m)+g(m)],m∈0,1,…,M-1
(3)
其中,min和max分別為求極小和極大值函數。
腐蝕和膨脹運算等價于離散函數在滑動濾波窗(g(m))內的最小值和最大值濾波。根據腐蝕和膨脹運算,x(n)關于g(m)的開運算和閉運算如下:
開運算(x°g)(n)=(xΘg)⊕g
閉運算(x·g)(n)=(x⊕g)Θg
(4)
形態學開閉運算是一對對偶變換,其中形態開運算可以抑制信號中的正脈沖噪聲,而形態閉運算可以抑制信號中的負脈沖噪聲[14]。
形態學濾波的效果由變換形式和結構元素共同決定。由于形態開、閉運算存在統計偏倚現象,導致單獨使用時濾波效果不夠好,因此通常采用開-閉或者閉-開組合濾波方法。結構元素的選取包括結構元素的形狀、寬度和高度等[15]。為達到理想的信噪分離效果,進行形態學濾波時,需要根據信號的幾何特征預先定義合適的結構元素對信號進行匹配。
李晉等[16]基于數學形態學的思想,研究了傳統形態濾波和廣義形態濾波的大地電磁強干擾分離方法,并在形態濾波的基礎上,研究了Top-hat變換、中值濾波和信號子空間增強的大地電磁二次信噪分離方法。仿真和實驗結果表明,基于數學形態學的大地電磁強干擾分離方法,對大地電磁強干擾中的大尺度干擾和基線漂移有較好的抑制作用,能夠有效改善MT數據質量。
進行MT阻抗估計時,當噪聲滿足獨立Gauss分布時最小二乘估計為最優無偏估計[17]。為了得到更可靠的阻抗估計,人們又相繼提出了更優的Robust阻抗估計法和遠參考法。遠參考大地電磁法可消除非相關噪聲,其基本思想就是在同步的情況下,利用與測點相隔一定距離的遠參考點的磁信號作為測點處的磁信號估算張量阻抗。當遠參考點和測點間距達到一定程度,且兩點位的噪聲特性非相關的情況下,遠參考處理能夠有效提高張量阻抗計算精度。這兩種方法在目前的大地電磁阻抗估計法中相對比較成熟,得到了廣泛地應用[18]。以上幾種阻抗估計方法的使用前提均是假設其殘差服從高斯分布,然而在實際電磁法勘探中,這一假設并不一定成立。
Shalivahan等[19]采用相干加權估計、抗差M估計方法和視電阻率加權估計處理10-3Hz~103Hz的寬頻MT觀測數據,結果表明三種方法聯合處理能夠有效抑制死區頻帶內的異常干擾。Antoine Saucier等[20]針對諧波干擾,在Nyman-Gaiser估計的基礎上改進諧波干擾模型,并利用頻域迭代法估計諧波基頻和各次諧波頻率,然后利用最小二乘法估計各諧波分量的振幅和相位。相比傳統最小二乘法,該方法的數據處理效率更高。但是該方法要求諧波基頻預估值與真實值相差不大且基頻波動范圍較小。Kappler等[21]為消除MT觀測資料中的脈沖和階躍噪聲,利用站間數據不一致性識別干擾成分,并利用維納濾波進行誤差補償。該方法在多臺站聯合觀測中表現出良好的應用效果。Ge Shuangchao等[22]利用全通濾波器減去梳狀濾波器方式,設計了一種全通帶諧波陷波器,可以一次性濾除觀測數據頻帶內的工頻干擾基頻及所有奇次諧波干擾。該濾波器適合作MT觀測系統的前端濾波。李晉等[23-24]為了分離出微弱的大地電磁有用信號,提出基于遞歸分析和聚類的大地電磁信噪辨識及分離方法。利用該方法對大地電磁信號和強于擾的本質特征進行了分析,并研究了基于信噪辨識的礦集區大地電磁噪聲壓制方法。仿真和實測實驗表明該方法,對于低頻段的緩變化信息識別具有較好的效果。燕歡等[25]針對礦集區典型的強干擾類型,引入子空間增強方法對大地電磁觀測數據中的大尺度強干擾進行去噪處理。針對礦集區MT觀測常見的大尺度方波和充放電三角波干擾,通過構造協方差矩陣和特征值分解等方法將噪聲子空間和信號子空間分離,并重構原始信號來實現信噪分離。
Mallat等[26-27]在小波分析的基礎上,提出信號在一定條件下可以用超完備字典進行稀疏表示,并在此思想上引入匹配追蹤算法。Candès等[28]提出壓縮感知技術后,信號的稀疏表示作為壓縮感知的重要組成部分,得到了系統性的研究和發展。利用冗余字典對信號進行稀疏分解,將信號的特征信息在變換域直觀呈現出來,然后采用合適的手段進行信噪分離,使得稀疏分解方法成為一種新的信號去噪技術[29],在波達方向估計[30]、音頻信號中風聲[31]和隨機噪聲[32]等干擾處理、圖像去噪[33]、數據采集系統的重尾脈沖干擾抑制[34]等多領域得到了廣泛應用。
壓縮感知和信號的稀疏分解在地球物理領域的應用,主要集中在地震勘探和合成孔徑雷達成像等方面。Liu Wei等[35]綜合利用壓縮感知和小波分析處理地震數據,通過梯度投影稀疏重構算法,將噪聲分解問題轉換為L1范數正則匹配問題,進而對隨機噪聲進行了有效抑制;Liu Sicong等[36]針對車載通信中的窄帶和脈沖噪聲,搭建了基于壓縮感知的時頻測量正交頻分復用框架,利用先驗稀疏度自適應匹配追蹤對上述兩種非高斯噪聲進行重構,提高了車載通信的效率;祁銳[37]針對塊稀疏信號重構問題,提出基于BAOMP 算法的分布式壓縮感知重構算法,該算法的原子選擇準則類似于正交匹配追蹤算法,不同之處在于每次選擇固定個數的塊指標,從而減少迭代次數,提高了算法的重構性能;Mohammad Amir Nazari Siahsar等[38]基于矩陣低秩和稀疏分解的方法處理地震數據中的非相干隨機噪聲,通過對原始觀測數據進行稀疏時頻表示將地震有效信號分解為低秩和稀疏分量,并利用匹配追蹤算法進行降噪處理,仿真和實測結果驗證了算法的有效性。
在電磁勘探領域,已經有學者嘗試采用稀疏表示或匹配追蹤技術對MT觀測資料進行處理,湯井田等[39]構造了諧波、方波和尖沖干擾原子庫,利用子空間追蹤和正交匹配追蹤相結合的方法處理MT觀測數據中的強干擾成分。由于原子庫冗余度較高,該算法的缺點是耗時較長,算法效率不夠高。為了改進匹配追蹤在大地電磁強干擾分離中的計算效率,李晉等[40]引入匹配追蹤和遺傳算法處理MT觀測資料中的大尺度強干擾時間序列。實驗結果證實,該方法能從過完備原子庫中快速、自適應地搜尋最優匹配原子重構大地電磁有用信號。但是在重構過程中會造成有用信號的損失,另外該算法對于稀疏度的自動獲取的研究還不夠完善。李廣等[41]針對音頻大地電磁觀測資料中的近場源干擾問題,首先利用形態濾波處理觀測數據中的大尺度強干擾,然后利用壓縮感知重構算法分離尖沖噪聲。該方法在強干擾區應用效果良好,但是當干擾信號強度較弱時,算法性能會隨干擾強度衰減而降低。另外為提高算法效率和信噪分離效果,亟需構建更加合理的冗余字典庫和開發更智能的匹配算法。湯井田等[42]通過分別設計與方波信號、脈沖信號以及諧波信號匹配,而對有用的大地電磁信號不敏感的稀疏表示字典,利用IOMP算法進行重構,有效地抑制了非平穩信號中的方波噪聲和脈沖噪聲。湯井田等[42]在分析噪聲的規律與特征基礎上,提出一種新的AMT數據處理方法。利用移不變稀疏編碼字典學習方法從觀測數據中自主學習人文噪聲的特征結構,構建冗余字典,然后利用學習到的冗余字典對觀測數據進行信噪分離,最后用去噪后的數據替換原始數據計算視電阻率和相位曲線。野外試驗結果表明,所述方法能夠快速、準確地分離出AMT信號中的人文干擾,保留有用信號,顯著改善AMT數據質量。葉華等[43]提出了基于稀疏表示與粒子群優化算法的非平穩信號去噪方法,以沖擊原子作為稀疏表示基,構建了僅對人文噪聲敏感的冗余字典;采用粒子群優化算法選擇最優原子,提高了匹配追蹤算法的效率。
Hilbert-Huang變換關鍵步驟是EMD分解,該過程需進行循環篩選,算法效率低,且易出現模態混疊現象。由于MT法觀測數據量較大,故該方法難以滿足實時性需求。所以為提高工作效率,需要對該方法進行優化處理。
小波變換是最經典的非平穩信號處理方法,在諸多領域取得較好的效果,但小波濾波過分依賴小波基函數的選取,如果小波基函數選擇不合適或閾值設置不合理,可能會損失部分有用信號,且小波基函數選定后,在分解與重構過程中將無法改變,缺乏自適應性。
經驗模態分解類的信號處理方法,對于諧波類的干擾具有良好的抑制效果,能夠較好地處理工頻電力諧波等干擾,但對于方波、充放電噪聲以及脈沖噪聲不甚理想,且存在模態混疊和端點效應問題。形態濾波法能夠壓制大尺度的強人文噪聲,但是對于脈沖噪聲效果不佳,經過形態濾波處理后,常常殘留大量的脈沖噪聲,此外,形態濾波容易損傷低頻有用信號。然而,礦集區的噪聲干擾情況復雜,當采集的大地電磁原始數據中的噪聲干擾形態特征不規則時,信號與噪聲出現混淆,信噪辨識的精度將有所下降。
信號稀疏表示常用算法屬于貪婪算法,由于原子庫冗余度較高,算法耗時較長,算法效率不夠高、容易陷入局部最優原子。另外,由于實測數據中的噪聲形態復雜多樣,預先設定的冗余字典無法與各類噪聲完全匹配,從而不可避免的割裂噪聲的完整結構,并損失部分有用信號。
近年來,仿生智能優化算法在圖像降噪處理、機械振動噪聲識別等多個研究領域得到廣泛關注。該類算法是研究者通過分析生物集群活動機理、覓食現象等而提出的用以解決復雜問題的新方法,代表性算法包括蟻群算法、螢火蟲算法、蝙蝠算法等。基于智能算法的信號特征提取、信噪分離等,成為各類信號處理領域的重要發展趨勢。通過系統分析和建模,構建各類不同噪聲模型的超完備冗余字典庫,基于智能算法實現干擾噪聲的快速自適應提取,是提高MT觀測資料信噪比的一種新思路。在進行此類問題研究時,需要構造合理的尋優模型、研究合理的參數調整判據,保證在干擾信號精確重構的同時最大限度保留有效MT信號,最終實現強干擾下微弱MT信號的有效提取。