基于數(shù)字遞歸陷波的多通道瞬變電磁法周期噪聲去除研究

張文偉，底青云，耿啟立，雷達(dá)，王中興，繆佳佳

(1.中地裝(北京)科學(xué)技術(shù)研究院有限公司，北京 100011; 2.中國(guó)地質(zhì)裝備集團(tuán)有限公司，北京 100102; 3.中國(guó)科學(xué)院地質(zhì)與地球物理研究所 中國(guó)科學(xué)院頁(yè)巖氣與地質(zhì)工程重點(diǎn)實(shí)驗(yàn)室，北京 100029; 4.中國(guó)科學(xué)院 地球科學(xué)研究院，北京 100029; 5.中國(guó)科學(xué)院大學(xué)，北京 100049)

0 引言

多通道瞬變電磁法(Multi-transient electromagnetic method, MTEM)由Wright和 Ziolkowski在2002年提出，同時(shí)測(cè)量發(fā)射電流和響應(yīng)電壓，通過反褶積得到大地脈沖響應(yīng)，利用大地脈沖響應(yīng)進(jìn)行進(jìn)一步數(shù)據(jù)處理[1-3]。多通道瞬變電磁法可以直接探測(cè)油氣與監(jiān)測(cè)儲(chǔ)藏，在陸地與淺海油氣探測(cè)中均取得了成功應(yīng)用[4-6]。近年來，國(guó)內(nèi)的許多學(xué)者從方法原理、數(shù)據(jù)處理與解釋、儀器研制與系統(tǒng)集成應(yīng)用等多個(gè)方面對(duì)多通道瞬變電磁法進(jìn)行了深入研究[7-13]。

本文開展多通道瞬變電磁法反褶積后陷波研究，設(shè)計(jì)數(shù)字遞歸陷波器去除大地脈沖響應(yīng)包含的周期噪聲。首先，對(duì)傳統(tǒng)的零點(diǎn)極點(diǎn)設(shè)計(jì)數(shù)字遞歸陷波器的方法進(jìn)行了修正，采用陷波帶寬作為設(shè)計(jì)參數(shù)，其物理意義明確，便于分析問題。文中推導(dǎo)了基于陷波帶寬計(jì)算數(shù)字遞歸陷波器濾波系數(shù)的公式，并與其他設(shè)計(jì)方法進(jìn)行了對(duì)比，說明修正的設(shè)計(jì)方法正確。其次，對(duì)數(shù)字遞歸陷波器的特征屬性進(jìn)行了簡(jiǎn)要闡述。同時(shí)，系統(tǒng)地分析了陷波帶寬、初始條件、及大地脈沖響應(yīng)本身的波形對(duì)數(shù)字遞歸陷波器的瞬態(tài)響應(yīng)的壓制效果及陷波結(jié)果的影響。針對(duì)第三類初始條件，討論了初始值個(gè)數(shù)的選擇對(duì)陷波結(jié)果的影響。最后，應(yīng)用反褶積后數(shù)字遞歸陷波方法去除實(shí)測(cè)數(shù)據(jù)計(jì)算的大地脈沖響應(yīng)中的工頻噪聲，效果良好。

1 數(shù)字遞歸陷波器設(shè)計(jì)

1.1 陷波器基本原理

理想單頻點(diǎn)固定型陷波器幅頻響應(yīng)為

(1)

其中：Ω為頻率，Ω0為陷波頻率，該頻率幅值響應(yīng)為0，其他頻率的幅值響應(yīng)為1，特別地，直流幅值與Nyquist頻率幅值均為1。理想的陷波器不可實(shí)現(xiàn)，實(shí)際陷波器分為通帶、阻帶與過渡區(qū)三個(gè)頻率范圍。通帶內(nèi)幅值基本沒有衰減，幅值近似為1；阻帶范圍內(nèi)幅值衰減非常大，阻帶內(nèi)最大的衰減，即陷波頻率的幅值稱為陷波器的深度；通帶與阻帶之間稱為過渡帶，幅值衰減介于二者之間。一個(gè)二階模擬陷波器的系統(tǒng)函數(shù)[21-22]為

