西安电子科技大学，陕西 西安

收稿日期：2021年11月19日；录用日期：2021年12月16日；发布日期：2021年12月23日

摘要

对于空间高速运动目标比如导弹来说，通过天基网络对目标进行光学定位，可以很好地改善地基观测系统定位精度不高的问题。在本文的研究中先对单星的可观测性能进行了分析，然后基于天基观测系统中的双星目标几何定位得到目标的位置信息，再使用无迹卡尔曼滤波算法(UKF)对目标进行定位跟踪。从仿真结果分析可知，无迹卡尔曼滤波算法在对超高速运动的目标进行定位跟踪时有更高的稳定性和优越性，可使得对目标的定位误差达到100 m/1000 Km以下。

关键词

天基光学，空间目标，单星观测，双星定位，无迹卡尔曼滤波

Research on Satellite Optical Positioning Based on Space-Based Observation Network

Hao Fu, Xiaoning Fu

Xidian University, Xi’an Shaanxi

Received: Nov. 19^th, 2021; accepted: Dec. 16^th, 2021; published: Dec. 23^rd, 2021

ABSTRACT

For space high-speed moving targets such as missiles, optical positioning of targets through space-based network can well improve the problem of low positioning accuracy of ground-based observation system. In this paper, the observability of a single satellite is analyzed, and then the position information of the target is obtained based on the geometric positioning of the double satellite target in the space-based observation system, and then the unscented Kalman filter algorithm (UKF) is used to locate and track the target. The simulation results show that the unscented Kalman filter algorithm has higher stability and superiority in positioning and tracking the ultra-high-speed moving target, and the positioning error of the target can reach less than 100 m/1000 km.

Keywords:Space Based Optics, Space Target, Single Satellite Observation, Double Satellite Positioning, Unscented Kalman Filter

Copyright © 2021 by author(s) and Hans Publishers Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY 4.0).

http://creativecommons.org/licenses/by/4.0/

1. 引言

随着各国向太空发射的卫星越来越多，利用天基网络观测系统对目标进行检测已成为当前发展的趋势，且相对于地基观测系统观测时间短，定位精度不高，天基网络光学观测系统受到的条件限制较少，探测范围大、定位精度高，可以用来实现特定目标的跟踪定位。

目前美国大力发展的空间目标监视系统(SBSS)对运动目标具有成像效果好、定轨精度高、监视能力强、处理速度快等优点 [1]；还有由天基红外预警系统(SBIRS)中低轨部分发展而来的空间跟踪与监视系统(STSS)和精确跟踪与监视系统(PTSS)目前已经成功实现部分在轨运行，其余部分也在不断完善 [2]。

在天基网络观测系统中，最基础也是最重要的部分就是分析系统的定位精度。由于光学传感器只能获得光学卫星与观测目标之间的相对角度信息，比如方位角和俯仰角，故使用单颗光学卫星来观测无法获取目标在实际空间内所处的具体位置 [3]。目前天基跟踪定位主要是通过结合多星光学成像信息的天基组网的方式，在星载光电平台上装上空间光学相机，通过控制光学相机的姿势对目标进行持续的跟踪定位 [4]。与直接将相机安装到卫星上相比，使用星载平台来观测目标可以减少对卫星本体的控制，降低目标观测的难度。

本文通过对理论分析以及仿真实验相结合，先后对单星的观测性能，双星的定位过程和定位误差进行了相应的分析研究，以及通过无迹卡尔曼滤波算法对目标进行跟踪定位，达到提高精度，降低误差的目的。

2. 单星观测模型

天基观测系统目标定位示意图如图1所示，是一种由多颗卫星在单个或者多个轨道平面实现特点需求的卫星网络系统，将观测平台装载到每颗卫星上，如图2所示。当目标出现在光学成像视场内时，可通过图像处理算法获取成像仪于目标之间的方位角和俯仰角，还可以通过闭环控制系统来驱动星载平台的转动跟踪和监视目标运动。

图1. 卫星网络目标定位示意图

图2. 星载观测平台

2.1. 卫星观测性能

对卫星的可观测性能分析是对卫星无源定位技术研究的基础，如果系统不可观测时，对目标定位的确定值就无从谈起。因此需要对卫星测角无源定位系统的可观测性进行分析，为接下来研究无源定位算法提供理论依据 [5]。卫星的观测模型如图3所示。

图3. 单星观测示意图

1) 可观测性理论

假定有无噪声的连续系统如下：

$\begin{array}{l}X\left(t\right)=F\left(t\right)X\left(t\right)+G\left(t\right)U\left(t\right)\\ Z\left(t\right)=H\left(t\right)X\left(t\right)\end{array}$ (1)

