傳感器故障后多變量經驗小波變換多點預測*

李春祥， 張佳麗

(上海大學土木工程系 上海，200444)

引 言

目前超高層建筑，特別是600 m以上超高層建筑需要安裝結構健康監測(structural health monitoring，簡稱SHM)系統，通過對結構響應等結構系統特性分析來監測結構損傷或退化。SHM系統使用的傳感器屬于精密測量儀器，使用環境惡劣、操作不當和安裝不穩定等因素都會導致故障發生。據資料顯示，最嚴重的風災往往由颶風和雷暴產生。在風災發生時，一旦傳感器發生故障，對數據記錄造成缺失，后果難以挽回。另外，在結構振動的主動控制中，整個建筑物表面所受的力需要同時得知，若多個傳感器同時發生故障，單點預測不能滿足計算需要，因此多點同步預測模型的建立十分必要。

要同步恢復缺失信號，必須要考慮同步分解多變量信號。隨著多變量信號分析的需求增加，Rehman等[1]提出MEMD算法來代替傳統經驗模態分解(empirical mode decomposition, 簡稱EMD)這種單通道信號分解方法。該算法能同步處理安放在不同位置的傳感器采集來的多變量信號，保證了固有模態函數(intrinsic mode function，簡稱IMF)在數量和尺度上的統一。近年來，MEMD常被用于機械狀態的監測和故障的診斷，極大地提高了多變量信號分解的準確性，同時降低了運算的復雜程度。Yong等[2]驗證了MEMD結合非局部均值 (non-local means，簡稱NL-means)降噪算法以及故障相關因素分析在滾動軸承故障診斷中的有效性。Huang等[3]提出一種部分噪聲輔助的多元經驗模態分解(partial noise assisted multivariate EMD，簡稱PNA-MEMD)算法，利用高頻窄帶的噪聲來取代傳統的白噪聲，獲得穩定效果的同時簡化了運算。熊炘等[4]利用MEMD分解多通道振動信號來融合識別齒輪箱齒面點蝕故障信號的多通道數據。王恒等[5]采用自適應噪聲輔助的多元經驗模態分解(noise assisted multivariate EMD，簡稱NA-MEMD)來提取故障特征。段若晨等[6]提出一種窄帶噪聲輔助多元經驗模態分解(narrowband noise assisted multivariate empirical mode decomposition，簡稱NNA-MEMD)算法，用來檢測換流變壓器用有載分接開關的機械狀態。從以上發展看出，NA-MEMD是一種常用多變量分析的算法，能一定程度上解決MEMD的模態混疊問題。考慮到經驗小波變換(empirical wavelet transform, 簡稱EWT)在單點信號分解中有很好的效果，利用EWT基本框架構建MEWT算法，對故障后的多變量信號進行同步恢復，同時采用MIMO策略代替傳統的滾動策略進行多步預測，提高精度的同時能獲得更長時間的預測值，并與NA-MEMD模型進行對比，說明MEWT模型的優越性。

1 MEWT-KELM-CS-MIMO多點同步多步預測模型

1.1 EWT基本原理

近年來，在非線性、非平穩信號去噪方面，小波變換(wavelet transform, 簡稱WT)和EMD等都取得了一定效果，但WT在強噪聲情況下去噪效果會退化，EMD存在數學理論缺失、對噪聲和取樣敏感的問題[7]。Gilles等[8-9]基于小波變換和經驗模態分解存在的問題，提出了經驗小波變換。EWT能夠通過完全自適應小波基提取信號的固有模態，與經典小波變換一樣具有完備的理論基礎，可以顯著降低EMD類分解方法存在的模態混疊現象。

EWT的步驟總結如下：

(1)

2) 根據Meyer小波的構造方法構造一系列經驗小波。對于任意的?n>0，經驗尺度函數和經驗小波函數分別為

(2)

(3)

(4)

(5)

3) 重建信號。重建的序列和經驗模態為

(6)

(7)

(8)

1.2 MEWT技術框架

筆者提出多變量經驗小波變換的概念，通過將經驗小波變換分解后的模態進行相空間重構，模態數不等用零矩陣補齊后再重構來實現同步預測。

這里以三點信號為例來展示MEWT技術框架。首先，對每個信號通過EWT分解后的模態進行以下相空間重構，即將一維時間序列轉換成矩陣的形式。第i個模態的訓練集的時間序列{xi1,xi2,…,xin} 進行嵌入維度為d的相空間重構，具體的輸入和輸出樣本對為

設3個信號分解后最大模態數為m，若EWT分解后的模態數不足m的信號用零矩陣來補齊缺失的模態。

