基于二叉樹的太陽系小天體數(shù)值歷表建立方法?

胡壽村季江徽

(1中國科學(xué)院紫金山天文臺(tái)南京210008)

(2中國科學(xué)院行星科學(xué)重點(diǎn)實(shí)驗(yàn)室南京210008)

(3中國科學(xué)院大學(xué)北京100049)

1 引言

小天體探測是21世紀(jì)全球深空探測的熱點(diǎn),歐美日目前都已經(jīng)對(duì)部分小天體(包括小行星、彗星和矮行星)展開過實(shí)質(zhì)性的探測活動(dòng),嫦娥二號(hào)也在2012年12月13日近距離飛越了4179 Toutatis小行星并取得重要的科學(xué)成果[1–3].目前新的小天體探測計(jì)劃(日本的Hayabusa-2、美國的OSIRIS-REx任務(wù))也在開展當(dāng)中,并且在2022年美國還將發(fā)射兩次小行星探測任務(wù)Lucy和Psyche,分別去探測6顆木星特洛伊小行星以及金屬型主帶小行星16 Psyche[4?5].

在小天體精確任務(wù)軌道設(shè)計(jì)中,需要計(jì)算小天體在完整力模型積分下的軌道.此時(shí)從效率和精度考慮,可以對(duì)事先積分好的軌道做擬合或插值來獲得高精度的光滑連續(xù)軌道,以計(jì)算任意時(shí)刻目標(biāo)天體的坐標(biāo).常用方法包括拉格朗日插值、Neville插值、切比雪夫多項(xiàng)式擬合、樣條插值等[6?7],其中切比雪夫多項(xiàng)式因具有較高的數(shù)值穩(wěn)定性而被廣泛應(yīng)用于構(gòu)建自然天體和人造衛(wèi)星的數(shù)值歷表[8?9].

目前已有多篇文獻(xiàn)介紹了切比雪夫多項(xiàng)式用于構(gòu)建大行星/月球和人造衛(wèi)星(如GPS衛(wèi)星)數(shù)值歷表的方法[10–12],其核心思想是將軌道做等區(qū)間分段后用一定階數(shù)的切比雪夫多項(xiàng)式進(jìn)行分段擬合,并將每一段擬合的切比雪夫多項(xiàng)式系數(shù)保存下來.在調(diào)用歷表時(shí),只需根據(jù)時(shí)間搜尋到對(duì)應(yīng)的系數(shù)段后即可反演出該時(shí)刻天體的軌道.然而與大行星/月球以及大部分人造衛(wèi)星不同,許多小天體的軌道偏心率較大,且存在近距離飛越大行星的可能,使得軌道隨時(shí)間的變化并不均勻.圖1給出了目前15630顆近地小行星和924顆彗星軌道偏心率的分布,分別有超過62.2%的近地小行星和超過81.7%的彗星軌道偏心率大于0.4[13]1.圖2給出了99942 Apophis近地小行星2000–2050年的軌道圖,其中左圖為日心平黃道系下Apophis的軌道,右圖為Apophis的地心距隨時(shí)間的變化.可以看到,雖然該小行星偏心率較小,但由于2029年4月13日Apophis將以31600 km的軌道高度近距離飛越地球,導(dǎo)致在此時(shí)間前后的軌道發(fā)生了較大變化[14?15].對(duì)于這部分小天體,它們不宜用傳統(tǒng)的等區(qū)間分段的方法做擬合或插值,否則容易造成較大的歷表誤差范圍.因此合理的方案是在保持?jǐn)M合階數(shù)不變的前提下在軌道變化較劇烈的位置處用更小的區(qū)間做擬合或插值.但這樣帶來的問題是如果區(qū)間分段不按照一定的規(guī)則,那么根據(jù)時(shí)間計(jì)算星歷時(shí)難以查找到其在歷表文件中所處的位置,因此需要一種新的數(shù)據(jù)結(jié)構(gòu)來組織區(qū)間分段信息以提高查找效率.

