皮纳卫星受限于体积、重量、功耗等方面的约束，无法采用传统大卫星多单机冗余备份的方式提高可靠性，一旦故障会直接影响到整个航天任务. 如何利用有限的资源，实现在轨自主故障诊断，是保证皮纳卫星长寿命、高可靠自主运行的关键.

在国外已公开的航天器在轨故障中，姿控分系统与电源分系统的故障比例最高^[1-2]. 姿控分系统包含较多敏感器和执行机构，内部集成多种算法，一旦出现故障，不仅影响卫星的姿态测量与控制，更可能使得整星失效，导致航天任务失败. 目前针对航天器姿控系统的故障检测技术主要可以划分为以下3类^[3-5]. 1）基于解析模型的方法，运用软件冗余代替硬件冗余得到残差，对残差进行处理，从而实现检测. 如Pirmoradi等^[6-7]利用卫星运动学与动力学模型，设计多组卡尔曼滤波滤波算法，通过对新息进行假设检验，实现单故障检测定位；Tudoroiu 等^[8]提出交互式多模型无极卡尔曼滤波算法，可以实现对执行机构的故障检测；贾庆贤等^[9]设计PD型学习观测器，估计卫星姿态，实现了单器件故障的检测与重构. 2）基于信号处理的检测方法，直接分析测量信号，提取频谱、幅值、方差等特征，实现故障检测，无需数学模型，适应性强. 如吴丽娜等^[10]利用离散小波变换对卫星姿控系统进行故障诊断，实现了陀螺突变故障的检测；陈婷艳等^[11]通过经验模态分解法（EMD）提取故障特征，实现航天器运动器件的故障检测；Gueddi等^[12-13]利用顶点主成分分析（VPCA），实现了执行器的故障检测与隔离. 3）基于人工智能的检测方法，无需精确模型，现阶段研究最多，也是今后自主检测发展的方向. 如樊久铭等^[14]基于模糊集理论，设计分布式模糊专家系统，初步验证了其在电源故障检测中的可行性；王展等^[15-17]分别利用神经网络，实现了器件故障类型的分类以及飞轮的故障检测. 综上所述，现有研究针对姿控系统的检测难以适应多器件故障的情况，偏向于单器件单故障类型的检测，适应的故障类型较少；人工智能的方法和部分阈值判断方法多数需要大量的历史数据以及专家经验，计算量巨大，不适用于皮纳卫星的在线实时诊断，难以满足皮纳卫星自主可靠运行的需求.

本文针对皮纳卫星姿控系统的自主故障检测问题，提出自上而下的分层检测方法. 将姿控系统划分为系统层和器件层，系统层进行全局监测，器件层完成故障定位，通过分层节省了不必要的计算资源，能够实现多器件多故障的检测.

皮纳卫星姿态确定与控制系统（attitude determination and control system，ADCS）主要由姿态敏感器、姿态控制器及中央处理器组成，完成姿态确定与控制. 图1给出姿态控制系统组成框图.

姿态确定系统（attitude determination system，ADS）是卫星实现特定指向、维持在轨稳定的重要前提. 常用姿态测量敏感器有磁强计、模拟太阳敏感器、数字太阳敏感器、星敏感器、红外地球敏感器及陀螺，ADS系统利用姿态敏感器获取磁场矢量、太阳矢量、恒星矢量、地球矢量及惯性角速度等测量信息，利用姿态确定算法实现对卫星在轨实时姿态的评估. 常用的姿态确定算法有双矢量算法和卡尔曼滤波. 姿态控制系统（attitude control system, ACS）是保证卫星姿态稳定与卫星功能实现的关键，常用的姿态控制器有反作用飞轮组、偏置轮、磁力矩器、推力器等. 控制系统根据姿态信息的估计值和参考值之间的偏差，实时解算所需的控制力矩，完成对系统的闭环控制.

ADCS故障主要可以分为敏感器故障、执行器故障以及算法异常. 根据故障时变特性，敏感器故障可以分为间歇性突变故障、偏差故障、缓变故障、输出卡死故障^[18]. 间歇性突变是由于环境突变或器件本身接触不良的原因导致器件出现瞬时突变，对系统的影响较小；偏差故障一般是由于元器件失效，导致信号出现一定偏差，通常表现为阶跃变化；缓变故障是由于器件老化导致性能退化或因噪声引起信号漂移等缓慢变化的故障；输出卡死故障则表现为输出信号某一时刻起保持常值. 敏感器故障会直接导致输出信号异常，使得接入该敏感器的定姿算法出错，导致控制力矩计算错误，在不进行故障隔离的情况下，故障影响将逐级扩大，极有可能使得卫星失控.

