CN103743563A - 一种基于温度数据的风机齿轮箱子空间故障预测方法 - Google Patents

Abstract

本发明公开了一种基于温度数据的风机齿轮箱子空间故障预测方法，包括以下步骤：A1温度数据的预处理；利用回归分析方法对温度数据进行单步预测，得到实际值和预测值之间的差值，称为残差，将残差作为随机状态空间模型的观测量Y;A2随机状态空间模型的识别；A3齿轮箱的故障预测；本发明通过分析温度数据所反映的齿轮箱的内在特征，可以在温度不高的故障初期阶段发出预警信号，以减轻齿轮箱的损伤，避免发生不可逆的故障。

Description

一种基于温度数据的风机齿轮箱子空间故障预测方法

技术领域

本发明涉及的是一种基于温度数据的风机齿轮箱子空间故障预测方法。

背景技术

在风电机组众多故障类型中，虽然齿轮箱的故障率相对较低，但其故障造成的停机时间是最长的、经济损失最大。目前针对风机齿轮箱的故障检测方法主要有：基于振动信号的方法、基于噪声信号的方法、油液分析法、声发射检测技术等，但现在风机上针对以上方法的传感器还未普及，现阶段国内的风电机组已经实现对齿轮箱的温度监测，因此，从现有技术和经济的角度考虑，直接利用温度信号来实现对齿轮箱的故障预测有其特定的优势。大多数故障的发生是一个渐变的过程，如轴承与齿轮的故障、润滑油不充足以及部件松动连接都会引起温度的异常升高，在故障发生前都会出现一些早期征兆，这些征兆与正常运行时的状态特征信号是有区别的。通过对实时的温度数据进行分析，在故障的早期阶段识别出这些异常特征，就有可能提前预测出齿轮箱未来的运行状态，并采取相应的应对措施。

发明内容

本发明所要解决的技术问题是针对现有技术的不足提供一种基于温度数据的风机齿轮箱子空间故障预测方法。

本发明的技术方案如下：

一种基于温度数据的风机齿轮箱子空间故障预测方法，包括以下步骤：

A1温度数据的预处理；利用回归分析方法对温度数据进行单步预测，得到实际值和预测值之间的差值，称为残差，将残差作为随机状态空间模型的观测量Y;

A11多元线性回归模型

多元线性回归模型的一般形式如下：

y＝β₀+β₁x₁+…+β_px_p+ε   (1)

式(1)中，β₀,β₁,…,β_p是未知参数，β₀为回归常数，β₁,…,β_p为回归系数；y为因变量；x₁,x₂,…,x_p为自变量，这里是与因变量相关的监测量；ε为随机误差。若已知参数的估计值

则可以实现温度的预测：

$\hat{y}={\widehat{\mathrm{\β}}}_0+{\widehat{\mathrm{\β}}}_1{x}_1+...+{\widehat{\mathrm{\β}}}_p{x}_p---\left(2\right)$

式(2)为经验回归方程，

称为y的预测值;

随机误差ε，常假定服从正态分布：E(ε)＝0，Var(ε)＝σ²;

A12参数估计

已知n组监测数据(x_i1,x_i2,…,x_ip;y_i)，样本预测的误差：

$y_i-{\hat{y}}_i=y_i-{\mathrm{\β}}_0-{\mathrm{\β}}_1{x}_{i1}-{\mathrm{\β}}_2{x}_{i2}...-{\mathrm{\β}}_p{x}_{\mathrm{ip}}---\left(3\right)$

采用最小二乘法估计多元线性回归模型的参数，即使式(4)取最小值时的解;

$Q=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{(y_i-{\hat{y}}_i)}^2=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{(y_i-{\mathrm{\β}}_0-{\mathrm{\β}}_1{x}_{i1}-{\mathrm{\β}}_2{x}_{i2}-...-{\mathrm{\β}}_p{x}_{\mathrm{ip}})}^2---\left(4\right)$

分别对β₀,β₁,…,β_p求偏导数，并令其等于零，然后联立求解即可求得回归参数的估计值

A13温度预测及残差求取

要预测k时刻的齿轮箱温度值

需要与齿轮箱温度相关的监测量作为自变量：k时刻之前两个时刻的环境温度（T_e(k-1)、T_e(k-2)）、k时刻之前两个时刻的齿轮箱油温度（T_o(k-1)、T_o(k-2)）、k时刻之前两个时刻的齿轮轴承温度（T_b(k-1)、T_b(k-2)）。则建立模型的自变量为：X_(k)＝(T_e(k-1),T_e(k-2),T_o(k-1),T_o(k-2),T_b(k-1),T_b(k-2));

齿轮箱温度T_k，采用齿轮箱油温度和齿轮轴承温度的加权和来表示，即T_k＝0.5T_o(k)+0.5T_b(k);

上述参数估计用到的(x_i1,x_i2,…,x_ip;y_i)已经知道（X_(k)、T_k），由最小二乘法估计出模型的参数那么齿轮箱温度的预测值为

温度的残差为

A2随机状态空间模型的识别

随机子空间的线性状态空间模型描述如下：