圖1 近地小行星(左)和彗星(右)的偏心率分布(數(shù)據(jù)來源于國際小行星中心1http://www.minorplanetcenter.net/iau/Ephemerides/EphemOrbEls.html)Fig.1 Eccentricity distribution for near-Earth asteroids(left)and comets(right)(data from IAU Minor Planet Center1)

本文引入計(jì)算機(jī)科學(xué)中的二叉樹結(jié)構(gòu)來表示不均勻分段的區(qū)間,在誤差較大處按照二分法對(duì)區(qū)間分段,直到滿足給定精度約束,然后調(diào)用星歷時(shí)在給定時(shí)刻可根據(jù)二分法查找來定位其系數(shù)段的位置.該工作可作為建立大量太陽系小天體數(shù)值歷表的統(tǒng)一方法,相比于傳統(tǒng)等時(shí)區(qū)間分段方法能夠有效壓縮數(shù)值歷表的存儲(chǔ)量,可應(yīng)用于大量目標(biāo)小天體的精確軌道設(shè)計(jì)中.

本文安排如下:首先,第2節(jié)介紹切比雪夫多項(xiàng)式擬合用于構(gòu)建大行星/月球數(shù)值歷表的基本方法;然后,第3節(jié)介紹二叉樹結(jié)構(gòu)的概念以及我們將其應(yīng)用于構(gòu)建小天體數(shù)值歷表的具體方法;第4節(jié)將通過對(duì)3顆大偏心率小天體和99942 Apophis小行星做數(shù)值歷表仿真來說明本文介紹的方法在小天體中的應(yīng)用.

圖2 99942 Apophis小行星2000—2050年的軌道(左)及地心距(右)Fig.2 The orbit(left panel)and geocentric distance(right panel)of 99942 Apophis from 2000 to 2050

2 建立數(shù)值歷表的基本方法

2.1 切比雪夫多項(xiàng)式

在數(shù)學(xué)上,Chebyshev多項(xiàng)式是以遞推方式定義的一系列多項(xiàng)式,計(jì)算效率很高,且具有良好的數(shù)值穩(wěn)定性[16?17].第1類Chebyshev多項(xiàng)式遞推公式如下:

第2類Chebyshev多項(xiàng)式遞推公式為:

以上(1)–(3)式中x的范圍均為?1≤x≤1.切比雪夫多項(xiàng)式具體應(yīng)用于數(shù)值歷表構(gòu)建的基本方法在文獻(xiàn)[10]中有較為詳細(xì)的敘述,這里只做簡略介紹.

假設(shè)擬合階數(shù)為D,記X和V分別為位置和速度矢量中的任一分量(采用位置和速度作為軌道參數(shù)),ck為對(duì)應(yīng)的切比雪夫多項(xiàng)式系數(shù),則X和V的計(jì)算公式為:

(4)式中需要將實(shí)際的時(shí)間引數(shù)(儒略日J(rèn)D)轉(zhuǎn)換到歸一化時(shí)間范圍[?1,1],即:

這樣,接下來的任務(wù)只需求出系數(shù)ck.

2.2 切比雪夫多項(xiàng)式系數(shù)計(jì)算

采用類似建立JPL(噴氣推進(jìn)實(shí)驗(yàn)室)歷表的方法,在各分段內(nèi)按均勻歸一化的時(shí)間步長1/4從1到?1(即1,3/4,1/2,1/4,···,?3/4,?1)進(jìn)行采樣. 這樣總共有18個(gè)采樣點(diǎn)(9個(gè)位置和9個(gè)速度),切比雪夫多項(xiàng)式系數(shù)滿足如下條件:

T是一個(gè)18×(D+1)的矩陣,因此階數(shù)最高不能超過17,如果D＜17,則(6)式是一個(gè)超定方程組,可以用最小二乘法來求解.但為了使得最后生成的歷表隨時(shí)間連續(xù),則需要保證在區(qū)間兩端連續(xù),因此在(6)式的基礎(chǔ)上還需要增加以下4個(gè)約束條件:

采用拉格朗日乘子法可以求解(6)式和(10)式.忽略過程,整理后得:

W為18×18的權(quán)重矩陣(對(duì)角矩陣).各對(duì)角元素分別對(duì)應(yīng)各采樣點(diǎn)的位置或速度權(quán)重(奇數(shù)點(diǎn)對(duì)應(yīng)位置,偶數(shù)點(diǎn)對(duì)應(yīng)速度),一般取位置權(quán)重為1,并可通過增加或減小速度權(quán)重值來減小或提高速度擬合精度.在文獻(xiàn)[10]中建議的取值為:

本文算例中W的取值與之相同.從(14)式中求出c即得到該分段內(nèi)的切比雪夫多項(xiàng)式系數(shù)(l可舍棄).

2.3 數(shù)值歷表構(gòu)建和讀取流程

對(duì)于大行星/月球和大部分人造地球衛(wèi)星,由于軌道偏心率不大,軌道變化均勻,傳統(tǒng)的方法就是直接對(duì)軌道做等區(qū)間分段(JPL大行星/月球歷表,GPS衛(wèi)星星歷表等都用這種方法),生成歷表文件的主要流程如下:

(1)獲得原始軌道數(shù)據(jù),給定擬合階數(shù)和區(qū)間長度;

(2)對(duì)整條軌道做等區(qū)間分段;

(3)對(duì)各區(qū)間做切比雪夫多項(xiàng)式擬合,并保存獲得的系數(shù)到歷表文件;

(4)保存階數(shù)、區(qū)間分段等信息到歷表文件頭.根據(jù)歷表文件計(jì)算t時(shí)刻軌道的流程如下:

(1)載入歷表文件到內(nèi)存中;

(2)根據(jù)時(shí)刻t和歷表文件頭信息查找到對(duì)應(yīng)系數(shù)段;

(3)根據(jù)系數(shù)反演計(jì)算t時(shí)刻的軌道.

以上方法應(yīng)用于軌道隨時(shí)間均勻變化的小天體也是合適的,但是正如前文所說,有相當(dāng)一部分小天體的軌道變化可能是較不均勻的,這導(dǎo)致等區(qū)間分段得到的擬合精度變化范圍較大.若在誤差較大處縮小分段區(qū)間,但如果沒有給定一定的規(guī)則,則調(diào)用歷表時(shí)不易根據(jù)時(shí)間反查到對(duì)應(yīng)的系數(shù)段位置.因此本文考慮采用二分法來對(duì)區(qū)間做分段,采用二叉樹結(jié)構(gòu)來存儲(chǔ)分段信息,這樣在給定時(shí)刻后可以用快速的二分法查找來定位系數(shù)位置,下一節(jié)將詳細(xì)介紹這種方法.

3 二叉樹結(jié)構(gòu)在數(shù)值歷表中的應(yīng)用

在計(jì)算機(jī)科學(xué)中,二叉樹是一種每個(gè)結(jié)點(diǎn)至多有兩個(gè)分支的樹結(jié)構(gòu).樹中每一個(gè)數(shù)據(jù)元素稱為結(jié)點(diǎn),如果該結(jié)點(diǎn)包含有子結(jié)點(diǎn),則稱之為分支結(jié)點(diǎn),否則稱為葉子結(jié)點(diǎn).本文考慮一種特殊的二叉樹結(jié)構(gòu),該二叉樹的結(jié)點(diǎn)要么是葉子結(jié)點(diǎn),要么是含有兩個(gè)子結(jié)點(diǎn)的分支結(jié)點(diǎn)(根據(jù)國外文獻(xiàn)的定義,這種樹結(jié)構(gòu)又稱為滿二叉樹).

為了闡述二叉樹在建立數(shù)值歷表中的應(yīng)用,圖3給出了一段軌道做分割后表示成二叉樹結(jié)構(gòu)的示意圖.

圖3 軌道段作不規(guī)則二分法分段后(左圖)用二叉樹結(jié)構(gòu)表示(右圖),帶下劃線的為葉子結(jié)點(diǎn)Fig.3 An orbit section is divided into irregular pieces(left panel)by dichotomy,and represented with binary tree(right panel),the underline symbol means leaf node

圖3中,各數(shù)字為軌道區(qū)間的編號(hào),其中帶下劃線的為不再繼續(xù)分段的軌道段(即葉子結(jié)點(diǎn)).以圖3為例,結(jié)合二叉樹結(jié)構(gòu)建立數(shù)值歷表的步驟可敘述如下:

(1)獲得原始軌道數(shù)據(jù),給定擬合階數(shù)和歷表擬合精度p;

(2)評(píng)估1的擬合誤差,由于誤差大于p,將1等分成2和7,記錄此次分段信息;

(3)評(píng)估2的擬合誤差,由于誤差大于p,繼續(xù)將2分割成3和4,記錄此次分段信息;

(4)評(píng)估3的擬合誤差,由于誤差小于p,保存擬合得到的系數(shù);

(5)評(píng)估4的擬合誤差,由于誤差大于p,繼續(xù)將4分割成5和6,記錄此次分段信息;

(6)評(píng)估5的擬合誤差,由于誤差小于p,保存擬合得到的系數(shù);

(7)評(píng)估6的擬合誤差,由于誤差小于p,保存擬合得到的系數(shù);

(8)評(píng)估7的擬合誤差,由于誤差小于p,保存擬合得到的系數(shù);

(9)保存整個(gè)軌道段的區(qū)間分段信息并建立區(qū)間與系數(shù)段索引的關(guān)系;

(10)保存階數(shù)、區(qū)間分段等信息到歷表文件頭,結(jié)束.

要計(jì)算某時(shí)刻t的星歷,此時(shí)可以利用二分法查找,例如假設(shè)t在6區(qū)間內(nèi),則程序?qū)凑找韵虏襟E查找(忽略其他步驟):

