基于正交解调的新型科氏流量计信号处理方法

A Novel Coriolis Flow Meter Signal Processing Method Based on Quadrature Demodulation

摘要

为实现科氏流量计实际应用中测量性能的最大化，提出一种基于相关运算和Hilbert变换的正交解调算法（CH-QD算法）.该算法利用自相关的同频余弦特性和Hilbert变换的相移特性产生正交解调所需的同频参考信号，有效解决了传统正交解调算法需预知信号频率、对频率检测算法的性能敏感等问题，提高了算法的检测精度，并实现了信号频率、相位差的联合估计.性能仿真测试表明，本文算法具有较高的检测精度、抗噪性和动态跟踪性能.在基于TMS320F28335DSP研制的科氏流量计变送器上实现本文算法，并进行性能分析和水流量标定实验.标定结果表明，在量程比20∶1的范围内，系统测量误差小于0.1%，重复性小于0.05%，具有较好的实用价值.

Abstract

In order to maximize the measurement performance of Coriolis flowmeter in practical application， a quadrature demodulation algorithm based on correlation operation and Hilbert transform （CH-QD algorithm） is proposed. The proposed algorithm uses the same-frequency cosine frequency characteristic of autocorrelation and the phase-shift characteristic of Hilbert transform to generate the reference signals required by quadrature demodulation， which effectively solves the problem that the traditional quadrature demodulation algorithm needs to predict the signal frequency and is sensitive to the performance of the frequency detection algorithm， which improves the detection accuracy of the algorithm， and realizes the joint estimation of signal frequency and phase difference. Simulation results show that the proposed algorithm has high detection accuracy， noise resistance and dynamic tracking performance. Based on TMS320F28335DSP， the algorithm is implemented on Coriolis flowmeter transmitter， and the performance analysis and water flow calibration experiments are carried out. The calibration results show that the measurement error of the system is less than 0.1% and the repeatability is less than 0.05% in the range of range ratio of 20∶1.

Graphical abstract

关键词

科里奥利质量流量计 / 信号处理方法 / 正交解调法 / Hilbert变换 / DSP

Key words

Coriolis mass flowmeter / signal processing methods / quadrature demodulation method / Hilbert transform / DSP

流体流量测量与人类的生产活动密切相关，精确的流量计量可以减少人力、物力和财力损耗.随着数字信号处理技术和仪器仪表的发展，科氏流量计（CMF）以其独特的测量优势广泛应用于石油化工、医疗制药、工业生产等诸多领域.

科氏流量计是一种利用流体流动产生的科里奥利力与质量流量成正比的原理而制成的直接式流量计，可以同时检测流体高精度的质量流量、密度等参数.当流体流经其以固有频率振动的测量管时，测量管两端由于科里奥利效应会产生不同的相位^［1］，传感器拾取该频率、相位差随流体密度和质量流量变化的振动信号，由变送器进行信号处理.由于质量流量与相位差成正比，密度与固有频率的平方成反比，因此，通过现代信号处理方法计算出两路传感器信号精确的频率和相位差，即可检测出流体的密度和质量流量.

多年来，众多学者投身于科氏流量计信号处理方法的研究中.文献［2］提出了计及负频率修正的离散时间傅里叶变换（DTFT）法，在计算离散傅里叶变换（DFT）时考虑负频率影响，提高了测量精度，但仍未从根本上解决频谱泄露问题.文献［3］采用四点拉格朗日插值的数字过零法，提高了零点检测的精度.文献［4］采用傅里叶级数拟合的过零检测法，但其最小二乘运算量大，后续仍需卡尔曼滤波校正.文献［5］将相关法与希尔伯特变换结合，提高了算法的抗干扰能力，但失去了希尔伯特变换的瞬时特性.文献［6-7］采用多次互相关、数据扩展和整周期化等手段抑制了非整周 期采样的影响.文献［8］的正交解调法（QD法）在批料流工况取得了较好的性能，但其使用的FIR陷波器在设计过程中需要预知陷波频率范围，解调信号的频率\omega一般设置为传感器的固有频率，因此该算法应用的前提是信号频率已知，且算法无法实现信号频率、相位差的联合估计.文献［9］提出将复带通滤波法（CBF法）应用到科氏流量计中，获得了快速的跟踪能力，但其并未针对小相位差（1°以下）和实际应用中流量检测情况开展实验.文献［10］的实复变换法将实信号转化为复信号，抑制了负频率分量对算法精度的影响.文献［11］采用Prony法进行参数估计，但复指数数量与流量计结构有关，因此其仅适用于特定类型流量计.文献［12］采用最小二乘法和曲线拟合的检测方法，在不降低精度或调整整采样周期的情况下，可以减少计算量，但其需经一段时间迭代才能达到所需的精度，且精度有限.文献［13-15］分别研究了锁相环法、QD法、DTFT法、随机梯度算法（SGA）、相关法和Hilbert变换法（HT法），验证了各个算法的优势.

