基于卡爾曼濾波的GNSS自動化監(jiān)測數(shù)據(jù)粗差分析

魏世玉，李 川

(1.重慶市地質(zhì)災(zāi)害自動化監(jiān)測工程技術(shù)研究中心(重慶地質(zhì)礦產(chǎn)研究院)，重慶 400042；2.重慶華地工程勘察設(shè)計院，重慶 400042)

基于卡爾曼濾波的GNSS自動化監(jiān)測數(shù)據(jù)粗差分析

魏世玉1,2，李 川1

(1.重慶市地質(zhì)災(zāi)害自動化監(jiān)測工程技術(shù)研究中心(重慶地質(zhì)礦產(chǎn)研究院)，重慶 400042；2.重慶華地工程勘察設(shè)計院，重慶 400042)

粗差探測一直是自動化變形監(jiān)測數(shù)據(jù)處理中的重點和難點，準(zhǔn)確有效的粗差探測是確保監(jiān)測工作能否反映變形體真實變形的關(guān)鍵。以GNSS自動化平面位移監(jiān)測數(shù)據(jù)為研究對象，利用卡爾曼濾波，建立監(jiān)測點坐標(biāo)序列粗差探測模型，對GNSS自動化監(jiān)測坐標(biāo)序列粗差探測過程進行了分析。通過工程實際監(jiān)測數(shù)據(jù)分析，驗證了模型的準(zhǔn)確性和有效性。

GNSS；自動化監(jiān)測；卡爾曼濾波；粗差探測

0 引言

在GNSS(Global Navigation Satellite System)自動化變形監(jiān)測中，由于受觀測條件的影響，如周跳、多路徑效應(yīng)以及接收機信號故障等，致使觀測數(shù)據(jù)中難免包含粗差，極大降低了觀測成果的準(zhǔn)確性和可靠性。粗差的存在給后續(xù)變形分析和解釋帶來了困難，甚至得出錯誤結(jié)論。因此，粗差探測工作是自動化變形監(jiān)測數(shù)據(jù)處理中的重要環(huán)節(jié)[1-4]。

對于粗差探測，目前普遍采用以均值漂移模型為基礎(chǔ)的粗差探測和以方差膨脹模型為基礎(chǔ)的抗差估計等方法[5-6]。這些方法具備嚴(yán)密的理論基礎(chǔ)，當(dāng)監(jiān)測網(wǎng)圖形強度較好且具備一定數(shù)量的多余觀測時，上述方法均能很好地消除粗差對觀測成果的影響。但在GNSS自動化變形監(jiān)測中，為了實時掌握變形體變形情況，通常需要提高監(jiān)測頻率，這就導(dǎo)致一次定位解算的觀測時間變短，多余觀測數(shù)也隨之減少。這時按上述方法很難對粗差進行準(zhǔn)確的定位和剔除。

卡爾曼濾波是20世紀(jì)60年代初由卡爾曼等人提出的一種遞推式濾波算法，其最大特點是能夠剔除系統(tǒng)中的隨機干擾，從而獲得逼近真實情況的有用信息。如果將粗差看作GNSS自動化監(jiān)測中的隨機干擾，便可利用卡爾曼濾波對其進行剔除，得到一組“干凈”的觀測值序列。因此，本文以GNSS自動化平面位移監(jiān)測數(shù)據(jù)為研究對象，利用卡爾曼濾波，建立監(jiān)測點坐標(biāo)序列粗差探測模型，對GNSS自動化監(jiān)測坐標(biāo)序列粗差探測過程進行詳細分析。最后通過實際監(jiān)測數(shù)據(jù)分析，驗證該模型的準(zhǔn)確性和有效性。

1 監(jiān)測點坐標(biāo)序列卡爾曼濾波模型

(1)

式中：Δt——相鄰時刻時間間隔。

由此可得系統(tǒng)狀態(tài)方程[7-11]：

(2)

記tk時刻有觀測向量Lk=Xk，則有觀測方程：

(3)

式中：Bk=[1 0]，為觀測矩陣；Δk——系統(tǒng)在tk時刻的觀測噪聲。

假設(shè)系統(tǒng)動態(tài)噪聲和觀測噪聲為互不相關(guān)的高斯白噪聲序列，即：

