小行星初軌計(jì)算的Fortran程序設(shè)計(jì)

陳姍姍,尹紅星

(山東大學(xué)威海分?？臻g科學(xué)與物理學(xué)院，山東 威海 264209)

小行星是早期太陽系殘留下來的物質(zhì)，對(duì)太陽系初期的演化研究有重要意義[1]。小行星初軌計(jì)算是研究小行星必不可少的環(huán)節(jié)，初期的有效預(yù)報(bào)是其發(fā)現(xiàn)的保證。目前共發(fā)現(xiàn)永久編號(hào)小行星21萬多顆。

國際小行星中心(MPC)[2]對(duì)小行星進(jìn)行編號(hào)確定工作，為方便觀測(cè)者觀測(cè)，開通了免費(fèi)的網(wǎng)上服務(wù)如小行星軌道預(yù)報(bào)、軌道根數(shù)下載等。意大利幾位科學(xué)家在網(wǎng)上提供的小行星軌道計(jì)算免費(fèi)軟件包OrbFit[3]是現(xiàn)在較常用的軟件包。國外有幾個(gè)網(wǎng)站[1]可以提供小行星歷表服務(wù)或軌道根數(shù)下載服務(wù)，如意大利Pisa大學(xué)的空間力學(xué)組在AstDys小行星動(dòng)力學(xué)網(wǎng)站[4]給出了永久編號(hào)小行星、有多次觀測(cè)記錄和單次觀測(cè)記錄小行星的軌道根數(shù)和其他固有根數(shù)的下載，美國洛威爾天文臺(tái)的Edward Bowell博士在網(wǎng)上免費(fèi)提供完整且每天更新的小行星高精度吻切軌道根數(shù)[5]，美國宇航局噴氣推進(jìn)實(shí)驗(yàn)室(JPL)免費(fèi)提供的HORIZONS[6]太陽系數(shù)據(jù)和歷表服務(wù)，可以通過遠(yuǎn)程登錄、電子郵件及在線使用方式獲取包括小行星在內(nèi)的太陽系內(nèi)天體的各種數(shù)據(jù)資源，包括軌道根數(shù)及星歷表計(jì)算。

小行星的運(yùn)動(dòng)研究是我國的傳統(tǒng)研究課題[7]。國家天文臺(tái)BATC項(xiàng)目組1995年5月開始實(shí)施大視場(chǎng)CCD小行星巡天計(jì)劃(SCAP)[8]，已建立了小行星發(fā)現(xiàn)數(shù)據(jù)庫，小行星巡天數(shù)據(jù)庫正在建設(shè)中。另外紫金山天文臺(tái)也開展小行星觀測(cè)工作，建立起包括觀測(cè)資料處理、初軌計(jì)算、攝動(dòng)計(jì)算、軌道改進(jìn)和位置預(yù)報(bào)等幾個(gè)系統(tǒng)的軟件包。南京大學(xué)天文系也有了規(guī)模較小的軟件包。

小行星的定軌工作離不開太陽系歷表的使用。自1984年起國際上普遍采用JPL提供的高精度歷表。紫金山天文臺(tái)現(xiàn)已在其網(wǎng)站[9]上提供電子天文歷表服務(wù)。我校于07年建成1m口徑的赤道式反射望遠(yuǎn)鏡[10]，并利用該望遠(yuǎn)鏡進(jìn)行小行星搜尋工作。

文章分3部分。首先簡要給出文中采用的兩種初軌計(jì)算方法的原理；然后給出兩種方法的實(shí)現(xiàn)及適用范圍；最后給出對(duì)Laplace方法流程的調(diào)整。

1 初軌計(jì)算原理

初軌計(jì)算方法很多，從根本上可以分為Laplace型和Gauss型[11]，Laplace型先設(shè)法從三次角觀測(cè)量求解某時(shí)刻的位置速度矢量，Gauss型則從三次角觀測(cè)量求解兩個(gè)時(shí)刻的位置矢量。下面簡要介紹文中使用的改進(jìn)的Laplace和Gauss方法的計(jì)算原理。

1.1 二體問題

小行星相對(duì)太陽的相對(duì)運(yùn)動(dòng)方程：

(1)

