徐科軍， 王 沁， 方 敏， 邵春莉

（1.合肥工業(yè)大學(xué)電氣與自動化工程學(xué)院，安徽合肥 230009；2.工業(yè)自動化安徽省工程技術(shù)研究中心，安徽合肥 230009）

基于數(shù)字濾波與自相關(guān)的渦街流量計抗強(qiáng)振動方法改進(jìn)

徐科軍1，2， 王 沁1， 方 敏1， 邵春莉1

（1.合肥工業(yè)大學(xué)電氣與自動化工程學(xué)院，安徽合肥 230009；2.工業(yè)自動化安徽省工程技術(shù)研究中心，安徽合肥 230009）

渦街流量計容易受到管道振動的干擾，對測量精度造成較大影響，特別是當(dāng)振動噪聲的功率大于渦街信號的功率時，采用通常的數(shù)字信號處理方法也無法得到正確的渦街頻率。為此，提出將數(shù)字濾波和自相關(guān)相結(jié)合的算法，運(yùn)用于渦街流量傳感器的信號處理，從含有強(qiáng)振動干擾的信號中識別出流量信息。針對該方法在實際實現(xiàn)中遇到的問題，對已有系統(tǒng)軟件進(jìn)行改進(jìn)和完善?；贛SP430F5418單片機(jī)，使用匯編語言編制實數(shù)FFT算法實現(xiàn)了自相關(guān)函數(shù)的運(yùn)算。

計量學(xué)；渦街流量計；數(shù)字濾波；自相關(guān)；強(qiáng)振動干擾；FFT；單片機(jī)

1 引 言

渦街流量計由于無機(jī)械可動部件，安裝方便，量程比較寬，并且可以測量氣體和液體介質(zhì)等優(yōu)點，被廣泛應(yīng)用于石油、化工等領(lǐng)域。但是，在實際應(yīng)用中，由于渦街流量計普遍采用應(yīng)力式傳感方式，極易受到管道振動的干擾，無法保證測量精度，有時還出現(xiàn)錯誤的結(jié)果［1］。采用數(shù)字信號處理方法可以從含有噪聲的信號中提取流量信息［2～6］。但是，當(dāng)振動噪聲的功率大于渦街信號的功率時，普通的數(shù)字信號處理方法也無法得到正確的流量信息。為了解決這個問題，文獻(xiàn)［7］提出了基于數(shù)字濾波和自相關(guān)算法相結(jié)合的抗強(qiáng)振動干擾處理方法，并用超低功耗單片機(jī)實現(xiàn)，取得較好的實驗結(jié)果。但是，要在實際應(yīng)用中取得好的應(yīng)用效果，還需要進(jìn)一步改進(jìn)和完善。為此，本文基于MSP430F5418單片機(jī)，編制匯編語言的實數(shù)FFT算法完成自相關(guān)函數(shù)的運(yùn)算，提高計算速度；改進(jìn)自相關(guān)衰減系數(shù)的計算方法，保證測量精度。

2 抗強(qiáng)振動干擾方法

2.1 基本原理

渦街流量信號在理想情況下為規(guī)則穩(wěn)定的正弦信號，渦街體積流量與該正弦信號頻率成正比。但是，在實際工作環(huán)境下，由于水泵或其它因素的影響均會引起管道內(nèi)流場不穩(wěn)定，因此實際的渦街流量信號并不是完全穩(wěn)定的正弦波，其頻率會因為干擾的影響存在晃動。而強(qiáng)振動信號一般由機(jī)泵或閥門開關(guān)等現(xiàn)場設(shè)備的機(jī)械振動引起，因此，其頻率為現(xiàn)場設(shè)備的振動頻率或其倍頻，一般較固定。文獻(xiàn)［8］通過對實驗現(xiàn)場采集的大量數(shù)據(jù)進(jìn)行分析和建模，以功率譜計算為基礎(chǔ)，定性地分析了渦街信號和振動信號的帶寬特征，得出渦街信號為寬帶信號，而振動信號為窄帶信號的結(jié)論。文獻(xiàn)［7］分別對強(qiáng)振動信號和渦街信號進(jìn)行自相關(guān)計算，得到不同的衰減特性。自相關(guān)計算可以反映數(shù)組序列不同時刻數(shù)據(jù)的相關(guān)性，因為強(qiáng)振動噪聲信號為窄帶信號，其頻率值波動范圍較小，而渦街信號為寬帶信號，頻率波動較大，因此強(qiáng)振動噪聲的相關(guān)性要優(yōu)于渦街信號，這一特征在自相關(guān)函數(shù)上表現(xiàn)為，渦街信號的自相關(guān)函數(shù)隨著時間的延時，其幅值有較大程度的衰減，見圖1；而強(qiáng)振動噪聲則衰減較小，見圖2。

