CN115453596A - 变分贝叶斯抗差自适应滤波方法、滤波器、设备及介质 - Google Patents

变分贝叶斯抗差自适应滤波方法、滤波器、设备及介质 Download PDF

Info

Publication number
CN115453596A
CN115453596A CN202211069615.0A CN202211069615A CN115453596A CN 115453596 A CN115453596 A CN 115453596A CN 202211069615 A CN202211069615 A CN 202211069615A CN 115453596 A CN115453596 A CN 115453596A
Authority
CN
China
Prior art keywords
covariance matrix
measurement
noise covariance
vehicle
vector
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
CN202211069615.0A
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.)
China University of Mining and Technology CUMT
Original Assignee
China University of Mining and Technology CUMT
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 China University of Mining and Technology CUMT filed Critical China University of Mining and Technology CUMT
Priority to CN202211069615.0A priority Critical patent/CN115453596A/zh
Publication of CN115453596A publication Critical patent/CN115453596A/zh
Pending legal-status Critical Current

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S19/00—Satellite radio beacon positioning systems; Determining position, velocity or attitude using signals transmitted by such systems
    • G01S19/38—Determining a navigation solution using signals transmitted by a satellite radio beacon positioning system
    • G01S19/39—Determining a navigation solution using signals transmitted by a satellite radio beacon positioning system the satellite radio beacon positioning system transmitting time-stamped messages, e.g. GPS [Global Positioning System], GLONASS [Global Orbiting Navigation Satellite System] or GALILEO
    • G01S19/42—Determining position
    • G01S19/45—Determining position by combining measurements of signals from the satellite radio beacon positioning system with a supplementary measurement
    • G01S19/47—Determining position by combining measurements of signals from the satellite radio beacon positioning system with a supplementary measurement the supplementary measurement being an inertial measurement, e.g. tightly coupled inertial
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S19/00—Satellite radio beacon positioning systems; Determining position, velocity or attitude using signals transmitted by such systems
    • G01S19/01—Satellite radio beacon positioning systems transmitting time-stamped messages, e.g. GPS [Global Positioning System], GLONASS [Global Orbiting Navigation Satellite System] or GALILEO
    • G01S19/13—Receivers
    • G01S19/21—Interference related issues ; Issues related to cross-correlation, spoofing or other methods of denial of service
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S19/00—Satellite radio beacon positioning systems; Determining position, velocity or attitude using signals transmitted by such systems
    • G01S19/01—Satellite radio beacon positioning systems transmitting time-stamped messages, e.g. GPS [Global Positioning System], GLONASS [Global Orbiting Navigation Satellite System] or GALILEO
    • G01S19/13—Receivers
    • G01S19/23—Testing, monitoring, correcting or calibrating of receiver elements
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S19/00—Satellite radio beacon positioning systems; Determining position, velocity or attitude using signals transmitted by such systems
    • G01S19/01—Satellite radio beacon positioning systems transmitting time-stamped messages, e.g. GPS [Global Positioning System], GLONASS [Global Orbiting Navigation Satellite System] or GALILEO
    • G01S19/13—Receivers
    • G01S19/35—Constructional details or hardware or software details of the signal processing chain
    • G01S19/37—Hardware or software details of the signal processing chain
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/10—Noise analysis or noise optimisation