$\left\{\begin{array}{c}{X}_{k+1}={\mathrm{AX}}_k+w_k\\ Y_k={\mathrm{CX}}_k+v_k\end{array}\right.---\left(5\right)$

方程(5)是随机状态空间模型，其中，X是状态量，一般不可测量，没有实际的物理意义，只是为了便于描述这个系统动态的数学对象；Y是观测量，在这里就是温度的残差；w是系统噪声，由于建模不精确和一些干扰造成的；v是测量噪声。这些量均为相应维数的列向量；A是系统矩阵；C是输出矩阵；

A21正交投影

定义由温度残差组成的分块Hankel矩阵Y：

其中，矩阵Y_p的维数q×N，矩阵Y_f的维数(p+1)×N，N值一般很大，除以

的意义是对测量值进行标准化处理;

对Y阵进行重新分块，如式(7)所示：

其中，

为Y_f的第一行移到了Y_p的末行之后的矩阵，

为Y_f去掉了第一行后的矩阵;

首先，将Y_f正交投影到Y_p的空间上，由正交投影的定义可以计算出O_i的值;

由于O_i一般是非常庞大的矩阵，在实际计算时，我们为了保证数值的稳定性，通常先对Y阵进行分块和LQ分解，得到一个稀疏的下三角阵，也对其进行分块，如式(9)所示，经过推导，投影矩阵O_i亦可以由式(10)求得；又根据随机子空间识别理论，投影矩阵O_i可分解为可观矩阵Γ_i和Kalman滤波状态序列的乘积。

$O_i=Y_f/Y_p=\left(\begin{array}{c}{L}_{21}\\ L_{31}\end{array}\right)Q_1={\mathrm{\Γ}}_i{\hat{X}}_i---\left(10\right)$

${\mathrm{\Γ}}_i=\left(\begin{array}{c}C\\ \mathrm{CA}\\ {\mathrm{CA}}^2\\ .\\ .\\ .\\ {\mathrm{CA}}^{i-1}\end{array}\right)---\left(11\right)$

同样的，对Y阵进行二次分块和LQ分解，并将

正交投影到

的空间上，得到了O_i-1的表达式(12);其中，Γ_i-1为Γ_i去掉最后一行CA^i-1后的矩阵;

$O_{i-1}=Y_f^{-}/Y_p^{+}\left(\begin{array}{cc}{L}_{31}& L_{32}\end{array}\right)\left(\begin{array}{c}{Q}_1\\ Q_2\end{array}\right)={\mathrm{\Γ}}_{i-1}{\hat{X}}_{i+1}---\left(13\right)$

A22奇异值分解

对式(10)中的L阵进行奇异值分解：

$\left(\begin{array}{c}{L}_{21}\\ L_{31}\end{array}\right)=\left(\begin{array}{cc}{U}_1& U_2\end{array}\right)\left(\begin{array}{cc}{S}_1& 0\\ 0& S_2\end{array}\right)\left(\begin{array}{c}{V_1}^T\\ {V_2}^T\end{array}\right)=U_1{S}_1{V_1}^T---\left(14\right)$

其中U、V为酉矩阵；S₁为对角阵，对角阵元素由大到小排列，系统的阶数为非零奇异值的个数，但在实际计算中由于噪声的影响，比较高次的奇异值不会完全等于0，而是和0比较接近的数，在这种情况下，通常奇异值会有一个跳跃，即σ_i+1<<σ_i，选择这个相对很小的奇异值σ_i+1，令这个奇异值近似等于0，并且从这个奇异值开始之后的值都近似处理为0即(σ_i+1≈σ_i+2≈…≈0);

根据奇异值分解(14)，得到可观矩阵Γ_i和状态量

的表达式：

${\mathrm{\&Gamma;}}_i=U_1{S}_1^{1/2}---\left(15\right)$

${\hat{X}}_i=S_1^{1/2}{V_1}^T{Q}_1---\left(16\right)$

再由式(13)，得到

的表达式，其中

为求Moore-Penrose逆;

A23估计系数矩阵A、C

已经估计出状态量根据随机状态空间模型，得到下面的方程;

$\left(\begin{array}{c}{\hat{X}}_{i+1}\\ Y_{i|i}\end{array}\right)=\left(\begin{array}{c}A\\ C\end{array}\right)\left({\hat{X}}_i\right)+\left(\begin{array}{c}{\mathrm{\&rho;}}_w\\ {\mathrm{\&rho;}}_v\end{array}\right)---\left(18\right)$

式中的Y_i|i就是之前出现过的Y_q+1|q+1，为了表示方便，这里对下标进行了改变;ρ_w、ρ_v是残差向量，与不相关;利用最小二乘思想，即令残差项最小时可以解出这个方程，得到A、C阵的估计值;

A3齿轮箱的故障预测

利用大量历史温度残差数据按照A2中所述的识别方法，得到很多组

阵(k＝1,2,…)，并求取

阵的特征值，在正常情况下，每一次计算所得特征值相差不大，也就是说系统的特征是基本一致的，均为正常状态;在这里定义参考特征值

它是历次

阵特征值的平均;

有了参考特征值，利用它与实时数据所求的特征值进行比较来判断齿轮箱的状态好坏;在这里定义离散度R作为评价指标，即由实时数据所求每个特征值与对应参考特征值之间距离的平均值，公式如下：

$R=\frac{1}{n}\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}\sqrt{{(x_i-x_{\mathrm{oi}})}^2+{(y_i-y_{\mathrm{oi}})}^2}---\left(20\right)$

其中，n为系统阶数；x_i、y_i为由实时温度残差数据求得的特征值的实部和虚部；x_oi、y_oi为参考特征值的实部和虚部;

每次计算都得到一个R值，取一段正常运行时间所求得的离散度值，计算其均值μ和根方差σ，定义预警门槛值R_k为：

R_k＝μ+3σ   (21)

有了门槛值R_k就可以通过监测R的变化情况，来实时评估齿轮箱的状态好坏，如果R越过所设定的门槛值，则说明齿轮箱有发生故障的风险。

