CN113311850A - 一种抗慢变干扰的相平面姿态控制方法 - Google Patents

一种抗慢变干扰的相平面姿态控制方法 Download PDF

Info

Publication number
CN113311850A
CN113311850A CN202110450138.1A CN202110450138A CN113311850A CN 113311850 A CN113311850 A CN 113311850A CN 202110450138 A CN202110450138 A CN 202110450138A CN 113311850 A CN113311850 A CN 113311850A
Authority
CN
China
Prior art keywords
axis
straight line
region
thruster
moment
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.)
Granted
Application number
CN202110450138.1A
Other languages
English (en)
Other versions
CN113311850B (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 Control Engineering
Original Assignee
Beijing Institute of Control Engineering
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 Control Engineering filed Critical Beijing Institute of Control Engineering
Priority to CN202110450138.1A priority Critical patent/CN113311850B/zh
Publication of CN113311850A publication Critical patent/CN113311850A/zh
Application granted granted Critical
Publication of CN113311850B publication Critical patent/CN113311850B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • 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/08—Control of attitude, i.e. control of roll, pitch, or yaw
    • G05D1/0808—Control of attitude, i.e. control of roll, pitch, or yaw specially adapted for aircraft
    • G05D1/0816—Control of attitude, i.e. control of roll, pitch, or yaw specially adapted for aircraft to ensure stability
    • G05D1/0833—Control of attitude, i.e. control of roll, pitch, or yaw specially adapted for aircraft to ensure stability using limited authority control

Landscapes

  • Engineering & Computer Science (AREA)
  • Computer Security & Cryptography (AREA)
  • Aviation & Aerospace Engineering (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Remote Sensing (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Automation & Control Theory (AREA)
  • Control Of Position, Course, Altitude, Or Attitude Of Moving Bodies (AREA)

Abstract

本发明提出一种抗慢变干扰的相平面姿态控制方法，包括建立干扰下的航天器动力学模型；针对由于航天器所受外部慢变干扰、航天器三轴动力学耦合及三轴控制输入耦合带来的总干扰角加速度，设计扩张状态观测器对其进行估计；设计三轴相平面分区；针对三轴相平面各个分区设计三轴相平面控制指令；航天器的推力器根据三轴相平面控制指令进行喷气。本发明方法实现了慢变干扰下的相平面控制补偿设计，能够有效提高系统对干扰的适应能力，兼顾闭环系统稳态性能，减少了喷气次数和燃料消耗。

Description

一种抗慢变干扰的相平面姿态控制方法
技术领域
本发明涉及一种抗慢变干扰的相平面姿态控制方法，特别涉及航天器在轨运行过程中慢变干扰的喷气补偿方法，属于姿态控制技术领域。
背景技术
航天器在轨运行中会受到多种外部慢变干扰，其中典型的慢变干扰包括：变轨过程中由于变轨发动机推力偏心和航天器的质心偏差产生的常值大干扰力矩；低轨道航天器受到的慢变气动干扰力矩；空间飞行器普遍受到的重力梯度力矩、太阳光压力矩等呈慢变特性的空间环境力矩。慢变干扰的长期积累会严重影响航天器姿态控制精度，导致稳态特性难以满足任务要求，喷气次数增多，燃料消耗增加。
发明内容
本发明的技术解决问题是：克服现有技术的不足之处，提供了一种抗慢变干扰的相平面姿态控制方法，通过对干扰的实时估计与补偿，实现了慢变干扰下的相平面控制补偿设计，能够有效提高系统对干扰的适应能力，兼顾闭环系统稳态性能，减少了喷气次数和燃料消耗。
本发明的技术解决方案是：
一种抗慢变干扰的相平面姿态控制方法，该方法所针对的航天器上配置有大推力器和小推力器，该具体步骤包括：
(1)建立慢变干扰下的航天器动力学模型；
(2)建立用于估计航天器三轴的总干扰角加速度的扩张状态观测器，并使用所建立的扩张状态观测器估计航天器三轴的总干扰角加速度；
(3)构建三轴的相平面，并将所构建的三轴相平面分别分为七个区；
(4)根据步骤(1)中所建立的慢变干扰下的航天器动力学模型、根据所述步骤(2)的扩张状态观测器估计出的三轴总干扰角加速度和步骤(3)的三轴相平面分区，得到每轴相平面的七个区的控制指令；
(5)航天器的推力器根据步骤(4)得到的三轴相平面控制指令进行喷气。
所述步骤(1)的建立的慢变干扰下的航天器动力学模型为：
Figure BDA0003038403670000021
其中，x1,x2,x3表示慢变干扰下的航天器的三轴姿态角，三轴姿态角为滚动姿态角、俯仰姿态角和偏航姿态角；
Figure BDA0003038403670000022
表示慢变干扰下的航天器的三轴姿态角速度，三轴姿态角速度为滚动姿态角速度、俯仰姿态角速度和偏航姿态角速度；
Figure BDA0003038403670000023
表示慢变干扰下的航天器的三轴姿态角加速度，三轴姿态角加速度为滚动姿态角加速度、俯仰姿态角加速度和偏航姿态角加速度；
Figure BDA0003038403670000024
Figure BDA0003038403670000025
表示航天器的转动惯量矩阵；
Figure BDA0003038403670000026
表示慢变干扰下的航天器的三轴姿态控制指令；
Figure BDA0003038403670000027
d1,d2,d3表示慢变干扰在航天器的三轴分量；
D＝[d1,d2,d3]T
Figure BDA0003038403670000031
表示由于航天器所受外部慢变干扰D、航天器三轴动力学耦合及三轴控制输入耦合带来的总干扰角加速度在三轴的分量；
Y＝[y1,y2,y3]T，其中，y1,y2,y3表示慢变干扰下的航天器的三轴姿态角的测量值；
V＝[v1,v2,v3]T，其中，v1,v2,v3表示慢变干扰下的航天器的三轴姿态角的测量噪声。
所述步骤(2)中，建立用于估计航天器滚动轴的总干扰角加速度的扩张状态观测器Σ1如下：
Figure BDA0003038403670000032
其中，z11，z21，z31为三个待求的未知函数；
Figure BDA0003038403670000033
为z11的时间导数，
Figure BDA0003038403670000034
为z2的时间导数，
Figure BDA0003038403670000035
为z3的时间导数；
β11为一常数，范围为h-1～10h-1；
β21为一常数，范围为0.01β11～0.1β11；
β31为一常数，范围为0.0001β11～0.001β11；
函数fal(e1,η,δ)＝e1s/δ1-η+|e1|ηsign(e1)(1-s)，s＝0.5(sign(e1+δ)-sign(e1-δ))，常数η∈(0,1)，常数δ>0取0.1h～10h，h表示数值离散步长，选取范围为0～0.01；
航天器滚动轴的总干扰角加速度的理论值
Figure BDA0003038403670000036
ad1的估计值
Figure BDA0003038403670000041
即求得的z31为航天器滚动轴的总干扰角加速度的估计值。
所述步骤(2)中，建立用于估计航天器俯仰轴的总干扰角加速度的扩张状态观测器Σ2如下：
Figure BDA0003038403670000042
其中，z12，z22，z32为三个待求的未知函数；
Figure BDA0003038403670000043
为z12的时间导数，
Figure BDA0003038403670000044
为z22的时间导数，
Figure BDA0003038403670000045
为z32的时间导数；
β12为一常数，范围为h-1～10h-1；
β22为一常数，范围为0.01β12～0.1β12；
β32为一常数，范围为0.0001β12～0.001β12；
函数fal(e2,η,δ)＝e2s/δ1-η+|e2|ηsign(e2)(1-s)，s＝0.5(sign(e2+δ)-sign(e2-δ))，常数η∈(0,1)，常数δ>0取0.1h～10h；
航天器俯仰轴的总干扰角加速度的理论值
Figure BDA0003038403670000046
ad2的估计值
Figure BDA0003038403670000047
即求得的z32为航天器俯仰轴的总干扰角加速度的估计值。
所述步骤(2)中，建立用于估计航天器偏航轴的总干扰角加速度的扩张状态观测器Σ3如下：
Figure BDA0003038403670000048
其中，z13，z23，z33为三个待求的未知函数；
Figure BDA0003038403670000051
为z13的时间导数，
Figure BDA0003038403670000052
为z23的时间导数，
Figure BDA0003038403670000053
为z33的时间导数；
β13为一常数，范围为h-1～10h-1；
β23为一常数，范围为0.01β13～0.1β13；
β33为一常数，范围为0.0001β13～0.001β13；
函数fal(e3,η,δ)＝e3s/δ1-η+|e3|ηsign(e3)(1-s)，s＝0.5(sign(e3+δ)-sign(e3-δ))，常数η∈(0,1)，常数δ>0取0.1h～10h；
航天器偏航轴的总干扰角加速度的理论值
Figure BDA0003038403670000054
ad3的估计值
Figure BDA0003038403670000055
即求得的z33为航天器偏航轴的总干扰角加速度的估计值。
6、根据权利要求1所述的一种抗慢变干扰的相平面姿态控制方法，其特征在于：
所述步骤(3)中，构建三轴的相平面，分别为：滚动轴相平面、俯仰轴相平面和偏航轴相平面，滚动轴相平面由滚动角x1和滚动角速度
Figure BDA00030384036700000517
构成，俯仰轴相平面由俯仰角x2和俯仰角速度
Figure BDA00030384036700000518
构成，偏航轴相平面由偏航角x3和偏航角速度
Figure BDA00030384036700000519
构成；
将滚动轴相平面分为七个区，七个区分别为：关机区S01，右侧大推力全开区S11，右侧小推力步进区S21，右侧小推力抗外滑区S31，左侧大推力全开区S11’，左侧小推力步进区S21’，左侧小推力抗外滑区S31’；
以x1＝0为纵坐标轴，
Figure BDA0003038403670000056
为横坐标轴；
所述S11区为直线
Figure BDA0003038403670000057
的右侧区域与直线
Figure BDA0003038403670000058
的上侧区域的重叠区域，以及直线
Figure BDA0003038403670000059
的上侧区域与直线x1＝xB的右侧区域的重叠区域，其中，直线
Figure BDA00030384036700000510
为右侧开关线，xD>0为死区阈值，
Figure BDA00030384036700000511
为所述小推力步进区的角速度最大值，K1>0为滚动轴相平面的开关线斜率绝对值，取值范围为
Figure BDA00030384036700000512
其中
Figure BDA00030384036700000513
为限速值，Tmin表示最小喷气长度；
所述S11’区为所述直线
Figure BDA00030384036700000514
的左侧区域与直线
Figure BDA00030384036700000515
下侧区域的重叠区域，以及直线
Figure BDA00030384036700000516
的下侧区域与直线x1＝-xB的左侧区域的重叠区域，其中直线
Figure BDA0003038403670000061
为左侧开关线；
所述S21区为直线
Figure BDA0003038403670000062
的右侧区域与直线
Figure BDA0003038403670000063
的上侧区域与直线
Figure BDA0003038403670000064
的下侧区域与直线x1＝xB的左侧区域的重叠区域，其中xB>0为所述大推力全开区的阈值；
所述S21’区为直线
Figure BDA0003038403670000065
的左侧区域与直线
Figure BDA0003038403670000066
下侧区域与直线
Figure BDA0003038403670000067
的上侧区域与直线x1＝-xB的右侧区域的重叠区域；
所述S31区为直线
Figure BDA00030384036700000625
的下侧区域与直线
Figure BDA0003038403670000068
的上侧区域与直线x1＝xD的右侧区域的重叠区域；
所述S31’区为直线
Figure BDA0003038403670000069
的上侧区域与直线
Figure BDA00030384036700000610
的下侧区域与直线x1＝-xD的左侧区域的重叠区域；
所述S01区为滚动轴相平面全区域除去所述的S11，S21，S31，S11’，S21’，S31’区以外的区域。
将俯仰轴相平面分为七个区，七个区分别为：关机区S02，右侧大推力全开区S12，右侧小推力步进区S22，右侧小推力抗外滑区S32，左侧大推力全开区S12’，左侧小推力步进区S22’，左侧小推力抗外滑区S32’；
以x2＝0为纵坐标轴，
Figure BDA00030384036700000611
为横坐标轴；
所述S12区为直线
Figure BDA00030384036700000612
的右侧区域与直线
Figure BDA00030384036700000613
的上侧区域的重叠区域，以及直线
Figure BDA00030384036700000614
的上侧区域与直线x2＝xB的右侧区域的重叠区域，其中，直线
Figure BDA00030384036700000615
为右侧开关线，xD>0为死区阈值，
Figure BDA00030384036700000616
为所述小推力步进区的角速度最大值，K2>0为俯仰轴相平面的开关线斜率绝对值，取值范围为
Figure BDA00030384036700000617
其中
Figure BDA00030384036700000618
为限速值，Tmin表示最小喷气长度；
所述S12’区为所述直线
Figure BDA00030384036700000619
的左侧区域与直线
Figure BDA00030384036700000620
下侧区域的重叠区域，以及直线
Figure BDA00030384036700000621
的下侧区域与直线x2＝-xB的左侧区域的重叠区域，其中直线
Figure BDA00030384036700000622
为左侧开关线；
所述S22区为直线
Figure BDA00030384036700000623
的右侧区域与直线
Figure BDA00030384036700000624
的上侧区域与直线
Figure BDA0003038403670000071
的下侧区域与直线x2＝xB的左侧区域的重叠区域，其中xB>0为所述大推力全开区的阈值；
所述S22’区为直线
Figure BDA0003038403670000072
的左侧区域与直线
Figure BDA0003038403670000073
下侧区域与直线
Figure BDA0003038403670000074
的上侧区域与直线x2＝-xB的右侧区域的重叠区域；
所述S32区为直线
Figure BDA0003038403670000075
的下侧区域与直线
Figure BDA0003038403670000076
的上侧区域与直线x2＝xD的右侧区域的重叠区域；
所述S32’区为直线
Figure BDA0003038403670000077
的上侧区域与直线
Figure BDA0003038403670000078
的下侧区域与直线x2＝-xD的左侧区域的重叠区域；
所述S02区为俯仰轴相平面全区域除去所述的S12，S22，S32，S12’，S22’，S32’区以外的区域。
将偏航轴相平面分为七个区，七个区分别为：关机区S03，右侧大推力全开区S13，右侧小推力步进区S23，右侧小推力抗外滑区S33，左侧大推力全开区S13’，左侧小推力步进区S23’，左侧小推力抗外滑区S33’；
以x3＝0为纵坐标轴，
Figure BDA0003038403670000079
为横坐标轴；
所述S13区为直线
Figure BDA00030384036700000710
的右侧区域与直线
Figure BDA00030384036700000711
的上侧区域的重叠区域，以及直线
Figure BDA00030384036700000712
的上侧区域与直线x3＝xB的右侧区域的重叠区域，其中，直线
Figure BDA00030384036700000713
为右侧开关线，xD>0为死区阈值，
Figure BDA00030384036700000714
为所述小推力步进区的角速度最大值，K3>0为偏航轴相平面的开关线斜率绝对值，取值范围为
Figure BDA00030384036700000715
其中
Figure BDA00030384036700000716
为限速值，Tmin表示最小喷气长度；
所述S13’区为所述直线
Figure BDA00030384036700000717
的左侧区域与直线
Figure BDA00030384036700000718
下侧区域的重叠区域，以及直线
Figure BDA00030384036700000719
的下侧区域与直线x3＝-xB的左侧区域的重叠区域，其中直线
Figure BDA00030384036700000720
为左侧开关线；
所述S23区为直线
Figure BDA00030384036700000721
的右侧区域与直线
Figure BDA00030384036700000722
的上侧区域与直线
Figure BDA00030384036700000723
的下侧区域与直线x3＝xB的左侧区域的重叠区域，其中xB>0为所述大推力全开区的阈值；
所述S23’区为直线
Figure BDA0003038403670000081
的左侧区域与直线
Figure BDA0003038403670000082
下侧区域与直线
Figure BDA0003038403670000083
的上侧区域与直线x3＝-xB的右侧区域的重叠区域；
所述S33区为直线
Figure BDA0003038403670000084
的下侧区域与直线
Figure BDA0003038403670000085
的上侧区域与直线x3＝xD的右侧区域的重叠区域；
所述S33’区为直线
Figure BDA0003038403670000086
的上侧区域与直线
Figure BDA0003038403670000087
的下侧区域与直线x3＝-xD的左侧区域的重叠区域；
所述S03区为偏航轴相平面全区域除去所述的S13，S23，S33，S13’，S23’，S33’区以外的区域。
所述步骤(4)中，关机区的相平面控制指令为：
针对滚动轴，在当前采样时刻tk，tk＝t0+kΔT，k＝0,1,2,...表示采样时刻，其中t0表示第一次采样时刻，当
Figure BDA0003038403670000088
位于所述S01区时，则关闭产生绕滚动轴力矩的推力器，不进行喷气；
针对俯仰轴，在当前采样时刻tk，当
Figure BDA0003038403670000089
位于所述S02区时，则关闭产生绕俯仰轴力矩的推力器，不进行喷气；
针对偏航轴，在当前采样时刻tk，当
Figure BDA00030384036700000810
位于所述S03区时，则关闭产生绕偏航轴力矩的大推力器和小推力器，不进行喷气；
所述步骤(4)中，右侧大推力全开区的相平面控制指令为：
针对滚动轴，在当前采样时刻tk，当
Figure BDA00030384036700000811
位于所述S11区时，则采用大推力器喷气，产生绕滚动轴力矩为-FM1，喷气长度为ΔT，其中FM1为所配置的大推力器喷气产生的绕滚动轴的力矩，ΔT表示采样控制周期；
针对俯仰轴，在当前采样时刻tk，当
Figure BDA00030384036700000812
位于所述S12区时，则采用大推力器喷气，产生绕俯仰轴力矩为-FM2，喷气长度为ΔT，其中FM2为所配置的大推力器喷气产生的绕俯仰轴的力矩；
针对偏航轴，在当前采样时刻tk，当
Figure BDA00030384036700000813
位于所述S13区时，则采用大推力器喷气，产生绕偏航轴力矩为-FM3，喷气长度为ΔT，其中FM3为所配置的大推力器喷气产生的绕偏航轴的力矩；
所述步骤(4)中，左侧大推力全开区的相平面控制指令为：
针对滚动轴，在当前采样时刻tk，当
Figure BDA0003038403670000091
位于所述S11’区时，则采用大推力器喷气，产生绕滚动轴力矩为FM1，喷气长度为ΔT；
针对俯仰轴，在当前采样时刻tk，当
Figure BDA0003038403670000092
位于所述S12’区时，则采用大推力器喷气，产生绕俯仰轴力矩为FM2，喷气长度为ΔT；
针对偏航轴，在当前采样时刻tk，当
Figure BDA0003038403670000093
位于所述S13’区时，则采用大推力器喷气，产生绕偏航轴力矩为FM3，喷气长度为ΔT；
所述步骤(4)中，右侧小推力步进区的相平面控制指令为：
(i)针对滚动轴，在当前采样时刻tk，当
Figure BDA0003038403670000094
位于所述S21：
若满足
Figure BDA0003038403670000095
且
Figure BDA0003038403670000096
时，则本周期关闭产生绕滚动轴力矩的小推力器和大推力器，不进行喷气，其中ε1为滚动轴干扰阈值，取值范围为(0,0.01Fm1/I11]；
若满足
Figure BDA0003038403670000097
且
Figure BDA0003038403670000098
时，则本周期采用小推力器喷气，产生绕滚动轴力矩为-Fm1，喷气长度TN为：
Figure BDA0003038403670000099
TN2∈[Tmin,0.8ΔT]，TN∈[Tmin,ΔT]，其中K11>0为滚动轴相平面小推力步进区干扰补偿系数，
Figure BDA00030384036700000910
为小推力角加速度参数，k2∈[0.5,0.9]为常数；
若满足
Figure BDA00030384036700000911
且
Figure BDA00030384036700000912
时，则本周期采用小推力器喷气，产生绕滚动轴力矩为-Fm1，喷气长度为Tmin；
若满足
Figure BDA00030384036700000913
且
Figure BDA00030384036700000914
时，则本周期采用小推力器喷气，产生绕滚动轴力矩为-Fm1，喷气长度TN为：TN＝TN1+TN2，
Figure BDA00030384036700000915
TN2＝Kj(|x1|-xe)/ac2，TN2∈[Tmin,0.5ΔT]，TN∈[Tmin,ΔT]，其中Kj＝((1-k2)ac2ΔT)/(xB-xe)为小推力步进区参数，xe>0为略小于xD的常数；
(ii)针对俯仰轴，在当前采样时刻tk，当
Figure BDA00030384036700000916
位于所述S22：
若满足
Figure BDA0003038403670000101
且
Figure BDA0003038403670000102
时，则本周期关闭产生绕俯仰轴力矩的小推力器和大推力器，不进行喷气，其中ε2为俯仰轴干扰阈值，取值范围为(0,0.01Fm2/I22]；
若满足
Figure BDA0003038403670000103
且
Figure BDA0003038403670000104
时，则本周期采用小推力器喷气，产生绕俯仰轴力矩为-Fm2，喷气长度TN为：TN＝TN1+TN2，
Figure BDA0003038403670000105
Figure BDA0003038403670000106
TN2∈[Tmin,0.8ΔT]，TN∈[Tmin,ΔT]，其中K12>0为俯仰轴相平面小推力步进区干扰补偿系数；
若满足
Figure BDA0003038403670000107
且
Figure BDA0003038403670000108
时，则本周期采用小推力器喷气，产生绕俯仰轴力矩为-Fm2，喷气长度为Tmin；
若满足
Figure BDA0003038403670000109
且
Figure BDA00030384036700001010
时，则本周期采用小推力器喷气，产生绕俯仰轴力矩为-Fm2，喷气长度TN为：TN＝TN1+TN2，
Figure BDA00030384036700001011
TN2＝Kj(|x2|-xe)/ac2，TN2∈[Tmin,0.5ΔT]，TN∈[Tmin,ΔT]；
(iii)针对偏航轴，在当前采样时刻tk，当
Figure BDA00030384036700001012
位于所述S23：
若满足
Figure BDA00030384036700001013
且
Figure BDA00030384036700001014
时，则本周期关闭产生绕偏航轴力矩的小推力器和大推力器，不进行喷气，其中ε3为偏航轴干扰阈值，取值范围为(0,0.01Fm3/I33]；
若满足
Figure BDA00030384036700001015
且
Figure BDA00030384036700001016
时，则本周期采用小推力器喷气，产生绕偏航轴力矩为-Fm3，喷气长度TN为：TN＝TN1+TN2，
Figure BDA00030384036700001017
TN2∈[Tmin,0.8ΔT]，TN∈[Tmin,ΔT]，其中K13>0为俯仰轴相平面小推力步进区干扰补偿系数；
若满足
Figure BDA00030384036700001018
且
Figure BDA00030384036700001019
时，则本周期采用小推力器喷气，产生绕偏航轴力矩为-Fm3，喷气长度为Tmin；
若满足
Figure BDA00030384036700001020
且
Figure BDA00030384036700001021
时，则本周期采用小推力器喷气，产生绕偏航轴力矩为-Fm3，喷气长度TN为：TN＝TN1+TN2，
Figure BDA00030384036700001022
TN2＝Kj(|x3|-xe)/ac2，TN2∈[Tmin,0.5ΔT]，TN∈[Tmin,ΔT]；
所述步骤(4)中，左侧小推力步进区的相平面控制指令为：
(i’)针对滚动轴，在当前采样时刻tk，当
Figure BDA0003038403670000111
位于所述S21’:
若满足
Figure BDA0003038403670000112
且
Figure BDA0003038403670000113
时，则本周期关闭产生绕滚动轴力矩的小推力器和大推力器，不进行喷气；
若满足
Figure BDA0003038403670000114
且
Figure BDA0003038403670000115
时，则本周期采用小推力器喷气，产生绕滚动轴力矩为Fm1，喷气长度TN＝TN1+TN2，
Figure BDA0003038403670000116
若满足
Figure BDA0003038403670000117
且
Figure BDA0003038403670000118
时，则本周期采用小推力器喷气，产生绕滚动轴力矩为Fm1，喷气长度为Tmin；
若满足
Figure BDA0003038403670000119
且
Figure BDA00030384036700001110
时，则本周期采用小推力器喷气，产生绕滚动轴力矩为Fm1，喷气长度TN为：TN＝TN1+TN2，
Figure BDA00030384036700001111
TN2＝Kj(|x1|-xe)/ac2，TN2∈[Tmin,0.5ΔT]，TN∈[Tmin,ΔT]；
(ii’)针对俯仰轴，在当前采样时刻tk，当
Figure BDA00030384036700001112
位于所述S22’：
若满足
Figure BDA00030384036700001113
且
Figure BDA00030384036700001114
时，则本周期关闭产生绕俯仰轴力矩的小推力器和大推力器，不进行喷气；
若满足
Figure BDA00030384036700001115
且
Figure BDA00030384036700001116
时，则本周期采用小推力器喷气，产生绕俯仰轴力矩为Fm2，喷气长度TN＝TN1+TN2，
Figure BDA00030384036700001117
若满足
Figure BDA00030384036700001118
且
Figure BDA00030384036700001119
时，则本周期采用小推力器喷气，产生绕俯仰轴力矩为Fm2，喷气长度为Tmin；
若满足
Figure BDA00030384036700001120
且
Figure BDA00030384036700001121
时，则本周期采用小推力器喷气，产生绕俯仰轴力矩为Fm2，喷气长度TN为：TN＝TN1+TN2，
Figure BDA00030384036700001122
TN2＝Kj(|x2|-xe)/ac2，TN2∈[Tmin,0.5ΔT]，TN∈[Tmin,ΔT]；
(iii’)针对偏航轴，在当前采样时刻tk，当
Figure BDA00030384036700001123
位于所述S23’：
若满足
Figure BDA00030384036700001124
且
Figure BDA00030384036700001125
时，则本周期关闭产生绕偏航轴力矩的小推力器和大推力器，不进行喷气；
若满足
Figure BDA0003038403670000121
且
Figure BDA0003038403670000122
时，则本周期采用小推力器喷气，产生绕偏航轴力矩为Fm3，喷气长度TN＝TN1+TN2，
Figure BDA0003038403670000123
若满足
Figure BDA0003038403670000124
且
Figure BDA0003038403670000125
时，则本周期采用小推力器喷气，产生绕偏航轴力矩为Fm3，喷气长度为Tmin；
若满足
Figure BDA0003038403670000126
且
Figure BDA0003038403670000127
时，则本周期采用小推力器喷气，产生绕偏航轴力矩为Fm3，喷气长度TN为：TN＝TN1+TN2，
Figure BDA0003038403670000128
TN2＝Kj(|x3|-xe)/ac2，TN2∈[Tmin,0.5ΔT]，TN∈[Tmin,ΔT]；
所述步骤(4)中，右侧小推力抗外滑区的相平面控制指令为：
(I)针对滚动轴，在当前采样时刻tk，当
Figure BDA0003038403670000129
位于所述S31：
若满足
Figure BDA00030384036700001210
时，则本周期关闭产生绕滚动轴力矩的小推力器和大推力器，不进行喷气；
若满足
Figure BDA00030384036700001211
时，则本周期采用小推力器喷气，产生绕滚动轴力矩为-Fm1，喷气长度TN为：TN＝TN1+TN2，TN1＝Tmin，
Figure BDA00030384036700001212
TN2∈[0,0.8ΔT]，TN∈[Tmin,ΔT]，其中K21为滚动轴小推力抗外滑区干扰补偿系数，满足0<K21<K11；
(II)针对俯仰轴，在当前采样时刻tk，当
Figure BDA00030384036700001213
位于所述S32：
若满足
Figure BDA00030384036700001214
时，则本周期关闭产生绕俯仰轴力矩的小推力器和大推力器，不进行喷气；
若满足
Figure BDA00030384036700001215
时，则本周期采用小推力器喷气，产生绕俯仰轴力矩为-Fm2，喷气长度TN为：TN＝TN1+TN2，TN1＝Tmin，
Figure BDA00030384036700001216
TN2∈[0,0.8ΔT]，TN∈[Tmin,ΔT]，其中K22为俯仰轴小推力抗外滑区干扰补偿系数，满足0<K22<K12；
(III)针对偏航轴，在当前采样时刻tk，当
Figure BDA00030384036700001217
位于所述S33：
若满足
Figure BDA0003038403670000131
时，则本周期关闭产生绕偏航轴力矩的小推力器和大推力器，不进行喷气；
若满足
Figure BDA0003038403670000132
时，则本周期采用小推力器喷气，产生绕偏航轴力矩为-Fm3，喷气长度TN为：TN＝TN1+TN2，TN1＝Tmin，
Figure BDA0003038403670000133
TN2∈[0,0.8ΔT]，TN∈[Tmin,ΔT]，其中K23为偏航轴小推力抗外滑区干扰补偿系数，满足0<K23<K13；
所述步骤(4)中，左侧小推力抗外滑区的相平面控制指令为：
(I’)针对滚动轴，在当前采样时刻tk，当
Figure BDA0003038403670000134
位于所述S31’：
若满足
Figure BDA0003038403670000135
时，则本周期关闭产生绕滚动轴力矩的小推力器和大推力器，不进行喷气；
若满足
Figure BDA0003038403670000136
时，则本周期采用小推力器喷气，产生绕滚动轴力矩为-Fm1，喷气长度TN：TN＝TN1+TN2，TN1＝Tmin，
Figure BDA0003038403670000137
TN2∈[0,0.8ΔT]，TN∈[Tmin,ΔT]；
(II’)针对俯仰轴，在当前采样时刻tk，当
Figure BDA0003038403670000138
位于所述S32’：
若满足
Figure BDA0003038403670000139
时，则本周期关闭产生绕俯仰轴力矩的小推力器和大推力器，不进行喷气；
若满足
Figure BDA00030384036700001310
时，则本周期采用小推力器喷气，产生绕俯仰轴力矩为-Fm2，喷气长度TN：TN＝TN1+TN2，TN1＝Tmin，
Figure BDA00030384036700001311
TN2∈[0,0.8ΔT]，TN∈[Tmin,ΔT]；
(III’)针对偏航轴，在当前采样时刻tk，当
Figure BDA00030384036700001312
位于所述S33’：
若满足
Figure BDA00030384036700001313
时，则本周期关闭产生绕偏航轴力矩的小推力器和大推力器，不进行喷气；
若满足
Figure BDA00030384036700001314
时，则本周期采用小推力器喷气，产生绕偏航轴力矩为-Fm3，喷气长度TN：TN＝TN1+TN2，TN1＝Tmin，
Figure BDA00030384036700001315
TN2∈[0,0.8ΔT]，TN∈[Tmin,ΔT]；
所述步骤(4)中，在计算右侧小推力步进区和左侧小推力步进区的相平面控制指令中的三轴相平面控制喷气长度之前，对三轴干扰角加速度估计进行限幅，具体限幅方法为：
针对滚动轴，当
Figure BDA0003038403670000141
时，则取
Figure BDA0003038403670000142
当
Figure BDA0003038403670000143
时，则取
Figure BDA0003038403670000144
其中限幅值Lmt1取0.1Fm1/I11～0.5Fm1/I11；
针对俯仰轴，当
Figure BDA0003038403670000145
时，则取
Figure BDA0003038403670000146
当
Figure BDA0003038403670000147
时，则取
Figure BDA0003038403670000148
其中限幅值Lmt2取0.1Fm2/I22～0.5Fm2/I22；
针对偏航轴，当
Figure BDA0003038403670000149
时，则取
Figure BDA00030384036700001410
当
Figure BDA00030384036700001411
时，则取
Figure BDA00030384036700001412
其中限幅值Lmt3取0.1Fm3/I33～0.5Fm3/I33；
在根据航天器三轴姿态角测量输出值设计所述扩张状态观测器之前，还对航天器三轴姿态角测量输出值进行了惯性滤波处理，具体滤波方法为：
针对滚动轴，设计惯性滤波形式为：
Figure BDA00030384036700001413
其中，参数α1>0表示滚动轴姿态滤波系数，yf1为滚动姿态角滤波输出值；
针对俯仰轴，设计惯性滤波形式为：
Figure BDA00030384036700001414
其中，参数α2>0表示俯仰轴姿态滤波系数，yf2为俯仰姿态角滤波输出值；
针对偏航轴，设计惯性滤波形式为：
Figure BDA00030384036700001415
其中，参数α3>0表示偏航轴姿态滤波系数，yf3为偏航姿态角滤波输出值。
所述步骤(5)中，航天器的推力器根据步骤(4)得到的三轴相平面控制指令进行喷气，具体为：
针对航天器的滚动轴，在当前采样时刻tk和下一个采样时刻tk+ΔT之间所构成的采样控制周期[tk,tk+ΔT]内，在当前采样时刻将滚动轴相平面控制指令中所采用的推力器开机，开机时间长度等于滚动轴相平面控制指令中的滚动轴喷气长度，然后将所采用的推力器关机；
针对航天器的俯仰轴，在当前采样时刻tk和下一个采样时刻tk+ΔT之间所构成的采样控制周期[tk,tk+ΔT]内，在当前采样时刻将俯仰轴相平面控制指令中所采用的推力器开机，开机时间长度等于俯仰轴相平面控制指令中的俯仰轴喷气长度，然后将所采用的推力器关机；
针对航天器的偏航轴，在当前采样时刻tk和下一个采样时刻tk+ΔT之间所构成的采样控制周期[tk,tk+ΔT]内，在当前采样时刻将偏航轴相平面控制指令中所采用的推力器开机，开机时间长度等于偏航轴相平面控制指令中的偏航轴喷气长度，然后将所采用的推力器关机。
本发明与现有技术相比的有益效果是：
本发明提出的一种抗慢变干扰的相平面控制方法，与目前常用的相平面控制设计方法相比，引入了根据实时干扰估计对相平面控制补偿量进行适应性调整的策略。通过快速准确估计干扰，进而根据干扰估计调整相平面控制补偿量，实现对干扰补偿控制的精细设计。因此，得到的新方法对干扰具有更强的适应性，同时兼顾闭环系统稳态性能，减少了喷气次数和燃料消耗。
附图说明
图1为本发明的抗慢变干扰的相平面控制方法流程图；
图2为本发明的第i轴(i＝1,2,3分别表示滚动轴、俯仰轴、偏航轴)相平面分区示意图；
图3为本发明控制方案对施加在三轴的外部慢变干扰力矩估计曲线；
图4为基于干扰角加速度估计器的相平面控制方法对施加在三轴的外部慢变干扰力矩估计曲线；
图5为本发明控制方案下的三轴姿态角变化曲线；
图6为基于干扰角加速度估计器的相平面控制方法下的三轴姿态角变化曲线；
图7为本发明控制方案下的三轴喷气控制力矩曲线；
图8为基于干扰角加速度估计器的相平面控制方法下的三轴喷气控制力矩曲线。
具体实施方式
下面结合附图对本发明的具体实施方式进行进一步的详细描述。
如图1所示，本发明的具体步骤实现如下(以下以航天器姿态稳定任务为例来说明方法的具体实现)：
(1)建立慢变干扰下的航天器动力学模型。
航天器姿态运动学方程为：
Figure BDA0003038403670000161
其中，
Figure BDA0003038403670000162
θ,ψ分别为航天器滚动角、俯仰角以及偏航角，ω1,ω2,ω3为航天器相对惯性空间的瞬时转速在本体系的三轴分量，ω0为轨道角速度。
航天器姿态动力学方程为
Figure BDA0003038403670000163
其中，
Figure BDA0003038403670000164
为航天器转动惯量阵，ω＝[ω1,ω2,ω3]T，u＝[u1,u2,u3]T和Td＝[Td1,Td2,Td3]T分别为控制力矩和外部慢变干扰力矩，
Figure BDA0003038403670000165
根据上述航天器姿态运动学方程可得：
Figure BDA0003038403670000166
其中，Cbo表示挠性航天器轨道坐标系到本体坐标系的坐标变换矩阵，具体表达式为：
Figure BDA0003038403670000167
函数矩阵
Figure BDA0003038403670000171
的具体表达式为：
Figure BDA0003038403670000172
根据上述航天器姿态运动学方程和动力学方程可得：
Figure BDA0003038403670000173
其中
Figure BDA0003038403670000174
表示函数矩阵
Figure BDA0003038403670000175
的逆矩阵，
Figure BDA0003038403670000176
表示函数矩阵
Figure BDA0003038403670000177
的时间导数，
Figure BDA0003038403670000178
表示Cbo的时间导数。
因此，可以得到慢变干扰下的航天器动力学模型为：
Figure BDA0003038403670000179
其中,y＝[y1,y2,y3]表示三轴姿态输出测量，采用数字太阳敏感器进行测量，v＝[v1,v2,v3]T为测量噪声，
Figure BDA00030384036700001710
和
Figure BDA00030384036700001711
分别为三轴相平面控制力矩指令，非线性函数
Figure BDA00030384036700001712
和
Figure BDA00030384036700001713
分别表示由于航天器所受外部慢变干扰Td、航天器三轴动力学耦合及三轴控制输入耦合带来的总干扰角加速度在三轴的分量，具体表达式为：
Figure BDA00030384036700001714
(2)设计扩张状态观测器估计三轴的总干扰角加速度。
下面说明针对第i，i＝1，2，3轴的扩张状态观测器的设计过程。基于第(1)步得到的慢变干扰下的航天器第i轴动力学模型可得：
Figure BDA0003038403670000181
其中，i＝1,2,3分别对应滚动轴、俯仰轴及偏航轴，adi表示航天器第i轴总干扰角加速度的理论值，即
Figure BDA0003038403670000182
x2＝θ，
Figure BDA0003038403670000183
x3＝ψ，
Figure BDA0003038403670000184
vi表示第i轴姿态角的测量噪声。
为减少测量噪声影响，采用惯性滤波方法进行处理：
Figure BDA0003038403670000185
其中，参数αi>0表示第i轴姿态滤波系数，yfi为第i轴姿态角测量yi经过惯性滤波后的输出结果。
设计扩张状态观测器的形式为
Figure BDA0003038403670000186
其中，函数
Figure BDA0003038403670000187
参数δ>0，η>0，β1i,β2i,β3i为扩张状态观测器增益，adi的估计值
Figure BDA0003038403670000188
即求得的z3i为航天器第i轴总干扰角加速度的估计值
Figure BDA0003038403670000189
(3)针对由三轴的姿态角及姿态角速度所构成的相平面，设计三轴的相平面分区。
以第i轴为例进行说明，第i轴姿态的相平面控制的分区呈中心对称，如图2所示，包括关机区S0i，右侧大推力全开区S1i，右侧小推力步进区S2i，右侧小推力抗外滑区S3i，右侧开关线为GCD，左侧大推力全开区S1i'，左侧小推力步进区S2i'，左侧小推力抗外滑区S3i'，左侧开关线为G'C'D'。
(4)针对相平面各分区设计三轴相平面控制指令。
在当前采样时刻tk，根据第i轴的总干扰角加速度估计值
Figure BDA0003038403670000191
以及姿态角xi和姿态角速度
Figure BDA0003038403670000192
构成的相点
Figure BDA0003038403670000193
作为第i轴的相平面控制输入，第i轴的相平面控制指令
Figure BDA0003038403670000194
具体如下表1所示。
Figure BDA0003038403670000195
Figure BDA0003038403670000201
其中，εi为第i,i＝1,2,3轴干扰阈值，取值范围为(0,0.01Fmi/Iii]，FMi和Fmi分别为所配置的大、小推力器产生的绕第i轴正向的力矩，Iii为转动惯量矩阵I的第i个对角元素。
具体地，在进行第i轴相平面控制指令中的喷气长度TN计算之前对
Figure BDA0003038403670000202
进行限幅，具体限幅方法为：当
Figure BDA0003038403670000203
时，取
Figure BDA0003038403670000204
当
Figure BDA0003038403670000205
时，取
Figure BDA0003038403670000206
其中限幅值Lmti选取为0.1Fmi/Iii～0.5Fmi/Iii。
(5)航天器的推力器根据步骤(4)得到的三轴相平面控制指令进行喷气。具体为：
针对航天器的滚动轴，在当前采样时刻tk和下一个采样时刻tk+ΔT之间所构成的采样控制周期[tk,tk+ΔT]内，在当前采样时刻将所采用的推力器开机，开机时间长度等于滚动轴喷气长度，然后将所采用的推力器关机；
针对航天器的俯仰轴，在当前采样时刻tk和下一个采样时刻tk+ΔT之间所构成的采样控制周期[tk,tk+ΔT]内，在当前采样时刻将所采用的推力器开机，开机时间长度等于俯仰轴喷气长度，然后将所采用的推力器关机；
针对航天器的偏航轴，在当前采样时刻tk和下一个采样时刻tk+ΔT之间所构成的采样控制周期[tk,tk+ΔT]内，在当前采样时刻将所采用的推力器开机，开机时间长度等于偏航轴喷气长度，然后将所采用的推力器关机。
本发明提出的抗慢变干扰相平面控制方法或系统工程易实现，同时扩张状态观测器以及相平面控制方法，都已在工程实践中得到了广泛验证，具有实用性。本方法可推广应用于采用喷气推力器完成姿态控制任务的各类航天器。
实施例：
通过方法一和方法二的对比仿真说明本发明优点。
方法一：本发明方法。具体设计如下：
(1)建立慢变干扰下的航天器动力学模型为
Figure BDA0003038403670000211
其中，航天器转动惯量矩阵为
Figure BDA0003038403670000212
外部慢变干扰取为
Figure BDA0003038403670000213
轨道角速度ω0＝0.0011rad/s，数字太阳敏感器的测量噪声v1,v2,v3表示均值为零方差为1×10-7rad的高斯白噪声。
(2)设计扩张状态观测器估计三轴的总干扰角加速度
先对三轴姿态测量进行惯性滤波
Figure BDA0003038403670000221
其中，yfi表示yi经过惯性滤波后的输出，三轴姿态惯性滤波系数取为α1＝4，α2＝4，α3＝4。
再设计三轴的扩张状态观测器为
Figure BDA0003038403670000222
其中，函数
Figure BDA0003038403670000223
扩张状态观测器数值离散步长h＝0.001，δ＝0.001，η＝0.5，三轴扩张状态观测器增益分别为β11＝1000,β21＝150,β31＝0.25，β12＝1000,β22＝20,β32＝1，β13＝1000,β23＝20,β33＝1，求得的z3i为航天器第i轴总干扰角加速度的估计值
Figure BDA0003038403670000224
即航天器第i轴总干扰角加速度的理论值
Figure BDA0003038403670000225
的估计值
Figure BDA0003038403670000226
(3)设计三轴相平面分区，具体相平面分区参数如下表2所示
Figure BDA0003038403670000227
(4)根据表1及三轴相平面分区，设计三轴相平面控制指令，涉及的三轴相平面控制参数如下表3所示。
Figure BDA0003038403670000228
Figure BDA0003038403670000231
(5)航天器的推力器根据上述步骤(4)得到的三轴相平面控制指令进行喷气。具体为：
针对航天器的滚动轴，在当前采样时刻tk和下一个采样时刻tk+ΔT之间所构成的采样控制周期[tk,tk+ΔT]内，在当前采样时刻将所采用的推力器开机，开机时间长度等于滚动轴喷气长度，然后将所采用的推力器关机；
针对航天器的俯仰轴，在当前采样时刻tk和下一个采样时刻tk+ΔT之间所构成的采样控制周期[tk,tk+ΔT]内，在当前采样时刻将所采用的推力器开机，开机时间长度等于俯仰轴喷气长度，然后将所采用的推力器关机；
针对航天器的偏航轴，在当前采样时刻tk和下一个采样时刻tk+ΔT之间所构成的采样控制周期[tk,tk+ΔT]内，在当前采样时刻将所采用的推力器开机，开机时间长度等于偏航轴喷气长度，然后将所采用的推力器关机。
方法二：《卫星姿态动力学与控制》(屠善澄主编，宇航出版社，2001)第443页的基于干扰角加速度估计器的相平面控制方法。具体设计方法如下：
(I)以俯仰轴为例，滚动轴和偏航轴的设计方法类似。采用干扰角加速度估计器估计俯仰轴的干扰角加速度，其中状态方程为
Figure BDA0003038403670000241
其中wd是方差不为零的虚拟白噪声，而量测方程为
Figure BDA0003038403670000242
式中ω0为轨道角速度，gθ为由陀螺得到的速率信号，
Figure BDA0003038403670000243
为由姿态测量系统提供的陀螺漂移估值，设置陀螺测量噪声为均值为零方差为1×10-7rad/s的高斯白噪声。
干扰角加速度估计器为
Figure BDA0003038403670000244
其中，干扰角速度估计器增益阵K＝[1000,10]T。
(II)采用和方法一同样的相平面分区和相平面控制参数。由于方法二没有在小推力抗外滑区对干扰进行补偿，因此在方法二中将三轴相平面控制中的小推力抗外滑区的喷气长度计算公式中的
Figure BDA0003038403670000245
i＝1,2,3置为0。
取系统仿真初值
Figure BDA0003038403670000246
单位为rad和rad/s。
根据仿真结果，得到如下结论：
1、本发明提出的抗慢变干扰相平面控制方案是可行的。图3给出了施加在三轴上的外部干扰力矩的估计情况，表明本方案可以快速准确地估计出外部慢变干扰。如图5给出了三轴姿态控制情况，表明本发明方案可以对干扰进行有效补偿，三轴姿态控制的稳态性能较好。图7给出了三轴喷气情况，表明在进行干扰补偿过程中俯仰轴和偏航轴喷气次数较少，没有造成明显的喷气频繁和对喷现象。
2、方法二难以对外部慢变干扰进行准确估计，导致闭环系统稳态性能难以达到满意效果。如图4所示，当施加在俯仰轴和偏航轴的外部慢变干扰较大时，干扰估计结果与理论值存在较大偏差。图6给出了在控制补偿作用下的三轴姿态变化情况，表明俯仰角和偏航角的稳态控制误差较大。图8给出了三轴喷气情况，表明在进行干扰补偿过程中，俯仰轴和偏航轴喷气频繁，消耗燃料较多。
本发明说明书中未作详细描述的内容属本领域技术人员的公知技术。

Claims (10)

1.一种抗慢变干扰的相平面姿态控制方法，其特征在于：该方法所针对的航天器上配置有大推力器和小推力器，该具体步骤包括：
(1)建立慢变干扰下的航天器动力学模型；
(2)建立用于估计航天器三轴的总干扰角加速度的扩张状态观测器，并使用所建立的扩张状态观测器估计航天器三轴的总干扰角加速度；
(3)构建三轴的相平面，并将所构建的三轴相平面分别分为七个区；
(4)根据步骤(1)中所建立的慢变干扰下的航天器动力学模型、根据所述步骤(2)的扩张状态观测器估计出的三轴总干扰角加速度和步骤(3)的三轴相平面分区，得到每轴相平面的七个区的控制指令；
(5)航天器的推力器根据步骤(4)得到的三轴相平面控制指令进行喷气。
2.根据权利要求1所述的一种抗慢变干扰的相平面姿态控制方法，其特征在于：
所述步骤(1)的建立的慢变干扰下的航天器动力学模型为：
Figure FDA0003038403660000011
其中，x1,x2,x3表示慢变干扰下的航天器的三轴姿态角，三轴姿态角为滚动姿态角、俯仰姿态角和偏航姿态角；
Figure FDA0003038403660000012
表示慢变干扰下的航天器的三轴姿态角速度，三轴姿态角速度为滚动姿态角速度、俯仰姿态角速度和偏航姿态角速度；
Figure FDA0003038403660000013
表示慢变干扰下的航天器的三轴姿态角加速度，三轴姿态角加速度为滚动姿态角加速度、俯仰姿态角加速度和偏航姿态角加速度；
X＝[x1,x2,x3]T，
Figure FDA0003038403660000021
Figure FDA0003038403660000022
表示航天器的转动惯量矩阵；
Figure FDA0003038403660000023
表示慢变干扰下的航天器的三轴姿态控制指令；
Figure FDA0003038403660000024
d1,d2,d3表示慢变干扰在航天器的三轴分量；
D＝[d1,d2,d3]T
Figure FDA0003038403660000025
表示由于航天器所受外部慢变干扰D、航天器三轴动力学耦合及三轴控制输入耦合带来的总干扰角加速度在三轴的分量；
Y＝[y1,y2,y3]T，其中，y1,y2,y3表示慢变干扰下的航天器的三轴姿态角的测量值；
V＝[v1,v2,v3]T，其中，v1,v2,v3表示慢变干扰下的航天器的三轴姿态角的测量噪声。
3.根据权利要求1所述的一种抗慢变干扰的相平面姿态控制方法，其特征在于：
所述步骤(2)中，建立用于估计航天器滚动轴的总干扰角加速度的扩张状态观测器Σ1如下：
Figure FDA0003038403660000031
其中，z11，z21，z31为三个待求的未知函数；
Figure FDA0003038403660000032
为z11的时间导数，
Figure FDA0003038403660000033
为z2的时间导数，
Figure FDA0003038403660000034
为z3的时间导数；
β11为一常数，范围为h-1～10h-1；
β21为一常数，范围为0.01β11～0.1β11；
β31为一常数，范围为0.0001β11～0.001β11；
函数fal(e1,η,δ)＝e1s/δ1-η+|e1|ηsign(e1)(1-s)，s＝0.5(sign(e1+δ)-sign(e1-δ))，常数η∈(0,1)，常数δ>0取0.1h～10h，h表示数值离散步长，选取范围为0～0.01；
航天器滚动轴的总干扰角加速度的理论值
Figure FDA0003038403660000035
ad1的估计值
Figure FDA0003038403660000036
即求得的z31为航天器滚动轴的总干扰角加速度的估计值。
4.根据权利要求1所述的一种抗慢变干扰的相平面姿态控制方法，其特征在于：
所述步骤(2)中，建立用于估计航天器俯仰轴的总干扰角加速度的扩张状态观测器Σ2如下：
Figure FDA0003038403660000037
其中，z12，z22，z32为三个待求的未知函数；
Figure FDA0003038403660000038
为z12的时间导数，
Figure FDA0003038403660000039
为z22的时间导数，
Figure FDA00030384036600000310
为z32的时间导数；
β12为一常数，范围为h-1～10h-1；
β22为一常数，范围为0.01β12～0.1β12；
β32为一常数，范围为0.0001β12～0.001β12；
函数fal(e2,η,δ)＝e2s/δ1-η+|e2|ηsign(e2)(1-s)，s＝0.5(sign(e2+δ)-sign(e2-δ))，常数η∈(0,1)，常数δ>0取0.1h～10h；
航天器俯仰轴的总干扰角加速度的理论值
Figure FDA0003038403660000041
ad2的估计值
Figure FDA0003038403660000042
即求得的z32为航天器俯仰轴的总干扰角加速度的估计值。
5.根据权利要求1所述的一种抗慢变干扰的相平面姿态控制方法，其特征在于：
所述步骤(2)中，建立用于估计航天器偏航轴的总干扰角加速度的扩张状态观测器Σ3如下：
Figure FDA0003038403660000043
其中，z13，z23，z33为三个待求的未知函数；
Figure FDA0003038403660000044
为z13的时间导数，
Figure FDA0003038403660000045
为z23的时间导数，
Figure FDA0003038403660000046
为z33的时间导数；
β13为一常数，范围为h-1～10h-1；
β23为一常数，范围为0.01β13～0.1β13；
β33为一常数，范围为0.0001β13～0.001β13；
函数fal(e3,η,δ)＝e3s/δ1-η+|e3|ηsign(e3)(1-s)，s＝0.5(sign(e3+δ)-sign(e3-δ))，常数η∈(0,1)，常数δ>0取0.1h～10h；
航天器偏航轴的总干扰角加速度的理论值
Figure FDA0003038403660000047
ad3的估计值
Figure FDA0003038403660000048
即求得的z33为航天器偏航轴的总干扰角加速度的估计值。
6.根据权利要求1所述的一种抗慢变干扰的相平面姿态控制方法，其特征在于：
所述步骤(3)中，构建三轴的相平面，分别为：滚动轴相平面、俯仰轴相平面和偏航轴相平面，滚动轴相平面由滚动角x1和滚动角速度
Figure FDA00030384036600000521
构成，俯仰轴相平面由俯仰角x2和俯仰角速度
Figure FDA00030384036600000522
构成，偏航轴相平面由偏航角x3和偏航角速度
Figure FDA00030384036600000523
构成；
将滚动轴相平面分为七个区，七个区分别为：关机区S01，右侧大推力全开区S11，右侧小推力步进区S21，右侧小推力抗外滑区S31，左侧大推力全开区S11’，左侧小推力步进区S21’，左侧小推力抗外滑区S31’；
以x1＝0为纵坐标轴，
Figure FDA0003038403660000051
为横坐标轴；
所述S11区为直线
Figure FDA0003038403660000052
的右侧区域与直线
Figure FDA0003038403660000053
的上侧区域的重叠区域，以及直线
Figure FDA0003038403660000054
的上侧区域与直线x1＝xB的右侧区域的重叠区域，其中，直线
Figure FDA0003038403660000055
为右侧开关线，xD>0为死区阈值，
Figure FDA0003038403660000056
为所述小推力步进区的角速度最大值，K1>0为滚动轴相平面的开关线斜率绝对值，取值范围为
Figure FDA0003038403660000057
其中
Figure FDA0003038403660000058
为限速值，Tmin表示最小喷气长度；
所述S11’区为所述直线
Figure FDA0003038403660000059
的左侧区域与直线
Figure FDA00030384036600000510
下侧区域的重叠区域，以及直线
Figure FDA00030384036600000511
的下侧区域与直线x1＝-xB的左侧区域的重叠区域，其中直线
Figure FDA00030384036600000512
为左侧开关线；
所述S21区为直线
Figure FDA00030384036600000513
的右侧区域与直线
Figure FDA00030384036600000514
的上侧区域与直线
Figure FDA00030384036600000515
的下侧区域与直线x1＝xB的左侧区域的重叠区域，其中xB>0为所述大推力全开区的阈值；
所述S21’区为直线
Figure FDA00030384036600000516
的左侧区域与直线
Figure FDA00030384036600000517
下侧区域与直线
Figure FDA00030384036600000518
的上侧区域与直线x1＝-xB的右侧区域的重叠区域；
所述S31区为直线
Figure FDA00030384036600000519
的下侧区域与直线
Figure FDA00030384036600000520
的上侧区域与直线x1＝xD的右侧区域的重叠区域；
所述S31’区为直线
Figure FDA0003038403660000061
的上侧区域与直线
Figure FDA0003038403660000062
的下侧区域与直线x1＝-xD的左侧区域的重叠区域；
所述S01区为滚动轴相平面全区域除去所述的S11，S21，S31，S11’，S21’，S31’区以外的区域。
7.根据权利要求6所述的一种抗慢变干扰的相平面姿态控制方法，其特征在于：
将俯仰轴相平面分为七个区，七个区分别为：关机区S02，右侧大推力全开区S12，右侧小推力步进区S22，右侧小推力抗外滑区S32，左侧大推力全开区S12’，左侧小推力步进区S22’，左侧小推力抗外滑区S32’；
以x2＝0为纵坐标轴，
Figure FDA0003038403660000063
为横坐标轴；
所述S12区为直线
Figure FDA0003038403660000064
的右侧区域与直线
Figure FDA0003038403660000065
的上侧区域的重叠区域，以及直线
Figure FDA0003038403660000066
的上侧区域与直线x2＝xB的右侧区域的重叠区域，其中，直线
Figure FDA0003038403660000067
为右侧开关线，xD>0为死区阈值，
Figure FDA0003038403660000068
为所述小推力步进区的角速度最大值，K2>0为俯仰轴相平面的开关线斜率绝对值，取值范围为
Figure FDA0003038403660000069
其中
Figure FDA00030384036600000610
为限速值，Tmin表示最小喷气长度；
所述S12’区为所述直线
Figure FDA00030384036600000611
的左侧区域与直线
Figure FDA00030384036600000612
下侧区域的重叠区域，以及直线
Figure FDA00030384036600000613
的下侧区域与直线x2＝-xB的左侧区域的重叠区域，其中直线
Figure FDA00030384036600000614
为左侧开关线；
所述S22区为直线
Figure FDA00030384036600000615
的右侧区域与直线
Figure FDA00030384036600000616
的上侧区域与直线
Figure FDA00030384036600000617
的下侧区域与直线x2＝xB的左侧区域的重叠区域，其中xB>0为所述大推力全开区的阈值；
所述S22’区为直线
Figure FDA00030384036600000618
的左侧区域与直线
Figure FDA00030384036600000619
下侧区域与直线
Figure FDA00030384036600000620
的上侧区域与直线x2＝-xB的右侧区域的重叠区域；
所述S32区为直线
Figure FDA00030384036600000621
的下侧区域与直线
Figure FDA00030384036600000622
的上侧区域与直线x2＝xD的右侧区域的重叠区域；
所述S32’区为直线
Figure FDA0003038403660000071
的上侧区域与直线
Figure FDA0003038403660000072
的下侧区域与直线x2＝-xD的左侧区域的重叠区域；
所述S02区为俯仰轴相平面全区域除去所述的S12，S22，S32，S12’，S22’，S32’区以外的区域。
8.根据权利要求7所述的一种抗慢变干扰的相平面姿态控制方法，其特征在于：
将偏航轴相平面分为七个区，七个区分别为：关机区S03，右侧大推力全开区S13，右侧小推力步进区S23，右侧小推力抗外滑区S33，左侧大推力全开区S13’，左侧小推力步进区S23’，左侧小推力抗外滑区S33’；
以x3＝0为纵坐标轴，
Figure FDA0003038403660000073
为横坐标轴；
所述S13区为直线
Figure FDA0003038403660000074
的右侧区域与直线
Figure FDA0003038403660000075
的上侧区域的重叠区域，以及直线
Figure FDA0003038403660000076
的上侧区域与直线x3＝xB的右侧区域的重叠区域，其中，直线
Figure FDA0003038403660000077
为右侧开关线，xD>0为死区阈值，
Figure FDA0003038403660000078
为所述小推力步进区的角速度最大值，K3>0为偏航轴相平面的开关线斜率绝对值，取值范围为
Figure FDA0003038403660000079
其中
Figure FDA00030384036600000710
为限速值，Tmin表示最小喷气长度；
所述S13’区为所述直线
Figure FDA00030384036600000711
的左侧区域与直线
Figure FDA00030384036600000712
下侧区域的重叠区域，以及直线
Figure FDA00030384036600000713
的下侧区域与直线x3＝-xB的左侧区域的重叠区域，其中直线
Figure FDA00030384036600000714
为左侧开关线；
所述S23区为直线
Figure FDA00030384036600000715
的右侧区域与直线
Figure FDA00030384036600000716
的上侧区域与直线
Figure FDA00030384036600000717
的下侧区域与直线x3＝xB的左侧区域的重叠区域，其中xB>0为所述大推力全开区的阈值；
所述S23’区为直线
Figure FDA00030384036600000718
的左侧区域与直线
Figure FDA00030384036600000719
下侧区域与直线
Figure FDA00030384036600000720
的上侧区域与直线x3＝-xB的右侧区域的重叠区域；
所述S33区为直线
Figure FDA00030384036600000721
的下侧区域与直线
Figure FDA00030384036600000722
的上侧区域与直线x3＝xD的右侧区域的重叠区域；
所述S33’区为直线
Figure FDA0003038403660000081
的上侧区域与直线
Figure FDA0003038403660000082
的下侧区域与直线x3＝-xD的左侧区域的重叠区域；
所述S03区为偏航轴相平面全区域除去所述的S13，S23，S33，S13’，S23’，S33’区以外的区域。
9.根据权利要求8所述的一种抗慢变干扰的相平面姿态控制方法，其特征在于：
所述步骤(4)中，关机区的相平面控制指令为：
针对滚动轴，在当前采样时刻tk，tk＝t0+kΔT，k＝0,1,2,...表示采样时刻，其中t0表示第一次采样时刻，当
Figure FDA0003038403660000083
位于所述S01区时，则关闭产生绕滚动轴力矩的推力器，不进行喷气；
针对俯仰轴，在当前采样时刻tk，当
Figure FDA0003038403660000084
位于所述S02区时，则关闭产生绕俯仰轴力矩的推力器，不进行喷气；
针对偏航轴，在当前采样时刻tk，当
Figure FDA0003038403660000085
位于所述S03区时，则关闭产生绕偏航轴力矩的大推力器和小推力器，不进行喷气；
所述步骤(4)中，右侧大推力全开区的相平面控制指令为：
针对滚动轴，在当前采样时刻tk，当
Figure FDA0003038403660000086
位于所述S11区时，则采用大推力器喷气，产生绕滚动轴力矩为-FM1，喷气长度为ΔT，其中FM1为所配置的大推力器喷气产生的绕滚动轴的力矩，ΔT表示采样控制周期；
针对俯仰轴，在当前采样时刻tk，当
Figure FDA0003038403660000087
位于所述S12区时，则采用大推力器喷气，产生绕俯仰轴力矩为-FM2，喷气长度为ΔT，其中FM2为所配置的大推力器喷气产生的绕俯仰轴的力矩；
针对偏航轴，在当前采样时刻tk，当
Figure FDA0003038403660000088
位于所述S13区时，则采用大推力器喷气，产生绕偏航轴力矩为-FM3，喷气长度为ΔT，其中FM3为所配置的大推力器喷气产生的绕偏航轴的力矩；
所述步骤(4)中，左侧大推力全开区的相平面控制指令为：
针对滚动轴，在当前采样时刻tk，当
Figure FDA0003038403660000091
位于所述S11’区时，则采用大推力器喷气，产生绕滚动轴力矩为FM1，喷气长度为ΔT；
针对俯仰轴，在当前采样时刻tk，当
Figure FDA0003038403660000092
位于所述S12’区时，则采用大推力器喷气，产生绕俯仰轴力矩为FM2，喷气长度为ΔT；
针对偏航轴，在当前采样时刻tk，当
Figure FDA0003038403660000093
位于所述S13’区时，则采用大推力器喷气，产生绕偏航轴力矩为FM3，喷气长度为ΔT；
所述步骤(4)中，右侧小推力步进区的相平面控制指令为：
(i)针对滚动轴，在当前采样时刻tk，当
Figure FDA0003038403660000094
位于所述S21：
若满足
Figure FDA0003038403660000095
且
Figure FDA0003038403660000096
时，则本周期关闭产生绕滚动轴力矩的小推力器和大推力器，不进行喷气，其中ε1为滚动轴干扰阈值，取值范围为(0,0.01Fm1/I11]；
若满足
Figure FDA0003038403660000097
且
Figure FDA0003038403660000098
时，则本周期采用小推力器喷气，产生绕滚动轴力矩为-Fm1，喷气长度TN为：TN＝TN1+TN2，
Figure FDA0003038403660000099
TN2∈[Tmin,0.8ΔT]，TN∈[Tmin,ΔT]，其中K11>0为滚动轴相平面小推力步进区干扰补偿系数，
Figure FDA00030384036600000910
为小推力角加速度参数，k2∈[0.5,0.9]为常数；
若满足
Figure FDA00030384036600000911
且
Figure FDA00030384036600000912
时，则本周期采用小推力器喷气，产生绕滚动轴力矩为-Fm1，喷气长度为Tmin；
若满足
Figure FDA00030384036600000913
且
Figure FDA00030384036600000914
时，则本周期采用小推力器喷气，产生绕滚动轴力矩为-Fm1，喷气长度TN为：TN＝TN1+TN2，
Figure FDA00030384036600000915
TN2＝Kj(|x1|-xe)/ac2，TN2∈[Tmin,0.5ΔT]，TN∈[Tmin,ΔT]，其中Kj＝((1-k2)ac2ΔT)/(xB-xe)为小推力步进区参数，xe>0为略小于xD的常数；
(ii)针对俯仰轴，在当前采样时刻tk，当
Figure FDA00030384036600000916
位于所述S22：
若满足
Figure FDA00030384036600000917
且
Figure FDA00030384036600000918
时，则本周期关闭产生绕俯仰轴力矩的小推力器和大推力器，不进行喷气，其中ε2为俯仰轴干扰阈值，取值范围为(0,0.01Fm2/I22]；
若满足
Figure FDA0003038403660000101
且
Figure FDA0003038403660000102
时，则本周期采用小推力器喷气，产生绕俯仰轴力矩为-Fm2，喷气长度TN为：TN＝TN1+TN2，
Figure FDA0003038403660000103
Figure FDA0003038403660000104
TN2∈[Tmin,0.8ΔT]，TN∈[Tmin,ΔT]，其中K12>0为俯仰轴相平面小推力步进区干扰补偿系数；
若满足
Figure FDA0003038403660000105
且
Figure FDA0003038403660000106
时，则本周期采用小推力器喷气，产生绕俯仰轴力矩为-Fm2，喷气长度为Tmin；
若满足
Figure FDA0003038403660000107
且
Figure FDA0003038403660000108
时，则本周期采用小推力器喷气，产生绕俯仰轴力矩为-Fm2，喷气长度TN为：TN＝TN1+TN2，
Figure FDA0003038403660000109
TN2＝Kj(|x2|-xe)/ac2，TN2∈[Tmin,0.5ΔT]，TN∈[Tmin,ΔT]；
(iii)针对偏航轴，在当前采样时刻tk，当
Figure FDA00030384036600001010
位于所述S23：
若满足
Figure FDA00030384036600001011
且
Figure FDA00030384036600001012
时，则本周期关闭产生绕偏航轴力矩的小推力器和大推力器，不进行喷气，其中ε3为偏航轴干扰阈值，取值范围为(0,0.01Fm3/I33]；
若满足
Figure FDA00030384036600001013
且
Figure FDA00030384036600001014
时，则本周期采用小推力器喷气，产生绕偏航轴力矩为-Fm3，喷气长度TN为：TN＝TN1+TN2，
Figure FDA00030384036600001015
TN2∈[Tmin,0.8ΔT]，TN∈[Tmin,ΔT]，其中K13>0为俯仰轴相平面小推力步进区干扰补偿系数；
若满足
Figure FDA00030384036600001016
且
Figure FDA00030384036600001017
时，则本周期采用小推力器喷气，产生绕偏航轴力矩为-Fm3，喷气长度为Tmin；
若满足
Figure FDA00030384036600001018
且
Figure FDA00030384036600001019
时，则本周期采用小推力器喷气，产生绕偏航轴力矩为-Fm3，喷气长度TN为：TN＝TN1+TN2，
Figure FDA00030384036600001020
TN2＝Kj(|x3|-xe)/ac2，TN2∈[Tmin,0.5ΔT]，TN∈[Tmin,ΔT]；
所述步骤(4)中，左侧小推力步进区的相平面控制指令为：
(i’)针对滚动轴，在当前采样时刻tk，当
Figure FDA0003038403660000111
位于所述S21’:
若满足
Figure FDA0003038403660000112
且
Figure FDA0003038403660000113
时，则本周期关闭产生绕滚动轴力矩的小推力器和大推力器，不进行喷气；
若满足
Figure FDA0003038403660000114
且
Figure FDA0003038403660000115
时，则本周期采用小推力器喷气，产生绕滚动轴力矩为Fm1，喷气长度TN＝TN1+TN2，
Figure FDA0003038403660000116
若满足
Figure FDA0003038403660000117
且
Figure FDA0003038403660000118
时，则本周期采用小推力器喷气，产生绕滚动轴力矩为Fm1，喷气长度为Tmin；
若满足
Figure FDA0003038403660000119
且
Figure FDA00030384036600001110
时，则本周期采用小推力器喷气，产生绕滚动轴力矩为Fm1，喷气长度TN为：TN＝TN1+TN2，
Figure FDA00030384036600001111
TN2＝Kj(|x1|-xe)/ac2，TN2∈[Tmin,0.5ΔT]，TN∈[Tmin,ΔT]；
(ii’)针对俯仰轴，在当前采样时刻tk，当
Figure FDA00030384036600001112
位于所述S22’：
若满足
Figure FDA00030384036600001113
且
Figure FDA00030384036600001114
时，则本周期关闭产生绕俯仰轴力矩的小推力器和大推力器，不进行喷气；
若满足
Figure FDA00030384036600001115
且
Figure FDA00030384036600001116
时，则本周期采用小推力器喷气，产生绕俯仰轴力矩为Fm2，喷气长度TN＝TN1+TN2，
Figure FDA00030384036600001117
若满足
Figure FDA00030384036600001118
且
Figure FDA00030384036600001119
时，则本周期采用小推力器喷气，产生绕俯仰轴力矩为Fm2，喷气长度为Tmin；
若满足
Figure FDA00030384036600001120
且
Figure FDA00030384036600001121
时，则本周期采用小推力器喷气，产生绕俯仰轴力矩为Fm2，喷气长度TN为：TN＝TN1+TN2，
Figure FDA00030384036600001122
TN2＝Kj(|x2|-xe)/ac2，TN2∈[Tmin,0.5ΔT]，TN∈[Tmin,ΔT]；
(iii’)针对偏航轴，在当前采样时刻tk，当
Figure FDA00030384036600001123
位于所述S23’：
若满足
Figure FDA00030384036600001124
且
Figure FDA00030384036600001125
时，则本周期关闭产生绕偏航轴力矩的小推力器和大推力器，不进行喷气；
若满足
Figure FDA00030384036600001126
且
Figure FDA00030384036600001127
时，则本周期采用小推力器喷气，产生绕偏航轴力矩为Fm3，喷气长度TN＝TN1+TN2，
Figure FDA0003038403660000121
若满足
Figure FDA0003038403660000122
且
Figure FDA0003038403660000123
时，则本周期采用小推力器喷气，产生绕偏航轴力矩为Fm3，喷气长度为Tmin；
若满足
Figure FDA0003038403660000124
且
Figure FDA0003038403660000125
时，则本周期采用小推力器喷气，产生绕偏航轴力矩为Fm3，喷气长度TN为：TN＝TN1+TN2，
Figure FDA0003038403660000126
TN2＝Kj(|x3|-xe)/ac2，TN2∈[Tmin,0.5ΔT]，TN∈[Tmin,ΔT]；
所述步骤(4)中，右侧小推力抗外滑区的相平面控制指令为：
(I)针对滚动轴，在当前采样时刻tk，当
Figure FDA0003038403660000127
位于所述S31：
若满足
Figure FDA0003038403660000128
时，则本周期关闭产生绕滚动轴力矩的小推力器和大推力器，不进行喷气；
若满足
Figure FDA0003038403660000129
时，则本周期采用小推力器喷气，产生绕滚动轴力矩为-Fm1，喷气长度TN为：TN＝TN1+TN2，TN1＝Tmin，
Figure FDA00030384036600001210
TN2∈[0,0.8ΔT]，TN∈[Tmin,ΔT]，其中K21为滚动轴小推力抗外滑区干扰补偿系数，满足0<K21<K11；
(II)针对俯仰轴，在当前采样时刻tk，当
Figure FDA00030384036600001211
位于所述S32：
若满足
Figure FDA00030384036600001212
时，则本周期关闭产生绕俯仰轴力矩的小推力器和大推力器，不进行喷气；
若满足
Figure FDA00030384036600001213
时，则本周期采用小推力器喷气，产生绕俯仰轴力矩为-Fm2，喷气长度TN为：TN＝TN1+TN2，TN1＝Tmin，
Figure FDA00030384036600001214
TN2∈[0,0.8ΔT]，TN∈[Tmin,ΔT]，其中K22为俯仰轴小推力抗外滑区干扰补偿系数，满足0<K22<K12；
(III)针对偏航轴，在当前采样时刻tk，当
Figure FDA00030384036600001215
位于所述S33：
若满足
Figure FDA00030384036600001216
时，则本周期关闭产生绕偏航轴力矩的小推力器和大推力器，不进行喷气；
若满足
Figure FDA0003038403660000131
时，则本周期采用小推力器喷气，产生绕偏航轴力矩为-Fm3，喷气长度TN为：TN＝TN1+TN2，TN1＝Tmin，
Figure FDA0003038403660000132
TN2∈[0,0.8ΔT]，TN∈[Tmin,ΔT]，其中K23为偏航轴小推力抗外滑区干扰补偿系数，满足0<K23<K13；
所述步骤(4)中，左侧小推力抗外滑区的相平面控制指令为：
(I’)针对滚动轴，在当前采样时刻tk，当
Figure FDA0003038403660000133
位于所述S31’：
若满足
Figure FDA0003038403660000134
时，则本周期关闭产生绕滚动轴力矩的小推力器和大推力器，不进行喷气；
若满足
Figure FDA0003038403660000135
时，则本周期采用小推力器喷气，产生绕滚动轴力矩为-Fm1，喷气长度TN：TN＝TN1+TN2，TN1＝Tmin，
Figure FDA0003038403660000136
TN2∈[0,0.8ΔT]，TN∈[Tmin,ΔT]；
(II’)针对俯仰轴，在当前采样时刻tk，当
Figure FDA0003038403660000137
位于所述S32’：
若满足
Figure FDA0003038403660000138
时，则本周期关闭产生绕俯仰轴力矩的小推力器和大推力器，不进行喷气；
若满足
Figure FDA0003038403660000139
时，则本周期采用小推力器喷气，产生绕俯仰轴力矩为-Fm2，喷气长度TN：TN＝TN1+TN2，TN1＝Tmin，
Figure FDA00030384036600001310
TN2∈[0,0.8ΔT]，TN∈[Tmin,ΔT]；
(III’)针对偏航轴，在当前采样时刻tk，当
Figure FDA00030384036600001311
位于所述S33’：
若满足
Figure FDA00030384036600001312
时，则本周期关闭产生绕偏航轴力矩的小推力器和大推力器，不进行喷气；
若满足
Figure FDA00030384036600001313
时，则本周期采用小推力器喷气，产生绕偏航轴力矩为-Fm3，喷气长度TN：TN＝TN1+TN2，TN1＝Tmin，
Figure FDA00030384036600001314
TN2∈[0,0.8ΔT]，TN∈[Tmin,ΔT]；
所述步骤(4)中，在计算右侧小推力步进区和左侧小推力步进区的相平面控制指令中的三轴相平面控制喷气长度之前，对三轴干扰角加速度估计进行限幅，具体限幅方法为：
针对滚动轴，当
Figure FDA0003038403660000141
时，则取
Figure FDA0003038403660000142
当
Figure FDA0003038403660000143
时，则取
Figure FDA0003038403660000144
其中限幅值Lmt1取0.1Fm1/I11～0.5Fm1/I11；
针对俯仰轴，当
Figure FDA0003038403660000145
时，则取
Figure FDA0003038403660000146
当
Figure FDA0003038403660000147
时，则取
Figure FDA0003038403660000148
其中限幅值Lmt2取0.1Fm2/I22～0.5Fm2/I22；
针对偏航轴，当
Figure FDA0003038403660000149
时，则取
Figure FDA00030384036600001410
当
Figure FDA00030384036600001411
时，则取
Figure FDA00030384036600001412
其中限幅值Lmt3取0.1Fm3/I33～0.5Fm3/I33；
在根据航天器三轴姿态角测量输出值设计所述扩张状态观测器之前，还对航天器三轴姿态角测量输出值进行了惯性滤波处理，具体滤波方法为：
针对滚动轴，设计惯性滤波形式为：
Figure FDA00030384036600001413
其中，参数α1>0表示滚动轴姿态滤波系数，yf1为滚动姿态角滤波输出值；
针对俯仰轴，设计惯性滤波形式为：
Figure FDA00030384036600001414
其中，参数α2>0表示俯仰轴姿态滤波系数，yf2为俯仰姿态角滤波输出值；
针对偏航轴，设计惯性滤波形式为：
Figure FDA00030384036600001415
其中，参数α3>0表示偏航轴姿态滤波系数，yf3为偏航姿态角滤波输出值。
10.根据权利要求1所述的一种抗慢变干扰的相平面姿态控制方法，其特征在于：
所述步骤(5)中，航天器的推力器根据步骤(4)得到的三轴相平面控制指令进行喷气，具体为：
针对航天器的滚动轴，在当前采样时刻tk和下一个采样时刻tk+ΔT之间所构成的采样控制周期[tk,tk+ΔT]内，在当前采样时刻将滚动轴相平面控制指令中所采用的推力器开机，开机时间长度等于滚动轴相平面控制指令中的滚动轴喷气长度，然后将所采用的推力器关机；
针对航天器的俯仰轴，在当前采样时刻tk和下一个采样时刻tk+ΔT之间所构成的采样控制周期[tk,tk+ΔT]内，在当前采样时刻将俯仰轴相平面控制指令中所采用的推力器开机，开机时间长度等于俯仰轴相平面控制指令中的俯仰轴喷气长度，然后将所采用的推力器关机；
针对航天器的偏航轴，在当前采样时刻tk和下一个采样时刻tk+ΔT之间所构成的采样控制周期[tk,tk+ΔT]内，在当前采样时刻将偏航轴相平面控制指令中所采用的推力器开机，开机时间长度等于偏航轴相平面控制指令中的偏航轴喷气长度，然后将所采用的推力器关机。
CN202110450138.1A 2021-04-25 2021-04-25 一种抗慢变干扰的相平面姿态控制方法 Active CN113311850B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110450138.1A CN113311850B (zh) 2021-04-25 2021-04-25 一种抗慢变干扰的相平面姿态控制方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110450138.1A CN113311850B (zh) 2021-04-25 2021-04-25 一种抗慢变干扰的相平面姿态控制方法

Publications (2)

Publication Number Publication Date
CN113311850A true CN113311850A (zh) 2021-08-27
CN113311850B CN113311850B (zh) 2024-03-19

Family

ID=77371023

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110450138.1A Active CN113311850B (zh) 2021-04-25 2021-04-25 一种抗慢变干扰的相平面姿态控制方法

Country Status (1)

Country Link
CN (1) CN113311850B (zh)

Citations (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2008143398A (ja) * 2006-12-12 2008-06-26 Toshiba Corp 飛翔体制御システムおよび飛翔体の飛行制御方法
CN105468007A (zh) * 2015-12-02 2016-04-06 哈尔滨工业大学 一种基于干扰观测器的挠性卫星轨迹线性化姿态控制方法
CN107450588A (zh) * 2017-09-29 2017-12-08 北京航空航天大学 一种挠性航天器姿态控制系统的强抗扰控制方法
CN108459610A (zh) * 2018-02-27 2018-08-28 北京控制工程研究所 一种着陆器动力下降液体晃动抑制方法
CN108501944A (zh) * 2018-05-14 2018-09-07 吕杉 汽车爆胎安全稳定控制方法
CN109856968A (zh) * 2018-09-19 2019-06-07 湘潭大学 一种基于相平面模糊自调整的倒立摆系统改进型自抗扰控制方法
CN110058520A (zh) * 2019-04-02 2019-07-26 清华大学 一种固定时间收敛输出反馈模型参考控制方法
CN110850887A (zh) * 2019-10-10 2020-02-28 南京航空航天大学 一种四旋翼无人机复合动态逆抗干扰姿态控制方法
CN112286042A (zh) * 2020-09-14 2021-01-29 西安工业大学 基于扩张状态观测器的光电稳瞄平台三步线性控制方法

Patent Citations (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2008143398A (ja) * 2006-12-12 2008-06-26 Toshiba Corp 飛翔体制御システムおよび飛翔体の飛行制御方法
CN105468007A (zh) * 2015-12-02 2016-04-06 哈尔滨工业大学 一种基于干扰观测器的挠性卫星轨迹线性化姿态控制方法
CN107450588A (zh) * 2017-09-29 2017-12-08 北京航空航天大学 一种挠性航天器姿态控制系统的强抗扰控制方法
CN108459610A (zh) * 2018-02-27 2018-08-28 北京控制工程研究所 一种着陆器动力下降液体晃动抑制方法
CN108501944A (zh) * 2018-05-14 2018-09-07 吕杉 汽车爆胎安全稳定控制方法
CN109856968A (zh) * 2018-09-19 2019-06-07 湘潭大学 一种基于相平面模糊自调整的倒立摆系统改进型自抗扰控制方法
CN110058520A (zh) * 2019-04-02 2019-07-26 清华大学 一种固定时间收敛输出反馈模型参考控制方法
CN110850887A (zh) * 2019-10-10 2020-02-28 南京航空航天大学 一种四旋翼无人机复合动态逆抗干扰姿态控制方法
CN112286042A (zh) * 2020-09-14 2021-01-29 西安工业大学 基于扩张状态观测器的光电稳瞄平台三步线性控制方法

Non-Patent Citations (5)

* Cited by examiner, † Cited by third party
Title
CARLO ZANONI 等: "《Estimation of the Electrostatic Effects in the LISA-Pathfinder Critical Test Mass Dynamics via the Method of Moments》", 《IEEE/ASME TRANSACTIONS ON MECHATRONICS》, pages 1 - 9 *
HONGHUA ZHANG 等: "《Attitudecontrolandsloshingsuppressionforliquid-filled spacecraftinthepresenceofsinusoidaldisturbance》", 《JOURNAL OFSOUNDANDVIBRATION》, pages 1 - 12 *
席露华 等: "《精确轨道航天器的姿态机动策略选择》", 《飞行器测控学报》, vol. 35, no. 1, pages 1 - 5 *
詹韬 等: "《基于ESO的再入飞行器姿态控制》", 《导航定位与授时》, vol. 1, no. 3, pages 1 - 7 *
郭振西 等: "《运载火箭无滚控装置的制导控制技术研究》", 《导弹与航天运载技术》, no. 343, pages 1 - 5 *

Also Published As

Publication number Publication date
CN113311850B (zh) 2024-03-19

Similar Documents

Publication Publication Date Title
CN108803649B (zh) 一种垂直起降重复使用运载器自抗扰滑模控制方法
CN109189085B (zh) 基于事件触发的航天器网络化系统姿态控制方法
CN109062052B (zh) 基于扩张状态观测器的四旋翼无人机积分滑模控制方法
CN107688295A (zh) 一种基于快速终端滑模的四旋翼飞行器有限时间自适应控制方法
CN110850887B (zh) 一种四旋翼无人机复合动态逆抗干扰姿态控制方法
CN109542112B (zh) 一种针对垂直起降可重复使用火箭返回飞行的固定时间收敛抗扰控制方法
CN112445234B (zh) 一种航天器的姿态控制方法和装置
CN113306747B (zh) 基于so(3)群的挠性航天器姿态稳定控制方法和系统
CN107402516B (zh) 基于联合执行机构的递阶饱和模糊pd姿态控制方法
CN104881553B (zh) 单滑块滚喷模式变质心飞行器模型及其结构布局参数的设计方法
CN111831002B (zh) 一种基于预设性能的超高声速飞行器姿态控制方法
CN112486193B (zh) 一种基于自适应增广控制理论的飞翼无人机三轴全权限控制方法
Huang et al. Finite-time control for a UAV system based on finite-time disturbance observer
CN114281092A (zh) 一种基于滑模干扰观测器的高超声速飞行器协调姿态控制方法
CN111258216A (zh) 一种适用于四旋翼飞行器的滑模重复控制器
CN109625329B (zh) 一种基于离散喷气的飞轮角动量自主卸载方法
CN113885543B (zh) 一种复杂海况干扰下的航行体自抗扰横滚控制器的设计方法
CN113311850B (zh) 一种抗慢变干扰的相平面姿态控制方法
CN111007867B (zh) 一种可预设调整时间的高超声速飞行器姿态控制设计方法
CN108459611A (zh) 一种近空间飞行器的姿态跟踪控制方法
CN116430828A (zh) 一种基于观测器的四旋翼故障容灾降级控制方法
CN114706413B (zh) 近地轨道微纳卫星变质心姿态控制方法及系统
CN116360258A (zh) 基于固定时间收敛的高超声速变形飞行器抗干扰控制方法
CN113238572B (zh) 基于预设性能控制的预设时间四旋翼无人机姿态跟踪方法
CN113534842B (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