广域电磁法(WFEM)是中南大学何继善院士发明的一种人工源频率域电磁测深方法，广域电磁法及其观测系统具有勘探深度大、观测范围广、工作效率高、适用性强等优势，目前已经在油气藏探测、页岩气勘探、金属矿探测、煤田采空区探测等方面取得了丰硕成果，已探明矿产资源潜在经济价值超过1100亿元(何继善，2010；Ren et al., 2013a, b；高金定，2015；张乔勋等，2017；Yang et al., 2017；凌帆等，2017；任政勇等，2017；Li et al., 2018).

FPGA(现场可编程门阵列)具有规模大、运算速度快、功耗低、可重配置等特点，用FPGA芯片来实现实时数字信号处理对提高电法仪器的运行速度、缩小仪器体积、降低功耗具有重要的实际意义(汤井田等，2015；Yan et al., 2015；柳建新等，2017；Auken et al., 2006).近年来，以FPGA芯片为代表的可编程逻辑器件开始逐步应用到电法仪器的研制中来(雷小青等，2010；李从庆等，2011；张启升等，2012；Ren et al., 2014a, b；任政勇等，2017).中南大学蒋奇云博士等人研制了以DSP(数字信号处理器)为核心的第一代广域电磁法接收机，在此基础上，近年来致力于以Microsemi公司SmartFusion2 SoC FPGA为核心的新一代广域电磁法接收机研制与应用，并且取得了丰硕的成果(蒋奇云，2010).由于FPGA等先进电子信息技术在地球物理仪器中应用的相对滞后，现有广域电磁法接收机中FPGA的应用仅仅停留在数据采集控制等层面，很少发挥其高速并行的优势来实现对采集到的电场信号进行FIR滤波、FFT变换等实时数字信号处理，导致接收机采集到的电场信号数据量庞大且大量冗余.

FFT(Fast Fourier Transform)快速傅里叶变换是频率域电磁法接收机数字信号处理关键技术之一，通过FFT快速傅里叶变换，可以提取某些关键频率点电场信号的振幅和相位信息，以便获取所需要的观测参数，从而来判断矿产资源类型.然而由于FFT算法运算量大、硬件实现相对复杂等特点，再加上广域电磁法发送机发送信号的特殊性，接收机中如何实现各个频组数据FFT实时变换是新一代高性能广域电磁法接收机研制的一个关键问题，本文提出的采取插值和抽取相结合的采样速率变换方法，将各频组FFT数据等效变换到了1024点，并在FPGA芯片上进行了仿真与测试，有效地解决了这个难题.

广域电磁法是相对于传统的可控音频大地电磁法(CSAMT)和磁偶源频率测深法(MELOS)方法提出来的.广域电磁法继承了CSAMT使用人工场源克服场源随机性的优点，也继承了MELOS方法非远区测量的优势，摒弃了CSAMT可控音频大地电磁法远区信号微弱的劣势，扩展了观测适用的范围，同时也摒弃了MELOS方法校正方法，保留了计算公式中的高次项.广域电磁测深法既不是沿用卡尼亚公式，也没有把非远区校正到远区，而是采用适合全域的公式计算视电阻率即广域视电阻率，这种方法大大拓展了人工源电磁法的观测范围，提高了观测速度、观测精度和野外测试的效率.

所谓“广域”是指突破“远区”的局限，在包括远区也包括非远区的广大地区进行测量，把电磁探测的观测范围扩大到包括非远区的广大区域.其算法是直接从电偶源全区精确公式(1)出发，即：

(1)

在(1)均匀大地表面上水平电偶极源的E_x的精确表达式的基础上，定义广域视电阻率为

(2)

式中，，.式(2)是广域电阻率的计算基础.在包括远区也包括部分非远区在内的广大区域进行测量，观测人工源电磁场的一个分量(而不是彼此正交的一组电、磁分量)，计算广域视电阻率值，这就是广域电磁法的基本原理.

现有广域电磁法发送机和接收机信号共分为0至5共6个频组，具体信号频率成份及采样频率如表 1所示.其中0频组、2频组和4频组包括了主频为2ⁿ Hz(-6≤n≤13，n为整数)共20个频率成份的信号.1频组、3频组和5频组包括了主频为0.75×2ⁿ Hz(-6≤n≤13，n为整数)共20个频率成份的信号.4频组和5频组0.015625 Hz和0.01171875 Hz两个低频信号发送了两倍强度的信息，0频组和1频组的采样频率为38400 Hz，2频组到5频组采样频率均为4800 Hz.

