基于空间目标异步观测的惯导误差快速确定.pdf

第32卷第1期 中国惯性技术学报 Vol.32 No.1 2024 年 01 月 Journal of Chinese Inertial Technology Jan.2024 收稿日期：收稿日期：2023-07-18；修回日期：修回日期：2023-12-10 基金项目：基金项目：国家自然科学基金（U21B6001）通讯作者：通讯作者：杨静（1975），女，副教授，硕士生导师，从事多源信息融合、故障诊断与容错、组合导航技术。文章编号：文章编号：1005-6734(2024)01-0016-11 doi.10.13695/ki.12-1222/o3.2024.01.003 基于空间目标异步观测的惯导误差快速确定基于空间目标异步观测的惯导误差快速确定 杨 静1，王 栋1，熊 凯2（1.北京航空航天大学 自动化科学与电气工程学院,北京 100191；2.北京控制工程研究所 空间智能控制技术重点实验室,北京 100094）摘要：摘要：针对惯性/天文组合导航系统中大的初始状态误差影响惯导误差收敛速度的问题，在星相机观测具有先验位置信息的有限空间目标的辅助下，提出一种基于目标-恒星角距异步测量的惯导误差在线快速确定方法。首先，在星相机光轴旋转角度和视场角大小受限的情况下，设计了通过异步照相观测方式获取有效空间目标参考信息的方案；其次，在利用惯导误差状态传播模型实现异步测量信息同步处理的基础上，构建基于空间目标与恒星之间角距的非线性最小二乘优化模型，避免了星相机的光轴扰动和安装误差对测量精度的影响；最后，基于高斯牛顿法设计了两轮迭代优化估计惯导位置误差和速度误差的方法。蒙特卡洛仿真结果表明，所提方法利用对空间目标和恒星的有限观测信息，可以有效估计惯导位置误差和速度误差，在初始位置误差约十千米量级的情况下，可以估计补偿约 97.73%的位置误差以及 66.25%的速度误差，优化求解误差参数的计算耗时为 0.0160 s。关关 键键 词：词：天文导航；异步观测；非线性最小二乘；高斯牛顿法 中图分类号：中图分类号：V249.32 文献标志码：文献标志码：A Fast determination of inertial navigation error based on asynchronous observation of space targets YANG Jing1,WANG Dong1,XIONG Kai2(1.School of Automation Science and Electrical Engineering,Beihang University,Beijing 100191,China;2.Science and Technology on Space Intelligent Control Laboratory,Beijing Institute of Control Engineering,Beijing 100094,China)Abstract:Aiming at the problem that large initial state error will affect the convergence rate of inertial navigation error in the inertial/celestial navigation system,assisted by star camera observing limited space target with prior position information,a kind of online fast determination method of inertial navigation error based on asynchronous measurement of target-star angular distance is proposed.Firstly,under the condition that rotation angle of optical axis and field of view of star camera are limited,the scheme of acquiring effective reference information of space target through asynchronous photographic observation is designed.Secondly,on the basis of synchronization processing of asynchronous measurement by exploiting state propagation model of inertial navigation error,nonlinear least square optimization model based on angular distance between target and star is constructed which can avoid influence of optical axis disturbance and installation error.Finally,a two-round iterative optimization method is designed to estimate the position error and velocity error of inertial navigation.Monte Carlo simulation results show that the proposed method can effectively estimate the position error and velocity error of inertial navigation by exploiting the limited observation information of space targets and stars,and can compensate the initial position error by 97.73%and the initial velocity error by 66.25%when the initial position error is about ten kilometers.The computation time of the optimization solution error parameter is 0.0160 s.第 1 期 杨静等：基于空间目标异步观测的惯导误差快速确定 17 Key words:celestial navigation;asynchronous observation;nonlinear least square;Gauss-Newton method 惯性导航系统自主性、抗干扰性较强且导航信息完整，但导航误差随时间积累1。作为一种重要的天文导航敏感器，星相机可以同时对恒星以及具有一定亮度和轨道先验信息的空间目标如人造卫星等进行照相观测2，通过星图匹配与图像处理提取恒星的星光矢量以及恒星与空间目标的角距，进而还可以采用轴系定位法或天文定位法得到空间目标的方向矢量，将其作为空间基准信息源可以辅助校正惯导解算误差，以满足高精度长航时的自主导航定位需求3。由于星相机成像曝光以及星图和目标的识别处理需要消耗一定的时间4，并且可观测恒星和空间目标的数量与光学镜头的视场角大小、探测灵敏度等性能指标以及被观测对象的分布和运动特性密切相关，导致基于照相观测星相机的天文导航系统的数据更新率较低5。此外，星相机提供的是与载体位置间接相关的测量信息，而对载体速度的可观性弱，速度误差难以被有效地估计出来。如何有效利用数量和精度有限的测量信息，实现对惯导误差参数的快速、准确估计是目前需要解决的问题。目前，常用的惯导误差估计方法分为两类：基于状态估计的方法和基于参数估计的方法。其中，基于状态估计的方法包括卡尔曼滤波6、扩展卡尔曼滤波7,8、无迹卡尔曼滤波9-11、容积卡尔曼滤波12,13等方法。文献14采用标准卡尔曼滤波器融合里程计、星敏感器和惯性器件的测量信息，实现动态的初始对准和误差校正。文献15和文献16将待估误差项扩维到状态向量中，然后采用不同的自适应容积卡尔曼滤波算法在线递推估计导航误差参数。文献17提出一种确定月球探测车初始位置和姿态以及惯性器件误差的方法，在月球车保持静止时利用恒星高度角、星光矢量以及水平速度误差建立测量模型，再采用无迹卡尔曼滤波算法估计并校正惯导误差。在以上基于状态估计的方法中，往往被估误差参数需要一定的时间才能稳定收敛，并且当初始误差参数较大时可能导致滤波收敛时间增长，严重情况下甚至可能造成滤波发散，进而难以快速有效地实现对误差参数的准确估计与补偿。此外，基于状态估计方法的估计精度还依赖于系统模型及噪声统计特性的准确程度。在进行参数估计时，最小二乘是广泛应用的准则之一。当待估参数与系统测量建立的是线性模型时，即可用线性最小二乘法获得误差参数的估计值，该方法无需测量噪声先验信息，实现简单，快速性较好18；当待估参数与系统测量是非线性关系时，可以基于最小二乘准则建立关于待估参数的非线性优化模型，通常采用高斯-牛顿19（Gauss Newton,GN）、列文伯格-马夸尔特20（Levenberg-Marquardt,LM）等非线性优化方法进行迭代求解，以保证参数估计的精度。在导弹、高超声速飞行器等应用平台，高速高动态特征将造成严重的气动光学效应，同时由于受到光学镜头观测视野的限制，获得的有效天文测量信息有限。针对大的初始位置和速度误差影响组合导航估计误差收敛速度，严重者甚至导致滤波发散的问题，本文结合异步测量同步处理思想，提出一种非线性最小二乘优化方法，通过融合有限的空间目标观测信息快速确定惯导误差参数。1 模型建立 1 模型建立 1.1 系统状态方程 由于初始导航参数误差、加速度计偏置、陀螺漂移等误差的存在，惯导解算中的积分运算将使导航误差随时间不断积累增大。通常，建立惯导小误差方程来描述载体状态误差和惯性器件主要误差分量的动态变化情况，并可用于实现对误差状态的递推估计。在发射点惯性系（简称发惯系）下采用导航状态误差和惯性器件的主要误差分量构成导航系统的状态向量为：T=liliccXVP(1)其中，为平台误差角向量；licV为速度误差向量，上标li表示发惯系，下标c表示载体坐标系；licP为位置误差向量；为三轴陀螺常值漂移误差向量；为三轴加速度计常值偏置误差向量；上标T表示转置。将天文/惯性组合导航系统的状态方程表示为：()()()()()tttttXFXGW(2)其中，t为时间；()tF为系统状态矩阵；()tG为噪声因数矩阵；()tW为系统噪声。将式(2)离散化可得：|11,1()kkk kkk ktXXXw(3)其中，下标k对应当前时间kt；,1k kw为 ttGW离散化后的过程噪声；|1k k表示从1kt到kt的状态转移矩阵，即|11,k kkkt t，可以按照式(4)计算得到。23|1112!3!k kTTTFFFI(4)其中，T为离散周期。通过式(3)可以建立不同时刻系统状态之间的转18 中国惯性技术学报 第 32 卷 换关系，并利用惯性传感器测量的角速度和加速度进行惯导解算。在此基础上，根据式(5)将1k 时刻的状态估计值1kX利用状态转移矩阵递推得到k时刻的状态估计值|1k kX以完成时间更新。|1|11k kk kkXX(5)其中，上标“”表示估计值。1.2 优化模型 构建准确的测量模型以有效利用测量信息是提升状态估计精度的一种重要途径。在基于照相观测星相机的天文/惯性组合导航系统中，关于空间目标的可用测量信息主要有两种：视线方向矢量测量以及与恒星间的角距测量。天文观测示意图如图1所示，在工作过程中，星相机同时对至少一颗空间目标和若干颗背景恒星进行照相观测，通过星图识别可获得恒星矢量以及空间目标的视线方向矢量，进而得到两者之间的角距。恒星恒星空间目标空间目标视线方向矢量视线方向矢量恒星矢量恒星矢量视场视场轨迹轨迹空间目标-恒星角距空间目标-恒星角距 图 1 天文观测示意图 Fig.1 Schematic diagram of celestial observation 目标视线方向矢量测量提供的是与载体位置直接相关的信息，信息利用率较高，据此建立的测量模型非线性度较小，但基于传统轴系定位法提取的目标方向矢量难以避免星相机轴扰动误差的影响，而坐标系转换的过程还会引入安装误差的影响。因此，当星相机存在较大的安装误差和光轴扰动误差时，基于目标视线方向矢量的测量模型精度有限。与之相比，基于角距的测量模型非线性度较高，但是由于角距信息不依赖于坐标系，故采用基于角距建立测量模型的方式可以有效避免星相机轴扰动和安装误差的影响。通常情况下，星相机视野中有足够数量的可观测恒星，对其进行照相观测可以提取较高精度的载体姿态四元数，以辅助校正惯导姿态误差，因此，基于卡尔曼滤波的状态估计方法可对姿态误差进行较准确的估计与补偿，并且由于有关载体姿态的测量信息（恒星矢量）比较充足，即使在较大的初始姿态误差下，对应姿态误差分量的滤波估计也可以快速收敛，无需在初始误差确定阶段对载体姿态误差进行优化与补偿。由于空间目标的数量和分布空间的限制，通常情况下星相机仅能观测到一颗可观性较好的目标。此外，通过观测空间目标得到的测量信息仅与载体位置间接相关，与载体速度无法建立函数关系。如果直接采用滤波方案对载体位置和速度误差进行估计，通常需要一定时间才能稳定收敛，并且当初始位置和速度误差参数较大时可能导致滤波收敛时间增长，严重情况下甚至可能造成滤波发散。因此，本文根据空间目标与恒星之间角距的测量模型，以导航位置和速度误差为待估参数，基于异步角距测量建立优化模型，通过非线性最小二乘方法迭代求解惯导初始误差，以保证后续滤波方案的快速收敛。建立待优化的目标函数F为：2,011()(|)2sik kmlik igk iins gk ijigFttt(6)其中，im为k it时刻测量得到的恒星矢量与空间目标方向矢量之间的角距数目；skt为启动照相观测获取角距测量辅助导航误差确定的起始时刻；kt为结束时刻；jt为待估计的导航误差所对应的时刻，该时刻的导航误差将被用来对所有异步测量进行时间对准的处理，,sjk k；k i表示k it时刻是否获取有效照相观测的指示变量，1k i时为存在有效观测，否则0k i；()gk it为星相机在k it时刻测量得到的空间目标与第g颗观测恒星之间的角距，即：T()()arccos()()k igk igk ik ittttLSL(7)其中，()gk itS为k it时刻第g颗观测恒星的星光方向矢量；()k itL为k it时刻空间目标的方向矢量；为角距测量误差。,(|)liins gk ijtt表示根据jt时刻的惯导误差向量解算得到的k it时刻第g颗观测恒星矢量和空间目标矢量间的角距估计值，下标ins表示该参数估计值是根据惯导解算结果计算得到的。具体计算方法如下：1）将jt时刻的惯导误差向量()jtX经过状态转移得到k it时刻的(|)k ijttX；2）利用状态向量中的位置误差分量(|)lick ijttP对k it时刻的惯导解算位置,()c insk itP的误差进行修正，可得载体位置：,|lic insk ic insk ick ijttttPPP(8)3）利用,()c insk itP和k it时刻的空间目标位置估计()litk itP计算目标视线方向单位矢量估计值：第 1 期 杨静等：基于空间目标异步观测的惯导误差快速确定 19 ,litk ic insk ik ilitk ic insk itttttPPLPP(9)4）在发惯系下，计算角距估计值：T,(|)arccos()liliiins gk ijigk ik ittttC SL(10)其中，右上标i表示地心惯性坐标系；liiC为从地心惯性坐标系到发惯系的转换矩阵，由载体发射点的初始经度0和初始纬度0L计算得到。由式(3)可推导出k it和jt时刻的状态向量满足：|,k ik i jjk i jXXw(11)其中，1|k i jj k i为jt到k it的状态转移矩阵；,k i jw为等效过程噪声。令|Pk i j为状态转移矩阵|k i j中与位置分量对应的第79行，再将其根据该列所对应的状态分量写成分块矩阵的形式：|Pppvppppk i jk i jk i jk i jk i jk i j(12)其中，上标、v、p、分别表示|Pk i j中与平台误差角、速度、位置、陀螺漂移和加速度计零偏对应的分量。根据式(11)和式(12)，进一步推导出由jt时刻的状态向量jX得到k it时刻位置误差估计为：|(,)()()()+()()(,)liPck ijk i jjjk ipplipvlipk i jcjk i jcjk i jjppPk i jjk i jjk i jjk ittt ttttttt tPXwPVw(13)式(13)给出了由jX的各个状态分量以及系统噪声解析表示的(|)lick ijttP，据此可以建立描述与空间目标的多个异步观测时刻相对应的位置误差与待估计jt时刻的导航误差之间关系的模型。在忽略各个状态分量之间相关性的前提下，由式(14)可得各分量协方差阵之间的关系。TT|TT|TT|()()()()()(,)pppppvpvpk ijk i jpjk i jk i jvjk i jppppk i jjk i jk i jjk i jppPPk i jjk i jk i jjk ik i jtttttttt tPPPPPPQ(14)其中，()jtP表示jt时刻物理量()的协方差阵；(,)jk it tQ表示由jt时刻到k it时刻的过程噪声协方差阵。由于在确定误差的短时段内，小的姿态误差和器件误差对系统的影响有限，而且期间获取的角距信息也难以提供足够的信息对这些误差分量进行有效估计，因此，在初始阶段仅考虑将主要导航误差分量()licjtP和()licjtV作为待确定参数j，其他分量均作为 等 效 模 型 误 差 进 行 处 理。利 用 式(10)建 立(|)lick ijttP和,(|)liins gk ijtt之间的关系，并将其代入式(6)建立非线性的优化模型。通过优化求解该模型可得参数估计值为：2,011argmin()(|)2sijk kmlijk igk iins gk ijigttt(15)2 2 高斯牛顿法高斯牛顿法 当目标函数是非线性最小二乘问题时，可以采用GN、LM等非线性优化方法进行求解。与传统牛顿法相比，GN算法利用目标函数的一阶泰勒展开项近似其二阶泰勒展开项，通过计算拟海塞矩阵降低了计算量，但当系数矩阵奇异或病态时存在增量不稳定甚至算法不收敛的问题。LM算法通过引入信赖域在一定程度上解决了系数矩阵的奇异和病态问题，但算法的收敛速度较慢。由于本文建立的优化模型具有较好的性质，通过合理设计空间目标的观测方式可以避免出现系数矩阵奇异或病态问题，因此，采用GN算法来实现优化模型的快速求解。下面给出GN算法的实现流程：1）设置迭代初始值0 x；2）对于第i次迭代，计算雅克比矩阵iJ x以及目标函数()if x，求解增量方程得到增量ix：TTiiiii J xJ xxJ xfx(16)3）如果ix小于增量阈值则停止迭代，否则更新迭代值：1iiixxx(17)4）如果已到达最大迭代次数或者目标函数值小于阈值则停止迭代，否则返回步骤2)。3 3 基于角距异步观测的惯导误差快速确定方法基于角距异步观测的惯导误差快速确定方法 本小节详细介绍本文方法各模块的实现流程。3.1 空间目标异步测量信息的采集 为了保证误差参数估计的有效性以及快速性，在天文导航启动初期，依靠载体姿态机动或通过旋转星相机安装平台的方式改变星相机光轴的方位角，使其指向不同范围的天区以保证待估参数具有较好的可观测度，并且在同一指向下星相机以较高的频率对视野内的空间目标和恒星连续进行多次观测以获取充足的异步观测信息。考虑到星相机的指向次数越多，导航误差的传播时间越长，异步测量信息的有效性下降，本文将光轴20 中国惯性技术学报 第 32 卷 的指向次数设置为两次，下面给出异步测量采集的具体实现方式：1）设置同一指向下高频观测频率以及观测次数。2）天文导航启动后，在星相机光轴的第一个指向下，以预设的观测频率对视野内的空间目标和恒星连续进行多次观测，通过成像曝光和星图识别等过程获得恒星的星光矢量和目标卫星的视线方向矢量，进而得到恒星和目标卫星之间的角距测量信息，并记录每个观测时刻时间更新过程中状态转移阵的变化。3）依靠载体姿态机动或通过旋转星相机安装平台的方式改变星相机的光轴指向。如果依靠载体的姿态机动来改变光轴指向，则光轴转角的大小取决于载体姿态机动的大小；如果通过旋转星相机安装平台的方式改变光轴指向，光轴转角可以设置在40 80 范围内，使得观测信息尽可能对状态形成较好的几何约束。4）在星相机光轴的第二个指向下，以预设的观测频率对视野内的空间目标和恒星连续进行多次观测，采集空间目标角距测量信息并记录每个观测时刻时间更新过程中状态转移阵的变化。综上，观测信息采集处理的时间示意图如图2所示。其中，1t为相邻两次惯导解算之间的时间间隔；2t为相邻两次时间更新之间的时间间隔；3t为同一指向下相邻两次天文观测的时间间隔；4t为星相机光轴改变指向的时间间隔。211tnt，322tnt，433tnt，1n、2n、3nN，N表示正整数。惯导解算时刻惯导解算时刻时间更新时刻天文观测时刻时间更新时刻天文观测时刻第一个指向下的观测第二个指向下的观测第一个指向下的观测第二个指向下的观测改变光轴转角改变光轴转角1t2t3t4t 图 2 观测信息采集处理时间示意图 Fig.2 Schematic diagram of time for observation information acquisition and processing 3.2 基于GN算法的优化模型求解 考虑到大的初始位置误差会导致较大的线性化误差，进而对模型线性化和观测同步处理过程产生严重影响。同时，还要避免不准确的速度误差估计对位置估计精度的影响，因此，采用了两轮优化的求解方案。在第一轮优化求解过程中，仅以位置误差作为待估参数建立优化模型，采用GN算法迭代计算出位置误差估计值；在第二轮优化求解时，利用第一轮估计结果修正位置误差，并建立以残余位置误差和速度误差作为待估参数的优化模型，从而对速度误差和残余位置误差进行同步求解。采用经过两轮估计得到的位置误差和速度误差对惯导误差进行修正，进而可以为后续的惯导解算或组合导航提供更准确的数据初值。算法的总体结构图如图3所示。第一轮优化建立以位置误差为待估参数的优化模型基于GN算法的优化求解第二轮优化建立以位置和速度误差为联合待估参数的优化模型基于GN算法的优化求解惯导解算时间更新空间目标角距异步测量获取天文观测解算位置、速度、姿态解算位置状态转移阵角距、星光矢量位置误差估计值位置和速度误差估计值IMU比力、角速度 图 3 算法总体结构图 Fig.3 The overall structure diagram of the algorithm 下面给出两轮优化求解的具体流程。1）第一轮优化 建立以位置误差为待估参数的优化模型，即待估参数向量j由jt时刻的位置误差()licjtP构成。由式(13)可得：|(|)(),(|)(),lipplick ijk i jcjlilick ijcjtttjkitttjki PPPP(18)基于式(10)建立如式(15)所示的优化目标，然后用GN算法对位置误差参数进行迭代估计，具体步骤为：设置导航误差参数的迭代初值,0j、增量变化阈值,minj、目标函数阈值minF以及最大迭代步数maxk；如果当前迭代步数k已达到最大迭代步数则停止迭代，否则进行步骤；记第k步迭代中参数估计值为,j k，由于第一轮 优 化 中 仅 以 位 置 误 差 作 为 待 估 参 数，即,(,)lij kcjt k P，(,)licjt kP表示第k步迭代得到的jt时刻的位置误差估计值。由式(13)可得：|(|,)(,),(|,)(,),lipplick ijk i jcjlilick ijcjtt kt kjkitt kt kjki PPPP(19)其中，(|,)lick ijtt kP为在第k步迭代中由(,)licjt kP经过状态转移到k it时刻位置误差估计值。k it时刻的惯导解算位置误差修正为：,(,)()(|,)lic insk ic insk ick ijtktttkPPP(20)第 1 期 杨静等：基于空间目标异步观测的惯导误差快速确定 21 其中，,(,)c insk itkP表示在第k步迭代中修正后的惯导解算位置。根据式(9)和式(10)依次确定k it时刻空间目标的单位方向矢量,k itkL以及其与第g颗观测恒星间的角距,(|,)liins gk ijtt k。据此，计算第k步迭代的目标函数值：2,011()()(|,)2sik kmlij kk igk iins gk ijigFttt k(21)如果目标函数值,()j kF 小于阈值minF则停止迭代，否则继续下一步骤；计算第k步迭代中的目标函数关于当前待估参数,j k 的雅克比矩阵,()j kJ。,()=jj kj kjFJ(22)构造正规增量方程：T,()()=()()j kj kj kj kj kFJJJ(23)其中，,j k 为第k步迭代的参数增量，若增量小于增量变化阈值,minj 则停止迭代，否则基于增量更新误差参数估计值,1j k。,1,j kj kj k (24)将,1j k 作为第1k次的参数估计值，返回步骤继续下一步迭代。2）第二轮优化 为了进一步提高位置误差的估计精度，并且利用载体位置变化中隐含的速度信息，实现对速度误差的估计。以第一轮优化结果作为位置误差初值，建立以位置误差()licjtP和速度误差()licjtV为联合待估参数j的优化模型。由式(13)可得：|(|)()(),(|)(),lipplipvlick ijk i jcjk i jcjlilick ijcjttttjkitttjkiPPVPP(25)将式(25)代入式(8)，再依次代入式(9)和式(10)，最终建立形如式(15)所示的优化模型，再采用与第一轮优化相同的方法进行处理。4 4 可观测性分析方法可观测性分析方法 系统的可观测性是指在有限的时间区间内根据测量信息推断系统状态的能力。通过对系统状态的可观测性进行表征、判定以及量化评估，可以为系统状态的准确完备估计提供方向。Fisher信息矩阵、可观测性矩阵的条件数和范数等均可以作为系统可观测性的量化指标。目前常用的可观测性分析方法包括基于估计误差方差阵的特征值特征向量分析法、基于Fisher信息矩阵的克拉美劳方差下界（Cramer-Rao Lower Bounds,CRLB）分析法、以及基于可观测性矩阵奇异值分解的分析方法。前两种方法仅适用于滤波框架下的可观测性分析，而本文是基于参数优化的思想确定惯导误差，因此采用基于可观测性矩阵奇异值分解的方法对误差可观测性分析。下面给出误差可观测性分析的具体实现过程。由于仅对位置误差和速度误差进行优化，此处仅以发惯系下载体三维位置误差和三维速度误差构成六维状态向量T XPV。在k时刻，空间目标与第i颗 恒 星 之 间 的 角 距 实 际 测 量 值 记 为ik，1,2sin，sn为可观恒星数量。空间目标与第i颗恒星之间的角距估计值为：T,1,1,1arccost kc kk kiik kkt kc kk kPPPSPPP (26)其中，ikS为k时刻第i颗恒星的星光矢量测量值；,t kP为k时刻空间目标位置估计值；,c kP为k时刻载体位置估计值；,1k kP为k时刻状态估计值,1k kX中前三维对应载体位置误差的状态分量。基于式(7)和式(26)，建立测量模型为：11,1,1,1,1ssk kkiikkk kkkkk kknnk kk Zh X(27)其中，,1k kk、分别为k时刻角距估计值和实际测量值向量；h为非线性的测量方程；k 为测量噪声。将式(27)在状态估计值,1k kX处进行泰勒展开，得到线性化的测量方程为：,1,1,1k kkk kkk kkXXhZh XXXX (28)将式(5)和式(28)改写成离散线性系统模型的形式：,11kk kkkkkkXXZH X (29)其中，状态向量nkXR；状态转移矩阵,1nk knR；测量向量mkZR；线性化后的测量阵nkmHR；6n为状态维数，smn为测量维数；k 为包含线性化误差的测量噪声。根据式(28)，测量阵的具体形式为：,1,1,1,1,k kk kk kk kk XXPPVVPPhhhHXPVhLLP0(30)22 中国惯性技术学报 第 32 卷 其中，L 为空间目标视线方向关于载体位置误差的非线性函数，具体形式见式(8)和式(9)。由式(30)可知，仅利用当前k时刻的角距测量信息，测量阵是不满秩的，速度误差完全不可观测。根据分段线性化的分析方法，将所有的异步测量信息均表示为同步处理时刻状态kX的函数：111,1111,1kkkkkk kkkkkkkZH XZHXZH X (31)系统的可观测性矩阵kO可表示为：11,111,1kkkk kkH OHH (32)由式(29)可知，在异步测量信息充足的情况下，通过状态转移矩阵的同步处理，kO可以满足列满秩的可观性条件，即可以通过测量值对状态进行求解。如果可观测性矩阵是列满秩的，则k时刻系统是可观测的，可以将系统可观测度定义为：1condkkO (33)其中，cond表示求条件数。当可观测性矩阵条件数较大时，说明其为病态矩阵，那么在相同测量误差条件下得到的状态估计误差较大，难以通过测量信息准确估计系统状态，系统可观性较差；反之，当kO的条件数较小时，则系统的可观测度较好。通过以上定义的可观测度仅能评估系统整体的可观测性，为了对每一维状态的可观测性进行分析，对可观测性矩阵进一步进行奇异值分解，再根据奇异值的大小进行分析。对于整个系统而言，奇异值越小则系统的可观测度越差。对于每一维状态而言，其所包含的最小奇异值在变量中所占成分越大，则该维状态的可观测度越差。5 5 仿真实验与分析仿真实验与分析 5.1 仿真验证 以安装有星相机的高速运载体为仿真研究对象，采用天文/惯性组合导航方式，在载体助推段结束后启动星相机参与导航。天文测量精度为3，星相机的观测频率为0.5 Hz，星相机安装误差和光轴扰动误差为3，星相机的光轴与安装基座法向轴的夹角为20。在每个观测时刻识别一个空间目标和五颗背景恒星，在同一光轴指向下以10 Hz连续观测10次，相对基座改变光轴指向1次。陀螺精度水平为0.02/h，加速度计精度水平为10 g。根据先验信息获得的空间目标的位置精度为10 m。载体的初始位置误差为1000 m，初始速度误差为100 m/s。若无特殊说明，以上为后续实验中的缺省仿真参数。在上述仿真条件下，进行50次蒙特卡洛仿真。采用迭代后的位置剩余误差和速度剩余误差的均方根统计值作为评估本文方法性能的指标，本文将剩余误差定义为传播到最后一个测量时刻的导航误差与经迭代优化后估计的导航误差参数之差。当最大迭代步数为10步时，两轮优化过程中的目标函数值的收敛曲线分别如图4和图5所示，不同阶段下的导航误差统计结果见表1。由图4可知，在以位置误差为待估参数的迭代优化过程中，目标函数值经过五步迭代后就趋于收敛。由图5可知，在以位置误差和速度误差为联合待估参数的迭代优化过程中，目标函数经过两步迭代后就趋于收敛，这展现了高斯牛顿法具有快速收敛的优势，且将第一轮优化得到的结果作为第二轮优化的位置误差初值，加快了第二轮优化的收敛速度。图 4 第一轮优化过程中目标函数收敛曲线 Fig.4 Objective function in the first round of optimization 图 5 第二轮优化过程中目标函数收敛曲线 Fig.5 Objective function in the second round of optimization 由表1可知，在初始误差较大（位置误差为十千米量级）的情况下，本文方法可以补偿大约97.73%的初始位置误差以及66.25%的初始速度误差。在计算机第 1 期 杨静等：基于空间目标异步观测的惯导误差快速确定 23 配置为Intel i7-6700、安装系统为win7版本、CPU主频为3.40 GHZ的条件下，优化求解误差参数的耗时为0.0160 s，这体现了所提方法具有一定的快速性。当同一指向下的观测次数降为四次时，优化求解误差参数的耗时降为0.0109 s，即异步测量信息数据量的减少有助于缩短求解误差参数的耗时。仿真统计结果表明，本文方法可以准确且快速地确定启动天文导航初期惯导解算的位置误差和速度误差，并且经过第二轮优化对速度误差的有效估计进一步提升了位置误差估计的准确程度。表 1 不同阶段的导航误差 Tab.1 Navigation errors under different stages 导航参数 位置/m 速度/(m/s)阶段 未补偿 第一轮 第二轮 未补偿第二轮误差 11247.5 312.2 255.1 173.4058.52补偿百分比-97.22%97.73%-66.25%5.2 误差可观测性分析 利用第3小节给出的基于可观测性矩阵奇异值分解的分析方法，从误差可观测性角度验证本文方法的有效性。基于5.1节的实验条件，分别在改变和固定星相机光轴指向情况下，将可观测性矩阵的秩以及可观测度随观测次数的变化情况统计在表2中，并将经过奇异值分解后可观测性矩阵的奇异值和对应的右奇异向量直方图绘制在图6和图7中，用于对每一维状态误差的可观测性改善情况进行分析。按照对应右奇异向量的奇异值从大到小的顺序对各子图进行排列，图中为对应的奇异值，横坐标16维依次对应x、y、z轴的位置误差和速度误差。表 2 不同观测方式下的可观测性分析 Tab.2 Observability analysis under different observation modes 星相机光轴 指向情况 观测次数 可观测性 矩阵的秩 可观测度 改变指向 1 4 0 2 6 2.10E-03 4 6 4.79E-03 6 6 6.04E-03 8 6 7.66E-03 10 6 8.34E-03 固定指向 1 4 0 2 6 7.03E-06 4 6 2.26E-05 6 6 3.51E-05 8 6 4.78E-05 10 6 6.46E-05 由表2可知，当观测次数为1时，无论是否改变星相机的光轴指向，可观测性矩阵都不是列满秩的，系统整体是不可观测的。当观测次数大于1时，可观测性矩阵是列满秩的，且系统可观测度随着观测次数的增加而不断改善。此外，与固定光轴指向相比，改变光轴指向下的可观测度大约提升了两个数量级。由图6可知，在固定光轴指向的情况下，奇异值由高到低依次对应的是y轴位置误差、y轴速度误差、x轴位置误差、x轴速度误差、z轴位置误差、z轴速度误差，这是由于从固定光轴指向得到的测量信息无法对每个轴向的状态误差形成完备的几何约束，从而导致各轴向的可观测度相差较大。与图6相比，在改变光轴指向的情况下，图7中y轴和z轴的位置误差及速度误差对应的奇异值明显增大，说明通过改变光轴指向可以对各轴状态误差形成比较完备的几何约束，因此可观测度明显提升。=2.7339e-0500.51=4.7128e-0600.51=3.17e-0600.51123456=4.555e-0712345600.51123456=2.2958e-0800.51123456=1.8726e-1