陳良捷，魏博文，戴國強(qiáng)，林太清

(1. 江西省水利科學(xué)院，江西 南昌 330029； 2. 南昌大學(xué) 工程建設(shè)學(xué)院，江西 南昌 330031)

大壩安全監(jiān)測是融合水工結(jié)構(gòu)學(xué)、儀器儀表學(xué)、計算統(tǒng)計學(xué)和網(wǎng)絡(luò)技術(shù)的“醫(yī)院”系統(tǒng)，廣泛應(yīng)用于大壩等水工建筑物的服役環(huán)境、變形效應(yīng)、滲流效應(yīng)、應(yīng)力應(yīng)變等數(shù)值量測[1-3]。原型監(jiān)測時序主要來源于監(jiān)測儀器對建筑物固定測點(diǎn)進(jìn)行一定頻次的自動化或人工觀測?；炷翂伪O(jiān)測項目主要分為5 類：環(huán)境監(jiān)測、變形監(jiān)測、滲流監(jiān)測、應(yīng)力應(yīng)變監(jiān)測、專項特殊監(jiān)測，其中變形監(jiān)測時序是壩體結(jié)構(gòu)在多元復(fù)雜環(huán)境共同驅(qū)動下產(chǎn)生變形效應(yīng)的直接表現(xiàn)，是混凝土壩服役健康診斷的重要基礎(chǔ)[4-5]。近年來，隨著計算機(jī)技術(shù)的發(fā)展及筑壩監(jiān)測硬件配置完備性的提高，大壩數(shù)據(jù)采集方式愈加多樣化、自動化，存儲機(jī)制與規(guī)模也逐步成熟龐大[6]。但由于設(shè)備更替修繕及人為管理疏漏等不確定因素，監(jiān)測時序中難免出現(xiàn)異常、空缺等現(xiàn)象，對海量數(shù)據(jù)集進(jìn)行科學(xué)高效的可靠性辨識及預(yù)處理是構(gòu)建大壩安全監(jiān)控體系的關(guān)鍵首步[7-10]。

國內(nèi)外學(xué)者針對大壩安全監(jiān)測數(shù)據(jù)集的可靠性辨識及異常突變識別技術(shù)方面展開了研究。鄭霞忠等[11]基于環(huán)境自變量與效應(yīng)因變量的關(guān)聯(lián)性約束機(jī)制提高了監(jiān)測異常數(shù)據(jù)清洗的準(zhǔn)確性，通過密度聚類算法識別時序中離群孤立點(diǎn)集，并依據(jù)關(guān)聯(lián)信號清洗規(guī)則有效辨識粗差部分，借助粒子群優(yōu)化的最小二乘支持向量機(jī)模型補(bǔ)齊數(shù)值空缺；西班牙學(xué)者Salazar 等[12]基于強(qiáng)化回歸樹模型計算結(jié)果與實際觀測之間的差異判斷疑似異常效應(yīng)序列，并輔以溫度場及靜力場耦合的有限元仿真驗證其可靠性，有效規(guī)避了因異常載荷而導(dǎo)致的虛假異?，F(xiàn)象；王麗蓉等[13]將大壩原始監(jiān)測時序繪制成數(shù)據(jù)過程線圖，基于卷積神經(jīng)網(wǎng)絡(luò)自主解析圖像特征的能力優(yōu)勢，識別序列點(diǎn)集的異常離群類型并搜尋具體位置；鄭東健等[14]為保護(hù)觀測資料樣本整體性并準(zhǔn)確反映異常干擾對各效應(yīng)分量的局部影響，構(gòu)建了指示變量數(shù)值模型以有效模擬各種系統(tǒng)干擾的影響量值。但是，目前現(xiàn)有的研究方案及成果未將效應(yīng)監(jiān)測時序的自身結(jié)構(gòu)特性融入其中，致使復(fù)雜混頻信號的內(nèi)部交互干擾成為異常突變識別工作的障礙。因此，如何開展科學(xué)可靠的監(jiān)測數(shù)據(jù)可靠性辨識及預(yù)處理工作亟需更深層次的研究分析。

