CN103453906A - 卫星轨道的预测方法 - Google Patents

卫星轨道的预测方法 Download PDF

Info

Publication number
CN103453906A
CN103453906A CN2013103456134A CN201310345613A CN103453906A CN 103453906 A CN103453906 A CN 103453906A CN 2013103456134 A CN2013103456134 A CN 2013103456134A CN 201310345613 A CN201310345613 A CN 201310345613A CN 103453906 A CN103453906 A CN 103453906A
Authority
CN
China
Prior art keywords
orbit
satellite
hours
state
calculating
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN2013103456134A
Other languages
English (en)
Other versions
CN103453906B (zh
Inventor
陈曦
张中华
葛宁
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Tsinghua University
Original Assignee
Tsinghua University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Tsinghua University filed Critical Tsinghua University
Priority to CN201310345613.4A priority Critical patent/CN103453906B/zh
Publication of CN103453906A publication Critical patent/CN103453906A/zh
Application granted granted Critical
Publication of CN103453906B publication Critical patent/CN103453906B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Landscapes

  • Position Fixing By Use Of Radio Waves (AREA)

Abstract

本发明提供了一种卫星轨道的预测方法，包括以下步骤：获得过去N小时的卫星事后精密轨道，获得卫星初始状态；计算卫星受力参数；根据所述的卫星受力参数计算未来M小时卫星第一预测轨道根数；根据所述过去N小时的卫星事后精密轨道计算轨道根数修正系数；根据所述轨道根数修正系数修正第一预测轨道根数，获得修正轨道根数；利用所述修正轨道根数计算轨道状态，获得第一预测轨道状态；根据所述过去N小时的卫星事后精密轨道计算卫星初始状态修正系数；以及根据初始状态修正系数修正第一预测轨道状态，获得预测的卫星轨道状态。本发明的技术可用于预测卫星轨道，满足导航定位类卫星的基于预测的方法实现实时定轨需求。

Description

卫星轨道的预测方法
技术领域
本发明涉及一种卫星轨道预测方法，尤其是一种适用于低轨卫星的精密轨道预测方法，属于航天轨道动力学技术领域。
背景技术
近年来，低轨卫星导航正引起越来越多地关注。与现有的全球导航卫星相比，低轨卫星导航具有三大主要优势：相同发射功率下，较低的路径损耗给地面用户带来高于20dB的接收功率增益；卫星视角变化快，有利于快速精确的解算载波相位整周模糊度；适应导航能力泛在化发展趋势。
低轨卫星导航的一个基本问题，是实时定轨问题。在低轨卫星实时定轨领域，最常见的手段是在低轨卫星上安装全球导航卫星系统(Global Navigation Satellite System, GNSS)接收机。全球导航卫星系统主要包括中国的北斗(BeiDou-2)、美国的全球定位系统(Global Positioning System, GPS)、俄罗斯的全球导航卫星系统(Global Navigation Satellite System, GLONASS)和欧洲的伽利略系统(Galileo)。全球导航卫星系统接收机工作的基本原理是：接收到导航卫星发送的无线电信号并提取伪距，并根据4个以上伪距计算自身在地理坐标系中的位置，常见的解算算法有最小二乘法和卡尔曼滤波法。
然而，这种方法的缺点是：全球导航卫星系统自身的实时轨道确定存在误差，时钟也存在误差，这导致低轨卫星星载接收机的实时定位精度大于6米，而低轨导航增强通常要求低轨卫星的实时轨道确定精度达到甚至优于1米，从而导致卫星轨道的预测不够精确。
发明内容
综上所述，确有必要提供一种具有较高预测精度的卫星精密轨道的预测方法。
一种卫星轨道的预测方法，主要包括以下步骤：步骤S10，获得过去N小时的卫星事后精密轨道，获得卫星初始状态；步骤S20，根据所述过去N小时的卫星事后精密轨道计算卫星受力参数；步骤S30，根据所述的卫星受力参数计算未来M小时卫星第一预测轨道根数；步骤S40，根据所述过去N小时的卫星事后精密轨道计算轨道根数修正系数；步骤S50，根据所述轨道根数修正系数修正第一预测轨道根数，获得修正轨道根数；步骤S60，利用所述修正轨道根数计算轨道状态，获得第一预测轨道状态；步骤S70，根据所述过去N小时的卫星事后精密轨道计算卫星初始状态修正系数；以及步骤S80，根据初始状态修正系数修正第一预测轨道状态，获得预测的卫星轨道状态。
本发明提供的卫星精密轨道的预测方法，利用事后精密轨道数据来预测未来轨道，综合利用了轨道动力学、事后精密轨道数据和低轨卫星轨道根数特征，其预测位置精度能达到米级。
附图说明
图1是本发明第一实施例提供的低轨卫星导航系统示意图。
图2是本发明提供的卫星轨道预测方法流程图。
图3本发明提供的卫星轨道预测方法应用于GRACE卫星的轨道预测位置误差。
图4本发明提供的卫星轨道预测方法应用于GRACE卫星的轨道预测速度误差。
图5是本发明第二实施例提供的卫星轨道预测方法。
具体实施方式
下面根据说明书附图并结合具体实施例对本发明的技术方案进一步详细表述。
请参阅图1，S1~S5为传统导航卫星，它们为中轨导航卫星，工作轨道高度一般为2万公里。现有一个低轨导航星L1，工作轨道不高于3000公里，能够发射导航信号。一用户终端U1通过接收低轨导航星L1和传统导航卫星S1~S5的信号实现定位。在这个系统中，为了提高导航精度，要求低轨导航星L1的精度越高越好。
请一并参阅图2，图2为本发明提供的卫星轨道预测方法流程图，主要包括以下步骤：
步骤S10，获得过去N小时的卫星事后精密轨道，获得卫星初始状态；
步骤S20，根据所述过去N小时的卫星事后精密轨道计算卫星受力参数；
步骤S30，根据所述的卫星受力参数计算未来M小时卫星第一预测轨道根数；
步骤S40，根据所述过去N小时的卫星事后精密轨道计算轨道根数修正系数；
步骤S50，根据所述轨道根数修正系数修正第一预测轨道根数，获得修正轨道根数；
步骤S60，利用所述修正轨道根数计算轨道状态，获得第一预测轨道状态；
步骤S70，根据所述过去N小时的卫星事后精密轨道计算卫星初始状态修正系数；以及
步骤S80，根据初始状态修正系数修正第一预测轨道状态，获得预测的卫星轨道状态。
在步骤S10中，所述“事后精密轨道”，指的是过去N小时的精密卫星轨道。所述卫星事后精密轨道可通过常规技术获得，典型精度为1~10厘米。具体的，在卫星过顶时，将过去N小时星上的全球导航卫星系统接收机的初始观测量下载到地面，使用地面根据全球导航卫星系统的轨道、钟差和电离层修正数据产品经过约化动力学等方法进行计算，得到事后精密轨道。
在步骤S20中，根据N小时的卫星事后精密轨道计算卫星受力参数的方法，可进一步包括以下步骤：
步骤S21，将卫星所受力假想为保守力，并使用重力模型表示该保守力；
步骤S22，根据N小时的卫星事后精密轨道计算重力模型参数修正量；
步骤S23，利用所述重力模型参数修正量修正重力模型，所获得的重力模型即所求卫星受力参数。
在步骤S22中，所述重力模型参数修正量包括余弦规范化球面谐波系数                                               
Figure 2013103456134100002DEST_PATH_IMAGE001
的修正量以及正弦规范化球面谐波系数
Figure 2013103456134100002DEST_PATH_IMAGE003
的修正量
Figure 2013103456134100002DEST_PATH_IMAGE004
。其中，
Figure 261056DEST_PATH_IMAGE001
和
Figure 920944DEST_PATH_IMAGE003
来源于重力场模型，在地心大地坐标系中的具体表达式为：
；
其中，
Figure 2013103456134100002DEST_PATH_IMAGE006
为地球摄动势能；
Figure 2013103456134100002DEST_PATH_IMAGE007
为地球引力常数；
Figure 2013103456134100002DEST_PATH_IMAGE008
为地球平均赤道半径；r为卫星事后精密轨道径矢；Pnm(
Figure 2013103456134100002DEST_PATH_IMAGE009
)为的n阶m次伴随勒让德多项式。在一般应用中，可使用重力模型获得规范化球面谐波系数的值、
Figure 896487DEST_PATH_IMAGE003
，如采用GRACE重力模型(GRACE Gravity Model 03, GGM03)参数。所述
Figure 949894DEST_PATH_IMAGE001
和
Figure 789674DEST_PATH_IMAGE003
的初始化方法可采用GGM03重力模型中的
Figure 270334DEST_PATH_IMAGE001
和作为卫星受力的初始参数，但并不限于所述方法。所述重力模型参数修正量为实际重力场规范化球面谐波系数与已知的规范化球面谐波系数的差值，包括前述余弦规范化球面谐波系数修正量
Figure 482189DEST_PATH_IMAGE002
和正弦规范化球面谐波系数修正量
Figure 442055DEST_PATH_IMAGE004
。
所述重力模型参数修正量
Figure 93616DEST_PATH_IMAGE002
和
Figure 184938DEST_PATH_IMAGE004
的计算可包括以下步骤：
步骤S221，将过去N小时中根据事后精密定轨结果所获得的任意时刻轨道状态定义为过去N小时中任意时刻的真实轨道状态；
步骤S222，根据所述真实轨道状态进行重力场反演，得到重力模型参数修正量
Figure 580147DEST_PATH_IMAGE002
和
Figure 394519DEST_PATH_IMAGE004
。
在步骤S222中，重力场反演的具体方法可包括如下步骤：
将过去N小时的真实轨道状态在动力学积分轨道处进行下述泰勒展开：
Figure 2013103456134100002DEST_PATH_IMAGE010
其中，
Figure 2013103456134100002DEST_PATH_IMAGE011
给定初始状态及初始重力场模型里预测的轨道状态；
Figure 2013103456134100002DEST_PATH_IMAGE012
初始轨道状态误差；
为状态向量对初始状态的偏导数，即状态转移矩阵；
Figure 2013103456134100002DEST_PATH_IMAGE014
为状态对重力场球面谐波系数的偏导数，即参数敏感矩阵。
其中，h即为待求重力场球面谐波系数构成的列向量，元素个数为nh个，
Figure 2013103456134100002DEST_PATH_IMAGE015
的排列顺序不限，举例如下：
Figure 2013103456134100002DEST_PATH_IMAGE016
。
h的初值
Figure 2013103456134100002DEST_PATH_IMAGE017
可以采用GRACE重力模型(GRACE Gravity Model 03, GGM03)参数。
所述状态转移矩阵、参数敏感矩阵满足如下关系
Figure 2013103456134100002DEST_PATH_IMAGE018
其中，
Figure 2013103456134100002DEST_PATH_IMAGE019
， 
Figure 2013103456134100002DEST_PATH_IMAGE020
。
若令
Figure 2013103456134100002DEST_PATH_IMAGE021
，则状态转移矩阵和参数敏感矩阵的联合微分可写为：
Figure 2013103456134100002DEST_PATH_IMAGE022
(1)
考虑到
Figure 2013103456134100002DEST_PATH_IMAGE023
中时间参数
Figure 2013103456134100002DEST_PATH_IMAGE024
连续且与
Figure 2013103456134100002DEST_PATH_IMAGE025
和
Figure 2013103456134100002DEST_PATH_IMAGE026
相关性小或者无关，可以交换微分的先后次序：
Figure 2013103456134100002DEST_PATH_IMAGE027
令
Figure 2013103456134100002DEST_PATH_IMAGE028
，
Figure 2013103456134100002DEST_PATH_IMAGE029
，
Figure 2013103456134100002DEST_PATH_IMAGE030
，上述变分方程可写为
(2)
其中
Figure 2013103456134100002DEST_PATH_IMAGE032
（N为重力场球谐系数的阶数）为球谐系数个数，由于卫星加速度由受力决定，可令。
将式
Figure 406862DEST_PATH_IMAGE010
写成矩阵形式为
Figure 2013103456134100002DEST_PATH_IMAGE034
计算出和
Figure 2013103456134100002DEST_PATH_IMAGE036
后，即根据所述反演变分方程形式(2)并利用最小二乘法解算地球重力场球面谐波系数修正量
Figure 2013103456134100002DEST_PATH_IMAGE037
。
所述反演变分方程形式（2）可通过微分法求解。
在步骤S23中，获得了
Figure 2013103456134100002DEST_PATH_IMAGE038
后，就可以获得反演后的
Figure 2013103456134100002DEST_PATH_IMAGE039
,即修正后的重力模型参数，即所述的卫星受力参数。
在步骤S30中，选择一个初始时间和历史精密轨道得到初始卫星状态，根据卫星初始状态和卫星受力参数来计算未来M小时卫星轨道根数（即第一预测轨道根数）。
所述第一预测轨道根数包括轨道倾角
Figure 2013103456134100002DEST_PATH_IMAGE040
和升交点赤经
Figure 2013103456134100002DEST_PATH_IMAGE041
等的计算方法可通过对下式进行积分获得：
Figure 2013103456134100002DEST_PATH_IMAGE043
Figure 2013103456134100002DEST_PATH_IMAGE044
Figure 2013103456134100002DEST_PATH_IMAGE045
Figure 2013103456134100002DEST_PATH_IMAGE046
其中，
Figure 2013103456134100002DEST_PATH_IMAGE048
摄动加速度分量，方向分别为径向背向地心为正，横向
Figure 2013103456134100002DEST_PATH_IMAGE050
速度增加方向为正，轨道面法线方向
Figure 2013103456134100002DEST_PATH_IMAGE051
；
Figure 2013103456134100002DEST_PATH_IMAGE052
真近点角；
Figure 2013103456134100002DEST_PATH_IMAGE053
平均角速度；
Figure 2013103456134100002DEST_PATH_IMAGE054
；
Figure 2013103456134100002DEST_PATH_IMAGE055
，升交点角距；
Figure 2013103456134100002DEST_PATH_IMAGE056
，其中
Figure 2013103456134100002DEST_PATH_IMAGE057
为卫星轨道半长轴，
Figure 2013103456134100002DEST_PATH_IMAGE058
为卫星离心率；
Figure 2013103456134100002DEST_PATH_IMAGE059
；
摄动源为非球形引力摄动，是为保守力，则摄动加速度为
Figure DEST_PATH_IMAGE060
，
Figure DEST_PATH_IMAGE061
，
Figure DEST_PATH_IMAGE062
至此，我们获得了未来M小时卫星第一预测轨道根数，包括第一预测轨道根数中的轨道倾角
Figure 169501DEST_PATH_IMAGE040
和第一预测轨道根数升交点赤经等。
在步骤S40中，将卫星所有受力归结为重力进行轨道预测仍然存在误差，主要是轨道倾角误差和升交点赤经误差。因此，我们需要进一步根据所述过去N小时的卫星事后精密轨道计算轨道根数修正量，并根据所述轨道根数修正系数修正轨道根数，获得修正后的修正轨道根数。具体地，轨道根数修正量包括轨道倾角修正量和升交点赤经修正量。
所述轨道倾角修正量可通过以下方法计算：
步骤S411，计算过去N小时内事后精密轨道状态对应的轨道根数的轨道倾角
Figure DEST_PATH_IMAGE063
；
步骤S412，将过去N小时的与过去N小时第一预测轨道根数的轨道倾角相减，获得
Figure DEST_PATH_IMAGE064
；
步骤S413，用多项式拟合
Figure 342228DEST_PATH_IMAGE064
的长期缓慢变化：，其中，
Figure DEST_PATH_IMAGE066
的具体值可以根据精度要求给定，是一个设计参数，
Figure DEST_PATH_IMAGE067
是过去N小时的起始时刻；
步骤S414，根据将过去N小时的的值使用最小二乘方法计算上式中的
Figure DEST_PATH_IMAGE068
；以及
步骤S415，根据
Figure 161465DEST_PATH_IMAGE068
计算未来M小时的轨道倾角修正量。
所述升交点赤经修正量可通过以下方法计算：
步骤S421，计算过去N小时内事后精密轨道状态对应的轨道根数的升交点赤经
Figure DEST_PATH_IMAGE069
；
步骤S422，将过去N小时的
Figure 497900DEST_PATH_IMAGE069
与过去N小时第一预测轨道根数的升交点赤经相减，获得
Figure DEST_PATH_IMAGE070
；
步骤S423，用多项式拟合
Figure 724482DEST_PATH_IMAGE069
的长期缓慢变化：
Figure DEST_PATH_IMAGE071
，其中，
Figure DEST_PATH_IMAGE072
的具体值可以根据精度要求给定，是一个设计参数，
Figure 966107DEST_PATH_IMAGE067
是过去N小时的起始时刻；
步骤S424，根据过去N小时的
Figure 821324DEST_PATH_IMAGE070
的值使用最小二乘方法计算上式中的
Figure DEST_PATH_IMAGE073
；
步骤S425，根据计算未来M小时的升交点赤经修正量。
在步骤S50中，所述修正后的修正轨道根数包括修正后的轨道倾角和修正后的升交点赤经，其中所述修正后的轨道倾角为
Figure DEST_PATH_IMAGE074
，修正后的轨道根数中的升交点赤经为
Figure DEST_PATH_IMAGE075
。
在步骤S60中，所述轨道状态包括卫星在地心坐标系下的三维位置和速度，所述第一预测轨道状态就可以通过上述修正后的轨道倾角和修正后的升交点赤经计算。
在步骤S70中，初始状态误差也是影响轨道预测的重要误差项。因此，我们需要进一步根据所述过去N小时的卫星事后精密轨道计算初始状态修正量，并根据初始状态修正量修正第一预测轨道状态，获得修正轨道状态即最终轨道状态。
所述初始误差修正量的计算方法如下：
根据位置与加速度的关系，由卫星初始状态造成的位置误差可以表示为
Figure DEST_PATH_IMAGE076
，其中，x为任意一轴的卫星实际轨道位置分量。、
Figure DEST_PATH_IMAGE078
由卫星初始状态误差决定，这说明在初始状态误差一定的情况下，某时刻有初始状态引起的位置预测误差与当前的位置为线性关系。若知道K（K大于等于2）个时刻对应的位置分量x对应预测位置误差
Figure DEST_PATH_IMAGE079
，则上式可写为
Figure DEST_PATH_IMAGE080
即可通过最小二乘法计算出
Figure 737645DEST_PATH_IMAGE077
和
Figure 512572DEST_PATH_IMAGE078
，进而可以由上式进而根据A、B计算未来M小时由初始状态引起的位置误差，即所述初始误差修正量。
步骤S80，将初始误差修正量加到第一预测轨道状态，就获得了修正轨道状态即最终轨道状态。
根据本发明的方法，可以利用前一天精密轨道数据，预测下一天的低轨卫星轨道。作为具体的应用，我们针对GRACE卫星的数据进行理论计算，I=J=6，得到的结果如图3和图4所示。其中图3的横坐标单位为分钟，纵坐标单位为米。可以看到，实际单轴24小时的预测误差最大约为1米。图4是速度预测误差，其中横坐标单位为分钟，纵坐标为单位为米/秒，可以看到，单轴速度预测速度误差最大约为0.06米/秒。
进一步，请参阅图5，图5为本发明第二实施例提供的卫星轨道预测方法的示意图，图5的情况出现在本发明的方法首次运行时。不失一般性，我们假设N=M。在时刻1，我们可以使用第一个N小时的事后精密轨道和重力模型做轨道状态预测；由于不存在第一个N小时的预测轨道，因此不能做根据本发明的修正计算，包括轨道根数修正和初始轨道状态。
根据本发明的轨道预测方法进行实时定轨的方法，包括以下步骤：
（1）在卫星过顶时下传过去N小时的全球导航卫星接收机观测量；
（2）根据本发明的轨道预测方法预测未来M小时的轨道状态,并将所预测的未来M小时的轨道状态在卫星过顶同时上传的卫星；
（3）在未来M小时使用根据当前时间所上传的预测的轨道状态，通过内插计算卫星的实时轨道位置和速度。
在实际应用中，卫星的平均过顶间隔约为12小时，考虑事后精密定轨的延迟，也可以考虑将时间分为若干12小时，即M=N=12小时，使用根据过去[2N,N]小时的事后精密轨道和本发明的轨道预测方法预测未来M小时的轨道状态（也可以是轨道根数）也能得到较好的预测结果，计算过程可以发生在两次卫星过顶之间。因此，本发明的上述实施例，是为了说明而不是限制本发明的方法。
本发明提供的卫星精密轨道的预测方法，利用事后精密轨道数据来预测未来轨道，综合利用了轨道动力学、事后精密轨道数据和低轨卫星轨道根数特征，其预测位置精度能达到米级，实现了24小时预测精度优于1米。本发明的技术可用于预测卫星轨道，满足导航定位类卫星的基于预测的方法实现实时定轨需求，从而能够解决低轨卫星的实时轨道计算问题。另外，所述卫星精密轨道的预测方法不仅可以用于低轨导航星，也可以用于其它需要实时精密确定轨道的卫星。
另外，本领域技术人员还可在本发明精神内作其它变化，当然这些依据本发明精神所作的变化，都应包含在本发明所要求保护的范围内。

Claims (13)

1.一种卫星轨道的预测方法，主要包括以下步骤：
步骤S10，获得过去N小时的卫星事后精密轨道，获得卫星初始状态；
步骤S20，根据所述过去N小时的卫星事后精密轨道计算卫星受力参数；
步骤S30，根据所述的卫星受力参数计算未来M小时卫星第一预测轨道根数；
步骤S40，根据所述过去N小时的卫星事后精密轨道计算轨道根数修正系数；
步骤S50，根据所述轨道根数修正系数修正第一预测轨道根数，获得修正轨道根数；
步骤S60，利用所述修正轨道根数计算轨道状态，获得第一预测轨道状态；
步骤S70，根据所述过去N小时的卫星事后精密轨道计算卫星初始状态修正系数；以及
步骤S80，根据初始状态修正系数修正第一预测轨道状态，获得预测的卫星轨道状态。
2.如权利要求1所述的卫星轨道的预测方法，其特征在于，所述卫星受力参数的计算方法包括以下步骤：
将卫星所受力假想为保守力，并使用重力模型表示该保守力；
根据N小时的卫星事后精密轨道计算重力模型参数修正量；
根据所述重力模型参数修正量修正重力模型，修正后的重力模型即所述卫星受力参数。
3.如权利要求2所述的卫星轨道的预测方法，其特征在于，所述重力模型参数修正量包括余弦规范化球面谐波系数                                               
Figure 2013103456134100001DEST_PATH_IMAGE001
的修正量
Figure 2013103456134100001DEST_PATH_IMAGE002
以及正弦规范化球面谐波系数
Figure 2013103456134100001DEST_PATH_IMAGE003
的修正量
Figure 2013103456134100001DEST_PATH_IMAGE004
。
4.如权利要求3为所述的卫星轨道的预测方法，其特征在于，
Figure 845087DEST_PATH_IMAGE001
和来源于重力场模型，在地心大地坐标系中的具体表达式为：
；
其中，
Figure 2013103456134100001DEST_PATH_IMAGE006
为地球摄动势能；
Figure 2013103456134100001DEST_PATH_IMAGE007
为地球引力常数；
Figure 2013103456134100001DEST_PATH_IMAGE008
为地球平均赤道半径；r为卫星事后精密轨道径矢；Pnm(
Figure 2013103456134100001DEST_PATH_IMAGE009
)为
Figure 643465DEST_PATH_IMAGE009
的n阶m次伴随勒让德多项式。
5.如权利要求4所述的卫星轨道的预测方法，其特征在于，所述重力模型参数修正量
Figure 852729DEST_PATH_IMAGE002
和
Figure 880728DEST_PATH_IMAGE004
的计算包括以下步骤：
将过去N小时中根据事后精密定轨结果所获得的任意时刻轨道状态定义为过去N小时中任意时刻的真实轨道状态；
根据所述真实轨道状态进行重力场反演，得到重力模型参数修正量
Figure 62311DEST_PATH_IMAGE002
和
Figure 330612DEST_PATH_IMAGE004
。
6.如权利要求5所述的卫星轨道的预测方法，其特征在于，所述重力场反演包括如下步骤：
将过去N小时的真实轨道状态在动力学积分轨道处进行下述泰勒展开：
Figure 2013103456134100001DEST_PATH_IMAGE010
，
其中，
Figure 2013103456134100001DEST_PATH_IMAGE011
给定初始状态及初始重力场模型里预测的轨道状态；
Figure 2013103456134100001DEST_PATH_IMAGE012
初始轨道状态误差；
Figure 2013103456134100001DEST_PATH_IMAGE013
为状态向量对初始状态的偏导数，即状态转移矩阵；
Figure 2013103456134100001DEST_PATH_IMAGE014
为状态对重力场球面谐波系数的偏导数，即参数敏感矩阵；
其中，h即为重力场球面谐波系数
Figure 2013103456134100001DEST_PATH_IMAGE015
构成的列向量，元素个数为nh个，h的初值
Figure 2013103456134100001DEST_PATH_IMAGE016
采用GRACE重力模型参数；
将式
Figure 669583DEST_PATH_IMAGE010
写成矩阵形式：
Figure 2013103456134100001DEST_PATH_IMAGE017
；
计算
Figure 2013103456134100001DEST_PATH_IMAGE018
和
Figure DEST_PATH_IMAGE019
，根据反演变分方程式：
Figure 2013103456134100001DEST_PATH_IMAGE020
；
并利用最小二乘法解算地球重力场球面谐波系数修正量
Figure 111060DEST_PATH_IMAGE002
和
Figure 147149DEST_PATH_IMAGE004
，其中
Figure DEST_PATH_IMAGE021
（N为重力场球谐系数的阶数）为球谐系数个数。
7.如权利要求1所述的卫星轨道的预测方法，其特征在于，所述第一预测轨道根数包括轨道倾角
Figure DEST_PATH_IMAGE022
和升交点赤经
Figure DEST_PATH_IMAGE023
，所述轨道倾角
Figure 147204DEST_PATH_IMAGE022
和升交点赤经
Figure 268743DEST_PATH_IMAGE023
的计算方法通过对下式进行积分获得：
Figure DEST_PATH_IMAGE025
Figure DEST_PATH_IMAGE026
Figure DEST_PATH_IMAGE027
Figure DEST_PATH_IMAGE028
Figure DEST_PATH_IMAGE029
；
其中，
Figure DEST_PATH_IMAGE030
摄动加速度分量，方向分别为径向
Figure DEST_PATH_IMAGE031
背向地心为正，横向
Figure DEST_PATH_IMAGE032
速度增加方向为正，轨道面法线方向
Figure DEST_PATH_IMAGE033
；
Figure DEST_PATH_IMAGE034
真近点角；
Figure DEST_PATH_IMAGE035
平均角速度；
Figure DEST_PATH_IMAGE036
；
Figure DEST_PATH_IMAGE037
，升交点角距；
Figure DEST_PATH_IMAGE038
，其中
Figure DEST_PATH_IMAGE039
为卫星轨道半长轴，
Figure DEST_PATH_IMAGE040
为卫星离心率；
Figure DEST_PATH_IMAGE041
。
8.如权利要求7所述的卫星轨道的预测方法，其特征在于，摄动源为非球形引力摄动，是为保守力，摄动加速度为：
Figure DEST_PATH_IMAGE042
，
Figure DEST_PATH_IMAGE043
，。
9.如权利要求8所述的卫星轨道的预测方法，其特征在于，所述轨道根数修正量包括轨道倾角修正量和升交点赤经修正量。
10.如权利要求9所述的卫星轨道的预测方法，其特征在于，所述轨道倾角修正量的计算包括以下步骤：
计算过去N小时内事后精密轨道状态对应的轨道根数的轨道倾角
Figure DEST_PATH_IMAGE045
；
将过去N小时的
Figure 9516DEST_PATH_IMAGE045
与过去N小时第一预测轨道根数的轨道倾角相减，获得；
用多项式拟合
Figure 480205DEST_PATH_IMAGE046
的长期缓慢变化：
Figure DEST_PATH_IMAGE047
，其中，
Figure DEST_PATH_IMAGE048
的具体值根据精度要求给定，是一个设计参数，是过去N小时的起始时刻；
根据将过去N小时的
Figure 339577DEST_PATH_IMAGE046
的值使用最小二乘方法计算上式中的
Figure DEST_PATH_IMAGE050
；以及
根据
Figure 761462DEST_PATH_IMAGE050
计算未来M小时的轨道倾角修正量。
11.如权利要求10所述的卫星轨道的预测方法，其特征在于，所述升交点赤经修正量的计算包括以下步骤：
计算过去N小时内事后精密轨道状态对应的轨道根数的升交点赤经
Figure DEST_PATH_IMAGE051
；
将过去N小时的
Figure 731692DEST_PATH_IMAGE051
与过去N小时第一预测轨道根数的升交点赤经相减，获得
Figure DEST_PATH_IMAGE052
；
用多项式拟合
Figure 742374DEST_PATH_IMAGE051
的长期缓慢变化：
Figure DEST_PATH_IMAGE053
，其中，
Figure DEST_PATH_IMAGE054
的具体值根据精度要求给定，是一个设计参数，
Figure 84231DEST_PATH_IMAGE049
是过去N小时的起始时刻；
根据过去N小时的
Figure 180363DEST_PATH_IMAGE052
的值使用最小二乘方法计算上式中的；
步骤S425，根据
Figure 705016DEST_PATH_IMAGE055
计算未来M小时的升交点赤经修正量。
12.如权利要求11所述的卫星轨道的预测方法，其特征在于，所述初始误差修正量由以下方法计算：由卫星初始状态造成的位置误差表示为
Figure DEST_PATH_IMAGE056
，其中x为任意一轴的卫星实际轨道位置分量;
Figure DEST_PATH_IMAGE057
、由卫星初始状态误差决定；通过最小二乘法计算出
Figure 898100DEST_PATH_IMAGE057
和
Figure 181550DEST_PATH_IMAGE058
，进而上式计算未来M小时由初始状态引起的位置误差
Figure DEST_PATH_IMAGE059
，即所述初始误差修正量。
13.如权利要求1所述的卫星轨道的预测方法，其特征在于，进一步包括以下步骤：
在卫星过顶时下传过去N小时的全球导航卫星接收机观测量；
预测卫星未来M小时的轨道状态,并将所预测的未来M小时的轨道状态在卫星过顶同时上传；
在未来M小时使用根据当前时间所上传的预测的轨道状态，通过内插计算卫星的实时轨道位置和速度。
CN201310345613.4A 2013-08-09 2013-08-09 卫星轨道的预测方法 Active CN103453906B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201310345613.4A CN103453906B (zh) 2013-08-09 2013-08-09 卫星轨道的预测方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201310345613.4A CN103453906B (zh) 2013-08-09 2013-08-09 卫星轨道的预测方法

Publications (2)

Publication Number Publication Date
CN103453906A true CN103453906A (zh) 2013-12-18
CN103453906B CN103453906B (zh) 2016-04-27

Family

ID=49736534

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201310345613.4A Active CN103453906B (zh) 2013-08-09 2013-08-09 卫星轨道的预测方法

Country Status (1)

Country Link
CN (1) CN103453906B (zh)

Cited By (11)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104615579A (zh) * 2014-12-30 2015-05-13 中国科学院数学与系统科学研究院 基于最大模分解的卫星轨道确定方法及装置
CN105044745A (zh) * 2015-07-15 2015-11-11 中国人民解放军理工大学 一种圆轨道低轨卫星过顶剩余可见时长预测方法
CN105573118A (zh) * 2015-12-16 2016-05-11 中国人民解放军国防科学技术大学 快速重访卫星轨道设计方法
CN105737834A (zh) * 2014-12-09 2016-07-06 上海新跃仪表厂 一种基于轨道平根数的相对导航鲁棒滤波方法
CN109000666A (zh) * 2018-06-05 2018-12-14 北京电子工程总体研究所 一种基于中心天体矢量观测的自主定轨方法及其系统
CN110068846A (zh) * 2019-04-30 2019-07-30 上海微小卫星工程中心 一种基于星载gnss接收机在星上自主确定轨道平根数的方法
CN110068845A (zh) * 2019-04-30 2019-07-30 上海微小卫星工程中心 一种基于平根数理论确定卫星理论轨道的方法
CN112013851A (zh) * 2020-06-30 2020-12-01 南京天际易达通信技术有限公司 一种卫星运控轨道计算方法
CN114118504A (zh) * 2020-08-27 2022-03-01 哈尔滨工业大学 一种卫星轨道预测方法及系统
CN115096319A (zh) * 2022-08-24 2022-09-23 航天宏图信息技术股份有限公司 一种基于光学测角数据的星链卫星初轨确定方法和装置
CN116886178A (zh) * 2023-09-06 2023-10-13 北京融为科技有限公司 轨道预报修正方法及装置

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20060026143A1 (en) * 2004-08-02 2006-02-02 Hirth Gerhard A System for querying databases
CN1923622A (zh) * 2006-07-07 2007-03-07 中国科学院力学研究所 一种卫星飞行参数实时预测方法
WO2009065206A1 (en) * 2007-11-19 2009-05-28 Rx Networks Inc. Distributed orbit modeling and propagation method for a predicted and real-time assisted gps system
CN101446628A (zh) * 2007-11-26 2009-06-03 联发科技股份有限公司 用于预测gnss卫星轨迹延伸数据的方法及装置
CN102636795A (zh) * 2012-04-27 2012-08-15 清华大学 一种多接收机网络化无线定位方法

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20060026143A1 (en) * 2004-08-02 2006-02-02 Hirth Gerhard A System for querying databases
CN1923622A (zh) * 2006-07-07 2007-03-07 中国科学院力学研究所 一种卫星飞行参数实时预测方法
WO2009065206A1 (en) * 2007-11-19 2009-05-28 Rx Networks Inc. Distributed orbit modeling and propagation method for a predicted and real-time assisted gps system
CN101446628A (zh) * 2007-11-26 2009-06-03 联发科技股份有限公司 用于预测gnss卫星轨迹延伸数据的方法及装置
CN102636795A (zh) * 2012-04-27 2012-08-15 清华大学 一种多接收机网络化无线定位方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
匡翠林等: "低轨卫星与GPS导航卫星联合定规研究", 《大地测量与地球动力学》 *
张雄等: "利用GPS伪距观测值确定低轨卫星轨道的探讨", 《地理空间信息》 *

Cited By (19)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105737834A (zh) * 2014-12-09 2016-07-06 上海新跃仪表厂 一种基于轨道平根数的相对导航鲁棒滤波方法
CN105737834B (zh) * 2014-12-09 2018-06-26 上海新跃仪表厂 一种基于轨道平根数的相对导航鲁棒滤波方法
CN104615579A (zh) * 2014-12-30 2015-05-13 中国科学院数学与系统科学研究院 基于最大模分解的卫星轨道确定方法及装置
CN105044745A (zh) * 2015-07-15 2015-11-11 中国人民解放军理工大学 一种圆轨道低轨卫星过顶剩余可见时长预测方法
CN105573118A (zh) * 2015-12-16 2016-05-11 中国人民解放军国防科学技术大学 快速重访卫星轨道设计方法
CN105573118B (zh) * 2015-12-16 2018-11-20 中国人民解放军国防科学技术大学 快速重访卫星轨道设计方法
CN109000666A (zh) * 2018-06-05 2018-12-14 北京电子工程总体研究所 一种基于中心天体矢量观测的自主定轨方法及其系统
CN109000666B (zh) * 2018-06-05 2021-11-23 北京电子工程总体研究所 一种基于中心天体矢量观测的自主定轨方法及其系统
CN110068845A (zh) * 2019-04-30 2019-07-30 上海微小卫星工程中心 一种基于平根数理论确定卫星理论轨道的方法
CN110068845B (zh) * 2019-04-30 2021-07-23 上海微小卫星工程中心 一种基于平根数理论确定卫星理论轨道的方法
CN110068846A (zh) * 2019-04-30 2019-07-30 上海微小卫星工程中心 一种基于星载gnss接收机在星上自主确定轨道平根数的方法
CN110068846B (zh) * 2019-04-30 2022-01-07 上海微小卫星工程中心 一种基于星载gnss接收机在星上自主确定轨道平根数的方法
CN112013851A (zh) * 2020-06-30 2020-12-01 南京天际易达通信技术有限公司 一种卫星运控轨道计算方法
CN114118504A (zh) * 2020-08-27 2022-03-01 哈尔滨工业大学 一种卫星轨道预测方法及系统
CN114118504B (zh) * 2020-08-27 2024-07-23 哈尔滨工业大学 一种卫星轨道预测方法及系统
CN115096319A (zh) * 2022-08-24 2022-09-23 航天宏图信息技术股份有限公司 一种基于光学测角数据的星链卫星初轨确定方法和装置
CN115096319B (zh) * 2022-08-24 2022-11-18 航天宏图信息技术股份有限公司 一种基于光学测角数据的星链卫星初轨确定方法和装置
CN116886178A (zh) * 2023-09-06 2023-10-13 北京融为科技有限公司 轨道预报修正方法及装置
CN116886178B (zh) * 2023-09-06 2024-01-19 北京融为科技有限公司 轨道预报修正方法及装置

Also Published As

Publication number Publication date
CN103453906B (zh) 2016-04-27

Similar Documents

Publication Publication Date Title
CN103453906B (zh) 卫星轨道的预测方法
Han et al. Integrated GPS/INS navigation system with dual-rate Kalman Filter
CN106405670B (zh) 一种适用于捷联式海洋重力仪的重力异常数据处理方法
CN106767787A (zh) 一种紧耦合gnss/ins组合导航装置
CN108344415B (zh) 一种组合导航信息融合方法
CN105371844B (zh) 一种基于惯性/天文互助的惯性导航系统初始化方法
Georgy et al. Vehicle navigator using a mixture particle filter for inertial sensors/odometer/map data/GPS integration
CN107015259B (zh) 采用多普勒测速仪计算伪距/伪距率的紧组合方法
CN103487820B (zh) 一种车载捷联/卫星紧组合无缝导航方法
CN103542854A (zh) 基于星载处理器的自主定轨方法
CN102508954A (zh) Gps/sins组合导航全数字仿真方法及装置
Bock et al. GPS single-frequency orbit determination for low Earth orbiting satellites
CN106997061B (zh) 一种基于扰动星间相对速度提高重力场反演精度的方法
CN103868514A (zh) 一种在轨飞行器自主导航系统
CN104048664A (zh) 一种导航卫星星座自主定轨的方法
CN103017787A (zh) 适用于摇摆晃动基座的初始对准方法
CN110727003A (zh) 一种北斗卫星导航系统的伪距仿真方法
Iqbal et al. Pseudoranges error correction in partial GPS outages for a nonlinear tightly coupled integrated system
CN101893712A (zh) 用于地球静止卫星精密定轨的选权拟合方法
Du et al. Experimental study on GPS non-linear least squares positioning algorithm
CN106643726B (zh) 一种统一惯性导航解算方法
CN102914781A (zh) 一种格洛纳斯卫星信号星历电文生成方法及装置
Mahmoud et al. Integrated INS/GPS navigation system
CN110554443A (zh) 基于载波相位观测值和点加速度法确定地球重力场的方法
Zhou Low-cost MEMS-INS/GPS integration using nonlinear filtering approaches

Legal Events

Date Code Title Description
C06 Publication
PB01 Publication
C10 Entry into substantive examination
SE01 Entry into force of request for substantive examination
C14 Grant of patent or utility model
GR01 Patent grant