張麗妹，高占寶，尹志兵

(1.北京航空航天大學(xué)自動(dòng)化科學(xué)與電氣工程學(xué)院，北京100191;2.中山市精威包裝機(jī)械有限公司，廣東中山528425)

0 引言

近年來(lái)國(guó)內(nèi)外稱重技術(shù)研究取得了顯著進(jìn)展，主要表現(xiàn)在稱重測(cè)量技術(shù)由靜態(tài)測(cè)量向動(dòng)態(tài)測(cè)量、在線測(cè)量和模型化方向發(fā)展［1］。所謂動(dòng)態(tài)稱重即稱重信號(hào)處于動(dòng)態(tài)變化之中，沒(méi)有達(dá)到穩(wěn)態(tài)，在動(dòng)態(tài)測(cè)量中較快地準(zhǔn)確地獲得稱重值是目前動(dòng)態(tài)稱重時(shí)急需解決的問(wèn)題。然而由于動(dòng)態(tài)稱重過(guò)程中機(jī)械振動(dòng)等影響，導(dǎo)致測(cè)量數(shù)據(jù)存在很大的噪聲，為了得到準(zhǔn)確的稱重值需要對(duì)測(cè)量信號(hào)進(jìn)行濾波處理。

目前常用于處理稱重?cái)?shù)據(jù)的濾波方法有中值濾波法、滑動(dòng)平均濾波法等傳統(tǒng)濾波方法，但是這類(lèi)傳統(tǒng)的濾波方法難以克服測(cè)量數(shù)據(jù)的非線性，對(duì)于處理強(qiáng)非平穩(wěn)數(shù)據(jù)難以得到較好的結(jié)果。為此本文在利用已有測(cè)量數(shù)據(jù)對(duì)系統(tǒng)建立模型的基礎(chǔ)上，采用基于模型的濾波方法對(duì)測(cè)量數(shù)據(jù)進(jìn)行處理。在此類(lèi)基于模型的濾波方法中卡爾曼濾波堪稱經(jīng)典，卡爾曼濾波是噪聲服從高斯分布的線性系統(tǒng)最優(yōu)濾波，對(duì)非高斯非線性系統(tǒng)，卡爾曼濾波不能達(dá)到最優(yōu)估計(jì)效果。粒子濾波是一種基于遞推貝葉斯估計(jì)的蒙特卡洛實(shí)現(xiàn)方法，在處理非線性非高斯系統(tǒng)問(wèn)題上有著明顯優(yōu)勢(shì)［2］。

對(duì)于卡爾曼濾波及擴(kuò)展卡爾曼濾波，需假設(shè)分布為高斯分布，若真實(shí)分布為非高斯分布，則單一高斯分布難以很好地描述，可考慮用高斯和分布來(lái)近似真實(shí)分布。為此本文提出了基于高斯和粒子濾波的處理方法，并對(duì)稱重系統(tǒng)進(jìn)行建模。，仿真結(jié)果證明采用高斯和濾波進(jìn)行處理的方法在提高稱重的快速性與精確性上有很好的效果。

1 粒子濾波原理

粒子濾波是遞推貝葉斯估計(jì)的一種蒙特卡洛實(shí)現(xiàn)方法，讓樣本粒子通過(guò)系統(tǒng)模型按時(shí)間順序向前傳播，得到各時(shí)刻系統(tǒng)的狀態(tài)樣本，從而求得系統(tǒng)狀態(tài)的后驗(yàn)概率密度函數(shù)。序貫蒙特卡洛方法的核心思想是用一組帶權(quán)隨機(jī)樣本集代表后驗(yàn)分布函數(shù)，并基于這些樣本計(jì)算概率值。

離散的動(dòng)態(tài)狀態(tài)空間模型可描述為

式中:xk∈Rn為k時(shí)刻狀態(tài)變量;zk∈Rm為k時(shí)刻得到的測(cè)量向量;wk，vk為過(guò)程噪聲與觀測(cè)噪聲;fk為量測(cè)函數(shù);hk為觀測(cè)函數(shù)。

