高熠飛,王建平,李林峰

(1.河海大學水文水資源學院，江蘇 南京 210098；2.南瑞集團(國網(wǎng)電力科學研究院)有限公司，江蘇 南京 211000)

時間序列泛指隨時間或空間有序變化的數(shù)據(jù)集合。水文行業(yè)中的流量、水位數(shù)據(jù)都是典型的時間序列數(shù)據(jù)。在水文時間序列中，存在2種情況：(a)數(shù)據(jù)受到新機制的作用,如突降大暴雨或是閘門啟閉造成數(shù)據(jù)點的跳變；(b)人為差錯或測量儀器故障導致錯數(shù)或漏數(shù)，以至于數(shù)據(jù)測量收集產(chǎn)生誤差。經(jīng)常會產(chǎn)生與其他數(shù)據(jù)在某種相似度量上存在預定義或顯著差異的數(shù)據(jù)[1]，一般將這種數(shù)據(jù)稱為水文序列中的異常值。

異常值對時間序列的模型識別、參數(shù)估計、診斷檢驗乃至預測都有重要的影響[2]，當前水情測報與水庫調度自動化系統(tǒng)在工作過程中需要大量調用并分析歷史時間序列，異常值的存在直接影響系統(tǒng)相關決策的正確性。此外，水文序列中的異常值可能潛藏著更有價值和意義的水文信息和知識[3]。因此，正確檢測出時間序列中的異常值具有十分重要的現(xiàn)實意義，可以減少異常值對數(shù)據(jù)分析的干擾，提高預報等活動的準確性，另外，還可以幫助分析異常值背后蘊藏的信息。

目前國內(nèi)外對于異常檢測的研究可以分為以下幾大類：

a. 基于預測模型的檢測方法?；陬A測的方法通過構建模型計算待檢測點的預測值，基于預測值與實測值的偏離程度進行異常值判別?；诰€性模型如自回歸模型、滑動平均模型[3-5]等的算法計算快捷，但線性模型適用于平穩(wěn)數(shù)據(jù)，對于非平穩(wěn)的水文數(shù)據(jù)預測準確率不高。非線性模型方法例如基于神經(jīng)網(wǎng)絡模型[6]或是支持向量機模型[7]的檢測算法，對水文數(shù)據(jù)的預測較準確，但在實際工程應用中會面臨閾值設置以及隨著檢測進行，閾值可能發(fā)生變化等問題，且非線性模型計算量大，耗費資源。

b. 基于距離的方法。通過設定距離函數(shù)，計算數(shù)據(jù)空間中數(shù)據(jù)點之間的距離，當某一個數(shù)據(jù)與其他數(shù)據(jù)距離較大時，認為該點數(shù)據(jù)異常[8]。付培國等[9]提出一種基于密度偏倚抽樣的異常檢測方法，首先利用基于密度偏倚的概率抽樣對待測數(shù)據(jù)進行概率抽樣，之后對抽樣數(shù)據(jù)計算局部異常系數(shù)，局部異常系數(shù)越大的數(shù)據(jù)越有可能異常。該方法改善了傳統(tǒng)基于距離的算法時間復雜度高的情況，處理大數(shù)據(jù)集時效果明顯，但前期需要設置大量參數(shù)，相當依賴使用者的經(jīng)驗，且該方法適合對歷史序列進行異常值挖掘，不適合作為實時異常值檢測方法使用。

c. 基于聚類的方法。將數(shù)據(jù)集分成若干類，不屬于任何類的或某個類中數(shù)據(jù)量特別少的點被判定為異常點[10]。此類方法面臨初始類采用隨機選取帶來的局部最優(yōu)問題。左進等[11]提出一種改進k均值聚類的異常檢測算法，通過計算所有數(shù)據(jù)點的緊密性，在數(shù)據(jù)緊密的地方均勻選擇k個初始中心，避免了陷入局部最優(yōu)，對于異常點的檢測效率具有明顯的提升，但基于聚類的方法依賴初始類個數(shù)的設置以及數(shù)據(jù)集中存在異常點的假定，且基于聚類的方法本質在于分類，用來進行異常點的診斷不夠高效。

