任曉芳，徐亮

?

利用輻射傳輸-擴散混合模型的DOT圖像重建算法

任曉芳，徐亮

摘 要：針對擴散光學(xué)層析成像中散射理論限制和逆向問題，提出了基于混合模型的擴散光學(xué)層析圖像重建算法，該算法將輻射傳輸?shù)仁胶蛿U散近似有效結(jié)合，是低擴散區(qū)域和三維熒光成像的有效果擴展。首先，研究了這種混合模型，通過最小化測量數(shù)據(jù)與混合模型解之間的最小二乘正則化估計吸收和散射分布，然后，證明了該模型可用于解決逆向問題，最后，仿真測試本文算法，給出低散射區(qū)域不同情況的重建，實驗結(jié)果表明，該算法可產(chǎn)生良好的重建效果。

關(guān)鍵詞：逆向問題；擴散光學(xué)層析圖像；輻射傳輸?shù)仁?；擴散近似；低散射區(qū)域

徐 亮（1977-），男，新疆工程學(xué)院，計算機工程系，講師，碩士，研究方向：圖像處理、機器學(xué)習(xí)等，烏魯木齊，830023

0 引言

擴散光學(xué)層析成像（Diffusion optical tomography, DOT）已在乳腺癌檢測、嬰幼兒腦組織血氧水平和功能性大腦活動研究等領(lǐng)域得到廣泛應(yīng)用[1-2]，其對于深度分辨率的影響[3]在醫(yī)學(xué)領(lǐng)域也具有很大的潛力。DOT中圖像重建問題是一種非線性病態(tài)求逆問題，因此，即使是測量或建模中的小誤差也會引起重建的大誤差[4-5]。為了克服該問題的病態(tài)性，必須使用正則化技術(shù)和貝葉斯估計[6]，且需要準(zhǔn)確描述介質(zhì)內(nèi)光傳播在計算上可行的正向模型[7]。

為了克服散射理論的限制，已經(jīng)開發(fā)出了組合擴散等式與其他模型的不同混合模型，這些模型包括：各階PN近似混合[8]，運用輻射傳輸?shù)仁脚c擴散等式[9]，輻射擴散模型[10]，可用在具有非散射區(qū)域的高度散射介質(zhì)、耦合蒙特卡羅散射方法和三階球諧函數(shù)（SP3)）作為正向模型。但是，這些模型無法解決逆向問題[11]。

本文提出一種混合模型的圖像重建算法，用于DOT求逆問題。主要目的是證明混合模型可用于逆向傳輸模型，開發(fā)混合模型作為DOT的光傳播模型，是對文獻[12]中三維熒光成像的擴展?；旌夏Ｐ椭?，用擴散等式近似無效的子域中的輻射傳輸?shù)仁浇９獾膫鞑?，擴散近似用在域中其他地方。DOT的逆向問題中，通過最小化測量數(shù)據(jù)和混合模型獲得數(shù)據(jù)之間的正則化最小二乘誤差估計物體內(nèi)吸收和散射分布。本文用設(shè)有在線搜索算法和估計參數(shù)積極性約束的高斯牛頓方法求解這個最小化問題，此外，運用數(shù)據(jù)和解空間的縮放，以提高優(yōu)化算法的性能。

1 正向問題及逆向問題

1.1 正向問題

DOT中正向問題是為了當(dāng)給定介質(zhì)光學(xué)性質(zhì)和輸入光源時計算可測量數(shù)據(jù)，生物介質(zhì)中光傳播通常通過可視作隨機算法的傳輸理論來建模，例如蒙特卡羅或確定性算法，基于用積分微分等式描述的粒子傳輸。

1.1.1 輻射傳輸?shù)仁?/p>

輻射傳輸?shù)仁绞菑V泛用于組織中光傳輸?shù)哪Ｐ蚚13]，RTE是傳輸?shù)仁降膯嗡俳疲虼怂僭O(shè)光子能量（或速度）在碰撞中不改變，介質(zhì)內(nèi)折射率是常量。

公式（1）中，c是介質(zhì)中光的速度，i是虛數(shù)單位，ω是輸入信號的角度調(diào)制頻率，和分別為介質(zhì)的散射和吸收系數(shù)，而且，是輻射，是散射項函數(shù)，是?內(nèi)光源。本文中沒有內(nèi)部光源，因此

1.1.2 擴散近似

擴散近似框架中，由下式近似輻射如公式（5）：

公式（7）中，q0( r )是?內(nèi)源，公式（7）是DA的頻域版本，邊界條件（3）不能以擴散近似的變量表示，相反，經(jīng)常向內(nèi)定向光子電流為零的近似代替它，而且，通過考慮介質(zhì)折射率和周圍介質(zhì)折射率之間的不匹配，可推導(dǎo)出一種羅賓類型的邊界條件，見文獻[1]中的例子，邊界條件可記作公式（8）：

