徐愛功,張建龍,隋 心,劉 偉,張兆南

(1.遼寧工程技術(shù)大學(xué) 測繪與地理科學(xué)學(xué)院，遼寧 阜新 123000；2.洛陽電子裝備試驗(yàn)中心，河南 洛陽 471000；3.東北煤田地質(zhì)局物探測量隊(duì)，遼寧 沈陽 110101)

(1.School of Geomatics, Liaoning Technical University, Fuxin 123000, China; 2.Luoyang Electronic Equipment Test Center, Luoyang 471000, China; 3.Physical Exploration and Surveying Team, Northeast Coalfield Geology Bureau, Shenyang 110101, China)

基于自回歸法的滑坡監(jiān)測數(shù)據(jù)Kalman濾波的應(yīng)用

徐愛功1,張建龍1,隋 心1,劉 偉2,張兆南3

(1.遼寧工程技術(shù)大學(xué) 測繪與地理科學(xué)學(xué)院，遼寧 阜新 123000；2.洛陽電子裝備試驗(yàn)中心，河南 洛陽 471000；3.東北煤田地質(zhì)局物探測量隊(duì)，遼寧 沈陽 110101)

介紹了卡爾曼濾波的基本模型，針對(duì)運(yùn)動(dòng)學(xué)方法以監(jiān)測點(diǎn)位置參數(shù)、速率參數(shù)為狀態(tài)向量，以加速度為噪聲向量建立觀測方程和狀態(tài)方程的過程比較復(fù)雜。把自回歸法引入狀態(tài)方程和觀測方程的建立中，并結(jié)合某露天礦滑坡動(dòng)態(tài)監(jiān)測數(shù)據(jù)分析了模型的應(yīng)用。結(jié)果表明，濾波值和觀測值之間差值保持在1～3 cm之間，濾波后圖像較原觀測值圖像更為光滑，表明所建模型是合適的，能夠反映滑坡動(dòng)態(tài)變化過程，從而為礦滑坡監(jiān)測預(yù)報(bào)提供數(shù)學(xué)工具。

自回歸模型；卡爾曼濾波；滑坡；預(yù)報(bào)

我國露天礦數(shù)量眾多，高陡邊坡穩(wěn)定性一直是影響露天礦生產(chǎn)安全的重要問題。據(jù)統(tǒng)計(jì)，不穩(wěn)定邊坡或具有潛在滑坡危險(xiǎn)的邊坡占邊坡總長度可達(dá)15%～20%，甚至更高[1]。由于滑坡具有時(shí)間上的突發(fā)性和后果上的嚴(yán)重性，因而加強(qiáng)對(duì)滑坡體的監(jiān)測，對(duì)滑坡作短期、中期和長期預(yù)報(bào)，具有重要的現(xiàn)實(shí)意義[2]。

Kalman濾波是美籍匈牙利數(shù)學(xué)家卡爾曼在濾波理論中引入狀態(tài)空間分析方法提出的一種處理動(dòng)態(tài)數(shù)據(jù)的方法，是一種線性無偏且以均方誤差最小為準(zhǔn)則的最優(yōu)估計(jì)算法,由于其具有可遞推實(shí)現(xiàn)，利于計(jì)算機(jī)運(yùn)算等優(yōu)點(diǎn)，已成為目前處理動(dòng)態(tài)觀測數(shù)據(jù)的一種常用方法?；率且粋€(gè)動(dòng)態(tài)發(fā)展過程，適合采用卡爾曼濾波處理觀測數(shù)據(jù)，并作出預(yù)報(bào)。文獻(xiàn)[3-6]在建立濾波模型時(shí)以監(jiān)測點(diǎn)的位置、速率和加速度等參數(shù)構(gòu)建狀態(tài)向量、觀測向量，這種方法能夠反映變形體的運(yùn)動(dòng)狀態(tài)，但是構(gòu)造過程較為復(fù)雜。本文將監(jiān)測點(diǎn)高程變化看做是一個(gè)依賴于時(shí)間的變量，監(jiān)測前后各期之間具有一定的相關(guān)性，利用自回歸法建立卡爾曼濾波的狀態(tài)方程和觀測方程，并在某露天礦滑坡動(dòng)態(tài)數(shù)據(jù)處理中的應(yīng)用做了研究。

1 Kalman濾波模型

1.1 Kalman濾波模型的建立

Kalman濾波方程的建立，需要已知系統(tǒng)模型和觀測模型，并且要給出噪聲的先驗(yàn)統(tǒng)計(jì)性質(zhì)，以及系統(tǒng)狀態(tài)的初始值。進(jìn)而根據(jù)建立的方程實(shí)時(shí)獲得系統(tǒng)的狀態(tài)變量和輸入信號(hào)的最優(yōu)估值。下面給出離散Kalman濾波模型的建立過程[7]：

