施工隧道數據處理中的可靠卡爾曼濾波①

張詩琪, 張晉豫

(北京交通大學 計算機與信息技術學院, 北京 100044)

隨著城市化的演進, 地鐵以快速便捷, 載客量大等優勢, 逐步發展為城市公共交通的首選. 作為城市發展的標志, 地鐵不僅為人民日常出行帶來了便利, 也為拉動經濟做出卓越的貢獻, 然而, 由于工程施工周期長, 風險系數高及地質多樣性、施工技術限制和過程規范性不足等因素會導致隧道結構發生形變, 而累積形變則會影響結構安全, 發生重大安全事故, 因此需要對施工期隧道的結構進行全方位監測, 以保障人民生命財產安全.

地下環境對通信網絡具有很強的限制, 且存在環境惡劣, 作業空間小的缺陷, 施工期隧道主要依賴高精度傳感器, 但傳感器所收集的數據并不能直接用于分析或安全評估, 主要原因如下: (1) 測量環境存在隨機噪聲; (2)施工現場存在電磁干擾; (3) 傳輸線存在信號衰減; (4) 電源電壓存在波動; (5) 傳感器本身存在硬件誤差. 因而需要對傳感器直接采集到的數據進行過濾和清洗, 去除測量之中存在的噪聲, 盡可能的“還原”出真實值, 以供后續的數據分析工作.

針對傳感器監測數據的去噪方法主要有時域濾波和頻域濾波[1]. 常用時域濾波方法如算數平均濾波、中值濾波、滑動平均濾波等都是設置一個滑動窗口, 并基于數據的統計特性進行平滑, 這些方法雖然原理簡單, 但計算量較大, 且在實時處理中若對數據一無所知則易發生嚴重偏移; 傳統的頻域濾波器 (高通/低通/帶通) 需要把信號從時域變換到頻域再設定相應閾值進行去噪, 這類方法通常無法還原信號序列的細節; 小波理論是時域-頻域濾波的一個新發展, 其克服了傅里葉變換的缺陷, 將信號經過小波變換后得到一系列分解的小波, 再通過分析小波系數對信號進行重構從而達到消噪. 但小波理論在分解過程中的基函數、閾值函數及每一層的閾值選取方面困難較大.

施工隧道結構安全監控數據具有強實時、無先驗知識及存在少量異常突變的特點, 以上方法難以直接使用. 卡爾曼濾波器(Kalman filter, KF)[2]作為一種針對傳感器數據的無監督濾波算法, 通過結合模型的估計值來動態調整觀測數據, 從而剔除隨機干擾、逼近真實值, 達到去除噪聲的目的. 該方法基于系統狀態空間模型, 僅需k -1時 刻的估計值及k 時刻的觀測值即可計算出k 時刻的估計值, 其遞歸特性也彰顯了KF適用于在線實時動態系統的數據處理. 另外, 由于KF算法針對線性系統進行設計, 在非線性系統方面, 多位學者進行了研究, 其中應用最為廣泛的當屬擴展卡爾曼濾波器(extended Kalman filter, EKF)[3]和無跡卡爾曼濾波器 (unscented Kalman filter, UKF)[4]. KF及其擴展算法自被提出起就廣泛應用于定位導航、目標檢測與追蹤、信號處理以及數據融合等多個領域.

1 基本原理與研究現狀

1.1 線性卡爾曼濾波器

1.1.1 基本原理

線性離散系統的系統方程可以描述為:

式(1)為狀態更新方程. xk為k 時刻信號的狀態變量, 即信號的真實值; Fk-1為 從k -1時 刻到k 時刻的轉移矩陣; uk-1和 Bk-1是系統外部已知的控制輸入及其對應的矩陣; wk-1表 示從k -1時 刻到k 時刻系統自身產生的過程噪聲. 式(2)為觀測方程. zk為k 時刻的實際觀測值;Hk為量測矩陣, 表示真實值與觀測值之間的變化關系;vk表 示k 時刻觀測時的隨機噪聲. 其中, KF對過程噪聲和觀測噪聲做了嚴格的假設: wk和vk均為高斯白噪聲,即它們分別服從均值為0, 協方差為Qk和Rk的高斯分布.

KF的原理即為通過不斷地預測和更新校正來去除噪聲、獲取信號的真實值, 算法實現主要可以分為2個過程: 預測和更新.