为了解决现有算法存在的弊端，进一步提高科氏流量计信号处理方法各方面的性能，本文提出基于相关运算和Hilbert变换的正交解调法（简称CH-QD法），详细阐述了算法的基本原理及实现步骤.将该算法与目前工业领域实用效果较好的几种算法进行精度、稳定性、抗噪性、动态跟踪性等性能测试，并在基于TMS320F28335DSP研制的数字科氏流量计变送器系统上实现算法，进行性能分析和水流量标定实验，以验证算法的有效性与实用性.

1 CH-QD算法

传统的正交解调法（QD法） ^［8］通过人为生成与待测信号同频率的正余弦信号，将两路传感器信号分别乘以正余弦信号，将信号频率变换到2倍频，再通过低通滤波器滤除高频成分，保留低频分量，最后利用低频分量解算信号相位差.该算法对谐波的抑制能力强，但其产生同频参考信号的前提是信号频率已知，而信号频率也是科氏流量计的待测量，因此，该算法通常需要搭配额外的频率检测算法共同使用.为解决上述问题并实现信号参数的联合估计，本文提出CH-QD算法.

首先对任意一路传感器信号进行自相关运算获得其自相关信号，作为传感器信号的同频余弦参考信号；将自相关信号经过Hilbert变换获得Hilbert变换信号，作为传感器信号的同频正弦参考信号；根据自相关信号与Hilbert变换信号之间的同频、正交关系，由复数运算和三角函数运算解算出信号的瞬时频率；两路传感器信号分别经过上述正弦、余弦参考信号解调后，信号的中心频率变换到原始频率的2倍频，同时信号中包含着携带信号相位信息的低频成分；采用多级低通滤波器滤除高频成分，保留低频成分；最后利用三角函数提取两路信号的相位差.

具体实现如下：

设两路传感器信号x_{\mathrm{1}}\left(n\right)和x_{\mathrm{2}}\left(n\right)为

\left.\begin{array}{c}{x}_{\mathrm{1}}\left(n\right)=A_{\mathrm{1}}\mathrm{s}\mathrm{i}\mathrm{n}\left(\mathrm{2}\mathrm{\pi }f/f_{\mathrm{s}}n+{\theta }_{\mathrm{1}}\right)+{\xi }_{\mathrm{1}}\left(n\right),\\ x_{\mathrm{2}}\left(n\right)=A_{\mathrm{2}}\mathrm{s}\mathrm{i}\mathrm{n}\left(\mathrm{2}\mathrm{\pi }f/f_{\mathrm{s}}n+{\theta }_{\mathrm{2}}\right)+{\xi }_{\mathrm{2}}\left(n\right).\end{array}\right\}(1)

其中：A_{\mathrm{1}}和A_{\mathrm{2}}为信号幅度；f为信号频率；f{}_{\mathrm{s}}为采样频率；{\theta }_{\mathrm{1}}和{\theta }_{\mathrm{2}}为信号相位；{\xi }_{\mathrm{1}}\left(n\right)和{\xi }_{\mathrm{2}}\left(n\right)为高斯白噪声；N为采样点数，n=1，…，N.

根据相关运算原理，对于周期为T的正弦功率信号的自相关运算如式（2）所示，由于高斯白噪声与不同时刻的信号之间不存在相关性，即在时间差\tau \ne \mathrm{0}的情况下，其自相关为0，因此x_{\mathrm{1}}\left(n\right)和x_{\mathrm{2}}\left(n\right)的自相关信号中残留的噪声分量为0.同时自相关信号会保留原信号的频率信息，丢失掉相位信息.因此，对其中一路传感器信号进行自相关运算，得到同频余弦信号x_{\mathrm{3}}\left(n\right).

\begin{array}{l}{R}_{xx}\left(\tau \right)=\frac{\mathrm{1}}{T}{\int }_{\mathrm{0}}^{T}A\mathrm{s}\mathrm{i}\mathrm{n}\left(wt+\theta \right)\cdot A\mathrm{s}\mathrm{i}\mathrm{n}\left[w\left(t-\tau \right)+\theta \right]\mathrm{d}t=\\ \frac{\mathrm{1}}{T}{\int }_{-\frac{T}{\mathrm{2}}}^{\frac{T}{\mathrm{2}}}-\frac{A^{\mathrm{2}}}{\mathrm{2}}\left[\begin{array}{l}\mathrm{c}\mathrm{o}\mathrm{s}\left(\mathrm{2}wt-w\tau +\mathrm{2}\theta \right)-\\ \mathrm{c}\mathrm{o}\mathrm{s}\left(w\tau \right)\end{array}\right]dt=\\ \frac{A^{\mathrm{2}}}{\mathrm{2}}\mathrm{c}\mathrm{o}\mathrm{s}\left(w\tau \right),\end{array}(2)

x_{\mathrm{3}}\left(n\right)=\frac{A_{\mathrm{1}}^{\mathrm{2}}}{\mathrm{2}}\mathrm{c}\mathrm{o}\mathrm{s}\left(\mathrm{2}\mathrm{\pi }f/f_{\mathrm{s}}n\right).(3)

由于Hilbert变换在时域上定义为x（t）与\frac{\mathrm{1}}{\mathrm{\pi }t}的卷积，在频域上可等效于一个幅频特性为\left|H\left(\omega \right)\right|=\mathrm{1}，相频特性为 的全通滤波器，因此对余弦信号x_{\mathrm{3}}\left(n\right)进行Hilbert变换，可以产生同频正弦信号：