重抽樣過(guò)程雖然在一定程度上解決了粒子退化問(wèn)題，但是重抽樣會(huì)帶來(lái)新的問(wèn)題。例如，重抽樣總是傾向于復(fù)制最合適的粒子，并將導(dǎo)致不合適的粒子子代為零。在極端情況下，這又將導(dǎo)致一種粒子退化，損失了粒子的多樣性。

2 高斯和粒子濾波

為了避免上述粒子貧化問(wèn)題，提出了高斯和粒子濾波方法，該方法結(jié)合高斯混合逼近后驗(yàn)分布的思想，并采用蒙特卡洛采樣和更新策略。

2.1 高斯混合逼近

可證明任何概率分布p(x)可由如下密度pG(x)無(wú)限逼近:

式中:G為混合部分的個(gè)數(shù);Φ={α1，…，αG，φ1，…，φG}為高斯混合模型的參數(shù)，且φi={μi，∑i};αi為混合權(quán)重;N(x;μi，∑i)表示均值為μi方差為∑i的正態(tài)分布。逼近誤差可為任意?。?］。

2.2 高斯和粒子濾波

在高斯和粒子濾波中，采用高斯混合模型來(lái)逼近后驗(yàn)分布。

在時(shí)間更新階段，假設(shè)時(shí)間的后驗(yàn)分布為p(xk|y0:k)并且可用高斯和來(lái)表達(dá)［5］，即

則預(yù)估分布為

積分可由高斯混合模型逼近為

在測(cè)量更新過(guò)程，當(dāng)新的測(cè)量值yk+1到達(dá)時(shí)，可通過(guò)貝葉斯估計(jì)算法更新后驗(yàn)密度［6］:

1)對(duì)i=1，…，G，從重要性分布中采樣M個(gè)樣本，記作。

2)i=1，…，G;j=1，…，M。計(jì)算相應(yīng)權(quán)重為

式中:π(xk|y0:k)為建議分布。當(dāng)假設(shè)先驗(yàn)分布為重要性分布時(shí)，權(quán)重計(jì)算如下:

3)由于

4)更新高斯混合分布的權(quán)重為

5)歸一化權(quán)重

在此基礎(chǔ)上可得后驗(yàn)分布為

則狀態(tài)xk+1的估計(jì)可計(jì)算如下

這種方式避免了粒子濾波的重采樣步驟，從而避免了粒子濾波的退化問(wèn)題，而且高斯和粒子濾波以高斯和分布逼近非高斯噪聲，較粒子濾波更精確。下面通過(guò)具體實(shí)例數(shù)據(jù)來(lái)驗(yàn)證這一點(diǎn)。

3 實(shí)例驗(yàn)證

為驗(yàn)證高斯和粒子濾波在動(dòng)態(tài)稱重?cái)?shù)據(jù)處理方面的效果，本文采集了大量動(dòng)態(tài)組合秤實(shí)測(cè)數(shù)據(jù)，并分別用EKF(擴(kuò)展卡爾曼濾波)、PF(粒子濾波)、GSPF(高斯和粒子濾波)對(duì)數(shù)據(jù)進(jìn)行處理，將處理結(jié)果進(jìn)行對(duì)比研究?，F(xiàn)以一組60 Hz 采樣頻率、標(biāo)稱稱重500 g的采集數(shù)據(jù)為例，說(shuō)明本文方法的優(yōu)點(diǎn)。

3.1 稱重系統(tǒng)建模

上述濾波方法需要對(duì)系統(tǒng)進(jìn)行建模。傳統(tǒng)的分析系統(tǒng)物理結(jié)構(gòu)進(jìn)行建模的方法，過(guò)程復(fù)雜且需要一定的專業(yè)知識(shí)，而系統(tǒng)的輸入輸出數(shù)據(jù)從一定程度上反映了系統(tǒng)的結(jié)構(gòu)，為此本文采用基于數(shù)據(jù)的建模方法［8］。由于真實(shí)稱重過(guò)程中會(huì)有物料下落的沖擊力帶來(lái)的強(qiáng)擾動(dòng)，所以采用負(fù)階躍測(cè)得數(shù)據(jù)對(duì)系統(tǒng)進(jìn)行建模，即在稱量盤(pán)上懸掛一定質(zhì)量重物，在某時(shí)刻剪斷繩子測(cè)量系統(tǒng)輸出，這樣就可以得到系統(tǒng)的負(fù)階躍響應(yīng)動(dòng)態(tài)校準(zhǔn)數(shù)據(jù)，進(jìn)而可得到相對(duì)精確的模型。