小行星和太陽組成的系統(tǒng)所受合外力為零、合力矩為零、系統(tǒng)內(nèi)只有保守力作用，分別應(yīng)用動(dòng)量守恒定律、角動(dòng)量守恒定律、能量守恒定律引入這6個(gè)常數(shù)，即為軌道根數(shù)。軌道根數(shù)能反映小行星的運(yùn)動(dòng)特性，通常選擇為軌道半長軸a、偏心率e、軌道傾角i、升交點(diǎn)黃經(jīng)Ω、近地(星)點(diǎn)幅角ω、平近點(diǎn)角M(近地點(diǎn)時(shí)刻τ)。對(duì)小行星運(yùn)動(dòng)的橢圓軌道，軌道根數(shù)具有明顯的幾何意義，其中a、e為描述軌道大小和形狀的參數(shù)，i、Ω、ω為描述軌道位置的參數(shù)，M(τ)給出小行星的具體位置。

1.2 幾何關(guān)系

(2)

1.3f、g形式

f、g級(jí)數(shù)表達(dá)式見文獻(xiàn)[11]。

1.4 改進(jìn)的Laplace方法主要迭代方程

由三次觀測(cè)的幾何關(guān)系得：

νixi=λizi+pi

νiyi=μizi+qi(i=1,2,3)

(4)

其中，pi=νiXi-λiZi；qi=νiYi-μiZi

(5)

式中λ，μ，ν及下面未注釋的量的意義見文獻(xiàn)[11]。將(3)式代入(4)式，整理得：

νifix-λifiz+νigix′-λigiz′=pi

νifiy-μifiz+νigiy′-μigiz′=qi(i=1,2,3)

(6)

迭代求解[11]即可求得中間時(shí)刻的位置速度值。迭代初值取上面f、g級(jí)數(shù)形式的第1項(xiàng)。

1.5 改進(jìn)的Gauss方法主要迭代方程

小行星在二體問題模型下的軌道是過日心的固定平面，并利用幾何關(guān)系整理得：

λ1n1ρ1-λ2ρ2+λ3n3ρ3=n1X1-X2+n3X3

μ1n1ρ1-μ2ρ2+μ3n3ρ3=n1Y1-Y2+n3Y3

ν1n1ρ1-ν2ρ2+ν3n3ρ3=n1Z1-Z2+n3Z3

(7)

n1、n3及下面各未標(biāo)注量的意義見文獻(xiàn)[12]。

上方程組消去n1ρ1,n3ρ3再利用f、g關(guān)系求出n1,n3的級(jí)數(shù)表達(dá)式，整理得：

(8)

為求得ρ2,r2，由矢量關(guān)系得ρ2,r2的另一關(guān)系式：

(9)

(10)

(11)

2 程序?qū)崿F(xiàn)

2.1 時(shí)間系統(tǒng)與坐標(biāo)系統(tǒng)

時(shí)間系統(tǒng)選用UT世界時(shí)，程序中通過美國海軍天文臺(tái)(USNO)的novas程序包[13]轉(zhuǎn)化為TDB儒略日形式。坐標(biāo)系統(tǒng)選用J2000赤道坐標(biāo)系和黃道坐標(biāo)系。

2.2 運(yùn)行環(huán)境

為統(tǒng)一形式，程序編寫使用f77版本，在linux系統(tǒng)的fedora(2.6.231-42.fcc8)下編譯運(yùn)行，編譯器為系統(tǒng)自帶的gcc(4.1.2)。

圖1 主程序流程圖Fig.1 Data flow in the main program

2.3 運(yùn)行準(zhǔn)備及要求

程序使用的是DE405星歷表[14]，該星歷表基于國際天文協(xié)會(huì)(IAU)所建立的ICRF[15]參考架。星歷表的使用借用了USNO的novas程序包和JPL提供的星歷使用程序selcon.f。使用前需安裝相應(yīng)的星歷文件及星歷表使用程序包。

程序用三次角觀測(cè)資料進(jìn)行初軌計(jì)算，運(yùn)行前建立觀測(cè)文件obs，格式與MPC相同。對(duì)小行星進(jìn)行位置預(yù)報(bào)需建立程序要求格式的pre文件。預(yù)報(bào)結(jié)果放在preorb文件中。

2.4 主要程序流程圖

圖1給出主程序流程圖，對(duì)較長時(shí)間間隔的三次觀測(cè)資料，可使用改進(jìn)的Laplace和Gauss方法，短時(shí)間間隔的觀測(cè)資料可使用文中調(diào)整的Laplace方法。圖2和圖3為改進(jìn)的Gauss和Laplace方法的程序流程圖。

