王蘇杭,王一棠,龐 勇,李一陽,宋學(xué)官

(大連理工大學(xué)機(jī)械工程學(xué)院,高性能精密制造全國重點實驗室,遼寧大連 116024)

0 引言

城市的地下空間被用于滿足日益增長的交通、能源和電信的基建需求,而隧道是常見的地下空間利用方式之一[1-2]。盾構(gòu)機(jī)是一種集光電、機(jī)械、液壓、計算機(jī)技術(shù)為一體的數(shù)字化大型工程機(jī)械裝備[3-4],是隧道掘進(jìn)盾構(gòu)施工的核心,具有高效、安全和環(huán)保等優(yōu)點,已廣泛應(yīng)用于國內(nèi)多個重點建設(shè)項目。為保障盾構(gòu)機(jī)安全、高效、平穩(wěn)運(yùn)行,其智能運(yùn)維系統(tǒng)將施工現(xiàn)場采集的實測掘進(jìn)數(shù)據(jù)進(jìn)行挖掘和分析,從而對盾構(gòu)機(jī)的性能進(jìn)行精準(zhǔn)預(yù)測。然而,由于施工現(xiàn)場環(huán)境復(fù)雜,實際采集的數(shù)據(jù)會存在異常值情況,嚴(yán)重降低數(shù)據(jù)的利用價值。

在隧道掘進(jìn)過程中,部署在盾構(gòu)機(jī)上的傳感器采集信號,通過數(shù)字轉(zhuǎn)換和相應(yīng)的通信技術(shù),存儲到計算機(jī)或者服務(wù)器,即可得到盾構(gòu)機(jī)的實測數(shù)據(jù)。傳感器按照一定的采樣頻率采集信號,所以盾構(gòu)機(jī)實測數(shù)據(jù)是嚴(yán)格的時間序列數(shù)據(jù)。時間序列數(shù)據(jù)的異常形式一般分為3類:點異常、模式異常和序列異常[5-6]。在盾構(gòu)機(jī)實測數(shù)據(jù)中,點異常出現(xiàn)較多,而模式異常和序列異常較少見,因此點異常是研究的重點。

當(dāng)前,面向時間序列異常值檢測算法模型一般有基于預(yù)測[7]、基于統(tǒng)計[8-9]、基于機(jī)器學(xué)習(xí)和深度學(xué)習(xí)[10]、基于窗口等方法[11-12]。基于窗口的方法通過設(shè)計滑動窗口的大小對時間序列進(jìn)行分割,對各子時間序列進(jìn)行特征提取,在各子時間序列中定位異常值點。與其他方法相比,該類方法易于理解且計算簡單,降低了時間復(fù)雜度,并且能將動態(tài)的數(shù)據(jù)量靜態(tài)化[13],因此得到廣泛的應(yīng)用。戴慧等[14]提出了一種基于滑動窗口的由粗到細(xì)的兩階段探測策略,采用子時間序列的均值和方差作為子時間序列的特征信息,實現(xiàn)了對某變電站各設(shè)備表采集數(shù)據(jù)的異常值檢測。高熠飛等[15]計算窗口內(nèi)子時間序列的中位數(shù)以及中位數(shù)絕對偏差作為子時間序列的特征信息,較好地消除了滑動窗口內(nèi)極端異常值對算法的影響,實現(xiàn)了對水文序列異常值的檢測。唐向紅等[16]計算子時間序列內(nèi)相鄰兩點的斜率與設(shè)定的斜率范圍進(jìn)行比較,實現(xiàn)了對電流信號的檢測。基于數(shù)據(jù)斜率的方式能較準(zhǔn)確描述子時間序列的結(jié)構(gòu)特征,但工程數(shù)據(jù)中異常數(shù)據(jù)分布隨機(jī)性大、波動性大,導(dǎo)致傳統(tǒng)的基于斜率的檢測方法對數(shù)據(jù)形態(tài)信息利用不充分,而置信區(qū)間能夠以一定的可靠程度估計總體數(shù)據(jù)所在的區(qū)間。田騰等[17]采用基于子時間序列斜率的置信區(qū)間的方式進(jìn)行特征提取,對盾構(gòu)機(jī)掘進(jìn)速度的異常值檢測取得了良好的效果。陸秋琴等[18]同樣采用斜率的置信區(qū)間距離半徑對污染物子序列特征進(jìn)行提取,實現(xiàn)了對傳感器數(shù)據(jù)采集中異常值的檢測。但其對子時間序列斜率的置信區(qū)間半徑計算上,采用了傳統(tǒng)方式,運(yùn)算效率有待改進(jìn)。因此,本文采用一種基于滑動窗口提取子時間序列特征的異常值檢測方法,采用快速計算的方式得到子時間序列斜率的置信距離半徑作為特征,能夠有效利用序列特征進(jìn)行異常識別,提高了運(yùn)算效率。同時對檢測出的異常值,提出了基于多子時間序列的填補(bǔ)方法。

