劉興亮,冉 典,余學(xué)祥
(1.安徽理工大學(xué) 測繪學(xué)院,安徽 淮南 232001)
自適應(yīng)Kalman濾波在地表沉降觀測中的應(yīng)用
劉興亮1,冉 典1,余學(xué)祥1
(1.安徽理工大學(xué) 測繪學(xué)院,安徽 淮南 232001)
地下采煤引起的地表沉陷是一個時間和空間的過程,據(jù)此提出了觀測站動態(tài)數(shù)據(jù)處理模型Kalman濾波和自適應(yīng)Kalman濾波,通過實例驗證了自適應(yīng)Kalman濾波比普通Kalman濾波在觀測站數(shù)據(jù)濾波和預(yù)測中具有優(yōu)越性。
動態(tài)數(shù)據(jù)處理;Kalman濾波;自適應(yīng)Kalman濾波
對地表移動觀測站數(shù)據(jù)建立預(yù)測模型的方法很多,如時間序列方法、頻譜分析方法、Kalman濾波方法等。前2種方法要求觀測數(shù)據(jù)較多,并且只適用于對一點進行預(yù)測,而Kalman濾波方法對初始數(shù)據(jù)個數(shù)要求較少,且適用于對整個觀測站或其中任一監(jiān)測點進行預(yù)測。下面以某地表移動觀測站為例,介紹Kalman濾波[1]和自適應(yīng)Kalman濾波[2]預(yù)測模型建立的原理和步驟。
假設(shè)觀測站由n個監(jiān)測點組成,監(jiān)測網(wǎng)高差觀測值個數(shù)為m。以測點的高程和高程速率為狀態(tài)向量。設(shè)點i在時刻t的位置向量為ξi(t),瞬時速率為λi(t),將瞬時加速率Ω(t)看成隨機干擾項,有以下微分關(guān)系式成立:
記i點的狀態(tài)向量為Xi(t),全網(wǎng)的狀態(tài)向量為X(t),則式(1)可寫成:
式(2)是一個常系數(shù)連續(xù)線性系統(tǒng)微分方程,采用拉普拉斯變換將其離散化[3],可得Kalman濾波的狀態(tài)方程為:
根據(jù)文獻[1]可得全網(wǎng)n個點的狀態(tài)方程為:
式(4)即為觀測站Kalman濾波的狀態(tài)方程。
將整個觀測站看作一個動態(tài)監(jiān)測系統(tǒng),而動態(tài)監(jiān)測系統(tǒng)由狀態(tài)方程和觀測方程[4]組成。當(dāng)以高差為觀測值時,某一觀測值Lij在第k+1次的觀測方程為:
Lij/k+1=-ξi/k+1+ξj/k+1-?tij/k+1λi/k+1+?tij/k+1λj/k+1+νij/k+1(5)式(5)中,?tij/k+1=tij/k+1-tk+1,而tij/k+1為高差Lij的觀測時刻,tk+1為第k+1次觀測時各高差觀測的中心時刻。對于移動觀測站而言,一般每1個月左右觀測一次(沉陷活躍期半個月左右觀測一次),?tij/k+1與兩次觀測的時間間隔?tk=tk+1-tk相比可忽略不記,則式(5)可以表示為:
整個觀測站的觀測方程[5]為:
式(4)和式(7)共同構(gòu)成了整個觀測站動態(tài)數(shù)據(jù)處理的Kalman濾波模型:
式中,Φk,k-1為k-1到k時刻的系統(tǒng)一步轉(zhuǎn)移矩陣;Γk,k-1為系統(tǒng)噪聲矩陣;Ωk-1為k-1時刻的系統(tǒng)噪聲;Bk為k時刻系統(tǒng)的觀測矩陣;Δk為k時刻系統(tǒng)的觀測噪聲;Xk為k時刻的系統(tǒng)待估狀態(tài)參數(shù);Lk為k時刻系統(tǒng)的觀測向量矩陣。
觀測站Kalman濾波模型由函數(shù)模型和隨機模型[6]組成,式(8)為函數(shù)模型,隨機模型為:
式中,DΩ(k)稱為系統(tǒng)動態(tài)噪聲方差陣;δkj為
根據(jù)廣義最小二乘估計原理,可利用高差觀測值Lk計算tk時刻狀態(tài)向量Xk的最佳估值X^(k,k)。用X(k,k)表示X^(k,k)。Kalman濾波方程為:
式中,
式中,X(k/k-1)為一步預(yù)測值;DX(k/k-1)為一步預(yù)測方差陣;Jk為狀態(tài)增益矩陣;Ek為預(yù)測新息或殘差。
一般的Kalman濾波器其預(yù)測結(jié)果容易發(fā)散,為克服發(fā)散現(xiàn)象,常采用自適應(yīng)濾波[7]方法。其構(gòu)建思想是:假定觀測值不存在粗差,存在標(biāo)準(zhǔn)的正態(tài)分布,當(dāng)動力模型存在異常誤差時,將動力模型信息作為一個整體,采用統(tǒng)一自適應(yīng)因子調(diào)整動力學(xué)模型信息對狀態(tài)參數(shù)的影響。按求條件極值的方法所確定的自適應(yīng)Kalman濾波公式如下:
式中,
Jk為新的增益矩陣:αk為自適應(yīng)調(diào)節(jié)因子:
式中,c0可取1.0~1.5;c1可取3.0~8.0。
式中,Xk為當(dāng)前參數(shù)的最小二乘解:
從自適應(yīng)濾波表達式可看出,自適應(yīng)Kalman濾波改進了標(biāo)準(zhǔn)Kalman濾波。隨著αk的變化,自適應(yīng)Kalman濾波在最小二乘、標(biāo)準(zhǔn)Kalman濾波和自身之間變化。
下面以某礦區(qū)的地表移動觀測站為例進行說明。選取走向線上具有典型代表性的2個監(jiān)測點ML05和ML16作為研究對象。ML05離開采工作面邊界較近,受開采影響較小,ML16在下沉曲線的拐點附近,其水平移動和傾斜最明顯。表1和表2分別給出了監(jiān)測點ML05和ML16采用最小二乘、標(biāo)準(zhǔn)Kalman濾波和自適應(yīng)Kalman濾波時的高程值及其差值;圖1和圖2為ML05和ML16高程平差值和Kalman及自適應(yīng)Kalman濾波值差值比較;表3和表4給出了平差值和采用標(biāo)準(zhǔn)Kalman濾波及自適應(yīng)Kalman濾波預(yù)測值;圖3、圖4為監(jiān)測點ML05、ML16高程平差值和標(biāo)準(zhǔn)Kalman、自適應(yīng)Kalman濾波預(yù)測值對比。這里給出了從第5期到第12期的結(jié)果,前4期主要用來確定濾波初始值。
表1 監(jiān)測點ML05高程濾波值比較
圖1 ML05點高程平差值和Kalman及自適應(yīng)Kalman濾波值差值比較
表2 監(jiān)測點ML16高程濾波值比較
圖2 ML16點高程平差值和Kalman及自適應(yīng)Kalman濾波值差值比較
表3 監(jiān)測點ML05高程預(yù)測值比較
圖3 監(jiān)測點ML05高程平差值和標(biāo)準(zhǔn)Kalman、
自適應(yīng)Kalman濾波預(yù)測值對比表4 監(jiān)測點ML16高程預(yù)測值比較
圖4 監(jiān)測點ML16高程平差值和標(biāo)準(zhǔn)Kalman、自適應(yīng)Kalman濾波預(yù)測值對比
對以上數(shù)據(jù)及圖表分析可得到以下結(jié)論:
1)對2個監(jiān)測點總共11期的觀測數(shù)據(jù),標(biāo)準(zhǔn)Kalman濾波前幾期濾波效果較差,后幾期濾波效果較好;而自適應(yīng)Kalman濾波在11期數(shù)據(jù)處理中濾波效果均很好,通過濾波差值中誤差可以看出自適應(yīng)Kalman濾波比標(biāo)準(zhǔn)Kalman濾波精度高。
2)對2個監(jiān)測點高程預(yù)測,標(biāo)準(zhǔn)卡爾曼濾波和自適應(yīng)卡爾曼濾波預(yù)測精度相當(dāng)。但對ML16點的預(yù)測精度低于ML05點,主要是ML16點受開采影響較大。ML16點第9期的預(yù)測偏差較大,可能是受地表突變影響。
[1] 崔希璋,於宗儔,陶本藻,等.廣義測量平差[M].武漢:武漢測繪科技大學(xué)出版社,2001
[2] 楊元喜,何海波,徐天河.論動態(tài)自適應(yīng)濾波[J].測繪學(xué)報,2001,30(4):293-298
[3] 馬攀,文鴻雁.離散卡爾曼濾波用于GPS動態(tài)變形數(shù)據(jù)處理[J].西安科技大學(xué)學(xué)報,2006,26(3):353-357
[4] 梅連友.卡爾曼濾波在滑坡監(jiān)測中的應(yīng)用[J].測繪工程,2004(3):13-15
[5] 余學(xué)祥,張華海.Kalman濾波在GPS監(jiān)測網(wǎng)中的應(yīng)用[J].工程勘察,2000(4):33-35
[6] 鄧自立.卡爾曼濾波與維納濾波[M].哈爾濱:哈爾濱工業(yè)大學(xué)出版社,2001
[7] 張海,常艷紅,車歡.基于GPS不同測量特性的自適應(yīng)卡爾曼濾波算法[J].中國慣性技術(shù)學(xué)報,2010,18(6):676-701
P258
B
1672-4623(2014)06-0091-03
10.3969/j.issn.1672-4623.2014.06.032
劉興亮,碩士,從事大地測量學(xué)與測繪工程研究。
2013-12-19。
項目來源:安徽高校省級自然科學(xué)研究重點資助項目(KJ2010A104);安徽省國土資源科技資助項目(2011-K-18)。