圖2 改進(jìn)的Gauss方法流程圖

圖3 改進(jìn)的Laplace方法流程圖，iter,imax,ft,hug為迭代控制量

2.5 主要子程序效果驗(yàn)證

首先對(duì)由某時(shí)刻位置速度矢量計(jì)算軌道根數(shù)及星歷表計(jì)算程序進(jìn)行了聯(lián)合測(cè)試，圖4為測(cè)試結(jié)果。由于初軌計(jì)算采用的是二體模型，而實(shí)際上在攝動(dòng)力作用下小行星的真實(shí)軌道是一族吻切軌道的包絡(luò)線，觀測(cè)點(diǎn)計(jì)算出的軌道與真實(shí)軌道相切，因此如圖出現(xiàn)觀測(cè)點(diǎn)附近精度高，距離半周期時(shí)偏離很大的現(xiàn)象，另外采用數(shù)據(jù)點(diǎn)少也是造成計(jì)算軌道與真實(shí)軌道偏離的原因。圖示兩方向最大偏移量小于2.0′,可以認(rèn)為兩個(gè)主要子程序效果良好。

對(duì)改進(jìn)的Laplace和Gauss方法進(jìn)行驗(yàn)證時(shí)采用了00014號(hào)小行星不同時(shí)間間隔的模擬觀測(cè)數(shù)據(jù)，驗(yàn)證結(jié)果見表1。由表中數(shù)據(jù)可見，對(duì)于不同時(shí)間間隔的角觀測(cè)數(shù)據(jù)，兩種方法計(jì)算效果不同，都有隨時(shí)間間隔縮小變壞的趨勢(shì)。即改進(jìn)的Laplace方法和Gauss方法對(duì)小行星進(jìn)行初軌計(jì)算時(shí)，只適用于弧段不是很短的情況，此時(shí)兩種方法效果接近，由于程序中Gauss方法考慮了光行差，且Laplace方法f、g階數(shù)只取了四階，高斯方法稍優(yōu)于拉普拉斯方法；對(duì)短時(shí)間隔的觀測(cè)數(shù)據(jù)(小時(shí)量級(jí)間隔)兩種方法計(jì)算的軌道根數(shù)都有較大偏離，因此無法進(jìn)行有效的軌道預(yù)報(bào)。

表1 初軌計(jì)算方法比較.表中模擬觀測(cè)數(shù)據(jù)是MPC小行星歷表服務(wù)生成的00014號(hào)小行星的相應(yīng)數(shù)據(jù)。每組數(shù)據(jù)給出MPC軌道根數(shù)文件中列出的軌道根數(shù)及改進(jìn)的Laplace和Gauss方法得出的軌道根數(shù)。其中MPC給出第6個(gè)軌道根數(shù)為平近點(diǎn)角，程序給出的是過近點(diǎn)時(shí)刻Table 1 Comparison of initial-orbit calculation methods. The simulated observation data are of Asteroid No.00014 generated with the Minor Planet & Comet Ephemeris Service(MPC) online. The table gives the orbital elements from Laplace method,Gaussian method,and MPC orbit database. The sixth orbit element from the MPC is the mean anomaly while the program gives the mean anomaly time.

續(xù)表

模擬觀測(cè)數(shù)據(jù)K08Y31HC200401010000004530570+2153450K08Y31HC200401030000004512180+2158080K08Y31HC200401050000004494380+2202340197V195V191VD39D39D39aeinodwm0(j2000)/t0mpc25854908016756029106908645520963059834601524(K08BU)gau264289220213260582886087052418385090245170953345543lap264266890213167682911287050568389406245170990843509模擬觀測(cè)數(shù)據(jù)K08Y31HC200401010000004530570+2153450K08Y31HC200401020000004521310+2155560K08Y31HC200401030000004512180+2158080197V195V191VD39D39D39aeinodwm0(j2000)/t0mpc25854908016756029106908645520963059834601524(K08BU)gau318423930400148880783487123627660508245114841533498lap318364340400061480756487125957659389245114904540179模擬觀測(cè)數(shù)據(jù)K08Y31HC200401010000004530570+2153450K08Y31HC200401010416704530350+2153500K08Y31HC200401010833304530120+2153560197V195V191VD39D39D39aeinodwm0(j2000)/t0mpc25854908016756029106908645520963059834601524(K08BU)gau153949250418554502620210157304489756245234382144170lap101329110032735500107111769518977016245265899135827

