CN103323209B - 基于双目立体视觉的结构模态参数识别系统 - Google Patents

Abstract

本发明公开了一种基于双目立体视觉识别结构模态参数的系统。其包括：搭建双目立体视觉系统，标定获取两个摄像单元的内外参数，进行三维立体校正；对被测结构表面建立易识别的关键点；加载冲击力让被测结构自由运动，获取同步的图像序列文件，依次对每一帧图像识别所有关键点的三维坐标，最终得到整个测试时间内所有关键点的三维运动数据；通过所有关键点的三维运动数据使用特征系统实现算法进行模态识别，获得被测结构的模态参数，输出结果。与传统技术相比，本发明能够非接触地或对被测结构施加较小影响地测量被测结构的模态参数，其成本低、速度快、精度高，尤其适合不方便布置传感器的大型质轻结构的模态参数测量。

Description

基于双目立体视觉的结构模态参数识别系统

技术领域

本发明属于图像视觉领域和动力学模态参数辨识的交叉领域，涉及双目立体视觉系统的建立、关键点的识别和三维运动数据的计算、结构模态参数的识别算法。

背景技术

结构的模态参数识别是结构动态性能设计、分析的重要手段，在线模态参数识别对于结构振动主动控制、损伤识别、故障预报及诊断、结构健康监测等具有重要的应用价值。现代工程结构的发展趋于轻量化、大型化，部分结构甚至工作在极端环境中，如高温、外太空等。传统基于接触式测量的振动测试由于存在影响轻质结构模态参数、不易在极端环境中测量、需大范围布置数据线路等问题，不再适用于某些现代工程结构的模态参数识别，如高温高转速下工作的涡轮叶片、空间站的大面积太阳翼。

相比于接触式测量，非接触测量对被测对象的附加影响极小，因此更适合于柔软、易变形结构的测量。视觉测量技术具有非接触、低成本、精度高、测量设备简单等突出特点，因此日益受到重视。鉴于此，本发明意在通过双目立体视觉系统自动识别轻质、极端环境下工作结构的模态参数，对结构主动振动控制、在线健康监测提供一种技术手段。

发明内容

本发明利用双摄像机拍摄待测结构的自由运动或受冲击载荷后的运动，识别图像中的关键点，计算得到关键点的三维运动数据，通过关键点的三维运动数据计算得到结构的模态，最终呈现给使用者。其成本低、速度快、精度高，对测量结构模态无影响或影响很小。

本发明通过如下的技术方案实现。

一种基于双目立体视觉的非接触式模态测量系统，包括摄像装置、数据采集和处理装置；

所述摄像装置包括两个摄像单元和一个触发信号发生装置，所述两个摄像单元基本以平行方式设置，左右相隔一定距离，固定在支撑单元上；两个摄像单元能够在触发信号发生装置的同步下对被测结构进行同步的稳定帧率的拍摄，从而生成同步图像序列文件；

所述数据采集和处理装置能够从所述摄像装置获得所述同步图像序列文件，并进行处理，获得所述被测结构的结构模态参数。

优选地，还包括对所述被测结构加载冲击载荷的装置。

优选地，所述摄像装置的所述两个摄像单元在所述被测结构加载冲击载荷之后同步地进行拍摄。

优选地，所述结构模态参数为各阶模态的无阻尼频率、阻尼比、阵型、模态置信度(MAC)、模态相位共线性(MPC)中的一项或多项。

一种使用根据以上技术方案之一所述的基于双目立体视觉的非接触式模态测量系统进行模态测量的方法，其特征在于，包括如下步骤：

步骤一：对于两个摄像单元进行标定，获得两个摄像单元各自的内参数和它们之间的外参数；通过三维立体校正，获得图像的消除畸变的行对齐映射和重投影矩阵。

步骤二：建立被测结构表面的关键点，加载冲击载荷让被测结构自由运动，同时使用标定完的双摄像单元摄像被测结构受冲击载荷后的自由运动，获得同步图像序列文件。

步骤三：处理摄像单元的同步图像序列文件，获得测试时间内，关键点的三维运动数据。

步骤四：通过所有关键点的三维运动数据，以特征系统实现算法识别被测结构的模态。

步骤五：根据识别结果，图形化显示各阶模态的无阻尼频率、阻尼比、阵型、模态置信度(MAC)、模态相位共线性(MPC)。

