CN118129767B - 一种动力学约束的低轨卫星gnss实时精密定轨方法及系统 - Google Patents

一种动力学约束的低轨卫星gnss实时精密定轨方法及系统 Download PDF

Info

Publication number
CN118129767B
CN118129767B CN202410552556.5A CN202410552556A CN118129767B CN 118129767 B CN118129767 B CN 118129767B CN 202410552556 A CN202410552556 A CN 202410552556A CN 118129767 B CN118129767 B CN 118129767B
Authority
CN
China
Prior art keywords
epoch
low
matrix
equation
orbit satellite
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.)
Active
Application number
CN202410552556.5A
Other languages
English (en)
Other versions
CN118129767A (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.)
Wuhan University WHU
Original Assignee
Wuhan University WHU
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 Wuhan University WHU filed Critical Wuhan University WHU
Priority to CN202410552556.5A priority Critical patent/CN118129767B/zh
Publication of CN118129767A publication Critical patent/CN118129767A/zh
Application granted granted Critical
Publication of CN118129767B publication Critical patent/CN118129767B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01C—MEASURING DISTANCES, LEVELS OR BEARINGS; SURVEYING; NAVIGATION; GYROSCOPIC INSTRUMENTS; PHOTOGRAMMETRY OR VIDEOGRAMMETRY
    • G01C21/00—Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00
    • G01C21/24—Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00 specially adapted for cosmonautical navigation
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01C—MEASURING DISTANCES, LEVELS OR BEARINGS; SURVEYING; NAVIGATION; GYROSCOPIC INSTRUMENTS; PHOTOGRAMMETRY OR VIDEOGRAMMETRY
    • G01C21/00—Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00
    • G01C21/10—Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00 by using measurements of speed or acceleration
    • G01C21/12—Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00 by using measurements of speed or acceleration executed aboard the object being navigated; Dead reckoning
    • 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