本文的異常值檢測算法是針對水情測報與水庫調度自動化系統(tǒng)實際運行中對水文流量數(shù)據(jù)的質量需求而研究的，因此算法不能占用主系統(tǒng)太多的資源。復雜的算法例如神經(jīng)網(wǎng)絡以及支持向量機，對單個點的準確性可能更高，但隨著數(shù)據(jù)量的增大，算法的時間復雜度較高，對資源的消耗較大，不適合水情測報與水庫調度自動化系統(tǒng)實時調用。

實際采集得到的流量數(shù)據(jù)由于風浪、船舶、發(fā)電負荷調整、泄洪閘門啟閉以及其他人類活動等影響，產(chǎn)生了較多的噪音數(shù)據(jù)，因此還需要算法具有較強的抗干擾能力?；诰€性模型預測的算法計算簡便，但流量數(shù)據(jù)局部波動以及水文數(shù)據(jù)非平穩(wěn)的特性使得模型預測結果的準確性難以保證。

綜上，本文考慮利用假設檢驗的統(tǒng)計學方法進行異常檢測，對于數(shù)據(jù)的局部波動從統(tǒng)計學層面進行判斷。統(tǒng)計學中常用的基于高斯分布的判斷方法，其參數(shù)均值以及標準差在序列中存在極端異常值時會產(chǎn)生較大偏差，不夠穩(wěn)定[12]，因此考慮選用魯棒性更好的中位數(shù)[13]、中位數(shù)絕對偏差來替代均值和方差，利用待測點前向滑動窗口中的中位數(shù)以及中位數(shù)絕對偏差作為參數(shù)回歸柯西分布，計算待測點的置信區(qū)間，根據(jù)待測點是否落在置信區(qū)間來判斷待測點是否異常。該方法利用待測點歷史數(shù)據(jù)的信息，可以自適應調整判斷區(qū)間，并且算法受窗口中可能存在異常值的影響小，具有較好的魯棒性，計算量小，適用于對大量水文數(shù)據(jù)的實時自動化檢測。

1 算 法 介 紹

1.1 相關定義

定義1水文時間序列。由觀測值及記錄時間組成的有序集合X={(x1,t1),(x2,t2),…,(xn,tn)}，令di=(xi,ti) 表示包含時刻以及觀測值的水文對象，其中，t表示當前水文對象的時間，x表示對應時間的值，比如水位值。

定義2滑動窗口?；瑒哟翱诔７譃閱芜呉约半p邊。雙邊窗口更適合對待檢測節(jié)點前驅和后繼窗口數(shù)據(jù)都已知的歷史數(shù)據(jù)進行挖掘，而單邊窗口則是忽略待測點以后的數(shù)據(jù)，只利用前驅窗口中的數(shù)據(jù)信息進行分析，而本文的期望是在工程應用中實時判斷數(shù)據(jù)異常與否，因此，單邊窗口明顯更適合。因此對于待檢測點di，假定窗口寬度為w個時間點，滑動窗口即是待測點di的前k時刻觀測值構成的集合，即

ηi={xi-k,…,xi-2,xi-1}

(1)

定義3中位數(shù)絕對偏差。描述序列各點與序列中位數(shù)的偏離程度。假設式(1)中ηi序列的中位數(shù)為mi，則計算出序列中每個點與中位數(shù)mi的殘差，生成殘差序列{|xi-k-mi|,…,|xi-2-mi|,|xi-1-mi|}，計算該殘差序列的中位數(shù)即原序列的中位數(shù)絕對偏差。