預測步是直接使用k -1時刻校正過的后驗估計值得到的, 其結果稱為 k時刻的先驗估計, 也可稱為預測值,計算公式為式(3)、式(4), 其中Pk|k-1是先驗誤差協方差.

更新步的主要思想為引入觀測值對預測步中的先驗估計進行修正, 其原理為式(5)-式(7). Kk為卡爾曼增益, 它是在誤差協方差最小的原則上求出來的, 作用為引入觀測值并動態調整估計值, Kk的大小決定了濾波結果的偏向性;則為結合了預測值和觀測值的k 時刻濾波后的真實值, 而 Pk則為校正后的誤差.

1.1.2 參數設置

首先針對系統方程的參數設置, Fk-1聯 系了k -1到k時刻的狀態變化過程, 當監測數據為一維標量值時,可設置此值為1[5]; 在系統沒有外部控制時, uk-1和Bk-1為0; 在傳感器所得觀測值與所求狀態值一致時,Hk為1.

wk-1和 vk的協方差矩陣Qk和Rk分別表示對模型和觀測的信任程度, 它們的值間接影響 Kk, 故對濾波結果精度影響較大; Rk可以從傳感器的精度參數得出, 而Qk則通常根據經驗設置.

1.2 非線性卡爾曼濾波器

非線性離散系統的系統方程可以描述為:

EKF是一種針對非線性系統常用的濾波算法, 其核心思想為通過一階泰勒展開將非線性的函數 f(x)和h(x)展開為線性函數, 從而轉入KF計算流程. 其預測和更新過程如下:

1.3 卡爾曼濾波器在隧道數據處理中的研究

在隧道變形監測方面, KF常常用于變形趨勢的預測. 陳冠宇等人[7]和陳大勇等人[8]直接將KF應用于隧道沉降數據, 實驗證明濾波前后的降噪效果明顯, 數據曲線更為平滑; Han等人[9]針對過程噪聲Q 和觀測噪聲R 結合滑動窗口根據極大似然估計推導出Q 和R 的計算公式, 算法可自適應調整 Q和R 的值, 提高了濾波器的精度和容錯性; 錢建國等人[10]和王賓賓等人[11]對比了經典KF和基于方差補償的自適應KF在沉降數據的去噪效果, 結果表明自適應算法的濾波和預測效果更好; Yi等人[12]針對有色噪聲通過Yule-Walker算法計算出自適應參數, 并借助KF結構提出一種在線去噪算法, 在混凝土結構監測數據上證明了算法的有效性.

本文通過研究經典卡爾曼濾波器, 發現其在抗野值, 噪聲估計以及非線性系統的適應性方面性能較差,從而導致濾波數據發生嚴重偏移, 針對以上問題, 提出了一種自適應噪聲優化的可靠卡爾曼算法, 通過對野值數據的合理補償及對噪聲的實時估計, 進一步優化了系統的建模精度, 獲得更準確的濾波結果.

2 可靠卡爾曼過濾算法設計

2.1 野值處理

在施工掘進過程中, 傳感器的傳輸線不可避免地會被破壞, 系統將此時的測量值定義為與正常測量范圍有巨大偏差的數值, 稱這樣的值為野值[13]. 在一個相對穩定的環境, KF可獲得很好的去噪效果, 但若系統中存在野值, 會嚴重影響濾波效果, 甚至導致濾波發散,從而影響后續數據分析.

2.1.1 野值檢測

根據野值的定義, 其檢測準則歸結到底即為合適的閾值設計. 常用的野值檢測方法如奈爾準則[14]等都需要提前獲知整個數據列的統計特性, 不適用于動態系統; 而根據統計特性設計的SD、MAD和IQR方法[15]在野值嚴重偏離正常數據時, 會直接影響計算結果從而使閾值失效.

在KF中, 將觀測值與預測值的差值稱為新息, 其定義為式(15). 觀測野值也正是通過影響新息值從而影響卡爾曼增益, 并最終使濾波結果出現巨大偏差, 因而從新息值角度選取合適閾值進行野值的檢測與判斷.

由于土壓力數據中野值與正常測量值差距懸殊,根據工程制定的報警閾值和傳感器自身量程, 當300kPa時, 判斷為野值.

2.1.2 野值補償

