摘 要：提出用频率、幅值和相位均按照随机游动模型变化的信号来描述科氏流量传感器的输出信号。采用具有跟踪频率变化能力的自适应陷波器，对信号进行滤波，以求其频率；采用自适应谱线增强器从含有噪声的数据中提取出所需要的信号；然后采用滑动Goertzel算法计算两路信号之间的相位差，并通过频率和相位差计算出时间差，求得质量流量。仿真结果表明所研究的方法是有效的。

关键字：科里奥利质量流量计 随机游动模型 自适应陷波器 自适应谱线增强 滑动Goertzel算法

    1 引言

    科里奥利质量流量计(以下简称为科氏流量计)，是一种可以直接实现质量流量测量并可同时获取流体密度值的流量计。它是当前发展最为迅速的流量计之一，具有广阔的应用前景。科氏流量计要求其信号处理电路精确地测量出来自两个流量传感器信号的相位差。目前广泛使用的国内外产品绝大部分采用基于模拟和数字电路的信号处理方法，如：放大、滤波、整形和计数等，但这种处理方法往往存在很多局限^[1]。随着数字信号处理技术和DSP芯片的发展，国内外的相关研究机构做了大量研究工作，希望采用数字信号处理技术来解决这些问题^[1~8]。但是这些研究基本上都是以频率、幅值和相位不随时间变化的信号模型作为研究对象，而实际的科氏流量计信号由于受到管内流体的影响，这些特征是随时间变化的，因此并不能较为真实地反映实际的情况。

    为了更为真实地描述科氏流量计的信号的实际特性，文中以频率、幅值和相位均按照随机游动模型(Random Walk Mode1)变化的信号模拟科氏流量计信号，并研究一种新的信号处理方法，从含有噪声干扰的信号中准确测得质量流量。首先采用基于ⅡR型自适应陷波器的自适应谱线增强器(ALE)将信号中噪声滤除，同时测得信号的频率，并跟踪其变化。然后采用滑动Goertzel算法，得到增强后的信号的相位差和时间差，从而测得质量流量。值得注意的是，由于许多实际信号都可以用该模型模拟，所以文中所介绍方法的应用领域并不仅局限于对科氏流量计信号的处理。

    2 信号模型

    根据科氏流量计的工作原理，其传感器的输出信号为正弦信号。由于受到管内流体特性，如流速，密度和流体脉动等因素变化的影响^[9~10]，信号的幅值、频率和相位也会相应地发生变化。这种变化通常是不可预见和随机的。同时由于传感器测量噪声和环境噪声的存在，所以定义如下模型来模拟真实的科氏流量计信号：

    y(n)=A(n)sin(ω(n)+φ(n))+ σ_ee(n)     (1)

    其中：e(n)是零均值，方差为1的白噪声。幅值A(n)、归一化频率ω(n)和相位φ(n)按照随机游动模型变化，即：

    A(n)=A(n-1)+σ_Ae_A(n)       (2)
    ω(n)=ω(n-1)+σ_ωe_ω(n)     (3)
    φ(n)=φ(n-1)+σ_φe_φ(n)    (4)

    其中：e_A(n)、e_ω(n)和e_φ(n)均为零均值，方差为l的白噪声，且e(n)、e_A(n)、e_ω(n)和e_φ(n)互不相关。需要指出的是，由于幅值、频率和相位应该随时问缓慢变化，因此在计算机仿真中σ_A、σ_ω和σ_Ф的取值与f_s有关系，f_s越大，σ_A、σ_ω和σ_φ应相应地取小一些，反之亦然。

    3 自适应ⅡR陷波器和谱线增强

    3.1 原理

    研究的是零极点约束的ⅡR陷波滤波器^[11~12]。零极点约束是指陷波器的零极点满足如下条件：对于每对零极点，零点在单位圆上，并且位于陷波频率处。极点则在单位圆内，且与零点同一角度，并尽可能的靠近零点，极点与原点的距离为ρ。由此可得陷波器的传递函数为：

        (5)

    其中：α_k=-2cosω_k，ω_k为陷波频率，n为陷井数。

    对于科氏流量计信号，n=1，则(5)式简化为：

        (6)

    假设信号为y(n)=Asin(ω₀n+φ)+e(n)，

    当α=-2cosω₀时，

    (1+αz^-1+z^-2)e(n)=(1+αz^-1+z^-2)y(n)      (7)

    成立，其中：z^-1为单位延迟因子，即：y(n)z^-1=y(n-1)。式(7)表明当陷波频率等于信号频率且ρ=1时，陷波器的输出中只含有白噪声e(n)，正弦信号被完全滤掉，陷波器相当于y(n)的一个白化滤波器。对于0<ρ<1，陷波器的输出为：

        (8)

    当通常信号的频率是未知的，因此需要对α进行估计。陷波器的输出误差为，根据递推预测误差理论^[16]，取代价函数

        (9)

    α的估计可表示为：

        (10)

    由于ρ趋近于1，，因此可由式(11)递推计算：

         (11)

    其中：

        (12)
        (13)

    P(n)可由式(14)递推计算：

        (14)

    λ(n)为遗忘因子：

    λ(n)=λ₀λ(n-1)+(1-λ₀)λ_∞     (15)

    在以上讨论中，ρ一直假设为定值。但在实际算法中，由于ρ决定每个陷井的带宽，在输入信号的先验知识未知的情况下，如果ρ非常趋近于1，即极点靠近零点，则陷井可能因为过窄而无法落在信号频率处，也就无法感知信号的存在。因此在陷波的开始阶段ρ应取得稍小一点，使算法收敛到期望的传递函数上，然后再增大ρ。为此ρ可改写为ρ(n)，

    ρ(n)=ρ₀ρ(n-1)+ (1-ρ₀)ρ_∞     (16)

    3.2 自适应算法的改进

    在文献[1]中使用

    λ(O)=0.95、λ₀=0.99、λ_∞=1
    ρ(0)=0.7、ρ_U=0.99、ρ_∞=0.995

    对时不变信号进行处理，取得很好的效果。但是将其应用到时变信号模型后，发现自适应算法无法对信号频率的变化进行跟踪。为此使用两种方法^[13~15]对算法进行改进，使其能够响应频率的缓慢变化。

    算法1 如前所述，ρ(n)决定陷波器的陷井带宽，在对频率固定的信号进行处理时，ρ应选的尽可能靠近1，以获得良好的滤波效果。但在频率时变的情况下，ρ如果过于靠近1，将使得陷波器在收敛后无法感应到频率的变化，从而丧失了跟踪的能力，因此需要适当地调小ρ。同时，遗忘因子λ=1意味着在T=1/(1-λ)个采样点内，信号频率可近似认为是常量。但在频率时变的情况下，该条件已不再满足，所以λ必须小于1。减小ρ和λ可以使陷波器具有了跟踪频率变化的能力，但却降低了陷波器的滤除噪声的能力。因此，ρ和λ的取值应综合考虑这两方面因素，文献[13~14]对此进行了详细的讨论，文中通过调整ρ_∞和λ_∞使ρ和λ满足要求。

    算法2 在算法1中只考虑了通过自适应算法调整α使得代价函数F(α)最小，而ρ的更新算法与F无关。为了使ρ也以自适应的方式更新^[15]，在本算法中将代价函数F(α)改写成F(α，ρ)，并使

        (17)

    与类似， 也可由以下递推公式进行计算：

        (18)

    其中：

          

    P_p(n)可由式(21)递推计算：

       (21)

    其中：λ_ρ为计算时的遗忘因子，应满足λ_ρ>λ以使ρ的自适应过程比信号频率即α的自适应过程稍慢。λ_ρ的取值与信号特征有关，需要在仿真中确定。由于λ的最优值可以认为与ρ的相等，所以此时λ可以按式(22)递推计算^[15]：

    λ(n)=0.995λ(n-1)+0.005ρ(n)     (22)

    4 相位差及时间差的计算

    采用傅立叶分析技术计算科氏流量计振动管信号的相位差。为了能够快速地计算出两信号的傅立叶系数，先使用矩形窗对经过自适应谱线增强的两路信号进行截取，然后用滑动Goertzel算法^[17~18]计算傅立叶系数并计算出相位差。

    4.1 滑动Goertzel算法

    常规Goertzel算法的传递函数如下^[17]：

        (23)

    其频率为等间隔分布，当信号的频率恰好对应于某个k值时，可以获得精确的结果。而当信号频率落在频率间隔内时，就会由于泄漏问题而产生较大误差。泄漏问题对于计算相位的影响要远远大于对功率谱的计算[7]。

    滑动Goertzel算法^[18]克服了常规Goertzel算法的这一缺点，它用实际的信号频率替换传递函数中的2πk/N，因此在信号频率不变的情况下，经过一定时间的收敛，可以精确地计算出信号的傅立叶系数，并且可以在每个采样点计算一次傅立叶系数。正是这一特性，使得滑动Goertzel算法在信号的幅值和相位随时间变化时也能够跟踪其变化。滑动Goertzel算法如图1所示。

    设固定频率信号y(n)的傅立叶展开为：y(n)=αsin(nω)+bcos(nω)，其变换为：

        (24)

    则：

        (25)

    V(z^-1)是中间变量v(n)的z变换，对其进行z反变换^[18]，得：

    

    其中：

    

    从式(26)~(29)可以看出在年n的值比较大时Δ₁和Δ₂可以忽略。经过计算可以得到：

    

    4.2 相位差及时间差的计算

    通过仿真发现当信号的频率随时间变化时滑动Goertzel算法并不能对频率的变化实施长时间的跟踪，而且由于陷波器存在一定的频率跟踪误差，使得信号相位差的计算误差较大。为此采用矩形窗对经过谱线增强信号进行截取，每个窗包含300个采样点。矩形窗宽度的选择是滑动Goertzel算法收敛时间与相位差计算精度之间的折衷，窗口过短可能使得滑动Goertzel算法并没有完全收敛，窗口过长则使相位计算精度下降。

    文献[18]中的滑动Goertzel算法只考虑了信号幅值和相位的变化，并没有考虑信号频率的变化。为了将该方法应用于文中所讨论的科氏流量计信号，将式(30)~(31)改写为：

   

    则相位差的计算公式为：

        (36)

    其中Ø(n)和Ø(n)分别为两路信号的相位估计。时间差的计算公式为：

    

    其中为信号频率的估计，f_s为采样频率。

    5 仿真结果

    为了检验上述方法的效果，使用Matlab®对算法进行了仿真。由于科氏流量计的信号频率通常在100±4Hz范围内变化，相位在±4°以内，因此仿真中选择10000个采样点，且：

    

    σ_e=0.2、σ_A=10^-3、σ_ω=10^-4、σ_φ=10^-5为了检验滑动Goertzel算法计算相位差的效果，在保持e、e_A、e_ω、e_φ以及e_Ø不变(这里的不变并不意味着这些量为常量，它们依旧是互不相关的白噪声)的情况下分别对φ(0)=π/180、Ø(0)=0.5π/18O和Ø(0)=0.1π/180进行仿真。为节约篇幅，在5.1节仅给出Ø(0)=0.5π/180时的时间差收敛波形。

    5.1 自适应陷波器的仿真结果

    对于第一种自适应算法，一些参数的初始值选择如下：

    
    λ(0)=0.9、λ0=0.9、λ∞=0.95
    ρ(0)=0.8、ρ0=0.8、ρ∞=0.98

    仿真结果如图2~图3所示。

    第二种自适应算法的仿真中，选择

     
   λ(0)=0.85、ρ(0)=0.8、λρ=0.999

    仿真结果如图4~图6所示。

    从图2和图4可以看出，算法1相对于算法2具有较快的收敛速度，而从图3和图5中又可以看出在跟踪频率变化的能力上算法2优于算法1。

    表1列出了实际信号频率的均值与其估计值的均值按下式计算：

        (38)

    以去除陷波器收敛过程的影响。和分别对应算法

