CN103884359B - 一种基于主元分析算法的卫星陀螺部件故障诊断方法 - Google Patents

Abstract

一种基于主元分析算法的卫星陀螺部件故障诊断方法，本发明涉及一种卫星陀螺部件的故障诊断方法，具体涉及一种基于巴特沃斯低通滤波器和主元分析算法的故障诊断方法。本发明是要解决现有主元分析算法进行故障诊断的不足，如出现误报现象的问题。步骤一：基于巴特沃斯低通滤波器对卫星陀螺部件的输出角速度数据X_p进行滤波预处理；步骤二：利用预处理后的卫星陀螺部件的输出角速度数据X，构建PCA数学模型；步骤三：利用步骤二中得到的PCA数学模型参数，采用平方预测误差SPE统计量对残差子空间中的过程数据进行检测；步骤四：检测到故障后根据过程变量贡献图诊断出故障位置。本发明应用于过程数据可观测的卫星陀螺部件运行过程故障诊断应用领域。

Description

一种基于主元分析算法的卫星陀螺部件故障诊断方法

技术领域

本发明涉及一种卫星陀螺部件的故障诊断方法，具体涉及一种基于巴特沃斯低通滤波器和主元分析算法的故障诊断方法。

背景技术

陀螺仪最早是用于航海导航，但随着科学技术的发展，它在航空航天事业中也得到广泛的应用。陀螺仪不仅可以作为指示仪表，还可以作为自动控制系统中的一个敏感元件，即可作为信号传感器。根据需要，陀螺仪能够提供准确的方位、水平、位置、速度和加速度等信号，以便驾驶员或用自动导航仪来控制飞机、舰炮或航天飞机等航行体按一定的航线飞行，而在导弹、卫星运载器或空间探测火箭等航行体的制导中，则直接利用这些信号完成航天体的姿态控制和轨道控制。此外陀螺仪在短期使用时积分误差小，具有较高的精度，且噪声也小，因此在航天器等的姿态测量中得到广泛的应用，是航天器上最为关键的分系统之一，陀螺仪一旦发生故障，将对航天任务的完成产生致命性的影响。

伴随着航天器空间应用的迅速发展，尤其是载人航天器、大型通信卫星、深空探测器等对高可靠性、长寿命日益增长的需求，对星载陀螺仪的性能和可靠性的要求也日益提高。研究卫星陀螺部件故障诊断技术，是发展系统重构和自主管理技术的基础，对保证我国空间任务顺利完成具有重要意义；同时，该技术作为地面试验故障诊断的手段，还可以保证地面试验的安全性，具有重要的现实意义和工程价值。

故障诊断技术已有30多年的发展历史，但作为一门综合性新学科——故障诊断学——还是这些年发展起来的。从不同的角度出发有多种故障诊断分类方法，如基于机理研究的诊断理论和方法、基于信号处理及特征提取的故障诊断方法、模糊诊断理论和方法等，分类方法众多，但可归纳为两类：1）基于非模型的故障诊断理论和方法，如信号空间特征、模态和信息处理方法的诊断理论与方法；2）基于系统数学模型和现代控制理论、方法的故诊断理论与方法。近年来随着相关领域的发展以及理论研究的深入，各种新的诊断方法层出不穷，以全新的角度对现有的故障诊断方法进行了重新分类，整体上可将其分为定量和定性分析两大类。其中，定量分析方法可细分为基于数据驱动和解析模型的方法，基于数据驱动的方法又进一步包括多元统计分析类方法、信息融合类方法等。

鉴于陀螺仪在航天任务中所承担的至关重要的作用，对陀螺仪的故障诊断显得尤为重要，本发明正是针对传统主元分析算法会出现故障误报的缺点，在原有主元分析算法的基础上，基于巴特沃斯低通滤波器对初始数据信息进行数据预处理，以解决传统主元分析算法的不足。

