楊軒，鄢建國，葉茂，金煒桐，曲春凱，劉素艷

（武漢大學(xué) 測繪遙感信息工程國家重點(diǎn)實(shí)驗(yàn)室，武漢 430079）

0 引 言

深空多種模式的探測如軌道環(huán)繞、探測器著陸、巡游、采樣、浮空探測等，可以獲取更為全面和更深層的科學(xué)探測數(shù)據(jù)[1]。從現(xiàn)實(shí)和長遠(yuǎn)來看，對(duì)深空的探測和開發(fā)具有十分重要的科學(xué)和經(jīng)濟(jì)意義。深空探測將是21世紀(jì)人類進(jìn)行空間資源開發(fā)與利用、空間科學(xué)與技術(shù)創(chuàng)新的重要途徑。隨著探測任務(wù)復(fù)雜性的提升，一個(gè)任務(wù)中會(huì)同時(shí)涉及多個(gè)探測目標(biāo)。多目標(biāo)的探測技術(shù)往往需要多個(gè)探測器協(xié)同工作。對(duì)多目標(biāo)同時(shí)定軌定位是保證其探測任務(wù)成功實(shí)施首先需要解決的問題。比較典型的多目標(biāo)探測任務(wù)包括地球衛(wèi)星任務(wù)GRACE（Gravity Recovery and Climate Experiment）和月球探測任務(wù)SELENE（SELenological and ENgineering Explorer）以及GRAIL（Gravity Recovery and Interior Laboratory）。GRACE與GRAIL任務(wù)中，通過在兩顆探測器之間進(jìn)行高精度的相對(duì)距離測定，為反演高精度的地球與月球重力場提供了必要的數(shù)據(jù)[2-3]。SELENE任務(wù)中，第1次使用了同波束VLBI（Very Long Baseline Interferometry）技術(shù)同時(shí)對(duì)兩顆小探測器Rstar與Vstar觀測，定軌精度達(dá)到10 m左右的水平；除此之外，第1次對(duì)兩顆軌道器成功實(shí)施了四程多普勒觀測，在國際上首次獲得了月球背面的重力場信號(hào)，提高了月球背面重力場的分辨率[4-6]。對(duì)于火星探測來說，以往大都是對(duì)單個(gè)目標(biāo)進(jìn)行研究，包括對(duì)探測器進(jìn)行定軌，或?qū)χ懫骶?xì)定位，沒有對(duì)探測器和著陸器之間的鏈路進(jìn)行研究。在不考慮火星大氣和電離層介質(zhì)影響、火星自轉(zhuǎn)模型誤差、火星重力場模型誤差等情況下，劉慶會(huì)等討論了利用同波束VLBI同時(shí)對(duì)火星探測器和火星車定位的可能性，通過弧段內(nèi)數(shù)小時(shí)的軌道器測距，以及火星車和軌道器的同波束VLBI時(shí)延數(shù)據(jù)實(shí)現(xiàn)了定位誤差幾百米的火星車絕對(duì)定位與位置誤差幾十米的軌道器定軌[1]。

本文在目前研究的基礎(chǔ)上，根據(jù)一種月球與火星探測多程微波測量鏈路，從幾何構(gòu)型出發(fā)，仿真分析了此多程鏈路對(duì)典型月球與火星探測器軌道定軌，著陸器定位的精度貢獻(xiàn)。文章首先給出了該多程微波測量鏈路模式的詳細(xì)數(shù)學(xué)模型；通過2個(gè)仿真算例，討論了多程多普勒模式相比于傳統(tǒng)雙程測量模式對(duì)月球與火星探測多目標(biāo)定軌定位的貢獻(xiàn)。

1 多程微波測量鏈路模型

與傳統(tǒng)雙程跟蹤模式的不同之處在于，此多程微波測量鏈路模型測量模式增加了火星著陸器與軌道器之間一條鏈路，如圖 1所示。地面觀測跟蹤站在Ti時(shí)刻發(fā)送一個(gè)上行信號(hào)給火星軌道器，軌道器收到信號(hào)后，在Sj時(shí)刻，立刻把信號(hào)轉(zhuǎn)發(fā)給火星著陸器，著陸器在Lk時(shí)刻收到信號(hào)后，向地球發(fā)射轉(zhuǎn)發(fā)信號(hào)，最后地球測站在Tn時(shí)刻收到下行信號(hào)。

