CN115327587A - 基于gnss定位信息的低轨卫星轨道误差修正方法及系统 - Google Patents

基于gnss定位信息的低轨卫星轨道误差修正方法及系统 Download PDF

Info

Publication number
CN115327587A
CN115327587A CN202210905416.2A CN202210905416A CN115327587A CN 115327587 A CN115327587 A CN 115327587A CN 202210905416 A CN202210905416 A CN 202210905416A CN 115327587 A CN115327587 A CN 115327587A
Authority
CN
China
Prior art keywords
low
orbit satellite
orbit
satellite
equation
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.)
Pending
Application number
CN202210905416.2A
Other languages
English (en)
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.)
Beijing Automation Control Equipment Institute BACEI
Original Assignee
Beijing Automation Control Equipment Institute BACEI
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 Beijing Automation Control Equipment Institute BACEI filed Critical Beijing Automation Control Equipment Institute BACEI
Priority to CN202210905416.2A priority Critical patent/CN115327587A/zh
Publication of CN115327587A publication Critical patent/CN115327587A/zh
Pending legal-status Critical Current

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S19/00—Satellite radio beacon positioning systems; Determining position, velocity or attitude using signals transmitted by such systems
    • G01S19/38—Determining a navigation solution using signals transmitted by a satellite radio beacon positioning system
    • G01S19/39—Determining a navigation solution using signals transmitted by a satellite radio beacon positioning system the satellite radio beacon positioning system transmitting time-stamped messages, e.g. GPS [Global Positioning System], GLONASS [Global Orbiting Navigation Satellite System] or GALILEO
    • G01S19/40—Correcting position, velocity or attitude
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S19/00—Satellite radio beacon positioning systems; Determining position, velocity or attitude using signals transmitted by such systems
    • G01S19/01—Satellite radio beacon positioning systems transmitting time-stamped messages, e.g. GPS [Global Positioning System], GLONASS [Global Orbiting Navigation Satellite System] or GALILEO
    • G01S19/13—Receivers
    • G01S19/14—Receivers specially adapted for specific applications
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S19/00—Satellite radio beacon positioning systems; Determining position, velocity or attitude using signals transmitted by such systems
    • G01S19/01—Satellite radio beacon positioning systems transmitting time-stamped messages, e.g. GPS [Global Positioning System], GLONASS [Global Orbiting Navigation Satellite System] or GALILEO
    • G01S19/13—Receivers
    • G01S19/35—Constructional details or hardware or software details of the signal processing chain
    • G01S19/37—Hardware or software details of the signal processing chain
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S19/00—Satellite radio beacon positioning systems; Determining position, velocity or attitude using signals transmitted by such systems
    • G01S19/38—Determining a navigation solution using signals transmitted by a satellite radio beacon positioning system
    • G01S19/39—Determining a navigation solution using signals transmitted by a satellite radio beacon positioning system the satellite radio beacon positioning system transmitting time-stamped messages, e.g. GPS [Global Positioning System], GLONASS [Global Orbiting Navigation Satellite System] or GALILEO
    • G01S19/42—Determining position

Abstract

本发明提供了一种基于GNSS定位信息的低轨卫星轨道误差修正方法及系统，该轨道误差修正方法包括：建立低轨卫星的运动方程；获取低轨卫星的TLE文件；利用GNSS接收机捕获卫星导航信号以获取本地位置矢量和本地速度矢量；利用低轨卫星接收机量测得到低轨卫星与低轨卫星接收机之间的多普勒频移量测值；根据运动方程建立以低轨卫星的位置矢量和速度矢量为状态量的状态方程；根据TLE文件、本地位置矢量、本地速度矢量以及多普勒频移量测值建立量测方程；基于状态方程和量测方程通过卡尔曼滤波对低轨卫星的轨道误差进行修正。应用本发明的技术方案，以解决现有技术中通过TLE文件计算得出的低轨卫星轨道信息误差大导致定位精度低的技术问题。

Description

