基于衍射理論的哈達瑪編碼光譜成像數據重構*

劉文龍 劉學斌 王爽 嚴強強?

1)(中國科學院西安光學精密機械研究所,中國科學院光譜成像技術重點實驗室,西安 710119)

2)(中國科學院大學,北京 100049)

哈達瑪變換光譜成像技術是一種在不犧牲分辨率的情況下,通過多通道復用的方法增加光學系統光通量,使得系統信噪比顯著高于傳統成像的新型計算光譜成像技術.本文針對系統成像過程中數字微鏡器件衍射造成的圖像降質問題,建立了基于標量衍射理論的哈達瑪編碼光譜成像圖譜數據退化模型,提出以Lucy-Richardson (L-R)算法為核心的重構數據修正算法,提高重構數據立方體的圖像質量和光譜精度.通過對成像過程的模擬和重構算法的仿真實驗,數字化生成了系統成像衍射降質的物理過程,并驗證了本文所采用修正方法的有效性.經過L-R 修正之后的復原數據立方體的光譜角距離評價結果為0.1296,圖像相似度評價因子優于0.85,相比修正前的重構數據質量有較大提升,表明算法對哈達瑪編碼光譜成像數據的重構及修正具有較好的效果.

1 引言

光譜成像作為一種有效信息獲取方式,由于其獲得的二維空間信息與一維光譜信息的數據立方體,能同時反映探測目標的幾何特征及光譜特征,更加準確地解譯探測目標的屬性,而被應用在醫學成像、環境遙感、天文觀測等眾多領域[1].孔徑編碼作為一種新型調制方式,代替了傳統成像光譜儀的狹縫[2],形成多路復用的多通道凝視光譜成像模式,使得系統能夠獲得高信噪比信號[3,4].孔徑編碼光譜成像的編碼方式可分為基于壓縮感知的稀疏編碼和基于完全采樣的完備編碼[5,6],相比于壓縮感知的稀疏編碼,完備編碼方案無需考慮信號的稀疏性,其數據重構算法簡單,信息幾乎無損失.其中,哈達瑪編碼從數學上被證明為完備編碼的最優方案之一.由于系統中調制編碼對象不同,可將哈達瑪編碼光譜成像技術(Hardmard transfer imaging spectrometer,HTIS)分為空間編碼和光譜編碼兩種技術方案,分別對應了單色散哈達瑪編碼光譜成像技術(SD-HTIS)和雙色散哈達瑪編碼光譜成像技術(DD-HTIS).

哈達瑪編碼光譜成像儀最為核心的器件是調制編碼模板,早期主要采用機械刻蝕的掩膜板.相比機械模板,DMD 具有更高的變換速度和分辨率.美國Sandia 國家實驗室第一次將數字微鏡器件(digtial micromirror devices,DMD)用于編碼成像系統中,驗證了采用DMD 器件的編碼光譜成像技術的可行性[7].美國杜克大學使用DMD 作為調制器件,研制了孔徑編碼快照成像系統,即壓縮感知成像系統,并進行了單色散和雙色散孔徑編碼光譜成像系統的數據重構[8,9].孔徑編碼光譜成像作為一種計算光譜成像技術,具有諸多傳統光譜成像技術不具備的優勢,但在數據使用過程中,極大程度上依賴于高精度的光譜和圖像重構算法.DMD 作為微納器件,對系統裝調精度有著極高的要求,另外,其自身較強的衍射效應,也給孔徑編碼成像系統帶來了能量偏移、響應不均勻、衍射降質等問題,造成重構數據立方體的圖像和光譜信號失真.因此,數據立方體的高精度重構是孔徑編碼光譜成像技術走向應用所需要解決的關鍵問題之一.Streeter 等[10]和Galvis 等[11]通過優化系統設計的方式來解決提高成像質量,提出采用照明光源的設計方案來修正基于DMD 的Hadamard變換成像光譜儀在采集過程中產生的能量偏移誤差.Stenven[12]和David[13]借助DMD 器件在實驗室搭建了一套可見光波段的孔徑編碼光譜成像系統,用于探索DMD 自身的衍射效應對系統成像的影響,并提出在光路前級添加第二塊DMD 作為衍射預補償的新思路,針對系統存在的問題,在系統設計和算法設計兩方面不斷演進解決方案.胡炳樑等[14,15]通過輻射標定的方式,結合后處理算法對系統編碼像素點錯位導致的能量偏移及系統存在的輻射響應不均勻進行了修正.Chi 等[16?18]認為光學系統的點擴散函數和光闌的非均勻照明是造成系統解碼錯誤的主要來源,通過對光學系統的編碼過程進行了標定,在此基礎上采用Gold 反卷積的方式對編碼矩陣進行了修正,消除了編碼數據相鄰空間的串擾.