針對混凝土重力壩變形監(jiān)測時序中異常突變值的識別難題，本文依據(jù)經(jīng)驗小波變換理論提取序列信號極大值點(diǎn)并自適應(yīng)地分割傅里葉頻譜，通過構(gòu)造正交小波濾波器組獲取信號的不同模態(tài)，進(jìn)而排除監(jiān)測時序中趨勢分量對異常突變診斷工作的干擾；以待測數(shù)據(jù)對象周圍領(lǐng)域的局部數(shù)據(jù)分布密度為檢測指標(biāo)，重點(diǎn)識別局部異常因子遠(yuǎn)大于1 的數(shù)據(jù)點(diǎn)，最后將計算結(jié)果通過統(tǒng)計箱形圖精準(zhǔn)找到對應(yīng)的數(shù)據(jù)序號以便剔除。

1 混凝土壩變形監(jiān)測時序的內(nèi)部信息特征

混凝土壩變形規(guī)律的復(fù)雜性導(dǎo)致其內(nèi)在的多元映射關(guān)系難以準(zhǔn)確描述，通?；陂L時間的監(jiān)測數(shù)據(jù)建立回歸模型[15-16]，根據(jù)各環(huán)境因子對因變量作用的顯著程度，從大到小依次逐個引入方程式，最大限度展示變形效應(yīng)的驅(qū)動來源。由于現(xiàn)場監(jiān)測技術(shù)的局限性及各類環(huán)境影響的強(qiáng)度差異，混凝土壩變形按照經(jīng)典分析理論可分為3 個部分：水壓分量 δH、溫度分量 δT和時效分量 δθ[17]，以混凝土重力壩為例，其擬合公式可表示為：

式中：H為水深；an為回歸模型水壓分量的擬合系數(shù)；b1i和b2i為回歸模型溫度分量的擬合系數(shù)；i為時間周期，i=1表示年周期，i=2表示半年周期；t為從建模資料序列到始測日的累計天數(shù)，t0為從建模資料序列的第1 天到始測日的累計天數(shù)；c1和c2為回歸模型時效分量的擬合系數(shù)；θ=t/100。

但是，現(xiàn)場觀測的項目種類有限且經(jīng)典統(tǒng)計模型產(chǎn)生的殘差序列 δR較大，無法解釋結(jié)構(gòu)變形效應(yīng)的全部因果關(guān)系，混凝土壩變形時序中各個組成分量大體趨勢如圖1 所示。從圖1 可以看出，水壓分量 δH、溫度分量 δT和時效分量δθ的趨勢曲線較為平穩(wěn)，而殘差序列 δR表現(xiàn)出強(qiáng)非線性、高頻率的不穩(wěn)定特征。因此，構(gòu)建混凝土壩變形監(jiān)測時序異常點(diǎn)診療體系的關(guān)鍵首步是解譯混頻信號的復(fù)雜多元內(nèi)部信息。

圖1 變形監(jiān)測數(shù)據(jù)各組成序列的近似趨勢Fig. 1 Approximate trend of each component sequence of the displacement monitoring data

2 基于EWT-LOF-BOXPLOT 的異常突變值診斷方案

2.1 經(jīng)驗小波變換

經(jīng)驗小波變換（Empirical Wavelet Transform, EWT）[18-19]是一種新型多分量信號處理方法，其主要思想是基于傅里葉支撐選擇1 組正交小波濾波器組，其中包括1 個低通濾波器和N?1個帶通濾波器，分別對應(yīng)趨勢分量和細(xì)節(jié)分量，隨后對調(diào)幅-調(diào)頻模態(tài)進(jìn)行希爾伯特變換，以獲取信號的瞬時頻率和瞬時幅值。該方法充分融合了小波變換（Wavelet Transform, WT）和經(jīng)驗?zāi)B(tài)分解（Empirical Mode Decomposition, EMD）的優(yōu)勢，適合處理非線性、非平穩(wěn)特性數(shù)據(jù)的分解難題[20-21]。

假設(shè)大壩原型變形監(jiān)測信號 δ(t)表達(dá)式為：

