浙江大学学报(工学版), 2022, 56(12): 2514-2522 doi: 10.3785/j.issn.1008-973X.2022.12.021

航空航天技术

多矢量信息下的星敏感器三轴旋转角求解方法

段辉,, 张志利, 周召发,, 徐泽乾, 赵芝谦, 王韶迪

火箭军工程大学 导弹工程学院,陕西 西安 710025

Three-axis rotation angle solving method of star sensor under multi-vector information

DUAN Hui,, ZHANG Zhi-li, ZHOU Zhao-fa,, XU Ze-qian, ZHAO Zhi-qian, WANG Shao-di

School of Missile Engineering, Rocket Force University of Engineering, Xi’an 710025, China

通讯作者: 周召发,男,教授, 博士. orcid.org/0000-0003-4262-9256. E-mail: 3199148797@qq.com

收稿日期: 2021-12-29  

基金资助: 航空科学基金资助项目(201808U8004)

Received: 2021-12-29  

Fund supported: 航空科学基金资助项目(201808U8004)

作者简介 About authors

段辉(1996—),男,博士生,从事组合导航研究.orcid.org/0000-0003-3886-654X.E-mail:1020423896@qq.com , E-mail:1020423896@qq.com

摘要

为了提高星敏感器的实时性能与解算精度,基于多矢量信息,提出四元数表示形式下的星敏感器三轴旋转角求解方法,从理论上对所提方法进行详细推导. 基于单个星光矢量在天球系与星敏系中的三维坐标,将方向余弦阵变换形式转变成四元数变换形式. 降次二次四元数变换形式,以便后续求解. 考虑不同星光矢量的权重,联立所有星光矢量的信息求解星敏感器三轴旋转角,给出解决三轴旋转角求解过程中的迭代不确定性和四元数矢量方向奇异性的具体方案. 将所提方法与传统方法的性能进行对比分析. 仿真结果表明,相比传统方法,所提方法具有更快的解算速度和更高的三轴旋转角求解精度.

关键词: 旋转角 ; 星敏感器 ; 四元数 ; 奇异性 ; 实时性

Abstract

In order to improve the real-time performance and solution accuracy of star sensor, based on multi-vector information, a quaternion representation method for solving three-axis rotation angle of star sensor was proposed, and the method was derived in detail theoretically. Based on the three-dimensional coordinates of a single starlight vector in celestial coordinate system and star sensor coordinate system, the direction cosine matrix transformation form was transformed into quaternion transformation form. The quadratic quaternion transformation form was reduced in order to facilitate the subsequent solution. Considering the weight of different starlight vectors, the information of all starlight vectors was combined to solve the three-axis rotation angle of the star sensor. A specific solution was given, aiming at the iteration uncertainty and the direction singularity of quaternion vector in the process of solving three-axis rotation angle, The performance of the proposed algorithm was compared with that of the traditional algorithm. Simulation results show that the proposed algorithm has faster solving speed and higher accuracy in solving three-axis rotation angle than the traditional algorithm.

Keywords: rotation angle ; star sensor ; quaternion ; singularity ; real time

PDF (989KB) 元数据 多维度评价 相关文章 导出 EndNote| Ris| Bibtex  收藏本文

本文引用格式

段辉, 张志利, 周召发, 徐泽乾, 赵芝谦, 王韶迪. 多矢量信息下的星敏感器三轴旋转角求解方法. 浙江大学学报(工学版)[J], 2022, 56(12): 2514-2522 doi:10.3785/j.issn.1008-973X.2022.12.021

DUAN Hui, ZHANG Zhi-li, ZHOU Zhao-fa, XU Ze-qian, ZHAO Zhi-qian, WANG Shao-di. Three-axis rotation angle solving method of star sensor under multi-vector information. Journal of Zhejiang University(Engineering Science)[J], 2022, 56(12): 2514-2522 doi:10.3785/j.issn.1008-973X.2022.12.021

求解星敏感器三轴旋转角主要依赖恒星星光矢量的天球系坐标、星敏系坐标. 其中矢量天球系坐标在星图识别成功后,由导航星库给出;矢量星敏系坐标由星图经星点提取后给出. 天球系坐标与对应的星敏系坐标组成该恒星的星光矢量对[1-2].

针对星敏感器的姿态确定问题,国内外学者基于不同的理论方法给出了相应的解决方案。Markley等[3]开发的奇异值分解(singular value decomposition, SVD)方法稳健且准确. 虽然SVD具有鲁棒的数值矩阵运算能力,但矩阵分解的时间消耗较大. Yang等[4]提出针对Wahba问题的解析方法,通过求解四次方程的根得到最大特征值. 该算法计算速度快,但可能会由于特征多项式没有实根导致无法得到正确的旋转四元数. Shuster等[5]提出的QUEST算法通过Cayley-Hamilton定理给出对应矩阵的闭型特征多项式. 数据矩阵的最大特征值非常接近1,因此基于牛顿算法,以1为初始值迭代,得到最大特征值. Mortari[6]提出采用ESOQ和ESOQ2方法计算最优旋转四元数. 其中ESOQ推导 $ {\boldsymbol{K}} $的闭型特征多项式( $ {\boldsymbol{K}} $为由权重、恒星星光矢量的天球系坐标与星敏系坐标构成的 ${\text{4}} \times {\text{4}}$矩阵), ESOQ2作为次优估计量,利用向量变换来迭代计算K的最大特征值,它的处理速度比ESOQ快. 由于损失函数可以完美地表示为单位四元数的二次函数,Davenport能够将Wahba问题转换为特征值-特征向量问题. 原因是对应矩阵为正定对称矩阵,因此单位四元数的二次函数为标准的凸函数,具有良好的全局极值性质,直接对单位四元数的二次函数求导即可得到最优解,与该最优解对应的就是矩阵的最小特征值. Horn等[7]指出,一般的特征值/特征向量算法计算出所有的特征值和特征向量,这些特征值对三轴旋转角求解不必要. 为了求出最大特征值及其特征向量,Horn等[7]将该问题简化为瑞利商问题进行求解. Farrell等[8]给出姿态估计损失函数的极性分解解,该解是Wahba问题的显式伪逆解,因此该方法实时性差,且鲁棒性弱.

随着我国航空航天事业以及军事装备的飞速发展,空间载体的自主导航问题成为研究热点. 相较于太阳敏感器、地平仪、磁力计等姿态测量设备,星敏感器测量精度高,抗干扰能力及自主导航能力强. 本研究就星敏感器设备的姿态解算问题给出既考虑解算精度、又保证实时性的方法. 将从四元数的角度出发,联立所有星光矢量的信息;考虑不同信息间的权重分布,将三轴旋转角的求解问题转换成求解构造出的数据矩阵的特征值特征向量问题;给出解决迭代过程中出现的不确定性和四元数向量奇异性的具体方案.

1. 星敏感器测量模型

恒星星光矢量在星敏感器电荷耦合器件(charge coupled device,CCD)上的成像过程如图1所示[9-13]. 在不考虑星敏感器测量噪声(如散粒噪声、读出噪声)的情况下,恒星星光矢量的理想映射模型为

图 1

图 1   星光矢量成像示意图

Fig.1   Schematic diagram of starlight vector imaging


$ {\boldsymbol{D}}_k^{{s}} = {\boldsymbol{C}}_{{i}}^{{s}}{\boldsymbol{D}}_k^{{i}}. $

式中: $ {\boldsymbol{C}}_{{i}}^{{s}} $为惯性坐标系i到星敏感器坐标系s的旋转矩阵; $ {\boldsymbol{D}}_k^{{i}} $为第 $k$个星光矢量在i中的坐标, $\parallel {\boldsymbol{D}}_k^{{i}}\parallel = 1$$ {\boldsymbol{D}}_k^{{s}} $为第k个星光矢量在s中的坐标, $\parallel {\boldsymbol{D}}_k^{{s}}\parallel = 1$.$k$颗恒星在CCD平面上的投影为 $p({x_k},{y_k})$,星敏感器的焦距为 $f$,星光矢量在s中的三维坐标为

$ {\boldsymbol{D}}_k^{{s}} = \left[ \begin{gathered} D_{{{x,}}k}^{{s}} \\ D_{{{y,}}k}^{{s}} \\ D_{{{z,}}k}^{{s}} \\ \end{gathered} \right] = \frac{1}{{\sqrt {{x_k}^2+{y_k}^2+{f^2}} }}\left[ \begin{gathered} - {x_k} \\ - {y_k} \\ \;\; f \\ \end{gathered} \right]. $

式中: $ {x_k} $$ {y_k} $为恒星在CCD平面上的二维坐标. $ {\boldsymbol{D}}_k^{{i}} $的信息无法从星敏感器的星图中直接获得,须经星图识别后由导航星表得到,计算式为[14-17]

$ {\boldsymbol{D}}_k^{{i}} = \left[ \begin{gathered} D_{{{x,}}k}^{{i}} \\ D_{{{y,}}k}^{{i}} \\ D_{{{z,}}k}^{{i}} \\ \end{gathered} \right] = \left[ \begin{gathered} \cos {\delta _k}\cos {\alpha _k} \\ \cos {\delta _k}\sin {\alpha _k} \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \sin {\delta _k} \\ \end{gathered} \right]. $

