(同濟(jì)大學(xué) 測繪與地理信息學(xué)院，上海 200092)

0 引言

在太陽質(zhì)心慣性系中，八大行星圍繞著太陽做近似圓周運(yùn)動[1]。但是由于地球自身有公轉(zhuǎn)、自轉(zhuǎn)、極移、章動、歲差等運(yùn)動現(xiàn)象，居住在地球上的人類是無法真實地觀察到這樣的運(yùn)動狀態(tài)的。因此，研究太陽系中各天體在地心參考系(包括地心慣性系和地心地固系)中實際的運(yùn)動規(guī)律具有重要的意義。文獻(xiàn)[2-7]只是從不同的層面模擬出個別天體在地心慣性系或地固系中的運(yùn)動規(guī)律，并未利用真實的天體坐標(biāo)對其進(jìn)行研究。本文正是基于DE430計算出天體的實際狀態(tài)向量，并利用C#和SharpGL對不同參考系中太陽系天體的實際運(yùn)動狀態(tài)和運(yùn)行軌跡進(jìn)行研究。

行星/月球(DE/LE)歷表是美國噴氣實驗室(JPL)于20世紀(jì)60年代開始研發(fā)的星歷表。該歷表被廣泛應(yīng)用于星際探測、太空導(dǎo)航、天體測量、引力定律的檢驗、空間大地測量、地球物理等諸多相關(guān)學(xué)科領(lǐng)域[11]。截止目前，JPL網(wǎng)站已經(jīng)發(fā)布了19份不同版本的DE/LE歷表，不同版本的歷表在其內(nèi)容和精度上都有一定的區(qū)別，各版JPL行星歷表的比較可參考JPL網(wǎng)站的介紹以及參考文獻(xiàn)[8-9]。

C#是由C與C++語言派生出來的，在吸收兩者的優(yōu)點(diǎn)的同時，彌補(bǔ)了兩者在界面開發(fā)等方面的缺點(diǎn)。與C和C++相比，C#的語法簡單易學(xué)，配置與調(diào)試簡單和開發(fā)周期短等優(yōu)點(diǎn)，如今被廣泛應(yīng)用于工程開發(fā)。SharpGL是進(jìn)行三維顯示的函數(shù)庫，是一種特殊的OpenGL庫，專門供C#調(diào)用[10]。與OpenGL相比，SharpGL在C#開發(fā)中更具吸引力，它可以通過控件直接將繪圖界面拖拽在窗體上，并進(jìn)一步設(shè)置其屬性事件，省去了一系列在OpenGL中的繁瑣程序。

1 DE星歷表

本文中采用的是DE430歷表，該歷表于2013年4月發(fā)布，包含章動和天平動。起始?xì)v元為1549-12-21(儒略日：2287184.5)，截止歷元為2650-11-25(儒略日：2688976.5)，涵蓋了約1000年的時長。

為了高精度地描述各天體的空間位置，反映其真實的運(yùn)動狀態(tài)，需利用DE歷表計算出天體的位置和速度。那么，如何根據(jù)歷表來獲取行星的狀態(tài)量呢？通?？衫肑PL網(wǎng)站上提供了Fortran版本的子程序進(jìn)行計算，其原理是通過切比雪夫插值多項式[8]插值得到。

Fortran版的PLEPH子程序可以計算出某一時刻目標(biāo)天體相對于參照天體的狀態(tài)量(6個參數(shù)包括位置和速度)。除Fortran版本外，還有其他網(wǎng)站上提供的C版本程序，大部分是由Fortran版本翻譯而成。PLEPH (tdb, npl, nctr, pv)子程序包含4個參數(shù)，其中tdb為儒略日，npl為行星的編號，nctr為中心行星編號(如果計算章動和天平動，中心天體編號為0)，具體行星編號可參見表1，pv為天體的狀態(tài)向量，包含天體的三維坐標(biāo)和速度，單位分別是au和au/day。調(diào)用時輸入4個參數(shù)便可得到相應(yīng)的天體運(yùn)動狀態(tài)向量。

表1 天體(包含章動和天平動)編號

2 坐標(biāo)系統(tǒng)轉(zhuǎn)換

本文除了研究太陽質(zhì)心系中天體的運(yùn)行規(guī)律外，還重點(diǎn)研究天體在地心參考系中的運(yùn)行規(guī)律，因為人類正是處于地心系中觀察各天體運(yùn)行狀態(tài)的。地心系包括地心慣性系和地心地固系，兩者的坐標(biāo)系原點(diǎn)都是地球的質(zhì)量中心，不同的是各自的坐標(biāo)軸指向，因此，實現(xiàn)地心慣性系與地固系轉(zhuǎn)換的關(guān)鍵問題就是求出兩者坐標(biāo)軸指向的旋轉(zhuǎn)矩陣。下面為地固系與地心慣性系的相互轉(zhuǎn)換關(guān)系[12-13]。