(1)發(fā)現(xiàn)t在區(qū)間1內(nèi),因?yàn)?為分支節(jié)點(diǎn),繼續(xù)查找;

(2)繼續(xù)判斷,發(fā)現(xiàn)t在2內(nèi),因?yàn)?為分支節(jié)點(diǎn),繼續(xù)查找;

(3)繼續(xù)判斷,發(fā)現(xiàn)t在4內(nèi),因?yàn)?為分支節(jié)點(diǎn),繼續(xù)查找;

(4)繼續(xù)判斷,發(fā)現(xiàn)t在6內(nèi),因?yàn)?為葉子節(jié)點(diǎn),可根據(jù)葉子節(jié)點(diǎn)保存的系數(shù)段索引提取相應(yīng)的切比雪夫多項(xiàng)式系數(shù),結(jié)束查找.

以上根據(jù)t查找系數(shù)段索引的次數(shù)最多為二叉樹的深度,由于查找運(yùn)算只是簡單的判斷程序,耗費(fèi)的計(jì)算機(jī)時(shí)間可忽略.查找結(jié)束后,即可根據(jù)提取到的切比雪夫多項(xiàng)式系數(shù)計(jì)算天體在該時(shí)刻的位置和速度.從以上敘述可以知道,我們的方法(下稱為新方法)與傳統(tǒng)方法的區(qū)別在于生成歷表時(shí)增加了一個(gè)區(qū)間分段的操作,而讀取歷表時(shí)增加了一個(gè)根據(jù)時(shí)間反查系數(shù)段位置的操作.最終生成的歷表文件不僅包含文件頭信息和切比雪夫多項(xiàng)式系數(shù),還要保存分段信息(但分段信息占用的存儲(chǔ)量非常小),相比于傳統(tǒng)方法復(fù)雜了一些,但流程是清晰的.

