CN104019831A - 基于emd和熵权的陀螺仪故障诊断方法 - Google Patents

Abstract

基于EMD和熵权的陀螺仪故障诊断方法，涉及一种陀螺仪的故障诊断方法，具体涉及一种基于经验模态分解和熵权的故障诊断方法，属于陀螺仪故障诊断领域。为解决针对现有的陀螺仪故障诊断方法的实用性和有效性受到限制，陀螺仪运行情况监测不足，单一信号故障检测的自适应性弱，检测过程繁琐的问题。本发明引入了一种经典的故障诊断方法—经验模态分解(EMD)方法，这是一种单一信号时域处理方法，即每次只能对一个轴上的陀螺角速度信号提取故障特征信息，然后创新性的使用熵权概念来实现故障诊断，这样可以更好的实现故障诊断，最后对所提出方法的有效性进行了仿真验证。本发明适用于卫星陀螺仪故障的诊断。

Description

基于EMD和熵权的陀螺仪故障诊断方法

技术领域

本发明涉及一种陀螺仪的故障诊断方法，具体涉及一种基于经验模态分解(EmpiricalMode Decomposition,EMD)和熵权的故障诊断方法，属于陀螺仪故障诊断领域。

背景技术

卫星的姿态是指卫星星体在轨道上运行时在所处的空间位置的状态，它一般相对于某一参考坐标系而言。在卫星应用任务中，通常要求卫星姿态在工作空间中应该保持高精度定向，或者能够按照既定要求在空间中完成相应的姿控。在卫星姿控过程中，卫星姿态运动学与卫星姿态动力学起着极为重要的作用，它是卫星姿态学中的重要理论基础，由它可以对被控对象进行数学建模，从而根据既定性能要求来设计卫星姿态控制系统，这也为卫星姿态控制系统的故障诊断提供了基础。因而卫星姿控系统最重要的任务就是实现姿态确定，进而实现姿态准确控制。

在姿控系统中，将陀螺仪与星敏感器两种传感器结合使用的方式在高精度姿态确定系统中经常使用，卫星整体的可靠性和姿态控制精度很大程度上受到陀螺仪的影响，因而对其进行故障诊断意义重大。对陀螺仪运行情况的监测主要有两种方式：一是采用某些机内检测设备来对陀螺仪的某些状态变量进行监测，例如电压、电流、温度等，但这仅能够检测出变量幅度变化较大的故障。二是利用具有冗余关系的姿态敏感器来实现故障诊断，这种方式也具有局限性，如果姿态敏感器之间不具有冗余关系，或者具有冗余关系但冗余部件不能正常工作，则这种方式不能使用。由于这两种方式各有利弊，因而实用性和有效性受到限制。

发明内容

本发明的目的是提出一种基于EMD和熵权的陀螺仪故障诊断方法，以解决针对现有的陀螺仪故障诊断方法的实用性和有效性受到限制，陀螺仪运行情况监测不足，单一信号故障检测的自适应性弱，检测过程繁琐的问题。

本发明为解决上述技术问题所采用的技术方案是：

本发明所述的基于EMD和熵权的陀螺仪故障诊断方法，是按照以下步骤实现的：

步骤一：对原始陀螺角速度输出数据X_p采用滑动窗的形式选取10s，基于巴特沃斯低通滤波器对数据进行滤波预处理后作为故障诊断用信号，记为X；

步骤二：利用EMD方法对步骤一中得到的故障诊断用信号X进行分解，获得各次IMF分量以及残差分量RES；

步骤三：选取步骤二中得到的一次IMF分量作为熵权计算的特征量，由能量熵公式计算得到选取的一次IMF分量的能量熵，进而求得熵权变化值，根据熵权变化阈值来对陀螺仪进行故障诊断。

本发明的有益效果是：

一、对陀螺部件的角速度输出信号，在原始经验模态分解算法的基础上，引入能量熵权的概念，可对阶跃、缓变、卡死等故障进行准确检测与识别。适用于单一信号的过程监测，且对于非平稳信号的故障诊断具有较强的优势。