现在风电场已经实现了对风机齿轮箱的在线温度监测，一般是当温度值越限以后发出警告通知运行人员，所设上限往往较高，例如当轴承温度高于80℃且持续60s则自动停机，可见这种应对方法具有一定的滞后性。另外当出现温度传感器故障或者风机不能自动切机时，会造成润滑油温度持续的升高，严重者甚至造成润滑油的自燃齿轮箱的烧毁。本发明通过分析温度数据所反映的齿轮箱的内在特征，可以在温度不高的故障初期阶段发出预警信号，以减轻齿轮箱的损伤，避免发生不可逆的故障。

附图说明

图1为齿轮箱正常状态的温度预测残差；

图2为齿轮箱异常状态的温度预测残差；

图3为齿轮箱正常状态的特征值；

图4为齿轮箱异常状态的特征值；

图5为正常状态离散度曲线；

图6为异常状态离散度曲线。

具体实施方式

以下结合具体实施例，对本发明进行详细说明。

实施例1

本发明技术方案的实现分为三个步骤：

一、利用回归分析方法对温度数据预处理；

二、利用子空间方法识别随机状态空间模型的参数；

三、实现对齿轮箱的故障预警。

1温度数据的预处理

子空间方法是一种时域分析方法，适合处理类似于振动信号那样的高频率且在某一数值上下波动的信号。利用回归分析方法对温度数据进行单步预测，得到实际值和预测值之间的差值，称为残差，将残差作为随机状态空间模型的观测量Y。

1.1多元线性回归模型

多元线性回归模型的一般形式如下：

y＝β₀+β₁x₁+…+β_px_p+ε   (1)

式(1)中，β₀,β₁,…,β_p是未知参数，β₀为回归常数，β₁,…,β_p为回归系数；y为因变量；x₁,x₂,…,x_p为自变量，这里是与因变量相关的监测量；ε为随机误差。若已知参数的估计值

则可以实现温度的预测：

$\hat{y}={\widehat{\mathrm{\&beta;}}}_0+{\widehat{\mathrm{\&beta;}}}_1{x}_1+...+{\widehat{\mathrm{\&beta;}}}_p{x}_p---\left(2\right)$

式(2)为经验回归方程，

称为y的预测值。

随机误差ε，常假定服从正态分布：E(ε)＝0，Var(ε)＝σ²;。

1.2参数估计

已知n组监测数据(x_i1,x_i2,…,x_ip;y_i)，样本预测的误差：

$y_i-{\hat{y}}_i=y_i-{\mathrm{\&beta;}}_0-{\mathrm{\&beta;}}_1{x}_{i1}-{\mathrm{\&beta;}}_2{x}_{i2}...-{\mathrm{\&beta;}}_p{x}_{\mathrm{ip}}---\left(3\right)$

采用最小二乘法估计多元线性回归模型的参数，即使式(4)取最小值时的解。

$Q=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{(y_i-{\hat{y}}_i)}^2=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{(y_i-{\mathrm{\&beta;}}_0-{\mathrm{\&beta;}}_1{x}_{i1}-{\mathrm{\&beta;}}_2{x}_{i2}-...-{\mathrm{\&beta;}}_p{x}_{\mathrm{ip}})}^2---\left(4\right)$

分别对β₀,β₁,…,β_p求偏导数，并令其等于零，然后联立求解即可求得回归参数的估计值

1.3温度预测及残差求取

要预测k时刻的齿轮箱温度值

需要与齿轮箱温度相关的监测量作为自变量：k时刻之前两个时刻的环境温度（T_e(k-1)、T_e(k-2)）、k时刻之前两个时刻的齿轮箱油温度（T_o(k-1)、T_o(k-2)）、k时刻之前两个时刻的齿轮轴承温度（T_b(k-1)、T_b(k-2)）。则建立模型的自变量为：X_(k)＝(T_e(k-1),T_e(k-2),T_o(k-1),T_o(k-2),T_b(k-1),T_b(k-2))。

齿轮箱温度T_k，采用齿轮箱油温度和齿轮轴承温度的加权和来表示，即T_k＝0.5T_o(k)+0.5T_b(k)。

上文参数估计用到的(x_i1,x_i2,…,x_ip;y_i)已经知道（X_(k)、T_k），由最小二乘法可以估计出模型的参数

那么齿轮箱温度的预测值为

温度的残差为 $e_k=T_k-{\hat{T}}_k.$

2随机状态空间模型的识别

随机子空间的线性状态空间模型描述如下：

