基于非延迟代价函数和t检验的瞬时频率估计方法

摘要

本发明公开了基于非延迟代价函数和t检验的瞬时频率估计方法。该方法采用短时傅里叶变换将原始信号转换为时频谱图，采用Canny检测算法获得多条脊带，采用t检验排除每条脊带的异常值，通过叠加构建一条具有完整清晰边缘的合成脊带，采用t检验排除合成脊带的异常值，计算合成脊带的均值曲线，对均值曲线进行平滑处理，计算该平滑均值曲线在95%置信水平上的置信区间，将平滑均值曲线及其置信区间映射到目标脊线上，得到目标脊线的参考线和局部搜索区间，采用非延迟代价函数提取目标脊线。本发明适合于估计复杂多分量变频信号的瞬时频率，克服了传统方法在机械振动信号瞬时频率估计中的缺陷，估计结果的准确度和精确度高，便于工程应用。

说明书

技术领域

本发明涉及旋转机械状态监测与故障诊断领域，具体涉及基于非延迟代价函数和t检验的瞬时频率估计方法。

背景技术

由于工作环境的复杂性，旋转机械经常工作在变速条件下。瞬时频率估计是评估旋转机械运行状态及进行故障诊断的重要前提。目前常用的瞬时频率估计方法是一步代价函数法（one-step cost function）。一步代价函数法能够在局部频率范围内搜索脊点，但是局部频率范围的中心点依赖于上一个脊点的位置，这导致一步代价函数存在着延迟。此外，局部频率范围的宽度是根据经验随意设置的，且在任何时刻的宽度都是固定的，不能随着时间发生变化，这导致一步代价函数缺乏足够的自适应性。上述缺陷导致一步代价函数法在估计瞬时频率时准确度和精确度较低。

发明内容

本发明要解决的问题是针对以上不足，提出基于非延迟代价函数和t检验的瞬时频率估计方法。与现有方法相比，本发明将映射后的平滑均值曲线作为目标脊线的参考线，将映射后的置信区间作为目标脊线的局部搜索区间，因此局部频率搜索范围的中心点不依赖上一个脊点的位置，没有任何延迟，局部频率搜索范围能够自动设定，搜索带宽能够随着时间的变化而自动变化，瞬时频率估计结果的准确度和精确度高。

为解决以上技术问题，本发明提供基于非延迟代价函数和t检验的瞬时频率估计方法，其特征在于，包括以下步骤：

步骤1：采用短时傅里叶变换算法将信号x（k）（k=1, 2, …,N）转换为时频谱图，N代表信号的长度；

步骤2：从时频谱图中选取一块具有较高信噪比的局部区域，采用Canny检测算法将该局部区域转换成二值图像，二值图像包含多条脊带；局部区域是指至少包含两条脊带，信噪比大于80dB的区域；

步骤3：采用t检验算法排除每条脊带上下边缘的异常值；

步骤4：将上述多条脊带按照相互之间的运动学比例关系叠加到其中一条轮廓最完整的脊带上，构建一条具有完整清晰边缘的合成脊带；运动学比例关系是指脊带所对应的机器部件之间的传动比；

步骤5：采用t检验算法排除上述合成脊带上下边缘的异常值；

步骤6：计算上述合成脊带的均值曲线，采用五点三次平滑算法对均值曲线进行平滑处理，得到平滑均值曲线，计算该平滑均值曲线在95%置信水平上的置信区间；

步骤7：将上述平滑均值曲线及其置信区间按照平滑均值曲线与待估计目标脊线之间的运动学比例关系映射到目标脊线上；

步骤8：将映射后的平滑均值曲线作为目标脊线的参考线，将映射后的置信区间作为目标脊线的局部搜索区间；

步骤9：采用非延迟代价函数在每个时刻所对应的局部搜索区间内搜索脊点，确定每个时刻所对应的瞬时频率，最后得到整个时间区间上的瞬时频率。

进一步地，所述步骤1中短时傅里叶变换算法包括以下步骤：

1）对信号x（k）进行短时傅里叶变换：

，

TF(t, f)代表信号x（k）的短时傅里叶变换结果，t代表时间因子，f代表尺度因子，函数w(z)代表自变量为z的窗口函数；

2）计算信号x（k）的时频谱：

，

spectrogram(t, f)代表x（k）的时频谱。

进一步地，所述步骤2中Canny检测算法包括以下步骤：

1) 采用高斯滤波器对图像f(x, y)进行平滑处理，消除图像中的噪声和无关细节：

，

，

g(x, y)代表平滑后的图像，G(x, y)代表二维高斯滤波器，x代表图像的时间点，y代表图像的频率点，符号*代表卷积计算，σ代表高斯标准差；本发明中，σ=1；

2) 计算g(x, y)强度梯度的幅值和方向

，

，

，

