基于梯度变换的太赫兹图像超分辨率重建

郭佑东，凌福日，姚建铨，3

（1.华中科技大学武汉光电国家研究中心，武汉430074；2.华中科技大学光学与电子信息学院，武汉430074；3.天津大学精密仪器与光电子工程学院，天津300072）

引 言

由于太赫兹波具有良好的穿透性、安全性等优良特性，太赫兹成像在医学影像、安全检查、环境监测等方面具有广泛的应用前景［1-6］。然而在实际应用中，受到成像信号频率、成像环境、探测器、器件制造工艺等因素的影响，利用实际太赫兹成像系统采集的太赫兹图像往往存在信噪比低，边缘模糊严重等现象［7］，这使得太赫兹原始图像的质量无法满足人们的正常视觉要求，以及研究人员对图像进行进一步研究的需要。提高太赫兹图像质量最直接的方法就是升级系统硬件，但该方法不仅有技术难度，还面临成本约束［8］。因此，研究进一步改善已有太赫兹图像质量的方法，能够在一定程度上弥补现有装置的不足。在不升级成像系统的前提下，超分辨率（super-resolution，SR）重建技术提供了一种从低分辨率（low-resolution，LR）图像恢复出相同场景的高分辨率（high-resolution，HR）图像的实用、有效的低成本解决方案［9-10］。该技术可以增加现有图像的空间分辨率，从而克服传感器和光学成像的固有分辨率限制［11-12］。近年来，研究人员持续开展了利用超分辨率重建技术对太赫兹图像进行重建的研究，包括基于多帧图像的凸集投影（projections onto convex sets，POCS）算法、迭代反向投影算法、基于单帧的Lucy-Richardson迭代、形态学算法、反卷积方法［13-16］等等。

SR技术是一个病态逆问题，与多帧图像相比，单帧图像仅含有少量信息，病态重建问题更为复杂。其核心在于建立从LR图像到HR图像的非线性映射，由于LR图像信息不充足，因此，需要根据先验信息建立正则条件来对重建过程进行约束［12］。基于梯度变换的图像超分辨率重建方法，主要思想是根据图像的梯度先验信息对重建过程进行约束来恢复图像。由于人类视觉对图像边缘的敏感性，2007年，FATTAL等人首次提出利用图像梯度轮廓的先验信息来描述图像边缘的空间分辨率［17］。2011年，SUN等人利用通用高斯函数对边缘梯度轮廓进行建模，从而首次提出了基于梯度变换（gradient field transformation，GFT）的图像超分辨率重建的方法［18］。2015年，NGUYEN等人利用三角模型以及混合高斯模型（triangle and Gaussian mixed model，TGM）来改进边缘拟合函数以描述图像边缘梯度轮廓的非对称性，但是该模型无法有效拟合较为复杂的边缘［19］。2018年，QIANG等人进一步提出自适应梯度变换方法，使用双边高斯模型拟合图像边缘梯度轮廓，同时利用“差函数”定义梯度变换来优化低分辨率图像的梯度图像，并以此作为正则项约束LR图像与HR图像之间的关系，进行超分辨率重建，锐化图像边缘［20］。由于太赫兹信号的波长属于亚毫米级，相较于可见光图像，太赫兹图像的边缘较为模糊。因此，基于梯度变换的重建方法为这一问题的改善提供了一种合适的解决方案。

本文中基于实际太赫兹透射式成像系统，提出利用有理分形插值结合基于梯度变换的超分辨率重建算法对该系统所采集的单帧太赫兹图像进行超分辨率重建，验证并扩展了基于梯度变换的重建方法的适用性。针对所采集的太赫兹图像中存在的对比度较低和空间噪声较大等影响其梯度图像质量的因素，引入基于空间信息熵的直方图匹配技术和双边滤波器减弱这些不良因素对梯度图像质量的影响，提高了重建算法的性能。经过重建之后，太赫兹图像的细节信息得到恢复，边缘得到锐化，图像质量得到提高；同时当成像信号频率较低时，该方法对图像质量的提升效果更加明显。

1 超分辨率重建模型

1.1 观测模型

超分辨率图像重建的观测模型如图1所示。图1中，PSF指点扩散函数（point spread function）。由于设备和环境的原因，成像系统获得的图像通常是高分辨率图像经过几何变形、模糊、降采样等退化过程之后形成的低分辨率图像［21］，该过程可以由下式表示：

式中，*表示卷积，ILR为LR图像，IHR为HR图像，D为退化矩阵，↓为降采样，n为成像过程中的各类噪声。

Fig.1 Observation model of super-resolution

