高空间分辨率高可见度的太赫兹光谱成像研究

褚致弘， 张逸竹， 曲秋红， 赵晋武， 何明霞*

1. 天津大学精密测试技术及仪器国家重点实验室， 天津 300072 2. 天津大学精密仪器与光电子工程学院， 天津 300072 3. 莱仪特太赫兹(天津)科技有限公司， 天津 300019

引 言

太赫兹波是一种频率范围在微波和红外之间的电磁波， 其振荡频率通常定义为0.1～10 THz[1]。 由于太赫兹波对非极性材料具有良好的穿透性、 物质在太赫兹波段的指纹谱特性以及光子能量低的特性， 使得太赫兹波适合应用在物质的鉴别、 安检安防、 非接触式无损检测等[2]。 太赫兹成像正是利用太赫兹波的独特性质与电磁波谱中其他波段的成像形成信息的互补， 如X射线、 紫外、 可见光、 红外、 微波[3]。 与红外等更高频段成像相比， 太赫兹波能透过塑料、 橡胶、 纸板等非极性材料进行层析成像， 探测内部物质[4-5]。 与微波毫米波等低频段成像相比， 太赫兹波的波长更短， 因此空间分辨率更高。

宽带太赫兹脉冲成像能够获得太赫兹频域光谱数据， 其除了具有二维的空间信息外， 还具有三维的光谱信息， 利用太赫兹的光谱数据能够对物质进行识别分类， 具有非常大的应用潜力。 太赫兹光谱成像可分为时域信息成像和频域光谱成像， 通常太赫兹光谱成像的极限分辨率体现在太赫兹频谱较高频区域光谱所成图像。 由于受到系统中激发太赫兹波的飞秒激光器、 斩波器、 锁相器等硬件的影响， 太赫兹频域中较高频段普遍存在能量低、 噪声大的问题， 直接导致了太赫兹较高频段的频谱成像图像模糊， 噪声干扰大， 影响图像的分辨率和对比度。 通过优化系统硬件降低太赫兹成像中的噪声问题成本高昂， 效果也很有限[6]。 因此， 采用图像去噪算法解决太赫兹频谱图像的噪声污染问题， 提高图像的边缘分辨率和整体质量尤为重要。

目前已知的太赫兹图像去噪方法大多采用传统图像去噪方法， 并且缺乏针对太赫兹频域光谱高频区域成像的去噪方法。 这类方法是基于空间域或变换域实现对太赫兹图像降噪。 例如基于空间域的非局部均值滤波、 中值滤波等方法实现对太赫兹图像去噪[7-8]， 其对图像像素点进行无差别处理会导致太赫兹图像损失边缘细节， 造成几何结构的破坏。 而小波变换这类基于变换域的太赫兹图像去噪算法[9-10]能够实现图像细节信息丰富的区域与较为平坦的区域区别处理， 但存在处理不同太赫兹图像时多尺度小波系数的最优化选择问题， 参数选择依靠经验， 处理后的灰度级偏离原始图像的灰度级， 影响在实际使用中的效率。 近年来深度学习技术的发展[11]， 越来越多的深度学习技术在图像去噪领域取得了良好的效果。 2017年Zhang等提出了DnCNN[12]， 其结合了ResNet中残差学习[13]、 批量标准化(BN)[14]、 修正线性单元(RELU)[15]的思想， 在图像盲去噪领域具有良好性能； 2019年Yuan等将残差学习的卷积神经网络应用在了高光谱图像去噪， 提出了HSID-CNN[16]； 基于深度学习的图像降噪技术相比于传统降噪算法能更好地将各像素点间的光谱信息相关联， 进而实现在降低图像噪声和保存图像边缘细节上取得更优的效果。

因此本文搭建了基于THz-TDS的太赫兹时域光谱成像系统， 实现了对高分辨率板成像。 并结合瑞利判据对成像系统的空间分辨率和景深进行标定， 找到了极限分辨率和最佳成像范围。 通过对采集得到的光谱数据进行多种方式成像分析， 研究高空间分辨率的太赫兹光谱成像算法。 最后训练出了针对太赫兹频域光谱信号中的高频区域所成图像的太赫兹深度去噪模型， 实现了在保留高频边缘细节特征的基础上对图像进行降噪， 并与传统太赫兹图像去噪算法的结果进行横向比较， 为太赫兹频域光谱高频区域的图像降噪提供了一种新的解决方案。

