CN111238489B - 一种低轨卫星大气阻力摄动建模和计算方法 - Google Patents

一种低轨卫星大气阻力摄动建模和计算方法 Download PDF

Info

Publication number
CN111238489B
CN111238489B CN202010202175.6A CN202010202175A CN111238489B CN 111238489 B CN111238489 B CN 111238489B CN 202010202175 A CN202010202175 A CN 202010202175A CN 111238489 B CN111238489 B CN 111238489B
Authority
CN
China
Prior art keywords
satellite
atmospheric resistance
perturbation
atmospheric
orbit
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
CN202010202175.6A
Other languages
English (en)
Other versions
CN111238489A (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.)
China Xian Satellite Control Center
Original Assignee
China Xian Satellite Control Center
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 China Xian Satellite Control Center filed Critical China Xian Satellite Control Center
Priority to CN202010202175.6A priority Critical patent/CN111238489B/zh
Publication of CN111238489A publication Critical patent/CN111238489A/zh
Application granted granted Critical
Publication of CN111238489B publication Critical patent/CN111238489B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

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
    • G05—CONTROLLING; REGULATING
    • G05B—CONTROL OR REGULATING SYSTEMS IN GENERAL; FUNCTIONAL ELEMENTS OF SUCH SYSTEMS; MONITORING OR TESTING ARRANGEMENTS FOR SUCH SYSTEMS OR ELEMENTS
    • G05B13/00—Adaptive control systems, i.e. systems automatically adjusting themselves to have a performance which is optimum according to some preassigned criterion
    • G05B13/02—Adaptive control systems, i.e. systems automatically adjusting themselves to have a performance which is optimum according to some preassigned criterion electric
    • G05B13/04—Adaptive control systems, i.e. systems automatically adjusting themselves to have a performance which is optimum according to some preassigned criterion electric involving the use of models or simulators
    • G05B13/042—Adaptive control systems, i.e. systems automatically adjusting themselves to have a performance which is optimum according to some preassigned criterion electric involving the use of models or simulators in which a parameter or coefficient is automatically adjusted to optimise the performance
    • G05B13/045—Adaptive control systems, i.e. systems automatically adjusting themselves to have a performance which is optimum according to some preassigned criterion electric involving the use of models or simulators in which a parameter or coefficient is automatically adjusted to optimise the performance using a perturbation signal
    • G—PHYSICS
    • G05—CONTROLLING; REGULATING
    • G05D—SYSTEMS FOR CONTROLLING OR REGULATING NON-ELECTRIC VARIABLES
    • G05D1/00—Control of position, course, altitude or attitude of land, water, air or space vehicles, e.g. using automatic pilots
    • G05D1/10—Simultaneous control of position or course in three dimensions

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Remote Sensing (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Automation & Control Theory (AREA)
  • General Physics & Mathematics (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Software Systems (AREA)
  • Medical Informatics (AREA)
  • Evolutionary Computation (AREA)
  • Astronomy & Astrophysics (AREA)
  • Health & Medical Sciences (AREA)
  • Artificial Intelligence (AREA)
  • Aviation & Aerospace Engineering (AREA)
  • Navigation (AREA)

Abstract

本发明公开了一种低轨卫星大气阻力摄动建模和计算方法，具体按照以下步骤实施：步骤1、建立卫星的Box‑Wing模型；步骤2、建立分段线性大气阻力摄动模型；步骤3、根据步骤1建立的Box‑Wing模型和步骤2建立的分段线性大气阻力摄动模型，计算t时刻折线形大气阻力摄动加速度。利用本发明的方法计算卫星的大气阻力摄动，能够克服传统的单一大气阻力因子难以精确描述大气变化引起的大气阻力摄动计算误差，有效提高定轨精度，为低轨卫星的精密定轨，特别是在磁暴等特殊条件下的精密定轨提供了技术支持。

Description

一种低轨卫星大气阻力摄动建模和计算方法
技术领域
本发明属于低轨卫星精密定轨技术领域，涉及一种低轨卫星大气阻力摄动建模和计算方法。
背景技术
卫星技术的广泛应用，对卫星轨道确定精度提出了越来越高的要求。卫星定轨精度主要取决于两方面的因素，一是测量数据的质量，二是采用的定轨策略和方法。在高精度测量系统中,卫星动力学模型精度往往难以满足精度需求，需要在精密定轨过程中增加估计大气阻力因子、太阳光压系数，或引入周期性的经验力，甚至采用简化动力学方法来补偿力模型的不足；另一方面，当跟踪设备误差较大时,又需要借助高精度的卫星动力学模型，在定轨过程中增加估计测量数据的设备系统差，对测量设备进行标校。因此，实际定轨中必须综合考虑测量精度和动力学模型精度。
低轨(Low Earth Orbit,LEO)卫星轨道运动受大气阻力摄动影响较大,因此在LEO精密定轨过程中,通常需要求解的参数包括卫星的初始位置和速度、大气阻力因子Cd，以及跟测站弧段相关的参数系统差。式(1)给出了传统大气阻力摄动加速度计算模型，主要参数包括大气阻力因子Cd，大气密度ρ、迎风面积A、卫星质量m、以及相对速度v。传统的定轨方法中，全弧段一般仅解算一个大气阻力因子Cd，用以补偿大气阻力摄动计算误差。这种方法仅适用于弧段适中且精度要求不高的情况，但对高精度LEO定轨任务由于目前的大气密度模型都不十分精确，特别是在太阳活动和地磁指数活动比较剧烈(即磁暴)的时候，大气阻力在整个弧段存在较大变化，用这种方法计算的大气阻力摄动精度较差，定轨精度受到很大限制。
Figure GDA0003871668710000021
以我国“XX四号”02星为例，该星位于200km高的近圆轨道上，大气阻力摄动是主要的定轨误差源。如果定轨弧段选取太长，则在积分运动方程时，大气阻力误差的累积效应将很明显；如果定轨弧段太短，则参加定轨的测量数据太少，同样难以精确定轨。因此，对该星的定轨弧长取为0.8～3天。“XX四号”02星采用GPS数据和USB数据定轨，并对测站USB外测设备进行标校。表1中给出了卫星在第90圈和91圈的两组实测数据算例，算例A采用1.5天数据，算例B采用0.8天数据，定轨期间发生了磁暴。
表1部分站USB测量系统差和随机差标校值(90和91圈算例)
Figure GDA0003871668710000022
表1结果显示，两组算例给出的活动站和渭南站USB测量系统差存在着明显差异，表明磁暴的发生给轨道摄动(主要是大气阻力摄动)的高精度建模和计算带来了困难，从而影响了测量系统的标校。因此，大气阻力摄动的高精度建模和计算(特别是在发生磁暴的情况下)，对于低轨卫星的精密定轨和一些科学应用，具有重要意义。
发明内容
本发明的目的是提供一种低轨卫星大气阻力摄动建模和计算方法，解决了现有技术中存在的单一大气阻力因子难以精确描述大气变化引起的大气阻力摄动计算误差的问题。
本发明所采用的技术方案是，一种低轨卫星大气阻力摄动建模和计算方法，具体按照以下步骤实施：
步骤1、建立卫星的Box-Wing模型；
步骤2、建立分段线性大气阻力摄动模型；
步骤3、根据步骤1建立的Box-Wing模型和步骤2建立的分段线性大气阻力摄动模型，计算t时刻折线形大气阻力摄动加速度。
本发明的特点还在于：
步骤1具体为，将卫星主体部分和太阳帆板部分简化为简单的几何模型盒和翼，即Box-Wing，卫星本体简化为长方体，太阳帆板简化为长方形的面板，根据地面标定的卫星轮廓的几何尺寸，卫星各表面的面积表示为：Ai(i＝1，…，nface)。
步骤2具体为，将整个定轨弧段[t0，tf]等分为n个弧段，则有n+1个节点，设各节点的大气阻力因子初值为Cd(ti)，(i＝1，2，…，n+1)，利用一次牛顿插值多项式，计算任意两个节点之间的折线型大气阻力因子Cd(t)：
Figure GDA0003871668710000031
步骤3中t时刻折线形大气阻力摄动加速度为：
Figure GDA0003871668710000032
其中，Ai为卫星第i个平面的面积；γi为卫星相对速度矢量v与第i个平面法向ni的夹角，ni由定义的卫星姿态确定；相对速度矢量
Figure GDA0003871668710000033
为卫星相对旋转大气的速度矢量，ω为地球的旋转角速度矢量，r、
Figure GDA0003871668710000034
分别为卫星位置、速度矢量，由卫星动力学积分得到。
步骤3还包括分段线性大气阻力因子估计。
分段线性大气阻力因子估计具体为：采用折线形大气阻力因子，实际的待估状态向量为：
Figure GDA0003871668710000041
卫星运动方程：
Figure GDA0003871668710000042
非线性观测方程：
Z＝h(Y,t) (5)
精密定轨的过程就是在给定初值Y0的条件下，利用最小二乘方法联合求解方程组(4)、(5)，得到新的状态估值
Figure GDA0003871668710000043
多次迭代直至收敛。
本发明的有益效果是：通过本发明一种低轨卫星大气阻力摄动建模和计算方法计算卫星的大气阻力摄动，能够克服传统的单一大气阻力因子难以精确描述大气变化引起的大气阻力摄动计算误差，有效提高定轨精度，为低轨卫星的精密定轨，特别是在磁暴等特殊条件下的精密定轨提供了技术支持。
附图说明
图1是本发明一种低轨卫星大气阻力摄动建模和计算方法的LEO卫星定轨计算流程图；
图2是本发明一种低轨卫星大气阻力摄动建模和计算方法的LEO卫星Box-Wing模型示意图。
图3是本发明一种低轨卫星大气阻力摄动建模和计算方法的折线形大气阻力因子示意图。
具体实施方式
下面结合附图和具体实施方式对本发明进行详细说明。
本发明一种低轨卫星大气阻力摄动建模和计算方法，具体按照以下步骤实施：
步骤1、建立卫星的Box-Wing模型；
低轨卫星大气阻力的计算与卫星的迎风面积相关，传统的卫星常数中一般给定一个参考迎风面积，然而，由于卫星在轨姿态存在变化，其真实迎风面积并不固定，基于卫星的结构分析，将卫星的主体部分和太阳帆板部分简化成简单的几何模型——盒和翼(“Box-Wing”)，如图2所示，即本体简化成长方体、太阳帆板简化为长方形的面板，因此，可以根据地面标定的卫星轮廓的几何尺寸，简单计算卫星各表面的面积Ai(i＝1，…，nface)
步骤2、建立分段线性大气阻力摄动模型；
将整个定轨弧段[t0，tf]等分为n个弧段，则有n+1个节点，给定各节点的大气阻力因子初值Cd(ti)，(i＝1，2，…，n+1)，利用一次牛顿插值多项式，计算任意两个节点之间的折线型大气阻力因子Cd(t)，如图3所示：
Figure GDA0003871668710000051
这样在整个弧段上计算的大气阻力因子呈分段线性，分段线性的大气阻力因子比直线型的大气阻力因子能够更加准确地反映大气密度快速波动引起的大气阻力变化；
步骤3、根据步骤1建立的Box-Wing模型和步骤2建立的分段线性大气阻力摄动模型，计算t时刻折线形大气阻力摄动加速度；
t时刻折线形大气阻力摄动加速度为：
Figure GDA0003871668710000052
其中，Ai为卫星第i个平面的面积；γi为卫星相对速度矢量v与第i个平面法向ni的夹角，ni由定义的卫星姿态确定；相对速度矢量
Figure GDA0003871668710000061
为卫星相对旋转大气的速度矢量，ω为地球的旋转角速度矢量，r、
Figure GDA0003871668710000062
分别为卫星位置、速度矢量，由卫星动力学积分得到；
卫星定轨过程是对状态参数迭代改进的过程，大气阻力摄动计算是其中的主要步骤之一，为了更加精确地计算大气阻力摄动加速度、估计卫星轨道，一般在定轨过程中将各节点的各节点的大气阻力因子Cd(ti)也作为待估参数进行估计，因此，精密定轨过程中的大气阻力摄动加速度计算过程实际上还包括了大气阻力因子Cd(ti)的估计过程，采用折线形大气阻力因子，实际的待估状态向量为：
Figure GDA0003871668710000063
卫星运动方程：
Figure GDA0003871668710000064
非线性观测方程：
Z＝h(Y，t) (5)
精密定轨的过程就是在给定初值Y0的条件下，利用最小二乘方法联合求解方程组(4)、(5)，得到新的状态估值
Figure GDA0003871668710000065
多次迭代直至收敛。
实施例
将本发明一种低轨卫星大气阻力摄动建模和计算方法具体与卫星的精密定轨过程相结合进行说明，如图1所示。
1、测量数据预处理
收集地面跟踪设备对低轨卫星的测量数据Z，对选取的定轨弧段上的测量数据进行预处理，剔除测量数据中的粗差和野值得到定轨观测值
Figure GDA0003871668710000066
2、设置初始概略轨道Y*、各节点的大气阻力因子Cd,i：
卫星的初始概略轨道由卫星的初定轨过程获得，一般精度在1km左右甚至更差，卫星的初始概略轨道用作精密定轨迭代过程的卫星轨道状态参数的初始值，利用各节点Cd,i的初值一般可设为2.0；
3、建立卫星轨道动力学模型
建立低轨卫星的轨道动力学模型，包括二体运动、地球非球形摄动、潮汐摄动、三体引力摄动、大气阻力摄动、太阳光压摄动、相对论效应摄动以及经验加速度摄动；其中大气阻力摄动按照本申请发明内容中方法建立，包括Box-Wing模型、分段线性阻力因子以及式(2)表示的各表面的阻力之和；
4、建立观测方程
根据观测数据类型建立相应的观测方程
Figure GDA0003871668710000071
5、逐一处理观测数据
(1)读入一条观测数据，记录观测值
Figure GDA0003871668710000072
和观测时刻ti；
(2)基于星轨道动力学模型和轨道初值积分到观测时刻ti，得到参考状态Y(ti)；
(3)由观测方程计算参考状态处的理论观测值h[Y(ti)]，并计算观测数据改进量
Figure GDA0003871668710000073
观测数据改进量又称为残差向量，包含了卫星轨道参数和待估参数的改进信息；
(4)计算相应的法方程系数
Figure GDA0003871668710000074
(5)重复(1)～(4)步，直到所有观测数据处理完毕，从而构成法方程：
Figure GDA0003871668710000081
6、轨道改进
利用最小二乘法对(7)式求解可解算出卫星轨道状态的改进量ΔY，则新的轨道状态和各节点的大气阻力因子Cd,i即为Y＝Y*+ΔY；
轨道改进的过程需要进行迭代，将Y作为新的参考轨道重复2～5步直到收敛。

Claims (1)

1.一种低轨卫星大气阻力摄动建模和计算方法，其特征在于，具体按照以下步骤实施：
步骤1、建立卫星的Box-Wing模型；具体为，将卫星主体部分和太阳帆板部分简化为简单的几何模型盒和翼，即Box-Wing，卫星本体简化为长方体，太阳帆板简化为长方形的面板，根据地面标定的卫星轮廓的几何尺寸，卫星各表面的面积表示为：Ai(i＝1，…，nface)；
步骤2、建立分段线性大气阻力摄动模型；具体为，将整个定轨弧段[t0，tf]等分为n个弧段，则有n+1个节点，设各节点的大气阻力因子初值为Cd(ti)，(i＝1，2，…，n+1)，利用一次牛顿插值多项式，计算任意两个节点之间的折线型大气阻力因子Cd(t)：
Figure FDA0003871668700000011
步骤3、根据步骤1建立的Box-Wing模型和步骤2建立的分段线性大气阻力摄动模型，计算t时刻折线形大气阻力摄动加速度；t时刻折线形大气阻力摄动加速度为：
Figure FDA0003871668700000012
其中，Ai为卫星第i个平面的面积；γi为卫星相对速度矢量v与第i个平面法向ni的夹角，ni由定义的卫星姿态确定；相对速度矢量
Figure FDA0003871668700000013
为卫星相对旋转大气的速度矢量，ω为地球的旋转角速度矢量，r、
Figure FDA0003871668700000014
分别为卫星位置、速度矢量，由卫星动力学积分得到；Cd(t)表示由(2)式计算的t时刻的大气阻力因子；ρ表示大气密度；m表示卫星质量；nface表示卫星面板个数；
步骤3还包括分段线性大气阻力因子估计，具体为：采用折线形大气阻力因子，实际的待估状态向量为：
Figure FDA0003871668700000021
r、
Figure FDA0003871668700000022
分别为卫星位置、速度矢量，上标T表示转置，Cd(t0)…Cd(tn+1)分别表示节点t0…tn+1待估计的大气阻力因子；
卫星运动方程：
Figure FDA0003871668700000023
非线性观测方程：
Z＝h(Y，t) (5)
精密定轨的过程就是在给定初始值
Figure FDA0003871668700000024
的条件下，利用最小二乘方法联合求解公式(4)和公式(5)，得到新的状态估计值
Figure FDA0003871668700000025
求解过程采用迭代算法，即将本次迭代估计值
Figure FDA0003871668700000026
作为新的初始值Y0代入公式(4)和公式(5)重新进行最小二乘估计，多次迭代直至收敛后，得到最终的状态估计
Figure FDA0003871668700000027
CN202010202175.6A 2020-03-20 2020-03-20 一种低轨卫星大气阻力摄动建模和计算方法 Active CN111238489B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010202175.6A CN111238489B (zh) 2020-03-20 2020-03-20 一种低轨卫星大气阻力摄动建模和计算方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010202175.6A CN111238489B (zh) 2020-03-20 2020-03-20 一种低轨卫星大气阻力摄动建模和计算方法

Publications (2)

Publication Number Publication Date
CN111238489A CN111238489A (zh) 2020-06-05
CN111238489B true CN111238489B (zh) 2022-11-29

Family

ID=70866073

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010202175.6A Active CN111238489B (zh) 2020-03-20 2020-03-20 一种低轨卫星大气阻力摄动建模和计算方法

Country Status (1)

Country Link
CN (1) CN111238489B (zh)

Families Citing this family (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112298614B (zh) * 2020-09-18 2021-09-14 中国人民解放军战略支援部队航天工程大学 一种推力在轨标定试验方法
CN112591150B (zh) * 2021-01-05 2022-09-13 成都天巡微小卫星科技有限责任公司 一种控制超低轨道卫星姿态的大气阻力矩补偿方法及系统
CN113591265B (zh) * 2021-05-28 2024-04-30 中国人民解放军63920部队 一种计算火星探测器大气阻力的方法
CN113218660B (zh) * 2021-06-14 2022-09-13 中国西安卫星测控中心 一种电推力矢量在轨标定方法
CN118413261A (zh) * 2024-03-20 2024-07-30 株洲太空星际卫星科技有限公司 卫星测控数传弧段自动计算方法、装置、设备及介质

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US5687084A (en) * 1992-05-26 1997-11-11 Microcosm, Inc. Satellite orbit maintenance system
CN109323698A (zh) * 2018-12-03 2019-02-12 西安四方星途测控技术有限公司 一种空间目标陨落多模型跟踪引导技术

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US5687084A (en) * 1992-05-26 1997-11-11 Microcosm, Inc. Satellite orbit maintenance system
CN109323698A (zh) * 2018-12-03 2019-02-12 西安四方星途测控技术有限公司 一种空间目标陨落多模型跟踪引导技术

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
HY-2卫星DORIS厘米级精密定轨;朱俊等;《宇航学报》;20130228;第34卷(第2期);第163-169页 *

Also Published As

Publication number Publication date
CN111238489A (zh) 2020-06-05

Similar Documents

Publication Publication Date Title
CN111238489B (zh) 一种低轨卫星大气阻力摄动建模和计算方法
CN110837094B (zh) 基于低轨卫星的无奇点20轨道根数拟合方法
CN110132261A (zh) 一种基于数值拟合的高精度星上轨道预报方法
CN109323698B (zh) 一种空间目标陨落多模型跟踪引导方法
CN111914411B (zh) 一种全姿态四轴转台框架角指令解算方法
CN111220347B (zh) 一种飞行器气动协调修正方法
CN101949703A (zh) 一种捷联惯性/卫星组合导航滤波方法
CN108959734B (zh) 一种基于实时递推太阳光压力矩辨识方法及系统
CN111766397B (zh) 一种基于惯性/卫星/大气组合的气象风测量方法
CN106840203A (zh) 惯导/气压高度计/gps组合导航系统中气压高度计校正方法
CN113074703A (zh) 多源卫星高度计的数据处理方法、装置、设备及存储介质
CN107101649A (zh) 一种空间飞行器制导工具在轨误差分离方法
CN110703355A (zh) 一种星载加速度计的校准方法及装置
CN108875174A (zh) 一种基于多段打靶法的不变拟周期轨道确定方法
CN116125503A (zh) 一种高精度卫星轨道确定及预报算法
CN110632636A (zh) 一种基于Elman神经网络的载体姿态估计方法
Beuselinck et al. Determination of attitude motion of the Foton M-3 satellite according to the data of onboard measurements of the Earth’s magnetic field
CN112393835B (zh) 一种基于扩展卡尔曼滤波的小卫星在轨推力标定方法
CN111308127B (zh) 一种基于大气物理机制的星载加速度计校准方法
CN111605736B (zh) 地月l2点转移轨道最优误差修正点选择方法
CN112327333A (zh) 一种卫星位置计算方法及装置
CN114771873A (zh) 一种超低轨卫星轨道自主精确维持方法
CN110006455A (zh) 用于冗余惯导系统中加速度计误差参数的快速标定方法
CN109100705A (zh) 星载激光测高仪在轨标定模型中权矩阵的确定方法
CN111737814B (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
CP02 Change in the address of a patent holder

Address after: No.462, East Xianning Road, Xi'an, Shaanxi 710043

Patentee after: CHINA XI'AN SATELLITE CONTROL CENTER

Address before: No.462, East Xianning Road, Xincheng District, Xi'an, Shaanxi 710043

Patentee before: CHINA XI'AN SATELLITE CONTROL CENTER

CP02 Change in the address of a patent holder