太赫茲時域脈沖層析成像的Matlab實(shí)現(xiàn)

唐粹偉, 丁芳媚, 郭田田，孫定仟, 李秋勇，周 俊

(電子科技大學(xué) 物理電子學(xué)院，四川 成都 610054)

太赫茲時域脈沖層析成像的Matlab實(shí)現(xiàn)

唐粹偉, 丁芳媚, 郭田田，孫定仟, 李秋勇，周 俊

(電子科技大學(xué) 物理電子學(xué)院，四川 成都 610054)

該文在透射和反射成像模式下，實(shí)現(xiàn)對物體二維成像的關(guān)鍵算法和技術(shù)研究。根據(jù)太赫茲成像儀采集到的時域太赫茲信號，利用Matlab編程實(shí)現(xiàn)二維成像的關(guān)鍵算法，實(shí)現(xiàn)對原有功能的補(bǔ)充、完善和擴(kuò)展。在實(shí)際的復(fù)雜應(yīng)用中，對復(fù)雜樣品進(jìn)行成像實(shí)驗(yàn)，并用所實(shí)現(xiàn)的層析成像技術(shù)進(jìn)行圖像重現(xiàn)，與實(shí)物進(jìn)行對比，驗(yàn)證程序的正確性和實(shí)用性。

太赫茲；層析；成像；Matlab

太赫茲波和其他波段的電磁輻射一樣可以用來對物體成像，而且根據(jù)太赫茲波的高透性、無損性以及大多物質(zhì)在太赫茲波段都有指紋譜等特性[1-2]，使太赫茲成像相比其他成像方式更具優(yōu)勢[3-5]。

太赫茲層析成像技術(shù)可以獲得樣本的二維橫截面分布圖像。目前對太赫茲層析成像技術(shù)的研究主要集中在重建圖像質(zhì)量的提高、成像時間的縮短以實(shí)現(xiàn)實(shí)時成像等方面[6-8]。目前在太赫茲層析成像中對投影數(shù)據(jù)的獲取主要通過點(diǎn)掃描方式獲得，需要花費(fèi)大量時間進(jìn)行數(shù)據(jù)測量[9-10]。通過選擇合適的圖像重建方法可以顯著減少所需投影方向數(shù)，從而明顯縮短太赫茲層析成像所花費(fèi)時間，具有很重要的意義[11,12]。

1 成像原理

太赫茲掃描成像實(shí)驗(yàn)主要是借助英國Teraview公司的TPS 3000型太赫茲光譜儀進(jìn)行的，成像模式包括透射模式和反射模式。反射式太赫茲脈沖成像系統(tǒng)的光路示意圖如圖1所示。系統(tǒng)主要包括飛秒激光器、太赫茲波產(chǎn)生部分、太赫茲波探測部分、延遲裝置及樣品二維平移裝置5個部分。

利用THz 源所產(chǎn)生的已知THz 波作為成像射線，經(jīng)目標(biāo)透射或反射后記錄了其譜信息(包括振幅、相位及復(fù)介電常數(shù)的空間分布信息)，這些信息經(jīng)過轉(zhuǎn)換得到相應(yīng)電信號響應(yīng)，隨后經(jīng)適當(dāng)?shù)男盘柗治鎏幚?，最終得到目標(biāo)的二維或三維圖像。

圖1 反射型太赫茲脈沖成像系統(tǒng)示意圖

太赫茲時域脈沖掃描成像的優(yōu)勢在于：可以通過對不同時間延遲的掃描, 在x-y平面上每一空間點(diǎn)上獲得帶有樣品信息的透射或反射時域波形(數(shù)據(jù)格式為三維格式:x×y×t)，其中t為時間軸。

2 各種成像顯示模式及Matlab實(shí)現(xiàn)

成像顯示模式(view mode)是指在計(jì)算機(jī)軟件中將采集的時域數(shù)據(jù)(由單點(diǎn)THz 波形組成)在整個C-Scan窗口呈現(xiàn)出來的形式。不同的顯示模式下，凸顯出樣品信息的重點(diǎn)也不一樣。

