近紅外腦功能DOT成像雙網(wǎng)格半三維算法

劉 明，孟 偉，周曉青，賈夢宇，趙會娟,

近紅外腦功能DOT成像雙網(wǎng)格半三維算法

劉 明1，孟 偉1，周曉青1，賈夢宇1，趙會娟1,2

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

采用基于光子輸運(yùn)模型的擴(kuò)散層析成像(DOT)實(shí)現(xiàn)在近紅外腦功能成像的諸多方案具有良好的定量研究潛力，但由于其逆問題固有的病態(tài)性，傳統(tǒng)的全三維重建存在精度和分辨率低以及計(jì)算量大等缺點(diǎn)，因此根據(jù)實(shí)際問題約束和簡化重建過程對改善成像性能具有重要意義．根據(jù)人體頭部的解剖結(jié)構(gòu)及近紅外光在人腦中的傳輸特性，提出了改善逆問題病態(tài)性的雙網(wǎng)格半三維圖像重建算法，即在正問題中采用密集網(wǎng)格剖分下的全三維有限元計(jì)算提高計(jì)算精度，在逆問題中采用稀疏網(wǎng)格剖分且在同一腦組織層進(jìn)行二維重建方法．模擬結(jié)果表明，對于單目標(biāo)重建，雙網(wǎng)格半三維重建算法與全三維重建算法相比，量化度提高30%、重建速度快4～12倍；對于雙目標(biāo)重建在信噪比30,dB以上的情況，半三維重建算法在CCS＝12.5,mm的空間分辨能力與全三維重建算法CCS＝20,mm的空間分辨能力相當(dāng)．

腦功能成像；半三維算法；有限元；雙網(wǎng)格

近紅外(near infrared，NIR)腦功能成像是一種研究大腦血液動力學(xué)變化的非侵入式技術(shù)．該技術(shù)的優(yōu)勢體現(xiàn)在：①其具有無創(chuàng)、無輻射、高的時(shí)間分辨率、價(jià)廉可移動等特點(diǎn)，對于大腦的思維和意識等活動研究有著重要意義；②因其可在自然條件下進(jìn)行監(jiān)測，NIR腦功能成像具有很好的運(yùn)動魯棒性，該技術(shù)對新生兒和嬰幼兒腦功能研究有著重要意義[1-2]．此外，近紅外光學(xué)成像技術(shù)還可用于診斷腦組織的某些生理異常，如血腫、腦室溢血、局部缺氧和缺血等．

目前商品化的近紅外腦功能成像儀主要基于拓?fù)涑上?optical topography，OT)原理[3-4]．該成像方式假設(shè)探測區(qū)域下深度和橫向方向的組織光學(xué)參數(shù)均勻分布，并利用修正的朗伯-比爾定律給出該區(qū)域組織的平均血氧濃度變化，因此OT可被認(rèn)為是2D成像．雖然基于OT理論的NIR腦功能成像儀具有圖像重建速度快的優(yōu)點(diǎn)，但是具有3方面局限：①NIR腦功能成像儀是透過頭皮、頭骨來探測發(fā)生在灰質(zhì)內(nèi)的腦血液動力學(xué)變化，而OT假設(shè)源-探測點(diǎn)之間的組織光學(xué)參數(shù)為均勻分布，勢必導(dǎo)致計(jì)算結(jié)果無法真實(shí)反映灰質(zhì)內(nèi)的血液動力學(xué)變化[5]；②由于OT不具有深度分辨能力，因此，心跳、血壓等正常生理活動將通過頭皮的血液循環(huán)間接影響NIR腦功能儀的準(zhǔn)確度[6]；③由于修正的朗伯-比爾定律不能準(zhǔn)確描述光在組織中傳播的真實(shí)過程，因此限制了NIR腦功能儀的空間分辨率和量化度的進(jìn)一步提高[7]．