式中：δn(t)為調(diào)幅-調(diào)頻信號；An(t)為調(diào)幅部分函數(shù)；Fn(t)為調(diào)頻部分函數(shù)。為完全自適應(yīng)分割監(jiān)測信號 δ(t)的頻譜，將傅里葉支撐區(qū)間[0, π]劃分為N個連續(xù)部分，即：

式中： ?i為不同子體部分，i=1,2,···,N； θi為不同子體部分的區(qū)間邊界，θ0=0且θN=π。設(shè)定以 θi為中心，寬度為2 ρ的過渡區(qū)間，為每個子體部分 ?i構(gòu)造以Meyer 小波函數(shù)為基礎(chǔ)的經(jīng)驗小波帶通濾波器。經(jīng)驗小波函數(shù)和經(jīng)驗尺度函數(shù)表達(dá)式為：

式中：λ(x)為區(qū)間[0, 1]內(nèi)滿足K階導(dǎo)的任意函數(shù)； ρi為過渡階段的寬度參數(shù)。兩者表達(dá)式分別為：

信號 δ(t)與經(jīng)驗小波函數(shù)fi(θ)的內(nèi)積為細(xì)節(jié)系數(shù)C(i,t)，信號 δ(t)與經(jīng)驗尺度函數(shù)g1(θ)的內(nèi)積為趨勢系數(shù)C(0,t)，具體表達(dá)式為：

式中：IFFT?1()為快速傅里葉逆變換；f?i(θ)和g?1(θ)分別為fi(θ)和g1(θ)的傅里葉變換；fi(θ)和g1(θ)分別為fi(θ)和g1(θ)的共軛復(fù)數(shù)。結(jié)合式（2）和（8），原信號的重構(gòu)表達(dá)式及單分量經(jīng)驗?zāi)B(tài)分別為：

式中： ?為卷積運(yùn)算；C?(0,t)和C?(i,t)分別為C(0,t)和C(i,t)的傅里葉變換形式。原始真實信號通過EWT 方法可自適應(yīng)地分解成許多具有緊湊支持頻譜的經(jīng)驗?zāi)B(tài)函數(shù)。

2.2 局部異常因子

局部異常因子（Local Outlier Factor, LOF）檢測方法是Breunig 等[22]于2000 年提出的基于密度的無監(jiān)督離群點(diǎn)識別工具。其主要手段是通過比較待測數(shù)據(jù)點(diǎn)與其領(lǐng)域數(shù)據(jù)點(diǎn)的聚集密度，以此密度比值為依據(jù)判斷該待測數(shù)據(jù)點(diǎn)是否為離群異常點(diǎn)。相比其他基于距離、統(tǒng)計分布的檢測方法，LOF 僅考慮待測數(shù)據(jù)對象周圍領(lǐng)域的局部數(shù)據(jù)分布，可避免其他領(lǐng)域外數(shù)據(jù)集的干擾[23-25]。

假設(shè)存在由j個i維數(shù)據(jù)點(diǎn)組成的樣本集合S，從中任意取兩個樣本點(diǎn)M=(xm1,xm2,···,xmi)與N=(xn1,xn2,···,xni)，其中m和n皆取自[1,j]，則兩點(diǎn)之間的歐幾里得距離表達(dá)式為：

定義Dk(M)為數(shù)據(jù)點(diǎn)M的第k距離，若Dk(M)=D(M,N)，則數(shù)據(jù)點(diǎn)N是指在樣本集合S中除了數(shù)據(jù)點(diǎn)M之外距離數(shù)據(jù)點(diǎn)M最近的第k個點(diǎn)。假設(shè)Ek(M)為數(shù)據(jù)點(diǎn)M的第k距離鄰域，則其表達(dá)式為：

由式（12）可知，數(shù)據(jù)點(diǎn)M的第k距離鄰域Ek(M)是指在樣本集合S中除了數(shù)據(jù)點(diǎn)M之外的所有距離數(shù)據(jù)點(diǎn)M小于或等于數(shù)據(jù)點(diǎn)M的第k距離Dk(M)的數(shù)據(jù)點(diǎn)集合，相關(guān)示意圖如圖2 所示。

Mk圖2 數(shù)據(jù)點(diǎn) 的第距離鄰域示意Fig. 2 The distance neighborhood diagram of data point k M