x_{\mathrm{4}}\left(n\right)=\frac{A_{\mathrm{1}}^{\mathrm{2}}}{\mathrm{2}}\mathrm{s}\mathrm{i}\mathrm{n}\left(\mathrm{2}\mathrm{\pi }f/f_{\mathrm{s}}n\right).(4)

由自相关信号x_{\mathrm{3}}\left(n\right)和Hilbert变换信号x_{\mathrm{4}}\left(n\right)构造复信号E\left(n\right)：

E\left(n\right)=x_{\mathrm{3}}\left(n\right)+\mathrm{j}{x}_{\mathrm{4}}\left(n\right)=\frac{A_{\mathrm{1}}^{\mathrm{2}}}{\mathrm{2}}{\mathrm{e}}^{\mathrm{j}\left(\mathrm{2}\mathrm{\pi }{f}_{\mathrm{0}}/f_{\mathrm{s}}n\right)}.(5)

复信号E\left(n\right)及其时移共轭信号E^{*}\left(n-\mathrm{1}\right)经过式（6），式（7）的数学运算和反三角函数运算，即可检测出信号的瞬时频率f_{\mathrm{0}}\left(n\right).

\begin{array}{l}Z\left(n\right)=E\left(n\right)\times E^{*}\left(n-\mathrm{1}\right)=\frac{A_{\mathrm{1}}^{\mathrm{2}}}{\mathrm{2}}{\mathrm{e}}^{\mathrm{j}\left(\frac{\mathrm{2}\mathrm{\pi }{f}_{\mathrm{0}}}{f_{\mathrm{s}}}\right)}=\\ Z_{\mathrm{r}}\left(n\right)+\mathrm{j}{Z}_{\mathrm{i}}\left(n\right),\end{array}(6)

f_{\mathrm{0}}\left(n\right)=\frac{f_{\mathrm{s}}}{\mathrm{2}\mathrm{\pi }}\times \mathrm{a}\mathrm{r}\mathrm{c}\mathrm{t}\mathrm{a}\mathrm{n}\frac{Z_{\mathrm{i}}\left(n\right)}{Z_{\mathrm{r}}\left(n\right)}.(7)

把两路信号x_{\mathrm{1}}\left(n\right)和x_{\mathrm{2}}\left(n\right)看作调频信号，用上述同频参考信号x_{\mathrm{3}}\left(n\right)和x_{\mathrm{4}}\left(n\right)分别对两路信号进行解调，以x_{\mathrm{1}}\left(n\right)为例，解调后的信号为x_{\mathrm{q}\mathrm{1}}\left(n\right)，x_{\mathrm{q}\mathrm{2}}\left(n\right)：

\begin{array}{l}{x}_{\mathrm{q}\mathrm{1}}\left(n\right)=x_{\mathrm{1}}\left(n\right)x_{\mathrm{4}}\left(n\right)=\\ \frac{A_{\mathrm{1}}^{\mathrm{3}}}{\mathrm{2}}\mathrm{s}\mathrm{i}\mathrm{n}\left(\mathrm{2}\mathrm{\pi }f/f_{\mathrm{s}}n+{\theta }_{\mathrm{1}}\right)\mathrm{s}\mathrm{i}\mathrm{n}\left(\mathrm{2}\mathrm{\pi }f/f_{\mathrm{s}}n\right)=\\ \frac{A_{\mathrm{1}}^{\mathrm{3}}}{\mathrm{4}}\left[\mathrm{c}\mathrm{o}\mathrm{s}\left({\theta }_{\mathrm{1}}\right)-\mathrm{c}\mathrm{o}\mathrm{s}\left(\mathrm{4}\mathrm{\pi }f/f_{\mathrm{s}}n+{\theta }_{\mathrm{1}}\right)\right]+{\tilde{\xi }}_{\mathrm{1}}\left(n\right),\end{array}(8)

\begin{array}{l}{x}_{\mathrm{q}\mathrm{2}}\left(n\right)=x_{\mathrm{1}}\left(n\right)x_{\mathrm{3}}\left(n\right)=\\ \frac{A_{\mathrm{1}}^{\mathrm{3}}}{\mathrm{2}}\mathrm{s}\mathrm{i}\mathrm{n}\left(\mathrm{2}\mathrm{\pi }f/f_{\mathrm{s}}n+{\theta }_{\mathrm{1}}\right)\mathrm{c}\mathrm{o}\mathrm{s}\left(\mathrm{2}\mathrm{\pi }f/f_{\mathrm{s}}n\right)=\\ \frac{A_{\mathrm{1}}^{\mathrm{3}}}{\mathrm{4}}\left[\mathrm{s}\mathrm{i}\mathrm{n}\left({\theta }_{\mathrm{1}}\right)+\mathrm{s}\mathrm{i}\mathrm{n}\left(\mathrm{4}\mathrm{\pi }f/f_{\mathrm{s}}n+{\theta }_{\mathrm{1}}\right)\right]+{\tilde{\xi }}_{\mathrm{2}}\left(n\right).\end{array}(9)

