考虑光行时校正的高精度全球导航卫星系统测速算法

摘要

全球导航卫星系统（GNSS）测速算法认为多普勒测量值为接收机与GNSS卫星的相对运动速度在视线方向上的投影，未考虑光行时的影响，存在一个微小的系统性偏差，测速精度不高，限制了其在精密测速领域中的应用。在传统的GNSS测速算法基础上，给出了一至三阶的积分多普勒插值算法，用瞬时多普勒替代了平均多普勒，对于低轨道地球卫星（LEO）用户，使用二阶插值积分多普勒的测速精度较使用平均多普勒的测速精度提升了三个数量级；另一方面考虑了光行时的影响，修正了因收发时刻不一致引入的多普勒偏差，并分别给出了地球惯性坐标系（ECI）和地球固连坐标系（ECEF）下的修正算法，实验结果表明，光行时修正算法降低了系统性偏差，三维速度误差的均方根（RMS）降低至3.3 mm·s^-1，能够满足高精度测速的应用需求。

关键词

全球导航卫星系统; 光行时; 多普勒测速; 积分多普勒; 拉格朗日插值

速度信息是描述物体运动状态的重要参数之一，全球导航卫星系统（ global navigation satellite system，GNSS） 测速技术具有高精度、高效率、简单方便等特点，这使其在航空飞行、智能导航、无人驾驶、海上交通等国民经济领域得到了广泛应用^［1］。

利用GNSS 进行测速主要有4种方式^［2-3］：①是利用定位结果通过位置差分的方法获取速度，即实时差分定位（real-time kinematic， RTK）；②是利用伪距观测值差分计算速度^［4］；③是利用原始多普勒观测值计算速度^［5-6］；④是利用积分多普勒即载波相位差分计算速度^［1］。近年来，众多学者对GNSS测速方法进行了深入研究：文献［7］分析了位置差分、多普勒和积分多普勒三种测速方法的精度与载体运动状态的关系，当全球定位系统（global positioning system， GPS）观测值的测量频率为1 Hz，位置差分方法在载体速度变化较大时测速精度较低；文献［4，8］提出了利用单频伪距差分测定载体速度的方法，通过历元间伪距差分可求得$\mathrm{\Delta }t$时间内的平均速度，伪距差分方法克服了积分多普勒方法中的周跳问题，但受伪距测量精度的影响，测速误差较大；文献［9］分析了多普勒测速法和位置差分测速法各自的优缺点，当载体做匀速运动时位置差分的精度较高，当载体速度变化较大时多普勒测速的精度较高。积分多普勒方法是对相邻历元间的载波相位观测值进行差分^［5-6，10］，它反映的是接收机在差分时间段内的平均多普勒^［11］。当接收机处于匀速运动或静止情况下，积分多普勒测速精度最高；但在非匀速运动时，积分多普勒测速精度会明显恶化。

为了将积分多普勒方法推广至载体的非匀速运动，本文提出了积分多普勒插值方法，采用拉格朗日（Lagrange）多项式插值方法^［12］，将积分多普勒测量值内插到测速历元时刻，测速结果即为历元时刻的瞬时速度，本文对积分多普勒的插值精度进行了评估分析，给出积分多普勒插值方法的使用建议。

此外，经过对测速结果进行分析，发现存在一个微小的系统性速度偏差，这是因为传统的GNSS测速算法认为影响多普勒测量值的主要因素包括：GNSS卫星和接收机的相对运动速度在视线方向上的投影、GNSS卫星钟速、接收机钟速，并未考虑光行时对多普勒测量值的影响^［3，5］。本文在传统的GNSS测速算法基础上，考虑了光行时的影响，修正了因收发时刻不一致引入的多普勒偏差，并给出了地球惯性坐标系（earth-centered inertial coordinate system， ECI）和地球固连坐标系（earth-centered earth-fixed coordinate system， ECEF）下的修正算法。

1 GNSS精密测速模型

1.1　积分多普勒插值算法

本文中积分多普勒的测量频度为1 Hz，即每秒1个测量值，设在$t_i$时刻的积分多普勒测量值为$f_i$，$f_i$为1 s内的平均多普勒，在1 s的时间尺度内，假设空间运动遵照匀加速直线运动的条件，即在每秒内，加速度为一定值，$f_i$可近似于$t_i-\frac{\mathrm{1}}{\mathrm{2}}$时刻的瞬时多普勒，则积分多普勒插值归结为以下数学问题：

在平面上有$(t_{\mathrm{1}}-\frac{\mathrm{1}}{\mathrm{2}},f_{\mathrm{1}})$，$(t_{\mathrm{2}}-\frac{\mathrm{1}}{\mathrm{2}},f_{\mathrm{2}})$，$(t_i-\frac{\mathrm{1}}{\mathrm{2}},f_i)$，…，$(t_n-\frac{\mathrm{1}}{\mathrm{2}},f_n)$共n个点，求插值函数$g(t)$，使其经过这n个点，这可通过拉格朗日插值方法进行求解，具体算法如下。

设$D_n=\left\{\mathrm{1,2},\cdots ,\mathrm{n}\right\}$，作n个多项式$p_j(t)$，其中$j\in D_n$，有：

