姚春霞 何其利 張錦 付天宇 吳朝 王山峰黃萬霞 袁清習(xí) 劉鵬 王研? 張凱?

1) (中國科學(xué)院高能物理研究所, 北京同步輻射裝置, X 射線光學(xué)與技術(shù)實(shí)驗(yàn)室, 北京 100049)

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

3) (中國科學(xué)技術(shù)大學(xué), 國家同步輻射實(shí)驗(yàn)室, 合肥 230029)

1 引 言

X 射線自被倫琴發(fā)現(xiàn)以來就被廣泛應(yīng)用于醫(yī)學(xué)影像診斷以及材料檢測等領(lǐng)域.但是傳統(tǒng)的X 射線成像技術(shù)主要是基于吸收襯度成像機(jī)制, 對由碳、氫、氧等元素構(gòu)成的弱吸收物質(zhì)的成像質(zhì)量不佳.為了克服上述困難, 發(fā)展基于X 射線相位襯度成像理論的新型X 射線成像技術(shù)成為近年來的研究熱點(diǎn).其中, X 射線光柵微分相位襯度成像技術(shù)由于對X 射線光源的相干性要求較低, 可以與實(shí)驗(yàn)室光源相結(jié)合實(shí)現(xiàn)弱吸收物質(zhì)的大視場、高分辨相位襯度成像, 因此被認(rèn)為是最有可能在臨床醫(yī)學(xué)中獲得應(yīng)用的一種X 射線相位襯度成像方法,近年來獲得了廣泛的發(fā)展[1?3].但是傳統(tǒng)X 射線光柵微分相位襯度成像技術(shù)需要采用分析光柵作為空間濾波器使用.一方面, 分析光柵會阻擋X 射線, 降低了整個成像系統(tǒng)X 射線的利用效率.另一方面, 分析光柵需要進(jìn)行步進(jìn)掃描來獲取樣品的多張投影圖像才能提取樣品的吸收、折射和散射信息, 因此存在樣品曝光時間長、輻射劑量高等問題.上述這些問題的存在都限制了X 射線光柵微分相位襯度成像技術(shù)的應(yīng)用范圍, 特別是生物醫(yī)學(xué)成像等對輻射劑量敏感的研究領(lǐng)域.因此如何解決上述X 射線光柵微分相位襯度成像技術(shù)發(fā)展過程中存在的X 射線的利用效率低和樣品輻射劑量高的問題, 從而實(shí)現(xiàn)快速、低劑量的X 射線光柵微分相位襯度成像技術(shù), 成為近年來X 射線相位襯度成像領(lǐng)域發(fā)展的重要方向.

