CN104048677A - 基于k-s分布性检验和hht的陀螺仪故障诊断方法 - Google Patents

Abstract

基于K-S分布性检验和HHT的陀螺仪故障诊断方法，本发明涉及一种陀螺仪的故障诊断方法。本发明是要解决现有螺仪故障诊断方法存在不足如产生虚假频率分量，并且故障诊断精度低的问题。步骤一：对原始陀螺角速度输出信号Xp采用EMD方法进行分解，获得不同频段IMF分量；步骤二：对步骤一中得到的不同频段的IMF分量利用K-S分布性检验方法进行相关性检验，判断不同频段的IMF分量是否是原始陀螺角速度输出信号的有效分量；步骤三：对步骤二中经过K-S方法检验过的IMF分量进行HHT变换，进而得到IMF分量的时频谱以及边际谱，结合时频谱上信号的能量与频率变化和边际谱上的信号频率分布判断系统运行过程中是否发生故障。本发明应用于信号处理领域。

Description

基于K-S分布性检验和HHT的陀螺仪故障诊断方法

技术领域

本发明涉及一种陀螺仪的故障诊断方法，具体涉及一种基于Kolmogorov-Smirnov(简称“K-S”)分布性检验与希尔伯特-黄变换(Hilbert-Huang Transform,HHT)的故障诊断方法。

背景技术

近些年，信号处理技术得到了不断的发展，并且由于信号处理学科的不断进步与完善，信号处理的方法也进行了不断的改进和更新。短时傅里叶变换、小波变换等算法将所能处理的对象信号从以前的平稳信号过渡到非平稳信号，但由于这些方法都是基于傅里叶变换的，因而当面对非平稳信号时，它们与傅里叶变换相关的一些不足就会展现出来，例如会产生虚假频率分量。Hilbert-Huang变换(HHT)在处理非平稳信号上表现的很出色，从被提出开始就得到了广泛的研究与应用。这种方法首先利用经验模态分解(EMD)算法对信号进行分解，得到各固有模态函数分量IMF与残差分量RES，然后将各固有模态函数分量IMF进行希尔伯特(Hilbert)变换，得出其时频谱和边际谱，结合时频谱与边际谱上的信号特征分析瞬时频率的特性，由于它可以摆脱傅里叶变换的缺点，因而在非线性、非平稳信号上得到了应用。

已有的基于EMD的故障诊断方法可有效诊断陀螺仪故障，但是在某些情况下，信号的某些时域特征信息可能不是很明显，因而会使得该方法的分析效果不好，如果转换到频域去分析信号特征可能会有所改善。

现有基于EMD的陀螺仪故障诊断方法存在不足，如针对时域特征不明显的诊断信号，则该方法进行故障诊断的精度较低。

发明内容

本发明是要解决现有螺仪故障诊断方法存在不足如产生虚假频率分量，并且故障诊断精度低的问题，而提供了基于K-S分布性检验和HHT的陀螺仪故障诊断方法。

基于K-S分布性检验和HHT的陀螺仪故障诊断方法按以下步骤实现：

步骤一：对原始陀螺角速度输出信号Xp采用EMD方法进行分解，获得不同频段IMF分量；

步骤二：对步骤一中得到的不同频段的IMF分量利用K-S分布性检验方法进行相关性检验，判断不同频段的IMF分量是否是原始陀螺角速度输出信号的有效分量；

步骤三：对步骤二中经过K-S方法检验过的IMF分量进行HHT变换，进而得到IMF分量的时频谱以及边际谱，结合时频谱上信号的能量与频率变化和边际谱上的信号频率分布判断系统运行过程中是否发生故障，并确定出故障发生的时刻。

发明效果：

本发明在已有基于EMD的陀螺仪故障诊断方法基础上引入了以EMD算法为核心的基于K-S分布性检验与希尔伯特-黄变换(HHT)方法，将K-S分布性检验与希尔伯特-黄变换(HHT)相结合，它在频域上对信号进行分析，并且是一种单信号处理方法，这样可以通过时、频信号分析方法的结合，更好地实现陀螺仪的故障诊断，这是首次将基于K-S分布性检验与HHT变换的方法应用于陀螺仪的故障诊断。

对陀螺角速度输出信号，在原始经验模态分解(EMD)算法的基础上，引入了K-S分布性检验与希尔伯特-黄变换，形成了以传统EMD算法为核心的基于K-S分布检验与希尔伯特-黄变换(HHT)的故障诊断方法，通过时、频信号分析方法的结合，可对阶跃、卡死、缓变等故障进行准确检测与识别，它在频域上对原始陀螺输出角速度数据进行分析，适用于单一信号的过程检测，且对于非平稳信号的故障诊断具有较强的优势，是首次创新性地将基于K-S分布性检验与HHT的方法应用在陀螺仪故障诊断上。