在判斷出野值后, 一種簡單的想法是直接用KF濾波值代替, 但這種方法存在著以下問題: (1)直接使用預測值, 使得觀測值喪失意義, 若野值連續長時間出現甚至可能導致濾波發散; (2) KF雖然能在一定程度上“拉回”野值, 但其作用有限, 使用預測值在后續工作中依然無法消除異常值帶來的影響.

由于數據帶有明顯的趨勢性, 常用的填充方法不再適用[16], 本文依數據特點選取結合滑動窗口的最小二乘法進行野值補償.

最小二乘法是一種經常用于工程領域的曲線擬合方法, 其根本思想為遵循最小化誤差平方原則(式(16)), 根據數據特點選取合適的基函數, 構造合適的擬合函數, 并代入已有的觀測數據并求解該函數的參數,從而得到最佳擬合方程. 其中 yi為已知的觀測數據,f(xi)為最佳擬合函數上的計算值.

通過分析數據可以發現: 如圖1所示, 在局部范圍內, 數據變化的趨勢基本為線性. 故本文采用線性最小二乘法進行計算, 選取基函數為1 ,x,x2,···,xn, 階數為1.

圖1 部分原始數據趨勢情況

由于系統為實時采集的動態系統, 且傳輸線在損壞后修復需要一定的時間, 故野值通常連續出現, 鑒于此, 使用滑動窗口進行局部擬合數據, 以更好地適應數據趨勢. 設置滑動窗口初始大小為0, 最大值為11, 結合滑動窗口的野值檢測與補償算法如算法1.

算法1. 野值檢測與補償算法zk size(window)1) 讀入 , 計算 .k ︿ek ︿ek＞300 2) 計算 時刻的新息值 , 若 則轉入3), 否則轉入6).zk size(window)=11 3) 是野值, 若 則轉入4), 否則轉入5).zk-10 zk-1 f(x)=akx+bk zk=f(11)4) 將 至 作為最小二乘的參數, 擬合函數為 ,, 轉入6).zk-size(window)+1 zk-1 f(x)=akx+bk zk=f(k)5) 將至作為最小二乘的參數, 擬合函數為 ,轉入6).size(window)=11 6) 若 , 窗口向后滑動, 否則轉入7).k 7) 若 為最后一個檢測時刻, 則結束算法, 否則轉入8).k++8) , 轉入1)進入下一輪迭代.

2.2 自適應噪聲的卡爾曼濾波器

由第1節可知, 濾波后真實值的校正信息主要由卡爾曼增益決定, 即歸根到底是由系統中的噪聲決定的, 在KF中, 不僅假定噪聲為高斯白噪聲, 而且需要在初始化時對噪聲協方差進行賦值. Q k 和 Rk的取值直接影響濾波效果, 若直接定值為常數, 存在以下缺陷:(1)實際工程往往無法得知較為準確的噪聲值; (2)動態系統的噪聲并非固定不變.

由此提出改進的卡爾曼濾波器, 遞推過程如下:

算法從自適應噪聲變化規律角度提出改進, Ek為第 k次測量的誤差, 即通過增益矯正過的后驗誤差;δ1(k)為預測誤差, 系統動態估計過程中, 過程噪聲主要由硬件誤差引起, 假設噪聲為白噪聲, 則噪聲協方差Qk可以描述為( xk-1-xk-2)2, 則第k 次的預測誤差構造為第k -1次 測量的誤差值向量與第k 次的噪聲向量的和;δ2(k) 為 測量噪聲, 在第k 次測量時, 其系統累積的測量噪聲可直接由第 k-1次測量的觀測值和真實值得出;Kk和 xk分別為校正增益和真實值, 與KF計算思想一脈相承.

在初值設定方面, 設 x0=z1, x-1=z-1=0, E0=0.

該可靠卡爾曼算法 (reliable Kalman filter, RKF)摒棄了傳統算法中在初始化時對噪聲的常數式假定,從噪聲產生的物理意義出發, Qk和Rk直接由上一時刻結果得到, 遵循KF的“一步調整”思想, 自適應地調整系統噪聲, 從而達到更精確的建模和估計.

3 實驗分析

本文數據來源為北京某在建地鐵線, 種類為拱肩土壓力, 共6 380條, 數據有以下特點: (1)前期為斷面區域施工期, 后期為基本穩定期, 整體呈先上升后下降最后逐漸平穩; (2)存在極少量野值, 占數據總量≤0.3%. 由于數據量較大, 圖形不易辨別, 取前100條數據進行研究, 其中第29-44條為野值.