圖1 多程微波測量鏈路模式示意圖Fig.1 The diagram of the multi-way Doppler tracking mode

與雙程跟蹤模式類似，多程微波測量鏈路測量模型也采用“回溯法”來求取計(jì)算值，即按照Tn-Lk-Sj-Ti的路徑。地面跟蹤站在Tn時(shí)刻接收到下行信號(hào)并記錄時(shí)間，其余的時(shí)刻需要進(jìn)行光行時(shí)解算依次獲取。下行信號(hào)到達(dá)Tn時(shí)刻記錄的為協(xié)調(diào)世界時(shí)UTC（Universal Time Coordinated）時(shí)標(biāo)，首先需要將該時(shí)刻轉(zhuǎn)換到太陽系質(zhì)心力學(xué)時(shí)TDB（Barycentric Dynamic Time）時(shí)刻，得到TDB（n）。之后通過光行時(shí)解算依次得到TDB（m）、TDB（k）、TDB（j）、TDB（i），最后在轉(zhuǎn)換為Ti處UTC時(shí)標(biāo)。整個(gè)計(jì)算過程是在BCRS（Barycentric Celestial Reference System）下進(jìn)行的，并考慮了廣義相對(duì)論時(shí)延。則距離模型可表示為[7]

其中：c表示光速，X（Sj）表示j時(shí)刻軌道器的位置狀態(tài)量；X（Ti）表示i時(shí)刻地面跟蹤站的位置狀態(tài)量；X（Lk）表示k時(shí)刻著陸器的位置；X（Tn）描述n時(shí)刻地面跟蹤站的位置；R1表示信號(hào)從地面深空站傳送到火星軌道器的幾何距離；R2表示信號(hào)從火星軌道器傳送到火星著陸器的幾何距離；R3表示信號(hào)從火星著陸器傳送到地面深空站的幾何距離；RLTnk與RLTij為由相對(duì)論引起的時(shí)間延遲，即Shapiro時(shí)延[8-9]，相關(guān)的公式可參考文獻(xiàn)[7]第8章的內(nèi)容，在此文中不再推導(dǎo)。所有的狀態(tài)量都是在太陽系質(zhì)心坐標(biāo)系（BCRS）下描述的。由于深空站記錄的觀測值是在UTC時(shí)間下給出的，所以需要進(jìn)行相應(yīng)的時(shí)標(biāo)改正，式1第4和第5項(xiàng)即UTC到TDB的時(shí)標(biāo)改正。

對(duì)于速度值的計(jì)算，一般按照積分多普勒的形式表示，在一個(gè)多普勒積分周期內(nèi)，起始時(shí)刻Ts和終止時(shí)刻Te各自對(duì)應(yīng)的距離差相對(duì)于時(shí)間的變化率。在Ts和Te時(shí)刻，分別由式（1）可得到Rs與Re，則速度大小可表示為

在Sj和Lk時(shí)刻，距離計(jì)算值對(duì)軌道器和著陸器向量的偏導(dǎo)數(shù)為

速度計(jì)算值對(duì)軌道器和著陸器向量的偏導(dǎo)數(shù)可由式（5）推導(dǎo)

其中：下標(biāo)s和e分別表示在多普勒積分起始時(shí)刻和終止時(shí)刻探測器的狀態(tài)。

2 多程微波測量模式仿真

本文的仿真計(jì)算平臺(tái)使用了武漢大學(xué)深空探測器精密定軌與重力場解算軟件系統(tǒng)（Wuhan University Deep-space Orbit Determination and Gravity Recovery System，WUDOGS）[10-12]，文獻(xiàn)[11]中詳細(xì)介紹了該系統(tǒng)，并且與國際上主流的深空探測器精密定軌系統(tǒng)GEODYN-II進(jìn)行了嚴(yán)格的交叉驗(yàn)證測試，表明其精度已初步達(dá)到GEODYN-II的精度水平。WUDOGS采用了經(jīng)典動(dòng)力學(xué)法對(duì)探測器精密定軌，使用加權(quán)最小二乘法解算待估參數(shù)得到單弧段解，之后聯(lián)合多個(gè)單弧段法方程求解得到全局參數(shù)，整個(gè)過程需要進(jìn)行多次迭代直到結(jié)果收斂。

