CN110133738A - 基于IpDFT的质子磁力仪自由感应衰减信号的频率估计方法 - Google Patents

Abstract

本发明公开了基于IpDFT的质子磁力仪自由感应衰减信号的频率估计方法。质子磁力仪式量子弱磁测量仪器史上最古老的仪器，由于其操作简单、稳定性好等优点，已经成为半个世纪以来使用最为广泛的地磁测量仪器之一。质子磁力仪中自由感应衰减信号频率测量的精度决定了导出磁场的精度。因此，对于高精度频率估计的研究具有十分重要的意义。本发明提供的方法提高了算法的估计精度，能在很短的时间内估计出信号频率，在低信噪比的环境下，仍有较高的估计精度。

Description

基于IpDFT的质子磁力仪自由感应衰减信号的频率估计方法

技术领域

本发明涉及地磁测量领域，特别是涉及电网用倒供电防护装置基于IpDFT的质子磁 力仪自由感应衰减信号的频率估计方法。

背景技术

质子磁力仪是量子弱磁测量仪器史上最古老的仪器，由于其操作简单、稳定性好等 优点，已经成为半个世纪以来使用最为广泛的地磁测量仪器之一。传统的质子磁力仪的分辨率是0.1nT，精度是1.0nT。悬垂式磁力仪可以获得更高的精度，精度为0.1nT，分辨 率为0.01nT。虽然相比于质子磁力仪，悬垂式磁力仪精度更高，更灵敏，功率更低，但 是它也有一些缺点。悬垂式磁力仪探头中自由基溶液的极化总需要射频激励才能得到自 由感应衰减(FID)信号。在一些勘探项目中，如近距离梯度测量，仪器会相互干扰，导 致磁场测量误差；当悬垂式磁力仪与其他系统一起工作时，也会影响整个系统，很难完 全屏蔽射频。另外，由于悬垂式磁力仪价格昂贵，在多点弱磁监测应用中，成本较高。 另一方面，由于原理不同，质子磁力仪没有仪器间相互干扰的问题，而且比悬垂式磁力 仪便宜。由于这些优点，质子磁力仪被广泛应用。但是随着勘探技术的进步，质子磁力 仪需要具有更低的功耗、更高的精度和更高的灵敏度。

质子在地磁场周围的拉摩尔旋进可以再测量线圈中产生一个FID信号。由于FID信号频率与磁场强度成正比，频率测量的精度决定了导出磁场的精度。因此，对于高精度 频率估计的研究具有十分重要的意义。目前已经提出的频率估计方法有自适应算法、基 于自适应分数谱图的时频分布法、基于自相关因子的算法、基于离散傅里叶变换(DFT) 的估计算法。但这些算法存在一些缺点，有的算法计算复杂度比较大，无法在短时间内 估计出信号频率，有的算法抗噪性能不好，在信噪比比较小时估计精度较差。目前，对 于FID信号频率估计主要存在两个难点，由于FID信号衰减快，频率测量时间有限，如 何在短时间内达到较高的测量精度仍然是一个有待解决的问题；进一步提高在低信噪比 的环境下频率的测量精度是另一个有待解决的问题。

发明内容

为了解决以上问题，本发明提供基于IpDFT的质子磁力仪自由感应衰减信号的频率 估计方法，该方法将信号频谱中的正、负序列都考虑在内进行计算，提高了算法的估计精度，能在很短的时间内估计出信号频率，在低信噪比的环境下，仍有较高的估计精度， 为达此目的，本发明提供基于IpDFT的质子磁力仪自由感应衰减信号的频率估计方法， 包括以下步骤：

步骤1：采集N点离散时间FID信号s(n)；

步骤2：获取信号s(n)的N点离散傅里叶变换序列，记做S(k)，其中因子因子是其共轭，因子λ＝e^-βe^jω，λ^*＝e^-βe^-jω是其共轭，

