CN103162755B - 一种基于改进自适应算法的科里奥利流量计信号跟踪方法 - Google Patents

Abstract

本发明涉及一种基于改进自适应算法的科里奥利流量计信号跟踪方法，解决现有信号跟踪方法难以长时间实时准确跟踪和测量信号频率的问题。本发明基于奇异谱分析方法和基于优化零极点配置滤波算法，通过对输入信号进行预处理，构建以

Description

一种基于改进自适应算法的科里奥利流量计信号跟踪方法

技术领域

本发明属于流体测量领域，具体涉及到一种基于改进自适应算法的科里奥利流量计信号跟踪方法。

背景技术

流体质量流量的测量与控制在许多工业领域有着重要而广泛的作用。质量流量测量包括直接和间接两种方式。由于技术方面的原因，传统的质量流量测量多采用间接方式。间接测量质量流量需要测压力、温度、体积流量、密度等，这使得流量计部件增多，从而使流量计成本高、精度低、故障率高，安装使用维护很不方便。目前在工业计量中，对于温度、压力、组分不断变化的气体介质来说，直接方式测量远比间接测量更具有实际价值。因此，直接测量流体的质量流量比间接测量质量流量有更大的优势。

科里奥利流量计已成为目前研究最多、最有前途的直接式质量流量仪器，市场需求量也越来越大。对科里奥利流量计的信号进行处理，关键在于能够精确的测算出来自两路传感器输出信号的频率和相位。而对于测量管输入段与输出段两路信号间相位差的精确测量则首先要求能够精确测得这两路信号的频率。因此，需要对传感器输出信号的频率进行实时、精确的跟踪和测量。

现有的基于自适应陷波器方法的频率测量可以根据被处理信号的特点，自动调整陷波器自身模型，使其陷波频率收敛到信号的基频处，并可由陷波器的参数求出基频，以实现频率的实时测量与跟踪，但存在难以兼顾频率跟踪的精度与跟踪的实时性的问题。另有基于格型IIR滤波器的自适应陷波器，其计算较为复杂，难以长时间持续跟踪信号频率的变化。也有采用的简化格型自适应陷波器，虽然计算简单、便于硬件实现，但其对于先验知识未知的信号，难以兼顾收敛速度和频率跟踪精度，在实际应用中一般并不单独使用。

发明内容

本发明的目的是为了解决上述问题，提出了一种基于改进自适应算法的科里奥利流量计信号跟踪方法，能够实时、精确跟踪和测量传感器输出信号频率。

本发明的技术方案：

一种基于改进自适应算法的科里奥利流量计信号跟踪方法，其特殊之处在于，包括以下步骤：

1】获得数字信号Y(n)；

2】对数字信号Y(n)进行信号预处理，得到数字信号X(n)，X(n)＝X(1)，X(2)…X(n)；信号预处理优选奇异谱分析的方法，具体包括步骤2.1】和2.2】；

2.1】对数字信号Y(n)进行分解：

构建行数为L，列数为K的矩阵Y₁。将Y(n)中的第i个元素Y(i)，第i+1个元素Y(i+1)，到第L+i-1个元素Y(L+i-1)依次赋值给Y₁的第i列所对应的L个元素，i的取值范围从1到K。

然后对矩阵Y₁进行奇异值分解(SVD)，得到特征矩阵U、V。

2.2】对数字信号Y(n)进行重构：

选取U、V中能够反映原数字信号中待测量信号特征的m个特征向量组生成重构矩阵RC；通过重构矩阵RC生成包含Y(n)中待测量信号特征的数字信号X(n)；

3】对数字信号X(n)进行滤波：

3.1】构建滤波器：

构建以

$H\left(z\right)=\frac{1+2k_0{z}^{-1}+z^{-2}}{1+2k_1{\mathrm{\ρz}}^{-1}+{\mathrm{\ρ}}^2{z}^{-2}}$

为传递函数的滤波器，给k₀、k₁、ρ赋初值，使输入滤波器的信号的频率落在滤波器的陷阱带宽内；

3.2】自适应调整滤波器模型，解算得到跟踪信号：

