鮑 園 王 研 朱佩平 吳自玉

1（中國科學技術(shù)大學 國家同步輻射實驗室 合肥 230029）

2（中國科學院高能物理研究所 北京 100049）

衍射增強成像中余弦曲線擬合方法的噪聲分析

鮑 園1,2王 研2朱佩平1,2吳自玉1,2

1（中國科學技術(shù)大學 國家同步輻射實驗室 合肥 230029）

2（中國科學院高能物理研究所 北京 100049）

衍射增強成像是一種功能強大的相位襯度成像技術(shù)，其信息分離研究一直都是研究熱點。提出了一種簡便的信息分離方法——余弦曲線擬合法(Cosine Fitting Radiography, CFR)，該方法具有簡便、樣品所受輻射劑量低等優(yōu)點。利用同步輻射實驗和蒙特卡羅仿真方法，從角度信號成像模型出發(fā)解釋了此方法，同時從噪聲傳遞角度分析了該方法的性能和給出定量的評價指標，建立了衍射增強成像和光柵微分干涉成像之間更直觀的聯(lián)系。對理解衍射增強成像系統(tǒng)的噪聲行為和提高信噪比具有一定的參考意義。

衍射增強成像，噪聲分析，余弦曲線擬合，蒙特卡羅仿真

衍射增強成像(Diffraction Enhanced Imaging, DEI)，也叫基于分析晶體的成像方法(Analyzerbased Imaging, ABI)，可以對主要由碳、氫、氧等輕元素構(gòu)成的軟組織成像，同時獲得吸收、折射和超小角散射信息。這種成像方法最初由Forster[1]在1980年提出，但在最近十幾年才成為相襯成像領(lǐng)域的研究熱點。衍射增強成像這個術(shù)語是由Chapman[2]在1997年引入，他將搖擺曲線線性近似，由兩個半腰像，提取出樣品的吸收像和折射像。隨后，DEI的發(fā)展主要沿著兩個方向進行：(1) 將搖擺曲線認為是角度信號響應(yīng)函數(shù)，探測器接收的像就可以表達成樣品角度函數(shù)和本征搖擺曲線的卷積形式；(2) 基于全波檢波法，對比加入樣品前后的搖擺曲線，分別計算其零階、一階、二階矩，得到相應(yīng)的吸收、折射和散射信息。Chapman的方法屬于前者，高斯擬合法[3]以及Oltulu[4]、Wernick[5]的方法屬于后者。

Chapman的方法采用線性近似，搖擺曲線在兩個腰位處的二階矩均為零，因此不能提取散射信息。為彌補這個缺點，Wernick等[5]提出多圖法(Multiple-image Radiography, MIR)，利用輻射輸運理論(Radiative Transport Theory)[6]建模X射線光場傳播，證明了多圖法可以真實地分離和提純物體的吸收、折射和散射信息，并且對大搖擺角處由散射引起的圖像退化有免疫作用。Rigon等[7-8]提出了一個可以替代多圖法的方法，利用三張圖提取三種信息，但是存在一些限制。胡春紅等[3]提出用高斯曲線近似搖擺曲線，但該方法需要較多的圖像采集次數(shù)。我們課題組根據(jù)Chapman的思路，將搖擺曲線視為角度信號響應(yīng)函數(shù)，并采用余弦函數(shù)擬合(Cosine Fitting Radiography, CFR)[9]本征搖擺曲線，該方法只需三張圖就可以獲得能夠與MIR相比較的結(jié)果。

DEI不僅為放射線照相術(shù)引入了折射像和散射像，而且觸發(fā)了其他相襯成像方法的發(fā)展。尤其難能可貴的是，作為一個經(jīng)典模型，它的結(jié)論可以擴展到所有基于準直器-分析器模型的成像方法中。雖然DEI是一種功能強大的相位襯度成像技術(shù)，但是仍然有很多問題亟待解決，比如與這個方法有關(guān)的定量性能分析，Brankov等[10]就分析了Chapman法的噪聲特性。與此同時，基于光柵干涉成像的噪聲分析也有大量工作[11-16]報道，Revol[11]發(fā)展了光柵微分干涉成像(Grating-based Phase Contrast Imaging, GDPCI)中隨機誤差的定量描述，Weber[14]引入矩陣形式的最小二乘方算法，同時計算吸收、折射、散射各自的方差以及它們之間的協(xié)方差?？紤]到衍射增強成像和光柵干涉成像有很多相似之處，比如：它們都描述物體的折射角，即相位的一階導數(shù)；散射像均來源于物體引起的可見度衰減；而且噪聲模型相近，因此我們可以借鑒光柵干涉成像中的思想，轉(zhuǎn)而分析衍射增強成像的噪聲特性。本文從角度信號響應(yīng)成像模型出發(fā)，全面分析了DEI中三種信息的噪聲特性，證明了我們提出的方法的正確性，為提高成像系統(tǒng)的靈敏度提供了定量的依據(jù)，同時也建立了DEI和GDPCI之間更直觀的聯(lián)系。