二、本发明可有效克服现有陀螺仪运行情况监测方法的不足，提高单一信号故障检测的自适应性，较好的克服了陀螺仪故障信号的非平稳性；

三、本发明采用能量熵权的概念并结合熵权变化曲线图，对陀螺仪的运行情况进行监测，可以保证在0.5s内检测出发生的故障，检测过程与现有方法相比提高了0.2～0.5s,能够较好地提高故障诊断能力。

四、本发明方法的实用性和有效性与现有方法相比适用范围更广。

附图说明

图1为本发明的流程图；

图2为本发明的EMD算法的流程图；

图3左右两图分别为本发明的滤波前后陀螺角速度输出信号对比图；

图4为本发明的正常陀螺角速度输出信号EMD分解结果图；

图5为本发明的发生故障时陀螺角速度输出信号EMD分解结果图；

图6为本发明的正常陀螺角速度输出信号一次IMF分量熵权变化曲线图；

图7为本发明的发生单阶跃故障时EMD分解结果图；

图8为本发明的发生单阶跃故障时一次IMF熵权变化曲线图；

图9为本发明的发生多阶跃故障时EMD分解结果图；

图10为本发明的发生多阶跃故障时一次IMF熵权变化曲线图；

图11为本发明的发生卡死故障时EMD分解结果图；

图12为本发明的发生卡死故障时一次IMF熵权变化曲线图；

图13为本发明的发生缓变故障时EMD分解结果图；

图14为本发明的发生缓变故障时一次IMF熵权变化曲线图；

具体实施方式

具体实施方式一：本实施方式所述的基于EMD和熵权的陀螺仪故障诊断方法，其特征在于所述方法是按以下步骤实现的：

步骤一：对原始陀螺角速度输出数据X_p采用滑动窗的形式选取10s，基于巴特沃斯低通滤波器对数据进行滤波预处理后作为故障诊断用信号，记为X；

步骤二：利用EMD方法对步骤一中得到的故障诊断用信号X进行分解，获得各次IMF分量以及残差分量RES；

步骤三：选取步骤二中得到的一次IMF分量作为熵权计算的特征量，由能量熵公式计算得到选取的一次IMF分量的能量熵，进而求得熵权变化值，根据熵权变化阈值来对陀螺仪进行故障诊断。

结合图1解释本实施方式。当陀螺发生故障时，陀螺角速度输出信号的熵权值将会出现较大的波动性。正常陀螺角速度输出信号EMD分解曲线与发生故障时陀螺角速度输出信号EMD分解结果曲线分别如图4和图5所示。由图中可以看到，当陀螺在7.8s时刻发生故障时，其前几个IMF分量的幅值甚至是频率通常会发生较大的变化，因而各个IMF分量的熵权值也必然会发生较大的波动，因而可以将原始信号经过EMD分解获得的IMF分量作为特征量计算熵权来实现对陀螺的故障诊断。

EMD算法的基本理论如下：

1998年，N.E.Huang等人提出了EMD算法，该方法是对信号局部特征进行分解的方法，它通过分析信号能量与频率随时间的变化，将信号分解为多个固有模态函数(IntrinsicMode Function,IMF)，本质上是对信号进行了平稳化处理。

EMD算法与传统方法不同，它的分解基底函数不需要使用固定形态窗口，而是由所处理的信号提取得到，即以IMF作为基底，因而信号不同，基底就有所不同。其中IMF必须满足条件：

(1)在整个信号范围上的极值点和过零点数量相等，或者只相差一个。

(2)在任意时刻，由极大值所形成的上包络线和由极小值所形成的下包络线的包络均值为零。

EMD算法的整体思路是利用信号的上下包络的均值来确定“瞬时平衡位置”，进而提取固有模态函数，并且方法基于如下假设：

(1)信号至少具有两个极值：极大值和极小值；

(2)信号的特征时间尺度是由极值点的时间间隔确定；

(3)如果信号没有极值而仅有拐点，可以通过微分、分解、再积分的方法获得IMF。

