一类改进的红外无损检测技术热波图像序列处理方法

陈炜耿，柴 琳，赵雅涵

(东南大学自动化学院，江苏 南京 210096)

在航空航天、汽车、铁路及船舶的生产制造及其产品的使用过程中，难免会产生各种各样的缺陷，这将严重影响产品的质量和使用性能。作为保证产品质量和使用安全的有效手段，无损检测(nondestructive testing，NDT)技术越来越受到相关技术研究人员的重视。无损检测是指在不损伤或基本不损伤被检试件的情况下，探测被检试件内部或表面的各种不同类型的缺陷，获得被检材料的质量信息，目前正广泛应用于金属、非金属、复合材料及其制品的缺陷检测。红外无损检测技术是一种新兴的无损检测手段，与超声检测、射线检测等传统常规无损检测技术相比，具有速度快、非接触、检测面积大和可远距离检测等优点。随着计算机视觉和图像处理技术的发展，红外无损检测技术获得了快速发展，已经被人们接受并且得到广泛应用[1-2]。

在红外无损检测技术中，脉冲热成像检测(pulsed thermography，PT)是研究最多和最成熟的方法，也是使用最广泛的方法。它利用短时而高能的脉冲激励注入被检试件，考察激励后试件表面的温度变化从而检测出有效的缺陷信息。Ludwig等[3]对脉冲热成像检测中传热模型进行了三维热扩散分析。Busse等[4]提出并发展的调制辐射检测技术，又被称为锁相热成像检测(lock-in thermography，LT)。它利用谐波调制激励源在被检试件内部产生正弦波形热波，然后采用数字锁相技术在激励周期内通过采集图像得到被测试件表面和内部的温度变化信号，并对各个像素点处的温度变化信号进行一维傅立叶变换，提取其相位和幅值信息，最后绘制出被测试件表面各点温度变化的相位图和幅值图。Swiderski[5]指出相位图不受表面发射率和光照非均匀辐射的影响。Mulaveesala等[6-7]指出调制频率限制了探测深度。Maldague等[8]提出脉冲相位热成像检测(pulsed phase thermography，PPT)，也称为脉冲相位法。脉冲相位法保留了PT和LT的优点，克服了PT对加热均匀性的严格要求和LT处理时间较长的局限性，其提取温度变换信号进行一维傅立叶变换后的最大相位值，绘制出被检试件的最大相位图。但脉冲相位法抑制加热不均现象的效果还不是很理想，也不能很好地抑制噪声，为此需要一种能进一步抑制加热不均和噪声、改善图像质量的算法。

1 红外图像序列处理算法

1.1 脉冲相位法

红外无损检测就是在对被测试件进行热激励的同时采集其表面热图，得到的是一个如图1所示的图像序列[9]。

图1 红外图像序列

红外图像可用矩阵T(n)表示：

(1)

式中：n为离散时间变量，表示红外图像帧数；Ti,j(n)为第n帧图像上第i行第j列像素点的温度值或灰度值。同一像素点处的灰度值序列Ti,j(n)可表示为该像素点处的离散时间序列，即Ti,j(0),Ti,j(1),Ti,j(2),…,Ti,j(N-1)表示第i行第j列像素点处的离散时间序列。

与锁相热成像法只能获得一种频率下的检测结果相比，脉冲相位法的检测结果包含许多不同频率的成分。脉冲相位法通过对图像序列中的离散时间序列进行一维离散傅立叶变换，可以获得多个频率下的检测结果，从检测结果中可以提取温度变化的相位信息和幅值信息。离散傅立叶变换的公式为：

jImi,j(k)k=0，1,…,N-1

(2)

幅值(A)公式和相位(Φ)公式为：

(3)

(4)

提取的相位信息和幅值信息可以用图2表示，因为获得的是多个频率下的检测结果，所以对应的相位信息和幅值信息也分别是一个图像序列。

与相位图相比，幅值图会受材料表面发射率和光照非均匀辐射影响，且相位图探测的深度是幅值图的两倍，因此一般选择相位图作为结果图。脉冲相位法将最大相位图作为结果图，最大相位图意即在每个像素点处选择最大相位值作为此坐标的像素值。由于只考虑最大相位值，而不是使用所有频率相位值的叠加，缺陷也就不会因温度扩散而展现出不同的表观尺寸。

图2 幅值图序列和相位图序列

1.2 红外图像序列预处理

一般而言，对被检试件加热前，总是假设被检试件处于热平衡状态。而由于周围环境噪声等因素的存在，实际上试件不可能处于绝对的热平衡状态。因此，在使用脉冲光源激励之前，首先用红外热像仪记录一幅背景热图，然后再触发脉冲光源激励并同时记录原始热图序列，并将该原始热图序列与背景热图相减，以降低背景噪声造成的影响[10]。

