一种测量风力机叶片相对变形的结构故障检测方法及系统

摘要

本发明公开了一种测量风力机叶片相对变形的结构故障检测方法及系统，该方法包括：计算变形区域相对于参考区域的旋转角度；生成变形区域参数；更新旋转角度后变形区域像素点坐标；计算参考区域与更新后的变形区域的像素点匹配相关系数；更新变形区域参数；更新像素点坐标；计算参考区域与二次更新后的变形区域的像素点匹配相关系数；比较所述二次更新后的变形区域的像素点匹配相关系数与预设阈值的大小；计算变形区域像素点的三维坐标和叶片的应变，判断和定位叶片的结构故障。采用本发明的测量风力机叶片相对变形的结构故障检测方法及系统，能够克服传统数字图像相关技术方法在大旋转角的匹配困难这一技术难题，适用于旋转叶片的测量。

说明书

技术领域

本发明涉及光电检测技术领域，特别是涉及一种测量风力机叶片相对变形的结构故障检测方法及系统。

背景技术

风力机叶片的各种故障在受力状态变化时有明显的外在表现，例如，表面裂纹张开、内部缺陷使表面鼓包等。风电叶片的许多故障只在叶片转动中能表现出来。目前，根据叶片的变形参数测量结果对叶片健康检测评估方法是基于传统的DIC方法。

数字图像相关技术(Digital Image Correlation，DIC)是一种基于现代数字图像处理和分析技术的新型先进光测技术，它通过分析变形前后物体表面的数字图像来获得被测物体表面的变形(位移和应变)信息，具有非接触、对环境和隔振要求低、测量范围广以及自动化程度高等优点。但是，传统的DIC方法对于旋转角大于7°的变形物体的图像匹配会失效，因此难以运用到旋转叶片的测量。

发明内容

本发明的目的是提供一种适用于任意旋转角度测量的测量风力机叶片相对变形的结构故障检测方法及系统。

为实现上述目的，本发明提供了如下方案：

一种测量风力机叶片相对变形的结构故障检测方法，包括：

获取参考区域与变形区域；所述参考区域为在变形前的叶片表面图像上选取的一块区域，所述参考区域包括若干个参考区域特征点；所述变形区域为对变形后的叶片表面图像进行特征点匹配，得到与所述参考区域匹配的变形区域，所述变形区域包括若干个与所述参考区域特征点对应的变形区域特征点；所述参考区域与所述变形区域横像素点个数与纵像素点个数均相同；

根据所述参考区域特征点与所述变形区域特征点，计算所述变形区域相对于所述参考区域的旋转角度θ；

将所述变形区域旋转θ角度后，将所述参考区像素点坐标(x_i,y_i)与对应的所述旋转θ角度后变形区域像素点坐标(x_i^*,y_i^*)做差，得到所述旋转θ角度后变形区域像素点相对于所述参考区域像素点的x方向的位移u和y方向的位移v；由所述旋转角度θ以及所述旋转θ角度后变形区域像素点相对于所述参考区域像素点的x方向的位移u和y方向的位移v生成变形区域参数

选取所述旋转θ角度后变形区域中任一像素点作为变形区基准点，根据所述变形区域参数以及所述旋转θ角度后变形区域中每一像素点到所述变形区基准点的距离(Δx^*,Δy^*)更新所述旋转θ角度后变形区域像素点坐标(x_i^*,y_i^*)，得到更新后的变形区域像素点坐标(x_i^*’,y_i^*’)；

根据更新后的变形区域像素点坐标(x_i^*’,y_i^*’)处的像素点的灰度值和所述参考区域像素点坐标(x_i,y_i)处的像素点的灰度值计算所述参考区域与所述更新后的变形区域的像素点匹配相关系数C(p)；

根据所述像素点匹配相关系数C(p)和所述变形区域参数更新所述变形区域参数得到更新后的变形区域参数

根据所述更新后的变形区域参数以及所述旋转θ’角度后变形区域中每一像素点到所述变形区基准点的距离(Δx^*’,Δy^*’)二次更新所述更新后的变形区域像素点坐标(x_i^*’,y_i^*’)，得到二次更新后的变形区域像素点坐标(x_i^*”,y_i^*”)；

根据所述二次更新后的变形区域像素点坐标(x_i^*”,y_i^*”)处的像素点的灰度值和所述参考区域像素点坐标(x_i,y_i)处的像素点的灰度值计算所述参考区域与所述二次更新后的变形区域的像素点匹配相关系数C’(p)；

