CN111623703A - 一种基于新型卡尔曼滤波的北斗变形监测实时处理方法 - Google Patents

一种基于新型卡尔曼滤波的北斗变形监测实时处理方法 Download PDF

Info

Publication number
CN111623703A
CN111623703A CN202010733758.1A CN202010733758A CN111623703A CN 111623703 A CN111623703 A CN 111623703A CN 202010733758 A CN202010733758 A CN 202010733758A CN 111623703 A CN111623703 A CN 111623703A
Authority
CN
China
Prior art keywords
time
matrix
real
variance
monitoring
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Pending
Application number
CN202010733758.1A
Other languages
English (en)
Inventor
梁晓东
雷孟飞
周俊华
杨振武
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Hunan Lianzhi Technology Co Ltd
Original Assignee
Hunan Lianzhi Technology Co Ltd
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Hunan Lianzhi Technology Co Ltd filed Critical Hunan Lianzhi Technology Co Ltd
Priority to CN202010733758.1A priority Critical patent/CN111623703A/zh
Publication of CN111623703A publication Critical patent/CN111623703A/zh
Pending legal-status Critical Current

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01B—MEASURING LENGTH, THICKNESS OR SIMILAR LINEAR DIMENSIONS; MEASURING ANGLES; MEASURING AREAS; MEASURING IRREGULARITIES OF SURFACES OR CONTOURS
    • G01B7/00—Measuring arrangements characterised by the use of electric or magnetic techniques
    • G01B7/16—Measuring arrangements characterised by the use of electric or magnetic techniques for measuring the deformation in a solid, e.g. by resistance strain gauge
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2218/00—Aspects of pattern recognition specially adapted for signal processing
    • G06F2218/02—Preprocessing

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Mathematical Physics (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Computational Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • Theoretical Computer Science (AREA)
  • Computing Systems (AREA)
  • Algebra (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Radar Systems Or Details Thereof (AREA)

Abstract

本发明提供一种基于新型卡尔曼滤波的北斗变形监测实时处理方法，包括：获取初始状态向量的最优估值等；基于卡尔曼滤波，得到状态向量的一步预测值和一步预测方差；获取实时监测数据；得到该时间段数据的标准方差和卡尔曼增益；计算状态向量的最优估值和最优估值的方差阵；构造滑动窗口残差向量并动态调整滑动窗口的大小，实时更新观测噪声方差阵；计算位移量并输出。本发明通过残差和标准方差的比值来实现粗差探测；通过实时缩放测量噪声方差阵的值来修复粗差，减弱粗差对观测结果的污染；在粗差探测的基础上，结合实时更新观测噪声方差阵和动态调整滑动窗口的大小，加速监测结果的收敛，迅速反应监测点的真实位移，能够满足变形的监测需求。

Description

一种基于新型卡尔曼滤波的北斗变形监测实时处理方法
技术领域
本发明涉及北斗变形监测技术领域，具体涉及一种基于新型卡尔曼滤波的北斗变形监测实时处理方法。
背景技术
北斗卫星导航系统（BeiDou Navigation Satellite System，BDS）为我国自行研制的全球卫星导航系统，随着北斗三代系统的逐步建立，北斗在变形监测、定位等领域应用更加广泛。
采用北斗进行变形监测时，由于周跳、多路径效应以及接收机信号噪声的影响，观测数据中包含粗差，除此之外，北斗解算的过程中受到噪声及解算算法的影响，使得解算结果包含有高频的随机误差，降低了观测结果的准确性和稳定性，从而为后续的数据分析带来困难。
变形监测数据处理的另一个难点在于：如何将真实位移和粗差进行识别。现实的监测中，监测点可能发生了真实位移，真实的位移能导致监测结果发生较大变化，因此，需要将粗差和真实位移进行区分，不能将真实位移当做粗差校正。
现有技术中对于变形监测数据处理有很多方法，比如普遍采用滑动平均算法、多项式拟合以及arma模型等方法，均不能解决上述技术问题，具体是：
1、滑动平均算法模型简单，计算方便，对于数据的完整性要求相对较低，但是滑动平均算法对于粗差只能削弱，对粗差和随机误差的处理效果与滑动窗口有关，另外对于真实变形的反应并不准确。
2、多项式拟合算法能够较好的处理粗差和随机误差，但是对于真实变形反应并不好，并且只能用于后处理算法中，对于实时监测并不适用。
3、arma模型涉及到数据平稳性检验以及p和q值的确定，算法比较复杂，并且同样不适用与实时监测。
综上所述，设计一种操作便捷、能够实现粗差探测以减弱粗差对观测结果的影响及能够真实区分监测点的真实位移变化的北斗变形监测实时处理方法具有重要意义。
发明内容
本发明提供一种基于新型卡尔曼滤波的北斗变形监测实时处理方法，技术方案如下：
一种基于新型卡尔曼滤波的北斗变形监测实时处理方法，包括以下步骤：
步骤一、获取N天的监测结果作为初始数据集，1≤N≤3；获取初始状态向量的最优估值
Figure DEST_PATH_IMAGE002
、初始状态向量的最优估值的方差阵
Figure DEST_PATH_IMAGE004
、初始动态噪声方差阵
Figure DEST_PATH_IMAGE006
、初始观测噪声方差阵
Figure DEST_PATH_IMAGE008
；
步骤二、基于卡尔曼滤波，通过表达式6）和表达式7）分别计算得到状态向量的一步预测值
Figure DEST_PATH_IMAGE010
和状态向量的一步预测方差
Figure DEST_PATH_IMAGE012
；
Figure DEST_PATH_IMAGE014
6）；
Figure DEST_PATH_IMAGE016
7）；
其中：k为历元数，k为大于等于1且小于等于M的自然数；A为转移矩阵，
Figure DEST_PATH_IMAGE018
，
Figure DEST_PATH_IMAGE020
为历元间的时间间隔，
Figure DEST_PATH_IMAGE022
为
Figure DEST_PATH_IMAGE024
的转置；
Figure DEST_PATH_IMAGE026
为状态向量在
Figure DEST_PATH_IMAGE028
时刻的最优估值，
Figure DEST_PATH_IMAGE028A
为第k-1历元的时刻，
Figure DEST_PATH_IMAGE030
，
Figure DEST_PATH_IMAGE032
为
Figure DEST_PATH_IMAGE028AA
时刻的三维坐标，
Figure DEST_PATH_IMAGE034
为
Figure DEST_PATH_IMAGE028AAA
时刻向
Figure DEST_PATH_IMAGE036
时刻变化的速度；
Figure DEST_PATH_IMAGE038
为在
Figure DEST_PATH_IMAGE028AAAA
时刻的状态向量的最优估值的方差阵；
Figure DEST_PATH_IMAGE040
为噪声矩阵，
Figure DEST_PATH_IMAGE042
，
Figure DEST_PATH_IMAGE044
为
Figure DEST_PATH_IMAGE040A
的转置；
Figure DEST_PATH_IMAGE046
为
Figure DEST_PATH_IMAGE028AAAAA
时刻的动态噪声，
Figure DEST_PATH_IMAGE048
为
Figure DEST_PATH_IMAGE028AAAAAA
时刻的动态噪声的方差阵；
步骤三、获取
Figure DEST_PATH_IMAGE028AAAAAAA
时刻至
Figure DEST_PATH_IMAGE036A
时刻的实时监测数据
Figure DEST_PATH_IMAGE050
；结合步骤二获得的状态向量的一步预测值
Figure DEST_PATH_IMAGE010A
，通过表达式8）得到
Figure DEST_PATH_IMAGE028AAAAAAAA
时刻至
Figure DEST_PATH_IMAGE036AA
时刻该时间段数据的标准方差
Figure DEST_PATH_IMAGE052
且通过表达式9）计算卡尔曼增益
Figure DEST_PATH_IMAGE054
：
Figure DEST_PATH_IMAGE056
8）；
Figure DEST_PATH_IMAGE058
9）；
其中：B为观测系数矩阵，
Figure DEST_PATH_IMAGE060
，
Figure DEST_PATH_IMAGE062
为B的转置；
Figure DEST_PATH_IMAGE064
为
Figure DEST_PATH_IMAGE036AAA
时刻的观测噪声方差阵；
步骤四、通过表达式10）和表达式11）分别计算
Figure DEST_PATH_IMAGE036AAAA
时刻状态向量的最优估值
Figure DEST_PATH_IMAGE066
和最优估值的方差阵
Figure DEST_PATH_IMAGE068
：
Figure DEST_PATH_IMAGE070
10）；
Figure DEST_PATH_IMAGE072
11）；
其中：
Figure DEST_PATH_IMAGE074
为4×4的单位矩阵；
构造滑动窗口残差向量并动态调整滑动窗口的大小，采用表达式13）实时更新观测噪声方差阵：
Figure DEST_PATH_IMAGE076
13）；
其中：
Figure DEST_PATH_IMAGE078
为观测噪声缩小因子，
Figure DEST_PATH_IMAGE080
；
Figure DEST_PATH_IMAGE082
为观测噪声方差阵系数；
步骤五、计算位移量并输出，具体是：监测点
Figure DEST_PATH_IMAGE036AAAAA
时刻的位移量
Figure DEST_PATH_IMAGE084
，其中：
Figure DEST_PATH_IMAGE086
通过
Figure DEST_PATH_IMAGE036AAAAAA
时刻状态向量的最优估值
Figure DEST_PATH_IMAGE088
获得，
Figure DEST_PATH_IMAGE090
；取k= k+1，返回步骤二。
以上技术方案中优选的，所述步骤四中构造滑动窗口残差向量并动态调整滑动窗口的大小以及实时更新观测噪声方差阵具体包括以下步骤：
步骤①、缓存最近的m个北斗结算结果构成观测值滑动窗口，
Figure DEST_PATH_IMAGE092
且m为整数；构建滑动窗口残差向量
Figure DEST_PATH_IMAGE094
；
步骤②、通过表达式12）计算观测数据的方差作为
Figure DEST_PATH_IMAGE096
时刻至
Figure DEST_PATH_IMAGE098
时刻时间段数据的标准方差
Figure DEST_PATH_IMAGE100
：
Figure DEST_PATH_IMAGE102
12）；
其中：
Figure DEST_PATH_IMAGE104
为m个实时数据的方差，
Figure DEST_PATH_IMAGE106
；
步骤③、根据2倍标准方差计算观测噪声方差阵系数
Figure DEST_PATH_IMAGE108
；
步骤④、调整观测噪声缩小因子
Figure DEST_PATH_IMAGE078A
的值：若
Figure DEST_PATH_IMAGE110
，则取
Figure DEST_PATH_IMAGE112
；若
Figure DEST_PATH_IMAGE114
，则取
Figure DEST_PATH_IMAGE116
；采用表达式13）实时更新观测噪声方差阵；
标准方差的滑动窗口m根据
Figure DEST_PATH_IMAGE078AA
的值动态调整：若
Figure DEST_PATH_IMAGE118
，则取
Figure DEST_PATH_IMAGE120
；若
Figure DEST_PATH_IMAGE122
，则取
Figure DEST_PATH_IMAGE124
，
Figure DEST_PATH_IMAGE126
代表
Figure DEST_PATH_IMAGE128
取整；返回步骤①。
以上技术方案中优选的，所述步骤一中：
初始状态向量的最优估值
Figure DEST_PATH_IMAGE130
的获得具体是：
Figure DEST_PATH_IMAGE132
，初始数据集中的监测数据总数为
Figure DEST_PATH_IMAGE134
，利用最小二乘计算
Figure DEST_PATH_IMAGE134A
个监测数据的最优估计值作为
Figure DEST_PATH_IMAGE136
，
Figure DEST_PATH_IMAGE138
，
Figure DEST_PATH_IMAGE140
为单位矩阵，
Figure DEST_PATH_IMAGE142
为
Figure DEST_PATH_IMAGE140A
的转置，
Figure DEST_PATH_IMAGE144
为观测值向量，P为观测值向量
Figure DEST_PATH_IMAGE144A
的协方差矩阵，数值根据北斗解算软件输出的Ratio值确定，
Figure DEST_PATH_IMAGE146
，R为
Figure DEST_PATH_IMAGE134AA
个监测数据的Ratio值矩阵，
Figure DEST_PATH_IMAGE148
为R中所有元素之和；
Figure DEST_PATH_IMAGE150
为初始速度，
Figure DEST_PATH_IMAGE152
，
Figure DEST_PATH_IMAGE154
为第1个监测数据；
初始状态向量的最优估值的方差阵
Figure DEST_PATH_IMAGE156
的获得具体是：
Figure DEST_PATH_IMAGE158
，
Figure DEST_PATH_IMAGE160
，
Figure DEST_PATH_IMAGE162
为初始数据集中的第j个监测数据，j为大于0的自然数；
Figure DEST_PATH_IMAGE164
；
初始动态噪声方差阵
Figure DEST_PATH_IMAGE166
的获得：
Figure DEST_PATH_IMAGE168
，初始观测噪声方差阵
Figure DEST_PATH_IMAGE170
为北斗解算后的残差方差值。
应用本发明的技术方案，具体是：
1、本发明采用新型的卡尔曼滤波通过残差和标准方差的比值来实现粗差探测；通过实时缩放测量噪声方差阵的值来修复粗差，减弱粗差对观测结果的污染。
2、本发明在粗差探测的基础上，结合实时更新观测噪声方差阵和动态调整滑动窗口的大小，可以加速监测结果的收敛，迅速反应监测点的真实位移，延迟时间较小，能够满足变形的监测需求。
除了上面所描述的目的、特征和优点之外，本发明还有其它的目的、特征和优点。下面将参照图，对本发明作进一步详细的说明。
附图说明
构成本申请的一部分的附图用来提供对本发明的进一步理解，本发明的示意性实施例及其说明用于解释本发明，并不构成对本发明的不当限定。在附图中：
图1是实施例中基于新型卡尔曼滤波的北斗变形监测实时处理方法的流程示意图；
图2是现有技术中普通卡尔曼滤波处理结果示意图；
图3是本实施例中新型卡尔曼滤波处理结果示意图。
具体实施方式
以下结合附图对本发明的实施例进行详细说明，但是本发明可以根据权利要求限定和覆盖的多种不同方式实施。
实施例：
一种基于新型卡尔曼滤波的北斗变形监测实时处理方法，具体是：
变形监测的结果为关于时间t的坐标序列，k为历元数（k>0），其中
Figure DEST_PATH_IMAGE036AAAAAAA
为第k历元的时刻，因此，
Figure DEST_PATH_IMAGE036AAAAAAAA
时刻的三维坐标
Figure DEST_PATH_IMAGE086A
采用表达式1）表示：
Figure DEST_PATH_IMAGE172
1）；
其中：
Figure DEST_PATH_IMAGE174
为
Figure DEST_PATH_IMAGE028AAAAAAAAA
时刻的三维坐标；
Figure DEST_PATH_IMAGE034A
和
Figure DEST_PATH_IMAGE176
分别为
Figure DEST_PATH_IMAGE028AAAAAAAAAA
时刻向
Figure 100002_DEST_PATH_IMAGE036AAAAAAAAA
变化的速度和加速度；
Figure DEST_PATH_IMAGE020A
为历元间的时间间隔，与接收机采样频率有关，本实施例中接收机采用15s采样间隔，
Figure 100002_DEST_PATH_IMAGE020AA
=15。
由于变形监测中监测物变化缓慢，所以在变形过程中只考虑
Figure DEST_PATH_IMAGE174A
和
Figure DEST_PATH_IMAGE034AA
两项，
Figure DEST_PATH_IMAGE176A
作为滤波中的动态噪声。
系统状态方程为表达式2）：
Figure DEST_PATH_IMAGE178
2）；
其中：
Figure DEST_PATH_IMAGE180
和
Figure DEST_PATH_IMAGE026A
为
Figure 100002_DEST_PATH_IMAGE036AAAAAAAAAA
时刻和
Figure DEST_PATH_IMAGE028AAAAAAAAAAA
时刻状态向量的最优估值；A为转移矩阵，
Figure DEST_PATH_IMAGE018A
；
Figure DEST_PATH_IMAGE040AA
为噪声矩阵，
Figure DEST_PATH_IMAGE042A
；
Figure DEST_PATH_IMAGE046A
为
Figure DEST_PATH_IMAGE028AAAAAAAAAAAA
时刻的动态噪声，
Figure DEST_PATH_IMAGE182
；
观测方程为表达式3）：
Figure DEST_PATH_IMAGE184
3）；
其中：
Figure DEST_PATH_IMAGE186
为
Figure 100002_DEST_PATH_IMAGE036AAAAAAAAAAA
时刻监测数据（即监测结果）；B为观测系数矩阵，
Figure DEST_PATH_IMAGE188
；
Figure DEST_PATH_IMAGE190
为观测噪声。
本实施例基于新型卡尔曼滤波的北斗变形监测实时处理方法，具体流程如图1所示，具体包括如下步骤：
第一步、获取N天的监测结果，得到初始状态向量的最优估值
Figure DEST_PATH_IMAGE002A
、初始状态向量的最优估值的方差阵
Figure DEST_PATH_IMAGE004A
、初始动态噪声方差阵
Figure DEST_PATH_IMAGE006A
、初始观测噪声方差阵
Figure DEST_PATH_IMAGE008A
；1≤N≤3（本实施例中取1天的监测结果进行计算）。
初始状态向量的最优估值
Figure DEST_PATH_IMAGE002AA
的获得具体是：
积累第一天的监测数据记为初始数据集，初始数据集中的监测数据的总数为
Figure 100002_DEST_PATH_IMAGE134AAA
，初始状态向量为
Figure DEST_PATH_IMAGE192
，则
Figure DEST_PATH_IMAGE132A
，
Figure DEST_PATH_IMAGE194
为初始三维坐标（即监测数据），
Figure DEST_PATH_IMAGE150A
为初始速度。利用最小二乘计算
Figure DEST_PATH_IMAGE134AAAA
个监测数据的最优估计值作为
Figure DEST_PATH_IMAGE196
，如表达式4）：
Figure DEST_PATH_IMAGE138A
4）；
其中：
Figure DEST_PATH_IMAGE140AA
为3×3的单位矩阵，
Figure DEST_PATH_IMAGE142A
为
Figure DEST_PATH_IMAGE198
的转置，
Figure 100002_DEST_PATH_IMAGE144AA
为观测值向量，P为观测值向量
Figure DEST_PATH_IMAGE144AAA
的协方差矩阵，数值根据北斗解算软件输出的Ratio值确定，
Figure DEST_PATH_IMAGE146A
，R为
Figure DEST_PATH_IMAGE134AAAAA
个监测数据的Ratio值矩阵，
Figure DEST_PATH_IMAGE148A
为R中所有元素之和；
Figure 100002_DEST_PATH_IMAGE150AA
通过表达式5）获得：
Figure DEST_PATH_IMAGE152A
5）；
其中：
Figure DEST_PATH_IMAGE154A
为初始状态向量计算后的第1个监测结果；
Figure DEST_PATH_IMAGE200
为初始三维坐标（即初始监测数据）；
Figure DEST_PATH_IMAGE202
为历元间的时间间隔。
初始状态向量的最优估值的方差阵
Figure DEST_PATH_IMAGE156A
的获得具体是：
Figure DEST_PATH_IMAGE158A
，其中：
Figure 100002_DEST_PATH_IMAGE160A
，
Figure 100002_DEST_PATH_IMAGE162A
为初始数据集中的第j个监测数据，j为大于0的自然数；根据误差传播定律，由初始速度的计算公式可以得到初始速度的初始方差，
Figure 100002_DEST_PATH_IMAGE164A
。
初始动态噪声方差阵
Figure DEST_PATH_IMAGE166A
的获得：
Figure DEST_PATH_IMAGE168A
，初始观测噪声方差阵
Figure DEST_PATH_IMAGE170A
为3×3矩阵，数值为北斗解算后的XYZ三个方向结果的残差方差以及协方差值。
第二步、基于卡尔曼滤波，通过表达式6）和表达式7）分别计算得到状态向量的一步预测值
Figure DEST_PATH_IMAGE010AA
和状态向量的一步预测方差
Figure DEST_PATH_IMAGE012A
；
Figure DEST_PATH_IMAGE014A
6）；
Figure DEST_PATH_IMAGE016A
7）；
其中：
Figure DEST_PATH_IMAGE090A
，
Figure 100002_DEST_PATH_IMAGE086AA
为时刻的三维坐标，
Figure DEST_PATH_IMAGE204
为
Figure 100002_DEST_PATH_IMAGE036AAAAAAAAAAAAA
时刻向
Figure DEST_PATH_IMAGE206
时刻变化的速度；
Figure DEST_PATH_IMAGE038A
为在
Figure DEST_PATH_IMAGE028AAAAAAAAAAAAA
时刻的状态向量的最优估值的方差阵；
Figure DEST_PATH_IMAGE022A
为
Figure DEST_PATH_IMAGE024A
的转置；将加速度作为随机扰动项，
Figure DEST_PATH_IMAGE048A
为
Figure DEST_PATH_IMAGE028AAAAAAAAAAAAAA
时刻的动态噪声的方差阵。
第三步、获取
Figure 100002_DEST_PATH_IMAGE028AAAAAAAAAAAAAAA
时刻至
Figure 100002_DEST_PATH_IMAGE036AAAAAAAAAAAAAA
时刻的实时监测数据
Figure DEST_PATH_IMAGE050A
；结合第二步获得的状态向量的一步预测值
Figure DEST_PATH_IMAGE010AAA
，通过表达式8）得到
Figure 100002_DEST_PATH_IMAGE028AAAAAAAAAAAAAAAA
时刻至
Figure 100002_DEST_PATH_IMAGE036AAAAAAAAAAAAAAA
时刻该时间段数据的标准方差
Figure DEST_PATH_IMAGE052A
：
Figure DEST_PATH_IMAGE056A
8）；
通过表达式9）计算卡尔曼增益
Figure DEST_PATH_IMAGE054A
：
Figure DEST_PATH_IMAGE058A
9）；
其中：
Figure DEST_PATH_IMAGE062A
为B的转置；
Figure DEST_PATH_IMAGE064A
为
Figure 100002_DEST_PATH_IMAGE036AAAAAAAAAAAAAAAA
时刻的观测噪声方差阵；
第四步、通过表达式10）和表达式11）分别计算
Figure DEST_PATH_IMAGE098A
时刻状态向量的最优估值
Figure DEST_PATH_IMAGE066A
和最优估值的方差阵
Figure DEST_PATH_IMAGE068A
：
Figure DEST_PATH_IMAGE070A
10）；
Figure DEST_PATH_IMAGE072A
11）；
其中：
Figure DEST_PATH_IMAGE074A
为4×4单位矩阵；
构造滑动窗口残差向量并动态调整滑动窗口的大小，实时更新观测噪声方差阵：
实时更新观测噪声方差阵具体包括以下步骤：
步骤①、缓存最近的m个北斗结算结果构成观测值滑动窗口，
Figure DEST_PATH_IMAGE092A
且m为整数，此处m的取值与接收机采样间隔有关，测试中采样间隔为15s时，
Figure 100002_DEST_PATH_IMAGE092AA
比较合适，即缓存数据的时间长度最小1h，最大12h，实际运算时m初始值取240，每增加一个监测结果，取m+1（解算结果没有累积够240个时不进行滤波处理）；当m达到2880个时，则暂时保持稳定不变；构建滑动窗口残差向量
Figure DEST_PATH_IMAGE094A
；
步骤②、通过表达式12）计算观测数据的方差作为
Figure 100002_DEST_PATH_IMAGE028AAAAAAAAAAAAAAAAA
时刻至
Figure DEST_PATH_IMAGE036AAAAAAAAAAAAAAAAA
时刻时间段数据的标准方差
Figure DEST_PATH_IMAGE104A
：
Figure DEST_PATH_IMAGE102A
12）；
其中：
Figure 100002_DEST_PATH_IMAGE104AA
为m个实时数据的方差，
Figure DEST_PATH_IMAGE106A
；
步骤③、根据2倍标准方差计算观测噪声方差阵系数
Figure DEST_PATH_IMAGE108A
；
步骤④、调整观测噪声缩小因子
Figure DEST_PATH_IMAGE078AAA
的值：若
Figure DEST_PATH_IMAGE110A
（表明系统稳定），则取
Figure DEST_PATH_IMAGE112A
，即不干涉
Figure 100002_DEST_PATH_IMAGE064AA
的值，也就是
Figure DEST_PATH_IMAGE064AAA
不做调整；若
Figure DEST_PATH_IMAGE114A
（表明模型发散），则取
Figure DEST_PATH_IMAGE116A
，此时可能是由粗差或监测点真实位移造成：
如果是由粗差造成，采用表达式13）实时更新观测噪声方差阵：
Figure DEST_PATH_IMAGE208
13）；
其中：
Figure DEST_PATH_IMAGE210
为观测噪声方差阵系数；
由于
Figure DEST_PATH_IMAGE210A
>1，所以
Figure DEST_PATH_IMAGE064AAAA
的值会增大，根据表达式9）可知
Figure DEST_PATH_IMAGE212
减小，因此，
Figure DEST_PATH_IMAGE214
会减小，相当于当前解算结果的权重，粗差并不会反映到滤波结果，下一个历元
Figure DEST_PATH_IMAGE216
则按照系统稳定处理。
如果是监测点位移引起的模型发散，随着监测结果的增多
Figure DEST_PATH_IMAGE218
，因此
Figure 100002_DEST_PATH_IMAGE114AA
，由于
Figure DEST_PATH_IMAGE114AAA
时，
Figure DEST_PATH_IMAGE220
，所以
Figure 100002_DEST_PATH_IMAGE078AAAA
值逐渐减小，
Figure DEST_PATH_IMAGE064AAAAA
逐渐减小，根据表达式9）可知
Figure DEST_PATH_IMAGE212A
增大，因此
Figure DEST_PATH_IMAGE222
会减大，即相当于当前解算结果的权重会增加，达到快速反应真实位移的效果。
标准方差的滑动窗口m根据
Figure 100002_DEST_PATH_IMAGE078AAAAA
的值动态调整：m的初始值为240；若
Figure DEST_PATH_IMAGE118A
，则增加一个解算结果，取
Figure DEST_PATH_IMAGE120A
，当m值达到最大值时，则不增加；若
Figure DEST_PATH_IMAGE122A
，则取
Figure DEST_PATH_IMAGE124A
，
Figure DEST_PATH_IMAGE126A
代表
Figure DEST_PATH_IMAGE128A
取整（即m的值为m的最大取值
Figure DEST_PATH_IMAGE224
与
Figure 100002_DEST_PATH_IMAGE078AAAAAA
的乘积取整后的整数）；返回步骤①。
第五步、计算位移量并输出，具体是：监测点
Figure DEST_PATH_IMAGE226
时刻的位移量
Figure DEST_PATH_IMAGE228
，其中：
Figure DEST_PATH_IMAGE230
通过
Figure DEST_PATH_IMAGE036AAAAAAAAAAAAAAAAAA
时刻状态向量的最优估值
Figure DEST_PATH_IMAGE088A
获得，
Figure DEST_PATH_IMAGE232
；取k= k+1，返回第二步。
利用设置的监测点进行实时解算实验，数据采样间隔为15s，实验中分别采用普通卡尔曼滤波和本实施例中新型卡尔曼滤波分别对实时解算数据进行解算，为了更好的模拟滤波对粗差和真正移动的识别情况，随机在原始结果2379历元和23679历元处加入30mm粗差，另外在53000历元后加入10mm的真实位移。
本实施例处理各个历元的实验数据得到新型卡尔曼滤波结果详见图3，与现有技术（普通卡尔曼滤波处理结果示意图如图2所示）比较可知：两种卡尔曼滤波算法对监测结果中的随机误差都有较好的剔除作用，滤波后的数据要平滑很多。图2普通卡尔曼滤波在两处粗差处的滤波结果呈现出较大跳动，并且后边的2000个历元数据也收到了影响；在53000历元处的真实位移并不能够较好反映，延迟大约4000个历元。图3的新型卡尔曼滤波在两处粗差处都没有受粗差影响，另外在53000历元处立即反映出真实的位移，延时时间几乎可以忽略。实例的测试说明新型卡尔曼滤波算法能区分粗差和真实位移，得到最优的滤波效果。
本实施例能用于实时变形监测，卡尔曼滤波满足动态系统实时处理数据的要求，能够对动态系统进行实时数据处理，本实施例在普通的卡尔曼滤波算法中加入观测噪声方差阵系数、观测噪声缩小因子以及标准方差的动态计算，用于北斗实时解算，动态调整观测噪声的方差值，防止滤波发散，对于实时数据中的粗差有抑制作用，能够过滤掉高频的随机噪声误差，另外也能够及时反映真实的位移情况，有效弥补普通卡尔曼滤波处理方法的缺点，实现更好地剔除数据粗差、识别监测结果中的真实位移以及平滑监测结果的效果，具有重要的应用价值。
以上所述仅为本发明的优选实施例而已，并不用于限制本发明，对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (3)

1.一种基于新型卡尔曼滤波的北斗变形监测实时处理方法，其特征在于，包括以下步骤：
步骤一、获取N天的监测结果作为初始数据集，1≤N≤3；获取初始状态向量的最优估值
Figure DEST_PATH_IMAGE002AAA
、初始状态向量的最优估值的方差阵
Figure DEST_PATH_IMAGE004AAAA
、初始动态噪声方差阵
Figure DEST_PATH_IMAGE006AAA
、初始观测噪声方差阵
Figure DEST_PATH_IMAGE008AAA
；
步骤二、基于卡尔曼滤波，通过表达式6）和表达式7）分别计算得到状态向量的一步预测值
Figure DEST_PATH_IMAGE010AAAA
和状态向量的一步预测方差
Figure DEST_PATH_IMAGE012AA
；
Figure DEST_PATH_IMAGE014AA
6）；
Figure DEST_PATH_IMAGE016AA
7）；
其中：k为历元数，k为大于等于1且小于等于M的自然数；A为转移矩阵，
Figure DEST_PATH_IMAGE018AA
，
Figure DEST_PATH_IMAGE020AA
为历元间的时间间隔，
Figure DEST_PATH_IMAGE022AA
为
Figure DEST_PATH_IMAGE024AA
的转置；
Figure DEST_PATH_IMAGE026AA
为状态向量在
Figure DEST_PATH_IMAGE028AAAAAAAAAAAAAAA
时刻的最优估值，
Figure DEST_PATH_IMAGE028AAAAAAAAAAAAAAAA
为第k-1历元的时刻，
Figure DEST_PATH_IMAGE030AA
，
Figure DEST_PATH_IMAGE032AA
为
Figure DEST_PATH_IMAGE028AAAAAAAAAAAAAAAAA
时刻的三维坐标，
Figure DEST_PATH_IMAGE034AAAAAA
为
Figure DEST_PATH_IMAGE028AAAAAAAAAAAAAAAAAA
时刻向
Figure DEST_PATH_IMAGE036AAAAAAAAA
时刻变化的速度；
Figure DEST_PATH_IMAGE038AAAAAAAA
为在
Figure DEST_PATH_IMAGE028AAAAAAAAAAAAAAAAAAA
时刻的状态向量的最优估值的方差阵；
Figure DEST_PATH_IMAGE040AAA
为噪声矩阵，
Figure DEST_PATH_IMAGE042AAA
，
Figure DEST_PATH_IMAGE044AA
为
Figure DEST_PATH_IMAGE040AAAA
的转置；
Figure DEST_PATH_IMAGE046AA
为
Figure DEST_PATH_IMAGE028AAAAAAAAAAAAAAAAAAAA
时刻的动态噪声，
Figure DEST_PATH_IMAGE048AA
为
Figure DEST_PATH_IMAGE028AAAAAAAAAAAAAAAAAAAAA
时刻的动态噪声的方差阵；
步骤三、获取
Figure DEST_PATH_IMAGE028AAAAAAAAAAAAAAAAAAAAAA
时刻至
Figure DEST_PATH_IMAGE036AAAAAAAAAA
时刻的实时监测数据
Figure DEST_PATH_IMAGE050AA
；结合步骤二获得的状态向量的一步预测值
Figure DEST_PATH_IMAGE010AAAAA
，通过表达式8）得到
Figure DEST_PATH_IMAGE028AAAAAAAAAAAAAAAAAAAAAAA
时刻至
Figure DEST_PATH_IMAGE036AAAAAAAAAAA
时刻该时间段数据的标准方差
Figure DEST_PATH_IMAGE052AA
且通过表达式9）计算卡尔曼增益
Figure DEST_PATH_IMAGE054AA
：
Figure DEST_PATH_IMAGE056AA
8）；
Figure DEST_PATH_IMAGE058AA
9）；
其中：B为观测系数矩阵，
Figure DEST_PATH_IMAGE060AA
，
Figure DEST_PATH_IMAGE062AA
为B的转置；
Figure DEST_PATH_IMAGE064AA
为
Figure DEST_PATH_IMAGE036AAAAAAAAAAAA
时刻的观测噪声方差阵；
步骤四、通过表达式10）和表达式11）分别计算
Figure DEST_PATH_IMAGE036AAAAAAAAAAAAA
时刻状态向量的最优估值
Figure DEST_PATH_IMAGE066AA
和最优估值的方差阵
Figure DEST_PATH_IMAGE068AAA
：
Figure DEST_PATH_IMAGE070AA
10）；
Figure DEST_PATH_IMAGE072AA
11）；
其中：
Figure DEST_PATH_IMAGE074AA
为4×4的单位矩阵；
构造滑动窗口残差向量并动态调整滑动窗口的大小，采用表达式13）实时更新观测噪声方差阵：
Figure DEST_PATH_IMAGE076AA
13）；
其中：
Figure DEST_PATH_IMAGE078AAAA
为观测噪声缩小因子，
Figure DEST_PATH_IMAGE080AAAA
；
Figure DEST_PATH_IMAGE082AA
为观测噪声方差阵系数；
步骤五、计算位移量并输出，具体是：监测点
Figure DEST_PATH_IMAGE036AAAAAAAAAAAAAA
时刻的位移量
Figure DEST_PATH_IMAGE084AA
，其中：
Figure DEST_PATH_IMAGE086AA
通过
Figure DEST_PATH_IMAGE036AAAAAAAAAAAAAAA
时刻状态向量的最优估值
Figure DEST_PATH_IMAGE088AA
获得，
Figure DEST_PATH_IMAGE090AA
；取k= k+1，返回步骤二。
2.根据权利要求1所述的基于新型卡尔曼滤波的北斗变形监测实时处理方法，其特征在于，所述步骤四中构造滑动窗口残差向量并动态调整滑动窗口的大小以及实时更新观测噪声方差阵具体包括以下步骤：
步骤①、缓存最近的m个北斗结算结果构成观测值滑动窗口，
Figure DEST_PATH_IMAGE092AA
且m为整数；构建滑动窗口残差向量
Figure DEST_PATH_IMAGE094AA
；
步骤②、通过表达式12）计算观测数据的方差作为
Figure DEST_PATH_IMAGE028AAAAAAAAAAAAAAAAAAAAAAAA
时刻至
Figure DEST_PATH_IMAGE036AAAAAAAAAAAAAAAA
时刻时间段数据的标准方差
Figure DEST_PATH_IMAGE096AA
：
Figure DEST_PATH_IMAGE098AAA
12）；
其中：
Figure DEST_PATH_IMAGE100AA
为m个实时数据的方差，
Figure DEST_PATH_IMAGE102AA
；
步骤③、根据2倍标准方差计算观测噪声方差阵系数
Figure DEST_PATH_IMAGE104AA
；
步骤④、调整观测噪声缩小因子
Figure DEST_PATH_IMAGE078AAAAA
的值：若
Figure DEST_PATH_IMAGE106AA
，则取
Figure DEST_PATH_IMAGE108AA
；若
Figure DEST_PATH_IMAGE110AA
，则取
Figure DEST_PATH_IMAGE112AA
；采用表达式13）实时更新观测噪声方差阵；
标准方差的滑动窗口m根据
Figure DEST_PATH_IMAGE078AAAAAA
的值动态调整：若
Figure DEST_PATH_IMAGE114AA
，则取
Figure DEST_PATH_IMAGE116AA
；若
Figure DEST_PATH_IMAGE118AA
，则取
Figure DEST_PATH_IMAGE120AA
，
Figure DEST_PATH_IMAGE122AA
代表
Figure DEST_PATH_IMAGE124AA
取整；返回步骤①。
3.根据权利要求1所述的基于新型卡尔曼滤波的北斗变形监测实时处理方法，其特征在于，所述步骤一中：
初始状态向量的最优估值
Figure DEST_PATH_IMAGE126AA
的获得具体是：
Figure DEST_PATH_IMAGE128AAAA
，初始数据集中的监测数据总数为
Figure DEST_PATH_IMAGE130AAAA
，利用最小二乘计算
Figure DEST_PATH_IMAGE130AAAAA
个监测数据的最优估计值作为
Figure DEST_PATH_IMAGE132AA
，
Figure DEST_PATH_IMAGE134AAA
，
Figure DEST_PATH_IMAGE136AAA
为单位矩阵，
Figure DEST_PATH_IMAGE138AAA
为
Figure DEST_PATH_IMAGE136AAAA
的转置，
Figure DEST_PATH_IMAGE140AAA
为观测值向量，P为观测值向量
Figure DEST_PATH_IMAGE140AAAA
的协方差矩阵，数值根据北斗解算软件输出的Ratio值确定，
Figure DEST_PATH_IMAGE142AA
，R为
Figure DEST_PATH_IMAGE130AAAAAA
个监测数据的Ratio值矩阵，
Figure DEST_PATH_IMAGE144AA
为R中所有元素之和；
Figure DEST_PATH_IMAGE146AA
为初始速度，
Figure DEST_PATH_IMAGE148AA
，
Figure DEST_PATH_IMAGE150AA
为第1个监测数据；
初始状态向量的最优估值的方差阵
Figure DEST_PATH_IMAGE004AAAAA
的获得具体是：
Figure DEST_PATH_IMAGE152AA
，
Figure DEST_PATH_IMAGE154AA
，
Figure DEST_PATH_IMAGE156AA
为初始数据集中的第j个监测数据，j为大于0的自然数；
Figure DEST_PATH_IMAGE158AA
；
初始动态噪声方差阵
Figure DEST_PATH_IMAGE160A
的获得：
Figure DEST_PATH_IMAGE162A
，初始观测噪声方差阵
Figure DEST_PATH_IMAGE164A
为北斗解算后的残差方差值。
CN202010733758.1A 2020-07-28 2020-07-28 一种基于新型卡尔曼滤波的北斗变形监测实时处理方法 Pending CN111623703A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010733758.1A CN111623703A (zh) 2020-07-28 2020-07-28 一种基于新型卡尔曼滤波的北斗变形监测实时处理方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010733758.1A CN111623703A (zh) 2020-07-28 2020-07-28 一种基于新型卡尔曼滤波的北斗变形监测实时处理方法

Publications (1)

Publication Number Publication Date
CN111623703A true CN111623703A (zh) 2020-09-04

Family

ID=72273137

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010733758.1A Pending CN111623703A (zh) 2020-07-28 2020-07-28 一种基于新型卡尔曼滤波的北斗变形监测实时处理方法

Country Status (1)

Country Link
CN (1) CN111623703A (zh)

Cited By (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112269192A (zh) * 2020-10-22 2021-01-26 云南航天工程物探检测股份有限公司 一种快速自适应的动态北斗监测实时解算去噪方法
CN112344847A (zh) * 2020-11-20 2021-02-09 中国有色金属长沙勘察设计研究院有限公司 一种地基合成孔径雷达数据降噪方法
CN112381309A (zh) * 2020-11-23 2021-02-19 珠江水利委员会珠江水利科学研究院 水库大坝安全监测预警方法、装置、系统及存储介质
CN113109848A (zh) * 2021-04-14 2021-07-13 长沙学院 一种载体低速旋转条件下的导航定位方法及其应用系统
CN114912551A (zh) * 2022-07-18 2022-08-16 中国铁路设计集团有限公司 面向桥梁变形监测的gnss和加速度计实时融合算法

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101750066A (zh) * 2009-12-31 2010-06-23 中国人民解放军国防科学技术大学 基于卫星定位的sins动基座传递对准方法
CN108317949A (zh) * 2018-02-07 2018-07-24 桂林电子科技大学 一种rtk高精度差分定位形变监测系统及方法
CN108548479A (zh) * 2018-04-16 2018-09-18 武汉大学 基于gnss的桥梁实时监测快速初始化方法

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101750066A (zh) * 2009-12-31 2010-06-23 中国人民解放军国防科学技术大学 基于卫星定位的sins动基座传递对准方法
CN108317949A (zh) * 2018-02-07 2018-07-24 桂林电子科技大学 一种rtk高精度差分定位形变监测系统及方法
CN108548479A (zh) * 2018-04-16 2018-09-18 武汉大学 基于gnss的桥梁实时监测快速初始化方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
雷孟飞 等: "自适应卡尔曼滤波在BDS变形监测数据处理中的应用", 《导航定位学报》 *

Cited By (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112269192A (zh) * 2020-10-22 2021-01-26 云南航天工程物探检测股份有限公司 一种快速自适应的动态北斗监测实时解算去噪方法
CN112269192B (zh) * 2020-10-22 2024-02-02 云南航天工程物探检测股份有限公司 一种快速自适应的动态北斗监测实时解算去噪方法
CN112344847A (zh) * 2020-11-20 2021-02-09 中国有色金属长沙勘察设计研究院有限公司 一种地基合成孔径雷达数据降噪方法
CN112381309A (zh) * 2020-11-23 2021-02-19 珠江水利委员会珠江水利科学研究院 水库大坝安全监测预警方法、装置、系统及存储介质
CN112381309B (zh) * 2020-11-23 2022-04-12 珠江水利委员会珠江水利科学研究院 水库大坝安全监测预警方法、装置、系统及存储介质
CN113109848A (zh) * 2021-04-14 2021-07-13 长沙学院 一种载体低速旋转条件下的导航定位方法及其应用系统
CN113109848B (zh) * 2021-04-14 2023-03-14 长沙学院 一种载体低速旋转条件下的导航定位方法
CN114912551A (zh) * 2022-07-18 2022-08-16 中国铁路设计集团有限公司 面向桥梁变形监测的gnss和加速度计实时融合算法

Similar Documents

Publication Publication Date Title
CN111623703A (zh) 一种基于新型卡尔曼滤波的北斗变形监测实时处理方法
CN106855628B (zh) 一种高动态卫星导航信号的快速捕获和跟踪系统和方法
CN108827322B (zh) 一种多星协同测向定位观测系统优化设计与评估方法
WO2024016369A1 (zh) 面向桥梁变形监测的gnss和加速度计实时融合算法
CN108120452B (zh) Mems陀螺仪动态数据的滤波方法
CN109031261B (zh) 一种时差估计方法及装置
CN110007299A (zh) 一种基于混合坐标伪谱技术的微弱目标检测跟踪方法
CN110231620A (zh) 一种噪声相关系统跟踪滤波方法
KR100951321B1 (ko) 파티클 필터 기반의 음향 센서를 이용한 3차원 공간에서의객체 추적 방법
CN109190647B (zh) 一种有源无源数据融合方法
CN116224320B (zh) 一种极坐标系下处理多普勒量测的雷达目标跟踪方法
CN111999750B (zh) 针对杆臂不准的实时单站周跳探测改进方法
Chen et al. IMM tracking of a 3D maneuvering target with passive TDOA system
CN106153046A (zh) 一种基于自适应Kalman滤波的陀螺随机噪声AR建模方法
CN115856963A (zh) 基于深度神经网络学习的高精度定位算法
CN112880659B (zh) 一种基于信息概率的融合定位方法
Liang et al. Multipath time delay estimation using higher order statistics
CN112241583A (zh) 一种最小化后验距离的传感器路径优化方法
CN109270344B (zh) 脉冲丢失下的相参脉冲信号频率估计方法
CN112747773A (zh) 基于Allan方差和随机多项式提高陀螺仪精度的方法
Yan et al. Optimal multirate filtering with its application in estimation of the current of a transformer
CN117908034B (zh) 基于自适应波束跟踪水下目标高稳健模基doa估计方法
CN116338755A (zh) 噪声积分数据融合定位方法、设备、系统及存储介质
CN113050036B (zh) 基于多谐振点传声器阵列的gis波束形成定位方法
CN106570330A (zh) 一种用于扩展目标跟踪的形态估计性能评估方法

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
RJ01 Rejection of invention patent application after publication
RJ01 Rejection of invention patent application after publication

Application publication date: 20200904