CN109063380B - 一种静止轨道电推进卫星故障检测方法及位置保持方法 - Google Patents

一种静止轨道电推进卫星故障检测方法及位置保持方法 Download PDF

Info

Publication number
CN109063380B
CN109063380B CN201811059085.5A CN201811059085A CN109063380B CN 109063380 B CN109063380 B CN 109063380B CN 201811059085 A CN201811059085 A CN 201811059085A CN 109063380 B CN109063380 B CN 109063380B
Authority
CN
China
Prior art keywords
orbit
satellite
perturbation
thrust
thruster
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
CN201811059085.5A
Other languages
English (en)
Other versions
CN109063380A (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.)
Beijing Institute of Technology BIT
Original Assignee
Beijing Institute of Technology BIT
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 Institute of Technology BIT filed Critical Beijing Institute of Technology BIT
Priority to CN201811059085.5A priority Critical patent/CN109063380B/zh
Publication of CN109063380A publication Critical patent/CN109063380A/zh
Application granted granted Critical
Publication of CN109063380B publication Critical patent/CN109063380B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Geometry (AREA)
  • General Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Control Of Position, Course, Altitude, Or Attitude Of Moving Bodies (AREA)

Abstract

本发明属于航天器轨道动力学与控制领域。本发明公开的一种静止轨道电推进卫星故障检测方法，实现方法如下，建立静止轨道卫星受环境摄动下的瞬时轨道运动模型；建立平均轨道运动模型；建立受控制力作用下的轨道运动模型；推导受所有力作用下的轨道运动模型；设计故障检测器，并验证所设计故障检测器的收敛性。本发明还公开的一种静止轨道电推进卫星故障模式的位置保持方法，在卫星故障检测方法基础上，计算推力器故障导致的轨道漂移；设计控制律，通过控制律修正推力器故障导致的轨道漂移；离散控制律，获得电推力器的开关序列；通过采用开关序列，实现静止轨道电推进卫星故障模式下的位置保持。本发明普遍适用于高轨及地球静止轨道电推进卫星。

Description

