高 峰，李 嬌，張麗敏，趙會娟

(天津大學(xué)精密儀器與光電子工程學(xué)院，天津 300072)

熒光分子層析(fluorescence molecular tomography，F(xiàn)MT)[1-3]是結(jié)合近紅外擴散熒光層析(fluorescence diffuse optical tomography，F(xiàn)DOT)理論和近紅外熒光探針技術(shù)以實現(xiàn)生物活體內(nèi)特異大分子生化過程的無損三維定量觀測方法，它可提供其他分子成像和現(xiàn)有平面熒光強度成像技術(shù)所不具備的超靈敏度和三維層析能力．FDOT作為 FMT的成像基礎(chǔ)技術(shù)已在連續(xù)光和頻域調(diào)制測量模式下獲得了原理性成功和初步應(yīng)用，而時域FDOT模式研究則處于基本測量技術(shù)和理論體系建立階段．由于時域 FDOT擁有多組分熒光產(chǎn)率和壽命圖像的同時重建能力和高信噪比完整測量信息的獲取能力，從而在小動物模型應(yīng)用中獲得了廣泛的重視．相對于熒光產(chǎn)率，熒光壽命客觀地提供了特異分子的生化微環(huán)境信息，對生物醫(yī)學(xué)基礎(chǔ)研究和臨床診斷均有重要的意義．

本文提出了基于透射解析模型的時域擴散熒光層析圖像重建方法．算法應(yīng)用廣義脈沖譜技術(shù)和歸一化玻恩比逆模型，其中，廣義脈沖譜技術(shù)通過拉普拉斯變換將耦合擴散方程從時域轉(zhuǎn)換到復(fù)頻域，可實現(xiàn)多參數(shù)、多組份熒光圖像的快速重建，而歸一化玻恩比逆模型表示則可免除系統(tǒng)耦合系數(shù)的標定及消除系統(tǒng)響應(yīng)的影響．本文采用外推條件下無限平板透射光學(xué)擴散方程解析解作為逆模型權(quán)重矩陣的計算基礎(chǔ)，由于求解熒光參數(shù)線性方程的過程的欠定性，將基于行操作的代數(shù)重建技術(shù)用于線性求逆過程，其應(yīng)用優(yōu)勢在于每次迭代只需單行運算與存儲，這有助于采用高密度空間剖分實現(xiàn)高分辨率的圖像重建．在算法驗證研究中，首先進行基于有限元正向模型數(shù)據(jù)的數(shù)值模擬研究，對算法的可行性及空間分辨率和魯棒性等關(guān)鍵性能方面做初步的驗證與評估；然后采用所提算法和多通道時間相關(guān)單光子計數(shù)(time-correlated single photon counting，TCSPC)測量系統(tǒng)對自制仿體進行圖像重建實驗，進一步驗證了重建算法的有效性和準確性．

1 算法原理

采用擴散方程作為光在組織體中的傳輸模型，熒光檢測量的數(shù)學(xué)表達式可由激發(fā)光和出射熒光之間耦合擴散方程得到[4]．另外采用了在解決以波動和擴散為主導(dǎo)現(xiàn)象的逆向問題中，已被證明的一種有效方法——廣義脈沖譜技術(shù)(generalized pulse spectrum technique，GPST)，通過拉氏變換將時域擴散方程轉(zhuǎn)化到復(fù)頻域，從而使原四維全時空問題降為與時間相關(guān)的三維問題[5-6]．

時域擴散熒光的逆向問題即時域 FDOT圖像重建，其目的是在背景光學(xué)參數(shù)分布已知的情況下，實現(xiàn)對熒光產(chǎn)率和熒光壽命空間分布的同時重建．即基于已得光子傳輸模型的正向算子，來求解所需的熒光參數(shù)分布，從而實現(xiàn)擴散熒光層析成像的圖像重建．具體過程如下：源在sr處接收在dr處的熒光光子密度可看作是所有熒光體元 dV在整個體積上的積分，從而熒光密度可表示為[7]

式中：Φm( rd, rs, q)為源位于 rs處在 rd處探測到的熒光光子密度的拉普拉斯變換值；Φx( r, rs, q)為源位于 rs處接收在r處的激發(fā)光光子密度的拉普拉斯變換值；Gm(rd, r, q)為源位于r處接收在 rd處的擴散方程的格林函數(shù)的拉普拉斯變換值[8]．

