開采沉陷監(jiān)測數(shù)據(jù)處理的分段卡爾曼濾波模型

郝登程，王國瑞，李培現(xiàn)，沈佳琦，曹禹錫，楊中輝

中國礦業(yè)大學(xué)(北京)地球科學(xué)與測繪工程學(xué)院，北京 100083

有效監(jiān)測礦區(qū)開采沉陷，能夠?qū)ΦV區(qū)經(jīng)濟(jì)發(fā)展、建筑物安全使用以及災(zāi)害預(yù)警預(yù)報起到指導(dǎo)作用[1]。隨著空間對地觀測技術(shù)的發(fā)展，具有高精度、高頻率、全天時等特點(diǎn)的GNSS技術(shù)逐漸取代傳統(tǒng)測量技術(shù)，在煤礦區(qū)沉陷監(jiān)測中得到了廣泛的應(yīng)用[2]，得到了大量的礦區(qū)沉降信息。但數(shù)據(jù)受多路徑效應(yīng)[3]、儀器穩(wěn)定性、外界環(huán)境[4-5]等因素影響，數(shù)據(jù)波動明顯，在長周期監(jiān)測中常存在較大偏差，影響決策及預(yù)警預(yù)報的可靠性[6]。與城市地表沉降監(jiān)測相比，礦區(qū)沉陷速度快，變化大，數(shù)據(jù)處理更為困難[7]。

1960年以來，卡爾曼濾波在電子、控制、導(dǎo)航、無線傳感器網(wǎng)絡(luò)等多個領(lǐng)域得到了廣泛應(yīng)用[8-9]?？柭鼮V波算法因其計算結(jié)果準(zhǔn)確、可靠等特點(diǎn)在軍事雷達(dá)跟蹤[10-11]、無線定位[12-13]等數(shù)據(jù)處理領(lǐng)域，發(fā)揮了巨大作用。隨著抗差卡爾曼濾波[14-15]、自適應(yīng)卡爾曼濾波[16-17]以及抗差自適應(yīng)卡爾曼濾波[18-19]等新算法的應(yīng)用，卡爾曼濾波在諸多研究領(lǐng)域發(fā)揮了越來越重要的作用。但因礦區(qū)沉陷速度快、變形量大等特點(diǎn)的影響，濾波結(jié)果偏差較大，在這方面的應(yīng)用較少。為提高濾波效果，本文將礦區(qū)沉陷監(jiān)測數(shù)據(jù)分為三段[20]，對礦區(qū)沉陷較為緩和的部分，采用標(biāo)準(zhǔn)卡爾曼濾波模型對數(shù)據(jù)進(jìn)行處理，并設(shè)定閾值，剔除含有較大粗差的數(shù)據(jù)。對礦區(qū)沉陷較大的部分，采用加入修正數(shù)卡爾曼濾波的方式對數(shù)據(jù)進(jìn)行處理[21]。以寧夏某礦區(qū)2015年2月至2020年5月間每小時采樣率的監(jiān)測數(shù)據(jù)處理為例，采用分段結(jié)合修正數(shù)的卡爾曼濾波方法，改善了計算結(jié)果的可靠性。

1 卡爾曼濾波基本原理及參數(shù)設(shè)置

1.1 基本原理

卡爾曼濾波是一種采用線性系統(tǒng)狀態(tài)方程，通過系統(tǒng)輸入輸出觀測數(shù)據(jù)，對系統(tǒng)狀態(tài)進(jìn)行最優(yōu)估計的算法。

卡爾曼濾波器狀態(tài)方程：

X(k+1)=φX(k)+ΓW(K)

(1)

卡爾曼濾波器觀測方程：

Y(k)=HX(k)+V(k)

(2)

式中，φ為狀態(tài)轉(zhuǎn)移矩陣；Γ為噪聲驅(qū)動矩陣；W(K)為過程噪聲；Y(k)為k時刻的觀測信號；H為觀測矩陣；V(k)為觀測噪聲；W(K)和V(k)分別為均值為零方差各為Q和R的高斯白噪聲[22]。

卡爾曼濾波器包含預(yù)測、更新兩個主要過程。預(yù)測過程通過時間更新建立當(dāng)前狀態(tài)的先驗估計，即對狀態(tài)變量與誤差協(xié)方差進(jìn)行估計，為下一時刻提供先驗估計值。更新過程是在當(dāng)前觀測量及先驗估計的基礎(chǔ)上，得出當(dāng)前狀態(tài)的濾波估計值[23]。

卡爾曼濾波具體實現(xiàn)過程[24]如下：

