趙新華，范振東，何 宇，查益華

（1.國網(wǎng)新源水電有限公司新安江水力發(fā)電廠，杭州311608；2.中國電建集團華東勘測設(shè)計研究院有限公司，杭州311122）

0 引言

大壩監(jiān)測數(shù)據(jù)是評估運行期大壩安全性態(tài)最直觀、最可靠的方法之一［1，3］，而自動化監(jiān)測因其具有快速、實時、在線監(jiān)測的優(yōu)勢在各大水電站大壩監(jiān)測中得到廣泛應(yīng)用，但是自動化監(jiān)測由于采集頻次高且容易受雷擊、電壓等影響，容易出現(xiàn)異常值和毛刺點［4］。若異常值量值與正常值相差較大，則往往會覆蓋監(jiān)測效應(yīng)量的變化規(guī)律。此外自動化監(jiān)測數(shù)據(jù)中異常值出現(xiàn)的位置、量值具有不確定性，人工刪除異常值存在工作量大、主觀性強等不足。目前，常用的大壩監(jiān)測數(shù)據(jù)異常值識別方法有邏輯判別法、統(tǒng)計判別法以及基于模型的方法［5，7］等，其中，邏輯判別法需要給定一個測值的邏輯合理范圍，不同監(jiān)測效應(yīng)量的合理范圍不同，不具有普遍性。統(tǒng)計判別法適用于荷載條件變化不大的監(jiān)測數(shù)據(jù)中，單純使用統(tǒng)計判別法的誤差檢測率往往不高?；谀Ｐ偷姆椒ㄐ韪鶕?jù)環(huán)境量等信息建立相應(yīng)的預(yù)測模型，以模型預(yù)測值和實際測值的殘差為檢測對象，根據(jù)統(tǒng)計學方法或者機器學習方法進行異常值檢測，該方法是一種有監(jiān)督學習方法，對粗差等異常值檢測的精度受限于所建模型精度，而所建模型精度往往又受異常值的影響，因此用于建模的樣本數(shù)據(jù)往往需人工審核去除異常值，沒有達到減少人工工作量的效果。本文提出一種基于數(shù)據(jù)重構(gòu)與孤立森林法［8］（Iso?lation Forest）的異常數(shù)據(jù)檢測方法，這是一種無監(jiān)督的學習方法，適用于連續(xù)型的時間序列數(shù)據(jù)。

孤立森林算法是由劉飛博士在陳開明、周志華教授指導(dǎo)下提出的一種基于集成學習的快速異常檢測方法，算法因具有高精準度、無監(jiān)督等優(yōu)點，被廣泛應(yīng)用于大數(shù)據(jù)的異常值檢測［9，10］。由于孤立森林算法不適用于有趨勢性變化的數(shù)據(jù)序列，因此需要先對數(shù)據(jù)序列進行分解與重構(gòu)，將趨勢項分離出來，剩余項作為孤立森林訓練和異常檢測的樣本。為此，本文利用小波變換［11］、快速傅里葉變換［12］、奇異譜法［13］對監(jiān)測數(shù)據(jù)進行分解與重構(gòu)，以分離出趨勢項，需要說明的是，這三種趨勢性分離方法適用于時間連續(xù)且不存在趨勢性誤差的時間序列，此外，時間序列樣本越長，分離出的趨勢項精度越高。為進一步消除異常值和毛刺點，經(jīng)孤立森林法進行異常值剔除后，對其余測點采用拉依達準則（Pauta criterion）進行異常值檢測。

1 大壩自動化監(jiān)測數(shù)據(jù)分解與重構(gòu)

大壩監(jiān)測自動化數(shù)據(jù)異常檢測的前提是將監(jiān)測數(shù)據(jù)的趨勢項分離出來，為此，本文利用小波變換、傅里葉變換和奇異譜法對監(jiān)測數(shù)據(jù)進行分解與重構(gòu)。

1.1 小波變換

對于一維數(shù)據(jù)樣本x(n)(n= 0，1，…，N- 1)，采用小波進行分解，x(n)可分離出高頻信號和低頻信號。設(shè)P0=x，則第j層上的一維離散小波可表示為：

式中：Pj為低頻序列；Qj為高頻序列；cAj為低頻系數(shù)；cDj為高頻系數(shù)，（j= 1，2，…，L，k= 0，1，2，…，N/2j- 1)；L為小波分解的層數(shù)；φ為尺度函數(shù)；ψ為小波基函數(shù)。式（1）可形象表述如圖1所示。

圖1 小波分解示意圖Fig.1 Wavelet decomposition diagram

小波重構(gòu)算法與分解算法相反［14］，小波分解出的高頻序列Qj可視為噪音，經(jīng)過多層分解后，可分離出表征時間序列數(shù)據(jù)變化規(guī)律的主分量，本文采用第L層低頻近似部分和高頻細節(jié)部分模擬數(shù)據(jù)樣本的趨勢項xc(n)，即：

1.2 傅里葉變換

對于一維數(shù)據(jù)樣本x(n)(n= 0，1，…，N- 1)，式（3）稱為x(t)的傅里葉變換。

