夏慧娟 吳衍青 張 磊 孫元鶴 邰仁忠

1（中國科學(xué)院上海應(yīng)用物理研究所 上海 201800）

2（中國科學(xué)院上海高等研究院 上海 201210）

3（中國科學(xué)院大學(xué) 北京 100049）

X射線成像技術(shù)在醫(yī)療、安檢、工業(yè)探傷、無損檢測等領(lǐng)域中具有舉足輕重的地位。由于X射線能量較高，所以實(shí)際應(yīng)用中往往利用X射線的熒光效應(yīng)將X射線轉(zhuǎn)化為可見光進(jìn)行間接成像。尤其對于高能X射線同步輻射光源，間接成像探測器比直接探測壽命長，接收能量范圍高，而且容易實(shí)現(xiàn)。但是由于在轉(zhuǎn)換過程中受到各種因素，尤其是高折射率的影響，X射線間接成像質(zhì)量有待提高。

在成像技術(shù)中，經(jīng)常通過圖像的反卷積算法來消除成像系統(tǒng)的各種不理想因素，如離焦、運(yùn)動(dòng)模糊、光學(xué)介質(zhì)畸變等不利影響。這些反卷積算法在X射線成像技術(shù)中也有不少應(yīng)用。這些算法分為非盲去卷積算法和盲去卷積算法兩大類［1-2］。前者更適合精確恢復(fù)圖像。目前常用的非盲去卷積算法有近鄰法、逆濾波法、維納濾波法、線性最小平方法等線性算法［3-4］；Jansson-Van Citter法、非線性最小平方法等約束迭代法，最大似然法、最大后驗(yàn)概率法等統(tǒng)計(jì)迭代方法；稀疏自適應(yīng)網(wǎng)絡(luò)、完全卷積學(xué)習(xí)網(wǎng)絡(luò)等新方法［5-7］。這些算法計(jì)算量普遍較大，且都需要事先估計(jì)成像系統(tǒng)的點(diǎn)擴(kuò)散函數(shù)（Point Spread Function，PSF）才能得到復(fù)原圖像；盲去卷積算法一般也需要預(yù)估PSF。PSF是衡量光學(xué)系統(tǒng)成像性能的數(shù)學(xué)工具，是恢復(fù)圖像過程中消除成像系統(tǒng)各種不利因素影響的關(guān)鍵。PSF估計(jì)準(zhǔn)確與否直接影響了圖像恢復(fù)的效果。其估計(jì)方法大致分為實(shí)驗(yàn)測量估計(jì)和計(jì)算構(gòu)建估計(jì)兩類。

在X射線閃爍體成像系統(tǒng)的PSF估計(jì)中常用的實(shí)驗(yàn)方法有傳統(tǒng)邊緣法、斜邊法、微孔測量法、微聚焦測量法等［8-9］。前兩種方法用于估計(jì)成像系統(tǒng)的噪聲圖像邊緣角度存在不準(zhǔn)確性。在微孔測量法中，微孔孔徑過大會(huì)導(dǎo)致PSF估計(jì)誤差大；過小的微孔只能在薄片上加工，其厚度不足以阻擋X射線的高能輻射，同樣導(dǎo)致較大誤差；微聚焦測量法要求將光斑聚焦到一個(gè)像素點(diǎn)以內(nèi)，一般都在1 μm以內(nèi)。對測量設(shè)備要求高，操作難度大；而且聚焦到微米級大小，光斑的能量會(huì)很大，探測器非常容易飽和，總之實(shí)現(xiàn)難度較大。一般來說，實(shí)測的PSF噪聲不易控制，這對圖像反卷積不利。在計(jì)算構(gòu)建估計(jì)方法中，首先構(gòu)建某種數(shù)學(xué)模型，然后根據(jù)各種形式的反饋信息修正該模型的參數(shù)，最終得到一個(gè)復(fù)原誤差可以接受的PSF模型。如基于小波分析的高斯PSF估計(jì)［10］，基于值平面上梯度矢量分布的PSF估計(jì)［11］等。總體來說，計(jì)算構(gòu)建類方法不會(huì)有噪聲等因素干擾，但目前該類方法普遍比較復(fù)雜，并且只適用于比較簡單成像系統(tǒng)的PSF獲取，如運(yùn)動(dòng)模糊圖像和散焦模糊圖像的PSF估計(jì)，對于復(fù)雜PSF的估計(jì)能力較弱，方法缺乏普適性。

