双精度浮点定义下的深空Doppler测速计算方法

摘要

本发明属于数据处理技术领域，公开了一种双精度浮点定义下的深空Doppler测速计算方法。该方法首先依据Doppler测量原理建立Doppler测量方程；然后用泰勒级数展开式来计算积分周期始末时刻的光行时差；最后将光行时差代入测量方程计算Doppler测速。本发明提高了深空探测Doppler测速的计算精度和效率。

说明书

技术领域

本发明属于数据处理技术领域，涉及一种双精度浮点定义下的深空Doppler(多普勒) 测速计算方法。

背景技术

目前深空探测中Doppler的计算通常采用《Formulation for observed and computed  values of Deep Space Network data types for navigation》一书中所述的方法。利用 该方法，若采用双精度浮点数计算光行时，由于舍入误差的影响，其计算精度约为皮秒左右(其中^[]表示取整，R为探测距离，单位为亿千米)，此精度在1亿千米以上的探 测任务Doppler测量数据的处理中，只能保证1mm/s(1秒积分周期)的计算精度，显然 不能满足当前0.1mm/s甚至0.02mm/s的精度要求。通常的做法是在计算光行时使用四精 度浮点数代替双精度浮点数，以提高计算精度，但该方法虽然提高了计算精度，却降低了 计算效率。

发明内容

本发明的目的在于提供一种双精度浮点定义下的深空Doppler测速计算方法，以提高 计算精度和计算效率。

为实现上述目的，本发明提供的双精度浮点定义下的深空Doppler测速计算方法包括 以下步骤：

<1>依据Doppler测量原理建立Doppler测量方程；

<2>用泰勒级数展开式来计算积分周期始末时刻的光行时差；

<3>将光行时差代入测量方程计算Doppler测速。

进一步的，所述Doppler测量方程如下：

$f_D=f_R-f_0=\frac{N_C}{\mathrm{\Delta T}}+\mathrm{\delta f}+v$

其中f_D表示Doppler频率偏移量，f_R表示接收端接收信号频率，f₀表示地面接收站 本振频率，N_C为一个积分周期内的周计数，ΔT为一个积分周期，δf为观测系统误差，v 为观测随机误差。

所述计算积分周期始末时刻光行时差的方法如下：

$\mathrm{\Delta \tau }={\tau }_{S,e}-{\tau }_{S,s}$

$=\frac{({\rho }_{\mathrm{SR},e}-{\rho }_{\mathrm{SR},s})}{c}+{\mathrm{RLT}}_e-{\mathrm{RLT}}_s$

$-\{{[{\mathrm{TDB}}_R-t\left({\tau }_R\right)]}_e-{[{\mathrm{TDB}}_R-t\left({\tau }_R\right)]}_s\}$

$+\{{[{\mathrm{TDB}}_S-t\left({\tau }_S\right)]}_e-{[{\mathrm{TDB}}_S-t\left({\tau }_S\right)]}_s\}$

其中，τ_S，e，τ_S，s分别表示星上信号发射端积分周期结束和起始时刻对应的原始， ρ_SR，e，ρ_SR，s分别表示积分结束和起始时刻对应的单程几何距离，RLT_e，RLT_s表示积分结束和 起始时刻引力弯曲，TDB表示太阳系质心动力学时，t(τ)表示对应的原时，[]_e，[]_s分别 对应积分结束和起始时刻；

其中：

${\rho }_{\mathrm{SR},e}-{\rho }_{\mathrm{SR},s}\approx -(\frac{{\partial \rho }_{\mathrm{SR},e}}{\partial {\overrightarrow{r}}_{\mathrm{SR},e}}\Delta {\overrightarrow{r}}_{\mathrm{SR}}+\frac{1}{2}{\Delta }^T{\overrightarrow{r}}_{\mathrm{SR}}\frac{{\partial }^2{\rho }_{\mathrm{SR},e}}{\partial {\overrightarrow{r}}_{\mathrm{SR},e}^2}\Delta {\overrightarrow{r}}_{\mathrm{SR}})$

$\Delta {\overrightarrow{r}}_{\mathrm{SR}}={[{\overrightarrow{r}}_{S/C}\left({\mathrm{TDB}}_{\mathrm{SC}}\right)-{\overrightarrow{r}}_R\left({\mathrm{TDB}}_R\right)]}_s-{[{\overrightarrow{r}}_{S/C}\left({\mathrm{TDB}}_{\mathrm{SC}}\right)-{\overrightarrow{r}}_R\left({\mathrm{TDB}}_R\right)]}_e$

$\frac{{\partial \rho }_{\mathrm{SR},e}}{\partial {\overrightarrow{r}}_{\mathrm{SR},e}}=\frac{{\overrightarrow{r}}_{\mathrm{SR},e}}{{\rho }_{\mathrm{SR},e}}$

