郭 明, 王仕興, 易國財, 何 可, 張振雄, 王緒本

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

0 引言

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

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

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

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

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

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

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

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

(4)

式中：U1和U2分別為極點p1、p2到點ω的矢量長度，可以表示為：

(5)

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

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

(6)

因此，當(dāng)

(7)

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

Δω≈(1-r)

(8)

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

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

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

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

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

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

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

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

2 模擬信號消噪效果

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

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

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

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

(9)

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

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

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

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

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

圖6 模擬含噪信號和消噪信號疊加后對比Fig.6 Comparison between the simulated noise signal and the denoised signal after superposition

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

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

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

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

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

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

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

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

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

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

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

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

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

圖11(b)為數(shù)字諧振器去除工頻干擾后的反演結(jié)果，測線0 m～25 m段和385 m～400 m段為沱江的兩岸，反演電阻率為相對高阻，與實際地質(zhì)情況相符合。25 m～385 m段為江面，反演電阻率呈現(xiàn)低-高分布，推測藍(lán)色低阻為江水，水深約為20 m；在深度20 m～60 m處，電阻率由淺層到深層逐漸升高，推測為基巖，巖性為砂泥巖互層。但是由于無人機半航空瞬變電磁法對低阻地質(zhì)體反應(yīng)敏感，所以實際水深可能不足20 m。圖12為該測線的地質(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ù)字諧振器消除半航空實測中工頻干擾的方法，討論了不同帶寬諧振器的消噪特點。

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

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