1.3 KELM基本概念

給定N個樣本{(xi,yi)|i=1,2,…,N}，其中xi=[xi1xi2… xin]T∈Rn，要尋找預測函數f:x→y，使得f(x)≈y。在隱層神經元數目為L的ELM中，這個預測函數的形式可表述為

(9)

在極限學習機(extremelearningmachine, 簡稱ELM)中，輸入權值和隱層節點偏置都是事先設定好的，激活函數也是選定的，所以ELM的訓練問題可歸結為求解輸出權值β的問題，建立如下的最優化問題

s.t.hT(xi)β=yi-εi(i=1,2,…,N)

(10)

其中：C為正則化參數，可以人為設定；εi為松弛變量，衡量實際值與預測值之間的誤差。

求解式(10)得到

(11)

其中：H=[h(x1) …h(xN)]T為隱含層輸出矩陣；Y=[y1… yN]T為輸出向量；IN為一個N×N維的單位矩陣。

類似于支持向量機(supportvectormachine, 簡稱SVM)等辦法，用核函數取代特征映射h(x)，定義核函數為k(x,y)=〈h(x),h(y)〉,這里〈·,·〉表示特征映射的內積。于是可以定義核矩陣Ω為

Ω=HHT：Ωi,j=h(xi)·h(xj)=k(xi,xj)

(12)

將核函數帶入式(9)和式(11)，得到KELM的預測函數為

(13)

研究顯示，KELM比起ELM有更好的泛化性能的同時需要更少的迭代參數[10]，KELM往往能和SVM達到相同的精度，但訓練預測花費的時間更少[11]，因此筆者選用KELM作為預測模型。徑向基核函數(radialbasiskernelfunction, 簡稱RBF)是適應性最好、使用最為普遍的一種核函數，其表達式為

K(x,xi)=exp(-‖x-xi‖2)/(2σ2)

(14)

其中：σ為RBF的寬度，也叫核參數。

1.4 CS算法概述

杜鵑搜索算法是文獻[12]提出的一種新興啟發算法,采用相關的Lévy飛行搜索機制。用杜鵑鳥的蛋來代表新的解，目的是使用新的和潛在的解來代替不那么好的解。該算法基于3個理想化規則：a. 每個杜鵑下一個蛋，堆放在一個隨機選擇的巢中；b. 最好的高品質的蛋巢將轉移到下一代；c. 巢數量固定，杜鵑的蛋被發現的概率為[0,1][13]。研究表明，杜鵑搜索比其他群體優化算法更有效。

筆者采用CS對KELM的正則化參數C和核函數參數σ進行優化，C的取值范圍設定為[2-8,28]，σ2的取值范圍設定為[0.01,100]，具體步驟如下：

1) 初始化參數C和σ，生成初始種群。

2) 將平均絕對誤差(mean absolute error，簡稱MAE)作為適應度函數計算每個鳥巢的適應度，求出種群最優位置。

3) 用xi(t+1)=xi(t)+α⊕levy(λ)更新鳥巢位置，其中：xi(t)代表第i個鳥巢在t代的鳥巢位置；α為步長；⊕代表點對點乘法；levy(λ)表示萊維(Lévy)隨機搜索路徑，且levy(λ)～u=t-λ(1<λ<3)。將現有鳥巢位置與上一代鳥巢位置進行對比，擇優作為當前最優位置。

4) 用隨機數r∈[0,1]與鳥巢主人發現外來鳥概率Pa對比，若r>Pa，則隨機地改變鳥巢位置，得到一組新的鳥巢位置。

5) 比較各鳥巢適應度值，更新當前鳥巢最優位置。

6) 滿足容許值停止迭代，否則重新執行步驟3。

7) 得到最優參數。

有研究表明，由于CS算法搜索過程采用Lévy飛行，短距離的探索與長距離探索相間，因此CS在迭代后期有更強的優化能力[14]。同時，有研究指出，CS算法和粒子群優化(particle swarm optimization, 簡稱PSO)算法都能收斂到全局最優，但是仍有機會陷入局部最優[15]，因此對每種情況運行10次求平均值得到最終結果[16]。

1.5 MIMO策略用于多步預測

