太赫兹波三维成像技术研究进展

王与烨，陈霖宇，徐德刚*，石 嘉，冯 华，姚建铨

(1.天津大学 精密仪器与光电子工程学院 激光与光电子研究所,天津 300072；2.天津大学 光电信息技术教育部重点实验室，天津 300072；3.陆军军医大学第一附属医院 神经外科，重庆 400038)

1 引 言

太赫兹(Terahertz，简称THz，1 THz=1012Hz)波段是指频率从100 GHz到10 THz，相应的波长从3 mm到30 μm，介于毫米波与红外光之间频谱范围相当宽的电磁波谱区域。由于该频段是宏观电子学向微观光子学过渡的频段，具有很多独特的性质，如许多生物大分子的振动与转动能级，以及很多凝聚态体系的声子和其他元激发都位于太赫兹波段，因此，可以通过特征共振与吸收对物质进行指纹识别。太赫兹波的光子能量(毫电子伏特)低，一般不会对生物组织与细胞造成损害，可以对生物体进行无损检测。其对许多在可见和红外波段不透明的非金属、非极性材料具有较强的穿透能力，可进行透射成像。因此，太赫兹技术在物理科学、生命科学、信息技术和国防科技等方面都具有广阔的应用前景。太赫兹波成像是太赫兹科学技术最主要的应用之一，当前已发展出了透射式[1]、反射式[2]、衰减全反射[3]、近场成像[4]等主动成像技术以及被动成像技术[5-6]，并在无损探伤[7]、爆炸物检测[8]、文物研究[9-10]、安检[11]、药品检测[12]等方面，尤其是生物医学方面如对皮肤癌[13]、脑胶质瘤[14]、乳腺癌[15]的检测应用方面取得了巨大的进展。但是，传统的太赫兹波二维成像技术只能反映样品的表面或整体信息，无法满足内部信息观测的需求。

三维成像技术是获取样品内部信息的有效手段之一，自提出以来，在微波、可见光、红外与太赫兹波等领域得到了验证。不同频段的电磁波谱各有其缺点：X射线光子能量较高，对生物体有电离辐射；红外光对样品的穿透深度较低；毫米波虽然透过率高，但其分辨率受到长波长的限制。而太赫兹波的穿透性好，分辨率能满足实际要求，因此利用太赫兹波进行三维成像有更好的效果。层析成像技术是一种主要的三维成像技术，其英文单词“Tomography”来源于希腊单词“tomos”与“graphia”的组合，分别是切面及描述的意思，因此，层析的意思就是对样品内部切面进行探测。自1997年太赫兹波层析成像首次成功实现以来[16]，太赫兹波三维成像技术得到了迅速的发展，一系列新的三维成像技术相继被提出。根据三维成像系统的结构以及原理，太赫兹波三维成像技术分为透射式和反射式两种方式。目前，较为成熟的透射式层析成像技术是太赫兹计算机辅助层析(Computed Tomography,CT)，其可以看作是X射线CT在电磁波段上的扩展。研究内容在理论上包括太赫兹波透射样品时发生的折射、反射以及衍射等过程模型的建立，基于不同模型对重建算法的优化改进等；在实验上包括系统成像速度的提高、突破瑞利长度对成像样品尺寸的限制、通过改进探测方式以降低光束折射与散射的影响等。反射式层析成像是根据太赫兹波在样品内部不同深度的反射信号传输延时不同，通过对样品内部反射信号进行处理得出其深度信息，从而实现层析成像。根据实现方式主要分为太赫兹飞行时间层析(THz Time-of-Flight Tomography,THz-TOF)、太赫兹光学相干层析(THz Optical Coherence Tomography,THz-OCT)、太赫兹调频连续波(THz Frequency Modulated Continuous Wave,THz-FMCW)雷达成像等。目前，国际上美国、法国、德国和日本等发达国家陆续开展了太赫兹波三维成像技术的研究，国内哈尔滨工业大学、首都师范大学、北京工业大学、电子科技大学和华中科技大学等单位也开展了相关方面的研究工作[17-21]。在三维成像采用的辐射源方面，飞行时间层析技术采用太赫兹时域光谱(Time-Domain Spectroscopy,TDS)系统，利用光电导天线或光整流法产生的脉冲太赫兹波作为光源。而早期的太赫兹CT实验也大多采用THz-TDS系统作为光源[22]，其优点是能够获取样品的光谱信息，实现光谱层析，但系统结构复杂且耗时。近年来，基于气体激光泵浦的连续太赫兹源、太赫兹参量振荡器、量子级联激光器、耿式振荡器等产生的太赫兹波也得到利用。在探测方面，多采用光电导取样以及电光取样的相干探测，也有课题组采用阵列式热释电探测器、肖特基二极管(Schottky-Barrier Diode,SBD)等探测器，而且采用太赫兹频率上转换的间接探测方法也有报道[23-24]。在成像方式方面，受限于探测器，绝大多数研究采用逐点扫描的方式，但为了提高成像速度，采用电控光采样、振镜扫描的成像系统以实现实时成像的研究也在逐步开展。

2 太赫兹计算机辅助层析

2.1 CT成像原理