在此假设基础上，Huang等人进一步指出，可以用EMD算法来筛选出信号的固有模态函数(IMF)，即EMD算法的核心是从信号中提取IMF分量的过程，称为“筛选过程”。给定信号x(t)∈R¹，其筛选过程如下：

(1)找出原始信号x(t)时间序列上的所有的局部极值点；

(2)分别由极大值和极小值构造生成x(t)的上包络函数和下包络函数，分别记为e_max(t)和e_min(t)；

(3)生成上、下包络函数的均值函数：

$m_1\left(t\right)=\frac{e_{\max }\left(t\right)+e_{\min }\left(t\right)}{2}---\left(1\right)$

(4)求得信号x(t)与包络均值函数m₁(t)的差值函数：

h₁(t)＝x(t)-m₁(t)   (2)

理想情况下，若h₁(t)满足IMF的两个条件，则h₁(t)为x(t)的第一个IMF分量；

(5)通常情况下，h₁(t)不满足IMF条件，此时将h₁(t)作为原始信号，重复(1)-(3)步骤可以得到该时的包络线均值函数m₁₁(t)，继而求得此时的差值函数h₁₁(t)＝h₁(t)-m₁₁(t)，继续判断h₁₁(t)是否满足IMF的条件，如果不满足则继续，循环到第k次时，得到的差值函数为h_1k(t)＝h_1(k-1)(t)-m_1k(t)，此时h_1k(t)满足IMF条件，第一个IMF分量：c₁(t)＝h_1k(t)。

为了使IMF分量能够尽可能多的反映原始信号的局部信息而又不进行额外的迭代循环，必须确定一个迭代停止准则。这里引入了标准差(SD)准则，SD的计算公式为：

$\mathrm{SD}=\frac{\underset{t=0}{\overset{T}{\mathrm{\Σ}}}{|h_{1(k-1)}\left(t\right)-h_{1k}\left(t\right)|}^2}{\underset{t=0}{\overset{T}{\mathrm{\Σ}}}{h}_{1(k-1)}^2\left(t\right)}---\left(3\right)$

SD的阈值一般取0.2～0.3，标准差值小于既定阈值，则筛选过程结束。

分量c₁(t)代表的是原信号x(t)中的最高频成分，因而残差函数r₁(t)中包含了低频成分：

r₁(t)＝x(t)-c₁(t)   (4)

对r₁(t)循环进行筛选，获得第二个IMF分量c₂(t)。如此进行n次，可以得到n个IMF分量，以及一个残差函数r_n(t)。由此，x(t)组成为：

$x\left(t\right)=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{c}_i\left(t\right)+r_n\left(t\right)---\left(5\right)$

其中r_n(t)代表了原始信号x(t)的平均趋势，而c₁(t),c₂(t),…,c_n(t)分别代表了x(t)从高到低不同频率段的成分，这些成分都是不同的，EMD算法的流程图如图2所示。

由于EMD算法法是一种自适应的时域信号处理方法，它将原始信号中的IMF分量作为基底，避免了因为基函数选取不同，而对信号分解产生影响。另一方面，前几个IMF分量包含原始信号中最重要、最显著的信息。鉴于此，本发明采用了EMD算法来对陀螺仪角速度信号进行故障特征分析。

具体实施方式二：本实施方式与具体实施方式一不同的是：步骤一的具体过程为：

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

步骤一二、利用巴特沃斯低通滤波器逼近方法将所得到的模拟低通滤波器的性能指标计算得到滤波器阶数，并通过表1，巴特沃斯低通滤波器系统函数一览表得到模拟低通滤波器的系统函数，以模拟低通滤波器的系统函数作为设计数字滤波器的样本；

步骤一三、利用双线性变换，将样本最终变换成所需的数字各型滤波器的系统函数b(z)为数字各型滤波器的系统函数的分子多项式，a(z)为数字各型滤波器的系统函数的分母多项式；

步骤一四、利用步骤一三中得到的b(z)，a(z)对卫星陀螺部件的角速度输出数据X_p进行滤波。

表1巴特沃斯低通滤波器系统函数一览表

