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

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

Instantaneous frequency estimation method based on non delay cost function and t test

The present invention discloses an instantaneous frequency estimation method based on non delay cost function and t test. This method uses the short-time Fourier transform of the original signal conversion for the spectrum, obtained using the Canny multi ridge detection algorithm, excluding outliers with each ridge by t test, to build a complete and clear edge of the ridge is synthesized with the superposition, combined with the exclusion of spinal anomalies using t test, calculation of the mean curve the synthesis of ridge belt, the mean curve smoothing, calculate the confidence interval of the mean smooth curve at 95% confidence level, will mean smooth curve and its confidence interval mapping to the target on the ridge, the ridge line of the target reference line and the local search interval, using non delay cost function to extract the target ridge. The invention is suitable for instantaneous frequency estimation in complex multi-component frequency signal, overcomes the defects of traditional method in instantaneous frequency estimation of mechanical vibration signal, accuracy and precision of the estimated results, convenient for engineering applications.

【技术实现步骤摘要】
基于非延迟代价函数和t检验的瞬时频率估计方法
本专利技术涉及旋转机械状态监测与故障诊断领域，具体涉及基于非延迟代价函数和t检验的瞬时频率估计方法。
技术介绍
由于工作环境的复杂性，旋转机械经常工作在变速条件下。瞬时频率估计是评估旋转机械运行状态及进行故障诊断的重要前提。目前常用的瞬时频率估计方法是一步代价函数法（one-stepcostfunction）。一步代价函数法能够在局部频率范围内搜索脊点，但是局部频率范围的中心点依赖于上一个脊点的位置，这导致一步代价函数存在着延迟。此外，局部频率范围的宽度是根据经验随意设置的，且在任何时刻的宽度都是固定的，不能随着时间发生变化，这导致一步代价函数缺乏足够的自适应性。上述缺陷导致一步代价函数法在估计瞬时频率时准确度和精确度较低。
技术实现思路
本专利技术要解决的问题是针对以上不足，提出基于非延迟代价函数和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)代表强度梯度的方向，gx(x,y)代表g(x,y)对x的偏导数，gy(x,y)代表g(x,y)对y的偏导数；3)采用非最大值镇压法消除虚假的边缘：在梯度方向θ(x,y)上，如果非零梯度值大于相邻的两个梯度值，则该非零梯度值保持不变，否则，该非零梯度值置零；4)采用大小不同的两个阈值对上述消除虚假边缘后的图像进行滤波，两个阈值记为T1和T2，T1<T2，由T1得到的图像记为I1，由T2得到的图像记为I2；本专利技术中，T1=0.0063，T2=0.0156；5)从I2中剔除不与强边缘连接的弱边缘，然后连接I1和I2中的边缘形成连续边缘。进一步地，所述步骤3中t检验算法包括以下步骤：1）估计信号xn（n=1,2,…,N）的标准差，，代表样本均值，σ代表样本标准差，N代表样本长度；2）如果，则剔除xn；tσ(N-1)代表标准差为σ，自由度为（N-1）的t分布。进一步地，所述步骤9中非延迟代价函数包括以下步骤：1）第k个时刻所对应的局部搜索区间FBk定义为，fk(pmc)代表映射后的平滑均值曲线在第k个时刻的值，代表映射后的平滑均值曲线置信区间在第k个时刻宽度的一半，m代表目标脊线的长度；2）第k个时刻所对应的非延迟代价函数CFk定义为：,，fk(i)代表在FBk范围内所取的频率值，TF(tk,fk)代表TF(t,f)在第k个时刻的值，tk代表t在第k个时刻的值，fk代表f在第k个时刻的值，ek代表权重因子。进一步地，相对误差≤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检验算法排除每条脊带上下边缘的异常值；步骤本文档来自技高网...

【技术保护点】
基于非延迟代价函数和t检验的瞬时频率估计方法，其特征在于，包括以下步骤：步骤1：采用短时傅里叶变换算法将信号x（k）（k=1, 2, …,N）转换为时频谱图，N代表信号的长度；步骤2：从时频谱图中选取一块具有较高信噪比的局部区域，采用Canny检测算法将该局部区域转换成二值图像，二值图像包含多条脊带；步骤3：采用t检验算法排除每条脊带上下边缘的异常值；步骤4：将上述多条脊带按照相互之间的运动学比例关系叠加到其中一条轮廓最完整的脊带上，构建一条具有完整清晰边缘的合成脊带；步骤5：采用t检验算法排除上述合成脊带上下边缘的异常值；步骤6：计算上述合成脊带的均值曲线，采用五点三次平滑算法对均值曲线进行平滑处理，得到平滑均值曲线，计算该平滑均值曲线在95%置信水平上的置信区间；步骤7：将上述平滑均值曲线及其置信区间按照平滑均值曲线与待估计目标脊线之间的运动学比例关系映射到目标脊线上；步骤8：将映射后的平滑均值曲线作为目标脊线的参考线，将映射后的置信区间作为目标脊线的局部搜索区间；步骤9：采用非延迟代价函数在每个时刻所对应的局部搜索区间内搜索脊点，确定每个时刻所对应的瞬时频率，最后得到整个时间区间上的瞬时频率。...

【技术特征摘要】
1.基于非延迟代价函数和t检验的瞬时频率估计方法，其特征在于，包括以下步骤：步骤1：采用短时傅里叶变换算法将信号x（k）（k=1,2,…,N）转换为时频谱图，N代表信号的长度；步骤2：从时频谱图中选取一块具有较高信噪比的局部区域，采用Canny检测算法将该局部区域转换成二值图像，二值图像包含多条脊带；步骤3：采用t检验算法排除每条脊带上下边缘的异常值；步骤4：将上述多条脊带按照相互之间的运动学比例关系叠加到其中一条轮廓最完整的脊带上，构建一条具有完整清晰边缘的合成脊带；步骤5：采用t检验算法排除上述合成脊带上下边缘的异常值；步骤6：计算上述合成脊带的均值曲线，采用五点三次平滑算法对均值曲线进行平滑处理，得到平滑均值曲线，计算该平滑均值曲线在95%置信水平上的置信区间；步骤7：将上述平滑均值曲线及其置信区间按照平滑均值曲线与待估计目标脊线之间的运动学比例关系映射到目标脊线上；步骤8：将映射后的平滑均值曲线作为目标脊线的参考线，将映射后的置信区间作为目标脊线的局部搜索区间；步骤9：采用非延迟代价函数在每个时刻所对应的局部搜索区间内搜索脊点，确定每个时刻所对应的瞬时频率，最后得到整个时间区间上的瞬时频率。2.根据权利要求1所述的基于非延迟代价函数和t检验的瞬时频率估计方法，其特征在于，所述步骤1中短时傅里叶变换算法包括以下步骤：1）对信号x（k）进行短时傅里叶变换：，TF(t,f)代表信号x（k）的短时傅里叶变换结果，t代表时间因子，f代表尺度因子，函数w(z)代表自变量为z的窗口函数；2）计算信号x（k）的时频谱：，spectrogram(t,f)代表x（k）的时频谱。3.根据权利要求1所述的基于非延迟代价函数和t检验的瞬时频率估计方法，其特征在于，所述步骤2中Canny检测算法包括以下步骤：1)采用高斯滤波器对图像f(x,y)进行平滑处理，消除图像中的噪声和无关细节：，，g(x,y)代表平滑后的图像，G(x,y)...

【专利技术属性】
技术研发人员：林近山，
申请(专利权)人：潍坊学院，
类型：发明
国别省市：山东,37