仿真過程中，真實(shí)的軌道由“標(biāo)準(zhǔn)”模型生成，此模型是事先選定并且具有目前最高的精度；通過真實(shí)軌道與相關(guān)多普勒模型，可仿真出多普勒觀測量；然后為“標(biāo)準(zhǔn)”模型添加噪聲和誤差用來模擬實(shí)際模型誤差，結(jié)合仿真的多普勒觀測量，進(jìn)行仿真定軌分析。在定軌計(jì)算中，軌道器6個(gè)初始軌道根數(shù)和每個(gè)弧段的測量偏差作為局部參數(shù)進(jìn)行估計(jì)，著陸器的3個(gè)位置參數(shù)則作為全局參數(shù)進(jìn)行估計(jì)。首先對(duì)每個(gè)弧段進(jìn)行精密定軌，解算局部待估參數(shù)，收斂后得到含有全局參數(shù)的法方程矩陣，之后融合多弧段方程進(jìn)行全局參數(shù)的求解，完成第1次全局迭代，之后經(jīng)過多次全局迭代直至全局參數(shù)的解算結(jié)果收斂，得到最終的著陸器位置，整個(gè)算法過程和重力場解算流程類似[13]。

考慮到觀測設(shè)備熱噪聲等，我們?cè)谀M觀測值上添加了標(biāo)準(zhǔn)差為1 mm/s的高斯白噪聲；并且考慮到現(xiàn)存模型誤差的影響，在仿真分析中，考慮了重力場誤差、坐標(biāo)轉(zhuǎn)換誤差、星歷誤差的影響，添加的誤差在每個(gè)算例中分別給出。

星上轉(zhuǎn)發(fā)延遲對(duì)多普勒計(jì)算值的影響在文獻(xiàn)[14]中進(jìn)行了詳細(xì)的研究。一般來說，星上轉(zhuǎn)發(fā)延遲的量級(jí)在微秒（μs）量級(jí)，如美國國家航空航天局（National Aeronautics and Space Administration，NASA）火星探測器“火星漫游者”（Mars Reconnaissance Orbiter，MRO）的轉(zhuǎn)發(fā)延遲為1.414 9 μs，歐洲空間局（European Space Agency，ESA）火星快車的轉(zhuǎn)發(fā)延遲為2.076 μs[14]。本文將此延遲設(shè)為10 μs，通過引入測量偏差來減弱這一延遲的影響。在本文仿真分析中沒有考慮包括火星大氣時(shí)延、地球電離層等的影響，這些在實(shí)際數(shù)據(jù)中可以作為改正項(xiàng)從數(shù)據(jù)預(yù)處理中剔除。

2.1 月球南極探測器與軌道器定位定軌

本節(jié)仿真討論了一個(gè)100 km高度的極軌圓軌道衛(wèi)星與位于月球南極的著陸器應(yīng)用此多程多普勒測量模式后能達(dá)到的精度水平。表 1為軌道器的開普勒軌道根數(shù)，著陸器位于月球南極附近（75°S，130°E）。仿真中考慮了包括月球660階非球形引力、各大天體的三體攝動(dòng)、固體潮攝動(dòng)、太陽光壓攝動(dòng)、地球扁率攝動(dòng)和月球扁率項(xiàng)間接攝動(dòng)以及相對(duì)論效應(yīng)的影響，具體描述由表 2給出。仿真解算中，模擬數(shù)據(jù)包括雙程多普勒與此多程多普勒，地面跟蹤站位于佳木斯與喀什，其中佳木斯深空站發(fā)射與接收多程微波多普勒信號(hào)，喀什深空站發(fā)射與接收雙程多普勒信號(hào)，采樣率均設(shè)置為10 s。

表1 月球探測器開普勒軌道根數(shù)Table 1 Kepler orbit elements of the Lunar spacecraft

表2 仿真分析中參數(shù)設(shè)置Table 2 A full description of the parameters used in the simulation of the Lunar spacecraft POD