狀態(tài)方程xk=Φk,k-1xk-1+Γk,k-1ωk-1.

(1)

觀測方程zk=Hkxk+vk.

(2)

式中:Φk,k-1為狀態(tài)轉(zhuǎn)移矩陣；Γk,k-1為狀態(tài)噪聲系數(shù)陣；Hk為觀測向量系數(shù)陣；ωk-1,vk分別表示狀態(tài)噪聲與觀測噪聲，且具有如下統(tǒng)計(jì)特性：

E(ωk-1)，E(vk)為ωk-1，vk的期望；Rωω(k,j)，Rvv(k,j)為ωk，vk的自協(xié)方差；Rωv(k,j)為ωk，vk的互協(xié)方差；E(x0)，var(x0)為已知的初值；R(x0,ωk-1)，R(x0,vk)為x0與ωk-1和vk的協(xié)方差。由于動(dòng)態(tài)噪聲和觀測噪聲均是零均值的白噪聲序列，且互不相關(guān),故以上建立的模型可以稱為完全不相關(guān)的白噪聲作用下的Kalman濾波。在滿足均方誤差最小的前提下可以得出以下濾波遞推公式[7]：

狀態(tài)估計(jì)校正

(3)

卡爾曼增益

(4)

誤差協(xié)方差預(yù)測

(5)

誤差協(xié)方差估計(jì)校正

Pk|k=(I-KkHk)Pk|k-1.

(6)

圖1 卡爾曼濾波示意圖

1.2 狀態(tài)方程和觀測方程的建立

在處理滑坡監(jiān)測數(shù)據(jù)時(shí)，狀態(tài)方程和觀測方程的建立是關(guān)鍵的一步，盡管Kalman濾波模型狀態(tài)方程和觀測方程的建立方法很多，但多數(shù)文獻(xiàn)是以監(jiān)測點(diǎn)的運(yùn)動(dòng)學(xué)模型為基礎(chǔ)建立的。例如，文獻(xiàn)[3-5]是以監(jiān)測點(diǎn)的位置參數(shù)(水平坐標(biāo)和高程)及三維坐標(biāo)速率作為狀態(tài)向量，以三維坐標(biāo)加速度為噪聲向量建立的。文獻(xiàn)[6]將某時(shí)刻變形量對(duì)上一時(shí)刻展開成泰勒級(jí)數(shù)取至三次方，仍然是采用運(yùn)動(dòng)學(xué)模型。該方法能夠反映變形體的運(yùn)動(dòng)狀態(tài)，但方程建立過程較為復(fù)雜。本文采用一種較為簡單的方法建立狀態(tài)方程和觀測方程，即自回歸法，避免了計(jì)算點(diǎn)位變化過程的速度和加速度，建立過程如下：

自回歸模型

xk=η1xk-1+η2xk-2+…+ηnxk-n+ωk-1.

(7)

η1,η2,…,ηn是n階加權(quán)系數(shù)，ωk-1是白噪聲序列[8]，將該模型引入Kalman濾波中。

令x1(k)=xk，x2(k)=xk-1…

xi(k)=xk-i+1，xn(k)=xk-n+1.

則狀態(tài)方程為

即xk+1=Φk+1,kxk+Γk+1,kωk.

其中:

觀測方程為

zk+1=[1 0 … 0]·[x1(k+1)x2(k+1) …

xn(k+1)]T+v(k+1).

即zk+1=Hk+1xk+1+vk+1.

其中：Hk+1=[1 0 … 0]。

至此，就完成了狀態(tài)方程和觀測方程的建立。

這里，對(duì)加權(quán)系數(shù)ηi(i=1,2,…,n)的求法通常有矩估計(jì)法、最小二乘法、最小方差法、極大似然估計(jì)法等[8]，通常采用最小二乘法，將式(7)改寫成矩陣形式為

即V=Aη-L.

(8)

1.3 狀態(tài)初值的確定

卡爾曼濾波是由一組遞推公式給出的，在濾波前需要給定初值x0和P0。x0可以根據(jù)前面幾組觀測值給定，P0可以由初始測量值經(jīng)過統(tǒng)計(jì)方式給出，通常P0應(yīng)為對(duì)角陣[9]。

2 實(shí)例分析

2.1 概 況