此时解调信号中包括低频分量、2倍频高频分量和残留干扰，经低通滤波滤去高频\mathrm{4}\mathrm{\pi }f/f_{\mathrm{s}}n+{\theta }_{\mathrm{1}}分量和残留干扰\tilde{\xi }\left(n\right)，得到同向分量I_{\mathrm{1}}和正交分量Q_{\mathrm{1}}：

\left.\begin{array}{c}{I}_{\mathrm{1}}=A_{\mathrm{1}}^{\mathrm{3}}/\mathrm{4}\cdot \mathrm{s}\mathrm{i}\mathrm{n}{\theta }_{\mathrm{1}},\\ Q_{\mathrm{1}}=A_{\mathrm{1}}^{\mathrm{3}}/\mathrm{4}\cdot \mathrm{c}\mathrm{o}\mathrm{s}{\theta }_{\mathrm{1}}.\end{array}\right\}(10)

以此类推，对信号x_{\mathrm{2}}\left(n\right)进行解调和低通滤波，得到同向分量I_{\mathrm{2}}和正交分量Q_{\mathrm{2}}：

\left.\begin{array}{c}{I}_{\mathrm{2}}=A_{\mathrm{1}}^{\mathrm{2}}{A}_{\mathrm{2}}/\mathrm{4}\cdot \mathrm{s}\mathrm{i}\mathrm{n}{\theta }_{\mathrm{2}},\\ Q_{\mathrm{2}}=A_{\mathrm{1}}^{\mathrm{2}}{A}_{\mathrm{2}}/\mathrm{4}\cdot \mathrm{c}\mathrm{o}\mathrm{s}{\theta }_{\mathrm{2}}.\end{array}\right\}(11)

利用式（12），式（13）进行三角函数运算得到同向分量I和正交分量Q，根据式（14）计算信号的实时相位差∆\phi.

\begin{array}{c}I=I_{\mathrm{1}}\cdot Q_{\mathrm{2}}-Q_{\mathrm{1}}\cdot I_{\mathrm{2}}=\\ \frac{A_{\mathrm{1}}^{\mathrm{3}}}{\mathrm{4}}\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_{\mathrm{1}}\cdot \frac{A_{\mathrm{1}}^{\mathrm{2}}{A}_{\mathrm{2}}}{\mathrm{4}}\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_{\mathrm{2}}-\frac{A_{\mathrm{1}}^{\mathrm{3}}}{\mathrm{4}}\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_{\mathrm{1}}\cdot \frac{A_{\mathrm{1}}^{\mathrm{2}}{A}_{\mathrm{2}}}{\mathrm{4}}\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_{\mathrm{2}}=\\ \frac{A_{\mathrm{1}}^{\mathrm{5}}{A}_{\mathrm{2}}}{\mathrm{16}}\mathrm{s}\mathrm{i}\mathrm{n}\left({\theta }_{\mathrm{1}}-{\theta }_{\mathrm{2}}\right),\end{array}(12)

\begin{array}{c}Q=Q_{\mathrm{1}}\cdot Q_{\mathrm{2}}+I_{\mathrm{1}}\cdot I_{\mathrm{2}}=\\ \frac{A_{\mathrm{1}}^{\mathrm{3}}}{\mathrm{4}}\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_{\mathrm{1}}\cdot \frac{A_{\mathrm{1}}^{\mathrm{2}}{A}_{\mathrm{2}}}{\mathrm{4}}\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_{\mathrm{2}}+\frac{A_{\mathrm{1}}^{\mathrm{3}}}{\mathrm{4}}\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_{\mathrm{1}}\cdot \frac{A_{\mathrm{1}}^{\mathrm{2}}{A}_{\mathrm{2}}}{\mathrm{4}}\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_{\mathrm{2}}=\\ \frac{A_{\mathrm{1}}^{\mathrm{5}}{A}_{\mathrm{2}}}{\mathrm{16}}\mathrm{c}\mathrm{o}\mathrm{s}\left({\theta }_{\mathrm{1}}-{\theta }_{\mathrm{2}}\right),\end{array}(13)

∆\phi ={\theta }_{\mathrm{1}}-{\theta }_{\mathrm{2}}=\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{t}\mathrm{a}\mathrm{n}\frac{I}{Q}.(14)

由于滤波器的设计直接影响到正交解调法的精度和动态跟踪性能，且不同类型传感器的信号频率之间存在差异，如果采用固定参数的滤波器难以同时满足通频带较窄、与各类传感器适配.综合考虑，本文采用4级无限脉冲响应（IIR）低通滤波的方案，滤波器系数LP1和LP2由式（15），式（16）设计，该系数只与设定的滤波器截止频率f_{\mathrm{c}}和信号的采样频率f_{\mathrm{s}}有关，设计简便，易于在嵌入式系统中实时实现.其中，前2个滤波器的默认截止频率选取为45.5 Hz，后2个滤波器的默认截止频率选取为10 Hz，首次采用CH-QD算法检测信号参数时采用默认截止频率，随后，根据检测到的传感器信号频率f_{\mathrm{0}}实时调整截止频率f_{\mathrm{c}}，前2个滤波器的截止频率设置为0.2f_{\mathrm{0}}，后2个滤波器的截止频率设置为0.05f_{\mathrm{0}}，以达到更好的滤波效果.