表1 的比较

ω(n)	(n)	(n)
0.31444	0.31448	0.31449

    表2列出了两种算法所得频率估计的方差σ²和均方误差(MSE)，其中MSE的计算公式如下：

        (39)

表2 的方差和MSE

	σ2	MSE
(n)	4.29×10-6	1.96×10-6
(n)	2.73×10-6	8.17×10-7

    从表1和表2中可以看出虽然算法1的均值比算法2的稍微接近真实值，但是算法2在方差和均方误差方面更具优势，尤其是算法2的MSE只相当于算法1的一半，说明算法2跟踪频率变化的能力比算法1强。

    5.2 滑动Goertzel算法的仿真结果

    选择一个矩形窗(第9601~9900采样点)来说明用滑动Goertzel算法计算相位差的结果。图7和图8分别为φ(0)=0.5π/18O时的两种算法下的时间差收敛波形。可以清楚地看到在相同条件下，由于算法2的频率跟踪误筹较小，时间差的计算精度也较高。

    下面给出仿真的数值结果，表3为三种φ(O)下分别使用两种陷波器算法所得到的相位差的均值和均方误差MSE以及实际信号的，单位为rad。表4为时间差的均值和均方误差MSE以及实际信号的Δt，单位为ms。考虑到滑动Goertzel算法的收敛问题，这里的均值和MSE的计算只使用了第9700~9900采样点。

