CN110222393A - 基于细粒度风电发电状态划分的风机叶片结冰异常监测方法 - Google Patents

Abstract

本发明公开了一种基于细粒度风电发电状态划分的风机叶片结冰异常检测方法。在揭示风机变工况运行时各变量关系的分段性特点后，针对性地提出了一种基于慢特征提取的细粒度状态划分建模方法。针对风力发电的动态变化特性，提出了一种动静态协同的监测方法。利用风场SCADA系统采集到的数据，用以上方法对风机建立各子状态的监测模型，并离线了验证本文提出的方法检测风机出力异常的效果。充分利用了风机运行时数据的动态特性，有效地提升了检测效果，有助于风场维护人员对叶片结冰的情况进行及时的诊断以及处理，从而保障了风力发电机组正常稳定的运行，同时提升了人员财产的安全保障系数。

Description

基于细粒度风电发电状态划分的风机叶片结冰异常监测方法

技术领域

本发明属于风力发电过程监测领域，特别是涉及一种基于细粒度风电发电状态划分的风 机叶片结冰异常检测方法。

背景技术

据行业统计，2018年1-9月，全国新增风电并网容量1261万千瓦，到9月底累计 风电并网容量达到1.76亿千瓦；1-9月，全国风电发电量2676亿千瓦时，同比增长26％； 平均利用小时数1565小时，同比增加178小时；1-9月，全国弃风电量222亿千瓦时，同 比减少74亿千瓦时。

与此同时，经济性仍是制约风电发展的重要因素。与传统的化石能源电力相比，风电的 发电成本仍比较高，补贴需求和政策依赖性较强，行业发展受政策变动影响较大。而在风电 项目的开发过程中，风机能否在运转时期发挥最佳性能是衡量风场投资成败的关键因素之一。 风机极易出现一些利用普通手段难以监测的细微故障，这些故障的发生往往会导致风力发电 机无法达到额定的工作点，例如风机叶片结冰。叶片结冰会导致风机产生无法并网的电能， 造成了资源的浪费甚至对电网的稳定运行造成干扰。结冰程度较为严重时会导致冰块脱落， 这往往会威胁到风场人员财产安全。同时因为风机结构众多，测点数很多，导致数据量十分 庞大，不利于对此类故障的分析以及监测。

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

现阶段，暂时没有人针对风机运行变工况的特点，以及数据的静态与动态特性，来研究 风机的状态监测，判断其结冰故障。当前对风机结冰等故障的判断，主要依靠风电操作员的 经验，或者通过二维数据特征，通常为即风电功率曲线，来简单进行判断。但是这样做出的 结冰故障判断通常只能在严重结冰判断出来，无法及时在结冰早期就做出判断。这可能是由 于先前的方法无法有效利用风机中大量其他测点带来的完整信息，且没有考虑风机运行变工 况的特点，以及数据的静态与动态特性，从而无法对风机的整体状态形成完成、准确、及时 的监测。

发明内容

本发明的目的在于针对现有技术的不足，提供一种基于细粒度风电发电状态划分的风机 叶片结冰异常检测方法。

本发明的目的通过以下技术方案实现：一种基于细粒度风电发电状态划分的风机叶片结 冰异常检测方法，本方法的基本分析和建模单元是风速片数据矩阵，输出是沿速度方向的状 态分类，并且不同的子状态沿速度方向是连续的。通过顺序地添加新的风速片以迭代更新可 变展开建模单元，基于其与参考分布相似度的变化来检查所得模型的相似特性的变化，若特 性变化超过设定阈值，则认为不属于同一子状态。对于得出的各状态分类，可根据其统计监 测量的数据分布得监测置信限。该算法实际上是通过考察模型的改变其对监测性能的影响来 识别不同状态的。该方法包括两个部分：

(1)基于细粒度风电发电状态划分获得离线子状态，包括以下子步骤：

(1.1)获得风场正常工作的I条历史数据并将其描述为一个二维矩阵X(I×J)。矩阵中的J个 可测变量为运行过程中被测量到的J个状态参数。

(1.2)准备风速片数据矩阵：将步骤(1.1)描述的二维矩阵X(I×J)中的风速变量为参考进 行离散化处理，得到K个原始风速片。每个风速片中有I_k(k＝1,2,...,K)个样本以及J个可测 变量，其中下标k表示原始风速片的序号，并有离散化后第k个原始风速片记为 X_k。对每一个原始风速片X_k进行各自单独的数据标准化处理，得到标准化后的风速片

其中，表示X_k的均值，D(X_k)表示X_k的方差。

(1.3)对风速片执行SFA算法并获得初始风速片模型：对的协方差矩阵进行奇异值分 解得到标准正交基所张成的矩阵U_k：

