基于摄影测量和ERA改进算法的叶片模态辨识方法

全伟铭; 郭迎福; 王文韫; 邹龙洲

doi:doi:10.3788/LOP202158.0411001

激光与光电子学进展, 2021, 58 (4): 0411001, 网络出版: 2021-02-04

基于摄影测量和ERA改进算法的叶片模态辨识方法下载： 875次

Modal Identification Method of Blade Based on Photogrammetry and Improved ERA

论文大纲

全伟铭 ¹郭迎福 ¹王文韫 ^2,*邹龙洲 ¹

作者单位

¹ 湖南科技大学机电工程学院, 湖南湘潭 411201

² 机械设备健康维护湖南省重点实验室, 湖南湘潭 411201

AI 词云图 AI一句话精读 AI短摘要

注：本部分内容由 AI 自动生成，请您知悉。

摘要

在风力机叶片表面排列编码目标,基于双目立体成像原理和图像特征匹配与跟踪方法,可以动态、实时地获得编码目标的空间运动信息。针对噪声污染影响模态辨识的问题,通过奇异值分解(SVD)和Cadzow算法构造了无噪声污染的Hankel矩阵,有效减少由噪声引起的虚假模态。针对模态定阶过程容易出现虚假模态的问题,提出一种改进的特征系统实现算法(ERA)。相比于传统的ERA模态辨识方法,该方法在固有频率的识别上提高了3.6%的精度,获取了更清晰的稳定图。

Abstract

Through arranging coded markers on the blade surface of a wind turbine, and based on binocular stereo imaging principle and image feature matching and tracking method, the spatial motion information of coded markers can be obtained dynamically and in real time. Aiming at the problem of noise pollution affecting modal identification, a singular value decomposition (SVD) and Cadzow algorithm was used to construct a Hankel matrix without noise pollution, which effectively reduces false modals caused by noise. In order to solve the problem that false modes are easy to appear in the process of modal order determination, an improved eigensystem realization algorithm (ERA) is proposed. Compared with the traditional ERA, the proposed method improves the accuracy of natural frequency identification by 3.6% and obtains a clearer stabilization diagram.

1 引言

风力机叶片是既昂贵又容易损坏的风电组件之一,随着风能产业的预期增长和未来对风能的依赖,迫切需要一种监测叶片运行状态的有效方法^[1]。运行模态分析是一种通过测量结构在运行状态下的动态响应,提取结构模态特征的技术,可作为检测叶片健康状况和定量评估结构运行状况的一种有效方法^[2-3]。结构响应的获取通常都是采用接触式的测量方法,如应变片、差分应变发生器和其他接触式传感器。然而,接触式的测量方法不仅容易受到传感器数量和采集通道的限制,而且测量系统的布置较为困难,使得风力机叶片的实时测量存在许多局限^[4-5]。摄影测量是一种通过分析照片,提取拍摄对象几何形状、位移和变形的学科。该方法不会对测试对象带来质量载荷和刚度变化,具有非接触、高精度、全场测量的优点,在大尺寸全场测量方面已有许多成熟的应用^[6-9]。本文在风力机叶片表面粘贴多个编码标志,通过双目相机的连续拍摄和摄影测量方法获取编码标志的三维空间运动信息,避免了因为传感器的附加质量和采集数据不足而引起的测量误差。

特征系统实现算法(ERA)是一种利用线性系统中Ho-Kalman的最小实现理论来识别系统模态参数的时域方法,该方法便于确定模态阶次且识别速度较快,同时对具有低频和模态密集特征的结构具有较强的辨识能力^[2]。Li等^[10]开发了一种基于ERA的数据驱动方法,并将其应用在机电振荡动态模式的提取中。Castellanos-Toro等^[11]使用智能手机获取用于模态分析的加速度数据,结合自然激发技术和ERA算法,成功辨识了400多座桥梁的固有频率和阻尼比。但由于脉冲响应中存在噪声或结构非线性较小的问题,模态辨识方法通常会出现虚假模态和错误定阶等问题,作为帮助区分物理极点和虚假极点的稳定图容易产生错误结果。针对数据噪声的问题,基于奇异值分解(SVD)的噪声消除算法相对于其他算法具有快捷且易于实现的优点。Jiang等^[12]将SVD用于重新分配连续小波变换的系数矩阵,以获得振动信号的去噪和更清晰的时频曲线。Golafshan等^[13]提出一种基于SVD和Hankel矩阵的球轴承人工故障去噪方法。本文采用SVD和Cadzow算法^[14]构造了无噪声污染的Hankel矩阵,并将其用于ERA算法,一定程度上减少了噪声对稳定图的影响。

