赵柯昕 甘庆波 刘 静

(1 中国科学院国家天文台北京100101)

(2 国家航天局空间碎片监测与应用中心北京100101)

(3 中国科学院大学北京100049)

1 引言

天基空间目标观测具有全天候、全天时、灵活机动等优势, 受到了国内外极大关注. 20世纪90年代至今, 已经发射了十余颗天基观测卫星, 如天基监测系统(Space-Based Surveillance System,SBSS)1、地球同步空间态势感知计划(Geosynchronous Space Situational Awareness Program,GSSAP)2和Sapphire卫星[1]等. 现阶段天基目标观测仍以光学载荷为主要观测设备, 测角资料为主要观测数据, 获得目标轨道参数为主要任务目标, 而初始轨道确定是实现这一目标的基础.

天基空间目标光学观测的初轨确定问题具有两个基本特征, 光学短弧段导致的几何约束不足问题和观测平台与空间目标动力学特征相似导致的奇点(平凡解)问题. 前者是初轨确定中普遍存在的问题, 而天基观测中由于相对运动速度更快导致观测弧段更短. 平凡解问题首先由甘庆波等人于2007年提出[2], 发现将类Laplace方法应用于天基光学观测的初轨确定中会出现计算结果收敛到观测平台自身轨道的问题, 并推导了平凡解存在的数学形式, 给出了初步解决方法. 2009年南京大学刘林等人确认了改进Laplace方法收敛到观测平台本身的问题[3]. 2018年, Doscher[4]指出针对近圆、共面轨道进行定轨时更容易收敛到平凡解.

理论上可以通过求解由角度测量数据构造的条件方程, 进而计算轨道根数. Gauss首次将初轨确定问题转换为求解非线性方程的问题, 提出了使用多组测角数据的最小二乘(least square, LS)法[5].Escobal[6]于1976年提出了双r迭代法(r为观测时刻目标地心距), 该方法可以适用于任何时间间隔的观测数据, 但存在收敛困难、鲁棒性较差等问题.Gooding[7–8]在其对Lambert问题研究的基础上, 提出了一种基于Lambert问题和Newton-Raphson法使用3次测角数据的初轨确定新方法, 通过迭代求解Lambert方程搜索探测目标的斜距. Vallado[9]于2010年评估了Gooding方法应用于天基空间目标观测任务中的性能, 发现该方法相当稳健, 可以解决大多数不同轨道构型的情况. 从算例来看Gooding方法对斜距的计算已经接近真实值. 1994年Dumoulin[10]提出了一种基于Stumpt函数和Lambert方程使用位置向量的迭代方法, 但该方法在倾角和斜距过小等特殊情况时会计算失败. 2008年陈务深等人将天基初轨确定问题建模为求解不等式约束的非线性最小二乘问题, 由此将该问题转换为计算最优解的数学问题[11–12]. 2017年章品等人提出了一种基于距离搜索的天基空间目标初轨确定方法[13]. Gooding[7–8]、Dumoulin[10]和章品等[13]的初轨计算方法都需要一个较合理的初始猜测值. 一旦猜测值与实际情况偏差过大, 通常会不收敛或收敛到观测平台轨道. 2010年Shefer[14]将利用两个观测时刻距离量进行初轨确定的基本问题转化为求解非线性方程的问题, 给出了椭圆、双曲线和抛物线3种情况下方程具体形式和初值计算方法, 并使用Newton-Raphson法计算轨道根数.2019年Kuznetsov[15]在Shefer[14]工作的基础上, 进一步推导出从测角数据求解斜距的非线性方程组,利用小行星实测数据进行验证, 取得了良好的效果.

本文分析了天基光学目标观测的初轨确定问题, 基于改进的Gauss方程, 提出了一种适用于天基光学观测数据的初轨确定新方法. 首先简要比较了天基观测和地基观测在观测构型上的异同, 分析了天基初轨确定过程中计算收敛到平凡解的数学原理. 从利用空间目标两次观测时刻的位置矢量求解半通径的问题出发, 通过改进经典Gauss方程, 将利用光学观测数据的初轨确定问题转换为了求解关于斜距的非线性方程组的问题, 推导了该方程组的解析形式. 而本质上该方程组与观测数据组数无关, 同样适用于多组观测资料的情况. 利用轨道能量约束减小了解的搜索区域, 消除了方程组的奇点. 使用天基目标观测的实测数据分析了本方法构建的非线性条件方程组解的性质. 最后利用天基低轨观测平台对低、中、高轨目标的仿真观测数据验证了本方法具有良好的性能.