其中，Λ_k是以对的协方差矩阵进行奇异值分解得到的所有标准化特征值为主对角线 的矩阵。

进一步得到白化矩阵Q_k：Q_k＝Λ_k ^-1/2U_k ^T

再对风速片进行白化，得到白化后的风速片Z_k：

对Z_k的一阶差分的协方差矩阵进行奇异值分解得到标准正交基所张成的矩阵P_k：

其中，Ω_k是以对的协方差矩阵进行奇异值分解得到的所有标准化特征值为主对角线 的矩阵。

进一步得投影矩阵W_k：W_k＝P_kQ_k

得到的Ω_k中标准化特征值数量为N，将特征值由小到大排序，保留前M个特征，得Ω_k,M， 剩余的为Ω_k,Me，Me＝N-M。对应的，保留按特征值大小排序的M个标准正交基，张成矩阵P_k,M。剩余的标准正交基张成矩阵P_k,Me。由于各风速片变量维度为J，其特征值数量相同，N＝J。

M决定方法如下：统计所有风速片低于平均变换速度的成分数量，采用出现次数最多的 成分的个数作为M。

把P_k,M、P_k,Me、Ω_k,M、Ω_k,Me作为监测模型，将风速片投影到特征空间，得风速片在特征空间的投影S_k，分为两部分计算：

S_k＝[S_k,d S_k,e]

S_k,d＝[S_k,1,S_k,2,...,S_k,M]＝P_k,MZ_k

S_k,e＝[S_k,1,S_k,2,...,S_k,Me]＝P_k,MeZ_k

其中，S_k,d包含M个慢特征[S_k,1,S_k,2,...,S_k,M]。S_k,e包含Me个剩余特征[S_k,1,S_k,2,...,S_k,Me]。

计算S_k,d中每条向量的T²监测统计值

其中，表示第k个风速片下第j条数据的T²监测统计值，s_k,d,j表示S_k,d的第j条向量。

统计离散化概率分布Pr_k,j，计算S_k,d的监测统计值与标准卡方分布的相似度Sim_k：

在标准卡方分布x轴上的与之间，取n个离散化段，n＝60；统计 在个离散化段中出现的概率Pr_k,j(j＝1,2,...,n)，得其概率分布，并统计标准卡方分布在 离散化段的概率分布为Pr_k,jr(j＝1,2,...,n)。

Sim_k的计算公式如下：

(1.4)将作为状态片Y_k,nu的初始风速片：Y_k,nu表示为Y_k,nu(Iu×J)，由当前 状态片中包含的风速片展开而成；I_u表示第u个风速片的样本数量；k为当前合成风速片中 初始风速片序数，最初等于1；nu为当前状态片内除初始风速片外的风速片个数，初始时， nu＝0。

向状态片Y_k,nu中逐一添加风速片：将nu加1，即nu＝nu+1，把风速片添加到状态 片Y_k,nu中，得到状态片Y_k,nu(Iu×J)，将状态片Y_k,nu按照(1.3)逐步计算得到状 态片的Z_k,nu、Ω_k,nu,M、Ω_k,nu,Me、P_k,nu,M和P_k,nu,Me，将P_k,nu,M、P_k,nu,Me作为状态片的监测模型，并代入状态片Y_k,nu中的所有风速片中，计算获得每个风速片在新状态片Y_k,nu中S_k,nu,u(u＝0,1,...,nu)，其中， S_k,nu,u＝[S_k,nu,u,dS_k,nu,u,e](u＝0,1,...,nu)，S_k,nu,u,d＝P_k,nu,MZ_k+u，S_k,nu,u,e＝P_k,nu,MeZ_k+u；S_k,nu,u表示使 用状态片Y_k,nu当前的监测模型对第k+u个风速片在特征空间的投影，表示第k+u个 风速片，Z_k+u表示第k+u个风速片经过如步骤(1.3)中所述白化处理后得到的矩阵。

进一步计算S_k,nu,u,d(u＝0,1,...,nu)每条向量的监测统计值统计计算获得每个风速片的监测统计值在状态片Y_k,nu中 与标准卡方分布的相似度Sim_k,nu,u(u＝0,1,...,nu)。

(1.5)将每个风速片在状态片Y_k,nu中的Sim_k,nu,u(u＝0,1,...,nu)与初始风速片 的Sim_k进行比较：若有连续三个u满足则把当前的nu作为nu^*， 继续执行步骤(1.6)；若不满足该条件，则重复执行步骤(1.4)。

其中，α是一个附加于Sim_k的常数，称为松弛因子。

