基于改進Kalman的傳感器數據加權處理算法

陳艷春，達鈺鵬

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

0 引 言

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

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

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

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

1 算法描述

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

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

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

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

(1)

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

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

2 算法流程

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

圖1 算法流程

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

(2)

(3)

std1=w1*x1+w2*x2

(4)

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

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

t=gbmean+K*(f-gbmean)

(5)

求出最優觀測值t，并計算Z時間段內各采樣點t值的標準差std3。

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

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

3 仿真數據驗證

3.1 數據準備

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

圖2 data2趨勢圖

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

圖3 etl2趨勢圖

數據表格式如表1所示。

表1 處理前數據

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

3.2 算法實現

在算法實現上，為了快速實現，采用python語言進行開發，在anaconda3環境下，使用了pandas、outliters、math包。其中格拉布斯法是outliters包的smirnov_grubbs函數實現。移動平均值使用pandas.dataframe.rolling函數計算，由于采用數據的短期內自相關而長期內正態分布的特點，根據參考文獻[10]的統計結果，選取15分鐘作為Z時間段長短，設置滑動時間窗口為最近15分鐘的采樣數據，移動平均值為f15，即N值為15，數據以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("含噪音數據.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值后，去除少量由于移動均值計算產生的缺失值，數據如表2所示。

表2 處理后數據

續表2

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

圖4 gbmean趨勢圖

圖5 t趨勢圖

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

3.3 結果驗證

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

圖6 濾波效果圖

然后用sklean.metrics包進行定量比較。計算etl2、gbmean、t與data2值的均方誤差(MSE)和平均絕對誤差(MAE)，MSE是最常用的回歸損失函數，計算方法是求預測值與真實值之間距離的平方和，MAE是目標值和預測值之差的絕對值之和的平均值。MSE會賦予異常點更大的權重，也就是說對預測誤差給予更大的權重，用MSE可以比較方法對異常值處理的好壞，而MAE值體現了誤差，顯示了方法處理的性能。源數據etl2的MSE約為119.45，f15值的MSE約為7.99，gbmean值的MSE約為8.00，t值的MSE約為6.98，t值的MSE相較于gbmean下降了12.73%；處理前數據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%，算法效果較為穩定。綜上可知，該算法性能和可靠性提升明顯。

3.4 對預測效果提升的驗證

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

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

表3 訓練特征

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

4 結束語

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