1 实验部分

1.1 太赫兹成像系统的搭建及其景深、 分辨率的标定

使用的太赫兹时域光谱成像系统由德国Menlo Systems公司生产的TeraSmart太赫兹时域光谱仪(THz-TDS)、 电控二维平移台、 手动一维平移台、 离轴抛物面镜(镀金膜)、 分束镜、 计算机等组成， 成像系统实物图如图1所示。

图1 太赫兹时域光谱成像系统示意图Fig.1 Schematic diagram of terahertz time-domain spectral imaging system

在进行太赫兹光谱成像实验前， 首先需要对成像系统的瑞利距离(共焦参数)进行计算， 并通过移动Z轴平移台对理论计算结果进行验证， 找出最佳成像范围。 太赫兹光谱成像系统分辨物体细节能力的重要表征指标参数就是空间分辨率， 在本套太赫兹成像系统中由光导天线法产生的太赫兹波经离轴抛物面镜后会在成像中心处聚焦， 但由于衍射的影响， 所成的像是一个椭圆形的艾里斑， 而不是理想状态下的一个“点”。 因此， 太赫兹光谱成像系统对物成像的艾里斑角半径决定了本系统的分辨极限， 而在实际应用中， 不同人对艾里斑的重叠有不同的评价， 主观区分的说服力明显不足， 为了科学的表征本成像系统的分辨率， 需要给出一个公认的标准， 瑞利判据就非常适合本套系统的分辨率表征。 对于缝隙型光阑， 其极限分辨率为两图样重叠区域的中点光强度约等于单个衍射图样中心光强度的81%(对于点状物光阑， 约为73.5%)。

在此基础上对系统进行分辨率标定， 采用了Edmund公司的USAF 1951高分辨率板分辨率测试目标板进行测定。 测试卡由若干元素组成， 每组元素分为三条水平线和三条竖直线组成， 定义这些线的空格分隔为线宽。

实验通过调节Z轴的手动平移台， 每次步进为0.5 mm进行二维扫描成像， 计算出最佳成像范围景深。 在最佳成像范围内得到的图像与高分辨率测试卡的查找表对比， 得出系统扫描成像能分辨出的最大线对数(Frequency)和最小线宽(Line Width)。 每毫米线对是成像分辨率的一种单位， 表征成像系统在一毫米内能分辨出多少对线， 而线宽就是每根线间的宽度， 其中线对数越多、 线宽越小代表成像系统的分辨率越高， 线对数与线宽能够相互转换。

通过对Z轴相对距离2～11 mm共19个位置进行太赫兹光谱反射式扫描成像， 二维扫描平移台X轴、Y轴电机步进均为0.1 mm， 像素点136×161共计21 896个。 由于空气中水分子对太赫兹信号质量的影响较大， 实验环境由聚乙烯薄膜封闭， 内部充干燥空气， 使实验环境的空气湿度恒为10%， 温度25 ℃。 对各点时域信号进行快速傅里叶变换(FFT)， 获得其频域光谱数据， 取各点的0～5 THz范围频谱总强度值进行成像， 并对每张图像中的各组条纹进行瑞利判据， 对本成像系统的最佳成像景深范围进行标定。 以高分辨率板的线宽表征Z轴平移台在不同位置处频域总强度成像的极限分辨率， 进而得出系统最佳成像景深为2 mm， 如图2(a)所示。 图2(b1)—(e2)分别展示的是Z轴平移台相对距离2， 11， 5.5和6.5 mm的高分辨率板组别(0, 5)的瑞利判据曲线及其对应的成像结果。 太赫兹图像的瑞利分辨极限判别方法以图2(d1)和(d2)为例， 共标出1—7共7个序号及所对应的点位其中1和7两点表示相对背景的总强度， 2， 4和6表示的是缝隙形光阑的峰， 3和5两点分别表示的是两缝隙光阑图样的重叠区中间的光强度， 按照瑞利判据对缝隙形光阑鞍-峰比的定义， 鞍与峰的比值在0.81处为极限分辨率， 数学表达式为

