CN113311456A - 一种基于卡尔曼滤波的qar数据噪声处理方法 - Google Patents

一种基于卡尔曼滤波的qar数据噪声处理方法 Download PDF

Info

Publication number
CN113311456A
CN113311456A CN202110538773.5A CN202110538773A CN113311456A CN 113311456 A CN113311456 A CN 113311456A CN 202110538773 A CN202110538773 A CN 202110538773A CN 113311456 A CN113311456 A CN 113311456A
Authority
CN
China
Prior art keywords
state
time
moment
matrix
covariance matrix
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.)
Granted
Application number
CN202110538773.5A
Other languages
English (en)
Other versions
CN113311456B (zh
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.)
Wuhan University WHU
Original Assignee
Wuhan University WHU
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 Wuhan University WHU filed Critical Wuhan University WHU
Priority to CN202110538773.5A priority Critical patent/CN113311456B/zh
Publication of CN113311456A publication Critical patent/CN113311456A/zh
Application granted granted Critical
Publication of CN113311456B publication Critical patent/CN113311456B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

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/393—Trajectory determination or predictive tracking, e.g. Kalman filtering
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01C—MEASURING DISTANCES, LEVELS OR BEARINGS; SURVEYING; NAVIGATION; GYROSCOPIC INSTRUMENTS; PHOTOGRAMMETRY OR VIDEOGRAMMETRY
    • G01C21/00—Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00
    • G01C21/20—Instruments for performing navigational calculations
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
    • Y02T90/00—Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation

Landscapes

  • Engineering & Computer Science (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Remote Sensing (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Automation & Control Theory (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • Feedback Control In General (AREA)
  • Indication And Recording Devices For Special Purposes And Tariff Metering Devices (AREA)
  • Measurement Of Velocity Or Position Using Acoustic Or Ultrasonic Waves (AREA)

Abstract

本发明提出了一种基于卡尔曼滤波的QAR数据噪声处理方法。本发明定义t时刻状态向量、t时刻状态转移矩阵、t时刻输入传递矩阵；根据t时刻飞机状态向量和t时刻实际状态转移协方差矩阵得到t+1时刻预测状态向量、t+1时刻飞机状态预测误差协方差矩阵、t+1时刻转移状态误差协方差矩阵；根据t+1时刻转移状态误差协方差矩阵得到t+1时刻卡尔曼系数；将t+1时刻预测状态向量与观测状态进行融合得到滤波后的t+1时刻状态。本发明优点在于，考虑到了在飞行过程中关于位置信息存在的测量误差，结合运动物理模型与测量数据有效地减小了这一误差，由于利用了更多信息，减小了单一传感器存在的噪声对结果产生的影响，从而得到了误差更小的位置数据。

Description

一种基于卡尔曼滤波的QAR数据噪声处理方法
技术领域
本发明涉及QAR数据预处理方法领域，具体涉及一种基于卡尔曼滤波的QAR数据噪声处理方法。
背景技术
利用QAR数据可以还原真实的飞行轨迹，但是由于GPS误差以及其他各种各样不可预知的影响因素的存在，QAR数据中总是存在部分噪声，如果不对轨迹数据加以处理，则可能会出现轨迹“打结”、轨迹“回跳”等现象，严重影响了轨迹的展示效果。而为了尽可能真实地还原飞机的飞行轨迹，提升轨迹展示效果，需要对QAR数据进行噪声处理。
卡尔曼滤波是一种利用观测值对预测值进行修正，来得到最优状态估计的算法，其实质是寻找满足均方误差最小条件下的状态估计值。卡尔曼滤波算法在实际应用中仅需要保存上一时刻的状态估计值，而无需保存所有的历史数据，借助系统的状态转移方程，利用上一时刻的状态估计值就可以对当前时刻的状态观测值进行修正，从而得到当前时刻的状态估计值，并以此为一个周期进行迭代，继续对下一时刻的状态值进行估计。卡尔曼递推算法计算的时间复杂度较低，所需要的内存空间较小，非常适用于实时连续系统。
发明内容
本发明的目的在于提供一种基于卡尔曼滤波的QAR数据噪声计算方法，该方法能够利用QAR数据中记录的丰富信息对飞行轨迹点进行有效的纠正，使纠正后的轨迹点更接近于真实的飞行轨迹点，具体步骤如下：
步骤1，定义t时刻状态向量、t时刻状态转移矩阵、t时刻输入传递矩阵；
步骤2，将t时刻飞机状态向量和t时刻实际状态转移协方差矩阵结合状态转移矩阵、输入传递矩阵进行计算，得到t+1时刻预测状态向量、t+1时刻飞机状态预测误差协方差矩阵、t+1时刻转移状态误差协方差矩阵。
步骤3，将t+1时刻状态向量观测值与t+1时刻观测值向量，根据t+1时刻转移状态误差协方差矩阵得到t+1时刻卡尔曼系数；
步骤4，将于步骤2得到的t+1时刻预测状态向量与观测状态通过t+1时刻卡尔曼系数进行融合得到滤波后的t+1时刻状态，并根据状态误差协方差矩阵得到t+1时刻实际状态转移协方差矩阵。
步骤5，重复步骤2到步骤4，如此循环往复直到所有时刻的状态完成处理，得到完整的航迹数据。
作为优选，步骤1所述t时刻状态向量为：
Figure BDA0003070930180000021
其中，Statet为t时刻的状态向量，t用来指代数据的序号，K代表最大序号，与采样次数即数据长度相等，lngt为t时刻飞机所在的经度，latt为时刻t飞机所在的纬度，vt为时刻t飞机相对于地面的速度，当t＝0时，Statet为初始状态向量；
步骤1所述t时刻状态转移矩阵为：
Figure BDA0003070930180000022
Ft表示t时刻的状态转移矩阵，当t＝0时，Ft为初始状态转移矩阵。其中，ωx t表示在t时刻把飞行距离转化成经度需要乘上的系数，ωy表示把飞行距离转化成纬度需要乘上的系数，θt为t时刻飞机的朝向角，Δt为QAR数据的采样间隔，上述ωx t，ωy为：
Figure BDA0003070930180000023
Figure BDA0003070930180000024
其中，Latt表示t时刻飞机所处的纬度，Earth-Radius_Long表示地球长半轴半径，Earth_Radius_Short表示地球短半轴半径。
步骤1所述t时刻输入传递矩阵为：
Figure BDA0003070930180000025
其中，Bt为t时刻的输入传递矩阵，ωx t表示在t时刻把飞行距离转化成经度需要乘上的系数，ωy表示把飞行距离转化成纬度需要乘上的系数，由于在对轨迹进行采样时Δt足够小，根据牛顿运动学定则，此时在做近似匀加速直线运动的状态转移方程为：
Xt+1＝Xt+(Vt*Δt+0.5*at*Δt2)*Sin(θt)*ωx
Yt+1＝Yt+(Vt*Δt+0.5*at*Δt2)*Cos(θt)*ωy
Vt+1＝Vt+at*Δt
其中，at表示t时刻飞机的水平加速度，Vt*Δt+0.5*at*Δt2表示Δt内飞机的水平位移，Vt为t时刻飞机的水平速度，Δt为离散时间的采样间隔；
作为优选，步骤2所述t+1时刻预测状态向量为：
Statet+1 -＝Ft*Statet+Bt*at+xt，0≤t≤K-1
式中，Ft为t时刻的状态转移矩阵，Bt为t时刻的输入传递矩阵，at为t时刻飞机的水平加速度，Ft*Statet表示飞机在两个时刻间的状态转移过程，Bt*at表示水平加速度对飞机状态的影响；xt为t时刻系统状态向量误差，该误差使得飞机动力装置对飞机的控制与记录数据之间存在偏差，符合正态分布xt～N(0，Qt)，Qt为t时刻状态向量误差的协方差矩阵；
步骤2所述t+1时刻飞机状态预测误差协方差矩阵为：
Cov(Statet+1)＝Cov(Ft*Statet)＝FtPtFt T，0≤t≤K-1
步骤2中所述t+1时刻转移状态误差协方差矩阵为：
Pt+1 -＝FtPtFt T+Qt，0≤t≤K-1
式中，Pt+1 -为t+1时刻的转移状态误差协方差矩阵，Pt为t时刻的状态误差协方差矩阵，Ft为t时刻的状态转移矩阵，Qt为t时刻状态向量误差的协方差矩阵；
作为优选，步骤3所述t+1时刻状态向量观测值为：
Zt+1＝H*Statet+1 -+R，0≤t≤K-1
其中，Zt+1为t+1时刻状态向量观测值，H为变换矩阵，R为观测误差协方差矩阵；
步骤3所述t+1时刻观测值向量为：
Figure BDA0003070930180000031
其中，Lngt+1为t+1时刻观测到的飞机所在经度，Latt+1为t+1时刻观测到的飞机所在纬度。
步骤3所述t+1时刻卡尔曼系数为：
Kt+1＝Pt+1 -*HT*(H*Pt+1-*HT+R)-1，0≤t≤K-1
其中，Pt+1 -为t+1时刻转移状态误差协方差矩阵，H为变换矩阵，R为观测误差协方差矩阵。卡尔曼系数的主要作用是衡量根据状态协方差矩阵Pt+1 -和观测状态协方差矩阵R的大小，决定卡尔曼滤波得到的估计值是更偏向于预测值还是更偏向观测值；
作为优选，步骤4所述滤波后的t+1时刻状态为：
Statet+1＝Statet+1 -+Kt+1*(Zt+1-H*Statet+1 -)，0≤t≤K-1
式中，Statet+1 -为t+1时刻的预测状态，Zt+1为t+1时刻的观测状态，H为变化矩阵，Kt+1为t+1时刻的卡尔曼系数。Zt+1-H*Statet+1 -为实际的观测与预测值之间的残差，其与Kt+1的乘积用来修正飞机位置的预测值。
步骤4所述t+1时刻实际状态转移协方差矩阵为：
Pt+1＝(I-Kt*H)*Pt+1 -，0≤t≤K-1
其中，I为单位矩阵，Kt为t时刻的卡尔曼系数，Pt+1 -为t+1时刻转移状态误差协方差矩阵，H为变换矩阵；随着迭代的进行，最佳状态估计值的噪声分布情况需要不断地进行更新；状态的不确定性随着迭代的进行总体上会逐渐变小，但是由于传递噪声的存在，状态的不确定性又会呈现出逐渐变大的趋势；卡尔曼滤波算法就是通过不断地调整从而在这一状态中达到一个平衡；
本发明可以纠正由于测量误差导致的QAR数据中的空间位置误差，过程中结合了预测值与观测值，经过滤波后的空间位置更具有连贯性，符合实际情况。
附图说明
图1：是本发明的基本方法流程示意图。
具体实施方式
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
下面结合图1介绍本发明的具体实施方式如下。
本发明实施例为一种基于卡尔曼滤波的QAR数据噪声处理方法，具体如下：
步骤1，定义t时刻状态向量、t时刻状态转移矩阵、t时刻输入传递矩阵；
步骤1所述t时刻状态向量为：
Figure BDA0003070930180000041
其中，Statet为t时刻的状态向量，t用来指代数据的序号，K代表最大序号，与采样次数即数据长度相等，lngt为t时刻飞机所在的经度，latt为时刻t飞机所在的纬度，vt为时刻t飞机相对于地面的速度，当t＝0时，Statet为初始状态向量，K＝7289；
步骤1所述t时刻状态转移矩阵为：
Figure BDA0003070930180000051
Ft表示t时刻的状态转移矩阵，当t＝0时，Ft为初始状态转移矩阵。其中，ωx t表示在t时刻把飞行距离转化成经度需要乘上的系数，ωy表示把飞行距离转化成纬度需要乘上的系数，θt为t时刻飞机的朝向角，Δt为QAR数据的采样间隔，K＝7289，上述ωx t，ωy为：
Figure BDA0003070930180000052
Figure BDA0003070930180000053
其中，Latt表示t时刻飞机所处的纬度，Earth_Radius-Long表示地球长半轴半径，单位为米，Earth_Radius_Short表示地球短半轴半径，单位为米。
步骤1所述t时刻输入传递矩阵为：
Figure BDA0003070930180000054
其中，Bt为t时刻的输入传递矩阵，ωx t表示在t时刻把飞行距离转化成经度需要乘上的系数，K＝7289，ωy表示把飞行距离转化成纬度需要乘上的系数，由于在对轨迹进行采样时Δt足够小，根据牛顿运动学定则，此时在做近似匀加速直线运动的状态转移方程为：
Xt+1＝Xt+(Vt*Δt+0.5*at*Δt2)*Sin(θt)*ωx
Yt+1＝Yt+(Vt*Δt+0.5*at*Δt2)*Cos(θt)*ωy
Vt+1＝Vt+at*Δt
其中，at表示t时刻飞机的水平加速度，Vt*Δt+0.5*at*Δt2表示Δt内飞机的水平位移，Vt为t时刻飞机的水平速度，Δt为离散时间的采样间隔；
步骤2，将t时刻飞机状态向量和t时刻实际状态转移协方差矩阵结合状态转移矩阵、输入传递矩阵进行计算，得到t+1时刻预测状态向量、t+1时刻飞机状态预测误差协方差矩阵、t+1时刻转移状态误差协方差矩阵。
步骤2所述t+1时刻预测状态向量为：
Statet+1 -＝Ft*Statet+Bt*at+xt，0≤t≤K-1
式中，Ft为t时刻的状态转移矩阵，Bt为t时刻的输入传递矩阵，at为t时刻飞机的水平加速度，Ft*Statet表示飞机在两个时刻间的状态转移过程，Bt*at表示水平加速度对飞机状态的影响；xt为t时刻系统状态向量误差，该误差使得飞机动力装置对飞机的控制与记录数据之间存在偏差，符合正态分布xt～N(0，Qt)，Qt为t时刻状态向量误差的协方差矩阵，K＝7289；
步骤2所述t+1时刻飞机状态预测误差协方差矩阵为：
Cov(Statet+1)＝Cov(Ft*Statet)＝FtPtFt T，0≤t≤K-1
步骤2中所述t+1时刻转移状态误差协方差矩阵为：
Pt+1 -＝FtPtFt T+Qt，0≤t≤K-1
式中，Pt+1 -为t+1时刻的转移状态误差协方差矩阵，Pt为t时刻的状态误差协方差矩阵，Ft为t时刻的状态转移矩阵，Qt为t时刻状态向量误差的协方差矩阵，K＝7289；
步骤3，将t+1时刻状态向量观测值与t+1时刻观测值向量，根据t+1时刻转移状态误差协方差矩阵得到t+1时刻卡尔曼系数；
步骤3所述t+1时刻状态向量观测值为：
Zt+1＝H*Statet+1 -+R，0≤t≤K-1
其中，Zt+1为t+1时刻状态向量观测值，H为变换矩阵，R为观测误差协方差矩阵，K＝7289；
步骤3所述t+1时刻观测值向量为：
Figure BDA0003070930180000061
其中，Lngt+1为t+1时刻观测到的飞机所在经度，Latt+1为t+1时刻观测到的飞机所在纬度，K＝7289。
步骤3所述t+1时刻卡尔曼系数为：
Kt+1＝Pt+1 -*HT*(H*Pt+1 -*HT+R)-1，0≤t≤K-1
其中，Pt+1 -为t+1时刻转移状态误差协方差矩阵，H为变换矩阵，R为观测误差协方差矩阵，K＝7289。卡尔曼系数的主要作用是衡量根据状态协方差矩阵Pt+1 -和观测状态协方差矩阵R的大小，决定卡尔曼滤波得到的估计值是更偏向于预测值还是更偏向观测值；
步骤4，将于步骤2得到的t+1时刻预测状态向量与观测状态通过t+1时刻卡尔曼系数进行融合得到滤波后的t+1时刻状态，并根据状态误差协方差矩阵得到t+1时刻实际状态转移协方差矩阵。
步骤4所述滤波后的t+1时刻状态为：
Statet+1＝Statet+1 -+Kalt+1*(Zt+1-H*Statet+1 -)，0≤t≤K-1
式中，Statet+1 -为t+1时刻的预测状态，Zt+1为t+1时刻的观测状态，H为变化矩阵，Kalt+1为t+1时刻的卡尔曼系数，Zt+1-H*Statet+1 -为实际的观测与预测值之间的残差，其与Kalt+1的乘积用来修正飞机位置的预测值，K＝7289。
步骤4所述t+1时刻实际状态转移协方差矩阵为：
Pt+1＝(I-Kt*H)*Pt+1 -，0≤t≤K-1
其中，I为单位矩阵，Kalt为t时刻的卡尔曼系数，Pt+1 -为t+1时刻转移状态误差协方差矩阵，H为变换矩阵，K＝7289；随着迭代的进行，最佳状态估计值的噪声分布情况需要不断地进行更新；状态的不确定性随着迭代的进行总体上会逐渐变小，但是由于传递噪声的存在，状态的不确定性又会呈现出逐渐变大的趋势；卡尔曼滤波算法就是通过不断地调整从而在这一状态中达到一个平衡；
步骤5，重复步骤2到步骤4，如此循环往复直到所有时刻的状态完成处理，得到完整的航迹数据。
应当理解的是，本说明书未详细阐述的部分均属于现有技术。
应当理解的是，上述针对较佳实施例的描述较为详细，并不能因此而认为是对本发明专利保护范围的限制，本领域的普通技术人员在本发明的启示下，在不脱离本发明权利要求所保护的范围情况下，还可以做出替换或变形，均落入本发明的保护范围之内，本发明的请求保护范围应以所附权利要求为准。

Claims (5)

1.一种基于卡尔曼滤波的QAR数据噪声处理方法，其特征在于，包括以下步骤：
步骤1，定义t时刻状态向量、t时刻状态转移矩阵、t时刻输入传递矩阵；
步骤2，将t时刻飞机状态向量和t时刻实际状态转移协方差矩阵结合状态转移矩阵、输入传递矩阵进行计算，得到t+1时刻预测状态向量、t+1时刻飞机状态预测误差协方差矩阵、t+1时刻转移状态误差协方差矩阵；
步骤3，将t+1时刻状态向量观测值与t+1时刻观测值向量，根据t+1时刻转移状态误差协方差矩阵得到t+1时刻卡尔曼系数；
步骤4，将于步骤2得到的t+1时刻预测状态向量与观测状态通过t+1时刻卡尔曼系数进行融合得到滤波后的t+1时刻状态，并根据状态误差协方差矩阵得到t+1时刻实际状态转移协方差矩阵；
步骤5，重复步骤2到步骤4，如此循环往复直到所有时刻的状态完成处理，得到完整的航迹数据。
2.根据权利要求1所述的基于卡尔曼滤波的QAR数据噪声处理方法，其特征在于，
步骤1所述t时刻状态向量为：
Figure FDA0003070930170000011
其中，Statet为t时刻的状态向量，t用来指代数据的序号，K代表最大序号，与采样次数即数据长度相等，lngt为t时刻飞机所在的经度，latt为时刻t飞机所在的纬度，vt为时刻t飞机相对于地面的速度，当t＝0时，Statet为初始状态向量；
步骤1所述t时刻状态转移矩阵为：
Figure FDA0003070930170000012
Ft表示t时刻的状态转移矩阵，当t＝0时，Ft为初始状态转移矩阵；其中，ωx t表示在t时刻把飞行距离转化成经度需要乘上的系数，ωy表示把飞行距离转化成纬度需要乘上的系数，θt为t时刻飞机的朝向角，Δt为QAR数据的采样间隔，上述ωx t，ωy为：
Figure FDA0003070930170000021
Figure FDA0003070930170000022
其中，Latt表示t时刻飞机所处的纬度，Earth_Radius_Long表示地球长半轴半径，Earth_Radius_Short表示地球短半轴半径；
步骤1所述t时刻输入传递矩阵为：
Figure FDA0003070930170000023
其中，Bt为t时刻的输入传递矩阵，ωx t表示在t时刻把飞行距离转化成经度需要乘上的系数，ωy表示把飞行距离转化成纬度需要乘上的系数，由于在对轨迹进行采样时Δt足够小，根据牛顿运动学定则，此时在做近似匀加速直线运动的状态转移方程为：
Xt+1＝Xt+(Vt*Δt+0.5*at*Δt2)*Sin(θt)*ωx
Yt+1＝Yt+(Vt*Δt+0.5*at*Δt2)*Cos(θt)*ωy
Vt+1＝Vt+at*Δt
其中，at表示t时刻飞机的水平加速度，Vt*Δt+0.5*at*Δt2表示Δt内飞机的水平位移，Vt为t时刻飞机的水平速度，Δt为离散时间的采样间隔。
3.根据权利要求1所述的基于卡尔曼滤波的QAR数据噪声处理方法，其特征在于，
步骤2所述t+1时刻预测状态向量为：
Statet+1 -＝Ft*Statet+Bt*at+xt，0≤t≤K-1
式中，Ft为t时刻的状态转移矩阵，Bt为t时刻的输入传递矩阵，at为t时刻飞机的水平加速度，Ft*Statet表示飞机在两个时刻间的状态转移过程，Bt*at表示水平加速度对飞机状态的影响；xt为t时刻系统状态向量误差，该误差使得飞机动力装置对飞机的控制与记录数据之间存在偏差，符合正态分布xt～N(0，Qt)，Qt为t时刻状态向量误差的协方差矩阵；
步骤2所述t+1时刻飞机状态预测误差协方差矩阵为：
Cov(Statet+1)＝Cov(Ft*Statet)＝FtPtFt T，0≤t≤K-1
步骤2中所述t+1时刻转移状态误差协方差矩阵为：
Pt+1 -＝FtPtFt T+Qt，0≤t≤K-1
式中，Pt+1 -为t+1时刻的转移状态误差协方差矩阵，Pt为t时刻的状态误差协方差矩阵，Ft为t时刻的状态转移矩阵，Qt为t时刻状态向量误差的协方差矩阵。
4.根据权利要求1所述的基于卡尔曼滤波的QAR数据噪声处理方法，其特征在于，
步骤3所述t+1时刻状态向量观测值为：
Zt+1＝H*Statet+1 -+R，0≤t≤K-1
其中，Zt+1为t+1时刻状态向量观测值，H为变换矩阵，R为观测误差协方差矩阵；
步骤3所述t+1时刻观测值向量为：
Figure FDA0003070930170000031
其中，Lngt+1为t+1时刻观测到的飞机所在经度，Latt+1为t+1时刻观测到的飞机所在纬度；
步骤3所述t+1时刻卡尔曼系数为：
Kt+1＝Pt+1 -*HT*(H*Pt+1 -*HT+R)-1，0≤t≤K-1
其中，Pt+1 -为t+1时刻转移状态误差协方差矩阵，H为变换矩阵，R为观测误差协方差矩阵；卡尔曼系数的主要作用是衡量根据状态协方差矩阵Pt+1 -和观测状态协方差矩阵R的大小，决定卡尔曼滤波得到的估计值是更偏向于预测值还是更偏向观测值。
5.根据权利要求1所述的基于卡尔曼滤波的QAR数据噪声处理方法，其特征在于，
步骤4所述滤波后的t+1时刻状态为：
Statet+1＝Statet+1 -+Kt+1*(Zt+1-H*Statet+1 -)，0≤t≤K-1
式中，Statet+1 -为t+1时刻的预测状态，Zt+1为t+1时刻的观测状态，H为变化矩阵，Kt+1为t+1时刻的卡尔曼系数；Zt+1-H*Statet+1 -为实际的观测与预测值之间的残差，其与Kt+1的乘积用来修正飞机位置的预测值；
步骤4所述t+1时刻实际状态转移协方差矩阵为：
Pt+1＝(I-Kt*H)*Pt+1 -，0≤t≤K-1
其中，I为单位矩阵，Kt为t时刻的卡尔曼系数，Pt+1 -为t+1时刻转移状态误差协方差矩阵，H为变换矩阵；随着迭代的进行，最佳状态估计值的噪声分布情况需要不断地进行更新；状态的不确定性随着迭代的进行总体上会逐渐变小，但是由于传递噪声的存在，状态的不确定性又会呈现出逐渐变大的趋势；卡尔曼滤波算法就是通过不断地调整从而在这一状态中达到一个平衡。
CN202110538773.5A 2021-05-18 2021-05-18 一种基于卡尔曼滤波的qar数据噪声处理方法 Active CN113311456B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110538773.5A CN113311456B (zh) 2021-05-18 2021-05-18 一种基于卡尔曼滤波的qar数据噪声处理方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110538773.5A CN113311456B (zh) 2021-05-18 2021-05-18 一种基于卡尔曼滤波的qar数据噪声处理方法

Publications (2)

Publication Number Publication Date
CN113311456A true CN113311456A (zh) 2021-08-27
CN113311456B CN113311456B (zh) 2022-08-16

Family

ID=77373571

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110538773.5A Active CN113311456B (zh) 2021-05-18 2021-05-18 一种基于卡尔曼滤波的qar数据噪声处理方法

Country Status (1)

Country Link
CN (1) CN113311456B (zh)

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20140139374A1 (en) * 2012-11-21 2014-05-22 Raytheon Company Kalman filtering with indirect noise measurements
CN105719314A (zh) * 2016-01-30 2016-06-29 西北工业大学 基于单应性估计和扩展卡尔曼滤波的无人机位置估计方法
CN107133635A (zh) * 2017-03-30 2017-09-05 南京航空航天大学 一种多源异构飞行事故航迹数据融合方法
CN109472418A (zh) * 2018-11-16 2019-03-15 西安电子科技大学 基于卡尔曼滤波的机动目标状态预测优化方法
CN110081878A (zh) * 2019-05-17 2019-08-02 东北大学 一种多旋翼无人机的姿态及位置确定方法

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20140139374A1 (en) * 2012-11-21 2014-05-22 Raytheon Company Kalman filtering with indirect noise measurements
CN105719314A (zh) * 2016-01-30 2016-06-29 西北工业大学 基于单应性估计和扩展卡尔曼滤波的无人机位置估计方法
CN107133635A (zh) * 2017-03-30 2017-09-05 南京航空航天大学 一种多源异构飞行事故航迹数据融合方法
CN109472418A (zh) * 2018-11-16 2019-03-15 西安电子科技大学 基于卡尔曼滤波的机动目标状态预测优化方法
CN110081878A (zh) * 2019-05-17 2019-08-02 东北大学 一种多旋翼无人机的姿态及位置确定方法

Non-Patent Citations (4)

* Cited by examiner, † Cited by third party
Title
WANG QING ET AL.: "Aerodynamic Modeling and Parameter Estimation from QAR Data of an Airplane Approaching a High-altitude Airport", 《CHINESE JOURNAL OF AERONAUTICS 25 (2012)》 *
刘家学等: "高阶无迹卡尔曼滤波算法在飞机定位中的应用", 《计算机应用与软件》 *
苏义鑫等: "基于自适应四元数卡尔曼滤波的姿态估计方法", 《武汉理工大学学报》 *
谷润平等: "QAR数据的数据融合算法", 《计算机系统应用》 *

Also Published As

Publication number Publication date
CN113311456B (zh) 2022-08-16

Similar Documents

Publication Publication Date Title
CN102435763B (zh) 一种基于星敏感器的航天器姿态角速度测量方法
CN109343550B (zh) 一种基于滚动时域估计的航天器角速度的估计方法
CN103940433B (zh) 一种基于改进的自适应平方根ukf算法的卫星姿态确定方法
CN105737823A (zh) 基于五阶ckf的gps/sins/cns组合导航方法
CN109443342A (zh) 新型自适应卡尔曼无人机姿态解算方法
CN107367941B (zh) 高超声速飞行器攻角观测方法
CN110375731A (zh) 一种混合交互式多模型滤波方法
CN103776449A (zh) 一种提高鲁棒性的动基座初始对准方法
CN102436437B (zh) 基于角速度的飞行器极限飞行时四元数傅里埃近似输出方法
CN113311456B (zh) 一种基于卡尔曼滤波的qar数据噪声处理方法
CN106931966B (zh) 一种基于泰勒高阶余项拟合的组合导航方法
CN102495829B (zh) 基于角速度的飞行器极限飞行时四元数沃尔什近似输出方法
CN102508819B (zh) 基于角速度的飞行器极限飞行时四元数勒让德近似输出方法
CN102519466A (zh) 基于角速度的欧拉角勒让德指数近似输出方法
CN102495830B (zh) 基于角速度的飞行器极限飞行时四元数Hartley近似输出方法
CN102445202B (zh) 一种刚体空间运动状态的拉盖尔输出方法
CN102495831B (zh) 基于角速度的飞行器极限飞行时四元数埃米特近似输出方法
Yu et al. Hybrid multi-frequency attitude estimation based on vision/inertial integrated measurement system
CN102323990B (zh) 一种刚体空间运动气动模型的建模方法
CN102506866B (zh) 基于角速度的飞行器极限飞行时四元数切比雪夫近似输出方法
Garcia et al. Adaptive and resilient flight control system for a small unmanned aerial system
CN102506865B (zh) 基于角速度的飞行器极限飞行时四元数多项式类近似输出方法
Zhu et al. Research on GFSINS/Star-sensor Integrated Attitude Estimation Algorithm Based on UKF.
CN102506864B (zh) 飞行器极限飞行时四元数任意步长正交级数近似输出方法
CN102508818B (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
GR01 Patent grant
GR01 Patent grant