(2)

式中b為陷波帶寬。如圖1，陷波帶寬是指在該頻率范圍內(nèi)功率衰減大于3 dB，表征陷波器的陷波性能。陷波帶寬越窄，陷波頻率之外的頻率受影響越小，因此，理論上陷波帶寬越窄越好。可以驗(yàn)證，式(2)中當(dāng)s=jΩ0時(shí)，H(s)=0，陷波頻率的幅值為零，此時(shí)，陷波器的深度為無窮大。

圖1 陷波器原理示意Fig.1 Sketch of notcher

本文主要討論數(shù)字遞歸陷波器，數(shù)字遞歸陷波器是對(duì)式(2)的離散實(shí)現(xiàn)。數(shù)字陷波器的陷波帶寬與陷波頻率，分別是用采樣率對(duì)模擬陷波頻率與陷波帶寬進(jìn)行規(guī)一化的結(jié)果。數(shù)字域頻率通常記為ω，對(duì)應(yīng)的數(shù)字域陷波頻率記為ω0，幅頻響應(yīng)為

(3)

數(shù)字域的陷波帶寬仍然表示陷波器功率衰減大于3 dB的頻率范圍，下文記為bw。事實(shí)上，上述數(shù)字域量均無量綱，而對(duì)應(yīng)模擬量有量綱，易于區(qū)分。因此，為了表述簡(jiǎn)潔，在不引起誤解的情況下，下文將數(shù)字域頻率簡(jiǎn)稱為頻率，數(shù)字域陷波頻率簡(jiǎn)稱為陷波頻率，數(shù)字域陷波帶寬簡(jiǎn)稱為陷波帶寬。

1.2 修正零點(diǎn)極點(diǎn)法

設(shè)計(jì)數(shù)字遞歸陷波器常用方法有雙線性變換法[21-22]、零點(diǎn)極點(diǎn)法[18]、全通濾波器法[23]。零點(diǎn)極點(diǎn)法思路與實(shí)現(xiàn)都非常簡(jiǎn)單，但是Strack的設(shè)計(jì)方法中計(jì)算濾波系數(shù)的參數(shù)物理意義不夠直觀，不便于參數(shù)選擇及問題分析。本小節(jié)先給出Strack的設(shè)計(jì)方法，然后推導(dǎo)基于陷波帶寬的修正零點(diǎn)極點(diǎn)法計(jì)算濾波系數(shù)的公式。

Strack設(shè)計(jì)方法直接考慮系統(tǒng)函數(shù)零點(diǎn)與極點(diǎn)的分布。首先，考慮只有一階零點(diǎn)的陷波器，零點(diǎn)對(duì)應(yīng)陷波頻率，其系統(tǒng)函數(shù)為[18]

(4)

令z=ejω，可以得到頻率響應(yīng)：

(5)

圖2 二階數(shù)字遞歸陷波器幾何解釋Fig.2 Geometric interpretation of 2nd-order digital recursive notcher

(6)

如果極點(diǎn)P與零點(diǎn)B足夠近，則當(dāng)頻點(diǎn)A不是在零點(diǎn)附近，幅值近似為1，當(dāng)頻點(diǎn)A為零點(diǎn)時(shí)，則幅值為0，這恰好就是所需要的陷波器。需要注意的是，為了使系統(tǒng)穩(wěn)定，極點(diǎn)P必須在單位圓內(nèi)，因此極點(diǎn)P在單位圓內(nèi)且離零點(diǎn)B足夠近。

考慮兩個(gè)特殊頻點(diǎn)，直流(1，0)與Nyquist頻率(-1，0)，二者的幅值響應(yīng)相等，均近似為1，有

(7)

式(7)左端表示直流(1，0)的幅值響應(yīng)，右端表示Nyquist頻率(-1，0)的幅值響應(yīng)。如果設(shè)zn=(α,β)，zp=(x,y)，并設(shè)