定義Dk(M,L)為數(shù)據(jù)點(diǎn)L與數(shù)據(jù)點(diǎn)M之間的第k可達(dá)距離，其實際意義是指當(dāng)數(shù)據(jù)點(diǎn)L位于數(shù)據(jù)點(diǎn)M的第k距離鄰域內(nèi)時，Dk(M,L)只能等于數(shù)據(jù)點(diǎn)M的第k距離Dk(M)。當(dāng)數(shù)據(jù)點(diǎn)L位于數(shù)據(jù)點(diǎn)M的第k距離鄰域外時，Dk(M,L)為兩點(diǎn)之間的歐幾里得距離D(M,L)，數(shù)學(xué)表達(dá)式為：

設(shè)定ρk(M)為數(shù)據(jù)點(diǎn)M的局部可達(dá)密度，可理解為在數(shù)據(jù)點(diǎn)M的第k距離鄰域Ek(M)中除了數(shù)據(jù)點(diǎn)M外的所有點(diǎn)與數(shù)據(jù)點(diǎn)M之間的平均第k可達(dá)距離的倒數(shù)，表達(dá)式為：

由此可得在數(shù)據(jù)點(diǎn)M的第k距離鄰域Ek(M)中除了數(shù)據(jù)點(diǎn)M外所有點(diǎn)的局部可達(dá)密度，將其分別與數(shù)據(jù)點(diǎn)M自身局部可達(dá)密度的比值疊加求平均，即得數(shù)據(jù)點(diǎn)M的局部異常因子，其表達(dá)式為：

若Lk(M)值接近1，則數(shù)據(jù)點(diǎn)M的第k距離鄰域Ek(M)內(nèi)的數(shù)據(jù)點(diǎn)的密度相差不大，數(shù)據(jù)點(diǎn)M與其鄰域內(nèi)其他數(shù)據(jù)點(diǎn)屬于同簇；若Lk(M)值遠(yuǎn)小于1，則數(shù)據(jù)點(diǎn)M的密度比其鄰域內(nèi)其他數(shù)據(jù)點(diǎn)的密度大，數(shù)據(jù)點(diǎn)M可被判定為密集點(diǎn)；若Lk(M)值遠(yuǎn)大于1，則數(shù)據(jù)點(diǎn)M的密度比其鄰域內(nèi)其他數(shù)據(jù)點(diǎn)的密度小，數(shù)據(jù)點(diǎn)M可被判定為異常點(diǎn)。

2.3 離群點(diǎn)檢測模型的建成方案

混凝土重力壩長期深處于靜力場和溫度場中，其變形效應(yīng)也備受水壓、溫度等周期性影響，因此變形監(jiān)測數(shù)據(jù)存在一定的趨勢分量，對LOF 異常點(diǎn)檢測工作產(chǎn)生較大影響。為了確保檢測模型的精準(zhǔn)運(yùn)行，采用EWT 信號處理方法對原型監(jiān)測數(shù)據(jù)進(jìn)行分解，去除變形監(jiān)測時序的運(yùn)行趨勢。異常值屬于短時突變值，用此檢測方案不但不會遺漏任何離群點(diǎn)信息，還有利于消除多變工況對異常數(shù)據(jù)點(diǎn)檢測的影響，提高離群點(diǎn)檢測模型的可靠性。

為了進(jìn)一步精準(zhǔn)確定數(shù)據(jù)點(diǎn)SLOF的判斷閾值，且避免人為設(shè)定帶來的諸多不確定因素，本方案采取可反映數(shù)據(jù)點(diǎn)分散情況的箱形圖。箱形圖于1977 年由美國統(tǒng)計學(xué)家John Tukey 發(fā)明，能顯示出一組數(shù)據(jù)的最大值、最小值、中位數(shù)及上下四分位數(shù)，既不需要對待檢數(shù)據(jù)進(jìn)行服從特定分布的事先假定操作，也沒有任何限制性要求，可真實直觀地表現(xiàn)數(shù)據(jù)形狀的原生面貌。

根據(jù)處理多種數(shù)據(jù)試算經(jīng)驗，箱形圖已經(jīng)定義了異常值識別標(biāo)準(zhǔn)：

