多星载光学传感器系统误差极大似然配准算法

0 引 言

弹道导弹的助推段红外辐射特征显著,不难被下视的红外预警卫星发现,但处于中段飞行的弹道导弹目标温度偏低,现代光学传感器(主要是红外传感器)只能在深空冷背景下对其进行探测,因此,利用低轨卫星作为载体,采用临边扫描方式探测中段飞行弹头的方法是目前的发展趋势,如美国在研的空间跟踪监视系统(space tracking and surveillance system, STSS)[1-2]。一方面，由于光学传感器采用被动探测方式,只能获取目标的视线,即相对于传感器的方位角和俯仰角量测,没有直接的距离量测,难以精确估计目标的状态参数[3-4]。另一方面,由于低轨星载光学传感器受视场限制,且目标与卫星间的相对位置关系快速变化,需要多传感器间的组网协同以实现对目标的快速定位与连续跟踪[5-8]。实际工程应用中,光学传感器获取的视线量测不可避免地包含系统误差,因此为提高多传感器对目标的融合精度,必须首先进行多传感器的系统误差配准[9-13]。

国内外学者对误差配准技术虽然已经展开了深入研究,但主要集中在地基和机载传感器的误差配准方面,而关于星载光学传感器配准的研究相对较少。文献[14-16]以STSS系统为背景,基于恒星观测与其已知位置之间的差值,构建偏差滤波器和非线性状态估计器,实现了传感器指向误差和目标状态的实时分离估计。文献[17-18]在文献[16]的基础上,将影响传感器视线的误差因素分别等效为姿态失配角和视线指向偏移角等偏差角,通过建立偏差角的状态转移模型和参考点的量测模型,采用非线性滤波器实现了对指向误差的在线估计与补偿。上述配准方法属于在线估计算法,能够实现系统误差的实时估计,但前提是已知参考点的先验信息,且需要增加额外的滤波器以实现目标状态估计。文献[19]又提出了基于迭代最小二乘(iterative least squares, ILS)的离线批处理误差估计算法,但仅研究了两部传感器情况下的配准问题。对于多传感器(传感器数量大于2)情况,文献[20-21]针对地基多传感器,提出极大似然配准(maximum likehood registration, MLR)算法,其不受传感器数量与体制的约束,利用多部传感器的观测信息,精确地实现了系统误差和目标轨迹的联合估计。文献[22]又将MLR算法应用到了多空基平台传感器的误差配准场景。考虑到机动平台姿态对传感器观测的影响,文献[23]提出平台姿态时变条件下的MLR算法,实现了传感器量测和姿态角系统误差的同时估计。文献[24]在上述文献基础上,推导了无源传感器观测量对目标状态的雅克比矩阵,进一步实现了缺失距离条件下多无源传感器的MLR,但未考虑到平台定位和姿态角测量等随机误差对配准性能产生的影响。

本文以STSS系统为背景,对低轨星座协同探测弹道目标情况下的多光学传感器配准问题进行了研究。由于低轨星座在轨道空间中高速运动,因此在考虑卫星轨道和姿态动态变化的条件下,建立了星载光学传感器的观测模型。通过一阶Taylor近似法对观测过程中涉及到的测量误差进行分析,推导出目标状态的误差协方差与卫星轨道定向、姿态角测量和传感器测量等随机误差的关系式,并基于视线交叉得到了观测在状态空间中的近似投影。在此基础上,将MLR算法扩展到多星载光学传感器的配准,提出多星载光学传感器系统误差MLR算法。该算法实现了多传感器系统误差与目标轨迹的融合估计,并通过引入各类测量误差的先验信息对目标状态的误差协方差进行修正,以提高算法的配准性能。

1 系统模型

定义地心惯性(earth-centered inertial,ECI)坐标系为公共参考系,即J2000平赤道平春分点坐标系,并且定义卫星的轨道坐标系与本体坐标系,以描述观测过程中低轨卫星的轨道与姿态变化[25]。