比较所述二次更新后的变形区域的像素点匹配相关系数C’(p)与预设阈值的大小；若所述二次更新后的变形区域的像素点匹配相关系数C’(p)小于所述阈值，将所述更新后的变形区域的像素点匹配相关系数C(p)更新为二次更新后的变形区域的像素点匹配相关系数C’(p)，并返回步骤“根据所述像素点匹配相关系数C(p)和所述变形区域参数更新所述变形区域参数得到更新后的变形区域参数否则，判断是否将所述旋转θ’角度后变形区域中所有像素点作为变形区基准点，若没有将所述旋转θ’角度后变形区域中所有像素点作为变形区基准点，更新所述变形区基准点，计算所述旋转θ’角度后变形区域中每一像素点到所述变形区更新后的基准点的距离(Δx^*’,Δy^*’)，根据所述更新后的变形区域参数以及所述旋转θ’角度后变形区域中每一像素点到所述变形区更新后的基准点的距离(Δx^*’,Δy^*’)二次更新所述变形区域参数中的u’,v’参数，得到二次更新后的变形区域参数将所述变形区域参数更新为二次更新后的变形区域参数并返回步骤“根据所述变形区域参数更新所述旋转θ角度后变形区域像素点坐标(x_i^*,y_i^*)，得到更新后的变形区域像素点坐标(x_i^*’,y_i^*’)”，若将所述旋转θ’角度后变形区域中所有像素点作为变形区基准点，则根据所述二次更新后的变形区域像素点坐标(x_i^*”,y_i^*”)计算所述旋转θ’角度后变形区域像素点的三维坐标；

根据所述旋转θ’角度后变形区域像素点的三维坐标和与所述旋转θ’角度后变形区域像素点对应的所述参考区域像素点的三维坐标计算所述变形后的叶片的应变，并根据所述变形后的叶片的应变判断和定位叶片的结构故障。

可选的，所述根据所述参考区域特征点与所述变形区域特征点，计算所述变形区域相对于所述参考区域的旋转角度θ，具体包括：

获取所述参考区域中两个参考区特征点坐标K(x_i,y_i)和K(x_i+1,y_i+1)，其中，i＝1,…,n，n≥2；

在所述变形区域中利用尺度不变特征变化算法匹配得到与所述参考区域中两个参考区特征点对应的两个变形区特征点坐标K^*(x_i^*,y_i^*)和K^*(x_i+1^*,y_i+1^*)，其中，i＝1,…,n，n≥2；

根据如下公式计算所述变形区域相对于所述参考区域的旋转角度θ：

可选的，根据如下公式计算得到所述旋转θ角度后变形区域像素点坐标(x_i^*,y_i^*)：

其中，(x_i,y_i)为参考区像素点坐标，(Δx,Δy)为参考区域中每一像素点到参考区基准点的距离，所述参考区基准点为参考区域中任一像素点。

可选的，所述根据所述变形区域参数以及所述旋转θ角度后变形区域中每一像素点到所述变形区基准点的距离(Δx^*,Δy^*)更新所述旋转θ角度后变形区域像素点坐标(x_i^*,y_i^*)，得到更新后的变形区域像素点坐标(x_i^*’,y_i^*’)，具体包括：

根据如下公式更新所述旋转θ角度后变形区域像素点坐标(x_i^*,y_i^*)：

x_i^*'＝x₀+u+Δx^*cosθ+Δy^*sinθ

y_i^*'＝y₀+v-Δx^*sinθ+Δy^*cosθ

其中，(x₀,y₀)为参考区基准点的坐标，所述参考区基准点为参考区域中任一像素点。

可选的，所述根据更新后的变形区域像素点坐标(x_i^*’,y_i^*’)处的像素点的灰度值和所述参考区域像素点坐标(x_i,y_i)处的像素点的灰度值计算所述参考区域与所述更新后的变形区域的像素点匹配相关系数C(p)，具体包括：

根据如下公式计算所述参考区域与所述更新后的变形区域的像素点匹配相关系数C(p)：

其中，M为所述参考区域横像素中间像素点所在序号，所述参考区域横像素个数为奇数，f(x_i,y_j)为所述参考区域像素点的灰度值，为所述参考区域像素点的平均灰度值，为所述所述旋转θ角度后变形区域像素点的灰度值，为所述旋转θ角度后变形区域像素点的平均灰度值。

可选的，所述根据所述像素点匹配相关系数C(p)和所述变形区域参数更新所述变形区域参数得到更新后的变形区域参数具体包括：

根据如下公式更新所述变形区域参数

其中，▽C(p)为C(p)的雅克比矩阵，▽▽C(p)为C(p)的海森矩阵。

可选的，根据所述更新后的变形区域参数以及所述旋转θ’角度后变形区域中每一像素点到所述变形区基准点的距离(Δx^*’,Δy^*’)二次更新所述更新后的变形区域像素点坐标(x_i^*’,y_i^*’)，得到二次更新后的变形区域像素点坐标(x_i^*”,y_i^*”)，具体包括：