圖1 渦街信號自相關(guān)衰減函數(shù)圖

圖2 強(qiáng)振動信號自相關(guān)衰減函數(shù)圖

為此，文獻(xiàn)［7］提出了基于數(shù)字濾波和自相關(guān)相結(jié)合的渦街流量傳感器抗強(qiáng)干擾信號處理方法，具體處理過程為：首先對采集到的時域信號進(jìn)行頻譜分析，按照一定的信噪比及幅值門限確定頻率個數(shù)，可能為0，可能為1，也可能大于1。當(dāng)頻率個數(shù)為0時，系統(tǒng)輸出0表示沒有流量。當(dāng)頻率個數(shù)為1時，會有3種情況：（1）只有渦街流量；（2）只有噪聲；（3）渦街流量信號與噪聲信號同頻?？梢郧笤盘栄舆t一段時間后的自相關(guān)函數(shù)在半個周期內(nèi)的幅值與在時間為0時的自相關(guān)函數(shù)的值的比值，稱該比值為衰減系數(shù)，當(dāng)衰減系數(shù)小于設(shè)定的閾值時，即只有渦街流量信號或渦街流量信號與噪聲信號同頻，輸出渦街流量頻率；當(dāng)衰減系數(shù)大于等于閾值時，即只有噪聲，輸出0表示沒有流量。當(dāng)頻率的個數(shù)大于1時，要先采用數(shù)字帶通濾波器對信號進(jìn)行濾波，將濾波以后得到的相應(yīng)頻率的信號進(jìn)行自相關(guān)計算，然后計算0.5 s附近的正弦信號波峰值與自相關(guān)函數(shù)的第一點的比值，然后判斷不同頻率點信號所得到的衰減系數(shù)的大小，其中衰減系數(shù)最小的頻率點對應(yīng)的信號即為渦街信號。

2.2 幾點局限

該方法根據(jù)振動信號和渦街信號自相關(guān)函數(shù)的衰減程度的不同，作為去除噪聲信號的依據(jù)，在一定程度上解決了抗強(qiáng)振動干擾的問題。但是，由于為了滿足低功耗、兩線制過程控制儀表的要求，無法選用DSP去實現(xiàn)算法，只能選用低功耗單片機(jī)，而單片機(jī)內(nèi)存單元和計算速度都很有限，所以，該方法在實時實現(xiàn)過程中還存在以下3方面的問題：

（1）該算法中采用了數(shù)字帶通濾波器，對原始的渦街流量傳感器輸出信號進(jìn)行濾波，得到相應(yīng)頻率點的信號。在Matlab仿真中，設(shè)定的帶通濾波器帶寬較窄，僅為［－5 Hz，5 Hz］。由于單片機(jī)的內(nèi)存單元有限，只能進(jìn)行2048點的濾波，而數(shù)字帶通濾波器的帶寬較窄，因此會有較長時間的收斂過程，如圖3和圖4所示。

圖3 帶通濾波后振動信號

當(dāng)渦街信號頻率和強(qiáng)振動信號頻率值很靠近時，經(jīng)過帶寬很窄的數(shù)字帶通濾波器濾波后，信號會出現(xiàn)較大幅度的波動；并且經(jīng)過濾波后的振動信號以及渦街信號在前300個點均為收斂過程，同時，整個范圍內(nèi)幅值也存在失真的情況，這對下一步的自相關(guān)計算造成不良影響。