步骤3：获取到频谱中k₁，k₂，k₃处三个值最大的频谱值S(k₁)、S(k₂)和S(k₃)， 其中k₁，k₂，k₃的值可以通过粗估计获得，k₂＝k₀，是l₀整数部分，可以通过寻找信号 频谱的谱峰所在位置，即得到k₂的值，k₁＝k₀-1，k₃＝k₀+1。

步骤4：根据插值傅里叶变换，利用S(k₁)、S(k₂)和S(k₃)三点频谱值得到信号频 率和阻尼因子的估计值。

作为本发明进一步改进，所述步骤4利用三点频谱值对信号频率和阻尼因子进行估 计的方法具体包括以下步骤：

步骤41：，因为有关系式(λ+λ^*)＝2e^-βcos(ω₀)，λλ^*＝e^-2β，所以解出λ的值就可 以得到信号阻尼因子β和频率ω₀的估计值，令(λ+λ^*)＝x，λ·λ^*＝y，通过对S(k₁)、 S(k₂)和S(k₃)表达式进行一系列的化简运算可以得到一元二次方程组

其中[]^*表示取共轭，公式(1)中各系数的值可以通过下面的公式求得

步骤42：求解公式(1)所示一元二次方程组，可以得到

其中代表x的估计值，代表y的估计值。

步骤43：根据关系式x＝2e^-βcos(ω₀)，y＝e^-2β，通过公式(3)最终可以得到信号频率ω₀和阻尼因子β的估计值，为

代表ω₀的估计值，代表β的估计值。

本申请工作原理如下：

在本发明中，利用插值傅里叶变换(IpDFT)来估计FID信号的频率。由于快速傅里叶变换(FFT)的存在，使得该算法计算量小，效率高，通过插值来解决离散傅里叶变换(DFT)中存在的栅栏效应问题。通过将信号的正、负频率都考虑在内进行计算，解决了 DFT频谱泄露的问题，提高了频率估计的精度，尤其在获取到的信号时间较短时。该算 法的稳定性、抗噪性能和计算复杂度要优于现有的同类频率估计算法。

本申请有益效果如下：

与现有的技术相比，本发明以下优点：1.利用信号的傅里叶变换序列进行插值运算， 降低了计算复杂度。2.将信号的正、负频率都考虑在内，提高了信号频率估计精度，能够 在较短时间内估计出信号频率。3.充分考虑了输入信号与输出信号的噪声，抗噪性能好。

附图说明

图1为无噪情况下，当信号频率变化时信号频率和阻尼因子估计值的均方误差图。其中 图(a)是频率估计值的均方误差随l₀的变化情况，图(b)为阻尼因子估计值的均方误差随l₀的变化情况。

图2为无噪情况下，当信号阻尼因子变化时信号频率和阻尼因子估计值的均方误差图。 其中图(a)是频率估计值的均方误差随β的变化情况，图(b)为阻尼因子估计值的均方误差随β的变化情况。

图3为有噪情况下，信噪比SNR＝40dB时，当信号频率变化时信号频率和阻尼因子估计 值的均方误差图。其中图(a)是频率估计值的均方误差随l₀的变化情况，图(b) 为阻尼因子估计值的均方误差随l₀的变化情况。

图4为有噪情况下，信噪比SNR＝40dB时，，当信号阻尼因子变化时信号频率和阻尼因子 估计值的均方误差图。其中图(a)是频率估计值的均方误差随β的变化情况，图(b)为阻尼因子估计值的均方误差随β的变化情况

图5为当信噪比SNR变化时信号频率和阻尼因子估计值的均方误差图。其中图(a)是频 率估计值的均方误差随SNR的变化情况，图(b)为阻尼因子估计值的均方误差随SNR的变化情况。

具体实施方式

下面结合附图与具体实施方式对本发明作进一步详细描述：

本发明提供基于IpDFT的质子磁力仪自由感应衰减信号的频率估计方法，该方法将 信号频谱中的正、负序列都考虑在内进行计算，提高了算法的估计精度，能在很短的时间内估计出信号频率，在低信噪比的环境下，仍有较高的估计精度。