$$p_j(t)=\prod _{i\in B_j}\frac{t-(t_i-\frac{\mathrm{1}}{\mathrm{2}})}{(t_j-\frac{\mathrm{1}}{\mathrm{2}})-(t_i-\frac{\mathrm{1}}{\mathrm{2}})}$$

（1）

式中：$p_j(t)$为n-1次多项式，$B_j=\{i\left|i\ne j,i\in D_n\right.\}$，且满足$\forall m\in B_j$，有$p_j(t_m-\frac{\mathrm{1}}{\mathrm{2}})=\mathrm{0}$，并且$p_j(t_j-\frac{\mathrm{1}}{\mathrm{2}})=\mathrm{1}$。

最后可得拉格朗日插值公式^［12］：

$$g(t)={\displaystyle \sum _{j=\mathrm{1}}^n}{f}_j{p}_j(t)$$

（2）

（1）当n=3时，二阶拉格朗日插值公式可写为

$$\begin{array}{l}{g}^{(\mathrm{2})}(t)=\frac{[t-(t_{\mathrm{2}}-\frac{\mathrm{1}}{\mathrm{2}})][t-(t_{\mathrm{3}}-\frac{\mathrm{1}}{\mathrm{2}})]}{[(t_{\mathrm{1}}-\frac{\mathrm{1}}{\mathrm{2}})-(t_{\mathrm{2}}-\frac{\mathrm{1}}{\mathrm{2}})][(t_{\mathrm{1}}-\frac{\mathrm{1}}{\mathrm{2}})-(t_{\mathrm{3}}-\frac{\mathrm{1}}{\mathrm{2}})]}{f}_{\mathrm{1}}+\\ \frac{[t-(t_{\mathrm{1}}-\frac{\mathrm{1}}{\mathrm{2}})][t-(t_{\mathrm{3}}-\frac{\mathrm{1}}{\mathrm{2}})]}{[(t_{\mathrm{2}}-\frac{\mathrm{1}}{\mathrm{2}})-(t_{\mathrm{1}}-\frac{\mathrm{1}}{\mathrm{2}})][(t_{\mathrm{2}}-\frac{\mathrm{1}}{\mathrm{2}})-(t_{\mathrm{3}}-\frac{\mathrm{1}}{\mathrm{2}})]}{f}_{\mathrm{2}}+\\ \frac{[t-(t_{\mathrm{1}}-\frac{\mathrm{1}}{\mathrm{2}})][t-(t_{\mathrm{2}}-\frac{\mathrm{1}}{\mathrm{2}})]}{[(t_{\mathrm{3}}-\frac{\mathrm{1}}{\mathrm{2}})-(t_{\mathrm{1}}-\frac{\mathrm{1}}{\mathrm{2}})][(t_{\mathrm{3}}-\frac{\mathrm{1}}{\mathrm{2}})-(t_{\mathrm{2}}-\frac{\mathrm{1}}{\mathrm{2}})]}{f}_{\mathrm{3}}\end{array}$$

（3）

式中：$t_{\mathrm{1}}$、$t_{\mathrm{2}}$、$t_{\mathrm{3}}$依次间隔1 s，当$t=t_{\mathrm{3}}$时，式（3）可简化为

$$g^{(\mathrm{2})}(t)=\frac{\mathrm{3}}{\mathrm{8}}{f}_{\mathrm{1}}-\frac{\mathrm{5}}{\mathrm{4}}{f}_{\mathrm{2}}+\frac{\mathrm{15}}{\mathrm{8}}{f}_{\mathrm{3}}$$

（4）

（2）当n=2时，不加推导，直接给出一阶拉格朗日插值公式：

$$g^{(\mathrm{1})}(t)=-\frac{\mathrm{1}}{\mathrm{2}}{f}_{\mathrm{1}}+\frac{\mathrm{3}}{\mathrm{2}}{f}_{\mathrm{2}}$$

（5）

（3）同理，当n=4时，三阶拉格朗日插值公式为：

$$g^{(\mathrm{3})}(t)=-\frac{\mathrm{5}}{\mathrm{16}}{f}_{\mathrm{1}}+\frac{\mathrm{21}}{\mathrm{16}}{f}_{\mathrm{2}}-\frac{\mathrm{35}}{\mathrm{16}}{f}_{\mathrm{3}}+\frac{\mathrm{35}}{\mathrm{16}}{f}_{\mathrm{4}}$$

（6）

以上给出了一至三阶的积分多普勒插值公式，将积分多普勒测量值内插到测速历元时刻，消除因测量时刻不准而引入的速度偏差。

不失一般性，若需要更高阶的积分多普勒插值公式，则参照上述一至三阶插值公式的推导过程，对公式（2）进行变换，得到对应的高阶积分多普勒插值公式。

1.2　光行时校正算法

GNSS信号在空中传播的测距方程可以表达为式（7）：^［13］

$$\begin{array}{l}\rho ={\rho }_{\mathrm{r}}^{\mathrm{s}}(t)+c\mathrm{\delta }{t}_{\mathrm{r}}(t)-c\mathrm{\delta }{T}_{\mathrm{s}}(t-\frac{\rho }{c})+\mathrm{\delta }{\rho }_{\mathrm{i}\mathrm{o}\mathrm{n}}+\\ \mathrm{\delta }{\rho }_{\mathrm{t}\mathrm{r}\mathrm{o}\mathrm{p}}+{\epsilon }_{\mathrm{\rho }}(t)\end{array}$$

