一种提高gnss实时监测数据精度的算法
【专利摘要】本发明公开了一种提高GNSS实时监测数据精度的算法，采用卡尔曼滤波器通过利用状态方程和量测方程来估计出整个系统在每一个釆样间隔的所有状态，通过对各个状态的统计与预测，反复更新，最终得到滤波结果。在计算机中其滤波方程的计算更新过程包含了以下两个方面：量测方程的更新和时间方程的更新，其基本原理为用前一次的估计值递推当前时刻的状态值，然后用当前时刻的量测值修正状态值，可以大大提高GNSS实时监测数据的精度。
【专利说明】
一种提高GNSS实时监测数据精度的算法
技术领域
[0001] 本发明涉及计算机数据处理算法领域，尤其涉及一种提高GNSS实时监测数据精度 的算法。
【背景技术】
[0002] 卡尔曼(Kalman)滤波是由卡尔曼(Kalman R E)博士为满足计算机用于人造地球 卫星定轨和导航等计算的要求，在1960年提出的一种新的线性滤波模型，通常称它为卡尔 曼滤波方法。卡尔曼滤波是一种具有无偏性的递推线性最小方差估计，即估计误差的均值 或数学期望为零。在计算方法上，卡尔曼滤波采用递推形式，即在t-Ι时刻估值的基础上，利 用t时刻的观测值，递推得倒t时刻的状态估值。由于一次仅处理一个时刻的观测值，无需存 储先前的观测数据，因而计算量大大减少。
[0003] 卡尔曼滤波能有效的消除噪声干扰，获得逼近真实情况的有用信息。这种滤波方 法很适合处理动态系统的观测数据，因而广泛应用于动态测量。但是卡尔曼滤波是一个连 续系统，通过离散化后，才能让计算机实现卡尔曼滤波的模拟，这里通过一种简单方法的是 卡尔曼滤波离散化。然而这种滤波离散化并不能有效地提高GNSS实时监测数据的精度。因 此，要对该滤波进行改进，以提高GNSS实时监测数据的精度。
【发明内容】
[0004] 本发明主要解决的技术问题是提供一种提高GNSS实时监测数据精度的算法，采 用卡尔曼滤波器通过利用状态方程和量测方程来估计出整个系统在每一个釆样间隔的所 有状态，通过对各个状态的统计与预测，反复更新，最终得到滤波结果，在计算机中其滤波 方程的计算更新过程包含了以下两个方面:量测方程的更新和时间方程的更新，其基本原 理为用前一次的估计值递推当前时刻的状态值，然后用当前时刻的量测值修正状态值，可 以大大提高GNSS实时监测数据的精度。
[0005] 为解决上述技术问题，本发明采用的一个技术方案是:提供了一种提高GNSS实时 监测数据精度的算法，采用卡尔曼滤波器通过利用状态方程和量测方程来计算出整个系统 在每一个釆样间隔的所有状态，通过对各个状态的统计与预测，反复更新，最终得到滤波结 果，
[0006] 其中，连续系统的离散化后的状态方程和量测方程分别为：
[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)系统的状态向量； Wk-1为tk-i时刻(nXl)动态噪声向量；?k/k-i为时间tk-iStk的系统状态转移矩阵(nXn);H k 为tk时刻的(mXn)量测矩阵；刻系统的噪声矩阵;Lk为tk时刻的(mXn)对系统 的观测向量;V k为tk时刻的(m X 1)维量测噪声向量；
[0010] 同时，所述状态方程和量测方程中的UPV满足：
[0012]式（1-3)中，正定矩阵Rk为量测噪声方差序列的方差阵；非负定矩阵Qk为系统噪声 序列的方差阵;是满足如下条件：
[0014] 根据式（1-1)、（1-2)、（1-3)和（1-4)所描述的条件即可由Lk和XH求得xk，g卩k时刻 的最优估计值。
[0015] 在本发明一个较佳实施例中，所述卡尔曼滤波器的卡尔曼滤波方程中状态预报的 公式为：
[0016] = φ?/?-ι^Α-ι (1-5),
[0017] 式（1-5)中，是Χη的卡尔曼滤波估值，li/iM是通过Χη而算得的一步预测，即 通过k-Ι时刻以及之前时刻的量测值得到的。
[0018] 在本发明一个较佳实施例中，所述卡尔曼滤波器的卡尔曼滤波方程中状态协方差 阵预报的公式为：
[0019] 乃-1 =Φ?'?-4 乃-1--必--t (.1-6)，.
[0020] 由式（1-6)，可知Pk/k-i是Pk-i加上噪声方差得到的。
[0021 ]在本发明一个较佳实施例中，所述卡尔曼滤波器的卡尔曼滤波方程中状态估计的 公式为：
[0022] 之-之 -丑it 之n) :(1-7),
[0023] 式(1-7)用来计算估值Xk，Xk是根据测量值Lk和一步预测计算出来的，式中，
[0024] ^HkXkl]i_^ = HkXk + V^-HkXkl]i^ = HkXk ik-\ + K (1-11)，
[0025] 其中，名^^是一步预测误差，它的计算公式为：
[0026] =Xj. - Xk /k-l (1-12):
[0027] 为量测值Hk的一步预测，而私-/4之,^即为量测的一步预测误差，由式 (1-11)可知，/4 -/4?·Μ4是由尤和Vk计算得到；由式（1-以)可知，矣是由尤i和知 得到，所以称巧-为新息；
[0028] Kk称为滤波增益矩阵，选取Kk的准则就是使&均方误差阵最小。
[0029] 在本发明一个较佳实施例中，所述卡尔曼滤波器的卡尔曼滤波方程中状态协方差 估计的公式为：
[0030] Pk=(I-KkHk)Pk/k-1 (1-8)
[0031] 或
[0032] h.=、卜 - K ky + K 丨'RkKl (1 -9),
[0033] 协方差误差估计的公式用式（1-8)和式（1-9)，式（1-8)的计算量较小，但式（1-9) 能始终保证算出的Pk是对称的，P k的对称性是使卡尔曼滤波稳定工作的前提，所以虽然式 (1-8)计算量小，但为了避免Pk的不对称造成的系统不稳定，通常使用式(1-9)。
[0034] 在本发明一个较佳实施例中，所述卡尔曼滤波器的卡尔曼滤波方程中滤波增益方 程的公式为：
[0035] Kt ^ +Rk)1
[0036] 式(1-10)中，Pk/H称为一步预测均方差阵或一步预测误差方差阵：
[0037] (143)j
[0038] 由式（1-10)得到Rk和Kk成反比：如果大，Kk值就小；若Rk小，Kk值就大，Pk-i是戈的 均方误差阵，即：
[0039] 乃-i =丑{尤-π 足―ι}
[0040] 式中足尤说―1 = A -之《-1，为之/4-!的估计误差。
[0041] 本发明的有益效果是:本发明的提高GNSS实时监测数据精度的算法，采用卡尔曼 滤波器通过利用状态方程和量测方程来估计出整个系统在每一个釆样间隔的所有状态，通 过对各个状态的统计与预测，反复更新，最终得到滤波结果，在计算机中其滤波方程的计算 更新过程包含了以下两个方面:量测方程的更新和时间方程的更新，其基本原理为用前一 次的估计值递推当前时刻的状态值，然后用当前时刻的量测值修正状态值，可以大大提高 GNSS实时监测数据的精度。
【具体实施方式】
[0042]下面将对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施 例仅是本发明的一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通 技术人员在没有做出创造性劳动前提下所获得的所有其它实施例，都属于本发明保护的范 围。
[0043]本发明实施例包括：
[0044] -种提高GNSS实时监测数据精度的算法，其特征在于，采用卡尔曼滤波器通过利 用状态方程和量测方程来计算出整个系统在每一个釆样间隔的所有状态，通过对各个状态 的统计与预测，反复更新，最终得到滤波结果，
[0045] 其中，连续系统的离散化后的状态方程和量测方程分别为：
[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)系统的状态向量； Wk-1为tk-1时刻(nXl)动态噪声向量；?k/k-1为时间tk-1至tk的系统状态转移矩阵(nXn);Hk 为tk时刻的(mXn)量测矩阵；rk-Atk-i时刻系统的噪声矩阵;Lk为tk时刻的(mXn)对系统 的观测向量;Vk为tk时刻的(m X 1)维量测噪声向量；
[0049]同时，所述状态方程和量测方程中的Wk-dPV满足：
[0051]式（1-3)中，正定矩阵Rk为量测噪声方差序列的方差阵；非负定矩阵Qk为系统噪声 序列的方差阵;是满足如下条件：
[0053] 根据式（1-1)、（1-2)、（1-3)和（1-4)所描述的条件即可由Lk和Xh求得X k，即k时刻 的最优估计值。
[0054]卡尔曼滤波器的卡尔曼滤波方程如下所示：
[0055] 状态预报：尤1 = φ?'/卜 1 尤-1 (1-5);
[0056] 状态协方差阵预报:= φ*-】巧-1 + rViUk 0-6),
[0057] 状态估计:尤=之/卜].(1 .-.7);
[0058] 状态协方差估计:Pk=(I-KkHk)Pk/k-1 (1-8)；
[0059] 或私二(I.. - κ!βdpi"人1 - κιβJ + ΚΑ'Κ!' (1-9):
[0060] 滤波增益方程:& +氏-)-1 (140)。
[0061 ]在式（1-5)中，是Xh的卡尔曼滤波估值，是通过Xh而算得的一步预测， 通过k-Ι时刻以及之前时刻的量测值得到的。
[0062] 式(1-7)用来计算估值Xk。它是根据测量值Lk和一步预测，计算出来的，式中，
[0063] Hk ~HkXklk_, = //,X, + Vk ~HkXklk_, = HkXk<k_{ + Vk (1_11)?
[0064] 其中是一步预测误差，它的计算公式为：
[0065] 文kik-\=Xk-hi (1-12)；
[0066] 巧名^为量测值Hk的一步预测，而尽-//乂^即为量测的一步预测误差，由式 (1 -11)可知，巧-是由和vk计算得到，由式（1 -12)可知，f是由i;/iM和父!^ 得到，所以称//d ,为新息；
[0067] Kk称为滤波增益矩阵。选取Kk的准则就是使&均方误差阵最小。
[0068] 式(1-10)中，Pk/H称为一步预测均方差阵或一步预测误差方差阵：
[0070]由式α-io)很容易得到Rk和Kk成反比：如果大，Kk值就小;若Rk小，Kk值就大。Pk-!是 的均方误差阵，即：
[0072]式中足尤/i-1= ? - 尤/Μ，:.为尤的估计误差。
[0073]由式（1-6)可知Pk/k-^Pk-加上噪声方差得到的。
[0074]协方差误差估计可以用式（1-8)和式（1-9)，式（1-8)的计算量较小，但式（1-9)能 始终保证算出的Pk是对称的。Pk的对称性是使卡尔曼滤波稳定工作的前提，所以虽然式（1-8)计算量小，但为了避免P k的不对称造成的系统不稳定，我们还是常使用式(1-9)。
[0075]实施例：以一个矿区变形监测为例：
[0076]当矿区开采工作进行到一定的阶段，井下岩体的移动和破坏逐步波及到地表，使 因采动影响的地表从原有的标高向下沉降，从而在采空区上方形成一个比采空区大得多的 沉陷区域。沉陷区域改变了原有的地表形态，引起地表标高、水平位置发生变化，常用的定 量表述这种变化的指标有:下沉、水平移动、倾斜、曲率、水平变形、下沉速度。
[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表 示两次观测的间隔天数。
[0095] 综上所述，本发明的提高GNSS实时监测数据精度的算法，采用卡尔曼滤波器通过 利用状态方程和量测方程来估计出整个系统在每一个釆样间隔的所有状态，通过对各个状 态的统计与预测，反复更新，最终得到滤波结果，在计算机中其滤波方程的计算更新过程包 含了以下两个方面:量测方程的更新和时间方程的更新，其基本原理为用前一次的估计值 递推当前时刻的状态值，然后用当前时刻的量测值修正状态值，可以大大提高GNSS实时监 测数据的精度。
[0096] 以上所述仅为本发明的实施例，并非因此限制本发明的专利范围，凡是利用本发 明说明书内容所作的等效结构或等效流程变换，或直接或间接运用在其它相关的技术领 域，均同理包括在本发明的专利保护范围内。
【主权项】
1. 一种提高GNSS实时监测数据精度的算法，其特征在于，采用卡尔曼滤波器通过利用 状态方程和量测方程来计算出整个系统在每一个釆样间隔的所有状态，通过对各个状态的 统计与预测，反复更新，最终得到滤波结果， 其中，连续系统的离散化后的状态方程和量测方程分别为： 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)系统的状态向量;Wk-i为 tk-i时刻(nXl)动态噪声向量；?k/k-1为时间tk-1至tk的系统状态转移矩阵(nXn);Hk为tk时 亥IJ的(m X η)量测矩阵；Γ k-Atk-i时刻系统的噪声矩阵;Lk为tk时刻的(m X η)对系统的观测 向量;Vk为tk时刻的(m X 1)维量测噪声向量； 同时，所述状态方程和量测方程中的Wk-dPV满足：式（1-3)中，正定矩阵Rk为量测噪声方差序列的方差阵;非负定矩阵Qk为系统噪声序列 的方差阵;是满足如下条件：根据式（1-1)、（1-2)、（1-3)和（1-4)所描述的条件即可由Lk和Xh求得Xk，即k时刻的最 优估计值。2. 根据权利要求1所述的提高GNSS实时监测数据精度的算法，其特征在于，所述卡尔曼 滤波器的卡尔曼滤波方程中状态预报的公式为：(1-5)， 式（1-5)中，是乂!^的卡尔曼滤波估值，是通过Xh而算得的一步预测，8卩 通过k-1时刻以及之前时刻的量测值得到的。3. 根据权利要求2所述的提高GNSS实时监测数据精度的算法，其特征在于，所述卡尔曼 滤波器的卡尔曼滤波方程中状态协方差阵预报的公式为：(1-6)， 由式（1-6)，可知Pk/k-i是Pk-i加上噪声方差得到的。4. 根据权利要求3所述的提高GNSS实时监测数据精度的算法，其特征在于，所述卡尔曼 滤波器的卡尔曼滤波方程中状态估计的公式为：(1-7), 式（1-7)用来计算估值Xk，Xk是根据测量值Lk和一步预测义^^计算出来的，式中，(1-11), 其中，是一步预测误差，它的计算公式为:1 技为量测值恥的一步预测，而/4 -馬尤/w即为量测的一步预测误差，由式（1 _ 11)可知，是由尤ijPVk计算得到；由式（1-12)可知，戈由fjn^PXk得 到，所以称/-/,-% A's 4为新息； Kk称为滤波增益矩阵，选取Kk的准则就是使fA均方误差阵最小。5. 根据权利要求4所述的提高GNSS实时监测数据精度的算法，其特征在于，所述卡尔曼 滤波器的卡尔曼滤波方程中状态协方差估计的公式为： Pk=(I_KkHk)Pk/k-1 (1-8) 或 ' '' ......... fV …，l ，l 'l V1-7/，. 协方差误差估计的公式用式（1-8)和式（1-9)，式（1-8)的计算量较小，但式（1-9)能始 终保证算出的Pk是对称的，Pk的对称性是使卡尔曼滤波稳定工作的前提，所以虽然式（1-8) 计算量小，但为了避免Pk的不对称造成的系统不稳定，通常使用式(1-9)。6. 根据权利要求5所述的提高GNSS实时监测数据精度的算法，其特征在于，所述卡尔曼 滤波器的卡尔曼滤波方程中滤波增益方程的公式为：(1-10)， 式(1-10)中，Pk/k-i称为一步预测均方差阵或一步预测误差方差阵：(1-13),由式（1-10)得到Rk和Kk成反比：如果大，Kk值就小;若Rk小，K k值就大，？1{-1是爲_1的均方误 差阵，即： ^ (1-14)， 式中足Λ w =? -名^，为矣/M的估计误差。
【文档编号】G01S19/37GK105866807SQ201610202484
【公开日】2016年8月17日
【申请日】2016年4月5日
【发明人】卢康
【申请人】南信大影像技术工程（苏州）有限公司