CN104537251B - 一种风机叶片的冲击载荷识别方法 - Google Patents

Abstract

本发明涉及一种风机叶片的冲击载荷识别方法，该方法基于共轭梯度最小二乘迭代算法，包括以下步骤：1)采用锤击法测量冲击载荷施加的风机叶片位置点与风机叶片的响应测点间的频响函数，通过快速傅里叶逆变换得到单位脉冲响应函数，进而解卷积获得传递矩阵；2)采用加速度传感器测量由作用于风机叶片的冲击载荷产生的加速度响应信号；3)初始化共轭梯度最小二乘迭代算法的迭代步数、待识别载荷、响应残差和迭代搜索方向；4)更新当前迭代步数下的迭代步数、待识别载荷、响应残差、共轭梯度系数和迭代搜索方向；5)判断是否满足启发式收敛终止准则；6)获得正则化的待识别载荷。

Description

一种风机叶片的冲击载荷识别方法

技术领域

本发明属于机械系统状态监测领域，具体涉及一种风机叶片的冲击载荷识别方法，该方法基于共轭梯度最小二乘迭代算法。

背景技术

现如今复合材料已经被广泛应用在风力发电的风机叶片上，由于风机运行以及维护过程中，不可避免的遭受到风沙、飞鸟、冰雹、维修工具等外来物的冲击，并且发生频率较高，冲击损伤积累下来会对风机叶片的完整性以及承载能力造成安全隐患。因此，识别出风机叶片的冲击载荷，对风机的正常运行和维护具有重要意义，同时对于叶片的安全性能的评估也具有指导意义。风机叶片的健康监测分为主动监测和被动监测，主动监测技术是利用激振器在结构中激发出信号，由传感器采集结构的响应信号，进而确定风机叶片的状态。被动监测系统中不存在激振器，仅采用传感器对结构的变化进行检测，再根据传感器采集的信息结合建模技术、信号处理方法确定风机叶片状态。冲击载荷识别技术属于新兴的被动监测技术。然而，由于测量点的不可达难以对作用于结构的外载荷特别是冲击载荷直接测量或计算。

载荷识别技术分为频域法和时域法。时域载荷识别方法不同于传统的频域法，无须在各个频点对频响函数矩阵求逆，而是根据载荷与系统传递函数之间的卷积关系，通过解卷积获取载荷的时域离散信号，能够处理瞬态信号，因而在冲击载荷识别方面应用较多。载荷识别作为一个典型的反问题，不满足Hadamard三个适定性条件(连续、稳定、唯一)，而需要研究正则化方法克服载荷识别的不适定性。

发明内容

本发明目的是提供一种识别精度优于经典的Landweber迭代正则化算法的基于共轭梯度最小二乘迭代正则化算法的载荷识别方法，用于解决当前风机叶片冲击载荷识别中存在的不适定性问题。

为达到上述目的，本发明采用的技术方案是，一种风机叶片的冲击载荷识别方法，该方法基于共轭梯度最小二乘迭代算法，其特征在于，将求解反问题的共轭梯度最小二乘迭代正则化算法，应用在复合材料风机叶片冲击载荷识别领域，并给出了迭代收敛的启发式终止准则。该方法具体包括如下步骤：

1)测量频响函数和计算传递矩阵。采用锤击法(一种试验模态测试方法)测量冲击载荷施加的风机叶片位置点与风机叶片的加速度响应测点间的频响函数H(ω)，通过快速傅里叶逆变换(InverseFastFourierTransform，IFFT)得到单位脉冲响应函数h(t)，进而解卷积获得传递矩阵H。其中，ω表示圆频率变量，t表示时间变量；

2)施加冲击载荷和采集加速度响应。采用加速度传感器测量由作用于风机叶片的冲击载荷产生的加速度响应信号y；

3)初始化迭代步数m、待识别载荷f₀、响应残差r₀和迭代搜索方向d₀：

m＝0，f₀＝0，r₀＝y-Hf₀，d₀＝H^Tr₀(1)

4)更新当前迭代步数m下的迭代步长α_m、待识别载荷f_m、响应残差r_m、共轭梯度系数β_m和迭代搜索方向d_m：

${\mathrm{\α}}_m=\frac{{\left|\right|H^T{r}_{m-1}\left|\right|}_2^2}{{\left|\right|{\mathrm{Hd}}_{m-1}\left|\right|}_2^2}---\left(2\right)$

f_m＝f_m-1+α_md_m-1(3)

r_m＝r_m-1-a_mHd_m(4)

${\mathrm{\β}}_m=\frac{{\left|\right|H^T{r}_m\left|\right|}_2^2}{{\left|\right|{\mathrm{Hr}}_m\left|\right|}_2^2}---\left(5\right)$

