基于自適應卡爾曼濾波的地表移動變形預報研究

陳長坤,焦寶文,喬方,石長偉,肖明,蘇迪

(1.安徽理工大學 測繪學院,安徽 淮南 232001； 2.廣東省重工建筑設計院有限公司,廣東 廣州 510700)

0 引 言

卡爾曼濾波通過建立狀態方程和量測方程描述系統的動態過程,依據濾波增益矩陣的變化,從量測數據中定量識別和提取有效信息,修正狀態參量,無須存儲各個不同時刻的觀測數據,便于實時數據處理以及預報[1-2]。因此,卡爾曼濾波廣泛應用于各種工程變形、滑坡監測、大壩監測等動態測量系統中[3-8]。卡爾曼濾波的應用需要動態系統的數學模型和噪聲的先驗知識,但在許多條件下,它們是未知的或近似已知的,即需研究動態系統的數學模型和系統中相關噪聲的統計學特性。動態系統的數字模型和相關噪聲統計特性兩者的先驗信息一般總是存在各種誤差,自適應卡爾曼濾波正是為了克服這一缺點而提出的。它可以解決濾波的發散問題,從而減弱模型誤差的影響,使濾波結果更接近于真實值[9-10],從而更好地預報各種工程變形例如地表移動變形監測[11]。

本文利用GNSS CORS地表移動自動化監測系統[12]對地表移動變形進行實時監測和分析。GNSS CORS在變形監測中具有實時、高效、高精度、自動化程度高的優點,現在已經廣泛應用在地鐵、大壩、橋梁以及滑坡監測、礦山開采沉陷等不同領域。對GNSS CORS連續運行實時監測站監測點的位置坐標序列進行實時濾波預報，并對比標準卡爾曼濾波預報值、自適應卡爾曼濾波預報值、實測值，對GNSS CORS地表移動實時監測進行實時預報分析。

1 標準卡爾曼濾波模型

1.1 卡爾曼濾波狀態方程

設某一GNSS CORS連續運行監測網由n個監測站構成,可測量得到監測站監測點的3維位置坐標序列(本文直接得到3維位置坐標序列是適用于礦區的BJ-54坐標系下的高斯平面直角坐標和正常高),將t時刻(歷元)的3維位置坐標及其速率構成狀態向量。GNSS CORS連續運行監測站i在時刻(歷元)t的位置坐標為ξi(t),瞬時速率為λi(t),瞬時加速率為Ωi(t),可將瞬時加速率看作一種隨機干擾,則ξi(t)、λi(t)、Ωi(t)有以下微分關系:

(1)

記i點的狀態向量為

=Xi(t)Yi(t)Zi(t)

則該監測點的狀態方程為

(2)

把n個監測點的狀態方程組合可得全網的狀態方程為

Xk+1=Φk+1,kXk+Γk+1,kΩk.

(3)

1.2 卡爾曼濾波觀測方程

測量獲得第i個GNSS CORS監測站監測點的3維位置序列Li,k+1為觀測值,某一監測站監測點在第k+1歷元的觀測方程為

Li,k+1=ξi,k+1+Δti,k+1λi,k+1+Δi,k+1,

(4)

式中:Δti,k+1=ti,k+1-tk+1,ti,k+1為Lij的觀測時刻;tk+1為本次k+1歷元觀測的中心時刻,在地表移動變形觀測的過程中,短時間內(數個歷元,在本文中可以忽略不計)地表的移動變形速率可以忽略不計,可得全網的觀測方程為

Lk+1=Bk+1Xk+1+Δk+1.

(5)

狀態方程和觀測方程兩者共同組成GNSS CORS連續運行實時監測網的卡爾曼濾波模型:

Xk=Φk,k-1Xk-1+Γk,k-1Ωk-1,

Lk=BkXk+Δk,

(6)

式中:Φk,k-1為k-1到k歷元的系統一步轉移矩陣;Γk,k-1為系統噪聲矩陣;Ωk-1為k-1歷元的系統噪聲;Bk為k歷元系統的觀測矩陣;Δk為k歷元系統的觀測噪聲;Xk為k歷元系統待估狀態參數;Lk為k歷元系統觀測向量矩陣。Xk,Lk均為GNSS監測站監測點的3維位置和速度向量。

1.3 卡爾曼濾波算法

卡爾曼濾波的函數模型由式(6)給出,其隨機模型為

E(Ωk)=0,

E(Δk)=0,

Cov(Ωk,Ωj)=DΩ(k)δkj,

Cov(Δj,Δj)=DΔ(k)δkj,

