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

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

時域光學(xué)拓撲成像方法的模擬和實驗研究

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

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

針對現(xiàn)有的以修正的朗伯-比爾定律(modified Lambert-Beer law，MLBL)為理論模型的光學(xué)拓撲方法存在量化度不準確的問題，提出了一種基于最小二乘擬合算法的時域光學(xué)拓撲成像方法．該方法利用時域擴散方程在半無限空間條件下的解析解和測量曲線實現(xiàn)最佳匹配以計算相鄰的源點與探測點之間的待測組織吸收系數(shù)變化值．對該時域光學(xué)拓撲成像方法進行了一系列的模擬驗證，并以基于時間相關(guān)單光子計數(shù)(time-correlated single photon counting，TCSPC)技術(shù)的多通道時間分辨測量系統(tǒng)為平臺進行了實驗驗證．結(jié)果表明：所提出的最小二乘擬合算法總體上優(yōu)于傳統(tǒng)的基于經(jīng)驗值光程的MLBL方法．

擴散方程；飛行時間；時間相關(guān)單光子計數(shù)；光學(xué)拓撲；腦功能成像

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

由于近紅外光在腦部組織穿透深度的限制，一般采用拓撲的方式研究大腦皮層應(yīng)激反應(yīng)過程．光學(xué)拓撲成像即是利用反射測量給出待測組織體表面下淺層的光學(xué)參數(shù)二維變化．與其他成像方式相比，光學(xué)拓撲成像具有靈活可變的源-探測陣列排布；由于只選擇測量與源相鄰的探測點處的反射光流，并將其與中間采樣點相關(guān)聯(lián)，因而具有較高的信噪比．

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

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

1 近紅外時域光學(xué)拓撲成像方法

1.1 基于時域擴散方程解析解的最小二乘擬合算法

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

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

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

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

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

在后面的驗證過程中，基于時域擴散方程解析解的最小二乘法用LSM表示．

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

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

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

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

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

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

2 近紅外光學(xué)拓撲成像方法模擬驗證

本文利用蒙特卡洛(Monte-Carlo，MC)模擬產(chǎn)生正向數(shù)據(jù)[2]．MC模擬的數(shù)值模型為80,mm×80 mm×30,mm 的平板．在圖1(a)所示的坐標系下，(-7.5,mm，-3.75,mm)位置處有一深度為20,mm的圓柱形凹槽，直徑D分別為 6,mm、7.5,mm、10,mm、14,mm．采用21源-4探測器配置，其中D1、D2、D3、D4為探測點位置，由此形成32個源-探測器對，對應(yīng)32個直接測量采樣點和17個間接采樣點，如圖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的增大，每種方法的量化度都會相應(yīng)提高．在所有情況下，按照量化度從大到小排列，依次是基于經(jīng)驗值光程的MLBL方法、基于計算值光程的 MLBL方法、基于時域擴散方程解析解的最小二乘擬合算法．但當(dāng)D大于源和探測器距離ρ時，基于經(jīng)驗值光程的 MLBL方法在量化度方面會出現(xiàn)過估計問題，出現(xiàn)偏差．而此時，時域測量方法的計算結(jié)果更加理想．

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

3 近紅外光學(xué)拓撲方法實驗驗證

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

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

(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)用于切換源通道，實現(xiàn)了發(fā)射光源在 16個不同的入射源點間進行切換．光開關(guān)通過FC(ferrule connection)接口與計算機進行連接和通信．用戶通過計算機VB界面控制源通路的選擇．

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

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

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

圖4 雙波長時域光學(xué)拓撲成像系統(tǒng)原理Fig.4 Schematic of the dual-wavelength，time domain optical topography system

3.2 仿體驗證

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

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

圖5 實驗示意Fig.5 Experimental setup

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

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

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

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

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

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

4 結(jié) 語

本文提出了一種在時域測量數(shù)據(jù)和時域擴散方程在平面半無限空間條件下的解析解之間尋求最佳匹配的時域光學(xué)拓撲成像方法，并進行了一系列的模擬和仿體實驗，從而得到待測組織光學(xué)參數(shù)變化的2D圖像．仿體實驗在實驗室已搭建的基于 TCSPC的多通道時域測量系統(tǒng)平臺上進行．模擬和實驗結(jié)果表明，該方法在量化度上優(yōu)于傳統(tǒng)的基于經(jīng)驗值光程的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］ 徐可欣，高 峰，趙會娟. 生物醫(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，屈軍樂. 高級時間相關(guān)單光子計數(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ā)展計劃(863計劃)資助項目(2009AA02Z413)；國家自然科學(xué)基金資助項目(81271618)；教育部高等學(xué)校博士學(xué)科點專項科研基金資助項目(20120032110056).

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

高 峰，gaofeng@tju.edu.cn.

時間：2013-11-18.

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