d_m＝H^Tr_m+β_md_m-1(6)

其中，||·||₂表示向量的二范数；

5)判断当前迭代步数m下的待识别载荷f_m是否满足启发式收敛终止准则η_m：

η_m＝|p′_m(0)|^1/2||y-Hf_m||₂(7)

其中，p_m为Ritz残差多项式，在共轭梯度最小二乘正则化算法过程中计算而得：

$\begin{array}{cc}\left|p_0^{\′}\left(0\right)\right|=0,& \left|p_1^{\′}\left(0\right)\right|={\mathrm{\α}}_0,\\ \left|p_{m+1}^{\′}\left(0\right)\right|=\left|p_m^{\′}\left(0\right)\right|+{\mathrm{\α}}_m+{\mathrm{\α}}_m\frac{{\mathrm{\β}}_m}{{\mathrm{\α}}_{m-1}}\left(\right|p_m^{\′}\left(0\right)|-|p_{m-1}^{\′}\left(0\right)\left|\right),& m\≥1\end{array}---\left(8\right)$

6)若启发式收敛准则有最小值，则终止迭代计算，获得重构冲击载荷的最优迭代正则化解f_m；否则，迭代步数m＝m+1，返回步骤4)继续迭代计算，直到终止迭代计算，获得重构冲击载荷的最优迭代正则化解f_m。

本发明与现有技术相比具有下列优点：

1.该共轭梯度最小二乘(ConjugateGradientLeastSquares，CGLS)迭代正则化算法不同于截断奇异值分解、Tikhonov正则化方法，具有不涉及传递矩阵求逆运算和不需要明确正则化参数的优点，由于迭代过程即是正则化的过程，可广泛应用在反问题求解中；

2.本发明给出了启发式迭代收敛终止准则并与共轭梯度最小二乘迭代正则化算法结合可高精度地重构作用于风机叶片的冲击载荷；

3.与传统的Landweber迭代正则化算法相比，CGLS迭代算法收敛速度快、识别结果精度高、耗时少。

附图说明

图1是本发明方法流程图；

图2是实施例复合材料风机叶片载荷识别装置示意图；

图3是传感器测量的冲击信号，其中，3(a)力信号，3(b)加速度信号；

图4是启发式终止准则在CGLS正则化算法中应用结果，其中，4(a)冲击事件1，4(b)冲击事件2；

图5是载荷识别结果，其中，5(a)冲击事件1，5(b)冲击事件2；

图6是CGLS和Landweber正则化迭代算法随迭代步数的相对误差，其中，6(a)冲击事件1，6(b)冲击事件2。

具体实施方式

下面结合附图1-6及一个具体实施例对本发明作进一步描述，应该强调的是，下述说明仅仅是示例性的，而本发明的应用对象不局限下述示例。

图1是本发明完成的一种风机叶片的冲击载荷识别方法的流程图，，该方法基于共轭梯度最小二乘迭代算法，可知该方法根据测量的频响函数和加速度信号，将共轭梯度最小二乘正则化迭代算法结合启发式收敛准则，克服了载荷识别中存在的病态问题，达到了载荷识别的目的，具体步骤如下：

1)测量频响函数和计算传递矩阵。采用锤击法测量冲击载荷施加的风机叶片位置点与风机叶片的加速度响应测点间的频响函数H(ω)，通过快速傅里叶逆变换(InverseFastFourierTransform，IFFT)得到单位脉冲响应函数h(t)，进而解卷积获得传递矩阵H。其中，ω表示圆频率变量，t表示时间变量；其中，所述锤击法为本领域常用的试验模态测试方法；

11)试验示意图如图2所示，选用的复合材料风机叶片长度为702mm、最大宽度为112mm和厚度为11mm。为了模拟风机叶片真实的安装条件，用悬臂梁支座的固支段加紧风机叶片根部0～60mm部位，型号为PCB333B50加速度传感器布置在距固支端根部120～160mm位置。模拟冲击力作用的位置在距离叶片根部70％处，已有资料表明该部位易受到冲击载荷作用。在无噪声的实验室环境下，采用型号PCB086C02的脉冲力锤(锤头顶部嵌有力传感器)，选择不锈钢锤帽，重复敲击作用点五次，同时由LMSSCADASIII数据采集系统同步记录冲击力和加速度信号，五次冲击载荷作用点到加速度测点间的频响函数为H₁(ω)，H₂(ω)，H₃(ω)，H₄(ω)和H₅(ω)，由LMSIMPACT模块计算得到其平均值为H(ω)；