1.1 轨道运动模型

对于低轨卫星,其在轨道空间中的运动特性近似服从二体运动规律。假设在ECI坐标系下的卫星位置矢量为r=[x,y,z]T,加速度矢量为并假定地球为WGS-84标准的椭球体,保留引力二阶带谐项,则其轨道运动满足如下微分方程[26]：

(1)

式中,为卫星到地心的距离;μ为地心引力常数;aJ2为二阶带谐项引起的瞬时摄动,表达式为

式中,J2=1.082 63×10-3;Re=6 378 137 m为地球赤道半径。

在初值给定的情况下,通过四阶Runge-Kutta方法对式(1)积分,可得任意时刻的卫星运动状态X,即位置r、速度由卫星的运动状态可求得与之相对应的3个角度参数Φ=[Ω,i,u]T,即升交点赤经、轨道倾角和纬度辐角,其共同确定了轨道在空间中的定向,表达式分别为

Ω=arctan(-Wx/Wy)

(2)

(3)

(4)

式中,W为单位面积速度矢量,可由位置和速度矢量计算得到。

1.2 传感器观测模型

对于星载光学传感器,观测模型描述了目标空间位置与视线角度的转换关系。设在ECI坐标系下N个传感器对同一目标进行连续同步观测,并已知低轨卫星的运动状态。在k时刻,卫星n的状态矢量为目标位置矢量为rk,则ECI坐标系下的相对位置矢量为

Δrk,n=rk-rk,n

(5)

定义卫星轨道坐标系的原点位于卫星质心,zo、xo轴在轨道平面内分别指向地心和运动方向,yo轴与xo、zo轴构成右手系。设卫星轨道的3个角度参数量测为Φn,卫星轨道测量噪声为其为零均值白色高斯过程,协方差为根据ECI坐标系到轨道坐标系的转换关系可得目标在卫星轨道坐标系下的位置矢量为

(6)

式中,为ECI坐标系到轨道平面坐标系的坐标变换;为常值矩阵,分别表示为

式中,R1(·)、R3(·)为绕x轴和z轴的旋转矩阵。

定义固连于卫星的本体坐标系,原点位于卫星质心,xs轴沿卫星纵轴指向前,zs轴在纵对称面内指向下,ys轴垂直于纵对称面指向右,构成右手系。设卫星n的姿态角量测为Θn=[ϑn,γn,ξn]T,姿态角测量噪声为依次绕坐标轴旋转相应的姿态角,可将转换到本体坐标系下

(7)

式中,ϑn、ξn、γn为描述卫星三轴姿态的偏航、俯仰和横滚姿态角;为零均值白色高斯噪声,协方差为为绕y轴的旋转矩阵。

设非线性观测函数为hn(·),则传感器n包含观测误差的观测量为

(8)

式中,zk,n=[θk,n,φk,n]T为观测矢量；θk,n、φk,n分别为视线方位角和俯仰角。

非线性观测函数的表达式为

联合式(6)～式(8),可得包含系统误差的观测方程为

(9)

式中,为传感器系统误差矢量;为测量噪声矢量,其为零均值白色高斯过程,协方差为且不同传感器间的相互独立。

2 星载光学传感器测量误差分析

2.1 测量误差传递

在对系统误差和目标轨迹进行估计前,首先要将传感器的观测矢量投影到目标状态空间,由式(9)可得映射关系

(10)

式中,为观测均值矢量。利用式(10),可将观测矢量zk,n投影到目标状态空间。

(11)

根据式(11),目标状态误差主要是在卫星轨道定向误差、姿态角测量误差和传感器测量误差的综合作用下,经非线性量测转换而产生。在已知系统误差bn的情况下,目标状态xk,n可近似为一个均值为rk的高斯过程,利用一阶Taylor近似对非线性量测转换进行线性化处理,可近似求得其误差协方差矩阵。

