基于遗传算法的极短弧定轨∗

李鑫冉 王 歆

(1中国科学院紫金山天文台南京210008)(2中国科学院空间目标与碎片观测重点实验室南京210008)(3中国科学院大学北京100049)

基于遗传算法的极短弧定轨∗

李鑫冉1,2,3†王 歆1,2‡

(1中国科学院紫金山天文台南京210008)
(2中国科学院空间目标与碎片观测重点实验室南京210008)
(3中国科学院大学北京100049)

空间目标的巡天观测获取了海量的极短弧观测数据,而经典初轨计算方法对于极短弧几乎不能获得合理的结果.将初轨计算问题转换为两个三变量的分层优化问题,采用遗传算法,针对具体问题选择了优化变量以及相应的遗传操作,建立了一种极短弧初轨计算方法.基于实测资料的数值实验表明,方法可为后续工作提供有效的初值.

航天器,天体力学,方法:数值

1 引言

短弧定轨,也称初轨计算或轨道计算,是指没有先验信息的情况下,根据单站单圈的少量观测数据计算二体意义下的轨道.测角资料的轨道计算是天体力学中的经典问题,两百多年来实测数据的轨道计算方法研究从未停止过,并且随着新的天体以及新的观测手段的出现,不断涌现出新的需求和困难[1].

经典的轨道计算方法是由Gauss和Laplace针对小行星的轨道计算问题建立的,人造卫星上天以后有了较大发展,Escobal对其进行了详细的介绍[2].Ta ff对各种方法进行了分类比较[3].我国学者也做了较多的研究,主要集中在改进Laplace方法方面,并结合我国自主的实测资料情况做了一些改进[4-10].

初轨计算问题始终没有得到满意的解决,根本上是由其求解方程的本征病态引起的,短弧定轨问题是不适定的.面对这种困难,为了能够获得较好的结果,只能依赖采集尽量长的弧段,对于光学观测5～6 min的弧段才足以获得较好的结果[11],在实践中一般要求最少连续采集3 min以上的弧段.

随着空间目标数量的不断增加,近年来空间目标观测中新出现了针对空域的观测模式,通过对空域的持续观测,记录下经过该空域的所有目标的方位信息.这种方式改变了原来一台设备只能同时对单个目标进行跟踪的情况,对同时经过空域的目标可以全部进行监视.这种方式提高了观测效率,大幅度提高了观测数量,但与此同时也带来了新的问题,即使采用大视场望远镜,对于低轨道空间目标能够采集到的弧长也是十分有限的,例如中国科学院空间目标与碎片光学观测网采用的超大视场望远镜采集到的弧段大多不超过30 s,相当多的仅有10余秒.近年来国际上许多大型巡天项目也都将空间目标作为其观测对象,这些项目采集了海量数据,极短弧定轨问题就成了有效利用数据的关键[12].

经典方法对于如此短弧段的初轨计算几乎无法适用,主要体现在计算不收敛或者得到完全不合理(unrealistic)的解(轨道半长径小于地球赤道半径).这样短的弧段被称为“极短弧”(TSA,Too Short Arc),以区别传统意义上的短弧[13−14].极短弧的具体弧长尚无明确定义,文献[11]认为2 min以下弧段为极短弧,而文献[14]则采用了1 min的弧段,一般认为极短弧指采用经典方法很难获得合理解的弧段,而具体弧长可以根据具体对象确定.

除了经典计算方法外,轨道计算中还有一类方法为优选法,采用优选法可克服经典方法中迭代不收敛的现象.但常应用在一些特殊情况,特别是可转换为一维优选问题的情况,如圆轨道的初轨计算等.对于多维情况,计算方法复杂不实用,应用极少.另一方面这些方法通常需要一个较为准确的初值,而初值选取本身就是一个初轨计算问题,对于极短弧轨道计算问题也不适用[1].

Ansalone等[14]提出将遗传算法用于解决极短弧定轨问题,由此可以避免迭代不收敛以及传统优选法需要较为准确的初值和多维问题求解复杂的问题.他们针对天基模拟资料采用双r优选法,构造了遗传算法的一个具体实现,但其采用的遗传运算与一般规律很不一致.Hinagawa等[15]将该方法用于同步卫星的一般初轨计算,并对部分遗传运算环节做了调整.

本文也采用遗传算法用于极短弧定轨问题,采用了完全不同于现有方法的变量和运算构造方法,实现了一种具体方法,并结合中国科学院空间目标与碎片光学观测网的实际数据情况做了参数选择和计算验证.

2 遗传算法

