适应多测站滤波定轨系统的可观测度分析方法

,,,

(1.航天工程大学 研究生院，北京 101416; 2.航天工程大学 光电装备系，北京 101416)

0 引言

可观测性这一概念是由Kalman为了解决线性系统的相关问题而提出的，如果系统的状态能被过去的观测唯一确定，则该系统为可观测的[1]。但航天工程领域，各种系统均为非线性的，例如自主导航及滤波定轨系统。对于非线性时变系统的可观测性，还没有一个统一而严密的定义。

Lee和Dunn提出的李函数准则是目前非线性系统可观测性判定的常用准则[2]。仅判断系统的可观测性在实际应用中作用并不大，所以引入可观测度的定义，能将系统的可观测性大小定量的体现出来。针对这一问题，国内外先后提出多种非线性系统的可观测度的判断方法，其中可观测矩阵SVD分解[3-4]和基于条件数的可观测度计算方法[5-8]最为常见，然而这两种方案均未考虑量测噪声对系统可观测性造成的影响，李恒年教授及孙仲康教授在文献 [9-10]中均提出了将可观测矩阵和量测噪声协方差矩阵结合的可观测度判定方法，所求得的可观测度受噪声影响较大。

在多站单测量体制的实际应用条件下，对滤波定轨系统进行可观测性分析。针对该种弱观测条件下可观测矩阵条件数过大时可观测度判断不准确的这一情况，分析了误差形成的原因并提出了相关解决方案。建立卫星动力学状态方程与多站测速的量测方程，并求出其观测矩阵，采用容积卡尔曼滤波(CKF)算法作为定轨算法，利用Matlab软件对多测站滤波定轨系统的可观测度进行仿真。通过比较不同测站组合情况下定轨精度与可观测度的关系。仿真结果表明，改进后可观测度判断方法相较于传统的基于条件数判断方法，能更准确判断系统的定轨精度与可观测度的对应关系。

1 系统描述

常用非线性滤波算法如扩展卡尔曼滤波(EKF)、无迹卡尔曼滤波(UKF)以及容积卡尔曼滤波(CKF)，适用于非线性系统的最优状态估计。在卫星轨道估计中常采用非线性滤波算法进行实时轨道估计，是一种在时域上对轨道状态矢量进行最优估计的方法，需要对系统建立状态方程与量测方程，利用状态方程进行时间更新同时对利用量测方程进行量测更新，最终给出最优估计。对系统进一步进行可观测性分析时，根据控制理论的要求，需要求出非线性系统的状态转移方程与观测方程。

1.1 状态方程

根据卫星在轨道运行的动力学规律，建立卫星的轨道动力学模型。地球是一个形状、质量分布皆非均匀的扁球体，对航天器的引力需考虑地球非球形引力，大气阻力和太阳光压、潮汐等其他摄动力。非球形引力项中最大的是J2项，数量级远远大于其摄动力，在计算中仅考虑该项。在地心惯性坐标系中，单颗卫星的运动方程可以表示为引力、地心距离与卫星运动矢量的微分方程。在J2000坐标系下建立状态方程，如下所示：

(1)

在计算机仿真分析中，连续系统无法进行时域上的状态估计，于是将式(1)的连续方程通过数值计算方法离散化，常用方法有欧拉法、龙格库塔离散方法等，再对系统进行泰勒一阶展开，变为(2)所示的离散非线性时变系统：

Xk+1=Φk+1,kXk

(2)

式中,

(3)

Xk+1为k时刻的航天器的状态量，Xk为k时刻的状态量，Φk+1,k为k时刻到k+1时刻的状态转移矩阵。

1.2 量测方程

量测过程是通过地面站对航天器进行测量，常用的测量元有测距、测速及测角。在实际应用中，相较于包含完整测元的体制，仅采用单测角或单测速的系统轨道估计精度较低，更具有研究意义，因此本文采用单种测速值的测量系统，对其进行可观测性研究。

每个测站给出的径向速度量测值通过航天器的状态矢量表示出，由此构成量测方程组。多个测站对同一航天器进行观测获取多个测速值，在单站的情况下进行扩维运算，构成多站仅测速系统的量测方程：

(4)

利用泰勒一阶展开对量测方程进行处理，对系统状态求偏导，获得系统观测矩阵[11]：

(5)

量测方程建立在站心坐标系中，状态方程建立在J2000坐标系中，进行滤波运算时需要提前将数据进行坐标系转换，在同一坐标系下进行轨道估计。

