一种基于音频的转速识别方法与流程

本申请涉及一种基于音频的转速识别方法，适用于转速识别的技术领域。

背景技术：

在旋转机械中，对于转速的识别对于生产和维护来说具有重要的意义。例如，对于风机叶片来说，无论是从风机安全运行方面还是从发电综合性能上来看，叶片转速都是风电机组最基本也是较为重要的控制量，从其他角度来看，出于对叶片运行状况的监控，叶片转速可以为后续的数据分析算法提供可靠输入，是重要的参考量。但是现实中，出于硬件条件或安全性、隐私性考虑，并不一定可以实时通过scada(数据采集与监视控制)系统采集到叶片实时转速，此时需要通过其他方法计算得到近似的实时转速，目前自动计算叶片实时转速的方法还是空白。

公开号为cn110658519a的中国专利中，通过对外辐射源雷达接收到的直达波信号和雷达回波信号进行预处理，获取距离谱数据，再根据自相关谱的谱峰时间间隔估计叶片转速。该专利中，计算叶片转速为中间步骤，在实施时需要对风电机组加装外辐射源雷达，根据两类波形信号进而计算谱数据，在工程实施上较为复杂，在多风机条件下如果外辐射源雷达系统参考通道中出现多径干扰，会影响谱数据的计算精度进而影响叶片转速的计算准确度。

公开号为cn110470265a的中国专利中，基于微波传感器对动叶片叶尖间隙与定时值的实时在线测量，可计算获得动叶片叶尖间隙、振动及应力等关键参数，进而通过相位峰值与转速计算模块计算相位峰值和叶片转速。该专利中需要加装精密的传感设备，具有一定的复杂度。

上述方法中，都需要加装额外的传感设备，具有一定的复杂度和高成本。计算转速为中间结果，难以为后续分析过程提供精准输入。

技术实现要素：

本申请的目的是提供一种基于音频的转速识别方法，其能够根据声音计算转速，能够脱离对数据采集和监控系统的依赖，可独立运行、效率高、运算准确，可以和后续算法配合，为后续算法输入转速信息。

本申请涉及一种基于音频的转速识别方法，包括以下步骤：

(1)通过音频采集设备采集音频数据；

(2)将采集到的音频数据转换为时频域矩阵，再转换为对数能量谱；

(3)通过滤波截取中频段的信号，将此频段范围内的对数能量谱按时域方向求均值得到中频能量向量；

(4)通过中值滤波对得到的中频能量向量进行降噪处理，选定零值位置并计算过零次数，根据过零次数估算转速；

(5)通过时间序列分析得到信号的自相关函数，通过定位自相关函数的最高点来确定信号的周期，根据周期来估算转速；

(6)计算功率谱密度，通过定位功率谱密度的峰值来确定信号的周期，根据周期来估算转速；

(7)将上述步骤(4)、步骤(5)和步骤(6)中估算得到的转速求平均值，作为最终确定的转速数值。

其中，在步骤(4)中，以中频能量向量均值为起点，上下一个标准差为波动范围定位零值位置；在步骤(5)中，对自相关函数值的峰值定位方法为：通过识别信号中的任一样本点，其直接的临近点取值均低于该点，定位到自相关函数中所有的极大值点；在步骤(7)中，对根据过零次数、自相关函数和功率谱密度确定的转速分别计算相对误差，若任两个计算得到的相对误差超过设定误差阈值，则不输出计算结果。

本申请仅需加装音频采集装置，就可以通过采集声音信号，记录风机叶片扫风声音，然后通过至少两种算法迅速计算当前叶片的实时转速，互相之间进行验证并对计算错误的样本进行过滤，整体资源消耗极低，可以在边缘设备运行。本申请还可以在无法接入scada系统时，作为补充手段迅速准确计算出风机当前的实时转速，并可以作为信息输入后续需要转速信息的故障检测算法中。

附图说明

图1为本申请实施例中的叶片扫风声音信号的示意图。

图2为本申请实施例中得到的对数能量谱的示意图。

图3为本申请实施例中截取的频率范围。

图4为本申请实施例中得到的中频能量特征。

图5为本申请实施例中对中频能量特征向量进行平滑去噪后确定的零值位置。

图6为本申请实施例中对滤波后信号计算滞后阶数的自相关函数的示意图。

图7为本申请实施例中根据极值点坐标得到滤波后信号周期的示意图。

图8为本申请实施例中对滤波后信号计算功率谱函数确定功率谱密度最大值及其坐标。

图9为本申请的转速识别方法的步骤示意图。

具体实施方式

为使本申请的目的、技术方案和优点更加清楚明白，下文中将结合附图对本申请的实施例进行详细说明。需要说明的是，在不冲突的情况下，本申请中的实施例及实施例中的特征可以相互任意组合。

参见图9，根据本申请的一种基于音频的转速识别方法，包括以下步骤：

(1)通过音频采集设备采集音频数据；

(2)然后将采集到的音频数据采用例如短时傅里叶变换的方法转换为时频域矩阵，再转换为对数能量谱，此时矩阵内元素为分贝形式表示；

(3)通过滤波截取中频段的信号，将此频段范围内对数能量谱按时域方向求均值得到中频能量向量；这是因为传感器接收风机叶片扫风声音信号为周期性信号，为了避免低频噪声和高频故障特征对信号质量的影响，可以通过带通滤波截取中频段的信号，通常为2048-4096hz之间，将此频段范围内对数能量谱按时域方向求均值得到中频能量向量；

(4)计算信号的过零次数

首先通过中值滤波对上述中频能量向量进行降噪处理，通过选定合适的邻域大小可以去除部分离群值对算法稳定性的影响。