實(shí)現(xiàn)上述快速、低劑量X 射線光柵微分相位襯度成像技術(shù)的發(fā)展需求可以通過兩種途徑來實(shí)現(xiàn).第一種途徑是改造光柵的結(jié)構(gòu), 避免分析光柵的步進(jìn)掃描過程從而減少樣品的輻射損傷; 第二種途徑是發(fā)展基于免分析光柵X 射線相位襯度成像系統(tǒng)的一次曝光相位襯度成像技術(shù).在對第一種途徑的探索中, Endrizzi 等[4]、Fu 等[5]、Wei 等[6]以及Ge 等[7]分別通過改造光柵的結(jié)構(gòu)實(shí)現(xiàn)了樣品一次曝光成像.但是這種方法僅僅是避免了分析光柵的步進(jìn)掃描過程, 減少了樣品的輻射損傷.其仍需要使用分析光柵, 因此無法解決分析光柵對X 射線的吸收導(dǎo)致的整個成像系統(tǒng)X 射線利用率低的問題.此外經(jīng)過改造的光柵結(jié)構(gòu)還降低了成像裝置在縱向上的分辨率.而基于免分析光柵X 射線相位襯度成像系統(tǒng)的一次曝光相位襯度成像技術(shù)則完全避免了分析光柵的使用.其不需要分析光柵進(jìn)行步進(jìn)掃描來進(jìn)行樣品的數(shù)據(jù)采集, 只需要進(jìn)行樣品的一次投影數(shù)據(jù)采集就能夠提取樣品的吸收、折射和散射信息.該方法可以有效地提高整個成像系統(tǒng)X 射線利用率的同時解決傳統(tǒng)X 射線光柵微分相位襯度成像技術(shù)中樣品曝光時間長、輻射劑量高的問題.因此通過第二條途徑實(shí)現(xiàn)上述快速、低劑量的X 射線光柵微分相位襯度成像技術(shù)的發(fā)展需求更為有效.Balles[8]、Marathe[9]以及Wen[10]、Bennett[11]等驗(yàn)證了第二種途徑的可行性,此后基于免分析光柵X 射線相位襯度成像系統(tǒng)的一次曝光成像理論的發(fā)展受到了廣泛的關(guān)注.但是目前基于該成像理論的樣品信息提取算法還不完善.例如Balles 等[8]和Marathe 等[9]使用的相位步進(jìn)信息提取算法(phase-stepping algorithm, PS算法)雖然能夠得到準(zhǔn)確的樣品吸收、折射和散射信息, 但PS 算法需要滿足光柵自成像期等于探測器像素尺寸的整數(shù)倍的限制條件.然而這一限制條件在實(shí)際實(shí)驗(yàn)過程中是難以嚴(yán)格滿足的, 其光柵自成像周期和探測器像素尺寸之間的差異會嚴(yán)重影響樣品信息提取的精度, 進(jìn)而影響X 射線相位襯度成像中定量化分析的準(zhǔn)確性和有效性; 而Wen等[10]和Bennett 等[11]所使用的樣品信息提取算法是對單張投影圖直接進(jìn)行空間諧波分析, 與PS算法相比較, 該算法會顯著降低樣品圖像信息的空間分辨率.

針對基于免分析光柵X 射線相位襯度成像系統(tǒng)中一次曝光成像理論在樣品信息提取算法上存在的上述問題, 本文開展了基于免分析光柵X 射線相位襯度成像系統(tǒng)的樣品信息提取算法的研究,提出一種基于超定方程理論的最小二乘法的樣品信息提取算法.該算法的適用性并不受光柵周期和探測器單個像素尺寸的匹配性限制.此外, 該算法在免分析光柵X 射線相位襯度成像系統(tǒng)中可以獲得與傳統(tǒng)PS 算法具有相同空間分辨率的樣品的吸收、折射和散射圖像信息.

2 免分析光柵一次曝光相位襯度成像理論與樣品信息提取算法

2.1 基于免分析光柵X 射線相位襯度成像系統(tǒng)的一次曝光相位襯度成像原理

基于免分析光柵X 射線相位襯度成像系統(tǒng)裝置如圖1 所示, 主要包括滿足一定相干條件的X 射線源、相位光柵 G1(π 相位光柵或 π /2 相位光柵)、以及探測器像素尺寸為 b ×b 的高分辨探測器.本文中為便于描述, 在免分析光柵X 射線相位襯度成像系統(tǒng)裝置中相位光柵 G1選取為 π 相位光柵.如圖1(a)所示, 由于Talbot 效應(yīng)相位光柵 G1在Talbot 距離D 處會形成 G1光柵的自成像條紋其中自成像條紋的周期p 為 π 相位光柵 G1周期p1的一半[12].高分辨探測器則位于 G1的自成像位置處, 可以直接對自成像條紋進(jìn)行圖像采集.

圖1 基于免分析光柵X 射線相位襯度成像系統(tǒng)的裝置成像示意圖Fig.1.Schematic imaging diagram of the phase contrast imaging system without analyser grating.

在實(shí)驗(yàn)過程中, 如圖2 所示當(dāng)高分辨率探測器的像素尺寸b 等于 p /N 時(N 為正整數(shù)), 可使用高分辨探測器分別采集一張樣品的投影圖像和背景的投影圖像.此后將得到的樣品投影圖像沿 xg方向每N 列抽取一列像素組成一張新投影圖像, 則可獲得N 張新投影圖像, 分別定義為 I1, I2, ···,IN.當(dāng)N 足夠大時, 上述抽取方法獲得的N 張新投影圖像中相同像素位置處的光強(qiáng)可以組成圖2所示的光強(qiáng)變化曲線.而當(dāng) N =4 時可以獲得I1(黑色), I2(紅色), I3(綠色), I4(黃色)共四張新投影圖像.

