陳紀(jì)輝 王峰 理玉龍 張興 姚科 關(guān)贊洋 劉祥明

1) (中國(guó)工程物理研究院,激光聚變研究中心,綿陽(yáng) 621900)

2) (復(fù)旦大學(xué),核物理與離子束應(yīng)用教育部重點(diǎn)實(shí)驗(yàn)室,上海 200433)

1 引言

慣性約束聚變(inertial confinement fusion,ICF)是實(shí)現(xiàn)可控核聚變反應(yīng)的一條重要途徑[1].對(duì)ICF 中高溫等離子體的X 射線輻射的空間分辨研究,可以提供ICF 靶內(nèi)由于做功和能量(電子熱傳導(dǎo)、輻射熱傳導(dǎo))輸運(yùn)導(dǎo)致的流體狀態(tài)的空間演化信息[2,3].目前,用于X 射線成像的設(shè)備主要有:X 射線針孔相機(jī)[4]、KB 反射鏡[5,6]、Wolter 顯微鏡[7]以及球面彎晶[8].上述成像診斷儀器采用的都是直接成像方式,即通過(guò)光透射、反射、衍射等傳輸過(guò)程,在探測(cè)平面直接反映成像目標(biāo)的光輻射空間分布的高分辨二維圖像.然而,通過(guò)這種方式采集到的二維圖像丟失了沿探測(cè)方向上的強(qiáng)度分布信息,其包含的光強(qiáng)分布圖像實(shí)際上是在一定軸距范圍內(nèi)光強(qiáng)的沿軸積分疊加.綜上所述,現(xiàn)有的診斷手段還無(wú)法通過(guò)單次測(cè)量實(shí)現(xiàn)對(duì)射線源的層析成像.

全息技術(shù)是實(shí)現(xiàn)對(duì)物體三維成像的一條有效途徑.相干全息術(shù)利用相干性極強(qiáng)的光源(通常是激光)作為參考光,在照明物體后,將物光波前的振幅和相位以干涉條紋的形式記錄為全息圖.之后再通過(guò)適當(dāng)?shù)姆绞?從這幅包含物光“全部信息”的全息圖中重建出物體的三維形貌[9].傳統(tǒng)的相干全息術(shù)對(duì)于照明光源的相干性要求極高,而非相干光源在自然界中更為普遍,這在很大程度上限制了其應(yīng)用范圍.1961 年,Mertz 和Young[10]將全息技術(shù)推廣至非相干光領(lǐng)域.他們將針孔相機(jī)中的透光小孔替換為菲涅耳波帶片,使得不同位置的光透過(guò)波帶片后在底片上產(chǎn)生放大倍數(shù)不同、中心位置不同的波帶片形貌投影.這些疊加投影形成的全息圖反映了原光源的空間信息,通過(guò)適當(dāng)?shù)闹亟ǚ椒ň湍苓€原出其在各個(gè)位置的光強(qiáng)分布.在非相干全息理論被提出后不久,核醫(yī)學(xué)[11,12]、ICF 成像診斷[13,14]以及天文學(xué)[15]等領(lǐng)域都采用了這種具有高收光效率的“編碼成像”方式.近年來(lái),非相干全息技術(shù)被廣泛應(yīng)用于無(wú)透鏡成像領(lǐng)域,并衍生出了菲涅耳相機(jī)的理念[16,17].2020 年,Wu 等[18]在反演菲涅耳相機(jī)底片(非相干全息圖)的過(guò)程中引入了信號(hào)傳輸中壓縮感知優(yōu)化算法,獲得了高分辨率、低背景噪聲的二維重建圖像.2023 年,Soltau 等[19]使用菲涅耳波帶片作為編碼孔徑完成了高光通量的X 射線全場(chǎng)熒光成像.值得注意的是,到目前為止的非相干全息成像(編碼成像)研究大都集中于二維光源圖像的采集與重建,而早期利用Gabor 波帶片進(jìn)行三維編碼成像的嘗試,受到了非聚焦像的嚴(yán)重干擾[20,21].

