对一种月球与火星探测多程微波测量链路定轨定位的数值模拟初步分析

杨轩，鄢建国，叶茂，金炜桐，曲春凯，刘素艳

（武汉大学 测绘遥感信息工程国家重点实验室，武汉 430079）

0 引 言

深空多种模式的探测如轨道环绕、探测器着陆、巡游、采样、浮空探测等，可以获取更为全面和更深层的科学探测数据[1]。从现实和长远来看，对深空的探测和开发具有十分重要的科学和经济意义。深空探测将是21世纪人类进行空间资源开发与利用、空间科学与技术创新的重要途径。随着探测任务复杂性的提升，一个任务中会同时涉及多个探测目标。多目标的探测技术往往需要多个探测器协同工作。对多目标同时定轨定位是保证其探测任务成功实施首先需要解决的问题。比较典型的多目标探测任务包括地球卫星任务GRACE（Gravity Recovery and Climate Experiment）和月球探测任务SELENE（SELenological and ENgineering Explorer）以及GRAIL（Gravity Recovery and Interior Laboratory）。GRACE与GRAIL任务中，通过在两颗探测器之间进行高精度的相对距离测定，为反演高精度的地球与月球重力场提供了必要的数据[2-3]。SELENE任务中，第1次使用了同波束VLBI（Very Long Baseline Interferometry）技术同时对两颗小探测器Rstar与Vstar观测，定轨精度达到10 m左右的水平；除此之外，第1次对两颗轨道器成功实施了四程多普勒观测，在国际上首次获得了月球背面的重力场信号，提高了月球背面重力场的分辨率[4-6]。对于火星探测来说，以往大都是对单个目标进行研究，包括对探测器进行定轨，或对着陆器精细定位，没有对探测器和着陆器之间的链路进行研究。在不考虑火星大气和电离层介质影响、火星自转模型误差、火星重力场模型误差等情况下，刘庆会等讨论了利用同波束VLBI同时对火星探测器和火星车定位的可能性，通过弧段内数小时的轨道器测距，以及火星车和轨道器的同波束VLBI时延数据实现了定位误差几百米的火星车绝对定位与位置误差几十米的轨道器定轨[1]。

本文在目前研究的基础上，根据一种月球与火星探测多程微波测量链路，从几何构型出发，仿真分析了此多程链路对典型月球与火星探测器轨道定轨，着陆器定位的精度贡献。文章首先给出了该多程微波测量链路模式的详细数学模型；通过2个仿真算例，讨论了多程多普勒模式相比于传统双程测量模式对月球与火星探测多目标定轨定位的贡献。

1 多程微波测量链路模型

与传统双程跟踪模式的不同之处在于，此多程微波测量链路模型测量模式增加了火星着陆器与轨道器之间一条链路，如图 1所示。地面观测跟踪站在Ti时刻发送一个上行信号给火星轨道器，轨道器收到信号后，在Sj时刻，立刻把信号转发给火星着陆器，着陆器在Lk时刻收到信号后，向地球发射转发信号，最后地球测站在Tn时刻收到下行信号。

图1 多程微波测量链路模式示意图Fig.1 The diagram of the multi-way Doppler tracking mode

与双程跟踪模式类似，多程微波测量链路测量模型也采用“回溯法”来求取计算值，即按照Tn-Lk-Sj-Ti的路径。地面跟踪站在Tn时刻接收到下行信号并记录时间，其余的时刻需要进行光行时解算依次获取。下行信号到达Tn时刻记录的为协调世界时UTC（Universal Time Coordinated）时标，首先需要将该时刻转换到太阳系质心力学时TDB（Barycentric Dynamic Time）时刻，得到TDB（n）。之后通过光行时解算依次得到TDB（m）、TDB（k）、TDB（j）、TDB（i），最后在转换为Ti处UTC时标。整个计算过程是在BCRS（Barycentric Celestial Reference System）下进行的，并考虑了广义相对论时延。则距离模型可表示为[7]

其中：c表示光速，X（Sj）表示j时刻轨道器的位置状态量；X（Ti）表示i时刻地面跟踪站的位置状态量；X（Lk）表示k时刻着陆器的位置；X（Tn）描述n时刻地面跟踪站的位置；R1表示信号从地面深空站传送到火星轨道器的几何距离；R2表示信号从火星轨道器传送到火星着陆器的几何距离；R3表示信号从火星着陆器传送到地面深空站的几何距离；RLTnk与RLTij为由相对论引起的时间延迟，即Shapiro时延[8-9]，相关的公式可参考文献[7]第8章的内容，在此文中不再推导。所有的状态量都是在太阳系质心坐标系（BCRS）下描述的。由于深空站记录的观测值是在UTC时间下给出的，所以需要进行相应的时标改正，式1第4和第5项即UTC到TDB的时标改正。