执行机构故障主要考虑反作用轮，最常见的为反作用轮力矩输出异常，这是由于轴承温度增加、润滑不好的原因，使得摩擦力矩增大、反作用轮增益效率下降导致；空转故障是指反作用轮由于某些因素不能响应正常的控制力矩指令，转速几乎维持不变（在摩擦力作用下缓慢减少），输出力矩几乎为零；停转故障是指反作用轮输出力矩在产生一个巨大反向扰动后转速迅速变为零的情况. 执行机构故障时，不影响敏感器正常工作，控制效率下降.

算法异常指参数传递错误、无数据、异常中断等软件方面导致的问题，与敏感器和执行机构无关，暂不考虑该类异常情况.

提出的故障检测方案采用分层方式，将姿控系统划分为系统层和器件层. 系统层利用非线性观测器残差，评估卫星是否存在异常，该层检测到故障，则进入器件层算法，否则结束该轮检测. 器件层利用定姿算法中卡尔曼新息评估使用的敏感器工作状态，通过对采样信号进行小波变换，结合假设检验实现敏感器的故障检测与定位. 在陀螺可用的前提下，利用数字陀螺估计误差实现执行机构故障检测，完成所有器件的故障检测.

因可能存在器件故障在系统层反映不明显的情况，当系统层残差判断连续N次正常时，进行一次器件层检测，以降低漏检率. 对于定姿新息异常但小波分析正常的情况，汇总检测结果后，切换至另一组敏感器进行定姿，结合2组检测结果可以定位故障器件. 分层检测方案如图2所示.

系统层的检测是构建一个全局观测器，通过对残差进行阈值判断，确定是否有器件故障. 在无故障时，通过该层的检测，可以略过器件层的运算，节省计算资源.

卫星运动学与动力学方程如下：

(1) ${{{\dot q}}_{{\rm{bi}}}} = {{\varOmega }}\left( {{{{w}}_{{\rm{bi}}}}} \right){{{q}}_{{\rm{bi}}}}/2,$

(2) ${{\dot w}} = {{{J}}^{ - 1}}\left\{ {{{Jw}} \times {{w}} + {{{C}}_{\rm{m}}}{{h}} \times {{w}} + {{{T}}_{\rm{c}}} - {{{C}}_{\rm{m}}}{{\dot h}} + {{{T}}_{\rm{d}}}} \right\}.$

式中： ${{J}}$为卫星转动惯量， ${{{C}}_{\rm{m}}}$为飞轮安装矩阵， ${{{T}}_{\rm{c}}}$为控制力矩， ${{{T}}_{\rm{d}}}$为干扰力矩， ${{h}}$为飞轮产生的角动量，并有

${{\varOmega }}\left( {{{{w}}_{{\rm{bi}}}}} \right) = \left[ {\begin{array}{*{20}{c}} { - [{{w}}_{{\rm{bi}}}^ \times ]}&{{{{w}}_{{\rm{bi}}}}}\\ { - {{w}}_{{\rm{bi}}}^{\rm{T}}}&0 \end{array}} \right].$

取状态变量 ${{x}} = {[\begin{array}{*{20}{c}}{{{{w}}^{\rm{T}}}}&{{{{q}}^{\rm{T}}}}\end{array}]^{\rm{T}}}$，以陀螺输出和定姿所得姿态四元数作为量测，则敏感器和执行器信息均可以接入观测器，能够实现姿控系统的全局故障检测.

考虑如下非线性系统：

(3) $\left.\begin{array}{l} {{\dot x}} = {{Ax}} + g\left( {{x}} \right) + {{Bu}} + {{B}}{{{f}}_{\rm{a}}} + {{{\delta }}_{\rm{s}}}{\rm{ + }}{{{T}}_{\rm{d}}},\\ {{y}} = {{Cx}}{\rm{ + }}{{{f}}_{\rm{s}}}{\rm{ + }}{{{v}}_{\rm{s}}}. \end{array}\right\}$