式中: ${\alpha _k}$${\delta _k}$分别为第 $k$颗恒星的星光矢量在天球坐标系下的赤经、赤纬. 视场内 $n$颗恒星的矢量信息对为

$\left. \begin{array}{l} {{\boldsymbol{D}}}_{k}^{s}={\left[{D}_{x,k}^{s},{D}_{y,k}^{s},{D}_{z,k}^{s}\right]}^{\text{T}}\text{,}k = 1,2,\cdots ,n \text{;}\\ {{\boldsymbol{D}}}_{k}^{i}={\left[{{{D}}}_{x,k}^{i},{D}_{y,k}^{i},{D}_{z,k}^{i}\right]}^{\text{T}}\text{,}k = 1,2,\cdots ,n . \end{array} \right\} $

可以得到

$ \left[ \begin{gathered} {\kern 1pt} D_{{{x}},1}^{{s}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} D_{{{x}},2}^{{s}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \cdots {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} D_{{{x}},n}^{{s}} \\ D_{{{y}},1}^{{s}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} D_{{{y}},2}^{{s}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \cdots {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} D_{{{y}},n}^{{s}} \\ {\kern 1pt} D_{{{z}},1}^{{s}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} D_{{{z}},2}^{{s}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \cdots {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} D_{{{z}},n}^{{s}} \\ \end{gathered} \right] = {\boldsymbol{C}}_{{i}}^{{s}}\left[ \begin{gathered} {\kern 1pt} D_{{{x}},1}^{{i}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} D_{{{x}},2}^{{i}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \cdots {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} D_{{{x}},n}^{{i}} \\ D_{{{y,}}1}^{{i}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} D_{{{y}},2}^{{i}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \cdots {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} D_{{{y}},n}^{{i}} \\ {\kern 1pt} D_{{{z}},1}^{{i}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} D_{{{z}},2}^{{i}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \cdots {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} D_{{{z}},n}^{{i}} \\ \end{gathered} \right]. $

$ \left. \begin{array}{l} {\boldsymbol{V}}=\left[{{\boldsymbol{D}}}_{1}^{{i}}{,}{{\boldsymbol{D}}}_{2}^{{i}}{,}\cdots , {{\boldsymbol{D}}}_{n}^{{i}}\right]{,}\\ {\boldsymbol{W}}=\left[{{\boldsymbol{D}}}_{1}^{{s}}{,}{{\boldsymbol{D}}}_{2}^{{s}}{,}\cdots, {{\boldsymbol{D}}}_{n}^{{s}}\right]. \end{array} \right\} $

考虑星敏感器散粒噪声、读出噪声的影响,将式(5)简写为

$ {\boldsymbol{W}} = {\boldsymbol{C}}_{{i}}^{{s}}{\boldsymbol{V}}+{\boldsymbol{\delta }}. $

式中: ${\boldsymbol{\delta }}$为总的合成噪声,通常被看作是均值为0的高斯噪声. 利用不同的三轴旋转角求解算法可以得到旋转矩阵 $ {\boldsymbol{C}}_{{i}}^{{s}} $(下文简写为 ${\boldsymbol{C}}$)和旋转四元数 $ {\boldsymbol{q}} $.

2. 三轴旋转角求解算法推导

旋转角求解的核心是利用所有的矢量信息对求解出最优旋转矩阵 ${\boldsymbol{C}}$. 方法是利用已知信息,在某种最优准则下,将噪声对旋转矩阵计算的影响降到最低。该准则为测量误差 $ {{\boldsymbol{e}}_k} $模的加权平方和最小[18-21]. 对于单个的矢量信息对,可以得到

$ \sqrt{{a}_{k}}({\boldsymbol{C}}{{\boldsymbol{D}}}_{k}^{{i}}-{{\boldsymbol{D}}}_{k}^{{s}})={{\boldsymbol{e}}}_{k}{;}k = 1,2,\cdots ,n . $

式中: $ {a_k} $为不同矢量信息对的权重; $ {{\boldsymbol{e}}_k} $为第 $ k $个残余误差向量,在理想情况下, $ {{\boldsymbol{e}}_k} = {\boldsymbol{0}} $.

Wahba问题构造的目标优化函数为

$ \begin{split} {J^*}({\boldsymbol{C}}) =& 0.5 \sum\nolimits_{k = 1}^n {{a_i}} {\left| {{\boldsymbol{D}}_k^{{s}} - {\boldsymbol{CD}}_k^{{i}}} \right|^2} = \\& 0.5 \sum\nolimits_{k = 1}^n {{{\left| {\sqrt {{a_k}} ({\boldsymbol{D}}_k^{{s}} - {\boldsymbol{CD}}_k^{{i}})} \right|}^2}} . \end{split} $

式中: $ {a_i} $为非负加权因子,由传感器的测量精度决定,满足 ${a_1}{\text{+}}{a_2}{\text{+}} \cdots {\text{+}}{a_n} = 1$. 设某一观测矢量对为

$ \left. \begin{array}{l} {{\boldsymbol{D}}}_{k}^{{s}}={\left[{D}_{{x},k}^{{s}},{D}_{{y},k}^{{s}},{D}_{{z},k}^{{s}}\right]}^{{{\rm{T}}}}{,}\\ {{\boldsymbol{D}}}_{k}^{{i}}={\left[{D}_{{x},k}^{{i}},{D}_{{y},k}^{{i}},{D}_{{z},k}^{{i}}\right]}^{{{\rm{T}}}}. \end{array} \right\} $

化简 ${\boldsymbol{CD}}_k^{{i}}$,得到

$ {\boldsymbol{CD}}_k^{{i}} = D_{{{x}},k}^{{i}}{{\boldsymbol{C}}_1}+D_{{{y}},k}^{{i}}{{\boldsymbol{C}}_2}+D_{{{z}},k}^{{i}}{{\boldsymbol{C}}_3}. $

式中: $ {{\boldsymbol{C}}}_{1}、{{\boldsymbol{C}}}_{2}、{{\boldsymbol{C}}}_{3} $${\boldsymbol{C}}$的列. ${\boldsymbol{C}}$与对应四元数 ${\boldsymbol{q}} = {\left[ {{q_0},{\boldsymbol{q}}_{\rm{v}}^{\text{T}}} \right]^{\text{T}}} = {\left[ {{q_0},{q_1},{q_2},{q_3}} \right]^{\text{T}}}$具有如下关系:

$ \left. \begin{array}{l} {\boldsymbol{C}} = ({q_0}^2 - {{\boldsymbol{q}}_{\rm{v} }^{\text{T}}}{{\boldsymbol{q}}_{\rm{v}}})I+2{{\boldsymbol{q}}_{\rm{v}}}{{\boldsymbol{q}}_{\rm{v}}^{\text{T}}}+2{q_0}({{\boldsymbol{q}}_{\rm{v}}} \times ), \\ \left(\boldsymbol{q}_{\rm{v}} \times\right)=\left[\begin{array}{ccc} 0 & -q_3 & q_2 \\ q_3 & 0 & -q_1 \\ -q_2 & q_1 & 0 \end{array}\right] . \end{array} \right\} $

展开式(12),得到

$ \begin{split} {\boldsymbol{C}} = \left[ \begin{gathered} q_0^2+q_1^2 - q_2^2 - q_3^2 {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} 2{q_1}{q_2}+{\kern 1pt} 2{q_0}{q_3} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} 2{q_1}{q_3}{\kern 1pt} - 2{q_0}{q_2} \\ 2{q_1}{q_2} - 2{q_0}{q_3} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} q_0^2 - q_1^2+q_2^2 - q_3^2{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} 2{q_2}{q_3}+{\kern 1pt} 2{q_0}{q_1} \\ 2{q_1}{q_3}{\kern 1pt} +2{q_0}{q_2} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} 2{q_2}{q_3} - 2{q_0}{q_1} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} q_0^2 - q_1^2 - q_2^2+q_3^2 \\ \end{gathered} \right] . \end{split}$

$ \begin{split} {{\boldsymbol{C}}_1} = & \left[ \begin{gathered} q_0^2+q_1^2 - q_2^2 - q_3^2 \\ 2{q_1}{q_2} - 2{q_0}{q_3} \\ 2{q_0}{q_2}+2{q_1}{q_3} \\ \end{gathered} \right] = \\ & \left[ \begin{gathered} \;\;\; {q_0}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {q_1}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} - {q_2}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} - {q_3} \\ - {q_3}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {q_2}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {q_1}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} - {q_0} \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {q_2}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {q_3}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {q_0}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {q_1} \\ \end{gathered} \right]{\kern 1pt} \left[ \begin{gathered} {q_0} \\ {q_1} \\ {q_2} \\ {q_3} \\ \end{gathered} \right] = {{\boldsymbol{P}}_1}{\boldsymbol{q}}. \end{split}$

同理,得到 $ {{\boldsymbol{C}}_2} = {{\boldsymbol{P}}_2}{\boldsymbol{q}} $$ {{\boldsymbol{C}}_3} = {{\boldsymbol{P}}_3}{\boldsymbol{q}} $,其中