图2 太赫兹成像系统景深及分辨率标定示意图(a)： 太赫兹成像系统景深标定图； (b1)： Z轴2 mm处的瑞利判据曲线； (b2)： Z轴2 mm处的频谱总强度成像； (c1)： Z轴11 mm处的瑞利判据曲线； (c2)： Z轴11 mm处的频谱总强度成像； (d1)： Z轴5.5 mm处的瑞利判据曲线； (d2)： Z轴5.5 mm处的频谱总强度成像； (e1)： Z轴6.5 mm处的瑞利判据曲线； (e2)： Z轴6.5 mm处的频谱总强度成像， 以上实验均在标准高分辨率板组别(0, 5)处进行Fig.2 Depth of field and resolution calibration diagram of terahertz imaging system(a)： Calibration map of depth of field for terahertz imaging system; (b1): Rayleigh criterion curve at 2 mm Z-axis; (b2): Total spectral intensity imaging at 2 mm Z-axis;(c1): Rayleigh criterion curve at 11 mm Z-axis;(c2): Total spectral intensity imaging at 11 mm Z-axis; (d1): Rayleigh criterion curve at 5.5 mm Z-axis; (d2): Total spectral intensity imaging at 5.5 mm Z-axis; (e1): Rayleigh criterion curve at 6.5 mm Z-axis; (e2): Total spectral intensity imaging at 6.5 mm Z-axis; All the above experiments were carried out at the standard high-resolution plate group (0,5)

(1)

式(1)中，a为本系统所计算瑞利判据的鞍-峰比，P1—P7分别表示点1—7的频域光谱强度值。

通过计算太赫兹频谱总强度图像的各组别鞍-峰比a判断出各幅图像极限分辨率所对应的组别， 对照高分辨率板的测试查找表获取其线宽， 即极限空间分辨率。

通过瑞利判据对成像系统的景深进行标定可以得出系统对Z轴平移台的相对距离较为敏感， 后续在扫描样品的时候可以根据待测样品的厚度等信息， 通过调节Z轴平移台使系统达到硬件上的成像最佳状态。

1.2 太赫兹光谱成像流程

基于THz-TDS的二维扫描反射式光谱成像主要分为时域信息成像和频域光谱成像。 在时域中主要有最大峰值成像、 最小峰值成像、 飞行时间差成像等， 在频域中主要有相位差成像、 频谱总强度成像、 某一频点成像等。

本节以成像系统最佳景深范围内的Z轴相对距离6 mm所扫描数据为例， 分别采用时域最大峰值、 频谱总强度成像、 频域某特定频点成像， 并对结果进行比对分析， 研究太赫兹光谱成像中高空间分辨率的成像算法。

图3为太赫兹成像系统Z轴相对距离6 mm时的二维扫描光谱数据及其所成图像。 图3(a)所示的是高分辨率板的缝隙型光阑处太赫兹频谱曲线和背景处太赫兹频谱曲线， 分别取其中较低频段处(0.89 THz)和较高频段处(1.94 THz)成像， 所选频段均在光阑材料和背景材料的幅值有明显区分。 图3(b)为太赫兹时域最大峰值成像， 最高分辨率为397 μm； 图3(c)为频域总强度值成像， 最高分辨率为315 μm； 图3(d)为频率为0.89 THz处幅值成像， 最高分辨率为397 μm； 图3(e)为频率为1.94 THz处幅值成像， 最高分辨率为157 μm。 利用瑞利判据应用在太赫兹波成像的空间分辨率计算式(1)进行理论值计算， 以1.94 THz处为例， 其中λ≈154.5 μm(f为1.94 THz)，l(反射焦距)=50.8 mm，D(离轴抛物面镜直径)=50.8 mm， 计算结果约为188.49 μm； 频率为0.89 THz成像的理论空间分辨率为411.14 μm。 理论计算结果与实测结果基本一致。

图3 标准分辨率板的太赫兹频谱及不同方式所成图像(a): 分辨率板的缝隙型光阑处和背景处的太赫兹频谱曲线； (b): 太赫兹时域最大峰值成像； (c): 频域总强度成像； (d): 频率为0.89 THz处幅值成像； (e): 频率为1.94 THz处幅值成像, 以上实验均在Z轴相对距离6 mm处进行Fig.3 Terahertz spectrum of standard resolution plates and images produced in different ways(a): Terahertz spectrum curve at the slit aperture of the resolution panel and at the background; (b): THz time-domain maximum peak imaging; (c): Total intensity imaging in the frequency domain; (d): Amplitude imaging at 0.89 THz; (e): Amplitude imaging at the frequency of 1.94 THz； All the above experiments were performed at a relative distance of 6mm in the Z-axis