Landscapes

  • Engineering & Computer Science (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Remote Sensing (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Signal Processing (AREA)
  • Automation & Control Theory (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • Astronomy & Astrophysics (AREA)
  • Control Of Position, Course, Altitude, Or Attitude Of Moving Bodies (AREA)

Abstract

本发明公开了一种动力学约束的低轨卫星GNSS实时精密定轨方法及系统，根据低轨卫星运动特性，建立连续历元卫星坐标与瞬时加速度间关系的动力学约束，利用地球重力场模型获得满足精度要求的瞬时加速度，将动力学约束虚拟观测方程形成法方程，与GNSS观测方程形成的法方程叠加，利用序贯最小二乘方法更新法方程，逐历元求解待估参数，得到高精度的定轨结果。本发明相比简化动力学定轨方法算力要求更低，同时定轨精度和稳定性优于未引入动力学信息的运动学方法。本发明解决了现有低轨卫星实时精密定轨时轨道精度、稳定性和算力消耗难以兼顾的问题，具有理论严密、模型简单、易于实现等特点，能够服务于低轨导航增强、自动驾驶等多个领域。

Description

一种动力学约束的低轨卫星GNSS实时精密定轨方法及系统
技术领域
本发明属于卫星导航技术领域，具体涉及一种动力学约束的低轨卫星GNSS实时精密定轨方法及系统。
背景技术
低轨卫星（Low Earth Orbiter，LEO）是一类重要的人造卫星，在通信、导航、遥感和其他空间科学研究中具有广泛应用，GNSS是其精密轨道确定的主要手段。面向未来L4、L5级自动驾驶等高精度实时应用场景时，低轨导航增强系统对卫星轨道精度和稳定性提出了更高的要求。然而，受限于定轨模式、GNSS精密产品可获取性、卫星平台可分配算力等因素，低轨卫星在GNSS实时精密定轨中存在轨道精度、稳定性和算力消耗之间难以兼顾的矛盾，制约了卫星的高精度服务性能。
低轨卫星GNSS实时精密定轨通常采用运动学定轨和简化动力学定轨两种方式。实时运动学定轨模型简单，求解方便，算力要求低，但定轨结果受制于GNSS卫星几何构型，轨道稳定性较差。实时简化动力学定轨引入低轨卫星动力学信息，模型强度更强，轨道稳定性更好，但相应的模型和求解方法更复杂，算力要求更高。在低轨卫星星上设备算力限制下，使用上述两种方法难以有效获得准确稳定的低轨卫星精密轨道。因此，发展一种精度高、稳定性强、算力要求低的动力学约束低轨卫星实时精密定轨方法，既能利用低轨卫星运动特性增强定轨模型，又能满足星上设备算力的限制，对未来自动驾驶及其他高精度位置服务具有重要意义。
发明内容
针对现有技术的不足，本发明提供一种动力学约束的低轨卫星GNSS实时精密定轨方法，既能提升低轨卫星轨道精度和稳定性，又能降低设备算力要求，包括以下步骤：
步骤1，构建能够同时估计连续三个历元低轨卫星坐标的序贯最小二乘参数估计器，作为动力学约束的低轨卫星GNSS实时定轨参数估计器；
步骤2，利用低轨卫星运动特性，建立连续三个历元位置与瞬时加速度的关系，得到动力学约束；
步骤3，利用地球重力模型获取高精度瞬时加速度；
步骤4，将动力学约束带入步骤1构建的低轨卫星GNSS实时定轨参数估计器，求解高精度低轨卫星坐标。
而且，所述步骤1中动力学约束的低轨卫星GNSS实时精密定轨首先需要构造适配的参数估计器，对于任意历元，低轨卫星GNSS观测方程描述为：
（1）
式中，下标表示第历元，表示第历元的观测值，表示第历元的设计矩阵，表示第历元的待估参数矩阵，、、表示第历元低轨卫星坐标，、、表示GNSS卫星到低轨卫星向量关于、、的方向余弦，表示接收机钟差参数，...表示无电离层组合模糊度参数，表示观测噪声。
联立第、、历元的观测方程、、，统一表示为观测方程：
（2）
式中，、、表示第、、历元的GNSS卫星到低轨卫星的方向余弦，、、、表示第、、历元低轨卫星的坐标和接收机钟差，...表示无电离层组合模糊度，表示第、、历元观测噪声，矩阵化后表示为、和。
观测值减去参数初值的OMC矩阵表示为：
（3）
式中，、、表示第、、历元低轨卫星的观测值，、、表示第、、历元低轨卫星初始坐标，、、表示第、、历元设计矩阵，基于低轨卫星初始坐标计算获得。
对式（3）根据最小二乘准则形成法方程，求解待估参数向量：
（4）
式中，权矩阵由设定的观测方程先验精度确定，、表示由、、三个历元对应的参数组成的法方程矩阵和残差矩阵。
待估参数向量中包含第历元的低轨卫星坐标、、，从而实现低轨卫星的精密定轨。
利用序贯最小二乘方法更新法方程，逐历元求解待估参数，具体流程为：①第历元，由式（4）形成第历元法方程，估计第历元低轨卫星坐标及其他参数；②第历元，将法方程中第历元对应低轨卫星坐标和接收机钟差参数通过参数消去方法消除，第和历元对应的法方程信息移至、历元对应位置，空出的历元对应位置填入历元信息，形成第历元法方程，估计第历元低轨卫星坐标及其他参数；③序贯执行步骤②，实现低轨卫星实时定轨。
步骤①、②具体过程推导如下：
①第历元法方程矩阵和残差矩阵写作：
（5）
（6）
其中：
（7）
（8）
（9）
（10）
（11）
（12）
式（7）-式（10）中，元素表示的方差，元素表示 的协方差，，；式（11）-式（12）中，元素表示的残差，。
式（5）为式（4）中法方程矩阵的展开形式，其中历元对应低轨卫星坐标、、及接收机钟差的法方程矩阵表示为，如式（7）所示，包含参数自身法方程系数项与参数间关系项；前序两个历元法方程矩阵、替换式（7）下标就能获得；不同历元参数间关系项，以式（8）为例，替换下标能够全部获得；无电离层组合模糊度参数对应法方程矩阵为，如式（9）所示；无电离层组合模糊度参数与其他参数间关系对应法方程矩阵为，如式（10）所示，下标替换能够获得全部项；式（6）为式（4）中残差矩阵的展开形式，其中历元对应低轨卫星坐标、、及接收机钟差的残差矩阵表示为，如式（11）所示；前序两个历元残差矩阵、替换式（11）下标就能获得，无电离层组合模糊度参数对应残差矩阵为，如式（12）所示；
②第历元，将法方程矩阵和残差矩阵中第历元对应低轨卫星坐标、、和接收机钟差参数通过参数消去方法消除，然后将第和历元对应的法方程信息移至、历元对应位置，并填入第历元信息，得历元的法方程矩阵和残差矩阵：
（13）
（14）
其中，法方程矩阵中元素和残差矩阵中元素描述为：
（15）
（16）
式中，为法方程矩阵中元素，为残差矩阵中元素。
而且，所述步骤2中基于牛顿运动方程，低轨卫星连续历元的平均加速度与相应历元卫星坐标描述为：
（17）
（18）
式中，为低轨卫星第历元的平均加速度；为常系数矩阵；为历元间采样间隔；、、分别为、、历元卫星坐标，为便于推导统一使用惯性系下坐标表示。
将低轨卫星的平均加速度写作多个历元瞬时加速度多项式组合的形式，即：
（19）
式中，为加权函数，为多项式系数矩阵，为低轨卫星瞬时加速度。
多项式阶数与选用的历元数有关，对于、、连续三个历元，使用三阶多项式将系数矩阵推导得到常系数矩阵，即：
（20）
与历元平均加速度相关的瞬时加速度表示为、、，进而式（19）写作：
（21）
由式（17）和式（21），、、连续三个历元平均加速度与瞬时加速度的关系推导为：
（22）
为满足实时解算要求，将、、换元为、、，则式（22）写为：
（23）
由此建立连续三个历元位置与瞬时加速度的关系，即为本发明所述动力学约束。
将动力学约束利用广义最小二乘原理写作虚拟观测值，计算公式为：
（24）
上述公式（17）-（24）均在惯性系下推导，若最终待求低轨卫星坐标为地固系下坐标，需使用对应历元地固系到惯性系的旋转矩阵进行转换，将、、历元的旋转矩阵表示为、和，则：
（25）
为便于下文推导，设虚拟观测方程地固系下设计矩阵为，即：
（26）
式中，为常系数矩阵。
而且，所述步骤3中对于任意历元低轨卫星先验初始坐标，首先将其转换至球坐标：
（27）
式中，是低轨卫星到地心的距离，是极角，是方位角。
在球坐标系中，地球重力场模型由球谐系数和表示，利用的模型阶次指取相应阶数的球谐系数进行引力势计算，引力势的球谐展开公式如下：
（28）
式中，是引力势，是低轨卫星到地心的距离，是极角，是方位角，是引力常数，是地球质量，是地球平均半径，是相关的Legendre函数，和表示其系数，，。
此时低轨卫星受到的引力为引力势的梯度，通过引力势对、、的单位向量、、求偏导获得：
（29）
式中，表示对引力势求偏导。
由于引力，为低轨卫星质量，由此能够计算得到任意历元低轨卫星瞬时加速度；利用式（27）-（29）计算得到、、历元瞬时加速度值、和后，式（24）中仅剩坐标参数为未知参数，后续通过参数估计方法求解。
将步骤2、3中动力学约束虚拟观测方程形成法方程，与步骤1中GNSS观测方程形成的法方程叠加，最终参数估计结果即为施加动力学约束的低轨卫星GNSS实时定轨结果，具体推导过程如下：
基于式（24）-（26）动力学约束形成的虚拟观测方程构建法方程，结合式（4）得到第历元法方程矩阵和残差矩阵：
（30）
式中，为第历元动力学约束方程形成的法方程矩阵，为第历元动力学约束方程形成的残差矩阵，为设计矩阵，由虚拟观测方程的精度决定，为观测值减去参数初值的OMC矩阵。
结合式（3），在此处表示为：
（31）
式中，虚拟观测值设为0；表示地固系下坐标初值；表示瞬时加速度初值，通过式（27）-（29）计算获得，计算时所需低轨卫星坐标即为；历元间采样间隔和常系数矩阵已知，由此获得虚拟观测方程的OMC矩阵。
第历元，动力学约束形成的虚拟观测方程法方程矩阵和残差矩阵由（30）和（31）推导得到，表示为：
（32）
（33）
其中：
（34）
（35）
（36）
式中，表示历元对应低轨卫星坐标、、的法方程矩阵， 替换式（34）标获得，不同历元参数间关系项以式（35）为例，替换下标就能全部获得，表示历元对应低轨卫星坐标、、的残差矩阵，式（34）和式（35）中元素表示的方差，元素表示和的协方差，，，式（36）中元素表示的残差，。
而且，所述步骤4中由法方程叠加原理，将式（13）与式（32）、式（14）与式（33）叠加，最终获得本发明所述动力学约束解，如式（37）所示：
（37）
其中，最终待估参数由叠加后的法方程矩阵和叠加后的残差矩阵计算得到：
（38）
（39）
式中，表示GNSS观测方程对应法方程矩阵中的元素，表示第历元对应的法方程矩阵及其与其他元素间关系，表示动力学约束方程的法方程矩阵中的元素，表示GNSS观测方程的残差矩阵中的元素，表示第历元GNSS观测方程的残差矩阵中对应元素，表示动力学约束方程的残差矩阵中的元素。
本发明还提供一种动力学约束的低轨卫星GNSS实时精密定轨系统，用于实现如上所述的一种动力学约束的低轨卫星GNSS实时精密定轨方法。
而且，包括处理器和存储器，存储器用于存储程序指令，处理器用于调用存储器中的存储指令执行如上所述的一种动力学约束的低轨卫星GNSS实时精密定轨方法。
或者，包括可读存储介质，所述可读存储介质上存储有计算机程序，所述计算机程序执行时，实现如上所述的一种动力学约束的低轨卫星GNSS实时精密定轨方法。
与现有技术相比，本发明具有如下优点：
本发明将观测方程建立在无电离层组合观测值上，通过低轨卫星运动特性建立连续历元卫星坐标与瞬时加速度间关系，利用阶次数较低的地球重力场模型获得满足精度要求的瞬时加速度，无需估计大气阻力、太阳光压、地球反照压、经验力等非保守力参数，相比简化动力学定轨方法算力要求更低，同时定轨精度和稳定性优于未引入动力学信息的运动学方法。
附图说明
为了更清楚地说明本发明或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
图1是本发明实施例基于动力学约束的低轨卫星GNSS实时精密定轨方法的流程图。
图2是本发明实施例重力场模型阶数与低轨卫星瞬时加速度精度关系示意图。
图3是本发明实施例四颗卫星在不同约束强度下的轨道误差3DRMS示意图。
具体实施方式
为使本发明的目的、技术方案和优点更加清楚，下面将结合本发明中的附图和实施例，对本发明中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
实施例1
如图1所示，本发明实施例提供一种动力学约束的低轨卫星GNSS实时精密定轨方法，包括以下几个步骤：
步骤1，构建能够同时估计连续三个历元低轨卫星坐标的序贯最小二乘参数估计器，作为动力学约束的低轨卫星GNSS实时定轨参数估计器。
动力学约束的低轨卫星GNSS实时精密定轨首先需要构造适配的参数估计器。对于任意历元，低轨卫星GNSS观测方程可描述为：
（1）
式中，下标表示第历元；表示第历元的观测值，为消除电离层影响通常采用无电离层（Ionosphere-Free，IF）组合观测值；表示第历元的设计矩阵；表示第历元的待估参数矩阵；、、表示第历元低轨卫星坐标；、、表示GNSS卫星到低轨卫星向量关于、、的方向余弦；表示接收机钟差参数；...表示无电离层组合模糊度参数；表示观测噪声。
联立第、、历元的观测方程、、，统一表示为观测方程：
（2）
式中，、、表示第、、历元的GNSS卫星到低轨卫星的方向余弦，、、、表示第、、历元低轨卫星的坐标和接收机钟差，...表示无电离层组合模糊度，表示第、、历元观测噪声，矩阵化后表示为、和。
观测值减去参数初值的OMC（Observation Minus Calculation）矩阵可表示为：
（3）
式中，、、表示第、、历元低轨卫星的观测值，、、表示第、、历元低轨卫星初始坐标，本实施例通过伪距单点定位方式求解；、、表示第、、历元设计矩阵，基于低轨卫星初始坐标计算获得。
根据最小二乘准则形成法方程，求解待估参数向量：
（4）
式中，权矩阵由设定的观测方程先验精度确定，、表示由、、三个历元对应的参数组成的法方程矩阵和残差矩阵。
待估参数向量中包含第历元的低轨卫星坐标、、，从而实现低轨卫星的精密定轨。
为实现低轨卫星实时定轨，利用序贯最小二乘方法更新法方程，逐历元求解待估参数。流程为：①第历元，由式（4）形成第历元法方程，估计第历元低轨卫星坐标及其他参数；②第历元，将法方程中第历元对应低轨卫星坐标和接收机钟差参数通过参数消去方法消除，第和历元对应的法方程信息移至、历元对应位置，空出的历元对应位置填入历元信息，形成第历元法方程，估计第历元低轨卫星坐标及其他参数；③序贯执行步骤②，实现低轨卫星实时定轨。
步骤①、②具体过程推导如下：
①第历元法方程矩阵和残差矩阵写作：
（5）
（6）
其中：
（7）
（8）
（9）
（10）
（11）
（12）
式（7）-式（10）中，元素表示的方差，元素表示 的协方差，，；式（11）-式（12）中，元素表示的残差，。
式（5）为式（4）中法方程矩阵的展开形式，其中历元对应低轨卫星坐标、、及接收机钟差的法方程矩阵表示为，如式（7）所示，包含参数自身法方程系数项与参数间关系项；前序两个历元法方程矩阵、替换式（7）下标就能获得；不同历元参数间关系项，以式（8）为例，替换下标能够全部获得；无电离层组合模糊度参数对应法方程矩阵为，如式（9）所示；无电离层组合模糊度参数与其他参数间关系对应法方程矩阵为，如式（10）所示，下标替换能够获得全部项。
式（6）为式（4）中残差矩阵的展开形式，其中历元对应低轨卫星坐标、、及接收机钟差的残差矩阵表示为，如式（11）所示；前序两个历元残差矩阵、替换式（11）下标就能获得，无电离层组合模糊度参数对应残差矩阵为，如式（12）所示。
②第历元，将法方程矩阵和残差矩阵中第历元对应低轨卫星坐标、、和接收机钟差参数通过参数消去方法消除，然后将第和历元对应的法方程信息移至、历元对应位置，并填入第历元信息，得历元的法方程矩阵和残差矩阵：
（13）
（14）
其中，法方程矩阵中元素和残差矩阵中元素描述为：
（15）
（16）
式中，为法方程矩阵中元素，为残差矩阵中元素。
步骤2，利用低轨卫星运动特性，建立连续三个历元位置与瞬时加速度的关系，得到动力学约束。
基于牛顿运动方程，低轨卫星连续历元的平均加速度与相应历元卫星坐标可描述为：
（17）
（18）
式中，为低轨卫星第历元的平均加速度；为常系数矩阵；为历元间采样间隔；、、分别为、、历元卫星坐标，为便于推导统一使用惯性系（CRS）下坐标表示。
低轨卫星的平均加速度也可以写作多个历元瞬时加速度多项式组合的形式，表示为：
（19）
式中，为加权函数，为多项式系数矩阵，为低轨卫星瞬时加速度。
多项式阶数与选用的历元数有关，对于、、连续三个历元，使用三阶多项式，依据文献（A technique for modeling the Earth’s gravity field on thebasis of satellite accelerations）中方法将系数矩阵推导得到常系数矩阵，即：
（20）
与历元平均加速度相关的瞬时加速度表示为、、，进而式（19）可以写作：
（21）
由式（17）和式（21），、、连续三个历元平均加速度与瞬时加速度的关系可推导为：
（22）
为满足实时解算要求，将、、换元为、、，则式（22）可写为：
（23）
由此建立连续三个历元位置与瞬时加速度的关系，即为本发明所述动力学约束。
将动力学约束利用广义最小二乘原理写作虚拟观测值，计算公式为：
（24）
上述公式（17）-（24）均在惯性系（CRS）下推导，若最终待求低轨卫星坐标为地固系（TRS）下坐标，需使用对应历元地固系到惯性系（TRS to CRS）的旋转矩阵进行转换。的求解方法可参照公开文献（The IERS Conventions (2010): Reference systems andnew models）。将、、历元的旋转矩阵表示为、和，则：
（25）
为便于下文推导，设虚拟观测方程地固系（TRS）下设计矩阵为，即：
（26）
式中，为常系数矩阵。
步骤3，利用地球重力模型获取高精度瞬时加速度。
本实施例利用伪距单点定位的方式求得低轨卫星任意历元的先验初始坐标。对于任意历元低轨卫星先验初始坐标，首先将其转换至球坐标：
（27）
式中，是低轨卫星到地心的距离，是极角，是方位角。
在球坐标系中，地球重力场模型由球谐系数和表示，利用的模型阶次指取相应阶数的球谐系数进行引力势计算。引力势的球谐展开公式如下：
（28）
式中，是引力势，是低轨卫星到地心的距离，是极角，是方位角，是引力常数，是地球质量，是地球平均半径，是相关的Legendre函数，和表示其系数，，。
此时低轨卫星受到的引力为引力势的梯度，通过引力势对、、的单位向量、、求偏导获得：
（29）
式中，表示对引力势求偏导。
由于引力，为低轨卫星质量，由此能够计算得到任意历元低轨卫星瞬时加速度。
利用式（27）-（29）计算得到、、历元瞬时加速度值、和后，式（24）中仅剩坐标参数为未知参数，后续通过参数估计方法求解。
步骤4，将动力学约束带入步骤1构建的参数估计器，求解高精度低轨卫星坐标。
将步骤2、3中动力学约束虚拟观测方程形成法方程，与步骤1中GNSS观测方程形成的法方程叠加，最终参数估计结果即为施加动力学约束的低轨卫星GNSS实时定轨结果。具体推导过程如下。
①基于式（24）-（26）动力学约束形成的虚拟观测方程构建法方程。结合式（4）可以得到第历元法方程矩阵和残差矩阵：
（30）
式中，为第历元动力学约束方程形成的法方程矩阵，为第历元动力学约束方程形成的残差矩阵，为设计矩阵，由虚拟观测方程的精度决定，为观测值减去参数初值的OMC矩阵。
结合式（3），在此处可表示为：
（31）
式中，虚拟观测值设为0；表示地固系下坐标初值，本实施例通过伪距单点定位方式获得；表示瞬时加速度初值，通过式（27）-（29）计算获得，计算时所需低轨卫星坐标即为；历元间采样间隔和常系数矩阵已知，由此获得虚拟观测方程的OMC矩阵。
第历元，动力学约束形成的虚拟观测方程法方程矩阵和残差矩阵由（30）和（31）推导得到，表示为：
（32）
（33）
其中：
（34）
（35）
（36）
式中，表示历元对应低轨卫星坐标、、的法方程矩阵， 替换式（34）标获得，不同历元参数间关系项以式（35）为例，替换下标就能全部获得，表示历元对应低轨卫星坐标、、的残差矩阵，式（34）和式（35）中元素表示的方差，元素表示和的协方差，，，式（36）中元素表示的残差，，这里是为了区分动力学约束方程与GNSS观测方程差异，将各项方差、协方差、残差带有不同标注符号。
②由法方程叠加原理，将式（13）与式（32）、式（14）与式（33）叠加，最终获得本发明所述动力学约束解，如式（37）所示：
（37）
其中，最终待估参数由叠加后的法方程矩阵和叠加后的残差矩阵计算得到：
（38）
（39）
式中，表示GNSS观测方程对应法方程矩阵中的元素，表示第历元对应的法方程矩阵及其与其他元素间关系，表示动力学约束方程的法方程矩阵中的元素，表示GNSS观测方程的残差矩阵中的元素，表示第历元GNSS观测方程的残差矩阵中对应元素，表示动力学约束方程的残差矩阵中的元素。
为验证本发明所提方法的有效性，利用GRACE-C、SWARM-A、SWARM-B、SWARM-C四颗卫星2019年1月2日的星载观测数据进行测试。重力场模型不同阶数与低轨卫星瞬时加速度精度间关系如图2所示。从图2中可以看出，选用30阶和60阶次重力场模型可满足和动力学约束精度要求。利用本发明方法对四颗卫星进行测试，定轨残差如图3所示，图中、、分别表示未施加动力学约束、施加约束强度为和三种方案。求解过程中为满足动力学约束精度要求，和两种方案分别使用了30阶和60阶重力场模型。由图3可以看出，施加动力学约束后轨道定轨残差明显变小，整体精度和稳定性相比未约束解（）存在明显提升；同时，三种方案中动力学约束解轨道精度和稳定性最好。上述实验证明了本发明方法在实测低轨卫星任务中的有效性。
实施例2
基于同一发明构思，本发明还提供一种动力学约束的低轨卫星GNSS实时精密定轨系统，包括处理器和存储器，存储器用于存储程序指令，处理器用于调用存储器中的程序指令执行如上所述的一种动力学约束的低轨卫星GNSS实时精密定轨方法。
实施例3
基于同一发明构思，本发明还提供一种动力学约束的低轨卫星GNSS实时精密定轨系统，包括可读存储介质，所述可读存储介质上存储有计算机程序，所述计算机程序执行时，实现如上所述的一种动力学约束的低轨卫星GNSS实时精密定轨方法。
具体实施时，本发明技术方案提出的方法可由本领域技术人员采用计算机软件技术实现自动运行流程，实现方法的系统装置例如存储本发明技术方案相应计算机程序的计算机可读存储介质以及包括运行相应计算机程序的计算机设备，也应当在本发明的保护范围内。
本文中所描述的具体实施例仅仅是对本发明精神作举例说明。本发明所属技术领域的技术人员可以对所描述的具体实施案例，做各种各样的修改或补充或采用类似的方式替代，但并不会偏离本发明的精神或者超越所附权利要求书所定义的范围。

Claims (10)

1.一种动力学约束的低轨卫星GNSS实时精密定轨方法，其特征在于，包括以下步骤：
步骤1，构建能够同时估计连续三个历元低轨卫星坐标的序贯最小二乘参数估计器，作为动力学约束的低轨卫星GNSS实时定轨参数估计器；
对于任意历元，低轨卫星GNSS观测方程描述为：
（1）
式中，下标表示第历元，表示第历元的观测值，表示第历元的设计矩阵，表示第历元的待估参数矩阵，、、表示第历元低轨卫星坐标，、、表示GNSS卫星到低轨卫星向量关于、、的方向余弦，表示接收机钟差参数，...表示无电离层组合模糊度参数，n表示无电离层组合模糊度参数的个数，表示观测噪声；
联立三个历元低轨卫星GNSS观测方程，根据最小二乘准则形成法方程，利用序贯最小二乘方法更新法方程，逐历元求解待估参数，具体流程为：①第历元，形成第历元法方程，估计第历元低轨卫星坐标及其他参数；②第历元，将法方程中第历元对应低轨卫星坐标和接收机钟差参数通过参数消去方法消除，第和历元对应的法方程信息移至、历元对应位置，空出的历元对应位置填入历元信息，形成第历元法方程，估计第历元低轨卫星坐标及其他参数；③序贯执行步骤②，实现低轨卫星实时定轨；
步骤2，利用低轨卫星运动特性，建立连续三个历元位置与瞬时加速度的关系，得到动力学约束；
基于牛顿运动方程，建立低轨卫星连续历元的平均加速度与相应历元卫星坐标的关系，将低轨卫星的平均加速度表示为多个历元瞬时加速度多项式组合的形式，进而得到连续三个历元位置与瞬时加速度的关系，即动力学约束，然后利用最小二乘算法对连续三个历元位置与瞬时加速度的关系建立动力学约束虚拟观测方程；
步骤3，利用地球重力模型获取高精度瞬时加速度；
将历元低轨卫星先验初始坐标转换至球坐标，通过地球重力场模型获得引力势，对球坐标的单位向量求偏导后获得低轨卫星受到的引力，根据牛顿第二定理和低轨卫星质量求得任意历元低轨卫星瞬时加速度；
步骤4，将动力学约束带入步骤1构建的低轨卫星GNSS实时定轨参数估计器，求解高精度低轨卫星坐标；
联合步骤2中动力学约束虚拟观测方程和步骤3中瞬时加速度，形成当前历元的法方程，然后与步骤1中GNSS观测方程形成的法方程叠加，最终参数估计结果即为施加动力学约束的低轨卫星GNSS实时定轨结果。
2.如权利要求1所述的一种动力学约束的低轨卫星GNSS实时精密定轨方法，其特征在于：步骤1中联立第、、历元的观测方程、、，统一表示为观测方程：
（2）
式中，、、表示第、、历元的GNSS卫星到低轨卫星的方向余弦，、、、表示第、、历元低轨卫星的坐标和接收机钟差，...表示无电离层组合模糊度，表示第、、历元观测噪声，矩阵化后表示为、和；
观测值减去参数初值的OMC矩阵表示为：
（3）
式中，、、表示第、、历元低轨卫星的观测值，、、表示第、、历元低轨卫星初始坐标， 、、表示第、、历元设计矩阵，基于低轨卫星初始坐标计算获得；
对式（3）根据最小二乘准则形成法方程，求解待估参数向量：
（4）
式中，权矩阵由设定的观测方程先验精度确定，、表示由、、三个历元对应的参数组成的法方程矩阵和残差矩阵；
待估参数向量中包含第历元的低轨卫星坐标、、，从而实现低轨卫星的精密定轨。
3.如权利要求2所述的一种动力学约束的低轨卫星GNSS实时精密定轨方法，其特征在于：步骤1中利用序贯最小二乘方法更新法方程，逐历元求解待估参数，具体流程为：①第历元，由式（4）形成第历元法方程，估计第历元低轨卫星坐标及其他参数；②第历元，将法方程中第历元对应低轨卫星坐标和接收机钟差参数通过参数消去方法消除，第和历元对应的法方程信息移至、历元对应位置，空出的历元对应位置填入历元信息，形成第历元法方程，估计第历元低轨卫星坐标及其他参数；③序贯执行步骤②，实现低轨卫星实时定轨；
步骤①、②具体过程推导如下：
①第历元法方程矩阵和残差矩阵写作：
（5）
（6）
其中：
（7）
（8）
（9）
（10）
（11）
（12）
式（7）-式（10）中，元素表示的方差，元素表示 的协方差，， ；式（11）-式（12）中，元素表示的残差，；
式（5）为式（4）中法方程矩阵的展开形式，其中历元对应低轨卫星坐标、、及接收机钟差的法方程矩阵表示为，如式（7）所示，包含参数自身法方程系数项与参数间关系项；前序两个历元法方程矩阵、替换式（7）下标就能获得；不同历元参数间关系项，以式（8）为例，替换下标能够全部获得；无电离层组合模糊度参数对应法方程矩阵为，如式（9）所示；无电离层组合模糊度参数与其他参数间关系对应法方程矩阵为，如式（10）所示，下标替换能够获得全部项；式（6）为式（4）中残差矩阵的展开形式，其中历元对应低轨卫星坐标、、及接收机钟差的残差矩阵表示为，如式（11）所示；前序两个历元残差矩阵、替换式（11）下标就能获得，无电离层组合模糊度参数对应残差矩阵为，如式（12）所示；
②第历元，将法方程矩阵和残差矩阵中第历元对应低轨卫星坐标、、和接收机钟差参数通过参数消去方法消除，然后将第和历元对应的法方程信息移至、历元对应位置，并填入第历元信息，得历元的法方程矩阵和残差矩阵：
（13）
（14）
其中，法方程矩阵中元素和残差矩阵中元素描述为：
（15）
（16）
式中，为法方程矩阵中元素，为残差矩阵中元素。
4.如权利要求1所述的一种动力学约束的低轨卫星GNSS实时精密定轨方法，其特征在于：步骤2中基于牛顿运动方程，低轨卫星连续历元的平均加速度与相应历元卫星坐标描述为：
（17）
（18）
式中，为低轨卫星第历元的平均加速度；为常系数矩阵；为历元间采样间隔；、、分别为、、历元卫星坐标，统一使用惯性系下坐标表示；
将低轨卫星的平均加速度写作多个历元瞬时加速度多项式组合的形式，即：
（19）
式中，为加权函数，为多项式系数矩阵，为低轨卫星瞬时加速度；
多项式阶数与选用的历元数有关，对于、、连续三个历元，使用三阶多项式将系数矩阵推导得到常系数矩阵，即：
（20）
与历元平均加速度相关的瞬时加速度表示为、、，进而式（19）写作：
（21）
由式（17）和式（21），、、连续三个历元平均加速度与瞬时加速度的关系推导为：
（22）
为满足实时解算要求，将、、换元为、、，则式（22）写为：
（23）
由此建立连续三个历元位置与瞬时加速度的关系，即为步骤2所述动力学约束。
5.如权利要求4所述的一种动力学约束的低轨卫星GNSS实时精密定轨方法，其特征在于：步骤2中将动力学约束利用广义最小二乘原理写作虚拟观测值 ，计算公式为：
（24）
上述公式（17）-（24）均在惯性系下推导，若最终待求低轨卫星坐标为地固系下坐标，需使用对应历元地固系到惯性系的旋转矩阵进行转换，将 、、历元的旋转矩阵表示为、和，则：
（25）
设虚拟观测方程地固系下设计矩阵为，即：
（26）
式中，为常系数矩阵。
6.如权利要求1所述的一种动力学约束的低轨卫星GNSS实时精密定轨方法，其特征在于：步骤3中对于任意历元低轨卫星先验初始坐标，首先将其转换至球坐标：
（27）
式中，是低轨卫星到地心的距离，是极角，是方位角；
在球坐标系中，地球重力场模型由球谐系数和表示，利用的模型阶次指取相应阶数的球谐系数进行引力势计算，引力势的球谐展开公式如下：
（28）
式中，是引力势，是低轨卫星到地心的距离，是极角，是方位角，是引力常数，是地球质量，是地球平均半径，是相关的Legendre函数，和表示其系数，，；
此时低轨卫星受到的引力为引力势的梯度，通过引力势对、、的单位向量、、求偏导获得：
（29）
式中，表示对引力势求偏导；
由于引力，为低轨卫星质量，由此能够计算得到任意历元低轨卫星瞬时加速度；利用式（27）-（29）计算得到、、历元瞬时加速度值、和后，式（24）中仅剩坐标参数为未知参数，后续通过参数估计方法求解。
7.如权利要求3所述的一种动力学约束的低轨卫星GNSS实时精密定轨方法，其特征在于：步骤4中将步骤2、3中动力学约束虚拟观测方程形成法方程，与步骤1中GNSS观测方程形成的法方程叠加，最终参数估计结果即为施加动力学约束的低轨卫星GNSS实时定轨结果，具体推导过程如下：
基于式（24）-（26）动力学约束形成的虚拟观测方程构建法方程，结合式（4）得到第历元法方程矩阵和残差矩阵：
（30）
式中，为第历元动力学约束方程形成的法方程矩阵，为第历元动力学约束方程形成的残差矩阵，为设计矩阵，由虚拟观测方程的精度决定，为观测值减去参数初值的OMC矩阵；
结合式（3），在此处表示为：
（31）
式中，虚拟观测值设为0；表示地固系下坐标初值；表示瞬时加速度初值，通过式（27）-（29）计算获得，计算时所需低轨卫星坐标即为；历元间采样间隔和常系数矩阵已知，由此获得虚拟观测方程的OMC矩阵；
第历元，动力学约束形成的虚拟观测方程法方程矩阵和残差矩阵由（30）和（31）推导得到，表示为：
（32）
（33）
其中：
（34）
（35）
（36）
式中，表示历元对应低轨卫星坐标、、的法方程矩阵， 替换式（34）下标获得，不同历元参数间关系项以式（35）为例，替换下标就能全部获得，表示历元对应低轨卫星坐标、、的残差矩阵，式（34）和式（35）中元素表示的方差，元素表示和的协方差，，，式（36）中元素表示的残差，。
8.如权利要求7所述的一种动力学约束的低轨卫星GNSS实时精密定轨方法，其特征在于：步骤4中由法方程叠加原理，将式（13）与式（32）、式（14）与式（33）叠加，最终获得步骤4所述动力学约束解，如式（37）所示：
（37）
其中，最终待估参数由叠加后的法方程矩阵和叠加后的残差矩阵计算得到：
（38）
（39）
式中，表示GNSS观测方程对应法方程矩阵中的元素，，；表示第历元对应的法方程矩阵及其与其他元素间关系，，；表示动力学约束方程的法方程矩阵中的元素，，；表示GNSS观测方程的残差矩阵中的元素，；表示第历元GNSS观测方程的残差矩阵中对应元素，表示动力学约束方程的残差矩阵中的元素，。
9.一种动力学约束的低轨卫星GNSS实时精密定轨系统，其特征在于，包括处理器和存储器，存储器用于存储程序指令，处理器用于调用存储器中的程序指令执行如权利要求1-8任一项所述的一种动力学约束的低轨卫星GNSS实时精密定轨方法。
10.一种动力学约束的低轨卫星GNSS实时精密定轨系统，其特征在于，包括可读存储介质，所述可读存储介质上存储有计算机程序，所述计算机程序执行时，实现如权利要求1-8任一项所述的一种动力学约束的低轨卫星GNSS实时精密定轨方法。
CN202410552556.5A 2024-05-07 2024-05-07 一种动力学约束的低轨卫星gnss实时精密定轨方法及系统 Active CN118129767B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202410552556.5A CN118129767B (zh) 2024-05-07 2024-05-07 一种动力学约束的低轨卫星gnss实时精密定轨方法及系统

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202410552556.5A CN118129767B (zh) 2024-05-07 2024-05-07 一种动力学约束的低轨卫星gnss实时精密定轨方法及系统

Publications (2)

Publication Number Publication Date
CN118129767A CN118129767A (zh) 2024-06-04
CN118129767B true CN118129767B (zh) 2024-07-19

Family

ID=91248079

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202410552556.5A Active CN118129767B (zh) 2024-05-07 2024-05-07 一种动力学约束的低轨卫星gnss实时精密定轨方法及系统

Country Status (1)

Country Link
CN (1) CN118129767B (zh)

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108871348A (zh) * 2018-05-08 2018-11-23 中国人民解放军国防科技大学 一种利用天基可见光相机的低轨卫星自主定轨方法
CN112129300A (zh) * 2020-09-16 2020-12-25 武汉大学 位置间动力学约束的低轨卫星星载gnss精密定轨方法及系统

Family Cites Families (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
KR100809425B1 (ko) * 2006-09-29 2008-03-05 한국전자통신연구원 ＧＰＳ와 Ｇａｌｉｌｅｏ 데이터를 이용한 정밀 궤도결정시스템 및 그 방법
CN106338753B (zh) * 2016-09-22 2019-03-12 北京航空航天大学 一种基于地面站/星间链路/gnss联合测量的地球同步轨道星座定轨方法
CN110058287B (zh) * 2019-05-16 2022-03-15 北京合众思壮科技股份有限公司 一种低轨卫星定轨方法、装置及系统
CN111947667B (zh) * 2020-06-24 2022-08-12 火眼位置数智科技服务有限公司 一种基于运动学和动力学组合的低轨卫星实时高精度定轨方法

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108871348A (zh) * 2018-05-08 2018-11-23 中国人民解放军国防科技大学 一种利用天基可见光相机的低轨卫星自主定轨方法
CN112129300A (zh) * 2020-09-16 2020-12-25 武汉大学 位置间动力学约束的低轨卫星星载gnss精密定轨方法及系统

Also Published As

Publication number Publication date
CN118129767A (zh) 2024-06-04

Similar Documents

Publication Publication Date Title
CN100501331C (zh) 基于x射线脉冲星的导航卫星自主导航系统与方法
CN104075715B (zh) 一种结合地形和环境特征的水下导航定位方法
CN101216319B (zh) 基于联邦ukf算法的低轨卫星多传感器容错自主导航方法
CN107421550B (zh) 一种基于星间测距的地球-Lagrange联合星座自主定轨方法
Tscherning et al. A comparison of methods for computing gravimetric quantities from high degree spherical harmonic expansions
CN106871928A (zh) 基于李群滤波的捷联惯性导航初始对准方法
CN110514203B (zh) 一种基于isr-ukf的水下组合导航方法
CN102928858B (zh) 基于改进扩展卡尔曼滤波的gnss单点动态定位方法
CN112129300B (zh) 位置间动力学约束的低轨卫星星载gnss精密定轨方法及系统
CN106597507A (zh) Gnss/sins紧组合滤波平滑的高精度快速算法
CN108844536A (zh) 一种基于量测噪声协方差矩阵估计的地磁导航方法
CN109541647B (zh) 基于半天球格网点模型的gnss多路径效应修正方法
CN104048664A (zh) 一种导航卫星星座自主定轨的方法
CN112161632B (zh) 一种基于相对位置矢量测量的卫星编队初始定位方法
CN108508463B (zh) 基于Fourier-Hermite正交多项式扩展椭球集员滤波方法
CN116125503A (zh) 一种高精度卫星轨道确定及预报算法
CN110554443A (zh) 基于载波相位观测值和点加速度法确定地球重力场的方法
Xu et al. Orbit determination and thrust force modeling for a maneuvered GEO satellite using two-way adaptive Kalman filtering
CN118129767B (zh) 一种动力学约束的低轨卫星gnss实时精密定轨方法及系统
CN105527639A (zh) 一种基于平滑与外推的卫星定位方法
CN114740541B (zh) 基于主从星测速模式的小行星重力场反演方法及系统
Snay Introducing two spatial reference frames for regions of the Pacific Ocean
CN105652333A (zh) 一种低低跟踪重力测量卫星四点三线模型及其建立方法
CN112731504B (zh) 对月球探测器自主定轨的方法及装置
CN110579784B (zh) 基于卫星组合导航系统的卫星自主导航方法

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
GR01 Patent grant
GR01 Patent grant