基于改進Kalman的傳感器數(shù)據(jù)加權(quán)處理算法

陳艷春，達鈺鵬

(1.石家莊鐵道大學，河北 石家莊 050043；2.河北省人力資源和社會保障廳信息中心，河北 石家莊 050071)

0 引 言

任何傳感器的測量都是帶有誤差的，誤差產(chǎn)生的原因既有設備本身的問題，在數(shù)據(jù)采集過程中，如受傳感器老化，也有轉(zhuǎn)換器以及無線電傳輸過程中的干擾，使得接收數(shù)據(jù)中經(jīng)常會產(chǎn)生異常跳變點，這種偏離的數(shù)據(jù)點被稱為異常值[1]。異常值分為兩類，一種是孤立的虛假異常值，這類異常值一般是孤立出現(xiàn)，屬于噪音的一種；另一種是真實的異常值，這類異常值一般連續(xù)出現(xiàn)，反映觀測對象發(fā)生異常變化。

在物聯(lián)網(wǎng)監(jiān)測中，為了兼顧計算成本和計算量，常常采用SPC、PCA等方法對建筑等監(jiān)測對象進行異常識別[2-3]。此類基于統(tǒng)計的方法可以基于歷史數(shù)據(jù)實現(xiàn)對異常的識別，其理論基礎是經(jīng)過實踐檢驗的[4]。但這類算法對輸入數(shù)據(jù)的準確性要求較高，直接使用未經(jīng)處理的原始數(shù)據(jù)極易發(fā)生誤報。同時，因為使用場景的不同，很多傳感器的觀測誤差存在周期性變化，處理算法如果采用靜態(tài)參數(shù)，則無法擬合周期性變化。

為了解決實時監(jiān)控中由于虛假異常值出現(xiàn)產(chǎn)生的誤報問題，常用的手段有兩個[5-6]。一是在同一位置部署多個同類型傳感器同時采樣，利用傳感器誤差近似正態(tài)分布的特點，利用統(tǒng)計學方法中的拉伊達法或格拉布斯法剔除虛假異常點后計算平均值，可以實現(xiàn)對噪音和孤立虛假異常點的較好過濾，當傳感器數(shù)量較多時使用拉伊達法，較少時使用格拉布斯法。這類基于統(tǒng)計理論的方法的優(yōu)勢是計算方法簡便，通用性強，魯棒性好，在系統(tǒng)邊緣節(jié)點便可部署使用，天然適合分布式部署，但缺點是需要在單個位置部署大量傳感器，系統(tǒng)硬件部署和維護成本較高。二是利用各種濾波手段對單個傳感器數(shù)據(jù)進行實時修正。例如小波分析、Kalman濾波等，但這類濾波手段也存在著一些不足。以Kalman濾波算法為例，Kalman濾波是一種利用線性系統(tǒng)狀態(tài)方程，通過系統(tǒng)輸入輸出觀測數(shù)據(jù)，對系統(tǒng)狀態(tài)進行最優(yōu)估計的算法，本質(zhì)是利用兩個正態(tài)分布的融合仍是正態(tài)分布這一特性迭代最優(yōu)估計值。Kalman濾波雖然性能優(yōu)越且得到大量實踐證明其正確性，存在的不足是：由于無法確定測量過程中的系統(tǒng)噪聲和量測噪聲特性，只是試驗中給定了Q和R的噪聲參數(shù)，而Kalman濾波是基于精確數(shù)學模型遞推的過程。隨著測量時刻的不斷增加，根據(jù)遞推公式求得的P(k|k)會逐漸趨于0或者某一穩(wěn)態(tài)的常數(shù)，但是最優(yōu)觀測值X(k|k)與實際數(shù)據(jù)的差距越來越大，這種情況下，Kalman濾波器的預測和估計的功能逐漸喪失。而且由于濾波的實現(xiàn)平臺在計算式計算誤差的不斷累計傳遞也可能會使濾波器出現(xiàn)發(fā)散的現(xiàn)象[7]。同時，這類濾波算法參數(shù)調(diào)試復雜，而隨著傳感器設備老化，設備誤差會變化，相應參數(shù)也需要做出調(diào)整，需要長期跟蹤濾波效果，適時調(diào)整算法參數(shù)，增加了日常維護人員的工作負擔[8]。

基于以上，該文提出一種結(jié)合以上兩種手段的傳感器實時數(shù)據(jù)處理算法，通過高頻率采樣，將單傳感器單位時間內(nèi)多次采樣值看作為多傳感器的一次采樣值，用統(tǒng)計方法剔除虛假異常值后的觀測值，再利用Kalman濾波處理求得最優(yōu)估計值，較好地解決了虛假異常值產(chǎn)生的誤報問題。

1 算法描述