结合图3解释本实施方式。其它步骤与具体实施方式一相同。

具体实施方式三：本实施方式与具体实施方式一或二不同的是：步骤一二所述的样本如下：

巴特沃斯低通滤波器系统函数式中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₀表示分母多项式的系数。其它步骤与具体实施方式一至三之一相同。

具体实施方式四：本实施方式与具体实施方式一至三之一不同的是：步骤一三所述数字各型滤波器包括：低通、高通、带通或带阻四个类型数字滤波器。其它步骤与具体实施方式一至三之一相同。

具体实施方式五：本实施方式与具体实施方式一至四之一不同的是：步骤一四所述的滤波方式为：采用matlab中自带的函数X＝filter(b(z),a(z),X_p)进行滤波；其中，X为滤波之后的数据。其它步骤与具体实施方式一至四之一相同。

具体实施方式六：本实施方式与具体实施方式一至五之一不同的是：步骤二的具体过程为：

步骤二一、找出步骤一中得到的数据X时间序列上的所有的局部极值点；

步骤二二、分别由极大值和极小值构造生成X的上包络函数和下包络函数，分别记为e_max(t)和e_min(t)；

步骤二三、生成上、下包络函数的均值函数：基于公式

$m_1\left(t\right)=\frac{e_{\max }\left(t\right)+e_{\min }\left(t\right)}{2}---\left(1\right)$

得到上、下包络函数的均值函数m₁(t)；

步骤二四、基于公式

h₁(t)＝x(t)-m₁(t)   (2)

求得信号X与包络均值函数m₁(t)的差值函数：

h₁(t)＝X-m₁(t)   (6)

步骤二五、当h₁(t)不满足IMF条件，此时将h₁(t)作为原始信号，重复(1)-(3)步骤得到该时的包络线均值函数m₁₁(t)，继而求得此时的差值函数h₁₁(t)＝h₁(t)-m₁₁(t)，继续判断h₁₁(t)是否满足IMF的条件，如果不满足则继续，循环到第k次时，得到的差值函数为h_1k(t)＝h_1(k-1)(t)-m_1k(t)，此时h_1k(t)满足IMF条件，一次IMF分量：c₁(t)＝h_1k(t)；

根据式

$\mathrm{SD}=\frac{\underset{t=0}{\overset{T}{\mathrm{\Σ}}}{|h_{1(k-1)}\left(t\right)-h_{1k}\left(t\right)|}^2}{\underset{t=0}{\overset{T}{\mathrm{\Σ}}}{h}_{1(k-1)}^2\left(t\right)}---\left(3\right)$

计算各次IMF标准差值(SD)，并与既定阈值比较，若标准差值小于既定阈值，则筛选过程结束；

分量c₁(t)代表的是原信号X中的最高频成分，因而残差函数r₁(t)中包含了低频成分：

r₁(t)＝X-c₁(t)   (7)

对r₁(t)循环进行筛选，获得二次IMF分量c₂(t)；如此进行n次，可以得到n个IMF分量，以及一个残差函数r_n(t)；由此，X组成为：

$X=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{c}_i\left(t\right)+r_n\left(t\right)---\left(8\right)$

其中r_n(t)代表了原始信号X的平均趋势，而c₁(t),c₂(t),…，c_n(t)分别代表了X从高到低不同频率段的成分。这些成分都是不同的。其它步骤与具体实施方式一至五之一相同。

具体实施方式七：本实施方式与具体实施方式一至六之一不同的是：步骤三的具体过程为：

步骤三一、计算步骤二中得到的一次IMF分量的熵权变化曲线图；首先将步骤二中得到的一次IMF分量作为特征量，对其分段，将数据分成10段，即m＝10，每段数据长度n＝40；然后基于公式

$H_j=-\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{p}_i\mathrm{lo}{g}_a{p}_i,j=\mathrm{1,2},...m---\left(10\right)$

其中，H_j为第j个数据段的能量熵，p_i＝E_i/E，E_i是一次IMF信号的瞬时能量，而E为一次IMF信号总能量，n为每段数据长度，m为原始数据的分段总数；

