-
双程多普勒测速(以下简称双程测速)是行星探测器精密定轨中最为常用的观测模式之一，也是无线电科学实验(radio science experiment)中最主要的数据源[1-2]。
随着深空网结构的优化、频段的升级和服务能力的提高，加之先进的校正技术[3]和对来自于不同误差源噪声的详细建模[4]，双程测速的测量精度越来越高。国外的部分探测任务诸如“卡西尼-惠更斯”、“朱诺”和欧空局火星快车号(Mars express, MEX)等在其特定的任务时段，双程测速精度均可以达到10-2 ~10-3 mm/s[5-8]。
精密定轨的精度同时也与模型理论值的计算精度有关。传统的双程测速的理论值由多普勒频移周计数始点的双程平均距离到周计数终点的双程平均距离之平均变化率计算[9-11]。在行星际的探测任务中，由于探测器距离地面测站遥远(距地球1亿km以上), 且计算机浮点数模型的精度有限，距离量和必要的距离差分量在计算机中只能精确表达至毫米量级，从而导致测速值的计算精度无法满足高精度定轨的需要。
在国外的研究中，文献[12]使用7节点的高斯求积公式数值求解距离差分量。对于典型的30 s多普勒计数间隔，该方法最终可以将双程测速值的计算精度提高至2×10-3 mm/s。但是该文献中并没有给出其详细的解法与计算公式。文献[3]基于Widrow提出的量化误差模型[13]对数值噪声进行建模并给出其统计特性，但并没有定量地讨论减少计算误差的策略。在国内的研究中，曹建峰等提出了一种间接使用中心天体、航天器的速度与加速度进行多普勒建模的方法[14]，但相关文献仍旧较少。
与此同时，国际上很多优秀定轨软件对国内用户限制[15-17]，无法与其进行更加详细的交叉对比验证。一种通用的方法是将定轨软件中的浮点数均改为四精度存储，但是由于定轨软件架构复杂，更改浮点变量的精度操作繁琐，且会极大增加内存占用和时间复杂度[3]。
本文建立行星际双程测速模型(以下简称星际模型)，高精度地计算双程测速的理论值，减少传统双程测速模型(以下简称传统模型)在计算过程中丢失的有效位，并将该模型在武汉大学自主研发的深空探测器精密定轨与重力场解算软件系统(Wuhan University deep-space orbit determination and gravity recovery system, WUDOGS)中进行了实现[18-19]，并以MEX探测任务为背景，利用该软件进行仿真测试，从计算精度和模拟定轨结果两个方面验证该模型的优越性。这一工作可以应用到后续高精度的火星、小行星等探测任务中。
HTML
-
Te3和Ts3分别表示多普勒周计数的结束时刻和起始时刻，若多普勒测速的时标位于Te3，则相应的载波信号被探测器在Te2时刻转发，地面测站在Te1时刻发射；同样，Ts3时刻对应的载波信号被探测器在Ts2时刻转发，地面测站在Ts1时刻发射。忽略探测器应答机的转发时延，则双程测速的理论计算值ρ(Te3)可以表示为：
式中, ρg(Te3)为几何双程测速值; ρr(Te3)为广义相对论时延等效距离变率改正; ρa(Te3)为大气时延等效距离变率改正; ρt(Te3)为时间尺度不同造成的等效距离变率改正。且各项满足：
式中，ρ2g和ρ4g为下行几何距离；ρ1g和ρ3g为上行几何距离；c为光速；ρ4r、ρ2r、ρ3r和ρ1r为广义相对论时延的等效距离改正；ρ4a、ρ2a、ρ3a和ρ1a为大气介质时延的等效距离改正；ΔTe和ΔTs代表收发信号时刻(barycentric dynamical time, TDB)和世界标准时间(coordinated universal time, UTC)两种时间尺度之差; tc为多普勒积分间隔。
-
模型的计算精度(model accuracy)至少需要比观测量的随机噪声水平小一个量级[20]。以MEX任务为例，10 s采样率的双程测速理论值的模型计算精度至少需要达到0.005 mm/s的量级才能够满足高精度定轨的需要。而根据IEEE754浮点数标准，最常用的双精度浮点数在计算机内存中分配的二进制数位为64位，1位存储符号位，11位存储指数位，52位存储有效数位，相当于十进制的[lg(252+1)]=15位有效数字。有效数位不能表达的部分按照最近舍入法则进行舍入或者截断[21]。
根据式(1)和式(2)，若想计算ρ(Te3)，首先需要计算7个差分值：ρ4g－ρ2g、ρ3g－ρ1g、ρ4r－ρ2r、ρ3r－ρ1r、ρ4a－ρ2a、ρ3a－ρ1a和cΔTe－cΔTs，其中误差最大的是前两个差分值的计算。以火星任务为例，由于地面测站至探测器的距离至少达到1011 m的量级[19]，ρ4g、ρ2g、ρ3g和ρ1g的整数部分至少占用12个十进制有效位，小数部分至多占用3个十进制有效位，计算机至多能有效表达至毫米量级，将两数直接相减，则ρ4g－ρ2g和ρ3g－ρ1g两个差分值也至多能有效表达至毫米量级；若多普勒积分间隔为10 s，则ρg(Te3)至多能有效表达至0.1 mm/s量级。这样即使式(1)中的后3项可以在计算机中表达的十分精确，计算值ρ(Te3)也只能有效表达至0.1 mm/s量级，简单的采用式(1)和(2)，显然无法满足行星探测器高精度定轨的需求。
-
星际双程测速模型的实质在于规避两个量级在1011及以上的大数相减的情况，从而高精度地计算差分值ρ4g－ρ2g和ρ3g－ρ1g。
-
以ρ4g－ρ2g为例，给出如下高精度的计算方法。
按照式(3)计算ρ2g和ρ4g[14]：
式中，XM(Te2)和XM(Ts2)为Te2和Ts2时刻中心天体质心相对于太阳系质心的位置；XS(Te2)和XS(Ts2)为Te2和Ts2时刻探测器相对于中心天体质心的位置；XB(Te3)和XB(Ts3)为Te3和Ts3时刻地月质心相对于太阳系质心的位置；XE(Te3)和XE(Ts3)为Te3和Ts3时刻地球质心相对地月质心的位置；XA(Te3)和XA(Ts3)为Te3和Ts3时刻地面测站的地心天球坐标。
将向量用坐标分量形式展开，引入以下表示：ΔXS(ΔXS, ΔYS, ΔZS)、ΔXM(ΔXM, ΔYM, ΔZM)、ΔXB2(ΔXB2, ΔYB2, ΔZB2)、ΔXE2(ΔXE2、ΔYE2, ΔZE2)和ΔXA2(ΔXA2、ΔYA2, ΔZA2)。对于ΔXS和ΔXM的坐标分量，分别由对应的Te2和Ts2时刻的坐标分量相减得到；对于ΔXB2、ΔXE2、ΔXA2则分别由对应的Te3和Ts3时刻的坐标分量相减得到。进一步令：
可以得到：
由于距离较远，两时刻距离之间的差分值相对于距离本身为小量，故在Ts3时刻处，将式(5)沿Ts3时刻对应的下行支路方向进行泰勒展开。
令：
则Ts3时刻对应的下行支路径向方向可以用向量D2的方向向量α2(a2, b2, c2)表示。
令f1、f2和f3分别代表Ts3时刻f沿α2方向的一阶、二阶和三阶方向导数值。则式(5)可以写成如下形式：
考虑到实际求解精度，式(7)中忽略了比三阶高的无穷小O(|D2|3)项。
式(7)中3个方向导数的计算公式为：
式中，w2=a2Rx2+b2Ry2+c2Rz2，代入式(7)，即可算出ρ4g－ρ2g的值。
对于测距上行支路，将式(3)至式(8)相应的下行量替换为上行量即可计算出ρ3g－ρ1g的值。
-
在计算几何距离差分值的过程中涉及到不同时刻天体星历差分值的计算。对于火星、地球和地月质心本身而言，其距离太阳系质心遥远，若直接将两时刻对应的位置相减，同样会引起上述提到的数值误差。因此，本文采用切比雪夫多项式差分的方式计算ΔXM、ΔXB和ΔXE。
根据历表插值的方法[22]，两个时刻同一天体位置之差可以表示为：
式中, N代表切比雪夫多项式的阶数; τ为按照对应时刻所在的时间区间标准化至(-1, 1)范围内的时刻; Ck(τ)为τ对应的切比雪夫多项式系数; Pi, k(τ)为τ对应的各坐标分量的切比雪夫多项式。
令$\left\{ \begin{array}{l} \Delta {C_k} = {C_k}\left( {{\tau _2}} \right) - {C_k}\left( {{\tau _1}} \right)\\ \Delta {P_{i,k}} = {P_{i,k}}\left( {{\tau _2}} \right) - {P_{i,k}}\left( {{\tau _1}} \right) \end{array} \right.$, 则有:
式中，ΔPi, k按照如下递推公式计算：
式中, k=2, 3…N-1。
当t2和t1两个时刻位于同一区间内，ΔCkPi, k(τ2)=0。式(10)可以写为：
至此，从式(3)至式(12)，严格避免了两个数量级在1011及以上的大数字相减，ρ4g－ρ2g和ρ3g－ρ1g的值可以通过该高精度算法计算得出。
-
不考虑式(1)中后3项具体计算细节，给出星际模型计算ρg(Te3)的详细步骤，如图 1所示。
Figure 1. Detailed Calculation Steps of Interplanetary Model
3.1. 几何距离差分值的高精度计算方法
3.2. 天体星历差分值的计算
3.3. 星际双程测速模型计算ρg(Te3)的具体步骤
-
双程测速模型计算值的具体策略描述如下：
1) 利用SPICE[23]读取某时刻MEX探测器的位置和速度作为初始轨道根数；
2) 给定力模型、积分弧段和积分步长，积分出MEX探测器的轨道；
3) 将积分出的轨道分别代入传统和星际双程测速模型，生成两模型的计算值。
为了定量分析两种模型的计算精度，需要一个“准确”的理论值作为参考。其中，美国戈达德宇航中心的GEODYN-Ⅱ软件系统广泛地应用于行星探测器的精密定轨、测地学参数的估计、卫星轨道预报、仪器校准等深空探测领域，经受过数次行星际实测任务的考验[15, 18-19, 24]，可以认为GEODYN-Ⅱ的双程测速模型可靠准确，基本参数设置参照文献[19]，将GEODYN-Ⅱ生成的计算值作为“准确”的理论计算值。
将两模型的计算值分别与GEODYN-Ⅱ生成的计算值作差，得到计算值的残差。图 2分别给出了两模型的计算值残差图，表 1给出了残差的统计信息，RMS表示均方根误差。由于测试中均使用同一组轨道且策略严格保持一致，所以残差可以直接显示出两模型对于计算机浮点误差的敏感程度。对于传统模型，残差均值为-0.007 005 mm/s，超过计算值精度需要的0.005 mm/s，计算机的浮点误差最大达到0.111 22 mm/s，超过了MEX 10 s采样率的测量精度，残差的RMS为0.040 993 mm/s，这表明残差围绕均值在十分接近于测量精度的范围内分布，传统测速模型对浮点误差具有较强的敏感性。这也同样对应于浮点数运算规则的特性：在两个相近的大数相减时，有效位数将会大量缺失，而最近舍入法则会让这种数值噪声呈现一定的随机性[3]。对于星际模型，计算机浮点误差的量级最大为5×10-3 mm/s，均值为-0.003 43 mm/s，而RMS比均值还要小一个量级。在现有的观测精度下，该量级的数值噪声完全可以忽略不计，这是因为星际模型规避了对浮点误差最具敏感性的大数相减运算。距离差分值通过在多普勒计数起始点沿支路方向泰勒展开，从而转化为“小项”相加得到。星际模型能够将计算值的精度提高2倍，在一定程度上保证了计算值的“纯净度”。
Figure 2. Residual Plot of Calculation Values Using Two Models
模型 计算值残差/mm·s-1 最大值 均值 RMS 传统双程测速 0.111 22 -0.007 005 0.040 993 星际双程测速 5.000 0×10-3 -3.429 86×10-3 9.419 85×10-4 Table 1. Statistical Information of Calculation Residuals
为定量分析数值误差对定轨的影响，使用GEODYN-Ⅱ生成一组3 d弧段的模型计算值作为输入到WUDOGS中的真实观测值，在初始轨道根数的位置分量上分别加上-100 m、+100 m、+100 m的偏差，分别采用两种模型进行模拟定轨，得到最后一次迭代的观测值残差图，如图 3所示。表 2给出了两种模型模拟定轨的初轨改正数，表 3和表 4分别给出了两种模型的模拟定轨结果与给定初始轨道根数的差值以及最后一次迭代的观测值残差统计信息。
Figure 3. Post-fit Residuals of Orbit Determination Using Two Models
模型 ΔPx/m ΔPy/m ΔPz/m ΔVx/m·s-1 ΔVy/m·s-1 ΔVz/m·s-1 传统模型 100.097 80 -100.104 90 -99.798 70 2.67×10-5 -6.32×10-7 -1.29×10-5 星际模型 99.999 98 -99.999 98 -100.000 03 -4.89×10-9 2.58×10-11 2.29×10-9 Table 2. Initial Orbit Corrections of POD Using Two Models, Respectively
模型 ΔPx/m ΔPy/m ΔPz/m ΔVx/m·s-1 ΔVy/m·s-1 ΔVz/m·s-1 传统模型 0.097 809 -0.104 888 0.201 331 2.670×10-5 -6.317×10-7 -1.287×10-5 星际模型 -1.764×10-5 1.934×10-5 -3.640×10-5 -4.890×10-5 2.581×10-11 2.294×10-9 Table 3. Initial Orbit Differences Between POD Simulation Results and Given Results Using Two Models, Respectively
模型 最后一次迭代观测值残差均值 最后一次迭代观测值残差RMS 传统模型 -0.004 1 0.003 1 星际模型 4.589 5×10-8 4.768 8×10-7 Table 4. Statistical Results of POD Measurement Residuals During the Last Iteration Using Two Models, Respectively/mm·s-1
由于模拟定轨中没有引入其他形式的误差，理论上，若计算不存在浮点误差，最后一次迭代的观测值残差应接近于0，且初轨位置分量的改正数应接近于+100 m、-100 m和-100 m，但是数值噪声的存在会造成与理论结果的不一致性。据表 3和表 4显示，传统模型的定轨结果与给定初始轨道位置的差值在分米量级，最后一次迭代收敛后残差的均值和RMS分别为-0.004 1 mm/s和0.003 1 mm/s，星际模型的定轨结果与给定初始轨道位置的差值在10-5 m量级，最后一次迭代收敛后残差的均值和RMS分别为4.589 5×10-8 mm/s和4.768 8×10-7 mm/s，且从图 3中可以看出，残差量级最大不超过2×10-6 mm/s，可以认为星际模型中的浮点误差成分极其微弱，在定轨过程中不会引入过大的数值噪声，即星际模型能够更好地提取观测值中的有用信息，使残差成分在最小二乘的意义下更“纯净”。
-
针对距离相对遥远的行星际精密定轨任务，本文建立星际双程测速模型，高精度地计算了两个几何距离差分值，得到如下结论：
1) 星际双程测速模型能避免两个相近的大数相减，从而高精度地计算多普勒周计数终点和始点上行与下行几何距离之间的差分值，避免传统双程测速模型计算机浮点数误差过大的问题。
2) 各大天体的星历在短时间间隔内的差分值需由切比雪夫差分多项式及其相应系数的线性组合计算，同样可以避免计算机浮点数误差过大的问题，是星际双程测速模型的一个必要步骤。
3) 由于计算机浮点数模型本身无法精确表达大部分实数，星际双程测速模型只能尽可能地减弱计算机浮点误差的影响，并不能完全消除。根据本文的仿真结果，星际双程测速模型将模型计算值在计算机中表达的精度提高2个量级，避免定轨过程中引入过大的数值误差，可以为后续行星际深空探测任务的高精度定轨提供参考。
DownLoad:
