杜荣华, 张 翔, 廖文和

(南京理工大学机械工程学院, 江苏 南京 210094)

0 引 言

仅测角相对导航并不是一个新概念,最早应用于19世纪科学家通过望远镜观测天体并计算其轨道的研究中。当前,仅测角相对导航已经在航海[1]、深空探测[2]、卫星轨道确定[3]以及编队飞行[4-5]任务中得到了广泛应用和研究。

仅测角相对导航的概念相当简单,即通过单个相机在一段时间内测量主动星与空间目标之间的视线角(line of sight,LoS)就可以确定其相对运动状态[6]。目前,常用来测量空间物体之间相对运动状态的有源主动式传感器(例如微波雷达、激光雷达等),由于其功耗高、体积大和造价贵,使得其在实际工程应用中受到很大限制,而光学相机等无源被动式传感器(例如红外相机、可见光相机等)在这一方面具有很大的优势[7-9],因为可应用于各种卫星间距,并且对卫星质量、功率的设计影响很小[10],非常适合微小卫星平台的应用。实际上,许多航天器配备了星敏感器,如果方向适当,可用于跟踪其视野范围内的空间目标,并根据LoS观测值执行相对导航任务[11]。其中,仅测角相对导航技术代表了多种高级分布式空间系统的明确使能技术,包括自主交会和对接、空间态势感知、分布式孔径科学以及非合作目标的在轨维修等任务[12]。此外,仅测角相对导航技术在天基反导监测以及空间碎片探测和清理等领域也有较大的应用前景[13-14]。

在过去十年内,已经有在轨项目利用仅测角相对导航技术进行了非合作目标交会接近的演示验证。瑞典的原型研究仪器和太空任务技术进步(prototype research instruments and space mission technology advancement,PRISMA)计划是使用仅测角相对导航技术实现非合作目标远距离自主交会的典型代表。2011年4月,PRISMA主承包商OHB-Sweden在自主交会(autonomous rendezvous, ARV)实验中实现了从30 km至50 m的基于LoS的自主会合[15-17]。2011年8月,在PRISMA运行即将结束时,德国宇航中心(Germany Space Operation Center,GSOC)操作追踪星Mango从60 km抵近目标星TanDEM-X至4 km处,抵近过程仍然使用基于LoS的相对导航滤波器[18-21]。2012年4月,GSOC使用全球定位系统和光学导航进行了高级交会论证(简称为ARGON)实验[22],在ARGON实验中,地面研发的专用飞行动力学系统可对星上收集的图像进行处理,并根据单目相机测得的LoS矢量估计追踪星和空间目标之间的相对轨道信息,进而实施一系列的轨道机动,在5天内完成了从30 km至3 km的抵近操作。2016年11月,GSOC进行了自主视觉接近导航和目标识别(autonomous vision approach navigation and target identification,AVANTI)实验[23],与ARGON实验相比,AVANTI实验中的图像处理、相对导航以及轨道机动规划均为星载自主实现,而ARGON实验中需要将追踪星拍摄的图像数据传回地面,在地面处理完成后再将轨道机动指令发送到卫星上。AVANTI实验共进行了两次基于仅测角相对导航的抵近操作,第一次抵近操作中将两颗卫星距离从13 km减少至1 km,然后依靠空气阻力将卫星间距拉大至2.5 km处,在第二抵近操作中将两颗卫星的距离减小至50 m以内[23]。

由于较差的动态可观性,AGRON和AVANTI两个实验都使用了先验初始解来初始化仅测角相对导航滤波器,即北美航空航天防御司令部(简称为NORAD)提供的两行根数(two line elements,TLE)和地面雷达测轨信息来确保仅测角相对导航滤波器的收敛[22-23]。从严格意义上讲,整个相对导航过程并不是完全自主的。为了提高仅测角相对导航系统的自主性,研究仅测角相对导航初始相对轨道确定(initial relative orbit determination, IROD)方法是实现非合作目标完全自主抵近的关键。然而,由于单目相机深度信息的缺失,仅测角相对导航存在星间距离不可观测的问题[24-25]。目前,已经有许多学者研究了提高仅测角相对导航可观性的方法,例如轨道机动法[26-27]、相机偏置法[28-29]、非线性动力学法[30-31]和多星多传感器法[32-33],其中相机偏置法只能应用于近距离范围内,而轨道机动法会增加卫星的燃耗,多星多传感器法则存在多卫星多传感器配置复杂的缺陷。