$\frac{{\partial }^2{\rho }_{\mathrm{SR},e}}{\partial {\overrightarrow{r}}_{\mathrm{SR},e}^2}=-\frac{1}{{\rho }_{\mathrm{SR},e}}\left(\begin{array}{ccc}\frac{X_{\mathrm{SR},e}^2}{{\rho }_{\mathrm{SR},e}^2}-1& \frac{X_{\mathrm{SR},e}{Y}_{\mathrm{SR},e}}{{\rho }_{\mathrm{SR},e}^2}& \frac{X_{\mathrm{SR},e}{Z}_{\mathrm{SR},e}}{{\rho }_{\mathrm{SR},e}^2}\\ \frac{X_{\mathrm{SR},e}{Y}_{\mathrm{SR},e}}{{\rho }_{\mathrm{SR},e}^2}& \frac{Y_{\mathrm{SR},e}^2}{{\rho }_{\mathrm{SR},e}^2}-1& \frac{Y_{\mathrm{SR},e}{Z}_{\mathrm{SR},e}}{{\rho }_{\mathrm{SR},e}^2}\\ \frac{X_{\mathrm{SR},e}{Z}_{\mathrm{SR},e}}{{\rho }_{\mathrm{SR},e}^2}& \frac{Y_{\mathrm{SR},e}{Z}_{\mathrm{SR},e}}{{\rho }_{\mathrm{SR},e}^2}& \frac{Z_{\mathrm{SR},e}^2}{{\rho }_{\mathrm{SR},e}^2}-1\end{array}\right)$

其中，X，Y，Z分别表示飞行器在广义相对论太阳系质心坐标系下的位置分量，下标S 表示发射端，R表示接收端，s表示积分开始时刻，e表示积分结束时刻。

本发明通过泰勒级数展开来计算积分周期始末两个时刻的光行时差，避免了直接使用 差分距离观测量获取Doppler、损失计算精度的缺点，消除了在双精度浮点定义时因字长 不足所产生的舍入误差对计算精度的影响，达到四精度浮点数计算的精度，保证了计算效 率。

附图说明

图1火星探测器使用常规方法与通过本发明计算的Doppler进行比较结果。

具体实施方式

本发明提供的双精度浮点定义下的深空Doppler测速计算方法包括以下步骤：

1.依据Doppler测量原理建立Doppler测量方程

积分Doppler的观测量是地面站接收信号频率相对标准频率的偏移量f_D，通过测量一 个积分周期内的周计数得到，即

$f_D=f_R-f_0=\frac{N_C}{\mathrm{\Delta T}}+\mathrm{\delta f}+v---\left(1\right)$

其中f_R表示接收端接收信号频率，f₀表示地面接收站本振频率，N_C为一个积分周期内的 周计数，ΔT为一个积分周期，δf为观测系统误差，v为观测随机误差。

地面接收站的一个积分周期ΔT内累积的Doppler周计数N_C的变化可以推导出 Doppler观测量f_D。对于给定的地面测站，连续的Doppler观测有连续的积分周期。积分 周期可以短到0.1s，也可以长到半天(43200s)，典型的积分周期是几秒到几千秒的区间。

对于N_C有

$N_C={\int }_{t_s}^{t_e}(f_R-f_0)\mathrm{dt}={\int }_{t_s}^{t_e}{f}_R\mathrm{dt}-f_0\mathrm{\Delta T}---\left(2\right)$

$={\int }_{{\tau }_s}^{{\tau }_e}{f}_S\mathrm{d\tau }-f_0\mathrm{\Delta T}$

其中f₀为地面接收站的本振频率，小写的下标s，e表示积分起始和结束，t表示地面测站 时间，一般采用协调世界时UTC时间系统，τ表示信号发射端的原时。

考虑比较理想的状态，即发射标准频率f_S是一个常数，则有：

f_RΔT＝f_S(τ_e-τ_s)＝f_SΔτ                           (3)

为了方便区分，ΔT为接收端一个积分周期，Δτ为接收端一个积分周期对应的发射端的原 时差。

2.用泰勒级数展开式来计算积分周期始末时刻的光行时差

考虑到深空探测中光行时一般在广义相对论太阳系质心坐标系下建立，原时差可由下 公式计算：

$\mathrm{\Delta \tau }={\tau }_{S,e}-{\tau }_{S,s}$

$=\frac{({\rho }_{\mathrm{SR},e}-{\rho }_{\mathrm{SR},s})}{c}+{\mathrm{RLT}}_e-{\mathrm{RLT}}_s$$\left(4\right)$

$-\{{[{\mathrm{TDB}}_R-t\left({\tau }_R\right)]}_e-{[{\mathrm{TDB}}_R-t\left({\tau }_R\right)]}_s\}$

$+\{{[{\mathrm{TDB}}_S-t\left({\tau }_S\right)]}_e-{[{\mathrm{TDB}}_S-t\left({\tau }_S\right)]}_s\}$

