【技术领域】
本发明提供一种多模型自校准秩滤波方法，属于非高斯鲁棒滤波技术领域。
背景技术：
针对工程中常见的系统状态方程受未知输入影响的问题，文献“自校准kalman滤波方法[j].航空动力学报.2014,29(06):1363-1368”提出了一种自校准卡尔曼滤波方法(self-calibrationkalmanfilter,skf)，该方法在依照原始状态方程进行迭代运算的同时，对未知输入项进行估计，从而使未知输入的影响自动得到补偿。在此基础上，研究人员又先后发展了自校准扩展卡尔曼滤波方法(self-calibrationextendedkalmanfilter,sekf)和自校准无迹卡尔曼滤波方法(self-calibrationunscentedkalmanfilter,sukf)，它们将自校准技术推广到了非线性领域。
但是由于系统不确定因素的存在，未知输入也有为零的可能。在这种情况下，自校准卡尔曼滤波方法由于在先验估计中引入了对未知输入项的估计，尽管该估计值很小，其滤波精度仍不及没有考虑未知输入影响的标准卡尔曼滤波方法。为了进一步提升自校准卡尔曼滤波在未知输入为零段的滤波精度，研究人员基于多模型估计理论，提出了多模型自校准卡尔曼滤波方法(multiple-modelself-calibrationkalmanfilter，mskf)。
特别地，对于非线性系统，多模型自校准扩展卡尔曼滤波(multiple-modelself-calibrationextendedkalmanfilter，msekf)和多模型自校准无迹卡尔曼滤波方法(multiple-modelself-calibrationunscentedkalmanfilter,msukf)也被开发出来，但是它们只适用于服从高斯分布的系统，对于非高斯分布的非线性系统未知输入存在为零段与非零段的情况便不再适用。
技术实现要素：
本发明的目的是提供一种多模型自校准秩滤波方法(multiple-modelself-calibrationrankfilter,msrf)，它通过将多模型估计理论引入到自校准秩滤波方法(self-calibrationrankfilter,srf)中，将mskf的应用范围拓展到了非高斯、非线性领域。其同时采用rf与srf进行计算，实时更新二者权重，进而得到状态估计。
本发明一种多模型自校准秩滤波方法，它包含以下六个步骤：
步骤一：建立系统基本方程
多模型自校准秩滤波采用自校准秩滤波与秩滤波两种方法进行运算，故系统包含两个状态方程，第一个为含有未知输入项的状态方程，第二个为标准的非线性状态方程，其具体表达式为
zk＝hk(xk)+vk(3)
式中，xk表示系统的状态向量，
式中，cov[·]为协方差，e[·]为数学期望，δkj为δ函数，当k＝j时，δkj＝1，当k≠j时，δkj＝0；
步骤二：对系统进行滤波初始化
设定状态估计与估计误差方差矩阵的初始值为
同时，为了完成两模型估计结果的融合，还需要设定两种模型的概率初始值
pr(1|z3)＝pr(2|z3)＝0.5(7)
以及用于迭代计算的概率初始值prmax和prmin；初始化prmax和prmin的原因如下：
在多模型估计运算过程中，一些模型会由于对应的概率逐渐趋近为零而被淘汰，故参与运算的模型数量n在不断减小，这会降低系统对复杂环境的适应能力；针对本文方法只选取两个动力学模型且采用最高概率法选取先验估计值，只需要定性分析两种模型的概率大小而不需要精确计算概率值的特点，多模型自校准卡尔曼滤波方法不再使用计算得到的条件概率值进行迭代，而是设定两个确定的概率初始值prmin和prmax＝1-prmin；在每一步滤波之前，通过比较上一步概率结果的大小将其分别赋给两个模型，并以它们为初始值更新当前时刻的模型概率；由于prmin并不是概率下限那样的极小值，因此可以保证概率恢复的速度，从而使卡尔曼滤波的实时性得到保证；
步骤三：对系统进行时间更新
设k-1时刻的状态估计值和误差方差矩阵分别为
当k＝1,2时
状态一步预测值
xk/k-1,i＝f(χk-1,i)(9)
式中，n为状态向量xk的维度，标准正态偏量
一步预测误差方差矩阵
式中，ω为协方差权重系数，其计算公式为
当k＞2时
状态一步预测值
一步预测误差方差矩阵
式中
j＝argmaxjpr(j|zk)(15)
其中，函数argmax[f(x)]返回当f(x)最大时x的值，并且
在上述计算过程中，式(18)和式(21)提供自校准状态方程的一步预测结果，式(20)和式(21)提供标准状态方程计算的一步预测结果，然后通过式(15)中概率大小的比对，由式(13)完成最终的一步预测值筛选；概率大小的计算依托于贝叶斯原理，比较的是在当前时刻量测值已知的情况下，两种模型的条件概率，如式(23)和式(24)所示；
步骤四：迭代变量更新
在步骤三中有很多中间变量需要实时更新，因此有必要得到它们的递推公式，进而保证整个滤波过程的顺利进行；这些迭代变量包括：
各模型量测更新
式中
各模型条件概率重置
pr(j|zk)＝prmax(32)
pr[(3-j)|zk]＝prmin(33)；
步骤五：量测更新
首先，基于状态一步预测
然后，根据{χk/k-1,i}计算量测估计值
zk/k-1,i＝h(χk/k-1,i)(35)
由式(35)和式(36)计算量测误差方差矩阵pzz
由式(34)、式(35)和式(36)计算误差协方差矩阵pxz
由式(37)和式(38)计算增益矩阵
则，由量测更新得到最终的状态估计值和状态估计误差方差矩阵
步骤六：迭代计算
根据k时刻的状态估计值
通过步骤一到步骤六，本发明充分利用了秩滤波和自校准秩滤波两种方法的计算结果，依托基于贝叶斯原理的多模型估计理论，自动区分未知输入为零段与非零段，从而可以精确地选择其中最合适的结果作为自身的先验估计。最重要的一点在于，本发明是针对非高斯、非线性系统开发的，相较于其他方法适用范围更广。
本发明的有益效果包括
(1)将多模型估计理论引入到自校准秩滤波方法中，推导得到了多模型自校准秩滤波方法的完整过程。
(2)可以自动识别未知输入为零段与非零段，分别发挥秩滤波方法与自校准秩滤波方法各自的优势。
(3)进一步提高了系统受未知输入影响时的滤波精度，同时提高了滤波过程的稳定性，增强了系统的鲁棒性。
(4)可以处理不服从高斯分布的非线性系统，适用范围得以扩展，符合工程应用需求。
【附图说明】
图1为本发明方法流程示意图。
图2为步骤三时间更新流程示意图。
发明中序号、符号、代号说明如下：
xk：系统的状态向量
zk：系统量测向量
f(·)、h(·)：非高斯、非线性向量函数
bk：未知输入
wk：系统噪声向量
vk：量测噪声向量
qk：系统噪声向量方差矩阵
rk：量测噪声向量方差矩阵
cov[·]：协方差
e[·]：数学期望
δkj：δ函数
{χi}：秩采样点集
pr(1|zk)、pr(2|zk)：迭代计算的概率初始值
prmax、prmin：用于迭代计算的概率初始值
pk/k-1：一步预测误差方差矩阵
ω：协方差权重系数
argmax[f(x)]：函数最大值返回函数
pzz:量测误差方差矩阵
pxz:误差协方差矩阵
【具体实施方式】
下面结合附图对本发明作详细说明。
本发明提出了一种多模型自校准秩滤波方法，其流程图如图1所示，时间更新流程图如图2所示，它包括以下六个步骤：
步骤一：建立系统基本方程
式中，xk表示系统的状态向量，
式中，cov[·]为协方差，e[·]为数学期望，δkj为δ函数，当k＝j时，δkj＝1，当k≠j时，δkj＝0。
步骤二：滤波初始化
设定状态估计与估计误差方差矩阵的初始值
设定两种模型的概率初始值
pr(1|z3)＝pr(2|z3)＝0.5(47)
以及用于迭代计算的概率初始值prmax和prmin。
步骤三：时间更新
当k＝1,2时
状态一步预测值
xk/k-1,i＝f(χk-1,i)(49)
式中，n为状态向量xk的维度，标准正态偏量
一步预测误差方差矩阵
式中，ω为协方差权重系数，其计算公式为
当k＞2时
状态一步预测值
一步预测误差方差矩阵
式中
j＝argmaxjpr(j|zk)(55)
其中，函数argmax[f(x)]返回当f(x)最大时x的值，并且
其时间更新流程示意图见图2所示；
步骤四：迭代变量更新
各模型量测更新
式中
各模型条件概率重置
pr(j|zk)＝prmax(72)
pr[(3-j)|zk]＝prmin(73)
步骤五：量测更新
进行重采样，得到新的秩采样点集{χk/k-1,i}
然后，根据{χk/k-1,i}计算量测估计值
zk/k-1,i＝h(χk/k-1,i)(75)
计算量测误差方差矩阵pzz
计算误差协方差矩阵pxz
计算增益矩阵
状态估计值
步骤六：迭代计算
k＝k+1(83)
重复步骤三、步骤四和步骤五，直至滤波过程结束。