若記分段數(shù)為P,則可以計(jì)算最終生成的歷表存儲(chǔ)的系數(shù)個(gè)數(shù)N為

通過下文的數(shù)值仿真可以看到,由于通過新方法減小了N的值,使得對(duì)于部分小天體可以大大減小歷表文件的大小,這對(duì)需要生成大量小天體數(shù)值歷表的情況(例如需要對(duì)大量候選目標(biāo)小行星計(jì)算發(fā)射窗口)是很重要的.

4 數(shù)值仿真

考 慮3顆 近 地 小 行 星4179 Toutatis、 3200 Phaethon、 99942 Apophis和 彗 星2P/Encke(其中,Toutatis、Phaethon和Encke均為大偏心率小天體,而Apophis將于2029年4月13日近距離飛越地球),其各自的軌道信息見表1的第2到4行.這里將采用本文介紹的新方法建立各自的數(shù)值歷表,取歷表的時(shí)間跨度為5 yr(2027–2032年),容許的最大位置誤差為1 km,擬合階數(shù)為10,生成的數(shù)值歷表信息如表1所示,其中第5行為區(qū)間分段的范圍,第6行為總的區(qū)間數(shù).可以看到:對(duì)于偏心率越大的天體,區(qū)間的變化范圍越大,而Apophis飛越地球?qū)е碌能壍雷兓沟米钚^(qū)間只有0.11 d.圖4以Phaethon和Apophis為例演示了區(qū)間分段對(duì)應(yīng)的軌道位置,顯然Phaethon的分段較密處發(fā)生在近日點(diǎn)前后,而Apophis發(fā)生在飛越地球前后.