(1.6)将序号为k至k+nu^*-1的状态片作为一个离线子状态；假设目前得到的是第a个离线 子状态Y_a＝Y_k,nu*-1，则该离线子状态的监测模型为P_a,M＝P_v,nu*-1,M和P_a,Me＝P_v,nu*-1,Me，且该子状 态的Z_a＝Z_k,nu*-1、Ω_a,M＝Ω_k,nu*-1,M、Ω_a,Me＝Ω_k,nu*-1,Me，进一步得S_a,d＝P_a,MZ_a，S_a,e＝P_a,MeZ_a。 计算该离线子状态Y_k,nu*-1的四个监测统计量并根据高斯核密度估计确定 这四个监测统计量对应的置信限(置信度α一般取0.95)ST1_a，ST2_a，DY1_a，DY2_a。

四个监测统计量的计算与服从分布如下：

T_a ²＝s_a,d,i ^Ts_a,d,i，服从分布：

T_a,e ²＝s_a,e,i ^Ts_a,e,i，服从分布：

S_a ²＝s_a,d,i ^TΩ_a,M ^-1s_a,d,i，服从分布：

S_a,e ²＝s_a,e,i ^TΩ_a,Me ^-1s_a,e,i，服从分布：

其中，s_a,d,i为S_a,d中的第i条向量，s_a,e,i为S_a,e中的第i条向量。表示自由度为p的卡方 分布，F_p,q表示自由度为p,q的方差比分布。

(1.7)将风速片作为新状态片的初始风速片，进行下一个离线子状态划分；新状态片 Y_k,nu中：k＝k+nu*，nu归零。递归地重复步骤(1.2)～(1.6)获得全部离线子状态。

(2)获取风场工作实时监测数据并利用步骤(1)状态划分的离线子状态对风场工作实时数 据划分状态，获得风机实时运行状态；具体包括如下步骤：

(2.1)获得数据并计算监测统计量：获取风场工作实时监测数据，包括当前时间及其之前的 Win条时序数据，每条数据包括J个状态参数。将监测数据标准化后根据风速划分成与步骤(1.7) 中获得的离线子状态相对应的λ个实时监测子状态，形成实时监测状态片Ym_v(v＝1,2,...,λ)。 根据步骤(1.3)中的相关公式计算得到Ym_v对应的Zm_v(v＝1,2,...,λ)，并根据各自对应的离线 子状态的P_a,M、P_a,Me、Ω_a,M、Ω_a,Me，计算得到Sm_v,d＝P_a,MZm_v、Sm_v,e＝P_a,MeZm_v，进一步按 照步骤(1.5)的相关公式计算得到当前实时监测状态片Ym_v(v＝1,2,...,λ)的T_v ²(v＝1,2,...,λ)、 T_v,e ²(v＝1,2,...,λ)、S_v ²(v＝1,2,...,λ)、S_v,e ²(v＝1,2,...,λ)监测统计量。

(2.2)比较实时监测统计量与离线监测阈值：根据对应的离线子状态中的控制限ST1_a，ST2_a， DY1_a，DY2_a判断实时监测子状态的监测统计量T_v ²(v＝1,2,...,λ)、T_v,e ²(v＝1,2,...,λ)、 S_v ²(v＝1,2,...,λ)、S_v,e ²(v＝1,2,...,λ)中每个数据是否超限，分别统计本次监测中各监测统计量 对应的超限个数TC_N、TeC_N、SC_N、SeC_N，及超限比例TC、TeC、SC、SeC， 其中TC、TeC作为静态指标，SC、SeC作为动态指标。

(2.3)判断风机实时运行状态：当指标超过预设的比例，这里采用10％，即判定指标异常。 如果动态指标正常，静态指标异常，则判断风机工况发生变动，但工作状态正常，判断其状 态为State1；如果动态指标正常，静态指标正常，则判断风机工况不变动，且工作状态正常， 判断其状态为State2；如果动态指标异常，静态指标也异常，则判断风机工况发生变动，且 工作状态异常，判断其状态为State3；如果动态指标异常，静态指标正常，则判断风机工况 不变动，但工作状态异常，判断其状态为State4。

本发明的有益效果在于：本发明方法利用风速作为主要关注变量，在揭示风机变工况运 行时的各变量关系分段性的特点和规律后，通过不同风速将风电数据划分为不同状态，详细 的划分了风力发电机组在不同风速下的数据特性，然后根据得到的不同的数据特性，建立不 同状态下的监测模型，提高了风电数据利用效率。该方法创新的提出了以风速为划分依据且 后续动静态协同监控的想法。有效利用了风机中大量其他测点带来的完整信息，且考虑风机 运行变工况的特点，以及数据的静态与动态特性，细化的区分了不同风速下的风机运行情况， 并动静协同，从而对风机的整体状态形成完整、准确、及时的监测。从而保证了风力发电机 组的安全可靠的运行。