x=γα

(8)

代入式(7)，則有

(9)

由于零點(diǎn)zn=(α,β)在單位圓上，從而有α2+β2=1，結(jié)合式(8)、(9)可以得到

x2+y2=2γ-1 。

(10)

將式(8)、(9)、(10)代入式(7)右端項(xiàng)，得到

(11)

實(shí)際應(yīng)用的陷波器通常為實(shí)系數(shù)陷波器，因此零點(diǎn)與極點(diǎn)均以共軛對(duì)出現(xiàn)，系統(tǒng)函數(shù)式(6)進(jìn)一步改寫成

(12)

α=cosω0。

(13)

陷波器按式(12)設(shè)計(jì)只有一個(gè)待定參數(shù)γ，但是，其物理意義并不直觀。為此，考慮頻率響應(yīng)，將z=ejω代入式(12)，整理得

(14)

(15)

(cosω-α)2γ2=(1-γ)2sin2ω。

(16)

(17)

利用和差化積公式，有

(18)

結(jié)合式(13，17，18)，則陷波帶寬(記為bw)為

(19)

由式(11)知γ<1，所以

(20)

結(jié)合式(12)、(13)、(20)，則可以根據(jù)陷波頻率ω0及陷波帶寬bw，得到陷波器的系統(tǒng)函數(shù)與濾波系數(shù)。為了便于后面分析，將式(12)改寫成通用二階濾波器形式：

(21)

其中，a1、a2與b1由陷波頻率與陷波帶寬決定：

(22)

(23)

(24)

通過以上推導(dǎo)，修正的數(shù)字遞歸陷波器設(shè)計(jì)方法(式(21)～式(24))可以直接由陷波帶寬計(jì)算濾波系數(shù)。陷波帶寬有明確的物理意義，是指在該頻率范圍內(nèi)陷波器的功率衰減大于3 dB，表征了陷波器的陷波性能。理論上講，陷波帶寬越窄越好，這樣對(duì)信息損害小，對(duì)噪聲抑制作用強(qiáng)。下文將利用式(21)～式(24)設(shè)計(jì)數(shù)字遞歸陷波器，并使用陷波帶寬分析相關(guān)問題。

1.3 數(shù)字遞歸陷波器設(shè)計(jì)實(shí)例

如果已知陷波頻率為0.25π，陷波帶寬為 0.02π，根據(jù)式(21)～式(24)計(jì)算濾波系數(shù)，并與雙線變換(BL)及全通濾波法(AP)進(jìn)行對(duì)比(表 1)，可以看出本文修正的零點(diǎn)極點(diǎn)法與雙線性變換得到的濾波系數(shù)完全一致。圖3為三種方法得到的陷波器頻譜響應(yīng)(橫軸是歸一化數(shù)字域頻率)，可見三者頻譜響應(yīng)基本重合。

表1 二階數(shù)字遞歸陷波器濾波系數(shù)Table 1 Coefficients of 2nd-order digital recursive notcher

2 數(shù)字遞歸陷波器的特征屬性

2.1 陷波頻率和陷波帶寬

根據(jù)式(21)可以得到時(shí)間域數(shù)字遞歸陷波器：

yn=-a1yn-1-a2yn-2+b1xn+b2xn-1+b1xn-2，

(25)

由第1節(jié)知，式(25)中的濾波系數(shù)由陷波頻率、陷波帶寬決定，陷波頻率與陷波帶寬是陷波器的首要特征屬性。

2.2 初始條件