$ \left. \begin{array}{l} \begin{split} & {{{{\boldsymbol{P}}}}_2} = \left[ \begin{array}{l} {q_3}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \; \; {q_2}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {q_1}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {q_0}\\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {q_0}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} - {q_1}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {q_2}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} - {q_3}\\ - {q_1}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} - {q_0}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {q_3}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \, {q_2} \end{array} \right]{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} , \\& {{{{\boldsymbol{P}}}}_3} = \left[ \begin{array}{l} - {q_2}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \;\; {q_3}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} - {q_0}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {q_1}\\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {q_1}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {q_0}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \; {q_3}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {q_2}\\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {q_0}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} - {q_1}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} - {q_2}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \;\; {q_3} \end{array} \right]. \end{split} \end{array} \right\}$

$ {{\boldsymbol{C}}}_{1}、{{\boldsymbol{C}}}_{2}、{{\boldsymbol{C}}}_{3} $代入式(11),得到

$ \begin{split} {\boldsymbol{CD}}_k^{{i}} =& D_{{{x}},k}^{{i}}{{\boldsymbol{C}}_1}+D_{{{y}},k}^{{i}}{{\boldsymbol{C}}_2}+D_{{{z}},k}^{{i}}{{\boldsymbol{C}}_3} = \\& (D_{{{x}},k}^{{i}}{{\boldsymbol{P}}_1}+D_{{{y}},k}^{{i}}{{\boldsymbol{P}}_2}+D_{{{z}},k}^{{i}}{{\boldsymbol{P}}_3}){\boldsymbol{q}} = {{\boldsymbol{P}}_{D_{k}}}{\boldsymbol{q}}{\kern 1pt} . \end{split}$

$ {{\boldsymbol{P}}_{D_{k}}} = D_{{{x}},k}^{{i}}{{\boldsymbol{P}}_1}+D_{{{y}},k}^{{i}}{{\boldsymbol{P}}_2}+D_{{{z}},k}^{{i}}{{\boldsymbol{P}}_3}. $

得到 $ \sqrt {{a_k}} {{\boldsymbol{P}}_{D_{k}}}{\boldsymbol{q}} - \sqrt {{a_k}} {\boldsymbol{D}}_k^{{s}} = {{\boldsymbol{e}}_k} $,假设此时无测量误差,即

$ \sqrt {{a_k}} {{\boldsymbol{P}}_{D_k}}{\boldsymbol{q}} = \sqrt {{a_k}} {\boldsymbol{D}}_k^{{s}}. $

解耦 $ {{\boldsymbol{P}}_{D_k}}{\boldsymbol{q}} $,分离2个矩阵变量. 可知

$ \begin{split} & \left( {\sqrt {{a_k}} {{\boldsymbol{P}}_{D_k}}} \right) {\left( {\frac{{\sqrt {{a_k}} {{\boldsymbol{P}}_{D_k}}}}{{{a_k}}}} \right)^{\rm{T}}} = {\left( {D_{{{x}},k}^{{i}}} \right)^2}{{\boldsymbol{P}}_1}{{\boldsymbol{P}}_1^{\rm{T}}} + {\left( {D_{{{y}},k}^{{i}}} \right)^2}{{\boldsymbol{P}}_2}{{\boldsymbol{P}}_2^{\rm{T}}}+ \\& {\left( {D_{{{z}},k}^{{i}}} \right)^2}{{\boldsymbol{P}}_3}{{\boldsymbol{P}}_3^{\rm{T}}}+{\kern 1pt} D_{{{x}},k}^{{i}}D_{{{y}},k}^{{i}}\left( {{{\boldsymbol{P}}_1}{{\boldsymbol{P}}_2^{\rm{T}}}+{{\boldsymbol{P}}_2}{{\boldsymbol{P}}_1^{\rm{T}}}} \right)+ \\& D_{{{x}},k}^{{i}}D_{{{z}},k}^{{i}}\left( {{{\boldsymbol{P}}_1}{{\boldsymbol{P}}_3^{\rm{T}}}+{{\boldsymbol{P}}_3}{{\boldsymbol{P}}_1^{\rm{T}}}} \right)+{\kern 1pt} {\kern 1pt} D_{{{y}},k}^{{i}}D_{{{z}},k}^{{i}}\left( {{{\boldsymbol{P}}_2}{{\boldsymbol{P}}_3^{\rm{T}}}+{{\boldsymbol{P}}_3}{{\boldsymbol{P}}_2^{\rm{T}}}} \right). \end{split} $

$ \begin{split} &\boldsymbol{P}_1 \boldsymbol{P}_2^{\mathrm{T}}+\boldsymbol{P}_2 \boldsymbol{P}_1^{\mathrm{T}}=\\& \left[\begin{array}{cccc} q_0 & q_1 & -q_2 & -q_3 \\ -q_3 & q_2 & q_1 & -q_0 \\ q_2 & q_3 & q_0 & q_1 \end{array}\right]\left[\begin{array}{ccc} q_3 & q_0 & -q_1 \\ q_2 & -q_1 & -q_0 \\ q_1 & q_2 & q_3 \\ q_0 & -q_3 & q_2 \end{array}\right]+\\& \left[\begin{array}{cccc} q_3 & q_2 & q_1 & q_0 \\ q_0 & -q_1 & q_2 & -q_3 \\ -q_1 & -q_0 & q_3 & q_2 \end{array}\right]\left[\begin{array}{ccc} q_0 & -q_3 & q_2 \\ q_1 & q_2 & q_3 \\ -q_2 & q_1 & q_0 \\ -q_3 & -q_0 & q_1 \end{array}\right]=\boldsymbol{0} . \end{split} $

同理,可知

$ \begin{array}{l}{{\boldsymbol{P}}}_{1}{{\boldsymbol{P}}}_{1}^{\text{T}}={{\boldsymbol{P}}}_{2}{{\boldsymbol{P}}}_{2}^{\text{T}}={{\boldsymbol{P}}}_{3}{{\boldsymbol{P}}}_{3}^{\text{T}}={\boldsymbol{I}}\text{,}\\ {{\boldsymbol{P}}}_{1}{{\boldsymbol{P}}}_{3}^{\text{T}}+{{\boldsymbol{P}}}_{3}{{\boldsymbol{P}}}_{1}^{\text{T}}={{\boldsymbol{P}}}_{2}{{\boldsymbol{P}}}_{3}^{\text{T}}+{{\boldsymbol{P}}}_{3}{{\boldsymbol{P}}}_{2}^{\text{T}}={\boldsymbol{0}}.\end{array} $

得到

$ \left( {\sqrt {{a_k}} {{\boldsymbol{P}}_{D_k}}} \right){\left( {\frac{{\sqrt {{a_k}} {{\boldsymbol{P}}_{D_k}}}}{{{a_k}}}} \right)^{\text{T}}} = {\boldsymbol{I}}. $

矩阵 $ \sqrt {{a_k}} {{\boldsymbol{P}}_{D_k}} $的右逆矩阵为 $ {({{(\sqrt {{a_k}} {{\boldsymbol{P}}_{D_k}})} \mathord{\left/ {\vphantom {{(\sqrt {{a_k}} {{\boldsymbol{P}}_{D_k}})} {{a_k}}}} \right. } {{a_k}}})^{\text{T}}} $. 式(18)转换为

