Application of trust region method in infrared image sequence processing

红外无损检测技术具有适用范围广、检测效率高、非接触测量、且无需耦合等优点，在近30年来得到了迅速发展，并广泛应用于航空、航天、汽车、军工、新材料研究、核工业及电力系统等各个领域^[1]。它通过使用外部激励源对被测试件进行主动式加热，使试件内部缺陷以表面温度分布异常的形式表现出来，通过红外热像仪实时记录这种表面温度分布随时间变化的情况，将其转化为图像序列，再利用多种分析处理方法提取温度异常区域，实现对试件中缺陷检测的目的。然而在红外图像的成像过程中，受到加热不均匀、周围环境红外辐射等因素的干扰，导致红外图像序列存在着噪声大、对比度低等问题，严重影响着缺陷的检测与识别。目前，红外图像序列处理方法主要包括单帧图像处理和多帧图像处理两大类^[2-4]。其中，单帧图像的处理方法主要包括：对比度调整、图像增强、非均匀性校正等方法。然而，单帧图像反映的是某一固定时刻的温度分布情况，由于内部缺陷深度不同，单帧图像不足以反映全部缺陷的信息。因此，多帧图像处理技术逐渐成为红外图像序列的研究热点之一。基于多帧图像的处理方法非常丰富，包括：多项式拟合法、脉冲相位法、正则化法、相关系数法、最大温差对比度法等。这些方法综合考虑了图像序列中不同时刻的温度信息，因此，处理效果往往优于单帧图像处理效果。然而这些处理方法也存在着一些不足，例如：多项式拟合法对噪声敏感；脉冲相位法的频率需人为选择；正则化法无法消除受热不均匀带来的背景噪声；相关系数法和最大温差法，需要选择无缺陷区域作为参考点。为了弥补上述方法存在的不足，文中提出基于信赖域反射（Trust Region Reflective，TRR）算法去除背景噪声，结合主成分分析（Principal Component Analysis，PCA）算法对去除噪声后的图像序列进行缺陷特征信息提取，并采用基于区域生长的图像分割技术对缺陷区域进行分割的红外图像序列处理方法。

红外热波无损检测中，通常采用卤素光源进行激励，光源辐射到试件表面的能量存在差异，使得获取的红外图像序列存在不均匀的背景噪声，增加了缺陷检测的困难程度。通过采用特定的图像处理算法，可消除红外图像序列中不均匀的背景，提高红外图像的信噪比。

红外图像序列背景去除算法的基本思路是根据试件表面温度分布不均匀的情况，选定合适的曲面对其加以拟合，然后将拟合后的曲面从原始图像中减去，从而消除背景噪声。传统的去背景算法多基于线性搜索法，然而红外图像数据呈非线性分布，进行线性搜索时存在收敛速度慢、易陷入局部最优解等问题，导致拟合背景中存在局部过拟合现象而误将缺陷信息引入拟合背景中。TRR算法是求解非线性最优化问题的典型算法之一，相比于传统的线性搜索算法，TRR算法具有鲁棒性强，收敛速度快的优点，且该算法的全局收敛特性^[5]有效避免了缺陷信息被引入拟合背景之中。因此，非常适合非线性曲面的热不均匀背景拟合。背景温度拟合函数如公式（1）所示^[6-7]：

式中：（x,y）为像素点位置坐标；a,b,c,d,e,f,g,h,i,j,k为待拟合系数。

为方便公式表达，对待拟合系数和像素点位置坐标进行向量化描述，记系数向量<$ a,b,c,d,e,f,g,h,i,j,k$>为β，坐标向量<$ x,y$>为X，因此，公式（1）可简记为f (X;β)，采用TRR算法可计算出系数向量β。

TRR算法通过引入信赖域，即当前解β_i的一个邻域N_i（由当前步长S_i、解β_i和信赖域半径r_i描述，${N_i} = \left\{ {S{}_i \in {P_i}\left| {\left\| {{S_i} - {\beta _i}} \right\| \leqslant {r_i}} \right.} \right\}$），在迭代过程中，将整个求解域划分为多个信赖域子问题，并在每一次迭代过程中，根据最小值条件得到一个试探步长S_i，该步长接受与否取决于特定的评价标准。如果当前步长被接受，则将当前解β_i更新为β_i+S_i，同时保持当前信赖域半径r_i；否则，保持当前解β_i不变，并缩小信赖域半径r_i。如此，解或信赖域半径通过迭代过程进行更新，直至满足容差收敛条件▽F(β_i)<tol。误差函数F(β)定义如公式（2）所示：

式中：Y_i为X_i所对应位置处的温度数据；n为拟合点个数。

TRR算法通过解决信赖域子问题，求解出相应的步长S_i，并对当前解β_i和信赖域N_i进行更新。信赖域子问题可用公式（3）来描述：