對(duì)于上述算法應(yīng)建立狀態(tài)空間模型，采用損失函數(shù)及最終預(yù)報(bào)誤差來(lái)衡量模型的精度，以選擇模型的階次。本文比較了二階、三階模型，各自的損失函數(shù)及最終預(yù)報(bào)誤差值如表1所示。由表1 可以看出三階模型較二階模型精度提高很多，且兩個(gè)指標(biāo)數(shù)值已經(jīng)很小，即模型較精確，考慮到計(jì)算復(fù)雜度問(wèn)題無(wú)需采用更高階模型。

表1 損失函數(shù)與最終預(yù)報(bào)誤差比較

用此方法建立的狀態(tài)空間模型如下:

其中，

Xk，Xk-1∈R3分別為k，k-1 時(shí)刻狀態(tài)變量，Yk∈R 為k時(shí)刻觀測(cè)向量，Uk∈R 為k時(shí)刻輸入向量，ek-1，ek∈R 為過(guò)程噪聲及觀測(cè)噪聲，均服從伽馬分布。

3.2 仿真結(jié)果

為比較濾波效果，本文采用三種方式來(lái)衡量濾波效果:圖形法直觀定性比較、均方根誤差定量比較、尺度函數(shù)定量比較。

圖形法可以直觀地看到各方法的處理效果，圖形結(jié)果如圖1、圖2所示。圖1 為三種濾波算法的結(jié)果圖，可以看到PF 和GSPF 與真實(shí)狀態(tài)很接近;圖2 為圖1 的局部細(xì)節(jié)圖。由圖2 可以直觀看出PF 與GSPF較之EKF 方法降噪效果要好很多，EKF 的處理效果最差，GSPF 的處理效果最好。

圖1 三種濾波算法結(jié)果

圖2 濾波結(jié)果細(xì)節(jié)

為了定量衡量各方法的優(yōu)劣，采用均方根誤差表示真實(shí)狀態(tài)與狀態(tài)估計(jì)的差值［9］:

式中:n為觀測(cè)點(diǎn)數(shù);括號(hào)內(nèi)為在t觀測(cè)點(diǎn)處的絕對(duì)誤差。三種方法的均方根誤差比較見(jiàn)表1。由表1 可以看出，EKF 的均方根誤差較其他兩個(gè)大很多，這是由于初始部分大噪聲導(dǎo)致的，這一點(diǎn)在圖1 中可以清楚地看到，而PF 與GSPF 兩方法的誤差結(jié)果都較小，GSPF的均方根誤差最小，估計(jì)狀態(tài)與真實(shí)狀態(tài)最接近，效果最好，這與圖1、圖2 得出的結(jié)果一致，從數(shù)量上更清楚地看到各方法效果的差距。

表2 均方根誤差比較結(jié)果

上面兩種方法是從濾波效果角度在定性與定量上比較了三種方法的優(yōu)劣，得出了GSPF 方法最優(yōu)的結(jié)論。然而從實(shí)際工程角度，濾波目的是要快速地得到精度較高的稱重?cái)?shù)據(jù)，即提高動(dòng)態(tài)稱重的速度與精度，為此需要得到能表征速度與精度的量。尺度函數(shù)是為衡量動(dòng)態(tài)稱重效果而編寫(xiě)的函數(shù)，用來(lái)在數(shù)量上衡量稱重的速度與精度，分別用穩(wěn)定時(shí)間與穩(wěn)定誤差來(lái)表示。穩(wěn)定時(shí)間為從開(kāi)始測(cè)量到數(shù)據(jù)進(jìn)入穩(wěn)定域的時(shí)間，穩(wěn)定誤差為穩(wěn)定域內(nèi)數(shù)據(jù)的分散程度，穩(wěn)定時(shí)間越小表示稱重速度越快，穩(wěn)定誤差越小表示稱重精度越高。三種方法的對(duì)比結(jié)果如表2。由表2 可以看出EKF 的速度慢、精度差，而PF 的速度太慢，只有GSPF 的速度快且精度高，對(duì)動(dòng)態(tài)稱重?cái)?shù)據(jù)處理有絕對(duì)的優(yōu)勢(shì)。