考慮在某一歷元時刻t下，協(xié)議地球坐標(biāo)系(CTS)到J2000.0協(xié)議地心慣性系(CIS)的旋轉(zhuǎn)，依次對協(xié)議地球系進(jìn)行極移、地球自轉(zhuǎn)、章動和歲差旋轉(zhuǎn)變換可以轉(zhuǎn)換至J2000.0協(xié)議慣性系。綜上所述，從地固系到J2000.0慣性系的坐標(biāo)轉(zhuǎn)換表達(dá)為：

(1)

由于旋轉(zhuǎn)矩陣為正交矩陣，正交矩陣的逆矩陣是自身的轉(zhuǎn)置矩陣，所以要想實現(xiàn)逆變，即J2000.0慣性系到地固系的轉(zhuǎn)換，可表達(dá)為：

(2)

式中，U為極移矩陣，它將協(xié)議歷元平地球坐標(biāo)系轉(zhuǎn)換至瞬時極地球坐標(biāo)系；S為地球的自轉(zhuǎn)矩陣，它將瞬時極地球坐標(biāo)系轉(zhuǎn)換至真天球坐標(biāo)系；N為章動矩陣，它將真天球坐標(biāo)系轉(zhuǎn)換至t歷元時刻的測瞬平天球坐標(biāo)系；P為歲差矩陣，它將t歷元時刻的測瞬平天球坐標(biāo)系轉(zhuǎn)換到J2000.0地心慣性系。PNSU具體的計算可參考文獻(xiàn)[12-13]。

3 系統(tǒng)設(shè)計與實現(xiàn)

3.1 功能設(shè)計

上文中提到，居住在地球上的人類是無法真實的看到八大行星圍繞著太陽做近似圓周運(yùn)動情景的，那是因為地球自身有公轉(zhuǎn)、自轉(zhuǎn)、極移、章動、歲差等運(yùn)動現(xiàn)象，如果在以地球為參考原點(diǎn)的坐標(biāo)系中，太陽和其他行星的運(yùn)動狀態(tài)是怎樣的呢？

以地球質(zhì)心為原點(diǎn)，如果僅考慮地球公轉(zhuǎn)，可以采用地心慣性系，如果進(jìn)一步考慮地球自轉(zhuǎn)、極移、章動、歲差等現(xiàn)象，可以采用地固系。實際上居住在地球上的人類就是處在地固坐標(biāo)系中觀察各個行星的運(yùn)動狀態(tài)的。為觀察太陽和八大行星的運(yùn)行規(guī)律，本文作者擬開發(fā)設(shè)計一套太陽系空間顯示系統(tǒng)，該系統(tǒng)的主要的設(shè)計功能有：行星的坐標(biāo)計算、天體的顯示、軌道的繪制、天體運(yùn)動的動態(tài)顯示、不同參考系的轉(zhuǎn)換以及觀察角度的變換。圖1為太陽系空間顯示系統(tǒng)設(shè)計開發(fā)界面。

圖1 太陽系空間顯示系統(tǒng)設(shè)計運(yùn)行界面

3.2 結(jié)構(gòu)設(shè)計

上文中已經(jīng)給出了太陽系空間顯示系統(tǒng)的具體設(shè)計功能，那么如何實現(xiàn)這些功能呢？本文將采用逆向思維法一步步進(jìn)行分析，給出詳細(xì)的設(shè)計思路，簡略結(jié)構(gòu)圖見圖2。

圖2 軟件設(shè)計結(jié)構(gòu)圖

1)天體的呈現(xiàn):要在屏幕上呈現(xiàn)出各天體的信息，需要考慮各天體的大小、形狀、表面、位置等信息。首先考慮到太陽和八大行星的體積比數(shù)量級太大，如太陽和地球的體積比約為30萬，為了在屏幕上展示行星，無法按照實際體積比設(shè)定天體的求半徑大小，因此可根據(jù)天體的體積比適當(dāng)確定天體的大?。浑m然各天體的形狀不是球體，但是在屏幕坐標(biāo)系下，各個天體的體積被適當(dāng)縮小，完全可以用球體代替各個天體的形狀；此外，各個天體的表面顏色也有不同，如太陽是橘黃色的、地球是藍(lán)色的，如何將真實的行星表面信息反應(yīng)在模擬的天體上呢？這里可以采用紋理貼圖的方式將不同的行星圖片粘貼在相應(yīng)的球體上；最后考慮如何描述行星的位置，可以采用DE歷表獲取行星的狀態(tài)向量。

2)軌道的繪制:軌道是由線形成的，而軌跡線是由天體的運(yùn)動軌跡點(diǎn)連接形成的，因此只要獲得不同歷元下天體的軌跡點(diǎn)的三維坐標(biāo)，而且這些軌跡點(diǎn)足夠密集，通過點(diǎn)連成線，即可形成天體的運(yùn)行軌道。