一种静止轨道电推进卫星故障检测方法及位置保持方法
技术领域
本发明涉及一种静止轨道电推进卫星故障检测方法及位置保持方法,尤其涉及静止轨道电推进卫星的推力故障检测、位置保持等重要技术，属于航天器轨道动力学与控制领域。
背景技术
静止轨道(Geostationary Orbit,GEO)指轨道高度为35786km、倾角为零的特殊地球同步轨道。在静止轨道上运行的航天器，通常为高成本、大功率、长寿命、高精度的大型复杂航天器，具有覆盖面积大，且相对于地面静止的特点，在通信、导航、预警、遥感、气象、数据中继等军用和民用领域发挥着重要的作用，经济价值非常之高。
GEO卫星的轨道运动主要受地球非球形摄动、日月引力摄动、太阳光压摄动三种主要摄动的影响。在无控情况下，GEO卫星会在摄动力作用下逐渐漂出定点窗口，一旦离开定点窗口，GEO卫星就会失效成为空间碎片，不仅无法实现自身正常的通信导航等功能，还会危及附近经度的其他卫星。GEO卫星的位置保持技术是指在特定的时刻实施一定的轨道机动使卫星保持在定点经度附近的一类控制方法，该方法能够有效抑制高轨摄动因素引起的漂移现象，提升GEO卫星的在轨运行稳定性。
目前GEO卫星的位置保持方法，按推进能源的不同分为两类，一类是基于化学推进的位置保持，另一类是基于电推进的位置保持。电推进系统最显著的优点是高比冲，即等量燃料产生的推力更大，能够有效节省燃料，从而减少发射成本、增加有效载荷的质量。因此，本发明针对的对象是搭载电推进系统的静止轨道卫星。
根据美国国家航空航天局(NASA)提供的电推力器测试报告，以及由佐治亚理工大学Saleh教授团队提供的近15年电推进任务统计报告，电推进系统在卫星在轨服务期间几乎无法避免故障情况。推进系统的故障将直接导致航天器部分或完全失效，不仅无法完成既定空间任务，还会威胁到空间环境的安全。因此，故障检测系统是保证卫星正常在轨运行的必备环节。
故障模式下的位置保持策略包含两个部分：第一，检测推力器的故障状态，跟踪推力器输出推力的变化；第二，关闭故障推力器，并消除推力器故障引起的轨道漂移。当前文献仅针对完备模式下的位置保持策略或是仅针对故障状态已知的位置保持策略进行了研究。文献(M.Guelman,Geostationary satellites autonomous closed loop stationkeeping,Acta Astronautica,97,2014,pp.9–15.)研究了以当前时刻轨道信息为输入的闭环位置保持策略，但前提是推力器正常工作。文献(J.Zhang,S.Zhao,Z.Zhou,X.Li,OptimalStation Keeping for XIP Thrusters in Failure Mode under Eclipse Constraints,Journal of Aerospace Engineering,29(6),2016,pp.1-11.)研究了考虑日蚀情况下的最优燃料位置保持策略，但假设推力器故障状态已知。目前的位置保持策略没有考虑推力器故障检测环节，更缺乏一套具有故障检测和漂移消除两大功能的故障模式位置保持方案。因此本发明旨在提出一种静止轨道电推进卫星故障检测方法及位置保持方法，该方法不仅能精确检测推力器的故障状态，而且在有效识别故障状态之后，能够快速消除由推力器故障引起的轨道漂移，从而实现故障模式下的位置保持。
发明内容
本发明公开的一种静止轨道电推进卫星故障检测方法要解决的技术问题是：提供一种静止轨道电推进卫星故障检测方法，实现静止轨道电推进卫星故障检测，能够有效跟踪任意形式的推力变化，普遍适用于高轨及地球静止轨道电推进卫星。
进一步的，在本发明公开的一种静止轨道电推进卫星故障检测方法基础上，本发明还公开的一种静止轨道电推进卫星故障模式的位置保持方法要解决的技术问题是：提供一种静止轨道电推进卫星故障模式的位置保持方法，实现静止轨道电推进卫星在故障模式下的位置保持，能够消除任意形式的推力变化引起的轨道漂移，普遍适用于高轨及地球静止轨道电推进卫星。
本发明的目的通过以下技术方案实现。
本发明公开的一种静止轨道电推进卫星故障检测方法，通过Lagrange行星运动方程建立静止轨道卫星受环境摄动下的瞬时轨道运动模型；通过平均法建立静止轨道卫星受环境摄动下的平均轨道运动模型；通过Gauss摄动方程和推力器布局建立静止轨道卫星受控制力作用下的轨道运动模型；推导静止轨道卫星受所有力作用下的轨道运动模型；设计故障检测器，并验证所设计故障检测器的收敛性；
本发明还公开的一种静止轨道电推进卫星故障模式的位置保持方法，在所述的一种静止轨道电推进卫星故障检测方法基础上，还包括如下步骤：通过滤波方法计算推力器故障导致的轨道漂移；设计控制律，通过控制律修正推力器故障导致的轨道漂移；离散控制律，获得电推力器的开关序列；通过采用开关序列，实现静止轨道电推进卫星故障模式下的位置保持。
本发明公开的一种静止轨道电推进卫星故障检测方法，包括如下步骤：
步骤一：通过Lagrange行星运动方程建立静止轨道卫星受环境摄动下的瞬时轨道运动模型，所述的环境摄动包括地球非球形摄动RENP、日月引力摄动RLSP、太阳光压摄动RSRP。
步骤1.1：建立无奇点轨道要素；
相对于经典轨道要素{a,e,i,Ω,ω,M}，无奇点轨道要素{a,ex,ey,ix,iy,λ}，能避免轨道倾角i和偏心率e极小情况下的奇异问题。无奇点轨道要素和经典轨道要素的映射关系如下所示：
Figure BDA0001796600240000031
其中，a是半场轴，e是偏心率，i是轨道倾角，ω是近地点角距，Ω是升交点赤经，M是平近点角，θG是格林尼治恒星时角。
步骤1.2：建立Lagrange行星运动方程；
静止轨道电推进卫星的轨道运动主要受三种环境摄动影响。三种自然摄动包括地球非球形摄动RENP、日月引力摄动RLSP和太阳光压摄动RSRP。
Lagrange行星运动方程用于描述静止轨道卫星在三种自然摄动作用下的轨道运动。采用Lagrange行星运动方程并略去二阶小量后的轨道运动模型如下所示：
Figure BDA0001796600240000041
其中，n是卫星轨道运动的平均角速率，R是摄动势函数，Ωe是地球自转角速率。
步骤1.3：建立地球非球形摄动势函数；
静止轨道非球形摄动(Earth’s Non-spherical Perturbation,ENP)势函数具有如下形式：
Figure BDA0001796600240000042
其中，μ是地球引力常数，r是航天器与地心的距离，Re是地球半径，Jn是带谐项系数，Jnm是田谐项系数，Pn(φ),Pnm(φ)是勒让德多项式，φ是地心纬度，λnm是引力场系数。
Pn(φ),Pnm(φ)的表达式如公式(4)所示，Pn(φ)视为Pn0(φ)
Figure BDA0001796600240000043
作为优选，步骤1.3所述的地球非球形摄动势函数优选仅包含前四阶带谐项J2,J3,J4和前三阶田谐项J22,J31,J33的摄动势函数，所述摄动势函数具体推导方法如下：
包含前四阶带谐项J2,J3,J4和前三阶田谐项J22,J31,J33的摄动势函数如公式(5)所示：
Figure BDA0001796600240000051
Figure BDA0001796600240000052
因为静止轨道的轨道倾角较小，所以地心纬度φ满足如下条件：
Figure BDA0001796600240000053
将公式(7)带入公式(5)和公式(6)，获得考虑静止轨道特点的摄动势函数如公式(8)所示：
Figure BDA0001796600240000054
步骤1.4：建立日月引力摄动势函数；
日月引力摄动(Luni-Solar Perturbation,LSP)的势函数具有如下形式：
Figure BDA0001796600240000055
下标k代表第三引力体，k＝s,m，s代表太阳，m代表月亮；μk是第三引力体的引力常数，rk是第三引力体与地心的距离，θk是rk和r的夹角。
步骤1.5：建立太阳光压摄动势函数；
太阳光压摄动(Solar Radiation Perturbation,SRP)的势函数具有如下形式：
Figure BDA0001796600240000061
CR是光压系数，S是卫星正对太阳的面积，m是卫星质量，P是单位面积的太阳辐射压，x,y,z是卫星相对地心的位置矢量在地心惯性系下的三个分量，xs,ys,zs是太阳相对地心的位置矢量在地心惯性系下的三个分量。
步骤1.6：将公式(3)、(9)、(10)分别带入公式(2)推导三种摄动作用下的瞬时轨道要素变化率：
由于公式(3)带入公式(10)后的表达式非常复杂，仅给出如公式(8)所示的优选的非球形摄动势函数带入公式(10)的结果，即非球形摄动(ENP)引起的瞬时轨道要素变化率如下：
Figure BDA0001796600240000062
其中，ac是平半长轴，取值与静止轨道卫星的定点经度有关，nc是ac对应的轨道运动平均角速率。
Figure BDA0001796600240000063
日月引力摄动(LSP)引起的瞬时轨道要素变化率如下：
Figure BDA0001796600240000071
其中，xk,yk,zk是第三引力体在地心惯性坐标下的三个分量。
太阳光压摄动(SRP)引起的瞬时轨道要素变化率如下：
Figure BDA0001796600240000072
步骤1.7：联立公式(11)、(13)、(14)，建立静止轨道卫星受环境摄动下的瞬时轨道运动模型，如公式(15)所示：
Figure BDA0001796600240000073
其中，
Figure BDA0001796600240000074
表示瞬时轨道要素变化率，
Figure BDA0001796600240000075
步骤二：在步骤一建立的静止轨道卫星受环境摄动下的瞬时轨道运动模型基础上，通过平均法计算三种环境摄动引起的平均轨道要素和短周期项，建立静止轨道卫星受环境摄动下的平均轨道运动模型。
步骤2.1：分别推导三种摄动作用下的平均轨道要素和短周期项；
分析公式(11)、(13)、(14)各项的周期性，定义周期大于一个静止轨道轨道周期的项为长周期项，周期小于一个静止轨道轨道周期的项为短周期项；并对公式(11)、(13)、(14)按照公式(16)求平均，即能够推导出三种摄动作用下平均轨道要素变化率
Figure BDA0001796600240000081
按照公式(17)能够推导三种摄动作用下短周期项变化率
Figure BDA0001796600240000082
Figure BDA0001796600240000083
Figure BDA0001796600240000084
积分短周期项变化率
Figure BDA0001796600240000085
能够获得短周期项(σsp)ENP,(σsp)LSP,(σsp)SRP。
三种环境摄动作用下的平均轨道要素变化率和短周期项的表达式如下：
非球形摄动(ENP)引起的平均轨道要素变化率和短周期项如下：
Figure BDA0001796600240000086
Figure BDA0001796600240000091
其中，
Figure BDA0001796600240000092
是半长轴和平经度的短周期修正项；
Figure BDA0001796600240000093
日月引力摄动(LSP)引起的平均轨道要素变化率和短周期项如下：
Figure BDA0001796600240000094
Figure BDA0001796600240000101
其中，CΘ-LSP＝cos(ωk+Mk)cosΩk，SΘ-LSP＝sin(ωk+Mk)sinΩkcosik。
太阳光压摄动(SRP)引起的平均轨道要素变化率和短周期项如下：
Figure BDA0001796600240000102
Figure BDA0001796600240000103
步骤2.2：联立公式(18)-(23)，获得静止轨道卫星受环境摄动下的平均轨道运动模型如公式(24)。
Figure BDA0001796600240000104
步骤三：通过Gauss摄动方程和推力器布局建立静止轨道卫星受控制力作用下的轨道运动模型，并将模型线性化；
步骤3.1：建立Gauss摄动方程；
除三种环境摄动外，静止轨道卫星的轨道运动还受到推力器输出推力的英雄。Gauss摄动方程用于描述静止轨道卫星在电推进系统提供的小推力作用下的轨道运动。采用Gauss摄动方并略去二阶小量后的轨道运动模型如下所示：
Figure BDA0001796600240000111
T是推力器输出的推力大小；dR,dT,dN是推力器在卫星质心轨道坐标系下沿径向、切向和法向的推力分量与T的比值；u是电推力器开关常数，u＝0or 1，0表示推力器关闭，1表示推力器开启；m是静止轨道卫星的总质量；V是静止轨道卫星的速度；l是静止轨道卫星的平赤经，l＝Ω+ω+M。
步骤3.2：将推力器按推力器布局投影到卫星质心轨道坐标系并带入到Gauss摄动方程。
步骤3.3：建立静止轨道卫星受控制力作用下的轨道运动模型，如公式(26)所示：
Figure BDA0001796600240000121
其中，Ti表示第i个推力器输出的推力大小，ui是第i个推力器的开关常数，
Figure BDA0001796600240000122
分别表示第i个推力器在卫星质心轨道坐标系下沿径向、切向和法向的推力分量与Ti的比值。
步骤3.4：线性化建立的静止轨道卫星受控制力作用下的轨道运动模型；
将公式(26)写成矩阵形式，如公式(27)所示：
Figure BDA0001796600240000123
其中fEP(σ)是关于当前轨道六要素的非线性函数矩阵，具有6×3维；u是3×1维的常值控制向量，表征推力大小在质心轨道坐标系径、切、法向的分量，
Figure BDA0001796600240000124
将公式(27)在期望的卫星定点位置σd按泰勒级数展开并保留一阶项，如公式(28)所示：
Figure BDA0001796600240000125
其中，因为
Figure BDA0001796600240000126
和
Figure BDA0001796600240000127
是仅取决于期望定点位置的常值矩阵，因此B(σ)是关于σ的线性函数矩阵，公式(28)是静止轨道卫星受控制力作用下的线性化轨道运动模型。
步骤四：基于步骤二建立的静止轨道卫星受环境摄动下的平均轨道运动模型和步骤三建立的静止轨道卫星受控制力作用下的轨道运动模型，推导静止轨道卫星受所有力作用下的轨道运动模型。
联立步骤2.2的公式(24)和步骤3.4的公式(28)，获得静止轨道卫星受所有力作用下的轨道运动模型如公式所示。
Figure BDA0001796600240000131
其中，B(σ)是关于σ的线性函数矩阵，w和v分别表示表示系统噪声和量测噪声。
步骤五：设计故障检测器，并验证所设计故障检测器的收敛性；
步骤5.1：设计故障检测器，通过故障检测器跟踪推力器故障引起的推力变化。
推力器故障会引起静止轨道卫星轨道漂移，如果不能及时检测故障状态并且消除故障引起的漂移，一旦漂出定点窗口，静止轨道卫星将无法维持正常的通信、导航功能。导致推力器故障主要有两方面的原因，一是自身系统老化；二是外界因素的影响，包括非合作目标的碰撞、干扰。电推力器在理想工作状态输出常值推力，因此，其故障形式主要表现为实际输出推力低于额定推力。
当推力器处于故障状态时，实际输出推力由已知变为未知。离散并反解如公式(29)所示的状态方程能够获得实际推力u与状态量
Figure BDA0001796600240000132
的映射关系，其中
Figure BDA0001796600240000133
属于超前状态信息。
Figure BDA0001796600240000134
引入辅助变量ξ用于观测故障状态下的实际推力u。
Figure BDA0001796600240000135
引入线性观测器用于描述估计推力
Figure BDA0001796600240000136
的递推关系。
Figure BDA0001796600240000137
设计系数矩阵Γ和Ψ实现估计推力
Figure BDA0001796600240000141
对实际推力u的有效跟踪，设计的参数如下：
Figure BDA0001796600240000142
将公式(31)和公式(32)带入公式(30)，获得辅助变量ξ的递推表达式，如公式(34)所示。
Figure BDA0001796600240000143
当且仅当Γ＝ΨB-1，辅助变量ξ的递推只取决于前一时刻的轨道信息，不需要对未来轨道信息进行预测，如公式(35)所示。
Figure BDA0001796600240000144
公式(31)和公式(35)共同组成所设计故障检测器的动力学方程，通过故障检测器能够跟踪推力器故障引起的推力变化。
步骤5.2：验证所设计故障检测器的收敛性；
定义如公式(36)所示的推力估计误差：
Figure BDA0001796600240000145
将公式(31)和公式(35)带入公式(36)：
Figure BDA0001796600240000146
在公式(37)右端增加uk构造ek，获得推力估计误差的递推表达式：
Figure BDA0001796600240000151
当k趋于无穷时，ek趋于零。因此，推力估计误差是渐进稳定的，且最终收敛于零。当Ψ＜I，即λi＜1,i＝1,2,3时，所设计的故障检测器具有良好的局部性能，即能够验证所设计的故障检测器能够使推力估计误差渐进稳定，且最终收敛于零。
在所述的一种静止轨道电推进卫星故障检测方法基础上，本发明还公开的一种静止轨道电推进卫星故障模式的位置保持方法，包括所述的一种静止轨道电推进卫星故障检测方法的步骤一至五，还包括步骤六至步骤九。
步骤六：通过滤波方法计算推力器故障导致的轨道漂移；
步骤五已验证估计推力
Figure BDA0001796600240000152
能够精确跟踪故障状态下的实际推力u，因此，估计推力
Figure BDA0001796600240000153
能够视为已知输入带入如公式(29)所示的静止轨道卫星的系统方程。针对输入已知的非线性系统，应用滤波方法即能够估计故障状态下的平均轨道要素
Figure BDA0001796600240000154
当平均轨道要素估计值
Figure BDA0001796600240000155
与平均轨道要素实际值
Figure BDA0001796600240000156
的误差趋近于零时，当前时刻q的平均轨道要素估计值
Figure BDA0001796600240000157
与理想平均轨道要素
Figure BDA0001796600240000158
之间的偏差即为推力故障持续期间引起的轨道漂移。
步骤六所述的滤波方法包括无迹卡尔曼滤波方法、扩展卡尔曼滤波。
作为优选，步骤六所述的滤波方法优选无迹卡尔曼滤波方法，针对形如公式(39)的系统方程，无迹卡尔曼滤波算法(UKF)的流程如下：
Figure BDA0001796600240000159
流程1：计算采样点集。
根据UT(unscented transform)变换选取采样点：
Figure BDA00017966002400001510
其中，Sx是Cholesky因子，初值为
Figure BDA00017966002400001511
函数chol(·)是Cholesky分解函数，Zi是第i个采样点对应的向量空间。
流程2：状态更新。
Figure BDA0001796600240000161
Figure BDA0001796600240000162
其中，函数qr(·)为QR分解函数，函数cholupdate(·)为Cholesky分解一次更新。
Figure BDA0001796600240000163
为协方差加权的权值，
Figure BDA0001796600240000164
为均值加权的权值。
流程3：量测更新：
Figure BDA0001796600240000165
其中，操作符A/B为矩阵相除，通过求解方程AATx＝ATB获得。
步骤七：设计控制律，通过控制律修正推力器故障导致的轨道漂移。
步骤七所述控制律包括开环控制律和闭环控制律，开环控制律基于最优控制理论，闭环控制律基于反馈信息。
作为优选，步骤七所述控制律优选反馈控制律，步骤七具体实现方法如下：
将步骤六获得的轨道漂移作为初始偏差，使用如下形式控制律修正偏差：
Figure BDA0001796600240000166
Figure BDA0001796600240000167
是期望轨道面沿地心轨道坐标系法向的单位向量，
Figure BDA0001796600240000168
是当前轨道面沿地心轨道坐标系径向的单位向量，
Figure BDA0001796600240000171
是当前时刻航天器的速度矢量。c1,c2是增益，
Figure BDA0001796600240000172
ac＝as+ka(a-ad)，as是标准静止轨道的半长轴，ka是设计的常值系数，ad是期望半长轴。
先给出由轨道要素描述的
Figure BDA0001796600240000173
在地心惯性系(ECI)的投影。
Figure BDA0001796600240000174
Figure BDA0001796600240000175
Figure BDA0001796600240000176
其中，
Figure BDA0001796600240000177
Figure BDA0001796600240000178
通过对公式(48)进行周期性分析，将包含sinkl,coskl且k≥2的项视为短周期项，去除短周期项后的表达式如下：
Figure BDA0001796600240000179
通过对公式(47)进行周期性分析，将包含sin2M,cos2M的项视为短周期项，去除短周期项后的表达式如下：
Figure BDA0001796600240000181
其中，u＝Ω+ω，函数ysg＝fsg(xsg),ycg＝fcg(xcg)表示(xsg,ysg)和(xcg,ycg)具有相同的周期性。
将公式(47)、(49)、(50)带入公式(44)，获得控制律在地心惯性系下的表达式，如公式(51)所示。
Figure BDA0001796600240000182
将公式(51)所示的
Figure BDA0001796600240000183
投影到航天器质心轨道坐标系，即得到能够直接带入如公式(29)所示系统方程的控制律。
Figure BDA0001796600240000184
其中，
Figure BDA0001796600240000185
步骤八：离散步骤七所述的控制律，获得电推力器的开关序列；
公式(52)给出了连续型的控制律，但由于电推力器输出常值推力且具有开关特性，因此需要根据推力器的性能参数将连续型的控制律离散化。
通过计算推力器常值推力产生的加速度arate，设计以arate为阈值的函数，如公式(53)：
Figure BDA0001796600240000191
联立公式(52)和公式(53)，即能够获得推力器的开关序列。
步骤九：通过采用步骤八所述的开关序列，能够消除任意形式的推力变化引起的轨道漂移，实现静止轨道电推进卫星故障模式下的位置保持。
有益效果
1、本发明公开的一种静止轨道电推进卫星故障检测方法，所有运行在高轨及静止轨道的航天器须考虑三种环境摄动，推导的平均轨道运动模型考虑了三种主要环境摄动的影响，且推导的公式不受摄动势函数阶数和次数的限制，因此适用于所有运行在高轨及静止轨道的航天器，通过该方法对静止轨道卫星进行控制，为静止轨道卫星在轨运行的稳定性和安全性提供了理论依据和技术支撑。
2、本发明公开的一种静止轨道电推进卫星故障检测方法，通过设计辅助变量和观测器能够精确估计故障模式下的推力，且估计的推力只取决于前一时刻的轨道信息，不需要对未来轨道信息进行预测，能够简化模型、提升计算效率。
3、本发明公开的一种静止轨道电推进卫星故障检测方法，通过设计辅助变量和观测器能够精确估计故障模式下的推力，进一步的，通过合理选取检测器的参数，能够有效跟踪任意形式的推力变化，从而有效应对卫星在轨运行过程中的推力故障情况。
4、本发明还公开的一种静止轨道电推进卫星故障模式的位置保持方法，在所述的一种静止轨道电推进卫星故障检测方法基础上实现，因此，具有所述的一种静止轨道电推进卫星故障检测方法的优点，此外，还包括步骤六至步骤九；通过滤波方法计算推力器故障导致的轨道漂移；设计控制律，通过控制律修正推力器故障导致的轨道漂移；离散所设计的控制律，获得电推力器的开关序列；通过采用所设计的开关序列，能够消除任意形式的推力变化引起的轨道漂移，从而实现静止轨道电推进卫星故障模式下的位置保持。
5、本发明还公开的一种静止轨道电推进卫星故障模式的位置保持方法，通过合理选取采样点并应用滤波方法，和建立的动力学模型具有良好的适配性，能够有效估计故障模式下的轨道漂移。
6、本发明还公开的一种静止轨道电推进卫星故障模式的位置保持方法，通过设计控制律并对其进行离散处理，能够获得电推进系统的推力开关序列，普遍适用于所有搭载电推进系统的静止轨道卫星的位置保持任务。
7、本发明还公开的一种静止轨道电推进卫星故障模式的位置保持方法，通过设计反馈控制律，并合理选取控制器参数，能够消除任意形式的推力变化引起的轨道漂移，相比正常模式下的位置保持策略精度有所提高。
附图说明
图1为一种静止轨道电推进卫星故障检测方法的流程图；
图2为一种静止轨道电推进卫星故障模式的位置保持方法的流程图；
图3为故障模式下实现位置保持的推力器布局；
图4为所设计故障检测器对故障推力的跟踪情况；
图5为故障模式下滤波算法估计的平均轨道要素与实际轨道要素的偏差变化，其中：图5(a)为半长轴偏差的变化、图5(b)为平经度偏差的变化、图5(c)为偏心率矢量x分量偏差的变化、图5(d)为偏心率矢量y分量偏差的变化、图5(e)为倾角矢量x分量偏差的变化、图5(f)为倾角矢量y分量偏差的变化；
图6为所设计故障模式位置保持控制律在切向和法向提供的控制加速度；
图7为所设计故障模式位置保持控制律在切向提供的控制加速度及其对应的推力器开关序列所提供的控制加速度；
图8为无控、理想、受控三种工况下的轨道要素变化，其中：图8(a)为半长轴变化、图8(b)为平经度变化、图8(c)为偏心率矢量x分量的变化、图8(d)为偏心率矢量y分量的变化、图8(e)为倾角矢量x分量的变化、图8(f)为倾角矢量y分量的变化。
具体实施方式
为了更好地说明本发明的目的和优点，下面结合附图对本发明的具体实施方式作进一步详细的说明。
实施例1：
如图1所示，本实施例公开的一种静止轨道电推进卫星故障检测方法，为验证该方法，首先，选取运行在地球静止轨道上的卫星作为主要研究对象。仿真的起始时刻为2020年1月1日0时0分0秒(世界时)，卫星的基本参数如下表所示。
表1卫星参数
Figure BDA0001796600240000201
Figure BDA0001796600240000211
步骤一：通过Lagrange行星运动方程建立静止轨道卫星受环境摄动下的瞬时轨道运动模型，所述的环境摄动包括地球非球形摄动RENP、日月引力摄动RLSP、太阳光压摄动RSRP。；
经典轨道要素在描述小倾角和小偏心率的GEO卫星的轨道运动时，容易出现奇异问题。为避免奇异，引入无奇点轨道要素，如公式(1)所示。静止轨道电推进卫星所受主要环境摄动有三种，分别是地球非球形摄动RENP、日月引力摄动RLSP、太阳光压摄动RSRP。用略去二阶小量的Lagrange行星运动方程如公式(2)所示，描述环境摄动力和控制力作用下的轨道运动模型。
用于计算地球非球形摄动RENP的地球引力场系数如下表所示：
表2地球引力场主项系数
Figure BDA0001796600240000212
用于计算日月引力摄动RLSP、太阳光压摄动RSRP的太阳、月亮平均轨道运动模型如下。
太阳平均轨道运动模型如公式(54)所示：
Figure BDA0001796600240000221
月亮平均轨道运动模型如公式(55)所示：
Figure BDA0001796600240000222
其中，T是当前时刻对应的儒略日，
Figure BDA0001796600240000223
将三种环境摄动的势函数带入Lagrange行星运动方程，即能够建立静止轨道卫星受环境摄动下的瞬时轨道运动模型，如公式(15)所示。
步骤二：在步骤一建立的静止轨道卫星受环境摄动下的瞬时轨道运动模型基础上，通过平均法计算三种环境摄动引起的平均轨道要素和短周期项，建立静止轨道卫星受环境摄动下的平均轨道运动模型。
建立的静止轨道卫星受环境摄动下的平均轨道运动模型如公式(24)所示。因为推导的平均轨道运动模型考虑了三种主要环境摄动的影响，且推导的公式不受摄动势函数阶数和次数的限制，因此适用于所有运行在高轨及静止轨道的航天器，
步骤三：通过Gauss摄动方程和推力器布局建立静止轨道卫星受控制力作用下的轨道运动模型，并将模型线性化；
实例中的卫星搭载电推进系统，用于实现故障模式下的位置保持。推力器布局如图3所示。引入卫星质心轨道坐标系RTN用于描述推力器布局，
Figure BDA0001796600240000224
轴由卫星质心指向地心，
Figure BDA0001796600240000225
轴与
Figure BDA0001796600240000226
轴相互垂直，指向卫星速度的正方向，
Figure BDA0001796600240000227
轴垂直于轨道面，方向满足右手定则。四个推力器正装在
Figure BDA0001796600240000228
轴和
Figure BDA0001796600240000229
轴的正负方向上，沿顺时针方向依次记作N(a)、E(b)、S(c)、W(d)。
推力器参数如下表所示：
表3推力器参数
Figure BDA0001796600240000231
建立的静止轨道卫星受控制力作用下的线性化轨道运动模型如公式(28)。
步骤四：基于步骤二建立的静止轨道卫星受环境摄动下的平均轨道运动模型和步骤三建立的静止轨道卫星受控制力作用下的轨道运动模型，推导静止轨道卫星受所有力作用下的轨道运动模型。
建立的静止轨道卫星受所有力作用下的轨道运动模型如公式(29)所示。
本实例设置系统噪声w为高斯白噪声，方差为Q；量测噪声v为高斯白噪声，方差为R。
Q＝diag(1×10-8,1×10-16,1×10-16,1×10-16,1×10-16,1×10-14) (56)
R＝diag(500.4,8.966×10-12,8.966×10-12,6.126×10-10,6.126×10-10,1.013×10-4) (57)
步骤五：设计故障检测器，并验证所设计故障检测器的收敛性；
所设计故障检测器由公式(31)和公式(35)共同组成。
对比公式(30)和公式(35)，实际推力u的计算不仅依赖当前时刻的轨道信息
Figure BDA0001796600240000233
还依赖未来的轨道信息
Figure BDA0001796600240000232
但估计的推力只取决于前一时刻的轨道信息，不需要对未来轨道信息进行预测，因此所设计的故障检测不仅能够精确跟踪故障状态下的实际推力u，还能够简化模型、提升计算效率。
本实例采用的故障检测器的参数如下：
Ψ＝diag[λ,λ,λ] (58)
λ分别取λ＝1.35，λ＝0.85和λ＝0.55做了三组对比仿真，仿真结果如图4所示。
图4为所设计故障检测器对故障推力的跟踪情况，红色实线表示实际推力，在第100s和第200s由于故障分别衰减常值推力的20％，在第300s时推力器关闭。绿色虚线、蓝色实线、黑色实线分别代表λ＝1.35，λ＝0.85和λ＝0.55时，所设计故障检测器估计的推力。通过对比能够发现，故障引起的推力变化在10s内即能够实现精确跟踪，且λ＜1时具有良好的局部性能，且越趋近1时跟踪精度越高。本实例中模拟的推力故障情况是间隔一段时间的两次推力衰减且推力衰减幅度递增，但所有可能出现的故障情况都能够拆分成单次推力变化的组合。因此，通过合理选取故障检测器的参数，能够有效跟踪任意形式的推力变化，从而有效应对卫星在轨运行过程中的推力故障情况。
实施例2：
如图2所示，本实施例在所述的一种静止轨道电推进卫星故障检测方法基础上，还公开的一种静止轨道电推进卫星故障模式的位置保持方法，包括所述的一种静止轨道电推进卫星故障检测方法的步骤一至五，还包括步骤六至步骤九，为验证该方法，前五步选取的参数与实施例1一致，步骤六至步骤九如下：
步骤六：通过滤波方法计算推力器故障导致的轨道漂移；
作为优选，步骤六所述的滤波方法优选无迹卡尔曼滤波方法，步骤具体实现方法如公式(39)-(43)所示。
本实例设定状态协方差矩阵的初值如公式(59)所示：
P＝diag(1×103,1×10-5,1×10-5,2×10-3,2×10-3,3×10-2) (59)
将步骤五所设计的故障检测器估计的推力
Figure BDA0001796600240000241
带入如公式(29)的系统方程，并应用上述的无迹卡尔曼滤波方法，估计推力故障导致的轨道漂移。仿真结果如图5所示。
图5表示采用无迹卡尔曼律方法估计的轨道要素与实际轨道要素的偏差。红色实线表示推力器正常工作时的偏差变化，偏差不恒等于零是因为系统误差和观测误差始终存在。蓝色实线表示单次推力故障引起的偏差变化，故障情况如下：推力器自第0s开始正常工作，在第100s时发生推力器故障，故障导致输出推力衰减至峰值的80％，从故障发生到主动关闭故障推力器总共持续100s。绿色实线表示两次推力故障引起的偏差变化，故障情况如下：推力器自第0s开始正常工作，在第100s时发生第一次推力器故障，故障导致推力器输出衰减至峰值的80％，在第200s发生第二次推力器故障，故障导致推力器输出衰减至峰值的60％，从故障发生到主动关闭故障推力器总共持续200s。从检测到故障出现到主动关闭故障推力器的总时长定义为故障持续时间，本例中分明设定100s和200s作为故障持续时间的上限，一旦故障持续时间超过设定时间，故障推力器将会被主动关闭。蓝色实线表示的单次推力故障引起的六个轨道要素偏差均在第1852s收敛，绿色实线表示的两次推力故障引起的六个轨道要素偏差均在第2432s收敛。由此能够分别计算在第1852s和第2432s时卫星所处位置与期望定点窗口的轨道漂移。蓝色实线所表示的单次推力故障引起的轨道要素偏差在1852s收敛，此时轨道要素估计值
Figure BDA0001796600240000251
与平均轨道要素实际值
Figure BDA0001796600240000252
的偏差近似为零，推力故障持续期间引起的轨道漂移为此时平均轨道要素估计值
Figure BDA0001796600240000253
与理想平均轨道要素
Figure BDA0001796600240000254
其值如下：
Figure BDA0001796600240000255
在本实例的两种故障情况下轨道要素的估计误差都能较快且较精确地收敛，证明选用的滤波方法和建立的动力学模型具有良好的适配性，能够有效估计故障模式下的轨道漂移。
步骤七：设计控制律，通过控制律修正推力器故障导致的轨道漂移；
公式(60)表示推力故障持续期间引起的轨道漂移，将其值作为初始偏差，使用如公式(52)所示的反馈控制律修正偏差。在所设计的控制律中，常值系数ka＝2000。仿真结果如图6所示。
步骤八：离散步骤七所述的控制律，获得电推力器的开关序列；
因为电推力器的输出推力通常为常值且具有开关特性，所以需要将所设计的连续控制律离散以获得实际推力器的开关序列。根据表2提供的推力器的性能参数，又因为四个推力器分别正装在切向和法向轴的正负方向，因此选取如下加速度阈值arate＝200mN/2000kg＝1×10-4m/s2。联立公式(53)和公式(52)，即能够获得推力器的开关序列。仿真结果如图7所示。
图6表示步骤七所设计的控制律在卫星质心轨道坐标系切向和法向提供的控制加速度。图7表示离散处理后获得的推力开关序列。图7所示矩形波的正向信号由安装在切向轴，推力方向为正向的推力器提供；矩形波的负向信号由安装在切向轴，推力方向为负向的推力器提供；对称安装的推力器在控制周期内交替运行，这种开关特性同样适用于对称安装在法向轴的推力器上。图6所示的切向加速度和径向加速度在一个控制周期内具有重叠部分，由此表明单个周期内存在两个推力器同时开机的情况。电推进卫星的寿命周期内，为实现燃料最省的目的，工程中通常使用“单个时刻仅有一个推力器开机”的约束条件，但在故障模式下，为保证卫星始终位于定点窗口，首要目标是在最短的时间内消除故障引起的轨道漂移，推力器的开机约束能够适当放宽，因此本发明设计的故障模式开机策略是合理的。因为所设计的控制律都通过设置阈值离散为电推力器的开关序列，因此本发明提出的方法普遍适用于所有搭载电推进系统的静止轨道卫星的位置保持任务
步骤九：通过采用步骤八所述的开关序列，能够消除任意形式的推力变化引起的轨道漂移，实现静止轨道电推进卫星故障模式下的位置保持。
将步骤八获得的推力器开关序列带入到如公式(29)所示的系统动力学方程，即能够获得卫星的轨道要素在受控情况下的变化规律，仿真结果如图8所示。
图8表示初始偏差一定的无控、理想、受控三种工况下的对比仿真结果。蓝色实现表示无控情况的轨道要素变化情况，黑色实现表示理想静止轨道卫星的轨道要素变化情况，绿色实线表示采用本发明设计的开关序列后的轨道要素变化情况。在收敛时间方面，六个轨道要素均能在2天(两个控制周期)内收敛。偏心率矢量(ex,ey)在当前算例中的收敛时间为0.5天，倾角矢量(ix,iy)在当前算例中的收敛时间为1天，具有良好的性能。在收敛精度方面，平半长轴
Figure BDA0001796600240000261
收敛至100米内，平经度
Figure BDA0001796600240000262
收敛至10-3度，平偏心率矢量
Figure BDA0001796600240000263
收敛至10-5，平倾角矢量
Figure BDA0001796600240000264
收敛至10-5。表给出了本发明设计的位保策略同主要参考文献的位保策略控制精度的对比。方案1代表本发明的结果；方案2代表文献(M.Guelman,Geostationary satellites autonomous closed loop station keeping,ActaAstronautica,97,2014,pp.9–15.)的结果；方案3代表文献(J.Zhang,S.Zhao,Z.Zhou,X.Li,Optimal Station Keeping for XIP Thrusters in Failure Mode under EclipseConstraints,Journal of Aerospace Engineering,29(6),2016,pp.1-11.)的结果。因此，本实施例公开的一种静止轨道电推进卫星故障模式的位置保持方法(方案1)，相较推力器正常工作的位置保持策略(方案2)和已知一对推力器故障的位置保持策略(方案3)，控制精度均有1到2个数量级的提升，具有更优的控制性能。
表4主流控制方案的精度对比
Figure BDA0001796600240000265
综上所述，以上仅为本发明的较佳实施例而已，并非用于限定本发明的保护范围。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (8)