表3 稱重效果對(duì)比

4 總結(jié)

為了提高動(dòng)態(tài)稱重的快速性與準(zhǔn)確性，提高生產(chǎn)效率，針對(duì)動(dòng)態(tài)稱重過(guò)程的強(qiáng)非平穩(wěn)性及過(guò)程噪聲的復(fù)雜性，在對(duì)系統(tǒng)建模的基礎(chǔ)上，本文提出了高斯和粒子濾波方法對(duì)動(dòng)態(tài)稱重?cái)?shù)據(jù)進(jìn)行濾波處理。仿真結(jié)果證明高斯和粒子濾波方法可以很好地濾除測(cè)量噪聲，從而快速的得到較精確的測(cè)量值，并且與EKF 及PF方法作對(duì)比證明處理結(jié)果優(yōu)于以上兩種方法。

由于PF 及GSPF 方法是基于蒙特卡洛的統(tǒng)計(jì)分析計(jì)算方法，在計(jì)算量上要比EKF 大得多，但是由于現(xiàn)在計(jì)算機(jī)處理速度的提高，計(jì)算復(fù)雜度已不再是阻礙粒子濾波發(fā)展的障礙，因此粒子濾波系列算法近年來(lái)被廣泛應(yīng)用于處理非線性非高斯濾波問(wèn)題。在此，通過(guò)仿真證明了該算法在動(dòng)態(tài)稱重?cái)?shù)據(jù)處理方面的優(yōu)勢(shì)，要將此方法應(yīng)用于生產(chǎn)過(guò)程中還需要做一些工作，可用并行方式利用FPGA(現(xiàn)場(chǎng)可編程門(mén)陣列)來(lái)在線實(shí)現(xiàn)該算法。

［1］Nledzwlecki M，Wasilewski A.Application of adaptive filtering to dynamic weighing of vehicles［J］.Control Eng.Practice，1996，5(4):635-644.

［2］朱林富，張三同.基于粒子濾波的故障診斷方法研究［D］.北京:北京交通大學(xué)，2010.

［3］Liu Qinming，Dong Ming，Peng Ying.A novel method for online health prognosis of equipment based on hidden semi- Markov model using sequential Monte Carlo methods［J］.Mechanical Systems and Signal Processing，2012，5.

［4］Shaodong YANG，Desheng WEN，Jing SUN，Junyong MA.Gaussian Sum Particle Filter for Spacecraft Attitude Estimation［C］//2010 the 2nd International Conference on Signal Processing Systems.Dalian:IEEE，2010，3:566-570.

［5］Kotecha J H，Djuric P M.Gaussian Sum Particle Filering［J］.IEEE Transactions on Signal Processing，2003，51(10):2602-2612.

［6］Oshman Y，Carmi A.Attitude Estimation from Vector Observation Using a Genetic-Algorithm-Embedded Quaternion Particle Filter［J］.Journal of Guidance，Control and Dynamics，2006.7(29):879-891.

［7］Kotecha J H，Djuric P M.Gaussian sum particle filtering for dynamic state space models［C］.//Proceedings of ICASSP-2001.Salt Lake City:IEEE，2001，5:3465-3468.

［8］Chaochao Chen，George Vachtsevanos，Marcos E.Orchard.Machine remaining useful life prediction:An integrated adaptive neuro-fuzzy and high-order particle filtering approach［J］.Mechanical Systems and Signal Processing，2012，28:597-607.

［9］Wahyu Caesarendra，Gang Niu，Bo-Suk Yang.Machine condition prognosis based on sequential Monte Carlo method［J］.Expert Systems with Applications，2010，37:2412-2420.