其中，D(Ωk)和D(Δk)分別為系統(tǒng)動態(tài)噪聲和觀測噪聲方差陣，δkj為Kronecker函數(shù)。

由廣義最小二乘原理可推導(dǎo)監(jiān)測點坐標(biāo)序列卡爾曼濾波遞推公式[8,11]：

(1)狀態(tài)向量一步預(yù)報值及其方差矩陣

(4)

(5)

(2)狀態(tài)向量濾波值及其方差矩陣

(6)

(7)

式中：Jk——濾波增益矩陣，具體形式為：

2 濾波模型初始值的確定

濾波啟動前，初始狀態(tài)往往難以精確確定，只能對其進行估算。但若與實際偏差較大，可能導(dǎo)致錯誤的濾波估計甚至引起發(fā)散。因此，合理確定初始值十分重要。

為近似求得初始狀態(tài)向量，分別將監(jiān)測點X、Y方向坐標(biāo)序列中前N個無粗差的數(shù)據(jù)進行二次多項式擬合:

(8)

同時得擬合殘差中誤差：

擬合系數(shù)協(xié)方差矩陣：

則由協(xié)方差傳播律得：

(9)

2.2 觀測噪聲方差陣D(Δk) 和系統(tǒng)噪聲方差陣D(Ωk)

對于觀測噪聲方差陣D(Δk)，主要由觀測值的數(shù)據(jù)處理方法來確定，載波相位相對定位中，觀測噪聲方差陣可由基線解算后的殘差向量進行估算。對于系統(tǒng)噪聲方差陣D(Ωk)，若將監(jiān)測點坐標(biāo)變化的加速度看作隨機擾動項，有[8]：

(10)

3 粗差探測

當(dāng)對系統(tǒng)狀態(tài)進行一步預(yù)報后，得一步預(yù)報殘差：

(11)

(12)

得預(yù)報殘差：

(13)

(14)

并從tk時刻起，繼續(xù)進行濾波。

圖1 卡爾曼濾波粗差探測流程Fig.1 Process of gross error detection using kalman filter

4 實驗分析

為驗證模型的有效性和準(zhǔn)確性，以重慶市三峽庫區(qū)某滑坡GNSS自動化位移監(jiān)測為分析對象，其監(jiān)測點的布設(shè)見圖2。

圖2 GNSS監(jiān)測點布設(shè)示意圖Fig.2 The layout diagram of GNSS monitoring point

選取同一剖面監(jiān)測點GNSS2(位于主滑面坡頂)和GNSS4(位于主滑面坡底)連續(xù)1個月(2015年5月)的監(jiān)測數(shù)據(jù)為分析對象，觀測時段長為2小時(即連續(xù)觀測2小時的數(shù)據(jù)作為一次靜態(tài)基線解算結(jié)果)，共359組有效坐標(biāo)序列數(shù)據(jù)。

為分析濾波模型對粗差探測的敏感性，選取坡頂觀測條件良好，且距基準(zhǔn)站較近的GNSS2在X軸和Y軸方向的坐標(biāo)序列(圖3)，經(jīng)分析發(fā)現(xiàn)，該序列無粗差情況。利用前5個坐標(biāo)分量序列分別進行二次擬合，以確定初始狀態(tài)，從序列第6個數(shù)據(jù)開始進行預(yù)報與濾波，得其在X軸和Y軸方向預(yù)報殘差序列見圖4。

對上述無粗差坐標(biāo)序列預(yù)報殘差統(tǒng)計分析得X方向預(yù)報殘差中誤差為0.65 mm，Y方向則為0.82 mm。由此可見，該濾波模型對粗差探測具有較高的靈敏性。

圖3 GNSS2坐標(biāo)序列(無粗差)Fig.3 Coordinates series of GNSS2(without gross error)

圖4 GNSS2預(yù)報殘差序列Fig. 4 Predicted residual series of GNSS2

