CN109886930B - 基于变化率与温度差异的热图像缺陷特征提取方法 - Google Patents

Abstract

本发明公开了一种基于变化率与温度差异的热图像缺陷特征提取方法，通过从涡流脉冲热图像中选出最大温度值像素点，根据最大温度值像素点得到变换步长，然后根据选定的步长提取代表性的瞬态热响应曲线。之后，再利用本专利所提出效率改进基于温度变化率与温度差异性的分类方法分类。最后采用Canny算子进行边缘轮廓提取，从而提取出涡流脉冲热图像的缺陷特征。

Description

基于变化率与温度差异的热图像缺陷特征提取方法

技术领域

本发明属于缺陷检测技术领域，更为具体地讲，涉及一种基于变化率与温度差异的热图像缺陷特征提取方法。

背景技术

红外热成像技术已广泛应用于汽车工业、造船业、石油化工工业以及航空航天领域，它能够有效的用于对缺陷的无损检测技术。其优势在于无需直接接触待检测试件，检测时间短等。

红外技术大体上可以分为两类，即主动加热和被动加热。对于主动加热来说，要求人为的将能量或者热量给予待检测试件。针对加热的方式不同，主动加热可以分为：光激励、电磁激励以及机械激励。光激励又可以分为脉冲激励和幅值调制激励两种，前者称为光脉冲热像法，后者称为光锁相热像法；电磁激励称为涡流热像法；而机械激励可以分为超声脉冲和幅值调制超声，前者称为超声红外热像法，后者称为超声红外锁相热成像法。对于被动加热，热量变化产生于待检测试件自身。

近年来，涡流等表面检测检测新技术得到飞速发展。其不损伤本体、快速高效等特性，能有效地解决传统无损检测方法存在劳动强度大、周期长、效率低、安全性差等问题，实现大面积快速检测、节省大量人力物力。涡流红外检测是基于电磁学中的涡流现象与焦耳热现象，运用高速高分辨率红外热像仪获取温度场分布，并通过对红外热图像序列的分析处理来检测结构缺陷及材料电磁热特性变化。其检测结果为图像，直观易懂，单次检测面积大，效率高，检测时无需接触被测件表面，同时可利用涡流效应检测表面及近表面缺陷，可检测更深层缺陷,这些都是这种检测方法的优势。

根据电磁感应定律，当通入高频的交变电流的感应线圈靠近导体试件时，在试件的表面会感生出涡流，如果被测件中有缺陷，涡流将被迫绕过缺陷，改变其流向，这将使得被测件内部涡流密度发生变化。由焦耳定律可知，涡流在被测件中转换成焦耳热，导致被测件中产生的热量不均匀，从而产生高温区和低温区，由于温度的差异性，高温区热量通过热传导向低温区传递，导致被测件不同区域温度发生变化，通过红外热像仪采集试件温度的变化过程，然后将采集的热图像视频交给计算机进行分析处理，来获取被测件相关信息，实现缺陷的定性与定量检测。

在现有技术中，名称为《脉冲涡流红外热图像的特征提取方法》的专利中，利用了步长搜索的方法进行缺陷特征的提取，在此之后采用了传统的模糊C均值的方法将瞬态热响应曲线分类。在该算法中，通过聚类中心和隶属度函数将曲线分类，由其目标函数可知，其分类原理是将样本与聚类中心之间的距离最小化，即每个样本的温度与聚类中心的温度差值最小化。然而这种方法在面对加热不均匀的情况时，就会对分类的准确性产生很大影响。因此，我们提出一种新的分类方法，不仅仅将温度的差异性作为分类的因素，同时还考虑的温度的变化率，即单位时间内温度的变化快慢。从而，使得分类结果更加准确。另外在本发明中，通过算法寻找中心点来初始化聚类中心，从而对原算法进行效率上改进。最后利用Canny算子提取缺陷边缘，使得结果可视化。

发明内容

本发明的目的在于克服现有技术的不足，提供一种基于变化率与温度差异的热图像缺陷特征提取方法，具有适应性强，准确性高，可视效果好等优点。

为实现上述发明目的，本发明一种基于变化率与温度差异的热图像缺陷特征提取方法，包括以下步骤：

