基于最小二乘反演的半航空瞬變電磁純二次場數(shù)據(jù)運(yùn)動(dòng)噪聲去除方法

楊洋， 張衡， 周長宇， 陳成棟， 孫懷鳳*

1 山東大學(xué)巖土與結(jié)構(gòu)工程研究中心， 濟(jì)南 250061 2 山東大學(xué)地球電磁探測研究所， 濟(jì)南 250061 3 山東省工業(yè)技術(shù)研究院先進(jìn)勘探與透明城市協(xié)同創(chuàng)新中心， 濟(jì)南 250061

0 引言

半航空瞬變電磁法(SATEM)最早是由Nabighian(1988)提出的一種地球物理勘探方法，它采用地面發(fā)射、空中接收的工作模式，這種模式結(jié)合了地面瞬變電磁法(TEM)大功率發(fā)射和航空瞬變電磁法(AEM)空中快速勘探的優(yōu)勢，克服了在復(fù)雜地形條件下傳統(tǒng)地面勘探手段無法有效開展和航空電磁勘探不適用且成本高的問題(嵇艷鞠等，2013, 2014).該方法具有觀測純二次場、適應(yīng)性強(qiáng)、靈敏度高和分辨率強(qiáng)的特點(diǎn)(孫懷鳳等，2019；張向陽等，2019)，能夠在我國中西部地區(qū)，尤其是地形復(fù)雜地區(qū)提供一種新的勘察技術(shù)手段(圖1).

圖1 半航空瞬變電磁勘探Fig.1 Semi-airborne transient electromagnetic exploration

由于接收線圈與飛行器在空中通過軟連接固定，導(dǎo)致飛行過程中接收線圈一直處于非穩(wěn)定運(yùn)動(dòng)狀態(tài)，從而產(chǎn)生運(yùn)動(dòng)噪聲，成為半航空瞬變電磁法影響最嚴(yán)重的噪聲之一.線圈運(yùn)動(dòng)噪聲是由接收線圈切割地磁場，導(dǎo)致線圈內(nèi)部磁通量變化而形成的感應(yīng)電動(dòng)勢，該噪聲頻率較低，并且伴隨著測量工作一直存在，是航空電磁勘探的主要噪聲源(殷長春，2018；毛鑫鑫等，2021).

國內(nèi)外學(xué)者在線圈運(yùn)動(dòng)噪聲的分析和去除方面取得了很大進(jìn)展，Lane等(2000)和Buselli等(1998)分別從時(shí)間域和頻率域分析了運(yùn)動(dòng)噪聲特點(diǎn)并得出運(yùn)動(dòng)噪聲主要能量集中在低頻域，頻率主要集中于0～1000 Hz.Munkholm(1997)利用三分量磁場之和在主場方向上的投影與運(yùn)動(dòng)噪聲之間耦合為最小原理抑制運(yùn)動(dòng)噪聲.Davis等(2006)考慮了線圈運(yùn)動(dòng)引起的系統(tǒng)幾何參數(shù)變化的影響，并設(shè)計(jì)濾波器進(jìn)行硬件補(bǔ)償去除運(yùn)動(dòng)噪聲.此外，還有一部分學(xué)者采用數(shù)據(jù)分解和擬合去除運(yùn)動(dòng)噪聲，Lemire(2001)采用樣條插值和拉格朗日優(yōu)化算法，去除航空電磁數(shù)據(jù)中的基線漂移；Kass等(2010)在反演運(yùn)算中利用主成分壓制噪聲；朱凱光和李楠(2009)研究了航空電磁探測數(shù)據(jù)中的多種噪聲及壓制技術(shù)；尹大偉等(2013)提出多項(xiàng)式擬合法去除運(yùn)動(dòng)噪聲；李肅義等(2013)在低空電磁數(shù)據(jù)處理中采用小波變換去除運(yùn)動(dòng)噪聲；Ji等(2016)采用小波閥值法去除時(shí)間域電磁數(shù)據(jù)中的運(yùn)動(dòng)噪聲，提高深部異常分辨率；劉富波等(2017)采用總體經(jīng)驗(yàn)?zāi)B(tài)分解法實(shí)現(xiàn)對運(yùn)動(dòng)噪聲的分離與剔除.隨著深度學(xué)習(xí)的興起，Xue等(2020)采用字典學(xué)習(xí)算法去除時(shí)間域航空電磁數(shù)據(jù)中的運(yùn)動(dòng)噪聲.由上述研究可知，航空電磁數(shù)據(jù)線圈運(yùn)動(dòng)噪聲去除研究較為充分.與航空電磁法類似，SATEM在野外數(shù)據(jù)采集過程中，會(huì)產(chǎn)生較大的運(yùn)動(dòng)噪聲，然而大部分勘探儀器只在關(guān)斷之后進(jìn)行數(shù)據(jù)采集，獲得的數(shù)據(jù)是非全時(shí)的，從而導(dǎo)致許多去噪方法無法有效去除線圈運(yùn)動(dòng)噪聲.