近年來，應(yīng)用擴(kuò)散層析成像(diffuse OT，DOT)理論進(jìn)行腦功能成像的研究得到了越來越多的關(guān)注[8-10]. 在反射型DOT中，光在腦內(nèi)傳播是三維的，全三維重建可以得到任意深度內(nèi)腦血氧的變化，但是由于成像區(qū)域大而導(dǎo)致圖像重建速度過慢．同時(shí)，由于DOT重建中正問題和逆問題的計(jì)算中需要對成像空間進(jìn)行離散，全三維重建將導(dǎo)致測量數(shù)據(jù)遠(yuǎn)遠(yuǎn)少于重建參數(shù)個(gè)數(shù)，因此逆問題的欠定性將十分嚴(yán)重．

對于實(shí)際的腦功能成像，人頭可分為頭皮、頭骨、灰質(zhì)和白質(zhì)4部分．由于頭皮層較薄且頭骨對光的吸收大，因此可將頭皮和頭骨設(shè)為主要采用頭骨的光學(xué)參數(shù)的均一的組織層．腦功能引起的腦血氧的變化主要發(fā)生在灰質(zhì)中，由于探測系統(tǒng)信噪比的限制，源-探距離是有限的(一般為3,cm左右)，此時(shí)NIR光只能到達(dá)的有限深度基本在灰質(zhì)表層，因此可忽略灰質(zhì)內(nèi)光學(xué)參數(shù)沿深度的變化，同時(shí)將白質(zhì)設(shè)為均一的組織層．綜上所述，同時(shí)考慮到光子在非均勻組織中的三維傳輸，本文提出半三維圖像重建算法，以彌補(bǔ)二維拓?fù)涑上窈腿S圖像重建的不足.

1 理論與方法

DOT系統(tǒng)可采用時(shí)域、頻域、連續(xù)光(穩(wěn)態(tài))3種測量模式．由于連續(xù)光模式提供了簡單的系統(tǒng)、快速的數(shù)據(jù)獲取方式和較高的信噪比，因此在腦功能成像研究中被廣泛地采用．本文對光在頭內(nèi)傳輸?shù)恼龁栴}采用穩(wěn)態(tài)擴(kuò)散方程并利用有限元法進(jìn)行求解[11]．在圖像重建中采用非線性最優(yōu)化Newton-Raphson算法．在忽略腦功能活動引起的約化散射系數(shù)(sμ′)的變化時(shí)，逆問題求解簡化為計(jì)算灰質(zhì)內(nèi)吸收系數(shù)(aμ)的微擾aμδ．采用迭代策略使吸收系數(shù)由相對粗略的初始估計(jì)逐步逼近真值[12-13]，即

式中：χ為邊界光子密度測量值；F(·)為正問題求解的邊界光子密度值；Jμa為關(guān)于μa的Jacobian矩陣；上標(biāo)i代表第i次迭代；δμi為第i次迭代的吸收系數(shù)

a微擾．

在全三維擴(kuò)散層析成像圖像重建中，要求重建組織體的某一個(gè)橫斷面甚至整個(gè)三維空間的光學(xué)參數(shù)圖像，此時(shí)關(guān)于μa的Jacobian矩陣元素計(jì)算式為

式中：j為有限元形函數(shù)體元編號；Ω為成像區(qū)域；Γ和Φ為光子能量密度；ξd為探測器所在位置，d＝1，2，…，D；ζs為第s個(gè)光源所在位置，s＝1，2，…，S；c為光在組織中的傳播速度；r為區(qū)域內(nèi)坐標(biāo)矢量；uxyz(r)為有限元形函數(shù)．Γ(ξd,ζs)為邊界光子能量密度，其中α是與邊界條件有關(guān)的系數(shù)，在Robin條件下，fR為擴(kuò)散傳輸內(nèi)反射系數(shù)．

