一种含未知输入和非高斯量测噪声的随机系统滤波器

摘要

本发明涉及一种含未知输入和非高斯量测噪声的随机系统滤波器。首先，针对含未知输入/干扰的离散时间线性随机系统，基于状态方程完成状态的预测估计；其次，针对实际系统的量测噪声通常含有较多野值，噪声的概率分布往往具有尖峰肥尾等强非高斯特性，在求解干扰估计与状态滤波时基于Huber函数代替传统的均方误差准则构造指标函数；再次，通过最小化指标函数，基于不动点迭代法解算出干扰估计与状态的滤波估计；最后，将状态预测、干扰估计、状态滤波按时刻递推，给出滤波器的设计流程。本发明可推广应用于组合导航、目标跟踪、信号处理等领域，解决实际系统的高精度状态估计问题。

权利要求

1.一种含未知输入和非高斯量测噪声的随机系统滤波器，其特征在于，包括以下步骤：第一步，基于系统状态方程进行状态预测，得到状态的预测估计与相应的协方差矩阵；

第二步，基于Huber函数构造干扰与状态联合估计的性能指标函数；

第三步，通过最小化第二步中的性能指标函数，基于矩阵运算和不动点迭代法解算出干扰估计与状态的滤波估计，以及相应的协方差矩阵；

第四步，将前述步骤中的状态预测、干扰估计、状态滤波按时刻递推，最终得到滤波器；

所述第一步，状态预测过程如下：(1)一类含有未知输入/干扰的离散时间线性随机系统，其状态空间模型表示如下：其中， 为系统状态变量， 为量测输出变量， 为系统的未知输入，即动态特性未知的干扰， 和 分别表示过程噪声和量测噪声，为互不相关的零均值白噪声且方差为正定阵Qk和Rk，过程噪声ωk认为是高斯分布，而量测噪声υk具有强非高斯特性，Ak，Gk，Ck和Hk为已知的时变矩阵，(Ak,Ck)可观且Hk列满秩，系统初始状态x0服从高斯分布 且与噪声均不相关；

(2)基于状态方程，根据前一时刻得到的干扰和状态的估计，得到当前时刻状态的预测估计、预测误差与协方差矩阵分别如下三个公式所示：在初始k＝0时刻，状态预测直接由系统的初始条件给定，即所述第二步，基于Huber函数的性能指标函数构造过程如下：基于量测噪声的非高斯特性，使用Huber函数构造相应的性能指标函数，Huber函数等价于残差的l1范数，对于噪声野值具有抑制能力，定义 其第i个分量表示为ξk[i]，构造出如下的性能指标函数：其中，ρ(·)为Huber函数，定义如下：其中τ是函数自变量，γ为一个可调节的参数；

所述第三步，干扰估计与状态滤波过程如下：(1)基于第二步构造的性能指标函数，将其最小化可以得到干扰的估计与状态的滤波估计，即需要求解如下方程：

定义函数ψ(τ)＝ρ′(τ)/τ，同时定义矩阵Ψk＝diag{ψ(ξk[1]),…,ψ(ξk[m])},经过相

关矩阵运算，方程的解表示为：

(2)Mk、Kk均与Ψk有关，Ψk是xk和dk的函数，上式是关于xk和dk的不动点方程，基于不动点迭代法，得到下述迭代过程：

其中，t为迭代时标， 和 由前述定义结合 给出，迭代初值可以取为 迭代终止条件可直接设定迭代次数或者给出判断准则及阈值；迭代终止后，记 其中tf是最后一次迭代的时标，则干扰估计与状态滤波表示为：同时，相应的估计误差与协方差矩阵分别如下两式所示：其中 是状态滤波误差与干扰估计误差之间的协方差阵， 是 的转置， 是干扰估计误差的协方差阵。

2.根据权利要求1所述的含未知输入和非高斯量测噪声的随机系统滤波器，其特征在于：所述第四步，滤波器的按时刻递推过程如下：在每一时刻首先进行状态预测获取状态的预测估计 及相应的协方差 接着进行干扰估计与状态滤波获取 及相应的协方差 进而递推至下一时刻，重复这些步骤，最终即得到了滤波器的完整工作流程。

说明书

