基于過零檢測(cè)的相位差計(jì)算方法的改進(jìn)

徐長(zhǎng)利　李波　鐘其水　鐘驍　曾學(xué)兵

摘要：針對(duì)科里奧利質(zhì)量流量計(jì)相位差計(jì)算中實(shí)時(shí)性較差和精度不高的問題，在過零檢測(cè)的基礎(chǔ)上對(duì)相位差計(jì)算方法進(jìn)行改進(jìn)。首先采用改進(jìn)的自適應(yīng)格型陷波器（Adaptive Lattice Notch Filter，ALNF）對(duì)信號(hào)進(jìn)行頻率跟蹤，通過快速傅里葉變換（Fast Fourier Transform，F(xiàn)FT）進(jìn)行頻率預(yù)估縮短ALNF收斂時(shí)間，引入反饋因子提高ALNF收斂精度;然后針對(duì)離散信號(hào)無法求解零點(diǎn)的問題，采用傅里葉級(jí)數(shù)擬合的方式求解信號(hào)的零點(diǎn);最后采用卡爾曼濾波器對(duì)誤差進(jìn)行校正，從而提高相位差的計(jì)算精度。

關(guān)鍵詞：自適應(yīng)格型濾波器;傅里葉級(jí)數(shù)擬合;過零檢測(cè);卡爾曼濾波器;相位差

中圖分類號(hào)：TH814 ? ?文獻(xiàn)標(biāo)志碼：A ? ?文章編號(hào)：1671-0797（2022）02-0035-04

DOI：10.19514/j.cnki.cn32-1628/tm.2022.02.010

0 ? ?引言

在許多工業(yè)領(lǐng)域當(dāng)中，流體的流量都有著十分重要的應(yīng)用，其中科里奧利質(zhì)量流量計(jì)（簡(jiǎn)稱“科氏力流量計(jì)”）的應(yīng)用尤為廣泛。該流量計(jì)主要是由一次儀表和二次儀表組成。

（1）一次儀表主要由測(cè)量管（U型管、直管等）、激振器和傳感器等組成。在進(jìn)行流量測(cè)量時(shí)，流體會(huì)流過測(cè)量管，在科里奧利力的作用下，測(cè)量管會(huì)發(fā)生扭曲，通過傳感器采集會(huì)得到兩路具有相位差的正弦信號(hào)，計(jì)算兩路信號(hào)之間的相位差便可進(jìn)一步轉(zhuǎn)化為流體的流量。

（2）二次儀表主要是數(shù)字信號(hào)處理部分，用來計(jì)算流量計(jì)的相位差。

迄今為止，有許多學(xué)者對(duì)科氏力流量計(jì)相位差計(jì)算方法進(jìn)行了深入而廣泛的研究，主要有復(fù)系數(shù)濾波算法[1-2]、計(jì)及負(fù)頻率的影響的DTFT算法（Sliding Discrete Time Fourier Transform，SDTFT）[3]、希爾伯特變換（Hilbert

Transform，HT）方法[4-5]以及基于相關(guān)法計(jì)算方法[6]等。其中，復(fù)系數(shù)濾波算法是采用復(fù)帶通濾波器（Complex Bandpass Filtering，CBF）、復(fù)陷波濾波器（Complex Notch Filtering，CNF）以及它們的組合來求解相位差，其計(jì)算量較大，收斂時(shí)間較長(zhǎng)。SDTFT雖然克服了希爾伯特變換方法易受噪聲影響的缺點(diǎn)以及相關(guān)法需要整周期采樣的缺點(diǎn)，但其為了防止數(shù)值溢出，采用滑動(dòng)窗口的方式減少了每次計(jì)算的采樣點(diǎn)數(shù)，導(dǎo)致精度有所缺失。

為解決上述問題，本文采用平穩(wěn)單相流信號(hào)模型作為科氏力流量計(jì)的信號(hào)模型。針對(duì)流量計(jì)信號(hào)處理中的實(shí)時(shí)性和精度問題，提出以過零檢測(cè)為基礎(chǔ)，利用改進(jìn)的自適應(yīng)格型陷波器進(jìn)行頻率跟蹤，在計(jì)算出結(jié)果之后，采用卡爾曼濾波器進(jìn)行誤差校正，最終實(shí)現(xiàn)相位差計(jì)算精度的提高。