1.一种静止轨道电推进卫星故障检测方法，其特征在于：包括如下步骤，
步骤一：通过Lagrange行星运动方程建立静止轨道卫星受环境摄动下的瞬时轨道运动模型，所述的环境摄动包括地球非球形摄动RENP、日月引力摄动RLSP、太阳光压摄动RSRP；
步骤二：在步骤一建立的静止轨道卫星受环境摄动下的瞬时轨道运动模型基础上，通过平均法计算三种环境摄动引起的平均轨道要素和短周期项，建立静止轨道卫星受环境摄动下的平均轨道运动模型；
步骤三：通过Gauss摄动方程和推力器布局建立静止轨道卫星受控制力作用下的轨道运动模型，并将模型线性化；
步骤四：基于步骤二建立的静止轨道卫星受环境摄动下的平均轨道运动模型和步骤三建立的静止轨道卫星受控制力作用下的轨道运动模型，推导静止轨道卫星受所有力作用下的轨道运动模型；静止轨道卫星受所有力作用下的轨道运动模型如公式(25)所示；
Figure FDA0003858699380000011
其中，B(σ)是关于σ的线性函数矩阵，w和v分别表示表示系统噪声和量测噪声；
步骤五：设计故障检测器，并验证所设计故障检测器的收敛性；
步骤五具体实现方法如下，
步骤5.1：设计故障检测器，通过故障检测器跟踪推力器故障引起的推力变化；
当推力器处于故障状态时，实际输出推力由已知变为未知；离散并反解如公式(25)所示的状态方程获得实际推力u与状态量
Figure FDA0003858699380000012
的映射关系，其中
Figure FDA0003858699380000013
属于超前状态信息；
Figure FDA0003858699380000014
引入辅助变量ξ用于观测故障状态下的实际推力u；
Figure FDA0003858699380000021
引入线性观测器用于描述估计推力
Figure FDA0003858699380000022
的递推关系；
Figure FDA0003858699380000023
设计系数矩阵Γ和Ψ实现估计推力
Figure FDA0003858699380000024
对实际推力u的有效跟踪，设计的参数如下：
Figure FDA0003858699380000025
将公式(27)和公式(28)带入公式(26)，获得辅助变量ξ的递推表达式，如公式(30)所示；
Figure FDA0003858699380000026
当且仅当Γ＝ΨB-1，辅助变量ξ的递推只取决于前一时刻的轨道信息，不需要对未来轨道信息进行预测，如公式(31)所示；
Figure FDA0003858699380000027
公式(27)和公式(31)共同组成所设计故障检测器的动力学方程，通过故障检测器跟踪推力器故障引起的推力变化；
步骤5.2：验证所设计故障检测器的收敛性；
定义如公式(32)所示的推力估计误差：
Figure FDA0003858699380000028
将公式(27)和公式(31)带入公式(32)：
Figure FDA0003858699380000031
在公式(33)右端增加uε构造eε，获得推力估计误差的递推表达式：
Figure FDA0003858699380000032
当ε趋于无穷时，eε趋于零；因此，推力估计误差是渐进稳定的，且最终收敛于零；当Ψ＜I，即γ1＜1,γ2＜1,γ3＜1时，所设计的故障检测器具有良好的局部性能，即验证所设计的故障检测器使推力估计误差渐进稳定，且最终收敛于零。
2.如权利要求1所述的一种静止轨道电推进卫星故障检测方法，其特征在于：步骤一具体实现方法如下，
步骤1.1：建立无奇点轨道要素；
相对于经典轨道要素{a,e,i,Ω,ω,M}，无奇点轨道要素{a,ex,ey,ix,iy,λ}，能避免轨道倾角i和偏心率e极小情况下的奇异问题；无奇点轨道要素和经典轨道要素的映射关系如下所示：
Figure FDA0003858699380000033
其中，a是半场轴，e是偏心率，i是轨道倾角，ω是近地点角距，Ω是升交点赤经，M是平近点角，θG是格林尼治恒星时角；
步骤1.2：建立Lagrange行星运动方程；
静止轨道电推进卫星的轨道运动受三种环境摄动影响；三种环境摄动包括地球非球形摄动RENP、日月引力摄动RLSP和太阳光压摄动RSRP；
Lagrange行星运动方程用于描述静止轨道卫星在三种环境摄动作用下的轨道运动；采用Lagrange行星运动方程并略去二阶小量后的轨道运动模型如下所示：
Figure FDA0003858699380000041
其中，nsat是卫星轨道运动的平均角速率，R是摄动势函数，Ωe是地球自转角速率；
步骤1.3：建立地球非球形摄动势函数；
地球非球形摄动(Earth’s Non-sphericalPerturbation,ENP)势函数具有如下形式：
Figure FDA0003858699380000042
其中，μ是地球引力常数，r是航天器与地心的距离，Re是地球半径，Jn是带谐项系数，JNQ是田谐项系数，PN(φ)、PNQ(φ)是勒让德多项式，φ是地心纬度，λNQ是引力场系数；
PN(φ)、PNQ(φ)的表达式如公式(4)所示，PN(φ)视为PN0(φ)
Figure FDA0003858699380000051
步骤1.4：建立日月引力摄动势函数；
日月引力摄动(Luni-Solar Perturbation,LSP)的势函数具有如下形式：
Figure FDA0003858699380000052
下标k代表第三引力体，k＝s,m，s代表太阳，m代表月亮；μk是第三引力体的引力常数，rk是第三引力体与地心的距离，θk是rk和r的夹角；
步骤1.5：建立太阳光压摄动势函数；
太阳光压摄动(SolarRadiationPerturbation,SRP)的势函数具有如下形式：
Figure FDA0003858699380000053
CR是光压系数，S是卫星正对太阳的面积，Mass是卫星质量，P是单位面积的太阳辐射压，x,y,z是卫星相对地心的位置矢量在地心惯性系下的三个分量，xs,ys,zs是太阳相对地心的位置矢量在地心惯性系下的三个分量；
步骤1.6：将公式(3)、(5)、(6)分别带入公式(2)推导三种摄动作用下的瞬时轨道要素变化率：
由于公式(3)带入公式(6)后的表达式非常复杂，仅给出优选的非球形摄动势函数带入公式(6)的结果，即非球形摄动(ENP)引起的瞬时轨道要素变化率如下：
Figure FDA0003858699380000061
其中，ac是平半长轴，取值与静止轨道卫星的定点经度有关，nc是ac对应的轨道运动平均角速率；
Figure FDA0003858699380000062
日月引力摄动(LSP)引起的瞬时轨道要素变化率如下：
Figure FDA0003858699380000063
其中，xk,yk,zk是第三引力体在地心惯性坐标下的三个分量；
太阳光压摄动(SRP)引起的瞬时轨道要素变化率如下：
Figure FDA0003858699380000071
步骤1.7：联立公式(7)、(9)、(10)，建立静止轨道卫星受环境摄动下的瞬时轨道运动模型，如公式(11)所示：
Figure FDA0003858699380000072
其中
Figure FDA0003858699380000073
表示瞬时轨道要素变化率，
Figure FDA0003858699380000074
3.如权利要求1或2所述的一种静止轨道电推进卫星故障检测方法，其特征在于：步骤二具体实现方法如下，
步骤2.1：分别推导三种摄动作用下的平均轨道要素和短周期项；
分析公式(7)、(9)、(10)各项的周期性，定义周期大于一个静止轨道轨道周期的项为长周期项，周期小于一个静止轨道轨道周期的项为短周期项；并对公式(7)、(9)、(10)按照公式(12)求平均，即推导出三种摄动作用下平均轨道要素变化率
Figure FDA0003858699380000075
按照公式(13)推导三种摄动作用下短周期项变化率
Figure FDA0003858699380000076
Figure FDA0003858699380000077
Figure FDA0003858699380000078
积分短周期项变化率
Figure FDA0003858699380000079
获得短周期项(σsp)ENP,(σsp)LSP,(σsp)SRP；
三种环境摄动作用下的平均轨道要素变化率和短周期变化项的表达式如下：
地球非球形摄动(ENP)引起的平均轨道要素变化率和短周期变化项如下：
Figure FDA0003858699380000081
Figure FDA0003858699380000082
其中，
Figure FDA0003858699380000083
是半长轴和平经度的短周期修正项；
Figure FDA0003858699380000091
日月引力摄动(LSP)引起的平均轨道要素变化率和短周期变化项如下：
Figure FDA0003858699380000092
Figure FDA0003858699380000093
其中，CΘ-LSP＝cos(ωk+Mk)cosΩk，SΘ-LSP＝sin(ωk+Mk)sinΩkcosik；太阳光压摄动(SRP)引起的平均轨道要素变化率和短周期变化项如下：
Figure FDA0003858699380000101
Figure FDA0003858699380000102
步骤2.2：联立公式(14)-(19)，获得静止轨道卫星受环境摄动下的平均轨道运动模型如公式(20)；
Figure FDA0003858699380000103
4.如权利要求3所述的一种静止轨道电推进卫星故障检测方法，其特征在于：步骤三具体实现方法如下，
步骤3.1：建立Gauss摄动方程；
除三种环境摄动外，静止轨道卫星的轨道运动还受到推力器输出推力的影响；Gauss摄动方程用于描述静止轨道卫星在电推进系统提供的小推力作用下的轨道运动；采用Gauss摄动方并略去二阶小量后的轨道运动模型如下所示：
Figure FDA0003858699380000111
T是推力器输出的推力大小；dR,dT,dN是推力器在卫星质心轨道坐标系下沿径向、切向和法向的推力分量与T的比值；u是电推力器开关常数，u＝0或1，0表示推力器关闭，1表示推力器开启；m是静止轨道卫星的总质量；V是静止轨道卫星的速度；l是静止轨道卫星的平赤经，l＝Ω+ω+M；
步骤3.2：将推力器按推力器布局投影到卫星质心轨道坐标系并带入到Gauss摄动方程；
步骤3.3：建立静止轨道卫星受控制力作用下的轨道运动模型，如公式(22)所示：
Figure FDA0003858699380000112
其中，Th表示第h个推力器输出的推力大小，uh是第h个推力器的开关常数，
Figure FDA0003858699380000121
分别表示第h个推力器在卫星质心轨道坐标系下沿径向、切向和法向的推力分量与Th的比值；
步骤3.4：线性化建立的静止轨道卫星受控制力作用下的轨道运动模型；
将公式(22)写成矩阵形式，如公式(23)所示：
Figure FDA0003858699380000122
其中fEP(σ)是关于当前轨道六要素的非线性函数矩阵，具有6×3维；u是3×1维的常值控制向量，表征推力大小在质心轨道坐标系径、切、法向的分量，
Figure FDA0003858699380000123
将公式(23)在期望的卫星定点位置σd按泰勒级数展开并保留一阶项，如公式(24)所示：
Figure FDA0003858699380000124
其中，因为
Figure FDA0003858699380000125
和
Figure FDA0003858699380000126
是仅取决于期望定点位置的常值矩阵，因此B(σ)是关于σ的线性函数矩阵，公式(24)是静止轨道卫星受控制力作用下的线性化轨道运动模型。
5.一种静止轨道电推进卫星故障模式的位置保持方法，其特征在于：包括如权利要求1、2、3或4或所述的一种静止轨道电推进卫星故障检测方法的步骤一至五，还包括步骤六至步骤九，
步骤六：通过滤波方法计算推力器故障导致的轨道漂移；
步骤七：设计控制律，通过控制律修正推力器故障导致的轨道漂移；
步骤八：离散步骤七所述的控制律，获得电推力器的开关序列；
步骤九：通过采用步骤八所述的开关序列，消除任意形式的推力变化引起的轨道漂移，实现静止轨道电推进卫星故障模式下的位置保持。
6.如权利要求5所述的一种静止轨道电推进卫星故障模式的位置保持方法，其特征在于：步骤六具体实现方法如下，
步骤五已验证估计推力
Figure FDA0003858699380000131
精确跟踪故障状态下的实际推力u，因此，估计推力
Figure FDA0003858699380000132
视为已知输入带入如公式(25)所示的静止轨道卫星的系统方程；针对输入已知的非线性系统，应用滤波方法即估计故障状态下的平均轨道要素
Figure FDA0003858699380000133
当平均轨道要素估计值
Figure FDA0003858699380000134
与平均轨道要素实际值
Figure FDA0003858699380000135
的误差趋近于零时，当前时刻q的平均轨道要素估计值
Figure FDA0003858699380000136
与理想平均轨道要素
Figure FDA0003858699380000137
之间的偏差即为推力故障持续期间引起的轨道漂移。
7.如权利要求6所述的一种静止轨道电推进卫星故障模式的位置保持方法，其特征在于：步骤七具体实现方法如下，
步骤七所述控制律包括开环控制律和闭环控制律，开环控制律基于最优控制理论，闭环控制律基于反馈信息；当步骤七所述控制律选反馈控制律，步骤七具体实现方法如下，
将步骤六获得的轨道漂移作为初始偏差，使用如下形式控制律修正偏差：
Figure FDA0003858699380000138
Figure FDA0003858699380000139
是期望轨道面沿地心轨道坐标系法向的单位向量，
Figure FDA00038586993800001310
是当前轨道面沿地心轨道坐标系径向的单位向量，
Figure FDA00038586993800001311
是当前时刻航天器的速度矢量；c1,c2是增益，
Figure FDA00038586993800001312
ac＝as+ka(a-ad)，as是标准静止轨道的半长轴，ka是设计的常值系数，ad是期望半长轴；
先给出由轨道要素描述的
Figure FDA00038586993800001313
在地心惯性系(ECI)的投影；
Figure FDA0003858699380000141
Figure FDA0003858699380000142
Figure FDA0003858699380000143
其中，
Figure FDA0003858699380000144
Figure FDA0003858699380000145
通过对公式(39)进行周期性分析，将包含sinkl,coskl且k≥2的项视为短周期项，去除短周期项后的表达式如下：
Figure FDA0003858699380000146
通过对公式(38)进行周期性分析，将包含sin2M,cos2M的项视为短周期项，去除短周期项后的表达式如下：
Figure FDA0003858699380000151
其中，β＝Ω+ω，函数ysg＝fsg(xsg),ycg＝fcg(xcg)表示(xsg,ysg)和(xcg,ycg)具有相同的周期性；
将公式(38)、(40)、(41)带入公式(35)，获得控制律在地心惯性系下的表达式，如公式(42)所示；
Figure FDA0003858699380000152
将公式(42)所示的
Figure FDA0003858699380000153
投影到航天器质心轨道坐标系，即得到直接带入如公式(25)所示系统方程的控制律；
Figure FDA0003858699380000154
其中，
Figure FDA0003858699380000155
8.如权利要求7所述的一种静止轨道电推进卫星故障模式的位置保持方法，其特征在于：步骤八具体实现方法如下，公式(43)给出了连续型的控制律，但由于电推力器输出常值推力且具有开关特性，因此需要根据推力器的性能参数将连续型的控制律离散化；
通过计算推力器常值推力产生的加速度frate，设计以frate为阈值的函数，如公式(44)：
Figure FDA0003858699380000161
联立公式(43)和公式(44)，即获得推力器的开关序列。
CN201811059085.5A 2018-09-12 2018-09-12 一种静止轨道电推进卫星故障检测方法及位置保持方法 Active CN109063380B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201811059085.5A CN109063380B (zh) 2018-09-12 2018-09-12 一种静止轨道电推进卫星故障检测方法及位置保持方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201811059085.5A CN109063380B (zh) 2018-09-12 2018-09-12 一种静止轨道电推进卫星故障检测方法及位置保持方法