其中:DΩ(k)為系統的動態噪聲方差陣,是非負定方差陣;DΔ(k)為系統的觀測噪聲方差陣,Δ為正定方差陣;δkj為Kronecker函數:

動態噪聲與觀測噪聲完全不相關,即:

Cov(Ωk,Δj)=0

用標準卡爾曼濾波方程處理GNSS監測站各歷元位置序列,其濾波方程為

X(k/k)=X(k/(k-1))+Jk

[Lk-BkX(k/(k-1))] ,

(7)

DX(k/k)=(I-JkBk)DX(k/(k-1)) .

(8)

X(k/k)為狀態濾波方程;DX(k/k)濾波誤差協方差陣,其中:

X(k/(k-1))=Φk,k-1X((k-1)/(k-1))

(9)

(10)

(11)

式中:X(k/(k-1))為一步預測值;Jk為濾波增益矩陣;Lk-BkX(k/(k-1))為預測殘差;X(k/(k-1))為預測方差陣;DX(k/(k-1))為預報誤差協方差陣。

式(7)和式(8)即是GNSS CORS連續運行監測站監測點的三維位置序列的標準卡爾曼濾波的遞推計算公式,在確定初始狀態后即可求得監測點三維位置的濾波值,X(k/k)即預報值。

2 方差補償自適應卡爾曼濾波模型

自適應卡爾曼濾波就是在對觀測數據進行遞推濾波的同時,不斷地對未知的或不確定的模型參數和噪聲的統計特性進行適當的估計和修正,以減小模型誤差,使濾波結果更接近于真實值[11]。方差補償自適應卡爾曼濾波是利用預測殘差對動態噪聲的協方差向量進行修正,計算出更接近實際的狀態向量。它的公式推導過程如下:

由第1章中標準卡爾曼濾波方程式(7)～(11)

假定{Ωk}和{Δk}為正態序列,X0為正態向量,{Ωk}和{Δk}為互不相關的零均值白噪聲序列。

定義i步的預測殘差為

Vk+i=Lk+i-L(k+i)/k,

(12)

式中:Lk+i、L(k+i)/k為第k+i歷元的觀測值和它的最佳預測值;Vk+i為預測殘差。

Lk+i=Bk+iΦ(k+i)/kXk+Δk+i,

則Vi+1的方差陣Dvv為

(13)

記

Bk+iΦ(k+i)/kΓγ/(γ-1)=A(k+i,γ)=[ahj(k+i,γ)],

式中,r=1,…,N;k=1,…,n,上標k+i,r表示與k+i,r有關。假定DΩr-1Ωr-1在觀測時間段tk+1,tk+2,…,tk+N上為常值對角陣,即

(14)

并記

有

記

(15)

其中ηk+i為零均值隨機變量,r=1,…,N.

令

(16)

又記

E=[Ek+1,…,Ek+N]T,

η=[ηk+1,…,ηk+N]T,

(17)

則有

E=AdiagDΩΩ+η.

(18)

式(18)為關于diagDΩΩ的線性方程組。當N≥r時,有唯一解。記diagDΩΩ的LS估計為

diagDΩΩ=(ATA)-1ATE.

(19)

根據上述各式求得任意時間段長度的DΩΩ,并作為動態噪聲協方差陣的實時估計。

由此得出的協方差可以對模型進行實時改正。自適應卡爾曼濾波在進行濾波的同時,能夠實時地按照相應的自適應方法對模型進行修改,有效地降低濾波過程中出現的發散現象。

3 算例分析

3.1 試驗簡介

為了驗證本文提出的自適應卡爾曼濾波預報的正確性,以淮南礦區朱集東礦1222(1)工作面上的GNSS CORS自動化監測系統進行試驗。該系統由1個基準站(CZJDM)和2個連續運行實時監測站(CORS1和CORS2)組成。

本文選用GNSS CORS1、CORS2連續運行監測站在2017年6月的監測數據進行橫向分析,選取兩個實時監測站每天1:00～2:00時段各歷元的標準卡爾曼濾波預報值和自適應卡爾曼濾波預報值(共60組數據,每組數據3600個歷元)進行濾波預報精度評價。

選取GNSS CORS1連續運行監測站2017年6月18日凌晨1:00～2:00的數據(觀測數據采樣率為1 s,共3600個歷元)進行縱向分析,分別采用標準卡爾曼濾波算法和自適應卡爾曼濾波算法對監測點的3維位置序列進行濾波預報,并獲取該實時監測站的實測3維位置序列,然后分別對比分析標準卡爾曼濾波、自適應卡爾曼濾波預報值與實測值的偏差。