1 ? ?相位差計(jì)算改進(jìn)方法與步驟

科氏力流量計(jì)在理想狀態(tài)下的兩路輸出信號(hào)均為頻率和幅值不變的正弦信號(hào)，它們之間的相位差反映流體的流量。

科氏力流量計(jì)在穩(wěn)定工作處狀態(tài)下，其信號(hào)主要受環(huán)境噪聲干擾，體現(xiàn)在信號(hào)的幅度發(fā)生變化，因此采用平穩(wěn)單相流信號(hào)進(jìn)行仿真[3]，信號(hào)模型如下：

y（n）=Asin（ωn+φ）+Be（n） ? ? （1）

式中：e（n）為均值為0、方差為1的高斯白噪聲;B主要用來控制噪聲的幅度，進(jìn)而控制信噪比;A為幅值，視為常量;ω為科氏力流量計(jì)的振動(dòng)頻率;φ為信號(hào)的原始相位。

相位差的計(jì)算流程如圖1所示，首先利用FIR低通濾波器進(jìn)行初步降噪。經(jīng)過降噪之后采用基于FFT和反饋因子改進(jìn)的自適應(yīng)格型陷波器進(jìn)行頻率跟蹤以及進(jìn)一步的降噪，跟蹤的頻率用于過零檢測(cè)中進(jìn)行時(shí)間差向相位差的轉(zhuǎn)化。

然后采用傅里葉級(jí)數(shù)擬合進(jìn)行零點(diǎn)的求解，將兩路信號(hào)零點(diǎn)之間的時(shí)間差ΔT轉(zhuǎn)換為相位差，相位差ΔP計(jì)算如下[7]：

ΔP=360°×ΔT×f ? ? （2）

式中：f為跟蹤后的頻率。

最后采用卡爾曼濾波器進(jìn)行相位差誤差的校正。

1.1 ? ?基于FFT頻率預(yù)估和反饋因子改進(jìn)的自適應(yīng)格型陷波器

ALNF是一種可以根據(jù)信號(hào)的特點(diǎn)進(jìn)行遞歸收斂并估計(jì)信號(hào)頻率的濾波器。其主要采用格型IIR陷波器實(shí)現(xiàn)，主要由全極點(diǎn)和全零點(diǎn)兩個(gè)格型濾波器級(jí)聯(lián)而成，其傳遞函數(shù)為：

H（z）= ? ? ? ? ? ?（3）

式中：k0為權(quán)系數(shù)，用于計(jì)算信號(hào)頻率;ρ為偏置參數(shù)，決定陷波帶寬。

此種形式是將零點(diǎn)固定在單位圓上，k0采用Burg算法進(jìn)行自適應(yīng)調(diào)整[3]。

ALNF雖然能夠跟蹤頻率，但是由于其算法存在迭代遞歸過程，在計(jì)算時(shí)收斂過程較長(zhǎng)，并且當(dāng)頻率發(fā)生改變時(shí)會(huì)浪費(fèi)較長(zhǎng)的時(shí)間進(jìn)行頻率的自適應(yīng)跟蹤。因此，需要采用FFT進(jìn)行信號(hào)頻率的預(yù)估，然后再進(jìn)行自適應(yīng)跟蹤濾波，這樣可以縮短收斂過程;同時(shí)引入反饋因子Fb，反饋因子采用輸入信號(hào)與輸出信號(hào)的互相關(guān)，選擇合適的閾值之后，如果互相關(guān)函數(shù)最大值小于該閾值，那么表示跟蹤效果差，這時(shí)需要立刻調(diào)整陷波的帶寬ρ，對(duì)頻率進(jìn)行跟蹤。

其改進(jìn)方式如圖2所示。

設(shè)輸入信號(hào)為s（n）=x（n）+e（n），其中x（n）為原始信號(hào)，e（n）為高斯白噪聲;（n+m）為跟蹤降噪之后的信號(hào)，其中m為信號(hào)延時(shí)。則反饋因子Fb為：

Fb=max[R（m）]=maxs（n）（n+m） ? ? ? ?（4）

當(dāng)輸出信號(hào)（n+m）與輸入信號(hào)s（n）相關(guān)性很低時(shí)，需要通過調(diào)整ρ來調(diào)節(jié)陷波器帶寬，進(jìn)而重新跟蹤頻率。