为使表达简明,省略传感器下标n。考虑Φ上存在测量误差δΦ,对旋转矩阵在均值处进行一阶Taylor级数展开,可得

(12)

式中,为关于Φ的梯度矩阵;等式右侧第二项为对轨道误差δΦ的全微分,可表示为

同理,在均值处的一阶Taylor展开形式为

(13)

其右侧第二项为关于姿态误差δΘ的全微分,即

ϑ

对于观测函数的反函数h-1(·),其在观测均值处的一阶Taylor近似表达式为

(14)

式中,h-1(·)关于观测误差δz的全微分为

zkh-1(zk-b)·δzk=

将上述3次变换的一阶Taylor展开式代入式(11),忽略位置误差和二阶及二阶以上的高阶项,则目标状态误差为

JΦ·δΦ+JΘ·δΘ+Jzk·δz

(15)

式中,JΦ、JΘ和Jzk为各类测量误差到目标状态的误差传递矩阵,即雅克比矩阵,式(15)的矢量紧凑形式可表示为

(16)

假设各量测量之间具有独立性,则各测量误差互不相关,因此目标状态的误差协方差为

(17)

如果式(12)～式(14)中的Taylor展开式保留到二阶项或二阶以上高阶项,则可得到二阶或更高阶的目标状态误差

(18)

文献[27-28]对不同阶数的Taylor近似精度进行了理论与仿真分析,结果表明包含二阶项的Taylor近似与一阶近似相比具有更高精度,而二阶以上Taylor展开的近似精度并没有明显的提高,因此一般不采用高阶Taylor近似。虽然利用二阶Taylor展开对非线性量测转换进行近似具有更高精度,但其计算难度高且计算量很大,一般情况下通常采用一阶Taylor近似[29]。

2.2 状态空间中的投影近似

由于光学传感器缺少距离维信息,导致观测在目标状态空间中无法获得唯一描述,即式(11)中的xk,n不唯一,进而无法获取目标状态误差协方差。解决该问题的一种可行方法是通过旋转变换获得目标视线量测在ECI坐标系下的投影,利用投影矢量交叉对目标位置进行近似估计,再结合卫星位置和目标位置估计结果,得到距离信息,并对观测矢量进行扩维处理。

假设某时刻传感器n的观测矢量zn=[θn,φn]T,则视线量测为Ln=[cos θncos φn, sin θncos φn,sin φn]T,结合卫星的轨道量测和姿态角量测,可得视线量测在ECI坐标系下的投影为

(19)

根据目标与卫星在ECI坐标系下的相对位置关系,可得

(20)

式中,‖r-rn‖为目标到卫星n的距离。将式(20)的第3个表达式分别代入前2个表达式,整理可得单颗卫星的定位方程为

Anr=an

(21)

式中

由于式(21)中方程个数小于未知量个数,因此单星条件下的定位方程为欠定方程。考虑将N颗卫星的定位方程组合在一起,令则定位方程为

Ar=a

(22)

由于式(22)是超定的,可利用最小二乘法对目标位置进行近似,得到目标位置估计[30]为

(23)

结合已知的卫星位置,可得距离项估计为通过对原始观测进行扩维处理,得到将其代入式(11)后,即可获得状态空间中的投影近似。

3 多星载光学传感器系统误差MLR算法

3.1 非随机误差的极大似然估计

考虑用观测集Z={zk;k=1,2,…,K}来估计非随机系统误差为当前观测集对应的时刻,zk为k时刻系统的观测矢量集

(24)

其矢量函数形式为

(25)

系统误差的极大似然估计量就是使似然函数p(Z|b)极大化的b值,即

(26)

由于目标位置rk未知,在b和rk均未知的情况下,仅能得到似然函数p(Z|X,b),其中目标轨迹X={rk;k=1,2,…,K}。通过联合极大化关于X和b的似然函数p(Z|X,b),依据噪声序列在时域内相互独立,可得X和b的极大似然估计量为

