一种基于脉冲频率特征模式识别的螺旋桨空化程度识别方法与流程

本发明属于信号处理和特征提取领域，尤其是涉及一种基于脉冲频率特征模式识别的螺旋桨空化程度识别方法。

背景技术：

空化是螺旋桨运行中的常见问题。空化的产生和发展不仅影响流道内速度分布，使螺旋桨工况变坏、效率降低，而且影响其动态响应，长时间的空化还可能严重损伤叶轮等过流部件。寻找适当的故障诊断、识别方法，对于有效地控制螺旋桨空化具有重要意义。

空化状态识别是螺旋桨状态监测的难点之一，目前信号处理领域常用的旋转机械故障检测方法主要有傅立叶变换、短时傅立叶变换、小波变换、第二代小波变换和多小波变换等，可以说都是基于内积原理的特征波形基函数信号分解，旨在灵活运用与特征波形相匹配的基函数去更好地处理信号，提取故障特征，从而实现故障诊断。

但是，现有技术中存在以下缺点和不足：

傅立叶变换、短时傅立叶变换、小波变换、第二代小波变换和多小波变换等故障检测的方法都建立在假设信号是平稳信号的基础上，而现实中往往是非平稳信号，从而这些检测方法都有不合理的地方，不合实际。另一方面，这些传统检测方法由于理论上的限制，很难检测到旋转机械的一些重要特征，如叶片通过频率bpf、叶片比频率brf等，有很大的局限性。

不同空化状态下，空泡的行为可以看作产生不同频段的脉冲信号，而在时域信号中这些脉冲彼此堆叠，且可能淹没在噪声中。

因此，在低信噪比条件下分离和提取脉冲频率就成为了螺旋桨空化状态识别的基础。

技术实现要素：

本发明提供了一种基于脉冲频率特征模式识别的螺旋桨空化程度识别方法，能将螺旋桨在不同空化状态下的统计特征表现出来，得到状态的估计更加准确，对进一步的信号处理和故障诊断都具有现实的指导意义。

一种基于脉冲频率特征模式识别的螺旋桨空化程度识别方法，包括：

(1)采集不同空化状态下的水下螺旋桨噪声信号；

(2)将采集的噪声信号导入程序，用快速循环平稳特征函数计算，得到循环密度谱；

(3)将得到的循环密度谱进行归一化后得到循环相干谱，然后进一步积分平均构建对数坐标下的增强包络谱；

(4)根据得到的增强包络谱判断特征频率，选择其对应时间周期的整数倍，对源数据进行改进的时域平均；

(5)对改进的时域平均后的信号进行集合经验模态分解eemd，得到对应的固有模态函数imf；

(6)采用恒虚警率检测统计不同空化状态、不同阶数的固有模态函数的脉冲频率；

(7)将上述脉冲频率作为特征矩阵，经bp神经网络训练识别后得到空化状态的判断。

相较于传统的方法，本发明中，基于循环平稳的周期估计会使时域平均的效果更准确，同时，eemd的结果可以减少不同频段的脉冲的堆叠，更准确的得到不同空化状态下的脉冲频率。

本发明的方法能将螺旋桨在不同空化状态下的脉冲频率表现出来，得到的状态估计更加贴近螺旋桨空化噪声的本质，通过得到的脉冲频率，可以实现初步的状态识别和故障诊断。

步骤(2)中，所述的循环平稳特征函数为：

其中，α为循环频率、f为频谱频率；x为待测信号；x为信号x的频谱；x^*表示的x共轭复数；e(·)为数学期望；其中，x的调幅模型的数学表达为：

其中，ai为各特征频率对应的幅度；αi为特征频率的2倍；t表示时间；n表示数目；v(t)表示背景噪声。

步骤(3)中，所述循环相干谱的数学表达式为：

其中，为有气泡工况对应的循环相干谱，为有气泡工况对应的循环密度谱，为有气泡工况对应下循环频率为0的循环密度谱。

构建对数坐标下的增强包络谱具体步骤为：

(3-1)计算增强包络谱各个循环频率对应的函数值；所述的增强包络谱的数学表达式为：

