高　凡　彭碧波　鐘　敏　許厚澤

1　中國(guó)科學(xué)院測(cè)量與地球物理研究所大地測(cè)量與地球動(dòng)力學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，武漢市徐東大街340號(hào)，430077 2　中國(guó)科學(xué)院大學(xué)，北京市玉泉路甲19號(hào)，100049

?

基于帶諧重力場(chǎng)的測(cè)高衛(wèi)星軌道設(shè)計(jì)

高凡1,2彭碧波1鐘敏1許厚澤1

1中國(guó)科學(xué)院測(cè)量與地球物理研究所大地測(cè)量與地球動(dòng)力學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，武漢市徐東大街340號(hào)，430077 2中國(guó)科學(xué)院大學(xué)，北京市玉泉路甲19號(hào)，100049

推導(dǎo)了一組測(cè)高衛(wèi)星軌道設(shè)計(jì)的理論公式，介紹了相關(guān)方法，并編寫了一套軟件。分別以激光測(cè)高衛(wèi)星ICESat和雷達(dá)測(cè)高衛(wèi)星HY-2A為例，基于31階地球帶諧重力場(chǎng)，給出了其平均軌道參數(shù)與密切軌道參數(shù)。軌道模擬顯示，衛(wèi)星經(jīng)過(guò)各自的回歸周期之后，其星下點(diǎn)均能精確回到初始位置，并且長(zhǎng)半軸、軌道傾角、偏心率和近地點(diǎn)幅角均沒有發(fā)生長(zhǎng)期和長(zhǎng)周期的變化，符合參考軌道要求。

測(cè)高衛(wèi)星；重復(fù)地面軌跡；凍結(jié)軌道；HY-2A

衛(wèi)星測(cè)高是一種先進(jìn)的衛(wèi)星大地測(cè)量手段，主要分為雷達(dá)測(cè)高和激光測(cè)高2種模式。在雷達(dá)模式下，其觀測(cè)直徑為1.2～6 km[1]，而對(duì)于激光模式，其地表光斑的直徑僅為大約70 m[2]。它們的星下點(diǎn)需要按所設(shè)計(jì)的軌跡周期性運(yùn)動(dòng)，因此，嚴(yán)格重復(fù)地面軌跡軌道在測(cè)高衛(wèi)星軌道設(shè)計(jì)中顯得尤為重要。另外，為了維持該模式，衛(wèi)星在軌運(yùn)行中長(zhǎng)半軸、偏心率、軌道傾角和近地點(diǎn)幅角這4個(gè)軌道根數(shù)不能發(fā)生長(zhǎng)期與長(zhǎng)周期的變化。

針對(duì)不同的衛(wèi)星任務(wù)，國(guó)際上一些學(xué)者對(duì)此問題作了深入研究。1990年，美國(guó)宇航局(NASA)的Carl[3]針對(duì)GeoSat衛(wèi)星的嚴(yán)格重復(fù)地面軌跡模式，基于地球重力場(chǎng)J2項(xiàng)，提出計(jì)算衛(wèi)星軌道平均長(zhǎng)半軸的方法。1995年，德克薩斯大學(xué)空間研究中心(CSR)的Samsung[4]基于31階地球重力場(chǎng)，系統(tǒng)地對(duì)NASA地球科學(xué)部的激光測(cè)高衛(wèi)星系統(tǒng)(GLAS)的軌道進(jìn)行了設(shè)計(jì)、模擬與分析。1997年，西班牙薩拉戈薩大學(xué)的Antonio等[5]提出直接在地固系下求出周期運(yùn)行軌道的方法。1999年，西班牙皇家天文臺(tái)的Martin[6]在Antonio基礎(chǔ)上，基于較高階的帶諧重力場(chǎng)給出了計(jì)算重復(fù)地面軌跡軌道模式下的方法。2008年，Martin[7]給出基于高階重力場(chǎng)的一種快速設(shè)計(jì)重復(fù)地面軌跡軌道的方法，并對(duì)Topex/Poseidon衛(wèi)星127圈/10 d模式下的軌道進(jìn)行分析。 2008年，CSR的Nadege[8]對(duì)NASA的ICESat衛(wèi)星在3種嚴(yán)密重復(fù)地面軌跡模式下的地面軌跡子周期、軌道機(jī)動(dòng)以及交叉點(diǎn)進(jìn)行了非常詳細(xì)的描述與分析。以上工作主要是來(lái)自美國(guó)的NASA與法國(guó)的CNES[9]，其目的是在滿足工程需求的同時(shí)，減少軌道機(jī)動(dòng)次數(shù)，節(jié)省燃料，延長(zhǎng)衛(wèi)星壽命。我國(guó)第一顆測(cè)高衛(wèi)星——HY-2A作為一顆試驗(yàn)衛(wèi)星已經(jīng)成功在軌運(yùn)行[10,11]，鑒于我國(guó)未來(lái)還有更多的測(cè)高衛(wèi)星任務(wù)[12-14]，積極開展此類衛(wèi)星的軌道設(shè)計(jì)研究具有較大的工程價(jià)值。

