CN110273818B - 一种基于轴变换粗细度分类的风机叶片结冰故障监测方法 - Google Patents

Abstract

本发明公开了一种基于轴变换粗细度分类的风机叶片结冰故障监测方法。本发明方法创新的将SCADA系统中的风速变量作为参考进行离散化，并以风速作为节点进行划数据片的操作，提出了“风速片”的概念，将大量的SCADA数据细化为不同的风速下的情况，从而有效地解决了风电数据量过大难以有效的对故障进行建模并进行及时的监测反馈的问题，有助于风场维护人员对叶片结冰的情况进行及时的诊断以及处理，从而保障了风力发电机组正常稳定的运行，同时提升了人员财产的安全保障系数。

Description

一种基于轴变换粗细度分类的风机叶片结冰故障监测方法

技术领域

本发明属于风力发电过程监测领域，特别是涉及一种基于轴变换粗细度分类的风机叶片结冰故障监测方法。

背景技术

自21世纪以来，全球人口、经济持续增长，世界能源需求增长强劲，油气资源竞争激烈，生态环境压力变大，全球气候变化备受关注。能源领域具有投资大、周期长、关联多、惯性强的发展规律。

未来二三十年，将是能源生产消费方式和能源结构调整变革的关键时期。风力和太阳能等可再生能源将快速增长，从而形成天然气、石油、煤炭、核能、可再生能源为五大支柱的能源新格局。其中，风能作为一种由空气流动所产生的能量，在我国有着广阔的开发前景。目前，根据国际上对风能资源技术开发量的评价指标，考虑了自然因素和政策因素的限制后，我国陆地70m高度层年平均风功率密度达到300W/m²以上的风能资源技术可开发量为2.6TW，70m高度层年平均风功率密度达到200W/m²以上的风能资源技术可开发量为3.6TW。不仅陆地，海上风能也有极大的开发潜力，我国沿海水深不超过50M的海上风力发电实际可装机容量约为500GW。在如此大的发展背景下，风机的装机量保持了每年10％—16％的增长率。风机的快速发展导致了故障的发生数量也大量增加，但由于风机的发电原理以及设计等原因，导致其都处于一种低温高风速的运行环境中，这也就使得风机极易出现一些利用普通手段难以监测的细微故障，这些故障的发生往往会导致风力发电机无法达到额定的工作点，例如风机叶片结冰，叶片结冰会导致风机产生无法并网的电能，造成了资源的浪费甚至对电网的稳定运行造成干扰，更有甚者，结冰程度较为严重时会导致冰块脱落，这往往会威胁到风场人员财产，同时因为风机的发电结构较为精细，因而需要较多的测点，从而导致了数据量的剧增，对此类故障的分析以及监测造成了负面影响。

综上所述，风力发电产业是一个发展极为迅速的新型产业，其工作环境复杂多变，并且由于其需要并入电网运行，因而对发电机运行的稳定性要求较高，并且需要及时的人员维护。由此可知，风力发电过程发生故障的机率较高，若不及时检测并排除故障就会导致风场无法正常运行，随之而来的就是一大笔资金使得工程重新正常运行。例如：2016年2月16日和2月20日，大唐河北乌登山风场和山西偏关水泉风场接连发生风机倒塔，造成设备损坏的事件。河北乌登山风场110号风机倒塔事故的直接原因是叶片质量出现问题，在运行中开裂，气流不平稳而引起风机的剧烈摇晃，山西偏关水泉风场项目14号风机倒塔事故的直接原因则是风机振动值严重超标，造成法兰疲劳开裂导致的风机倒塔。

针对风机的数据量大，数据特性不稳定导致的工况特征难以确定的情况，还未能有前人做出有效的处理，现阶段的风机叶片结冰等故障的监测主要依靠人工观察，或是通过二维的数据，即观察功率在对应的风速点是否达标从而判断是否结冰等其他故障，但这样往往忽略了其他变量数据带来的有效信息，从而难以有效的及时监测到风机的初期故障。

发明内容

本发明的目的在于利用风速作为主要关注变量，并通过不同风速将风电数据划分为不同工况，详细的划分了风力发电机组在不同风速下的数据特性，然后根据得到的不同的数据特性，建立不同工况下的监测模型，提高了风电数据利用效率，并提供了一种基于轴变换粗细度分类的风机叶片结冰故障监测方法。