其中，为有气泡工况对应的循环相干谱；

(3-2)将函数值通过取10的对数计算得到声压级，根据得到的对数函数值范围，设置取值区间，将剩余的对数函数值赋值为对应的最值；

(3-3)根据对应的坐标点和函数值，构建对数坐标下的增强包络谱。

步骤(4)中，所述改进的时域平均的公式如下：

其中，x(n)为以时间间隔δt离散采样得到的时间信号，n为平均的周期段数目，m为一个周期中的采样点数，y(n)为改进的时域平均后得到的信号。

步骤(4)中，判断特征频率的方法采用步骤(3)中增强包络谱中幅值最高的线谱对应的频率；然后，选择其对应时间周期的整数倍作为改进的时域平均的长度。时间周期长度选择需满足两方面的需求：一是循环频率分辨率的要求，其值δα≈1/t，其中t＝m/fs，为时间周期，fs为采样频率，例如要求循环频率分辨率达到0.1hz时，所需的t大约为10s；二是平均段数的要求，在满足分辨率要求和计算效率允许的条件下，n越大降噪越明显。

步骤(5)中，集合经验模态分解eemd的计算步骤如下：

(5-1)向信号中加入正态分布白噪声；

(5-2)找出加入正态分布白噪声后的信号序列x(t)的所有极大值点和极小值点，将其用三次样条函数分别拟合为原序列的上下包络线；上下包络线的均值为m1(t)；将原数据序列减去m1(t)可得到一个减去低频的新序列h1(t)，h1(t)＝x(t)－m1(t)；

将差值与imf条件进行判断是否继续分解重复上述过程，得到第一个本征模函数分量c1(t)，它表示信号数据序列最高频率的成分；

(5-3)用x(t)减去c1(t)，得到一个去掉高频成分的新数据序列r1(t)；对r1(t)加入正态分布的白噪声后再进行上述步骤(5-2)中的分解，得到第二个本征模函数分量c2(t)；如此重复直到最后一个数据序列rn(t)不可被分解；

(5-4)重复步骤(5-1)至(5-3)，每次加入新的白噪声序列；

(5-5)将每次得到的固有模态函数imf集成均值作为最终结果。

imf条件如下：

1.信号的极值点(极大值或极小值)数目和过零点数目相等或最多相差一个；

2.由局部极大值构成的上包络线和由局部极小值构成的下包络线的平均值为零。

步骤(6)中，所述的恒虚警率检测采用上升沿计数法，上升沿判定条件为：原状态为下降，且对于时间序列的连续两点，第一点值小于阈值且第二点大于阈值；下降沿判定条件为：原状态为上升，且对于时间序列的连续两点，第一点值大于阈值且第二点小于阈值。

自适应门限检测器作为cfar检测的一种，原理是通过临近距离的单元功率来估计待检测单元杂波的平均功率，将估计的杂波平均功率乘以一个系数作为检测门限。

假设噪声服从瑞利分布，其概率密度函数为:

其中，pn＝2b²为噪声的平均功率。

则，虚警率为：

其中，η＝-pnlnpf为门限值。

因此可得，当噪声的平均功率确定时，乘以一系数得到的门限值可以保证恒虚警率，从而自适应的实现了较低虚警的统计检测。

上升沿计数法的具体步骤如下：

(6-1)输入单位时间的固有模态函数imf信号，根据信号的均值和方差设定恒定虚警率大小，由自适应门限检测器得到阈值；

(6-2)对信号前两点依次进行上升沿和下降沿判定，满足上升沿判定的，即对脉冲计数器进行“加一”操作；

(6-3)将判定的点数逐渐后移，直至遍历输入信号的所有点。

步骤(7)中，所述bp神经网络的构建和训练由matlab的工具包完成，将特征矩阵分为训练验证和测试后，分别做归一化处理；采用3层网络，隐含层节点数采用默认的10，输出层节点数为4，对应4种空化状态。4种空化状态分别为梢涡空化、片空化、泡空化、导叶空化。

与现有技术相比，本发明具有以下有益效果：