针对模态定阶过程容易出现虚假模态的问题,模态振幅相干系数、模态相位共线性指标和一致模态等指标在ERA稳定图中被引入,用于区分真实模态和虚假模态^[15-17]。Zhang等^[15]使用EEMD方法降低数据噪声,根据奇异熵增量谱确定系统阶数,再通过设定频率和阻尼比的阈值进一步消除虚假模态,实现屋顶溢流水电站动力装置的模态参数识别。Moaveni等^[16]通过对采样频率的合理选择和有效特征值的筛选,改善了ERA算法在系统阶数未知和初始条件不为零时的辨识结果,但是该方法需多次试算来确定适当的采样频率和特征值,计算量较大。Chiang等^[17]通过建立平稳响应的相关函数矩阵,研究与系统响应相关的正弦值的分布来确定系统的阶次,但结果对虚假极点剔除的程度不够明显。本文提出一种基于多次重构Hankel矩阵和层次聚类分析的ERA方法,该方法通过传统的ERA以及SVD去噪后的ERA分别求解不同维度的Hankel矩阵的模态极点,引入模态相似度准则对得到的模态极点进行初步判定,再通过层次聚类分析进一步剔除虚假模态,获取效果更好的稳定图,实现了风力机叶片模态参数更准确的辨识。

2 双目摄影测量原理

图1为双目立体成像模型,该模型的图像像素坐标系(o--uv)描述了像素在图像中的位置,其中,单位为pixel;成像平面坐标系(o--xy)是指像素坐标变换到实际物理尺寸的坐标系,其中坐标原点O为图像的主点,x和y轴分别平行于图像像素坐标系(o--uv)的u和v轴;摄像机坐标系(O_C--X_CY_CZ_C)的坐标原点O_C为摄像机的光学中心,Z_C轴与光轴重合,X_C和Y_C轴分别平行于成像平面坐标系的x和y轴。

图 1. 双目立体成像模型

Fig. 1. Binocular stereo imaging model

下载图片查看所有图片

为了描述世界坐标系中的任意一点坐标P(x_w,y_w,z_w)的转换关系,现定义:下标l和r分别代表左、右摄像机,且左摄像机坐标系(O_Cl--X_ClY_ClZ_Cl)和三维世界坐标系重合。具体转换步骤如下:右摄像机坐标系(O_Cr--X_CrY_CrZ_Cr)通过旋转矩阵R₀绕X_Cr、Y_Cr和Z_Cr轴分别旋转α、β和γ角,得到与左摄像机坐标系平行的坐标系(O_k--X_kY_kZ_k),再通过平移矩阵T(t_x,t_y,t_z)的转换,即可实现世界坐标系与右摄像机坐标系之间的转换。根据小孔成像原理,得到摄像机坐标系与成像平面坐标系之间的转换关系,最后通过转换坐标原点和单位尺寸,生成图像像素坐标。世界坐标系中的任意一点P(x_w,y_w,z_w)在三维世界坐标系、摄像机坐标系、图像平面坐标系和图像像素坐标系的空间变换关系为