本发明的目的通过以下技术方案实现：一种基于轴变换粗细度分类的风机叶片结冰故障监测方法，该方法包括以下步骤：

(1)离线建立监测模型，该步骤主要由以下子步骤实现：

(1.1)获取历史数据：获取风场的I个可用历史数据，每个历史数据具有J个可测变量。将I个历史数据描述为一个二维矩阵X_I(I×J)。

(1.2)将二维矩阵X_I(I×J)进行离散化处理，具体为：将二维矩阵X_I(I×J)中的“风速”变量单独提取出来，按照等宽法进行离散化处理；按照“风速”变量的离散结果对X_I(I×J)进行离散化处理，得到K个初始风速片，每个初始风速片中有I_k(k＝1，2…K)个样本以及J个可测变量，其中下标k表示属于第k个初始风速片，并有

离散化后第k个初始风速片记为X_I，k。

(1.3)提出一种步进式粗分类方法，根据变量相关特性的不同，将整个运行过程自动识别划分不同风速段，该步骤通过以下子步骤实现：

(1.3.1)风速片的标准化：对每一个初始风速片X_I，k都进行各自单独的数据标准化处理得到新的经处理后的风速片：

其中

表示X_I，k的均值，D(X_I，k)表示X_I，k的方差。

(1.3.2)对风速片

建立其PCA模型：

其中T_k(I_k×C)，

分别表示该风速片的主元空间以及投影矩阵，C表示保留的主元的个数。E_k表示第k个风速片的残差空间，t_k，i表示T_k的第i列向量，p_k，i表示

的第i行向量。

(1.3.3)计算单个控制限：计算每一个经过PCA建模后风速片的残差空间的SPE指标：

SPE_k，i＝e_k，i ^Te_k，i(i＝1，2…I_k)

其中e_k，i是第k个风速片的第i个样本的残差向量。因此，在第k个风速片上，存在I_k个SPE指标，利用该I_k个SPE指标确定第k个风速片的控制限Ctr_k。

(1.3.4)风速片步进式分类：步进式合并风速片，直到满足截止条件时，停止合并，形成一个SSPP段。

其中，风速片合并形成风速段，具体合并方式为：按照变量展开方式加入下一个风速片

使得样本数量增加I_k+1个。

所述截止条件为：若新加进的第b+n个风速片导致风速段中连续三个风速片的新控制指标Ctr_k，new都大于α×Ctr_k，则完成划分，并认为b～b+n-1片初始风速片为一个SSPP段。所述的Ctr_k，new通过如下方法计算得到：

依据步骤1.3.2和1.3.3分别计算由第b～b+n个风速片组成的风速段的PCA模型，并利用该风速段的PCA投影矩阵分别对该风速段中所包含的所有风速片进行投影处理，并分别得到各个风速片的残差空间。随后，再根据(1.3.3)所述，重新计算出每个对应风速片的SPE指标以及控制限Ctr_k，new。

(1.3.5)将第b+n片风速片作为新开始的风速片，然后重复步骤(1.3.3)～(1.3.4)，直至没有数据剩余。并记录每一个SSPP段中风速变量所处的范围，并记第m个SSPP段的风速范围为：Rwind_m。

(1.4)完成基于SSPP分段结果的GMM聚类划分：完成SSPP分段后，设S_m代表第m个SSPP段的二维数据矩阵，对S_m再次进行聚类，将其划分为G个高斯元数据类，该聚类步骤主要通过以下子步骤实现：

(1.4.1)对SSPP段进行数据白化：首先对S_m进行标准化，并进行PCA白化处理得到

保留部分主元个数使

含有S_m的90％的信息，同时得到对应SSPP段的PCA白化投影矩阵P_m，w。

(1.4.2)确定高斯混合模型GMM的初始参数：GMM如下：

其中，p(x)表示概率密度函数，π_k，θ_k，δ_k分别表示高斯混合模型的第k个高斯元的权重系数、均值以及方差。G代表一共有G个高斯元，G的值一般利用先验知识即可确定。N表示正态分布密度函数；随后，可随机确定一个π_k，θ_k，δ_k的初始参数用于后续的EM算法迭代中。

(1.4.3)利用EM迭代求得精确的参数π_k，θ_k，δ_k，该迭代主要通过以下子步骤实现：

(1.4.3.1)随机设置初始的参数π_k，0，θ_k，0，δ_k，0。其中下标0代表迭代次数。