1、本发明通过改进的时域平均方法，能够将信号的周期包络信息提取一部分，并剔除其影响；其中的周期包络信息由循环平稳分析获得，相比传统检测方式，更加接近实际信号，从而解决低信噪比的问题。

2、本发明采用集合经验模态分解(eemd)，在emd的基础上增加了白噪声和平均的步骤，解决了emd带来的混叠问题；然后，利用恒虚警检测，能具有鲁棒性地统计不同imf的脉冲频率，反映不同空化阶段的特征，最终得到空化状态的估计更加准确。

附图说明

图1为本发明一种基于脉冲频率特征模式识别的螺旋桨空化程度识别方法的流程示意图；

图2为本发明实施例螺旋桨的imf2～imf11及余差图；

图3为本发明实施例螺旋桨的imf3脉冲频率图；

图4为本发明实施例螺旋桨的imf4脉冲频率图；

图5为本发明实施例螺旋桨的imf5脉冲频率图；

图6为本发明实施例螺旋桨的imf6脉冲频率图；

图7为本发明实施例螺旋桨的imf7脉冲频率图；

图8为本发明实施例螺旋桨的imf8脉冲频率图；

图9为本发明实施例bp神经网络训练结果图。

具体实施方式

下面结合附图和实施例对本发明做进一步详细描述，需要指出的是，以下所述实施例旨在便于对本发明的理解，而对其不起任何限定作用。

如图1所示，一种基于脉冲频率特征模式识别的螺旋桨空化程度识别方法，包括以下步骤：

s01，使用水听器采集水下螺旋桨的噪声。

s02，在程序中设定好相应的参数，将采集到的信号导入到程序中，计算循环密度谱：

其中：α为循环频率、f为频谱频率；x为待测信号；x为信号x的频谱；x^*表示的x共轭复数，e为数学期望。

其中，x的调幅模型的数学表达为：

其中：ai为各特征频率对应的幅度；αi为特征频率的2倍；t为表示时间；n表示数目。

s03，由s02中循环平稳特征函数计算得到的循环密度函数，根据下述公式计算增强包络谱各个循环频率对应的函数值：

s04，将函数值通过取10的对数计算等得到声压级，根据得到的对数函数值范围，设置最值限制，根据对应的坐标点和函数值，构建对数坐标下的增强包络谱。

s05，根据得到的增强包络谱，将最高峰对应的频率估计为特征频率，用其对应的时间周期，选择合适的整数倍，对源数据进行改进的时域平均。改进的时域平均的输出如下：

其中，x(n)为以时间间隔δt离散采样得到的时间信号，n为平均的周期段数目，m为一个周期中的采样点数，y(n)为改进的时域平均后得到的信号。

时间周期的估计，来自于步骤s04中的增强包络谱对应的线谱，一般采用幅值最高的线谱对应的频率作为特征频率。

然后，选择估计的特征频率的整数倍作为改进的时域平均的长度。时间周期长度的选择需要满足两方面的需求：一时循环频率分辨率的要求，其值δα≈1/t，其中t＝m/fs，为时间周期，fs为采样频率，例如要求循环频率分辨率达到0.1hz时，所需的t大约为10s；二是平均段数的要求，在满足分辨率要求和计算效率允许的条件下，n越大降噪越明显。

s06，对改进的时域平均后的信号y(n)，进行集合经验模态分解(eemd)，得到对应的固有模态函数(imf)。

imf条件如下：

1.信号的极值点(极大值或极小值)数目和过零点数目相等或最多相差一个；

2.由局部极大值构成的上包络线和由局部极小值构成的下包络线的平均值为零。

eemd的计算步骤如下：

s06-1，向信号加入正态分布白噪声；

s06-2，找出加入正态分布白噪声后的信号序列x(t)的所有极大值点和极小值点，将其用三次样条函数分别拟合为原序列的上和下包络线；上下包络线的均值为m1(t)；将原数据序列减去m1(t)可得到一个减去低频的新序列，即h1(t)＝x(t)－m1(t)；