根据如下公式更新所述更新后的变形区域像素点坐标(x_i^*’,y_i^*’)：

其中，(x₀,y₀)为参考区基准点的坐标，u_x，u_y为u的一阶导数，v_x，v_y为v的一阶导数，u_xx，u_xy，u_yy为u的二阶导数，v_xx，v_xy，v_yy为v的二阶导数。

可选的，根据所述二次更新后的变形区域像素点坐标(x_i^*”,y_i^*”)处的像素点的灰度值和所述参考区域像素点坐标(x_i,y_i)处的像素点的灰度值计算所述参考区域与所述二次更新后的变形区域的像素点匹配相关系数C’(p)，具体包括：

根据如下公式计算所述参考区域与所述二次更新后的变形区域的像素点匹配相关系数C’(p)：

其中，M为所述参考区域横像素中间像素点所在序号，所述参考区域横像素个数为奇数，f(x_i,y_j)为所述参考区域像素点的灰度值，为所述参考区域像素点的平均灰度值，为所述所述旋转θ角度后变形区域像素点的灰度值，为所述旋转θ角度后变形区域像素点的平均灰度值。

可选的，所述根据所述更新后的变形区域参数以及所述旋转θ’角度后变形区域中每一像素点到所述变形区更新后的基准点的距离(Δx^*’,Δy^*’)二次更新所述变形区域参数中的u’,v’参数，得到二次更新后的变形区域参数具体包括：

根据如下公式二次更新所述变形区域参数中的u’,v’参数：

u”＝u'+x₀-x₁+Δx^*'cosθ'+Δy^*'sinθ'

v₁”＝v'+y₀-y₁-Δx^*'sinθ'+Δy^*'cosθ'

其中，(x₀,y₀)为参考区基准点的坐标，所述参考区基准点为参考区域中任一像素点，(x₁,y₁)为参考区更新后的基准点的坐标，所述参考区更新后的基准点与所述变形区更新后的基准点位置对应。

可选的，所述根据所述旋转θ’角度后变形区域像素点的三维坐标和与所述旋转θ’角度后变形区域像素点对应的所述参考区域像素点的三维坐标计算所述变形后的叶片的应变，具体包括：利用坐标网格的方法计算所述变形后的叶片的应变。

本发明还提供一种测量风力机叶片相对变形的结构故障检测系统，包括：

参考区域与变形区域获取模块，用于获取参考区域与变形区域；所述参考区域为在变形前的叶片表面图像上选取的一块区域，所述参考区域包括若干个参考区域特征点；所述变形区域为对变形后的叶片表面图像进行特征点匹配，得到与所述参考区域匹配的变形区域，所述变形区域包括若干个与所述参考区域特征点对应的变形区域特征点；所述参考区域与所述变形区域横像素点个数与纵像素点个数均相同；

参考区域的旋转角度θ计算模块，用于根据所述参考区域特征点与所述变形区域特征点，计算所述变形区域相对于所述参考区域的旋转角度θ；

变形区域参数生成模块，用于将所述变形区域旋转θ角度后，将所述参考区像素点坐标(x_i,y_i)与对应的所述旋转θ角度后变形区域像素点坐标(x_i^*,y_i^*)做差，得到所述旋转θ角度后变形区域像素点相对于所述参考区域像素点的x方向的位移u和y方向的位移v；由所述旋转角度θ以及所述旋转θ角度后变形区域像素点相对于所述参考区域像素点的x方向的位移u和y方向的位移v生成变形区域参数

旋转θ角度后变形区域像素点坐标更新模块，用于选取所述旋转θ角度后变形区域中任一像素点作为变形区基准点，根据所述变形区域参数以及所述旋转θ角度后变形区域中每一像素点到所述变形区基准点的距离(Δx^*,Δy^*)更新所述旋转θ角度后变形区域像素点坐标(x_i^*,y_i^*)，得到更新后的变形区域像素点坐标(x_i^*’,y_i^*’)；

参考区域与更新后的变形区域的像素点匹配相关系数计算模块，用于根据更新后的变形区域像素点坐标(x_i^*’,y_i^*’)处的像素点的灰度值和所述参考区域像素点坐标(x_i,y_i)处的像素点的灰度值计算所述参考区域与所述更新后的变形区域的像素点匹配相关系数C(p)；

变形区域参数更新模块，用于根据所述像素点匹配相关系数C(p)和所述变形区域参数更新所述变形区域参数得到更新后的变形区域参数

变形区域像素点坐标二次更新模块，用于根据所述更新后的变形区域参数以及所述旋转θ’角度后变形区域中每一像素点到所述变形区基准点的距离(Δx^*’,Δy^*’)二次更新所述更新后的变形区域像素点坐标(x_i^*’,y_i^*’)，得到二次更新后的变形区域像素点坐标(x_i^*”,y_i^*”)；