2)施加冲击载荷和采集加速度响应。采用加速度传感器测量由作用于风机叶片的冲击载荷产生的加速度响应信号y；

21)采用软质锤帽，使用冲击力锤连续敲击风机叶片易受冲击位置两次(如图2(a)所示)，并同时由LMSSCADASIII数据采集以2048Hz的采样频率，同步记录加速度信号和冲击载荷信号；传感器测量的加速度信号和冲击载荷信号如图3所示，来自冲击力锤的力信号作为载荷识别的比较对象；注意该步骤施加冲击力的作用点与步骤1)测量频响函数的作用点一致，同时加速度位置始终保持不变。其中，软质锤帽的载荷激励频带小于硬质锤帽的激励频带。

22)冲击试验时，冲击力锤连续敲击风机叶片点两次，分别标记为冲击事件1和冲击事件2，且第二次冲击在前次冲击彻底衰减后再施加；

23)在冲击载荷识别中，为了应用本发明的迭代正则化算法分别识别每次冲击载荷，对每次冲击事件截取的数据长度为2050、冲击持续时间为1s；

3)初始化迭代步数m、待识别载荷f₀、响应残差r₀和迭代搜索方向d₀：

m＝0，f₀＝0，r₀＝y-Hf₀，d₀＝H^Tr₀(1)

31)本实例的迭代总步数设置为1000步，计算环境为Win732、内存4G和CPUi5-3450；

4)更新当前迭代步数m下的迭代步长α_m、待识别载荷f_m、响应残差r_m、共轭梯度系数β_m和迭代搜索方向d_m：

${\mathrm{\α}}_m=\frac{{\left|\right|H^T{r}_{m-1}\left|\right|}_2^2}{{\left|\right|{\mathrm{Hd}}_{m-1}\left|\right|}_2^2}---\left(2\right)$

f_m＝f_m-1+α_md_m-1(3)

r_m＝r_m-1-a_mHd_m-1(4)

${\mathrm{\β}}_m=\frac{{\left|\right|H^T{r}_m\left|\right|}_2^2}{{\left|\right|{\mathrm{Hr}}_m\left|\right|}_2^2}---\left(5\right)$

d_m＝H^Tr_m+β_md_m-1(6)

5)判断当前迭代步数m下的待识别载荷f_m是否满足启发式收敛终止准则η_m：

η_m＝|p′_m(0)|^1/2||y-Hf_m||₂(7)

其中，p_m为Ritz残余多项式。在共轭梯度最小二乘算法过程中可以计算而得：

$\begin{array}{cc}\left|p_0^{\′}\left(0\right)\right|=0,& \left|p_1^{\′}\left(0\right)\right|={\mathrm{\α}}_0,\\ \left|p_{m+1}^{\′}\left(0\right)\right|=\left|p_m^{\′}\left(0\right)\right|+{\mathrm{\α}}_m+{\mathrm{\α}}_m\frac{{\mathrm{\β}}_m}{{\mathrm{\α}}_{m-1}}\left(\right|p_m^{\′}\left(0\right)|-|p_{m-1}^{\′}\left(0\right)\left|\right),& m\≥1\end{array}---\left(8\right)$

6)若启发式收敛准则有最小值，则终止迭代计算，获得重构冲击载荷的最优迭代正则化解f_m；否则，迭代步数m＝m+1，返回步骤4)继续迭代计算，直到终止迭代计算，获得重构冲击载荷的最优迭代正则化解f_m。

61)图4是CGLS迭代算法的启发式终止准则的收敛图。从图4可知，冲击事件1和冲击事件2分别需要255和334个迭代步数获得最优迭代正则化解f_m，记为f₂₅₅和f₃₃₄；

62)图5是两次冲击事件应用CGLS迭代算法的载荷识别结果，从图5可知，利用CGLS迭代算法重构的两次冲击事件的正则化结果f₂₅₅和f₃₃₄与实测的载荷结果一致，特别是在冲击力极值点附近；

7)将基于共轭梯度最小二乘迭代正则化算法(CGLS)的载荷识别结果与实测载荷(Exact)、基于Landweber迭代正则化算法(LW)的载荷识别结果比较，比较结果如图5、图6和表1。

71)为了定量评价两种迭代算法所识别载荷的精度，定义载荷识别精度的相对误差为

$\frac{{\left|\right|f_{\mathrm{exact}}-f_{\mathrm{identified}}\left|\right|}_2}{{{\left|\right|f}_{\mathrm{real}}\left|\right|}_2}\×100\%---\left(6\right)$

其中，f_exact和f_identified分别是力传感器实测的冲击载荷和应用迭代算法重构的冲击载荷。