\mathrm{L}\mathrm{P}\mathrm{1}=\frac{f_{\mathrm{s}}}{\left(f_{\mathrm{s}}+\mathrm{2}\mathrm{\pi }{f}_{\mathrm{c}}\right)},(15)

\mathrm{L}\mathrm{P}\mathrm{2}=\mathrm{1}-\mathrm{L}\mathrm{P}\mathrm{1}.(16)

当对第k个采样数据滤波时，输入x\left(k\right)的滤波输出y\left(k\right)即为

y\left(k\right)=\mathrm{L}\mathrm{P}\mathrm{1}\cdot y\left(k-\mathrm{1}\right)+\mathrm{L}\mathrm{P}\mathrm{2}\cdot x\left(k\right).(17)

算法原理框图见图1，与传统QD法相比，本文算法无需提前预知信号频率，实现了CMF所需的频率、相位差的联合估计，同时滤波器参数的设计简单有效，多级滤波在保证实时性的基础上进一步增强了算法的抗噪性.

2 仿真验证

工业应用中，科氏流量计两路传感器信号在各个稳定流量点上可视为频率相同、相位差与流量成正比的正弦波，信号中普遍夹杂着50 Hz工频、随机噪声、二次和三次固有频率谐波的干扰.为此建立基本信号模型：

\left.\begin{array}{c}{x}_{\mathrm{1}}\left(n\right)=\mathrm{s}\mathrm{i}\mathrm{n}\left(\mathrm{2}\mathrm{\pi }{f}_{\mathrm{0}}/f_{\mathrm{s}}n+{\theta }_{\mathrm{1}}\right)+e_{\mathrm{1}}\left(n\right),\\ x_{\mathrm{2}}\left(n\right)=\mathrm{s}\mathrm{i}\mathrm{n}\left(\mathrm{2}\mathrm{\pi }{f}_{\mathrm{0}}/f_{\mathrm{s}}n+{\theta }_{\mathrm{2}}\right)+e_{\mathrm{2}}\left(n\right).\end{array}\right\}（18）

其中：{\theta }_{\mathrm{1}}和{\theta }_{\mathrm{2}}是信号相位；e_{\mathrm{1}}\left(n\right)和e_{\mathrm{2}}\left(n\right)为干扰项，包含相对幅度10%的工频，二次、三次谐波和信噪比可调的高斯白噪声.不同类型科氏流量计传感器的固有频率和扭转角存在差异，一般固有频率在60~1 000 Hz范围内，最大量程质量流量对应的相位差小于4°.为了适配各类传感器、保证采样信号的还原度，进一步提升仪表的测量精度，本文的采样频率f_{\mathrm{s}}选取为100 kHz，并模拟现有一款U型传感器的振动特性，将仿真信号基频设定为84.5 Hz，幅值1 V，采样点数N=8 192.

信噪比为10 dB时，上述信号经过滤波后，输出信号的频谱见图2.本文对解调信号采用逐级滤波的方案，前2次滤波可以很大限度地消除高频成分，后2次滤波除进一步消除高频成分外，还可以滤除低频干扰，改善信号质量.

2.1 算法抗噪性能测试

为验证算法的抗噪性能，基于基本信号模型式（18），{\theta }_{\mathrm{1}}和{\theta }_{\mathrm{2}}分别设置为30°和30.2°，在0~80 dB不同信噪比下，对本文算法与QD法、CBF法、HT法和DTFT法进行500次独立的仿真实验，获得每次实验中逐点频率和相位差结果的均方误差MSE，500次实验的平均均方误差MSE如图3，图4所示.

QD法中先利用常用的自适应滤波器检测并跟踪信号频率，据此生成参考信号，并采用常用的椭圆低通IIR滤波器进行解调信号的滤波；HT法中采用常用的加窗函数法来抑制其端点效应，并舍弃检测结果中两端各100个畸变的数据点；DTFT法中采用重叠汉宁窗的方法来抑制其频率泄漏问题.为模拟科氏流量计信号处理流程，不失一般性，先由FIR滤波器对信号进行初步降噪，以保证与实际应用情况相符.

图3频率检测结果表明：HT法在信噪比较低的条件下稳定性较差，随着信噪比增大逐渐趋于稳定，但与CH-QD算法相比，始终波动更剧烈；CBF法自身的滤波特性使其几乎不受信噪比影响，在低信噪比环境下频率检测MSE仍然小于10^-3，较CH-QD算法大10倍；而CH-QD算法内部的相关运算和多级滤波环节使其在低信噪比下频率检测MSE依旧小于10^-4，算法对噪声的敏感度较低，频率估计的稳定性更高、抗噪性更强.

