星光/惯性组合导航系统由于其精度高、自主性强等优点，在军用领域有着广泛的应用需求^[1-3]。在战时状态，这一类组合导航系统与卫星/惯性组合导航系统相比，在安全性方面更具优势^[4-5]。

在星光/惯性组合导航中，惯导与星跟踪器之间一般采取固联安装，无需考虑隔离平台运动。惯导向星跟踪器传递导航信息包括协调世界时、惯导姿态以及地理经纬度等信息^[2]，星跟踪器借助辅助信息完成星图曝光、星点识别、矢量计算等，之后组合滤波解算出载体在导航系的姿态和位置等。实际上，惯导存在姿态误差，惯导和星跟踪器之间的安装阵由于力热变化也存在漂移，惯导的姿态误差和安装阵在组合导航算法中属于待估计的参数，对这些参数的估计是组合导航算法中必须要面临的问题。在全天时星光惯性组合导航中，由于大气层内的高背景杂光的影响，为了提高信噪比满足全天时的探测需求，一般选择小视场长焦折反式光学系统^[6]，如美国诺斯罗普公司用于RC-135侦察机、B-2轰炸机的典型产品LN-120G视场只有6″。为了提高白昼观星的探测概率，组合导航系统一般自带摆镜，采用主动寻星的方式工作，每次只跟踪单颗恒星。本文的研究工作就是基于这一类组合导航方式。

Wang^[7]、王可东^[8]和付建楠^[9]等对星光惯性组合导航的系统误差进行了建模，但没有考虑安装矩阵的影响。张鹏^[10]、杨波^[11]、张磊^[12]和林星辰^[13]将安装矩阵作为系统的状态观测量，基于卡尔曼滤波方法进行统一求解，但是由于惯导漂移随时间累积的特点，滤波算法收敛较慢，无法满足星跟踪器恒星识别的实时要求。李新鹏等^[14]提出了一种基于四元数自适应卡尔曼滤波的星敏感器安装矩阵在轨实时校准方法，但是该方法依赖于多星敏感器之间安装矩阵的在轨互标，对惯导的姿态精度有一定要求，工程实用性不足。焦宏伟等^[15]基于导航星矢量与地心矢量夹角观测值与真值相等的条件建立了状态模型和测量模型，对安装矩阵的求解假定了组合导航系统的位置不变，因而只适合估计惯导的初始对准偏差。Ning等^[16]通过分析星敏感器的观测数据得到存在紧耦合的安装误差，本文也基于观测数据得到安装阵的迭代计算初值。张广军等^[17]给出了一种基于Kalman滤波的在轨安装阵标定方法，与本文累积多颗星的方法有很大不同。Kim等^[18]考虑了时间同步对组合导航算法的影响，本文算法未涉及这一方面，主要是通过算法之前的硬同步实现。熊智等^[19]提出了地理系下耦合位置误差的机载惯性/星光组合导航滤波方法，该方法基于地理系下姿态量测方程线性化，其算法应用场合与本文方法有很多不同。Jam等^[20]提出了一种针对CCD(Charge Coupled Device)星敏感器的组合导航算法，该算法不具有全天时的应用能力。

本文针对全天时星跟踪器相对捷联惯导的安装矩阵在线快速估计方法开展了研究，主要用于消除安装阵的漂移来提升组合导航系统的精度。由于星跟踪器视场小，每次寻星的视场内只出现一颗恒星，而要保证能估计出多个未知参数，需要将多颗恒星的信息转移至同一时刻同一坐标系进行解算。分析了惯性天文组合导航的系统误差、安装阵与星跟踪器测量信息之间的关系，通过采用Levenberg-Marquat (L-M)算法来保证实时求解。实际的跑车试验结果证明了这种方法具有良好的工程实用性。

1 星光惯性组合导航安装阵误差模型

对于组合导航系统，常用的坐标系定义有：① J2000.0惯性参考系C_i；②地心地固坐标系C_e；③地理坐标系C_n，采用东北天坐标系；④惯导本体系C_g；⑤星跟踪器测量系C_s。