(1.4.3.2)E-step：假设存在另一个潜变量Y＝[1，2…，G]，并假设现有的数据样本为X，第i个样本记为x_i，则Y的分布函数

如下：

其中

表示潜变量y在第i个样本以及第k个高斯元下的密度函数。其中

表示潜变量Y中的第i个样本在第k个高斯元下对应的潜变量，

表示在第k个高斯元中取到样本i以及潜变量样本i的概率，其式如下：

其中，X表示样本，θ，δ分别表示均值和方差向量，β表示协方差矩阵，J为数据表中的变量个数。

(1.4.3.3)M-step：潜变量Y可与现有变量X一同构建出一个q函数，其形式如下：

其中下标j代表迭代的次数，若j＝0则代表初值。θ_k，j+1表示第j+1次迭代的第k个高斯元的均值，δ_k，j+1表示第j+1次迭代的第k个高斯元的方差；β_k表示第k个高斯元的协方差。

随后对该q函数分别求π_j+1，θ_j+1，b_j+1的偏导：

分别求得π_j+1，θ_j+1，δ_j+1后再次带入步骤(1.4.3.2)中，按照E-step、M-step进行不断迭代，迭代结束后，记第m个SSPP段中的第i个高斯元数据为S_m，i。

(1.4.4)不同工况阈值的确定，该步骤由以下子步骤实现：

(1.4.4.1)确定每个离线样本的SPE监测阈值以及T²阈值：利用步骤(1.3.3)所述的方法确定SPE监测阈值，T²统计量的计算公式如下：

其中S表示特征值对角矩阵，t_k为T_k的行向量。

随后，每个样本都可获得一个T²统计量，从而得到T²监测阈值。将第m个SSPP段中第i个高斯元数据的SPE与T²监测阈值分别记为MSPE_m，i，MT² _m，i，则有：

MSPE_m＝[MSPE_m，1 MSPE_m，2 … MSPE_m，G]

MT² _m＝[MT² _m，1 MT² _m，2 … MT² _m，G]；

并以高斯元权重为离线建模系数：WPm＝[π_m，1 π_m，2 … π_m，G]。

(1.4.4.2)计算第m个SSPP段中的加权监测阈值：加权阈值计算公式如下：

WMSPE_m＝MSPE_m·WP_m ^T

WMT² _m＝MT² _m·WP_m ^T

获得对应SSPP段的监测阈值后，即可进行在线监测。

(2)在线监测主要由以下步骤实现：

(2.1)获得数据并确定对应工况：风场在线运行时，获得实时数据，记获得的实时数据为

将风速变量单独取出并记作

按照(1.3.5)中各个SSPP段的风速范围，确定该样本点的风速属于哪一个SSPP段。

(2.2)若有该数据属于第m个SSPP段，则利用该SSPP段中每一个高斯元分别对实时数据

进行标准化处理：

其中

为S_m，j的均值，D(S_m，j)为第S_m，j的方差。

为经过第j个高斯元处理后的数据。

现假设该SSPP段被GMM算法分为了G类，由此可以得到：

(2.3)利用GMM模型计算样本的后验概率：在步骤(1.4.3)得到高斯混合模型后，可以利用GMM模型计算新来的样本数据属于第y_i个高斯元的后验概率，此概率由以下公式求得：

(2.4)计算实时数据的统计量：分别对G个数据类{S_m，1，S_m，2，…，S_m，G}进行PCA建模，得到对应的投影矩阵，记作{P_m，1，P_m，2，…，P_m，G}，按照步骤1.3.2、1.3.3以及1.4.4.1进行计算出G个SPE统计量{SPE_m，1，SPE_m，2，…，SPE_m，G}与G个T²监测统计量{T² _m，1，T² _m，2，…，T² _m，G}，将m个SSPP段的每个新来的样本的的后验概率作为监测权重系数，记作MP_m＝[p(Y＝y_m，1|x)p(Y＝y_m，2|x)…p(Y＝y_m，G|x)]，其中x为

利用P_m，w进行白化处理得到的数据；计算出该SSPP段的MP_m。最后，由下式计算出实时数据的加权监测统计量：

OLSPE_m＝[SPE_m，1 SPE_m，2 … SPE_m，G]×MP_m ^T

OLT² _m＝[T² _m，1 T² _m，2 … T² _m，G]×MP_m ^T

(3)比较实时数据的加权监测统计量与离线的监测阈值：