針對非全時(shí)SATEM數(shù)據(jù)，本文提出一種基于最小二乘的非全時(shí)SATEM數(shù)據(jù)運(yùn)動(dòng)噪聲去除方法.通過分析電磁數(shù)據(jù)衰減曲線的特點(diǎn)，將非全時(shí)電磁數(shù)據(jù)延拓至全時(shí)長，基于傅里葉級數(shù)構(gòu)造衰減晚期數(shù)據(jù)的超定線性方程組，通過最小二乘反演求解運(yùn)動(dòng)噪聲，并將所獲得運(yùn)動(dòng)噪聲從原數(shù)據(jù)中剔除.本文提出的運(yùn)動(dòng)噪聲去除方法得到了實(shí)際驗(yàn)證，對非全時(shí)SATEM數(shù)據(jù)預(yù)處理體現(xiàn)了明顯處理效果.

1 瞬變電磁純二次場數(shù)據(jù)運(yùn)動(dòng)特征分析

半航空瞬變電磁系統(tǒng)工作時(shí)，地面的發(fā)射源以一定電流，向地下發(fā)射一次脈沖場.當(dāng)激勵(lì)源的電流在某一關(guān)斷時(shí)間不為零的瞬間突然關(guān)斷后，地下探測區(qū)域的介質(zhì)將被激勵(lì)起呈環(huán)狀穿透的包含地電信息的二次場，以維持在斷開電流前的穩(wěn)定電磁場(一次場)，二次場的變化特性與地下探測介質(zhì)的電性分布密切相關(guān)，在一次場的關(guān)斷間隙觀測二次場的衰減特征.

圖2a給出實(shí)測的雙極性矩形發(fā)射波形，圖2b給出對應(yīng)的理想二次場衰減曲線，圖2c給出對應(yīng)的實(shí)測二次場衰減曲線.在一次場的關(guān)斷間隙，二次場隨時(shí)間呈指數(shù)衰減，可以根據(jù)其特性將其衰減曲線劃分為三個(gè)階段：早期、中期、晚期，各階段信號的特點(diǎn)如下：早期，信號強(qiáng)、能量大，衰減很快，有效電磁響應(yīng)信號和噪聲重疊；中期，電磁響應(yīng)信號逐漸減小，強(qiáng)度變?nèi)?，噪聲影響逐漸變大；晚期，有效電磁響應(yīng)信號很小，能量微弱，幾乎被噪聲淹沒.

在實(shí)測過程中，當(dāng)一次脈沖場突然關(guān)斷后，二次場的產(chǎn)生需要一定的過程.如圖2b所示，實(shí)測每個(gè)半周期數(shù)據(jù)二次場尖峰前存在能量增強(qiáng)段，而最前端低能量平滑區(qū)為噪聲部分，與衰減曲線晚期信號特點(diǎn)相似，以下統(tǒng)稱為晚期數(shù)據(jù).

圖3a給出雙極性矩形波的全時(shí)理論響應(yīng)，圖3b給出雙極性矩形波的非全時(shí)理論響應(yīng)，圖3c給出廣西某地實(shí)測非全時(shí)半航空瞬變電磁響應(yīng)剖面.從圖3c能夠很容易地看出每個(gè)半周期衰減曲線的起點(diǎn)和終點(diǎn)之間存在明顯的階躍，這是由于只在關(guān)斷后接收數(shù)據(jù)、供電時(shí)段不采集信號，造成所觀測到的數(shù)據(jù)是非全時(shí)的.同時(shí)由于線圈通過軟連接懸掛于無人機(jī)下方，無人機(jī)飛行時(shí)線圈姿態(tài)變化所引起運(yùn)動(dòng)噪聲較為明顯.在進(jìn)行地球物理反演之前，需要從觀測到的電磁數(shù)據(jù)中去除線圈運(yùn)動(dòng)噪聲.

圖2 半航空瞬變電磁二次場衰減曲線(a) 發(fā)射波形; (b) 理想衰減曲線; (c) 實(shí)測衰減曲線.Fig.2 Decay curve of semi-airborne transient electromagnetic secondary field(a) Emission waveform;(b) Ideal decay curve; (c) Measured decay curve.

2 基于最小二乘反演的運(yùn)動(dòng)噪聲去除方法

2.1 純二次場數(shù)據(jù)的全時(shí)延拓

如前所述，很多瞬變電磁儀器供電時(shí)并不采集信號，只采集關(guān)斷后信號，因此造成瞬變電磁時(shí)域數(shù)據(jù)并非全時(shí)的.線圈姿態(tài)隨時(shí)間連續(xù)變化，因此，每兩個(gè)相鄰半周期衰減曲線起點(diǎn)和終點(diǎn)之間存在明顯的階躍，導(dǎo)致運(yùn)動(dòng)噪聲以及有效電磁數(shù)據(jù)的連續(xù)性消失.本文首先將數(shù)據(jù)延拓至全時(shí)長，使各點(diǎn)對應(yīng)到二次場響應(yīng)的真實(shí)時(shí)刻.

