CN103033827B - 一种卫星位置解算方法 - Google Patents

Abstract

本发明公开了一种卫星位置解算方法，包括：确定观测时刻t_e，从导航电文中获取参考时刻t_b的卫星位置参数，包括位置X₀＝{x₀,y₀,z₀}，速度V₀＝{v_x0，v_y0，v_z0}，计算t_b加速度a₀＝{a_x,a_y,a_z}，设置积分步长h，0＜h＜200，且为整数；确定需要积分的时间长度t_L＝t_e-t_b，并计算需要积分的次数当N＞4，利用线性单步方法计算前4个步长时刻的位置X₀、X₁、X₂、X₃和速度V₀、V₁、V₂、V₃；根据X₀、X₁、X₂、X₃和速度V₀、V₁、V₂、V₃，计算剩余步长时刻的卫星位置、速度及加速度信息，得到观测时刻t_e的卫星位置。本发明每一步积分用到前四步的速度和位置信息，每步积分都是首先进行预测，然后进行校正，完成每步积分只需要计算两次加速度值；通过较多地利用前面几步的已知信息，实现计算量小、精度高的卫星位置解算。

Description

一种卫星位置解算方法

技术领域

本发明涉及卫星定位技术领域，特别是涉及一种卫星位置解算方法。

背景技术

研究GLONASS（GLOBAL NAVIGATION SATELLITE SYSTE，格洛纳斯）卫星系统定位，尤其是多系统组合定位有重要实用价值，如何根据GLONASS广播星历进行卫星位置计算是组合定位的重要一环。GLONASS每30min更换一组新的星历参数，采用数值积分的方法求解某一时刻t_e的卫星速度和位置时，要以最新星历的参考时刻t_b为参考时刻，并以该时刻卫星的坐标和速度为起算数据，对加速度和速度进行积分，其积分通式如（1-1）所示：