作为本发明一种实施例，本发明提供基于IpDFT的质子磁力仪自由感应衰减信号的 频率估计方法，具体实施方式如下；

质子磁力仪产生的FID信号的离散时间表达式为

由上式可知，FID信号也就是带阻尼的实正弦信号，其中β表示信号的阻尼因子，A是信号幅度，ω₀＝2πl₀/N＝2π(k₀+δ₀)/N是信号频率，l₀表示获取的正弦信号的周期数， k₀是其整数部分，δ₀(||δ₀||≤0.5)是它的小数部分，N为采样点数，表示信号初始相位。 该信号N点DFT表达式为

其中S(k)表示第k点的频谱值，k∈{0,1,2,…,N-1}，通过化简合并得到因子因子是其共轭，因子λ＝e^-βe^jω， λ^*＝e^-βe^-jω是其共轭，通过计算可以得到关系式(λ+λ^*)＝2e^-βcos(ω₀)， λλ^*＝e^-2β。该发明利用两点处的频谱值来计算得到信号频率的估计值，通过公式(2)可以 得到频谱中值最大的三条谱线的表达式S(k₁)、S(k₂)和S(k₃)，令(λ+λ^*)＝x，λ·λ^*＝y， 有如下关系式

其中k₂＝k₀，是l₀整数部分，可以通过寻找信号频谱的谱峰所在位置，即得到k₀的值，k₁＝k₀-1，k₃＝k₀+1。公式(3)和(4)可以组成如下矩 阵

通过上述矩阵可以得到如下两个关系式

同理根据公式(4)和(5)可以得到以下关系式

观察分析上述四个公式，式(7)和式(9)消去v+v^*可得

同理式(8)和式(10)消去-(vλ^*+v^*λ)可以得到

上述两个公式中只有x、y两个未知数，因此可以组成一元二次方程组来求解x、y的值。

其中各系数的值分别为

通过计算得知a₁/a₂≈b₁/b₂≈c₁/c₂，即式(13)中两个方程近似为同一方程，因此不能 通过方程组(13)求得x和y的唯一解。因为这两个等式为复值，可以利用等号左右两边的实部、虚部相等得到一个新的一元二次方程组为

其中Re[]表示取实部，Im[]取虚部，上述方程组的解为

根据关系式x＝(λ+λ*)＝2e^-βcos(ω₀)，y＝λλ^*＝e-^2β可以得到阻尼因子β和信号频率ω₀的估计值为

至此，通过公式(17)即可得到FID信号阻尼因子和频率的估计值。

通过仿真实验来验证该方法的估计性能，并使用均方误差(MSE)这一性能指标来反映算法的估计情况。在仿真中设定信号幅度A＝1，信号初始相位在[0,2π)范围内随 机获取，采样点数及DFT点数N＝256。

首先分析无噪情况下算法的估计情况。图1展示了当β＝10^-4，频率l₀在0.5到4的范围内变化时，算法估计精度的变化情况。图(a)是频率估计值的均方误差随 l₀的变化情况，图(b)为阻尼因子估计值的均方误差随l₀的变化情况，从两图 中可以看出无论是对信号频率估计还是对阻尼因子的估计，该算法的估计精确，在-300dB 左右，此时估计误差主要来自于计算机软件本身的误差。且l₀变化时，算法估计精度基本 不变，估计性能稳定。如图2所示是在无噪声的情况下各算法的估计随β的变化情况， 其中图(a)是频率估计值的均方误差随β的变化情况，图(b)为阻尼因子估 计值的均方误差随β的变化情况。令l₀＝1.2，β在10^-4到10^-1范围内变化，从 两图中可以看出随着阻尼因子β的增大，该算法对频率和阻尼因子的估计精度都逐渐变 差，但变化不是特别明显。