考慮到目前研究中的模型存在誤差，為盡可能模擬實(shí)際情況，我們?cè)谠虑蛑亓瞿Ｐ偷南禂?shù)上添加了10倍的標(biāo)準(zhǔn)差。JPL發(fā)布的DE系列歷表是目前使用最為廣泛的行星歷表，由于激光測月數(shù)據(jù)的加入，使得其月球軌道精度處于亞米級(jí)，月球天平動(dòng)的精度可達(dá)到幾個(gè)毫角秒[15-17]，因此在月球位置三個(gè)分量上各添加1 m的偏差，為3個(gè)歐拉角添加10 mas的偏差。

此外，由于初軌也可能存在誤差，因此在軌道器初始位置向量X、Y、Z三個(gè)方向上各加100 m的偏差，著陸器分別在在月球主軸坐標(biāo)系X、Y、Z三個(gè)方向上添加100 m的偏差。定軌過程中初軌方向誤差增加為1 000 m仍然可以收斂，不影響本文的結(jié)論。由于是仿真實(shí)驗(yàn)，真實(shí)軌道與真實(shí)的著陸器坐標(biāo)已知，這里對(duì)定軌定位的精度評(píng)價(jià)通過重建軌道與真實(shí)軌道的差異、解算的坐標(biāo)與真實(shí)著陸器坐標(biāo)的差異給出。圖 2為雙程多普勒與多程多普勒觀測量占比統(tǒng)計(jì)圖。

圖2 月球探測器算例雙程多普勒與多程多普勒觀測量占比統(tǒng)計(jì)圖Fig.2 The observable number of the two-way Doppler and the multi-way Doppler in the Lunar spacecraft POD experiment

圖2中，可以發(fā)現(xiàn)多程多普勒觀測量比雙程多普勒觀測量少得多，大約只有雙程多普勒的5%左右，這是因?yàn)閷?shí)施此多程多普勒比雙程多普勒需要更嚴(yán)苛的條件，即地面跟蹤站–軌道器、軌道器–著陸器、著陸器–地面跟蹤站之間同時(shí)可見，由于具有100 km的軌道高度，每次軌道器過境，與著陸器通視的時(shí)間只有大約10 min，因此限制了此多程多普勒觀測量的生成。表 3～4為仿真定軌與著陸器定位的結(jié)果。

表3列出了只使用雙程多普勒數(shù)據(jù)模擬定軌與使用雙程多普勒結(jié)合此多程多普勒數(shù)據(jù)聯(lián)合定軌的結(jié)果。結(jié)合圖 2雙程多普勒與多程多普勒觀測量占比統(tǒng)計(jì)圖，發(fā)現(xiàn)由于觀測條件的限制，弧段3、弧段6和弧段7中并不含有此多程多普勒數(shù)據(jù)，得到了與只含有雙程多普勒數(shù)據(jù)一致的結(jié)果，在只使用雙程多普勒算例中，弧段5未能得到收斂解，猜測原因是因?yàn)榇嘶《蝺?nèi)數(shù)據(jù)的不連續(xù)，需要解算的測量偏差達(dá)到8個(gè)，導(dǎo)致了法方程求逆的不穩(wěn)定?；《?、2、4和弧段5得益于此多程多普勒數(shù)據(jù)的加入，軌道初軌與真實(shí)軌道初軌的偏差均得到了不同程度的減小，提升的精度最大超過一個(gè)量級(jí)。選擇了一個(gè)有代表性的弧段（弧段2），對(duì)這個(gè)弧段的重建軌道精度進(jìn)行分析。圖 3（a）為只使用雙程多普勒數(shù)據(jù)模擬定軌得到的重建軌道，圖 3（b）為雙程多普勒數(shù)據(jù)與此多程多普勒數(shù)據(jù)聯(lián)合定軌得到的重建軌道。很明顯，與只使用雙程多普勒數(shù)據(jù)相比，添加此多程多普勒數(shù)據(jù)后，精密定軌得到的軌道與真實(shí)軌道的差值減小了一個(gè)量級(jí)以上，軌道三維精度最大可達(dá)厘米級(jí)。軌道切向（T）與法向（N）方向的精度得到了提升，而在軌道徑向（R），由于精度已足夠高，未能得到顯著的提升。這一結(jié)果與文獻(xiàn)[18]仿真分析的四程多普勒定軌精度類似。

表3 月球探測器算例重建軌道與真實(shí)軌道初軌差值Table 3 Initial orbit position difference between reconstructed orbits and true one for the Lunar spacecraft POD experiment m

表4 著陸器的迭代計(jì)算過程Table 4 Iterative calculation of lander position m

進(jìn)一步的，考慮到在仿真數(shù)據(jù)中包含有歷表誤差、月球天平動(dòng)誤差以及星上轉(zhuǎn)發(fā)延遲的影響，我們考慮了若不使用測量偏差來吸收這些誤差，可以得到軌道精度的情況。表 5比較了這兩種不同定軌設(shè)置的結(jié)果。

從表 5可以看出，在含有歷表誤差、月球天平動(dòng)誤差以及星上轉(zhuǎn)發(fā)延遲的影響下，若不解算測量偏差，影響一般在厘米到米之間，得益于激光測月數(shù)據(jù)的加入，使得雙程多普勒數(shù)據(jù)處理過程中可忽略這些誤差的影響[19]；但是此多程多普勒模型中著陸器與軌道器之間存在一條鏈路，因此對(duì)月球天平動(dòng)誤差以及著陸器位置誤差特別敏感，在計(jì)算的過程中必須加入測量偏差以保證結(jié)果的正常收斂。

圖3 弧段2重建軌道與真實(shí)軌道的差異Fig.3 The differences of arc 2 between true orbit and reconstructed orbits

圖4 火星探測器算例雙程多普勒與多程多普勒觀測量占比統(tǒng)計(jì)圖Fig.4 The observable number of the two-way Doppler and the multi-way Doppler in the Mars spacecraft POD experiment

圖5 弧段3重建軌道與真實(shí)軌道的差異Fig.5 The differences of arc 3 between true orbit and reconstructed orbits

2.2 火星軌道器著陸器定軌定位

此算例中應(yīng)用了一個(gè)大偏心率小傾角橢圓軌道，表 6列出了開普勒軌道根數(shù)。仿真中著陸器的位置選在了亞馬遜平原（Amazonis Planitia）（30°N，165°W），模擬定軌參數(shù)見表 7，仿真解算的數(shù)據(jù)的設(shè)置與2.1節(jié)保持一致。

與月球算例相似，考慮到模型誤差的存在，我們?cè)诨鹦侵亓瞿Ｐ蚆RO120D上添加了10倍的中誤差；依據(jù)文獻(xiàn)[20]，為火星歷表添加1 km的偏差；依據(jù)文獻(xiàn)[21]，為火星固聯(lián)坐標(biāo)系至慣性系的轉(zhuǎn)換參數(shù)添加相應(yīng)的誤差。同樣的，在軌道器初始位置向量X、Y、Z三個(gè)方向上各加100 m的偏差，著陸器分別在火星固連坐標(biāo)系X、Y、Z三個(gè)方向上添加1 000 m的偏差。圖 4為此算例中雙程多普勒與多程多普勒觀測量占比統(tǒng)計(jì)圖。

表5 雙程多普勒解算測量偏差與不解算測量偏差結(jié)果比較Table 5 The comparison between solving the measurement bias and not solving the measurement bias for two-way Doppler m

由于仿真中使用了一個(gè)大偏心率小傾角橢圓軌道，此多程多普勒數(shù)據(jù)量相比100 km的月球極軌情況增加到大約占雙程多普勒數(shù)據(jù)的15%左右。表 8列出了只使用雙程多普勒數(shù)據(jù)模擬定軌與使用雙程多普勒結(jié)合此多程多普勒數(shù)據(jù)聯(lián)合定軌的結(jié)果。表 9給出了著陸器定位結(jié)果的收斂過程。