$\left\{\begin{array}{c}{X}_{k+1}={\mathrm{AX}}_k+w_k\\ Y_k={\mathrm{CX}}_k+v_k\end{array}\right.---\left(5\right)$

方程(5)是随机状态空间模型，其中，X是状态量，一般不可测量，没有实际的物理意义，只是为了便于描述这个系统动态的数学对象；Y是观测量，在这里就是温度的残差；w是系统噪声，由于建模不精确和一些干扰造成的；v是测量噪声。这些量均为相应维数的列向量。A是系统矩阵；C是输出矩阵。

2.1正交投影

定义由温度残差组成的分块Hankel矩阵Y：

其中，矩阵Y_p的维数q×N，矩阵Y_f的维数(p+1)×N，N值一般很大，除以的意义是对测量值进行标准化处理。

对Y阵进行重新分块，如式(7)所示：

其中，

为Y_f的第一行移到了Y_p的末行之后的矩阵，为Y_f去掉了第一行后的矩阵。

首先，将Y_f正交投影到Y_p的空间上，由正交投影的定义可以计算出O_i的值。

由于O_i一般是非常庞大的矩阵，在实际计算时，我们为了保证数值的稳定性，通常先对Y阵进行分块和LQ分解，得到一个稀疏的下三角阵，也对其进行分块，如式(9)所示，经过推导，投影矩阵O_i亦可以由式(10)求得；又根据随机子空间识别理论，投影矩阵O_i可分解为可观矩阵Γ_i和Kalman滤波状态序列的乘积。

$O_i=Y_f/Y_p=\left(\begin{array}{c}{L}_{21}\\ L_{31}\end{array}\right)Q_1={\mathrm{\&Gamma;}}_i{\hat{X}}_i---\left(10\right)$

${\mathrm{\&Gamma;}}_i=\left(\begin{array}{c}C\\ \mathrm{CA}\\ {\mathrm{CA}}^2\\ .\\ .\\ .\\ {\mathrm{CA}}^{i-1}\end{array}\right)---\left(11\right)$

同样的，对Y阵进行二次分块和LQ分解，并将

正交投影到

的空间上，得到了O_i-1的表达式(12)。其中，Γ_i-1为Γ_i去掉最后一行CA^i-1后的矩阵。

$O_{i-1}=Y_f^{-}/Y_p^{+}\left(\begin{array}{cc}{L}_{31}& L_{32}\end{array}\right)\left(\begin{array}{c}{Q}_1\\ Q_2\end{array}\right)={\mathrm{\&Gamma;}}_{i-1}{\hat{X}}_{i+1}---\left(13\right)$

2.2奇异值分解

对式(10)中的L阵进行奇异值分解：

$\left(\begin{array}{c}{L}_{21}\\ L_{31}\end{array}\right)=\left(\begin{array}{cc}{U}_1& U_2\end{array}\right)\left(\begin{array}{cc}{S}_1& 0\\ 0& S_2\end{array}\right)\left(\begin{array}{c}{V_1}^T\\ {V_2}^T\end{array}\right)=U_1{S}_1{V_1}^T---\left(14\right)$

其中U、V为酉矩阵；S₁为对角阵，对角阵元素由大到小排列，系统的阶数为非零奇异值的个数，但在实际计算中由于噪声的影响，比较高次的奇异值可能不会完全等于0，而是和0比较接近的数，在这种情况下，通常奇异值会有一个跳跃，即σ_i+1<<σ_i，那么我们选择这个相对很小的奇异值σ_i+1，令这个奇异值近似等于0，并且从这个奇异值开始之后的值都近似处理为0即(σ_i+1≈σ_i+2≈…≈0)，这样近似以后，不仅确定了系统的阶数n，而且也降低了矩阵的维数，减少了计算量。

根据奇异值分解(14)，我们就可以得到可观矩阵Γ_i和状态量

的表达式：

Γ_i＝U₁S₁ ^1/2   (15)

${\hat{X}}_i=S_1^{1/2}{V_1}^T{Q}_1---\left(16\right)$

再由式(13)，我们可以得到

的表达式，其中

为求Moore-Penrose逆。

2.3估计系数矩阵A、C

我们已经估计出状态量

根据随机状态空间模型，可以得到下面的方程。

$\left(\begin{array}{c}{\hat{X}}_{i+1}\\ Y_{i|i}\end{array}\right)=\left(\begin{array}{c}A\\ C\end{array}\right)\left({\hat{X}}_i\right)+\left(\begin{array}{c}{\mathrm{\&rho;}}_w\\ {\mathrm{\&rho;}}_v\end{array}\right)---\left(18\right)$

式中的Y_i|i就是之前出现过的Y_q+1|q+1，为了表示方便，这里对下标进行了改变。ρ_w、ρ_v是残差向量，与

不相关。利用最小二乘思想，即令残差项最小时可以解出这个方程，得到A、C阵的估计值。

3齿轮箱的故障预测

利用大量历史温度残差数据按照2中所述的识别方法，可以得到很多组

阵(k＝1,2,…)，并求取阵的特征值，在正常情况下，每一次计算所得特征值相差几部不大，也就是说系统的特征是基本一致的，均为正常状态。在这里定义参考特征值

它是历次

阵特征值的平均。

有了参考特征值，我们就可以利用它与实时数据所求的特征值进行比较来判断齿轮箱的状态好坏。在这里定义离散度R作为评价指标，即由实时数据所求每个特征值与对应参考特征值之间距离的平均值，公式如下：

$R=\frac{1}{n}\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}\sqrt{{(x_i-x_{\mathrm{oi}})}^2+{(y_i-y_{\mathrm{oi}})}^2}---\left(20\right)$

其中，n为系统阶数；x_i、y_i为由实时温度残差数据求得的特征值的实部和虚部；x_oi、y_oi为参考特征值的实部和虚部。

每次计算都可以得到一个R值，我们取一段正常运行时间所求得的离散度值，计算其均值μ和根方差σ，定义预警门槛值R_k为：

R_k＝μ+3σ   (21)

有了门槛值R_k就可以通过监测R的变化情况，来实时评估齿轮箱的状态好坏，如果R越过所设定的门槛值，则说明齿轮箱有发生故障的风险。

实施例2

一、利用回归分析方法对温度数据预处理；

1.求取回归参数的估计值

利用一段时间的环境温度T_e，齿轮箱油温度T_o，齿轮轴承温度T_b估计回归模型参数。假定k时刻的温度值T_k与前2个时刻的环境温度值（T_e(k-1)、T_e(k-2)）、k时刻之前两个时刻的齿轮箱油温度（T_o(k-1)、T_o(k-2)）有关、k时刻之前两个时刻的齿轮轴承温度（T_b(k-1)、T_b(k-2)）有关。T_k是齿轮箱油温度和齿轮轴承温度的平均值。

若已经获得n组这样的监测数据(T_ie1,T_ie2,T_io3,T_io4,T_ib5,T_ib6;T_i)，i＝1,2,…,n，则回归模型为

