自適應(yīng)Kalman濾波在地表沉降觀測(cè)中的應(yīng)用

劉興亮，冉 典，余學(xué)祥

（1.安徽理工大學(xué) 測(cè)繪學(xué)院，安徽 淮南 232001）

自適應(yīng)Kalman濾波在地表沉降觀測(cè)中的應(yīng)用

劉興亮1，冉 典1，余學(xué)祥1

（1.安徽理工大學(xué) 測(cè)繪學(xué)院，安徽 淮南 232001）

地下采煤引起的地表沉陷是一個(gè)時(shí)間和空間的過(guò)程，據(jù)此提出了觀測(cè)站動(dòng)態(tài)數(shù)據(jù)處理模型Kalman濾波和自適應(yīng)Kalman濾波，通過(guò)實(shí)例驗(yàn)證了自適應(yīng)Kalman濾波比普通Kalman濾波在觀測(cè)站數(shù)據(jù)濾波和預(yù)測(cè)中具有優(yōu)越性。

動(dòng)態(tài)數(shù)據(jù)處理；Kalman濾波；自適應(yīng)Kalman濾波

對(duì)地表移動(dòng)觀測(cè)站數(shù)據(jù)建立預(yù)測(cè)模型的方法很多，如時(shí)間序列方法、頻譜分析方法、Kalman濾波方法等。前2種方法要求觀測(cè)數(shù)據(jù)較多，并且只適用于對(duì)一點(diǎn)進(jìn)行預(yù)測(cè)，而Kalman濾波方法對(duì)初始數(shù)據(jù)個(gè)數(shù)要求較少，且適用于對(duì)整個(gè)觀測(cè)站或其中任一監(jiān)測(cè)點(diǎn)進(jìn)行預(yù)測(cè)。下面以某地表移動(dòng)觀測(cè)站為例，介紹Kalman濾波[1]和自適應(yīng)Kalman濾波[2]預(yù)測(cè)模型建立的原理和步驟。

1 Kalman濾波狀態(tài)方程和觀測(cè)方程

假設(shè)觀測(cè)站由n個(gè)監(jiān)測(cè)點(diǎn)組成，監(jiān)測(cè)網(wǎng)高差觀測(cè)值個(gè)數(shù)為m。以測(cè)點(diǎn)的高程和高程速率為狀態(tài)向量。設(shè)點(diǎn)i在時(shí)刻t的位置向量為ξi（t），瞬時(shí)速率為λi（t），將瞬時(shí)加速率Ω（t）看成隨機(jī)干擾項(xiàng)，有以下微分關(guān)系式成立：

記i點(diǎn)的狀態(tài)向量為Xi（t），全網(wǎng)的狀態(tài)向量為X（t），則式（1）可寫(xiě)成：

式（2）是一個(gè)常系數(shù)連續(xù)線性系統(tǒng)微分方程，采用拉普拉斯變換將其離散化[3]，可得Kalman濾波的狀態(tài)方程為：

根據(jù)文獻(xiàn)[1]可得全網(wǎng)n個(gè)點(diǎn)的狀態(tài)方程為：

式（4）即為觀測(cè)站Kalman濾波的狀態(tài)方程。

將整個(gè)觀測(cè)站看作一個(gè)動(dòng)態(tài)監(jiān)測(cè)系統(tǒng)，而動(dòng)態(tài)監(jiān)測(cè)系統(tǒng)由狀態(tài)方程和觀測(cè)方程[4]組成。當(dāng)以高差為觀測(cè)值時(shí)，某一觀測(cè)值Lij在第k+1次的觀測(cè)方程為：

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的觀測(cè)時(shí)刻，tk+1為第k+1次觀測(cè)時(shí)各高差觀測(cè)的中心時(shí)刻。對(duì)于移動(dòng)觀測(cè)站而言，一般每1個(gè)月左右觀測(cè)一次（沉陷活躍期半個(gè)月左右觀測(cè)一次），?tij/k+1與兩次觀測(cè)的時(shí)間間隔?tk=tk+1-tk相比可忽略不記，則式（5）可以表示為：

