黃昆學，常曉濤，孫　玉

(1. 武漢大學測繪學院，湖北 武漢 430079； 2. 國家測繪地理信息局衛(wèi)星測繪應(yīng)用中心，北京 100830；

3. 代爾夫特理工大學，荷蘭 代爾夫特)

Using Lunar Prospector Doppler Data to Recover and Analyze

Lunar Gravity Field

HUANG Kunxue，CHANG Xiaotao，SUN Yu

?

利用Lunar Prospector Doppler數(shù)據(jù)解算并分析月球重力場

黃昆學1，常曉濤2，孫玉3

(1. 武漢大學測繪學院，湖北 武漢 430079； 2. 國家測繪地理信息局衛(wèi)星測繪應(yīng)用中心，北京 100830；

3. 代爾夫特理工大學，荷蘭 代爾夫特)

Using Lunar Prospector Doppler Data to Recover and Analyze

Lunar Gravity Field

HUANG Kunxue，CHANG Xiaotao，SUN Yu

摘要：Doppler技術(shù)可高精度測定衛(wèi)星之間的距離變化率，可以和其他技術(shù)聯(lián)合用于月球衛(wèi)星定軌。本文利用Lunar Prospector月球項目4個月的Doppler跟蹤數(shù)據(jù)，解算了75階次月球重力場模型。首先介紹了數(shù)據(jù)預(yù)處理方法、力模型的使用及參數(shù)設(shè)置等關(guān)鍵技術(shù)，然后把解算的重力場模型與同階次的其他月球重力場模型進行比較，可以看出低階位系數(shù)較LP75G月球重力場模型稍差，但高階位系數(shù)精度有所提高，可能是僅利用了低軌單星數(shù)據(jù)的原因。

關(guān)鍵詞：Lunar Prospector；月球重力場；Doppler；數(shù)據(jù)解算

Lunar Prospector(LP)是美國NASA的月球探索任務(wù)，衛(wèi)星于1998年1月6日發(fā)射升空，經(jīng)過一系列調(diào)整之后進入月球極地近圓軌道。LP繞月軌道被設(shè)計為過兩極的圓軌道,軌道高度設(shè)計為100 km,運行速度約為5517 km/h,運行周期為118 min。LP的設(shè)計初衷是用于探測、研究月球全球構(gòu)成等，重點是探測南極水冰和分析月球的重力場及磁場。為完成不同的科學任務(wù),LP任務(wù)分成3個階段。正常任務(wù)階段從衛(wèi)星進入月球軌道至1999年1月，歷時1年。其軌道為極地圓軌道，平均高度為100 km[1]。之后，LP進入短暫的過渡階段，平均軌道高度為40 km。過渡階段的主要任務(wù)是驗證軌道預(yù)測的能力，為其后的擴展任務(wù)階段做準備。LP擴展任務(wù)階段開始于1999年1月29日，直到7月31日撞擊月球表面為止，平均軌道高度約為30 km。采用較低的環(huán)月軌道是為了得到高精度、高分辨率的月球重力場模型。

美國JPL已經(jīng)利用LP數(shù)據(jù)研制了一系列高精度、高階次月球重力場模型，包括LP75D、LP75G、LP100J、LP100K、LP150和LP165等[2-3]。此外，也有學者利用GEODYN和SOLVE軟件處理了公布在PDS上的LP Doppler跟蹤數(shù)據(jù)，從而實踐了月球重力場解算技術(shù)流程。文獻[4]利用正常任務(wù)階段數(shù)據(jù)模制了75階月球重力場模型。也有學者利用擴展任務(wù)階段數(shù)據(jù)，以LP100J為初始重力場模型解算了100階月球重力場模型。但上述文獻主要介紹了基本解算原理和方法，均未涉及LP數(shù)據(jù)預(yù)處理等具體信息。本文詳盡介紹了數(shù)據(jù)預(yù)處理方法、力模型的使用及參數(shù)設(shè)置等關(guān)鍵技術(shù)，然后對解算的重力場模型進行了比較分析。