在本文提出的近紅外擴(kuò)散光半三維成像算法中，在正問題中，采用三維有限元法計(jì)算擴(kuò)散方程，得到光子密度三維空間分布；在逆問題中，考慮在有限的源-探距離下漫反射光測量穿透深度的限制，重建過程采用半三維方式，在頭部每一解剖結(jié)構(gòu)層內(nèi)均假設(shè)光學(xué)參數(shù)為二維分布，即該層沿深度方向光學(xué)參數(shù)相等，沿重建面方向光學(xué)參數(shù)為二維分布．由于NIR腦功能成像儀探測的是發(fā)生在灰質(zhì)內(nèi)的血氧變化，因此可進(jìn)一步假設(shè)頭皮、頭骨和白質(zhì)分別是均一物質(zhì)，如圖1所示．基于上述思想，近紅外擴(kuò)散半三維Jacobian矩陣元素計(jì)算方法可表示為

layer平面內(nèi)新的有限元形函數(shù)．

圖1 人頭簡化3層模型與正、逆問題中的網(wǎng)格剖分Fig.1 Simplified view of three-lay head model，dense mesh for forward problem and sparse mesh for inverse problem

比較式(2)和式(3)可見，與全三維算法相比，采用半三維算法可減少逆問題重建節(jié)點(diǎn)數(shù)目，從而縮小了Jacobian矩陣規(guī)模，節(jié)約了計(jì)算時(shí)間．同時(shí)，由于重建節(jié)點(diǎn)數(shù)目變少，逆問題的病態(tài)性減弱，有望提高重建的量化度．

在近紅外擴(kuò)散光半三維重建算法中，對式(1)進(jìn)行有限元離散化，計(jì)算過程中采用雙網(wǎng)格剖分方式，如圖1所示．在正問題中應(yīng)用密集網(wǎng)格(見圖1(b))的剖分計(jì)算節(jié)點(diǎn)光子能量密度，使漫射方程光子能量密度的解更精確．在光學(xué)參數(shù)重建過程中，根據(jù)半三維算法假設(shè)應(yīng)用稀疏網(wǎng)格(見圖1(c))進(jìn)行計(jì)算，有利于減少重建中過多的自由度．

老齡化為養(yǎng)老保障帶來巨大壓力。隨著人口平均預(yù)期壽命的上升，越來越多的人開始享受養(yǎng)老保障的福利，但上交老年保險(xiǎn)的比例不是很高，導(dǎo)致保險(xiǎn)行業(yè)盈利受阻，甚至負(fù)盈利[5]。人口老齡化會對醫(yī)療保險(xiǎn)和老年保險(xiǎn)形成沖擊，不能有效解決這個(gè)問題將給我國經(jīng)濟(jì)社會帶來負(fù)面影響。

采用稀疏網(wǎng)格分布，aμJ矩陣元素的上、下標(biāo)分別表示源-探位置和各節(jié)點(diǎn)的坐標(biāo)，則其Jacobian矩陣可以寫為

式中X、Y、Z為3個(gè)坐標(biāo)方向的坐標(biāo)最大值．

考慮到半三維假設(shè)利用圖1(a)中3層簡化模型，則相應(yīng)于式(4)同一層內(nèi)的Jacobian矩陣元素可以表示為

同時(shí)，針對半三維重建，式(1)中的吸收系數(shù)微擾δμa可以用Ωxy平面內(nèi)吸收系數(shù)微擾δμa代替．

2 模擬計(jì)算結(jié)果與分析