（7）

式中：$\rho$为接收机测距值；${\rho }_{\mathrm{r}}^{\mathrm{s}}(t)$为星地真实距离；$c$为光速；$\mathrm{\delta }{T}_{\mathrm{s}}(t-\frac{\rho }{c})$和$\mathrm{\delta }{t}_{\mathrm{r}}(t)$分别为卫星在发射时刻和接收机在接收时刻的钟差；$\mathrm{\delta }{\rho }_{\mathrm{i}\mathrm{o}\mathrm{n}}$和$\mathrm{\delta }{\rho }_{\mathrm{t}\mathrm{r}\mathrm{o}\mathrm{p}}$分别为电离层和对流层的延迟；${\epsilon }_{\mathrm{\rho }}(t)$是测距噪声。

下面分别给出了ECEF和ECI坐标系下考虑光行时校正的GNSS测速算法。

1.2.1　ECEF坐标系校正算法

在ECEF地心地固坐标系中，星地距离${\rho }_{\mathrm{r}}^{\mathrm{s}}(t)$可表示为^［14］

$${\rho }_{\mathrm{r}}^{\mathrm{s}}(t)=\left|{\mathit{C}}_{\mathrm{e}}^{\mathrm{I}}(t-\frac{\rho }{c}){\mathit{r}}_{\mathrm{e}}^{\mathrm{s}}(t-\frac{\rho }{c})-{\mathit{C}}_{\mathrm{e}}^{\mathrm{I}}(t){\mathit{r}}_{\mathrm{e}}^{\mathrm{r}}(t)\right|$$

（8）

式中：${\mathit{r}}_{\mathrm{e}}^{\mathrm{s}}(t-\frac{\rho }{c})$和${\mathit{r}}_{\mathrm{e}}^{\mathrm{r}}(t)$分别为GNSS卫星在发射时刻的ECEF坐标系和接收机在接收时刻的ECEF坐标系的位置矢量；${\mathit{C}}_{\mathrm{e}}^{\mathrm{I}}(t-\frac{\rho }{c})$和${\mathit{C}}_{\mathrm{e}}^{\mathrm{I}}(t)$分别为GNSS卫星和接收机的ECEF坐标系至ECI坐标系的坐标旋转矩阵。

将式（8）代入式（7）中，对时间求导，并忽略电离层和对流层的变化率对GNSS测速的影响，可得：

$$\begin{array}{l}\dot{\rho }={\mathit{\mu }}_{\mathrm{a}\mathrm{s}}^{\mathrm{I}}\{(\mathrm{1}-\frac{\dot{\rho }}{c})[{\dot{\mathit{C}}}_{\mathrm{e}}^{\mathrm{I}}(t-\frac{\rho }{c}){\mathit{r}}_{\mathrm{e}}^{\mathrm{s}}(t-\frac{\rho }{c})+\\ {\mathit{C}}_{\mathrm{e}}^{\mathrm{I}}(t-\frac{\rho }{c}){\mathit{v}}_{\mathrm{e}}^{\mathrm{s}}(t-\frac{\rho }{c})]-[{\dot{\mathit{C}}}_{\mathrm{e}}^{\mathrm{I}}(t){\mathit{r}}_{\mathrm{e}}^{\mathrm{r}}(t)+{\mathit{C}}_{\mathrm{e}}^{\mathrm{I}}(t){\mathit{v}}_{\mathrm{e}}^{\mathrm{r}}(t)]\}+\\ c\mathrm{\delta }{\dot{t}}_r(t)-c\mathrm{\delta }{\dot{T}}_{\mathrm{s}}(t-\frac{\rho }{c})(\mathrm{1}-\frac{\dot{\rho }}{c})+{\dot{\epsilon }}_{\mathrm{\rho }}(t)\end{array}$$

（9）

式中：${\mathit{\mu }}_{\mathrm{a}\mathrm{s}}^{\mathrm{I}}$为ECI惯性系下距离矢量的方向余弦函数；$\dot{\rho }$为距离变化率测量值；${\dot{\mathit{C}}}_{\mathrm{e}}^{\mathrm{I}}(t-\frac{\rho }{c})$和${\dot{\mathit{C}}}_{\mathrm{e}}^{\mathrm{I}}(t)$分别为GNSS卫星和接收机的ECEF至ECI的坐标旋转矩阵的变化率；${\mathit{v}}_{\mathrm{e}}^{\mathrm{s}}(t-\frac{\rho }{c})$和${\mathit{v}}_{\mathrm{e}}^{\mathrm{r}}(t)$分别为GNSS卫星在发射时刻的ECEF坐标系和接收机在接收时刻的ECEF坐标系的速度矢量；$\delta {\dot{T}}_{\mathrm{s}}(t-\frac{\rho }{c})$和$\delta {\dot{t}}_{\mathrm{r}}(t)$分别为卫星在发射时刻的和接收机在接收时刻的钟速；${\dot{\epsilon }}_{\mathrm{\rho }}(t)$是距离变化率噪声。

将速度矢量投影在视线方向上，等式两边都有$\dot{\rho }$，将式（9）整理后得：