阜新礦業(yè)集團(tuán)白音華某露天礦位于內(nèi)蒙古西烏珠穆沁旗境內(nèi)，在不斷開挖過程中邊坡受采礦爆破影響強(qiáng)烈，加之大型機(jī)械作業(yè)震動(dòng)及長期雨水侵蝕，以及不斷堆積擴(kuò)大的排土場對(duì)滑坡體壓力不斷增大[10]，穩(wěn)定性逐漸變差，最終導(dǎo)致滑坡的發(fā)生。

通過監(jiān)測獲取的數(shù)據(jù)采用Kalman濾波模型處理。濾波模型是通過已知觀測值建立的，模型系數(shù)矩陣在模型確定的前提下即可得到。本文首先采用某監(jiān)測點(diǎn)前30期的高程觀測值(見表1)求解自回歸方程，建立基于自回歸法的Kalman濾波模型，通過分析確定轉(zhuǎn)移矩陣及觀測矩陣，建立狀態(tài)方程和觀測方程。再對(duì)后面10期觀測數(shù)據(jù)進(jìn)行濾波，檢驗(yàn)?zāi)Ｐ偷恼_性及濾波效果。

表1 某監(jiān)測點(diǎn)各期觀測值高程

2.2 模型建立

Γ=[1 0 0]T，H=[1 0 0].

2.3 程序設(shè)計(jì)

由于Matlab具有變量定義簡單、編程效率高、作圖方便等優(yōu)點(diǎn)[11]，故對(duì)上述建立的卡爾曼濾波采用Matlab進(jìn)行程序設(shè)計(jì)，Kalman濾波程序如下：

clc;

close all;

z1=xlsread(′data.xls′,′sheet1′,′A1:A30′); %觀測值

p0=eye(3); %初值方差陣賦初值

p=p0;

phik=[0.4049,0.2994,0.2957;1,0,0;0,1,0]; %狀態(tài)轉(zhuǎn)移矩陣

Qk=eye(3); %狀態(tài)噪聲矩陣

Hk=[1,0,0]; %觀測向量系數(shù)矩陣

Rk=1; %觀測向量噪聲矩陣

x0=[954.615,954.575,954.525]′;

Xk=x0; %狀態(tài)矩陣賦初值

for k=1:30

Pk=phik*p*phik′+Qk; %誤差協(xié)方差

Kk=Pk* Hk′*inv(Hk*Pk*Hk′+Rk); %增益矩陣

Zk=z1(k);

Xk=phik*Xk+Kk*(Zk-Hk*phik*Xk); %最優(yōu)估計(jì)值

x1(k)=Xk(1);

p=(eye(3)-Kk*Hk)*Pk; %誤差協(xié)方差估計(jì)校正

end

2.4 數(shù)據(jù)分析

由圖2、圖3可以看出濾波值和原始測量值保持一致的下降趨勢，且最大偏差為6.28 cm，是由于第16和21兩次觀測值可能含有較大誤差，大部分差值保持在1～4 cm。濾波后圖像比原測量曲線要平滑，說明模型建立的是合理的，能夠達(dá)到消除測量中干擾誤差的目的，經(jīng)過預(yù)測、修正所得到的濾波最優(yōu)估值能夠更好地反映滑坡動(dòng)態(tài)變化過程。

圖2 觀測值與濾波值對(duì)比

圖3 觀測值與濾波值的偏差

通過以上建立的濾波模型再對(duì)未來10期(31～40期)觀測數(shù)據(jù)做出預(yù)測，初值的確定采用以上數(shù)據(jù)最后3期，即x0=[953.873,953.850,953.730]。觀測值與估計(jì)值如表2所示。觀測值與估計(jì)值對(duì)比和差值如圖4、圖5所示?？梢钥闯?，僅預(yù)測初始2期偏差較大，這是由于遞推過程初始值選取與實(shí)際觀測值之間有一定的偏差造成的，之后幾期偏差都在1～3 cm，說明初值的選取僅會(huì)影響前面較少幾期的預(yù)測結(jié)果，而不會(huì)對(duì)之后幾期有影響；濾波估計(jì)值呈一致的下降趨勢，而原始觀測值由于存在觀測誤差并沒有呈一致下降趨勢，說明濾波后能夠較好反映滑坡動(dòng)態(tài)發(fā)展規(guī)律，濾波對(duì)滑坡的預(yù)測效果良好。

表2 未來10期高程觀測值與估計(jì)值對(duì)照

圖4 未來10期觀測值與估計(jì)值對(duì)比

圖5 未來10期觀測值與估計(jì)值的差值

3 結(jié) 論

通過以上模型的建立及數(shù)據(jù)處理結(jié)果分析可以得出以下結(jié)論：

