在实际导航系统中，单个传感器可能会出现故障或零位漂移等特性导致无法正常工作. 例如，GPS因受遮挡等影响无法接收定位信号^[1-2]，惯性导航传感器由于物理上的限制而存在漂移误差^[3]，视觉传感器由于特征点被遮挡无法观测等. 因此，为了保证导航系统的正常运行，需要利用2种或2种以上的观测方式，构成多传感器组合导航系统. Bian等^[4]提出了一种INS/GPS自适应组合导航系统，利用GPS作为惯性导航系统（inertial navigation system，INS）的辅助观测系统，同时设计自适应卡尔曼滤波（adaptive Kalman filte，AKF），用于解决GPS时变的观测噪声. Gao等^[5]使用GPS和激光雷达(LiDAR)对INS进行周期性修正，设计了基于INS/GPS/LiDAR的组合导航系统. Li等^[6]为了保证火星探测器实现精确和安全的火星着陆，设计了基于微型相干高度计、测速仪（miniature coherent altimeter velocimeter，MCAV）和INS的组合导航系统，该系统通过MCAV对惯性偏差和漂移进行校正，提高了组合导航的性能.

组合导航系统中由于观测数据数量大，存在冗余现象，需要设计相应的信息融合算法对所有观测数据进行融合以获得最优状态估计值. 根据对多传感器观测信息的处理方法不同，可以将信息融合系统的体系结构分成2种，分别是集中式融合结构、分布式融合结构和混合式融合结构. 集中式融合算法是将所有观测信息整合，然后将其共同作为观测输入量，用于进行状态的量测更新. 该方法已经比较成熟，也有较多研究成果^[7-9]. 集中式融合算法具有估计精度高的特点，但是随着观测量的增加，观测矩阵维度增大，会带来巨大的计算负担，因此，近些年，分布式融合算法得到了学者们的广泛关注. 分布式融合算法对每个传感的观测信息进行单独处理，然后将处理后的结果通过适当的方法融合. 常用的融合方式有线性加权融合^[10-11]、协方差交互融合^[12-13]等，混合式融合算法是将集中式融合和分布式融合相结合的方法, 可以继承前2种融合算法的优势，但计算负担大，实时性低.

在航天任务（例如：轨道的确定和自主导航任务）中，通常存在非高斯噪声^[14]. Student’s t分布噪声作为非高斯噪声的一种特殊形式，由于其属于广义高斯分布，更适合于模拟非高斯噪声^[15]. 针对Student’s t分布噪声，Huang等^[16]提出了一个基于Student’s t的滤波框架. 与高斯滤波（Gaussian filter，GF）框架类似，Student’s t的滤波框架需要计算Student’s t的加权积分函数. 因此，Huang等^[17-18]又提出了采用不同算法计算加权积分的多种滤波算法. 此外，Roth等^[19]针对过程噪声和观测噪声均是Student’s t分布的问题，提出了基于Student’s t的鲁棒滤波器.

综上所述，在导航任务中，如果传感器的观测噪声为非高斯噪声，则需要导航系统中的信息融合算法对噪声具有一定的鲁棒性. 本文考虑系统噪声为Student’s t分布的特殊情况，以已有的Student’s t的滤波框架为基础，利用线型加权融合的方法，设计对Student’s t分布噪声具有鲁棒性的多传感器信息融合算法，以提高组合导航系统的导航精度和对特殊情况的适应能力.本文基于组合导航系统的运动模型和观测模型，选取观测噪声为满足Student’s t分布的噪声，设计基于Student’s t分布的无迹四元数局部滤波算法；设计最优权重计算方法及线型加权信息融合算法；利用仿真分析验证所设计算法的有效性和鲁棒性.

考虑航天器运动系统存在Student’s t分布噪声的情况，采用四元数描述的航天器运动学模型如下：

(1) ${\dot {{q}}} = \frac{1}{2}{\left[ {{{\omega }},0} \right]^{\rm{T}}} \otimes {{q}} = \frac{1}{2}{{\varOmega }}({{\omega }}){{q}}.$

式中： ${{q}}$为四元数向量； ${{\omega }} = {\left[ {\begin{array}{*{20}{c}} {{\omega _1},}&{{\omega _2},}&{{\omega _3}} \end{array}} \right]^{\rm{T}}}$为相对惯性系下的角速度在体坐标系下的表示；

(2) ${{\bf\varOmega }}({{\omega }}) = \left[ {\begin{array}{*{20}{c}} { - \left[ {{{\omega }} \times } \right]}&{{\omega }} \\ { - {{{\omega }}^{\rm{T}}}}&0 \end{array}} \right],$

其中， $\left[ {{{\omega }} \times } \right]$为 ${{\omega }}$的反对称矩阵，

(3) $\left[ {{{\omega }} \times } \right] = \left[ {\begin{array}{*{20}{c}} 0&{ - {\omega _3}}&{{\omega _2}} \\ {{\omega _3}}&0&{ - {\omega _1}} \\ { - {\omega _2}}&{{\omega _1}}&0 \end{array}} \right];$

四元数乘法 $ \otimes $采用Shuster的定义^[20]，

(4) ${{q}}'' = {{{{q}}' \otimes {{q}}}} = \left[ {\begin{array}{*{20}{c}} {{{q'\!\!}_4}{{\rho }} + {q_4}{{\rho '}} - {{\rho '}} \times {{\rho }}} \\ {{{q'\!\!}_4}{q_4} - {{{{\rho '}}}^{\rm{T}}}{{\rho }}} \end{array}} \right].$

