基于盲去卷积算法的太赫兹时域光谱图像处理方法

刘 拓,王洪金, 何赟泽, TOMASZ Chady, 许 健, 李 响

(1.湖南大学 电气与信息工程学院，长沙 410082；2.天津大学 精密仪器与光电子工程学院，天津 300000；3.波兰西波莫瑞工业大学 电气工程系，什切青 70-310)

将太赫兹波用于树脂玻璃材料的无损检测研究，所研究的问题是如何利用图像处理方法提高太赫兹成像分辨率。太赫兹波束广泛应用于非金属材料以及非极性材料中，具有宽频带、低能量、高穿透等特性，在无损检测领域有着重要的应用[1]。太赫兹波束的高穿透特性在同等条件下比涡流携带了更多信息；而其光子低能的特性，又使其比X射线更难影响材料内部的组织结构[2]。太赫兹时域光谱(Terahertz Time-domain Spectroscopy, THz-TDS)则是材料与太赫兹脉冲相互作用后，测量得到的太赫兹场强随着时间变动的曲线，其信噪比较高、稳定性好，可以在室温下工作[3]。在医学领域，太赫兹图像可以对微波、红外、可见光、紫外线或X射线检测图像所获得的信息进行互补[4]；在安防领域，适用于安检的太赫兹设备已经逐步应用在车站、机场等，与太赫兹设备相匹配的算法也在不断地出现，不仅可以节省人力资源，还能降低安检成本[5]；在军事领域，太赫兹合成孔径雷达在重量、体积、功耗等方面具有较大的优势，更适合于极大带宽信号的产生与处理，更容易实现高空间分辨力的成像[6]。但另一方面，太赫兹成像系统的低分辨率一直是其应用于无损检测或其他领域的巨大障碍，大量相关的研究工作都致力于提高太赫兹系统的成像质量[7]。笔者以树脂玻璃为试验材料，研究了一种基于盲去卷积算法的图像处理方法，通过图像灰度范围调整、盲去卷积、滤波去噪等步骤处理其太赫兹时域光谱，以提高其成像的分辨率，进一步拓展了太赫兹在无损检测领域的应用范围。

1 试验系统介绍

使用Picometrix公司的型号为T-ray 4000的脉冲太赫兹时域频谱系统进行逐点扫描成像。太赫兹时域频谱系统工作原理示意如图1所示。通过飞秒激光发射的脉冲分为两支，其中一支通过光纤作用于太赫兹发射器中的太赫兹发射材料，促使其发射频率约为1.8 THz的太赫兹波束。该系统为反射式，即太赫兹源和太赫兹探测器位于试验材料的同一侧，波束从太赫兹源出发，因其波长在毫米到亚毫米级别，太赫兹波在大多数绝缘体中可无损失传播或存在衰减传播。当太赫兹波遇到具有不同介电常数的两个物体表面时，会在边界界面的上下表面分别发生反射，根据分界界面与被测物体表面的距离，太赫兹场强峰值会在皮秒级的时域频谱中发生变化，从而得知分界界面的存在。通过对被检测物体表面进行逐点扫描，最终能实现被检测物体内部的分层成像。太赫兹时域频谱系统外观如图2所示。

图2 太赫兹时域频谱系统外观

树脂玻璃刻槽能够在平底刻槽的底端形成一个良好的反射界面，因此文中采用树脂玻璃刻槽来标定基于太赫兹时域频谱的成像系统的空间分辨率。试验材料为树脂玻璃，厚度为5.08 mm。在其表面从左到右共设置12组人工缺陷，这些表面人工缺陷是利用数控机床铣削加工的矩形槽。对于该类型的样本，目的不是模拟真正的缺陷，而是用于操作校准、构建校准曲线、确认探头的设置或验证算法。在试验中，此类缺陷用于评估所用算法对最大分辨率或灵敏度的处理效果。每一组缺陷之间的距离为10 mm，每组两个缺陷之间的距离逐渐增大，且所有缺陷的深度都是1 mm。缺陷的宽度为2.05 mm，缺陷的长度为5 mm。每组每两个缺陷的距离逐渐增加(见表1)。试验材料结构示意如图3所示。

表1 每组缺陷之间的距离 mm

图3 试验材料结构示意

2 盲去卷积算法

实际的太赫兹成像结果可以看成是材料的样品函数与太赫兹波束的点扩散函数(Point Spread Function, PSF)的卷积，成像过程如式(1)所示[8]

T(x,y)=I(x,y)*F(x,y)+N(x,y)

(1)

式中：T(x,y)为观测图像；I(x,y)为真实图像函数；F(x,y)为太赫兹波束的点扩散函数；N(x,y)为图像的噪声函数。

显然，当太赫兹波束的点扩散函数已知时，可以将其代入式(1)，由于T(x,y)也是已知量，所以可以反导出I(x,y)的数值。在太赫兹图像中，I(x,y)是太赫兹波在物体中某一时刻的传播方程在垂直于传播路径平面上的投影[7]。