主元分析的基本思想就是在尽量保持初始数据信息的原则下，对由相互之间存在相关性的一系列过程变量所组成的数据集合进行降维处理，进而获得相互之间互不相关的特征信号的过程，即用较少的主元信号表征过程数据矩阵的动态变化，换言之，主元分析将原始数据空间划分为两个子空间——主元特征信号子空间(Principal ComponentSubspace,PCS)和残差子空间(Residual Subspace,RS)，主元特征信号子空间包含了过程主要的、大部分的过程特征信息，而残差子空间刻画了过程中次要的、极小部分的过程动态信息。主元分析就是利用原始数据在这两个子空间的投影来表示数据发生的变化。

现有主元分析算法进行故障诊断存在不足，如出现误报现象。

发明内容

本发明是要解决现有主元分析算法进行故障诊断的不足，如出现误报现象的问题，而提供了一种基于巴特沃斯低通滤波器与主元分析（PCA）的故障检测和诊断方法。

基于主元分析算法的卫星陀螺部件故障诊断方法按以下步骤实现：

步骤一：基于巴特沃斯低通滤波器对卫星陀螺部件的输出角速度数据X_p进行滤波预处理；

步骤二：利用预处理后的卫星陀螺部件的输出角速度数据X，构建PCA数学模型；

步骤三：利用步骤二中得到的PCA数学模型参数，采用平方预测误差SPE统计量对残差子空间中过程数据进行检测；

步骤四：检测到故障后根据过程变量贡献图诊断出故障位置。

发明效果：

对陀螺输出角速度数据进行滤波，去除测量噪声的影响，避免故障误报现象的发生。适用于多个过程变量且相互耦合的过程建模。

1）本发明提出的基于巴特沃斯低通滤波器的主元分析方法可有效克服传统主元分析算法中的数据误报现象，提高故障检测的精度，在阶跃和斜坡故障中，精度分别提高了23.7%和10%。

2）本发明采用平方预测误差(SPE)统计量并结合过程变量贡献图，简化了检测过程，并能较好地完成故障诊断。

附图说明

图1为本发明的流程图；

图2为本发明实施例中的正常过程监测SPE统计量图；

图3为本发明实施例中的陀螺仪X轴发生阶跃故障时未对陀螺输出角速度数据滤波时的SPE统计量图；

图4为本发明实施例中的陀螺仪X轴发生阶跃故障时对陀螺输出角速度数据滤波之后的SPE统计量图；

图5为本发明实施例中的陀螺仪X轴发生斜坡故障时未对陀螺输出角速度数据滤波时的SPE统计量图；

图6为本发明实施例中的陀螺仪X轴发生斜坡故障时对陀螺输出角速度数据滤波之后的SPE统计量图；

图7为本发明实施例中的发生阶跃故障时的过程变量贡献图；

图8为本发明实施例中的斜坡故障时的过程变量贡献图。

具体实施方式

具体实施方式一：本实施方式的基于主元分析算法的卫星陀螺部件故障诊断方法按以下步骤实现：

步骤一：基于巴特沃斯低通滤波器对卫星陀螺部件的输出角速度数据Xp进行滤波预处理；

步骤二：利用预处理后的卫星陀螺部件的输出角速度数据X，构建PCA数学模型；

步骤三：利用步骤二中得到的PCA数学模型参数，采用平方预测误差SPE统计量对残差子空间中的过程数据进行检测；

步骤四：检测到故障后根据过程变量贡献图诊断出故障位置。

具体实施方式二：本实施方式与具体实施方式一不同的是：所述步骤一中基于巴特沃斯低通滤波器对卫星陀螺部件的输出角速度数据X_p进行滤波预处理具体为：

一、模型定义

速率陀螺仪的模型表示为：

w_out=w_in+D+N

其中，w_out为速率陀螺仪的测量输出，w_in为卫星实际角速度，D为随机漂移，N为测量噪声；

巴特沃斯低通滤波器幅度平方函数|H_a(jΩ)|²定义为

${\left|H_a\left(\mathrm{j\Ω}\right)\right|}^2=\frac{1}{1+{\left(\frac{\mathrm{\Ω}}{{\mathrm{\Ω}}_c}\right)}^{2N}}$

式中N为正整数代表滤波器的阶次，Ω表示模拟角频率，Ω_c称为截止频率，j表示虚数单位，a表示为模拟滤波器，H(·)表示模拟滤波器系统函数；

