李能菲,孫宇松,黃見

(1安徽職業(yè)技術(shù)學院機電工程學院,安徽 合肥 230011;2中國科學院合肥物質(zhì)科學研究院安徽光學精密機械研究所,中國科學院大氣光學重點實驗室,安徽 合肥 230031;3中國科學技術(shù)大學研究生院科學島分院,安徽 合肥 230026;4先進激光技術(shù)安徽省實驗室,安徽 合肥 230031)

0 引言

關(guān)聯(lián)成像是一種新型的計算成像技術(shù),其與傳統(tǒng)面陣成像的區(qū)別在于關(guān)聯(lián)成像的成像器件為無空間分辨率能力的點探測器,由于點探測器相對于面陣CCD/CMOS探測器具有光譜選擇范圍大、量子效率高等優(yōu)勢,使得關(guān)聯(lián)成像在現(xiàn)有面陣成像無法工作的非可見光波段具有潛在的成像優(yōu)勢,近年來成為光學、電子學和計算機科學等眾多交叉學科的研究熱點。關(guān)聯(lián)成像起源于糾纏光子光源[1],因此部分研究人員也將關(guān)聯(lián)成像稱為量子成像。后續(xù)的研究表明經(jīng)典光源[2?4]也可以實現(xiàn)關(guān)聯(lián)成像。在傳統(tǒng)的關(guān)聯(lián)成像系統(tǒng)中,光源發(fā)出的光被分成兩路,一路稱為參考光路,另一路稱為探測光路。參考光路的光不與成像物體接觸,直接被面陣探測器接收;探測光路的光與成像物體相互作用后被一點探測器收集,聯(lián)合參考光路記錄的光強分布與探測光路記錄的光強經(jīng)關(guān)聯(lián)運算可重建出物體圖像。研究發(fā)現(xiàn),使用空間光調(diào)制器產(chǎn)生光強空間分布預(yù)知的調(diào)制光源可以省略掉參考光路[5],大大簡化了關(guān)聯(lián)成像系統(tǒng),推動了關(guān)聯(lián)成像進一步向?qū)嵱没较虬l(fā)展。目前,關(guān)聯(lián)成像在多光譜成像[6?8]、紅外成像[9]、太赫茲成像[10,11]、氣體檢測[12]、偏振成像[13]、三維成像[14?16]以及目標跟蹤[17?19]等領(lǐng)域都展現(xiàn)了其廣闊的應(yīng)用前景。

不同于傳統(tǒng)面陣成像,關(guān)聯(lián)成像通過一系列被掩膜圖案調(diào)制的光源照射場景,無空間分辨能力的點探測器記錄相應(yīng)的光強,通過聯(lián)合光強與掩膜圖案做關(guān)聯(lián)運算來重建場景圖像,關(guān)聯(lián)成像的這種成像機制決定了其是一種以犧牲時間分辨率換取空間分辨率的成像技術(shù),成像空間分辨率越高,所需的調(diào)制散斑越多,使得采樣時間長、成像效率低。因此,如何在保持較高成像質(zhì)量的基礎(chǔ)上、在低采樣數(shù)下獲取更高的空間分辨率與幀頻成為關(guān)聯(lián)成像走向?qū)嶋H應(yīng)用必須突破的關(guān)鍵科學問題。文獻[7,8]和文獻[20]采用相互正交的隨機編碼復(fù)用分別實現(xiàn)了多光譜關(guān)聯(lián)成像、多物體成像與加密,由于編碼信息空間分布是隨機的,使得在圖像重構(gòu)時需利用壓縮感知優(yōu)化算法來對隨機編碼對應(yīng)的欠定方程組進行求解,導(dǎo)致圖像復(fù)原時間消耗大大增加。