在時間序列預測中，往往希望了解未來一段時間的數值或者趨勢，為了獲得更長時間的預測值而增大序列的間隔時間往往導致信號信息丟失和預測值不能用于實時調度的問題，多步預測就成為一種重要方法。因為采用多輸入單輸出函數，先前的多步預測策略，例如滾動策略、直接策略和直接滾動策略被認為是單輸出的策略，影響多步預測精度的主要因素有誤差的積累、準確性的降低和不確定性的增加[17]。MIMO策略能避免上述單輸出策略引起的未來預測值之間的聯系被丟失的情況，提高準確性，降低不確定性，同時MIMO策略一次輸出所有步長預測值，消除了滾動法中的誤差積累現象，因此MIMO 策略相比于其他多步預測策略有更高的預測精度。近年來，多篇文獻也反應MIMO策略相比于其他多步預測策略有更高的預測精度[17-19]。另外，與直接策略以及直接滾動策略需要構建多個模型相比，MIMO策略僅需一個模型，建模簡便。

MIMO策略僅需通過時間序列[x1,…,xN]來訓練一個多輸出模型F

[xt+H,…,xt+1]=F(xt,…,xt-d+1)+w

t∈{d,…,N-H}

(15)

以矩陣的形式表現更為直觀，MIMO策略訓練模型時輸入矩陣為

輸出矩陣為

(16)

筆者按照訓練集占比75%左右的標準，選取前800個點作為訓練集，即N取800，預測后200個點。在超前3步和超前6步的預測中，H分別取3和6，每次預測得到后一步長所有的預測值，循環得到后200個點的預測值。

1.6 信號恢復流程

筆者創新性地提出用MEWT-KELM-CS-MIMO模型來對多點缺失數據進行同步恢復，多步預測得到未來一段時間的預測值。整個流程框架如圖1所示。具體步驟如下：

1) 按MEWT的技術框架，將多點中斷信號前一定長度的時間序列進行同步相空間重構后組合成新的矩陣；

2) 將重構后的矩陣分為訓練集和測試集；

3) 歸一化處理能提高預測的精度和收斂的速度，將新組合成的矩陣歸一化到[-1,1]；

4) 徑向基核函數極限學習機用來訓練和預測，其中正則化參數和核參數用CS算法優化。

5) 將各模態的預測值相加得到預測結果，計算評價指標評價模型；

6) MIMO策略用來實現模型的1步預測、3步預測和6步預測，并與滾動策略對比。所提出的模型還與EWT-KELM-CS單點模型、NA-MEMD-KELM-

圖1 信號恢復流程圖Fig.1 Flowchart of signal recovery

CS多點同步模型通過評價指標進行對比，最終證明MEWT-KELM-CS-MIMO多點同步模型的優越性；

7) 多步預測結果可用于傳感器信號恢復，也可在結構振動控制的主動控制中提前計算結構下一步荷載。

2 基于實測風壓數據的模型驗證

對某辦公樓樓頂砌筑的矩形結構在2012年11月23日測得的實測風壓數據[20]進行研究。該測試結構位于辦公樓樓頂，視野開闊，當天風向為東北風，風力為3至4級，實測方案和實測平面布置圖如圖2所示。在結構AB墻表面，沿豎向每隔21 cm布置1#～5#風壓傳感器，DA墻表面，沿豎向每隔21 cm布置6#～10#風壓傳感器，1#和6#風壓傳感器距離結構頂面18 cm，同時，所有傳感器距角A水平距離均為23 cm。可知，6#～10#風壓傳感器位于迎風面，1#～5#風壓傳感器位于背風面，筆者取其中2#～4#和7#～9#傳感器的實測數據來分析。

現場采樣頻率為20 Hz，取1 000個數據點用于模型的訓練和測試，為了獲得較長的時間序列，采樣點之間間隔取0.8 s。

圖2 實測方案和平面布置圖Fig.2 Field measurement scheme and plane layout

2.1 迎風面模型驗證

假設迎風面3點數據同時缺失，用MEWT-KELM-CS-MIMO模型對這3點數據同步恢復。7#，8#和9#風壓數據用EWT分解后產生的模態數分別為11，12和11，嵌入維度d取10，MEWT確定模態數為12，不足12個模態的用零矩陣補齊后重構為新矩陣，之后將其分為訓練集和測試集兩部分。訓練集取790個30維向量，測試集取200個30維向量。3種常見模型和所提出的MEWT-KELM-CS-MIMO分別預測后對比，采用以下4個指標評價4種模型的預測效果。

平均絕對誤差

(17)

均方根誤差

(18)

相關系數

(19)

平均絕對相對誤差

(20)

由于篇幅限制，8#超前1步、3步和6步的預測結果及誤差分布圖如圖3所示。三點的各模型預測性能評價指標如表1所示。

圖3 8#預測結果和誤差分布圖Fig.3 Forecasting results and Error distribution diagram of 8#