$\left\{\begin{array}{c}{T}_1={\mathrm{\&beta;}}_0+{\mathrm{\&beta;}}_1{T}_{1e1}+{\mathrm{\&beta;}}_2{T}_{1e2}+{\mathrm{\&beta;}}_3{T}_{1o3}+{\mathrm{\&beta;}}_4{T}_{1o4}+{\mathrm{\&beta;}}_5{T}_{1b5}+{\mathrm{\&beta;}}_6{T}_{1b6}+{\mathrm{\&epsiv;}}_1\\ T_2={\mathrm{\&beta;}}_0+{\mathrm{\&beta;}}_1{T}_{2e1}+{\mathrm{\&beta;}}_2{T}_{2e2}+{\mathrm{\&beta;}}_3{T}_{2o3}+{\mathrm{\&beta;}}_4{T}_{2o4}+{\mathrm{\&beta;}}_5{T}_{2b5}+{\mathrm{\&beta;}}_6{T}_{2b6}+{\mathrm{\&epsiv;}}_2\\ .\\ .\\ .\\ T_n={\mathrm{\&beta;}}_0+{\mathrm{\&beta;}}_1{T}_{\mathrm{ne}1}+{\mathrm{\&beta;}}_2{T}_{\mathrm{ne}2}+{\mathrm{\&beta;}}_3{T}_{\mathrm{no}3}+{\mathrm{\&beta;}}_4{T}_{\mathrm{no}4}+{\mathrm{\&beta;}}_5{T}_{\mathrm{nb}5}+{\mathrm{\&beta;}}_6{T}_{\mathrm{nb}6}+{\mathrm{\&epsiv;}}_n\end{array}\right.$

令

则上式简写为：y＝Xβ+ε

令 $Q=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{(T_i-{\hat{T}}_i)}^2=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{(T_i-{\mathrm{\&beta;}}_0-{\mathrm{\&beta;}}_1{T}_{\mathrm{ie}1}-{\mathrm{\&beta;}}_2{T}_{\mathrm{ie}2}-{\mathrm{\&beta;}}_3{T}_{\mathrm{ib}3}-{\mathrm{\&beta;}}_4{T}_{\mathrm{ib}4}-{\mathrm{\&beta;}}_5{T}_{\mathrm{ib}5}-{\mathrm{\&beta;}}_6{T}_{\mathrm{ib}6})}^2$ 取最小值时的回归系数值即为所求估计参数：

$({\widehat{\mathrm{\&beta;}}}_0,{\widehat{\mathrm{\&beta;}}}_1,{\widehat{\mathrm{\&beta;}}}_2,{\widehat{\mathrm{\&beta;}}}_3,{\widehat{\mathrm{\&beta;}}}_4,{\widehat{\mathrm{\&beta;}}}_5,{\widehat{\mathrm{\&beta;}}}_6)=(-\mathrm{0.0338,0.0051},-\mathrm{0.0073,0.4356,0.0646,0.4102,0.0892})$

2.求温度残差

由于已经得到回归系数β，则k时刻的温度预测值为

温度残差为

即为表1中所列的值。

二、利用子空间方法识别随机状态空间模型的参数；

取Y阵的维数是12×4000，Y_p、Y_f维数都是6×4000。

经过两次分块后分别进行LQ分解

投影矩阵为：

$O_i=Y_f/Y_p=\left(\begin{array}{c}{L}_{21}\\ L_{31}\end{array}\right)Q_1={\mathrm{\&Gamma;}}_i{\hat{X}}_i---\left(5\right)$

$O_{i-1}=Y_f^{-}/Y_p^{+}\left(\begin{array}{cc}{L}_{31}& L_{32}\end{array}\right)\left(\begin{array}{c}{Q}_1\\ Q_2\end{array}\right)={\mathrm{\&Gamma;}}_{i-1}{\hat{X}}_{i+1}---\left(6\right)$

对式(5)中的L进行奇异值分解，得下式(7)，同时确定了系统的阶数为5阶。

$\left(\begin{array}{c}{L}_{21}\\ L_{31}\end{array}\right)=\left(\begin{array}{cc}{U}_1& U_2\end{array}\right)\left(\begin{array}{cc}{S}_1& 0\\ 0& S_2\end{array}\right)\left(\begin{array}{c}{V_1}^T\\ {V_2}^T\end{array}\right)=U_1{S}_1{V_1}^T---\left(7\right)$

可得到可观矩阵和状态序列的估计值为

Γ_i＝U₁S₁ ^1/2   (8)

${\hat{X}}_i=S_1^{1/2}{V_1}^T{Q}_1---\left(9\right)$

利用最小二乘法得到系统矩阵A和输出矩阵C的估计值：

由大量历史数据可以求得一组与该齿轮箱系统特征对应的参考特征值：

参考特征值 ${\mathrm{\&lambda;}}_o=\left(\begin{array}{c}0.6893\\ -0.7015\\ 0.9878\\ -0.2709-j0.9007\\ -0.2709+j0.9007\end{array}\right)$

同时我们利用实时监测的温度残差数据，根据上述过程每次计算都可以得到一组相应的新特征值，将其与参考特征值进行比较，即可以得到系统状态是否正常。表2，表3列出了正常和异常状态时所求特征值。

三、实现对齿轮箱的故障预警。

已知参考特征值λ_o和每次实时计算出的特征值λ_i，根据定义的离散度公式：

$R=\frac{1}{n}\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}\sqrt{{(x_i-x_{\mathrm{oi}})}^2+{(y_i-y_{\mathrm{oi}})}^2}---\left(12\right)$

即可计算出表4中的值。门槛值R_k＝0.0191

由图5可知，在正常运行状态时，所计算出的离散度值，均在所定义的门槛值R_k以下。对于图6，离散度曲线在第15个点时，达到了门槛值附近，经过一段时间的振荡后，完全的越过了门槛值。

