主成分分析在航空瞬變電磁去噪中的應(yīng)用

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

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

0 前言

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

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

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

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

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

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

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

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

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

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

2 L曲線法

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

3 實(shí)驗(yàn)結(jié)果及分析

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

3.1 模擬數(shù)據(jù)實(shí)驗(yàn)

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

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

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

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

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

圖3 成分趨勢(shì)對(duì)比Fig．3 Contrastive plots of component trend

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

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

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

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

第一主成分代表衰減曲線的主要趨勢(shì)，因此其反映了地下介質(zhì)分布的主要信息。現(xiàn)在主要判斷其他成分到底反映的是地下介質(zhì)信息還是噪聲？為此，分別合并第一主成分和其他成分，如果某成分代表有效信號(hào)，那么該成分和第一主成分合并后所獲得的曲線趨勢(shì)不同于第一主成分；如果某成分代表噪聲，那么合并后的曲線趨同于第一主成分，但是其中會(huì)有個(gè)別點(diǎn)發(fā)生跳動(dòng)，這是因?yàn)樵肼暡豢赡芨淖冋麄€(gè)曲線趨勢(shì)，而只是影響曲線的局部形狀。

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

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

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

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

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

從上述分析可以看到，分別應(yīng)用基于L曲線的信噪分離準(zhǔn)則和基于曲線趨勢(shì)對(duì)比的分離準(zhǔn)則所獲得的結(jié)果是一致的。重構(gòu)成分1、成分2、成分3獲得的單點(diǎn)衰減曲線如圖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í)測(cè)數(shù)據(jù)實(shí)驗(yàn)

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

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

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

4 結(jié)論

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

參考文獻(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é)報(bào)：信息科學(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] 呂東偉, 毛立峰. 時(shí)間域航空電磁資料小波去噪方法研究[C]. 第九屆中國(guó)國(guó)際地球物理電磁學(xué)術(shù)討論會(huì)論文, 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] 夏江海.奇異值分解在位場(chǎng)資料中的應(yīng)用[J]. 物探化探計(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é)報(bào), 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.