2.1 時間切片/頻域切片

從成像數(shù)據(jù)文件中提取出來的時域信號為512×N的矩陣，矩陣每一列的數(shù)據(jù)是每個像素點(diǎn)的時域信號離散值(512個)，而N代表有N個像素點(diǎn)。時域切片成像模式(time slice)是在每一個像素點(diǎn)選取同一時間延遲位置的信息幅值的大小，進(jìn)行二維圖像顯示。

頻域切片成像模式(freq slice)，是對頻域數(shù)據(jù)而言的。對于每一個像素點(diǎn)的時域數(shù)據(jù)，對其進(jìn)行傅里葉變換即可得到該點(diǎn)采集到的信號的頻域數(shù)據(jù)。同樣，在每一個像素點(diǎn)選取同一頻域位置的信號幅值進(jìn)行二維顯示。

2.2 最大/最小飛行時間

對于x-y平面上的每一個像素點(diǎn)，此模式(time of flight(max))下會呈現(xiàn)最大信號所對應(yīng)的延遲時間大小的圖像。最大信號值處的時間坐標(biāo)就是最大飛行時間。該模式下圖像中每個像素點(diǎn)顯示值即為該時間坐標(biāo)值。通常情況下，由掃描物體前表面反射的信號最強(qiáng)因而會被認(rèn)為是最大值，樣品中表面最高的像素點(diǎn)處信號最先被反射回，因此飛行時間最小；而表面最低的像素點(diǎn)處返回最大值會需要更長的時間。如果表面平坦，這時最強(qiáng)信號對應(yīng)的激光飛行時間近似為常值。這時，圖像為一個無明顯變化的圖樣。

2.3 最大/最小峰值

最大/最小峰值(Max/Min peak size)模式下會呈現(xiàn)每個像素點(diǎn)在單點(diǎn)波形中的最大/小信號值。選取出每個像素點(diǎn)時域信號的最大/小值，即可得到1×N(N為像素點(diǎn)個數(shù))的序列，將此數(shù)據(jù)用二維圖展示出來。一般說來，在吸收度小的像素點(diǎn)上太赫茲波衰減較小，因而，此類像素點(diǎn)中最大峰值較其他像素點(diǎn)較大，反之則較小。

2.4 峰峰值切片

峰峰值切片模式(peak-to-peak in slice)，在每一個像素點(diǎn)都選取相同的時間區(qū)間，將區(qū)間內(nèi)的最大、最小時間信號提取出來相減得到該切片內(nèi)的差值(峰峰值)。將N個像素點(diǎn)的峰峰值結(jié)果呈現(xiàn)在x-y二維圖像中即是峰峰值切片成像效果圖。這種模式下，可以去除空氣/物體界面上的強(qiáng)反射波，而只關(guān)心樣品內(nèi)部的信息。

2.5 峰峰值

峰峰值模式(peak-to-peak)，是在每個像素點(diǎn)展示整個時域波形上最大值與最小值的差值。在每一個像素點(diǎn)都提取出其時域信號的最大、最小值相減得到該峰峰值。將N個像素點(diǎn)的峰峰值結(jié)果呈現(xiàn)在x-y二維圖像中即是峰峰值成像效果圖。因此，在這個模式下每個像素點(diǎn)都展示其最大值與最小值之間的差值。這種模式下，可以清晰展現(xiàn)樣品的輪廓。

2.6 微分/積分

一階微分信號表現(xiàn)每個時域信號在不同時間點(diǎn)變化的大小，即其變化率，能突出表現(xiàn)信號變化較大的部分。該模式(1st deriv time slice)下選擇一個時間點(diǎn)，計(jì)算出每個像素點(diǎn)該時間位置下信號的一階微分，微分值近似于該時間點(diǎn)和下一個時間點(diǎn)間的差值。微分信號表現(xiàn)每個時域信號在不同時間點(diǎn)變化的大小，即其變化率，能突出表現(xiàn)在某一時間切片點(diǎn)信號變化較大的部分。