遗传算法(Genetic Algorithm,GA)是模拟达尔文生物进化论的自然选择的计算模型,通过模拟自然进化优胜劣汰的过程搜索最优解,是一种启发式方法.遗传算法自上世纪70年代由Holland和De Jong提出以来,得到了蓬勃发展,应用广泛,特别是算法的数学原理已得到充分研究,其收敛性和全局搜索能力已得到证明[16].

遗传算法将问题的解向量编码为染色体,解向量中的每个元素称为染色体上的基因.算法模拟生物进化的过程,在解的取值范围内随机产生G个染色体作为初始种群(population),对每个染色体进行适应度( fitness)评估,以衡量解的优劣性.通过选择(selection),交叉(crossover),变异(mutation)遗传运算,使染色体得到进化,获得适应度更好的新种群,通过不断进化从而求解出问题的最佳近似解.其计算流程见图1.

遗传算法中种群中个体的数目称为种群数,种群数越大其搜索能力也越强,但计算效率也越低.用来进行适应度评估的函数称为适值函数,函数值称为适值,即适应度越好,越接近最佳近似解.选择操作是为了从当代种群中选出优秀的个体,使其可以作为父代繁衍子代,将优秀的基因传递下去.选择依据适应度的大小,适应度越高,被选中的概率越大.交叉和变异是遗传算法中最重要的两个运算步骤,交叉操作通过对两个父代染色体进行组合从而产生新的子代个体,交叉操作使得优秀的基因有机会组合在一起,从而获得更优的个体.对于通过选择操作得到的子代,以一定的交叉概率Pc再进行交叉操作,交叉操作决定了遗传算法的全局搜索能力,一般取值大多在0.6～0.9之间.单纯通过交叉操作无法产生新的基因,遗传算法通过变异操作来产生新的基因,以实现较强的局部搜索能力.对于经过选择和交叉后的子代各染色体的基因以变异概率Pm进行变异操作,与生物进化类似,变异发生概率较低,一般取值在0.01～0.1之间.经过选择、交叉和变异后完成一代的进化.如此迭代进化,遗传算法将收敛到全局最优解.对于实际计算总需要一个终止条件,一般当最优解的适应度达到给定阈值或者进化代数(迭代次数)达到指定的最大代数作为终止条件.

对于不同问题,遗传算法流程是一致的,重点在于结合具体问题特点设计各种运算,关键在于避免种群出现早衰,即使得解陷入局部极值,而不是全局最优解.

图1 遗传算法的流程图Fig.1 The flowchart of genetic algorithm

3 初轨计算的遗传算法

3.1 变量和编码的选择

对于优选问题,首先需要明确优化变量.对于轨道计算问题最容易想到的是以轨道根数作为优选变量,同时优选6个变量使得参数空间维数太高,不利于求解.基于这个原因,文献[14]将问题转换为两个变量的优化问题,优选变量选择为观测首末时刻的斜距.这种选择必须将优选变量结合首末观测资料共同得到轨道,计算受到首末观测资料精度的制约,特别当首末资料出现较大误差时,无法得到合理的结果.因此我们选择和观测量不相关的历元时刻t0的3个Kepler根数(a,e,M0)作为优选变量,这样参数空间维数不太高,同时解不和具体观测资料相关,便于后续资料处理.

遗传算法常采用二进制编码,对于二进制编码,交叉、变异算子比较容易设计,但对于轨道计算这种对于精度要求较高的连续问题,采用二进制编码将使编码长度过长,因此本文采用了实数编码,这样无需再做转换[17].

3.2 初始种群的生成和终止条件

为了生成初始种群,首先根据一些已知信息定义3个优选变量的值域,有a∈[al,au], e∈[el,eu],M0∈[Ml,Mu].不同于传统优选方法需要较为准确的初值,取值范围并无特殊要求,可根据已知信息确定,没有太多信息时,可将范围取得大一些,这是容易做到的.初始种群采用均匀随机采样方式建立,对于每个优化变量,由

可得到一个取值,其中X为任一优化变量,λ为[0,1)上均匀分布的随机数.对于3个优化变量分别重复G次上述过程,即可得到种群数为G的初始种群{xi,i=1,2,···,G}.

为了简化计算流程,选择迭代次数达到最大进化代数K为终止条件.

3.3 适值函数

令已知一组观测量{ti,αi,δi,i=1,2,···,N},其中(αi,δi)为赤经和赤纬.根据历元时刻t0的(a,e,M0),对于任意时刻ti,可得到Mi:

其中Zi为ti时刻目标的天顶距.对任意一对观测时间(tk,tj),tk＞tj,即可得到(rk,rj)和相应的(fk,fj).定义适值函数:

对于每个个体xi,其适值为:

显然,对于轨道计算问题,适值越小越优.

适值除了表明解的优劣,在遗传算法中最主要的作用是决定每个个体被选择进入子代的概率,因此需要对适值进行一定的变换.对于讨论的问题,使得适值越小的个体有越大的选择概率,这个过程称为适值的标定(Scaling)[18]7-21.

遗传算法中全局搜索能力主要取决于种群的多样性,基于适值F的选择方式,容易在进化前期使得一些相对较好的解成为超强个体,被大量选择进入子代,降低种群的多样性;而在进化后期由于种群相对集中,个体间适值差异较小,从而变为随机搜索.本文选择基于适值的秩作为选择操作的依据,即标定后的适值和原适值取值无关,只与其排序有关.对于每代种群,对每个个体计算其适值Fi后,从小到大排序,令个体xi的排名为ri.显然最优解的ri=1,最劣解的ri=G.标定后的适值取

其中c为常数,用于维持种群的多样性,可根据具体问题调整.c＜1时可减小各适值间的差异,增加样本多样性,使种群中较差个体的选择概率不会太低;而c＞1时则增大了各适值间的差异,使种群中较差个体的选择概率更低.由于标定后的适值和具体取值无关,显然在进化的各阶段都有助于保持种群的多样性.

3.4 选择操作

选择操作的原则是适值越大的被选择的概率越大,对于个体xi,选择概率Pi取为:

可见经过适值标定后,最优个体的选择概率是确定的,从而可使不太优的解也有一定的被选择的概率.选择操作根据选择概率可采用轮盘赌(Roullete wheel)算法实现,考虑采用一般轮盘赌算法对于一次具体选择操作,随机因素会使得较优个体没有被选择,为了避免这种情况,采用随机均匀采样方法(Stochastic uniform sampling)实现选择操作.随机均匀采样方法等价于只随机转动一次轮盘,之后的转动按照预设的规则进行.

初轨计算问题是要获得最优解,并不需要种群中所有个体都达到最优,因此在选择操作后引入精英保留策略,即对于进化中最优的Ne个个体总是选择进入子代,采用精英保留策略后,每代的最优个体不会发生衰退.

3.5 交叉操作

由于采用了实数编码,无法直接使用类似二进制编码的交叉操作,文献[2]仿照二进制单点交叉方法,通过直接交换两个父代相关解向量元素的方法实现交叉.本文采用了算术交叉的方法[18]42-44,令和分别是父代(k−1代)选为交叉的两个个体,则子代(k代)的两个个体为:

其中,α为[0,1)上均匀分布的随机数.对于M0∈[0,2π],为角度量,因此可有两个方向进行交叉,约定总是在劣弧(＜π)上进行交叉.从上式可看出,采用的交叉运算为凸变换,即

交叉新生成的子代一定满足约束条件,避免了子代个体有效性的判断.

3.6 变异操作

遗传算法中变异算子决定了解的局部搜索能力,文献[14]采用了高斯变异.高斯变异的变异范围较难确定,文献[14]是直接给定了10 km作为变异范围.本文为了减少参数,采用了动态变异方法[18]42-44,令xi,j表示第i个个体的第j个基因,变异操作为:

其中β为[0,1)上均匀分布的随机数,决定变异发生的方向.和分别表示第j个基因的上下界.Δ函数具体形式为:

γ为[0,1)上均匀分布的随机数,k和K分别是当前进化代数和最大进化代数,b为常数,主要调整全局搜索能力.从式中可看出Δ(k,x)随k增大逐渐趋于0,这样在进化早期有较大的搜索范围,而到后期搜索范围较小,增强细调能力.显然当b＜1时倾向于全局开拓能力,而b＞1时则更倾向于尽快进入局部搜索.同样变异算子也可以确保变异后代符合约束条件.

3.7 (i,Ω,ω)的求解

上述各节实现了一个完整的遗传算法,经过求解后可得到t0时刻的(a,e,M0),由(5)式可得到ti时刻的地心位置矢量ri,在此基础上求解完整的6个轨道根数已无任何困难．

仍可采用和上述遗传算法相同过程求解(i,Ω,ω).求解中(a,e,M0)作为已知量, (i,Ω,ω)作为优化变量,仍是三变量优化问题.由(a,e,M0,i,Ω,ω)按照二体问题模型可计算ti时刻的,只需将适值函数变为:

由于已有一定信息,初始种群也无需在全区域搜索,根据每对(ri,rj)可得到一组(i,Ω),取其平均可得到参考的作为初值范围选取依据,值域可取

4 数值验证