3 Laplace方法流程的調(diào)整

為了使程序?qū)崿F(xiàn)用短時(shí)間間隔的角觀測(cè)資料進(jìn)行軌道預(yù)報(bào)，對(duì)流程相對(duì)簡單的Laplace方法進(jìn)行調(diào)整。

3.1 調(diào)整原理

圖5 幾何關(guān)系驗(yàn)證Fig.5 Test chart for geometric positions.The vertical coordinates are for the differences of the three components of e+-s，respectively.The comparisons are made to the data from DE405 covering 3000 days starting from March 1,2009

圖6 p、q關(guān)系驗(yàn)證

由圖5可見，三方向的差值最大為10e-14量級(jí)，可以認(rèn)為幾何關(guān)系是成立的。但圖6中p、q的差值在10e-7量級(jí)，相比幾何關(guān)系的差值很大，對(duì)p、q對(duì)求解的影響進(jìn)行了測(cè)試(表2)。由表中數(shù)據(jù)可見，p、q在此量級(jí)上的線性變化對(duì)結(jié)果影響不大，非線性變化使求解值有較大波動(dòng)，方程組對(duì)p、q變化敏感。p、q關(guān)系本身的差別造成了求解結(jié)果的較大偏離，因而直接求解不能得出很好的結(jié)果，為此求解之前需要對(duì)p、q進(jìn)行調(diào)整。

表2 p、q變化對(duì)結(jié)果影響.表中求解方程為Laplace方法的主要方程(6)式，

3.2 流程調(diào)整方法

增加調(diào)整項(xiàng)：考慮到p、q共6個(gè)量，如果每個(gè)量進(jìn)行10次調(diào)整，將要進(jìn)行一百萬次迭代計(jì)算過程，一次運(yùn)行要幾個(gè)小時(shí)。程序采用折中辦法，只對(duì)p做調(diào)整，q調(diào)整取決于p，以減少程序運(yùn)行時(shí)間。

判斷p、q調(diào)整是否良好的條件：

條件1：比較由計(jì)算得到的赤經(jīng)赤緯與觀測(cè)值的差別是否足夠小；

條件2：由結(jié)果r、v計(jì)算出軌道根數(shù)，剔除明顯不合適的調(diào)整項(xiàng)(如a<0，e<0,e>1)。單用上面兩個(gè)條件求出的結(jié)果計(jì)算軌道根數(shù)并進(jìn)行軌道預(yù)報(bào)不能得到預(yù)期結(jié)果；

條件3：加入軌道預(yù)報(bào)項(xiàng)，挑選預(yù)報(bào)值與觀測(cè)值最接近的調(diào)整計(jì)算結(jié)果。

3.3 調(diào)整后的Laplace程序流程及預(yù)報(bào)效果

圖7 調(diào)整后的Laplace流程圖.虛線框內(nèi)為所加調(diào)整部分

程序分別用火星星歷表計(jì)算值和山東大學(xué)威海天文臺(tái)數(shù)據(jù)進(jìn)行了測(cè)試。表3、表4給出兩測(cè)試結(jié)果，校天文臺(tái)望遠(yuǎn)鏡視場(chǎng)為6′，程序可以給出3天內(nèi)的有效預(yù)報(bào)值，指導(dǎo)后續(xù)觀測(cè)。

表3 火星數(shù)據(jù)測(cè)試結(jié)果.表中第一行為日期，取該天0.5h、1.0h、1.5h的火星星歷計(jì)算值為模擬觀測(cè)數(shù)據(jù)，進(jìn)行10天的預(yù)報(bào)，下面對(duì)應(yīng)列給出程序預(yù)報(bào)值與星歷表計(jì)算值在赤經(jīng)赤緯方向上的差值

表4 校天文臺(tái)數(shù)據(jù)測(cè)試結(jié)果.表中給出我校天文臺(tái)1m口徑望遠(yuǎn)鏡的一組觀測(cè)數(shù)據(jù)，給出所列時(shí)刻程序預(yù)報(bào)結(jié)果與MPC歷表生成服務(wù)相應(yīng)結(jié)果的差值