圖4 帶通濾波后渦街信號

（2）自相關(guān)函數(shù)可以描述隨機(jī)過程中任意2個不同時刻之間相關(guān)性，反映函數(shù)在2個不同時刻的狀態(tài)之間的統(tǒng)計關(guān)聯(lián)程度。自相關(guān)估計一般需要點數(shù)較多，在Matlab中使用了16384個點來進(jìn)行仿真，但是，在系統(tǒng)實時實現(xiàn)過程中單片機(jī)內(nèi)存的限制，只能進(jìn)行2048個點的自相關(guān)計算，點數(shù)減少，從而降低自相關(guān)計算的準(zhǔn)確性。此外，用來計算自相關(guān)的序列是經(jīng)過數(shù)字濾波得到的，而當(dāng)數(shù)字濾波的結(jié)果不能真實反映相應(yīng)頻率點信號時，也會對自相關(guān)函數(shù)造成影響。圖5和圖6分別為振動信號和渦街信號的2048個點的自相關(guān)函數(shù)圖。由圖5可以看出，振動信號的自相關(guān)函數(shù)也存在一定的衰減，而由6可以看出，渦街信號的自相關(guān)函數(shù)的衰減曲線存在波動。

圖5 振動信號自相關(guān)函數(shù)

（3）該方法采用無偏估計進(jìn)行自相關(guān)計算：

圖6 渦街信號自相關(guān)函數(shù)

式中：R為數(shù)據(jù)個數(shù)；r為延遲時間；x（n）為濾波后的原始信號；n＝0，1，2，…，R－r－1。

采用式（1）的方法乘法次數(shù)多，運(yùn)算量較大。因此，為了節(jié)省計算時間，提高系統(tǒng)實時性，原來的系統(tǒng)中根據(jù)振動自相關(guān)函數(shù)和渦街信號自相關(guān)函數(shù)的特點進(jìn)行了計算上的簡化，因為在仿真中振動信號和渦街信號自相關(guān)函數(shù)的衰減規(guī)律都是單向衰減，并且振動信號自相關(guān)函數(shù)衰減較慢，因此我們一般選擇0.5 s延時附近的峰值點與自相關(guān)函數(shù)的第一個峰值點進(jìn)行比較，從而計算出振動信號和渦街信號的衰減系數(shù)，然后進(jìn)行比較。即按照式（1）只計算了1個點的自相關(guān)系數(shù)。但是，經(jīng)過分析，由圖5和圖6所示振動信號或者渦街信號的自相關(guān)函數(shù)并不完全是單向衰減的，其原因是濾波后信號的失真以及自相關(guān)點數(shù)不夠多，因此原來軟件中只計算0.5 s延時處的一點的自相關(guān)系數(shù)，會造成結(jié)果的誤判。

3 方法的改進(jìn)

針對上述局限性，對基于單傳感器的抗強(qiáng)干擾系統(tǒng)軟件進(jìn)行了相應(yīng)的改進(jìn)。由于MSP430F5418單片機(jī)的內(nèi)存為16 kB，為了能夠在單片機(jī)上實現(xiàn)抗強(qiáng)干擾算法，自相關(guān)計算點數(shù)的上限為2048個點。因此不再僅僅計算0.5 s延時處的自相關(guān)系數(shù)，而是將2048個點的自相關(guān)系數(shù)都計算出來，然后在找出第一個最大衰減的點數(shù)，根據(jù)這一點與第1個峰值點的衰減比值來消除噪聲。

3.1 采用實數(shù)FFT實現(xiàn)自相關(guān)計算

基于MSP430單片機(jī)用實數(shù)FFT方法實現(xiàn)自相關(guān)計算；自相關(guān)函數(shù)的無偏估計公式如式（1）所示，如果計算過程中R和r均較大，則計算過程非常繁瑣，占用時間較長，因此很難在單片機(jī)中進(jìn)行應(yīng)用。自相關(guān)算法還有另一種估計方法即有偏估計方法：

式中~Rxx的點數(shù)為2R－1，對于規(guī)定的R值，當(dāng)r遠(yuǎn)小于R時，~Rxx估計值才越接近真值。對~Rxx求傅里葉變換，則有