上式中， $X\left(t\right)$ 为状态向量； $F\left(t\right)$ 为状态矩阵； $G\left(t\right)$ 、 $U\left(t\right)$ 分别为控制矩阵和控制矢量； $Z\left(t\right)$ 、 $H\left(t\right)$ 分别为观测矢量和观测矩阵。

求解上述微分方程，得

$X\left(t\right)=\Phi \left(t,t_0\right)X\left(t_0\right)+C\left(t\right)$ (2)

其中， $\Phi \left(t,t_0\right)$ 为系统状态转移矩阵，且有如下关系：

$\begin{array}{l}\frac{\partial \Phi \left(t,t_0\right)}{\partial t}=F\left(t\right)\Phi \left(t,t_0\right)\\ \Phi \left(t,t_0\right)=I_n\\ C\left(t\right)={\displaystyle {\int }_{t_0}^t\Phi }\left(t,\tau \right)G\left(\tau \right)U\left(\tau \right)\text{d}\tau \end{array}$ (3)

上式中， $I_n$ 为n阶单位矩阵。

有如下定理：

$M={\displaystyle {\int }_{t_0}^{t_1}{\Phi }^{\text{T}}}\left(t,t_0\right)H^{\text{T}}\left(t\right)H\left(t\right)\Phi \left(t,t_0\right)\text{d}t$ (4)

则系统在 $t\in \left[t_0,t_1\right]$ 的可观测性充要条件为观测矩阵M正定。

考虑到在计算或判断观测矩阵M正定性会遇到困难，则上述定理有如下等价形式：

系统在 $t\in \left[t_0,t_1\right]$ 的可观测性充要条件为：对 $\forall t\in \left[t_0,t_1\right]$，有 $H\left(t\right)\Phi \left(t,t_0\right)Y\equiv 0$，则 $Y=0$。

2) 观测卫星系统建模

根据二体运动规律以及万有引力定律可假设 $r\left(t\right)={\left[\begin{array}{lll}x\left(t\right)\hfill & y\left(t\right)\hfill & z\left(t\right)\hfill \end{array}\right]}^{\text{T}}$ 为卫星在地心惯性坐标系下的位置矢量，则卫星的运动微分方程为：

$a=-\mu \frac{r}{R^3}$ (5)

式中，a为卫星的加速度矢量即 $a={\left[\begin{array}{lll}\ddot{x}\left(t\right)\hfill & \ddot{y}\left(t\right)\hfill & \ddot{z}\left(t\right)\hfill \end{array}\right]}^{\text{T}}$ ；R为卫星与地心的质心距； $\mu$ 为开普勒常数。

取观测卫星位置矢量和速度矢量作为状态变量 [6]，即 $X\left(t\right)={\left[\begin{array}{ll}{r}^{\text{T}}\hfill & {\stackrel{\dot{}}{r}}^{\text{T}}\hfill \end{array}\right]}^{\text{T}}$，其中：r为位置矢量 $r\left(t\right)={\left[\begin{array}{lll}x\left(t\right)\hfill & y\left(t\right)\hfill & z\left(t\right)\hfill \end{array}\right]}^{\text{T}}$，$\stackrel{\dot{}}{r}$ 为速度矢量 $r\left(t\right)={\left[\begin{array}{lll}\stackrel{\dot{}}{x}\left(t\right)\hfill & \stackrel{\dot{}}{y}\left(t\right)\hfill & \stackrel{\dot{}}{z}\left(t\right)\hfill \end{array}\right]}^{\text{T}}$，则卫星运动状态方程为：

$\begin{array}{l}r=b\left(\cos E\left(t\right)-e\right)P+b\sqrt{1-e^2}\sin E\left(t\right)Q\\ \stackrel{\dot{}}{r}=\frac{\sqrt{\mu b}}{r}\left(-\sin E\left(t\right)P+\sqrt{1-e^2}\cos E\left(t\right)Q\right)\end{array}$ (6)

式中，b为卫星轨道的长半轴；e为轨道偏心率；i为轨道倾角； $\Omega$ 为升交点赤经； $\omega$ 为近地点幅角； $E\left(t\right)$ 为当前t时刻卫星的偏近角 [7]；P、Q分别为

$\begin{array}{l}P=\left[\begin{array}{c}\cos \Omega \cos \omega -\sin \Omega \sin \omega \cos i\\ \sin \Omega \cos \omega +\cos \Omega \sin \omega \cos i\\ \sin \omega \sin i\end{array}\right]\\ Q=\left[\begin{array}{c}-\cos \Omega \sin \omega -\sin \Omega \cos \omega \cos i\\ \sin \Omega \sin \omega +\cos \Omega \cos \omega \cos i\\ \cos \omega \sin i\end{array}\right]\end{array}$ (7)

3) 观测方程

