基于数字递归陷波的多通道瞬变电磁法周期噪声去除研究

张文伟，底青云，耿启立，雷达，王中兴，缪佳佳

(1.中地装(北京)科学技术研究院有限公司，北京 100011; 2.中国地质装备集团有限公司，北京 100102; 3.中国科学院地质与地球物理研究所 中国科学院页岩气与地质工程重点实验室，北京 100029; 4.中国科学院 地球科学研究院，北京 100029; 5.中国科学院大学，北京 100049)

0 引言

多通道瞬变电磁法(Multi-transient electromagnetic method, MTEM)由Wright和 Ziolkowski在2002年提出，同时测量发射电流和响应电压，通过反褶积得到大地脉冲响应，利用大地脉冲响应进行进一步数据处理[1-3]。多通道瞬变电磁法可以直接探测油气与监测储藏，在陆地与浅海油气探测中均取得了成功应用[4-6]。近年来，国内的许多学者从方法原理、数据处理与解释、仪器研制与系统集成应用等多个方面对多通道瞬变电磁法进行了深入研究[7-13]。

本文开展多通道瞬变电磁法反褶积后陷波研究，设计数字递归陷波器去除大地脉冲响应包含的周期噪声。首先，对传统的零点极点设计数字递归陷波器的方法进行了修正，采用陷波带宽作为设计参数，其物理意义明确，便于分析问题。文中推导了基于陷波带宽计算数字递归陷波器滤波系数的公式，并与其他设计方法进行了对比，说明修正的设计方法正确。其次，对数字递归陷波器的特征属性进行了简要阐述。同时，系统地分析了陷波带宽、初始条件、及大地脉冲响应本身的波形对数字递归陷波器的瞬态响应的压制效果及陷波结果的影响。针对第三类初始条件，讨论了初始值个数的选择对陷波结果的影响。最后，应用反褶积后数字递归陷波方法去除实测数据计算的大地脉冲响应中的工频噪声，效果良好。

1 数字递归陷波器设计

1.1 陷波器基本原理

理想单频点固定型陷波器幅频响应为

(1)

其中：Ω为频率，Ω0为陷波频率，该频率幅值响应为0，其他频率的幅值响应为1，特别地，直流幅值与Nyquist频率幅值均为1。理想的陷波器不可实现，实际陷波器分为通带、阻带与过渡区三个频率范围。通带内幅值基本没有衰减，幅值近似为1；阻带范围内幅值衰减非常大，阻带内最大的衰减，即陷波频率的幅值称为陷波器的深度；通带与阻带之间称为过渡带，幅值衰减介于二者之间。一个二阶模拟陷波器的系统函数[21-22]为

(2)

式中b为陷波带宽。如图1，陷波带宽是指在该频率范围内功率衰减大于3 dB，表征陷波器的陷波性能。陷波带宽越窄，陷波频率之外的频率受影响越小，因此，理论上陷波带宽越窄越好。可以验证，式(2)中当s=jΩ0时，H(s)=0，陷波频率的幅值为零，此时，陷波器的深度为无穷大。

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

本文主要讨论数字递归陷波器，数字递归陷波器是对式(2)的离散实现。数字陷波器的陷波带宽与陷波频率，分别是用采样率对模拟陷波频率与陷波带宽进行规一化的结果。数字域频率通常记为ω，对应的数字域陷波频率记为ω0，幅频响应为

(3)

数字域的陷波带宽仍然表示陷波器功率衰减大于3 dB的频率范围，下文记为bw。事实上，上述数字域量均无量纲，而对应模拟量有量纲，易于区分。因此，为了表述简洁，在不引起误解的情况下，下文将数字域频率简称为频率，数字域陷波频率简称为陷波频率，数字域陷波带宽简称为陷波带宽。

1.2 修正零点极点法

设计数字递归陷波器常用方法有双线性变换法[21-22]、零点极点法[18]、全通滤波器法[23]。零点极点法思路与实现都非常简单，但是Strack的设计方法中计算滤波系数的参数物理意义不够直观，不便于参数选择及问题分析。本小节先给出Strack的设计方法，然后推导基于陷波带宽的修正零点极点法计算滤波系数的公式。