$$\begin{array}{l}\dot{\rho }\left(\mathrm{1}+\frac{{\mathit{\mu }}_{\mathrm{a}\mathrm{s}}^{\mathrm{I}}}{c}[{\dot{\mathit{C}}}_{\mathrm{e}}^{\mathrm{I}}(t-\frac{\rho }{c}){\mathit{r}}_{\mathrm{e}}^{\mathrm{s}}(t-\frac{\rho }{c})+{\mathit{C}}_{\mathrm{e}}^{\mathrm{I}}(t-\frac{\rho }{c}){\mathit{v}}_{\mathrm{e}}^{\mathrm{s}}(t-\frac{\rho }{c})]-\delta {\dot{T}}_{\mathrm{s}}(t-\frac{\rho }{c})\right)=\\ {\mathit{\mu }}_{\mathrm{a}\mathrm{s}}^{\mathrm{I}}\{{\dot{\mathit{C}}}_{\mathrm{e}}^{\mathrm{I}}(t-\frac{\rho }{c}){\mathit{r}}_{\mathrm{e}}^{\mathrm{s}}(t-\frac{\rho }{c})+{\mathit{C}}_{\mathrm{e}}^{\mathrm{I}}(t-\frac{\rho }{c}){\mathit{v}}_{\mathrm{e}}^{\mathrm{s}}(t-\frac{\rho }{c})-\\ [{\dot{\mathit{C}}}_{\mathrm{e}}^{\mathrm{I}}(t){\mathit{r}}_{\mathrm{e}}^{\mathrm{r}}(t)+{\mathit{C}}_{\mathrm{e}}^{\mathrm{I}}(t){\mathit{v}}_{\mathrm{e}}^{\mathrm{r}}(t)]\}+c\mathrm{\delta }{\dot{t}}_{\mathrm{r}}(t)-c\mathrm{\delta }{\dot{T}}_{\mathrm{s}}(t-\frac{\rho }{c})+{\dot{\epsilon }}_{\mathrm{\rho }}(t)\end{array}$$

（10）

为了得到ECEF坐标系下的测速方程，设定ECI惯性坐标系与接收时刻的ECEF坐标系重合，即，${\mathit{C}}_{\mathrm{e}}^{\mathrm{I}}(t)=\mathit{I}$，并将${\dot{\mathit{C}}}_{\mathrm{\epsilon }}^{\mathrm{{\rm I}}}(t-\frac{\rho }{c})={\mathit{C}}_{\mathrm{e}}^{\mathrm{I}}(t-\frac{\rho }{c}){\mathit{\Omega }}_{\mathrm{i}\mathrm{e}}^{\mathrm{e}}$和${\dot{\mathit{C}}}_{\mathrm{e}}^{\mathrm{I}}(t)={\mathit{C}}_{\mathrm{e}}^{\mathrm{I}}(t){\mathit{\Omega }}_{\mathrm{i}\mathrm{e}}^{\mathrm{e}}$代入公式（10）中，则式（10）可以简化为公式（11）：^［14］

$$\begin{array}{l}\dot{\rho }\left(\mathrm{1}+\frac{{\mathit{\mu }}_{\mathrm{a}\mathrm{s}}^{\mathrm{I}}}{c}[{\mathit{C}}_{\mathrm{e}}^{\mathrm{I}}(t-\frac{\rho }{c})({\mathit{\Omega }}_{\mathrm{i}\mathrm{e}}^{\mathrm{e}}{\mathit{r}}_{\mathrm{e}}^{\mathrm{s}}(t-\frac{\rho }{c})+{\mathit{v}}_{\mathrm{e}}^{\mathrm{s}}(t-\frac{\rho }{c}))]-\delta {\dot{T}}_{\mathrm{s}}(t-\frac{\rho }{c})\right)-\\ {\mathit{\mu }}_{\mathrm{a}\mathrm{s}}^{\mathrm{I}}[{\mathit{C}}_{\mathrm{e}}^{\mathrm{I}}(t-\frac{\rho }{c})({\mathit{\Omega }}_{\mathrm{i}\mathrm{e}}^{\mathrm{e}}{\mathit{r}}_{\mathrm{e}}^{\mathrm{s}}(t-\frac{\rho }{c})+{\mathit{v}}_{\mathrm{e}}^{\mathrm{s}}(t-\frac{\rho }{c}))-{\mathit{\Omega }}_{\mathrm{i}\mathrm{e}}^{\mathrm{e}}{\mathit{r}}_{\mathrm{e}}^{\mathrm{r}}(t)]+c\delta {\dot{T}}_{\mathrm{s}}(t-\frac{\rho }{c})=\\ -{\mathit{\mu }}_{\mathrm{a}\mathrm{s}}^{\mathrm{I}}{\mathit{v}}_{\mathrm{e}}^{\mathrm{r}}(t)+c\delta {\dot{t}}_{\mathrm{r}}(t)+{\dot{\epsilon }}_{\mathrm{\rho }}(t)\end{array}$$

（11）