附图说明：

图1是一种基于细粒度风电发电状态划分的风机叶片结冰异常检测方法的流程图，(a) 为离线建模过程流程图，(b)为在线监测流程图；

图2是本发明方法的细粒度状态划分展示图，为功率曲线的分类展示；

图3是本发明方法的在线监测图示，(a)正常状态段监测示意图，(b)异常状态段监测 示意图。

具体实施方式

针对风机有大量测点、变工况运行、数据的静态动态特性不一的情况，提出了一种细粒 度风电发电状态划分和风机叶片结冰异常监测方法。该方法在揭示风机变工况运行时的各变 量关系的分段性特点后，针对性地提出了一种基于慢特征提取的细粒度风电发电状态划分和 动静协同风机叶片结冰异常监测方法，实现离线的风机出力性能评价与刻画。其特点在于把 时序数据转化为风速片数据，然后提取慢特征进行状态划分建模。针对风力发电的动态变化 特点，提出了一种动静态协同的监测方法，对风机运行状态进行在线监测，通过判断风机当 前运行异常的情况，可以描述风机出力性能下降的情况。其特点在于监测时将数据对应到之 前划分的子模型中去，同时考虑静态与动态监测统计量来进行监测。下面结合附图及具体实 例，对本发明作进一步详细说明。

本发明以开源数据集，天池大赛的风机结冰故障预测数据集为例。包括如下变量：风速 (wind_speed)、发电机转速(generator_speed)、网侧有功功率/kw(power)、对风角/°(wind_direction)、25秒平均风向角(wind_direction_mean)、偏航位置(yaw_position)、偏 航速度(yaw_speed)、叶片1角度(pitch1_angle)、叶片2角度(pitch2_angle)、叶片3角度 (pitch3_angle)、变桨电机1温度(pitch1_moto_tmp)、变桨电机2温度(pitch2_moto_tmp)、 变桨电机3温度(pitch3_moto_tmp)、x方向加速度(acc_x)、y方向加速度(acc_y)、环境 温度(environment_tmp)、机舱温度(int_tmp)、ng5 1温度(pitch1_ng5_tmp)、ng5 2温度 (pitch2_ng5_tmp)、ng5 3温度(pitch3_ng5_tmp)、ng5 1充电器直流电流(pitch1_ng5_DC)、 ng5 2充电器直流电流(pitch2_ng5_DC)、ng5 3充电器直流电流(pitch3_ng5_DC)。

由先验知识，结合皮尔逊相关系数，可以判定pitch1_angle叶片1角度、pitch2_angle叶 片2角度、pitch3_angle叶片3角度是强相关的，可以保留其均值为pitch_angle_mean，并删 去原来的这三个变量。

pitch1_speed叶片1速度、pitch2_speed叶片2速度、pitch3_speed叶片3速度也为强相关， 可以保留其均值为pitch_speed_mean，并删去原来的这三个变量。

pitch1_moto_tmp变桨电机1温度、pitch2_moto_tmp变桨电机2温度、pitch3_moto_tmp 变桨电机3温度也为强相关，可以保留其均值为pitch_moto_tmp_mean，并删去原来的这三个 变量。

共有20个变量。因此，本实施例中提到J为20，历史数据I为350255。

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

如图1所示，本发明是一种基于细粒度风电发电状态划分的风机叶片结冰异常检测方法， 包括以下两个部分：

本发明的目的通过以下技术方案实现：一种基于细粒度风电发电状态划分的风机叶片结 冰异常检测方法，该方法包括两个部分：

(1)基于细粒度风电发电状态划分获得离线子状态，包括以下子步骤：

(1.1)获得风场正常工作的I条历史数据并将其描述为一个二维矩阵X(I×J)。矩阵中的J个 可测变量为运行过程中被测量到的J个状态参数。

(1.2)准备风速片数据矩阵：将步骤(1.1)描述的二维矩阵X(I×J)中的风速变量为参考进 行离散化处理，得到K个原始风速片。每个风速片中有I_k(k＝1,2,...,K)个样本以及J个可测 变量，其中下标k表示原始风速片的序号，并有离散化后第k个原始风速片记为 X_k。对每一个原始风速片X_k进行各自单独的数据标准化处理，得到标准化后的风速片

其中，表示X_k的均值，D(X_k)表示X_k的方差。

(1.3)对风速片执行SFA算法并获得初始风速片模型：对的协方差矩阵进行奇异值分 解得到标准正交基所张成的矩阵U_k：

其中，Λ_k是以对的协方差矩阵进行奇异值分解得到的所有标准化特征值为主对角线 的矩阵。

进一步得到白化矩阵Q_k：Q_k＝Λ_k ^-1/2U_k ^T