Strack设计方法直接考虑系统函数零点与极点的分布。首先，考虑只有一阶零点的陷波器，零点对应陷波频率，其系统函数为[18]

(4)

令z=ejω，可以得到频率响应：

(5)

图2 二阶数字递归陷波器几何解释Fig.2 Geometric interpretation of 2nd-order digital recursive notcher

(6)

如果极点P与零点B足够近，则当频点A不是在零点附近，幅值近似为1，当频点A为零点时，则幅值为0，这恰好就是所需要的陷波器。需要注意的是，为了使系统稳定，极点P必须在单位圆内，因此极点P在单位圆内且离零点B足够近。

考虑两个特殊频点，直流(1，0)与Nyquist频率(-1，0)，二者的幅值响应相等，均近似为1，有

(7)

式(7)左端表示直流(1，0)的幅值响应，右端表示Nyquist频率(-1，0)的幅值响应。如果设zn=(α,β)，zp=(x,y)，并设

x=γα

(8)

代入式(7)，则有

(9)

由于零点zn=(α,β)在单位圆上，从而有α2+β2=1，结合式(8)、(9)可以得到

x2+y2=2γ-1 。

(10)

将式(8)、(9)、(10)代入式(7)右端项，得到

(11)

实际应用的陷波器通常为实系数陷波器，因此零点与极点均以共轭对出现，系统函数式(6)进一步改写成

(12)

α=cosω0。

(13)

陷波器按式(12)设计只有一个待定参数γ，但是，其物理意义并不直观。为此，考虑频率响应，将z=ejω代入式(12)，整理得

(14)

(15)

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

(16)

(17)

利用和差化积公式，有

(18)

结合式(13，17，18)，则陷波带宽(记为bw)为

(19)

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

(20)

结合式(12)、(13)、(20)，则可以根据陷波频率ω0及陷波带宽bw，得到陷波器的系统函数与滤波系数。为了便于后面分析，将式(12)改写成通用二阶滤波器形式：

(21)

其中，a1、a2与b1由陷波频率与陷波带宽决定：

(22)

(23)

(24)

通过以上推导，修正的数字递归陷波器设计方法(式(21)～式(24))可以直接由陷波带宽计算滤波系数。陷波带宽有明确的物理意义，是指在该频率范围内陷波器的功率衰减大于3 dB，表征了陷波器的陷波性能。理论上讲，陷波带宽越窄越好，这样对信息损害小，对噪声抑制作用强。下文将利用式(21)～式(24)设计数字递归陷波器，并使用陷波带宽分析相关问题。

1.3 数字递归陷波器设计实例

如果已知陷波频率为0.25π，陷波带宽为 0.02π，根据式(21)～式(24)计算滤波系数，并与双线变换(BL)及全通滤波法(AP)进行对比(表 1)，可以看出本文修正的零点极点法与双线性变换得到的滤波系数完全一致。图3为三种方法得到的陷波器频谱响应(横轴是归一化数字域频率)，可见三者频谱响应基本重合。

表1 二阶数字递归陷波器滤波系数Table 1 Coefficients of 2nd-order digital recursive notcher

2 数字递归陷波器的特征属性

2.1 陷波频率和陷波带宽

根据式(21)可以得到时间域数字递归陷波器：

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

(25)

由第1节知，式(25)中的滤波系数由陷波频率、陷波带宽决定，陷波频率与陷波带宽是陷波器的首要特征属性。

2.2 初始条件

仔细观察式(25)可以发现，输出yn不仅与输入有关，还与之前的输出yn-1,yn-2有关。特别地，y0与y-1,y-2有关。y-1,y-2表示零时刻之前陷波器的输出，称为初始条件。尽管初始条件的选择是任意的，但是初始条件对陷波结果影响巨大。表2总结了三类常用的初始条件，其中第一类初始条件初始输出信号为0；第二类初始条件采用输入信号作为初始输出信号，即Strack[18]使用的初始条件；第三类初始条件采用正交投影算子[20]估计输入包含的真实信号分量(表2中s0与s1)作为初始输出。如果y0,y1,…,yM-1均采用第三类初始条件估计得到，则称初始值的个数为M。表2为M=2的特殊情况，下文在没有特殊说明的情况下，M默认等于2。与前两类不同，第三类初始条件直接估计真实信号分量作为初始输出，估计的信号分量是有效数据，因此初始值的个数不会减少有效数据长度。