將響應(yīng)剖面各樣本點(diǎn)根據(jù)其所屬的半周期進(jìn)行分割，假設(shè)共a個(gè)樣本點(diǎn)，b個(gè)半周期，則每個(gè)半周期有a/b個(gè)樣本點(diǎn).各個(gè)樣本點(diǎn)之間時(shí)間間隔一致，每兩個(gè)相鄰半周期之間存在一個(gè)半周期的時(shí)間間隔.假設(shè)一個(gè)半周期時(shí)間間隔為單位“1”，將實(shí)測第1、2、3,…,b個(gè)半周期樣本數(shù)據(jù)分別移動(dòng)到第1、3、5,…,2b-1半周期位置處，第2、4、6,…,2b位置不做任何處理.

圖4a給出圖3c實(shí)測響應(yīng)剖面延拓至全時(shí)長的結(jié)果圖.已知廣西某地半航空實(shí)測數(shù)據(jù)共9600個(gè)樣本點(diǎn)，32個(gè)半周期，每個(gè)半周期含有300個(gè)樣本點(diǎn).因此，將第1～300、第301～600,…、第9301～9600樣本點(diǎn)分別移動(dòng)到第1～300、第601～900,…、第18901～19200樣本點(diǎn)位置.

從圖3c和圖4a的比較中可以明顯看出，不同半周期衰減曲線之間的連續(xù)性已經(jīng)恢復(fù).而且，在線圈運(yùn)動(dòng)噪聲去除方面，連續(xù)數(shù)據(jù)明顯優(yōu)于不連續(xù)數(shù)據(jù).上述全時(shí)長延拓方式，符合噪聲的連續(xù)性規(guī)律，能夠獲得更好的線圈運(yùn)動(dòng)噪聲擬合結(jié)果.需要注意的是，如果采集的半航空瞬變電磁數(shù)據(jù)為全時(shí)長數(shù)據(jù)，則無需進(jìn)行上述延拓，同時(shí)，本文所提出去噪方法同樣適用于全時(shí)數(shù)據(jù).

圖3 半航空瞬變電磁雙極性矩形波響應(yīng)剖面(a) 全時(shí)理論響應(yīng); (b) 非全時(shí)理論響應(yīng); (c) 實(shí)測響應(yīng)剖面.Fig.3 Semi-airborne transient electromagnetic bipolar rectangular wave response profile(a) Full-time theoreticalresponse; (b) Non-full-time theoretical response; (c) Measured response profile.

圖4 數(shù)據(jù)預(yù)處理(a) 非全時(shí)瞬變電磁數(shù)據(jù)延拓至全時(shí)長; (b) 非全時(shí)瞬變電磁衰減晚期數(shù)據(jù).Fig.4 Data preprocessing(a) Extension of non-full-time TEM data to full-time; (b) Late data of non-full-time transient electromagnetic decay curve.

特別強(qiáng)調(diào)的是，每個(gè)半周期二次場早、中期的衰減特性會(huì)對擬合運(yùn)動(dòng)噪聲的發(fā)展趨勢產(chǎn)生影響.如圖4b所示，基于瞬變電磁晚期數(shù)據(jù)有效信號很小、幾乎被噪聲淹沒的特點(diǎn)，從圖4a中剔除每個(gè)半周期衰減早、中期的數(shù)據(jù)，選擇衰減晚期的數(shù)據(jù)點(diǎn)作為噪聲在時(shí)間域的已知點(diǎn)，重建全時(shí)長的運(yùn)動(dòng)噪聲.

2.2 基于最小二乘反演的運(yùn)動(dòng)噪聲去除方法

線圈運(yùn)動(dòng)噪聲是由線圈姿態(tài)變化引起磁通量改變而產(chǎn)生的感應(yīng)電動(dòng)勢，線圈的姿態(tài)變化類似于單擺運(yùn)動(dòng)，具有一定的周期性特征，而傅里葉級數(shù)對表示周期性信號或者近似周期信號具有天然的優(yōu)勢(Bracewell,1966; Yang et al.,2018).因此，本文利用傅里葉正交基來重建全時(shí)長的運(yùn)動(dòng)噪聲.更具體地說，通過傅里葉級數(shù)構(gòu)造超定線性方程組來擬合運(yùn)動(dòng)噪聲.

通過離散傅立葉逆變換構(gòu)建針對離散信號的傅立葉級數(shù)，離散傅里葉逆變換公式為：

(1)

(2)