1)本发明提出的K-S分布性检验方法，对采用传统EMD方法得到的不同频段的IMF分量进行相关性检验，剔除与原始信号分布相关性较小的分量，从而可在最大程度上确保得到的IMF分量是原始信号的真实分量；

2)本发明采用希尔伯特-黄变换(HHT)，在频域上对信号进行分析，并采用时频谱与边际谱，根据时频谱上信号的能量与频率变化和边际谱上的信号频率分布来判断系统运行过程中是否发生故障，并确定出故障发生的时刻，简化了检测过程，并能较好地完成故障诊断。

附图说明

图1为本发明的流程图；

图2为本发明的正常陀螺角速度输出信号的合成时频谱；

图3为本发明的正常陀螺角速度输出信号的合成边际谱；

图4为本发明的发生单间歇故障时陀螺角速度输出信号的合成时频谱；

图5为本发明的发生单间歇故障时陀螺角速度输出信号的合成边际谱；

图6为本发明的发生多间歇故障时陀螺角速度输出信号的合成时频谱；

图7为本发明的发生多间歇故障时陀螺角速度输出信号的合成边际谱；

图8为本发明的发生持续性故障时陀螺角速度输出信号的合成时频谱；

图9为本发明的发生持续性故障时陀螺角速度输出信号的合成边际谱；

图10为本发明的发生卡死故障时陀螺角速度输出信号的合成时频谱；

图11为本发明的发生卡死故障时陀螺角速度输出信号的合成边际谱；

图12为本发明的发生缓变故障时陀螺角速度输出信号的合成时频谱；

图13为本发明的发生缓变故障时陀螺角速度输出信号的合成边际谱。

具体实施方式

具体实施方式一：本实施方式的基于K-S分布性检验和HHT的陀螺仪故障诊断方法按以下步骤实现：

步骤一：对原始陀螺角速度输出信号Xp采用EMD方法进行分解，获得不同频段IMF分量；

步骤二：对步骤一中得到的不同频段的IMF分量利用K-S分布性检验方法进行相关性检验，判断不同频段的IMF分量是否是原始陀螺角速度输出信号的有效分量；

步骤三：对步骤二中经过K-S方法检验过的IMF分量进行HHT变换，进而得到IMF分量的时频谱以及边际谱，结合时频谱上信号的能量与频率变化和边际谱上的信号频率分布判断系统运行过程中是否发生故障，并确定出故障发生的时刻。

本实施方式效果：

本实施方式在已有基于EMD的陀螺仪故障诊断方法基础上引入了以EMD算法为核心的基于K-S分布性检验与希尔伯特-黄变换(HHT)方法，将K-S分布性检验与希尔伯特-黄变换(HHT)相结合，它在频域上对信号进行分析，并且是一种单信号处理方法，这样可以通过时、频信号分析方法的结合，更好地实现陀螺仪的故障诊断，这是首次将基于K-S分布性检验与HHT变换的方法应用于陀螺仪的故障诊断。

对陀螺角速度输出信号，在原始经验模态分解(EMD)算法的基础上，引入了K-S分布性检验与希尔伯特-黄变换，形成了以传统EMD算法为核心的基于K-S分布检验与希尔伯特-黄变换(HHT)的故障诊断方法，通过时、频信号分析方法的结合，可对阶跃、卡死、缓变等故障进行准确检测与识别，它在频域上对原始陀螺输出角速度数据进行分析，适用于单一信号的过程检测，且对于非平稳信号的故障诊断具有较强的优势，是首次创新性地将基于K-S分布性检验与HHT的方法应用在陀螺仪故障诊断上。

1)本实施方式提出的K-S分布性检验方法，对采用传统EMD方法得到的不同频段的IMF分量进行相关性检验，剔除与原始信号分布相关性较小的分量，从而可在最大程度上确保得到的IMF分量是原始信号的真实分量；

2)本实施方式采用希尔伯特-黄变换(HHT)，在频域上对信号进行分析，并采用时频谱与边际谱，根据时频谱上信号的能量与频率变化和边际谱上的信号频率分布来判断系统运行过程中是否发生故障，并确定出故障发生的时刻，简化了检测过程，并能较好地完成故障诊断。

具体实施方式二：本实施方式与具体实施方式一不同的是：所述步骤一中对原始陀螺角速度输出信号Xp采用EMD方法进行分解，获得不同频段IMF分量具体为：

(一)利用matlab极值函数找出原始陀螺角速度输出信号Xp时间序列上的所有的局部极值点；

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