為驗(yàn)證本文所發(fā)展的雙網(wǎng)格半三維算法(簡稱半三維算法)并與全三維算法進(jìn)行對比，采用圖2所示模型進(jìn)行模擬驗(yàn)證．設(shè)模型體積為60,mm ×60,mm× 50,mm ，第1層厚度為5,mm，用來模擬頭皮和頭骨，背景光學(xué)吸收系數(shù)aμ設(shè)為0.01,mm-1，約化散射系數(shù)sμ′設(shè)為1,mm-1[14]；第2層厚度為10,mm，用來模擬灰質(zhì)背景，aμ設(shè)為0.02,mm-1，sμ′設(shè)為1,mm-1；第3層厚度為35,mm，用來模擬白質(zhì)，其光學(xué)參數(shù)與灰質(zhì)相同．設(shè)16個(gè)光源點(diǎn)和16個(gè)探測點(diǎn)以4×4方式均勻分布在頂平面上的中心，相鄰最近兩點(diǎn)之間相距10,mm，且探測點(diǎn)和源點(diǎn)重合．每個(gè)光源分別作用時(shí)，其余15個(gè)探測點(diǎn)同時(shí)接收對應(yīng)位置的漫射光．重建過程應(yīng)用Matlab和COMSOL Mutiphysics軟件混合編程環(huán)境實(shí)現(xiàn)，重建目標(biāo)為z＝5,mm處灰質(zhì)層橫切面圖像．在本文中，產(chǎn)生光子密度“測量值”和圖像重建中的正問題所用的網(wǎng)格剖分相同．

圖2 數(shù)值仿真所用的重建模型Fig.2Reconstruction models with one target and withtwo targets for numerical simulation

2.1 重建尺寸精度、量化度和速度分析

由于正問題的計(jì)算精度和網(wǎng)格剖分密度有關(guān)，為了獲得合理的剖分密度，通過COMSOL Mutiphysics軟件對圖1(b)所示正問題模型采用不同剖分方式．在剖分過程中按照COMSOL Mutiphysics軟件提供的“單元尺寸參數(shù)”為標(biāo)準(zhǔn)確定剖分網(wǎng)格疏密程度，分別產(chǎn)生采用極密集剖分時(shí)的35,010個(gè)節(jié)點(diǎn)以及采用密集剖分時(shí)的9,137個(gè)節(jié)點(diǎn)．采用圖2(a)所示的單目標(biāo)模型比較全三維算法和本文發(fā)展的雙網(wǎng)格半三維算法的重建尺寸精度、量化度和速度．其中目標(biāo)體為位于模型第2層中部的圓柱，半徑為5,mm，吸收系數(shù)分別設(shè)為0.03、0.04、0.05,mm-1，用以模擬發(fā)生在灰質(zhì)中的腦血氧變化．

計(jì)算微擾時(shí)正則化方法選用代數(shù)重建技術(shù)(ART)，松弛因子取為0.1[15]．在重建計(jì)算過程中通過迭代誤差判斷探測點(diǎn)上重建值,sdχ與“測量值”Fs,d(μa)的接近程度，誤差定義為ε=0.03的迭代次數(shù)為迭代截止限制次數(shù)，其中newε為本次迭代誤差，preε為前一次迭代誤差．經(jīng)多次模擬發(fā)現(xiàn)，在極密集網(wǎng)格剖分方式下，全三維算法迭代截止限制次數(shù)在20～30次之間．雙網(wǎng)格半三維算法迭代截止限制次數(shù)在10～15次之間．

為了進(jìn)一步比較所獲得的圖像，采用量化度(QR)、重建誤差(QE)和半高寬(FWHM)對獲得的結(jié)果進(jìn)行評估．量化度定義為重建目標(biāo)吸收系數(shù)和背景吸收系數(shù)之差與真實(shí)值之比[16]；重建誤差定義為重建吸收系數(shù)最大值與目標(biāo)吸收系數(shù)最大值之差的絕對值與目標(biāo)吸收系數(shù)之比；半高寬定義為某一方向上重建圖像的形貌曲線的半高寬度．

圖3(a)、(c)、(e)分別為密集網(wǎng)格剖分下對不同目標(biāo)體采用全三維算法的重建圖像，圖3(b)、(d)、(f)分別為密集網(wǎng)格剖分下對不同目標(biāo)體采用雙網(wǎng)格半三維算法重建的圖像，其中虛線為目標(biāo)的真實(shí)大小及位置．圖4為全三維和半三維算法對z＝5,mm處重建的圖像在y＝30,mm時(shí)沿x方向重建值的形貌曲線，兩種算法及真實(shí)目標(biāo)分別用不同線型表示．圖5和圖6分別為極密集網(wǎng)格剖分下的形貌曲線及重建圖像．