為驗證模型探測粗差的有效性和準(zhǔn)確性，選取坡底觀測條件相對較差的GNSS4在X軸方向的坐標(biāo)序列進行分析(圖5)。該月監(jiān)測期間，GNSS4監(jiān)測點上接收機存在蓄電池?zé)o法供電的情況，導(dǎo)致某些觀測時段觀測值數(shù)量不足，另外該監(jiān)測點位于斜坡底部，衛(wèi)星信號質(zhì)量在某些時段可能會較差，因此解算出的坐標(biāo)序列中極有可能存在粗差。利用前5個無粗差坐標(biāo)序列進行二次擬合，以確定濾波初始狀態(tài)，并取3倍擬合殘差中誤差作為粗差探測閾值(μ=3.5 mm,μ′=4.5 mm),從序列第6個數(shù)據(jù)開始進行預(yù)報與濾波，得其在X軸方向預(yù)報殘差序列見圖6。

圖5 GNSS4坐標(biāo)序列(南北方向)Fig.5 Coordinates series of GNSS4

圖6 GNSS4預(yù)報殘差序列(南北方向)Fig.6 Predicted residual series of GNSS4

表1 坐標(biāo)序列粗差探測結(jié)果

由上述分析可知，該模型能夠有效探測出GNSS坐標(biāo)序列中不同大小的粗差。當(dāng)發(fā)現(xiàn)某一時刻觀測值存在粗差時，采用(14)給出的插值法進行插補后，繼續(xù)進行遞推和濾波，得GNSS4監(jiān)測點濾波后X軸方向坐標(biāo)序列見圖7。

圖7 GNSS4濾波后坐標(biāo)序列(南北方向)Fig.7 Coordinates series of GNSS4 after filtered

為檢驗濾波后坐標(biāo)序列中粗差是否完全剔除，最后對濾波后的坐標(biāo)序列再次進行遞推預(yù)報與濾波，得其預(yù)報殘差序列見圖8。

圖8 剔除粗差后預(yù)報殘差序列Fig.8 Predicted residual series after removing the gross error

上述殘差序列中，個別值在3～5 mm之間，這是由于監(jiān)測點的異動變形引起，而非粗差所致，由此證明了該模型的準(zhǔn)確性。

5 結(jié)語

本文將卡爾曼濾波應(yīng)用于GNSS自動化監(jiān)測粗差探測，通過理論分析和實際監(jiān)測數(shù)據(jù)驗證，得出以下幾點結(jié)論：

(1)卡爾曼濾波是一種對動態(tài)系統(tǒng)進行實時狀態(tài)估計的有效方法，系統(tǒng)數(shù)學(xué)模型和噪聲統(tǒng)計模型的合理建立是系統(tǒng)準(zhǔn)確狀態(tài)估計的關(guān)鍵。

(3)變形監(jiān)測數(shù)據(jù)粗差探測的關(guān)鍵在于如何區(qū)分粗差與變形。本文利用粗差僅在某一時刻具有突變特征而變形往往伴隨后續(xù)數(shù)據(jù)整體發(fā)生偏移，有效地剔除了GNSS坐標(biāo)序列中粗差。

[1] 宋迎春,朱建軍,陳正陽,等.Kalman濾波中測量粗差的探測與修復(fù)[J]. 工程數(shù)學(xué)學(xué)報,2008,25(3):437-441.SONGYingchun，ZHUJianjun，CHENZhengyang，etal.Thedetectionandcorrectionofgrosserrorinobservedvaluesforthekalmanfilter[J].ChineseJournalofEngineeringMathematics，2008,25(3):437-441.

[2] 劉超云,尹小波,張彬. 基于Kalman濾波數(shù)據(jù)融合技術(shù)的滑坡變形分析與預(yù)測[J].中國地質(zhì)災(zāi)害與防治學(xué)報,2015,26(4):30-35.LIUChaoyun，YINXiaobo，ZHANGBin.AnalysisandpredictionoflandslidedeformationsbasedondatafusiontechnologyofKalman-filter[J].TheChineseJournalofGeologicalHazardandControl,2015,26(4):30-35.

[3] 劉偉平,郝金明,汪平,等.Kalman濾在周跳探測與修復(fù)中的應(yīng)用[J]. 大地測量與地球動力學(xué),2009(6):101-103.LIUWeiping,HAOJinming,WANGPing,etal.StudyoncycleslipdetectionandrepairbyuseofKalmanFiltering[J].JournalofGeodesyandGeodynamics, 2009(6):101-103