1 盾構(gòu)機(jī)實測數(shù)據(jù)點異常處理方法

所提方法分為異常值檢測和修正2部分。其中異常值檢測部分有2個階段,第1階段是異常子時間序列檢測,采用滑動窗口的方法將未處理的盾構(gòu)機(jī)實測數(shù)據(jù)劃分成若干個子時間序列,隨后提取子時間序列的特征,若該特征超過設(shè)定的閾值,則認(rèn)為該序列為異常子時間序列;第2階段是采用局部離群因子(local outlier factor,LOF)算法檢測異常子時間序列中的異常值。異常值修正部分是對異常子序列中正常值點采用最小二乘法建立線性回歸模型,預(yù)測異常值點數(shù)值,最終異常值點數(shù)值取多個異常子時間序列預(yù)測值的均值。

1.1 異常子時間序列檢測

一條時間長度為n的時間序列表示為

X(t)=(x(t1),x(t2),…,x(tn))

式中x(ti)(1≤i≤n)為ti時刻記錄的數(shù)據(jù),采集時間ti是嚴(yán)格遞增的。

采用長度為w(w?n)的滑動窗口對該時間序列進(jìn)行等長度分割,步長為r,連續(xù)滑動(n-w)/r次,則會形成((n-w)/r)+1個長度為w的子時間序列。

1.1.1 直接方法提取子時間序列特征