优选地，所述步骤一通过标定获得两个摄像单元分别的内参数矩阵、畸变向量和摄像单元间的外参数矩阵。

优选地，所述步骤二中的关键点的图案为黑白同心圆环。

附图说明

图1是根据本发明的结构模态参数识别系统示意图；

图2是根据本发明的结构模态参数识别系统流程图；

图3是获取关键点三维运动数据处理流程图；

图4是关键点形状示意图；

图5是特征系统实现算法计算流程图。

具体实施方式

为了更加明晰的阐述本发明的技术方案和内容，下面结合附图对本发明做进一步详细描述。

如图1所示，本发明的结构模态参数识别系统包括摄像装置、数据采集和处理装置。所述摄像装置包括两个摄像单元，基本以平行方式设置，左右相距90mm，固定在支撑单元上。

图2示出的是根据本发明进行结构模态参数识别的流程图。其中的具体流程包括如下步骤。

步骤1，搭建双目立体视觉系统，进行标定，获得两个摄像单元的内参数和相互的外参数；通过三维立体校正，获得图像的消除畸变的行对齐映射和重投影矩阵。

本发明中使用的摄像单元均满足针孔相机模型。记Q为某点在摄像单元坐标系中的三维坐标；为它在此摄像单元上的理想成像位置(齐次坐标表示)，x_p和y_p是图像上的x坐标和y坐标(x轴为水平方向，向右为正，y轴为竖直方向，向上为正)；为考虑成畸变后真实的位置，x_d和y_d是图像上的x坐标和y坐标；M为内参数矩阵，f_x，f_y，c_x，c_y为摄像单元的x向焦距、y向焦距和光轴中心坐标；F为畸变向量；T，R为右摄像单元相对左摄像单元的平移向量和旋转向量；s是一个尺度比例。

$Q=\left[\begin{array}{c}X\\ Y\\ Z\end{array}\right],{\tilde{x}}_p=\left[\begin{array}{c}{x}_p\\ y_p\\ 1\end{array}\right],{\tilde{x}}_d=\left[\begin{array}{c}{x}_d\\ y_d\end{array}\right],M=\left[\begin{array}{ccc}{f}_x& 0& c_x\\ 0& f_y& c_y\\ 0& 0& 1\end{array}\right]$

F＝[k₁k₂p₁p₂k₃]

T＝[t_xt_yt_z]^T

R＝[R₁R₂R₃]

这些参数间满足关系式：

${\tilde{x}}_p=sMQ$

$\left[\begin{array}{c}{x}_p\\ y_p\end{array}\right]=(1+k_1{r}^2+k_2{r}^4+k_3{r}^6)\left[\begin{array}{c}{x}_d\\ y_d\end{array}\right]+\left[\begin{array}{c}2p_1{x}_d{y}_d+p_2(r^2+2{x_d}^2)\\ 2p_2{x}_d{y}_d+p_1(r^2+2{y_d}^2)\end{array}\right]$

标定的目的是获得两个摄像单元分别的内参数矩阵M、畸变向量F和摄像单元间的外参数矩阵T，R。本发明采用由张正友提出的棋盘格标定方法。通过变换棋盘的位置和角度拍摄20对图像，完成后进行标定，获得内外参数：M_left，M_rjght，F_left，，F_right，R，T。

标定完成后，两摄像单元的相对位置、焦距等不再改变，否则需要重新进行标定。

标定完成后，进行三维立体校正，使两个摄像单元成像平面落在同一个平面，并完全地前向平行。具体实施使用Bouguet立体校正算法，使两图像中的每一副重投影次数最小化，同时也使得观测面积最大化，输出得到消除畸变的行对齐图像。通过Bouguet算法，最终获得左右图像的图像行对齐映射和重投影矩阵Q。

$Q=\left[\begin{array}{cccc}1& 0& 0& -c_x\\ 0& 1& 0& -c_y\\ 0& 0& 0& f\\ 0& 0& -\frac{1}{T_x}& 0\end{array}\right]$

步骤2，建立待测结构上的关键点，利用双目立体摄像系统拍摄待测结构加载冲击载荷后的自由运动，获得同步图像序列。

如果被测结构上有明显的可识别的关键点，则不需要专门建立关键点；若无明显可识别关键点，则粘贴易识别的关键点。在本实施方案中，采用粘贴易识别的黑白同心圆环，参见附图4，圆环尺寸视被测结构大小确定，外环直径与内环直径之比为2。关键点数量按结构复杂程度选取，一般沿两个方向各取10个左右。

关键点图案在本实施方案中采用黑白同心圆环，但本发明不局限于使用此特征图像，也可采用其它易识别的特征图案，或是直接识别待测结构上易识别的特征点。

拍摄图像序列时，保证在整个测试过程中所有关键点均在两摄像单元的拍摄区域内，并尽量让关键点正对镜头。

步骤3，处理摄像单元的同步图像序列文件，获得测试时间内，关键点的三维运动数据。

对同一时刻的左右两张图片，进行步骤3.1至3.7的处理，获得关键点集的屏幕坐标

$p_{left}^{\left(i\right)}={\left[\begin{array}{cc}{x}_{left}^{\left(i\right)}& y_{left}^{\left(i\right)}\end{array}\right]}^T,p_{right}^{\left(i\right)}={\left[\begin{array}{cc}{x}_{right}^{\left(i\right)}& y_{right}^{\left(i\right)}\end{array}\right]}^T$

上标i表示第i个关键点，两个分量为关键点在校正图像中的x坐标、y坐标。

步骤3.1，对两张图片使用图像行对齐映射处理得到消除畸变的行对齐图像。

步骤3.2，对两张图片使用高斯模糊去除图像中的噪点。

步骤3.3，对两张图片使用二值阀值化，将图像转化为黑白图像。

步骤3.4，求取二值阀值化后的图像的轮廓。

步骤3.5，分别对两张图片的所有轮廓执行剔除操作。

(1)周长条件：轮廓的周长∈[L_min，L_max]

(2)面积条件：轮廓的实际面积∈[S_min，S_max]

(3)面积周长比条件：轮廓的实际面积比周长＜ratio

(4)拟合误差条件：(轮廓的实际面积-轮廓拟合成椭圆的面积)/轮廓拟合面积∈[-e_e，e_e]

(5)当前轮廓数量如果小于2倍于需要的关键点数，则暂停计算，提示用户更改参数，重新进行剔除。

(6)匹配轮廓对，当轮廓的圆心相距不超过d_D时，匹配为同一对轮廓。

(7)计算匹配的轮廓对数量。如果轮廓对数量不等于关键点数量，则暂停计算，提示用户更改参数，重新进行剔除。

(8)对匹配出来的轮廓对求外圆的圆心坐标。

(9)与上一帧的位置匹配，保存结果。

以上参数L_min，L_max，S_min，S_max，ratio，e_e，d_D由用户根据实验情况设定。

步骤3.6，对左右两张图片的轮廓对匹配，计算得到所有关键点集的屏幕坐标。

步骤3.7，针对关键点集的轮廓，使用Zernike矩对轮廓边界进行亚像素定位，之后重新进行椭圆拟合，计算拟合椭圆的中心，更新相应的屏幕坐标。

传统的轮廓检测算法只能检测出像素级边缘，因而使得最终得到的椭圆中心定位精度较差；采用Zernike模板矩，可以计算轮廓的亚像素边界，椭圆中心定位精度较高。其步骤为：

计算三个不同阶次的矩，分别为：

$A_{00}=\∫{\∫}_{x^2+y^2\≤1}f(x,y)dxdy$

$A_{11}=\∫{\∫}_{x^2+y^2\≤1}f(x,y)(x-yj)dxdy$

$A_{20}=\∫{\∫}_{x^2+y^2\≤1}f(x,y)(x^2+y^2-1)dxdy$

$\mathrm{\φ}={\tan }^{-1}\left(\frac{\mathrm{Im}\left(A_{11}\right)}{\mathrm{Re}\left(A_{11}\right)}\right)$

$l=\frac{A_{20}}{A_{11}^{\′}}$

轮廓的亚像素边界点为：

$\left[\begin{array}{c}{x}_{sub}\\ y_{sub}\end{array}\right]=\left[\begin{array}{c}{x}_o\\ y_o\end{array}\right]+l\left[\begin{array}{c}cos\mathrm{\φ}\\ sin\mathrm{\φ}\end{array}\right]$