仔細(xì)觀察式(25)可以發(fā)現(xiàn)，輸出yn不僅與輸入有關(guān)，還與之前的輸出yn-1,yn-2有關(guān)。特別地，y0與y-1,y-2有關(guān)。y-1,y-2表示零時(shí)刻之前陷波器的輸出，稱為初始條件。盡管初始條件的選擇是任意的，但是初始條件對(duì)陷波結(jié)果影響巨大。表2總結(jié)了三類常用的初始條件，其中第一類初始條件初始輸出信號(hào)為0；第二類初始條件采用輸入信號(hào)作為初始輸出信號(hào)，即Strack[18]使用的初始條件；第三類初始條件采用正交投影算子[20]估計(jì)輸入包含的真實(shí)信號(hào)分量(表2中s0與s1)作為初始輸出。如果y0,y1,…,yM-1均采用第三類初始條件估計(jì)得到，則稱初始值的個(gè)數(shù)為M。表2為M=2的特殊情況，下文在沒有特殊說明的情況下，M默認(rèn)等于2。與前兩類不同，第三類初始條件直接估計(jì)真實(shí)信號(hào)分量作為初始輸出，估計(jì)的信號(hào)分量是有效數(shù)據(jù)，因此初始值的個(gè)數(shù)不會(huì)減少有效數(shù)據(jù)長(zhǎng)度。

圖3 數(shù)字遞歸陷波器的頻譜Fig.3 Spectrum of digital recursive notcher

表2 二階數(shù)字遞歸陷波器初始條件Table 2 Initial conditions of 2nd-order digital recursive notcher

2.3 瞬態(tài)響應(yīng)

瞬態(tài)響應(yīng)是指單頻點(diǎn)正弦波作為輸入(即輸入只包含噪聲，信號(hào)分量為零)時(shí)數(shù)字遞歸陷波器的輸出響應(yīng)。Pei等[20]與Piskorowski[24-25]討論了初始條件對(duì)瞬態(tài)響應(yīng)的影響，并分析了瞬態(tài)響應(yīng)對(duì)心電圖噪聲去除的破壞性。大地脈沖響應(yīng)與心電脈沖都是瞬態(tài)信號(hào)，在受瞬態(tài)響應(yīng)影響這個(gè)方面應(yīng)具有相似性，因此，本文主要研究大地脈沖響應(yīng)噪聲去除受瞬態(tài)響應(yīng)的影響。

2.4 相位

數(shù)字遞歸陷波器的相位響應(yīng)是非線性的，如果用于去除大地脈沖響應(yīng)信號(hào)的周期噪聲，會(huì)造成時(shí)間延遲與相位畸變。為此，本文采用正反兩次濾波的方式消除非線性相位[18, 26]。

3 數(shù)字遞歸陷波效果的影響因素

3.1 陷波帶寬及初始條件對(duì)瞬態(tài)響應(yīng)的影響

事實(shí)上，瞬態(tài)響應(yīng)除了受初始條件影響，還受陷波帶寬影響。下面采用合成數(shù)據(jù)，對(duì)比不同陷波帶寬，不同初始條件的瞬態(tài)響應(yīng)。輸入信號(hào)為50 Hz的正弦波，采樣率為16 384 Hz，因此，數(shù)字域陷波頻率為ω0=50×2π/16 384。分別選取3個(gè)陷波帶寬(ω0/2、ω0/10、ω0/50)應(yīng)用三類初始條件，得到對(duì)應(yīng)的瞬態(tài)響應(yīng)結(jié)果如圖4所示。圖4a為第一類初始條件下，數(shù)字遞歸陷波器的瞬態(tài)響應(yīng)，當(dāng)陷波帶寬為ω0/2時(shí)，瞬態(tài)響應(yīng)在約0.1 s時(shí)衰減殆盡；當(dāng)陷波帶寬增加到ω0/10，瞬態(tài)響應(yīng)在約0.3 s時(shí)衰減殆盡；當(dāng)陷波帶寬進(jìn)一步增加到ω0/50，瞬態(tài)響應(yīng)需要更長(zhǎng)時(shí)間才能衰減殆盡。圖4b為第二類初始條件下，不同陷波帶寬條件下，數(shù)字遞歸陷波器的瞬態(tài)響應(yīng)，規(guī)律與第一類相同，隨著陷波帶寬的減小，瞬態(tài)響應(yīng)衰減時(shí)間越長(zhǎng)。圖4c為第三類初始條件的數(shù)字遞歸陷波器的瞬態(tài)響應(yīng)，可以看出，盡管隨著陷波帶寬減小瞬態(tài)響應(yīng)衰減時(shí)間變長(zhǎng)，但是瞬態(tài)響應(yīng)顯著受到抑制，開始時(shí)刻就明顯小于前兩類初始條件，即使陷波帶寬達(dá)到ω0/50，瞬態(tài)響應(yīng)也比較小，并且能夠快速衰減。

