CN106643944A - 一种科里奥利质量流量计高精度频率跟踪及估计方法 - Google Patents

Abstract

本发明提供了一种科里奥利质量流量计高精度频率跟踪及估计方法，涉及数字信号处理领域，采用改进的分段相位差频率估计方法对科里奥利质量流量计输出信号进行解算，以得到与流体流量质量和密度相关的瞬时频率，本发明采用了数字信号处理方法估计CMF谐振频率，不受模拟电路时延和RLC参数影响，具有灵活、高速与高精度的优点，无需使用自适应线性增强信号等复杂算法，在处理数据的复杂度上相比于基于AFF和SGA的科氏质量流量计数字信号处理系统有明显改善，正确稳定跟踪到CMF输出频率精度达到99.999％。

Description

一种科里奥利质量流量计高精度频率跟踪及估计方法

技术领域

本发明涉及数字信号处理领域，尤其是科里奥利质量流量计的频率跟踪及估计方法。

背景技术

科里奥利质量流量计(Coriolis Mass Flowmeter，CMF)是一种基于谐振原理和科氏效应(流体在直线运动时产生与质量流量成正比的科里奥利力原理)，并能以极高的测量精度和可靠性实现精确测量流体流量质量的高端仪器仪表。目前，国内市场上销售的CMF大多采用模拟驱动方式和模拟信号处理方法，易受外界震动和流体流速不稳定的影响，使得工业现场测量质量流量精度低于实验时标定的精度。模拟电路驱动和处理方法是对来自CMF输出的两路频率相同、相位不同的正弦信号进行滤波和过零点检测，将其转换成两路方波信号，再利用方波信号的上升沿或下降沿来计数解算信号的频率。模拟驱动及信号处理的缺点在于驱动增益有限、频率跟踪缓慢，不能解决两相流、批料流发生时流量管固有频率及阻尼比剧烈变化情况下的流量管停振而无法跟踪到频率的问题。

《面向时变的科里奥利质量流量计信号的处理方法研究与实现》公开了一种科里奥利质量流量计的信号处理方法，提出将多抽一滤波器、自适应格型陷波滤波器和负频率修正的滑动DTFT递推算法组合起来作为一套科里奥利流量计信号处理方法，该方法采用两级多抽一滤波器，同时采用30阶FIR低通滤波器，虽然能较好的保证线性相位特性，但是其数据运算量大、处理实时性不强，只能作为一种计算机仿真验证方法来研究，很难运用在存储容量有限的嵌入式数字信号处理器上。针对CMF输出两路信号的高精度频率算法估计和跟踪表述，该方法对于其采用的自适应格型陷波滤波器的参数描述模糊，虽然提到使用Burg算法进行自适应调整，却没有详细说明参数整定的具体方法。

合肥工业大学徐科军等的中国专利“基于AFF和SGA的科氏质量流量计数字信号处理系统，CN1467485”采用自适应陷波器(Adaptive Funnel Shaped Filter，AFF)，以自适应线性增强信号、跟踪和测量信号频率，能兼顾频率跟踪能力及跟踪精度，但问题在于该方案是基于定点型DSP(Digital Signal Processor，DSP)开发的算法。参考德州仪器TMS320系列DSP数据手册以及文献《DSP基本体系结构和特点》可知，定点型DSP相比于浮点型DSP，定点DSP的字长每增加1bit，动态范围扩大6dB，在处理低信噪比信号的场合，需要对数据不断地移位定标或者截尾。但是频繁地移位定标要耗费大量的程序空间和执行时间，截尾处理则会使计算的精度劣化。32bit浮点运算DSP的动态范围可以达到1536dB，这不仅扩大了数据动态范围，还因为大大减少了定标，移位和溢出检查，从而提高了运算精度，可节省大量运算时间和存储空间。由于浮点DSP的浮点运算用硬件来实现，可以在单周期内完成，因而其处理速度远高于定点DSP，这一优点在实现高精度复杂算法时尤为突出，为复杂算法的实时处理提供了保证。

西安东风机电有限公司的中国专利“一种基于谐波分解的科里奥利质量流量计数字信号解算装置及方法，CN1837758”提到其在DSP中采用特殊设计的梳状滤波器实现带通滤波来抑制直流和工频干扰，且其所述的滤波器可以覆盖其使用的CMF装置的全部工作频率。该专利采用FIR滤波器来实现其所述功能，且要求中心频率是固定的，若采用该专利所述方法实现CMF高精度频率估计，所要面对的问题在于解决其所述的覆盖所有CMF装置的全部工作频率，此算法设计复杂，虽然可以起到稳定滤波器特性的效果，但需要以付出大量的运算空间和运算时间为代价，其使用的最小二乘法对信号进行谐波分解，虽然不受非整周期截断的影响，但这种方法计算量大，实时跟踪频率的效果不佳。