本文采用微孔-反卷積測量法測量PSF，使用直徑為20 μm的微孔進(jìn)行測量，并用抗噪聲效果較好的Lucy-Richardson算法進(jìn)行反卷積［12］，得到信噪比理想的PSF。由于該算法不能滿足PSF中央?yún)^(qū)域的高精度要求，我們利用非相干PSF光信息分布的非定域性重建了PSF中央?yún)^(qū)域［13］。本文以一種X射線新型閃爍體成像系統(tǒng)為例，利用重建后的PSF獲得了較為理想的圖像反卷積效果，驗(yàn)證了該方法的有效性。

1 理論分析

1.1 非相干光PSF測量中的問題

通過有限離散化處理之后，PSF可以看成一個(gè)幾何點(diǎn)（δ函數(shù)）在其像平面上許多有限尺寸像點(diǎn)的疊加。

由式（1）可知，已知微孔的物分布和像分布后，通過反卷積得到我們需要的PSF［14］。因此，測量PSF實(shí)驗(yàn)方案是采用大小合適的微孔，成像之后做反卷積，得到PSF。這應(yīng)該是一種比較可行的方法。

常用的反卷積函數(shù)包括維納濾波器（Wiener）復(fù)原，正則化（Regularized）濾波器復(fù)原和Lucy-Richardson方法復(fù)原。在知道清晰圖像和噪聲頻率特性的情況下，維納濾波器復(fù)原被看作是一種有效的圖像復(fù)原技術(shù)。正則化濾波器也是在知道噪聲信息的情況下有效的反卷積方法。但是，我們實(shí)驗(yàn)中圖像頻率和噪聲信息都是未知的，而Lucy-Richardson濾波器方法假設(shè)噪聲服從泊松分布，基于貝葉斯理論使產(chǎn)生圖像的似然性達(dá)到最大值，不需要已知噪聲信息和圖像頻率特征，保留了原始數(shù)據(jù)的計(jì)數(shù)性和非負(fù)性，抗噪聲性能比較好，每次迭代的結(jié)果都是非負(fù)的，并且對PSF中的小錯(cuò)誤具有很強(qiáng)的魯棒性。因此我們選擇Lucy-Richardson方法復(fù)原［15］。

Lucy-Richardson方法雖然具有很多優(yōu)點(diǎn)，但是它也存在一定的誤差，特別是區(qū)域小，邊緣銳利的情況下，會(huì)存在很大的誤差［16］。在實(shí)驗(yàn)中測量的PSF，由于-1級衍射光受非線性效應(yīng)影響小，使-1級信息失真小，另外-1級區(qū)域面積大，強(qiáng)度遠(yuǎn)小于0級信息的強(qiáng)度，用Lucy-Richardson方法做恢復(fù)，誤差小。Lucy-Richardson方法廣泛應(yīng)用，大部分應(yīng)用情況下，該方法都可以滿足精度要求，而我們研究中，0級包含圖像細(xì)節(jié)信息，因此我們對0級的精度要求更高。

1.2 CTF-PSF和OTF-PSF的關(guān)系及信息冗余

由衍射光學(xué)的知識(shí)可知，非相干成像的光學(xué)傳遞函數(shù)（Optical Transfer Function，OTF）等于同等透鏡條件下相干衍射受限系統(tǒng)的相干傳遞函數(shù)（Coherent Transfer Function，CTF）的自相關(guān)歸一化函數(shù)［17］。PSF的傅里葉變換就是光學(xué)系統(tǒng)的傳遞函數(shù)。CTF-PSF是指在相干光照明下成像系統(tǒng)的PSF；OTF-PSF是指在非相干光照明下成像系統(tǒng)對應(yīng)的PSF。因此：

式中：h(x0，y0；xi，yi)是由菲涅爾衍射公式計(jì)算而來的 CTF。其中：d0為物距，di為像距；(x0，y0)為物平面上一點(diǎn)坐標(biāo)，(xi，yi)為像平面上一點(diǎn)坐標(biāo)；dU′(x0，y0；x，y)是光波經(jīng)過孔徑函數(shù) P(x，y)，焦距為透鏡后的單位脈沖響應(yīng)的復(fù)振幅。將式（3）帶入式（2）中并化簡可以得到式（4）。