式中ω為圓周頻率。

又由于2個R點的序列進(jìn)行線性卷積，將得到一個長為2R－1點的序列，如果用DFT來實現(xiàn)線性卷積，則需要將序列補(bǔ)零延長至2R－1點，從而得到序列x2R（n）。由于x2R（n）的后R－1點均為0，因此，式（3）可以進(jìn)行如下變換：

設(shè)m＝r＋n，且0≤m≤2R－1，則式（4）右邊可以寫成如下所示：

式中X2R（ejω）為序列x2R（n）的傅里葉變換。

由式（3）和式（5）可知：

式（6）說明，序列的有偏自相關(guān)估計和該序列的功率譜是一對傅里葉變換。而傅里葉變換可以通過FFT實現(xiàn)，因此，在單片機(jī)中通過FFT計算實現(xiàn)2048個點的自相關(guān)計算。

在以基于MSP430F5418單片機(jī)為核心的硬件系統(tǒng)上完成軟件研制，因此算法的實時實現(xiàn)需要合理分配單片機(jī)資源。單片機(jī)內(nèi)存為16kB，采樣數(shù)組為2400個點，數(shù)據(jù)類型為整型，占用2 B，保存數(shù)組需要約4kB的內(nèi)存，因此自相關(guān)計算最多只能選擇2048個點。為了能通過FFT實現(xiàn)，首先需要將數(shù)據(jù)補(bǔ)零至4096個點，數(shù)據(jù)類型為整型，需要占用內(nèi)存空間8 kB。

實現(xiàn)過程中利用實數(shù)FFT來進(jìn)行自相關(guān)計算，而實數(shù)FFT是將序列按照奇偶分別存放于奇數(shù)序列和偶數(shù)序列。設(shè)存放原來序列的數(shù)組為x（n），長度為2N，奇數(shù)項h（r）＝x（2r＋1），其中r＝0，1，…，N－1，奇數(shù)項作為FFT計算的虛部；偶數(shù)項g（r）＝x（2r），其中r＝0，1，…，N，偶數(shù)項作為FFT計算的實部。具體實現(xiàn)過程如下：

（1）將經(jīng)過數(shù)字濾波后的2048個點存放于x（n）中，該序列長度為2048點，然后通過補(bǔ)零的方法將序列x（n）補(bǔ)到4096個點，為了節(jié)省單片機(jī)內(nèi)存，補(bǔ)出的0不是直接放于x（n）中，而是按照奇偶項直接分配到h（r）＝x（2r＋1）和g（r）＝x（2r）中，h（r）為2048點，后1024個點為0，同理g（r）也為2048個點，后1024個點也為0。

（2）h（r）和g（r）序列就緒后，進(jìn)行實數(shù)FFT計算，由于h（r）和g（r）的總點數(shù)為4096個點，因此完成的實際是4096個點的FFT計算，設(shè)得到的傅里葉變換序列為X（k），因為X（k）是關(guān)于中心共軛對稱的，因此只保存前2048個點的實部和虛部，后半部分根據(jù)其共軛對稱的關(guān)系求出，其中序列h（r）和g（r）中分別存放著前2048個點的虛部和實部。同時求出X（k）的共軛序列X*（k）。

（3）根據(jù)第（2）歩計算X（k）和X*（k）的乘積，他們?yōu)楣曹椥蛄?，因此乘積為實數(shù)，又因為X（k）后2048個點與前2048個點為對稱共軛的關(guān)系，因此在這個過程中可以計算出4096個點的乘積值，按順序存放于h（r）和g（r）序列中。

（4）對式（3）中求得的乘積進(jìn)行IFFT計算，其計算公式為：

式中N＝2 048為自相關(guān)實際點數(shù)。其中利用FFT的子程序完成IFFT的計算。

3.2 優(yōu)化自相關(guān)衰減系數(shù)計算方法