计算机辅助层析成像的概念是由挪威物理学家Abel在1826年最先提出的，1917年，奥地利数学家Radon发展了Abel的思想并建立了CT的数学理论基础，提出了Radon变换的概念。由于20世纪70年代初拉东的成果才被人们发现，用于重建图像的代数算法是由美国物理学家Cormack于1963年确立的。1972年，英国EMI公司工程师Hounsfield研制了第一台用于临床的X射线CT扫描装置[25]。如今，X射线CT的发展已经经历了笔束平移旋转扫描、窄扇束平移旋转扫描、宽扇束射线源与探测器同步旋转扫描、宽扇束射线源旋转扫描、电子束扫描、螺旋扇束扫描、螺旋锥束扫描七代。太赫兹CT的概念来源于X射线CT，其基本原理与X射线CT相似，即：样品置于一个二维移动平台带动的旋转平台上进行透射式扫描，样品每旋转一个角度dθ则进行一次线投影扫描，理论上共旋转180°就可以完成一次样品切面的CT扫描，得到的数据在以横坐标为投影角度，纵坐标为投影位置的强度图上表示，称为正弦图。以Shepp-Logan头部模型为例，图1(a)为模型的切面图，其模拟投影的正弦图如图1(b)所示，以若干个投影角度对头部模型进行投影，将得到的数据通过一定的算法处理，可以重建出样品一个切面的分布信息，再经过纵轴的平移对多个切面成像，最后利用三维建模即可得到样品的空间分布信息。

图1 (a)Shepp-Logan头部模型；(b)模型的投影正弦图 Fig.1 (a)Shepp-Logan phantom; (b)Projection sinogram of Shepp-Logan

上述投影成像过程可以用拉东变换来解释，拉东变换就是实际的射线投影，即样品的某个物理参量(如吸收系数、折射率)在不同角度的线积分。以二维拉东变换为例，在θ角度线扫描获得的数据Rθ(ρ)可以表示为样品分布函数f(x,y)的拉东变换，即在投影路径上的积分：

(1)

其中，ρ为线扫描的位置，δ是狄拉克δ函数，满足：

(2)

为了获得样品的空间分布函数f(x,y)，需要对得到的数据Rθ(ρ)进行逆拉东变换，也就是数据重建。最常用的重建算法是基于傅里叶中心切片定理发展出的滤波反投影算法(Filtered Back Projection,FBP)。在傅里叶中心切片定理中，线投影函数Rθ(ρ)的傅里叶变换为：

(3)

对应的分布函数f(x,y)的二维傅里叶变换为：

(4)

其沿着θ角度的值与公式(3)等价。

进一步，将公式(3)代入公式(4)，可以得到样品空间分布函数的表达式为：

|ω|e2πiω(xcosθ+ysinθ)dωdθ.

(5)

其中，|ω|表示由笛卡尔坐标系变换到极坐标系时需要引入的滤波，其可以提高内部成像质量并消除伪影。对于点扫描成像系统，样品空间分布函数可以用离散形式表示为：

(6)

其中：

A(θ,ρ),(i, j)为每个像素及投影线间的权矩阵，|v|表示斜坡滤波器。

公式(5)与(6)所表示的即为FBP算法，这是目前效率最高、应用最广的解析重建算法，被大多数的X射线CT仪器所采用。但是这种算法对投影数较敏感，在投影数少的情况下，重建图像会出现伪影。除了解析算法，还有基于迭代的重建算法，如联合代数重建算法(Simultaneous Algebraic Reconstruction Technique,SART)、有序子集最大期望值法(Ordered Subsets Expectation Maximization,OSEM)等，这些算法在低投影数的情况下比解析算法效果更好，但它们需要大量的计算。

在太赫兹CT成像过程中，根据使用的太赫兹源类型，可以将其分为连续波太赫兹CT和脉冲太赫兹CT。对于连续波太赫兹CT，由朗伯-比尔定律(Lambert-Beer Law)可知：

(7)

其中，I0与I分别表示无样品时和经过样品投影后的光强。根据公式(5)可以得到样品吸收系数的分布，即α(x,y)=f(x,y)。

当采用基于THz-TDS系统的太赫兹脉冲进行层析成像时，不仅能获得光波的强度变化，还能得到太赫兹波的相位信息。通过傅立叶变换，可以在频域对样品折射率进行分析。投影后太赫兹波的相移Δφ满足：

(8)

其中，c表示空气中的光速，ω表示频率。根据公式(5)可以得到折射率实部的分布为nω(x,y)=fω(x,y)=1，由振幅衰减可以得到折射率虚部即吸收系数随频率与空间位置的分布。由此可见，相比于X射线CT，太赫兹CT可以获取更多的样品信息。

2.2 研究进展

2002年，B.Ferguson等人首次将太赫兹波应用到CT成像中，利用基于钛宝石飞秒激光器泵浦的THz-TDS系统对聚乙烯材料进行成像，飞秒激光入射到光电导天线产生脉冲太赫兹波，使用光谱仪和CCD相机来探测恢复太赫兹信号[22]，成像系统以及对火鸡骨的成像结果如图2所示。太赫兹CT的首次成功实现预示着其在无损探伤、生产过程诊断及监控等方面的巨大潜力。由于超短太赫兹脉冲在频域具有宽带性，采用THz-TDS系统能够实现太赫兹光谱层析。2010年，德国A.Brahm等人采用THz-TDS系统对装有葡萄糖和乳糖的聚苯乙烯泡沫进行光谱层析成像，通过与数据库的对比分析，可以分辨出材料中葡萄糖和乳糖的位置[26]。同年，日本Advantest公司E.Kato等人利用两个高功率超短脉冲光纤激光器及光电导开关来实现太赫兹脉冲的产生与探测。由于采用异步采样技术取代了传统的机械式光学延时线，系统的扫描频率可达250 Hz。实验中对装有不同浓度的乳糖(L)、酪氨酸(T)以及乳糖酪氨酸混合粉末(LT)的聚乙烯圆柱状物进行光谱层析成像，结合物质的折射率补偿方法，实现了对高折射率样品聚乙烯(n=1.53)的层析成像。通过对不同太赫兹频率的图像重建，可以分辨出不同物质的位置及浓度的差别，如图3所示[27]。该公司已将这种太赫兹CT技术做成了实用化的产品[28]。太赫兹光谱层析的成功实现极大地增加了CT成像所能获取的信息量，其能够对样品的内部成分进行精确识别及定位，这是太赫兹CT相比于X射线CT的一大优势。