該文提出的算法，根據(jù)實踐中觀測值短期內(nèi)存在變化自相關(guān)性，而長期內(nèi)正態(tài)分布的特點[9]，將單位時間段內(nèi)傳感器監(jiān)測數(shù)據(jù)看作一個線性定常系統(tǒng)。假設在單位時間段內(nèi)所觀測的值是恒定的，可將單傳感器多次采樣值看作為多傳感器的一次采樣值，則該時間段內(nèi)傳感器的觀測誤差和噪音都是近似正態(tài)分布的，那么基于該值計算的移動平均值也是近似正態(tài)分布的，這樣的情況下是可以通過Kalman進行數(shù)據(jù)融合的。

首先，利用統(tǒng)計方法剔除虛假異常值影響，再用移動平均值降低隨機噪音的影響?；瑒哟翱诜绞将@取數(shù)據(jù)，設采樣傳感器數(shù)為N，當N>100時，采用拉伊達法即3西格瑪法，計算觀測值Xi與平均值X殘差的絕對值Vi，當Vi大于3倍標準差std即(|Xi-X|≥3*std)時，認為該觀測值為異常值，剔除所有異常值后計算剩余值的平均值gbmean。當N≤100時，采用格拉布斯法，根據(jù)n和顯著性水平a計算g(n,a)，將觀測值從小到大排序，計算最大值和最小值各自與平均值的殘差的絕對值Vi，當Vi>g(n,a)*std，判斷其為異常值剔除，并重新計算剩余觀測值的均值和std，重復以上步驟直至沒有出現(xiàn)異常值，而后計算剩余值的平均值gbmean[10-11]。格拉布斯法相較于拉伊達法更嚴謹準確，但由于需要排序和反復計算均值、標準差，當N較大時計算量較大，而當N>100時候，其結(jié)果和拉伊達法接近，為簡化計算可使用拉伊達法剔除異常值[12]。

然后，利用Kalman濾波計算最優(yōu)估計值。以gbmean作為經(jīng)驗值，移動平均值f作為觀測值，二值單位時間內(nèi)各自標準差作為其觀測誤差，而后根據(jù)上一次濾波時f值和gbmean值各自與t值的絕對誤差進行加權(quán)融合出新的誤差std1和std2，這樣的誤差結(jié)合傳感器自身最近誤差以及與最優(yōu)估計值誤差，當出現(xiàn)孤立異常值時更相信gbmean值，當出現(xiàn)連續(xù)異常值時更相信f值。

實時計算Kalman增益，根據(jù)一元卡爾曼濾波增益系數(shù)計算公式(1)：

(1)

計算出最優(yōu)估計值t。文中f值和gbmean值的誤差并非如在經(jīng)典Kalman濾波中是估計值，而是基于單位時間段值滑動數(shù)據(jù)計算而來的，避免估計中錯誤累積導致的離散問題。在計算出估計值后，結(jié)合歷史t值，基于SPC法，根據(jù)該時間段t值的采樣次數(shù)，利用格拉布斯法查表或拉伊達法確定控制限，進行異常值判斷。

之所以選擇f值和gbmean值進行融合，一是為了解決缺失值問題，由于各種因素，傳感器數(shù)據(jù)出現(xiàn)缺失值是十分常見的，移動均值可以比較好地解決常見線性系統(tǒng)的缺失值填充問題。二是為了更好地識別連續(xù)出現(xiàn)的異常值。如果直接使用帶有噪音的原始觀測值，f值雖然可以改善但仍無法完全剔除虛假異常值影響，很容易出現(xiàn)誤報；而僅使用gbmean值，由于會將第一個出現(xiàn)的真實異常值認為是虛假異常值，需連續(xù)出現(xiàn)的異常值影響總體方差后才能體現(xiàn)異常情況，故對于真實發(fā)生的異常反應具有滯后性。該文使用了Kalman法將f值和gbmean值進行融合，當異常值連續(xù)出現(xiàn)時能夠較快反映異常情況，在解決孤立虛假異常值的基礎上，實現(xiàn)對異常情況的較快反應。

2 算法流程

以圖1為例，說明此算法的計算流程。

圖1 算法流程

第一步：獲取Z時間內(nèi)的所有采樣值，先計算Z時間內(nèi)采樣數(shù)據(jù)的移動平均值f，以其標準差為其觀測誤差x1，上一次預測絕對誤差為x2。根據(jù)式(2)～式(4)[13]：

(2)

(3)

std1=w1*x1+w2*x2

(4)

計算其誤差std1。而后，如采樣次數(shù)N≤100時，根據(jù)設定的顯著性水平a，利用格拉布斯法剔除異常值后計算移動平均值gbmean；當N>100時，利用拉伊達法剔除異常值后計算移動平均值gbmean，以其標準差為其觀測誤差x3，上一次預測絕對誤差為x4，計算其誤差std2。