假設(shè)等式中f(x)為運(yùn)動(dòng)噪聲的連續(xù)采樣，則F(y)為運(yùn)動(dòng)噪聲對應(yīng)的離散頻率域系數(shù)，但事實(shí)上只有部分運(yùn)動(dòng)噪聲是準(zhǔn)確知曉的，即瞬變電磁晚期數(shù)據(jù)，而在早期，有效信號與運(yùn)動(dòng)噪聲重疊是無法直接分離的，運(yùn)動(dòng)噪聲一般以中低頻為主，能量亦主要集中于中低頻段，因此刻畫運(yùn)動(dòng)噪聲并不需要全部頻率域系數(shù)，只需要部分頻率域系數(shù)即可，尤其是低頻段對應(yīng)頻率域系數(shù)，如果噪聲在時(shí)間域已知點(diǎn)的數(shù)目大于頻率域未知系數(shù)的數(shù)目，則可以通過構(gòu)建超定方程組求解噪聲的頻率域系數(shù).

不妨設(shè)運(yùn)動(dòng)噪聲在頻率域中有兩個(gè)期望頻率(實(shí)際上可能有更多頻率)，即四個(gè)頻率域系數(shù)，則可以通過選擇時(shí)間域中任意四個(gè)值準(zhǔn)確已知的點(diǎn)來構(gòu)造一個(gè)適定的線性方程組求得唯一解，即：

(3)

式中，右邊紅色變量為期望頻率，即運(yùn)動(dòng)噪聲對應(yīng)頻率.左邊藍(lán)色變量為在時(shí)間域中所選中數(shù)據(jù)點(diǎn).

由于實(shí)際數(shù)據(jù)中除運(yùn)動(dòng)噪聲外仍存在高斯噪聲等其他噪聲，所以需要在時(shí)間域中選擇更多的點(diǎn)構(gòu)造一個(gè)超定的線性方程組求解運(yùn)動(dòng)噪聲.在這樣的情況下，假設(shè)選擇了六個(gè)時(shí)間域數(shù)據(jù)點(diǎn)，有四個(gè)待求的頻率系數(shù)，超定的線性方程組如式(4)所示：

(4)

將式(4)改寫為一個(gè)矩陣向量公式，可得：

Ax=b，

(5)

其中：

(6)

將(5)式兩邊左乘A的轉(zhuǎn)置矩陣AT，得到：

ATAx=ATb，

(7)

在式(7)兩邊左乘(ATA)-1，利用最小二乘求解期望頻率系數(shù)向量x，得到：

x=(ATA)-1ATb.

(8)

然后，將求解得到的期望頻率系數(shù)向量x中各個(gè)元素放到在頻率域?qū)?yīng)的位置，利用離散傅里葉逆變換反演得到全時(shí)長的時(shí)間域運(yùn)動(dòng)噪聲，并將其從對應(yīng)位置的非全時(shí)電磁數(shù)據(jù)中去掉以達(dá)到去除運(yùn)動(dòng)噪聲的目的.公式(9)給出運(yùn)動(dòng)噪聲構(gòu)造方法，其中，我們利用4個(gè)不同頻率系數(shù)構(gòu)造整個(gè)時(shí)間序列噪聲，即：

(9)

更一般的，當(dāng)實(shí)測瞬變電磁數(shù)據(jù)均值明顯偏移二次場響應(yīng)能量為0位置時(shí)，可以在公式(6)中矩陣A最右側(cè)加入一列一階勒讓德多項(xiàng)式，修正曲線的偏移，以得到更準(zhǔn)確的結(jié)果.

一階勒讓德多項(xiàng)式可表示為：

p1(x)=x，

(10)

式中，x∈[-1,1]，x個(gè)數(shù)等于N，各點(diǎn)間隔2∕(N-1)，則求解矩陣可表示為：

(11)

式中，p1(m1),…,p1(m6)代表樣本點(diǎn)在時(shí)間域?qū)?yīng)位置的一階勒讓德多項(xiàng)式，P1(k)代表一階勒讓德多項(xiàng)式系數(shù).

由于數(shù)據(jù)本身是非全時(shí)的，無法直接通過傅里葉變換獲得運(yùn)動(dòng)噪聲期望頻譜.為得到合理的期望頻率，本文首先通過勒讓德多項(xiàng)式對延拓至全時(shí)長的非全時(shí)數(shù)據(jù)進(jìn)行擬合，對擬合曲線進(jìn)行頻譜分析，從而估計(jì)運(yùn)動(dòng)噪聲對應(yīng)期望頻率.勒讓德多項(xiàng)式階數(shù)的選擇，則以衰減晚期數(shù)據(jù)去除多項(xiàng)式擬合曲線后剩余噪聲平均能量小于去噪前相對于均值的平均能量的10%為準(zhǔn).由于運(yùn)動(dòng)噪聲能量主要集中在低頻域(1 Hz～1 kHz)，且隨著頻率的增加噪聲幅值迅速衰減(Buselli et al.,1998)，因此，本文將基于勒讓德多項(xiàng)式擬合基線得到的頻譜進(jìn)行能量累加，累加范圍為1 Hz～1 kHz.基于大量實(shí)際數(shù)據(jù)測試的統(tǒng)計(jì)結(jié)果，最終選擇以80%累加能量和為閾值確定期望頻率最大值，以此確定實(shí)測數(shù)據(jù)的運(yùn)動(dòng)噪聲期望頻率范圍，進(jìn)而基于期望頻率建立針對運(yùn)動(dòng)噪聲頻率系數(shù)的超定方程組，以求解完整運(yùn)動(dòng)噪聲.

