一种音频大地电磁数据工频噪声干扰的滤波方法与流程

1.本发明属于地球物理勘探领域的数据处理技术，具体涉及一种音频大地电磁数据工频噪声干扰的滤波方法。


背景技术：

2.通常，音频大地电磁仪在野外工作时均是同时采集测点上的电场信号和磁场信号，采集的频率范围大致是1～10000hz，不同仪器采集的频率范围略有不同。在音频大地电磁仪采集的这个信号频段内，非常容易受人为干扰的影响，尤其是工业用电的影响。然而，随着科技的进步和人为活动范围的日渐扩大，即便是山区也往往存在着较强的人为干扰，特别是1万伏、3.5万伏、11万伏等高压输电线路，对天然电磁场信号干扰非常严重。这些人为干扰严重降低了50hz及其倍频段的数据质量，甚至导致在高压线附近区域采集的数据无法使用。目前，世界上最先进和成熟的电磁仪器，内部已经集成了软、硬两方面的工频陷波技术，在采集方法上也采用了降低噪声干扰的布设装置，如v8多功能电法仪，在数据采集过程中采用电磁分离和远参考方法，但在强干扰地区很难克服工作中高压线附近带来的干扰问题，从而造成在50hz及其倍频段频点附近视电阻率及相位值发生畸变，影响最终的反演结果。


技术实现要素：