针对上述各种方法存在的问题,本文拟在无机动操纵观测弧内,引入非线性相对运动模型,该非线性相对运动模型基于相对轨道根数(relative orbit elements, ROE)建立追踪星和空间目标之间的相对运动状态,并将星间距离和相对轨道形状进行解耦。然后,在线性理论获得的共线性解附近系统地改变星间距离大小,并执行一系列最小二乘拟合,随后利用二分法或牛顿迭代法快速在全局范围内寻找最小拟合残差的最优解。最后,在4种不同的轨道场景中对提出的方法进行了半物理仿真校验,仿真结果证明所提出的方法可以较快地获得高精度的仅测角相对导航IROD的解。

1 问题描述

仅测角相对导航的目的在于推导出空间目标和主动星之间的相对运动状态,该相对运动状态对应于ti时刻获得的一组LoS观测值{ui:i=1, 2,…,k},这些观测值是指向空间目标的单位LoS矢量。一般情况下,仅测角相对导航的问题与系统的可观性密切相关。Woffinden和Geller[34]已经证明,在线性化的测量模型和相对运动模型的假设下,仅测角相对导航系统是不可观测的。

首先,令x(t)表示空间目标和主动星在时间t处的相对运动状态,则线性化的相对运动模型的演化方程可表示为

x(t)=Φ(t,t0)x(t0)

(1)

式中,Φ(t,t0)是状态转移矩阵。

假设相对位置矢量r(t)与相对运动状态x(t)之间具有线性相关性,其关系可表示为

r(t)=L(t)x(t)

(2)

一组给定的LoS观测值{ui:i=1, 2,…,k}和相对位置矢量r(t)之间满足以下关系:

ui×r(ti)=0,i∈[1,k]

(3)

将式(1)和式(2)代入式(3)中,得到

ui×(L(ti)Φ(ti,t0)x(t0))=0,i∈[1,k]

(4)

从式(4)中可以看出,如果x(t0)=x0是式(4)的解,则缩放解φx0也是式(4)的解,从而导致一组给定的LoS观测值具有无限个解,造成系统不可观测,这个被称为Woffinden困境[35]。因此,为了通过仅测角相对导航获得空间目标和主动星之间的相对运动状态,必须解决Woffinden困境,提高仅测角相对导航系统的可观性。

2 算法说明

由第1节的Woffinden困境可知,采用线性化的相对运动模型和测量模型将导致仅测角相对导航系统不可观测,因此本文引入了非线性相对运动模型。首先,为捕获线性模型所忽略的细微差异,重新建立相对运动状态更新方程,可记为

x(t)=f(t,x(t0))

(5)

式中,f是非线性函数,对应于包含J2摄动的相对运动方程的数值积分。

测量模型h可表示为

(6)

式中,r(t,x0)为t时刻空间目标相对于主动星在直角笛卡尔坐标系中的相对位置矢量。可令r(t,x0)=g(x(t)),其中g一般为非线性转换函数。

假设ui=[ux(ti),uy(ti),uz(ti)]是由相机采集的图像获得的空间目标在ti时刻的LoS单位矢量,该LoS单位矢量也可用一组角度测量值(即方位角α和俯仰角ε)进行表示,其中i=1, 2,…,k,表示第i个测量时刻。理想情况下,矢量ui与h(ti,x0)在方向上是平行的,即

ui×h(ti,x0)=0,i∈[1,k]

(7)

由于ui中含有LoS矢量测量误差,h(ti,x0)含有建模误差,两者实际有偏差,导致ui×h≠0。因此,IROD问题可以描述为最优化问题:

(8)

式中,J为损失函数,表示前k个时刻LoS矢量误差的平方和。

由于仅测角导致的弱可观性,如何设计优化算法保证式(8)收敛到全局最小值是找到精确IROD解的关键。在线性化的相对运动模型中,空间目标相对于主动星之间的相对运动状态可以表示为

x(t)=Φ(t,t0)x0

(9)

假设相对位置矢量r(t,x0)和x(t)之间具有线性相关性,可记为

r(t,x0)=C(t)x(t)

(10)

式中,C(t)为r(t)与x(t)之间的转换矩阵。

对于k个LoS矢量测量值ui,由式(7)可构造以下线性方程组:

ui×(C(ti)Φ(ti,t0)x0)=0,i∈[1,k]

(11)

(13)