图3 数字递归陷波器的频谱Fig.3 Spectrum of digital recursive notcher

表2 二阶数字递归陷波器初始条件Table 2 Initial conditions of 2nd-order digital recursive notcher

2.3 瞬态响应

瞬态响应是指单频点正弦波作为输入(即输入只包含噪声，信号分量为零)时数字递归陷波器的输出响应。Pei等[20]与Piskorowski[24-25]讨论了初始条件对瞬态响应的影响，并分析了瞬态响应对心电图噪声去除的破坏性。大地脉冲响应与心电脉冲都是瞬态信号，在受瞬态响应影响这个方面应具有相似性，因此，本文主要研究大地脉冲响应噪声去除受瞬态响应的影响。

2.4 相位

数字递归陷波器的相位响应是非线性的，如果用于去除大地脉冲响应信号的周期噪声，会造成时间延迟与相位畸变。为此，本文采用正反两次滤波的方式消除非线性相位[18, 26]。

3 数字递归陷波效果的影响因素

3.1 陷波带宽及初始条件对瞬态响应的影响

事实上，瞬态响应除了受初始条件影响，还受陷波带宽影响。下面采用合成数据，对比不同陷波带宽，不同初始条件的瞬态响应。输入信号为50 Hz的正弦波，采样率为16 384 Hz，因此，数字域陷波频率为ω0=50×2π/16 384。分别选取3个陷波带宽(ω0/2、ω0/10、ω0/50)应用三类初始条件，得到对应的瞬态响应结果如图4所示。图4a为第一类初始条件下，数字递归陷波器的瞬态响应，当陷波带宽为ω0/2时，瞬态响应在约0.1 s时衰减殆尽；当陷波带宽增加到ω0/10，瞬态响应在约0.3 s时衰减殆尽；当陷波带宽进一步增加到ω0/50，瞬态响应需要更长时间才能衰减殆尽。图4b为第二类初始条件下，不同陷波带宽条件下，数字递归陷波器的瞬态响应，规律与第一类相同，随着陷波带宽的减小，瞬态响应衰减时间越长。图4c为第三类初始条件的数字递归陷波器的瞬态响应，可以看出，尽管随着陷波带宽减小瞬态响应衰减时间变长，但是瞬态响应显著受到抑制，开始时刻就明显小于前两类初始条件，即使陷波带宽达到ω0/50，瞬态响应也比较小，并且能够快速衰减。

由以上分析可看出，陷波带宽对瞬态响应影响较大，随着陷波带宽的减小，瞬态响应增强且衰减时间变长。同时，不同的初始条件对瞬态响应的压制效果有明显差异，在相同陷波带宽条件下，第一类与第二类初始条件的瞬态响应明显比第三类初始条件强，且衰减时间更长，第三类初始条件对数字递归陷波器的瞬态响应抑制作用明显，这与Pei等[20]的结果一致。

3.2 不同初始条件下大地脉冲响应周期噪声的去除效果

从前节可以看出， 数字递归陷波器本身的瞬态响应需要经过一段时间才能衰减殆尽。大地脉冲响应是一种瞬态信号，通常持续时间为几毫秒到几秒，因此，可以推测大地脉冲响应中周期噪声去除效果受瞬态响应影响尤其明显。本节分析初始条件对大地脉冲响应周期噪声去除的影响。

a—第一类初始条件；b—第二类初始条件；c—第三类初始条件a—first kind of initial condition；b—second kind of initial condition； c—third kind of initial condition图4 数字递归陷波器瞬态响应Fig.4 Transient response of digital recursive notcher