步骤3.8，根据 $Q\left[\begin{array}{c}{x}_{left}^{\left(i\right)}\\ y_{left}^{\left(i\right)}\\ x_{right}^{\left(i\right)}-x_{left}^{\left(i\right)}\\ 1\end{array}\right]=\left[\begin{array}{c}{X}^{\left(i\right)}\\ Y^{\left(i\right)}\\ Z^{\left(i\right)}\\ W^{\left(i\right)}\end{array}\right],$ 计算得到所有关键点集的三维坐标： $\left[\begin{array}{ccc}\frac{X^{\left(i\right)}}{W^{\left(i\right)}}& \frac{Y^{\left(i\right)}}{W^{\left(i\right)}}& \frac{Z^{\left(i\right)}}{W^{\left(i\right)}}\end{array}\right];$

步骤3.9，输出测试时间区间内的所有关键点的三维运动数据。

步骤4，根据被测结构的关键点三维运动数据使用特征系统实现算法方法识别结构的模态，输出并保存结果。

对于n维线性时不变系统，离散时间状态空间方程为：

x(k+1)＝Ax(k)+Bu(k)

y(k)＝Cx(k)

其中x为n×1状态向量，u为m×1控制向量，y为p×1观测向量，A为n×n系统矩阵，B为n×m维控制矩阵，C为p×n维观测矩阵对于脉冲响应：

Y(k)＝CA^k-1B，k＝1，2，3，...

Hankel矩阵定义为

步骤4.1，由双目视觉系统测量得到的运动数据，给定块矩阵行数r和列数s，分别构造零阶汉克尔矩阵H_rs(0)和一阶汉克尔矩阵H_rs(1)；

步骤4.2，对H_rs(0)进行奇异值分解(SVD)，H_rs(0)＝PDQ^T；

步骤4.3，根据给定阀值ε得到系统的阶次n，则

$H_{rs}\left(0\right)={\mathrm{PDQ}}^T=P\left[\begin{array}{cc}{D}_n& 0\\ 0& 0\end{array}\right]Q^T,$ D_n＝diag(d₁d₂…d_n)

步骤4.4，根据 $A=D_n^{-1/2}{P}_n^T{H}_{rs}\left(1\right)Q_n{D}_n^{-1/2},B=D_n^{1/2}{Q_n}^T{E}_m,C=E_p^T{P}_n{D}_n^{1/2}$ 得到系统最小实现，其中P_n和Q_n分别为P和Q的前n列，

$E_m^T=\left[\begin{array}{cc}{I}_m& 0_{m\×\left(ms-m\right)}\end{array}\right],E_p^T=\left[\begin{array}{cc}{I}_p& 0_{p\×\left(pr-p\right)}\end{array}\right]$

步骤4.5，对系统矩阵A进行特征值分解A＝VDV^-1。D＝diag(z₁z₂…z_n)，为特征值对角阵；为特征向量矩阵；

步骤4.6，计算各阶模态参数：无阻尼频率ω_i、阻尼比振型、模态置信度(MAC)、模态相位共线性(MPC)；

${\mathrm{\λ}}_i=\frac{ln\left(z_i\right)}{T},i=1,2,...,n$

${\mathrm{\ω}}_i=\sqrt{\mathrm{Re}{\left({\mathrm{\λ}}_i\right)}^2+Im{\left({\mathrm{\λ}}_i\right)}^2}$

${\mathrm{MAC}}_i=\frac{\left|q_i^{*}\·{\stackrel{\‾}{q}}_i^{*}\right|}{\sqrt{\left|{\stackrel{\‾}{q}}_i^{*}\·{\stackrel{\‾}{q}}_i^{*}\right|\left|q_i^{*}\·q_i^{*}\right|}}$

对c_i进行归一化，使||c_i||₂＝1。设θ₀是振幅最大的点对应的相位，对c_i的每个分量旋转角度-θ₀，

${c_i}^{*}={\left(\begin{array}{cccc}{c}_{i1}{e}^{-{\mathrm{\θ}}_{0}j}& c_{i2}{e}^{-{\mathrm{\θ}}_{0}j}& \mathrm{...}& c_{ip}{e}^{-{\mathrm{\θ}}_{0}i}\end{array}\right)}^T$

${\mathrm{MPC}}_i=\frac{\left|\right|c_i^{*}|{|}_2}{\left|\right|c_i|{|}_2}=\left|\right|c_i^{*}|{|}_2$

二范数归一化振型，

步骤4.7，输出并保存结果。

步骤5：根据用户选择显示各阶模态无阻尼频率、阻尼比、阵型、MAC、MPC。

