郭佑?xùn)|，凌福日，姚建銓，3

（1.華中科技大學(xué)武漢光電國(guó)家研究中心，武漢430074；2.華中科技大學(xué)光學(xué)與電子信息學(xué)院，武漢430074；3.天津大學(xué)精密儀器與光電子工程學(xué)院，天津300072）

引 言

由于太赫茲波具有良好的穿透性、安全性等優(yōu)良特性，太赫茲成像在醫(yī)學(xué)影像、安全檢查、環(huán)境監(jiān)測(cè)等方面具有廣泛的應(yīng)用前景［1-6］。然而在實(shí)際應(yīng)用中，受到成像信號(hào)頻率、成像環(huán)境、探測(cè)器、器件制造工藝等因素的影響，利用實(shí)際太赫茲成像系統(tǒng)采集的太赫茲圖像往往存在信噪比低，邊緣模糊嚴(yán)重等現(xiàn)象［7］，這使得太赫茲原始圖像的質(zhì)量無法滿足人們的正常視覺要求，以及研究人員對(duì)圖像進(jìn)行進(jìn)一步研究的需要。提高太赫茲圖像質(zhì)量最直接的方法就是升級(jí)系統(tǒng)硬件，但該方法不僅有技術(shù)難度，還面臨成本約束［8］。因此，研究進(jìn)一步改善已有太赫茲圖像質(zhì)量的方法，能夠在一定程度上彌補(bǔ)現(xiàn)有裝置的不足。在不升級(jí)成像系統(tǒng)的前提下，超分辨率（super-resolution，SR）重建技術(shù)提供了一種從低分辨率（low-resolution，LR）圖像恢復(fù)出相同場(chǎng)景的高分辨率（high-resolution，HR）圖像的實(shí)用、有效的低成本解決方案［9-10］。該技術(shù)可以增加現(xiàn)有圖像的空間分辨率，從而克服傳感器和光學(xué)成像的固有分辨率限制［11-12］。近年來，研究人員持續(xù)開展了利用超分辨率重建技術(shù)對(duì)太赫茲圖像進(jìn)行重建的研究，包括基于多幀圖像的凸集投影（projections onto convex sets，POCS）算法、迭代反向投影算法、基于單幀的Lucy-Richardson迭代、形態(tài)學(xué)算法、反卷積方法［13-16］等等。

SR技術(shù)是一個(gè)病態(tài)逆問題，與多幀圖像相比，單幀圖像僅含有少量信息，病態(tài)重建問題更為復(fù)雜。其核心在于建立從LR圖像到HR圖像的非線性映射，由于LR圖像信息不充足，因此，需要根據(jù)先驗(yàn)信息建立正則條件來對(duì)重建過程進(jìn)行約束［12］?；谔荻茸儞Q的圖像超分辨率重建方法，主要思想是根據(jù)圖像的梯度先驗(yàn)信息對(duì)重建過程進(jìn)行約束來恢復(fù)圖像。由于人類視覺對(duì)圖像邊緣的敏感性，2007年，F(xiàn)ATTAL等人首次提出利用圖像梯度輪廓的先驗(yàn)信息來描述圖像邊緣的空間分辨率［17］。2011年，SUN等人利用通用高斯函數(shù)對(duì)邊緣梯度輪廓進(jìn)行建模，從而首次提出了基于梯度變換（gradient field transformation，GFT）的圖像超分辨率重建的方法［18］。2015年，NGUYEN等人利用三角模型以及混合高斯模型（triangle and Gaussian mixed model，TGM）來改進(jìn)邊緣擬合函數(shù)以描述圖像邊緣梯度輪廓的非對(duì)稱性，但是該模型無法有效擬合較為復(fù)雜的邊緣［19］。2018年，QIANG等人進(jìn)一步提出自適應(yīng)梯度變換方法，使用雙邊高斯模型擬合圖像邊緣梯度輪廓，同時(shí)利用“差函數(shù)”定義梯度變換來優(yōu)化低分辨率圖像的梯度圖像，并以此作為正則項(xiàng)約束LR圖像與HR圖像之間的關(guān)系，進(jìn)行超分辨率重建，銳化圖像邊緣［20］。由于太赫茲信號(hào)的波長(zhǎng)屬于亞毫米級(jí)，相較于可見光圖像，太赫茲圖像的邊緣較為模糊。因此，基于梯度變換的重建方法為這一問題的改善提供了一種合適的解決方案。

本文中基于實(shí)際太赫茲透射式成像系統(tǒng)，提出利用有理分形插值結(jié)合基于梯度變換的超分辨率重建算法對(duì)該系統(tǒng)所采集的單幀太赫茲圖像進(jìn)行超分辨率重建，驗(yàn)證并擴(kuò)展了基于梯度變換的重建方法的適用性。針對(duì)所采集的太赫茲圖像中存在的對(duì)比度較低和空間噪聲較大等影響其梯度圖像質(zhì)量的因素，引入基于空間信息熵的直方圖匹配技術(shù)和雙邊濾波器減弱這些不良因素對(duì)梯度圖像質(zhì)量的影響，提高了重建算法的性能。經(jīng)過重建之后，太赫茲圖像的細(xì)節(jié)信息得到恢復(fù)，邊緣得到銳化，圖像質(zhì)量得到提高；同時(shí)當(dāng)成像信號(hào)頻率較低時(shí)，該方法對(duì)圖像質(zhì)量的提升效果更加明顯。

1 超分辨率重建模型

1.1 觀測(cè)模型

超分辨率圖像重建的觀測(cè)模型如圖1所示。圖1中，PSF指點(diǎn)擴(kuò)散函數(shù)（point spread function）。由于設(shè)備和環(huán)境的原因，成像系統(tǒng)獲得的圖像通常是高分辨率圖像經(jīng)過幾何變形、模糊、降采樣等退化過程之后形成的低分辨率圖像［21］，該過程可以由下式表示：

式中，*表示卷積，ILR為L(zhǎng)R圖像，IHR為HR圖像，D為退化矩陣，↓為降采樣，n為成像過程中的各類噪聲。

Fig.1 Observation model of super-resolution

超分辨率圖像重建是圖1所示圖像退化過程的逆過程。該過程是一個(gè)不適定問題，需要引入外部約束和先驗(yàn)信息進(jìn)行優(yōu)化求解。由LR圖像估計(jì)HR圖像的表達(dá)式為：

式中，IHR為重建后的圖像，U（IHR）為正則項(xiàng)，參量β為加權(quán)平衡參量。減小β將減少正則項(xiàng)對(duì)圖像重建的作用，當(dāng)β=0時(shí)，該式退化為最大似然估計(jì)?；谔荻茸儞Q的重建方法將理想圖像的梯度先驗(yàn)信息作為正則項(xiàng)（即（2）式中的 U（IHR））來約束（2）式，達(dá)到銳化圖像邊緣的目的。

1.2 基于梯度變換的超分辨率重建

基于梯度變換的超分辨率重建方法框架如圖2所示。該方法首先對(duì)低分辨率太赫茲圖像進(jìn)行插值，初步提高圖像質(zhì)量，消除鋸齒狀邊緣。插值后的太赫茲圖像（Iu）的細(xì)節(jié)信息包含在其梯度圖像（▽Iu）當(dāng)中。然后提取Iu的梯度圖像，通過梯度變換使該梯度圖像更加尖銳，變換后的梯度圖像（▽I^HR）將包含更多的細(xì)節(jié)信息。利用▽I^HR作為正則項(xiàng)并根據(jù)最大后驗(yàn)概率

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

（maximum a posteriori，MAP）方法約束圖像的重建過程可以獲得高分辨率圖像（IHR），達(dá)到銳化圖像邊緣的目的。▽越尖銳，IHR的邊緣也越尖銳，同時(shí)也將包含更多的細(xì)節(jié)信息。受到太赫茲圖像對(duì)比度較低以及空間噪聲的影響，▽Iu的質(zhì)量無法滿足有效進(jìn)行梯度變換的要求。本文中引入基于空間信息熵的直方圖匹配技術(shù)和雙邊濾波器（bilateral filter，BF）對(duì)Iu進(jìn)行處理，優(yōu)化其梯度圖像。

1.2.1 圖像插值 受到成像系統(tǒng)效率的限制，太赫茲原始圖像中的像素點(diǎn)較少，邊緣不連續(xù)。從太赫茲原始圖像中提取的梯度圖像質(zhì)量較差，無法對(duì)該梯度圖像進(jìn)行梯度變換。因此，本文中首先利用有理分形插值對(duì)原始圖像進(jìn)行處理，采用具有較高插值精度的有理插值函數(shù)和能夠較好地描述圖像紋理細(xì)節(jié)的分形插值方法，擴(kuò)大圖像規(guī)模的同時(shí)，能夠更好地保持圖像邊緣，從而提高梯度圖像的質(zhì)量。分形插值通過一個(gè)特定的迭代函數(shù)系統(tǒng)生成分形插值函數(shù)而實(shí)現(xiàn)，假設(shè)原始圖像所在的平面區(qū)域?yàn)?Ω=I×J=［a，b］×［c，d］，其中I，J分別表示原始圖像水平方向和垂直方向像素坐標(biāo)的集合，原始圖像的像素點(diǎn)為｛（xi，yj，zi，j），i=1，…，N；j=1，…，M｝，記區(qū)間 Ii=［xi，xi+1］，Ji=［yj，yj+1］，其中 i∈Γ=｛1，…，N-1｝，分形插值函數(shù)的表達(dá)式為：

式中，φi（x）：I→Ii，φj（y）：J→Jj為插值區(qū)間映射，F(xiàn)i，j：Ω→R為各插值點(diǎn)的灰度值映射，R為插值后圖像的區(qū)域，si，j為尺度因子。有理分形插值將（3）式中的qi，j（x，y）定義為：

式中，Pi，j（φi（x），φj（y））和 Bi，j（x，y）分別為雙變量有理插值函數(shù)和有理擾動(dòng)基函數(shù)［22］。

1.2.2 圖像增強(qiáng) 受到成像信號(hào)強(qiáng)度的影響，太赫茲原始圖像的對(duì)比度較低，從插值后的太赫茲圖像中提取的梯度圖像質(zhì)量較低，對(duì)該梯度圖像直接進(jìn)行梯度變換將無法得到理想的結(jié)果。因此，本文中利用基于空間信息熵的直方圖匹配技術(shù)對(duì)插值后的太赫茲圖像進(jìn)行增強(qiáng)，提高對(duì)比度，優(yōu)化梯度圖像?？紤]到圖像灰度級(jí)的空間分布，將圖像劃分為多個(gè)子區(qū)域，并在每個(gè)子區(qū)域根據(jù)下式進(jìn)行圖像增強(qiáng)處理，可以實(shí)現(xiàn)圖像的全局增強(qiáng)。

式中，k為子區(qū)域內(nèi)的灰度級(jí)，Ik，e為利用（5）式對(duì)灰度級(jí)k進(jìn)行映射后在經(jīng)過增強(qiáng)的圖像中的灰度級(jí)，F(xiàn)k為灰度級(jí)映射函數(shù)。fl為一個(gè)關(guān)于l的離散函數(shù)，根據(jù)圖像的空間信息熵計(jì)算［23］。Imin和Imax分別為輸入圖像最小和最大的灰度值。經(jīng)過全局增強(qiáng)之后，再對(duì)圖像進(jìn)行離散余弦變換并對(duì)其變換系數(shù)根據(jù)下式進(jìn)行加權(quán)實(shí)現(xiàn)局部增強(qiáng)，進(jìn)一步提升圖像的對(duì)比度。

式中，K為子區(qū)域內(nèi)的灰度級(jí)數(shù)量，參量γ用于平衡全局和局部增強(qiáng)的程度，其取值范圍為，本文中該參量取0.25以保證圖像的增強(qiáng)效果，同時(shí)不會(huì)由于過度增強(qiáng)圖像而放大其中的噪聲。

1.2.3 雙邊濾波 由于成像過程中太赫茲信號(hào)存在衰減且不同采樣點(diǎn)采集到的信號(hào)衰減程度不同，同時(shí)受到成像系統(tǒng)其它噪聲（如探測(cè)器隨機(jī)噪聲、掃描系統(tǒng)電磁噪聲等）的影響，太赫茲圖像中的空間噪聲將對(duì)梯度圖像的質(zhì)量產(chǎn)生較大影響。因此，本文中采用雙邊濾波器來消除空間噪聲的影響，該方法的優(yōu)點(diǎn)在于其能夠在保持圖像邊緣的前提下對(duì)圖像進(jìn)行平滑，不會(huì)破壞圖像的邊緣梯度輪廓從而不影響后續(xù)的梯度變換，雙邊濾波器［24］的表達(dá)式為：

式中，W=3σq，i和 j為圖像中像素點(diǎn)的坐標(biāo)；σq和 σr為幾何擴(kuò)散和廣度擴(kuò)散，由實(shí)驗(yàn)確定，分別用于控制空間域和值域權(quán)值的衰減［24-25］，本文中取 σq=3，σr=20。

1.2.4 梯度變換 自適應(yīng)梯度變換（adaptive gradient transform，AGT）利用雙邊高斯模型擬合圖像邊緣［20］，構(gòu)造“差函數(shù)”形式的梯度變換優(yōu)化梯度圖像，使其更加尖銳，自適應(yīng)梯度變換表達(dá)式為：

式中，σLR和σHR分別為L(zhǎng)R圖像的梯度圖像和變換后梯度圖像中梯度輪廓的銳度（擬合高斯模型的標(biāo)準(zhǔn)差），n為梯度的方向。σLR和σHR之間的關(guān)系［18］為：

式中，參量σLR利用LR圖像中每個(gè)像素點(diǎn)及其八鄰域內(nèi)的信息計(jì)算［20］。

1.2.5 圖像重建 根據(jù)MAP方法，利用優(yōu)化后的梯度圖像（▽I^HR）作為正則約束條件代入（2）式可以獲得高分辨率太赫茲圖像（IHR）［10，18］，其表達(dá)式為：

式中，▽IHR為重建后的圖像的梯度圖像。本文中，參量β取值為0.15?？紤]到太赫茲光束的強(qiáng)度服從高斯分布，因此重建時(shí)所采用的退化矩陣D為高斯矩陣。利用梯度下降算法［18-19］優(yōu)化（10）式可以獲得高分辨率太赫茲圖像。

2 實(shí)驗(yàn)結(jié)果與分析

2.1 圖像采集

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

實(shí)驗(yàn)中采用的太赫茲透射式成像系統(tǒng)（Zomega Terahertz Corporation（USA））如圖3所示。其中M1～M10為反射鏡，OAP1和OAP2為離軸拋物面鏡（off axis parabolic mirror，OAP），L1和 L2為聚乙烯透鏡。飛秒激光器產(chǎn)生的激光光束被偏振分束器（polarized beam splitter，PBS）分為抽運(yùn)光和探測(cè)光。半波板（halfwave-plate，HWP）用于改變抽運(yùn)光和探測(cè)光的比例。延線由兩個(gè)反射鏡構(gòu)成，用于調(diào)節(jié)抽運(yùn)光和探測(cè)光之間的時(shí)間延遲。抽運(yùn)光經(jīng)過衰減器，同時(shí)被反射鏡反射后照射到砷化鎵（GaAs）光電導(dǎo)天線上產(chǎn)生太赫茲波。探測(cè)光被用于探測(cè)太赫茲波，氧化銦錫（indium tin oxide，ITO）導(dǎo)電玻璃和平衡探測(cè)器（balance detector，BD）共同構(gòu)成成像系統(tǒng)的探測(cè)裝置。測(cè)量過程中，系統(tǒng)的溫度維持在21℃（294K）。成像系統(tǒng)工作頻率范圍為0.1THz～1.0THz。樣品為一個(gè)金屬制斬波器，其半徑為5cm，原始圖像的大小為55pixle×55pixle，插值后的太赫茲圖像的大小為440pixle×440pixle。樣品被放入遮蔽物中并被放置于2維電動(dòng)平移臺(tái)上，通過控制程序移動(dòng)平移臺(tái)實(shí)現(xiàn)對(duì)樣品的逐點(diǎn)掃描。掃描過程中，探測(cè)器對(duì)每一個(gè)掃描點(diǎn)測(cè)量太赫茲光譜，再將提取的太赫茲光譜按掃描次序排列成一幅太赫茲圖像。

探測(cè)器在非目標(biāo)區(qū)域采集到的時(shí)域信號(hào)和頻譜分別如圖4a、圖4b所示。在放置遮蔽物后，太赫茲信號(hào)會(huì)存在一定的時(shí)延和衰減。經(jīng)過傅里葉變換之后，對(duì)頻域信號(hào)進(jìn)行濾波即可獲得單頻信號(hào)。本文中使用0.25THz，0.50THz，0.75THz圖像來分析成像信號(hào)頻率對(duì)重建性能的影響。

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

2.2 梯度圖像提取

本文中采用離散梯度算子［-1/2，0，1/2］和［-1/2，0，1/2］T與圖像進(jìn)行卷積，分別獲得每個(gè)像素點(diǎn)水平方向和垂直方向的梯度。以0.75THz圖像為例，圖5a為太赫茲原始圖像的梯度圖像，圖5b為插值后的太赫茲圖像的梯度圖像，圖5c為對(duì)插值后的太赫茲圖像直接進(jìn)行增強(qiáng)后（未引入雙邊濾波）提取的梯度圖像，圖5d為對(duì)插值后的太赫茲圖像進(jìn)行增強(qiáng)和雙邊濾波后提取的梯度圖像，圖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可知，經(jīng)過插值之后，部分圖像細(xì)節(jié)得到恢復(fù)，邊緣變得連續(xù)；由圖5c可知，直接對(duì)經(jīng)過插值的太赫茲圖像進(jìn)行增強(qiáng)雖然可以提高其邊緣梯度輪廓的幅值，但也會(huì)放大圖像中的噪聲，從而降低梯度圖像的質(zhì)量；由圖5d可知，引入基于空間信息熵的直方圖匹配技術(shù)和雙邊濾波器對(duì)經(jīng)過插值的太赫茲圖像進(jìn)行處理后，其梯度圖像的質(zhì)量大幅提高，從而可以使得梯度變換的效果更加理想；另外，從圖5e中可以看出，梯度變換可以有效壓縮梯度輪廓的范圍（白色線條更細(xì)），提高梯度的幅值。

2.3 重建結(jié)果分析

經(jīng)過插值處理之后，太赫茲圖像的分辨率得到提高（如圖6b、圖6g和圖6l所示），部分細(xì)節(jié)信息更加清晰。相比于0.25THz的太赫茲信號(hào)，頻率為0.50THz和0.75THz的信號(hào)的光斑較小，因此0.50THz和0.75THz圖像經(jīng)過插值后邊緣更清晰。總體而言，插值處理雖然恢復(fù)了太赫茲圖像中的部分細(xì)節(jié)信息，增強(qiáng)了圖像邊緣的連續(xù)性，但插值并不能解決由于衍射極限的影響而產(chǎn)生的模糊問題。因此，在對(duì)圖像進(jìn)行插值之后，還需要對(duì)圖像進(jìn)行重建，銳化圖像邊緣。

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迭代和基于梯度變換的方法對(duì)本實(shí)驗(yàn)中所采用的成像系統(tǒng)獲得的太赫茲圖像進(jìn)行復(fù)原的結(jié)果如圖6所示。對(duì)經(jīng)過插值的太赫茲圖像采用基于梯度變換的方法重建之后，其邊緣更加尖銳，同時(shí)圖像中的空間噪聲并沒有因?yàn)檫吘変J化而被放大且沒有出現(xiàn)Lucy-Richardson迭代算法處理后出現(xiàn)的振鈴現(xiàn)象（如圖6d、圖6i、圖6n所示）。其中0.75THz圖像邊緣最尖銳，這是因?yàn)?幅圖像中獲取0.75THz圖像的太赫茲信號(hào)頻率最高，波長(zhǎng)最短，光斑最小。但由于其信號(hào)強(qiáng)度較弱，0.75THz圖像的對(duì)比度不及0.25THz和0.50THz圖像。由圖6e可見，0.25THz圖像邊緣銳化的效果最明顯，這是由于受到衍射效應(yīng)的影響，和經(jīng)過插值的0.50THz，0.75THz圖像相比，插值后的 0.25THz圖像的邊緣較為模糊，因此，該方法在對(duì)0.25THz圖像質(zhì)量的提升更加明顯。

本文中采用邊緣強(qiáng)度（ICV）和平均梯度（Idef）兩種評(píng)價(jià)標(biāo)準(zhǔn)［26］定量評(píng)價(jià)圖像重建的效果。兩種評(píng)價(jià)標(biāo)準(zhǔn)的計(jì)算式分別為：

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

經(jīng)過插值的太赫茲圖像和重建后的圖像的邊緣強(qiáng)度和平均梯度計(jì)算結(jié)果如表1所示。由表1中的數(shù)據(jù)可知，該重建方法可以有效提高太赫茲圖像的邊緣強(qiáng)度和平均梯度。其中該方法對(duì)0.25THz圖像的邊緣強(qiáng)度的提升效果最為明顯。0.50THz圖像經(jīng)過重建之后，邊緣強(qiáng)度和平均梯度最高，這是因?yàn)楹?.25THz圖像相比，其成像信號(hào)頻率較高，光斑較小，邊緣較為清晰。和0.75THz圖像相比，雖然0.50THz圖像成像信號(hào)頻率較低，但其信號(hào)強(qiáng)度更大，因此圖像的對(duì)比度更好。

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

綜上所述，利用基于梯度變換的超分辨率重建方法可以有效恢復(fù)低分辨率太赫茲圖像中的細(xì)節(jié)信息，銳化圖像邊緣，重建高分辨率太赫茲圖像。由于需要提取較為清晰的梯度圖像進(jìn)行梯度變換，該方法對(duì)于空間噪聲較小以及圖像對(duì)比度較高的太赫茲原始圖像將有更好的重建效果。

3 結(jié) 論

本文中針對(duì)實(shí)際太赫茲透射式成像系統(tǒng)采集到的原始圖像，采用有理分形插值結(jié)合基于梯度變換的超分辨率重建方法進(jìn)行了圖像重建，利用基于空間信息熵的直方圖匹配技術(shù)和雙邊濾波器優(yōu)化梯度輪廓，提高了重建性能，分析了成像信號(hào)頻率和強(qiáng)度對(duì)重建性能的影響。結(jié)果表明，該方法可以有效恢復(fù)太赫茲圖像當(dāng)中的細(xì)節(jié)信息，銳化圖像邊緣，提高圖像質(zhì)量且不會(huì)出現(xiàn)振鈴現(xiàn)象。對(duì)經(jīng)過插值的太赫茲圖像采用基于梯度變換的超分辨率重建方法處理之后，0.25THz，0.50THz和0.75THz圖像的邊緣強(qiáng)度分別提高了169%，116%和104%，平均梯度分別提高了16%，28%和24%。