$\left|H_a\left(\mathrm{j\Ω}\right)\right|=\frac{1}{\sqrt{2}}$

${\mathrm{\δ}}_1=20\left|1g\frac{H_a\left(j0\right)}{H_a\left({\mathrm{j\Ω}}_c\right)}\right|=3\mathrm{dB}$

δ₁表示通带允许的最大衰减，所以又称Ω_c为巴特沃斯低通滤波器的3分贝带宽；

二、设计巴特沃斯低通滤波器及基于巴特沃斯低通滤波器进行滤波预处理：

1）将给定的低通数字滤波器的性能指标按照双线性变换规则转换成相应的模拟低通滤波器的性能指标；其中，所述性能指标包括通带截止频率w_p、阻带截止频率w_s、通带最大衰减α_p与阻带最小衰减α_s；

2）利用巴特沃斯低通滤波器逼近方法将所得到的模拟低通滤波器的性能指标进行设计并查表得到模拟低通滤波器的系统函数，以模拟低通滤波器的系统函数作为设计数字滤波器的“样本”；

3）利用双线性变换，将“样本”最终变换成所需的数字各型滤波器的系统函数其中，所述数字各型滤波器包括低通、高通、带通和带阻四个类型数字滤波器；

4）利用步骤3）中得到的b(z)，a(z)对卫星陀螺部件的输出角速度数据X_p进行滤波，即采用matlab中自带的函数X=filter(b(z),a(z),X_p)；其中，X为滤波之后的数据。

所述步骤二2）中的利用巴特沃斯低通滤波器逼近方法将所得到的模拟低通滤波器的性能指标进行设计并查表得到模拟低通滤波器的系统函数，以模拟低通滤波器的系统函数作为设计数字滤波器的“样本”如下：

巴特沃斯低通滤波器系统函数式中b(s)为分子多项式，a(s)为分母多项式，具体表示为：a(s)=s^N+a_N-1s^N-1+...+a₂s²+a₁s+1(a₀=a_N=1),s表示复变量，N表示巴特沃斯模拟低通滤波器的阶次,a_N,a_N-1,a_N-2...a₂,a₁,a₀表示分母多项式的系数

本实施方式中，假设测量噪声为白噪声，假设随机漂移的导数也是白噪声。测量噪声的水平将直接影响到最终主元模型的准确性，如出现故障误报现象，因而需要采用巴特沃斯低通滤波器对陀螺输出数据进行预处理。

其它步骤及参数与具体实施方式一相同。

具体实施方式三：本实施方式与具体实施方式一或二不同的是：

所述步骤二中利用预处理后的卫星陀螺部件的输出角速度数据X构建PCA数学模型具体为：

一、首先对滤波预处理后的角速度数据X采用以下方法进行标准化处理：

${X1}_{\mathrm{ij}}=\frac{X_{\mathrm{ij}}-\stackrel{\‾}{X}\left(j\right)}{s\left(j\right)},\mathrm{1,2},...n,j=\mathrm{1,2}...m$

其中，X∈R^n×m，R表示实数，n为陀螺仪输出数据即观测数据样本采样点个数，m为每一采样点过程变量个数，X(j)表示X的第j列，表示X(j)的均值，s(j)表示X(j)的标准差；

记标准化处理后的样本数据阵为X1；

二、求样本数据阵的协方差矩阵V

$V=\frac{1}{n-1}{X1}^{T}X1$

三、对协方差矩阵V进行特征值分解

V=PΛP^T

其中，Λ是协方差V的对角矩阵，其中Λ的主对角线元素为幅值递减的非负实特征值(λ₁≥λ₂≥...≥λ_m≥0)，P为负荷矩阵，且是标准正交矩阵，即：

p_ip_j=0(i≠j)

p_ip_j=1(i=j)

p_i是特征值λ_i对应的标准化特征向量；

四、采用主元分析法确定主元个数d

采用主元累计方差贡献率百分比法(CPV)来确定主元个数，在求主元时，数据的降维过程即为使数据在它的协方差阵所对应的特征向量上变化的过程；

在这引出方差贡献率的概念：