滤波器零点到原点与X轴的夹角α为cos^-1(-k₀)，滤波器极点到原点与X轴的夹角β为cos^-1(-k₁)，也就是说滤波器系数-k₀为滤波器零点到原点与X轴的夹角α的余弦，滤波器系数-k₁为滤波器极点到原点与X轴的夹角β的余弦；每次进行迭代运算时，固定滤波器零点，使滤波器极点以变角度的方式趋近于滤波器零点：

当输入第一个信号X(1)，求得滤波器的系数k₁，k₀＝k₁，根据公式f＝arccos(-k₀)解算得到第一个跟踪信号f₁，

当输入第二个信号X(2)，求得滤波器的系数k₁ ⁽¹⁾，当k₁ ⁽¹⁾>k₁，调整k₀使k₀<k₁ ⁽¹⁾，当k₁ ⁽¹⁾<k₁，调整k₀使得k₀>k₁ ⁽¹⁾，根据公式f＝arccos(-k₀)解算得到第二个跟踪信号f₂，

当输入第三个信号X(3)，求得滤波器的系数k₁ ⁽²⁾，当k₁ ⁽²⁾>k₁ ⁽¹⁾，调整k₀使k₀<k₁ ⁽²⁾，当k₁ ⁽²⁾<k₁ ⁽¹⁾，调整k₀使得k₀>k₁ ⁽²⁾，根据公式f＝arccos(-k₀)解算得到第三个跟踪信号f₃，

……，

当输入第n个信号X(n)时，求得滤波器的系数k₁ ^(n-1)，当k₁ ^(n-1)>k₁ ^{(n- 2)}，调整k₀使k₀<k₁ ^(n-1)；当k₁ ^(n-1)<k₁ ^(n-2)，调整k₀使k₀>k₁ ^(n-1)，根据公式f＝arccos(-k₀)解算得到第n个跟踪信号f_n，……；

4】曲线拟合：将步骤3.2】中一定时间内解算的所有跟踪信号进行曲线拟合获得最终跟踪信号；

5】依次输出或显示最终跟踪信号。

上述步骤3.2】调整k₀时，根据

γ(n)＝γ₀γ(n-1)+(1-γ₀)γ_∞,k₀＝k₁[1-γ(n)tan^-1(Δf)]进行调整。

步骤4】中曲线拟合的方法优先选择线性最小二乘拟合的方法。

步骤1】的数字信号Y(n)由以下方法得到：从科里奥利流量计的振动管传感器上获得原始频率信号，对原始频率信号进行A/D采样，获得数字信号Y(n)。

本发明带来的有益效果是：

本发明提出了一种采用基于奇异谱分析方法的科里奥利流量计信号预处理方法，并通过建立预估计信号频率变化趋势与自适应滤波器模型之间的联系，提出了一种基于优化零极点配置滤波算法的信号跟踪方法。该方法能够在自适应滤波的过程中根据频率的变化趋势动态地调节陷波器的模型，从而使陷波器更好的跟踪传感器输出信号的频率。理论分析和实际现场测试都表明该方法与传统方法相比在频率跟踪方面具有更好的收敛性和实时性。实现了在工业现场中对科里奥利流量计传感器信号进行长时间、高精度的跟踪。

附图说明

图1为α<β时，陷波器的幅频特性；

图2为α>β时，陷波器的幅频特性；

图3为相同信噪比条件下实际信号与仿真信号频谱的比较；

图4为基于格型自适应算法的信号跟踪方法与本发明信号跟踪方法的频率解算结果对比；

图5为基于格型自适应算法的信号跟踪方法与本发明信号跟踪方法的频率跟踪结果对比。

具体实施方式

1.基于奇异谱分析的信号重构方法：

奇异谱方法(SSA)是一种适用于研究周期振荡行为的分析方法,它是一种从时间序列的动力重构出发,并与经验正交函数(EOF)相联系的统计技术,是EOF分解的一种特殊应用。