图2 (a)基于THz-TDS系统的太赫兹CT实验装置图；(b)火鸡骨的实物图与太赫兹三维成像图[22] Fig.2 (a)Experimental setup for THz CT based on THz-TDS; (b)Physical map and 3D image of a turkey bone[22]

图3 装有乳糖和酪氨酸的聚苯乙烯层析成像(a)样品示意图和(b)层析成像图[27] Fig.3 (a)Schematic representation and (b)tomographic image of the cylindrical polystyrene with lactose and tyrosine[27]

当采用THz-TDS系统进行点扫描光谱层析时，通常需要一个样品二维移动平台与一个转动平台。特别是基于机械延迟线的THz-TDS系统普遍采用机械平移台所获取的光学延迟线进行采样探测。该方法虽然可以实现很长的精确时间延迟，但由于其机械惯性而不能实现快速扫描，而且现有的锁相放大器由于硬件本身的原因，在处理速度上难以突破瓶颈。显然，多个移动平台的使用无法满足一些快速实时测量场合的要求。成像时间过长成为限制太赫兹CT扩大其应用范围的主要因素之一。

2007年，澳大利亚X.X.Yin等人采用二维离散小波变换原理导出的逆拉东算法对感兴趣区域(Region of Interest,ROI)的投影数据进行处理，这种算法的优良局部特性可以在较少投影角的情况下实现对样品的局部成像，有效地减少了基于THz-TDS系统的CT成像时间[29]。第二年，他们又和英国剑桥大学的研究人员合作，基于量子级联激光器进行了局部感兴趣区域的图像三维重建[30]。2011年，日本大阪大学T.Yasui与法国波尔多大学E.Abraham等人将非共线二维自由空间电光采样(2D Free-Space Electro-Optic Sampling,2D-FSEOS)技术与实时线投影相结合，有效避免了机械延迟线和样品扫描所消耗的时间。对于360°旋转的样品，在0.6°旋转角步长下获得600个线投影数据，每一线投影所需时间为10 ms，即仅需6 s就可以完成一个切面的光谱层析成像，实现了快速的太赫兹光谱层析成像，成像系统如图4所示[31]。2013年，他们进一步采用飞秒激光泵浦LN晶体，产生高强度的太赫兹波输出，对四根金属棒、装牙签的塑料盒以及一个胶囊进行太赫兹CT成像。对于横向20 mm包含100个像素点，纵向包含116个像素点，旋转步长5°的样品，成像时间仅为6 min[32]。如果采用更快速的位移台，成像时间可以进一步缩短。2017年，北京工业大学B.Li等人采用扩束的太赫兹光斑透射样品，采用面阵式热释电探测器采集投影数据，结合角谱衍射传播算法抑制了太赫兹波在样品外部的衍射效应，提高了成像速度[19]。

图4 快速三维太赫兹光谱层析实验示意图[31] Fig.4 Schematic of fast three-dimensional terahertz spectral tomography experiment[31]

图5 高功率连续波量子级联激光器的太赫兹CT成像结果[33] Fig.5 THz CT imaging results of high-power continuous wave quantum cascade laser[33]

随着太赫兹技术的发展，并且考虑到基于THz-TDS系统的成像系统较为复杂、成像速度慢，研究人员采用不同种类的太赫兹源进行了相关实验，并且太赫兹CT的应用研究也在不断开展。2006年，英国剑桥大学K.L.Nguyen等人首次采用频率为2.9 THz的高功率连续波量子级联激光器对聚苯乙烯泡沫样品及聚四氟乙烯样品进行吸收系数的太赫兹CT成像，实验结果如图5所示[33]。2012年，哈尔滨工业大学Q.Li等人采用频率为2.52 THz的连续太赫兹气体激光器对不同形状的聚苯乙烯泡沫进行成像，使用结合了高斯低通滤波运算的修正迭代算法以及数学形态学操作对数据进行处理，利用仅为FBP算法所需投影数的四分之一，即可获得高质量的3D成像，该方法在节省成像时间的同时保证了成像质量[17]。2013年，日本M.Suga等人引入了基于马赫曾德尔干涉仪(Mach-Zehnder interferometer)的相位调制技术，利用频率为540 GHz的连续倍频太赫兹源进行CT成像，通过获取透射光在每一点的相移信息实现了对样品折射率的重建，有效避免了吸收系数重建时由于折射率不匹配导致的折射、反射、散射等损耗而引起的图像边缘缺陷[34]。2014年，日本T.Kashiwagi等人采用由Bi2Sr2CaCu2O8+δ(Bi2212)本征约瑟夫森结构成的高温超导太赫兹辐射源进行了太赫兹CT成像。实验中采用频率为440 GHz的太赫兹波对内置玩偶的鸡蛋壳以及内含种子的心豆藤进行了CT成像。实验结果可以清晰显示出内部物体的位置与特征。这项研究有望促进太赫兹CT在生物医学中的应用[35]。2016年，日本名古屋大学S.R.Tripathi等人基于窄线宽、可调谐太赫兹参量源与太赫兹频率上转换探测技术对铅笔和内部有缺陷的塑料物体进行CT成像，由于该方法是将太赫兹频率的探测转换到近红外波段探测，在1.5 THz附近探测动态范围可达90 dB，实验中采用频率为1.5 THz的太赫兹CT成像可以很好地反映出物体内部信息及缺陷位置[23]。2017年，上海微系统与信息技术研究所周涛等人报道了基于单载流子光电二极管(Uni-Traveling-Carrier Photodiode,UTC-PD)产生的90～140 GHz低相干太赫兹波[36]以及德国P. Hillger报道了基于低成本的SiGe-HBT源产生的490 GHz太赫兹波进行的CT成像实验[37]，这些研究为太赫兹CT成像应用走向低成本、小型化奠定了基础。