综上所述，对于高精度估计科里奥利质量流量计的输出信号频率，且保持实时跟踪，需要采用合适的频率估计跟踪算法，尽量优化算法流程，在充分发挥硬件性能的基础上，保证频率估计的精度和算法的快速收敛性，并且不增加算法的时间复杂度和空间复杂度。

发明内容

为了克服现有技术的不足，本发明为提高科里奥利质量流量计输出信号频率估计的精度，保证跟踪频率的稳定性，提出采用改进的分段相位差频率估计方法对科里奥利质量流量计输出信号进行解算，以得到与流体流量质量和密度相关的瞬时频率。

本发明解决其技术问题所采用的技术方案包括下述步骤：

步骤一、根据CMF一次仪表有效输出信号确定贝塞尔滤波器的通带截止频率、通带最大衰减、阻带截止频率、阻带最小衰减四个技术指标，将四个技术指标输入电路设计软件，得到用于科里奥利质量流量计二次仪表的滤波器电路图；

步骤二、按照步骤一生成的滤波器电路图，在CMF二次仪表的模-数转换器

(Analog-to-Digital Converter，AD)前端搭建有源模拟贝塞尔低通滤波器；

步骤三、用一个具有4通道同步采样同步转换功能的AD转换器将CMF一次仪表的输出信号转换为数字信号；

步骤四、对于步骤三中AD转换器转换后输出的数字信号，用离散正弦信号x[n]表示：

式(1)中A为信号幅度，ω₀为信号数字角频率，φ为信号初始相位，n＝0,1,2…N-1，取x[n]的前N/2个点离散傅里叶变换(Discrete Fourier Transform，DFT)，变换后的频域谱为X[k]，则

其中k＝0，1，2…N/2-1，f_s为信号采样频率，设k_max为离散信号快速傅里叶变换(Fast Fourier Transformation，FFT)处理后的峰值谱位置，式(2)中的表示为(k_max+δ)，则X[k]写为：

则该正弦信号的频率ω₀表示为如下形式：

提取步骤三AD转换器转换的数字信号的4通道中的一路，截取一组长度为2N-1的输出数组，以标志位x[i]分成前后两个子段，即{x(n)，-N+1≤n≤i}与{x′(n)，i≤n≤N-1}，得到该同一路序列前后两子段的频偏估计：

其中为前向FFT序列的峰值谱的相位谱，为后向FFT序列的峰值谱的相位谱；

步骤五、使用解卷绕补偿2π角度使回归线性相位，综合(4)、(5)式得到：

公式(6)中为信号的无偏估计频率。

本发明的有益效果是采用了数字信号处理方法估计CMF谐振频率，不受模拟电路时延和RLC参数影响，具有灵活、高速与高精度的优点，本发明涉及的CMF高精度频率跟踪方法无需使用自适应线性增强信号等复杂算法，在处理数据的复杂度上相比于基于自适应漏斗形滤波器(Adaptive Funnel Shaped Filter，AFF)和滑动Goertzel算法(SlidingGoertzel Algorithm，SGA)的科氏质量流量计数字信号处理系统有明显改善，正确稳定跟踪到CMF输出频率精度达到99.999％。

附图说明

图1是本发明涉及的科里奥利质量流量计高精度频率跟踪方法流程图。

图2是本发明的低通滤波器的幅频特性曲线图。

图3是本发明的低通滤波器的群延迟图。

具体实施方式

下面结合附图和实施例对本发明进一步说明。

如图1所示，本发明涉及的CMF高精度频率跟踪方法需要经过确定滤波器性能指标、设计滤波器、AD同步采集转换信号、解算CMF输出频率和补偿CMF输出频率精度等五大部分组成，其具体实施步骤如下：

步骤一、根据CMF一次仪表有效输出信号确定贝塞尔滤波器的通带截止频率、通带最大衰减、阻带截止频率、阻带最小衰减四个技术指标，将四个技术指标输入德州仪器FilterPro Desktop滤波器设计软件，设计一个贝塞尔滤波器电路，用于CMF原始输出信号滤波，得到用于科里奥利质量流量计二次仪表的滤波器电路图；