為了滿足渦街流量計的低功耗要求，選擇MSP430單片機(jī)作為處理器。但是，目前430系列單片機(jī)最大內(nèi)存僅為18 kB，仍然不能實現(xiàn)多點數(shù)的數(shù)字濾波器和多點數(shù)的自相關(guān)計算，因此根據(jù)圖（5）和圖（6）的特性，在自相關(guān)衰減計算中做了相應(yīng)的調(diào)整，以保證算法的有效性。雖然少點數(shù)的自相關(guān)特性有些變化，但是渦街信號的自相關(guān)函數(shù)仍然比振動信號衰減得更明顯，由于濾波造成的失真，使得渦街信號自相關(guān)函數(shù)會先出現(xiàn)快速衰減的趨勢，然后又再上升；而振動信號的自相關(guān)函數(shù)在越接近N的位置也會衰減比較快。因此不再選擇0.5 s的位置來比較二者的衰減系數(shù)。對于渦街信號，通過逐點尋找，找出第一個衰減最大的正弦峰值點；對于振動信號仍采用此方法，如果自相關(guān)函數(shù)的衰減趨勢為單向的，則仍選擇1024個點附近的正弦峰值點。

4 實驗結(jié)果

分別在重慶川儀自動化股份有限公司和合肥工業(yè)大學(xué)實驗室的振動平臺上做了抗強(qiáng)振動干擾的實驗測試，驗證在存在強(qiáng)振動干擾的情況下改進(jìn)方法和系統(tǒng)的有效性。通過鼓風(fēng)機(jī)吹風(fēng)模擬氣體流量，由振動實驗平臺模擬強(qiáng)振動干擾，并施加于渦街流量計上。分別選擇膜盒式和懸臂梁式渦街傳感器進(jìn)行實驗。

在重慶川儀自動化股份有限公司，針對DN50口徑膜盒式渦街傳感器進(jìn)行抗強(qiáng)振動的實驗。為了驗證抗強(qiáng)干擾算法的有效性，在實驗過程中使振動信號的功率大于渦街信號的功率。以振動信號頻率為110 Hz，渦街信號為100 Hz時的功率譜圖為例，從圖7中可以看出振動信號的功率比渦街信號的大。該圖是用數(shù)字存儲示波器采集渦街流量傳感器輸出信號，用Matlab分析的結(jié)果。在實驗過程中，選擇了3個不同的流量點，分別大約為101 Hz，130 Hz以及148 Hz，在這3個流量點下，通過改變振動臺的振動頻率以及振動加速度來進(jìn)行考核。當(dāng)振動頻率較低時，其振感很強(qiáng)，儀表液晶在振動的作用下其讀數(shù)有些模糊。用基于MSP430F5418單片機(jī)的抗強(qiáng)振動型渦街流量計實時測量和處理，實驗結(jié)果表明，改變振動頻率和加速度，儀表均能顯示正確的渦街頻率，如表1所示。

圖7 膜盒式渦街傳感器輸出信號功率譜

在合肥工業(yè)大學(xué)實驗室搭建的振動實驗平臺上，采用DN50口徑懸臂梁式渦街傳感器進(jìn)行抗強(qiáng)振動實驗。實驗過程中，當(dāng)振動實驗臺加速度為0.8g及以上時，其振動強(qiáng)度能夠?qū)⒎庞谡駝优_表面的小工件慢慢振落下去。同樣地，選擇3個流量點分別大約為101 Hz、88 Hz和121 Hz，圖8所示為用數(shù)字存儲示波器采集渦街流量傳感器輸出信號，用Matlab分析的結(jié)果。用基于MSP430F5418單片機(jī)的抗強(qiáng)振動型渦街流量計實時測量和處理，實驗結(jié)果如表2所示，可見，儀表在不同頻率振動的影響下，均能正確顯示渦街信號頻率。

表1 DN50膜盒式渦街傳感器抗強(qiáng)振動實驗結(jié)果

圖8 懸臂梁式渦街傳感器輸出信號功率譜

表2 DN50懸臂梁式渦街傳感器抗強(qiáng)振動實驗結(jié)果

5 結(jié)束語