表3 相位差均值和MSE

		(算法1)	(算法2)	MSE(算法1)	MSE(算法2)
φ(0)=1.50π/180	0.04968	0.05049	0.04939	1.14×10-6	5.55×10-7
φ(0)=1.25π/180	0.04095	0.04165	0.04073	8.09×10-7	3.73×10-7
φ(0)=1.00π/180	0.03539	0.03480	0.03507	6.40×10-7	3.37×10-7
φ(0)=0.90π/180	0.02874	0.02927	0.02860	4.46×10-7	1.81×10-7
φ(0)=0.75π/180	0.02350	0.02397	0.02340	3.24×10-7	1.21×10-7
φ(0)=0.60π/180	0.01826	0.01866	0.01820	2.23×10-7	7.33×10-8
φ(0)=0.50π/180	0.01794	0.01756	0.01777	2.21×10-7	9.05×10-8
φ(0)=0.40π/180	0.01128	0.01159	0.01127	1.21×10-7	3.08×10-8
φ(0)=0.25π/180	0.00605	0.00629	0.00607	6.76×10-8	1.44×10-8
φ(0)=0.10π/180	0.00397	0.00377	0.00393	5.28×10-8	1.13×10-8

表4 时间差均值和MSE

		(算法1)	Δt(算法2)	MSE(算法1)	MSE(算法2)