圖3 采用密集網(wǎng)格下全三維算法與雙網(wǎng)格半三維算法獲得的重建圖像Fig.3Reconstructed images by using 3D and dual mesh semi-3D algorithms under finer mesh for FEM calculation

對于不同吸收系數(shù)的目標(biāo)體，采用密集網(wǎng)格與極密集網(wǎng)格剖分下全三維與半三維算法重建所得的結(jié)果如表1所示．從重建尺寸FWHM的角度，由于目標(biāo)的真實(shí)直徑為10,mm，半三維算法重建出的目標(biāo)尺寸更接近真實(shí)值；從量化度QR角度，無論密集剖分還是極密集剖分，應(yīng)用雙網(wǎng)格半三維算法得到的量化度均高于全三維算法；從量化誤差QE角度，在極密集網(wǎng)格下半三維算法的量化誤差明顯優(yōu)于三維算法；從計(jì)算時(shí)間角度，在密集節(jié)點(diǎn)數(shù)下，采用半三維算法迭代1次的時(shí)間約為57,s，而采用全三維算法約需91,s；在極密集節(jié)點(diǎn)數(shù)下，采用半三維算法完成一次迭代的時(shí)間約為2,min，而采用全三維算法約需12,min．達(dá)到相同迭代截止限制條件，全三維算法的迭代次數(shù)約為雙網(wǎng)格半三維算法的2倍．因此在極密集網(wǎng)格剖分下，雙網(wǎng)格半三維算法將重建速度提高了約12倍．本模擬結(jié)果表明：對于大量有限元節(jié)點(diǎn)正問題的圖像重建，半三維算法在量化度、計(jì)算時(shí)間和重建尺寸精度上明顯優(yōu)于全三維算法；對于半三維重建算法，采用極密集網(wǎng)格剖分比密集網(wǎng)格剖分獲得更高質(zhì)量的圖像，因此在下文中均采用極密集網(wǎng)格剖分.

圖4 采用密集網(wǎng)格下全三維算法與雙網(wǎng)格半三維算法獲得的形貌曲線Fig.4 Profile by using 3D and dual mesh semi-3D algorithms under finer mesh for FEM calculation

圖5 采用極密集網(wǎng)格下全三維與雙網(wǎng)格半三維算法獲得形貌曲線Fig.5Profile by using 3D and dual mesh semi-3D algorithms under extra fine mesh for FEM calculation

圖6 采用極密集網(wǎng)格下全三維與雙網(wǎng)格半三維算法獲得的重建圖像Fig.6 Reconstructed images by using 3D and dual mesh semi-3D algorithms under extra fine mesh for FEM calculation

表1 密集網(wǎng)格與極密集網(wǎng)格剖分下全三維與雙網(wǎng)格半三維算法比較Tab.1 Comparison of 3D and dual mesh semi-3D reconstruction algorithms with finer mesh and extra fine mesh

2.2 空間分辨率與抗噪性分析

為進(jìn)一步對比全三維算法和雙網(wǎng)格半三維算法的空間分辨能力和抗噪性，采用圖2(b)所示的雙目標(biāo)模型，兩目標(biāo)半徑均為5,mm，吸收系數(shù)均為0.03,mm-1．

考慮到實(shí)際應(yīng)用中光電倍增管等探測器會引起熱噪聲和散粒噪聲等，向測量信號加入噪聲形成新的“測量值”，如式(6)所示．

式中：Γ為無噪聲的“測量值”；SNR為信噪比；Gnoise為高斯白噪聲．

改變目標(biāo)之間的中心距(CCS)使兩目標(biāo)不斷靠近，令CCS分別為20.0、15.0、12.5,mm，即xy面上雙目標(biāo)最近點(diǎn)距離分別為10.0、5.0、2.5,mm，對模型進(jìn)行圖像重建，為了定量評估重建圖像的空間分辨性，定義，即R＝1時(shí)，兩物體可被完全區(qū)分開，一般人眼可將兩物體分開時(shí)的R＞0.1．