由于方位角 $\beta$ 和俯仰角 $\epsilon$ 为卫星对目标定位主要的观测信息，则有

$\begin{array}{l}\tan \left(\beta \left(t\right)\right)=x^{\prime }/y^{\prime }\\ \tan \left(\epsilon \left(t\right)\right)=z^{\prime }/\sqrt{{x^{\prime }}^2+{y^{\prime }}^2}\end{array}$ (8)

其中， $\left(x^{\prime },y^{\prime },z^{\prime }\right)$ 表示所测目标在以卫星为坐标原点的观测坐标系中的位置。

4) 可观测性分析

结合轨道根数相互关系以及二体运动相关理论，整理并化简得

$X\left(t\right)=\left[\begin{array}{ll}{I}_{3}f\hfill & I_{3}g\hfill \\ I_3\stackrel{\dot{}}{f}\hfill & I_3\stackrel{\dot{}}{g}\hfill \end{array}\right]\left[\begin{array}{l}r\left(t_0\right)\hfill \\ \stackrel{\dot{}}{r}\left(t_0\right)\hfill \end{array}\right]=\left[\begin{array}{ll}{I}_{3}f\hfill & I_{3}g\hfill \\ I_{3}f\hfill & I_{3}g\hfill \end{array}\right]X\left(t_0\right)$ (9)

$f,g$ 分别为

$\begin{array}{l}f=\frac{\cos \left(E\left(t\right)-E\left(t_0\right)\right)-e\cos E\left(t_0\right)}{1-e\cos E\left(t_0\right)}\\ g=\frac{b^3}{\sqrt{\mu }}\sin \left(E\left(t\right)-E\left(t_0\right)\right)-\frac{b^3}{\sqrt{\mu }}\left[e\left(\sin E\left(t\right)-\sin E\left(t_0\right)\right)\right]\end{array}$ (10)

上式中， $b,g$ 为固定值。

由式(8)可得

$\begin{array}{l}0=x^{\prime }\left(t\right)\cos \beta \left(t\right)-y^{\prime }\left(t\right)\sin \beta \left(t\right)\\ 0=-x^{\prime }\left(t\right)\sin \beta \left(t\right)\sin \epsilon \left(t\right)+z^{\prime }\left(t\right)\cos \epsilon \left(t\right)-y^{\prime }\left(t\right)\cos \beta \left(t\right)\sin \epsilon \left(t\right)\end{array}$ (11)

将 $X\left(t\right)={\left[\begin{array}{ll}{r}^{\text{T}}\hfill & {\stackrel{\dot{}}{r}}^{\text{T}}\hfill \end{array}\right]}^{\text{T}}$ 作为状态量，式(11)可写为矩阵形式，

$\left[\begin{array}{c}0\\ 0\end{array}\right]=\left[\begin{array}{cc}\cos \beta \left(t\right)& -\sin \beta \left(t\right)\sin \epsilon \left(t\right)\\ -\sin \beta \left(t\right)& -\cos \beta \left(t\right)\sin \epsilon \left(t\right)\\ 0& \cos \epsilon \left(t\right)\\ 0& 0\\ 0& 0\\ 0& 0\end{array}\right]\left[\begin{array}{l}{x}^{\prime }\left(t\right)\hfill \\ y^{\prime }\left(t\right)\hfill \\ z^{\prime }\left(t\right)\hfill \\ {\stackrel{\dot{}}{x}}^{\prime }\left(t\right)\hfill \\ {\stackrel{\dot{}}{y}}^{\prime }\left(t\right)\hfill \\ {\stackrel{\dot{}}{z}}^{\prime }\left(t\right)\hfill \end{array}\right]$ (12)

将上式整理得

$Z\left(t\right)={\left[\begin{array}{cc}\cos \beta \left(t\right)& -\sin \beta \left(t\right)\sin \epsilon \left(t\right)\\ -\sin \beta \left(t\right)& -\cos \beta \left(t\right)\sin \epsilon \left(t\right)\\ 0& \cos \epsilon \left(t\right)\\ 0& 0\\ 0& 0\\ 0& 0\end{array}\right]}^{\text{T}}X\left(t\right)$ (13)

根据上式(9)和(13)可得

$\Phi \left(t,t_0\right)=\left[\begin{array}{ll}{I}_{3}f\hfill & I_{3}g\hfill \\ I_3\stackrel{\dot{}}{f}\hfill & I_3\stackrel{\dot{}}{g}\hfill \end{array}\right]$ (14)

$H\left(t\right)=\left[\begin{array}{cc}\cos \beta \left(t\right)& -\sin \beta \left(t\right)\sin \epsilon \left(t\right)\\ -\sin \beta \left(t\right)& -\cos \beta (t)\sin \epsilon \left(t\right)\\ 0& \cos \epsilon \left(t\right)\\ 0& 0\\ 0& 0\\ 0& 0\end{array}\right]$ (15)