M(x, y)代表强度梯度的幅值，θ(x, y)代表强度梯度的方向，g_x(x,>y(x,>

3) 采用非最大值镇压法消除虚假的边缘：

在梯度方向θ(x, y)上，如果非零梯度值大于相邻的两个梯度值，则该非零梯度值保持不变，否则，该非零梯度值置零；

4) 采用大小不同的两个阈值对上述消除虚假边缘后的图像进行滤波，两个阈值记为T₁和T₂，T₁<T₂，由T₁得到的图像记为I₁，由T₂得到的图像记为I₂；本发明中，T₁=0.0063，T₂=0.0156；

5) 从I₂中剔除不与强边缘连接的弱边缘，然后连接I₁和I₂中的边缘形成连续边缘。

进一步地，所述步骤3中t检验算法包括以下步骤：

1）估计信号x_n（n=1,>

，

代表样本均值，σ代表样本标准差，N代表样本长度；

2）如果，则剔除x_n；t_σ(N-1)代表标准差为σ，自由度为（N-1）的t分布。

进一步地，所述步骤9中非延迟代价函数包括以下步骤：

1）第k个时刻所对应的局部搜索区间FB_k>

，

f_k(pmc)代表映射后的平滑均值曲线在第k个时刻的值，代表映射后的平滑均值曲线置信区间在第k个时刻宽度的一半，m代表目标脊线的长度；

2）第k个时刻所对应的非延迟代价函数CF_k定义为：

,

，

f_k(i)>k范围内所取的频率值，TF(t_k,>k)代表TF(t,>k代表t在第k个时刻的值，f_k代表f在第k个时刻的值，e_k代表权重因子。

进一步地，相对误差≤0.682%，平均相对误差≤0.066%。

本发明采用以上技术方案，与现有技术相比，本发明具有以下优点：

1) 本发明具有实时性：本发明以合成脊带平滑均值曲线映射为参考线，能够即时确定当前时刻局部频率搜索范围的中心点，避免了对前一个脊点的依赖，消除了时间延迟，具有实时性。

2) 本发明具有自适应性：本发明利用合成脊带平滑均值曲线置信区间映射所提供的局部范围，能够自适应确定每个时刻所对应的局部频率搜索范围，搜索带宽能够随着时间的变化而自动变化，不需要凭借经验设置搜索带宽，从而消除了由于人为原因而产生的误差。

3) 实验结果表明：由本发明得到的瞬时频率估计值与实测值之间的最大相对误差为0.682%，平均相对误差为0.066%；与一步代价函数法的结果相比，最大相对误差降低95.84%，平均相对误差降低96.92%。

下面结合附图和实施例对本发明做进一步说明。

附图说明

附图1为本发明实施例中基于非延迟代价函数和t检验的瞬时频率估计方法的流程图；

附图2为本发明实施例中行星齿轮箱振动信号；

附图3为本发明实施例中行星齿轮箱振动信号的时频谱；

附图4为本发明实施例中从时频谱图上选取的具有高信噪比的局部区域；

附图5为本发明实施例中由Canny算法检测到的局部图像区域的边缘；

附图6为本发明实施例中采用t检验算法消除每条脊带异常点后的结果；

附图7为本发明实施例中利用脊带之间的运动学比例关系叠加而成的合成脊带（最下层的脊带即为合成脊带）；

附图8为本发明实施例中采用t检验算法消除合成脊带异常点后的结果；

附图9为本发明实施例中合成脊带的均值平滑曲线及其95%置信区间；

附图10为本发明实施例中映射得到的均值平滑曲线及其置信区间；

附图11为本发明实施例中瞬时频率估计值。

具体实施方式

实施例，如图1所示，基于非延迟代价函数和t检验的瞬时频率估计方法，包括以下步骤：

步骤1：采用短时傅里叶变换算法将信号x（k）（k=1, 2, …,N）转换为时频谱图，N代表信号的长度；

步骤2：从时频谱图中选取一块具有较高信噪比的局部区域，采用Canny检测算法将该局部区域转换成二值图像，二值图像包含多条脊带；局部区域是指至少包含两条脊带，信噪比大于80dB的区域；

步骤3：采用t检验算法排除每条脊带上下边缘的异常值；

步骤4：将上述多条脊带按照相互之间的运动学比例关系叠加到其中一条轮廓最完整的脊带上，构建一条具有完整清晰边缘的合成脊带；运动学比例关系是指脊带所对应的机器部件之间的传动比；

步骤5：采用t检验算法排除上述合成脊带上下边缘的异常值；

步骤6：计算上述合成脊带的均值曲线，采用五点三次平滑算法对均值曲线进行平滑处理，得到平滑均值曲线，计算该平滑均值曲线在95%置信水平上的置信区间；

步骤7：将上述平滑均值曲线及其置信区间按照平滑均值曲线与待估计目标脊线之间的运动学比例关系映射到目标脊线上；

步骤8：将映射后的平滑均值曲线作为目标脊线的参考线，将映射后的置信区间作为目标脊线的局部搜索区间；

步骤9：采用非延迟代价函数在每个时刻所对应的局部搜索区间内搜索脊点，确定每个时刻所对应的瞬时频率，最后得到整个时间区间上的瞬时频率。

步骤1中短时傅里叶变换算法包括以下步骤：