再对风速片进行白化，得到白化后的风速片Z_k：

对Z_k的一阶差分的协方差矩阵进行奇异值分解得到标准正交基所张成的矩阵P_k：

其中，Ω_k是以对的协方差矩阵进行奇异值分解得到的所有标准化特征值为主对角线 的矩阵。

进一步得投影矩阵W_k：W_k＝P_kQ_k

得到的Ω_k中标准化特征值数量为N，将特征值由小到大排序，保留前M个特征，得Ω_k,M， 剩余的为Ω_k,Me，Me＝N-M。对应的，保留按特征值大小排序的M个标准正交基，张成矩阵P_k,M。剩余的标准正交基张成矩阵P_k,Me。由于各风速片变量维度为J，其特征值数量相同，N＝J。

M决定方法如下：统计所有风速片低于平均变换速度的成分数量，采用出现次数最多的 成分的个数作为M。

把P_k,M、P_k,Me、Ω_k,M、Ω_k,Me作为监测模型，将风速片投影到特征空间，得风速片在特征空间的投影S_k，分为两部分计算：

S_k＝[S_k,d S_k,e]

S_k,d＝[S_k,1,S_k,2,...,S_k,M]＝P_k,MZ_k

S_k,e＝[S_k,1,S_k,2,...,S_k,Me]＝P_k,MeZ_k

其中，S_k,d包含M个慢特征[S_k,1,S_k,2,...,S_k,M]。S_k,e包含Me个剩余特征[S_k,1,S_k,2,...,S_k,Me]。

计算S_k,d中每条向量的T²监测统计值

其中，表示第k个风速片下第j条数据的T²监测统计值，s_k,d,j表示S_k,d的第j条向量。

统计离散化概率分布Pr_k,j，计算S_k,d的监测统计值与标准卡方分布的相似度Sim_k：

在标准卡方分布x轴上的与之间，取n个离散化段，n＝60；统计 在个离散化段中出现的概率Pr_k,j(j＝1,2,...,n)，得其概率分布，并统计标准卡方分布在 离散化段的概率分布为Pr_k,jr(j＝1,2,...,n)。

Sim_k的计算公式如下：

(1.4)将作为状态片Y_k,nu的初始风速片：Y_k,nu表示为Y_k,nu(Iu×J)，由当前 状态片中包含的风速片展开而成；I_u表示第u个风速片的样本数量；k为当前合成风速片中 初始风速片序数，最初等于1；nu为当前状态片内除初始风速片外的风速片个数，初始时， nu＝0。

向状态片Y_k,nu中逐一添加风速片：将nu加1，即nu＝nu+1，把风速片添加到状态 片Y_k,nu中，得到状态片Y_k,nu(Iu×J)，将状态片Y_k,nu按照(1.3)逐步计算得到状 态片的Z_k,nu、Ω_k,nu,M、Ω_k,nu,Me、P_k,nu,M和P_k,nu,Me，将P_k,nu,M、P_k,nu,Me作为状态片的监测模型，并代入状态片Y_k,nu中的所有风速片中，计算获得每个风速片在新状态片Y_k,nu中S_k,nu,u(u＝0,1,...,nu)，其中， S_k,nu,u＝[S_k,nu,u,dS_k,nu,u,e](u＝0,1,...,nu)，S_k,nu,u,d＝P_k,nu,MZ_k+u，S_k,nu,u,e＝P_k,nu,MeZ_k+u；S_k,nu,u表示使 用状态片Y_k,nu当前的监测模型对第k+u个风速片在特征空间的投影，表示第k+u个 风速片，Z_k+u表示第k+u个风速片经过如步骤(1.3)中所述白化处理后得到的矩阵。

进一步计算S_k,nu,u,d(u＝0,1,...,nu)每条向量的监测统计值统计计算获得每个风速片的监测统计值在状态片Y_k,nu中 与标准卡方分布的相似度Sim_k,nu,u(u＝0,1,...,nu)。

(1.5)将每个风速片在状态片Y_k,nu中的Sim_k,nu,u(u＝0,1,...,nu)与初始风速片 的Sim_k进行比较：若有连续三个u满足则把当前的nu作为nu^*， 继续执行步骤(1.6)；若不满足该条件，则重复执行步骤(1.4)。

其中，α是一个附加于Sim_k的常数，称为松弛因子。

(1.6)将序号为k至k+nu^*-1的状态片作为一个离线子状态；假设目前得到的是第a个离线 子状态Y_a＝Y_k,nu*-1，则该离线子状态的监测模型为P_a,M＝P_v,nu*-1,M和P_a,Me＝P_v,nu*-1,Me，且该子状 态的Z_a＝Z_k,nu*-1、Ω_a,M＝Ω_k,nu*-1,M、Ω_a,Me＝Ω_k,nu*-1,Me，进一步得S_a,d＝P_a,MZ_a，S_a,e＝P_a,MeZ_a。 计算该离线子状态Y_k,nu*-1的四个监测统计量并根据高斯核密度估计确定 这四个监测统计量对应的置信限(置信度α一般取0.95)ST1_a，ST2_a，DY1_a，DY2_a。