对于速度值的计算，一般按照积分多普勒的形式表示，在一个多普勒积分周期内，起始时刻Ts和终止时刻Te各自对应的距离差相对于时间的变化率。在Ts和Te时刻，分别由式（1）可得到Rs与Re，则速度大小可表示为

在Sj和Lk时刻，距离计算值对轨道器和着陆器向量的偏导数为

速度计算值对轨道器和着陆器向量的偏导数可由式（5）推导

其中：下标s和e分别表示在多普勒积分起始时刻和终止时刻探测器的状态。

2 多程微波测量模式仿真

本文的仿真计算平台使用了武汉大学深空探测器精密定轨与重力场解算软件系统（Wuhan University Deep-space Orbit Determination and Gravity Recovery System，WUDOGS）[10-12]，文献[11]中详细介绍了该系统，并且与国际上主流的深空探测器精密定轨系统GEODYN-II进行了严格的交叉验证测试，表明其精度已初步达到GEODYN-II的精度水平。WUDOGS采用了经典动力学法对探测器精密定轨，使用加权最小二乘法解算待估参数得到单弧段解，之后联合多个单弧段法方程求解得到全局参数，整个过程需要进行多次迭代直到结果收敛。

仿真过程中，真实的轨道由“标准”模型生成，此模型是事先选定并且具有目前最高的精度；通过真实轨道与相关多普勒模型，可仿真出多普勒观测量；然后为“标准”模型添加噪声和误差用来模拟实际模型误差，结合仿真的多普勒观测量，进行仿真定轨分析。在定轨计算中，轨道器6个初始轨道根数和每个弧段的测量偏差作为局部参数进行估计，着陆器的3个位置参数则作为全局参数进行估计。首先对每个弧段进行精密定轨，解算局部待估参数，收敛后得到含有全局参数的法方程矩阵，之后融合多弧段方程进行全局参数的求解，完成第1次全局迭代，之后经过多次全局迭代直至全局参数的解算结果收敛，得到最终的着陆器位置，整个算法过程和重力场解算流程类似[13]。

考虑到观测设备热噪声等，我们在模拟观测值上添加了标准差为1 mm/s的高斯白噪声；并且考虑到现存模型误差的影响，在仿真分析中，考虑了重力场误差、坐标转换误差、星历误差的影响，添加的误差在每个算例中分别给出。

星上转发延迟对多普勒计算值的影响在文献[14]中进行了详细的研究。一般来说，星上转发延迟的量级在微秒（μs）量级，如美国国家航空航天局（National Aeronautics and Space Administration，NASA）火星探测器“火星漫游者”（Mars Reconnaissance Orbiter，MRO）的转发延迟为1.414 9 μs，欧洲空间局（European Space Agency，ESA）火星快车的转发延迟为2.076 μs[14]。本文将此延迟设为10 μs，通过引入测量偏差来减弱这一延迟的影响。在本文仿真分析中没有考虑包括火星大气时延、地球电离层等的影响，这些在实际数据中可以作为改正项从数据预处理中剔除。

2.1 月球南极探测器与轨道器定位定轨

本节仿真讨论了一个100 km高度的极轨圆轨道卫星与位于月球南极的着陆器应用此多程多普勒测量模式后能达到的精度水平。表 1为轨道器的开普勒轨道根数，着陆器位于月球南极附近（75°S，130°E）。仿真中考虑了包括月球660阶非球形引力、各大天体的三体摄动、固体潮摄动、太阳光压摄动、地球扁率摄动和月球扁率项间接摄动以及相对论效应的影响，具体描述由表 2给出。仿真解算中，模拟数据包括双程多普勒与此多程多普勒，地面跟踪站位于佳木斯与喀什，其中佳木斯深空站发射与接收多程微波多普勒信号，喀什深空站发射与接收双程多普勒信号，采样率均设置为10 s。

表1 月球探测器开普勒轨道根数Table 1 Kepler orbit elements of the Lunar spacecraft

表2 仿真分析中参数设置Table 2 A full description of the parameters used in the simulation of the Lunar spacecraft POD

考虑到目前研究中的模型存在误差，为尽可能模拟实际情况，我们在月球重力场模型的系数上添加了10倍的标准差。JPL发布的DE系列历表是目前使用最为广泛的行星历表，由于激光测月数据的加入，使得其月球轨道精度处于亚米级，月球天平动的精度可达到几个毫角秒[15-17]，因此在月球位置三个分量上各添加1 m的偏差，为3个欧拉角添加10 mas的偏差。

