0 引言
随着各国航天技术迅速发展，地球周围空间产生了包括航天器、火箭末级和空间碎片等在内的大量人造物体，对空间目标态势的感知和掌握，是人类不断开拓视野、不断进步的过程[1-2]。空间态势感知包括对空间目标的探测、跟踪、识别，对空间事件的评估和预报等[3],空间态势感知是人类开展航天活动的基础。空间态势感知按照信息源的部署位置划分，主要分为地基和天基两类，与地基相比，天基不会受到天气、环境和地缘条件的限制，可长时间对空间目标进行观测。同时，诸多实践证明，天基观测平台对于监视空间目标很有优势，能够很好地弥补地基观测设备的不足[4]。因此，天基态势感知受到了广泛关注。
空间目标跟踪是天基空间态势感知的关键技术之一。目前，空间目标跟踪主要分为测距跟踪和测角跟踪[3]。测距跟踪主要是采用有源、主动式的工作方式，利用微波雷达、激光雷达等传感器获取到的目标测距信息进行跟踪定轨。测距跟踪具有较高的跟踪精度，但传感器功耗较大，且有源主动工作方式使其隐蔽性较差。测角跟踪是空间目标中更常见的跟踪方式，利用光电传感器获取到的测角信息进行定轨，采用无源、被动式的工作方式。此外，为了进一步提升空间目标跟踪性能，获取更加精确的状态估计结果，可以使用多颗观测卫星对同一个目标进行跟踪[5],并使用非线性的滤波算法[6]，进一步提升空间目标跟踪能力。
空间目标跟踪系统的可观测性是能通过滤波估计获得较好跟踪精度的先决条件。因此，研究跟踪系统可观测性有助于选择相应的策略提升系统的可观测性，进一步改善系统的跟踪性能。本文针对传统空间目标跟踪精度不高的问题，提出一种基于可观测性分析的高精度空间目标跟踪方法。首先，基于Fisher信息矩阵法定义了一种新型可观性指标，并以空间目标跟踪为背景，针对采用天基雷达，天基红外或地基雷达所得到的不同量测模型系统分别给出了可观性指标值，并选出一种可观性指标值最大的空间目标量测模型。接着，基于无迹卡尔曼滤波和选择的量测模型，实现空间目标跟踪。本文基于可观性指标筛选出观测值较高的空间目标跟踪测量方案，进而提高空间目标跟踪精度。
1 空间目标动力学模型与测量模型
本文采用地心固连坐标系OXYZ，大地坐标系和传感器本体坐标系，分别如下所示。
1) 地心固连坐标系OXYZ坐标原点O位于地心，OX轴为格林尼治子午线与赤道面正向交线的方向，OZ轴垂直于赤道面并指向地球的北极，OY轴由OZ轴与OX轴根据右手螺旋定律确定。
2) 大地坐标系包括大地纬度Be、经度Le以及高度He。大地纬度Be为赤道平面与地球椭球面法线的夹角，向北为正，大地经度Le为在赤道面内从格林尼治向东度量的角度，大地高度He为从地球椭球表面量起的法向距离。
3) 传感器本体坐标系OSXSYSZS的坐标原点OS为传感器的中心位置，OSXS轴为地心O与传感器的中心OS连线的延长线，OSZS轴为在OSXS轴与地球自转轴组成的平面内并垂直于OSXS轴，OSYS由OSZS轴与OSXS轴根据右手螺旋定律确定。
在OXYZ坐标系下，假设空间目标的位置r和速度v分别为
(1)
式中：x，y，z表示在OXYZ坐标系3个方向的位置；vx，vy，vz表示在OXYZ坐标系3个方向的速度矢量。
根据文献[7]，可得在OXYZ坐标系下，考虑J2项摄动的影响，空间目标的轨道动力学方程为
(2)
(3)
式中：μe为地球引力常数，且μe=3.986 006 4×1014 m3/s2；ωe为非惯性坐标系的角速度，根据文献[7]中的WGS-84模型，ωe=7.292 115×10-5 rad/s；J2为参考地球模型的J2项摄动常数，且J2=1.082 626 836×10-3；r为空间目标到地心的距离；Re为地球赤道半径且Re=6.378 14×106 m；ce为
(4)
选取状态变量X=[rT, vT]T，且假设太阳光压力摄动等未建模摄动的影响为零均值高斯白噪声，根据式(2)和式(3)，得到空间目标的连续状态方程为
(5)
式中：w为未建模系统误差，是零均值高斯白噪声，且E(wwT)=Q。
根据式(5)，空间目标的连续状态方程在变量Xk-1处采用四阶Runge Kutta方法展开，得到如下空间目标的非线性离散状态方程：
Xk=f(Xk-1)+wk-1
(6)
式中：w为离散化未建模系统误差，且
假设卫星i (i=1, 2) 在OXYZ坐标系下的位置矢量为ri(i=1, 2)，卫星的纬度和经度分别为φ和θ，则空间目标在OSXSYSZS坐标系下的位置矢量rSi为
rSi=MOXYZ→OSXSYSZS(ri-rS) (i=1, 2)
(7)
式中，MOXYZ→OSXSYSZS为OXYZ坐标系到OSXSYSZS坐标系的转移矩阵，且
(8)
图1表示传感器本体坐标系。如图 1所示，rSi =(xSi, ySi, zSi)表示空间目标在OSXSYSZS坐标系下的位置矢量，α和β分别为方位角和俯仰角，且
(9)
图1 传感器本体坐标系
本文是基于天基雷达传感器，在双星观测情况下，采用测角测量和测距测量，根据式(9)，考虑到测角误差和测距误差，测量方程为
Z=g(X)+v=
[α1, β1, rS1, α2, β2, rS2]T+v
(10)
式中：rSi(i=1,2)表示位置矢量rSi(i=1,2)的模；v为测量误差，是零均值高斯白噪声，且E(vvT)=R；Z为某一时刻空间目标在两颗卫星的传感器本体坐标系下测角和测距的测量值。
根据式(10)，得到非线性离散观测状态方程：
Zk=g(Xk)+vk
(11)
式中,vk为离散化观测误差，且
2 无迹卡尔曼滤波
根据式(6)与式(11)，可得空间目标的双星天基雷达测量的离散化数学模型为
(12)
由于当系统非线性程度较高时，传统卡尔曼滤波及其改进型扩展卡尔曼滤波算法精度会降低甚至会引发滤波发散。因此本文采用Julier和Uhlman[8]提出的无迹卡尔曼滤波方法。相比传统卡尔曼滤波及其改进型[9]，UKF法不需要解析求解雅克比矩阵，且易工程实现。而且，由于UKF法采用确定性采样的方法，其在处理非线性估计时具有较高精度[10]。
无迹卡尔曼滤波过程如下：
1) 初始条件设置
(13)
式中：α*为刻度因子常数；β*为第二刻度因子常数；κ为第三刻度因子常数，κ在高斯分布的情况下最优值为2；L为状态的总维数；
和
分别为计算均值和方差的权重。
初始均值
和初始方差P0表示为
(14)
式中：X0和
表示状态变量的初始值及其期望值；P0表示状态变量的初始协方差矩阵。
2) 采样点确定
式中：
表示协方差矩阵Pk-1的Cholesky分解，即
(16)
3) 计算时间更新
(17)
(18)
(19)
4) 计算观测更新
(20)
(21)
(22)
(23)
5) 计算滤波增益和滤波更新
(24)
(25)
(26)
式中：
为状态估计的后验均值；Pk为状态估计的后验协方差矩阵。
上述为UKF滤波算法的详细步骤，该算法针对非线性高斯系统，采用无迹变换对称采样的方式得到的状态估计值能够精确到Tayor展开级数3阶[11]。针对任意非线性系统，状态估计值逼近精度到Tayor展开级数2阶[11]。本文测量方程的非线性程度较高，因此采用UKF算法进行空间目标的状态估计。
3 Fisher信息矩阵法
空间目标跟踪系统可观性是能通过滤波估计获得较好的跟踪精度的先决条件。因此，研究跟踪系统可观性有助于选择相应的策略提升系统的可观性，进一步改善系统导航性能。非线性系统可观性计算矩阵判断方法比线性系统更为复杂，本文通过线性化假设近似计算非线性系统可观性矩阵的Fisher信息矩阵法[12]。
首先对系统状态方程(6)和观测方程(11)进行线性化，得到线性化的系统模型：
(27)
式中，Φk和Hk分别表示状态转移矩阵和观测转移矩阵，且
(28)
(29)
线性化后的信息的递推更新过程为
(30)
由于Fisher信息矩阵的逆代表着无偏估计的所能达到的误差下线，因此可以利用该矩阵来衡量跟踪系统的状态估计性能。通常，可以使用Fisher信息矩阵的行列式或迹来度量系统可观性，行列式或迹的值越大，空间目标跟踪系统可观性更好。本文基于Fisher信息矩阵的迹，定义Ω表示可观测性因子，选取所有离散时刻中Fisher信息矩阵的迹的最大值来定义可观测性因子Ω的值，并利用该参数描述空间目标跟踪系统的可观测性，且
Ω=logζ1{max[det(F)]-ζ2}
(31)
式中,det(F) 表示Fisher信息矩阵的行列式，ζ1， ζ2 表示可观测性因子参数，且ζ1=10，ζ2=9.996 999 71×1047。
4 仿真与分析
对近地轨道飞行卫星编队进行仿真，考虑到地球J2项摄动的影响，以地球WGS-84模型为参考对象，采用如式(6)所示的动力学模型，比较三种不同的测量模型：①双星均携带天基红外传感器；②一星携带天基雷达传感器一星携带天基红外传感器；③双星均携带天基雷达传感器。基于双星编队对空间目标的跟踪数据，采用UKF算法实现空间目标在OXYZ坐标中精确定位。本文所采用的参考卫星，伴随卫星和空间目标的初始状态矢量如表 1所示。
表1 初始状态矢量
类别x/my /mz /m参考卫星3.23×1053.51×1066.54×106伴随卫星5.18×1053.21×1066.79×106空间目标2.30×10-92.42×1066.65×106类别vx/(m·s-1)vy/(m·s-1)vz /(m·s-1)参考卫星-6.78×103-3.20×1031.30×103伴随卫星-6.35×103-3.11×1031.56×103空间目标-5.76×103-4.54×1031.64×103
本文初始状态估计误差ΔX，初始协方差矩阵P0，初始Fisher信息矩阵F0和系统状态方程噪声协方差矩阵Q如下所示：
(32)
(33)
F0=inv(P0)=
(34)
(35)
测量模型1为双星均携带天基红外传感器，则测量模型1为
Z1=g1(X)+v1=[α1, β1, α2, β2]T+v1
(36)
式中：v1为测量误差，是零均值高斯白噪声，且
为某一时刻空间目标在两颗卫星的传感器本体坐标系下测角的测量值。
测量模型1的测量方程噪声协方差矩阵R1为
(37)
测量模型2为一星携带天基雷达传感器一星携带天基红外传感器，则测量模型2为
Z2=g2(X)+v2=
[α1, β1, α2, β2, rS2]T+v2
(38)
式中：v2为测量误差，是零均值高斯白噪声，且
为某一时刻空间目标在两颗卫星的传感器本体坐标系下测角和测距的测量值。
测量模型2的测量方程噪声协方差矩阵R2为
(39)
测量模型3为双星均携带天基雷达传感器，测量模型3如式(10)所示。测量模型3的测量方程噪声协方差矩阵R为
(40)
基于上述参数和式(31)，得出三种不同的测量模型情况下可观测性因子值如表2所示。测量模型3的可观测性因子值最高，测量模型1的可观测性因子值最低。因此双星编队上携带天基雷达传感器时，可观测性因子值最好，此时该空间目标跟踪系统的跟踪性能最好。图2表示三种不同的测量模型情况下空间目标跟踪精度。在图2中，横坐标为时间，纵坐标为位置估计标准差以10为底的对数。如图2所示，测量模型3的空间目标跟踪精度最高，测量模型1的空间目标跟踪精度最低，因此，双星编队上携带天基雷达传感器时，空间目标跟踪精度最高，与采用可观测性因子的分析结果一致，验证了可观测性因子的有效性。
表2 三种不同的测量模型情况下可观测性因子值
类别模型1模型2模型3可观测因子值38.5338.7038.87
图2 三种不同的测量模型情况下空间目标跟踪精度
5 结束语
针对传统空间目标跟踪精度不高的问题，提出一种可观测性分析的高精度空间目标跟踪方法。首先，基于Fisher信息矩阵法定义了一种新型可观性指标，并以空间目标跟踪为背景，针对采用天基雷达，天基红外或地基雷达所得到的不同量测模型系统分别给出了可观性指标值，并选出一种可观性指标值最大的空间目标量测模型。接着，基于无迹卡尔曼滤波和选择的量测模型，实现高精度的空间目标跟踪。相对于仅天基红外，天基红外与天基雷达相结合两种方案，仅天基雷达的空间目标跟踪测量方案具有较大的可观测性指标值，进而具有较高目标跟踪精度。最后，数值仿真验证了该方法的有效性。
[1] 王兴龙, 蔡亚星, 陈士明, 等. 多源信息融合在空间态势感知领域的应用与发展[J]. 航天返回与遥感,2021,42(1):11-20.
[2] 黄晓斌，张燕，肖锐，等. 空间目标的雷达定轨实时识别问题研究[J]. 雷达科学与技术, 2021,19(1):63-68.
[3] 杜扬钦. 基于双目视觉卫星编队的空间态势感知关键技术研究[D]. 浙江: 浙江大学,2020.
[4] 胡云鹏, 黎克波, 陈磊. 面向空间态势感知的天基可见光空间目标自主跟踪方法[J]. 中国科学：技术科学,2021,51(4):424-434.
[5] 姚培东, 李星秀, 吴盘龙, 等. 双星编队对空间非合作目标的仅测角跟踪方法[J]. 导航与控制,2020,19(2):27-35.
[6] 何山, 吴盘龙, 恽鹏, 等. 一种面向临近空间高超声速再入滑翔目标跟踪算法[J]. 宇航学报,2020,41(5):553-559.
[7] 秦政. 基于低轨红外星座的中段目标跟踪技术研究[D]. 长沙：国防科技大学，2020.
[8] JULIER S J,UHLMANN J K. Unscented Filtering and Nonlinear Estimation[J]. Proceedings of the IEEE,2004,92(3):401-422.
[9] 王祎, 杏建军, 于洋, 等. 编队卫星相对位置测量信息自主确定月球引力场的方法[J]. 计算机工程与科学,2016,38(1):125-130.
[10] 欧阳威. 火星探测器自主导航方案设计与分析[D]. 长沙：国防科技大学，2017.
[11] SIMON D. Optimal State Estimation:Kalman, H infini-ty, and Nonlinear Approaches [M]. USA:John Wiley & Sons,2006.
[12] SIMANDL M, KRALOVEC J, TICHAVSKY P. Filtering, Predictive, and Smoothing Cramer-Rao Bounds for Discrete-Time Nonlinear Dynamic Systems[J]. Automatica,2001,37(11):1703-1716.