(1)、将涡流脉冲热图像用三维矩阵S表示，其中，S(i,j,:)表示三维矩阵S的第i行和第j列，第三个维度表示时间；

(2)、从三维矩阵S选出像素值最大的点S(I_zz,J_zz,T_zz)，其中，I_zz、J_zz和T_zz分别表示最大像素值点的行对应值、列对应值和时间对应值；

(3)、设定K个温度阈值T(m)，m＝1,2,…,K，将最大像素值点S(I_zz,J_zz,T_zz)所在行进行温度划分，得到K+1个数据块，S^k(m,n,:)表示第k个数据块在m行n列的瞬态热响应值；

在第k个数据块中，找到温度最大值点，记为

设置第k个数据块的温度阈值THRE_CL^k，计算距离温度最大值点

最近的温度点

间的相关度Re，其中，j＝1,2,…,N，N表示最大像素值点所在行中像素点的总个数；

再判断Re是否小于THRE_CL^k，如果Re≥THRE_CL^k，则继续计算下一个距离次近的温度点间的相关度，直到得到Re＜THRE_CL^k时，计算结束，然后统计Re≥THRE_CL^k的温度点个数，记为CL^k，最后将CL^k最为第k个数据块的列步长；

(4)、设定P个温度阈值T(p)，p＝1,2,…,P，将最大像素值点S(I_zz,J_zz,T_zz)所在列进行温度划分，得到P+1个数据块；

在第

个数据块中，找到温度最大值点，记为

设置第

个数据块的温度阈值

计算距离温度最大值点

最近的温度点

间的相关度

其中，i＝1,2,…,M，M表示最大像素值点所在列中像素点的总个数；

再判断

是否小于

如果

则继续计算下一个距离次近的温度点间的相关度，直到得到

时，计算结束，然后统计

的温度点个数，记为

最后将

最为第

个数据块的行步长；

(5)、分块分步计算每一个温度点的瞬态热响应；

(5.1)、将最大瞬态热响应值存储在X(:,1)中，然后计算S^k(i,j,:)与X(:,1)间的相关度

(5.2)、设置阈值DD，集合X(:,g)；如果

则将S^k(i,j,:)作为一个新特征存储在X(:,g)中；否则，令

继续计算下一个与X(:,1)的相关度；如果i＞M，则令i＝i-M，j＝j+CL^k，即变化到第j+CL^k列进行计算，如果j＞N,则瞬态热响应的计算过程完毕；

(6)、将集合X(:,g)中的像素点分为L类；

(6.1)、对集合X(:,g)中的像素点进行初始化聚类；

(6.1.1)、计算每两个像素点之间的距离：

其中，

和

分别表示第

个像素点和第

个像素点的像素值，

g为集合X(:,g)中像素点数量；

(6.1.2)、初始迭代计数器i'＝1，i'∈[1，L]；设置半径th；

(6.1.3)、以每一个像素点

为中心，以半径th画圆，然后统计每个圆内像素点数量，找到包含像素点数量最多的圆，以该圆的圆心为第一个初始化聚类中心，记为V₁；

(6.1.4)、以V₁为中心，在集合X(:,g)中将

的像素点去除，然后令i'＝i'+1，在

的剩余像素点中，重新设置半径th，并继续重复步骤(6.1.3)，得到第二个初始化聚类中心，记为V₂；然后以此类推，当迭代计数器i'的计数值达到上限时，得到第二个初始化聚类中心，记为V_L；