参考区域与二次更新后的变形区域的像素点匹配相关系数计算模块，用于根据所述二次更新后的变形区域像素点坐标(x_i^*”,y_i^*”)处的像素点的灰度值和所述参考区域像素点坐标(x_i,y_i)处的像素点的灰度值计算所述参考区域与所述二次更新后的变形区域的像素点匹配相关系数C’(p)；

第一判断模块，用于比较所述二次更新后的变形区域的像素点匹配相关系数C’(p)与预设阈值的大小；若所述二次更新后的变形区域的像素点匹配相关系数C’(p)小于所述阈值，将所述更新后的变形区域的像素点匹配相关系数C(p)更新为二次更新后的变形区域的像素点匹配相关系数C’(p)，并将指令发送至变形区域参数更新模块；否则，将指令发送至第二判断模块；

第二判断模块，用于判断是否将所述旋转θ’角度后变形区域中所有像素点作为变形区基准点，若没有将所述旋转θ’角度后变形区域中所有像素点作为变形区基准点，则将指令发送至变形区域参数二次更新模块，否则将指令发送至三维坐标计算模块；

所述变形区域参数二次更新模块，用于更新所述变形区基准点，计算所述旋转θ’角度后变形区域中每一像素点到所述变形区更新后的基准点的距离(Δx^*’,Δy^*’)，根据所述更新后的变形区域参数以及所述旋转θ’角度后变形区域中每一像素点到所述变形区更新后的基准点的距离(Δx^*’,Δy^*’)二次更新所述变形区域参数中的u’,v’参数，得到二次更新后的变形区域参数将所述变形区域参数更新为二次更新后的变形区域参数并将指令发送至旋转θ角度后变形区域像素点坐标更新模块；

所述三维坐标计算模块，用于根据所述二次更新后的变形区域像素点坐标(x_i^*”,y_i^*”)计算所述旋转θ’角度后变形区域像素点的三维坐标；

变形后的叶片的应变计算模块，用于根据所述旋转θ’角度后变形区域像素点的三维坐标和与所述旋转θ’角度后变形区域像素点对应的所述参考区域像素点的三维坐标计算所述变形后的叶片的应变，并根据所述变形后的叶片的应变判断和定位叶片的结构故障。

与现有技术相比，本发明的有益效果是：

本发明提出了一种测量风力机叶片相对变形的结构故障检测方法及系统，考虑了变形区域相对于所述参考区域的旋转角度θ，克服了传统DIC方法在大旋转角的匹配困难这一技术难题并获得了叶片在旋转状态下的相对应变，适用于旋转叶片的测量。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。

图1为本发明实施例一中测量风力机叶片相对变形的结构故障检测方法流程图；

图2为本发明实施例一中应变计算中虚拟德罗内三角网格构建示意图；

图3为本发明实施例一中应变计算中空间三角形网格坐标变换示意图；

图4为本发明实施例二中测量风力机叶片相对变形的结构故障检测系统结构图；

图5为本发明实施例三中含缺陷叶片沿径向方向的应变变化曲线；

图6为本发明实施例四中5kw小型风电叶片其它位置相对水平位置在径向的应变变化曲线。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

本发明的目的是提供一种适用于任意旋转角度测量的测量风力机叶片相对变形的结构故障检测方法及系统。

为使本发明的上述目的、特征和优点能够更加明显易懂，下面结合附图和具体实施方式对本发明作进一步详细的说明。

实施例一

图1为本发明实施例中基于测量风力机叶片相对变形的结构故障检测方法流程图，如图1所示，本发明提供的测量风力机叶片相对变形的结构故障检测方法，包括：

步骤101：获取参考区域与变形区域。

参考区域为在变形前的叶片表面图像上选取的一块区域，所述参考区域包括若干个参考区域特征点；所述变形区域为对变形后的叶片表面图像进行特征点匹配，利用尺度不变特征变换算法得到与所述参考区域匹配的变形区域，所述变形区域包括若干个与所述参考区域特征点对应的变形区域特征点；所述参考区域与所述变形区域横像素点个数与纵像素点个数均相同。

步骤102：计算变形区域相对于参考区域的旋转角度θ。

获取所述参考区域中两个参考区特征点坐标K(x_i,y_i)和K(x_i+1,y_i+1)，其中，i＝1,…,n，n≥2；

在所述变形区域中利用尺度不变特征变化算法匹配得到与所述参考区域中两个参考区特征点对应的两个变形区特征点坐标K^*(x_i^*,y_i^*)和K^*(x_i+1^*,y_i+1^*)，其中，i＝1,…,n，n≥2；

根据如下公式计算所述变形区域相对于所述参考区域的旋转角度θ：