(3)

(4)

(5)

K(k+1)=P(k+1|k)HT
[HP(k+1|k)HT+R]-1

(6)

(4) 一步預(yù)測協(xié)方差矩陣：協(xié)方差矩陣預(yù)測是利用k時刻的協(xié)方差矩陣P(k|k)通過狀態(tài)轉(zhuǎn)移矩陣φ、噪聲驅(qū)動矩陣Γ及過程噪聲W(K)中誤差Q，來預(yù)測k+1時刻協(xié)方差P(k+1|k)。

P(k+1|k)=φP(k|k)φT+ΓQΓT

(7)

(5) 協(xié)方差矩陣更新：根據(jù)預(yù)測協(xié)方差P(k+1|k)，計算更新協(xié)方差陣P(k+1|k+1)。

P(k+1|k+1)=[In-K(k+1)H]P(k+1|k)

(8)

(9)

卡爾曼濾波存在狀態(tài)預(yù)測式(3)與狀態(tài)更新式(4)兩個主要過程。式(6)～式(9)的增益矩陣、協(xié)方差陣預(yù)測與更新，為計算狀態(tài)更新值提供計算依據(jù)。利用增益矩陣與協(xié)方差矩陣不斷調(diào)整預(yù)測值與觀測值間的關(guān)系，提高濾波值準(zhǔn)確性?？柭鼮V波方程是一組遞推公式 ，可不斷循環(huán)進(jìn)行預(yù)測與更新。在求解狀態(tài)向量時，不需要保存大量的歷史觀測數(shù)據(jù)；當(dāng)獲得新的量測值時，可求得新的狀態(tài)向量濾波值，以達(dá)到濾波效果[25]?？柭鼮V波流程如圖1所示。

圖1 卡爾曼濾波流程

1.2 參數(shù)設(shè)置

一般假設(shè)觀測噪聲V(k)和過程噪聲W(K)為高斯白噪聲，方差分別是R和Q。

觀測噪聲V(k)取值一般與傳感器的測量精度有關(guān)，其方差R可通過對數(shù)據(jù)長期的觀測統(tǒng)計得出。R值影響濾波收斂的快慢，取值越小收斂越快，取值越大收斂越慢，進(jìn)而影響濾波效果。對沉陷監(jiān)測數(shù)據(jù)分段回歸分析，提取沉陷監(jiān)測數(shù)據(jù)開始階段與衰退階段中穩(wěn)定數(shù)據(jù)，計算中誤差的平均值作為R值。

過程噪聲W(K)與外界自然因素有關(guān)，其方差Q影響因素復(fù)雜，難以準(zhǔn)確獲取，通過不斷調(diào)整來確定對于某過程的最優(yōu)Q值[5]。

2 卡爾曼濾波程序開發(fā)

Python編程語言是一種應(yīng)用廣泛的高級編程語言，由Guido van Rossum在1989年創(chuàng)建。本研究程序開發(fā)過程中使用到了tkinter、matplotlib、numpy、xlwt、xlrd等Python庫[26]。

2.1 數(shù)據(jù)處理流程

采用回歸分析算法對沉陷監(jiān)測數(shù)據(jù)進(jìn)行分割，計算K值。將K值與0比較，規(guī)定K值連續(xù)小于0的地表移動階段為活躍階段，活躍階段開始前為開始階段，活躍階段結(jié)束之后為衰退階段。

對沉陷較為平緩階段的監(jiān)測數(shù)據(jù)采用標(biāo)準(zhǔn)卡爾曼濾波處理，通過設(shè)定閾值剔除粗大誤差。對沉陷較快階段的監(jiān)測數(shù)據(jù)做一次差分，使用卡爾曼濾波對差分?jǐn)?shù)據(jù)濾波處理，利用得到的濾波值作為修正數(shù)對沉陷監(jiān)測數(shù)據(jù)濾波處理。將分段后的數(shù)據(jù)濾波結(jié)果拼接，得到最終計算結(jié)果。數(shù)據(jù)處理流程如圖2所示。

圖2 沉陷監(jiān)測數(shù)據(jù)處理流程

2.2 參數(shù)R值計算流程

對于R值可通過觀測儀器的精度獲得，也可以通過提取監(jiān)測數(shù)據(jù)中穩(wěn)定監(jiān)測數(shù)據(jù)段的中誤差進(jìn)行估計。具體計算過程如下：

(1) 對數(shù)據(jù)進(jìn)行分割，分割后每份數(shù)據(jù)集中含有指定數(shù)量的數(shù)據(jù)。