然而，太赫兹时域频谱在绝缘介质中的传播时空方程极为复杂，其不仅与光程有关，而且与太赫兹的波长有关。对于宽谱太赫兹时域频谱仪得到的数据，使用物理扩散模型解逆卷积则会非常复杂，产生不必要的莫尔条纹[7]，因此采用盲去卷积算法进行处理。

盲去卷积算法的基本原理是：根据最大似然图像模型以及其稀疏特性[9-10]，不断迭代解卷积核，从而达到在不需要知道PSF具体数据的情况下，通过迭代以得到去卷积后的理想图像与PSF具体的数值。该算法基于这样一种假设：成像的光子满足泊松分布，因此理想图像将会是一个满足泊松分布的模型。通过原先输入图像的先验概率，以及模糊图像的先验概率，运用大似然估计获得原始清晰图像的后验概率，通过这种方式让复原图像尽可能接近原始图像[11]。假定理想图像各像素点独立，可以联立所有像素点构造关于所求图像、所求PSF与当前图像的似然函数。该似然函数的值是从理想图像推导出的当前图像的概率，当此函数取得最大值时，就说明当前求出的理想图像和当前求出的PSF能够最大可能地从这一理想图像推导出当前的图像。也就是说，能从当前图像最大程度地复原出理想的图像。具体的处理步骤如下所述。

(1) 读取图像，并对图像进行预处理。因为PSF收敛时，图像一般已经卷积过度，故需要提前设定迭代次数。

(2) 由于PSF未知，保守地估计其为一个所有值均相等的矩阵，以此矩阵作为初始PSF，以预处理后的图像作为初始图像。

(3) 使用最大似然估计的方法，即在频率域中构造卷积核作为PSF，使得当前图像可以反向寻找出理想图像的概率最大。

(4) 使用此PSF，反向迭代出新的当前图像。

(5) 使用新的当前图像，反向迭代出新的PSF。

(6) 重复第(4),(5)步的操作，直到达到预先设置的次数为止，输出最终的PSF结果与图像，完成整个盲去卷积的过程。

将这些步骤以流程图的形式表示，如图4所示。

图4 盲去卷积算法流程

3 成像结果与处理

由于太赫兹波的特殊物理性质，其在刻槽分界表面反射时会产生明显的瑞利衍射。为了能够直观地观测这些现象对太赫兹检测结果的影响,检测结果均以B扫描结果呈现(见图5)。

图5 材料的B扫描结果

图5中，横坐标为当前探头的水平位置，且只截取了与缺陷处相关的坐标，纵坐标为时间轴，且只截取了与材料有关的时间范围。此外，其为一幅伪彩色图，颜色坐标代表太赫兹的场强，为相对单位。A处的蓝色线段倾斜的原因是材料没有被水平放置，而B区域的一对色块对应于材料中的一对缺陷。色块成对出现，且从左至右，每对色块之间的距离越来越大，所有的色块都在一条时间坐标的直线上，这也与材料中缺陷的特征相吻合。C区域所在直线为材料的下表面，在垂直方向上，这条红线浅色的部分都与中间的红色色块一一对应。浅色的部分其实是该缺陷边缘周围的太赫兹衍射造成的，波束在入射到缺陷表面周围的下表面时，会对缺陷产生影响，即虽然浅色部分对应的区域实际上是空气而不是表面，却依然有场强存在。如果太赫兹系统的定向性足够优越，那么图5中底部红线所在区域只会出现断续的红色与背景色。

分别取横坐标47.5，45.5，153 mm，将3条曲线的谷值对齐，观察试验材料对应的A扫描结果(见图6)。太赫兹探测器在每个点记录从太赫兹脉冲发出后320 ps内返回的太赫兹波在探测器内引起的光电流所代表的电场强度波形图，如图6所示。仪器输出的太赫兹场强度用太赫兹脉冲的场强作为参考，以相对场强的方式进行记录，因而以相对单位进行表示。

图6 试件2在不同x位置时的A扫描结果

图6的横坐标为时间(进行了一定的截取),纵坐标为太赫兹场强。紧邻场强波谷的波峰的物理意义为太赫兹波束从材料的上表面反射，而右起的波峰物理意义为太赫兹波束从材料的下表面反射。

黑色、红色、黄色曲线分别对应试验材料中远离缺陷处、缺陷缝隙处与缺陷处。由于黑色曲线远离缺陷，太赫兹波衍射的影响较小，位于两个波峰之间，可以理解为波束穿过了试验材料；红色曲线在缺陷缝隙处，与缺陷的距离接近，因此曲线在两个表面峰值之间出现了衍射造成的峰值；黄色曲线没有下表面，只有缺陷表面，同样是由于衍射，在下表面的位置出现了一个本来不应存在的波峰。

以图5中蓝色线段为基准，将图像对齐后仅裁剪缺陷部分的图像。又由于图像处理算法中一般不支持小于零的像素值，故亦需要调整图像的灰度，将灰度转化为0到255，调整后的图像如图7所示。

