針對微尺寸X 射線源的非相干全息層析成像*

陳紀輝 王峰 理玉龍 張興 姚科 關贊洋 劉祥明

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

2) (復旦大學,核物理與離子束應用教育部重點實驗室,上海 200433)

1 引言

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

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

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

2 非相干全息層析成像機制

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

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

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

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

其中,r*是經過投影放大后的波帶片最內環半徑:

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

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

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

2.2 壓縮感知重建模型

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

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

3 模擬成像與重建結果

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

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

圖3 模擬成像系統二維側視圖Fig.3.2D lateral view of the simulative imaging system.

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

為了定量表征壓縮感知重建后的圖像質量,使用峰值信噪比(peak signal-to-noise ratio,PSNR,單位dB)作為評價指標,其表達式定義為

其中,Omax是圖像單個像素點的最大數值.MSE是均方誤差,用于反映原始數據和重建數據之間的差異,其表達式可以寫作:

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

圖5 PSNR 與重建迭代次數間的關系Fig.5.Relationship between PSNR and reconstruction iterations.

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

圖6 使用不同最內環半徑波帶片進行模擬成像的分辨水平對比.在距離波帶片54,58 和62 mm 平面上的二維重建結果,以及在58 mm 平面內虛線標識區域的一維光強分布情況 (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 結論

本文研究了非相干全息成像機制以及三維圖像重建算法.為了驗證該技術的可行性,數值模擬了微尺寸體光源的非相干全息成像過程,并使用壓縮感知算法對其產生的全息圖進行逐層反解,得到了體光源距離探測器不同深度位置上的二維光強分布情況.結果表明: 在文中所例舉的參數體系下,非相干全息技術能夠對微尺寸體光源進行層析成像;相較于傳統的背傳輸算法和維納濾波算法,基于壓縮感知算法的重建結果能有效的消除來自不同物距上光源的非聚焦贗像干擾.綜上,該項技術有望用于對ICF 中黑腔等離子體分布進行層析成像,完善ICF 領域中現有的成像診斷手段.