在星光惯性组合导航中，定义t=n时刻的惯导的姿态真值和测量值分别为C_gn、Ĉ_gn，t=n时刻相对于t=n-1时刻对应的姿态真值变化和测量值变化分别为C_Rn、Ĉ_Rn，定义惯导的初始姿态误差${\mathit{\boldsymbol{ \boldsymbol{\varDelta} }}_{{\rm{e0}}}}$，有

$ {\mathit{\boldsymbol{\hat C}}_{{\rm{g0}}}} = {\mathit{\boldsymbol{ \boldsymbol{\varDelta} }}_{{\rm{e0}}}}{\mathit{\boldsymbol{C}}_{{\rm{g0}}}} $	（1）

$ \begin{array}{*{20}{l}} {{\mathit{\boldsymbol{C}}_{{\rm{g}}n}} = {\mathit{\boldsymbol{C}}_{{\rm{R}}n}}{\mathit{\boldsymbol{C}}_{{\rm{g}}(n - 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} {\kern 1pt} {\mathit{\boldsymbol{C}}_{{\rm{R}}n}}{\mathit{\boldsymbol{C}}_{{\rm{R}}(n - 1)}}{\mathit{\boldsymbol{C}}_{{\rm{R}}(n - 2)}} \cdots {\mathit{\boldsymbol{C}}_{{\rm{R1}}}}{\mathit{\boldsymbol{C}}_{{\rm{g0}}}}} \end{array} $

（2）

$ \begin{array}{l} {{\mathit{\boldsymbol{\hat C}}}_{{\rm{g}}n}} = {{\mathit{\boldsymbol{\hat C}}}_{{\rm{R}}n}}{{\mathit{\boldsymbol{\hat C}}}_{{\rm{g}}(n - 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} {\kern 1pt} {\kern 1pt} {{\mathit{\boldsymbol{\hat C}}}_{{\rm{R}}n}}{{\mathit{\boldsymbol{\hat C}}}_{{\rm{R}}(n - 1)}}{{\mathit{\boldsymbol{\hat C}}}_{{\rm{R}}(n - 2)}} \cdots {{\mathit{\boldsymbol{\hat C}}}_{{\rm{R1}}}}{{\mathit{\boldsymbol{\hat C}}}_{{\rm{g0}}}} \end{array} $

（3）

定义${\mathit{\boldsymbol{ \boldsymbol{\varDelta} }}_{{\rm{R}}n}}$为测量值Ĉ_Rn到其真实值C_Rn的误差阵，有

$ {\mathit{\boldsymbol{\hat C}}_{{\rm{R}}n}} = {\mathit{\boldsymbol{ \boldsymbol{\varDelta} }}_{{\rm{R}}n}}{\mathit{\boldsymbol{C}}_{{\rm{R}}n}} $

考虑惯导短时精度高且采样周期短(本项目所用惯导采样周期5 ms)，有${\mathit{\boldsymbol{ \boldsymbol{\varDelta} }}_{{\rm{R}}n}}$，则

$ \begin{array}{*{20}{c}} {{{\mathit{\boldsymbol{\hat C}}}_{{\rm{g}}n}}\mathit{\boldsymbol{\hat C}}_{{\rm{g}}(n - 1)}^{\rm{T}} = {{\mathit{\boldsymbol{\hat C}}}_{{\rm{R}}n}} = {\mathit{\boldsymbol{ \boldsymbol{\varDelta} }}_{{\rm{R}}n}}{\mathit{\boldsymbol{C}}_{{\rm{R}}n}} \buildrel\textstyle.\over= }\\ {{\mathit{\boldsymbol{C}}_{{\rm{R}}n}} = {\mathit{\boldsymbol{C}}_{{\rm{g}}n}}\mathit{\boldsymbol{C}}_{{\rm{g}}(n - 1)}^{\rm{T}}} \end{array} $

（4）

式(4)说明惯导t=n时刻相对t=n-1时刻姿态真值变化量可以用对应的测量值变化量代替。

根据式(4)且考虑到本项目所用10颗星的周期不超过30 s，以试验中时漂为0.01 (°)/h的惯导为例，30 s内等效的时漂误差为0.3″，可以忽略不计。因此，有

$ \begin{array}{l} {{\mathit{\boldsymbol{\hat C}}}_{{\rm{g}}n}}\mathit{\boldsymbol{\hat C}}_{{\rm{g1}}}^{\rm{T}} = {{\mathit{\boldsymbol{\hat C}}}_{{\rm{g}}n}}\mathit{\boldsymbol{\hat C}}_{{\rm{g}}(n - 1)}^{\rm{T}}{{\mathit{\boldsymbol{\hat C}}}_{{\rm{g}}(n - 1)}}\mathit{\boldsymbol{\hat C}}_{{\rm{g}}(n - 2)}^{\rm{T}} \cdots {{\mathit{\boldsymbol{\hat C}}}_{{\rm{g2}}}}\mathit{\boldsymbol{\hat C}}_{{\rm{g1}}}^{\rm{T}} \buildrel\textstyle.\over= \\ {\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} {\mathit{\boldsymbol{C}}_{{\rm{g}}n}}\mathit{\boldsymbol{C}}_{{\rm{g}}(n - 1)}^{\rm{T}}{\mathit{\boldsymbol{C}}_{{\rm{g}}(n - 1)}}\mathit{\boldsymbol{C}}_{{\rm{g}}(n - 2)}^{\rm{T}} \cdots {\mathit{\boldsymbol{C}}_{{\rm{g2}}}}\mathit{\boldsymbol{C}}_{{\rm{gl}}}^{\rm{T}} = {\mathit{\boldsymbol{C}}_{{\rm{g}}n}}\mathit{\boldsymbol{C}}_{{\rm{g1}}}^{\rm{T}} \end{array} $

（5）

定义t=n时刻的星跟踪器测量系C_s到惯导本体系C_g下的投影矩阵A_m为安装阵，定义t=n时刻的星跟踪器的姿态阵为C_sn，则有

$ {\mathit{\boldsymbol{C}}_{{\rm{s}}n}} = {\mathit{\boldsymbol{A}}_{\rm{m}}}{\mathit{\boldsymbol{C}}_{{\rm{g}}n}} $

利用式(5)得到

$ {\mathit{\boldsymbol{C}}_{{\rm{s}}n}}\mathit{\boldsymbol{C}}_{{\rm{sl}}}^{\rm{T}} = {\mathit{\boldsymbol{A}}_{\rm{m}}}{\mathit{\boldsymbol{C}}_{{\rm{g}}n}}\mathit{\boldsymbol{C}}_{{\rm{gl}}}^{\rm{T}}\mathit{\boldsymbol{A}}_{\rm{m}}^{\rm{T}} = {\mathit{\boldsymbol{A}}_{\rm{m}}}{\mathit{\boldsymbol{\hat C}}_{{\rm{g}}n}}\mathit{\boldsymbol{\hat C}}_{{\rm{gl}}}^{\rm{T}}\mathit{\boldsymbol{A}}_{\rm{m}}^{\rm{T}} $

（6）

式(6)说明在{T_j}时间内星跟踪器的测量星矢量可以转换到同一时刻。

对于星跟踪器而言，假设对于某颗恒星S进行观测，设其在J2000.0惯性坐标系和星跟踪器测量系下的矢量分别为V_sⁱ、V_s^s，则有

$ \mathit{\boldsymbol{V}}_{\rm{s}}^{\rm{s}} = {\mathit{\boldsymbol{A}}_{\rm{m}}}{\mathit{\boldsymbol{C}}_{{\rm{ge}}}}{\mathit{\boldsymbol{C}}_{{\rm{ei}}}}\mathit{\boldsymbol{V}}_{\rm{s}}^{\rm{i}} = {\mathit{\boldsymbol{A}}_{\rm{m}}}{\mathit{\boldsymbol{C}}_{{\rm{gn}}}}{\mathit{\boldsymbol{C}}_{{\rm{ne}}}}{\mathit{\boldsymbol{C}}_{{\rm{ei}}}}\mathit{\boldsymbol{V}}_{\rm{s}}^{\rm{i}} $	（7）

式中：C_ne为当地地理坐标系C_n在地心地固坐标系C_e下的投影矩阵；C_ei为地心地固坐标系C_e在J2000.0惯性系C_i下的投影矩阵；Ĉ_ge、C_ei为经过计算可以得到的已知量。

考虑N颗星的星点序列{T_j}，j=1, 2, …, N，定义惯导T_j时刻的姿态误差阵为${\mathit{\boldsymbol{ \boldsymbol{\varDelta} }}_{{\rm{g}}j}}$，根据式(7)得到某颗星某时刻的测量矢量为

$ \mathit{\boldsymbol{\hat V}}_{{\rm{s}}j}^{\rm{s}} = {\mathit{\boldsymbol{A}}_{\rm{m}}}{\mathit{\boldsymbol{C}}_{{\rm{ge}}j}}{\mathit{\boldsymbol{C}}_{{\rm{ei}}j}}\mathit{\boldsymbol{V}}_{{\rm{s}}j}^{\rm{i}} = {\mathit{\boldsymbol{A}}_{\rm{m}}}{\mathit{\boldsymbol{ \boldsymbol{\varDelta} }}_{{\rm{g}}j}}{\mathit{\boldsymbol{\hat C}}_{{\rm{ge}}j}}{\mathit{\boldsymbol{C}}_{{\rm{ei}}j}}\mathit{\boldsymbol{V}}_{{\rm{s}}j}^{\rm{i}} $

（8）

式中，未知量只有安装阵A_m和${\mathit{\boldsymbol{ \boldsymbol{\varDelta} }}_{{\rm{g}}j}}$。定义惯导坐标系下从T₁时刻到T_N时刻的姿态真实值和实测值的变化量分别为${\mathit{\boldsymbol{ \boldsymbol{\varDelta} }}_{{\rm{1}}N}}$、${\mathit{\boldsymbol{ \boldsymbol{\hat \varDelta} }}_{1N}}$，则有

$ \begin{array}{l} {\mathit{\boldsymbol{ \boldsymbol{\varDelta} }}_{1N}} = ({\mathit{\boldsymbol{ \boldsymbol{\varDelta} }}_{{\rm{g1}}}}{{\mathit{\boldsymbol{\hat C}}}_{{\rm{gel}}}}{\mathit{\boldsymbol{C}}_{{\rm{eil}}}}){({\mathit{\boldsymbol{ \boldsymbol{\varDelta} }}_{{\rm{g}}N}}{{\mathit{\boldsymbol{\hat C}}}_{{\rm{ge}}N}}{\mathit{\boldsymbol{C}}_{{\rm{ei}}N}})^\mathit{\boldsymbol{T}}} = \\ {\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} \begin{array}{*{20}{l}} {{\mathit{\boldsymbol{ \boldsymbol{\varDelta} }}_{{\rm{gl}}}}({{\mathit{\boldsymbol{\hat C}}}_{{\rm{gel}}}}{\mathit{\boldsymbol{C}}_{{\rm{eil}}}}){{({{\mathit{\boldsymbol{\hat C}}}_{{\rm{ge}}N}}{\mathit{\boldsymbol{C}}_{{\rm{ei}}N}})}^{\rm{T}}}\mathit{\boldsymbol{ \boldsymbol{\varDelta} }}_{{\rm{g}}N}^{\rm{T}} = }\\ {{\mathit{\boldsymbol{ \boldsymbol{\varDelta} }}_{{\rm{g1}}}}{\mathit{\boldsymbol{ \boldsymbol{\varDelta} }}_{1N}}\mathit{\boldsymbol{ \boldsymbol{\varDelta} }}_{{\rm{g}}N}^{\rm{T}}} \end{array} \end{array} $

式中：${\mathit{\boldsymbol{ \boldsymbol{\varDelta} }}_{\rm{1}}}_N$主要与T_{1和TN时刻的姿态信息有关，而与中间过程无关。因此，有}

$ \begin{array}{l} ({\mathit{\boldsymbol{A}}_{\rm{m}}}{\mathit{\boldsymbol{ \boldsymbol{\varDelta} }}_{{\rm{gl}}}}{{\mathit{\boldsymbol{\hat C}}}_{{\rm{gel}}}}{\mathit{\boldsymbol{C}}_{{\rm{ei1}}}})\mathit{\boldsymbol{V}}_{{\rm{s}}N}^{\rm{i}} = \\ {\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} ({\mathit{\boldsymbol{A}}_{\rm{m}}}{\mathit{\boldsymbol{ \boldsymbol{\varDelta} }}_{{\rm{gl}}}}){{\mathit{\boldsymbol{ \boldsymbol{\hat \varDelta} }}}_{1N}}{{\mathit{\boldsymbol{\hat C}}}_{{\rm{ge}}N}}{\mathit{\boldsymbol{C}}_{{\rm{ei}}N}}\mathit{\boldsymbol{V}}_{{\rm{s}}N}^{\rm{i}} = \\ {\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} ({\mathit{\boldsymbol{A}}_{\rm{m}}}{\mathit{\boldsymbol{ \boldsymbol{\varDelta} }}_{{\rm{g1}}}}){{\mathit{\boldsymbol{ \boldsymbol{\hat \varDelta} }}}_{1N}}\mathit{\boldsymbol{ \boldsymbol{\varDelta} }}_{{\rm{g}}N}^{\rm{T}}\mathit{\boldsymbol{A}}_{\rm{m}}^{\rm{T}}\mathit{\boldsymbol{\hat V}}_{{\rm{s}}N}^{\rm{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} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \begin{array}{*{20}{l}} {({\mathit{\boldsymbol{A}}_{\rm{m}}}{\mathit{\boldsymbol{ \boldsymbol{\varDelta} }}_{{\rm{g1}}}}){{\mathit{\boldsymbol{ \boldsymbol{\hat \varDelta} }}}_{1N}}\mathit{\boldsymbol{ \boldsymbol{\varDelta} }}_{{\rm{g}}N}^{\rm{T}}{\mathit{\boldsymbol{ \boldsymbol{\varDelta} }}_{{\rm{g1}}}}{{({\mathit{\boldsymbol{A}}_{\rm{m}}}{\mathit{\boldsymbol{ \boldsymbol{\varDelta} }}_{{\rm{g1}}}})}^{\rm{T}}}\mathit{\boldsymbol{\hat V}}_{{\rm{s}}N}^{\rm{s}} = }\\ {({\mathit{\boldsymbol{A}}_{\rm{m}}}{\mathit{\boldsymbol{ \boldsymbol{\varDelta} }}_{{\rm{g1}}}})({{\mathit{\boldsymbol{\hat C}}}_{{\rm{gel}}}}{\mathit{\boldsymbol{C}}_{{\rm{eil}}}})\mathit{\boldsymbol{V}}_{{\rm{s}}N}^{\rm{i}}} \end{array} \end{array} $

（9）

这里，令

$ {\mathit{\boldsymbol{A}}_{{\rm{m1}}}} = {\mathit{\boldsymbol{A}}_{\rm{m}}}{\mathit{\boldsymbol{ \boldsymbol{\varDelta} }}_{{\rm{g1}}}},{\mathit{\boldsymbol{ \boldsymbol{\varDelta} }}_{\rm{g}}} = \mathit{\boldsymbol{ \boldsymbol{\varDelta} }}_{{\rm{g}}N}^{\rm{T}}{\mathit{\boldsymbol{ \boldsymbol{\varDelta} }}_{{\rm{g1}}}} $

表示${\mathit{\boldsymbol{ \boldsymbol{\varDelta} }}_{{\rm{g1}}}}$的影响已分别投影到A_m和$\mathit{\boldsymbol{ \boldsymbol{\varDelta} }}_{{\rm{g}}N}^{\rm{T}}$中。设A_m1和${\mathit{\boldsymbol{ \boldsymbol{\varDelta} }}_{\rm{g}}}$分别由绕x、y、z轴的旋转角α、β、γ和θ、δ、φ表示，有

$ \left\{ {\begin{array}{*{20}{l}} {{\mathit{\boldsymbol{A}}_{{\rm{m1}}}} = {\mathit{\boldsymbol{R}}_2}(\beta ){\mathit{\boldsymbol{R}}_1}(\alpha ){\mathit{\boldsymbol{R}}_3}(\gamma )}\\ {{\mathit{\boldsymbol{ \boldsymbol{\varDelta} }}_{\rm{g}}} = {\mathit{\boldsymbol{R}}_2}(\delta ){\mathit{\boldsymbol{R}}_1}(\theta ){\mathit{\boldsymbol{R}}_3}(\varphi )} \end{array}} \right. $	（10）

式中: R₁、R₂、R₃分别为绕x、y、z轴的方向余弦矩阵。基于式(9)，得到转移到同一时刻T₁的星点信息满足：

$ \left\{ {\begin{array}{*{20}{l}} {\mathit{\boldsymbol{\hat V}}_{{\rm{s1}}}^{\rm{s}} = {\mathit{\boldsymbol{A}}_{{\rm{m1}}}}({{\mathit{\boldsymbol{\hat C}}}_{{\rm{gel}}}}{\mathit{\boldsymbol{C}}_{{\rm{eil}}}})\mathit{\boldsymbol{V}}_{{\rm{s1}}}^{\rm{i}}}\\ {{\mathit{\boldsymbol{A}}_{{\rm{m1}}}}{{\mathit{\boldsymbol{ \boldsymbol{\hat \varDelta} }}}_{12}}{\mathit{\boldsymbol{ \boldsymbol{\varDelta} }}_{\rm{g}}}\mathit{\boldsymbol{A}}_{{\rm{m1}}}^{\rm{T}}\mathit{\boldsymbol{\hat V}}_{{\rm{s2}}}^{\rm{s}} = {\mathit{\boldsymbol{A}}_{{\rm{m1}}}}({{\mathit{\boldsymbol{\hat C}}}_{{\rm{gel}}}}{\mathit{\boldsymbol{C}}_{{\rm{ei1}}}})\mathit{\boldsymbol{V}}_{{\rm{s2}}}^{\rm{i}}}\\ \vdots \\ {{\mathit{\boldsymbol{A}}_{{\rm{m1}}}}{{\mathit{\boldsymbol{ \boldsymbol{\hat \varDelta} }}}_{1N}}{\mathit{\boldsymbol{ \boldsymbol{\varDelta} }}_{\rm{g}}}\mathit{\boldsymbol{A}}_{{\rm{m1}}}^{\rm{T}}\mathit{\boldsymbol{\hat V}}_{{\rm{s}}N}^{\rm{s}} = {\mathit{\boldsymbol{A}}_{{\rm{m1}}}}({{\mathit{\boldsymbol{\hat C}}}_{{\rm{gel}}}}{\mathit{\boldsymbol{C}}_{{\rm{eil}}}})\mathit{\boldsymbol{V}}_{{\rm{s}}N}^{\rm{i}}} \end{array}} \right. $

（11）

当N=1时，${\mathit{\boldsymbol{ \boldsymbol{\hat \varDelta} }}_{1N}}$和${\mathit{\boldsymbol{ \boldsymbol{\varDelta} }}_{\rm{g}}}$均退化为单位矩阵，方程组最后一行退化为第1行。上面的非线性方程组可整理成F(x)=0，x为待求的未知量，其大小为3×N，未知数如式(11)中所示，有6个，计算A_m1和${\mathit{\boldsymbol{ \boldsymbol{\varDelta} }}_{\rm{g}}}$即是求解非线性方程组。由于存在多个相互耦合的隐式变量，又考虑到实时性的计算需求，适于选择基于L-M算法的安装阵在线估计。

2 基于L-M算法的安装阵在线估计

L-M算法是一种非线性最小二乘优化算法，采用变步长的“下山”寻优策略，同时具有最优梯度法和高斯牛顿法的优点。考虑到叉乘操作满足:

$ (\mathit{\boldsymbol{A}}{\mathit{\boldsymbol{V}}_i}) \times (\mathit{\boldsymbol{A}}{\mathit{\boldsymbol{V}}_j}) = {\mathit{\boldsymbol{V}}_i} \times {\mathit{\boldsymbol{V}}_j} $

式中：A为3×3的矩阵；V_i、V_j代表 3×1的列矢量。分别计算T₁和T_j时刻2颗星在惯性系下的乘积，结合式(10)，得到3×(N-1)的方程组。其中，第j(j=2, 3, …, N)颗星对应的方程为

$ \begin{array}{l} \mathit{\boldsymbol{F}}(\mathit{\boldsymbol{x}}){|_j} = {\rm{arcsin}}\left\| {\mathit{\boldsymbol{\hat V}}_{{\rm{s1}}}^{\rm{s}} \times ({\mathit{\boldsymbol{A}}_{{\rm{m1}}}}{{\mathit{\boldsymbol{ \boldsymbol{\hat \varDelta} }}}_{1j}}{\mathit{\boldsymbol{ \boldsymbol{\varDelta} }}_{\rm{g}}}\mathit{\boldsymbol{A}}_{{\rm{m1}}}^{\rm{T}}\mathit{\boldsymbol{\hat V}}_{{\rm{s}}j}^{\rm{s}})} \right\| - \\ {\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} {\rm{arcsin}}\left\| {\mathit{\boldsymbol{\hat V}}_{{\rm{s1}}}^{\rm{i}} \times \mathit{\boldsymbol{\hat V}}_{{\rm{s}}j}^{\rm{i}}} \right\| = 0 \end{array} $

（12）

上述3×(N-1)的方程组通过求偏导计算雅克比矩阵不断迭代直至满足收敛条件为止，迭代最终的结果即为待求的6参数。设

$ \mathit{\boldsymbol{P}} = \left[ {\begin{array}{*{20}{l}} \alpha &\beta &\gamma &\theta &\delta &\varphi \end{array}} \right] $

为待求的参数。6个未知数，需要12个方程式进行求解，为保证解算精度，一般取N=10。

L-M算法对应的迭代式为

$ {\rm{d}}\mathit{\boldsymbol{P}} = - {(\mathit{\boldsymbol{J}}_k^{\rm{T}}{\mathit{\boldsymbol{J}}_{\rm{k}}} + {\mu _k}\mathit{\boldsymbol{I}})^{ - 1}}\mathit{\boldsymbol{J}}_k^{\rm{T}}{\mathit{\boldsymbol{F}}_k} $	（13）

式中：dP为待求参数每一步的更新量：P_k为第k步迭代时P的取值，F_k=F(P_k)为P_k处的函数值F(P_k)；μ_k为阻尼因子，J_k为P_k处F(P_k)的雅克比矩阵：

$ {\mathit{\boldsymbol{J}}_k} = \frac{{\partial {\mathit{\boldsymbol{F}}_k}}}{{\partial \mathit{\boldsymbol{P}}}} $	（14）

式(14)的求解这里不便列出。具体的迭代过程如下：

步骤1    给定初值P₀，选择初始阻尼因子μ₀以及步长因子r。

步骤2  计算F_k、J_k，j=0。

步骤3  计算dP。

步骤4  分别计算F^T(P_k+dP)F(P_k+dP)及F^T(P_k)F(P_k)，两者须满足：

$ {\mathit{\boldsymbol{F}}^{\rm{T}}}({\mathit{\boldsymbol{P}}_k} + {\rm{d}}\mathit{\boldsymbol{P}})\mathit{\boldsymbol{F}}({\mathit{\boldsymbol{P}}_k} + {\rm{d}}\mathit{\boldsymbol{P}}) < {\mathit{\boldsymbol{F}}^{\rm{T}}}({\mathit{\boldsymbol{P}}_k})\mathit{\boldsymbol{F}}({\mathit{\boldsymbol{P}}_k}) $	（15）

步骤5  若式(15)成立, 且j=0，则P_k+1=P_k+dP，_μk+1=_μk/r，转步骤3，否则j ≠ 0，转步骤7。

步骤6  若式(15)不成立，则_μk+1=_μkr，转步骤3。

步骤7  如果P_k满足精度，则结束；否则n=n+1，P_k=P_k+1，_μk=_μk+1，回到步骤2。

图 1为安装阵求解过程的示意图，Δ为F对应的阈值，n_max为最大迭代次数。

3 试验数据及分析

本文的全天时星跟踪器与某航天总体单位的激光惯导系统一起工作，针对的应用背景为军用长航时飞机平台和船舶平台。

为验证该系统的性能开展了地面跑车试验，场景有公路、街区和高速公路等，都在白昼环境下进行，如图 2所示。

图 2中，惯导上面通过转接板与全天时星跟踪器捷联安装。在30 s的组合导航周期内，实际的全天时星跟踪器能在白昼条件下平均3 s自动搜索到当前时刻距离当前方位最近的1颗恒星，30 s对应10颗星。具体搜星方式是星跟踪器利用自带的摆镜提供俯仰方向20°左右范围的扫描，同时星跟踪器依靠底部固连的惯导提供偏航方向360°的全周扫描，偏航方向扫描一周对应一个环形的天区条带，扫描一周平均能搜索到20~23颗星，30 s对应扫描半周。星跟踪器为保证星点跟踪的准确性，星表内的星点角距都大于全天时星跟踪器的视场，以此为基础可保证同一时刻出现在视场内的最多只有1颗星。

在实际应用基于L-M算法在线估计安装阵之前，将观测星曝光时刻通过线性插值方式统一到惯导发给星跟踪器的数据的对应时刻，以减少时标不统一带来的测量误差。星跟踪器在累积10颗星形成星点队列后，开始执行基于L-M的算法。每执行一次L-M算法后，按照先进先出的原则，更新一次星点队列。

试验条件如下：

1) 捷联惯导精度：等效陀螺漂移0.01 (°)/h，等效加速度零偏10^-5g。

2) 星跟踪器瞬时视场为2°，观测灵敏度极限为白昼3.5Mv(视星等)，地面标定精度为3″(1σ)；

3) 跑车试验时间约10 h。

实际的跑车试验分别使用了固定安装阵和本文在线估计安装矩阵的方法。定义星跟踪器的单星投影误差为多个观测星矢量转移到同一时刻捷联惯导本体坐标系中的误差，其主要包括安装矩阵误差和星跟踪器的测量误差，采用星跟踪器的单星投影误差衡量星跟踪器的输出精度。

图 3显示了实际12 h跑车试验中星跟踪器固定安装阵时的单星投影误差曲线，此图是中值滤波滤掉野值和跳点后的结果。图 3中，单星x方向投影误差为7.251 3″(1σ)，y方向投影误差为7.323 4″(1σ)，均大于星跟踪器的地面标定精度3″(1σ)一倍以上，且存在周期性波动，并随时间而发散。图 3显示，如果不进行安装阵的在线估计，单星定位精度将无法得到保证。

图 4为实际试验中安装阵3个轴角实时估计后的曲线。与图 4的安装阵的实时估计相对应，图 5为实时估计出的惯导3轴姿态误差曲线。

图 5显示，在组合导航的过程中，虽然惯导误差在逐渐放大，但并不影响安装阵的实时估计。如图 6所示，也不影响对单星投影误差的修正。

图 6为使用实时估计安装阵方法后星跟踪器的单星投影误差曲线。其中，单星x方向投影误差为3.132 7″(1σ)，y方向投影误差为3.718 3″(1σ)，与星跟踪器地面标定精度一致。

对比图 3和图 6，使用该方法后全天时星跟踪器提供的指向精度提高了1倍，提供的单星指向精度与星跟踪器的单星指向精度3″接近，基本上消除了安装阵和惯导姿态误差对星跟踪器输出测星精度的影响。

图 7为对应的安装阵3个轴角的收敛次数统计曲线。迭代次数最多为38次，最少为2次，平均迭代9.46次，平均计算时间为3.8 ms < 5 ms。

图 8为组合导航静态条件下导航1 h的对比试验数据，实线为安装阵固定时的结果，虚线为本文基于L-M算法估计安装阵后的结果。组合导航算法在静态条件下使用本文算法前后的位置误差分别为优于519 m/h和优于192 m/h。

4 结论

1) 针对全天时星光惯性组合导航系统，提出了一种基于L-M算法的在线实时估计安装阵方法，能有效消除安装阵偏移和惯导姿态误差累积对星跟踪器输出精度的影响。

2) 跑车试验结果表明采用该算法后，惯导的输出测星精度由7.251 3″(1σ)变为优于3.8″(1σ)，在线估计平均时间 < 5 ms，满足星光惯性组合导航的实时性需求。

3) 经过在线估计安装阵后，本文方法能够将静态条件下组合导航的定位精度提升1倍以上。动态跑车的试验结果可以满足总体技术指标要求。