2 可观测性分析

可观性及可控性是线性系统的一种定性分析方法，指系统内的状态是否可以由输出和输入进行分析，对于系统的运行有重要研究意义。在控制理论研究中，线性系统可观测性有明确的判断准则，如格拉姆矩阵判据。任何线性系统的定理推广到非线性系统时，复杂程度大大增加。因为非线性系统中某一时刻的状态无法由初始状态推导，非线性系统可观测性分析没有统一的判断准则。

2.1 非线性系统的可观测性分析

非线性系统可通过泰勒展开等近似为线性化系统，再采用线性系统的判断准则来判断，但该种方法对牺牲了对系统的准确描述使计算变得简明。而直接对非线性系统的可观测性判定在工程运用中，通常采用由Lee和Dunn提出的李函数定理[2]：

对非线性系统:

(6)

y(t)=h(x(t),t)

(7)

如果在t∈[t0,t1]上，对凸集S∈Rn上的任意x0，x1∈S，h(x(t;x),t)=h(x(t;x1),t)， 都有x1=x0，则系统在凸集S上是完全可观测的。线性系统在S集上的可观测性可以通过对初始状态的可观性分析推导出来，与其状态向量的变化无关。非线性系统的可观性不能由初始状态推导，因此考虑由状态方程f(·)和量测方程h(·)进行推导。

定义：

(8)

(9)

式中，H(t)、Φ(t,t0)为h(·)和f(·)的雅各比矩阵。若M(x0)是正定的，则系统在凸集S上是完全可观测的。

将该定理应用在非线性离散时变系统中，定义可观测性矩阵Γ为：

(10)

式中，N表示观测次数，n表示状态维数。

如果rankΓ(k0,k0+N-1)=n，即在N次观测中，可观测矩阵Γ的秩等于n，则系统在S上是完全可观测的。将式(3)、(5)带入(10)可以得到可观测矩阵，对其求秩可以对系统的可观测性进行判断，通过求取特征值或进行奇异值分解，都可以对系统的可观测度进行判断。

系统中的状态噪声与量测噪声对系统的可观测性也存在影响，对系统本身结构的可观测性描述可以不考虑噪声的影响。考虑噪声时的系统可观测性及可观测度判断有其他的判断方法。

2.2 可观测度分析方法

可观测性定性的反应出系统的实际运行能力，为了更清楚反应系统可观测性的大小，引入可观测度的概念，对系统的可观测性给出一种定量的分析方法。

矩阵的条件数定义为最大奇异值和最小奇异值的比值，可观测矩阵的条件数反映可观测矩阵的病态程度，可以理解为可观测矩阵各行列之间的相关性，相关性强表明量测值中获取的相似信息量多，有效信息量少；相关性弱则表明量测值中获取的有效信息量多。这样能直观的反映出算法的稳定性、收敛性及收敛速度，将条件数的大小作为可观测度的判断依据。

对可观测矩阵Γ进行式(11)所示的奇异值分解：

Γ=U∑V

(11)

式中，U、V为奇异值分解后的酉矩阵，∑为对角矩阵，对角线上的元素为可观测矩阵Γ的奇异值。条件数cond(Γ)为最大奇异值和最小奇异值的商。

(12)

在实际应用及仿真中发现，式(11)、(12)直接计算出的条件数并不适用于多站单测量体制的滤波定轨算法。传统的条件数计算中将不同测站获取的同种体制的量测值直接带入可观测矩阵Γ计算，所得到的可观测度仅仅是量测信息本身之间数值大小上的相关性，当两个处于不同位置的测站所获取的量测值在仅仅在数值上相似时，条件数会猛然增大，可知其并不能正确体现不同测站所包含有效信息的多少，定轨精度的好坏与所得到的条件数大小并不存在对应关系。

(13)

(14)

3 仿真分析

为了研究可观测度是否能作为滤波定轨算法精度的参考依据，采用基于容积卡尔曼滤波的滤波定轨算法，在卫星工具包(Satellite Tool Kit，STK)中建立场景模型，卫星模型采用太阳同步轨道卫星(Sun-synchronous satellite)，轨道估计数据在轨道估计算法中采用高精度轨道预报(high precision orbit propagation, HPOP)，在条件数计算中采用仅考虑J2项的轨道数据，更利于状态方程的雅各比矩阵的计算同时仿真出的曲线平滑，最终利用MATLAB进行算法验证。

