CN115857529A - 航天器姿态控制系统的执行器故障重构方法 - Google Patents

航天器姿态控制系统的执行器故障重构方法 Download PDF

Info

Publication number
CN115857529A
CN115857529A CN202211550334.7A CN202211550334A CN115857529A CN 115857529 A CN115857529 A CN 115857529A CN 202211550334 A CN202211550334 A CN 202211550334A CN 115857529 A CN115857529 A CN 115857529A
Authority
CN
China
Prior art keywords
spacecraft
fault
state
time
matrix
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
CN202211550334.7A
Other languages
English (en)
Other versions
CN115857529B (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.)
Harbin Institute of Technology
Original Assignee
Harbin Institute of Technology
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 Harbin Institute of Technology filed Critical Harbin Institute of Technology
Priority to CN202211550334.7A priority Critical patent/CN115857529B/zh
Publication of CN115857529A publication Critical patent/CN115857529A/zh
Application granted granted Critical
Publication of CN115857529B publication Critical patent/CN115857529B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
    • Y02T90/00—Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation

Landscapes

  • Control Of Position, Course, Altitude, Or Attitude Of Moving Bodies (AREA)

Abstract

航天器姿态控制系统的执行器故障重构方法,属于航空航天飞行控制及故障诊断领域。解决了现有航天器姿态控制系统执行器故障诊断、辨识准确率低的问题。本发明建立航天器姿态控制系统的运动学模型和带有模型不确定性、空间环境扰动、执行机构故障、传感器故障的航天器姿态控制系统动力学模型；设计滑模控制器，获取航天器姿态的控制力矩；对扰动和故障参数进行状态增广和估计，获取满足卡尔曼滤波理论的状态空间方程模型；采用双重自适应强跟踪平方根容积卡尔曼滤波算法和LSTM神经网络，构建FDIR框架；将航天器姿态的控制力矩和测量的航天器的角速度输入至FDIR框架，对故障参数进行估计。本发明适用于航天器故障估计。

Description