計(jì)算出信號(hào)之間的相關(guān)性之后，利用相關(guān)函數(shù)的最大值作為反饋因子Fb，并設(shè)定閾值Th=8，小于閾值Th則認(rèn)為信號(hào)之間相關(guān)性低，相關(guān)性低的情況則需要減小ρ來加快信號(hào)的收斂。

ALNF改進(jìn)前后對(duì)比如圖3所示。

由圖3對(duì)比分析可知，經(jīng)過FFT頻率預(yù)估和反饋因子改進(jìn)之后，信號(hào)的收斂區(qū)間減小，達(dá)到了改進(jìn)的目的，可以用來提高相位差計(jì)算的實(shí)時(shí)性。

1.2 ? ?基于傅里葉級(jí)數(shù)擬合的過零檢測(cè)

過零檢測(cè)的基本原理是通過提取周期性的零點(diǎn)來得到左右信號(hào)之間的時(shí)間差，進(jìn)而獲取相位差，其具有無收斂過程、運(yùn)算量小和響應(yīng)速度快的優(yōu)點(diǎn)，但是在零點(diǎn)的求解過程中極易受到噪聲的干擾。

在經(jīng)過兩次降噪之后，信號(hào)的信噪比有了較好的提升，因此可以采用過零檢測(cè)的方式實(shí)現(xiàn)信號(hào)相位差的計(jì)算。

傅里葉級(jí)數(shù)擬合是以三角函數(shù)表示的無窮級(jí)數(shù)，比較適合用來擬合周期性函數(shù)。

在信號(hào)滿足狄利克雷條件的情況下，函數(shù)f（t）可以表示為：

f（t）=a0+[akcos（kωt）+bksin（kωt）] （5）

式中：a0、ak、bk為傅里葉系數(shù)。

在實(shí)際應(yīng)用中，無法采用無限個(gè)諧波疊加，因此需要選取合適的諧波數(shù)量進(jìn)行擬合：

Y（t）=a0+[akcos（kωt）+bksin（kωt）] （6）

寫成矩陣形式為AX=Y，其中：

X=[a0，a1，b1，a2，b2，…，am，bm]T

Y=[Y（t1），Y（t2），…，Y（tn）]T

A= 1 ? cos（ωt1） … ?cos（mωt1） ?sin（mωt1） ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 1 ? cos（ωtn） ?… ? cos（mωtn） ? sin（mωtn）

采用最小二乘法求解為：

X=（ATA）-1ATY ? （7）

據(jù)此可以求出傅里葉級(jí)數(shù)擬合方程的系數(shù)，進(jìn)而求出信號(hào)零點(diǎn)。

1.3 ? ?基于卡爾曼濾波器的誤差校正

在平穩(wěn)單相流信號(hào)條件下，相位差保持在一個(gè)恒定的值，經(jīng)過計(jì)算求解之后的相位差通常是在真實(shí)值附近上下波動(dòng)，因此可以采用卡爾曼濾波器對(duì)所求解的誤差值進(jìn)行預(yù)測(cè)估計(jì)，使其在真實(shí)值附近減小波動(dòng)。卡爾曼濾波主要公式如下[8]：

狀態(tài)先驗(yàn)估計(jì)：

k | k-1=Φk，k-1k-1 | k-1 ? （8）

方差先驗(yàn)估計(jì)：

Pk | k-1=Φk，k-1Pk-1 | k-1Φk，k-1 ?T ? ?+？祝k，k-1Qk，k-1？祝k，k-1 ?T ? ? ? ? ? （9）

增益矩陣：

Kk=Pk | k-1HkT（HkPk | k-1HkT+Rk）-1 ? ?（10）

狀態(tài)估計(jì)校正（最優(yōu)估計(jì)）：

k | k=k | k-1+Kk（zk-Hkk | k-1） ?（11）

方差估計(jì)校正：

Pk | k=（I-KkHk）Pk | k-1 ? （12）

卡爾曼濾波器主要包括兩個(gè)過程：預(yù)測(cè)和校正。狀態(tài)先驗(yàn)估計(jì)和方差先驗(yàn)估計(jì)為時(shí)間更新方程，用于預(yù)測(cè)部分;增益矩陣、狀態(tài)估計(jì)校正和方差估計(jì)校正為狀態(tài)更新方程，用于校正部分。通過這兩組方程便可以進(jìn)行誤差校正。