将差值与imf条件进行判断是否继续分解重复上述过程，这样就得到第一个本征模函数分量c1(t)，它表示信号数据序列最高频率的成分；

s06-3，用x(t)减去c1(t)，得到一个去掉高频成分的新数据序列r1(t)；对r1再进行上述s06-2中的分解，得到第二个本征模函数分量c2(t)；如此重复直到最后一个数据序列rn(t)不可被分解；

s06-4，重复步骤s06-1至s06-3，每次加入新的白噪声序列；

s06-5，将每次得到的imf集成均值作为最终结果。

s07，根据s06得到的不同空化状态、不同阶数的imf，采用上升沿计数法统计脉冲频率，对比得到不同特征对应的脉冲频率。

自适应门限检测器作为cfar检测的一种，原理是通过临近距离的单元功率来估计待检测单元杂波的平均功率，将估计的杂波平均功率乘以一个系数作为检测门限。

上升沿计数法步骤如下：

s07-1，输入单位时间的imf信号，根据信号的均值和方差设定恒定虚警率大小，由自适应门限检测器得到阈值；

s07-2，对信号前两点依次进行上升沿和下降沿判定，满足上升沿判定的，即对脉冲计数器进行“加一”操作；

s07-3，将判定的点数逐渐后移，直至遍历输入信号的所有点。

s08，将上述脉冲频率作为特征矩阵，经bp神经网络训练后得到空化状态的判断。

bp神经网络的构建和训练由matlab的工具包neuralnetfitting完成，将特征矩阵分为训练验证和测试后，分别做归一化处理。采用3层网络，隐含层节点数采用默认的10，输出层节点数为4，对应4种空化状态。其余设置如训练函数等均选择默认设置。最后输入若干信号特征验证空化类别。

为了具体表现本方法在不同空化状态下螺旋桨噪声检测领域的优势和特征，采用十组不同空化状态下的七叶螺旋桨噪声进行分析。

七叶桨的转速为21转每秒左右，其某一状态下不同阶数的imf如图2所示，从上到下为imf2～imf11(对应上文中的c1～c10)及余差。可以看到不同阶数的imf均具有不同频段的脉冲。imf3脉冲频率(数量)图如图3所示，梢涡空化阶段脉冲数量明显较大，对应可能的高频特征。imf4脉冲频率图如图4所示，泡空化阶段脉冲数量明显较大，对应可能的中高频特征。imf5脉冲频率图如图5所示，导叶空化阶段脉冲数量明显较大。imf6脉冲频率图如图6所示，不同空化阶段差别不大。imf7脉冲频率图如图7所示，从梢涡到导叶空化阶段脉冲数量逐渐增长，其中泡空化阶段脉冲数量明显下降，尤其是第七个阶段，泡发展阶段。imf8脉冲频率图如图8所示，结果类似imf7，泡发展阶段脉冲数量最少。

之后，将以上的特征矩阵作为输入，得到bp神经网络的分类结果。300组数据的分组和训练结果如图9所示，其中测试集的正确率可达94％，可以认为基本达成了正确分类。

通过图3～图9，可以得到如下结论：

1.随imf阶数增大，即频率降低，脉冲频率的优势阶段逐渐向空化发展后期移动；

2.imf7的结论非常有典型性，整体线性增长，其泡空化阶段的降低可能代表此阶段的片空化被抑制，气泡溃灭数量减小；

3.bp神经网络的测试结果，验证了上述对于脉冲频率的特征对于空化状态识别的有效性；

4.由于eemd的性质，并非每阶imf的脉冲都一定有明确的物理意义，反应的统计特性未必准确。

相较于传统的方法，基于循环平稳的周期估计会是时域平均的效果更准确，而整数倍的处理使得分辨率得到了保证，同时，从处理结果的验证上来看imf7等典型特征，可以有效的判断螺旋桨的不同空化状态。

以上所述的实施例对本发明的技术方案和有益效果进行了详细说明，应理解的是以上所述仅为本发明的具体实施例，并不用于限制本发明，凡在本发明的原则范围内所做的任何修改、补充和等同替换，均应包含在本发明的保护范围之内。