第二步：移動窗口取數(shù)據(jù)進行滾動計算，實時更新f值和gbmean值及其各自誤差，根據(jù)式(1)計算K值，進而根據(jù)Kalman校正式(5)：

t=gbmean+K*(f-gbmean)

(5)

求出最優(yōu)觀測值t，并計算Z時間段內(nèi)各采樣點t值的標準差std3。

第三步：計算控制限，當N>100時，Vi控制限為3*std3;當N≤100時，需根據(jù)格拉布斯法則計算g(n,a)，Vi控制限為g(n,a)*std3。超出控制限的為異常值。

整個算法中需要設定的參數(shù)只有N以及使用格拉布斯法時的顯著性水平a，也只需存在N個f值、N個gbmean值和N個t值共3N個值，所需存儲量和計算量較小，調(diào)整參數(shù)難度較低，可在物聯(lián)網(wǎng)終端或邊緣服務器進行數(shù)據(jù)處理，實現(xiàn)分布式部署。

3 仿真數(shù)據(jù)驗證

3.1 數(shù)據(jù)準備

該文所采用的數(shù)據(jù)源于石家莊站鐵路站房2013年7月至2014年7月數(shù)據(jù)中的M16BL-1傳感器數(shù)據(jù)data2，在anaconda3環(huán)境下進行數(shù)據(jù)分析。為了驗證算法的濾波性能，在已有數(shù)據(jù)基礎上，將采樣密度調(diào)整至每分鐘一次，結(jié)合原始值利用pandas中線性插值法填充調(diào)整后產(chǎn)生的缺失值，生成數(shù)據(jù)525 583條，填充后數(shù)據(jù)折線趨勢圖如圖2所示。

圖2 data2趨勢圖

而后先利用numpy.randon.normal函數(shù)生成均值為0，標準差為8的正態(tài)分布隨機噪聲，再利用numpy.randon.randint函數(shù)生成1 000個[50,100]和1 000個[-100,-50]的隨機異常值點，這些異常值隨機分布，有孤立有連續(xù)，data2值加上噪音和異常值點的數(shù)據(jù)為etl2值，也就是要處理的數(shù)據(jù)值，折線趨勢圖如圖3所示。

圖3 etl2趨勢圖

數(shù)據(jù)表格式如表1所示。

表1 處理前數(shù)據(jù)

當然，為了更加直觀地體現(xiàn)算法處理效果，增加的噪音和異常值較為明顯，實際異常值不會這么多。

3.2 算法實現(xiàn)

在算法實現(xiàn)上，為了快速實現(xiàn)，采用python語言進行開發(fā)，在anaconda3環(huán)境下，使用了pandas、outliters、math包。其中格拉布斯法是outliters包的smirnov_grubbs函數(shù)實現(xiàn)。移動平均值使用pandas.dataframe.rolling函數(shù)計算，由于采用數(shù)據(jù)的短期內(nèi)自相關(guān)而長期內(nèi)正態(tài)分布的特點，根據(jù)參考文獻[10]的統(tǒng)計結(jié)果，選取15分鐘作為Z時間段長短，設置滑動時間窗口為最近15分鐘的采樣數(shù)據(jù)，移動平均值為f15，即N值為15，數(shù)據(jù)以dataframe格式存儲為df1，顯著性水平為0.99，初始預測誤差為15，gbmean值和t值計算代碼如下：

Import pandas aspd

Import matplotlib.pyplot as plt

From sklean.metrics import mean_squared_error

From sklean.metrics import mean_absolute_error

Import numpy as np

Form outliers import smirnov_grubbs as grubbs

Import math

df1=pd.read_csv("含噪音數(shù)據(jù).csv")

df1['gbmean']=None

df1['t']=None

df1['vi1']=None

df1.loc[14,'vi1']=15

df1.loc[14,'vi2']=15

for i in range(15,len(df1)):

s=grubbs.test(df1.loc[i-14:i,'etl2'],0.99)

df1.loc[i,'gbmean']=s.mean()

std1=df1.loc[i-14:i,'f15'].std()

x1=df1.loc[i-1,'vi1']

x2=df1.loc[i-1,'vi2']

w1=cou(std1,x1)

w2=1-w1

std2=s.std()

w3=cou(std2,x2)

w4=1-w3

std1=std1*w1+w2*x1

std2=std2*w3+w4*x2

stdx2=math.pow(std2,2)/(math.pow(std2,2)+math.pow(std1,2))

f1.loc[i,'t']=df1.loc[i,'gbmean']+stdx2*(df1.loc[i,'f15']-df1.loc[i,'gbmean'])

df1.loc[i,'vi2']=abs(df1.loc[i,'t']-df1.loc[i,'gbmean'])

df1.loc[i,'vi1']=abs(df1.loc[i,'t']-df1.loc[i,'f15'])