根据可观测性理论，若有矢量 $Y={\left[y_1,y_2,y_3,y_4,y_5,y_6\right]}^{\text{T}}$ 使得 $H\left(t\right)\Phi \left(t,t_0\right)Y=0$，则有下式：

$\begin{array}{l}\left(f{y}_1+g{y}_4\right)\cos \beta \left(t\right)-\left(f{y}_2+g{y}_3\right)\sin \beta \left(t\right)=0\\ -\left(f{y}_1+g{y}_4\right)\sin \beta \left(t\right)\sin \epsilon \left(t\right)-\left(f{y}_2+g{y}_5\right)\cos \beta \left(t\right)\sin \epsilon \left(t\right)+\left(f{y}_3+g{y}_6\right)\cos \epsilon \left(t\right)=0\end{array}$ (16)

因正弦函数和余弦函数是两个线性独立的函数，则由式(15)得：

$\{\begin{array}{l}f{y}_1+g{y}_4=0\\ f{y}_2+g{y}_5=0\\ f{y}_3+g{y}_6=0\end{array}$ (17)

下面对上式中形如 $a_{1}f+a_{2}g=0$ 的式子进行讨论：

1) 若 $a_1\ne 0$ 且 $a_2\ne 0$，则有 $f/g=-a_1/a_2$ 为常数，但由式(10)可得， $f,g$ 都为时间t的函数，且 $f,g$ 之间不存在固定的倍数关系，因而 $f/g$ 不等于常数，故假设不成立。

2) 若 $a_1,a_2$ 只有一个为0，不妨设 $a_1$ 为0， $a_2\ne 0$，则有 $a_{2}g=0$，根据式(10)，g为t的函数且不恒等于0，故可知假设不成立，同理亦可证明 $a_1\ne 0,a_2=0$ 不成立。

则只有当 $a_1=0$ 且 $a_2=0$ 时， $a_{1}f+a_{2}g=0$ 成立，则由式(6)得， $Y=0$。

特别需要说明的是，虽然在上述定性分析过程中忽略了观测误差以及系统噪声对卫星可观测性分析的影响，这是因为系统噪声和观测误差的存在只对目标定位跟踪的精度有影响，而与系统的可观测性分析无关。综上可得，卫星无源测角系统是完全可观测的，上述定理论证成立。

2.2. 常用坐标系及其转换

2.2.1. 常用坐标系

本文主要介绍了四种坐标系，分别为：大地坐标系 $\left\{B,L,H\right\}$，地固坐标系 $\left\{X_F,Y_F,Z_F\right\}$，地惯坐标系 $\left\{X_I,Y_I,Z_I\right\}$ 和轨道坐标系 $\left\{X_O,Y_O,Z_O\right\}$ [7]。具体定义如表1所示。

表1. 常见坐标系及其定义

2.2.2. 坐标系转换

1) 大地坐标系与地固坐标系的转换公式如下：

$\left[\begin{array}{c}{X}_F\\ Y_F\\ Z_F\end{array}\right]=\left[\begin{array}{c}\left(N+H\right)\cos B\cos L\\ \left(N+H\right)\cos B\sin L\\ \left[N\left(1-e_0^2\right)+H\right]\sin B\end{array}\right]$ (18)

其中， $N=a/\sqrt{1-e_0^2{\sin }^{2}B}$ 为卯酉圈半径； $e_0=\sqrt{a^2-b^2}/a$ 为地球第一偏心率；a为地球椭球长半周；b为地球椭球短半轴。

2) 地心惯性坐标系和轨道坐标系的转换：

由地心惯性系依次绕Z轴旋转Ω，绕X轴旋转 $-\left(\pi /2-i\right)$，绕Y轴旋转 $-\left(\pi /2+\upsilon \right)$，然后进行平移，即可得到：

$\left[\begin{array}{c}{X}_O\\ Y_O\\ Z_O\end{array}\right]=\left[\begin{array}{c}{X}_{OI}\\ Y_{OI}\\ Z_{OI}\end{array}\right]+R_Z\left(\Omega \right)R_X\left(-\left(\pi /2-i\right)\right)R_Y\left(-\left(\pi /2+\upsilon \right)\right)\left[\begin{array}{c}{X}_I\\ Y_I\\ Z_I\end{array}\right]$ (19)

其中 $\left(X_{OI},Y_{OI},Z_{OI}\right)$ 为地心惯性系坐标原点的位置；Ω为升交点赤经；i为轨道倾角； $\upsilon$ 为纬度幅角，其中 $\upsilon =\omega +f$，$\omega$ 为近地点角距，f为真近角点。

