時(shí)域光學(xué)拓?fù)涑上穹椒ǖ哪M和實(shí)驗(yàn)研究

高 峰，段臨晶，李 嬌，張 偉，易 茜，武林會(huì)，王 欣

(1.天津大學(xué)精密儀器與光電子工程學(xué)院，天津 300072；2.天津市生物醫(yī)學(xué)檢測技術(shù)與儀器重點(diǎn)實(shí)驗(yàn)室，天津 300072)

時(shí)域光學(xué)拓?fù)涑上穹椒ǖ哪M和實(shí)驗(yàn)研究

高 峰1,2，段臨晶1，李 嬌1,2，張 偉1，易 茜1，武林會(huì)1，王 欣1

(1.天津大學(xué)精密儀器與光電子工程學(xué)院，天津 300072；2.天津市生物醫(yī)學(xué)檢測技術(shù)與儀器重點(diǎn)實(shí)驗(yàn)室，天津 300072)

針對現(xiàn)有的以修正的朗伯-比爾定律(modified Lambert-Beer law，MLBL)為理論模型的光學(xué)拓?fù)浞椒ù嬖诹炕炔粶?zhǔn)確的問題，提出了一種基于最小二乘擬合算法的時(shí)域光學(xué)拓?fù)涑上穹椒ǎ摲椒ɡ脮r(shí)域擴(kuò)散方程在半無限空間條件下的解析解和測量曲線實(shí)現(xiàn)最佳匹配以計(jì)算相鄰的源點(diǎn)與探測點(diǎn)之間的待測組織吸收系數(shù)變化值．對該時(shí)域光學(xué)拓?fù)涑上穹椒ㄟM(jìn)行了一系列的模擬驗(yàn)證，并以基于時(shí)間相關(guān)單光子計(jì)數(shù)(time-correlated single photon counting，TCSPC)技術(shù)的多通道時(shí)間分辨測量系統(tǒng)為平臺(tái)進(jìn)行了實(shí)驗(yàn)驗(yàn)證．結(jié)果表明：所提出的最小二乘擬合算法總體上優(yōu)于傳統(tǒng)的基于經(jīng)驗(yàn)值光程的MLBL方法．

擴(kuò)散方程；飛行時(shí)間；時(shí)間相關(guān)單光子計(jì)數(shù)；光學(xué)拓?fù)洌荒X功能成像

近紅外光譜成像技術(shù)(near-infrared spectroscopy，NIRS)以其安全無損、非電離輻射等優(yōu)點(diǎn)而在醫(yī)學(xué)成像領(lǐng)域得到廣泛關(guān)注．選擇近紅外“治療窗口”范圍內(nèi)的波長，不僅可以實(shí)現(xiàn)對活體組織的無損傷、實(shí)時(shí)地探測和成像，還可以直接同時(shí)提供腦部氧合血紅蛋白和還原血紅蛋白的濃度變化信息[1]．