Landscapes

  • Engineering & Computer Science (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Remote Sensing (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • Signal Processing (AREA)
  • Theoretical Computer Science (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Geometry (AREA)
  • General Engineering & Computer Science (AREA)
  • Position Fixing By Use Of Radio Waves (AREA)

Abstract

本发明公开了一种变分贝叶斯抗差自适应滤波方法、滤波器、设备及介质，涉及GNSS/INS组合导航定位及滤波器领域，包括：获取车载全球导航卫星系统GNSS接收机以及车载惯性传感器INS的测量信息，并基于测量信息确定量测向量和状态向量；基于量测向量和状态向量，采用自适应学生T分布方法重构厚尾量测噪声以重构量测噪声协方差阵，基于变分贝叶斯方法估计系统噪声协方差阵；基于上述重构后的协方差阵，确定计算后验PDF的自由分解近似时所需要的参数；根据参数，采用定点迭代方式确定最终定位结果；所述最终定位结果包括最终的量测噪声协方差阵、系统噪声协方差阵和状态估计值。本发明能够提高在不精确噪声干扰下的适应能力。

Description

变分贝叶斯抗差自适应滤波方法、滤波器、设备及介质
技术领域
本发明涉及GNSS/INS组合导航定位及滤波器领域，特别是涉及一种变分贝叶斯抗差自适应滤波方法、滤波器、电子设备及计算机可读存储介质。
背景技术
在GNSS/INS组合导航系统定位过程中，滤波器的性能将对最终定位精度取决定作用，随着滤波算法的不断发展，不仅促使多系统组合精密单点定位技术在各行业都有着广泛的应用，也更保证了GNSS/INS组合导航系统长时序复杂环境中的导航精度。目前在GNSS/INS组合导航系统定位技术领域中通常采用卡尔曼滤波器进行估计，卡尔曼滤波是线性高斯状态空间模型的一种最优状态估计，在线性高斯模型条件下算法的优势发挥到极致，然而卡尔曼滤波本身极其依赖于准确的模型与噪声统计。由于高精度的结果需要以精确的随机模型以及观测噪声为前提，在车辆载体复杂的运动状态下，这些噪声往往是未知甚至随时间不断变化。此外，当经过城市区域等高遮挡环境时，由于对卫星信号的遮挡也会导致多路径以及数据缺失，使观测质量大幅下降，这也就导致噪声协方差阵的先验知识并不准确，从而产生大量估计误差，甚至滤波发散。
传统变分贝叶斯自适应滤波器针对不精确噪声有着很好的估计效果，然而在实际应用中由于观测环境变化与卫星数变化导致量测噪声协方差阵无法有效继承，从而导致定位性能下降，另一方面，车辆等经过不同环境受到树木房屋遮挡以及湖面多路径影响都会导致严重的噪声干扰，使观测噪声明显偏离正态分布，且高机动载体下系统噪声往往未知，这些都将严重影响GNSS/INS组合导航系统的定位性能。
发明内容
基于此，本发明的目的是提供一种可应用于GNSS/INS组合导航系统定位中系统噪声与量测噪声不准确时的变分贝叶斯抗差自适应滤波方法、滤波器、电子设备及计算机可读存储介质。
为实现上述目的，本发明提供了如下方案：
第一方面，本发明提供了一种变分贝叶斯抗差自适应滤波方法，包括：
获取车载全球导航卫星系统GNSS接收机以及车载惯性传感器INS的测量信息，并基于所述测量信息确定量测向量和状态向量；
基于所述量测向量和所述状态向量，采用自适应学生T分布方法重构厚尾量测噪声以重构量测噪声协方差阵，基于变分贝叶斯方法估计系统噪声协方差阵；
基于重构后的测噪声协方差阵和估计的系统噪声协方差阵，确定计算后验PDF的自由分解近似时所需要的参数；所述参数为预测协方差阵、辅助随机变量、尺度矩阵、自由度参数和状态向量；
根据所述参数，采用定点迭代方式确定最终定位结果；所述最终定位结果包括最终的量测噪声协方差阵、系统噪声协方差阵和状态估计值。
可选地，所述测量信息至少包括车载GNSS接收机提供的伪距观测和载波相位观测值、以及车载INS提供的车辆运动载体的线加速度信息和角加速度信息；所述基于所述测量信息确定量测向量和状态向量，具体包括：
对车载GNSS接收机与车载INS获得的测量信息进行处理分别获得车载GNSS接收机以及车载INS解算的位置信息与速度信息；
将GNSS接收机解算的位置信息和速度信息，与车载INS解算的位置信息和速度信息作差，并将差值作为量测向量，且在传递进滤波中估计出最终的状态向量。
可选地，所述根据所述参数，采用定点迭代方式确定最终的量测噪声协方差阵、系统噪声协方差阵和状态估计值，具体包括：
确定相邻两次迭代对应的状态估计值的差值是否小于设定阈值或者当前迭代次数是否达到预设最大迭代次数；
当相邻两次迭代对应的状态估计值的差值小于设定阈值或者当前迭代次数达到预设最大迭代次数，则将当前迭代次数对应的量测噪声协方差阵、系统噪声协方差阵和状态估计值确定为最终定位结果；
当相邻两次迭代对应的状态估计值的差值大于或者等于设定阈值且当前迭代次数未达到预设最大迭代次数时，基于当前迭代次数对应的量测噪声协方差阵、系统噪声协方差阵和状态估计值，计算下一迭代次数对应的当前迭代次数对应的量测噪声协方差阵、系统噪声协方差阵和状态估计值，并返回步骤确定相邻两次迭代对应的状态估计值的差值是否小于设定阈值或者当前迭代次数是否达到预设最大迭代次数。
第二方面，本发明提供了一种变分贝叶斯抗差自适应滤波器，包括：
获取模块，用于获取车载全球导航卫星系统GNSS接收机以及车载惯性传感器INS的测量信息，并基于所述测量信息确定量测向量和状态向量；
重构模块，用于基于所述量测向量和所述状态向量，采用自适应学生T分布方法重构厚尾量测噪声以重构量测噪声协方差阵，基于变分贝叶斯方法估计系统噪声协方差阵；
参数确定模块，用于基于重构后的测噪声协方差阵和估计的系统噪声协方差阵，确定计算后验PDF的自由分解近似时所需要的参数；所述参数为预测协方差阵、辅助随机变量、尺度矩阵、自由度参数和状态向量；
定位结果确定模块，用于根据所述参数，采用定点迭代方式确定最终定位结果；所述最终定位结果包括最终的量测噪声协方差阵、系统噪声协方差阵和状态估计值。
可选地，所述测量信息至少包括车载GNSS接收机提供的伪距观测和载波相位观测值、以及车载INS提供的车辆运动载体的线加速度信息和角加速度信息；在所述基于所述测量信息确定量测向量和状态向量方面，所述获取模块具体包括：
对车载GNSS接收机与车载INS获得的测量信息进行处理分别获得车载GNSS接收机以及车载INS解算的位置信息与速度信息；
将GNSS接收机解算的位置信息和速度信息，与车载INS解算的位置信息和速度信息作差，并将差值作为量测向量，且在传递进滤波中估计出最终的状态向量。
可选地，所述定位结果确定模块，具体包括：
条件判断单元，用于确定相邻两次迭代对应的状态估计值的差值是否小于设定阈值或者当前迭代次数是否达到预设最大迭代次数；
定位结果确定单元，用于当相邻两次迭代对应的状态估计值的差值小于设定阈值或者当前迭代次数达到预设最大迭代次数，则将当前迭代次数对应的量测噪声协方差阵、系统噪声协方差阵和状态估计值确定为最终定位结果；
返回单元，用于当相邻两次迭代对应的状态估计值的差值大于或者等于设定阈值且当前迭代次数未达到预设最大迭代次数时，基于当前迭代次数对应的量测噪声协方差阵、系统噪声协方差阵和状态估计值，计算下一迭代次数对应的当前迭代次数对应的量测噪声协方差阵、系统噪声协方差阵和状态估计值，并返回条件判断单元。
第三方面，本发明提供了一种电子设备，包括存储器及处理器，所述存储器用于存储计算机程序，所述处理器运行所述计算机程序以使所述电子设备执行根据第一方面所述的变分贝叶斯抗差自适应滤波方法。
第四方面，本发明提供了一种计算机可读存储介质，其存储有计算机程序，所述计算机程序被处理器执行时实现如第一方面所述的变分贝叶斯抗差自适应滤波方法。
根据本发明提供的具体实施例，本发明公开了以下技术效果：
基于现状，本发明构建了一个新的滤波器，基于GNSS/INS组合导航定位领域中观测噪声粗差所固有的厚尾特性，通过采用自适应学生T分布对厚尾量测噪声进行建模，更新量测噪声协方差阵，同时采用变分贝叶斯方法估计未知的系统噪声协方差阵并进行迭代以达到最优估计，这种新的滤波器具有很优秀的鲁棒性，以此提高在不精确噪声干扰下的适应能力，是一种在复杂环境下保证良好估计性能的有效方法。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。
图1为本发明变分贝叶斯抗差自适应滤波方法的流程图；
图2为本发明一种针对复杂环境的变分贝叶斯抗差自适应滤波方法流程图；
图3为本发明变分贝叶斯抗差自适应滤波器的结构图。
具体实施方式
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
为使本发明的上述目的、特征和优点能够更加明显易懂，下面结合附图和具体实施方式对本发明作进一步详细的说明。
本发明的目的在于提供一种可应用于GNSS/INS组合导航系统定位中系统噪声与量测噪声不准确时的变分贝叶斯抗差自适应滤波方法及滤波器。
实现本发明目的的技术方案为：获取观测数据后，基于学生T分布修正厚尾量测噪声协方差阵，并采用变分贝叶斯更新未知的系统噪声同时进行定点迭代，以此达到最优估计，此发明可提升滤波的适应性以及精度，从而在复杂环境下获得更好的定位性能。
实施例一
如图1所示，本发明实施例提供的变分贝叶斯抗差自适应滤波方法，具体包括以下步骤：
步骤100：获取车载全球导航卫星系统GNSS接收机以及车载惯性传感器INS的测量信息，并基于所述测量信息确定量测向量和状态向量。
所述测量信息包括车载GNSS接收机提供的伪距观测和载波相位观测值、以及车载INS提供的车辆运动载体的线加速度信息和角加速度信息。其中车载GNSS接收机接收数据时受到的随机干扰称为量测噪声，在不受粗差影响下符合正态分布，然而实际观测环境下极易受到外界环境干扰从而产生粗差；车载INS中加速度计以及陀螺仪等仪器原因产生的随机误差称为系统噪声，系统噪声呈正态分布且不易受外界环境干扰，然而系统噪声与车载仪器相关，这也导致了系统噪声往往是未知且时变的，这都会使获取的测量信息并不准确。
通过对车载GNSS接收机与车载INS获得的测量信息进行处理分别获得车载GNSS接收机以及车载INS解算的位置信息与速度信息，并将两者间的差值作为量测向量zk，传递进滤波中从而估计出最终的状态向量xk，即车辆运动载体的位置、速度与姿态等定位结果。
针对上述系统噪声与量测噪声不准确的问题，需在步骤200中进行重新建模。
步骤200：基于所述量测向量和所述状态向量，采用自适应学生T分布方法重构厚尾量测噪声以重构量测噪声协方差阵，基于变分贝叶斯方法估计系统噪声协方差阵。
首先将厚尾量测噪声建模为学生T分布，即St(zk；Hkxk,Rk,vk)，定义修正后的量测噪声协方差阵为
Figure BDA0003829214800000061
其中St()为学生T后验概率密度(PDF)；zk表示量测向量；Hkxk为均值向量；Rk为尺度矩阵。初始值由位置与速度的验后协方差组成；vk为自由度参数；λk为辅助随机变量。然而此后验PDF无法直接求得，但学生T分布可以表示为无限多高斯分布的混合，因此可采用下式进行重新表示：
Figure BDA0003829214800000062
其中，N()为高斯后验PDF，G()为伽马后验PDF。
其次，由于实际工程应用中，系统噪声与量测噪声不同，系统噪声不易被环境因素所影响且定义修正后的系统噪声协方差阵与真实值并非简单的线性关系，此时构建为学生T分布时会由于先验分布与真实的噪声分布存在明显差异导致定位精度下降，因此这里依旧构建为高斯分布，并采用变分贝叶斯VB方法去联合估计未知的系统噪声，即N(xk；xk|k-1,Pk|k-1)，其中xk|k-1为预测状态向量；Pk|k-1为预测协方差阵，由INS机械编排组成。
步骤300：基于重构后的测噪声协方差阵和估计的系统噪声协方差阵，确定计算后验PDF的自由分解近似时所需要的参数。
为联合求取量测噪声与系统噪声的后验PDF需要对状态向量xk、预测协方差阵Pk|k-1、辅助随机变量λk、尺度矩阵Rk以及自由度参数vk进行联合估计，即p(Θk|z1:k)，其中Θk{xk,Pk|k-1,λk,Rk,vk}。然而上述联合后验PDF并无精确的解析解，需要利用变分贝叶(VB)方法找出联合后验PDF的自由分解近似，即p(Θk∣z1:k)≈q(xk)q(Pk|k-1)q(λk)q(Rk)q(vk)。
其中q()表示近似后验PDF，最优解可表示为
Figure BDA0003829214800000071
其中E[]代表期望；φ代表Θ中任意一元素；c为常数。由于q(xk)、q(Pk|k-1)、q(λk)、q(Rk)与q(vk)之间相互耦合，在步骤400中需求取各元素近似后验PDF的同时需要进行定点迭代，以保证各元素后验PDF均达到最优解，并得到修正后的量测噪声Rk、系统噪声
Figure BDA0003829214800000072
与状态估计值
Figure BDA0003829214800000073
步骤400：根据所述参数，采用定点迭代方式确定最终定位结果；所述最终定位结果包括最终的量测噪声协方差阵、系统噪声协方差阵和状态估计值。具体如下：
步骤A：确定相邻两次迭代对应的状态估计值的差值是否小于设定阈值或者当前迭代次数是否达到预设最大迭代次数；
步骤B：当相邻两次迭代对应的状态估计值的差值小于设定阈值或者当前迭代次数达到预设最大迭代次数，则将当前迭代次数对应的量测噪声协方差阵、系统噪声协方差阵和状态估计值确定为最终定位结果；
步骤C：当相邻两次迭代对应的状态估计值的差值大于或者等于设定阈值且当前迭代次数未达到预设最大迭代次数时，基于当前迭代次数对应的量测噪声协方差阵、系统噪声协方差阵和状态估计值，计算下一迭代次数对应的当前迭代次数对应的量测噪声协方差阵、系统噪声协方差阵和状态估计值，并返回步骤确定相邻两次迭代对应的状态估计值的差值是否小于设定阈值或者当前迭代次数是否达到预设最大迭代次数。
1)预测协方差阵
令φ＝Pk|k-1，可得出Pk|k-1于第i+1次时的近似后验PDFq(i+1)(Pk|k-1)，由于Pk|k-1是未知的，因此需要选择共轭先验分布，因为共轭性可以保证后验分布与先验分布有相同的函数形式，这里采用逆Wishart分布作为Pk|k-1的共轭先验分布，定义其自由度参数与逆尺度矩阵分别为
Figure BDA0003829214800000081
和
Figure BDA0003829214800000082
即：
Figure BDA0003829214800000083
自由度参数以及逆尺度矩阵的表达式分别如下：
Figure BDA0003829214800000084
其中uk|k-1为自由度参数初始迭代值，这里设为n+τ+1，其中n为状态估计向量维数，τ为调节因子设为4；Uk|k-1为尺度矩阵初始迭代值，设为τPk|k-1；Ak为中间变量；
Figure BDA0003829214800000085
为估计误差协方差阵；xk|k为次优估计值。将得到的自由度参数
Figure BDA0003829214800000086
与逆尺度矩阵
Figure BDA0003829214800000087
求解后修正系统噪声协方差阵
Figure BDA0003829214800000088
即
Figure BDA0003829214800000089
用于第5步中状态向量的计算。
2)辅助随机变量
令φ＝λk，代入公式中，q(λk)可被计算为
Figure BDA00038292148000000810
其中
Figure BDA00038292148000000811
为形状参数；
Figure BDA00038292148000000812
为速率参数；m为量测噪声向量维数；tr()为矩阵的迹；E(i)[vk]为自由度参数的期望值，初次迭代时设为5；Bk为中间变量。可以表示为如下：
Figure BDA0003829214800000091
将得到的形状参数与速率参数计算λk期望值
Figure BDA0003829214800000092
参与第3步尺度矩阵的计算。
3)尺度矩阵
令φ＝Rk，代入公式中，q(Rk)可被计算为
Figure BDA0003829214800000093
其中
Figure BDA0003829214800000094
为自由度参数；
Figure BDA0003829214800000095
为逆尺度矩阵。将得到的自由度参数与逆尺度矩阵用于更新尺度矩阵Rk逆矩阵的期望值
Figure BDA0003829214800000096
结合公式即可得出本次迭代后修正的量测噪声协方差阵
Figure BDA0003829214800000097
用于第5步中状态向量的计算。
4)自由度参数
令φ＝vk，代入公式中，q(vk)可被计算为
Figure BDA0003829214800000098
其中
Figure BDA0003829214800000099
为形状参数；
Figure BDA00038292148000000910
为速率参数；辅助随机变量对数的期望
Figure BDA00038292148000000911
将得到的形状参数与速率参数用于更新自由度参数vk的期望
Figure BDA00038292148000000912
并用于该历元下一次辅助随机变量的迭代当中。
5)状态向量
令φ＝xk，代入公式中，将第1步与第3步得到的修正后的系统噪声协方差阵
Figure BDA00038292148000000913
与量测噪声协方差阵Rk代入，q(xk)可被计算为
Figure BDA0003829214800000101
其中Kk为卡尔曼增益矩阵；Hk为转换矩阵。将得到的次优估计
Figure BDA0003829214800000102
与估计误差协方差阵
Figure BDA0003829214800000103
返回第1步重新修正预测协方差阵。
当两次迭代中步骤400第5步状态估计值xk|k之差小于1e-16时或达到预设最大迭代次数时，输出当前迭代次数时步骤400第3步量测噪声协方差阵Rk、步骤400第5步系统噪声协方差阵
Figure BDA0003829214800000104
以及精确的状态估计值
Figure BDA0003829214800000105
从而获得较高精度的定位结果。
有益效果
本发明提出一种变分贝叶斯自适应滤波算法，针对复杂环境下系统噪声时变未知以及量测噪声呈现厚尾分布的问题，通过基于学生T分布重新构造量测噪声协方差阵，并采用变分贝叶斯估计未知系统噪声的方式相结合，通过定点迭代机制进行优化，以此提高目标状态的跟踪精度，于复杂环境下可有效提高滤波的鲁棒性。本发明相比于现有的滤波算法，构造了自适应学生T分布来修正量测噪声协方差阵，在变分贝叶斯框架下修正未知的系统噪声协方差阵并定点迭代优化估计结果，得到一种新的自适应滤波算法，对于滤波器领域具有一定的实际工程意义。
实施例二
图2是本发明一种针对复杂环境的变分贝叶斯抗差自适应滤波方法流程图，通过构造学生T分布修正量测噪声矩阵并采用变分贝叶斯估计未知系统噪声进行迭代优化从而达到目标效果，获取观测信息后，各部分具体实施细节如下：
1.进行更新并初始化
进行时间更新：
Figure BDA0003829214800000111
进行系统噪声协方差阵更新：为获取Pk∣k-1的先验信息，首先将Pk∣k-1的均值选择为一步预测误差协方差阵Pk∣k-1，构造
Figure BDA0003829214800000112
因此自由度参数uk|k-1＝n+τ+1，τ为调节参数，一般为4，n为状态估计向量维数；逆尺度矩阵初始值为Uk|k-1＝τPk|k-1。
进行量测噪声协方差阵更新：构造
Figure BDA0003829214800000113
计算Tk，其中自由度参数tk＝n+ρ+1，ρ为调节参数，一般为4，m为量测向量的维数；逆尺度矩阵初始值为Tk＝ρRk；Rk为名义量测噪声矩阵。定义修正后的量测噪声协方差阵为
Figure BDA0003829214800000114
因此需要计算初始期望，分别为：
Figure BDA0003829214800000115
2.循环迭代修正系统噪声协方差阵与量测噪声协方差阵
1)令φ＝Pk|k-1代入到公式
Figure BDA0003829214800000116
中可计算得出：
Figure BDA0003829214800000117
其中
Figure BDA0003829214800000118
将log q(i+1)(Pk|k-1)更新为逆WishartPDF，自由度参数与逆尺度矩阵分别表示为：
Figure BDA0003829214800000119
Figure BDA00038292148000001110
2)令φ＝λk，代入公式
Figure BDA00038292148000001111
中，可以得到：
Figure BDA00038292148000001112
其中vk为尺度矩阵Rk的自由度参数，Bk可表示为：
Figure BDA0003829214800000121
求取期望，可得
Figure BDA0003829214800000122
利用log q(i+1)(λk)，将q(i+1)(λk)更新为伽马PDF，形状参数设为
Figure BDA0003829214800000123
速率参数设为
Figure BDA0003829214800000124
得到：
Figure BDA0003829214800000125
3)令φ＝Rk，代入公式
Figure BDA0003829214800000126
中，可以得到：
Figure BDA0003829214800000127
其中
Figure BDA0003829214800000128
将q(i+1)(Rk)更新为逆WishartPDF，得到自由度参数
Figure BDA0003829214800000129
与逆尺度矩阵
Figure BDA00038292148000001210
Figure BDA00038292148000001211
4)令φ＝vk，代入公式
Figure BDA00038292148000001212
中，可以得到：
Figure BDA00038292148000001213
其中logΓ(0.5vk)可以采用Stirling近似求得，即：
logΓ(0.5vk)≈(0.5vk-0.5)log(0.5vk)-0.5vk (21)
利用上式，将q(vk)更新为伽马PDF，形状参数设为
Figure BDA00038292148000001214
速率参数设为
Figure BDA00038292148000001215
得到
Figure BDA00038292148000001216
6)更新修正后的量测噪声协方差阵：
Figure BDA0003829214800000131
其中
Figure BDA0003829214800000132
E[λk]采用公式得出的
Figure BDA0003829214800000133
公式得出的
Figure BDA0003829214800000134
与
Figure BDA0003829214800000135
用于计算E(i+1)[vk]参与下次迭代，即
Figure BDA0003829214800000136
7)令φ＝xk代入到公式
Figure BDA0003829214800000137
中可计算得出：
Figure BDA0003829214800000138
其中
Figure BDA0003829214800000139
经过i+1次迭代后，修正的一步预测误差协方差阵
Figure BDA00038292148000001310
另外可以从公式中看出xk的后验分布仍是均值为
Figure BDA00038292148000001311
协方差为
Figure BDA00038292148000001312
的高斯分布，其均值与协方差可分别表示如下：
Figure BDA00038292148000001313
3.进行迭代终结判断
当两次迭代中状态估计值
Figure BDA00038292148000001314
之差小于1e-16时或达到预设最大迭代次数时迭代终止，即：
Figure BDA00038292148000001315
或i+1＞imax迭代终止，输出xk|k与Pk|k参与下一历元迭代。
本发明提供的一种可应用于GNSS/INS组合导航系统定位过程中在复杂环境下系统噪声与量测噪声不准确时的变分贝叶斯抗差自适应滤波方法。针对未知的系统噪声以及现有的变分贝叶斯滤波中量测噪声协方差阵无法有效继承且呈现厚尾分布的问题，提出了通过学生T分布修正量测噪声协方差阵，在变分贝叶斯框架下估计未知的系统噪声并通过迭代机制进行联合优化达到最优估计，以此提高滤波的适应性以及精度，从而在复杂环境下获得更好的定位性能。本发明相比于传统变分贝叶斯滤波算法，采取将量测噪声协方差阵建模为学生T分布而非高斯分布，并在变分贝叶斯框架下迭代优化目标状态估计结果，得到一种新的自适应滤波算法，具有一定的实际工程意义
实施例三
如图3所示，本发明提供了一种变分贝叶斯抗差自适应滤波器，包括：
获取模块10，用于获取车载全球导航卫星系统GNSS接收机以及车载惯性传感器INS的测量信息，并基于所述测量信息确定量测向量和状态向量；
重构模块20，用于基于所述量测向量和所述状态向量，采用自适应学生T分布方法重构厚尾量测噪声以重构量测噪声协方差阵，基于变分贝叶斯方法估计系统噪声协方差阵；
参数确定模块30，用于基于重构后的测噪声协方差阵和估计的系统噪声协方差阵，确定计算后验PDF的自由分解近似时所需要的参数；所述参数为预测协方差阵、辅助随机变量、尺度矩阵、自由度参数和状态向量；
定位结果确定模块40，用于根据所述参数，采用定点迭代方式确定最终定位结果；所述最终定位结果包括最终的量测噪声协方差阵、系统噪声协方差阵和状态估计值。
所述测量信息至少包括车载GNSS接收机提供的伪距观测和载波相位观测值、以及车载INS提供的车辆运动载体的线加速度信息和角加速度信息；在所述基于所述测量信息确定量测向量和状态向量方面，所述获取模块具体包括：
对车载GNSS接收机与车载INS获得的测量信息进行处理分别获得车载GNSS接收机以及车载INS解算的位置信息与速度信息；
将GNSS接收机解算的位置信息和速度信息，与车载INS解算的位置信息和速度信息作差，并将差值作为量测向量，且在传递进滤波中估计出最终的状态向量。
所述定位结果确定模块40，具体包括：
条件判断单元，用于确定相邻两次迭代对应的状态估计值的差值是否小于设定阈值或者当前迭代次数是否达到预设最大迭代次数；
定位结果确定单元，用于当相邻两次迭代对应的状态估计值的差值小于设定阈值或者当前迭代次数达到预设最大迭代次数，则将当前迭代次数对应的量测噪声协方差阵、系统噪声协方差阵和状态估计值确定为最终定位结果；
返回单元，用于当相邻两次迭代对应的状态估计值的差值大于或者等于设定阈值且当前迭代次数未达到预设最大迭代次数时，基于当前迭代次数对应的量测噪声协方差阵、系统噪声协方差阵和状态估计值，计算下一迭代次数对应的当前迭代次数对应的量测噪声协方差阵、系统噪声协方差阵和状态估计值，并返回条件判断单元。
实施例四
本发明实施例提供一种电子设备包括存储器及处理器，该存储器用于存储计算机程序，该处理器运行计算机程序以使电子设备执行实施例一的变分贝叶斯抗差自适应滤波方法。
可选地，上述电子设备可以是服务器。
另外，本发明实施例还提供一种计算机可读存储介质，其存储有计算机程序，该计算机程序被处理器执行时实现实施例一的变分贝叶斯抗差自适应滤波方法。
本说明书中各个实施例采用递进的方式描述，每个实施例重点说明的都是与其他实施例的不同之处，各个实施例之间相同相似部分互相参见即可。对于实施例公开的系统而言，由于其与实施例公开的方法相对应，所以描述的比较简单，相关之处参见方法部分说明即可。
本文中应用了具体个例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处。综上所述，本说明书内容不应理解为对本发明的限制。