(2) 對分割后的每份數(shù)據(jù)回歸分析。

(3) 通過對每份數(shù)據(jù)進(jìn)行F檢驗，挑選出符合K=0要求的數(shù)據(jù)。

(4) 對符合要求的數(shù)據(jù)集數(shù)據(jù)分別計算誤差及其平均值，該平均值設(shè)為R值。

3 工程應(yīng)用

3.1 研究區(qū)域概況

寧東煤炭基地是黃河中上游流域重要的煤炭生產(chǎn)基地。近年來，礦區(qū)受開采影響區(qū)域沉陷面積3 485 km2，土地沙化、水土流失嚴(yán)重。為掌握該地區(qū)開采沉陷情況，采用GNSS對該區(qū)域進(jìn)行監(jiān)測。某GNSS監(jiān)測點(diǎn)位于工作面正上方，自2015年2月到2020年5月對該點(diǎn)監(jiān)測，每隔1 h獲取1次數(shù)據(jù)，數(shù)據(jù)監(jiān)測量為43 556個。該監(jiān)測點(diǎn)地表沉陷速度快、下沉劇烈，沉陷監(jiān)測數(shù)據(jù)受外界影響大，其中含有粗差較大的偏離數(shù)據(jù)。原始數(shù)據(jù)的下沉過程曲線如圖3所示。

圖3 原始數(shù)據(jù)下沉過程曲線

3.2 數(shù)據(jù)分割

利用回歸分析的方法對沉陷監(jiān)測數(shù)據(jù)進(jìn)行數(shù)據(jù)分割。通過與0值比較發(fā)現(xiàn)，沉陷變化快的數(shù)據(jù)段位于0值界限以下，將監(jiān)測數(shù)據(jù)分割為三部分，分割點(diǎn)為4 521，6 761。K值計算結(jié)果如圖4所示。

圖4 K值計算結(jié)果

第一部分有監(jiān)測數(shù)據(jù)4 521個，該段累積沉降122 mm，平均每小時下沉量0.03 mm，監(jiān)測數(shù)據(jù)變化小，沉降較穩(wěn)定。數(shù)據(jù)受粗差影響嚴(yán)重且較為發(fā)散。分割后第一部分?jǐn)?shù)據(jù)下沉過程曲線如圖5所示。

圖5 分割后第一部分?jǐn)?shù)據(jù)下沉過程曲線

第二部分有監(jiān)測數(shù)據(jù)2 240個，該段累積沉降約2 335 mm，平均每小時下沉量為1.04 mm，監(jiān)測數(shù)據(jù)變化大，沉降速度快。分割后第二部分?jǐn)?shù)據(jù)下沉過程曲線如圖6所示。

圖6 分割后第二部分?jǐn)?shù)據(jù)下沉過程曲線

第三部分有監(jiān)測數(shù)據(jù)36 795個，該段數(shù)據(jù)累積沉降約565 mm，平均每小時下沉量為0.01 mm，監(jiān)測數(shù)據(jù)變化小，沉降較穩(wěn)定。數(shù)據(jù)存在粗差影響大且較為發(fā)散。分割后第三部分?jǐn)?shù)據(jù)下沉過程曲線如圖7所示。

圖7 分割后第三部分?jǐn)?shù)據(jù)下沉過程曲線

3.3 R值計算

將每500個數(shù)據(jù)為一組進(jìn)行劃分。通過F檢驗提取穩(wěn)定數(shù)據(jù)，計算結(jié)果如圖8所示，共有2 500個5組數(shù)據(jù)符合要求。求得5組數(shù)據(jù)中誤差平均值為0.386 mm，R值為0.386。

圖8 F檢驗圖

3.4 數(shù)據(jù)濾波處理

使用卡爾曼濾波模型對第一部分?jǐn)?shù)據(jù)進(jìn)行濾波處理，設(shè)定閾值來剔除粗差數(shù)據(jù)對濾波的影響。

計算卡爾曼濾波模型中預(yù)測值與實際觀測值之差，將其與閾值進(jìn)行比較，大于閾值的點(diǎn)認(rèn)為觀測值偏離較大，使用預(yù)測值代替實際觀測值，以此實現(xiàn)粗差數(shù)據(jù)剔除。

第一部分?jǐn)?shù)據(jù)濾波結(jié)果如圖9所示，濾波后能夠有效剔除粗差，降低數(shù)據(jù)的離散，濾波效果明顯。

圖9 第一部分?jǐn)?shù)據(jù)卡爾爾曼濾波圖