定義4水文時間序列異常。根據(jù)待定點di所對應的滑動窗口ηi中數(shù)據(jù)計算出窗口內(nèi)數(shù)據(jù)的中位數(shù)mi以及中位數(shù)絕對偏差δi，代入式(2)，計算ti時刻對應的累積分布函數(shù)值，記為C(xi)，若C(xi)在置信區(qū)間之外，則判定xi為異常點.其中C(xi)計算公式如下:

(2)

定義5基于中位數(shù)與中位數(shù)絕對偏差的柯西分布。柯西分布，又稱柯西-洛倫茲分布，其概率密度函數(shù)為

(3)

式中：x0——分布峰值位置的位置參數(shù)，在本文中用中位數(shù)進行估計；γ——最大值1/2處1/2寬度的尺度參數(shù)，本文用中位數(shù)絕對偏差進行估計。

1.2 方法對比及選擇

為了直觀地看出基于高斯分布的算法在遇到極端異常值干擾時的情況，選取一段存在6個異常值的2 h流量數(shù)據(jù)進行試驗。

圖1是一段包含異常點的水文流量時間序列集X={d1,d2,…,dn}的流量-時間關系曲線，可以看到，流量數(shù)據(jù)中有較明顯的6個偏差很大的異常點。

圖1 包含異常點的流量數(shù)據(jù)Fig.1 Hydrological flow data containing outliers

設定滑動窗口η寬度為10個時間點。對于每一時刻的待測點di都有其對應的滑動窗口ηi。對于每個待測點對應的滑動窗口內(nèi)的數(shù)據(jù)進行分析，分別繪制其Q-Q圖(quantile-quantile plot)，其橫軸為水文流量值的分位值，縱軸為標準正態(tài)分布的分位值。若樣本數(shù)據(jù)近似于正態(tài)分布，Q-Q圖上的點會排列在直線y=σx+μ附近，其中σ表示該組數(shù)據(jù)的方差，μ表示均值。點的散布與直線越接近，就越接近正態(tài)分布[14]。

雖然根據(jù)隨機選取部分滑動窗口內(nèi)數(shù)據(jù)的Q-Q圖實例，可以看出，窗口寬度較小時，窗口內(nèi)數(shù)據(jù)可基本看成符合正態(tài)分布。

分別計算窗口內(nèi)數(shù)據(jù)的平均值、中位數(shù)、標準差以及中位數(shù)絕對偏差，得到序列{xmean}、{xmedian}、{xSD}、{xMAD}。其中，均值序列與中位數(shù)序列的對比見圖2，標準差序列與中位數(shù)絕對偏差序列的對比見圖3。

圖2 滑動窗口內(nèi)數(shù)據(jù)的均值與中位數(shù)對比Fig.2 Comparison of mean and median data in sliding windows

圖3 滑動窗口內(nèi)數(shù)據(jù)的標準差與中位數(shù)絕對偏差對比Fig.3 Comparison of standard deviation and median absolute deviation of data in sliding windows

由圖2、圖3可以明顯看出，當滑動窗口內(nèi)包含異常值時，均值與標準差均會發(fā)生窗口寬度的跳變，然而同時刻的中位數(shù)以及中位數(shù)絕對偏差卻幾乎不受窗口中異常值的影響，顯得很穩(wěn)定。

基于正態(tài)分布的檢測算法在窗口中出現(xiàn)異常值，尤其是偏差較大的異常值時，其參數(shù)均值、標準差會受到很大的影響，極易造成誤判，魯棒性差。所以每次檢測到異常值后都需要及時進行修正來減少異常值的影響，然后利用修正值進行下一個時間點的檢測。修正一般分為人工修正或是自動修正[3]。基于人工的修正對操作人員的專業(yè)知識要求較高，且實時性難以保證。自動修正一般是基于預測，例如通過前向k-近鄰窗口內(nèi)的數(shù)據(jù)得出預測值[3]，以此代替異常值，但以水情測報與水庫調度自動化系統(tǒng)中最常用的水庫出庫流量數(shù)據(jù)為例，數(shù)據(jù)局部波動很大，這種情況下預測模型的準確性難以保證，因此對后續(xù)判別的準確性也會帶來影響。因此，期望能夠在檢測到異常值的情況下，即使異常值偏差很大，也可以在不對異常值做任何操作的情況下繼續(xù)準確地進行下一個點的檢測。