3.基于此，有必要针对音频大地电磁数据受工频噪声干扰的问题，提供一种音频大地电磁数据工频噪声干扰的滤波方法。
4.为了实现上述目的，本发明提供如下技术方案：
5.一种音频大地电磁数据工频噪声干扰的滤波方法，包括如下步骤：对有工频噪声干扰的音频大地电磁数据作带阻滤波，得到无工频噪声干扰数据。
6.进一步地，所述滤波方法具体包括如下步骤：
7.(1)判断音频大地电磁数据的原始时间序列文件的信号道功率谱曲线是否存在工频噪声干扰，得到第一判断结果；
8.(2)所述第一判断结果为是，对存在工频噪声干扰的信号道作带阻滤波，得到滤波后的数据；
9.(3)将滤波后的数据保存在一维数组中；
10.(4)对将滤波后的数据进行格式处理，得到格式文件。
11.进一步地，所述步骤(1)具体包括如下步骤：读取采样率为150hz的原始时间序列文件中的数据，对数据的每个信号道均做功率谱分析，得到每个信号道的功率谱曲线，根据每个信号道的功率谱曲线判断所述信号道是否存在工频噪声干扰，如果信号道的功率谱曲线出现毛刺信号或增益异常，则判断该信号道存在工频噪声干扰；如果信号道的功率谱曲线是一条近似平行横轴且连续光滑的曲线，则判断该信号道不存在工频噪声干扰。
12.进一步地，所述步骤(2)具体包括如下步骤：设定凯瑟窗参数，调用matlab中的凯
瑟窗函数构建凯瑟窗winfun，winfun＝(kaiser(n,β))'；将二维数组存储区中要滤波的信号道数据与winfun作卷积运算，得到滤波后的数据，完成对存在工频噪声干扰的信号道的滤波。
13.进一步地，步骤(2)中以fir数字滤波器作带阻滤波。
14.进一步地，所述凯瑟窗参数的设定包括：定频率带阻滤波器阻带低频边缘ω(p)、带阻滤波器阻带高频边缘ω(s)和衰减系数as，所述ω(p)、ω(s)和as分别设定为40hz～50hz，50hz～60hz和0～100。
15.进一步地，所述步骤(2)第一判断结果为是，对存在工频噪声干扰的信号道作带阻滤波前，还包括如下步骤：读取采样率为2400hz的原始时间序列文件的每个信号道的数据并保存在二维数组存储区中对应的行数中。
16.进一步地，所述读取采样率为2400hz的原始时间序列文件的每个信号道的数据，具体包括如下步骤：跳过采样率为2400hz的原始时间序列文件的标签(tag)信息，直接读取采样率为2400hz的原始时间序列文件的每个信号道的数据，将每个信号道的数据保存在二维数组存储区中对应的行数中。
17.进一步地，所述二维数组存储区的分配方法如下：二维数组存储区＝zeros(采样道数，采样率
×
采集总时间)，其中，zeros指的是二维数组存储区初始化全部全部置0，所述采样道数是指二维数组存储区的行数，所述采样率
×
采集总时间为二维数组存储区的列数，即每行数据的总长度，所述采样率、采样道数和采集总时间通过读取采样率为2400hz的原始时间序列文件的标签(tag)信息得到。
18.进一步地，所述标签(tag)信息为采样率为2400hz的原始时间序列文件的每个信号道的二进制数据前32个字节。
19.进一步地，所述一维数组的数据长度与所述二维数组的列数相同。
20.进一步地，步骤(4)具体包括如下步骤：将滤波后数据转换成v8数据处理软件ssmt2000能识别，且以ts3为后缀的文件。
21.进一步地，步骤(4)具体包括如下步骤：将二维数组中被滤波的信号道数据替换成一维数组数据，然后将替换后的二维数组数据转换成二进制数据，文件保存后缀为ts3，文件名与采样率为150hz的原始时间序列文件相同。
22.本发明的有益技术效果：
23.本发明提供的音频大地电磁数据工频噪声干扰的滤波方法，能够有效既可以滤除工频噪声干扰，获取了高质量的原始数据，提高了工作效率，减少了野外工作量，可以更好地为地质勘探工作服务。
附图说明
24.图1为采样率为150hz的受50hz工频噪声干扰的原始时间序列文件ts4的信号电道ex原始时间序列信号；
25.图2为采样率为150hz的受50hz工频噪声干扰的原始时间序列文件ts4的信号电道ey原始时间序列信号；
26.图3为采样率为150hz的受50hz工频噪声干扰的原始时间序列文件ts4的信号电道ex功率谱曲线；
27.图4为采样率为150hz的受50hz工频噪声干扰的原始时间序列文件ts4的信号电道ey功率谱曲线；
28.图5为采样率为2400hz的受50hz工频噪声干扰的原始时间序列文件ts3的信号电道功率谱曲线；
29.图6为采样率为2400hz的受50hz工频噪声干扰的原始时间序列文件ts3经去噪后的信号电道功率谱曲线；
30.图7为采样率为2400hz的受50hz工频噪声干扰的原始时间序列文件ts3的视电阻率曲线图；
31.图8为采样率为2400hz的受50hz工频噪声干扰的原始时间序列文件ts3的相位曲线图；
32.图9为采样率为2400hz的受50hz工频噪声干扰的原始时间序列文件ts3经去噪后的视电阻率曲线图；
33.图10为采样率为2400hz的受50hz工频噪声干扰的原始时间序列文件ts3经去噪后的相位曲线图。
具体实施方式
34.下面结合实施例对本发明作进一步详细描述。
35.实施例1
36.v8多功能电法仪采集主机包括两个信号电道(ex、ey)和三个信号磁道(hx、hy、hz)，v8多功能电法仪采集盒子仅包括两个信号电道(ex、ey)信号。
37.本实施例针对v8多功能电法仪采集的音频大地电磁原始时间序列信号受50hz工频噪声干扰的问题，提供一种音频大地电磁数据工频噪声干扰的滤波方法，该滤波方法采用作带阻滤波，滤除信号电道(ex或ey)或磁道(hx或hy)的工频噪声，包括如下步骤：
38.步骤一，判断音频大地电磁数据的原始时间序列文件的信号道功率谱曲线是否存在噪声干扰，得到第一判断结果，具体包括如下步骤：
39.tsviewer软件打开v8采集主机的采样率150hz的4245102a.ts4文件，读取4245102a.ts4文件中的数据，对数据的每个信号道均做功率谱分析，得到每个信号道的功率谱曲线，根据每个信号道的功率谱曲线判断所述信号道是否存在工频噪声干扰，
40.参见图3
‑
4，信号道ex和ey的功率谱曲线在45hz～55hz之间出现明显的毛刺噪声干扰，在50hz处增益为1db，明显大于
‑
3db，则判断信号道ex和ey存在工频噪声干扰，需要做滤波处理。
41.步骤二，在matlab软件中，读取4245102a.ts3文件的标签(tag)信息，标签(tag)信息包括采样率2400hz、采样道数4和采集总时间1285秒；根据采样率、采样道数和采集总时间分配二维数组存储区，分配方法如下：二维数组存储区＝zeros(采样道数，采样率
×
采集总时间)，其中zeros指的是二维数组存储区初始化全部全部置0，采样道数是指二维数组存储区的行数，采样率
×
采集总时间为二维数组存储区的列数，即每行数据的总长度，二维数组存储区的行数为4，分别为信号道ex、信号道ey、信号道hx和信号道hy，列数为1285
×
2400。
42.步骤三，在matlab软件中，跳过4245102a.ts3文件的标签(tag)信息，直接读取每
个信号道的数据，将每个信号道的数据保存在步骤二的二维数组存储区中对应的行数中,ex信号道数据保存在第一行，ey信号道数据保存在第二行，hx信号道数据保存在第三行，hy信号道数据保存在第四行。
43.步骤四，设定凯瑟窗参数，调用matlab中的凯瑟窗函数构建凯瑟窗winfun，winfun＝(kaiser(n,β))'；将二维数组第一行数据(信号道ex)与winfun作卷积运算，完成对ex的滤波。
44.凯瑟窗函数公式如下：
[0045][0046]
i0是修正后的零阶贝赛尔(bessel)函数，β为依赖于滤波器长度n的参数，给定频率带阻滤波器阻带低频边缘ω(p)、带阻滤波器阻带高频边缘ω(s)和衰减系数as，则β和n的取值公式如下：
[0047][0048][0049]
δω＝ω(p)
‑
ω(s)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(4)
[0050]
ω(p)设置为45hz之间，ω(s)设置为55hz之间；as设置为80。利用公式(2)和公式(3)分别计算β和滤波器长度n，当n计算为偶数时，将其加1转换成奇数。
[0051]
步骤六，将滤波后的信号道ex数据保存在长度为1285
×
2400的一维数组中。
[0052]
参见图5
‑
6，滤波后的功率谱曲线在45hz到55hz之间的毛刺信号消失，在50hz处增益明显减小到
‑
4.8db。
[0053]
采用与信号道ex相同的参数设置对信号道ey滤波，滤波后的数据保存在另一个长度为1285
×
2400的一维数组中。
[0054]
步骤七，将二维数据存储区中的ex和ey道数据分别替换成滤波后一维数组中的数据，将被替换后二维数组数据转换成ssmt2000软件能识别，以ts3为后缀，文件名为4245102a的文件。
[0055]
对滤波后的数据采用ssmt2000软件处理，滤波前与滤波后处理结果见图7
‑
10。图7
‑
8中，滤波前在50hz处视电阻率和相位均有明显的跳点，说明是噪声干扰引起的。图9
‑
10中，滤波后在50hz处视电阻率和相位跳点消失，曲线变得相对连续，说明滤波后噪声得以消除。
[0056]
以上所述实施例仅表达了本发明的几种实施方式，其描述较为具体和详细，但并不能因此而理解为对本发明专利范围的限制。应当指出的是，对于本领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干变形和改进，这些都属于本发明的保护范围。因此，本发明专利的保护范围应以所附权利要求为准。