(三)生成Xp的上包络函数和下包络函数的均值函数：基于公式(1)得到Xp的上包络函数和下包络函数的均值函数m₁(t)；

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

(四)基于公式(2)求得信号Xp与包络均值函数m₁(t)的差值函数：

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

(五一)筛选获得第一个IMF分量：

通常情况下，h₁(t)不满足IMF条件，将h₁(t)作为原始信号，重复(一)～(三)步骤得到Xp的包络均值函数m₁₁(t)，求得差值函数h₁₁(t)＝h₁(t)-m₁₁(t)，判断h₁₁(t)是否满足IMF的条件，如果不满足则继续重复(一)～(三)步骤，重复到第k次时，得到的差值函数为h_1k(t)＝h_1(k-1)(t)-m_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)$

式中，h_1k(t)为第k次筛选第1个IMF分量时得到的结果，h_1(k-1)(t)为第k-1次筛选第1个IMF分量时得到的结果，T为所选数据的时间长度，k为第k次筛选过程；

根据公式(3)依次计算两个筛选出来的数据，得到其标准差值并与既定阈值作比较，若标准差值小于既定阈值，则筛选过程结束，确定第一个IMF分量，即c₁(t)＝h_1k(t)，如果标准差值不小于既定阈值，则继续筛选到标准差值小于既定阈值为止；

(五二)获得不同频段IMF分量：

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

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

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

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

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

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

具体实施方式三：本实施方式与具体实施方式一或二不同的是：所述步骤二中得到的不同频段的IMF分量利用K-S分布性检验方法进行相关性检验，判断不同频段的IMF分量是否是原始陀螺角速度输出信号的有效分量：

(一)K-S分布性检验：

对于N点时间序列y(n)＝{y₁,y₂,…,y_n}，定义其累积分布函数：

$E_i=\frac{n\left(i\right)}{N}---\left(6\right)$

其中，n(i)是将包含N个数据的数据样本首先进行升序排序，得到的在样本总体中数据值小于y(i)的样本数；

设f(x)和r(x)分别为两个信号的累积分布函数，定义f(x)和r(x)在同一数据点处的最大差异D：

$D=\underset{-\∞\≤x\≤+\∞}{\max }|f\left(x\right)-r\left(x\right)|---\left(7\right)$

进一步定义两组数据样本的相似概率prob(D)如下：

$\mathrm{prob}\left(D\right)=Q_{\mathrm{ks}}\left(\mathrm{\λ}\right)=2\underset{j=1}{\overset{\∞}{\mathrm{\Σ}}}{(-1)}^{j-1}{e}^{-2j^2{\mathrm{\λ}}^2}---\left(8\right)$

其中，Q_ks(λ)为K-S概率分布函数，N₁,N₂分别为两序列含有的样本总数，j为虚数单位；

通过式(8)，当λ→0时，Q_KS→1；当λ→∞时，Q_KS→0；得到结论：如果两组信号的累积分布函数上相似，二者的相似概率趋于1，如果累积分布函数不相似，二者的相似概率趋于0；

(二)判断不同频段的IMF分量是否是原始陀螺角速度输出信号的有效分量：

(1):求取累积分布函数：基于公式(6)求取原始陀螺角速度输出信号Xp与步骤一中得到的不同频段的IMF分量的累计分布函数；

(2):求累积分布函数的最大差异：将原始陀螺角速度输出信号Xp作为参考信号，将步骤一中得到的不同频段的IMF分量与参考信号做比较，基于公式(7)，求取所有IMF分量与原始信号Xp的累计分布函数之间的最大差异，记为D₁,D₂,…D_n；

(3):求相似概率值：基于公式(8)求取不同频段的IMF分量与原始陀螺角速度输出数据信号之间的相似概率值。

Kolmogorov-Smirnov分布性检验法的基本理论：在统计领域中有一类问题是假设检验问题，它通常包含两种形式：一种是完全未知总体分布，第二种是已知总体分布，未知所含参数。因而提出了统计假设来了解总体的某些特性，其中的假设可以来自于理论的分析，也可以来自于实际问题的观测。通过已知的样本的信息来对假设进行判断，其中参数假设检验是对总体参数提出的假设，除此之外为非参数假设检验。

通常总体分布的类型在参数假设检验中被假定已知并满足正态分布。但在实际对信号分析时，信号总体服从的分布一般是未知的，这时需要根据获得的样本来对信号总体进行假设检验，即进行非参数假设检验。在文中用到的Kolmogorov-Smirnov分布性检验法即是非参数假设检验的一种。

Kolmogorov-Smirnov分布拟合优度检验法，简称K-S分布性检验法，用于判断两个样本概率分布函数差异，描述两个独立统计样本的相关性。