在实际应用方面，太赫兹波由于具有安全性与穿透性，很适合用于研究文物以及材料的无损探伤。2011年，美国D.J.Roth等人对太空飞船外部使用的起隔热保护的塑料泡沫材料进行太赫兹CT成像，以探测其内部的缺陷及孔洞。以X射线CT图像为标准，实验结果表明太赫兹CT能够正确反映样品结构，如图6(a)所示。这展示了太赫兹CT在航天无损探伤应用的潜力[38]。同年，法国波尔多大学课题组采用耿氏二极管产生的110 GHz太赫兹波对古埃及陶器进行了太赫兹CT成像，如图6(b)所示，通过CT成像能够得出陶器内物质的大致形状、体积及吸收系数[39]。第二年，他们采用相同装置对干燥的人类腰椎骨进行CT成像，发现太赫兹CT可以区分硬质骨和松质骨，并取得了优良的成像结果[40]，如图6(c)所示。这些研究对于考古学与人类学具有重要的意义。

图6 太赫兹CT在无损检测方面的应用。(a)NASA隔热保护材料[38]；(b)古埃及陶器[39]；(c)人类腰椎骨[40] Fig.6 Applications of Terahertz computed tomography in nondestructive testing. (a)NASA thermal insulation system material[38]; (b)ancient Egyptian sealed pottery[39]; (c)human lumbar vertebra bone[40]

在上述实验中，大多采用FBP算法对图像进行重建，但该算法忽略了太赫兹光束本身的光学特性，太赫兹波在样品中的传播过程与X射线并不相同，具体表现为：(1)X射线不管是在表面还是内部都沿着直线传输，几乎没有反射与折射。但太赫兹波在空气及样品中传播时折射率有较大的差异，太赫兹波在样品表面的折射以及反射现象很严重，当其以不同入射角入射时，其反射率不同且由折射引起的偏转也不同，这对图像的重建和光束的探测提出了较高的要求；(2)太赫兹光束适用于高斯光束模型，其与样品的相互作用不能简单地用拉东变换理论模型中的线积分来描述；(3)由于太赫兹波波长较长，散射、衍射、边缘效应等现象会对成像结果造成一定影响。基于上述原因，在太赫兹CT的前期研究中发现，对复杂且折射率较高的样品(如火鸡骨、聚四氟乙烯)成像时，其内部结构并不能展现出来，如图2(b)、图5(b)所示。这都是由于边界的强烈反射与折射所造成的。由此可见，鉴于太赫兹波本身的光束特性，在对物体吸收系数和折射率重建时，太赫兹波探测器本身的限制和样品的折射、反射损耗等问题会在很大程度上影响成像质量。

图7 (a)干燥鸡腿骨照片(黑线表示断面区域)；(b)断面照片；(c)强度断面重建；(d)多峰值平均断面重建[44] Fig.7 (a)Digital photograph of the dried chicken femur (the horizontal black line indicates the region of the cross-section); (b)digital photograph of the sectioned chicken femur; (c)Strength section reconstruction; (d)cross-section reconstruction of multi-peak average[44]

图8 (a)聚合物颅骨模型照片，使用优化的(b)FBP；(c)SART；(d)OSEM算法重建的三维结构[46] Fig.8 (a)Photograph of the polymer skull, 3D visualization of the optimized (a)FBP; (b)SART; (c)OSEM results[46]

图9 (a)喷头照片；(b)FBP三维重建结果；(c)SART三维重建结果；(d)ML-TR三维重建结果[47] Fig.9 (a)Photograph of nozzle; (b)3D visualizations of FBP; (c)SART; (d)ML-TR results[47]