整個(gè)觀測(cè)站的觀測(cè)方程[5]為：

式（4）和式（7）共同構(gòu)成了整個(gè)觀測(cè)站動(dòng)態(tài)數(shù)據(jù)處理的Kalman濾波模型：

式中，Φk,k-1為k-1到k時(shí)刻的系統(tǒng)一步轉(zhuǎn)移矩陣；Γk,k-1為系統(tǒng)噪聲矩陣；Ωk-1為k-1時(shí)刻的系統(tǒng)噪聲；Bk為k時(shí)刻系統(tǒng)的觀測(cè)矩陣；Δk為k時(shí)刻系統(tǒng)的觀測(cè)噪聲；Xk為k時(shí)刻的系統(tǒng)待估狀態(tài)參數(shù)；Lk為k時(shí)刻系統(tǒng)的觀測(cè)向量矩陣。

2 觀測(cè)站Kalman濾波模型

觀測(cè)站Kalman濾波模型由函數(shù)模型和隨機(jī)模型[6]組成，式（8）為函數(shù)模型，隨機(jī)模型為：

式中，DΩ（k）稱為系統(tǒng)動(dòng)態(tài)噪聲方差陣；δkj為

根據(jù)廣義最小二乘估計(jì)原理，可利用高差觀測(cè)值Lk計(jì)算tk時(shí)刻狀態(tài)向量Xk的最佳估值X^（k，k）。用X（k，k）表示X^（k，k）。Kalman濾波方程為：

式中，

式中，X（k/k-1）為一步預(yù)測(cè)值；DX（k/k-1）為一步預(yù)測(cè)方差陣；Jk為狀態(tài)增益矩陣；Ek為預(yù)測(cè)新息或殘差。

3 自適應(yīng)Kalman濾波

一般的Kalman濾波器其預(yù)測(cè)結(jié)果容易發(fā)散，為克服發(fā)散現(xiàn)象，常采用自適應(yīng)濾波[7]方法。其構(gòu)建思想是：假定觀測(cè)值不存在粗差，存在標(biāo)準(zhǔn)的正態(tài)分布，當(dāng)動(dòng)力模型存在異常誤差時(shí)，將動(dòng)力模型信息作為一個(gè)整體，采用統(tǒng)一自適應(yīng)因子調(diào)整動(dòng)力學(xué)模型信息對(duì)狀態(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)濾波表達(dá)式可看出，自適應(yīng)Kalman濾波改進(jìn)了標(biāo)準(zhǔn)Kalman濾波。隨著αk的變化，自適應(yīng)Kalman濾波在最小二乘、標(biāo)準(zhǔn)Kalman濾波和自身之間變化。

4 實(shí)例分析

下面以某礦區(qū)的地表移動(dòng)觀測(cè)站為例進(jìn)行說(shuō)明。選取走向線上具有典型代表性的2個(gè)監(jiān)測(cè)點(diǎn)ML05和ML16作為研究對(duì)象。ML05離開(kāi)采工作面邊界較近，受開(kāi)采影響較小，ML16在下沉曲線的拐點(diǎn)附近，其水平移動(dòng)和傾斜最明顯。表1和表2分別給出了監(jiān)測(cè)點(diǎn)ML05和ML16采用最小二乘、標(biāo)準(zhǔn)Kalman濾波和自適應(yīng)Kalman濾波時(shí)的高程值及其差值；圖1和圖2為ML05和ML16高程平差值和Kalman及自適應(yīng)Kalman濾波值差值比較；表3和表4給出了平差值和采用標(biāo)準(zhǔn)Kalman濾波及自適應(yīng)Kalman濾波預(yù)測(cè)值；圖3、圖4為監(jiān)測(cè)點(diǎn)ML05、ML16高程平差值和標(biāo)準(zhǔn)Kalman、自適應(yīng)Kalman濾波預(yù)測(cè)值對(duì)比。這里給出了從第5期到第12期的結(jié)果，前4期主要用來(lái)確定濾波初始值。