圖2 由采集到的投影圖像組成的N 張新投影圖像中同一像素點(diǎn)的光強(qiáng)隨X 射線的偏折角 ψ =xg/D 的變化曲線(角度曲線)Fig.2.Plot of the new images’ intensity oscillation (shifting curve) of single pixel as a function of the deviation angle ψ =xg/D.

這些新投影圖像等效于傳統(tǒng)的含分析光柵的X 射線光柵微分相位襯度成像裝置中, 使用占空比為 1 /N 的分析光柵沿 xg方向做步進(jìn)掃描運(yùn)動, 并以 p /N 為步進(jìn)間距進(jìn)行N 次采集得到的投影圖像.因此傳統(tǒng)的含分析光柵的X 射線光柵微分相位襯度成像技術(shù)和圖1 所示的免分析光柵的相位襯度成像技術(shù)的原理是類似的, 都是對自成像條紋進(jìn)行N 點(diǎn)采樣.根據(jù)上述分析, 則圖2 中有樣品時的新投影圖像的光強(qiáng)可以表示為[13,14]:

其中 ? 表示卷積; f (ψ;x,y) 表示樣品對X 射線的角度調(diào)制函數(shù), 包含吸收、折射和散射效應(yīng);是無樣品時的N 張新投影圖像的光強(qiáng)所形成的背景角度曲線, η =1, 2, ··· ,N 表示一組新投影圖像的序號;與傳統(tǒng)X 射線光柵微分相位襯度成像系統(tǒng)中通過吸收光柵做步進(jìn)掃描所獲得的系統(tǒng)位移曲線類似; Iη(ψ;x,y) 是有樣品時的N 張新投影圖像中相同像素點(diǎn) (x,y) 上的光強(qiáng)變化曲線.此處與 Iη(ψ;x,y) 均可以采用余弦函數(shù)來擬合[14?16].可表示為

其中 ψ0表示背景角度曲線的初始角度;=(Smax+Smin)/2 表示背景角度曲線的平均值; Smax和 Smin分別表示背景角度曲線的最大值和最小值;V0=(Smax?Smin)/(Smax+Smin)表示背景角度曲線的可見度.而樣品對X 射線的角度調(diào)制函數(shù)f(ψ;x,y)可表示為[17]

其中M 代表樣品對X 射線的吸收; θ 代表X 射線經(jīng)過樣品后的折射角; σ2表示X 射線經(jīng)過樣品后散射角度高斯分布的方差.根據(jù)(1)式—(3)式, 則免分析光柵X 射線相位襯度成像系統(tǒng)中探測器上的樣品光強(qiáng)方程可表示為

在實(shí)際實(shí)驗(yàn)過程中, 相位光柵自成像周期p 是難以滿足為探測器像素尺寸b 的N 倍的要求的,因此為了保證樣品信息提取的精度, 需要在角度曲線余弦擬合中考慮探測器的像素尺寸b 和光柵自成像周期p 之間的差值, 這里令探測器像素尺寸b 與光柵自成像的周期的 1 /N 的差值所對應(yīng)的角度曲線上的相位差為 α , 即

2.2 信息提取算法

將(5)式代入成像方程式(4)可以得到N 張樣品光強(qiáng)圖的表達(dá)為

令

以及

聯(lián)立(6)式中的N 個方程, 則(6)式又可以寫為

進(jìn)而簡寫為

其中,

(10)式中的A 是由采樣點(diǎn)數(shù)N 以及探測器的像素尺寸b 以及相位光柵自成像的周期p 決定的常數(shù)矩陣.由于 N ＞3 , 故簡化后的樣品投影圖像的表達(dá)式(10)式是一個超定線性方程組, 為了解(10)式,可以將其轉(zhuǎn)換為一個最小二乘優(yōu)化問題:

其解為

(14)式中的 p inv(A) 表示矩陣A 的偽逆.將該解中的 x1, x2和 x3代入(8)式可以得到:

同理, 當(dāng)沒有放置樣品時, 可以將探測器上獲得的背景投影圖像分為N 張光強(qiáng)圖像:

其解為

由于無樣品時吸收 Mbg, 折射 θbg和散射 σbg都為0,故可得

因此, 結(jié)合背景扣除[18,19], 并將(7)式分別代入(15)式和(18)式得扣除背景后樣品的吸收、折射和散射信息提取公式分別為

其中 a rctan 是 tan 的反三角函數(shù).

下文中為描述方便, 將上述樣品信息提取算法稱為基于超定方程理論的最小二乘法的樣品信息提取算法(least-squares algorithm, LS 算法).而傳統(tǒng)的光柵微分相位襯度成像中N 步等間距平移分析光柵的相位步進(jìn)法對應(yīng)的信息提取公式為[15,19,20]

其中

3 仿真實(shí)驗(yàn)與結(jié)果

3.1 系統(tǒng)參數(shù)設(shè)計(jì)

模擬實(shí)驗(yàn)裝置如圖1 所示, X 射線的能量為E =12 keV.選用的 π 相位光柵周期 p1為4.80 μm,探測器放置在一階Talbot 距離D 為27.88 mm處.實(shí)驗(yàn)中設(shè)計(jì)了一個圓柱體和三張厚度不同的紙作為實(shí)驗(yàn)樣品, 如圖3 所示.其中圓柱體的材質(zhì)為PMMA, 其高為 r =0.61 mm , 對X 射線具有吸收和折射效應(yīng), 其復(fù)折射率虛部 β 為 1.79×10?9、實(shí)部 δ 為 1.85×10?6; 而紙張對X 射線只有散射效應(yīng),其散射寬度與紙的層數(shù)c 有關(guān), c 層紙對X 射線的散射角度方差為 c σ2, 實(shí)驗(yàn)中 σ2取 1.67×10?11rad2.模擬實(shí)驗(yàn)過程為: X 射線穿過圓柱體時, 吸收效應(yīng)會引起X 射線光強(qiáng)的衰減, 同時折射效應(yīng)會使得X 射線偏離原來的入射方向, 折射角為 θ.隨后穿過圓柱體的X 射線會照射在紙上, 此時X 射線會以折射角為中心發(fā)生小角度的前向高斯散射, 此后小角度前向散射光會被成像探測器上的像素單元所接收.當(dāng)X 射線穿過不同厚度的紙張時, 小角度前向高斯散射光的散射寬度不同(散射寬度為因此會造成探測器上單個像素單元接收到的X 射線光子數(shù)發(fā)生變化, 即X 射線光強(qiáng)發(fā)生變化.上述樣品仿真模型的具體描述可參見參考文獻(xiàn)[17,21].

圖3 模擬所用的樣品示意圖.樣品由直徑為0.461 mm PMMA 圓柱狀和PMMA 后不同厚度(0—3 層)的紙組成,其中PMMA 圓柱狀對X 光具有吸收和折射效應(yīng), 而紙張對X 光只有散射效應(yīng)Fig.3.Schematic diagram of the simulation sample.The simulation sample has a 0.46 mm diameter PMMA cylinder that combined refraction and absorption effects.The PMMA cylinder overlie paper layers (0?3 layers) that exhibit the scattering effects.

在上述模擬過程中, 假設(shè)樣品厚度遠(yuǎn)小于成像距離, 則可使用基于菲涅爾衍射理論來計(jì)算X 射線經(jīng)過樣品及 G1光柵后在探測器上的光強(qiáng)分布圖像[19,22,23].此外, 為了比較探測器的像素尺寸與光柵自成像周期的匹配差對樣品信息提取精度的影響, 可以分別使用兩個不同像素尺寸的X 射線高分辨成像探測器來對模擬樣品進(jìn)行數(shù)據(jù)采集.這里假設(shè)兩個探測器的總視場均為 r ×r , 且探測器的像素尺寸分別為 b1=0.60 μm (α =0)和 b2=0.66 μm(α /=0).模擬得到的樣品投影圖像如圖4(a),(c)所示.其中圖4(a)對應(yīng)的探測器像素尺寸為0.60 μm,圖4(c)對應(yīng)的探測器像素尺寸為0.66 μm, 從圖4(a)右下角的放大圖圖4(b)可以看到, 當(dāng)探測器的像素尺寸與光柵自成像的周期相匹配時, 探測器采集到的樣品投影圖像也為周期性分布, 且投影圖像的周期與光柵自成像周期相等.而從圖4(c)右下角的放大圖圖4(d)可以觀察到, 當(dāng)探測器的像素尺寸與光柵自成像周期不匹配時, 匹配差的存在使投影圖中出現(xiàn)了莫爾條紋.

使用(19)式描述的基于超定方程理論的樣品信息提取算法將圖3 中所示的不同探測器采集到的樣品投影圖像以 N ·b 為抽取周期進(jìn)行列像素抽取, 則可以分別重排為 N =4 張樣品投影圖像.如圖5 所示, 其中探測器的像素尺寸為0.60 μm 時采集的樣品投影圖像所組合出的4 張新的投影圖像如圖5(a)—(d)所示, 分別為樣品處于角度曲線的谷位、上坡位、峰位、下坡位時的投影圖像.從圖5(b)和圖5(d)中可以觀察到, 兩幅圖像在圓柱體的左右兩邊緣位置分別出現(xiàn)亮暗的現(xiàn)象, 其代表了折射角方向的不同, 該現(xiàn)象和傳統(tǒng)含分析光柵的X 射線光柵微分相位襯度成像裝置中分析光柵做4 次步進(jìn)掃描后在位移曲線上坡、下坡位得到的樣品投影圖像是吻合的[14].此外圖5(a)和圖5(c)分別為樣品位于角度曲線的谷位和峰位時的投影圖像, 此時由于角度曲線是對稱分布的, 因此圖5(a)和圖5(c)樣品邊界處的光強(qiáng)分布也是對稱的.

圖4 探測器上的樣品投影圖 (a) 探測器像素尺寸為0.60 μm (α = 0)時模擬樣品的投影圖像; (b) 圖(a)的局部(紅框內(nèi))放大圖; (c) 探測器像素尺寸為0.66 μm (α ≠ 0)時模擬樣品的投影圖像; (d) 圖(c)的局部放大圖Fig.4.The projective images of the simulation sample:(a) The projective image of the simulation sample with pixel size of 0.60 μm (α = 0); (b) local enlarged drawing of Fig.(a); (c) the projective image of the simulation sample with pixel size of 0.66 μm (α ≠ 0); (d) local enlarged drawing of Fig.(c).

當(dāng)采用探測器的像素尺寸為0.66 μm 時, 采集到的樣品投影數(shù)據(jù)經(jīng)過抽取和重排后, 組合出新的4 張樣品投影圖像如圖5(e)—(h)所示, 其4 張投影圖像在 xg方向上均存在周期性的明暗條紋,這是因?yàn)楫?dāng)探測器的像素尺寸與光柵自成像周期為非整數(shù)倍時, 按周期抽取后的同一投影圖像中的多列像素不再分布在角度曲線的同一位置, 從而出現(xiàn)了莫爾條紋.

