基于短样本频谱校正旋转机械振动信号盲分离方法及装置的制造方法
【专利摘要】本发明公开了一种基于短样本频谱校正的旋转机械振动信号盲分离方法及装置，方法包括：获取合并后的频率序列中满足给定一个阈值，都存在且仅存在一个下标满足第一表达式的频率，则该频率对应一个候选有效模式；将候选有效模式组成的集合作为输入，识别候选有效模式集合中的单源模式；对单源模式进行基于DB指数的K?均值聚类，获取混合矩阵的估计；通过混合矩阵的估计、观测信号恢复源信号。装置包括：振动传感器、模数转化器、微控制器和输出驱动及其显示模块，经过DSP芯片的内部处理，得到混合矩阵的估计；最后借助输出驱动及其显示模块显示恢复的源信号。本发明解决了旋转机械振动信号在样本数目很少的超定情况下的盲源分离问题。
【专利说明】
基于短样本频谱校正旋转机械振动信号盲分离方法及装置
技术领域
[0001] 本发明涉及旋转机械振动信号盲分离领域，尤其涉及一种基于短样本频谱校正的 旋转机械振动信号盲分离方法及装置，本发明解决了旋转机械振动信号在样本数目很少的 超定情况下的盲源分离问题。
【背景技术】
[0002] 作为大型机械系统中最主要的机械类型，旋转机械在现代工程应用中占有很大比 重。由于现代机械系统的巨型化、精细化、自动化的发展趋势，旋转机械系统常常要工作在 恶劣的环境中，因而导致故障的时常发生。旋转机械出现故障，往往导致整个机械系统的破 坏或者是工作性能的降低，比如生产质量、操作安全性等等。
[0003] 盲源分离(Blind Source Separation,BSS)是指在未知混合系统的信道信息的情 况下，从观测中恢复源信号的问题。盲源分离问题广泛地应用于机械诊断信号处理[14]、语 音识别[5]和无线通信 [6]等领域。众所周知，旋转部件(齿轮、轴承等等)是现代大型机械系统 中最常见、最关键的部件[7]。然而由于现场诸多因素的影响（例如:将多个小型马达固定在 同一结构上，或者多个机械故障同时发生等等），从传感器获得的观测信号，往往不能反映 特定旋转部件的真实工作状态。与此同时，在工程实践中，这些观测信号往往会被环境干扰 (如环境噪声，其他机械系统等等）。因此BSS可以且很有必要作为提取不同源信号特征的预 处理方法 [8]，为机械故障监测与诊断提供基础。
[0004] 通常情况下，由于机械中多个转子的同时工作，旋转机械的振动信号往往由多个 周期谐波成分构成。不同的机械故障类型，产生的振动信号有着不同的频谱特征。例如:转 子不对中所产生的振动信号中，主要的特征是存在明显的二次谐波 [9];齿轮在齿轮座中的 松动常常产生高于10次的谐波成分(有时甚至达到20次谐波）；油膜振荡故障总是生成一些 接近半谐波的次谐波成分 [1()]。
[0005] 为了减小故障事故所产生的损失，在现场故障检测与诊断中迫切需要提高响应速 度，从而及时发现机械故障。快速响应的最直接解决方法是利用尽量短时间的观测信号来 进行故障检测与诊断，即利用短样本信号进行现场故障检测与诊断。
[0006] 然而，作为预处理手段，现有的BSS方法很难满足短样本的要求。例如:旋转机械故 障诊断中最主流的BSS算法独立分量分析（Independent Component Analysis，ICA)[n]。许 多基于ICA的方法[12，13]，及基于改进ICA的方法(如二阶ICA [14]，非线性ICA[15]，kernel ICA
[3] 等等）已经引入到故障识别及分析中。
[0007] 发明人在实现本发明的过程中，发现现有技术中至少存在以下缺点和不足：
[0008] ICA在仅仅提供短样本观测时，常常产生不确定的解。这是由于，ICA是一种基于对 四阶累积量峭度所衍生的一个目标函数的最优化算法。作为随机信号的四阶累积量，峭度 的计算需消耗大量样本以保证其稳定及精确。其他基于统计学的方法，例如:基于四阶张量 分析的F00BI(Fourth Order Cumulant-based Blind Identification)1^16]，同样会在短样 本情况下呈现较差性能。
【发明内容】
[0009] 本发明提供了一种基于短样本频谱校正的旋转机械振动信号盲分离方法及装置， 本发明仅需给定较短的观测样本，即可获得较好的源信号恢复效果，有助于提高旋转机械 监测与诊断系统的快速响应，详见下文描述：
[0010] 一种基于短样本频谱校正的旋转机械振动信号盲分离方法，所述方法包括以下步 骤：
[0011] 获取合并后的频率序列中满足给定一个阈值，都存在且仅存在一个下标满足第一 表达式的频率，则该频率对应一个候选有效模式；
[0012] 将候选有效模式组成的集合作为输入，识别候选有效模式集合中的单源模式；
[0013] 对单源模式进行基于DB指数的K-均值聚类，获取混合矩阵的估计；
[0014] 通过混合矩阵的估计、观测信号恢复源信号。
[0015] 所述第一表达式具体为：
[0016] ^?)\<s
[0017] 其中，U频率估计集;乙为频率;e为阈值。
[0018] 所述将候选有效模式组成的集合作为输入，识别候选有效模式集合中的单源模式 具体为：
[0019] 1)单源模式的幅值向量与混合矩阵的某个列向量平行；
[0020] 2)单源模式的相位向量具有相位一致性，即相位向量的元素满足下式：
[0022]其中，l<r，l<M，r辛1;|>0为阈值;C为排列组合运算;I.为第p*单源模式相位向 量中第r个元素;为第!/单源模式相位向量中第1个元素为单源模式的下标;M为观测 信号的个数。所述对单源模式进行基于DB指数的K-均值聚类，获取混合矩阵的估计具体为： [0023]将聚类数目初始化为1 = 2,利用K-均值聚类法将单源模式集聚成I类，并计算DB指 数；
[0024] 令I = 1+1，并再次聚类，求DB指数;直到寻找到I满足DKDh，且DKDI+1，则I可作为 源信号数目N的估计，而聚类中心所组成的矩阵，作为混合矩阵的估计。
[0025] 所述获取合并后的频率序列具体为：
[0026] 对每个观测信号分别做频谱校正，获取若干个频率估计集合；
[0027] 将频率估计集合进行混合，按照升序排列，获取合并后的频率序列。
[0028] 所述对每个观测信号分别做频谱校正，获取若干个频率估计集合具体为：
[0029] 对观测信号做加汉宁窗L点DFT变换，得到观测频谱；
[0030] 收集观测频谱的所有峰值的位置，计算峰值及其次峰的比值，获取中间参数：
[0031] 根据所述中间参数估计频率偏移:通过所述频率偏移获取频率估计、幅值估计及 相位估计。
[0032] -种基于短样本频谱校正的旋转机械振动信号盲分离装置，所述装置包括:振动 传感器、模数转化器、微控制器和输出驱动及其显示模块，
[0033]所述振动传感器采集到的多路观测信号x(t)经过所述模数转化器采样得到样本 序列，以并行数字输入的形式进入所述DSP芯片，经过所述DSP芯片的内部算法处理，得到混 合矩阵的估计;最后借助所述输出驱动及其显示模块显示恢复的源信号。
[0034]所述微控制器具体为DSP芯片。
[0035] 本发明提供的技术方案的有益效果是:第一、与其他盲分离方法(如fast-ICA)相 比，本发明能够在短样本观测的条件下实现高精度的源信号恢复，而短样本观测十分符合 实际当中机械检测与诊断的快速响应的要求。
[0036] 第二、本发明所引入的频谱校正技术十分擅长于谐波信息的提取，因而尤其适合 于旋转机械振动信号的处理，因为旋转机械振动信号的故障往往与谐波信息相关。
[0037] 第三、本发明所提出的BSS方法可以通过改进的K-均值聚类法精确估计出源信号 数目，而实际应用于旋转机械信号处理时，源信号数目也常常未知，因而本发明与实际应用 更加符合。
[0038]第四、不同于基于统计特性的方法(如€&^-10六^0081，等等），本发明不需随机初 始化过程，也不需多次迭代，因而具有较高的稳定性，以及较低的复杂性，使得旋转机械振 动信号的分析更加可靠，更加有效。
【附图说明】
[0039]图1为基于短样本频谱校正的旋转机械振动信号盲分离方法的流程图；
[0040] 图2为基于短样本频谱校正的旋转机械振动信号盲分离方法的详细流程图；
[0041] 图3为fast-ICA合成信号长样本(L = 400)恢复信号的示意图；
[0042]图4为本发明合成信号长样本(L = 400)恢复信号的示意图；
[0043]图5为fast-ICA合成信号短样本(L = 70)恢复信号的示意图；
[0044] 图6为本发明合成信号短样本(L = 70)恢复信号的示意图；
[0045] 图7为fast-ICA实际振动信号长样本(L = 400)恢复信号的示意图；
[0046] 图8为本发明实际振动信号长样本(L = 400)恢复信号的示意图；
[0047]图9为fast-ICA实际振动信号短样本(L = 70)恢复信号的示意图；
[0048]图10为本发明实际振动信号短样本(L = 70)恢复信号的示意图；
[0049] 图11为基于短样本频谱校正的旋转机械振动信号盲分离装置的结构示意图；
[0050] 图12为DSP内部的流程图。
【具体实施方式】
[0051] 为使本发明的目的、技术方案和优点更加清楚，下面对本发明实施方式作进一步 地详细描述。
[0052]为了解决【背景技术】中存在的问题和不足，本发明实施例提出了短样本情况下的基 于短样本频谱校正的旋转机械信号盲分离方法及装置，结合了频谱校正、相位一致性准则 及改进的K-均值聚类法。同时本发明也能够工作在较长样本情形。
[0053] 实施例1
[0054] 一种基于短样本频谱校正的旋转机械振动信号盲分离方法，参见图1，该方法包括 以下步骤：
[0055] 101:获取合并后的频率序列中满足给定一个阈值，都存在且仅存在一个下标满足 第一表达式的频率，则该频率对应一个候选有效模式；
[0056] 102:将候选有效模式组成的集合作为输入，识别候选有效模式集合中的单源模 式；
[0057] 103:对单源模式进行基于DB指数的K-均值聚类，获取混合矩阵的估计；
[0058] 104:通过混合矩阵的估计、观测信号恢复源信号。
[0059]其中，步骤101中的第一表达式具体为：
[0060] |U 卜
[0061 ]其中，为频率估计集;尤为频率;e为阈值。
[0062]其中，步骤102中的将候选有效模式组成的集合作为输入，识别候选有效模式集合 中的单源模式具体为：
[0063] 1)单源模式的幅值向量与混合矩阵的某个列向量平行；
[0064] 2)单源模式的相位向量具有相位一致性，即相位向量的元素满足下式：
[0066]其中，l<r，l<M，r辛1;|>0为阈值;C为排列组合运算；为第！/单源模式相位向 量中第r个元素；^:为第！/单源模式相位向量中第1个元素;为单源模式的下标;M为观测 信号的个数。其中，步骤103中的对单源模式进行基于DB指数的K-均值聚类，获取混合矩阵 的估计具体为：
[0067]将聚类数目初始化为1 = 2,利用K-均值聚类法将单源模式集聚成I类，并计算DB指 数；
[0068] 令I = 1+1，并再次聚类，求DB指数;直到寻找到I满足DKDh，且DKDI+1，则I可作为 源信号数目N的估计，而聚类中心所组成的矩阵，作为混合矩阵的估计。
[0069] 其中，步骤101中的获取合并后的频率序列具体为：
[0070] 对每个观测信号分别做频谱校正，获取若干个频率估计集合；
[0071 ]将频率估计集合进行混合，按照升序排列，获取合并后的频率序列。
[0072]进一步地，对每个观测信号分别做频谱校正，获取若干个频率估计集合具体为： [0073]对观测信号做加汉宁窗L点DFT变换，得到观测频谱；
[0074] 收集观测频谱的所有峰值的位置，计算峰值及其次峰的比值，获取中间参数：
[0075] 根据中间参数估计频率偏移:通过频率偏移获取频率估计、幅值估计及相位估计。
[0076] 综上所述，本方法可准确提取不同源信号的谐波特征参数(频率、幅值和相位），在 短样本所导致的FFT分辨率较低的情况下，频谱校正有效的抵御了谱泄漏和栅栏效应(为本 领域技术人员所公知的技术术语，本发明实施例对此不做赘述）的影响。同时，一个由频率 分辨、候选有效模式选择、单源频率判别准则的步骤构成的频率筛选流程可以排除谐波成 分间的相互干扰。更进一步，该方法通过引入改进的K-均值算法，可以对源信号的数目进行 估计，从而不需对源信号的数目的先验知识。
[0077] 实施例2
[0078]下面结合具体的数学公式、图2对实施例1中的方案进行详细的描述，并给出本方 法的详细操作流程，详见下文：
[0079] 一、盲源分离模型
[0080] 1.1、时域模型
[0081] 考虑N个未知源和M路观测。假定所研究的旋转机械结构具有高度的刚性，且尺寸 较小，因而其传播时延可以忽略[1]。在此情况下，旋转机械振动信号混合系统可以被当作即 时混合，其表达式可以表示为：
[0082] x(t) =As(t)+n(t) (1)
[0083] 其中，8(1：) = [81(1：)，'"，811(1：)，'"，841：)]1'为源信号矢量，叉（1：) = |^1(1：)，'"，叉111 (t)，…，XM(t)]T为观测信号矢量，11(〇 = [111(0，112(0，~，1?(0]7是噪声矢量4为混合矩 阵。短样本BSS所要完成的任务是:在不知道混合矩阵A的情况下，仅从少量样本的观测信号 Xl(t)~XM(t)中恢复源信号Sl(t)~SN(t)。
[0084] 根据N和M大小的相对关系，BSS问题可分为2个类型，超定或者正定情形(NSM)，及 欠定情形(N>M)。本发明实施例仅针对超定情形进行分析和描述。
[0085] 因为振动信号起源于机械部件的旋转，第n(n=l，2,…，N)个源信号~(〇可看作 一系列谐波信号的组合：
[0086] 以= + (2) 产1
[0087] 其中，？"是频率成分数，(^,1)，匕,1)，011,1)分别是第1 1个源信号，第?个频率成分的幅值、 频率、以及相位参数。
[0088]基于上述现有的时域模型，本发明实施例旨在提出一种BSS算法，仅利用少量观测 样本，估计混合矩阵A并恢复源信号S1(t)~sN(t)。同时，值得注意的是，在工业应用中，源信 号的数目往往是未知的，本发明实施例同时也对源信号的数目进行估计。
[0089] 1.2、基于谐波的BSS模型
[0090] 联立式（1)和（2)，可以发现，若观测信号Xm(t)可以进一步与谐波参数关联，则可 通过参数c n,p，fn,p，9n,p估计出混合矩阵A。
[0091] 因为实信号包含共辄的两边带谱线，将源信号Sn(t)重写为：
[0092] ~,(，) =元(〇 + ;:(，） (3)
[0093] 其中A(〇为信号的右边带所对应的时域信号;<(〇:为信号的左边带所对应的时域 信号，且满足如下的表达式：
(4)
[0095] 当某频率成分fn,P远离直流成分，仅考虑其中一个边带谱线的时域信号(此处仅考 虑信号右边带时域信号1(0)时，就可包含全部的谱信息。再与式（1)联立，可得BSS的频域 模型：
[0096] X(/) = AS(/) (5)
[0097] 其中，句/)为观测信号x(t)的频谱的右半边带JK)为源信号频谱的右半边带;f 为源信号及观测信号在频域的自变量。易知，复指数信号的傅立叶变换是一个冲激函数'8 (?)'。因此，第n(l<n<N)个源信号s n(t)的右半边带无⑴的傅立叶变换表达式为：
[0098]文(/)=疋|>"./(/ 心 (6) 产1
[0099]其中，S( ?)为冲激函数。
[0100]将混合矩阵A表示为列向量的形式[m，…，aN]，将式(6)代入式(5)中可得：
[0102]其中，之(/)为第m(l彡m彡M)个观测信号xm(t)的右半谱。
[0103]为了估计出混合矩阵A的各列ai~aN，考虑这样一种特殊的频率成分//，它仅仅出 现在一个源信号（不失一般性，假设出现在第n个源信号&(〇)，而其他源信号不含有该频 率成分//，即该频率成分4满足如下的要求： \fp- i{fw_p,P = l---Pw},n=l---,N,W^n
[0105] 其中，&为除了第n个源信号&⑴以外的、第泛个源信号中的第p个频率成分;Pn* 第n个源信号的频率成分数目；P n为第i个源信号的频率成分数目。
[0106] 则对于频率成分//'式(7)变为
[0108] 其中，为右半谱成分的谱峰心对应的频域向量；足4/丨.)(19KM)为第m个观 测信号xm(t)的右半谱的谱峰心的频谱值。从式(9)可以看出，与频率成分/；+对应的频域向 量、与混合矩阵向量 &"平行。因此，只要收集到充分多的单源模式./；*，即可估计混合矩 阵A的所有列向量。
[0109] 1.3、短样本盲分离的困难
[0110] 注意到式（7)是BSS的理想傅立叶变换模型，其中的频率f?是一个连续的变量。然 而，在实际中理想的傅立叶变换是不可实现的，因为这需要无穷多的样本。在实际应用当 中，理想傅立叶变换被L点(L是观测信号样点的数目）的离散傅立叶变换(Digital Fourier Transform，DFT)代替。在DFT中，频率f仅仅允许为频率分辨率A f( A f = fs/L，fs为采样速 率)的整数倍。
[0111 ]然而，第n个源信号的频率fn,P往往不会正好等于DFT频率分辨率A f = fs/L的整数 倍，这将导致式(7)的冲激函数S(f-fn,p)不能够获得理想的取样结果。这个偏差同样反映在 观测信号^⑴的DFT频谱之(/#)上，即出现谱泄漏现象。
[0112] 不失一般性，将第n个源信号~(〇的频率成分fn,P定义为频率分辨率A f的整数倍 kn,P和小数倍Sn,p的和，表达式如下所示：
[0113] fn>P = (kn>P+8n>P) Af,kn>PGZ+;8pG(-0.5,0.5] (10)
[0114] 其中，Z+为正整数。
[0115] 当样本长度L变小，DFT频率分辨率A f更大，DFT频谱变得更粗糙，即栅栏效应更严 重。由于栅栏效应的作用，式（10)被忽略的小数项'Sn, P A f'的值将会变得更大，这将导致 DFT谱峰的谱线位置(即DFT频率估计值，又.,,=)与实际频率的偏差会更大。
[0116] 为了克服这一难题，本发明实施例引入频谱校正以得到精确的谐波参数估计。
[0117] 二、基于频谱校正的BSS
[0118] 2.1、频谱校正
[0119] 为了去除短样本所带来的频谱泄漏和栅栏效应所导致的巨大偏差的影响，本发明 实施例引入了比值频谱校正法[17]。
[0120] 比值频谱校正法利用谱峰及其相邻最大的谱线对频偏进行估计，再利用频偏对频 率成分的频率、幅值和相位进行校正。对各路观测信号141〇， 111=1，2，-^进行频谱校正，每 个观测信号都可得到校正后的谐波参数，即：频率估计集、幅度估计集相位估 计集{U。
[0121] 得到校正的谐波参数后，对于一个频率成分，如式(8)和式(9)可知，仅当该频 率成分只属于一个源信号，不属于其他源信号时，其谐波参数可用来估计A的列向量。因此， 在估计混合矩阵之前，必须筛选出单源模式。
[0122] 2.2、单源成分筛选
[0123] 单源成分筛选分为三个步骤:频率合并、候选有效模式选择、单源模式判别准则。
[0124] 2.2.1、频率合并
[0125] 值得注意的是，受噪声、干扰及校正法本身精度的影响，即使是同一个频率成分， 其从所有观测中经过频谱校正估计的结果仍会存在轻微的不同。因此，需通过频率合并步 骤将同一个频率成分的不同估计进行合并，以得到更加精确的频率估计，并且组合出每个 频率成分的幅值向量和相位向量。若将M个观测的频率估计集合合并为一个大集合，并按升 序排列，则前面所提的同一个频率成分的不同估计由于仅存在很小的差异，必定倾向于形 成一个簇。假定总共有P个簇形成(其中每个簇对应一个合并的频率成分），不失一般性，将 第p个簇的元素定义为tt, j =】工…，( r p是该簇的元素数目）。则该r 元素应由其均 值7；来代替，以达到合并的效果，其均值表达式为：
(ID
[0127] 2.2.2、候选有效模式选择
[0128] 理论上，只要混合矩阵A的元素中不含零，从式（1)中可以看出，所有源信号的所有 频率成分应同时包含在M路观测信号中。因而，那些不被所有观测信号包含的频率成分往往 为伪频率成分，很可能由噪声或频率成分间的干扰产生。由大量实验可发现，这种伪频率成 分往往幅值很小，且不与混合矩阵A的列向量平行。为了排除伪频率成分，必须选择(在一定 误差范围内）被所有观测信号包含的频率作为后续步骤的候选模式。
[0129] 在实际操作时，假定一小的阈值￡>0,对于经过合并后的某频率成分叉，若对每个 观测信号^⑴的频率估计集|.H U =U,…，M )，存在且仅存在一个下标瓦，满足
[0130] \f"uAi-fr\<￡ (12)
[0131] 则某频率成分乙可被当作候选有效频率。同时，根据式(9)，可利用候选模式在M个 观测信号内的幅值及相位参数构建如下式所示的候选有效模式向量巧 ：
[0132] ip= (13) _"A/，内广 _
[0133] 其中，戶为候选有效模式的数目。经过候选有效频率的选择后，合并频率的数目下 降到歹个。
[0134] 2.2.3、单源模式判别准则
[0135] 在旋转机械故障分析中，常常会发生多个源信号同时包含同一个频率成分（即重 叠频率）的情况。显然，这种频率成分不能满足式(8)，从而，无法得出式(9)的结论，即该种 频率成分将不与混合矩阵A的列向量平行，无法用来估计混合矩阵A。因而，重叠频率应被去 除，以保证混合矩阵A的估计精度。
[0136] 假定在候选有效频率中，/#是单源的频率成分，即频率成分满足式(8)，由式(9) 可以看出，M/j的所有项的相位均来自于同一个频率成分//的相位，因而M/j各项的相 位必定具有一致性，即在一定误差范围内均相等，且频率成分.//对应的模式向量^组成的 模式集#/，￥= 1，2，~，1为单源模式集"为单源模式的下标'为单源模式的数目）。
[0137] 由以上分析可以总结出单源的频率成分的两个特性：
[0138] 1)单源的频率成分的幅值向量与混合矩阵A的某个列向量平行。
[0139] 2)单源的频率成分的相位向量具有相位一致性，即相位向量的元素满足下式：
(14)
[0141 ] 其中，1彡r，1彡M，r乒1，D0为小阈值，C为排列组合运算。
[0142] 2.3、基于DB指数的源数及混合矩阵估计
[0143] 若源信号的数目已知，则可以利用K-均值聚类法对上节所得单源模式集 / =1，…，#1进行聚类来估计混合矩阵A。然而在实际工程应用中，源信号的数目'N'往往是未 知的。为了在源信号数目未知的情形下估计混合矩阵A，并且恢复源信号，本发明实施例提 出了基于DB指数[18]和1(-均值聚类法的聚类技术，在未知源信号数目前提下，寻找聚类模式 的自然划分。
[0144] 首先将聚类数目设为I，然后对单源模式集进行K-均值聚类法，将 其聚成I簇(^匕为聚类所得的第i(i = l，2,…，I)个簇，则
[0145] C； = ' ri = 1,2,? ? | (15)
[0146] 其中，n为聚类后的单源模式的下标，I；为第i个簇的ri个单源模式，Ri为第i个簇 的单源模式个数。(^满足如下表达式：
[0147] 1丨(、,=[z, =!.丄..../) '丨 （16)
[0148] (17) 2=1
[0149] 其中，U为求并集操作。
[0150] 对于某一特定的聚类结果(聚类数目为I)，DB指数可以评价该结果的合适程度，其 衡量方法独立于聚类算法，因而可以作为一个通用的衡量指标。DB指数可以定义为：
(IS)
[0152]其中，0:为聚类数目等于I时的聚类结果的DB指数;代表第i簇和第j簇之间的 相似度指标;Gi及Gj代表第i簇及第j簇之间的离散度指标。Gi的表达式为
(19)
[0154] 其中，ci为第i簇的聚类中心，Gj的表达式同理，本发明实施例对此不做赘述。
[0155] 其中，Mi，j的表达式为
[0156] Mi, j = | | Ci-Cj | (20)
[0157] 其中，Cj为第j簇的聚类中心。
[0158] 通过上述的论述，易知，同一簇内的元素，离散度越小，代表聚集的程度越高，则聚 类效果越好。不同簇之间的距离越大，代表相似度越低，则聚类效果越好。因而，针对DB指数 的表达式，可以得出，DB指数最小所对应的聚类的结果性能最好。由此可知，源信号数目的 数估计可通过对数据进行多次聚类，求其DB指数，其中DB指数的极小值对应的聚类的簇数 既可作为源信号数目的估计，而该聚类的各簇的聚类中心所组成的矩阵A，则可作为混合 矩阵A的估计。
[0159] 2.4、源信号的恢复
[0160] 针对超定情形，在得到混合矩阵的估计1后，源信号的恢复可以通过以下表达式 得到
[0161] s(f) = A^x.(f) (21)
[0162] 其中，表示求混合矩阵A的伪逆操作。
[0163] 3、本发明实施例提供的基于频谱校正的短样本旋转机械信号盲分离的详细流程 如下，下面结合图2进行描述：
[0164] 3.1频谱校正步骤
[0165] Step 1:对观测信号Xm(t)，m=l，…，M做加汉宁窗L点DFT变换，得到观测频谱Xm (k),k = 0,l,---,L；
[0166] Step 2:收集观测频谱Xm(k)的所有峰值的位置。对谱峰，仏(Pm是第m个 观测的谱峰数目），计算峰值尤)及其次峰的比值~a .
[0168] 其中，总为谱峰左侧的谱线；+1)为谱峰右侧的谱线。
[0169] 则可得到中间参数~1
[0170] ?mA=(2-^J/(l + ^J： (23)
[0171] Step 3:根据中间参数~a估计频率偏移4,$ : r ,, 卜w若丨碌味~+1)| 削
[0172] 4而= 丨. 丨丨， 丨 (24) .卜％，若|碌中未，1)|
[0173] 则频率估计也，?为：
[0174] +4,^)2^ (25)
[0175] Step 4:幅值估计及相位估计九.~则分别为：
[0176] 4^ =2^,a, (26)
[0177] 之a =3_[尤取,,J] 一喊 i(￡-l)/￡ (27)
[0178] 其中，angle( ?)代表取相角运算。经过校正，得到频率估计集合{.元、幅值估计 集合、以及相位估计集合化M丨w =
[0179] 具体实现时，上述对频谱进行校正的步骤为本领域技术人员所公知，本发明实施 例对此不做赘述。
[0180] 3.2、单源模式筛选步骤
[0181] Step 5:将M个频率集合进行混合，按照升序排列，得到P个频率簇， 其中第p(p = l，…，P)个簇由差别较小的频率估计结果组成^，0=1，2，-,1^:，1\为该簇的 元素个数，则合并后的频率为式（11)所示，收集所有合并频率，则可得到合并后的频率序列 {Jp'P = c
[0182] Step 6:找出合并后的频率序列试，# = 1，…，4中满足给定一个小阈值￡>0,对于每 个m都存在且仅存在一个下标瓦w满足式（12)的频率，则该频率为候选有效频率，对应一个 候选有效模式％，如式(13)所示，其中P为候选有效模式数目。
[0183] Step 7:将候选有效模式集合K)作为输入，按照式（⑷识别中的单源模式模 式并记为^，/=1, ????，}，即满足上式(14)的即为单源模式。
[0184] 3.3基于DB指数的K-均值聚类及源信号的恢复步骤
[0185] Step 8:将聚类数目初始化为1 = 2,利用K-均值聚类法将单源模式集聚成I类， 并按照式（18)~式(20)计算其对应的DB指数D:。
[0186] Step 9:令1 = 1+1，并再次聚类，求DB指数Didtep 10:重复Step 9直到寻找到I满 SDKDpi，且口工③^(若I = 2，只需满SDKDm )，贝丨JI可作为源信号数目N的估计，而聚类中心 所组成的矩阵，可以作为混合矩阵的估计A。
[0187] Step 11:源信号的恢复可以通过表达式(21)得到。
[0188] 综上所述，本方法可准确提取不同源信号的谐波参数，有效的抵御了谱泄漏和栅 栏效应；同时，本方法提供的一个由频率分辨、候选有效模式选择、单源频率判别准则的步 骤构成的频率筛选流程可以排除谐波成分间的相互干扰;更进一步，该方法通过引入改进 的K-均值算法，可以对源信号的数目进行估计，从而不需对源信号的数目的先验知识。
[0189] 实施例3
[0190] 下面结合具体的计算结果、附图对实施例1-2中的方法进行可行性验证，详见下文 描述：
[0191] 四、实验
[0192] 4.1、合成信号仿真
[0193] 考虑3X2混合系统 '0.59 0.89'
[0194] A= 0.97 0.47 (28； 0.59 0.53
[0195] 源信号表达式为 i', {t) = cos (2^-100?+ tcJ 9) + l,4cos(2w200? + 2yr / 9) + I.85gos(2^400/ + 7T / 3) .(29)
[0196] (r) = 1.7cos(2^r50f + ^/ 36) + 0.8cos(2^80r + 2^- /9) +1,2cos(2^800f + 5ni 12)
[0197] 采样速率匕=2000泡，考虑4种样本长度(1 = 400,200，100,70)。引入€&8卜1〇八算 法作为比较算法，并引入恢复信号与源信号的相关系数，即恢复相关系数作为性能衡量指 标。因为fast-ICA算法基于信号的统计特性，所以每次运行结果有差异。为了能更好的测得 其性能，将实验重复1〇〇〇次，记录其中成功次数，和成功的实验的恢复相关系数的平均值。
[0198] 其中一次的运行结果展示在图3~图6中，1000次实验的成功次数，及成功试验平 均恢复相关系数如表1。
[0199] 表1合成信号成功次数及平均恢复相关系数
[0201] 如图3和图4所示，在长样本(L = 400)的情形下，fast-ICA和本发明都能够精确的 恢复源信号。然而，当样本数目减少到L = 70,即短样本的情形下，由图5和图6可以得出， fast-ICA的恢复信号出现了明显的失真，而本发明仍然精确的恢复了源信号。
[0202] 从表1中可知，随着样本的数目减少，fast-ICA的成功恢复次数急剧下降，且成功 实验的平均恢复相关系数也轻微下降，而本方法的性能始终保持在相关系数为1的理想情 形，且所有实验都成功的进行。因而，在处理合成信号方面，本方法在短样本情况下，其性能 明显优于fast-ICA。
[0203] 4.2、实际振动信号
[0204] 实验
[0205] 仍然利用上述的混合系统，而信号则是两个实际的机械故障信号。Sl(t)为工频 89.6853Hz的不平衡故障信号，s 2(t)为工频102.8811Hz转子不对中故障信号。考虑不同的 样本长度(L = 400，200，100，70 )，同样在样本长度下进行1000次重复试验，其中一次的运行 结果展示在图7~图10中，1000次实验的成功次数，及成功试验平均恢复相关系数如表2。
[0206] 表2实际振动信号成功次数及平均恢复相关系数
[0208] 从图7和图8中可知，在长样本(L = 400)情形下，fas t-ICA和本方法都可精确恢复 源信号。而对于短样本情形(L = 70)，图9中可以清晰地看出，fast-ICA的恢复信号存在明显 的失真，而本方法仍然能够精确恢复源信号，见图10。从表2中可以看出，随着样本数目从 400降到70,本方法的性能均强于fas t-ICA。更进一步，当样本数目减小到一定程度时， fast-ICA的成功次数出现急剧下降的情况，而本方法性能较为稳定。在处理实际振动信号 方面，本方法在短样本情况下，其性能仍然明显优于fast-ICA。
[0209] 实施例4
[0210] 本发明实施例提供了一种与实施例1至2中的方案对应的装置，详见下文描述:硬 件实施图如图11所示，将从振动传感器采集到的多路观测信号x(t)经过A/D(模数转化器） 采样得到样本序列x(n)，以并行数字输入的形式进入DSP芯片，经过DSP芯片的内部算法处 理，得到混合矩阵的估计;最后借助输出驱动及其显示模块显示恢复的源信号。
[0211] 其中，图10的DSP(Digital Signal Processor，数字信号处理器）为核心器件，在 信号参数估计过程中，完成如下主要功能：
[0212] (1)调用核心算法，完成信号的FFT变换，频谱校正，单源成分筛选(包含频率合并、 候选成分筛选和单源成分识别)得到信号的单源模式集；
[0213] ⑵利用基于DB指数的改进K-均值聚类法对单源模式集进行聚类，得到源数目及 混合矩阵的估计，并通过观测信号和混合矩阵估计的伪逆求得源信号的估计，并将结果输 出至驱动和显示模块；
[0214] DSP器件的内部程序流程如图12所示。本发明将所提出的算法植入DSP器件内，基 于此完成高精度、低复杂度、高效的旋转机械振动信号盲分离。
[0215] 图12流程分为如下几个步骤：
[0216] 1)首先需根据旋转机械故障诊断的实时性要求，设置信号的采样点数L;
[0217] 2)其次，CPU主控器从I/O端口读采样数据，进入内部RAM;
[0218] 3)最终，按图1和图2本发明的处理过程进行源信号的恢复，并将其通过外部显示 装置进行显示。
[0219] 本发明实施例对各器件的型号除做特殊说明的以外，其他器件的型号不做限制， 只要能完成上述功能的器件均可。
[0220] 参考文献
[0221] [l]Gelle G,Colas M,Serviere C.Blind source separation:a tool for rotating machine monitoring by vibrations analysis?[J].Journal of sound and vibration,2001,248(5):865-85.
[0222] [2]Sadhu A，Hazra B.A Novel Damage Detection Algorithm using Time-Series Analysis-Based Blind Source Separation[J].Shock and Vibration,2013,20 (3):423-38.
[0223] [3]Li Z，Yan X,Tian Z，et al.Blind vibration component separation and nonlinear feature extraction applied to the nonstationary vibration signals for the gearbox multi-fault diagnosis[J].Measurement,2013,46(1):259-71.
[0224] [4]Cui L,ffu C，Ma C，et al.Diagnosis of roller bearings compound fault using underdetermined blind source separation algorithm based on null-Space pursuit[J].Shock and Vibration,2015.
[0225] [5]Bao G，Ye Z,Xu X,et al.A compressed sensing approach to blind separation of speech mixture based on a two-layer sparsity model[J].IEEE Transactions on Audio Speech and Language Processing，2013,21(5):899_906.
[0226] [6]Sha Z~C,Huang Z_T，Zhou Y_Y，et al.Frequency-hopping signals sorting based on underdetermined blind source separation[J].IET Communications，2013，7 (14):1456-64.
[0227] [7]Lei Y，Lin J，He Z,et al.A review on empirical mode decomposition in fault diagnosis of rotating machinery[J].Mechanical Systems and Signal Processing,2013,35(1):108-26.
[0228] [8]Serviere C,Fabry P.Blind source separation of noisy harmonic signals for rotating machine diagnosis[J]?Journal of Sound and Vibration, 2004,272(l):317-39.
[0229] [9]Lees A.Misalignment in rigidly coupled rotors[J].Journal of Sound and Vibration,2007,305(1):261-71.
[0230] [10]Chen Y,Du R,Qu L.Fault features of large rotating machinery and diagnosis using sensor fusion[J].Journal of Sound and Vibration，1995，188(2): 227-42.
[0231] [11]Hyvrinen A，0ja E.A fast fixed-point algorithm for independent component analysis[J].Neural computation,1997,9(7):1483-92.
[0232] [12]Mcneill S I,Zimmerman D C.Relating independent components to free-vibration modal responses[J]?Shock and Vibration，2010，17(2):161-70?
[0233] [13]Zuo M J，Lin J，Fan X.Feature separation using ICA for a onedimensional time series and its application in fault detection[J].Journal of Sound and Vibration，2005,287(3):614_24.
[0234] [14]Ypma A,Pajunen P.Rotating machine vibration analysis with second-order independent component analysis proceedings of the Proceedings of the First International Workshop on Independent Component Analysis and Signal Separation,ICA,F,1999[C].Citeseer.
[0235] [15]Roan M，Erling J，Sibul L.A new,non-linear,adaptive,blind source separation approach to gear tooth failure detection and analysis[J] .Mechanical Systems and Signal Processing，2002，16(5):719-40?
[0236] [16]De Lathauwer L，Castaing J，Cardoso J-F.Fourth-order cumulant-based blind identification of underdetermined mixtures[J]?IEEE Transactions on Signal Processing，2007,55(6):2965-73.
[0237] [17]Agrez D.Improving phase estimation with leakage minimization[J] .IEEE Transactions on Instrumentation and Measurement，2005，54(4):1347-53?
[0238] [18]Davies D L,Bouldin D ff.A cluster separation measure[J].IEEE Transactions on Pattern Analysis and Machine Intelligence,1979,(2):224-7.
[0239] 本领域技术人员可以理解附图只是一个优选实施例的示意图，上述本发明实施例 序号仅仅为了描述，不代表实施例的优劣。
[0240] 以上所述仅为本发明的较佳实施例，并不用以限制本发明，凡在本发明的精神和 原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。
【主权项】
1. 一种基于短样本频谱校正的旋转机械振动信号盲分离方法，其特征在于，所述方法 包括W下步骤： 获取合并后的频率序列中满足给定一个阔值，都存在且仅存在一个下标满足第一表达 式的频率，则该频率对应一个候选有效模式； 将候选有效模式组成的集合作为输入，识别候选有效模式集合中的单源模式； 对单源模式进行基于DB指数的K-均值聚类，获取混合矩阵的估计； 通过混合矩阵的估计、观测信号恢复源信号。2. 根据权利要求1所述的一种基于短样本频谱校正的旋转机械振动信号盲分离方法， 其特征在于，所述第一表达式具体为：其中，/5。,为频率估计集;又为频率;e为阔值。3. 根据权利要求1所述的一种基于短样本频谱校正的旋转机械振动信号盲分离方法， 其特征在于，所述将候选有效模式组成的集合作为输入，识别候选有效模式集合中的单源 模式具体为： 1) 单源模式的幅值向量与混合矩阵的某个列向量平行； 2) 单源模式的相位向量具有相位一致性，即相位向量的元素满足下式：其中，1《'，1《1^声1成〇为阔值;(：为排列组合运算;<^,,,;为第1/单源模式相位向量中 第r个元素;为第l/单源模式相位向量中第1个元素单源模式的下标;M为观测信号 的个数。4. 根据权利要求1所述的一种基于短样本频谱校正的旋转机械振动信号盲分离方法， 其特征在于，所述对单源模式进行基于DB指数的K-均值聚类，获取混合矩阵的估计具体为： 将聚类数目初始化为I = 2，利用K-均值聚类法将单源模式集聚成I类，并计算DB指数； 令I = 1+1，并再次聚类，求DB指数;直到寻找到I满足Di<Di-i，且Di<Di+i，则I可作为源信 号数目N的估计，而聚类中屯、所组成的矩阵，作为混合矩阵的估计。5. 根据权利要求1所述的一种基于短样本频谱校正的旋转机械振动信号盲分离方法， 其特征在于，所述获取合并后的频率序列具体为： 对每个观测信号分别做频谱校正，获取若干个频率估计集合； 将频率估计集合进行混合，按照升序排列，获取合并后的频率序列。6. 根据权利要求5所述的一种基于短样本频谱校正的旋转机械振动信号盲分离方法， 其特征在于，所述对每个观测信号分别做频谱校正，获取若干个频率估计集合具体为： 对观测信号做加汉宁窗L点DFT变换，得到观测频谱； 收集观测频谱的所有峰值的位置，计算峰值及其次峰的比值，获取中间参数： 根据所述中间参数估计频率偏移:通过所述频率偏移获取频率估计、幅值估计及相位 估计。7. -种用于实施权利要求1-6中任一权利要求所述的一种基于短样本频谱校正的旋转 机械振动信号盲分离方法的装置，其特征在于， 所述装置包括:振动传感器、模数转化器、微控制器和输出驱动及其显示模块， 所述振动传感器采集到的多路观测信号X(t)经过所述模数转化器采样得到样本序列， W并行数字输入的形式进入所述DSP忍片，经过所述DSP忍片的内部算法处理，得到混合矩 阵的估计;最后借助所述输出驱动及其显示模块显示恢复的源信号。8.根据权利要求7所述的装置，其特征在于，所述微控制器具体为DSP忍片。
【文档编号】G01H17/00GK105910701SQ201610240265
【公开日】2016年8月31日
【申请日】2016年4月15日
【发明人】黄翔东, 靳旭康, 马欣
【申请人】天津大学