$a_i=\frac{{\mathrm{\λ}}_i}{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\λ}}_i}$

a_i为第i主元t_i的方差贡献率，λ_i为第i个主元ti的方差，称是前d个主元的累计方差贡献率，因而降维后的主元个数d就可以通过累计方差贡献率来确定，一般情况下选取85%或者90%以上即可满足要求，原来的m维空间就变成了d维空间，d≤m；

五、求得分矩阵T

选取负荷矩阵P∈R^m×d的各列，使其与前d个特征值的负荷向量相对应，那么X1到低维的投影就包含在了得分矩阵T中：

T=X1·P

最终得到的PCA模型如下： $X1={\mathrm{TP}}^T+E=\underset{i=1}{\overset{d}{\mathrm{\Σ}}}{t}_i{p}_i^T+E$

其中，E是X1的残差矩阵用来表示噪声，t、p分别代表矩阵T、P的列向量，T为主元模型得分矩阵，P为主元负荷矩阵，i取值从1～d。

本实施方式中，标准化的目的是使得样本点集合重心与坐标原点重合，压缩处理则可消除因为量纲不同而引起的虚假变化信息，从而使分析结果更为合理。

本实施方式中，步骤四中的主元分析法的理论推导

设X_p是一个有n个样本、m个变量组成的数据矩阵，X_p=(X_p(i,j）)_n×m=[X_p(1),X_p(2),...,X_p(m)],i=1,2,...n,j=1,2,...m

假设原始数据矩阵Xp经滤波预处理和标准化处理后记为X1（即均值E(X1(j))=0;方差var(X1(j))=1），现需求一个综合变量t₁，t₁是X1(1),...,X1(m)的线性组合，即t₁=X1·p₁，||p₁||=1；现使t₁能够携带尽可能多的原变异信息，也就是要求t₁的方差取最大值，t₁的方差为： $\mathrm{var}\left(t_1\right)=\frac{1}{n}{\left|\right|t_1\left|\right|}^2=\frac{1}{n}{p}_1^T\×{X1}^{T}X1p_1=p_1^T{V}_{p_1}$ 令 $V=\frac{1}{n}{X1}^{T}X1$ 是X1的协方差矩阵。

将以上问题转化为数学表达式的形式，即求

由拉格朗日算法，设λ₁是拉格朗日系数，令分别求L对λ₁和p1的偏导，并使其为零，即：

$\frac{\∂L}{{\∂\mathrm{\λ}}_1}=-(p_1^T{p}_1-1)=0$

$\frac{\∂L}{\∂p_1}=2{\mathrm{Vp}}_1-2{\mathrm{\λ}}_1{p}_1=0$

由此可得Vp₁=λ₁p₁，即p₁是V的一个标准化特征向量，对应的特征值为λ₁。

对于目标函数：即p₁所对应的特征值λ₁为t₁方差的最大值，换句话说，也就是p₁是矩阵V的最大特征值λ₁所对应的标准化特征向量，并把p₁称为主轴，t₁=Xp₁称为第一主成分。

与以上过程类似，可求出X1的第j主轴p_j，即协方差矩阵V的第j个特征值所对应的标准化特征向量，第j主成分t_j为：

t_j=X1·p_j

$\mathrm{var}\left(t_j\right)=\frac{1}{n}{p}_j^T{X1}^{T}X{1p}_j={\mathrm{\λ}}_j$

故有var(_t1)≥var(t₂)...≥var(t_d)，用数据变异大小来表示原始数据，则第一主元t₁携带的信息量最大，t₂次之，以此类推。

我们定义如下变量：

PCA的负荷矩阵为：

P_h=[p₁,p₂,...,p_m]=[P,P_r]；

P=[p₁,p₂,...,p_d]其构成的子空间称为主元子空间(PCS)，P∈R^m×d;

P_r=[p_d+1,p_d+2,...,p_m],其构成的子空间称为残差子空间(RCS)，P_r∈R^m×(m-d)；

PCA主成分得分矩阵为：

T_h=[T₁,T₂,...,T_m]=[T,T_r]；

T=[T₁,T₂,...,T_d]，前d个主元构成的得分矩阵，T∈R^n×d；

T_r=[T_d+1,T_d+2,...,T_m]，后m-d个主元构成的矩阵，T_r∈R^n×(m-d)。

则数据矩阵X1可表示为： $X1={\mathrm{TP}}^T+E=\underset{i=1}{\overset{d}{\mathrm{\Σ}}}{t}_i{p}_i^T+E$

令

称为X1的估计值，称为X1的残差，则X1可表示为：即主元分析将原始数据分解到两个子空间中。

其它步骤及参数与具体实施方式一或二相同。

具体实施方式四：本实施方式与具体实施方式一至三之一不同的是：所述步骤三中利用步骤二中得到的PCA数学模型参数，采用平方预测误差SPE统计量对残差子空间中过程数据进行检测具体为：

利用步骤二中得到的PCA模型参数，采用平方预测误差（SPE）统计量在残差子空间中检测运行的过程数据，采用SPE统计量的控制门限值判断运行过程中是否发生故障，当SPE值过大时，说明过程中出现了故障；

SPE统计量在任意时刻的值是一个标量，第i个样本点的SPE值为：

${\mathrm{SPE}}_i=e_i{e}_i^T=X1\left(i\right)(I-{\mathrm{PP}}^T)X1{\left(i\right)}^T$

式中，X1(i)是标准化之后第i个采样点数据，e_i是残差E的第i行，P=[p₁,p₂,...,p_d]，I是单位矩阵；SPE统计量表征了数据阵X1对主元模型的偏离程度；

当检验水平为α时，SPE统计量的控制门限值可按下式计算：

${\mathrm{SPE}}_a={\mathrm{\θ}}_1[\frac{h_o{c}_a\sqrt{{2\mathrm{\θ}}_2}}{{\mathrm{\θ}}_1}+1+\frac{{\mathrm{\θ}}_2{h}_0(h_0-1)}{{\mathrm{\θ}}_1^1}]1{/}_{h0}$

其中c_α是与上限(1-α)×100%相应的标准正态偏差，λ_j是与数据协方差第j个负荷向量相关的特征值，所述特征值按数值大小降序排列，d是主元个数，m是过程变量个数；

$h_0=1-\frac{{2\mathrm{\θ}}_1{\mathrm{\θ}}_3}{{3\mathrm{\θ}}_2^2},{\mathrm{\θ}}_l=\underset{j=d+1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\λ}}_j^l,l=\mathrm{1,2,3},$