式中：積分號前面的相位因子不影響最終的衍射強(qiáng)度分布，在討論成像場強(qiáng)分布的時(shí)候可以近似忽略。h(x0，y0；xi，yi)做自相關(guān)并歸一化即可得到 OTF。因此我們提出了冗余信息提取的方法，冗余信息是指光傳輸過程中出現(xiàn)其他系統(tǒng)的多余信息，通過提取多余的信息，利用非相干光照明條件下光信息的非定域性質(zhì)，修正某一區(qū)域扭曲的信息。在本研究的例子中，非相干光PSF中-1級衍射區(qū)域中存在0級光信息，且誤差小。我們將設(shè)法將該0級信息提出，重新構(gòu)造0級衍射區(qū)域。利用這一點(diǎn)，將非相干衍射光場轉(zhuǎn)換到相干衍射光場，并在理論上與其對應(yīng)的相干光照明下的光場信息進(jìn)行對比修正，再轉(zhuǎn)回非相干光場，并用實(shí)際光場修正計(jì)算衍射光場；此過程不斷循環(huán)迭代，直至得到正確的0級信息。

2 X射線閃爍體成像系統(tǒng)

以一種新型的X射線閃爍成像系統(tǒng)為例，對閃爍晶體表面修飾二維光柵陣列，可以突破全內(nèi)反射的限制，提高光學(xué)信息的提取，特別是高空間頻率的光學(xué)信息［18-19］。提取抑制高頻噪聲，大大提高信噪比。然而，提取的信息會(huì)失真，導(dǎo)致多個(gè)虛擬物點(diǎn)，并且具有不同的空間分布，使最終的圖像失焦、分散和退化，因此在沒有圖像恢復(fù)的情況下不利于成像。圖1（a）是成像系統(tǒng)的對照組的模型圖。X射線照射到樣品上，X射線通過閃爍體后轉(zhuǎn)換成可見光，最后由CCD成像。圖1（b）是成像系統(tǒng)的實(shí)驗(yàn)組的模型圖，兩者的唯一區(qū)別是實(shí)驗(yàn)組的閃爍體的出射面修飾了二維光柵陣列。

圖1 對照組(a)和實(shí)驗(yàn)組(b)設(shè)備模塊圖Fig.1 Device module diagram of imaging system(a)Control group,(b)Experimental group

2.1 新型閃爍體成像系統(tǒng)信息冗余

在二維光柵閃爍成像系統(tǒng)中，非相干系統(tǒng)中PSF的-1級衍射包含相干傳遞函數(shù)的-1級衍射和0級衍射信息。將實(shí)驗(yàn)組和對照組測得的微孔進(jìn)行反卷積得到PSF，即對應(yīng)于OTF的PSF，去掉第0階中心。以理論結(jié)構(gòu)PSF的0級衍射為初始值。為了實(shí)驗(yàn)的方便，我們將PSF簡化為一個(gè)對稱函數(shù)，使模型的自相關(guān)運(yùn)算等價(jià)于卷積運(yùn)算。構(gòu)建點(diǎn)擴(kuò)散函數(shù)計(jì)算過程中，迭代次數(shù)200次，計(jì)算時(shí)間大約20 min，迭代之前與迭代之后相同位置的差值的誤差在10-5。不斷迭代，不斷添加正確的信息，直到收斂得到正確的PSF。流程圖如圖2所示。

圖2 冗余信息法構(gòu)建PSF的流程圖Fig.2 Flow chart of constructing point spread function with redundant information method

2.2 新型閃爍體成像系統(tǒng)PSF及圖像復(fù)原結(jié)果

實(shí)驗(yàn)數(shù)據(jù)由上海光源BL13WI線站獲得。實(shí)驗(yàn)時(shí)X射線能量15 keV。成像探測器的CCD像素陣列2 048×2 048，5倍鏡頭下像素分辨率為1.3 μm×1.3 μm，數(shù)值孔徑NA=0.5。對照組閃爍體為YAG：Ce，晶體尺寸?25.4 mm×0.2 mm，雙面拋光，折射率1.82，中心發(fā)光波長550 nm；實(shí)驗(yàn)組閃爍體為同樣YAG：Ce閃爍體，但在其表面修飾加工一層二維光柵，周期300 nm，直徑180 nm。所用實(shí)驗(yàn)樣品為脫水斑馬魚樣本，微孔大小20 μm。小孔的成像結(jié)果、反卷積PSF和冗余信息法得到的PSF在圖3中顯示。

圖4是分別用實(shí)驗(yàn)20 μm小孔，反卷積得到的PSF和信息冗余法得到的PSF對斑馬魚進(jìn)行圖像復(fù)原得到的效果圖。