從表 8中，除了弧段1和弧段7的定軌精度有小幅降低外，其余5個(gè)弧段定軌的精度有了不同程度的提升，最大超過一個(gè)量級(jí)。在統(tǒng)計(jì)意義上仍然說明多程多普勒數(shù)據(jù)對(duì)火星探測器精密定軌的貢獻(xiàn)。從上面7個(gè)弧段中，我們也選擇了一個(gè)有代表性的弧段（弧段3）進(jìn)行詳細(xì)分析。圖 5（a）為只使用雙程多普勒數(shù)據(jù)模擬定軌得到的重建軌道，圖 5（b）為雙程多普勒數(shù)據(jù)與此多普勒數(shù)據(jù)聯(lián)合定軌得到的重建軌道。加入此多程多普勒數(shù)據(jù)后，軌道切向（T）與法向（N）方向的精度得到了大幅度的提升，在徑向（R）方向上，其精度并未得到顯著提升，這是因?yàn)檐壍缽较蛳啾确ㄏ蚝颓邢蚪８_，本身的誤差占比較小的緣故。從表 9中可知，對(duì)著陸器定位的模擬表明，與月球南極著陸器算例類似，收斂速度較快，只需3次全局迭代便可計(jì)算出米級(jí)精度的坐標(biāo)。

表6 火星探測器軌道根數(shù)Table 6 Kepler orbit elements of the Mars spacecraft

表7 模擬定軌參數(shù)設(shè)置Table 7 A full description of the parameters used in the simulation of the Mars spacecraft POD

表8 重建軌道與真實(shí)軌道初軌差值Table 8 Initial orbit position difference between reconstructed orbits and true one m

表9 著陸器的迭代計(jì)算過程Table 9 Iterative calculation of lander position m

由于添加了著陸器與軌道器之間的鏈路，本文研究的多程多普勒模式在軌道面與著陸器–軌道器方向上新增了一個(gè)約束，這一約束與軌道切向（T）與法向（N）均形成了一個(gè)小角度，為軌道切向與法向提供了額外的約束，彌補(bǔ)了雙程多普勒的不足。

3 結(jié) 論

本文對(duì)一種月球與火星探測多程微波測量鏈路進(jìn)行了初步的數(shù)值模擬分析，結(jié)果表明，結(jié)合傳統(tǒng)雙程多普勒數(shù)據(jù)與此多程多普勒數(shù)據(jù)進(jìn)行精密定軌，重建軌道的精度可在月球探測任務(wù)中穩(wěn)定在幾米左右，在火星探測任務(wù)中可達(dá)到數(shù)10 m精度的水平，著陸器位置精度可分別穩(wěn)定在分米和米級(jí)。

由于文章在模擬過程中考慮的是較為理想的狀態(tài)，仿真數(shù)據(jù)除地球與火星、月球遮擋外，不存在間斷的情況，保證了充足的數(shù)據(jù)量；除此之外，添加的模型誤差，如重力場誤差、歷表誤差等，都屬于形式誤差，有可能會(huì)存在低估這些誤差的情況；太陽光壓等其他攝動(dòng)力的誤差也未考慮，因此在實(shí)際應(yīng)用中，多程多普勒模型的精度可能會(huì)受到一定的限制。

本文仿真分析工作在火星著陸器定位模式上做了初步嘗試，有望對(duì)我國火星探測任務(wù)提供一定程度的參考。

[1]劉慶會(huì)，吳亞軍，黃勇，等.基于同波束VLBI 的火星車測定位技術(shù)[J].中國科學(xué)：物理學(xué) 力學(xué) 天文學(xué)，2015，45（9）：1-8.LIU Q H，WU Y J，HUANG Y，et al.Mars rover positioning technology based on same-beam VLBI[J].SCIENTIA SINICA Physica，Mechanica & Astronomica，2015，45（9）：1-8.

[2]ZHENG W，XU H Z，ZHONG M，et al.Precise recovery of the Earth’s gravitational field with GRACE：intersatellite range-rate interpolation approach[J].IEEE Geoscience and Remote Sensing Letters，2012，9（3）：422-426.

[3]ZUBER M T，SMITH D E，WATKINS M M，et al.Gravity field of the Moon from the gravity recovery and interior laboratory（GRAIL）mission[J].Science，2013，339（6120）：668-671.

[4]GOOSSENS S，MATSUMOTO K，ISHIHARA Y，et al.Results for orbit determination of the three satellites of Kaguya[J].Journal of the Geodetic Society of Japan，2009，55（2）：255-268.