步骤二、按照步骤一生成的滤波器电路图，在CMF二次仪表的模-数转换器(Analog-to-Digital Converter，AD)前端搭建有源模拟贝塞尔低通滤波器，用于滤除叠加在CMF一次仪表输出信号上的高频噪声，同时利用贝塞尔滤波器优良的线性相位特性保证信号相位失真程度最低；

为衰减高频信号，在模-数转换器前端设计一个低通滤波器，考虑到高精度频率估计算法处理的CMF输出信号要求相位不失真，故所采用的滤波器必须保持最佳的线性相位特性，本发明采用一个有源模拟贝塞尔低通滤波器，可以在CMF输出信号到达ADC之前，消除叠加在信号上的高频噪声，且保证信号相位失真程度最低。

步骤三、用一个具有4通道同步采样同步转换功能的AD转换器将CMF一次仪表的输出信号转换为数字信号；

将CMF一次仪表输出信号送入DSP处理之前需要转换成数字信号，本发明采用的具有4通道同步采样同步转换功能的AD转换器，实现了对同步性的要求。

步骤四、对于AD输出的信号，本发明用离散正弦信号表示为：

取x[n]谱的前N/2点离散傅里叶变换为X[k]：

其中k＝0，1，2…N/2-1，f_s为信号采样频率，设k_max为FFT的峰值谱位置，式(2)中的可表示为(k_max+δ)，则X[k]为：

该正弦信号的频率ω₀为：

提取步骤三AD转换器转换的数字信号的4通道中的一路，截取一组长度为2N-1的输出数组，以标志位x[i]分成前后两个子段，即{x(n)，-N+1≤n≤i}与{x′(n)，i≤n≤N-1}，得到则该同一路序列前后两子段的频偏估计：

其中为前向FFT序列的峰值谱的相位谱，为后向FFT序列的峰值谱的相位谱。

步骤五、反正切计算得到的相位需要回归到线性相位特性，使用解卷绕补偿2π角度使得回归线性相位，综合(4)、(5)式得到：

公式(6)中为信号的无偏估计频率。

实施例：

步骤一、针对国内某厂家生产的四个通用系列型号的CMF流量管，其CMF信号输出频率的范围介于40-240Hz，CMF一次仪表输出信号频率除共振频率之外，还包含有各种高频、低频噪声信号；此外，CMF输出的信号共振频率还与流体的性质、流速和密度状态有关，当流量管处于两相流或者批料流状态时，流量管的固有频率及阻尼比都会发生剧烈变化，在处理CMF一次仪表输出的信号时需对信号进行预处理，本发明使用一个低通滤波器滤除CMF一次仪表输出信号的高频噪声。本发明使用的模拟低通滤波器技术指标确定为：通带截止频率Ω_p＝400Hz，通带最大衰减α_p＝3dB，阻带截止频率Ω_s＝600Hz，阻带最小衰减α_s＝-30dB。

步骤二、本发明采用两级多路反馈的3阶贝塞尔低通滤波器来实现CMF原始输出信号的原始滤波。本发明设计的低通滤波器的通带增益为1V，滤波器由RC有源电路(Resistence-Capacitance电路，以下简称RC电路)组成，RC电路前一部分为一个1阶低通贝塞尔滤波器，RC电路后一部分为一个2阶低通贝塞尔滤波器，两个贝塞尔低通滤波器级联成一个3阶贝塞尔低通滤波器。如图2所示，是本发明得到的三阶模拟贝塞尔低通滤波器的幅频特性曲线，在信号频率为255.44Hz处，信号已衰减-3.03dB，可见滤波器对240Hz以上的频率能较好的抑制，图3是滤波器的群延迟图，频率小于257.62Hz时，滤波器能保持较好的线性相位特性，群延迟时间为690μs。

步骤三、本发明基于CMF输出两路信号同步性的要求采用了一片具有4通道同步采样同步转换功能的AD转换器，实现对CMF一次仪表输出的A、B两路信号的同步采样。根据奈奎斯特采样定律，AD采样频率必须设置为阻带截止频率的2倍以上，故本发明使用采样频率为4687.5Hz，以此来满足抗混叠滤波以及FFT运算要求。

步骤四、离散时间情况下的正弦信号一般表达式为：