此外，由于初轨也可能存在误差，因此在轨道器初始位置向量X、Y、Z三个方向上各加100 m的偏差，着陆器分别在在月球主轴坐标系X、Y、Z三个方向上添加100 m的偏差。定轨过程中初轨方向误差增加为1 000 m仍然可以收敛，不影响本文的结论。由于是仿真实验，真实轨道与真实的着陆器坐标已知，这里对定轨定位的精度评价通过重建轨道与真实轨道的差异、解算的坐标与真实着陆器坐标的差异给出。图 2为双程多普勒与多程多普勒观测量占比统计图。

图2 月球探测器算例双程多普勒与多程多普勒观测量占比统计图Fig.2 The observable number of the two-way Doppler and the multi-way Doppler in the Lunar spacecraft POD experiment

图2中，可以发现多程多普勒观测量比双程多普勒观测量少得多，大约只有双程多普勒的5%左右，这是因为实施此多程多普勒比双程多普勒需要更严苛的条件，即地面跟踪站–轨道器、轨道器–着陆器、着陆器–地面跟踪站之间同时可见，由于具有100 km的轨道高度，每次轨道器过境，与着陆器通视的时间只有大约10 min，因此限制了此多程多普勒观测量的生成。表 3～4为仿真定轨与着陆器定位的结果。

表3列出了只使用双程多普勒数据模拟定轨与使用双程多普勒结合此多程多普勒数据联合定轨的结果。结合图 2双程多普勒与多程多普勒观测量占比统计图，发现由于观测条件的限制，弧段3、弧段6和弧段7中并不含有此多程多普勒数据，得到了与只含有双程多普勒数据一致的结果，在只使用双程多普勒算例中，弧段5未能得到收敛解，猜测原因是因为此弧段内数据的不连续，需要解算的测量偏差达到8个，导致了法方程求逆的不稳定。弧段1、2、4和弧段5得益于此多程多普勒数据的加入，轨道初轨与真实轨道初轨的偏差均得到了不同程度的减小，提升的精度最大超过一个量级。选择了一个有代表性的弧段（弧段2），对这个弧段的重建轨道精度进行分析。图 3（a）为只使用双程多普勒数据模拟定轨得到的重建轨道，图 3（b）为双程多普勒数据与此多程多普勒数据联合定轨得到的重建轨道。很明显，与只使用双程多普勒数据相比，添加此多程多普勒数据后，精密定轨得到的轨道与真实轨道的差值减小了一个量级以上，轨道三维精度最大可达厘米级。轨道切向（T）与法向（N）方向的精度得到了提升，而在轨道径向（R），由于精度已足够高，未能得到显著的提升。这一结果与文献[18]仿真分析的四程多普勒定轨精度类似。

表3 月球探测器算例重建轨道与真实轨道初轨差值Table 3 Initial orbit position difference between reconstructed orbits and true one for the Lunar spacecraft POD experiment m

表4 着陆器的迭代计算过程Table 4 Iterative calculation of lander position m

进一步的，考虑到在仿真数据中包含有历表误差、月球天平动误差以及星上转发延迟的影响，我们考虑了若不使用测量偏差来吸收这些误差，可以得到轨道精度的情况。表 5比较了这两种不同定轨设置的结果。

从表 5可以看出，在含有历表误差、月球天平动误差以及星上转发延迟的影响下，若不解算测量偏差，影响一般在厘米到米之间，得益于激光测月数据的加入，使得双程多普勒数据处理过程中可忽略这些误差的影响[19]；但是此多程多普勒模型中着陆器与轨道器之间存在一条链路，因此对月球天平动误差以及着陆器位置误差特别敏感，在计算的过程中必须加入测量偏差以保证结果的正常收敛。

图3 弧段2重建轨道与真实轨道的差异Fig.3 The differences of arc 2 between true orbit and reconstructed orbits

图4 火星探测器算例双程多普勒与多程多普勒观测量占比统计图Fig.4 The observable number of the two-way Doppler and the multi-way Doppler in the Mars spacecraft POD experiment

图5 弧段3重建轨道与真实轨道的差异Fig.5 The differences of arc 3 between true orbit and reconstructed orbits

2.2 火星轨道器着陆器定轨定位

此算例中应用了一个大偏心率小倾角椭圆轨道，表 6列出了开普勒轨道根数。仿真中着陆器的位置选在了亚马逊平原（Amazonis Planitia）（30°N，165°W），模拟定轨参数见表 7，仿真解算的数据的设置与2.1节保持一致。