3 合成模型算例與方法測試

針對仿真數(shù)據(jù)進(jìn)行去噪算法測試，假設(shè)仿真信號采樣頻率為30000 Hz，周期樣本點(diǎn)個(gè)數(shù)為600個(gè)，信號峰峰值為11.4343 mV，共32個(gè)周期,19200個(gè)采樣點(diǎn).

圖5a為全時(shí)仿真原始信號，圖5b為模擬運(yùn)動(dòng)噪聲信號，將圖5a、b合成得到圖5c所示的含運(yùn)動(dòng)噪聲全時(shí)仿真信號.進(jìn)而，將圖5c中第1、3、5,…,63半周期數(shù)據(jù)摘取出來，拼接成如圖5d所示的非全時(shí)瞬變電磁數(shù)據(jù)仿真信號，其中，圖5a所示全時(shí)仿真原始信號為雙極性矩形波正演得到的全時(shí)理論響應(yīng)，圖5b所示模擬運(yùn)動(dòng)噪聲信號由Matlab編程隨機(jī)產(chǎn)生的未知周期、未知幅值的三角函數(shù)構(gòu)成.

圖5 含運(yùn)動(dòng)噪聲非全時(shí)瞬變電磁信號仿真(a) 仿真原始信號; (b) 模擬運(yùn)動(dòng)噪聲; (c) 全時(shí)含噪仿真信號; (d) 非全時(shí)含噪仿真信號.Fig.5 Simulation of non-full-time transient electromagnetic signal with motion noise(a) Simulated original signal; (b) Simulated motion noise; (c) Full-time noisy simulated signal; (d) Non-full-time noisy simulated signal.

圖6a給出將圖5d非全時(shí)含噪仿真信號延拓至全時(shí)長的結(jié)果.圖7a給出仿真信號衰減曲線.其中，第1～10、51～300樣本點(diǎn)有效電磁響應(yīng)信號很幾乎被噪聲淹沒，為晚期數(shù)據(jù)，將其從圖6a中剔除，結(jié)果見圖6b.

圖6 勒讓德多項(xiàng)式擬合非全時(shí)含噪仿真信號運(yùn)動(dòng)噪聲(a) 含噪仿真信號延拓至全時(shí)長; (b) 含噪仿真信號晚期數(shù)據(jù); (c) 勒讓德多項(xiàng)式擬合基線; (d) 去噪后信號.Fig.6 Motion noise of non-full-time noisy simulation signal fitted by Legendre polynomials(a) Noisy simulation signal extended to full-time; (b) Late data of noisy simulation signal; (c) Baseline of Legendre polynomials fitting; (d) Signal after denoising.

圖7b給出選用階數(shù)1～10階的勒讓德多項(xiàng)式擬合后剩余噪聲平均能量值.去噪前相對于均值的平均能量為15.76 mV，勒讓德多項(xiàng)式最高階為6階時(shí)，剩余噪聲平均能量小于10%，因此采用6階勒讓德多項(xiàng)式初步擬合全時(shí)運(yùn)動(dòng)噪聲基線，從原始數(shù)據(jù)中直接剔除基線，由圖6d可見，去噪后相鄰半周期之間仍存在明顯的階躍，表明利用勒讓德多項(xiàng)式整體擬合方式的去噪效果并不理想.然而，如圖6c所示，擬合基線與原始數(shù)據(jù)發(fā)展趨勢幾乎一致，因此可以估計(jì)運(yùn)動(dòng)噪聲頻譜.實(shí)測數(shù)據(jù)也將采用上述方式進(jìn)行頻譜分析.短時(shí)傅里葉變換得到如圖8a所示頻譜，80%累加能量和(1 Hz～1 kHz)界限為82 Hz.

圖7 仿真信號衰減曲線和勒讓德多項(xiàng)式階數(shù)選取(a) 仿真信號衰減曲線; (b) 勒讓德多項(xiàng)式階數(shù)-剩余噪聲平均能量.Fig.7 Decay curve of simulated original signal and selection of the order of Legendre polynomials(a) Decay curve of simulated original signal; (b) The order of Legendre polynomials-average energy of residual noise.

圖8 勒讓德多項(xiàng)式擬合基線頻譜分析(a) 頻譜; (b) 1 Hz～1 kHz能量累加.Fig.8 Spectral analysis of baseline fitted by Legendre polynomials(a) Spectral; (b) Accumulation of energy from 1 Hz to 1 kHz.

采用0～82 Hz、頻率間隔為30000/(2×9600)Hz的傅里葉正交基擬合運(yùn)動(dòng)噪聲基線，系數(shù)矩陣A為8320×132的超定方程組.如圖9b所示，用非全時(shí)仿真信號減去擬合運(yùn)動(dòng)噪聲，得到圖9c所示的去噪數(shù)據(jù)，各半周期信號峰峰值能量高度一致，晚期數(shù)據(jù)平滑.模擬運(yùn)動(dòng)噪聲與擬合運(yùn)動(dòng)噪聲基線對比結(jié)果如圖9a所示，二者發(fā)展趨勢幾乎重疊.

圖9 傅里葉正交基及小波變換擬合非全時(shí)含噪仿真信號運(yùn)動(dòng)噪聲(a) 仿真運(yùn)動(dòng)噪聲與擬合運(yùn)動(dòng)噪聲基線; (b) 傅里葉正交基擬合基線; (c) 去除傅里葉正交基基線的信號; (d) 非全時(shí)仿真信號差分; (e) 小波變換擬合基線; (f) 去除小波變換基線的信號.Fig.9 Fitting motion noise of non-full-time noisy simulation signal by Fourier orthogonal basis and wavelet transform(a) Simulated motion noise and fitted motion noise baseline; (b) Baseline fitted by Fourier orthogonal basis; (c) The signal removed baseline fitted by Fourier orthogonal basis; (d) Non-full-time simulation signal differential; (e) Baseline fitted by wavelet transform; (f) The signal removed baseline fitted by wavelet transform.

為了對比不同方法的去噪效果，本文采用了小波分解的方式進(jìn)行去噪(Wang et al., 2013)，在應(yīng)用中采用了db5和分解層數(shù)9對仿真信號進(jìn)行處理，去噪結(jié)果如圖9f所示.運(yùn)動(dòng)噪聲在一部分樣本點(diǎn)被很好的抑制，但同時(shí)能夠看到，在某些位置去噪結(jié)果仍然不甚理想，特別是在每個(gè)周期兩端樣本點(diǎn)存在較明顯的跳躍.

利用小波分解對該信號處理效果不理想的原因，主要來自兩方面：(1)運(yùn)動(dòng)噪聲在時(shí)間上是相對連續(xù)的，但該信號為純二次場信號，其在時(shí)間上不連續(xù)的，因此不同周期對應(yīng)衰減信號之間是存在一定階躍的，如果直接將該信號進(jìn)行小波分解，將會(huì)人為引入新的噪聲，對處理效果造成影響； (2)小波分解是針對整個(gè)信號進(jìn)行的，由于本身瞬變電磁信號的頻譜相對較寬，低頻運(yùn)動(dòng)噪聲與瞬變電磁信號中均存在低頻成分，在頻率域兩者頻譜存在重疊，在利用小波分解去除低頻運(yùn)動(dòng)噪聲的同時(shí)，易對瞬變電磁信號中低頻成分產(chǎn)生影響，進(jìn)而造成去噪后的信號失真.相比小波分解方法，本文所提出方法，第一步：將純二次場信號進(jìn)行時(shí)間延拓，將信號恢復(fù)到其準(zhǔn)確時(shí)間軸上，以解決運(yùn)動(dòng)噪聲時(shí)間連續(xù)性問題；第二步：避開瞬變電磁信號幅值較大位置，截取不同周期信號中的晚期部分，由于其瞬變電磁有效信號幾乎衰減殆盡，其幾乎全為噪聲，基于該部分信號，通過最小二乘法擬合獲得全時(shí)間對應(yīng)的運(yùn)動(dòng)噪聲，避免瞬變電磁信號本身對基線去除的影響.

圖10給出采用本文方法含噪仿真信號所有半周期去噪前、后疊加衰減曲線與仿真原始信號單一半周期衰減曲線對比結(jié)果.由圖可以看出，加入運(yùn)動(dòng)噪聲后仿真信號衰減曲線明顯向上偏移，使用上述方法進(jìn)行去噪后，衰減曲線與原始仿真信號幾乎重疊，運(yùn)動(dòng)噪聲被有效去除，驗(yàn)證了本文所提出方法針對仿真信號的有效性.

圖10 仿真原始信號、仿真含噪信號與去噪后信號衰減曲線對比Fig.10 Comparison of the decay curve of simulated original signal, simulated noisy signal and denoised signal

4 實(shí)測數(shù)據(jù)去噪

為了進(jìn)一步驗(yàn)證本文方法對非全時(shí)半航空電磁數(shù)據(jù)運(yùn)動(dòng)噪聲去除的有效性，我們對廣西某地半航空瞬變電磁法實(shí)測非全時(shí)純噪聲數(shù)據(jù)以及非全時(shí)電磁數(shù)據(jù)進(jìn)行處理.

4.1 實(shí)測非全時(shí)純噪聲數(shù)據(jù)去噪實(shí)例