這里與傳統(tǒng)方法生成的歷表做一比較,分別取4顆小天體的總分段數(shù)為表1中第6行的值(即分別取為20、53、44和32),采用均勻區(qū)間分段,生成的數(shù)值歷表的位置、速度誤差與新方法的對(duì)比如圖5所示.圖5的結(jié)果表明:對(duì)于3顆大偏心率小天體Toutatis、Phaethon和Encke,傳統(tǒng)方法的擬合誤差在近日點(diǎn)附近迅速放大,遠(yuǎn)日點(diǎn)和近日點(diǎn)的位置誤差分別可以相差2個(gè)、8個(gè)和6個(gè)量級(jí),而Apophis由于在飛越地球后軌道有較大的改變,導(dǎo)致位置誤差變化范圍也可以達(dá)到8個(gè)量級(jí).因此,傳統(tǒng)的等區(qū)間分段法確實(shí)不適用于建立這些小天體的數(shù)值歷表,否則為了達(dá)到1 km精度,根據(jù)表1第5行的結(jié)果,它們分別需要取相等的分段區(qū)間為28.5 d、1.8 d、0.11 d和7.1 d.相比于新方法,這會(huì)導(dǎo)致系數(shù)的存儲(chǔ)量迅速增大.表1第7行給出了在相同的1 km位置精度下,新方法相對(duì)于傳統(tǒng)方法所節(jié)省的存儲(chǔ)空間,顯然偏心率越大的小天體節(jié)省的空間越可觀,而Apophis的軌道特性甚至可以導(dǎo)致節(jié)省99.7%的存儲(chǔ)量.

表1 4顆小天體軌道生成的數(shù)值歷表信息Table 1 The information of numerical ephemerides for the four minor bodies

圖4 左圖為3200 Phaethon小行星在一個(gè)軌道周期內(nèi)的軌道分段,右圖為99942 Apophis在2029年1月—2030年3月的軌道分段,其中方框內(nèi)為Apophis飛越地球前后30 d的軌道段,圖中以紅色和藍(lán)色來區(qū)分相鄰分段.Fig.4 The demonstration of orbit division for 3200 Phaethon in one orbit period(left panel)and 99942 Apophis from 2029 Jan.to 2030 Mar.(right panel).The square frame in the right panel contains Apophis orbit in 30 d before and after the Earth- flyby.The red and blue colors are used to distinguish the adjoining sections.

對(duì)以上4顆小天體,若取不同的擬合階數(shù)D(D∈[6,16]),各自的分段數(shù)P和系數(shù)個(gè)數(shù)N隨D的變化如圖6所示.顯然隨著D的增大,P值逐漸減小,但對(duì)于不同的天體,最小N值對(duì)應(yīng)的D值也不同.為了最大程度減小系數(shù)個(gè)數(shù),對(duì)不同天體應(yīng)取不同的擬合階數(shù)D.從圖6結(jié)果來看,D從10到13取值并無顯著差異,但從計(jì)算效率角度考慮,D取值一般不應(yīng)過大.文獻(xiàn)[10]也給出了JPL歷表中各大行星和月球擬合采用的階數(shù),可供參考.

圖5 傳統(tǒng)方法(藍(lán)色)與新方法(橙色)生成的歷表的位置和速度誤差的對(duì)比Fig.5 The comparison of position and velocity errors of ephemerides for old method(blue)and new method(orange)

圖6 對(duì)這4顆小天體,取1 km的最大允許位置誤差,分段數(shù)P和系數(shù)個(gè)數(shù)N隨擬合階數(shù)D的變化.Fig.6 For the four minor bodies,assuming the allowable position error is 1 km,the figure shows the variation of section number P and coefficient number N with fitting order D.

5 總結(jié)