本节所研究的系统为含有Student’s t分布噪声的惯导和多个星敏感器组合导航系统，导航系统框图具体如图1所示.

假设陀螺安装在航天器上，且安装方向与航天器本体坐标系重合，则陀螺的输出模型表示如下：

(5) ${{\tilde \omega }} = {{\omega }} + {{\beta }} + {{{\eta }}_v},\;\;{{\dot \beta }} = {{{\eta }}_u}.$

式中： ${{\tilde \omega }}$为陀螺实际输出； ${{\beta }}$为陀螺漂移； ${{{\eta }}_v}$和 ${{{\eta }}_u}$为Student’s t分布噪声，即

(6) $p\left( {{{{\eta }}_v}} \right) = {\rm{St}} \left( {{{{\eta }}_v};0,{{{Q}}_v},{v_v}} \right),$

(7) $p\left( {{{{\eta }}_u}} \right) = {\rm{St}} \left( {{{{\eta }}_u};0,{{{Q}}_u},{v_u}} \right),$

其中， ${{{Q}}_v}$和 ${{{Q}}_u}$为方差， ${v_v}$和 ${v_u}$为自由度，v为漂移误差，u为漂移速率坡度。.

如果星敏感器的安装方向与航天器本体坐标系完全重合，则星光矢量在体系下的观测方程为

(8) ${{b}} = {{A}}_I^b\left( {{q}} \right){{r}} + {{l}}.$

式中： ${{b}}$为星敏感器输出向量； ${{r}}$为星光矢量在惯性坐标系下的单位矢量方向，可通过星历表查询； ${{A}}_I^b$为惯性系到星体系的坐标转换矩阵（I 为惯性系，b 为星体坐标系； ${{l}}$为敏感器的矢量测量误差，假设该噪声为Student’s t分布噪声，即

$p\left( {{{{l}}^{\left( i \right)}}} \right) = {\rm{St}} \left( {{{{l}}^{\left( i \right)}};0,{{{R}}^{\left( i \right)}},v_l^{\left( i \right)}} \right).$

其中， ${{{R}}^{\left( i \right)}}$为第 $i$个观测噪声方差， $v_l^{\left( i \right)}$为自由度.

在 $k$时刻，基于四元数的矢量观测模型为

(9) $\begin{split} {{{z}}_k} = & \left[ {\begin{array}{*{20}{c}} {{{{b}}_1}},&{{{{b}}_2}},& \cdots ,&{{{{b}}_{\rm o}}} \end{array}} \right]_k^{\rm{T}} + \left[ {\begin{array}{*{20}{c}} {{{{v}}_1}},&{{{{v}}_2}},& \cdots , &{{{{v}}_{\rm o}}} \end{array}} \right]_k^{\rm{T}}= \\ & \left[ {\begin{array}{*{20}{c}} {{{A}}({{q}}){{{r}}_1}},&{{{A}}({{q}}){{{r}}_2}},& \cdots ,&{{{A}}({{q}}){{{r}}_o}} \end{array}} \right]_k^{\rm{T}} + {{l}} . \end{split}\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!$

式中： ${{{b}}_{\rm o}}$和 ${{{r}}_{\rm o}}$分别为第 ${\rm o}$个参考矢量在体系和惯性系下的分量， ${{A}}({{q}})$为姿态转移矩阵，其四元数形式^[21]如下：

(10) ${{A}}\left( {{q}} \right) = \left( {q_4^2 - {{\left\| {{\rho }} \right\|}^2}} \right){{{I}}_{3 \times 3}} - 2{q_4}\left[ {{{\rho }} \times } \right] + 2{{\rho }}{{{\rho }}^{\rm{T}}}.$

假设1　系统的初始状态 ${{{x}}_0}$和系统过程噪声 ${{{\eta }}_v}$和 ${{{\eta }}_u}$及各个观测噪声 ${{l}}_k^{\left( i \right)}$之间互不相关，且初始状态 ${{{x}}_0}$满足Student’s t分布，其概率密度函数（probability density function，PDF）如下：

(12) $p\left( {{{{x}}_0}} \right) = {\rm{St}} \left( {{{{x}}_0};{{\hat{ {x}}}_{0|0}},{{{P}}_{0|0}},{v_x}} \right).$

针对系统中存在Student’s t分布噪声的情况，文献[16]中提出基于Student’s t分布的滤波框架. 为方便设计，提出以下2个假设.

假设2　系统状态和观测的联合预测PDF满足Student’s t分布，即

(13) $\begin{array}{l} \!\!\!\!\!\!\!\!\!\!\!\! \! p\left( {{{{x}}_k},{{{z}}_k}|{{{Z}}_{k - 1}}} \right) = \\\!\!\!\!\!\!\!\!\!\!\!\! \! {\rm{St}} \left( {\left[ {\begin{array}{*{20}{c}} {{{{x}}_k}} \\ {{{{z}}_k}} \end{array}} \right];\left[ {\begin{array}{*{20}{c}} {{{\hat{{ x}}}_{k|k - 1}}} \\ {{{\hat{{ z}}}_{k|k - 1}}} \end{array}} \right],\left[ {\begin{array}{*{20}{c}} {{{{P}}_{k|k - 1}}}&{{{P}}_{k|k - 1}^{xz}} \\ {{{\left( {{{P}}_{k|k - 1}^{xz}} \right)}^{\rm{T}}}}&{{{P}}_{k|k - 1}^{zz}} \end{array}} \right],{v_x}} \right) . \end{array} $

假设3　当前状态和一步预测状态向量的联合PDF满足Student’s t分布，即

