一种广域电磁法工频滤波与实时计算技术

本发明涉及地球物理数据处理领域，更具体地涉及广域电磁法接收机采集数据的工频滤波和实时fft计算。

背景技术：

1、传统电磁法勘探主要缺陷有：理论方法上采用平面波近似公式，探测深部信息不准确；探测技术上单次探测数据量少，大量谐波数据被丢弃，勘测效率不高、分辨率较差等；仪器装备上采集信号信噪比不高，转置笨重成本高等。而广域电磁法建立了全新的曲面波电磁勘探理论，以严格的电磁波方程为基础，不做简化，将电流-频率-位置-大地电阻率-磁导率等全部信息考虑在内，只用测量电分量或磁分量中的一个，采用迭代方式计算出广域视电阻率，实现全息电磁法勘探；探测技术上广域电磁法与伪随机电法相结合，实现了一次可发送和接收多个主频信号，能够反映不同深度的地电信息，提高了勘测效率；目前何继善院士团队基于广域电磁法理论，已经设计研发出高精度广域电磁法接收机，测量精度在1％以内，模数转换分辨率达到了24bit。

2、目前广域电磁法接收机主控芯片采用的是xilinx的zynq系列芯片，主要包含arm和fpga两部分，资源丰富，功能强大。广域电磁法接收基于zynq芯片设计，实现了对广域电磁法电场信号的高精度采集，其主要思路为：传感器输入模拟信号后经过前置放大滤波和噪声压制后，经过a/d高精度模数转换输入到zynq芯片的fpga端，在fpga内实现对信号采集的控制以及输入信号的串并转换和缓冲，通过axi总线传输至arm的ddr存储器中储存。数据的处理一般将数据读入到pc电脑中进行后续处理。

3、fpga仅仅完成了广域电磁法电场信号的接收控制与传输，不能完全发挥fpga内部更多逻辑资源的作用。fpga作为可编程逻辑器件，基于ram实现，非常适合并行计算。fpga在数字信号处理领域也有着广泛的应用。相比于电脑上的软件处理，基于fpga设计广域电磁法电场信号的工频滤波和fft计算，延迟能够控制在us级，达到实时计算的效果。

技术实现思路

1、本发明提出一种广域电磁法工频滤波与实时计算技术，解决了广域电磁法电场信号采集数据中工频噪声影响大和不能实时fft计算的问题。

2、为实现上述目的，本发明通过以下技术方案予以实现：本发明提出一种广域电磁法工频滤波与实时计算技术，包括：基于现场可编程门阵列fpga芯片实现，由fir滤波器、巴特沃斯带阻滤波器、插值抽取模块、fft计算模块组成，包括：

3、一套fir高通滤波器参数和多套fir低通抗混叠滤波器参数，根据带计算广域电磁法电场信号频率特征，并基于凯瑟窗函数设计，由matlab的filter designer工具计算得出滤波器系数，计算出的fir滤波器系数数符合对称性，减小设计复杂度。

4、一套巴特沃斯50hz工频带阻滤波器参数，在50hz附近设置滤波器的通带低频截止频率、通带高频截止频率、阻带低频截止频率、阻带高频截止频率、通带最大衰减、阻带最小衰减等参数，在matlab上由buttord算法计算出高阶巴特沃斯带阻滤波器的阶数和截止频率，接着由butter函数求出巴特沃斯带阻滤波器系数，并转换为二阶级联形式，更容易在fpga上实现。

5、插值和抽取模块根据广域电磁法电场信号频率特征和a/d采样芯片采样率特征，设计插值与抽取倍数。为减少信号损失，先进行插值，fir低通滤波除去镜像后，再进行抽取，再通过fir低通滤波去除高频噪声混叠。经过插值抽取后的信号变换采样率同样满足奈奎斯特采样定理。

6、经过插值和抽取后的广域电磁法电场信号每次计算的fft点数统一变换为2048点，通过fpga上的fft ip核实时计算。