这样可以知道如果两个时间序列类型不同，那么其累积分布函数差异将非常明显，可以根据这来判断两个序列的分布是否相似。

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

具体实施方式四：本实施方式与具体实施方式一至三之一不同的是：所述步骤三：对步骤二中经过K-S方法检验过的IMF分量进行HHT变换，进而得到IMF分量的时频谱以及边际谱，结合时频谱上信号的能量与频率变化和边际谱上的信号频率分布判断系统运行过程中是否发生故障，并确定出故障发生的时刻：

一、Hilbert-Huang时频谱和边际谱：

利用EMD将原始信号x(t)进行分解得到不同频段的IMF分量，各频段IMF分量都是单分量信号，接着对IMF分量进行Hilbert变换，具体为：

${\hat{c}}_i\left(t\right)=\frac{1}{\mathrm{\π}}\underset{-\∞}{\overset{+\∞}{\∫}}\frac{c_i\left(\mathrm{\τ}\right)}{t-\mathrm{\τ}}\mathrm{d\τ}---\left(13\right)$

式中，c_i(τ),i＝1,2,…n为EMD分解得到的不同频段的IMF分量，为Hilbert变换后的不同频段的IMF分量；

随后构造解析信号如下:

其中A_i(t)和分别表示z_i(t)的瞬时幅值和相位，各IMF分量的瞬时频率定义为：

那么原信号x(t)表示为：

用时间和瞬时频率来描述信号x(t)的幅值，Hilbert时频谱定义为：

$H(\mathrm{\ω},t)=\mathrm{Re}\left[\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{A}_i\left(t\right)e^{j{\∫\mathrm{\ω}}_i\left(t\right)\mathrm{dt}}\right]---\left(17\right)$

由信号的时频谱可以进一步求得其边际谱，边际谱为：

$h\left(\mathrm{\ω}\right)={\∫}_0^{T}H(\mathrm{\ω},t)\mathrm{dt}---\left(18\right)$

二、具体步骤如下：

(1)：Hilbert变换：基于式(13)，计算步骤二中经过K-S分布性检验方法得到的IMF分量的Hilbert变换信号，记作：并构造解析信号，如式(14)所示；

(2):求取信号瞬时频率：对Step(1)中得到的不同频段IMF分量的解析信号z_i(t),i＝1,2,…n，基于公式(15)，求取瞬时频率，记作：ω₁(t),ω₂(t),…,ω_n(t)；

(3):求时频谱以及边际谱并进行故障诊断：基于(1)与(2)中计算得到的不同频段IMF分量的解析信号与瞬时频率，分别基于式(17)与式(18)求取原始陀螺角速度输出信号的时频谱与边际谱。

本实施方式中，HHT变换基本理论

1、Hilbert变换

Hilbert变换被广泛应用在信号瞬时特性(瞬时振幅、瞬时相位、瞬时频率)以及非线性与非平稳性分析等方面，是信号分析中的重要工具。对于给定的连续时间信号x(t)，其Hilbert变换定义为：

$\tilde{x}\left(t\right)=\frac{1}{\mathrm{\π}}{\∫}_{-\∞}^{+\∞}\frac{x\left(\mathrm{\τ}\right)}{t-\mathrm{\τ}}\mathrm{d\τ}=\frac{1}{\mathrm{\π}}{\∫}_{-\∞}^{+\∞}-\frac{x(t-\mathrm{\τ})}{\mathrm{\τ}}\mathrm{d\τ}=x\left(t\right)*\frac{1}{\mathrm{\πt}}---\left(9\right)$

由上式可以看出的卷积，因而可以将Hilbert看成是x(t)通过一个线性时不变系统的输出，该系统的单位脉冲响应为：

以原始信号x(t)为实部，Hilbert变换后的信号为虚部，构造解析信号：

$Z\left(t\right)=x\left(t\right)+j\tilde{x}\left(t\right)---\left(10\right)$

Z(t)可进一步表示成：

$Z\left(t\right)=x\left(t\right)+j\tilde{x}\left(t\right)=a\left(t\right)e^{\mathrm{j\θ}\left(t\right)}---\left(11\right)$

其中，a(t)和θ(t)分别表示Z(t)的瞬时幅值和相位。

2、瞬时频率