1）对信号x（k）进行短时傅里叶变换：

，

TF(t, f)代表信号x（k）的短时傅里叶变换结果，t代表时间因子，f代表尺度因子，函数w(z)代表自变量为z的窗口函数；

2）计算信号x（k）的时频谱：

，

spectrogram(t, f)代表x（k）的时频谱。

步骤2中Canny检测算法包括以下步骤：

1) 采用高斯滤波器对图像f(x, y)进行平滑处理，消除图像中的噪声和无关细节：

，

，

g(x, y)代表平滑后的图像，G(x, y)代表二维高斯滤波器，x代表图像的时间点，y代表图像的频率点，符号*代表卷积计算，σ代表高斯标准差；本发明中，σ=1；

2) 计算g(x, y)强度梯度的幅值和方向

，

，

，

M(x, y)代表强度梯度的幅值，θ(x, y)代表强度梯度的方向，g_x(x,>y(x,>

3) 采用非最大值镇压法消除虚假的边缘：

在梯度方向θ(x, y)上，如果非零梯度值大于相邻的两个梯度值，则该非零梯度值保持不变，否则，该非零梯度值置零；

4) 采用大小不同的两个阈值对上述消除虚假边缘后的图像进行滤波，两个阈值记为T₁和T₂，T₁<T₂，由T₁得到的图像记为I₁，由T₂得到的图像记为I₂；本发明中，T₁=0.0063，T₂=0.0156；

5) 从I₂中剔除不与强边缘连接的弱边缘，然后连接I₁和I₂中的边缘形成连续边缘。

步骤3中t检验算法包括以下步骤：

1）估计信号x_n（n=1,>

，

代表样本均值，σ代表样本标准差，N代表样本长度；

2）如果，则剔除x_n；t_σ(N-1)代表标准差为σ，自由度为（N-1）的t分布。

步骤9中非延迟代价函数包括以下步骤：

1）第k个时刻所对应的局部搜索区间FB_k>

，

f_k(pmc)代表映射后的平滑均值曲线在第k个时刻的值，代表映射后的平滑均值曲线置信区间在第k个时刻宽度的一半，m代表目标脊线的长度；

2）第k个时刻所对应的非延迟代价函数CF_k定义为：

,

，

f_k(i)>k范围内所取的频率值，TF(t_k,>k)代表TF(t,>k代表t在第k个时刻的值，f_k代表f在第k个时刻的值，e_k代表权重因子。

利用风机涡轮行星齿轮箱振动数据对本发明所述算法的性能进行验证。

振动数据从靠近行星齿轮系的齿轮箱外壳上采集，数据长度N=2736825，采样频率fs= 5000 Hz。

采集到的行星齿轮箱振动数据如图2所示。

采用短时傅里叶变换算法将图2所示的行星齿轮箱振动数据转换为时频谱图，得到的时频谱图如图3所示。

从图3所示的时频谱图上选取具有高信噪比的局部区域，得到的局部区域如图4所示。

采用Canny检测算法对如图4所示的局部区域进行边缘检测，得到的图像边缘如图5所示。

采用t检验算法消除图5中各条脊带的异常点，得到的结果如图6所示。

按照脊带之间的运动学比例关系将各条脊带叠加到其中一条轮廓最完整的脊带上，所构建的合成脊带如图7所示（最下层的脊带即为合成脊带）。

采用t检验算法消除合成脊带的异常点，结果如图8所示。

计算合成脊带的平滑均值曲线及其95%的置信区间，结果如图9所示。

按照平滑均值曲线与目标脊线之间的运动学比例关系将平滑均值曲线及其置信区间映射到目标脊线上，结果如图10所示。

采用非延迟代价函数搜索目标脊线的脊点，得到的瞬时频率曲线如图11所示。

经多次实验表明，由本发明得到的瞬时频率估计值与实测值之间的最大相对误差为0.682%，平均相对误差为0.066%，而采用一步代价函数法获得的瞬时频率估计值与实测值之间的最大相对误差为16.39%，平均相对误差为2.14%，本发明最大相对误差降低95.84%，平均相对误差降低96.92%。

根据实验结果，分析后认为：

1) 传统的一步代价函数在确定当前搜索区间的中心点时需要依赖上一个脊点的位置，存在着时间延迟现象，本发明利用映射后的平滑均值曲线作为参考线，可以即时确定当前搜索区间的中心位置，完全不依赖上一个脊点，因此具有实时性。

2) 传统的的一步代价函数法缺乏自适应性，需要人为设置搜索区间，且搜索宽度是固定的，因而不可避免地带来误差，本发明利用映射后的平滑均值曲线置信区间来自动确定局部搜索区间，搜索带宽能够随着时间的变化而自动变化，不需要人工参与，因此具有自适应性。

3) 与传统的一步代价函数法相比，本发明精确度和准确度高。

本领域技术人员应该认识到，上述的具体实施方式只是示例性的，是为了使本领域技术人员能够更好的理解本发明内容，不应理解为是对本发明保护范围的限制，只要是根据本发明技术方案所作的改进，均落入本发明的保护范围。