式中：${\mathit{\Omega }}_{\mathrm{i}\mathrm{e}}^{\mathrm{e}}$为地球自转角速度矢量，令${\dot{\rho }}_{\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{r}}=\dot{\rho }\left(\mathrm{1}+\frac{{\mathit{\mu }}_{\mathrm{a}\mathrm{s}}^{\mathrm{I}}}{c}[{\mathit{C}}_{\mathrm{e}}^{\mathrm{I}}(t-\frac{\rho }{c})({\mathit{\Omega }}_{\mathrm{i}\mathrm{e}}^{\mathrm{e}}{\mathit{r}}_{\mathrm{e}}^{\mathrm{s}}(t-\frac{\rho }{c})+{\mathit{v}}_{\mathrm{e}}^{\mathrm{s}}(t-\frac{\rho }{c}))]-\delta {\dot{T}}_{\mathrm{s}}(t-\frac{\rho }{c})\right)$，在对距离变化率测量值$\dot{\rho }$进行光行时校正后，则测速方程与经典的GNSS测速方程在形式上完全一致，经典的GNSS测速算法在对公式（7）进行求导时，并未考虑$t-\frac{\rho }{c}$关于$t$的隐函数求导，即未考虑收发时刻不一致对测速方程的影响，这是引起光行时测速误差的本质原因。

令${\mathit{v}}_{\mathrm{e}}={\left[{\mathit{v}}_{\mathrm{e}}^{\mathrm{r}}{(t)}^{\mathrm{T}},c\delta {\dot{t}}_{\mathrm{r}}(t)\right]}^{\mathrm{T}}$为待估状态参数，通过对公式（11）进行线性化，并运用最小二乘方法^［15］，解算出接收机在ECEF坐标系下的速度和钟速，该速度即为考虑了光行时校正的ECEF坐标系下的速度。速度求解过程与经典GNSS测速方程的最小二乘求解方法完全相同，这里不做赘述。

1.2.2　ECI坐标系校正算法

在ECI惯性坐标系中，星地距离${\rho }_{\mathrm{r}}^{\mathrm{s}}(t)$可表示为

$${\rho }_{\mathrm{r}}^{\mathrm{s}}(t)=\left|{\mathit{r}}_{\mathrm{I}}^{\mathrm{s}}(t-\frac{\rho }{c})-{\mathit{r}}_{\mathrm{I}}^{\mathrm{r}}(t)\right|$$

（12）

式中：${\mathit{r}}_{\mathrm{I}}^{\mathrm{s}}(t-\frac{\rho }{c})$和${\mathit{r}}_{\mathrm{I}}^{\mathrm{r}}(t)$分别为ECI坐标系下卫星在发射时刻和接收机在接收时刻的位置矢量。

ECI坐标系下测速方程的推导过程通过类比ECEF坐标系的测速方程给出，假设发射时刻$t-\frac{\rho }{c}$的与接收时刻$t$的ECEF坐标系都是ECI坐标系，即${\mathit{C}}_{\mathrm{e}}^{\mathrm{I}}(t-\frac{\rho }{c})=\mathit{I}$，${\mathit{C}}_{\mathrm{e}}^{\mathrm{I}}(t)=\mathit{I}$，${\dot{\mathit{C}}}_{\mathrm{e}}^{\mathrm{I}}(t-\frac{\rho }{c})=\mathrm{0}$，${\dot{\mathit{C}}}_{\mathrm{e}}^{\mathrm{I}}(t)=\mathrm{0}$，则式（10）可以简化为

$$\begin{array}{l}\dot{\rho }\left[\mathrm{1}+\frac{{\mathit{\mu }}_{\mathrm{a}\mathrm{s}}^{\mathrm{I}}}{c}{\mathit{v}}_{\mathrm{I}}^{\mathrm{s}}(t-\frac{\rho }{c})-\mathrm{\delta }{\dot{T}}_{\mathrm{s}}(t-\frac{\rho }{c})\right]-{\mathit{\mu }}_{\mathrm{a}\mathrm{s}}^{\mathrm{I}}{\mathit{v}}_{\mathrm{I}}^{\mathrm{s}}(t-\frac{\rho }{c})+\\ c\mathrm{\delta }{\dot{T}}_{\mathrm{s}}(t-\frac{\rho }{c})=-{\mathit{\mu }}_{\mathrm{a}\mathrm{s}}^{\mathrm{I}}{\mathit{v}}_{\mathrm{I}}^{\mathrm{r}}(t)+c\mathrm{\delta }{\dot{t}}_{\mathrm{r}}(t)+{\dot{\epsilon }}_{\mathrm{\rho }}(t)(\mathrm{13})\end{array}$$

式中：${\mathit{\mu }}_{\mathrm{a}\mathrm{s}}^{\mathrm{I}}$为ECI坐标系下距离矢量的方向余弦函数；${\mathit{v}}_{\mathrm{I}}^{\mathrm{s}}(t-\frac{\rho }{c})$和${\mathit{v}}_{\mathrm{I}}^{\mathrm{r}}(t)$分别为ECI坐标系下GNSS卫星在发射时刻和接收机在接收时刻的速度矢量；令${\dot{\rho }}_{\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{r}}=\dot{\rho }\left[\mathrm{1}+\frac{{\mathit{\mu }}_{\mathrm{a}\mathrm{s}}^{\mathrm{I}}}{c}{\mathit{v}}_{\mathrm{I}}^{\mathrm{s}}(t-\frac{\rho }{c})-\delta {\dot{T}}_{\mathrm{s}}(t-\frac{\rho }{c})\right]$，在对距离变化率测量值$\dot{\rho }$进行光行时校正后，即修正了光行时引入的多普勒偏差。