由以上分析可看出，陷波帶寬對(duì)瞬態(tài)響應(yīng)影響較大，隨著陷波帶寬的減小，瞬態(tài)響應(yīng)增強(qiáng)且衰減時(shí)間變長(zhǎng)。同時(shí)，不同的初始條件對(duì)瞬態(tài)響應(yīng)的壓制效果有明顯差異，在相同陷波帶寬條件下，第一類與第二類初始條件的瞬態(tài)響應(yīng)明顯比第三類初始條件強(qiáng)，且衰減時(shí)間更長(zhǎng)，第三類初始條件對(duì)數(shù)字遞歸陷波器的瞬態(tài)響應(yīng)抑制作用明顯，這與Pei等[20]的結(jié)果一致。

3.2 不同初始條件下大地脈沖響應(yīng)周期噪聲的去除效果

從前節(jié)可以看出， 數(shù)字遞歸陷波器本身的瞬態(tài)響應(yīng)需要經(jīng)過一段時(shí)間才能衰減殆盡。大地脈沖響應(yīng)是一種瞬態(tài)信號(hào)，通常持續(xù)時(shí)間為幾毫秒到幾秒，因此，可以推測(cè)大地脈沖響應(yīng)中周期噪聲去除效果受瞬態(tài)響應(yīng)影響尤其明顯。本節(jié)分析初始條件對(duì)大地脈沖響應(yīng)周期噪聲去除的影響。

a—第一類初始條件；b—第二類初始條件；c—第三類初始條件a—first kind of initial condition；b—second kind of initial condition； c—third kind of initial condition圖4 數(shù)字遞歸陷波器瞬態(tài)響應(yīng)Fig.4 Transient response of digital recursive notcher

圖5中所加入的正弦噪聲為50 Hz，最大值與脈沖響應(yīng)最大值相等，初始相位為π/4；采樣率為16 384Hz。對(duì)測(cè)試信號(hào)進(jìn)行數(shù)字遞歸陷波，陷波頻率為ω0=2π×50/16 384，使用3種陷波帶寬：ω0/2，ω0/5，ω0/50，分別采用三類初始條件，結(jié)果分別如圖6所示。第一類初始條件(圖6a)，當(dāng)陷波帶寬為ω0/2時(shí)，正弦噪聲受到很大程序的抑制，但是殘余量仍然較大，尤其在脈沖響應(yīng)的峰值附近較強(qiáng)，到0.05 s后，正弦噪聲基本被消除;減小陷波帶寬(ω0/10)，陷波效果反而變差，這是由于陷波帶寬變小引起瞬態(tài)響應(yīng)增強(qiáng);當(dāng)陷波帶寬減小到ω0/50時(shí)，陷波器近乎失效。第二類初始條件(圖6b)，當(dāng)陷波帶寬為ω0/2，陷波結(jié)果略優(yōu)于第一類初始條件，正弦噪聲在大地脈沖響應(yīng)的峰值附近仍有大量殘余，減小陷波帶寬，與第一類初始條件類似，由于瞬態(tài)響應(yīng)增強(qiáng)，陷波不能恢復(fù)原始信號(hào)。第三類初始條件(圖6c)，當(dāng)陷波帶寬為ω0/2時(shí)，正弦噪聲也有大量殘余,當(dāng)減小陷波帶寬時(shí)，正弦噪聲殘余量減小，當(dāng)陷波帶寬為ω0/10時(shí)，噪聲基本得到壓制；進(jìn)一步減小陷波帶寬至ω0/50，陷波結(jié)果與原始信號(hào)完全重合，噪聲被近乎完全壓制,這是由于陷波帶寬減小，陷波器抑制噪聲能力增強(qiáng)；同時(shí)，第三類初始條件能夠有效壓制陷波器本身的瞬態(tài)響應(yīng)，從而可以最大限度地消除正弦噪聲。

