CN113311803B - 一种基于核主元分析的在轨航天器飞轮故障检测方法 - Google Patents

一种基于核主元分析的在轨航天器飞轮故障检测方法 Download PDF

Info

Publication number
CN113311803B
CN113311803B CN202110535087.2A CN202110535087A CN113311803B CN 113311803 B CN113311803 B CN 113311803B CN 202110535087 A CN202110535087 A CN 202110535087A CN 113311803 B CN113311803 B CN 113311803B
Authority
CN
China
Prior art keywords
flywheel
data
matrix
fault
spacecraft
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202110535087.2A
Other languages
English (en)
Other versions
CN113311803A (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.)
Beihang University
Original Assignee
Beihang University
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 Beihang University filed Critical Beihang University
Priority to CN202110535087.2A priority Critical patent/CN113311803B/zh
Publication of CN113311803A publication Critical patent/CN113311803A/zh
Application granted granted Critical
Publication of CN113311803B publication Critical patent/CN113311803B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G05—CONTROLLING; REGULATING
    • G05B—CONTROL OR REGULATING SYSTEMS IN GENERAL; FUNCTIONAL ELEMENTS OF SUCH SYSTEMS; MONITORING OR TESTING ARRANGEMENTS FOR SUCH SYSTEMS OR ELEMENTS
    • G05B23/00—Testing or monitoring of control systems or parts thereof
    • G05B23/02—Electric testing or monitoring
    • G05B23/0205—Electric testing or monitoring by means of a monitoring system capable of detecting and responding to faults
    • G05B23/0218—Electric testing or monitoring by means of a monitoring system capable of detecting and responding to faults characterised by the fault detection method dealing with either existing or incipient faults
    • G05B23/0221—Preprocessing measurements, e.g. data collection rate adjustment; Standardization of measurements; Time series or signal analysis, e.g. frequency analysis or wavelets; Trustworthiness of measurements; Indexes therefor; Measurements using easily measured parameters to estimate parameters difficult to measure; Virtual sensor creation; De-noising; Sensor fusion; Unconventional preprocessing inherently present in specific fault detection methods like PCA-based methods

Landscapes

  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • Automation & Control Theory (AREA)
  • Testing Of Devices, Machine Parts, Or Other Structures Thereof (AREA)

Abstract

本发明公开一种基于核主元分析的在轨航天器飞轮故障检测方法，针对飞轮不同的控制模式、检测样式，通过传感器相应地采集飞轮转速、力矩等实时状态信息，运用小波阈值法对其降噪，再进行同源扩维；通过核主元分析方法对数据进行高维映射，并根据历史正常数据得到反映偏离训练模型程度的T2、SPE控制限，进而得到综合指标故障控制限；计算实时数据的T2、SPE统计量，进而得到综合指标，根据综合指标是否超限判断有无故障。本发明对部件级的飞轮单机故障和多机故障都能达到较高的准确率，并能定位到故障飞轮；本发明在检测指标展示上更直观，利用工程实践；本发明可扩展性强，具有较大的应用前景。

Description