\{\begin{array}{l} [\begin{array}{l} X_{C} \\ Y_{C} \\ Z_{C} \end{array}] = R_{r} \cdot [\begin{array}{l} x_{w} \\ y_{w} \\ z_{w} \end{array}] + T \\ x = f \frac{X_{C}}{Z_{C}} \\ y = f \frac{Y_{C}}{Z_{C}} \\ u_{p} = u_{0} + \frac{x_{p}}{d_{x}} \\ v_{p} = v_{0} + \frac{y_{p}}{d_{y}} \end{array}, (1)

式中:R_r为右摄像机坐标系(O_Cr--X_CrY_CrZ_Cr)依次绕X_Cr、Y_Cr和Z_Cr轴分别旋转α、β和γ角的旋转正交矩阵;T为平移矢量;f为摄像机镜头焦距;d_x和d_y分别表示每个像素在x轴和y轴的物理尺寸(单位:mm/pixel)。根据(1)式可以得到左摄像机图像像素坐标到三维世界坐标的映射关系为

Z_{Cl} [\begin{array}{l} u_{pl} \\ v_{pl} \\ 1 \end{array}] = [\begin{array}{l} 1 / d_{x} & 0 & u_{0 l} \\ 0 & 1 / d_{y} & v_{0 l} \\ 0 & 0 & 1 \end{array}] [\begin{array}{l} f_{l} & 0 & 0 \\ 0 & f_{l} & 0 \\ 0 & 0 & 1 \end{array}] [\begin{array}{l} x_{w} \\ y_{w} \\ z_{w} \end{array}] 。 (2)

根据(2)式可以得到右摄像机图像像素坐标到三维世界坐标的映射关系为

\begin{array}{l} Z_{Cr} [\begin{array}{l} u_{pr} \\ v_{pr} \\ 1 \end{array}] = [\begin{array}{l} 1 / d_{x} & 0 & u_{0 r} \\ 0 & 1 / d_{y} & v_{0 r} \\ 0 & 0 & 1 \end{array}] \times \\ [\begin{array}{l} f_{r} & 0 & 0 \\ 0 & f_{r} & 0 \\ 0 & 0 & 1 \end{array}] [\begin{array}{l} r_{1} & r_{2} & r_{3} & t_{x} \\ r_{4} & r_{5} & r_{6} & t_{y} \\ r_{7} & r_{8} & r_{9} & t_{z} \end{array}] [\begin{array}{l} x_{w} \\ y_{w} \\ z_{w} \\ 1 \end{array}] 。 (3) \end{array}

采用张正友的棋盘格标定方法^[18],得到摄像机的内部参数(焦距f和畸变系数s)和外部参数(旋转矩阵R₀和平移矢量T),联立(2)式和(3)式即可求解编码标志的三维坐标为

\{\begin{array}{l} \begin{array}{l} x_{w} = z_{w} x_{l} / f_{l} \\ y_{w} = z_{w} y_{l} / f_{l} \end{array} \\ z_{w} = \frac{f_{l} (f_{r} t_{x} - y_{r} t_{z})}{y_{r} (r_{7} x_{l} + r_{8} y_{l} + f_{l} r_{9}) - f_{r} (r_{4} x_{l} + r_{5} y_{l} + f_{l} r_{6})} \end{array}, (4)

\{\begin{array}{l} x_{l} = x_{pl} \cdot d_{x} \\ y_{l} = x_{pl} \cdot d_{y} \\ x_{r} = x_{pr} \cdot d_{x} \\ y_{r} = x_{pr} \cdot d_{y} \end{array} 。 (5)

本文的编码标志由中心圆、环形编码带以及对角扇形组成,编码标志的坐标定义为对角扇形的角点(中心圆圆心)。编码标志采用12位Schneider编码方案,根据环形编码带的位置、大小和数量识别编码标志的编码代号^[19]。在中心坐标的提取上,首先通过Canny算子得到中心圆的轮廓,然后基于最小二乘椭圆拟合法求解圆心坐标,同时结合Hough变换对扇形角点进行检测定位,最后根据圆心坐标和扇形角点坐标的差值,判定坐标提取的准确性。

3 模态参数的识别

3.1 特征系统实现算法

对于一个n维线性时不变振动系统,在时间离散域上的状态方程和输出方程可以表示为

\{\begin{array}{l} x (k + 1) = Ax (k) + Bu (k) \\ y (k) = Cx (k) \end{array}, (6)

式中:A为系统转移矩阵,A∈ℝ²ⁿ^×2ⁿ;B为输入矩阵,B∈ℝ²ⁿ^×^l,其中,l为输入的数量;C为观测矩阵,C∈ℝ^m^×2ⁿ,其中,m为输出的数量;x(k)、u(k)、y(k)分别为状态向量、输入向量和输出向量;k为离散时间。系统的响应数据可以表示为

Y (k) = C \cdot A^{k} \cdot B 。 (7)

根据系统响应数据构造广义Hankel矩阵为

\begin{array}{l} H (k - 1) = \\ [\begin{array}{l} h (k) & h (k + 1) & \dots & h (k + s - 1) \\ h (k + 1) & h (k + 2) & \dots & h (k + s) \\ ⋮ & ⋮ & ⋮ \\ h (k + r - 1) & h (k + r) & \dots & h (k + r + s - 2) \end{array}] \\ = P \cdot A^{k - 1} \cdot Q, (8) \end{array}

式中:P为能观矩阵,P∈ℝ^rm^×^rs;Q为能控矩阵,Q∈ℝ²ⁿ^×^rs;rs为系统的阶次。当k=1时,对H(0)进行奇异值分解可得

\{\begin{array}{l} H (0) = LΛ R^{T} \\ Λ = diag (λ_{1}, λ_{2}, \dots, λ_{n}, λ_{n + 1}, \dots) \end{array}, (9)

式中:L和R分别为H(0)的左、右奇异值向量矩阵,L∈ℝ^rm^×2ⁿ,R∈ℝ²ⁿ^×^sl;Λ为奇异值对角阵,Λ∈ℝ²ⁿ^×2ⁿ;λ_n为奇异值。通过设定阈值ε截断奇异值矩阵Λ,进而获得系统的最小实现阶次n',截断条件为

\{\begin{array}{l} \frac{λ_{n'}}{λ_{1}} > ε \\ \frac{λ_{n' + 1}}{λ_{1}} \leq ε \end{array} 。 (10)

系统的最小实现[A',B',C']为

\{\begin{array}{l} A' = Λ_{n'}^{- \frac{1}{2}} {L^{T}}_{n'} H (1) R_{n'} Λ_{n'}^{- \frac{1}{2}} \\ B' = Λ_{n'}^{\frac{1}{2}} {R^{T}}_{n'} {E^{T}}_{N_{i}} \\ C' = Λ_{n'}^{\frac{1}{2}} L_{n'} E_{N_{o}} \end{array}, (11)

式中:Λ_n'为Λ截断后的特征值所构成的对角阵,Λ_n'=diag(λ₁,λ₂,…,λ_n');L_n'和R_n'分别为L和R矩阵的前n'列; $E_{N_{i}}$ = $[\begin{array}{l} I_{N_{i}} & 0_{N_{i}} & \dots & 0_{N_{i}} \end{array}]$ ,I_Ni为N_i×N_i的单位矩阵; $E_{N_{o}}$ = $[\begin{array}{l} I_{N_{o}} & 0_{N_{o}} & \dots & 0_{N_{o}} \end{array}]$ ,I_No为N_o×N_o的单位矩阵。对A'进行特征值分解可解得系统的固有频率ω_i、阻尼比ξ_r及模态矩阵Φ为

\{\begin{array}{l} λ_{r} = \frac{In (z_{i})}{Δt} = λ_{r}^{Re} + j λ_{r}^{Im} \\ ω_{n} = \sqrt[]{(λ_{r}^{Re})^{2} + (λ_{r}^{Im})^{2}} \\ ξ_{r} = \frac{λ_{r}^{Re}}{\sqrt[]{(λ_{r}^{Re})^{2} + (λ_{r}^{Im})^{2}}} \\ Φ = C \times Ψ \end{array} 。 (12)

3.2 系统响应的去噪与重构

当数据采集系统用于实验测量时,被测信号不可避免地被噪声污染。在模态分析过程中,如果数据噪声污染严重,与某些微弱信号相关的真实模态很可能被噪声所掩盖,模态分析算法生成的稳定图容易出现大量的虚假模态,影响参数的辨识。为剔除原始数据中的噪声,在特征系统实现算法的基础上,采用奇异值分解的性质来估计初始Hankel矩阵的秩,对得到的奇异值进行筛选和置零操作,生成去除系统噪声的近似矩阵,最后通过Cadzow算法重构为Hankel矩阵的形式,实现系统数据的去噪与重构。

定义由离散脉冲响应函数构成的Hankel矩阵为H(0),可以写为

H (0) = [\begin{array}{l} h (1) & h (2) & \dots & h (s) \\ h (2) & h (3) & \dots & h (s + 1) \\ ⋮ & ⋮ & ⋮ \\ h (r) & h (r + 1) & \dots & h (r + s - 1) \end{array}] 。 (13)

由于测量环境的影响,此时Hankel矩阵H(0)不仅包含有用的数据信息还包含数据噪声,奇异值分解后可以写为

H_{0} = {\overset{̅}{H}}_{0} + N_{0} = U_{r \times r} S_{r \times s} {V^{T}}_{s \times s}, (14)

式中: ${\overset{̅}{H}}_{0}$ 和N₀分别表示无噪声信号与噪声信号;U和V是正交矩阵;S是其对角元素为σ_i的实对角矩阵,σ_i为H₀的奇异值。一般来说,非零奇异值的个数就是矩阵的秩,但实际测量的数据不可避免地存在随机噪声,超出矩阵秩的奇异值可能不会变为零,而是变得非常小,为了更好地表述,将(14)式写为

H_{0} = {\overset{̅}{H}}_{0} + N_{0} = [\begin{array}{l} U_{g} & U_{0} \end{array}] [\begin{array}{l} S_{g} & 0 \\ 0 & S_{0} \end{array}] [\begin{array}{l} V_{g}^{T} \\ V_{0}^{T} \end{array}], (15)

式中:S_g为大于噪声阈值ε的奇异值σ_i(i=1,2,…,g),代表没有噪声的信号部分;S₀为小于噪声阈值ε的奇异值σ_i(i=g+1,g+2,…,L),代表具有噪声的信号部分,其中,L=min(r,s)。通过绘制归一化奇异值图谱的方法来确定大于噪声阈值ε的奇异值个数g,根据(10)式即可确定合适的噪声阈值大小ε。把S₀中的奇异值全部置为0后,此时得到的矩阵 ${\dot{H}}_{0}$ 可以近似为H₀去除数据噪声后的矩阵,即

\{\begin{array}{l} {\dot{H}}_{0} = UΣ V^{T} \\ Σ = [\begin{array}{l} S_{g} & 0 \\ 0 & 0 \end{array}] \end{array}, (16)

式中:Σ为消除较小的奇异值S₀后的奇异值矩阵。需要注意的是,此时的 ${\dot{H}}_{0}$ 不再具有Hankel矩阵的结构,为了使其恢复为Hankel矩阵的结构,通过Cadzow算法将 ${\dot{H}}_{0}$ 中反对角线上的所有元素,替换为沿反对角线元素的算术平均值,计算过程如下式所示。重构后的Hankel矩阵 ${\dot{H}}_{1}$ 的每个元素可表示为

\dot{h} (i) = \frac{1}{β - α + 1} \overset{β}{\sum_{j = α}} {\dot{H}}_{i - j + 1, j}, (17)

式中:α=max(1,i-r+1);β=min(s,i)。

3.3 多次重构Hankel矩阵的虚假极点剔除法

使用特征系统实现算法识别模态参数时,如果选用的阶次低于系统真实的阶次,有可能会遗漏真实模态,如果选用的阶次高于系统真实的模态阶次,则识别的稳定图会呈现较多的虚假模态,对真实模态造成干扰。

针对该问题,定义初始的Hankel矩阵H₀为参考矩阵,以降噪后重构的Hankel矩阵 ${\dot{H}}_{1}$ 为标准矩阵,采用传统的ERA算法分别获取这两个矩阵的模态参数,并设为参考极点组K₀和标准极点组K₁。需要注意的是,这里标准极点组K₁是标准矩阵 ${\dot{H}}_{1}$ 生成的极点。以标准组极点为标准,选取标准组与参考组在同阶次情况下固有频率的最近极点,两极点间距离 ${r_{ij}}^{[3]}$ 的计算公式为

r_{ij} = \frac{W_{f}}{d f} \frac{|f_{i} - f_{j}|}{\max (f_{i}, f_{j})} + \frac{W_{ξ}}{d ξ} \frac{|ξ_{i} - ξ_{j}|}{\max (ξ_{i}, ξ_{j})}, (18)

式中:f_i和ξ_i为标准组极点对应的固有频率和阻尼比;f_j和ξ_j为参考组极点对应的固有频率和阻尼比;W_f和W_ξ为固有频率和阻尼比的权重;df和dξ为固有频率和阻尼比的容差。若r_ij≤1,则判定极点为稳定极点;反之,则为不稳定极点。通过多次改变初始Hankel矩阵的行列数和模态相似度判定,生成一个稳定模态集合,对该稳定的模态集合进行层次聚类分析,进一步剔除虚假模态,对剔除虚假模态后的结果进行均值处理,得到最终辨识的模态参数,具体的算法流程如图2所示,其中虚线标记的为传统的ERA算法。

图 2. 改进的ERA算法流程图

Fig. 2. Flow chart of improved ERA

下载图片查看所有图片

4 实验与结果分析

为了验证该方法的有效性及优越性,基于摄影测量技术对某型1.5 kW风力机叶片开展了测振实验,图3(a)为双目高速摄像测量仪,它由两个德国Optronis高速工业相机,一个LED直流光源和采集存储硬件组成,其中相机镜头的分辨率和焦距分别为4080 pixel×3072 pixel和35 mm。图3(b)为实验的约束状态。图3(c)为风力机叶片表面9号编码标志的识别效果图。图3(d)为待测量的风力机叶片,该叶片的材料由尼龙和玻璃纤维复合组成,直径为890 mm。在本次实验中,相机和风力机叶片的距离为1.5 m,叶片上共贴有15个尺寸为20 mm×20 mm的编码标志。

图 3. 摄影测量实验图。(a) 双目摄像测量仪; (b) 实验约束状态; (c) 9号编码标志识别图; (d) 风力机叶片

Fig. 3. Experiment configuration. (a) Binocular camera measuring instrument; (b) experimental constraint state; (c) identification of No.9 coded marker; (d) wind turbine blade

下载图片查看所有图片

以37、39和41号编码标志重构坐标系,并定义39号编码标志为坐标原点,以37和39号编码标志的连线为X轴,39和41号编码标志的连线为Y轴,以XY平面的法线方向为Z轴。设置相机的拍摄帧率为120 Hz,在Z方向的敲击激励下,连续采集1000幅图像,通过对编码标志的匹配跟踪,实现风力机叶片在X、Y、Z三个方向上振动数据的测量。其中图4为左右相机的标定参数和结果。图5为在敲击激励下,所有编码标志在X、Y、Z三个方向的振动响应测量结果。

图 4. 相机标定结果

Fig. 4. Camera calibration result

下载图片查看所有图片

表 2. 叶片固有频率识别误差

Table 2. Natural frequency identification error of blades unit: %

Item	Traditional ERA				Improved ERA				B&Khammer
Item	Xdirection	YDirection	Zdirection	Averageerror	Xdirection	YDirection	Zdirection	Averageerror	B&Khammer
First-order	0.18	3.66	4.77	2.87	1.04	2.01	1.12	1.39
Second-order	0.79	1.80	0.43	3.02	2.01	1.85	1.96	1.94
Third-order	5.18	5.20	2.62	4.3	1.99	2.04	2.14	2.05

查看所有表

表 1. 叶片固有频率识别结果

Table 1. Natural frequency identification results of blades unit: Hz

Item	Traditional ERA			Improved ERA			B&Khammer
Item	Xdirection	YDirection	Zdirection	Xdirection	YDirection	Zdirection	B&Khammer
First-order	6.59	6.36	6.28	6.67	6.73	6.67	6.61
Second-order	22.35	22.57	22.08	22.62	22.58	22.60	22.17
Third-order	39.87	39.86	40.95	41.21	41.19	41.15	42.05

查看所有表

图 5. 编码标志的振动位移曲线。(a) X方向; (b) Y方向; (c) Z方向

Fig. 5. Vibration displacement curves of coded marker. (a) X direction; (b) Y direction; (c) Z direction

下载图片查看所有图片

为避免生成过多的虚假极点,设置模型迭代阶次为20,并将编码标志的振动响应数据导入ERA辨识算法中,以频率为横坐标,绘制系统的稳定图并求解模态参数。图6分别列举了在敲击激励下,传统的ERA算法辨识的稳定图,由图可知,尽管Z方向振动数据的振幅较X方向和Y方向更明显,但其生成的模态却近似相同。图7为重构无噪声污染的Hankel矩阵 ${\dot{H}}_{1}$ 后生成的ERA稳定图。通过对比图6和图7可知,对初始Hankel矩阵进行SVD分解去噪,可以减少大部分由数据噪声生成的孤立点,并可使部分虚假极点向主模态靠近,如图6圈标记所示。但由于数学模态的存在,部分虚假模态仍然无法完全剔除,如图7圈标记所示。对多次重构Hankel矩阵后所得到的稳定极点进行层次聚类分析,以阻尼和频率的形式所表示的聚类结果如图8所示。进一步剔除虚假模态后的稳定图如图9所示。需要注意的是,在本次识别过程中,Hankel矩阵重构的次数为25,固有频率和阻尼比的权重设置为0.2,两者的容差分别设置为0.01和0.04。其中,图10为SVD分解去噪时所绘制的归一化奇异值图。在三个方向数据的模态辨识中,保留的奇异值个数分别为29、34和28,选取的噪声阈值分别为0.2%、0.5%和0.2%。

图 6. 传统的ERA稳定图。(a) X方向; (b) Y方向; (c) Z方向

Fig. 6. Stabilization diagrams for traditional ERA. (a) X direction; (b) Y direction; (c) Z direction

下载图片查看所有图片

图 7. SVD去噪声后的ERA稳定图。 (a) X方向; (b) Y方向; (c) Z方向

Fig. 7. Stabilization diagrams after noise removal by SVD. (a) X direction; (b) Y direction; (c) Z direction

下载图片查看所有图片

图 8. 阻尼频率关系图。(a) X方向; (b) Y方向; (c) Z方向

Fig. 8. Diagrams of damping and frequency. (a) X direction; (b) Y direction; (c) Z direction

下载图片查看所有图片

图 9. 改进后的ERA稳定图。 (a) X方向; (b) Y方向; (c) Z方向

Fig. 9. Stabilization diagrams for improved ERA. (a) X direction; (b) Y direction; (c) Z direction

下载图片查看所有图片

图 10. 归一化的奇异值。 (a) X方向; (b) Y方向; (c) Z方向

Fig. 10. Normalized singular values. (a) X direction; (b) Y direction; (c) Z direction

下载图片查看所有图片

为了验证本文方法的准确性和有效性,采用B&K公司的8206-002力锤、4514B-001加速度传感器和3053-B-120采集卡对某型1.5 kW风力机叶片开展模态力锤实验,首先使用PULSE软件的MTC Hammer模块对实验对象进行建模与信号采集,最后通过Reflex后处理模块进行模态参数的辨识,并取10次实验结果的平均值作为最终得到的结果。图11为B&K仪器辨识的稳定图。

图 11. B&K仪器识别的稳定图

Fig. 11. Stabilization diagram for B&K instrument

下载图片查看所有图片

表1列出了ERA方法辨识的前三阶固有频率。表2列出了叶片模态频率识别的误差,需要注意的是,各阶固有频率都是取间隔理论极点为1.5 Hz的全部实际极点的平均值作为最终结果,各点的具体取值是通过Matlab中的数据游标工具得到。计算方法为

{\overset{̅}{f}}_{e} = \frac{\overset{n}{\sum_{i = 1}} f_{i, e}}{n}, f_{i, e} \in (f_{e}^{*} - ε, f_{e}^{*} + ε), (19)

式中:f_i_,_e为稳定图上实际生成的极点值; $f_{e}^{*}$ 为理论极点值; ${\overset{̅}{f}}_{e}$ 为实际极点的平均值,其中e代表阶数;ε=1.5 Hz。表2中各阶固有频率的相对误差,都是以B&K实验的结果作为真实值,计算公式为

δ_{e} = \frac{|{\overset{̅}{f}}_{e} - f_{b, e}|}{f_{b, e}}, (20)

式中:f_b_,_e为锤击法实验得到的结果;δ_e为固有频率的相对误差。由表1可知,虽然编码标志在X、Y、Z三个方向的振幅不同,生成的稳定图也不同,但辨识的结果和B&K力锤实验测得的结果相近。由表2可知,改进的ERA方法比传统的ERA方法辨识的结果更可靠,各阶固有频率的相对误差均小于3%,三个方向数据辨识的各阶固有频率的平均误差在2%左右,其中,Z方向数据辨识的第一阶固有频率相较传统的ERA方法提高了3.6%的精度。实验结果表明,基于摄影测量获取的编码标志的振动响应数据可靠,基于奇异值分解去噪与多次重构Hankel矩阵的ERA改进算法能够准确地识别模态参数,获得识别效果更好的稳定图。

5 结论

针对传统接触式测量方法在大尺寸结构应用中的局限性,本文通过摄影测量的方法获取编码标志在空间的运动信息,克服了接触式测量方法采集点数量不够以及传感器附加质量改变结构特性的问题。针对模态计算过程中容易出现虚假模态的问题,提出了一种改进的ERA方法,获取了更清晰和更精确的稳定图,实现了风力机叶片模态参数更准确的辨识。为了证明该方法的有效性,对某型1.5 kW风力机叶片进行实验验证。通过分析发现,风力机叶片在X、Y、Z三个方向的振动数据,均可单独用于辨识模态参数。数据噪声会导致主模态发生偏移,并生成散乱分布的虚假极点,虚假模态虽然不能完全消除,但是可以采用层次聚类分析进行剔除。实验结果表明,摄影测量方法可以准确地获取风力机叶片表面编码标志的空间运动信息,改进的ERA方法识别的各阶固有频率的相对误差均小于3%,三个方向数据辨识的各阶固有频率的平均误差在2%左右,其中使用Z方向数据辨识的模态参数的改进效果最明显,研究成果为大型结构的模态测试和叶片的状态监测提供了一种参考方法。

参考文献

[1] Li D S. Ho S C M, Song G B, et al. A review of damage detection methods for wind turbine blades[J]. Smart Materials and Structures, 2015, 24(3): 033001.

[2] Juang J N, Pappa R S. An eigensystem realization algorithm for modal parameter identification and model reduction[J]. Journal of Guidance, Control, and Dynamics, 1985, 8(5): 620-627.

[3] 张永祥, 刘心, 褚志刚, 等. 基于随机子空间法的模态参数自动提取[J]. 机械工程学报, 2018, 54(9): 187-194.

Zhang Y X, Liu X, Chu Z G, et al. Autonomousmodal parameter extraction based on stochastic subspace identification[J]. Journal of Mechanical Engineering, 2018, 54(9): 187-194.

[4] Ozbek M, Meng F Z, Rixen D J. Challenges in testing and monitoring the in-operation vibration characteristics of wind turbines[J]. Mechanical Systems and Signal Processing, 2013, 41(1/2): 649-666.

[5] Wang W Y, Chen A H. Target-less approach of vibration measurement with virtual points constructed with cross ratios[J]. Measurement, 2020, 151: 107238.

[6] 杨谢柳, 尹晨宇, 方素平, 等. 基于全站仪的大型三维形貌摄影测量方法[J]. 激光与光电子学进展, 2020, 57(10): 101505.

Yang X L, Yin C Y, Fang S P, et al. Three-dimensional measurement of large-scale objects using photogrammetry based on total station[J]. Laser & Optoelectronics Progress, 2020, 57(10): 101505.

[7] 王君, 董明利, 李巍, 等. 大型槽式太阳能反射镜面摄影测量方法[J]. 激光与光电子学进展, 2018, 55(5): 051204.

Wang J, Dong M L, Li W, et al. Photogrammetric method for large trough solar reflector[J]. Laser & Optoelectronics Progress, 2018, 55(5): 051204.

[8] 严俊, 叶南, 李廷成, 等. 无编码点的工业摄影测量技术的研究及实现[J]. 光学学报, 2019, 39(10): 1015002.

Yan J, Ye N, Li T C, et al. Research and implementation of industrial photogrammetry without coded points[J]. Acta Optica Sinica, 2019, 39(10): 1015002.

[9] 解则晓, 王晓东, 宫韩磊. 基于双目视觉的薄壁零件圆孔轮廓测量[J]. 中国激光, 2019, 46(12): 1204004.

Xie Z X, Wang X D, Gong H L. Round-hole profile measurement of thin-walled parts based on binocular vision[J]. Chinese Journal of Lasers, 2019, 46(12): 1204004.

[10] Li X, Jiang T, Yuan H Y, et al. An eigensystem realization algorithm based data-driven approach for extracting electromechanical oscillation dynamic patterns from synchrophasor measurements in bulk power grids[J]. International Journal of Electrical Power & Energy Systems, 2020, 116: 105549.

[11] Castellanos-Toro S, Marmolejo M, Marulanda J, et al. Frequencies and damping ratios of bridges through operational modal analysis using smartphones[J]. Construction and Building Materials, 2018, 188: 490-504.

[12] Jiang Y, Kan J, Zhang Z, et al. Machinefault diagnosis based on reassigned wavelet scalogram and SVD[J]. International Journal of Advancements in Computing Technology, 2012, 4(20): 10-21.

[13] GolafshanR, Yuce Sanliturk K.SVD and Hankel matrix based de-noising approach for ball bearing fault detection and its assessment using artificial faults[J]. Mechanical Systems and SignalProcessing, 2016, 70/71: 36- 50.

[14] Cadzow J A. Signal enhancement-a composite property mapping algorithm[J]. IEEE Transactions on Acoustics, Speech, and Signal Processing, 1988, 36(1): 49-62.

[15] Zhang Y, Lian J J, Liu F. Modal parameter identification for a roof overflow powerhouse under ambient excitation[J]. Water Science and Engineering, 2016, 9(1): 67-80.

[16] Moaveni B, Masoumi Z. Modifying the ERA and fast ERA to improve operational performance for structural system identification[J]. Mechanical Systems and Signal Processing, 2019, 120: 664-692.

[17] Chiang D Y, Lin C S. Identification of modal parameters from ambient vibration data using eigensystem realization algorithm with correlation technique[J]. Journal of Mechanical Science and Technology, 2010, 24(12): 2377-2382.

[18] Zhang Z. A flexible new technique for camera calibration[J]. IEEE Transactions on Pattern Analysis and Machine Intelligence, 2000, 22(11): 1330-1334.

[19] 张宗华, 王森, 王宇莹, 等. 基于编码环带径向直线拟合的圆环编码标志点中心提取方法[J]. 激光与光电子学进展, 2020, 57(7): 071203.

Zhang Z H, Wang S, Wang Y Y, et al. Method for extracting the center of circular coded target based on radial straight line fittingof circular coding[J]. Laser & Optoelectronics Progress, 2020, 57(7): 071203.

3.3 多次重构Hankel矩阵的虚假极点剔除法

4 实验与结果分析

5 结论

全伟铭, 郭迎福, 王文韫, 邹龙洲. 基于摄影测量和ERA改进算法的叶片模态辨识方法[J]. 激光与光电子学进展, 2021, 58(4): 0411001. Weiming Quan, Yingfu Guo, Wenyun Wang, Longzhou Zou. Modal Identification Method of Blade Based on Photogrammetry and Improved ERA[J]. Laser & Optoelectronics Progress, 2021, 58(4): 0411001.

基于摄影测量和ERA改进算法的叶片模态辨识方法下载： 875次

1 引言

2 双目摄影测量原理

图 1. 双目立体成像模型

Fig. 1. Binocular stereo imaging model

3 模态参数的识别

3.1 特征系统实现算法

3.2 系统响应的去噪与重构

3.3 多次重构Hankel矩阵的虚假极点剔除法

图 2. 改进的ERA算法流程图

Fig. 2. Flow chart of improved ERA

4 实验与结果分析

图 3. 摄影测量实验图。(a) 双目摄像测量仪; (b) 实验约束状态; (c) 9号编码标志识别图; (d) 风力机叶片

Fig. 3. Experiment configuration. (a) Binocular camera measuring instrument; (b) experimental constraint state; (c) identification of No.9 coded marker; (d) wind turbine blade

图 4. 相机标定结果

Fig. 4. Camera calibration result

表 2. 叶片固有频率识别误差

Table 2. Natural frequency identification error of blades unit: %

表 1. 叶片固有频率识别结果

Table 1. Natural frequency identification results of blades unit: Hz

图 5. 编码标志的振动位移曲线。(a) X方向; (b) Y方向; (c) Z方向

Fig. 5. Vibration displacement curves of coded marker. (a) X direction; (b) Y direction; (c) Z direction

图 6. 传统的ERA稳定图。(a) X方向; (b) Y方向; (c) Z方向

Fig. 6. Stabilization diagrams for traditional ERA. (a) X direction; (b) Y direction; (c) Z direction

图 7. SVD去噪声后的ERA稳定图。 (a) X方向; (b) Y方向; (c) Z方向

Fig. 7. Stabilization diagrams after noise removal by SVD. (a) X direction; (b) Y direction; (c) Z direction

图 8. 阻尼频率关系图。(a) X方向; (b) Y方向; (c) Z方向

Fig. 8. Diagrams of damping and frequency. (a) X direction; (b) Y direction; (c) Z direction

图 9. 改进后的ERA稳定图。 (a) X方向; (b) Y方向; (c) Z方向

Fig. 9. Stabilization diagrams for improved ERA. (a) X direction; (b) Y direction; (c) Z direction

图 10. 归一化的奇异值。 (a) X方向; (b) Y方向; (c) Z方向

Fig. 10. Normalized singular values. (a) X direction; (b) Y direction; (c) Z direction

图 11. B&K仪器识别的稳定图

Fig. 11. Stabilization diagram for B&K instrument

5 结论

Article Outline

关于本站 Cookie 的使用提示

全站搜索

基于摄影测量和ERA改进算法的叶片模态辨识方法 下载： 875次

1 引言

2 双目摄影测量原理

图 1. 双目立体成像模型

Fig. 1. Binocular stereo imaging model

3 模态参数的识别

3.1 特征系统实现算法

3.2 系统响应的去噪与重构

3.3 多次重构Hankel矩阵的虚假极点剔除法

图 2. 改进的ERA算法流程图

Fig. 2. Flow chart of improved ERA

4 实验与结果分析

图 3. 摄影测量实验图。(a) 双目摄像测量仪; (b) 实验约束状态; (c) 9号编码标志识别图; (d) 风力机叶片

Fig. 3. Experiment configuration. (a) Binocular camera measuring instrument; (b) experimental constraint state; (c) identification of No.9 coded marker; (d) wind turbine blade

图 4. 相机标定结果

Fig. 4. Camera calibration result

表 2. 叶片固有频率识别误差

Table 2. Natural frequency identification error of blades unit: %

表 1. 叶片固有频率识别结果

Table 1. Natural frequency identification results of blades unit: Hz

图 5. 编码标志的振动位移曲线。(a) X方向; (b) Y方向; (c) Z方向

Fig. 5. Vibration displacement curves of coded marker. (a) X direction; (b) Y direction; (c) Z direction

图 6. 传统的ERA稳定图。(a) X方向; (b) Y方向; (c) Z方向

Fig. 6. Stabilization diagrams for traditional ERA. (a) X direction; (b) Y direction; (c) Z direction

图 7. SVD去噪声后的ERA稳定图。 (a) X方向; (b) Y方向; (c) Z方向

Fig. 7. Stabilization diagrams after noise removal by SVD. (a) X direction; (b) Y direction; (c) Z direction

图 8. 阻尼频率关系图。(a) X方向; (b) Y方向; (c) Z方向

Fig. 8. Diagrams of damping and frequency. (a) X direction; (b) Y direction; (c) Z direction

图 9. 改进后的ERA稳定图。 (a) X方向; (b) Y方向; (c) Z方向

Fig. 9. Stabilization diagrams for improved ERA. (a) X direction; (b) Y direction; (c) Z direction

图 10. 归一化的奇异值。 (a) X方向; (b) Y方向; (c) Z方向

Fig. 10. Normalized singular values. (a) X direction; (b) Y direction; (c) Z direction

图 11. B&K仪器识别的稳定图

Fig. 11. Stabilization diagram for B&K instrument

5 结论

Article Outline

相关论文

相关资讯

关于本站 Cookie 的使用提示

全站搜索

基于摄影测量和ERA改进算法的叶片模态辨识方法下载： 875次