$ \frac{{\sqrt {{a_k}} }}{{{a_k}}}{\boldsymbol{P}}_{D_k}^{\text{T}}\sqrt {{a_k}} {\boldsymbol{D}}_k^{\text{s}} = {\boldsymbol{q}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \Rightarrow {a_k}{\boldsymbol{P}}_{D_k}^{\text{T}}{\boldsymbol{D}}_k^{\text{s}} = {a_k}{\boldsymbol{q}}. $

将式(23)左侧展开, 得到

$ \begin{split} &{a_k}{\boldsymbol{P}}_{D_k}^{{{\rm{T}}}}{\boldsymbol{D}}_k^{{s}} = {a_k}\left[ {D_{{{x}},k}^{{i}}{\boldsymbol{P}}_1^{{{\rm{T}}}}+D_{{{y}},k}^{{i}}{\boldsymbol{P}}_2^{{{\rm{T}}}}+D_{{{z}},k}^{{i}}{\boldsymbol{P}}_3^{{{\rm{T}}}}} \right]{\boldsymbol{D}}_k^{{s}}= \\ &{a_k}\left[ {D_{{{x}},k}^{{i}}\left[ \begin{array}{l} {\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{q_0}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt} - {q_3}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{q_2}\\ {\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt} {q_1}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt} {q_2}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt} {q_3}\\ - {q_2}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt} {q_1}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt} {q_0}\\ - {q_3}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt} - {q_0}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt} {q_1} \end{array} \right]+D_{{{y}},k}^{{i}}\left[ \begin{array}{l} {q_3}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt} \,\,\,\, {q_0}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt} - {q_1}\\ {\kern 1pt} {q_2}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt} - {q_1}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt} - {q_0}\\ {\kern 1pt} {q_1}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{q_2}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt} {q_3}\\ {\kern 1pt} {q_0}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt} - {q_3}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{q_2} \end{array} \right]+D_{{{z}},k}^{{i}}\left[ \begin{array}{l} - {q_2}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt} {q_1}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{q_0}\\ {\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{q_3}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{q_0}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt} - {q_1}\\ - {q_0}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{q_3}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt} - {q_2}\\ {\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{q_1}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{q_2}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{\kern 1pt}{q_3} \end{array} \right]} \right]\left[ \begin{array}{l} D_{{{x}},k}^{{s}}\\ D_{{{y}},k}^{{s}}\\ D_{{{z}},k}^{{s}} \end{array} \right]= \\ &\begin{array}{l} {{a_k}\left[ {\begin{array}{*{20}{c}} {D_{x,k}^iD_{x,k}^s + D_{y,k}^iD_{y,k}^s + D_{z,k}^iD_{z,k}^s} & {D_{z,k}^iD_{y,k}^s - D_{y,k}^iD_{z,k}^s} & { - D_{z,k}^iD_{x,k}^s + D_{x,k}^iD_{z,k}^s} & {D_{y,k}^iD_{x,k}^s - D_{x,k}^iD_{y,k}^s} \\ * & {D_{x,k}^iD_{x,k}^s - D_{y,k}^iD_{y,k}^s - D_{z,k}^iD_{z,k}^s} & {D_{y,k}^iD_{x,k}^s + D_{x,k}^iD_{y,k}^s} &{D_{z,k}^iD_{x,k}^s + D_{x,k}^iD_{z,k}^s} \\ * & * & { - D_{x,k}^iD_{x,k}^s + D_{y,k}^iD_{y,k}^s - D_{z,k}^iD_{z,k}^s} &{D_{z,k}^iD_{y,k}^s + D_{y,k}^iD_{z,k}^s} \\ * & * & * & { - D_{x,k}^iD_{x,k}^s - D_{y,k}^iD_{y,k}^s + D_{z,k}^iD_{z,k}^s} \end{array}} \right] \left[ {\begin{array}{*{20}{l}} {{q_0}}\\ {{q_1}}\\ {{q_2}}\\ {{q_3}} \end{array}} \right]{{. }}} \end{array}\\[-26pt] \end{split} $

式中: $ * $为对称矩阵的对称元素. 取 $ {a_k}{\boldsymbol{P}}_{D_k}^{\text{T}}{\boldsymbol{D}}_k^{\text{s}} = {a_k}{\boldsymbol{q}} $的第1行,

$ {\left[ {\begin{array}{*{20}{c}} {a_k}\left( {D_{{ {x}},k}^{ {i}}D_{{ {x}},k}^{ {s}}+D_{{ {y}},k}^{ {i}}D_{{ {y}},k}^{ {s}}+D_{{ {z}},k}^{ {i}}D_{{ {z}},k}^{ {s}}} \right) \\ {a_k}\left( {D_{{ {z}},k}^{ {i}}{\kern 1pt} D_{{ {y}},k}^{ {s}} - D_{{ {y}},k}^{ {i}}D_{{ {z}},k}^{ {s}}} \right) \\ {a_k}\left( { - D_{{ {z}},k}^{ {i}}D_{{ {x}},k}^{ {s}}+D_{{ {x}},k}^{ {i}}D_{{ {z}},k}^{ {s}}} \right){\kern 1pt} \\ {a_k}\left( {D_{{ {y}},k}^{ {i}}D_{{ {x}},k}^{ {s}} - D_{{ {x}},k}^{ {i}}D_{{ {y}},k}^{ {s}}} \right){\kern 1pt} {\kern 1pt} {\kern 1pt} \\ \end{array}} \right]^{ {{\rm{T}}}}}\left[ \begin{gathered} {q_0} \\ {q_1} \\ {q_2} \\ {q_3} \\ \end{gathered} \right] = {a_k}{q_0}. $

$ n $个星光矢量信息对求和,得到

$ {\left[ {\begin{array}{*{20}{c}} \sum\limits_{k = 1}^n {{a_k}\left( {D_{{ {x}},k}^{ {i}}D_{{ {x}},k}^{ {s}}+D_{{ {y}},k}^{ {i}}D_{{ {y}},k}^{ {s}}+D_{{ {z}},k}^{ {i}}D_{{ {z}},k}^{ {s}}} \right)} \\ \sum\limits_{k = 1}^n {{a_k}\left( {D_{{ {z}},k}^{ {i}}{\kern 1pt} D_{{ {y}},k}^{ {s}} - D_{{ {y}},k}^{ {i}}D_{{ {z}},k}^{ {s}}} \right)} {\kern 1pt} {\kern 1pt} {\kern 1pt} \\ \sum\limits_{k = 1}^n {{a_k}\left( { - D_{{ {z}},k}^{ {i}}D_{{ {x}},k}^{ {s}}+D_{{ {x}},k}^{ {i}}D_{{ {z}},k}^{ {s}}} \right)} {\kern 1pt} {\kern 1pt} {\kern 1pt} \\ \sum\limits_{k = 1}^n {{a_k}\left( {D_{{ {y}},k}^{ {i}}D_{{ {x}},k}^{ {s}} - D_{{ {x}},k}^{ {i}}D_{{ {y}},k}^{ {s}}} \right)} {\kern 1pt} {\kern 1pt} {\kern 1pt} \\ \end{array}} \right]^{ {{\rm{T}}}}}\left[ \begin{gathered} {q_0} \\ {q_1} \\ {q_2} \\ {q_3} \\ \end{gathered} \right] = {q_0}. $

$ {a_k}{\boldsymbol{P}}_{D_k}^{\text{T}}{\boldsymbol{D}}_k^{{s}} = {a_k}{\boldsymbol{q}} $的每行求和,令