式中：X(k)稱為x(n)的第k層頻譜函數(shù)。

時間序列x(n)經(jīng)過傅里葉變換后失去時間特性，即X(k)僅具有頻率特性，X(k)由x(n)在時間序列上的特性決定。

本文以高斯函數(shù)濾波［15］為低通濾波器，以通過高斯濾波器后輸出的時域信號模擬數(shù)據(jù)樣本的趨勢項xc(n)。

1.3 奇異譜

奇異譜分析是由Broomhead 等［16］提出的一種適用于非線性、非平穩(wěn)信號的分解與重構(gòu)方法，具有計算工作量小的優(yōu)點［17］。

對于一維數(shù)據(jù)樣本x(n)(n= 0，1，…，N- 1)，假設(shè)嵌入的維數(shù)為m，時間延遲為τ，則數(shù)據(jù)樣本嵌入m×l維相空間為：

式中：k= 1，2，…，l；l=N-(m- 1)τ；X=[X1，X2，…，Xl]代表m×l維相空間的l個映射點。

令C為X的m×m維協(xié)方差矩陣，則：

通過對協(xié)方差矩陣C進行奇異值分解，可得到一組非負的奇異值（ei，i= 1，2，…，m），按e1≥e2≥…≥em≥0 排列組成奇異譜。其中大的奇異值對應(yīng)信號中的有效成分，小的奇異值可視為信號中的噪聲。

ek對應(yīng)的特征向量Ek稱為經(jīng)驗正交函數(shù)（EOF），第k個主分量定義為樣本序列x(n)在Ek上的正交投影系數(shù)：

當Ek和已知時，重構(gòu)的樣本信號如下：

本文采用第1主分量重構(gòu)的樣本信號模擬數(shù)據(jù)樣本的趨勢項xc(n)。

2 孤立森林算法

孤立森林（iForest）算法的提出是基于異常點的兩個特征［18］：分布稀疏以及離密度高的群體較遠，具有上述特征的點也被稱為容易被孤立的點。針對一組連續(xù)型的數(shù)據(jù)集，孤立森林算法的核心在于隨機進行采樣并構(gòu)造一定數(shù)量的iTree，由這些iTree組成一個iForest。構(gòu)造iForest的主要步驟如下：

（1）本文以分離趨勢項后的剩余項的絕對值即|x(n) -xc(n)|為訓練集，越靠近0 則該數(shù)據(jù)點是正常點的概率越大，為此建立一組用于訓練的樣本集g={g1，g2，…，g2N}，gi=-gi+N= |x(i) -xc(i)|，從樣本集中隨機選擇m個樣點作為子采樣集D={d1，d2，…，dm}，作為樹的根節(jié)點。

（2）從當前子采樣集中隨機一個分裂點p，p介于當前子采樣集中最小值到最大值之間。

（3）對子采樣集的每個數(shù)據(jù)di，若di

（4）重復(fù)步驟（2）、（3），不斷構(gòu)造新的左、右子樹，直至滿足下列條件之一：①D中只有一個數(shù)值或者N個一樣的數(shù)值，無法進一步劃分；②樹的高度達到限值。

（5）重復(fù)上述步驟，直至iTree 的數(shù)目達到限值，由這些iTree組成一個孤立森林。

對任何一個查詢數(shù)據(jù)x，通過循環(huán)一次孤立森林里所有樹，可以計算該查詢數(shù)據(jù)在每一顆iTree的路徑長度h(x)，進而可以計算在該孤立森林里路徑長度的期望E[h(x)]。對于提供含有m個樣本數(shù)量的子樣集D，由于iTree 的結(jié)構(gòu)與二叉樹等價，因此其搜索路徑長度平均值c(m)等價于二叉樹中失敗查詢的路徑長度［18］，文獻［19］給出了二叉樹中失敗查詢的路徑長度：

式中：H（m-1）為調(diào)函數(shù)，可以被估計為H（m-1）=ln（m-1）+γ，γ為歐拉常數(shù)；c（m）作用為標準化查詢數(shù)據(jù)x的路徑長h（x）。

查詢數(shù)據(jù)x的異常分數(shù)s計算公式如下：

查詢數(shù)據(jù)x根據(jù)如下準則進行異常檢測：當E[h(x)]→c(m)時，s→0.5，不能判斷查詢數(shù)據(jù)x是否為異常；當E[h(x)]→0 時，s→1 時，被檢測為異常點；當E[h(x)]→m- 1時，s→0時，被檢測為正常點。

拉依達準則通過計算標準偏差判斷測值偏離整體的程度，當偏離程度超出允許區(qū)間時，認為該測值為粗大誤差，需要注意的是使用拉依達準則需要足夠多的樣本，當樣本數(shù)較少時，其剔除粗大誤差的準確性會降低。為進一步消除粗大誤差，經(jīng)孤立森林法進行異常值剔除后，對其余測點，若剩余項小于指定值，由于剩余項為無趨勢性的平穩(wěn)時間序列，可通過拉依達準則方法對異常值進一步進行檢測。綜合上述理論，基于數(shù)據(jù)重構(gòu)和孤立森林法的大壩自動化監(jiān)測數(shù)據(jù)異常檢測的流程如圖2。