1 余弦曲線擬合方法的基本原理

衍射增強成像原理如圖1所示。

圖1 衍射增強成像裝置及其搖擺曲線Fig.1 Experimental set up for diffraction enhanced imaging and its rocking curve.

余弦曲線擬合方法的思路是：將角度信號響應(yīng)函數(shù)視為余弦曲線，并用其擬合本征搖擺曲線：

式中，IθA為探測器每個像素上接收到的光子數(shù)；I0為入射光平均強度（光通量）；A為探測器面積；t為采集一幅圖時樣品的曝光時間；表示余弦函數(shù)在一個周期內(nèi)的平均值；表示本征搖擺曲線的等效可見度；),(yxV為樣品搖擺曲線的可見度。為保證擬合結(jié)果的唯一性，只用峰位和腰位來擬合，效果如圖2所示。圖2中，上坡位(Up Load)和下坡位(Down Load)分別對應(yīng)搖擺曲線兩個半腰位置，峰位(Peak)對應(yīng)搖擺曲線最大值。

圖2 余弦函數(shù)擬合效果Fig.2 Fit of the intrinsic RC by the cosine function.

將樣品放置在兩塊晶體的光路中，物體對光線的偏折等價于分析晶體的偏轉(zhuǎn)：

式中，(,)V x y為樣品搖擺曲線的等效可見度；0V為無樣品時的可見度?？梢姸葟?V變?yōu)?,)V x y的退化可看作是搖擺曲線對于散射信號的響應(yīng)，其關(guān)系為：

2 噪聲特性分析

我們考慮兩種主要的噪聲源，分別是光源光子數(shù)的波動和分析晶體的位置偏差。由此引起的總的誤差可表示為各自誤差之和：

2.1光源光子數(shù)的波動

由于光子數(shù)的波動，最終提取的吸收系數(shù)、折射角以及散射方差也會圍繞其理論值波動。不失一般性，我們假設(shè)光源光子數(shù)波動符合泊松分布，運用誤差傳遞理論[17]，計算這三種信息各自的方差：

從式(8)-(10)可以看出，折射像和散射像的噪聲受接收到的光子數(shù)、樣品的吸收以及搖擺曲線的等效可見度影響，而吸收像的噪聲只與接收到的光子數(shù)有關(guān)。

如果折射角足夠小，式(9)、(10)可分別簡化為：

其中，V′=1+4/V+3/V2。

2.2分析晶體的位置偏差

DEI中采集圖像需要將分析晶體放置在不同的角度位置，然而由于機械誤差和抖動，實際的位置是圍繞理想位置分布的。這些位置偏差最終會導致提取信息的誤差。假設(shè)分析晶體在上坡位、峰位和下坡位的方向位置偏差分別為UθΔ、PθΔ和DθΔ，同樣運用誤差傳遞理論，可以得到：

可以發(fā)現(xiàn)，上坡位和下坡位的偏差是等權(quán)重的，且峰位的偏差可以忽略。因此，采樣時降低了分析晶體在峰位處的位置精度要求。

信噪比定義為信號變化系數(shù)的倒數(shù)，用信號的平均值與標準方差之比表示。最終，可以得到吸收、折射和散射信息的信噪比為：

其中，Itotal=I0At是入射光子數(shù)。

3 仿真和實驗

實驗在北京同步輻射光源4W1A形貌成像站的微分相襯成像裝置上進行。樣品是圓柱型PMMA和不同層數(shù)的濾紙，PMMA是主要的折射材料，濾紙是主要的散射材料。兩塊Si(111)晶體分別用作單色器和分析器。單色器晶體選擇的能量為15keV，分析晶體的轉(zhuǎn)動精度為0.1 mrad，曝光時間從40ms以5ms為步長增加到100ms，用像素為7.4μm×7.4μm的探測器檢測光子信號。與此同時，運用基于Geant4[18-19]的蒙特卡羅仿真，設(shè)置與實驗一致的條件，模擬粒子與材料的真實作用過程。Geant4是當前模擬粒子輸運過程的主流開源工具，可以很好地描述X射線的粒子性。基于斯涅耳定律引入折射效應(yīng)，使Geant4具有相襯成像能力[20]。構(gòu)建的樣品模型如圖3所示，圓柱是直徑3mm的PMMA，濾紙分為7個不同的高度代表不同的散射能力，緊貼放置于圓柱體之后。