以上分析表明，對(duì)于大地脈沖響應(yīng)這類瞬態(tài)信號(hào)，陷波帶寬較大時(shí)，正弦噪聲殘余量大，對(duì)早期信號(hào)恢復(fù)影響較大。減小陷波帶寬，如果初始條件選擇不當(dāng)，由于瞬態(tài)響應(yīng)的原因，陷波結(jié)果反而變差。

圖5 電阻率25 Ω·m均勻半空間偏移距1 000 m處大地脈沖響應(yīng)及添加50 Hz正弦噪聲后信號(hào)Fig.5 Earth impulse response of 25 Ω·m homogenous half space at an offset of 1 000 m and contaminated signal by 50 Hz sinusoidal noise

a—第一類初始條件；b—第二類初始條件；c—第三類初始條件a—first kind of initial condition；b—second kind of initial condition； c—third kind of initial condition圖6 大地脈沖響應(yīng)周期噪聲去除受瞬態(tài)響應(yīng)的影響Fig.6 Influence of transient response on removal of periodical noise in earth impulse response

第三類初始條件能夠有效壓制數(shù)字遞歸陷波器的瞬態(tài)響應(yīng)，從而可以實(shí)現(xiàn)窄帶陷波，正弦噪聲能夠最大限度被消除。

3.3 第三類初始條件初始值個(gè)數(shù)對(duì)陷波結(jié)果的影響

看另一個(gè)例子。如圖7所示，考慮兩類典型的非均勻半空間大地脈沖響應(yīng)，與由解析解求得的均勻半空間大地脈沖響應(yīng)不同，由于數(shù)值計(jì)算本身的局限性，這兩個(gè)大地脈沖響應(yīng)的早期信號(hào)都有缺失，起始時(shí)刻的信號(hào)值不為零，圖7a中的大地脈沖響應(yīng)包含空氣波，圖7b中的大地脈沖響應(yīng)不含空氣波。采用50 Hz正弦信號(hào)模擬工頻噪聲，添加噪聲后的信號(hào)分別如圖7a與7b的虛線所示。采樣率為16 000 Hz，陷波頻率為ω0=2π×50/16 000，陷波帶寬依次取ω0/2，ω0/5，ω0/50。圖7c～7e顯示了利用數(shù)字遞歸陷波器去除圖7a中含噪大地脈沖響應(yīng)中的正弦噪聲的結(jié)果；圖7f～7h為圖7b中含噪大地脈沖響應(yīng)的陷波結(jié)果。從圖中可以看出，兩個(gè)信號(hào)的陷波結(jié)果顯示的規(guī)律一致，當(dāng)使用第一類(圖7c與圖7f)與第二類(圖7d與圖7g)初始條件時(shí)，瞬態(tài)響應(yīng)未得到有效壓制，減小陷波帶寬， 瞬態(tài)響應(yīng)增強(qiáng)，陷波效果反而變差，這與前文結(jié)果一致。不同之處，第三類(圖7c與圖7h)初始條件也未能有效消除正弦噪聲。

進(jìn)一步研究發(fā)現(xiàn)，應(yīng)用第三類初始條件的陷波結(jié)果受初始值的個(gè)數(shù)影響顯著。為此，固定陷波帶寬為ω0/50，初始值個(gè)數(shù)分別取不同的值，陷波結(jié)果展示在圖8中。圖8a中，當(dāng)初始值個(gè)數(shù)為20時(shí)，瞬態(tài)響應(yīng)仍然沒有得到壓制；當(dāng)初始值個(gè)數(shù)增加到200個(gè)，瞬態(tài)響應(yīng)基本被抑制，正弦噪聲得到消除。繼續(xù)增加初始值個(gè)數(shù)，陷波結(jié)果變化不顯著。圖8b中，當(dāng)初始值個(gè)數(shù)為20時(shí)，瞬態(tài)響應(yīng)也非常強(qiáng)，正弦噪聲殘余量大；當(dāng)初始值個(gè)數(shù)增加到200時(shí)，仍然有少量噪聲殘留；進(jìn)一步增加初始值的個(gè)數(shù)，正弦噪聲基本被消除。