其中，τ_S，e，τ_S，s分别表示星上信号发射端积分周期结束和起始时刻对应的原始，ρ_SR，e，ρ_SR，s分 别表示积分结束和起始时刻对应的单程几何距离，RLT_e，RLT_s表示积分结束和起始时刻引 力弯曲，TDB表示太阳系质心动力学时，t(τ)表示对应的原时，[]_e，[]_s分别对应积分结 束和起始时刻。

对于深空探测(如火星探测)而言，(4)式右端第1项的数值约在200-1300s之间， 第2、3项为相对论效应项，数值约为10^-6s量级，而第4、5项的数值在10^-9s量级。因此， 在双精度浮点定义下，按照式(4)计算光行时，右端第一项引入的舍入误差约为3×10^-12s， 对应8.4GHz的X波段的Doppler测量误差为1mm/s(1s积分周期)，远低于当前的测量精 度。因此，为了满足双精度浮点系统的精度要求，需将上式右端第一项计算改为以下：

${\rho }_{\mathrm{SR},e}-{\rho }_{\mathrm{SR},s}\approx -(\frac{{\partial \rho }_{\mathrm{SR},e}}{\partial {\overrightarrow{r}}_{\mathrm{SR},e}}\Delta {\overrightarrow{r}}_{\mathrm{SR}}+\frac{1}{2}{\Delta }^T{\overrightarrow{r}}_{\mathrm{SR}}\frac{{\partial }^2{\rho }_{\mathrm{SR},e}}{\partial {\overrightarrow{r}}_{\mathrm{SR},e}^2}\Delta {\overrightarrow{r}}_{\mathrm{SR}})---\left(5\right)$

这里，

$\Delta {\overrightarrow{r}}_{\mathrm{SR}}={[{\overrightarrow{r}}_{S/C}\left({\mathrm{TDB}}_{\mathrm{SC}}\right)-{\overrightarrow{r}}_R\left({\mathrm{TDB}}_R\right)]}_s-{[{\overrightarrow{r}}_{S/C}\left({\mathrm{TDB}}_{\mathrm{SC}}\right)-{\overrightarrow{r}}_R\left({\mathrm{TDB}}_R\right)]}_e$

$\frac{{\partial \rho }_{\mathrm{SR},e}}{\partial {\overrightarrow{r}}_{\mathrm{SR},e}}=\frac{{\overrightarrow{r}}_{\mathrm{SR},e}}{{\rho }_{\mathrm{SR},e}}$

$\frac{{\partial }^2{\rho }_{\mathrm{SR},e}}{\partial {\overrightarrow{r}}_{\mathrm{SR},e}^2}=-\frac{1}{{\rho }_{\mathrm{SR},e}}\left(\begin{array}{ccc}\frac{X_{\mathrm{SR},e}^2}{{\rho }_{\mathrm{SR},e}^2}-1& \frac{X_{\mathrm{SR},e}{Y}_{\mathrm{SR},e}}{{\rho }_{\mathrm{SR},e}^2}& \frac{X_{\mathrm{SR},e}{Z}_{\mathrm{SR},e}}{{\rho }_{\mathrm{SR},e}^2}\\ \frac{X_{\mathrm{SR},e}{Y}_{\mathrm{SR},e}}{{\rho }_{\mathrm{SR},e}^2}& \frac{Y_{\mathrm{SR},e}^2}{{\rho }_{\mathrm{SR},e}^2}-1& \frac{Y_{\mathrm{SR},e}{Z}_{\mathrm{SR},e}}{{\rho }_{\mathrm{SR},e}^2}\\ \frac{X_{\mathrm{SR},e}{Z}_{\mathrm{SR},e}}{{\rho }_{\mathrm{SR},e}^2}& \frac{Y_{\mathrm{SR},e}{Z}_{\mathrm{SR},e}}{{\rho }_{\mathrm{SR},e}^2}& \frac{Z_{\mathrm{SR},e}^2}{{\rho }_{\mathrm{SR},e}^2}-1\end{array}\right)$

其中，X，Y，Z分别表示飞行器在广义相对论太阳系质心坐标系下的位置分量，下标S表示 发射端，R表示接收端，s表示积分开始时刻，e表示积分结束时刻。

如果Doppler积分周期相对较短(＜100s)，可以忽略式中相对论修正项的影响，由此引 入的误差不会超过0.3mm，可以优于1ps的计算精度要求。

3.将光行时差代入测量方程计算Doppler测速

将上述(4)-(5)式代入公式(1)即可得到深空Doppler测量值。该计算方法对差分单程 距离的计算使用了泰勒展开，避开了直接使用差分距离观测量获取Doppler，损失计算精 度的缺点，实现了使用双精度浮点数进行高精度的Doppler计算，保证了计算精度，提高 了计算效率。

本发明已经成功应用于火星探测器的轨道计算中。图1为使用双精度常规方法进行计 算与使用本发明计算的残差(观测值与计算值的差)比较。图1上半部分为常规方法计算结 果，下半部分为使用该专利方法的计算结果，可以看到使用该方法计算的Doppler残差基 本为白噪声，Doppler的计算值精度提高一个量级。