[5] 朱建軍,曾卓喬. 污染誤差模型下的測量數(shù)據(jù)處理理論[J]. 測繪學(xué)報,1999(3):215-220.ZHUJianjun,ZENGZhuoqiao.TheTheoryofsurveyingadjustmentundercontaminatederrormodel[J].ActaGeodaeticaetCartographicaSinica，1991(3):215-220.

[6] 宋力杰,楊元喜. 均值漂移模型粗差探測法與LEGE法的比較[J]. 測繪學(xué)報,1999(4):295-300.SONGLijie，YANGYuanxi.ComparisonbetweendatasnoopingandLEGE[J].ActaGeodaeticaetCartographicaSinica,1999(4):295-300.

[7] 黑君淼.基于Kalman濾波的滑坡監(jiān)測數(shù)據(jù)處理與分析[D].安徽理工大學(xué)，碩士學(xué)位論文，2011.HEIJunmiao.DataprocessingandanalysisformonitoringlandslidedeformationbasedonKalmanfiltering[D].AnhuiUniversityofScienceandTechnology，Master’sDissertation，2011.

[8] 黃聲享,尹暈暉,蔣征,等. 變形監(jiān)測數(shù)據(jù)處理[M]. 武漢：武漢大學(xué)出版社, 2010.HUANGShengxiang,YINYunhui,JIANGZheng,etal.Deformationmonitoringdataprocessing[M].Wuhan:WuhanUniversityPress,2010.

[9] 張月超,陳義. 基于改進Kalman濾波算法的粗差修正及應(yīng)用[J]. 測繪與空間地理信息,2013,36(12):257-259. ZHANG Yuechao，CHEN Yi.Gross error revision of filter algorithm and its application in improved kalman[J].Geomatics & Spatial Information Technology，2013,36(12):257-259.

[10] 許國輝,張新長. 卡爾曼濾波模型粗差的探測及其在施工變形測量中的應(yīng)用[J]. 中山大學(xué)學(xué)報(自然科學(xué)版),2003,42(3):89-91. XU Guohui，ZHANG Xinchang.Gross error detection of kalman filtering model and its application in construction deformation measurements[J].Journal of Sun Yat-sen University(JCR Science Edition)，2003,42(3):89-91.

[11] 邱衛(wèi)寧,陶本藻,姚宜斌,等. 測量數(shù)據(jù)處理理論與方法[M].武漢：武漢大學(xué)出版社, 2008. QIU Weining，TAO Benzao，YAO Yibin, et al. Theories and methods of measuring data processing [M].Wuhan: Wuhan University Press,2008.

The gross error detection for GNSS automatic monitoring data based on kalman filter

WEI Shiyu1,2，LI Chuan1

(1.ChongqingEngineeringResearchCenterofAutomaticMonitoringforGeologicalHazards(ChongqingInstituteofGeologyandMineralResources),Chongqing400042,China; 2.ChongqingHuadiEngineeringSurveyandDesignInstitute,Chongqing400042,China)

Gross error detection has been the key points and difficulties in automatic deformation monitoring data processing, the accurately and effectively gross error detection is the key to make sure that the monitoring data can reflect the real deformation. This paper selected the plane displacement monitoring data as the research object，established the gross error detection model for the coordinates series of monitoring point using kalman filter, and analyzed the progress of gross error detection. Finally, through the actual monitoring data analysis, it verified the accuracy and effectiveness of the model.

GNSS；automatic monitoring；Kalman filter；gross error detection

10.16031/j.cnki.issn.1003-8035.2017.01.23

2016-03-30;

2016-05-23

基于北斗定位的地質(zhì)災(zāi)害監(jiān)測設(shè)備研發(fā)與示范(CQGT-KJ-2014045)；地質(zhì)災(zāi)害多維度空間形變監(jiān)測儀研發(fā)與示范(CQGT-KJ-2014044)

魏世玉(1987-)，男，云南保山人，碩士，工程師，主要研究方向為地質(zhì)災(zāi)害智能監(jiān)測。E-mail:1832382150@qq.com

P694

A

1003-8035(2017)01-0146-05