Claims (8)

1.一种变分贝叶斯抗差自适应滤波方法，其特征在于，包括：
获取车载全球导航卫星系统GNSS接收机以及车载惯性传感器INS的测量信息，并基于所述测量信息确定量测向量和状态向量；
基于所述量测向量和所述状态向量，采用自适应学生T分布方法重构厚尾量测噪声以重构量测噪声协方差阵，基于变分贝叶斯方法估计系统噪声协方差阵；
基于重构后的测噪声协方差阵和估计的系统噪声协方差阵，确定计算后验PDF的自由分解近似时所需要的参数；所述参数为预测协方差阵、辅助随机变量、尺度矩阵、自由度参数和状态向量；
根据所述参数，采用定点迭代方式确定最终定位结果；所述最终定位结果包括最终的量测噪声协方差阵、系统噪声协方差阵和状态估计值。
2.根据权利要求1所述的一种变分贝叶斯抗差自适应滤波方法，其特征在于，所述测量信息至少包括车载GNSS接收机提供的伪距观测和载波相位观测值、以及车载INS提供的车辆运动载体的线加速度信息和角加速度信息；所述基于所述测量信息确定量测向量和状态向量，具体包括：
对车载GNSS接收机与车载INS获得的测量信息进行处理分别获得车载GNSS接收机以及车载INS解算的位置信息与速度信息；
将GNSS接收机解算的位置信息和速度信息，与车载INS解算的位置信息和速度信息作差，并将差值作为量测向量，且在传递进滤波中估计出最终的状态向量。
3.根据权利要求1所述的一种变分贝叶斯抗差自适应滤波方法，其特征在于，所述根据所述参数，采用定点迭代方式确定最终的量测噪声协方差阵、系统噪声协方差阵和状态估计值，具体包括：
确定相邻两次迭代对应的状态估计值的差值是否小于设定阈值或者当前迭代次数是否达到预设最大迭代次数；
当相邻两次迭代对应的状态估计值的差值小于设定阈值或者当前迭代次数达到预设最大迭代次数，则将当前迭代次数对应的量测噪声协方差阵、系统噪声协方差阵和状态估计值确定为最终定位结果；
当相邻两次迭代对应的状态估计值的差值大于或者等于设定阈值且当前迭代次数未达到预设最大迭代次数时，基于当前迭代次数对应的量测噪声协方差阵、系统噪声协方差阵和状态估计值，计算下一迭代次数对应的当前迭代次数对应的量测噪声协方差阵、系统噪声协方差阵和状态估计值，并返回步骤确定相邻两次迭代对应的状态估计值的差值是否小于设定阈值或者当前迭代次数是否达到预设最大迭代次数。
4.一种变分贝叶斯抗差自适应滤波器，其特征在于，包括：
获取模块，用于获取车载全球导航卫星系统GNSS接收机以及车载惯性传感器INS的测量信息，并基于所述测量信息确定量测向量和状态向量；
重构模块，用于基于所述量测向量和所述状态向量，采用自适应学生T分布方法重构厚尾量测噪声以重构量测噪声协方差阵，基于变分贝叶斯方法估计系统噪声协方差阵；
参数确定模块，用于基于重构后的测噪声协方差阵和估计的系统噪声协方差阵，确定计算后验PDF的自由分解近似时所需要的参数；所述参数为预测协方差阵、辅助随机变量、尺度矩阵、自由度参数和状态向量；
定位结果确定模块，用于根据所述参数，采用定点迭代方式确定最终定位结果；所述最终定位结果包括最终的量测噪声协方差阵、系统噪声协方差阵和状态估计值。
5.根据权利要求4所述的一种变分贝叶斯抗差自适应滤波器，其特征在于，所述测量信息至少包括车载GNSS接收机提供的伪距观测和载波相位观测值、以及车载INS提供的车辆运动载体的线加速度信息和角加速度信息；在所述基于所述测量信息确定量测向量和状态向量方面，所述获取模块具体包括：
对车载GNSS接收机与车载INS获得的测量信息进行处理分别获得车载GNSS接收机以及车载INS解算的位置信息与速度信息；
将GNSS接收机解算的位置信息和速度信息，与车载INS解算的位置信息和速度信息作差，并将差值作为量测向量，且在传递进滤波中估计出最终的状态向量。
6.根据权利要求4所述的一种变分贝叶斯抗差自适应滤波器，其特征在于，所述定位结果确定模块，具体包括：
条件判断单元，用于确定相邻两次迭代对应的状态估计值的差值是否小于设定阈值或者当前迭代次数是否达到预设最大迭代次数；
定位结果确定单元，用于当相邻两次迭代对应的状态估计值的差值小于设定阈值或者当前迭代次数达到预设最大迭代次数，则将当前迭代次数对应的量测噪声协方差阵、系统噪声协方差阵和状态估计值确定为最终定位结果；
返回单元，用于当相邻两次迭代对应的状态估计值的差值大于或者等于设定阈值且当前迭代次数未达到预设最大迭代次数时，基于当前迭代次数对应的量测噪声协方差阵、系统噪声协方差阵和状态估计值，计算下一迭代次数对应的当前迭代次数对应的量测噪声协方差阵、系统噪声协方差阵和状态估计值，并返回条件判断单元。
7.一种电子设备，其特征在于，包括存储器及处理器，所述存储器用于存储计算机程序，所述处理器运行所述计算机程序以使所述电子设备执行根据权利要求1至3中任一项所述的变分贝叶斯抗差自适应滤波方法。
8.一种计算机可读存储介质，其特征在于，其存储有计算机程序，所述计算机程序被处理器执行时实现如权利要求1至3中任一项所述的变分贝叶斯抗差自适应滤波方法。
CN202211069615.0A 2022-09-02 2022-09-02 变分贝叶斯抗差自适应滤波方法、滤波器、设备及介质 Pending CN115453596A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202211069615.0A CN115453596A (zh) 2022-09-02 2022-09-02 变分贝叶斯抗差自适应滤波方法、滤波器、设备及介质

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202211069615.0A CN115453596A (zh) 2022-09-02 2022-09-02 变分贝叶斯抗差自适应滤波方法、滤波器、设备及介质

