一种数据去噪降维的一次风机故障预警方法

1.本发明涉及风机故障预警的技术领域，尤其涉及利用数据去噪降维的聚类型mset一次风机故障预警方法。


背景技术：

2.一次风机通常由电动机、叶轮、机壳、进气箱、集流器、基座、调节控制系统等组成，目的是用来输送和干燥煤粉，并送入喷燃器，同时提供燃料燃烧初期所需的氧气。
3.一次风机作为火电厂中的重要辅机设备，运行工况复杂多变，环境恶劣，属于火电厂中易发生故障的设备，一旦发生故障会导致机组负荷返回，甚至主燃料量跳闸，影响机组安全、稳定运行，造成电厂经济损失。
4.为此提出了利用振动信号监测实现故障预警，但是一次风机测点信号众多，且信号间存在互相耦合，信号间的特性随电厂负荷变化而不同。传统的利用振动信号监测实现故障预警的方法因未充分利用一次风机信息而导致故障预警准确率低。尤其是在新能源渗透率不断增加的背景下，火电机组负荷经常大范围调整，一次风机运行工况更加复杂，故障预警难度进一步增大。
5.为了克服上述问题，提出了基于多元状态估计技术(multivariate state estimation technique，mset)的故障预警方法，充分利用测点信号的综合信息实现故障预警。但任存在以下问题：1)一次风机运行环境恶劣，测量信号中的高频干扰信号无法避免；2)一次风机监测点数目众多，传统的mset未能考虑数据间的耦合关系，导致训练矩阵维度高，计算量大；3)随着电厂信息化程度不断提高，积累了海量运行数据，传统的mset预警方法记忆矩阵构造方法未能充分考虑海量数据中关键特性的提取，导致记忆矩阵覆盖面有限，维度高，计算复杂。


技术实现要素：