图5中所加入的正弦噪声为50 Hz，最大值与脉冲响应最大值相等，初始相位为π/4；采样率为16 384Hz。对测试信号进行数字递归陷波，陷波频率为ω0=2π×50/16 384，使用3种陷波带宽：ω0/2，ω0/5，ω0/50，分别采用三类初始条件，结果分别如图6所示。第一类初始条件(图6a)，当陷波带宽为ω0/2时，正弦噪声受到很大程序的抑制，但是残余量仍然较大，尤其在脉冲响应的峰值附近较强，到0.05 s后，正弦噪声基本被消除;减小陷波带宽(ω0/10)，陷波效果反而变差，这是由于陷波带宽变小引起瞬态响应增强;当陷波带宽减小到ω0/50时，陷波器近乎失效。第二类初始条件(图6b)，当陷波带宽为ω0/2，陷波结果略优于第一类初始条件，正弦噪声在大地脉冲响应的峰值附近仍有大量残余，减小陷波带宽，与第一类初始条件类似，由于瞬态响应增强，陷波不能恢复原始信号。第三类初始条件(图6c)，当陷波带宽为ω0/2时，正弦噪声也有大量残余,当减小陷波带宽时，正弦噪声残余量减小，当陷波带宽为ω0/10时，噪声基本得到压制；进一步减小陷波带宽至ω0/50，陷波结果与原始信号完全重合，噪声被近乎完全压制,这是由于陷波带宽减小，陷波器抑制噪声能力增强；同时，第三类初始条件能够有效压制陷波器本身的瞬态响应，从而可以最大限度地消除正弦噪声。

以上分析表明，对于大地脉冲响应这类瞬态信号，陷波带宽较大时，正弦噪声残余量大，对早期信号恢复影响较大。减小陷波带宽，如果初始条件选择不当，由于瞬态响应的原因，陷波结果反而变差。

图5 电阻率25 Ω·m均匀半空间偏移距1 000 m处大地脉冲响应及添加50 Hz正弦噪声后信号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 大地脉冲响应周期噪声去除受瞬态响应的影响Fig.6 Influence of transient response on removal of periodical noise in earth impulse response

第三类初始条件能够有效压制数字递归陷波器的瞬态响应，从而可以实现窄带陷波，正弦噪声能够最大限度被消除。

3.3 第三类初始条件初始值个数对陷波结果的影响

看另一个例子。如图7所示，考虑两类典型的非均匀半空间大地脉冲响应，与由解析解求得的均匀半空间大地脉冲响应不同，由于数值计算本身的局限性，这两个大地脉冲响应的早期信号都有缺失，起始时刻的信号值不为零，图7a中的大地脉冲响应包含空气波，图7b中的大地脉冲响应不含空气波。采用50 Hz正弦信号模拟工频噪声，添加噪声后的信号分别如图7a与7b的虚线所示。采样率为16 000 Hz，陷波频率为ω0=2π×50/16 000，陷波带宽依次取ω0/2，ω0/5，ω0/50。图7c～7e显示了利用数字递归陷波器去除图7a中含噪大地脉冲响应中的正弦噪声的结果；图7f～7h为图7b中含噪大地脉冲响应的陷波结果。从图中可以看出，两个信号的陷波结果显示的规律一致，当使用第一类(图7c与图7f)与第二类(图7d与图7g)初始条件时，瞬态响应未得到有效压制，减小陷波带宽， 瞬态响应增强，陷波效果反而变差，这与前文结果一致。不同之处，第三类(图7c与图7h)初始条件也未能有效消除正弦噪声。

进一步研究发现，应用第三类初始条件的陷波结果受初始值的个数影响显著。为此，固定陷波带宽为ω0/50，初始值个数分别取不同的值，陷波结果展示在图8中。图8a中，当初始值个数为20时，瞬态响应仍然没有得到压制；当初始值个数增加到200个，瞬态响应基本被抑制，正弦噪声得到消除。继续增加初始值个数，陷波结果变化不显著。图8b中，当初始值个数为20时，瞬态响应也非常强，正弦噪声残余量大；当初始值个数增加到200时，仍然有少量噪声残留；进一步增加初始值的个数，正弦噪声基本被消除。

a—包含空气波的大地脉冲响应与含噪信号；b—不包含空气波的大地脉冲响应与含噪信号；c～e——为(a)中含噪信号的陷波结果；f～h—为(b)中含噪信号的陷波结果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 非均匀半空间大地脉冲响应周期噪声去除Fig.7 Removal of periodic noise in impulse response of inhomogeneous earth

a—图7a中含噪大地脉冲响应的陷波结果；b—图7b中含噪大地脉冲响应的陷波结果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 第三类初始条件初始值个数对陷波结果的影响Fig.8 Influence of different number of initial values of the third kind of initial condition on notch result