式中： $g\left( {{x}} \right)$为非线性部分； ${{u}}$为控制力矩； ${{{f}}_{\rm{a}}}$为执行器故障力矩； ${{{f}}_{\rm{s}}}$为敏感器等效故障矢量； ${{{\delta}} _{\rm{s}}}$为建模误差； ${{{T}}_{\rm{d}}}$为干扰力矩； ${{{v}}_{\rm{s}}}$为敏感器噪声； ${{y}}$为敏感器测量输出，在正常状态下，敏感器输出信息可以看作由真值与零均值高斯噪声构成.

假设系统中（A，C）可观测，且满足如下条件.

a）g（x）为Lipschitz函数，有

(4) $\begin{array}{l} \left\| {\left. {g\left( {{{{x}}_1}} \right) - g\left( {{{{x}}_2}} \right)} \right\|} \right. \leqslant r\left\| {\left. {{{{x}}_1} - {{{x}}_2}} \right\|} \right.; \;\;\;r > 0. \end{array}$

b）干扰力矩有界，即存在 $\eta > 0$满足：

(5) $\left\| {\left. {{{{T}}_{\rm{d}}}} \right\|} \right. \leqslant \eta .$

设计观测器为

(6) $\left.\begin{array}{l} {{\dot {\hat x}}} = {{A\hat x}} + g\left( {{{\hat x}}} \right) + {{Bu}} + {{L}}\left( {{{y}} - {{\hat y}}} \right),\\ {{\hat y}} = {{C\hat x}}. \end{array}\right\}$

式中： $\hat {{x}}$表示观测器估计状态； ${{L}}$为增益矩阵，并能保证（A−LC）是稳定矩阵.

定义状态误差为 ${{e}} = {{x}} - {{\hat x}}$，输出误差 ${{\tilde y}} = {{y}} - {{\hat y}}$，则由式（1）、（2）可得状态误差方程：

(7) $\left.\begin{array}{l} {{\dot e}} = ({{A}} - {{LC}}){{e}} + g\left( {{x}} \right) - g\left( {{{\hat x}}} \right) + {{B}}{{{f}}_{\rm{a}}} + \\ \quad\;\;{{{\delta }}_{\rm{s}}} + {{{T}}_{\rm{d}}}{\rm{ + }}{{L}}{{{f}}_{\rm{s}}}{\rm{ + }}{{L}}{{{v}}_{\rm{s}}},\\ {{\tilde y}} = {{Ce}} + {{{f}}_{\rm{s}}}{\rm{ + }}{{{v}}_{\rm{s}}}. \end{array}\right\}$

考虑建模无误差，且未出现器件故障的情况，式（7）可以简化为

(8) $\left.\begin{array}{l} {{\dot e}} = ({{A}} - {{LC}}){{e}} + g\left( {{x}} \right) - g\left( {{{\hat x}}} \right) + {{{T}}_{\rm{d}}}{\rm{ + }}{{L}}{{{v}}_{\rm{s}}},\\ {{\tilde y}} = {{Ce}}{\rm{ + }}{{{v}}_{\rm{s}}}. \end{array}\right\}$

引理1^[19] 　给定常数μ>0和对称正定矩阵Q，可以推出

(9) $2{{{x}}^{\rm{T}}}{{y}} \leqslant \frac{1}{\mu }{{{x}}^{\rm{T}}}{{Qx}} + \mu {{{y}}^{\rm{T}}}{{{Q}}^{ - 1}}{{y}};\;\;\;\;{{x}},{{y}} \in {{\bf R}^n}.$

引理2^[20] 　给定对称正定矩阵P，可以推出

(10) $2{{{e}}^{\rm{T}}}{{P}}\left[ {g\left( {{x}} \right) - g\left( {{{\hat x}}} \right)} \right] \leqslant {r^2}{{{e}}^{\rm{T}}}{{PPe}} + {{{e}}^{\rm{T}}}{{e}}.$

定义如下变量：

(11) ${{E}} = 2{\left( {{{A}} - {{LC}}} \right)^{\rm{T}}}{{P}} + {r^2}{{PP}} + {{I}} + \mu {{P}}{{{Q}}^{ - 1}}{{P}}.$

若满足E为负定矩阵，则观测器状态估计误差可以渐进收敛.

证明：

选取李雅普诺夫函数为

(12) $V = {{{e}}^{\rm{T}}}{{Pe}}.$

结合式（6）~（8），微分式可以表示为