一种含未知输入和非高斯量测噪声的随机系统滤波器

技术领域

[0001] 本发明涉及一种含未知输入/干扰和非高斯量测噪声的随机系统滤波器，针对系统所受到的未知输入或未知动态干扰，采用状态预测、干扰估计、状态滤波的递推滤波器结
构，同时考虑到量测噪声的非高斯特性，基于Huber函数构造干扰与状态估计的性能指标函
数，增强了滤波器对于量测野值的鲁棒性，本发明可应用于飞行器、船舶、车辆等运动体的
组合导航系统，也可推广应用于目标跟踪、信号处理等相关领域。

背景技术

[0002] 飞行器、船舶、车辆等运动体依靠导航系统实时获取自身的运动信息和姿态信息，而滤波算法又是导航系统信息处理的核心所在。传统的以卡尔曼滤波为典型代表的滤波算
法在实际工程应用中已经取得了很大的成功，发挥了积极的作用。然而，运动体的任务需求
日益提高，所面临的环境也日益复杂多样，这些都对导航系统的实时性、可靠性和精确性提
出了进一步的挑战。为了提高导航系统的性能，以便能够快速、可靠、精准地获取和估计运
动信息和姿态信息，一方面需要从硬件角度提高传感器性能或者研制新的智能感知机制，
另一方面则要从软件角度提升滤波算法的适应性和鲁棒性。特别是在硬件资源有限的条件
下，后者更是具有重要意义。

[0003] 实际系统总是或多或少受到未知输入或未知干扰的影响。以惯性导航系统为例，陀螺仪的漂移、加速度计的常值偏置、外部环境的振动、安装误差等都可以视为干扰。因此，
研究存在未知输入/干扰下的滤波方法即抗干扰滤波，具有积极的理论意义和工程应用价
值。作为一种典型的抗干扰滤波方法，未知输入滤波无需对干扰信号的特性做出假设，能够
在干扰信号完全未知的情形下给出状态的最优估计，一般为最小方差无偏估计。未知输入
滤波最初由美国学者Kitanidis于上世纪八十年代提出，经过三十余年的发展，已建立起了
一套较为成熟的理论体系。同时，这些理论成果也成功应用在了导航定位、环境监测、故障
检测等相关领域。例如，专利1(火星动力下降段非线性三步滤波方法，申请号为
201310341821.7)利用未知输入滤波技术对火星探测器动力下降段的位置和速度进行精确
估计；专利2(一种基于智能手机移动检测的城市快速路交通状态估计方法，申请号为
201610056587.7)将未知输入滤波应用于城市快速路交通状态估计问题当中，实现交通密
度和边界流量的同步估计。

[0004] 然而，现有的研究成果都是基于系统的过程噪声和量测噪声均服从高斯分布的假设，在最小方差意义下给出状态的最优估计。事实上，非高斯噪声广泛存在于实际系统当
中。当随机变量的概率分布具有偏斜非对称、尖峰肥尾等强非高斯特性时，一般不能够再近
似为高斯分布，见专利3(一种基于负熵的非高斯线性随机系统滤波方法，申请号为
201410124598.5)。特别是对于量测噪声，其概率分布经常出现肥尾现象，这意味着量测噪
声取得极端异常值的可能性很大。或者说，系统的量测通道经常受到野值的影响。在这种情
况下，现有滤波器的估计性能会变差，甚至有可能导致滤波发散。因此，有必要研究一种针
对非高斯随机系统的未知输入滤波器。

发明内容

[0005] 本发明的技术解决问题是：针对现有的未知输入滤波方法局限于高斯系统的不足，提供一种适合于非高斯系统的新型随机抗干扰滤波器，解决复杂系统的高精度状态估
计问题。

[0006] 本发明的技术解决方案为：一种含未知输入和非高斯量测噪声的随机系统新型抗干扰滤波器，其特征在于，包括以下步骤：首先，基于状态方程得到状态的预测估计与预测
误差的协方差矩阵；其次，基于Huber函数构造干扰与状态联合估计的性能指标函数；进而，
通过最小化性能指标函数，借助矩阵相关运算和不动点迭代法给出干扰估计与状态的滤波
估计，以及相应的估计误差和协方差矩阵；最后，总结归纳出包括状态预测、干扰估计、状态
滤波在内的递推滤波器设计流程；具体步骤如下：