Publications (2)

Publication Number Publication Date
CN109063380A CN109063380A (zh) 2018-12-21
CN109063380B true CN109063380B (zh) 2023-02-17

Family

ID=64760325

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201811059085.5A Active CN109063380B (zh) 2018-09-12 2018-09-12 一种静止轨道电推进卫星故障检测方法及位置保持方法

Country Status (1)

Country Link
CN (1) CN109063380B (zh)

Families Citing this family (10)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109977463A (zh) * 2019-02-15 2019-07-05 南京航空航天大学 一种使用缩比模型测量大型装备红外辐射特征的相似实验方法
CN110209190B (zh) * 2019-03-01 2022-05-20 苏州纳飞卫星动力科技有限公司 一种卫星标称轨道无偏飞行控制的方法
CN110053788B (zh) * 2019-03-15 2022-05-13 中国西安卫星测控中心 一种考虑复杂摄动的星座长期保持控制频次估计方法
CN111989265B (zh) * 2019-11-26 2022-08-12 中国科学院微小卫星创新研究院 超低轨道卫星轨道自主维持方法
CN111854765B (zh) * 2020-06-08 2022-04-26 中国人民解放军战略支援部队航天工程大学 一种中轨道导航卫星轨道长期预报方法
CN112036037B (zh) * 2020-08-31 2022-09-02 北京理工大学 一种倾斜地球同步轨道的长期演化快速分析方法
CN112298614B (zh) * 2020-09-18 2021-09-14 中国人民解放军战略支援部队航天工程大学 一种推力在轨标定试验方法
CN112278330B (zh) * 2020-09-27 2022-02-01 北京控制工程研究所 一种基于星时驱动的电推进位置保持方法
CN113218660B (zh) * 2021-06-14 2022-09-13 中国西安卫星测控中心 一种电推力矢量在轨标定方法
CN114771873B (zh) * 2022-03-24 2024-05-03 北京控制工程研究所 一种超低轨卫星轨道自主精确维持方法

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103472820A (zh) * 2013-09-18 2013-12-25 哈尔滨工业大学 一种基于偏最小二乘算法的推进系统故障诊断方法
CN104898645A (zh) * 2015-04-30 2015-09-09 北京控制工程研究所 一种卫星故障检测隔离恢复策略及策略动态调整方法
CN105353621A (zh) * 2015-11-30 2016-02-24 北京控制工程研究所 一种地球静止轨道卫星电推力器故障模式推力分配方法
CN108490963A (zh) * 2018-02-08 2018-09-04 中国空间技术研究院 全电推进卫星电推力器故障模式下的位置保持方法及系统

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103472820A (zh) * 2013-09-18 2013-12-25 哈尔滨工业大学 一种基于偏最小二乘算法的推进系统故障诊断方法
CN104898645A (zh) * 2015-04-30 2015-09-09 北京控制工程研究所 一种卫星故障检测隔离恢复策略及策略动态调整方法
CN105353621A (zh) * 2015-11-30 2016-02-24 北京控制工程研究所 一种地球静止轨道卫星电推力器故障模式推力分配方法
CN108490963A (zh) * 2018-02-08 2018-09-04 中国空间技术研究院 全电推进卫星电推力器故障模式下的位置保持方法及系统