4 結(jié) 論

初軌計(jì)算原理清楚，方法很多但根本上可以分成Laplace型和Gauss型。Fortran程序?qū)崿F(xiàn)了改進(jìn)的Laplace和Gauss方法，給出其實(shí)現(xiàn)流程。對(duì)兩程序的測(cè)試選用了00014號(hào)小行星，該小行星屬于主帶小行星，有代表性地給出兩種方法的使用范圍：初軌計(jì)算屬于短弧定軌，但弧段不宜太短，對(duì)兩天間隔以上的三次角觀測(cè)資料可以進(jìn)行有效的軌道預(yù)報(bào)；對(duì)更短弧段的觀測(cè)資料由于兩種方法計(jì)算出的軌道根數(shù)偏離太大，不能進(jìn)行有效的預(yù)報(bào)以指導(dǎo)觀測(cè)。

為了可以利用一天內(nèi)的三次角觀測(cè)數(shù)據(jù)進(jìn)行軌道預(yù)報(bào)指導(dǎo)后續(xù)觀測(cè)，對(duì)流程較簡單的Laplace方法進(jìn)行了調(diào)整。通過對(duì)該方法主要原理的驗(yàn)證，找出問題的所在： Laplace方法主要迭代方程組對(duì)方程組右矩陣的變化較敏感，而方程右矩陣本身的不確定度導(dǎo)致無法求得理想的位置速度矢量，因而相應(yīng)軌道根數(shù)有較大偏離，也就無法進(jìn)行有效的軌道預(yù)報(bào)。因此程序考慮增加主要迭代方程右矩陣的調(diào)整項(xiàng)和反饋項(xiàng)，經(jīng)調(diào)整后可以指導(dǎo)我校天文臺(tái)的小行星觀測(cè)，初步實(shí)現(xiàn)預(yù)報(bào)功能。

要改善初軌計(jì)算結(jié)果需要考慮更多因素如進(jìn)行觀測(cè)資料預(yù)處理[16]。對(duì)多次觀測(cè)可以挑選觀測(cè)誤差較小的觀測(cè)值，或者將統(tǒng)計(jì)原理應(yīng)用于初軌計(jì)算，提高初軌計(jì)算的精度。

[1] 朱進(jìn)，高建，關(guān)敏，等.小行星的搜尋和定軌[J].云南天文臺(tái)臺(tái)刊, 2002(3):17-20.

ZHU Jin,GAO Jian,GUAN Min,et al.Asteroid Searching and Orbit Determination[J].Publications of Yunnan Observatory,2002(3):17-20.

[2] http://www.cfa.harvard.edu/iau/MPC.html.

[3] http://adams.dm.unipi.it/～orbmaint/orbfit/.

[4] http://hamilton.dm.unipi.it/astdys/index.php?pc=4.

[5] ftp://ftp.lowell.edu/pub/elgb/astorb.html.

[6] http://ssd.jpl.nasa.gov/?horizons.

[7] 王綬琯.20世紀(jì)中國學(xué)術(shù)大典/天文學(xué)[M].福州：福建教育出版社，2002:75-77.

[8] http://batc.bao.ac.cn/work/jobs/work/work10.htm.

[9] http://159.226.71.170/dianzili.htm.

[10] http://astro.wh.sdu.edu.cn/apparatus/ShowArticle.asp?ArticleID=4.

[11] 劉林，胡松杰，王歆.航天動(dòng)力學(xué)引論[M].南京：南京大學(xué)出版社,2006:45-52.

[12] 易照華.天體力學(xué)基礎(chǔ)[M].南京：南京大學(xué)出版社,1993: 96-120.

[13] http://aa.usno.navy.mil/software/novas/novas_f/novasf_intro.html.

[14] ftp://ssd.jpl.nasa.gov/pub/eph/planets/usrguide.

[15] 喬書波，李金嶺，孫付平，等. ICRF的現(xiàn)狀分析及未來的發(fā)展[J].天文學(xué)進(jìn)展,2007(2):147-160.

QIAO Shu-bo,LI Jin-ling,SUN Fu-ping,et al.Analysis of Present Status of ZCRF and ZCRF in Future[J].Progress in Astronmy,2007(2):149-160.

[16] 王威，于志堅(jiān).航天器軌道確定——模型與算法[M].北京：國防工業(yè)出版社,2007. 87-115.