CN112393835B - 一种基于扩展卡尔曼滤波的小卫星在轨推力标定方法 - Google Patents
一种基于扩展卡尔曼滤波的小卫星在轨推力标定方法 Download PDFInfo
- Publication number
- CN112393835B CN112393835B CN202011212314.XA CN202011212314A CN112393835B CN 112393835 B CN112393835 B CN 112393835B CN 202011212314 A CN202011212314 A CN 202011212314A CN 112393835 B CN112393835 B CN 112393835B
- Authority
- CN
- China
- Prior art keywords
- orbit
- time
- satellite
- earth
- thrust
- 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
Links
Images
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01L—MEASURING FORCE, STRESS, TORQUE, WORK, MECHANICAL POWER, MECHANICAL EFFICIENCY, OR FLUID PRESSURE
- G01L5/00—Apparatus for, or methods of, measuring force, work, mechanical power, or torque, specially adapted for specific purposes
- G01L5/13—Apparatus for, or methods of, measuring force, work, mechanical power, or torque, specially adapted for specific purposes for measuring the tractive or propulsive power of vehicles
- G01L5/133—Apparatus for, or methods of, measuring force, work, mechanical power, or torque, specially adapted for specific purposes for measuring the tractive or propulsive power of vehicles for measuring thrust of propulsive devices, e.g. of propellers
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01L—MEASURING FORCE, STRESS, TORQUE, WORK, MECHANICAL POWER, MECHANICAL EFFICIENCY, OR FLUID PRESSURE
- G01L25/00—Testing or calibrating of apparatus for measuring force, torque, work, mechanical power, or mechanical efficiency
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01M—TESTING STATIC OR DYNAMIC BALANCE OF MACHINES OR STRUCTURES; TESTING OF STRUCTURES OR APPARATUS, NOT OTHERWISE PROVIDED FOR
- G01M15/00—Testing of engines
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/10—Geometric CAD
- G06F30/15—Vehicle, aircraft or watercraft design
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
- G06F2119/14—Force analysis or force optimisation, e.g. static or dynamic forces
Landscapes
- Physics & Mathematics (AREA)
- Engineering & Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Geometry (AREA)
- Combustion & Propulsion (AREA)
- Chemical & Material Sciences (AREA)
- Theoretical Computer Science (AREA)
- Mathematical Optimization (AREA)
- Mathematical Analysis (AREA)
- Pure & Applied Mathematics (AREA)
- Computer Hardware Design (AREA)
- Evolutionary Computation (AREA)
- General Engineering & Computer Science (AREA)
- Computational Mathematics (AREA)
- Aviation & Aerospace Engineering (AREA)
- Automation & Control Theory (AREA)
- Ocean & Marine Engineering (AREA)
- Navigation (AREA)
Abstract
本发明公开了一种基于扩展卡尔曼滤波的小卫星在轨推力标定方法,根据小卫星在轨运行情况,建立基于高斯变分方程的小卫星动力学模型,小卫星动力学模型中考虑包括J2摄动、太阳光压和第三体引力的多种摄动力,基于小卫星动力学模型,进行扩展卡尔曼滤波的算法设计,在滤波器内对状态预测进行更新。对小推力大小和其本体坐标系下的方向进行在轨标定。本发明可以在轨标定多种摄动力下小卫星的推力,保证了小推力机动的有效性。
Description
技术领域
本发明属于航天器轨道运动技术领域,具体涉及一种基于扩展卡尔曼滤波的小卫星在轨推力标定方法。
背景技术
近几年,随着空间资源的不断减少,空间卫星越来越趋向于小型化,小卫星具有体积小、质量轻、反应快速、性价比高等优点,随着小卫星任务数量的增加和小卫星技术的不断提高,小卫星的发展已进入到一个高速阶段。
由于小卫星本身质量小,能够携带的燃料较少,因此为了实现小卫星在轨的快速机动能力,需要采用小推力的推进形式。小推力既可用于小卫星的姿态控制、轨道调整和机动等,也能进行小卫星的编队任务,因此小推力机动在小卫星任务上的应用前景很广阔,其研究工作越来越受到重视。
传统的推力标定方法,一般在地面对推力发动机进行推力测试,通过对推力参数的测量,以得到的实验数据对发动机进行改进。但是,由于太空环境温度、压力等的特殊性,地面测量技术手段有限,很难实现小卫星推力的高精度标定。为了保证小推力机动的有效性,需要对其进行在轨标定。现存的在轨标定方法,仅仅考虑了存在推力的情况,考虑的摄动有限,造成所标定的推力精度有限。
发明内容
为解决现有技术中存在的问题,本发明的目的是提供一种基于扩展卡尔曼滤波的小卫星在轨推力标定方法,使用该方法能够将小卫星在轨运动的动力学模型进行状态扩维,利用扩展卡尔曼滤波器对卫星轨道信息进行预测更新,对小推力大小和其本体坐标系下的方向进行在轨标定的仿真验证,实现小卫星的小推力在轨标定。
为了实现上述目的,本发明采用的技术方案如下:
一种基于扩展卡尔曼滤波的小卫星在轨推力标定方法,包括如下过程:
建立地球惯性坐标系和卫星轨道坐标系,小卫星在以地球为中心的椭圆轨道上运行;
结合地球惯性坐标系、卫星轨道坐标系、J2摄动、太阳光压、第三体引力和轨道六要素建立基于高斯变分方程的小卫星动力学模型;
基于所述小卫星动力学模型,进行扩展卡尔曼滤波的算法设计,得到扩展卡尔曼滤波器,在扩展卡尔曼滤波器内对小卫星的状态预测进行更新,对小推力大小和小卫星本体坐标系下的方向进行在轨标定的仿真验证,实现小卫星的小推力在轨标定。
优选的:地球惯性坐标系EXYZ中,坐标系原点在地心E,X轴指向春分点,Z轴垂直于地球赤道平面向上,Y轴满足右手法则;
卫星轨道坐标系Txyz中,坐标系原点在卫星质心T,z轴沿卫星矢径方向,指向地心,y轴垂直于轨道平面,x轴沿速度方向,遵循右手法则。
优选的,小卫星动力学模型如下:
其中,FT+FJ2+FSRP+FM=[Fx,Fy,Fz]T,分别FT为小推力,FJ2为J2摄动力,FSRP为太阳光压摄动力,FM为第三体摄动力,为表示方便,上述力在下文中均由该力的加速度表示,即aT为小推力,aJ2为J2摄动力,aSRP为太阳光压摄动力,aM为第三体摄动力;a为轨道半长轴,e为偏心率,i为轨道倾角,f为真近点角,Ω为升交点赤经,ω为近地点幅角,h为角动量,p为半通径,r为位置矢量的大小,t为时间。
式中,表示地球惯性系下太阳到小卫星的位置矢量,AU表示地球到太阳的距离,P表示单位AU距离下的太阳光压力,ν表示地球的阴影函数,Cr表示辐射压力系数,m表示小卫星的质量,A表示小卫星的横截面积,上标ECI是指该矢量在惯性系下表示,如太阳光压摄动力aSRP在惯性系下表示为下文中的上标ECI均为此意。
式中,rM,r分别表示地球惯性系下第三体的位置矢量和卫星的位置矢量,其中,所述第三体包括太阳或者月球,AU表示地球到太阳的距离,G为万有引力常量,M表示太阳或者月球的质量,m表示卫星的质量,A表示卫星的横截面积。
优选的,所述扩展卡尔曼滤波器如下:
其中,为k时刻状态最优估值,为k+1时刻状态最优估值,是由用状态模型推算获得的k+1时刻的预测状态量,是离散化的系统方程,目的是得到k+1时刻的状态预测值,既是时间t的非线性函数,又是k时刻状态最优估值的非线性函数;Γk为k时刻的系统干扰矩阵,W为系统干扰量,Kk+1为k+1时刻的增益矩阵,Zk+1为k+1时刻的量测值,h为角动量,Pk+1|k为k时刻转移至k+1时刻的预测误差方差矩阵,Hk+1为k+1时刻的量测矩阵,Rk+1为k+1时刻的量测噪声方差矩阵,Φk+1,k为系统从k时刻转移至k+1时刻的状态转移矩阵,Pk为k时刻的滤波(估计)误差方差矩阵,Qk为k时刻系统干扰误差矩阵,I为单位矩阵,上标T是指该矩阵的转置,即为Hk+1的转置,为Φk+1,k的转置,为Γk的转置。
与现有技术相比,本发明具有以下有益效果:
本发明根据小卫星在轨运行的特点,考虑了J2摄动、太阳光压和第三体引力,建立了基于高斯变分方程的小卫星动力学模型,小卫星动力学模型含有与卫星轨道相关的六要素,更容易描述轨道的特征,同时六要素也容易建立与时间的联系,能够刻画出某时间点的轨道要素,进而可以表示出推力大小。利用扩展卡尔曼滤波器对状态信息进行预测,考虑了轨道运动过程中的噪声,更贴近实际工程问题,利用仿真标定出的推力大小,限制了误差范围,更加精确。综上所述,本发明可以在轨标定多种摄动力下小卫星的推力,保证了小推力机动的有效性。
附图说明
图1为本发明提出的建立地球惯性坐标系EXYZ示意图;
图2为本发明提出的卫星轨道坐标系Txyz示意图,用于描述小卫星在轨运动情况;
图3为本发明提出的卫星小推力的示意图,用于描述推力的方向;
图4为本发明提出的小推力在轨标定流程示意图,用来描述整个小卫星运行过程中在轨标定推力的流程。
具体实施方式
下面通过附图和实施例,对本发明的技术方案做进一步的详细描述。
本发明基于扩展卡尔曼滤波的小卫星在轨推力标定方法,包括以下步骤:
步骤一、建立地球惯性坐标系EXYZ和卫星轨道坐标系Txyz,小卫星在以地球为中心的椭圆轨道上运行;
步骤二、建立基于高斯变分方程的小卫星动力学模型,小卫星动力学模型中考虑包括J2摄动、太阳光压、第三体引力的多种摄动力以及和轨道六要素;
步骤三、基于小卫星动力学模型,进行扩展卡尔曼滤波的算法设计,在滤波器内对状态预测进行更新;
步骤四、对小推力大小和其本体坐标系下的方向进行在轨标定,满足既定标定误差;具体的,对小推力大小和其本体坐标系下的方向在3小时内进行在轨标定,要求标定误差小于5%。
参见图1,建立的地球惯性坐标系EXYZ具体细节如下:
地球惯性坐标系EXYZ,坐标系原点在地心E,X轴指向春分点,Z轴垂直于地球赤道平面向上,Y轴满足右手法则。
参见图2,建立的卫星轨道坐标系Txyz具体细节如下:
卫星轨道坐标系Txyz,坐标系原点在卫星质心T,z轴沿卫星矢径方向,指向地心,y轴垂直于轨道平面,x轴沿速度方向,遵循右手法则;
参见图3,小卫星在近地球轨道上运行,机动方式为小推力,同时在轨运行过程中,小卫星也会受到J2摄动、太阳光压摄动、第三体引力摄动,小卫星的星上定轨方式基于全球导航卫星系统(GNSS)的信号,全球导航卫星系统输出轨道信息,并存在一定大小的误差。
参见图4,建立如下的基于高斯变分方程的动力学模型:
常用于确定轨道的六个轨道根数,称之为轨道六要素,分别为轨道半长轴a、偏心率e、轨道倾角i、真近点角f、升交点赤经Ω和近地点幅角ω。由于存在摄动力和控制力,高斯变分方程来描述轨道要素随时间的变化,具体的,小卫星动力学模型如下:
其中,h为角动量,p半通径,r为位置矢量的大小,t为时间,摄动力和控制力的统一形式
F=FT+FJ2+FSRP+FM=[Fx,Fy,Fz]T
其中FT、FJ2、FSRP、FM分别为小推力、J2摄动力、太阳光压摄动力和第三体摄动力。
在地球惯性系下的太阳光压和太阳月球等第三体引力引起的摄动力为
式中,rM,r分别表示地球惯性系下太阳到卫星的位置矢量、第三体(太阳或者月球)的位置矢量和卫星的位置矢量,AU表示地球到太阳的距离,P表示单位AU距离下的太阳光压力,ν表示地球的阴影函数,G为万有引力常量,M表示太阳或者月球的质量,Cr表示辐射压力系数,m表示卫星的质量,A表示卫星的横截面积。
扩展卡尔曼滤波器算法模型建立过程如下:
针对非线性系统的扩展卡尔曼滤波器(EKF)的基本思想是将非线性系统展开成泰勒级数的形式,得到非线性系统的线性化模型,再利用递推方程进行系统的状态估计。
其中,Xk+1为k+1时刻的状态量,f(Xk,k)为系统状态方程,目的是由k时刻的状态量得到k+1时刻的状态量,Γk为k时刻的系统干扰矩阵,Wk为k时刻的系统干扰量,Zk+1为k+1时刻的量测值,h(Xk+1,k+1)为量测方程,目的是得到k+1时刻的量测,Vk+1为k+1时刻的量测噪声。
可得非线性系统的最优估计状态方程和量测方程为
令非线性系统标称状态和最优估值之间存在偏差
将状态方程和量测方程在估值附近进行一阶泰勒级数展开,写成标准差分方程的形式
然后推导出误差状态的卡尔曼滤波方程为
式中,为k时刻误差状态量,δXk+1为k+1时刻误差状态量,是由误差状态用状态模型推算获得的k+1时刻的预测误差状态量,δZk+1为k+1时刻误差量测值,Φk+1,k为系统从k时刻转移至k+1时刻的状态转移矩阵,Hk+1为k+1时刻的量测矩阵,Pk+1|k为k时刻转移至k+1时刻的预测误差方差矩阵,Qk为k时刻系统干扰误差矩阵,I为单位矩阵。
由于初始时刻的状态偏差最优估计和预测值均为零,因此对应的扩展卡尔曼滤波方程为
基于上述方程在扩展卡尔曼滤波器内对状态预测进行更新。
利用EKF对小推力的大小及方向的在轨标定的过程如下:
首先,将笛卡尔动力学模型进行状态扩维,即包含带估计的参数:
式中,为卫星状态的二阶导,代表其加速度,μ为地心引力常数,r为卫星位置矢量,r为卫星位置矢量的大小,aT为作用在卫星上的小推力,aJ2为J2摄动力,aSRP为太阳光压摄动力,aM为第三体摄动力,为卫星质量的变化,为小推力的变化,参见图3,和θ为定义小推力方向的角,和分别为该角度的变化。
对一个测量周期内的卫星扩张状态进行预测,写出预测方程;然后根据EKF基本原理中的上述扩展卡尔曼滤波方程分别写出状态估值方程、滤波增益方程、一步预测均方误差方程和估计均方误差方程。
根据扩展卡尔曼滤波器的算法,对小推力大小和其本体坐标系下的方向在3小时内进行在轨标定的仿真验证,要求标定误差小于5%。
以上内容仅为说明本发明的技术思想,不能以此限定本发明的保护范围,凡是按照本发明提出的技术思想,在技术方案基础上所做的任何改动,均落入本发明权利要求书的保护范围之内。
Claims (2)
1.一种基于扩展卡尔曼滤波的小卫星在轨推力标定方法,其特征在于,包括如下过程:
建立地球惯性坐标系和卫星轨道坐标系,小卫星在以地球为中心的椭圆轨道上运行;
结合地球惯性坐标系、卫星轨道坐标系、J2摄动、太阳光压、第三体引力和轨道六要素建立基于高斯变分方程的小卫星动力学模型;
基于所述小卫星动力学模型,进行扩展卡尔曼滤波的算法设计,得到扩展卡尔曼滤波器,在扩展卡尔曼滤波器内对小卫星的状态预测进行更新,对小推力大小和小卫星本体坐标系下的方向进行在轨标定的仿真验证,实现小卫星的小推力在轨标定;
地球惯性坐标系EXYZ中,坐标系原点在地心E,X轴指向春分点,Z轴垂直于地球赤道平面向上,Y轴满足右手法则;
卫星轨道坐标系Txyz中,坐标系原点在卫星质心T,z轴沿卫星矢径方向,指向地心,y轴垂直于轨道平面,x轴沿速度方向,遵循右手法则;
小卫星动力学模型如下:
其中,FT+FJ2+FSRP+FM=[Fx,Fy,Fz]T,分别FT为小推力,FJ2为J2摄动力,FSRP为太阳光压摄动力,FM为第三体摄动力;a为轨道半长轴,e为偏心率,i为轨道倾角,f为真近点角,Ω为升交点赤经,ω为近地点幅角,h为角动量,p为半通径,r为位置矢量的大小,t为时间;
式中,rM,r分别表示地球惯性系下第三体的位置矢量和卫星的位置矢量,其中,所述第三体包括太阳或者月球,AU表示地球到太阳的距离,G为万有引力常量,M表示太阳或者月球的质量,m表示卫星的质量,A表示卫星的横截面积。
2.根据权利要求1所述的一种基于扩展卡尔曼滤波的小卫星在轨推力标定方法,其特征在于,所述扩展卡尔曼滤波器如下:
其中,为k时刻状态最优估值,为k+1时刻状态最优估值,是由用状态模型推算获得的k+1时刻的预测状态量,是离散化的系统方程,目的是得到k+1时刻的状态预测值,既是时间t的非线性函数,又是k时刻状态最优估值的非线性函数;Γk为k时刻的系统干扰矩阵,W为系统干扰量,Kk+1为k+1时刻的增益矩阵,Zk+1为k+1时刻的量测值,h为角动量,Pk+1|k为k时刻转移至k+1时刻的预测误差方差矩阵,Hk+1为k+1时刻的量测矩阵,Rk+1为k+1时刻的量测噪声方差矩阵,Φk+1,k为系统从k时刻转移至k+1时刻的状态转移矩阵,Pk为k时刻的滤波(估计)误差方差矩阵,Qk为k时刻系统干扰误差矩阵,I为单位矩阵。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202011212314.XA CN112393835B (zh) | 2020-11-03 | 2020-11-03 | 一种基于扩展卡尔曼滤波的小卫星在轨推力标定方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202011212314.XA CN112393835B (zh) | 2020-11-03 | 2020-11-03 | 一种基于扩展卡尔曼滤波的小卫星在轨推力标定方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN112393835A CN112393835A (zh) | 2021-02-23 |
CN112393835B true CN112393835B (zh) | 2021-12-21 |
Family
ID=74598026
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN202011212314.XA Active CN112393835B (zh) | 2020-11-03 | 2020-11-03 | 一种基于扩展卡尔曼滤波的小卫星在轨推力标定方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN112393835B (zh) |
Families Citing this family (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN113602534B (zh) * | 2021-06-26 | 2023-02-28 | 山东航天电子技术研究所 | 一种微型电推进推力大小的在轨标定方法 |
Family Cites Families (6)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
GB2388578A (en) * | 2002-05-18 | 2003-11-19 | Astrium Ltd | Satellite launch vehicles and a method of launching satellites into orbit |
CN101381004B (zh) * | 2008-08-20 | 2010-11-10 | 南京航空航天大学 | 基于大气阻力的微小卫星编队飞行控制方法及控制装置 |
CN103940431B (zh) * | 2014-04-11 | 2016-08-10 | 北京空间飞行器总体设计部 | 基于gnss精密定轨的圆轨道切向小推力在轨标定方法 |
CN108061660B (zh) * | 2017-10-23 | 2019-09-17 | 上海卫星工程研究所 | 基于线振动测量的卫星发动机在轨推力实时标定方法 |
CN108454886B (zh) * | 2018-01-09 | 2019-11-12 | 北京控制工程研究所 | 一种电推进系统毫牛级推力在轨标定方法 |
CN109190155B (zh) * | 2018-07-25 | 2022-09-16 | 西北工业大学 | 一种采用电推进/太阳帆推进的混合连续小推力轨道设计方法 |
-
2020
- 2020-11-03 CN CN202011212314.XA patent/CN112393835B/zh active Active
Also Published As
Publication number | Publication date |
---|---|
CN112393835A (zh) | 2021-02-23 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN110487301B (zh) | 一种雷达辅助机载捷联惯性导航系统初始对准方法 | |
CN110794863B (zh) | 一种控制性能指标可定制的重型运载火箭姿态控制方法 | |
Olds et al. | IRVE-3 post-flight reconstruction | |
Yuan | JPL level-2 processing standards document for level-2 product release 06 | |
CN108548542B (zh) | 一种基于大气阻力加速度测量的近地轨道确定方法 | |
CN104792340A (zh) | 一种星敏感器安装误差矩阵与导航系统星地联合标定与校正的方法 | |
CN111552003B (zh) | 基于球卫星编队的小行星引力场全自主测量系统及方法 | |
CN110316402A (zh) | 一种编队控制模式下的卫星姿态控制方法 | |
Roh et al. | Dynamic accuracy improvement of a MEMS AHRS for small UAVs | |
CN112393835B (zh) | 一种基于扩展卡尔曼滤波的小卫星在轨推力标定方法 | |
Srivastava et al. | Attitude determination and control system for a leo debris chaser small satellite | |
CN112817233B (zh) | 一种基于迭代学习控制的小天体探测器绕飞段轨道跟踪控制方法 | |
CN105068425A (zh) | 一种适用于敏捷卫星姿态确定的状态反馈鲁棒非脆弱控制方法 | |
CN108871312B (zh) | 一种重力梯度仪及星敏感器的联合定姿方法 | |
Hong et al. | Application of EKF for missile attitude estimation based on “SINS/CNS” integrated guidance system | |
Prieto et al. | A drag free control based on model predictive techniques | |
CN111649738B (zh) | 微重力场下的加速度计初始姿态解算方法 | |
Xiong et al. | Online calibration research on the lever arm effect for the hypersonic vehicle | |
CN112498747A (zh) | 深空探测器推力矢量与目标加速度偏差角计算方法及系统 | |
Xie et al. | Autonomous guidance, navigation, and control of spacecraft | |
Kinzie et al. | Dual quaternion-based dynamics and control for gravity recovery missions | |
Zhang et al. | A Systematic Approach for Calibrating the Inertial Sensor of Taiji-1 | |
CN113282097B (zh) | 一种geo博弈航天器相对位置非球形摄动误差的控制方法 | |
CN116300417B (zh) | 一种大型分布式空间望远镜主次镜编队控制方法 | |
CN116817932B (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 |