一、數(shù)據(jù)處理方法和流程

本文所采用的數(shù)據(jù)來自美國PDS的LP軌道數(shù)據(jù)文件(ODF)。該數(shù)據(jù)文件包含美國深空網(wǎng)(DSN)對LP的Doppler跟蹤觀測數(shù)據(jù)，包括距離變率(測速)和距離(測距)兩種數(shù)據(jù)。距離變率數(shù)據(jù)采樣間隔為10 s，精度較高(約為0.3 mm/s)；距離數(shù)據(jù)采樣間隔為33 s，精度平均約為0.5 m。由于最終要使用GEODYN Ⅱ軟件處理觀測數(shù)據(jù)，因此需要將ODF中的數(shù)據(jù)記錄轉(zhuǎn)換為GEODYN Ⅱ可識別的格式。首先根據(jù)數(shù)據(jù)格式說明文件TRK-2-18將ODF由二進制文件轉(zhuǎn)為ASCII文件。ODF中的距離變率數(shù)據(jù)以Hz為單位,距離數(shù)據(jù)則以range unit(RU)為單位，為了使其符合GEODYN Ⅱ?qū)τ跀?shù)據(jù)格式的特殊要求，需要按下式將二者單位分別轉(zhuǎn)換為m/s和m

Doprate=0.5×C/(T×Fref)×Dopfreq

(1)

Doprange=0.5×C/(FRU×Fref)×Doprange

(2)

式中，C表示光速；T為轉(zhuǎn)換器上下行信號比，S波段應(yīng)取為240/221；Fref為基準頻率，可由數(shù)據(jù)文件中讀取，其值接近于2093 MHz；FRU取值為0.5。由于ODF實際上為Ramp格式數(shù)據(jù)，其基準頻率通常是隨時間線性變化的。但在實際計算過程中，注意到絕大部分情況下基準頻率的變率實際上為零，因此可以將其作為普通的Doppler數(shù)據(jù)來處理，即基準頻率為常量。

本文研究所采用的衛(wèi)星定軌及大地參數(shù)解算軟件為GEODYN Ⅱ和SOLVE。GEODYN Ⅱ軟件主要由3部分構(gòu)成，即TDF、IIS和IIE。TDF是格式轉(zhuǎn)換器，能夠?qū)⒛撤N數(shù)據(jù)記錄轉(zhuǎn)換為GEODYN II可處理的格式。但由ODF轉(zhuǎn)換得到的ASCII數(shù)據(jù)記錄格式仍然不能被TDF識別，因此需要利用程序再次進行轉(zhuǎn)換。經(jīng)過TDF轉(zhuǎn)化后，觀測數(shù)據(jù)文件的格式變?yōu)镚EODYN II binary格式，該格式可以被GEODYN II所識別。根據(jù)情況進行弧段劃分(Arc1，Arc2，…，Arcn)。由于月球背面存在重力信息空白，軌道弧段不宜過長，一般取為2 d。還要注意同一弧段中不應(yīng)包括對衛(wèi)星的人工操作，即衛(wèi)星機動的時間段應(yīng)排除在連續(xù)弧段之外。上述數(shù)據(jù)預(yù)處理過程如圖1所示。SOLVE嚴格意義上應(yīng)該屬于GEODYN Ⅱ的一部分，其主要功能是進行矩陣求逆運算。圖1中的ReadODF與Trans2TDF為數(shù)據(jù)格式轉(zhuǎn)換程序。

圖1　LP 數(shù)據(jù)預(yù)處理流程

二、月球重力場模型解算

本文解算使用1998年1月11日至5月1日正常任務(wù)階段約4個月的Doppler跟蹤雙程、三程數(shù)據(jù)，所使用的LP數(shù)據(jù)稍多于LP75G(3個月)。解算流程如圖2所示，圖2中Ematrix表示法方程矩陣。所采用的力學模型、星歷及參數(shù)設(shè)置等見表1。用于解算GLGM2的數(shù)據(jù)主要來自軌道高度約為200km的Clementine任務(wù)。以該模型作為初始重力場模型對LP進行軌道積分時，顯然精度不夠，即使積分弧段為2d仍然會導致部分弧段不收斂。對于不收斂的弧段，先縮短其積分長度直至收斂[5-6]。利用所有積分弧段解算重力場模型，將更新后的重力場模型作為先驗重力場模型替換GLGM2[7-12]。更新后的重力場模型包含LP探測的重力場信息，可使所有2d積分弧段均收斂。