四个监测统计量的计算与服从分布如下：

T_a ²＝s_a,d,i ^Ts_a,d,i，服从分布：

T_a,e ²＝s_a,e,i ^Ts_a,e,i，服从分布：

S_a ²＝s_a,d,i ^TΩ_a,M ^-1s_a,d,i，服从分布：

S_a,e ²＝s_a,e,i ^TΩ_a,Me ^-1s_a,e,i，服从分布：

其中，s_a,d,i为S_a,d中的第i条向量，s_a,e,i为S_a,e中的第i条向量。表示自由度为p的卡方 分布，F_p,q表示自由度为p,q的方差比分布。

(1.7)将风速片作为新状态片的初始风速片，进行下一个离线子状态划分；新状态片 Y_k,nu中：k＝k+nu*，nu归零。递归地重复步骤(1.2)～(1.6)获得全部离线子状态。

本实施例中，分为了7个子状态，如图2所示。

(2)获取风场工作实时监测数据并利用步骤(1)状态划分的离线子状态对风场工作实时数 据划分状态，获得风机实时运行状态；具体包括如下步骤：

(2.1)获得数据并计算监测统计量：获取风场工作实时监测数据，包括当前时间及其之前的 Win条时序数据，每条数据包括J个状态参数。将监测数据标准化后根据风速划分成与步骤(1.7) 中获得的离线子状态相对应的λ个实时监测子状态，形成实时监测状态片Ym_v(v＝1,2,...,λ)。 根据步骤(1.3)中的相关公式计算得到Ym_v对应的Zm_v(v＝1,2,...,λ)，并根据各自对应的离线 子状态的P_a,M、P_a,Me、Ω_a,M、Ω_a,Me，计算得到Sm_v,d＝P_a,MZm_v、Sm_v,e＝P_a,MeZm_v，进一步按 照步骤(1.5)的相关公式计算得到当前实时监测状态片Ym_v(v＝1,2,...,λ)的T_v ²(v＝1,2,...,λ)、 T_v,e ²(v＝1,2,...,λ)、S_v ²(v＝1,2,...,λ)、S_v,e ²(v＝1,2,...,λ)监测统计量。

(2.2)比较实时监测统计量与离线监测阈值：根据对应的离线子状态中的控制限ST1_a，ST2_a， DY1_a，DY2_a判断实时监测子状态的监测统计量T_v ²(v＝1,2,...,λ)、T_v,e ²(v＝1,2,...,λ)、 S_v ²(v＝1,2,...,λ)、S_v,e ²(v＝1,2,...,λ)中每个数据是否超限，分别统计本次监测中各监测统计量 对应的超限个数TC_N、TeC_N、SC_N、SeC_N，及超限比例TC、TeC、SC、SeC， 其中TC、TeC作为静态指标，SC、SeC作为动态指标。

(2.3)判断风机实时运行状态：当指标超过预设的比例，这里采用10％，即判定指标异常。 如果动态指标正常，静态指标异常，则判断风机工况发生变动，但工作状态正常，判断其状 态为State1；如果动态指标正常，静态指标正常，则判断风机工况不变动，且工作状态正常， 判断其状态为State2；如果动态指标异常，静态指标也异常，则判断风机工况发生变动，且 工作状态异常，判断其状态为State3；如果动态指标异常，静态指标正常，则判断风机工况 不变动，但工作状态异常，判断其状态为State4。

(2.4)判断得到结果后，即可指导风机的运行与故障维修。

选取风机正常工作的某段时间段和风机异常工作，即叶片结冰的某段时间段，对监测效 果进行验证。若State3到达高位，对应严重异常，若State4到达高位，对应微小异常。其中 Abnormal到达高位，对应人工标注的结冰异常。如图3(a)所示，正常段监测时会出现少量 误报，在可接受范围之内；如图3(b)所示，异常段监测时，能够提前预测结冰异常的出现， 对风机的整体状态形成了完整、准确、及时的监测。有助于指导风机的故障维修，从而保障 风力发电机组安全可靠运行。

本方法创新的提出了以风速为划分依据且后续动静态协同监控的想法。有效利用了风机 中大量其他测点带来的完整信息，且考虑风机运行变工况的特点，以及数据的静态与动态特 性，细化的区分了不同风速下的风机运行情况，并动静协同，从而对风机的整体状态形成完 整、准确、及时的监测。从而保证了风力发电机组的安全可靠的运行。

