一種提高gnss實時監(jiān)測數(shù)據(jù)精度的算法
【專利摘要】本發(fā)明公開了一種提高GNSS實時監(jiān)測數(shù)據(jù)精度的算法，采用卡爾曼濾波器通過利用狀態(tài)方程和量測方程來估計出整個系統(tǒng)在每一個釆樣間隔的所有狀態(tài)，通過對各個狀態(tài)的統(tǒng)計與預(yù)測，反復(fù)更新，最終得到濾波結(jié)果。在計算機中其濾波方程的計算更新過程包含了以下兩個方面：量測方程的更新和時間方程的更新，其基本原理為用前一次的估計值遞推當(dāng)前時刻的狀態(tài)值，然后用當(dāng)前時刻的量測值修正狀態(tài)值，可以大大提高GNSS實時監(jiān)測數(shù)據(jù)的精度。
【專利說明】
一種提高GNSS實時監(jiān)測數(shù)據(jù)精度的算法
技術(shù)領(lǐng)域
[0001] 本發(fā)明涉及計算機數(shù)據(jù)處理算法領(lǐng)域，尤其涉及一種提高GNSS實時監(jiān)測數(shù)據(jù)精度 的算法。
【背景技術(shù)】
[0002] 卡爾曼(Kalman)濾波是由卡爾曼(Kalman R E)博士為滿足計算機用于人造地球 衛(wèi)星定軌和導(dǎo)航等計算的要求，在1960年提出的一種新的線性濾波模型，通常稱它為卡爾 曼濾波方法。卡爾曼濾波是一種具有無偏性的遞推線性最小方差估計，即估計誤差的均值 或數(shù)學(xué)期望為零。在計算方法上，卡爾曼濾波采用遞推形式，即在t-Ι時刻估值的基礎(chǔ)上，利 用t時刻的觀測值，遞推得倒t時刻的狀態(tài)估值。由于一次僅處理一個時刻的觀測值，無需存 儲先前的觀測數(shù)據(jù)，因而計算量大大減少。
[0003] 卡爾曼濾波能有效的消除噪聲干擾，獲得逼近真實情況的有用信息。這種濾波方 法很適合處理動態(tài)系統(tǒng)的觀測數(shù)據(jù)，因而廣泛應(yīng)用于動態(tài)測量。但是卡爾曼濾波是一個連 續(xù)系統(tǒng)，通過離散化后，才能讓計算機實現(xiàn)卡爾曼濾波的模擬，這里通過一種簡單方法的是 卡爾曼濾波離散化。然而這種濾波離散化并不能有效地提高GNSS實時監(jiān)測數(shù)據(jù)的精度。因 此，要對該濾波進行改進，以提高GNSS實時監(jiān)測數(shù)據(jù)的精度。

【發(fā)明內(nèi)容】