圖3 小孔的成像結(jié)果(a)，反卷積PSF(b)，冗余信息法得到的PSF(c)Fig.3 Imaging results of pinholes(a),deconvolution PSF(b),PSF obtained by redundant information method(c)

在無參考圖像質(zhì)量評價(jià)當(dāng)中，常用的方法是在待評價(jià)圖像的基礎(chǔ)上通過低通濾波器構(gòu)建一個(gè)足夠模糊圖像（幾乎損失所有高頻圖像紋理細(xì)節(jié)）作為參考圖像；假設(shè)它與真實(shí)圖像做相同濾波后的圖像相同；然后對評價(jià)圖像和參考圖像進(jìn)行各種對比運(yùn)算，以運(yùn)算結(jié)果評價(jià)待評價(jià)圖像的質(zhì)量，Q值算法就是其中的一種。Q值算法建立在矩陣的奇異值分解理論（Singular Value Decomposition，SVD）基礎(chǔ)上。SVD在圖像處理、圖像評價(jià)和模式識(shí)別等領(lǐng)域廣泛應(yīng)用。同一幅圖像，經(jīng)某種處理前后，圖像Q值的相對大小具有很強(qiáng)的評價(jià)意義。若是處理過后，Q值增大，表示該處理促使圖像失真度增加，圖像質(zhì)量變差；若是處理過后，Q值減少，則表示該處理降低了圖像的失真程度，尤其是低頻部分，這一般是圖像的主要信息［20］。表1是斑馬魚的Q值變化，表1中Qorg、Qre分別表示對照圖的Q值和復(fù)原圖的Q值。ΔQ表示復(fù)原圖相對于對照圖的Q值變化率。

圖4 斑馬魚對照組圖片(a)，實(shí)驗(yàn)20 μm恢復(fù)斑馬魚的細(xì)節(jié)結(jié)果圖(b)，反卷積得到的PSF恢復(fù)斑馬魚的細(xì)節(jié)結(jié)果圖(c)，冗余信息提取法得到的PSF恢復(fù)斑馬魚的細(xì)節(jié)結(jié)果圖(d)Fig.4 The picture of zebrafish control group(a),image restored by experiment 20 μm(b),image restored by point spread function obtained by deconvolution(c),image restored by point spread function obtained by redundant information extraction method(d)

表1 采用三種復(fù)原方法復(fù)原斑馬魚前后Q值的變化Table 1 Changes of Q value before and after zebrafish restoration by three restoration methods

與對照組圖像相比，20 μm圓孔作為PSF、實(shí)驗(yàn)反卷積PSF和冗余信息法PSF的復(fù)原圖像噪聲水平都有所下降，但是由Q值評價(jià)結(jié)果可以看出，使用20 μm小孔復(fù)原和實(shí)驗(yàn)反卷積得到的PSF復(fù)原圖像Q值下降比較少，而冗余信息法復(fù)原圖像Q值則下降很多。所以冗余信息法做圖像復(fù)原處理降低了圖像的失真程度，尤其是低頻部分，圖像的主體信息。信息冗余法恢復(fù)的圖像信噪比提高更加明顯；而且得益于二維光柵的作用，多提取的光子對圖像高低頻信息起到了增強(qiáng)的作用，由于實(shí)驗(yàn)影響因素繁多，所以很難反卷積得到十分精確的系統(tǒng)PSF，故復(fù)原誤差較大，而冗余信息提取法可以有效地避免這一點(diǎn)，達(dá)到良好的復(fù)原效果，消除了二維光柵對成像造成的不利影響，保證對照組圖像的基本特征的前提下，增強(qiáng)了圖像的高低頻細(xì)節(jié)，提高了圖像的信噪比。因此冗余信息法復(fù)原的圖像高低頻細(xì)節(jié)都更加清晰，圖像質(zhì)量更高。

3 結(jié)語

本文以一種新型X射線閃爍體成像系統(tǒng)為例，采用信息冗余提取方法構(gòu)建PSF，成功恢復(fù)了受二維光柵影響而退化的成像。信息冗余提取方法利用非相干光照明下光學(xué)傳遞函數(shù)-1級衍射光信息冗余來提取0級光信息，不斷迭代加入正確信息，最終獲得系統(tǒng)的PSF，可以恢復(fù)復(fù)雜光學(xué)系統(tǒng)的成像。本研究從理論和實(shí)驗(yàn)上說明了信息冗余提取方法的有效性，且具有較強(qiáng)的適用性。