Claims (1)

1.基于细粒度风电发电状态划分的风机叶片结冰异常检测方法，其特征在于，包括如下步骤：

(1)基于细粒度风电发电状态划分获得离线子状态，包括以下子步骤：

(1.1)获得风场正常工作的I条历史数据并将其描述为一个二维矩阵X(I×J)。矩阵中的J个可测变量为运行过程中被测量到的J个状态参数。

(1.2)准备风速片数据矩阵：将步骤(1.1)描述的二维矩阵X(I×J)中的风速变量为参考进行离散化处理，得到K个原始风速片。每个风速片中有I_k(k＝1,2,...,K)个样本以及J个可测变量，其中下标k表示原始风速片的序号，并有离散化后第k个原始风速片记为X_k。对每一个原始风速片X_k进行各自单独的数据标准化处理，得到标准化后的风速片

其中，表示X_k的均值，D(X_k)表示X_k的方差。

(1.3)对风速片执行SFA算法并获得初始风速片模型：对的协方差矩阵进行奇异值分解得到标准正交基所张成的矩阵U_k：

其中，Λ_k是以对的协方差矩阵进行奇异值分解得到的所有标准化特征值为主对角线的矩阵。

进一步得到白化矩阵Q_k：Q_k＝Λ_k ^-1/2U_k ^T

再对风速片进行白化，得到白化后的风速片Z_k：

对Z_k的一阶差分的协方差矩阵进行奇异值分解得到标准正交基所张成的矩阵P_k：

其中，Ω_k是以对的协方差矩阵进行奇异值分解得到的所有标准化特征值为主对角线的矩阵。

进一步得投影矩阵W_k：W_k＝P_kQ_k

得到的Ω_k中标准化特征值数量为N，将特征值由小到大排序，保留前M个特征，得Ω_k,M，剩余的为Ω_k,Me，Me＝N-M。对应的，保留按特征值大小排序的M个标准正交基，张成矩阵P_k,M。剩余的标准正交基张成矩阵P_k,Me。由于各风速片变量维度为J，其特征值数量相同，N＝J。

M决定方法如下：统计所有风速片低于平均变换速度的成分数量，采用出现次数最多的成分的个数作为M。

把P_k,M、P_k,Me、Ω_k,M、Ω_k,Me作为监测模型，将风速片投影到特征空间，得风速片在特征空间的投影S_k，分为两部分计算：

S_k＝[S_k,d S_k,e]

S_k,d＝[S_k,1,S_k,2,...,S_k,M]＝P_k,MZ_k

S_k,e＝[S_k,1,S_k,2,...,S_k,Me]＝P_k,MeZ_k

其中，S_k,d包含M个慢特征[S_k,1,S_k,2,...,S_k,M]。S_k,e包含Me个剩余特征[S_k,1,S_k,2,...,S_k,Me]。

计算S_k,d中每条向量的T²监测统计值

其中，表示第k个风速片下第j条数据的T²监测统计值，s_k,d,j表示S_k,d的第j条向量。

统计离散化概率分布Pr_k,j，计算S_k,d的监测统计值与标准卡方分布的相似度Sim_k：

在标准卡方分布x轴上的与之间，取n个离散化段，n＝60；统计在个离散化段中出现的概率Pr_k,j(j＝1,2,...,n)，得其概率分布，并统计标准卡方分布在离散化段的概率分布为Pr_k,jr(j＝1,2,...,n)。

Sim_k的计算公式如下：

(1.4)将作为状态片Y_k,nu的初始风速片：Y_k,nu表示为Y_k,nu(Iu×J)，由当前状态片中包含的风速片展开而成；I_u表示第u个风速片的样本数量；k为当前合成风速片中初始风速片序数，最初等于1；nu为当前状态片内除初始风速片外的风速片个数，初始时，nu＝0。

向状态片Y_k,nu中逐一添加风速片：将nu加1，即nu＝nu+1，把风速片添加到状态片Y_k,nu中，得到状态片Y_k,nu(Iu×J)，将状态片Y_k,nu按照(1.3)逐步计算得到状态片的Z_k,nu、Ω_k,nu,M、Ω_k,nu,Me、P_k,nu,M和P_k,nu,Me，将P_k,nu,M、P_k,nu,Me作为状态片的监测模型，并代入状态片Y_k,nu中的所有风速片中，计算获得每个风速片在新状态片Y_k,nu中S_k,nu,u(u＝0,1,...,nu)，其中，S_k,nu,u＝[S_k,nu,u,dS_k,nu,u,e](u＝0,1,...,nu)，S_k,nu,u,d＝P_k,nu,MZ_k+u，S_k,nu,u,e＝P_k,nu,MeZ_k+u；S_k,nu,u表示使用状态片Y_k,nu当前的监测模型对第k+u个风速片在特征空间的投影，表示第k+u个风速片，Z_k+u表示第k+u个风速片经过如步骤(1.3)中所述白化处理后得到的矩阵。