圖3 蒙特卡羅仿真樣品模型Fig.3 Sample construction in Monte Carlo simulations.

按照CFR方法，得到的仿真和實驗結(jié)果如圖4。

可以發(fā)現(xiàn)，隨著紙張厚度增加，吸收和散射逐漸增強，但折射像沒有受到明顯影響。分別計算圖4方框內(nèi)部以及沿虛線的噪聲標準差進行比較，將折射像和散射像的噪聲偏差分別以曝光時間和等效可見度為因變量進行關(guān)聯(lián)并擬合，結(jié)果見圖5。

圖4 實驗提取的吸收(a)、折射(b)和散射(c)方差像以及蒙特卡羅模擬得到的對應(yīng)結(jié)果(d-f)Fig.4Absorption (a), refraction (b) and scattering (c) images extracted from experimental data, while (d-f) are those of Monte Carlo simulations.

圖5 折射角的標準差和散射方差的標準差隨時間和搖擺曲線等效可見度的變化關(guān)系(a)和(b)表示折射角的結(jié)果，(c)和(d)表示散射方差像的結(jié)果Fig.5 Standard deviations σθxand σplotted vs. the exposure time and the equivalent visibility of the RC. (a) and (b) show the results of the refraction-angle images, (c) and (d) show the results of scattering images

根據(jù)羅斯判據(jù)[21]，信噪比至少達到5的時候，才有可能100%的區(qū)分圖像細節(jié)?；诜抡娼Y(jié)果，用每個像素的平均光子數(shù)描述探測器平面的平均強度，計算出圓柱邊界處的信噪比隨平均強度變化的關(guān)系如圖6所示。

圖6 吸收、折射和散射像在邊界處的信噪比隨平均光子數(shù)的變化趨勢Fig.6 SNRs at the boundary of the sample vs. the mean intensity in the detector plane.

模擬結(jié)果顯示，信噪比隨探測器平面的平均強度而提高。相比吸收像，即使探測器平面的光子數(shù)僅為100，折射像和散射像的在樣品邊界處信噪比也是足夠的。體現(xiàn)了折射像和散射像在表達由輕元素組成的物質(zhì)的邊界時具有優(yōu)勢。

4 討論

文中折射定義為相鄰面積元之間的密度不均勻引起的相干散射，其對應(yīng)的物理圖像是單一方向光線入射在樣品上，雖然相鄰面積元的出射光線方向有可能不同，但是每個面積元的出射光線只有一個方向；小角散射定義為面積元內(nèi)部的密度不均勻引起的相干散射，其對應(yīng)的物理圖像是單一方向光線入射在一個面積元上，出射光線有多個方向。由此可知，在衍射增強成像中折射和小角散射的區(qū)分不是絕對的而是相對的，兩者在一定條件下可以相互轉(zhuǎn)化。當成像分辨率降低時，樣品的面積元增大，一部分折射轉(zhuǎn)化為小角散射；而當成像分辨率提高時，樣品面積元減小，一部分小角散射轉(zhuǎn)化為折射。

基于CFR計算出DEI的噪聲特性，與之前用最大似然估計(Maximum Likelihood Estimation, MLE)[10]所得的結(jié)果一致。值得注意的是，這個結(jié)果也類似于用相步進方法[14]或者正反像投影方法[22]在光柵微分相襯成像中的結(jié)論。因為這兩種相襯成像方法的共性在于提取的相襯信息均是折射角。關(guān)于吸收、折射和散射方差這三種信息之間的協(xié)方差可以通過協(xié)方差矩陣來計算，衍射增強CT系統(tǒng)中的噪聲特性將作為下一步的工作進行開展。此外，文中所采用的噪聲模型考慮了光源的光子數(shù)波動和分析晶體的位置偏差。實際上，其他因素諸如暗電流噪聲、晶體缺陷和熱不穩(wěn)定性也在影響結(jié)果，但他們都太小以致可以忽略。

