一種基于數(shù)字諧振器重構(gòu)半航空瞬變電磁信號(hào)壓制工頻干擾的方法

郭 明, 王仕興, 易國(guó)財(cái), 何 可, 張振雄, 王緒本

(成都理工大學(xué) 地球勘探與信息技術(shù)教育部重點(diǎn)實(shí)驗(yàn)室，成都 610059)

0 引言

半航空瞬變電磁法(SATEM)，是一種融合了地面TDEM和航空VTEM的優(yōu)點(diǎn)的新勘探方法，該方法一般由地面長(zhǎng)接地導(dǎo)線源和懸掛在無(wú)人機(jī)上的接收機(jī)組成，其不僅具有較高的空間分辨率，還適合在河流、沼澤及其他地面工作難以展開(kāi)的地區(qū)工作。該方法最先由國(guó)外學(xué)者M(jìn)ogi[1-2]提出，將基于直升機(jī)的低空電磁探測(cè)系統(tǒng)成功應(yīng)用于Mount Bandai火山結(jié)構(gòu)探查；目前國(guó)內(nèi)研究熱點(diǎn)主要體現(xiàn)在理論研究[3-4]和儀器研究方面[5-6]，但實(shí)際勘探項(xiàng)目的相關(guān)應(yīng)用成果仍發(fā)表較少[7-8]。

半航空瞬變電磁法采用無(wú)人機(jī)搭載接收線圈進(jìn)行連續(xù)的數(shù)據(jù)采集，相對(duì)地面物探方法而言具有快速高效等優(yōu)點(diǎn)，但是其觀測(cè)二次場(chǎng)晚期信號(hào)較弱，容易受噪聲影響。當(dāng)前半航空瞬變電磁數(shù)據(jù)處理方法研究熱點(diǎn)，主要集中在小波多分辨率分析[9]和小波閾值消噪方面[10]，這些方法一般用于處理動(dòng)態(tài)噪聲、白噪聲以及天電噪聲等，而對(duì)工頻干擾消噪的文章較少。然而工頻干擾作為一種城市地區(qū)常見(jiàn)的電磁干擾，將會(huì)嚴(yán)重影響半航空瞬變電磁勘探的數(shù)據(jù)質(zhì)量，限制該方法在城市地下空間勘探中的應(yīng)用，瞬變電磁法中常采用雙極性同步采樣抑制工頻干擾，即發(fā)射機(jī)采用雙極性發(fā)射，在接收機(jī)中或者軟件中，把疊加之后的正極性信號(hào)與負(fù)極性信號(hào)相減，即可得兩次測(cè)量之和，以此達(dá)到消除外部工頻周期性噪聲的目的。理論上當(dāng)發(fā)射機(jī)周期為工頻周期的偶數(shù)倍時(shí)，在接收到的正負(fù)信號(hào)同時(shí)間道處，工頻信號(hào)相等，正、負(fù)信號(hào)相減后，即可消除工頻信號(hào)干擾[11]，但對(duì)于野外實(shí)測(cè)工頻信號(hào)，該方法僅能壓制部分干擾，難以達(dá)到半航空瞬變電磁法反演的數(shù)據(jù)質(zhì)量要求。所以筆者結(jié)合實(shí)際勘探結(jié)果，提出一種數(shù)字諧振器重構(gòu)信號(hào)去除工頻干擾的方法。

這里主要以無(wú)人機(jī)半航空瞬變電磁法在長(zhǎng)江某地的成功勘探成果為例，研究了該方法在城市地下空間勘探中工頻干擾對(duì)勘探數(shù)據(jù)的影響，以及數(shù)字諧振器重構(gòu)信號(hào)消除工頻干擾的效果。

1 半航空數(shù)字諧振器設(shè)計(jì)

Dupis等[12]對(duì)直流電力線工頻干擾及其諧波進(jìn)行了頻譜分析，結(jié)果如圖1所示。為研究工頻干擾的特點(diǎn)，筆者先利用采樣頻率為32 kfs的半航空瞬變電磁儀，采集該地區(qū)地面強(qiáng)工頻干擾的數(shù)據(jù)(圖2(a))，然后對(duì)實(shí)測(cè)數(shù)據(jù)進(jìn)行頻譜分析，分析結(jié)果如圖2(b)所示，該地區(qū)的強(qiáng)工頻干擾數(shù)據(jù)頻譜主要為50 Hz、250 Hz、350 Hz和550 Hz。

圖1 直流電力線工頻干擾頻譜Fig.1 Dc power line interference spectrum