将实测载荷作为参考，可以计算每个迭代步的相对误差，由图6可知两次载荷识别结果的最小相对误差的迭代步数分别是282和374。从表1可知，两次冲击事件的正则化相对误差分别为36.63％和32.14％，迭代过程中最小相对误差分别为36.60％和31.82％。可知启发式终止准则确定的正则化解，非常接近于迭代过程中最小相对误差对应的最优解，因此将启发式终止准则用来确定CGLS算法的收敛步数是可行的。

表1CGLS和Landweber迭代算法对比结果

进一步，比较CGLS和Landweber两种迭代正则化算法的载荷识别性能。从图6可知，在1000步内，CGLS相对误差快速下降然后缓慢上升，符合“半收敛”迭代算法收敛特点，而Landweber相对误差一直缓慢下降，且在1000步内尚未收敛。实际上，本实施例中的两次冲击载荷识别应用Landweber迭代正则化算法在30000迭代步数以内，尚未获得最小的相对误差。可知Landweber相对CGLS迭代算法收敛速度非常慢，且在相同的计算迭代次数，前者仅需后者一半的计算时间。从图5中可知，Landweber迭代算法相对误差比较大，主要在于不能精确的识别冲击加载区的大幅值冲击力。

通过上述的两次冲击载荷识别实例分析可知：共轭梯度最小二乘迭代正则化算法结合启发式终止准则可高精度重构作用于复合材料风机叶片的冲击载荷；所识别的冲击载荷的结果可以为风机叶片是否收到致命冲击损伤提供依据；与传统的Landweber迭代正则化算法相比，CGLS迭代算法收敛速度快、识别结果精度高、耗时少。

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims (2)

1.一种风机叶片的冲击载荷识别方法，该方法基于共轭梯度最小二乘迭代算法，包括如下步骤：

1)测量频响函数和计算传递矩阵，采用锤击法测量冲击载荷施加的风机叶片位置点与风机叶片的加速度响应测点间的频响函数H(ω)，通过快速傅里叶逆变换得到单位脉冲响应函数h(t)，进而解卷积获得传递矩阵H，其中，ω表示圆频率变量，t表示时间变量；

2)施加冲击载荷和测量加速度响应，采用加速度传感器测量由作用于风机叶片的冲击载荷产生的加速度响应信号y；

3)初始化迭代步数m、待识别载荷f₀、响应残差r₀和迭代搜索方向d₀：

m＝0，f₀＝0，r₀＝y-Hf₀，d₀＝H^Tr₀；

4)更新当前迭代步数m下的迭代步长α_m、待识别载荷f_m、响应残差r_m、共轭梯度系数β_m和迭代搜索方向d_m：

${\mathrm{\α}}_m=\frac{\left|\right|H^T{r}_{m-1}|{|}_2^2}{\left|\right|{\mathrm{Hd}}_{m-1}|{|}_2^2};$

f_m＝f_m-1+α_md_m-1；

r_m＝r_m-1-a_mHd_m-1；

${\mathrm{\β}}_m=\frac{\left|\right|H^T{r}_m|{|}_2^2}{\left|\right|{\mathrm{Hr}}_{m-1}|{|}_2^2};$

d_m＝H^Tr_m+β_md_m-1；

其中，||·||₂表示向量的二范数；

5)判断当前迭代步数m下的待识别载荷f_m是否满足启发式收敛终止准则η_m：

η_m＝|p′_m(0)|^1/2||y-Hf_m||₂；

其中，p_m为Ritz残差多项式，可在上述迭代过程中计算得到：

$\begin{array}{cc}\left|p_0^{\′}\left(0\right)\right|=0,& \left|p_1^{\′}\left(0\right)\right|={\mathrm{\α}}_0,\\ \left|p_{m+1}^{\′}\left(0\right)\right|=\left|p_m^{\′}\left(0\right)\right|+{\mathrm{\α}}_m+{\mathrm{\α}}_m\frac{{\mathrm{\β}}_m}{{\mathrm{\α}}_{m-1}}\left(\right|p_m^{\′}\left(0\right)|-|p_{m-1}^{\′}\left(0\right)\left|\right),& m\≥1\end{array};$

6)若启发式收敛准则有最小值，则终止迭代计算，获得重构冲击载荷的最优迭代正则化解f_m；否则，迭代步数m＝m+1，返回步骤4)继续迭代计算，直到终止迭代计算，获得重构冲击载荷的最优迭代正则化解f_m。

2.根据权利要求1所述的风机叶片的冲击载荷识别方法，其特征在于，所述的步骤6)中，共轭梯度最小二乘迭代算法的迭代步数m起到了正则化参数的作用，迭代过程即是正则化的过程。