3)動畫顯示:以動態(tài)的方式展示天體的運(yùn)動，首先考慮的是如何讓天體自己動起來，如果以一定的時間間隔觸發(fā)一次事件，這個事件包括重畫計算機(jī)屏幕和重新設(shè)定新的天體坐標(biāo)，那么設(shè)置適當(dāng)?shù)臅r間間隔后觸發(fā)該事件，天體的位置不同，便實現(xiàn)了天體的運(yùn)動。C#中的定時器(timer)便有這樣的功能，它可以設(shè)定一定的時間間隔然后不斷觸發(fā)事件[14]。

4)坐標(biāo)系轉(zhuǎn)換:坐標(biāo)系轉(zhuǎn)換主要涉及太陽質(zhì)心系、地球質(zhì)心系和地固系的轉(zhuǎn)換。太陽質(zhì)心系中的天體坐標(biāo)可利用DE歷表得到。地球質(zhì)心系中的天體坐標(biāo)可通過簡單的坐標(biāo)系平移轉(zhuǎn)換得到，即利用天體在太陽質(zhì)心系中的坐標(biāo)減去中心天體在太陽質(zhì)心系中坐標(biāo)得到。地固系中的天體坐標(biāo)需要考慮地球自轉(zhuǎn)、極移、章動和歲差，構(gòu)造旋轉(zhuǎn)矩陣可地心慣性系轉(zhuǎn)換到地固系。

5)視角變換:不同的方向和角度下觀察到天體的運(yùn)動狀態(tài)也是不同的，如從Z軸上方沿下方看，看到的是XOY平面的天體的運(yùn)行狀態(tài)。實現(xiàn)此功能首先要設(shè)定一組眼睛、目標(biāo)和上方向的向量，利用此向量通過矩陣變換可以得到不同視角下的天體運(yùn)動狀態(tài)。

6)其他功能:除上述主要功能外，為提高軟件的人機(jī)交互性能，軟件設(shè)計功能還包括天體演示過程暫停與播放功能、行星是否顯示功能、軌跡是否顯示功能以及幫助等其他功能。

3.4 系統(tǒng)實現(xiàn)

1)定義繪制行星的方法，可以設(shè)置行星的紋理、坐標(biāo)、旋轉(zhuǎn)角度。

public void drawplanet(string bmppath,floatsize,int number, float rot, float x,floaty,float z)

{

OpenGL gl = openGLControl.OpenGL;

gl.LoadIdentity(); //加載單位矩陣

gl.Translate(x, y, z); //行星平移坐標(biāo)

gl.Rotate(rotation[number], 0.0f, 0.0f, 1.0f); //饒Z軸旋轉(zhuǎn)

texture.Create(gl,bmppath); //設(shè)置紋理

gl.Sphere(ptrBall, size, 100, 100); //繪制星球

rotation[number] += rot;

}

2)繪制行星軌道，每天畫一個點(diǎn)，將其連接起來。Pleph(jd,target, center, R)方法可以得到某一歷元下目標(biāo)行星相對于參照行星的相對狀態(tài)量。

gl.Begin(OpenGL.GL_LINE_STRIP); //開始繪制線

for (int i=0; j < 365 * jplyear; i++)

{

jd = 2415020.5 + i; //改變?nèi)迓匀?/p>

ClassPLEPH.Pleph(jd,target, center, R); //獲得天體的相對坐標(biāo)

gl.Vertex(R[0], R[1], R[2]); //繪制行星的位置點(diǎn)

}

gl.End(); //結(jié)束繪制

3)設(shè)置定時器，定時器每觸發(fā)一次，儒略日增加一日，便重新繪制一次行星，行星的坐標(biāo)與上次不同，實現(xiàn)了行星的運(yùn)動。

private void timer1_Tick(object sender, EventArgs e)

{

double jd;

OpenGL gl = openGLControl.OpenGL;

gl.Clear(OpenGL.GL_COLOR_BUFFER_BIT | OpenGL.GL_DEPTH_BUFFER_BIT);

//清除顏色緩沖以及深度緩沖

jd = 2415080.5 + ii ; //儒略日變化

ClassPLEPH.Pleph(jd, 1, 12, PO1); //新儒略日的行星相對坐標(biāo)

drawplanet(bmppath[0], 0.1f, 0, 0.103f, (float)PO1[0], (float)PO1[1], (float)PO1[2]);

//繪制行星

ii++

}