一种基于核主元分析的在轨航天器飞轮故障检测方法
【技术领域】
本发明针对近地轨道航天器在执行惯性稳定任务时，提供了一种基于核主元分析的在轨航天器飞轮故障检测方法，本发明属于航天器执行机构故障检测领域。
【背景技术】
自人类进入太空以来，航天器事故就屡见不鲜。据统计，1990年至2001年仅美欧日加等国发射的764个航天器中，就有121个发生故障；1990-2006年期间仅商业卫星和科学卫星故障造成的经济损失就高达44亿美元。近年来，控制分系统故障占故障总数的约30％。典型案例有1996年美国GPS BⅡ-7卫星由于姿态反作用飞轮故障，卫星完全失控；2001年意大利、荷兰联合研制的BeppoSAX卫星由于陀螺仪故障，最终于2003年坠毁等。为减少由此带来的巨大损失，航天器故障诊断技术应运而生。
美俄(苏)等主要航天强国自上世纪六七十年代就已经开始进行相关领域研究，并在诸多航天工程中采用了以状态监测为主的故障诊断方法。NASA已形成完善的航天器故障诊断体系，并将其归类为集成健康管理系统范围，其研究的感应监控系统(InductiveMonitoring System，IMS)通过聚类方法自动对数据进行分类，实验表明，复盘哥伦比亚号航天飞机失事前的数据后，IMS能较航天飞机控制中心更早发现故障，该系统还于2009年应用于国际空间站的管理。我国于2003年研制出了智能导弹武器综合测试与诊断系统，但直到2014年才在西安某卫星测控基地宣告成立首个航天器在轨故障诊断与维修实验室，且国内在线故障诊断大多针对部件级单故障，部件级多故障和系统故障则由测控站处理。
传统的基于知识的故障诊断方法过于依赖现有知识，对未知故障效果不佳，且对故障量级难以估计。而基于解析模型的故障诊断方法对模型精确性要求较高，不适用于高维和非线性等难以得到精确数学模型或是计算太过繁琐的系统，且对未知故障的识别能力不足，模型参数的实时确定也存在较大难度。基于数据驱动的故障诊断方法无须充分了解参数间的复杂关系，能够显著降低对解析模型的依赖，提高航天器在轨自主运维能力，实现相对智能的故障诊断。其所需的先验信息，几乎全部来自系统离线训练的数据，且数据模型还可根据系统状态进行灵活调整。核主元分析就是数据驱动方法的典型代表。
综上，现有故障诊断方法存在四个方面的不足，一是对未知故障的独立研判能力不足。二是对海量数据的在线分析能力不足。三是对故障数据的深度提炼能力不足。四是对复杂故障的精确辨识能力不足。
【发明内容】
本发明的目的在于针对上述不足，针对执行惯性稳定任务的卫星，提出一种基于核主元分析的在轨航天器飞轮故障检测方法，通过本发明提供的方法，可以有效利用反作用飞轮组的状态测量信息等多维数据，实现噪声抑制和降低信息的不确定性，实现航天器在线故障检测，以提高航天器自主运维能力，为后续实现故障重构和容错控制打下坚实基础。
本发明针对使用飞轮的近地轨道航天器在执行惯性稳定任务时，基于刚体动力学和核主元分析等基本理论，提出一种飞轮故障检测方法，从而为航天器在线自主运维提供有效支撑。
针对上述问题，本发明技术方案如下：
针对飞轮不同的控制模式、检测样式，通过传感器相应地采集飞轮转速、力矩等实时状态信息，运用小波阈值法对其降噪，再进行同源扩维；通过核主元分析方法对数据进行高维映射，并根据历史正常数据得到反映偏离训练模型程度的T2、SPE控制限，进而得到综合指标故障控制限；计算实时数据的T2、SPE统计量，进而得到综合指标，根据综合指标是否超限判断有无故障。
其中，T2统计量，又称霍特林(Hoteling)统计量，一般用于反映样本在主元子空间内偏离主元模型的程度。其检测原理为当测试样本得到的T2统计量值大于通过训练样本计算所得的T2控制限，则可能发生故障。
其中，SPE统计量，即平方预测误差统计量(Squared Prediction Error)，又称Q统计量，主要用于反映原始数据在残差子空间中偏离主元模型的程度。其检测原理为当测试样本得到的SPE统计量值大于通过训练样本计算所得的SPE控制限，则可能发生故障。
具体操作步骤为：
步骤1：应用本发明基于如下假设
为简化航天器控制系统设计，做出如下假设1～3。
假设1：航天器本体系视为主轴坐标系，飞轮组为航天器执行机构，其转动惯量相对航天器本体转动惯量为小量。
假设2：控制任务过程中，航天器姿态角和角速度始终处于小角度范围内。
假设3：环境干扰力矩不是本方法研究重点，做整体简化，通过轨道频率ω0为基频的时变正弦函数乘以三轴对应的干扰力矩系数进行模拟。
为得到满足故障检测方法要求的数据输入，做出如下假设4～6。
假设4：飞轮状态信息(包括输出转速、输出力矩、指令力矩、指令转速等)均可得到。
假设5：测量噪声为高斯噪声，服从正态分布。
假设6：故障诊断方法开始运行时，航天器已经处于稳定状态，在此之前，飞轮组不发生故障，飞轮的每个故障均持续一定时间。
步骤2：输入数据预处理。具体步骤如下：
步骤2.1判断飞轮控制模式
根据假设4，算法所需的关键信息可得：若飞轮为力矩控制模式，则采集每个飞轮的指令力矩和输出力矩，求得力矩差作为算法初始输入；若飞轮为转速控制模式，则采集指令转速和输出转速，求得转速差作为算法初始输入。
步骤2.2判断算法检测模式
若检测模式为单一输入，则对单个飞轮逐一进行故障检测，算法每次运行时将经降噪处理的单个飞轮差值数据直接作为下步输入；若算法检测模式为综合输入，对飞轮组整体同时进行故障检测，则将经降噪处理的单个飞轮差值数据按其对应的飞轮序号合并作为下步输入，其数据维度等于飞轮个数。
步骤2.3采集运行数据
根据假设6，航天器已经处于稳定状态。若当前时刻小于预设时刻T，算法不运行；若当前时刻大于等于T，则对此时刻至T-n+1时刻的飞轮数据进行采样，n为采样窗口大小，n＝m+s，s为前推的数据点数，m为算法单次运行时输入的样本点数，后续步骤将进行说明，若当前时刻等于T，则此时采集的数据为原始训练数据，若当前时刻大于T，则此时采集的数据为原始测试数据。
步骤2.4小波阈值法实时降噪
利用小波阈值法对初始输入的每一维度逐个进行降噪，小波基选取dbN小波或symN小波，小波分解层数最大不超过：
Figure GDA0003145226630000051
其中，
Figure GDA0003145226630000052
意为向下取整，n为样本长度。小波阈值采用通用阈值：
Figure GDA0003145226630000053
其中，
Figure GDA0003145226630000054
为噪声的标准差，n为样本长度。
步骤2.5数据同源扩维
根据飞轮控制模式和算法检测模式，将相应传感器(转速模式下为转速传感器，力矩模式下为力矩传感器；单一输入时为待检的单个飞轮的上述传感器，综合输入时为飞轮组所有飞轮的上述传感器)第t时刻采集的数据，前推s时刻，建立增广矩阵，扩展到其t-s时刻的数据，其中Xt＝[xt1 xt2 … xtr]，r为算法当前模式下须检测的飞轮个数。具体如下：
Figure GDA0003145226630000055
步骤3：实时运行数据故障检测。具体步骤如下：
步骤3.1：输入数据标准化
采取先中心化后单位化的方法进行标准化，使训练数据和测试数据的各数据列均服从期望为0、方差为1的标准正态分布。如下式：
Figure GDA0003145226630000061
其中，
Figure GDA0003145226630000062
为标准化后的训练样本矩阵第j列，
Figure GDA0003145226630000063
为原训练样本第j个变量所对应的列，
Figure GDA0003145226630000064
为
Figure GDA0003145226630000065
的列期望，
Figure GDA0003145226630000066
为
Figure GDA0003145226630000067
的列标准差。其中，所述的训练数据和测试数据两种数据均由航天器运行时产生；其中，预先采集正常工况下的一段数据作为训练数据，而测试数据则是当前待检测的实时数据。
测试样本标准化所用的列期望和列方差均来自训练样本，即：
Figure GDA0003145226630000068
其中，
Figure GDA0003145226630000069
为标准化后的测试样本矩阵第j列，
Figure GDA00031452266300000610
为原测试样本第j个变量所对应的列，
Figure GDA00031452266300000611
为
Figure GDA00031452266300000612
的列期望，
Figure GDA00031452266300000613
为
Figure GDA00031452266300000614
的列标准差。
步骤3.2：计算中心化核矩阵并特征分解
选取径向基核函数RBF作为核函数，核矩阵K第i行第j列的元素Kij的求法为：
Figure GDA00031452266300000615
求取原始训练样本核矩阵K0时，xi,yj均来自标准训练样本矩阵，且均为列向量；求取原始测试样本核矩阵K1时，xi来自标准测试样本矩阵，yj来自标准训练样本矩阵，m为数据行数。
对原始训练样本核矩阵和原始测试样本核矩阵都进行中心化：
Ktr＝K0-ImK0-K0Im+ImK0Im (7)
其中，Ktr是中心化后的训练样本核矩阵，
Figure GDA0003145226630000071
即所有元素均为
Figure GDA0003145226630000072
的m×m维方阵。
对测试样本有：
Kte＝K1-IpK0-K1Im+IpK0Im (8)
其中，Kte是中心化后的训练样本核矩阵，K0为m×m维，K1为p×m维，p为测试样本数，
Figure GDA0003145226630000073
即所有元素均为
Figure GDA0003145226630000074
的p×m维方阵。
Ktr做特征分解得到主元特征值矩阵Λ0和对应的特征向量矩阵N0，求得新的主元特征值矩阵
Figure GDA0003145226630000075
步骤3.3：主元选取及归一化特征向量
累计方差贡献率阈值设为85％。计算方法如下所示：
Figure GDA0003145226630000076
其中，k为保留的主元个数，λi为当前已保留的主元对应的特征值，λj为所有的特征值。考虑到运算量及精度要求，仅对保留的主元特征向量做标准化处理即可，由此得到负载矩阵
Figure GDA0003145226630000077
以及新的特征向量矩阵
Figure GDA0003145226630000081
步骤3.4：计算故障控制限和故障统计量
T2控制限
Figure GDA0003145226630000082
的定义式为：
Figure GDA0003145226630000083
其中，Fα(k,n-k)为显著性水平α下，服从自由度为k和n-k的F分布，k为保留的主元个数。
测试样本核矩阵的T2统计量的定义式为：
Figure GDA0003145226630000084
其中，Ti 2为第i个时刻的测试样本的T2统计量，ti＝Kte(i)P＝[t1 t2… tk]，即其在主元子空间前k个特征向量方向上的得分。Kte(i)为经中心化处理后的第i个时刻的测试样本核矩阵，P为m*k阶的负载矩阵。
SPE控制限SPEUCL的定义式为：
Figure GDA0003145226630000085
其中，
Figure GDA0003145226630000086
其中a，b分别为正常训练数据的SPE均值和方差，式中的
Figure GDA0003145226630000087
指置信度为α，自由度为h的χ2分布。
测试样本的SPE统计量的定义式为：
Figure GDA0003145226630000088
其中，SPEd为经中心化处理后的第d个时刻数据的测试样本核矩阵的SPE统计量，tdi和tdj分别为该数据在主元子空间所有特征向量方向上和前k个特征向量方向上的得分向量，ed是残差向量。
综合指标计算公式为：
Figure GDA0003145226630000091
其中，Φ是正定对称矩阵，
Figure GDA0003145226630000092
Λ-1为矩阵Λ的逆，
Figure GDA0003145226630000093
和
Figure GDA0003145226630000094
即相应的SPE和T2统计量控制限。控制限计算公式为：
Figure GDA0003145226630000095
其中，
Figure GDA0003145226630000096
S＝cov(Ktr)，即Ktr的协方差矩阵。
当
Figure GDA0003145226630000097
时，则认为有故障发生。
步骤3.5：计算故障贡献率。若对飞轮组整体同时作故障检测，则须进行此步骤，若对单个飞轮逐一作故障检测，则省略此步骤。
对多个飞轮的状态数据同时进行故障检测，通过计算其贡献率确定故障部位，一般认为，故障贡献率大的变量偏离正常模型的程度越大，但并非贡献率小的变量不存在故障。
取归一化因子为v＝[v1 v2 … vm]＝[1 1 … 1]，则第i个变量对T2统计量的贡献率
Figure GDA0003145226630000098
为：
Figure GDA0003145226630000099
Figure GDA00031452266300000910
代表对矩阵
Figure GDA00031452266300000911
求vi的偏导，其偏导矩阵对应的第p行第q列元素为：
Figure GDA00031452266300000912
第i个变量SPE统计量的贡献率
Figure GDA0003145226630000108
为：
Figure GDA0003145226630000102
带入
Figure GDA0003145226630000103
即可得
Figure GDA0003145226630000104
由于引入同源扩维环节，须对同类传感器的同源数据进行合并。计算公式为：
Figure GDA0003145226630000105
即当第i个变量不能整除传感器数(飞轮个数)r时，其故障贡献率累加到飞轮j对应的变量贡献率上，而j等于i除r取余；当第i个变量能够整除传感器数(飞轮个数)r时，其故障贡献率则累加到飞轮r对应的变量贡献率上。
Figure GDA0003145226630000106
的计算与之类似。
为便于工程实现，减轻因比对多个统计量贡献率而造成的不便，同样应将单一的T2和SPE统计量故障贡献率融合到一起，得到综合贡献率：
Figure GDA0003145226630000107
此即为综合指标的综合贡献率，从而最终确定故障飞轮。
本发明针对使用飞轮的执行惯性稳定控制任务的近地轨道航天器提出了一种基于核主元分析的在轨航天器飞轮故障检测方法，本发明主要优势体现在：1)本发明提出的故障检测方法，使得航天器在执行惯性稳定控制任务时，能够在线对飞轮运行状态进行监测，对部件级的飞轮单机故障和多机故障都能达到较高的准确率，并能定位到故障飞轮；
2)本发明提出的综合指标和综合故障贡献率计算方法融合了传统的T2和SPE两种单一统计量的特点，在检测指标展示上更直观，利用工程实践；
3)本发明支持逐一对单个飞轮和同时对飞轮组整体进行故障检测，后续可以加入其它状态测量数据以进一步提高故障检测准确度，可扩展性强，具有较大的应用前景。
【附图说明】
图1为训练集处理流程示意图。
图2为测试集处理流程示意图。
图3为故障贡献率计算流程示意图。
图4为故障检测方法整体处理流程示意图。
【具体实施方式】
下面以采用三正交飞轮组为执行机构的某型号航天器为例，具体说明本发明的实施流程。
飞轮故障建模公式为：
Tm(tf)＝emfTmt(tf)+fa (21)
tf为故障时刻，Tm(tf)为此时飞轮电机实际产生的电磁力矩，emf∈[0,1]为飞轮当前的乘性有效系数，Tmt(tf)为tf时刻飞轮电机正常情况下应产生的电磁力矩，fa为其当前的加性偏差。
选取飞轮常见的增益故障、偏差故障和空转三种故障进行仿真，将故障对应设置于飞轮1、2、3上，具体故障参数如表1所示：
表1飞轮故障参数
Figure GDA0003145226630000121
1、应用本发明基于如下假设
按照步骤1中做出假设。
基于假设1～2进行航天器控制系统建模，航天器动力学模型可简化为：
Figure GDA0003145226630000122
其中，Ib为计入飞轮质量后的整个航天器的常值惯量矩阵，ωbi为航天器在本体系下的绝对角速度，hw为整个飞轮系统的相对角动量，Tw为反作用飞轮输出的总控制力矩，Td为外干扰力矩。
四元数运动学方程为：
Figure GDA0003145226630000123
其中，q＝[q1 q2 q3]T，四元数矢量为Q＝[q0 q]T。
基于假设1，指令控制力矩可简化为：
Figure GDA0003145226630000124
其中，kp和kd分别是PD控制的比例系数和微分系数，qbid和ωbid是由星敏感器和陀螺分别测得后，经EKF算法滤波得到。
三正交飞轮组操纵律为：
Figure GDA0003145226630000131
其中，
Figure GDA0003145226630000132
为各个飞轮角加速度组成的列向量，Cw＝CIw仍为安装矩阵，C＝diag(c1 c2 c3)为飞轮组的常值安装矩阵，Iw＝diag(Iw1 Iw2 Iw3)，即反作用飞轮轴向转动惯量对角阵，
Figure GDA0003145226630000133
为Cw的逆矩阵。
飞轮电机基本数学模型为：
Figure GDA0003145226630000134
其中，Ω为飞轮转速，U为飞轮电机电枢电压，Kt为电机的转矩系数，R为电机的电枢电阻，Jw为单个飞轮转动惯量，Ke为电机的反电动势系数，
Figure GDA0003145226630000135
为电机增益系数，
Figure GDA0003145226630000136
为电机时间常数。
力矩模式下，指令输入为指令力矩，单个飞轮输出力矩为：
Figure GDA0003145226630000137
其中，Twr为飞轮输出力矩，Tm为飞轮电磁力矩，Tf为电机摩擦力矩，Twc为飞轮指令力矩，
Figure GDA0003145226630000138
K为电流反馈增益。
转速模式下，指令输入为指令转速，具体参考Bialke B提出的ITHACO Type A型飞轮高精度模型。
基于假设3，取轨道干扰力矩基频γ0＝0.0011rad/s，对环境干扰力矩整体进行简化建模如下：
Figure GDA0003145226630000139
基于假设4，假设所采用的飞轮力矩传感器测量误差是均值为0，方差为1×10-8的高斯噪声(转速模式下，假设飞轮转速传感器测量误差是均值为0，方差为1×10-6的高斯噪声)。
基于假设5，引入姿态测量误差和环境干扰力矩，由陀螺测得的航天器本体系相对惯性系的绝对角速度ωbim为：
ωbim＝ωbi+b+nb (29)
其中，b为陀螺常值漂移，取[1.2217e-4 1.2217e-4 1.2217e-4]，nb为陀螺随机噪声，取[1.7453e-5 1.7453e-5 1.7453e-5]。
由星敏感器测得的航天器本体系相对惯性系的四元数矢量Qbim为：
Figure GDA0003145226630000141
上式中，ns为星敏噪声，取[1e-41e-41e-4]。
2：输入数据预处理。具体步骤如下：
2.1判断飞轮控制模式
根据假设4，算法所需的关键信息可得；此处选择力矩模式，采集每个飞轮的指令力矩和输出力矩，求得力矩差作为算法初始输入。
2.2判断算法检测模式
此处选择单一输入，即对单个飞轮逐一进行故障检测，算法每次运行时将经降噪处理的单个飞轮指令力矩与输出力矩之差作为下步输入。
2.3采集运行数据
根据假设6，航天器已经处于稳定状态。取T＝100s，若当前时刻小于T，算法不运行；若当前时刻大于等于T，则对当前时刻至T-n+1时刻的飞轮数据进行采样，n为采样窗口大小，n＝m+s，s为前推的数据点数，m为单次运行时输入算法的训练样本点数，检测单一飞轮时取m＝100，(检测飞轮组整体时取m＝200)，若当前时刻等于T，则此数据为原始训练数据，其后续数据流向和处理步骤参照图1；若当前时刻大于T，则此数据为原始测试数据，其后续数据流向和处理步骤参照图2。
2.4小波阈值法实时降噪
利用小波阈值法对其进行预处理，本实施例中，小波基函数选取sym6小波，分解层数为1，采用通用阈值：
Figure GDA0003145226630000151
其中，
Figure GDA0003145226630000152
为噪声的标准差，n为数据长度，也即步骤2.2中的采样窗口n。
2.5数据同源扩维
检测单一飞轮时将第t时刻采集的数据，前推s时刻，取s＝2(检测飞轮组整体时取s＝0)，即扩展到t-2时刻的数据，其中Xt＝[xt1 xt2 … xtr]，r为算法当前模式下须检测的飞轮个数，本例中，单一输入时r＝1，综合输入时r＝3。
同源扩维的具体公式如下：
Figure GDA0003145226630000153
本例中对单一飞轮故障检测，则同源扩维后，t时刻算法输入的数据集为：
Figure GDA0003145226630000154
3、故障检测方法设计。具体步骤如下：
3.1：输入数据标准化
采取先中心化后单位化的方法进行标准化，使各数据列服从期望为0，方差为1的标准正态分布。如下式：
Figure GDA0003145226630000161
其中，
Figure GDA0003145226630000162
为标准训练样本矩阵第j列，
Figure GDA0003145226630000163
为原训练样本第j个变量所对应的列，
Figure GDA0003145226630000164
为
Figure GDA0003145226630000165
的列期望，
Figure GDA0003145226630000166
为
Figure GDA0003145226630000167
的列标准差。
测试样本标准化所用的列期望和列方差均来自训练样本，即：
Figure GDA0003145226630000168
其中，
Figure GDA0003145226630000169
为标准测试样本矩阵第j列，
Figure GDA00031452266300001610
为原测试样本第j个变量所对应的列，
Figure GDA00031452266300001611
为
Figure GDA00031452266300001612
的列期望，
Figure GDA00031452266300001613
为
Figure GDA00031452266300001614
的列标准差，w为数据列数，本例中为3(综合输入时由于未扩展维数，因而数据列数也为3)。
3.2：计算中心化核矩阵并特征分解
选取径向基核函数RBF作为核函数，核矩阵K第i行第j列的元素Kij的求法为：
Figure GDA00031452266300001615
求取原始训练样本核矩阵K0时，xi,yj均来自标准训练样本矩阵，且均为列向量；求取原始测试样本核矩阵K1时，xi来自标准测试样本矩阵，yj来自标准训练样本矩阵，m为数据行数，本例中为m＝n+s＝102(综合输入时m＝200)，核参数σ＝100。
对训练样本矩阵和测试样本矩阵都进行中心化：
Ktr＝K0-ImK0-K0Im+ImK0Im (37)
其中，Ktr是中心化后的训练样本核矩阵，
Figure GDA0003145226630000171
即所有元素均为
Figure GDA0003145226630000172
的m×m维方阵。
对测试样本有：
Kte＝K1-IpK0-K1Im+IpK0Im (38)
其中，Kte是中心化后的训练样本核矩阵，K0为m×m维，K1为p×m维，p为测试样本数，
Figure GDA0003145226630000173
即所有元素均为
Figure GDA0003145226630000174
的p×m维方阵。
Ktr做特征分解得到主元特征值矩阵Λ0和对应的特征向量矩阵N0，求得新的主元特征值矩阵
Figure GDA0003145226630000175
3.3：主元选取及归一化特征向量
累计方差贡献率阈值设为85％。计算方法如下所示：
Figure GDA0003145226630000176
其中，k为保留的主元个数，λi为当前已保留的主元对应的特征值，λj为所有主元的特征值。
仅对保留的特征向量做标准化处理，由此得到特征空间的负载矩阵
Figure GDA0003145226630000181
特征向量矩阵
Figure GDA0003145226630000182
3.4：计算故障控制限和故障统计量
T2控制限
Figure GDA0003145226630000183
的定义式为：
Figure GDA0003145226630000184
其中，Fα(k,n-k)为显著性水平α下，服从自由度为k和n-k的F分布，k为保留的主元个数。
测试样本核矩阵的T2统计量的定义式为：
Figure GDA0003145226630000185
其中，Ti 2为第i个时刻的测试样本的T2统计量，ti＝Kte(i)P＝[t1 t2… tk]，为该数据在主元子空间前k个特征向量方向上的得分。Kte(i)为经中心化处理后的第i个时刻的测试样本核矩阵，P为m*k阶的负载矩阵。
SPE控制限SPEUCL的定义式为：
Figure GDA0003145226630000186
其中，
Figure GDA0003145226630000187
其中a，b分别为正常训练数据的SPE均值和方差，式中的
Figure GDA0003145226630000188
指置信度为α，自由度为h的χ2分布。
测试样本的SPE统计量的定义式为：
Figure GDA0003145226630000189
其中，SPEs为经中心化处理后的第s个时刻的测试样本核矩阵的SPE统计量，tsi和tsj分别为该数据在主元子空间所有特征向量方向上和前k个特征向量方向上的得分向量，es是残差向量。
参考PCA方法得到KPCA的综合指标：
Figure GDA0003145226630000191
其中，Φ是正定对称矩阵，
Figure GDA0003145226630000192
Λ-1为主元特征向量矩阵Λ的逆，
Figure GDA0003145226630000193
和
Figure GDA0003145226630000194
即相应的SPE和T2统计量控制限。控制限计算公式为：
Figure GDA0003145226630000195
其中，
Figure GDA0003145226630000196
S＝cov(Ktr)，即Ktr的协方差矩阵。
当
Figure GDA0003145226630000197
时，则认为有故障发生。若对单个飞轮进行故障检测，不必进行步骤3.5，在此可根据综合指标是否超限判断该飞轮是否发生故障，若认为有故障则可以直接定位。
在此直接给出本专业常用的算法性能量化指标计算方法：
误报率计算公式为：
Figure GDA0003145226630000198
其中，FP为检测有故障但实际无故障的误报点数，N为实际无故障的点数，TN为检测无故障且实际无故障的正确点数。
漏报率计算公式为：
Figure GDA0003145226630000199
其中，FN为检测无故障但实际有故障的漏报点数，P为实际有故障的点数，TP为检测有故障且实际有故障的正确点数。
正确率即为所有检测结果与实际相符的点占总检测点数的比例，计算公式：
Figure GDA00031452266300001910
根据公式(46)至公式(48)得到检测性能如表2所示：
表2单一输入时综合指标性能(力矩模式)
Figure GDA0003145226630000201
3.5：计算故障贡献率
对多个飞轮的状态数据同时进行故障检测，通过计算其贡献率确定故障部位，一般认为，故障贡献率大的变量偏离正常模型的程度越大。
第i个变量对T2统计量的贡献率
Figure GDA0003145226630000202
为：
Figure GDA0003145226630000203
Knew(p·),
Figure GDA0003145226630000204
分别代表矩阵Knew的第p行和第q行，其偏导为：
Figure GDA0003145226630000205
第i个变量SPE统计量的贡献率
Figure GDA0003145226630000206
为：
Figure GDA0003145226630000207
带入
Figure GDA0003145226630000211
即可得
Figure GDA0003145226630000212
从而最终确定故障飞轮。
由于引入同源扩维环节，须对同类传感器的同源数据进行合并。计算公式为：
Figure GDA0003145226630000213
即当第i个变量不能整除传感器数(飞轮个数)r时，其故障贡献率累加到飞轮j对应的变量贡献率上，而j等于i除r取余；当第i个变量能够整除传感器数(飞轮个数)r时，其故障贡献率则累加到飞轮r对应的变量贡献率上。
Figure GDA0003145226630000214
的计算与之类似。
为便于工程实现，减轻因比对多个统计量贡献率而造成的不便，同样应将单一的T2和SPE统计量故障贡献率融合到一起，求得综合指标的故障贡献率：
Figure GDA0003145226630000215
此即为综合指标的综合贡献率，具体步骤参照图3，检测性能如表3所示：
表3综合输入时综合指标性能(力矩模式)
Figure GDA0003145226630000216
由故障参数可知，200s时仅飞轮1发生故障，300s时飞轮1、2发生故障，400s时飞轮1、2、3均发生故障，计算得到相应时刻的故障贡献率，如下表4所示：
表4综合输入时故障综合贡献率(力矩模式)
Figure GDA0003145226630000221
可知，300s时飞轮1的故障综合贡献率最大，最有可能发生故障，但这并不意味着其他飞轮未发生故障，仅代表此时飞轮1的数据偏离正常模型的程度最大，故障量对其造成的影响最为显著，飞轮2的故障综合贡献率同样不可忽视。400s时的情况与之类似。
力矩模式下，飞轮伺服系统跟踪的输入信号为指令力矩，而转速模式下，输入信号为指令转速，而姿态控制器输出的是三轴指令力矩，显然转速模式对航天器平台的控制不如力矩模式快速精确。转速模式下，若某飞轮发生灾难性故障，例如空转、卡死等，并持续较长时间，则本方法正确率将显著下降，同时误报率和漏报率将明显提高。不难理解，这是由于飞轮实际转速需要一定时间跟踪上指令转速。
转速模式下故障检测步骤与力矩模式相同，故障检测方法整体处理流程参照图4。为避免故障设置过大造成转速模式下系统直接失控，飞轮3的乘性有效系数改为emf＝0.1的增益故障，其余故障参数已在前文给出，在此直接附上结果，单一输入时如表5和表6所示。
表5单一输入时综合指标性能(转速模式)
Figure GDA0003145226630000231
表6综合输入时综合指标性能(转速模式)
Figure GDA0003145226630000232
综合输入时故障综合贡献率如表7所示：
表7综合输入时故障综合贡献率(转速模式)
Figure GDA0003145226630000233
其结果分析如力矩模式。
本方法可以较好克服传统方法对解析模型和知识的依赖，更多借助航天器在轨运行产生的实时数据进行故障检测，显著提高其自主运维能力，并能够达到较高的正确率，同时将误报率和漏报率控制在较低水平。