本文在前人工作的基礎(chǔ)上,介紹了一種基于非相干全息成像機(jī)制與壓縮感知重建算法,針對(duì)微尺寸X 射線源的層析成像技術(shù);并通過(guò)模擬三維光源的照明全息記錄和層析重建闡明了該診斷手段的應(yīng)用價(jià)值.模擬結(jié)果顯示了在選用的實(shí)驗(yàn)參數(shù)體系下,非相干全息成像系統(tǒng)能對(duì)軸向尺寸為16 mm的光源進(jìn)行層析成像,以4 mm 為步長(zhǎng)的重建圖像中幾乎不存在來(lái)自其他層面信號(hào)的干擾.

2 非相干全息層析成像機(jī)制

2.1 非相干全息成像理論模型

非相干全息的成像過(guò)程如圖1 所示,空間體光源可以看作是數(shù)個(gè)點(diǎn)光源的集合,位于不同空間位置的點(diǎn)光源獨(dú)立透過(guò)波帶片,并在探測(cè)器上留下不同中心位置和放大倍數(shù)的波帶片投影(當(dāng)像距遠(yuǎn)小于波帶片焦距時(shí),不考慮衍射帶來(lái)的影響),各個(gè)投影疊加形成的非相干全息圖由下式表示:

圖1 非相干全息成像示意圖Fig.1.Procedure of incoherent holography.

其中,h是探測(cè)器上獲取的光強(qiáng)分布;(x′,y′,z′)表示物方區(qū)域的坐標(biāo);t'是波帶片投影在探測(cè)面上的強(qiáng)度分布;o是待測(cè)體光源的光強(qiáng)分布,M是成像縮放倍數(shù),且該光源延軸向的固有長(zhǎng)度為z0.這里考慮將正弦波帶片作為“編碼孔徑”,其在探測(cè)面上的投影可以寫(xiě)為

其中,r*是經(jīng)過(guò)投影放大后的波帶片最內(nèi)環(huán)半徑:

式中,d是波帶片距離探測(cè)器之間的距離(像距),z是光源與波帶片之間的距離(物距),r1是波帶片孔徑的最內(nèi)環(huán)半徑.根據(jù)(1)式的表述,對(duì)于二維光源,非相干全息圖可以視為光源平面上的光強(qiáng)分布與該物距下產(chǎn)生的波帶片投影的卷積;而三維光源產(chǎn)生的全息圖則是通過(guò)所有物距平面產(chǎn)生的二維全息圖累加形成.通常,在空間域較為復(fù)雜的卷積和退卷積運(yùn)算可以轉(zhuǎn)化為在頻域中更為簡(jiǎn)潔的表述,即

式中,H是h的二維傅里葉變換,而T和O分別是t和o在給定物距上的二維傅里葉變換.使用傳輸算子A表征物空間光強(qiáng)分布與非相干全息圖之間的映射:

其中,F-1代表二維傅里葉逆變換算符,Ti代表(4)式中T沿Z方向的第i個(gè)分量.Fi意為取三維矩陣中沿Z方向第i個(gè)分量進(jìn)行二維傅里葉變換.因此,h與o之間的運(yùn)算關(guān)系可以表示為

2.2 壓縮感知重建模型

上述三維非相干全息成像過(guò)程可以看作是一個(gè)對(duì)數(shù)據(jù)的壓縮采樣過(guò)程,即使用具有(Nx×Ny)個(gè)探測(cè)單元的探測(cè)器去采集體量為(Nx×Ny×Nz)的目標(biāo)信號(hào).因此,在反解二維全息圖重建三維光源光強(qiáng)分布的過(guò)程中,不可避免地會(huì)遇到方程個(gè)數(shù)遠(yuǎn)小于未知數(shù)的情況.根據(jù)壓縮感知理論[22,23],這種求解“欠定方程組”的問(wèn)題需要引入正則條件來(lái)對(duì)方程的解空間范圍施加約束,繼而求出既滿足全息映射過(guò)程,又符合正則約束的最優(yōu)解.本文三維圖像信號(hào)通過(guò)最小化目標(biāo)函數(shù)重建:

其中,第1 項(xiàng)為殘差項(xiàng),表征反解出的三維圖像經(jīng)過(guò)非相干全息投影后與實(shí)際采集到的全息圖之間的差異;第2 項(xiàng)為正則項(xiàng),對(duì)于微尺寸X 射線源的三維圖像,使用它的全變分范數(shù)(total variation norm,TV norm)作為正則項(xiàng)來(lái)表征圖像信號(hào)的稀疏性;β為正則參數(shù),用于平衡目標(biāo)函數(shù)中殘差項(xiàng)和正則項(xiàng)的比重.采用兩步迭代收縮閾值算法[24](TWIST)對(duì)(7)式進(jìn)行求解,以獲得不同物距層面上的光強(qiáng)分布.

3 模擬成像與重建結(jié)果

為了驗(yàn)證非相干全息層析成像技術(shù)的可行性,模擬了針對(duì)微尺寸X 射線源的成像過(guò)程.這里假設(shè)待測(cè)信號(hào)是分布在三維空間的I,C,F 三個(gè)獨(dú)立的字母形貌光源,其軸向間距均為8 mm,水平間距為0.20 mm;每個(gè)獨(dú)立字母光源占據(jù)的平面空間約為0.16 mm × 0.20 mm;在距離字母光源“C”50 mm 的位置上放置一片正弦波帶片作為收光孔徑,其最內(nèi)環(huán)半徑為0.08 mm,環(huán)帶數(shù)為256 環(huán),波帶片總半徑為1.28 mm;探測(cè)器位于距離波帶片遠(yuǎn)離光源一面500 mm 的位置上,用于采集字母光源產(chǎn)生的非相干全息圖.其中,探測(cè)器的探測(cè)面尺寸為16 mm × 16 mm,像素點(diǎn)數(shù)目為2000 × 2000,單位像素點(diǎn)尺寸為8 μm × 8 μm;字母光源輻射的X 射線能量設(shè)定為15 keV,上述波帶片對(duì)該能量X 射線的焦距約為7.75×104mm.在本模擬成像系統(tǒng)的物像距和波帶片參數(shù)下,其衍射效應(yīng)可忽略不計(jì),符合非相干全息的幾何光學(xué)成像機(jī)制;模擬成像的三維示意圖以及在上述參數(shù)下產(chǎn)生的非相干全息圖如圖2 所示,圖3 是模擬成像系統(tǒng)二維側(cè)視圖.

圖2 針對(duì)微尺寸X 射線源的成像示意圖與二維全息圖Fig.2.Imaging system for microscale X-ray source and the corresponding 2D hologram.

圖3 模擬成像系統(tǒng)二維側(cè)視圖Fig.3.2D lateral view of the simulative imaging system.