Publications (1)

Publication Number Publication Date
CN115453596A true CN115453596A (zh) 2022-12-09

Family

ID=84300478

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202211069615.0A Pending CN115453596A (zh) 2022-09-02 2022-09-02 变分贝叶斯抗差自适应滤波方法、滤波器、设备及介质

Country Status (1)

Country Link
CN (1) CN115453596A (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116338573A (zh) * 2023-03-30 2023-06-27 中国矿业大学 一种封闭空间考虑噪声误差特性的无人系统定位方法

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116338573A (zh) * 2023-03-30 2023-06-27 中国矿业大学 一种封闭空间考虑噪声误差特性的无人系统定位方法
CN116338573B (zh) * 2023-03-30 2023-12-22 中国矿业大学 一种封闭空间考虑噪声误差特性的无人系统定位方法

Similar Documents

Publication Publication Date Title
CN109738917B (zh) 一种北斗变形监测中的多路径误差削弱方法及装置
Gao et al. Adaptive Kalman filtering with recursive noise estimator for integrated SINS/DVL systems
CN107110650B (zh) 在能观测性方面受约束的导航状态的估计方法
CN109477900A (zh) 全球导航卫星系统接收器中用于模糊度解算的频率间偏差的估算
CN111156987A (zh) 基于残差补偿多速率ckf的惯性/天文组合导航方法
CN110567454B (zh) 一种复杂环境下sins/dvl紧组合导航方法
CN108827310A (zh) 一种船用星敏感器辅助陀螺仪在线标定方法
CN110567455B (zh) 一种求积更新容积卡尔曼滤波的紧组合导航方法
CN116358566B (zh) 一种基于抗差自适应因子的粗差探测组合导航方法
CN115453596A (zh) 变分贝叶斯抗差自适应滤波方法、滤波器、设备及介质
CN116007620A (zh) 一种组合导航滤波方法、系统、电子设备及存储介质
CN114323007A (zh) 一种载体运动状态估计方法及装置
Psiaki Kalman filtering and smoothing to estimate real-valued states and integer constants
Han et al. Precise positioning with machine learning based Kalman filter using GNSS/IMU measurements from android smartphone
CN113671551B (zh) Rtk定位解算方法
CN104407366B (zh) 一种对伪距进行平滑处理的方法
CN113434806A (zh) 一种抗差自适应多模型滤波方法
CN114894222B (zh) Imu-gnss天线的外参数标定方法和相关方法、设备
CN114019954B (zh) 航向安装角标定方法、装置、计算机设备和存储介质
CN116182847A (zh) 车载捷联组合导航方法、装置、电子设备及存储介质
Avzayesh et al. Improved-Performance Vehicle’s State Estimator Under Uncertain Model Dynam
CN114485651B (zh) Fls与ukf相结合的组合导航方法、装置、存储介质及设备
CN110716219A (zh) 一种提高定位解算精度的方法
CN111123323B (zh) 一种提高便携设备定位精度的方法
US11519730B2 (en) Method for determining the position and orientation of a vehicle

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