Claims (4)

1.一种基于核主元分析的在轨航天器飞轮故障检测方法，其特征在于：包括如下步骤：
步骤1：首先做出如下假设：
简化航天器控制系统设计，做出如下假设1～3：
假设1：航天器本体系视为主轴坐标系，飞轮组为航天器执行机构，其转动惯量相对航天器本体转动惯量为小量；
假设2：控制任务过程中，航天器姿态角和角速度始终处于小角度范围内；
假设3：环境干扰力矩不是本方法研究重点，做整体简化，通过轨道频率ω0为基频的时变正弦函数乘以三轴对应的干扰力矩系数进行模拟；
为得到满足故障检测方法要求的数据输入，做出如下假设4～6：
假设4：飞轮状态信息均可得到；
假设5：测量噪声为高斯噪声，服从正态分布；
假设6：故障诊断方法开始运行时，航天器已经处于稳定状态，在此之前，飞轮组不发生故障，飞轮的每个故障均持续一定时间；
步骤2：输入数据预处理；具体步骤如下：
步骤2.1判断飞轮控制模式
根据假设4，若飞轮为力矩控制模式，则采集每个飞轮的指令力矩和输出力矩，求得力矩差作为初始输入；若飞轮为转速控制模式，则采集指令转速和输出转速，求得转速差作为初始输入；
步骤2.2判断检测模式
若检测模式为单一输入，则对单个飞轮逐一进行故障检测，每次运行时将经降噪处理的单个飞轮差值数据直接作为下步输入；若检测模式为综合输入，对飞轮组整体同时进行故障检测，则将经降噪处理的单个飞轮差值数据按其对应的飞轮序号合并作为下步输入，其数据维度等于飞轮个数；
步骤2.3采集运行数据
根据假设6，航天器已经处于稳定状态；若当前时刻小于预设时刻T，故障检测方法不运行；若当前时刻大于等于T，则对此时刻至T-n+1时刻的飞轮数据进行采样，n为采样窗口大小，n＝m+s，s为前推的数据点数，m为算法单次运行时输入的样本点数，若当前时刻等于T，则此时采集的数据为原始训练数据，若当前时刻大于T，则此时采集的数据为原始测试数据；
步骤2.4小波阈值法实时降噪
利用小波阈值法对初始输入的每一维度逐个进行降噪，小波基选取dbN小波或symN小波，小波分解层数最大不超过：
Figure FDA0003069287290000021
其中，
Figure FDA0003069287290000022
意为向下取整，n为样本长度；小波阈值采用通用阈值：
Figure FDA0003069287290000023
其中，
Figure FDA0003069287290000024
为噪声的标准差，n为样本长度；
步骤2.5数据同源扩维
根据飞轮控制模式和检测模式，将相应传感器第t时刻采集的数据，前推s时刻，建立增广矩阵，扩展到其t-s时刻的数据，其中Xt＝[xt1 xt2 … xtr]，r为算法当前模式下须检测的飞轮个数；
Figure FDA0003069287290000031
步骤3：实时运行数据故障检测；具体步骤如下：
步骤3.1：输入数据标准化
采取先中心化后单位化的方法进行标准化，使训练数据和测试数据的各数据列均服从期望为0、方差为1的标准正态分布；如下式：
Figure FDA0003069287290000032
其中，
Figure FDA0003069287290000033
为标准化后的训练样本矩阵第j列，
Figure FDA0003069287290000034
为原训练样本第j个变量所对应的列，
Figure FDA0003069287290000035
为
Figure FDA0003069287290000036
的列期望，
Figure FDA0003069287290000037
为
Figure FDA0003069287290000038
的列标准差；
测试样本标准化所用的列期望和列方差均来自训练样本，即：
Figure FDA0003069287290000039
其中，
Figure FDA00030692872900000310
为标准化后的测试样本矩阵第j列，
Figure FDA00030692872900000311
为原测试样本第j个变量所对应的列，
Figure FDA00030692872900000312
的列期望，
Figure FDA00030692872900000313
的列标准差；
步骤3.2：计算中心化核矩阵并特征分解
选取径向基核函数RBF作为核函数，核矩阵K第i行第j列的元素Kij的求法为：
Figure FDA00030692872900000314
求取原始训练样本核矩阵K0时，xi,yj均来自标准训练样本矩阵，且均为列向量；求取原始测试样本核矩阵K1时，xi来自标准测试样本矩阵，yj来自标准训练样本矩阵，m为数据行数；
对原始训练样本核矩阵和原始测试样本核矩阵都进行中心化：
Ktr＝K0-ImK0-K0Im+ImK0Im (7)
其中，Ktr是中心化后的训练样本核矩阵，
Figure FDA0003069287290000041
即所有元素均为
Figure FDA0003069287290000042
的m×m维方阵；
对测试样本有：
Kte＝K1-IpK0-K1Im+IpK0Im (8)
其中，Kte是中心化后的训练样本核矩阵，K0为m×m维，K1为p×m维，p为测试样本数，
Figure FDA0003069287290000043
即所有元素均为
Figure FDA0003069287290000044
的p×m维方阵；
Ktr做特征分解得到主元特征值矩阵Λ0和对应的特征向量矩阵N0，求得新的主元特征值矩阵
Figure FDA0003069287290000045
步骤3.3：主元选取及归一化特征向量
累计方差贡献率阈值设为85％；计算方法如下所示：
Figure FDA0003069287290000046
其中，k为保留的主元个数，λi为当前已保留的主元对应的特征值，λj为所有的特征值；考虑到运算量及精度要求，仅对保留的主元特征向量做标准化处理即可，由此得到负载矩阵
Figure FDA0003069287290000051
以及新的特征向量矩阵
Figure FDA0003069287290000052
步骤3.4：计算故障控制限和故障统计量
T2控制限
Figure FDA0003069287290000053
的定义式为：
Figure FDA0003069287290000054
其中，Fα(k,n-k)为显著性水平α下，服从自由度为k和n-k的F分布，k为保留的主元个数；
测试样本核矩阵的T2统计量的定义式为：
Figure FDA0003069287290000055
其中，Ti 2为第i个时刻的测试样本的T2统计量，ti＝Kte(i)P＝[t1 t2 … tk]，即其在主元子空间前k个特征向量方向上的得分；Kte(i)为经中心化处理后的第i个时刻的测试样本核矩阵，P为m*k阶的负载矩阵；
SPE控制限SPEUCL的定义式为：
Figure FDA0003069287290000056
其中，
Figure FDA0003069287290000057
其中a，b分别为正常训练数据的SPE均值和方差，式中的
Figure FDA0003069287290000058
指置信度为α，自由度为h的χ2分布；
测试样本的SPE统计量的定义式为：
Figure FDA0003069287290000059
其中，SPEd为经中心化处理后的第d个时刻数据的测试样本核矩阵的SPE统计量，tdi和tdj分别为该数据在主元子空间所有特征向量方向上和前k个特征向量方向上的得分向量，ed是残差向量；
综合指标计算公式为：
Figure FDA0003069287290000061
其中，Φ是正定对称矩阵，
Figure FDA0003069287290000062
Λ-1为矩阵Λ的逆，
Figure FDA0003069287290000063
和
Figure FDA0003069287290000064
即相应的SPE和T2统计量控制限；控制限计算公式为：
Figure FDA0003069287290000065
其中，
Figure FDA0003069287290000066
S＝cov(Ktr)，即Ktr的协方差矩阵；
当
Figure FDA0003069287290000067
时，则认为有故障发生。
2.根据权利要求1所述的一种基于核主元分析的在轨航天器飞轮故障检测方法，其特征在于：若对飞轮组整体同时作故障检测，则所述的故障检测方法，进一步包括：步骤3.5：计算故障贡献率：
对多个飞轮的状态数据同时进行故障检测，通过计算其贡献率确定故障部位，一般认为，故障贡献率大的变量偏离正常模型的程度越大，但并非贡献率小的变量不存在故障；
取归一化因子为v＝[v1 v2 … vm]＝[1 1 … 1]，则第i个变量对T2统计量的贡献率
Figure FDA0003069287290000068
为：
Figure FDA0003069287290000069
Figure FDA00030692872900000610
代表对矩阵
Figure FDA00030692872900000611
求vi的偏导，其偏导矩阵对应的第p行第q列元素为：
Figure FDA0003069287290000071
第i个变量SPE统计量的贡献率
Figure FDA0003069287290000072
为：
Figure FDA0003069287290000073
带入
Figure FDA0003069287290000074
即可得
Figure FDA0003069287290000075
由于引入同源扩维环节，须对同类传感器的同源数据进行合并；计算公式为：
Figure FDA0003069287290000076
即当第i个变量不能整除传感器数r时，其故障贡献率累加到飞轮j对应的变量贡献率上，而j等于i除r取余；当第i个变量能够整除传感器数r时，其故障贡献率则累加到飞轮r对应的变量贡献率上；
Figure FDA0003069287290000077
的计算与之类似；
为便于工程实现，减轻因比对多个统计量贡献率而造成的不便，同样应将单一的T2和SPE统计量故障贡献率融合到一起，得到综合贡献率：
Figure FDA0003069287290000078
此即为综合指标的综合贡献率，从而最终确定故障飞轮。
3.根据权利要求1所述的一种基于核主元分析的在轨航天器飞轮故障检测方法，其特征在于：步骤2.5所述的相应传感器，具体是：转速模式下为转速传感器，力矩模式下为力矩传感器；单一输入时为待检的单个飞轮的上述传感器，综合输入时为飞轮组所有飞轮的上述传感器。
4.根据权利要求1所述的一种基于核主元分析的在轨航天器飞轮故障检测方法，其特征在于：步骤3.1所述的训练数据和测试数据两种数据均由航天器运行时产生；其中，预先采集正常工况下的一段数据作为训练数据，而测试数据则是当前待检测的实时数据。
CN202110535087.2A 2021-05-17 2021-05-17 一种基于核主元分析的在轨航天器飞轮故障检测方法 Active CN113311803B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110535087.2A CN113311803B (zh) 2021-05-17 2021-05-17 一种基于核主元分析的在轨航天器飞轮故障检测方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110535087.2A CN113311803B (zh) 2021-05-17 2021-05-17 一种基于核主元分析的在轨航天器飞轮故障检测方法