航天器姿态控制系统的执行器故障重构方法
技术领域
本发明属于航空航天飞行控制及故障诊断领域。
背景技术
由于航天器在恶劣的太空环境中长期运行，始终受到空间环境的各种干扰的作用以及高低温交替、电磁辐射、空间粒子照射等多方面因素，使得航天器各分系统、元器件的可靠度降低，从而难免会出现各种故障，造成不可估量且难以挽回的损失。姿态控制系统(Attitude Control System,ACS)是航天器中最关键的子系统之一，其安全性和可靠性对于确保任务的顺利完成至关重要。航天器姿态控制系统的故障一般是由其部件、执行器和传感器的故障引起的，如姿控发动机、反作用轮、陀螺仪和星敏感器。目前，航天器故障诊断通常依靠地面站遥测数据阈值告警，速度慢，准确率低。随着航天器变得越来越复杂和任务时间延长，对自主故障诊断的要求越来越高。所以为传感器、执行器或系统组件开发一种智能的FDIR框架是必要的，也是实现航天器容错控制(FTC)和健康管理系统(PHM)的重要前提。
故障检测、隔离和重构(Fault Detection,Isolation,and Reconstruction,FDIR)，FDIR方案旨在检测和隔离故障的存在和位置，准确估计故障的大小，FDIR框架旨在检测和隔离故障的存在和位置，准确估计故障的大小。目前，解决FDIR问题的方法主要有两种，即基于模型的故障位置估计方法和基于数据的故障位置估计方法。基于数据的故障位置估计，方法需要布置大量的传感器并收集全生命周期的数据。然而，由于航天器的特殊性，通常很难安装各种传感器并收集大量各种类型的故障数据。基于模型的方法充分利用了系统内现有的少量传感器和已知的解析冗余，具有更可靠、易于实现、节省设备成本和空间等优点。
由于基于模型的故障估计方法要求系统模型绝对精确，但实际模型中包含着空间扰动、未建模动态、系统参数摄动等因素，制约着高性能FDIR的设计，现有基于模型的FDIR方法存在四个方面的不足，仍然面临着诸多挑战。第一，如何解决执行器故障种类较多且随机不确定，且多个执行器同时故障的诊断问题。第二，如何解决模型不确定性对传统基于模型方法精度的影响。第三，如何提升FDIR的速度，以满足故障诊断的实时性要求。第四，基于模型的方法通常假设测量信号是准确的，但传感器同样有可能发生故障，出现测量不确定性，如何保证测量信号不准确时的FDIR结果仍然满足性能要求是当今的难题。
发明内容
本发明目的是为了解决现有航天器姿态控制系统执行器故障诊断、辨识准确率低的问题，提供了航天器姿态控制系统的执行器故障重构方法。
本发明所述航天器姿态控制系统的执行器故障重构方法，其特征在于，具体过程为：
步骤一、建立航天器姿态控制系统的运动学模型和带有模型不确定性、空间环境扰动、执行机构故障、传感器故障的航天器姿态控制系统动力学模型；
步骤二、利用航天器姿态控制系统的运动学模型和所述航天器姿态控制系统动力学模型，设计滑模控制器，获取航天器姿态的控制力矩；
对所述航天器姿态控制系统动力学模型进行离散化，同时对扰动和故障参数进行状态增广和估计，获取满足卡尔曼滤波理论的状态空间方程模型；
步骤三、基于满足卡尔曼滤波理论的状态空间方程模型，采用双重自适应强跟踪平方根容积卡尔曼滤波算法和LSTM神经网络，构建FDIR框架；
步骤四、将步骤二获取航天器姿态的控制力矩和测量的航天器的角速度输入至FDIR框架，对故障参数进行估计。
进一步地，本发明中，步骤一中，航天器姿态控制系统的运动学模型为：
Figure BDA0003980746460000021
ω＝ωbo+Cboωoi+b+η
其中，q＝[q0,qv]T为姿态四元数，描述航天器自身固定坐标系
Figure BDA0003980746460000022
相对于惯性坐标系FI的方向，满足||q||2＝1，其中q0为q的标量部分，qv为q的矢量部分；ω为航天器自身固定坐标系/>
Figure BDA0003980746460000023
相对于惯性坐标系/>
Figure BDA0003980746460000024
的角速度，ω＝[ω1,ω2,ω3]T，其中ω1,ω2,ω3分别代表航天器在自身坐标系中三轴的角速度；Cbo表示航天器自身固定坐标系/>
Figure BDA0003980746460000025
到轨道坐标系/>
Figure BDA0003980746460000026
对应的旋转矩阵，/>
Figure BDA0003980746460000027
ωbo为航天器自身固定坐标系/>
Figure BDA0003980746460000028
相对于轨道坐标系/>
Figure BDA0003980746460000029
的角速度，ωoi表示轨道坐标系/>
Figure BDA00039807464600000210
相对于惯性坐标系/>
Figure BDA00039807464600000211
的角速度；b和η分别是角速率陀螺漂移和测量噪声，上角T代表矩阵的转置；I3为3×3的单位矩阵，/>
Figure BDA00039807464600000212
是q的一阶导数，/>
Figure BDA00039807464600000213
为qv的叉乘矩阵。
进一步地，本发明中，步骤一中，带有模型不确定性、空间环境扰动、执行机构故障、传感器故障的航天器姿态控制系统动力学模型为：
Figure BDA00039807464600000214
其中J0+ΔJ∈R3×3，J0+ΔJ表示实际转动惯量矩阵，包括标称部分J0和不确定部分ΔJ；通过进一步地整理及细化得到：
Figure BDA0003980746460000031
Figure BDA0003980746460000032
Figure BDA0003980746460000033
J0+ΔJ＝Ja
Figure BDA0003980746460000034
其中，D0是航天器执行器的标称安装矩阵；ΔD代表航天器执行器实际安装与标称安装之间的误差矩阵；τ＝[τ1,τ2,…τa]T∈Ra×1是执行器的输出力矩，τi是第i个执行器的输出力矩，其中i＝1,2,…a，a表示航天器姿态控制系统执行器个数；uc为待设计的控制器输出力矩；sat(·)表示饱和函数；Λ＝diag[Λ1,…,Λa]∈Ra×a是执行器健康因子矩阵，并满足0≤Λi≤1，Λi＝1表示第i个执行器完全健康，0＜Λi＜1表示执行器部分失效，Λi＝0表示完全失效故障；Td为环境扰动力矩；d为包含不确定性和内扰动以及环境扰动力矩在内的总扰动，
Figure BDA0003980746460000035
是ω的一阶导数，/>
Figure BDA0003980746460000036
是中间变量，ω×为ω的叉乘矩阵。
进一步地，本发明中，步骤二中，获取的航天器姿态的控制力矩的过程为：
步骤二一、设航天器姿态的姿态角控制目标为qd，姿态角控制目标为ωd；
利用航天器姿态的姿态角控制目标为qd，计算姿态误差四元数qe＝[qe0,qev]T，qe0为误差四元数标量部分，qev为误差四元数矢量部分，计算公式为：
Figure BDA0003980746460000037
为四元数乘法运算；
利用姿态角控制目标为ωd，计算角速度误差ωe＝[ωe1,ωe2,ωe3]T，ωe1,ωe2,ωe3分别为航天器在自身坐标系中三轴的角速度误差，计算公式为：ωe＝ω-ωd；
步骤二二、利用角速度误差和姿态误差四元数设计滑模控制器；利用滑膜控制器获取航天器姿态的控制力矩，具体为：
Figure BDA0003980746460000038
Figure BDA0003980746460000039
Figure BDA00039807464600000310
其中k1,k2,k3是三个正标量，s为滑模面，s＝ωe+k1qev，sat(·)为饱和函数，sgn(·)为符号函数，τmax是单个执行器的最大输出力矩，τi为第i个执行机构的输出力矩，
Figure BDA0003980746460000041
是qev的一阶导数，其中/>
Figure BDA0003980746460000042
为ωe的叉乘矩阵，us是中间变量。
进一步地，本发明中，步骤二中，满足卡尔曼滤波理论的状态空间方程模型包括：扰动观测器模型和故障观测器模型；
所述扰动观测器模型：
Figure BDA0003980746460000043
扰动观测器模型中，前两个是状态方程，最后一个是观测方程；
所述故障观测器模型：
Figure BDA0003980746460000044
故障观测器模型中，前两个是状态方程，最后一个是观测方程；
式中，ωk+1,ωk分别代表航天器自身固定坐标系
Figure BDA0003980746460000045
相对于惯性坐标系/>
Figure BDA0003980746460000046
的角速度ω在第k+1和k时刻的值；/>
Figure BDA0003980746460000047
Figure BDA0003980746460000048
代表航天器自身固定坐标系/>
Figure BDA0003980746460000049
相对于惯性坐标系/>
Figure BDA00039807464600000410
的角速度ω的过程噪声；yk代表在第k时刻的航天器自身固定坐标系/>
Figure BDA00039807464600000411
相对于惯性坐标系/>
Figure BDA00039807464600000412
的角速度的测量值；vk代表航天器自身固定坐标系/>
Figure BDA00039807464600000413
相对于惯性坐标系/>
Figure BDA00039807464600000414
的角速度的测量噪声，在传统的卡尔曼滤波当中，/>
Figure BDA00039807464600000415
需要满足统计特性恒定不变的白噪声序列，τk为执行器输出控制力矩τ在k时刻的值，T离散化采样时间；Λk+1表示k+1时刻的健康因子矩阵，Λk表示k时刻的健康因子矩阵，/>
Figure BDA00039807464600000416
是健康因子矩阵的随机有走噪声，dk+1和dk分别是包含不确定性和内外扰动在内的总扰动在k+1时刻和k时刻的值；/>
Figure BDA00039807464600000417
为k时刻总扰动dk的随机游走噪声。/>
Figure BDA00039807464600000418
分别为k时刻需要估计并替换补偿的健康因子矩阵和总扰动。
进一步地，本发明中，步骤三、构建FDIR框架的过程为：
步骤三一，利用扰动观测器模型，采用双重自适应强跟踪平方根容积卡尔曼滤波算法，对航天器系统初始运行h个时刻的总扰动进行状态估计，获得扰动估计值序列
Figure BDA0003980746460000051
其中h表示系统初始运行时间结束时刻；
步骤三二、采用LSTM神经网络利用扰动估计值序列
Figure BDA0003980746460000052
预测航天器系统未来f个时刻的扰动状态，获得预测值序列/>
Figure BDA0003980746460000053
其中f＜h/2；
步骤三三、利用预测值序列
Figure BDA0003980746460000054
对故障观测器模型中当前运行m个时刻获得的总扰动[dk]m进行替换，其中，m小于或等于f；其具体替换方法为按照时间顺序一一对应进行替换；
步骤三四、故障观测器模型采用双重自适应强跟踪平方根容积卡尔曼滤波算法对系统的健康因子状态进行估计得到
Figure BDA0003980746460000055
(此时的/>
Figure BDA0003980746460000056
为故障重构结果，)直至获得p个时刻的故障估计值：/>
Figure BDA0003980746460000057
并判断/>
Figure BDA0003980746460000058
中故障估计值随时间增加是否趋于平稳，若是，将扰动观测器模型中的Λk替换为/>
Figure BDA0003980746460000059
完成一次FDIR框架的构建，否则，执行步骤三五；
步骤三五、令扰动观测器模型采用双重自适应强跟踪平方根容积卡尔曼滤波算法再次对系统运行h个时刻的扰动进行状态估计，获得扰动估计值序列
Figure BDA00039807464600000510
返回执行步骤三二，其中，r为返回次数。
进一步地，本发明中，步骤四中，对故障参数进行估计的过程为：
步骤四一、令
Figure BDA00039807464600000511
作为状态空间方程变量，对所述状态空间方程变量的初始状态估计值和状态估计误差协方差矩阵进行初始化；
步骤四二、计算状态估计误差协方差矩阵的平方根Sk-1|k-1：
步骤四三、利用公式：
Figure BDA00039807464600000512
/>
构造k-1时刻容积点、采用扰动观测器模型和故障观测器模型中的状态方程传播容积点，获得k时刻状态变量的预测值；利用k时刻的状态预测值，计算状态误差协方差预测矩阵的平方根；
其中，Xi,k-1|k-1表示第k-1时刻的最优估计值对应的容积点，Sk-1|k-1为k-1时刻状态误差协方差预测矩阵的平方根，
Figure BDA0003980746460000061
为第k-1时刻的状态预测值；ζi表示第i个容积点对应的权重，/>
Figure BDA0003980746460000062
n为xk的维数，所述xk为/>
Figure BDA0003980746460000063
或/>
Figure BDA0003980746460000064
所述/>
Figure BDA0003980746460000065
和/>
Figure BDA0003980746460000066
维数相同，[l]i是点集[l]的第i列，[l]表示为：
Figure BDA0003980746460000067
步骤四三、利用：
Figure BDA0003980746460000068
构造k时刻容积点、利用扰动观测器模型和故障观测器模型中测量方程传播容积点，获得k时刻的量测变量的预测值；利用k时刻的量测预测值，计算新息协方差矩阵的平方根；
其中，Xi,k|k-1表示第k时刻的最优估计值对应的容积点，Sk|k-1为k时刻状态误差协方差预测矩阵的平方根，
Figure BDA0003980746460000069
为k-1时刻对k时刻的状态量预测值；
步骤四四、利用
Figure BDA00039807464600000610
Figure BDA00039807464600000611
计算状态量与量测量的互协方差矩阵的平方根；其中，其中X1,k|k-1表示k-1时刻状态量的第1个容积点，X2n,k|k-1表示表示k-1时刻状态量的第2n个容积点，；Pxz,k|k-1表示第k时刻的状态量与量测量的互协方差矩阵,
Figure BDA00039807464600000612
其中，Z1,k|k-1表示k-1时刻状态量的第1个容积点对k时刻的量测预测值，X2n,k|k-1表示k-1时刻状态量的第2n个容积点对k时刻的量测预测值，
Figure BDA00039807464600000613
表示第k-1时刻对k时刻的量测预测值，/>
Figure BDA00039807464600000614
为第k-1时刻对k时刻的状态量预测值；
步骤四五、利用k时刻的状态变量的测量值和量测变量的预测值计算新息，利用OCSVM对新息进行数据特征判断，当新息的样本数据与OCSVM的特征数据相同时，执行步骤四六，否则，执行步骤四七；
步骤四六、计算满足卡尔曼滤波理论的状态空间方程模型中观测方程的多重渐消因子Ψ，利用状态误差协方差预测矩阵的平方根和状态量与量测量的互协方差矩阵的平方根对双重自适应强跟踪平方根容积卡尔曼滤波器进行构建；执行步骤四八；
步骤四七、计算满足卡尔曼滤波理论的状态空间方程模型中测量方程STF的多重渐消因子B，利用新息协方差矩阵的平方根和状态量与量测量的互协方差矩阵的平方根对双重自适应强跟踪平方根容积卡尔曼滤波算法进行构建；执行步骤四八；
步骤四八、计算双重自适应强跟踪平方根容积卡尔曼滤波算法的增益矩阵，对故障参数进行估计；
其中，所述双重自适应强跟踪平方根容积卡尔曼滤波算法的增益矩阵Kk：
Figure BDA0003980746460000072
故障参数估计值：
Figure BDA0003980746460000071
Figure BDA0003980746460000073
为k时刻的故障参数的最优估计值,zk表示第k时刻的状态变量xk的测量值。
本发明针对一类航天器姿态控制系统提出一种新型双重自适应强跟踪平方根容积卡尔曼滤波器(Dual Adaptive Strong Tracking SCKF,DA-ST-SCKF)用于实现系统状态和未知参数的高精度和快速的联合估计，能够提升传统卡尔曼滤波对模型不确定性的自适应能力，同时加入了对测量不确定性的自适应能力，且精度和速度均得到提升。以及基于DA-ST-SCKF和LSTM神经网络的双观测器FDIR框架，进一步地降低了模型不确定性对故障重构精度的影响。实现模型不确定和测量不确定性同时存在的情况下的高精度，快速的故障重构。本发明相较于其他成果，考虑的问题更加全面，覆盖范围广，适应性更强，且性能上得到了提升。能够提升航天器自主状态感知能力，为主动容错和健康管理打下了坚实基础。
附图说明
图1是本发明所述航天器姿态控制系统的执行器故障重构方法的流程图
图2是FDIR框架示意图；
图3是航天器姿态控制系统三轴所受总扰动信号曲线图；
图4是基于双重自适应强跟踪平方根容积卡尔曼滤波算法(DA-ST-SCKF)对φ轴总扰动的估计；
图5是基于双重自适应强跟踪平方根容积卡尔曼滤波算法(DA-ST-SCKF)对θ轴总扰动的估计；
图6是基于双重自适应强跟踪平方根容积卡尔曼滤波算法(DA-ST-SCKF)对
Figure BDA0003980746460000082
轴总扰动的估计；
图7是基于双重自适应强跟踪平方根容积卡尔曼滤波算法(DA-ST-SCKF)对三轴总扰动的估计误差；
图8是传感器出现故障故障前后的φ轴的卡尔曼滤波新息；
图9是传感器出现故障故障前后的θ轴的卡尔曼滤波新息；
图10是传感器出现故障故障前后的
Figure BDA0003980746460000081
轴的卡尔曼滤波新息；
图11是基于OCSVM对新息样本进行训练和测试的结果对比图；其中SVM Score表示数据分类分数；
图12是基于OCSVM方法对新息分类的结果图；
图13是传统阈值方法对新息分类的结果图；
图14是FDIR框架对姿态控制反作用轮1故障时的估计结果；
图15是FDIR框架对姿态控制反作用轮2故障时的估计结果；
图16是FDIR框架对姿态控制反作用轮3故障时的估计结果。
具体实施方式
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动的前提下所获得的所有其他实施例，都属于本发明保护的范围。
需要说明的是，在不冲突的情况下，本发明中的实施例及实施例中的特征可以相互组合。
具体实施方式一：下面结合图1说明本实施方式，本实施方式所述航天器姿态控制系统的执行器故障重构方法，具体过程为：
步骤一、建立航天器姿态控制系统的运动学模型和所述航天器姿态控制系统动力学模型；
步骤二、利用航天器姿态控制系统的运动学模型和所述航天器姿态控制系统动力学模型，设计滑模控制器，获取航天器姿态的控制力矩；
对所述航天器姿态控制系统动力学模型进行离散化，同时对扰动和故障参数进行状态增广和估计，获取满足卡尔曼滤波理论的状态空间方程模型；
步骤三、基于满足卡尔曼滤波理论的状态空间方程模型，采用双重自适应强跟踪平方根容积卡尔曼滤波算法(DA-ST-SCKF)和LSTM神经网络，构建FDIR框架；
步骤四、将步骤二获取航天器姿态的控制力矩和测量的航天器的角速度输入至FDIR框架，对故障参数进行估计。
进一步地，本实施方式中，步骤一中，航天器姿态控制系统的运动学模型为：
Figure BDA0003980746460000091
ω＝ωbo+Cboωoi+b+η
其中，q＝[q0,qv]T为姿态四元数，描述航天器自身固定坐标系
Figure BDA0003980746460000092
相对于惯性坐标系FI的方向，满足||q||2＝1，其中q0为q的标量部分，qv为q的矢量部分；ω为航天器自身固定坐标系/>
Figure BDA0003980746460000093
相对于惯性坐标系/>
Figure BDA0003980746460000094
的角速度，ω＝[ω1,ω2,ω3]T，其中ω1,ω2,ω3分别代表航天器在自身坐标系中三轴的角速度；Cbo表示航天器自身固定坐标系/>
Figure BDA0003980746460000095
到轨道坐标系
Figure BDA0003980746460000096
对应的旋转矩阵，/>
Figure BDA0003980746460000097
ωbo为航天器自身固定坐标系/>
Figure BDA0003980746460000098
相对于轨道坐标系/>
Figure BDA0003980746460000099
的角速度，ωoi表示轨道坐标系/>
Figure BDA00039807464600000916
相对于惯性坐标系/>
Figure BDA00039807464600000910
的角速度；b和η分别是角速率陀螺漂移和测量噪声，上角T代表矩阵的转置；I3为3×3的单位矩阵，/>
Figure BDA00039807464600000911
是q的一阶导数，/>
Figure BDA00039807464600000917
为qv的叉乘矩阵。
进一步地，本实施方式中，步骤一中，带有模型不确定性、空间环境扰动、执行机构故障、传感器故障的航天器姿态控制系统动力学模型为：
Figure BDA00039807464600000912
其中J0+ΔJ∈R3×3，J0+ΔJ表示实际转动惯量矩阵，包括标称部分J0和不确定部分ΔJ；通过进一步地整理及细化得到：
Figure BDA00039807464600000913
/>
Figure BDA00039807464600000914
Figure BDA00039807464600000915
J0+ΔJ＝Ja
Figure BDA0003980746460000105
其中，D0是航天器执行器的标称安装矩阵；ΔD代表航天器执行器实际安装与标称安装之间的误差矩阵；τ＝[τ1,τ2,…τa]T∈Ra×1是执行器的输出力矩，τi是第i个执行器的输出力矩，其中i＝1,2,…a，a表示航天器姿态控制系统执行器个数；uc为待设计的控制器输出力矩；sat(·)表示饱和函数；Λ＝diag[Λ1,…,Λa]∈Ra×a是执行器健康因子矩阵，并满足0≤Λi≤1，Λi＝1表示第i个执行器完全健康，0＜Λi＜1表示执行器部分失效，Λi＝0表示完全失效故障；Td为环境扰动力矩；d为包含不确定性和内扰动以及环境扰动力矩在内的总扰动，
Figure BDA0003980746460000104
是ω的一阶导数，/>
Figure BDA0003980746460000103
是中间变量，ω×为ω的叉乘矩阵。
本实施方式中，总干扰dk由模型不确定性和内外部干扰组成。模型不确定性主要是指系统模型与实际系统的差异，包括为了计算方便而使用简化建模，导致的未建模动态；系统参数理论值与真实值之间的误差；由于机械精度和装配工艺的限制，发射过程中的强烈冲击和振动造成的执行器和传感器的安装偏差；系统参数可能会随着系统老化或故障等而出现。外部干扰包括重力梯度力矩、辐射光压、气动力矩、剩磁力矩、非轴向力矩和内部微振动。
姿态角速率通常由陀螺仪测量，它可能会被各种误差干扰，例如常值偏差、失准、比例因子误差和噪声。具有高保真度的角速率陀螺测量模型可以建模为：
ωmk＝ωk+bk+gk+mak+vk
其中bk为k时刻陀螺漂移率偏差；vk是k时刻陀螺测量随机噪声；gk和mak分别是k时刻陀螺比例因子和陀螺的未对准误差，ωmk为陀螺的量测值即角速度实际测量值，ωk为k时刻航天器自身固定坐标系
Figure BDA0003980746460000102
相对于惯性坐标系/>
Figure BDA0003980746460000101
的角速度。
本实施方式中，角速率测量随机噪声源主要包括角白噪声、随机游走和速率随机游走，由于它们具有随机特性，它们难以被星上姿态确定算法标定。而bk,gk,mak可以通过姿态确定算法在线估计和补偿，因此此处不考虑。异常测量噪声源主要包括两个方面，一是陀螺仪内部元件由于生产过程和空间辐射而老化或失效，二是由于空间环境辐射干扰等外部问题。大量实验和理论分析表明，这些因素会导致异常漂移和测量噪声的增加，不再满足白噪声条件，可以通过地面或在轨校准来补偿漂移，因此此处主要考虑的异常噪声以增加的白噪声、有色噪声和测量异常值的形式出现。
进一步地，本实施方式中，步骤二中，获取的航天器姿态的控制力矩uc的过程为：
步骤二一、设航天器姿态的姿态角控制目标为qd，姿态角控制目标为ωd；
利用航天器姿态的姿态角控制目标为qd，计算姿态误差四元数qe＝[qe0,qev]T，qe0为误差四元数标量部分，qev为误差四元数矢量部分，计算公式为：
Figure BDA0003980746460000111
为四元数乘法运算；
利用姿态角控制目标为ωd计算角速度误差ωe＝[ωe1,ωe2,ωe3]T，ωe1,ωe2,ωe3分别为航天器在自身坐标系中三轴的角速度误差，计算公式为：ωe＝ω-ωd；
步骤二二、利用角速度误差和姿态误差四元数设计滑模控制器；利用滑模控制器获取航天器姿态的控制力矩，具体为：
Figure BDA0003980746460000112
Figure BDA0003980746460000113
Figure BDA0003980746460000114
其中k1,k2,k3是三个正标量，s为滑模面，s＝ωe+k1qev，sat(·)为饱和函数，sgn(·)为符号函数，τmax是单个执行器的最大输出力矩，τi为第i个执行机构的输出力矩，
Figure BDA0003980746460000115
是qev的一阶导数，其中/>
Figure BDA0003980746460000116
为ωe的叉乘矩阵，us是中间变量。
进一步地，本实施方式中，步骤二中，满足卡尔曼滤波理论的状态空间方程模型包括：扰动观测器模型和故障观测器模型；
所述扰动观测器模型：
Figure BDA0003980746460000117
扰动观测器模型中，前两个是状态方程，最后一个是观测方程；
所述故障观测器模型：
Figure BDA0003980746460000121
故障观测器模型中，前两个是状态方程，最后一个是观测方程；
式中，ωk+1,ωk分别代表航天器自身固定坐标系
Figure BDA0003980746460000122
相对于惯性坐标系/>
Figure BDA0003980746460000123
的角速度ω在第k+1和k时刻的值；/>
Figure BDA0003980746460000124
Figure BDA0003980746460000125
代表航天器自身固定坐标系/>
Figure BDA0003980746460000126
相对于惯性坐标系/>
Figure BDA0003980746460000127
的角速度ω的过程噪声；yk代表在第k时刻的航天器自身固定坐标系
Figure BDA0003980746460000128
相对于惯性坐标系/>
Figure BDA0003980746460000129
的角速度的测量值；vk代表航天器自身固定坐标系/>
Figure BDA00039807464600001210
相对于惯性坐标系/>
Figure BDA00039807464600001211
的角速度的测量噪声，在传统的卡尔曼滤波当中，/>
Figure BDA00039807464600001212
需要满足统计特性恒定不变的白噪声序列，τk为执行器输出控制力矩τ在k时刻的值，T离散化采样时间；Λk+1表示k+1时刻的执行器健康因子矩阵，Λk表示k时刻的执行器健康因子矩阵，/>
Figure BDA00039807464600001213
是健康因子矩阵的随机有走噪声，dk+1和dk分别是包含不确定性和内外扰动在内的总扰动在k+1时刻和k时刻的值；/>
Figure BDA00039807464600001214
为k时刻总扰动dk的随机游走噪声。/>
Figure BDA00039807464600001215
分别为k时刻需要估计并替换补偿的健康因子矩阵和总扰动。
进一步地，本实施方式中，步骤三、构建FDIR框架的过程为：
步骤三一，利用扰动观测器模型，采用双重自适应强跟踪平方根容积卡尔曼滤波算法，对航天器系统初始运行h个时刻的总扰动进行状态估计，获得扰动估计值序列
Figure BDA00039807464600001216
其中h表示系统初始运行时间结束时刻；
步骤三二、采用LSTM神经网络利用扰动估计值序列
Figure BDA00039807464600001217
预测航天器系统未来f个时刻的扰动状态，获得预测值序列/>
Figure BDA00039807464600001218
其中f＜h/2；
步骤三三、利用预测值序列
Figure BDA00039807464600001219
对故障观测器模型中当前运行m个时刻获得的总扰动[dk]m进行替换，其中，m小于或等于f；其具体替换方法为按照时间顺序一一对应进行替换；
步骤三四、故障观测器模型采用双重自适应强跟踪平方根容积卡尔曼滤波算法对系统的健康因子状态进行估计得到
Figure BDA0003980746460000131
(此时的/>
Figure BDA0003980746460000132
为故障重构结果，)直至获得p个时刻的故障估计值：/>
Figure BDA0003980746460000133
并判断/>
Figure BDA0003980746460000134
中故障估计值随时间增加是否趋于平稳，若是，将扰动观测器模型中的Λk替换为/>
Figure BDA0003980746460000135
完成一次FDIR框架的构建，否则，执行步骤三五；
步骤三五、令扰动观测器模型采用双重自适应强跟踪平方根容积卡尔曼滤波算法再次对系统运行h个时刻的扰动进行状态估计，获得扰动估计值序列
Figure BDA0003980746460000136
返回执行步骤三二，其中，r为返回次数。
本实施方式具体的FDIR方案结构如图所示。DO(扰动观测器)使用控制输入和角速度测量来估计总干扰dk。然后存储历史观测数据
Figure BDA0003980746460000137
周期性地抽取部分历史数据作为神经网络样本数据，引入长短时记忆网络(Long Short Term Memory Network,LSTM)进行训练、预测。
LSTM是循环神经网络(Recurrent Neural Network,RNN)的一个分支，具有保留之前输入数据信息的记忆功能，该网络的优点是可以更好地存储和访问历史数据信息。在LSTM模型中，包含神经元的输入以及控制输入、输出和记忆信息的三个门控单元。这三个门分别是输入门、输出门和遗忘门，这三个门控单元分别控制变量的输入、输出和记忆。
本实施方式中
Figure BDA0003980746460000138
时间序列数据集分为两部分。选取前68％的数据作为训练集，训练LSTM模型；最后32％的数据是测试集，测试网络的实际学习效果，迭代次数为80次。然后将训练数据输入LSTM模型，通过Adam优化算法调整模型参数，根据设定的迭代次数进行训练，不断降低损失函数。最后，输入测试数据，预测结果并进行评估。
得到未来一段时间内的d预测值后，将预测值作为当前时间总扰动信号，对FO中的数学模型进行补偿，以减少模型的不确定性。FO(故障观测器)利用控制输入和角速度测量计算执行器故障位置和大小，然后将重构结果代入DO(扰动观测器)的数学模型，实现执行器故障与干扰的解耦，提高故障重构的准确性。
进一步地，本实施方式中，步骤四中，对故障参数进行估计的过程为：
步骤四一、令
Figure BDA0003980746460000139
作为状态空间方程变量，对所述状态空间方程变量的初始状态估计值和状态估计误差协方差矩阵进行初始化；
步骤四二、采用Cholesky分解法，计算状态估计误差协方差矩阵的平方根Sk-1|k-1：
步骤四三、时间更新，利用公式：
Figure BDA0003980746460000141
构造k-1时刻容积点、采用扰动观测器模型和故障观测器模型中的状态方程传播容积点，利用传播后的容积点获得k-1时刻对k时刻的状态预测值；利用k-1时刻对k时刻的状态预测值，计算预测误差协方差矩阵Pk|k-1及其平方根Sxx,k|k-1；
其中，Xi,k-1|k-1表示第k-1时刻的状态方程对应的第i个容积点，Sk-1|k-1为第k-1时刻状态误差协方差预测矩阵的平方根，
Figure BDA0003980746460000142
为第k-1时刻的状态估计值；ζi表示第i个基本容积点矩阵，/>
Figure BDA0003980746460000143
n为xk的维数，此处n＝6，所述xk为/>
Figure BDA0003980746460000149
或/>
Figure BDA00039807464600001410
所述/>
Figure BDA00039807464600001411
和/>
Figure BDA00039807464600001412
维数相同，[l]i是点集[l]的第i列，[l]表示为：
Figure BDA0003980746460000144
步骤四四、测量更新，利用公式：
Figure BDA0003980746460000145
再次构造k-1时刻容积点、利用扰动观测器模型和故障观测器模型中测量方程传播容积点，用传播后的容积点获得k-1时刻对k时刻的量测预测值；利用k-1时刻对k时刻的量测预测值，计算新息协方差矩阵Pzz,k|k-1及其平方根Szz,k|k-1；
其中，Xi,k|k-1表示表示k-1时刻测量方程的第i个容积点，Sxx,k|k-1为k时刻预测误差协方差矩阵的平方根，
Figure BDA0003980746460000146
为k-1时刻对k时刻的状态量预测值；
步骤四五、利用
Figure BDA0003980746460000147
Figure BDA0003980746460000148
计算状态量与量测量的互协方差矩阵；其中，其中X1,k|k-1表示测量方程k-1时刻状态量的第1个容积点，X2n,k|k-1表示测量方程k-1时刻状态量的第2n个容积点，；Pxz,k|k-1表示第k时刻的互协方差矩阵；
Figure BDA0003980746460000151
其中，Z1,k|k-1表示k-1时刻状态量的第1个容积点对k时刻的量测预测值，Z2n,k|k-1表示k-1时刻状态量的第2n个容积点对k时刻的量测预测值，
Figure BDA0003980746460000152
表示第k-1时刻对k时刻的量测预测值，/>
Figure BDA0003980746460000153
为第k-1时刻对k时刻的状态量预测值；
步骤四六、利用k时刻的状态变量的测量值和量测变量的预测值计算新息，利用OCSVM对新息进行数据特征判断，当新息的样本数据与OCSVM的特征数据相同时，执行步骤四六，否则，执行步骤四七；
步骤四七、计算满足卡尔曼滤波理论的状态空间方程模型中观测方程的多重渐消因子Ψ，利用状态误差协方差预测矩阵的平方根和状态量与量测量的互协方差矩阵的平方根对双重自适应强跟踪平方根容积卡尔曼滤波器进行构建；执行步骤四八；
步骤四八、计算满足卡尔曼滤波理论的状态空间方程模型中测量方程STF的多重渐消因子B，利用新息协方差矩阵的平方根和状态量与量测量的互协方差矩阵的平方根对双重自适应强跟踪平方根容积卡尔曼滤波算法进行构建；执行步骤四八；
步骤四九、计算双重自适应强跟踪平方根容积卡尔曼滤波算法的增益矩阵，对故障参数进行估计；
其中，所述双重自适应强跟踪平方根容积卡尔曼滤波算法的增益矩阵Kk：
Figure BDA0003980746460000154
故障参数估计值：
Figure BDA0003980746460000155
Figure BDA0003980746460000156
为k时刻的故障参数的最优估计值,zk表示第k时刻的状态变量xk的测量值。
本实施方式中，参数的估计过程的原理具体为：
第1步：初始化
使用初始状态估计值
Figure BDA0003980746460000157
和状态估计误差协方差矩阵P0|0初始化滤波器，设定系统噪声和量测噪声所满足的初始协方差矩阵Qk,Rk，假设Qk,Rk保持恒定，不随时间k变化。以及最大估计时刻k。设定x0为xk在k＝0时刻的初始值，P0|0为误差协方差矩阵Pk|k的初始值。并根据P0|0计算S0|0，计算方式为S0|0＝chol{P0|0}，其中chol是取Cholesky分解结果的下三角矩阵的运算。
第2步：计算误差协方差矩阵的平方根Sk-1|k-1；
Sk-1|k-1表示k-1时刻的误差协方差矩阵的平方根
Sk-1|k-1＝chol{Pk-1|k-1}
其中chol是取Cholesky分解结果的下三角矩阵的运算；Pk-1|k-1为k-1时刻的误差协方差矩阵；
第3步：时间更新
生成容积点、点权重并传播容积点并估计预测状态
Figure BDA0003980746460000161
Figure BDA0003980746460000162
Figure BDA0003980746460000163
其中Xi,k-1|k-1和
Figure BDA0003980746460000164
分别表示k-1时刻状态量的第i个容积点及其对k时刻的状态预测值，其中/>
Figure BDA0003980746460000165
为第k-1时刻对k时刻的状态量预测值，f(·)具体指代状态空间方程中的非线性函数；/>
Figure BDA0003980746460000166
这里的[l]i是点集[l]的第i列，它表示以下点集：
Figure BDA0003980746460000167
然后计算预测误差协方差矩阵平方根因子Sxx,k|k-1及预测误差协方差矩阵Pk|k-1：
Figure BDA0003980746460000168
Figure BDA0003980746460000169
Figure BDA00039807464600001610
其中SQk-1是系统过程噪声Qk-1在k-1时刻的平方根，qr(·)是分解函数获得矩阵的平方根因子，
Figure BDA0003980746460000171
为容积点的权重；Pk|k-1为预测误差协方差矩阵；Sxx,k|k-1为预测误差协方差矩阵的平方根因子。
第4步：测量更新
获取容积点并通过测量函数传播容积点，然后得到预测测量值：
Figure BDA0003980746460000172
Zi,k|k-1＝h(Xi,k|k-1)
Figure BDA0003980746460000173
其中Xi,k|k-1,Zi,k|k-1分别表示k-1时刻状态量的第i个容积点及其对k时刻的量测预测值其中，
Figure BDA0003980746460000174
为第k-1时刻对k时刻的状态量预测值；h(Xi,k|k-1)为状态空间方程中的测量函数h(Xi,k|k-1)＝Xi,k|k-1+vk；/>
Figure BDA0003980746460000175
表示第k-1时刻对k时刻的量测预测值；
计算新息协方差矩阵平方根Szz,k|k-1及新息协方差矩阵Pzz,k|k-1：
Figure BDA0003980746460000176
Figure BDA0003980746460000177
Figure BDA0003980746460000178
其中，Szz,k|k-1表示第k时刻的新息协方差矩阵的平方根；
Figure BDA0003980746460000179
表示第k时刻的量测噪声矩阵Rk的平方根，即/>
Figure BDA00039807464600001710
Pzz,k|k-1为第k时刻的新息协方差矩阵；
计算互协方差矩阵的平方根
Figure BDA00039807464600001711
Figure BDA00039807464600001712
Figure BDA00039807464600001715
为加权中心矩阵；Pxz,k|k-1表示第k时刻的状态量与量测量的互协方差矩阵；
第5步：计算新息，利用OCSVM对新息进行分类：
Figure BDA00039807464600001713
rk表示第k时刻的新息；zk表示第k时刻的状态变量xk的测量值；
Figure BDA00039807464600001714
表示第k时刻的量测预测值。
假设在算法运行的一段初始时间内，传感器无故障，即量测噪声矩阵Rk未发生变化，也即rk为正常状态。基于此假设，在这段时间(h个时刻)内会积累新息rk的历史数据序列，记为[rk]h。基于[rk]h作为OCSVM的训练集S＝{([rk]h,t)|t＝1,k＝1,2…N}，其中N表示训练样本的数量，并且训练样本标签t设置为1(正常数据)，在此处，我们将由于模型不确定性导致的新息rk作为正常数据并抽取部分作为训练样本，利用训练样本对OCSVM进行训练，由测量不确定性(传感器故障)导致的新息rk作为异常数据。获取OCSVM的特征数据。
第6步：计算模型STF的多重渐消因子Ψ，对双重自适应强跟踪平方根容积卡尔曼滤波算法进行构建；
Figure BDA0003980746460000181
Figure BDA0003980746460000182
Mk＝Pzz,k|k-1+Nk-Vk
ψk＝max[1,tr(Nk)/tr(Mk)]
Figure BDA0003980746460000183
其中ρ＝0.95，σk(υ)是对角权重矩阵σk的第υ个对角元素，pυ,υ是对角矩阵Pk|k第υ行和υ列的元素；tr(·)代表对矩阵的求迹运算；ψk为模型STF单一修正因子；Pk|k-1为状态预测误差协方差矩阵，Qk为k时刻的过程噪声矩阵；σk为多重渐消因子的权重矩阵，σk(υ)为σk的第υ个对角元素；
计算多重渐消因子Ψ：
Ψ＝diag([ψk·σk(1),…,ψk·σk(n)])
引入多重渐消因子Ψ后，对状态预测误差协方差矩阵Pk|k-1修正如下：
Pk|k-1＝Ψ(Pk|k-1-Qk)+Qk
第7步：计算测量STF的多重渐消因子B，对双重自适应强跟踪平方根容积卡尔曼滤波器进行构建；
NRk＝Vk-Pzz,k|k-1+Rk
MRk＝Rk
βk＝max[1,tr(NRk)/tr(MRk)]
βk为测量STF单一修正因子；
进一步使用多重渐消因子STF将单个标量渐消因子替换为多个渐消因子，使不同的状态分量具有不同的调节能力，提高了滤波器的整体滤波性能。进一步设计权重系数如下：
Figure BDA0003980746460000191
其中，Rυ,υ为测量噪声矩阵Rk的第υ行，第υ列的元素。κk为多重渐消因子的权重矩阵，κk(υ)为κk的第υ个对角元素
将多重渐消因子校正矩阵B代入，获得修正后的量测噪声矩阵Rk
B＝diag([βk·κk(1),…,βk·κk(n)])
Rk＝BRk
第8步、计算卡尔曼增益矩阵，并利用测量预测值和卡尔曼增益矩阵计算状态量估计值，再次执行第3步和第4步，得到互协方差矩阵和新息协方差矩阵的平方根系数，即Pxz,k|k-1和Sk|k-1，并计算卡尔曼滤波增益矩阵：
Figure BDA0003980746460000192
其中Kk为k时刻的卡尔曼增益矩阵，
计算状态量最优估计值
Figure BDA0003980746460000193
其中，
Figure BDA0003980746460000194
为k时刻的状态量最优估计值；zk表示第k时刻的状态变量xk的测量值；
Figure BDA0003980746460000195
表示第k-1时刻对k时刻的量测预测值；
并更新状态误差协方差矩阵的平方根：
Figure BDA0003980746460000196
Figure BDA0003980746460000197
至此完成一个时刻的滤波过程其中，Pk|k为误差协方差矩阵；
对状态协方差矩阵进行更新，k时刻+1，重新执行第2步。
仿真实施例：
结合图2-图16来说明；通过仿真研究来验证所提出的FDIR方案的有效性和改进。初始姿态角选择为[0,0,0]Tdeg，初始角速度选择为ω0＝[0.4,-0.2,0.3]Tdeg/s。航天器在近地轨道运行，轨道角速度选择为ωoi＝0.00192rad/s。转动惯量不确定性参数为ΔJ11＝0.001,ΔJ22＝0.001,ΔJ33＝0.001。J0被选择为：
Figure BDA0003980746460000201
假设有4个反作用轮(Reaction Wheel,RW)用于航天器姿态控制。对执行器饱和进行定义，即每个执行器的最大输出力矩约束为0.5Nm。航天器标称安装矩阵RW可表示为：
Figure BDA0003980746460000202
假设，在不丧失一般性的情况下，每次机动只启用三组RW，使得D0保持精确已知。因为备份RW只有在接近完全失效或影响安全性和稳定性的情况下才会被激活。
ωt＝[0.4sin(0.05t),0.4sin(0.05t),0.4sin(0.05t)]T
航天器姿态机动利用以下控制目标，使角速度在小范围内以正弦信号的形式变化，使姿态控制系统得到充分的激励，从而能够提高对扰动和故障的识别效果。对于包含上述所有不确定性的总扰动d，设为：
d＝10-4·[d1,d2,d3]T
d1＝sin(0.0683t)sin(0.201t)+7cos(0.024t)+0.001cos(200.5t)
d2＝1.5sin(0.0683t)sin(0.035t)+0.4sin(0.038t)+0.001cos(200.5t)-0.7e0.0025t
d3＝2sin(0.0683t)+4sin(0.022t)cos(0.033t)+0.001cos(200.5t)+e0.002t
其中φ轴d1表示稳定但随机的环境扰动；θ轴d2和ψ轴d3代表RW不平衡动力学引起的扰动环境和高频微振动扰动的逐渐恶化。
利用不同的扰动场景验证了滤波算法的准确性和快速性，以及LSTM神经网络滚动预测补偿性能。图2是本发明所述FDIR框架图，图3是航天器姿态控制系统三轴所受总扰动信号曲线图，图4-7所示基于DA-ST-SCKF的总扰动估计结果及估计误差图。
为了验证DA-ST-SCKF对传感器故障的鲁棒性，在仿真过程中在550秒后加入异常测量噪声ωm(k)＝ω(k)+vf(k)。
vf(k)＝5v(k),550≤t＜730
vf(k)＝10v(k)+v(k-50)+rand,730≤t≤1000
其中ωm(k)为实测角速度，v(k)设为方差为0.0001的白噪声序列，vf为异常白噪声或彩色噪声；rand表示测量野值。具体传感器异常下的滤波新息变化如图8-10，基于不同方法对异常新息的判别结果如图11-13，基于DA-ST-SCKF和新型FDIR框架的故障重构结果如图14-16。
虽然在本文中参照了特定的实施方式来描述本发明，但是应该理解的是，这些实施例仅仅是本发明的原理和应用的示例。因此应该理解的是，可以对示例性的实施例进行许多修改，并且可以设计出其他的布置，只要不偏离所附权利要求所限定的本发明的精神和范围。应该理解的是，可以通过不同于原始权利要求所描述的方式来结合不同的从属权利要求和本文中所述的特征。还可以理解的是，结合单独实施例所描述的特征可以使用在其他所述实施例中。