本文提出了一種余弦編碼復(fù)用高空間分辨率關(guān)聯(lián)成像技術(shù),通過構(gòu)造多個低空間分辨率的余弦編碼散斑復(fù)用為高空間分辨率調(diào)制散斑,對成像物體進行調(diào)制照明,單像素探測器接收被調(diào)制物體信號,由迭代算法復(fù)原出成像目標低空間分辨率的混疊圖像;鑒于編碼信息所特有的確定性頻譜結(jié)構(gòu),利用數(shù)字圖像處理解碼重構(gòu)出多個低空間分辨率物體圖像,進而拼接為高空間分辨率目標圖像。理論介紹了余弦編碼復(fù)用高空間分辨率關(guān)聯(lián)成像技術(shù)的實現(xiàn)方法,數(shù)值仿真驗證了此方法的有效性。

1 傳統(tǒng)關(guān)聯(lián)成像方法

圖1為關(guān)聯(lián)成像系統(tǒng)示意圖,投影器件產(chǎn)生變化的調(diào)制光源Si(x,y),對成像物體O(x,y)進行調(diào)制照明,調(diào)制照明光源與成像物體相互作用產(chǎn)生的后向散射信號被光學透鏡匯聚后,由無空間分辨能力的單像素探測器收集,從而實現(xiàn)光電轉(zhuǎn)換,數(shù)據(jù)采集單元同步采集相應(yīng)的電信號,并存儲在計算機中供后續(xù)圖像復(fù)原使用。設(shè)與調(diào)制信號Si(x,y)對應(yīng)的探測光強為Ii,則Ii可表示為

圖1 關(guān)聯(lián)成像系統(tǒng)示意圖Fig.1 Schematic diagram of the ghost imaging system

式中:L為調(diào)制次數(shù)。對于空間分辨率為N×M的物體成像,一般情況下全采樣時L=NM。對(1)式求解常用的方法是迭代算法,即

傳統(tǒng)關(guān)聯(lián)成像對空間分辨率為N×M的物體關(guān)聯(lián)成像時,全采樣需要N×M個調(diào)制散斑,散斑數(shù)據(jù)量隨著空間分辨率提高呈平方關(guān)系增加。如,對128×128分辨率物體成像需要調(diào)制散斑數(shù)為16384,假設(shè)調(diào)制頻率為20 kHz,則全采樣一幅圖像的時間至少為0.8 s;當對256×256分辨率物體成像,需要調(diào)制散斑數(shù)為65536,在同樣的調(diào)制頻率下,則全采樣一幅圖像的時間至少為3.2 s,成像幀頻大幅降低。可以看出,成像物體的空間分辨率越高,需要的調(diào)制散斑越多,相對應(yīng)采樣次數(shù)就越多,進而不可避免地增加了成像時間。

2 余弦編碼復(fù)用高空間分辨率關(guān)聯(lián)成像方法

本研究提出采用多個余弦編碼的低空間分辨率復(fù)用調(diào)制散斑來實現(xiàn)N×M高空間分辨率的關(guān)聯(lián)成像。為便于表述,這里假設(shè)將N×M空間分辨率的物體等分為2部分,即P1:(1:N,1:N)和P2:(1:N,N+1:M),其中M=2N。通過獲取2個(1:N,1:N)低空間分辨率的圖像來實現(xiàn)對N×M高空間分辨率物體的成像。構(gòu)造兩個N×N二維余弦結(jié)構(gòu)編碼矩陣,分別對應(yīng)被成像物體的兩個部分,記為EP1和EP2,即

式中:a為編碼矩陣的對比度,b為編碼矩陣的偏置,fx1和fy1為編碼矩陣EP1的頻率,fx2和fy2為編碼矩陣EP2的頻率,φ0為編碼矩陣的初始相位。將余弦編碼矩陣和基于Hadamard基生成的圖案Si(x,y)相融合,可生成兩個空間分辨率大小為N×N的調(diào)制散斑EP1·Si(x,y)和EP2·Si(x,y),符號·表示點乘運算。將這兩個調(diào)制散斑橫向排列,可構(gòu)成空間分辨率為N×M的調(diào)制散斑,該高空間分辨率的調(diào)制散斑同時對成像物體的P1和P2部分調(diào)制照明,單像素探測器收集空間分辨率為N×M物體的反射或者透射信號Ui,可表示為