Matlab實(shí)現(xiàn)的具體步驟是，首先選取兩個相鄰時間切片點(diǎn)的時域數(shù)據(jù)，讓兩個時間點(diǎn)的時域信號幅值相減，即得到所要的一階微分信號。之后，將得到的一階微分信號進(jìn)行二維顯示即可。

積分切片模式(integrated slice)，在每一個像素點(diǎn)選取相同的時間區(qū)間，將區(qū)間內(nèi)的所有時間信號相加得到積分切片值。將1×N個積分結(jié)果呈現(xiàn)在x-y二維圖像中即是積分成像效果圖。該模式相當(dāng)于將樣品物理剖分成小切片，然后對小切片進(jìn)行成像。這種成像模式的優(yōu)勢在于能在不對樣品進(jìn)行解剖切片的前提下，依次觀察樣品某一深度某一部分的狀態(tài)。

2.7 時域/頻域差分切片

對于每一個像素點(diǎn)的時域信號，所選時間區(qū)間的結(jié)束時間位置的信號值減掉起始時間位置處的信號幅值即為所需差值(diff time slice)，將每個像素點(diǎn)處的不同差值組合起來呈現(xiàn)于x-y平面得到二維圖像。Matlab實(shí)現(xiàn)的具體步驟是，首先選取任意兩個時間切片點(diǎn)(自己指定)的時域數(shù)據(jù)，讓兩個時間點(diǎn)的時域信號幅值相減，即得到所要的時域差分信號。之后，將得到的時域差分信號進(jìn)行二維顯示即可。

同樣，對于每一個像素點(diǎn)的頻域信號，所選頻率區(qū)間的結(jié)束頻率位置的信號值減去起始頻率位置處的信號幅值即為所需差值(diff freq slice)，將每個像素點(diǎn)處的不同差值組合起來呈現(xiàn)于x-y平面得到二維圖像。Matlab實(shí)現(xiàn)的具體步驟是，首先選取任意兩個頻率切片點(diǎn)(自己指定)的頻域數(shù)據(jù)，讓兩個時間點(diǎn)的頻域信號幅值相減，即得到所要的頻域差分信號。之后，將得到的頻域差分信號進(jìn)行二維顯示即可。

3 俯視圖和剖面圖成像效果

不同的模式下的數(shù)據(jù)進(jìn)行進(jìn)一步的插值函數(shù)的處理后，即可得到所需的俯視圖和剖面圖的成像效果圖。除此之外，對比度的調(diào)節(jié)功能可以手動調(diào)節(jié)，得到最好最清楚的成像效果。

俯視圖就是x-y切面的圖像。所測的封轉(zhuǎn)芯片的最小飛行時間的效果圖如圖2所示。

圖2 最小飛行時間成像效果圖

剖面圖可以手動選擇x-t平面或y-t的切面圖像。x-t切面成像效果圖如圖3所示，y-t切面成像效果圖如圖4所示。

圖3 水平切面成像效果圖

圖4 豎直切面成像效果圖

4 GUI設(shè)計(jì)

利用Matlab自帶的GUI界面設(shè)計(jì)功能，結(jié)合前文所述的不同顯示模式的原理，加入了時間段選擇控件、選擇成像點(diǎn)位置控件和控制成像圖像對比度的控件，并完成不同控件之間的關(guān)聯(lián)，完成軟件的設(shè)計(jì)，結(jié)果如圖5所示。

圖5 GUI界面圖

5 結(jié)束語

本文提取每一個數(shù)據(jù)點(diǎn)的時域波形，首先在時域上通過Matlab編程實(shí)現(xiàn)各種成像功能，包括幅值成像、最大/最小值成像、峰峰值成像、最長/最短飛行時間成像和時間片成像等等多種成像顯示模式。然后，對時域波形進(jìn)行傅里葉變換，得到頻域波形，在頻域上通過Matlab編程實(shí)現(xiàn)頻域幅值成像。

但是，整個程序中還是缺少對實(shí)驗(yàn)中數(shù)據(jù)進(jìn)行函數(shù)處理的功能，如濾波器處理、幅值成像、相位成像和單頻率點(diǎn)成像等，還需繼續(xù)加強(qiáng)。