(6.1.5)、将所有的初始化聚类中心用集合表示为：V₀＝(V₁,V₂,…V_i',…,V_L)；

(6.2)、对初始化聚类后的集合X(:,g)中的像素点分为L类；

(6.2.1)、设置聚类数目L，L满足：2≤L≤n；并初始化聚类中心V₀，初始化迭代次数c＝0；设定终止迭代条件阈值ε；

(6.2.2)、计算影响因子η_i'k'；

其中，x_k'(t)表示第t帧时第k'个像素点的像素值，V_i'(t)表示第t帧时第i'个聚类中心的像素值，b表示帧数；

(6.2.3)、利用公式

计算隶属度矩阵U；

其中，i'＝1,2,…,L，c∈L，d_n'k'＝||x_k'-V_i'||，n'＝i',j'，d_n'k'表示第k'个像素点与第i'次聚类中心V_i'的欧氏距离，x_k'表示第k'个像素点的像素值；τ为常数；u_i'k'表示第k'个像素点隶属于第i'类的程度；

(6.2.4)、更新聚类中心V_i'；

其中，g表示集合X(:,g)中的像素点总个数；

表示第k'个像素点的热响应值；

(6.2.5)、如果迭代次数到达最大值L或者前后两次聚类中心之差绝对值小于ε，则算法结束，并输出隶属度矩阵U和聚类中心V，再进入步骤(6.2.6)；否则，令c＝c+1，返回步骤(6.2.2)；

(6.2.6)、利用隶属度最大化准则对所有像素点去模糊化，得到每个像素点所属类别，即M_k'＝arg_i'max(u_i'k')；

(7)、对三维矩阵S进行降维处理；

(7.1)、计算第i'个类别中所有温度点瞬态响应的均值MCⁱ'；

(7.2)、计算MCⁱ'对应的瞬态响应值与第i^*个类别中第j^*个温度点瞬态响应值

间的相关度，记为

其中，i^*＝1,2…,L,i'＝1,2…,L,i^*≠i'，j^*＝1,2,…,K^*，K^*表示第i^*个类别中温度点的个数；

对第i^*个类别中得到的

求和，得到

再从所有的

中选出最大的

并记为

最后将

存在二维矩阵Y中；

(8)、将三维矩阵S变换为二维矩阵O，再对二维矩阵O和Y进行线性变换，即：

其中，

是Y的伪逆矩阵；

(9)、采用Canny算子算法对矩阵R进行特征提取；

(9.1)、选取一高斯滤波器

利用高斯滤波器对矩阵R进行平滑处理，即对矩阵R中每个像素点进行卷积运算：g₁(x,y)＝h(x,y,σ)*R(x,y)，R(x,y)表示矩阵R中坐标为(x,y)的像素点的像素值；

(9.2)、利用一阶偏导的有限差分法计算g₁(x,y)的梯度幅值；

(9.2.1)、计算梯度幅值G(x,y)：G_x＝g₁(x,y)-g₁(x+1,y+1)，G_y＝g₁(x+1,y)-g₁(x,y+1)，G(x,y)＝|G_x|+|G_y|，其中，G_x代表在X轴方向梯度，G_y代表在Y轴方向梯度；

(9.2.2)、计算幅值

其中，

是以像素点(x,y)为中心8邻域像素点的连线，G_x,y是以像素点(x,y)的梯度方向的直线，G_a(x,y)为其两者交点的幅值；

(9.3)、幅值G(x,y)同G_a(x,y)进行比较，若G(x,y)>G_a(x,y)，则保留G(x,y)的值，否则将此时的幅值设置为0；然后将保留的幅值对应的像素点进行非极大值抑制，得到图像G₂；

(9.4)、对非极大值抑制后的图像G₂与预设的高阈值H_th和低阈值L_th进行判断；

如果图像G₂中某一像素点g₂(x,y)的梯度幅值超过高阈值H_th，则将像素点g₂(x,y)记为边缘像素点；

如果图像G₂中某一像素点g₂(x,y)的梯度幅值低于低阈值L_th，则将像素点g₂(x,y)删除；

如果图像G₂中某一像素点g₂(x,y)的梯度幅值介于高阈值H_th和低阈值L_th之间，则判断像素点g₂(x,y)的8领域空间内是否存像素点的梯度幅值高于高阈值H_th，若存在，则保留像素点g₂(x,y)，并记为边缘像素点；否则将像素点g₂(x,y)删除；最终得到一幅显现缺陷特征的图像。

本发明的发明目的是这样实现的：

本发明基于变化率与温度差异的热图像缺陷特征提取方法，通过从涡流脉冲热图像中选出最大温度值像素点，根据最大温度值像素点得到变换步长，然后根据选定的步长提取代表性的瞬态热响应曲线。之后，再利用本专利所提出效率改进基于温度变化率与温度差异性的分类方法分类。最后采用Canny算子进行边缘轮廓提取，从而提取出涡流脉冲热图像的缺陷特征。

同时，本发明基于变化率与温度差异的热图像缺陷特征提取方法。还具有以下有益效果：

(1)、本发明在分类时，通过选取最大密度圆的圆心来初始聚类中心，相比于传统分类方法随机产生初始聚类中心，本发明的方法更加接近真实的聚类中心，所以效率上会大大提高。

(2)。本发明在分类时，相比于传统模糊C均值算法，不仅仅考虑了温度的差异性，同时还考虑的温度变化率。充分的减小了对待检测试件加热不均匀所产生的误差；

(3)、本发明采用行列变步长搜索实现了高效提取试件中的缺陷信息，并准确的刻画缺陷轮廓，弥补了传统方法对于缺陷提取上的一些不足。

附图说明

图1是本发明基于变化率与温度差异的热图像缺陷特征提取方法流程图；

图2是试件的示意图；

图3是试件实际情况下在不同时刻的瞬态热响应曲线；

图4是试件通过本发明在不同时刻的瞬态热响应曲线；

图5是试件通过本发明进行缺陷特征提取的结果图；

图6是试件在不同时刻的混叠向量对应曲线；

图7是试件在T1时刻三种情况曲线归一化后对比图；

图8是试件在T2时刻三种情况曲线归一化后对比图；

图9是试件在T3时刻三种情况曲线归一化后对比图；

图10是试件通过本发明和ICA算法提取的缺陷特征对比图；

图11是本发明算法处理时间、传统算法与ICA算法处理时间对比图。

具体实施方式

下面结合附图对本发明的具体实施方式进行描述，以便本领域的技术人员更好地理解本发明。需要特别提醒注意的是，在以下的描述中，当已知功能和设计的详细描述也许会淡化本发明的主要内容时，这些描述在这里将被忽略。

实施例

图1是本发明基于变化率与温度差异的热图像缺陷特征提取方法流程图。

在本实施例中，本发明一种效率改进基于变化率与温度差异性的变步长涡流脉冲热图像缺陷特征提取方法，包括以下步骤：

S1、将涡流脉冲热图像用三维矩阵S表示，其中，S(i,j,:)表示三维矩阵S的第i行和第j列，第三个维度表示时间；

S2、从三维矩阵S选出像素值最大的点S(I_zz,J_zz,T_zz)，其中，I_zz、J_zz和T_zz分别表示最大像素值点的行对应值、列对应值和时间对应值；

S3、从大到小设定K个温度阈值T(m)，m＝1,2,…,K，将最大像素值点S(I_zz,J_zz,T_zz)所在行进行温度划分，得到K+1个数据块，S^k(m,n,:)表示第k个数据块在m行n列的瞬态热响应值；

在第k个数据块中，找到温度最大值点，记为

设置第k个数据块的温度阈值THRE_CL^k，计算距离温度最大值点

最近的温度点

间的相关度Re，再判断Re是否小于THRE_CL^k，如果Re≥THRE_CL^k，则继续计算下一个距离次近的温度点间的相关度，直到得到Re＜THRE_CL^k时，计算结束，然后统计Re≥THRE_CL^k的温度点个数，记为CL^k，最后将CL^k最为第k个数据块的列步长；

S4、从大到小设定P个温度阈值T(p)，p＝1,2,…,P，将最大像素值点S(I_zz,J_zz,T_zz)所在列进行温度划分，得到P+1个数据块；

在第

个数据块中，找到温度最大值点，记为

设置第

个数据块的温度阈值

计算距离温度最大值点

最近的温度点

间的相关度

再判断

是否小于

如果

则继续计算下一个距离次近的温度点间的相关度，直到得到

时，计算结束，然后统计

的温度点个数，记为

最后将

最为第

个数据块的行步长；

S5、分块分步计算每一个温度点的瞬态热响应；

S5.1、将最大瞬态热响应值存储在X(:,1)中，然后计算S^k(i,j,:)与X(:,1)间的相关度

S5.2、设置阈值DD，集合X(:,g)；如果

则将S^k(i,j,:)作为一个新特征存储在X(:,g)中；否则，令

继续计算下一个与X(:,1)的相关度；如果i＞M，则令i＝i-M，j＝j+CL^k，即变化到第j+CL^k列进行计算，如果j＞N,则瞬态热响应的计算过程完毕；

S6、将集合X(:,g)中的像素点分为L类

S6.1、对集合X(:,g)中的像素点进行初始化聚类；

S6.1.1、计算每两个像素点之间的距离：

其中，

和

分别表示第

个像素点和第

个像素点的像素值，

g为集合X(:,g)中像素点数量；

S6.1.2、初始迭代计数器i'＝1，i'∈[1，L]；设置半径th；

S6.1.3、以每一个像素点

为中心，以半径th画圆，然后统计每个圆内像素点数量，找到包含像素点数量最多的圆，以该圆的圆心为第一个初始化聚类中心，记为V₁；

S6.1.4、以V₁为中心，在集合X(:,g)中将

的像素点去除，然后令i'＝i'+1，在

的剩余像素点中，重新设置半径th，并继续重复步骤S6.1.3，得到第二个初始化聚类中心，记为V₂；然后以此类推，当迭代计数器i'的计数值达到上限时，得到第二个初始化聚类中心，记为V_L；

S6.1.5、将所有的初始化聚类中心用集合表示为：V₀＝(V₁,V₂,…V_i',…,V_L)；

S6.2、对初始化聚类后的集合X(:,g)中的像素点分为L类；

S6.2.1、设置聚类数目L，L满足：2≤L≤n；并初始化聚类中心V₀，初始化迭代次数c＝0；设定终止迭代条件阈值ε；

S6.2.2、计算影响因子η_i'k'；

其中，x_k'(t)表示第t帧时第k'个像素点的像素值，V_i'(t)表示第t帧时第i'个聚类中心的像素值，b表示帧数；

S6.2.3、利用公式

计算隶属度矩阵U；

其中，i'＝1,2,…,L，c∈L，d_n'k'＝||x_k'-V_i'||，n'＝i',j'，d_n'k'表示第k'个像素点与第i'次聚类中心V_i'的欧氏距离，x_k'表示第k'个像素点的像素值；τ为常数；u_i'k'表示第k'个像素点隶属于第i'类的程度；

S6.2.4、更新聚类中心V_i'；

其中，g表示集合X(:,g)中的像素点总个数；

表示第k'个像素点的热响应值；

S6.2.5、如果迭代次数到达最大值L或者前后两次聚类中心之差绝对值小于ε，则算法结束，并输出隶属度矩阵U和聚类中心V，再进入步骤S6.2.6；否则，令c＝c+1，返回步骤S6.2.2；

S6.2.6、利用隶属度最大化准则对所有像素点去模糊化，得到每个像素点所属类别，即M_k'＝arg_i'max(u_i'k')；

S7、对三维矩阵S进行降维处理；

S7.1、计算第i'个类别中所有温度点瞬态响应的均值MCⁱ'；

S7.2、计算MCⁱ'对应的瞬态响应值与第i^*个类别中第j^*个温度点瞬态响应值

间的相关度，记为

其中，i^*＝1,2…,L,i'＝1,2…,L,i^*≠i'，j^*＝1,2,…,K^*，K^*表示第i^*个类别中温度点的个数；

对第i^*个类别中得到的

求和，得到

再从所有的

中选出最大的

并记为

最后将

存在二维矩阵Y中；

S8、将三维矩阵S变换为二维矩阵O，再对二维矩阵O和Y进行线性变换，即：

其中，

是Y的伪逆矩阵；

S9、采用Canny算子算法对矩阵R进行特征提取；

S9.1、选取一高斯滤波器

利用高斯滤波器对矩阵R进行平滑处理，即对矩阵R中每个像素点进行卷积运算：g₁(x,y)＝h(x,y,σ)*R(x,y)，R(x,y)表示矩阵R中坐标为(x,y)的像素点的像素值；

S9.2、利用一阶偏导的有限差分法计算g₁(x,y)的梯度幅值；

S9.2.1、计算梯度幅值G(x,y)：G_x＝g₁(x,y)-g₁(x+1,y+1)，G_y＝g₁(x+1,y)-g₁(x,y+1)，G(x,y)＝|G_x|+|G_y|，其中，G_x代表在X轴方向梯度，G_y代表在Y轴方向梯度；这样通过一阶偏导的有限差分法可以检测出图像中的像素点是否是边缘像素点，也就是判断一个像素点是否在斜坡上；

S9.2.2、计算幅值G_a(x,y)：

在本实施例中，如图2所示，

是以像素点(x,y)为中心8邻域像素点的连线，即田子型连线；G_x,y是以像素点(x,y)的梯度方向的直线，即穿过田子型连线的直线；G_a(x,y)为其两者交点的幅值；

S9.3、对梯度幅值进行非极大值抑制，仅仅得到全局的梯度并不足以确定边缘，因此，为了确定边缘，必须保留局部梯度最大的点，且要抑制非极大值点。下面我们将幅值G(x,y)同G_a(x,y)进行比较，若G(x,y)>G_a(x,y)，则保留G(x,y)的值，否则将此时的幅值设置为0；然后将保留的幅值对应的像素点进行非极大值抑制，得到图像G₂ g₂(x,y)；

S9.4、对非极大值抑制后的图像G₂ g₂(x,y)与预设的高阈值H_th和低阈值L_th进行判断；

如果图像G₂中某一像素点g₂(x,y)的梯度幅值超过高阈值H_th，则将像素点g₂(x,y)记为边缘像素点；

如果图像G₂中某一像素点g₂(x,y)的梯度幅值低于低阈值L_th，则将像素点g₂(x,y)删除；

如果图像G₂中某一像素点g₂(x,y)的梯度幅值介于高阈值H_th和低阈值L_th之间，则判断像素点g₂(x,y)的8领域空间内是否存像素点的梯度幅值高于高阈值H_th，若存在，则保留像素点g₂(x,y)，并记为边缘像素点；否则将像素点g₂(x,y)删除；最终得到一幅显现缺陷特征的图像。

实验仿真

为了能够更好的加强试件缺陷的信息，本实施例针对图2(a)和图2(b)所示的试件1、试件2的圆孔型缺陷所采集到的热图像序列进行了处理。采用本发明中描述的首先确定变步长搜索，然后针对瞬态热响应曲线运用改进的分类算法分类，然后进行反混合矩阵伪逆矩阵向量处理，最终通过Canny算子处理得到边缘提取结果。

下面分别利用本发明所述方法和ICA算法对试件进行计算分析。

本发明所述方法首先采集试件的原始数据序列，热序列的采样时间是15秒，两个温度阈值分别设置为：T(1)＝29，T(2)＝30，Ref_cl＝0.993。大于阈值Ref_cl的瞬态热响的个数是10。REFR¹＝0.97,REFR²＝0.95,REFR³＝0.93。在进行初始化聚类中心步骤中，密度圆的半径th＝25。在实际情况下提取试件的瞬态热响应如图3(a)、图3(b)和图3(c)所示。在新算法情况下提取试件的瞬态热响应如图4(a)、图4(b)和图4(c)所示，然后再进行缺陷特征的提取，得到图5(a)、图5(b)和图5(c)所示的结果。在同样的基础上，ICA算法提取到试件的混叠向量如图6(a)、图6(b)和图6(c)所示。

通过选取的缺陷1位置的瞬态响应、混叠向量以及实际情况的比较，如图7所示，本发明所述方法以及ICA算法的峰值与曲线走势与实际情况基本相同，故本发文明所述方法同ICA一样可以提取相应的特征信息。

通过选取的缺陷2位置的瞬态响应、混叠向量以及实际情况的比较，如图8所示，本发明所述方法以与实际情况基本相同，然而ICA算法的结果却与实际结果存在差异性，因此本专利所述方法能够准确的提取相应的特征信息。

通过选取的缺陷周围区域的瞬态响应、混叠向量以及实际情况的比较，如图9所示，发明所述方法以及ICA算法的峰值与曲线走势与实际情况基本相同，故本发文明所述方法同ICA一样可以提取相应的特征信息。

最后，试件经过Canny算子处理之后，得到图10(a)所示的缺陷，而ICA算法处理之后，得到图10(b)所示的缺陷，通过对试件进行Canny算子处理后与ICA算法处理后的比较，新方法可以滤去更多的噪声，精确的提取缺陷轮廓，可视效果明显。

图11为效率改进算法、传统算法与ICA算法在处理相同数量样本时间长短的对比。通过对比得，改进算法的效率高于传统算法与ICA算法。

尽管上面对本发明说明性的具体实施方式进行了描述，以便于本技术领域的技术人员理解本发明，但应该清楚，本发明不限于具体实施方式的范围，对本技术领域的普通技术人员来讲，只要各种变化在所附的权利要求限定和确定的本发明的精神和范围内，这些变化是显而易见的，一切利用本发明构思的发明创造均在保护之列。

Claims (1)

1.一种基于变化率与温度差异的热图像缺陷特征提取方法，包括以下步骤：

(1)、将涡流脉冲热图像用三维矩阵S表示，其中，S(i,j,:)表示三维矩阵S的第i行和第j列，第三个维度表示时间；

(2)、从三维矩阵S选出像素值最大的点S(I_zz,J_zz,T_zz)，其中，I_zz、J_zz和T_zz分别表示最大像素值点的行对应值、列对应值和时间对应值；

(3)、设定K个温度阈值T(m)，m＝1,2,…,K，将最大像素值点S(I_zz,J_zz,T_zz)所在行进行温度划分，得到K+1个数据块，S^k(m,n,:)表示第k个数据块在m行n列的瞬态热响应值；

在第k个数据块中，找到温度最大值点，记为

设置第k个数据块的温度阈值THRE_CL^k，计算距离温度最大值点

最近的温度点

间的相关度Re，其中，j＝1,2,…,N，N表示最大像素值点所在行中像素点的总个数；

再判断Re是否小于THRE_CL^k，如果Re≥THRE_CL^k，则继续计算下一个距离次近的温度点间的相关度，直到得到Re＜THRE_CL^k时，计算结束，然后统计Re≥THRE_CL^k的温度点个数，记为CL^k，最后将CL^k最为第k个数据块的列步长；

(4)、设定P个温度阈值T(p)，p＝1,2,…,P，将最大像素值点S(I_zz,J_zz,T_zz)所在列进行温度划分，得到P+1个数据块；

在第

个数据块中，找到温度最大值点，记为

设置第

个数据块的温度阈值

计算距离温度最大值点

最近的温度点

间的相关度

其中，i＝1,2,…,M，M表示最大像素值点所在列中像素点的总个数；

再判断

是否小于

如果

则继续计算下一个距离次近的温度点间的相关度，直到得到

时，计算结束，然后统计

的温度点个数，记为

最后将

最为第

个数据块的行步长；

(5)、分块分步计算每一个温度点的瞬态热响应；

(5.1)、将最大瞬态热响应值存储在X(:,1)中，然后计算S^k(i,j,:)与X(:,1)间的相关度

(5.2)、设置阈值DD，集合X(:,g)；如果

则将S^k(i,j,:)作为一个新特征存储在X(:,g)中；否则，令

继续计算下一个与X(:,1)的相关度；如果i＞M，则令i＝i-M，j＝j+CL^k，即变化到第j+CL^k列进行计算，如果j＞N,则瞬态热响应的计算过程完毕；

(6)、将集合X(:,g)中的像素点分为L类；

(6.1)、对集合X(:,g)中的像素点进行初始化聚类；

(6.1.1)、计算每两个像素点之间的距离：

其中，

和

分别表示第

个像素点和第

个像素点的像素值，

g为集合X(:,g)中像素点数量；

(6.1.2)、初始迭代计数器i'＝1，i'∈[1，L]；设置半径th；

(6.1.3)、以每一个像素点

为中心，以半径th画圆，然后统计每个圆内像素点数量，找到包含像素点数量最多的圆，以该圆的圆心为第一个初始化聚类中心，记为V₁；

(6.1.4)、以V₁为中心，在集合X(:,g)中将

的像素点去除，然后令i'＝i'+1，在

的剩余像素点中，重新设置半径th，并继续重复步骤(6.1.3)，得到第二个初始化聚类中心，记为V₂；然后以此类推，当迭代计数器i'的计数值达到上限时，得到第二个初始化聚类中心，记为V_L；

(6.1.5)、将所有的初始化聚类中心用集合表示为：V₀＝(V₁,V₂,…V_i',…,V_L)；

(6.2)、对初始化聚类后的集合X(:,g)中的像素点分为L类；

(6.2.1)、设置聚类数目L，L满足：2≤L≤n；并初始化聚类中心V₀，初始化迭代次数c＝0；设定终止迭代条件阈值ε；

(6.2.2)、计算影响因子η_i'k'；

其中，x_k'(t)表示第t帧时第k'个像素点的像素值，V_i'(t)表示第t帧时第i'个聚类中心的像素值，b表示帧数；

(6.2.3)、利用公式

计算隶属度矩阵U；

其中，i'＝1,2,…,L，c∈L，d_n'k'＝||x_k'-V_i'||，n'＝i',j'，d_n'k'表示第k'个像素点与第i'次聚类中心V_i'的欧氏距离，x_k'表示第k'个像素点的像素值；τ为常数；u_i'k'表示第k'个像素点隶属于第i'类的程度；

(6.2.4)、更新聚类中心V_i'；

其中，g表示集合X(:,g)中的像素点总个数；

表示第k'个像素点的热响应值；

(6.2.5)、如果迭代次数到达最大值L或者前后两次聚类中心之差绝对值小于ε，则算法结束，并输出隶属度矩阵U和聚类中心V，再进入步骤(6.2.6)；否则，令c＝c+1，返回步骤(6.2.2)；

(6.2.6)、利用隶属度最大化准则对所有像素点去模糊化，得到每个像素点所属类别，即M_k'＝arg_i'max(u_i'k')；

(7)、对三维矩阵S进行降维处理；

(7.1)、计算第i'个类别中所有温度点瞬态响应的均值MCⁱ'；

(7.2)、计算MCⁱ'对应的瞬态响应值与第i^*个类别中第j^*个温度点瞬态响应值

间的相关度，记为

其中，i^*＝1,2…,L,i'＝1,2…,L,i^*≠i'，j^*＝1,2,…,K^*，K^*表示第i^*个类别中温度点的个数；

对第i^*个类别中得到的

求和，得到

再从所有的

中选出最大的

并记为

最后将

存在二维矩阵Y中；

(8)、将三维矩阵S变换为二维矩阵O，再对二维矩阵O和Y进行线性变换，即：

其中，

是Y的伪逆矩阵；

(9)、采用Canny算子算法对矩阵R进行特征提取；

(9.1)、选取一高斯滤波器

利用高斯滤波器对矩阵R进行平滑处理，即对矩阵R中每个像素点进行卷积运算：g₁(x,y)＝h(x,y,σ)*R(x,y)，R(x,y)表示矩阵R中坐标为(x,y)的像素点的像素值；

(9.2)、利用一阶偏导的有限差分法计算g₁(x,y)的梯度幅值；

(9.2.1)、计算梯度幅值G(x,y)：G_x＝g₁(x,y)-g₁(x+1,y+1)，G_y＝g₁(x+1,y)-g₁(x,y+1)，G(x,y)＝|G_x|+|G_y|，其中，G_x代表在X轴方向梯度，G_y代表在Y轴方向梯度；

(9.2.2)、计算幅值G_a(x,y)：

其中，

是以像素点(x,y)为中心8邻域像素点的连线，G_x,y是以像素点(x,y)的梯度方向的直线，G_a(x,y)为其两者交点的幅值；

(9.3)、幅值G(x,y)同G_a(x,y)进行比较，若G(x,y)>G_a(x,y)，则保留G(x,y)的值，否则将此时的幅值设置为0；然后将保留的幅值对应的像素点进行非极大值抑制，得到图像G₂；

(9.4)、对非极大值抑制后的图像G₂与预设的高阈值H-th和低阈值L-th进行判断；

如果图像G₂中某一像素点g₂(x,y)的梯度幅值超过高阈值H-th，则将像素点g₂(x,y)记为边缘像素点；

如果图像G₂中某一像素点g₂(x,y)的梯度幅值低于低阈值L-th，则将像素点g₂(x,y)删除；

如果图像G₂中某一像素点g₂(x,y)的梯度幅值介于高阈值H-th和低阈值L-th之间，则判断像素点g₂(x,y)的8领域空间内是否存像素点的梯度幅值高于高阈值H-th，若存在，则保留像素点g₂(x,y)，并记为边缘像素点；否则将像素点g₂(x,y)删除；最终得到一幅显现缺陷特征的图像。

Images