以上所述的实施例，只是本发明较优选的具体实施方式的一种，本领域的技术人员在本发明技术方案范围内进行的通常变化和替换都应包含在本发明的保护范围内。

Claims (7)

1.一种使用基于双目立体视觉的结构模态参数识别系统进行模态测量的方法，所述基于双目立体视觉的结构模态参数识别系统包括摄像装置、数据采集和处理装置；

所述摄像装置包括两个摄像单元和一个触发信号发生装置，所述两个摄像单元基本以平行方式设置，左右相隔一定距离，固定在支撑单元上；两个摄像单元能够在所述触发信号发生装置的同步下对被测结构进行同步拍摄，从而生成同步图像序列文件；

所述数据采集和处理装置能够从所述摄像装置获得所述同步图像序列文件，并进行处理，获得所述被测结构的结构模态参数；

其特征在于，所述方法包括如下步骤：

步骤一：对两个摄像单元进行标定，获得两个摄像单元各自的内参数和它们之间的外参数；通过三维立体校正，获得图像的消除畸变的行对齐映射和重投影矩阵；

步骤二：建立被测结构表面的关键点，加载冲击载荷让被测结构自由运动，同时使用标定完的双摄像单元摄像被测结构受冲击载荷后的自由运动，获得同步图像序列文件；

步骤三：处理摄像单元的同步图像序列文件，获得测试时间内，关键点的三维运动数据；

步骤四：通过所有关键点的三维运动数据，以特征系统实现算法识别被测结构的模态；

步骤五：根据识别结果，图形化显示各阶模态的无阻尼频率、阻尼比、阵型、模态置信度(MAC)和模态相位共线性(MPC)；

所述步骤一中，记Q为某点在摄像单元坐标系中的三维坐标；为它在此摄像单元上的以齐次坐标表示的理想成像位置，x_p和y_p是图像上的水平方向x坐标和竖直方向y坐标，其中水平方向以向右为正，竖直方向以向上为正；为考虑成畸变后真实的位置，x_d和y_d是图像上的x坐标和y坐标；M为内参数矩阵，f_x，f_y，c_x，c_y为摄像单元的x向焦距、y向焦距和光轴中心坐标；F为畸变向量；T，R为右摄像单元相对左摄像单元的平移向量和旋转向量；s是一个尺度比例：

$Q=\left[\begin{array}{c}X\\ Y\\ Z\end{array}\right],{\tilde{x}}_p=\left[\begin{array}{c}{x}_p\\ y_p\\ 1\end{array}\right],{\tilde{x}}_d=\left[\begin{array}{c}{x}_d\\ y_d\end{array}\right],M=\left[\begin{array}{ccc}{f}_x& 0& c_x\\ 0& f_y& c_y\\ 0& 0& 1\end{array}\right]$

F＝[k₁k₂p₁p₂k₃]

T＝[t_xt_yt_z]^T

R＝[R₁R₂R₃]

这些参数间满足关系式：

${\tilde{x}}_p=sMQ$

$\left[\begin{array}{c}{x}_p\\ y_p\end{array}\right]=(1+k_1{r}^2+k_2{r}^4+k_3{r}^6)\left[\begin{array}{c}{x}_d\\ y_d\end{array}\right]+\left[\begin{array}{c}2p_1{x}_d{y}_d+p_2\left(r^2+2{x_d}^2\right)\\ 2p_2{x}_d{y}_d+p_1\left(r^2+2{y_d}^2\right)\end{array}\right]$

通过标定获得内外参数：M_left，M_right，F_left，F_right，R，T。

2.根据权利要求1所述的方法，其特征在于，所述基于双目立体视觉的结构模态参数识别系统还包括对所述被测结构加载冲击载荷的装置。

3.根据权利要求2所述的方法，其特征在于，所述摄像装置的所述两个摄像单元在所述被测结构加载冲击载荷之后同步地以稳定帧率进行拍摄。

4.根据权利要求3所述的方法，其特征在于，所述结构模态参数为各阶模态的无阻尼频率、阻尼比、阵型、模态置信度(MAC)、模态相位共线性(MPC)中的一项或多项。

5.根据权利要求1所述的方法，其特征在于，所述步骤二中的关键点的图案为黑白同心圆环。

6.根据权利要求1所述的方法，其特征在于，所述步骤三中进行步骤3.1至3.7的处理，获得关键点集的屏幕坐标