1)將自回歸法引入狀態(tài)方程和觀測方程的建立過程中，使得濾波模型的建立過程簡單化，避免采用運(yùn)動(dòng)學(xué)方法的復(fù)雜性。

2)Kalman濾波把參數(shù)的估計(jì)和預(yù)報(bào)結(jié)合起來，有效消除觀測噪聲，通過參數(shù)的不斷修正，增強(qiáng)模型適應(yīng)不斷獲得的數(shù)據(jù)，可以實(shí)時(shí)計(jì)算濾波值，也可以對(duì)未來幾期做出預(yù)測，為預(yù)測滑坡的變形發(fā)展提供依據(jù)。

3)如果觀測不平穩(wěn)，一些觀測值有較大起伏或觀測值含有較大誤差，濾波偏差會(huì)較大；初值選取與實(shí)際觀測值有一定偏差時(shí)，預(yù)報(bào)結(jié)果前面較少幾期偏差會(huì)較大，所以應(yīng)用Kalman濾波也是有一定的局限性。

[1]姚潁康，周傳波，郭廖武，等.深凹露天礦山巖質(zhì)邊坡穩(wěn)定性預(yù)測[J].金屬礦山，2008(3):42-45.

[2]張正祿，黃全義，文鴻雁，等.工程的變形監(jiān)測分析與預(yù)報(bào)[M].北京：測繪出版社，2007.

[3]高雅萍，馮曉亮.GPS變形監(jiān)測動(dòng)態(tài)數(shù)據(jù)處理中卡爾曼濾波的應(yīng)用[J].海洋測繪，2006,26(4)：36-38.

[4]許國輝，余春林.卡爾曼濾波模型的建立及其在施工變形測量中的應(yīng)用[J].測繪通報(bào)，2004(4):22-23.

[5]陳帥，趙麗，李超，等.卡爾曼濾波在壩體變形監(jiān)測中的應(yīng)用[J].全球定位系統(tǒng)，2011(5):83-84.

[6]陸付民，王尚慶，李勁.離散卡爾曼濾波法在滑坡變形預(yù)測中的應(yīng)用[J].水利水電科技進(jìn)展，2009,29(4):6-9.

[7]劉勝，張紅梅.最優(yōu)估計(jì)理論[M].北京：科學(xué)出版社，2011：129-164.

[8]張樹京，齊立心.時(shí)間序列簡明教程[M].北京：清華大學(xué)出版社，北方交通大學(xué)出版社，2003：49-54.

[9]修延霞，候凱.卡爾曼濾波在大壩變形監(jiān)測中的應(yīng)用[J].城市勘測，2010(1):92-94.

[10]溫德娟，滕壽仁，楊子榮.白音華3號(hào)露天礦礦坑充水因素研究[J].煤炭技術(shù)，2006,25(7):97-98.

[11]劉衛(wèi)國.Matlab程序設(shè)計(jì)教程[M].2版.北京：中國水利水電出版社，2010.

[責(zé)任編輯：劉文霞]

Application of Kalman filtering to the landslide monitoring data procession based on auto-regression method

XU Ai-gong1， ZHANG Jian-long1， SUI Xin1， LIU Wei2， ZHANG Zhao-nan3

It introduces a basic model of Kalman filtering.The kinematic method establishes the observation equation and state equation based on location parameters and rate as the state vector, acceleration as the noise vector.The auto-regression method is introduced into the establishment of the state equation and observation equation, simply compared with the kinematic method.Taking one opencast mine landslide dynamic monitoring data as an example, it analyzes the model application.The results show that the deviation of filtering value and observed value is between 1 to 3cm, and the filtering image is smoother than the original observation image.It shows that the model is appropriate, and it can reflect the dynamic change process of landslide.So the method provides a powerful mathematical tool for this opencast mine landslide monitoring and prediction.

auto-regression method; Kalman filtering; landslide; forecast

2013-08-30

精密工程與工業(yè)測量國家測繪地理信息局重點(diǎn)實(shí)驗(yàn)室開放基金資助項(xiàng)目(PF2012-8)；現(xiàn)代城市測繪國家測繪地理信息局重點(diǎn)實(shí)驗(yàn)室開放課題資助項(xiàng)目(20111203W)

徐愛功(1963-)，男，教授，博士生導(dǎo)師.

P642

：A

：1006-7949(2014)09-0052-04

(1.School of Geomatics, Liaoning Technical University, Fuxin 123000, China; 2.Luoyang Electronic Equipment Test Center, Luoyang 471000, China; 3.Physical Exploration and Surveying Team, Northeast Coalfield Geology Bureau, Shenyang 110101, China)