超分辨率图像重建是图1所示图像退化过程的逆过程。该过程是一个不适定问题，需要引入外部约束和先验信息进行优化求解。由LR图像估计HR图像的表达式为：

式中，IHR为重建后的图像，U（IHR）为正则项，参量β为加权平衡参量。减小β将减少正则项对图像重建的作用，当β=0时，该式退化为最大似然估计。基于梯度变换的重建方法将理想图像的梯度先验信息作为正则项（即（2）式中的 U（IHR））来约束（2）式，达到锐化图像边缘的目的。

1.2 基于梯度变换的超分辨率重建

基于梯度变换的超分辨率重建方法框架如图2所示。该方法首先对低分辨率太赫兹图像进行插值，初步提高图像质量，消除锯齿状边缘。插值后的太赫兹图像（Iu）的细节信息包含在其梯度图像（▽Iu）当中。然后提取Iu的梯度图像，通过梯度变换使该梯度图像更加尖锐，变换后的梯度图像（▽I^HR）将包含更多的细节信息。利用▽I^HR作为正则项并根据最大后验概率

Fig.2 The framework of super-resolution reconstruction for terahertz imaging based on gradient field transform

（maximum a posteriori，MAP）方法约束图像的重建过程可以获得高分辨率图像（IHR），达到锐化图像边缘的目的。▽越尖锐，IHR的边缘也越尖锐，同时也将包含更多的细节信息。受到太赫兹图像对比度较低以及空间噪声的影响，▽Iu的质量无法满足有效进行梯度变换的要求。本文中引入基于空间信息熵的直方图匹配技术和双边滤波器（bilateral filter，BF）对Iu进行处理，优化其梯度图像。

1.2.1 图像插值 受到成像系统效率的限制，太赫兹原始图像中的像素点较少，边缘不连续。从太赫兹原始图像中提取的梯度图像质量较差，无法对该梯度图像进行梯度变换。因此，本文中首先利用有理分形插值对原始图像进行处理，采用具有较高插值精度的有理插值函数和能够较好地描述图像纹理细节的分形插值方法，扩大图像规模的同时，能够更好地保持图像边缘，从而提高梯度图像的质量。分形插值通过一个特定的迭代函数系统生成分形插值函数而实现，假设原始图像所在的平面区域为 Ω=I×J=［a，b］×［c，d］，其中I，J分别表示原始图像水平方向和垂直方向像素坐标的集合，原始图像的像素点为｛（xi，yj，zi，j），i=1，…，N；j=1，…，M｝，记区间 Ii=［xi，xi+1］，Ji=［yj，yj+1］，其中 i∈Γ=｛1，…，N-1｝，分形插值函数的表达式为：

式中，φi（x）：I→Ii，φj（y）：J→Jj为插值区间映射，Fi，j：Ω→R为各插值点的灰度值映射，R为插值后图像的区域，si，j为尺度因子。有理分形插值将（3）式中的qi，j（x，y）定义为：

式中，Pi，j（φi（x），φj（y））和 Bi，j（x，y）分别为双变量有理插值函数和有理扰动基函数［22］。

1.2.2 图像增强 受到成像信号强度的影响，太赫兹原始图像的对比度较低，从插值后的太赫兹图像中提取的梯度图像质量较低，对该梯度图像直接进行梯度变换将无法得到理想的结果。因此，本文中利用基于空间信息熵的直方图匹配技术对插值后的太赫兹图像进行增强，提高对比度，优化梯度图像。考虑到图像灰度级的空间分布，将图像划分为多个子区域，并在每个子区域根据下式进行图像增强处理，可以实现图像的全局增强。

式中，k为子区域内的灰度级，Ik，e为利用（5）式对灰度级k进行映射后在经过增强的图像中的灰度级，Fk为灰度级映射函数。fl为一个关于l的离散函数，根据图像的空间信息熵计算［23］。Imin和Imax分别为输入图像最小和最大的灰度值。经过全局增强之后，再对图像进行离散余弦变换并对其变换系数根据下式进行加权实现局部增强，进一步提升图像的对比度。

式中，K为子区域内的灰度级数量，参量γ用于平衡全局和局部增强的程度，其取值范围为，本文中该参量取0.25以保证图像的增强效果，同时不会由于过度增强图像而放大其中的噪声。

1.2.3 双边滤波 由于成像过程中太赫兹信号存在衰减且不同采样点采集到的信号衰减程度不同，同时受到成像系统其它噪声（如探测器随机噪声、扫描系统电磁噪声等）的影响，太赫兹图像中的空间噪声将对梯度图像的质量产生较大影响。因此，本文中采用双边滤波器来消除空间噪声的影响，该方法的优点在于其能够在保持图像边缘的前提下对图像进行平滑，不会破坏图像的边缘梯度轮廓从而不影响后续的梯度变换，双边滤波器［24］的表达式为：