$p_{left}^{\left(i\right)}={\left[\begin{array}{cc}{x}_{left}^{\left(i\right)}& y_{left}^{\left(i\right)}\end{array}\right]}^T,p_{right}^{\left(i\right)}={\left[\begin{array}{cc}{x}_{right}^{\left(i\right)}& y_{right}^{\left(i\right)}\end{array}\right]}^T$

上标i表示第i个关键点，两个分量为关键点在校正图像中的x坐标、y坐标；

所述步骤三进一步包括：

步骤3.1，对两张图片使用图像行对齐映射处理得到消除畸变的行对齐图像；

步骤3.2，对两张图片使用高斯模糊去除图像中的噪点；

步骤3.3，对两张图片使用二值阀值化，将图像转化为黑白图像；

步骤3.4，求取二值阀值化后的图像的轮廓；

步骤3.5，分别对两张图片的所有轮廓执行剔除操作；

(1)周长条件：轮廓的周长∈[L_min，L_max]

(2)面积条件：轮廓的实际面积∈[S_min，S_max]

(3)面积周长比条件：轮廓的实际面积比周长＜ratio

(4)拟合误差条件：(轮廓的实际面积-轮廓拟合成椭圆的面积)/轮廓拟合面积∈[-e_e，e_e]

(5)当前轮廓数量如果小于2倍于需要的关键点数，则暂停计算，提示用户更改参数，重新进行剔除；

(6)匹配轮廓对，当轮廓的圆心相距不超过d_D时，匹配为同一对轮廓；

(7)计算匹配的轮廓对数量，如果轮廓对数量不等于关键点数量，则暂停计算，提示用户更改参数，重新进行剔除；

(8)对匹配出来的轮廓对求外圆的圆心坐标；

(9)与上一帧的位置匹配，保存结果；

以上参数L_min，L_max，S_min，S_max，ratio，e_e，d_D由用户根据实验情况设定；

步骤3.6，对左右两张图片的轮廓对匹配，计算得到所有关键点集的屏幕坐标；

步骤3.7，针对关键点集的轮廓，使用Zernike矩对轮廓边界进行亚像素定位，之后重新进行椭圆拟合，计算拟合椭圆的中心，更新相应的屏幕坐标；

其中，使用Zernike模板矩计算亚像素边界的步骤为：

选取合适的模板圆半径N，计算三个不同阶次的矩，分别为：

$A_{00}=\∫{\∫}_{x^2+y^2\≤1}f\left(x,y\right)dxdy$

$A_{11}=\∫{\∫}_{x^2+y^2\≤1}f\left(x,y\right)\left(x-yj\right)dxdy$

$A_{20}=\∫{\∫}_{x^2+y^2\≤1}f\left(x,y\right)\left(x^2+y^2-1\right)dxdy$

其中x，y是以模板圆半径N归一化的像素坐标，f(x，y)为像素点(x，y)的灰度值；

$\mathrm{\φ}={\tan }^{-1}\left(\frac{\mathrm{Im}\left(A_{11}\right)}{\mathrm{Re}\left(A_{11}\right)}\right)$

$l=\frac{A_{20}}{A_{11}^{\′}}$

原轮廓上点(x₀，y₀)^T的亚像素边界点为：

$\left[\begin{array}{c}{x}_{sub}\\ y_{sub}\end{array}\right]=\left[\begin{array}{c}{x}_0\\ y_0\end{array}\right]+Nl\left[\begin{array}{c}\cos \mathrm{\φ}\\ \sin \mathrm{\φ}\end{array}\right]$

步骤3.8，根据 $Q\left[\begin{array}{c}{x}_{left}^{\left(i\right)}\\ y_{left}^{\left(i\right)}\\ x_{right}^{\left(i\right)}-x_{left}^{\left(i\right)}\\ 1\end{array}\right]=\left[\begin{array}{c}{X}^{\left(i\right)}\\ Y^{\left(i\right)}\\ Z^{\left(i\right)}\\ W^{\left(i\right)}\end{array}\right],$ 计算得到所有关键点集的三维坐标： $\left[\begin{array}{ccc}\frac{X^{\left(i\right)}}{W^{\left(i\right)}}& \frac{Y^{\left(i\right)}}{W^{\left(i\right)}}& \frac{Z^{\left(i\right)}}{W^{\left(i\right)}}\end{array}\right],$ 其中，Q为重投影矩阵，