2 天基光学观测的初轨确定问题

地基和天基光学观测的初始轨道确定问题的几何构型如图1所示. 图1中空间目标地心位置矢量r、观测平台地心位置矢量R以及空间目标相对观测平台斜距矢量ρ可以组成一个三角形. 下标S和G分别表示天基观测和地基观测的情况. 空间目标的地心矢量r可表示为:

图1 地基和天基光学观测的几何构型Fig.1 The geometrical configuration of ground-based and space-based optical surveillance

式中,L为观测视线方向的单位矢量,ρ为斜距. 观测视线方向的单位矢量为:

式中(α,δ)分别为空间目标在地心天球坐标系下的赤经和赤纬角度测量数据.

对(1)式进行求导, 得到目标地心位置矢量的一阶和二阶导数分别为:

地基观测与天基观测的几何构型基本相同, 不同之处在于观测平台的运动特征. 当观测平台位于地表时, 其遵循地表测站的运动规律.

式中,RS是观测平台的地心距. 联立(1)、(4)、(6)式后, 经典Laplace方法或Gauss方法在计算八次方程时会出现一个公因子(r-RS). 由此会出现一个平凡解r=RS. 将平凡解代入(1)式, 得到观测平台的轨道根数. 这是经典初轨确定方法收敛到平凡解的数学本质, 是观测平台和空间目标动力学相似造成的[2].

3 天基初轨确定新方法的原理和特点

定义在地心天球坐标系下, 3次观测时间tj=1,2,3, 对应的3次赤经、赤纬角度测量值(αj,δj)j=1,2,3、观测平台的地心位置矢量Rj=1,2,3、空间目标的斜距ρj=1,2,3与地心位置矢量rj=1,2,3和空间目标相对观测平台观测视线方向的单位矢量Lj=1,2,3. 从求解轨道半通径的问题出发, 为了避免经典Gauss方程出现奇点, 对其进行改进. 由此将使用天基光学测角数据的初轨确定问题转换为求解观测时刻斜距的非线性条件方程组的问题.

3.1 半通径计算方程

当观测时间间隔较小时, 可以认为3次观测时刻的空间目标地心位置矢量位于同一平面内[5]:

当已知t1和t2时刻的空间目标相对观测平台的斜距时, 可求出t2时刻的斜距ρ2:

各符号定义如下:

因此, 可以通过(1)式和角度测量数据得到空间目标在观测时刻的地心位置矢量(r1,r2,r3). 由于测量弧段较短, 地心位置矢量r1与r2、r1与r3、r2与r33者相互之间的夹角都小于π. 则半通径p可通过下式得到:

式中,r1、r2和r3分别是地心位置矢量r1、r2和r3的模. 现在考虑由空间目标在t1和t3时刻的地心位置矢量(r1,t1;r3,t3)求解半通径的问题. 位置矢量r1和r3可以和椭圆轨道围成一个扇形和一个三角形[6], 如图2所示, 其中点A与B分别为t3与t1时刻目标所在的位置, 角度ν1与ν3分别为t1与t3时刻目标的真近点角,点O为地心.此时位置矢量r1、r3与弧�AB可组成一个扇形, 其面积为AS. 而两次观测时刻的位置矢量与线段AB可组成一个三角形, 其面积为AT. 由图可知扇形与三角形的面积比是一个小值.

图2 空间目标位置矢量组成的扇形和三角形Fig.2 The triangle and sector formed by the position vectors of space target

引入辅助变量y表示扇形面积AS和三角形面积AT的比:

其中,τ13=t3-t1. 当计算出面积比y后, 可以较为方便地计算半通径. 为求解y引入经典Gauss方程[6], 如下(11)–(12)式:

式中

对于椭圆轨道, 函数X和x可表示为:

式中, dE=E3-E1,E1、E3为t1、t3观测时刻空间目标的偏近点角. 将(11)式代入(12)式:

(11)和(12)式需要迭代计算,Escobal[6]于1965年给出了一种迭代计算方法. 当地心位置矢量(r1;r3)已知时, 可计算得到参数m和l. 第1次循环中, 可以定y= 1, 函数x和X使用(14)和(13)式依次求出. 再利用(15)式求出新面积比y. 重复上述计算过程, 直至达到预定精度. 对于观测弧段较短, 即两次观测时刻的地心位置矢量(r1;r3)夹角较小时, 具有很快的收敛速度. 但当夹角较大或夹角等于π时无法求出面积比.

3.2 改进Gauss方程

为了避免除以零或无穷出现奇点, 将经典Gauss方程(11)式改写为[15–16]:

式中

当r1和r3重叠时,σ等于0, 但在天基观测中是很罕见的情况. 因此可以排除σ= 0的情况. 联立(11)和(12)式, 可得

式中

对于椭圆轨道, 函数X可表示为:

最终, 将求解半通径的问题转化为求解x的问题.

3.3 非线性条件方程

定义时间间隔τ12=t2-t1和τ23=t3-t2, 关于(t1,ρ1;t3,ρ3)的非线性条件方程可以将其代入(17)式得[15]:

式中

8

由(1)式和(9)式, 函数x12和x23可表示为:

由此构造出了含有t1和t3时刻斜距(ρ1,ρ3)两个未知数的两个方程, 求解该方程经过后续计算可得t2观测时刻斜距,由此可通过观测数据计算出3次观测时刻空间目标的地心位置矢量. 推导出的约束方程(21)式本质上与观测数据组数无关, 同样适用于多组观测资料的情况.

3.4 利用轨道能量约束缩小搜索范围

利用迭代方法搜索f1(ρ1,ρ3) = 0和f3(ρ1,ρ3)= 0的根时计算效率较低, 可以使用轨道能量约束将搜索范围控制在一个较小的区间. 空间目标一般在椭圆轨道上运行, 应满足由Milani等人提出的地心二体能量约束[17]:

式中,ε为轨道能量, ˙ρ是斜距ρ的一阶导数. 斜距ρ应为正值. 同时因为光学望远镜的能力限制,ρ应存在最大值ρmax. 因此, 斜距应满足以下约束条件:

3.5 非线性条件方程的奇点问题

在求解f1(ρ1,ρ3) = 0和f3(ρ1,ρ3) = 0时, 当p=0、f11= 0和f31= 0会导致方程(21)式无意义,从而导致计算失败. 通过分析, 当且仅当r1=-r2和r2=-r3时, 会导致f11和f31等于0, 但这在天基短弧光学观测时是不存在的. 仅剩下需要避免p=0这一种情况. 给定θ1、θ2分别作为地心位置矢量r1和r2、r2和r3的夹角, 由(9)式可得:

当位置矢量r1、r2和r3满足(25)式时会使p=0, 导致条件方程是奇异的, 该点需要去除.

3.6 利用Gibbs法或Herrick-Gibbs计算轨道根

数

迭代求解非线性条件方程(21)式, 得到t1和t3时刻的斜距ρ1、ρ3. 使用(8)式求得t2时刻的斜距ρ2.由此可得空间目标在3次观测时刻的位置矢量(t1,r1)、(t2,r2)、(t3,r3). 当t1和t3时刻的观测方向矢量夹角大于1°时, 采用Gibbs法计算初始轨道根数.当夹角小于1°时, Herrick-Gibbs方法具有较好的精度[6]. 得到了t2观测时刻位置和速度矢量, 经过转换即可得到经典Kepler轨道根数.

3.7 初轨确定方法流程

针对天基光学观测初轨确定问题, 提出了使用本文方法的初始轨道确定流程:

(1)根据先验信息给出t1和t3时刻的斜距初始猜测值ρ1和ρ3, 利用赤经、赤纬角度观测数据, 结合(8)式得到t2时刻的斜距ρ2;

