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

任曉芳，徐亮

?

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

任曉芳，徐亮

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

關鍵詞：逆向問題；擴散光學層析圖像；輻射傳輸等式；擴散近似；低散射區域

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

0 引言

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

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

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

1 正向問題及逆向問題

1.1 正向問題

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

1.1.1 輻射傳輸等式

輻射傳輸等式是廣泛用于組織中光傳輸的模型[13]，RTE是傳輸等式的單速近似，因此它假設光子能量（或速度）在碰撞中不改變，介質內折射率是常量。

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

1.1.2 擴散近似

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

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

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

1.1.3輻射傳輸-擴散模型

混合模型中，域?劃分成兩個不相連的子集?rte和?da，RTE用作子域?rte中的正向模型，其中DA的假設不可行，DA用作子域?da中的正向模型，等于剩余域。混合模型的形式如公式（10）、（11）、（12）：

1.2逆向問題

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

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

已知測量數據Z時。式（15）中，F是映射吸收和散射參數到數據的光傳輸的正向模型。而且，矩陣Le是一個加權矩陣，從統計學觀點看，可解釋為噪聲協方差逆矩陣的Cholesky因子，項是正則化懲罰函數，正則化可克服不穩定性，即由病態特性造成的問題，使用一種恰當的平滑先驗模型作為正則化函數[16]。

DOT中，通常測量的光強度的動態范圍非常大，必須縮放數據，以確保優化問題的數值穩定性，且解空間的轉換可用于組成正確的預處理，縮放數據空間如公式（16）：

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

2. 數值實現

2.1 混合模型的FE-近似

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

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

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

2.2 雅可比

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

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

3 實驗與分析

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

表1　背景介質和內含物的吸收系數μ(mm-1 )和散射系數μ(mm-1 )

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

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

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

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

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

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

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

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

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

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

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

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

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

4 總結

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

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

參考文獻

[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] 張芹芹,吳曉靜,朱思偉,等.譜域光學相干層析成像量化技術及其在生物組織定量分析中的應用[J].光學精密工程,2012, 34(6):1188-1193.

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

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

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

[6] 高應俊,金重星,林林,等.基于光學相干層析成像技術的強散射介質光學特性測量[J].光子學報,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] 王彩芳.擴散光學層析成像的EM重建算法[D].北京:北京大學, 2009.

[11] 秦轉萍,趙會娟,周曉青,等.有效探測區域的內窺式漫射層析成像圖像重建算法[J].光學學報,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].光學學報,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-），女，新疆工程學院，計算機工程系，講師，碩士，研究方向：圖像處理、機器學習，烏魯木齊，830023

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

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

中圖分類號：TP391

文獻標志碼：A