(14) $\begin{array}{l} \!\!\!\!\!\!\!\!\! p\left( {{{{x}}_k},{{{x}}_{k{\rm{ + }}1}}|{{{Z}}_k}} \right) = \\ \!\!\!\!\!\!\!\!\! {\rm{St}} \left( {\left[ {\begin{array}{*{20}{c}} {{{{x}}_k}} \\ {{{{x}}_{k{\rm{ + }}1}}} \end{array}} \right];\left[ {\begin{array}{*{20}{c}} {{{\hat{{ x}}}_{k|k}}} \\ {{{\hat{{ x}}}_{k{\rm{ + }}1|k}}} \end{array}} \right],\left[ {\begin{array}{*{20}{c}} {{{{P}}_{k|k}}}&{{{{P}}_{k,k + 1|k}}} \\ {{{P}}_{k,k + 1|k}^{\rm{T}}}&{{{{P}}_{k + 1|k}}} \end{array}} \right],{v_x}} \right). \end{array} $

基于假设2和3设计滤波算法如下.

1）系统状态一步预测及相关协方差计算.

(15) $\begin{split} {{\hat{{ x}}}_{k|k - 1}} = & \int_{{{\bf{R}}^n}} {{{{A}}_{k - 1}}{{{x}}_{k - 1}}} \times \\ & {\rm{St}} \left( {{{{x}}_{k - 1}};{{\hat{{ x}}}_{k - 1|k - 1}},{{{P}}_{k - 1|k - 1}},{v_3}} \right){\rm{d}} {{{x}}_{k - 1}}, \end{split} $

(16) $\begin{split} {{{P}}_{k|k - 1}} = & \frac{{{v_x} - 2}}{{{v_x}}}\int_{{{\bf{R}}^n}} {{{{A}}_{k - 1}}{{{x}}_{k - 1}}} {{x}}_{k - 1}^{\rm{T}}{{A}}_{k - 1}^{\rm{T}} \times \\ & {\rm{St}} \left( {{{{x}}_{k - 1}};{{\hat{{ x}}}_{k - 1|k - 1}},{{{P}}_{k - 1|k - 1}},{v_x}} \right){\rm{d}} {{{x}}_{k - 1}} - \\ & \frac{{{v_x} - 2}}{{{v_x}}}{{\hat{{ x}}}_{k|k - 1}}{\hat{{ x}}}_{k|k - 1}^{\rm{T}}{\rm{ + }}\frac{{{v_\eta }\left( {{v_x} - 2} \right)}}{{{v_x}\left( {{v_\eta } - 2} \right)}}{{{Q}}_{k - 1}}. \end{split} $

式中： ${\hat{{ x}}_{k|k - 1}}$和 ${{{P}}_{k|k - 1}}$为一步预测的状态向量和相应的协方差矩阵， ${{{Q}}_{k - 1}} = {\rm{diag}}\left( {{{{Q}}_v},{{{Q}}_u}} \right)$为 $k - 1$时刻惯导系统误差的方差矩阵， ${v_\eta } = {v_u} + {v_v}$.

2)系统状态的量测更新值及协方差计算.

(17) ${{\mathit{\Delta}} _k} = {\left[{{{\left( {{{{z}}_k} - {{\hat{{ z}}}_{k|k - 1}}} \right)}^{\rm{T}}}{{\left( {{{P}}_{k|k - 1}^{zz}} \right)}^{ - 1}}\left( {{{{z}}_k} - {{\hat{{ z}}}_{k|k - 1}}} \right)}\right]}^{1/2},\;\; \quad$

(18) ${{{K}}_k} = {{P}}_{k|k - 1}^{xz}{\left( {{{P}}_{k|k - 1}^{zz}} \right)^{ - 1}},\qquad\quad\qquad\qquad\qquad\;\;\;\;\;\;$

(19) ${{\hat{{ x}}}_{k|k}} = {{\hat{{ x}}}_{k|k - 1}} + {{{K}}_k}\left( {{{{z}}_k} - {{\hat{{ z}}}_{k|k - 1}}} \right).\quad\qquad\qquad\qquad\qquad$

(20) ${{{P}}_{k|k}} = \frac{{\left( {{v_x} - 2} \right)\left( {{v_x} + {\mathit{\Delta}} _k^2} \right)}}{{{v_x}\left( {{v_x} + m - 2} \right)}}\left( {{{{P}}_{k|k - 1}} - {{{K}}_k}{{P}}_{k|k - 1}^{zz}{{K}}_k^{\rm{T}} } \right).$

式中： $m$为观测向量维度， ${{\hat{{ x}}}_{k|k}}$和 ${{{P}}_{k|k}}$为后验估计值和相应的协方差矩阵， ${{\hat{{ z}}}_{k|k - 1}}$和 ${{P}}_{k|k - 1}^{zz}$为一步预测的观测向量和相应的协方差矩阵， ${{P}}_{k|k - 1}^{xz}$为状态和观测向量的互协方差矩阵.

(21) ${{\hat{{ z}}}_{k|k - 1}} = \int_{{{\rm{R}}^n}} {{{{H}}_k}{{{x}}_{k - 1}}} {\rm{St}} \left( {{{{x}}_k};{{\hat{{ x}}}_{k|k - 1}},{{{P}}_{k|k - 1}},{v_x}} \right){\rm{d}} {{{x}}_k},$