$ {\boldsymbol{E}} = \left[ \begin{gathered} \sum\limits_{k = 1}^n {{a_k}D_{{ {x}},k}^{ {i}}D_{{ {x}},k}^{ {s}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \sum\limits_{k = 1}^n {{a_k}D_{{ {x}},k}^{ {i}}D_{{ {y}},k}^{ {s}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} } \sum\limits_{k = 1}^n {{a_k}D_{{ {x}},k}^{ {i}}D_{{ {z}},k}^{ {s}}} } \\ \sum\limits_{k = 1}^n {{a_k}D_{{ {y}},k}^{ {i}}D_{{ {x}},k}^{ {s}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \sum\limits_{k = 1}^n {{a_k}D_{{ {y}},k}^{ {i}}D_{{ {y}},k}^{ {s}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} } \sum\limits_{k = 1}^n {{a_k}D_{{ {y}},k}^{ {i}}D_{{ {z}},k}^{ {s}}} } \\ \sum\limits_{k = 1}^n {{a_k}D_{{ {z}},k}^{ {i}}D_{{ {x}},k}^{ {s}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \sum\limits_{k = 1}^n {{a_k}D_{{ {z}},k}^{ {i}}{\kern 1pt} D_{{ {y}},k}^{ {s}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} } \sum\limits_{k = 1}^n {{a_k}D_{{ {z}},k}^{ {i}}D_{{ {z}},k}^{ {s}}} } \\ \end{gathered} \right] . $

得到

$ \begin{split} \\[-3pt] \left. \begin{array}{l} {\boldsymbol{ Qq}}={\boldsymbol{q}}\text{;}\\ {Q}_{11}={E}_{11}+{E}_{22}+{E}_{33}\text{,} {Q}_{31}=-{E}_{31}+{E}_{13}\text{,}\\ {Q}_{12}=-{E}_{23}+{E}_{32}\text{,} \qquad {Q}_{32}={E}_{21}+{E}_{12}\text{,}\\ {Q}_{13}={E}_{13}-{E}_{31}\text{,} \qquad \;\; {Q}_{33}=-{E}_{11}+{E}_{22}-{E}_{33}\text{,}\\ {Q}_{14}=-{E}_{12}+{E}_{21}\text{,} \qquad {Q}_{34}={E}_{32}+{E}_{23}\text{,}\\ {Q}_{21}=-{E}_{23}+{E}_{32}\text{,} \qquad {Q}_{41}={E}_{21}-{E}_{12}\text{,}\\ {Q}_{22}={E}_{11}-{E}_{22}-{E}_{33}\text{,} {Q}_{42}={E}_{31}+{E}_{13}\text{,}\\ {Q}_{23}={E}_{21}+{E}_{12}\text{,} \qquad \;\; {Q}_{43}={E}_{32}+{E}_{23}\text{,}\\ {Q}_{24}={E}_{31}+{E}_{13}\text{,} \qquad \;\; {Q}_{44}=-{E}_{11}-{E}_{22}+{E}_{33}. \end{array} \right\} \qquad \\[-70pt] \end{split} $

因为星敏感器在测量星光矢量时的系统偏差与随机噪声不可避免,所以星光矢量 ${\boldsymbol{D}}_k^{{s}} = [ D_{{{x}},k}^{{s}} , D_{{{y}},k}^{{s}}, D_{{{z}},k}^{{s}} ]^{\text{T}}$必定不是真值,因此 $ {\boldsymbol{Qq}}= {\boldsymbol{q}} $右侧须补充误差因子 $ \varepsilon $,即

$ {\boldsymbol{Qq}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} = \left( {1+\varepsilon } \right){\boldsymbol{q}}. $

可以看出, $ 1+\varepsilon $正好是数据矩阵 $ {\boldsymbol{Q}} $的特征值, $ {\boldsymbol{q}} $为特征值 $ 1+\varepsilon $对应的特征向量. 星敏感器的精度普遍为角秒级,因此误差因子 $ \varepsilon $非常小,由此,问题转变成求解数据矩阵 $ {\boldsymbol{Q}} $最接近1的特征值 $ {\lambda _{{\text{opt}}}} $$ {\lambda _{{\text{opt}}}} $对应的特征向量即为最优旋转四元数 $ {{\boldsymbol{q}}_{{\text{opt}}}} $. 求解 $ {\lambda _{{\text{opt}}}} $考虑使用牛顿-辛普森梯度迭代法,迭代初始值为1. 注意到,当以初始值为1进行迭代时,得到的 $ {\lambda _{{\text{opt}}}} $有可能是小于1的特征值 $ {\lambda _1} $,也可能是大于1的特征值 $ {\lambda _2} $,且 $ {\lambda _1} $$ {\lambda _2} $满足:

$ \begin{split} & {\boldsymbol{Q}}{{\boldsymbol{q}}}_{\lambda _1}={\lambda }_{1}{{\boldsymbol{q}}}_{\lambda _1}\text{,} {\boldsymbol{Q}}{{\boldsymbol{q}}}_{\lambda _2}={\lambda }_{2}{{\boldsymbol{q}}}_{\lambda _2}. \end{split} $

$ {\lambda _{{\text{opt}}}} $唯一,只能为 $ {\lambda _1} $$ {\lambda _2} $,因此以初始值为1进行迭代不可行. 若能证明数据矩阵 $ {\boldsymbol{Q}} $的特征值均大于1或者均小于1,则上述问题无需考虑.

证明  数据矩阵 $ {\boldsymbol{Q}} $的所有特征值均小于1.

$ {\boldsymbol{Qx}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} = \lambda {\boldsymbol{x}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \Rightarrow {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \left[ {{\boldsymbol{Q}} - \lambda {\boldsymbol{I}}} \right]{\kern 1pt} {\kern 1pt} {\boldsymbol{x}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} =\mathbf{ 0} . $

$ \left. \begin{array}{l} {\boldsymbol{Q}}=\left[{{\boldsymbol{Q}}}_{1}\text{,}{{\boldsymbol{Q}}}_{2}\text{,}{{\boldsymbol{Q}}}_{3}\text{,}{{\boldsymbol{Q}}}_{4}\right]\text{,}\\ {\boldsymbol{x}}= \;\, \left[{x}_{1}\text{,}{x}_{2}\text{,}{x}_{3}\text{,}{x}_{4}\right]{}^{\text{T}}. \end{array} \right\} $

$ \begin{split} & \left[1\text{,} 1\text{,}1\text{,}1\right]\left[{\boldsymbol{Q}}-\lambda I\right]{\boldsymbol{x}}=0 \Rightarrow \\& \qquad \lambda \left({x}_{1}+{x}_{2}+{x}_{3}+{x}_{4}\right)=\\&\qquad \left[1\text{,}1\text{,}1\text{,}1\right]{{\boldsymbol{Q}}}_{1}{x}_{1}+\left[1\text{,}1\text{,}1\text{,}1\right]{{\boldsymbol{Q}}}_{2}{x}_{2}+\\&\qquad \left[1\text{,}1\text{,}1\text{,}1\right]{{\boldsymbol{Q}}}_{3}{x}_{3}+\left[1\text{,}1\text{,}1\text{,}1\right]{{\boldsymbol{Q}}}_{4}{x}_{4}. \end{split} $

若存在特征值 $ {\lambda _i} > 1 $,则

$ \begin{split} & \left[1\text{,}1\text{,}1\text{,}1\right]{{\boldsymbol{Q}}}_{1}+\left[1\text{,}1\text{,}1\text{,}1\right]{{\boldsymbol{Q}}}_{2}+\\&\qquad \left[1\text{,}1\text{,}1\text{,}1\right]{{\boldsymbol{Q}}}_{3}+\left[1\text{,}1\text{,}1\text{,}1\right]{{\boldsymbol{Q}}}_{4} > 4. \end{split} $

计算得到

$ \begin{split} & \left[1 {,}1 {,}1 {,}1\right]{{\boldsymbol{Q}}}_{1}+\left[1 {,}1 {,}1 {,}1\right]{{\boldsymbol{Q}}}_{2}+\\&\qquad \left[1 {,}1 {,}1 {,}1\right]{{\boldsymbol{Q}}}_{3}+\left[1 {,}1 {,}1 {,}1\right]{{\boldsymbol{Q}}}_{4}=\\& {\displaystyle \sum _{k=1}^{n}4{a}_{k}{D}_{ {z},k}^{ {i}}{D}_{ {y},k}^{ {s}}}+{\displaystyle \sum _{k=1}^{n}4{a}_{k}{D}_{ {x},k}^{ {i}}{D}_{ {z},k}^{ {s}}}+{\displaystyle \sum _{k=1}^{n}4{a}_{k}{D}_{ {y},k}^{ {i}}{D}_{ {x},k}^{ {s}}} > 4\\& \iff {\displaystyle \sum _{k=1}^{n}{a}_{k}{D}_{ {z},k}^{ {i}}{D}_{ {y},k}^{ {s}}} + {\displaystyle \sum _{k=1}^{n}{a}_{k}{D}_{ {x},k}^{ {i}}{D}_{ {z},k}^{ {s}}}+{\displaystyle \sum _{k=1}^{n}{a}_{k}{D}_{ {y},k}^{ {i}}{D}_{ {x},k}^{ {s}}} > 1. \end{split} $

$ {a_k}D_{{ {z}},k}^{ {i}}{\kern 1pt} D_{{ {y}},k}^{ {s}}+{a_k}D_{{ {x}},k}^{ {i}}D_{{ {z}},k}^{ {s}}+{a_k}D_{{ {y}},k}^{ {i}}D_{{ {x}},k}^{ {s}} $中的 $ D_{{ {z}},k}^{ {i}}{\kern 1pt} $$ D_{{ {y}},k}^{ {s}} $$ D_{{ {x}},k}^{ {i}} $$ D_{{ {z}},k}^{ {s}} $$ D_{{ {y}},k}^{ {i}} $$ D_{{ {x}},k}^{ {s}} $看成自变量,构造拉格朗日函数 $ f $,令 $ {c_1} $$ {c_2} $为拉格朗日乘子,达到

$ \left. \begin{array}{l} f= {D}_{ {z},k}^{ {i}}{D}_{ {y},k}^{ {s}}+{D}_{ {x},k}^{ {i}}{D}_{ {z},k}^{ {s}}+{D}_{ {y},k}^{ {i}}{D}_{ {x},k}^{ {s}}+\\ \qquad {c}_{1}\left(1-{\left({D}_{ {x},k}^{ {i}}\right)}^{2}-{\left({D}_{ {y},k}^{ {i}}\right)}^{2}-{\left({D}_{ {z},k}^{ {i}}\right)}^{2}\right)+\\ \qquad {c}_{2}\left(1-{\left({D}_{ {x},k}^{ {s}}\right)}^{2}-{\left({D}_{ {y},k}^{ {s}}\right)}^{2}-{\left({D}_{ {z},k}^{ {s}}\right)}^{2}\right) {;}\\ {{\rm{s}}} {.{\rm{t}}}. \;\;\; {\left({D}_{ {x},k}^{ {i}}\right)}^{2}+{\left({D}_{ {y},k}^{ {i}}\right)}^{2}+{\left({D}_{ {z},k}^{ {i}}\right)}^{2}=1 {,}\\ \;\;\;\;\;\;\;\;{\left({D}_{ {x},k}^{ {s}}\right)}^{2}+{\left({D}_{ {y},k}^{ {s}}\right)}^{2}+{\left({D}_{ {z},k}^{ {s}}\right)}^{2}=1. \end{array} \right\} $

$ f' = 0 $,得到

$ {D}_{ {x},k}^{ {s}}={D}_{ {y},k}^{ {i}} {,}{D}_{ {y},k}^{ {s}}={D}_{ {z},k}^{ {i}} {,}{D}_{ {z},k}^{ {s}}={D}_{ {x},k}^{ {i}}. $

$ \begin{split} & {\left( {{a_k}D_{{ {z}},k}^{ {i}}D_{{ {y}},k}^{ {s}}+{a_k}D_{{ {x}},k}^{ {i}}D_{{ {z}},k}^{ {s}}+{a_k}D_{{ {y}},k}^{ {i}}D_{{ {x}},k}^{ {s}}} \right)_{\max }} = \\& \qquad {a_k}\left( {{{\left( {D_{{ {x}},k}^{ {i}}} \right)}^2}+{{\left( {D_{{ {y}},k}^{ {i}}} \right)}^2}+{{\left( {D_{{ {z}},k}^{ {i}}} \right)}^2}} \right) = {a_k}. \end{split} $

得到

$ \begin{split} & {\left( {\sum\limits_{k = 1}^n {{a_k}D_{{ {z}},k}^{ {i}}{\kern 1pt} D_{{ {y}},k}^{ {s}}} +\sum\limits_{k = 1}^n {{a_k}D_{{ {x}},k}^{ {i}}D_{{ {z}},k}^{ {s}}} +\sum\limits_{k = 1}^n {{a_k}D_{{ {y}},k}^{ {i}}D_{{ {x}},k}^{ {s}}} } \right)_{\max }} = \\& \qquad \quad {a_1}+{a_2}+{a_3}+ \cdots {a_n} = 1. \\[-5pt] \end{split} $

式(39)与式(35)矛盾,因此 $ {\boldsymbol{Q}} $的所有特征值都必定小于1,得证.

证明结果表明,可以基于牛顿-辛普森梯度迭代法,以迭代初始值为1求解 $ {\boldsymbol{Q}} $的最接近1的特征值. 实际上, $ {\boldsymbol{Q}} $最接近1的特征值即最大特征值. 为了求解 $ {\boldsymbol{Q}} $最接近1的特征值,须得到关于特征值 $ \lambda $的函数. $ {\boldsymbol{Q}} $的特征多项式为

$ f\left( \lambda \right) = \det \left( {{\boldsymbol{Q}} - \lambda {\boldsymbol{I}}} \right) = {\lambda ^4}+{r_1}{\lambda ^2}+{r_2}\lambda +{r_3}. $

$ \left. \begin{array}{l} \begin{array}{l}{r}_{1} = - 2\left( {E}_{11}^{2} + {E}_{12}^{2} + {E}_{13}^{2} + {E}_{21}^{2} + {E}_{22}^{2} + {E}_{23}^{2} + {E}_{31}^{2} + {E}_{32}^{2} + {E}_{33}^{2} \right) \text{,}\\ {r}_{2} = 8( {E}_{13}{E}_{22}{E}_{31} - {E}_{12}{E}_{23}{E}_{31} - {E}_{13}{E}_{21}{E}_{32} + {E}_{11}{E}_{23}{E}_{32} + \\ \qquad {E}_{12}{E}_{21}{E}_{33} - {E}_{11}{E}_{22}{E}_{33} )\text{,}\\ {r}_{3}=\mathrm{det} \;\; {\boldsymbol{Q}} .\end{array} \end{array} \right\} $

在星敏感器输出星光矢量的天球系坐标和星敏系坐标后,可以并行计算得到 $ {r_1} $$ {r_2} $$ {r_3} $. 特征多项式 $ f\left( \lambda \right) $的导数为

$ f'\left( \lambda \right) = 4{\lambda ^3}+2{r_1}\lambda +{r_2}. $

经过6次迭代,使所得特征值达到符合要求的精度. 求出特征值 $ {\lambda _{{\text{opt}}}} $后,计算 $ {\lambda _{{\text{opt}}}} $对应的特征向量 $ {{\boldsymbol{q}}}_{{\rm{opt}}}={\left[{q}_{0},{q}_{1},{q}_{2},{q}_{3}\right]}^{\text{T}} $,方法如下.

$ \begin{split} &{\left[\boldsymbol{Q}-\lambda_{\mathrm{opt}} \boldsymbol{I}\right] \boldsymbol{q}_{\mathrm{opt}}^{\prime}=\mathbf{0 } \Rightarrow} \\ &\boldsymbol{Q}-\lambda_{\mathrm{opt}} \boldsymbol{I}=\left[\begin{array}{llll} a_{11}&a_{12}&a_{13}&a_{14} \\ a_{21}&a_{22}&a_{23}&a_{24} \\ a_{31}&a_{32}&a_{33}&a_{34} \\ a_{41}&a_{42}&a_{43}&a_{44} \end{array}\right]= \\ &{\left[ {\begin{array}{*{20}{c}} Q_{11}-\lambda&Q_{12}&Q_{13}&Q_{14} \\ Q_{21}&Q_{22}-\lambda&Q_{23}&Q_{24} \\ Q_{31}&Q_{32}&Q_{33}-\lambda&Q_{34} \\ Q_{41}&Q_{42}&Q_{43}&Q_{44}-\lambda \end{array}} \right]} . \end{split} $

实施初等行变换,将 $ \left[ {{\boldsymbol{Q}} - {\lambda _{{\text{opt}}}}{\boldsymbol{I}}} \right]{{\boldsymbol{q'}}_{{\text{opt}}}}= \mathbf{0 } $变换为

$ \left. \begin{array}{l} \begin{split} & {a}_{1,1}{{q}^{\prime }}_{0}+{a}_{1,2}{{q}^{\prime }}_{1}+{a}_{1,3}{{q}^{\prime }}_{2}+{a}_{1,4}{{q}^{\prime }}_{3}=0\text{,}\\& {b}_{2,2}{{q}^{\prime }}_{1}+{b}_{2,3}{{q}^{\prime }}_{2}+{b}_{2,4}{{q}^{\prime }}_{3}=0\text{,}\\& {c}_{3,3}{{q}^{\prime }}_{2}+{c}_{3,4}{{q}^{\prime }}_{3}=0\text{,} {d}_{4,4}{{q}^{\prime }}_{3}=0\text{.} \end{split} \end{array} \right\} $

只有 $ {d_{4,4}} $=0,才能使得矩阵 $ {\boldsymbol{Q}} - {\lambda _{{\text{opt}}}}{\boldsymbol{I}} $非满秩,特征向量 ${{\boldsymbol{q}}'_{{\text{opt}}}}$非零. 令 $ {q'_3} = - 1 $,求解 ${{{\boldsymbol{q}}}}'_{{\rm{opt}}}= {[{{q} }'_{0},{{q} }'_{1},{{q} }'_{2},{{q} }'_{3}]}^{\text{T}}$的解析表达式为

$ \left. \begin{array}{l}{{q}^{\prime }}_{0}=\dfrac{{a}_{14}}{{a}_{11}}-\dfrac{{a}_{12}{b}_{24}}{{a}_{11}{b}_{22}}+\dfrac{{a}_{12}{b}_{23}{c}_{34}}{{a}_{11}{b}_{22}{c}_{33}}-\dfrac{{a}_{13}{c}_{34}}{{a}_{11}{c}_{33}}\text{,}\\ {{q}^{\prime }}_{1}=\dfrac{{b}_{24}}{{b}_{22}}-\dfrac{{b}_{23}{c}_{34}}{{b}_{22}{c}_{33}}\text{,}{{q}^{\prime }}_{2}=\dfrac{{c}_{34}}{{c}_{33}}\text{,}{{q}^{\prime }}_{3}=-1. \end{array} \right\} $

其中

$ \left. \begin{array}{l} {b}_{22}={a}_{22}-\dfrac{{a}_{21}{a}_{12}}{{a}_{11}}\text{,}{b}_{33}={a}_{33}-\dfrac{{a}_{31}{a}_{13}}{{a}_{11}}\text{,}\\ {b}_{23}={a}_{23}-\dfrac{{a}_{21}{a}_{13}}{{a}_{11}}\text{,}{b}_{34}={a}_{34}-\dfrac{{a}_{31}{a}_{14}}{{a}_{11}}\text{,}\\ {b}_{24}={a}_{24}-\dfrac{{a}_{21}{a}_{14}}{{a}_{11}}\text{,}{c}_{33}={b}_{33}-\dfrac{{b}_{32}{b}_{23}}{{b}_{22}}\text{,}\\ {b}_{32}={a}_{32}-\dfrac{{a}_{31}{a}_{12}}{{a}_{11}}\text{,}{c}_{34}={b}_{34}-\dfrac{{b}_{32}{b}_{24}}{{b}_{22}}. \end{array} \right\} $

由于特征向量 $ {{\boldsymbol{q}}'_{{\rm{opt}}}} $的模一定不为1,在对其进行归一化处理后得到最优旋转四元数 $ {{\boldsymbol{q}}_{{\rm{opt}}}} $

$ \left. \begin{array}{l} \left|{{{\boldsymbol{q}}}}'_{{\rm{opt}}}\right|=\sqrt{{{q}^{\prime }}_{0}^{2}+{{q}^{\prime }}_{1}^{2}+{{q}^{\prime }}_{2}^{2}+{{q}^{\prime }}_{3}^{2}}\text{,}\\ {{\boldsymbol{q}}}_{{\rm{opt}}}=\dfrac{{{{\boldsymbol{q}}}}'_{{\rm{opt}}}}{\left|{{{\boldsymbol{q}}} }'_{{\rm{opt}}}\right|}. \end{array} \right\} $

在式(45)中,若令 $ {q'_3} = 1 $,则可以得到 $ {\boldsymbol{q}}_{{\rm{opt}}}^{''}= {\left[{{q}}'_{0},{{q}}'_{1},{{q}}'_{2},{{q}}'_{3}\right]}^{\text{T}} $的另一种解析表达式为

$ \left. \begin{array}{l} {{q}}'_{0}=-\left(\dfrac{{a}_{14}}{{a}_{11}}-\dfrac{{a}_{12}{b}_{24}}{{a}_{11}{b}_{22}}+\dfrac{{a}_{12}{b}_{23}{c}_{34}}{{a}_{11}{b}_{22}{c}_{33}}-\dfrac{{a}_{13}{c}_{34}}{{a}_{11}{c}_{33}}\right)\text{,}\\ {{q} }'_{1}=-\left(\dfrac{{b}_{24}}{{b}_{22}}-\dfrac{{b}_{23}{c}_{34}}{{b}_{22}{c}_{33}}\right)\text{,}{{q}^{\prime }}_{2}=-\left(\dfrac{{c}_{34}}{{c}_{33}}\right)\text{,}{{q}^{\prime }}_{3}=1. \end{array} \right\} $

$ {\boldsymbol{q}}_{{\rm{opt}}}^{'}= - {\boldsymbol{q}}_{{\rm{opt}}}^{''}$,进而可以得到

$ {{\boldsymbol{q}}_{{\text{opt}}}} = \frac{ { {\boldsymbol{q}}_{{\rm{opt}}}^{''} }}{{\left| { {\boldsymbol{q}}_{{\rm{opt}}}^{''} } \right|}}. $

求出的2个最优旋转四元数互为负向量. 用旋转四元数的轴角表示角度,得到

$ \left. \begin{array}{l} \begin{split} & {\boldsymbol{q}}=\mathrm{cos}\dfrac{\phi }{2}+{\boldsymbol{u}}\mathrm{sin}\dfrac{\phi }{2}\Rightarrow \left[{\boldsymbol{u}},\phi \right]\text{,}\\& -{\boldsymbol{q}}=-\mathrm{cos}\dfrac{\phi }{2}-{\boldsymbol{u}}\mathrm{sin}\dfrac{\phi }{2}=\\& \mathrm{cos}\dfrac{2\text{π}-\phi }{2}+\left(-{\boldsymbol{u}}\right)\mathrm{sin}\dfrac{2\text{π}-\phi }{2}\Rightarrow \left[-{\boldsymbol{u}},2\text{π}-\phi \right]. \end{split} \end{array} \right\} $

即旋转四元数 $ {\boldsymbol{q}} $的含义是绕 $ {\boldsymbol{u}} $旋转 $ \phi $角度;旋转四元数 $ - {\boldsymbol{q}} $的含义是绕 $ - {\boldsymbol{u}} $旋转 $ 2{\text{π }} - \phi $角度. 可以看出,二者含义相同, $ {\boldsymbol{q}} $$ - {\boldsymbol{q}} $代表同样的旋转变化. 从四元数的方向余弦阵形式考虑,即将 $ {\boldsymbol{q}} $$ - {\boldsymbol{q}} $分别代入式(13),容易得出 $ {\boldsymbol{q}} $$ - {\boldsymbol{q}} $对应的方向余弦阵相等. 式(50)从方向余弦阵的角度推导出 $ {\boldsymbol{q}} $$ - {\boldsymbol{q}} $都能够得到正确的姿态. 可知, $ {q'_3} $的取值不影响星敏感器三轴旋转角的取值.

3. 仿真验证及分析

利用数值仿真对比验证三轴旋转角求解算法和传统的三轴旋转角求解算法的姿态解算精度. 仿真步骤如下.

1)随机生成单位正交矩阵. 2)随机生成三维矢量 ${{{{\boldsymbol{D}}}}}_{k}^{i}={\left[{{{D}}}_{x,k}^{i}, {{{D}}}_{y,k}^{i}, \;{{{D}}}_{z,k}^{i}\right]}^{\text{T}},(k=1,2,\cdots ,m)$,其中 $ {\boldsymbol{D}}_k^i $为第 $ k $个星光矢量在天球坐标系中的三维坐标, $ m = 15 $. 星敏感器存在系统噪声与随机噪声影响,同时星敏感器测量星光矢量在星敏坐标系中的三维坐标时存在噪声误差. 第 $ k $个星光矢量在星敏坐标系中的噪声误差向量 $ {{\boldsymbol{\varepsilon}} }_{k}={\left[{\varepsilon }_{x,k},{\varepsilon }_{y,k},{\varepsilon }_{z,k}\right]}^{\text{T}} $,令 $ D_{j,k}^i $服从均值为0、标准差为20的高斯分布; $ {\varepsilon _{j,k}} $服从均值为0,标准差为0.001的高斯分布. 3)计算 $ {{\boldsymbol{D}}}_{k}^{s}={\left[{D}_{x,k}^{s}, {D}_{y,k}^{s},{D}_{z,k}^{s}\right]}^{\text{T}} $${\boldsymbol{D}}_k^s = {\boldsymbol{C}}{\boldsymbol{D}}_k^i+{{\boldsymbol{\varepsilon }}_k},$ $ (k = 1, 2, \cdots , m) $. 4)分别采用SVD算法、LS算法、QUEST算法和本研究算法求解三轴旋转角. 其中SVD算法、LS算法先求解最优旋转矩阵,再得到三轴旋转角;QUEST算法、本研究算法先求解最优旋转四元数,再得到三轴旋转角. 5)计算4种算法的三轴旋转角误差. 6)计算4种算法的耗时. 7)重复步骤1)~6),共计 $ n = 1 \; 000 $次. 仿真使用Matlab 2018b数值计算软件在计算机上进行,系统配置为Inter i7-8550U 2.0GHz CPU,32G RAM[22-26].