Claims (9)

1.航天器姿态控制系统的执行器故障重构方法，其特征在于，具体过程为：
步骤一、建立航天器姿态控制系统的运动学模型和带有模型不确定性、空间环境扰动、执行机构故障、传感器故障的航天器姿态控制系统动力学模型；
步骤二、利用航天器姿态控制系统的运动学模型和所述航天器姿态控制系统动力学模型，设计滑模控制器，获取航天器姿态的控制力矩；
对所述航天器姿态控制系统动力学模型进行离散化，同时对扰动和故障参数进行状态增广和估计，获取满足卡尔曼滤波理论的状态空间方程模型；
步骤三、基于满足卡尔曼滤波理论的状态空间方程模型，采用双重自适应强跟踪平方根容积卡尔曼滤波算法和LSTM神经网络，构建FDIR框架；
步骤四、将步骤二获取航天器姿态的控制力矩和测量的航天器的角速度输入至FDIR框架，对故障参数进行估计。
2.根据权利要求1所述的航天器姿态控制系统的执行器故障重构方法，其特征在于，步骤一中，航天器姿态控制系统的运动学模型为：
Figure FDA0003980746450000011
ω＝ωbo+Cboωoi+b+η
其中，q＝[q0,qv]T为姿态四元数，描述航天器自身固定坐标系
Figure FDA0003980746450000012
相对于惯性坐标系FI的方向，满足||q||2＝1，其中q0为q的标量部分，qv为q的矢量部分；ω为航天器自身固定坐标系/>
Figure FDA0003980746450000013
相对于惯性坐标系/>
Figure FDA0003980746450000014
的角速度，ω＝[ω1,ω2,ω3]T，其中ω1,ω2,ω3分别代表航天器在自身坐标系中三轴的角速度；Cbo表示航天器自身固定坐标系/>
Figure FDA0003980746450000015
到轨道坐标系/>
Figure FDA0003980746450000016
对应的旋转矩阵，/>
Figure FDA0003980746450000017
ωbo为航天器自身固定坐标系/>
Figure FDA0003980746450000018
相对于轨道坐标系/>
Figure FDA0003980746450000019
的角速度，ωoi表示轨道坐标系/>
Figure FDA00039807464500000110
相对于惯性坐标系/>
Figure FDA00039807464500000111
的角速度；b和η分别是角速率陀螺漂移和测量噪声，上角T代表矩阵的转置；I3为3×3的单位矩阵，/>
Figure FDA00039807464500000112
是q的一阶导数，/>
Figure FDA00039807464500000113
为qv的叉乘矩阵。
3.根据权利要求2所述的航天器姿态控制系统的执行器故障重构方法，其特征在于，步骤一中，带有模型不确定性、空间环境扰动、执行机构故障、传感器故障的航天器姿态控制系统动力学模型为：
Figure FDA00039807464500000114
其中J0+ΔJ∈R3×3，J0+ΔJ表示实际转动惯量矩阵，包括标称部分J0和不确定部分ΔJ；
进一步地：
Figure FDA0003980746450000021
Figure FDA0003980746450000022
Figure FDA0003980746450000023
J0+ΔJ＝Ja
Figure FDA0003980746450000024
/>
其中，D0是航天器执行器的标称安装矩阵；ΔD代表航天器执行器实际安装与标称安装之间的误差矩阵；τ＝[τ1,τ2,…τa]T∈Ra×1是执行器的输出力矩，τi是第i个执行器的输出力矩，其中i＝1,2,…a，a表示航天器姿态控制系统执行器个数；uc为待设计的控制器输出力矩；sat(·)表示饱和函数；Λ＝diag[Λ1,…,Λa]∈Ra×a是执行器健康因子矩阵，并满足0≤Λi≤1，Λi＝1表示第i个执行器完全健康，0＜Λi＜1表示执行器部分失效，Λi＝0表示完全失效故障；Td为环境扰动力矩；d为包含不确定性和内扰动以及环境扰动力矩在内的总扰动，
Figure FDA0003980746450000025
是ω的一阶导数，/>
Figure FDA0003980746450000026
是中间变量，ω×为ω的叉乘矩阵。
4.根据权利要求3所述的航天器姿态控制系统的执行器故障重构方法，其特征在于，步骤二中，获取的航天器姿态的控制力矩的过程为：
步骤二一、设航天器姿态的姿态角控制目标为qd，姿态角控制目标为ωd；
利用航天器姿态的姿态角控制目标为qd，计算姿态误差四元数qe＝[qe0,qev]T，qe0为误差四元数标量部分，qev为误差四元数矢量部分，计算公式为：
Figure FDA0003980746450000027
Figure FDA0003980746450000028
为四元数乘法运算；
利用姿态角控制目标为ωd，计算角速度误差ωe＝[ωe1,ωe2,ωe3]T，ωe1,ωe2,ωe3分别为航天器在自身坐标系中三轴的角速度误差，计算公式为：ωe＝ω-ωd；
步骤二二、利用角速度误差和姿态误差四元数设计滑模控制器；利用滑模控制器获取航天器姿态的控制力矩uc；具体为：
Figure FDA0003980746450000031
Figure FDA0003980746450000032
Figure FDA0003980746450000033
其中k1,k2,k3是三个正标量，s为滑模面，s＝ωe+k1qev，sat(·)为饱和函数，sgn(·)为符号函数，τmax是单个执行器的最大输出力矩，τi为第i个执行机构的输出力矩，
Figure FDA0003980746450000034
是qev的一阶导数，其中/>
Figure FDA0003980746450000035
为ωe的叉乘矩阵，us是中间变量。
5.根据权利要求1所述的航天器姿态控制系统的执行器故障重构方法，其特征在于，步骤二中，满足卡尔曼滤波理论的状态空间方程模型包括：扰动观测器模型和故障观测器模型。
6.根据权利要求5所述的航天器姿态控制系统的执行器故障重构方法，其特征在于，所述扰动观测器模型：
Figure FDA0003980746450000036
扰动观测器模型中，前两个是状态方程，最后一个是观测方程；
式中，
Figure FDA0003980746450000037
ωk+1,ωk分别代表航天器自身固定坐标系/>
Figure FDA0003980746450000038
相对于惯性坐标系/>
Figure FDA0003980746450000039
的角速度ω在第k+1和k时刻的值；/>
Figure FDA00039807464500000310
代表航天器自身固定坐标系/>
Figure FDA00039807464500000311
相对于惯性坐标系/>
Figure FDA00039807464500000312
的角速度ω的过程噪声；yk代表在第k时刻的航天器自身固定坐标系/>
Figure FDA00039807464500000313
相对于惯性坐标系/>
Figure FDA00039807464500000314
的角速度的测量值；vk代表航天器自身固定坐标系/>
Figure FDA00039807464500000315
相对于惯性坐标系/>
Figure FDA00039807464500000316
的角速度的测量噪声，在传统的卡尔曼滤波当中，/>
Figure FDA00039807464500000317
vk需要满足统计特性恒定不变的白噪声序列，dk+1和dk分别是包含不确定性和内外扰动在内的总扰动在k+1时刻和k时刻的值；/>
Figure FDA00039807464500000318
为k时刻总扰动dk的随机游走噪声，/>
Figure FDA00039807464500000319
为k时刻需要估计并替换补偿的健康因子矩阵。
7.根据权利要求6所述的航天器姿态控制系统的执行器故障重构方法，其特征在于，所述故障观测器模型：
Figure FDA0003980746450000041
故障观测器模型中，前两个是状态方程，最后一个是观测方程；
式中，
Figure FDA0003980746450000042
τk为执行器输出控制力矩τ在k时刻的值，T离散化采样时间；Λk+1表示k+1时刻的健康因子矩阵，Λk表示k时刻的健康因子矩阵，/>
Figure FDA0003980746450000043
是健康因子矩阵的随机游走噪声，/>
Figure FDA0003980746450000044
为k时刻要估计并替换补偿的总扰动。
8.根据权利要求7所述的航天器姿态控制系统的执行器故障重构方法，其特征在于，步骤三、构建FDIR框架的过程为：
步骤三一，利用扰动观测器模型，采用双重自适应强跟踪平方根容积卡尔曼滤波算法，对航天器系统初始运行h个时刻的总扰动进行状态估计，获得扰动估计值序列
Figure FDA0003980746450000045
其中h表示系统初始运行时间结束时刻；
步骤三二、采用LSTM神经网络利用扰动估计值序列
Figure FDA0003980746450000046
预测航天器系统未来f个时刻的扰动状态，获得预测值序列/>
Figure FDA0003980746450000047
其中f＜h/2；
步骤三三、利用预测值序列
Figure FDA0003980746450000048
对故障观测器模型中当前运行m个时刻获得的总扰动[dk]m进行替换，其中，m小于或等于f；其具体替换方法为按照时间顺序一一对应进行替换；
步骤三四、故障观测器模型采用双重自适应强跟踪平方根容积卡尔曼滤波算法对系统的健康因子状态进行估计得到
Figure FDA0003980746450000049
(此时的/>
Figure FDA00039807464500000410
为故障重构结果，)直至获得p个时刻的故障估计值：/>
Figure FDA00039807464500000411
并判断/>
Figure FDA00039807464500000412
中故障估计值随时间增加是否趋于平稳，若是，将扰动观测器模型中的Λk替换为/>
Figure FDA00039807464500000413
完成一次FDIR框架的构建，否则，执行步骤三五；
步骤三五、令扰动观测器模型采用双重自适应强跟踪平方根容积卡尔曼滤波算法再次对系统运行h个时刻的扰动进行状态估计，获得扰动估计值序列
Figure FDA00039807464500000414
令/>
Figure FDA00039807464500000415
返回执行步骤三二，其中，r为返回次数。
9.根据权利要求8所述的航天器姿态控制系统的执行器故障重构方法，其特征在于，步骤四中，对故障参数进行估计的过程为：
步骤四一、令
Figure FDA0003980746450000051
作为状态空间方程变量，对所述状态空间方程变量的初始状态估计值和状态估计误差协方差矩阵进行初始化；
步骤四二、计算状态估计误差协方差矩阵的平方根Sk-1|k-1：
步骤四三、利用公式：
Figure FDA0003980746450000052
构造k-1时刻容积点、采用扰动观测器模型和故障观测器模型中的状态方程传播容积点，获得k时刻状态变量的预测值；利用k时刻的状态预测值，计算状态误差协方差预测矩阵的平方根；
其中，Xi,k-1|k-1表示第k-1时刻的最优估计值对应的容积点，Sk-1|k-1为k-1时刻状态误差协方差预测矩阵的平方根，
Figure FDA0003980746450000053
为第k-1时刻的状态预测值；ζi表示第i个容积点对应的权重，/>
Figure FDA0003980746450000054
n为xk的维数，所述xk为/>
Figure FDA0003980746450000055
或/>
Figure FDA0003980746450000056
所述/>
Figure FDA0003980746450000057
和/>
Figure FDA0003980746450000058
维数相同，[l]i是点集[l]的第i列，[l]表示为：
Figure FDA0003980746450000059
步骤四三、利用：
Figure FDA00039807464500000510
构造k时刻容积点、利用扰动观测器模型和故障观测器模型中测量方程传播容积点，获得k时刻的量测变量的预测值；利用k时刻的量测预测值，计算新息协方差矩阵的平方根；
其中，Xi,k|k-1表示第k时刻的最优估计值对应的容积点，Sk|k-1为k时刻状态误差协方差预测矩阵的平方根，
Figure FDA00039807464500000511
为k-1时刻对k时刻的状态量预测值；
步骤四四、利用
Figure FDA00039807464500000512
Figure FDA00039807464500000513
计算状态量与量测量的互协方差矩阵的平方根；其中，其中X1,k|k-1表示k-1时刻状态量的第1个容积点，X2n,k|k-1表示表示k-1时刻状态量的第2n个容积点；Pxz,k|k-1表示第k时刻的状态量与量测量的互协方差矩阵。
Figure FDA0003980746450000061
其中，Z1,k|k-1表示k-1时刻状态量的第1个容积点对k时刻的量测预测值，X2n,k|k-1表示k-1时刻状态量的第2n个容积点对k时刻的量测预测值，
Figure FDA0003980746450000062
表示第k-1时刻对k时刻的量测预测值，/>
Figure FDA0003980746450000063
为第k-1时刻对k时刻的状态量预测值；/>
步骤四五、利用k时刻的状态变量的测量值和量测变量的预测值计算新息，利用OCSVM对新息进行数据特征判断，当新息的样本数据与OCSVM的特征数据相同时，执行步骤四六，否则，执行步骤四七；
步骤四六、计算满足卡尔曼滤波理论的状态空间方程模型中观测方程的多重渐消因子Ψ，利用状态误差协方差预测矩阵的平方根和状态量与量测量的互协方差矩阵的平方根对双重自适应强跟踪平方根容积卡尔曼滤波器进行构建；执行步骤四八；
步骤四七、计算满足卡尔曼滤波理论的状态空间方程模型中测量方程STF的多重渐消因子B，利用新息协方差矩阵的平方根和状态量与量测量的互协方差矩阵的平方根对双重自适应强跟踪平方根容积卡尔曼滤波算法进行构建；执行步骤四八；
步骤四八、计算双重自适应强跟踪平方根容积卡尔曼滤波算法的增益矩阵，对故障参数进行估计；
其中，所述双重自适应强跟踪平方根容积卡尔曼滤波算法的增益矩阵Kk：
Figure FDA0003980746450000064
故障参数估计值：
Figure FDA0003980746450000065
Figure FDA0003980746450000066
为k时刻的故障参数的最优估计值,zk表示第k时刻的状态变量xk的测量值。/>
CN202211550334.7A 2022-12-05 2022-12-05 航天器姿态控制系统的执行器故障重构方法 Active CN115857529B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202211550334.7A CN115857529B (zh) 2022-12-05 2022-12-05 航天器姿态控制系统的执行器故障重构方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202211550334.7A CN115857529B (zh) 2022-12-05 2022-12-05 航天器姿态控制系统的执行器故障重构方法

