武 瑩， 陸從德， 杜興忠， 余小東

(1. 成都理工大學(xué) 信息科學(xué)與技術(shù)學(xué)院，成都 610059；2. 成都理工大學(xué) 地球物理學(xué)院，成都 610059；3.中國水電顧問集團 貴陽勘測設(shè)計研究院，貴陽 550081)

0 前言

航空瞬變電磁法是一種快速地球物理勘查方法，因其勘探效率高，成本相對較低，可大面積勘探等優(yōu)勢，而廣泛用于地質(zhì)填圖、礦產(chǎn)資源勘查、油氣勘查，以及水文、工程、環(huán)境監(jiān)測等各個領(lǐng)域[1]。為了避免激發(fā)源的影響，航空瞬變電磁法常常對二次場數(shù)據(jù)進(jìn)行處理和解釋，以獲得地下介質(zhì)分布信息。關(guān)于航空瞬變電磁二次場信號：①由于該信號是一種寬頻帶信號，因此易受到多種噪聲的影響，例如隨機噪聲、天電噪聲和人文噪聲等；②由于二次場信號能量較弱，導(dǎo)致中晚期的信號被噪聲嚴(yán)重污染，使得處理和解釋結(jié)果不可靠或者不可信。為了解決上述問題，設(shè)計出保幅去噪方法就顯得尤為重要[2-3]。

關(guān)于航空瞬變電磁去噪研究，已經(jīng)有許多學(xué)者提出了不同的去噪方法，①最小二乘濾波和頻率域的線性濾波法可以減弱隨機噪聲的影響[4-5]，但是前者在有效信號和噪聲相關(guān)時不能獲得較好的效果，而后者往往存在吉布斯現(xiàn)象，即所謂的邊緣效應(yīng)；②裁剪法、中值濾波法、小波閾值濾波法可以有效地抑制天電噪聲[6-9]，但是裁剪法主要針對疊前數(shù)據(jù)，中值濾波在抑制具有不同時寬的天電噪聲時性能較差(這是因為所取窗口的大小形狀會對濾波效果造成較大影響)，小波閾值去噪法對于幅值變化不大的信號有很好的效果，但是對于幅值變化較大的航空瞬變電磁信號不能獲得較佳的去噪性能；③錐形疊加法以及局部噪聲預(yù)測濾波法、自適應(yīng)濾波方法等基于預(yù)測或遞歸濾波的去噪算法，對于抑制航空瞬變電磁數(shù)據(jù)中的人文噪聲有很好的效果[6,10-12]，但是錐形疊加法主要用于疊前數(shù)據(jù)，預(yù)測濾波算法要求很多的輸入?yún)?shù)，其實現(xiàn)較為復(fù)雜，很難應(yīng)用于工程實踐，自適應(yīng)濾波方法具有普適性，但是其去噪效果有待改進(jìn)。在實際的航空瞬變電磁勘查中，影響有效信號的噪聲種類很多，往往是復(fù)雜的和不可預(yù)測的，而上述各種方法都只能處理單一噪聲，不能同時抑制天電噪聲和人文噪聲。從目前有關(guān)研究文獻(xiàn)可知，天電噪聲和人文噪聲是航空瞬變電磁噪聲的主要來源[3]，為此尋找一種計算簡單又同時可以抑制航空瞬變電磁主要噪聲的去噪方法，就成為了提高處理和解釋精度的關(guān)鍵內(nèi)容之一。