图4相位差检测结果表明，5种相位差估计算法的MSE均随信噪比的增大而呈减小的趋势，并在信噪比增大到某个程度后趋于稳定.从稳定之前来看，CD-QD算法的MSE较其他算法小2个数量级，在低信噪比环境下的稳定性优于其他算法；从趋于稳定时的信噪比值来看，CH-QD法在信噪比50 dB时就可以达到其最佳运算的稳定状态，而其他算法需在信噪比60 dB时才能达到稳定状态.因此，CH-QD算法对信噪比的要求较低，在较低的信噪比下就可以达到算法的最优估计状态，证明其具有较佳的稳定性和抗噪性.

2.2 算法相位差检测精度测试

为验证算法的精度，以信号模型式（18）作为信号源，在30 dB随机噪声下，相位差范围为0.001°~4°，对上述5种算法进行500次仿真实验，不同相位差下算法相位差检测结果的相对误差和均方误差MSE的均值如图5，图6所示.

由图可知，5种方法的相对误差均随着相位差的增大而呈现减小的趋势，并在相位差增大到某一特定值后保持稳定.其中，CH-QD算法的相对误差始终小于其他算法，在相位差为0.01°时仍能保持在0.04%以下，并且由相位差大小导致的相对误差波动较小，算法对相位差大小的敏感度较低.

5种方法的MSE均随着相位差的增大而呈现增大的趋势.其中，CH-QD算法的MSE明显小于其他算法，在小相位差的检测上优势更为显著，而小相位差检测正是科氏流量计低流速测量所面临的难点.因此，CH-QD算法在不同相位差下的检测精度、稳定性均优于其他算法，算法的精度可以满足科氏流量计的测量需求.

2.3 算法相位差检测动态跟踪性能测试

在流体流速变化等条件下，科氏流量计传感器信号受流体特性和流体状态的影响，其幅度、频率和相位差均随时间变化呈现小范围、不规则的变化，表现为基本信号模型式（18）的幅度、频率和相位参数在小幅度范围内的瞬时波动，即为时变信号.这就要求信号处理方法在保证测量精度和稳定性的前提下，能够及时跟踪信号参数的变化，具有较小的运算量和较快的响应速度.并且，由于信号频率和相位同时变化，依赖于频率检测精度的相位差检测算法就不再适用.因此，目前许多信号处理方法均无法达到时变信号的测量需求.

为模拟时变信号特征，文献［16］提出使用随机游动模型（RWM）生成时变信号，但Li等^［9］指出该模型并未对参数的瞬时变化率采取任何限制，而实际CMF的参数受机械惯性、有限流体速度等因素限制不会产生快速的变化.因此，结合实际信号的研究，本文采用其提出的修正的随机游动模型（MRWM）来表征时变信号，如式（19）所示.由边界限制幅度、频率和相位差生成随机序列A\left(n\right)，f\left(n\right)，∆\theta \left(n\right)后，通过低通滤波以限制其变化率，由滤波后的参数A_{\mathrm{l}}\left(n\right)，f_{\mathrm{l}}\left(n\right)，∆{\theta }_{\mathrm{l}}\left(n\right)生成两路时变信号，产生的信号参数见图7.

\left.\begin{array}{l}A\left(n\right)=\frac{A_{\mathrm{m}\mathrm{i}\mathrm{n}}+A_{\mathrm{m}\mathrm{a}\mathrm{x}}}{\mathrm{2}}+\frac{A_{\mathrm{m}\mathrm{a}\mathrm{x}}-A_{\mathrm{m}\mathrm{i}\mathrm{n}}}{\mathrm{2}}\cdot e_A\left(n\right),\\ f\left(n\right)=\frac{f_{\mathrm{m}\mathrm{i}\mathrm{n}}+f_{\mathrm{m}\mathrm{a}\mathrm{x}}}{\mathrm{2}}+\frac{f_{\mathrm{m}\mathrm{a}\mathrm{x}}-f_{\mathrm{m}\mathrm{i}\mathrm{n}}}{\mathrm{2}}\cdot e_f\left(n\right),\\ ∆\theta \left(n\right)=\frac{∆{\theta }_{\mathrm{m}\mathrm{i}\mathrm{n}}+∆{\theta }_{\mathrm{m}\mathrm{a}\mathrm{x}}}{\mathrm{2}}+\frac{∆{\theta }_{\mathrm{m}\mathrm{a}\mathrm{x}}-∆{\theta }_{\mathrm{m}\mathrm{i}\mathrm{n}}}{\mathrm{2}}\cdot e_{∆\theta }\left(n\right),\\ A_{\mathrm{l}}\left(n\right)=h_A\left(n\right)\mathrm{*}A\left(n\right),\\ f_{\mathrm{l}}\left(n\right)=h_f\left(n\right)\mathrm{*}f\left(n\right),\\ ∆{\theta }_{\mathrm{l}}\left(n\right)=h_{∆\theta }\left(n\right)\mathrm{*}∆\theta \left(n\right),\\ y_{\mathrm{1}}\left(n\right)=A_{\mathrm{l}}\left(n\right)\mathrm{s}\mathrm{i}\mathrm{n}\left[\mathrm{2}\mathrm{\pi }{f}_{\mathrm{l}}\left(n\right)/f_{\mathrm{s}}n+∆{\theta }_{\mathrm{l}}\left(n\right)/\mathrm{2}\right]+\\ {\sigma }_{\mathrm{e}\mathrm{1}}\cdot e_{\mathrm{1}}\left(n\right),\\ y_{\mathrm{2}}\left(n\right)=A_{\mathrm{l}}\left(n\right)\mathrm{s}\mathrm{i}\mathrm{n}\left[\mathrm{2}\mathrm{\pi }{f}_{\mathrm{l}}\left(n\right)/f_{\mathrm{s}}n-∆{\theta }_{\mathrm{l}}\left(n\right)/\mathrm{2}\right]+\\ {\sigma }_{\mathrm{e}\mathrm{2}}\cdot e_{\mathrm{2}}\left(n\right).\end{array}\right\}（19）