比较式(8)与式(13)的寻优过程,可见寻优搜索空间从6个维度的x0减少到只有一个维度的φ。图1给出了该IROD研究方法的流程图。需要指出的是,缩放解φx0中φ的上限和下限取决于主动星和空间目标之间的距离范围。对于本研究的应用场景,主动星和空间目标的间隔大约介于1～100 km之间,通过在1～100 km范围内改变比例因子φ,就可以找到使拟合残差均方根σ(φ)最小的比例因子φ,即为仅测角相对导航IROD的全局最优解。

图1 仅测角相对导航IROD方法流程框图

证明σ(φ)为凸函数

实际上,可以通过使用一组ROE或一组笛卡尔坐标系的坐标分量来参数化相对运动状态x。在这种情况下,相对运动模型f描述了相对运动状态x的时间演变,而函数g将相对运动状态x映射为直角笛卡尔坐标系中的坐标分量。相对运动状态x的建模可采用曲线型Hill-Clohessy-Wiltshire(HCW)模型[36]和Gim-Alfriend模型[37]。后者基于ROE,但提供了一个线性映射L,以将ROE转换为曲线笛卡尔坐标系表示的相对运动状态(即x(t)=L(t)O(t)ω(t,t0)O(t)L(t)x0=Φ((t,t0)x0),其中ω为基于平均ROE的状态转移矩阵,O为平均轨道根数和密切轨道根数之间的转换矩阵),可以根据Brouwer[38]和Lyddane[39]中提出的理论进行平均轨道根数至密切轨道根数的一阶非线性映射。

根据式(3),可以描述一般非线性的测量方程为

ui×r(ti)=ui×(g(ti,f(ti,x0)))=0,i∈[1,k]

(14)

进一步简化为

ui×(g(ti,f(ti,x0)))=0,i∈[1,k]

(15)

假设曲线型相对运动状态在轨道坐标系中进行表示,并且分量按径向-切向-法向序列排序。对于正在考虑的问题(即相对运动主要为沿轨道方向的星间距离),直线相对位置r可以由曲线相对位置rcur得到

(16)

式中,R是圆形轨道的半径;rcur,2表示相对位置沿轨道方向的分量,且R≥rcur,2,所以式(16)可以简化为

(17)

把状态转移矩阵拆分成两块,得到Φ=[Φ1-3,1-6,Φ4-6,1-6]T。让Φ6,1-6定义矩阵Φ的第6行,且r=g(x)=g(Φx0)=Φ1-3,1-6x0。因此,式(15)可以采取以下形式进行表示:

(18)

因此,对式(15)的一般非线性方程,该问题已简化为二次函数的最小化问题,为简单起见,引入与卫星轨道曲率模型相对应的二次部分:

(19)

ui×(Φ1-3,1-6x0)=0,i∈[1,k]

(20)

并扩展成

(21)

式中,Φ1-3,1-6=[Φ1-3,1-5,Φ1-3,6]已拆分为多个块。

累加k次测量后,获得线性系统(由表示线性问题的下标L标识):

(22)

其五维解可由最小二乘法给出:

(23)

0,i∈[1,k]

(24)

使用式(21)中采用的相同方法将问题限制在五维空间中:

(25)

(26)

在累加k个测量值之后,再次获得另一个线性系统:

(27)

Γ(φ)T(I-P(φ))Γ(φ)

(28)

为简化起见,引入投影矩阵P(φ)=A(φ)(A(φ)T·A(φ))-1A(φ)T,得出任意固定为φ时损失函数J的最小二乘意义的解析公式,但其凸性仍有待证明。可以通过考虑以下因素来简化式(25)以便近似得

(29)

(30)

进而矢量Γ可以表达成以下形式:

Γ= sign(φ)(Γ1+φΓ2)

(31)

(32)

(33)

3 数值分析

3.1 相对运动状态参数化

引入一组无量纲的ROE来进行相对运动状态的参数化[40]:

δα=[δa,δex,δey,δix,δiy,δλ]T

(34)

式中,δa是无量纲的相对半长轴;δλ代表相对平均经度;δe=[δex,δey]T和δi=[δix,δiy]T分别为相对偏心矢量和相对倾斜矢量。

与笛卡尔坐标系表示的相对运动状态参数相比,通过ROE进行参数化具有以下优点:可以快速了解相对运动的几何形状,如图2所示,描述了局部笛卡尔坐标系(径向、切向和法向或RTN)中表示的相对运动状态,其单位矢量定义如下:eR和eN分别与主动星的绝对位置和轨道角动量方向对齐,eT与eR和eN满足右手定则。

