一種基于數字諧振器重構半航空瞬變電磁信號壓制工頻干擾的方法

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

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

0 引言

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

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

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

1 半航空數字諧振器設計

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

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

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

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

其系統函數為：

(1)

可以看出，此系統的幅度特性在ω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)

該系統在任意頻點的幅度特性可以寫成式(4)。

(4)

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

(5)

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

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

(6)

因此，當

(7)

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

Δω≈(1-r)

(8)

通過以上諧振器的系統函數分析可知，設計一個數字諧振器的主要步驟如下：

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

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

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

3)根據式(1)設計諧振器。

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

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

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

2 模擬信號消噪效果

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

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

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

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

(9)

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

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

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

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

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

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

3 野外實測數據消噪

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

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

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

圖8 工區測線布置圖Fig.8 Work area survey line layou

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

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

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

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

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

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

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

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

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

圖12 地質解釋圖Fig.12 Geological interpretation graphics

4 結論

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

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

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