采集连续6个数据的SPE与T²监测量：

A.当SPE与T²监测量均连续超限时，则判断叶片结冰。

B.当SPE与T²监测量中的其中一个均连续超限时，则判断风电机组未正常工作。

C.当SPE与T²监测量均不存在连续超限时，则判断风电机组正常工作。

D.当SPE与T²均连续等于监测阈值时，则需要等待更多的数据进入再做进一步判断。

进一步地，步骤(1.3.3)中，利用KDE方法以99％或者95％的置信区间确定出第k个风速片的控制限Ctr_k。

进一步地，步骤(1.4.4.1)中，同利用KDE方法并取99％或者95％的置信区间，获得T²监测阈值。

这里以第一个高斯元数据

为例，按照步骤1.3.2，1.3.3以及1.4.4.1进行计算：

SPE₁＝e₁ ^Te₁

t₁表示

经过投影矩阵P₁ ^T投影后的主元空间，T² ₁表示第一个高斯元的T²统计量，E₁表示第一个高斯元的残差空间，SPE₁表示残差监测统计量。

本发明的有益效果在于：本发明针对风力发电机这样的大数据量，难分工况的情况，提出了一种基于轴变换粗细度分类的风机叶片结冰故障监测方法。该方法将风速看作离散化的主要变量，然后递进式的对风力发电机组的数据进行了划分，完成了基于风速为横轴的工况划分。然后，针对每个不同的工况，建立了相应的监测模型。该方法创新的提出了以风速为划分依据的想法，细化的区分了不同风速下的风机运行情况，解决了风力发电机组数据难以进行结冰故障监测的问题。从而保证了风力发电机组的安全可靠的运行，同时达到其及时监测结冰等故障的需求。

附图说明：

图1是本发明基于轴变换粗细度分类的风机叶片结冰故障监测方法的流程图，(a)为离线建模过程流程图，(b)为在线监测流程图；

图2是本发明方法的SSPP分类展示图，(a)为功率曲线的分类展示；

图3是本发明方法的诊断图示，(a)离线检测示意图，(b)在线监测示意图；

图4是监测结果图。

具体实施方式

下面结合附图及具体实例，对本发明作进一步详细说明。

本发明以开源数据集为例。该例数据包中可用样本有393886个，变量共有28个。

应该理解，本发明不止局限于上述实例的风力发电过程，凡是熟悉本领域的技术人员在不违背本发明的前提下还可以做出等同变型或替换，这些等同的变型或替换均包含在本申请权利要求所限定的范围内。

如图1所示，本发明是一种基于轴变换粗细度分类的风机叶片结冰故障监测方法，包括以下步骤：

(1)离线建立监测模型，该步骤主要由以下子步骤实现：

(1.1)获取历史数据：设某风场建设有SCADA系统，该系统一共有J个可测变量，以及I个可用历史数据，并且每个历史数据都有对应的人工标注是否为故障的标签。随后可以将这些一共有I个时刻的历史数据，描述为一个二维矩阵X_I(I×J)，该矩阵中的J个可测变量为运行过程中可以被测量到的J个状态参数，包括风速、功率、电机转速、桨距角、桨速等。本实例中，变量参数为28，可用历史数据数目为393886。

(1.2)将二维矩阵X_I(I×J)进行离散化处理：得到了历史数据X_I后，将其中的“风速”变量单独提取出来，按照等宽法进行离散化处理；按照“风速”变量的离散结果对X_I(I×J)进行离散化处理，得到K个数据片，即将数据以“风速”变量为参考进行离散化处理，因其是以“风速”变量作为参考进行离散化，因而下称其为“初始风速片”，每个初始风速片中有I_k(k＝1,2…K)个样本以及J个可测变量，其中下标k表示属于第k个初始风速片，并有

离散化后第k个初始风速片记为X_I,k。

(1.3)进行步进式分类，该步骤通过以下子步骤实现：

(1.3.1)风速片的标准化：对每一个初始风速片X_I,k都进行各自单独的数据标准化处理得到新的经处理后的风速片(下称风速片)：

其中

表示X_I，k的均值，D(X_I，k)表示X_I，k的方差。

(1.3.2)建立PCA模型：对风速片

建立其PCA模型：

其中T_k(I_k×C)，

分别表示该风速片的主元空间以及投影矩阵；E_k表示第k个风速片的残差空间，t_k，i表示T_k的第i列向量，p_k，i表示