本文采用早期Gabor 波帶片編碼成像實(shí)驗(yàn)用到的背傳輸算法[19],維納濾波算法[21]以及前文理論部分介紹的壓縮感知模型三種重建算法求解全息圖.圖4 展示了對(duì)本次模擬生成的非相干全息圖的層析重建結(jié)果,可知上述3 種重建方法都能在確定的物距上反解出相應(yīng)的字母光源,這說(shuō)明非相干全息技術(shù)具有三維成像的潛力.對(duì)比這3 種算法的重建結(jié)果可以發(fā)現(xiàn): 背傳輸算法計(jì)算出光源信號(hào)的邊緣還存在著大量串聯(lián)噪聲,這會(huì)在實(shí)際應(yīng)用中對(duì)信號(hào)識(shí)別造成影響;維納濾波算法雖然能在一定程度上識(shí)別出待測(cè)圖像信號(hào),但由于重建出的圖像對(duì)比度過(guò)低,使得目標(biāo)源強(qiáng)分布混雜于背景噪聲之中.不能有效地凸顯重建有效信號(hào)的強(qiáng)度特征;而基于壓縮感知算法的重建結(jié)果顯示: 其非信號(hào)區(qū)域的串聯(lián)噪聲顯著降低,光源物距平面上的光強(qiáng)分布基本上都集中在信號(hào)區(qū)域.除了存在光源信號(hào)的物距平面,我們還額外計(jì)算了距字母光源前后4 mm平面上的光源分布.從圖4(b)和圖4(c)的重建結(jié)果可以看出,利用背傳輸算法和維納濾波算法計(jì)算得到的非信號(hào)平面上,仍然存在著大量來(lái)自其他平面上光源信號(hào)的非聚焦贗像.如前所述,在壓縮感知重建的模型((7)式)中,引入TV norm 作為正則項(xiàng),以此來(lái)表征重建圖像在梯度域的稀疏度.由于上述非聚焦贗像無(wú)法在特定的物平面內(nèi)重建成為具有明顯邊界的清晰圖像,這些在梯度域中不稀疏的干擾項(xiàng),會(huì)在求解(7)式極小值的過(guò)程中不斷被消除.因此,通過(guò)壓縮感知算法的重建,非信號(hào)平面內(nèi)的背景噪聲會(huì)隨著算法迭代次數(shù)的增加而有效地被抑制.在圖4(d)—(f)中,分別展示了TWIST 算法迭代100 次、500 次與2000 次的重建結(jié)果.

為了定量表征壓縮感知重建后的圖像質(zhì)量,使用峰值信噪比(peak signal-to-noise ratio,PSNR,單位dB)作為評(píng)價(jià)指標(biāo),其表達(dá)式定義為

其中,Omax是圖像單個(gè)像素點(diǎn)的最大數(shù)值.MSE是均方誤差,用于反映原始數(shù)據(jù)和重建數(shù)據(jù)之間的差異,其表達(dá)式可以寫(xiě)作:

其中,o和o*分別是源強(qiáng)分布矩陣和重建后的強(qiáng)度分布矩陣,Nx代表圖像矩陣的總行數(shù),Ny代表圖像矩陣的總列數(shù),Nz代表圖像矩陣的總頁(yè)數(shù).由(9)式可以看出,MSE 值越小說(shuō)明重建后的圖像與源光強(qiáng)分布越接近.而基于MSE 計(jì)算得出的PSNR,其數(shù)值越大表明圖像的失真程度越小.圖5展示了TWIST 算法迭代次數(shù)與重建圖像PSNR值之間的關(guān)系.可以看出,PSNR 數(shù)值在算法運(yùn)行初期經(jīng)歷微弱的震蕩后快速上升,這表明圖像質(zhì)量在該迭代區(qū)間(Iterations ≤ 200)內(nèi)被迅速優(yōu)化,并在經(jīng)歷大量迭代次數(shù)后逐漸趨于平穩(wěn).

圖5 PSNR 與重建迭代次數(shù)間的關(guān)系Fig.5.Relationship between PSNR and reconstruction iterations.