式中:TA和TB分別對應(yīng)成像物體空間分辨率大小為N×N的P1部分和P2部分;T表示成像物體的混疊圖像,辨率為N×N,數(shù)值上T=EP1·TA+EP2·TB。鑒于Si(x,y)為基于Hadamard基的正交散斑,可通過線性迭代運算對(4)式中T進行求解,即

在獲取T后,如何從分辨率為N×N的混疊圖像T中解調(diào)重構(gòu)出TA和TB是本技術(shù)的核心。由于余弦編碼矩陣EP1和EP2具有特定的頻譜特性,可以利用傅里葉變換理論來處理。對混疊圖像T做傅里葉變換,將其從空間域轉(zhuǎn)換到傅里葉頻域,即

式中符號F表示對圖像做二維傅里葉變換。將歐拉公式cosx=(eix+e?ix)/2代入(6)式,由傅里葉移位定理可推導(dǎo)出

式中:FTA和FTB分別為低空間分辨率圖像TA和TB對應(yīng)的頻譜,fx和fy表示頻譜坐標。由(7)式可以看出,在頻域中由于TA和TB的編碼頻率組合(fxi,fyi)不同,TA和TB的頻譜信息被頻移到混疊圖像頻譜的不同高頻區(qū)域。因此,在混疊圖像頻譜圖中會出現(xiàn)兩對頻點(對稱性),對應(yīng)圖像TA和TB頻譜信息中的最大值。首先將混疊圖像頻譜中一對頻點平移到傅里葉域的中心,同時保持另外一對頻點信息位置不變,然后用低通濾波器適當?shù)靥崛∑揭坪蟮念l譜信息,最后對提取到的頻譜信息做二維傅里葉逆變換處理,可以得到TA和TB,即

式中:符號F?1表示對圖像做二維傅里葉逆變換,RA和RB分別為平移重組后TA和TB的頻譜圖,X表示低通濾波器,最終按照調(diào)制照明時散斑的橫向排列,將TA和TB融合為一空間分辨率為N×M的目標圖像。