式中： δ′′為待檢數(shù)據(jù)的異常值；Qu和Qd分別為待檢數(shù)據(jù)的上下四分位數(shù)；sIQR為待檢數(shù)據(jù)的四分位距，即sIQR=Qu?Qd。箱形圖對數(shù)據(jù)集異常值的識別標(biāo)準(zhǔn)是基于具有一定耐抗性的四分位數(shù)和四分位距，避免了異常值對自身產(chǎn)生的擾動影響，增強(qiáng)了整體模型的穩(wěn)定性和準(zhǔn)確性。

綜上所述，混凝土重力壩變形監(jiān)測時序的離群點(diǎn)檢測優(yōu)化模型建成方案如圖3 所示。

圖3 優(yōu)化診斷方案流程Fig. 3 Flow chart of optimized diagnosis scheme

3 工程應(yīng)用

某混凝土重力壩位于浙江省境內(nèi)，最大壩高115 m，壩頂全長466.5 m，水庫總庫容2.2×1010m3。該混凝土壩設(shè)有變形、滲流、應(yīng)力應(yīng)變、裂隙裂縫等監(jiān)測項目，其中位移監(jiān)測主要有正垂線（PL）、倒垂線（IP）、壩頂引張線（EX）及多點(diǎn)位移計（M）等，其位移監(jiān)測系統(tǒng)布置如圖4 所示。

圖4 大壩位移監(jiān)測系統(tǒng)布置示意（單位：m）Fig. 4 Schematic layout of dam displacement monitoring system (unit: m)

本文選取該混凝土壩第4 號壩段115 m 高程處的水平位移自動化監(jiān)測時序 δ115為研究對象，其時間區(qū)間為1998 年4 月20 日至2018 年4 月30 日，并規(guī)定以向下游方向位移為正，向上游方向位移為負(fù)。由于該混凝土壩監(jiān)測系統(tǒng)更新改造過數(shù)次且管理方案有改變，致使其時間序列步長各不相等且無規(guī)律可尋，經(jīng)整理后得其序列數(shù)據(jù)容量為757 組。

由于大壩在運(yùn)行期間會受到庫水位、溫度等工況變化且離群異常點(diǎn)具有短時突變特性，故大部分隱藏于數(shù)據(jù)集的高頻區(qū)間段內(nèi)。因此在對該混凝土重力壩位移監(jiān)測時序進(jìn)行LOF 異常點(diǎn)檢測之前，需利用EWT 信號處理方法去除監(jiān)測時序的趨勢分量，將分解出的趨勢分量與原始數(shù)據(jù)進(jìn)行對比，結(jié)果見圖5。

圖5 趨勢分量與原始數(shù)據(jù)對比Fig. 5 Comparison between trend component and original data

由圖5 可知，大部分?jǐn)?shù)據(jù)點(diǎn)位于監(jiān)測時序的趨勢曲線上。通過整理得757 組水平位移序列的平均值為0.475 mm，最大值為5.270 mm，最小值為?4.800 mm，區(qū)間跨度較大，且存在個別疑似異常值。

在通過EWT 信號處理方法對該位移監(jiān)測時序進(jìn)行分解處理后，人工剔除蘊(yùn)藏于原始信號中的趨勢分量，對應(yīng)每個數(shù)據(jù)序號重構(gòu)剩余的高頻分量，再對重構(gòu)后的位移效應(yīng)細(xì)節(jié)信號進(jìn)行異常點(diǎn)檢測以求出各數(shù)據(jù)點(diǎn)的SLOF，結(jié)果如圖6 所示。

圖6 LOF 異常點(diǎn)檢測結(jié)果Fig. 6 The detection result diagram of abnormal points

由圖6 可知，基于大壩位移時序細(xì)節(jié)信號計算出的SLOF大多數(shù)位于1 附近，說明大多數(shù)據(jù)點(diǎn)第k距離鄰域內(nèi)的數(shù)據(jù)密度相差不大，即這些數(shù)據(jù)點(diǎn)與其鄰域內(nèi)其他數(shù)據(jù)點(diǎn)屬于同簇。通過整理計算，該細(xì)節(jié)信號序列SLOF平均值為2.057，最大值為8.563，最小值為0.783，中位數(shù)為1.604，下四分位數(shù)（25%）為1.059，上四分位數(shù)（75%）為2.685?，F(xiàn)將其LOF 數(shù)據(jù)集作箱形圖以尋找異常值，其結(jié)果如圖7 所示。