(22) $\begin{split} {{P}}_{k|k - 1}^{zz} = & \frac{{{v_x} - 2}}{{{v_x}}}\int_{{{\rm{R}}^n}} {{{{H}}_k}{{{x}}_k}} {{x}}_k^{\rm{T}}{{H}}_k^{\rm{T}} \times \\ & {\rm{St}} \left( {{{{x}}_k};{{\hat{{ x}}}_{k|k - 1}},{{{P}}_{k|k - 1}},{v_x}} \right){\rm{d}} {{{x}}_k} - \\ & \frac{{{v_x} - 2}}{{{v_x}}}{{\hat{{ z}}}_{k|k - 1}}{\hat{{ z}}}_{k|k - 1}^{\rm{T}}{\rm{ + }}\frac{{{v_l}\left( {{v_x} - 2} \right)}}{{\left( {{v_l} - 2} \right){v_x}}}{{{R}}_k}, \end{split} $

(23) $\begin{split} {{P}}_{k|k - 1}^{xz} = & \frac{{{v_x} - 2}}{{{v_x}}}\int_{{{\rm{R}}^n}} {{{{x}}_k}} {{x}}_k^{\rm{T}}{{H}}_k^{\rm{T}} \times \\ & {\rm{St}} \left( {{{{x}}_k};{{\hat{{ x}}}_{k|k - 1}},{{{P}}_{k|k - 1}},{v_x}} \right){\rm{d}} {{{x}}_k} - \\ & \frac{{{v_x} - 2}}{{{v_x}}}{{\hat{{ x}}}_{k|k - 1}}{\hat{{ z}}}_{k|k - 1}^{\rm{T}}. \end{split}\;\;\;\;\;\;\;\; $

假设初始量 ${{{x}}_0} = {\left[ {\begin{array}{*{20}{c}} {\delta {{p}}_0^{\rm{T}} ,}&{{\hat{{ \beta} }}_0^{\rm{T}} } \end{array}} \right]^{\rm{T}} }$，其中， $\delta {{p}}$为误差修正罗德里格参数（error modified Rodrigues parameters，EMRPs），可以有效避免四元数自身的单位特性的限制. ${{{P}}_0}$、 ${{\hat{{ q}}}_0}$均已知，且 $k - 1$时刻的状态量估计值 ${{\hat{{ x}}}_{k - 1}}$和 ${{{P}}_{k - 1}}$已获取，利用无迹变换法则，可生成Sigma点如下：