θ表示(x_i+1-x_i,y_i+1-y_i)和(x_i+1^*-x_i^*,y_i+1^*-y_i^*)的向量夹角，其范围为0°到180°，而旋转角度θ是表示参考区域逆时针旋转到变形区域，其范围为0°到360°。为了正确的计算这个旋转角度θ，向量(x_i+1-x_i,y_i+1-y_i)和(x_i+1^*-x_i^*,y_i+1^*-y_i^*)的旋转方向需要被确定。为此，引进了向量(x_i+1-x_i,y_i+1-y_i)到(x_i+1^*-x_i^*,y_i+1^*-y_i^*)的向量积,其表达式如下:

假如是正值，则(x_i+1^*-x_i^*,y_i+1^*-y_i^*)是关于(x_i+1-x_i,y_i+1-y_i)逆时针旋转，于是θ＝θ；否则，θ＝360°-θ。对于一系列向量(x_i+1-x_i,y_i+1-y_i)和(x_i+1^*-x_i^*,y_i+1^*-y_i^*)，其最终的旋转角θ可以通过求θ的平均值获得。

步骤103：生成变形区域参数

首先根据如下公式计算得到所述旋转θ角度后变形区域像素点坐标(x_i^*,y_i^*)：

其中，(x_i,y_i)为参考区像素点坐标，(Δx,Δy)为参考区域中每一像素点到参考区基准点的距离，所述参考区基准点为参考区域中任一像素点。

然后将所述变形区域旋转θ角度后，将所述参考区像素点坐标(x_i,y_i)与对应的所述旋转θ角度后变形区域像素点坐标(x_i^*,y_i^*)做差，得到所述旋转θ角度后变形区域像素点相对于所述参考区域像素点的x方向的位移u和y方向的位移v；即u＝x_i-x_i^*，v＝y_i-y_i^*。

由所述旋转角度θ以及所述旋转θ角度后变形区域像素点相对于所述参考区域像素点的x方向的位移u和y方向的位移v生成变形区域参数

步骤104：更新旋转θ角度后变形区域像素点坐标(x_i^*,y_i^*)。

选取所述旋转θ角度后变形区域中任一像素点作为变形区基准点，根据所述变形区域参数以及所述旋转θ角度后变形区域中每一像素点到所述变形区基准点的距离(Δx^*,Δy^*)更新所述旋转θ角度后变形区域像素点坐标(x_i^*,y_i^*)，得到更新后的变形区域像素点坐标(x_i^*’,y_i^*’)。

更新所述旋转θ角度后变形区域像素点坐标(x_i^*,y_i^*)的具体公式为：

x_i^*'＝x₀+u+Δx^*cosθ+Δy^*sinθ

y_i^*'＝y₀+v-Δx^*sinθ+Δy^*cosθ

其中，(x₀,y₀)为参考区基准点的坐标，所述参考区基准点为参考区域中任一像素点。

步骤105：计算参考区域与更新后的变形区域的像素点匹配相关系数C(p)。

根据更新后的变形区域像素点坐标(x_i^*’,y_i^*’)处的像素点的灰度值和所述参考区域像素点坐标(x_i,y_i)处的像素点的灰度值计算所述参考区域与所述更新后的变形区域的像素点匹配相关系数C(p)，计算公式如下：

其中，M为所述参考区域横像素中间像素点所在序号，所述参考区域横像素个数为奇数，2M+1为所述参考区域横像素点个数，f(x_i,y_j)为所述参考区域像素点的灰度值，为所述参考区域像素点的平均灰度值，为所述所述旋转θ角度后变形区域像素点的灰度值，为所述旋转θ角度后变形区域像素点的平均灰度值。

步骤106：更新变形区域参数

根据所述像素点匹配相关系数C(p)和所述变形区域参数更新所述变形区域参数得到更新后的变形区域参数考虑旋转角的牛顿拉夫逊迭代算法，即NR迭代算法更新变形区域参数的具体公式为：

其中，▽C(p)为C(p)的雅克比矩阵，即Jacobian矩阵，▽▽C(p)为C(p)的海森矩阵，即Hessian矩阵。

在此迭代中，利用高精度的双三次插值算法：

其中，M为所述参考区域横像素中间像素点所在序号，所述参考区域横像素个数为奇数，2M+1为所述参考区域横像素点个数，f(x,y)为所述参考区域像素点的灰度值，为所述参考区域像素点的平均灰度值，g(x^**,y^**)为所述所述旋转θ角度后变形区域像素点的亮度值，为所述旋转θ角度后变形区域像素点的亮度值，a_mn为每个像素通过双三次插值算法拟合得到的系数值。

步骤107：更新像素点坐标(x_i^*’,y_i^*’)。