如果SPE_i≤SPEα，说明时刻i即第i个检测点的SPE统计值正常；反之，若SPE_i>SPE_α，说明SPE统计值异常，即检测到有故障发生，称该时刻为故障时刻；

当SPE统计量超过其门限值，判断过程中出现了故障，因而检测出故障后需要进行故障定位。

本实施方式中，

SPE统计量是模型外部变化的一种测度，表征了数据阵X1对主元模型的偏离程度，

无论是否发生故障，均可以随着时间的推移，按照本步骤来对过程数据进行检测。

其它步骤及参数与具体实施方式一至三之一相同。

具体实施方式五：本实施方式与具体实施方式一至四之一不同的是：所述步骤四中检测到故障后根据过程变量贡献图诊断出故障位置，具体为：

设第i时刻检测到故障，变量X1_(i)为故障时刻i各变量的采样值，变量为故障时刻i各变量的估计值，则预测误差（即残差E的第i行）为：

第j个过程变量在第i时刻的SPE统计量的贡献值为：

其中e_i,j为故障时刻i时的预测误差ei的第j列，根据得到的故障时刻i，遍历j得到各个过程变量的贡献值，较大的贡献值所对应的过程变量即为引起SPE统计量超出其控制门限值的原因，由此诊断出故障发生的位置。

过程变量贡献值的大小代表了变量引起过程数据异常程度的大小，过程变量贡献值较大者多是引起故障的原因，由过程变量贡献值可明显观测出异常变量。

其它步骤及参数与具体实施方式一至四之一相同。

仿真实验：

执行步骤一：对某卫星系统仿真运行，得到陀螺输出角速度值作为样本数据XP，其中XP由陀螺仪四个轴的输出角速度组成，共有1800组采样值。

