一种基于数字谐振器重构半航空瞬变电磁信号压制工频干扰的方法

郭 明, 王仕兴, 易国财, 何 可, 张振雄, 王绪本

(成都理工大学 地球勘探与信息技术教育部重点实验室，成都 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)将本文的数字谐振器技术应用于沱江半航空瞬变电磁实测数据，有效压制了工频干扰。去噪后的二次场衰减曲线用于反演后，得到的电阻率剖面更符合地质情况。