本发明涉及机械的故障诊断领域，尤其是涉及一种基于振动信号联合谱相关系数分析的轴流风机状态识别方法。
背景技术
轴流风机在长时间使用过程中会出现一些故障，目前基于轴流风机的振动信号的故障诊断方法主要有两类：时域方法和频域方法。
时域方法主要分为两类，第一类时域方法是时域波形分析，直接对振动时域信号进行整体观察和分析，主要关心振动信号峰值的大小、振动的整体幅度大小、振动信号峰值对应的时间、同一形状的波形重复出现的周期长短等。另一类时域方法是时域统计参数分析，通过截取一段振动时域信号，然后对其进行各类统计参数的计算，如平均值、峰值、峰峰值、标准差、方差、均方根等，通过与已有的标准对比，进而判断当前风机的振动是否处于正常范围内。
频域方法主要有快速傅里叶变换、短时傅里叶变换以及小波变换等。快速傅里叶变换是将信号的时域信息完全转换成频域信息的最常用方法，是一种快速高效的算法，能够较为完美地描述平稳信号的特征。而短时傅里叶变换是将时域信号加窗，将加窗之后的信号进行快速傅里叶变换，从而得到信号的时频联合分布。小波变换也是对信号进行加窗分析，但其窗函数形状可变，进而保证信号分析结果在低频处有较高的频率分辨率，在高频处有较高的时间分辨率，进而实现信号的自适应分析。
但以上方法中，时域波形分析的方式太过简单粗糙，难以有针对性地获取有用的风机状态信息；时域统计参数分析所获取的参数只能判断振动幅度的大小，不能进一步判断轴流风机的运行状态；频域方法能将振动信号中的时域信息转换成频域信息，能有效获取振动的主要频率成分，但是也需要通过频谱分析进一步推测轴流风机可能的故障原因，且需要检查风机零部件以验证推测，不能够直观有效地获取轴流风机的运行状态信息。
技术实现要素：
本发明提供了一种基于振动信号联合谱相关系数分析的轴流风机状态识别方法，能够对实时运行中的轴流风机的运行状态进行准确有效的识别判断，可广泛运用于轴流风机的实时监测、状态识别和故障诊断等领域。
一种基于振动信号联合谱相关系数分析的轴流风机状态识别方法，包括以下步骤：
(1)营造轴流风机的各种工况状态，分别采集振动时域信号序列，建立标准工况数据库，经过数据处理后得到标准工况幅值谱数据库和倒频谱数据库；
(2)采集实时运行的轴流风机振动信号，经过数据处理后得到轴流风机实时工况的振动信号幅值谱和倒频谱；
(3)将实时工况的振动信号幅值谱依次与标准工况幅值谱数据库中的每一组幅值谱进行相关性分析，得到一系列幅值谱相关系数；
(4)将实时工况的振动信号倒频谱依次与标准工况倒频谱数据库中的每一组倒频谱进行相关性分析，得到一系列倒频谱相关系数；
(5)将幅值谱相关系数与倒频谱相关系数进行坐标组合，构成一系列联合谱相关系数点；
(6)在直角坐标系中计算联合谱相关系数点与点(1，1)的绝对距离，根据绝对距离最小的联合谱相关系数点确定轴流风机当前的运行状态。
本发明的方法能够克服已有的时域方法、频域方法不能够直接有效地对轴流风机的运行状态进行识别诊断的缺陷，该方法简便易行，能够建立起轴流风机标准工况库，并对实时运行中的轴流风机的运行状态进行准确有效的识别判断。
步骤(1)中，所述的工况状态包括正常工况、转子不平衡工况和地脚螺栓松动工况；所述的转子不平衡工况通过在叶片根部增加配重设定若干组工况，所述的地脚螺栓松动工况根据不同地脚螺栓松动的组合形式设定若干组工况。
作为优选，步骤(1)中，对于同一种工况，以固定的采样频率重复采集三组振动时域信号序列。用这三组振动原始数据作为标准工况建立轴流风机标准工况数据库，并对其进行快速傅里叶变换得到相应的幅值谱，建立起轴流风机幅值谱数据库；对其进行倒频谱变换得到相应的倒频谱，建立起轴流风机倒频谱数据库。
每一组工况下重复采集三组数据，可以创建更加丰富有效的轴流风机标准工况数据库，并且可以有效避免实验偶然性带来的误差和不确定性，降低所建立的轴流风机标准工况数据库存在错误工况的可能性。
作为优选，步骤(1)和步骤(2)之间增加以下步骤：再次营造轴流风机的各种工况状态，令轴流风机在n种不同的工况下运行，在这n种不同工况下，以某一固定的采样频率重复测三组振动原始数据，以这三组振动原始数据作为验证工况，运用基于振动信号联合谱相关系数分析的轴流风机状态识别方法进行诊断验证，统计状态识别和故障诊断的成功率。
将验证正确的验证工况加入到轴流风机标准工况数据库中，并将其对应的幅值谱、倒频谱加入到轴流风机幅值谱数据库、轴流风机倒频谱数据库中。
在重复工况下采集三组数据，并用已有方法进行诊断验证，可以在已知具体工况的条件下，对已有方法的正确率进行全方位的验证。将验证正确的验证工况加入轴流风机标准工况数据库中，能够进一步丰富轴流风机标准工况数据库，实现标准工况的原始数据的自适应扩展。
步骤(1)和步骤(2)中，所述的数据处理具体为：对振动时域信号序列分别进行快速傅里叶变换和倒频谱变换，得到振动信号的幅值谱和倒频谱。其中，所述的倒频谱变换具体步骤如下：
(1-1)对振动时域信号序列进行快速傅里叶变换；
(1-2)对所得结果先取绝对值，然后取自然对数；
(1-3)在取自然对数基础上，进行快速傅里叶逆变换，得到振动时域信号序列的每个点所对应的实倒频谱值；
(1-4)以振动时域信号序列每个点的时间值为横坐标，以每个点所对应的实倒频谱值为纵坐标，绘制振动时域信号的倒频谱。
通过以上处理和分析，可以有效地提取出轴流风机振动信号中的主要调制频率成分，在结合轴流风机的机械运转规律的基础上，能够进一步有效识别轴流风机的轴频和叶频。
步骤(3)和步骤(4)中，关于相关系数的计算，可以分为三步，具体为：
(3-1)计算两个数组的标准差，计算公式如下：
式中，x代表第一个数组，y代表第二个数组；xi代表数组x中第i个元素，yi代表数组y中第i个元素；μx代表数组x的期望值，μy代表数组y的期望值；σx代表数组x的标准差，σy代表数组y的标准差；n代表数组x和数组y中的元素总数。
(3-2)计算两个数组之间的协方差，计算公式如下：
cov(x，y)＝e[(x-μx)(y-μy)]
式中，x代表第一个数组，y代表第二个数组；cov(x，y)代表数组x和数组y之间的协方差；符号e代表求期望值。
(3-2)计算两个数组之间的相关系数，计算公式如下：
式中，ρ代表数组x与数组y的相关系数。
步骤(5)的具体步骤为：在直角坐标系中，以幅值谱相关系数为横坐标值，以倒频谱相关系数为纵坐标值，进行联合谱相关系数点的绘制，依次产生一系列联合谱相关系数点。
步骤(6)的具体步骤为：
(6-1)依次比较每个联合谱相关系数点到点(1，1)之间的绝对距离，得到离点(1，1)绝对距离最小的一个联合谱相关系数点；
(6-2)根据找到联合谱相关系数点，确定工况编号，从而找到对应的标准工况；
(6-3)将这组标准工况认定为轴流风机当前的实时运行状态。
本发明提供的技术方案的有益效果如下：
(1)本发明提出了一种基于振动信号联合谱相关系数分析的轴流风机状态识别方法，该方法能够对轴流风机的运行状态进行有效的识别和判断。
(2)本发明能够建立起轴流风机典型工况的标准工况数据库，并在此基础上，建立对应的标准工况幅值谱数据库和标准工况倒频谱数据库。
(3)本发明运用曲线的相关系数分析方法，将实时工况的幅值谱(倒频谱)与标准工况幅值谱(倒频谱)数据库中的数据依次进行相关系数分析，从幅值谱(倒频谱)的角度出发，得到最接近的一组标准工况。
(4)本发明联合幅值谱和倒频谱的相关系数分析，幅值谱提供振动信号的频率成分信息，倒频谱提供振动信号的主要调制频率信息，从幅值谱和倒频谱两个角度观察，得到轴流风机最有可能的实时运行状态。
附图说明
图1为本发明一种基于振动信号联合谱相关系数分析的轴流风机状态识别方法的流程示意图；
图2为本发明实施例的实验系统结构示意图；
图3为本发明实施例的风机底座部位俯视图；
图4为本发明实施例轴流风机实时工况的振动时域信号序列；
图5为本发明实施例轴流风机实时工况的振动信号幅值谱；
图6为本发明实施例轴流风机实时工况的振动信号倒频谱；
图7为本发明实施例轴流风机实时工况的联合谱相关系数分布图。
具体实施方式
为了使本发明的目的、技术方案和优点更加清楚，下面将举实施例并结合附图对本发明作进一步地详细描述，但本发明并不限于下述的实施例。
如图1所示，一种基于振动信号联合谱相关系数分析的轴流风机状态识别方法，包括以下步骤：
1、营造轴流风机的各种工况状态，建立轴流风机标准工况数据库。
1-1、本发明在轴流风机上进行实验，令轴流风机在n种不同的工况下运行，通过人为营造各种工况，测得其振动数据进行分析。实验系统如图2所示，交变控制箱11主要用于控制轴流风机的运行与转速；出口管道13与风机12的出口紧密连接，管道上附有压力、温度计、阀门等；采集仪14负责将振动加速度传感器采集到的振动信号转化为数据存储起来，与它相连的pc利用远程桌面访问采集仪的主机来调取采集到的数据进行分析。加速度传感器15分别安装在风机12的垂直、水平方向和轴向位置，用来测定风机不同方向的加速度进行对照分析；压力和温度传感器16，用来测定风机故障运行时出口处的压力与温度，观察风机在故障运行时出口压力温度是否有较大的变化；出口开关17，通过旋转开关可以控制风机出口的开度。
分别令轴流风机在正常工况、转子不平衡工况、地脚螺栓松动工况下运行，其中正常工况是指在轴流风机机械状态良好的情况下，以2980r/min的额定转速运行。
转子不平衡工况是通过在叶片根部的固定螺栓上增加指定配重来实现的，增加的配重分别有10g、20g等2种情况，也是以2980r/min的额定转速运行，即转子不平衡工况可分为2种。
地脚螺栓松动工况是将特定部位的螺栓松动半圈进行实验，如图3所示，风机底座共有四枚地脚螺栓，分别为螺栓1、螺栓2、螺栓3和螺栓4，实验遍历螺栓松动的所有情况，一枚地脚螺栓松动时有螺栓1松动、螺栓2松动、螺栓3松动和螺栓4松动共4种情况。两枚地脚螺栓松动时有螺栓1、螺栓2松动，螺栓1、螺栓3松动，螺栓1、螺栓4松动，螺栓2、螺栓3松动，螺栓2、螺栓4松动和螺栓3、螺栓4松动共6种情况。三枚地脚螺栓松动时有螺栓1、螺栓2、螺栓3松动，螺栓1、螺栓2、螺栓4松动，螺栓1、螺栓3、螺栓4松动，螺栓2、螺栓3、螺栓4松动共4种情况。以及四枚螺栓全部松动的情况，也是以2980r/min的额定转速运行，即地脚螺栓松动共有15种情况。
将以上18种工况依次进行编号，正常工况为1号工况，转子不平衡(不平衡质量10g)为2号工况，转子不平衡(不平衡质量20g)为3号工况，螺栓1松动为4号工况，螺栓2松动为5号工况，螺栓3松动为6号工况，螺栓4松动为7号工况，螺栓1、螺栓2松动为8号工况，螺栓1、螺栓3松动为9号工况，螺栓1、螺栓4松动为10号工况，螺栓2、螺栓3松动为11号工况，螺栓2、螺栓4松动为12号工况，螺栓3、螺栓4松动为13号工况，螺栓1、螺栓2、螺栓3松动为14号工况，螺栓1、螺栓2、螺栓4松动为15号工况，螺栓1、螺栓3、螺栓4松动为16号工况，螺栓2、螺栓3、螺栓4松动为17号工况，螺栓1、螺栓2、螺栓3、螺栓4松动为18号工况。
1-2、根据实际需要选定采样频率fs，进而确定采样时间间隔ts，选取恰当的采样总时长t。在此条件下，对轴流风机的振动信号进行测取，共测取1组正常工况，2组转子不平衡工况和15组地脚螺栓松动工况，总计18组标准工况。振动信号采样频率与采样时间间隔的关系如下：
式中，fs为振动信号的采样频率；ts为振动信号的采样时间间隔；
振动信号采样总时长与信号频率分辨率的关系如下：
式中，f0为振动信号的频率分辨率；n为振动信号的采样总点数；t为振动信号的采样总时长。
1-3、将轴流风机在以上18组标准工况下所测取的振动时域信号按照一定的命名规则进行命名，建立轴流风机标准工况数据库。
正常工况命名为normal-x，x代表第x次测试，如normal-2代表正常工况第2次实验测取的振动时域信号序列。
转子不平衡工况命名为imbalance-x-y，x代表不平衡的质量，y代表第y次测试，如imbalance-10g-3代表转子不平衡工况，其不平衡质量为10g，第三次实验测取的振动时域信号序列。
地脚螺栓松动工况命名为boltlooseness-x-y-z，x代表松动的地脚螺栓的数量，y代表松动的地脚螺栓的位置，z代表第z次测试，如boltlooseness-3-123-4代表地脚螺栓松动工况，其松动的地脚螺栓的数量是3个，松动位置是1号、2号和3号，第4次实验测取的振动信号序列。
2、建立轴流风机标准工况幅值谱数据库和倒频谱数据库。
2-1、在得到以上18种标准工况的振动时域序列的基础上，对其进行快速傅里叶变换，得到其对应的振动信号幅值谱，建立轴流风机标准工况幅值谱数据库。
matlab(用于算法开发、数据可视化、数据分析以及数值计算的高级技术计算语言和交互式环境)中用以实现快速傅里叶变换的语句为：
x(f)＝fft(x(t))
式中，x(f)为振动信号时间序列经过快速傅里叶变换得到的幅值谱；x(t)为离散的振动时域序列；fft表示快速傅里叶变换。
2-2、在得到以上18种标准工况的振动时域序列的基础上，对其进行倒频谱变换，得到其对应的振动信号倒频谱，建立轴流风机标准工况倒频谱数据库。
matlab中用以实现倒频谱变换的语句为：
c(t)＝rceps(x(t))
式中：c(t)为振动信号时间序列经过倒频谱变换得到的倒频谱；x(t)为离散的振动时域序列；rceps表示倒频谱变换。
3、采集实时运行的轴流风机振动信号，经过数据处理后得到轴流风机实时工况的振动信号幅值谱和倒频谱。
3-1、令轴流风机在正常工况下以2980r/min的转速运行，即营造轴流风机的正常工况，以采样频率fs对实时运行中的轴流风机振动信号进行采集，获得实时工况的振动时域信号序列，如图4所示。
3-2、对实时工况的振动时域信号序列进行快速傅里叶变换，得到轴流风机实时工况的振动信号幅值谱，如图5所示。
3-3、对实时工况的振动时域信号序列进行倒频谱变换，得到轴流风机实时工况的振动信号倒频谱，如图6所示。
4、分别计算实时工况与标准工况的振动信号幅值谱相关系数ρxi和倒频谱相关系数ρyi。
4-1、将轴流风机实时工况的振动幅值谱依次与轴流风机标准工况幅值谱数据库中的第i(i＝1，2，3……n，n＝18)组标准工况的幅值谱进行相关性分析，得到相应的幅值谱相关系数ρxi；
对两个大小相同的一维数组(曲线)进行相似程度的分析，即相关性系数计算，可以分为三步：
第一步，计算两个数组的标准差，计算公式如下：
式中，x代表第1个数组，y代表第二个数组；xi代表数组x中第i个元素，yi代表数组y中第i个元素；μx代表数组x的期望值，μy代表数组y的期望值；σx代表数组x的标准差，σy代表数组y的标准差；n代表数组x和数组y中的元素总数。
第二步，计算两个数组之间的协方差，计算公式如下：
cov(x，y)＝e[(x-μx)(y-μy)]
式中，x代表第1个数组，y代表第二个数组；cov(x，y)代表数组x和数组y之间的协方差；符号e代表求期望值。
第三步，计算两个数组之间的相关系数，计算公式如下：
式中，ρ代表数组x与数组y的相关系数。
4-2、将轴流风机实时工况的振动倒频谱依次与轴流风机标准工况倒频谱数据库中的第i(i＝1，2，3……n，n＝18)组标准工况的倒频谱进行相关性分析，得到相应的倒频谱相关系数ρyi；
其中，相关系数的计算过程与步骤4-1中一致。
5、绘制联合谱相关系数点并计算与点(1，1)之间的绝对距离di。
5-1、在直角坐标系中，以幅值谱相关系数ρxi为横坐标值，以倒频谱相关系数ρyi为纵坐标值，进行联合谱相关系数点ui(ρxi，ρyi)的绘制，如图7所示，依次产生n(n＝18)个联合谱相关系数点。
5-2依次计算第i(i＝1，2，3……n，n＝18)个联合谱相关系数点与点(1，1)之间的绝对距离di；
联合谱相关系数点ui(ρxi，ρyi)，其到点(1，1)之间的绝对距离di按照以下公式计算：
6、根据绝对距离di，确定轴流风机当前的实时运行状态。
6-1依次比较联合谱相关系数点ui(ρxi，ρyi)到点(1，1)之间的绝对距离di，得到离点(1，1)绝对距离最小的一个联合谱相关系数点；
matlab中用以实现寻找数组元素最小值的语句为：
[dmin，index]＝min(di)；
式中，di代表联合谱相关系数点ui(ρxi，ρyi)到点(1，1)之间的绝对距离；dmin代表di构成的数组中的元素最小值；index代表di构成的数组中的元素最小值在数组中的位置；min代表求数组最小值的操作。
6-2、根据dmin找到对应的联合谱相关系数点umin，确定工况编号，从而找到dmin对应的标准工况；
6-3、将这组标准工况认定为轴流风机当前的实时运行状态。
以上所述仅为本发明的较佳实施举例，并不用于限制本发明，凡在本发明精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。