图3所示为采集的背景热图、某帧原始热图和减背景处理后图像。从图中可以看出，执行减背景操作后，图像中缺陷的对比度明显增大，有利于提高检测效果。若无特殊说明，下文出现的图像均指减背景处理后的图像。

对于红外图像中存在的噪声问题，在预处理阶段使用中值滤波进行去噪。中值滤波是基于排序统计理论的一种能够有效抑制噪声的非线性信号处理技术。中值滤波的原理是选择一个模板，然后对模板中像素的集合由小到大进行排序，再用集合中的中位数替代原像素值。当模板在图像中进行移动遍历后，就可以很好地对图像进行平滑处理。中值滤波对脉冲噪声有良好的滤除作用，因为中值滤波的输出结果是由模板范围内像素值的中间值决定的，所以对孤立的像素值远不如平均值敏感，从而可以消除孤立的噪点。中值滤波在滤除噪声的同时，能够保护信号的边缘，使之不被模糊。这些优良特性是线性滤波方法所不具有的。

图3 原始热图与减背景处理后热图比较

二维中值滤波的公式为：

(5)

式中：f(x,y)和g(x,y)分别为原始图像和滤波后图像；W为模板，通常为3×3，5×5等奇数区域，可以是不同形状；a和b为模板W中的数值。

在对图像序列进行减背景和中值滤波处理后，还需要对离散时间序列进行拟合。如前所述，红外图像序列可以看成是由全部像素点的离散时间序列组成的，每个像素点处都有一个离散时间序列。在这一步使用Savitzky-Golay滤波器进一步消除噪声和平滑曲线。

Savitzky-Golay滤波器[11]，通常简称为S-G滤波器，是一种在时域中基于局域多项式最小二乘法拟合的滤波方法，同时也是光谱分析中常用的预处理方法之一。使用S-G滤波器进行平滑滤波，可以提高光谱的平滑性，并降低噪声的干扰，其最大的特点在于在滤除噪声的同时可以确保信号的形状、宽度不变。S-G平滑滤波的效果，随着选取窗宽不同而不同，可以满足多种场合的需求。

(6)

式中：E为均方误差；x(n)为拟合前离散时间序列在第n点处的值。

为了最小化E，对E中多项式的每个参数求偏导并使其为零：

(7)

即:

(8)

式中：ar为多项式系数；nk+r为多项式自变量的幂次方。

(9)

需要拟合的点数M、多项式阶次N和待拟合序列x(n)均已知，故可以求得多项式系数，从而确定多项式。拟合时阶次N取得太高无法起到消除噪声的作用，故阶次N最好取4～6[12]。常用的拟合方法还有最小二乘法，其与Savitzky-Golay滤波拟合方法所使用的原理相同，都是最小化误差的平方和，区别主要在于最小二乘法是对整体数据进行拟合，而S-G滤波器是先将数据分成一个个区间，然后对每个区间的局部数据分别进行拟合。图4为Savitzky-Golay滤波拟合算法和最小二乘法算法对离散时间序列拟合的对比图，其中拟合阶次均取5，S-G滤波器的M取5。

图4 拟合算法对比图

图中点划线是原始离散时间序列曲线，虚线是最小二乘法拟合后的曲线，实线是S-G滤波拟合后的曲线。从图4可以看出，相比最小二乘法拟合，由于Savitzky-Golay滤波拟合算法每次是对局部离散时间序列进行拟合，其能最大程度地保持温度变化曲线，虽然最小二乘法也能正确拟合曲线的趋势，却不能很好地保持信号的形状和宽度。

1.3 改进的脉冲相位法

经过1.2节的图像预处理后，对图像序列使用脉冲相位法得到图2所示幅值图序列和相位图序列，接下来的步骤与传统脉冲相位法有所不同，不是绘制最大相位图，而是对相位图序列使用主成分分析(principal component analysis，PCA)法进行进一步处理。采集到的图像序列的离散时间序列是一个实数序列，由离散傅立叶变换的共轭对称性可知，其傅立叶变换左右对称，故在N张相位图序列中，只需要取前N/2张相位图用于后续处理[12]，为了方便叙述，仍记为N张。

主成分分析法顾名思义，就是找出数据中最主要的成分，用数据里最主要的成分来代替原始数据，通常用于萃取特征、减少冗余。它通过正交变换将一组可能存在相关性的变量转换为一组线性不相关的变量，转换后的这组变量叫主成分。第一个主成分应尽可能多地包含数据的信息，并且每个后续主成分都应尽可能多地包含剩余数据的信息。主成分分析通过使投影后数据的方差最大来寻找正交变换的坐标系来保留尽可能多的原始信息。它是一种非参数分析，具有唯一解，独立于任何关于数据概率分布的假设。