圖2 算法流程圖Fig.2 Algorithm flowchart

3 實例分析

為驗證算法的可行性，構(gòu)造一組含異常值的大壩變形時間序列，該序列從2015年1月1日至2018年12月31日，共有1 461 個數(shù)據(jù)測點，其中異常值91 個，占總測值的6.22%，見圖3。

分別用小波變換-孤立森林法、傅里葉變換-孤立森林法、奇異譜-孤立森林法對圖3 測值進行異常檢測，其中iTree 樣本容量為256，最大層高為8層，森林里共計有100顆iTree，最大迭代次數(shù)為50 次，異常分數(shù)判別閾值y選取將影響到算法的檢測率，若取值過大，則多數(shù)異常值不能檢測，反之，則誤判率較高，故本文按式（10）取值，閾值y先小后大。

圖3 含異常值時間序列Fig.3 Time series with outliers

式中：c1表示閾值初始值，實際取值時可先取0.50，根據(jù)實際情況進行調(diào)整；c2為參數(shù)，文中取1.0；n為當前迭代次數(shù)。

各算法檢測效果見表1、圖4 和圖5。由圖表可見，各算法的異常值檢測率均在95%以上，小波變換和傅里葉變換的檢測率大于奇異譜法。誤判率均在8%以下，奇異譜法的誤判率最低，誤判點主要出現(xiàn)在連續(xù)出現(xiàn)異常值的區(qū)域。由圖可見，對于具有趨勢性和周期性變化的時間序列測值，本文的三種方法均具有較高的檢測率，再經(jīng)拉依達準則判別后，基本能消除異常值帶來的影響。

圖4 各算法異常檢測結(jié)果Fig.4 Abnormal recognition curves of different algorithms

圖5 各算法與拉依達準則異常檢測結(jié)果Fig.5 Abnormal recognition curves of different algorithms and Pauta criterion

表1 各算法異常檢測成果表Tab.1 Abnormal recognition results of different algorithms

大壩自動化監(jiān)測數(shù)據(jù)異常值大小、數(shù)量以及分布規(guī)律要復(fù)雜得多，為進一步驗證本文算法的有效性，以兩個大壩工程自動化監(jiān)測數(shù)據(jù)進行異常檢測。其中A 水電站大壩位于福建閩江，最大壩高101 m，壩頂長783 m，對其繞壩滲流測點UP60 監(jiān)測數(shù)據(jù)進行異常值檢測。另一水電站大壩B 位于云南省南盤江，最大壩高96.5 m，壩頂長456.8 m，對壩體內(nèi)部沉降測點4 號自動化監(jiān)測數(shù)據(jù)進行異常值檢測，該數(shù)據(jù)無明顯異常點。

檢測成果見表2和圖6。A 水電站繞壩滲流測點UP60共有8 527 個數(shù)據(jù)，由于異常值存在，無法判斷繞壩滲流孔水位測值情況，剔除異常點后，該繞壩滲流孔測值在60～62.5 m 之間，無明顯趨勢性變化。B 水電站大壩內(nèi)部沉降變形測點4 號共有306 個數(shù)據(jù)，測點測值無明顯異常點，但由于測值初期，測值較少且變化較大，初期測值被誤判為異常點。通過數(shù)據(jù)重構(gòu)-孤立森林-拉依達準則法進行異常檢測后，基本能消除異常值帶來的影響且誤判率整體較低，能夠滿足工程實際應(yīng)用。

圖6 大壩自動化監(jiān)測測點異常檢測成果Fig.6 The results of abnormal detection of the dam automatic monitoring points

表2 典型測點異常檢測成果表Tab.2 The results of abnormal detection of the dam typical points

4 結(jié)論

本文利用小波變換、傅里葉變換和奇異譜法對時間序列數(shù)據(jù)進行分解與重構(gòu)，提取出趨勢項后通過孤立森里法對剩余項進行異常檢測，最后根據(jù)拉依達準則法進一步提出異常值。將組合數(shù)據(jù)重構(gòu)-孤立森林-拉依達準則法的組合算法用于大壩自動化監(jiān)測數(shù)據(jù)異常檢測中，得到以下結(jié)論。

（1）本文的組合方法是一種無監(jiān)督學習，除時間序列以外不需其他特征標簽，不需要對訓練樣本進行學習，擴大了其適用范圍，降低了工作量。

（2）通過數(shù)據(jù)重構(gòu)-孤立森林-拉依達準則法進行異常檢測后，基本能消除異常值帶來的影響且誤判率整體較低，能夠滿足工程實際應(yīng)用。

（3）誤判點主要出現(xiàn)在連續(xù)異常值出現(xiàn)區(qū)域的相鄰正常測點，這主要是因為小波變換、傅里葉變換和奇異譜法重構(gòu)分離出的趨勢項受連續(xù)異常值的影響，因此對于時間不連續(xù)、測值少且變化劇烈的時間序列適用性不強?！?/p>