由上述实验结果可知， 太赫兹频域光谱成像的空间分辨率整体优于太赫兹时域相关参数成像。 太赫兹频域较高频段数据所成图像的空间分辨率最高， 但所含噪声也很高， 整体信噪比低； 频域较低频段数据所成图像的整体噪声较少， 但空间分辨率较低。 因为随着频率升高， 太赫兹波长变短， 分辨率提高， 而受到THz-TDS的硬件限制， 如图3(a)所示信号随着频率升高会迅速衰减， 最终会淹没在系统噪声中。 因此， 为获取高空间分辨率高可见度的太赫兹图像需采用太赫兹频域较高频段光谱数据进行成像， 并需要降低图像的无用噪声。

1.3 深度卷积神经网络在太赫兹频谱成像的去噪算法研究

1.3.1 太赫兹频谱图像去噪网络结构

本文所设计的太赫兹频域图像深度去噪网络借鉴了ResNet中的残差学习的理论以及DnCNN的思想， 网络拟合的是残差而不是数据本身， 将模型的训练结果设置为训练集中的标准图像和经对标准图像加噪后的重建图像之间的残差图像， 并创新性地引入了成像系统内真实的太赫兹高频噪声。

本文所提出的太赫兹图像深度去噪网络是在DnCNN的基础上优化改进的。 网络共设计了20层， 激活函数全部采用ReLU， 其计算简单能节省模型的运算时间并能克服一定的梯度消失问题。 其中第1层为卷积+ReLU， 卷积核的大小为3×3×1， 其将图像映射成为64个特征图； 中间第2至19层为卷积+Bn+ReLU， 批量标准化(Bn)能加速网络收敛、 提高泛化能力， 其中每次的卷积核个数都为64个， 卷积核的大小为3×3； 最后1层为卷积层， 包含64个3×3大小的卷积核， 网络的最终输出为噪声图像中的噪声估计。

训练网络的损失函数数学表达式为

(2)

式(2)中： LossMSE为网络的损失函数， 反映的是训练集图像的残差与网络所预测的残差之间的均方根误差， 该值越小说明模型所预测的残差图像精度越高；f(y；Θ)表示网络所预测的残差值；y和x分别表示加噪图像和原始图像， (y-x)即表示图像的残差真实值。

1.3.2 训练数据集与实验环境

测试集为所搭建的太赫兹光谱成像系统对高分辨率板进行扫描得到的宽谱数据， 对该数据中频域较高频段图像进行盲去噪实验。 太赫兹频域深度去噪网络的训练和测试环境基于64位Windows10操作系统， 采用Python语言基于Tensorflow框架训练， 实验平台采用单块GeForce RTX 2080Ti显卡为模型训练进行加速。 网络的初始学习率设置为0.1， 每训练10个周期学习率降低10%， 共训练90个周期。

2 结果与讨论

模型训练完成后， 采用搭建的太赫兹光谱成像系统在Z轴相对距离6 mm处光谱成像数据进行测试， 取频域1.94 THz处太赫兹频谱数据测试。 为评价网络性能， 对比了四种传统的太赫兹图像去噪方法， 其结果如图4(a)—(f)。 其中图4(a2)表示成像系统对高分辨率板在Z轴相对距离6 mm处的频域1.94 THz原始光谱数据所成图像； 图4(b2)表示双边滤波(Bilateral filter)算法对太赫兹频谱图像去噪的结果， 其相较原图有一定的去噪能力， 并且保存了高频的边缘细节； 图4(c2)表示非局部均值(NL-means)去噪算法对原频谱数据的去噪结果， 随着方框半径的增大去噪效果会提高， 但边缘细节会越模糊； 图4(d2)表示原频谱数据经维纳滤波(Wiener filter)去噪算法处理后的成像结果； 图4(e2)表示中值滤波算法对原频谱数据的去噪结果， 邻域范围取14×14， 像素点周围被无差别处理， 边缘细节遭到破坏； 图4(f2)表示在提高对比度后经本算法处理后的重建图像效果； 图4(a1)—(f1)分别表示其所对应图像在组别(1, 3)处的瑞利判据。 与传统去噪图像对比可知， 经去噪算法重建后的图像可以看出， 在图像边缘细节没有遭到破坏的前提下去除了大多数噪点， 整体的对比度、 清晰程度得到大幅度改善， 组别(1, 3)的瑞利判据结果曲线也更加平滑。