余弦編碼復(fù)用高空間分辨率關(guān)聯(lián)成像的基本流程如圖2所示,余弦編碼矩陣與Hadamard基圖案復(fù)用來產(chǎn)生調(diào)制信息,按照與成像物體(圖2中object)空間結(jié)構(gòu)相同的方向?qū)φ{(diào)制信息進行排列,構(gòu)成高空間分辨率的調(diào)制照明散斑(圖2中的MP);照明散斑對成像物體進行調(diào)制,單像素探測器收集調(diào)制照明光與物體相互作用后的總能量,然后利用線性迭代算法復(fù)原出物體的混疊圖像T;對T做二維傅里葉變換,得到混疊圖像的頻譜圖F,對頻譜進行重組,從而獲得不同編碼信息對應(yīng)的低空間分辨率物體頻譜信息(圖2中RA和RB),利用低通濾波器(圖2中X)提取合適的頻譜信息,再進行傅里葉逆變換,從而實現(xiàn)對低空間分辨率物體圖像的重構(gòu)(圖2中TA和TB),最后按照調(diào)制照明散斑排列的順序獲得高空間分辨率物體圖像(圖2中RI)。

圖2 余弦編碼復(fù)用高分辨率關(guān)聯(lián)成像的基本流程Fig.2 Flow chart of the cosine encoded multiplexing high-resolution ghost imaging

3 仿真驗證

上節(jié)理論分析了余弦編碼復(fù)用高空間分辨率關(guān)聯(lián)成像方法,本節(jié)將通過數(shù)值仿真來驗證此方法的有效性。設(shè)成像物體的空間分辨率為256×512,按照所提方法可以將其分為兩個空間分辨率為256×256的子圖像來處理,其中P1對應(yīng)物體的(1:256,1:256)部分,P2對應(yīng)物體的(1:256,257:512)部分。因此,Hadamard基圖案和余弦編碼矩陣EP1與EP2的維數(shù)設(shè)置為256×256。其中,(3)式中編碼矩陣的對比度a和偏置b均設(shè)置為0.5,fx1=0,fy1=128,fx2=128,fy2=0,生成的編碼矩陣及其對應(yīng)的傅里葉頻譜如圖3所示?？梢钥闯鼍幋a矩陣EP1的頻譜在空間坐標(129,1)位置上有一處明顯的沖激點,編碼矩陣EP2的頻譜在空間坐標(1,129)位置上具有明顯沖激點,這種特征有利于后續(xù)的高空間分辨率圖像解碼重構(gòu)。

圖3 余弦編碼矩陣及其傅里葉頻譜。EP1與EP2為余弦編碼空間表現(xiàn)形式;F(EP1)和F(EP2)為編碼矩陣的傅里葉頻譜Fig.3 Cosine encoded matrices and their Fourier spectra.EP1 and EP2 demonstrate the spatial representations of the cosine encoded matrices.F(EP1)and F(EP2)demonstrate the Fourier spectra of the employed cosine encoded matrices

本方法在應(yīng)用傅里葉逆變換復(fù)原低空間分辨率子圖像過程中,選取的頻譜范圍對重構(gòu)的圖像質(zhì)量影響較大。這里選用理想低通濾波器來提取子圖像的頻譜,理想低通濾波器X(fx,fy)定義為

式中:R(fx,fy)為傅里葉域中空間坐標(fx,fy)到中心原點的距離;R0為濾波半徑,濾波半徑的大小表征了物體頻譜信息對應(yīng)的截至頻率。為定量評估重構(gòu)圖像的質(zhì)量,分別采用均方差(MSE,EMS)和結(jié)構(gòu)相似性指數(shù)[21](SSIM,MSSI)來對重構(gòu)圖像進行評估,分別表示為

式中:c(x,y)為空間位置(x,y)上的灰度值;o(x,y)表示仿真使用的原始圖像空間位置(x,y)上的灰度值,MSE大小反應(yīng)了復(fù)原圖像與原始圖像之間的誤差,其值越小表示復(fù)原圖像質(zhì)量越高,越接近原始圖像;μc和μo為重建圖像和原始圖像的像素點均值;σc和σo分別是重建圖像和原始圖像的標準差;σco是協(xié)方差;C1和C2是穩(wěn)定弱分母常數(shù),C1=(K1L)2,C2=(K2L)2,一般情況下K1=0.01,K2=0.03,L=255。SSIM衡量圖片的失真程度,當兩張圖像一模一樣時,SSIM的值等于1。

本文方法在實現(xiàn)對空間分辨率為256×512的物體成像時,全采樣所需的余弦編碼復(fù)用調(diào)制散斑數(shù)量為65536(等于256×256),65536個復(fù)用的調(diào)制散斑與空間分辨率為256×512的成像物體[如圖4(a)]相互作用后產(chǎn)生對應(yīng)的65536個強度值,聯(lián)合調(diào)制散斑和強度值,利用線性迭代算法復(fù)原出成像物體的混疊圖像,如圖4(b)所示?；殳B圖像的空間分辨率為256×256,混疊圖像在數(shù)學上是對高空間分辨率物體圖像抽樣所形成的低分辨率的圖像,其融合了編碼信息與空間強度信息。對混疊圖像做傅里葉變換可得其對應(yīng)的傅里葉頻譜,如圖4(c)。為更直觀地展示,將低頻信息平移到圖像的中心位置,即低頻信息在頻譜圖像的中間位置,高頻信息在頻譜圖像的邊緣位置。鑒于所構(gòu)造編碼矩陣具有的特定頻譜結(jié)構(gòu),低空間分辨率圖像的頻譜被頻移到混疊圖像傅里葉域的高頻區(qū)域,即圖4(c)空間坐標(129,1)和(129,256)以及(1,129)和(256,129)附近區(qū)域,將低空間分辨率對應(yīng)圖像的頻譜平移到頻譜圖像的中心位置,可得到兩個低空間分辨率圖像的頻譜圖,如圖4(d)、(e)所示。

圖4成像物體、復(fù)原的混疊圖像及其傅里葉頻譜。(a)空間分辨率為256×512成像物體圖像;(b)復(fù)原的256×256成像物體的混疊圖像;(c)混疊圖像的傅里葉頻譜;(d),(e)為對(c)頻移后對應(yīng)低空間分辨率子圖像的頻譜圖Fig.4 Picture of the imaged object,restored mixed image and its Fourier spectrum.(a)Image of the object with 256×512 resolution;(b)Reconstructed mixed image of the object with 256×256 resolution;(b)Fourier spectrum of the reconstructed mixed image;(d),(e)Fourier spectrum of the low spatial resolution image obtained through frequency shift from(c)

圖5 為余弦編碼復(fù)用高空間分辨率關(guān)聯(lián)成像仿真實驗結(jié)果。通過使用理想低通濾波器來分別提取圖4(d)、(e)低頻信息,然后應(yīng)用傅里葉逆變換來復(fù)原相應(yīng)的低空間分辨率物體圖像,最后將復(fù)原的兩個低空間分辨率圖像合成為一個高空間分辨率圖像。圖5(a)～(h)展示了濾波半徑分別為12、24、32、40、48、56、64、72 pixel時復(fù)原的結(jié)果。當濾波半徑較小時提取的頻譜信息較少,復(fù)原的圖像比較模糊;隨著濾波半徑的逐漸增大,提取的頻譜信息增大,復(fù)原圖像的MSE逐漸減小,誤差變小,SSIM逐漸增加,即復(fù)原的圖像和原圖的相似度提高;但是當濾波半徑大于56時[圖5(g)和(h)],可能是選取頻譜出現(xiàn)了混疊,導(dǎo)致MSE和SSIM逐漸變差?？傮w上,通過選擇合適的濾波半徑,余弦編碼復(fù)用能夠高效地實現(xiàn)對高空間分辨率物體成像。在高空間分辨率圖像復(fù)原時,本研究分別使用了迭代算法和數(shù)字圖像處理的相關(guān)方法,使得時間消耗大大降低;在主頻為2.3 GHz、i7-10875H處理器以及16 G內(nèi)存計算機上,應(yīng)用MATLAB R2015b實現(xiàn)256×512分辨率成像復(fù)原的時間消耗約為21.13 s,去除迭代算法復(fù)原混疊圖像的時間消耗,僅從混疊圖像中解碼重構(gòu)兩個低空間分辨率圖像時間消耗約為0.82 s。

圖5 利用理想低通濾波器復(fù)原的結(jié)果,濾波半徑分別為:(a)12 pixel;(b)24 pixel;(c)32 pixel;(d)40 pixel;(e)48 pixel;(f)56 pixel;(g)64 pixel;(h)72 pixelFig.5 Recovered results with ideal low-pass filters.The filter radius are(a)12 pixel;(b)24 pixel;(c)32 pixel;(d)40 pixel;(e)48 pixel;(f)56 pixel;(g)64 pixel;(h)72 pixel,respectively

利用多個低空間分辨率的余弦編碼散斑復(fù)用為高空間分辨率調(diào)制散斑對物體調(diào)制成像時,調(diào)制散斑數(shù)量與低空間分辨率的調(diào)制散斑相同,能夠大幅降低調(diào)制散斑數(shù)據(jù),以提高高空間分辨率關(guān)聯(lián)成像的成像效率。利用本方法對成像物體劃分的塊數(shù)增多,頻譜尺寸減小,低通濾波器的半徑受到進一步限制,復(fù)原的高空間分辨率圖像質(zhì)量會有一定程度的降低。下面將256×256空間分辨率的物體分別等分為2部分和4部分來處理,并定量計算重構(gòu)圖像的質(zhì)量(SSIM和MSE)。結(jié)果如圖6所示,其中圖6(a)為仿真使用的空間分辨率為256×256的測試圖像。當將256×256分辨率物體等分為4部分時,每部分的空間分辨率為128×128,此時全采樣所需的散斑數(shù)量為16384,為了防止頻譜重疊,最大的濾波半徑為32 pixel。圖6(b)～(e)分別是濾波半徑為8、16、24、32 pixel時的成像結(jié)果,可以看出隨著濾波半徑的增加,復(fù)用的圖像質(zhì)量逐漸增加,最大的SSIM為0.7112,最小的MSE為0.0032。相應(yīng)地,當將256×256分辨率物體等分為2部分時,每部分的空間分辨率為128×256,此時全采樣所需的散斑數(shù)量為32768,為了防止頻譜重疊,此時最大的濾波半徑為64 pixel。圖6(f)～(j)分別為濾波半徑為32、40、48、56、64 pixel時的成像結(jié)果,同樣可以看出,隨著濾波半徑的增加,復(fù)用的圖像質(zhì)量逐漸增加,即使在濾波半徑為32 pixel時,復(fù)原圖像的SSIM和MSE分別達到了0.7698和0.0025,比等分為4部分時的高。仿真表明,等分數(shù)量的增加雖然降低了采樣次數(shù),但也降低了高空間分辨率成像質(zhì)量,因此在應(yīng)用本研究所提出方法時需兼顧成像分辨率(劃分的塊數(shù))和成像質(zhì)量。

圖6 應(yīng)用所提出方法將成像物體等分為4份和2份時的成像對比。(a)仿真使用的測試圖像;(b)～(e)將成像物體等分為4份,濾波半徑分別為8、16、24、32 pixel時利用理想低通濾波器復(fù)原的結(jié)果;(f)～(j)將成像物體等分為2份,濾波半徑分別為32、40、48、56、64 pixel時利用理想低通濾波器復(fù)原的結(jié)果Fig.6 Imaging comparison when the imaging object is equally divided into 4 and 2 parts by using the proposed method.(a)Testing imaging object;(b)～(e)Recovered results using ideal low-pass filters with the filter radius at 8,16,24 and 32 pixel respectively when the imaging object is equally divided into 4 parts;(f)～(j)Recovered results with the filter radius at 32,40,48,56 and 64 pixel respectively using ideal low-pass filters when the imaging object is equally divided into 2 parts

4 討論與結(jié)論

關(guān)聯(lián)成像的機理決定了其成像效率和空間分辨率相互矛盾。橫向空間分辨率越高則所需的調(diào)制信息越多,在調(diào)制頻率一定的情況下,使得采集時間過長,無法滿足實際需求。本研究提出了一種余弦編碼復(fù)用高空間分辨率關(guān)聯(lián)成像技術(shù)。通過構(gòu)造多個低空間分辨率的余弦編碼散斑復(fù)用為高空間分辨率調(diào)制散斑對成像物體進行調(diào)制照明,聯(lián)合調(diào)制信號以及調(diào)制信號與成像物體相互作用后的強度值,利用線性迭代算法復(fù)原出成像場景的混疊圖像;應(yīng)用傅里葉變換將混疊圖像從空間域轉(zhuǎn)換到頻域,通過頻域中的頻移操作分別獲取低空間分辨率物體頻譜,然后利用理想低通濾波器來選取合適頻譜以重構(gòu)低空間分辨率物體圖像,最后按照低空間分辨率調(diào)制散斑的排列次序合成高空間分辨率物體圖像。數(shù)值仿真驗證了所提出方法的有效性,所實現(xiàn)的余弦編碼復(fù)用高空間分辨率關(guān)聯(lián)成像大幅降低了調(diào)制散斑數(shù)量,減少了在線采樣時間,同時在低空間分辨率圖像重構(gòu)時規(guī)避了傳統(tǒng)優(yōu)化算法,進一步降低了圖像重構(gòu)的時間消耗,在生物醫(yī)學等要求高空間分辨率且對時間苛刻的領(lǐng)域具有一定的應(yīng)用價值。