图7 材料在缺陷附近的灰度图像

图7中红色部位对应反射最为强烈的区域，即分界层表面。从图7可以看出，在太赫兹B扫描图中，分界上下出现明显的瑞利衍射现象；太赫兹的反射强度也并非迅速衰减，从而导致相近刻槽之间的边界模糊，左起第二组刻槽间的间隙不可分辨，第三组刻槽间的间隙勉强可以区别。由于只是改变了图像的灰度范围，故图像的形态没有改变，此外，很难忽略缺陷周围“X”状的深蓝色图形，而该图形由太赫兹波在有限尺寸透镜中的传播本质引起[12]。又由于纵坐标是时间轴，对应深度尺寸，无需进行去卷积，分别使用3×1、4×1与5×1大小的PSF对图7进行处理，迭代10次，结果如图8所示。

图8 材料使用不同大小PSF去卷积后的图像

从图8可以看出，3×1大小的初始PSF处理后的图像不如4×1大小初始PSF处理后的图像改善明显，而5×1大小的初始PSF去卷积已经过度，在缺陷的缝隙处出现了错误的深红色形状。因此，最终采用4×1大小初始PSF对图7进行处理。也可以通过预估，设置初始PSF为4×1大小，来观察图像在不同迭代次数下去卷积的结果(见图9)。

图9 材料经过不同次数迭代去卷积后的图像

从图9可以发现，去卷积过程中引入了不希望被看到的高频噪声。当迭代次数过少时，去卷积达不到应有的效果，而在迭代次数过多时，缺陷部分的图像则会发生混叠，高频噪声的影响也过大，难以被去除。因此，选择初始PSF大小为4×1，迭代次数为10次。去卷积后，还要对图像进行滤波操作，削弱或消除高频噪声的影响，去卷积后的图像滤波频谱如图10所示。

将去卷积后的图像进行快速傅里叶变换(FFT)， 图像中心的位置是低频部分，图像的4个角是高频部分。将图像频谱4个角处的值设置为0，仅仅保留图像中心处的十字部分，就达到了低通滤波的效果，将此频谱进行反FFT，即可得到最后的图像，将图7与处理后的最终图像进行对比(见图11)。

空间分辨率的比较是通过横坐标的变化来实现的，横坐标有变化就意味着分辨率有变化。由图11可知，在对图像进行处理后，缺陷部分的分辨率明显提高，每对缺陷的中间缝隙都更容易辨认，说明处理效果良好。证实了盲去卷积的算法可以提高分辨率、改善图像质量。

图10 去卷积后的图像滤波频谱

图11 材料处理前后的缺陷处图像

为了对图像的清晰程度进行评价和估计，也可以引入灰度平均梯度(GMG)这一指标[13]。对于M×N大小的图像，其定义如式(2)所示。

(2)

一般情况下，GMG与图像的清晰度成正相关，可以较好地体现出图像的边缘特性。使用式(2)分别求取处理前后图像的GMG，得到图像在进行处理之前的本指标值为3.567 0，而在处理之后的本指标值上升到了4.243 6。显然，经过这些处理步骤之后的图像边缘特性比处理之前的更好，处理过程有效地降低了图像的模糊程度，缩小了在太赫兹B扫描图像中可区分两区域间的间隙。此外，将缺陷所在时间剖面处的处理前和处理后的曲线进行对比(见图12)。

图12 材料处理前后的缺陷处时间剖面图像

图12是将图11上下两幅图像缺陷处的时间坐标剖开的结果，横坐标为探头水平位置，纵坐标为该处的灰度。且黑色曲线是图像处理之前的缺陷剖面，红色曲线是图像处理之后的缺陷剖面。需要说明的是，图12并非在图11上的某一固定时间坐标取值，而是依照以下原则：当横坐标对应的值为缺陷时，就直接取灰度峰值作为函数值；当横坐标对应的值没有缺陷时，就取其临近缺陷处的图像纵坐标作为纵坐标值，而取该处横纵坐标的灰度。

通过图12也可以清晰看到，红线的每对峰值相较于黑线的每对峰值更容易分辨，处理后的曲线峰值较大，而两峰之间鞍部处的谷值较小，尤其在左侧，每对缺陷的谷值跟峰值都拉开了更大的差距，这进一步证实了文中图像处理方法有着良好的效果。

4 结语

利用以盲去卷积为主的太赫兹时域光谱处理方法，不需要得到太赫兹波束的点扩散函数，只需要提供初始PSF的规模与迭代次数，具备良好的实用价值。将太赫兹时域光谱进行了灰度范围调整、去卷积处理、滤波去噪等操作，分别对处理前后的图像进行直观对比、处理前后的图像灰度平均梯度对比、处理前后的图像缺陷时间剖面处对比，证实了算法能够减弱太赫兹成像中由于瑞利散射引起的模糊，但是对于太赫兹传播中的x-t(空间-时间)耦合分离的作用尚有待提高。