選取GNSS CORS1、CORS2連續運行監測站2017年6月每周第一天1:00～5:00的監測數據進行內符合精度分析(共8組數據,即選取8個時間段,每組數據14400個歷元),分別采用標準卡爾曼濾波算法和自適應卡爾曼濾波算法對監測點的3維位置序列進行濾波預報,分別求出相應濾波值的內符合精度,并對比分析標準卡爾曼濾波、自適應卡爾曼濾波預報值的內符合精度。

3.2 數值結果與分析

橫向數據分析中,對各個預報時段所有觀測歷元的平面坐標中誤差和高程中誤差進行精度評定;通過該時段中每個歷元3維位置序列的預報偏差m求得該預報時段i整體的預報中誤差。

(20)

式中:mP為濾波預報的平面坐標中誤差；mH為濾波預報的高程中誤差；k為該預報時段的歷元個數；mx、my、mH分別為該預報時段歷元k的預報偏差。

表1 濾波預報平面精度信息統計

表2 濾波預報高程精度信息統計

從表1和表2中可以看出:對60個時段預報的精度進行對比分析,自適應卡爾曼濾波預報相較于標準卡爾曼濾波,其平面位置精度和高程位置精度得到大幅度提升(大約提高了2倍左右),同時精度分布得更為均勻,即自適應卡爾曼濾波預報后的成果更為穩定可靠。

縱向數據分析中,表3～表5以及圖1～圖3給出了本次試驗中標準卡爾曼濾波和自適應卡爾曼濾波在GNSS CORS連續運行監測站一個預報時段的平面X、Y坐標序列偏差和高程序列偏差信息。

表3 X坐標濾波預報偏差統計信息

從表3可以看出，對于平面位置X坐標序列，標準卡爾曼濾波預報最大偏差值3.21 mm，平均偏差值0.53 mm，偏差值小于0.5 mm的比率占57.67%；自適應卡爾曼濾波預報最大偏差值為0.80 mm，平均偏差值0.18 mm，預報偏差值相較于標準卡爾曼濾波平均縮減四分之一，偏差值小于0.5 mm的比率為99.94%.

表4 Y坐標濾波預報偏差統計信息

從表4可以看出，對于平面位置Y坐標序列。標準卡爾曼濾波預報最大偏差值為1.83 mm，平均偏差值0.33 mm，偏差值小于0.4 mm的比率占67.81%；自適應卡爾曼濾波預報最大值為0.80 mm，平均偏差值0.18 mm，預報偏差值相比較于標準卡爾曼濾波平均縮減四分之一，偏差值小于0.4 mm的比率為99.97%.

表5 高程濾波預報偏差統計信息

從表5看出，對于高程序列。標準卡爾曼濾波預報最大偏差值為8.27 mm，平均偏差值1.45 mm，偏差值小于1 mm的比率占43.39%；自適應卡爾曼濾波預報最大偏差值為1.29 mm，平均偏差值0.58 mm，預報偏差值相比較于標準卡爾曼濾波平均縮減四分之一，偏差值小于1 mm的比率為99.94%.由圖3可知高程序列偏差值分布更為均勻，用自適應卡爾曼濾波預報的成果穩定性更高。

3.3 內符合精度分析

對于某一時段,若觀測了k個歷元,取各歷元觀測值的平均值,即可求得該監測站的平均位置,并根據各觀測值與平均值之差,對該監測站的內符合精度進行評定,即

(21)

(22)

從表6可以看出2個GNSS CORS 連續運行監測站(共8個時段)的標準卡爾曼濾波預報值的平面內符合精度在4.3 ～6.2 mm之間,平均值5.4 mm,高程內符合精度在7.0 ～9.1 mm之間平均值為8.3 mm;從內符合精度可知,GNSS CORS連續運行監測站的平面位置精度滿足開采沉陷精度要求(根據《煤礦測量規程》,為了保證解算的開采沉陷的關鍵參數的精度,要求相鄰兩期間平面點位中誤差≤±20 mm,最弱的高程中誤差≤±10 mm,這就要求一次測量平面點位相對中誤差≤±14 mm,最弱點高程中誤差≤±7 mm),但是高程內符合精度比平面位置測量精度低一倍,標準卡爾曼濾波預報值大于7 mm,不滿足開采沉陷的精度要求。而自適應卡爾曼濾波預報值的平面內符合精度在3.2 ～4.3 mm之間,平均值為3.9 mm,高程內符合精度在5.9 ～7.4 mm之間,平均值6.6 mm;從內符合精度可知,GNSS CORS連續運行監測站的平面位置精度完全滿足開采沉陷精度要求,高程方向均基本滿足開采沉陷精度要求。且自適應卡爾曼濾波預報值相比較于標準卡爾曼濾波預報值,平面位置精度最低提高14.0%,最高提高37.9%,平均提高26.6%,高程方向精度最低提高10.0%,最高提高33%,平均提高20.1%.

