基于高斯和粒子濾波的動態稱重數據處理

張麗妹，高占寶，尹志兵

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

0 引言

近年來國內外稱重技術研究取得了顯著進展，主要表現在稱重測量技術由靜態測量向動態測量、在線測量和模型化方向發展［1］。所謂動態稱重即稱重信號處于動態變化之中，沒有達到穩態，在動態測量中較快地準確地獲得稱重值是目前動態稱重時急需解決的問題。然而由于動態稱重過程中機械振動等影響，導致測量數據存在很大的噪聲，為了得到準確的稱重值需要對測量信號進行濾波處理。

目前常用于處理稱重數據的濾波方法有中值濾波法、滑動平均濾波法等傳統濾波方法，但是這類傳統的濾波方法難以克服測量數據的非線性，對于處理強非平穩數據難以得到較好的結果。為此本文在利用已有測量數據對系統建立模型的基礎上，采用基于模型的濾波方法對測量數據進行處理。在此類基于模型的濾波方法中卡爾曼濾波堪稱經典，卡爾曼濾波是噪聲服從高斯分布的線性系統最優濾波，對非高斯非線性系統，卡爾曼濾波不能達到最優估計效果。粒子濾波是一種基于遞推貝葉斯估計的蒙特卡洛實現方法，在處理非線性非高斯系統問題上有著明顯優勢［2］。

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

1 粒子濾波原理

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

離散的動態狀態空間模型可描述為

式中:xk∈Rn為k時刻狀態變量;zk∈Rm為k時刻得到的測量向量;wk，vk為過程噪聲與觀測噪聲;fk為量測函數;hk為觀測函數。

重抽樣過程雖然在一定程度上解決了粒子退化問題，但是重抽樣會帶來新的問題。例如，重抽樣總是傾向于復制最合適的粒子，并將導致不合適的粒子子代為零。在極端情況下，這又將導致一種粒子退化，損失了粒子的多樣性。

2 高斯和粒子濾波

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

2.1 高斯混合逼近

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

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

2.2 高斯和粒子濾波

在高斯和粒子濾波中，采用高斯混合模型來逼近后驗分布。

在時間更新階段，假設時間的后驗分布為p(xk|y0:k)并且可用高斯和來表達［5］，即

則預估分布為

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

在測量更新過程，當新的測量值yk+1到達時，可通過貝葉斯估計算法更新后驗密度［6］:

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

2)i=1，…，G;j=1，…，M。計算相應權重為

式中:π(xk|y0:k)為建議分布。當假設先驗分布為重要性分布時，權重計算如下:

3)由于

4)更新高斯混合分布的權重為

5)歸一化權重

在此基礎上可得后驗分布為

則狀態xk+1的估計可計算如下

這種方式避免了粒子濾波的重采樣步驟，從而避免了粒子濾波的退化問題，而且高斯和粒子濾波以高斯和分布逼近非高斯噪聲，較粒子濾波更精確。下面通過具體實例數據來驗證這一點。

3 實例驗證

為驗證高斯和粒子濾波在動態稱重數據處理方面的效果，本文采集了大量動態組合秤實測數據，并分別用EKF(擴展卡爾曼濾波)、PF(粒子濾波)、GSPF(高斯和粒子濾波)對數據進行處理，將處理結果進行對比研究?，F以一組60 Hz 采樣頻率、標稱稱重500 g的采集數據為例，說明本文方法的優點。

3.1 稱重系統建模

上述濾波方法需要對系統進行建模。傳統的分析系統物理結構進行建模的方法，過程復雜且需要一定的專業知識，而系統的輸入輸出數據從一定程度上反映了系統的結構，為此本文采用基于數據的建模方法［8］。由于真實稱重過程中會有物料下落的沖擊力帶來的強擾動，所以采用負階躍測得數據對系統進行建模，即在稱量盤上懸掛一定質量重物，在某時刻剪斷繩子測量系統輸出，這樣就可以得到系統的負階躍響應動態校準數據，進而可得到相對精確的模型。

對于上述算法應建立狀態空間模型，采用損失函數及最終預報誤差來衡量模型的精度，以選擇模型的階次。本文比較了二階、三階模型，各自的損失函數及最終預報誤差值如表1所示。由表1 可以看出三階模型較二階模型精度提高很多，且兩個指標數值已經很小，即模型較精確，考慮到計算復雜度問題無需采用更高階模型。

表1 損失函數與最終預報誤差比較

用此方法建立的狀態空間模型如下:

其中，

Xk，Xk-1∈R3分別為k，k-1 時刻狀態變量，Yk∈R 為k時刻觀測向量，Uk∈R 為k時刻輸入向量，ek-1，ek∈R 為過程噪聲及觀測噪聲，均服從伽馬分布。

3.2 仿真結果

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

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

圖1 三種濾波算法結果

圖2 濾波結果細節

為了定量衡量各方法的優劣，采用均方根誤差表示真實狀態與狀態估計的差值［9］:

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

表2 均方根誤差比較結果

上面兩種方法是從濾波效果角度在定性與定量上比較了三種方法的優劣，得出了GSPF 方法最優的結論。然而從實際工程角度，濾波目的是要快速地得到精度較高的稱重數據，即提高動態稱重的速度與精度，為此需要得到能表征速度與精度的量。尺度函數是為衡量動態稱重效果而編寫的函數，用來在數量上衡量稱重的速度與精度，分別用穩定時間與穩定誤差來表示。穩定時間為從開始測量到數據進入穩定域的時間，穩定誤差為穩定域內數據的分散程度，穩定時間越小表示稱重速度越快，穩定誤差越小表示稱重精度越高。三種方法的對比結果如表2。由表2 可以看出EKF 的速度慢、精度差，而PF 的速度太慢，只有GSPF 的速度快且精度高，對動態稱重數據處理有絕對的優勢。

表3 稱重效果對比

4 總結

為了提高動態稱重的快速性與準確性，提高生產效率，針對動態稱重過程的強非平穩性及過程噪聲的復雜性，在對系統建模的基礎上，本文提出了高斯和粒子濾波方法對動態稱重數據進行濾波處理。仿真結果證明高斯和粒子濾波方法可以很好地濾除測量噪聲，從而快速的得到較精確的測量值，并且與EKF 及PF方法作對比證明處理結果優于以上兩種方法。

由于PF 及GSPF 方法是基于蒙特卡洛的統計分析計算方法，在計算量上要比EKF 大得多，但是由于現在計算機處理速度的提高，計算復雜度已不再是阻礙粒子濾波發展的障礙，因此粒子濾波系列算法近年來被廣泛應用于處理非線性非高斯濾波問題。在此，通過仿真證明了該算法在動態稱重數據處理方面的優勢，要將此方法應用于生產過程中還需要做一些工作，可用并行方式利用FPGA(現場可編程門陣列)來在線實現該算法。

［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］.北京:北京交通大學，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.