主成分分析法(PCA)能對含噪信號分析，達(dá)到信噪分離的目的。目前已有研究者應(yīng)用PCA法到地球物理數(shù)據(jù)處理中。夏江海[13]在位場資料處理中應(yīng)用奇異值分解有效抑制了隨機噪聲；王權(quán)海等[14]將PCA方法應(yīng)用于地震數(shù)據(jù)處理中，極大提高了地震數(shù)據(jù)的信噪比；Reninger等[15]把奇異值分解作為一種去噪工具,應(yīng)用到航空瞬變電磁數(shù)據(jù)處理中。在文獻(xiàn)[ 15]中，盡管奇異值分解法對天電噪聲和人文噪聲具有良好的去噪效果，但是其在確定信噪分離準(zhǔn)則時需要較為復(fù)雜的條件(例如飛行高度和飛行視頻等)，不易應(yīng)用到工程實踐中。作者同樣使用主成分分析方法對航空瞬變電磁數(shù)據(jù)進(jìn)行分析，但是在確定信噪分離準(zhǔn)則時使用較為簡單的L曲線法，不僅能獲得良好的去噪效果，而且易于嵌入到航空瞬變電磁處理與解釋軟件中，應(yīng)用于工程實踐。

1 基于主成分分析的航空瞬變電磁去噪方法

對于航空瞬變電磁數(shù)據(jù)來講，由于地下電性特征具有一定的分布規(guī)律，所以觀測到的電磁場數(shù)據(jù)隨空間、時間、頻率、極距和發(fā)收距的變化也遵循一定規(guī)律，這種規(guī)律性在數(shù)學(xué)上就是線性相關(guān)性?？梢哉J(rèn)為反映地下電性特征的數(shù)據(jù)是由數(shù)據(jù)矩陣中各列間相關(guān)性很高的部分組成，而數(shù)據(jù)中不相關(guān)的部分則認(rèn)為是各種噪聲和干擾[16]。這就是說，能量較大的主要成分代表地下電性特征，能量較小的成分代表了各種噪聲和干擾。實際上這也是在航空瞬變電磁去噪中應(yīng)用主成分分析法的依據(jù)。

設(shè)XXm×n表示航空瞬變電磁二次場，其中n表示測點，m表示時間道，于是獲得基于主成分分析的航空瞬變電磁去噪方法的計算步驟為：

(1)對XXm×n進(jìn)行歸一化獲得Xm×n。

(2)計算Xm×n的協(xié)方差矩陣Dm×m。

(3)對協(xié)方差矩陣Dm×m進(jìn)行特征值分解，即D·V=Λ·V。其中，Λ是特征值矩陣，Λ=diag(λ1,λ2,…,λm)，且λ1≥λ2≥…≥λm≥0；V是m×m的特征向量矩陣，且VVT=I。

(4)計算X在V中的投影，即Y=VTX，得到矩陣X的主成分y1、y2、…、ym。

2 L曲線法

L曲線法最先由Lawson等[17]提出，然后被Hansen[18]用來確定正則化法的參數(shù)，以便求解病態(tài)問題。該方法認(rèn)為：把殘差項作為橫坐標(biāo)，正則化項作為縱坐標(biāo)，然后通過L曲線的角點即可確定所要求的正則化參數(shù)。所獲得的曲線具有L-形狀，所以命名為L曲線法。實際上早在1989年Hansen[19]就已經(jīng)證明了Tikhonov正則化方法和奇異值分解(包括廣義奇異值分解)的關(guān)系，并認(rèn)為正則化參數(shù)控制了正則化項和殘差項的加權(quán)，而這些權(quán)值相當(dāng)于奇異值分解中的特征值，描述了其對各個特征分量的貢獻(xiàn)。在應(yīng)用PCA對信號進(jìn)行分析時，特征值表征了其相應(yīng)的特征向量(或稱為成分)對信號的貢獻(xiàn)程度，因此可以引入L曲線法作為信噪分離準(zhǔn)則。在基于L曲線的信噪分離準(zhǔn)則中，對特征值進(jìn)行降序排列，以特征值索引作為橫坐標(biāo)，特征值作為縱坐標(biāo)，那么所獲得的特征值曲線可以近似為L曲線。對于完全滿足L形狀的曲線，要確定其角點是比較容易的，但是對于不滿足L形狀的曲線，要確定其角點就比較困難。對于后一種情況，作者采用的方法是：首先求解相鄰兩個特征值的斜率，然后求前一斜率和后一斜率的比值，具有最大比值所對應(yīng)的特征點就是所要求的L曲線的角點。也就是說，從特征值的斜率曲線來看，具有最大突變的特征值就是信噪分離的基準(zhǔn)點。當(dāng)然，再結(jié)合衰減曲線趨勢對比法，可以獲得更為合理的結(jié)果。