3.2 樣品的信息提取

利用PS 以及LS 樣品信息提取算法, 對圖5中的兩組數(shù)據(jù)分別進(jìn)行樣品的吸收、折射和散射信息的提取.當(dāng)采用的探測器像素尺寸為0.60 μm時, 樣品信息提取結(jié)果如圖6 所示.其中, 圖6(a)—(c)為使用PS 算法提取得到的樣品的吸收、折射和散射信息, 圖6(d)—(f)為使用LS 算法提取得到的樣品的吸收、折射和散射信息.可以看到圖6(a)中由于圓柱體樣品是厚度不均勻的, 中心區(qū)域的光強(qiáng)吸收較大, 使得其在中心區(qū)域樣品的吸收信息更強(qiáng).而在折射信息圖6(b)中, 圓柱邊界區(qū)域的折射角大且兩側(cè)折射角方向相反, 使得折射圖像兩側(cè)一亮一暗.圖6(c)所示為樣品的散射信息, 由于紙的厚度從上到下依次變薄, 其散射系數(shù)也不斷減小,因此圖6(c)中可以看到從上到下, 圖像依次由亮變暗.

圖5 由原始的樣品投影圖組成的4 張新圖像 (a)—(d) 探測器像素尺寸為0.60 μm (α = 0)時所得到的4 張新投影圖;(e)—(h) 探測器像素尺寸為0.66 μm ()時所得到的4 張新投影圖Fig.5.Four new images extracted from one original projective image: (a)?d) Contain four new images with the pixel size of 0.60 μm (α = 0); (e)?(h) contain four new images with the pixel size of 0.66 μm (α ≠ 0).