基于GNSS定位信息的低轨卫星轨道误差修正方法及系统
技术领域
本发明涉及机会信号导航技术领域，尤其涉及一种基于GNSS定位信息的低轨卫星轨道误差修正方法及系统。
背景技术
当前，卫星导航应用领域广泛，涉及到国民生产生活的方方面面，同时也是各种精确制导武器的重要组成部分。传统导航卫星均分布在中高轨道上，信号到达地面时已经非常微弱，加上信号格式都是公开的，很容易遭受到干扰和欺骗。相比较之下，低轨卫星由于其轨道高度较低，信号到达地面的功率更强，且卫星的几何位置变化快，因此，目前各国争相竞争发射低轨卫星。
利用低轨卫星进行定位目前已经得到较为广泛的应用，比如美国铱星通信公司提供的卫星授时与定位系统(STL)已因具备全球定位系统(GPS)的备份功能而被使用。但是低轨卫星的主要通途是用来通信，并不会像GNSS(全球导航卫星系统)一样实时播发导航电文，而且对于非合作的卫星其通信内容和格式都是经过加密的，因此，得到低轨卫星轨道信息的方式大多依赖于北美航空司令部发布的两行根数文件(TLE文件)。虽然TLE文件每天都会进行更新，但是根据TLE文件计算出来的卫星轨道具有非常大的误差，该误差甚至是公里级别的，对定位结果的影响非常大，减小TLE文件带来的轨道误差有助于提高低轨卫星的定位精度。
发明内容
为了解决现有技术中存在的问题之一，本发明提供了一种基于GNSS定位信息的低轨卫星轨道误差修正方法及系统。
根据本发明的一方面，提供了一种基于GNSS定位信息的低轨卫星轨道误差修正方法，轨道误差修正方法包括：
建立低轨卫星的运动方程；
获取低轨卫星的TLE文件；
利用GNSS接收机捕获卫星导航信号以获取本地位置矢量和本地速度矢量；
利用低轨卫星接收机量测得到低轨卫星与低轨卫星接收机之间的多普勒频移量测值；
根据运动方程建立以低轨卫星的位置矢量和速度矢量为状态量的状态方程；
根据TLE文件、本地位置矢量、本地速度矢量以及多普勒频移量测值建立量测方程；
基于状态方程和量测方程通过卡尔曼滤波对低轨卫星的轨道误差进行修正。
进一步地，根据TLE文件、本地位置矢量、本地速度矢量以及多普勒频移量测值建立量测方程包括：
根据TLE文件、本地位置矢量和本地速度矢量预测得到低轨卫星与低轨卫星接收机之间的多普勒频移预测值；
根据多普勒频移量测值和多普勒频移预测值建立量测方程。
进一步地，根据TLE文件、本地位置矢量和本地速度矢量预测得到低轨卫星与低轨卫星接收机之间的多普勒频移预测值包括：
根据TLE文件预测得到低轨卫星的预测位置矢量和预测速度矢量；
根据预测位置矢量和本地位置矢量计算得到低轨卫星到低轨卫星接收机的相对位置矢量；
根据预测速度矢量和本地速度矢量计算得到低轨卫星到低轨卫星接收机的相对速度矢量；
根据相对位置矢量和相对速度矢量计算得到多普勒频移预测值。
进一步地，根据多普勒频移量测值和多普勒频移预测值建立量测方程包括：
计算多普勒频移预测值与多普勒频移量测值之间的多普勒频移差值；
建立以多普勒频移差值为量测量的量测方程。
进一步地，通过以下公式根据相对位置矢量和相对速度矢量计算得到多普勒频移预测值：
Figure BDA0003772175810000031
上式中，fr(k)表示第k个采样点的多普勒频移预测值，f0(k)表示第k个采样点低轨卫星的发射频率，c表示电磁波在空间中的传播速度，n(k)表示第k个采样点低轨卫星的发射频率的量测误差，ρ(k)表示第k个采样点低轨卫星到低轨卫星接收机的相对位置矢量，
Figure BDA0003772175810000032
表示第k个采样点低轨卫星到低轨卫星接收机的相对速度矢量，ρx(k)、ρy(k)和ρz(k)表示相对位置矢量ρ(k)在三方向上的分量，
Figure BDA0003772175810000033
和
Figure BDA0003772175810000034
表示相对速度矢量
Figure BDA0003772175810000035
在三方向上的分量。
进一步地，基于状态方程和量测方程通过卡尔曼滤波对低轨卫星的轨道误差进行修正包括：
基于状态方程和量测方程通过卡尔曼滤波得到低轨卫星的位置误差和速度误差；
根据位置误差和速度误差对预测位置矢量和预测速度矢量进行修正以得到低轨卫星的修正位置矢量和修正速度矢量；
根据修正位置矢量和修正速度矢量转换得到低轨卫星修正轨道误差后的轨道参数。
进一步地，低轨卫星的运动方程为：
Figure BDA0003772175810000041
上式中，x、y和z分别表示低轨卫星在惯性直角坐标系中的位置矢量在三方向上的分量，
Figure BDA0003772175810000042
和
Figure BDA0003772175810000043
分别表示低轨卫星在惯性直角坐标系中的加速度在x、y、z三方向上的分量，μ表示地球开普勒常数，Re表示地球赤道半径，J2表示地球形状力学因子，
Figure BDA0003772175810000044
表示低轨卫星惯性直角坐标系原点的距离。
进一步地，状态方程为：
Figure BDA0003772175810000045
上式中，X(k+1)表示第k+1个采样点的状态量，X(k)表示第k个采样点的状态量，X(t)表示t时刻的状态量，tk表示第k个采样点的时刻，
Figure BDA0003772175810000046
r(k)表示低轨卫星在第k个采样点在惯性直角坐标系中的位置矢量，
Figure BDA0003772175810000047
表示低轨卫星在第k个采样点在惯性直角坐标系中的速度矢量，r(k)＝[xk yk zk]，
Figure BDA0003772175810000048
xk、yk和zk分别表示低轨卫星在第k个采样点在惯性直角坐标系中的位置矢量在三方向上的分量，
Figure BDA0003772175810000049
和
Figure BDA00037721758100000410
分别表示低轨卫星在第k个采样点在惯性直角坐标系中的速度矢量在三方向上的分量，T表示采样间隔，I表示单位矩阵，F表示状态量X的非线性变换矩阵，ω(k)表示系统噪声矩阵。
进一步地，量测方程为：
Figure BDA0003772175810000051
上式中，Z(k)表示第k个采样点的量测量，Δfk表示多普勒频移差值，
Figure BDA0003772175810000052
表示低轨卫星到低轨卫星接收机的距离。
根据本发明的另一方面，提供了一种基于GNSS定位信息的低轨卫星轨道误差修正系统，误差修正系统包括：
数据读取单元，数据读取单元用于读取低轨卫星的TLE文件；
GNSS接收机，GNSS接收机用于捕获卫星导航信号以获取本地位置矢量和本地速度矢量；
低轨卫星接收机，低轨卫星接收机用于量测得到低轨卫星与低轨卫星接收机之间的多普勒频移量测值；
卡尔曼滤波修正单元，卡尔曼滤波单元用于根据建立的低轨卫星运动方程、状态方程和量测方程进行卡尔曼滤波以对低轨卫星的轨道误差进行修正，状态方程根据运动方程建立并且以低轨卫星的位置矢量和速度矢量为状态量，量测方程根据TLE文件、本地位置矢量、本地速度矢量以及多普勒频移量测值建立。
应用本发明的技术方案，提供了一种基于GNSS定位信息的低轨卫星轨道误差修正方法及系统，该方法通过建立低轨卫星运动方程，获取低轨卫星的TLE文件、低轨卫星与低轨卫星接收机之间的多普勒频移量测值以及GNSS定位信息，基于前述建立的运动方程和获取到的信息建立状态方程和量测方程，进而基于该状态方程和量测方程通过卡尔曼滤波对低轨卫星的轨道误差进行修正，能够显著减小低轨卫星的轨道误差，进而大幅提升低轨卫星的定位精度，保证定位误差收敛在一定范围内。
附图说明
所包括的附图用来提供对本发明实施例的进一步的理解，其构成了说明书的一部分，用于例示本发明的实施例，并与文字描述一起来阐释本发明的原理。显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
图1示出了根据本发明的具体实施例提供的基于GNSS定位信息的低轨卫星轨道误差修正方法的流程示意图；
图2示出了根据本发明的具体实施例提供的基于GNSS定位信息的低轨卫星轨道误差修正方法的原理示意图。
具体实施方式
需要说明的是，在不冲突的情况下，本申请中的实施例及实施例中的特征可以相互组合。下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。以下对至少一个示例性实施例的描述实际上仅仅是说明性的，决不作为对本发明及其应用或使用的任何限制。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
需要注意的是，这里所使用的术语仅是为了描述具体实施方式，而非意图限制根据本申请的示例性实施方式。如在这里所使用的，除非上下文另外明确指出，否则单数形式也意图包括复数形式，此外，还应当理解的是，当在本说明书中使用术语“包含”和/或“包括”时，其指明存在特征、步骤、操作、器件、组件和/或它们的组合。
除非另外具体说明，否则在这些实施例中阐述的部件和步骤的相对布置、数字表达式和数值不限制本发明的范围。同时，应当明白，为了便于描述，附图中所示出的各个部分的尺寸并不是按照实际的比例关系绘制的。对于相关领域普通技术人员已知的技术、方法和设备可能不作详细讨论，但在适当情况下，所述技术、方法和设备应当被视为授权说明书的一部分。在这里示出和讨论的所有示例中，任何具体值应被解释为仅仅是示例性的，而不是作为限制。因此，示例性实施例的其它示例可以具有不同的值。应注意到：相似的标号和字母在下面的附图中表示类似项，因此，一旦某一项在一个附图中被定义，则在随后的附图中不需要对其进行进一步讨论。
如图1所示，根据本发明的具体实施例提供了一种基于GNSS定位信息的低轨卫星轨道误差修正方法，轨道误差修正方法包括：
S1，建立低轨卫星的运动方程；
S2，获取低轨卫星的TLE文件；
S3，利用GNSS接收机捕获卫星导航信号以获取本地位置矢量和本地速度矢量；
S4，利用低轨卫星接收机量测得到低轨卫星与低轨卫星接收机之间的多普勒频移量测值；
S5，根据运动方程建立以低轨卫星的位置矢量和速度矢量为状态量的状态方程；
S6，根据TLE文件、本地位置矢量、本地速度矢量以及多普勒频移量测值建立量测方程；
S7，基于状态方程和量测方程通过卡尔曼滤波对低轨卫星的轨道误差进行修正。
其中，GNSS接收机依次捕获卫星导航信号，通过捕获、跟踪、帧同步、位同步等步骤解析出原始电文，通过最小二乘得到定位结果，也就是本地位置矢量和本地速度矢量。本发明实施例中，坐标系采用惯性直角坐标系，例如WGS84惯性坐标系。卡尔曼滤波方式根据实际情况进行选择，例如如图2所示，可以采用拓展卡尔曼滤波，状态量的滤波初值设置为TLE文件解算得到的卫星位置矢量和速度矢量，根据建立的状态方程和量测方程实时解算状态向量，即低轨卫星的位置矢量和速度矢量，并且只有在能够接收到低轨卫星信号(可见卫星的瞬时多普勒频移)时才进行卡尔曼滤波更新。同时，本发明状态量中不包含低轨卫星接收机的钟差和钟漂，时间以GNSS为基准，能够降低运算量以便于工程实现。
应用此种配置方式，提供了一种基于GNSS定位信息的低轨卫星轨道误差修正方法，该方法通过建立低轨卫星运动方程，获取低轨卫星的TLE文件、低轨卫星与低轨卫星接收机之间的多普勒频移量测值以及GNSS定位信息，基于前述建立的运动方程和获取到的信息建立状态方程和量测方程，进而基于该状态方程和量测方程通过卡尔曼滤波对低轨卫星的轨道误差进行修正，能够显著减小低轨卫星的轨道误差，进而大幅提升低轨卫星的定位精度，保证定位误差收敛在一定范围内。与现有技术相比，本发明的技术方案能够解决现有技术中通过TLE文件计算得出的低轨卫星轨道信息误差大导致定位精度低的技术问题。
本发明实施例中，状态量X的表达式为：
Figure BDA0003772175810000091
其中，r＝[x y z]和
Figure BDA0003772175810000092
分别代表WGS84惯性坐标系下低轨卫星的位置矢量和速度矢量。
低轨卫星的状态微分方程为：
Figure BDA0003772175810000093
其中，F表示状态量X的非线性变换矩阵。
根据万有引力定律，低轨卫星的运动方程为：
Figure BDA0003772175810000094
进一步地，低轨卫星的运动方程为考虑J2摄动项的方程，即低轨卫星的运动方程为：
Figure BDA0003772175810000095
上式中，x、y和z分别表示低轨卫星在惯性直角坐标系中的位置矢量在三方向上的分量，
Figure BDA0003772175810000096
和
Figure BDA0003772175810000097
分别表示低轨卫星在惯性直角坐标系中的加速度在x、y、z三方向上的分量，μ表示地球开普勒常数，Re表示地球赤道半径，J2表示地球形状力学因子，
Figure BDA0003772175810000101
表示低轨卫星惯性直角坐标系原点的距离。
状态预测方程为：
Figure BDA0003772175810000102
其中，
Figure BDA0003772175810000103
T表示采样间隔，tk表示采样点时刻。
状态转移矩阵离散化后可得：
Figure BDA0003772175810000104
经过线性化和离散化后得到的状态方程为：
X(k+1)＝Φ(k+1/k)·X(k)+ω(k)，
也就是，本发明实施例中，状态方程为：
Figure BDA0003772175810000105
上式中，X(k+1)表示第k+1个采样点的状态量，X(k)表示第k个采样点的状态量，X(t)表示t时刻的状态量，tk表示第k个采样点的时刻，
Figure BDA0003772175810000106
r(k)表示低轨卫星在第k个采样点在惯性直角坐标系中的位置矢量，
Figure BDA0003772175810000107
表示低轨卫星在第k个采样点在惯性直角坐标系中的速度矢量，r(k)＝[xk yk zk]，
Figure BDA0003772175810000108
xk、yk和zk分别表示低轨卫星在第k个采样点在惯性直角坐标系中的位置矢量在三方向上的分量，
Figure BDA0003772175810000109
和
Figure BDA00037721758100001010
分别表示低轨卫星在第k个采样点在惯性直角坐标系中的速度矢量在三方向上的分量，T表示采样间隔，I表示单位矩阵，F表示状态量X的非线性变换矩阵，ω(k)表示系统噪声矩阵，对应的协方差矩阵为Q(k)＝E[ω(k)ωT(k)]。
此外，本发明中，由于GNSS接收机与低轨卫星之间以及低轨卫星接收机与低轨卫星之间均存在相对运动，因此GNSS接收机与低轨卫星之间以及低轨卫星接收机与低轨卫星之间均存在多普勒频移。基于此，本发明实施例中，根据TLE文件、本地位置矢量、本地速度矢量以及多普勒频移量测值建立量测方程包括：
根据TLE文件、本地位置矢量和本地速度矢量预测得到低轨卫星与低轨卫星接收机之间的多普勒频移预测值；
根据多普勒频移量测值和多普勒频移预测值建立量测方程。
基于上述实施例，本发明实施例中，根据TLE文件、本地位置矢量和本地速度矢量预测得到低轨卫星与低轨卫星接收机之间的多普勒频移预测值包括：
根据TLE文件预测得到低轨卫星的预测位置矢量和预测速度矢量；
根据预测位置矢量和本地位置矢量计算得到低轨卫星到低轨卫星接收机的相对位置矢量；
根据预测速度矢量和本地速度矢量计算得到低轨卫星到低轨卫星接收机的相对速度矢量；
根据相对位置矢量和相对速度矢量计算得到多普勒频移预测值。
进一步地，作为本发明的一个具体实施例，通过以下公式根据相对位置矢量和相对速度矢量计算得到多普勒频移预测值：
Figure BDA0003772175810000111
上式中，fr(k)表示第k个采样点的多普勒频移预测值，f0(k)表示第k个采样点低轨卫星的发射频率，c表示电磁波在空间中的传播速度，n(k)表示第k个采样点低轨卫星的发射频率的量测误差，ρ(k)表示第k个采样点低轨卫星到低轨卫星接收机的相对位置矢量，
Figure BDA0003772175810000121
表示第k个采样点低轨卫星到低轨卫星接收机的相对速度矢量，ρx(k)、ρy(k)和ρz(k)表示相对位置矢量ρ(k)在三方向上的分量，
Figure BDA0003772175810000122
和
Figure BDA0003772175810000123
表示相对速度矢量
Figure BDA0003772175810000124
在三方向上的分量，其中，n(k)一般为白噪声，设均值为零，方差为σ2，k＝0,1,…,N，N表示采样点的数量，采样点也就是观测点。
此外，本发明实施例中，根据多普勒频移量测值和多普勒频移预测值建立量测方程包括：
计算多普勒频移预测值与多普勒频移量测值之间的多普勒频移差值；
建立以多普勒频移差值为量测量的量测方程。
作为本发明的一个具体实施例，量测量Z(k)的表达式为Z(k)＝Δfk＝fr(k)-fc(k)，其中，Δfk表示多普勒频移预测值与多普勒频移量测值之间的多普勒频移差值，也就是，多普勒频移预测值减去多普勒频移量测值，fc(k)表示多普勒频移量测值。基于该实施例，量测量Z(k)表示为状态量X的非线性函数为：
Z(k)＝H(X(k))+n(k)，
其中，H(X(k))表示量测量对状态量的非线性变换矩阵，n(k)对应的协方差矩阵为R(k)＝E[n(k)nT(k)]＝σ2。
进一步地，量测雅克比矩阵为：
Figure BDA0003772175810000125
其中的偏微分为：
Figure BDA0003772175810000126
Figure BDA0003772175810000131
其中
Figure BDA0003772175810000132
表示低轨卫星到低轨卫星接收机的距离。
经过线性化和离散化后得到的量测方程为：
Z(k)＝H(k+1/k)·X(k)+n(k)，
也就是，本发明实施例中，量测方程为：
Figure BDA0003772175810000133
上式中，Z(k)表示第k个采样点的量测量，Δfk表示多普勒频移差值，
Figure BDA0003772175810000134
表示低轨卫星到低轨卫星接收机的距离。
进一步地，将前述构建的离散化和线性化后的状态方程和量测方程带入卡尔曼滤波工时进行计算，从而对低轨卫星的轨道误差进行修正。本发明实施例中，基于状态方程和量测方程通过卡尔曼滤波对低轨卫星的轨道误差进行修正包括：
基于状态方程和量测方程通过卡尔曼滤波得到低轨卫星的位置误差和速度误差；
根据位置误差和速度误差对预测位置矢量和预测速度矢量进行修正以得到低轨卫星的修正位置矢量和修正速度矢量；
根据修正位置矢量和修正速度矢量转换得到低轨卫星修正轨道误差后的轨道参数。
其中，滤波估计步骤如下：
状态一步预测均方误差阵：
P(k+1/k)＝Φ(k+1/k)·P(k)·ΦT(k+1/k)+Q，
其中Q阵中的参数根据实际情况合理选择。
K(k+1)＝P(k+1/k)·H(k+1/k)·(H(k+1/k)·P(k)·HT(k+1/k)+R(k+1))，
状态估计：
Figure BDA0003772175810000141
状态估计均方误差阵：
P(k+1)＝(I-K(k+1)·H(k+1))·P(k+1/k)，
滤波估计过程为本领域技术人员所熟知，在此不再赘述，通过滤波估计即可得到低轨卫星的位置误差和速度误差，根据位置误差和速度误差对预测位置矢量和预测速度矢量进行修正就能够得到低轨卫星的修正位置矢量和修正速度矢量，再将修正位置矢量和修正速度矢量转换为轨道参数就得到了轨道误差修正后的轨道参数。本发明实施例中，轨道参数主要包括轨道倾角、升交点赤经、近地点角距、轨道半长轴、纬度幅角和真近点角这六个参数，轨道参数转换的具体过程如下：
卫星的面积积分公式
Figure BDA0003772175810000142
其中h为积分常矢量，其在惯性直角坐标系中的表达为h＝(hx,hy,hz)，则
Figure BDA0003772175810000143
其中，h表示积分常矢量h的长度，hx、hy和hz分别表示积分常矢量h在三方向上的分量。
在以h为z'轴，r方向为x'轴的轨道坐标系中有h＝(0,0,h)，再由惯性直角坐标系和轨道直角坐标系的转换关系可得：
Figure BDA0003772175810000144
上式中，i表示轨道倾角，Ω表示升交点赤经，两者通过以下公式即可得到：
Figure BDA0003772175810000151
Figure BDA0003772175810000152
进一步地，由
Figure BDA0003772175810000153
可得：
Figure BDA0003772175810000154
上式中，e表示轨道偏心率，其在惯性直角坐标系中的坐标为e＝(ex,ey,ez)，则
Figure BDA0003772175810000155
其中，ex、ey和ez分别表示轨道偏心率矢量在坐标轴三个方向投影。
接下来，可以计算近地点角距w，其计算公式为：
Figure BDA0003772175810000156
由圆锥曲线的公式可得轨道半长轴a的计算公式为：
Figure BDA0003772175810000157
再按下式计算纬度幅角u，其定义为在轨道平面内由升交点到低轨卫星位置相对地心的角距：
Figure BDA0003772175810000158
真近点角f的计算公式为：
f＝u-w，
此外，偏近点角与真近点角的关系为
Figure BDA0003772175810000159
再由开普勒方程可求平近点角M，平近点角是指在轨道平面内从近地点起算低轨卫星以平均角速度运行的角度，其计算公式为M＝E-e sin E，其中，E表示偏近点角。
根据本发明的另一方面，提供了一种基于GNSS定位信息的低轨卫星轨道误差修正系统，误差修正系统包括：
数据读取单元，数据读取单元用于读取低轨卫星的TLE文件；
GNSS接收机，GNSS接收机用于捕获卫星导航信号以获取本地位置矢量和本地速度矢量；
低轨卫星接收机，低轨卫星接收机用于量测得到低轨卫星与低轨卫星接收机之间的多普勒频移量测值；
卡尔曼滤波修正单元，卡尔曼滤波单元用于根据建立的低轨卫星运动方程、状态方程和量测方程进行卡尔曼滤波以对低轨卫星的轨道误差进行修正，状态方程根据运动方程建立并且以低轨卫星的位置矢量和速度矢量为状态量，量测方程根据TLE文件、本地位置矢量、本地速度矢量以及多普勒频移量测值建立。
关于该卫星轨道误差修正系统的相关示例性描述请参考前述关于卫星轨道误差修正方法的相关示例性描述，在此不再赘述，利用该系统能够显著减小低轨卫星的轨道误差，进而大幅提升低轨卫星的定位精度，保证定位误差收敛在一定范围内。
进一步地，基于本发明前述提出的低轨卫星轨道误差修正方法或系统，本发明还提供一种组合导航定位方法及系统。具体来说，在利用本发明前述提出的低轨卫星轨道误差修正方法或系统对轨道误差修正后，利用修正后的修正轨道参数以及根据TLE文件解算出来的原始轨道参数作为神经网络的训练集，通过训练得到低轨卫星轨道参数的预测模型，利用该预测模型对轨道参数进行预测以得到预测轨道参数，进而得到由低轨卫星的预测轨道参数构成的轨道数据库，从而在需要低轨卫星导航时直接调用该数据库中的预测轨道参数。具体请参考图2的实施例，组合导航定位系统包括GNSS接收机和低轨卫星接收机，在GNSS定位信息可用的情况下，组合导航定位系统工作在GNSS定位模式下，并且利用本发明前述提出的轨道误差修正方法或系统对低轨卫星的轨道误差进行实时修正，同时通过神经网络算法利用修正后的修正轨道参数和TLE文件解算出的原始轨道参数对低轨卫星的轨道参数进行实时预测以得到预测轨道参数(时间序列)，存储于轨道数据库中，一旦GNSS定位结果不可靠或GNSS接收机受到遮挡或者干扰而无法提取GNSS导航观测量，则组合导航定位系统切换至低轨卫星定位模式，从前述的轨道数据库中提取预测的当前时刻的低轨卫星轨道参数，通过预测轨道参数求解得到低轨卫星当前的位置矢量和速度矢量，再结合低轨卫星的多普勒观测量也即多普勒频移量测量进行定位，从而保证定位误差收敛在一定范围内。其中，神经网络算法根据实际情况选用，作为本发明的一个具体实施例，采用LSTM预测算法进行训练。该导航方法及对应的导航系统能够显著提升导航的精度和可靠性。
进一步地，运用LSTM算法进行低轨卫星轨道参数(时间序列)预测时，神经网络的输入为轨道六参数，包括根据TLE文件解算的原始轨道六参数和修正后的轨道六参数，其中原始轨道六参数作为先验信息。训练过程如下：首先初始化网络，设置偏差值，通过前向运算得到输出值，之后确定目标函数，计算实际值与估计值的误差，构造误差函数，最终根据梯度下降原则更新网络的权重和偏置。神经网络算法训练的过程为神经网络领域的技术人员所熟知，在此不再赘述。
综上所述，本发明提供了一种基于GNSS定位信息的低轨卫星轨道误差修正方法及系统，该方法通过建立低轨卫星运动方程，获取低轨卫星的TLE文件、低轨卫星与低轨卫星接收机之间的多普勒频移量测值以及GNSS定位信息，基于前述建立的运动方程和获取到的信息建立状态方程和量测方程，进而基于该状态方程和量测方程通过卡尔曼滤波对低轨卫星的轨道误差进行修正，能够显著减小低轨卫星的轨道误差，进而大幅提升低轨卫星的定位精度，保证定位误差收敛在一定范围内。与现有技术相比，本发明的技术方案能够解决现有技术中通过TLE文件计算得出的低轨卫星轨道信息误差大导致定位精度低的技术问题。
为了便于描述，在这里可以使用空间相对术语，如“在……之上”、“在……上方”、“在……上表面”、“上面的”等，用来描述如在图中所示的一个器件或特征与其他器件或特征的空间位置关系。应当理解的是，空间相对术语旨在包含除了器件在图中所描述的方位之外的在使用或操作中的不同方位。例如，如果附图中的器件被倒置，则描述为“在其他器件或构造上方”或“在其他器件或构造之上”的器件之后将被定位为“在其他器件或构造下方”或“在其他器件或构造之下”。因而，示例性术语“在……上方”可以包括“在……上方”和“在……下方”两种方位。该器件也可以其他不同方式定位(旋转90度或处于其他方位)，并且对这里所使用的空间相对描述作出相应解释。
此外，需要说明的是，使用“第一”、“第二”等词语来限定零部件，仅仅是为了便于对相应零部件进行区别，如没有另行声明，上述词语并没有特殊含义，因此不能理解为对本发明保护范围的限制。
以上所述仅为本发明的优选实施例而已，并不用于限制本发明，对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围内。