圖11a給出廣西某地半航空實(shí)測非全時(shí)純噪聲數(shù)據(jù)，即發(fā)射源未供電時(shí)，所采集純運(yùn)動(dòng)噪聲數(shù)據(jù)，其采樣頻率為30000 Hz,單個(gè)觀測點(diǎn)數(shù)據(jù)每個(gè)半周期樣本點(diǎn)個(gè)數(shù)為300，共32個(gè)半周期，9600個(gè)采樣點(diǎn).如圖11b所示，將觀測到的非全時(shí)數(shù)據(jù)延拓至全時(shí)長，將第301～600樣本點(diǎn)移動(dòng)至第601～900位置，第601～900樣本點(diǎn)移動(dòng)至第1201～1500位置，依次向后延拓.

圖11 勒讓德多項(xiàng)式擬合實(shí)測非全時(shí)純噪聲信號運(yùn)動(dòng)噪聲(a) 實(shí)測非全時(shí)純噪聲信號; (b) 非全時(shí)信號延拓至全時(shí)長; (c) 勒讓德多項(xiàng)式擬合基線.Fig.11 Motion noise of measured non-full-time pure noise signal fitted by Legendre polynomial(a) Measured non-full-time pure noisesignal; (b) Non-full-time signals extend to full-time; (c) Baseline fitted by Legendre polynomial.

圖12給出選用階數(shù)15～25階的勒讓德多項(xiàng)式擬合后剩余噪聲平均能量值.去噪前相對于均值的平均能量為15.132 mV，勒讓德多項(xiàng)式最高階為23階時(shí)，剩余噪聲平均能量小于10%，因此采用23階勒讓德多項(xiàng)式初步擬合全時(shí)運(yùn)動(dòng)噪聲基線，結(jié)果見圖11c.傅里葉變換得到如圖13a所示頻譜，80%累加能量和(1 Hz～1 kHz)界限為110 Hz，參見圖13b.

圖12 實(shí)測純噪聲數(shù)據(jù)勒讓德多項(xiàng)式階數(shù)選取Fig.12 Selection of the order of Legendre polynomials for measured pure noise data

圖13 勒讓德多項(xiàng)式擬合基線頻譜分析(a) 頻譜; (b) 1 Hz～1 kHz能量累加.Fig.13 Spectral analysis of baseline fitted by Legendre polynomials(a) Spectral; (b) Accumulation of energy from 1 Hz to 1 kHz.

選取每個(gè)半周期的第1～10、51～300位置的樣本點(diǎn)數(shù)據(jù)為衰減晚期數(shù)據(jù)，采用期望頻率為0～110 Hz、頻率間隔為30000/(9600×2)Hz的傅里葉正交基擬合運(yùn)動(dòng)噪聲基線，系數(shù)矩陣A為8320×160，即方程數(shù)目為8320個(gè)，待解未知數(shù)為160個(gè)，利用最小二乘法求解超定方程組可獲得160個(gè)待定系數(shù)，進(jìn)而用于構(gòu)建運(yùn)動(dòng)噪聲.

如圖14所示，運(yùn)動(dòng)噪聲擬合結(jié)果與原始數(shù)據(jù)的發(fā)展趨勢幾乎吻合，特別是衰減早、中期運(yùn)動(dòng)噪聲得到很好的擬合.消除運(yùn)動(dòng)噪聲之后，剩余噪聲在零值附近上下振蕩，無明顯波動(dòng)，表現(xiàn)了很好的去噪結(jié)果.

圖14 傅里葉正交基擬合非全時(shí)純噪聲數(shù)據(jù)的運(yùn)動(dòng)噪聲(a) 實(shí)測非全時(shí)純噪聲信號晚期數(shù)據(jù); (b) 傅里葉正交基擬合運(yùn)動(dòng)噪聲; (c) 去噪后信號.Fig.14 Motion noise of non-full-time pure noise data fitted by Fourier orthogonal basis(a) Late data of measured non-full-time pure noise signal; (b) Motion noise fitted by Fourier orthogonal basis; (c) Signal after denoising.

4.2 實(shí)測非全時(shí)瞬變電磁數(shù)據(jù)去噪實(shí)例

圖15a給出廣西某地半航空實(shí)測電磁數(shù)據(jù)，其采樣頻率為30000 Hz，單個(gè)觀測點(diǎn)數(shù)據(jù)每個(gè)半周期樣本點(diǎn)個(gè)數(shù)為300，共32個(gè)半周期，9600個(gè)采樣點(diǎn).參見圖15b，將觀測到的非全時(shí)數(shù)據(jù)延拓至全時(shí)長，將第301～600樣本點(diǎn)延拓至第601～900位置，第601～900樣本點(diǎn)延拓至第1201～1500位置，依次向后延拓.根據(jù)去噪前疊加衰減曲線，認(rèn)為每個(gè)半周期的第11～50位置為衰減早、中期數(shù)據(jù)，將其從圖15b中剔除，結(jié)果見圖15c.