的第i行向量。C可以根据经验进行确定，或是确定为对每个风速片提取其90％的信息后，保留的主元个数数目出现次数最多的数值。本实施例中，保留的主元个数C为9。

(1.3.3)计算单个控制限：计算每一个经过PCA建模后风速片的残差空间的SPE指标：

SPE_k，i＝e_k，i ^Te_k，i(i＝1，2…I_k)

其中e_k，i是第k个风速片的第i个样本的残差向量。因此，在第k个风速片上，存在I_k个SPE指标，可利用KDE方法以99％或者95％的置信区间确定出第k个风速片的控制限Ctr_k。本实施例中采用99％的置信区间来确定。

(1.3.4)风速片步进式分类：步进式合并风速片，直到满足截止条件时，停止合并，形成一个SSPP段。

其中，风速片合并形成合并风速段，具体合并方式为：按照变量展开方式加入下一个风速片

使得样本数量增加I_k+1个。例如将下一个风速片

按照变量展开方式，加入到上一个风速片

中，得到的样本数量I_k+I_k+1的合并风速段，记为

其中下标k～k+1表示该合并风速段

是由第k到第k+1片之间的所有风速片组成。

所述截止条件为：若新加进的第b+n个风速片导致风速段中连续三个风速片的新控制指标Ctr_k，new都大于α×Ctr_k，则完成划分，并认为b～b+n-1片初始风速片为一个SSPP段。所述的Ctr_k，new通过如下方法计算得到：

依据步骤1.3.2和1.3.3分别计算由第b～b+n个风速片组成的合并风速段的PCA模型，并利用该风速段的PCA投影矩阵分别还原该风速段中所包含的所有风速片，并分别得到各个风速片的残差空间。随后，再根据(1.3.3)所述，重新计算出每个对应风速片的SPE指标以及控制限Ctr_k，new。

本例中，如图2(a)的功率曲线图所示，原始数据一共被划分成了7个SSPP段。记录每一个SSPP段中风速变量所处的范围，并记第m个SSPP段的风速范围为：Rwind_m

(1.4)完成基于SSPP分段结果的GMM聚类划分：完成SSPP分段后，设S_m代表第m个SSPP段的数据，对数据S_m再次进行聚类，将其划分为G个高斯元数据类，本例中，为计算方便，认为每个SSPP段中的高斯元数据类均为2个，该聚类步骤主要通过以下子步骤实现：

(1.4.1)对SSPP段进行数据白化：首先对S_m进行标准化，并进行PCA白化处理得到

保留部分主元个数使

含有S_m的90％的信息，同时得到对应SSPP段的PCA白化投影矩阵P_m，w。本例中，白化保留的主元个数为9。

(1.4.2)确定GMM的初始参数：高斯混合模型如下：

其中π_k，θ_k，δ_k表示高斯混合模型的第k个高斯元的权重系数、均值以及方差。随后，可随机确定一个π_k，θ_k，δ_k的初始参数用于后续的EM算法迭代中。

(1.4.3)利用EM迭代求得精确的参数π_k，θ_k，δ_k，该迭代主要通过以下子步骤实现：

(1.4.3.1)随机设置初始的参数π_k，0＝[0.5，0.5]，θ_k，0＝0，δ_k，0＝1。其中下标0代表迭代次数。

(1.4.3.2)E-step：本例中潜变量Y＝[1，2]，并有的数据样本为X，第i个样本记为x_i，则可以推导出Y的分布函数

如下：

其中

表示潜变量y在第i个样本以及第k个高斯元下的密度函数；

表示在第k个高斯元中取到样本i以及潜变量样本i的概率，其式如下：

(1.4.3.3)M-step：潜变量Y可与现有变量X一同构建出一个q函数，其形式如下：

其中下标j代表迭代的次数，若j＝0则代表初值。

随后对该q函数分别求π_j+1，θ_j+1，δ_j+1的偏导：

分别求得π_j+1，θ_j+1，δ_j+1后再次带入步骤(1.4.3.2)中，本例中选则迭代1000次结束。并记第m个SSPP段中的第i个高斯元数据为S_m，i

(1.4.4)不同工况阈值的确定，该步骤由以下子步骤实现：

(1.4.4.1)确定每个离线样本的SPE监测阈值以及T²阈值：利用步骤(1.3.3)所述的方法确定SPE监测阈值，T²统计量的计算公式如下：