3 實驗結(jié)果及分析

為了驗證主成分分析方法在航空瞬變電磁去噪中的正確性，首先對正演模擬響應(yīng)曲線進(jìn)行主成分分析，確定有效信號的分離準(zhǔn)則；然后使用實測數(shù)據(jù)進(jìn)一步檢驗該方法的有效性和適用性。

3.1 模擬數(shù)據(jù)實驗

3.1.1 模擬數(shù)據(jù)的產(chǎn)生及主成分分析結(jié)果

MAXWELL軟件是一款由澳大利亞EMIT公司開發(fā)、基于Windows操作系統(tǒng)平臺的電磁正反演和處理軟件。本文的模擬數(shù)據(jù)由MAXWELL軟件中的有限元法正演得到，地電模型為層狀介質(zhì)。首先正演產(chǎn)生了15 s的模擬數(shù)據(jù)，同時也產(chǎn)生模擬的隨機噪聲和天電噪聲，然后通過五個周期疊加和抽道，分別獲得了75個測點的模擬數(shù)據(jù)和含噪數(shù)據(jù)，且每個測點的時間道數(shù)為20個。圖1(a)給出了單測點的模擬衰減曲線，圖1(b)所示為含噪的模擬衰減曲線。

根據(jù)上述的計算步驟，對圖1(b)所示的含噪電磁數(shù)據(jù)進(jìn)行主成分分析，獲得20個成分如圖2所示。

圖1 模擬電磁響應(yīng)衰減曲線Fig．1 The attenuation curves of EM response simulated(a)模擬響應(yīng)曲線； (b) 含噪的模擬響應(yīng)曲線

圖2 20個互不相關(guān)的成分Fig．2 Plots of 20 uncorrelated components

圖3 成分趨勢對比Fig．3 Contrastive plots of component trend

從圖2中可以看出，響應(yīng)曲線經(jīng)主成分分析后，把原始數(shù)據(jù)分解為20個互不相關(guān)的成分，而且每一個成分都是按照特征值的降序排列的。第一個成分比其他的成分都光滑，更接近原始的衰減曲線；隨著特征值的減小，所對應(yīng)的成分出現(xiàn)了不同程度的波動，說明在這些成分中存在很強的噪聲。

盡管從圖2中可以看到隨著特征值的減小，數(shù)據(jù)各個成分的變化趨勢，但不能明確地判別出反映地下信息的有效成分。為了能合理獲得反映地下介質(zhì)分布的有效成分，本文使用曲線趨勢對比法和L曲線法對分解后的成分進(jìn)行分析，從而分離出分別代表有效信號和噪聲的成分。

3.1.2 基于曲線趨勢對比的信噪分離準(zhǔn)則

首先用表示衰減曲線趨勢的第一主成分與其他的各個成分合并，然后與衰減曲線的主要趨勢(第一主成分)作對比[15]，如圖3所示。

第一主成分代表衰減曲線的主要趨勢，因此其反映了地下介質(zhì)分布的主要信息?，F(xiàn)在主要判斷其他成分到底反映的是地下介質(zhì)信息還是噪聲？為此，分別合并第一主成分和其他成分，如果某成分代表有效信號，那么該成分和第一主成分合并后所獲得的曲線趨勢不同于第一主成分；如果某成分代表噪聲，那么合并后的曲線趨同于第一主成分，但是其中會有個別點發(fā)生跳動，這是因為噪聲不可能改變整個曲線趨勢，而只是影響曲線的局部形狀。

利用上述的信噪分離準(zhǔn)則分析圖3，可以看出，成分1分別與成分2、成分3合并后的整體曲線變化趨勢不同于成分1，而成分4-成分20分別與成分1合并后的曲線除了在個別點有不同程度的跳變外，整體的趨勢基本與成分1一致，所以認(rèn)為成分1、成分2、成分3 是有效地質(zhì)成分，成分4-成分20 是噪聲成分。