(27)

式中,p(zk|rk,b)为观测矢量zk的概率密度函数。

3.2 目标轨迹估计

使得概率密度函数p(zk|rk,b)在目标状态空间中极大化的rk就是目标轨迹的极大似然估计,考虑到噪声序列在观测平台间具有相互独立性,可将p(zk|rk,b)表示为紧凑形式,即

p(zk,1,zk,2,…,zk,N|rk,b)=

(28)

其在目标状态空间中的近似形式可表示为

p(zk|rk,b)≈

(29)

式中,为xk,n误差协方差的逆。

令式(29)的极大化条件等价于使其右侧的指数函数中第一项等于零,即因此k时刻目标轨迹的极大似然估计量其表达式为

(30)

3.3 传感器系统误差估计

假设系统误差的迭代项为将迭代项代入式(11)可得

(31)

考虑bn上存在微小扰动,对式(11)在处进行一阶Taylor级数展开,则目标状态与系统误差的增量之间具有线性关系,即

(32)

式中,Jbn为关于bn在处的雅克比矩阵,可表示为

(33)

令k时刻目标的状态矢量集联合式(32)有

(34)

式中

Qk=blkdiag(-Jb1,-Jb2,…,-JbN)

其中,blkdiag(·)是由多个矩阵构成的分块对角形矩阵。

为获得误差矢量的极大似然估计,需将目标轨迹的极大似然估计量代入到似然函数p(Z|X,b),然后对进行极大化处理。将式(30)代入到式(29),则指数函数中第一项为零,进一步化简可得

(35)

式中

式中,右侧第二项是由N2个3×3的子矩阵构成的方阵,i、j分别代表子矩阵所在的行列位置。

结合独立性假设,似然函数可展开为

(36)

代入式(33)后,进一步化简为式(28)的形式。

(37)

式中,C与均为常量。同样,使式(37)极大化的条件等价于使其右侧的指数函数中第一项等于零,即因此系统误差的极大似然估计量即

(38)

3.4 算法执行步骤

本文算法主要包括误差估计和轨迹估计两个步骤。误差估计步骤中,利用当前系统误差的估计值b0对原始观测去偏化,根据视线交叉定位对缺失的距离项进行估计,并计算出观测矢量zk在目标状态空间的投影及误差协方差矩阵,使极大化后可得到一个新的误差估计轨迹估计步骤中,根据估计值的收敛性,计算出目标状态空间中对应的配准后的传感器观测,利用极大化似然函数p(zk|rk,b),可获得对目标的轨迹估计。图1总结了MLR算法的执行关键步骤。

图1 算法执行流程
Fig.1 Flow chart of the algorithm

4 仿真实验与性能分析

4.1 场景设置及评价指标

假设仿真场景由Walker低轨星座和1个弹道目标构成,以历元时间7 Jul 2018 01:50:00作为场景初始时刻[31]。Walker星座的构型参数为28/4/2/77.8°/1 600 km;目标的发射时刻为7 Jul 2018 02:07:00,初始位置为(140° E, 42°N),射向(84.2° E, 41.5°N)附近,射程约为4 541.8 km,具体场景如图2所示。根据卫星工具箱(satellite tool kit,STK)分析可知,星座中有3颗卫星可同时观测到目标,其编号分别为1、7、16,对应的轨道根数如表1所示。

图2 仿真场景示意
Fig.2 Sketch of simulation scenario

表1 历元时刻卫星的轨道根数

Table 1 Orbital elements of satellites at the epoch

假设星载光学传感器对目标同步观测,观测间隔为1 s,其观测性能参数如下:各传感器的系统误差分别为b1=[-3.5°,-1.5°]T,b2=[3°,1°]T,b3=[-2.5°,-1°]T;测量误差均为σθn=0.2°,σφn=0.15°。卫星姿态参数采用欧拉角描述法,相应的姿态角测量误差为σϑn= σγn=σξn=0.1°。各低轨卫星的轨道测量误差分别为σΩn=σin=0.01°,σun=0.02°。