(2)利用(1)式求解3次观测时刻空间目标地心位置矢量r1、r2和r3, 利用(9)式计算轨道半通径p;

(3)利用(23)式计算参数x12与x23, 由此构建非线性条件方程组(22)式;

(4)重复步骤(1)至步骤(3), 更新ρ1和ρ3, 直至(22)式中的f1(ρ1,ρ3)=0和f3(ρ1,ρ3)=0同时成立,利用(8)式得到ρ2;

(5)计算3次观测时刻的目标位置矢量, 根据观测视线矢量的夹角, 选择使用Gibbs法或Herrick-Gibbs法确定t2观测时刻目标的速度矢量, 经过转换获得目标的Kepler轨道根数.

4 数值验证

4.1 基于实测数据的定轨结果

使用光学观测卫星对空间目标的角度实测数据来分析非线性条件方程组根的特性. 经过平滑后, 得到了3次观测时刻的赤经α与赤纬δ角度测量数据和该观测卫星位置矢量(Rx,Ry,Rz), 如表1所示. 观测时刻用简化儒略日(Modified Julian Date,MJD)表示, 角度测量数据和卫星地心位置矢量位于J2000地心天球坐标系中, 测量精度为10′′.

表1 天基光学观测平台对空间目标3次观测数据Table 1 Three observation measurements of space targets by space-based optical observation platform

将表中数据代入(22)式, 得到了关于(ρ1,ρ3)的非线性方程组. 将f1(ρ1,ρ3) = 0和f3(ρ1,ρ3) = 0的曲线绘制于以第1、3观测时刻斜距(ρ1,ρ3)为横轴和纵轴的图中, 如图3和图4所示. 图中虚线为f1(ρ1,ρ3) = 0, 实线为f3(ρ1,ρ3) = 0. 图3为(22)式在8倍地球半径(ρ1,ρ3)∈[0,8]的范围内图像情况. 图4为图3中的交点在(ρ1,ρ3)∈[0,1]的放大情况. 地球半径定义为REarth= 6378.14 km. 两条曲线的交点认为是研究范围内(22)式的解.

图3 条件方程在8倍地球半径下的交点Fig.3 The intersection points of condition equation under eight times Earth radius

图4 条件方程在1倍地球半径下的交点Fig.4 The intersection points of condition equation under one time Earth radius

由图可以发现在8倍地球半径的搜索区域内,方程组仅存在1个交点, 即仅存在唯一解. 图中的阴影部分为不满足轨道能量约束的区域. 交点1所对应的第1、3次观测时刻的斜距量分别为0.274452REarth和0.258642REarth, 通过进一步求解可得该目标的Kepler轨道根数分别为:

其中,a、e、i、ω、Ω、M分别为轨道半长轴、偏心率、轨道倾角、近地点俯角、升交点赤经和平近点角.

4.2 不同观测情况的仿真验证与精度评估

为了进一步比较本文提出的天基初轨确定方法的性能, 使用天基低轨观测平台对低轨近圆轨道目标、中轨近圆轨道目标、大椭圆轨道目标和地球同步轨道目标分别进行仿真观测. 使用本文方法和经典Gauss法对仿真观测数据进行初始轨道计算, 对比两种方法在使用不同长度观测弧段时的位置计算误差Δr与速度计算误差Δv. 表2展示了MJD = 59410.166667时刻下, 所使用的观测平台、低轨(Low Earth Orbit, LEO)、中轨(Middle Earth Orbit, MEO)、大椭圆轨道(Highly Elliptical Orbit, HEO)与同步轨道(Geostationary Earth Orbit, GEO)目标的Kepler轨道根数.

表2 观测平台与目标的Kepler轨道根数Table 2 The Kepler orbit elements of observation platform and targets

观测平台和空间目标的轨道外推考虑了由美国国家地理空间情报局(US National Geospatial-Intelligence Agency, NGA)发布的EGM 2008 40阶×40阶的地球重力场模型和太阳、月球三体引力.天基观测平台的观测噪声主要是观测平台自身轨道误差、观测系统自身误差等因素引起的. 为了评估初轨确定方法, 观测数据增加5′′的Gauss随机噪声. 将t1时刻仿真目标的位置和速度矢量作为真值.