3.1.3 基于L曲線的信噪分離準(zhǔn)則

盡管基于曲線趨勢對比法可以區(qū)分地質(zhì)成分和噪聲成分，但是它主要是依靠人來判斷，因此具有一定的主觀性。為了更為客觀地、定量地判別地質(zhì)成分和噪聲成分，這里給出了基于L曲線的信噪分離準(zhǔn)則。同時，該準(zhǔn)則也可以和曲線趨勢對比法進(jìn)行相互驗證。

圖4(a)為含噪的模擬電磁響應(yīng)數(shù)據(jù)的L曲線(歸一化)，該曲線表征了特征值的變化情況。一般來講，含噪數(shù)據(jù)的L曲線可明顯分為兩個部分：對應(yīng)有效地質(zhì)成分的特征值幅度較大，衰減較快；而對應(yīng)噪聲成分的特征值幅度要小得多，衰減較慢且數(shù)值變化平穩(wěn)。根據(jù)這個特性，很容易確定L曲線的角點位置，如圖4(a)中第3個特征值。如果要定量地確定L曲線的角點位置，那么可以先求出相鄰兩個特征值的斜率，然后依次求出前一斜率和后一斜率的比值，最后確定出最大的斜率比值所對應(yīng)的特征值，其中斜率比值如圖4(b)所示。

由圖4(a)可見，第3個特征值之前的曲線幅值較大，變化較快；而其值之后的曲線幅值較小，變化緩慢，所以第3個特征值就是L曲線的角點位置。此外，在圖4(b)中，我們也可以發(fā)現(xiàn)，具有最大斜率比值所對應(yīng)的特征值正是第3個特征值。實際上，斜率比值反映了特征值曲線變化率的突變情況。所以應(yīng)用L曲線法可以定量判別地質(zhì)成分和噪聲成分。因此根據(jù)上述方法，我們判定成分1、成分2、成分3代表有效的地質(zhì)成分，而成分4-成分20主要反映噪聲信息。

從上述分析可以看到，分別應(yīng)用基于L曲線的信噪分離準(zhǔn)則和基于曲線趨勢對比的分離準(zhǔn)則所獲得的結(jié)果是一致的。重構(gòu)成分1、成分2、成分3獲得的單點衰減曲線如圖5所示。從圖5中可以看到，噪聲得到了較好地抑制。

圖4 含噪數(shù)據(jù)的L曲線Fig．4 L curve (a)含噪模擬數(shù)據(jù)的L曲線；(b)特征值斜率比值

圖5 去噪后的電磁響應(yīng)曲線Fig．5 The attenuation curve of EM response denoised

3.2 實測數(shù)據(jù)實驗

圖6為吊艙式直升機時間域電磁測量系統(tǒng)在某一勘探區(qū)的實測單點衰減曲線，其中該數(shù)據(jù)已經(jīng)進(jìn)行了疊加、抽道、歸一化處理。疊加和抽道主要是為了抑制隨機噪聲。從圖6中，還可以發(fā)現(xiàn)，盡管通過疊加抽道可以抑制大部分隨機噪聲，但是不能有效減弱天電噪聲和人文噪聲的影響。對圖6中的數(shù)據(jù)應(yīng)用主成分分析法處理后，獲得的去噪結(jié)果如圖7所示。對比分析圖6和圖7可以發(fā)現(xiàn)，經(jīng)過主成分分析去噪后衰減曲線變得較為光滑，這說明主成分分析法能有效抑制航空瞬變電磁數(shù)據(jù)中的天電噪聲和人文噪聲，且計算較為簡單。

圖6 實測電磁響應(yīng)衰減曲線Fig．6 A raw decay measured

圖7 去噪后的電磁響應(yīng)衰減曲線Fig．7 A denoised decay

4 結(jié)論