式中，W=3σq，i和 j为图像中像素点的坐标；σq和 σr为几何扩散和广度扩散，由实验确定，分别用于控制空间域和值域权值的衰减［24-25］，本文中取 σq=3，σr=20。

1.2.4 梯度变换 自适应梯度变换（adaptive gradient transform，AGT）利用双边高斯模型拟合图像边缘［20］，构造“差函数”形式的梯度变换优化梯度图像，使其更加尖锐，自适应梯度变换表达式为：

式中，σLR和σHR分别为LR图像的梯度图像和变换后梯度图像中梯度轮廓的锐度（拟合高斯模型的标准差），n为梯度的方向。σLR和σHR之间的关系［18］为：

式中，参量σLR利用LR图像中每个像素点及其八邻域内的信息计算［20］。

1.2.5 图像重建 根据MAP方法，利用优化后的梯度图像（▽I^HR）作为正则约束条件代入（2）式可以获得高分辨率太赫兹图像（IHR）［10，18］，其表达式为：

式中，▽IHR为重建后的图像的梯度图像。本文中，参量β取值为0.15。考虑到太赫兹光束的强度服从高斯分布，因此重建时所采用的退化矩阵D为高斯矩阵。利用梯度下降算法［18-19］优化（10）式可以获得高分辨率太赫兹图像。

2 实验结果与分析

2.1 图像采集

Fig.3 The THz transmission imaging systema—the illustration of THz imaging system b—the experimental setup of the THz imaging system

实验中采用的太赫兹透射式成像系统（Zomega Terahertz Corporation（USA））如图3所示。其中M1～M10为反射镜，OAP1和OAP2为离轴抛物面镜（off axis parabolic mirror，OAP），L1和 L2为聚乙烯透镜。飞秒激光器产生的激光光束被偏振分束器（polarized beam splitter，PBS）分为抽运光和探测光。半波板（halfwave-plate，HWP）用于改变抽运光和探测光的比例。延线由两个反射镜构成，用于调节抽运光和探测光之间的时间延迟。抽运光经过衰减器，同时被反射镜反射后照射到砷化镓（GaAs）光电导天线上产生太赫兹波。探测光被用于探测太赫兹波，氧化铟锡（indium tin oxide，ITO）导电玻璃和平衡探测器（balance detector，BD）共同构成成像系统的探测装置。测量过程中，系统的温度维持在21℃（294K）。成像系统工作频率范围为0.1THz～1.0THz。样品为一个金属制斩波器，其半径为5cm，原始图像的大小为55pixle×55pixle，插值后的太赫兹图像的大小为440pixle×440pixle。样品被放入遮蔽物中并被放置于2维电动平移台上，通过控制程序移动平移台实现对样品的逐点扫描。扫描过程中，探测器对每一个扫描点测量太赫兹光谱，再将提取的太赫兹光谱按扫描次序排列成一幅太赫兹图像。

探测器在非目标区域采集到的时域信号和频谱分别如图4a、图4b所示。在放置遮蔽物后，太赫兹信号会存在一定的时延和衰减。经过傅里叶变换之后，对频域信号进行滤波即可获得单频信号。本文中使用0.25THz，0.50THz，0.75THz图像来分析成像信号频率对重建性能的影响。

Fig.4 Terahertz signal and spectruma—terahertz signal b—terahertz spectrum

2.2 梯度图像提取

本文中采用离散梯度算子［-1/2，0，1/2］和［-1/2，0，1/2］T与图像进行卷积，分别获得每个像素点水平方向和垂直方向的梯度。以0.75THz图像为例，图5a为太赫兹原始图像的梯度图像，图5b为插值后的太赫兹图像的梯度图像，图5c为对插值后的太赫兹图像直接进行增强后（未引入双边滤波）提取的梯度图像，图5d为对插值后的太赫兹图像进行增强和双边滤波后提取的梯度图像，图5e为变换后的梯度图像。图5中白色线条为太赫兹图像的边缘梯度轮廓，白色线条越窄，越明亮，梯度的幅值越大，图像边缘越尖锐。

Fig.5 Gradient field（0.75THz）a—gradient field of original images b—gradient field of interpolated images c—gradient field of enhanced images（without filter） d—gradient field of images processed by enhancement and bilateral filtering e—transformed gradient field