Publications (2)

Publication Number Publication Date
CN113311803A CN113311803A (zh) 2021-08-27
CN113311803B true CN113311803B (zh) 2022-05-17

Family

ID=77373389

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110535087.2A Active CN113311803B (zh) 2021-05-17 2021-05-17 一种基于核主元分析的在轨航天器飞轮故障检测方法

Country Status (1)

Country Link
CN (1) CN113311803B (zh)

Families Citing this family (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114235005B (zh) * 2021-11-23 2023-08-29 北京航天控制仪器研究所 一种适用于长期加电下快速更新六项陀螺仪诸元方法
CN114239321A (zh) * 2022-01-10 2022-03-25 华东理工大学 一种基于大数据的炼油过程模式识别及优化方法
CN114923503B (zh) * 2022-04-11 2024-04-19 北京航空航天大学 一种基于主元分析的在轨航天器陀螺仪和星敏感器故障诊断方法
CN114877940A (zh) * 2022-05-16 2022-08-09 中山大学 一种用于卫星飞轮的在轨故障检测方法
CN115993075B (zh) * 2022-11-14 2023-07-11 南京航空航天大学 基于sslle和自适应阈值的导弹舵面故障检测方法
CN117422673B (zh) * 2023-10-11 2024-04-05 盐城市日盛机械有限公司 一种用于飞轮制造的飞轮加工缺陷检测方法

Family Cites Families (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107153748A (zh) * 2017-06-07 2017-09-12 北京信息科技大学 基于加权核主元分析（wkpca）的回转窑故障诊断方法
CN107608335B (zh) * 2017-09-14 2020-02-14 山东科技大学 无人机飞行控制系统故障检测与故障分离的数据驱动方法
CN110347138A (zh) * 2019-05-17 2019-10-18 贵州大学 自适应核主元分析的化工过程故障诊断方法
CN111144017A (zh) * 2019-12-30 2020-05-12 北京化工大学 一种基于ff-rvm的多时段间歇过程软测量建模方法

Also Published As

Publication number Publication date
CN113311803A (zh) 2021-08-27

Similar Documents

Publication Publication Date Title
CN113311803B (zh) 一种基于核主元分析的在轨航天器飞轮故障检测方法
Jiang et al. Parameter estimation-based fault detection, isolation and recovery for nonlinear satellite models
CN102175266B (zh) 一种运动体陀螺惯性组件的故障诊断方法
CN101697079B (zh) 用于航天器实时信号处理的盲系统故障检测与隔离方法
CN114923503B (zh) 一种基于主元分析的在轨航天器陀螺仪和星敏感器故障诊断方法
Cimpoesu et al. Fault detection and diagnosis using parameter estimation with recursive least squares
Wan et al. Real-time fault-tolerant moving horizon air data estimation for the reconfigure benchmark
Han et al. Quadratic-Kalman-filter-based sensor fault detection approach for unmanned aerial vehicles
CN110531737A (zh) 基于混合模型的卫星执行机构故障诊断方法、系统及介质
Qin et al. Sensor fault diagnosis of autonomous underwater vehicle based on LSTM
Mehra et al. Autonomous failure detection, identification and fault-tolerant estimation with aerospace applications
US9921306B1 (en) Active optimal reduced state estimator
Yuksel et al. An optimal sensor fusion method for skew-redundant inertial measurement units
Sarychev et al. Optimal regressors search subjected to vector autoregression of unevenly spaced TLE series
Wang et al. Fault detection of flywheel system based on clustering and principal component analysis
Song et al. A decoupling three-position calibration method based on observability analysis for SINS/CNS integrated navigation
Patton et al. Reliable fault diagnosis scheme for a spacecraft attitude control system
Fonod et al. Robust thruster fault diagnosis: Application to the rendezvous phase of the Mars sample return mission
Wang et al. A robust state estimation method for unknown, time-varying and featureless aircraft sensor failures
Li et al. Asymptotic local approach in fault detection based on predictive filters
Fabiani et al. A NLPCA hybrid approach for AUV thrusters fault detection and isolation
Cornejo et al. Model-based fault detection for the DELFI-N3XT attitude determination system
Lopez-Encarnacion et al. Model-based FDI for agile spacecraft with multiple actuators working simultaneously
Rahimi et al. Fault isolation of reaction wheels onboard 3-axis controlled in-orbit satellite using ensemble machine learning techniques
Song et al. Fault location and separation method of Distributed Inertial Measurement Units based on IAC

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