Claims (10)

1.一种基于GNSS定位信息的低轨卫星轨道误差修正方法，其特征在于，所述轨道误差修正方法包括：
建立所述低轨卫星的运动方程；
获取所述低轨卫星的TLE文件；
利用GNSS接收机捕获卫星导航信号以获取本地位置矢量和本地速度矢量；
利用低轨卫星接收机量测得到所述低轨卫星与所述低轨卫星接收机之间的多普勒频移量测值；
根据所述运动方程建立以所述低轨卫星的位置矢量和速度矢量为状态量的状态方程；
根据所述TLE文件、所述本地位置矢量、所述本地速度矢量以及所述多普勒频移量测值建立量测方程；
基于所述状态方程和所述量测方程通过卡尔曼滤波对所述低轨卫星的轨道误差进行修正。
2.根据权利要求1所述的轨道误差修正方法，其特征在于，根据所述TLE文件、所述本地位置矢量、所述本地速度矢量以及所述多普勒频移量测值建立量测方程包括：
根据所述TLE文件、所述本地位置矢量和所述本地速度矢量预测得到所述低轨卫星与所述低轨卫星接收机之间的多普勒频移预测值；
根据所述多普勒频移量测值和所述多普勒频移预测值建立量测方程。
3.根据权利要求2所述的轨道误差修正方法，其特征在于，根据所述TLE文件、所述本地位置矢量和所述本地速度矢量预测得到所述低轨卫星与所述低轨卫星接收机之间的多普勒频移预测值包括：
根据所述TLE文件预测得到所述低轨卫星的预测位置矢量和预测速度矢量；
根据所述预测位置矢量和所述本地位置矢量计算得到所述低轨卫星到所述低轨卫星接收机的相对位置矢量；
根据所述预测速度矢量和所述本地速度矢量计算得到所述低轨卫星到所述低轨卫星接收机的相对速度矢量；
根据所述相对位置矢量和所述相对速度矢量计算得到所述多普勒频移预测值。
4.根据权利要求3所述的轨道误差修正方法，其特征在于，根据所述多普勒频移量测值和所述多普勒频移预测值建立量测方程包括：
计算所述多普勒频移预测值与所述多普勒频移量测值之间的多普勒频移差值；
建立以所述多普勒频移差值为量测量的量测方程。
5.根据权利要求4所述的轨道误差修正方法，其特征在于，通过以下公式根据所述相对位置矢量和所述相对速度矢量计算得到所述多普勒频移预测值：
Figure FDA0003772175800000021
上式中，fr(k)表示第k个采样点的多普勒频移预测值，f0(k)表示第k个采样点低轨卫星的发射频率，c表示电磁波在空间中的传播速度，n(k)表示第k个采样点低轨卫星的发射频率的量测误差，ρ(k)表示第k个采样点低轨卫星到低轨卫星接收机的相对位置矢量，
Figure FDA0003772175800000022
表示第k个采样点低轨卫星到低轨卫星接收机的相对速度矢量，ρx(k)、ρy(k)和ρz(k)表示所述相对位置矢量ρ(k)在三方向上的分量，
Figure FDA0003772175800000031
和
Figure FDA0003772175800000032
表示所述相对速度矢量
Figure FDA0003772175800000033
在三方向上的分量。
6.根据权利要求5所述的轨道误差修正方法，其特征在于，基于所述状态方程和所述量测方程通过卡尔曼滤波对所述低轨卫星的轨道误差进行修正包括：
基于所述状态方程和所述量测方程通过卡尔曼滤波得到所述低轨卫星的位置误差和速度误差；
根据所述位置误差和所述速度误差对所述预测位置矢量和所述预测速度矢量进行修正以得到所述低轨卫星的修正位置矢量和修正速度矢量；
根据所述修正位置矢量和所述修正速度矢量转换得到所述低轨卫星修正轨道误差后的轨道参数。
7.根据权利要求6所述的轨道误差修正方法，其特征在于，所述低轨卫星的运动方程为：
Figure FDA0003772175800000034
上式中，x、y和z分别表示所述低轨卫星在惯性直角坐标系中的位置矢量在三方向上的分量，
Figure FDA0003772175800000035
和
Figure FDA0003772175800000036
分别表示所述低轨卫星在惯性直角坐标系中的加速度在x、y、z三方向上的分量，μ表示地球开普勒常数，Re表示地球赤道半径，J2表示地球形状力学因子，
Figure FDA0003772175800000037
表示低轨卫星惯性直角坐标系原点的距离。
8.根据权利要求7所述的轨道误差修正方法，其特征在于，所述状态方程为：
Figure FDA0003772175800000041
上式中，X(k+1)表示第k+1个采样点的状态量，X(k)表示第k个采样点的状态量，X(t)表示t时刻的状态量，tk表示第k个采样点的时刻，
Figure FDA0003772175800000042
r(k)表示所述低轨卫星在第k个采样点在惯性直角坐标系中的位置矢量，
Figure FDA0003772175800000043
表示所述低轨卫星在第k个采样点在惯性直角坐标系中的速度矢量，r(k)＝[xk yk zk]，
Figure FDA0003772175800000044
xk、yk和zk分别表示所述低轨卫星在第k个采样点在惯性直角坐标系中的位置矢量在三方向上的分量，
Figure FDA0003772175800000045
和
Figure FDA0003772175800000046
分别表示所述低轨卫星在第k个采样点在惯性直角坐标系中的速度矢量在三方向上的分量，T表示采样间隔，I表示单位矩阵，F表示状态量X的非线性变换矩阵，ω(k)表示系统噪声矩阵。
9.根据权利要求8所述的轨道误差修正方法，其特征在于，所述量测方程为：
Figure FDA0003772175800000047
上式中，Z(k)表示第k个采样点的量测量，Δfk表示所述多普勒频移差值，
Figure FDA0003772175800000048
表示所述低轨卫星到所述低轨卫星接收机的距离。
10.一种基于GNSS定位信息的低轨卫星轨道误差修正系统，其特征在于，所述误差修正系统包括：
数据读取单元，所述数据读取单元用于读取所述低轨卫星的TLE文件；
GNSS接收机，所述GNSS接收机用于捕获卫星导航信号以获取本地位置矢量和本地速度矢量；
低轨卫星接收机，所述低轨卫星接收机用于量测得到所述低轨卫星与所述低轨卫星接收机之间的多普勒频移量测值；
卡尔曼滤波修正单元，所述卡尔曼滤波单元用于根据建立的低轨卫星运动方程、状态方程和量测方程进行卡尔曼滤波以对所述低轨卫星的轨道误差进行修正，所述状态方程根据所述运动方程建立并且以所述低轨卫星的位置矢量和速度矢量为状态量，所述量测方程根据所述TLE文件、所述本地位置矢量、所述本地速度矢量以及所述多普勒频移量测值建立。
CN202210905416.2A 2022-07-29 2022-07-29 基于gnss定位信息的低轨卫星轨道误差修正方法及系统 Pending CN115327587A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210905416.2A CN115327587A (zh) 2022-07-29 2022-07-29 基于gnss定位信息的低轨卫星轨道误差修正方法及系统

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210905416.2A CN115327587A (zh) 2022-07-29 2022-07-29 基于gnss定位信息的低轨卫星轨道误差修正方法及系统