$Q=\left[\begin{array}{cccc}1& 0& 0& -c_x\\ 0& 1& 0& -c_y\\ 0& 0& 0& f\\ 0& 0& -1/t_x& 0\end{array}\right]$

其中t_x为右相机单元相对左相机单元的平移向量x值，其它参数均取自左相机，焦距f为(f_x+f_y)/2；

步骤3.9，输出测试时间区间内的所有关键点的三维运动数据。

7.根据权利要求1所述的方法，其特征在于，在所述步骤四中，

对于n维线性时不变系统，离散时间状态空间方程为：

x(k+1)＝Ax(k)+Bu(k)

y(k)＝Cx(k)

其中x为t×1阶状态向量，u为m×1阶控制向量，y为p×1阶观测向量，A为t×t阶系统矩阵，B为t×m阶控制矩阵，C为p×t阶观测矩阵。对于脉冲响应：

Y(k)＝CA^k-1B，k＝1，2，3，...

Hanke1矩阵定义为pr×ms阶矩阵，

步骤4.1，由双目视觉系统测量得到的运动数据，给定块矩阵行数r和列数s，分别构造零阶汉克尔矩阵H_rs(0)和一阶汉克尔矩阵H_rs(1)；

步骤4.2，对H_rs(0)进行奇异值分解(SVD)，H_rs(0)＝PDQ^T，其中P为pr×pr阶酉矩阵，D为半正定pr×ms阶对角矩阵，且对角线元素按从大到小排列，Q^T为ms×ms阶酉矩阵；

步骤4.3，根据给定阀值ε得到系统的阶次n，则

$H_{rs}\left(0\right)={\mathrm{PDQ}}^T\≈P\left[\begin{array}{cc}{D}_n& 0\\ 0& 0\end{array}\right]Q^T,D_n=diag\left(\begin{array}{cccc}{d}_1& d_2& ...& d_n\end{array}\right)$

其中D_n为大于阀值ε的n个对角线元素组成的方阵；

步骤4.4，根据 $A=D_n^{-1/2}{P}_n^T{H}_{rs}\left(1\right)Q_n{D}_n^{-1/2},B=D_n^{1/2}{Q_n}^T{E}_m,C=E_p^T{P}_n{D}_n^{1/2}$ 得到系统最小实现，其中P_n和Q_n分别为P和Q的前n列，

$E_m^T=\left[\begin{array}{cc}{I}_m& 0_{m\×(ms-m)}\end{array}\right],E_p^T=\left[\begin{array}{cc}{I}_p& 0_{p\×(pr-p)}\end{array}\right]$

步骤4.5，对矩阵A进行特征值分解，A＝VDV^-1，其中：D为特征值对角阵，D＝diag(z₁z₂…z_n)；V为特征向量矩阵，

步骤4.6，计算各阶模态参数：无阻尼频率ω_i、阻尼比振型、模态置信度(MAC)、模态相位共线性(MPC)；

${\mathrm{\λ}}_i=\frac{\ln \left(z_i\right)}{T},i=1,2,...,n$

${\mathrm{\ω}}_i=\sqrt{\mathrm{Re}{\left({\mathrm{\λ}}_i\right)}^2+\mathrm{Im}{\left({\mathrm{\λ}}_i\right)}^2}$

${\mathrm{MAC}}_i=\frac{|q_i^{*}\·{\stackrel{\‾}{q}}_i^{*}|}{\sqrt{|{\stackrel{\‾}{q}}_i^{*}\·{\stackrel{\‾}{q}}_i^{*}||q_i^{*}\·q_i^{*}|}}$

对c_i进行归一化，使||c_i||₂＝1，设θ₀是振幅最大的点对应的相位，对c_i的每个分量旋转角度-θ₀，

${c_i}^{*}={\left(\begin{array}{cccc}{c}_{i1}{e}^{-{\mathrm{\θ}}_{0}i}& c_{i2}{e}^{-{\mathrm{\θ}}_{0}i}& ...& c_{ip}{e}^{-{\mathrm{\θ}}_{0}i}\end{array}\right)}^T$

${\mathrm{MPC}}_i=\frac{\left|\right|c_i^{*}|{|}_2}{\left|\right|c_i|{|}_2}=\left|\right|c_i^{*}|{|}_2$

二范数归一化振型，

步骤4.7，输出并保存结果。