CN103676941B - 基于运动学和动力学模型的卫星控制系统故障诊断方法 - Google Patents

Abstract

基于运动学和动力学模型的卫星控制系统故障诊断方法，能够在硬件冗余条件不满足的情况下，实现动量轮、陀螺和星敏感器的故障分离。该方法基本流程为：首先，建立用于描述执行机构、惯性敏感器和光学敏感器之间输入输出关系的解析模型；然后，基于模型分别设计两个状态估计器，第1个状态估计器产生的残差对执行机构故障和惯性敏感器故障都敏感，不受光学敏感器故障的影响；而第2个状态估计器的残差对惯性敏感器故障和光学敏感器故障都敏感，不受执行机构故障的影响。这样，通过同时检测两个状态估计器的残差并与在线计算的阈值相比较，可以区分执行机构故障、惯性敏感器故障和光学敏感器故障。

Description

基于运动学和动力学模型的卫星控制系统故障诊断方法

技术领域

本发明属于卫星姿态控制领域，涉及一种基于运动学和动力学模型的卫星控制系统执行机构和敏感器故障诊断方法。

背景技术

现代航天器大多结构复杂，工作环境恶劣，运行时间长，在轨工作过程中极易发生故障。一旦发生故障，若不采取积极有效的措施进行应对，很可能导致航天任务失败，从而带来巨大的资源浪费和不利的社会影响。

作为构成卫星的各个分系统中最关键也是最复杂的一个分系统，卫星姿态控制系统的可靠性是星上多个分系统正常运行的基本保障。如遥测遥控分系统要实现与地面站的通信，要求天线指向正确的位置；能源分系统要保障正常的电力供应，要求太阳帆板跟踪太阳；对于对地观测卫星而言，有效载荷的正常工作要求卫星进行高精度高稳定度姿态控制。因此，卫星姿态控制系统能否正常工作往往关系到整个航天任务的成败。

卫星姿态控制系统由姿态敏感器、控制器、执行机构和卫星本体一起构成闭环控制回路。基于姿态敏感器测量和确定卫星相对于空间某些已知基准目标的方位；控制器对测得的信息进一步处理后确定卫星姿态，并根据所确定的姿态按满足设计要求的控制律发出指令，控制执行机构按控制指令产生所需的控制力矩，实现卫星姿态控制。卫星姿态控制系统的故障率较高且危害较大。

为了提高卫星姿态控制系统的可靠性，一方面，应提高其组成部件如各敏感器、执行机构和控制器等自身的可靠性，降低发生故障的可能性；另一方面，应对关键部件采用冗余设计，并通过故障诊断技术消除或减弱故障造成的不利影响。

卫星姿态敏感器包括光学敏感器和惯性敏感器，其中，光学敏感器有太阳敏感器、地球敏感器等，惯性敏感器以陀螺为代表。一般采用一致性检验的方法进行卫星姿态敏感器故障诊断，如通过比较两个光学敏感器输出的方式进行故障检测，但当硬件冗余条件不满足时（如无故障光学敏感器个数少于2个），不能采用这种方法进行故障诊断。为了解决硬件冗余条件不满足情况下卫星姿态控制系统执行机构与敏感器的故障诊断问题，本发明提出一种基于状态估计的执行机构、惯性敏感器和光学敏感器故障分离方法。

发明内容

本发明的技术解决问题是：克服现有技术的不足，提供了一种运动学和动力学模型的卫星姿态控制系统执行机构与敏感器的故障诊断方法，能够在硬件冗余条件不满足的情况下，实现执行机构、惯性敏感器和光学敏感器故障分离。

本发明的技术解决方案是：

基于运动学和动力学模型的卫星控制系统故障分离方法，步骤如下：

（1）建立卫星动力学子系统故障诊断模型和卫星运动学子系统故障诊断模型；

（2）根据步骤（1）中的动力学子系统故障诊断模型设计基于状态估计的残差生成器，计算残差评价值J_FD1,r(k)与相应的阈值J_FD1,th(k)；

（3）根据步骤（1）中的运动学子系统故障诊断模型设计基于状态估计的残差生成器，计算残差评价值J_FD2,r(k)与相应的阈值J_FD2,th(k)；

（4）如果J_FD1,r(k)>J_FD1,th(k)并且J_FD2,r(k)≤J_FD2,th(k)，则卫星姿态控制系统的执行机构发生了故障；如果J_FD1,r(k)>J_FD1,th(k)并且J_FD2,r(k)>J_FD2,th(k)，则卫星姿态控制系统的惯性敏感器发生了故障；如果J_FD1,r(k)≤J_FD1,th(k)并且J_FD2,r(k)>J_FD2,th(k)，则卫星姿态控制系统的光学敏感器发生了故障；否则，卫星姿态控制系统的执行机构和敏感器没有发生故障。