与月球算例相似，考虑到模型误差的存在，我们在火星重力场模型MRO120D上添加了10倍的中误差；依据文献[20]，为火星历表添加1 km的偏差；依据文献[21]，为火星固联坐标系至惯性系的转换参数添加相应的误差。同样的，在轨道器初始位置向量X、Y、Z三个方向上各加100 m的偏差，着陆器分别在火星固连坐标系X、Y、Z三个方向上添加1 000 m的偏差。图 4为此算例中双程多普勒与多程多普勒观测量占比统计图。

表5 双程多普勒解算测量偏差与不解算测量偏差结果比较Table 5 The comparison between solving the measurement bias and not solving the measurement bias for two-way Doppler m

由于仿真中使用了一个大偏心率小倾角椭圆轨道，此多程多普勒数据量相比100 km的月球极轨情况增加到大约占双程多普勒数据的15%左右。表 8列出了只使用双程多普勒数据模拟定轨与使用双程多普勒结合此多程多普勒数据联合定轨的结果。表 9给出了着陆器定位结果的收敛过程。

从表 8中，除了弧段1和弧段7的定轨精度有小幅降低外，其余5个弧段定轨的精度有了不同程度的提升，最大超过一个量级。在统计意义上仍然说明多程多普勒数据对火星探测器精密定轨的贡献。从上面7个弧段中，我们也选择了一个有代表性的弧段（弧段3）进行详细分析。图 5（a）为只使用双程多普勒数据模拟定轨得到的重建轨道，图 5（b）为双程多普勒数据与此多普勒数据联合定轨得到的重建轨道。加入此多程多普勒数据后，轨道切向（T）与法向（N）方向的精度得到了大幅度的提升，在径向（R）方向上，其精度并未得到显著提升，这是因为轨道径向相比法向和切向建模更精确，本身的误差占比较小的缘故。从表 9中可知，对着陆器定位的模拟表明，与月球南极着陆器算例类似，收敛速度较快，只需3次全局迭代便可计算出米级精度的坐标。

表6 火星探测器轨道根数Table 6 Kepler orbit elements of the Mars spacecraft

表7 模拟定轨参数设置Table 7 A full description of the parameters used in the simulation of the Mars spacecraft POD

表8 重建轨道与真实轨道初轨差值Table 8 Initial orbit position difference between reconstructed orbits and true one m

表9 着陆器的迭代计算过程Table 9 Iterative calculation of lander position m

由于添加了着陆器与轨道器之间的链路，本文研究的多程多普勒模式在轨道面与着陆器–轨道器方向上新增了一个约束，这一约束与轨道切向（T）与法向（N）均形成了一个小角度，为轨道切向与法向提供了额外的约束，弥补了双程多普勒的不足。

3 结 论

本文对一种月球与火星探测多程微波测量链路进行了初步的数值模拟分析，结果表明，结合传统双程多普勒数据与此多程多普勒数据进行精密定轨，重建轨道的精度可在月球探测任务中稳定在几米左右，在火星探测任务中可达到数10 m精度的水平，着陆器位置精度可分别稳定在分米和米级。

由于文章在模拟过程中考虑的是较为理想的状态，仿真数据除地球与火星、月球遮挡外，不存在间断的情况，保证了充足的数据量；除此之外，添加的模型误差，如重力场误差、历表误差等，都属于形式误差，有可能会存在低估这些误差的情况；太阳光压等其他摄动力的误差也未考虑，因此在实际应用中，多程多普勒模型的精度可能会受到一定的限制。

本文仿真分析工作在火星着陆器定位模式上做了初步尝试，有望对我国火星探测任务提供一定程度的参考。

[1]刘庆会，吴亚军，黄勇，等.基于同波束VLBI 的火星车测定位技术[J].中国科学：物理学 力学 天文学，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，VOKROUHLICKÝ D.Light-time computations for the BepiColombo radio science experiment[J].Celestial Mechanics and Dynamical Astronomy，2010，107（1-2）：285-298.

[10]叶茂，李斐，鄢建国，等.国内外深空探测器精密定轨软件研究综述及WUDOGS简介[J].飞行器测控学报，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数据处理[J].武汉大学学报信息科学版（已接收）.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].测绘学报，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，FOLKNER 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]曹建峰，张宇，胡松杰，等.嫦娥三号着陆器精确定位与精度分析[J].武汉大学学报（信息科学版），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，FOLKNER W M.An improved JPL Mars gravity field and orientation from Mars orbiter and lander tracking data[J].Icarus，2016（274）：253-260.