图2所示为噪声标准差 ${\sigma _{\rm{s}}} $=1.0×10−3的情况下,不同算法的三轴姿态角误差曲线. 图中,n为样本数,φ1为横滚角误差,φ2为俯仰角误差,φ3为方位角误差. 可以看出,SVD算法、QUEST算法和本研究算法的定姿误差分散程度较小,三轴姿态角误差处于−0.03°~0.03°,该误差主要是由较强的噪声引起,数值计算截断误差可以忽略. 相比之下,LS算法的定姿误差分散程度较大,三轴姿态角误差处于−0.06°~0.06°,该误差主要由较强的噪声和模型的精度引起. 如表1所示为星敏系矢量坐标受到不同噪声标准差影响时,本研究算法与3种传统算法的姿态角误差标准差. 表中,σA为姿态角误差标准差. 仿真数据计算可知,在σs=1.0×10−3的情况下,SVD算法、QUEST算法和本研究算法的姿态角误差标准差为0.027 51°;LS算法的姿态角误差标准差为0.037 94°.

图 2

图 2   1.0×10−3噪声标准差下不同算法的三轴姿态角误差曲线

Fig.2   Three-axis attitude angle error curve of different algorithms under 1.0×10−3 standard deviation of noise


表 1   不同噪声标准差下不同算法的姿态角误差标准差