首先对样本数据进行基于巴特沃斯低通滤波器的滤波处理，巴特沃斯滤波器的最优参数为：

fp=10hz通带边界频率fs=105hz阻带截止频率

A_P=1dB通带最大衰减A_s=20dB阻带最小衰减

F_s=1000hz采样频率

基于以上参数设计的巴特沃斯低通滤波器如下：

$H\left(z\right)=\frac{0.0101(z^2+2z+1)}{z^2-1.6961z+0.7365}$

滤波器阶数：N=2；3dB截止频率：217.0827hz

利用上述设计的巴特沃斯低通滤波器对样本数据XP滤波得到新的样本数据，记为X。

执行步骤二：首先对样本数据X进行标准化处理工作，得到新的样本数据，记为X1，然后求解PCA模型，得到相应主元的贡献率。如表1所示：

表1主元贡献率及累积方差贡献率

变量序号	特征值	主元贡献率（%）	累积贡献率（%）
				1	2.8458	57.27	57.27
2	1.3895	27.96	85.24
				3	0.6076	12.23	97.47

4

0.1259

2.53

100

取累积方差贡献率为95%，则由表1可知，前三个主元的贡献率累积和百分比为97.47%，因而选择主元数d=3。

执行步骤三：依据样本模型，求出平方预测误差的控制门限值为SPE_alfa=0.2008，对于运行过程采样得到的数据，选取其中的700组采样值，分别计算出每组采样值的SPE统计量值，通过过程数据的SPE统计量值是否超过控制门限值来判断运行过程中是否发生故障。

对于正常运行阶段，得到1800个采样值，由图2可知，在正常运行阶段，SPE统计量值均处在控制门限值以下，说明运行过程正常，没有故障发生。

为了验证改进PCA(即先滤波再采用传统PCA算法)的有限性，在第300～320采样区间人为的对陀螺仪X轴注入故障等级为1的阶跃干扰，滤波前后SPE统计量的变化如图3和图4所示，由图3可知，在0～300之间有多点超出了控制限，即出现了误报现象，因此常规主元分析方法存在不足；而采用巴特沃斯低通滤波器对样本数据进行滤波处理后，由图4可知，在第300个采样时刻之前，SPE统计量都在控制门限值下面，在300～320采样时刻的SPE值超过控制门限值，即此时SPE统计图可准确检测出故障，且精度提高了23.7%。

为了进一步验证改进PCA的有效性，在第300～320采样区间人为的对陀螺仪X轴注入故障值为1.155的斜坡干扰，滤波前后SPE统计量的变化如图5和图6所示，由图5可知，在0～300之间有多点超出了控制限，即出现了误报现象，由图6可知，在第300个采样时刻之前，SPE统计量都在控制门限值下面，在300～320采样时刻的SPE值超过控制门限值，即此时SPE统计图可准确检测出故障，且精度提高了10%。

执行步骤四：根据过程变量贡献值计算公式得到阶跃故障时的过程变量贡献图如图7所示，由图7可知，第一个变量的贡献值最高，而这个变量为陀螺仪的X轴，由此可以判断必然是陀螺仪X轴出现异常，即X轴是引起故障的原因，与预期的故障原因一致。

同理，根据过程变量贡献率计算公式得到斜坡故障时的过程变量贡献图如图8所示，由图8可知，第一个变量的贡献值最高，而这个变量为陀螺仪的X轴，由此可以判断必然是陀螺仪X轴出现异常，即X轴是引起故障的原因，与预期的故障原因一致。

本发明详细说明了如何应用巴特沃斯低通滤波器和主元分析方法进行故障检测与诊断。综合实施例的上述分析，对于卫星陀螺部件的过程检测和故障诊断，本发明的算法能够有效克服传统主元分析模型误报的缺点，可以有效地对故障进行检测与诊断，并能有效提高故障检测的灵敏度。

Claims (4)

1.一种基于主元分析算法的卫星陀螺部件故障诊断方法，其特征在于基于主元分析算法的卫星陀螺部件故障诊断方法按以下步骤实现：