[0007] 第一步，状态预测

[0008] (1)一类含有未知输入/干扰的离散时间线性随机系统，其状态空间模型表示如下：

[0009]

[0010] 其中， 为系统状态变量， 为量测输出变量， 为系统的未知输入，也可认为是动态特性未知的干扰， 和 分别表示过程噪声和量测噪
声，为互不相关的零均值白噪声且协方差分别为正定矩阵Qk和Rk，过程噪声ωk认为服从高
斯分布，而量测噪声υk为非高斯随机变量，尤其是具有尖峰肥尾的强非高斯特性，Ak，Gk，Ck
和Hk为已知的时变矩阵，假设(Ak,Ck)可观且Hk列满秩，系统初始状态x0服从高斯分布
且与噪声ωk和υk均不相关。

[0011] (2)基于状态方程，根据系统运行的前一时刻的干扰估计和状态滤波估计，可得到当前时刻的状态预测估计、预测误差、预测误差的协方差矩阵，分别如下三式所示：

[0012]

[0013]

[0014]

[0015] 需要指出的是，在k＝0时刻，状态预测不用以上公式，直接由初始条件给定，即

[0016] 第二步，基于Huber函数的性能指标函数构造

[0017] 考虑到量测噪声的尖峰肥尾非高斯特性，这里采用Huber函数构造相应的性能指标函数。当量测噪声含有野值时，Huber函数等价于残差的l1范数，因此对量测野值具有较
强的鲁棒性。定义 其第i个分量表示为ξk[i]。给出如下的性能指
标函数：

[0018]

[0019] 其中，ρ(·)为Huber函数，定义如下：

[0020]

[0021] 其中τ是函数自变量，γ为一个可调节的参数。

[0022] 第三步，干扰估计与状态滤波

[0023] (1)基于第二步构造的性能指标函数，将其最小化则可得到干扰的估计与状态的滤波估计，即估计问题转化为优化问题 为了求解这一优化问题，将指
标函数分别对xk和dk求导，可以得到以下方程组：

[0024]

[0025] 进一步，定义函数ψ(·)为ψ(τ)＝ρ′(τ)/τ，同时定义矩阵Ψk＝diag{ψ(ξk[1]),…,ψ(ξk[m])}， 因此，可以得到如下所示的方程组：

[0026]

[0027] 根据分块矩阵求逆公式、矩阵逆引理以及代数运算，可以求解上述方程组。为了表述简便，进一步定义矩阵
则方程组的解表示如下：

[0028]

[0029] (2)注意到在上述结果中，Mk，Kk均与Ψk有关，而Ψk是xk和dk的函数，故上式可看作是关于xk和dk的不动点方程。基于不动点迭代法，我们可以得到下述的迭代求解过程：

[0030]

[0031] 其中，t表示迭代时标， 和 由前述定义结合给出，迭代初值可以取为 迭代终止条件或者可以通过直接设定迭代
次数给出，或者可以给出相应判断准则及阈值。迭代终止后，记
其中tf表示最后一次迭代的时标。这样，
干扰的估计与状态的滤波估计表示为：

[0032]

[0033] 同时，我们可以得到相应的估计误差与协方差矩阵如下两式所示：

[0034]

[0035]

[0036] 第四步，滤波器按时刻递推

[0037] 基于前三个步骤的结果，在系统运行的每一时刻首先进行状态预测获取接着进行干扰估计与状态滤波获取 进而递
推至系统的下一时刻，重复这些步骤，即可得到滤波器的完整工作流程。

[0038] 本发明与现有技术相比的优点在于：针对同时含有未知输入/干扰和非高斯量测噪声的复杂随机系统，提出了一种新型随机抗干扰滤波方法，通过采用基于Huber函数的未
知输入滤波算法，相对于现有的基于最小方差准则的未知输入滤波方法而言，提高了滤波
器对量测野值的鲁棒性，增强了滤波器对受扰非高斯随机系统的适应性，使得滤波器能够
在复杂情况下实现高精度的状态估计。

附图说明

[0039] 图1为本发明一种含未知输入和非高斯量测噪声的随机系统滤波器的设计流程图；