(13) $\begin{split} {\dot V}& = {{{{\dot e}}}^{\rm{T}}}{{Pe}} + {{{e}}^{\rm{T}}}{{P\dot e}}=\\ &{\left[ {({{A}} - {{LC}}){{e}} + g\left( {{x}} \right) - g\left( {{{\hat x}}} \right) + {{{T}}_{\rm{d}}}{\rm{ + }}{{L}}{{{v}}_{\rm{s}}}} \right]^{\rm{T}}}{{Pe}}+\\ & {{{e}}^{\rm{T}}}{{P}}\left[ {({{A}} - {{LC}}){{e}} + g\left( {{x}} \right) - g\left( {{{\hat x}}} \right) + {{{T}}_{\rm{d}}}{\rm{ + }}{{L}}{{{v}}_{\rm{s}}}} \right]=\\ & {{{e}}^{\rm{T}}}\left[ {{{\left( {{{A}} - {{LC}}} \right)}^{\rm{T}}}{{P}} + {{P}}\left( {{{A}} - {{LC}}} \right)} \right]{{e}}+\\ & 2{\left[ {g\left( {{x}} \right) - g\left( {{{\hat x}}} \right)} \right]^{\rm{T}}}{{Pe}} + 2{\left( {{{{T}}_{\rm{d}}} + {{L}}{{{\nu }}_{\rm{s}}}} \right)^{\rm{T}}}{{Pe}}\leqslant\\ & {{{e}}^{\rm{T}}}\left[ {{{\left( {{{A}} - {{LC}}} \right)}^{\rm{T}}}{{P}} + {{P}}\left( {{{A}} - {{LC}}} \right)} \right]{{e}} + {r^2}{{{e}}^{\rm{T}}}{{PPe}}+\\ & {{{e}}^{\rm{T}}}{{e}} + \frac{1}{\mu }{\lambda _{\max }}\left( {{Q}} \right){\left\| {\left. {{{{T}}_{\rm{d}}} + {{L}}{{{\nu }}_{\rm{s}}}} \right\|} \right.^2} + \mu {{{e}}^{\rm{T}}}{{P}}{{{Q}}^{ - 1}}{{Pe}}\leqslant\\ & {{{e}}^{\rm{T}}}\left[ {{{\left( {{{A}} - {{LC}}} \right)}^{\rm{T}}}{{P}} + {{P}}\left( {{{A}} - {{LC}}} \right)} \right]{{e}} + {r^2}{{{e}}^{\rm{T}}}{{PPe}}+\\ & {{{e}}^{\rm{T}}}{{e}} + \frac{1}{\mu }{\lambda _{\max }}\left( {{Q}} \right){\left\| {\left. {{{{T}}_{\rm{d}}} + {{L}}{{{\nu }}_{\rm{s}}}} \right\|} \right.^2} + \mu {{{e}}^{\rm{T}}}{{P}}{{{Q}}^{ - 1}}{{Pe}}. \end{split}$

取 $\varsigma {\rm{ = }}{\lambda _{\max }}\left( {{Q}} \right){\left\| {\left. {{{{T}}_{\rm{d}}} + {{L}}{{{\nu }}_{\rm{s}}}} \right\|} \right.^2}/ \mu \geqslant 0$，将式（11）代入式（13），可以化简得到 $\dot V \leqslant - {\lambda _{\min }}\left( { - {{E}}} \right){\left\| {\left. {{e}} \right\|} \right.^2} + \varsigma $，则只需E负定，即可以使得状态误差一致有界，证毕.

假设系统无故障时，状态误差的阈值为τ，系统一旦出现故障，观测器状态误差发散，超出阈值，即系统层的检测条件可以设置为