由于近紅外光在腦部組織穿透深度的限制，一般采用拓?fù)涞姆绞窖芯看竽X皮層應(yīng)激反應(yīng)過程．光學(xué)拓?fù)涑上窦词抢梅瓷錅y量給出待測組織體表面下淺層的光學(xué)參數(shù)二維變化．與其他成像方式相比，光學(xué)拓?fù)涑上窬哂徐`活可變的源-探測陣列排布；由于只選擇測量與源相鄰的探測點(diǎn)處的反射光流，并將其與中間采樣點(diǎn)相關(guān)聯(lián)，因而具有較高的信噪比．

近紅外光學(xué)拓?fù)涑上窨刹捎?種測量模式：連續(xù)光、頻域和時(shí)域[2]．時(shí)域模式在信息完整性、數(shù)據(jù)靈活性、系統(tǒng)穩(wěn)定性以及隨之體現(xiàn)的成像質(zhì)量穩(wěn)健性、多參數(shù)重建等諸多關(guān)鍵性能上具有其他測量模式不可比擬的綜合優(yōu)勢[3]．

針對現(xiàn)有的光學(xué)拓?fù)涑上穹椒ɡ碚撃Ｐ秃唵巍⒘炕`差大等問題，筆者提出了一種利用時(shí)域擴(kuò)散方程在半無限空間條件下的解析解和測量曲線實(shí)現(xiàn)最佳匹配以計(jì)算相鄰的源點(diǎn)與探測點(diǎn)之間的待測組織吸收系數(shù)變化值的時(shí)域光學(xué)拓?fù)涑上穹椒ǎ疚膶υ摃r(shí)域光學(xué)拓?fù)涑上穹椒ㄟM(jìn)行了一系列的模擬實(shí)驗(yàn)驗(yàn)證，并以基于時(shí)間相關(guān)單光子計(jì)數(shù)(time-correlated single photon counting，TCSPC)技術(shù)[4]的多通道時(shí)間分辨測量系統(tǒng)為平臺(tái)進(jìn)行了實(shí)驗(yàn)驗(yàn)證．實(shí)驗(yàn)結(jié)果表明，相較于基于經(jīng)驗(yàn)值光程的修正的 Lambert-Beer定律(modified Lambert-Beer law，MLBL)方法的量化度過估計(jì)和成像區(qū)域不均勻現(xiàn)象，時(shí)域測量方法對異質(zhì)體的重建值和重建形狀更為理想．該時(shí)域光學(xué)拓?fù)涑上穹椒ú粌H可以靈活設(shè)計(jì)源-探測器排布陣列，亦可采取多波長測量模式，獲取生物組織更多組分信息，因而有望成為現(xiàn)有腦功能成像方法的有益補(bǔ)充．

1 近紅外時(shí)域光學(xué)拓?fù)涑上穹椒?/h2>

1.1 基于時(shí)域擴(kuò)散方程解析解的最小二乘擬合算法

對于光在低吸收、高散射組織體的傳播行為，擴(kuò)散方程應(yīng)用得較為廣泛．時(shí)域形式的擴(kuò)散方程表示如下：

反射測量方式下，源和探測器距離較小，相較于大腦輪廓，每對源-探測器下的組織體表面曲率很大，故本文將待測組織體近似為平面半無限空間．在平面半無限空間媒質(zhì)外推邊界條件下，組織體表面距源點(diǎn)ρ處的探測點(diǎn)反射光流量為

在實(shí)際臨床應(yīng)用中，研究者更關(guān)心光學(xué)參數(shù)的相對變化量．在此，本文定義待測組織在光學(xué)參數(shù)變化前的狀態(tài)為“rest”，光學(xué)參數(shù)變化后的狀態(tài)為“task”.待測組織從 rest狀態(tài)變化到 task狀態(tài)期間，約化散射系數(shù)改變量遠(yuǎn)小于吸收系數(shù)改變量可忽略不計(jì)[5]．結(jié)合式(2)，有如下關(guān)系式：

實(shí)際測量中，一般要考慮儀器響應(yīng)函數(shù)(instrument response function，IRF)的影響，俗稱“空測曲線”．對式(3)兩邊做卷積，有

綜上，需要應(yīng)用兩次非線性最小二乘擬合算法．對式(4)左邊應(yīng)用一次最小二乘擬合算法可以同時(shí)重建待測組織rest狀態(tài)的吸收系數(shù)μa0和散射系數(shù)對式(4)右邊應(yīng)用一次最小二乘擬合算法即可獲得Δμa．

在后面的驗(yàn)證過程中，基于時(shí)域擴(kuò)散方程解析解的最小二乘法用LSM表示．

1.2 基于計(jì)算值光程的修正的Lambert-Beer方法

對于光學(xué)厚層媒質(zhì)，由于多次散射的作用增加了光子在組織體中傳輸?shù)穆窂介L度，因此，被接受的光信號在介質(zhì)內(nèi)部的平均傳播距離要大于入射點(diǎn)與接收點(diǎn)之間的幾何距離ρ．光子在組織中的運(yùn)行路程定義為光程，用L表示[6]．MLBL定理可表示為