随后，每个样本都可获得一个T²统计量，同样利用KDE方法并取99％的置信区间从而得到T²监测阈值。将第m个SSPP段中第i个高斯元数据的SPE与T²监测阈值分别记为MSPE_m，i，MT² _m，i，则有：

MSPE_m＝[MSPE_m，1 MSPE_m，2 … MSPE_m，G]

MT² _m＝[MT² _m，1 MT² _m，2 … MT² _m，G]

并以高斯元权重为离线建模系数：WPm＝[π_m，1 π_m，2 … π_m，G]。

(1.4.4.2)计算第m个SSPP段中的加权监测阈值：加权阈值计算公式如下：

WMSPE_m＝MSPE_m·WP_m ^T

WMT² _m＝MT² _m·WP_m ^T

获得对应SSPP段的监测阈值后，即可进行在线监测。

(2)在线监测主要由以下步骤实现：

(2.1)获得数据并确定对应工况：风场在线运行时，获得实时数据，记获得的实时数据为

本例数据中包含有风速、功率、电机转速、风向等变量，将风速变量单独取出并记作

按照(1.3.5)中各个SSPP段的风速范围，确定该样本点的风速属于哪一个SSPP段。

(2.2)若有该数据属于第m个SSPP段，则利用该SSPP段中每一个高斯元分别对实时数据

进行标准化处理：

其中

为S_m，j的均值，D(S_m，j)为第S_m，j的方差。

为经过第j个高斯元处理后的数据。

现假设该SSPP段被GMM算法分为了G类，由此可以得到：

(2.3)利用GMM模型计算数据的后验概率：在步骤(1.4.3)得到的高斯混合模型后，可以利用GMM模型计算每个数据属于第y_i个高斯元的后验概率，此概率由以下公式求得：

(2.4)计算实时数据的统计量：分别对G个数据类{S_m，1，S_m，2，…，S_m，G}进行PCA建模，得到对应的投影矩阵，记作{P_m，1，P_m，2，…，P_m，G}，按照步骤1.3.2、1.3.3以及1.4.4.1进行计算出G个SPE统计量{SPE_m，1，SPE_m，2，…，SPE_m，G}与G个T²监测统计量{T² _m，1，T² _m，2，…，T² _m，G}，将m个SSPP段的每个新来的样本的后验概率作为监测权重系数，记作MP_m＝[p(Y＝y_m，1|x) p(Y＝y_m，2|x) … p(Y＝y_m，G|x)]，其中x为

利用P_m，w进行白化处理得到的数据，计算出该SSPP段的MP_m。最后，由下式计算出实时数据的加权监测统计量：

OLSPE_m＝[SPE_m，1 SPE_m，2 … SPE_m，G]×MP_m ^T

OLT² _m＝[T² _m，1 T² _m，2 … T² _m，G]×MP_m ^T

(3)比较实时数据的加权监测统计量与离线的监测阈值：

采集连续6个数据的SPE与T²监测量：

A.当SPE与T²监测量均连续超限时，则判断叶片结冰。

B.当SPE与T²监测量中的其中一个均连续超限时，则判断风电机组未正常工作。

C.当SPE与T²监测量均不存在连续超限时，则判断风电机组正常工作。

D.当SPE与T²均连续等于监测阈值时，则需要等待更多的数据进入再做进一步判断。

现取用21个故障样本点进行简单展示，如图4所示，可以明确的看到，监测结果图中所有样本点的监测值均高于其对应工况下的监测阈值。

Claims (4)

1.一种基于轴变换粗细度分类的风机叶片结冰故障监测方法，其特征在于，该方法包括以下步骤：

(1)离线建立监测模型，该步骤主要由以下子步骤实现：

(1.1)获取历史数据：获取风场的I个可用历史数据，每个历史数据具有J个可测变量；将I个历史数据描述为一个二维矩阵X_I(I×J)；

(1.2)将二维矩阵X_I(I×J)进行离散化处理，具体为：将二维矩阵X_I(I×J)中的“风速”变量单独提取出来，按照等宽法进行离散化处理；按照“风速”变量的离散结果对X_I(I×J)进行离散化处理，得到K个初始风速片，每个初始风速片中有I_k(k＝1，2...K)个样本以及J个可测变量，其中下标k表示属于第k个初始风速片，并有

离散化后第k个初始风速片记为X_I，k；