圖7和圖8分別為固定信噪比SNR＝40,dB下，CCS分別為20.0,mm和12.5,mm時(shí)全三維算法和雙網(wǎng)格半三維算法獲得的重建圖像和形貌曲線．其中圖8(a)、(b)分別為CCS＝20.0,mm和12.5,mm時(shí)全三維和半三維算法對z＝5,mm處重建的圖像在y＝30,mm時(shí)沿x方向重建值的形貌曲線，兩種算法及真實(shí)目標(biāo)分別用不同線型表示．

圖9和圖10分別為固定CCS＝20.0,mm下，SNR分別為35,dB和30,dB時(shí)全三維算法和雙網(wǎng)格半三維算法獲得的重建圖像和形貌曲線．其中圖

圖7 信噪比SNR＝40,dB，CCS分別為20.0,mm和12.5,mm時(shí)全三維和雙網(wǎng)格半三維算法獲得的重建圖像Fig.7Reconstructed images when SNR is 40,dB and CCS is 20.0，12.5,mm respectively with 3D and dual mesh semi-3D algorithms

10(a)、(b)分別為SNR＝35,dB和30,dB時(shí)全三維和半三維算法對z＝5,mm處重建的圖像在y＝30,mm時(shí)沿x方向重建值的形貌曲線．在不同SNR和CCS下，采用全三維與半三維算法重建所得空間分辨率R如表2所示．雙網(wǎng)格半三維算法的R值在相同CCS和不同SNR下都大于全三維算法的R值．在添加3個(gè)噪聲的情況下，噪聲增大，R降低，不易被分辨，并且影響重建目標(biāo)形貌．在CCS為12.5,mm時(shí)，全三維算法的R已經(jīng)小于0.1，即無法分辨雙目標(biāo)；而此時(shí)雙網(wǎng)絡(luò)半三維算法R≥0.3可對雙目標(biāo)進(jìn)行分辨．重建過程的抗噪性與正則化方法有關(guān)，ART方法松弛因子的大小取決于測量信號的信噪比．如在松弛因子為0.1時(shí)，SNR小于30,dB則圖像形貌較差．對比圖10形貌曲線，雙網(wǎng)格半三維算法在SNR＝35,dB時(shí)的形貌優(yōu)于SNR＝30,dB；而全三維算法形貌曲線變化不明顯．說明全三維算法在SNR為30,dB時(shí)抗噪性優(yōu)于雙網(wǎng)格半三維算法．結(jié)合算法的應(yīng)用場合，在較高信噪比(SNR≥35,dB)的儀器系統(tǒng)，選用雙網(wǎng)格半三維算法可獲得更高的空間分辨率．

圖8 信噪比SNR＝40,dB，CCS分別為20.0,mm和12.5,mm時(shí)全三維和雙網(wǎng)格半三維算法獲得的形貌曲線Fig.8 Profile when SNR is 40,dB and CCS is 20.0，12.5,mm respectively with 3D and dual mesh semi-3D algorithms

圖9 信噪比SNR分別為35、30,dB，CCS為20.0,mm時(shí)全三維和雙網(wǎng)格半三維算法獲得的重建圖像Fig.9 Reconstructed images when SNR is 35，30,dB respectively and CCS is 20.0,mm with 3D and dual mesh semi-3D algorithms

圖10 信噪比SNR分別為35、30,dB，CCS為20.0,mm時(shí)全三維和雙網(wǎng)格半三維算法獲得的形貌曲線Fig.10Profile when SNR is 35，30,dB respectively and CCS is 20.0,mm with 3D and dual mesh semi-3D algorithms

表2 全三維與雙網(wǎng)格半三維算法空間分辨率RTab.2Comparison of R between 3D and dual mesh semi-3D algorithm

3 結(jié) 語