光譜數據立方體的高精度重構是哈達瑪編碼光譜成像技術亟須解決的問題.在上述研究工作的基礎上,本文針對哈達瑪孔徑編碼光譜成像過程中系統衍射效應導致的點擴散彌散引起的數據降質問題[19,20],從成像機理上開展了衍射效應的正向傳輸過程的模擬.結合了光譜成像系統多通道編碼的思想和菲涅爾標量衍射效應對光學系統的降質的物理模型,開展了哈達瑪孔徑編碼光譜成像系統數據降質過程的仿真實驗.針對光學系統實際成像過程中的點擴散效應對數據重構的影響,提出了基于Lucy-Richardson (L-R)的單譜段圖像修復算法,其中,L-R 算法能夠按照泊松噪聲統計標準對給定點擴散函數的退化圖像進行反卷積迭代推演計算,充分考慮了信號的統計漲落特性,經過一定次數的迭代計算后可得到接近理想清晰圖像的最大似然估計值[21?23].通過對修正后的信號進行哈達瑪解碼運算,實現對編碼光譜數據立方體的高精度重構.

2 哈達瑪編碼光譜成像數據重構原理

2.1 成像原理和重構理論

相比傳統光譜成像技術而言,單色散哈達瑪編碼光譜成像技術在一次像面位置用多通道的編碼模板代替了狹縫,通過數字微鏡器件對空間信號進行調制,實現對信號的編碼.其光學系統的成像原理示意如圖1 所示,系統包括前置物鏡、編碼元件、準直光路和色散元件、二次成像鏡和探測器焦面.在HTIS 中,前置物鏡和一次像面相當于一臺普通的相機,獲取目標的全波段的圖像信息.一次像面位置的DMD 通過微鏡在“0”和“1”之間翻轉對目標的空間信息進行選通和關閉,實現信號的編碼.繼而,由系統中的準直鏡對編碼后的空間圖像進行準直,再通過色散光柵對各個譜段圖像進行重排,其中相鄰單譜段圖像編碼后的空間信號在光柵的作用下平移一個像素的位置.不同譜段的相鄰空間位置的信號在同一探測器像元上疊加,獲得最終的成像結果.

圖1 空間調制型哈達瑪編碼光譜成像系統光學結構示意圖Fig.1.Schematic diagram of SD-HTIS optic structure.

HTIS 對目標的圖像和光譜信息的編碼過程如圖2 所示.輸入目標的全譜段信號經過編碼模板和色散元件的調制作用,各個譜段信號在光譜維按照編碼規則進行疊加,通過對編碼模板多次變換能夠實現信號的完備編碼,使得信號疊加滿足向量方程.采用哈達瑪編碼,保證了信號矩陣的非奇異性,使多波段信號可以根據編碼模板信息和多次采樣壓縮信號在數學上能夠被有效求解.

在單色散哈達瑪編碼光譜成像技術中,每一次成像,相當于光學系統實現了一次相鄰空間位置的不同譜段信號的編碼疊加,獲得一組包含多個未知量的方程.可以表示為

其中n指探測器像元行位置,yn表示探測器像元響應,L表示編碼階數,i表示譜段數,hi表示對應i譜段信號的編碼,xn+1?i,i表示對應n+1–i像素對應空間位置處i譜段的信號.經過L次變換構成編碼向量空間記為H,編碼的光譜信號記為X,探測器接收到的信號記為Y,成像系統的數學模型可以表達為[24]

傳統哈達瑪矩陣并不適用于工程應用,一般在光學系統中采用Hadamrd 矩陣的變體Sylvester矩陣作為編碼矩陣,其特點是矩陣的不同行是由首行向量進行平移得到的[25,26].對于一個7 階哈達瑪編碼的光譜成像系統,若S循環矩陣的首行元素為[s1s2s3s4s5s6s7];xij為信號在i位置處j譜段的原始信號,其中j∈[1,2,3,···,7];表示i位置第r次成像得到的探測信號,其中r∈[1,2,3,···,7].基于該編碼矩陣的成像可以通過下式表示:

左移S矩陣具有以下特點:

1)Si=[si,si+1,···,sN,s1,···,si?1],其中i是指矩陣的行序列;

2)如果Sij=1,則=1;如果Sij=0,則=?1.且S和S′之間滿足S?1=,其中i表示S矩陣的行序號,j表示S矩陣的列序號,L表示編碼模板的階數,i,j=1,2,3,···,L.