根据所述更新后的变形区域参数以及所述旋转θ’角度后变形区域中每一像素点到所述变形区基准点的距离(Δx^*’,Δy^*’)二次更新所述更新后的变形区域像素点坐标(x_i^*’,y_i^*’)，得到二次更新后的变形区域像素点坐标(x_i^*”,y_i^*”)，二次更新所述更新后的变形区域像素点坐标(x_i^*’,y_i^*’)的具体公式为：

其中，(x₀,y₀)为参考区基准点的坐标，u_x，u_y为u的一阶导数，v_x，v_y为v的一阶导数，u_xx，u_xy，u_yy为u的二阶导数，v_xx，v_xy，v_yy为v的二阶导数。该等式与传统的二阶形函数不同是其包含了旋转角θ，可以看到在二阶形函数中总共有13个变量，即(u’,v’,u_x,u_y,v_x,v_y,u_xx,u_xy,u_yy,v_xx,v_xy,v_yy,θ’)。则NR迭代的初始值为(u’,v’,0,0,0,0,0,0,0,0,0,0,θ’)，利用该形函数，结合传统的NR迭代算法则可消除由大旋转与变形引起的去相关，从而能精确的获得变形物体的参数(u’,v’,u_x,u_y,v_x,v_y,u_xx,u_xy,u_yy,v_xx,v_xy,v_yy)。

步骤108：计算参考区域与二次更新后的变形区域的像素点匹配相关系数C’(p)。

根据所述二次更新后的变形区域像素点坐标(x_i^*”,y_i^*”)处的像素点的灰度值和所述参考区域像素点坐标(x_i,y_i)处的像素点的灰度值计算所述参考区域与所述二次更新后的变形区域的像素点匹配相关系数C’(p)，计算所述参考区域与所述二次更新后的变形区域的像素点匹配相关系数C’(p)的具体公式为：

其中，2M+1为所述参考区域横像素点个数，f(x_i,y_j)为所述参考区域像素点的灰度值，为所述参考区域像素点的平均灰度值，为所述所述旋转θ角度后变形区域像素点的灰度值，为所述旋转θ角度后变形区域像素点的平均灰度值。

步骤109：比较所述二次更新后的变形区域的像素点匹配相关系数C’(p)与预设阈值的大小。

比较所述二次更新后的变形区域的像素点匹配相关系数C’(p)与预设阈值的大小.

若所述二次更新后的变形区域的像素点匹配相关系数C’(p)小于所述阈值，将所述更新后的变形区域的像素点匹配相关系数C(p)更新为二次更新后的变形区域的像素点匹配相关系数C’(p)，并返回步骤106“根据所述像素点匹配相关系数C(p)和所述变形区域参数更新所述变形区域参数得到更新后的变形区域参数

若所述二次更新后的变形区域的像素点匹配相关系数C’(p)大于或等于所述阈值，执行步骤110，判断是否将所述旋转θ’角度后变形区域中所有像素点作为变形区基准点。

步骤110：判断是否将所述旋转θ’角度后变形区域中所有像素点作为变形区基准点。

若没有将所述旋转θ’角度后变形区域中所有像素点作为变形区基准点，更新所述变形区基准点，执行步骤111；若将所述旋转θ’角度后变形区域中所有像素点作为变形区基准点，则执行步骤112“根据所述二次更新后的变形区域像素点坐标(x_i^*”,y_i^*”)计算所述旋转θ’角度后变形区域像素点的三维坐标”。

步骤111：二次更新后的变形区域参数

计算所述旋转θ’角度后变形区域中每一像素点到所述变形区更新后的基准点的距离(Δx^*’,Δy^*’)，根据所述更新后的变形区域参数以及所述旋转θ’角度后变形区域中每一像素点到所述变形区更新后的基准点的距离(Δx^*’,Δy^*’)二次更新所述变形区域参数中的u’,v’参数，得到二次更新后的变形区域参数将所述变形区域参数更新为二次更新后的变形区域参数并返回步骤104“根据所述变形区域参数更新所述旋转θ角度后变形区域像素点坐标(x_i^*,y_i^*)，得到更新后的变形区域像素点坐标(x_i^*’,y_i^*’)”。其中，根据如下公式二次更新所述变形区域参数中的u’,v’参数：

u”＝u'+x₀-x₁+Δx^*'cosθ'+Δy^*'sinθ'

v₁”＝v'+y₀-y₁-Δx^*'sinθ'+Δy^*'cosθ'

其中，(x₀,y₀)为参考区基准点的坐标，所述参考区基准点为参考区域中任一像素点，(x₁,y₁)为参考区更新后的基准点的坐标，所述参考区更新后的基准点与所述变形区更新后的基准点位置对应。

步骤112：计算旋转θ’角度后变形区域像素点的三维坐标。