4)如要得到地心慣性系的天體坐標(biāo)，通過Pleph(jd,target, center, R)方法可將地球質(zhì)心設(shè)置為中心參照體，從而得到天體在地心慣性系中的坐標(biāo)，其原理是將各天體質(zhì)心在太陽質(zhì)心系中的狀態(tài)量減去地球質(zhì)心在太陽質(zhì)心系中的狀態(tài)量。把地心慣性系的坐標(biāo)轉(zhuǎn)換到地固系中，需要進(jìn)行極移、地球自轉(zhuǎn)、章動和歲差旋轉(zhuǎn)變換，作者定義的Frametrans(year,month,day,h,m,s,style, x0,y0,z0, x,y, z)方法可以實現(xiàn)地固系和地心系的相互轉(zhuǎn)換。

圖3為天體在太陽質(zhì)心系、地心系和地固系中的運(yùn)動狀態(tài)，由圖可知各天體在不同參考系下的運(yùn)行狀態(tài)和運(yùn)行軌跡是完全不同的。在太陽質(zhì)心慣性系中，八大行星圍繞太陽質(zhì)心做近似圓周運(yùn)動；而在地心慣性系中，由于地球的公轉(zhuǎn)現(xiàn)象，各天體(除太陽外)圍繞地球質(zhì)心運(yùn)動的同時還伴有回旋現(xiàn)象發(fā)生，且天體離地球的距離越近，回旋現(xiàn)象越明顯；在地心地固系中，正如人類在地球真實觀察到的樣子，各天體圍繞地心做螺旋曲線運(yùn)動，例如太陽以地球為參照中心，人類可觀察到太陽每天東升西落的現(xiàn)象。

圖3 不同參考系下天體的運(yùn)動

5)獲取鼠標(biāo)的屏幕坐標(biāo)，將其傳遞至函數(shù)gl.LookAt(eye[0], eye[1], eye[2], center[0], center[1], center[2], up[0], up[1], up[2])中作為參數(shù)值，從而實現(xiàn)視角變化，函數(shù)中的三組參數(shù)分別代表眼睛的位置、瞄準(zhǔn)的位置和向上的方向。該函數(shù)實際上定義了一個視圖矩陣，與當(dāng)前矩陣相乘實現(xiàn)視點(diǎn)轉(zhuǎn)換。

private void openGLControl_MouseDown(object sender, MouseEventArgs e)

{

//拾取鼠標(biāo)的屏幕坐標(biāo)

center[1] = e.X;

center[2] = e.Y;

}

除了利用鼠標(biāo)控制視點(diǎn)變換外，也可以通過鍵盤按鍵改變gl.LookAt的其他參數(shù)值實現(xiàn)場景的旋轉(zhuǎn)和平移等功能。圖4為不同視角下太陽系中天體的運(yùn)動狀態(tài)：

圖4 不同視點(diǎn)下天體的運(yùn)動

3.5 系統(tǒng)優(yōu)化

為了觀察較長時間天體運(yùn)動軌跡的變化情況，在繪制軌道的過程中，如設(shè)置繪制20年的軌道，一天繪制一個天體位置點(diǎn)，繪制八條軌道大約需要58 000個點(diǎn)，因此每次屏幕重畫過程中，軌道的繪制占用了較長的時間，導(dǎo)致程序運(yùn)行效率不高。同時還應(yīng)考慮到天體的運(yùn)動具有周期性，每天相同的時刻繪制天體的位置點(diǎn)可能會導(dǎo)致周期點(diǎn)位重合，造成一定的資源重復(fù)和浪費(fèi)。

在此情況下，考慮使用切比雪夫插節(jié)點(diǎn)，摒棄等距節(jié)點(diǎn)。使用切比雪夫插值點(diǎn)可以保留更多更有效的軌道信息，更加真實地反映出行星的運(yùn)動軌跡，避免了資源的重復(fù)浪費(fèi)。在任意區(qū)間[]上的n+1個切比雪夫節(jié)點(diǎn)可通過線性變換得到[15]:

(3)

4 結(jié)語

本文首先基于DE430得到太陽系各天體的空間狀態(tài)量，然后利用C#和SharpGL混合開發(fā)出太陽系三維空間顯示系統(tǒng)，動態(tài)顯示太陽和八大行星在不同參考系中的運(yùn)動狀態(tài)和運(yùn)行軌跡。顯示結(jié)果表明，不同參考系下天體的運(yùn)行狀態(tài)和軌跡是不同的，在研究天體時若選擇太陽質(zhì)心系，各天體的圍繞太陽做近似圓周運(yùn)動；若選擇地心慣性系，各天體圍繞地球運(yùn)動，除太陽外各天體還發(fā)生回旋現(xiàn)象，這是由于地球的公轉(zhuǎn)產(chǎn)生的；若選擇地固系，地球為參考原點(diǎn)且保持不動，各天體沿螺旋曲線繞地球進(jìn)行運(yùn)轉(zhuǎn)。由此可見，選擇不同的參考系對于研究目的有著極其重要的意義。