觀察到中位數(shù)以及中位數(shù)絕對偏差受到窗口中極端異常值的影響幾乎可以忽略，有很好的主體數(shù)據(jù)(或稱大部分數(shù)據(jù))的特性[15]。

利用中位數(shù)判別異常值的方法早有一定的研究[16]，主要是直接考慮中位數(shù)與待測數(shù)據(jù)之間的差異來判斷是否異常。但此種方法忽略了窗口中數(shù)據(jù)的局部波動情況，會使一些在合理波動內(nèi)的數(shù)值被誤判成異常值。因此，考慮在原方法的基礎上加入反映數(shù)據(jù)變化尺度的參數(shù)來進行綜合判斷。本文在中位數(shù)的基礎上加入中位數(shù)絕對偏差來進行綜合判斷。

鑒于柯西分布與正態(tài)分布的相似性，窗口內(nèi)數(shù)據(jù)也可近似看作符合柯西分布。有研究證明，用樣本的中位數(shù)來對柯西分布進行參數(shù)估計是一種相當優(yōu)良的估計[17]。有學者使用中位數(shù)以及中位數(shù)絕對偏差作為參數(shù)回歸柯西分布并進行了網(wǎng)絡流量數(shù)據(jù)的異常檢測，取得了很好的效果。因此本文利用滑動窗口內(nèi)數(shù)據(jù)的中位數(shù)以及中位數(shù)絕對偏差來回歸柯西分布進行水文序列的異常值檢測，這樣既能考慮實際水文流量數(shù)據(jù)的局部波動情況，還大幅度減小了檢測出的異常數(shù)據(jù)對后續(xù)檢測的影響。

1.3 算法

1.3.1 算法步驟