7、凯塞窗函数实现fir滤波器的方法如下：

8、理想低通滤波器可表示为：

9、

10、其中wc为滤波器低通截止频率，α为群延时。(1)式所对应的单位冲激响应为：

11、

12、由于上式冲激响应函数无限长，且非因果，需要使用窗函数进行截断。凯塞窗函数具有阻带衰减快，通带起伏小的优势，所以这里采用凯瑟窗函数设计fir滤波器，凯塞窗函数为：

13、

14、其中

15、

16、其中i0(x)表示零阶第一类贝塞尔函数。同样在matlab上的filter desiner工具上设置好滤波器参数后，便能直接求出基于凯塞窗的fir滤波器系数。

17、fir滤波器在fpga上实现，需要先求解出其滤波器系数，而matlab的filterdesigner工具上，选择基于凯塞窗的fir低通滤波器，设置好滤波器采样频率38400hz，可以根据截止频率8500hz和9500hz，通带最大衰减3db，阻带最小衰减50db等参数设计fir滤波器并自动计算出滤波器系数，对滤波器系数进行指定16位宽的量化，最后生成coe文件导入到fir ip核中，在fpga上实现fir滤波。

18、设计的工频滤波器滤除的频率在50hz左右，现将受工频影响较大的第2、3频组采样率降低为400hz，再设计7个参数：滤波器采样频率fs为400hz,通带低频截止信号频率fpl为48hz、通带高频截止信号频率fph为52hz、阻带低频截止信号频率fsl为49hz、阻带高频截止信号频率fsh为51hz、通带最大衰减ap为2db、阻带最小衰减as为16db。将截止信号频率转换为归一化角频率：

19、

20、将以上参数代入matlab中的buttord函数求出巴特沃斯带阻滤波器的阶数n为8和截止频率wn为[0.2421,0.2580]。利用n和wn求出巴特沃率零点z、极点p、增益g模式的传递函数h(z)：

21、

22、将上式分解为二阶级联形式：

23、

24、求出的系数a和b组成一个4行6列矩阵，增益g为0.9370。8阶巴特沃斯带阻滤波器可以拆分为4个二阶巴特沃斯带阻滤波器级联形式，其中每个二阶传递函数可看成是输入信号移位相乘的结果除去递归输出移位相乘的结果，存储系数g，a，b，在fpga实现二阶巴特沃斯带阻滤波器，通过级联多个二阶巴特沃斯带阻滤波器，实现高阶滤波器工频滤波。

25、插值与抽取模块，广域电磁法电场信号频点范围在0.01hz到8192hz频带上分布，共分为0频组、1频组、2频组、3频组、4频组、5频组共六个频组，每个频组包含7个频点，其中0频组、2频组、4频组主频成分按2nhz(-6≤n≤13)分布，1频组、3频组、5频组主频成分按0.75*2nhz(-6≤n≤13)分布。在信号接收过程中，0频组、1频组按38400hz频率采样接收，其他频组按4800hz频率采样接收。不同采样率下，各个主频进行fft计算最少采样点数也不同，而在进行fft计算时，则需要采样点数统一，以提高适用性与计算效率。这里采用先插值i倍提高采样率，采用fir低通滤波器去除多余镜像后然后抽取d倍降低采样率，将采样点数统一变换为2048点。利用fpga上的fft ip核计算2048点采样数据，在100mhz系统频率下，计算一次耗时在100us以内，满足实时计算的要求。

26、综上所述，与现有技术相比，本发明具有以下有益效果：

27、实现了广域电磁法接收机数字信号处理中的50hz工频干扰滤除，提高了原始数据信噪比；

28、基于fpga实现了实时fft计算，发挥了fpga并行处理信号的优势，相比软件处理速度更快，使广域电磁法接收机具有了同时“采集-计算-存储”的功能，减少了后期数据处理的复杂性，及时监控采集到的信号质量与频谱信息。