[0040] 图2为本发明一种含未知输入和非高斯量测噪声的随机系统滤波器的不动点迭代法过程示意图。

具体实施方式

[0041] 如图1所示，本发明具体实现步骤如下(以某一做近似匀加速直线运动的运动体为例来说明方法的具体实现)：

[0042] 1、状态预测

[0043] (1)当运动体做近似匀加速直线运动时，运动学方程如下所示：

[0044]

[0045] 量测方程为：

[0046] yk＝sk+dk+υk

[0047] 其中sk，vk和ak分别表示运动体在k时刻的位置、速度和加速度，量测信息yk由位置传感器提供，T为采样周期，dk表示未知输入/干扰， (i＝1,2,3)为相互独立
的高斯白噪声，υk表示零均值方差为Rk的噪声，其具有学生t分布、Laplace分布或者其他尖
峰肥尾分布的特性。将运动学方程和量测方程写成紧凑形式，则有：

[0048]

[0049] 其中，相关矩阵和向量定义如下：

[0050]

[0051] Ck＝[1 0 0]，Hk＝1

[0052] 显然，(Ak,Ck)满足可观性条件，Hk满足列满秩要求，ωk为零均值高斯向量方差阵为 初始状态 且与噪声均不相关。

[0053] (2)基于状态方程，根据运动体前一时刻的干扰估计和状态滤波，可得到当前时刻的状态预测、预测误差、预测误差的协方差矩阵分别如下三式所示：

[0054]

[0055]

[0056]

[0057] 需要指出的是，在k＝0时刻，状态预测不用以上公式，直接由初始条件给定，即

[0058] 2、基于Huber函数的性能指标函数构造

[0059] 考虑到量测噪声的尖峰肥尾非高斯特性，使用Huber函数构造相应的性能指标函数。当量测噪声含有野值时，Huber函数等价于残差的l1范数，因此对量测野值具有较强的
鲁棒性。定义 其第i个分量表示为ξk[i]。我们给出如下的性能指
标函数：

[0060]

[0061] 其中，ρ(·)为Huber函数，定义如下：

[0062]

[0063] 其中τ是函数自变量，γ为一个可调节的参数。

[0064] 3、干扰估计与状态滤波

[0065] (1)基于上一步构造的性能指标函数，将其最小化则可得到干扰估计与状态滤波，即需要求解优化问题： 将指标函数Jk对xk和dk分别求导，可以得到以下
方程组：

[0066]

[0067] 定义函数ψ(τ)＝ρ′(τ)/τ，同时定义矩阵Ψk＝diag{ψ(ξk[1]),…,ψ(ξk[m])}，则上述方程组具有如下形式：

[0068]

[0069] 根据分块矩阵求逆公式、矩阵逆引理以及相关代数运算，我们可以求解上述方程组。为了表述简便，进一步定义矩阵
则方程组的解可表示为：

[0070]

[0071] (2)在上述结果中，注意到Mk，Kk均与Ψk有关，而Ψk是xk和dk的函数，故上式可看作是关于xk和dk的不动点方程。基于不动点迭代法，我们可以得到下述迭代过程：

[0072]

[0073] 其中，t为迭代时标， 和 由前述定义结合给出，迭代初值可以取为 迭代终止条件可以通过直接设定迭代次数
给 出 ，或 者 也 可以 给 出 相 应 判断 准 则 及阈 值 。迭 代 终 止后 ，我 们 记
其中tf是最后一次迭代的时标，则干扰的
估计和状态的滤波估计表示为：

[0074]

[0075] 同时，可以计算得到干扰估计误差、状态滤波误差，以及相应的协方差矩阵，如下所示：

[0076]

[0077]

[0078] 4、滤波器按时刻递推

[0079] 基于前述步骤中的结果，在系统运行的每一时刻首先进行状态预测获取接着进行干扰估计与状态滤波获取 进而递
推至下一时刻，重复这些步骤，可得到滤波器的完整工作流程，即可给出复杂环境下运动体
的运动信息(包括位置、速度和加速度)的实时高精度估计。

[0080] 本发明说明书中未作详细描述的内容属于本领域专业技术人员公知的现有技术。