取x[n]谱的前N/2点离散时间信号的DFT观测为X[k]，A为CMF输出信号的幅度，工程上要求信号的幅度满足1Hz/1mV，为方便推理论证，本发明以幅度A表示，则X[k]可表示为：

本发明把该正弦信号的频率ω₀可表示为如下形式：

k_max为FFT的峰值谱位置，本发明通过提取一路转换后得到的长度为2N-1的输入数组，再截取一个随机标志位x[i]的前后两个子段，即{x(n)，-N+1≤n≤i}与{x′(n)，i≤n≤N-1}，该同一路序列前后两子段的频偏估计为：

其中为前向FFT序列的峰值谱的相位谱，为后向FFT序列的峰值谱的相位谱。

步骤五、本发明在求解到两子段的相位差之后，计算得到的相位 需要回归到线性相位特性，本发明使用解卷绕补偿2π角度使得回归线性相位，再综合(4)、(5)式，得到：

公式(6)中为信号的无偏估计频率。

Claims (1)

1.一种科里奥利质量流量计高精度频率跟踪及估计方法，其特征在于包括下述步骤：

步骤一、根据CMF一次仪表有效输出信号确定贝塞尔滤波器的通带截止频率、通带最大衰减、阻带截止频率、阻带最小衰减四个技术指标，将四个技术指标输入电路设计软件，得到用于科里奥利质量流量计二次仪表的滤波器电路图；

步骤二、按照步骤一生成的滤波器电路图，在CMF二次仪表的模-数转换器前端搭建有源模拟贝塞尔低通滤波器；

步骤三、用一个具有4通道同步采样同步转换功能的AD转换器将CMF一次仪表的输出信号转换为数字信号；

步骤四、对于步骤三中AD转换器转换后输出的数字信号，用离散正弦信号x[n]表示：

$x\[n\]=A\cos ({\mathrm{\ω}}_{0}n+\mathrm{\φ})=\frac{A}{2}{e}^{j\mathrm{\φ}}{e}^{{\mathrm{j\ω}}_{0}n}+\frac{A}{2}{e}^{-j\mathrm{\φ}}{e}^{-{\mathrm{j\ω}}_{0}n}---\left(1\right)$

式(1)中A为信号幅度，ω₀为信号数字角频率，φ为信号初始相位，n＝0,1,2…N-1，取x[n]的前N/2个点离散傅里叶变换，变换后的频域谱为X[k]，则

$X\[k\]=\frac{A}{2}\·e^{j\mathrm{\φ}}\·e^{j\mathrm{\π}\frac{N-1}{N}(\frac{{\mathrm{\ω}}_0}{2\mathrm{\π}}\·\frac{N}{f_s}-k)}\frac{\sin \mathrm{\π}(k-\frac{{\mathrm{\ω}}_0}{2\mathrm{\π}}\·\frac{N}{f_s})}{sin\[n\·(k-\frac{{\mathrm{\ω}}_0}{2\mathrm{\π}}\·\frac{N}{f_s})/N\]}---\left(2\right)$

其中k＝0，1，2…N/2-1，f_s为信号采样频率，设k_max为离散信号快速傅里叶变换处理后的峰值谱位置，式(2)中的表示为(k_max+δ)，则X[k]写为：

$X\[k\]=\frac{A}{2}\·e^{j\mathrm{\φ}}\·e^{j\mathrm{\π}\frac{N-1}{N}\left(k_{\max }+\mathrm{\δ}-k\right)}\frac{sin\mathrm{\π}\left(k-k_{\max }-\mathrm{\δ}\right)}{sin\[n\·\left(k-k_{\max }-\mathrm{\δ}\right)/N\]}---\left(3\right)$

则该正弦信号的频率ω₀表示为如下形式：

${\mathrm{\ω}}_0=(k_{\max }+\mathrm{\δ})\frac{2{\mathrm{\πf}}_s}{N},\left|\mathrm{\δ}\right|\≤0.5---\left(4\right)$

提取步骤三AD转换器转换的数字信号的4通道中的一路，截取一组长度为2N-1的输出数组，以标志位x[i]分成前后两个子段，即{x(n)，-N+1≤n≤i}与{x′(n)，i≤n≤N-1}，得到该同一路序列前后两子段的频偏估计：

其中为前向FFT序列的峰值谱的相位谱，为后向FFT序列的峰值谱的相位谱；

步骤五、使用解卷绕补偿2π角度使回归线性相位，综合(4)、(5)式得到：

公式(6)中为信号的无偏估计频率。