(24) $\!\!\!\!\!\!\!\!\!\!\!\begin{split}\!\!\! &{{{X}}_{i,k - 1}} = \left\{ {\begin{aligned} & {{{\hat{{ x}}}_{k - 1}},} \\ & {{{\hat{{ x}}}_{k - 1}} + {{\left[ {\frac{{{v_3}\left( {n + \kappa } \right)}}{{{v_3} - 2}}{{{P}}_{k - 1}}{{{e}}_i}} \right]}^{{\rm{1/2}}}},\;i = 1,{\rm{2}}, \cdots ,n;} \\ & {{{\hat{{ x}}}_{k - 1}} \!-\! {{\left[ {\frac{{{v_3}\left( {n + \kappa } \right)}}{{{v_3} - 2}}{{{P}}_{k - 1}}{{{e}}_i}} \right]}^{{\rm{1/2}}}},\;i \!=\! n \!+\! 1,n \!+\! {\rm{2}},\! \cdots\! ,2n.} \end{aligned}} \right. \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\\ & {w_i} = \left\{ {\begin{aligned} & {{\kappa / {\left( {n + \kappa } \right),\;i = 0;}}} \\ & {{{0.5} / {\left( {n + \kappa } \right){\rm{,}}\;i = 1,{\rm{2}}, \cdots ,n;{\rm{ }}}}} \\ & {{{0.5} / {\left( {n + \kappa } \right),\;i = n + 1,n + {\rm{2}}, \cdots ,2n.}}} \end{aligned}} \right. \\ \end{split} $

式中： ${{{X}}_{i,k - 1}} = {\left[ {\begin{array}{*{20}{c}} {\delta {{p}}_{i,k - 1}^{\rm{T}} ,}&{{\hat{{ \beta }}}_{i,k - 1}^{\rm{T}} } \end{array}} \right]^{\rm{T}} }$， ${w_i}$为无迹变换的权重系数， $\kappa $为自由参数.

1）系统状态一步预测及相关协方差计算.

由于EMRP在描述姿态时，容易产生奇异问题，通过式（25）将EMRPs转化为误差四元数 $\delta {{{q}}_{i,k - 1}}$:

(25) $\left. {\begin{split} & {\delta {q_4} = \frac{{ - a{{\left\| {\delta {{p}}} \right\|}^2} + b{{\left( {{b^2} + \left( {1 - {a^2}} \right){{\left\| {\delta {{p}}} \right\|}^2}} \right)}^{{1 / 2}}}}}{{{b^2} + {{\left\| {\delta {{p}}} \right\|}^2}}},}\\ & {\delta {{\rho }} = {b^{ - 1}}\left( {a + \delta {q_4}} \right)\delta {{p}}.} \end{split}} \right\}$

由 $\delta {{{q}}_{i,k - 1}}$和 ${{\hat{{ q}}}_{k - 1}}$可得

(26) ${{\hat{{ q}}}_{i,k - 1}} = \delta {{{q}}_{i,k - 1}} \otimes {{\hat{{ q}}}_{k - 1}}.$

通过式（5）对Sigma点传播，可得 ${{\hat{{ q}}}_{i,k{\rm{|}}k - 1}}$，式中角速率可通过下式获得：

(27) ${{{\hat \omega }}_{i,k - 1}} = {{{\tilde \omega }}_{k - 1}} - {{{\beta }}_{i,k - 1}}{\rm{,\;}}i = 0,2, \cdots ,2n.$

由 ${\hat{{ q}}}_{i,k}^ - $和 ${{\hat{{ q}}}_{k - 1}}$可得一步预测的误差四元数如下：

(28) $\delta {{{q}}_{i,k{\rm{|}}k - 1}} = {{\hat{{ q}}}_{i,k{\rm{|}}k - 1}} \otimes {\hat{{ q}}}_{k - 1}^{ - 1}.$

通过式（29）将误差四元数转换为EMRPs，则系统一步预测值及协方差如下：

(29) $\delta {{p}} = b\frac{{\delta {{\rho }}}}{{a + \delta {q_4}}}.$

(30) $\delta {{\hat{{ p}}}_{k|k - 1}} = \sum\nolimits_{i = 0}^{2n} {{w_i}} \delta {{{p}}_{i,k{\rm{|}}k - 1}}.$

(31) $\begin{split} {{{P}}_{p,k|k - 1}} = & \frac{{{v_x} - 2}}{{{v_x}}}\left[ {\sum\nolimits_{i = 0}^{2n} {{w_i}\delta {{{p}}_{i,k{\rm{|}}k - 1}}{{\left( {\delta {{{p}}_{i,k{\rm{|}}k - 1}}} \right)}^{\rm{T}}} - } } \right. \\ & \left. {\delta {{\hat {{p}}}_{k|k - 1}}\delta {{\hat{{ p}}}_{k|k - 1}} + \frac{{{v_\eta }}}{{{v_\eta } - 2}}{{{Q}}_{k - 1}}} \right]. \end{split} $

则一步预测状态为 ${{\hat{{ x}}}_{k|k - 1}} = {\left[ {\begin{array}{*{20}{c}} {\delta {{p}}_{k|k - 1}^{\rm{T}} ,}&{{\hat{{\beta }}}_{k - 1}^{\rm{T}} } \end{array}} \right]^{\rm{T}} }$，协方差矩阵为

${{{P}}_{k|k - 1}} = \left[ {\begin{array}{*{20}{c}} {{{{P}}_{p,k|k - 1}}}&0 \\ 0&{{{{P}}_{\beta ,k - 1}}} \end{array}} \right].$

2）系统状态的量测更新值及协方差计算.

以第 $j$个星敏感器为例，生成的观测Sigma传播点形式如下：

(32) ${{Z}}_{i,k{\rm{|}}k - 1}^{\left( j \right)} = {{{A}}^{\left( j \right)}}\left( {{{{\hat{{ q}}}}_{i,k{\rm{|}}k - 1}}} \right){{{r}}^{\left( j \right)}},\;{\rm{ }}i = 0,2, \cdots ,2n.$

式中： ${{{A}}^{\left( j \right)}}\left( \right)$的形式由式（11）获得，观测值的一步预测形式为

(33) ${\hat{{ z}}}_{k{\rm{|}}k - 1}^{\left( j \right)} = \sum\nolimits_{i = 0}^{2n} {{w_i}{{Z}}_{i,k{\rm{|}}k - 1}^{\left( j \right)}} .$

根据一步预测状态 ${{\hat{{ x}}}_{k|k - 1}}$ 和 ${{{P}}_{k|k - 1}}$ ，生成Sigma点如下：

(34) $\left. \begin{split} & {{{X}}_{i,k|k - 1}} = {{\hat{{ x}}}_{k|k - 1}},i = 0;\\ & {{{X}}_{i,k|k - 1}} = {{\hat{{ x}}}_{k|k - 1}} + {\left( {\frac{{{v_3}\left( {n + \kappa } \right)}}{{{v_3} - 2}}{{{P}}_{k|k - 1}}{{{e}}_i}} \right)^{{1 / 2}}},\\ & \qquad\quad i = 1,{\rm{2,}} \cdots ,n;\\ & {{{X}}_{i,k|k - 1}}= {{\hat{{ x}}}_{k|k - 1}} - {\left( {\frac{{{v_3}\left( {n + \kappa } \right)}}{{{v_3} - 2}}{{{P}}_{k|k - 1}}{{{e}}_i}} \right)^{{1 / 2}}},\\ & \qquad\quad i = n + 1,n + {\rm{2,}} \cdots ,2n. \end{split} \right\}$

然后，由Student’s t滤波框架计算观测协方差 ${{P}}_{yy,k|k - 1}^{\left( j \right)}$和互协方差 ${{P}}_{xy,k|k - 1}^{\left( j \right)}$如下：

(35) $\begin{split} {{P}}_{yy,k|k - 1}^{\left( j \right)} =& \frac{{{v_x} - 2}}{{{v_x}}}\left( {\sum\nolimits_{i = 0}^{2n} {{w_i}{{Z}}_{i,k{\rm{|}}k - 1}^{\left( j \right)}{{\left( {{{Z}}_{i,k{\rm{|}}k - 1}^{\left( j \right)}} \right)}^{\rm{T}}}} } \right. - \\ & \left. {{\hat{{ z}}}_{k|k - 1}^{\left( j \right)}{{\left( {{\hat{{ z}}}_{k|k - 1}^{\left( j \right)}} \right)}^{\rm{T}}} + \frac{{v_l^{\left( j \right)}}}{{v_l^{\left( j \right)} - 2}}{{R}}_k^{\left( j \right)}} \right). \end{split}$