Norden E.Huang等人指出：只有在模态函数中，计算瞬时频率才具有非常明显的物理意义。频率这一概念，在传统的傅里叶分析中是这样定义的：在整个过程数据中，具有恒定幅值的正弦或者余弦函数，其瞬时频率应该与正弦或者余弦函数有关。如果一个信号长度过短，少于一个波长，那么频率无从谈起，对于非平稳信号而言，因为其频率会变化较大，这时的瞬时频率显得特别重要。瞬时频率可以被理解成被研究信号的频率值用正弦波的局部最佳来逼近，如果被研究信号含有多个频段分量或者其信号频率会随着时间变化，那么这时，瞬时频率就不会是单值了。这样为了使瞬时频率更具有意义，需要首先将这些含有多频段的信号分解为单成分信号。

对于1中所述的Hilbert变换，是原始信号x(t)和1/πt的卷积，所以构造的解析信号Z(t)利用随时间变化的幅值和相位三角函数来局部最佳拟合信号x(t)，这时瞬时频率被定义为：

$\mathrm{\ω}\left(t\right)=\frac{\mathrm{d\θ}\left(t\right)}{\mathrm{dt}}---\left(12\right)$

由上式可以看出，信号的瞬时频率值在每一个时刻都是唯一的，并且瞬时频率与时间是单值映射的。

3、希尔伯特-黄(HHT)变换

HHT的基本思想源于信号瞬时频率的求解。关于信号瞬时频率的定义方法有多种，其中普遍接受的一种定义是对原始信号进行如1中所述的Hilbert变换，再由原始信号与Hilbert变换信号构成解析信号，采用如2中所述的瞬时频率求解方法，即构造的解析信号的相位部分的导数作为信号的瞬时频率。但是这种方法对于任何一个信号在每一时刻只能求得一个瞬时频率，因此当信号中包含多个频率成分时，就不能给出合理的结果。

因而HHT基于上述问题，将求解信号瞬时频率的过程分成两步：首先将原始信号经过经验模态分解(EMD)分解成不同频段的IMF分量和残差分量RES，然后再对每个IMF分量进行Hilbert变换，求瞬时频率。

从Hilbert边际谱中可以看到信号的幅值随频率变化的情况，描述了信号在概率上的累积幅值，并且是对信号中各不同频率成分幅值的一种量度。在Hilbert边际谱上，如果信号在某个频率处存在着幅值，那么说明原信号中含有这个频率成分。

对于正常陀螺角速度输出信号，其时频谱中，信号能量值分布差别不大(通过颜色值)，并且随着时间的变化，信号瞬时频率变化较为光滑连续；边际谱中，信号主要分布在低频。当陀螺发生故障时，陀螺角速度输出数据信号的时频谱中，信号能量值会发生跃变，同时会发生频率突变；边际谱中，故障相应频段上的幅值分量会增加。

因而通过正常与故障时的陀螺角速度输出数据信号的时频谱与边际谱的比较,综合考虑时频谱中能量跃变与频率突变的时间点或者时间段以及边际谱中信号幅值分量增大的时间点或者时间段可有效进行故障诊断。

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

具体实施方式五：本实施方式与具体实施方式一至四之一不同的是：所述步骤(五一)中IMF条件为：

(1)在整个所考虑信号范围上，存在极值点和过零点，并且极值点和过零点在数量上相等，或者只相差一个；

(2)在任意时刻，由极小值包络与极大值包络二者形成的包络均值为零。

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

具体实施方式六：本实施方式与具体实施方式一至五之一不同的是：所述步骤(五一)中所述SD的既定阈值取0.2～0.3。

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

仿真实验：

执行步骤一：卫星姿态控制系统各组成部分的模型在Matlab平台上Simulink环境下搭建闭环仿真系统，其中仿真步长为0.025s，对仿真系统运行，得到陀螺角速度输出信号X_p。对陀螺角速度输出数据X_p基于经验模态分解(EMD)方法进行分解,得到不同频段IMF分量以及残差分量RES。

执行步骤二：对步骤一中得到的各个不同频段的IMF分量利用K-S分布性检验方法进行检验，来判断各个不同频段的IMF分量是否是原始陀螺角速度输出信号的有效分量。各个不同频段的IMF分量与原始陀螺角速度输出信号X_p之间的相似概率值如表1所示：

表1正常情况下相似概率值计算

从表1可看到，正常情况下，得到各个IMF分量与原信号的分布相关性很强，可以认为是有效分量。

执行步骤三：对步骤二中得到的经过K-S方法检验过的IMF分量进行希尔伯特-黄变换，进而得到IMF分量的时频谱以及边际谱，结合时频谱上信号的能量与频率变化和边际谱上的信号频率分布来判断系统运行过程中是否发生故障，并确定出故障发生的时刻。

对步骤二中得到的各个IMF分量进行叠加，合成的信号时频谱与边际谱分别如图2、图3所示。在时频谱中看到正常状态下信号能量值分布差别不大(通过颜色值)，并且随着时间的变化，信号瞬时频率变化较为光滑连续，在边际谱中发现该陀螺信号主要分布在低频。

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

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

