![航天器多源信息融合自主导航技术](https://wfqqreader-1252317822.image.myqcloud.com/cover/473/32855473/b_32855473.jpg)
2.4 动态系统中的状态估计算法
第2.2和2.3节介绍的是一般的参数估计问题,而航天器自主导航系统中要估计的是动态系统中的状态,其随时间演化的方程由轨道动力学给出,因此在本节中对动态系统的估计算法进行介绍。
2.4.1 递归贝叶斯估计算法
图2-1给出了一阶隐式马尔科夫模型,在该模型中,状态量x不能直接观测到;z为x的函数,且能够被观测到,称为观测量。
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0080_0001.jpg?sign=1739341055-Kl5y4xTJ3Cv4hb5mgqC76Ts88TkseUfr-0-2f52b50ad55add1651403884cc42c5a0)
图2-1 一阶隐式马尔科夫模型
一阶隐式马尔科夫模型具有如下属性:
属性1:在给定tk-1时刻的状态xk-1的前提下,tk时刻状态xk(以及未来时刻的状态xk+1,xk+2,…)出现的概率和tk-1以前时刻的状态、观测量条件独立,也即
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0080_0002.jpg?sign=1739341055-drlwRTbRVYoNrT2pi2YdYDev94dKVnhZ-0-414066a6932c1b2918d9d623b0873677)
式中,。
同样,在给定当前时刻状态的前提下,过去的状态和未来的状态、观测量独立,也即
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0081_0001.jpg?sign=1739341055-5pgUjjlwKC0oosWuZNn6UdaV1UD7R0pe-0-4dd43138c6027e19e79cd8157eb5d26c)
属性2:当前时刻的观测量仅依赖于当前时刻的状态量,和其他时刻的状态条件独立,也即
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0081_0002.jpg?sign=1739341055-tLnvbfNDJCS58Cz8HChkzpK17bi6TVuq-0-830420616b9352c3ff6c60238b4e987e)
定理2.4:递归贝叶斯估计算法给出在获得了tk时刻以及tk时刻以前的所有观测数据的前提下,递归地计算的方法。具体如下:
(1)初始化:给出初始时刻的状态先验概率密度函数p(x0)。
(2)预测:给定动力学模型,在tk时刻xk预测的概率密度函数由Chapman-Kolmogorov公式计算得到,即
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0081_0004.jpg?sign=1739341055-g8zn3rvIRxkRoo7Fh7MA8zpsq1Jt9X5a-0-b184a988ae95fb3b0d3e7009e651d1cf)
(3)更新:在tk时刻,给定观测量zk,则xk的后验概率密度函数可以由贝叶斯法则给出,即
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0081_0005.jpg?sign=1739341055-lX6HvOA97nF5IsVQ1B6liL4bueuBCOQz-0-628c058774d00ba04d91b79e54049d0b)
式中,,为归一化的常数。
证明:
首先由贝叶斯法则可以得到
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0081_0007.jpg?sign=1739341055-Q3hQKEBZQUuz6JTvcrqDic6RAsYAiUAd-0-6361b649c3ff94b42a44fa2992d83b30)
式(2-63)右边的分子可以写作
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0081_0008.jpg?sign=1739341055-Nkniz6S73GI0yT8mepWBXi6wJjUwwJJJ-0-9ada3abff15180f4fb8c7d7eff5c632d)
式(2-64)利用到了如下公式
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0081_0009.jpg?sign=1739341055-81MKodS0UsEmw0foo9avkYSopyJWYFv0-0-40cb99b2a2e78651a476d2b6e31f489f)
利用马尔科夫属性2,式(2-64)可以写作
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0081_0010.jpg?sign=1739341055-cWEgyk6z82T7seRUXeulaGSTnixvGwLs-0-f6e23b63cf584e86ae7a94f0819d0f5c)
式(2-63)的分母可以写作
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0081_0011.jpg?sign=1739341055-7xwKPdYBzsZBQV8u1AGCu5aefkfE7iUD-0-de028b0bfbbbc7cec7d7781a9e94bd03)
将式(2-65)和式(2-66)代入式(2-63)可以得到
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0082_0001.jpg?sign=1739341055-FHQpLKzPHw0M8kwAxvCf1qi4ikaNSToA-0-4eef83de787ae29229df3800394b9435)
在式(2-67)中为归一化的系数,
和测量方程对应,为更新过程;
由下式给出
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0082_0005.jpg?sign=1739341055-nkPQew3lonTCIG5f9W1av32ORuP57igT-0-5e4f92f9a71430f1af204d5c3b0a81db)
式中和状态预测过程相对应。在推导该式的过程中,利用了马尔科夫属性1 [式(2-58)]。
式(2-67)给出了马尔科夫动态系统贝叶斯估计的一般算法,观察式(2-67)和式(2-68)可以知道为
的函数,因此可以用定理中的递归形式进行计算。
2.4.2 卡尔曼滤波算法
上一小节中给出的一阶马尔科夫动态系统递归贝叶斯估计算法,对于一般的系统很难得到解析解,如果要得到解析解,必须做出一些假设。本节中将介绍的卡尔曼滤波(Kalman Filter,KF)就是在对动态模型和观测模型做出特定假设基础上的一种解析的递归式贝叶斯估计算法。KF由卡尔曼(Kalman)于1960年首次提出,它是一种最优的递归式数据处理算法。由第2.2节可知,存在很多种最优判据,而在给定假设的前提下,无论是哪种最优判据KF均为最优。KF是将所有可用的测量数据、系统和测量设备的先验信息组合在一块对待估状态进行估计,可使得误差在某种统计意义下最小。尽管如此,KF不需要对以前时刻的所有测量数据进行存储,并在有新测量数据的时候进行重新处理,这一特性对于滤波算法的实际实施至关重要。
相对于仅考虑测量方程的最小二乘算法而言,KF的优势是能够实现对动态变化状态变量的估计;相对在频域空间设计的维纳滤波算法而言,KF使用状态空间法在时域空间设计,所采用的递推形式便于在计算机上实现。由于上述特点,KF理论一提出便立即受到工程界的重视,阿波罗飞船导航系统的设计是KF算法早期应用的成功案例。目前,KF理论作为一种重要的最优估计理论,被广泛地应用于各个领域,尤其是基于多源信息融合的组合导航系统的设计。
考虑如下的线性定常系统
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0083_0001.jpg?sign=1739341055-JXVtmCpSJRTgpnKlfRD0F0bv6mmfk0hd-0-42e92897ac206dd533d5c5f37b7c2e6b)
式中,x和z分别为状态和观测量;Fk,k-1为tk-1时刻至tk时刻的状态转移矩阵;Hk为测量矩阵;wk为系统噪声;νk为测量噪声。噪声wk和νk均为高斯噪声,均值为0,并且满足以下条件
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0083_0002.jpg?sign=1739341055-kkFzinyq3gh0kueWuunvVJduXATOdj1U-0-8479efcde4df81660d5871c18f726816)
式中,Qk为系统噪声方差阵;Rk为测量噪声方差阵;通常假设Qk和Rk为正定阵。
此外,状态的先验概率分布满足,其中
代表高斯概率密度函数(见附录D5.1)。
式(2-69)的概率密度模型可以写作(证明参见附录D5.2节)
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0083_0005.jpg?sign=1739341055-6MTQAkacMedulfgA0sNEGqBm1Qji2kSc-0-6b5eef86accb57ec4dc5ddb253fe24dc)
定理2.5:考虑如式(2-69)所示的系统模型,则定理2.4中的贝叶斯滤波方程具有解析的表达式,且
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0083_0006.jpg?sign=1739341055-NkZbJyMHw0P2NPraRQPFQIcJj1eixeMe-0-391929469ed2d9222bfa890859b6fc1c)
上述分布的参数由下面的滤波预测方程和更新方程给出。
滤波预测方程为
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0083_0007.jpg?sign=1739341055-8XBZ1lMDmfhZwGYVELiC77lCRqSRg4UN-0-1c330e7b88a1cfa6249899775a167112)
KF更新方程为
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0083_0008.jpg?sign=1739341055-rDIJdrFAmXpGbWMamDDwVl8oFqG6U0ik-0-65241046de439e4c98dd72b2559dbcb6)
其中
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0083_0009.jpg?sign=1739341055-P9GbzyJDed1VKpdDnXuYNlvKCG2cR6XN-0-c53f599a200f4cf2e43a65ef777db8df)
在证明这一定理前,首先给出联合高斯概率密度函数的两个重要属性(具体推导参见附录D5.2):
属性1(由联合概率密度函数到条件概率密度函数):如果x和z满足联合高斯分布,则仍为高斯概率密度函数,并由下式给出
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0083_0011.jpg?sign=1739341055-m8ruZgd2Wo0Y6XkaCbjFvQlQIClX87Dz-0-a64cbc05e971848a437b3c9db970c608)
其中
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0084_0001.jpg?sign=1739341055-WAYQp5YZGtfIQlI2ChIuRaDOCm6pkNOD-0-449699decb49fa3ffd9f7e2b7fde0208)
属性2(由条件概率密度函数到联合概率密度函数):对于线性方程
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0084_0002.jpg?sign=1739341055-aWMtVlnF35rIjeixBhnLtSBzpbhZwy8N-0-6c51c420a043c314e402f79458d53ca8)
如果
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0084_0003.jpg?sign=1739341055-uHBoO3dALTXH77CSaMUKDMJIO8PmK365-0-af350667816512e5408b138fe9b1fe48)
且。记
,则
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0084_0006.jpg?sign=1739341055-aUYhqzxmz6dmHi8mrBHZ8gFicOsezerk-0-d67dd23879e58e144b40ac31f35d1f9e)
接下来证明定理2.5,证明的主要流程如图2-2所示。
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0084_0007.jpg?sign=1739341055-NsordUcnDxtNysut73tIMv0BpA3ou5h0-0-2e9fa38f823bfbe87fb0e5d7e8b5bfe4)
图2-2 KF推导流程
具体步骤如下:
假设
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0084_0008.jpg?sign=1739341055-7SBwpwKVu4M0dxTqJ6tasjxuubF2dPIt-0-740e533a8a09565b0b24f51c8b1739f8)
由马尔科夫属性1和式(2-70)第一式可以得到
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0084_0009.jpg?sign=1739341055-3IdxlJfJsn6LvdTE0IsVeFg9oeYqVGMu-0-b2b3a692eb15378da724e7707e9af738)
式(2-80)和式(2-78)对应,式(2-81)和式(2-79)的第一式对应,则由式(2-79)的第二式可以得到
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0084_0010.jpg?sign=1739341055-G7PyYGL4byQCryIbBvLVjnOEgP7W1GAt-0-27b8bf70666ecb4b9542c94f38b8cbdb)
从而可以得到
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0085_0001.jpg?sign=1739341055-s9xjTP4VFTKQgU4S2ppzYZNYxhxTi1dl-0-a76be0a841d6f3e190b084d235edc552)
记
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0085_0002.jpg?sign=1739341055-iBQffv4eAeX6GU07GaC8L259A3N4Y95T-0-c50ae33c0595b0c0ac6c15cafa8b97a5)
则式(2-83)可以写成式(2-71)第一式的形式,其中
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0085_0003.jpg?sign=1739341055-RftEvvjxsrzryh4lDzpVDzjZlKp1RuTR-0-ae39cd410eb0f0dbf9ba966fd7dc1414)
由式(2-79)的第二式和式(2-70)可以得到
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0085_0004.jpg?sign=1739341055-pljKweGBOAPTUqHkLtvEfWeJmv1AEwZY-0-a5055fb7a5f20973c75bef055cc36f6e)
再由式(2-75)可以得到
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0085_0005.jpg?sign=1739341055-Pqbebfvk4QIOrrr8qZwyYa5Oesu6U2p8-0-0642690232c08ae5f201a70db71d60b0)
其中
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0085_0006.jpg?sign=1739341055-TR76WxBOoDL4Q9mpUV1rrJJZRLiKL7WU-0-435c00da709a6077bd9949bebfd9fbfe)
记
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0085_0007.jpg?sign=1739341055-MNusaWAwmM08FKcYbT5oL58i1ZmTNK8A-0-b1ab38acb1751c38d728ee986a147cf6)
则可以得到式(2-73)。从而得证。
应当指出的是,证明的过程中利用了假设式(2-80),这一假设可以利用归纳法证明成立。
注:
(1)根据MMSE估计的最优判据可知,Pk为最优估计。
(2)KF更新方程有着几种不同的形式,参见表2-2(与LMMSE估计的表2-1对应)。
表2-2 几种不同形式的KF更新方程
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0086_0001.jpg?sign=1739341055-JKO3poEhiRTwanN7TeYLTxnl4YblttoI-0-acb160645ad266590f47c48cba1366a5)
(3)如果初始状态为任意分布,均值为、方差为P0;过程噪声wk是均值为0、方差为Qk的独立过程;测量噪声νk也是均值为0、方差为Rk,且初始状态、过程噪声、测量噪声相互独立,则定理2.5给出的KF是最优线性无偏估计。
2.4.3 扩展卡尔曼滤波算法
Kalman最初提出的滤波基本理论只适用于线性系统,并且要求测量方程也必须是线性的。但是,在工程实践中遇到的物理系统,其数学模型往往是非线性的,如在航天器自主导航技术研究中经常用到的轨道动力学方程是非线性的,天文自主导航系统中对应的视线方向测量方程也是非线性的。为了将KF算法推广用于非线性系统,在KF基本理论提出后的10多年时间里,Kalman、Bucy和Sunahara等学者致力于研究KF在非线性系统中的应用,提出了扩展卡尔曼滤波(Extended Kalman Filter,EKF)算法。
EKF算法设计的基本思路是通过截取系统方程和测量方程中非线性函数泰勒级数展开式的一阶项,对系统进行线性化,然后将KF方程用于线性化模型以获得状态估计值。EKF算法是针对非线性系统进行状态估计最常用的算法,如美国的“星尘(STARDUST)”号探测器采用EKF处理光学成像测量信息,从而确定探测器的位置。基于EKF的姿态确定方法广泛应用于高精度卫星姿态控制系统,如美国的哈勃太空望远镜(H ST)和日本的先进陆地观测卫星(ALOS)等。
用于EKF算法研究的非线性系统模型为
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0087_0001.jpg?sign=1739341055-rZOMbTPEOFSH5PrBzt3a2fYuhS4Io8gi-0-e50545aa5221714a9aa8df922114b398)
式中,xk∈ℝn为状态变量;zk∈ℝm为观测量;f: ℝn→ℝn为状态转移函数;h: ℝn→ℝm为测量函数;wk和νk均为零均值白噪声,并且满足以下条件
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0087_0002.jpg?sign=1739341055-NZPFsCyvWZlxSFqu0cIegMlRbXN9758r-0-5d76c752e6ec617bff8f9fc13ad435bc)
对于式(2-89)所示的非线性系统,很难找到一种严格的递推滤波方法,通常采用近似方法来处理非线性滤波问题,应用比较广泛的近似方法是非线性系统的线性化。为了针对非线性系统应用KF方程,作如下基本假设:状态变量预测值与实际值之间的差能够用一个线性方程表示,该线性方程能够足够准确地对滤波器的实际误差特性给予描述。这个基本假设在工程实践中往往可以得到满足,描述预测值与实际值之差的线性方程称为线性干扰方程。常用的EKF算法是针对线性化后的模型设计的。
EKF算法方程描述如下,定义预测误差为
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0087_0003.jpg?sign=1739341055-drpCSFuUyXgnyOj1quT2QYsG0nRXKsAN-0-aeac7c28d76831c446285dfce1f3d674)
式中,表示状态变量的预测值。线性干扰方程的形式为
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0087_0005.jpg?sign=1739341055-J24oablOHnU83yOb24IGzCBOYIsZB9Xp-0-6a084a7f029d5bfeb27c2dd8234c83c7)
式中,称为雅克比矩阵;Δzk=
。
在线性干扰方程的基础上,仿照线性KF基本方程,不难给出对偏差Δxk进行估计的KF方程为
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0087_0008.jpg?sign=1739341055-tH1YYoFGGYba0xPR2O6Jd8vDzH7wfh5g-0-0ceb9345f158affdf6d759fa65274ea6)
值得注意的是,在应用过程中,往往不对偏差Δxk进行预测,而是直接设置状态偏差的一步预测值,此时
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0087_0010.jpg?sign=1739341055-FXti6vz9pHn5ptjR9CcALZRS2A755fMo-0-98e9479797ff74413a4237097573a40a)
一般通过系统方程式(2-89)的第一式对状态变量预测,即
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0088_0001.jpg?sign=1739341055-oZPGHtLB973WnsKaakiIXbNkbEnL2u0m-0-4e35a0a1392196f30bdec4fe05502cd4)
利用估计得到的状态偏差Δ^xk对预测值^xnk或^xk|k-1进行修正,得到
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0088_0002.jpg?sign=1739341055-z0eB1QoPmXuv77ltyFvmKTZXlCx0wFWL-0-05f95b6c7247966ff7016b88e080821e)
为了精确实现对状态变量的预测,针对式(2-95),人们提出了多种改进方法,如数值积分、求解微分方程等。相应的,为了提升状态变量修正的效能,针对式(2-96),也有不同的实现方法,如求解非线性优化问题等。在实际应用过程中,应综合考虑系统非线性程度、系统噪声和测量噪声影响的权重、精度要求,以及计算量等因素,设计或选择适当的状态预测和更新方法。在系统动态变化不明显或呈现周期性变化规律时,可采用常增益或周期增益滤波器,避免了对滤波增益阵Kk的递推解算,能够降低计算量、减小滤波算法的复杂程度。类似算法广泛用于基于陀螺和星敏感器的卫星姿态确定系统。
2.4.4 无迹卡尔曼滤波算法
1.问题描述
20世纪60年代提出的EKF将KF技术推广用于非线性系统,其基本思路是通过截取状态方程和测量方程中非线性函数泰勒级数展开式的一阶项,对系统进行线性化,然后将KF方程用于所获得的线性化模型以获得状态估计值。EKF算法历史悠久,原理直观,已被用于航天器定轨和飞行器导航等多个领域。EKF算法的主要问题在于,在线性化过程中忽略了非线性函数泰勒级数展开式的二阶项和其他高阶项;在初始误差较大或系统受到外界干扰的情况下,有可能引入较大的线性化误差,从而降低滤波精度。为了解决这一问题,人们提出了二阶滤波等改进方法,试图通过引入非线性函数泰勒级数展开式的高阶项来改善滤波性能。虽然采用二阶滤波可以提高滤波精度,但对于复杂的非线性系统而言,计算非线性函数的高阶导数往往比较困难。另外,基于贝叶斯估计理论和蒙特卡洛方法的PF算法也可以解决非线性系统的状态估计问题。但是,应用PF算法需要对大量粒子进行预测和更新,对于性能有限的星载计算机而言,计算负担较重。
20世纪90年代,关于无导数(Derivative Free)非线性滤波算法的研究取得了较大发展,其中比较著名的包括UKF、二阶插值滤波和容积卡尔曼滤波(Cubature Kalman Filter,CKF)算法等。UKF和二阶插值滤波算法分别采用无迹变换和二阶插值技术描述非线性函数泰勒级数展开式的二阶项和其他高阶项,基本思路是通过确定性方法选择少量样本描述状态变量的均值和方差特性,然后对各个样本分别进行非线性变换,并将它们通过非线性系统方程和测量方程后的样本均值和样本方差分别作为状态预测值、测量预测值和相应方差阵。UKF可以达到与二阶滤波相似的精度,而不必对非线性函数进行求导,应用更简便且其计算量与EKF处于同一个数量级,不会显著增大计算负担。CKF是近年来提出的一种新型非线性滤波算法,该算法基于球面径向规则设计,经过严格的数学证明,其逼近非线性变换后概率分布的精度优于UKF,并且能够解决UKF在处理高维非线性状态估计问题时滤波性能不佳甚至发散的问题。应当指出,UKF、二阶插值滤波和CKF在原理上都是基于一组加权样本点来逼近非线性状态的统计特性,且采样过程都是根据确定的数学表达式来实现的,因此,它们可以统一归类为确定采样型滤波器。本节主要以UKF算法为例展开论述。
EKF、UKF和PF算法的区别在于描述状态变量通过非线性系统后统计特性的方式不同,其中,EKF通过线性化技术实现对状态变量均值和方差的递推计算,UKF通过少量根据确定性方法选择的样本点描述随机变量通过非线性系统后的统计特性,而PF通过大量随机样本描述状态变量通过非线性系统后的分布,如图2-3所示。
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0089_0001.jpg?sign=1739341055-CkcqBin3EDMFXMxuod9rml4qNKlQ5Lz2-0-b2b465c22c3d2fc0d45c8f745962fbf0)
图2-3 常用非线性滤波算法实施状态预测的方式
常用滤波算法的性能特点如表2-3所示。
表2-3 常用滤波算法的性能特点
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0090_0001.jpg?sign=1739341055-as4JTtQZqi7szQq4Se1oL2qlmknat7wL-0-d8c5627b1c2fb0195941175bfabcef55)
为了满足航天器自主导航研究的要求,需要从精度和计算量两个方面着眼,选择适当的非线性滤波算法。从精度方面考虑,采用UKF可以获得高于EKF的估计精度;从计算量方面考虑,EKF的计算量与状态向量维数的3次方成正比;UKF算法的计算量约为EKF算法的3倍;而PF的计算量与所选择的粒子数的大小有关,通常可达EKF的数百倍甚至上千倍,对于性能有限的星载计算机而言,应用PF算法计算负担较重。因此,建议集中研究EKF或UKF算法,即采用线性化或无迹变换技术,计算状态变量通过非线性系统后的均值和方差,在滤波过程中实现对状态变量的预测。此外,应用UKF算法不必求解雅克比矩阵,适用于地球重力场导航等雅克比矩阵不易求解的情况。UKF算法在2007年发射的美国海军NPSAT 1(Naval Postgraduate School Satellite 1)卫星上得到应用,该卫星的姿态确定系统利用UKF算法处理三轴磁强计的测量信息,获得卫星姿态及其变化率的估计值。
2.UKF算法流程
无迹卡尔曼滤波最先由Julier提出,又经Merve扩展出了多种算法。与EKF不同,UKF中采用无迹变换(Unscented Transform,UT)技术代替线性化,以描述状态变量通过非线性系统方程或测量方程后的均值和方差。所谓无迹变换,指的是通过确定性采样方法在状态空间选择若干样本(即Sigma点),用这些样本的分布来描述状态变量的均值和方差特性。可以证明,对每个Sigma点分别进行非线性变换,得到一组新的Sigma点,则这些新样本的均值和方差能够以较高精度逼近真实的状态变量经过非线性变换后的均值和方差。图2-4所示为二维状态变量的无迹变换示意图。
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0091_0001.jpg?sign=1739341055-amdNfLX30SGzvnZv5sJFYOltvfGjnEJX-0-62fdfb02530329bef299bd5142808f7b)
图2-4 二维状态变量的无迹变换示意图
无迹变换的计算方法为:假设x的均值为,协方差为Pk,通过非线性函数x′=f(x)传播到x′,要计算x′的均值和协方差,需要构造一个由x的2n+1个采样点χi(称为Sigma点)组成的矩阵
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0091_0003.jpg?sign=1739341055-QWLmkr3E82aT99glcl0cuUqVshLPyvhQ-0-58a711b2c855220d575178601f15c33d)
式中,表示矩阵平方根的第i列;λ=α2(n+κ)-n,是一个标量。
常数α一般取小的正值(如10-4≤α≤1),控制西格玛点的分布状态。调节α可以减小非线性方程的高阶项影响。κ的具体取值虽然没有界限,但应确保矩阵(n+λ)Pk为半正定矩阵。对高斯分布的情况,当状态变量为单变量时,取κ=2;当状态变量为多变量时,取κ=3-n。
本节先给出UKF的算法流程,再说明该算法相对EKF的优势。
第1步:初始化,假定初始时刻未知状态的先验分布均值为,方差为P0。
第2步:选择Sigma点,给定上一步的状态估计值及其方差阵Pk-1,Sigma点可用下式选择
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0091_0007.jpg?sign=1739341055-zzKbXKVCNrPsW4fsAvXU6KEv56Znu1FV-0-a0029a66e456217ed3420d23af845eaf)
式中,可调参数a描述了Sigma点的散布,一般取为一个小正数;在一些特殊情况下,通过适当选择a的取值可以消除某些高次误差的影响。表示矩阵平方根的第i列,矩阵
是通过矩阵分解得到的,满足Pk-1=
。对于l维系统来说,需要选取2l+1个Sigma点。
第3步:预测,将各Sigma点分别代入状态转移函数进行计算,得到一组新的样本
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0092_0004.jpg?sign=1739341055-eyDRFMG76eYk1X5tCA1MjOE0P9i7KMke-0-2e2e38178a79ba04f5442ff33a6bcc90)
状态变量的预测值和方差阵可按下式计算
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0092_0005.jpg?sign=1739341055-WNeQlwnYwgH1RnYpaMWySISPeIQYtKKA-0-961652469b3e971d39d80a4292b048f0)
式中,权值为
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0092_0006.jpg?sign=1739341055-lSYOtTZE6C9ITneti0gganhKZ51i35bn-0-ec7950610d8cc595e14928ef4ec1cca6)
第4步:更新,根据预测的结果和Pk|k-1重新选取Sigma点,即
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0092_0008.jpg?sign=1739341055-m9V8F2KdLACjaxUonHGKwhjUazGzp4Ww-0-f2290c41698f386b85ee1accdba46a04)
将重新选取的Sigma点依次代入测量方程进行计算,得到
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0092_0009.jpg?sign=1739341055-lws0iWOP7iAlmoOYSXNVV6Mjv9nX23cV-0-0eadd7053cd88261eb221999996b0cbc)
观测量的预测均值及其方差Pzz可按下式计算
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0092_0011.jpg?sign=1739341055-mB89H9Yt6CF3nzxfaRHj98MfMwBJN5gv-0-3b47bef41fc1abfdfd65449d4700cf15)
状态和观测量的互协方差阵Pxz为
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0092_0012.jpg?sign=1739341055-ValeovkJi1BYzt4p4fDmZh1clLfilFMH-0-ea205ab90605bb81e83e8a4424eb8be7)
通过观测量对预测值进行修正,得到状态估计值及其方差阵为
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0092_0013.jpg?sign=1739341055-CODLWq6cIH9eFPBXiQXdAnojT5tp0JEl-0-18adfa2c401d35da12acff2d2e570498)
第2~4步反复迭代运行,就可以得到状态xk的估计值。
UKF与EKF的不同之处在于采用无迹变换技术代替线性化,计算状态变量通过非线性系统方程与测量方程后的均值和方差。可以证明,采用无迹变换技术可以精确描述非线性函数泰勒级数展开式的二次项,而EKF中所用的线性化技术仅考虑了泰勒级数展开式的一次项。相对EKF而言,UKF算法估计精度更高,并且不需要求解雅克比矩阵。
UKF算法的优势在于对预测分布均值的描述比EKF更准确。为了便于比较,先推导预测均值的表达式。将状态变量先验分布的均值和方差记为
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0093_0002.jpg?sign=1739341055-kT6DzEX8qkJDURTmkNetHhTwqsbeBrOk-0-7429fa0491272113f80e2f80e42e4daf)
则预测分布的均值为
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0093_0003.jpg?sign=1739341055-OX36FXEXzQgLgQ6qZwl07ibhQYvCYAEO-0-22137040e8f0d5a6fd8223199c4d3ab4)
式中,为估计误差。不难得到其统计特性为
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0093_0005.jpg?sign=1739341055-Pgb4r0h912PynbjllnO9PAZXxtlE0uXc-0-2d7b2facebdba851f1d9569b548e60c6)
将在
处展成泰勒级数,得到
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0093_0008.jpg?sign=1739341055-D8hCMaukemfzuKIG5YdmBi2AXfgSwUQn-0-a780a5c96db9ad6c979dfc0c8fc002fb)
式中,ϕi∈ℝl的第i个元素为1,其他元素均为0; fi为f的第i项,从而
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0093_0009.jpg?sign=1739341055-QEf0YpK0bzwBsMZKggbfASRZzCLKBk0V-0-192d219f4f1def1e33921ed9edd1db2a)
式(2-107)中给出了非线性函数泰勒级数展开式的二次项。
接下来计算UKF中通过无迹变换得到的样本均值,将式(2-98)代入式(2-99)第一式可得
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0093_0010.jpg?sign=1739341055-8ZuQjR4J460x1I7Fom8eNCtpIptrwrIB-0-e7e21487c23ce4295ad29dd3cdaf7551)
将和
在
处展成泰勒级数,得到
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0094_0001.jpg?sign=1739341055-DsI48wCmNPBS9QubPC5j0cZ2RszCROpD-0-bc3f189d1096fb056007bb304b677c7b)
整理上式得到
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0094_0002.jpg?sign=1739341055-9tY0sKRxmvdrOSbSyuzgDiPws9DuZZy5-0-d598c87b5bd1414c2dff78039eac7660)
对比式(2-107)和式(2-108)可知,UKF中的样本均值与E{xk
的泰勒级数展开式的一次项和二次项是相等的。由此可知,采用无迹变换技术得到的样本均值
以二次精度逼近预测分布均值
。同理,样本均值
将以二次精度逼近预测均值
。
回顾在EKF中对的近似方法,为了与UKF区别,用符号
表示EKF的预测值
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0094_0011.jpg?sign=1739341055-Pbf7VmtKzJGNGtGq7k2woEqJldNJRxAB-0-150a7a5ef791df3a32df9a4c4ba1265a)
观察式(2-107)和式(2-109),不难发现,EKF简单地忽略了非线性函数泰勒级数展开式的高次项,给出的预测值二次项为0,与预测均值E{xkz1∶k-1}的二次项不符合。可见,EKF的对状态变量均值的预测误差体现在泰勒级数展开式的二次项和其他高次项中,而采用UKF得到的预测误差体现在泰勒级数展开式的二次以上项中。因此,UKF对预测均值的描述更精确,理论上能够达到比EKF更高的估计精度。
2.4.5 约束卡尔曼滤波
1.问题描述
在状态估计问题中,系统的状态往往满足一些代数约束,如果设计估计算法忽略这种约束,则估计的性能未必最优,而且最终估计未必满足约束。在本节中将对考虑等式约束的卡尔曼滤波进行研究。考虑两种典型的等式约束:线性等式约束(Linear Equality Constraint,LEC)和二次型等式约束(Quadratic Equality Constraint,QEC)。其中线性等式约束可以写作
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0095_0001.jpg?sign=1739341055-bZwc0bKkPMdok3qdbqVYWMPrHLl3eAfT-0-8567ca28aa6f87962ce393e627421f38)
式中,D∈ℝs×n; d∈ℝs; s≤n为约束的维数。
通常情况下,可以假设D行满秩。如果不是行满秩,则意味着存在冗余的约束,这时候可以去掉冗余的约束使得D行满秩。
对于二次型等式约束可以写作
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0095_0002.jpg?sign=1739341055-5CmAeuOcMmloVFenjec2W6UFlq8iPYB9-0-0470859c497f2cfec41389e7521c6964)
式中,Ai=ATi∈ℝn×n; q为等式二次型约束的个数。
值得注意的是,对于一般的非线性等式约束模型有
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0095_0003.jpg?sign=1739341055-iE6EZqhRTfhm93l1bmnc7M2ERDUKLZFU-0-2397b2dcd5676984630dae2c60a7cee6)
其在参考点泰勒展开可以得到
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0095_0004.jpg?sign=1739341055-jmcEGf1tjeuGinp2UOsPj2FiFfiyBJsf-0-87adec463f490efab3741c5a6f3332bf)
式中,De为微分算子,有
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0095_0005.jpg?sign=1739341055-hRvqHVXhkd2j8xE55taPMDeQLKa3Z29N-0-1a94edc73576e0bc763f591494643405)
式中,ej为e的第j个元素,e为扰动量。如果仅考虑一阶项,则约束可以写作线性等式约束的形式;如果考虑二阶项,则约束可以近似成二次型等式约束的形式。
对于最小均方误差估计问题,给出状态的先验估计和方差阵
以及测量方程
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0095_0008.jpg?sign=1739341055-Nx9NzsTyLGlznk7kDWTCnutJeCSXzHQq-0-a3940ee08ecd42903a4abf55990c149b)
式中,z∈ℝm; H∈ℝm×n行满秩为测量方程。求后验估计,使目标函数
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0095_0010.jpg?sign=1739341055-aKn3x5wgogTP1jC5NxqHNG1kn3BjXX9V-0-8450002604275b8bc57d15aeae60e113)
极小,且满足约束方程式(2-112)。
2.状态线性等式约束滤波
状态线性等式约束滤波问题可以描述为:寻找最优估计∈ℝn使得目标函数式(2-116)最小,且满足式(2-110)给出的约束方程。
将写作先验估计
和观测量z的线性组合形式,即
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0096_0004.jpg?sign=1739341055-86ALhXhxJwuM0y8ob61poyExc6Hj6GgN-0-a49a891f941fa98dcbc13591870e8cef)
式中,K∈ℝn×m; N∈ℝn×n; n∈ℝn为待确定矩阵。在和η=0的情况下有
=x,因此
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0096_0007.jpg?sign=1739341055-497ogVF8j0OPmlvGXpzqjlEdUUieaBsy-0-6597906da538c8f24d292a419c14b758)
从而使K、N和n满足约束
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0096_0008.jpg?sign=1739341055-V3k6x7GAQrutGqldgEwseJow6PBM5xJu-0-fb6a97d2f0961947adebc71edfe75123)
将式(2-119)代入式(2-117)中可以得到
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0096_0009.jpg?sign=1739341055-iOPvyGCbWlY4l3260uctgFCDDfrQl70f-0-fb35268edc932f2048a9bf7b7e6d1483)
式中,。
对于非约束线性估计问题,K的最优值可以通过对式(2-116)求极值得到。这里考虑的是线性等式约束估计,则极值函数可以重新写作
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0096_0011.jpg?sign=1739341055-FR7gMCVV6mje4Wu7PCj2ptXaWaibT1k5-0-f0c3bd65f45a0c6cba33eb878501d79f)
假设测量噪声和先验估计不相关,将式(2-120)代入式(2-121)中可以得到
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0096_0012.jpg?sign=1739341055-3QCgxLVyj8NNUhhiw3sMjlVXworp8epe-0-f379e1bdc0febf2474d86d72265787a1)
式中,R为η的方差阵,P-为的方差阵。
利用最优必要条件可以得到
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0096_0014.jpg?sign=1739341055-TSbL5s3lNd0RJnbXwyqe3q2kzbfb2ywi-0-47d25f404ccaa12efbbf73c3046aabf8)
由式(2-123)第一式可以解得
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0096_0015.jpg?sign=1739341055-YvP6I4TDkLohBjjLtfEOwlnBwGb18jGx-0-ae836f318d4b97de50a4d7f0fc2256ac)
其中
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0096_0016.jpg?sign=1739341055-m8ebFNCHke6pxMfpAaQ24UpzgWcm8tEi-0-8c6b6a7ab356a06541d5a11b8c0ec6a1)
将式(2-124)代入式(2-123)第二式可得
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0097_0001.jpg?sign=1739341055-lFkToGLIkBN8I52TRSe0lbHuFJOj4QGY-0-dbfc4e1e4923e3fa1b183fa8f4e42c44)
经过简化可以解得
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0097_0002.jpg?sign=1739341055-bTeTMCJZoWmbh52DTrwgsbosyR7m191b-0-f25780ece6c52d01c5bdac8d50928cc7)
其中
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0097_0003.jpg?sign=1739341055-oscQK7wjh6nluf0a6BvpU8DadpeVZwoi-0-41d457e08c09d91649171683ae8bd3e4)
将式(2-127)代入式(2-124)中可以得到约束估计的卡尔曼滤波增益矩阵,然后利用下式对方差进行测量更新
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0097_0004.jpg?sign=1739341055-h4zSm3pwi6FD6WqnG7eY1zwC1TS78zol-0-2ed78bb774252414c7f547d658601f1d)
3.状态二次型等式约束滤波
状态二次型等式约束滤波问题可以描述为:寻找最优估计使得目标函数式(2-116)最小,且满足约束方程式(2-111),这里仅考虑单个约束的情况。和状态线性等式约束滤波类似,
可以写成式(2-120)的形式。
相应的,二次型状态等式约束滤波的极值函数可以写作
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0097_0007.jpg?sign=1739341055-wavw5bB2n7h2N2jGKH9rpCms5D4lrs5A-0-b34443db9d0ccec0e0e844c74b7c9adf)
式中,λ为拉格朗日乘子。假设测量噪声和先验估计不相关,将式(2-120)代入式(2-130)中可以得到
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0097_0008.jpg?sign=1739341055-q0HbjpMhd2kNVcXZ8jvdbCUSxbKIySjO-0-2844689d1b98c91fbe2b68b110663de4)
根据附录E的最优必要条件可以得到
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0097_0009.jpg?sign=1739341055-6gslIJICDp6sBhZKdOih5i2y8kUvMcc5-0-08723651725ad5d93a004ebac160cbe1)
上式第一式可以重新写作
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0097_0010.jpg?sign=1739341055-yMAikNwcPRHdp6doIbB8K5F29zC29sWo-0-ea69615f3838c05e09700ecaf8de1dd8)
式中,W由式(2-125)得到。
式(2-133)有K和λ两个未知参数需要求解。典型的求解方法就是从式(2-133)的第一式求解出K为λ的函数,然后代入式(2-133)的第二式中求解出λ。式(2-133)第一式为离散Sylvester方程。如果A=I,可以直接对矩阵求逆得到K的解。如果A为一般矩阵,K的解则没那么直观。下面进行具体分析。
对于m维空间,存在m-1个向量βi(i=1,2,…,m-1)使得βi⊥βj(i≠j)和ϵ⊥βi。因此,λϵϵTW-1的m个特征值为0,0,…,0,λε,其中ε由式(2-137)给出,相应的特征向量为Wβ1,Wβ2,…,Wβm-1,ϵ。假设存在m个αi∈ℝ,i=1,2,…,m,使得
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0098_0001.jpg?sign=1739341055-oLAZtJWCs1aQ6LTWSDabXsQJTvDubovt-0-0f013e6bf3266db1abd2195fdb4a5ae0)
因为W>O,式(2-134)左乘以ϵTW-1得到α1=0,然后左乘以βTjW-1(j=1,2,…,m-1)得到
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0098_0002.jpg?sign=1739341055-uASdGjVge4l8gD8d9qo6OWSQl3jT1X0N-0-7cb65987eba9206169344c4e2038ee34)
这意味着λϵϵTW-1的m个特征向量线性独立。换句话说,λϵϵTW-1可以对角化为
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0098_0003.jpg?sign=1739341055-lshBiECIAimS4pu0fMuxttfpxZu8awZ5-0-06fa71e43b4b7db4d721f3561aabc328)
式中,([γ 1,γ 2,…,γ m])=diag([0,0,…,0,λε]),V=[W β 1,Wβ2,…,Wβm-1,ϵ]。由于对称矩阵同样可以被正交矩阵对角化,因此A可以写作
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0098_0005.jpg?sign=1739341055-J0GvYFAHvsEsBKpPK5LSewddGYmEsJZG-0-7585eeda9cef94533ae3e8a0fa625bb4)
式中,U为正交矩阵;Ξ=diag([ξ1,ξ2,…,ξn])。将式(2-137)代入式(2-133)的第一式,并左乘以UT、右乘以V可以得到
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0098_0006.jpg?sign=1739341055-JU2vUsZMqT9pOejOTs41n9qf3EdyEr56-0-9ba6c85b83a414d992fe0f34d9b4f3a0)
式中,K~=UTKV;C~=UTCV;C=(P-HT-λA^x-ϵT)W-1。记:K~和C~的第i、j个元素分别为k~ij和c~ij,式(2-138)可以写作
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0098_0007.jpg?sign=1739341055-RHDEqPnVxkhNBUZpf0dVszksc2tRUarW-0-090493a3234099d790e7ab637a6c769c)
这意味着
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0098_0008.jpg?sign=1739341055-0vI3a04ou50jrPg54WEHOlR78wqrMK0K-0-be2bab087d3bc533e9e17d997bb7f1f4)
一旦计算出来,则有
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0098_0010.jpg?sign=1739341055-mieJgncRI5SnucfSzP5ep66tEJqnD0fn-0-a3758cdcdfaec7a465b5f62a1bd694ac)
对于状态二次型等式约束滤波,很难得到后验估计的误差协方差阵。但是基于扩展卡尔曼滤波的假设,可以将其近似为式(2-129)。
注意到c~ij和γj依赖于λ,因此需要先计算λ才能计算K。对式(2-133)的第一式右乘以ϵ可以得到
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0098_0011.jpg?sign=1739341055-0KUAkn4LJGzSiw10BEItFGrurIHkrjsF-0-18165008dc8187cb83a7a81e4e4e6c50)
式中,d=P-HT W-1ϵ。
则
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0098_0012.jpg?sign=1739341055-T8Cs6Q6pQ4SBNpQrgP460JVPr0mYyO5R-0-893ef6ac65a947d19dac16027b0192ac)
将式(2-143)代入式(2-133)的第二式可以得到
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0099_0001.jpg?sign=1739341055-q1Eh0n1PsqIaiGh4iuWxWf47mURwd9LL-0-b96ba74263fda4396caffae3c6e964b1)
对式(2-144)进一步化简可以得到
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0099_0002.jpg?sign=1739341055-RQmZOu8QWUUwq99r4r5Ozv5OhVTA6ehC-0-adca7e062cd4be2e97162948c49ba266)
式中,h=d~+; d~=UTd;
。
式(2-145)可以写成标量的形式
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0099_0005.jpg?sign=1739341055-nvyB4TSmhNmegSNzcSeb7CFokrGBMSxL-0-b80f8254ddf0237710ade6c739c490fb)
式中,。通常情况下,可以采用牛顿迭代法对
进行求解。但是这种方法只能求解出一个解,而实际上式(2-146)有2p个根,其中
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0099_0008.jpg?sign=1739341055-GUsdKeamORSaOKCv98QjTrmDXsL4V50W-0-559bba79c1b9083aa64c8a3e9d3095d5)
式中,q为A的不同特征值个数。因此牛顿法得到的解未必是最优解,而且牛顿法未必能够保证收敛,这就需要更为有效的方法去计算。令
为A的q个不同特征值,并定义多项式函数为
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0099_0011.jpg?sign=1739341055-fqwLGpyDDeimXOid51QJHq9Xe2RfGKeQ-0-1074526fb568570869fc29770aca4f4c)
式(2-148)可以重新写作
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0099_0012.jpg?sign=1739341055-5IlQyrLCxs2H5NXsexSF0i06H0Y7MMJK-0-334b2010c87414fe5d030126c155e187)
式中,可以由差分方法计算得到。注意到式(2-148)和式(2-146)有相同的解。构造式(2-149)的伴随矩阵为
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0099_0014.jpg?sign=1739341055-ULo1sIO68W3ov76tRtgmZXu3tC2x4leR-0-29e079bbcc4596a936bfa328a605c3d5)
G的特征值对应式(2-149)的解,通过求解G的特征值可以求解出式(2-149)的所有解。因为G有2p个特征值,必须确定最优的那个特征值。利用约束优化的二阶充分条件可以得到如下定理:
定理2.6:如果为G的特征值,且对于∀i∈[1,r],r=2,3,…n,有
,则
对应最优泛函指标。
证明:由附录E中的约束优化内容,可以将最优泛函式(2-131)可以重新写作
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0100_0004.jpg?sign=1739341055-xIqlrqy7gwrTSb5CeN8hdyRVdVW2ZGwo-0-da7562fc72f41ee814debb1f71c8e375)
其中
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0100_0005.jpg?sign=1739341055-YSzPphxIw2B7LIX73Zo3Ngn0ofXwdnNC-0-74e3b24e085fba1c9feca21aa299ba4c)
假设K∗为式(2-142)的解,λ∗为式(2-148)的解。根据附录E约束优化的内容可知,如果g(K∗)=0,且对于所有满足下式
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0100_0006.jpg?sign=1739341055-Z2Lqtkfa8JGpKhuAQIVPWyktrpCCrNQE-0-0686043a6f5afeec54aa240521e91d87)
的z都有
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0100_0007.jpg?sign=1739341055-iSXizcMln7RJRsljKsaG6Ffke84dsp4M-0-5d42df65efa9229b937937d4c9c0e165)
则K∗对应最优泛函。
式(2-153)和式(2-154)中,vec为向量化算子;⊗为Kronecker算子(相关定义见附录C)。
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0100_0008.jpg?sign=1739341055-GpA7zbpHy1MQ3HMiG7bbOjmgzp5NFHVp-0-ed4cc640a0b2b9af57fca02519730f2e)
为了判断式(2-154)的成立条件,构造矩阵
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0100_0009.jpg?sign=1739341055-oEIhZdn6cIazVktyCmsbVj3NHs2Nvs4m-0-5f8b17b00823dd4e6108656f2e423993)
式中,Sr为对称分块矩阵;Lrr(K∗,λ∗)为L(K∗,λ∗)的左上角r ×r矩阵;(Im⊗▽gT)r矩阵为Im⊗▽gT的前r列。对于所有的z·▽g(K∗)=0,当且仅当下面条件成立时,有Q(z)>0,即
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0100_0010.jpg?sign=1739341055-XhJvwIMrNVixTFCWbPhNCAZPOYQpIjAP-0-51b3664589dfb66b5e4edd11ba35035b)
因此问题的关键是判断式(2-157)的成立条件。为了推导方便,仅考虑标量测量数据的情形(m=1)。这样式(2-156)可以重新写作
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0100_0011.jpg?sign=1739341055-jf85tswDvkCY9iENBlwwIDuWN27sndJw-0-a072cd0f49a2780b1fdd915c83e9eccb)
Sr的行列式由下式得到
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0101_0001.jpg?sign=1739341055-BCMvHshmTOXNFHpSn3qnj8idZfII46Zx-0-a838e8066165399012b3c2815c90dcd2)
通常需要对上式进行计算得到det(Sr)的符号,但是det(Sr)<0的一个充分条件为
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0101_0002.jpg?sign=1739341055-JUpg4GIm7b3UUu4uV4H2kIEbFricUgs4-0-00e04ad5fad908e2592795ee258ca4bc)
这是由于det(I r+λε Ξ r)(I r+λε Ξ r)-1的第i个特征值为,从而得证。
注:
(1)状态二次型等式约束滤波的实施流程如图2-5所示。
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0102_0001.jpg?sign=1739341055-cXhaBzAVJX0nGMU0sXLLXyhFaKJ44AI1-0-c0104db4238304bdebc1c4d531662aca)
图2-5 状态二次型等式约束滤波的实施流程
(2)如果A=I,则二次型约束为范数约束。因为=1,q=1,则式(2-148)可以简化为
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0101_0005.jpg?sign=1739341055-O6e7QenUHvRjOVg8j2HAUaPCR2wmAnnI-0-2a6310b74a018e374da0ea7b847ddd60)
从而拉格朗日乘子的解可以写作
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0101_0006.jpg?sign=1739341055-536gwKWMsbXAOPB0ULouSBLaBV8sesf3-0-2661ce9665365240bec3772e713ee56e)
根据定理2.6可以得到上式中的“+”号对应极小泛函指标。
(3)定义Y∈ℝn ×n由下式给出
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0101_0007.jpg?sign=1739341055-rmwtIb4RUjf5pjCm4usN4UX6dUdqgCcM-0-92d0b18a73a3264ef9f8772ee457b513)
则变换x~Y=Yx ~可以将约束方程变换成
![](https://epubservercos.yuewen.com/6C85C2/17640067807576006/epubprivate/OEBPS/Images/figure_0101_0008.jpg?sign=1739341055-QlvkalAqP2TvRh6VCvw4dSmzGRFbXxmQ-0-7b6409cb9e98f1b1defa2ba17666bbda)
式中,的对角矩阵元素为0,1,-1。对于变量
,其二次型系数矩阵最多有两个不同的非零特征值,这意味着拉格朗日乘子最多有4个根。
(4)对于约束系统,约束滤波得到的估计更为合理。比如自旋卫星的姿态估计问题,由于缺少陀螺测量量,而且状态变化很快,常用Markley变量替代姿态四元素表征姿态。Markley变量中的状态量包含角动量在惯性系和本体系下的投影Li和Lb,由于两者是同一个变量在不同坐标系下的投影,满足约束二次型等式约束‖Li‖=‖Lb‖。如果不考虑该约束,得到的估计未必合理。