与公式

$q_j=\frac{H_j}{\underset{j=1}{\overset{m}{\mathrm{\Σ}}}{H}_j},j=\mathrm{1,2},\·\·\·,m---\left(11\right)$

计算能量熵与熵权值，并将各点的熵权值连接形成熵权变化曲线图；

步骤三二、计算熵权变化的阈值；熵权变化的阈值是根据正常无故障状态下的陀螺角速度信号的熵权波动性进行设定的；由于标准差可以描述变量的波动情况，因而将由陀螺角速度正常信号计算得到的熵权值的标准差作为熵权变化的阈值；

步骤三三、根据步骤三一与步骤三二得到的故障诊断用信号的熵权变化曲线图和熵权变化阈值，进行故障诊断；

若熵权变化曲线图中的某个时间点或者时间段的熵权值超过熵权变化阈值，则可判断该时间点或者时间段发生故障。其它步骤与具体实施方式一至六之一相同。

熵和熵权的基本理论如下：

熵是对体系混乱程度的一种度量，开始是由鲁道夫.克劳修斯(RudolfClausius)提出，并应用在热力学中，后在1948年，由克劳德.埃尔伍德.香农(ClaudeElwoodShannon)第一次将熵的概念引入到信息论中，用熵表示不确定性，并且在著作《AMathematicalTheoryofCommunication》中提出了建立在概率统计模型上的信息度量，首次使用信息熵来表示事物中含有的信息量，实现了对参量特征的定量描述。

假设有一个离散随机变量Xs，Xs＝{x₁,x₂,…x_n}，则它的信息熵定义：

$H-\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{p}_i\mathrm{lo}{g}_a{p}_i---\left(9\right)$

其中，p_i＝P(x_i)为x_i出现的概率，并且满足

信息熵是从宏观上来表示随机变量Xs的信息量，同时又是对随机变量Xs有序性和不确定性的一种度量，信息熵越大，信息量越多，Xs的不确定性越大。

类似于信息熵，将一次IMF分量数据分成若干段，每段数据长度相同，则每个分段数据的能量熵定义为：

$H_j=-\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{p}_i\mathrm{lo}{g}_a{p}_i,j=\mathrm{1,2},...m---\left(10\right)$

其中，H_j为第j个数据段的能量熵，p_i＝E_i/E，E_i是一次IMF信号的瞬时能量，而E为一次IMF信号总能量，n为每段数据长度，m为原始数据的分段总数。

熵权，即为各瞬时熵在系统总熵中所占的权重，可以描述每一时刻信号的变异程度。

每个分段数据的熵权值q_j定义为

$q_j=\frac{H_j}{\underset{j=1}{\overset{m}{\mathrm{\Σ}}}{H}_j},j=\mathrm{1,2},\·\·\·,m---\left(11\right)$

从能量熵的公式可以看到：

(1)如果某一时刻的能量熵值越小，说明此时刻信号的变异程度越大，所提供的能量越多，在整体上所起的影响越大，其熵权值越大；

(2)如果某一时刻的瞬时能量熵值越大，说明此时刻信号的变异程度越小，所提供的能量越少，在整体上所起的影响越小，其熵权值越小；

(3)在具体应用时，可根据各个时刻的变异程度，由瞬时熵来计算瞬时的熵权值。

仿真实验：

步骤一：卫星姿态控制系统各组成部分的模型在Matlab平台上Simulink环境下搭建闭环仿真系统，其中仿真步长为0.025s，对仿真系统运行，得到陀螺角速度输出信号X_p，并采用滑动窗的形式选取400个共10s长的陀螺角速度输出数据，仍记为X_p。

首先对陀螺角速度输出数据X_p采用巴特沃斯低通滤波器预处理，巴特沃斯滤波器的最优参数为：

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

利用上述设计的巴特沃斯低通滤波器对陀螺角速度输出数据X_p滤波得到新的样本数据，记为X。滤波前后陀螺角速度输出数据的对比图如图3所示。