模型1步預測3步預測6步預測MAERMSERMAPEMAERMSERMAPEMAERMSERMAPE7#A0.1320.1710.9990.1110.4380.5710.9920.3541.3301.5650.9401.278B0.8891.2560.9620.6761.7292.3730.8551.0782.1472.8810.7781.410C0.1940.2420.9990.1390.4710.5800.9920.4550.9771.2080.9651.031D0.1750.2210.9990.1380.4020.5190.9940.2490.7090.9250.9800.4688#A0.1560.2080.9990.0460.4690.6040.9910.1320.9101.1870.9650.285B0.8581.2110.9620.5041.6582.3100.8481.0592.0902.7620.7741.049C0.2080.2760.9980.0780.5420.7110.9870.1920.9481.2060.9610.311D0.2020.2680.9980.0730.3700.4880.9940.1200.7430.9780.9750.2979#A0.1220.1480.9990.0330.3920.4990.9920.1420.7130.9310.9720.215B0.7561.0650.9620.2141.4011.9140.8720.3941.6742.2630.8160.473C0.1860.2310.9980.0540.4580.5740.9890.1460.9291.1450.9560.315D0.1800.2230.9980.0510.3590.4690.9930.1150.6260.8490.9760.185

A為EWT-KELM-CS單點模型；B為NA-MEMD-KELM-CS；C為MEWT-KELM-CS；D為MEWT-KELM-CS-MIMO

2.2 背風面模型驗證

假設背風面3點數據同時缺失，用MEWT-KELM-CS-MIMO模型對這3點數據同步恢復。2#，3#和4#風壓數據用EWT分解后產生的模態數分別為9，8和8，嵌入維度d取10，MEWT確定模態數為9，不足9個模態的用零矩陣補齊后重構為新矩陣，之后將其分為訓練集和測試集兩部分。訓練集取790個30維向量，測試集取200個30維向量。3種常見模型和所提出的MEWT-KELM-CS-MIMO分別預測后對比，采用上述評價指標對各模型進行評價。3#超前1 步、3步和6步的預測結果及誤差分布圖如圖4所示。3點的各模型預測性能評價指標如表2所示。

由表2可知，在超前1步預測時，提出的模型接近EWT-KELM-CS單點模型，在超前3步預測和超前6步預測中，對于3#和4#，所提出的模型要優于其他3種模型，對于2#，模型D效果和模型A在超前3步預測時相差不多，但是優于模型B，C。在超前3步預測時，2#模型D相對于模型B，C的ρMAE達到78.4%，17.5%；3#模型D相對于模型A，B，C的ρMAE達到9.7%，76.0%，28.8%；4#模型D相對于模型A，B，C的ρMAE達到9.7%，71.3%，25.4%。在超前6步預測時，2#所提出的模型D相對于模型A，B，C的ρMAE達到30.8%，57.7%，33.2%；3#模型D相對于模型A，B，C的ρMAE達到33.7%，54.8%，35.8%；4#模型D相對于模型A，B，C的ρMAE達到20.4%，44.6%，33.3%。

2.3 雙面模型驗證

隨機選取迎風面和背風面的數據點用MEWT-KELM-CS-MIMO對數據進行同步預測，結果和單獨預測迎風面或單獨預測背風面結果相差不多。同時，從以上結果看出，3點同步預測結果中中間點和兩端點效果相近，反映出筆者采用的多點同步模型不受空間的限制，僅用各點以往的數據樣本來訓練模型，同步預測之后的信號能達到比空間點預測更高的精度。

圖4 3#預測結果和誤差分布圖Fig.4 Forecasting results and Error distribution diagram of 3#

模型1步預測3步預測6步預測MAERMSERMAPEMAERMSERMAPEMAERMSERMAPE2#A0.0890.1150.9990.3150.3130.4480.9930.5181.0991.3850.9302.493B0.7531.0640.9571.4741.5302.1080.8173.2951.7972.3930.7577.862C0.1210.1570.9990.3670.4010.5190.9911.6261.1391.4730.9235.576D0.1270.1620.9990.3320.3310.4540.9921.3870.7611.0970.9548.2363#A0.0960.1290.9990.7110.3720.4810.9912.7941.2071.5540.9068.253B0.7130.9830.9627.4681.4061.9660.83110.7821.7702.3390.7517.834C0.1490.1890.9991.1770.4730.5970.9863.8421.2461.5800.9077.846D0.1430.1810.9990.7600.3370.4690.9911.4840.8001.1640.9451.6634#A0.0970.1340.9990.0770.4230.5370.9890.2981.1821.4980.9131.069B0.7111.0320.9590.5181.3331.8520.8601.2501.6992.3100.7711.130C0.1410.1780.9990.0920.5120.6470.9840.3261.4111.7870.8781.396D0.1370.1750.9990.0700.3820.5170.9900.2720.9411.3440.9300.899