根据上述过程,采用MATLAB1http://www.mathworks.com编写了求解程序,并选取中国科学院空间目标与碎片光学观测网中的一圈实测资料进行了数值验证,资料精度优于5′′,采样频率1 Hz.为了便于比对,选取弧段开始的不同弧长进行了计算.历元选择为第1个观测时刻,考虑到空间目标探测中绝大多数目标都是近圆轨道(e＜0.003)[1],值域选择为a∈[1.03,4.15], e∈[0,0.003]以及M0∈[0,2π],对于a并无利用任何先验信息,取值范围已覆盖了整个中低轨道区域;初轨计算中一般轨道面相对较准,取δi=0.1◦,δΩ=1.5◦.

算法中的参数采用了常见取值,G=50,K=50,Ne=2,Pc=0.8,Pm=0.05,对于初轨计算更侧重全局搜索能力,因此取(9)式中c=0.5,(14)式中b=0.5.分别计算了10～60 s弧长的结果,图2和图3分别给出了一组计算过程中从初始种群开始到进化结束,每代最优个体的适值F和轨道半长径a的结果.从图中可看出收敛速度非常快,从第10代起,适值F已经很接近最小值;a的变化也类似,第10代以后也基本稳定在准确解附近的较小区域内,算法搜索效率较高.

图2 适值F的收敛过程Fig.2 The convergence process of the fitness value

表1给出了不同弧长下的计算结果,其中POD表示由多天多站精密定轨得到的结果,作为参考标准,Laplace为改进Laplace方法得到的结果,GA为本文方法得到的结果.从表中可见,对于Laplace方法,40 s弧长时虽能获得结果,但半长径偏差已近100 km,当弧长缩短到20 s时,已无法得到合理的结果.采用本文方法对于各种弧长均得到了合理的结果,在10 s弧长下结果依然较好.

图3 轨道半长径a的收敛过程Fig.3 The convergence process of the semi-major axis a

从遗传算法的计算过程可见,计算结果依赖随机数的生成.为了验证随机数对计算的影响,避免随机数对极短弧定轨结果产生普适的影响,采用不同的随机数种子进行了多组计算.图4为10 s弧长采用不同随机数种子计算100次得到的轨道半长径直方图.从结果看,采用不同随机数得到的结果有一定差异,但这并不是遗传算法不收敛的体现,而是初轨计算问题本身病态的体现.图4中a仍在合理区间内,远优于极短弧定轨本身能够达到的精度.文献[19-20]给出的35 s弧长的定轨精度只有79 km,而70 s弧段的定轨精度也仅为16 km.表2和表3分别给出了20 s和30 s弧长计算中的10组结果,也说明了这点. 10组计算结果有一定差异,但最终残差都非常小,相互差异小于0.1′′,在精度范围内这些解应当是等价的.对实测资料的短弧定轨问题而言解并不唯一,极短弧情况下更为突出,遗传算法的全局搜索能力使得不同组只是得到了不同的最优解.在实践中可进行多组计算,选择适值相对较小或者多组平均结果作为定轨结果.

为了评估初轨计算精度,文献[19-20]引入Bootstrap方法实现了不依赖其他信息的精度估计方法,方法结合经典迭代过程的特点,通过对残差的重采样构建伪观测量进行精度评估.而遗传算法本身的特点使得残差重采样方法不再有效,因此采用参数化的Bootstrap方法对计算精度进行估计[21],在原始观测资料上根据测量误差分布来构建伪观测量用于精度评估.对观测资料增加5′′的随机误差后采用本文方法定轨,采用不同随机数种子对上述过程重复100次.由于原始观测资料中已有误差,增加5′′随机误差后得到的伪观测量误差应大于5′′,因此结果可作为5′′误差资料定轨精度的上界.由于不同随机数种子的使用,结果中包含了资料误差和随机数对求解过程的两重影响因素,可全面反映求解的实际精度.由于定轨精度与弧长密切相关,选择了最短的10 s弧长进行数值计算,得到的均值和标准差列于表4.结果表明采用本文方法定轨结果的精度特征与经典方法相仿,轨道面精度略优,沿迹量精度略差.对于10 s弧长的定轨结果,半长径a精度优于10 km,该精度已可为后续轨道改进以及轨道识别提供一个有效的初值.

表1 各弧长下的定轨结果Table 1 The results of orbit determination for various arc lengths

5 结论与讨论

综上所述,本文方法克服了经典方法对于极短弧定轨的一些不足,能够应用于极短弧定轨,计算无需先验初值,具有普适性.本文采用的遗传运算和参数选择与遗传算法的常规选取保持一致,保证了其求解的稳定性.