由图5a、图5b可知，经过插值之后，部分图像细节得到恢复，边缘变得连续；由图5c可知，直接对经过插值的太赫兹图像进行增强虽然可以提高其边缘梯度轮廓的幅值，但也会放大图像中的噪声，从而降低梯度图像的质量；由图5d可知，引入基于空间信息熵的直方图匹配技术和双边滤波器对经过插值的太赫兹图像进行处理后，其梯度图像的质量大幅提高，从而可以使得梯度变换的效果更加理想；另外，从图5e中可以看出，梯度变换可以有效压缩梯度轮廓的范围（白色线条更细），提高梯度的幅值。

2.3 重建结果分析

经过插值处理之后，太赫兹图像的分辨率得到提高（如图6b、图6g和图6l所示），部分细节信息更加清晰。相比于0.25THz的太赫兹信号，频率为0.50THz和0.75THz的信号的光斑较小，因此0.50THz和0.75THz图像经过插值后边缘更清晰。总体而言，插值处理虽然恢复了太赫兹图像中的部分细节信息，增强了图像边缘的连续性，但插值并不能解决由于衍射极限的影响而产生的模糊问题。因此，在对图像进行插值之后，还需要对图像进行重建，锐化图像边缘。

Fig.6 The experimental results of terahertz image reconstructiona—original image（0.25THz） b—interpolated image（0.25THz） c—reconstructed image（0.25THz，Wiener） d—reconstructed image（0.25THz，Lucy-Richardson） e—reconstructed image（0.25THz，proposed） f—original image（0.50THz） g—interpolated image（0.50THz） h—reconstructed image（0.50THz，Wiener） i—reconstructed image（0.50THz，Lucy-Richardson） j—reconstructed image（0.50THz，proposed） k—original image（0.75THz） l—interpolated image（0.75THz） m—reconstructed image（0.75THz，Wiener） n—reconstructed image（0.75THz，Lucy-Richardson） o—reconstructed image（0.75THz，proposed）

利用有理分形插值，Wiener滤波，Lucy-Richardson迭代和基于梯度变换的方法对本实验中所采用的成像系统获得的太赫兹图像进行复原的结果如图6所示。对经过插值的太赫兹图像采用基于梯度变换的方法重建之后，其边缘更加尖锐，同时图像中的空间噪声并没有因为边缘锐化而被放大且没有出现Lucy-Richardson迭代算法处理后出现的振铃现象（如图6d、图6i、图6n所示）。其中0.75THz图像边缘最尖锐，这是因为3幅图像中获取0.75THz图像的太赫兹信号频率最高，波长最短，光斑最小。但由于其信号强度较弱，0.75THz图像的对比度不及0.25THz和0.50THz图像。由图6e可见，0.25THz图像边缘锐化的效果最明显，这是由于受到衍射效应的影响，和经过插值的0.50THz，0.75THz图像相比，插值后的 0.25THz图像的边缘较为模糊，因此，该方法在对0.25THz图像质量的提升更加明显。

本文中采用边缘强度（ICV）和平均梯度（Idef）两种评价标准［26］定量评价图像重建的效果。两种评价标准的计算式分别为：

式中，▽2为拉普拉斯算子，M和N表示图像的大小。

经过插值的太赫兹图像和重建后的图像的边缘强度和平均梯度计算结果如表1所示。由表1中的数据可知，该重建方法可以有效提高太赫兹图像的边缘强度和平均梯度。其中该方法对0.25THz图像的边缘强度的提升效果最为明显。0.50THz图像经过重建之后，边缘强度和平均梯度最高，这是因为和0.25THz图像相比，其成像信号频率较高，光斑较小，边缘较为清晰。和0.75THz图像相比，虽然0.50THz图像成像信号频率较低，但其信号强度更大，因此图像的对比度更好。

Table 1 I CV and I def of interpolated and reconstructed images

综上所述，利用基于梯度变换的超分辨率重建方法可以有效恢复低分辨率太赫兹图像中的细节信息，锐化图像边缘，重建高分辨率太赫兹图像。由于需要提取较为清晰的梯度图像进行梯度变换，该方法对于空间噪声较小以及图像对比度较高的太赫兹原始图像将有更好的重建效果。

3 结 论

本文中针对实际太赫兹透射式成像系统采集到的原始图像，采用有理分形插值结合基于梯度变换的超分辨率重建方法进行了图像重建，利用基于空间信息熵的直方图匹配技术和双边滤波器优化梯度轮廓，提高了重建性能，分析了成像信号频率和强度对重建性能的影响。结果表明，该方法可以有效恢复太赫兹图像当中的细节信息，锐化图像边缘，提高图像质量且不会出现振铃现象。对经过插值的太赫兹图像采用基于梯度变换的超分辨率重建方法处理之后，0.25THz，0.50THz和0.75THz图像的边缘强度分别提高了169%，116%和104%，平均梯度分别提高了16%，28%和24%。