式中：A為光密度；I0為入射到組織表面的平均光強(qiáng)；I1為探測器在組織表面探測到的平均光強(qiáng)；ε為待測組分的摩爾吸光系數(shù)；C為待測組分的濃度；B為路徑修正因子；G為與μs相關(guān)、與μa無關(guān)的強(qiáng)散射衰減．B和G可認(rèn)為是常數(shù)．

當(dāng)一定光強(qiáng)入射到組織體，設(shè)在 rest、task狀態(tài)下，出射光強(qiáng)分別為Irest和Itask．由MLBL可知，待測組織task狀態(tài)相對于rest狀態(tài)的光密度變化為

由式(6)可知，時(shí)域和連續(xù)波模式的差別在于對光程的確定方法不同．時(shí)域模式可以利用光子在組織體中的平均飛行時(shí)間〈t〉來確定如式(7)所示．連續(xù)波測量模式一般取經(jīng)驗(yàn)值來近似平均光程，如式(8)所示．

在后面的驗(yàn)證過程中，基于計(jì)算值光程的MLBL計(jì)算結(jié)果用飛行時(shí)間的符號＜t＞表示；基于經(jīng)驗(yàn)值光程的 MLBL計(jì)算結(jié)果用經(jīng)驗(yàn)值(empirical value，EV)表示．

2 近紅外光學(xué)拓?fù)涑上穹椒M驗(yàn)證

本文利用蒙特卡洛(Monte-Carlo，MC)模擬產(chǎn)生正向數(shù)據(jù)[2]．MC模擬的數(shù)值模型為80,mm×80 mm×30,mm 的平板．在圖1(a)所示的坐標(biāo)系下，(-7.5,mm，-3.75,mm)位置處有一深度為20,mm的圓柱形凹槽，直徑D分別為 6,mm、7.5,mm、10,mm、14,mm．采用21源-4探測器配置，其中D1、D2、D3、D4為探測點(diǎn)位置，由此形成32個(gè)源-探測器對，對應(yīng)32個(gè)直接測量采樣點(diǎn)和17個(gè)間接采樣點(diǎn)，如圖1(b)所示．源-探測器最小距離ρ=7.5,mm，成像區(qū)域范圍為 30,mm×30,mm．?dāng)?shù)值模型的背景光學(xué)參數(shù)：吸收系數(shù)約化散射系數(shù)；異質(zhì)體光學(xué)參數(shù)：吸收系數(shù)，約化散射系數(shù)和背景保持一致．

利用前面提及的 3種方法處理模擬的 MC數(shù)據(jù)，重建結(jié)果如圖2和圖3所示．

圖1 5×5源-探測器陣列Fig.1 5×5,source-detector array

圖2 5×5源-探測器陣列模擬重建結(jié)果(D=6 mm，7.5 mm)Fig.2 Reconstructed images under a 5×5 source-detector array(D=6 mm，7.5 mm)

圖3 5×5源-探測器陣列模擬重建結(jié)果(D=10 mm，14 mm)Fig.3 Reconstructed images under a 5×5 source-detector array(D=10 mm，14 mm)

隨著異質(zhì)體直徑 D的增大，每種方法的量化度都會(huì)相應(yīng)提高．在所有情況下，按照量化度從大到小排列，依次是基于經(jīng)驗(yàn)值光程的MLBL方法、基于計(jì)算值光程的 MLBL方法、基于時(shí)域擴(kuò)散方程解析解的最小二乘擬合算法．但當(dāng)D大于源和探測器距離ρ時(shí)，基于經(jīng)驗(yàn)值光程的 MLBL方法在量化度方面會(huì)出現(xiàn)過估計(jì)問題，出現(xiàn)偏差．而此時(shí)，時(shí)域測量方法的計(jì)算結(jié)果更加理想．

對于目標(biāo)體的重建結(jié)果形狀方面，基于經(jīng)驗(yàn)值光程的 MLBL方法與實(shí)際偏差最大，經(jīng)常會(huì)出現(xiàn)異質(zhì)體中間值估計(jì)過高而邊緣值估計(jì)過低的現(xiàn)象．相對而言，時(shí)域測量方法的重建結(jié)果中此現(xiàn)象有所改善，其中基于時(shí)域擴(kuò)散方程解析解的最小二乘擬合算法對異質(zhì)體的形狀重建最符合實(shí)際．因此，當(dāng)D超過源和探測器之間的距離ρ時(shí)，時(shí)域測量方法對異質(zhì)體的重建值和重建形狀更為理想．