圖2　利用GEODYN II和SOLVE解算月球重力場位系數(shù)的具體流程及所需文件

ftn01行星星歷文件ftp:∥navigator.jpl.nasa.govJPLftn02集成了地磁、太陽通量等文件ftp:∥ftp.ngdc.noaa.govNOAAftn12重力場模型GLGM2(70階)PDSLemoineetal.1997ftn05控制文件太陽光壓模型:簡單柱狀模型多體攝動模型:太陽及其他主要模型精密初始軌道:利用spicetoolkit軟件提取并可同時進行坐標框架轉(zhuǎn)換 PDSNAIFDSN測站坐標:DSN16、24、27、54、61、66、34、42、46Folkner1997測站速度模型:ITRF93框架Folkner1997測站坐標改正:固體潮、海潮、極潮ftp:∥geodesy.gsfc.nasa.govLOVE:h2=0.609;k2=0.0852

續(xù)表1

首先，根據(jù)所選初始力模型(月球重力場模型、三體攝動模型、太陽光壓模型等)對LP軌道進行數(shù)值積分；然后便可得到線性觀測方程，該觀測方程建立了觀測數(shù)據(jù)(距離、距離變率數(shù)據(jù))與重力場位系數(shù)、初始狀態(tài)向量等未知參數(shù)之間的函數(shù)關(guān)系，并可形成法方程矩陣(Ematrix)；最終利用SOLVE對該法方程矩陣求逆，便可得到待求參數(shù)[13-14]。

三、模型分析

LPD75的模型階方差曲線與LP75G比較接近(如圖3所示)，但在低階次(7～20階)出現(xiàn)了明顯誤差，這在其誤差階方差曲線中也體現(xiàn)了出來。主要原因是由于GLGM2及LP75G都融合了多個歷史探月任務(wù)的數(shù)據(jù)，而本文解算的LPD75僅利用了LP的Doppler跟蹤數(shù)據(jù)，所恢復的月球重力場模型低階位系數(shù)精度較差。圖4顯示了LPD75和LP75G模型與月球地形模型STM359_grid-03的相關(guān)性。LPD75與地形的相關(guān)性在前50階次與LP75G相比總體上略差，50階次以上則優(yōu)于LP75G。顯然，軌道較低(正常任務(wù)階段月100 km)的LP任務(wù)對月球重力場短波信息較為敏感，所增加的一個月的LP數(shù)據(jù)，有助于提高模型高階位系數(shù)的精度，但對50階以下的重力場位系數(shù)沒有改進。

圖3　模型階方差及誤差階方差

由LPD75解算的月球自由空氣重力異常 (如圖5(a)所示)、月球大地水準面起伏(如圖6(a)所示)與LP75G(如圖5(b)、圖6(b)所示)可以看出，兩個模型的差別主要集中在月球背面和極區(qū)。LPD75的高頻噪聲明顯少于LP75G。

圖4　重力場模型與地形模型相關(guān)性

圖5　月表自由空氣重力異常及其差值

圖6　月球大地水準面起伏及其差值

四、結(jié)束語

本文基于LP Doppler數(shù)據(jù)，利用GEODYN Ⅱ和SOVLE軟件解算了月球重力場模型，得到了合理可靠的結(jié)果，驗證了數(shù)據(jù)處理、軟件使用的正確性。文中重點說明了其他學者忽略的LP數(shù)據(jù)預(yù)處理方法和過程，詳細闡述了月球重力場解算的流程，給出了必要文件及參數(shù)設(shè)置信息。除LP數(shù)據(jù)預(yù)處理過程外，文中所述GEODYN Ⅱ與SOLVE的解算過程均為通用過程，對于其他任務(wù)和類型的數(shù)據(jù)，只需要進對預(yù)處理過程和控制文件(ftn05)作相應(yīng)改動即可，可對今后利用GEODYN Ⅱ和SOLVE軟件進行月球重力場模型解算提供有借鑒。

