浅谈全相位时移相位差法在O...auser型磁力仪中的应用_任振华.pdf

|117科技论坛0 引言Overhauser 磁力仪是一种基于 Overhauser 效应原理的弱磁场测量仪器1，因其灵敏度高、功耗低等特点，被广泛应用于军事和各类工程领域。核的欧沃豪斯（Nuclear Overhauser Effect）效应是Overhauser 于 1953 年在电子自旋和核自旋的样品中首先发现的。当把特殊的含有不成对电子的液体与氢原子结合起来、并在射频（RF）磁场的极化下便产生 Overhauser 效应，即不成对的电子将其自身的极化强度传输给氢原子后，就产生了很强的进动信号，在地磁场的作用下，可在感应线圈上产生自由感应衰减（free induction decay,FID）信号，通过测量此信号的频率，便可得到地磁场强度，此种磁力仪称为 Overhause 磁力仪2。由于 FID 信号衰减过快，持续时间较短，通常在几百毫秒左右，为测量其频率，通过傅里叶变换得到其幅度谱，即可将要寻找的频率值锁定在一个范围内，但受限于信号时间长度，这个范围难以缩小。为提升测量的分辨率，可通过一些 FFT 插值算法去实现优化，本文即将全相位时移相位差法引入 FID 信号频率测量中。1 全相位时移相位差法原理全相位快速傅里叶变换（All-phase FFT，ApFFT）算法是一种新型谱分析方法。相比于传统的 FFT，ApFFT 具有优良的抑制频谱泄漏和抗噪声干扰的性能3，且各条谱线的相位值与频率偏离值无关，即“相位不变性”4。数据预处理和 FFT 计算构成 ApFFT 的两部分，前者相对复杂，图 1 示意了 ApFFT 这一流程，具体步骤如下5：(1)构造长为（2N-1）的卷积窗 Wc，对中心采样点x(0)前后的(2N-1)个数据进行加权；(2)原序列间隔 N 提取元素，构造长度为 N 的预处理新序列；(3)新序列进行 FFT 运算。Z-1Z-1Z-1Z-1+x(-2)x(-1)x(0)x(1)x(2)Wc(-2)Wc(-1)Wc(0)Wc(1)Wc(2)FFTy(-2)y(-1)y(0)Y(2)Y(1)Y(0)图 1 N 为 3 的 APFFT 谱分析框图可 认 为 对 于 长 为 2N-1 时 间 序 列 x(n)，其 中，-(N-1)n(N-1)，中的一点x(0)，存在且只存在N个包含该点的 N 维向量：x0=x(0),x(1),x(N-1)T x1=x(-1),x(0),x(N-2)T 式(1)xn-1=x(-N+1),x(-N+2),x(0)T通过循环移位，将(1)式中各向量样本点 x(0)移至各向量的首位，得到新的 NN 维向量：x0=x(0),x(1),x(N-1)T x1=x(0),x(1),x(-1)T 式(2)xn-1=x(0),x(-N+1),x(-1)T式(2)各向量相加求均值，即可得到全相位数据向量：xAP=1NNx(0),(N-1)x(1)+x(-N+1),.,x(N-1)+(N-1)x(-1)T 式(3)由 DFT 的移位性质可知，式(2)与式(1)的离散傅里叶变换()()2ikijNiXkXk e=和()()2ikijNiXkXk e=存在如下关系：浅谈全相位时移相位差法在 Overhauser型磁力仪中的应用任振华1，陈乾1，李斌2（1.西安华舜测量设备有限责任公司，陕西西安，710065；2.西北工业大学航海学院，陕西西安，710072）摘要：Overhauser型磁力仪要测量拉莫尔信号的频率，受限于信号采集时间，普通的FFT算法难以应用，而全相位快速傅里叶变换（ApFFT）算法有效地拓展了信号频谱分析方法的应用范围，本文将此方法引入Overhauser型磁力仪频率测量中，实现了高精度的测量，同时又没有增加硬件成本。关键词：拉莫尔频率；全相位FFT；频率测量；OverhauserDOI:10.16589/11-3571/tn.2023.06.004118|电子制作 2023 年 3 月科技论坛 ()()2ikijNiXkXk e=,0,1i kN 式(4)()()()()000112a0011222002()2()112002()2()22()2()22111111111sin ()1sinNNkijNpiiiiNNknkijjNNink ik nNNjjjNNinjkjkjjkjkNNjXkXkXk eNNx ni eeNeeeNeeeNeekeN =2()kN 式(5)由(5)式可知，0为 ApFFT 谱的相位值，即 x(0)的理论相位值，并且不受k的影响。ApFFT 谱分析具有优良的抑制频谱泄露和“相位不变性”的性质5，但频率还需要进一步校正。假设一个单频复指数信号如(6)式所示。00()(2/)()11jnjn Nx nAeAeNnN+=+式(6)式中0为信号数字频率，可变换为频率2/N的倍，则 x(n)的归一化（除以 N）FFT 谱计算如(7)式所示。00000212012()2()/2()/0()()()()/()/()/1()()11sin()sin()/nNknjjjNNiNjkjjjk n NjkNijkjkjkjjkNjkNjkNNjkNAX keeeNAAeeeeNNeAeeeeNeeekAeNkN =+=式(7)其中，0,1kN。对比(5)式和(7)式，可知：在 k=k0 点，FFT 计算给出振幅谱 a1=|X(k0)|，相位谱 p1=ang|X(k0)|，ApFFT 计算给出振幅谱 a1=|Xap(k0)|，相位谱 p1=ang|Xap(k0)|，则有：频率校正值：0121()ppkNN=+式(8)振幅校正值：2(1)2aAa=式(9)此方法利用不变的相位信息进行频率和振幅校正，可取得较高的精度。2 磁力仪的硬件组成地磁场是各种磁场成分的加，全球范围内地磁场值大小范围为 20000nT120000nT6，其对应的拉莫尔频率范围约为 0.85kHz。图 2 为 Overhauser 磁力仪整体结构框图，以主控CPU 为核心，完成信号采集、计算和上传。上位机主控CPU探头激励电路滤波放大电路AD采样电路外围接口电路切换电路电源管理 图 2 测量装置连接示意图系统采用 12V 锂电池供电，通过 DC/DC 进行电压转换，射频激励采用固定频率的有源晶振和功率放大器实现，通过开关电路控制直流激励，如图 3 所示。供电电池电源管理固定频率有源晶振功率放大电路射频输出直流激励脉冲输出开关电路直流电源图 3 激励电源管理探头包括激励和感应两组线圈，激励线圈使探头溶液中质子自旋发生极化，感应线圈接收进动信号产生的感应电压。激励电路包括射频激励和直流激励两部分，射频激励频率几十 MHz 左右，与拉莫尔信号频率相差较大，有利于信号采集滤波电路的设计。滤波放大电路包括探头配谐和信号放大两部分，配谐部分采用扫描办法锁定信号大致的频率范围，经过后续多级放大，使信号达到几百毫伏的级别，如图 4 所示。图 4 测量装置连接示意图AD 采样电路采用 TI 的芯片，型号为 ADS131A02，分辨率 24 位，最高采样速率 128kSPS，满足采样定理的要求。主控 CPU 采用 STM32L47x 系列芯片，其配有 FPU（浮点运算单元），给编程提供了很大的方便，同时具有多级休眠功能，为设备适应多种应用场景提供支持。|119科技论坛3 磁力仪的软件实现软件采用 C 语言编写，框架基于有限状态机的模式，工作流程如图 5 所示。开始外设初始化（RTC实时时钟、FatFS文件系统、USB、USART等）更新传感器工作参数，设定工作定时器进入系统状态机返回 图 5 软件流程图系统状态机包括五个状态，分别为：命令处理、地磁测量、数据处理、休眠、文件读取，如图 6 所示。命令处理状态，解码上位机命令，并根据内容更新相关工作参数，例如测量周期、直流激励脉冲时长等。地磁测量状态，包括探头极化和信号测量两个部分，软件负责控制各步骤的时序。探头极化通过定时器定时中断功能实现，确保时序控制的精确性，探头极化牵涉到射频激励与直流激励，前者可持续工作，后者有一定的时序要求，如图 7 所示。信号测量包括信号的窄带滤波和宽带放大。首先根据目标拉莫尔信号的频率范围，调整匹配网络的中心频率，通过窄带滤波的信号提升信噪比后，进入三级放大电路，幅值提升至 ADC 可检测的范围。激励线圈工作时长测量周期直流激励切换开关防抖延时接收线圈工作时长射频激励 图 7 探头工作时序数据处理状态，将 AD 转换得到的数据，应用全相位时移相位差法，估计出其频率值，图 8 示意其流程，然后根据拉莫尔信号的频率值与地磁场的关系，计算出地磁场强度，并基于频谱主峰幅值，给出信号质量评价，最终结果保存为文件，若需要则还可以通过串口传输至上位机。休眠状态，若工作周期时间较长，则通过 RTC 闹钟周期唤醒系统进行一次工作，然后再进入休眠；若暂时不需要设备工作，也可发送命令强制系统休眠，再次发送任何命令都会将唤醒系统。文件读取状态，系统检测到USB插入时，即进入此状态，运行文件读取、删除等任务，为实现更快的操作速度，系统在这一状态中仅执行USB相关的任务，不响应上位机的命令；当系统检测到 USB 拔出时，则退出此状态，返回上一状态。开始数据序列进行FFT运算，找出谱峰频率、幅值卷积窗与数据序列构造ApFFT预处理新序列ApFFT预处理新序列进行FFT运算根据频率校正公式估算频率结束 图 8 频率估算流程4 实验结果与总结为验证本文提出的测频方法的效果，用信号源 DG1032Z 生成按固定斜率衰减的正弦波信号，在 800Hz5kHz 范围内，分别对 0.8kHz、3kHz 和 5kHz 频率的信号进行测量，每 1 秒测量1个结果，连续进行100次，得到第一组数据；然后分别在原频率点的基础上，将信号频率增地磁测量数据处理文件读取命令处理测量完成默认进入处理完毕收到命令否USB插入USB拔出USB插入USB拔出USB插入USB拔出休眠处理完成闹钟唤醒外部I/O唤醒休眠命令休眠模式？是图 6 系统状态机120|电子制作 2023 年 3 月科技论坛加 0.001Hz，重复进行测量，得到第二组数据；最后将测量的时间结果转换为频率，实测数据如图 9 图 11 所示。0102030405060708090100799.977799.9772799.9774799.9776799.9778799.978测量点数/个频率/Hz信号频率800Hz及800.001Hz测量值分布图 信号频率800Hz信号频率800.001Hz 图 9 信号频率为 800Hz 及 800.001Hz 的测量结果01020304050607080901002999.9022999.9032999.9042999.9052999.9062999.9072999.9082999.9092999.912999.9112999.912测量点数/个频率/Hz信号频率3000Hz及3000.001Hz 测量值分布图 信号频率3000Hz信号频率3000.001Hz 图 10 信号频率为 3000Hz 及 3000.001Hz 的测量结果01020304050607080901004999.8354999.8364999.8374999.8384999.8394999.844999.8414999.8424999.8434999.8444999.845测量点数/个频率/Hz信号频率5000Hz及5000.001Hz测量值分布图 信号频率5000Hz信号频率5000.001Hz 图 11 信号频率为 5000Hz 及 5000.001Hz 的测量结果观察上述结果，虽然波形有很多“毛刺”，但基本可分辨出某一频率点 f 和 f+0.001Hz 的大小关系，总体上波形跳变范围随着信号频率增加而呈现出变大的趋势。测试数据的统计分析结果如表 1 所示。分析表 1 规律，测量误差均为负值，可加入补偿值提升测量准确性；随着待测频率变大，测量结果的标准差也增加，与测量装置的分辨率固定有关。测量结果的偏差和标准偏差上限分别为 0