主成分分析是一種統(tǒng)計分析方法，從統(tǒng)計角度能對有效信號和噪聲進(jìn)行分離。從上述的實驗結(jié)果及分析可知，主成分分析能有效抑制航空瞬變電磁法中的天電噪聲和人文噪聲。此外，主成分分析主要包含分解和重構(gòu)兩個部分，在分解時是對整個勘探區(qū)的數(shù)據(jù)進(jìn)行分解，所獲得特征向量比較穩(wěn)定；而在重構(gòu)時所要計算的成分大為減少，所以基于主成分的航空瞬變電磁去噪方法不僅算法穩(wěn)定，而且計算效率較高。

參考文獻(xiàn)：

[1] 雷棟, 胡祥云, 張素芳. 航空電磁法發(fā)展現(xiàn)狀[J]. 地質(zhì)找礦論叢, 2006, 21(1): 42-44.

[2] MCCRACKEN K G, PIK J P, HARRIS R W. Noise in EM exploration systems[J]. Exploration geophysics, 1984, 15(3): 169-174.

[3] SZARKA L. Geophysical aspects of man-made electromagnetic noise in the earth—review[J]. Surveys in geophysics, 1988, 9(3-4): 287-318.

[4] 陶本藻. 論最小二乘擬合[J]. 武漢大學(xué)學(xué)報：信息科學(xué)版, 1980(1): 27-34.

[5] AL-FOUZAN F, WILLIAM HARBERT, ROBERT DILMORE，et al. Methods for Removing Signal Noise from Helicopter Electromagnetic Survey Data[J]. Mine Water and the Environment, 2004, 23(1): 28-33.

[6] MACNAE J C, LAMONTAGNE Y, WEST G F. Noise processing techniques for time- domain EM systems[J]. Geophysics, 1984, 49(7): 934-948.

[7] SCOTT R, ATKINS F, HARPER P V. Median window filter as a smoothing-edge preserving technique[J]. Journal of nuclear medicine, 1978, 19(6):749.

[8] BUSELLI G, CAMERON M. Robust statistical methods for reducing sferics noise contaminating transient electromagnetic measurements[J]. Geophysics, 1996, 61(6):1633-1646.

[9] 呂東偉, 毛立峰. 時間域航空電磁資料小波去噪方法研究[C]. 第九屆中國國際地球物理電磁學(xué)術(shù)討論會論文, 2009.

[10] SPIES B R. Local noise prediction filtering for central induction transient electromagnetic sounding[J]. Geophysics, 1988, 53(8): 1068-1079.

[11] OLSEN K B, HOHMANN G W. Adaptive noise cancellation for time-domain EM data[J]. Geophysics, 1992, 57(3):466-469.

[12] DOUGLAS C, NYMAN, JAMES E G, et al. Adaptive rejection of high-line contamination[C]. 1983 SEG annual meeting, 1983.

[13] 夏江海.奇異值分解在位場資料中的應(yīng)用[J]. 物探化探計算技術(shù), 1989, 11(2): 93-98.

[14] 王權(quán)海,苗放.提高地震數(shù)據(jù)信噪比的PCA方法[J]. 地球物理學(xué)進(jìn)展, 2011, 26 (3): 1039-1044.

[15] RENINGER P A, MARTELET G. Singular value decomposition as a denoising tool for airborne time domain electromagnetic data[J]. Journal of Applied Geophysics, 2011, 75(2) : 264-276.

[16] 黃皓平.電磁法數(shù)據(jù)處理的奇異值分解法[J]. 地球物理學(xué)報, 1991, 34(5): 644-650.

[17] LAWSON C L, HANSON R J. Solving least squares problems[M]. Englewood Cliffs: Prentice-Hall, 1974.

[18] HANSEN P C. Analysis of discrete ill-posed problems by means of the L-curve[J]. SIAM review, 1992, 34(4): 561-580.

[19] HANSEN P C. Regularization, GSVD and truncaed GSVD[J]. BIT numerical mathematics, 1989, 29(3): 491-504.