(14) $ {\left\| {{{\tilde y}}} \right\|} \left\{ {\begin{array}{*{20}{c}} \leqslant \tau,\;\,{\text{系统正常}};\\ > \tau,\;\,{\text{系统故障}}. \end{array}} \right.$

式中： ${{\tilde y}}$为量测值与观测器估计值的差值.

器件层检测主要实现故障定位. 卫星卡尔曼定姿算法可以参见文献[7]，本文不再赘述. 根据文献[21]可知，卡尔曼滤波在收敛后新息有界，一旦器件出现故障，会使得算法中的参数不匹配，新息超出阈值，对相关敏感器进行小波分析，可以定位故障源.

在卫星稳定后，每个控制周期所需的控制力矩较小，使得飞轮转速变化小，其发生增益下降故障时对转速的影响不大，小波系数变化不明显，无法检测出故障；在飞轮发生空转故障后，转速几乎保持不变，小波系数将趋于零，难以通过阈值判断检测出该类故障. 基于扩展卡尔曼滤波算法 (EKF) 与动力学模型设计的数字动力学陀螺，将飞轮故障转化为输出力矩异常问题，由新息变化反映执行机构故障，从而实现敏感器与执行机构故障隔离与定位.

小波多分辨分析的基本思想是将主体空间用相互正交的尺度空间与小波空间表示. 尺度空间是对信号作概貌近似，表示信号的发展趋势，反映信号的低频部分；小波空间对信号作细节近似，反映信号的高频部分，涵盖信号的故障信息^[22-23].

考虑到短支集小波在处理间断问题时更有效，选择支集长度为6的dB3作为小波基. 为了简化计算，采用离散二进小波变换：

${W_{2{{m}}}} = \left\langle {\left. {f\left( t \right),{\psi _{2{{m}}}}\left( k \right)} \right\rangle } \right. = {2^{ - {{m}}}}\int\limits_{ - \infty }^{ + \infty } {f\overline {\psi \left( {{2^{{{ - m}}}}t - n} \right)} } {\rm d}t.$

为了便于卫星在线检测，采用适当宽度的时间窗截取敏感器采样信号，进行2层小波分解. 对每层高频小波系数进行小波重构，通过重构信号中的极值确定器件状态，重构算法如下：

$f\left( t \right) = \sum\limits_{m \in Z} {{W_{2{{m}}}}f\left( k \right)} * {\psi _{2{{m}}}}\left( t \right).$

令 ${{{\chi}} _{{j}}}$为第j层高频小波系数重构信号，在正常采样信号中，各时间窗内 ${{{\chi}} _{{j}}}$极值存在较大差异，选用全局阈值进行判断极易造成误判，根据时间窗采用不同的阈值，极大地增加了工作量. 根据文献[24]可知，当原始信号可以看作由真值和零均值高斯噪声构成时，高频重构信号 ${{{\chi}} _{{j}}}$满足高斯分布. 由式（3）可知，敏感器输出信号满足上述条件，由此可以提出如下假设检验：

$\left. \begin{array}{l} {H_0}:{{{\chi }}_{{j}}}\text{服从高斯分布} \Rightarrow \text{正常};\\ {H_1}:{{{\chi }}_{{j}}}\text{不服从高斯分布} \Rightarrow \text{故障}. \end{array} \right\}$

重构信号数据点数较少，采用Anderson-Darling（AD）检验重构信号的正态性. AD检验将样本数据的经验累积分布函数与假设数据呈高斯分布时理论的分布进行比较，得到统计量，与同类型的KS检验相比，AD检验在小样本条件下，具有良好的检验性能. AD检验统计量可以按式（14）计算：

(15) ${A^2} = - n - \sum\limits_{i = 1}^n {\frac{{2i - 1}}{n}\left[ {\ln \left( {\varPhi \left( {w_{{j}}^{{i}}} \right)} \right) - \ln \left( {\varPhi \left( {w_{{j}}^{{{n}} + 1 - {{i}}}} \right)} \right)} \right]} .$

式中：n为数据量个数，Φ为序列假设的高斯分布函数.

根据所选取的显著性水平α求得临界值，若A²小于临界值则接受零假设，即重构信号符合高斯分布，器件正常，否则在显著性水平下α拒绝零假设.

卫星动力学方程包含角速度信息，通过结合部分矢量敏感器信息，可以实现卫星的无陀螺定姿^[25-26]. 本文基于动力学方程设计EKF算法，以陀螺采样值作为量测，得到角速度的最优估计，整个滤波过程即为数字动力学陀螺. 由文献[21]可知，在无故障且算法收敛时，新息有界. 在执行机构故障后，状态递推方程与状态均方误差阵不再匹配，将导致算法发散，新息超出阈值.

卫星所受的外部干扰力矩主要包括重力梯度力矩、气动力矩、辐射力矩及磁力矩，可以将外部干扰合力矩简化为

(16) ${{{T}}_{\rm{d}}} = {{{T}}_{{\rm{dcon}}}} + {{\upsilon }}.$