3.1 野值對濾波結果的影響

在傳輸線被破壞時, 傳感器收集的數據顯示為1 500 kPa, 為極度偏離正常值的野值. 雖然野值出現的概率非常低, 但若不加處理直接使用KF, 觀測值對濾波結果的影響依然很大. 經實驗, 原始數據直接使用KF, 此時RMSE高達195.484 6; 由于野值數量極低, 若直接刪除, RMSE則驟降至4.989 0, 彰顯了野值處理的重要性.

為了驗證本文所提野值補償方法的有效性, 首先將第11至20條數據設置為干擾野值(1 500 kPa), 分別對其使用本文方法及常用的均值補償方法, 結果對比如表1所示, 使用滑動最小二乘法補償后數據的最大誤差為3.20, 而傳統均值補償后最大誤差為10.40, 這也證明了統計特性方式的補償方法難以用于趨勢性數據.

表1 不同方法野值補償情況(kPa)

由此, 將此方法應用于實際野值數據, 圖2為補償數據前后的數據對比. 圖2(a)為補償效果圖, 補償前后數據的方差從171 483.809 6下降至197.485 3, 數據趨勢明顯更加平穩; 圖2(b)中的實線即為補償后的數據,可以看出數據整體呈先上升后下降, 且局部還有不少波動, 非線性特征明顯. 在補償效果具體細節方面, 表2為野值補償前后的數據值, 可以看出, 結合滑動窗口的最小二乘擬合補償的數據基本符合原始數據段遞增的變化規律. 從工程角度來看, 野值出現的原因是傳感器傳輸線被破壞, 即施工操作較劇烈, 經與工程進度記錄相比較, 此時間段內斷面挖掘頻繁, 土壓力理應愈來愈大, 與補償值趨勢一致.

表2 野值補償詳情(kPa)

圖2 野值補償前后數據對比

遵循第1節中初值的設置規范以及KF的基本原理可以得到: 在處理一維數據時, Qk和Rk退化為常數Q ,R. KF的濾波結果與初值設置息息相關. KF需要的初值有4個: P0, Q, R 和 x0, 其中, x0通常設置為第一個監測值[5], 即 x0=358; P0的大小對結果影響不大, 故采用常規方式將其設置為平穩段數據(第3968-6380號數據)的方差值, 值為8.222 6; Q和R 為過程噪聲和測量噪聲, 為可調超參數. 表3為調節 Q 值時的濾波結果對比,其中計算所用到的R 為9: 當用于工程項目時, R通常由儀器本身決定, 根據傳感器出廠說明可知土壓力計測量范圍為1 000 kPa, 精度等級為0.3%, 即可計算出測量噪聲協方差為9. Q為系統噪聲協方差, 即建立的狀態空間模型和實際值的差距. 由表3數據變化情況可以看出, Q越小, 濾波后曲線越平滑, 距離測量值越遠;Q 越大, 則越接近觀測曲線. 在Q 值的選取上通常依賴于經驗值, 本系統為標量輸入, 即狀態轉換過程確定,此時Q 的取值越小越好, 這樣可以保證結果收斂快; 而若 Q 值過大時, 則有可能收斂過慢甚至發散. 由于在土壓力濾波方面目前國內外研究情況較少, 且土體運動情況的局限性較強, Q的取值沒有過多可參考的經驗值, 本文采用常用的標量卡爾曼 Q 值設定經驗, 將Q 值固定在[0.01, 1]區間, 并將MAE作為重要參考進行Q值選取: MAE為平均絕對誤差, 若MAE<儀器絕對誤差, 則選取此值作為 Q 值. 由儀器精度等級可以得出絕對誤差為3 kPa, 即選取MAE<3 kPa時的 Q值作為過程噪聲協方差, 由表3可得 Q =0.55. 另外, 為了驗證測量噪聲 R對濾波結果的影響, 將R 作為變量進行調參實驗, 結果為表4, 可以看出R 對濾波結果的調控與Q 相反: R 越小越接近觀測值, R越大則越信任估計值.

表3 Q值對濾波結果的影響

表4 R值對濾波結果的影響