φ(0)=1.50π/180	0.07887	0.08093	0.07908	5.97×10-6	1.21×10-6
φ(0)=1.25π/180	0.06502	0.06676	0.06521	4.21×10-6	8.27×10-7
φ(0)=1.00π/180	0.05623	0.05558	0.05589	1.08×10-6	6.97×10-7
φ(0)=0.90π/180	0.04562	0.04692	0.04578	2.26×10-6	4.18×10-7
φ(0)=0.75π/180	0.03731	0.03842	0.03746	1.61×10-6	2.88×10-7
φ(0)=0.60π/180	0.02900	0.02992	0.02914	1.08×10-6	1.84×10-7
φ(0)=0.50π/180	0.02850	0.02805	0.02832	3.80×10-7	1.86×10-7
φ(0)=0.40π/180	0.01791	0.01858	0.01804	5.36×10-7	8.67×10-8
φ(0)=0.25π/180	0.00960	0.01008	0.00971	2.60×10-7	4.51×10-8
φ(0)=0.10π/180	0.00632	0.00603	0.00626	1.10×10-7	2.51×10-8

    由此可以看出，与算法1相比，算法2除了收敛时稍长这个缺点外，其余方面均比算法1优越。

    6 结束语

    根据科氏流量计的实际工作原理，以频率、幅值和相位均按照随机游动模型变化的时变信号取代先前研究[1~8]中所普遍使用的频率、幅值和相位均固定的时不变信号来更为精确地模拟科氏流量计信号，并研究了具有跟踪频率变化能力的自适应陷波器，以及基于这种陷波器的自适应谱线增强器将混叠在信号中的噪声中滤除，同时对信号的时变频率进行估计。然后采用滑动Goertzel算法计算增强后的信号的傅立叶系数，以求得科氏流量计两路信号的相位差，并进一步地计算它们的时间差。从仿真的结果可以看出频率估计、相位差以及时间差计算都具有较高的精度，获得了令人满意的效果。

    参考文献

    1 徐科军，倪伟.一种科里奥利质量流量计的信号处理方法.计量学报，2001，22(4)：254~257.
    2 Yoshimura Hiroyuki.Phase difference measuring appa-ratus and flowmeter thereof.European Patent Applica-tion。EP 0702212A2，20.03.1996.
    3 H.V.Derby。T.Bose.S.Rajan.Method and apparatus for adaptive line enhancement in Coriolis mass flow meter measurement. US Patent No.5555190，Sep.10，1996.
    4 Yoshimura Hiroyuki.Phase difference measuring apparatus for measuring phase difference between input signals. European Patent Application，EP 0791807A2，1997.
    5 徐科军，姜汉科，苏建徽，等.科氏流量计信号处理中频率跟踪方法的研究.计量学报，1999，19(4)：304~307.
    6 徐科军，于翠欣，姜汉科，等.Research on signal processing of eoriolis mass flowmeter.ICEM I’99，Harbin，China，August，18~21，1999，电子测量与仪器学报，1999，13(增刊)：835~841.
    7 倪伟.用于科里奥利质量流量计的相位差测量方法与装置[D].合肥：合肥工业大学自动化研究所，1999.
    8 于翠欣.科里奥利质量流量计数字信号处理系统的研制[D].合肥：合肥工业大学自动化研究所，2000.
    9 R.Cheesewright，C.Clark. The effect of flow pulsations on Coriolis mass flow meters.Journal of Fluidsand Structures，1998，12：1025~1039.
    10 R.Cheesewright，C.Clark，D.Bisset.Understanding the experimental response of Coriolis mass flow meters to flow pulsations.Flow M easurement and Instrumen-tation，1999，10(4)：207~215.
    11 A.Nehorai.A minimal parameter adaptive notch-filter with Constrained poles and zeros. IEEE Trans.Acoust，Speech Signal Processing，1985，ASSP-33：983~996.
    12 P.Stoica.A.Nehorai. Performance analysis of an adaptive notch filter with constrained poles and zeros.IEEE Trans. Acoust， Speech Signal Processing，1988，36：911~919.
    13 P.Hande1. A.Nehorai. Tracking analysis of a constrained adaptive notch filter. Proc.IEEE ICASSP 1991，Toronto，May 1991，3209~3212.
    14 P.Handel，A.Nehorai.Tracking analysis of an adaptive notch filter with constrained poles and zeros.IEEE Trans.Signal Processing，1994，42：281~291.
    15 M.Dragosevic.S.Stankovic. An adaptive notch filter with improved tracking properties.IEEE Trans.Signal Processing，1995，43：2068~2078.
    16 L.Ljung.System identification：theory for the user (Second Edition).Prentice Hall PTR，1999，北京：清华大学出版社，2002.
    17 S.K.Mitra. Digital signal processing：a computer-based approach (Second Edition).New York，NY ：McGraw-Hill，2001，北京：清华大学出版社，2001.
    18 J.Chieharo.M.Kilani.A sliding Goertzel algorithm.Signal Processing，1996，52：283~297.