Publications (2)

Publication Number Publication Date
CN115857529A true CN115857529A (zh) 2023-03-28
CN115857529B CN115857529B (zh) 2023-12-22

Family

ID=85670015

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202211550334.7A Active CN115857529B (zh) 2022-12-05 2022-12-05 航天器姿态控制系统的执行器故障重构方法

Country Status (1)

Country Link
CN (1) CN115857529B (zh)

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116923730A (zh) * 2023-07-24 2023-10-24 哈尔滨工业大学 一种具有自调节预设性能约束的航天器姿态主动容错控制方法
CN117109571A (zh) * 2023-10-25 2023-11-24 北京控制工程研究所 导航误差快速收敛方法、装置、电子设备及存储介质

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20120053780A1 (en) * 2010-08-31 2012-03-01 Seoul National University R&Db Foundation Fault Detector and Fault Detection Method for Attitude Control System of Spacecraft
CN107703742A (zh) * 2017-09-06 2018-02-16 南京邮电大学 一种柔性航天器传感器故障调节方法
CN110412997A (zh) * 2019-07-22 2019-11-05 中国人民解放军军事科学院国防科技创新研究院 一种基于神经网络的航天器姿控喷管故障智能诊断系统及方法
US20200266743A1 (en) * 2019-02-14 2020-08-20 The Board Of Trustees Of The University Of Alabama Systems, methods and devices for neural network control for ipm motor drives
CN114879515A (zh) * 2022-07-12 2022-08-09 南京航空航天大学 基于学习神经网络的航天器姿态重构容错控制方法

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20120053780A1 (en) * 2010-08-31 2012-03-01 Seoul National University R&Db Foundation Fault Detector and Fault Detection Method for Attitude Control System of Spacecraft
CN107703742A (zh) * 2017-09-06 2018-02-16 南京邮电大学 一种柔性航天器传感器故障调节方法
US20200266743A1 (en) * 2019-02-14 2020-08-20 The Board Of Trustees Of The University Of Alabama Systems, methods and devices for neural network control for ipm motor drives
CN110412997A (zh) * 2019-07-22 2019-11-05 中国人民解放军军事科学院国防科技创新研究院 一种基于神经网络的航天器姿控喷管故障智能诊断系统及方法
CN114879515A (zh) * 2022-07-12 2022-08-09 南京航空航天大学 基于学习神经网络的航天器姿态重构容错控制方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
ZE YANG: "IAR-STSCKF-Based_Fault_Diagnosis_and_Reconstruction_for_Spacecraft_Attitude_Control_Systems", 《IEEE TRANSACTIONS ON INSTRUMENTATION AND MEASUREMENT》, pages 1 - 12 *
慕晓冬: "基于注意力残差网络的航天器测控系统故障诊断", 《仪器仪表学报》, vol. 43, no. 9, pages 81 - 87 *

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116923730A (zh) * 2023-07-24 2023-10-24 哈尔滨工业大学 一种具有自调节预设性能约束的航天器姿态主动容错控制方法
CN117109571A (zh) * 2023-10-25 2023-11-24 北京控制工程研究所 导航误差快速收敛方法、装置、电子设备及存储介质
CN117109571B (zh) * 2023-10-25 2024-01-23 北京控制工程研究所 导航误差快速收敛方法、装置、电子设备及存储介质