步骤二：利用EMD方法对步骤一中得到的滤波预处理后的故障诊断用信号X进行分解，得到多个不同频段IMF分量以及残差分量RES。如图4所示。

步骤三：基于步骤二中得到的各次IMF分量，选取一次IMF分量作为熵权计算的特征量，由能量熵公式可以计算得到选取的一次IMF分量的能量熵，进而求得熵权变化值，各点的熵权变化值如表2所示：

表2正常信号的分段数据熵权值表

将各点的熵权值连接形成曲线，如图6所示。然后根据既定的熵权变化阈值来对陀螺进行故障诊断。

经过计算，可以得到图6中正常情况时陀螺角速度输出信号熵权的标准差为0.2024，这表示了正常状态下熵权值的一个波动大小，由熵权原理可以得知，当过程中发生故障时，则此刻的熵权波动性应该大于0.2024，从而将熵权值0.2024作为故障诊断时的熵权阈值。

为了验证本发明所提出方法的有效性，这里以偏航轴上陀螺角速度输出信号为例，对其中几种常见故障，如阶跃故障、缓变故障、卡死故障，进行故障诊断验证，具体的应用过程如下。

1)情况1：发生阶跃故障

当卫星姿态控制系统的偏航轴陀螺出现发生阶跃故障情况时，又可具体分为两种情况，一种是在过程中仅仅发生了一次故障，称为单故障情况；另一种是在过程中频繁的出现了多次单故障情况，这称为多故障情况。

(1)发生单故障情况

此种情况下，卫星姿态控制系统中的偏航轴陀螺在t＝1s时刻发生了阶跃故障，故障大小为2×10^-4rad/s，约为50倍噪声均方差值。

故障诊断用信号选取了400个共10s长的陀螺角速度数据，经过EMD方法分解得到其一次IMF分量，如图7所示，然后对分量进行分段，每隔40个数据共1s时间长度计算其能量熵和熵权值，然后作出其整个过程的熵权变化曲线，各分段数据的熵权值如表3所示，熵权变化曲线如图8所示。

表3单阶跃故障的分段数据熵权值表

结果分析：在图7中可以看到发生单阶跃故障后陀螺角速度输出信号的一次IMF分量发生了同样的阶跃变化，但是效果并不是非常明显，对这个过程中的熵权变化曲线图8中可以看到，在区间1s～2s这个时间段上对应的熵权值明显超过了阈值，从而判断出该轴发生了故障。

(2)发生多故障情况

在此情况下，以过程中发生两次阶跃故障为例对方法进行验证，卫星姿态控制系统中的偏航轴陀螺在t＝1s，t＝5s时刻发生了阶跃故障，故障值为2×10^-4rad/s，约为50倍噪声均方差值。

这里故障诊断用信号同样也选取了400个共10s长的陀螺角速度数据，经过EMD方法对信号进行分解，得到的结果如图9所示，对其中的一次IMF分量进行分段处理，每隔40个数据共1s时间长度计算其能量熵和熵权值，然后作出其整个过程的熵权变化曲线，各分段数据的熵权值如表4所示，熵权变化曲线图如图10所示。

表4多阶跃故障的分段数据熵权值表

结果分析，在图9中可以看到发生单阶跃故障后陀螺角速度输出信号的一次IMF分量发生了同样的阶跃变化但是效果并不是十分明显，对这个过程中的熵权变化曲线图10中可以看到，在区间1s～2s和5s～6s这两个时间段上对应的熵权值明显超过了阈值，从而判断出该轴发生了多次故障。

2)情况2：发生卡死故障

当陀螺在运行过程中因某些原因，如堵转、断电等因素，会导致此刻陀螺角速度输出信号为零，发生这种情况时，使卫星姿态控制系统中偏航轴上陀螺在t＝7.8s时刻输出信号为零。

同样的，这里故障诊断用信号也选取了400个共10s长的陀螺角速度数据，经过EMD方法对信号进行分解，得到的结果如图11所示，对其中的一次IMF分量进行分段处理，每隔40个数据共1s时间长度计算其能量熵和熵权值，然后作出其整个过程的熵权变化曲线，各分段数据的熵权值如表5所示，熵权变化曲线图如图12所示。