本文基于二叉樹結(jié)構(gòu)提出了一種通用的小天體數(shù)值歷表構(gòu)建方法.該方法繼承了傳統(tǒng)的分段軌道切比雪夫多項(xiàng)式擬合方案,但將傳統(tǒng)的等區(qū)間分段修改為非等區(qū)間分段,即依據(jù)區(qū)間內(nèi)的實(shí)際擬合精度來調(diào)整分段區(qū)間,使得最終所有區(qū)間內(nèi)的精度都滿足給定要求.該方法采用二分法做區(qū)間分段,最后采用二叉樹結(jié)構(gòu)來組織所有的區(qū)間分段信息,這樣在根據(jù)時(shí)間來反查對(duì)應(yīng)的切比雪夫多項(xiàng)式系數(shù)時(shí)能夠采用快速的二分法搜索進(jìn)行查找,具有較高的查找效率,解決了非等區(qū)間分段下的系數(shù)搜索問題.由于該方法能保證只在擬合誤差較大的位置(一般為軌道劇烈變化的位置)采用較小的區(qū)間分段,而其他位置仍可采用較大的區(qū)間,從而能減小歷表文件的大小而不影響整個(gè)歷表的擬合精度.此優(yōu)點(diǎn)在需要構(gòu)建大量小天體數(shù)值歷表時(shí)更能夠明顯體現(xiàn)出來.如果要將多顆小天體的歷表合并成一個(gè)歷表文件,只需在文件頭中多增加相應(yīng)信息,具體步驟不贅述.本文提供的方法尤其適用于構(gòu)建大偏心率或可能近距離飛越大行星的小天體數(shù)值歷表.文中以3顆大偏心率小天體4179 Toutatis、3200 Phaethon和2P/Encke以及將近距離飛越地球的99942 Apophis小行星為例,通過數(shù)值仿真證實(shí)了該方法的有效性.另外,由于該方法沒有對(duì)天體的軌道性質(zhì)做任何要求,因此既可作為構(gòu)建太陽系小天體數(shù)值歷表的通用方法,甚至也可用于構(gòu)建大偏心率人造衛(wèi)星的數(shù)值歷表.

[1]Huang J C,Ji J H,Ye P J,et al.NatSR,2013,3:3411

[2]Zhao Y H,Ji J H,Huang J C,et al.MNRAS,2015,450:3620

[3]Jiang Y,Ji J H,Huang J C,et al.NatSR,2015,5:16029

[4]Lauretta D S,Team O R.An Overview of the OSIRIS-REx Asteroid Sample Return Mission.43rd Lunar and Planetary Science Conference,Woodlands,March 19–23,2012

[6]魏二虎,柴華.全球定位系統(tǒng),2006,31:13

[7]趙齊樂,劉經(jīng)南,葛茂榮,等.武漢大學(xué)學(xué)報(bào)(信息科學(xué)版),2006,31:879

[8]陳劉成,賈曉林,莫中秋.天文學(xué)進(jìn)展,2006,24:167

[9]張如偉,劉根友.大地測量與地球動(dòng)力學(xué),2008,28:115

[10]Newhall X X.Numerical Representation of Planetary Ephemerides.Proceedings of the 109th Colloquium of the International Astronomical Union.Gaithersburg,Maryland,July 27–29,1989:305

[11]Standish E M,Williams J G.Orbital Ephemerides of the Sun,Moon,and Planets//Seidelmann P K.Explanatory Supplement to the Astronomical Almanac.Mill Valley,CA:University Science Books,1992:279

[12]Folkner W M,Williams J G,Boggs D H,et al.The Planetary and Lunar Ephemerides DE430 and DE431.Interplanetary Network Progress Report:IPN Progress Report 42-196.US:California Institute of Technology,2014

[13]Malhotra R,Wang X Y.MNRAS,2016,265:4381

[14]Manca F,Sicoli P,Test A.MmSAI,2016,87:72

[15]Wlodarczyk I.Proceedings of the Polish Astronomical Society,2016,3:108

[16]Rivlin T J.The Chebyshev Polynomials.New York:John Wiley&Sons,1974

[17]林成森.數(shù)值計(jì)算方法:下冊(cè).2版.北京:科學(xué)出版社,2005