科里奥利流量计（以下简称科氏流量计）具有测量精度高、能直接得到流体的质量流量等优点，在天然气计量领域应用广泛[1-2]。相位差算法是科氏流量计实现流量准确计量的关键之一，目前，相位差算法主要有离散时间傅里叶变换（DTFT）法、相关法、正交解调法、希尔伯特变换法等[3-6]。DTFT法考虑负频率的影响，计算出各采样点处的傅里叶系数，计量精度高，但需要知道信号频率，且存在迭代步骤，实时性不足。相关法利用噪声与噪声，噪声与信号互不相关的特性计算相位差，抗干扰能力强，但在非整周期采样时会引入较大误差。正交解调法通过构造同频正、余弦信号与测量信号进行解调，再用低通滤波后的分量计算出相位差，对低通滤波器的设计要求高。希尔伯特变换法没有迭代步骤，计算量小，适用于时变信号的相位差估计，但存在端点效应，降低了相位差计算精度。本研究通过分析科氏流量计输出信号的非线性特征，提出了一种基于粒子滤波算法的科氏流量计信号处理方法。
科氏流量计的测量原理是科里奥利效应，通过设计不同弯曲形状的测量管产生科里奥利效应以实现质量流量的测量，以常见的U型管科氏流量计来说明其测量原理（见图1）[7]。如图1所示，在激振线圈的作用下，U型管绕OO轴上下摆动，当没有流体通过时，测量管两路输出信号的频率和相位均相同。当有流体通过时，流体在直管段运动时不断碰撞管壁，管壁会受到流体的反作用力，这个反作用力就是科里奥利力（以下简称科氏力）。
在测量管两侧，$ {F}_{{\mathrm{c}},1} $与$ {F}_{{\mathrm{c}},2} $大小相等，方向相反。科氏力的计算公式见式（1）。
式中：Fc为科氏力，N；m为流体质量，kg；v为流体流速，m/s；$\omega$为角频率，rad/s。
测量管直管段所受科氏力变化量计算公式如式（2）所示。
式中：$\Delta {F}_{{\mathrm{c}}} $为科氏力变化量，N；$ {q}_{m} $为质量流量，kg/h，$ {q}_{m}=\rho vS $，$ \rho $为流体密度，kg/m3，S为测量管截面积，m2；$ \Delta L $为直管段上一单位长度，m。
测量管直管段一单位微小长度的扭矩，其计算公式见式（3）。
式中：$ {\mathrm{d}}M $为直管段上一单位微小长度的扭矩，N∙m；r为弯管段半径，m；dFc为直管段一单位微小长度上所受科氏力，N；dL为直管段一单位微小长度，m。
对$ \mathrm{d}M $进行积分变换，可得到整个直管段受到的扭矩，其计算公式见式（4）。
式中：M为整个直管段受到的扭矩，N·m；L为测量长度，m。
在该扭矩作用下，U型管产生一微小的扭转角的计算公式见式（5）。
式中：$ {K}_{{\mathrm{s}}} $为测量管角弹性模量，Pa/rad；θ为扭转角，rad。
由式（4）和式（5）可得式（6）。
在形成$ \theta $的时间内，θ计算公式见式（7）。
式中：$\Delta t $为形成θ的时间，s；$ {v}_{{\mathrm{l}}} $为流量管在振动方向上的线速度，$ {v}_{{\mathrm{l}}} $=ωL，m/s。
由式（6）和式（7）联立得出质量流量计算公式，见式（8）。
由式（8）可知，质量流量与时间差成正比。
粒子滤波适用于强非线性系统，广泛应用于目标定位、视觉跟踪等领域。与卡尔曼滤波一样，粒子滤波算法使用状态空间模型来描述待测系统[8]，其方程式见式（9）与式（10）。
式中$ {:X}_{k}\mathrm{为}k $时刻的状态估计；f为状态函数；$ {Z}_{k} $为$ k $时刻的测量值；h为观测函数；$ {W}_{k-1} $和$ {V}_{k} $分别为测量系统的过程噪声和测量噪声；Xk-1为在k−1时刻的状态向量。
根据科氏流量计的信号特点，建立科氏流量计测量系统的状态空间模型，其状态矢量、状态转移方程和观测方程如式（11）~式（14）所示[9-10]。
式中：Xi|k为在时间为k时关于随机变量Xi的条件概念分布；$ {{\textit{z}}}_{k} $为传感器测量值；$ {A}_{k} $、$ {\theta }_{k} $分别为信号幅值和相位；$ {\omega}_{k} $为数字角频率，$ \omega_k=2\mathrm{\pi}f_k/f\mathrm{_s} $，fk为信号在k时刻的频率，fs为采样频率，kHz；$ F $为系统的状态转移矩阵。
粒子滤波算法是一种基于蒙特卡洛仿真的近似贝叶斯滤波算法[11]，其核心思想是用一些离散随机采样点来近似系统随机变量的概率密度函数，以样本均值代替积分运算，从而获得状态的最小方差估计。重采样是粒子滤波算法的关键环节，在保证粒子数量不变的情况下，通过抛弃权重低的粒子，多复制权重高的粒子而防止粒子退化，以提高算法的估计精度。
系统重采样算法由Doucet提出[12]，它与分层采样算法相似，算法步骤如下：
1） 将（0,1]划分为n个连续的互补区间，即：$ （\mathrm{0,1}]=（\mathrm{0,1}/n]\cup \ldots \cup ((n-1)/n,1] $。
2） 对每个子区间独立同分布采样得到$ {U}_{i} $，其计算公式见式（15）。
式中：Ui为重采样过程中用于选择粒子的随机数；$ U\left(\left(0,\dfrac{1}{n}\right]\right) $为区间$ \left(0,\dfrac{1}{n}\right] $上的均匀分布。
3） Ii为粒子的索引，用于标识n特定的粒子，令$ {I}_{i}=\mathrm{c}\mathrm{d}\mathrm{f}\left({u}_{i}\right) $，其中cdf是权值集合$ \{{w}_{i}{\}}_{i=1}^{n} $的累积分布函数。即对于$ u\in \left(\displaystyle\sum _{j=1}^{i-1}{w}^{j},\displaystyle\sum _{j=1}^{i}{w}^{j}\right] $，$ \mathrm{c}\mathrm{df}\left(u\right)=i $。设$ \gamma \left(i\right)={\gamma }_{i} $满足函数映射，$ \gamma :\{\mathrm{1,2},\dots, m\}\to X $，则$ {\gamma }_{i} $可以表示为$ \gamma \mathrm{c}\mathrm{d}\mathrm{f}\left({u}_{i}\right) $。
4） 初始化权值$ {w}_{j}=1/n $。记{$ {v}_{i}{\}}_{i=1}^{n} $为重采样后对应粒子复制数目的集合，其中$ {v}_{i} $表示重采样前的第$ i $个粒子在重采样后被复制的数目，$ {0 < v}_{i}\le m $。
标准粒子滤波算法通过重采样算法来解决序贯重要性采样存在的粒子退化问题，利用重要性概率密度函数$ q({x_k}|{{\textit{z}}_{1:k}}) $来抽取样本。具体步骤如下[13]：
1） 根据重要性概率密度函数抽取n个新粒子,其具体抽取方法见式（16），式（16）表示每个新粒子xi是通过抽样重要性概率密度函数q（xk|xk-1,${\textit{z}}_k $）得到的，其中i取1~n。
式中：$x_{i,k}\sim q $表示每个新粒子x是通过抽样条件概念密度函数q得到的；${q(x_{i,k}|x_{i,k - 1},{{\textit{z}}_k})} $为用于生成新粒子的概率密度函数，作为抽取粒子的依据；$ x_{i,k} $为得到的新粒子；${\textit{z}}_k $为k时刻的观测值。
2） 为每一个粒子计算权重，其计算公式见式（17）。
式中：$ w_{i,k} $为粒子权值；$ p({{\textit{z}}_k}|x_{i,k}) $为观测量，是${\textit{z}}_k $时刻的似然概率密度函数；$ p(x_{i,k}|x_{i,k - 1}) $为具有先验性质的系统状态转移概率密度函数。
3） 归一化权重，见式（18）。
4） 计算有效粒子数，见式（19）。
式中：Neff为有效粒子数。
5） 判断是否进行重采样，见式（20）。
式中：Nth为规定阈值。
若Neff小于阈值Nth，表示出现了粒子退化，须重新采样更新粒子集。
6） 得到当前时刻的状态和方差估计，返回步骤2进行下一时刻的更新，见式（21）与式（22）。
式中：$ {\hat x_k} $和$ {p_k} $分别为k时刻的状态和方差估计；T为状态转移矩阵。
在流量测量过程中，科氏流量计振动信号的幅值、频率和相位参数会发生微小的变化。为了模拟实际情况下的科氏流量计输出信号，在Matlab中使用如下的随机游动模型进行相位差算法仿真[14]，见式（23）~式（26）。
式中：y（n）为信号，V；A（n）为幅值，V；ω（n）为角频率，rad/s；φ（n）为相位，rad；$ e\left(n\right) $为在时间点n处的高斯白噪声，V；$ {e}_{A}\left(n\right) $、$ {e}_{\omega}\left(n\right) $、$ {e}_{\varphi }\left(n\right) $分别为影响幅值、角频率和相位的高斯白噪声分量；$ {\sigma }_{A} $、$ {\sigma }_{\omega} $、$ {\sigma }_{\varphi } $为调节系数，用来改变相应参数的游动幅度。
随机游动模型的仿真参数设置为：$ A\left(0\right)=10\;{\mathrm{V}} $，$\varphi(0)=0 $，$ \omega\left(0\right)=0.314 \;2 $，$ {f}_{{\mathrm{s}}}=20 \;\mathrm{kH}\mathrm{z} $，$ {\sigma }_{A}={10}^{-3} $，$ {\sigma }_{\omega}= {10}^{-5}，{\sigma }_{\varphi }= {10}^{-5} $ 。使用希尔伯特变换法计算结果作为本方法的对照，图2为两种方法的相位差估计曲线，可看出粒子滤波算法经过短时间的收敛后追踪到真实的相位差，能较好地反映相位差的变化情况。而希尔伯特变换法的计算结果围绕真实值上下浮动，且在端点效应的影响下，相位差计算结果在两端点处发散，降低了相位差的计算精度。图3为两种方法采样点处相位差的相对误差，可以看出粒子滤波算法相对误差更小，估计结果更接近真实值。
在不同初始相位差下，使用相对误差均值和均方误差来比较两种计算方法的性能。其中，均方误差按式（27）计算。
式中：SME为均方误差；$ \hat{\theta }（i） $为相位差计算值；$ \theta （i）$为相位差真实值。
由于粒子滤波算法存在一定的收敛过程，故均取50点以后的数据计算。
表1所列为不同初始条件下相位差算法的测量误差。由表1可知，粒子滤波算法的相对误差均值小于希尔伯特变换法，均方误差也降低了一个数量级，说明粒子滤波算法的相位差计算精度更高。
科氏流量计用于高压天然气的流量计量尚处于探索阶段，测试方法还不完善。目前，用于天然气计量的科氏流量计大多采用水作为介质进行测试验证，因此，本研究以水为介质进行相关实验研究，在水流量标准装置上开展算法测试实验。以称量法流量标准装置提供的流量作为参考流量，装置电子秤的测量范围为0~1 000 kg，精度为±0.01 kg，不确定度为0.05%（k=2）。通过高速数据采集卡采集科氏流量计的原始振动信号，采样频率为20 kHz，每次使用20 000个数据进行计算。
实验过程为：开启给水泵，进行水流量循环，通过调节阀将流量调至所需的流量点，待流量稳定后启动称量装置以获得参考流量，同时采集科氏流量计的原始振动信号；切换流量点，采集新流量点下科氏流量计的原始振动信号；在平稳流量下，利用累计流量和时间差的线性关系，通过线性回归得到拟合直线，进而得到参考时间差（Δt参考），并与采用希尔伯特变换法和粒子滤波算法计算出的时间差（分别以Δt希尔伯特和Δt粒子滤波表示）进行比较（见表2）。
从表2可以看出，与希尔伯特变换法相比，粒子滤波算法得到的时间差更接近参考值，说明基于粒子滤波的计算方法可有效地应用于科氏流量计实测振动信号的处理，可提升时间差或者相位差的测量精度。
本研究提出了一种基于粒子滤波的科氏流量计相位差算法。根据科氏流量计输出信号的特点，建立了测量系统的状态方程，采用粒子滤波算法对状态值进行迭代更新，获得每个采样点处的实时相位差。仿真分析结果表明，与希尔伯特变换法计算的相位差相比，粒子滤波算法的精度显著提升。通过实验系统采集了科氏流量计原始实际振动信号，处理结果表明，基于粒子滤波的计算方法可有效地应用于科氏流量计实测振动信号的处理，并在一定程度上提升了时间差或者相位差的测量精度。
2024, Vol. 53 Issue (6): 109-113 2024, Vol. 53 Issue (6): 109-113