Publications (1)

Publication Number Publication Date
CN115327587A true CN115327587A (zh) 2022-11-11

Family

ID=83919216

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210905416.2A Pending CN115327587A (zh) 2022-07-29 2022-07-29 基于gnss定位信息的低轨卫星轨道误差修正方法及系统

Country Status (1)

Country Link
CN (1) CN115327587A (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116552812A (zh) * 2023-04-12 2023-08-08 四川大学 一种电推进geo卫星自学习定轨方法

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116552812A (zh) * 2023-04-12 2023-08-08 四川大学 一种电推进geo卫星自学习定轨方法
CN116552812B (zh) * 2023-04-12 2024-01-23 四川大学 一种电推进geo卫星自学习定轨方法

Similar Documents

Publication Publication Date Title
Fang et al. Predictive iterated Kalman filter for INS/GPS integration and its application to SAR motion compensation
CN104406605B (zh) 机载多导航源综合导航仿真系统
US6424914B1 (en) Fully-coupled vehicle positioning method and system thereof
US6496778B1 (en) Real-time integrated vehicle positioning method and system with differential GPS
Gebre-Egziabher Design and performance analysis of a low-cost aided dead reckoning navigator
Hasan et al. A review of navigation systems (integration and algorithms)
CN110779518B (zh) 一种具有全局收敛性的水下航行器单信标定位方法
CN104678408A (zh) 星载导航接收机授时方法和授时型星载导航接收机以及星载导航应用系统
CN111380518B (zh) 一种引入径向速度的sins/usbl紧组合导航定位方法
CN115356754A (zh) 一种基于gnss和低轨卫星的组合导航定位方法
CN111965685B (zh) 一种基于多普勒信息的低轨卫星/惯性组合导航定位方法
CN113203418B (zh) 基于序贯卡尔曼滤波的gnssins视觉融合定位方法及系统
RU2314553C1 (ru) Система оценки точностных характеристик бортовой радиолокационной станции
Gebre-Egziabher et al. MAV attitude determination by vector matching
CN115327587A (zh) 基于gnss定位信息的低轨卫星轨道误差修正方法及系统
Ifqir et al. Fault tolerant multi-sensor data fusion for autonomous navigation in future civil aviation operations
Iiyama et al. Terrestrial GPS time-differenced carrier-phase positioning of lunar surface users
CN109471102B (zh) 一种惯组误差修正方法
Spiridonov et al. Orbit determination methods for leo satellites from probabilistic analysis, circular motion model and single pass doppler measurements
Ward Spacecraft attitude estimation using GPS: Methodology and results
İpek Satellite orbit estimation using kalman filters
CN113324539A (zh) 一种sins/srs/cns多源融合自主组合导航方法
Runnalls et al. Terrain-referenced navigation using the IGMAP data fusion algorithm
CN113551669A (zh) 基于短基线的组合导航定位方法及装置
JPH01500714A (ja) 分散型カルマンフィルタ

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination