吳佳琛，曹良才

（清華大學精密儀器系精密測試技術及儀器國家重點實驗室，北京 100084）

0 引言

基于透鏡的成像系統(tǒng)通過設計光學鏡組，將物空間發(fā)出的光線會聚到像空間，由像感器直接記錄下物體的圖像。為了提高成像質量，光學鏡組通常采用多片透鏡來校正像差，因此光學鏡組占據(jù)成像系統(tǒng)大部分體積和重量，限制了成像系統(tǒng)的小型化和輕量化。基于編碼掩膜的無透鏡成像將具有特定圖案的掩膜代替透鏡置入成像系統(tǒng)內，入射光受到掩膜圖案的調制，在像感器上形成編碼圖像，最后通過算法從編碼圖像中恢復出原始圖像，相較于透鏡成像系統(tǒng)具有結構簡單、易于構建的特點。編碼掩膜無透鏡成像打破了物點到像點一一對應的成像方式，將成像的重心由硬件轉移到算法上，并有可能進一步提取深度、波長等多維度的物體信息［1，2］。如何建立原始圖像與編碼圖像之間的聯(lián)系，并通過求解逆問題重建圖像則是編碼掩膜無透鏡成像中的關鍵問題。

早期編碼掩膜成像技術主要應用于高能射線成像中，也稱為編碼孔徑成像。這是由于在高能射線波段中，各種材料的折射率都近似等于1，并有很高的吸收率，無法采用常規(guī)光學元件進行折射成像。而無需光學元件的小孔成像技術具有極低的通光量，需要長時間曝光才能獲取清晰圖像。1961 年，由MERTZ L 和YOUNG N 提出了波帶片編碼成像技術（Zone Plate Coded Imaging，ZPCI）［3］，用波帶片代替透鏡實現(xiàn)了編碼成像。1968 年，DICKE R［4］和ABLES J［5］對小孔成像原理進行拓展，分別獨立提出了用于X 射線和伽馬射線成像的隨機孔徑陣列。之后FENIMORE E E 等提出的均勻冗余陣列（Uniformly Redundant Array，URA）［6］，以及GOTTESMAN S R 等提出的改進的均勻冗余陣列（Modified Uniformly Redundant Array，MURA）［7］已成功應用于工業(yè)伽馬射線相機和高能天文望遠鏡系統(tǒng)［8］。這些方法主要基于幾何光學模型，未充分考慮入射光透過掩膜版所產(chǎn)生的衍射效應。在可見光波段，由于波長和掩膜版圖案特征尺寸相當，衍射效應不能被忽視，因而限制了其應用推廣。為此，研究人員設計了各種類型的編碼掩膜，如可分離掩膜［9-11］、散射屏［12，13］、菲涅爾波帶片［14-16］等，并提出了相應的重建算法以提高成像系統(tǒng)的魯棒性。

菲涅爾孔徑編碼成像技術將波帶片編碼成像技術拓展至可見光領域，其原理是利用菲涅爾波帶片將目標物體編碼成與同軸全息圖具有相同形式的編碼圖像，因此可借鑒數(shù)字全息成像算法來重建圖像。SHIMANO T 等提出采集至少四幅不同相位的波帶片編碼圖像用來實現(xiàn)無孿生像的重建［14］。本文作者提出全變差正則化方法［15］和深度學習方法［16］實現(xiàn)了菲涅爾孔徑編碼成像的單次采集成像。菲涅爾孔徑編碼成像分辨率與波帶片最外環(huán)寬度相當，而波帶片最外環(huán)寬度隨著菲涅爾波帶片的半徑增加而縮小，因此大幅面的像感器可以獲得高分辨率的圖像，但是對像感器成本提出了更高的要求。

為了在不損失圖像分辨率的情況下降低硬件成本，可以使用多個小尺寸像感器來代替大尺寸像感器進行圖像采集。由于多個像感器只能獲得部分測量值，因此需要通過壓縮感知技術來實現(xiàn)菲涅爾孔徑（Fresnel Zone Aperture，F(xiàn)ZA）編碼成像。壓縮感知是DONOHO D 等數(shù)學家于2006 年左右提出的一種新的信息獲取指導理論［17-19］，能夠通過稀疏采樣，在遠低于奈奎斯特采樣頻率條件下，通過尋找欠定線性系統(tǒng)稀疏解的方式，獲得更好的圖像重構效果。壓縮采樣技術已成功應用于磁共振成像［20］、相位恢復［21，22］、熒光顯微［23］和全息成像［24，25］中。

本文基于全變差正則化重建算法，驗證了FZA 編碼成像系統(tǒng)在波帶片常數(shù)小于某一閾值時具有壓縮重建能力，提出一種基于部分采樣的FZA 編碼圖像的壓縮成像模型，能夠實現(xiàn)在部分測量數(shù)據(jù)缺失的情況下重建高質量圖像，并對不同采樣數(shù)據(jù)量下重建圖像質量進行了討論，分別針對矩形采樣和輻射線采樣的兩種數(shù)據(jù)采集方式進行了測試分析。

1 菲涅爾孔徑編碼成像模型

1.1 編碼圖像的記錄

透光率連續(xù)變化的波帶片被稱為伽柏波帶片（Gabor Zone Plate，GZP），其振幅傳遞函數(shù)可表示為

式中，(x，y)為空間坐標，r1為波帶片常數(shù)，表示波帶片最內圈環(huán)帶的半徑。但現(xiàn)有制造工藝難以加工正弦透過率型波帶片。通常采用具有二值透過率的菲涅爾波帶片（Fresnel Zone Plate，F(xiàn)ZP）代替GZP 實現(xiàn)相應的功能。因為FZP 在編碼掩膜成像中充當編碼孔徑的功能，也稱之為菲涅爾孔徑（FZA）。FZA 編碼成像原理如圖1 所示，物體到FZA 的距離為z1，F(xiàn)ZA 到像感器距離為z2。物體受非相干光照明，其透射或反射的光線經(jīng)FZA 調制后，由像感器記錄下編碼圖像，最后利用優(yōu)化算法計算重建出原始圖像。

圖1 FZA 編碼成像示意圖Fig.1 Principle diagram of FZA coded imaging

對于編碼掩膜成像而言，像感器采集的圖像可以表示為物體的像和掩膜版投影的卷積，即

式中，符號*表示卷積算符；O(x，y)表示待復原的物體圖像；T(x，y)表示掩膜版投影強度分布，當z1?z2時，T(x，y)與掩膜版透過率函數(shù)等價；e(x，y)表示成像系統(tǒng)中各種因素引入的噪聲，包括光電探測器噪聲、環(huán)境光噪聲，以及衍射效應引起的誤差等。將T(x，y)中的余弦項用復指數(shù)的形式表達，即

1.2 反向傳播重建

若要重建原始物體圖像，可以數(shù)值模擬相干光對編碼圖像的衍射過程。采用角譜法計算衍射光場，可得重建光場的表達式為

2 菲涅爾孔徑編碼圖像的壓縮重建方法

2.1 壓縮重建的RIP 條件

CANDèS E 和TAO T 提出的約束等距性（Restricted Isometry Property，RIP）性質是判斷傳感矩陣能否實現(xiàn)信號壓縮采樣的一個重要標準［26］。RIP 性質的表達式為

式中，θ為任意稀疏度為K的信號，δ為常數(shù)且δ∈(0，1)，A為傳感矩陣。若傳感矩陣A滿足式（6），則可以稱傳感矩陣A滿足RIP 性質，此時優(yōu)化算法能夠以高概率精確恢復原始信號。CANDèS E 和TAO T 還證明了獨立同分布的高斯隨機測量矩陣可以成為普適的壓縮感知傳感矩陣［26］。但在實際中判斷矩陣是否滿足RIP 性質被證明是十分困難的［27］，而且RIP 只是一個矩陣是否可以作為觀測矩陣的充分不必要條件，在實際應用中，人們更關心所設計的測量矩陣在能夠準確恢復信號的條件下，所需的測量數(shù)量以及可以恢復的信號的稀疏度。因此，可將高斯隨機測量矩陣對信號的恢復能力作為參照進行類比。本文采用正弦型波帶片的一維徑向函數(shù)對應的卷積矩陣作測量矩陣，對長度N= 2 048，稀疏度K= 200 的信號進行重構，用來驗證菲涅爾孔徑編碼成像系統(tǒng)是否滿足RIP 條件。模擬實驗中分別測試了波帶片常數(shù)為0.8 mm、0.7 mm、0.6 mm、0.5 mm 時的重構情況，信號和波帶片的采樣間隔為10 μm。其函數(shù)圖像如圖2（a）所示，信號重構結果如圖2（b）所示，同時給出高斯隨機傳感矩陣作為參考?？梢钥吹絩1越小，重構誤差越小，當r1= 0.5 mm時，重構性能與高斯隨機傳感矩陣幾乎一致。即波帶片常數(shù)越小，在壓縮重建中對信號的重構誤差也就越小。

圖2 FZA 傳感矩陣與高斯隨機傳感矩陣對信號的恢復能力比較Fig.2 Comparison of signal recovery ability between FZA sensing matrix and Gaussian random sensing matrix

2.2 全變差正則化

通過構造優(yōu)化目標函數(shù)實現(xiàn)編碼圖像的壓縮重建。若以矩陣向量相乘的形式抽象出菲涅爾孔徑編碼成像的數(shù)學模型，則有

假設待復原圖像的像素數(shù)為Nx×Ny=Nxy。其中u∈RNxy表示目標圖像，f∈RNxy表示編碼圖像，c和e分別表示常數(shù)和噪聲，F(xiàn)∈CNxy×Nxy表示二維傅里葉變換矩陣，F(xiàn)-1∈CNxy×Nxy則是對應的逆傅里葉變換矩陣；H∈CNxy×Nxy是對角矩陣，其對角線上的元素是菲涅爾衍射傳遞函數(shù)的離散采樣值；Re{·}表示取實部操作；K算子則是結合了F-1HF以及取實部操作，代表了整個成像系統(tǒng)的前向模型。在菲涅爾孔徑編碼成像應用中，物函數(shù)u始終為實函數(shù)，且菲涅爾衍射傳遞函數(shù)H為徑向對稱函數(shù)，則Re{(F-1HF)u}等價于(F-1Re{H}F)u，即K算子為線性算子，因此可以用線性回歸模型相關的優(yōu)化算法對原始圖像進行求解。

由于孿生像和原始像以及二者之和都能滿足式（7），因此圖像重建具有多個解，屬于不適定性問題，需要施加正則化項使解保持唯一且穩(wěn)定。孿生像本質上是原始圖像的離焦圖像，聚焦圖像的邊緣清晰銳利，邊緣以外的區(qū)域則過渡平滑，而離焦圖像在圖像邊緣會產(chǎn)生衍射條紋，并且衍射條紋隨著傳播距離增大擴散到整個圖像。反映在圖像的梯度上，聚焦圖像的梯度大部分趨于零值，呈現(xiàn)稀疏性，而離焦圖像則是非稀疏的。所以聚焦圖像的全變差要比離焦圖像的全變差小得多。圖3 給出了離焦圖像和聚焦圖像梯度域的圖像及相應的直方圖分布，其中聚焦圖像在梯度域具有明顯的稀疏特性。因此，可以采用全變差正則化實現(xiàn)抑制孿生像的效果，全變差正則化項可寫為

圖3 聚焦圖像和離焦圖像的梯度域稀疏性比較Fig.3 Comparison of sparsity in gradient domain between focused and defocused images

式中，D表示差分算子，m、n表示二維離散圖像的像素索引。

2.3 優(yōu)化目標函數(shù)

根據(jù)上述分析，設計如下優(yōu)化目標函數(shù)

式中，Ω表示采樣區(qū)域的指標集，‖ ·‖Ω表示在僅在指標集Ω上計算l2范數(shù)，τ＞0 為正則化參數(shù)。值得注意的是，在誤差項Ku-f前添加一個差分算子D，能夠有效消除編碼圖像中常數(shù)項的干擾，提高成像質量。

式（9）可以通過交替方向乘子法（Alternating Direction Method of Multipliers，ADMM）進行求解。ADMM 方法可將復雜的問題分解成若干個易于求解的子問題，縮小了問題的規(guī)模，降低了求解難度。首先將式（9）改寫為等價的約束優(yōu)化問題，即

式中，Dh和Dv分別為水平方向和垂直方向的差分算子，bh=Dhf，bv=Dvf。該問題的增廣拉格朗日函數(shù)定義為

式中，y1、y2、y3、y4為尺度對偶變量，μ＞0、η＞0 為懲罰參數(shù)。通過給定中間變量wh、wv、zh、zv和尺度對偶變量y1、y2、y3、y4的初值，ADMM 算法在第k次迭代過程中依次求解如下子優(yōu)化問題

并更新對偶變量

式中，β∈(0，( 5 +1) 2 )用來控制更新的步長。通過交替迭代優(yōu)化u、wh、wv、zh、zv、y1、y2、y3、y4最終恢復原始圖像u。

3 重建結果與分析

3.1 數(shù)值重建結果

在數(shù)值模擬實驗中，原始圖像尺寸為256×256 像素，像素間隔為10 μm；掩膜版圖案為菲涅爾波帶片，波帶片常數(shù)為0.25 mm，掩膜版尺寸為512×512 像素；完整的編碼圖像尺寸為768×768 像素。取完整編碼圖像中心的256×256 像素部分作為像感器所采集的圖像，現(xiàn)根據(jù)該采集圖像，通過減少該采集圖像的數(shù)據(jù)量，測試壓縮重建算法對原始圖像恢復的性能。首先采用矩形采樣區(qū)域對編碼圖像進行采樣，采樣比例以256×256 像素的采集圖像作為100%的測量數(shù)據(jù)進行計算。由于重建圖像的像素值分布范圍與原圖不盡相同，因此采用相關系數(shù)（Correlation Coefficient，CC）對重建圖像質量進行評價。對于兩幅圖像A和B，相關系數(shù)定義為

為了進一步分析采樣數(shù)據(jù)量對重建圖像質量之間的相關性，對分類主觀圖像質量（Categorical Subjective Image Quality，CSIQ）數(shù)據(jù)庫［28］中30 幅圖像進行測試，每幅圖像均采用ADMM 算法迭代200 次，最終重建圖像的相關系數(shù)分布如圖4（c）所示。圖中箱線圖框內的線條表示樣本中位數(shù)，框的下邊緣和上邊緣分別表示下四分位數(shù)和上四分位數(shù)，即數(shù)據(jù)從小到大排序后處于25%和75%位置上的值，二者之差為四分位距（Inter-Quartile Range，IQR）。圓圈表示離群值，即超出上下四分位數(shù)1.5 倍IQR 的值，豎線的上下兩端分別表示除去離群值之后的最大值和最小值。最后用折線圖連接了每種采樣數(shù)據(jù)量下相關系數(shù)的均值，當采樣數(shù)據(jù)量小于50%，圖像重建質量隨著采樣數(shù)據(jù)量的減少而迅速下降，在采樣數(shù)據(jù)量僅為6.3%的情況下，重建圖像相關系數(shù)均值降至0.30。

圖4 基于矩形采樣區(qū)域的壓縮重建結果Fig.4 Compressive reconstruction results based on rectangular sampling area

實際上，矩形采樣并沒有充分考慮編碼圖像在頻譜域中的能量分布情況。由于像感器對斜入射光線不敏感，實際視場角被限制在小范圍內，也就是波帶片投影的位移遠小于波帶片的直徑，此時，像感器的各個像素接收到光強僅來自于波帶片投影對應的局部小區(qū)域的疊加，因此，編碼圖像整體上呈現(xiàn)和波帶片類似的頻率分布，其特點為圖像中心低頻成分居多，圖像邊緣高頻成分居多。矩形采樣只對圖像中心部分進行采樣，會導致重建圖像的高頻信息缺失而變得模糊不清。隨著衍射距離的增大，菲涅爾衍射圖像逐漸向夫瑯禾費衍射圖像轉變，而夫瑯禾費衍射圖像和原始圖像的頻譜分布具有相同的形式［29］，因此，菲涅爾衍射圖像處于空域圖像向頻譜圖像轉化的中間狀態(tài)，與原始圖像的頻譜能量分布具有一定相似性。由于大多數(shù)自然圖像的頻譜能量都集中在低頻，所以應該在靠近中心的地方進行更密集的采樣，以匹配能量分布。采用如圖5（a）所示輻射線采樣圖案，輻射線由圖像中心向外發(fā)散，并均勻遍布整幅圖像，輻射線數(shù)量從左到右分別為64、48、32，相應的采樣數(shù)據(jù)量分別為完整編碼圖像的13.9%、10.7%、7.3%。這種采樣方式可以保證中心和邊緣的圖像信息都能被采集到，并且在圖像中心處得到更多的采樣數(shù)據(jù)，可通過線陣相機來實現(xiàn)。類似的采樣方案已應用于相位恢復［21］、計算機斷層掃描［30］、核磁共振成像［31］等領域。相對于矩形采樣而言，輻射線采樣能夠以更少的采樣數(shù)據(jù)量重建出高質量的圖像。如圖5（b）所示，輻射線采樣僅用10.7%的采樣數(shù)據(jù)重建的圖像與矩形采樣使用56.3%的采樣數(shù)據(jù)重建的圖像有相同水平的重建質量，重建圖像相關系數(shù)均為0.89。從圖5（c）中也可看出，基于輻射線采樣的圖像重建質量隨著采樣數(shù)據(jù)量的減少下降緩慢，即使在采樣數(shù)據(jù)量僅為5.7%的極端情況下，重建圖像相關系數(shù)均值仍有0.77。

圖5 基于輻射線采樣區(qū)域的壓縮重建結果Fig.5 Compressive reconstruction results based on radial line sampling area

3.2 實驗重建結果

實驗中采用液晶顯示屏（Liquid Crystal Display，LCD）顯示待成像的圖像，LCD 被放置在距離無透鏡相機約30 cm 的地方。目標圖像為一含有“HOLOLAB”字樣的矩形標志，其顯示在屏幕中的尺寸為130 mm×140 mm。無透鏡相機所采用像感器為QHY163M，像素間隔為3.75 μm。菲涅爾波帶片放置在距離像感器3 mm 處，取其中心2 048×2 048 像素作為采集到的完整的編碼圖像。分別使用完整的編碼圖像和不同采樣模式進行圖像重建。使用圖6（a）中的采樣模式對測量數(shù)據(jù)進行重建，對應的重建結果如圖6（b）所示，所展示的圖像為重建圖像中心500×500 像素部分。結果表明，在采集數(shù)據(jù)量僅為7.3%的情況下，該方法仍能有效識別出字母。通過圖6（c）字母“O”的橫截面強度可以看出，在三種不同采樣模式下，圖像對比度并未見明顯下降，驗證了該方法具有良好的圖像邊緣保持特性。矩形采樣模式可由多塊小尺寸面陣像感器拼接實現(xiàn)，輻射線采樣模式可由單個線陣像感器多次旋轉采集，或者采用多個線陣像感器拼接實現(xiàn)?？紤]到實際裝配的可行性，也可采用中心放置面陣像感器，周圍放置線陣像感器的混合排布方式，并對裝配完成后像感器陣列進行一次標定，減小裝配誤差對圖像重建的影響。

圖6 實驗測量圖像的壓縮重建結果Fig.6 Compressive reconstruction results of experimental measurement images

4 結論

本文提出了一種用于菲涅爾孔徑編碼成像的壓縮采樣方法，構建了從不完全測量中恢復原始圖像的壓縮重建模型。分析了菲涅爾孔徑編碼壓縮采樣模型的不相關特性，比較了菲涅爾孔徑編碼的測量矩陣和隨機高斯測量矩陣對信號的恢復能力，驗證了一維情況下波帶片常數(shù)小于0.5 mm 時，重構性能與高斯測量矩陣幾乎一致。提出并推導了基于ADMM 的壓縮重建方法，對不同采樣數(shù)據(jù)量下的矩形采樣模式和輻射線采樣模式進行了定量分析，驗證了輻射線采樣模式相比矩形采樣模式具有更高的圖像采樣效率。實驗表明，該方法僅通過7.3%的測量數(shù)據(jù)就可以獲得質量良好的圖像，可以用于菲涅爾孔徑編碼像感器陣列成像。所提出的方法能夠在保持重建圖像質量的同時，大幅降低測量數(shù)據(jù)，驗證了基于菲涅爾孔徑編碼成像構建多像感器架構的可行性。