此外,从图2可以看出[35],相对运动的大小可以通过ROE的尺寸(即轨道半长轴a)来描述[41]。

图2 笛卡尔坐标系和ROE表示的相对运动

这组ROE特别适用于仅测角相对导航问题,因为弱可观测的星间距离几乎与aδλ分量是相等的,从而将相对运动的形状和星间距离从几何上进行了解耦。

3.2 相对运动模型

需要建立主动星和空间目标相对运动的精确模型,以捕获系统微弱的非线性效应,提高仅测角相对导航系统的可观性。本文使用第3.1节中介绍的ROE(即x=δα)对相对运动模型进行参数化。在这种情况下,g不再是线性的。为了获得笛卡尔坐标系中的相对位置,首先需要将平均ROE转换为密切ROE,从而恢复J2摄动引起的长期和短期效应,然后将密切相对轨道根数映射到笛卡尔坐标系表示的相对位置,并且允许在相对运动模型中引入其他经验参数,例如反映差动阻力扰动效应的参数[42-43]。

在数学上,该模型通过状态转移矩阵Φ将时间t处的相对运动状态δα(t)与时间t0处的相对运动状态δα(t0)相关联,即

δα(t)=Φ(t,t0)δα(t0)

(35)

(36)

(37)

(38)

3.3 IROD的快速求解方法

由于使用ROE进行相对运动状态的参数化,因此矩阵C表示的是ROE和直角笛卡尔坐标系表示的相对位置之间的转换关系[35]:

C(t)=

(39)

式中,u代表主动星的纬度幅角;i代表主动星的轨道倾角。

在最小二乘拟合残差意义上对应的特定值S=aδλ=φ的线性解可推导如下:

(40)

Γ=-SA6

(41)

(42)

(43)

式中,κi代表κ的第i个分量。

3.3.1 二分法

3.3.2 牛顿迭代法

牛顿迭代法的思想是将非线性函数(原方程)线性化(切线方程),以线性方程的解逐步逼近非线性方程的解,其原理图如图3所示。

图3 牛顿迭代法原理图

其原理公式为

(44)

牛顿迭代法的具体步骤如下:

4 仿真验证

4.1 地面半物理仿真平台搭建

在完成快速仅测角相对导航IROD算法的设计后,搭建验证该算法的地面半物理仿真平台,其结构框图如图4所示。

图4 地面半物理仿真平台结构框图

该地面半物理仿真平台主要包括以下部分。

(1)轨道/姿态发生器

该模块内部装有Matlab/Simulink软件和STK软件,主要用于生成主动星和空间目标的轨道/姿态数据。

(2)星光星图模拟计算机

该计算机作为星光星图模拟器,内部装有导航星库、星图生成软件等,主要用来模拟恒星和空间目标在星图中的分布情况及星光等级。

(3)星图识别和目标检测计算机

该计算机内部装有星图识别软件、星体质心定位软件和空间目标自主检测软件等,主要功能是在原始图像中找到空间目标,并输出空间目标相对于摄像机的LoS测量信息。

(4)仅测角相对导航计算机

该计算机的主要功能是完成仅测角相对导航算法的实时运算,并输出空间目标相对于主动星的相对运动状态信息。

(5)实时监视显示计算机

该计算机主要用于监视和显示轨道/姿态发生器模拟的主动星和空间目标的轨道/姿态数据以及仅测角相对导航计算机的输出结果等。

搭建的实际地面半物理仿真平台如图5所示。在原有设备的基础上,配备了三轴精密转台和三轴直线导轨等机械微调装置,可调整星图显示装置和星图采集模块同轴。

图5 地面半物理仿真平台

4.2 实验测试与结果分析

在搭建的地面半物理仿真平台上对该方法进行校验,主要仿真参数如表1所示。

表1 主要仿真参数

该仿真校验在4种不同的轨道场景中进行,包括远距离保持点相对运动轨迹(场景1)、自然漂移相对运动轨迹(场景2)、远近距切换保持点相对运动轨迹(场景3)和主从随动相对运动轨迹(场景4)。这4个轨道场景刚好包含了主动星从远距离交会至抵近空间目标的各个过程,其相对运动轨迹如图6所示。出于模拟的需要,主动星和空间目标的参考运动使用20×20阶的重力场模型的数值积分传播得到,其中包括第三体引力、太阳辐射压力和大气阻力扰动等,因此可以得到LoS测量的参考值以及主动星和空间目标之间无机动的观测弧,最终根据主动星和空间目标模拟的状态创建一组观测值。

