CN107229060A - 一种基于自适应滤波的gps测量数据处理方法 - Google Patents

一种基于自适应滤波的gps测量数据处理方法 Download PDF

Info

Publication number
CN107229060A
CN107229060A CN201710497032.0A CN201710497032A CN107229060A CN 107229060 A CN107229060 A CN 107229060A CN 201710497032 A CN201710497032 A CN 201710497032A CN 107229060 A CN107229060 A CN 107229060A
Authority
CN
China
Prior art keywords
mrow
msub
mover
time
target
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
CN201710497032.0A
Other languages
English (en)
Other versions
CN107229060B (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.)
Beijing Technology and Business University
Original Assignee
Beijing Technology and Business University
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 Beijing Technology and Business University filed Critical Beijing Technology and Business University
Priority to CN201710497032.0A priority Critical patent/CN107229060B/zh
Publication of CN107229060A publication Critical patent/CN107229060A/zh
Application granted granted Critical
Publication of CN107229060B publication Critical patent/CN107229060B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

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/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

Landscapes

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

Abstract

本发明提供一种基于自适应滤波的GPS测量数据处理方法，该方法是一种基于二阶自适应模型的改进的卡尔曼滤波器的方法。首先，采用了递归卡尔曼滤波对目标进行实时跟踪；其次，运用二阶自适应模型来抓取目标的运动特性，同时也去除估计过程中的有色噪音；最后，为了求出测量噪音方差真值，引入遗忘因子，并采用自适应指数平滑的误差补偿方法，使其测量噪音方差快速收敛于真值，但是，因为在收敛过后的值抖动较大，所以又引入滑动窗口，并且运用一种自适应平滑滤波的误差补偿方法，此方法收敛速度很慢，但收敛完后十分稳定，所以可以使其测量误差方差一直在真值附近平滑的波动，从而使跟踪效果更加精确，解决了GPS测量数据的高精度在线去噪问题。

Description

