专利名称：基于乘子法的自发荧光断层成像重建方法
技术领域：
本发明属于光学分子影像领域，涉及自发荧光断层成像重建方法，尤其是一种基于乘子法的自发荧光断层成像重建方法。
背景技术：
光学分子影像技术是在分子细胞水平上实现生物体生理、病理变化的实时、无创、在体成像。荧光断层成像作为一种重要的光学分子影像成像模态，具有高灵敏度、低费用、无放射性等诸多优点，已经在新型药物研发，在体观测肿瘤发生、发展、转移机理等多个领域获得了广泛的应用。由于在近红外和可见光波长范围里，光子在生物体 组织中传播时会发生复杂的散射，因而二维平面自发荧光成像往往不能提供自发荧光光源的准确位置，而且具有无法获得荧光光源深度信息的固有缺陷。自发荧光断层成像重建技术通过测量边界表面的突光光强分布，实现了突光光源的三维精确定位。自发荧光断层成像重建是一个典型的病态问题。而且通常，表面荧光光强非常微弱，导致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|其中，w表示光源的分布信息，λ表示正则化参数。步骤102 :使用对偶变换将步骤I中所述的无约束条件最优化问题模型进行转化，得到相应的对偶模型max is (α, V) = - 1| - Φ If + ^ || ra If - C (^) subject to v_AT a =0该对偶模型是一个约束条件最优化问题模型，目标解向量a的维数要远小于w的维数，这将会显著的提高成像的效率。由于步骤101中所提出的最优化问题目标函数明显是不可微分的，这使得现有很多基于L2范数的荧光光源重建方法不能直接用来求解该模型。针对上述问题，本发明使用对偶变换，将步骤101中的非约束条件最优化问题模型进行转化，并得到与之等价的约束条件最优化问题模型。原模型中的目标变量w和对偶模型中的目标变量α的维数分别是M和N，M和N分别对应于系统矩阵的维数，还分别对应了表面荧光可测量结点的维数和成像区域所有结点的维数。显然，运用对偶变换将明显降低模型目标变量的维数。同时，转化所得到的约束条件最优化问题模型也为下一步使用增广拉格朗日函数技术奠定了基础。步骤103 :使用增广拉格朗日函数函数技术，建立步骤2所述对偶模型的增广拉格朗日函数L(a, v. w, μ) = E(a, r) - wT(ATa — v) — γ||^Γα - 其中μ为惩罚向量,w为拉格朗日乘子。在步骤102所得出的对偶最优化目标函数maxEU，v)的基础上添加增广拉格朗日惩罚项-就可以得到模型的增广拉格朗日函数；当μ =0时，得到的是模型的拉格朗日函数。此时，表示荧光光源分布信息的W成为了拉格朗日乘子。步骤104 :将V表示为α的削顶函数= CM；； + #-’增广拉格朗日函数的最
H.
大值函数可以化简为
2
I Il112 TfwmaxL(a,w,μ) = max--p-Φ η--ShrinkΛ(Α a-\——)
211丨丨2 2η 2从而降低了模型求解的复杂度。将步骤104中的增广拉格朗日乘子函数看作是V的函数并展开为向量的形式，由于最大化函数中的每一项都能够被展开为N维，并且都包含Vj，由此可以得到
JiTL(a, v, W，μ) = --||α — Φ 『+-|| m||2 —乞(令(vf — (-+^4 )2 + S^vj))
Z2,/=i LJU其中Vj表示向量的第j个元素。结合指示函数€的性质，vU)的最大值可以用包含α的削顶函数1= + 表示，并得到新的增广拉格朗日函数
β
I2H2max L(a, w, μ) = max-—1| - Φ一.Shrink^ (Ara + —)从而降低了模型求解的复杂度。其中，ShirnkA (w)是软阈值函数，其定义为
I IW,shirnkA(w) = (max(|%|-2,0)|^|), (, = 1，...，《)，削顶函数的定义为CLa O’) := (min |> I，λ)步骤105 :使用融合积极集策略的截断牛顿法求解增广拉格朗日函数的最大值maxLa ( a , w, u )。本发明的特性决定了重建荧光光源所需的计算量主要集中在求解步骤104所述 增广拉格朗日函数的最大值，所以在本步骤中，通过综合使用牛顿法和积极集策略提高最 大值的求解效率。下面对本步骤所涉及的子步骤进行详细说明步骤201 :计算积极集A+，该矩阵是系统矩阵A的子矩阵，由
权利要求
1.一种基于乘子法的自发荧光断层成像重建方法，包括步骤 SI使用有限元分析方法扩散方程进行离散化，基于LI范数的惩罚项建立无约束条件最优化问题模型； S2得到所述无约束条件最优化问题模型的对偶模型； S3建立所述对偶模型的增广拉格朗日函数； S4 :简化增广拉格朗日函数的最大值函数； S5使用截断牛顿法求解增广拉格朗日函数的最大值； S6将增广拉格朗日函数的梯度作为目标向量的最快下降方向对目标向量进行更新； S7更新惩罚向量； S8计算目标函数值J(W),如果
2.根据权利要求1所述的方法，其特征在于所述无约束条件最优化问题模型由下式表达
3.根据权利要求2所述的方法，其特征在于所述对偶模型由下式表达
4.根据权利要求1所述的方法，其特征在于所述增广拉格朗日函数由下式表达
5.根据权利要求4所述的方法，其特征在于通过削顶函数技术将ν表示成α的函数
6.根据权利要求1所述的方法，其特征在于所述使用截断牛顿法求解增广拉格朗日函数的最大值包括 计算积极集A+ ； 集合积极集计算出增广拉格朗日函数的海森矩阵和雅克比矩阵 通过使用预共轭梯度法计算
7.根据权利要求6所述的方法，其特征在于所述雅克比矩阵和海森矩阵由下式计算得到Y2aLia) = Γ+μ, A+A ▼鄭、=α-Φπ! +ΑΕΤλ (wk + μ(Α α) 其中ShirnkA (w)为标准软阈值函数，其表达式为Shirnki (w) = (max(|w 1-2,0) τ-^τ) (/ = 1，…，w) '1 Id 。
8.根据权利要求7所述的方法，其特征在于选取海森矩阵的对角元素得到预共轭梯度算子P: P = diagiV1!⑷)=diag(—Im —凡4^) O
9.根据权利要求6所述的方法，其特征在于积极集A+是系统矩阵A的子矩阵，包含了由 J+ = Ue {1,2, , η} I Xk+ μ kAT α > λ μ |J 所得到的积极集。
10.根据权利要求1所述的方法，其特征在于使用软阈值函数更新目标向量。
全文摘要
一种基于乘子法的自发荧光断层成像重建方法，使用有限元分析方法扩散方程进行离散化，基于L1范数的惩罚项建立无约束条件最优化问题模型；得到所述无约束条件最优化问题模型的对偶模型；建立所述对偶模型的增广拉格朗日函数；简化增广拉格朗日函数的最大值函数；使用截断牛顿法求解增广拉格朗日函数的最大值；将增广拉格朗日函数的梯度作为目标向量的最快下降方向对目标向量进行更新；更新惩罚向量；计算目标函数值J(w)，如果||(J(w)k-J(w)k-1)||/||Φm||≥tol为真，计算k＝k+1并跳至步骤S4，否则，结束计算，其中，tol为目标函数的收敛效率阈值。本发明能够快速地在较大成像区域内得到准确可靠的光源分布信息，除正则化参数的以外其他参数都可以实现自适应调整，从而提高了成像的鲁棒性。
文档编号A61B5/00GK102988026SQ20121052358
公开日2013年3月27日 申请日期2012年12月7日 优先权日2012年12月7日
发明者田捷, 郭伟, 杨鑫 申请人:中国科学院自动化研究所