根據(jù)NIR腦功能成像特點(diǎn)，本文提出了近紅外擴(kuò)散光雙網(wǎng)格半三維重建算法．模擬結(jié)果表明，本算法與全三維算法相比，由于對腦功能問題進(jìn)行合理的假設(shè)，重建過程中采用雙網(wǎng)格，在逆問題重建過程中有效地降低了重建節(jié)點(diǎn)數(shù)目，量化度明顯優(yōu)于全三維算法；同時(shí)節(jié)約了計(jì)算開銷，提高了運(yùn)算速度，對于網(wǎng)格節(jié)點(diǎn)數(shù)目越大的圖像重建效果越明顯，在正問題采用極密集網(wǎng)格剖分條件下本方法比全三維算法運(yùn)算速度提高了約12倍．對雙目標(biāo)重建和添加噪聲結(jié)果分析表明，半三維算法能夠?qū)﹄p目標(biāo)CCS≥12.5 mm、信噪比30,dB以上情況良好分辨．為了進(jìn)一步驗(yàn)證本研究方法的實(shí)用性，實(shí)際腦組織功能成像實(shí)驗(yàn)將在后續(xù)進(jìn)行研究．

［1］ Hebden J C，Gibson A，Yusof R M，et al. Threedimensional optical tomography of the premature infant brain[J]. Physics in Medicine and Biology，2002，47：4155-4166.

［2］ Koizumi H，Yamamoto T，Maki A，et al. Optical topography：Practical problems and new applications[J]. Applied Optics，2003，42(16)：3054-3062.

［3］ Huppert T J，Diamond S F，F(xiàn)ranceschini M A，et al. HomER：A review of time-series analysis methods for near-infrared spectroscopy of the brain[J]. Applied Optics，2009，48(10)：280-298.

［4］ Hoshi Y，Huang J，Kohri S，et al. Recognition of human emotions from cerebral blood flow changes in the frontal region：A study with event related near-infrared spectroscopy[J]. Journal of Neuroimaging，2011，21(2)：94-101.

［5］ Yamada T，Umeyama S，Matsuda K，et al. Multidistance probe arrangement to eliminate artifacts in functional near-infrared spectroscopy[J]. Journal of Biomedical Optics，2009，14(6)：064034.

［6］ N?si T，M?ki H，Hiltunen P，et al. Effect of taskrelated extracerebral circulation on diffuse optical tomography：Experimental data and simulations on the forehead[J]. Biomedical Optics Express，2013，4(3)：412-426.

［7］ Virtanen J，Noponen T，Meril?inen P，et al. Comparison of principal and independent component analysis in removing extracerebral interference from near-infrared spectroscopy signals[J]. Journal of Biomedical Optics，2009，14(5)：054032.

［8］ Franceschini M A，Joseph D K，Huppert T J，et al. Diffuse optical imaging of the whole head[J]. Journal of Biomedical Optics，2006，11(5)：054007.

［9］ Niu H，Guo P，Jiang T，et al. Improving image quality of diffuse optical tomography with a projection-errorbased adaptive regularization method[J]. Optics Express，2008，16(17)：12423-12434.

［10］ Perdue K L，F(xiàn)ang Q，Diamond S G.Quantitative assessment of diffuse optical tomography sensitivity to the cerebral cortex using a whole-head probe[J]. Physics in Medicine and Biology，2012，57(10)：2857-2872.

［11］ Arridge S R. Optical tomography in medical imaging[J]. Inverse Problems，1999，15(2)：41-93.

［12］ 秦轉(zhuǎn)萍，趙會娟，楊彥雙. 內(nèi)窺式近紅外頻域漫射層析成像圖像重構(gòu)算法[J]. 天津大學(xué)學(xué)報(bào)，2012，45(5)：423-429.

Qin Zhuanping，Zhao Huijuan，Yang Yanshuang. Image reconstruction method for endoscopic near-infrared frequency-domain diffuse optical tomography[J]. Journal of Tianjin University，2012，45(5)：423-429(in Chinese).