式中： ${{{T}}_{{\rm{dcon}}}}$为常值干扰力矩； ${{\upsilon }} $为未知干扰力矩，以零均值高斯白噪声表示.

取状态变量Δx=Δw，则将式（2）线性化可得

(17) $\Delta {{\dot x}} = {{F}}\Delta {{x}} + {{B}}{{{T}}_{\rm{c}}} + {{G\upsilon }}.$

式中：Δx为角速度估计误差， ${{B}} = {{G}} = {{{J}}^{ - 1}}$， ${{F}} = - {{{J}}^{ - 1}}{{C}} + {{{J}}^{ - 1}}{\left( {{{{C}}_{\rm{m}}}{{h}}} \right)^ \times }$，

${{C}} = \left[ {\begin{array}{*{20}{c}} {{J_{{{xz}}}}{w_{{y}}} - {J_{{{xy}}}}{w_{{z}}}}&{{J_{{{xz}}}}{w_{{x}}} - {J_{{y}}}{w_{{z}}} + 2{J_{{{yz}}}}{w_{{y}}} + {J_{{z}}}{w_{{z}}}}&{{J_{{z}}}{w_{{y}}} - {J_{{y}}}{w_{{y}}} - 2{J_{{{yz}}}}{w_{{z}}} - {J_{{{xy}}}}{w_{{x}}}}\\ {{J_{{x}}}{w_{{z}}} - 2{J_{{{xz}}}}{w_{{x}}} - {J_{{{yz}}}}{w_{{y}}} - {J_{{z}}}{w_{{z}}}}&{{J_{{{xy}}}}{w_{{z}}} - {J_{{{yz}}}}{w_{{x}}}}&{{J_{{x}}}{w_{{x}}} + {J_{{{xy}}}}{w_{{y}}} + 2{J_{{{xz}}}}{w_{{z}}} - {J_{{z}}}{w_{{x}}}}\\ {2{J_{{{xy}}}}{w_{{x}}} - {J_{{x}}}{w_{{y}}} + {J_{{y}}}{w_{{y}}} + {J_{{{yz}}}}{w_{{z}}}}&{{J_{{y}}}{w_{{x}}} - 2{J_{{{xy}}}}{w_{{y}}} - {J_{{{xz}}}}{w_{{z}}} - {J_{{x}}}{w_{{x}}}}&{{J_{{{yz}}}}{w_{{x}}} - {J_{{{xz}}}}{w_{{y}}}} \end{array}} \right].$

考虑陀螺仪模型^[27]为

(18) ${{{w}}_{\rm{m}}}\left( t \right) = {{w}}\left( t \right) + {{b}}\left( t \right) + {{{n}}_{\rm{A}}}\left( t \right).$

(19) ${{\dot b}}\left( t \right) = {{{n}}_{\rm{R}}}\left( t \right).$

式中： ${{{w}}_{\rm{m}}}$为陀螺输出值， ${{w}}$为卫星实际角速度在陀螺输入轴的投影， ${{b}}\left( t \right)$为陀螺等效常值漂移， ${{{n}}_{\rm{A}}}$、 ${{{n}}_{\rm{R}}}$为不相关的零均值白噪声.

根据卡尔曼滤波方程进行角速度估计. 式（20）称为卡尔曼预估方程，式（21）称为卡尔曼校正方程.

(20) $\begin{array}{l} \Delta {{{x}}_{{{k}} + 1|{{k}}}} = {{{\varPhi }}_{{k}}}\Delta {{{x}}_{{k}}}.\\ {{{P}}_{{{k}} + 1|{{k}}}} = {{{\varPhi }}_{{k}}}{{{P}}_{{k}}}{{\varPhi }}_{{k}}^{\rm{T}} + {{{G}}_{{k}}}{{{Q}}_{{k}}}{{G}}_{{k}}^{\rm{T}}. \end{array} $

式中： ${{{\varPhi }}_{{k}}} = {{{I}}_{6 \times 6}} + {{F}}T$，其中T为递推周期.