A為EWT-KELM-CS單點模型；B為NA-MEMD-KELM-CS；C為MEWT-KELM-CS；D為MEWT-KELM-CS-MIMO

3 基于實測強非平穩風速數據的模型驗證

Derecho是一種典型的風暴，通常可以存活超過8 h以上，其特點在于具有大面積超過65節(1節=1.852 km/h)的持續強風，移速較高，破壞力很強。2002年5月20日至7月15日德克薩斯理工大學大氣科學系在瑞茜技術中心由北向南布置了7個便攜式塔，測得了Derecho水平風速數據。塔間間隔為263 m，其中塔1和塔7最高觀測點位置高為3 m，塔2，3，5，6最高觀測點位置高為10 m，塔4最高觀測點位置高為15 m。選取塔4上高度為3.96, 6.10以及10.06 m處的水平風速數據來驗證多點同步模型，將3.96 m處的點記為1#，6.10 m處的點記為2#，10.06 m處的點記為3#。各塔布置圖及觀測點位置如圖5所示。原始采樣頻率為2 Hz，筆者采用1 Hz ，各點實測數據如圖6所示。以3#為例，運用增廣的迪基-福勒檢驗法(augmented dickey-fuller test，簡稱ADF)檢驗時間序列平穩性。若存在單位根，則為非平穩時間序列；否則為平穩時間序列。風速樣本檢驗值為-2.820，大于1%，5%顯著性水平下的臨界值-3.437和-2.864，存在單位根的原假設成立，所以樣本為非平穩時間序列。2#超前1步、3步和6步的預測結果及誤差分布圖如圖7所示。2#的各模型預測性能評價指標如表3所示。

圖5 各塔布置圖及觀測點位置Fig.5 The layout of towers and the location of observation points

圖6 各點實測數據圖Fig.6 Original wind speed time series

由表3可得，在超前1步預測和超前3步預測中，提出的模型效果接近EWT-KELM-CS單點模型，優于其他2種模型。但在超前6步預測時，模型D相對于模型A，B，C的ρMAE達到15.1%，62.3%，28.0%。

圖7 2#預測結果和誤差分布圖Fig.7 Forecasting results and Error distribution diagram of 2#

模型1步預測3步預測6步預測MAERMSERMAPEMAERMSERMAPEMAERMSERMAPE2#A0.0190.0260.9990.0010.1170.1530.9970.0080.5220.6640.9510.037B0.4550.5750.9690.0320.7670.9780.9060.0551.1761.4950.8020.086C0.1020.1270.9980.0070.3270.4340.9790.0230.6150.7860.9390.043D0.0970.1230.9990.0070.2280.3370.9880.0160.4430.5850.9690.031

A為EWT-KELM-CS單點模型；B為NA-MEMD-KELM-CS；C為MEWT-KELM-CS；D為MEWT-KELM-CS-MIMO

4 結 論

1) MEWT-KELM-CS-MIMO的多點同步多步預測模型預測精度高，與EWT-KELM-CS單點模型接近，滿足精度要求，同時發展多點同步多步預測模型能大大提升工程應用效率。

2) MEWT-KELM-CS-MIMO模型的預測精度高于MEWT-KELM-CS模型，說明多步預測時MIMO策略優于滾動策略，與國際文獻結論一致。

3) 隨著步數的增加，MEWT-KELM-CS-MIMO模型精度可以超過EWT-KELM-CS單點模型，可以實現更長時間預測的同時保證精度。多點同步多步預測模型能提高計算的效率，比單點模型具有更大的工程應用價值。

4) 傳統使用NA-MEMD對多維信號進行分解的方法精度較低，因為加入噪聲的MEMD不能完全解決模態混疊問題，采用有良好數學基礎的EWT能很好地解決這個問題，但是無法保證多點分解后模態數一致。筆者創新性采用零矩陣來補齊模態數，提出MEWT概念，得到更高精度的同步預測結果，證明了MEWT結合傳統的KELM適用于多點同步預測。

5) 若傳感器徹底損壞，無法帶入新鮮樣本，只能恢復后面一個步長的數據，要想恢復很長一段時間的數據，可以將采樣間隔加大或將預測值不斷滾動迭代實時分解后預測下一步長。同時多點同步多步預測模型在結構振動控制的主動控制中應用前景廣闊，可以提前計算結構整個面的受力，讓結構提前知道下一步的反應，減小時滯，提高效率。