本文首先依據(jù)重復(fù)地面軌跡軌道與凍結(jié)軌道的特征，利用拉格朗日行星方程分析平均帶諧重力場(chǎng)對(duì)各軌道根數(shù)的攝動(dòng)，優(yōu)化平均軌道元素確定的公式，使之便于程序計(jì)算。然后，分別對(duì)短周期項(xiàng)中的一階項(xiàng)與余項(xiàng)進(jìn)行修正，并以ICESat與HY-2A為例給出31階帶諧重力場(chǎng)下符合條件的平均軌道元素和密切軌道元素。

1　重復(fù)地面軌跡軌道與凍結(jié)軌道的特征

所謂重復(fù)地面軌跡軌道是指：在慣性系下，衛(wèi)星圍繞地球運(yùn)動(dòng)了N圈，在地固系下，衛(wèi)星的升交點(diǎn)軸(不連續(xù)的)旋轉(zhuǎn)了D圈之后，兩者正好重復(fù)通過(guò)地球上固定的同一子午面，其中N和D是兩個(gè)互質(zhì)的正整數(shù)。

為了使軌道維持上述重復(fù)地面軌跡軌道模式，要求衛(wèi)星在軌運(yùn)行中，其長(zhǎng)半軸、偏心率、軌道傾角和近地點(diǎn)幅角均不產(chǎn)生長(zhǎng)期和長(zhǎng)周期的變化。

2　平均軌道根數(shù)的確定

衛(wèi)星在軌運(yùn)行時(shí)，重力場(chǎng)是最主要的攝動(dòng)源，其對(duì)軌道根數(shù)的攝動(dòng)分為長(zhǎng)期項(xiàng)、長(zhǎng)周期項(xiàng)和短周期項(xiàng)[12-13]。只考慮長(zhǎng)期項(xiàng)和長(zhǎng)周期項(xiàng)影響而得到的軌道根數(shù)，稱為平均軌道根數(shù)。本節(jié)首先介紹分離長(zhǎng)期項(xiàng)和長(zhǎng)周期項(xiàng)的方法，并推導(dǎo)出經(jīng)正則化的平均帶諧重力場(chǎng)表達(dá)式，以便于程序計(jì)算。其次，利用拉格朗日行星方程，給出平均帶諧重力場(chǎng)對(duì)各個(gè)軌道根數(shù)攝動(dòng)的表達(dá)式。最后，提出一種較為簡(jiǎn)潔的確定平均長(zhǎng)半軸的方法，并介紹確定平均偏心率的方法。

2.1平均帶諧重力場(chǎng)對(duì)各軌道根數(shù)的攝動(dòng)分析

平均帶諧重力場(chǎng)的表達(dá)式如下：

(1)

將式(1)代入拉格朗日行星方程[15]可得：

(2)

(3)

(4)

(5)

(6)

(7)

2.2平均軌道根數(shù)的確定