假设图像序列中每张图像的分辨率为320×256，即行数为256，列数为320。图像序列既包含空间信息又包含时间信息，可以看成是三维数据，将其转化为二维数据矩阵就有两种方式[13]，一种是行数为256×320，列数为帧数N，记为A1；一种是行数为帧数N，列数为256×320，记为A2。由此可知A1和A2互为转置，本文使用A1作为后续步骤的输入矩阵。图5是对相位图序列进行主成分分析的算法流程。

下面对每一步分别进行说明：

1)对所有数据进行中心化，得到行数为256×320，列数为帧数N的矩阵X1。

X1=[A1-Amean]=[(I1-Imean),(I2-Imean),(I3-Imean),…,(IN-Imean)]

(10)

图5 算法流程图

2)计算协方差矩阵，令

(11)

其中S1是一个行数和列数都为256×320的对称阵，即行和列均为81 920的方阵，不利于计算。若令：

(12)

则S2是一个行和列均为N的方阵，其大小取决于帧数，在绝大多数情况下，N肯定是远小于81 920的，故使用S2进行接下来的计算。

3)对协方差矩阵进行特征值分解，特征值分解公式为：

(13)

式中：Di是对角阵，其元素是以降序排列的Si的特征值；Ui为特征向量矩阵。Ui的每一列都是Si的特征向量，排列顺序与Di中的特征值对应。由A1和A2互为转置可得：

(14)

基于式(10)、(12)、(13)和(14)，可以推出U1和U2的关系为:

(15)

由此可知计算S1还是S2并没有本质的区别，但与计算S1相比，计算S2容易得多。

4)将特征向量还原成结果图，令：

EigImage=(A1-Amean)·EigVec2

(16)

其中EigVec2为U2的第2列，即第二大的特征值所对应的特征向量，则EigImage为一个行数为256×320、列数为1的列向量。将其还原成行数为256、列数为320的图像，此图像即为所求结果图。之所以取第二大的特征值所对应的特征向量是因为第一大特征值所对应的特征向量体现的是温度下降的平均趋势，而第二大的特征值所对应的特征向量才体现热量对比信息[14]。尽管后面的其他特征向量也包含了一些热量对比信息，但在引入的同时会引入噪声，故只取第二大的特征值所对应的特征向量。

与传统的脉冲相位法采用最大相位图作为结果图不同，改进的脉冲相位法对得到的相位图序列采用主成分分析法，进行了进一步的处理。通过主成分分析，图像的加热不均现象得到进一步抑制，同时也再次滤除了噪声。由前述可知，用于主成分分析的相位图序列帧数仅为原始热波图像序列帧数的一半，与传统的脉冲相位法需要计算每个像素点处最大的灰度值相比，虽然多执行了主成分分析的运算操作，但实际增加的执行时间仍在可接受范围内。

2 实验结果

实验所用的红外图像序列使用FLIR公司的紧凑型多功能红外热像仪FLIR A35采集，图像分辨率为320 × 256，图像的采集频率为60Hz，热灵敏度为50mK，采集的图像序列总共800张，其中第1张为背景热图。图6为第240帧热图，图7为脉冲相位法处理后的结果图，图8为用本文改进后算法处理的结果图。

图6 第240帧热图

图7 脉冲相位法处理后的结果图

对比图6、图7和图8可知，减背景处理后的图像(图6)能观察到的缺陷数量有限，缺陷特征也不明显。经过脉冲相位法处理后的结果图，抑制了原始热图中加热不均的现象，能明显观察到原始热图难以观察到的左边第2列的缺陷，图中的缺陷特征也更加明显。不严重的加热不均可以看成是由描述热量变化的增益函数和描述加热均匀的均匀场所构成，在用脉冲相位法计算相位值时，式(4)中虚部除以实部相当于把增益抵消，消除了一部分加热不均，同理也能消除表面辐射的变化。经过本文改进后的算法处理后的结果图，加热不均的现象基本消失，能观察到前2张图像观察不到的最左边第3列缺陷，图像中的缺陷特性也是3张图中最明显的。一般而言，信号中的时不变成分，如增益项，不会被映射到主成分中，因为增益项体现于式(13)中的对角阵D2的特征值中，特征向量中一般再无公约数，故主成分分析法能够进一步将体现加热不均的增益项消除。此外，主成分分析使投影后数据方差最大的过程也是进一步消除噪声的过程，由式(13)可知，数据投影后的协方差矩阵是对角阵，消除了噪声。实验结果证明了本文算法能有效去除图像中加热不均现象和噪声的干扰，提高信噪比，凸显缺陷信息，有利于提高检测效果。

图8 本文算法处理后的结果图

3 结束语

本文改进的基于脉冲相位法的图像序列处理算法，用于对红外图像序列的无损检测，能够有效消除图像中加热不均和噪声的影响，提高缺陷的对比度。加热不均和噪声是影响红外无损检测效果的最主要的不利因素，改进算法通过使用主成分分析进一步抑制了加热不均和噪声的不良影响，有效地增强了图像的质量，为下一步缺陷的定量分析打下了基础。