(36) $\begin{split} {{P}}_{xy,k|k - 1}^{\left( j \right)} =& \frac{{{v_x} - 2}}{{{v_x}}}\left[ {\sum\nolimits_{i = 0}^{2n} {{w_i}{{{X}}_{i,k|k - 1}}{{\left( {{{Z}}_{i,k{\rm{|}}k - 1}^{\left( j \right)}} \right)}^{\rm{T}}}} } \right. - \\& \left. {{{\hat{{x}}}_{k|k - 1}}{{\left( {{\hat{{ z}}}_{k|k - 1}^{\left( j \right)}} \right)}^{\rm{T}}}} \right]. \end{split}$

式中： ${{R}}_k^{\left( j \right)}$为第 $k$时刻第 $j$个星敏感器的观测误差的方差矩阵.

最后，计算滤波增益矩阵 ${{K}}_k^{\left( j \right)}$、 后验估计值 ${{\hat{{ x}}}_{k|k}}$及相关协方差矩阵 ${{{P}}_{k|k}}$：

(37) ${{K}}_k^{\left( j \right)} = {{P}}_{xy,k|k - 1}^{\left( j \right)}{\left( {{{P}}_{yy,k|k - 1}^{\left( j \right)}} \right)^{ - 1}},$

(38) ${\hat{{ x}}}_{k|k}^{\left( j \right)} = {{\hat{{ x}}}_{k|k - 1}} + {{K}}_k^{\left( j \right)}\left( {{{z}}_k^{\left( j \right)} - {\hat{{ z}}}_{k|k - 1}^{\left( j \right)}} \right),$

(39) $ {{P}}_{k|k}^{\left( j \right)} = \frac{{\left( {{v_x} - 2} \right)\left( {{v_x} + {{\left( {\varDelta _k^{\left( j \right)}} \right)}^2}} \right)}}{{{v_x}\left( {{v_x} + m - 2} \right)}} \times \left( {{{{P}}_{k|k - 1}} - {{K}}_k^{\left( j \right)}{{P}}_{yy,k|k - 1}^{\left( j \right)}{{\left( {{{K}}_k^{\left( j \right)}} \right)}^{\rm{T}}}} \right). $

式中： $\varDelta _k^{\left( j \right)} = {\left[{{{\left( {{{z}}_k^{\left( j \right)}\! -\! {\hat{{ z}}}_{k|k \!- \!1}^{\left( j \right)}} \right)}^{\rm{T}}}{{\left( {{{P}}_{yy,k|k \!-\! 1}^{\left( j \right)}} \right)}^{\! -\! 1}}\left( {{{z}}_k^{\left( j \right)} \!-\! {\hat{{ z}}}_{k|k \!-\! 1}^{\left( j \right)}} \right)}\right]} ^{1/2}$.

由式（38）求得 ${\hat{{ x}}}_{k|k}^{\left( j \right)} = {\left[ {\begin{array}{*{20}{c}}{{{\left( {\delta {{p}}_{k|k}^{\left( j \right)}} \right)}^{\mathop{\rm T}\nolimits} },}&{{{\left( {{{\hat \beta }}_k^{\left( j \right)}} \right)}^{\mathop{\rm T}\nolimits} }}\end{array}} \right]^{\mathop{\rm T}\nolimits} }$，则可将 $\delta {{p}}_{k|k}^{\left( j \right)}$代入式（25）求得误差四元数 $\delta {{q}}_{k|k}^{\left( j \right)}$，接着计算 $k$时刻的四元数估计值 ${\hat{{ q}}}_k^{\left( j \right)}$如下：

(40) ${\hat{{ q}}}_k^{\left( j \right)} = \delta {{q}}_{k|k}^{\left( j \right)} \otimes {{\hat{{ q}}}_{k - 1}}.$

当多个星敏感器对同一目标进行测量时，会获得多个观测值 ${{y}}_k^{\left( e \right)}$， $e = 1,2,\cdots,o$. 根据每个观测值进行局部滤波可得到 $o$个估计值 ${\hat{{ x}}}_{k|k}^{\left( e \right)}$，采用分布式融合算法对所有估计值进行线性加权融合. 因此，权重系数的选取是分布式融合算法精确估计的关键. 本节基于最小线性方差原则，对融合权重系数进行选取.

信息融合后估计结果可表示为

(41) ${{\hat{{ x}}}_{k|k}} = {a_1}{\hat{{ x}}}_{k|k}^{\left( 1 \right)} + {a_2}{\hat{{ x}}}_{k|k}^{\left( 2 \right)} + \cdots + {a_o}{\hat{{ x}}}_{k|k}^{\left( o \right)}.$

式中： ${a_i}\left( {i = 1,2,\cdots,o} \right)$为权重因子，通过使代价函数 ${{J}}$最小，可计算权重因子的值. 代价函数的形式如下：

(42) ${{J}} = {\rm{tr}} \left( {{P}} \right).$

其中， ${\mathop{\rm tr}\nolimits} \left( \right)$表示矩阵的迹， ${{P}} = E\left[ {{\tilde{{ x}}}{{\tilde{{ x}}}^{\rm{T}}}} \right]$为估计误差的协方差矩阵，其中估计误差

(43) ${\tilde{ x}} = {{{x}}_k} - {{\hat{{ x}}}_{k|k}}.$

为保证融合结果的无偏性，则权重因子应满足：

(44) ${a_1} + {a_2} + \cdots +{a_o} = 1.$

由式（41）可得 $k$时刻的融合估计误差如下：

(45) ${{\tilde{ x}}_{k|k}} = {{{x}}_k} - {{\hat{{ x}}}_{k|k}} = \sum\limits_{i = 1}^o {{a_e}\left( {{{{x}}_k} - {\hat{{ x}}}_{k|k}^{\left( e \right)}} \right)} = \sum\limits_{e = 1}^o {{a_i}{\tilde{ x}}_{k|k}^{\left( e \right)}}. $