式（13）即为ECI坐标系下的GNSS测速方程，可通过最小二乘方法求解${\mathit{v}}_{\mathrm{I}}={\left[{\mathit{v}}_{\mathrm{I}}^{\mathrm{r}}{(t)}^{\mathrm{T}},c\delta {\dot{t}}_{\mathrm{r}}(t)\right]}^{\mathrm{T}}$，该速度和钟差为考虑了光行时校正的ECI坐标系下的速度和钟差。

2 数据分析与测试验证

2.1　实验场景

为了对考虑光行时校正的GNSS测速算法进行验证，本文搭建了全实物仿真测试系统，在sprient GNSS信号模拟器中建立了鸿雁卫星运动场景，鸿雁为低地球轨道（low earth orbit， LEO）卫星， GNSS星座设置为单GPS系统，信号频率为GPS L1 C/A，鸿雁卫星的轨道高度为1 175 km，轨道倾角为86.5°，鸿雁卫星的平均运动速度为7.36 km·s^-1，平均向心加速度为7.38 m·s^-2。sprient模拟器播发出GPS信号，并被自研的星载GNSS接收机进行接收与捕跟处理，实验时长2 h，可见卫星数目8~12颗，观测值采样频率1 Hz，星载GNSS接收机输出GNSS星历、L1 C/A伪距及积分多普勒等测试数据，进行事后数据处理，并与sprient模拟器输出的积分多普勒真值、LEO卫星速度真值进行比对，评估各种算法的处理精度，在此基础上，将各种算法在星载GNSS接收机中进行实现，接收机在线运行时长3 h，测试验证各种算法的实时处理精度，整个处理过程如图1所示。

图1  GNSS测速算法仿真测试验证系统

Fig. 1  Simulation and test verification system of GNSS velocimetry algorithm

2.2　积分多普勒插值精度分析

为了考核一至三阶拉格朗日插值公式的插值精度，在第2.1节LEO运行场景下进行仿真分析，LEO卫星的速度动态大，仿真结论更具说服力。

积分多普勒插值精度评估的主要步骤如下：

（1）首先通过仿真得到不同GNSS卫星与LEO卫星在历元时刻$t_i$的真距$\rho (t_i)$和瞬时多普勒真值$f(t_i)$；

（2）利用$t_i$时刻的真距$\rho (t_i)$和$t_{i+\mathrm{1}}$时刻的真距$\rho (t_{i+\mathrm{1}})$通过历元间差分得到$t_{i+\mathrm{1}}$时刻的积分多普勒仿真结果$f_{i+\mathrm{1}}$，即：

$$f_{i+\mathrm{1}}=-\frac{\rho (t_{i+\mathrm{1}})-\rho (t_i)}{(t_{i+\mathrm{1}}-t_i)\lambda }$$

（14）

式中：$\lambda$为GNSS信号的波长，对于L1信号，其值约为0.190 3 m；

（3）利用一至三阶拉格朗日插值公式分别计算$t_{i+\mathrm{1}}$时刻的瞬时多普勒插值结果$\widehat{f}(t_{i+\mathrm{1}})$，并与真值$f(t_{i+\mathrm{1}})$进行比较，评估插值精度。

$$\delta f(t_{i+\mathrm{1}})=\widehat{f}(t_{i+\mathrm{1}})-f(t_{i+\mathrm{1}})$$

（15）

图2与图3分别给出了LEO卫星与GPS卫星之间的L1 C/A积分多普勒测量值不插值以及一至三阶拉格朗日插值方法的精度评估结果。表1为积分多普勒不同插值方法的多普勒精度评估结果。

图2  LEO卫星积分多普勒不插值的多普勒误差分析（不同线条代表不同的GPS卫星）

Fig. 2  Doppler error analysis of integral Doppler non interpolation for LEO satellite (different lines represent different GPS satellites)

图3  LEO卫星积分多普勒一至三阶插值方法的多普勒误差分析（不同线条代表不同的GPS卫星）

Fig. 3  Doppler error analysis of first - to third-order integrated Doppler interpolation method for LEO satellite (different lines represent different GPS satellites)

表 1  积分多普勒不同插值方法的多普勒精度评估结果

Tab. 1  Doppler accuracy evaluation results of different interpolation methods for integral Doppler ( Hz )

方法	积分多普勒精度分析
	均值	标准差	均方根
原始积分多普勒	12.440 9	8.169 7	14.883 6
一阶Lagrange插值	2.8329×10-4	0.012 7	0.012 7
二阶Lagrange插值	3.0311×10-5	0.001 7	0.001 7
三阶Lagrange插值	3.1101×10-5	0.001 9	0.001 9

在原始积分多普勒测量值精度一定的情况下，积分多普勒的插值精度直接决定了用户的GNSS测速精度，为了对不同类型用户的积分多普勒插值效果进行对比，本文对地球静止轨道（geostationary orbit， GEO）卫星与GPS卫星之间、地面静止站与GPS卫星之间的积分多普勒测量值不插值（图4）以及采用一至三阶拉格朗日插值方法的插值精度进行了分析，结果如表2所示。