(21) $\left.\begin{array}{l} {{{K}}_{{{k}} + 1}} = {{{P}}_{{{k}} + 1|{{k}}}}{{H}}_{{{k}} + 1}^{\rm{T}}{\left( {{{{H}}_{{{k}} + 1}}{{{P}}_{{{k}} + 1|{{k}}}}{{H}}_{{{k}} + 1}^{\rm{T}} + {{{R}}_{{{k}} + 1}}} \right)^{ - 1}},\\ \Delta {{{x}}_{{{k}} + 1}} = {{{\varPhi }}_{{k}}}\Delta {{{x}}_{{k}}} + {{{K}}_{{{k}} + 1}}\Delta {{{y}}_{{{k}} + 1}},\\ {{{P}}_{{{k}} + 1}} = \left( {{{{I}}_{6 \times 6}} - {{{K}}_{{{k}} + 1}}{{{H}}_{{{k}} + 1}}} \right){{{P}}_{{{k}} + 1|{{k}}}}. \end{array}\right\}$

计算执行机构正常时EKF算法收敛后新息 $\Delta {{{y}}_{{{k}} + 1}}$模的最大值 $\sigma $，则可以设置如下检测条件：

(22) $\left. {\begin{array}{*{20}{c}} {\left\| {\left. {\Delta {{{y}}_{{{k}} + 1}}} \right\|} \right. \leqslant \sigma ,\;\;\text{执行器正常}},\\ {\left\| {\left. {\Delta {{{y}}_{{{k}} + 1}}} \right\|} \right. > \sigma ,\;\;\text{执行器故障}}. \end{array}} \right\}$

为了验证分层算法的准确性与实时性，基于浙江大学皮星二号相关参数的设计仿真场景，实现对算法多器件故障检测效果的评估.

卫星采用太阳同步轨道，初始姿态角和角速度分别设置为0.1°与0.1°/s；敏感器模型参数如表1所示.

卫星转动惯量取为

${{J}} = \left[ {\begin{array}{*{20}{c}} {0.578\;6}&{0.000\;3}&{0.007\;3} \\ {0.000\;3}&{0.581\;0}&{0.000\;3} \\ {0.007\;3}&{0.000\;3}&{0.238\;3} \end{array}} \right].$

计算得到利普希茨参数r=1.001，非线性观测器的相关参数分别如下所示：

${{P}} = 1.143\;6 {{{I}}_{7 \times 7}},$

${{L}} = 1.093\;0 {{{I}}_{7 \times 7}}.$

式中：I_7×7为7阶单位矩阵.

取显著性水平α=0.05，查表可得对应的AD检验临界值A₀=0.732 9.

考虑到卫星在轨运行中，陀螺仪可能因为硬件接触不良，发生间歇性突变故障；星敏感器在长时间使用后会由于器件老化，使得输出信号存在固定偏差；模拟太阳敏感器的太阳能电池片可能出现部分失效的情况；磁强计的磁敏电阻老化，会使得输出信号呈比例增大；供电模块的异常有可能使得飞轮无法通电，导致停转. 设置2种故障场景，如表2所示.

仿真中各故障设计如下.