式中：g_i为F(β_i)的梯度矩阵；H_i为F(β_i)的Hessian矩阵；S_i表示第i次迭代试探步长。为提高计算效率，加快算法的收敛速度，考虑采用F(β_i)的梯度方向▽F(β_i)，记为向量${p_1}$，以及Gauss-Newton搜索方向，记为${p_2}$，生成解平面P_i。${p_2}$可通过解方程（4）得到。

采用TRR算法确定系数向量β_best的计算流程图如图1所示。

Figure 1.  Flowchart for solving coefficient vector β_best by TRR algorithm

PCA是基于二阶统计特性的一种典型数据压缩算法。经TRR去背景算法预处理后的图像序列中，缺陷信号的能量将远大于噪声信号的能量。在PCA算法处理过程中，能量越大的信号对应的奇异值也越大，因此，可通过去除较小奇异值，滤除能量较小的噪声信号。然后，根据有效信号子空间重构图像矩阵，达到缺陷信息提取的目的。因此，采用PCA方法对红外图像序列进行后处理能够显著提高图像的信噪比^[8]。

每一帧红外图像可由一个M×N维的矩阵${A_i}$表示：

式中：p表示图像总帧数，图像序列可由一个M×N×p维的矩阵X表示。为了方便使用PCA算法，需将图像序列转化为MN×p维的矩阵，因此，每一帧图像数据需进行向量化处理并作为矩阵X中的一列，记为X_i：

对数据进行中心化，可使数据尽可能分布在原点附近，该过程可以提高数据的可视性和算法的效率，具体过程为将矩阵X的每列数据减去所有列的平均值。在数据中心化后，下一步可通过计算得到协方差矩阵C：

式中：$ \mathop { X}\limits^ - = \frac{1}{p}\sum\limits_{i = 1}^p {{{{X}}_i}}$。

PCA需对协方差矩阵C进行奇异值分解，并得到其特征向量。分解式如公式（7）所示：

式中：矩阵S为对角矩阵，其对角线元素是矩阵C的奇异值；矩阵U由矩阵S对角线元素所对应奇异值的左特征向量构成，矩阵V则相应地由其右特征向量构成。矩阵U的每一列向量用U_i来表示，即：

通过对矩阵U中每一列向量U_i进行向量化逆过程可重构包含特征信息的各个子空间图像矩阵K_i：

红外图像在经过特征提取后仍然会存在一些较小的噪声干扰，为进一步突出显示缺陷信息，可根据图像中包含的缺陷特征对缺陷区域进行分割。区域生长法是一种典型的区域分割算法，其思想是根据事先定义的生长准则，将像素或子域组合为更大区域，基本方法是将与“种子点”性质相似的邻域像素归并到种子点上，从而形成生长区域^[9]。在处理后的红外图像中，选择缺陷中心的像素点作为区域生长算法的初始种子点。其次，需要确定一个生长准则，为了避免发生欠分割与过分割现象，算法上采用基于4连接域的区域生长方式，同时定义基于灰度差的相似性测度来判断每个位置（$ x,y$）处的属性Q，属性Q定义如下：

式中：T_max为图像最大灰度值；α取典型值0.02；I_seed（x,y）为种子点灰度均值；I（x,y）为待判断点灰度值。

令f($ x,y$)表示图像数据阵列，区域生长算法可作如下描述：

（1） 选定初始种子点，并将其属性标记为1，其他点标记为0。

（2） 对种子点4连接域上点的灰度值，由相似性测度判断其属性Q，若满足Q=1的条件，则标记该点属性为1，否则为0。

（3） 将种子点所有的4连接域上属性为1的点添加到种子点中，重复步骤(2)至(3)，直到区域停止生长。

（4） 采用半径为2的圆形结构元素进行先开后闭的形态学运算消除噪点，平滑边界；标记各个连通区域，得到区域生长后的分割图像。

实验采用2个输出功率为1 kW的卤素灯做为热源，以长脉冲方式对试件进行激励，激励时长为6 s。并对降温过程的红外图像进行采集，总采集时间为60 s（共采集300帧图像，并取前100帧进行处理）。红外热像仪采用CEDIP Jade MWIR系列，探测波段为3~5 μm，焦平面阵列大小为320×240像素，噪声等效温差为室温20 mK，并设定热像仪采样频率为5 Hz。

实验试件为增强型碳纤维复合板（Carbon Fiber Reinforced Polymer，CFRP），厚度为11 mm，规格为250 mm×250 mm，为模拟缺陷情况，试件上共有3种不同直径（10 mm，14 mm 和 20 mm）的平底孔，每种直径平底孔有8个深度，深度范围为1 ~5 mm。试件平底孔分布的示意图如图2所示。

Figure 2.  Schematic diagram of test specimen