公式中，Is是源位置εj???處的擴散邊界電流，γn是維度依賴常量，取值為γ2= 1/ π和γ3=1/4，ξ是管理邊界??上內(nèi)反射的參數(shù)[14]，在折射率匹配的情況下，ξ=1，如公式（9）：

1.1.3輻射傳輸-擴散模型

混合模型中，域?劃分成兩個不相連的子集?rte和?da，RTE用作子域?rte中的正向模型，其中DA的假設(shè)不可行，DA用作子域?da中的正向模型，等于剩余域?；旌夏Ｐ偷男问饺绻剑?0）、（11）、（12）：

1.2逆向問題

DOT的逆向問題中，基于物體表面上光傳輸測量值Z估計物體的光學(xué)參數(shù)。本文中，待估計光學(xué)參數(shù)為介質(zhì)內(nèi)吸收和散射系數(shù)(μa,μs )，在一個離散框架中考慮逆向問題的解，離散化域?為K個不相交元素?k，依據(jù)下式表示吸收和散射系數(shù)如公式（13）、（14）：

針對DOT逆向問題的傳統(tǒng)算法是基于非線性最小二乘方法，這種算法中，使用單一數(shù)據(jù)獲取估計介質(zhì)光學(xué)屬性的絕對值，正則化非線性最小二乘問題是為了估計吸收和散射分布?，其最小化函數(shù)為公式（15）：

已知測量數(shù)據(jù)Z時。式（15）中，F(xiàn)是映射吸收和散射參數(shù)到數(shù)據(jù)的光傳輸?shù)恼蚰Ｐ汀６?，矩陣Le是一個加權(quán)矩陣，從統(tǒng)計學(xué)觀點看，可解釋為噪聲協(xié)方差逆矩陣的Cholesky因子，項是正則化懲罰函數(shù)，正則化可克服不穩(wěn)定性，即由病態(tài)特性造成的問題，使用一種恰當(dāng)?shù)钠交闰災(zāi)Ｐ妥鳛檎齽t化函數(shù)[16]。

DOT中，通常測量的光強度的動態(tài)范圍非常大，必須縮放數(shù)據(jù)，以確保優(yōu)化問題的數(shù)值穩(wěn)定性，且解空間的轉(zhuǎn)換可用于組成正確的預(yù)處理，縮放數(shù)據(jù)空間如公式（16）：

因此使用幅值和相位對數(shù)作為數(shù)據(jù)，而且，解空間中，用平均值縮放吸收值和散射值為公式（17）：

2. 數(shù)值實現(xiàn)

2.1 混合模型的FE-近似

本文使用有限元方法(FEM)對混合模型的數(shù)值近似，推導(dǎo)混合模型的FE-近似，實現(xiàn)參考文獻[12-17]。其結(jié)果是，獲得下列矩陣等式如公式（19）：

公式中，A0、A1、A2、A3和A4是公式（21）~（25）：

混合模型的FE解通過求解式（23）獲得，因此，作為正向問題的解，獲得RTE子域中空間和角度離散化節(jié)點中輻射α以及DA子域節(jié)點中光子密度a。出射度，即可測量值，可使用式（4）計算為公式（26）：

2.2 雅可比

高斯牛頓算法中，每輪迭代中必須計算雅可比矩陣，考慮縮放因素，縮放雅可比形式為公式（27）~（28）：

公式中，對應(yīng)于元素?k的雅可比Jμa第k列由逐列向量化獲得公式（29）~（33）：

3 實驗與分析

使用二維仿真測試混合模型的可行性，半徑為20mm的圓域?，域邊界上，設(shè)置32個等間距源和檢測器，考慮3種不同的內(nèi)部結(jié)構(gòu)：具有吸收內(nèi)含物和散射內(nèi)含物的高度散射介質(zhì)（情況1），具有兩個低散射內(nèi)含物的高度散射介質(zhì)（情況2），低散射層內(nèi)部具有吸收內(nèi)含物和散射內(nèi)含物的高度散射介質(zhì)（情況3）。內(nèi)含物的半徑為4mm，低散射層的寬度為0.5mm，仿真的吸收和散射系數(shù)如表1所示：

表1　背景介質(zhì)和內(nèi)含物的吸收系數(shù)μ(mm-1 )和散射系數(shù)μ(mm-1 )

最小化函數(shù)式（15）重建吸收和散射分布，為了吸收和散射的表示，劃分域為個三角形，在分段常量基（13）-（14）中表示吸收和散射系數(shù)，吸收和散射的離散化網(wǎng)格如圖1所示：