图4  GEO卫星及地面站积分多普勒不插值的多普勒误差分析（不同线条代表不同的GPS卫星）

Fig. 4  Doppler error analysis of integral Doppler non-interpolation for GEO satellite and ground station (different lines represent different GPS satellites)

表 2  不同类型用户积分多普勒不同插值方法的多普勒插值精度分析

Tab. 2  Doppler accuracy evaluation results of different interpolation methods for integral Doppler of different users ( Hz )

方法	积分多普勒插值精度（均方根）
	LEO	GEO	地面站
原始积分多普勒	14.883 6	1.256 8	0.223 0
一阶Lagrange插值	0.012 7	0.005 4	0.001 1
二阶Lagrange插值	0.001 7	0.008 9	0.001 8
三阶Lagrange插值	0.001 9	0.014 7	0.003 0

由表2可知，对于不同类型的用户，积分多普勒插值方法与不插值相比，能够显著提高积分多普勒的测量精度，其中LEO卫星的提升量级最大，且有如下结论：

（1）LEO卫星的运动速度变化快，插值算法对动态的补偿效果最明显，认为加加速度保持不变或变化较小，因此采用二阶及以上的插值算法为佳，二阶插值算法与不插值相比，积分多普勒插值精度提升将近4个数量级。

（2）GEO卫星或者地面站的运动速度变化相对较小，插值算法对动态的补偿效果较弱，认为加速度保持不变或变化较小，因此采用一阶及以上的插值算法为佳，一阶插值算法与不插值相比，积分多普勒插值精度提升2个数量级以上。

在上述不同插值方法下积分多普勒插值误差分析（图2、图3）的基础上，对不同积分多普勒插值方法对LEO卫星GNSS测速精度的影响进行评估，分析结果如图5、表3及图6所示。

图5  积分多普勒不插值时测速误差分析

Fig. 5  Velocity error analysis of integral Doppler non-interpolation

表 3  积分多普勒不同插值方法的测速精度评估结果

Tab. 3  Velocity Accuracy evaluation results of different interpolation methods for integral Doppler ( mm·s^-1 )

方法	测速精度分析
	均值			标准差			均方根
	Vx	Vy	Vz	Vx	Vy	Vz	V3D
原始积分多普勒	41.8	409.6	-123.0	3 522.7	1 433.6	3 447.2	5 150.6
一阶Lagrange插值	-0.266 9	-0.107 6	0.002 5	3.0	1.6	3.2	4.7
二阶Lagrange插值	0.047 6	0.007 4	-0.004 2	2.0	1.4	2.3	3.3
三阶Lagrange插值	0.045 2	0.008 5	-0.008 5	3.0	2.2	3.5	5.1

图6  积分多普勒一至三阶插值方法时测速误差分析

Fig. 6  Velocity error analysis of first - to third-order integrated Doppler interpolation method

由以上分析可得出如下结论：

（1）当不进行积分多普勒插值时，积分多普勒精度很差，三维速度误差的均方根（root mean square， RMS）达到14.883 6 Hz，受此影响，三维速度误差RMS达到5.150 6 m·s^-1，充分说明原始积分多普勒测速方法并不适用于速度快速变化的运动载体。

（2）一至三阶拉格朗日插值方法对积分多普勒测量精度和测速精度都有改善效果，一阶拉格朗日插值方法受积分多普勒插值方法误差较大的影响，RMS为0.012 7Hz，测速误差存在一个系统性偏差，二阶、三阶插值方法能有效消除1阶插值方法中的系统性偏差，积分多普勒插值RMS精度达到mHz；

（3）二阶、三阶插值方法的测速精度RMS分别为3.3和5.1 mm·s^-1，与原始积分多普勒测速方法相比，性能提升了三个数量级，且二阶插值方法的测速精度与三阶插值方法的测速精度相当，这与第2.2节LEO卫星积分多普勒插值精度的提升密切相关，从算法的复杂性考虑，二阶插值方法工程实现更为简单，因此本文推荐使用二阶插值方法作为积分多普勒测速的插值算法。

2.3　测速误差分析

由1.2节的理论分析过程可知，无论是ECEF坐标系还是ECI坐标系，考虑光行时校正的多普勒测速算法，都是对距离变化率$\dot{\rho }$进行光行时修正，得到${\dot{\rho }}_{\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{r}}$，基于${\dot{\rho }}_{\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{r}}$的GNSS测速方程跟经典GNSS测速方程在形式上保持一致，方程求解方法也完全相同。

2.3.1　ECEF坐标系的测速结果

不考虑光行时校正时，距离变化率采用$\dot{\rho }$，本文称之为经典测速算法；当考虑光行时校正时，距离变化率采用${\dot{\rho }}_{\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{r}}$，称之为修正算法，图7和图8给出了ECEF坐标系下两种算法的GNSS测速结果。

图7  ECEF坐标系下经典算法测速误差分析

Fig. 7  Velocity error analysis of traditional classical algorithm in ECEF coordinate system

图8  ECEF坐标系下修正算法测速误差分析

Fig. 8  Velocity error analysis of modified algorithm in ECEF coordinate system

（1）经典算法（图7）。

（2）修正算法（图8）。

在ECEF坐标系下，LEO卫星与GPS卫星之间的L1 C/A距离变化率的光行时修正分析结果如图9所示（不同的GPS卫星）。