其中：e_A，e_f，e_{∆\theta }，e_{\mathrm{1}}(n)，e_{\mathrm{2}}(n)是非相关的白噪声序列；{\sigma }_{\mathrm{e}\mathrm{1}}，{\sigma }_{\mathrm{e}\mathrm{2}}是噪声增益因子，本文设置为10^-5；h_A\left(n\right)，h_f\left(n\right)，h_{∆\theta }\left(n\right)是低通滤波器.边界限制分别设置为A_{\mathrm{m}\mathrm{i}\mathrm{n}}=0.95 V，A_{\mathrm{m}\mathrm{a}\mathrm{x}}=1.05 V，f_{\mathrm{m}\mathrm{i}\mathrm{n}}=80 Hz，f_{\mathrm{m}\mathrm{a}\mathrm{x}}=89 Hz，∆{\theta }_{\mathrm{m}\mathrm{i}\mathrm{n}}=0°，∆{\theta }_{\mathrm{m}\mathrm{a}\mathrm{x}}=0.4°.

采用上述MRWM模型对5种算法进行测试，逐点相位差跟踪结果见图8，跟踪结果与真实频率、真实相位差之间的MSE见表1.结果表明，HT法的检测结果在实际相位差上下持续、大幅度波动，基本无法使用；QD法受其算法特性的影响，实时检测中存在明显的延迟，动态跟踪性能较差，对于相位差的微小波动敏感性低，不适合应用于时变信号条件下；DTFT法具有一定的响应速度和跟踪性能，但仍存在检测偏差较大的情况；本文提出的CH-QD法的动态跟踪性能和检测精度均表现良好，更适用于作为非平稳流况下CMF的信号处理方法.

3 系统实现

以DSP+ARM双核架构研制CMF变送器，选用TI公司的TMS320F28335DSP作为主控芯片搭建数字信号处理模块，控制AD7606芯片以100 kHz采样率完成传感器双极性信号16位同步采样，并以STM32F103为核心设计驱动模块，提供传感器充足的振动能量.整套系统硬件主要包括信号调理和AD采集模块、驱动模块、数字信号处理模块、人机交互模块和温度补偿模块.首先，ADC采集的两路传感器信号经过一个带宽较大的FIR带通滤波器滤波后，由信号处理算法计算信号的频率、相位差等关键参数，结合温度补偿、仪表系数计算出流体的瞬时质量流量、累计质量流量、密度等信息.与此同时，变送器完成数据传输、人机交互等仪表功能.

在所研制的变送器系统下，同等条件下上述5种算法的运算时间见表2.其中，CBF法的运算速度较快，其大部分运算都集中在滤波和复数运算环节；QD法的运算量集中在低通IIR滤波环节，运算时间较长；HT法内部不包含任何滤波环节，计算量小，运算速度最快，但其精度、稳定性和抗噪性较差，不宜用于时变信号条件下；DTFT法需采用重叠汉宁窗抑制频率泄漏问题导致其每个采样点都需要计算窗口内数据的DTFT以保证精度，因此计算量最大；CH-QD法的运算环节主要包括相关运算、Hilbert变换、解调和滤波.较为复杂，但其运算速度仍取决于滤波环节.而本文所采用的多级滤波方案极大地降低了运算量，使得算法运算速度较快，完全能够满足科氏流量计实时测量的需求.

4 标定实验

为了验证算法在实际工况条件下的效果，选用科隆公司生产的DN15微弯型科氏流量计传感器，满管固有频率为430 Hz，匹配本文所研制的变送器，采用称重法进行水流量标定实验，实验装置主要由PLC控制柜、水箱、水泵、待测流量计、称重装置、换向器等组成，如图9所示.整个装置通过WinCC监控软件对各部分器件进行控制，并进行自动标定与数据记录.其中，电子秤的精度为±0.02%，不确定度为0.012%，装置的扩展不确定度为0.023 2%.

实验过程中，先将换向器调到1路，调节手动调节阀使流量达到最大，冲水一段时间，以消除流体挂壁现象和温差引起的误差.关闭上游阀和下游阀，使传感器处于满管零流量状态，运行程序一段时间，直至零点稳定后进行零点校正.换向器切换至2路开始标定，打开阀门、调节手动调节阀至标定的流速，在平稳流速下标定5个流量点，流量范围为1.83～37.41 kg·min^-1，每个流量点测试3组数据.水经由水箱、水泵、科氏流量计流向称重装置，科氏流量计变送器将检测结果以PWM脉冲的形式上传至PLC控制柜，由WinCC监控软件进行累计和标定.