為了更好的分析圖6 中的樣品信息提取的準(zhǔn)確性, 沿圖6 中虛線所示方向, 分別獲取了樣品的吸收信息、折射信息和散射信息的強(qiáng)度分布曲線.如圖6(g)—(i)所示, 當(dāng)探測器的像素尺寸與光柵自成像的周期相匹配時, PS 算法與LS 算法所得的吸收、折射和散射基本一致, 且和相應(yīng)的理論值能夠較好地吻合.

當(dāng)模擬中所用的探測器的像素尺寸為0.66 μm時, 樣品信息提取結(jié)果如圖7 所示.其中, 圖7(a)—(c)使用PS 算法提取得到的樣品的吸收、折射和散射信息, 圖7(d)—(f)為使用LS 算法提取得到的樣品的吸收、折射和散射信息.圖7(a)—(f)與圖6(a)—(f)基本一致, 不過當(dāng)使用PS 算法進(jìn)行樣品信息提取時, PS 算法提取的樣品吸收、折射和散射圖像中出現(xiàn)了與光柵方向平行的條狀偽影, 如圖7(a)—(c)所示.同樣的, 圖7(a)—(f)中虛線所示的樣品的吸收信息、折射信息和散射信息的強(qiáng)度分布曲線的比較結(jié)果如圖7(g)—(i)所示, 可以發(fā)現(xiàn)當(dāng)探測器的像素尺寸與光柵自成像的周期不匹配時, PS 算法提取的吸收、折射和散射像中出現(xiàn)了與光柵方向平行的條狀偽影, 所得到的吸收、折射和散射信息與理論值差異較大, 而LS 算法所得的吸收、折射與散射信息則與理論值吻合較好.另外, 由圖7(h)和圖7(i)圖中可以看到, 折射信息和散射信息在圓柱樣品邊界區(qū)域的提取值與理論值之間存在比較明顯的偏差.

圖6 探測器像素尺寸為0.60 μm 時模擬樣品的信息提取結(jié)果 (a)—(c) PS 算法時提取的吸收、折射和散射信息; (d)—(f) LS算法時提取的吸收、折射和散射信息; (g)—(i) 在虛線位置處(PS 算法(綠色)和LS 算法(紫色))的提取的吸收、折射和散射信息與理論值(藍(lán)色)的對比圖Fig.6.Sample information extracted with pixel size of 0.60 μm: (a)?(c) Depict the absorption, refraction and scattering information of the simulated sample obtained by PS algorithm; (d)?(f) depict the absorption, refraction and scattering information obtained by LS algorithm; (g)?(i) are profiles of the absorption, refraction and scattering images extracted by LS (purple) and PS(green) at the dotted lines in (a)?(f), as well as the theoretical values (blue).