表1.正常和异常状态温度残差对比

以上分别是正常状态和异常状态的温度残差值，由表1可以看出，异常状态的温度残差值只是略微比正常状态的残差值要大，从图1、图2的温度残差图也可以看出，在最后的一段时间内，异常状态的齿轮箱温度残差出现了不同于正常状态时的增大。

表2.正常状态的特征值

序号	λ1	λ2	λ3	λ4、λ5实部	λ4、λ5虚部
						1	0.71425	-0.68690	0.98726	-0.27041	0.89830
2	0.71357	-0.70127	0.98817	-0.27932	0.90976
						3	0.71282	-0.69986	0.98843	-0.28053	0.91148
4	0.71483	-0.70256	0.98772	-0.28011	0.91042
						5	0.71972	-0.70323	0.98703	-0.27917	0.91114

6	0.72390	-0.69344	0.98602	-0.27183	0.90346
						7	0.71611	-0.69281	0.98755	-0.26714	0.90443
8	0.71606	-0.69424	0.98712	-0.26432	0.90432
						9	0.71578	-0.69042	0.98598	-0.26883	0.90151
10	0.71646	-0.68964	0.98586	-0.26609	0.89852
						11	0.71648	-0.68948	0.98613	-0.26683	0.89833
12	0.71630	-0.68986	0.98617	-0.26674	0.89842
						13	0.70725	-0.68290	0.98776	-0.26978	0.89804
14	0.70866	-0.68084	0.98757	-0.26866	0.89552
						15	0.70881	-0.68256	0.98762	-0.27028	0.89544
16	0.70195	-0.67893	0.98883	-0.26868	0.89948
						17	0.70193	-0.67777	0.98793	-0.26835	0.89861
18	0.70185	-0.67795	0.98924	-0.26845	0.89893
						19	0.70170	-0.67818	0.98920	-0.26842	0.89888
20	0.72176	-0.69848	0.99233	-0.27601	0.90859
						21	0.72209	-0.69068	0.99557	-0.27301	0.90931
22	0.71900	-0.69643	0.99779	-0.27943	0.90706
						23	0.71407	-0.69255	0.99490	-0.27943	0.90952
24	0.72115	-0.69588	0.99593	-0.27783	0.91059
						25	0.72117	-0.69663	0.99396	-0.27528	0.90831
26	0.72261	-0.69437	0.99174	-0.27080	0.90836
						27	0.71902	-0.69625	0.99600	-0.27343	0.90984
28	0.71696	-0.69493	0.99538	-0.27741	0.90747
						29	0.72255	-0.69658	0.99681	-0.27762	0.90648
30	0.72021	-0.69032	0.98875	-0.27203	0.90658

31	0.71639	-0.69323	0.99125	-0.27141	0.91085
						32	0.71512	-0.69275	0.99468	-0.27414	0.91032
33	0.71873	-0.69325	0.98891	-0.27798	0.91136
						34	0.71631	-0.69534	0.99084	-0.27964	0.90844
35	0.71997	-0.69130	0.99125	-0.27215	0.91093
						36	0.72013	-0.69848	0.99680	-0.27468	0.90944
37	0.71421	-0.69398	0.99210	-0.27105	0.90706
						38	0.72089	-0.69482	0.99333	-0.27992	0.90874
39	0.72159	-0.69633	0.99293	-0.27468	0.90964
						40	0.72246	-0.69944	0.99816	-0.27768	0.91044
41	0.71471	-0.69319	0.99430	-0.27225	0.90928
						42	0.71957	-0.69625	0.99464	-0.27862	0.90670

注：λ₄、λ₅分别为实部±j虚部

表3.异常状态的特征值

序号	λ1	λ2	λ3	λ4、λ5实部	λ4、λ5虚部
						1	0.71512	-0.65650	0.98403	-0.16001	0.85705
2	0.71268	-0.65630	0.98460	-0.15839	0.85880
						3	0.70720	-0.65460	0.98464	-0.16710	0.85752
4	0.70113	-0.65035	0.98576	-0.16386	0.85591
						5	0.70306	-0.65760	0.98593	-0.16868	0.85769
6	0.69748	-0.65533	0.98703	-0.16834	0.86017
						7	0.69537	-0.65331	0.98741	-0.16745	0.85930
8	0.69714	-0.65526	0.98756	-0.16908	0.86038
						9	0.69945	-0.65324	0.98708	-0.16933	0.86166
10	0.71150	-0.65799	0.98306	-0.16943	0.85269

11	0.69513	-0.65147	0.98970	-0.16045	0.85059
						12	0.69506	-0.65362	0.98919	-0.16237	0.85220
13	0.70995	-0.65750	0.98725	-0.16604	0.84567
						14	0.71310	-0.66725	0.98615	-0.17178	0.84481
15	0.72321	-0.67449	0.98167	-0.17169	0.84680
						16	0.72721	-0.67689	0.98303	-0.17621	0.83006
17	0.71882	-0.66274	0.98603	-0.18066	0.83475
						18	0.70669	-0.65446	0.98772	-0.17632	0.83699
19	0.71065	-0.64335	0.98821	-0.16813	0.82779
						20	0.67717	-0.61013	0.98851	-0.13306	0.83538
21	0.65615	-0.60099	0.99112	-0.12782	0.84389
						22	0.65870	-0.60422	0.99071	-0.13099	0.84021
23	0.67884	-0.60809	0.98654	-0.12418	0.82464
						24	0.68548	-0.60706	0.98705	-0.10830	0.81907
25	0.68284	-0.60721	0.98721	-0.10716	0.82143
						26	0.69067	-0.62130	0.98787	-0.11433	0.82443