$\left\{\begin{array}{c}{V}_{\mathrm{te}}=V_b+{\∫}_{\mathrm{tb}}^{\mathrm{te}}\mathrm{adt}\\ X_{\mathrm{te}}=X_b+{\∫}_{\mathrm{tb}}^{\mathrm{te}}\mathrm{Vdt}\end{array}\right.---(1-1)$

其中，V、X表示速度和位置，a表示加速度。目前GLONASS卫星轨道数值积分方法主要有龙格-库塔法（以下简称R-K法），其中以四阶龙格-库塔法最为常用。但是，龙格-库塔法在计算某一时刻卫星位置和速度时只用到前一步的信息，为了保证定位精度，需要多次重新计算多个点处的函数值，计算量较大。

使用R-K方法每进行一步积分，即由t_i时刻的速度和位置（V_i，X_i）计算计算t_i+1时刻的速度和位置坐标（V_i+1X_i+1），都要按照公式（1-2）进行计算，其中h＝t_i+1-t_i为步长，a₁,a₂,a₃,a₄为加速度中间变量值，V₁，V₂，V₃，V₄为速度中间变量值，为日月摄动加速度。

$\left\{\begin{array}{c}{a}_1=f(t_i,X_i,V_i,\stackrel{\·\·}{X})\\ V_1=V_i\\ a_2=f(t_i+h/2,X_i+V_{1}h/2,V_i+a_{1}h/2,\stackrel{\·\·}{X})\\ V_2=V_i+a_{1}h/2\\ a_3=f(t_i+h/2,X_i+V_{2}h/2,V_i+a_{2}h/2,\stackrel{\·\·}{X})\\ V_3=V_i+a_{2}h/2\\ a_4=f(t_i+h,X_i+V_{3}h,V_i+a_{3}h,\stackrel{\·\·}{X})\\ V_4=V_i+a_{3}h\\ V_{i+1}=V_i+h(a_1+2a_2+2a_3+a_4)/6\\ X_{i+1}=X_i+h(V_1+2V_2+2V_3+V_4)/6\end{array}\right.---(1-2)$

使用R-K法每计算一步积分，即由t_i时刻的速度和位置（V_i，X_i）计算t_i+1时刻的速度和位置坐标（V_i+1X_i+1）时，只用到前一步的信息V_i和X_i，为了提高运算精度，重新计算了4次加速度值a₁,a₂,a₃,a₄，4次速度值V₁，V₂，V₃，V₄。由于进行轨道积分所花费的代价主要在于求解加速度，因此，使用R-K方法进行GLONASS轨道积分存在运算量较大的缺陷。

发明内容

本发明要解决的技术问题是提供一种卫星位置解算方法，用以解决现有技术使用R-K方法进行轨道积分存在运算量较大的问题。

为解决上述技术问题，本发明提供一种卫星位置解算方法，包括：

确定观测时刻t_e，从导航电文中获取参考时刻t_b的卫星位置参数，包括位置X₀＝{x₀,y₀,z₀}，速度V₀＝{v_x0，v_y0，v_z0}，计算t_b加速度a₀＝{a_x,a_y,a_z}，设置积分步长h，0＜h＜200，且为整数；确定需要积分的时间长度t_L＝t_e-t_b，并计算需要积分的次数

当N＞4，利用线性单步方法计算前4个步长时刻的位置X₀、X₁、X₂、X₃和速度V₀、V₁、V₂、V₃；

根据X₀、X₁、X₂、X₃和速度V₀、V₁、V₂、V₃，计算剩余步长时刻的卫星位置、速度及加速度信息，得到观测时刻t_e的卫星位置。

进一步，通过公式（2）计算t_b加速度a₀＝{a_x,a_y,a_z}；

$\left\{\begin{array}{c}{a}_x=-\frac{\mathrm{\μ}}{r^3}x-\frac{3}{2}{J}_0^2\·\frac{\mathrm{\μ}{a}_e^2}{r^5}\·x\·(1-\frac{{5z}^2}{r^2})+{\mathrm{\ω}}^{2}x+2\mathrm{\ω}\·v_y+\stackrel{\·\·}{x}\\ a_y=-\frac{\mathrm{\μ}}{r^3}y-\frac{3}{2}{J}_0^2\·\frac{\mathrm{\μ}{a}_e^2}{r^5}\·y\·(1-\frac{{5z}^2}{r^2})+{\mathrm{\ω}}^{2}y-2\mathrm{\ω}\·v_x+\stackrel{\·\·}{y}\\ a_z=-\frac{\mathrm{\μ}}{r^3}z-\frac{3}{2}{J}_0^2\·\frac{\mathrm{\μ}{a}_e^2}{r^5}\·z\·(3-\frac{{5z}^2}{r^2})+\stackrel{\·\·}{z}\end{array}\right.---\left(2\right)$

其中，地球引力常量μ＝398600.44km³/s²；卫星到地球质心的距离地球的长半轴a_e＝6378.136km；地球重力位第二带谐系数 地球自转速度ω＝0.00007292115rad/s；为日月摄动加速度在x、y、z三个方向的分量。

进一步，根据X₀、X₁、X₂、X₃和速度V₀、V₁、V₂、V₃，计算剩余步长时刻的卫星位置、速度及加速度信息，得到观测时刻t_e的卫星位置，具体包括：

按照公式（3）计算t_i+1时刻卫星速度和位置的预测值之后，按照公式（4）对预测值进行校正，得到t_i+1时刻卫星位置X_i+1＝{x_i+1，y_i+1,z_i+1}和速度V_i+1＝{v_x，i+1，v_y，i+1，v_z，i+1}；

$\left\{\begin{array}{c}\\ a_i=f(t_i,X_i,V_i,\stackrel{\·\·}{X})\\ {\stackrel{\‾}{V}}_{i+1}=V_i+\frac{h}{24}(55a_i-59a_{i-1}+37a_{i-2}-9a_{i-3})\\ {\stackrel{\‾}{X}}_{i+1}=X_i+\frac{h}{24}(55V_i-59V_{i-1}+37V_{i-2}-9V_{i-3})\end{array}\right.---\left(3\right)$

$\left\{\begin{array}{c}\\ a_{i+1}=f(t_i+h,{\stackrel{\‾}{X}}_{i+1},{\stackrel{\‾}{V}}_{i+1},\stackrel{\·\·}{X})\\ V_{i+1}=V_i+\frac{h}{720}(251a_{i+1}+646a_i-264a_{i-1}+106a_{i-2}-19a_{i-3})\\ X_{i+1}=X_i+\frac{h}{720}[251V_{i+1}+646V_i-264V_{i-1}+106V_{i-2}-19V_{i-3}]\end{array}\right.---\left(4\right)$

判断是否i＜N，如果是，则置X_i＝X_i+1,V_i＝V_i+1，积分计算下一个时刻的位置、速度；如果否得到t_N时刻的卫星位置、速度及加速度信息；

根据t_N时刻的卫星位置、速度及加速度信息，利用线性单步方法得到观测时刻t_e时刻的卫星位置。

进一步，当N≤4，利用线性单步方法，计算观测时刻的卫星位置。

本发明有益效果如下：

本发明每一步积分用到前四步的速度和位置信息，每步积分都是首先进行预测，然后进行校正，完成每步积分只需要计算两次加速度值；通过较多地利用前面几步的已知信息，实现计算量小、精度高的卫星位置解算。

附图说明

图1是本发明实施例中一种卫星位置解算方法的流程图。

具体实施方式

以下结合附图以及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不限定本发明。

如图1所示，本发明实施例涉及一种卫星位置解算方法，具体包括如下步骤：

a、确定观测时刻t_e，从导航电文中获取参考时刻t_b的卫星位置参数，包括位置X₀＝{x₀,y₀,z₀}，速度V₀＝{v_x0，v_y0，v_z0}，按照公式（2）计算该时刻加速度a₀＝{a_x,a_y,a_z}，设置积分步长h（0＜h＜200，且为整数）；确定需要积分的时间长度t_L＝t_e-t_b，并计算需要积分的次数

$\left\{\begin{array}{c}{a}_x=-\frac{\mathrm{\μ}}{r^3}x-\frac{3}{2}{J}_0^2\·\frac{\mathrm{\μ}{a}_e^2}{r^5}\·x\·(1-\frac{{5z}^2}{r^2})+{\mathrm{\ω}}^{2}x+2\mathrm{\ω}\·v_y+\stackrel{\·\·}{x}\\ a_y=-\frac{\mathrm{\μ}}{r^3}y-\frac{3}{2}{J}_0^2\·\frac{\mathrm{\μ}{a}_e^2}{r^5}\·y\·(1-\frac{{5z}^2}{r^2})+{\mathrm{\ω}}^{2}y-2\mathrm{\ω}\·v_x+\stackrel{\·\·}{y}\\ a_z=-\frac{\mathrm{\μ}}{r^3}z-\frac{3}{2}{J}_0^2\·\frac{\mathrm{\μ}{a}_e^2}{r^5}\·z\·(3-\frac{{5z}^2}{r^2})+\stackrel{\·\·}{z}\end{array}\right.---\left(2\right)$

其中，地球引力常量μ＝398600.44km³/s²；卫星到地球质心的距离地球的长半轴a_e＝6378.136km；地球重力位第二带谐系数 地球自转速度ω＝0.00007292115rad/s；为日月摄动加速度在x、y、z三个方向的分量。

b、若N≤4则用线性单步方法（如R-K法），计算观测时刻的卫星位置。若N＞4，前4个步长时刻的位置和速度X₀、X₁、X₂、X₃与V₀、V₁、V₂、V₃用线性单步方法（如R-K法）计算，其他值按照步骤c计算。

c、按照公式（3）计算t_i+1时刻卫星速度和位置的预测值之后按照公式（4）对预测值进行校正，得到t_i+1时刻卫星位置X_i+1＝{x_i+1，y_i+1,z_i+1}和速度V_i+1＝{v_x，i+1，v_y，i+1,v_z，i+1}。

$\left\{\begin{array}{c}\\ a_i=f(t_i,X_i,V_i,\stackrel{\·\·}{X})\\ {\stackrel{\‾}{V}}_{i+1}=V_i+\frac{h}{24}(55a_i-59a_{i-1}+37a_{i-2}-9a_{i-3})\\ {\stackrel{\‾}{X}}_{i+1}=X_i+\frac{h}{24}(55V_i-59V_{i-1}+37V_{i-2}-9V_{i-3})\end{array}\right.---\left(3\right)$

$\left\{\begin{array}{c}\\ a_{i+1}=f(t_i+h,{\stackrel{\‾}{X}}_{i+1},{\stackrel{\‾}{V}}_{i+1},\stackrel{\·\·}{X})\\ V_{i+1}=V_i+\frac{h}{720}(251a_{i+1}+646a_i-264a_{i-1}+106a_{i-2}-19a_{i-3})\\ X_{i+1}=X_i+\frac{h}{720}[251V_{i+1}+646V_i-264V_{i-1}+106V_{i-2}-19V_{i-3}]\end{array}\right.---\left(4\right)$

d、判断是否i＜N，如果是，则置X_i＝X_i+1,V_i＝V_i+1，积分计算下一个时刻的位置、速度，当i＝N-1时，得到接近t_e的t_N时刻的卫星位置；如果否，即i＝N，得到t_N时刻的卫星位置、速度及加速度信息；

e、根据t_N时刻的卫星位置、速度及加速度信息，用线性单步方法（如R-K法）得到观测时刻t_e时刻的卫星位置。

以下是采用本专利方法的具体实例：

选取2009年8月29日23:45至2009年8月31日16:45每间隔半小时共82个时刻，GLONASS 2号卫星的广播星历作为实验数据。

分别使用本发明实施例方法与R-K方法，以2009年8月29日23:45至2009年8月31日15:45每半小时间隔的参考时间向前轨道积分30分钟，得到2009年8月30日00:15至2009年8月31日16:15半小时间隔的共81个时刻的，GLONASS 2号卫星位置的X、Y、Z坐标，再将计算结果与广播星历里给出的卫星坐标进行比较，得到两者差值的统计结果，如表1所示（30min轨道积分结果与星历所给坐标之差的绝对值，单位/m）。

分别使用本发明实施例方法与R-K方法，以2009年8月29日23:45至2009年8月31日15:15每半小时间隔的参考时间向前轨道积分60分钟，得到2009年8月30日00:45至2009年8月31日16:15半小时间隔的共80个时刻的GLONASS 2号卫星位置的X、Y、Z坐标，再将计算结果与广播星历里给出的卫星坐标进行比较，得到两者差值的统计结果，如表2所示（60min轨道积分结果与星历所给坐标之差的绝对值，单位/m）。

表1

表2

（1）向前积分30min，选取不同步长对两种算法的误差影响不大，两种算法误差都比较小，但本发明实施例算法比R-K算法误差的均值要小。向前积分60min，本发明实施例算法优势更加明显，如以200为步长积分时，R-K算法距离误差均值为8.9717，而本发明实施例算法为8.6893。

（2）随着积分时间延长，两种算法的误差随之增大，但本发明实施例算法的均值误差小于R-K算法，且随着积分时间越长越明显，说明本发明实施例算法稳定性优于R-K算法。

（3）采用R-K方法进行一步积分，需要计算4次三维方向的加速度；采用本发明实施例方法进行一次运算，只需要计算2次；由于轨道积分运算的代价主要取决于计算加速度的次数，所以本发明实施例算法的运算量约为R-K算法的50%。

尽管为示例目的，已经公开了本发明的优选实施例，本领域的技术人员将意识到各种改进、增加和取代也是可能的，因此，本发明的范围应当不限于上述实施例。

Claims (2)

1.一种卫星位置解算方法，其特征在于，包括：

确定观测时刻t_e，从导航电文中获取参考时刻t_b的卫星位置参数，包括位置X₀＝{x₀,y₀,z₀}，速度V₀＝{v_x0，v_y0，v_z0}，计算t_b加速度a₀＝{a_x,a_y,a_z}，设置积分步长h，0<h<200，且为整数；确定需要积分的时间长度t_L＝t_e-t_b，并计算需要积分的次数

当N>4，利用线性单步方法计算前4个步长时刻的位置X₀、X₁、X₂、X₃和速度V₀、V₁、V₂、V₃；

根据X₀、X₁、X₂、X₃和速度V₀、V₁、V₂、V₃，计算剩余步长时刻的卫星位置、速度及加速度信息，得到观测时刻t_e的卫星位置，具体包括：

按照公式(3)计算t_i+1时刻卫星速度和位置的预测值之后，按照公式(4)对预测值进行校正，得到t_i+1时刻卫星位置X_i+1＝{x_i+1,y_i+1,z_i+1}和速度V_i+1＝{v_x,i+1,v_y,i+1,v_z,i+1}；

$\left\{\begin{array}{c}{a}_i=f(t_i,X_i,V_i,\stackrel{\&CenterDot;\&CenterDot;}{X})\\ {\stackrel{\&OverBar;}{V}}_{i+1}=V_i+\frac{h}{24}(55a_i-59a_{i-1}+37a_{i-2}-9a_{i-3})\\ {\stackrel{\&OverBar;}{X}}_{i+1}=X_i+\frac{h}{24}(55V_i-59V_{i-1}+37V_{i-2}-9V_{i-3})\end{array}\right.---\left(3\right)$

$\left\{\begin{array}{c}{a}_{i+1}=f(t_i+h,{\stackrel{\&OverBar;}{X}}_{i+1},{\stackrel{\&OverBar;}{V}}_{i+1},\stackrel{\&CenterDot;\&CenterDot;}{X})\\ V_{i+1}=V_i+\frac{h}{720}(251a_{i+1}+646a_i-264a_{i-1}+106a_{i-2}-19a_{i-3})\\ X_{i+1}=X_i+\frac{h}{720}[251V_{i+1}+646V_i-264V_{i-1}+106V_{i-2}-19V_{i-3}]\end{array}\right.---\left(4\right)$

判断是否i<N，如果是，则置X_i＝X_i+1,V_i＝V_i+1，积分计算下一个时刻的位置、速度；如果否得到t_N时刻的卫星位置、速度及加速度信息；

根据t_N时刻的卫星位置、速度及加速度信息，利用线性单步方法得到观测时刻t_e时刻的卫星位置：

其中，通过公式(2)计算t_b加速度a₀＝{a_x,a_y,a_z}；

$\left\{\begin{array}{c}{a}_x=-\frac{\mathrm{\&mu;}}{r^3}x-\frac{3}{2}{J}_0^2\&CenterDot;\frac{\mathrm{\&mu;}{a}_e^2}{r^5}\&CenterDot;x\&CenterDot;(1-\frac{5z^2}{r^2})+{\mathrm{\&omega;}}^{2}x+2\mathrm{\&omega;}\&CenterDot;v_y+\stackrel{\&CenterDot;\&CenterDot;}{x}\\ a_y=-\frac{\mathrm{\&mu;}}{r^3}y-\frac{3}{2}{J}_0^2\&CenterDot;\frac{\mathrm{\&mu;}{a}_e^2}{r^5}\&CenterDot;y\&CenterDot;(1-\frac{5z^2}{r^2})+{\mathrm{\&omega;}}^{2}y-2\mathrm{\&omega;}\&CenterDot;v_x+\stackrel{\&CenterDot;\&CenterDot;}{y}\\ a_z=-\frac{\mathrm{\&mu;}}{r^3}z-\frac{3}{2}{J}_0^2\&CenterDot;\frac{\mathrm{\&mu;}{a}_e^2}{r^5}\&CenterDot;z\&CenterDot;(3-\frac{5z^2}{r^2})+\stackrel{\&CenterDot;\&CenterDot;}{z}\end{array}\right.---\left(2\right)$

其中，地球引力常量μ＝398600.44km³/s²；卫星到地球质心的距离地球的长半轴a_e＝6378.136km；地球重力位第二带谐系数 地球自转速度ω＝0.00007292115rad/s；为日月摄动加速度在x、y、z三个方向的分量。

2.如权利要求1所述的卫星位置解算方法，其特征在于，当N≤4，利用线性单步方法，计算观测时刻的卫星位置。