4.2.1 LEO平台观测LEO目标的情况

图5为两种方法使用不同长度的LEO目标观测数据计算得到的初始轨道和该目标位置、速度矢量的残差. 观测弧段长度从4 s至30 s, 每秒进行一次计算. 图中实线和虚线分别为本文方法和经典Gauss方法初始轨道计算结果的残差随观测弧段长度的变化趋势. 从图中可以发现, 当观测弧段较短且小于10 s时本文方法依旧有良好的收敛性, 可以收敛到真解附近. 当观测弧长大于10 s时本文方法和经典Gauss法具有相近的计算精度.

图5 观测弧长为4–30 s时LEO目标初轨计算结果Fig.5 The IOD (initial orbit determination) result of the LEO target with observational time intervals from 4 to 30 s

4.2.2 LEO平台观测MEO目标的情况

图6为两种方法使用不同长度的MEO目标观测数据计算得到的初始轨道和该目标位置、速度矢量的残差. 观测弧段长度由30 s至300 s, 每间隔5 s进行一次计算. 可以发现当观测弧段长度小于60 s时本文方法具有更好的精度. 观测弧段较长时本文方法拥有和Gauss法相似的精度.

图6 观测弧长为30–300 s时MEO目标初轨计算结果Fig.6 The IOD result of the MEO target with observational time intervals from 30 to 300 s

4.2.3 LEO平台观测HEO目标的情况

图7为两种方法使用不同长度的HEO目标观测数据计算得到的初始轨道和该目标位置、速度矢量的残差. 观测弧段长度由90 s至300 s, 每间隔5 s进行一次计算. 可以发现当观测弧段较短且小于100 s或观测弧段长度大于250 s时, 相比于经典Gauss法本文方法具有更好的精度.

图7 观测弧长为90–300 s时HEO目标初轨计算结果Fig.7 The IOD result of the HEO target with observational time intervals from 90 to 300 s

4.2.4 LEO平台观测GEO目标的情况

图8为两种方法使用不同长度的GEO目标观测数据计算得到的初始轨道和该目标位置、速度矢量的残差. 使用观测弧段长度为60–360 s, 每间隔10 s进行一次计算. 可以发现此时对于GEO目标本文方法和经典Gauss法具有相似的精度.

图8 观测弧长为60–360 s时GEO目标初轨计算结果Fig.8 The IOD result of the GEO target with observational time intervals from 60 to 360 s

4.2.5 轨道根数计算结果

选取不同目标不同长度的观测数据, 使用本文方法进行轨道确定, 并计算经典轨道根数. LEO目标的观测弧长为15 s, MEO目标的观测弧长为120 s, HEO目标的观测弧长为160 s, GEO目标的观测弧长为180 s时, 并对观测数据增加5′′的随机误差后, 使用本方法计算得到的不同目标轨道根数如表3所示. 观测弧长没有特殊选择, 仅是典型的观测值. 将表3中的计算结果与表2中的真值进行比较可知, 在不同的观测弧长下, LEO目标、MEO目标、HEO目标和GEO目标的半长轴误差分别约为50 km、15 km、55 km和2 km.

表3 空间目标初轨确定计算结果Table 3 The IOD results of space targets

5 总结与讨论

天基空间目标观测是未来空间态势感知领域的一个重要发展方向. 本文研究了天基观测数据的初始轨道确定问题, 并提出了一种基于改进Gauss方程的初轨确定计算方法. 与其他方法相比, 本方法的关键是构造观测时刻空间目标斜距量的非线性方程组. 文中利用3组角度测量数据推导了该方程组的解析形式,也可适用于多组观测资料的情况.使用天基光学观测卫星对LEO目标的实际测量数据分析了方程组解的形式. 最后利用天基光学观测平台分别对低、中、高轨和大椭圆轨道目标进行了仿真观测, 得到不同长度的角度测量数据. 通过该仿真测量数据, 分析了本文方法与经典Gauss方法对不同弧长数据的计算结果, 发现对于较短的观测弧段, 本文方法具有更好的精度, 由此表明了本文提出的天基空间目标光学观测初轨确定方法具有一定的优越性.