图9  ECEF坐标系下距离变化率的光行时修正量级分析（不同线条代表不同的GPS卫星）

Fig. 9  Correction analysis of light travel time for pseudo velocity in ECEF coordinate system (different lines represent different GPS satellites)

由以上分析结果可得如下结论：

（1）经典算法由于未考虑光行时的影响，测速结果存在非零均值的系统性误差，三维速度误差RMS为13.1 mm·s^-1，且测速误差会有间歇性的跳变，该跳变与几何精度因子（geometric dilution of precision， GDOP）的波动基本一致，表明该跳变是由新的GNSS卫星信号进入锁定或已有GNSS卫星信号失锁引起的；

（2）修正算法修正了光行时引入的多普勒偏差，消除了系统性速度误差，三维速度误差RMS大幅降低至约3.3 mm·s^-1，降幅达74.8 %；

（3）基于${\dot{\rho }}_{\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{r}}$的修正算法相比于基于$\dot{\rho }$的经典算法，两者的距离变化率差异最大值为0.025 m·s^-1，所有距离变化率差异的RMS值为0.009 8 m·s^-1，这是引起经典算法中测速系统误差的主要原因，也解释了测速误差在切星时发生跳变的原因。

2.3.2　ECI坐标系的测速结果

本文还对ECI坐标系下经典算法以及修正算法的测速精度进行了数值实验，并和ECEF坐标系的测速精度进行了对比，分析结果如表4所示，在ECI坐标系下，LEO卫星与GPS卫星之间的L1 C/A距离变化率的光行时修正量如图10所示。

表 4  ECI / ECEF坐标系下经典算法与修正算法的测速精度对比

Tab. 4  Comparison of velocity accuracy of traditional classical algorithm and modified algorithm in ECI/ECEF coordinate system ( mm·s^-1 )

算法类型	坐标系	测速精度分析
		均值			标准差			均方根
		Vx	Vy	Vz	Vx	Vy	Vz	V3D
经典算法	ECEF	-1.1	4.0	0.4	6.6	7.5	7.4	13.1
	ECI	-1.1	4.0	0.4	6.6	7.5	7.4	13.1
修正算法	ECEF	0	0	0.0	2.0	1.4	2.3	3.3
	ECI	0	0	0.0	2.0	1.4	2.3	3.3

图10  ECI坐标系下距离变化率光行时的修正量级分析（不同线条代表不同的GPS卫星）

Fig. 10  Correction analysis of light travel time for pseudo velocity in ECI coordinate system (different lines represent different GPS satellites)

由以上分析结果可知：

（1）ECI坐标系下经典算法和修正算法的测速精度分别与ECEF坐标系下的测速精度完全相同，由图10可以看出，ECI坐标系下距离变化率的修正量与ECEF坐标系下的修正量完全一致，这是ECI坐标系的测速精度与ECEF坐标系的测速精度相同的本质原因；

（2）光行时校正算法的精度与坐标系无关，只是在ECEF坐标系和ECI坐标系有不同的方程表达形式。

2.4　GNSS测速算法试验验证

本节将GNSS测速经典算法和修正算法在星载GNSS接收机中进行实现，对两种算法的实时处理性能进行试验验证，实测结果如图11所示。

图11  ECEF坐标系下经典算法和修正算法的测速误差实测结果

Fig. 11  Real-time processing results of velocity measurement error using classic algorithm and modified algorithm in ECEF coordinate system

从图11中可以看出，经典算法的实时测速结果存在非零均值的系统性误差，三维速度误差RMS为13.2 mm·s^-1，修正算法的实时测速结果有效消除了系统性速度误差，三维速度误差RMS大幅降低至约3.3 mm·s^-1，这与第2.3.1节事后数据处理方法的结论一致，充分验证了光行时修正算法的优异性能。

3 结 论

为了解决基于积分多普勒的经典GNSS测速算法存在厘米级系统性偏差的问题，本文在经典的GNSS测速算法基础上，提出了考虑光行时校正的GNSS测速算法，消除了系统性偏差，提升了测速精度，主要包括以下方面：

（1）本文给出了一至三阶积分多普勒的插值公式，采用拉格朗日插值方法，将积分多普勒测量值内插到测速历元时刻，用瞬时积分多普勒替代平均积分多普勒，提升了瞬时积分多普勒的精度。对于LEO卫星用户，以二阶插值方法的插值性能最优，RMS为0.001 7 Hz，本文拟推荐使用二阶插值方法作为LEO卫星积分多普勒测速的插值算法；

（2）本文提出了考虑光行时校正的GNSS测速算法，并分别给出了地球惯性坐标系（ECI）和地球固连坐标系（ECEF）下的光行时修正算法，通过对距离变化率进行修正，在方程形式上与经典测速算法保持了一致，实验结果表明，光行时校正算法消除了经典测速算法的系统性偏差，修正后的三维速度精度可达毫米级，且ECEF坐标系下的测速精度与ECI坐标系下的测速精度完全相同。

作者贡献声明

韩星远：提出数学模型，算法设计，论文撰写。

吴 昊：数值仿真分析。

杨力强：文献调研，论文修改。

张园园：绘制图形，数据汇总。

贺一峰：指导论文撰写，修改论文框架。