Tab.1  Standard deviation of attitude angle error for different algorithms with different noise standard deviations

算法 σA/(°)
σs=1.0×10−2 σs=1.0×10−3 σs=1.0×10−4 σs=1.0×10−5
SVD 2.071×10−1 2.751×10−2 3.084×10−3 4.055×10−4
LS 3.679×10−1 3.794×10−2 4.286×10−3 4.622×10−4
QUEST 2.071×10−1 2.751×10−2 3.084×10−3 4.055×10−4
本研究 2.071×10−1 2.751×10−2 3.084×10−3 4.055×10−4

新窗口打开| 下载CSV


不同的噪声强度对上述4种算法的处理速度影响不大,因此以σs= $ 1.0 \times {\text{1}}{{\text{0}}^{{{ - 3}}}} $为例,分析不同算法的处理速度,如图3所示. 可以看出,本研究算法的耗时总体趋势在t= $ 4.0 \times {\text{1}}{{\text{0}}^{{{ - 5}}}} $ s上下波动,计算速度最快;QUEST、SVD定姿算法的耗时总体趋势的均值在 $ 7.0 \times {\text{1}}{{\text{0}}^{{{ - 5}}}} $s上下波动;LS算法的耗时总体趋势在 $ 1.0 \times {\text{1}}{{\text{0}}^{{{ - 4}}}} $ s上下波动,计算速度最慢. 如表2所示为星敏系矢量坐标受到不同噪声标准差影响时,本研究算法与3种传统算法的平均耗时. 在σs= $ 1.0 \times {\text{1}}{{\text{0}}^{{{ - 3}}}} $的情况下,本研究算法、QUEST算法、SVD算法、LS算法的平均耗时分别为 $ 3.322 \times {\text{1}}{{\text{0}}^{{{ - 5}}}} $$ 5.671 \times {\text{1}}{{\text{0}}^{{{ - 5}}}} $$ 5.971 \times {\text{1}}{{\text{0}}^{{{ - 5}}}} $$ 9.219 \times {\text{1}}{{\text{0}}^{{{ - 5}}}} $ s. 可以看出,σs= $ 1.0 \times {\text{1}}{{\text{0}}^{{{ - 3}}}} $的情况下,本研究算法的计算速度最快,LS算法的计算速度最慢、计算效率最低.

图 3

图 3   1.0×10−3噪声标准差下不同算法的时间消耗曲线

Fig.3   Time consumption curve of different algorithms under 1.0×10−3 standard deviation of noise