所述步骤（1）中建立卫星动力学子系统故障诊断模型和卫星运动学子系统故障诊断模型，卫星动力学子系统故障诊断模型如下所示：

$\left\{\begin{array}{c}{x}_g(k+1)=A_{g,d}{x}_g\left(k\right)+B_{g.d}u\left(k\right)\\ y_g\left(k\right)=C_g{x}_g\left(k\right)\end{array}\right.$

其中，k表示离散的时间，x_g(k)表示第k步的动力学子系统状态量卫星本体角速度，A_g,d表示动力学子系统状态矩阵，B_g,d表示动力学子系统输入矩阵，u(k)表示动力学子系统输入量，y_g(k)表示动力学子系统观测量，C_g表示动力学子系统观测矩阵；

卫星运动学子系统故障诊断模型如下所示：

$\left\{\begin{array}{c}q(k+1)=\stackrel{\‾}{\mathrm{\Ω}}\left(\mathrm{\ω}\left(k\right)\right)q\left(k\right)\\ s_{b,m}\left(k\right)=\frac{\tilde{q}\left(k\right)\⊗s_i\left(k\right)\⊗q\left(k\right)+n_s\left(k\right)}{\left|\right|\tilde{q}\left(k\right)\⊗s_i\left(k\right)\⊗q\left(k\right)+n_s\left(k\right)\left|\right|}\end{array}\right.$

其中，q(k)表示第k步的运动学子系统状态量姿态四元数，为q(k)的共轭，ω(k)表示运动学子系统输入量，表示运动学子系统的与输入量相关的时变状态矩阵，s_b,m(k)表示第k步的运动学子系统观测量，s_i(k)表示与s_b,m(k)相对应的已知量，n_s(k)表示服从高斯正态分布的测量噪声，为n_s(k)的方差，符号“～”、分别表示四元数的共轭与四元数间的乘法；

所述步骤（2）中根据步骤（1）中的动力学子系统故障诊断模型设计基于状态估计的残差生成器，计算残差评价值与相应的阈值具体为：

基于状态估计的动力学子系统故障诊断残差生成器为：

$\left\{\begin{array}{c}{\hat{x}}_g(k+1)=A_{g,d}{\hat{x}}_g\left(k\right)+B_{g,d}u\left(k\right)+L_g(y_g\left(k\right)-C_g{\hat{x}}_g\left(k\right))\\ r_{W+G}\left(k\right)=y_g\left(k\right)-C_g{\hat{x}}_g\left(k\right)\end{array}\right.$

其中，表示动力学子系统状态量的估计值，L_g表示动力学子系统残差生成器的反馈增益阵，r_W+G(k)表示动力学子系统的残差；

动力学子系统的残差评价值通过如下方式计算：

J_FD1,r(k)=||r_W+G(k)||

其中，J_FD1,r(k)表示动力学子系统的残差评价值；

与动力学子系统的残差评价值相对应的阈值计算方式为：

J_FD1,th(k)=3·trace(σ_a)

其中，J_FD1,th(k)表示动力学子系统的阈值，σ_a表示观测量y_g(k)的测量均方差。

所述步骤（3）中根据步骤（1）中的运动学子系统故障诊断模型设计基于状态估计的残差生成器，计算残差评价值与相应的阈值具体为：

运动学子系统的残差值计算方式如下：

$r_{G+S}\left(k\right)=s_{\mathrm{bm}}\left(k\right)-\stackrel{\widetilde{^}}{q}\left(k\right)\⊗s_i\left(k\right)\⊗\hat{q}\left(k\right)$

其中，表示由无迹卡尔曼滤波算法得到的运动学子系统第k步的状态量估计值，r_G+S(k)表示运动学子系统的残差；

运动学子系统的残差评价值计算方式如下：

J_FD2,r(k)=||r_G+S(k)||

其中，J_FD2,r(k)表示运动学子系统的残差评价值；

与运动学子系统的残差评价值相对应的阈值计算方式如下：

$J_{\mathrm{FD}2,\mathrm{th}}\left(k\right)=3\sqrt{\mathrm{trace}\left(P_{k+1}^{\mathrm{ss}}\left(k\right)\right)}$

其中，J_FD2,th(k)表示运动学子系统的阈值，表示无迹卡尔曼滤波算法中新息的估计方差。

本发明与现有技术相比的有益效果是：

国内卫星多采用基于硬件冗余的卫星姿态敏感器故障诊断方式，但当硬件冗余条件不满足时（如光学敏感器个数少于2个，陀螺个数少于4个），不能采用这种方法进行故障诊断。针对这一问题，本发明提出一种基于状态估计的卫星姿态控制系统输入输出型故障分离方法，该方法能够在硬件冗余条件不满足的情况下，实现执行机构、惯性敏感器和光学敏感器的故障分离。

附图说明

图1为本发明流程图；

图2为无故障情况下故障检测曲线；

图3为动量轮恒值偏差故障情况下故障检测曲线；

图4为动量轮卡死故障情况下故障检测曲线；

图5为动量轮缓变偏差故障情况下故障检测曲线；

图6为陀螺漂移增大故障情况下故障检测曲线；

图7为陀螺精度降低故障情况下故障检测曲线；

图8为星敏感器恒值偏差故障情况下故障检测曲线；

图9为星敏感器精度降低故障情况下故障检测曲线；

具体实施方式

下面结合附图对本发明的具体实施方式进行进一步的详细描述。

如图1所示，本发明提出一种基于运动学和动力学模型的卫星控制系统敏感器和执行机构故障分离方法，步骤如下：

（1）建立用于描述星载的执行机构、惯性敏感器和光学敏感器之间输入输出关系的故障诊断模型，即卫星动力学子系统故障诊断模型和卫星运动学子系统故障诊断模型。假设卫星处于三轴稳定的工作模式下，此时卫星的姿态角和姿态角速度都是小量，

对卫星姿态动力学公式在平衡位置近似线性化、离散化，得到动力学子系统故障诊断模型形式如下所示：

x_g(k+1)=A_g,dx_g(k)+B_g,du(k)

y_g(k)=C_gx_g(k)

其中，k表示离散的时间，x_g(k)表示第k步的动力学子系统状态量卫星本体角速度，A_g,d表示动力学子系统状态矩阵，B_g,d表示动力学子系统输入矩阵，u(k)表示动力学子系统输入量，y_g(k)表示动力学子系统观测量，C_g表示动力学子系统观测矩阵；

卫星运动学子系统故障诊断模型形式如下所示：

$q(k+1)=\stackrel{\‾}{\mathrm{\Ω}}\left(\mathrm{\ω}\left(k\right)\right)q\left(k\right)$

$s_{b,m}\left(k\right)=\frac{\tilde{q}\left(k\right)\⊗s_i\left(k\right)\⊗q\left(k\right)+n_s\left(k\right)}{\left|\right|\tilde{q}\left(k\right)\⊗s_i\left(k\right)\⊗q\left(k\right)+n_s\left(k\right)\left|\right|}$

其中，符号“～”、分别表示四元数的共轭与四元数间的乘法，q(k)表示第k步的运动学子系统状态量姿态四元数，ω(k)表示运动学子系统输入量，表示运动学子系统的与输入量相关的时变状态矩阵，

$\stackrel{\‾}{\mathrm{\Ω}}\left(\mathrm{\ω}\right)=\left[\begin{array}{cccc}\cos \left(\frac{\left|\right|\mathrm{\ω}\left|\right|\mathrm{\Δt}}{2}\right)& -\frac{{\mathrm{\ω}}_1}{\left|\right|\mathrm{\ω}\left|\right|}\sin \left(\frac{\left|\right|\mathrm{\ω}\left|\right|\mathrm{\Δt}}{2}\right)& -\frac{{\mathrm{\ω}}_2}{\left|\right|\mathrm{\ω}\left|\right|}\sin \left(\frac{\left|\right|\mathrm{\ω}\left|\right|\mathrm{\Δt}}{2}\right)& -\frac{{\mathrm{\ω}}_3}{\left|\right|\mathrm{\ω}\left|\right|}\sin \left(\frac{\left|\right|\mathrm{\ω}\left|\right|\mathrm{\Δt}}{2}\right)\\ \frac{{\mathrm{\ω}}_1}{\left|\right|\mathrm{\ω}\left|\right|}\sin \left(\frac{\left|\right|\mathrm{\ω}\left|\right|\mathrm{\Δt}}{2}\right)& \cos \left(\frac{\left|\right|\mathrm{\ω}\left|\right|\mathrm{\Δt}}{2}\right)& \frac{{\mathrm{\ω}}_3}{\left|\right|\mathrm{\ω}\left|\right|}\sin \left(\frac{\left|\right|\mathrm{\ω}\left|\right|\mathrm{\Δt}}{2}\right)& -\frac{{\mathrm{\ω}}_2}{\left|\right|\mathrm{\ω}\left|\right|}\sin \left(\frac{\left|\right|\mathrm{\ω}\left|\right|\mathrm{\Δt}}{2}\right)\\ \frac{{\mathrm{\ω}}_2}{\left|\right|\mathrm{\ω}\left|\right|}\sin \left(\frac{\left|\right|\mathrm{\ω}\left|\right|\mathrm{\Δt}}{2}\right)& -\frac{{\mathrm{\ω}}_3}{\left|\right|\mathrm{\ω}\left|\right|}\sin \left(\frac{\left|\right|\mathrm{\ω}\left|\right|\mathrm{\Δt}}{2}\right)& \cos \left(\frac{\left|\right|\mathrm{\ω}\left|\right|\mathrm{\Δt}}{2}\right)& \frac{{\mathrm{\ω}}_1}{\left|\right|\mathrm{\ω}\left|\right|}\sin \left(\frac{\left|\right|\mathrm{\ω}\left|\right|\mathrm{\Δt}}{2}\right)\\ \frac{{\mathrm{\ω}}_3}{\left|\right|\mathrm{\ω}\left|\right|}\sin \left(\frac{\left|\right|\mathrm{\ω}\left|\right|\mathrm{\Δt}}{2}\right)& \frac{{\mathrm{\ω}}_2}{\left|\right|\mathrm{\ω}\left|\right|}\sin \left(\frac{\left|\right|\mathrm{\ω}\left|\right|\mathrm{\Δt}}{2}\right)& -\frac{{\mathrm{\ω}}_1}{\left|\right|\mathrm{\ω}\left|\right|}\sin \left(\frac{\left|\right|\mathrm{\ω}\left|\right|\mathrm{\Δt}}{2}\right)& \cos \left(\frac{\left|\right|\mathrm{\ω}\left|\right|\mathrm{\Δt}}{2}\right)\end{array}\right]$

ω=[ω₁ω₂ω₃]^T量，Δt表示采样时间，s_b,m(k)表示第k步的运动学子系统观测量，s_i(k)表示与s_b,m(k)相对应的已知量，表示测量噪声，符号“～”、分别表示四元数的共轭与四元数间的乘法；

（2）根据步骤（1）中的动力学子系统故障诊断模型设计基于状态估计的残差生成器，计算残差评价值与相应的阈值；

基于状态估计的动力学子系统故障诊断残差生成器形式为：

${\hat{x}}_g(k+1)=A_{g,d}{\hat{x}}_g\left(k\right)+B_{g,d}u\left(k\right)+L_g(y_g\left(k\right)-C_g{\hat{x}}_g\left(k\right))$

$r_{W+G}\left(k\right)=y_g\left(k\right)-C_g{\hat{x}}_g\left(k\right)$

其中，表示动力学子系统状态量的估计值，L_g表示动力学子系统残差生成器的反馈增益阵，r_W+G(k)表示动力学子系统的残差；

动力学子系统的残差评价值计算方式如下所示：

J_FD1,r(k)=||r_W+G(k)||

其中，J_FD1,r(k)表示动力学子系统的残差评价值；

与动力学子系统的残差评价值相对应的阈值计算方式如下所示：

J_FD1,th(k)=3·trace(σ_a)

其中，J_FD1,th(k)表示动力学子系统的阈值，σ_a表示观测量y_g(k)的测量均方差；

（3）根据步骤（1）中的运动学子系统故障诊断模型设计基于状态估计的残差生成器，计算残差评价值与相应的阈值；

基于状态估计的运动学子系统故障诊断残差生成器利用了无迹卡尔曼滤波（UKF）算法，该算法的计算过程如下所示：

UKF算法初始化

给定

赋值 ${\hat{x}}_0^{+}=\left[\begin{array}{c}0\\ b_0^{+}\end{array}\right],$ a，λ，f=2(a+1)， $\mathrm{\δp}=f\frac{\mathrm{\δ}\stackrel{\→}{q}}{a+{\mathrm{\δq}}_0};$

${\stackrel{\‾}{Q}}_k=\frac{\mathrm{\Δt}}{2}\left[\begin{array}{cc}{\mathrm{\σ}}_a^2-\frac{1}{6}{\mathrm{\σ}}_b^2\mathrm{\Δt}& 0_{3\×3}\\ 0_{3\×3}& {\mathrm{\σ}}_b^2\end{array}\right];$

UKF算法迭代过程

${\mathrm{\σ}}_k\←2\mathrm{ncolumnsfrom}\±\sqrt{(n+\mathrm{\λ})(P_k^{+}+{\stackrel{\‾}{Q}}_k)}$

${\mathrm{\χ}}_k\left(0\right)={\hat{x}}_k^{+}$

${\mathrm{\χ}}_k\left(i\right)={\mathrm{\σ}}_k\left(i\right)+{\hat{x}}_k^{+}$

$\mathrm{\δ}{\hat{q}}_k^{+}\left(i\right)=\left[\begin{array}{c}{\mathrm{\δq}}_{0,k}^{+}\left(i\right)\\ \mathrm{\δ}{\stackrel{\→}{q}}_k^{+}\left(i\right)\end{array}\right]$

${\mathrm{\δq}}_{0,k}^{+}\left(i\right)=\frac{-a{\left|\right|{\mathrm{\χ}}_k^{\mathrm{\δp}}\left(i\right)\left|\right|}^2+f\sqrt{f^2+(1-a^2){\left|\right|{\mathrm{\χ}}_k^{\mathrm{\δp}}\left(i\right)\left|\right|}^2}}{f^2+{\left|\right|{\mathrm{\χ}}_k^{\mathrm{\δp}}\left(i\right)\left|\right|}^2},i=1,...,12$

$\mathrm{\δ}{\stackrel{\→}{q}}_k^{+}\left(i\right)=(a+{\mathrm{\δq}}_{0,k}^{+}\left(i\right)){\mathrm{\χ}}_k^{\mathrm{\δp}}\left(i\right)/f,i=1,...,12$

${\hat{q}}_k^{+}\left(0\right)={\hat{q}}_k^{+}$

${\hat{q}}_k^{+}\left(i\right)=\mathrm{\δ}{\hat{q}}_k^{+}\left(i\right)\⊗{\hat{q}}_k^{+},i=1,...,12$

${\mathrm{\δq}}_{k+1}^{-}\left(i\right)={\hat{q}}_{k+1}^{-}\left(i\right)\⊗{\left({\hat{q}}_{k+1}^{-}\left(0\right)\right)}^{-1},i=\mathrm{0,1},...,12$

${\mathrm{\χ}}_{k+1}^{\mathrm{\δp}}\left(0\right)=0$

${\mathrm{\χ}}_{k+1}^{\mathrm{\δp}}\left(i\right)=f\frac{\mathrm{\δ}{\stackrel{\→}{q}}_{k+1}^{-}\left(i\right)}{a+\mathrm{\δ}{q}_{0,k+1}^{-}\left(i\right)},i=1,...,12$

${\mathrm{\χ}}_{k+1}^b\left(i\right)={\mathrm{\χ}}_k^b\left(i\right),i=\mathrm{0,1},...,12$

${\hat{x}}_{k+1}^{-}=\frac{1}{n+\mathrm{\λ}}(\mathrm{\λ}{\mathrm{\χ}}_{k+1}\left(0\right)+\frac{1}{2}\underset{i=1}{\overset{2n}{\mathrm{\Σ}}}{\mathrm{\χ}}_{k+1}\left(i\right))$

$P_{k+1}^{-}=\frac{1}{n+\mathrm{\λ}}[\mathrm{\λ}({\mathrm{\χ}}_{k+1}\left(0\right)-{\hat{x}}_{k+1}^{-}){({\mathrm{\χ}}_{k+1}\left(0\right)-{\hat{x}}_{k+1}^{-})}^T$

$+\frac{1}{2}\underset{i=1}{\overset{2n}{\mathrm{\Σ}}}({\mathrm{\χ}}_{k+1}\left(i\right)-{\hat{x}}_{k+1}^{-}){({\mathrm{\χ}}_{k+1}\left(i\right)-{\hat{x}}_{k+1}^{-})}^T]+{\stackrel{\‾}{Q}}_k$

${\mathrm{\γ}}_{k+1}\left(i\right)=h({\mathrm{\χ}}_{k+1}\left(i\right),k)={\left[\begin{array}{c}A\left({\hat{q}}^{-}\left(i\right)\right)s_1\\ A\left({\hat{q}}^{-}\left(i\right)\right)s_2\end{array}\right]}_{k+1},i=\mathrm{0,1},...,12$

${\hat{y}}_{k+1}^{-}=\frac{1}{n+\mathrm{\λ}}({\mathrm{\λ\γ}}_{k+1}\left(0\right)+\frac{1}{2}\underset{i=1}{\overset{2n}{\mathrm{\Σ}}}{\mathrm{\γ}}_{k+1}\left(i\right))$

$P_{k+1}^{\mathrm{yy}}=\frac{1}{n+\mathrm{\λ}}[\mathrm{\λ}({\mathrm{\γ}}_{k+1}\left(0\right)-{\hat{y}}_{k+1}^{-}){({\mathrm{\γ}}_{k+1}\left(0\right)-{\hat{y}}_{k+1}^{-})}^T$

$+\frac{1}{2}\underset{i=1}{\overset{2n}{\mathrm{\Σ}}}({\mathrm{\γ}}_{k+1}\left(i\right)-{\hat{y}}_{k+1}^{-}){({\mathrm{\γ}}_{k+1}\left(i\right)-{\hat{y}}_{k+1}^{-})}^T]$

$P_{k+1}^{\mathrm{ss}}=P_{k+1}^{\mathrm{yy}}+R_{k+1}$

$P_{k+1}^{\mathrm{xy}}=\frac{1}{n+\mathrm{\λ}}[\mathrm{\λ}({\mathrm{\χ}}_{k+1}\left(0\right)-{\hat{x}}_{k+1}^{-}){({\mathrm{\γ}}_{k+1}\left(0\right)-{\hat{y}}_{k+1}^{-})}^T$

$+\frac{1}{2}\underset{i=1}{\overset{2n}{\mathrm{\Σ}}}({\mathrm{\χ}}_{k+1}\left(i\right)-{\hat{x}}_{k+1}^{-}){({\mathrm{\γ}}_{k+1}\left(i\right)-{\hat{y}}_{k+1}^{-})}^T]$

$K_k=P_k^{\mathrm{xy}}{\left(P_k^{\mathrm{ss}}\right)}^{-1}$

${\hat{x}}_k^{+}={\hat{x}}_k^{-}+K_k({\tilde{y}}_k-{\hat{y}}_k^{-})$

$P_k^{+}=P_k^{+}-K_k{P}_k^{\mathrm{ss}}{K}_k^T$

${\mathrm{\δq}}_{0,k+1}^{+}\left(i\right)=\frac{-a{\left|\right|\mathrm{\δ}{\hat{p}}_{k+1}^{+}\left|\right|}^2+f\sqrt{f^2+(1-a^2){\left|\right|\mathrm{\δ}{\hat{p}}_{k+1}^{+}\left|\right|}^2}}{f^2+{\left|\right|\mathrm{\δ}{\hat{p}}_{k+1}^{+}\left|\right|}^2}$

$\mathrm{\δ}{\stackrel{\→}{q}}_{k+1}^{+}\left(i\right)=(a+\mathrm{\δ}{q}_{0,k+1}^{+}\left(i\right))\mathrm{\δ}{\hat{p}}_{k+1}^{+}/f$

$\mathrm{\δ}{\hat{q}}_{k+1}^{+}=\left[\begin{array}{c}\mathrm{\δ}{q}_{0,k+1}^{+}\\ \mathrm{\δ}{\stackrel{\→}{q}}_{k+1}^{+}\end{array}\right]$

${\hat{q}}_{k+1}^{+}=\mathrm{\δ}{\hat{q}}_{k+1}^{+}\⊗{\hat{q}}_{k+1}^{-}\left(0\right)$

置零，准备下一循环

运动学子系统的残差值计算方式如下所示：

$r_{G+S}\left(k\right)=s_{\mathrm{bm}}\left(k\right)-\stackrel{\widetilde{^}}{q}\left(k\right)\⊗s_i\left(k\right)\⊗\hat{q}\left(k\right)$

其中，表示由UKF算法得到的运动学子系统第k步的状态量估计值，r_G+S(k)表示运动学子系统的残差；

运动学子系统的残差评价值计算方式如下所示：

J_FD2,r(k)=||r_G+S(k)||

其中，J_FD2,r(k)表示运动学子系统的残差评价值；

与运动学子系统的残差评价值相对应的阈值计算方式如下所示：

$J_{\mathrm{FD}2,\mathrm{th}}\left(k\right)=3\sqrt{\mathrm{trace}\left(P_{k+1}^{\mathrm{ss}}\left(k\right)\right)}$

其中，J_FD2,th(k)表示运动学子系统的阈值，表示UKF算法中新息的估计方差；

（4）如果J_FD1,r(k)>J_FD1,th(k)并且J_FD2,r(k)≤J_FD2,th(k)，则卫星姿态控制系统的执行机构发生了故障；如果J_FD1,r(k)>J_FD1,th(k)并且J_FD2,r(k)>J_FD2,th(k)，则卫星姿态控制系统的惯性敏感器发生了故障；如果J_FD1,r(k)≤J_FD1,th(k)并且J_FD2,r(k)>J_FD2,th(k)，则卫星姿态控制系统的光学敏感器发生了故障；否则，卫星姿态控制系统的执行机构和敏感器没有发生故障。

下面，以某低轨卫星姿态控制系统为诊断对象，通过仿真实例验证本发明所述方法的有效性。卫星执行机构包括三个沿卫星本体轴方向安装的动量轮、一组测量三轴角速度的陀螺和沿滚动轴、俯仰轴方向安装的星敏感器。故障诊断模型和残差生成器参数如下。

$A_{g,d}=\left[\begin{array}{ccc}1& 0& 7\×{10}^{-5}\\ 0& 1& 0\\ -1.429\×{10}^4& 0& 1\end{array}\right]$

$B_{g,d}=\left[\begin{array}{ccc}0.0001& 0& 5\×{10}^{-9}\\ 0& 0.0002& 0\\ -7.143\×{10}^{-9}& 0& 0.0001429\end{array}\right]$

$C_g=\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right]$

$L_g=\left[\begin{array}{ccc}0.9& 0& 0\\ 0& 0.9& 0\\ 0& 0& 0.9\end{array}\right]$

对以下8种情况进行仿真研究，分别是：

（a）无故障情况；

（b）动量轮故障1，俯仰轴动量轮输出力矩出现恒值偏差，偏差值为-0.05Nm；

（c）动量轮故障2，俯仰轴动量轮卡死，输出力矩等于外部扰动力矩；

（d）动量轮故障3，俯仰轴动量轮输出力矩出现缓变偏差，偏差值由0逐渐增长到-0.05Nm；

（e）陀螺故障4，俯仰轴陀螺漂移增大，偏差值为0.0001rad/s；

（f）陀螺故障5，俯仰轴陀螺测量精度降低，测量噪声标准差由10^-6rad/s变为10^-4rad/s；

（g）星敏感器故障6，滚动轴方向星敏感器测量输出出现偏差，偏差值为-0.0001，约20″；

（h）星敏感器故障7，滚动轴方向星敏感器测量精度逐渐降低，测量噪声标准差由10^-5逐渐增长到10^-3。

分别对以上8种情况进行仿真，卫星控制系统姿态动力学子系统和运动学子系统的故障检测曲线如图2-图9所示。图中，上半部分为动力学子系统的故障检测残差评价值与阈值曲线，下半部分为运动学子系统的故障检测残差评价值与阈值曲线。

由图2可以看到，无故障发生时，两组残差评价值均没有超过各自的阈值，没有出现误报；由图3、图4和图5可以看到，在动量轮输出力矩出现恒值偏差、卡死和输出力矩出现缓变偏差3种故障情况下，第一组残差评价值在故障发生后的一定时间内均超过了相应的阈值，而第二组残差评价值均没有超过各自的阈值，说明第一组残差对动量轮故障敏感，而第二组残差对动量轮故障解耦；由图6、图7可以看到，在陀螺漂移增大和测量精度降低2种故障情况下，两组残差评价值在故障发生后的一定时间内均超过了相应的阈值，说明两组残差均对陀螺故障敏感；由图8、图9可以看到，在星敏感器测量输出出现偏差和测量精度逐渐降低2种故障情况下，第一组残差评价值均没有超过各自的阈值，而第二组残差评价值在故障发生后的一定时间内均超过了相应的阈值，说明第一组残差对星敏感器故障解耦，而第二组残差对星敏感器故障敏感。

这样，通过综合分析两组残差评价值是否超过各自的阈值，便可以根据步骤（4）中的故障分离逻辑，初步实现动量轮、陀螺和星敏感器的故障分离，为进一步实现故障的准确定位提供重要的信息。

仿真研究表明，所提方法可以实现对执行机构、惯性敏感器和光学敏感器故障的分离。相对于现有的基于硬件冗余的敏感器故障诊断方式，该方法对硬件冗余要求低，易于星上实现。本项专利的主要技术内容可用于各类卫星姿态控制系统故障诊断，具有广阔的应用前景。

本发明说明书中未作详细描述的内容属于本领域专业技术人员的公知技术。

Claims (3)

1.基于运动学和动力学模型的卫星控制系统故障诊断方法，其特征在于步骤如下：

(1)建立卫星动力学子系统故障诊断模型和卫星运动学子系统故障诊断模型；

(2)根据步骤(1)中的动力学子系统故障诊断模型设计基于状态估计的残差生成器，计算残差评价值J_FD1,r(k)与相应的阈值J_FD1,th(k)；

(3)根据步骤(1)中的运动学子系统故障诊断模型设计基于状态估计的残差生成器，计算残差评价值J_FD2,r(k)与相应的阈值J_FD2,th(k)；

(4)如果J_FD1,r(k)>J_FD1,th(k)并且J_FD2,r(k)≤J_FD2,th(k)，则卫星姿态控制系统的执行机构发生了故障；如果J_FD1,r(k)>J_FD1,th(k)并且J_FD2,r(k)>J_FD2,th(k)，则卫星姿态控制系统的惯性敏感器发生了故障；如果J_FD1,r(k)≤J_FD1,th(k)并且J_FD2,r(k)>J_FD2,th(k)，则卫星姿态控制系统的光学敏感器发生了故障；否则，卫星姿态控制系统的执行机构和敏感器没有发生故障；

所述步骤(1)中建立卫星动力学子系统故障诊断模型和卫星运动学子系统故障诊断模型，卫星动力学子系统故障诊断模型如下所示：

$\left\{\begin{array}{c}{x}_g(k+1)=A_{g,d}{x}_g\left(k\right)+B_{g,d}u\left(k\right)\\ y_g\left(k\right)=C_g{x}_g\left(k\right)\end{array}\right.$

其中，k表示离散的时间，x_g(k)表示第k步的动力学子系统状态量卫星本体角速度，A_g,d表示动力学子系统状态矩阵，B_g,d表示动力学子系统输入矩阵，u(k)表示动力学子系统输入量，y_g(k)表示动力学子系统观测量，C_g表示动力学子系统观测矩阵；

卫星运动学子系统故障诊断模型如下所示：

$\left\{\begin{array}{c}q(k+1)=\stackrel{\‾}{\mathrm{\Ω}}\left(\mathrm{\ω}\left(k\right)\right)q\left(k\right)\\ s_{b,m}\left(k\right)=\frac{\tilde{q}\left(k\right)\⊗s_i\left(k\right)\⊗q\left(k\right)+n_s\left(k\right)}{||\tilde{q}\left(k\right)\⊗s_i\left(k\right)\⊗q\left(k\right)+n_s\left(k\right)||}\end{array}\right.$

其中，q(k)表示第k步的运动学子系统状态量姿态四元数，为q(k)的共轭，ω(k)表示运动学子系统输入量，表示运动学子系统的与输入量相关的时变状态矩阵，s_b,m(k)表示第k步的运动学子系统观测量，s_i(k)表示与s_b,m(k)相对应的已知量，n_s(k)表示服从高斯正态分布的测量噪声，为n_s(k)的方差，符号“～”、分别表示四元数的共轭与四元数间的乘法。

2.根据权利要求1所述的基于运动学和动力学模型的卫星控制系统故障诊断方法，其特征在于：所述步骤(2)中根据步骤(1)中的动力学子系统故障诊断模型设计基于状态估计的残差生成器，计算残差评价值与相应的阈值具体为：

基于状态估计的动力学子系统故障诊断残差生成器为：

$\left\{\begin{array}{c}{\hat{x}}_g(k+1)=A_{g,d}{\hat{x}}_g\left(k\right)+B_{g,d}u\left(k\right)+L_g(y_g\left(k\right)-C_g{\hat{x}}_g\left(k\right))\\ r_{W+G}\left(k\right)=y_g\left(k\right)-C_g{\hat{x}}_g\left(k\right)\end{array}\right.$

其中，表示动力学子系统状态量的估计值，L_g表示动力学子系统残差生成器的反馈增益阵，r_W+G(k)表示动力学子系统的残差；

动力学子系统的残差评价值通过如下方式计算：

J_FD1,r(k)＝||r_W+G(k)||

其中，J_FD1,r(k)表示动力学子系统的残差评价值；

与动力学子系统的残差评价值相对应的阈值计算方式为：

J_FD1,th(k)＝3·trace(σ_a)

其中，J_FD1,th(k)表示动力学子系统的阈值，σ_a表示观测量y_g(k)的测量均方差。

3.根据权利要求1所述的基于运动学和动力学模型的卫星控制系统故障诊断方法，其特征在于：所述步骤(3)中根据步骤(1)中的运动学子系统故障诊断模型设计基于状态估计的残差生成器，计算残差评价值与相应的阈值具体为：

运动学子系统的残差值计算方式如下：

$r_{G+S}\left(k\right)=s_{bm}\left(k\right)-\tilde{\hat{q}}\left(k\right)\⊗s_i\left(k\right)\⊗\hat{q}\left(k\right)$

其中，表示由无迹卡尔曼滤波算法得到的运动学子系统第k步的状态量估计值，r_G+S(k)表示运动学子系统的残差；

运动学子系统的残差评价值计算方式如下：

J_FD2,r(k)＝||r_G+S(k)||

其中，J_FD2,r(k)表示运动学子系统的残差评价值；

与运动学子系统的残差评价值相对应的阈值计算方式如下：

$J_{FD2,th}\left(k\right)=3\sqrt{trace\left(P_{k+1}^{ss}\right(k\left)\right)}$

其中，J_FD2,th(k)表示运动学子系统的阈值，表示无迹卡尔曼滤波算法中新息的估计方差。