太赫茲脈沖成像技術(shù)除了普通成像應(yīng)用，更可用于獲取物質(zhì)光譜、折射率、幅值、相位及樣品厚度等信息。太赫茲脈沖成像技術(shù)在無損探測、安檢、材料分析、醫(yī)療等方面用途廣泛，是未來國防、科研工作中需要提前占領(lǐng)的制高點(diǎn)。

[1] MITTLEMAN D M, GUPTA M, NEELAMANI R, et al.Recent advances in terahertz imaging[J].Applied Physics B, 1999, 68(6):1085-1094.

[2] HUNSCHE S, KOCH M, BRENER I, et al.THz near-field imaging[J].Optics Communications, 1997, 150(16):22-26.

[3] Chen Y, Huang S, Pickwell-MacPherson E.Frequency-wavelet domain deconvolution for terahertz reflection imaging and spectroscopy[J].Optics Express, 2010, 18(2): 1177-1190.

[4]Guo T T, DING Z L.Study on algorithms of terahertz pulsed reflection imaging[J].Millimeter Waves & Thz Technology Workshop,2013(1):1-2.

[5] 張振偉, 崔偉麗, 張 巖,等.太赫茲成像技術(shù)的實(shí)驗(yàn)研究[J].紅外與毫米波學(xué)報(bào), 2006(3)：217-220.

[6] SIEBERT K J, QUAST H, LEONHARDT R, et al.Continuous-wave all-optoelectronic terahertz imaging[J].Applied Physics Letters, 2002, 80(16):3003-3005.

[7] MICHAN S P, ABBOTT D, MUNCH J, et al.Analysis of system trade-offs for terahertz imaging[C]// Asia Pacific Symposium on Microelectronics and MEMS.International Society for Optics and Photonics, 1999:503-514.

[8] GASS J, DAKOFF A, KIM M K.Phase imaging without 2π ambiguity by multiwavelength digital holography[J].Optics letters, 2003, 28(13): 1141-1143.

[9] ZHANG L, ZHONG H, ZHANG Y, et al.Terahertz wave focal-plane multiwavelength phase imaging[J].JOSA A, 2009, 26(5): 1187-1190.

[10] FEDERICI J F, GARY D, BARAT R, et al.Terahertz imaging using an interferometric array[C]// Defense and Security.International Society for Optics and Photonics, 2005:2477-2479.

[11] HOSHINA H, NAKAJIMA S, YAMSHITA M, et al.Terahertz imaging diagnostics of the cancer tissues with chemometrics technique[C]// Infrared Millimeter Waves and 14th International Conference on Teraherz Electronics, 2006:195-195.

[12] KNOBLOCH P, SCHMALSTIEG K, KOCH M, et al.THz imaging of histo-pathological samples[C]// European Conference on Biomedical Optics.International Society for Optics and Photonics, 2001:239-245.

Implementation of Terahertz Time-Domain Tomography Using Matlab

TANG Cuiwei, DING Fangmei, GUO Tiantian, SUN Dingqian, LI Qiuyong, ZHOU Jun

(School of Physical Electronic,University of Electronic Science and Technology of China， Chengdu 610054，China)

In the transmission and reflection imaging modes, we realize the key algorithms and techniques of 2D imaging of objects.According to time-domain Terahertz signals collected by THz imaging instrument, we use MATLAB programming to realize key algorithm for 2D imaging to accomplish the improvement and expansion added to the original function.In the complex practical applications, imaging experiments are performed on complex samples, and the image is reconstructed by using the tomography technology.Compared with the real,the correctness and practicability of the program are verified.

terahertz; tomography; imaging; Matlab

2015-05-14；修改日期:2016-12-02

中央高?；究蒲袠I(yè)務(wù)費(fèi)(ZYGX2015J040)。

唐粹偉(1993-)，女，本科，電子信息科學(xué)與技術(shù)專業(yè)。

周俊，副教授，zhoujun123@uestc.edu.cn

TN722.5+8

A

10.3969/j.issn.1672-4550.2017.01.041