基于梯度變換的太赫茲圖像超分辨率重建

郭佑東，凌福日，姚建銓，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%。