以上分析表明，大地脉冲响应的波形对数字递归陷波器的瞬态响应有影响。特别地，对于部分大地脉冲响应， 第三类初始条件如果初始值个数选择太少，瞬态响应不能被抑制，正弦噪声不能被消除；通过增大初始值的个数，瞬态响应能够被抑制，从而大地脉冲响应中包含的周期噪声能够被有效去除。

4 野外实测数据

将数字递归陷波方法应用于野外实测数据中。发射偶极子位于0～240 m，接收偶极分别位于1 160～1 200 m，偏移距为发射与接收偶极子中点距离，为1 060 m。发射电流±25 A、12阶，基频为512 Hz的m序列电流，发射电流实际波形如图9a所示，接收电压波形如图9b所示。采用Wiener滤波法[27]估计的大地脉冲响应如图9c中的虚线所示。

从接收电压和估计的大地脉冲响应均可以看出，工频干扰非常强，大地脉冲响应完全被湮没。采用式(21～24)设计数字递归陷波器，其中，陷波频率ω0=2π×50/16 000，陷波带宽取ω0/30，初始条件采用第三类初始条件，初始值的个数取1 600，对Wiener滤波估计的大地脉冲响应(图9c虚线)进行陷波，结果如图9c的实线所示。可以看出，周期工频噪声基本受到压制。

为了进一步验证本文方法的有效性，处理一组共发射点道集的实测数据，发射位于480～720 m处，发射电流为m序列，接收排列2 420～2 700 m，点距为40 m，共8道。首先，采用Wiener滤波估计大地脉冲响应(如图10中的虚线所示)；然后应用数字递归陷波法去除大地脉冲响应中周期噪声(去噪后的脉冲响应如图10的实线所示)，陷波频率为ω0=2π×50/16 000，陷波带宽取ω0/30，初始条件采用第三类初始条件，初始值的个数取1 600。从图中可以看出，大地脉冲响应中周期噪声能够被完全压制。

5 结论

数字递归陷波器的滤波系数可以由陷波频率与陷波带宽直接计算，修正的零点极点法设计的数字递归陷波器正确。陷波带宽对瞬态响应影响显著，一般情况下，陷波带宽越窄，瞬态响应越强。如果初始条件选择不恰当，陷波带宽越小，由于瞬态响应的原因，陷波结果反而越差。第三类初始条件能显著压制瞬态响应，可以实现窄带陷波，更好地压制周期噪声。同时，大地脉冲响应波形对陷波结果也有影响，当大地脉冲响应早期信息缺失而初始值不为零时，通过增加第三类初始条件初始值的个数，可以改善陷波结果。实测数据处理结果表明，数字递归陷波可以去除大地脉冲响应的周期性工频噪声。

a——发射电流; b—接收电压; c—直接计算的大地脉冲响应(虚线)与去除工频噪声后的大地脉冲响应(实线)a—source current; b—received voltage; c—direct estimated earth impulse response (dashed line) and denoised earth impulse response (solid line)图9 MTEM野外实测数据周期噪声去除Fig.9 Removal of periodic noise in MTEM real field data

虚线为计算的大地脉冲响应，实线为去除工频噪声后的大地脉冲响应dashed line is evaluated earth impulse response and solid line denoised earth impulse response 图10 MTEM野外实测数据周期噪声去除Fig.10 Removal of periodic noise in MTEM real field data

虽然理论上陷波带宽越窄越好，但是，实测数据中包含的工频噪声的频率并非绝对不变，加之有限字长的精度影响，陷波带宽不应无限小。当工频为50 Hz，陷波带宽为陷波频率的1/50时，陷波带宽对应1 Hz，其物理意义表示49～51 Hz范围内陷波器功率衰减超过3 dB。处理实测数据时，应根据需要选择的合适陷波带宽，一般情况下，1 Hz精度已经能满足要求。第三类初始条件中初始值个数的具体值可以通过尝试选择，当继续增加对陷波结果影响不大时，就没有必要继续增加。

致谢:感谢自然资源航空物探遥感中心栾晓东博士为本文提供非均匀半空间大地脉冲响应正演数据，感谢中国科学院地质与地球物理研究所真齐辉博士对本文提出的修改意见。特别感谢审稿专家及编辑对本文提出的宝贵修改意见。