选定零值位置用于计算过零次数，通常在稳定的扫风情况下，零值位置的定位比较灵活，为了避免叶片间不对称扫风形态的影响，通常以中频能量向量均值为起点，上下一个标准差为波动范围定位零值位置。

本实施例中，设滤波后能量信号为x，x为周期性信号，则零值位置t可表示为：

t＝mean(x)+std(x)

对信号x进行去零值处理，减去零值后得到信号x_可表示为：

x_＝x-t

过零次数的计算方法为，将序列x_分解为由x_第二个元素至末尾的子序列x_1和由x_起始至倒数第二个元素的子序列x_0，将两个子序列相乘后计算小于零元素的个数，即为信号x_的过零次数。这是因为，如果一个元素相邻两点乘积小于0，即相邻两点分别位于正负位置，说明此处经过了0点。由于叶片扫风信号的特殊性，每一支叶片转动一次会产生两个过零点，对于三支叶片风机来说，此时根据声音计算得到的估计转速rs可表示为：

rs＝sum((x1×x0)＜0)//6

(5)计算自相关函数

通过时间序列分析得到信号x的自相关函数，其中自相关函数acf用来度量时间序列中每隔k个时间单位(xt和xt-k)的观测值之间的相关，可使用不同滞后长度观测值的子序列之间的相关系数表示。因此，在选定一定的滞后阶数k后，会计算k个不同滞后长度的观测值子序列与原始序列的相关系数，得到长度为k的acf序列，一个序列和其不同滞后阶数k的历史观测值之间的相似程度可表示为：

在自相关函数值第一次达到最高点时，代表此时原序列经过了一个扫风周期，因此通过定位自相关函数最高点即可确定原始信号x的周期。对自相关函数值的峰值定位方法为：通过识别信号中的任一样本点，其直接的临近点取值均低于该点定位到自相关函数中所有的极大值点；可以设置0.8为阈值，仅保留超过阈值的最大极大值点，为最终定位的第一次最高点坐标h，如图6所示。h代表着信号x中一支叶片和临近的一支叶片扫风之间的间隔周期，对于三支叶片风机来说此时根据声音计算得到的估计转速rs可表示为：

其中，length(x)表示信号的长度；

(6)计算功率谱密度

同样基于信号x，对其进行离散傅里叶变换(dft),序列x(n),n＝1,2,…,n-1的dft是一个复数序列x(f)：

其中n表示傅里叶变换的点数，k表示傅里叶变换的第k个频谱，k/n表示每个系数捕捉到的频率。基于dft计算功率谱密度p，可表示为：

通过定位功率谱密度的峰值即可定位到原始信号的周期，功率谱密度中每个元素表示在频率k/n的能量，也即周期n/k的能量，因此功率谱密度峰值所在频率的倒数即为原始信号的周期。

对于三支叶片风机来说，此时根据声音计算得到的估计转速rs可表示为：

其中，length(x)表示信号的长度，argmax(p)表示功率谱密度最大值对应的坐标。

(7)最终确认转速

设根据过零次数、自相关函数和功率谱密度定位的转速分别为rs1，rs2，rs3，则首先对rs1，rs2，rs3计算相对误差，以5％为预先设定的误差阈值，对每两组转速计算误差百分比，若有一次计算误差超过5％则代表此组数据转速计算有误，不输出计算结果。在保证组间相对误差都小于5％时，以rs1，rs2，rs3的均值输出作为最终的估计转速，可表示为：

rs＝mean(rs1，rs2，rs3)

本申请针对已有的依据声音进行故障诊断，且依赖实时转速输入的算法，提出了无需接入scada系统即可计算转速的方法，能够降低整体系统的复杂度，使故障检测系统能够脱离风机自身数据采集和监控系统独立运行。

实施例

图1所示为通过麦克风采集风机运行时长为1分钟的叶片扫风声音信号。对信号进行短时傅里叶变换，fft点数选择8192，交叠长度选择7168，得到stfts谱图矩阵，根据对数变换公式20log10(stfts)对谱图矩阵计算得到对数能量谱，如图2所示。如图3中虚线所示，截取频率范围在2048hz-4096hz的信号，按照时域方向计算对数能量均值，得到如图4所示的中频能量特征。使用中值滤波，邻域长度选择为11，对中频能量特征向量进行平滑去噪。如图5所示，根据滤波后信号x_均值确定零值位置，计算此时信号x_过零次数为66，则根据过零次数得到的叶片转速为66/6＝11rpm。对滤波后信号x_计算滞后阶数1000的自相关函数，如图6所示。以0.8为阈值，通过峰值定位算法确定超过最大值0.8的极值点，坐标为86，如图7所示。根据极值点坐标得到滤波后信号周期为86，信号长度为3001，此时得到转速为3001/86//3＝11。对滤波后信号x_计算功率谱函数，信号时长为60秒，长度为3001，采样率为len(x_)/60＝3001/60≈50hz，分别得到功率谱密度pxx_den和对应的频率值f，如图8所示。确定功率谱密度最大值及其坐标，最大值所在坐标为argmax(pxx_den)＝35，滤波后信号x_在第35个元素时达到最大能量，在f中对应的频率为0.583，则信号x_的周期为1/f[argmax(pxx_den)]*(len(x_)/60)＝0.583*50＝85.74285714285715，此时叶片转速为3001/85.74285714285715//3＝11。由于根据过零次数、自相关函数和功率谱计算得到的转速都为11，则最终转速为11rpm。

虽然本申请所揭露的实施方式如上，但所述的内容只是为了便于理解本申请而采用的实施方式，并非用以限定本申请。任何本申请所属技术领域内的技术人员，在不脱离本申请所揭露的精神和范围的前提下，可以在实施的形式上及细节上作任何的修改与变化，但本申请的专利保护范围，仍须以所附的权利要求书所界定的范围为准。