表5卡死故障的分段数据熵权值表

结果分析：在图11中可以看到发生卡死故障后陀螺角速度输出信号的一次IMF分量发生了突变，但是效果并不是十分明显，从这个过程中的熵权变化曲线图12中可以看到，在区间7s～8s时间段上对应的熵权值明显超过了阈值，可以检测出过程中发生的故障。

3)情况3：发生缓变故障

有时候由于陀螺上某些部件的原因可能会导致陀螺在运行过程中发生缓变漂移故障，这种情况下，卫星姿态控制系统偏航轴上的陀螺在t＝8s时刻发生缓变漂移故障，其斜率为2.5×10^-5rad/s²。

对该陀螺角速度信号作EMD分解，得到的分解结果如图13所示。然后对其进行分段处理，每40个数据作为一段分别计算能量熵与熵权，得到在这个过程中的熵权变化曲线，各分段数据的熵权值如表6所示，熵权变化曲线图如图14所示。

表6缓变故障的分段数据熵权值表

结果分析：在图13中可以看到发生缓变故障后陀螺角速度输出信号的一次IMF分量发生了不十分明显的变化，从这个过程中熵权变化曲线图14中可以看到，在区间8s～10s时间段上对应的熵权值明显超过了阈值，说明方法能够可以检测出过程中发生的故障。

本发明详细说明了如何应用经验模态分解算法与熵权进行陀螺仪故障诊断。综合实施例的上述分析，对于陀螺仪的故障诊断，本发明的算法能够有效克服现有陀螺仪运行情况监测方法的不足，可以有效地对阶跃、缓变、卡死等故障进行诊断。

Claims (7)

1.基于EMD和熵权的陀螺仪故障诊断方法，其特征在于所述方法是按以下步骤实现的：

步骤一：对原始陀螺角速度输出数据X_p采用滑动窗的形式选取10s，基于巴特沃斯低通滤波器对数据进行滤波预处理后作为故障诊断用信号，记为X；

步骤二：利用EMD方法对步骤一中得到的故障诊断用信号X进行分解，获得各次IMF分量以及残差分量RES；

步骤三：选取步骤二中得到的一次IMF分量作为熵权计算的特征量，由能量熵公式计算得到选取的一次IMF分量的能量熵，进而求得熵权变化值，根据熵权变化阈值来对陀螺仪进行故障诊断。

2.根据权利要求1所述的基于EMD和熵权的陀螺仪故障诊断方法，其特征在于步骤一的具体过程为：

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

步骤一二、利用巴特沃斯低通滤波器逼近方法将所得到的模拟低通滤波器的性能指标计算得到滤波器阶数，并通过巴特沃斯低通滤波器系统函数一览表得到模拟低通滤波器的系统函数，以模拟低通滤波器的系统函数作为设计数字滤波器的样本；

步骤一三、利用双线性变换，将样本最终变换成所需的数字各型滤波器的系统函数b(z)为数字各型滤波器的系统函数的分子多项式，a(z)为数字各型滤波器的系统函数的分母多项式；

步骤一四、利用步骤一三中得到的b(z)，a(z)对卫星陀螺部件的角速度输出数据X_p进行滤波。

3.根据权利要求2所述的基于EMD和熵权的陀螺仪故障诊断方法，其特征在于步骤一二所述的样本如下：

巴特沃斯低通滤波器系统函数式中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₀表示分母多项式的系数。

4.根据权利要求3所述的基于EMD和熵权的陀螺仪故障诊断方法，其特征在于步骤一三所述数字各型滤波器包括：低通、高通、带通或带阻四个类型数字滤波器。

5.根据权利要求4所述的基于EMD和熵权的陀螺仪故障诊断方法，其特征在于步骤一四所述的滤波方式为：采用matlab中自带的函数X＝filter(b(z),a(z),X_p)进行滤波；其中，X为滤波之后的数据。

6.根据权利要求5所述的基于EMD和熵权的陀螺仪故障诊断方法，其特征在于步骤二的具体过程为：