參考文獻：

[1]BINDER A B. Lunar Prospector: Overview [J]. Science, 1998, 281(5382): 1475-1476.

[2]KONOPLIV A S, BANERDT W B, SJOGREN W L. Venus Gravity: 180th Degree and Order Model [J]. Icarus, 1999, 139(1): 3-18.

[3]KONOPLIV A S, ASMAR S W, CARRANZA E, et al. Recent Gravity Models as a Result of the Lunar Prospector Mission [J]. Icarus, 2001, 150(1): 1-18.

[4]GOOSSENS S, MATSUMOTO K. Lunar Satellite Orbit Determination Analysis and Quality Assessment from Lunar Prospector Tracking Data and Selene Simulations [J]. Advances in Space Research, 2007, 40(1): 43-50.

[5]豐海，李建成，李大煒，等.月球Airy均衡狀態(tài)與月殼厚度估計[J]. 測繪學報,2012 ,41 (4): 543-548.

[6]PING J S, HUANG Q, YAN J G, et al. Lunar Topographic Model Cltm-s01 from Chang’e-1 Laser ALtimeter [J]. Science in China Series G: Physics Mechanics and Astronomy, 2009, 52(7): 1105-1114.

[7]羅志才,李瓊,鐘波. 利用GRACE時變重力場反演黑河流域水儲量變化[J]. 測繪學報，2012,41 (5): 676-681.

[8]GUO J Y, SUN Y, CHANG X T,et al. Lunar Deflections of the Vertical and Their Distribution [J]. Natural Science, 2011, 3(5): 339-343.

[9]ARAKI H, TAZAWA S, NODA H, et al. Lunar Global Shape and Polar Topography Derived from Kaguya-LALT Laser Altimetry [J]. Science, 2009, 323(5916): 897-900.

[10]李森，安智明. 克里金模型在區(qū)域重力數(shù)據(jù)插值中的應(yīng)用[J]. 測繪通報, 2013 (10): 63-66.

[11]MAZARICO E, LEMOINE F, HAN S C et al. GLGM-3: A Degree-150 Lunar Gravity Model from the Historical Tracking Data of NASA Moon Orbiters [J]. Journal of Geophysical Research, 2010(115): E05001.

[12]FOLKNER W M, WILLIAMS J G, BOGGS D H. The Planetary and Lunar Ephemeris de 421 [J]. IPN Progress Report, 2008:172-178.

[13]曹金國，王建斌，戴山嶺. 重力場擾動位泰勒展開式低階項應(yīng)用研究[J].測繪通報，2014(1):54-55.

[14]GOOSSENS S, MATSUMOTO K. Lunar Satellite Orbit Determination Analysis and Quality Assessment from Lunar Prospector Tracking Data and Selene Simulations [J]. Advances in Space Research, 2007, 40(1): 43-50.

引文格式： 黃昆學，常曉濤，孫玉. 利用Lunar Prospector Doppler數(shù)據(jù)解算并分析月球重力場[J].測繪通報，2016(1)：19-22.DOI:10.13474/j.cnki.11-2246.2016.0005.

作者簡介：黃昆學(1985—)，男，博士生，主要研究方向為月球重力場解算。E-mail：shulang2003@163.com

基金項目：國家自然科學基金 (41204007); 國家高技術(shù)研究發(fā)展計劃(2013AA122502); 國家重點基礎(chǔ)研究發(fā)展計劃(2013CB733302); 中國地質(zhì)調(diào)查局地質(zhì)調(diào)查工作項目(12120113019100)

收稿日期：2014-11-25； 修回日期： 2015-10-16

中圖分類號：P228

文獻標識碼：B

文章編號：0494-0911(2016)01-0019-04