a—包含空氣波的大地脈沖響應(yīng)與含噪信號(hào)；b—不包含空氣波的大地脈沖響應(yīng)與含噪信號(hào)；c～e——為(a)中含噪信號(hào)的陷波結(jié)果；f～h—為(b)中含噪信號(hào)的陷波結(jié)果a—impulse response with air wave and contaminated signal by 50 Hz sinusoidal noise; b—impulse response without air wave and contaminated signal by 50 Hz sinusoidal noise; c～e— are notch result of noised signal in (a); f～h— are notch result of noised signal in (b) 圖7 非均勻半空間大地脈沖響應(yīng)周期噪聲去除Fig.7 Removal of periodic noise in impulse response of inhomogeneous earth

a—圖7a中含噪大地脈沖響應(yīng)的陷波結(jié)果；b—圖7b中含噪大地脈沖響應(yīng)的陷波結(jié)果a—notch result of the noised earth impulse response in figure 7a; b—notch result of the noised earth impulse response in figure 7b圖8 第三類初始條件初始值個(gè)數(shù)對(duì)陷波結(jié)果的影響Fig.8 Influence of different number of initial values of the third kind of initial condition on notch result

以上分析表明，大地脈沖響應(yīng)的波形對(duì)數(shù)字遞歸陷波器的瞬態(tài)響應(yīng)有影響。特別地，對(duì)于部分大地脈沖響應(yīng)， 第三類初始條件如果初始值個(gè)數(shù)選擇太少，瞬態(tài)響應(yīng)不能被抑制，正弦噪聲不能被消除；通過增大初始值的個(gè)數(shù)，瞬態(tài)響應(yīng)能夠被抑制，從而大地脈沖響應(yīng)中包含的周期噪聲能夠被有效去除。

4 野外實(shí)測(cè)數(shù)據(jù)

將數(shù)字遞歸陷波方法應(yīng)用于野外實(shí)測(cè)數(shù)據(jù)中。發(fā)射偶極子位于0～240 m，接收偶極分別位于1 160～1 200 m，偏移距為發(fā)射與接收偶極子中點(diǎn)距離，為1 060 m。發(fā)射電流±25 A、12階，基頻為512 Hz的m序列電流，發(fā)射電流實(shí)際波形如圖9a所示，接收電壓波形如圖9b所示。采用Wiener濾波法[27]估計(jì)的大地脈沖響應(yīng)如圖9c中的虛線所示。

從接收電壓和估計(jì)的大地脈沖響應(yīng)均可以看出，工頻干擾非常強(qiáng)，大地脈沖響應(yīng)完全被湮沒。采用式(21～24)設(shè)計(jì)數(shù)字遞歸陷波器，其中，陷波頻率ω0=2π×50/16 000，陷波帶寬取ω0/30，初始條件采用第三類初始條件，初始值的個(gè)數(shù)取1 600，對(duì)Wiener濾波估計(jì)的大地脈沖響應(yīng)(圖9c虛線)進(jìn)行陷波，結(jié)果如圖9c的實(shí)線所示?？梢钥闯?，周期工頻噪聲基本受到壓制。