图4 1.94 THz成像经不同算法的去噪效果对比及其组别(1, 3)处瑞利判据(a1): 原始图像瑞利判据； (a2): 原始数据成像结果； (b1): 双边滤波图像瑞利判据； (b2): 双边滤波去噪结果； (c1): 非局部均值滤波瑞利判据； (c2): 非局部均值滤波去噪结果； (d1): 维纳滤波瑞利判据； (d2): 维纳滤波去噪结果； (e1): 中值滤波瑞利判据； (e2): 中值滤波去噪结果； (f1): 本文算法瑞利判据； (f2): 本文算法去噪结果, 以上实验均在频率为1.94 THz处进行， 瑞利判据均在组别(1, 3)处展示Fig.4 Depth of field and resolution calibration diagram of terahertz imaging system(a1): Rayleigh criterion of original image; (a2): Raw data imaging result; (b1): Rayleigh criterion for Bilateral filtering; (b2): Denoising result of Bilateral filtering; (c1): Rayleigh criterion for NL-means filtering; (c2): Denoising result of NL-means filtering; (d1): Rayleigh criterion of Wiener filtering; (d2): Denoising result of Wiener filtering; (e1): Rayleigh criterion of Median filtering; (e2): Denoising result of Median filtering;(f1): Rayleigh criterion of the proposed algorithm; (f2): Denoising result of the proposed algorithm； All the above experiments were performed at a frequency of 1.94 THz, and the Rayleigh criterion was demonstrated at group (1,3)

由于本研究的去噪图像是真实的太赫兹频谱图像， 没有与其对应的干净图像作为参考， 因此对太赫兹频谱图像经不同去噪算法处理后的重建图像进行评价时需要引入无参考图像评价体系。 无参考图像评价分为主观评价和客观评价指标， 选取了瑞利判据(Rayleigh Criterion)鞍-峰比表征图像中缝隙型光阑边缘细节的噪声程度、 梯度； 用对比度(CON)表征图像整体的清晰度、 渐变层次。 表1给出的是经不同算法处理后的太赫兹频谱图像与原图的图像质量客观评价对比表， 经分析可知去噪算法在边缘细节的梯度、 对比度均优于传统太赫兹图像的去噪算法， 能够在保留局部细节特征的基础上， 去除高频噪声， 降噪效果显著。

表1 不同算法对太赫兹频谱图像去噪前后的图像质量评价对比Table 1 Comparison of image quality evaluation of terahertz spectrum image before and after denoising by different algorithms

3 结 论

搭建了三维可移动式太赫兹时域光谱成像系统， 对高分辨率板进行的太赫兹反射式成像， 分别在时域、 频域进行重建成像， 标定了成像系统的景深。 通过不同方式的太赫兹成像， 研究了高空间分辨率的太赫兹成像方式， 发现在太赫兹频域较高频段成像的空间分辨率最高， 极限空间分辨率可达157 μm， 景深在2 mm左右。

针对太赫兹频域较高频段的图像整体信噪比低、 对比度低的特点， 在基于DnCNN网络设计思想的基础上， 在训练集中引入成像系统内真实的“太赫兹残差噪声”， 优化网络深度， 提出了太赫兹深度去噪网络。 通过与传统太赫兹图像去噪算法进行对比， 分别从主观和客观两个方面评价了不同算法对太赫兹图像去噪后的结果。 实验结果表明， 本文算法能够在保留图像边缘细节特征的基础上去除高频噪声， 并提高了太赫兹图像整体的清晰度， 在瑞利判据鞍-峰比和对比度的图像客观评价指标中均优于传统太赫兹去噪算法， 为太赫兹频域较高频段的图像去噪提供了一种新的思路。