进一步计算S_k,nu,u,d(u＝0,1,...,nu)每条向量的监测统计值统计计算获得每个风速片的监测统计值在状态片Y_k,nu中与标准卡方分布的相似度Sim_k,nu,u(u＝0,1,...,nu)。

(1.5)将每个风速片在状态片Y_k,nu中的Sim_k,nu,u(u＝0,1,...,nu)与初始风速片的Sim_k进行比较：若有连续三个u满足则把当前的nu作为nu^*，继续执行步骤(1.6)；若不满足该条件，则重复执行步骤(1.4)。

其中，α是一个附加于Sim_k的常数，称为松弛因子。

(1.6)将序号为k至k+nu^*-1的状态片作为一个离线子状态；假设目前得到的是第a个离线子状态则该离线子状态的监测模型为和且该子状态的进一步得S_a,d＝P_a,MZ_a，S_a,e＝P_a,MeZ_a。计算该离线子状态的四个监测统计量并根据高斯核密度估计确定这四个监测统计量对应的置信限(置信度α一般取0.95)ST1_a，ST2_a，DY1_a，DY2_a。

四个监测统计量的计算与服从分布如下：

T_a ²＝s_a,d,i ^Ts_a,d,i，服从分布：

T_a,e ²＝s_a,e,i ^Ts_a,e,i，服从分布：

S_a ²＝s_a,d,i ^TΩ_a,M ^-1s_a,d,i，服从分布：

S_a,e ²＝s_a,e,i ^TΩ_a,Me ^-1s_a,e,i，服从分布：

其中，s_a,d,i为S_a,d中的第i条向量，s_a,e,i为S_a,e中的第i条向量。表示自由度为p的卡方分布，F_p,q表示自由度为p,q的方差比分布。

(1.7)将风速片作为新状态片的初始风速片，进行下一个离线子状态划分；新状态片Y_k,nu中：k＝k+nu*，nu归零。递归地重复步骤(1.2)～(1.6)获得全部离线子状态。

(2)获取风场工作实时监测数据并利用步骤(1)状态划分的离线子状态对风场工作实时数据划分状态，获得风机实时运行状态；具体包括如下步骤：

(2.1)获得数据并计算监测统计量：获取风场工作实时监测数据，包括当前时间及其之前的Win条时序数据，每条数据包括J个状态参数。将监测数据标准化后根据风速划分成与步骤(1.7)中获得的离线子状态相对应的λ个实时监测子状态，形成实时监测状态片Ym_v(v＝1,2,...,λ)。根据步骤(1.3)中的相关公式计算得到Ym_v对应的Zm_v(v＝1,2,...,λ)，并根据各自对应的离线子状态的P_a,M、P_a,Me、Ω_a,M、Ω_a,Me，计算得到Sm_v,d＝P_a,MZm_v、Sm_v,e＝P_a,MeZm_v，进一步按照步骤(1.5)的相关公式计算得到当前实时监测状态片Ym_v(v＝1,2,...,λ)的T_v ²(v＝1,2,...,λ)、T_v,e ²(v＝1,2,...,λ)、S_v ²(v＝1,2,...,λ)、S_v,e ²(v＝1,2,...,λ)监测统计量。

(2.2)比较实时监测统计量与离线监测阈值：根据对应的离线子状态中的控制限ST1_a，ST2_a，DY1_a，DY2_a判断实时监测子状态的监测统计量T_v ²(v＝1,2,...,λ)、T_v,e ²(v＝1,2,...,λ)、S_v ²(v＝1,2,...,λ)、S_v,e ²(v＝1,2,...,λ)中每个数据是否超限，分别统计本次监测中各监测统计量对应的超限个数TC_N、TeC_N、SC_N、SeC_N，及超限比例TC、TeC、SC、SeC，其中TC、TeC作为静态指标，SC、SeC作为动态指标。

(2.3)判断风机实时运行状态：当指标超过预设的比例，这里采用10％，即判定指标异常。如果动态指标正常，静态指标异常，则判断风机工况发生变动，但工作状态正常，判断其状态为State1；如果动态指标正常，静态指标正常，则判断风机工况不变动，且工作状态正常，判断其状态为State2；如果动态指标异常，静态指标也异常，则判断风机工况发生变动，且工作状态异常，判断其状态为State3；如果动态指标异常，静态指标正常，则判断风机工况不变动，但工作状态异常，判断其状态为State4。