圖15 勒讓德多項(xiàng)式擬合實(shí)測非全時(shí)電磁信號運(yùn)動(dòng)噪聲(a) 實(shí)測非全時(shí)瞬變電磁信號; (b) 非全時(shí)信號延拓至全時(shí)長; (c) 晚期數(shù)據(jù); (d) 勒讓德多項(xiàng)式擬合基線.Fig.15 Motion noise of measured non-full-time electromagnetic signal fitted by Legendre polynomial(a) Measured non-full-time transient electromagnetic signals; (b) Non-full-time signals extend to full-time; (c) Late data; (d) Baseline fitted by Legendre polynomial.

圖16給出選用階數(shù)5～15階的勒讓德多項(xiàng)式擬合后剩余噪聲平均能量值.去噪前相對于均值的平均能量為19.08 mV，勒讓德多項(xiàng)式最高階為10階時(shí)，剩余噪聲平均能量小于10%，因此采用10階勒讓德多項(xiàng)式初步擬合全時(shí)運(yùn)動(dòng)噪聲基線，結(jié)果如圖15d所示.傅里葉變換得到如圖17a所示頻譜，80%累加能量和(1 Hz～1 kHz)界限為132 Hz.

圖16 實(shí)測非全時(shí)數(shù)據(jù)勒讓德多項(xiàng)式最高階數(shù)選取Fig.16 Selection of the highest order of Legendre polynomials for measured non-full-time data

圖17 勒讓德多項(xiàng)式擬合基線頻譜分析(a) 頻譜; (b) 1 Hz～1 kHz能量累加.Fig.17 Spectral analysis of baseline fitted by Legendre polynomials(a) Spectral; (b) Accumulation of energy from 1 Hz to 1 kHz.

采用期望頻率為0～132 Hz、頻率間隔為30000/(2×9600)Hz的傅里葉正交基擬合運(yùn)動(dòng)噪聲，可以獲得系數(shù)矩陣A為8320×176的超定方程組.

如圖18所示，擬合基線與原始數(shù)據(jù)的發(fā)展趨勢幾乎吻合，運(yùn)動(dòng)噪聲得到很好的擬合.消除運(yùn)動(dòng)噪聲之后，二次場早期脈沖具有高度一致性，晚期數(shù)據(jù)平滑性也很好.傅里葉正交基擬合方法可以真實(shí)反映線圈運(yùn)動(dòng)噪聲.

圖18 傅里葉正交基擬合實(shí)測非全時(shí)電磁信號運(yùn)動(dòng)噪聲(a) 實(shí)測非全時(shí)信號與傅里葉正交基擬合基線; (b) 去噪后信號.Fig.18 Motion noise of measured non-full-time electromagnetic signal fitted by Fourier orthogonal basis(a) Measured non-full-time signal and baseline fitted by Fourier orthogonal basis; (b) Signal after denoising.

圖19給出非全時(shí)實(shí)測瞬變電磁信號所有半周期去噪前、后疊加衰減曲線對比結(jié)果，可以看出運(yùn)動(dòng)噪聲去除后疊加衰減曲線整體向下偏移，與圖10所示仿真信號一致，運(yùn)動(dòng)噪聲被有效去除，驗(yàn)證了第2節(jié)所述方法的有效性.

圖19 非全時(shí)實(shí)測瞬變電磁信號所有半周期去噪前、后疊加衰減曲線對比Fig.19 Comparison of superimposed decay curves before and after all half-period denoising of non-full-time measured transient electromagnetic signals

5 結(jié)論

(1)針對半航空瞬變電磁數(shù)據(jù)中的運(yùn)動(dòng)噪聲進(jìn)行去除時(shí)，由于運(yùn)動(dòng)噪聲是連續(xù)的，而采集信號只含關(guān)斷后數(shù)據(jù)，是非連續(xù)的，應(yīng)先將非全時(shí)數(shù)據(jù)投影為全時(shí)長數(shù)據(jù)，以此作為前提進(jìn)行運(yùn)動(dòng)噪聲的剔除.

(2)由于半航空瞬變電磁衰減曲線早、中期數(shù)據(jù)衰減特性會(huì)對運(yùn)動(dòng)噪聲的擬合產(chǎn)生較大影響，且易產(chǎn)生畸變，而晚期數(shù)據(jù)及脈沖前端部分?jǐn)?shù)據(jù)衰減信號極小，幾乎只含噪聲，因此通過擬合晚期數(shù)據(jù)及脈沖前端部分?jǐn)?shù)據(jù)，通過建立超定方程組求解運(yùn)動(dòng)噪聲傅里葉級數(shù)對應(yīng)系數(shù)，能夠獲得重建后的全時(shí)長運(yùn)動(dòng)噪聲.

(3)通過對仿真數(shù)據(jù)、廣西某地實(shí)測數(shù)據(jù)進(jìn)行運(yùn)動(dòng)噪聲剔除，去噪結(jié)果表明本文提出的基于最小二乘的半航空瞬變電磁運(yùn)動(dòng)噪聲去除方法，能夠獲得較好去噪效果，是一種針對半航空瞬變電磁運(yùn)動(dòng)噪聲的有效去噪方法.