哈達瑪編碼光譜信號的重構是對線性方程的求解.考慮到S矩陣逆變換的特點,光譜信號可通過以下方程進行重構:

2.2 基于標量衍射理論的信號修正理論

相比傳統的掃描型光譜儀,哈達瑪編碼光譜成像技術通過多通道編碼原理,將各個譜段的信號在探測器上進行疊加,使得信號具有較高的信噪比.這一理想編碼過程可以用如下公式表示:

其中,系統光譜譜段數和編碼模板階數含義一致,用L表示,Mdmd(i)為λ(i) 波段對應的DMD 模板的掩膜,Iin(λ)為輸入目標的單個譜段信號,Mdet(i)為探測器對應的掩膜,Is(λ) 表示經過DMD 調制的λ波段的編碼圖像,Id表示探測器上接收的多個譜段編碼圖像的疊加.

編碼成像的原理是通過空間光調制器件使得探測器上接收到的信號被打開或者徹底關斷,從而通過“1”和“0”實現對信號的編碼調制.但實際成像過程中,由于DMD 的微鏡片尺寸很小,在調制編碼過程中產生的衍射效應使系統點擴散函數彌散較為嚴重,使得編碼過程中探測器上接收到編碼為“0”的信息受到編碼為“1”的信息的干擾,這對于編碼圖像造成較大的影響,編碼光譜數據立方體的復原結果和原始數據比較有極大的降質,造成探測目標的信息失真.針對這一問題,基于標量衍射理論對系統光路傳輸進行推演,建立系統的成像降質過程的數理模型.

在哈達瑪編碼光譜成像系統中,(x,y)目標的λ波段經過z1距離傳輸并經過前置鏡相位調制后的衍射光場U(u,v,x,y,λ)可以表示為

其中k=2π/λ,A(λ)為λ波段對應信號強度,P1(x,y)表示前置物鏡孔徑,f1為物鏡焦距.U(u,v,x,y,λ)經過z2距離正向傳播至(s,t)平面的光場分布U2(s,t,x,y,λ)為

U2經過DMD 調制后在(w,r)平面上的光場分布表達如下:

其中T(s,t,λ) 表征DMD 的數理模型,d為單個微鏡邊長,則微鏡的數學表達如下:

哈達瑪孔徑編碼圖像是基于單波段信號的編碼疊加實現的,光柵的分光作用通過不同波段的像對應的幾何位置進行表征,即不同波段的圖像在(l,p)平面上的位置受到波長的影響,最終的(l,p)平面上的光場U4(l,p,x,y,λ)為

式中T3(w,r,λ)和P3(w,r)分別表示(w,r)平面位置處透鏡的透射調制函數和孔徑函數.定義(l,p)平面處的透鏡的調制函數和孔徑函數分別為T4(w,r,λ)和P4(w,r),最終探測器像面上的光場復振幅Ud(g,h,x,y,λ)的分布為

記哈達瑪編碼光譜成像系統單個譜段成像過程中的點擴散調制函數為Fpsf,

通過上述模型可以計算物點在探測器面對應的點擴散函數,進而通過點擴散函數和對應物點強度信息的卷積計算,獲得目標編碼調制的單譜段衍射圖像:

此時,探測器上得到降質后的編碼圖像Id可表示為

基于標量衍射模型的哈達瑪編碼光譜成像系統的構建,為系統成像的理論分析奠定了基礎.從物理光學的角度進行分析,當光學系統中的孔徑較小時,會產生較為嚴重的衍射效應,即點擴散函數受到系統中微小孔徑衍射效應的影響.在基于DMD 的孔徑編碼光學系統中,系統成像過程中將DMD 衍射帶來的降質因子通過點擴散函數和編碼調制信號的卷積向系統的終端信號進行了傳遞,從而導致了成像信息頻率混疊和圖像降質.因此,要消除衍射效應對系統復原光譜數據的影響,需要對受衍射影響的圖像數據進行頻域上的分離,在此基礎上,對單譜段圖像信號進行修正.通過對(15)式進行逆矩陣操作,即可實現單譜段信號分離,獲得單譜段圖像信息.

(16)式從數學層面上解釋了單波段圖像降質過程,本質是衍射效應下生成的點擴散函數和真實信號強度信息卷積.而通過卷積對灰度圖像進行修正的研究已有諸多成果.L-R 圖像修復算法是一種基于卷積因子盲估計的退化圖像反卷積復原優化算法.L-R 算法認為,如果模糊圖像g是由原始圖像f和模糊因子卷積得到,假定圖像各個點之間相互獨立,根據泊松分布統計模型,g和f的條件概率分布可用以下公式表示:

那圖像中g和f應滿足條件:

其中g=h ?f,g表示采集信號,f表示輸入信號,h表示信號的傳遞函數.(17)式中p(g|f) 的約束條件為(18)式,為獲得最大概率,將(17)式兩邊對數化后進行梯度運算:

由(19)式解得

最終得到的優化迭代函數為

上述優化過程能夠實現對降質圖像的修復,通過L-R 算法實現了對編碼單譜段圖像的修正,讓其更加接近理想成像結果.在此基礎上,將修復后的圖像重新進行壓縮編碼,獲得修正后的編碼圖像.

通過上述過程對完備編碼的多幀圖像進行修正,根據(4)式對光譜數據進行重構,獲得目標的數據立方體.

3 實驗結果和討論

3.1 降質分析

上述數理模型從理論上對哈達瑪編碼成像過程及數據降質問題進行了解釋說明,并從機理上提出了數據修正算法.本節在系統數理模型的基礎上設定哈達瑪編碼光譜成像系統的參數,并根據參數進行成像系統的仿真實驗.實驗參數如表1 所列.

表1 哈達瑪編碼光譜成像系統的參數表Table 1.Parameters of the Hardmard transfer imaging spectrometer system.

表1 中,f1為物鏡焦距,f2為準直鏡焦距,f3為成像鏡焦距,z1=1000 m,z2=f1,z3=f2,z4=f2+f3,z5=f3;D1為物鏡孔徑,D2為準直鏡大小,D3為成像鏡大小,D4為探測器像面大小;單個DMD 微鏡的尺寸為10.8 μm,單個探測器像元大小為30 μm,光譜儀的光譜范圍為3.7—4.8 μm,分為7 個譜段.

首先,在不考慮成像系統像差和成像衍射等因素影響的情況下,開展哈達瑪編碼成像仿真工作,以AVIRIS 數據為仿真輸入數據,仿真7 階哈達瑪編碼光譜成像過程,并針對10 dB 噪聲的信號進行重構,結果如圖3 中的偽彩色圖所示,其中圖3(a)為原始輸入數據,圖3(b)為編碼得到的多幀圖像,圖3(c)為重構的數據立方體.

圖3 AVIRIS 數據仿真生成及復原光譜數據 (a)原始輸入數據;(b)編碼生成的多幀圖像數據;(c)復原后的光譜數據立方體Fig.3.Simulation of AVIRIS data and the restoration spectral data:(a) Original input data;(b) multi-frame image data generated for encoding;(c) restored spectrum data cube.

通過圖像相似度因子對重構圖像進行評價,將重構數據立方體的單譜段圖像和原始輸入圖像進行比較,兩者的相似度因子的評價結果如表2 所列.

取原始數據和復原數據的物點光譜曲線進行比較,原始數據和復原光譜數據同一目標的光譜曲線進行評價,光譜角距離評價因子的均值為0.006,反演數據和原始數據幾乎一致.上述仿真結果證明,在不存在衍射和系統像差的情況下,僅由系統噪聲對數據重構產生的影響較小.

實際成像過程中,哈達瑪編碼光譜成像過程中受到系統中微納光學器件的影響,造成圖像結果降質,根據預設的光學系統參數,在衍射成像模型的基礎上開展仿真實驗,通過(13)式計算獲得目標點在探測器表面的點擴散函數,圖4(a)表示得到的像點擴散函數形狀;圖4(b)表示沿橫向過像點中心點位置的歸一化的強度分布曲線,其中橫向表示像素點序號,縱向表示各個像素點上的信號強度和輸入信號強度的比值.

歸一化的系統點擴散分布如圖4 所示,可以看出,受衍射效應的影響,系統的能量彌散到多個探測器像元之上,使得生成的壓縮編碼圖像質量下降.根據仿真實驗結果直接重構得到的數據立方體的光譜角距離達到了0.502,相比系統噪聲帶來的影響,數字微鏡器件衍射對數據重構的影響更為嚴重.因此,減小和消除系統的衍射降質影響是哈達瑪編碼光譜成像數據重構急需解決的問題.

圖4 歸一化的系統點擴散函數分布圖 (a)像點的點擴散函數圖;(b)歸一化的像點強度分布曲線Fig.4.Normalized system point spread distribution diagram:(a) Point spread function diagram of image spot;(b) distribution curve of normalized image point intensity.

3.2 衍射圖像修復仿真實驗與數據分析