圖7 箱形圖檢測結(jié)果及異常值對應(yīng)的序號位置Fig. 7 Boxplot results and number positions for exceptions

在箱形圖7 中，上下兩條末端橫線并非準(zhǔn)確的1.5 倍箱子長度，而是不超過該長度的最遠(yuǎn)值，即表示該批數(shù)據(jù)正常值的分布區(qū)間。根據(jù)箱形圖7 顯示的數(shù)據(jù)序號可對原始數(shù)據(jù)進(jìn)行對應(yīng)剔除，得出異常突變值共32 組，優(yōu)化處理后得監(jiān)測數(shù)據(jù)725 組。

為了驗證該基于EWT-LOF-BOXPLOT 的異常突變值診斷方案對監(jiān)測數(shù)據(jù)集優(yōu)化的有效性和實用性，將優(yōu)化前后及單一LOF 方法檢測處理的數(shù)據(jù)分別代入經(jīng)典多元回歸模型中參與計算，得到的模型誤差擬合結(jié)果如圖8 所示；并引入統(tǒng)計學(xué)領(lǐng)域中衡量計算精準(zhǔn)度的評判指標(biāo)：平均絕對誤差EMA、均方誤差EMS、平均絕對百分誤差EMAP及決定系數(shù)R2，計算結(jié)果如表1 所示。

表1 計算精度的評判指標(biāo)Tab. 1 Evaluation index of calculation accuracy

圖8 優(yōu)化前后的數(shù)據(jù)擬合結(jié)果Fig. 8 Data fitting results before and after optimization

綜合對比圖8 與表1 可見，通過經(jīng)驗小波變換去除監(jiān)測時序的運(yùn)行趨勢，使樣本數(shù)據(jù)集的LOF 識別檢測更加精確，提高異常突變值的診斷效力；監(jiān)測樣本數(shù)據(jù)集經(jīng)過基于EWT-LOF-BOXPLOT 的異常突變值診斷方案優(yōu)化后，參與回歸擬合計算的平均絕對誤差、均方誤差、平均絕對百分誤差皆小于優(yōu)化前及單一LOF 檢測識別，且R2值更接近于1，這表明該方案優(yōu)化處理后的混凝土壩變形監(jiān)測時序更適用于回歸建模，有助于后續(xù)研究的監(jiān)控模型建立及反演反饋分析。

4 結(jié) 語

（1）混凝土壩變形時序因監(jiān)測設(shè)備更新改造或人為管理疏漏等不確定因素而產(chǎn)生異常突變點(diǎn)值，為提高原型觀測數(shù)據(jù)的真實性及計算可靠性，提出了一種針對混凝土壩變形監(jiān)測時序中異常突變值的優(yōu)化診斷方案。

（2）考慮到變形監(jiān)測時序是由多元環(huán)境組合驅(qū)動的綜合表現(xiàn)，其內(nèi)部信息特征復(fù)雜且存有多頻信號干擾，采用經(jīng)驗小波變換理論和箱形圖統(tǒng)計手段對局部異常因子計算模式進(jìn)行優(yōu)化，有效規(guī)避了單一檢測模型在進(jìn)行離群點(diǎn)識別時的判診局限。

（3）工程實例分析表明，相比于經(jīng)典局部異常因子計算機(jī)制，本文提出的異常點(diǎn)診療體系的檢測能力更優(yōu)，且建模方法簡便高效，加以修改可推廣應(yīng)用于大壩其他監(jiān)測效應(yīng)時序的優(yōu)化診斷。但是，由于本方案是基于效應(yīng)樣本數(shù)據(jù)集進(jìn)行自體內(nèi)部運(yùn)算，未引入周邊環(huán)境資料的考證，故對具有調(diào)節(jié)性能的水庫大壩存在一定局限性。