利用两个相机图像相关匹配得到的像素坐标与双相机系统的标定参数即可计算出叶片变形前后的三维坐标，二次更新后的变形区域像素点坐标为(x_i^*”,y_i^*”)，假设某个像素点的三维世界坐标为(X_w,Y_w,Z_w)则可通过以下公式计算得到：

其中，(x_i1^*”,y_i1^*”)和(x_i2^*”,y_i2^*”)为分别为相机1和2的二次更新后的变形区域像素点坐标，表示投影矩阵，它由相机的内外参决定。在上式中，世界坐标(X_w,Y_w,Z_w)能被这四个等式利用最小二乘法精确求解。

步骤113：计算变形后的叶片的应变，判断和定位叶片的结构故障。

图2为本发明实施例应变计算中虚拟德罗内Delaunay三角网格构建示意图，如图2所示，利用坐标网格的方法计算所述变形后的叶片的应变。该方法首先在参考状态中待测点P与其附近的4个点构建4个虚拟德罗内Delaunay三角网格。

图3为应变计算中空间三角形网格坐标变换示意图，对于其中一个三角网格，如图3(a)所示，假设该三角网格由P(x₀,y₀,z₀)，P1(x₁,y₁,z₁)和P2(x₂,y₂,z₂)组成。在计算之前需要将三角形网格的三维空间坐标进行坐标转换，将其转换到如图3(b)所示的平面正交直角坐标系中，其转换后的坐标公式如下：

其中，v₁，v₂,v’₁，v’₂分别为三角形各边长所对应的向量。

根据如下公式计算上式中坐标定义面内变形梯度张量F：

F₁₁X₁+F₁₂Y₁＝X′₁,F₁₁X₂+F₁₂Y₂＝X′₂

F₂₁X₁+F₂₂Y₁＝Y′₁,F₂₁X₂+F₂₂Y₂＝Y′₂

根据有限变形理论，面内梯度张量F和柯西-格林张量C的具有关系：

C＝F^TF

进一步根据柯西-格林应变张量C与拉格朗日应变张量E的关系C＝I+2E可求出E,其中I为单位矩阵。定义：

真实应变可由拉格朗日应变表示：

其中ε₁，ε₂分别为最大和最小主应变。根据如下公式计算其等效应变：

同理计算4个三角形网格所对应的应变，并取其平均值作为P点最终的应变。选取其中任意状态作为参考状态，并求出其他所有状态相对该状态的下的应变分布，提取其应变云图，并将应变云图覆盖在叶片的参考图中。从而可以通过应变云图判断集中应变，进一步提取叶片沿着径向方向的应变变化，结合云图判断应变集中现象。假如某处在旋转过程中始终存在应变集中，则可判断该处存在裂纹或者鼓包等缺陷。

实施例二

图4为本发明实施例中测量风力机叶片相对变形的结构故障检测系统结构图，如图4所示，本发明提供的测量风力机叶片相对变形的结构故障检测系统，包括：

参考区域与变形区域获取模块201，用于获取参考区域与变形区域；所述参考区域为在变形前的叶片表面图像上选取的一块区域，所述参考区域包括若干个参考区域特征点；所述变形区域为对变形后的叶片表面图像进行特征点匹配，得到与所述参考区域匹配的变形区域，所述变形区域包括若干个与所述参考区域特征点对应的变形区域特征点；所述参考区域与所述变形区域横像素点个数与纵像素点个数均相同；

参考区域的旋转角度θ计算模块202，用于根据所述参考区域特征点与所述变形区域特征点，计算所述变形区域相对于所述参考区域的旋转角度θ；

变形区域参数生成模块203，用于将所述变形区域旋转θ角度后，将所述参考区像素点坐标(x_i,y_i)与对应的所述旋转θ角度后变形区域像素点坐标(x_i^*,y_i^*)做差，得到所述旋转θ角度后变形区域像素点相对于所述参考区域像素点的x方向的位移u和y方向的位移v；由所述旋转角度θ以及所述旋转θ角度后变形区域像素点相对于所述参考区域像素点的x方向的位移u和y方向的位移v生成变形区域参数

像素点坐标更新模块204，用于选取所述旋转θ角度后变形区域中任一像素点作为变形区基准点，根据所述变形区域参数以及所述旋转θ角度后变形区域中每一像素点到所述变形区基准点的距离(Δx^*,Δy^*)更新所述旋转θ角度后变形区域像素点坐标(x_i^*,y_i^*)，得到更新后的变形区域像素点坐标(x_i^*’,y_i^*’)；

像素点匹配相关系数第一计算模块205，用于根据更新后的变形区域像素点坐标(x_i^*’,y_i^*’)处的像素点的灰度值和所述参考区域像素点坐标(x_i,y_i)处的像素点的灰度值计算所述参考区域与所述更新后的变形区域的像素点匹配相关系数C(p)；