图10 (a)修正带照片；(b)内部机械构造；(c)塑料外表面[48] Fig.10 (a)Photograph of correction tape; (b)Internal mechanical structure(ruban(blue) and head of the roller(green); (c)external surface of plastics[48]

图11 使用(a)传统的算法；(b)修正的算法重建的天然软木塞二维断面图像[50] Fig.11 Reconstructed 2D tomographic images of the natural cork stopper (a)without and (b)with correction algorithm applied[50]

为了提高重构图像质量，研究人员主要在太赫兹光束特性变换、重建算法优化以及图像处理等方面进行了大量研究。在太赫兹光束特性变换方面，2011年，英国K.B.Ozanyan等人在THz-TDS系统信号光的数据采集前采用由两个小孔光阑组成的准直系统对太赫兹波进行空间低频滤波，用来滤除折射及散射的光[41]。该方法有效地避免了由于光束聚焦而引起的对样品尺寸的限制。2014年，以色列A.Bitman等人使用锥透镜将高斯光束变成宽带贝塞尔光束进行CT成像，由于贝塞尔光束具有零级光束无衍射、传播过程中分布图案不变的特性，相比于高斯光束可以获得更好的成像质量[42]。2014年，德国A.Brahm等人利用THz-TDS系统深入地研究了不同材料、不同形状样品进行CT投影时的光学效应，同时采用Zemax软件对投影进行理论模拟。结果表明，实验观察到的折射效应与理论计算能够很好的吻合。需特别指出的是，当太赫兹波的聚焦光斑与样品尺寸接近时，CT成像过程中的折射、衍射、边缘效应、噪声等因素带来的影响将会增加，传统的经典光学模拟将不再适用[43]，建议在理论上发展新的物理近似方法以处理衍射和脉冲在样品边缘处的分裂，实验中采用多通道探测器以克服折射效应。在重建算法优化方面，法国波尔多大学的研究人员进行了大量的工作，2010年，E.Abraham等人采用THz-TDS系统对高折射率样品特氟龙圆柱进行CT成像，利用多峰值平均的处理方法，有效降低了由于强烈折射导致的多峰值现象的影响。图7为他们采用该方法对干燥鸡腿骨的成像结果。对比可知，多峰值平均的方法可以更准确地反映样品结构和折射率的变化[44]。2011年，B.Recur等人研究了FBP、SART、OSEM这3种不同重建方法的优缺点，分别从强度、对比度以及几何复原情况等方面比较了不同的投影角度数对物体的断面重建图像质量的影响[45]。2012年，他们使用高斯光束模型对这3种CT重建算法进行优化，并采用结构相似性(structural similarity,SSIM)作为评价成像的参数，实验中基于耿氏二极管产生240 GHz的连续波对插有金属杆的聚合物颅骨模型进行层析成像。结果表明，优化的重建算法可以减少金属杆的图像散布，得到的边界体积也更接近实际情况，如图8所示[46]。2014年，B.Recur等人对X射线CT中用到的极大似然层析(Maximum Likelihood for TRansmission tomography,ML-TR)方法进行了研究，通过结合高斯光束模型，并考虑空背景扫描和暗视场的测量强度分布，对ML-TR中的有序凸子集算法(Ordered Subsets Convex,OSC)进行优化。在对喷头样品的CT三维重建中，这种算法相比于传统的FBP和SART算法，减少了噪声和条形伪影，可以更准确地反映样品内外部的结构，如图9所示[47]。但恢复图形中有一些空洞缺陷，推测是因为实验系统自身条件所引起的。随后，他们进一步提出了一系列太赫兹CT图像处理方法，包括采用OSC算法进行重建，采用K均值(K-means)算法对感兴趣区域进行自动分割，采用连通区域标记算法对组分进行标记，三维渲染，图像骨架化算法以及口径跟踪对图像进行细化处理，实现了高质量的三维成像。图10为他们对修正带样品的成像以及处理结果[48]。这些研究推动了太赫兹CT在工业无损探伤领域的应用。美国S.Mukherjee等人采用THz-TDS系统对圆柱形有机塑料棒进行了层析成像，根据菲涅尔公式理论计算了光束在样品表面的反射，同时使用BEAM4软件模拟了光束偏移导致的影响，进一步对重建算法进行了校正。将这种校正方法再次应用在软木塞上，相比他们之前采用的方法[49]，可以消除边界效应并能显示样品内部的缺陷结构，如图11所示[50]。2016年，德国J. Tepe等人报道了一种新型、高效的太赫兹CT重建算法，这种算法对代数重建技术(Algebraic Reconstruction Technique,ART)进行了修正，考虑了折射、反射对成像所带来的影响，基于光程差和透射损耗重建了复折射率。这种方法的优点是可以避免求解亥姆霍兹偏微分方程及波动方程，从而省去了大量计算，但并未考虑衍射效应及光束的高斯特性[51]。

3 反射式太赫兹层析

基于透射式的太赫兹CT技术并不适用于所有样品和环境，对于吸收系数高的样品，太赫兹波在样品内部具有一定的穿透深度却无法完全穿透样品，特别是在样品无法旋转等特殊情况下，太赫兹CT难以进行探测。反射式层析技术可以对太赫兹波穿透深度范围内的样品进行成像，有效地解决上述问题。相比于透射式的太赫兹CT技术，反射式太赫兹层析技术由于避免了样品旋转而具有更快的成像速度。

3.1 太赫兹飞行时间层析

当太赫兹脉冲入射到样品中后，会在不同深度发生反射，根据反射脉冲峰值的时间延迟可以得出折射率的深度分布信息，由此可以重建出样品折射率的三维分布。这种层析技术有3个前提假设：(1)样品的散射及衍射忽略不计；(2)样品内部的多次反射忽略不计；(3)每一层的折射率均匀不变。因此，该方法适用于结构层次分明的样品。特别是，由于飞行时间层析的纵向分辨率取决于太赫兹波的脉冲宽度，其通常为几十微米的量级，相比于透射式太赫兹层析，太赫兹飞行时间层析技术最大的优势是纵向分辨率极高。因此在水污染检测[52]、艺术品成像[53]、缺陷探测[54]、药片包层分析[55]等方面有着广泛的应用。飞行时间层析是最早在太赫兹层析中使用的技术[16]。2001年，美国J.L.Johnson等人将迈克尔逊干涉仪的结构引入飞行时间层析系统，通过聚焦产生的古依相移实现相干相消，去除背景信号，从而提高样品信号的对比度，将纵向分辨率提高到了12.5 μm[56]。该方法对于具有单层结构的样品能实现有效的测量。但对于多层结构，特别是层厚度在10 μm以下的许多工业产品或生物组织(如皮肤)，层析技术的纵向分辨率有待进一步提高。2009年，日本J.Takayanagi等人基于脉宽为17 fs的全光纤飞秒激光器泵浦DAST晶体，产生宽带太赫兹波。通过压窄泵浦脉冲宽度并结合高斯窗口的反卷积信号处理技术来获得单峰分布的太赫兹短脉冲，从而提高纵向分辨率。该系统成功地对叠加的三张纸及包含仅2 μm厚的GaAs薄层的半导体样品进行了层析成像[57]，其实验系统及对样品的层析成像结果如图12所示。为了实现样本的精确测量，不仅需要提高纵向分辨率，还要解决信号在多层结构内部的多次反射问题。2016年，韩国H.Park等人通过建立分层反射模型来估计上层反射信号，消除不同层间反射信号的干涉，提高了成像质量[58]。此外，对于背向反射信号难以分析的情况，可以采用反射式计算机层析的方法，通过旋转样品，得到样品的边界分布[59]。

图12 (a)基于光纤激光器的高分辨率太赫兹飞行时间层析实验装置；(b)三张纸的三维层析成像结果；(c)半导体样品示意图；(d)GaAs薄层层析成像图[57] Fig.12 (a)Experimental set up for the fiber-laser, high-resolution THz time-of-flight tomography system; (b)3D terahertz-tomography image of three sheets of paper; (c)schematic diagram of the semiconductor sample; (d)GaAs layer tomography[57]

在提高成像速度方面，2005年，日本大阪大学T.Yasui与法国波尔多大学E. Abraham等人采用结合线投影的非共线2D-FSEOS太赫兹成像技术，首次实现了对涂有醇酸漆膜的运动物体的实时层析成像，并且，在对湿润漆膜干燥过程的动态成像中，成功地观测到了漆膜缩水的过程[60]。2012年，韩国标准科学研究院D.S.Yee等人采用电控光采样(Electronically-Controlled Optical Sampling,ECOPS)技术对软盘进行飞行时间层析成像。这种技术能够以1 kHz的扫描频率测量太赫兹波形以获得轴向数据，在对横向面积为100 mm×100 mm的样品成像时，仅需80 s，极大地提高了成像速度。实验中分别使用太赫兹飞行时间层析以及超声波层析对含有人为内部缺陷的玻璃纤维增强聚合物(Glass-Fiber-Reinforced Polymer,GFRP)样品进行层析成像。相比于超声波层析成像，太赫兹飞行时间层析能更好地反映样品纵向结构及内部缺陷，在深度结构检测时由于其低损耗而更具优势。2015年，他们进一步采用扫描振镜与太赫兹f-θ聚焦透镜偏转太赫兹光束，取代移动样品的机械扫描，从而提高成像速度。实验中对体积为100 mm×100 mm×3 mm的GFRP样品进行包含200×200×725个体素的层析成像仅需40 s[61-62]。也有研究者采用压缩感知(Compressed Sensing,CS)的方法来提高成像速度，并实现了对生物样品的层析成像[63]。

3.2 太赫兹光学相干层析

光学相干层析成像(Optical Coherence Tomography,OCT)是在光学低相干干涉的基础上发展起来的一种光学成像技术。1991年，美国麻省理工学院D.Huang等人首次报道了光学相干层析技术[64]。他们基于光纤耦合的迈克尔逊干涉仪并使用超辐射发光二极管发出的波长为830 nm低相干性红外光对人体视网膜和冠状动脉壁进行活体成像，通过移动参考镜的位置来调节参考光的光程，从而控制样品不同深度的反射光与参考光干涉，通过干涉强度图像获取样品的深度信息。由于参考光与样品反射光的光程差在相干长度范围内才会发生干涉，因此其纵向分辨率取决于光的相干长度，在成像光源的光谱为高斯型的情况下，纵向分辨率可以表示为：

(9)

其中，λ0是激光光源的中心波长，Δλ是激光光源光谱带宽的半峰全宽(Full Width Half Maximum,FWHM)值，n是样品的折射率。

当前，OCT技术已经在生物组织成像领域如牙科、皮肤科、眼科等有很大的发展，但由于其通常采用800、1 000、1 300 nm等波长附近的近红外光作为光源，探测深度只有2～3 mm，应用受到极大的限制。考虑到太赫兹波具有良好穿透能力，利用太赫兹波进行OCT成像可以获得更高的穿透深度[65]。2012年，日本大阪大学T.Isogawa等人首次将太赫兹波应用在时域光学相干层析中，采用一个放大受激发射(Amplified Stimulated Emission,ASE)光源以及一个单载流子光电二极管产生中心频率为330 GHz，光谱带宽为100 GHz的宽带低相干太赫兹波，其纵向分辨率为1mm，对塑料材料的探测深度可达10 mm[66]，实验系统及其对塑料样品的成像结果如图13所示。2013年，Z.G.Zuo等人采用中压汞灯作为辐射源进行OCT成像，其输出光谱覆盖1～20 THz频率范围，具有极短的相干长度。实验中采用平滑处理、小波降噪以及Ge片对探测光谱进行滤波等方法提高信噪比，达到了100 μm以下的纵向分辨率[67]。

图13 (a)太赫兹时域光学相干层析实验装置；(b)样品空间分布示意图；(c)层析成像结果[66] Fig.13 (a)Experimental set up for THz time-domain optical coherent tomography; (b)schematic of the sample spatial distribution; (c)tomographic images of the objects[66]

上述太赫兹OCT技术均为时域OCT(Time Domain Optical Coherence Tomography,TDOCT)，属于第一代OCT技术，由于需要参考镜的机械移动，成像速度较慢。另一种OCT技术为傅里叶域OCT(Fourier Domain Optical Coherence Tomography,FDOCT)，FDOCT通过记录低相干干涉的光谱信号，利用傅立叶分析实现深度信息的并行获取，因此不需要参考镜的机械移动，在不降低信噪比的情况下可以实现高速成像[67]。根据干涉光谱信号的获取方式的不同，FDOCT大致有两种实现方法：谱域OCT(Spectral Domain Optical Coherence Tomography,SDOCT)和扫频OCT(Swept Source Optical Coherence Tomography,SSOCT)。当前报道的太赫兹FDOCT多采用SSOCT。

SSOCT使用波长扫描光源，设扫频源输出的功率谱为S(k)，探测到的干涉光谱条纹可表示为：

(10)

其中，Pr表示参考镜反射的光功率，P0为照射到样品上的光功率，r(z)和φ(z)分别表示样品反射系数的幅度和相位，Γ(z)为光源相干函数，k(t)表示波数随时间变化的函数，z表示样品深度坐标。公式(10)右边第一项为直流项，第二项为样品各层返回的光强和样品各层之间的干涉光强，称为自相关项，前两项代表不需要的背景和噪声。第三项为参考光与样品各层反射光之间的干涉光强，即扫频光学相干层析系统探测到的有效干涉光谱信号。根据维纳-辛钦定理，信号的自相关函数和功率谱密度是一对傅里叶变换对，因此对干涉光谱数据进行从波数空间到深度空间的傅立叶变换，即可恢复样品各层信号和参考镜信号的自相关信号，实现从波数空间到深度空间的转换，获得深度分辨的样品信息[68]。对干涉光谱信号进行傅立叶变换之后的表达式为：

FT-1[I(k)]=FT-1[S(k)]⊗{[δ(z)]+

(11)

其中，S(k)为光源功率谱，AC项为自相关项经过快速傅里叶变换(Fast Fourier Transformation,FFT)后的结果。扫频OCT的纵向分辨率与成像系统中所使用的激光光源的中心频率和光谱带宽有关，其理论值与传统时域OCT相同。

2013年，日本大阪大学T.Isogawa等人将太赫兹时域OCT进一步扩展到了太赫兹扫频OCT，利用一个频率为193.7～194 THz的扫频激光器以及一个194.2 THz的单频激光器，并采用UTC-PD模块光混频产生250～450 GHz的扫频太赫兹波，实现太赫兹扫频OCT[69]。利用SBD探测干涉光谱，纵向分辨率达到0.92 mm，且成像耗时仅为该小组之前报道的太赫兹时域OCT的一半[66]。第二年，他们改用193.4～193.8 THz的扫频激光器产生400～800 GHz的扫频太赫兹波，从而进一步地将纵向分辨率提高到了0.61 mm。此外，他们采用电子元件经过一次倍频以及两次三倍频后产生400～780 GHz太赫兹波进行扫频太赫兹OCT，纵向分辨率为0.84 mm。虽然采用电学方法的系统纵向分辨率不如采用光学方法的系统，但通过电学方法产生的太赫兹波功率更高，有利于系统的搭建和实际应用[70]。在提高成像速度方面，2016年，韩国标准科学研究院D.S.Yee等人采用波长为1 544～1 558 nm快速扫频源与波长为1 545 nm的分布反馈激光二极管(Distributed Feedback Laser Diode,DFB-LD)组成光拍频源，基于光混频原理产生快速扫频太赫兹信号，通过相干零差探测取代以往的自由空间干涉仪来实现纵向扫描，并采用振镜及太赫兹f-θ透镜实现横向快速扫描，当扫频速度为1 kHz时，仅需100s即可完成纵向5次平均、截面100×100像素的GFRP样品成像[71]。

3.3 太赫兹调频连续波雷达成像

太赫兹FMCW雷达成像是利用外差探测的原理测量光波传输延时从而得出样品距离或深度。光波由分束器分为信号光及参考光，信号光经样品反射回探测器，参考光直接入射探测器，信号光波和参考光波的频率均被调制而随时间变化，因此两束光在探测处混频后产生拍频信号fb，其满足下式[24]：

(12)

其中，fSR为扫描速率，Δt为信号光传输延时，Δf为频率调制范围，Ts为扫频时间，R为样品距离(深度)。因此样品距离可以由探测到的拍频fb得出：

(13)

且其纵向分辨率可表示为：

(14)

2009年，德国Synview公司基于电子学器件研发了一套结构简单、紧凑快速的太赫兹三维成像系统，当采用中心频率为300 GHz(带宽90 GHz)的太赫兹源时，空间分辨率为1 mm。可实现多个样品表面距离或多层间距大于1 mm的分辨。采用该系统对藏有手枪的假人的成像结果如图14所示[72]。由公式(14)可知，太赫兹FMCW技术的纵向分辨率由频率范围决定，与毫米波相比，其具有更高的分辨率，有望应用于高精度测距、相位灵敏探测、远距离安检等实际场合中[73-74]。

图14 (a)藏有手枪的假人；(b)相应的太赫兹图像[72] Fig.14 (a)Dummy with a pistol and its (b)corresponding THz image[72]

4 其它太赫兹三维成像技术

其它太赫兹三维成像技术如太赫兹断层合成层析(THz Tomosynthesis,THz-TS)、太赫兹衍射层析(THz Diffraction Tomography,THz-DT)、太赫兹菲涅尔透镜层析(THz Tomography with Fresnel Lens)等，适合于某些特殊情况下的成像。这些技术由于各种因素，在提出之后并没有太大的研究进展。还有一些在其它波段广泛使用的三维成像技术，但在太赫兹波段的研究较少，如太赫兹三维数字全息(THz 3D holography)、合成孔径雷达技术(Synthetic Aperture Radar,SAR)。这一节将对它们做简单介绍，不再赘述。

4.1 太赫兹断层合成层析

断层合成层析最早于20世纪30年代应用在X射线成像中[75]，是一种类似于计算机辅助层析的技术。不同的是，这种技术只需要几个投影角度就可以完成重建，其投影角度通常在-50°到50°之间[76]，成像速度较快。2009年，日本N.Sunaguchi等人采用频率为540 GHz的连续太赫兹固态倍频源对3个位于便签纸不同深度的字母“T”“H”“Z”进行TS成像。采用反投影的方法及维纳滤波对数据进行处理，能较好地重建3个字母的空间分布[77]。其成像原理如图15所示。这种技术可以看作一种丢失信息的CT技术，其分辨率不高但成像速度更快。相比于太赫兹CT，这种技术还有一个巨大的优势，正如上文所提到的，太赫兹CT受到严重的反射、折射影响，尤其是对于宽且扁平的板状样品，在边缘处入射角度很大，信号衰减严重，很难获得完整的投影数据，成像过程还会有许多冗余信息，而太赫兹TS技术由于投影时的入射角较小，不会受到折射、反射现象的影响。太赫兹TS技术可以作为太赫兹CT的一种替代手段，对于需要局部成像、快速成像或是无法完全旋转样品的情况，是一种十分有效的方法。

图15 断层合成层析实验原理 [77] Fig.15 Experimental principle of tomosynthesis[77]

4.2 太赫兹衍射层析

太赫兹衍射层析是利用太赫兹平面波射入样品后的衍射场分布来提取样品折射率的分布信息。根据傅立叶衍射理论，当平面波照射在样品上时，测量平面上前向散射场的傅立叶变换与样品函数的三维频域傅立叶变换在移位半球上的值成正比。由此可以将测量平面的分布与样品的分布函数联系起来，通过算法重建样品的三维分布。2004年，S.Wang等人首次采用基于THz-TDS系统的衍射层析对聚乙烯样品进行成像[78]，其实验系统如图16所示。相比于太赫兹CT，太赫兹衍射层析的成像速度更快，但由于缺乏有效的重建算法，它的成像质量较差。

图16 (a)太赫兹衍射层析实验装置；(b)样品；(c)样品重建图像[78] Fig.16 (a)Experimental setup of terahertz diffraction tomography; (b)sample; (c)reconstruction of the sample[78]

4.3 太赫兹菲涅尔透镜层析

菲涅尔透镜由法国物理学家A.Fresnel发明，其镜片一面为平面，另一面由从小到大的同心圆组成。得益于现在的微加工技术，大孔径、短焦距的菲涅尔透镜制造成为可能，其为太赫兹成像提供了光束质量较好的聚焦器件，且结构紧凑、重量轻。菲涅尔透镜的焦距具有与入射光波频率成正比的特性，当改变入射光波频率时，可以实现样品不同深度信息的成像。2003年，S.Wang等人首次进行了太赫兹菲涅尔透镜层析[79]，采用THz-TDS系统产生的宽带太赫兹脉冲透射样品，通过改变频率对空间不同位置分布的3个字母进行成像，如图17所示。由于太赫兹菲涅尔透镜层析的纵向分辨率受太赫兹频率测量分辨率、光束聚焦焦深(即两倍的瑞利长度)、非傍轴光束的焦散曲线等因素的限制，其分辨率较差，为5 mm左右。此外，由于菲涅尔透镜本身的衍射特性，基于菲涅尔透镜的三维成像技术对吸收系数的分布比对折射率的分布更为敏感，且不包含光谱信息。

图17 对“T”“H”“Z”字母掩模板的菲涅尔透镜层析成像[78] Fig.17 Fresnel lens tomographic imaging using ‘T’, ‘H’ and ‘Z’ masks[78]

4.4 太赫兹三维数字全息

数字全息技术是随着计算机的迅猛发展而出现的，通常使用CCD相机记录光的相位及强度，再经过计算机重建样品信息。由于太赫兹波的波长较长，其成像质量的提高受限于衍射。而太赫兹全息成像不仅可以记录每个频率下对应的相位与振幅信息，还包含其他的散射信息。因此可以通过对光束传播行为的数值计算获得任意平面上的光场分布，从而消除了衍射的影响[80]。典型的太赫兹三维数字全息实验系统如图18所示[78]。

图18 太赫兹三维数字全息实验装置[78] Fig.18 Experimental setup for three-dimensional THz digital holography[78]

4.5 合成孔径雷达

合成孔径雷达是一种二维成像系统，其发射线性调频信号实现距离向的高分辨率，并利用回波的多谱勒相位变化合成等效的大孔径天线来获取高的方位分辨率。2010年，R.Heremans等人将合成孔径雷达成像技术应用在太赫兹波段[81]。目前合成孔径雷达正经历着从二维成像向三维成像的发展，采用干涉合成孔径雷达(Interferometric Synthetic Aperture Radar,InSAR)等方法可以实现对目标的三维成像，以合成孔径雷达复数据提取的相位信息为信息源来获取三维信息[82]。

5 结束语

本文介绍了基于透射式以及反射式的几种太赫兹波三维成像技术。每种技术适用于不同的实验环境及样品，且成像时间与分辨率也各有不同，在实际应用中，应该根据具体的要求选择合适的三维成像技术。总体而言，透射式太赫兹CT能够获取样品的光谱信息实现光谱层析，相对成熟的算法与技术保证了其优良的成像质量，且系统结构简单适用范围广，因此成为最受关注的太赫兹波三维成像技术。尽管已有公司开发出了实用化的成像仪器，但目前的太赫兹CT仍处于研究发展阶段，有许多问题亟待解决：(1)缺乏专门针对太赫兹CT的重建算法；(2)成像速度慢；(3)太赫兹光束的反射、折射以及散射现象降低了成像质量，尤其是太赫兹波在透射样品后会产生严重的光束偏移，难以应用于高折射率样品，这是限制太赫兹CT发展的重要因素。反射式层析技术具有更高的纵向分辨率，基于THz-TDS系统的飞行时间层析是最早使用的太赫兹层析技术，具有极高的纵向分辨率，但其系统复杂且昂贵，难以实用化。基于连续波的太赫兹OCT是近年来兴起的太赫兹层析成像技术，其系统结构简单、紧凑且便宜，具有广泛的应用价值。但其横向分辨率受到衍射极限的限制，纵向分辨率与光谱范围成正比。因此，反射式层析技术更深入的应用与发展依赖于宽光谱，需要继续对短脉冲太赫兹源以及用于提高信噪比并消除光束在样品内部多重反射现象的更加有效的信号处理算法进行研究。反射式层析成像技术适用于难以透射的样品，其与透射式太赫兹CT技术相互补充，使得太赫兹三维成像技术展现出了具有巨大的科学研究价值和广阔的应用前景。未来的太赫兹三维成像技术发展受到太赫兹源、探测器及相关功能器件的影响，包括高功率、可调谐、稳定性好的太赫兹源，高灵敏度、动态范围大的太赫兹波探测器以及损耗低的波导、偏振、聚焦等器件。随着成像系统的优化改进以及重建算法的突破性进展，太赫兹三维成像技术将朝着高速、高分辨率、小型化、实用化的研究方向发展。