表1 監(jiān)測(cè)點(diǎn)ML05高程濾波值比較

圖1 ML05點(diǎn)高程平差值和Kalman及自適應(yīng)Kalman濾波值差值比較

表2 監(jiān)測(cè)點(diǎn)ML16高程濾波值比較

圖2 ML16點(diǎn)高程平差值和Kalman及自適應(yīng)Kalman濾波值差值比較

表3 監(jiān)測(cè)點(diǎn)ML05高程預(yù)測(cè)值比較

圖3 監(jiān)測(cè)點(diǎn)ML05高程平差值和標(biāo)準(zhǔn)Kalman、

自適應(yīng)Kalman濾波預(yù)測(cè)值對(duì)比表4 監(jiān)測(cè)點(diǎn)ML16高程預(yù)測(cè)值比較

圖4 監(jiān)測(cè)點(diǎn)ML16高程平差值和標(biāo)準(zhǔn)Kalman、自適應(yīng)Kalman濾波預(yù)測(cè)值對(duì)比

5 結(jié) 語(yǔ)

對(duì)以上數(shù)據(jù)及圖表分析可得到以下結(jié)論：

1）對(duì)2個(gè)監(jiān)測(cè)點(diǎn)總共11期的觀測(cè)數(shù)據(jù)，標(biāo)準(zhǔn)Kalman濾波前幾期濾波效果較差，后幾期濾波效果較好；而自適應(yīng)Kalman濾波在11期數(shù)據(jù)處理中濾波效果均很好，通過(guò)濾波差值中誤差可以看出自適應(yīng)Kalman濾波比標(biāo)準(zhǔn)Kalman濾波精度高。

2）對(duì)2個(gè)監(jiān)測(cè)點(diǎn)高程預(yù)測(cè)，標(biāo)準(zhǔn)卡爾曼濾波和自適應(yīng)卡爾曼濾波預(yù)測(cè)精度相當(dāng)。但對(duì)ML16點(diǎn)的預(yù)測(cè)精度低于ML05點(diǎn)，主要是ML16點(diǎn)受開(kāi)采影響較大。ML16點(diǎn)第9期的預(yù)測(cè)偏差較大，可能是受地表突變影響。

[1] 崔希璋,於宗儔,陶本藻,等．廣義測(cè)量平差[M]．武漢:武漢測(cè)繪科技大學(xué)出版社,2001

[2] 楊元喜,何海波,徐天河．論動(dòng)態(tài)自適應(yīng)濾波[J]．測(cè)繪學(xué)報(bào),2001,30(4):293-298

[3] 馬攀,文鴻雁．離散卡爾曼濾波用于GPS動(dòng)態(tài)變形數(shù)據(jù)處理[J]．西安科技大學(xué)學(xué)報(bào),2006,26(3):353-357

[4] 梅連友．卡爾曼濾波在滑坡監(jiān)測(cè)中的應(yīng)用[J]．測(cè)繪工程,2004(3):13-15

[5] 余學(xué)祥,張華海．Kalman濾波在GPS監(jiān)測(cè)網(wǎng)中的應(yīng)用[J]．工程勘察,2000(4):33-35

[6] 鄧自立．卡爾曼濾波與維納濾波[M]．哈爾濱:哈爾濱工業(yè)大學(xué)出版社,2001

[7] 張海,常艷紅,車歡．基于GPS不同測(cè)量特性的自適應(yīng)卡爾曼濾波算法[J]．中國(guó)慣性技術(shù)學(xué)報(bào),2010,18(6):676-701

P258

B

1672-4623（2014）06-0091-03

10.3969/j.issn.1672-4623.2014.06.032

劉興亮，碩士，從事大地測(cè)量學(xué)與測(cè)繪工程研究。

2013-12-19。

項(xiàng)目來(lái)源：安徽高校省級(jí)自然科學(xué)研究重點(diǎn)資助項(xiàng)目（KJ2010A104）；安徽省國(guó)土資源科技資助項(xiàng)目（2011-K-18）。