圖1　凍結(jié)偏心率的確定(a=6 971 519.180,i=94)Fig.1　The determination of frozen eccentricity (a=6 971 519.180,i=94)

根據(jù)重復(fù)地面軌跡軌道的特性有：

(8)

(9)

(10)

式(10)兩邊都包含有a，其初值?。?/p>

(11)

表1　ICESat和HY-2A的平均軌道根數(shù)

表中的D、N和i等相關(guān)信息來(lái)自文獻(xiàn)[2]和文獻(xiàn)[17]。

3　密切軌道根數(shù)的確定

§2.2中得到的平均軌道根數(shù)，需要經(jīng)過(guò)短周期項(xiàng)修正轉(zhuǎn)化成密切軌道元素，才能產(chǎn)生實(shí)際應(yīng)用價(jià)值。本節(jié)將介紹修正短周期項(xiàng)的方法。

3.1修正短周期攝動(dòng)項(xiàng)中一階項(xiàng)的影響

密切軌道根數(shù)與平均軌道根數(shù)的關(guān)系如下：

(12)

(13)

利用拉格朗日行星方程，可以求出q*中的各個(gè)分量[4]。

表2給出了ICESat經(jīng)一階短周期項(xiàng)改正的初始軌道和由該初始軌道在31×0的GGM03帶諧重力場(chǎng)下經(jīng)8階Guass-Jackson積分器積分[18]得到的最終軌道。

表2　ICESat的平均軌道根數(shù)+短周期項(xiàng)一階項(xiàng)

表2中的B和L分別表示衛(wèi)星的大地緯度與大地經(jīng)度。從這兩項(xiàng)可以看出，加上一階短周期攝動(dòng)項(xiàng)之后，地面軌跡并不能閉合，在地面上緯度方向和經(jīng)度方向分別會(huì)產(chǎn)生大約-122 332 m 和-14 297 m的距離之差。

3.2修正其他短周期項(xiàng)的影響

從表2可以看出，始末軌道根數(shù)之間的差別主要是在近地點(diǎn)幅角ω與真近點(diǎn)角f 上。為了消除這一差別，需要繼續(xù)調(diào)整初始軌道元素，使之滿足下列方程：

(14)

(15)

式中,ω(q0,tf)、f(q0,tf)是在31×0 GGM03帶諧重力場(chǎng)下由8階Gauss-Jackson 積分器得出，ω0為初始的近地點(diǎn)幅角，f0為初始的真近點(diǎn)角。對(duì)于式(14)和式(15)，只有兩個(gè)方程但是有7個(gè)未知數(shù)，在此只選擇調(diào)整a和e，那么有：

f(a0,e0,tf)-f0=0

式中，a 的單位取地球赤道平均半徑為1。

(16)

式中，δa0、δe0可分別取值為1 m與1×10-6，ff為f(a0,e0,tf)，ωf為ω(a0,e0,tf)，式(16)中偏導(dǎo)數(shù)的計(jì)算方法參見文獻(xiàn)[4]。

表3顯示，兩顆衛(wèi)星經(jīng)過(guò)特定時(shí)間都能精確地回到初始位置，并且初始軌道根數(shù)與最終軌道根數(shù)之差很小。

圖2～5分別給出了HY-2A在193圈/14 d模式下，采樣率為60 s時(shí)密切軌道根數(shù)中的長(zhǎng)半軸、軌道傾角、偏心率和近地點(diǎn)幅角隨時(shí)間變化的情況，白線為其平均軌道元素值。圖2、圖3顯示長(zhǎng)半軸和軌道傾角都在其平均軌道根數(shù)附近小范圍均勻變化，圖4、圖5顯示偏心率和近地點(diǎn)幅角都集中在其平均軌道元素附近變化。以上結(jié)果表明，圖中各個(gè)軌道根數(shù)均沒有產(chǎn)生長(zhǎng)期和長(zhǎng)周期的變化，符合凍結(jié)軌道要求。

表3　ICESat與HY-2A的密切軌道元素