計算出gbmean值、t值后，去除少量由于移動均值計算產(chǎn)生的缺失值，數(shù)據(jù)如表2所示。

表2 處理后數(shù)據(jù)

續(xù)表2

利用matplotlib包進行數(shù)據(jù)可視化，趨勢圖如圖4和圖5所示。

圖4 gbmean趨勢圖

圖5 t趨勢圖

計算出t值，每次計算最近15個t值的均值xt和方差stdt，當N=15,顯著性水平為0.99時，最優(yōu)估計值的殘差g(n,a)=g(15,9)=3.292，當|t-xt|>3.292*stdt時，判斷發(fā)生異常。

3.3 結(jié)果驗證

首先，利用matplotlib包進行數(shù)據(jù)可視化，比較將f15、gbmean、t和真值data2值放在一張圖里進行對比，定性比較處理效果，由于數(shù)據(jù)集較大，圖6是整體濾波效果。

圖6 濾波效果圖

然后用sklean.metrics包進行定量比較。計算etl2、gbmean、t與data2值的均方誤差(MSE)和平均絕對誤差(MAE)，MSE是最常用的回歸損失函數(shù)，計算方法是求預測值與真實值之間距離的平方和，MAE是目標值和預測值之差的絕對值之和的平均值。MSE會賦予異常點更大的權(quán)重，也就是說對預測誤差給予更大的權(quán)重，用MSE可以比較方法對異常值處理的好壞，而MAE值體現(xiàn)了誤差，顯示了方法處理的性能。源數(shù)據(jù)etl2的MSE約為119.45，f15值的MSE約為7.99，gbmean值的MSE約為8.00，t值的MSE約為6.98，t值的MSE相較于gbmean下降了12.73%；處理前數(shù)據(jù)etl2的MAE約為7.04，f15值的MAE約為2.11，gbmean值的MAE約為2.23，t值的MAE約為2.03，t值的MAE相較于gbmean下降了8.72%。

同時計算各采樣點f15值、t值和gbmean值相較于真值data2絕對殘差的標準差，f15值殘差標準差約為1.88，t值殘差標準差約為1.68，gbmean值殘差標準差約為1.74，t值殘差標準差相較于gbmean值殘差標準差降低了3.18%，算法效果較為穩(wěn)定。綜上可知，該算法性能和可靠性提升明顯。

3.4 對預測效果提升的驗證

為了進一步證明數(shù)據(jù)處理對預測結(jié)果的影響，該文以工業(yè)界常用的lightgbm算法為例，對處理前后的數(shù)據(jù)進行模型訓練和預測。LigthGBM是決策樹預測模型(GBDT)的一種，由微軟提供，它和XGBoost一樣是對GBDT的高效實現(xiàn)，原理上它和GBDT及XGBoost類似，都采用損失函數(shù)的負梯度作為當前決策樹的殘差近似值，去擬合新的決策樹?；贖istogram的決策樹算法，更低的內(nèi)存占用和更快的處理速度，基于OpenMP多線程加速和基于OpenCL的異構(gòu)加速，可以利用多種硬件加速學習過程，應用范圍更廣[14]。

首先，利用上文的處理算法對數(shù)據(jù)進行處理；然后，對于周期性較強的時間序列問題，可以通過將時間序列問題轉(zhuǎn)化為回歸問題進行預測，這里將時間這一特征進行特征工程處理，分解出表3的10個特征。

表3 訓練特征

然后，分別使用原始數(shù)據(jù)etl2和處理后數(shù)據(jù)xmean1進行l(wèi)ightgbm模型訓練，這里由于已經(jīng)是回歸問題，故采用隨機抽取的方式取75%的數(shù)據(jù)為訓練集，25%的數(shù)據(jù)為測試集，以MSE值作為訓練依據(jù)，以ETL2值訓練的模型MSE值為118.447，以xmean1值訓練的模型MSE值為7.762 15，數(shù)據(jù)處理后模型預測的MSE值下降了93.44%，性能提升明顯。

4 結(jié)束語

該算法結(jié)合了現(xiàn)有物聯(lián)網(wǎng)數(shù)據(jù)處理中的兩種處理方法，將多傳感器數(shù)據(jù)測量方法應用于單傳感器數(shù)據(jù)處理，兼顧了成本和準確性，可以解決實際工作中物聯(lián)網(wǎng)大量傳感器高采樣次數(shù)，傳感器誤差存在周期性變化的場景，實現(xiàn)動態(tài)剔除異常值、實時濾波和實時計算控制限。相較于傳統(tǒng)Kalman濾波法，調(diào)整參數(shù)少，維護成本低；相較于神經(jīng)網(wǎng)絡等機器學習方法，計算量和存儲數(shù)據(jù)量小，適合實時數(shù)據(jù)處理和分布式進行部署。