在檢測量處理方面，為了重建熒光參數(shù)，采用了歸一化玻恩比[7]，即用在檢測器處探測到的熒光光流量除以相應(yīng)的激發(fā)光的光流量，這樣可以有效地消除不同的源和探測器對計算引起的不便．由式(1)推導(dǎo)及檢測量的處理后得

在無限平板透射模型中， Φx( ri)為外推邊界條件下的拉普拉斯變換后的擴散方程的解析解[8]，即

式中：r±,m為根據(jù)外推邊界近似中對稱的正負源位置；zb為外推邊界的坐標；z0= 1 /μs′為單一的近似各向同性源入射光線位置；為有效衰減系數(shù)．

由式(3)中的離散過程導(dǎo)致了所需求解的體元上的熒光參數(shù)個數(shù)遠遠多于測量數(shù)據(jù)個數(shù)，使得對體元熒光參數(shù)的求解過程變?yōu)榍范▎栴}，從而意味著測量數(shù)據(jù)上的微小變化可能引起重建圖像的完全變異．且線性方程組的解易受到噪聲干擾，因此很難直接用矩陣求逆的方法得到．在此情形下，只能按照某種近似原則將原問題的直接求解轉(zhuǎn)化為求原問題的穩(wěn)定的、合理的近似解，該過程稱為正規(guī)化過程．本文采用了代數(shù)重建技術(shù)(algebraic reconstruction technique，ART)[9]，它是根據(jù)求解線性方程組思想，以解決離散化數(shù)字圖像重建為目的，通過對矩陣進行逐行計算，從而有效快速地對方程進行求解[10]．為了同時求出熒光產(chǎn)率af()ημr及熒光壽命()τr，本文還引入了一對實變換因子對：．其中分別為背景在激發(fā)光和熒光波長下的吸收系數(shù)為背景的熒光壽命．由式(1)，可以得出其熒光參數(shù)為

2 模擬研究的重建結(jié)果討論及分析

本文在模擬研究中應(yīng)用了同課題組的有限元方法對非均勻目標體產(chǎn)生正向問題的數(shù)據(jù)，再將所得數(shù)據(jù)代入設(shè)計的算法中進行圖像重建．圖 1(a)所示為設(shè)計的非均勻目標體，通過在正向數(shù)據(jù)中混入白噪聲及改變目標體的中心間距(center-to-center separation，CCS)來對算法進行噪聲魯棒性和空間分辨率的研究．首先討論噪聲對圖像重建的影響，結(jié)果表明圖像的抗噪性受到q取值的影響：當時，如圖 1(b)所示，重建圖像的熒光壽命可以承受 50,dB以上的噪聲，抗噪聲的能力較差；而產(chǎn)率幾乎不受噪聲影響．當時，如圖1(c)所示，重建的熒光壽命可以承受 35,dB以上的噪聲，有較好的抗噪能力；而產(chǎn)率相對的抗噪聲性能變差．且熒光壽命和產(chǎn)率的量化值都變大，從而影響了重建算法的分辨率．由以上的重建結(jié)果可以看出，q值增大可以提高結(jié)果中熒光壽命的抗噪性，但也影響了成像的質(zhì)量，犧牲了圖像的分辨率．

圖1 模型的模擬結(jié)果Fig.1 Simulated results of model

對于算法空間分辨能力的討論結(jié)果，如圖 1(d)和(e)所示，隨著 CCS的減小，在熒光產(chǎn)率和熒光壽命的重建中，兩個目標體之間可分度均逐漸下降．當CCS大于 10,mm 時，熒光產(chǎn)率可以分辨；而當 CCS大于14,mm時，熒光壽命可以分辨．

由上述結(jié)果可以看出不同熒光壽命和熒光產(chǎn)率區(qū)域可以明顯區(qū)分，重建圖像質(zhì)量較好．同時值得注意的是，熒光產(chǎn)率的重建深度及量化精度要比熒光壽命的好些．對重建算法的空間分辨率進行測試，結(jié)果表明熒光產(chǎn)率比壽命好，壽命在深度信息上有偏差．