3. 双星定位模型

在天基观测系统中，目前常用的无源定位算法一般有三种：基于奇异值分解的最小二乘算法，总体最小二乘算法以及几何定位算法 [8]。对于由两颗光学卫星组成的观测体系，使用几何定位算法不仅定位精度较高，而且计算速度最快，可以获得一个较优解。

3.1. 双星定位基本原理

在使用单颗卫星对空间目标进行观测时，只能得到目标的方位信息即方位角e和俯仰角a，如图4所示，而无法直接得到目标的具体位置信息。

图4. 目标观测角度示意图

其中观测角度信息定义如下：

方位角e：卫星的视线矢量在面 $X_O{Y}_O{O}_1$ 内的投影与坐标轴 $X_O{O}_1$ 的夹角，由 $X_O{O}_1$ 指向 $O_1{T}^{\prime }$ 为正，其取值范围是 $\left[0,\pi \right]$。

俯仰角a：卫星的视线矢量在面 $X_O{Y}_O{O}_1$ 内的投影 $O_1{T}^{\prime }$ 与视线 $O_{1}T$ 的夹角，由 $O_1{T}^{\prime }$ 指向 $O_{1}T$ 为正，其取值范围是 $\left[-\frac{\pi }{2},\frac{\pi }{2}\right]$。

通过方位角e和俯仰角a可以计算出在轨道坐标系下所测目标的视线单位矢量如下：