步骤一：基于巴特沃斯低通滤波器对卫星陀螺部件的输出角速度数据X_p进行滤波预处理；具体为：

一、模型定义

速率陀螺仪的模型表示为：

w_out＝w_in+D+N

其中，w_out为速率陀螺仪的测量输出，w_in为卫星实际角速度，D为随机漂移，N为测量噪声；

巴特沃斯低通滤波器幅度平方函数|H_a(jΩ)|²定义为

$|H_a\left(j\mathrm{\Ω}\right){|}^2=\frac{1}{1+{\left(\frac{\mathrm{\Ω}}{{\mathrm{\Ω}}_c}\right)}^{2N}}$

式中N为正整数代表滤波器的阶次，Ω表示模拟角频率，Ω_c称为截止频率，j表示虚数单位，a表示为模拟滤波器，H(·)表示模拟滤波器系统函数；

当Ω＝Ω_c时，有

$\left|H_a\left(j\mathrm{\Ω}\right)\right|=\frac{1}{\sqrt{2}}$

${\mathrm{\δ}}_1=20\left|\lg \frac{H_a\left(j0\right)}{H_a\left({\mathrm{j\Ω}}_c\right)}\right|=3dB$

δ₁表示通带允许的最大衰减，所以又称Ω_c为巴特沃斯低通滤波器的3分贝带宽；

二、设计巴特沃斯低通滤波器及基于巴特沃斯低通滤波器进行滤波预处理：

1)将给定的低通数字滤波器的性能指标按照双线性变换规则转换成相应的模拟低通滤波器的性能指标；其中，所述性能指标包括通带截止频率w_p、阻带截止频率w_s、通带最大衰减α_p与阻带最小衰减α_s；

2)利用巴特沃斯低通滤波器逼近方法将所得到的模拟低通滤波器的性能指标进行设计并查表得到模拟低通滤波器的系统函数，以模拟低通滤波器的系统函数作为设计数字滤波器的“样本”；

3)利用双线性变换，将“样本”最终变换成所需的数字各型滤波器的系统函数其中，所述数字各型滤波器包括低通、高通、带通和带阻四个类型数字滤波器；

利用步骤3)中得到的b(z)，a(z)对卫星陀螺部件的输出角速度数据X_p进行滤波，即采用matlab中自带的函数X＝filter(b(z),a(z),X_p)；其中，X为滤波之后的数据；

步骤二：利用预处理后的卫星陀螺部件的输出角速度数据X，构建PCA数学模型；

步骤三：利用步骤二中得到的PCA数学模型参数，采用平方预测误差SPE统计量对残差子空间中的过程数据进行检测；

步骤四：检测到故障后根据过程数据变量贡献图诊断出故障位置。

2.根据权利要求1所述的一种基于主元分析算法的卫星陀螺部件故障诊断方法，其特征在于所述步骤二中利用预处理后的卫星陀螺部件的输出角速度数据X，构建PCA数学模型具体为：

一、首先对滤波预处理后的角速度数据X，采用以下方法进行标准化处理：

$X{1}_{ij}=\frac{X_{ij}-\stackrel{\‾}{X}\left(j\right)}{s\left(j\right)},i=1,2,...n,j=1,2,...m$

其中，X∈R^n×m，R表示实数，n为陀螺仪输出数据即观测数据样本采样点个数，m为每一采样点过程变量个数，X(j)表示X的第j列，表示X(j)的均值，s(j)表示X(j)的标准差；

记标准化处理后的数据阵为X1；

二、求样本的协方差矩阵V

$V=\frac{1}{n-1}X{1}^{T}X1$

三、对协方差矩阵V进行特征值分解

V＝PΛP^T

其中，Λ是协方差V的对角矩阵，其中Λ的主对角线元素为幅值递减的非负实特征值，即λ₁≥λ₂≥…≥λ_m≥0，P为负荷矩阵，且是标准正交矩阵，即：

p_ip_j＝0(i≠j)

p_ip_j＝1(i＝j)

p_i是特征值λ_i对应的单位化特征向量；

四、采用主元分析法确定主元个数d

