小行星初軌計算的Fortran程序設計

陳姍姍,尹紅星

(山東大學威海分校空間科學與物理學院，山東 威海 264209)

小行星是早期太陽系殘留下來的物質，對太陽系初期的演化研究有重要意義[1]。小行星初軌計算是研究小行星必不可少的環節，初期的有效預報是其發現的保證。目前共發現永久編號小行星21萬多顆。

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

小行星的運動研究是我國的傳統研究課題[7]。國家天文臺BATC項目組1995年5月開始實施大視場CCD小行星巡天計劃(SCAP)[8]，已建立了小行星發現數據庫，小行星巡天數據庫正在建設中。另外紫金山天文臺也開展小行星觀測工作，建立起包括觀測資料處理、初軌計算、攝動計算、軌道改進和位置預報等幾個系統的軟件包。南京大學天文系也有了規模較小的軟件包。

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

文章分3部分。首先簡要給出文中采用的兩種初軌計算方法的原理；然后給出兩種方法的實現及適用范圍；最后給出對Laplace方法流程的調整。

1 初軌計算原理

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

1.1 二體問題

小行星相對太陽的相對運動方程：

(1)

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

1.2 幾何關系

(2)

1.3f、g形式

f、g級數表達式見文獻[11]。

1.4 改進的Laplace方法主要迭代方程

由三次觀測的幾何關系得：

νixi=λizi+pi

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

(4)

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

(5)

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

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

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

(6)

迭代求解[11]即可求得中間時刻的位置速度值。迭代初值取上面f、g級數形式的第1項。

1.5 改進的Gauss方法主要迭代方程

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

λ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及下面各未標注量的意義見文獻[12]。

上方程組消去n1ρ1,n3ρ3再利用f、g關系求出n1,n3的級數表達式，整理得：

(8)

為求得ρ2,r2，由矢量關系得ρ2,r2的另一關系式：

(9)

(10)

(11)

2 程序實現

2.1 時間系統與坐標系統

時間系統選用UT世界時，程序中通過美國海軍天文臺(USNO)的novas程序包[13]轉化為TDB儒略日形式。坐標系統選用J2000赤道坐標系和黃道坐標系。

2.2 運行環境

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

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

2.3 運行準備及要求

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

程序用三次角觀測資料進行初軌計算，運行前建立觀測文件obs，格式與MPC相同。對小行星進行位置預報需建立程序要求格式的pre文件。預報結果放在preorb文件中。

2.4 主要程序流程圖

圖1給出主程序流程圖，對較長時間間隔的三次觀測資料，可使用改進的Laplace和Gauss方法，短時間間隔的觀測資料可使用文中調整的Laplace方法。圖2和圖3為改進的Gauss和Laplace方法的程序流程圖。

圖2 改進的Gauss方法流程圖

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

2.5 主要子程序效果驗證

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

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

表1 初軌計算方法比較.表中模擬觀測數據是MPC小行星歷表服務生成的00014號小行星的相應數據。每組數據給出MPC軌道根數文件中列出的軌道根數及改進的Laplace和Gauss方法得出的軌道根數。其中MPC給出第6個軌道根數為平近點角，程序給出的是過近點時刻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.

續表

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

3 Laplace方法流程的調整

為了使程序實現用短時間間隔的角觀測資料進行軌道預報，對流程相對簡單的Laplace方法進行調整。

3.1 調整原理

圖5 幾何關系驗證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關系驗證

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

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

3.2 流程調整方法

增加調整項：考慮到p、q共6個量，如果每個量進行10次調整，將要進行一百萬次迭代計算過程，一次運行要幾個小時。程序采用折中辦法，只對p做調整，q調整取決于p，以減少程序運行時間。

判斷p、q調整是否良好的條件：

條件1：比較由計算得到的赤經赤緯與觀測值的差別是否足夠小；

條件2：由結果r、v計算出軌道根數，剔除明顯不合適的調整項(如a<0，e<0,e>1)。單用上面兩個條件求出的結果計算軌道根數并進行軌道預報不能得到預期結果；

條件3：加入軌道預報項，挑選預報值與觀測值最接近的調整計算結果。

3.3 調整后的Laplace程序流程及預報效果

圖7 調整后的Laplace流程圖.虛線框內為所加調整部分

程序分別用火星星歷表計算值和山東大學威海天文臺數據進行了測試。表3、表4給出兩測試結果，校天文臺望遠鏡視場為6′，程序可以給出3天內的有效預報值，指導后續觀測。

表3 火星數據測試結果.表中第一行為日期，取該天0.5h、1.0h、1.5h的火星星歷計算值為模擬觀測數據，進行10天的預報，下面對應列給出程序預報值與星歷表計算值在赤經赤緯方向上的差值

表4 校天文臺數據測試結果.表中給出我校天文臺1m口徑望遠鏡的一組觀測數據，給出所列時刻程序預報結果與MPC歷表生成服務相應結果的差值

4 結 論

初軌計算原理清楚，方法很多但根本上可以分成Laplace型和Gauss型。Fortran程序實現了改進的Laplace和Gauss方法，給出其實現流程。對兩程序的測試選用了00014號小行星，該小行星屬于主帶小行星，有代表性地給出兩種方法的使用范圍：初軌計算屬于短弧定軌，但弧段不宜太短，對兩天間隔以上的三次角觀測資料可以進行有效的軌道預報；對更短弧段的觀測資料由于兩種方法計算出的軌道根數偏離太大，不能進行有效的預報以指導觀測。

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

要改善初軌計算結果需要考慮更多因素如進行觀測資料預處理[16]。對多次觀測可以挑選觀測誤差較小的觀測值，或者將統計原理應用于初軌計算，提高初軌計算的精度。

[1] 朱進，高建，關敏，等.小行星的搜尋和定軌[J].云南天文臺臺刊, 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世紀中國學術大典/天文學[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] 劉林，胡松杰，王歆.航天動力學引論[M].南京：南京大學出版社,2006:45-52.

[12] 易照華.天體力學基礎[M].南京：南京大學出版社,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的現狀分析及未來的發展[J].天文學進展,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] 王威，于志堅.航天器軌道確定——模型與算法[M].北京：國防工業出版社,2007. 87-115.