当偏航轴陀螺发生阶跃故障情况时，又可具体分为两种情况：间歇性阶跃故障与持续性阶跃故障。下面将分别对这两种情况进行方法的验证。

(1)发生间歇性故障情况

当发生间歇性阶跃故障时，又具体分为发生的是单故障情况还是多故障情况。当发生单故障时，这里在原始陀螺角速度输出信号的数据区间t＝15s～16s上添加故障值为0.04rad/s的阶跃故障，并作为故障诊断用信号。

首先利用EMD算法对故障诊断用信号进行分解，得到不同频段的IMF分量及残差分量RES，然后基于K-S分布性检验方法计算IMF分量与故障诊断用信号之间的相关性，得到其相似概率值如表2所示。

表2单间歇故障相似概率值计算

结果分析，由表2可以看到各IMF分量与故障诊断用信号之间的相关性很强，可认为是有效分量。由各IMF分量叠加合成的信号的时频谱与边际谱分别如图4、图5所示。通过图4时频谱可以看到，在t＝15s～16s区间内信号发生了明显的频率突变，并且其能量值也发生了跃变；而在图5边际谱中看到，与正常状态下相比，信号在低频段区域，相对较高一些的频段上增加了分量，由此断定信号在该过程中发生了故障。

当发生多故障时，以在过程中发生两次故障为例，这里在原始陀螺角速度输出信号的数据区间t＝6s～7s、t＝15s～16s上添加故障值为0.04rad/s的阶跃故障，并作为故障诊断用信号。

首先利用EMD算法对故障诊断用信号进行分解，得到不同频段的IMF分量及残差分量RES，然后基于K-S分布性检验方法计算IMF分量与故障诊断用信号之间的相关性，得到其相似概率值如表3所示。

表3多间歇性故障相似概率值计算

结果分析，由表3可以看到，各分量与故障诊断用信号之间在概率分布上的相关性很强，这样由各分量叠加合成信号的时频谱与边际谱如图6、图7所示。从时频谱中看到在t＝6s～7s，t＝15s～16s区间上，信号的能量值发生了突变，并且其频率发生大幅度突变，而其它区间上的信号分布较为光滑稳定，分析其边际谱，在低频段，相对较高一些的频率区间上增加了信号分量，由此断定发生了故障。

(2)持续性故障情况

当发生持续性阶跃故障时，这里在原始陀螺角速度输出信号的数据点t＝16s时刻添加故障值为0.04rad/s的阶跃故障，并作为故障诊断用信号。

首先利用EMD算法对故障诊断用信号进行分解，，得到不同频段的IMF分量及残差分量RES，然后基于K-S分布性检验方法计算IMF分量与故障诊断用信号之间的相关性，得到其相似概率值如表4所示

表4持续故障时相似概率值计算

结果分析，由表4可以看到，各分量与故障诊断用信号之间在概率分布上的相关性很强，由这几个分量合成的信号的时频谱和边际谱分别如表4所示。在时频谱的t＝16s时刻频率发生了突变，时刻前后的信号能量值也发生了很大的变化，边际谱较高频段上增加了信号值，因而判断此刻发生了故障。

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

当陀螺在运行过程中遇到特殊情况而发生卡死故障时，陀螺的角速度信号输出为零，这里设置在原始陀螺角速度输出信号的数据点t＝15s时刻陀螺发生卡死故障，并作为故障诊断用信号。

首先利用EMD算法对故障诊断用信号进行分解，，得到不同频段的IMF分量及残差分量RES，然后基于K-S分布性检验方法计算IMF分量与故障诊断用信号之间的相关性，得到其相似概率值如表5所示

表5卡死故障相似概率值计算

结果分析，由表5可以看出，各IMF分量与故障诊断用信号之间的相关性较小，但统计学上认为相关性大于0.5即可认为相关性较强。由这几个分量得到的合成信号的时频谱与边际谱分别如表4所示。从时频谱可以看到，当在t＝15s时刻发生卡死故障时，信号在频率上发生了突变，时刻前后的信号能量值也发生了很大的变化，在边际谱上看到在低频段，相对较高一些的频段上增加了分量，因而判断此刻发生了故障。

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

当陀螺在运行过程中发生缓变漂移故障的原因主要分为内因和外因，内因主要是由于陀螺本身内部结构由于生产时的工艺水平不足而形成的各种干扰力矩；外因可能较多，如陀螺运行转速未在理想转速情况下运行、测量时的温度影响、工作场合不满足水平等因素，这些都会导致其获得的角速度信号发生缓变漂移。在这里设置在原始陀螺角速度输出信号的数据点t＝15.5s时刻陀螺发生缓变故障，缓变速率为1.5×10^-4rad/s²，并作为故障诊断用信号。