27	0.68998	-0.62146	0.98669	-0.11294	0.82777
						28	0.69420	-0.62459	0.98161	-0.10467	0.82495
29	0.69188	-0.63276	0.98461	-0.10670	0.82967
						30	0.70025	-0.64094	0.98272	-0.10605	0.83253
31	0.69807	-0.65680	0.98622	-0.11306	0.85032
						32	0.71421	-0.66557	0.98134	-0.12296	0.84578
33	0.70176	-0.66436	0.98315	-0.12255	0.85036
						34	0.70460	-0.66562	0.98278	-0.12323	0.84932
35	0.70628	-0.66469	0.98147	-0.12206	0.84778

36	0.68948	-0.66988	0.98622	-0.11940	0.84722
						37	0.68956	-0.67457	0.98627	-0.12059	0.84508
38	0.68168	-0.66290	0.98575	-0.12861	0.83344
						39	0.68385	-0.65637	0.98499	-0.13147	0.82407
40	0.68469	-0.64743	0.98663	-0.13902	0.81266
						41	0.68247	-0.64819	0.98686	-0.13928	0.81610
42	0.68211	-0.64233	0.98766	-0.13673	0.81829

由表2可以看出，正常状态时系统矩阵的特征值，基本都在参考特征值附近变化，而且变化较小；而由表3，当齿轮箱异常时，由随机子空间算法识别的系统矩阵特征值变化较剧烈，且呈现偏离参考特征值的趋势，可见齿轮箱的状态出现了劣化的趋势。

表4.正常和异常状态离散度对比

由表4可以看出，当齿轮箱正常运行时，所得有关特征值的离散度值较小，且变化波动较小，而当齿轮箱处于异常状态时，一开始其离散度值也变化较小，但是随着图6中对应的曲线，离散度值变化越来越大，并且越过了所设定的门槛值，并且越过门槛值以后的离散度值的波动程度，明显较正常状态时的更剧烈。

应当理解的是，对本领域普通技术人员来说，可以根据上述说明加以改进或变换，而所有这些改进和变换都应属于本发明所附权利要求的保护范围。

Claims (1)

1.一种基于温度数据的风机齿轮箱子空间故障预测方法，其特征在于，包括以下步骤：

A1温度数据的预处理；利用回归分析方法对温度数据进行单步预测，得到实际值和预测值之间的差值，称为残差，将残差作为随机状态空间模型的观测量Y;

A11多元线性回归模型

多元线性回归模型的一般形式如下：

y＝β₀+β₁x₁+…+β_px_p+ε   (1)

式(1)中，β₀,β₁,…,β_p是未知参数，β₀为回归常数，β₁,…,β_p为回归系数；y为因变量；x₁,x₂,…,x_p为自变量，这里是与因变量相关的监测量；ε为随机误差。若已知参数的估计值

则可以实现温度的预测：

$\hat{y}={\widehat{\mathrm{\&beta;}}}_0+{\widehat{\mathrm{\&beta;}}}_1{x}_1+...+{\widehat{\mathrm{\&beta;}}}_p{x}_p---\left(2\right)$

式(2)为经验回归方程，

称为y的预测值;

随机误差ε，常假定服从正态分布：E(ε)＝0，Var(ε)＝σ²;

A12参数估计

已知n组监测数据(x_i1,x_i2,…,x_ip;y_i)，样本预测的误差：

$y_i-{\hat{y}}_i=y_i-{\mathrm{\&beta;}}_0-{\mathrm{\&beta;}}_1{x}_{i1}-{\mathrm{\&beta;}}_2{x}_{i2}...-{\mathrm{\&beta;}}_p{x}_{\mathrm{ip}}---\left(3\right)$

采用最小二乘法估计多元线性回归模型的参数，即使式(4)取最小值时的解;

$Q=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{(y_i-{\hat{y}}_i)}^2=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{(y_i-{\mathrm{\&beta;}}_0-{\mathrm{\&beta;}}_1{x}_{i1}-{\mathrm{\&beta;}}_2{x}_{i2}-...-{\mathrm{\&beta;}}_p{x}_{\mathrm{ip}})}^2---\left(4\right)$

分别对β₀,β₁,…,β_p求偏导数，并令其等于零，然后联立求解即可求得回归参数的估计值

A13温度预测及残差求取

要预测k时刻的齿轮箱温度值

需要与齿轮箱温度相关的监测量作为自变量：k时刻之前两个时刻的环境温度（T_e(k-1)、T_e(k-2)）、k时刻之前两个时刻的齿轮箱油温度（T_o(k-1)、T_o(k-2)）、k时刻之前两个时刻的齿轮轴承温度（T_b(k-1)、T_b(k-2)）。则建立模型的自变量为：X_(k)＝(T_e(k-1),T_e(k-2),T_o(k-1),T_o(k-2),T_b(k-1),T_b(k-2));

齿轮箱温度T_k，采用齿轮箱油温度和齿轮轴承温度的加权和来表示，即T_k＝0.5T_o(k)+0.5T_b(k);

上述参数估计用到的(x_i1,x_i2,…,x_ip;y_i)已经知道（X_(k)、T_k），由最小二乘法估计出模型的参数

那么齿轮箱温度的预测值为

温度的残差为 $e_k=T_k-{\hat{T}}_k;$

A2随机状态空间模型的识别

随机子空间的线性状态空间模型描述如下：