对模型中随机放置的两个测站和随机放置的三个测站在不同几何布设的情况下计算优化方法所得的条件数，同时进行定轨算法仿真以得出不同测站组合下的定轨误差结果，滤波步长为1 s，条件数计算时长为100 s，定轨算法计算时长为170 s，在仿真模型中F1、F2、F3、F4分别代表四个随机测站，四个测站的经纬度如表1所示。

表1 四个测站的经纬度

仿真初始轨道值在模型真实轨道值的基础上位置量加入1 km的误差，速度量加上100 m/s的误差，参考真实初始轨道值为：

初始协方差矩阵为：

状态噪声协方差矩阵为：

选择通过均方根误差判断定轨精度，根据数据的量级分为速度均方根误差和位置均方根误差。

Eposition=

(15)

Evelocity=

(16)

图1 双测站的条件数变化曲线

图2 双测站轨道估计位置均方根误差

图3 双测站轨道估计速度均方根误差

图4 三测站的条件数变化曲线

图5 三测站轨道估计位置均方根误差

图6 三测站轨道估计速度均方根误差

由图1和图4可以看出，由于系统是时变的，所以系统的条件数也是时变的，结合系统的条件数与图2、图3以及图4、图5的定轨精度误差结果可以看出，条件数越小，定轨精度越高，算法的收敛速度更快。为了更直观的分析条件数与定轨精度之间的关系，如表2所示，将传统方法所求的各系统的条件数均值及改进方法所求的条件数均值与均方误差值进行对比，比较不同方法求得的条件数大小与系统定轨精度的关系。

表2 定轨精度与两种可观测度方法所求得条件数的比较

图7 两种方法所得条件数与估计精度的变化趋势

由表2和图7可以看出，传统方法所求得的条件数的大小无法准确反映系统的定轨精度，并不适用于条件数大，量测值之间存在不同数值意义的弱观测系统。改进方法将不同测站坐标作为向量进行投影处理，代入可观测矩阵所求条件数，仿真结果数值上大于传统方法，但是对系统可观测性的判断准确度高。因此，实验结果表明改进方法所求得的条件数大小能正确反映出定轨精度的变化，当条件数较大时，系统的定轨精度较差；条件数较小时，系统的定轨精度较高。

4 结论

研究了基于多站测速的滤波定轨算法的可观测性及可观测度判断方法。传统可观测度判定方法中条件数的计算方法其结论对于多站单测量值的弱观测系统不适用，改进方法将不同测站获取的量测值看作矢量进行坐标系转换，将测站坐标值带入计算，给出一种适应多站组合测量弱观测条件下的基于条件数的可观测度判定方法。数值仿真结果验证了改进方法相对于传统方法，有更好的估计效果，更能准确作为系统可观测度判定的方法。

参考文献：

[1] Kalman R E. A New Approach to Linear Filtering and Prediction Problems.[J]. J.basic Eng.trans.ASME, 1960; 82D(1):35-45.

[2] Lee T, Dunn K, Chang C. On observability and unbiased estimation of nonlinear systems[J]. System Modeling and Optimization, 1982: 258-266.

[3] 周卫东, 蔡佳楠, 孙 龙,等. GPS/SINS超紧组合导航系统可观测性分析[J]. 北京航空航天大学学报, 2013, 39(9):1157-1162.

[4] 黄翔宇, 崔平远, 崔祜涛. 深空自主导航系统的可观性分析[J]. 宇航学报, 2006, 27(3):332-337.

[5] 宁晓琳, 房建成. 航天器自主天文导航系统的可观测性及可观测度分析[J]. 北京航空航天大学学报, 2005, 31(6):673-677.

[6] 刘 准, 陈 哲. 条件数在系统可观测性分析中的应用研究[J]. 系统仿真学报, 2004, 16(7): 1552-1555.

[7] 曲 毅, 刘 忠, 薛 锋. 基于TDOA的双站自适应滤波算法及可观测性研究[J]. 系统仿真学报, 2007, 19(22):5222-5225.

[8] 常晓华, 崔平远, 崔祜涛. 一种深空自主导航系统可观测性分析方法[J]. 哈尔滨工业大学学报, 2010, 42(11):1681-1685.

[9] 黄 普, 钱 山, 李恒年. 空间目标地基雷达跟踪可观测性分析[J]. 宇航动力学学报, 2015, 5(3).

[10] 孙仲康,周一宇,何黎星.单/多基地有源无源定位技术[M].北京：国防工业出版社,1996.

[11] 李恒年. 航天测控最优估计方法[M]. 北京：国防工业出版社, 2015.