表 2   不同噪声标准差下不同算法的耗时表

Tab.2  Time-consuming table for different algorithms with different noise standard deviations

算法 t/s
σs=1×10−2 σs=1×10−3 σs=1×10−4 σs=1×10−5
SVD 6.021×10−5 5.971×10−5 6.153×10−5 6.568×10−5
LS 9.354×10−5 9.219×10−5 9.612×10−5 9.674×10−5
QUEST 5.721×10−5 5.671×10−5 5.953×10−5 5.368×10−5
本研究 3.284×10−5 3.322×10−5 3.268×10−5 3.109×10−5

新窗口打开| 下载CSV


星敏感器是嵌入式设备,其运行环境也是嵌入式环境. 为了检验本研究算法在工程实践的应用价值,在嵌入式平台上对本研究算法与对比算法进行测试. 嵌入式设备型号为DSP6748,星敏感器的视场为 ${18^ \circ } \times {18^ \circ }$,焦距为42 mm,分辨率为 $2 \;048 \times 2 \;048$,像元尺寸为11 μm,探测灵敏度为 $6.5$ ${\text{mV}}$. 3倍标准差下星敏感器的精度为 $3''$. 取多次测试结果的平均值作为算法的最终耗时,SVD算法、LS算法、QUEST算法和本研究算法在嵌入式环境下的耗时测试结果分别为2.274×10−3、3.519×10−3、2.012×10−3、1.458×10−3 s. 可以看出,算法在嵌入式设备上的测试结果与在计算机上的仿真结果结论一致,本研究算法在计算效率上有一定的提高.

仿真与试验验证可知,在旋转角精度方面,本研究算法和SVD算法、QUEST算法的精度相当,LS算法稍差. 实际上,前3种算法满足Wahba问题的极限精度,由于算法运用的数学手段不同,导致三者在算法形式上有所差异,但结果一致,都是对三轴旋转角的最优估计. 在耗时方面,本研究算法的计算效率最高,QUEST算法和SVD算法的计算效率居中,LS算法的计算效率最低,结果体现出本研究算法在算法性能上的进步. 试验结果表明,与3种传统算法相比,保持相同最高估计精度的情况下,本研究算法在计算效率上有一定的提高.

4. 结 语

本研究算法基于单个星光矢量在天球坐标系与星敏坐标系中的三维坐标推导得到四元数变换形式,充分利用伪逆矩阵的性质对二次四元数矩阵方程进行降次,并将最优旋转四元数的求解转变为求解矩阵的最大特征值及其对应的特征向量问题. 本研究从理论上证明了特征值的取值范围,解决了三轴旋转角求解过程中迭代问题的前提条件;详细分析了四元数矢量方向问题,排除了方向奇异性带来的三轴旋转角求解问题. 对本研究算法与传统算法的性能仿真与试验对比分析结果表明,本研究算法和SVD算法、QUEST算法的三轴旋转角求解精度最高,本研究算法的计算效率较高、综合性能较好.

参考文献

CHANG L B, QIN F J, LI A

A novel backtracking scheme for attitude determination-based initial alignment

[J]. IEEE Transactions on Automation Science and Engineering, 2015, 12 (1): 384- 390

DOI:10.1109/TASE.2014.2346581      [本文引用: 1]

DU S, GAO Y

Inertial aided cycle slip detection and identification for integrated PPP GPS and INS

[J]. Sensors, 2012, 12 (11): 14344- 14362

DOI:10.3390/s121114344      [本文引用: 1]

MARKLEY F L, MORTARI D

How to estimate attitude from vector observations

[J]. Advance Astronomy Science, 2000, 103: 1979- 1996

[本文引用: 1]

YANG Y, ZHOU Z

An analytic solution to Wahba’s problem

[J]. Aerospace Science and Technology, 2013, 30 (1): 46- 49

DOI:10.1016/j.ast.2013.07.002      [本文引用: 1]

SHUSTER M D, OH S D

Three-axis attitude determination from vector observations

[J]. Journal of Guidance Control, 1981, 4 (1): 70- 77

DOI:10.2514/3.19717      [本文引用: 1]

MORTARI D

ESOQ: a closed-form solution to the Wahba problem

[J]. The Journal of the Astronomy Sciences, 1997, 45: 195- 204

[本文引用: 1]

HORN R A, JOHNSON C R. Matrix analysis [M]. Cambridge: Cambridge University Press, 1985: 57-73.

[本文引用: 2]

FARRELL J L, STUELPNAGEL J C, WESSNER R H, et al

A least squares estimate of spacecraft attitude

[J]. SIAM Review, 1966, 8 (3): 384- 386

DOI:10.1137/1008080      [本文引用: 1]

ZHANG Y J, ZHENG M T, XIONG J X, et al

On-orbit geometric calibration of ZY-3 three-linear array imagery with multistrip data sets

[J]. IEEE Transactions on Geoscience and Remote Sensing, 2014, 52 (1): 224- 234

DOI:10.1109/TGRS.2013.2237781      [本文引用: 1]

RATNAWEERA A, HALGAMUGE S K, WATSON H C

Self-organizing hierarchical particle swarm optimizer with time-varying acceleration coefficients

[J]. IEEE Transactions on Evolutionary Computation, 2004, 8 (3): 240- 255

DOI:10.1109/TEVC.2004.826071     

WU J, ZHOU Z B, GAO B, et al

Fast linear quaternion attitude estimator using vector observations

[J]. IEEE Transactions on Automation Science and Engineering, 2018, 15 (1): 307- 319

DOI:10.1109/TASE.2017.2699221     

LUO L, HUANG Y L, ZHANG Z, et al

A position loci-based in-motion initial alignment method for low-cost attitude and heading reference system

[J]. IEEE Transactions on Instrumentation and Measurement, 2021, 70: 7500618

ZHOU X R, XU X, YAO Y Q, et al

A robust quaternion Kalman filter method for MIMU/GPS in-motion alignment

[J]. IEEE Transactions on Instrumentation and Measurement, 2021, 70: 8503109

[本文引用: 1]

FRASER C S, HANLEY H B

Bias compensation in rational functions for Ikonos satellite imagery

[J]. Photogrammetric Engineering and Remote Sensing, 2003, 69 (1): 53- 57

DOI:10.14358/PERS.69.1.53      [本文引用: 1]

WU Y H, GAO Y, LIN J W, et al

Low-cost, high-performance monocular vision system for air bearing table attitude determination

[J]. Journal of Spacecraft and Rockets, 2014, 51 (1): 66- 75

DOI:10.2514/1.A32465     

WU J

Real-time magnetometer disturbance estimation via online nonlinear programming

[J]. IEEE Sensors Journal, 2019, 19 (12): 4405- 4411

DOI:10.1109/JSEN.2019.2901925     

WU J, ZHOU Z B, FOURATE H, et al

Generalized linear quaternion complementary filter for attitude estimation from multisensor observations: an optimization approach

[J]. IEEE Transactions on Automation Science and Engineering, 2019, 16 (3): 1330- 1343

DOI:10.1109/TASE.2018.2888908      [本文引用: 1]

ZHANG G, JIANG Y H, LI D, et al

In-orbit geometric calibration and validation of ZY-3 linear array sensors

[J]. The Photogrammetric Record, 2014, 29 (145): 68- 88

DOI:10.1111/phor.12052      [本文引用: 1]

LEFFERTS E J, MARKLEY F L, SHUSTER M D

Kalman filtering for spacecraft attitude estimation

[J]. Journal of Guidance, 1982, 5 (5): 417- 429

DOI:10.2514/3.56190     

LI J C, GAO W, ZHANG Y, et al

Gradient descent optimization-based self-alignment method for stationary SINS

[J]. IEEE Transactions on Instrumentation and Measurement, 2019, 68 (9): 3278- 3286

DOI:10.1109/TIM.2018.2878071     

MORTARI D, ROMOLI A. StarNav III: a three field of view star tracker [C]// IEEE Aerospace Conference Proceeding. Big Sky: IEEE, 2002: 47-67.

[本文引用: 1]

DOUIK A, LIU X, BALLAL T, et al

Precise 3-D GNSS attitude determination based on Riemannian manifold optimization algorithms

[J]. IEEE Transactions on Signal Processing, 2020, 68: 284- 299

DOI:10.1109/TSP.2019.2959226      [本文引用: 1]

COLE C L, CRASSIDIS J L

Fast star-pattern recognition using planar triangles

[J]. Journal of Guidance, Control, and Dynamics, 2006, 29 (1): 64- 71

LIEBE C C

Star trackers for attitude determination

[J]. IEEE Aerospace and Electronic Systems Magazine, 1995, 10 (6): 10- 16

DOI:10.1109/62.387971     

KUDVA P, THROCKMORTON A

Attitude determination studies for the earth observation system AM1 (EOS-AM1) mission

[J]. Journal of Guidance, Control, and Dynamic, 1996, 19 (6): 1326- 1331

DOI:10.2514/3.21789     

CANDAN B, SOKEN H E

Robust attitude estimation using IMU-only measurements

[J]. IEEE Transactions on Instrumentation and Measurement, 2021, 70: 9512309

[本文引用: 1]

/