一种基于熵权法的风机叶片故障预警方法

技术领域

本发明涉及风机故障诊断技术领域，具体而言，涉及一种基于熵权法的风机叶片故障预警方法。

背景技术

风机长期运行在露天的环境中，叶片由于较恶劣的天气环境常常发生腐蚀、开裂等故障，如果不能够及时的将故障排除，不但影响风机的发电效率，长期的运行会造成叶片的故障程度加深，甚至影响叶片的使用寿命，而随着人工智能在工业各方面的渗透，针对叶片的故障诊断涌现出多种智能算法，主要以神经网络等建立整个风场的故障识别模型为主，但是这种算法需要叶片的信号声音库足够完善，故障类型足够多，考虑到实际情况，完整数据集的建立过程困难较大，且这类方法多处于实验室阶段，实际应用案例极少。

发明内容

本发明目的在于提出一种基于熵权法的风机叶片故障预警方法，考虑到一系列实际问题，采用风机叶片正常时的信号，建立不同风速、不同风机的叶片正常时的最大可预警指标的门槛值矩阵。在对风机进行诊断预警时，计算该风机的预警指标值，并通过对比分析该风机的相同风速下的门槛值，来对风机作出是否预警的指示。进一步而言，利用小波变换从风机叶片正常的声音信号中提取能量特征矩阵，并基于熵权法建立故障预警指标的门槛值矩阵，将计算得到的待测风机的预警指标与相对应的门槛值进行比较分析，并改变小波分解的层数重复上述步骤，最后综合分析实现待测风机是否进行故障预警的任务。

本发明提供的技术方案，核心构思是通过待测信号与参考信号计算预警指标，并构建预警指标的门槛值矩阵，最后通过对比分析预警指标及其门槛值的大小综合判断风机叶片是否故障。在本发明的多个具体实施例中所述的方法包括，

步骤1.建立一数据存储单元，所述存储单元定期存储和更新一风机组在无故障情况下提取的所有风机叶片声音信号；所述声音信号包括所述风机组所有的风机处于不同风速情况下发出的声音；提取所述声音信号作为所述风机组中所有风机叶片故障预警时的备用参考信号。

步骤2.从所述风机组某风机中提取待测信号X和参考信号Y，并分别对所述待测信号X和所述参考信号Y进行小波分解并提取各自能量系数特征属性值，各个所述能量系数特征属性值组成能量系数特征属性特征维，进而分别构建待测信号特征矩阵A和参考信号特征矩阵B。

步骤3.分别提取所述待测信号特征矩阵A和所述参考信号特征矩阵B的第n列同一位置的所述能量系数特征属性值，组成m×2的特征矩阵Q且共计K个，其中n＝{1，2，3，……K}，K为所述小波分解的层数，m为所述信号划分的段数；利用熵权法分别计算K个所述特征矩阵Q中第一维和第二维所述能量系数特征维的特征权重值，分别为b

步骤4.构建所述风机组风机叶片故障预警的指标a，用以定量描述叶片的运行状态，公式如下：

步骤5.根据所述待测风机的编号和风速，遍历所述门槛值矩阵alfa，确定所述待测风机的预警门槛值a

步骤6.调整所述小波分解层数K值大小及其变化次数，并将每次计算结果记录在所述风机组中风机叶片状态信息表Tab中，判断所述风机组中风机的运行状况。

进一步的，所述步骤2中所述待测信号X和所述参考信号Y设置一定数目的采样点作为固定帧长，定义为wlen，则待测信号X和所述参考信号Y分别被分成m段信号，即

进一步的，所述步骤2中所述小波分解，利用小波能量公式计算所述参考信号Y和所述待测信号X分解后的能量系数，即E

进一步的，所述步骤3分别计算K个所述特征矩阵Q中第一维和第二维所述特征属性的特征权重值采用的熵权法，即首先根据所述特征矩阵Q中各个所述特征属性值进行归一化即

进一步的，所述步骤4中所述构建叶片故障预警指标a即

进一步的，所述步骤4所述的门槛值矩阵alfa，分别从三个时间段的声音信号存储单元中提取同风速、同一个风机的声音信号，且相互之间所述风机叶片故障预警指标a按照步骤2到4进行计算，分别用a

进一步的，所述步骤5中所述待测风机的预警门槛值a

进一步的，根据所述归一化矩阵Q′中的能量系数通过计算公式