首先利用EMD算法对故障诊断用信号进行分解，得到不同频段的IMF分量及残差分量RES，然后基于K-S分布性检验方法计算IMF分量与故障诊断用信号之间的相关性，得到其相似概率值如表6所示

表6缓变故障相似概率值计算

结果分析，由表6可以看到各个IMF分量与原信号之间的相关性很强，所以认为是有效分量。由各分量合成的信号的时频谱与边际谱分别如表5所示。通过时频谱可以看到，当在t＝15.5s时刻发生缓变故障时，信号频率上发生了突变，时刻前后的信号能量值也发生了很大的变化，在边际谱上看到在较高频段上增加了信号值，因而判断此刻发生了故障。

本发明详细说明了如何应用K-S分布性检验与HHT进行陀螺仪故障诊断。综合实施例的上述分析，对于陀螺仪的故障诊断，本发明的算法能够有效克服传统的基于经验模态分解(EMD)方法的不足，可以有效地对阶跃、缓变、卡死等故障进行诊断。

Claims (6)

1.基于K-S分布性检验和HHT的陀螺仪故障诊断方法，其特征在于基于K-S分布性检验和HHT的陀螺仪故障诊断方法按以下步骤实现：

步骤一：对原始陀螺角速度输出信号Xp采用EMD方法进行分解，获得不同频段IMF分量；

步骤二：对步骤一中得到的不同频段的IMF分量利用K-S分布性检验方法进行相关性检验，判断不同频段的IMF分量是否是原始陀螺角速度输出信号的有效分量；

步骤三：对步骤二中经过K-S方法检验过的IMF分量进行HHT变换，进而得到IMF分量的时频谱以及边际谱，结合时频谱上信号的能量与频率变化和边际谱上的信号频率分布判断系统运行过程中是否发生故障，并确定出故障发生的时刻。

2.根据权利要求1所述的基于K-S分布性检验和HHT的陀螺仪故障诊断方法，其特征在于所述步骤一中对原始陀螺角速度输出信号Xp采用EMD方法进行分解，获得不同频段IMF分量具体为：

(一)利用matlab极值函数找出原始陀螺角速度输出信号Xp时间序列上的所有的局部极值点；

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

(三)生成Xp的上包络函数和下包络函数的均值函数：基于公式(1)得到Xp的上包络函数和下包络函数的均值函数m₁(t)；

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

(四)基于公式(2)求得信号Xp与包络均值函数m₁(t)的差值函数：

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

(五一)筛选获得第一个IMF分量：

通常情况下，h₁(t)不满足IMF条件，将h₁(t)作为原始信号，重复(一)～(三)步骤得到Xp的包络均值函数m₁₁(t)，求得差值函数h₁₁(t)＝h₁(t)-m₁₁(t)，判断h₁₁(t)是否满足IMF的条件，如果不满足则继续重复(一)～(三)步骤，重复到第k次时，得到的差值函数为h_1k(t)＝h_1(k-1)(t)-m_1k(t)；

确定迭代停止准则，引入标准差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)$

式中，h_1k(t)为第k次筛选第1个IMF分量时得到的结果，h_1(k-1)(t)为第k-1次筛选第1个IMF分量时得到的结果，T为所选数据的时间长度，k为第k次筛选过程；

根据公式(3)依次计算两个筛选出来的数据，得到其标准差值并与既定阈值作比较，若标准差值小于既定阈值，则筛选过程结束，确定第一个IMF分量，即c₁(t)＝h_1k(t)，如果标准差值不小于既定阈值，则继续筛选到标准差值小于既定阈值为止；

(五二)获得不同频段IMF分量：

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

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

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

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

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

3.根据权利要求2所述的基于K-S分布性检验和HHT的陀螺仪故障诊断方法，其特征在于所述步骤二中得到的不同频段的IMF分量利用K-S分布性检验方法进行相关性检验，判断不同频段的IMF分量是否是原始陀螺角速度输出信号的有效分量：

(一)K-S分布性检验：

对于N点时间序列y(n)＝{y₁,y₂,…,y_n}，定义其累积分布函数：

$E_i=\frac{n\left(i\right)}{N}---\left(6\right)$

其中，n(i)是将包含N个数据的数据样本首先进行升序排序，得到的在样本总体中数据值小于y(i)的样本数；

设f(x)和r(x)分别为两个信号的累积分布函数，定义f(x)和r(x)在同一数据点处的最大差异D：

$D=\underset{-\∞\≤x\≤+\∞}{\max }|f\left(x\right)-r\left(x\right)|---\left(7\right)$