3 近紅外光學(xué)拓?fù)浞椒▽?shí)驗(yàn)驗(yàn)證

3.1 多通道時(shí)間分辨測量系統(tǒng)

為了進(jìn)一步驗(yàn)證提及的兩種時(shí)域算法的可行性和有效性，本文利用實(shí)驗(yàn)室現(xiàn)有的基于 TCSPC的多通道時(shí)間分辨測量系統(tǒng)[7]進(jìn)行了一系列的仿體驗(yàn)證實(shí)驗(yàn)．系統(tǒng)原理如圖 4所示，包括光源部分、光纖傳輸部分、探測部分和計(jì)算機(jī)．

(1)光源部分．為了更好地反映組織體中血紅蛋白的濃度變化并有效區(qū)分氧合血紅蛋白和還原血紅蛋白，本系統(tǒng)采用兩套發(fā)射波長位于近紅外波段分別為830,nm和780,nm的皮秒半導(dǎo)體激光器系統(tǒng)(PDL 828Sepia Ⅱ，德國 PicoQuant Inc)作為發(fā)射光源．每套激光器系統(tǒng)由激光頭和激光器控制模塊組成．兩套系統(tǒng)激光頭(型號分別為 LHD-P-830和 LHD-P-780，德國 PicoQuant Inc)的發(fā)射波長分別為 830,nm和780,nm，它們的激光器控制模塊相同．激光器的重復(fù)頻率為 80,MHz，功率可調(diào)，峰值功率為 300～400 μW．本系統(tǒng)采用了波分復(fù)用器(WDM，加拿大 OZ Optics Ltd)將兩路激光信號耦合成一路激光信號．1∶16光開關(guān)用于切換源通道，實(shí)現(xiàn)了發(fā)射光源在 16個(gè)不同的入射源點(diǎn)間進(jìn)行切換．光開關(guān)通過FC(ferrule connection)接口與計(jì)算機(jī)進(jìn)行連接和通信．用戶通過計(jì)算機(jī)VB界面控制源通路的選擇．

(2)光纖傳輸部分．包括16孔洞光纖架和16根傳輸光纖．光纖架由鋼質(zhì)材料制成，其表面鍍有一層黑色的氧化膜以防止對光線的反射．光纖架用于固定源、探測光纖以實(shí)現(xiàn)光源從不同源點(diǎn)入射，探測器從不同的探測點(diǎn)接收反射光．為了減少位于組織體上探頭的數(shù)目以及精確地獲得反射光信息，傳輸光纖采用同軸光纖結(jié)構(gòu)，即源光纖位于中心，探測光纖位于外環(huán)．系統(tǒng)包括 16根源光纖，其中 4根同時(shí)用作探測光纖．入射光纖芯徑為 62.5 μm/0.22，探測光纖芯徑為500 μm/0.37．

(3)探測部分．光子經(jīng)過組織體的吸收和散射作用，有的被組織體吸收，有的溢出組織體表面．溢出組織體表面的光子被4根探測光纖所接收，經(jīng)過準(zhǔn)直器(F230FC-B，Thorlabs Inc)和調(diào)節(jié)光強(qiáng)的濾光輪(FW102B，Thorlabs Inc)注入 4個(gè)光電倍增管(photomultiplier tube，PMT)，PMT輸出的光子脈沖通過TCSPC(SPC-134，德國Becker& Hickl Inc)進(jìn)行計(jì)數(shù)，得到出射光子的時(shí)間擴(kuò)展曲線．需要注意的是，由于 PMT為微弱光測量儀器，環(huán)境日光強(qiáng)度較大會(huì)損壞PMT，故本實(shí)驗(yàn)操作環(huán)境必須為暗室．

(4)計(jì)算機(jī)．集成控制各模塊工作，進(jìn)行數(shù)據(jù)處理和圖像重建．

圖4 雙波長時(shí)域光學(xué)拓?fù)涑上裣到y(tǒng)原理Fig.4 Schematic of the dual-wavelength，time domain optical topography system

3.2 仿體驗(yàn)證

本文采用固-液混合仿體，固態(tài)模擬均勻背景，液態(tài)模擬異質(zhì)體．仿體是材質(zhì)為聚甲醛、高H'=50 mm、直徑D'=110,mm的圓柱．在圓柱體正中心，有一個(gè)深20,mm，直徑D分別為10,mm、20,mm的圓柱形凹槽．凹槽的直徑和位置代表異質(zhì)體的大小和位置，如圖 5所示．在凹槽中加入按一定比例配制的Intralipid和印度墨水的混合液，其中 Intralipid用來模擬散射介質(zhì)，印度墨水用來模擬吸收介質(zhì)[8]．相應(yīng)實(shí)驗(yàn)重建結(jié)果如圖6～圖11所示．

實(shí)驗(yàn)結(jié)果表明，3種方法的重建結(jié)果量化值與前面提到過的模擬結(jié)果相似．即按照從大到小排列依次是基于經(jīng)驗(yàn)值光程的 MLBL方法、基于計(jì)算值光程的 MLBL方法、基于時(shí)域擴(kuò)散方程解析解的最小二乘擬合算法．但當(dāng)異質(zhì)體的直徑D=ρ=10,mm的時(shí)候，基于經(jīng)驗(yàn)值光程的 MLBL方法就已經(jīng)開始出現(xiàn)過估計(jì)問題，即在量化度方面會(huì)出現(xiàn)偏差．分析其原因，可能是由于實(shí)驗(yàn)測量過程中的基底噪聲與其他非線性噪聲的影響造成的．而重建結(jié)果的形狀情況亦與模擬結(jié)果相似，即基于經(jīng)驗(yàn)值光程的 MLBL方法亦會(huì)出現(xiàn)異質(zhì)體中間值估計(jì)過高而邊緣值估計(jì)過低的現(xiàn)象．綜上所述，實(shí)驗(yàn)結(jié)果均與模擬結(jié)果相符合，證明了當(dāng)異質(zhì)體尺寸 D接近或超過源和探測器之間的距離ρ 時(shí)，時(shí)域測量方法對異質(zhì)體的重建值和重建形狀更為理想．

圖5 實(shí)驗(yàn)示意Fig.5 Experimental setup

圖6 仿體實(shí)驗(yàn)結(jié)果(D＝10,mm，)Fig.6 Phantom experimental results(D＝10,mm，)

圖7 仿體實(shí)驗(yàn)結(jié)果(D＝10,mm，)Fig.7 Phantom experimental results(D＝10,mm，)

圖8 仿體實(shí)驗(yàn)結(jié)果(D＝10,mm，)Fig.8 Phantom experimental results(D＝10,mm，)

圖9 仿體實(shí)驗(yàn)結(jié)果(D＝20,mm，)Fig.9 Phantom experimental results(D＝20,mm，)

圖10 仿體實(shí)驗(yàn)結(jié)果(D＝20,mm，)Fig.10 Phantom experimental results(D＝20,mm，)

圖11 仿體實(shí)驗(yàn)結(jié)果(D＝20,mm，)Fig.11Phantom experimental results(D＝20,mm，)

4 結(jié) 語

本文提出了一種在時(shí)域測量數(shù)據(jù)和時(shí)域擴(kuò)散方程在平面半無限空間條件下的解析解之間尋求最佳匹配的時(shí)域光學(xué)拓?fù)涑上穹椒?，并進(jìn)行了一系列的模擬和仿體實(shí)驗(yàn)，從而得到待測組織光學(xué)參數(shù)變化的2D圖像．仿體實(shí)驗(yàn)在實(shí)驗(yàn)室已搭建的基于 TCSPC的多通道時(shí)域測量系統(tǒng)平臺(tái)上進(jìn)行．模擬和實(shí)驗(yàn)結(jié)果表明，該方法在量化度上優(yōu)于傳統(tǒng)的基于經(jīng)驗(yàn)值光程的MLBL方法．