由于当 $i \ne j$时，任意 ${\tilde{ x}}_{k|k}^{\left( i \right)}$和 ${\tilde{ x}}_{k|k}^{\left( j \right)}$不相关，因此由 ${{{P}}_{k|k}}$定义可得

(46) ${{{P}}_{k|k}} = E\left[ {{{{\tilde{ x}}}_{k|k}}{{\left( {{{{\tilde{ x}}}_{k|k}}} \right)}^{\rm{T}}}} \right] = \sum\nolimits_{e = 1}^o {a_e^2{{P}}_{k|k}^{\left( e \right)}} .$

则将 ${{{P}}_{k|k}}$代入式（42）可得

(47) ${{J}} = \sum\nolimits_{e = 1}^o {a_e^2 {\rm{tr}} \left( {{{P}}_{k|k}^{\left( e \right)}} \right)} .$

需要在式（44）为限制条件下，获得使 ${{J}}$最小的权重因子 ${a_i}$，因此，根据拉格朗日乘子法，可得辅助函数如下：

(48) $F = \sum\nolimits_{e = 1}^o {a_e^2{\rm{tr}} \left( {{{P}}_{k|k}^{\left( e \right)}} \right)} + \lambda \left( {\sum\nolimits_{e = 1}^o {{a_e}} - 1} \right).$

令

(49) $\frac{{\partial F}}{{\partial {a_e}}} = 2{a_e}{\rm{tr}} \left( {{{P}}_{k|k}^{\left( e \right)}} \right) + \lambda = 0,\;i = 1,2,\cdots,o.$

式（49）等式两侧同乘 $\prod_{{j \!= \!1},\;{j \!\ne \!i}}^o {{\mathop{\rm tr}\nolimits} \left( {{{P}}_{k|k}^{\!\left( j \right)}} \right)} ,{\rm{ }}j \!= \!1,2,\! \cdots \!,$ $o $，可得

(50) $2{a_e}\prod\limits_{e = 1}^o {{\rm{tr}}\left( {{ P}_{k|k}^{\left( e \right)}} \right)} + \lambda \;\prod\limits_{j = 1,\;j \ne i}^o {\;{\rm{tr}}\left( {{ P}_{k|k}^{\left( j \right)}} \right) = 0.} $

令 $e = 1,2,\cdots,o$对式（50）进行加和可得

(51) $\lambda = - 2\prod\limits_{e = 1}^o {{\rm{tr}}\left( {{ P}_{k|k}^{\left( e \right)}} \right)} {\left[ {\sum\limits_{e = 1}^o {\;\;\prod\limits_{j = 1,\;j \ne i}^o {\;{\rm{tr}}\left( {{ P}_{k|k}^{\left( j \right)}} \right)} } } \right]^{ - 1}}.$

将式（51）代入式（49）可得

(52) ${a_e} = \frac{1}{{{\rm{tr}}\left( {{{P}}_{k|k}^{\left( e \right)}} \right)}}{\left[ {\sum\limits_{e = 1}^o {\frac{1}{{{\rm{tr}}\left( {{{P}}_{k|k}^{\left( e \right)}} \right)}}} } \right]^{ - 1}}.$

定义变量如下：

(53) $\zeta = {\left( {\frac{1}{{{\rm{tr}}\left( {{{P}}_{k|k}^{\left( 1 \right)}} \right)}} + \frac{1}{{{\rm{tr}}\left( {{{P}}_{k|k}^{\left( 2 \right)}} \right)}} + \cdots \frac{1}{{{\rm{tr}}\left( {{{P}}_{k|k}^{\left( o \right)}} \right)}}} \right)^{ - 1}}.$

则式（52）可改写为

(54) ${a_e} ={\zeta }/{{{\rm{tr}}\left( {{{P}}_{k|k}^{\left( e \right)}} \right)}}.$

最后，将式（54）代入式（41）可得融合状态估计值如下：

(55) $\begin{split} {{\hat{{ x}}}_{k|k}} = & {a_1}{\hat{{ x}}}_{k|k}^{\left( 1 \right)} + {a_2}{\hat{{ x}}}_{k|k}^{\left( 2 \right)} + \cdots + {a_m}{\hat{{ x}}}_{k|k}^{\left( o \right)}{\rm{ }} = \\ & {\rm{ }}\frac{\zeta }{{{\rm{tr}}\left( {{{P}}_{k|k}^{\left( 1 \right)}} \right)}}{\hat{{ x}}}_{k|k}^{\left( 1 \right)} + \frac{\zeta }{{{\rm{tr}}\left( {{{P}}_{k|k}^{\left( 2 \right)}} \right)}}{\hat{{ x}}}_{k|k}^{\left( 2 \right)} + \cdots +\frac{\zeta }{{{\rm{tr}}\left( {{{P}}_{k|k}^{\left( o \right)}} \right)}}{\hat{{ x}}}_{k|k}^{\left( o \right)} \end{split} .$

将式（54）代入式（46）可得估计误差协方差矩阵如下：

(56) $\begin{split} {{{P}}_{k|k}} = & {\left( {\dfrac{\zeta }{{{\rm{tr}}\left( {{{P}}_{k|k}^{\left( 1 \right)}} \right)}}} \right)^2}{{P}}_{k|k}^{\left( 1 \right)} + {\left( {\dfrac{\zeta }{{{\rm{tr}}\left( {{{P}}_{k|k}^{\left( 2 \right)}} \right)}}} \right)^2}{{P}}_{k|k}^{\left( 2 \right)} + \cdots + \\ & {\left( {\dfrac{\zeta }{{{\rm{tr}}\left( {{{P}}_{k|k}^{\left( o \right)}} \right)}}} \right)^2}{{P}}_{k|k}^{\left( o \right)}. \\ \end{split} $