對第二部分沉陷監(jiān)測數(shù)據(jù)做一次差分，利用卡爾曼濾波模型對第二部分差分?jǐn)?shù)據(jù)進(jìn)行濾波處理，得到差分?jǐn)?shù)據(jù)濾波值?？柭鼮V波降低了差分?jǐn)?shù)據(jù)的離散度，結(jié)果如圖10所示。

圖10 第二部分?jǐn)?shù)據(jù)差分值卡爾曼濾波圖

在卡爾曼濾波模型預(yù)測值計算中，加入差分?jǐn)?shù)據(jù)濾波值作為修正數(shù)對第二部分?jǐn)?shù)據(jù)卡爾曼濾波處理，結(jié)果如圖11所示。

圖11 第二部分?jǐn)?shù)據(jù)卡爾爾曼濾波圖

觀察第二部分?jǐn)?shù)據(jù)濾波細(xì)節(jié)(圖12)可知，濾波值能夠很好地反映沉陷監(jiān)測變化曲線，且降低沉陷監(jiān)測數(shù)據(jù)的離散，濾波效果好。

圖12 第二部分?jǐn)?shù)據(jù)卡爾爾曼濾波圖細(xì)節(jié)

使用卡爾曼濾波模型對第三部分?jǐn)?shù)據(jù)進(jìn)行濾波處理，并設(shè)定閾值剔除粗差數(shù)據(jù)對濾波的影響，結(jié)果如圖13所示?？梢姡瑢Υ植畹奶蕹Ч^好，濾波效果較理想。

圖13 分割后第三部分?jǐn)?shù)據(jù)卡爾曼濾波圖

對三部分處理后的數(shù)據(jù)進(jìn)行拼接得到最終的數(shù)據(jù)處理結(jié)果，如圖14所示。

圖14 卡爾曼濾波數(shù)據(jù)處理結(jié)果

由處理結(jié)果知，通過分段卡爾曼濾波處理消除了較大粗差對數(shù)據(jù)的影響，降低了數(shù)據(jù)的離散，濾波后的結(jié)果與實測結(jié)果過程曲線吻合，濾波效果好。

3.5 結(jié)果分析

對于沉陷緩慢的監(jiān)測數(shù)據(jù)，通過設(shè)定閾值使用標(biāo)準(zhǔn)卡爾曼濾波模型處理，能有效剔除粗大誤差和降低數(shù)據(jù)發(fā)散。

對于下沉速度較快的數(shù)據(jù)，采用加入修正數(shù)(前一時刻下沉量濾波值)的卡爾曼濾波模型處理，解決了監(jiān)測下沉引起的濾波失真問題，有效改善了濾波結(jié)果。標(biāo)準(zhǔn)卡爾曼濾波計算結(jié)果與加入修正數(shù)卡爾曼濾波計算處理結(jié)果的對比，如圖15所示。

圖15 卡爾曼濾波數(shù)據(jù)處理結(jié)果對比

使用標(biāo)準(zhǔn)卡爾曼濾波處理沉陷變化快的數(shù)據(jù)時，由于曲線存在非線性下沉，沉陷數(shù)據(jù)在一步預(yù)測過程中的預(yù)測值始終位于實際觀測值的上方，致使卡爾曼濾波模型濾波值與實際沉陷觀測值之間存在向上偏離的問題，因此濾波結(jié)果不能有效反應(yīng)沉陷趨勢，濾波效果差。加入修正數(shù)后，濾波值能夠有效反映沉陷曲線變化趨勢，濾波值與實測值吻合好，結(jié)果更加準(zhǔn)確可靠。

4 結(jié) 論

(1) 針對長周期、高頻率的GNSS地表沉陷監(jiān)測數(shù)據(jù)存在沉降變化大的問題，采用Python編程語言、用分段方法進(jìn)行卡爾曼濾波處理。采用回歸分析將監(jiān)測數(shù)據(jù)分為三部分，針對沉陷監(jiān)測數(shù)據(jù)較為穩(wěn)定的下沉啟動階段和下沉衰退階段，使用標(biāo)準(zhǔn)卡爾曼濾波模型及設(shè)定閾值的方法處理，針對地表沉陷活躍段濾波結(jié)果向上偏離的問題，使用差分濾波后的值作為修正數(shù)對監(jiān)測數(shù)據(jù)進(jìn)行濾波處理，取得了較好的濾波效果。

(2) 以寧東某礦5年的GNSS監(jiān)測數(shù)據(jù)為研究對象，改善的卡爾曼濾波方法降低了噪聲和地表快速下沉趨勢的影響，提高了結(jié)果的精度和可靠性。