圖3為處理野值前后使用KF的對比. 圖3(a)為直接使用KF時的濾波效果, 可以看出含野值的數據對濾波結果的影響巨大, 未經野值處理直接使用KF時,檢測值遠遠超過正常的400 kPa. 另外, 由于KF的估計值與觀測值相互交織的特性, 使得野值造成的巨大誤差會長久地在系統中累積, 影響時間較長. 通過分析整體6 380條數據發現: 第29-44條數據為野值, 而僅16個野值在此之后持續影響近200條數據才通過觀測值將真實值“拉回”; 圖3(b)為最小二乘法補償后使用KF的濾波結果, 可以看出, KF濾波后數據明顯平緩, 達到去噪效果, 但同時其在局部細節方面表現很差, 例如第40-100條數據明顯的趨勢波動在濾波后直接被平滑,抹去了變化痕跡, 與實際土體運動情況差距過大.

圖3 野值補償前后KF濾波效果

雖然整體效果不算理想, 但在濾波效果方面, 直接使用KF濾波時MAE = 80.1058, MSE = 38214.2423,RMSE = 195.4846, 而補償野值濾波后MAE = 2.9761,MSE = 15.8522, RMSE = 3.9815, 有了明顯的提升.

3.2 不同濾波器對濾波結果的影響

由圖3(b)可以看出, 由于原始數據的非線性較強,濾波后數據與原始數據差異較大, 而這在結構安全監測方面有較大的影響, KF有可能將本應該出現報警的值平滑為安全值, 從而導致嚴重的后果; 也有可能將未達到告警閾值的數據平滑至報警邊緣, 從而誤報警.

分別使用KF, EKF和本文所提的RKF對數據進行濾波實驗. 其中EKF實驗中使用Matlab進行傅里葉逼近進行擬合, 擬合函數為:f(x)=377.8-15.65cos(0.08265x)-11.01sin(0.08265x). 擬合后與原始數據的對比如圖4所示, 模型擬合優度R2為0.931 8.

圖4 傅里葉擬合情況

使用不同濾波器的濾波效果如圖5所示, 具體結果如表5所示. 結合圖表可以看出, 3種算法濾波后數據都較濾波前更加平滑, 即都具備降噪能力, 而濾波精度則可由濾波值與觀測值的接近程度來描述. 由于數據段為非線性, KF濾波結果與實測值差距過大, 尤其在反復波動的情況下, KF無法跟隨信號細節, 這也突出了KF在處理非線性問題時的缺陷; EKF在高度擬合的情況下, 大體上能夠跟隨真實數據變化趨勢, 但由于系統非線性程度較高, 在僅進行一次泰勒展開的情況下, RMSE依然較大, 說明EKF在面對強非線性系統時的不足. 實際上, 高精度擬合在工程中難以實現,從圖4中可以看出, 即使是在高精度擬合的情況下, 擬合函數依然很難跟隨觀測曲線, 第40-80條數據中的3次劇烈變化都被抹去, 而這可能在力學分析中有重要的意義; 另外, 整體數據變化趨勢伴隨施工過程的進行而呈明顯分段, 顯然不能用同一個函數擬合所有數據,這也大大增加了EKF使用的難度; 而RKF直接動態跟蹤噪聲, 突破了強非線性的限制, 從圖5中也可以看出RKF跟蹤觀測曲線的效果最好, 即精度最高. 表5是3種濾波器的濾波效果對比, 最終結果顯示: 由于系統的強非線性, KF和EKF的濾波精度不夠高; 而對于RKF來說, 即使是在強非線性環境, 算法依然可以得到很好的濾波效果, 比KF和EKF更加接近觀測值, 其RMSE = 0.706 8, 比EKF提高了74.53%, 進一步提高了濾波精度.

圖5 不同濾波器濾波結果對比

表5 不同濾波器濾波結果對比

4 結論與展望

本文在經典卡爾曼濾波器的基礎上, 針對其無法有效處理野值及應用于強非線性系統的缺陷, 從實時動態跟蹤系統噪聲的角度出發, 將新息值作為野值檢測指標, 結合滑動窗口使用線性最小二乘方法進行野值補償, 并將原本需要確切知曉并以常數表示的噪聲Qk和Rk用上一時刻的硬件誤差及測量誤差作為替代,通過融入最終誤差進行分析和計算, 實現了系統噪聲的自適應估計, 從而解決野值補償及精確追蹤信號細節的問題. 經實際工程數據的實驗分析可知, 本文算法較經典算法有了明顯的提升, 具有一定的實用價值.