目标星轨道高度为350 km，轨道角速率为0.065 6 deg/s. 假设航天器的体坐标系和轨道坐标系是完全重合的. 陀螺输出的随机漂移标准差 ${\sigma _\upsilon }$和漂移速率斜坡标准差 ${\sigma _u}$分别为 ${\sigma _\upsilon } = 1.812 \times $ ${10^{ - 5}}{{\;\deg } / {{{\rm{s}} ^{{0.5}}}}} $和 ${\sigma _u} = 1.812 \times {10^{ - 8}}{\rm{ }}\;{{{\deg } / {\rm{s}} }^{{1.5}}}$. 陀螺采样周期设为1 s. 本文仿真时间设为5个轨道周期，为验证所提出算法在长时间的跟踪过程中，仍可以保持稳定的估计结果和较好的收敛性. 假设有3个星敏感器同时测量，采样周期均为10 s. 测量标准差分别为 ${\sigma _{{{{\xi }}^{(1)}}}} = 10''$， ${\sigma _{{{{\xi }}^{(2)}}}} = 12''$和 ${\sigma _{{{{\xi }}^{(3)}}}} = 14''$ . 初值为 $\delta {{p}}_0^{\rm{T}} = \left[ {\begin{array}{*{20}{c}} {10{\rm{ }}\deg },\;& { - 10{\rm{ }}\deg },\;& {20{\rm{ }}\deg } \end{array}} \right]$和 ${{\beta }}_0^{\rm{T}} = 0$，初始方差分别为 ${{{P}}_{\delta {p_0}}} = {\left( {5{\rm{ }}\deg } \right)^2}{{{I}}_{3 \times 3}}$ 和 ${{{P}}_{{\beta _0}}} = {\left( {1{\rm{ }}\deg } \right)^2}{{{I}}_{3 \times 3}}$ . 利用所提出的Student’s t滤波框架下的无迹四元数分布式融合滤波算法（STUDF）先对3组观测值分别进行局部滤波，然后利用所设计分布式融合算法进行融合. 为验证所提出算法的高效性，将所提出算法同文献[16]中的鲁棒Student's t 无迹滤波（robust Student's t unscented filter，RSTUF）及容积四元数分布融合滤波算法（cubature quaternion distribution fusion, CQDF）进行仿真对比。

如图2~9所示分别为目标航天器的角度估计误差和误差范数，以及惯导漂移的估计误差和误差范数. 其中，α 为航向角，θ 为俯仰角, γ 为滚转角.

从图2~7可以看出，由于受到非高斯噪声的影响，高斯滤波算法的假设条件无法得到满足，CQDF的估计误差最大，角度估计误差都无法收敛到−0.02 ~0.02. RSTUF由于对Student’s t分布噪声具有一定的鲁棒性，角度估计误差可以收敛到−2×10⁻³~2×10⁻³，但是波动较大. 这是由于单个滤波器可能因为某时刻观测误差较大的影响导致估计精度下降. STUDF算法的估计误差远小于RSTUF算法，并且曲线波动小.证明所设计算法对Student’s t分布噪声具有鲁棒性. 由于采用多传感器测量时，不同的传感器估计精度不同，多个传感器信息融合滤波，可以通过高精度估计结果弥补观测误差大造成的低精度估计结果，证明了多传感器融合具有互补性.

从图8和9可以看出，CQDF由于受到非高斯噪声的影响，角度误差2范数无法收敛到10⁻²以下. 其他2种算法的误差2范数可以达到10⁻⁴ ~10⁻²，但是STUDF算法的估计误差范数小于RSTUF算法，且收敛速度更快.证明所设计融合算法对Student’s t分布噪声具有鲁棒性，且多传感器融合可以提高系统容错性，保证估计结果稳定，不易产生较大波动. 另外，单个RSTUF的平均仿真运行时间为0.149 8 s，而所设计的融合算法STUDF的平均仿真运行时间为0.204 9 s. 融合算法耗时增加了25%，并未增加过多的计算时间，这也是分布式融合算法的优势，保证了所设计算法在组合导航系统中的实时性.

本文考虑多传感器组合导航系统存在观测噪声为Student’s t分布的情况，设计了Student’s t滤波框架下的无迹四元数融合滤波算法. 为验证所设计算法的性能，以目标的四元数姿态模型为仿真模型，将所设计融合滤波算法同已有的RSTUF及CQDF算法进行滤波性能对比. 仿真结果表明，当噪声为Student’s t分布时，传统的高斯滤波算法，即CQDF算法已无法正常运行，甚至存在滤波结果发散的情况. 单个滤波器RSTUF虽然对非高斯噪声具有鲁棒性，但仍存在由于单个传感器噪声偏大而导致的结果波动大，收敛速度慢的问题. 所设计的STUDF即对Student’s t分布噪声具有鲁棒性，可以精确估计导航系统的目标姿态信息，并利用多传感器信息融合技术的互补性和容错性，减小了估计结果的波动，并且提高了收敛速度. 仿真时间对比结果证明，所设计融合算法并未增加过多计算量，保证了导航系统的实时性. 本文假设多传感器为同步观测的理想情况，但实际的多传感器组合导航系统中还存在传感器观测频率不一致的问题，因此需要在后续研究中提出相应的解决方案.