［1］ Zhao Huijuan，Gao Feng. Time-resolved diffuse optical tomographic imaging for the provision of both anatomical and functional information about bio1ogical tissue [J]. Applied Optics，2005，44(10)：1905-1916.

［2］ 徐可欣，高 峰，趙會(huì)娟. 生物醫(yī)學(xué)光子學(xué)[M]. 北京：科學(xué)出版社，2011.

Xu Kexin，Gao Feng，Zhao Huijuan. Biomedical Photonics[M]. Beijing：Science Press，2011(in Chinese).

［3］ Correia T，Lloyd-Fox S，Everdell N，et al. Threedimensional optical topography of brain activity in in fants watching videos of human movement [J]. Physics in Medicine and Biology，2012，57(5)：1135-1146.

［4］ Becker W，屈軍樂. 高級時(shí)間相關(guān)單光子計(jì)數(shù)技術(shù)[M]. 北京：科學(xué)出版社，2009.

Becker W，Qu Junle. The Technique About Advanced Time-Correlated Single Photon Counting [M]. Beijing：Science Press，2009(in Chinese).

［5］ Taroni P，Bassi A，Spinelli L，et al. Time-resolved diffuse optical spectroscopy：A differential absorption approach [J]. Appl Spectrosc，2010，64(11)：1220-1226.

［6］ Liebert Adam，Kacprzak Michal，Maniewski Roman. Time-resolved reflectometry and spectroscopy for assessment of brain perfusion and oxygenation[J]. Biocybernetics and Biomedical Engineering，2007，27(2)：237-266.

［7］ Contini D，Torricelli A，Pifferi A，et al. Multichannel time-resolved tissue oximeter for functional imaging of the brain[J]. IEEE Transactions on Instrumentation and Measurement，2006，55(1)：85-90.

［8］ Gao Feng，Zhao Huijuan，Tanikawa Yukari，et al. Optical tomographic mapping of cerebral hemodynamic by means of time-domain detection：Methodology and phantom validation[J]. Physics in Medicine and Biology，2004，49(6)：1055-1078.

（責(zé)任編輯：趙艷靜）

Simulation and Experimental Investigation on Time Domain Optical Topography

Gao Feng1,2，Duan Linjing1，Li Jiao1,2，Zhang Wei1，Yi Xi1，Wu Linhui1，Wang Xin1
(1. School of Precision Instrument and Opto-Electronics Engineering，Tianjin University，Tianjin 300072，China；2. Tianjin Key Laboratory of Biomedical Detecting Techniques and Instruments，Tianjin 300072，China)

To cope with the low quantification in the established optical topography that originates from the excessively simplified computation model based on the modified Lambert-Beer’s Law(MLBL)，we proposed a leastsquares fitting scheme for time domain optical topography that seeks data matching between the time-resolved measurement and model prediction calculated by analytically solving the time-domain diffusion equation in semi-infinite geometry. A series of experiments were also conducted，including simulation and phantom verification．The phantom verification was performed with multi-channel time-resolved measurement system based on TCSPC．Our simulative and phantom experiments demonstrate that the proposed curve-fitting method is on the whole，superior to the conventional MLBL-based one in quantitative performance.

diffusion equation；flight time；time-correlated single photon counting；optical topology；functional cerebral imaging

Q63

：A

：0493-2137(2014)09-0829-07

10.11784/tdxbz201305068

2013-05-21；

2013-09-24.

國家高技術(shù)研究發(fā)展計(jì)劃(863計(jì)劃)資助項(xiàng)目(2009AA02Z413)；國家自然科學(xué)基金資助項(xiàng)目(81271618)；教育部高等學(xué)校博士學(xué)科點(diǎn)專項(xiàng)科研基金資助項(xiàng)目(20120032110056).

高 峰（1963— ），男，教授.

高 峰，gaofeng@tju.edu.cn.

時(shí)間：2013-11-18.

http://www.cnki.net/kcms/detail/12.1127.N.20131118.1627.001.html.