一种基于自适应滤波的GPS测量数据处理方法
技术领域
本发明涉及测量信号处理技术领域，尤其涉及一种基于自适应滤波的GPS测量数据处理方法。
背景技术
现如今，卫星定位系统得到了普遍的关注，发展迅速，广泛地应用于航天、航空、航海各军事领域以及交通、测量等民用领域。但是，仅仅靠一般的GPS测量数据进行跟踪的精度有限，无法满足实际需求；而精度较高的GPS系统的价格极其昂贵，因此，处理带有噪音的精度有限的GPS测量数据，提高其跟踪精度非常重要。
为了提高跟踪精度，一个性能优良的处理GPS测量数据的滤波算法是必要的，目前已经有很多的算法，如三维匹配滤波，基于Hough变化的滤波算法，动态规划算法和粒子滤波算法。首先，三维匹配滤波算法适合于数字处理，且可用于多目标的情况，但必须已知目标速度且该速度必须为匀速，当目标速度未知或者当目标速度机动时，算法性能严重下降，所以，在实际系统中，面对强机动性且带有有色噪音的GPS数据，并不适用。其次，Hough变化虽然在传统方法的基础上有所改进，但是同样在面临目标轨迹不是直线时，性能会明显下降。然后，对于动态规划算法而言，有很多学者提出了各种各样的改进算法，如2008年Buzzi等人提出的基于广义似然的动态规划算法，可以把目标跟踪问题看成一个路径搜索问题，但是并没有考虑到实际问题中有色噪音的存在，只是做了仿真实验，而2010年Orlando等人提出了一种基于空间自适应处理的动态规划算法，此后，也有学者提出了基于极大似然-概率数据互联的动态规划算法等等，但是，无论哪种动态规划算法对于全局的最优解都需要全局的数据，不能实现实时数据处理。最后，粒子滤波的算法虽然也可以针对于强机动性的目标，但是，其需要样本数量巨大，并且算法复杂，耗时长，目标的机动性越强，算法的复杂度越高，不利于快速跟踪目标。
由此可见，提出一种可用于实际系统中有高强机动性且含有有色噪音的GPS测量数据的快速处理算法成为亟待解决的课题。
发明内容
本发明的目的在于提供一种基于自适应滤波的GPS测量数据处理方法，其解决了GPS测量数据的高精度在线去噪问题，从而很大程度上增加了跟踪的精度。
本发明涉及一种基于自适应滤波的GPS测量数据处理方法，利用多个自适应模型实时修正参数，解决了GPS测量数据的高精度在线去噪问题，从而很大程度上增加了跟踪的精度。该方法是一种基于二阶自适应模型的改进的卡尔曼滤波器的方法。首先，采用了递归卡尔曼滤波对目标进行实时跟踪；其次，运用二阶自适应模型来抓取目标的运动特性，同时也去除估计过程中的有色噪音；最后，为了求出测量噪音方差真值，引入遗忘因子，并采用自适应指数平滑的误差补偿方法，使其测量噪音方差快速收敛于真值，但是，因为在收敛过后的值抖动较大，所以又引入滑动窗口，并且运用一种自适应平滑滤波的误差补偿方法，此方法收敛速度很慢，但收敛完后十分稳定，所以可以使其测量误差方差一直在真值附近平滑的波动，从而使跟踪效果更加精确。
具体实现步骤如下：
步骤1：目标运动状态和系统自适应参数初始化；
步骤2：建立具有系统自适应参数的二阶自适应模型，求出符合目标运动特性的状态转移矩阵Φ(k+1,k)、输入矩阵U(k)以及有色过程噪音方差Q(k)；其中，系统自适应参数包括机动频率α和变化率方差
步骤3：根据二阶自适应模型和目标初始状态对目标状态进行一步预测，得出下一时刻目标的状态预测值
步骤4：根据二阶自适应模型和过程噪音计算向前一步估计方差P(k+1|k)；
步骤5：引入遗忘因子dk，dk＝(1-b)/(1-bk+1)，b∈(0.1)；
步骤6：根据遗忘因子dk、测量数据y(k+1)和状态预测值计算出测量噪音的均值；
步骤7：根据测量数据y(k+1)和状态预测值计算出估计过程中的残差e(k)；
步骤8：利用求得的测量噪音均值、残差e(k)以及向前一步估计方差P(k+1|k)来计算出测量噪音方差；
步骤9：根据向前一步估计方差P(k+1|k)和测量噪音方差求出卡尔曼滤波增益K(k+1)；
步骤10：根据状态预测值测量数据y(k+1)和卡尔曼滤波增益对目标状态进行更新，得到目标状态估计值
步骤11：通过卡尔曼滤波增益和向前一步估计方差P(k+1|k)得到更新后的估计方差P(k+1|k+1)；
步骤12：根据目标状态估计值计算出目标的速度均值及速度估计值；
步骤13：利用Yule-Walker方法，通过目标状态估计值更新系统自适应参数，进而更新步骤2中的二阶自适应模型；
步骤14：重复步骤2至步骤13，直到处理完GPS测量数据第100个值时，进入步骤15；
步骤15：从k＝100开始，依次进行步骤2、步骤3、步骤4，再通过测量数据y(k+1)和状态预测值计算出估计过程中的残差e(k)；
步骤16：引入滑动窗口；
步骤17：通过滑动窗口和残差计算出窗口内残差均值
步骤18：通过残差均值和向前一步估计方差求出测量噪音R(k)；
步骤19：根据向前一步估计方差和测量噪音求出卡尔曼滤波增益K(k+1)；
步骤20：根据向前一步状态预测值、测量数据和卡尔曼滤波增益更新状态估计值
步骤21：通过卡尔曼滤波增益和向前一步估计方差，得到更新后的估计方差P(k+1|k+1)；
步骤22：利用步骤12和步骤13的方法，通过目标状态估计值更新自适应当前模型参数；
步骤23：重复步骤15至步骤22，直至所有测量数据全部执行完毕，则结束。
在上述技术方案的基础上，本发明还可以做如下改进。
进一步地，步骤1中目标运动状态和系统自适应参数初始化包括：
设置初值，从k＝0开始，k为采样时刻；
目标运动状态的初始值
估计方差的初始值P(0|0)＝P0；
二阶自适应模型参数初值：α(0)＝α0、
过程噪音方差初值Q(0)＝Q0；
测量噪音均值初值：
测量噪音方差初值：R(0)＝R0。
进一步地，步骤6的计算公式如下：
表示从0时刻开始至k时刻的测量噪音均值，表示从0时刻开始至k+1时刻的测量噪音均值，dk为遗忘因子，y(k+1)为测量数据，表示k时刻预测目标在k+1时刻的状态，k为采样时刻，|表示条件操作符，H(k+1)表示k+1时刻的测量矩阵。
进一步地，步骤7的计算公式如下：
ek+1为估计过程中的残差，y(k+1)为测量数据，表示k时刻预测目标在k+1时刻的状态，k为采样时刻，|表示条件操作符，H(k+1)表示k+1时刻的测量矩阵，表示从0时刻开始至k时刻的测量噪音均值。
进一步地，步骤8的计算公式如下：
表示k+1时刻的测量噪音方差，表示k时刻的测量噪音方差，dk为遗忘因子，ek+1表示估计过程中的残差，P(k+1|k)表示向前一步估计方差，H(k+1)表示k+1时刻的测量矩阵，HT(k+1)表示k+1时刻的测量矩阵的转置，表示k时刻预测目标在k+1时刻的状态，k为采样时刻，|表示条件操作符。
进一步地，步骤9的计算公式如下：
K(k+1)表示卡尔曼滤波增益，P(k+1|k)表示向前一步估计方差，|表示条件操作符，表示k+1时刻的测量噪音方差，H(k+1)表示k+1时刻的测量矩阵，HT(k+1)表示k+1时刻的测量矩阵的装置。
进一步地，步骤10的计算公式如下：
表示状态估计值，y(k+1)为测量数据，表示k时刻预测目标在k+1时刻的状态，k为采样时刻，|表示条件操作符，K(k+1)表示卡尔曼滤波增益，H(k+1)表示k+1时刻的测量矩阵。
进一步地，步骤12中根据目标状态估计值计算出目标的速度均值及速度估计值的具体步骤如下；
步骤12.1：利用下式计算目标速度均值：
其中为0至k+1时刻速度均值，为k时刻目标的状态估计值的第二行值，k为采样时刻；
步骤12.2：按照下式获取系统k时刻和k+1时刻的加速度估计值
其中为k时刻状态估计的第二行值，为k+1时刻状态估计的第二行值。
进一步地，步骤13中利用Yule-Walker方法，通过目标状态估计值更新系统自适应参数的具体方法如下；
根据目标速度的估计值对系统自适应参数进行修正，根据采样时刻k值的大小，选择修正系统自适应参数α和若k小于等于4进入步骤13.1，若k大于4进入步骤13.2；
13.1当采样时刻k小于等于4时，按下式计算取系统自适应参数α和α＝α0,其中α0为系统自适应参数α的初值，
如果则取
如果则取
如果则取(0,10]之间的任意数，
其中，为k时刻目标加速度估计值，π为圆周率，取为3.14，vM为正的常数，取为3，v-M为与vM绝对值相等的负常数，取为-3；
13.2当采样时刻k大于4时，利用Yule-Walker方法，按下式计算系统自适应参数α和
rk+1(1)为k+1时刻目标速度向前一步相关函数，rk(1)为k时刻目标速度向前一步相关函数，和分别为k时刻和k+1时刻目标速度估计值；rk+1(0)为k+1时刻目标速度自相关函数，rk(0)为k时刻目标速度自关函数；
根据系统方程得到目标运动速度，方程满足如下一阶马尔科夫随机序列：
其中为k+1时刻的速度，为k时刻的速度，β为离散后速度随机序列的机动频率，wv(k)为零均值白噪声离散序列，方差为其中为零均值白噪声w(t)的方差，β与α的关系为β＝e-αT；x(k+1)为k+1时刻目标的状态向量，k为取样时刻，Φ(k+1,k)表示状态转移矩阵，x(k)为k时刻目标的状态向量，U(k)表示输入矩阵，为0时刻开始至k时刻目标的速度均值，w(k)表示过程噪音，其均值为0，方差为Q(k)，Φ(k+1,k)、U(k)和Q(k)中含有目标机动频率α和变化率方差随着系统自适应参数的变化而变化。
一阶马尔科夫时间加速度序列满足以下参数关系：
其中，α与β分别为加速度的机动频率及其离散化后加速度序列的机动频率，自适应参数α和可按照下式计算得到：
其中，ln为取以e为底的对数计算；α和为系统自适应参数，T为采样间隔。
进一步地，步骤18的计算公式如下：
其中k为采样时刻，ek+1表示估计过程中的残差，表示k+1时刻的残差均值，H(k+1)表示k+1时刻的测量矩阵，P(k+1|k)表示向前一步估计方差，|表示条件操作符，H(k+1)表示k+1时刻的测量矩阵，N为残差的个数。
本发明的有益效果是：
本发明的基于自适应滤波的GPS测量数据处理方法，可以很好的处理实际系统中具有强机动性和有色噪音的GPS测量数据，具有以下优点：
1.针对GPS测量数据实时处理的问题，利用卡尔曼递归估计算法实现对GPS数据的实时处理。
2.针对GPS测量数据机动性强的问题，利用二阶自适应模型，来抓取目标的运动性能，与此同时，去除估计过程中的有色噪音。
3.针对GPS测量数据的测量噪音未知的问题，引入了两种自适应滤波方法一起估计测量噪音的真值，既提高了精度，又提高了运行速度。
4.针对平滑滤波处理方法的收敛极慢的问题，首先采用自适应指数平滑的方式使其测量噪音快速收敛，即引入了遗忘因子，通过自适应指数平滑滤波的误差补偿方法，求出测量噪音。
5.针对自适应指数平滑滤波处理测量方差所得的值不稳定的问题，在得到快速收敛的测量方差值过后，采用一种运算小稳定性高的自适应平滑滤波算法，即引入滑动窗口，利用一种自适应平滑滤波的误差补偿方法，将测量方差稳定收敛于真值。
附图说明
图1为本发明一实施例的基于自适应滤波的GPS测量数据处理方法的流程图；
图2为测量数据的波形图；
图3为估计数据的波形图；
图4为测量值和估计值对比示意图。
具体实施方式
为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例对本发明的基于自适应滤波的GPS测量数据处理方法进行进一步详细说明。需要说明的是，在不冲突的情况下，以下各实施例及实施例中的特征可以相互组合。应当理解，此处所描述的具体实施例仅用于解释本发明，并不用于限定本发明。
参照图1，本发明一实施例的基于自适应滤波的GPS测量数据处理方法包括以下步骤：
步骤S1：目标运动状态和系统自适应参数初始化；
步骤S2：建立具有系统自适应参数的二阶自适应模型，求出符合目标运动特性的状态转移矩阵Φ(k+1,k)、输入矩阵U(k)以及有色过程噪音方差Q(k)；其中，系统自适应参数包括机动频率α和变化率方差
步骤S3：根据二阶自适应模型和目标初始状态对目标状态进行一步预测，得出下一时刻目标的状态预测值
步骤S4：根据二阶自适应模型和过程噪音计算向前一步估计方差P(k+1|k)；
步骤S5：引入遗忘因子dk，dk＝(1-b)/(1-bk+1)，b∈(0.1)；
步骤S6：根据遗忘因子dk、测量数据y(k+1)和状态预测值计算出测量噪音的均值；
步骤S7：根据测量数据y(k+1)和状态预测值计算出估计过程中的残差e(k)；
步骤S8：利用求得的测量噪音均值、残差e(k)以及向前一步估计方差P(k+1|k)来计算出测量噪音方差；
步骤S9：根据向前一步估计方差P(k+1|k)和测量噪音方差R(k)求出卡尔曼滤波增益K(k+1)；
步骤S10：根据状态预测值测量数据y(k+1)和卡尔曼滤波增益对目标状态进行更新，得到目标状态估计值
步骤S11：通过卡尔曼滤波增益和向前一步估计方差P(k+1|k)得到更新后的估计方差P(k+1|k+1)；
步骤S12：根据目标状态估计值计算出目标的速度均值及速度估计值；
步骤S13：利用Yule-Walker方法，通过目标状态估计值更新系统自适应参数，进而更新步骤S2中的二阶自适应模型；
步骤S14：重复步骤S2至步骤S13，直到处理完GPS测量数据第100个值时，进入步骤S15；
步骤S15：从k＝100开始，依次进行步骤S2、步骤S3、步骤S4，再通过测量数据y(k+1)和状态预测值计算出估计过程中的残差e(k)；
步骤S16：引入滑动窗口；
步骤S17：通过滑动窗口和残差计算出窗口内残差均值
步骤18：通过残差均值和向前一步估计方差求出测量噪音R(k)；
步骤S19：根据向前一步估计方差和测量噪音求出卡尔曼滤波增益K(k+1)；
步骤S20：根据向前一步状态预测值、测量数据和卡尔曼滤波增益更新状态估计值
步骤S21：通过卡尔曼滤波增益和向前一步估计方差，得到更新后的估计方差P(k+1|k+1)；
步骤S22：利用步骤S12和步骤S13的方法，通过目标状态估计值更新自适应当前模型参数；
步骤S23：重复步骤S15至步骤S22，直至所有测量数据全部执行完毕，则结束。
以上各步骤具体如下：
步骤S1：设置初值，从k＝0开始，k为采样时刻；
目标运动状态的初始值
估计方差的初始值P(0|0)＝P0；
二阶自适应模型参数初值：α(0)＝α0、
过程噪音方差初值Q(0)＝Q0；
测量噪音均值初值：
测量噪音方差初值：R(0)＝R0。
本实例中：x0＝0；P0＝100*eye(3)；α0＝1/20；Q0＝anyvalue；r0＝anyvalue；R0＝anyvalue；
步骤S2：建立二阶自适应模型：
为了更好的抓取目标的运动特性，设系统速度为非零均值的时间相关随机过程：其中为速度均值，v(t)为零均值指数相关有色噪音模型，其相关函数为：
式中：Rv(τ)表示相关函数，表示速度的方差，α为机动频率，反应数据的波动特性。
接着对有色噪音v(t)做白化处理，得到：
式中：表示有色噪音v(t)的一阶导数；w(t)为零均值白噪音，方差为
由和得到目标运动的连续方程：
其中
以周期T采样，离散化后系统目标运动满足以下离散时间方程：
其为三维状态列向量，分别表示位移、速度和加速度，x(k+1)为k+1时刻目标的状态向量，k为取样时刻，Φ(k+1,k)表示状态转移矩阵，x(k)为k时刻目标的状态向量，U(k)表示输入矩阵，为0时刻开始至k时刻目标的速度均值，w(k)表示过程噪音，其均值为0，方差为Q(k)，Φ(k+1,k)、U(k)和Q(k)中含有目标机动频率α和变化率方差随着系统自适应参数的变化而变化。
则：
其中：
步骤S3：根据自适应当前统计模型和系统状态初值计算下一时刻的状态预测值
表示k时刻预测目标在k+1时刻的状态，k为采样时刻，|表示条件操作符，Φ(k+1,k)表示状态转移矩阵，表示目标在k时刻的状态估计值，特别的在k＝0时刻为初值U(k)表示输入矩阵，为速度均值。
步骤S4：基于自适应当前统计模型和过程噪音计算向前一步估计方差P(k+1|k)：
P(k+1|k)＝Φ(k+1,k)P(k|k)ΦT(k+1,k)+Q(k)
式中，P(k+1|k)表示向前一步估计方差，即k时刻预测目标在k+1时刻的状态协方差，k为采样时刻，|表示条件操作符；Φ(k+1,k)表示状态转移矩阵，U(k)表示输入矩阵，P(k|k)表示k时刻目标的状态协方差的估计值，Q(k)为过程噪音协方差。
步骤S5：引入遗忘因子dk，dk＝(1-b)/(1-bk+1)，b∈(0.1)，在本次实验中b的值取为0.8；
步骤S6：根据遗忘因子dk，测量数据y(k+1)和一步预测值计算测量噪音的均值：
式中，表示从0时刻开始至k+1时刻的测量噪音均值；dk为遗忘因子，表示从0时刻开始至k时刻的测量噪音均值，y(k+1)为测量数据，表示k时刻预测目标在k+1时刻的状态，k为采样时刻，|表示条件操作符；H(k+1)表示k+1时刻的测量矩阵。
步骤S7：通过测量数据y(k+1)和一步预测值计算估计过程中的残差e(k)：
式中，ek+1为估计过程中的残差，y(k+1)为测量数据，表示k时刻预测目标在k+1时刻的状态，k为采样时刻，|表示条件操作符；H(k+1)表示k+1时刻的测量矩阵，表示从0时刻开始至k时刻的测量噪音均值。
步骤S8：利用求得的测量噪音均值和残差以及向前一步估计方差来计算出测量噪音方差；
式中，表示k+1时刻的测量噪音方差，dk为遗忘因子，表示k时刻的测量噪音方差，ek+1表示估计过程中的残差，P(k+1|k)表示向前一步估计方差，H(k+1)表示k+1时刻的测量矩阵，HT(k+1)表示k+1时刻的测量矩阵的转置，表示k时刻预测目标在k+1时刻的状态，k为采样时刻，|表示条件操作符。
步骤S9：根据向前一步估计方差和测量噪音方差求出卡尔曼滤波增益K(k+1)；
式中，K(k+1)表示卡尔曼滤波增益，P(k+1|k)表示向前一步估计方差，|表示条件操作符，表示k+1时刻的测量噪音方差，H(k+1)表示k+1时刻的测量矩阵，HT(k+1)表示k+1时刻的测量矩阵的装置。
步骤S10：根据向前一步状态预测值，测量数据和卡尔曼滤波增益更新状态估计值
表示状态估计值，y(k+1)为测量数据，表示k时刻预测目标在k+1时刻的状态，k为采样时刻，|表示条件操作符，K(k+1)表示卡尔曼滤波增益，H(k+1)表示k+1时刻的测量矩阵。
步骤S11：通过卡尔曼滤波增益和向前一步估计方差，得到更新后的估计方差P(k+1|k+1)；
步骤S12：根据目标状态估计值计算目标速度均值及速度估计值，具体如下：
利用下式计算目标速度均值：
其中为0至k+1时刻速度均值，为k时刻目标的状态估计值的第二行值，k为采样时刻；
按照下式获取系统k时刻和k+1时刻的加速度估计值
其中，为k时刻状态估计的第二行值，为k+1时刻状态估计的第二行值；
步骤S13：利用Yule-Walker方法，通过状态估计值更新自适应当前模型参数，具体如下：
根据目标速度的估计值对系统自适应参数进行修正，根据采样时刻k值的大小，选择修正系统自适应参数α和若k小于等于4进入步骤13.1，若k大于4进入步骤13.2。
13.1当采样时刻k小于等于4时，因为采样数据较少，采用以下方法求取系统自适应参数α和
α＝α0,其中α0为系统自适应参数α的初值，
如果则取
如果则取
如果则取(0,10]之间的任意数，
其中，为k时刻目标加速度估计值，π为圆周率，取为3.14，vM为正的常数，取为3，v-M为与vM绝对值相等的负常数，取为-3；
13.2当采样时刻k大于4时，利用Yule-Walker方法，按下式计算系统自适应参数α和
其中，rk+1(1)为k+1时刻目标速度向前一步相关函数，rk(1)为k时刻目标速度向前一步相关函数，和分别为k时刻和k+1时刻目标速度估计值；rk+1(0)为k+1时刻目标速度自相关函数，rk(0)为k时刻目标速度自关函数；
根据系统方程得到目标运动速度，方程满足如下一阶马尔科夫随机序列：
其中，为k+1时刻的速度，为k时刻的速度，β为离散后速度随机序列的机动频率，wv(k)为零均值白噪声离散序列，方差为其中为零均值白噪声w(t)的方差，β与α的关系为β＝e-αT；
一阶马尔科夫时间加速度序列满足以下参数关系：
其中，α与β分别为加速度的机动频率及其离散化后加速度序列的机动频率，自适应参数α和可按照下式计算得到：
其中，ln为取以e为底的对数计算；α和为系统自适应参数，T为采样间隔；
步骤S14：重复步骤S2至步骤S13，直到处理完GPS测量数据第100个值时，进入步骤S15。
步骤S15：从k＝100开始，依次进行步骤S2、步骤S3、步骤S4，再通过测量数据y(k+1)和一步预测值计算估计过程中的残差e(k)；
式中，e(k+1)为估计过程中的残差，y(k+1)为测量数据，H(k+1)表示k+1时刻的测量矩阵，表示k时刻预测目标在k+1时刻的状态，k为采样时刻，|表示条件操作符。
步骤S16：引入滑动窗口，需取一个较大的数，这个窗口的大小对结果影响不大，本实验中滑动窗口长度取100。
步骤S17：优化残差值，利用平滑滤波的原理得到残差均值
其中，表示k+1时刻的残差均值，k为采样时刻，表示k时刻的残差均值，e(k+1)表示估计过程中的残差，N为残差的个数。
步骤S18：通过残差均值和向前一步估计方差求出测量噪音R(k+1)：
其中，e(k+1)表示估计过程中的残差，表示k+1时刻的残差均值，H(k+1)表示k+1时刻的测量矩阵，P(k+1|k)表示向前一步估计方差，|表示条件操作符，k为采样时刻；H(k+1)表示k+1时刻的测量矩阵，N为残差的个数，残差的个数与样本中数据的数量相等。
步骤S19：根据向前一步估计方差和测量噪音求出卡尔曼滤波增益K(k+1)：
K(k+1)＝P(k+1|k)HT(k+1)[H(k+1)P(k+1|k)HT(k+1)+R]-1
K(k+1)为k+1时刻的卡尔曼滤波增益，k为采样时刻，P(k+1|k)表示向前一步估计方差，|表示条件操作符；H(k+1)表示k+1时刻的测量矩阵，HT(k+1)表示k+1时刻的测量矩阵的装置，R为高斯测量白噪音的方差。
步骤S20：根据向前一步状态预测值、测量数据和卡尔曼滤波增益更新状态估计值
其中，表示状态估计值，y(k+1)为测量数据，表示k时刻预测目标在k+1时刻的状态，k为采样时刻，|表示条件操作符；K(k+1)表示k+1时刻的卡尔曼滤波增益，H(k+1)表示k+1时刻的测量矩阵。
步骤21：通过卡尔曼滤波增益和向前一步估计方差，得到更新后的估计方差：
P(k+1|k+1)＝[I-K(k+1)H(k+1)]P(k+1|k)
其中，I是3维矩阵，P(k+1|k+1)表示k+1时刻的目标状态估计方差，k为采样时刻，K(k+1)为k+1时刻的卡尔曼滤波增益，P(k+1|k)表示向前一步估计方差，|表示条件操作符，H(k+1)表示k+1时刻的测量矩阵。
步骤S22：利用步骤S12和步骤S13的方法，通过状态估计值更新自适应当前模型参数：
步骤S23：重复步骤S15至步骤S22，直至所有测量数据全部执行完毕，则结束。
实施例1：
本发明针对GPS测量数据进行在线去噪，选用7000个跟踪点，具体数据见下表1，时间步长为0.01s，其中X,Y轴分别为横轴坐标值和纵轴坐标值。
表1
序号 t(s) X轴(m) Y轴(m) 序号 t(s) X轴(m) Y轴(m)
1 0 0 0 4263 42.620 -31.2230 81.5673
2 0.010 -81.4130 -80.0434 4264 42.630 -37.1058 85.9250
3 0.020 -95.0864 -84.9268 4265 42.640 -36.4918 83.8843
4 0.030 -81.8948 -95.7984 4266 42.650 -35.8594 87.9904
5 0.040 -85.3587 -85.5665 4267 42.660 -31.5501 82.2402
6 0.050 -65.0065 -75.4497 4268 42.670 -39.8242 76.7223
7 0.060 -78.7037 --74.1371 … … … …
8 0.070 -78.4814 -89.6396 6996 99.950 86.5073 -96.0795
… … … … 6997 99.960 100.1355 -94.6065
4260 42.590 -37.2217 78.2453 6998 99.970 95.4480 -93.0502
4261 42.600 -37.1711 85.4215 6999 99.980 85.8574 -102.1952
4262 42.610 -43.3685 75.4372 7000 99.990 97.7219 -92.6440
按上述具体实施方式的步骤进行仿真验证，验证结果见附图2-附图4。在附图2-附图4中：横坐标为X轴,表示横轴坐标值，单位为m(米)；纵坐标为Y轴，表示纵轴坐标值，单位为m(米)。附图2是测量数据的波形图，从附图2中可以看出含很大的噪音。附图3是估计数据的波形图，也就是在过滤掉噪音之后的波形图。附图4为测量值和估计值对比示意图，即附图2中测量数据的波形与附图3中估计数据的波形的对比。在附图4中，测量数据波形中部的浅色连续实线为估计数据的波形。从附图4可以看出，采用本发明的基于自适应滤波的GPS测量数据处理方法，可以很好的处理实际系统中具有强机动性和有色噪音的GPS测量数据，解决了GPS测量数据的高精度在线去噪问题，从而很大程度上增加了跟踪的精度。
以上各实施例的基于自适应滤波的GPS测量数据处理方法是一种基于二阶自适应模型的改进的卡尔曼滤波器的方法。首先，采用了递归卡尔曼滤波对目标进行实时跟踪；其次，运用二阶自适应模型来抓取目标的运动特性，同时也去除估计过程中的有色噪音；最后，为了求出测量噪音方差真值，引入遗忘因子，并采用自适应指数平滑的误差补偿方法，使其测量噪音方差快速收敛于真值，但是，因为在收敛过后的值抖动较大，所以又引入滑动窗口，并且运用一种自适应平滑滤波的误差补偿方法，此方法收敛速度很慢，但收敛完后十分稳定，所以可以使其测量误差方差一直在真值附近平滑的波动，从而使跟踪效果更加精确，
以上基于自适应滤波的GPS测量数据处理方法，利用多个自适应模型实时修正参数，解决了GPS测量数据的高精度在线去噪问题，从而很大程度上增加了跟踪的精度。
以上所述仅为本发明的较佳实施例，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (10)