实验虽采用对称双灯激励方式，加热时试件表面在一定程度上仍然会存在受热不均匀的现象，表现为距离灯源较近的区域温度偏高，较远区域温度偏低，从而在试件表面形成不均匀的受热背景。该背景会掩盖试件红外图像中存在的缺陷信息。利用TRR算法对背景噪声进行曲面拟合，并在原始图像中将其去除，可使被掩盖的缺陷信息显现出来。图3为采用TRR算法对原始红外图像序列中第10帧图像进行背景拟合得到的背景拟合曲面（TRR算法经11次迭代，计算时间为1.3 s；传统的线性搜索方法经37次迭代，计算时间为2.9 s），相应的拟合函数系数向量β如表1所示。从图3中可见试件上下两侧表面温度存在着一定的温差，这种不均匀的受热现象使得存在于下方较深的缺陷信息被掩盖。

Figure 3.  Fitted background image with 10^th frame thermographic data of the original thermal image sequence

图4显示了TRR算法处理前后的第10帧红外图像第195行图像数据（该行缺陷较深，缺陷信号受背景噪声影响最严重，其位置由图2中红色虚线标出）。可以看出：经TRR去背景算法处理后，红外图像中不均匀的受热背景得到了一定程度的抑制，处理前在曲线中被掩盖的表征缺陷信息的异常峰值信号在处理后的曲线中显现了出来。为对两组数据进行评价，定义（10）所示的峰值信噪比（Peak Signal-to-Noise Ratio，PSNR）指标为：

Figure 4.  Comparison graph of image data in line 195 before and after TRR algorithm processing

式中：T_max为最大温度数据幅值；T_nc为无缺陷区域温度数据平均值；S_nc为无缺陷区域数据标准差。原始图像数据峰值信噪比为15.41 dB，经TRR算法处理后的图像数据峰值信噪比为22.42 dB。对比结果表明了TRR去背景算法的有效性。

PCA算法可提取整个图像序列中包含的缺陷信息，达到进一步去除图像噪声，增强图像对比度的目的。图5（a）是直接对原始图像序列进行PCA处理得到的结果图（该图为所有主成分中缺陷显示效果最好的第二主成分图）；图5（b）则是对原始图像序列进行背景去除后再采用PCA处理得到的结果图（该图为所有主成分中缺陷显示效果最好的第二主成分图）。对比图5（a）和图5（b）可以看出，结合PCA处理技术，去背景算法可以显著提高图像的信噪比和缺陷的检出率。

Figure 5.  Images of experimental results

图5（c）为图5（b）中第195行图像数据（与图4位置对应）。从图中可见：与图4相比，图5（c）中的噪声得到了良好的抑制。

通过引入信噪比（Signal-to-Noise Ratio, SNR）指标，可客观地对图像质量进行评价，SNR计算式如公式（11）所示^[10]：

式中：$\mathop {{T_{\rm{c}}}}\limits^ - $和$\mathop {{T_{\rm nc}}}\limits^ - $分别为缺陷区域和无缺陷区域像素值的平均值；S_nc为无缺陷区域像素值的标准方差。图5（a）的信噪比为7.21，而图5（b）的信噪比则为14.54。缺陷的检出率也随之提高。图5（a）中只有17个缺陷显现了出来，缺陷检出率为70.8%；而在图5（b）中全部24个缺陷被检测了出来，缺陷检出率为100%。

为了进一步消除噪声干扰，增强图像数据的可视性，采用区域生长算法将缺陷区域从背景中分离出来。图6（a）为对图5（b）进行分割后得到的二值图像。结果图表明：通过采用基于灰度差相似性测度的生长准则，区域生长算法可准确地分割出各个缺陷区域。图6（b）为图5（b）中图像数据的三维显示图，图6（c）则为分割后缺陷位置处图像数据的三维显示图。

Figure 6.  Results of defect segmentation

图6（a）中带有缺陷的位置和形态信息，而图6（b）、（c）中还带有缺陷信号的幅度信息。实现准确的缺陷区域分割，有利于对缺陷区域进行后续的定量研究。实验结果图表明，对经TRR和PCA算法处理后的图像，采用基于区域生长的图像分割方法可准确的分割出缺陷区域。

在红外热波无损检测图像序列处理过程中，针对图像存在不均匀的加热背景导致缺陷检出率降低的问题，文中提出了采用TRR去背景算法进行图像预处理，并结合PCA和区域生长算法实现缺陷检测及分割的方法。为检验方法的有效性，对含有人工缺陷CFRP试件的红外图像序列进行了处理。实验结果表明：TRR去背景算法可去除红外图像序列中不均匀的加热背景；PCA算法可对红外图像序列进行缺陷信息提取，从而提高了图像的信噪比；为方便后续对缺陷的定量研究，采用区域生长算法对缺陷区域进行了分割。处理结果表明，对经过去背景和PCA算法处理后的红外图像，采用基于区域生长的图像分割算法能够较准确地分割出各个缺陷区域。结合使用以上图像处理技术，可使图像信噪比由7.21提高至14.54，进而将缺陷检出率从70.8%提高至100%。