基于乘子法的自发荧光断层成像重建方法技术

一种基于乘子法的自发荧光断层成像重建方法，使用有限元分析方法扩散方程进行离散化，基于L1范数的惩罚项建立无约束条件最优化问题模型；得到所述无约束条件最优化问题模型的对偶模型；建立所述对偶模型的增广拉格朗日函数；简化增广拉格朗日函数的最大值函数；使用截断牛顿法求解增广拉格朗日函数的最大值；将增广拉格朗日函数的梯度作为目标向量的最快下降方向对目标向量进行更新；更新惩罚向量；计算目标函数值J(w)，如果||(J(w)k-J(w)k-1)||/||Φm||≥tol为真，计算k＝k+1并跳至步骤S4，否则，结束计算，其中，tol为目标函数的收敛效率阈值。本发明专利技术能够快速地在较大成像区域内得到准确可靠的光源分布信息，除正则化参数的以外其他参数都可以实现自适应调整，从而提高了成像的鲁棒性。

【技术实现步骤摘要】

本专利技术属于光学分子影像领域，涉及自发荧光断层成像重建方法，尤其是一种。
技术介绍
光学分子影像技术是在分子细胞水平上实现生物体生理、病理变化的实时、无创、在体成像。荧光断层成像作为一种重要的光学分子影像成像模态，具有高灵敏度、低费用、无放射性等诸多优点，已经在新型药物研发，在体观测肿瘤发生、发展、转移机理等多个领域获得了广泛的应用。由于在近红外和可见光波长范围里，光子在生物体 组织中传播时会发生复杂的散射，因而二维平面自发荧光成像往往不能提供自发荧光光源的准确位置，而且具有无法获得荧光光源深度信息的固有缺陷。自发荧光断层成像重建技术通过测量边界表面的突光光强分布，实现了突光光源的三维精确定位。自发荧光断层成像重建是一个典型的病态问题。而且通常，表面荧光光强非常微弱，导致CCD探测器所采集到的信号混有大量干扰噪声，从而进一步增加了荧光光源重建的复杂性，限制了自发荧光成像技术在实际中的应用。研发一种能够准确、快速、鲁棒的荧光光源重建技术一直是国际上一个研究热点，但是到目前为止，还没有一种公认的技术能够实现这一目标。现有的重建方法大都是将荧光光源重建看作是求解最优化问题。正则化方法是一个很通用的方法。其中，最常见的方法是使用基于L2范数的正则化方法，但是这类方法得到的重建结果往往过于平滑，精度也不够理想。考虑到在实际应用中，如肿瘤的早期检测和观察，荧光光源分布具有稀疏性分布的特点，因此，将稀疏性特点作为先验信息融入到优化问题模型中可以增强荧光光源重建的质量。因此基于LI范数的正则化方法能够很好地体现光源的稀疏性特点，但这会造成最优化问题目标函数不可微分，这就使得现有的很多用于L2范数正则化方法的光源重建方法不能直接用来求解基于LI范数的光源重建问题。除此之外，尽管现在已经存在一些可以被用来求解基于LI范数的正则化问题的数学方法，如迭代收缩方法，内点法等等，但是，自发荧光断层成像问题有病态性强，以及系统矩阵为实数稠密矩阵的特点，这使得以上的方法存在成像效率低，参数选择困难，甚至只能得到局部最优解而不能够得到准确的重建结果。
技术实现思路
本专利技术的目的是提供一种。为实现上述目的，一种，包括步骤SI使用有限元分析方法扩散方程进行离散化，基于LI范数的惩罚项建立无约束条件最优化问题模型；S2得到所述无约束条件最优化问题模型的对偶模型；S3建立所述对偶模型的增广拉格朗日函数；S4 :简化增广拉格朗日函数的最大值函数；S5使用截断牛顿法求解增广拉格朗日函数的最大值；S6将增广拉格朗日函数的梯度作为目标向量的最快下降方向对目标向量进行更新;S7更新惩罚向量；S8计算目标函数值J(W)，如果I I (JWk-J(W)H) I /I |Φ I≥tol为真，计算k=k+Ι并跳至步骤S4，否则，结束计算，其中，tol为目标函数的收敛效率阈值。本专利技术能够快速地在较大成像区域内得到准确可靠的光源分布信息。并且，本专利技术除正则化参数的以外其他参数都可以实现自适应调整，从而提高了成像的鲁棒性。附图说明 图1，本专利技术方法的流程图。图2，小鼠躯干的非匀质模型示意图以及对应的四面体网格。图3，非匀质模型体内光源位置示意图以及表面测量值。图4,两种方法的成像结果示意图(a) (b) (c)依次为Tikhonov方法(全域成像)、Tikhonov方法(局部成像)以及本专利技术所得到的成像结果。左侧一列为使用iso-surface的三维视图，右侧一列为在最大值点处进行切片得到的荧光光强分布图。图5，两个光源成像结果示意图左侧为使用iso-surface的三维视图，右侧为在最大值点处进行切片得到的荧光光强分布图。图6,三个光源成像结果示意图左侧为使用iso-surface的三维视图,右侧为在最大值点处进行切片得到的荧光光强分布图。具体实施例方式本专利技术所提出的方法充分考虑了自发荧光断层成像病态性强，系统矩阵为稠密矩阵的特点，在构建和求解拉格朗日乘子函数的过程中运用了对偶技术、自适应的迭代收缩策略、以及预算子共轭梯度法等技术。因此，本专利技术不但能够实现基于LI范数的光源重建，更重要的是降低了对参数选择的敏感度和提高了重建的效率。下面将结合附图详细描述本专利技术的重建方法，应指出的是，所描述的实施例仅旨在便于对本专利技术的理解，而对其不起任何限定作用。图1是本专利技术的重建方法的总体流程。下面对本专利技术的重建方法所涉及的关键步骤进行逐一详细说明，具体形式如下面所述步骤101 :使用有限元分析方法将扩散方程进行离散化，结合成像区域的四面体网格剖分，得到光源分布和边界测量数据之间的线性关系，然后使用基于LI范数的惩罚项建立起无约束条件最优化问题模型min J(w) = -^Aw - Φ I +A||w| 22其中A是系统矩阵，Φ"1为通过CCD相机所得到的被观测物体表面的光强值，λ是正则化参数，J(w)为最小二乘最优化问题的目标函数，w表示光源的分布信息，其初始值设定为零；尽管辐射传输方程能够对光子在生物组织中的传输进行精确地描述，但是该方程是一个包含多个变量的微分-积分方程，直接对该方程进行求解非常复杂。巨大的计算量限制了该方程在实际情况中的应用，所以在实际应用大都是使用该方程的简化和近似模型。由于在近红外和可见光波长范围内，大多数生物体组织都具有高散射体吸收的特性，即生物体对光子的散射作用远远大于吸收作用，所以可以使用扩散方程代替辐射传输方程来描述光子在生物体组织中的传输-V<D(x,2)VΦ(χ,λ)) + μα(χ,λ)Φ{χ,I) = S(x,2)(χ e Ω)其中Ω表示整个成像区域，Φ (r)表光流密度，S(r)表示光源强度，μ a(r)是组织吸收系数，μ ’s(r)表示约化组织反射系数，D(r) = 1/3(μ3ω + μ \(r))表示光学扩算系数。 在实际应用过程中，边界上光学信号的采集工作都是在暗箱中进行的，所以假定没有光子通过边界0Ω进入成像区域Ω，这符合使用Robin边界条件的要求Φ(χ, λ) + 2Α(χ; η, η,) (χ,λ)(\(χ}'^Φ(χ, 2))=0(χ e 5Ω)结合有限元分析方法，我们可以得到扩散方程及其边界条件的弱解形式f (D(x, λ)(νφ(χ, 2))(V￥(x, λ)) + μα (χ, Λ)Φ(χ, λ)Ψ(χ, λ)) χ J Ω+I ——-——Φ(χ λ)Ψ(χ,X)dx = f S(x，Ζ)Ψ(χ，Jl)dx其中，Ψ(χ，λ)为基函数。根据对成像区域进行四面体剖分得到的网格，结合对应的基函数，可以将上述的弱解形式离散化，得到下面的等式ΡΦ = FS其中，P是正定矩阵，Φ表示边界结点上光强，F表示光源加权矩阵，S表示成像区域内的光强分布。由于Φ中包含不可测量结点，所以去除掉F1F中与之对应的行向量，并得到AS = Φ111其中A e Rmxn表示系统矩阵，Φ111表示测量值。由于荧光断层成像是一个有严重病态性的逆问题，现有技术都是需要融入一个惩罚项函数作为正则化项使得求解过程变得更加稳定，本专利技术方法采用基于LI范数的惩罚项将光源的稀疏性特点作为先验知识融入到光源的重建中，使得荧光断层成像重建问题转化为具有如下形式的优化问题模型mill J(w) = —||,4w—Φ || +2||w|其中本文档来自技高网...

【技术保护点】
一种基于乘子法的自发荧光断层成像重建方法，包括步骤：S1使用有限元分析方法扩散方程进行离散化，基于L1范数的惩罚项建立无约束条件最优化问题模型；S2得到所述无约束条件最优化问题模型的对偶模型；S3建立所述对偶模型的增广拉格朗日函数；S4：简化增广拉格朗日函数的最大值函数；S5使用截断牛顿法求解增广拉格朗日函数的最大值；S6将增广拉格朗日函数的梯度作为目标向量的最快下降方向对目标向量进行更新；S7更新惩罚向量；S8计算目标函数值J(w)，如果||(J(w)k?J(w)k?1)||/||Φm||≥tol为真，计算k＝k+1并跳至步骤S4，否则，结束计算，其中，tol为目标函数的收敛效率阈值。

【技术特征摘要】

【专利技术属性】
技术研发人员：田捷，郭伟，杨鑫，
申请(专利权)人：中国科学院自动化研究所，
类型：发明
国别省市：