表 1 (Table 1)

表 1 广域电磁法接收机信号特点 Table 1 Characteristics of receiver of Wide-Field Electromagnetic Method (WFEM) 频组名 信号频率成份/Hz 采样频率fs/Hz 频率分辨率Δf/Hz 最少采样点数N 0 128，256, 512, 1024 2048, 4096, 8192 38400 128 300 1 96, 192, 384，768 1536, 3072, 6144 38400 96 400 2 1, 2, 4, 8 16, 32, 64 4800 1 4800 3 0.75, 1.5, 3 6, 12, 24, 48 4800 0.75 6400 4 0.015625，0.015625 0.03125, 0.0625 0.125, 0.25, 0.5 4800 0.015625 307200 5 0.01171875，0.01171875 0.0234375，0.046875 0.09375，0.1875，0.375 4800 0.01171875 409600

表 1 广域电磁法接收机信号特点 Table 1 Characteristics of receiver of Wide-Field Electromagnetic Method (WFEM)

广域电磁法接收机要求能够检测出0频组到5频组40个频率点信号的振幅和相位信息，如前所述，各频组采样频率已经固定，待检测信号频率也是确定的，根据频率分辨率的概念，各频组频率分辨率Δf至少应该达到该频组最低频率信号的频率大小，如表 1所示，同时据此计算出各频组最少采样点数N.为了不引起频谱泄露，实际采样点数应为表 1所计算出来的采样点数的整数倍，亦即采样时间应为各频组最低频率信号周期的整数倍.

进一步观察表 1可以看出，满足频率分辨率要求的各频组FFT变换的点数均不相同，而且均不是2^M(M为正整数)点数；而且低频组信号的点数至少要达到409600点，而即便是目前世界上最先进的FFT IP核也仅仅只支持2¹⁷=131072点FFT变换，况且即便是目前世界上最先进的FPGA芯片资源难以满足这么多点FFT变换的需求，如何用FPGA实现达到广域电磁法要求的FFT实时变换成了当前以FPGA为核心的新一代广域电磁法接收机研制中数字信号处理的关键技术问题.

有理数I/D采样率转换技术首先对输入序列x(n)按整数因子I倍插值，然后再对插值输出的序列按整数因子D倍数抽取，达到按有理数因子I/D的采样率转换.输入序列x(n)经过I倍插值及h_I(n)数字低通滤波器以后得到序列x_I(n)，再经过h_D(n)数字低通滤波器后进行D倍抽取，得到输出序列x_ID(n).h_I(n)数字低通滤波器将零值样点变成插值样点，h_D(n)数字低通滤波器主要用作防混叠失真，都是工作在If_s采样频率上，它们可以合并为一个数字低通滤波器h(n)，它逼近理想的低通数字滤波器的特性，表达式为

(3)

其中，，.

本文利用插值和抽取相结合的有理数I/D采样率转换技术，将各频组FFT变换点数统一变换到1024点，提出了满足广域电磁法接收机要求的等长度FFT变换技术方案，具体如表 2所示.

表 2 (Table 2)

表 2 广域电磁法接收机信号等长度FFT变换实现方案 Table 2 Realization scheme of signal equal-length FFT transform for the WFEM receiver 频组名 原采样频率fs/Hz 原频率分辨率Δf/Hz 最少采样点数N 实际取采样点数N′ 实际频率分辨率Δf′/Hz 插值系数I 抽取系数D 等效采样频率f′s/Hz 实际FFT点数 0 38400 128 300 4800 8 16 75 8192 1024 1 38400 96 400 3200 12 8 25 12288 1024 2 4800 1 4800 4800 1 16 75 1024 1024 3 4800 0.75 6400 6400 0.75 4 25 768 1024 4 4800 0.015625 307200 307200 0.015625 1 300 16 1024 5 4800 0.01171875 409600 409600 0.01171875 1 400 12 1024

表 2 广域电磁法接收机信号等长度FFT变换实现方案 Table 2 Realization scheme of signal equal-length FFT transform for the WFEM receiver