本文提出的基于柯西分布的流量時間序列異常值檢測方法的核心思想是：讀取水文流量時間序列；設置滑動窗口的寬度為N；將di的前N時刻的數(shù)據(jù)的觀測值{xi-N,xi-N+1,…,xi-1}放入窗口中；計算出前向窗口內(nèi)數(shù)據(jù)的中位數(shù)、中位數(shù)絕對偏差；在預設置信度p的前提下將待檢測點觀測值xi代入基于柯西分布的異常值判斷函數(shù)E(x)進行判斷，對數(shù)據(jù)狀態(tài)進行標記(正?；蚴钱惓?；將窗口向后滑動一個時間點并進行下一個數(shù)據(jù)的判斷。

1.3.2 標準化函數(shù)

(4)

式中：m——待測點對應滑動窗口內(nèi)數(shù)據(jù)的中位數(shù)。

1.3.3 異常值判斷函數(shù)

假定當前數(shù)據(jù)di的觀測值xi、對應的滑動窗口ηi及其中位數(shù)xmediani、中位數(shù)絕對偏差xMADi均已得到。di與滑動窗口ηi內(nèi)的數(shù)據(jù)可近似認為同分布。

根據(jù)式(4)首先對待測數(shù)據(jù)用標準化函數(shù)N(x)進行處理，設置信度為p，則異常值判斷函數(shù)E(x)定義如下(其中E(x)=1代表通過檢測，E(x)=0代表數(shù)據(jù)存在異常)：

(5)

1.3.4 參數(shù)選擇

在讀取數(shù)據(jù)后，對di進行異常檢測的第一步就是確定滑動窗口ηi的寬度N。若窗口寬度過寬，則窗口內(nèi)的數(shù)據(jù)不一定符合柯西分布的假定。又因為基于柯西分布的方法依賴窗口中數(shù)據(jù)的中位數(shù)以及中位數(shù)絕對偏差，若窗口寬度太窄，則窗口內(nèi)數(shù)據(jù)的代表性不夠，會導致判斷出現(xiàn)誤差。

判斷函數(shù)E(x)中，置信度p的影響很大。p太大，會出現(xiàn)漏檢現(xiàn)象；反之p若是一味求小，則會出現(xiàn)誤檢測，對于水情測報與水庫調度自動化系統(tǒng)的正常工作會造成很大的影響。

定義檢出率αDR、誤檢率αFR。分別表示正確檢測出的異常值占所有異常值的比例、被誤檢測為異常值的正常值占所有正常值的比例。

對于待檢測的水文序列，檢測前要確定2個參數(shù)，即選擇滑動窗口的寬度N以及置信度p。窗口寬度太小，窗口內(nèi)數(shù)據(jù)缺乏代表性。但窗口寬度太大，窗口內(nèi)數(shù)據(jù)就無法服從柯西分布的假設。p越大，判斷函數(shù)允許范圍越小，但可能會造成誤檢；p越小，發(fā)生誤檢的概率越低，但發(fā)生漏檢的概率也隨之提升。因此需要從檢出率以及誤檢率來綜合考慮。首先選取一段人工標記過異常值的歷史序列，對算法的參數(shù)選擇實際是選擇使得檢出率與誤檢率的比值αDR/αFR最大的一組N、p組合。因此從(N,p)=(1,75%)開始，首先固定p，考慮N,設置增量為1,計算每個(N,75%)對應的模型在歷史數(shù)據(jù)集上的αDR/αFR值，當比值連續(xù)出現(xiàn)降低時停止增加N。接著設置p的增量為1%，計算每一個p對應每個N模型的αDR/αFR值直到p=100%，取其中比值最大的一組N、p作為后續(xù)檢測的參數(shù)。由于算法檢測到異常值后無需進行修改數(shù)據(jù)的操作即可繼續(xù)進行檢測，單次計算量小，且滑動窗口內(nèi)的數(shù)據(jù)較多時，窗口內(nèi)數(shù)據(jù)不再服從柯西分布，檢出率開始下降，因此參數(shù)選擇使用窮舉法能在實現(xiàn)參數(shù)自動化設置的同時不會帶來很大的算法時間復雜度。

2 實 例 分 析

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

選擇三峽出庫流量數(shù)據(jù)來進行驗證。采用2012-05-18至2012-07-23的實測小時流量數(shù)據(jù)。為了更直觀地顯示，將具體時間(日期)用時間點表示。

從圖4可以看出，流量數(shù)據(jù)的局部波動劇烈，異常值信息可能隱藏在波動中。

圖4 三峽原始小時出庫流量數(shù)據(jù)集Fig.4 Originally and hourly outbound flow data set of Three Gorges

2.2 結果分析

不同參數(shù)下算法的檢測結果如表1所示。

表1 不同參數(shù)下的檢測結果

評估結果表明，本文的算法能夠對流量數(shù)據(jù)中的異常點進行較準確的判別，盡管出庫流量數(shù)據(jù)局部波動劇烈，但在窗口寬度以及置信度選取合理的情況下，檢出率以及誤檢率均比較理想，說明算法具有可行性。并且發(fā)現(xiàn)，當檢測到偏差較大的異常值時，算法仍然能對接下來的數(shù)據(jù)進行正確檢測，證明本算法具有較好的魯棒性，可以保證實際工程應用中檢測系統(tǒng)的平穩(wěn)運行。同時本算法計算量小，隨著數(shù)據(jù)量的增大，算法時間復雜度線性增加，面對海量數(shù)據(jù)也能進行快速檢測。因此，本算法很適合水情測報與水庫調度自動化系統(tǒng)中實時快速檢測的應用場景需求。

表2 不同算法在本文數(shù)據(jù)集上的檢測結果對比

圖5是窗口寬度為15個時間點、置信度為95%時本算法的異常檢測結果，大部分點被判定為正常，也有部分點在置信區(qū)間以外(根據(jù)本算法判斷函數(shù)，即判定為異常點)。

圖5 (N, p)=(15, 95%)條件下本算法的檢測結果Fig.5 Results of the algorithm when (N, p) = (15, 95%)

本算法以前向窗口中的中位數(shù)以及中位數(shù)絕對偏差作為判斷參數(shù)，無法準確描述窗口中數(shù)據(jù)的趨勢特征，雖然使用滑動窗口一定程度上減弱了趨勢項的影響，但趨勢的存在還是造成了一些漏檢、誤檢。

2.3 算法對比與分析

為了分析本文所提算法相較于傳統(tǒng)的基于正態(tài)分布的算法在水文數(shù)據(jù)檢測下的優(yōu)勢，將本文算法與常用的基于自回歸模型的算法[4]以及基于滑動窗口預測[3]的算法進行對比，在本文數(shù)據(jù)集上的檢測效果如表2所示。

基于自回歸模型的算法在擬合回歸時采用相鄰幾個時間點的數(shù)據(jù)，由于流量數(shù)據(jù)局部波動較大，模型擬合的準確性難以保證，在本文數(shù)據(jù)集上的試驗也證實了這一現(xiàn)象，且該算法基于高斯分布判斷異常，異常點的出現(xiàn)會使臨近時間點的均值方差發(fā)生偏差，影響到對異常值后續(xù)時刻的數(shù)據(jù)檢測，因此該算法在本文數(shù)據(jù)集檢出率較低且出現(xiàn)了較多的誤檢?；诨瑒哟翱陬A測的算法檢測到異常值后要將異常值進行修正，并且利用修正值繼續(xù)進行異常值檢測，根據(jù)試驗發(fā)現(xiàn)，線性模型對這種劇烈波動的數(shù)據(jù)預測適用性較差，導致檢測出現(xiàn)了誤判，發(fā)生誤判后，模型又使用預測值錯誤地修改了誤判點的數(shù)據(jù)，導致后續(xù)的檢測連續(xù)出現(xiàn)誤檢，在本文數(shù)據(jù)集的檢測效果很不理想。

最后繪制這3種算法在本文數(shù)據(jù)集上的操作特性曲線[19]。判斷某個算法的優(yōu)劣主要是看曲線下方面積的大小，面積越大檢測效果越好。本文算法和其他算法在本文數(shù)據(jù)集上檢測效果的操作特性曲線對比如圖6所示。

圖6 不同方法在本文數(shù)據(jù)集下的檢測結果ROC對比Fig.6 Comparison of ROC results of different algorithms in the data set of this study

由圖6可以看出，本文所提算法，在局部波動較強烈的流量數(shù)據(jù)集上的檢測效果明顯優(yōu)于另外2種算法，且更加穩(wěn)定，表明該算法可以應對這種場景下水情測報與水庫調度自動化系統(tǒng)對異常值檢測的要求。

3 結 語

針對實測水文序列中數(shù)據(jù)存在異常跳變點的問題，基于滑動窗口的異常檢測算法，提出了基于柯西分布的方法來檢測流量時間序列中的點異常。采用中位數(shù)以及中位數(shù)絕對偏差來代替均值以及標準差進行回歸。即使數(shù)據(jù)局部波動劇烈或是滑動窗口中出現(xiàn)極端異常值，模型都保持穩(wěn)定，且檢測到異常值后無需對異常值進行處理即可繼續(xù)對下一個值進行檢測。采用三峽水庫小時出庫流量數(shù)據(jù)進行試驗，并通過不同算法的檢測結果進行對比，驗證了本文算法在這種數(shù)據(jù)條件下的可行性以及魯棒性。未來需要考慮對滑動窗口內(nèi)的數(shù)據(jù)進行趨勢影響消除，降低水文數(shù)據(jù)的趨勢對算法準確性的影響。