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

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

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

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

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

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

太赫茲；層析；成像；Matlab

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

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

1 成像原理

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

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

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

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

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

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

2.1 時(shí)間切片/頻域切片

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

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

2.2 最大/最小飛行時(shí)間

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

2.3 最大/最小峰值

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

2.4 峰峰值切片

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

2.5 峰峰值

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

2.6 微分/積分

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

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

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

2.7 時(shí)域/頻域差分切片

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

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

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

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

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

圖2 最小飛行時(shí)間成像效果圖

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

圖3 水平切面成像效果圖

圖4 豎直切面成像效果圖

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

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

圖5 GUI界面圖

5 結(jié)束語(yǔ)

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

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

太赫茲脈沖成像技術(shù)除了普通成像應(yīng)用，更可用于獲取物質(zhì)光譜、折射率、幅值、相位及樣品厚度等信息。太赫茲脈沖成像技術(shù)在無(wú)損探測(cè)、安檢、材料分析、醫(yī)療等方面用途廣泛，是未來(lái)國(guó)防、科研工作中需要提前占領(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ù)專(zhuān)業(yè)。

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

TN722.5+8

A

10.3969/j.issn.1672-4550.2017.01.041