Also Published As

Publication number Publication date
CN115857529B (zh) 2023-12-22

Similar Documents

Publication Publication Date Title
Wang et al. Active fault-tolerant control for a quadrotor helicopter against actuator faults and model uncertainties
Wang et al. Multivariate regression-based fault detection and recovery of UAV flight data
CN115857529B (zh) 航天器姿态控制系统的执行器故障重构方法
CN102037421B (zh) 辅助检测航空器的硬着陆的辅助方法
Lu et al. Adaptive three-step Kalman filter for air data sensor fault detection and diagnosis
Aboutalebi et al. A novel sensor fault detection in an unmanned quadrotor based on adaptive neural observer
Lu et al. Nonlinear aircraft sensor fault reconstruction in the presence of disturbances validated by real flight data
Sushchenko et al. Processing of redundant information in airborne electronic systems by means of neural networks
Guo et al. Fault detection and isolation for Unmanned Aerial Vehicle sensors by using extended PMI filter
Liu et al. A fault-tolerant attitude estimation method for quadrotors based on analytical redundancy
Cao et al. Anti-disturbance fault tolerant initial alignment for inertial navigation system subjected to multiple disturbances
Han et al. Quadratic-Kalman-filter-based sensor fault detection approach for unmanned aerial vehicles
CN115877717B (zh) 一种基于自抗扰控制的飞行器容错控制结构与控制方法
Wang et al. Hybrid modeling based double-granularity fault detection and diagnosis for quadrotor helicopter
Thanaraj et al. Actuator fault detection and isolation on multi-rotor UAV using extreme learning neuro-fuzzy systems
Park et al. Data-driven fault detection and isolation of system with only state measurements and control inputs using neural networks
Rudin et al. A sensor fault detection for aircraft using a single Kalman filter and hidden Markov models
Fontanella et al. Improving inertial attitude measurement performance by exploiting MEMS gyros and neural thermal calibration
Chen et al. Multi-frequency data fusion for attitude estimation based on multi-layer perception and cubature kalman filter
Liu et al. Sensor fault diagnosis of GPS/INS tightly coupled navigation system based on state chi-square test and improved simplified fuzzy ARTMAP neural network
Khan et al. Active fault tolerant flight control system design
CN113359809B (zh) 一种基于rbfnn辅助的桥梁检测无人机自主定位方法
Liu et al. LGC-Net: A lightweight gyroscope calibration network for efficient attitude estimation
CN115218927A (zh) 基于二次卡尔曼滤波的无人机imu传感器故障检测方法
Li et al. Asymptotic local approach in fault detection based on predictive filters

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