说明书一种动态信号分析方法及装置
技术领域
本发明涉及信号分析技术领域，更具体的说，是涉及一种动态信号分析方法及装置。
背景技术
在旋转机械故障诊断以及生物医学信号处理等领域，一般通过对系统产生的动态信号进行分析，进而对系统的健康状况进行评估。
由于系统产生的动态信号常常具有非平稳特性，在对此类信号进行分析时，单独的时域分析方法以及单独的频域分析方法，都不能提供足够的非平稳信息来对系统的健康状况进行有效评估。
在现有技术中，一般采用传统的短时傅里叶变换、连续小波变换或者维格纳—威尔分布等时频分析方法，对系统产生的动态信号进行分析。通过上述任意一种时频分析方法，都能同时反映系统产生的动态信号在时域和频域上丰富的非平稳信息，用以对系统的健康状况进行有效评估。
但是，一方面由于系统本身的复杂性，其产生的动态信号又表现为非线性。而传统的时频分析方法在对系统产生的动态信号进行分析的过程中，并未考虑动态信号的非线性。另一方面，由于各种原因，系统产生的动态信号往往含有大量的噪声。而传统的时频分析方法在对系统产生的动态信号进行分析的过程中，并未考虑噪声的影响。
由上可知，现有的时频分析方法，并不能有效分析系统产生的动态信号，进而对系统的健康状况的评估造成了不良影响。
发明内容
有鉴于此，本发明提供了一种动态信号分析方法及装置，以克服现有技术中由于采用传统的时频分析方法分析系统产生的动态信号时，不能考虑动态信号的非线性以及噪声影响，对系统的健康状况的评估造成不良影响的问题。
为实现上述目的，本发明提供如下技术方案：
一种动态信号分析方法，包括：
对接收到的动态信号进行时频分析，生成与所述动态信号相对应的具有预设嵌入维数的高维时频分析信号；
对所述时频分析信号进行流形学习，生成具有预设本征维数的低维流形结构信号，所述预设本征维数小于所述预设嵌入维数；
对所述流形结构信号进行转换，生成包含多维数据的时频流形结构信号，所述时频流形结构信号中的第一维数据为所述动态信号的时频流形特征数据。
优选地，所述对接收到的动态信号进行时频分析，生成与所述动态信号相对应的高维时频分析信号的过程包括：
接收所述动态信号；
根据预设嵌入维数，将所述动态信号在相空间中重构，生成具有所述预设嵌入维数的第一信号，所述预设嵌入维数至少为12维；
根据预设时频分析方法，对所述第一信号整个频带范围内的所有时间点数进行时频分析，生成具有所述预设嵌入维数的第二信号，所述预设时频分析方法包括短时傅里叶变换、连续小波变换或维格纳—威尔分布变换；
将所述第二信号进行转换，生成具有所述预设嵌入维数的所述高维时频分析信号。
优选地，所述对所述时频分析信号进行流形学习，生成具有所述预设本征维数的低维流形结构信号的过程包括：
根据局部线性嵌入算法、等距映射算法或局部切空间排列算法对所述时频分析信号进行流形学习，生成具有所述预设本征维数的低维流形结构信号。
优选地，当所述预设时频分析方法为短时傅里叶变换时，所述根据预设时频分析方法，对所述第一信号整个频带范围内的时间点数进行时频分析，生成具有所述预设嵌入维数的第二信号的过程还包括：
在所述第一信号整个频带范围内通过增加短时傅里叶变换窗函数的移动步长减少时间点数，并对所述第一信号整个频带范围内的部分时间点数进行时频分析，生成具有所述预设嵌入维数的第二信号。
优选地，当所述预设时频分析方法为短时傅里叶变换、连续小波变换方法或维格纳—威尔分布变换方法时，所述根据预设时频分析方法，对所述第一信号整个频带范围内的时间点数进行时频分析，生成具有所述预设嵌入维数的第二信号的过程还包括：
在所述接收到的动态信号整个频带范围内选取故障频带；
对所述第一信号故障频带范围内的所有时间点数进行时频分析，生成具有所述预设嵌入维数的第二信号。
优选地，根据局部切空间排列算法对所述高维时频分析信号进行流形学习，生成具有所述预设本征维数的低维流形结构信号的过程包括：
计算所述时频分析信号的局部信息；
根据所述局部信息计算所述时频分析信号的排列矩阵；
根据所述排列矩阵计算具有所述预设本征维数的全局坐标，所述全局坐标即为具有所述预设本征维数的低维流形结构信号。
优选地，所述计算所述时频分析信号的局部信息的过程包括：
在所述时频分析信号中选取待处理元素；
选取与所述待处理元素最近的预设个数的元素，所述预设个数大于所述预设嵌入维数；
将所述预设个数的元素进行中心化处理得到中心化矩阵；
计算所述中心化矩阵的右奇异向量，并按照右奇异值从大到小的顺序，选取所述预设本征维数个数右奇异值相对应的右奇异向量，作为所述时频分析信号的局部信息。
优选地，所述根据所述局部信息计算所述时频分析信号的排列矩阵的过程包括：
根据所述时频分析信号以及所述中心化矩阵计算得到第一矩阵；
根据所述时频分析信号的局部信息计算得到第二矩阵；
根据第一矩阵以及第二矩阵计算所述时频分析信号的排列矩阵。
一种动态信号分析装置，包括：
时频分析单元：用于对接收到的动态信号进行时频分析，生成与所述动态信号相对应的具有预设嵌入维数的高维时频分析信号；
流形学习单元：用于对所述高维时频分析信号进行流形学习，生成具有预设本征维数的低维流形结构信号，所述预设本征维数小于所述预设嵌入维数；
流形转换单元：用于对所述流形结构信号进行转换，生成包含多维数据的时频流形结构信号，所述时频流形结构信号中的第一维数据为所述动态信号的时频流形特征数据。
优选地，所述时频分析单元包括：
接收子单元：用于接收所述动态信号；
相空间重构子单元：用于根据预设嵌入维数，将所述动态信号在相空间中重构，生成具有所述预设嵌入维数的第一信号，所述预设嵌入维数至少为12维；
时频分析子单元：用于根据预设时频分析方法，对所述第一信号整个频带范围内的所有时间点数进行时频分析，生成具有所述预设嵌入维数的第二信号，所述预设时频分析方法包括短时傅里叶变换、连续小波变换或维格纳—威尔分布变换；
转换子单元：用于将所述第二信号进行转换，生成具有所述预设嵌入维数的所述时频分析信号。
优选地，所述流形学习单元包括：
局部信息计算子单元：用于计算所述时频分析信号的局部信息；
排列矩阵计算子单元：用于根据所述局部信息计算所述时频分析信号的排列矩阵；
全局坐标计算子单元：用于根据所述排列矩阵计算具有所述预设本征维数的全局坐标，所述全局坐标即为具有所述预设本征维数的低维流形结构信号。
优选地，当所述预设时频分析方法为短时傅里叶变换时，所述时频分析子单元包括：
第一计算模块：在所述第一信号整个频带范围内通过增加短时傅里叶变换窗函数的移动步长减少时间点数；
第二计算模块：对所述第一信号整个频带范围内的部分时间点数进行时频分析，生成具有所述预设嵌入维数的第二信号。
优选地，当所述预设时频分析方法为短时傅里叶变换、连续小波变换方法或维格纳—威尔分布变换方法时，所述时频分析子单元包括：
故障频带选择模块：用于在所述接收到的动态信号整个频带范围内选取故障频带；
第三计算模块：用于对所述第一信号故障频带范围内的所有时间点数进行时频分析，生成具有所述预设嵌入维数的第二信号。
优选地，所述局部信息计算子单元包括：
待处理元素选取模块：用于在所述时频分析信号中选取待处理元素；
元素选取模块：用于选取与所述待处理元素最近的预设个数的元素，所述预设个数大于所述预设嵌入维数；
中心化处理模块：用于将所述预设个数的元素进行中心化处理得到中心化矩阵；
局部信息计算模块：用于计算所述中心化矩阵的右奇异向量，并按照右奇异值从大到小的顺序，选取所述预设本征维数个数右奇异值相对应的右奇异向量，作为所述时频分析信号的局部信息。
优选地，所述排列矩阵计算子单元包括：
第一矩阵计算模块：用于根据所述时频分析信号以及所述中心化矩阵计算得到第一矩阵；
第二矩阵计算模块：用于根据所述时频分析信号的局部信息计算得到第二矩阵；
第三矩阵计算模块：用于根据第一矩阵以及第二矩阵计算所述时频分析信号的排列矩阵。
经由上述的技术方案可知，与现有技术相比，本发明公开了一种动态信号分析方法及装置。首先，对接收到的动态信号进行时频分析，生成与动态信号相对应的具有预设嵌入维数的高维时频分析信号；然后，对时频分析信号进行流形学习，生成具有预设本征维数的低维流形结构信号，预设本征维数小于所述预设嵌入维数；最后，对流形结构信号进行转换，生成包含多维数据的时频流形结构信号，时频流形结构信号中的第一维数据为动态信号的时频流形特征数据。通过对接收的动态信号进行时频分析、流形学习以及转换，能够充分考虑动态信号的非线性以及信号中噪声的影响，使得最终的时频流形特征数据能够作为评估系统的健康状况的准确依据。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据提供的附图获得其他的附图。
图1为本发明实施例一公开的一种动态信号分析方法流程图；
图2为本发明实施例二公开的对动态信号时频分析的方法流程图；
图3为本发明实施例二公开的局部切空间排列算法流程图；
图4为本发明实施例三示例一对无故障阶段的齿轮健康信号分析图；
图5为本发明实施例三示例一对有严重磨损故障的齿轮振动信号分析图；
图6为本发明实施例三示例二对有外圈故障的列车轴承声音信号的分析图；
图7为本发明实施例三示例二对有内圈故障的列车轴承声音信号的分析图；
图8为本发明实施例三示例二对有滚子故障的列车轴承声音信号的分析图；
图9为本发明实施例三示例三对球轴承外圈故障的特征分析图；
图10为本发明实施例三示例三对球轴承滚珠故障的特征分析图；
图11为本发明实施例三示例三对球轴承内圈故障的特征分析图；
图12为本发明实施例四公开的动态信号分析装置结构示意图；
图13为本发明实施例五公开的时频分析单元结构示意图；
图14为本发明实施例五公开的流形学习单元结构示意图；
图15为本发明实施例五公开的时频分析子单元结构示意图；
图16为本发明实施例五公开的另一种时频分析子单元结构示意图；
图17为本发明实施例五公开的局部信息计算子单元结构示意图；
图18为本发明实施例五公开的排列矩阵计算子单元结构示意图。
具体实施方式
为了引用和清楚起见，下文中使用的技术名词的说明、简写或缩写总结如下：
流形学习：是一种非线性数据降维及信息挖掘方法，可用于系统状态特征高维数据的非线性特征挖掘，尤其是在旋转机械故障诊断中。高维系统状态参数可以是相空间重构的时域信号，也可以是时域和频域统计特征参数。基于流形学习的流形特征经过研究验证具有更佳的非线性效果，表明了流形学习在故障诊断中挖掘非线性结构的优势。但是，这些研究方法中高维数据本身携带的状态信息主要考虑了时域或频域特征等平稳信息，因而其流形特征是在平稳信息基础上提取的，也可以说目前主要使用平稳流形来表征系统状态，缺乏对反映系统非平稳状态的非平稳信息的考虑。
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
由背景技术可知，现有的时频分析方法，并不能有效分析系统产生的动态信号，进而对系统的健康状况的评估造成了不良影响。
因此，本发明公开了一种动态信号分析方法及装置。首先，对接收到的动态信号进行时频分析，生成与动态信号相对应的具有预设嵌入维数的高维时频分析信号；然后，对时频分析信号进行流形学习，生成具有预设本征维数的低维流形结构信号，预设本征维数小于所述预设嵌入维数；最后，对流形结构信号进行转换，生成包含多维数据的时频流形结构信号，时频流形结构信号中的第一维数据为动态信号的时频流形特征数据。通过对接收的动态信号进行时频分析、流形学习以及转换，能够充分考虑动态信号的非线性以及信号中噪声的影响，使得最终的时频流形特征数据能够作为评估系统的健康状况的准确依据。
上述动态信号分析方法的具体执行过程以及动态信号分析装置的具体构成将通过以下实施例进行详细说明。
实施例一
请参阅附图1，为本发明实施例一公开的一种动态信号分析方法流程图。该动态信号分析方法具体包括：
步骤101：对接收到的动态信号进行时频分析，生成与所述动态信号相对应的具有预设嵌入维数的高维时频分析信号；
首先，利用接收到的动态信号重构高维相空间，以充分显示嵌入在信号中的内在动态流形；需要说明的是，此处所说的高维相空间的维数为预设嵌入维数，理论上说，维数越大，计算精度也越大，但是相应的计算量也会增大，因此，本实施例中，预设嵌入维数可根据具体的实际应用进行设置。
其次，计算每一维信号的时频分布，以表示相空间中的非平稳信息，将每一维信号的时频分布的矩阵转化为列向量，组成高维时频分析信号，该时频分析信号与接收到的动态信号相对应，且维数为预设嵌入维数，能够充分表示原始动态信号的非平稳信息。
步骤102：对所述高维时频分析信号进行流形学习，生成具有预设本征嵌入维数的流形结构信号，所述预设本征维数小于所述预设嵌入维数；
接收的动态信号的流形结构揭示了系统内在的动态信息，它只存在于高维相空间中的低维子空间，可看成是动态信号的骨架。而动态信号的干扰噪声存在于相空间中的所有维数中，且是一种随机信息。通过对所述高维时频分析信号进行流形学习，得到所述低维流形结构信号，能够保留所述动态信号的骨架，去除随机白噪声。
步骤103：对所述流形结构信号进行转换，生成包含多维数据的时频流形结构信号，所述时频流形结构信号中的第一维数据为所述动态信号的时频流形特征数据。
需要说明的是，这里所说的转换只是对计算而得的流形结构信号表现形式的一种转换，是为了能够方便分析信号，进而确定系统的健康情况。
本发明实施例一公开的动态信号分析方法，首先，对接收到的动态信号进行时频分析，生成与动态信号相对应的具有预设嵌入维数的高维时频分析信号；然后，对时频分析信号进行流形学习，生成具有预设本征维数的低维流形结构信号，预设本征维数小于预设嵌入维数；最后，对流形结构信号进行转换，生成包含多维数据的时频流形结构信号，时频流形结构信号中的第一维数据为动态信号的时频流形特征数据。通过对接收的动态信号进行时频分析、流形学习以及转换，能够充分考虑动态信号的非线性以及信号中噪声的影响，使得最终的时频流形特征数据能够作为评估系统健康状况的准确依据。
在实施例一的基础上，本发明还公开了实现动态信号分析的具体实现过程，下面将通过以下实施例详细说明。
实施例二
本实施例为对实施例一中公开的一种动态信号分析方法的进一步描述。该方法具体包括：
步骤101：对接收到的动态信号进行时频分析，生成与所述动态信号相对应的具有预设嵌入维数的高维时频分析信号。
请参阅附图2，步骤101具体可分为以下几步：
步骤201：接收所述动态信号。
接收测得的非平稳动态信号，即时间序列x(t)=[x1,x2,...,xN]。需要说明的是，该动态信号为一维信号。
步骤202：根据预设嵌入维数，将所述动态信号在相空间中重构，生成具有所述预设嵌入维数的第一信号，所述预设嵌入维数至少为12维。
将一维信号x(t)重构在多维相空间中。公式为：
上式中，m是嵌入维数，τ是时间延迟，为相空间中具有所述预设嵌入维数的第一信号的第i个数据。为了使时频流形具有较高的时间分辨率，这里取τ=1。记n=N‑m+1，则上式可以表示为如下形式：
其中，它表示相空间中的第j维时间序列。嵌入维数m应足够大，一般可取m≥12。
步骤203：根据预设时频分析方法，对所述第一信号整个频带范围内的所有时间点数进行时频分析，生成具有所述预设嵌入维数的第二信号，所述预设时频分析方法包括短时傅里叶变换、连续小波变换或维格纳—威尔分布变换。
对相空间中具有所述预设嵌入维数的第一信号的每一维时间序列分别进行时频分析，生成具有所述预设嵌入维数的第二信号，即获取m维时频分布。
时频分布TFD(t,f)表示以时间t和频率f为变量的二维联合分布函数，利用它可以求某一确定的时间和频率范围内的能量分布，从而反映出信号在时频域上的非平稳信息。常用的时频分析方法有短时傅里叶变换，连续小波变换，和维格纳‑威尔分布。本发明采用其中的任何一种方法均可在下一步中实现时频流形的提取。
通过三种传统的时频分析方法，的时频分布函数可以用统一表示。它又可以表示为如下形式：
上式中，Ct,f表示信号在点（t，f）处的能量，L表示频率点的个数，是矩阵TFDj的第i列。当相空间中的每一维信号都转化为时频分布之后，m维的时频分布数据就构造出来了，它是一个m×L×n的3维数据矩阵。系统的内在时频结构，即时频流形，就嵌入其中，并可通过流形学习提取出来。
步骤204：将所述第二信号进行转换，生成具有所述预设嵌入维数的所述时频分析信号。
由于流形学习算法的输入是一个多维的数据集，这个数据集可以看成是一个矩阵，所以需要将每一维的时频分布矩阵转化成向量。这可以通过首尾连接时频分布矩阵的每一列来实现。即上式可以改写成如下向量：
最终，一个新的m维数据就形成了，即：
TFI=[TF1 TF2 … TFm]T=[z1 z2 … zL×n]
上式中，zi(∈Rm)表示m维空间中的一个点。矩阵TFI就是具有所述预设嵌入维数的时频分析信号，即m维的时频分析信号，同时，该矩阵也是下一步中流形学习算法的输入数据。
步骤102：对所述高维时频分析信号进行流形学习，生成具有所述预设本征维数的低维流形结构信号，所述预设本征维数小于所述预设嵌入维数。
从m维的时频分析信号中提取时频流形特征。信号的时频流形揭示了系统内在的动态信息，它只存在于高维相空间中的低维子空间，可看成是信号时频结构的骨架。而信号的干扰噪声存在于相空间中的所有维数中，且是一种随机信息，通过在m维的时频分析信号中进行流形学习，稳固的骨架结构能够被保持，而随机白噪声则会被去除。常用的流形学习算法有局部线性嵌入（LLE），等距映射（IsoMap），局部切空间排列（LTSA）。本发明采用局部切空间排列作为流形学习算法，因为它对参数具有鲁棒性，而且在主流形重构上具有优势。
局部切空间排列的算法认为流形是局部线性的。因而从数据的高维空间到它的局部切空间存在一个线性映射，从数据的低维流形到这个切空间也存在一个线性映射。局部切空间排列的主要思想是，首先利用各个数据点的k点近邻域近似构造其局部切空间，然后通过一个转换矩阵把邻域点在切空间中的局部坐标映射成d（小于m）维全局坐标，这里的d即为所述预设本征维数。从局部切空间坐标到低维全局表示的过程中，有一个仿射误差需要被最小化。而这个最小化过程可以通过求解对应于一个排列矩阵B的d个最小非零特征值的特征向量来实现。下面结合附图3简单介绍局部切空间排列的算法。
将前一个步骤中构造出的m维时频分析信号TFI作为输入数据，局部切空间排列的算法主要包含以下三步：
步骤301：计算所述时频分析信号的局部信息。
对于TFI中的每一点zi，选择其包含自身在内的k个最近邻点组成邻域矩阵Zi。中心化Zi，使之成为其中是Zi的均值向量，ek是包含k个1的向量。然后计算中心化矩阵的前d个最大右奇异向量，组成局部切空间矩阵Vi，即获得所述时频分析信号的局部信息。
步骤302：根据所述局部信息计算所述时频分析信号的排列矩阵。
选取0‑1选择矩阵Si，并通过Vi计算矩阵Wi。最后通过矩阵Si和Wi构造出所述时频分析信号的排列矩阵B。
步骤303：根据所述排列矩阵计算具有所述预设本征维数的全局坐标，所述全局坐标即为具有所述预设本征维数的低维流形结构信号。
计算排列矩阵B的前d+1个最小特征值对应的特征向量。这样就得到了d维全局坐标TFo(∈Rd×(L×n))。它里面的各个行向量对应于矩阵B的第2个到第（d+1）个最小特征值。
输出矩阵TFo=[U2 U3 … Ud+1]T即为提取出来的流形结构信号，它与输入矩阵TFI=[TF1 TF2 … TFm]T具有相同的数据点数，只是维数由m维减少为d维。
步骤103：对所述流形结构信号进行转换，生成包含多维数据的时频流形结构信号，所述时频流形结构信号中的第一维数据为所述动态信号的时频流形特征数据。
通过把TFo的第i维数据Ui+1反变换为矩阵，第i维数据流形就转化成了包含多维数据的时频流形结构信号，记为TFMi。本发明把时频流形结构信号的第一个时频结构，TFM1，作为所分析信号的时频流形特征。最近邻域点数k略大于嵌入维数m。
需要说明的是，本实施例中所采用的流形学习的算法为局部切空间排列算法，但是，采用等距映射算法或局部线性嵌入算法，同样能够实现本发明对动态信号进行分析的目的。对此，本实施例不再进行详细说明。
本实施例公开的动态信号分析方法，通过对时频分析以及流形学习的具体步骤的进一步描述，能够更充分考虑动态信号的非线性以及信号中噪声的影响，使得最终的时频流形特征数据能够作为评估系统的健康状况的准确依据。
实施例三
基于以上实施例二的描述，可以推出本发明的计算复杂性主要体现在流形学习的输入数据的总量（m×L×n）上。取较小的嵌入维数m、频率点数L和时间点数n可以提高计算效率。对于给定数据长度的信号，在保证流形提取效果的前提下，本发明分别从以下两个角度来减少数据量。
1，在所述第一信号整个频带范围内通过增加短时傅里叶变换窗函数的移动步长来减少时间点数。
对相空间中具有所述预设嵌入维数的第一信号的每一维时间序列分别进行短时傅里叶变换，生成具有所述预设嵌入维数的第二信号，即获取m维时频分布。
相空间中第j维信号的短时傅里叶变换结果为：
h(t)是位于t=0，f=0处的短时分析窗。短时傅里叶变换的频谱图（Spectrogram）即为其相应的时频分布函数即：
通常，短时分析窗在运算时一次只移动一个时间点，得到的时频分布在时间轴上仍然有n个点。当窗函数的移动步长为r（>1）时，得到的时频分布的时间点数为n′=n/r，这样就大大减少了时频分布矩阵的尺寸，此时，上式可以展开为如下形式：
上式中，Ct,f表示信号在点（t，f）处的能量，L表示频率点的个数，是矩阵TFDj的第i列。当相空间中的每一维信号都转化为时频分布之后，m维的时频分布数据就构造出来了，它是一个m×L×n′的3维数据矩阵。机械故障的内在时频结构，即时频流形，就嵌入其中，并可通过流形学习提取出来。
此方法仅适用于利用短时傅里叶变换进行时频分析的时频流形特征提取。因为短时傅里叶变换是分析局部时间区域内的频谱特性；而连续小波变换和维格纳‑威尔分布强调每一时刻点的频率分布，增加窗函数的移动步长会丢失一些时刻点的频率信息。采用该方法，即对所述第一信号整个频带范围内的部分时间点数进行时频分析，生成具有所述预设嵌入维数的第二信号，可以使本发明适用于分析健康或有局部故障的机械动态信号。
需要说明的是，其他步骤与实施例二中的相关步骤相类似，具体可参见实施例二中的相关描述。
为了进一步体现计算复杂度处理后，本方案的优势，下面将通过以下示例进行详细说明。
示例一
齿轮箱健康状况特征分析
实验是对一个汽车变速齿轮箱进行疲劳试验。该齿轮箱具有5个前进档和一个倒档。测试齿轮是变速器的第三档驱动齿轮，试验时其啮合频率和旋转频率分别为500Hz和20Hz。在整个试验过程中，齿轮箱的振动信号是通过固定在外壳上的加速度传感器和数据采集箱采集并存储的，采样频率为3kHz。测试齿轮在试验之初的一段时间内被认为是无故障阶段，在发生断齿之前的一段时间内被认为是严重磨损故障阶段，这两个阶段的振动数据被用于本发明的实验验证。分析中短时傅里叶变换的窗函数的移动步长为5。
首先对无故障阶段的健康信号进行分析。图4(a)显示了所分析信号的波形图，从图中可以看出，信号淹没在噪声之中，难以从时域波形中判断齿轮箱的健康状况。图4(b)显示了该信号的功率谱图，图中500Hz啮合频率处有一个明显的尖峰，但并不知道该频率发生的时刻，因而也难以判断齿轮箱是否有故障。这两个图都不能表示信号完整的非平稳信息。图4(c)显示了信号通过短时傅里叶变换获得的时频分布图，它能够表示信号在时频域上完整的非平稳信息，但依然受噪声的干扰。采用本发明的方法，信号的时频流形如图4(d)所示，它比图4(c)中的时频分布图更加清晰，去除了噪声的干扰，反映了信号内在的时频结构。从图4(d)中可以看出，信号只包含测试齿轮的啮合频率成分，且该频率发生在该阶段的每一个时刻，因而可以判断该阶段的齿轮箱是健康的。
对于有严重磨损故障的齿轮振动信号，其时域波形如图5(a)所示，它也包含了大量的噪声，导致不能判断故障的存在。图5(b)显示了信号的功率谱图，它包含了测试齿轮的啮合频率，同时在210Hz到350Hz的范围内有明显的幅值，暗示了信号的奇异性，但由于缺乏时域信息，所以难以判断故障的物理特性。图5(c)显示了短时傅里叶变换的频谱图，它结合了时域和频域上的信息，反映了信号频率随时间的变化，但背景噪声的存在干扰了对故障特征的判断。最后用本发明中的方法来分析该信号，其结果如图5(d)所示。从图中可以看出，时频流形的结果比图5(c)中的时频分布要清晰很多。图5(d)所示的时频流形里只有500Hz和210‑350Hz的频带内有特征。500Hz的频率发生在每一个时刻，对应于健康信号的啮合频率；210‑350Hz的频带内周期性出现能量脉冲，反映出测试齿轮局部故障的周期性脉冲特性。故障脉冲的重复频率为1/0.05s=20Hz,恰好等于测试齿轮的旋转频率。所以，时频流形特征能够用来判断齿轮箱局部故障的存在。
从本示例中可以看出，采用本发明提取的时频流形特征带有信号本质流形的特点，它去除了噪声的干扰，更便于观察齿轮箱的健康状况。
示例二
列车轴承局部故障特征分析
典型的轴承局部故障（如局部的剥落、腐蚀、划痕等）常常发生在轴承的外圈、内圈和滚子上，可分别在这三个部位设置局部故障来测试轴承的振动信号。试验通过对国内常用的列车轴承进行振动测试，采集其声压信号进行故障特征分析。列车轴承为单排圆柱滚子轴承，型号为NJ(P)3226X1。局部故障是用线切割分别在轴承的外圈、内圈和滚子上设置凹痕。声压信号是通过麦克风和数据采集箱采集并存储的，采样频率为12.5kHz。试验时轴的旋转速度为1430rpm，测试轴承的承载为30kN。列车轴承的故障参数和运动学特征如表1所示。
表1列车轴承的故障参数及运动学特征
采用本发明中的方法对这三类有局部故障的声音信号进行分析，取短时傅里叶变换的窗函数的移动步长为5。图6显示了对有外圈故障的声音信号的分析结果。图6(a)为声音信号的波形图，图6(b)为该信号的功率谱图，图6(c)为短时傅里叶变换的频谱图，图6(d)为信号的时频流形图。从图6中可以看出，信号的时域波形因受噪声的影响而不能判断周期性故障脉冲的存在；功率谱图因缺乏频率随时间的变化规律，不能观察故障的物理机理；短时傅里叶变换的频谱图虽然能展示信号完整的非平稳信息，但在很宽的频率范围内存在大量的噪声干扰；时频流形只保留了反映故障流形的时频结构，去除了大量的噪声，提高了信噪比和时频分辨率。时频流形中出现了周期性的能量聚集区，反映了故障脉冲的存在。故障脉冲的时间间隔约为7.21ms，说明列车轴承的局部故障发生在外圈。
图7显示了对有内圈故障的声音信号的分析结果。图7(a)为声音信号的波形图，图7(b)为该信号的功率谱图，图7(c)为短时傅里叶变换的频谱图，图7(d)为信号的时频流形图。从图7中可以看出，时域波形里的故障脉冲淹没在噪声之中，不能作为故障存在的依据；功率谱图因缺乏频率随时间的变化规律，不能观察故障的物理机理；短时傅里叶变换的频谱图虽然能展示信号完整的非平稳信息，但在很宽的频率范围内存在大量的噪声干扰；时频流形能够去除大量的噪声，提高故障信息的信噪比和时频分辨率。时频流形里故障脉冲的时间间隔约为5.13ms，说明列车轴承的局部故障发生在内圈。
图8显示了对有滚子故障的声音信号的分析结果。图8(a)为声音信号的波形图，图8(b)为该信号的功率谱图，图8(c)为短时傅里叶变换的频谱图，图8(d)为信号的时频流形图。从图8中可以看出，信号的时域波形和功率谱图都不能判断故障的存在；短时傅里叶变换的频谱图里也存在大量的非相关噪声的干扰；时频流形去除了大部分的噪声，保留了故障脉冲的能量聚集区，提高了信噪比和时频分辨率。时频流形里故障脉冲的时间间隔约为7.27ms，说明列车轴承的局部故障出现在滚子上。
从本示例中可以看出，本发明所采用的方法能够从轴承声音信号中提取出反映局部故障流形特点的时频结构，从而更方便的判断局部故障发生的位置。
2，在不减少时频分布时间点数的情况下，仅在所述第一信号整个频带范围内选取故障频带进行时频流形特征的提取。
当机械设备发生局部故障时，其动态信号表现为周期性脉冲形式，这些脉冲的能量主要集中在时频分布的某个频率带范围内，选取出这个频率带进行流形学习，可以排除其它频带内非相关分量的干扰，并直接提取出故障脉冲的时频流形特征。该方法适用于三种常用的时频分析方法，但是只能分析有局部故障的机械动态信号。下面介绍一种故障频带自动选取算法。
对于给定的时频分布TFD(t,f)∈RL×N，设其频带序列为F=[f1,f2,…,fL]，频率步长为fstep，它的行向量的和函数为：
上式是关于频率f的函数。它的全局最大值出现在故障频带的中央频率fm处。设阈值为1/e，则故障频带的下限fl和上限fh分别取SMR（fm）向左和向右减少1/e倍的频率。也就是说，上下限频率处的能量与故障中心带的能量关系为：SMR2(fl)=SMR2(fh)=SMR2(fm)/e2，其中，f1≤fl<fm<fh≤fL，这样就把时频分布限定在仅包含故障信息的窄频带序列F′=[fl,fl+1,…,fm,…,fh]内，它的频率点数为L′=(fh‑fl)/fstep+1<<L，这样本发明的计算效率就大大提高了。
有关于如何采用上述故障频带选取方法实现本发明的目的，将通过以下描述进行详细说明。
对相空间中的每一维时间序列分别进行时频分析，获取m维时频分布。常用的时频分析方法有短时傅里叶变换，连续小波变换，和维格纳‑威尔分布。本发明采用其中的任何一种方法均可在下一步中实现时频流形的提取。为了消除维格纳‑威尔分布的交叉项，本发明采用平滑伪维格纳‑威尔分布。
对于相空间中第j维信号其短时傅里叶变换结果为：
短时傅里叶变换的频谱图（Spectrogram）即为其相应的时频分布函数即：
的连续小波变换结果为：
上式中，a为尺度因子，b为平移因子，为小波母函数。CWT(a,b)表示在点(a,b)的小波系数。平移因子b对应于时间点t，而尺度因子a与频率f的对应关系为：
其中fs为采样频率，fc为母小波中心频率。这样，CWTj(a,b)可以转化为以t和f为变量的小波系数函数CWTj(t,f)。小波变换的尺度谱（Scalogram）（以t和f为坐标）即为其相应的时频分布函数即：
的平滑伪维格纳‑威尔分布结果为：
上式中，g(t)和h(t)分别为时间平滑窗和时域上的频率平滑窗。平滑伪维格纳‑威尔分布的绝对值即为其相应的时频分布函数即：
的以上三种时频分布函数可以用统一表示，它是一个L×n的矩阵，包含了所有频带范围内的非平稳信息。为了提取出故障频带内的非平稳信息，采用以上故障频带自动选取算法。
对于接收的动态信号的时频分布TFD(t,f),设其频带序列为F=[f1,f2,…,fL]，频率步长为fstep，它的行向量的和函数为：
上式是关于频率f的函数。它的全局最大值出现在故障频带的中央频率fm处。设阈值为1/e，则故障频带的下限fl和上限fh分别取SMR(fm)向左和向右减少1/e倍的频率。也就是说，上下限频率处的能量与故障中心带的能量关系为：SMR2(fl)=SMR2(fh)=SMR2(fm)/e2，其中f1≤fl<fm<fh≤fL，这样就把时频分布限定在仅包含故障信息的窄频带序列F′=[fl,fl+1,…,fm,…,fh]内，它的频率点数为L′=(fh‑fl)/fstep+1<<L，这样本发明的计算效率就大大提高了。
故障频带选出之后，对于在整个频带范围内的时频分布选择故障频带F′对应的时频分布，用TFDj(t,f)表示，它又可以展开为如下形式：
上式中，是矩阵TFDj的第i列，其数据长度为L′。当相空间中的每一维信号都转化为时频分布之后，m维的时频分布数据就构造出来了，它是一个m×L′×n的3维数据矩阵。采用首尾连接各维时频分布矩阵每一列的方式来满足流形学习算法对输入数据的要求，即上式可以改写为如下列向量：
最终，一个新的m维数据就形成了，即：
TFI=[TF1 TF2 … TFm]T=[z1 z2 … zL′×n]
上式中，zi(∈Rm)表示m维空间中的一个点。矩阵TFI就是下一步中流形学习算法的输入数据。
需要说明的是，其他步骤与实施例二中的相关步骤相类似，具体可参见实施例二中的相关描述。
为了进一步体现计算复杂度处理后，本方案的优势，下面将通过以下示例进行详细说明。
示例三
球轴承局部故障特征分析
试验数据来自美国凯斯西储大学轴承数据中心。测试轴承为深沟球轴承，其型号为6205‑2RS JEM SKF。试验通过电切割的方式分别在外圈、内圈和滚珠上设置单点故障，故障尺寸、测试时轴承的旋转速度、故障频率和周期如表2所示。采样频率为12kHz。
表2深沟球轴承的故障参数及运动学特征
图9显示了轴承外圈故障的特征分析结果，其时频分析方法为短时傅里叶变换，频带选取时阈值取1/3。图9(a)为该振动信号的波形图，图9(b)为该信号的功率谱图，图9(c)为短时傅里叶变换的频谱图，图9(d)为提取出的时频流形图。从图9可见，受噪声的影响，时域波形和频域功率谱都不能观察出信号有周期性故障特征的存在，而时频流形能在时频域中反映出故障脉冲完整的非平稳信息，且与短时傅里叶变换的频谱图相比，去除了噪声的干扰，具有更高的信噪比和时频分辨率。时频流形里故障脉冲的时间间隔约为9.71ms，说明轴承的局部故障发生在外圈。
图10显示了轴承滚珠故障的特征分析结果,其时频分析方法为连续小波变换，频带选取时阈值取1/2。图10(a)为该振动信号的波形图，图10(b)为该信号的功率谱图，图10(c)为连续小波变换的频谱图，图10(d)为提取出的时频流形图。从图10中可以看出，信号的时域波形和功率谱图都不能诊断出局部故障的存在。连续小波变换的频谱图中显示了在3.3kHz左右的频带范围内有故障脉冲的存在，但是频带内的噪声降低了故障脉冲的时频分辨率。信号的时频流形特征只保留了故障信息的内在时频结构，强调了故障脉冲特性，因而具有更高的信噪比和时频分辨率。故障脉冲的时间间隔约为7.28ms，说明轴承的局部故障发生在滚珠上。
图11显示了轴承内圈故障的特征分析结果。图11(a)为该振动信号的波形图，图11(b)为该信号的功率谱图，图11(c)为维格纳‑威尔分布图，图11(d)为提取出的时频流形图。从图11中可以看出，与维格纳‑威尔分布图相比，信号的时频流形特征能排除故障脉冲之间的噪声，提高信号时频结构的信噪比和时频分辨率。故障脉冲的时间间隔约为6.34ms，说明轴承的局部故障发生在内圈。
本示例说明，不管采用哪一种时频分析方法，本发明都能提取出故障脉冲的时频流形结构，从更好地判断出轴承局部故障发生的位置。
本实施例通过考虑机械系统非平稳表示的非线性特性，为系统状态模式提供了一种体现流形特点的新的时频表示方法，构造出一种新的特征标签—时频流形。新的特征不仅携带有非平稳信息，而且揭示了系统状态的内在流形结构，因而反映了机械系统状态模式的本质。这正是现有技术所缺乏的。与现有时频表示方法相比，本发明更具噪声干扰小、时频分辨率高的优点。
上述本发明公开的实施例中详细描述了方法，对于本发明的方法可采用多种形式的装置实现，因此本发明还公开了一种装置，下面给出具体的实施例进行详细说明。
实施例四
请参阅附图12，为本发明实施例四公开的一种动态信号分析装置具体结构示意图。该动态信号分析装置具体包括：
时频分析单元1：用于对接收到的动态信号进行时频分析，生成与所述动态信号相对应的具有预设嵌入维数的高维时频分析信号；
流形学习单元2：用于对所述时频分析信号进行流形学习，生成具有所述预设本征维数的低维流形结构信号，所述预设本征维数小于所述预设嵌入维数；
流形转换单元3：用于对所述流形结构信号进行转换，生成包含多维数据的时频流形结构信号，所述时频流形结构信号中的第一维数据为所述动态信号的时频流形特征数据。
有关于该装置中上述各个单元的具体功能及其实现过程请参考实施例一中的相关描述，本实施例不再详细说明。
通过采用本实施例中公开的动态信号分析装置，对接收的动态信号进行时频分析、流形学习以及转换，能够充分考虑动态信号的非线性以及信号中噪声的影响，使得最终的时频流形特征数据能够作为评估系统健康状况的准确依据。
基于上述实施例四，本发明还对上述装置进行了进一步详细说明，下面将通过以下实施例进行详细描述。
实施例五
请参阅附图13，为本实施例五公开的动态信号处理装置中时频分析单元的具体结构，该时频分析单元具体包括：
接收子单元11：用于接收所述动态信号；
相空间重构子单元12：用于根据预设嵌入维数，将所述动态信号在相空间中重构，生成具有所述预设嵌入维数的第一信号，所述预设嵌入维数至少为12维；
时频分析子单元13：用于根据预设时频分析方法，对所述第一信号整个频带范围内的所有时间点数进行时频分析，生成具有所述预设嵌入维数的第二信号，所述预设时频分析方法包括短时傅里叶变换、连续小波变换或维格纳—威尔分布变换；
转换子单元14：用于将所述第二信号进行转换，生成具有所述预设嵌入维数的所述时频分析信号。
请参阅附图14，为本实施例五公开的动态信号分析装置中，流形学习单元的具体结构示意图，该流形学习单元包括：
局部信息计算子单元21：用于计算所述时频分析信号的局部信息；
排列矩阵计算子单元22：用于根据所述局部信息计算所述时频分析信号的排列矩阵；
全局坐标计算子单元23：用于根据所述排列矩阵计算具有所述预设本征维数的全局坐标，所述全局坐标即为具有所述预设本征维数的低维流形结构信号。
请参阅附图15，为本实施例五公开的动态信号分析单元中时频分析子单元的具体结构示意图，当所述预设时频分析方法为短时傅里叶变换时，该时频分析子单元具体包括：
第一计算模块131：在所述第一信号整个频带范围内通过增加短时傅里叶变换窗函数的移动步长减少时间点数；
第二计算模块132：对所述第一信号整个频带范围内的部分时间点数进行时频分析，生成具有所述预设嵌入维数的第二信号。
需要说明的是，通过增加短时傅里叶变换窗函数的移动步长减少时间点数与对时间点数进行时频分析是同时进行的，在本实施例中，只是将时频分析子单元的上述两个功能分别通过不同的模块进行展现，但是两个模块执行其对应功能是没有先后顺序之分的，而是同时进行。
请参阅附图16，为本实施例五公开的动态信号分析单元中时频分析子单元的另一种具体结构示意图，当所述预设时频分析方法为短时傅里叶变换、连续小波变换方法或维格纳—威尔分布变换方法时，该时频分析子单元具体包括：
故障频带选择模块141：用于在所述接收到的动态信号整个频带范围内选取故障频带；
第三计算模块142：用于对所述第一信号故障频带范围内的所有时间点数进行时频分析，生成具有所述预设嵌入维数的第二信号。
需要说明的是，在对第一信号进行时频分析的时候，是只针对故障频带而言的，而故障频带的选取，是在接收到动态信号之后，对接收到的动态信号进行相空间重构之前就确定的，也就是说，在确定好动态信号的故障频带以后，在对动态信号进行相空间重构的过程与前面任意一个实施例中的相关描述是相同的，只是在时频分析的时候，只选取故障频带对应的时频分布。
请参阅附图17，为本实施例五公开的动态信号分析装置中，局部信息计算子单元的具体结构示意图，该局部信息计算子单元包括：
待处理元素选取模块211：用于在所述时频分析信号中选取待处理元素；
元素选取模块212：用于选取与所述待处理元素最近的预设个数的元素，所述预设个数大于所述预设嵌入维数；
中心化处理模块213：用于将所述预设个数的元素进行中心化处理得到中心化矩阵；
局部信息计算模块214：用于计算所述中心化矩阵的右奇异向量，并按照右奇异值从大到小的顺序，选取所述预设本征维数个数右奇异值相对应的右奇异向量，作为所述时频分析信号的局部信息。
请参阅附图18，为本实施例五公开的动态信号分析装置中，排列矩阵计算子单元的具体结构示意图，该排列矩阵计算子单元包括：
第一矩阵计算模块221：用于根据所述时频分析信号以及所述中心化矩阵计算得到第一矩阵；
第二矩阵计算模块222：用于根据所述时频分析信号的局部信息计算得到第二矩阵；
第三矩阵计算模块223：用于根据第一矩阵以及第二矩阵计算所述时频分析信号的排列矩阵。
需要说明的是，本实施例中公开的装置对应于前面实施例中所描述的方法，有关于整个装置，以及装置中各个子单元以及模块的功能及其实现过程具体请参见前述实施例中的相关描述，本实施例中不再进行详细说明。
采用本实施例公开的动态信号分析装置对动态信号进行分析的过程中，通过考虑机械系统非平稳表示的非线性特性，为系统状态模式提供了一种体现流形特点的新的时频表示方法，构造出一种新的特征标签—时频流形。新的特征不仅携带有非平稳信息，而且揭示了系统状态的内在流形结构，因而反映了机械系统状态模式的本质。这正是现有技术所缺乏的。与现有时频表示方法相比，本发明更具噪声干扰小、时频分辨率高的优点。
综上所述：
本发明公开了一种动态信号分析方法及装置。首先，对接收到的动态信号进行时频分析，生成与动态信号相对应的具有预设嵌入维数的高维时频分析信号；然后，对时频分析信号进行流形学习，生成具有预设本征维数的低维流形结构信号，预设本征维数小于所述预设嵌入维数；最后，对流形结构信号进行转换，生成包含多维数据的时频流形结构信号，时频流形结构信号中的第一维数据为动态信号的时频流形特征数据。通过对接收的动态信号进行时频分析、流形学习以及转换，能够充分考虑动态信号的非线性以及信号中噪声的影响，使得最终的时频流形特征数据能够作为评估系统的健康状况的准确依据。
本说明书中各个实施例采用递进的方式描述，每个实施例重点说明的都是与其他实施例的不同之处，各个实施例之间相同相似部分互相参见即可。对于实施例公开的装置而言，由于其与实施例公开的方法相对应，所以描述的比较简单，相关之处参见方法部分说明即可。
结合本文中所公开的实施例描述的方法或算法的步骤可以直接用硬件、处理器执行的软件模块，或者二者的结合来实施。软件模块可以置于随机存储器（RAM）、内存、只读存储器（ROM）、电可编程ROM、电可擦除可编程ROM、寄存器、硬盘、可移动磁盘、CD‑ROM、或技术领域内所公知的任意其它形式的存储介质中。
对所公开的实施例的上述说明，使本领域专业技术人员能够实现或使用本发明。对这些实施例的多种修改对本领域的专业技术人员来说将是显而易见的，本文中所定义的一般原理可以在不脱离本发明的精神或范围的情况下，在其它实施例中实现。因此，本发明将不会被限制于本文所示的这些实施例，而是要符合与本文所公开的原理和新颖特点相一致的最宽的范围。