${f_{{\rm{gyro}}\_x}} = \left\{ \begin{array}{l} 0,\;\,t < 51 ; \\ {0.02^ \circ },\;\,t = 51; \\ 0,\;\,t > 51 ; \\ \end{array} \right.\begin{array}{*{20}{c}} {}&{{f_{{\rm{gyro}}\_y}} = 0};&{{f_{{\rm{gyro}}\_z}} = 0} \end{array}.$

${f_{{\rm{star}}\_x}} = \left\{ \begin{array}{l} 0,\;\,t < 51; \\ {1^ \circ },\;\,t \geqslant 51; \\ \end{array} \right.\begin{array}{*{20}{c}} {}&{{f_{{\rm{star}}\_y}} = 0};&{{f_{{\rm{star}}\_z}} = 0} \end{array}.$

${{{s}}_{{\rm s}{\rm{un}}}} = \left\{ \begin{array}{l} {{{s}}_{{\rm s}{\rm{un}}}},\;\;t < 112; \\ 0,\;\;t \geqslant 112. \\ \end{array} \right.$

${w_{{\rm{wheel1}}}} = \left\{ \begin{array}{l} {w_{{\rm{wheel1}}}},\;\;t < 112; \\ 0,\;\;t \geqslant 112. \\ \end{array} \right.$

${{{s}}_{{\rm{Magm}}}} = \left\{ \begin{array}{l} {{{s}}_{{\rm{Magm}}}},\;\;t < 112; \\ 1.2{{{s}}_{{\rm{Magm}}}},\;\;t \geqslant 112. \end{array} \right.$

式中： ${f_{{\rm{gyro}}\_x}}$表示陀螺x轴突变故障； ${f_{{\rm{star}}\_x}}$表示星敏的滚转角测量存在偏差故障； ${{{s}}_{{\rm s}{\rm{un}}}}$、 ${{{s}}_{{\rm{Magm}}}}$分别表示太敏、磁强计输出信号； ${w_{{\rm{wheel1}}}}$为1号飞轮反馈的角速度信号.

场景1中的系统层检测结果如图3所示. 图中，r₁为观测器残差模值. 由图3可知，r₁在第51 s超出阈值，且持续增大，表明系统中存在故障，进一步对器件层进行检测.

器件层检测结果如图4所示. 图中，r₂为定姿算法的新息模值，r₃为数字陀螺的新息模值，s_x、w_x分别为星敏滚转角、陀螺滚转角速度信号中高频部分经小波重构得到的信号. 由图4可知，r₂从第51 s起新息发散，表明陀螺和星敏至少有一个发生故障. 分别对星敏和陀螺从当前时刻往前取32个采样数据进行小波分析，高频部分的重构信号均有较明显的突变. 计算AD统计量得到AD_star=5.077 9，AD_gyro=1.979 6，两者均大于临界值，说明星敏和陀螺在第51 s时均出现故障，与场景1设置的故障情况相符.

由于此时陀螺采样值异常，导致数字陀螺新息超出阈值，无法判断执行机构的工作状态. 根据器件检测结果，立即切换至备陀螺，采用太敏、磁强计定姿，可以看到数字陀螺新息在短暂超出阈值后趋于正常，表明执行机构无故障.

场景2中，假设星敏不可用，采用陀螺、太敏及磁强计定姿. 系统层检测结果如图5所示，在第112 s时，观测器残差超出阈值，进入器件层检测.

场景2下的器件层检测结果如图6所示. 图中，sun_x、m_x分别为太敏、磁强计X轴采样信号中高频部分经小波重构得到的信号. 此时，定姿算法新息超出阈值，表明至少存在一个敏感器故障. 对陀螺、磁强计、太敏采样数据进行小波分析，计算出统计量分别为：AD_Magm=6.944 9，AD_gyro=0.514 1，AD_sun=1.866 0，可知在第112 s时陀螺仪正常，太敏与磁强计均故障，与场景2设置的故障情况相符. 根据数字陀螺仿真结果可知，在第115 s时新息超出阈值，对陀螺进行小波分析，计算得到AD_gyro=0.251 2，小于临界值，表明陀螺仪正常，则可以推出执行机构故障，至此故障全部检出.

为了进一步验证方案的有效性，以磁强计为例，验证系统对器件不同故障模式的检测性能. 分别对输出卡死、恒偏差、恒增益故障进行300次仿真验证，与常规的小波阈值检测方法对比，此处选取局部阈值（即对不同时间窗内的高频小波系数选取对应的阈值）. 2种方案的检测统计结果如表3所示.

根据表3可以计算各故障下算法检测的准确率、误检率以及漏检率，结果如表4所示. 对比2种方式可知，分层检测的准确率有一定程度的提高，误检率大幅度降低. 常规的阈值检测几乎无法检出输出卡死故障，这是由于该类故障发生后，小波系数将逐渐减小，一直在阈值内，但因其高频小波系数不再服从正态分布，可以由假设检验判断出故障.

本文针对皮纳卫星在轨故障检测，提出分层故障检测方案. 系统层监测全局故障状况，器件层完成故障定位与隔离. 相比于常规检测方法，层次的划分能够有效地降低误检率，提高检测准确性，并且能够适应更多的故障类型，可以实现多器件同时故障的检测. 系统层的监测能够在未发生故障时避免进行不必要的运算，节省计算资源. 在器件层中，利用假设检验分析小波重构系数，无需利用历史信息计算阈值，适用于历史数据少的情况. 通过仿真分析，该方案能够实现星敏、陀螺等器件同时故障的实时检测与定位，并且对磁强计突变、恒偏差、恒增益、输出卡死类故障的检测准确率达到92%以上，误检率低于2%，对皮纳卫星的在线故障诊断具有一定的借鉴价值. 由于小波分析对缓变类故障的检测效率低，当多个器件发生缓变故障时，提出的方案尚无法有效地实现故障定位. 如何提高对缓变故障的检测率，是今后研究的重点.