[5]LIU Q，KIKUCHI F，MATSUMOTO K，et al.Same-beam VLBI observations of SELENE for improving lunar gravity field model[J].Radio Science，2010，45（1）：1-16.

[6]NAMIKI N，IWATA T，MATSUMOTO K，et al.Farside gravity field of the Moon from four-way Doppler measurements of SELENE（Kaguya）[J].Science，2009，323（5916）：900-905.

[7]MOYER T D.Formulation for observed and computed values of deep space network data types for navigation[M].USA：John Wiley & Sons，2005.

[8]SHAPIRO I I.Fourth test of general relativity[J].Physical Review Letters，1964，13（26）：789.

[9]TOMMEI G，MILANI A，VOKROUHLICKY D.Light-time computations for the BepiColombo radio science experiment[J].Celestial Mechanics and Dynamical Astronomy，2010，107（1-2）：285-298.

[10]葉茂，李斐，鄢建國，等.國內(nèi)外深空探測器精密定軌軟件研究綜述及WUDOGS簡介[J].飛行器測控學(xué)報(bào)，2017，36（1）：45-55.YE M，LI F，YAN J G，et al.Review of deep space orbit determination software and introduction to WUDOGS[J].Journal of Spacecraft TT&C Technology，2017，36（1）：45-55.

[11]YAN J G，YANG X，YE M，et al.Independent Mars spacecraft precise orbit determination software development and its application[J].Astrophysics and Space Science，2017，362，（47）：181-190.

[12]楊軒，鄢建國，葉茂，等.火星探測器精密定軌軟件研制及MEX數(shù)據(jù)處理[J].武漢大學(xué)學(xué)報(bào)信息科學(xué)版（已接收）.YANG X，YAN J G，YE M，et al.，Development of precise orbit determination software and data processing for MEX[J].Geomatics and Information Science of Wuhan University，accepted.

[13]鄢建國，李斐，平勁松.基于MGS測圖段部分弧段的精密定軌及火星重力場模型解算[J].測繪學(xué)報(bào)，2010，39（5）：484-490.YAN J G，LI F，PING J S.Precision orbit determination of mgs mapping phase arcs and martian gravity field model solution[J].Acta Geodaetica et Cartographica Sinica，2010，39（5）：484-490.

[14]BERTONE S，PONCIN-LAFITTE C L，ROSENBLATT P，et al.Impact analysis of the transponder time delay on radio-tracking observables[J].Advances in Space Research，2018，61（1）：89-96.

[15]RAMBAUX N，WILLIAMS J G.The Moon’s physical librations and determination of their free modes[J].Celestial Mechanics and Dynamical Astronomy，2011，109（1）：85-100.

[16]WILLIAMS J G，SLADE M A，ECKHARDT D H，et al.Lunar physical librations and laser ranging[J].Earth，Moon，and Planets，1973，8（4）：469-483.

[17]WILLIAMS，J G，BOGGS D H，F(xiàn)OLKNER W M.DE430 lunar orbit，physical librations and surface coordinates[R].USA：JPL Interoffice Memorandum， Jet Propulsion Laboratory，California Institute of Technology，Pasadena，California，2013，19.

[18]LI F，YE M，YAN J G，et al.A simulation of the four-way lunar lander–orbiter tracking mode for the Chang’E-5 mission[J].Advances in Space Research，2016，57（11）：2376-2384.

[19]曹建峰，張宇，胡松杰，等.嫦娥三號(hào)著陸器精確定位與精度分析[J].武漢大學(xué)學(xué)報(bào)（信息科學(xué)版），2016，41（2）：274-278.CAO J F，ZHANG Y，HU S J，et al.An analysis of precise positioning and accuracy of the CE-3 lunar lander soft landing[J].Geomatics and Information Science of Wuhan University，2016，41（2）：274-278.

[20]FOLKNER W M，WILLIAMS J G，BOGGS D H，et al.The planetary and lunar ephemerides DE430 and DE431[R].[S.l]：Interplanetary Network Progress Report，2014，196：1-81.

[21]KONOPLIV A S，PARK R S，F(xiàn)OLKNER W M.An improved JPL Mars gravity field and orientation from Mars orbiter and lander tracking data[J].Icarus，2016（274）：253-260.