[0004] 本發(fā)明主要解決的技術(shù)問題是提供一種提高GNSS實時監(jiān)測數(shù)據(jù)精度的算法，采 用卡爾曼濾波器通過利用狀態(tài)方程和量測方程來估計出整個系統(tǒng)在每一個釆樣間隔的所 有狀態(tài)，通過對各個狀態(tài)的統(tǒng)計與預(yù)測，反復(fù)更新，最終得到濾波結(jié)果，在計算機中其濾波 方程的計算更新過程包含了以下兩個方面:量測方程的更新和時間方程的更新，其基本原 理為用前一次的估計值遞推當(dāng)前時刻的狀態(tài)值，然后用當(dāng)前時刻的量測值修正狀態(tài)值，可 以大大提高GNSS實時監(jiān)測數(shù)據(jù)的精度。
[0005] 為解決上述技術(shù)問題，本發(fā)明采用的一個技術(shù)方案是:提供了一種提高GNSS實時 監(jiān)測數(shù)據(jù)精度的算法，采用卡爾曼濾波器通過利用狀態(tài)方程和量測方程來計算出整個系統(tǒng) 在每一個釆樣間隔的所有狀態(tài)，通過對各個狀態(tài)的統(tǒng)計與預(yù)測，反復(fù)更新，最終得到濾波結(jié) 果，
[0006] 其中，連續(xù)系統(tǒng)的離散化后的狀態(tài)方程和量測方程分別為：
[0007] Xk=?k/k-iXk-i+r k-1ffk-1 (1-1)，
[0008] Lk = HkXk+Vk (1-2)，
[0009] 式（1 -1)和（1 -2)中，Xk和Xk-i分別為tk時刻和tk-i時刻的（η X 1)系統(tǒng)的狀態(tài)向量； Wk-1為tk-i時刻(nXl)動態(tài)噪聲向量；?k/k-i為時間tk-iStk的系統(tǒng)狀態(tài)轉(zhuǎn)移矩陣(nXn);H k 為tk時刻的(mXn)量測矩陣；刻系統(tǒng)的噪聲矩陣;Lk為tk時刻的(mXn)對系統(tǒng) 的觀測向量;V k為tk時刻的(m X 1)維量測噪聲向量；
[0010] 同時，所述狀態(tài)方程和量測方程中的UPV滿足：
[0012]式（1-3)中，正定矩陣Rk為量測噪聲方差序列的方差陣；非負定矩陣Qk為系統(tǒng)噪聲 序列的方差陣;是滿足如下條件：
[0014] 根據(jù)式（1-1)、（1-2)、（1-3)和（1-4)所描述的條件即可由Lk和XH求得xk，g卩k時刻 的最優(yōu)估計值。
[0015] 在本發(fā)明一個較佳實施例中，所述卡爾曼濾波器的卡爾曼濾波方程中狀態(tài)預(yù)報的 公式為：
[0016] = φ?/?-ι^Α-ι (1-5),
[0017] 式（1-5)中，是Χη的卡爾曼濾波估值，li/iM是通過Χη而算得的一步預(yù)測，即 通過k-Ι時刻以及之前時刻的量測值得到的。
[0018] 在本發(fā)明一個較佳實施例中，所述卡爾曼濾波器的卡爾曼濾波方程中狀態(tài)協(xié)方差 陣預(yù)報的公式為：
[0019] 乃-1 =Φ?'?-4 乃-1--必--t (.1-6)，.
[0020] 由式（1-6)，可知Pk/k-i是Pk-i加上噪聲方差得到的。
[0021 ]在本發(fā)明一個較佳實施例中，所述卡爾曼濾波器的卡爾曼濾波方程中狀態(tài)估計的 公式為：
[0022] 之-之 -丑it 之n) :(1-7),
[0023] 式(1-7)用來計算估值Xk，Xk是根據(jù)測量值Lk和一步預(yù)測計算出來的，式中，
[0024] ^HkXkl]i_^ = HkXk + V^-HkXkl]i^ = HkXk ik-\ + K (1-11)，
[0025] 其中，名^^是一步預(yù)測誤差，它的計算公式為：
[0026] =Xj. - Xk /k-l (1-12):
[0027] 為量測值Hk的一步預(yù)測，而私-/4之,^即為量測的一步預(yù)測誤差，由式 (1-11)可知，/4 -/4?·Μ4是由尤和Vk計算得到；由式（1-以)可知，矣是由尤i和知 得到，所以稱巧-為新息；
[0028] Kk稱為濾波增益矩陣，選取Kk的準(zhǔn)則就是使&均方誤差陣最小。
[0029] 在本發(fā)明一個較佳實施例中，所述卡爾曼濾波器的卡爾曼濾波方程中狀態(tài)協(xié)方差 估計的公式為：
[0030] Pk=(I-KkHk)Pk/k-1 (1-8)
[0031] 或
[0032] h.=、卜 - K ky + K 丨'RkKl (1 -9),
[0033] 協(xié)方差誤差估計的公式用式（1-8)和式（1-9)，式（1-8)的計算量較小，但式（1-9) 能始終保證算出的Pk是對稱的，P k的對稱性是使卡爾曼濾波穩(wěn)定工作的前提，所以雖然式 (1-8)計算量小，但為了避免Pk的不對稱造成的系統(tǒng)不穩(wěn)定，通常使用式(1-9)。
[0034] 在本發(fā)明一個較佳實施例中，所述卡爾曼濾波器的卡爾曼濾波方程中濾波增益方 程的公式為：
[0035] Kt ^ +Rk)1
[0036] 式(1-10)中，Pk/H稱為一步預(yù)測均方差陣或一步預(yù)測誤差方差陣：
[0037] (143)j
[0038] 由式（1-10)得到Rk和Kk成反比：如果大，Kk值就小；若Rk小，Kk值就大，Pk-i是戈的 均方誤差陣，即：
[0039] 乃-i =丑{尤-π 足―ι}
[0040] 式中足尤說―1 = A -之《-1，為之/4-!的估計誤差。
[0041] 本發(fā)明的有益效果是:本發(fā)明的提高GNSS實時監(jiān)測數(shù)據(jù)精度的算法，采用卡爾曼 濾波器通過利用狀態(tài)方程和量測方程來估計出整個系統(tǒng)在每一個釆樣間隔的所有狀態(tài)，通 過對各個狀態(tài)的統(tǒng)計與預(yù)測，反復(fù)更新，最終得到濾波結(jié)果，在計算機中其濾波方程的計算 更新過程包含了以下兩個方面:量測方程的更新和時間方程的更新，其基本原理為用前一 次的估計值遞推當(dāng)前時刻的狀態(tài)值，然后用當(dāng)前時刻的量測值修正狀態(tài)值，可以大大提高 GNSS實時監(jiān)測數(shù)據(jù)的精度。
【具體實施方式】
[0042]下面將對本發(fā)明實施例中的技術(shù)方案進行清楚、完整地描述，顯然，所描述的實施 例僅是本發(fā)明的一部分實施例，而不是全部的實施例。基于本發(fā)明中的實施例，本領(lǐng)域普通 技術(shù)人員在沒有做出創(chuàng)造性勞動前提下所獲得的所有其它實施例，都屬于本發(fā)明保護的范 圍。
[0043]本發(fā)明實施例包括：
[0044] -種提高GNSS實時監(jiān)測數(shù)據(jù)精度的算法，其特征在于，采用卡爾曼濾波器通過利 用狀態(tài)方程和量測方程來計算出整個系統(tǒng)在每一個釆樣間隔的所有狀態(tài)，通過對各個狀態(tài) 的統(tǒng)計與預(yù)測，反復(fù)更新，最終得到濾波結(jié)果，
[0045] 其中，連續(xù)系統(tǒng)的離散化后的狀態(tài)方程和量測方程分別為：
[0046] Xk=?k/k-iXk-i+rk-iffk-i (1-1),
[0047] Lk = HkXk+Vk (1-2)，
[0048] 式（1 -1)和（1 -2)中，Xk和Xk-i分別為tk時刻和tk-i時刻的（η X 1)系統(tǒng)的狀態(tài)向量； Wk-1為tk-1時刻(nXl)動態(tài)噪聲向量；?k/k-1為時間tk-1至tk的系統(tǒng)狀態(tài)轉(zhuǎn)移矩陣(nXn);Hk 為tk時刻的(mXn)量測矩陣；rk-Atk-i時刻系統(tǒng)的噪聲矩陣;Lk為tk時刻的(mXn)對系統(tǒng) 的觀測向量;Vk為tk時刻的(m X 1)維量測噪聲向量；
[0049]同時，所述狀態(tài)方程和量測方程中的Wk-dPV滿足：
[0051]式（1-3)中，正定矩陣Rk為量測噪聲方差序列的方差陣；非負定矩陣Qk為系統(tǒng)噪聲 序列的方差陣;是滿足如下條件：
[0053] 根據(jù)式（1-1)、（1-2)、（1-3)和（1-4)所描述的條件即可由Lk和Xh求得X k，即k時刻 的最優(yōu)估計值。
[0054]卡爾曼濾波器的卡爾曼濾波方程如下所示：
[0055] 狀態(tài)預(yù)報：尤1 = φ?'/卜 1 尤-1 (1-5);
[0056] 狀態(tài)協(xié)方差陣預(yù)報:= φ*-】巧-1 + rViUk 0-6),
[0057] 狀態(tài)估計:尤=之/卜].(1 .-.7);
[0058] 狀態(tài)協(xié)方差估計:Pk=(I-KkHk)Pk/k-1 (1-8)；
[0059] 或私二(I.. - κ!βdpi"人1 - κιβJ + ΚΑ'Κ!' (1-9):
[0060] 濾波增益方程:& +氏-)-1 (140)。
[0061 ]在式（1-5)中，是Xh的卡爾曼濾波估值，是通過Xh而算得的一步預(yù)測， 通過k-Ι時刻以及之前時刻的量測值得到的。
[0062] 式(1-7)用來計算估值Xk。它是根據(jù)測量值Lk和一步預(yù)測，計算出來的，式中，
[0063] Hk ~HkXklk_, = //,X, + Vk ~HkXklk_, = HkXk<k_{ + Vk (1_11)?
[0064] 其中是一步預(yù)測誤差，它的計算公式為：
[0065] 文kik-\=Xk-hi (1-12)；
[0066] 巧名^為量測值Hk的一步預(yù)測，而盡-//乂^即為量測的一步預(yù)測誤差，由式 (1 -11)可知，巧-是由和vk計算得到，由式（1 -12)可知，f是由i;/iM和父!^ 得到，所以稱//d ,為新息；
[0067] Kk稱為濾波增益矩陣。選取Kk的準(zhǔn)則就是使&均方誤差陣最小。
[0068] 式(1-10)中，Pk/H稱為一步預(yù)測均方差陣或一步預(yù)測誤差方差陣：
[0070]由式α-io)很容易得到Rk和Kk成反比：如果大，Kk值就小;若Rk小，Kk值就大。Pk-!是 的均方誤差陣，即：
[0072]式中足尤/i-1= ? - 尤/Μ，:.為尤的估計誤差。
[0073]由式（1-6)可知Pk/k-^Pk-加上噪聲方差得到的。
[0074]協(xié)方差誤差估計可以用式（1-8)和式（1-9)，式（1-8)的計算量較小，但式（1-9)能 始終保證算出的Pk是對稱的。Pk的對稱性是使卡爾曼濾波穩(wěn)定工作的前提，所以雖然式（1-8)計算量小，但為了避免P k的不對稱造成的系統(tǒng)不穩(wěn)定，我們還是常使用式(1-9)。
[0075]實施例：以一個礦區(qū)變形監(jiān)測為例：
[0076]當(dāng)?shù)V區(qū)開采工作進行到一定的階段，井下巖體的移動和破壞逐步波及到地表，使 因采動影響的地表從原有的標(biāo)高向下沉降，從而在采空區(qū)上方形成一個比采空區(qū)大得多的 沉陷區(qū)域。沉陷區(qū)域改變了原有的地表形態(tài)，引起地表標(biāo)高、水平位置發(fā)生變化，常用的定 量表述這種變化的指標(biāo)有:下沉、水平移動、傾斜、曲率、水平變形、下沉速度。
[0077] (1)、m次觀測時η點的下沉公式為：
[0078] = H，>, ~H.，" - mm ( [ " 1 5) ?
[0079] 式中，Wn表示n點的下沉至;分別為首次和m次觀測時n點的高程。
[0080] ⑵、相鄰兩點間的傾斜公式為：
[0082] 式中，1η/η+ι表不η號點至n+1號點的水平距離;W n+i，Wn分別表不n+1號點和η號點的 下沉值。
[0083] (3)、η號點附近的曲率，即η-1號點至η+1號點之間的曲率：
[0085] 式中，in+l/n、in/n-l分別表不η+1號點至η號點和η號點至η-1號點的傾斜，ln+1/n、 1 η/η-ι分別表不η+1號點至η號點和η號點至η-1號點的水平距離。
[0086] (4)、η號點的水平移動公式：
[0087] Un=Lnm-Ln〇，mm (1-18)，
[0088] 式中，Un表示η號點的水平移動;Lnm、LnQ分別表示m次觀測時和首次觀測時η號點至 觀測線控制點的水平距離，用點間距離累加求得。
[0089] (5 )、η號點至η +1號點間的水平變形公式：
[0091 ] 式中，（ln+i/n)o、（ln+i/n)m分別表不n+1號點至η號點在首次觀測時和m次觀測時的水 平距離。
[0092] (6)、n號點的下沉速度公式：
[0094] 式中，Wnm-1、Wnm分別表不m-1次和m次觀測時（即前、后兩次觀測)η點的下沉值;t表 示兩次觀測的間隔天數(shù)。
[0095] 綜上所述，本發(fā)明的提高GNSS實時監(jiān)測數(shù)據(jù)精度的算法，采用卡爾曼濾波器通過 利用狀態(tài)方程和量測方程來估計出整個系統(tǒng)在每一個釆樣間隔的所有狀態(tài)，通過對各個狀 態(tài)的統(tǒng)計與預(yù)測，反復(fù)更新，最終得到濾波結(jié)果，在計算機中其濾波方程的計算更新過程包 含了以下兩個方面:量測方程的更新和時間方程的更新，其基本原理為用前一次的估計值 遞推當(dāng)前時刻的狀態(tài)值，然后用當(dāng)前時刻的量測值修正狀態(tài)值，可以大大提高GNSS實時監(jiān) 測數(shù)據(jù)的精度。
[0096] 以上所述僅為本發(fā)明的實施例，并非因此限制本發(fā)明的專利范圍，凡是利用本發(fā) 明說明書內(nèi)容所作的等效結(jié)構(gòu)或等效流程變換，或直接或間接運用在其它相關(guān)的技術(shù)領(lǐng) 域，均同理包括在本發(fā)明的專利保護范圍內(nèi)。
【主權(quán)項】
1. 一種提高GNSS實時監(jiān)測數(shù)據(jù)精度的算法，其特征在于，采用卡爾曼濾波器通過利用 狀態(tài)方程和量測方程來計算出整個系統(tǒng)在每一個釆樣間隔的所有狀態(tài)，通過對各個狀態(tài)的 統(tǒng)計與預(yù)測，反復(fù)更新，最終得到濾波結(jié)果， 其中，連續(xù)系統(tǒng)的離散化后的狀態(tài)方程和量測方程分別為： Xk= ?k/k-lXk-l+ Γ k-lffk-1 (卜1)， Lk = HkXk+Vk (1-2)， 式（1 -1)和（1 -2)中，Xk和Xh分別為tk時刻和tk-i時刻的(η X 1)系統(tǒng)的狀態(tài)向量;Wk-i為 tk-i時刻(nXl)動態(tài)噪聲向量；?k/k-1為時間tk-1至tk的系統(tǒng)狀態(tài)轉(zhuǎn)移矩陣(nXn);Hk為tk時 亥IJ的(m X η)量測矩陣；Γ k-Atk-i時刻系統(tǒng)的噪聲矩陣;Lk為tk時刻的(m X η)對系統(tǒng)的觀測 向量;Vk為tk時刻的(m X 1)維量測噪聲向量； 同時，所述狀態(tài)方程和量測方程中的Wk-dPV滿足：式（1-3)中，正定矩陣Rk為量測噪聲方差序列的方差陣;非負定矩陣Qk為系統(tǒng)噪聲序列 的方差陣;是滿足如下條件：根據(jù)式（1-1)、（1-2)、（1-3)和（1-4)所描述的條件即可由Lk和Xh求得Xk，即k時刻的最 優(yōu)估計值。2. 根據(jù)權(quán)利要求1所述的提高GNSS實時監(jiān)測數(shù)據(jù)精度的算法，其特征在于，所述卡爾曼 濾波器的卡爾曼濾波方程中狀態(tài)預(yù)報的公式為：(1-5)， 式（1-5)中，是乂!^的卡爾曼濾波估值，是通過Xh而算得的一步預(yù)測，8卩 通過k-1時刻以及之前時刻的量測值得到的。3. 根據(jù)權(quán)利要求2所述的提高GNSS實時監(jiān)測數(shù)據(jù)精度的算法，其特征在于，所述卡爾曼 濾波器的卡爾曼濾波方程中狀態(tài)協(xié)方差陣預(yù)報的公式為：(1-6)， 由式（1-6)，可知Pk/k-i是Pk-i加上噪聲方差得到的。4. 根據(jù)權(quán)利要求3所述的提高GNSS實時監(jiān)測數(shù)據(jù)精度的算法，其特征在于，所述卡爾曼 濾波器的卡爾曼濾波方程中狀態(tài)估計的公式為：(1-7), 式（1-7)用來計算估值Xk，Xk是根據(jù)測量值Lk和一步預(yù)測義^^計算出來的，式中，(1-11), 其中，是一步預(yù)測誤差，它的計算公式為:1 技為量測值恥的一步預(yù)測，而/4 -馬尤/w即為量測的一步預(yù)測誤差，由式（1 _ 11)可知，是由尤ijPVk計算得到；由式（1-12)可知，戈由fjn^PXk得 到，所以稱/-/,-% A's 4為新息； Kk稱為濾波增益矩陣，選取Kk的準(zhǔn)則就是使fA均方誤差陣最小。5. 根據(jù)權(quán)利要求4所述的提高GNSS實時監(jiān)測數(shù)據(jù)精度的算法，其特征在于，所述卡爾曼 濾波器的卡爾曼濾波方程中狀態(tài)協(xié)方差估計的公式為： Pk=(I_KkHk)Pk/k-1 (1-8) 或 ' '' ......... fV …，l ，l 'l V1-7/，. 協(xié)方差誤差估計的公式用式（1-8)和式（1-9)，式（1-8)的計算量較小，但式（1-9)能始 終保證算出的Pk是對稱的，Pk的對稱性是使卡爾曼濾波穩(wěn)定工作的前提，所以雖然式（1-8) 計算量小，但為了避免Pk的不對稱造成的系統(tǒng)不穩(wěn)定，通常使用式(1-9)。6. 根據(jù)權(quán)利要求5所述的提高GNSS實時監(jiān)測數(shù)據(jù)精度的算法，其特征在于，所述卡爾曼 濾波器的卡爾曼濾波方程中濾波增益方程的公式為：(1-10)， 式(1-10)中，Pk/k-i稱為一步預(yù)測均方差陣或一步預(yù)測誤差方差陣：(1-13),由式（1-10)得到Rk和Kk成反比：如果大，Kk值就小;若Rk小，K k值就大，？1{-1是爲(wèi)_1的均方誤 差陣，即： ^ (1-14)， 式中足Λ w =? -名^，為矣/M的估計誤差。
【文檔編號】G01S19/37GK105866807SQ201610202484
【公開日】2016年8月17日
【申請日】2016年4月5日
【發(fā)明人】盧康
【申請人】南信大影像技術(shù)工程（蘇州）有限公司