技术领域
[0001] 本
发明涉及一种静止轨道电推进卫星故障检测方法及位置保持方法,尤其涉及静止轨道电推进卫星的推
力故障检测、位置保持等重要技术，属于
航天器轨道动力学与控制领域。
背景技术
[0002] 静止轨道(Geostationary Orbit,GEO)指轨道高度为35786km、倾
角为零的特殊地球同步轨道。在静止轨道上运行的航天器，通常为高成本、大功率、长寿命、高
精度的大型复杂航天器，具有
覆盖面积大，且相对于地面静止的特点，在通信、导航、预警、遥感、气象、数据中继等军用和民用领域发挥着重要的作用，经济价值非常之高。
[0003] GEO卫星的轨道运动主要受地球非球形摄动、日月引力摄动、太阳光压摄动三种主要摄动的影响。在无控情况下，GEO卫星会在摄动力作用下逐渐漂出定点窗口，一旦离开定点窗口，GEO卫星就会失效成为空间碎片，不仅无法实现自身正常的通信导航等功能，还会危及附近经度的其他卫星。GEO卫星的位置保持技术是指在特定的时刻实施一定的轨道机动使卫星保持在定点经度附近的一类控制方法，该方法能够有效抑制高轨摄动因素引起的漂移现象，提升GEO卫星的在轨运行
稳定性。
[0004] 目前GEO卫星的位置保持方法，按推进
能源的不同分为两类，一类是基于化学推进的位置保持，另一类是基于电推进的位置保持。电推进系统最显著的优点是高比冲，即等量
燃料产生的推力更大，能够有效节省燃料，从而减少发射成本、增加有效
载荷的
质量。因此，本发明针对的对象是搭载电推进系统的静止轨道卫星。
[0005] 根据美国国家航空航天局(NASA)提供的电推力器测试报告，以及由佐治亚理工大学Saleh教授团队提供的近15年电推进任务统计报告，电推进系统在卫星在轨服务期间几乎无法避免故障情况。推进系统的故障将直接导致航天器部分或完全失效，不仅无法完成既定空间任务，还会威胁到空间环境的安全。因此，故障检测系统是保证卫星正常在轨运行的必备环节。
[0006] 故障模式下的位置保持策略包含两个部分：第一，检测推力器的故障状态，
跟踪推力器
输出推力的变化；第二，关闭故障推力器，并消除推力器故障引起的轨道漂移。当前文献仅针对完备模式下的位置保持策略或是仅针对故障状态已知的位置保持策略进行了研究。文献(M.Guelman,Geostationary satellites autonomous closed loop station keeping,Acta Astronautica,97,2014,pp.9–15.)研究了以当前时刻轨道信息为输入的闭环位置保持策略，但前提是推力器正常工作。文献(J.Zhang,S.Zhao,Z.Zhou,X.Li,Optimal Station Keeping for XIP Thrusters in Failure Mode under Eclipse Constraints,Journal of Aerospace Engineering,29(6),2016,pp.1-11.)研究了考虑日蚀情况下的最优燃料位置保持策略，但假设推力器故障状态已知。目前的位置保持策略没有考虑推力器故障检测环节，更缺乏一套具有故障检测和漂移消除两大功能的故障模式位置保持方案。因此本发明旨在提出一种静止轨道电推进卫星故障检测方法及位置保持方法，该方法不仅能精确检测推力器的故障状态，而且在有效识别故障状态之后，能够快速消除由推力器故障引起的轨道漂移，从而实现故障模式下的位置保持。
发明内容
[0007] 本发明公开的一种静止轨道电推进卫星故障检测方法要解决的技术问题是：提供一种静止轨道电推进卫星故障检测方法，实现静止轨道电推进卫星故障检测，能够有效跟踪任意形式的推力变化，普遍适用于高轨及地球静止轨道电推进卫星。
[0008] 进一步的，在本发明公开的一种静止轨道电推进卫星故障检测方法
基础上，本发明还公开的一种静止轨道电推进卫星故障模式的位置保持方法要解决的技术问题是：提供一种静止轨道电推进卫星故障模式的位置保持方法，实现静止轨道电推进卫星在故障模式下的位置保持，能够消除任意形式的推力变化引起的轨道漂移，普遍适用于高轨及地球静止轨道电推进卫星。
[0009] 本发明的目的通过以下技术方案实现。
[0010] 本发明公开的一种静止轨道电推进卫星故障检测方法，通过Lagrange行星运动方程建立静止轨道卫星受环境摄动下的瞬时轨道运动模型；通过平均法建立静止轨道卫星受环境摄动下的平均轨道运动模型；通过Gauss摄动方程和推力器布局建立静止轨道卫星受控制力作用下的轨道运动模型；推导静止轨道卫星受所有力作用下的轨道运动模型；设计故障检测器，并验证所设计故障检测器的收敛性；
[0011] 本发明还公开的一种静止轨道电推进卫星故障模式的位置保持方法，在所述的一种静止轨道电推进卫星故障检测方法基础上，还包括如下步骤：通过滤波方法计算推力器故障导致的轨道漂移；设计控制律，通过控制律修正推力器故障导致的轨道漂移；离散控制律，获得电推力器的
开关序列；通过采用开关序列，实现静止轨道电推进卫星故障模式下的位置保持。
[0012] 本发明公开的一种静止轨道电推进卫星故障检测方法，包括如下步骤：
[0013] 步骤一：通过Lagrange行星运动方程建立静止轨道卫星受环境摄动下的瞬时轨道运动模型，所述的环境摄动包括地球非球形摄动RENP、日月引力摄动RLSP、太阳光压摄动RSRP。
[0014] 步骤1.1：建立无奇点轨道要素；
[0015] 相对于经典轨道要素{a,e,i,Ω,ω,M}，无奇点轨道要素{a,ex,ey,ix,iy,λ}，能避免轨道倾角i和偏心率e极小情况下的奇异问题。无奇点轨道要素和经典轨道要素的映射关系如下所示：
[0016]
[0017] 其中，a是半场轴，e是偏心率，i是轨道倾角，ω是近地点角距，Ω是升交点赤经，M是平近点角，θG是格林尼治
恒星时角。
[0018] 步骤1.2：建立Lagrange行星运动方程；
[0019] 静止轨道电推进卫星的轨道运动主要受三种环境摄动影响。三种自然摄动包括地球非球形摄动RENP、日月引力摄动RLSP和太阳光压摄动RSRP。
[0020] Lagrange行星运动方程用于描述静止轨道卫星在三种自然摄动作用下的轨道运动。采用Lagrange行星运动方程并略去二阶小量后的轨道运动模型如下所示：
[0021]
[0022] 其中，n是卫星轨道运动的平均角速率，R是摄动势函数，Ωe是地球自转角速率。
[0023] 步骤1.3：建立地球非球形摄动势函数；
[0024] 静止轨道非球形摄动(Earth’s Non-spherical Perturbation,ENP)势函数具有如下形式：
[0025]
[0026] 其中，μ是地球引力常数，r是航天器与地心的距离，Re是地球半径，Jn是带谐项系数，Jnm是田谐项系数，Pn(φ),Pnm(φ)是勒让德多项式，φ是地心纬度，λnm是引力场系数。
[0027] Pn(φ),Pnm(φ)的表达式如公式(4)所示，Pn(φ)视为Pn0(φ)
[0028]
[0029] 作为优选，步骤1.3所述的地球非球形摄动势函数优选仅包含前四阶带谐项J2,J3,J4和前三阶田谐项J22,J31,J33的摄动势函数，所述摄动势函数具体推导方法如下：
[0030] 包含前四阶带谐项J2,J3,J4和前三阶田谐项J22,J31,J33的摄动势函数如公式(5)所示：
[0031]
[0032]
[0033] 因为静止轨道的轨道倾角较小，所以地心纬度φ满足如下条件：
[0034]
[0035] 将公式(7)带入公式(5)和公式(6)，获得考虑静止轨道特点的摄动势函数如公式(8)所示：
[0036]
[0037] 步骤1.4：建立日月引力摄动势函数；
[0038] 日月引力摄动(Luni-Solar Perturbation,LSP)的势函数具有如下形式：
[0039]
[0040] 下标k代表第三引力体，k＝s,m，s代表太阳，m代表月亮；μk是第三引力体的引力常数，rk是第三引力体与地心的距离，θk是rk和r的夹角。
[0041] 步骤1.5：建立太阳光压摄动势函数；
[0042] 太阳光压摄动(Solar Radiation Perturbation,SRP)的势函数具有如下形式：
[0043]
[0044] CR是光压系数，S是卫星正对太阳的面积，m是卫星质量，P是单位面积的
太阳辐射压，x,y,z是卫星相对地心的位置矢量在地心惯性系下的三个分量，xs,ys,zs是太阳相对地心的位置矢量在地心惯性系下的三个分量。
[0045] 步骤1.6：将公式(3)、(9)、(10)分别带入公式(2)推导三种摄动作用下的瞬时轨道要素变化率：
[0046] 由于公式(3)带入公式(10)后的表达式非常复杂，仅给出如公式(8)所示的优选的非球形摄动势函数带入公式(10)的结果，即非球形摄动(ENP)引起的瞬时轨道要素变化率如下：
[0047]
[0048] 其中，ac是平半长轴，取值与静止轨道卫星的定点经度有关，nc是ac对应的轨道运动平均角速率。
[0049]
[0050] 日月引力摄动(LSP)引起的瞬时轨道要素变化率如下：
[0051]
[0052] 其中，xk,yk,zk是第三引力体在地心惯性坐标下的三个分量。
[0053] 太阳光压摄动(SRP)引起的瞬时轨道要素变化率如下：
[0054]
[0055] 步骤1.7：联立公式(11)、(13)、(14)，建立静止轨道卫星受环境摄动下的瞬时轨道运动模型，如公式(15)所示：
[0056]
[0057] 其中，表示瞬时轨道要素变化率，
[0058] 步骤二：在步骤一建立的静止轨道卫星受环境摄动下的瞬时轨道运动模型基础上，通过平均法计算三种环境摄动引起的平均轨道要素和短周期项，建立静止轨道卫星受环境摄动下的平均轨道运动模型。
[0059] 步骤2.1：分别推导三种摄动作用下的平均轨道要素和短周期项；
[0060] 分析公式(11)、(13)、(14)各项的周期性，定义周期大于一个静止轨道轨道周期的项为长周期项，周期小于一个静止轨道轨道周期的项为短周期项；并对公式(11)、(13)、(14)按照公式(16)求平均，即能够推导出三种摄动作用下平均轨道要素变化率按照公式(17)能够推导三种摄动作用下短周期项变化率
[0061]
[0062]
[0063] 积分短周期项变化率 能够获得短周期项(σsp)ENP,(σsp)LSP,(σsp)SRP。
[0064] 三种环境摄动作用下的平均轨道要素变化率和短周期项的表达式如下：
[0065] 非球形摄动(ENP)引起的平均轨道要素变化率和短周期项如下：
[0066]
[0067]
[0068] 其中， 是半长轴和平经度的短周期修正项；
[0069]
[0070] 日月引力摄动(LSP)引起的平均轨道要素变化率和短周期项如下：
[0071]
[0072]
[0073] 其中，CΘ-LSP＝cos(ωk+Mk)cosΩk，SΘ-LSP＝sin(ωk+Mk)sinΩkcosik。
[0074] 太阳光压摄动(SRP)引起的平均轨道要素变化率和短周期项如下：
[0075]
[0076]
[0077] 步骤2.2：联立公式(18)-(23)，获得静止轨道卫星受环境摄动下的平均轨道运动模型如公式(24)。
[0078]
[0079] 步骤三：通过Gauss摄动方程和推力器布局建立静止轨道卫星受控制力作用下的轨道运动模型，并将模型线性化；
[0080] 步骤3.1：建立Gauss摄动方程；
[0081] 除三种环境摄动外，静止轨道卫星的轨道运动还受到推力器输出推力的英雄。Gauss摄动方程用于描述静止轨道卫星在电推进系统提供的小推力作用下的轨道运动。采用Gauss摄动方并略去二阶小量后的轨道运动模型如下所示：
[0082]
[0083] T是推力器输出的推力大小；dR,dT,dN是推力器在卫星质心轨道
坐标系下沿径向、切向和法向的推力分量与T的比值；u是电推力器开关常数，u＝0or 1，0表示推力器关闭，1表示推力器开启；m是静止轨道卫星的总质量；V是静止轨道卫星的速度；l是静止轨道卫星的平赤经，l＝Ω+ω+M。
[0084] 步骤3.2：将推力器按推力器布局投影到卫星质心轨道坐标系并带入到Gauss摄动方程。
[0085] 步骤3.3：建立静止轨道卫星受控制力作用下的轨道运动模型，如公式(26)所示：
[0086]
[0087] 其中，Ti表示第i个推力器输出的推力大小，ui是第i个推力器的开关常数，分别表示第i个推力器在卫星质心轨道坐标系下沿径向、切向和法向的推力分量与Ti的比值。
[0088] 步骤3.4：线性化建立的静止轨道卫星受控制力作用下的轨道运动模型；
[0089] 将公式(26)写成矩阵形式，如公式(27)所示：
[0090]
[0091] 其中fEP(σ)是关于当前轨道六要素的非线性函数矩阵，具有6×3维；u是3×1维的常值控制向量，表征推力大小在质心轨道坐标系径、切、法向的分量，
[0092]
[0093] 将公式(27)在期望的卫星定点位置σd按泰勒级数展开并保留一阶项，如公式(28)所示：
[0094]
[0095] 其中，因为 和 是仅取决于期望定点位置的常值矩阵，因此B(σ)是关于σ的线性函数矩阵，公式(28)是静止轨道卫星受控制力作用下的线性化轨道运动模型。
[0096] 步骤四：基于步骤二建立的静止轨道卫星受环境摄动下的平均轨道运动模型和步骤三建立的静止轨道卫星受控制力作用下的轨道运动模型，推导静止轨道卫星受所有力作用下的轨道运动模型。
[0097] 联立步骤2.2的公式(24)和步骤3.4的公式(28)，获得静止轨道卫星受所有力作用下的轨道运动模型如公式所示。
[0098]
[0099] 其中，B(σ)是关于σ的线性函数矩阵，w和v分别表示表示系统噪声和量测噪声。
[0100] 步骤五：设计故障检测器，并验证所设计故障检测器的收敛性；
[0101] 步骤5.1：设计故障检测器，通过故障检测器跟踪推力器故障引起的推力变化。
[0102] 推力器故障会引起静止轨道卫星轨道漂移，如果不能及时检测故障状态并且消除故障引起的漂移，一旦漂出定点窗口，静止轨道卫星将无法维持正常的通信、导航功能。导致推力器故障主要有两方面的原因，一是自身系统老化；二是外界因素的影响，包括非合作目标的碰撞、干扰。电推力器在理想工作状态输出常值推力，因此，其故障形式主要表现为实际输出推力低于额定推力。
[0103] 当推力器处于故障状态时，实际输出推力由已知变为未知。离散并反解如公式(29)所示的
状态方程能够获得实际推力u与状态量 的映射关系，其中 属于超前状态信息。
[0104]
[0105] 引入辅助变量ξ用于观测故障状态下的实际推力u。
[0106]
[0107] 引入线性观测器用于描述估计推力 的递推关系。
[0108]
[0109] 设计系数矩阵Γ和Ψ实现估计推力 对实际推力u的有效跟踪，设计的参数如下：
[0110]
[0111] 将公式(31)和公式(32)带入公式(30)，获得辅助变量ξ的递推表达式，如公式(34)所示。
[0112]
[0113] 当且仅当Γ＝ΨB-1，辅助变量ξ的递推只取决于前一时刻的轨道信息，不需要对未来轨道信息进行预测，如公式(35)所示。
[0114]
[0115] 公式(31)和公式(35)共同组成所设计故障检测器的动力学方程，通过故障检测器能够跟踪推力器故障引起的推力变化。
[0116] 步骤5.2：验证所设计故障检测器的收敛性；
[0117] 定义如公式(36)所示的推力估计误差：
[0118]
[0119] 将公式(31)和公式(35)带入公式(36)：
[0120]
[0121] 在公式(37)右端增加uk构造ek，获得推力估计误差的递推表达式：
[0122]
[0123] 当k趋于无穷时，ek趋于零。因此，推力估计误差是渐进稳定的，且最终收敛于零。当Ψ＜I，即λi＜1,i＝1,2,3时，所设计的故障检测器具有良好的局部性能，即能够验证所设计的故障检测器能够使推力估计误差渐进稳定，且最终收敛于零。
[0124] 在所述的一种静止轨道电推进卫星故障检测方法基础上，本发明还公开的一种静止轨道电推进卫星故障模式的位置保持方法，包括所述的一种静止轨道电推进卫星故障检测方法的步骤一至五，还包括步骤六至步骤九。
[0125] 步骤六：通过滤波方法计算推力器故障导致的轨道漂移；
[0126] 步骤五已验证估计推力 能够精确跟踪故障状态下的实际推力u，因此，估计推力能够视为已知输入带入如公式(29)所示的静止轨道卫星的系统方程。针对输入已知的非线性系统，应用滤波方法即能够估计故障状态下的平均轨道要素 当平均轨道要素估计值 与平均轨道要素实际值 的误差趋近于零时，当前时刻q的平均轨道要素估计值 与理想平均轨道要素 之间的偏差即为推力故障持续期间引起的轨道漂移。
[0127] 步骤六所述的滤波方法包括无迹卡尔曼滤波方法、
扩展卡尔曼滤波。
[0128] 作为优选，步骤六所述的滤波方法优选无迹卡尔曼滤波方法，针对形如公式(39)的系统方程，无迹卡尔曼滤波
算法(UKF)的流程如下：
[0129]
[0131] 根据UT(unscented transform)变换选取采样点：
[0132]
[0133] 其中，Sx是Cholesky因子，初值为 函数chol(·)是Cholesky分解函数，Zi是第i个采样点对应的向量空间。
[0134] 流程2：状态更新。
[0135]
[0136]
[0137] 其中，函数qr(·)为QR分解函数，函数cholupdate(·)为Cholesky分解一次更新。为协方差加权的权值， 为均值加权的权值。
[0138] 流程3：量测更新：
[0139]
[0140] 其中，操作符A/B为矩阵相除，通过求解方程AATx＝ATB获得。
[0141] 步骤七：设计控制律，通过控制律修正推力器故障导致的轨道漂移。
[0142] 步骤七所述控制律包括开环控制律和闭环控制律，开环控制律基于最优控制理论，闭环控制律基于反馈信息。
[0143] 作为优选，步骤七所述控制律优选反馈控制律，步骤七具体实现方法如下：
[0144] 将步骤六获得的轨道漂移作为初始偏差，使用如下形式控制律修正偏差：
[0145]
[0146] 是期望轨道面沿地心轨道坐标系法向的单位向量， 是当前轨道面沿地心轨道坐标系径向的单位向量， 是当前时刻航天器的速度矢量。c1,c2是增益，ac＝as+ka(a-ad)，as是标准静止轨道的半长轴，ka是设计的常值系数，ad是期望半长轴。
[0147] 先给出由轨道要素描述的 在地心惯性系(ECI)的投影。
[0148]
[0149]
[0150]
[0151] 其中，
[0152]
[0153] 通过对公式(48)进行周期性分析，将包含sinkl,coskl且k≥2的项视为短周期项，去除短周期项后的表达式如下：
[0154]
[0155] 通过对公式(47)进行周期性分析，将包含sin2M,cos2M的项视为短周期项，去除短周期项后的表达式如下：
[0156]
[0157] 其中，u＝Ω+ω，函数ysg＝fsg(xsg),ycg＝fcg(xcg)表示(xsg,ysg)和(xcg,ycg)具有相同的周期性。
[0158] 将公式(47)、(49)、(50)带入公式(44)，获得控制律在地心惯性系下的表达式，如公式(51)所示。
[0159]
[0160] 将公式(51)所示的 投影到航天器质心轨道坐标系，即得到能够直接带入如公式(29)所示系统方程的控制律。
[0161]
[0162] 其中，
[0163] 步骤八：离散步骤七所述的控制律，获得电推力器的开关序列；
[0164] 公式(52)给出了连续型的控制律，但由于电推力器输出常值推力且具有开关特性，因此需要根据推力器的性能参数将连续型的控制律离散化。
[0165] 通过计算推力器常值推力产生的
加速度arate，设计以arate为
阈值的函数，如公式(53)：
[0166]
[0167] 联立公式(52)和公式(53)，即能够获得推力器的开关序列。
[0168] 步骤九：通过采用步骤八所述的开关序列，能够消除任意形式的推力变化引起的轨道漂移，实现静止轨道电推进卫星故障模式下的位置保持。
[0169] 有益效果
[0170] 1、本发明公开的一种静止轨道电推进卫星故障检测方法，所有运行在高轨及静止轨道的航天器须考虑三种环境摄动，推导的平均轨道运动模型考虑了三种主要环境摄动的影响，且推导的公式不受摄动势函数阶数和次数的限制，因此适用于所有运行在高轨及静止轨道的航天器，通过该方法对静止轨道卫星进行控制，为静止轨道卫星在轨运行的稳定性和安全性提供了理论依据和技术
支撑。
[0171] 2、本发明公开的一种静止轨道电推进卫星故障检测方法，通过设计辅助变量和观测器能够精确估计故障模式下的推力，且估计的推力只取决于前一时刻的轨道信息，不需要对未来轨道信息进行预测，能够简化模型、提升计算效率。
[0172] 3、本发明公开的一种静止轨道电推进卫星故障检测方法，通过设计辅助变量和观测器能够精确估计故障模式下的推力，进一步的，通过合理选取检测器的参数，能够有效跟踪任意形式的推力变化，从而有效应对卫星在轨运行过程中的推力故障情况。
[0173] 4、本发明还公开的一种静止轨道电推进卫星故障模式的位置保持方法，在所述的一种静止轨道电推进卫星故障检测方法基础上实现，因此，具有所述的一种静止轨道电推进卫星故障检测方法的优点，此外，还包括步骤六至步骤九；通过滤波方法计算推力器故障导致的轨道漂移；设计控制律，通过控制律修正推力器故障导致的轨道漂移；离散所设计的控制律，获得电推力器的开关序列；通过采用所设计的开关序列，能够消除任意形式的推力变化引起的轨道漂移，从而实现静止轨道电推进卫星故障模式下的位置保持。
[0174] 5、本发明还公开的一种静止轨道电推进卫星故障模式的位置保持方法，通过合理选取采样点并应用滤波方法，和建立的动力学模型具有良好的适配性，能够有效估计故障模式下的轨道漂移。
[0175] 6、本发明还公开的一种静止轨道电推进卫星故障模式的位置保持方法，通过设计控制律并对其进行离散处理，能够获得电推进系统的推力开关序列，普遍适用于所有搭载电推进系统的静止轨道卫星的位置保持任务。
[0176] 7、本发明还公开的一种静止轨道电推进卫星故障模式的位置保持方法，通过设计反馈控制律，并合理选取
控制器参数，能够消除任意形式的推力变化引起的轨道漂移，相比正常模式下的位置保持策略精度有所提高。
附图说明
[0177] 图1为一种静止轨道电推进卫星故障检测方法的
流程图；
[0178] 图2为一种静止轨道电推进卫星故障模式的位置保持方法的流程图；
[0179] 图3为故障模式下实现位置保持的推力器布局；
[0180] 图4为所设计故障检测器对故障推力的跟踪情况；
[0181] 图5为故障模式下滤波算法估计的平均轨道要素与实际轨道要素的偏差变化，其中：图5(a)为半长轴偏差的变化、图5(b)为平经度偏差的变化、图5(c)为偏心率矢量x分量偏差的变化、图5(d)为偏心率矢量y分量偏差的变化、图5(e)为倾角矢量x分量偏差的变化、图5(f)为倾角矢量y分量偏差的变化；
[0182] 图6为所设计故障模式位置保持控制律在切向和法向提供的控制加速度；
[0183] 图7为所设计故障模式位置保持控制律在切向提供的控制加速度及其对应的推力器开关序列所提供的控制加速度；
[0184] 图8为无控、理想、受控三种工况下的轨道要素变化，其中：图8(a)为半长轴变化、图8(b)为平经度变化、图8(c)为偏心率矢量x分量的变化、图8(d)为偏心率矢量y分量的变化、图8(e)为倾角矢量x分量的变化、图8(f)为倾角矢量y分量的变化。
具体实施方式
[0185] 为了更好地说明本发明的目的和优点，下面结合附图对本发明的具体实施方式作进一步详细的说明。
[0187] 如图1所示，本实施例公开的一种静止轨道电推进卫星故障检测方法，为验证该方法，首先，选取运行在地球静止轨道上的卫星作为主要研究对象。仿真的起始时刻为2020年1月1日0时0分0秒(世界时)，卫星的基本参数如下表所示。
[0188] 表1卫星参数
[0189]
[0190]
[0191] 步骤一：通过Lagrange行星运动方程建立静止轨道卫星受环境摄动下的瞬时轨道运动模型，所述的环境摄动包括地球非球形摄动RENP、日月引力摄动RLSP、太阳光压摄动RSRP。；
[0192] 经典轨道要素在描述小倾角和小偏心率的GEO卫星的轨道运动时，容易出现奇异问题。为避免奇异，引入无奇点轨道要素，如公式(1)所示。静止轨道电推进卫星所受主要环境摄动有三种，分别是地球非球形摄动RENP、日月引力摄动RLSP、太阳光压摄动RSRP。用略去二阶小量的Lagrange行星运动方程如公式(2)所示，描述环境摄动力和控制力作用下的轨道运动模型。
[0193] 用于计算地球非球形摄动RENP的地球引力场系数如下表所示：
[0194] 表2地球引力场主项系数
[0195]
[0196] 用于计算日月引力摄动RLSP、太阳光压摄动RSRP的太阳、月亮平均轨道运动模型如下。
[0197] 太阳平均轨道运动模型如公式(54)所示：
[0198]
[0199] 月亮平均轨道运动模型如公式(55)所示：
[0200]
[0201] 其中，T是当前时刻对应的儒略日，
[0202] 将三种环境摄动的势函数带入Lagrange行星运动方程，即能够建立静止轨道卫星受环境摄动下的瞬时轨道运动模型，如公式(15)所示。
[0203] 步骤二：在步骤一建立的静止轨道卫星受环境摄动下的瞬时轨道运动模型基础上，通过平均法计算三种环境摄动引起的平均轨道要素和短周期项，建立静止轨道卫星受环境摄动下的平均轨道运动模型。
[0204] 建立的静止轨道卫星受环境摄动下的平均轨道运动模型如公式(24)所示。因为推导的平均轨道运动模型考虑了三种主要环境摄动的影响，且推导的公式不受摄动势函数阶数和次数的限制，因此适用于所有运行在高轨及静止轨道的航天器，
[0205] 步骤三：通过Gauss摄动方程和推力器布局建立静止轨道卫星受控制力作用下的轨道运动模型，并将模型线性化；
[0206] 实例中的卫星搭载电推进系统，用于实现故障模式下的位置保持。推力器布局如图3所示。引入卫星质心轨道坐标系RTN用于描述推力器布局， 轴由卫星质心指向地心，轴与 轴相互垂直，指向卫星速度的正方向， 轴垂直于轨道面，方向满足右手定则。四个推力器正装在 轴和 轴的正负方向上，沿顺
时针方向依次记作N(a)、E(b)、S(c)、W(d)。
[0207] 推力器参数如下表所示：
[0208] 表3推力器参数
[0209]
[0210] 建立的静止轨道卫星受控制力作用下的线性化轨道运动模型如公式(28)。
[0211] 步骤四：基于步骤二建立的静止轨道卫星受环境摄动下的平均轨道运动模型和步骤三建立的静止轨道卫星受控制力作用下的轨道运动模型，推导静止轨道卫星受所有力作用下的轨道运动模型。
[0212] 建立的静止轨道卫星受所有力作用下的轨道运动模型如公式(29)所示。
[0213] 本实例设置系统噪声w为高斯白噪声，方差为Q；量测噪声v为高斯白噪声，方差为R。
[0214] Q＝diag(1×10-8,1×10-16,1×10-16,1×10-16,1×10-16,1×10-14) (56)[0215] 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)
[0216] 步骤五：设计故障检测器，并验证所设计故障检测器的收敛性；
[0217] 所设计故障检测器由公式(31)和公式(35)共同组成。
[0218] 对比公式(30)和公式(35)，实际推力u的计算不仅依赖当前时刻的轨道信息 还依赖未来的轨道信息 但估计的推力只取决于前一时刻的轨道信息，不需要对未来轨道信息进行预测，因此所设计的故障检测不仅能够精确跟踪故障状态下的实际推力u，还能够简化模型、提升计算效率。
[0219] 本实例采用的故障检测器的参数如下：
[0220] Ψ＝diag[λ,λ,λ] (58)
[0221] λ分别取λ＝1.35，λ＝0.85和λ＝0.55做了三组对比仿真，仿真结果如图4所示。
[0222] 图4为所设计故障检测器对故障推力的跟踪情况，红色实线表示实际推力，在第100s和第200s由于故障分别衰减常值推力的20％，在第300s时推力器关闭。绿色虚线、蓝色实线、黑色实线分别代表λ＝1.35，λ＝0.85和λ＝0.55时，所设计故障检测器估计的推力。通过对比能够发现，故障引起的推力变化在10s内即能够实现精确跟踪，且λ＜1时具有良好的局部性能，且越趋近1时跟踪精度越高。本实例中模拟的推力故障情况是间隔一段时间的两次推力衰减且推力衰减幅度递增，但所有可能出现的故障情况都能够拆分成单次推力变化的组合。因此，通过合理选取故障检测器的参数，能够有效跟踪任意形式的推力变化，从而有效应对卫星在轨运行过程中的推力故障情况。
[0223] 实施例2：
[0224] 如图2所示，本实施例在所述的一种静止轨道电推进卫星故障检测方法基础上，还公开的一种静止轨道电推进卫星故障模式的位置保持方法，包括所述的一种静止轨道电推进卫星故障检测方法的步骤一至五，还包括步骤六至步骤九，为验证该方法，前五步选取的参数与实施例1一致，步骤六至步骤九如下：
[0225] 步骤六：通过滤波方法计算推力器故障导致的轨道漂移；
[0226] 作为优选，步骤六所述的滤波方法优选无迹卡尔曼滤波方法，步骤具体实现方法如公式(39)-(43)所示。
[0227] 本实例设定状态协方差矩阵的初值如公式(59)所示：
[0228] P＝diag(1×103,1×10-5,1×10-5,2×10-3,2×10-3,3×10-2) (59)
[0229] 将步骤五所设计的故障检测器估计的推力 带入如公式(29)的系统方程，并应用上述的无迹卡尔曼滤波方法，估计推力故障导致的轨道漂移。仿真结果如图5所示。
[0230] 图5表示采用无迹卡尔曼律方法估计的轨道要素与实际轨道要素的偏差。红色实线表示推力器正常工作时的偏差变化，偏差不恒等于零是因为系统误差和观测误差始终存在。蓝色实线表示单次推力故障引起的偏差变化，故障情况如下：推力器自第0s开始正常工作，在第100s时发生推力器故障，故障导致输出推力衰减至峰值的80％，从故障发生到主动关闭故障推力器总共持续100s。绿色实线表示两次推力故障引起的偏差变化，故障情况如下：推力器自第0s开始正常工作，在第100s时发生第一次推力器故障，故障导致推力器输出衰减至峰值的80％，在第200s发生第二次推力器故障，故障导致推力器输出衰减至峰值的60％，从故障发生到主动关闭故障推力器总共持续200s。从检测到故障出现到主动关闭故障推力器的总时长定义为故障持续时间，本例中分明设定100s和200s作为故障持续时间的上限，一旦故障持续时间超过设定时间，故障推力器将会被主动关闭。蓝色实线表示的单次推力故障引起的六个轨道要素偏差均在第1852s收敛，绿色实线表示的两次推力故障引起的六个轨道要素偏差均在第2432s收敛。由此能够分别计算在第1852s和第2432s时卫星所处位置与期望定点窗口的轨道漂移。蓝色实线所表示的单次推力故障引起的轨道要素偏差在1852s收敛，此时轨道要素估计值 与平均轨道要素实际值 的偏差近似为零，推力故障持续期间引起的轨道漂移为此时平均轨道要素估计值 与理想平均轨道要素 其值如下：
[0231]
[0232] 在本实例的两种故障情况下轨道要素的估计误差都能较快且较精确地收敛，证明选用的滤波方法和建立的动力学模型具有良好的适配性，能够有效估计故障模式下的轨道漂移。
[0233] 步骤七：设计控制律，通过控制律修正推力器故障导致的轨道漂移；
[0234] 公式(60)表示推力故障持续期间引起的轨道漂移，将其值作为初始偏差，使用如公式(52)所示的反馈控制律修正偏差。在所设计的控制律中，常值系数ka＝2000。仿真结果如图6所示。
[0235] 步骤八：离散步骤七所述的控制律，获得电推力器的开关序列；
[0236] 因为电推力器的输出推力通常为常值且具有开关特性，所以需要将所设计的连续控制律离散以获得实际推力器的开关序列。根据表2提供的推力器的性能参数，又因为四个推力器分别正装在切向和法向轴的正负方向，因此选取如下加速度阈值arate＝200mN/2000kg＝1×10-4m/s2。联立公式(53)和公式(52)，即能够获得推力器的开关序列。仿真结果如图7所示。
[0237] 图6表示步骤七所设计的控制律在卫星质心轨道坐标系切向和法向提供的控制加速度。图7表示离散处理后获得的推力开关序列。图7所示矩形波的正向
信号由安装在切向轴，推力方向为正向的推力器提供；矩形波的负向信号由安装在切向轴，推力方向为负向的推力器提供；对称安装的推力器在控制周期内交替运行，这种开关特性同样适用于对称安装在法向轴的推力器上。图6所示的切向加速度和径向加速度在一个控制周期内具有重叠部分，由此表明单个周期内存在两个推力器同时开机的情况。电推进卫星的寿命周期内，为实现燃料最省的目的，工程中通常使用“单个时刻仅有一个推力器开机”的约束条件，但在故障模式下，为保证卫星始终位于定点窗口，首要目标是在最短的时间内消除故障引起的轨道漂移，推力器的开机约束能够适当放宽，因此本发明设计的故障模式开机策略是合理的。因为所设计的控制律都通过设置阈值离散为电推力器的开关序列，因此本发明提出的方法普遍适用于所有搭载电推进系统的静止轨道卫星的位置保持任务
[0238] 步骤九：通过采用步骤八所述的开关序列，能够消除任意形式的推力变化引起的轨道漂移，实现静止轨道电推进卫星故障模式下的位置保持。
[0239] 将步骤八获得的推力器开关序列带入到如公式(29)所示的系统动力学方程，即能够获得卫星的轨道要素在受控情况下的变化规律，仿真结果如图8所示。
[0240] 图8表示初始偏差一定的无控、理想、受控三种工况下的对比仿真结果。蓝色实现表示无控情况的轨道要素变化情况，黑色实现表示理想静止轨道卫星的轨道要素变化情况，绿色实线表示采用本发明设计的开关序列后的轨道要素变化情况。在收敛时间方面，六个轨道要素均能在2天(两个控制周期)内收敛。偏心率矢量(ex,ey)在当前算例中的收敛时间为0.5天，倾角矢量(ix,iy)在当前算例中的收敛时间为1天，具有良好的性能。在收敛精度方面，平半长轴 收敛至100米内，平经度 收敛至10-3度，平偏心率矢量 收敛至10-5 -5，平倾角矢量 收敛至10 。表给出了本发明设计的位保策略同主要参考文献的位保策略控制精度的对比。方案1代表本发明的结果；方案2代表文献(M.Guelman,Geostationary satellites autonomous closed loop station keeping,Acta Astronautica,97,2014,pp.9–15.)的结果；方案3代表文献(J.Zhang,S.Zhao,Z.Zhou,X.Li,Optimal Station Keeping for XIP Thrusters in Failure Mode under Eclipse Constraints,Journal of Aerospace Engineering,29(6),2016,pp.1-11.)的结果。因此，本实施例公开的一种静止轨道电推进卫星故障模式的位置保持方法(方案1)，相较推力器正常工作的位置保持策略(方案2)和已知一对推力器故障的位置保持策略(方案3)，控制精度均有1到2个数量级的提升，具有更优的控制性能。
[0241] 表4主流控制方案的精度对比
[0242]
[0243] 综上所述，以上仅为本发明的较佳实施例而已，并非用于限定本发明的保护范围。凡在本发明的精神和原则之内，所作的任何
修改、等同替换、改进等，均应包含在本发明的保护范围之内。