首先，构建行数为L，列数为K的矩阵Y₁。将Y(n)中的第i个元素Y(i)，第i+1个元素Y(i+1)，到第L+i-1个元素Y(L+i-1)赋值给Y₁的第i列所对应的L个元素，i的取值范围从1到K。

然后对矩阵Y₁进行奇异值分解(SVD)：构建矩阵S＝XX^T，求得S的特征向量矩阵U。构建矩阵V＝X^TU。选取U、V中能够反映原数字信号中待测量信号特征的m个特征向量组U_m、V_m生成重构矩阵RC＝U_mV_m ^T。

最后通过重构矩阵RC生成包含Y(n)中待测量信号特征的数字信号X(n)。

2.基于优化零极点配置的自适应IIR滤波器模型的建立：

格型IIR陷波器是由Cho等人提出的，它由两个格型滤波器级联而成，是一种极联的全极点和全零点格型陷波器，其传递函数如下：

式(1)中的分子部分为整个格型滤波器提供一对共轭零点相当于MA模型；分母部分为整个格型滤波器提供一对共轭极点：相当于AR模型。对于二阶IIR陷波器，设零点到原点与X轴的夹角为α，极点到原点与X轴的夹角为β。当α<β时，陷波器的幅频特性(采样频率f_s＝5kHz)如图1所示，此时在陷波频率附近低频信号的衰减要小于高频信号的衰减；当α>β时，陷波器的幅频特性(采样频率f_s＝5kHz)如图2所示，此时在陷波频率附近低频信号的衰减要大于高频信号的衰减。

而式(1)中零点到原点与X轴的夹角α为：cos^-1(-k₀)，极点到原点与X轴的夹角β为：又因为格型IIR陷波器的参数ρ趋近于1但小于1，所以α始终小于β，这就限制了陷波器对于变化频率的跟踪能力。设陷波器零点到原点与X正半轴的夹角为α，极点到原点与X正半轴的夹角为β。如果能够在频率跟踪的过程中动态调节α角与β角的关系，则可以使陷波器更好的跟踪传感器输出信号的频率。因此，构建基于优化零极点配置的自适应IIR滤波器模型如下：

式(2)中零点到原点与X轴的夹角α为：cos^-1(-k₀)，极点到原点与X轴的夹角β为：cos^-1(-k₁)。因为通常信号的频率是未知的，所以需要对k₀进行估计。由陷波器的原理可知，陷波器的输出误差为：ε（n)＝e(n)-δ(n)，根据递推预测误差理论，取代价函数则k₀的估计可表示为：在跟踪传感器输出信号频率的过程中，如果检测到频率呈上升趋势，则调整k₀使陷波频率附近对于高频的衰减变小，而对于低频的衰减增大；反之，如果检测到频率呈下降趋势，则调整k₀使陷波频率附近对于低频的衰减变小，而对于高频的衰减增大。基于此，得到如下自适应更新公式：

ρ(n)＝ρ_nρ(n-1)+(1-ρ_n)ρ_∞

λ(n)＝λ_nλ(n-1)+(1-λ₀)λ_∞

Ψ(n)＝-X(n-1)+ρ(n)e(n-1)

P(n)＝P(n-1)/[λ(v)+Ψ(n)²P(n-1)]

W(n+1)＝W(n)+2μe(n)d(n)，

根据上述公式求得k₁后，本发明按照公式γ(n)＝γ₀γ(n-1)+(1-γ₀)γ_∞，k₀＝k₁[1-γ(n)tan^-1(Δf)]对k₀进行调整，体现了自适应变角度的过程，能够实时、精确跟踪信号频率。