$\left\{\begin{array}{c}{X}_{k+1}={\mathrm{AX}}_k+w_k\\ Y_k={\mathrm{CX}}_k+v_k\end{array}\right.---\left(5\right)$

方程(5)是随机状态空间模型，其中，X是状态量，一般不可测量，只是为了便于描述这个系统动态的数学对象；Y是观测量，在这里就是温度的残差；w是系统噪声，由于建模不精确和一些干扰造成的；v是测量噪声。这些量均为相应维数的列向量；A是系统矩阵；C是输出矩阵；

A21正交投影

定义由温度残差组成的分块Hankel矩阵Y：

其中，矩阵Y_p的维数q×N，矩阵Y_f的维数(p+1)×N，N值一般很大，除以

的意义是对测量值进行标准化处理;

对Y阵进行重新分块，如式(7)所示：

其中，为Y_f的第一行移到了Y_p的末行之后的矩阵，

为Y_f去掉了第一行后的矩阵;

首先，将Y_f正交投影到Y_p的空间上，由正交投影的定义可以计算出O_i的值;

由于O_i一般是非常庞大的矩阵，在实际计算时，我们为了保证数值的稳定性，通常先对Y阵进行分块和LQ分解，得到一个稀疏的下三角阵，也对其进行分块，如式(9)所示，经过推导，投影矩阵O_i亦可以由式(10)求得；又根据随机子空间识别理论，投影矩阵O_i可分解为可观矩阵Γ_i和Kalman滤波状态序列的乘积。

$O_i=Y_f/Y_p=\left(\begin{array}{c}{L}_{21}\\ L_{31}\end{array}\right)Q_1={\mathrm{\&Gamma;}}_i{\hat{X}}_i---\left(10\right)$

${\mathrm{\&Gamma;}}_i=\left(\begin{array}{c}C\\ \mathrm{CA}\\ {\mathrm{CA}}^2\\ .\\ .\\ .\\ {\mathrm{CA}}^{i-1}\end{array}\right)---\left(11\right)$

同样的，对Y阵进行二次分块和LQ分解，并将正交投影到

的空间上，得到了O_i-1的表达式(12);其中，Γ_i-1为Γ_i去掉最后一行CA^i-1后的矩阵;

$O_{i-1}=Y_f^{-}/Y_p^{+}\left(\begin{array}{cc}{L}_{31}& L_{32}\end{array}\right)\left(\begin{array}{c}{Q}_1\\ Q_2\end{array}\right)={\mathrm{\&Gamma;}}_{i-1}{\hat{X}}_{i+1}---\left(13\right)$

A22奇异值分解

对式(10)中的L阵进行奇异值分解：

$\left(\begin{array}{c}{L}_{21}\\ L_{31}\end{array}\right)=\left(\begin{array}{cc}{U}_1& U_2\end{array}\right)\left(\begin{array}{cc}{S}_1& 0\\ 0& S_2\end{array}\right)\left(\begin{array}{c}{V_1}^T\\ {V_2}^T\end{array}\right)=U_1{S}_1{V_1}^T---\left(14\right)$

其中U、V为酉矩阵；S₁为对角阵，对角阵元素由大到小排列，系统的阶数为非零奇异值的个数，但在实际计算中由于噪声的影响，比较高次的奇异值不会完全等于0，而是和0比较接近的数，在这种情况下，通常奇异值会有一个跳跃，即σ_i+1<<σ_i，选择这个相对很小的奇异值σ_i+1，令这个奇异值近似等于0，并且从这个奇异值开始之后的值都近似处理为0即（σ_i+1≈σ_i+2≈…≈0）;

根据奇异值分解(14)，得到可观矩阵Γ_i和状态量

的表达式：

${\mathrm{\&Gamma;}}_i=U_1{S}_1^{1/2}---\left(15\right)$

${\hat{X}}_i=S_1^{1/2}{V_1}^T{Q}_1---\left(16\right)$

再由式(13)，得到的表达式，其中

为求Moore-Penrose逆;

A23估计系数矩阵A、C

已经估计出状态量

根据随机状态空间模型，得到下面的方程;

$\left(\begin{array}{c}{\hat{X}}_{i+1}\\ Y_{i|i}\end{array}\right)=\left(\begin{array}{c}A\\ C\end{array}\right)\left({\hat{X}}_i\right)+\left(\begin{array}{c}{\mathrm{\&rho;}}_w\\ {\mathrm{\&rho;}}_v\end{array}\right)---\left(18\right)$

式中的Y_i|i就是之前出现过的Y_q+1|q+1，为了表示方便，这里对下标进行了改变;ρ_w、ρ_v是残差向量，与

不相关;利用最小二乘思想，即令残差项最小时可以解出这个方程，得到A、C阵的估计值;

A3齿轮箱的故障预测

利用大量历史温度残差数据按照A2中所述的识别方法，得到很多组

阵（k＝1,2…），并求取

阵的特征值，在正常情况下，每一次计算所得特征值相差不大，也就是说系统的特征是基本一致的，均为正常状态;在这里定义参考特征值

它是历次

阵特征值的平均;

有了参考特征值，利用它与实时数据所求的特征值进行比较来判断齿轮箱的状态好坏;在这里定义离散度R作为评价指标，即由实时数据所求每个特征值与对应参考特征值之间距离的平均值，公式如下：

$R=\frac{1}{n}\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}\sqrt{{(x_i-x_{\mathrm{oi}})}^2+{(y_i-y_{\mathrm{oi}})}^2}---\left(20\right)$

其中，n为系统阶数；x_i、y_i为由实时温度残差数据求得的特征值的实部和虚部；x_oi、y_oi为参考特征值的实部和虚部;

每次计算都得到一个R值，取一段正常运行时间所求得的离散度值，计算其均值μ和根方差σ，定义预警门槛值R_k为：

R_k＝μ+3σ   (21)

有了门槛值R_k就可以通过监测R的变化情况，来实时评估齿轮箱的状态好坏，如果R越过所设定的门槛值，则说明齿轮箱有发生故障的风险。

Images