根據(jù)現(xiàn)場實驗，對數(shù)字濾波和自相關(guān)算法在單片機(jī)實現(xiàn)過程中存在的局限性進(jìn)行分析，通過示波器采樣數(shù)據(jù)在Matlab中的分析結(jié)果與單片機(jī)計算過程中得到的結(jié)果進(jìn)行比較，分析得出：由于單片機(jī)內(nèi)存較小，只能做2048點的濾波計算，因此濾波結(jié)果存在失真；單片機(jī)資源有限，選擇了2048點自相關(guān)計算，無法達(dá)到Matlab的仿真效果。以現(xiàn)有單片機(jī)資源為考慮因素，進(jìn)行了2方面的改進(jìn)：一方面利用實數(shù)FFT實現(xiàn)自相關(guān)計算，將2048點的自相關(guān)計算快速計算出，節(jié)省計算時間。另一方面改進(jìn)了自相關(guān)衰減系數(shù)計算方法，從而保證算法計算的準(zhǔn)確性。

致謝：感謝重慶川儀自動化股份有限公司提供實驗條件和渦街流量計一次儀表。

［1］ Miau J J，Hu C C，Chou J H.Response of a vortex flowmeter to impulsive vibrations［J］.FlowMeasurementandInstrumentation，2000，11：41－49.

［2］ Schlatter G L，Barrett W D，Waers J F，etal．Signal processing method and apparatus for flowmeters：WO，1990004230［P］.1990－04－19.

［3］ Clark D W.Designing phase-locked loops for instrumentation applications［J］.Measurement，2002，32：205－227.

［4］ Clark D W，Ghaoud T.A dual phase locked loop for vortex flow metering［J］.FlowMeasurementand Instrumentation，2003，14：1－11.

［5］ 曾憲俊，徐科軍，周楊，等.低功耗、兩線制渦街流量計數(shù)字信號處理系統(tǒng)［J］.儀表技術(shù)與傳感器，2009，（1）：76－79.

［6］ 徐科軍，劉三山，劉家祥，等.基于數(shù)字頻譜與模擬帶通濾波組的兩線制渦街流量計［J］.計量學(xué)報，2011，32（1）：25－30.

［7］ Xu K J，Luo Q L，Wang G.Frequency-feature based anti-strong disturbance signal processing method and system for vortex flowmeterwith single sensor［J］.Review ofScientificInstuments，2010，7：075104－7.

［8］ 徐科軍，羅清林，朱永強(qiáng)，等.渦街流量傳感器信號的幅值和頻率特征［J］.計量學(xué)報，2010，31（6）：504－508.

Im provements of Anti-strong-vibration Method Combing Digital Filter w ith Auto-correlation for Vortex Flowmeter

XU Ke-jun1，2， WANG Qin1， FANG Min1SHAO Chun-li1
（1.School of Electrical and Automation Engineering，Hefei University of Technology，Hefei，Anhui230009，China；
2.Engineering Technology Research Center of Industrial Automation，Anhui Province，Hefei，Anhui230009，China）

Vortex flowmeter is susceptible to the interference of pipe vibration，which affects on the measurement accuracy greatly.Adopting normal digital signal processingmethods cannotobtain the accurate vortex frequencieswhen the power of vibration noise is larger than that of the vortex signal.Therefore a digital signal processing technique combining digital filter with autocorrelation was proposed to apply the signal processing of vortex flow sensor for extracting the flow rate information from the signal containing strong vibration interference.Aiming at the implementing problems of thismethod in practice application，the existing system software is improved.The autocorrelation function is calculated by using real number FFT algorithm with assemble language based on MSP430F5418 MCU in order to speed up computation.

Metrology；Vortex flowmeter；Digital filter；Autocorrelation；Strong vibration interference；FFT；MCU

TB937

A

1000-1158（2014）05-0449-06

10．3969／j．issn．1000-1158．2014．05．09

2012-08-10；

2013-01-23

國家自然科學(xué)基金（61074164）

徐科軍（1956－），男，江蘇無錫人，合肥工業(yè)大學(xué)教授、博士生導(dǎo)師，主要從事傳感器技術(shù)、自動化儀表和數(shù)字信號處理研究。dsplab＠hfut.edu.cn