一种去除时间域航空电磁数据中工频干扰的方法与流程

1.本发明涉及时间域航空电磁数据处理技术领域，尤其涉及一种去除时间域航空电磁数据中工频干扰的方法。


背景技术：

2.时间域航空电磁信号在采集过程中易受到来自各方面的干扰，导致信号淹没在噪声中而无法正确解译出地下介质分布信息。电磁勘探测区内难免会存在高压输电和变电装置对采集的信号产生干扰，因此合理正确的对其进行去除是解译地下介质构造的重要前提。我国的电网标准频率是50hz，习惯上称为工频，在地球物理探测中将其对信号产生的干扰称为工频干扰。
3.传统的阻带滤波方法主要有两种，一种是在频率域中进行阻带滤波后再转回时间域，一种是在时间域中进行陷波滤波。前者在应对海量数据时运算速度较慢，后者只能对50hz频率干扰进行去除，而无法对其谐波进行压制。文件一：王保利.工频干扰信号的时频域压制方法研究[j].中国煤炭地质，2015,4(27)：73
‑
77。公开了时频域的阻带滤波方法对工频干扰进行压制，在其计算过程中，将时域信号利用傅里叶变换转到频率域，在频率域对数据同时进行多次窄带滤波，然后将滤波结果再转换到时间域，从而实现对工频干扰的去除。文件二：steven w smith.digital signal processing[m].2003.library of congress cataloging in publication.公开了陷波滤波器的方法对含有的工频干扰噪声进行去除。陷波滤波器是通过对输入信号在某一频率点快速衰减，从而达到对该信号频点进行抑制的一种滤波器，属于带阻滤波器的一种，其阻带非常狭窄。
[0004]
文件一的多次阻带滤波方法主体思路基于频率域设计，且该方法受限于计算速度不能应用于海量数据的实时处理。
[0005]
文件一的方案是采用了时
‑
频和频
‑
时变换的阻带滤波方法。其将数据中工频干扰信号表示为：
[0006][0007]
其中，a
i
、f
i
、φ
i
分别为第i个频率工频干扰信号的振幅、频率、相位，t为时间，n为单频干扰信号的个数。根据有效信号的变化情况，将单道信号变换到频率域，来提取时频域的信号
[0008][0009]
表示沿时间轴取中值，通过在频率域进行滤波后，通过反傅里叶变换得到去噪后的时域信号：
[0010][0011]
其中x(t)表示原始信号，表示去噪后的信号，ifft表示反傅里叶变换。如果数据量大，将数据量裁剪，分块进行计算。
[0012]
利用该方案对地震数据进行了工频干扰的去除，但其去除方法采用了时间域数据转换到频率域，对频率域数据进行工频干扰去除后，再将频率域数据转换到时间域中，该方法可以有效的压制工频干扰。但此种方法应对少量数据时速度尚可，当应对时间域航空电磁海量数据时，区段式的时
‑
频和频
‑
时转换导致了处理速度慢，且远不能达到实时处理的要求。
[0013]
文件二由于陷波器的幅频特性只对工频干扰的50hz频率噪声进行去除，不能对50hz谐频干扰进行去除，从而导致原数据曲线形态发生改变，破坏了原数据衰减规律；因此陷波滤波器在工频干扰去除方面，存在着不足，需要寻找其他更好的方法对工频干扰进行去除。
[0014]
陷波滤波器是通过对输入信号在某一频率点快速衰减，从而达到对该信号频点进行抑制的一种滤波器。陷波滤波器是无限冲击响应(iir)数字滤波器，属于带阻滤波器的一种，其阻带非常狭窄，理想的陷波滤波器的频率响应为：
[0015][0016]
其中，ω0是陷波频率。在陷波频率点的信号频率响应为零，其他频点的响应值为1。滤波器中的3db频带的取值影响陷波的阻带宽度，3db频带取值越小，陷波滤波器越接近公式(4)给出结果，其滤除陷波频点的效果越好。
[0017]
为此，我们利用陷波滤波器去除50hz工频干扰，给出z域中陷波滤波器的二阶传递函数：
[0018][0019]
其中k为滤波器的增益，由幅度响应来确定，r为陷波滤波器的极半径，ω0＝2πf0，f0是陷波频点，将系统转换到时域后，再进行逆变换可得到对时间域数据进行处理的方法：
[0020]
y(n)＝a0x(n) a1x(n
‑
1) a2x(n
‑
2) b1y(n
‑
1) b2y(n
‑
2)
ꢀꢀꢀ
(6)
[0021]
a0＝k，a1＝
‑
2kcosω0，a2＝k，b1＝2rcosω0b2＝
‑
r2ꢀꢀꢀ
(7)
[0022]
式中x(n)可视为时间序列输入信号，y(n)为输出信号。
[0023]
该方案中详细的给出了陷波滤波器在时间域数据中的应用，但陷波滤波器存在的问题在于仅能对50hz频率进行窄带滤波，不能对50hz谐波进行压制，因此不能完全达到去除工频干扰的目的。