(1.3)提出一种步进式粗分类方法，根据变量相关特性的不同，将整个运行过程自动识别划分不同风速段，该步骤通过以下子步骤实现：

(1.3.1)风速片的标准化：对每一个初始风速片X_I，k都进行各自单独的数据标准化处理得到新的经处理后的风速片：

其中

表示X_I，k的均值，D(X_I，k)表示X_I，k的方差；

(1.3.2)对风速片

建立其PCA模型：

其中T_k(I_k×C)，

分别表示该风速片的主元空间以及投影矩阵，C表示保留的主元的个数，E_k表示第k个风速片的残差空间，t_k，i表示T_k的第i列向量，p_k，i表示

的第i行向量；

(1.3.3)计算单个控制限：计算每一个经过PCA建模后风速片的残差空间的SPE指标：

SPE_k，i＝e_k，i ^Te_k，i(i＝1，2...I_k)

其中e_k，i是第k个风速片的第i个样本的残差向量；因此，在第k个风速片上，存在I_k个SPE指标，利用该I_k个SPE指标确定第k个风速片的控制限Ctr_k；

(1.3.4)风速片步进式分类：步进式合并风速片，直到满足截止条件时，停止合并，形成一个SSPP段；

其中，风速片合并形成风速段，具体合并方式为：按照变量展开方式加入下一个风速片

使得样本数量增加I_k+1个；

所述截止条件为：若新加进的第b+n个风速片导致风速段中连续三个风速片的新控制指标Ctr_k，new都大于α×Ctr_k，则完成划分，并认为b～b+n-1片初始风速片为一个SSPP段；所述的Ctr_k，new通过如下方法计算得到：

依据步骤1.3.2和1.3.3分别计算由第b～b+n个风速片组成的风速段的PCA模型，并利用该风速段的PCA投影矩阵分别对该风速段中所包含的所有风速片进行投影处理，并分别得到各个风速片的残差空间；随后，再根据(1.3.3)所述，重新计算出每个对应风速片的SPE指标以及控制限Ctr_k，new；

(1.3.5)将第b+n片风速片作为新开始的风速片，然后重复步骤(1.3.3)～(1.3.4)，直至没有数据剩余；并记录每一个SSPP段中风速变量所处的范围，并记第m个SSPP段的风速范围为：Rwind_m；

(1.4)完成基于SSPP分段结果的GMM聚类划分：完成SSPP分段后，设S_m代表第m个SSPP段的二维数据矩阵，对S_m再次进行聚类，将其划分为G个高斯元数据类，该聚类步骤主要通过以下子步骤实现：

(1.4.1)对SSPP段进行数据白化：首先对S_m进行标准化，并进行PCA白化处理得到

保留部分主元个数使

含有S_m的90％的信息，同时得到对应SSPP段的PCA白化投影矩阵P_m，w；

(1.4.2)确定高斯混合模型GMM的初始参数：GMM如下：

其中，p(x)表示概率密度函数，π_k，θ_k，δ_k分别表示高斯混合模型的第k个高斯元的权重系数、均值以及方差，G代表一共有G个高斯元，G的值一般利用先验知识即可确定，N表示正态分布密度函数；随后，可随机确定一个π_k，θ_k，δ_k的初始参数用于后续的EM算法迭代中；

(1.4.3)利用EM迭代求得精确的参数π_k，θ_k，δ_k，该迭代主要通过以下子步骤实现：

(1.4.3.1)随机设置初始的参数π_k，0，θ_k，0，δ_k，0，其中下标0代表迭代次数；

(1.4.3.2)E-step：假设存在另一个潜变量Y＝[1，2...，G]，并假设现有的数据样本为X，第i个样本记为x_i，则Y的分布函数

如下：

其中

表示潜变量y在第i个样本以及第k个高斯元下的密度函数，其中

表示潜变量Y中的第i个样本在第k个高斯元下对应的潜变量，

表示在第k个高斯元中取到样本i以及潜变量样本i的概率，其式如下：

其中，X表示样本，θ，δ分别表示均值和方差向量，β表示协方差矩阵，J为数据表中的变量个数；

(1.4.3.3)M-step：潜变量Y可与现有变量X一同构建出一个q函数，其形式如下：

其中下标j代表迭代的次数，若j＝0则代表初值，θ_k，j+1表示第j+1次迭代的第k个高斯元的均值，δ_k，j+1表示第j+1次迭代的第k个高斯元的方差；β_k表示第k个高斯元的协方差；