遗传算法是基于种群的算法,其计算效率较经典方法大大降低,当测量数据弧段较长时,经典方法已可获得很好的结果,遗传算法毫无优势.但对于极短弧而言,计算效率上的差异就当前计算能力而言已可忽略,在经典方法失效时,仍可获得有效的定轨结果,为之后的工作提供初值.

本文将定轨问题转换为分层优化问题,分解为两个三参数优化问题,在优化变量仅增加一个的基础上,解不再依赖具体的观测量,提高了方法对于资料误差和野值的容错能力,也便于各种资料处理方法的应用.

此外,在确定(a,e,M0)时放弃了对轨道面的约束,可使得对残差较为敏感的沿迹方向相关量的计算精度有所提高,由于初轨计算中得到较为准确的a是关键,这种处理对于定轨结果的后续应用显然是有帮助的.

图4 不同随机数种子下的半长径的分布Fig.4 The histogram of the semi-major axis with di ff erent random seeds

表2 采用不同随机数种子下20 s弧长的轨道根数Table 2 The results of 20-second arc length with di ff erent random seeds

在初轨实践中固定根数等处理都是十分有效的手段[22],在原有方法中引入这些约束需要在前几次迭代中得到合理的结果或者较好的初值,这些对于极短弧定轨是很难满足的,而采用遗传算法后这些处理都可以很自然地融入其中,没有任何困难,大大方便了实际使用.

表3 采用不同随机种子下30 s弧长的轨道根数Table 3 The results of 30-second arc length with di ff erent random seeds

表4 轨道根数的均值和标准差Table 4 The mean and standard errors of elements

[1]吴连大.人造卫星和空间碎片的轨道和探测.北京:中国科学技术出版社,2011:158-191

[2]Escobal P R.Methods of Orbit Determination.New York:John Wiley&Sons,1965:239-289

[3]Ta ffL G.AJ,1984,89:1426

[4]贾沛璋,吴连大.天文学报,1997,38:353

[5]陆本魁,戎鹏志,吴建民,等.宇航学报,1997,18:1

[6]贾沛璋,吴连大.天文学报,1998,39:337

[7]Jia P Z,Wu L D.ChA&A,1999,23:384

[8]陆本魁,李剑峰,马静远,等.宇航学报,1999,20:14

[9]刘林,王歆.天文学报,2003,47:175

[10]Liu L,Wang X.ChA&A,2003,27:335

[11]Vallado D A,Carter S S.JAnSc,1998,46:195

[12]Milani A,Knezevic Z.CeMDA,2005,92:1

[13]Tommei G,Milani A,Rossi A.CeMDA,2007,97:289

[14]Ansalone L,Curti F.AdSpR,2013,52:477

[15]Hinagawa H,Yamaoka H,Hanada T.AdSpR,2014,53:532

[16]De Jong K A.Evolutionary Computation—A Uni fied Approach.Cambridge:MIT Press,2006:26-27

[17]Michalewicz Z.Genetic Algorithm+Data Structures=Evolution Programs.Berlin:Springer,1992: 97-105

[18]玄光男,程润伟.遗传算法与工程优化.北京:科学出版社,2000

[19]王歆.天文学报,2013,54:73

[20]Wang X.ChA&A,2013,37:338

[21]Efron B,Tibshirani R.An Introduction to the Bootstrap.New York:Chapman&Hall,1993:53-56

[22]贾沛璋,吴连大.紫金山天文台台刊,1997,16:175

Genetic Algorithm for Initial Orbit Determination with Too Short Arc

LI Xin-ran1,2,3WANG Xin1,2
(1 Purple Mountain Observatory,Chinese Academy of Sciences,Nanjing 210008)
(2 Key Laboratory for Space Object and Debris Observation,Purple Mountain Observatory,Chinese Academy of Sciences,Nanjing 210008)
(3 University of Chinese Academy of Sciences,Beijing 100049)

The sky surveys of space objects have obtained a huge quantity of tooshort-arc(TSA)observation data.However,the classical method of initial orbit determination(IOD)can hardly get reasonable results for the TSAs.The IOD is reduced to a two-stage hierarchical optimization problem containing three variables for each stage. Using the genetic algorithm,a new method of the IOD for TSAs is established,through the selection of optimizing variables as well as the corresponding genetic operator for speci fic problems.Numerical experiments based on the real measurements show that the method can provide valid initial values for the follow-up work.

space vehicles,celestial mechanics,methods:numerical

P135;

:A

10.15940/j.cnki.0001-5245.2016.01.007

2015-07-16收到原稿,2015-08-21收到修改稿∗国家自然科学基金项目(11373072)资助

†lixr@pmo.ac.cn

‡wangxin@pmo.ac.cn