3 實驗及分析

3.1 實驗材料及設(shè)備

3.1.1 測量系統(tǒng)

實驗研究基于本實驗室所發(fā)展的多通道TCSPC皮秒時間分辨測量系統(tǒng)，采用透射檢測模式．測量過程為：皮秒脈沖激光器發(fā)出波長為660,nm的激發(fā)光，經(jīng)直徑為 62.5,μm、數(shù)值孔徑為 0.22的光纖導(dǎo)出，之后1∶16光開關(guān)將激發(fā)光源依次導(dǎo)入16個源位置．對于每次源入射，4個 4∶1光開關(guān)分4次、每次選擇 4個檢測點做并行TCSPC 檢測，其中4個并行檢測點被置于相鄰區(qū)域，探測光纖的直徑為 500,μm、數(shù)值孔徑為 0.37．隨后切換源位置重復(fù)上述測量，直至 16個源位置“掃描”完畢．其中，在測量熒光時，接收光纖需經(jīng)過帶通濾波器將激發(fā)光濾除．實驗測量系統(tǒng)示意如圖2所示．

圖2 TCSPC測量系統(tǒng)示意Fig.2 Measurement system of TCSPC

3.1.2 實驗仿體選用固體模型作為背景，其主要組成成分是聚甲醛(Polyformaldehyde)．如圖3所示，仿體為100,mm×25,mm×70,mm的長方體，吸收系數(shù)μa(r)=0.003,8 mm-1，約化散射系數(shù)(實驗中所有光學(xué)參數(shù)均基于反射模式的光學(xué)參數(shù)時域測量方法測得).根據(jù)經(jīng)驗選擇其熒光差率 ημaf(r)＝0.000,01 mm-1，熒光壽命 τ(r)＝10,000,ps．在本實驗中，設(shè)計了2個固態(tài)仿體：仿體 1為背景中含有一個目標體，仿體 2背景中則包含了兩個相同幾何尺寸的目標體且其中心距為 20,mm．目標體是直徑為 5,mm、高為60,mm的圓柱形孔，由1% Intralipid 溶液和Cy5.5熒光染料(激發(fā)和發(fā)射峰值分別在 670,nm 和 710,nm)的混合溶液填充．其吸收系數(shù)μa(r)=0.001,7,mm-1，約化散射系數(shù)目標體的熒光產(chǎn)率和壽命為需重建的參數(shù)．在重建過程中只需在視場范圍(field of view，F(xiàn)OV)進行重建．由于正向模型計算采用了外推邊界條件，因此重建區(qū)域為 64,mm×29,mm×60,mm的長方體．光源的位置在外推邊界條件下，分布于 Z＝2,mm 的 XY平面，視場范圍為21,mm×21,mm，光源間隔為 7,mm×7,mm，共有 16個光源．16個探測器的位置分布于Z＝27,mm的XY平面，與光源位置相對．目標體中心坐標(X＝32,mm，Z＝15,mm)．離散化后體積元為V＝L×W×H＝1,mm×1,mm×1,mm的立方體．整個計算模型可離散成111,360個立方體單元．圖 3(a)為仿體、計算模型及坐標系示意．在 ART計算過程中迭代次數(shù)為15，松弛因子選為0.5．

圖3 計算和實驗?zāi)Ｐ褪疽釬ig.3 Calculating and experimental model

3.2 實驗測量及結(jié)果分析

在實驗測量時，TCSPC各部分的設(shè)定如下：TAC范圍為70.01×10-9；CFD最低限定為-15.69；PMT增益為 80%，共有 4,096個采樣點，采樣時間間隔為17.1,ps．在整個測量過程中，整個系統(tǒng)在黑暗狀態(tài)下，以防止外部光源影響．激發(fā)光和熒光測量過程中積分時間均固定在 10,s，改變激發(fā)光的光強．圖 4(a)顯示了第 1個光源激發(fā)時，16個探測通道測得的激發(fā)光(EXTdata)和熒光(EMSdata)的歸一化時間擴展曲線．可以看出，熒光信號較激發(fā)光有一定的時間延遲，且有一定程度的展寬．