采用主元累计方差贡献率百分比法CPV，确定主元个数，在求主元时，数据的降维过程即为使数据在它的协方差阵所对应的特征向量上变化的过程；

在这引出方差贡献率的概念：

$a_i=\frac{{\mathrm{\λ}}_i}{\underset{i=1}{\overset{m}{\Σ}}{\mathrm{\λ}}_i}$

a_i为第i主元t_i的方差贡献率，λ_i为第i个主元t_i的方差，称是前d个主元的累计方差贡献率，因而降维后的主元个数d就可以通过累计方差贡献率来确定，原来的m维空间变成了d维空间，d≤m；

五、求得分矩阵T

选取负荷矩阵P∈R^m×d的列，使其与前d个特征值的负荷向量相对应，那么X1到低维的投影就包含在了得分矩阵T中：

T＝X1·P

最终得到的PCA模型如下：

其中，E是X1的残差矩阵用来表示噪声，t、p分别代表矩阵T、P的列向量，i取值从1～d。

3.根据权利要求1所述的一种基于主元分析算法的卫星陀螺部件故障诊断方法，其特征在于所述步骤三中利用步骤二中得到的PCA数学模型参数，采用平方预测误差SPE统计量对残差子空间中检测过程数据进行检测具体为：

利用步骤二中得到的PCA模型参数，采用平方预测误差SPE统计量在残差子空间中检测运行的过程数据，采用SPE统计量的控制界限值判断运行的过程是否发生故障，当SPE值过大时，说明过程中出现了故障并及时检测；

SPE统计量在任意时刻的值是一个标量，第i个样本点的SPE值为：

${\mathrm{SPE}}_i=e_i{e}_i^T=X1\left(i\right)(I-{\mathrm{PP}}^T)X1{\left(i\right)}^T$

式中，X1(i)是标准化之后第i个采样点数据，e_i是残差E的第i行，P＝[p₁,p₂,...,p_d]，I是单位矩阵；SPE统计量表征了数据阵X1对主元模型的偏离程度；

当检验水平为α时，SPE统计量的控制门限值可按下式计算：

${\mathrm{SPE}}_{\mathrm{\α}}={\mathrm{\θ}}_1{\[\frac{h_o{c}_{\mathrm{\α}}\sqrt{2{\mathrm{\θ}}_2}}{{\mathrm{\θ}}_1}+1+\frac{{\mathrm{\θ}}_2{h}_0(h_0-1)}{{\mathrm{\θ}}_1^2}\]}^{1/h_0}$

其中c_α是与上限(1-α)×100％相应的标准正态偏差，λ_j是与数据协方差第j个负荷向量相关的特征值，所述特征值按数值大小降序排列，d是主元个数，m是过程变量个数；

$\begin{array}{cc}{h}_0=1-\frac{2{\mathrm{\θ}}_1{\mathrm{\θ}}_3}{3{\mathrm{\θ}}_2^2}& {\mathrm{\θ}}_l=\underset{j=d+1}{\overset{m}{\Σ}}{\mathrm{\λ}}_j^l,l=1,2,3,\end{array}$

如果SPE_i≤SPE_α，说明时刻i即第i个检测点的SPE统计值正常；反之，若SPE_i>SPE_α，说明SPE统计值异常，即检测到有故障发生，称该时刻为故障时刻；

当SPE统计量超过其门限值，判断过程中出现了故障，因而检测出故障后需要进行故障定位。

4.根据权利要求1所述的一种基于主元分析算法的卫星陀螺部件故障诊断方法，其特征在于所述步骤四中检测到故障后根据过程数据变量贡献图诊断出故障位置，具体为：

设第i时刻检测到故障，变量X1_(i)为故障时刻i各变量的采样值，变量为故障时刻i各变量的估计值，则预测误差即残差E的第i行为：

第j个过程变量在第i时刻的SPE统计量的贡献值为：其中e_i,j为故障时刻i时的预测误差e_i的第j列，根据得到的故障时刻i，遍历j得到各个过程变量的贡献值，较大的贡献值所对应的过程变量即为引起SPE统计量超出其控制门限值的原因，由此诊断出故障发生的位置。