进一步定义两组数据样本的相似概率prob(D)如下：

$\mathrm{prob}\left(D\right)=Q_{\mathrm{ks}}\left(\mathrm{\λ}\right)=2\underset{j=1}{\overset{\∞}{\mathrm{\Σ}}}{(-1)}^{j-1}{e}^{-2j^2{\mathrm{\λ}}^2}---\left(8\right)$

其中，Q_ks(λ)为K-S概率分布函数，N₁,N₂分别为两序列含有的样本总数，j为虚数单位；

通过式(8)，当λ→0时，Q_KS→1；当λ→∞时，Q_KS→0；得到结论：如果两组信号的累积分布函数上相似，二者的相似概率趋于1，如果累积分布函数不相似，二者的相似概率趋于0；

(二)判断不同频段的IMF分量是否是原始陀螺角速度输出信号的有效分量：

(1):求取累积分布函数：基于公式(6)求取原始陀螺角速度输出信号Xp与步骤一中得到的不同频段的IMF分量的累计分布函数；

(2):求累积分布函数的最大差异：将原始陀螺角速度输出信号Xp作为参考信号，将步骤一中得到的不同频段的IMF分量与参考信号做比较，基于公式(7)，求取所有IMF分量与原始信号Xp的累计分布函数之间的最大差异，记为D₁,D₂,…D_n；

(3):求相似概率值：基于公式(8)求取不同频段的IMF分量与原始陀螺角速度输出数据信号之间的相似概率值。

4.根据权利要求3所述的基于K-S分布性检验和HHT的陀螺仪故障诊断方法，其特征在于所述步骤三：对步骤二中经过K-S方法检验过的IMF分量进行HHT变换，进而得到IMF分量的时频谱以及边际谱，结合时频谱上信号的能量与频率变化和边际谱上的信号频率分布判断系统运行过程中是否发生故障，并确定出故障发生的时刻：

一、Hilbert-Huang时频谱和边际谱：

利用EMD将原始信号x(t)进行分解得到不同频段的IMF分量，各频段IMF分量都是单分量信号，接着对IMF分量进行Hilbert变换，具体为：

${\hat{c}}_i\left(t\right)=\frac{1}{\mathrm{\π}}\underset{-\∞}{\overset{+\∞}{\∫}}\frac{c_i\left(\mathrm{\τ}\right)}{t-\mathrm{\τ}}\mathrm{d\τ}---\left(13\right)$

式中，c_i(τ),i＝1,2,…n为EMD分解得到的不同频段的IMF分量，为Hilbert变换后的不同频段的IMF分量；

随后构造解析信号如下:

其中A_i(t)和分别表示z_i(t)的瞬时幅值和相位，各IMF分量的瞬时频率定义为：

那么原信号x(t)表示为：

用时间和瞬时频率来描述信号x(t)的幅值，Hilbert时频谱定义为：

$H(\mathrm{\ω},t)=\mathrm{Re}\left[\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{A}_i\left(t\right)e^{j{\∫\mathrm{\ω}}_i\left(t\right)\mathrm{dt}}\right]---\left(17\right)$

由信号的时频谱可以进一步求得其边际谱，边际谱为：

$h\left(\mathrm{\ω}\right)={\∫}_0^{T}H(\mathrm{\ω},t)\mathrm{dt}---\left(18\right)$

二、具体步骤如下：

(1)：Hilbert变换：基于式(13)，计算步骤二中经过K-S分布性检验方法得到的IMF分量的Hilbert变换信号，记作：并构造解析信号，如式(14)所示；

(2):求取信号瞬时频率：对Step(1)中得到的不同频段IMF分量的解析信号z_i(t),i＝1,2,…n，基于公式(15)，求取瞬时频率，记作：ω₁(t),ω₂(t),…,ω_n(t)；

(3):求时频谱以及边际谱并进行故障诊断：基于(1)与(2)中计算得到的不同频段IMF分量的解析信号与瞬时频率，分别基于式(17)与式(18)求取原始陀螺角速度输出信号的时频谱与边际谱。

5.根据权利要求4所述的基于K-S分布性检验和HHT的陀螺仪故障诊断方法，其特征在于所述步骤(五一)中IMF条件为：

(1)在整个所考虑信号范围上，存在极值点和过零点，并且极值点和过零点在数量上相等，或者只相差一个；

(2)在任意时刻，由极小值包络与极大值包络二者形成的包络均值为零。

6.根据权利要求5所述的基于K-S分布性检验和HHT的陀螺仪故障诊断方法，其特征在于所述步骤(五一)中所述SD的既定阈值取0.2～0.3。