圖7 探測器像素尺寸為0.66 μm 時模擬用樣品的信息提取結(jié)果 (a)—(c) PS 算法時提取的吸收、折射和散射信息; (d)—(f) LS算法時提取的吸收、折射和散射信息; (g)—(i) 在虛線位置處(PS 算法(綠色)和LS 算法(紫色))的提取的吸收、折射和散射信息與理論值(藍(lán)色)的對比圖Fig.7.Sample information extracted with pixel size of 0.66 μm: (a)?(c) Depict the absorption, refraction and scattering information of the simulated sample obtained by PS algorithm; (d)?(f) depict the absorption, refraction and scattering information obtained by LS algorithm; (g)?(i) are profiles of the absorption, refraction and scattering images extracted by LS (purple) and PS(green) at the dotted lines in Fig.(a)?(f), as well as the theoretical values (blue).

為了進(jìn)一步對LS 算法和PS 算法兩種算法進(jìn)行定量比較, 這里使用平均絕對誤差(mean absolute error, MAE)來評估樣品信息提取值e 與理論值 e0的誤差情況, 其平均絕對誤差可表示為

其中 m ×n 為樣品信息值的總像素數(shù).根據(jù)(22)式計(jì)算出圓柱樣品區(qū)域內(nèi)的MAE 誤差的結(jié)果如表1所列.在探測器像素尺寸與光柵自成像的周期匹配時, LS 算法樣品信息提取值的MAE 誤差與PS 算法樣品信息提取值的MAE 誤差基本相等, 且誤差均較小, 說明兩種算法提取的樣品信息值和理論值均能較好地吻合.當(dāng)探測器的像素尺寸與光柵自成像的周期不匹配時, 探測器的像素尺寸為0.66 μm時, 兩種算法樣品信息提取值的MAE 誤差均有所增加, 但是LS 算法樣品信息提取值的MAE 誤差比PS 算法樣品信息提取值的MAE 誤差小將近1 個量級, 說明在周期不匹配的情況下LS 算法比PS 算法提取得到的樣品信息更加準(zhǔn)確.

表1 PS 和LS 算法的理論值和提取值的平均絕對誤差Table 1.Mean absolute error of theoretical and extracted information by PS and LS algorithm.

此外, 圖6(a),(c),(d),(f)和圖7(a),(c),(d),(f)中樣品邊界處均出現(xiàn)了兩條縱向的明暗相間的條紋狀偽影, 這主要是由兩個方面的原因造成的: 1) 本文采用了菲涅耳傳播理論對X 射線穿過樣品的物理過程進(jìn)行模擬, 隨著X 射線傳播距離的增加,在Talbot 距離D 處樣品邊界處出現(xiàn)了相位的二階導(dǎo)數(shù)信息, 即衍射效應(yīng), 而LS 算法或PS 算法無法提取出樣品的相位二階導(dǎo)數(shù)信息, 因此引起了樣品邊界處的明暗相間的條紋狀偽影; 2) 在樣品信息提取算法中, 需要對角度曲線 S′(ψ) 進(jìn)行余弦擬合,而實(shí)際角度曲線在上述位置存在較強(qiáng)高次諧波, 并非完全可以用單一余弦函數(shù)來表述, 故擬合過程存在誤差, 這同樣會使得樣品邊界部分信息提取出現(xiàn)上述明暗相間的條紋狀偽影[17,24].

為了進(jìn)一步研究探測器的像素尺寸與光柵自成像周期不匹配時對信息提取的準(zhǔn)確性的影響.在圖1 所示成像裝置中換用多個不同像素尺寸(0.54,0.57, 0.60, 0.63, 0.66 μm)的探測器進(jìn)行樣品的投影數(shù)據(jù)的獲取, 此后利用PS 算法以及LS 算法分別提取樣品的吸收、折射和散射信息, 并按照圖6和圖7 中虛線位置處提取樣品的吸收信息、折射信息和散射信息的強(qiáng)度分布曲線, 其結(jié)果如圖8所示.