估计算法的抗噪性能是很重要的性能之一，接下来分析在噪声背景下算法的估计情 况。当SNR＝40dB，β＝10^-4，l₀在0.5到4范围内变化时，算法的估计性能随l₀的变化 情况由图3所示。图(a)是频率估计值的均方误差随l₀的变化情况，图(b) 为阻尼因子估计值的均方误差随l₀的变化情况。从两图中可以看出即使在有噪 声存在的情况下，该算法对信号频率和阻尼因子仍能较为准确的进行估计，虽然当l₀＜1.5 时估计精度会稍有下降，这是因为此时正、负频率分布很近，会对估计算法造成一定的 影响，但这种情况下该算法的估计仍较为准确，可以满足估计要求。图4(a)和4(b) 分别展示了当SNR＝40dB，l₀＝1.2，β在10^-4到10^-1范围内变化时，各算法频率估计性能 和阻尼因子估计性能的变化情况。从图中可以看出，无论是还是都 随着β的增大而增大，即各个算法的估计精度都在变差。这是因为当阻尼较大时，采样 到的信号随着n的增大，信号能量减小速度很快，后面的采样点会几乎淹没在噪声中， 使得各算法估计不再准确。

各算法估计性能随SNR的变化情况如图5所示，仿真中令l₀＝1.2，β＝10^-4，SNR 在0dB到50dB的范围内变化。图(a)是频率估计值的均方误差随SNR的变化 情况，图(b)为阻尼因子估计值的均方误差随SNR的变化情况。如图所示， 该算法的估计性能随着SNR的增大而逐渐变好，说明该算法完全将负频率频谱泄露的影 响完全消除。即使在SNR很小的情况下，估计的均方误差仍较小，对信号频率和阻尼因 子仍有较为精确的估计，说明该算法有很好的抗噪性能。

以上所述，仅是本发明的较佳实施例而已，并非是对本发明作任何其他形式的限制， 而依据本发明的技术实质所作的任何修改或等同变化，仍属于本发明所要求保护的范围。

Claims (2)

1.基于IpDFT的质子磁力仪自由感应衰减信号的频率估计方法，其特征在于：包括以下步骤：

步骤1：采集N点离散时间FID信号s(n)；

步骤2：获取信号s(n)的N点离散傅里叶变换序列，记做S(k)，其中因子因子是其共轭，因子λ＝e^-βe^{j ω}，λ^*＝e^-βe^-jω是其共轭，

步骤3：获取到频谱中k₁，k₂，k₃处三个值最大的频谱值S(k₁)、S(k₂)和S(k₃)，其中k₁，k₂，k₃的值可以通过粗估计获得，k₂＝k₀，是l₀整数部分，可以通过寻找信号频谱的谱峰所在位置，即得到k₂的值，k₁＝k₀-1，k₃＝k₀+1。

步骤4：根据插值傅里叶变换，利用S(k₁)、S(k₂)和S(k₃)三点频谱值得到信号频率和阻尼因子的估计值。

2.根据权利要求1所述的基于IpDFT的质子磁力仪自由感应衰减信号的频率估计方法，其特征在于：所述步骤4利用三点点频谱值对信号频率和阻尼因子进行估计的方法具体包括以下步骤：

步骤41：，因为有关系式(λ+λ^*)＝2e^-βcos(ω₀)，λλ^*＝e^-2β，所以解出λ的值就可以得到信号阻尼因子β和频率ω₀的估计值，令(λ+λ^*)＝x，λ·λ^*＝y，通过对S(k₁)、S(k₂)和S(k₃)表达式进行一系列的化简运算可以得到一元二次方程组

其中[]^*表示取共轭，公式(1)中各系数的值可以通过下面的公式求得

步骤42：求解公式(1)所示一元二次方程组，可以得到

其中代表x的估计值，代表y的估计值。

步骤43：根据关系式x＝2e^-βcos(ω₀)，y＝e^-2β，通过公式(3)最终可以得到信号频率ω₀和阻尼因子β的估计值，为

代表ω₀的估计值，代表β的估计值。