圖4(b)和(c)分別為仿體 1和仿體 2的熒光產(chǎn)率和壽命空間分布的重建結(jié)果．由重建圖像的XY平面的剖面圖看出算法對熒光產(chǎn)率圖像進行了很好的重建，重建目標的形狀、位置、大小均與實際目標相符．重建目標的熒光壽命形狀和位置與實際情況較符合，但體積要比實際目標大，與模擬實驗結(jié)果也是相符的．特別在仿體 2的重建圖像中可以很好分辨兩個目標體，從而體現(xiàn)了算法良好的分辨力．兩個仿體的成功重建效果說明了算法具有可靠性且適用于不同目標體的成像．但在重建圖像中目標邊界處出現(xiàn)了偽逆，尤其熒光壽命圖像較明顯，我們認為此現(xiàn)象是由于目標與背景交接處的光學(xué)特性(主要是熒光特性)變換較大，使得熒光擴散光方程產(chǎn)生了誤差．因此，由重建圖像及實驗討論可知，實驗結(jié)果很好地驗證了本文提出的透射解析模型時域擴散熒光層析重建算法可行性，且與模擬結(jié)果相似．

圖4 測量和重建結(jié)果Fig.4 Measurement and reconstruction results

4 結(jié) 語

本文提出了基于透射解析模型時域擴散熒光層析重建方法．應(yīng)用有限元方法提供正向問題的數(shù)據(jù)，在數(shù)值模擬的基礎(chǔ)上評價了算法的空間分辨率及魯棒性，驗證了算法具有強度測量成像模式無可比擬的空間分辨率和多參數(shù)成像功能．另外采用此重建方法對自行制作非均勻仿體進行圖像重建，結(jié)果中對于仿體的熒光產(chǎn)率和熒光壽命的良好重建，驗證了本重建算法的可行性和實用前景．

［1］ Ntziachristos V，Tung C H，Bremer C，et al. Fluorescence molecular tomography resolves protease activity in vivo [J]. Nature Medicine，2002，8(7)：757-760.

［2］ 徐可欣，高 峰，趙會娟．生物醫(yī)學(xué)光子學(xué)[M].北京：科學(xué)出版社，2007.Xu Kexin，Gao Feng，Zhao Huijuan. Biomedical Photonics [M]. Beijing：Science Press，2007(in Chinese).

［3］ Cherry S R. In vivo molecular and genomic imaging：New challenges for imaging physics[J]. Physics in Medicine and Biology，2004，49：13-48.

［4］ O’Leary M A，Boas D A，Li X D. Fluorescence lifetime imaging in turbid media[J]. Optics Letters，1996，21(2)：158-160.

［5］ Gao Feng，Zhao Huijuan，Yamada Y，et al. Timeresolved diffuse optical tomography using a modified generalized pulse spectrum technique [J]. IEICE Transactions on Information and Systems，2002，E85-D (1)：133-142.

［6］ Gao Feng，Zhao Huijuan，Yamada Y．A linear，featured-data scheme for image reconstruction in timedomain fluorescence molecular tomography[J]. Optics Express，2006，14(16)：7109-7124.

［7］ Ntziachristos V，Weissleder R. Experimental threedimensional fluorescence reconstruction of diffuse media by use of a normalized Born approximation[J]. Optics Letters，2001，26(12)：893-895.

［8］ Andrea C D，Spinelli L，Comelli D，et al. Localization and quantification of fluorescent inclusions embedded in a turbid medium[J]. Physics in Medicine and Biology，2005，50：2313-2327.

［9］ 赫爾曼 G T. 由投影重建圖像：CT的理論基礎(chǔ) [M].北京：科學(xué)出版社，1985：25-38.Herman G T. Image Reconstruction from Projections：The Fundamentals of Computerized Tomography [M].Beijing：Science Press，1985：25-38(in Chinese).

［10］ Kienle A，Patterson M S. Improved solutions of the steady-state and the time-resolved diffusion equations for reflectance from a semi-infinite turbid medium [J].Journal of the Optical Society of America，1997，A14：246-254.

［11］ Kumar A T N，Skoch J，Bacskai J．Fluorescent-lifetimebased tomography for turbid media[J]. Optics Letters，2005，30(24)：3347-3349.