［13］ 金 蒙，高 峰，馬文娟，等. 基于自然邊界條件求解輻射傳輸方程的方法[J]. 天津大學(xué)學(xué)報(bào)，2011，44(9)：816-822.

Jin Meng，Gao Feng，Ma Wenjuan，et al. A finite difference method to solve radiative transfer equation with natural boundary condition[J]. Journal of Tianjin University，2011，44(9)：816-822(in Chinese).

［14］ Tian Fenghua，Niu Haijing，Khan Bilal，et al. Enhanced functional brain imaging by using adaptive filtering and a depth compensation algorithm in diffuse optical tomography[J]. IEEE Transactions on Medical Imagine，2011，30(6)：1239-1251.

［15］ 馬文娟，高 峰，張 偉，等. 基于輻射傳輸方程三階簡化球諧近似模型的時(shí)域熒光擴(kuò)散層析圖像重建方法[J]. 光學(xué)學(xué)報(bào)，2011，31(5)：0510001.

Ma Wenjuan，Gao Feng，Zhang Wei，et al. Image reconstruction method of time-domain fluorescence diffuse optical tomography based on the third-order simplified spherical harmonics approximation to radiative transfer equation[J]. Acta Optica Sinica，2011，31(5)：0510001(in Chinese).

［16］ 周曉青，王 倩，范 穎，等. 面向大數(shù)據(jù)量的擴(kuò)散光層析成像快速重建算法[J]. 光學(xué)學(xué)報(bào)，2013，33(4)：0417001.

Zhou Xiaoqing，Wang Qian，F(xiàn)an Ying，et al. Fast reconstruction scheme of diffuse optical tomography for dense sampling dataset[J]. Acta Optica Sinica，2013，33(4)：0417001(in Chinese).

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

Dual Mesh Semi-3D Algorithm Towards NIR Diffuse Optical Tomographic Imaging of Brain Function

Liu Ming1，Meng Wei1，Zhou Xiaoqing1，Jia Mengyu1，Zhao Huijuan1,2
(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)

Diffuse optical tomography(DOT)based on the photon transport model provides good quantitative potential for near infrared cerebral functional imaging. Due to the inherent ill-posed inverse problem，the disadvantages of the traditional 3D reconstruction algorithm are low-precision，low-resolution and a big burden in calculation time. Therefore，the simplification and constraint of the reconstruction procedure according tothe practical situation have significant meaning. Based on the anatomic structure of head and the diffusive natureof near infrared light propagation in the brain，the semi-3D image reconstruction algorithm was proposed to alleviate the ill-posed problem. Namely，the dense mesh of 3D finite element method is used to improve the calculation accuracy in forward problem，and the sparse mesh of 2D image reconstruction algorithm is used in inverse problem by assuming that the optical properties in grey matter are invariable along the depth. Simulation results show that，for the single target sample，the quantitativeness ratio from the semi-3D reconstruction algorithm is 30%higher than that from the 3D reconstruction algorithm，and the reconstruction speed by using the semi-3D reconstruction algorithm is 4—12,times faster than that by using the 3D reconstruction algorithm. Additionally，for the two targets sample，the spatial resolution of the semi-3D reconstruction algorithm at CCS of 12.5,mm is almost the same as that ofthe 3D reconstruction algorithm at CCS of 20,mm with the value of SNR≥30.

cerebral functional imaging；semi-3D reconstruction algorithm；finite element method；dual mesh

TK448.21；Q63

A

0493-2137(2015)01-0062-08

10.11784/tdxbz201307070

2013-07-22；

2013-12-10.

國家自然科學(xué)基金資助項(xiàng)目(81101106，61108081，81271618)；教育部高等學(xué)校博士學(xué)科點(diǎn)專項(xiàng)科研基金資助項(xiàng)目(20120032110056).

劉 明（1981— ），男，博士研究生.

趙會娟，huijuanzhao@tju.edu.cn.

時(shí)間：2014-01-07.

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