进一步的，所述步骤2中所述待测信号X为采集所述风机组中某风机在一固定时刻风机叶片发出的一段声音信号；所述参考信号Y为遍历所述数据存储单元中的所述备用参考信号，提取与所述待测信号X对应的风机在同一风速下的声音信号，且所述待测信号X和所述参考信号Y两种信号的保持时间和帧长一致。

进一步的，所述步骤2针对所述待测信号X和所述参考信号Y两种信号中同一段声音信号分别进行K层小波分解，K是人为设定的一个常数，并利用小波能量公式提取分解后各自的能量系数，各自生成对应的待测信号特征维X1和参考信号特征维Y1；剩余段所述声音信号做上述同样的处理，各自生成待测信号特征矩阵A和参考信号特征矩阵B且大小均为m×K，其中m为所述特征矩阵行数，即声音的分段数，K为所述特征矩阵的列数，即小波分解的层数。

本发明技术方案带来的技术效果的一个方面在于，规避了采用以神经网络等建立整个风场故障识别模型过程中建立完整的数据集的困难度，提高了风机叶片故障预警的精确度和稳定性，相比于建立故障识别模型方法而言，操作简便，计算量较小，不要求完备的故障叶片数据库，更贴合于实际工程中的应用。

附图说明

图1为本发明一种基于熵权法的风机叶片故障预警方法的流程图；

图2为本发明一个实施例中风机K层小波分解波形分析示意图。

具体实施方式

实施例1

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

具体实现方法如下所述：

如图1所示，为本发明一种基于熵权法的风机叶片故障预警方法的流程图，其步骤包括：

步骤1.针对一个风场，建立风机叶片声音信号的数据存储单元，采集的声音包括风场所有风机不同风速尺度下的声音信号，将这些声音信号作为风机叶片故障预警时的备用参考信号，这个存储单元需要每隔一定的时间进行全部更新，时间段以月为单位来计量；

步骤2.在一固定时刻，对一风机的叶片进行故障诊断时，采集叶片在该时刻发出的一段声音信号作为待测信号X，设置一定数目的采样点作为固定帧长，定义为wlen，则待测信号可以被分成m段声音信号，具体的计算公式如下：

m是式(1)中结果的整数部分，wlen的采样点个数要求至少包含风机旋转一周声音采样点个数，低速轴的转速为w，单位是r/min，声音的采样频率是f，单位是Hz，则wlen需满足以下关系式：

在参考信号存储单元中寻找与待测信号同台风机，同样风速下的声音信号，且参考信号与待测信号保持时间、帧长一致；

针对两种信号的每段声音信号进行K层小波分解，K是人为设定的一个常数，并利用小波能量公式提取分解后的能量系数，计算公式如下：

E

其中x

于是每段声音信号都可以用能量系数生成的向量来表示，则对于待测信号和参考信号，可以通过上述的变换提取出能够表征两种信号本质的能量特征矩阵A和B，其中两个矩阵大大小都为m×K，即矩阵的行数为声音的分段数，矩阵的列为小波分解的层数；

步骤3.分别提取所述待测信号特征矩阵A和所述参考信号特征矩阵B的第n列同一位置的所述能量系数特征属性值，组成m×2的特征矩阵Q且共计K个，其中n＝{1，2，3，……K}，K为所述小波分解的层数，m为所述信号划分的段数；利用熵权法分别计算K个所述特征矩阵Q中第一维和第二维所述能量系数特征维的特征权重值，分别为b

所述步骤3分别计算K个矩阵Q中第一维和第二维的特征权重值采用的所述熵权法，即为根据所述特征矩阵Q中各个所述能量系数值进行归一化，即

步骤4.构建所述风机组风机叶片故障预警的指标a，用以定量描述叶片的运行状态，公式如下：

步骤5.利用声音信号的存储单元，构建预警指标a的门槛值矩阵alfa，

建立三个时间段上的声音信号的存储单元，分别从三个存储单元中取出同风速下，同一个风机的声音信号，相互之间的a指标按照步骤2到4进行计算，针对这个风机固定工况下可以计算出三个指标的值，分别用a1、a2、a3表示，则该个风机在这一风速下的故障预警指标为取三个指标中的最大值。针对不同的风机做同样的处理即可以得到预警指标a的门槛值矩阵alfa：

式(5)中，u是某一风场的风机的编号总数，v是风场不同的风速尺度数目，i的取值范围是[1，u]间的整数，j的取值范围是[1，v]的整数，