Non-Patent Citations (4)

* Cited by examiner, † Cited by third party
Title
Geostationary satellites autonomous closed loop station keeping;Mauricio M.Guelman;《Acta Astronautica》;20140531;全文 *
Geostationary station-keeping with electric propulsion in full and failure modes;Lincheng Li 等;《Acta Astronautica》;20190314;全文 *
地球静止卫星电推进轨道保持策略优化;蒯政中 等;《中国空间科学技术》;20180625;第38卷(第3期);全文 *
基于GEO卫星的小推力推进器构型设计与轨道转移设计研究;邵珠军;《中国优秀硕士学位论文全文数据库》;20170331;全文 *

Also Published As

Publication number Publication date
CN109063380A (zh) 2018-12-21

Similar Documents

Publication Publication Date Title
CN109063380B (zh) 一种静止轨道电推进卫星故障检测方法及位置保持方法
Subbarao et al. Nonlinear control of motion synchronization for satellite proximity operations
Abdelrahman et al. Sigma-point Kalman filtering for spacecraft attitude and rate estimation using magnetometer measurements
Bak et al. Autonomous attitude determination and control system for the Orsted satellite
Li et al. Geostationary station-keeping with electric propulsion in full and failure modes
Malyshev et al. Orbital corrections of space vehicles while performing dynamic operations
Varma Control of satellites using environmental forces: aerodynamic drag/solar radiation pressure
Yoo et al. Attitude control system design & verification for CNUSAIL-1 with solar/drag sail
Marchetti et al. Electric propulsion and controller design for drag-free spacecraft operation
Ghisi et al. Drag-free attitude and orbit control system performance of ESA’s GOCE mission during low orbit operations and de-orbiting
Allasio et al. GOCE mission: design phases and in-flight experiences
Ceriotti et al. Hybrid solar sail and SEP propulsion for novel Earth observation missions
Boskovic et al. A globally stable scheme for spacecraft control in the presence of sensor bias
Adnane et al. Reliable Kalman filtering for satellite attitude estimation under gyroscope partial failure
Moorthy et al. Extended lifetime of CubeSats in the lower thermosphere with active attitude control
Romanazzo et al. In-orbit experience with the drag-free attitude and orbit control system of ESA’s gravity Mission GOCE
Rigatos et al. A nonlinear optimal control method for attitude stabilization of micro-satellites
Imani et al. Optimal sliding mode control for spacecraft formation flying
Xing et al. An efficient momentum dumping method through an alternative sun pointing strategy for small Near Equatorial Orbit satellite
Henry et al. Sliding-Mode Control for On-Orbit Rendezvous with a Fleeing Passive Target on a Circular Capture Trajectory
Moorthy et al. Extended Orbital Flight of a CubeSat in the Lower Thermosphere with Active Attitude Control
Harinath et al. A Novel Technique for Reference Attitude Generation in Inclined Orbit Constellation
Mashtakov et al. Study of the accuracy provided by small satellite attitude determination & control system
Robertson Control system design for spacecraft formation flying: theory and experiment
Hamzah et al. Unmodeled disturbances torque exerted on RazakSAT’s attitude during sun tracking mode

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