圖1　光學(xué)參數(shù)的離散化（左圖），用在情況1和情況3中的正向模型的FE-離散化（中圖），情況2（右圖）中混合模型的RTE子域用黑色表示，DA子域用灰色表示

圖1（a）重建中使用的正向模型是RTE和混合模型，F(xiàn)EM用于數(shù)值實現(xiàn)，F(xiàn)E-網(wǎng)格的空間離散化包括5853個節(jié)點和11177個三角形元素，混合模型中，空間域劃分為RTE 和DA子域。情況1和情況3中，RTE子域包括的元素在距邊界5mm距離內(nèi)，情況2中，RTE子域包含的元素在距邊界12mm距離內(nèi)。以這樣一種方式選擇RTE子域，假設(shè)已知一些低散射區(qū)域位置的先驗知識，選擇覆蓋的區(qū)域以及近邊界區(qū)域的子域，剩余元素包括在DA子域內(nèi)，正向模型離散化和RTE、DA子域的劃分如圖1（b）和圖1圖（c）所示。RTE的角度離散化和混合模型中RTE部分包括32個等間距角度方向。

重建的吸收和散射如圖2至圖4所示：

圖2　具有一個吸收內(nèi)含物和一個散射內(nèi)含物的域內(nèi)吸收（第一行）和散射（最后一行）分布（情況1），圖像從左到右：仿真的分布(左列)、使用RTE獲得的重建（第二列）、混合模型作為正向模型（第三列）

圖3　具有兩個低散射內(nèi)含物的域內(nèi)吸收（第一行）和散射（最后一行）分布（情況2），圖像從左到右：仿真的分布（左列）、使用RTE獲得的重建（第二列）、混合模型作為正向模型（第三列）

圖4　低散射層內(nèi)具有一個內(nèi)含物和一個散射內(nèi)含物的域內(nèi)吸收（第一行）和散射（最后一行）分布（情況3），圖像從左到右：仿真的分布（左列）、使用RTE獲得的重建（第二列）、混合模型作為正向模型（第三列）

為了比較估計參數(shù)的準(zhǔn)確度，計算相對估計誤差。因此，取域?內(nèi)31428個點，計算吸收和散射的相對估計誤差如公式（34）、（35）：

從圖2至圖4可以看出，混合模型獲得的重建類似于RTE獲得的重建。情況1和2中，可很好區(qū)分內(nèi)含物，RTE和混合模型估計具有相同準(zhǔn)確性的參數(shù)，情況3中，很難區(qū)分位于低散射層內(nèi)部的吸收內(nèi)含物，無關(guān)正向模型，另一方面，仍然能良好區(qū)分高度散射內(nèi)含物，兩種模型給出相同幅值的相對誤差。情況3代表一種具有挑戰(zhàn)性的成像情況，其中測量僅攜帶低散射層內(nèi)部的非常有限的信息，這是由于光子的物理行為往往沿低散射區(qū)域運行。

除了估計的光學(xué)參數(shù)，本文還比較了估計的正向解，使用與參數(shù)估計誤差情況中相同的31428個點計算估計的正向數(shù)據(jù)的幅值和相位偏移對數(shù)的相對誤差，使用等式如公式（36）、（37）：

幅值和相位的估計的光子密度對數(shù)的相對誤差如表2所示：

表2　幅值△In|Ф|(%)和相位偏移△arg|Ф|(%)對數(shù)正向解的相對估計誤差

從表2可以看出，本文提出的混合模型在3種情況下的相對誤差均小于RTE，表明了本文模型具有更小的重建誤差，有助于提升重建效果。

4 總結(jié)

本文使用混合輻射傳輸-擴散模型作為光傳輸?shù)恼蚰Ｐ颓蠼釪OT的圖像重建問題，圖像重建基于正則化最小二乘算法解決高斯牛頓算法，用FEM進行數(shù)值求解混合模型，仿真測試該結(jié)果表明混合模型可用作擴散光學(xué)層析成像，而且可以成功求解逆向問題。從重建質(zhì)量看，與使用完全RTE求解器得到的重構(gòu)差不多。是三維熒光成像的擴展。

本文方法仍有改進空間，可選擇不同的算法進行正則化，如全變差正則化或稀疏正則化等，雅克比迭代也可以選擇更為優(yōu)良的自適應(yīng)步長，這將是下一步研究重點。

參考文獻

[1] Ferradal S L, Eggebrecht A T, Hassanpour M, et al. Atlas-based head modeling and spatial normalization for high-density diffuse optical tomography: In vivo validation against fMRI[J]. Neuroimage, 2014, 85(4): 117-126.

[2] 張芹芹,吳曉靜,朱思偉,等.譜域光學(xué)相干層析成像量化技術(shù)及其在生物組織定量分析中的應(yīng)用[J].光學(xué)精密工程,2012, 34(6):1188-1193.

[3] 郭昕,王向朝,步鵬,等.樣品散射對頻域光學(xué)相干層析成像光譜形狀和深度分辨率的影響[J].光學(xué)學(xué)報,2014,(1):181-186.

[4] 金重星.基于光學(xué)相干層析成像技術(shù)的生物組織散射特性的研究[D].廣州:暨南大學(xué),2011.

[5] Bal G. Inverse transport theory and applications [J]. Inverse Problems, 2009, 25(5): 1024-1029.

[6] 高應(yīng)俊,金重星,林林,等.基于光學(xué)相干層析成像技術(shù)的強散射介質(zhì)光學(xué)特性測量[J].光子學(xué)報,2011,17(1): 98-102.

[7] Tarvainen T, Kolehmainen V, Arridge S R, et al. Image reconstruction in diffuse optical tomography using the coupled radiative transport-diffusion model[J]. Journal of Quantitative Spectroscopy and Radiative Transfer, 2011, 112(16): 2600-2608.

[8] Surya Mohan P, Tarvainen T, Schweiger M, et al. Variable order spherical harmonic expansion scheme for the radiative transport equation using finite elements[J]. Journal of Computational Physics, 2011, 230(19): 7364-7383.

[9] Gorpas D, Yova D, Politopoulos K. A three-dimensional finite elements approach for the coupled radiative transfer equation and diffusion approximation modeling in fluorescence imaging[J]. Journal of Quantitative Spectroscopy and Radiative Transfer, 2010, 111(4): 553-568.

[10] 王彩芳.擴散光學(xué)層析成像的EM重建算法[D].北京:北京大學(xué), 2009.

[11] 秦轉(zhuǎn)萍,趙會娟,周曉青,等.有效探測區(qū)域的內(nèi)窺式漫射層析成像圖像重建算法[J].光學(xué)學(xué)報,2011, (24)11:548-554.

[12] Gorpas D, Yova D, Politopoulos K. A three-dimensional finite elements approach for the coupled radiative transfer equation and diffusion approximation modeling in fluorescence imaging [J]. Journal of Quantitative Spectroscopy and Radiative Transfer, 2010, 111(4): 553-568.

[13] Louedec K, Urban M. Ramsauer approach for light scattering on nonabsorbing spherical particles and application to the Henyey-Greenstein phase function[J]. Applied optics, 2012, 51(32): 7842-7852.

[14] Kaipio J, Somersalo E. Statistical and computational inverse problems[M]. New York: Springer, 2005.

[15] Habermehl C, Holtze S, Steinbrink J, et al. Somatosensory activation of two fingers can be discriminated with ultrahigh-density diffuse optical tomography[J]. Neuroimage, 2012, 59(4): 3201-3211.

[16] 馬文娟,高峰,張偉,等.基于輻射傳輸方程三階簡化球諧近似模型的時域熒光擴散層析圖像重建方法[J].光學(xué)學(xué)報,2011,19(5): 541-547.

[17] Tarvainen T, Vauhkonen M, Kolehmainen V, Kaipio J. Finite element model for the coupled radiative transfer equation and diffusion approximation [J]. Int J Numer Meth Eng 2006, 65(3):383-405.

Image Reconstruction Algorithm of DOT Based on Radiation Transport-Diffusion Hybrid Model

Ren Xiaofang, Xu Liang
(Department of Computer Engineering, Xinjiang Institute of Engineering, Urumqi 830023, China)

Abstract:For the limitations of scattering theory and the inverse problem in diffuse optical tomography, a hybrid model to reconstruct the diffusion tomographic images is proposed. This algorithm combines the radiation transport equation and diffusion approximation effectively. It is an extension of low diffusion regions and three-dimensional fluorescence imaging. Firstly this paper researches the hybrid model, absorption and scattering distributions are estimated by minimizing a regularized least-squares error between the measured data and solution of the coupled model. Then proves the model can be used to solve the inverse problem. Finally gives the simulation to test the algorithm and the reconstruction results of different low diffusion regions. Experimental results show the algorithm can produce good reconstruction results.

Key words:Inverse Problems; Diffusion Optical Tomography Image; Radiation Transport Equation; Diffusion Approximation; Low Diffusion Regions

收稿日期：（2015.06.27）

作者簡介：任曉芳（1979-），女，新疆工程學(xué)院，計算機工程系，講師，碩士，研究方向：圖像處理、機器學(xué)習(xí)，烏魯木齊，830023

基金項目：新疆維吾爾自治區(qū)自然科學(xué)基金項目（No.2013211A031）；新疆工程學(xué)院人文基金項目（2014xgy311612）

文章編號：1007-757X(2016)02-0020-05

中圖分類號：TP391

文獻標(biāo)志碼：A