6.本发明目的是提供一种含数据去噪降维的聚类型mset一次风机故障预警方法，利用测量数据间的关联性，降低了模型的复杂度和计算量，从而提高了模型运算效率，有效解决了干扰信号对故障预警准确性的影响。
7.一种数据去噪降维的一次风机故障预警方法，包括如下步骤：
8.步骤一：从电厂运行历史库中，获取与一次风机运行相关的测点历史数据；
9.步骤二：采用离散小波变换滤除步骤一测点历史数据中的高频成分；
10.步骤三：采用主成分分析对步骤二滤波后的数据进行降维；
11.步骤四：对步骤三降维后的数据基于k-means 构建记忆矩阵d；
12.步骤五：利用步骤四的记忆矩阵计算偏离度阈值；
13.步骤六：一次风机故障预警。
14.优选的是，本发明步骤一中从电厂运行历史库中，获取与一次风机运行相关的测点历史数据，所述测点包括：机组负荷、一次风机出口风压、一次风机电流、轴承温度1、轴承
温度2、轴承振动1、轴承振动2、电机非轴伸端轴承温度、电机轴伸端轴承温度、电机线圈a相温度1、电机线圈a相温度2、电机线圈b相温度1、电机线圈b相温度2、电机线圈c相温度1、电机线圈c相温度2。
15.优选的是，本发明步骤二中采用离散小波变换滤除步骤一测点历史数据中的高频成分；具体过程为：
16.步骤21，信号小波分解；选择小波基函数和小波分解层数，按照mallat算法计算小波系数w
l,k
，其中，l表示小波分解层数，k表示平移因子；
17.步骤22，阈值处理；在小波域内，信号有用部分小波系数大，而噪声部分小波系数小，通过设置阈值的方式，若小波系数绝对值小于阈值，则令对应的小波系数为零，得到下式
[0018][0019]
式中，θ表示阈值；
[0020]
阈值θ采用
‘
sqtwolog’规则计算得到，如下式
[0021][0022]
式中，n表示历史数据个数；
[0023]
步骤23，信号重构；根据阈值出力后得到的出力小波系数通过小波逆变换实现信号重构，获得各测点历史数据滤除高频噪声后的同长度有用信号。
[0024]
优选的是，本发明步骤三中采用主成分分析对步骤二滤波后的数据进行降维；
[0025]
步骤31，将经过步骤二滤波后得到的有用信号，表示为样本矩阵x
[0026][0027]
式中，n表示历史数据个数，m表示与一次风机运行相关的测点数，i为测点数据位置索引，j为测点信号索引；
[0028]
步骤32，样本矩阵x进行z-score标准化处理
[0029][0030]
其中为测点运行历史数据均值，sj为测点运行历史数据均方差；
[0031][0032]
[0033]
其中n为历史数据个数；
[0034]
z-score标准化后的样本有用数据集合矩阵记为z；
[0035][0036]
步骤33，计算样本有用数据集合矩阵z的相关系数矩阵r；根据下式计算样本有用数据集合矩阵z中任意两列序列zc＝[z
1c
ꢀ…ꢀzic
ꢀ…ꢀznc
]
t
与zd＝[z
1d
ꢀ…ꢀzid
ꢀ…ꢀznd
]
t
间相关系数；
[0037][0038]
式中，z
ic
、z
id
分别表示序列zc、zd中第i个位置对应的元素，中第i个位置对应的元素，为序列zc的平均值，为序列zd的平均值；
[0039]
由相关系数构造的相关系数矩阵r的表达式如下
[0040][0041]
步骤34，根据矩阵运算，计算相关系数矩阵r的特征值λ和特征向量v；
[0042]
步骤35，计算各参数变量的贡献率αj，对特征值λj从大到小进行排序,使得λ1≥λ2≥
…
≥λm，则各特征值对应的各参数变量的贡献率计算公式为：
[0043][0044]
步骤36，根据累积贡献率确定主元信号，主元信号是样本有用数据集合矩阵z信号的一种线性组合，主元信号中包含了样本有用数据集合矩阵z信号的大量特征，通过累积贡献率确定需要保留的主元信号数，并通过矩阵变换将样本有用数据集合矩阵z变为压缩后的样本有用数据矩阵z’，具体步骤如下：
[0045]
①
令p＝1，根据下式计算第一个主元的累积贡献率
[0046][0047]
②
令p＝p 1，根据上式计算前p个主元的累计贡献率；
[0048]
③
判断η
p
≥95％是否成立，如成立，则p为保留主元数，否则，跳转到步骤
②
；
[0049]
④
由前p个相关系数矩阵r的特征向量构造变换矩阵p，p＝[v1,...,v
p
]，由z’＝zp，将n
×
m的样本有用数据集合矩阵z压缩为n
×
p的样本有用数据矩阵z’。
[0050]
优选的是，本发明步骤四中对步骤三降维后的数据基于k-means 构建记忆矩阵d，具体过程为：
[0051]
步骤41，由样本有用数据矩阵z’组成初始记忆矩阵d
ini
；
[0052][0053]
步骤42，确定第一个质心，从初始记忆矩阵d
ini
中随机选择一行测点样本作为第一个质心c1，计算初始记忆矩阵d
ini
中任一行t∈[1,2,
…
,n]组成的测点样本z
t
＝[z’t1
,
…
,z’tp
]到c1的距离，记为d(z
t
,c1)，令质心计算器s＝1；
[0054]
步骤43，确定第二个质心，将初始记忆矩阵d
ini
中任一行t∈[1,2,
…
,n]组成的测点样本z
t
为第二个质心的概率记为p
t,2
，根据下式计算p
t,2
后，按照所得概率分布从初始记忆矩阵d
ini
中随机选择一行测点样本第二个质心c2,令质心计算器s＝2；
[0055][0056]
步骤44，确定第s个质心；令质心计算器s＝s 1,计算初始记忆矩阵d
ini
中任一行t∈[1,2,
…
,n]组成的测点样本z
t
至任一质心c
ss
，ss∈[1,2,
…
,s-1]的距离，并将样本分配给距离最近的质心c
ss
；样本分配完成后，将质心c
ss
所包含的样本集合记为c
ss
；将初始记忆矩阵d
ini
中任一行t∈[1,2,
…
,n]组成的测点样本z
t
为第s个质心的概率记为p
t,s
，根据下式计算p
t,s
后，按照所得概率分布从初始记忆矩阵d
ini
中随机选择一行测点样本第s个质心cs；
[0057][0058]
式中，{h；zh∈c
ss
}表示任意属于集合c
ss
的样本，h为集合c
ss
的样本索引；
[0059]
步骤45，重复步骤44，直到s＝nd，nd为记忆矩阵d保留的样本数；
[0060]
步骤46，计算初始记忆矩阵d
ini
中剩余行至任一质心c
ds
，ds∈[1,2,
…
,nd]的距离，并将其分配至距离最短的质心；用c’ds
表示第ds个质心所含样本集合，计算c’ds
各列均值，并将其赋值给c
ds
；
[0061]
步骤47，重复步骤46，直到质心分配不变，或达到最大迭代次数，由各个质心组成记忆矩阵d
[0062]
d＝[c1,c2,
…
,c
nd-1
,c
nd
]
t
。
[0063]
优选的是，本发明步骤五中利用步骤四的记忆矩阵计算偏离度阈值；具体过程为：
[0064]
从初始记忆矩阵d
ini
中选取nl个样本组成训练矩阵l，根据下式计算训练样本估计矩阵l
est
；
[0065]
[0066]
式中，表示矩阵间欧式距离计算符；
[0067]
计算训练样本估计矩阵l
est
与训练矩阵l间的偏离度ε
l
，称为训练偏离度
[0068][0069]
式中，nl为训练矩阵的样本数；w
pi
为第pi维主元的权重系数；l
pi
、l
est,pi
分别表示l、l
est
的第pi维元素的对应值；
[0070]
用tt表示训练偏离度的滑动窗口索引，根据滑动窗口宽度ws，则tt窗口的训练偏离度均值为序列ε
l
从tt至ws tt元素的平均值，如下式
[0071][0072]
将训练偏离度均值序列记为从中选出其最大值根据下式计算偏离度阈值e
an
[0073][0074]
其中δ为故障预警系数，δ≥0.9。
[0075]
优选的是，本发明步骤六中的一次风机故障预警，具体过程为：
[0076]
一次风机测量数据经离散小波变换、主成分分析处理后，形成观测矩阵，记为x
obs
；根据下式计算一次风机运行状态估计矩阵x
est
；
[0077][0078]
式中，表示矩阵间欧式距离计算符；
[0079]
计算状态估计矩阵与观测矩阵x
obs
间的偏离度ε
x
，称为观测偏离度
[0080][0081]
式中，ll为观测矩阵中的样本数；w
pi
为第pi维主元的权重系数；x
obs,pi
、x
est,pi
为x
obs
、x
est
的第pi维元素的对应值；
[0082]
用dd表示观测偏离度的滑动窗口索引，根据滑动窗口宽度ws，则dd窗口的观测偏离度均值为观测偏离度ε
x
从dd至ws dd元素的平均值，如下式
[0083][0084]
将观测偏离度均值序列记为从中选出其最大值
[0085]
根据与偏离度阈值e
an
间关系确定是否发出预警；若则发出故障预警。
[0086]
本发明采用上述技术方案，与现有技术相比具有如下优点：
[0087]
(1)采用离散小波变换((discrete wavelet transform,dwt)滤除测量信号中的高频干扰信号，解决现场测量信号干扰问题，有效解决了干扰信号对预警准确性的影响；
[0088]
(2)采用主成分分析(principal component analysis,pca)分析一次风机测量信号间关系，实现一次风机相关测量信号降维，有效利用了测量数据间的关联性，降低了模型的复杂度和计算量；
[0089]
(3)采用k-means 聚类方法构造mset记忆矩阵，充分挖掘了正常数据运行特征，有效降低了数据的冗杂性，增大了记忆矩阵覆盖面，提高了模型运算效率；
[0090]
(4)该一次风机故障预警方法可提前3小时35分钟实现预警，为运行检修人员故障排查提供了充足的时间。
附图说明
[0091]
图1是本发明基于多元状态估计的一次风机故障预警方法流程图。
[0092]
图2是本发明获得的原始信号及采用dwt滤除高频干扰后的有用信号示意图。
[0093]
图3是本发明主元个数与累计贡献率之间的关系图。
[0094]
图4是本发明获得的故障残差变化曲线示意图。
具体实施方式
[0095]
下面结合附图对本发明的技术方案进行详细说明：
[0096]
如图1所示，本发明的一次风机故障预警方法，包括如下步骤：
[0097]
步骤一：从电厂运行历史库中，获取与一次风机运行相关的测点历史数据，测点历史数据包括：机组负荷、一次风机出口风压、一次风机电流、轴承温度1、轴承温度2、轴承振动1、轴承振动2、电机非轴伸端轴承温度、电机轴伸端轴承温度、电机线圈a相温度1、电机线圈a相温度2、电机线圈b相温度1、电机线圈b相温度2、电机线圈c相温度1、电机线圈c相温度2等。
[0098]
步骤二：采用dwt滤除测点历史数据序列中高频成分。以第m个测点历史序列,滤波过程为：
[0099]
步骤21，信号小波分解。选择小波基函数和小波分解层数，按照mallat算法计算小波系数w
l,k
，其中，l表示小波分解层数，k表示平移因子。小波基函数选为“db1”小波基函数，分解层数为4层。
[0100]
步骤22，阈值处理。在小波域内，信号有用部分小波系数大，而噪声部分小波系数小，因此，可通过设置阈值的方式，若小波系数绝对值小于阈值，则令对应的小波系数为零，如下式
[0101][0102]
式中，θ表示阈值。
[0103]
阈值θ采用
‘
sqtwolog’规则计算得到，如下式
[0104]
[0105]
式中，n表示历史数据个数。
[0106]
步骤23，信号重构。根据阈值出力后得到的出力小波系数通过小波逆变换实现信号重构，获得各测点历史数据滤除高频噪声后的同长度有用信号。各测点历史数据及采用dwt滤除高频干扰后的有用信号如图2所示。
[0107]
步骤三：采用pca实现滤波后的有用信号降维，具体步骤为：
[0108]
步骤31，将步骤一中所述的测点信号经步骤二滤波后得到的有用信号表示为样本矩阵x
[0109][0110]
式中，n表示历史数据个数，m表示与一次风机运行相关的测点数，i为测点数据位置索引，j为测点信号索引。
[0111]
步骤32，对样本矩阵x进行z-score标准化处理
[0112][0113]
其中为测点运行历史数据均值，sj为测点运行历史数据均方差，表示为
[0114][0115][0116]
z-score标准化后的样本有用数据集合矩阵记为z
[0117][0118]
步骤33，计算样本有用数据集合矩阵z的相关系数矩阵r。根据下式计算样本有用数据集合矩阵z中任意两序列zc＝[z
1c
l z
ic
l z
nc
]
t
与zd＝[z
1d
l z
id
l z
nd
]
t
间相关系数：
[0119][0120]
式中，为序列zc的平均值，为序列zd的平均值；
[0121]
由相关系数构造相关系数矩阵r如下
[0122][0123]
步骤34，根据矩阵运算，计算相关系数矩阵r的特征值λ和特征向量v。
[0124]
步骤35，计算各参数变量的贡献率αj，对特征值λj从大到小进行排序，使得λ1≥λ2≥
…
≥λm，则各特征值对应的各参数变量的贡献率计算公式为：
[0125][0126]
步骤36，根据累积贡献率确定主元信号，主元信号是样本有用数据集合矩阵z信号的一种线性组合，主元信号中包含了样本有用数据集合矩阵z信号的大量特征，通过累积贡献率确定需要保留的主元信号数，并通过矩阵变换将样本有用数据集合矩阵z变为压缩后的样本有用数据矩阵z’，具体步骤如下：
[0127]
①
令p＝1，根据下式计算第一个主元的累积贡献率
[0128][0129]
②
令p＝p 1，根据上式计算前p个主元的累计贡献率；
[0130]
③
判断η
p
≥95％是否成立，如成立，则p为保留主元数，否则，跳转到步骤
②
；
[0131]
④
由前p个相关系数矩阵r的特征向量构造变换矩阵p，p＝[v1,...,v
p
]，由z’＝zp，将n
×
m的样本有用数据集合矩阵z压缩为n
×
p的样本有用数据矩阵z’。
[0132]
采用pca提取15维测量数据中的主元，累计贡献率随主元个数的变化规律如图3所示。从图3的主元个数与累计贡献率关系图可看出，主元数为5时，主元累计贡献率大于95％。
[0133]
⑤
由前5个相关系数矩阵r的特征向量构造变换矩阵p，p＝[v1,...,v5]，由z’＝zp，将n
×
15矩阵z压缩为n
×
5矩阵z’。
[0134]
步骤四：设置聚类中心数为2000，基于k-means 构建记忆矩阵d，步骤为：
[0135]
步骤41，由样本有用数据矩阵z’组成初始记忆矩阵d
ini
[0136][0137]
步骤42，确定第一个质心，从初始记忆矩阵d
ini
中随机选择一行测点样本作为第一个质心c1，计算初始记忆矩阵d
ini
中任一行t∈[1,2,
…
,n]组成的测点样本z
t
＝[z’t1
,
…
,z’tp
]到c1的距离，记为d(z
t
,c1)，令质心计算器s＝1；
[0138]
步骤43，确定第二个质心，将初始记忆矩阵d
ini
中任一行t∈[1,2,
…
,n]组成的测
点样本z
t
为第二个质心的概率记为p
t,2
，根据下式计算p
t,2
后，按照所得概率分布从初始记忆矩阵d
ini
中随机选择一行测点样本第二个质心c2,令质心计算器s＝2；
[0139][0140]
步骤44，确定第s个质心；令质心计算器s＝s 1,计算初始记忆矩阵d
ini
中任一行t∈[1,2,
…
,n]组成的测点样本z
t
至任一质心c
ss
，ss∈[1,2,
…
,s-1]的距离，并将样本分配给距离最近的质心c
ss
；样本分配完成后，将质心c
ss
所包含的样本集合记为c
ss
；将初始记忆矩阵d
ini
中任一行t∈[1,2,
…
,n]组成的测点样本z
t
为第s个质心的概率记为p
t,s
，根据下式计算p
t,s
后，按照所得概率分布从初始记忆矩阵d
ini
中随机选择一行测点样本第s个质心cs；
[0141][0142]
式中，{h；zh∈c
ss
}表示任意属于集合c
ss
的样本，h为集合c
ss
的样本索引。
[0143]
步骤45，重复步骤44，直到s＝nd，nd为记忆矩阵d保留的样本数；
[0144]
步骤46，计算初始记忆矩阵d
ini
中剩余行至任一质心c
ds
，ds∈[1,2,
…
,nd]的距离，并将其分配至距离最短的质心；用c’ds
表示第ds个质心所含样本集合，计算c’ds
各列均值，并将其赋值给c
ds
；
[0145]
步骤47，重复步骤46，直到质心分配不变，或达到最大迭代次数，由各个质心组成记忆矩阵d
[0146]
d＝[c1,c2,
…
,c
nd-1
,c
nd
]
t
。
[0147]
步骤五：利用步骤四的记忆矩阵计算偏离度阈值；具体过程为：
[0148]
从初始记忆矩阵d
ini
中选取nl个样本组成训练矩阵l，根据下式计算训练样本估计矩阵l
est
；
[0149][0150]
式中，表示向量间欧式距离计算符。
[0151]
计算训练样本估计矩阵l
est
与训练矩阵l间的偏离度ε
l
，称为训练偏离度
[0152][0153]
式中，nl为训练矩阵的样本数；w
pi
为第pi维主元的权重系数；l
pi
、l
est,pi
分别表示l、l
est
的第pi维元素的对应值；
[0154]
用tt表示训练偏离度的滑动窗口索引，根据滑动窗口宽度ws，则tt窗口的训练偏离度均值为序列ε
l
从tt至ws tt元素的平均值，如下式
[0155][0156]
将训练偏离度均值序列记为从中选出其最大值
根据下式计算偏离度阈值e
an
[0157][0158]
其中δ为故障预警系数，δ≥0.9。
[0159]
步骤六：一次风机故障预警，具体过程为：
[0160]
一次风机测量数据经离散小波变换、主成分分析处理后，形成观测矩阵，记为x
obs
；根据下式计算一次风机运行状态估计矩阵x
est
；
[0161][0162]
式中，表示矩阵间欧式距离计算符；
[0163]
计算状态估计矩阵与观测矩阵x
obs
间的偏离度ε
x
，称为观测偏离度
[0164][0165]
式中，ll为观测矩阵中的样本数；w
pi
为第pi维主元的权重系数；x
obs,pi
、x
est,pi
为x
obs
、x
est
的第pi维元素的对应值；
[0166]
用dd表示观测偏离度的滑动窗口索引，根据滑动窗口宽度ws，则dd窗口的观测偏离度均值为观测偏离度ε
x
从dd至ws dd元素的平均值，如下式
[0167][0168]
将观测偏离度均值序列记为从中选出其最大值
[0169]
根据与偏离度阈值e
an
间关系确定是否发出预警；若则发出故障预警。
[0170]
在某故障情况下，基于mset的一次风机残差变化趋势如图4所示，从图可知从180点开始，残差开始成增大趋势，当到245点时超过阈值，因此，提前发出故障预警。经查验，轴承振动为2.65mm/s，相比于达到3.5mm/s的报警提前了3小时35分钟，有效实现了故障预警，为运行检修人员故障排查提供了充足的时间。