1.一种基于自适应滤波的GPS测量数据处理方法，其特征在于，包括以下步骤：
步骤1：目标运动状态和系统自适应参数初始化；
步骤2：建立具有系统自适应参数的二阶自适应模型，求出符合目标运动特性的状态转移矩阵Φ(k+1,k)、输入矩阵U(k)以及有色过程噪音方差Q(k)；其中，系统自适应参数包括机动频率α和变化率方差
步骤3：根据二阶自适应模型和目标初始状态对目标状态进行一步预测，得出下一时刻目标的状态预测值
步骤4：根据二阶自适应模型和过程噪音计算向前一步估计方差P(k+1|k)；
步骤5：引入遗忘因子dk，dk＝(1-b)/(1-bk+1)，b∈(0.1)；
步骤6：根据遗忘因子dk、测量数据y(k+1)和状态预测值计算出测量噪音的均值；
步骤7：根据测量数据y(k+1)和状态预测值计算出估计过程中的残差e(k)；
步骤8：利用求得的测量噪音均值、残差e(k)以及向前一步估计方差P(k+1|k)来计算出测量噪音方差；
步骤9：根据向前一步估计方差P(k+1|k)和测量噪音方差求出卡尔曼滤波增益K(k+1)；
步骤10：根据状态预测值测量数据y(k+1)和卡尔曼滤波增益对目标状态进行更新，得到目标状态估计值
步骤11：通过卡尔曼滤波增益和向前一步估计方差P(k+1|k)得到更新后的估计方差P(k+1|k+1)；
步骤12：根据目标状态估计值计算出目标的速度均值及速度估计值；
步骤13：利用Yule-Walker方法，通过目标状态估计值更新系统自适应参数，进而更新步骤2中的二阶自适应模型；
步骤14：重复步骤2至步骤13，直到处理完GPS测量数据第100个值时，进入步骤15；
步骤15：从k＝100开始，依次进行步骤2、步骤3、步骤4，再通过测量数据y(k+1)和状态预测值计算出估计过程中的残差e(k)；
步骤16：引入滑动窗口；
步骤17：通过滑动窗口和残差计算出窗口内残差均值
步骤18：通过残差均值和向前一步估计方差求出测量噪音R(k)；
步骤19：根据向前一步估计方差和测量噪音求出卡尔曼滤波增益K(k+1)；
步骤20：根据向前一步状态预测值、测量数据和卡尔曼滤波增益更新状态估计值
步骤21：通过卡尔曼滤波增益和向前一步估计方差，得到更新后的估计方差P(k+1|k+1)；
步骤22：利用步骤12和步骤13的方法，通过目标状态估计值更新自适应当前模型参数；
步骤23：重复步骤15至步骤22，直至所有测量数据全部执行完毕，则结束。
2.根据权利要求1所述的基于自适应滤波的GPS测量数据处理方法，其特征在于，步骤1中目标运动状态和系统自适应参数初始化包括：
设置初值，从k＝0开始，k为采样时刻；
目标运动状态的初始值
估计方差的初始值P(0|0)＝P0；
二阶自适应模型参数初值：α(0)＝α0、
过程噪音方差初值Q(0)＝Q0；
测量噪音均值初值：
测量噪音方差初值：R(0)＝R0。
3.根据权利要求1所述的基于自适应滤波的GPS测量数据处理方法，其特征在于，步骤6的计算公式如下：
<mrow> <msub> <mover> <mi>r</mi> <mo>^</mo> </mover> <mrow> <mi>k</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> <mo>=</mo> <mrow> <mo>(</mo> <mn>1</mn> <mo>-</mo> <msub> <mi>d</mi> <mi>k</mi> </msub> <mo>)</mo> </mrow> <msub> <mover> <mi>r</mi> <mo>^</mo> </mover> <mi>k</mi> </msub> <mo>+</mo> <msub> <mi>d</mi> <mi>k</mi> </msub> <mrow> <mo>(</mo> <mi>y</mi> <mo>(</mo> <mrow> <mi>k</mi> <mo>+</mo> <mn>1</mn> </mrow> <mo>)</mo> <mo>-</mo> <mi>H</mi> <mo>(</mo> <mrow> <mi>k</mi> <mo>+</mo> <mn>1</mn> </mrow> <mo>)</mo> <mover> <mi>x</mi> <mo>^</mo> </mover> <mo>(</mo> <mrow> <mi>k</mi> <mo>+</mo> <mn>1</mn> <mo>|</mo> <mi>k</mi> </mrow> <mo>)</mo> <mo>)</mo> </mrow> </mrow>
表示从0时刻开始至k时刻的测量噪音均值，表示从0时刻开始至k+1时刻的测量噪音均值，dk为遗忘因子，y(k+1)为测量数据，表示k时刻预测目标在k+1时刻的状态，k为采样时刻，|表示条件操作符，H(k+1)表示k+1时刻的测量矩阵。
4.根据权利要求1所述的基于自适应滤波的GPS测量数据处理方法，其特征在于，步骤7的计算公式如下：
<mrow> <msub> <mi>e</mi> <mrow> <mi>k</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> <mo>=</mo> <mi>y</mi> <mrow> <mo>(</mo> <mi>k</mi> <mo>+</mo> <mn>1</mn> <mo>)</mo> </mrow> <mo>-</mo> <mi>H</mi> <mrow> <mo>(</mo> <mi>k</mi> <mo>+</mo> <mn>1</mn> <mo>)</mo> </mrow> <mover> <mi>x</mi> <mo>^</mo> </mover> <mrow> <mo>(</mo> <mi>k</mi> <mo>+</mo> <mn>1</mn> <mo>|</mo> <mi>k</mi> <mo>)</mo> </mrow> <mo>-</mo> <msub> <mover> <mi>r</mi> <mo>^</mo> </mover> <mi>k</mi> </msub> </mrow>
ek+1为估计过程中的残差，y(k+1)为测量数据，表示k时刻预测目标在k+1时刻的状态，k为采样时刻，|表示条件操作符，H(k+1)表示k+1时刻的测量矩阵，表示从0时刻开始至k时刻的测量噪音均值。
5.根据权利要求1所述的基于自适应滤波的GPS测量数据处理方法，其特征在于，步骤8的计算公式如下：
<mrow> <msub> <mover> <mi>R</mi> <mo>^</mo> </mover> <mrow> <mi>k</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> <mo>=</mo> <mrow> <mo>(</mo> <mn>1</mn> <mo>-</mo> <msub> <mi>d</mi> <mi>k</mi> </msub> <mo>)</mo> </mrow> <msub> <mover> <mi>R</mi> <mo>^</mo> </mover> <mi>k</mi> </msub> <mo>+</mo> <msub> <mi>d</mi> <mi>k</mi> </msub> <mrow> <mo>(</mo> <msub> <mi>e</mi> <mrow> <mi>k</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> <msubsup> <mi>e</mi> <mrow> <mi>k</mi> <mo>+</mo> <mn>1</mn> </mrow> <mi>T</mi> </msubsup> <mo>-</mo> <mi>H</mi> <mo>(</mo> <mrow> <mi>k</mi> <mo>+</mo> <mn>1</mn> </mrow> <mo>)</mo> <mi>P</mi> <mo>(</mo> <mrow> <mi>k</mi> <mo>+</mo> <mn>1</mn> <mo>|</mo> <mi>k</mi> </mrow> <mo>)</mo> <msup> <mi>H</mi> <mi>T</mi> </msup> <mo>(</mo> <mrow> <mi>k</mi> <mo>+</mo> <mn>1</mn> </mrow> <mo>)</mo> <mo>)</mo> </mrow> </mrow>
表示k+1时刻的测量噪音方差，表示k时刻的测量噪音方差，dk为遗忘因子，ek+1表示估计过程中的残差，P(k+1|k)表示向前一步估计方差，H(k+1)表示k+1时刻的测量矩阵，HT(k+1)表示k+1时刻的测量矩阵的转置，表示k时刻预测目标在k+1时刻的状态，k为采样时刻，|表示条件操作符。
6.根据权利要求1所述的基于自适应滤波的GPS测量数据处理方法，其特征在于，步骤9的计算公式如下：
<mrow> <mi>K</mi> <mrow> <mo>(</mo> <mi>k</mi> <mo>+</mo> <mn>1</mn> <mo>)</mo> </mrow> <mo>=</mo> <mi>P</mi> <mrow> <mo>(</mo> <mi>k</mi> <mo>+</mo> <mn>1</mn> <mo>|</mo> <mi>k</mi> <mo>)</mo> </mrow> <msup> <mi>H</mi> <mi>T</mi> </msup> <mrow> <mo>(</mo> <mi>k</mi> <mo>+</mo> <mn>1</mn> <mo>)</mo> </mrow> <msup> <mrow> <mo>&amp;lsqb;</mo> <mi>H</mi> <mrow> <mo>(</mo> <mi>k</mi> <mo>+</mo> <mn>1</mn> <mo>)</mo> </mrow> <mi>P</mi> <mrow> <mo>(</mo> <mi>k</mi> <mo>+</mo> <mn>1</mn> <mo>|</mo> <mi>k</mi> <mo>)</mo> </mrow> <msup> <mi>H</mi> <mi>T</mi> </msup> <mrow> <mo>(</mo> <mi>k</mi> <mo>+</mo> <mn>1</mn> <mo>)</mo> </mrow> <mo>+</mo> <msub> <mover> <mi>R</mi> <mo>^</mo> </mover> <mrow> <mi>k</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> <mo>&amp;rsqb;</mo> </mrow> <mrow> <mo>-</mo> <mn>1</mn> </mrow> </msup> </mrow>
K(k+1)表示卡尔曼滤波增益，P(k+1|k)表示向前一步估计方差，|表示条件操作符，表示k+1时刻的测量噪音方差，H(k+1)表示k+1时刻的测量矩阵，HT(k+1)表示k+1时刻的测量矩阵的装置。
7.根据权利要求1所述的基于自适应滤波的GPS测量数据处理方法，其特征在于，步骤10的计算公式如下：
<mrow> <mover> <mi>x</mi> <mo>^</mo> </mover> <mrow> <mo>(</mo> <mi>k</mi> <mo>+</mo> <mn>1</mn> <mo>|</mo> <mi>k</mi> <mo>+</mo> <mn>1</mn> <mo>)</mo> </mrow> <mo>=</mo> <mover> <mi>x</mi> <mo>^</mo> </mover> <mrow> <mo>(</mo> <mi>k</mi> <mo>+</mo> <mn>1</mn> <mo>|</mo> <mi>k</mi> <mo>)</mo> </mrow> <mo>+</mo> <mi>K</mi> <mrow> <mo>(</mo> <mi>k</mi> <mo>+</mo> <mn>1</mn> <mo>)</mo> </mrow> <mo>&amp;lsqb;</mo> <mi>y</mi> <mrow> <mo>(</mo> <mi>k</mi> <mo>+</mo> <mn>1</mn> <mo>)</mo> </mrow> <mo>-</mo> <mi>H</mi> <mrow> <mo>(</mo> <mi>k</mi> <mo>+</mo> <mn>1</mn> <mo>)</mo> </mrow> <mover> <mi>x</mi> <mo>^</mo> </mover> <mrow> <mo>(</mo> <mi>k</mi> <mo>+</mo> <mn>1</mn> <mo>|</mo> <mi>k</mi> <mo>)</mo> </mrow> <mo>&amp;rsqb;</mo> </mrow> 2
表示状态估计值，y(k+1)为测量数据，表示k时刻预测目标在k+1时刻的状态，k为采样时刻，|表示条件操作符，K(k+1)表示卡尔曼滤波增益，H(k+1)表示k+1时刻的测量矩阵。
8.根据权利要求1所述的基于自适应滤波的GPS测量数据处理方法，其特征在于，步骤12中根据目标状态估计值计算出目标的速度均值及速度估计值的具体步骤如下；
步骤12.1：利用下式计算目标速度均值：
<mrow> <mover> <mi>v</mi> <mo>&amp;OverBar;</mo> </mover> <mrow> <mo>(</mo> <mi>k</mi> <mo>+</mo> <mn>1</mn> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mn>1</mn> <mrow> <mi>k</mi> <mo>+</mo> <mn>1</mn> </mrow> </mfrac> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>k</mi> <mo>=</mo> <mn>0</mn> </mrow> <mi>k</mi> </munderover> <mover> <mover> <mi>x</mi> <mo>&amp;CenterDot;</mo> </mover> <mo>^</mo> </mover> <mrow> <mo>(</mo> <mi>k</mi> <mo>+</mo> <mn>1</mn> <mo>|</mo> <mi>k</mi> <mo>+</mo> <mn>1</mn> <mo>)</mo> </mrow> </mrow>
其中为0至k+1时刻速度均值，为k时刻目标的状态估计值的第二行值，k为采样时刻；
步骤12.2：按照下式获取系统k时刻和k+1时刻的加速度估计值
<mrow> <mover> <mi>v</mi> <mo>^</mo> </mover> <mrow> <mo>(</mo> <mi>k</mi> <mo>)</mo> </mrow> <mo>=</mo> <mover> <mover> <mi>x</mi> <mo>&amp;CenterDot;</mo> </mover> <mo>^</mo> </mover> <mrow> <mo>(</mo> <mi>k</mi> <mo>|</mo> <mi>k</mi> <mo>)</mo> </mrow> </mrow>
<mrow> <mover> <mi>v</mi> <mo>^</mo> </mover> <mrow> <mo>(</mo> <mi>k</mi> <mo>+</mo> <mn>1</mn> <mo>)</mo> </mrow> <mo>=</mo> <mover> <mover> <mi>x</mi> <mo>&amp;CenterDot;</mo> </mover> <mo>^</mo> </mover> <mrow> <mo>(</mo> <mi>k</mi> <mo>+</mo> <mn>1</mn> <mo>|</mo> <mi>k</mi> <mo>+</mo> <mn>1</mn> <mo>)</mo> </mrow> </mrow>
其中为k时刻状态估计的第二行值，为k+1时刻状态估计的第二行值。
9.根据权利要求1所述的基于自适应滤波的GPS测量数据处理方法，其特征在于，步骤13中利用Yule-Walker方法，通过目标状态估计值更新系统自适应参数的具体方法如下；
根据目标速度的估计值对系统自适应参数进行修正，根据采样时刻k值的大小，选择修正系统自适应参数α和若k小于等于4进入步骤13.1，若k大于4进入步骤13.2；
13.1当采样时刻k小于等于4时，按下式计算取系统自适应参数α和α＝α0,其中α0为系统自适应参数α的初值，
如果则取
如果则取
如果则取(0,10]之间的任意数，
其中，为k时刻目标加速度估计值，π为圆周率，取为3.14，vM为正的常数，取为3，v-M为与vM绝对值相等的负常数，取为-3；
13.2当采样时刻k大于4时，利用Yule-Walker方法，按下式计算系统自适应参数α和
<mrow> <msub> <mi>r</mi> <mrow> <mi>k</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> <mrow> <mo>(</mo> <mn>1</mn> <mo>)</mo> </mrow> <mo>=</mo> <msub> <mi>r</mi> <mi>k</mi> </msub> <mrow> <mo>(</mo> <mn>1</mn> <mo>)</mo> </mrow> <mo>+</mo> <mfrac> <mn>1</mn> <mrow> <mi>k</mi> <mo>+</mo> <mn>1</mn> </mrow> </mfrac> <mo>&amp;lsqb;</mo> <mover> <mi>v</mi> <mo>^</mo> </mover> <mrow> <mo>(</mo> <mi>k</mi> <mo>+</mo> <mn>1</mn> <mo>)</mo> </mrow> <mover> <mi>v</mi> <mo>^</mo> </mover> <mrow> <mo>(</mo> <mi>k</mi> <mo>)</mo> </mrow> <mo>-</mo> <msub> <mi>r</mi> <mi>k</mi> </msub> <mrow> <mo>(</mo> <mn>1</mn> <mo>)</mo> </mrow> <mo>&amp;rsqb;</mo> </mrow>
<mrow> <msub> <mi>r</mi> <mrow> <mi>k</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> <mrow> <mo>(</mo> <mn>0</mn> <mo>)</mo> </mrow> <mo>=</mo> <msub> <mi>r</mi> <mi>k</mi> </msub> <mrow> <mo>(</mo> <mn>0</mn> <mo>)</mo> </mrow> <mo>+</mo> <mfrac> <mn>1</mn> <mrow> <mi>k</mi> <mo>+</mo> <mn>1</mn> </mrow> </mfrac> <mo>&amp;lsqb;</mo> <mover> <mi>v</mi> <mo>^</mo> </mover> <mrow> <mo>(</mo> <mi>k</mi> <mo>+</mo> <mn>1</mn> <mo>)</mo> </mrow> <mover> <mi>v</mi> <mo>^</mo> </mover> <mrow> <mo>(</mo> <mi>k</mi> <mo>+</mo> <mn>1</mn> <mo>)</mo> </mrow> <mo>-</mo> <msub> <mi>r</mi> <mi>k</mi> </msub> <mrow> <mo>(</mo> <mn>0</mn> <mo>)</mo> </mrow> <mo>&amp;rsqb;</mo> </mrow>
rk+1(1)为k+1时刻目标速度向前一步相关函数，rk(1)为k时刻目标速度向前一步相关函数，和分别为k时刻和k+1时刻目标速度估计值；rk+1(0)为k+1时刻目标速度自相关函数，rk(0)为k时刻目标速度自关函数；
根据系统方程得到目标运动速度，方程满足如下一阶马尔科夫随机序列：
<mrow> <mover> <mi>v</mi> <mo>^</mo> </mover> <mrow> <mo>(</mo> <mi>k</mi> <mo>+</mo> <mn>1</mn> <mo>)</mo> </mrow> <mo>=</mo> <mi>&amp;beta;</mi> <mover> <mi>v</mi> <mo>^</mo> </mover> <mrow> <mo>(</mo> <mi>k</mi> <mo>)</mo> </mrow> <mo>+</mo> <msup> <mi>w</mi> <mi>v</mi> </msup> <mrow> <mo>(</mo> <mi>k</mi> <mo>)</mo> </mrow> </mrow>
其中为k+1时刻的速度，为k时刻的速度，β为离散后速度随机序列的机动频率，wv(k)为零均值白噪声离散序列，方差为其中为零均值白噪声w(t)的方差，β与α的关系为β＝e-αT；x(k+1)为k+1时刻目标的状态向量，k为取样时刻，Φ(k+1,k)表示状态转移矩阵，x(k)为k时刻目标的状态向量，U(k)表示输入矩阵，为0时刻开始至k时刻目标的速度均值，w(k)表示过程噪音，其均值为0，方差为Q(k)，Φ(k+1,k)、U(k)和Q(k)中含有目标机动频率α和变化率方差随着系统自适应参数的变化而变化。
一阶马尔科夫时间加速度序列满足以下参数关系：
<mfenced open = "" close = ""> <mtable> <mtr> <mtd> <mrow> <mi>&amp;beta;</mi> <mo>=</mo> <mfrac> <mrow> <msub> <mi>r</mi> <mrow> <mi>k</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> <mrow> <mo>(</mo> <mn>1</mn> <mo>)</mo> </mrow> </mrow> <mrow> <msub> <mi>r</mi> <mrow> <mi>k</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> <mrow> <mo>(</mo> <mn>0</mn> <mo>)</mo> </mrow> </mrow> </mfrac> </mrow> </mtd> <mtd> <mrow> <msubsup> <mi>&amp;delta;</mi> <mrow> <mi>v</mi> <mi>w</mi> </mrow> <mn>2</mn> </msubsup> <mo>=</mo> <msub> <mi>r</mi> <mrow> <mi>k</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> <mrow> <mo>(</mo> <mn>0</mn> <mo>)</mo> </mrow> <mo>-</mo> <msub> <mi>&amp;alpha;r</mi> <mrow> <mi>k</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> <mrow> <mo>(</mo> <mn>1</mn> <mo>)</mo> </mrow> </mrow> </mtd> </mtr> </mtable> </mfenced>
其中，α与β分别为加速度的机动频率及其离散化后加速度序列的机动频率，自适应参数α和可按照下式计算得到：
<mrow> <mi>&amp;alpha;</mi> <mo>=</mo> <mo>-</mo> <mfrac> <mrow> <mi>ln</mi> <mi> </mi> <msub> <mi>r</mi> <mrow> <mi>k</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> <mrow> <mo>(</mo> <mn>1</mn> <mo>)</mo> </mrow> <mo>-</mo> <mi>ln</mi> <mi> </mi> <msub> <mi>r</mi> <mrow> <mi>k</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> <mrow> <mo>(</mo> <mn>0</mn> <mo>)</mo> </mrow> </mrow> <mi>T</mi> </mfrac> </mrow>
<mrow> <msubsup> <mi>&amp;delta;</mi> <mi>v</mi> <mn>2</mn> </msubsup> <mo>=</mo> <mfrac> <mrow> <msub> <mi>r</mi> <mrow> <mi>k</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> <mrow> <mo>(</mo> <mn>0</mn> <mo>)</mo> </mrow> <mo>-</mo> <msub> <mi>&amp;alpha;r</mi> <mrow> <mi>k</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> <mrow> <mo>(</mo> <mn>1</mn> <mo>)</mo> </mrow> </mrow> <mrow> <mn>1</mn> <mo>-</mo> <msup> <mrow> <mo>(</mo> <mfrac> <mrow> <msub> <mi>r</mi> <mrow> <mi>k</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> <mrow> <mo>(</mo> <mn>1</mn> <mo>)</mo> </mrow> </mrow> <mrow> <msub> <mi>r</mi> <mrow> <mi>k</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> <mrow> <mo>(</mo> <mn>0</mn> <mo>)</mo> </mrow> </mrow> </mfrac> <mo>)</mo> </mrow> <mn>2</mn> </msup> </mrow> </mfrac> </mrow>
其中，ln为取以e为底的对数计算；α和为系统自适应参数，T为采样间隔。
10.根据权利要求1所述的基于自适应滤波的GPS测量数据处理方法，其特征在于，步骤18的计算公式如下：
<mrow> <mi>&amp;Delta;</mi> <mi>R</mi> <mrow> <mo>(</mo> <mi>k</mi> <mo>+</mo> <mn>1</mn> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mn>1</mn> <mrow> <mi>N</mi> <mo>-</mo> <mn>1</mn> </mrow> </mfrac> <mrow> <mo>(</mo> <mi>e</mi> <mo>(</mo> <mrow> <mi>k</mi> <mo>+</mo> <mn>1</mn> </mrow> <mo>)</mo> <mo>-</mo> <mover> <mi>e</mi> <mo>&amp;OverBar;</mo> </mover> <mo>(</mo> <mrow> <mi>k</mi> <mo>+</mo> <mn>1</mn> </mrow> <mo>)</mo> <mo>)</mo> </mrow> <msup> <mrow> <mo>(</mo> <mi>e</mi> <mo>(</mo> <mrow> <mi>k</mi> <mo>+</mo> <mn>1</mn> </mrow> <mo>)</mo> <mo>-</mo> <mover> <mi>e</mi> <mo>&amp;OverBar;</mo> </mover> <mo>(</mo> <mrow> <mi>k</mi> <mo>+</mo> <mn>1</mn> </mrow> <mo>)</mo> <mo>)</mo> </mrow> <mi>T</mi> </msup> <mo>-</mo> <mfrac> <mn>1</mn> <mi>N</mi> </mfrac> <mrow> <mo>(</mo> <mi>H</mi> <mo>(</mo> <mrow> <mi>k</mi> <mo>+</mo> <mn>1</mn> </mrow> <mo>)</mo> <mi>P</mi> <mo>(</mo> <mrow> <mi>k</mi> <mo>+</mo> <mn>1</mn> <mo>|</mo> <mi>k</mi> </mrow> <mo>)</mo> <mi>H</mi> <msup> <mrow> <mo>(</mo> <mrow> <mi>k</mi> <mo>+</mo> <mn>1</mn> </mrow> <mo>)</mo> </mrow> <mi>T</mi> </msup> <mo>)</mo> </mrow> </mrow>
<mrow> <mi>R</mi> <mrow> <mo>(</mo> <mi>k</mi> <mo>+</mo> <mn>1</mn> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mrow> <mi>N</mi> <mo>-</mo> <mn>1</mn> </mrow> <mi>N</mi> </mfrac> <mi>R</mi> <mrow> <mo>(</mo> <mi>k</mi> <mo>)</mo> </mrow> <mo>+</mo> <mi>&amp;Delta;</mi> <mi>R</mi> <mrow> <mo>(</mo> <mi>k</mi> <mo>+</mo> <mn>1</mn> <mo>)</mo> </mrow> </mrow>
其中k为采样时刻，e(k+1)表示估计过程中的残差，表示k+1时刻的残差均值，H(k+1)表示k+1时刻的测量矩阵，P(k+1|k)表示向前一步估计方差，|表示条件操作符，H(k+1)表示k+1时刻的测量矩阵，N为残差的个数。
CN201710497032.0A 2017-06-26 2017-06-26 一种基于自适应滤波的gps测量数据处理方法 Active CN107229060B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201710497032.0A CN107229060B (zh) 2017-06-26 2017-06-26 一种基于自适应滤波的gps测量数据处理方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201710497032.0A CN107229060B (zh) 2017-06-26 2017-06-26 一种基于自适应滤波的gps测量数据处理方法

Publications (2)

Publication Number Publication Date
CN107229060A true CN107229060A (zh) 2017-10-03
CN107229060B CN107229060B (zh) 2019-12-03

Family

ID=59935504

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201710497032.0A Active CN107229060B (zh) 2017-06-26 2017-06-26 一种基于自适应滤波的gps测量数据处理方法

Country Status (1)

Country Link
CN (1) CN107229060B (zh)

Cited By (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108646277A (zh) * 2018-05-03 2018-10-12 山东省计算中心（国家超级计算济南中心） 基于抗差自适应与扩展卡尔曼滤波融合的北斗导航方法
CN109147390A (zh) * 2018-08-20 2019-01-04 浙江工业大学 一种基于量化自适应卡尔曼滤波的车辆轨迹跟踪方法
CN109240171A (zh) * 2018-10-25 2019-01-18 中核新科（天津） 精密机械制造有限公司 用于表面处理工艺的工作液主盐浓度监测及补偿装置和方法
CN110032711A (zh) * 2019-04-22 2019-07-19 中南大学 一种基于参数动态调节的快速卡尔曼在线检测滤波方法
CN110069748A (zh) * 2019-05-05 2019-07-30 北京无线电测量研究所 实时确定数据方差方法及装置
CN110850363A (zh) * 2019-10-22 2020-02-28 南京大学 一种基于实时定位轨迹数据进行动态滤波优化的方法
CN111683134A (zh) * 2020-06-04 2020-09-18 勇鸿(重庆)信息科技有限公司 基于区块链技术的分布式车联网数据传输系统及方法
CN115859007A (zh) * 2023-02-17 2023-03-28 广东石油化工学院 一种石化仪表采样数据滑窗约束容错滤波降噪方法和装置
CN115980803A (zh) * 2023-03-17 2023-04-18 北京航空航天大学 基于双频码伪距和载波相位观测量进行伪距平滑的方法

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104112079A (zh) * 2014-07-29 2014-10-22 洛阳理工学院 一种模糊自适应变分贝叶斯无迹卡尔曼滤波方法
CN105549049A (zh) * 2015-12-04 2016-05-04 西北农林科技大学 一种应用于gps导航的自适应卡尔曼滤波算法

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104112079A (zh) * 2014-07-29 2014-10-22 洛阳理工学院 一种模糊自适应变分贝叶斯无迹卡尔曼滤波方法
CN105549049A (zh) * 2015-12-04 2016-05-04 西北农林科技大学 一种应用于gps导航的自适应卡尔曼滤波算法

Non-Patent Citations (4)

* Cited by examiner, † Cited by third party
Title
JIN XUE-BO ET AL.: "Maneuvering target tracking by adaptive statistics model", 《SCIENCEDIRECT》 *
YI SHENG-LUN ET AL.: "An Improved Online Denoising Algorithm Based on Kalman Filter and Adaptive Current Statistics Model", 《RESEARCHGATE》 *
李理敏 等: "基于自适应扩展卡尔曼滤波的载波跟踪算法", 《航空学报》 *
胡洪涛 等: "基于’当前’ 统计模型的模糊自适应跟踪算法", 《系统仿真学报》 *

Cited By (11)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108646277A (zh) * 2018-05-03 2018-10-12 山东省计算中心（国家超级计算济南中心） 基于抗差自适应与扩展卡尔曼滤波融合的北斗导航方法
CN109147390A (zh) * 2018-08-20 2019-01-04 浙江工业大学 一种基于量化自适应卡尔曼滤波的车辆轨迹跟踪方法
CN109147390B (zh) * 2018-08-20 2020-06-02 浙江工业大学 一种基于量化自适应卡尔曼滤波的车辆轨迹跟踪方法
CN109240171A (zh) * 2018-10-25 2019-01-18 中核新科（天津） 精密机械制造有限公司 用于表面处理工艺的工作液主盐浓度监测及补偿装置和方法
CN110032711A (zh) * 2019-04-22 2019-07-19 中南大学 一种基于参数动态调节的快速卡尔曼在线检测滤波方法
CN110032711B (zh) * 2019-04-22 2022-07-12 中南大学 一种基于参数动态调节的快速卡尔曼滤波的在线检测方法
CN110069748A (zh) * 2019-05-05 2019-07-30 北京无线电测量研究所 实时确定数据方差方法及装置
CN110850363A (zh) * 2019-10-22 2020-02-28 南京大学 一种基于实时定位轨迹数据进行动态滤波优化的方法
CN111683134A (zh) * 2020-06-04 2020-09-18 勇鸿(重庆)信息科技有限公司 基于区块链技术的分布式车联网数据传输系统及方法
CN115859007A (zh) * 2023-02-17 2023-03-28 广东石油化工学院 一种石化仪表采样数据滑窗约束容错滤波降噪方法和装置
CN115980803A (zh) * 2023-03-17 2023-04-18 北京航空航天大学 基于双频码伪距和载波相位观测量进行伪距平滑的方法

Also Published As

Publication number Publication date
CN107229060B (zh) 2019-12-03

Similar Documents

Publication Publication Date Title
CN107229060B (zh) 一种基于自适应滤波的gps测量数据处理方法
CN108682023A (zh) 基于Elman神经网络的紧耦合无迹卡尔曼跟踪滤波算法
CN108759833B (zh) 一种基于先验地图的智能车辆定位方法
CN104567871B (zh) 一种基于地磁梯度张量的四元数卡尔曼滤波姿态估计方法
CN109211276A (zh) 基于gpr与改进的srckf的sins初始对准方法
CN111985093A (zh) 一种带噪声估计器的自适应无迹卡尔曼滤波状态估计方法
CN104121907B (zh) 一种基于平方根容积卡尔曼滤波器的飞行器姿态估计方法
CN111156987A (zh) 基于残差补偿多速率ckf的惯性/天文组合导航方法
CN114609912A (zh) 基于伪线性最大相关熵卡尔曼滤波的仅测角目标追踪方法
CN107169478A (zh) 一种自适应在线滤波方法
CN103308896A (zh) 一种适于非引擎机动目标的高精度跟踪方法
CN111983927A (zh) 一种新型的最大协熵椭球集员滤波方法
CN111291319B (zh) 一种应用于非高斯噪声环境下的移动机器人状态估计方法
CN102706345A (zh) 一种基于衰减记忆序贯检测器的机动目标跟踪方法
CN110989341B (zh) 一种约束辅助粒子滤波方法及目标跟踪方法
Zhao et al. A filtering approach based on MMAE for a SINS/CNS integrated navigation system
CN111189454A (zh) 基于秩卡尔曼滤波的无人车slam导航方法
CN117419723A (zh) 一种基于因子图的交互模型自适应组合导航方法
CN114445459B (zh) 基于变分贝叶斯理论的连续-离散最大相关熵目标跟踪方法
CN110007298B (zh) 一种目标超前预测跟踪方法
CN114061592B (zh) 基于多模型的自适应鲁棒auv导航方法
CN115828533A (zh) 一种基于Student’s t分布的交互多模型鲁棒滤波方法
Zarei et al. Performance improvement for mobile robot position determination using cubature Kalman filter
CN111998854B (zh) 基于Cholesky分解计算的精确扩展Stirling插值滤波方法
CN114637956A (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