从表 2可以看出，通过合理选取各频组实际采样点数N′和各频组插值I/抽取D系数，在各频组等效采样频率f'_s满足奈奎斯特采样定理和频率分辨率的前提下，各频组实际FFT变换的点数均变换成了1024点，实现了各频组的等长度FFT变换，降低了低频组对硬件的要求，同时还变换成了2^M(M为正整数)点数的形式，适合FPGA实现.

由于采集到的电场信号受到很严重的工频及其高次谐波干扰，对采集到信号直接做抽取操作，对各频率成份的振幅和相位均有非常大的影响，因此对广域电磁法接收机实际采集到的信号必须滤除这些干扰后才能进行插值或抽取操作.以广域电磁法接收机第3频组实际采集到的电场信号为例来测试等长度FFT变换技术在广域电磁法接收机数字信号处理中的可行性，信号处理步骤如下：

(1) 取第3频组第一低频(0.75 Hz)采样周期6400点原始数据.

(2) 原始数据送入100 Hz的低通数字滤波器，有效滤除高频干扰.

(3) 按照表 2提出的等长度FFT变换技术实现方案，对FIR滤波器滤波输出的6400点信号先进行4倍插值，然后再进行25倍抽取，将信号长度变换成1024点.

(4) 分别对100 Hz低通滤波输出的6400点信号、4倍插值后的25600点信号以及25倍抽取后的1024点信号进行FFT快速傅里叶变换.

(5) N从0开始，分别提取N=1、2、4、8、16、32、64点幅值(模值)和相位信息(分别对应着0.75 Hz、1.5 Hz、3 Hz、6 Hz、12 Hz、24 Hz、48 Hz频谱成份)，并计算各频率成份的振幅.

图 1是信号尺度变换时域波形图，图 2是0~50 Hz有用信号频段范围局部放大的频谱图.从图 2可以看出，6400点信号经过FIR低通滤波后进行的尺度变换对0~50 Hz有用信号影响非常小.

图 1

Fig. 1

图 1 信号尺度变换时域波形图 Fig. 1 Time domain waveforms of signal-scale transformation

图 2

Fig. 2

图 2 信号尺度变换频域波形图(局部放大) Fig. 2 Frequency domain waveforms of signal-scale transformation(local amplification)

为了进一步定量分析方案的可行性，本文计算了变换前后各有用频率成份的振幅和相位等信息，如表 3所示.可以看出，经过4倍插值和25倍抽取后，还是能较好地从变换后的信号中有效提取出各频率成份的振幅和相位信息，因此，本文所提出的等尺度FFT变换技术是完全可行的.在实际应用中，还可以通过进一步提高FIR滤波器的阶数等措施，进一步加大滤波器阻带范围内衰减，从而进一步减小由于高频干扰引起的误差.

表 3 (Table 3)

表 3 尺度变换过程中各频率成份信息(振幅单位nV) Table 3 Information of frequency components in the process of scale transformation 频率/Hz 100 HzFIR低通滤波后信号(6400点) 4倍插值以后信号(25600点) 25倍抽取以后信号(1024点) 幅值/nV 相位/rad 振幅/nV 幅值/nV 相位/rad 振幅/nV 幅值/nV 相位/rad 振幅/nV 0.75 282.9 -0.7913 88.4 1030.6 -0.7923 88.3 45.2 -0.7891 88.2 1.5 225.7 -0.2014 70.5 901.3 -0.2018 70.4 36.6 -0.1976 71.6 3 98.9 1.5801 30.9 395.8 1.5841 30.9 15.7 1.5501 30.7 6 69.7 -0.0785 21.8 277.0 -0.0790 21.6 12.5 -0.1280 24.5 12 77.3 1.7522 24.2 309.6 1.7572 24.2 12.5 1.7369 24.5 24 76.2 0.2213 23.8 303.4 0.2224 23.7 13.0 0.2587 25.3 48 65.2 -0.0821 20.4 259.1 -0.0826 20.2 11.1 -0.0590 21.7

表 3 尺度变换过程中各频率成份信息(振幅单位nV) Table 3 Information of frequency components in the process of scale transformation

利用Altera公司开发的FFT快速傅里叶变换的IP核，定制1024点的实时FFT处理模块，构建FFT信号处理模型，如图 3所示.