數(shù)字諧振器是一種只通過(guò)信號(hào)單頻成分的系統(tǒng)，可以設(shè)計(jì)一種僅能通過(guò)如圖2(b)中頻率的數(shù)字諧振器，將原始數(shù)據(jù)與諧振輸出的數(shù)據(jù)相減作為消噪后的數(shù)據(jù)，達(dá)到消除原始信號(hào)中工頻干擾的目的。

數(shù)字諧振器是一個(gè)二階濾波器，也是種特殊的雙極點(diǎn)濾波器。該濾波器有一對(duì)共軛極點(diǎn)re±jω0,r接近于1，幅度特性在ω0附近最大，相當(dāng)于在該頻率發(fā)生了諧振，故稱為數(shù)字諧振器。數(shù)字諧振器非常適合頻帶非常窄、難以用常規(guī)IIR濾波器實(shí)現(xiàn)的帶通濾波器。數(shù)字諧振器的零點(diǎn)有兩種放置方法，一種是一個(gè)零點(diǎn)放置在原點(diǎn),另一種是兩個(gè)零點(diǎn)分別放置在z=1和z=-1處，筆者設(shè)計(jì)的濾波器零點(diǎn)在原點(diǎn),一對(duì)共軛極點(diǎn)為re±jω0的數(shù)字濾波器[13]。

其系統(tǒng)函數(shù)為：

(1)

可以看出，此系統(tǒng)的幅度特性在ω0附近取最大值,選取b0使|H(ejω0)=1|,則將z=ejω0帶入得式(2)。

(2)

從而：

b0=(1-r)|(1-re-j2ω0)|=

(1-r)|(1-rcos 2ω0+jrsin 2ω0)|=

(3)

該系統(tǒng)在任意頻點(diǎn)的幅度特性可以寫(xiě)成式(4)。

(4)

式中：U1和U2分別為極點(diǎn)p1、p2到點(diǎn)ω的矢量長(zhǎng)度，可以表示為：

(5)

當(dāng)U1(ω)、U2(ω)取最小值時(shí)，該系統(tǒng)具有最大幅度，我們求其最小值，得式(6)。

4r(1+r2)sinω0cosω+4r2sin 2ω]=0

(6)

因此，當(dāng)

(7)

時(shí)，幅度取最大值，此時(shí)的頻譜值為諧振器的精確諧振頻率。由此可以看出，如果兩個(gè)極點(diǎn)非常接近單位圓，則ω0≈ωr，可以證明其3dB帶寬為：

Δω≈(1-r)

(8)

通過(guò)以上諧振器的系統(tǒng)函數(shù)分析可知，設(shè)計(jì)一個(gè)數(shù)字諧振器的主要步驟如下：

1)根據(jù)實(shí)際數(shù)據(jù)要求的帶寬Δω得到諧振器的r值。

2)根據(jù)r值和諧振頻率ωr得到ω0。

圖2 地面強(qiáng)工頻干擾數(shù)據(jù)及頻譜Fig.2 Ground power frequency interference data and spectrum(a)工頻干擾數(shù)據(jù);(b)工頻干擾頻譜

3)根據(jù)式(1)設(shè)計(jì)諧振器。

筆者先設(shè)計(jì)了頻率為25 Hz，帶寬分別為1 Hz、2 Hz、3 Hz和4 Hz的數(shù)字諧振器，并繪制了對(duì)應(yīng)的頻率響應(yīng)曲線(圖3)，用來(lái)討論不同帶寬諧振器的特性。

圖3 不同帶寬的諧振器頻率響應(yīng)Fig.3 Frequency response of resonators with different bandwidths(a)幅頻特性;(b)相頻特性

理論上諧振器的帶寬越窄，消噪效果越明顯，但是從圖2(b)可以看出實(shí)測(cè)信號(hào)中工頻干擾不是標(biāo)準(zhǔn)的50 Hz，諧振器去噪時(shí)頻帶過(guò)窄導(dǎo)致去噪效果不佳，頻帶過(guò)寬會(huì)導(dǎo)致TDEM時(shí)間域衰減曲線畸變。

2 模擬信號(hào)消噪效果

為驗(yàn)證該方法的去噪效果，筆者正演模擬了發(fā)射頻率為25 Hz、采樣頻率為32 kfs，時(shí)間長(zhǎng)度為1 s的半航空瞬變電磁信號(hào)，將該信號(hào)作為有效信號(hào)s(n)，見(jiàn)圖4(a)。

圖4 半航空瞬變電磁響應(yīng)Fig.4 Semi-airbrone transient electromagnetic respons(a)理論信號(hào);(b)模擬含噪信號(hào)

