陳艷春，達(dá)鈺鵬

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

0 引 言

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

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

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

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

1 算法描述

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

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

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

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

(1)

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

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

2 算法流程

以圖1為例，說(shuō)明此算法的計(jì)算流程。

圖1 算法流程

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

(2)

(3)

std1=w1*x1+w2*x2

(4)

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

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

t=gbmean+K*(f-gbmean)

(5)

求出最優(yōu)觀測(cè)值t，并計(jì)算Z時(shí)間段內(nèi)各采樣點(diǎn)t值的標(biāo)準(zhǔn)差std3。

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

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

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

3.1 數(shù)據(jù)準(zhǔn)備

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

圖2 data2趨勢(shì)圖

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

圖3 etl2趨勢(shì)圖

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

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

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

3.2 算法實(shí)現(xiàn)

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

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'])

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

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

續(xù)表2

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

圖4 gbmean趨勢(shì)圖

圖5 t趨勢(shì)圖

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

3.3 結(jié)果驗(yàn)證

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

圖6 濾波效果圖

然后用sklean.metrics包進(jìn)行定量比較。計(jì)算etl2、gbmean、t與data2值的均方誤差(MSE)和平均絕對(duì)誤差(MAE)，MSE是最常用的回歸損失函數(shù)，計(jì)算方法是求預(yù)測(cè)值與真實(shí)值之間距離的平方和，MAE是目標(biāo)值和預(yù)測(cè)值之差的絕對(duì)值之和的平均值。MSE會(huì)賦予異常點(diǎn)更大的權(quán)重，也就是說(shuō)對(duì)預(yù)測(cè)誤差給予更大的權(quán)重，用MSE可以比較方法對(duì)異常值處理的好壞，而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%。

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

3.4 對(duì)預(yù)測(cè)效果提升的驗(yàn)證

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

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

表3 訓(xùn)練特征

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

4 結(jié)束語(yǔ)

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