$\{\begin{array}{l}{v}_{c1}={\left[\cos e_1\cos a_1,\cos a_1\sin e_1,\sin a_1\right]}^{\text{T}}\\ v_{c2}={\left[\cos a_2\cos e{}_2,\cos a_2\sin e{}_2,\sin a_2\right]}^{\text{T}}\end{array}$ (20)

双星定位过程如图5所示，由于观测目标在空间中存在一定的体积，导致卫星观测系统中两颗卫星的观测视线可能不能在空间中交于一点，而是产生一定的偏离 [4]。因此在计算目标实际位置时，目标在三维空间中的具体位置可以通过两颗卫星观测视线的公垂线中点来表示，如图6所示，这样可以尽量避免目标自身体积带来的误差。双星几何定位具体流程如图7所示。

图5. 目标双星定位模型

图6. 目标点求解示意图

图7. 双星几何定位流程图

图5中，C代表观测卫星1在地惯坐标系下的空间位置 $r_{s1}$，D表示观测卫星2在地惯坐标系下的空间位置 $r_{s2}$，CA与DB表示卫星1和卫星2的观测视线，观测视线的单位矢量即为观测矢量 $v_1,v_2$，AB则为视线CA和视线DB的公垂线。

假设 $C\left(x_1,y_1,z_1\right),C\left(x_2,y_2,z_2\right),v_1\left(a_1,b_1,c_1\right),v_2\left(a_2,b_2,c_2\right)$，可得：

${\theta }_1=\arccos \left[\frac{a_1\left(x_2-x_1\right)+b_1\left(y_2-y_1\right)+c_1\left(z_2-z_1\right)}{\sqrt{{\left(x_2-x_1\right)}^2+{\left(y_2-y_1\right)}^2+{\left(z_2-z_1\right)}^2}}\right]$ (21)

${\theta }_2=\arccos \left[\frac{a_2\left(x_1-x_2\right)+b_2\left(y_1-y_2\right)+c_2\left(z_1-z_2\right)}{\sqrt{{\left(x_2-x_1\right)}^2+{\left(y_2-y_1\right)}^2+{\left(z_2-z_1\right)}^2}}\right]$ (22)

$\theta =\arccos \left(a_1{a}_2+b_1{b}_2+c_1{c}_2\right)$ (23)

根据空间异面直线定理：

$\left|CA\right|=\frac{cos{\theta }_1+cos{\theta }_{2}cos\theta }{{\sin }^2\theta }\ast \sqrt{{\left(x_2-x_1\right)}^2+{\left(y_2-y_1\right)}^2+{\left(z_2-z_1\right)}^2}$ (24)

$\left|DB\right|=\frac{cos{\theta }_2+cos{\theta }_{1}cos\theta }{{\sin }^2\theta }\ast \sqrt{{\left(x_2-x_1\right)}^2+{\left(y_2-y_1\right)}^2+{\left(z_2-z_1\right)}^2}$ (25)

则空间目标在地惯坐标系下的位置为：

$r=\frac{1}{2}\left[r_{s1}+r_{s2}+\left|CA\right|\ast v_1+\left|DB\right|\ast v_2\right]$ (26)

3.2. 双星定位仿真

STK (Satellite Tool Kit)是由美国Analytical Graphics公司开发的一款在航天领域处于领先地位的商业分析软件 [9]。使用STK可以很方便的建立Walker星座来搭建天基观测网络来模拟目标定位。在该实验中，设置Walker星座的轨道数量为3条，每个轨道上的卫星数量为8颗，相邻轨道卫星相位为1，将星座高度设置为1000 KM。

在天基观测网络目标定位过程，选取第一条轨道上的第4颗卫星和第二条轨道上的第2颗卫星来分析双星定位，选择观测时间段为16:20:00~16:35:00，即双星观测到目标在空间中的飞行时间为900 s。具体过程如图8所示。

图8. STK双星定位示意图

根据天基光学传感器的观测视线测量误差 [10] 和观测卫星自身的位置误差 [11]，设置观测卫星自身的位置误差为16 m，速度误差为0.5 m/s，观测角度误差为(0.5 ± 0.05) e⁻³ rad。对与定位目标的均方根误差可以使用蒙特卡罗法来计算，其计算公式如下：

$RMSE=\sqrt{\frac{1}{N}{\displaystyle \underset{i=1}{\overset{N}{\sum }}\left[{\left(x_r-x\right)}^2+{\left(y_r-y\right)}^2+{\left(z_r-z\right)}^2\right]}}$ (27)

其中N为蒙特卡罗打靶仿真次数，根据 [12] 可知，在本次仿真中设为180次即可满足相应精度； $\left(x_r,y_r,z_r\right)$ 是定位算法计算得到的目标位置； $\left(x,y,z\right)$ 是目标的真实位置；经过Matlab仿真，可以得到目标的定位精度误差，如图9所示。

从图中可以看出，在本次仿真中，双星对于目标的定位精度平均处于280 m/1000 Km附近，考虑到卫星到目标的平均观测距离为3000~4000 Km左右，其实际定位误差在1000 m左右，x轴方向定位精度较高，误差在100 m/1000 Km左右，y轴方向次之，定位误差大约为180 m/1000 Km，而z轴方向的定位误差大约为220 m/1000 Km。

图9. 目标定位的均方根误差

3.3. 定位误差分析

目标的定位精度受很多因素的影响，观测系统的最终定位误差来源种类繁多。对于观测系统来讲，系统自身的静态指向误差、动态稳定误差、观测系统与卫星之间的对准误差等等；从卫星的角度来讲，卫星的自身位置和卫星观测目标时的观测角度误差都将影响定位精度；而且目标定位跟踪是根据卫星星载平台上的光学相机来队目标进行观测，其对目标进行观测时也存在诸多误差，比如卫星姿态角引起的误差，相机指向误差等。且通过以上仿真发现，定位误差也与卫星所处的位置有关，当目标与两颗卫星的相对位置大致处于一条直线上时，形成的观测模型将会影响对目标的定位计算，降低目标的定位精度；同时，当目标分别与两卫星的距离接近时，定位精度相较于距离较大时有所提高。

4. 无迹卡尔曼滤波

4.1. 卡尔曼滤波

卡尔曼滤波(KF)是一种基于目标上一时刻的状态估计值和观测系统的观测数据，通过系统的状态方程递推估计目标下一时刻的状态的线性最小方差估计方法。卡尔曼滤波可以通过递推计算，计算起来也非常简单方便，其对目标的状态也可进行试试更新和处理，非常便于计算机编程实现，因而其在工程上应用尤其广泛。在本文的研究中，目标的信息是通过天基观测网络中的双星定位得到的，且目标的运动是一种连续的状态，故可以通过该时刻之前的测量信息来递推该时刻目标的运动状态。

系统的状态方程和量测方程经离散化可表示如下：

$\{\begin{array}{l}{X}_k={\Phi }_{k/k-1}{X}_{k-1}+G_{k-1}{W}_{k-1}\\ Z_k=H_k{X}_k+V_k\end{array}$ (28)

其中， $X_k$ 为目标时刻的状态变量； ${\Phi }_{k/k-1}$ 为 $k-1$ 到k时刻的系统一步转移矩阵； $G_{k-1}$ 为系统噪声转移矩阵； $W_{k-1}$ 为 $k-1$ 时刻的系统噪声； $W_{k-1}$ 本质是均值为零；方差为 $Q_k$ 的高斯白噪声。 $Z_k$ 为k时刻的观测向量； $H_k$ 为k时刻的量测噪声矩阵； $V_k$ 为k时刻的量测噪声； $V_k$ 也是零均值的高斯白噪声，其方差为 $R_k$。

KF算法计算步骤如下：

步骤1：目标状态预测：

${\hat{X}}_{k/k-1}={\Phi }_{k/k-1}{\hat{X}}_{k-1/k-1}$ (29)

步骤2：计算卡尔曼增益：

$K_k=P_{k/k-1}{H}_k^{\text{T}}{\left(H_k{P}_{k/k-1}{H}_k^{\text{T}}+R_k\right)}^{-1}$ (30)

步骤3：状态估计方程：

${\hat{X}}_{k/k}={\hat{X}}_{k/k-1}+K_k\left(Z_k-H_k{\hat{X}}_{k/k-1}\right)$ (31)

步骤4：预测协方差：

$P_{k/k-1}={\Phi }_{k/k-1}{P}_{k-1}{\Phi }_{k/k-1}^{\text{T}}+G_{k-1}{Q}_{k-1}{G}_{k-1}^{\text{T}}$ (32)

步骤5：估计协方差更新：

$P_{k/k}=\left(I-K_k{H}_k\right)P_{k/k-1}$ (33)

式中， ${\hat{X}}_{k/k}$ 是目标k时刻的状态估计值； ${\hat{X}}_{k/k-1}$ 是由 $k-1$ 时刻的估计值； ${\hat{X}}_{k-1}$ 对下一时刻 ${\widehat{X}}_k$ 递推得到的预测值； $P_{k/k-1}$ 为预测值 ${\hat{X}}_{k/k-1}$ 的协方差阵。

4.2. 无迹卡尔曼滤波

基于线性回归的思想，无迹卡尔曼滤波(UKF)是通过UT变换构造Sigma点然后结合KF的基本框架来更新目标的状态和协方差 [13]。无迹卡尔曼滤波不用计算雅可比矩阵且无需忽略高阶项，所以精度会更高。并且UKF是在原状态分布中按照一定的规则采样，使这些点的均值以及协方差和原状态分布中的点相同，再之后将这些点的集合带入到非线性状态方程中，求得经过变换后的均值和方差。具体步骤如图10所示。

其非线性系统可以表示为：

$\{\begin{array}{l}{X}_{k+1}=F\left(X_k\right)+w_k\\ Z_k=H\left(X_k\right)+v_k\end{array}$ (34)

其中， $X_{k+1}$ 为状态向量； $Z_k$ 为观测量； $w_k,v_k$ 分别为系统过程噪声以及观测噪声； $w_k,v_k$ 都是均值为零的高斯白噪声；其对角方差阵对应为 $Q_k,R_k$。

UKF的基本步骤如下：

步骤1：初始化：

状态向量初值： ${\hat{x}}_0=E\left[x_0\right]$

协方差初始估计值： $P_0=E\left[\left(x_0-{\hat{x}}_0\right){\left(x_0-{\hat{x}}_0\right)}^{\text{T}}\right]$

步骤2：采样点计算：

使用UT变换后，2n + 1个对称的Sigma采样点如下：

$\{\begin{array}{l}{\chi }_{0,k-1}={\hat{x}}_{k-1}\hfill \\ {\chi }_{i,k-1}={\hat{x}}_{k-1}+{\left(\sqrt{\left(n+\lambda \right)P_{k-1}}\right)}_i,i=1~n\hfill \\ {\chi }_{i,k-1}={\hat{x}}_{k-1}-{\left(\sqrt{\left(n+\lambda \right)P_{k-1}}\right)}_{i-n},i=n+1~2n\hfill \end{array}$ (35)

其中 ${\hat{x}}_{k-1}$ 和 $P_{k-1}$ 分别为 $k-\text{1}$ 时刻的状态估计值和协方差矩阵； $\lambda ={\alpha }^2\left(n+k\right)-n$，$\alpha$ 决定采样点的分布，其取值范围是： $0\le \alpha \le 1$ ；n为状态向量的维数；k一般大于0来保证方差阵的半正定性。

步骤3：状态方程玉协方差方程更新：

$\{\begin{array}{l}{\chi }_{k/k-1}^i={\Phi }_{k/k-1}{\chi }_{k-1/k-1}^i\\ {\hat{X}}_{k/k-1}={\displaystyle \underset{i=0}{\overset{2n}{\sum }}{W}_i^m}{\chi }_{k/k-1}^i\\ P_{k/k-1}={\displaystyle \underset{i=0}{\overset{2n}{\sum }}{W}_i^m}\left[{\chi }_{k/k-1}^i-{\hat{X}}_{k/k-1}\right]{\left[{\chi }_{k/k-1}^i-{\hat{X}}_{k/k-1}\right]}^{\text{T}}+Q_k\\ Z_{k/k-1}^i=h\left({\chi }_{k/k-1}^i\right)\\ {\hat{Z}}_{k/k-1}={\displaystyle \underset{i=0}{\overset{2n}{\sum }}{W}_i^m}{Z}_{k/k-1}^i\end{array}$ (36)

计算权值为：

$W_0^m=\frac{\lambda }{\lambda +n},W_i^m=\frac{1}{2\left(\lambda +n\right)},i=1,2,\cdots ,2n$ (37)

步骤4：量测更新方程：

系统测量预测协方差为：

$\{\begin{array}{l}{P}_{{\hat{z}}_k{\hat{z}}_k}={\displaystyle \underset{i=0}{\overset{2n}{\sum }}{W}_i^c}\left[Z_{k/k-1}^i-{\hat{Z}}_{k/k-1}\right]{\left[Z_{k/k-1}^i-{\hat{Z}}_{k/k-1}\right]}^{\text{T}}+R_k\\ P_{x_{k,}{z}_k}={\displaystyle \underset{i=0}{\overset{2n}{\sum }}{W}_i^c}\left[{\chi }_{k/k-1}^i-{\hat{X}}_{k/k-1}\right]{\left[Z_{k/k-1}^i-{\hat{Z}}_{k/k-1}\right]}^{\text{T}}\end{array}$ (38)

权值：

$W_0^c=\frac{\lambda }{\lambda +n}+1-{\alpha }^2+\beta ,W_i^c=\frac{1}{2\left(\lambda +n\right)},i=1,2,\cdots ,2n,\beta \ge 0$ (39)

滤波增益：

$K_k=P_{x_k{z}_k}{P}_{{\hat{z}}_k{\hat{z}}_k}^{-1}$ (40)

状态估计值与下一时刻协方差矩阵：

$\{\begin{array}{l}{\hat{x}}_k={\hat{x}}_{k/k-1}+K_k\left(Z_k-{\hat{Z}}_{k/k-1}\right)\hfill \\ P_k=P_{k/k-1}-K_k{P}_{{\hat{z}}_k{\hat{z}}_k}{K}_k^{\text{T}}\hfill \end{array}$ (41)

图10. 无迹卡尔曼滤波算法示意图

4.3. 仿真实验与跟踪定位

4.3.1. 仿真实验

为了控制误差来源，依旧设置观测卫星自身的位置误差为16 m，导弹目标的速度误差为0.5 m/s，光学卫星的观测角度误差为(0.5 ± 0.05) e⁻³ rad。通过STK仿真软件设置导弹发射点为美国关岛军事基地(E13, N114deg)，落点为朝鲜平壤(E39, N125)。选取卫星观测网络中的第二条轨道上的低2颗卫星和第三条轨道上的第8颗卫星进行定位跟踪，选取双星共同观测到导弹的时间段：17:44:00~17:47:00，共计3分钟.如图11所示。

图11. 双星导弹跟踪定位

4.3.2. 仿真结果分析

分别使用线性卡尔曼滤波和无迹卡尔曼滤波对目标定位跟踪，在50次蒙特卡洛仿真下的对目标定位跟踪的位置误差和定位精度如图12所示：

图12. 滤波算法跟踪定位轨迹图

图13. 滤波算法跟踪定位误差

从仿真图13中可以看出，无迹卡尔曼滤波算法明显的减小了跟踪的定位的误差，并能够保持较好的跟踪性能。其x轴方向的平均定位误差在450 m左右；y轴方向的平均定位误差在250 m附近；z轴方向的平均定位误差大约为50 m；总体位置的平均定位误差大约在400~500 m范围内；考虑到卫星到被测目标的距离大致为5000 Km，故通过无迹卡尔曼滤波算法对目标进行定位可以使误差保持在100 m/1000 Km以下。

5. 结论

本文首先介绍了单星的可观测性能以及常用的坐标系及其变换，然后通过STK仿真软件搭建了基于Walker卫星星座的天基观测网络，基于星座中的双星目标几何定位得到目标的位置信息，再使用无迹卡尔曼滤波算法(UKF)对目标进行定位跟踪，可以使得对目标的定位误差达到100 m/1000 Km以下。并且从仿真结果分析可知，无迹卡尔曼滤波算法在对超高速运动的目标进行定位跟踪时有更高的稳定性和优越性。

基于天基观测系统的空间目标定位检测还需进一步探索，在使用双星定位过程中可以发现，随着目标与卫星的相对运动，刚开始选取的两颗卫星在一段时间之后并不是最理想的观测卫星，会随着时间的推移出现观测位置和角度更好的卫星。如果能够实现对定位目标的跨卫星观测域跟踪定位，及时更换观测效果更好的卫星来定位跟踪，将会进一步地降低定位误差，提高定位精度。

文章引用

付 豪,付小宁. 基于天基观测网络的卫星光学定位研究
Research on Satellite Optical Positioning Based on Space-Based Observation Network[J]. 计算机科学与应用, 2021, 11(12): 2953-2970. https://doi.org/10.12677/CSA.2021.1112299

参考文献