余弦曲線擬合方法還有一個優(yōu)勢：可以用三塊分析晶體，分別定位在角度信號響應(yīng)曲線的上坡腰位、峰位和下坡腰位，同步輻射單色X射線光束一次掃描樣品，可以同時拍攝到上坡腰位像、峰位像和下坡腰位像，從而利用余弦曲線擬合方法進行信息分離。這個設(shè)想為醫(yī)學和生物研究者帶來了新的機遇，有望在降低輻射損傷和提高成像效率方面取得進步。

5 結(jié)語

基于角度信號成像理論，分析了余弦擬合算法的噪聲特性。通過蒙特卡羅模擬和實驗驗證了此算法和噪聲結(jié)果的正確性。隨著通量的增加、等效可見度的提高，噪聲引起的圖像退化將迅速減小。在樣品邊界處，折射角信號和散射方差信號的信噪比高于吸收系數(shù)信號，探測器平面的平均光子數(shù)只要達到100，便可滿足羅斯判據(jù)。本文的噪聲分析工作建立了DEI和GDPCI之間更直觀的聯(lián)系，并且廣泛適用于角度信號成像模型下的各種噪聲測量和分析。

1 Forster E, Goetz K, Zaumseil P. Double crystal diffractometry for the characterization of targets for laser fusion experiments[J]. Crystal Research and Technology, 1980, 15(8): 937-945. DOI: 10.1002/crat.19800150812

2 Chapman D, Thomlinson W, Johnston R, et al. Diffraction enhanced X-ray imaging[J]. Physics in Medicine and Biology, 1997, 42(11): 2015. DOI: 10.1088/0031-9155/ 42/11/001

3 Hu C, Zhang L, Li H, et al. Comparison of refraction information extraction methods in diffraction enhanced imaging[J]. Optics Express, 2008, 16(21): 16704-16710. DOI: 10.1364/OE.16.016704

4 Oltulu O. A unified approach to X-ray absorptionrefraction-extinction contrast with diffraction-enhanced imaging[D]. Illinois Institute of Technology, UMI Number: 3087854, 2003

5 Wernick M N, Wirjadi O, Chapman D, et al. Multiple-image radiography[J]. Physics in Medicine and Biology, 2003, 48(23): 3875-3895. DOI: 10.1088/ 0031-9155/48/23/006

6 Khelashvili G, Brankov J G, Chapman D, et al. A physical model of multiple-image radiography[J]. Physics in Medicine and Biology, 2006, 51(2): 221-236. DOI: 10.1088/0031-9155/51/2/003

7 Rigon L, Besch H J, Arfelli F, et al. A new DEI algorithm capable of investigating sub-pixel structures[J]. Journal of Physics D: Applied Physics, 2003, 36(10A): 107-112. DOI: 10.1088/0022-3727/36/10A/322

8 Rigon L, Arfelli F, Menk R H. Three-image diffraction enhanced imaging algorithm to extract absorption, refraction, and ultrasmall-angle scattering[J]. Applied Physics Letters, 2007, 90(11): 114102. DOI: 10.1063/ 1.2713147

9 趙雪嬌, 張凱, 洪友麗, 等. 衍射增強成像提取多種信息的簡便方法研究[J]. 物理學報, 2013, 62(12): 124202. DOI: 10.7498/aps.62.124202 ZHAO Xuejiao, ZHANG Kai, HONG Youli, et al. A simple method of extracting multiple-information with diffraction enhanced imaging[J]. Acta Physica Sinica, 2013, 62(12): 124202. DOI: 10.7498/aps.62.124202

10 Brankov J G, Iz-Herranz S A, Wernick M N. Noise analysis for diffraction enhanced imaging[A]. 2004 2ndIEEE International Symposium on Biomedical Imaging: Macro to Nano[C]. New York: IEEE, 2004, Vols 1 and 2: 1428-1431

11 Revol V, Kottler C, Kaufmann R, et al. Noise analysis of grating-based X-ray differential phase contrast imaging[J]. Review of Scientific Instruments, 2010, 81(7): 073709. DOI: 10.1063/1.3465334

12 Tang X Y, Yang Y, Tang S J. Characterization of imaging performance in differential phase contrast CT compared with the conventional CT: spectrum of noise equivalent quanta NEQ(k)[J]. Medical Physics, 2012, 39(7): 4467-4482. DOI: 10.1118/1.4730287

13 Diemoz P, Bravin A, Langer M, et al. Analytical and experimental determination of signal-to-noise ratio and figure of merit in three phase-contrast imaging techniques[J]. Optics Express, 2012, 20(25): 27670-90. DOI: 10.1364/OE.20.027670