步骤6.根据待测风机的编号和风速情况，从alfa中确定该个风机的预警门槛值，当计算出来的a＜a

步骤7.改变小波分解层数K值，则步骤2到步骤6中涉及到的小波分解层数必须保持一致，通过设置K值的变化次数，即可以实现整个流程的循环的次数，将每次结果都记录在叶片状态信息表Tab中，如果表中的预警次数大于正常次数，则该个风机预警叶片故障，否则显示叶片正常，从而实现了风机叶片的故障预警。

实施例2

本实施例结合现场一实例进一步阐述一种基于熵权法的风机叶片故障预警方法的实施过程及注意事项。

某一风机现场共计5台风机组，现需要对上述5台风机组中所有的风机利用熵权法建立风机叶片故障预警系统，用于现场风机设备的管理。其具体方法如下：

1.设在某一风速下提取待测声音信号Y和参考声音信号X且将待测声音信号Y和参考声音信号X分成4段进行3层小波分解，故测得参考声音信号Y和待测声音信号X如下所示：

所述待测参考声音信号

2.分别提取所述待测信号特征矩阵A和所述参考信号特征矩阵B的第n列所述能量系数并组合成m×2的矩阵Q，其中n＝{1,2,3，……K}，共计K个矩阵Q，此处m为4，K为3，故所述特征矩阵Q分别为：

构建叶片故障预警的指标a，用以定量描述叶片的运行状态，计算公式如下所示：

3.构建预警指标a的门槛值矩阵alfa。建立三个时间段上的声音信号的存储单元并采取同风速下、同一个风机的声音信号，预警指标a按照上述步骤进行计算，针对该风机固定工况下可以计算出三个指标的值，分别用a1、a2、a3表示，则该个风机在这一风速下的故障预警指标为取三个指标中的最大值。具体过程如下，假设该风机为1号风机，分别在时间段1、时间段2和时间段3等三个时间段，在风机无故障情况下和实际运行状态下各提取了一段声音信号，且此段声音信号均被分成4段进行3层小波分解，进行一系列运算得出相对应的预警指标a1、a2、a3。测得参考声音信号Y和待测声音信号X如下所示：

上述各组声音信号对应的Q矩阵即分别提取参考声音信号X和待测声音信号Y同一列能量系数，组合成Q矩阵，如下：

对上述各个Q矩阵进行熵权法归一化处理，得出各自的归一化矩阵，具体如下：

根据计算公式

故障预警的指标a，用以定量描述叶片的运行状态，计算公式如下所示：

其他风机也做上述运算，由于实例上述已设定对现场5台风机利用熵权法建立风机叶片故障预警系统，用于现场风机设备的管理。设在某一风速下提取参考声音信号Y和待测声音信号X且将参考声音信号Y和待测声音信号X分成4段进行3层小波分解，且对风机在三个速度段进行取样，故v为3，u为5，具体如下：

根据本实施例第2步计算所述风机组风机叶片的故障预警指标0.45和上表中风机组各个风机在不同风速下的故障预警实际值做比较，即根据待测风机的编号和风速情况，从alfa中确定该个风机的预警门槛值，当计算出来的a

最后尝试调整K值大小和次数，即小波分解的层数，将每次结果都记录在叶片状态信息表Tab中，如下表所示，用以风机组叶片故障的预警。

实施例3

本实施例和实施例2的区别在于利用熵权法对归一化矩阵Q′中的特征属性权重进行二次熵权，比较下最终结果和实施例1之间的差异和联系。下表为实施例2求得的归一化矩阵Q′中的特征权重：

根据所述风机样本特征属性各个特征维的权重计算公式

此时对于故障预警的指标a利用公式

根据本实施例所述风机组风机叶片的故障预警指标0.08，参照表中风机组各个风机在不同风速下的故障预警预警实际值做比较，即根据待测风机的编号和风速情况，从alfa中确定该个风机的预警门槛值，当计算出来的a

最后尝试调整K值大小和次数，即小波分解的层数，将每次结果都记录在叶片状态信息表Tab中，如下表所示，用以风机组叶片故障的预警。

通过本实施例的计算结果可知，对基于熵权法提取所述风机组所有风机叶片能量系数矩阵特征维的二次权重相对于实施例2结果一致，即采用熵权法求取所述归一化矩阵特征属性权重或者对归一化矩阵特征维的二次权重均可以对风机组叶片故障进行预警判断。

实施例4