图6 4种不同轨道场景的相对运动轨迹图

此外,为了评估快速仅测角相对导航IROD算法的运行时间,仅测角相对导航计算机选用的是常规配置的Intel Core i5系列计算机。测试结果如表2～表3和图7～图14所示,为该算法在4种不同轨道场景中的性能。其中,表2为不同方法单次运算时耗;表3为采样周期为30 s,采样次数为5次时,不同方法的百分比误差;图7～图14为不同采样周期和观测次数情况下,各方法在不同轨道场景中的IROD性能。

图7 不同观测次数下的IROD方法的性能(场景1,T=10 s)

图8 采样周期下的IROD方法的性能(场景1, 观测次数为5)

图9 不同观测次数下的IROD方法的性能(场景2, T=10 s)

图10 不同采样周期下的IROD方法的性能(场景2,观测次数为5)

图11 不同观测次数下的IROD方法的性能(场景3, T=10 s)

图12 不同采样周期下的IROD方法的性能(场景3,观测次数为5)

图13 不同观测次数下的IROD方法的性能(场景4, T=10 s)

图14 不同采样周期下的IROD方法的性能(场景4,观测次数为5)

表3 IROD方法的百分比误差

从表2中可以看出,在相同轨道场景中,批量搜索法单次运算时耗大于牛顿迭代法和二分法,并且牛顿迭代法单次运算时耗最低。在不同轨道场景中,场景4的单次运算时耗明显高于其他3种轨道场景,其次是场景2,场景1和场景3两种轨道场景的单次运算时间相差不大,这是由于场景1和场景3的轨道场景相似,只是距离不同导致的。

表2 单次运算时耗

为了评估算法性能,定义均方根误差、标准差和百分比误差分别为

(45)

(46)

(47)

本文将相对平均经度aδλ的均方根误差Mx、标准差σx和百分比误差Px作为评价算法性能的指标,因为仅测角相对导航的不可观测性主要集中在该参数上[45]。

分析表3、图7、图9、图11和图13,在采样周期相同但观测次数不同的情况下,在相同轨道场景中,牛顿迭代法的IROD性能最好,其相对定轨误差可达所测距离的6%左右,其次是二分法和牛顿迭代法。此外,不同的观测次数也会影响算法的IROD性能,但并不是次数越多性能越好,在观测次数大于一定数量后,算法性能趋于稳定。

分析表3、图8、图10、图12和图14,在观测次数相同但采样周期不同的情况下,采样周期会影响算法的IROD性能,其中采样周期越大,算法的IROD性能越好。这是由于大的采样周期可以增加系统的非线性和扰动效应,并抵消相机噪声产生的影响,从而增加了仅测角相对导航的可观性。比较不同轨道场景中算法的IROD性能,在场景4的轨道场景中,各方法的IROD性能明显低于其他3种轨道场景,其相对定轨误差都在所测距离的30%以上。这是由于场景4的系统可观性较差导致的,这说明主动星和空间目标之间的相对运动轨迹也会影响所提算法的IROD性能。

5 结 论

本文针对仅测角相对导航IROD问题,提出了一种快速仅测角相对导航IROD方法,并通过地面半物理仿真平台测试了该方法在4种不同轨道场景中的性能,得到以下结论。

(1)采用基于ROE建立的非线性相对动力学模型,通过最小二乘拟合残差的方法可以获得较高精度的仅测角相对导航IROD的解,且性能受采样周期和观测次数的影响。增大采样周期可以提高所提方法的IROD的性能,这是由于大的采样周期可以增加系统的非线性和扰动效应,并抵消相机噪声产生的影响,从而提高系统的可观性。增加观测次数也可以提高所提方法的IROD的性能,但并不是次数越多性能越好,在观测次数大于一定数量后,算法性能趋于稳定。因此,在实际工程应用中,应该合理选择采样周期和观测次数。

(2)提出的二分法和牛顿迭代法可以快速地获得仅测角相对导航IROD的解。但是,所提出的方法在主从随动轨道场景中的IROD性能明显低于其他3种轨道场景,这是由于主从随动轨道场景中的系统可观性较差所导致的,所以未来将继续开展仅测角相对导航可观性最优的轨迹设计的研究。