技术实现要素：

[0024]
为了解决以上问题，本发明提供一种去除时间域航空电磁数据中工频干扰的方法，它提出了iir型梳状滤波器方法，是一种可快速、高效地对海量时间域航空电磁数据进行实时处理的工频干扰压制方法。
[0025]
本发明的技术方案是：提供一种去除时间域航空电磁数据中工频干扰的方法，包括以下步骤：
[0026]
首先，观察时间域航空电磁数据中是否含有工频干扰噪声；
[0027]
其次，设计对时间域航空电磁数据中工频干扰进行去除的滤波器，该滤波器需应
对海量数据时的计算速度和精度；所述滤波器为iir型梳状滤波器；
[0028]
最后，利用设计好的滤波器对时间域航空电磁数据进行滤波处理；滤波后，观察干扰去除效果，并抽取部分数据进行频谱对比分析。
[0029]
进一步的，掌握时间域航空电磁数据特点后，观察数据趋势，若数据整体趋势呈现正/余弦变化，变化的主频率周期为50hz。
[0030]
进一步的，所述iir型梳状滤波器通过3db带宽来限制阻带的宽窄度：其系统函数为：
[0031]
其中
[0032]
将其进行z变换为：
[0033]
y(n)＝a[x(n)
‑
x(n
‑
n)] by(n
‑
n)，
ꢀꢀꢀ
(9)
[0034]
进一步的，所述iir型梳状滤波器的输出信号不仅取决于输入信号x(n)和x(n
‑
n)，同时与输出信号y(n
‑
n)也相关。
[0035]
与现有技术相比，本发明的有益效果是：不仅可以高效的在时间域中直接去除工频干扰中50hz主频噪声，且能将其谐波噪声去除，去除后的数据保持了原有数据的幅值，并且保证了数据变化趋势的正确性，本发明的算法可应用于时间域航空电磁连续采样海量数据的实时处理中。
附图说明
[0036]
下面根据图进一步对本发明加以说明:
[0037]
图1是本发明航空时间域航空电磁原始数据中工频干扰噪声形态图；
[0038]
图2是本发明iir型梳状滤波器的信号流图；
[0039]
图3是本发明iir型梳状滤波器特性曲线图；
[0040]
图4是本发明iir型梳状滤波器滤波效果图；
[0041]
图5是本发明iir型梳状滤波器幅频特性图；
[0042]
图6是本发明iir型梳状滤波器滤波效果图；
[0043]
图7是本发明iir型梳状滤波器幅频特性图；
具体实施方式
[0044]
下面结合图对本发明作进一步详细的说明，需要说明的是，图仅用于解释本发明，是对本发明实施例的示意性说明，而不能理解为对本发明的限定。
[0045]
在本发明的描述中，需要说明的是，除非另有明确的规定和限定，术语“相连”、“连接”应做广义理解，例如，可以是固定连接，也可以是可拆卸连接，或一体地连接；可以是机械连接，也可以是电连接；可以是直接相连，也可以通过中间媒介间接相连。对于本领域的普通技术人员而言，可以具体情况理解上述术语在本发明中的具体含义。
[0046]
工频干扰对时间域航空电磁实测数据影响较大。工频干扰指的是电力系统运行过程中对信号引起的一种干扰，其利用电磁感应电流和接地漏电电流等途径，存在于交流高压输电、变电和用电设施及大型建筑物周边。工频干扰的频率根据不同国家和地区交流电频率而定，中国大陆地区为50hz及其谐波构成。时间域航空电磁勘探过程中测区内难免会
存在高压输电和变电装置，因此在响应中会出现50hz及其谐波振荡，导致二次场数据曲线形态中混入了正/余弦信号。在进行时间域航空电磁的数据处理中，工频干扰则被视为电磁信号的噪声予以去除。正确、精准地对工频干扰进行去除，可以提高数据质量，是正确反映地下电性构造的必要手段。
[0047]
梳状滤波器是由许多按一定频率间隔相同排列的通带和阻带，只允许某些特定频率范围的信号通过的一种滤波器。其特性曲线像梳子一样，故称为梳状滤波器。梳状滤波器不仅可以去除某一固定频率干扰，同时可去除该频率产生的k次谐波，即用来抑制周期性的噪声。梳状滤波器属于数字滤波器的一种。
[0048]
如图1、图2、图3、图4、图5、图6、图7所示，一种去除时间域航空电磁数据中工频干扰的方法，包括以下步骤：
[0049]
首先，观察时间域航空电磁数据中是否含有工频干扰噪声。
[0050]
如图1所示，掌握时间域航空电磁数据特点后，观察数据趋势，若数据整体趋势呈现正/余弦变化，变化的主频率周期为50hz。例如，数据基频为12.5hz，一个周期中含有4个周期的工频干扰，即12.5hz*4＝50hz，可以判断出数据的干扰是电塔等人文设施产生的工频干扰，黑色曲线给出了工频干扰的正/余弦变化特征。
[0051]
其次，设计对时间域航空电磁数据中工频干扰进行去除的滤波器，该滤波器需应对海量数据时的计算速度和精度；所述滤波器为iir型梳状滤波器；
[0052]
iir型梳状滤波器通过3db带宽来限制阻带的宽窄度，其系统函数为：
[0053]
其中
[0054]
将其进行z变换为：
[0055]
y(n)＝a[x(n)
‑
x(n
‑
n)] by(n
‑
n)，
ꢀꢀꢀ
(9)
[0056]
如图2所示，从公式(9)中可以明确的看出，iir型梳状滤波器的输出信号不仅取决于输入信号x(n)和x(n
‑
n)，同时与输出信号y(n
‑
n)也相关。
[0057]
如图3所示，通过对时间域函数特性曲线进行分析，图中曲线所示，可见其阻带可以很好的限制50hz及其k次谐波的通过，因此选择iir型梳状滤波器更为适合去除工频干扰。
[0058]
最后，利用设计好的滤波器对时间域航空电磁数据进行滤波处理。滤波后，观察干扰去除效果，并抽取部分数据进行频谱对比分析。
[0059]
实施例1
[0060]
如图4、图5所示，选取固定翼时间域航空电磁系统(iftem
‑
ii)在吉林省大安市龙沼镇地区的实测飞行数据进行研究。系统数据采样率为100khz，基频为12.5hz，发射波形采用半正弦波，供电时间为4.15ms，抽取一部分实测数据进行滤波展示。
[0061]
选取一定周期数据对其进行梳状滤波和幅频特性的展示。原始数据(original data)由于存在50hz工频干扰导致曲线形态发生异常，当进行梳状滤波后曲线形态回归正常，为了更好的理解滤波效果，对滤波后数据进行傅里叶变换，观察其频谱特征。
[0062]
可以看出滤波器不仅可以很好的滤除50hz噪声且较明显的滤除了50hz谐波产生的干扰，很好的保持了原有数据的幅值大小，因此利用iir型梳状滤波器可以达到预期的效果。
[0063]
实施例2
[0064]
如图6、图7所示，选取固定翼时间域航空电磁系统(iftem
‑
i)在吉林省哈尔滨市宾县水库测区的实测飞行数据进行研究。系统数据采样率为100khz，基频为25hz，发射波形采用半正弦波，供电时间为4.15ms，抽取一部分实测数据进行滤波展示。
[0065]
选取一个周期数据进行梳状滤波和幅频特性的展示。可以看出原始数据(original data)由于存在50hz工频干扰导致曲线形态发生异常，当进行梳状滤波后曲线形态回归正常，为了更好的展示滤波效果，对滤波后数据进行傅里叶变换，观察其频谱特征。可以看出滤波器不仅可以很好的滤除50hz噪声且较明显的滤除了50hz谐波产生的干扰，且很好的保持了原有数据的幅值大小，因此利用iir型梳状滤波器可以达到预期的效果。
[0066]
本发明可以直接去除时间域航空电磁数据中的50hz噪声；进行处理时数据不需要时
‑
频和频
‑
时转换；去除时间域航空电磁数据中的50hz谐波噪声，如100hz，150hz，200hz
…
n*50hz；去噪后保持原有数据曲线的正确变化趋势；去噪后保持原有数据幅频特性的幅值大小；可应用于海量数据的实时处理；本发明的算法可应用于时间域航空电磁连续采样海量数据的实时处理中。本发明的iir型梳状滤波器不仅可以高效的去除工频干扰中50hz主频噪声，且能将其谐波噪声去除，去除后的数据保持了原有数据的幅值，并且保证了数据变化趋势的正确性。
[0067]
以上所述为本发明的实施例，并不用于限制本发明，对于本领域的技术人员来说，本发明可以有各种改进和变化。凡在本发明的精神和原则之内所作的任何修改、等同替换、改进等均应含在本发明的权利要求范围之内。