由于工頻干擾不是標(biāo)準(zhǔn)的50 Hz，筆者用時(shí)間為1 s的50 Hz正弦波和平頂窗函數(shù)相乘，生成一個(gè)頻率為[50-p，50+p]的噪音信號(hào)模擬工頻干擾e(n)，其中噪音信號(hào)每個(gè)周期的頻移為p=5的定值；最后在有效信號(hào)s(n)上疊加模擬工頻干擾e(n)生成模擬含噪信號(hào)x(n)(圖4(b))。

為了找到寬帶合適的數(shù)字諧振器，達(dá)到消除噪聲，且對(duì)有效數(shù)據(jù)影響較小的效果，分別繪制了在雙對(duì)數(shù)坐標(biāo)系中，不同寬帶的消噪后的疊加數(shù)據(jù)曲線(圖5)，來(lái)對(duì)比不同寬帶的消噪效果，并引入相對(duì)擬合誤差作為衡量標(biāo)準(zhǔn)，本文擬合誤差(Rms)定義如下：

(9)

式中:si為理論數(shù)據(jù);xi為消噪后數(shù)據(jù);N為采樣時(shí)間個(gè)數(shù)。

圖5中消噪后曲線和原始數(shù)據(jù)曲線形態(tài)差異較小，根據(jù)相對(duì)擬合誤差(表1)可以得出，當(dāng)數(shù)字諧振器帶寬為2 Hz時(shí)，該消噪方法對(duì)有效數(shù)據(jù)影響較小。

圖5 不同帶寬數(shù)字諧振器消噪效果Fig.5 Noise elimination effect of digital resonator with different bandwidth

表1 不同帶寬數(shù)字諧振器消噪后信號(hào)相對(duì)擬合誤差

圖6中消噪選擇帶寬為2 Hz數(shù)字諧振器，原始數(shù)據(jù)和消噪后數(shù)據(jù)基本重合。從圖6可以看到，使用數(shù)字諧振器可以較好地壓制工頻干擾。

圖6 模擬含噪信號(hào)和消噪信號(hào)疊加后對(duì)比Fig.6 Comparison between the simulated noise signal and the denoised signal after superposition

3 野外實(shí)測(cè)數(shù)據(jù)消噪

此次長(zhǎng)江水域勘探目的是探明河流水深以及地下地質(zhì)結(jié)構(gòu)分布，該段江面寬度約為370 m，水流湍急，常規(guī)物探方法難以實(shí)施，所以選擇半航空瞬變電磁法作為勘探手段。通過(guò)對(duì)不同時(shí)間的場(chǎng)源進(jìn)行正演模擬(圖7)可知，場(chǎng)源在平行與線源方向變化速度較慢，垂直于場(chǎng)源的方向場(chǎng)源變化速度較快。

圖7 工區(qū)場(chǎng)源的正演模擬Fig.7 Forward simulation of field sources in the work area

根據(jù)對(duì)長(zhǎng)導(dǎo)線源的場(chǎng)源擴(kuò)散的特點(diǎn)，為保證采集數(shù)據(jù)質(zhì)量，此次勘探工作測(cè)線布置選擇垂直于江面布設(shè)長(zhǎng)度約500 m長(zhǎng)導(dǎo)線源，測(cè)線平行于線源，測(cè)線偏移距為200 m(圖8)。發(fā)射基頻為25 Hz、占空比50%的雙極方波電流，電流強(qiáng)度為10 A，接收機(jī)采樣頻率選擇32 kfs。

圖8 工區(qū)測(cè)線布置圖Fig.8 Work area survey line layou

此次勘探中數(shù)據(jù)整體質(zhì)量較好，但測(cè)線東南段接近岸邊的位置采集到的數(shù)據(jù)受到了50 Hz工頻信號(hào)及其諧波干擾的干擾，這些干擾可能由于公路上高壓輸電線引起。

利用本文的數(shù)字諧振器方法處理含噪數(shù)據(jù)，能有效地校正工頻畸變，然后將雙極性同步采樣所得信號(hào)進(jìn)行反向疊加，就可以有效地壓制工頻信號(hào)干擾。

通過(guò)圖9對(duì)比發(fā)現(xiàn)，數(shù)字諧振器消噪方法在損失較小有效信號(hào)的情況下，能在很大程度地壓制工頻干擾，為后期的反演解釋提供可靠的數(shù)據(jù)。需要注意的是，數(shù)字諧振器濾波過(guò)程中會(huì)有明顯的邊界效應(yīng)，處理數(shù)據(jù)時(shí)需要先對(duì)含噪的數(shù)據(jù)進(jìn)行延拓，然后再消噪處理。