根據理論降質模型,首先對編碼圖像進行分解,得到圖像的多個譜段的編碼信息,通過L-R 算法對分解得到的單譜段圖像進行修復.由于光學系統的點擴散函數更接近于高斯分布,因此修復算法中用于逆卷積的核函數采用高斯核,本實驗中采用的高斯濾波器維度為11×11,經過優化迭代計算,最終采用的濾波器標準值為0.730.經過上述單譜段圖像的修復之后,重新對各個譜段信息進行壓縮計算,獲得修正后的編碼圖像如圖5 所示.和修正前的圖像相比,修正后圖像的對比度和圖像的細節信息得到較大提升,圖5 中紅色框圖標記位置兩者對比十分明顯.

圖5 Hadamard 編碼單譜段圖像 (a)衍射降質圖像;(b)復原修正圖像Fig.5.Encoded single-spectrum image:(a) Image with the influence of diffraction;(b) image optimized by the L-R repair algorithm.

以理想編碼圖像為基準,對受衍射影響的編碼圖像和經過修正的編碼圖像進行比較,相似度因子評價的結果如表3 所列.相比受衍射影響的圖像,修正后的數據更加接近理想的編碼圖像.這一結果證明了L-R 算法能夠有效實現對降質數據的修正,提升數據重構的質量.

表3 受衍射影響的編碼圖像和經過修正的編碼圖像的相似度因子評價Table 3. Similarity factor evaluation of uncorrected coded image and corrected encoded image.

在此基礎上,通過(4)式實現對數據立方體的重構,圖6 分別給出了原始輸入數據立方體、未經過衍射修正的數據重構的立方體和修正后的數據重構的數據立方體.

圖6(b)和圖6(c)直觀地反映了修正前后重構的數據立方體的質量.對上述復原光譜數據立方體的準確性進行評價,將數據立方體的各個譜段和輸入光譜圖像進行比較,相似度評價結果如表4 所列.

表4 未修正SSIM 與修正后SSIM 的相似度評價Table 4. Similarity evaluation of uncorrected SSIM and corrected SSIM.

隨機選取圖6 數據立方體中3 組光譜曲線進行比較,結果如圖7 所示,其中橫坐標位光譜譜段數,縱坐標表示數字信號的灰度值.可以看出,修正后的光譜數據和原始輸入數據更加接近,而未經過修正的重構光譜數據和原始輸入偏差較大.以輸入數據立方體為對照數據,分別對直接重構的數據立方體和經過修正后復原的數據立方體包含的全部光譜信息的光譜角距離進行評價,其中直接對衍射編碼復原譜角距離均值為0.508,而經過L-R 修正之后的復原得到的數據立方體的光譜角距離均值為0.130.

圖6 不同譜段光譜圖像 (a)原始輸入和理想復原圖像;(b)經過衍射效應影響復原的多譜段圖像;(c)經過L-R 修正后復原的多譜段圖像Fig.6.Spectral images of different bands:(a) Original input and ideal restoration image;(b) multi-spectral image restored by diffraction effect;(c) multi-spectral image restored after L-R correction.

上述實驗說明,空間調制型哈達瑪編碼光譜成像系統中,衍射是造成重構數據失真的主要因素.在哈達瑪編碼光譜成像系統衍射模型的基礎上,通過L-R 算法能夠有效地實現信號的修正,提高重構信號的質量.

4 結論

本文針對多通道孔徑編碼光譜成像技術中的數據衍射降質問題,結合哈達瑪編碼光譜成像的機理和標量衍射理論,建立了哈達瑪編碼光譜成像的數理模型,通過對模型的分析,發現哈達瑪編碼光譜成像系統中的衍射效應是導致系統重構數據退化嚴重的主要因素,并構建了基于點擴散函數的哈達瑪編碼光譜數據的退化模型.在此基礎上,設計了一種基于L-R 算法的數據逆卷積重構算法,該算法結合了光譜成像系統的多通道編碼的思想和光學系統實際成像過程中的點擴散響應.算法首先編碼壓縮的圖像進行逆向分解,獲得多個譜段的編碼圖像,針對不同譜段的圖像使用L-R 算法對降質圖像進行迭代優化,降低了衍射降質因子對圖譜數據的影響,從而提高了重構光譜數據立方體的精度.該方法從數理層面上開展了對系統重構數據降質因素的分析,并通過降質模型,實現了對復原圖譜數據的修正.這項工作提高了采用DMD 作為編碼器件的哈達瑪光譜成像系統的數據立方體重構精度,為孔徑編碼光譜成像在多學科的應用研究上發揮重要作用.