式中，ρ(n)为去偏置参数；λ(n)为遗忘因子；γ(n)为校准系数；Ψ(n)为梯度参数，其中X(n)＝[X(n-10，...，X(n-M)]^T为滤波器输入；P(n)为协方差参数；e(n)为误差信号，其中d(n)为所期望的响应；W＝[W₁，W₂，W₃，...W_m]为权向量，μ为收敛参数；Δf为自适应迭代后估算频率的变化量。当Δf→0时，Δf与tan^-1(Δf)呈线性变化关系；当Δf→∞时，从而限制自适应迭代的最大步长以保证陷波器的稳定性。经过每次迭代使得k₁逐渐趋于k₀，当估计频率趋于稳定时，即Δf→0时，α→β，k₀→k₁，此时滤波器性能达到最优的状态。

另外，为保证陷波器的稳定性和平滑性，还需加入如下判断环节和平滑过程：

$k_0\left(n\right)=\left\{\begin{array}{cc}-1& k_0\left(n\right)\&le;-1\\ k_0\left(n\right)& -1\&le;k_0\left(n\right)\&le;1\\ 1& k_0\left(n\right)\&GreaterEqual;1\end{array}\right.$

k₀(n)＝0.5k₀(n)+0.5k₀(n-1)

3.上述方法跟踪科里奥利流量计的频率信号，具体步骤如下：

1】科里奥利流量计内的传感器敏感流体的振动频率信号获得原始信号，对原始信号进行A/D采样，获得采样频率为f_s的数字信号Y(n)。

2】对数字信号Y(n)进行基于奇异谱分析(SSA)方法的信号预处理：

2.1】对数字信号Y(n)进行分解：

构建行数为L，列数为K的矩阵Y₁。将Y(n)中的第i个元素Y(i)，第i+1个元素Y(i+1)，到第L+i-1个元素Y(L+i-1)赋值给Y₁的第i列所对应的L个元素，i的取值范围从1到K。

然后对矩阵Y₁进行奇异值分解(SVD)，得到矩阵U、V。

2.2】对数字信号Y(n)进行重构：

选取U、V中能够反映原数字信号中包含待测量信号特征的m个特征向量组生成重构矩阵RC。通过重构矩阵RC生成包含待测量信号特征的数字信号X(n)。

3】对数字信号X(n)进行滤波：

3.1】构建滤波器：

构建以

$H\left(z\right)=\frac{1+2k_0{z}^{-1}+z^{-2}}{1+2k_1{\mathrm{\&rho;z}}^{-1}+{\mathrm{\&rho;}}^2{z}^{-2}}$

为传递函数的滤波器，给k₀、k₁、ρ赋初值，使输入滤波器的信号的频率落在滤波器的陷阱带宽内；

3.2】自适应调整滤波器模型：

每次进行迭代运算时，固定滤波器零点，使滤波器极点以变角度的方式趋近于滤波器零点：

滤波器零点到原点与X轴的夹角α为cos^-1(-k₀)，滤波器极点到原点与X轴的夹角β为cos^-1(-k₁)，也就是说滤波器系数-k₀为滤波器零点到原点与X轴的夹角α的余弦，滤波器系数-k₁为滤波器极点到原点与X轴的夹角β的余弦；设输入信号X(n)＝X(1)，X(2)…X(n)

当输入第一个信号X(1)，求得滤波器的系数k₁，k₀＝k₁，根据公式f＝arccos(-k₀)解算得到第一个跟踪信号f₁，

当输入第二个信号X(2)，求得滤波器的系数k₁ ⁽¹⁾，当k₁ ⁽¹⁾>k₁，调整k₀使k₀<k₁ ⁽¹⁾，当k₁ ⁽¹⁾<k₁，调整k₀使得k₀>k₁ ⁽¹⁾，根据公式f＝arccos(-k₀)解算得到第二个跟踪信号f₂，

当输入第三个信号X(3)，求得滤波器的系数k₁ ⁽²⁾，当k₁ ⁽²⁾>k₁ ⁽¹⁾，调整k₀使k₀<k₁ ⁽²⁾，当k₁ ⁽²⁾<k₁ ⁽¹⁾，调整k₀使得k₀>k₁ ⁽²⁾，根据公式f＝arccos(-k₀)解算得到第三个跟踪信号f₃，

……，

当输入第n个信号X(n)时，求得滤波器的系数k₁ ^(n-1)，当k₁ ^(n-1)>k₁ ^{(n- 2)}，调整k₀使k₀<k₁ ^(n-1)；当k₁ ^(n-1)<k₁ ^(n-2)，调整k₀使k₀>k₁ ^(n-1)，根据公式f＝arccos(-k₀)解算得到第n个跟踪信号f_n，……，

4】曲线拟合：将一定时间内解算的所有跟踪信号进行曲线拟合获得最终跟踪信号；优先选择线性最小二乘拟合方法，最小二乘拟合方法较一般处理中对多个频率求平均从而获得最终频率的方法更能趋近于实际频率的真实值。

5】依次输出或显示最终跟踪信号。

步骤3.2】调整k₀时，根据

γ(n)＝γ₀γ(n-1)+(1-γ₀)γ_∞,k₀＝k₁[1-γ(n)tan^-1(Δf)]进行调整。

4.实际现场试验的结果与分析

在一般计算机仿真过程中使用的仿真信号为正弦波加高斯白噪声。然而在相同信噪比条件下，对比实际流量计传感器输出信号的频谱和仿真信号的频谱(如图3所示)可以发现，实际信号在中心频率附近的噪声要高于仿真信号的噪声，而仿真信号在远离中心频率处的噪声要高于实际信号的噪声，所以如果使用正弦波加高斯白噪声作为仿真信号进行试验是不能完全反映实际情况的。为了能够准确检验本发明方法在频率跟踪和解算方面的收敛性和快速跟踪特性，本发明通过对实际流量计传感器输出信号的频率进行跟踪和解算，从而验证本发明方法的可行性。

为验证本发明方法解算频率的收敛性，设计实验如下：通过泵将液体送入流量计中，待检测到通过流量计的瞬时质量流量趋于稳定后，通过DSP对A/D采样的信号进行频率解算，检测算法解算频率的收敛范围，即可验证算法解算频率的收敛性。设置A/D采样频率f_s为5kHz，每0.1秒输出一个解算频率，测试10秒，测试结果如图4所示。从图4中可以看出，在传感器输出信号频率趋于稳定后，应用本发明方法解算频率的收敛范围为0.0005Hz，满足实际应用要求。而应用格型自适应方法解算频率的收敛范围为0.0015Hz。

表1给出了图4中应用两种不同方法解算频率结果的均方误差，从中可以看出本发明在频率跟踪方面的收敛性要明显优于基于格型自适应算法的信号跟踪方法。

表1

方法	频率的均方误差
		基于格型自适应算法的信号跟踪方法	1.7619×10-7
本发明的信号跟踪方法	2.1452×10-8

为验证本发明在频率变化时的快速跟踪特性，设计实验如下：通过泵将水送入流量计中，待检测到通过流量计的瞬时质量流量趋于稳定后，向流量计中送入水和酒精的混合液，使通过流量计的液体的密度迅速发生变化，则其间传感器输出信号的频率也将迅速发生变化。通过DSP对A/D采样的信号进行频率解算，根据两种算法解算变化频率的结果是否分别进入其收敛范围，判定算法跟踪上新频率所用的时间，即可验证算法在频率变化时的快速跟踪特性。设置A/D采样频率f_s为5kHz，每0.1秒输出一个解算频率，测试10秒，测试结果如图5所示。

从图5中可以看出，当传感器输出频率发生突变后，本发明跟踪上了新的频率(频率解算结果进入0.0005Hz的收敛范围)所用时间明显快于基于格型自适应算法的信号跟踪方法跟踪上新的频率(频率解算结果进入0.0015Hz的收敛范围)所用时间。

综上所述，本发明与基于格型自适应算法的信号跟踪方法相比不仅在信号频率的解算方面具有更好的收敛性，而且在信号频率的跟踪上具有更好的实时性。

Claims (5)

1.一种基于改进自适应算法的科里奥利流量计信号跟踪方法，其特征在于，包括以下步骤：

1】获得数字信号Y(n)；

2】对数字信号Y(n)进行信号预处理后获得数字信号X(n)，X(n)＝X(1)，X(2)…X(n)；

3】对数字信号X(n)进行滤波：

3.1】构建滤波器：

构建以

$H\left(z\right)=\frac{1+2k_0{z}^{-1}+z^{-2}}{1+2k_1{\mathrm{\&rho;z}}^{-1}+{\mathrm{\&rho;}}^2{z}^{-2}}$

为传递函数的滤波器，给k₀、k₁、ρ赋初值，使输入滤波器的信号的频率落在滤波器的陷阱带宽内；

3.2】自适应调整滤波器模型，解算得到跟踪信号：

滤波器零点到原点与X轴的夹角α为cos^-1(-k₀)，滤波器极点到原点与X轴的夹角β为cos^-1(-k₁)，也就是说滤波器系数-k₀为滤波器零点到原点与X轴的夹角α的余弦，滤波器系数-k₁为滤波器极点到原点与X轴的夹角β的余弦；每次进行迭代运算时，固定滤波器零点，使滤波器极点以变角度的方式趋近于滤波器零点：

当输入第一个信号X(1)，求得滤波器的系数k₁，k₀＝k₁，根据公式f＝arccos(-k₀)解算得到第一个跟踪信号f₁，

当输入第二个信号X(2)，求得滤波器的系数k₁ ⁽¹⁾，当k₁ ⁽¹⁾>k₁，调整k₀使k₀<k₁ ⁽¹⁾，当k₁ ⁽¹⁾<k₁，调整k₀使得k₀>k₁ ⁽¹⁾，根据公式f＝arccos(-k₀)解算得到第二个跟踪信号f₂，

当输入第三个信号X(3)，求得滤波器的系数k₁ ⁽²⁾，当k₁ ⁽²⁾>k₁ ⁽¹⁾，调整k₀使k₀<k₁ ⁽²⁾，当k₁ ⁽²⁾<k₁ ⁽¹⁾，调整k₀使得k₀>k₁ ⁽²⁾，根据公式f＝arccos(-k₀)解算得到第三个跟踪信号f₃，

……，

当输入第n个信号X(n)时，求得滤波器的系数k₁ ^(n-1)，当k₁ ^(n-1)>k₁ ^{(n- 2)}，调整k₀使k₀<k₁ ^(n-1)；当k₁ ^(n-1)<k₁ ^(n-2)，调整k₀使k₀>k₁ ^(n-1)，根据公式f＝arccos(-k₀)解算得到第n个跟踪信号f_n，……；

4】曲线拟合：将步骤3.2】中一定时间内解算的所有跟踪信号进行曲线拟合获得最终跟踪信号；

5】依次输出或显示最终跟踪信号。

2.根据权利要求1所述的基于改进自适应算法的科里奥利流量计信号跟踪方法，其特征在于：所述信号预处理方法是奇异谱分析的方法，具体步骤是：

2.1】对数字信号Y(n)进行分解：

构建行数为L，列数为K的矩阵Y₁：将Y(n)中的第i个元素Y(i)，第i+1个元素Y(i+1)，到第L+i-1个元素Y(L+i-1)依次赋值给Y₁的第i列所对应的L个元素，i的取值范围从1到K；

然后对矩阵Y₁进行奇异值分解，得到特征矩阵U、V；

2.2】对数字信号Y(n)进行重构：

选取U、V中能够反映原数字信号Y(n)中待测量信号特征的m个特征向量组生成重构矩阵RC，通过重构矩阵RC生成包含Y(n)中待测量信号特征的数字信号X(n)。

3.根据权利要求1所述的基于改进自适应算法的科里奥利流量计信号跟踪方法，其特征在于，步骤3.2】调整k₀时，根据γ(n)＝γ₀γ(n-1)+(1-γ₀)γ_∞，k₀＝k₁[1-γ(n)tan^-1(Δf)]进行调整。

4.根据权利要求1所述的基于改进自适应算法的科里奥利流量计信号跟踪方法，其特征在于，所述曲线拟合的方法是线性最小二乘拟合。

5.根据权利要求1所述的基于改进自适应算法的科里奥利流量计信号跟踪方法，其特征在于，步骤1】的数字信号Y(n)由以下方法得到：从科里奥利流量计的振动管传感器上获得原始频率信号，对原始频率信号进行A/D采样，获得数字信号Y(n)。