随后对该q函数分别求π_j+1，θ_j+1，b_j+1的偏导：

分别求得π_j+1，θ_j+1，δ_j+1后再次带入步骤(1.4.3.2)中，按照E-step、M-step进行不断迭代，迭代结束后，记第m个SSPP段中的第i个高斯元数据为S_m，i；

(1.4.4)不同工况阈值的确定，该步骤由以下子步骤实现：

(1.4.4.1)确定每个离线样本的SPE监测阈值以及T²阈值：利用步骤(1.3.3)所述的方法确定SPE监测阈值，T²统计量的计算公式如下：

其中S表示特征值对角矩阵，t_k为T_k的行向量；

随后，每个样本都可获得一个T²统计量，从而得到T²监测阈值；将第m个SSPP段中第i个高斯元数据的SPE与T²监测阈值分别记为MSPE_m，i，MT² _m，i，则有：

MSPE_m＝[MSPE_m，1 MSPE_m，2 ... MSPE_m，G]

MT² _m＝[MT² _m，1 MT² _m，2 ... MT² _m，G]；

并以高斯元权重为离线建模系数：WPm＝[π_m，1 π_m，2 ... π_m，G]；

(1.4.4.2)计算第m个SSPP段中的加权监测阈值：加权阈值计算公式如下：

WMSPE_m＝MSPE_m·WP_m ^T

WMT² _m＝MT² _m·WP_m ^T

获得对应SSPP段的监测阈值后，即可进行在线监测；

(2)在线监测主要由以下步骤实现：

(2.1)获得数据并确定对应工况：风场在线运行时，获得实时数据，记获得的实时数据为

将风速变量单独取出并记作

按照(1.3.5)中各个SSPP段的风速范围，确定该样本点的风速属于哪一个SSPP段；

(2.2)若有该数据属于第m个SSPP段，则利用该SSPP段中每一个高斯元分别对实时数据

进行标准化处理：

其中

为S_m，j的均值，D(S_m，j)为第S_m，j的方差，

为经过第j个高斯元处理后的数据；

现假设该SSPP段被GMM算法分为了G类，由此可以得到：

(2.3)利用GMM模型计算样本的后验概率：在步骤(1.4.3)得到高斯混合模型后，可以利用GMM模型计算新来的样本数据属于第y_i个高斯元的后验概率，此概率由以下公式求得：

(2.4)计算实时数据的统计量：分别对G个数据类{S_m，1，S_m，2，...，S_m，G}进行PCA建模，得到对应的投影矩阵，记作{P_m，1，P_m，2，...，P_m，G}，按照步骤1.3.2、1.3.3以及1.4.4.1进行计算出G个SPE统计量{SPE_m，1，SPE_m，2，...，SPE_m，G}与G个T²监测统计量{T² _m，1，T² _m，2，...，T² _m，G}，将m个SSPP段的每个新来的样本的的后验概率作为监测权重系数，记作MP_m＝[p(Y＝y_m，1|x) p(Y＝y_m，2|x) ... p(Y＝y_m，G|x)]，其中x为

利用P_m，w进行白化处理得到的数据；计算出该SSPP段的MP_m；最后，由下式计算出实时数据的加权监测统计量：

OLSPE_m＝[SPE_m，1 SPE_m，2 ... SPE_m，G]×MP_m ^T

OLT² _m＝[T² _m，1 T² _m，2 ... T² _m，G]×MP_m ^T

(3)通过比较实时数据的加权监测统计量与离线的监测阈值实现监测。

2.根据权利要求1所述的监测方法，其特征在于，采集连续6个数据的SPE与T²监测量：

A.当SPE与T²监测量均连续超限时，则判断叶片结冰；

B.当SPE与T²监测量中的其中一个均连续超限时，则判断风电机组未正常工作；

C.当SPE与T²监测量均不存在连续超限时，则判断风电机组正常工作；

D.当SPE与T²均连续等于监测阈值时，则需要等待更多的数据进入再做进一步判断。

3.根据权利要求1所述的方法，其特征在于，步骤(1.3.3)中，利用KDE方法以99％或者95％的置信区间确定出第k个风速片的控制限Ctr_k。

4.根据权利要求1所述的方法，其特征在于，步骤(1.4.4.1)中，同利用KDE方法并取99％或者95％的置信区间，获得T²监测阈值。

Images