利用辐射传输-扩散混合模型的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