变形区域参数更新模块206，用于根据所述像素点匹配相关系数C(p)和所述变形区域参数更新所述变形区域参数得到更新后的变形区域参数

变形区域像素点坐标二次更新模块207，用于根据所述更新后的变形区域参数以及所述旋转θ’角度后变形区域中每一像素点到所述变形区基准点的距离(Δx^*’,Δy^*’)二次更新所述更新后的变形区域像素点坐标(x_i^*’,y_i^*’)，得到二次更新后的变形区域像素点坐标(x_i^*”,y_i^*”)；

像素点匹配相关系数第二计算模块208，用于根据所述二次更新后的变形区域像素点坐标(x_i^*”,y_i^*”)处的像素点的灰度值和所述参考区域像素点坐标(x_i,y_i)处的像素点的灰度值计算所述参考区域与所述二次更新后的变形区域的像素点匹配相关系数C’(p)；

第一判断模块209，用于比较所述二次更新后的变形区域的像素点匹配相关系数C’(p)与预设阈值的大小；若所述二次更新后的变形区域的像素点匹配相关系数C’(p)小于所述阈值，将所述更新后的变形区域的像素点匹配相关系数C(p)更新为二次更新后的变形区域的像素点匹配相关系数C’(p)，并将指令发送至变形区域参数更新模块206；否则，将指令发送至第二判断模块210；

第二判断模块210，用于判断是否将所述旋转θ’角度后变形区域中所有像素点作为变形区基准点，若没有将所述旋转θ’角度后变形区域中所有像素点作为变形区基准点，则将指令发送至变形区域参数二次更新模块211，否则将指令发送至三维坐标计算模块212；

变形区域参数二次更新模块211，用于更新所述变形区基准点，计算所述旋转θ’角度后变形区域中每一像素点到所述变形区更新后的基准点的距离(Δx^*’,Δy^*’)，根据所述更新后的变形区域参数以及所述旋转θ’角度后变形区域中每一像素点到所述变形区更新后的基准点的距离(Δx^*’,Δy^*’)二次更新所述变形区域参数中的u’,v’参数，得到二次更新后的变形区域参数将所述变形区域参数更新为二次更新后的变形区域参数并将指令发送至旋转θ角度后变形区域像素点坐标更新模块104；

三维坐标计算模块212，用于根据所述二次更新后的变形区域像素点坐标(x_i^*”,y_i^*”)计算所述旋转θ’角度后变形区域像素点的三维坐标；

变形后的叶片的应变计算模块213，用于根据所述旋转θ’角度后变形区域像素点的三维坐标和与所述旋转θ’角度后变形区域像素点对应的所述参考区域像素点的三维坐标计算所述变形后的叶片的应变，并根据所述变形后的叶片的应变判断和定位叶片的结构故障。

实施例三

将本发明用于一含裂纹缺陷的风扇叶片变形测量，并诊断其裂纹位置。根据本发明提供的测量方法，测量过程为：

对一铁质风扇叶片进行测量，将叶片背面锯出一条深度约为0.5mm的裂纹以模拟叶片裂纹缺陷。实验前将其表面喷洒黑白散斑颗粒，安装好后将风扇转轴固定在实验平台上，将双相机系统固定在叶片前端，调整相机的位置，确保叶片在旋转过程中的清晰成像，并标定好系统的内外参数。

将叶片在水平位置作为参考位置，并拍摄其图像。随后将叶片逆时针旋转大约20°的方位作为目标位置，并在目标位置叶片尖端的正面利用螺旋测微计分别加载至10mm、20mm、30mm、40mm使其产生变形，同时拍摄其图像。

提取沿径向方向的应变变形曲线，如图5所示，可看到在某处出现了明显的应变集中，根据这一特性同样可判断该处存在裂纹。

实施例四

将本发明用于一5kW的小型风力发电叶片(直径约5m)的运行过程中的相对变形测量，并进行故障分析。根据本发明提供的方法，测量过程为：

实验前用自来水将叶片表面清洗干净，晾干，然后在叶片表面制作对比度较好的黑白散斑图案。

摆放好相机系统的位置并确保风力发电机叶片清晰成像，两相机同时记录下风力发电机叶片的旋转过程，将叶片逆时针旋转90°、180°、270°、360°，拍摄其图像。本实验图像采集的频率为10Hz，叶片的转速约为0.134r/s。

提取叶片沿径向的应变变化曲线，如图6所示，同样也没有观察到应变集中或者突变的情况。于是可以判断叶片不存在裂纹或者鼓包等故障。

本文中应用了具体个例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处。综上所述，本说明书内容不应理解为对本发明的限制。