2 ? ?仿真分析

由于科氏力流量計(jì)的相位差通常在4°～10°，因此在此次仿真中選擇了5 000個(gè)采樣點(diǎn)用來模擬信號(hào)，且fs=16 000 Hz，A=2.5，f0=188.64 Hz，Phase difference=

4.01°，SNR=38.96 dB。

在進(jìn)行計(jì)算時(shí)考慮DSP處理器一次處理的采樣點(diǎn)數(shù)不能太多，因此5 000個(gè)采樣點(diǎn)每次會(huì)處理1 000個(gè)，將1 000個(gè)采樣點(diǎn)內(nèi)的零點(diǎn)求出之后進(jìn)行時(shí)間差的計(jì)算，然后將該范圍內(nèi)的零點(diǎn)求平均：

Δt=Δti ? ?（13）

式中：Δt為獲得的時(shí)間差輸出。

最后通過卡爾曼濾波器進(jìn)行校正，校正后的結(jié)果及誤差如圖4所示。

經(jīng)過卡爾曼濾波器校正之后的均方誤差如表1所示。

從仿真結(jié)果中可以看出，由于過零檢測(cè)無收斂過程，因此可以很好地滿足實(shí)時(shí)性的要求，同時(shí)經(jīng)過卡爾曼濾波之后，所求的相位差的均方誤差要遠(yuǎn)小于濾波之前的均方誤差，提高了相位差計(jì)算結(jié)果的精度。

3 ? ?結(jié)語

本文基于平穩(wěn)單相流信號(hào)模型，提出了改進(jìn)的科氏力流量計(jì)相位差的計(jì)算方法。首先采用FFT預(yù)估信號(hào)頻率，減少了自適應(yīng)格型濾波器的收斂時(shí)間，同時(shí)引入了反饋因子解決了因頻率變化導(dǎo)致的自適格型濾波器無法收斂的問題;然后采用傅里葉級(jí)數(shù)擬合的方式進(jìn)行零點(diǎn)的求解;最后采用卡爾曼濾波器進(jìn)行誤差的校正。從仿真結(jié)果可以看出，卡爾曼濾波器很好地修正了相位差計(jì)算的誤差，使精度得到了提高。

[參考文獻(xiàn)]

[1] LI M，HENRY M.Complex signal processing for Coriolis mass flow metering in two-phase flow[J] . Flow Measurement and Instrumentation ，2018 ，

64：104-115.

[2] LI M，HENRY M.Complex bandpass filtering for Coriolis mass flow meter signal processing[C]//

IECON 2016-42nd Annual Conference of the IEEE Industrial Electronics Society，2016.

[3] 李葉，徐科軍，朱志海，等.面向時(shí)變的科里奧利質(zhì)量流量計(jì)信號(hào)的處理方法研究與實(shí)現(xiàn)[J].儀器儀表學(xué)報(bào)，2010，31（1）：8-14.

[4] 黃丹平，汪俊其，于少東，等.基于小波變換和改進(jìn)Hilbert變換對(duì)科氏質(zhì)量流量計(jì)信號(hào)處理[J].中國(guó)測(cè)試，2016，42（6）：37-41.

[5] 楊輝躍，涂亞慶，張海濤，等.一種基于SVD和Hilbert變換的科氏流量計(jì)相位差測(cè)量方法[J].儀器儀表學(xué)報(bào)，2012，33（9）：2101-2107.

[6] 涂亞慶，沈廷鰲，李明，等.基于多次互相關(guān)的非整周期信號(hào)相位差測(cè)量算法[J].儀器儀表學(xué)報(bào)，2014，35（7）：1578-1585.

[7] HOU Q L，XU K J，F(xiàn)ANG M，et al.Development of Coriolis mass flow meter with digital drive and signal processing technology[J].ISA Trans-

actions，2013，52（5）：692-700.

[8] 彭丁聰.卡爾曼濾波的基本原理及應(yīng)用[J].軟件導(dǎo)刊，2009，8（11）：32-34.

收稿日期：2021-11-09

作者簡(jiǎn)介：徐長(zhǎng)利（1999—），男，山東臨沂人，碩士研究生，研究方向：智能制造。

通信作者：李波（1975—），男，四川人，博士，教授，研究方向：智能制造。