本实施例对某一风机组中所有风机进行的K层小波分解得到一系列离散平稳的小波进行波形分析,此实施例中会涉及Matlab软件的使用。

在对离散型波形分析所用到的函数有swt小波分解和iswt小波重构。假设某一风机场风机组共计5台风机，对所有风机提取组件所述风机特征矩阵D，并利用熵权法求取所述风机的特征属性特征权重值和一系列风机叶片故障预警值进行3层小波分解，对得出的波形进一步分析，和上述实施例2和实施例3中所述风机叶片预警值Tab表做横向比较。具体结果如下：

首先使用命令行实现一维分析，即

1.加载信号(所述风机组某一风机某段声音信号)load noisdopp；

2.设置变量s＝noisdopp；对于SWT变换，如果需要在第k层分解信号,k∈[1,2,……K]，那么原始信号需要能够平分成2^k份。所以如果原始信号的长度不满足要求，需要使用Signal Extension GUI工具或使用wextend函数来扩展它。

3.执行一层平稳小波分解[swa,swd]＝swt(s,1,'db1')；函数执行将产生1层近似和细节的系数，二者信号的长度相等，这也是平稳小波和普通小波两者不同的地方。

4.显示近似、细节小波系数

Subplot(1,2,1),plot(swa)；title('Approximation cfs')

subplot(1,2,2),plot(swd)；title('Detail cfs')

5.使用平稳小波逆变换恢复原始信号

A0＝iswt(swa,swd,'db1')；

重构的误差为err＝norm(s-A0),err＝2.1450e-014

6.从系数中构建近似和细节，输入如下代码

nulcfs＝zeros(size(swa))；

A1＝iswt(swa,nulcfs,'db1')；

D1＝iswt(nulcfs,swd,'db1')；

显示结果如下

Subplot(1,2,1),plot(A1)；title('Approximation A1')；

Subplot(1,2,2),plot(D1)；title('Detail D1')；

7.执行多层平稳小波分解

使用db1小波做3层小波分解[swa,swd]＝swt(s,3,'db1')；产生的近似系数在swa中，细节系数在swd中，而且它们有相同的长度。

8.显示近似和细节的系数

kp＝0；

for i＝1:3

subplot(3,2,kp+1),plot(swa(i,:))；

title(['Approx.cfs level',num2str(i)])

subplot(3,2,kp+2),plot(swd(i,:))；

title(['Detail cfs level',num2str(i)])

kp＝kp+2；

end

9.从系数中重建第3层的近似

mzero＝zeros(size(swd))；

A＝mzero；

A(3,:)＝iswt(swa,mzero,'db1')；

10.从系数中重建细节

D＝mzero；

for i＝1:3

swcfs＝mzero；

swcfs(i,:)＝swd(i,:)；

D(i,:)＝iswt(mzero,swcfs,'db1')；

End

11.从第3层近似和第2、3层细节中重建第1、2层近似

重建2、3层的近似

A(2,:)＝A(3,:)+D(3,:)

A(1,:)＝A(2,:)+D(2,:)；

显示第1、2、3层的近似和细节

kp＝0；

for i＝1:3

subplot(3,2,kp+1),plot(A(i,:))；

title(['Approx.level',num2str(i)]；

subplot(3,2,kp+2),plot(D(i,:))；

title(['Detail level',num2str(i)])

kp＝kp+2；

end

执行后，得到波形如图2所示，根据1、2、3层波形的分析，提取某个风机在某一风速下的趋势，其离散频次和离散幅度代表着再次瞬间风机叶片的稳定状态，映射于上述实施例2和实施例3中所述风机叶片预警值Tab表记录的相关信息，最终结果毫无二致。

实施例5

本实施例将对所述风机组中所有风机提取的所述待测信号特征矩阵A和所述参考信号特征矩阵B的第n列所述能量系数并组合成m×2的矩阵Q，第一特征维和第二特征维对应的参考声音信号Y或者待测声音信号X的先后顺序是否影响最终判断结果做进一步说明。同样假设某一风机场某风机组所有风机共计5台，对5台风机所有的风机叶片提取的特征矩阵做3层小波分解处理，其结果如下：

所述待测声音信号

此时待测信号特征矩阵特征维为特征矩阵Q的第一特征维，参考信号特征矩阵特征维为特征矩阵Q的第二特征维。而实际在求取故障预警的指标a，用以定量描述叶片的运行状态，计算公式如下所示：