14 Weber T, Bartl P, Durst J, et al. Measurements and simulations analysing the noise behaviour of grating-based X-ray phase-contrast imaging[J]. Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, 2010, 648: S273-S5. DOI: 10.1016/ j.nima.2010.11.060

15 Rizzi J, Mercere P, Idir M, et al. X-ray phase contrast imaging and noise evaluation using a single phase grating interferometer[J]. Optics Express, 2013, 21(14): 17340-17351. DOI: 10.1364/OE.21.017340

16 Li K, Bevins N, Zambelli J, et al. Fundamental relationship between the noise properties of grating-based differential phase contrast CT and absorption CT: theoretical framework using a cascaded system model and experimental validation[J]. Medical Physics, 2013, 40(2): 021908-021922. DOI: 10.1118/1.4788647

17 Chen G H, Zambelli J, Li K, et al. Scaling law for noise variance and spatial resolution in differential phasecontrast computed tomography[J]. Medical Physics, 2011,38: 584-588. DOI: 10.1118/1.3533718

18 Agostinelli S, Allison J, Amako K A, et al. GEANT4-a simulation toolkit[J]. Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, 2003,506(3): 250-303. DOI: 10.1016/S0168-9002(03)01368-8

19 Allison J, Amako K, Apostolakis J, et al. Geant4 developments and applications[J]. IEEE Transactions on Nuclear Science, 2006,53(1): 270-278. DOI: 10.1109/ TNS.2006.869826

20 Wang Z, Huang Z, Zhang L, et al. Implement X-ray refraction effect in Geant4 for phase contrast imaging[A]. Bo Y. 2009 IEEE Nuclear Science Symposium and Medical Imaging Conference (NSS/MIC 2009)[C]. Piscataway: IEEE, 2009: 2395-2398. DOI: 10.1109/ NSSMIC.2009.5402180

21 Bushberg J T, Boone J M. The essential physics of medical imaging[M]. Philadelphia, Pa: Lippincott Williams & Wilkins, 2012: 86-92. DOI: 10.1148/radiol. 14144045

22 Zhu P, Zhang K, Wang Z, et al. Low-dose, simple, and fast grating-based X-ray phase-contrast imaging[J]. Proceedings of the National Academy of Sciences, 2010,107(31): 13576-13581. DOI: 10.1073/pnas.1003198107

Noise analysis of cosine fitting radiography for diffraction enhanced imaging

BAO Yuan1,2WANG Yan2ZHU Peiping1,2WU Ziyu1,2

1(National Synchrotron Radiation Laboratory,University of Science and Technology of China,Hefei 230029,China)
2(Institute of High Energy Physics,Chinese Academy of Sciences,Beijing 100049,China)

Background:Diffraction enhanced imaging (DEI) is a powerful phase contrast imaging technique.Purpose:This study aims to investigate the noise properties of a new simple cosine fitting radiography for DEI based on the angular signal response imaging prototype.Methods:The rocking curve can be treated as a angular signal response and a simple multi-information retrieval algorithm based on the cosine function fitting. A comprehensive analysis of noise properties in DEI, based on a noise model including the fluctuation of the photon source, and a bias of the analyzer crystal are also given. The validations have been performed with Monte Carlo simulations based on the Geant4 toolkit combined with the refractive process of X-ray and synchrotron radiation experimental data.Results:The cosine fitting radiography (CFR) has been verified by Monte Carlo simulations and experimental data which are in good agreement with each other. The noise penalty is drastically reduced with high photon flux, high visibility and high angular precision.Conclusions:Based on the angular signal response function prototype, we analyze and calculate the noise properties of DEI. In addition, the analytical method can build a strong connection between DEI and grating-based phase contrast imaging (GDPCI) and widely suitable for all kinds of measurement noises in angular signal response imaging prototype. The analysis is useful to understand the noise characteristics of DEI images and improve the signal-noise-ratio of the extracted information for biomedical imaging and material science.

Diffraction enhanced imaging, Noise analysis, Cosine fitting radiography, Monte Carlo simulations

TL99

10.11889/j.0253-3219.2015.hjs.38.090101

基礎(chǔ)研究國家重點項目(No.2012CB825801)、創(chuàng)新研究團體基金項目(No.11321503)、國家自然科學基金(No.11205157、No. 11205189、No.11375225、No.U1332109)資助

鮑園，男，1990年出生，2011年畢業(yè)于安徽大學，現(xiàn)為碩士研究生

朱佩平，E-mail: zhupp@ihep.ac.cn；吳自玉，E-mail: wuzy@ustc.edu.cn

2014-08-28，

2015-03-10

CLCTL99