圖9 實(shí)測(cè)數(shù)據(jù)和去噪后數(shù)據(jù)對(duì)比Fig.9 Comparison between measured data and denoised data

圖10為半航空瞬變電磁法在沱江測(cè)量的原始數(shù)據(jù)經(jīng)過(guò)數(shù)據(jù)處理方法得到的dbdt多測(cè)道圖，該剖面曲線由等時(shí)間間隔抽15道數(shù)據(jù)組成。圖10(a)為僅使用常規(guī)數(shù)據(jù)處理方法生成的dbdt多測(cè)道圖，在測(cè)線的330 m～400 m處由于受到岸邊工頻信號(hào)影響，剖面曲線的數(shù)據(jù)出現(xiàn)明顯的異常。異常主要體現(xiàn)為不同時(shí)間道的電磁感應(yīng)數(shù)據(jù)發(fā)生了交叉，且部分?jǐn)?shù)據(jù)出現(xiàn)了負(fù)值。圖10(b)為對(duì)野外數(shù)據(jù)先利用數(shù)字諧振器濾波去除工頻干擾后，再進(jìn)行常規(guī)處理，最后生成的dbdt多測(cè)道圖，該曲線的成功的去掉了工頻干擾的影響。

圖10 消噪前后dbdt多測(cè)道圖Fig.10 DBDT multi-track diagram before and after noise eliminationa(a)原始數(shù)據(jù)dbdt多測(cè)道圖;(b)消噪后dbdt多測(cè)道圖

為對(duì)比野外數(shù)據(jù)工頻干擾去除前后對(duì)反演結(jié)果的影響，利用半航空瞬變電磁自適應(yīng)正則化-阻尼最小二乘法反演方法進(jìn)行反演解釋[14],繪制了該測(cè)線段電阻率斷面圖。圖11(a)為數(shù)字諧振器濾波前的電阻率斷面圖，在圖中330 m～400 m處出現(xiàn)了明顯的異常，對(duì)比該處各測(cè)點(diǎn)的實(shí)測(cè)二次場(chǎng)衰減曲線和反演擬合曲線發(fā)現(xiàn)，這些測(cè)點(diǎn)的曲線擬合差很大，這樣的反演結(jié)果不可信，且對(duì)于出現(xiàn)負(fù)值的測(cè)點(diǎn)基本無(wú)法反演。

圖11 消噪前后反演結(jié)果Fig.11 Inversion results before and after noise reduction(a)原始數(shù)據(jù)反演結(jié)果;(b)諧振器消噪后數(shù)據(jù)的反演結(jié)果

圖11(b)為數(shù)字諧振器去除工頻干擾后的反演結(jié)果，測(cè)線0 m～25 m段和385 m～400 m段為沱江的兩岸，反演電阻率為相對(duì)高阻，與實(shí)際地質(zhì)情況相符合。25 m～385 m段為江面，反演電阻率呈現(xiàn)低-高分布，推測(cè)藍(lán)色低阻為江水，水深約為20 m；在深度20 m～60 m處，電阻率由淺層到深層逐漸升高，推測(cè)為基巖，巖性為砂泥巖互層。但是由于無(wú)人機(jī)半航空瞬變電磁法對(duì)低阻地質(zhì)體反應(yīng)敏感，所以實(shí)際水深可能不足20 m。圖12為該測(cè)線的地質(zhì)解釋結(jié)果，工區(qū)地質(zhì)概況屬于第四系更新統(tǒng)藍(lán)家坡組，上部為黃褐色砂質(zhì)亞粘土，下部為黃褐色礫石層，基巖為白堊系下統(tǒng)七曲寺組，以砂泥巖互層為主。

圖12 地質(zhì)解釋圖Fig.12 Geological interpretation graphics

4 結(jié)論

1)根據(jù)半航空瞬變電磁數(shù)據(jù)和工頻干擾的頻譜特性，筆者提出了一種基于數(shù)字諧振器消除半航空實(shí)測(cè)中工頻干擾的方法，討論了不同帶寬諧振器的消噪特點(diǎn)。

2)通過(guò)對(duì)比合成數(shù)據(jù)不同帶寬的數(shù)字諧振器效果，帶寬為2 Hz的效果較好，即不嚴(yán)重影響二次場(chǎng)衰減曲線，又能有效壓制50Hz工頻干擾。

3)將本文的數(shù)字諧振器技術(shù)應(yīng)用于沱江半航空瞬變電磁實(shí)測(cè)數(shù)據(jù)，有效壓制了工頻干擾。去噪后的二次場(chǎng)衰減曲線用于反演后，得到的電阻率剖面更符合地質(zhì)情況。