假設(shè)時間序列X(t)的某一個子時間序列表示為Xj(1≤j

(1)

(2)

(3)

這樣就可計算第j個子時間序列斜率的置信區(qū)間半徑為

(4)

(5)

(6)

式中:Z為滿足N(0,1)的正態(tài)分布隨機(jī)變量;α為置信水平,本文取α=0.05。

1.1.2 快速方法計算子時間序列斜率均值和均方差

子時間序列特征值提取的關(guān)鍵在于相鄰兩點斜率的求值k(i),以及后續(xù)斜率的均值和均方差的計算。通常情況下時間序列會有較大的長度n,從而就會有(n-w+1)·(w-1)次斜率的計算量,以及(n-w+1)個滑動窗口的斜率均值和均方差的計算量,因此研究如何提高數(shù)據(jù)的運(yùn)算效率是非常有必要的。

第j個子時間序列任意兩點之間斜率的和與第(j+1)個子時間序列任意兩點之間斜率的和之間的關(guān)系為

(7)

故第(j+1)個子時間序列只需計算k(j+w)即可。

(8)

(9)

式(9)中的分子,由式(8)可得:

(10)

將式(10)帶入式(9)可得:

(11)

式(7)、式(8)和式(11)給出了前后2個子時間序列斜率、斜率均值和斜率均方差之間的關(guān)系,通過前一個子時間序列的斜率、斜率均值和斜率均方差計算當(dāng)前子時間序列的斜率、斜率均值和斜率均方差,會比直接計算更迅速。由表1可知,所提出的快速計算方法與滑動窗口長度w無關(guān),且明顯優(yōu)于直接計算方法。

表1 直接計算和快速計算的時間復(fù)雜度

得到子時間序列的特征,即斜率的置信區(qū)間半徑,若該半徑大于閾值γ,則認(rèn)為該子時間序列為異常子時間序列,其中包含異常值,具體的異常值信息需進(jìn)一步分析。

根據(jù)以上描述,異常子時間序列檢測的步驟可總結(jié)如下:

輸入:長度為n的時間序列,滑動窗口長度w,異常子時間序列判斷閾值γ。

輸出:異常子時間序列和正常子時間序列。

(1)使用長度為w的滑動窗口對時間序列X(t)進(jìn)行等長度劃分,步長為1,得到若干個子時間序列Xj(1≤j≤n-w+1);

(4)將子時間序列斜率的置信區(qū)間半徑與閾值γ進(jìn)行比較,初步確定含有異常值的異常子時間序列。

1.2 LOF算法檢測異常值

針對每個檢測出的異常子時間序列,采用LOF算法[19-20]計算每個數(shù)據(jù)點的離群因子,以此來識別異常值。

設(shè)異常子時間序列為Xl(1≤l≤n-w+1),LOF算法首先計算每個數(shù)據(jù)點的可達(dá)距離,通過可達(dá)距離計算數(shù)據(jù)點的局部可達(dá)密度,進(jìn)一步計算出數(shù)據(jù)點的局部離群因子。數(shù)據(jù)點的局部離群因子越大于1,說明該數(shù)據(jù)點的密度比其鄰域內(nèi)點的平均密度小,該數(shù)據(jù)點越可能是一個異常值點,反之,則說明該數(shù)據(jù)點是正常點。在本方法中,計算出異常子數(shù)據(jù)序列中每個數(shù)據(jù)點的局部離群因子,采用離群因子中最大值的1/2作為閾值,若一個數(shù)據(jù)點的離群因子大于閾值,則為異常值點,反之,則為正常點。

點p,o∈Xl,采用式(12)計算出點p到點o的第k可達(dá)距離,即:

RDk(p,o)=max{k-distance(o),dist(p,o)}

(12)

式中:k-distance(o)為點o的K-鄰近距離,是指距待測樣本點o的最近幾個樣本中,第k個最近的樣本點跟待測檢測點o之間的距離;dist(p,o)為點p到點o的歐式距離。

(13)

式中N為數(shù)據(jù)點的維度。

由式(14)計算點p的局部可達(dá)密度為

(14)

式中|Ndistance(p)|為到點p的距離小于k-distance(p)的集合。

由式(15)計算點p的局部離群因子為

(15)

根據(jù)以上描述,LOF算法檢測異常值步驟可總結(jié)如下:

輸入:異常子時間序列Xl。

輸出:異常子時間序列中的異常值。

(1)使用式(15)計算異常子時間序列中每個數(shù)據(jù)點的局部離群因子;

(2)選取局部離群因子中最大值的1/2作為閾值ρ;

(3)將每個數(shù)據(jù)點的離群因子與閾值ρ進(jìn)行比較,若大于閾值,則該數(shù)據(jù)點為異常值,否則為正常值。

1.3 基于多子時間序列的異常值填補(bǔ)

最小二乘法[21](least square method,LSM)是通過最小化模型輸出與真實值之間的誤差平方和計算樣本點的最佳匹配模型,在線性回歸中,就是試圖尋找到一條直線,使得所有樣本點到直線上的歐式距離之和最小。

在1.2節(jié)中LOF算法從異常子時間序列中檢測出具體的異常值點,對異常子時間序列中正常值點采用最小二乘法,建立線性回歸模型,預(yù)測異常值點位置的數(shù)值。本方法中滑動窗口步長r=1,所以異常值點存在于多個子時間序列中,記錄異常值點在盾構(gòu)機(jī)實測數(shù)據(jù)中的索引。若異常值填補(bǔ)數(shù)據(jù)集中無此索引,則創(chuàng)建該索引和記錄該子時間序列中預(yù)測異常值點的數(shù)值;否則,直接記錄該索引位置點的預(yù)測值。所有子時間序列檢測完后,每個異常值點的最終預(yù)測值取其在多個異常子時間序列中的預(yù)測值均值。

根據(jù)以上描述,基于多子時間序列的異常值填補(bǔ)步驟可總結(jié)如下:

輸入:異常子時間序列Xl和異常子時間序列中的異常值。

輸出:異常值填補(bǔ)數(shù)據(jù)集。

(1)對輸入的異常子時間序列中的正常值采用最小二乘法預(yù)測異常值點的數(shù)值;

(2)記錄異常值點在盾構(gòu)機(jī)實測數(shù)據(jù)中的索引,若異常值填補(bǔ)數(shù)據(jù)集中無該異常值點索引,則創(chuàng)建該索引和記錄預(yù)測后的數(shù)值,否則,直接記錄該索引位置點的預(yù)測值;

(3)若該子時間序列是最后子時間序列,計算異常值填補(bǔ)數(shù)據(jù)集中每個異常值點的多個子時間序列預(yù)測數(shù)值的均值,作為該異常值點最終的填補(bǔ)數(shù)據(jù);否則,返回1.1節(jié)異常子時間序列檢測中步驟(4)。

完整的盾構(gòu)機(jī)實測數(shù)據(jù)點異常處理方法流程圖如圖1所示。

圖1 盾構(gòu)機(jī)實測數(shù)據(jù)點異常處理方法流程圖

1.4 評價指標(biāo)

為了評價所提出方法的性能,點異常值檢測和異常值填補(bǔ)采用不同的指標(biāo)衡量。點異常值檢測采用召回率R和準(zhǔn)確率P,其數(shù)值越大,說明方法檢測性能越好,計算公式下:

(16)

(17)

異常值填補(bǔ)采用決定系數(shù)R2(coefficient of determination),其數(shù)值越接近1,表明填充性能越好,計算公式如下:

(18)

2 試驗結(jié)果及其分析

2.1 盾構(gòu)機(jī)實測數(shù)據(jù)基本情況及試驗設(shè)置

為了驗證所提方法的有效性和實用性,通過所收集的盾構(gòu)機(jī)實測數(shù)據(jù)集進(jìn)行檢驗。該數(shù)據(jù)來源于某城市地鐵修建中的隧道施工標(biāo)段,此隧道前段和后端選用盾構(gòu)掘進(jìn)施工,中段的硬巖階段選用傳統(tǒng)的礦山法施工。在盾構(gòu)掘進(jìn)施工過程中,隧道掘進(jìn)裝備為土壓平衡式盾構(gòu)機(jī),其主要系統(tǒng)參數(shù)如表2所示。

表2 盾構(gòu)機(jī)主要系統(tǒng)參數(shù)

搭載在該盾構(gòu)機(jī)的數(shù)據(jù)采集系統(tǒng)采樣頻率為1次/s。在掘進(jìn)過程中,傳感器采集了刀盤系統(tǒng)、土倉系統(tǒng)、推進(jìn)系統(tǒng)等子系統(tǒng)運(yùn)行參數(shù)的實測數(shù)據(jù)。

從盾構(gòu)機(jī)實測數(shù)據(jù)集中選取88環(huán)段分段運(yùn)行數(shù)據(jù)進(jìn)行試驗驗證。88環(huán)段包含1 605個數(shù)據(jù)樣本,選取盾構(gòu)機(jī)推力作為分析對象。在源數(shù)據(jù)中隨機(jī)插入20個異常值點,異常值位置依次是6、37、120、221、361、405、506、607、685、749、802、900、1 006、1 094、1 212、1 249、1 324、1 427、1 518、1 595。

本文使用的計算機(jī)配置為Windows10操作系統(tǒng),處理器為Intel(R) Core(TM) i7-6700 CPU @ 3.40 GHz,編程語言為Python3.6.13。

2.2 盾構(gòu)機(jī)實測數(shù)據(jù)驗證

在試驗前對源數(shù)據(jù)采用最小-最大規(guī)范化處理,公式如下:

(19)

式中:n為盾構(gòu)機(jī)實測數(shù)據(jù)的數(shù)量;min(X)為盾構(gòu)機(jī)實測數(shù)據(jù)的最小值;max(X)為盾構(gòu)機(jī)實測數(shù)據(jù)的最大值。

滑動窗口長度為7,步長為1,計算各子時間序列斜率置信區(qū)間半徑,最優(yōu)閾值設(shè)為γ=0.04,如圖2所示。將置信區(qū)間半徑大于閾值的子時間序列設(shè)為異常子時間序列,如圖3所示。每個被檢測出的異常子時間序列,采用LOF算法識別點異常值,如圖4所示。

圖2 子時間序列斜率置信區(qū)間半徑

圖3 檢測出的異常子時間序列

圖4 檢測出的異常值點

用準(zhǔn)確率、召回率指標(biāo)和執(zhí)行時間衡量不同算法的檢測結(jié)果,將AD-SV-LOF、AD-SV-GG、AD-Variance-LOF、原LOF算法與本文方法進(jìn)行比較。AD-SV-LOF方法與本文方法都是使用子時間序列斜率置信區(qū)間半徑作為該子時間序列的特征,但在計算過程中,AD-SV-LOF方法采用直接計算的方式,本文方法采用快速計算的方式。AD-SV-GG方法在異常子識別序列識別階段與本方法相同,但在異常值識別中采用了Gath-Geva聚類算法,AD-Variance-LOF方法在對異常子時間序列識別中選擇傳統(tǒng)的置信區(qū)間閾值判斷。采用10次重復(fù)實驗,各項指標(biāo)取其均值,結(jié)果如表3所示。

表3 異常檢測結(jié)果比較

在識別異常子時間序列時,本文方法采用子時間序列斜率的置信區(qū)間半徑作為子時間序列的特征,能夠較好地反映該子時間序列的變化特征,而AD-Variance-LOF方法采用子時間序列的均值和方差描述該子時間序列的結(jié)構(gòu)特征變化情況,不能很好反映子時間序列特征,因此,在召回率R、準(zhǔn)確率P評價指標(biāo)數(shù)值上低于本文方法。AD-SV-GG方法在異常值檢測階段采用Gath-Geva聚類算法,Gath-Geva算法的結(jié)果會受隨機(jī)生成的初始矩陣影響[17],只能給出是不是異常值點的判斷,相比較而言,基于密度的LOF算法通過給每個數(shù)據(jù)點都分配一個依賴于鄰域密度的離群因子,進(jìn)而判斷該數(shù)據(jù)點是否為離群點,因此本文方法的評價指標(biāo)優(yōu)于AD-SV-GG方法。相對于傳統(tǒng)的LOF算法,本文方法在檢測過程中考慮到了時間序列數(shù)據(jù)的結(jié)構(gòu)變化特征,先識別異常子時間序列,再使用LOF算法,因此本文方法的準(zhǔn)確率和召回率均優(yōu)于原LOF算法。

AD-Variance-LOF方法對當(dāng)前窗口子時間序列進(jìn)行初步斷定時,采用的子時間序列是均值和方差,該方法簡單且運(yùn)行成本低,所以執(zhí)行時間小于本文方法。AD-SV-LOF在提取子時間序列特征時均采用了直接方法,而本文方法采用了快速方法計算子時間序列的均值和方差,快速計算方法時間復(fù)雜度與滑動窗口長度無關(guān),隨著數(shù)據(jù)量的增大,優(yōu)勢越大,所以執(zhí)行時間比AD-SV-LOF方法少了約0.02 s,但在實際施工過程中,采集的盾構(gòu)機(jī)運(yùn)行數(shù)據(jù)量巨大,本文方法的優(yōu)勢會進(jìn)一步得到體現(xiàn)。

綜合以上結(jié)果分析,本文方法能夠更好地反映窗口內(nèi)子時間序列的形態(tài)信息以及內(nèi)在變化特征,因此本文方法在召回率R、準(zhǔn)確率P評價指標(biāo)數(shù)值上均優(yōu)于其余 4 種方法。盡管計算成本相比較于AD-Variance-LOF方法有所提升,考慮到本文方法具有較好的異常識別率,且對子時間序列特征采用了快速計算方式,提高了計算效率,故所提出的方法依舊具有較好的實用性。

最后采用基于多子時間序列對異常值進(jìn)行填補(bǔ),如圖5所示,填補(bǔ)結(jié)果的決定系數(shù)R2大于0.9,達(dá)到了預(yù)期。

圖5 多子時間序列異常值填補(bǔ)

為進(jìn)一步驗證本文所提方法在盾構(gòu)機(jī)實測數(shù)據(jù)的適用性,選取盾構(gòu)機(jī)其他重要的2個參數(shù),即推進(jìn)速度和刀盤扭矩分別進(jìn)行試驗驗證。推力、推進(jìn)速度和刀盤轉(zhuǎn)矩是3種不同的物理量,且具有不同的運(yùn)行趨勢。同樣的,對這2種參數(shù)隨機(jī)插入異常值點,用本文方法進(jìn)行檢測和修正試驗結(jié)果如表4所示。由表4數(shù)據(jù)可知,本文方法在盾構(gòu)機(jī)其他參數(shù)的處理上均取得良好結(jié)果,滿足了工程數(shù)據(jù)要求。

表4 盾構(gòu)機(jī)其他參數(shù)處理結(jié)果

3 結(jié)束語

本文針對盾構(gòu)機(jī)實測數(shù)據(jù)中產(chǎn)生的點異常值進(jìn)行處理,提出了一種基于滑動窗口的時間序列異常檢測方法,通過快速方法計算子時間序列的斜率置信區(qū)間半徑并識別異常子時間序列,再使用LOF算法對異常子時間序列進(jìn)行異常值檢測。試驗結(jié)果表明:與以子時間序列均值和方差信息提取子時間序列特征的檢測方法和原LOF算法相比,本文方法提高了檢測精度;與直接計算子時間斜率置信區(qū)間半徑的方法相比,本文方法提高了運(yùn)算效率。針對檢測出的異常值點,采用基于多子時間序列對異常值進(jìn)行填補(bǔ),試驗結(jié)果表明,填補(bǔ)效果達(dá)到了預(yù)期。本文的點異常值檢測和修正方法能夠滿足工程數(shù)據(jù)的處理要求,但對點異常值的產(chǎn)生原因尚未研究,因此下一步將結(jié)合盾構(gòu)機(jī)實測數(shù)據(jù)的空間相關(guān)性對點異常值產(chǎn)生原因進(jìn)行研究。