為了進(jìn)一步驗(yàn)證本文方法的有效性，處理一組共發(fā)射點(diǎn)道集的實(shí)測(cè)數(shù)據(jù)，發(fā)射位于480～720 m處，發(fā)射電流為m序列，接收排列2 420～2 700 m，點(diǎn)距為40 m，共8道。首先，采用Wiener濾波估計(jì)大地脈沖響應(yīng)(如圖10中的虛線所示)；然后應(yīng)用數(shù)字遞歸陷波法去除大地脈沖響應(yīng)中周期噪聲(去噪后的脈沖響應(yīng)如圖10的實(shí)線所示)，陷波頻率為ω0=2π×50/16 000，陷波帶寬取ω0/30，初始條件采用第三類初始條件，初始值的個(gè)數(shù)取1 600。從圖中可以看出，大地脈沖響應(yīng)中周期噪聲能夠被完全壓制。

5 結(jié)論

數(shù)字遞歸陷波器的濾波系數(shù)可以由陷波頻率與陷波帶寬直接計(jì)算，修正的零點(diǎn)極點(diǎn)法設(shè)計(jì)的數(shù)字遞歸陷波器正確。陷波帶寬對(duì)瞬態(tài)響應(yīng)影響顯著，一般情況下，陷波帶寬越窄，瞬態(tài)響應(yīng)越強(qiáng)。如果初始條件選擇不恰當(dāng)，陷波帶寬越小，由于瞬態(tài)響應(yīng)的原因，陷波結(jié)果反而越差。第三類初始條件能顯著壓制瞬態(tài)響應(yīng)，可以實(shí)現(xiàn)窄帶陷波，更好地壓制周期噪聲。同時(shí)，大地脈沖響應(yīng)波形對(duì)陷波結(jié)果也有影響，當(dāng)大地脈沖響應(yīng)早期信息缺失而初始值不為零時(shí)，通過增加第三類初始條件初始值的個(gè)數(shù)，可以改善陷波結(jié)果。實(shí)測(cè)數(shù)據(jù)處理結(jié)果表明，數(shù)字遞歸陷波可以去除大地脈沖響應(yīng)的周期性工頻噪聲。

a——發(fā)射電流; b—接收電壓; c—直接計(jì)算的大地脈沖響應(yīng)(虛線)與去除工頻噪聲后的大地脈沖響應(yīng)(實(shí)線)a—source current; b—received voltage; c—direct estimated earth impulse response (dashed line) and denoised earth impulse response (solid line)圖9 MTEM野外實(shí)測(cè)數(shù)據(jù)周期噪聲去除Fig.9 Removal of periodic noise in MTEM real field data

虛線為計(jì)算的大地脈沖響應(yīng)，實(shí)線為去除工頻噪聲后的大地脈沖響應(yīng)dashed line is evaluated earth impulse response and solid line denoised earth impulse response 圖10 MTEM野外實(shí)測(cè)數(shù)據(jù)周期噪聲去除Fig.10 Removal of periodic noise in MTEM real field data

雖然理論上陷波帶寬越窄越好，但是，實(shí)測(cè)數(shù)據(jù)中包含的工頻噪聲的頻率并非絕對(duì)不變，加之有限字長(zhǎng)的精度影響，陷波帶寬不應(yīng)無限小。當(dāng)工頻為50 Hz，陷波帶寬為陷波頻率的1/50時(shí)，陷波帶寬對(duì)應(yīng)1 Hz，其物理意義表示49～51 Hz范圍內(nèi)陷波器功率衰減超過3 dB。處理實(shí)測(cè)數(shù)據(jù)時(shí)，應(yīng)根據(jù)需要選擇的合適陷波帶寬，一般情況下，1 Hz精度已經(jīng)能滿足要求。第三類初始條件中初始值個(gè)數(shù)的具體值可以通過嘗試選擇，當(dāng)繼續(xù)增加對(duì)陷波結(jié)果影響不大時(shí)，就沒有必要繼續(xù)增加。

致謝:感謝自然資源航空物探遙感中心欒曉東博士為本文提供非均勻半空間大地脈沖響應(yīng)正演數(shù)據(jù)，感謝中國(guó)科學(xué)院地質(zhì)與地球物理研究所真齊輝博士對(duì)本文提出的修改意見。特別感謝審稿專家及編輯對(duì)本文提出的寶貴修改意見。