本文主要以系统误差估计的误差均值和误差标准差作为性能评价指标,用于校验所提算法对星载光学传感器的配准性能。同时,为评估算法理论上能达到的最佳估计精度,以文献[20]定义的克拉美罗下界 (Cramer-Rao lower bound,CRLB)为指标,该指标是一个保守值,但对于数量不大的蒙特卡罗仿真,可以很好地预估算法最优估计性能,其定义式为

(39)

式中,xk为目标的真实值rk而非估计值。

4.2 实验结果与性能分析

实验首先通过进行单次配准仿真,对所提算法的正确性进行验证,计算得到的3颗低轨卫星的系统误差估计结果如图3所示。由图3可知,对于任意一颗卫星,算法经过3～5次迭代运算后,对系统误差的估计值均能稳定收敛于预设值附近,估计精度基本达到95 %以上,表明算法在多颗卫星观测的环境下具有良好的误差估计性能。然而,在迭代次数到达一定数量之后,对误差估计精度的提高能力有限。在此基础上,图4给出了200次蒙特卡罗仿真得到的目标轨迹融合结果,可以看出算法对目标轨迹的估计接近于真实轨迹,表明算法在距离缺失条件下仍能较好地实现对目标轨迹的融合估计。

图3 系统误差估计结果
Fig.3 Results of systemic error estimation

通过100次蒙特卡罗仿真,统计得到在迭代次数为5、采样点数为200～400环境下算法对系统误差估计的误差均值和误差标准差,分别如图5和图6所示。由图5可知,在采样点数较少时,算法对不同卫星的误差估计仍存在一定程度的偏差;当采样点数大于310时,误差均值维持在0.1°以下,并且随着采样点数的增加,均趋近于0°,表明得到的系统误差估计是渐进无偏估计量。由图6可知,对于各个传感器的系统误差估计,估计精度均稳定收敛于验证了算法的有效性。

图4 轨迹估计结果
Fig.4 Results of trajectory estimation

图5 系统误差估计的误差均值
Fig.5 Error mean of systemic error estimation

图6 系统误差估计的误差标准差
Fig.6 Error standard deviation of systemic error estimation

为比较算法的配准性能,基于上述仿真场景,在相同条件下对文献[24]提出的适用于被动传感器的MLR算法(算法1)和本文算法(算法2)的配准性能进行对比分析。图7给出了在迭代次数为5、采样点数为400条件下对卫星1的蒙特卡罗仿真统计结果。由图7可知,算法2在方位偏差估计和俯仰偏差估计上的误差标准差均低于算法1,方位偏差估计精度提高了约12%,俯仰偏差估计精度提高了约70%。综上可知,通过引入卫星轨道定向误差、姿态角测量误差和传感器测量误差的先验信息对目标状态的误差协方差进行修正,能有效降低量测中包含的随机误差对系统误差配准造成的不利影响,进一步提高了配准精度。

图7 算法性能对比
Fig.7 Performance comparison between algorithms

5 结 论

对低轨星座协同探测弹道目标情况下的多光学传感器配准问题进行了研究,提出了多星载光学传感器系统误差MLR算法。通过引入量测转换所涉及的测量误差的先验信息对目标状态的误差协方差进行修正,并基于视线交叉获得观测在状态空间中的近似投影,从而将MLR算法扩展至多星载光学传感器的配准应用中。仿真结果表明了该算法的有效性,且配准性能优于传统的被动传感器MLR算法。由于算法本身是利用期望极大化迭代实现系统误差估计的,实时性较差,同时未对观测异步以及多目标等复杂情形进行讨论,这也是多星载光学传感器的配准下一步亟待解决的问题。