圖2　長(zhǎng)半軸隨時(shí)間變化情況Fig.2　Semi-major history of HY-2A

圖3　軌道傾角隨時(shí)間變化情況Fig.3　Inclination history of HY-2A

圖4　偏心率隨時(shí)間變化情況Fig.4　Eccentricity history of HY-2A

圖5　近地點(diǎn)幅角隨時(shí)間變化情況Fig.5　Argument of perigee history of HY-2A

4　結(jié)　語(yǔ)

測(cè)高衛(wèi)星的軌道設(shè)計(jì)需要滿足兩個(gè)條件：嚴(yán)密重復(fù)地面軌跡軌道和凍結(jié)軌道。圍繞這個(gè)問題，本文系統(tǒng)介紹了基于帶諧重力場(chǎng)的測(cè)高衛(wèi)星軌道設(shè)計(jì)的方法與步驟，并以激光測(cè)高衛(wèi)星ICESat和雷達(dá)測(cè)高衛(wèi)星HY-2A為例對(duì)各個(gè)步驟作詳細(xì)說(shuō)明。

本文主要完成了以下工作：1)推導(dǎo)出經(jīng)正則化的平均帶諧重力場(chǎng)表達(dá)式，并利用拉格朗日行星方程，給出平均帶諧重力場(chǎng)對(duì)各個(gè)軌道根數(shù)攝動(dòng)的表達(dá)式；2)介紹了確定平均偏心率的方法，并提出一種能滿足重復(fù)地面軌跡軌道的條件和確定平均長(zhǎng)半軸的簡(jiǎn)潔方法；3)介紹了修正短周期項(xiàng)攝動(dòng)的方法，并確定了密切軌道參數(shù)；4)在GGM03中的帶諧重力場(chǎng)下，利用Gauss-Jackson積分器，通過(guò)軌道模擬驗(yàn)證了所求得密切軌道符合嚴(yán)密重復(fù)地面軌跡軌道和凍結(jié)軌道的要求。結(jié)果表明，本文的算法以及自主編寫的軟件精確可靠，具有一定的理論價(jià)值和較大的工程應(yīng)用價(jià)值。

下一步的工作是將更多的攝動(dòng)因素納入到軌道設(shè)計(jì)的范圍，如田諧重力場(chǎng)、日月引力、大氣阻力和太陽(yáng)光壓等，以及針對(duì)未來(lái)的衛(wèi)星任務(wù)，積極開展更為細(xì)致的工作，如對(duì)更長(zhǎng)回歸周期的軌道進(jìn)行設(shè)計(jì)，對(duì)其地面軌跡的子周期進(jìn)行分析等。

[1]Seeber G.Satellite Geodesy: Foundations, Methods, and Applications [M].Berlin: Walter de Gruyter,2003

[2]Waleed A, Zwally H J,Rorbit B, et al.The ICESat-2 Laser Altimetry Mission[J]. Proceedings of the IEEE,2010,98(5):735-751

[3]Carl A W.A Prograde Geosat Exact Repeat Mission?[J].Journal of the Astronautical Sciences, 1991,39(1):313-326

[4]Samsung L.Orbit Analysis and Maneuver Design for the Geoscience Laser Altimeter System [D]. Austin: The University of Texas, 1995

[5]Antonio E, Martine L.Frozen Orbits in the Tesseral Artificial Satellite Theory[C]. The 12th International Symposium on Space Flight Dynamics,1997

[6]Martin L.Searching for Repeating Ground Track Orbits: A Systematic Approach [J].The Journal of the Astronautical Sciences,1999,47(3):177-188

[7]Martin L, Ryan P R.Fast Design of Repeat Ground Track Orbits in High-fidelity Geopotentials[J]. The Journal of the Astronautical Sciences, 2008, 56(3):311-324.

[8]Nadege P.Mission Design Concepts for Repeat Groundtrack Orbits and Application to the ICESat Mission[D]. Austin: The University of Texas,2008