如圖8 所示, 當(dāng)探測器的像素尺寸與光柵自成像的周期0.60 μm 之間的偏差越大, 即在 α 的絕對值越大時, PS 算法提取的吸收、折射和散射像中得到的與光柵方向平行的條狀偽影就越強(qiáng).而LS算法的提取結(jié)果受到的探測器的像素尺寸的影響相對較小.考慮到兩種算法都是假設(shè)搖擺曲線為余弦曲線, 并進(jìn)行N 點(diǎn)采樣后進(jìn)行擬合.因此, 其信息分離準(zhǔn)確性很大程度上取決于其余弦假設(shè)的可靠性.由于PS 算法要求擬合所用余弦函數(shù)周期必須等于實(shí)際角度曲線的周期, 在免分析光柵相位襯度成像系統(tǒng)中PS 算法擬合角度曲線所用的余弦函數(shù)的周期是N 個探測器像素寬度之和 N ·b , 當(dāng)α=0時, 擬合角度曲線所用的余弦函數(shù)的周期N ·b與實(shí)際角度曲線周期p 相等, 故PS 算法與LS 算法都能得到精確的信息提取結(jié)果, 但時, 實(shí)際角度曲線的周期p 不再等于N 個探測器像素寬度之和 N ·b , 即擬合所用余弦函數(shù)周期與角度曲線周期不匹配, 所以會影響擬合的可靠性,并且隨著 α 絕對值的增大, 余弦函數(shù)周期與角度曲線周期的差異越大, 提取得到的樣品信息的失真更加嚴(yán)重.而LS 算法中余弦函數(shù)的周期可以是任意的, 不需要是一組探測器的像素尺寸之和 N ·b , 因此隨著像素大小的變化, 能夠較為穩(wěn)定的分離得到吸收、折射和散射像.

圖8 不同像素尺寸的探測器提取樣品的吸收、折射和散射信息 (a)—(c) 利用PS 算法提取的樣品吸收、折射和散射信息的強(qiáng)度曲線; (d)—(f) 利用LS 算法提取的樣品吸收、折射和散射信息的強(qiáng)度曲線Fig.8.The absorption, refraction, and scattering information of the simulated sample with different pixel sizes: (a)?(c) Depict the profiles of extracted absorption, refraction and scattering images with different pixel sizes by PS algorithm; (d)?(f) depict the profiles of extracted absorption, refraction and scattering images with different pixel sizes by LS algorithm.

此外, 在上述免分析光柵一次曝光相位襯度成像系統(tǒng)中, LS 算法與PS 算法對單張投影圖像的處理方法相同, 均是將探測器采集到的單張?jiān)纪队皥D按像素抽取為N 張新投影圖像, 故LS 算法與PS 算法所得的樣品吸收、折射和散射圖像信息的空間分辨率是相同的, 均為探測器采集到的單張投影圖空間分辨率的 1 /N 倍.

4 結(jié) 論

基于免分析光柵X 射線相位襯度成像系統(tǒng)的一次曝光成像理論由于不需要使用分析光柵, 從而避免了分析光柵對X 射線的衰減, 提升了整個成像系統(tǒng)中X 射線的利用效率, 同時一次曝光也減小了對樣品的輻射損傷, 因此基于免分析光柵X 射線相位襯度成像系統(tǒng)的一次曝光成像理論的發(fā)展是近年來光柵成像技術(shù)發(fā)展的方向, 但是在該成像理論中傳統(tǒng)的相位步進(jìn)算法(PS 算法)對光柵周期與探測器像素尺寸的匹配提出了嚴(yán)格的限制條件, 而空間諧波分析算法會嚴(yán)重降低樣品吸收、折射和散射圖像信息的空間分辨.這些問題的存在限制了一次曝光成像理論的應(yīng)用和發(fā)展.針對上述問題本文提出了適用于上述成像系統(tǒng)的基于超定方程理論的最小二乘法的樣品信息提取算法(LS 算法).該算法可以獲得與傳統(tǒng)PS 算法具有相同空間分辨率的樣品吸收、折射和散射信息, 此外在光柵周期與探測器像素尺寸存在匹配偏差時依然可以得到可靠的信息提取結(jié)果.因此, LS 算法的使用有助于降低上述成像系統(tǒng)對實(shí)驗(yàn)裝置精度的要求, 增加成像裝置的實(shí)用性, 從而使得快速、低劑量光柵相襯成像技術(shù)在臨床醫(yī)學(xué)、材料等領(lǐng)域更具實(shí)用前景.