图 3

Fig. 3

图 3 1024点FFT快速傅里叶变换模型 Fig. 3 Model of 1024-point FFT

图 4是1024点广域电磁法接收机实测信号FFT变换幅频特性曲线，上图为matlab理论仿真结果，下图为图 3模型计算输出的1024点复数取模输出的结果，可以看出，经过对第3频组6400点原始数据进行4倍插值和25倍抽取变换成的1024点数据，等效的采样频率也变成了768 Hz.两个幅频特性曲线频谱结构几乎一致，幅度上前者是后者的2¹¹=4096倍，这是因为exp指数为-11，实际等效为output×2^-(-11)=output×2¹¹= output×4096.

图 4

Fig. 4

图 4 1024点广域电磁法接收机实测信号FFT幅频特性曲线 Fig. 4 FFT amplitude-frequency curves of 1024-point WFEM receiver

表 4是提取的0.75 Hz、1.5 Hz、3 Hz、6 Hz、12 Hz、24 Hz和48 Hz七个频率点信号的振幅和相位信息.左侧是3频组经过插值\抽取变换后的1024点FFT变换理论计算值，这和表 3是一致的，右侧是图 3 FFT模型计算结果，可以看出1024点信号理论计算值和FFT模型计算结果基本一致.

表 4 (Table 4)

表 4 1024点广域电磁法接收机实测信号各频率成份信息表 Table 4 Frequency component information of measured signal by 1024-point WFEM receiver 频率/Hz 1024点信号理论计算值 FFT模型计算结果 幅值/nV 相位/rad 振幅/nV 幅值/nV 相位/rad 振幅/nV 0.75 2959023.0 -0.7891 88.2 1446.8 -0.7893 88.3 1.5 2401315.3 -0.1976 71.6 1174.0 -0.1981 71.7 3 1028354.1 1.5501 30.7 503.1 1.5489 30.7 6 820597.5 -0.1280 24.5 401.3 -0.1275 24.5 12 822254.5 1.7369 24.5 402.6 1.7380 24.6 24 848830.2 0.2587 25.3 415.0 0.2608 25.3 48 726473.7 -0.0590 21.7 354.6 -0.0593 21.6

表 4 1024点广域电磁法接收机实测信号各频率成份信息表 Table 4 Frequency component information of measured signal by 1024-point WFEM receiver

将模型转换成RTL寄存器传输级硬件描述语言，用Modelsim软件进行了时序仿真，用Quartus Ⅱ软件中进行了编译及时序分析.结果表明：在Altera公司Cyclone Ⅱ系列EP2C35F484C8芯片上实现1024点FFT快速傅里叶变换最高响应频率达到了93.84 MHz，在Cyclone Ⅱ系列EP2C35F484C8芯片上完成1024点数据输入、FFT快速傅里叶转换以及结果输出最快只需要约34 μs，因此完全可以满足广域电磁法接收机数字信号FFT变换实时处理的需求.

在EP2C35F484C8型FPGA芯片上，完成1024点FFT快速傅里叶变换，内部逻辑单元只消耗了17%，存储单元只消耗了51%，嵌入式硬件乘法器只消耗了34%.因此，在一块FPGA芯片上，完全可以考虑根据需要有选择性的将FIR数字滤波器、IIR工频滤波器、LMS自适应滤波器、FFT快速傅里叶变换以及Nios Ⅱ嵌入式软核处理器等模块全部集成在同一块FPGA芯片上实现，实现SOPC可编程的片上系统广域电磁法接收机.

针对广域电磁法接收机数字信号FFT实时变换的需求，提出了通过插值和抽取相结合、适合广域电磁法接收机数字信号处理的等长度FFT变换实现技术，并以广域电磁法接收机实际电场信号进行了测试.结果表明：在对广域电磁法接收机实际信号进行FIR低通滤波后，能够通过插值和抽取相结合的方法，将各个频组FFT变换的长度统一变换成1024点，解决了各个频组FFT变换长度不统一、不适合FPGA实现的技术难题.最后利用Altera公司的FFT IP核，实现了1024点FFT变换并进行了测试，完全能够满足广域电磁法接收机FFT实时实现的需求.对于基于FPGA的广域电磁法接收机数据采集模块与信号处理模块的信息协同处理问题，值得进一步研究.