摘 要：提出用頻率、幅值和相位均按照隨機游動模型變化的信號來描述科氏流量傳感器的輸出信號。采用具有跟蹤頻率變化能力的自適應陷波器，對信號進行濾波，以求其頻率；采用自適應譜線增強器從含有噪聲的數據中提取出所需要的信號；然后采用滑動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.