從以上分析中可以看出在GNSS CORS連續運行監測站的位置坐標序列的預報中,其自適應卡爾曼濾波預報值相較標準卡爾曼濾波預報精度得到大幅提高,自適應卡爾曼濾波更能真實體現地表移動變形動態數據的實時性;對于標準卡爾曼濾波預報偏離實測值大,并且分布不均勻,自適應卡爾曼濾波預報偏差大幅度降低,預報偏差分布更為均勻,內符合精度也得到提高,使濾波過后的每個歷元的測量成果更為穩定和可靠,這為更進一步獲得穩定的移動變形量提供了堅實基礎。

表6 濾波內符合精度

4 結束語

本文中的自適應卡爾曼濾波算法,對GNSS CORS連續運行實時監測站各測點的3維位置序列進行預報,根據該研究區域的實測數據表明,相比于標準卡爾曼濾波算法,其預報值相較于實測值的偏差值縮小約四分之一,預報精度提高了約2倍,同時分布更為均勻,表明本文采用的自適應卡爾曼濾波對GNSS變形監測具有很好的剔除噪聲作用;平面位置內符合精度和高程內符合精度均提高了20%左右,自適應卡爾曼濾波預報的成果更為穩定和可靠,為獲得準確的實時變形量提供基礎。

但是本文中的自適應卡爾曼算法存在一定問題,還有待進一步改進。由圖3可知雖然自適應卡爾曼濾波預報偏差值相比與標準卡爾曼濾波偏差值明顯減小,但是自適應卡爾曼濾波預報偏差值最終趨近于某一值,而不是隨著歷元增加預報偏差值變更小。在某些歷元標準卡爾曼濾波的預報偏差值會小于自適應卡爾曼濾波,這正是標準卡爾曼濾波預報的隨機性導致的。本文的自適應卡爾曼濾波算法對于GNSS CORS連續運行實時監測站監測的地表移動變形監測的精度是足夠了。

另外,由于該自適應卡爾曼濾波算法通過建立狀態方程和量測方程來描述系統的動態過程,依據濾波增益矩陣的變化,從測量數據中定量識別和提取有效信息,修正狀態參量,以此對實時監測數據進行處理及預報。但由于GNSS CORS連續運行實時監測站本身存在系統誤差,或者是其他因素導致的粗差,該模型并不能探測并剔除粗差,這也是后續需要考慮的問題。

[1] 崔希璋,於宗儔,陶本藻,等.廣義測量平差[M].2版.武漢:武漢大學出版社,2009.

[2] 尹暉.時空變形分析與預報的理論和方法[M].北京:測繪出版社，2002.

[3] 余學祥,呂偉才.GPS監測網動態數據處理抗差Kalman濾波模型[J].中國礦業大學學報,2000,29(6):553-557.

[4] 余學祥,呂偉才.抗差卡爾曼濾波模型及其在GPS監測網中的應用[J].測繪學報,2001,30(1);27-31

[5] LU W C,XU S Q.Reasearch on Kahnan filtering algorithm for deformation informantion series of similar single-difference model[J].Journal of China University of Mining and Technology,2004,14(2):189-194.

[6] 余學祥,呂偉才.GPS監測網動態數據處理抗差Kalrnan濾波模型[J].中國礦業大學學報,2000,29(6):553-557.

[7] 楊元喜.動態Kalman濾波模型誤差的影響[J].測繪科學,2006,31(1):17-18

[8] 尹航,朱紀洪,周池軍,等.基于Kalman預報觀測器的增量動態逆控制[J].清華大學學報(自然科學版),2011.

[9] 賈志軍,單甘霖,程興亞,等.GPS動態定位中的自適應擴展卡爾曼濾波算法[J].軍械工程學院學報,2001(2):39-43.

[10] 高雅萍,張勤.方差補償自適應卡爾曼濾波在GPS滑坡監測中的應用研究[J].水土保持研究,2008(12): 150-152.

[11] 張福榮.自適應卡爾曼濾波在變形監測數據處理中的應用研究[D].西安:長安大學，2009.

[12] 余學祥.煤礦開采沉陷自動化監測系統[M].北京:測繪出版社，2014.