除了對(duì)非相干全息技術(shù)應(yīng)用于層析成像的可行性研究,本文還研究了不同參數(shù)的波帶片對(duì)成像分辨水平的影響.與波帶片聚焦成像類(lèi)似,非相干全息成像的水平分辨率也與波帶片的最內(nèi)環(huán)半徑成反比(在一定環(huán)帶數(shù)目的條件下).只是在非相干全息成像中用于參考的是波帶片投影的最內(nèi)環(huán)半徑,它不僅與波帶片編碼孔徑相關(guān),還受到物像距關(guān)系的調(diào)制.圖6 展示了選用兩種不同最內(nèi)環(huán)半徑波帶片對(duì)前述的字母光源進(jìn)行非相干全息模擬成像的重建結(jié)果對(duì)比.基于圖5 所示的圖像質(zhì)量隨迭代次數(shù)變化關(guān)系,為了使用較好且穩(wěn)定的圖像重建質(zhì)量進(jìn)行對(duì)比,兩種參數(shù)的波帶片重建結(jié)果均經(jīng)過(guò)TWIST 算法的2000 次迭代計(jì)算得出.重建結(jié)果包含兩個(gè)部分;一部分是光源在物距為54,58和62 mm 平面處的二維重建結(jié)果,另一部分是位于58 mm 平面上虛線標(biāo)識(shí)區(qū)域內(nèi)的一維光強(qiáng)分布情況.用作對(duì)比的兩種波帶片的總環(huán)帶數(shù)均為256 環(huán),最內(nèi)環(huán)半徑分別為0.08 和0.06 mm.通過(guò)對(duì)比圖6(a)和圖6(b)的重建結(jié)果可以發(fā)現(xiàn),在選用最內(nèi)環(huán)半徑為0.06 mm 的波帶片進(jìn)行成像后,在光源信號(hào)平面內(nèi)的分辨率得到了明顯提高,具體表現(xiàn)在一維光強(qiáng)分布中的信號(hào)上升和下降沿所經(jīng)過(guò)的距離變窄.其中,圖6(a)的虛線區(qū)域由上至下共經(jīng)歷4 個(gè)邊沿,其距離分別為80,64,64 和72 μm,而在圖6(b)中相應(yīng)的邊沿距離為40,40,32 和40 μm.但值得注意的是,使用小內(nèi)環(huán)半徑的波帶片后在光源信號(hào)平面的非信號(hào)區(qū)域和非信號(hào)平面內(nèi)重建產(chǎn)生了來(lái)自其他平面光源的非聚焦像,這意味著其軸向分辨率降低,非相干全息成像的重建算法無(wú)法正確復(fù)現(xiàn)以4 mm 作為軸間距的層析影像.因此,非相干全息成像的軸向和水平面內(nèi)分辨率相互制約,在實(shí)際應(yīng)用中需要根據(jù)需求選取特定參數(shù)的波帶片作為編碼孔徑.

圖6 使用不同最內(nèi)環(huán)半徑波帶片進(jìn)行模擬成像的分辨水平對(duì)比.在距離波帶片54,58 和62 mm 平面上的二維重建結(jié)果,以及在58 mm 平面內(nèi)虛線標(biāo)識(shí)區(qū)域的一維光強(qiáng)分布情況 (a) 0.08 mm;(b) 0.06 mmFig.6.Comparisons of resolution level in simulative imaging when applied FZP with different innermost radius.The 2D reconstruction result at the following objective depth: 54,58 and 62 mm,and the 1D intensity distribution of the dotted line area in 58 mm plane: (a) 0.08 mm;(b) 0.06 mm.

4 結(jié)論

本文研究了非相干全息成像機(jī)制以及三維圖像重建算法.為了驗(yàn)證該技術(shù)的可行性,數(shù)值模擬了微尺寸體光源的非相干全息成像過(guò)程,并使用壓縮感知算法對(duì)其產(chǎn)生的全息圖進(jìn)行逐層反解,得到了體光源距離探測(cè)器不同深度位置上的二維光強(qiáng)分布情況.結(jié)果表明: 在文中所例舉的參數(shù)體系下,非相干全息技術(shù)能夠?qū)ξ⒊叽珞w光源進(jìn)行層析成像;相較于傳統(tǒng)的背傳輸算法和維納濾波算法,基于壓縮感知算法的重建結(jié)果能有效的消除來(lái)自不同物距上光源的非聚焦贗像干擾.綜上,該項(xiàng)技術(shù)有望用于對(duì)ICF 中黑腔等離子體分布進(jìn)行層析成像,完善ICF 領(lǐng)域中現(xiàn)有的成像診斷手段.