[9]Martin L.SADSaM: A Software Assistant for Designing Satellite Missions[R].CNES Report No. DTS/MPI/MS/MN/99-053,1999

[10]張慶君, 張健, 張歡,等. 海洋二號(hào)衛(wèi)星工程研制及在軌運(yùn)行簡(jiǎn)介[J]. 中國(guó)工程科學(xué), 2013(7):12-18(Zhang Qingjun,Zhang Jian,Zhang Huan, et al. The Study of HY-2A Satellite Engineering Development and In-orbit Movement [J]. Engineering Sciences, 2013(7): 12-18)

[11]Gao F,Peng B B,Zhang Y, et al. Analysis of HY-2A Precise Orbit Determination Using DORIS[J]. Journal of Advances in Space Research,2015,55(5):1 394-1 404

[12]Bao L F, Xu H Z, Li Z C. Towards A 1 mGal Accuracy and 1 min Resolution Altimetry Gravity Field [J]. Journal of Geodesy,2013, 87: 961-969

[13]鮑李峰, 許厚澤. 雙星伴飛衛(wèi)星測(cè)高模式及其軌道設(shè)計(jì)[J]. 測(cè)繪學(xué)報(bào), 2014,43 (7): 661-667(Bao Lifeng, Xu Houze.Twin-Satellites Altimetry Mode and Its Orbit Design[J]. Acta Geodaetica et Cartographica Sinica, 2014, 43 (7): 661-667)

[14]李洋,張潤(rùn)寧. 高度計(jì)測(cè)距精度對(duì)沿軌跡重力異常反演的影響[J]. 測(cè)繪學(xué)報(bào)，2015, 44(4): 363-369(Li Yang, Zhang Running．Inference of Altimeter Accuracy on Along-Track Gravity Anomaly Recovery[J]. Acta Geodaetica et Cartographica Sinica,2015, 44(4): 363-369)

[15]Escobal P R. Methods of Orbit Determination [M].Florida: Krieger Publishdhing Company, 1975.

[16]Tapley B, Ries J, Bettadpur S, et al. The GGM03 Mean Earth Gravity Model from GRACE[C]. AGU Fall Meet, 2007

[17]Zhang Y, Vincent T. Input Data for HY-2 Precise Orbit Determination[R]. NAOAS and CNES,2011

[18]Matthew M B, Liam M H. Implementation of Gauss-Jackson Integration for Orbit Propagation[J]. Journal of the Astronautical Sciences, 2004,52(3):331-357

About the first author:GAO Fan, PhD candidate,majors in orbit design and precise orbit determination, E-mail: gaofan02006@126.com.

Altimetry Satellite Orbit Design Based on Zonal Geopotential

GAOFan1,2PENGBibo1ZHONGMin1XUHouze1

1State Key Laboratory of Geodesy and Earth’s Dynamics, Institute of Geodesy and Geophysics, CAS,340 Xudong Street, Wuhan 430077,China 2University of Chinese Academy of Sciences, A19 Yuquan Road, Beijing 100049,China

The orbits of altimetry satellites should be repeat groundtrack and frozen. We derive an algorithm of the orbit design based on zonal geopotential field and develop the software by ourselves. This paper takes ICESat and HY-2A as examples, computing their mean orbit elements and osculation orbit elements. The orbit generation results show that the satellites can get a nice closure after a full cycle and there is no secular and long period change for the semi-major, eccentricity, inclination and the argument of perigee. The orbit can meet the requirement of altimetry satellites and can be taken as a reference orbit.

altimetry satellite; repeat groundtrack; frozen orbit; HY-2A

National Natural Science Foundation of China, No. 41131067.

2015-09-30

高凡，博士生，主要從事衛(wèi)星軌道設(shè)計(jì)與精密定軌研究，E-mail: gaofan02006@126.com。

10.14075/j.jgg.2016.10.004

1671-5942(2016)010-0859-05

P228

A

項(xiàng)目來(lái)源：國(guó)家自然科學(xué)基金(41131067)。