标定结果如表3所示，在量程比20∶1范围内，本文研制的变送器系统搭配所提出算法实现的质量流量测量误差小于0.1%，重复性小于0.05%，符合国家工业0.1级科氏流量计标准，具有较好的实用价值.

5 结 论

为了提高科氏流量计信号处理方法各方面的性能，本文提出CH-QD算法进行信号频率、相位差检测，克服了传统QD法需预知信号频率、对频率波动敏感等问题，所采用的滤波方案极大程度上提高了算法的运算速度.算法性能测试分析结果表明，本文算法在低信噪比、小相位差和时变信号下均表现突出，具有较高的检测精度、稳定性、抗噪性和动态跟踪性能.

在基于TMS320F28335DSP研制的科氏流量计变送器上实现本文算法，并进行水流量标定实验，标定结果表明，在量程比20∶1的范围内，系统测量误差优于0.1%，重复性优于0.05%，整套变送器系统性能良好，具有较高的实用价值.

参考文献

[1] Romano P. Coriolics mass flow rate meter having a substantially increased noise immunity： US Patent， No 4934196［P］. 1990-06-19.

[2] Hou Q L， Xu K J. A DSP-based signal processing method and system for CMF［J］. Measurement，2013，46（7）：2184–2192.

[3] 刘陈慈，徐科军，黄雅. 尿素机中高频微弯型科氏质量流量计信号处理［J］. 计量学报，2022，43（6）：767-776.

[4] Liu Chen-ci， Xu Ke-jun， Huang Ya. Signal processing of high-frequency micro-bending Coriolis mass flowmeter in urea machine ［J］. Journal of Metrology，2022，43 （6）：767-776.

[5] 徐长利，李波，钟其水，等. 基于过零检测的相位差计算方法的改进［J］. 机电信息，2022（2）：35-38.

[6] Xu Chang-li， Li Bo， Zhong Qi-shui，et al. Improvement of phase difference calculation method based on zero-crossing detection ［J］. Electromechanical Information，2022（2）：35-38.

[7] Liu Y J， Wang L Y， Luo Q，et al. A novel method of Coriolis flowmeter phase difference measurement based on improved correlation analysis method［J］. Flow Measurement and Instrumentation，2021， 80：1-11.

[8] Shen T A， Huang S Y， Chen P，et al. A novel phase difference measurement method for Coriolis mass flowmeter based on correlation theory［J］. Energies，2022，15（3710）：1-11.

[9] 沈廷鳌，涂亚庆，李明，等. 数据延拓式相关的相位差测量方法及验证［J］. 仪器仪表学报，2014，35（6）：1331-1337.

[10] Shen Ting-ao， Tu Ya-qing， Li Ming，et al. Phase difference measurement method and verification of data continuation correlation ［J］. Journal of Instrument，2014，35（6）：1331-1337.

[11] 张伦，徐科军，徐浩然，等. 面向批量流测量的科氏质量流量计正交解调方法实现［J］. 计量学报，2020，41（7）：808-815.

[12] Zhang Lun， Xu Ke-jun， Xu Hao-ran，et al. Implementation of orthogonal demodulation method for Coriolis mass flowmeter for batch flow measurement ［J］. Journal of Metrology，2020，41（7）：808-815.

[13] Li M， Henry M. Complex bandpass filtering for coriolis mass flow meter signal processing［C］// 42nd Annual Conference of the IEEE Industrial Electronics Society（IECON）. Florence，2016：4952-4957.

[14] Chen P， Tu Y Q， Li M，et al. A real-to-complex conversion phase difference estimation method for Coriolis mass flowmeter signal［C］ // 2019 International Conference on Communications， Information System and Computer Engineering （CISCE）. Haikou，2019：280-284.

[15] Semenov A， Fedosov I. Novel prony-based algorithm for estimating oscillation parameters of Coriolis flowmeter at two-phase flow［C］// XXIX International Scientific Symposium Metrology and Metrology Assurance （MMA）. Sozopol，2019：1-6.

[16] Lei H. An improved phase difference detection method for a Coriolis flowmeter［J］. Measurement，2021，172：1-7.

[17] Li M， Henry M. Signal processing methods for Coriolis mass flow metering in two-phase flow conditions［C］// IEEE International Conference on Industrial Technology（ICIT）. Taipei，2016：690-695.

[18] Li C H， Sun L J. Improvement of signal processing in Coriolis mass flowmeters for gas-liquid two-phase flow［J］. Frontiers of Information Technology and Electronic Engineering，2021，22（2）：272-286.

[19] Chen N， Xiao Y， Liu Z，et al. Research on Coriolis mass flowmeter signal processing and driving method［C］ // 6th International Conference on Computing Methodologies and Communication （ICCMC）. Erode，2022：589-592.

[20] Tu Y Q， Yang H Y， Zhang H T，et al. CMF signal processing method based on feedback corrected ANF and Hilbert transformation［J］. Measurement Science Review，2014，14（1）：41-47.

基金资助

辽宁省科学技术计划项目(2021JH1/10400011)

辽宁省科学技术计划项目(2020JH2/10100031)