步骤二一、找出步骤一中得到的数据X时间序列上的所有的局部极值点；

步骤二二、分别由极大值和极小值构造生成X的上包络函数和下包络函数，分别记为e_max(t)和e_min(t)；

步骤二三、生成上、下包络函数的均值函数：基于公式

$m_1\left(t\right)=\frac{e_{\max }\left(t\right)+e_{\min }\left(t\right)}{2}---\left(1\right)$

得到上、下包络函数的均值函数m₁(t)；

步骤二四、基于公式

h₁(t)＝x(t)-m₁(t)   (2)

求得信号X与包络均值函数m₁(t)的差值函数：

h₁(t)＝X-m₁(t)   (6)

步骤二五、当h₁(t)不满足IMF条件，此时将h₁(t)作为原始信号，重复(1)-(3)步骤得到该时的包络线均值函数m₁₁(t)，继而求得此时的差值函数h₁₁(t)＝h₁(t)-m₁₁(t)，继续判断h₁₁(t)是否满足IMF的条件，如果不满足则继续，循环到第k次时，得到的差值函数为h_1k(t)＝h_1(k-1)(t)-m_1k(t)，此时h_1k(t)满足IMF条件，一次IMF分量：c₁(t)＝h_1k(t)；

根据式

$\mathrm{SD}=\frac{\underset{t=0}{\overset{T}{\mathrm{\Σ}}}{|h_{1(k-1)}\left(t\right)-h_{1k}\left(t\right)|}^2}{\underset{t=0}{\overset{T}{\mathrm{\Σ}}}{h}_{1(k-1)}^2\left(t\right)}---\left(3\right)$

计算各次IMF标准差值(SD)，并与既定阈值比较，若标准差值小于既定阈值，则筛选过程结束；

分量c₁(t)代表的是原信号X中的最高频成分，因而残差函数r₁(t)中包含了低频成分：

r₁(t)＝X-c₁(t)   (7)

对r₁(t)循环进行筛选，获得二次IMF分量c₂(t)；如此进行n次，可以得到n个IMF分量，以及一个残差函数r_n(t)；由此，X组成为：

$X=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{c}_i\left(t\right)+r_n\left(t\right)---\left(8\right)$

其中r_n(t)代表了原始信号X的平均趋势，而c₁(t),c₂(t),…,c_n(t)分别代表了X从高到低不同频率段的成分。

7.根据权利要求6所述的基于EMD和熵权的陀螺仪故障诊断方法，其特征在于步骤三的具体过程为：

步骤三一、计算步骤二中得到的一次IMF分量的熵权变化曲线图；首先将步骤二中得到的一次IMF分量作为特征量，对其分段，将数据分成10段，即m＝10，每段数据长度n＝40；然后基于公式

$H_j=-\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{p}_i\mathrm{lo}{g}_a{p}_i,j=\mathrm{1,2},...m---\left(10\right)$

其中，H_j为第j个数据段的能量熵，p_i＝E_i/E，E_i是一次IMF信号的瞬时能量，而E为一次IMF信号总能量，n为每段数据长度，m为原始数据的分段总数；

与公式

$q_j=\frac{H_j}{\underset{j=1}{\overset{m}{\mathrm{\Σ}}}{H}_j},j=\mathrm{1,2},\·\·\·,m---\left(11\right)$

计算能量熵与熵权值，并将各点的熵权值连接形成熵权变化曲线图；

步骤三二、计算熵权变化的阈值；熵权变化的阈值是根据正常无故障状态下的陀螺角速度信号的熵权波动性进行设定的；由于标准差可以描述变量的波动情况，因而将由陀螺角速度正常信号计算得到的熵权值的标准差作为熵权变化的阈值；

步骤三三、根据步骤三一与步骤三二得到的故障诊断用信号的熵权变化曲线图和熵权变化阈值，进行故障诊断；

若熵权变化曲线图中的某个时间点或者时间段的熵权值超过熵权变化阈值，则可判断该时间点或者时间段发生故障。