專利名稱：基于乘子法的自發(fā)熒光斷層成像重建方法
技術領域：
本發(fā)明屬于光學分子影像領域，涉及自發(fā)熒光斷層成像重建方法，尤其是一種基于乘子法的自發(fā)熒光斷層成像重建方法。
背景技術：
光學分子影像技術是在分子細胞水平上實現(xiàn)生物體生理、病理變化的實時、無創(chuàng)、在體成像。熒光斷層成像作為一種重要的光學分子影像成像模態(tài)，具有高靈敏度、低費用、無放射性等諸多優(yōu)點，已經(jīng)在新型藥物研發(fā)，在體觀測腫瘤發(fā)生、發(fā)展、轉移機理等多個領域獲得了廣泛的應用。由于在近紅外和可見光波長范圍里，光子在生物體 組織中傳播時會發(fā)生復雜的散射，因而二維平面自發(fā)熒光成像往往不能提供自發(fā)熒光光源的準確位置，而且具有無法獲得熒光光源深度信息的固有缺陷。自發(fā)熒光斷層成像重建技術通過測量邊界表面的突光光強分布，實現(xiàn)了突光光源的三維精確定位。自發(fā)熒光斷層成像重建是一個典型的病態(tài)問題。而且通常，表面熒光光強非常微弱，導致CCD探測器所采集到的信號混有大量干擾噪聲，從而進一步增加了熒光光源重建的復雜性，限制了自發(fā)熒光成像技術在實際中的應用。研發(fā)一種能夠準確、快速、魯棒的熒光光源重建技術一直是國際上一個研究熱點，但是到目前為止，還沒有一種公認的技術能夠實現(xiàn)這一目標?，F(xiàn)有的重建方法大都是將熒光光源重建看作是求解最優(yōu)化問題。正則化方法是一個很通用的方法。其中，最常見的方法是使用基于L2范數(shù)的正則化方法，但是這類方法得到的重建結果往往過于平滑，精度也不夠理想。考慮到在實際應用中，如腫瘤的早期檢測和觀察，熒光光源分布具有稀疏性分布的特點，因此，將稀疏性特點作為先驗信息融入到優(yōu)化問題模型中可以增強熒光光源重建的質量。因此基于LI范數(shù)的正則化方法能夠很好地體現(xiàn)光源的稀疏性特點，但這會造成最優(yōu)化問題目標函數(shù)不可微分，這就使得現(xiàn)有的很多用于L2范數(shù)正則化方法的光源重建方法不能直接用來求解基于LI范數(shù)的光源重建問題。除此之外，盡管現(xiàn)在已經(jīng)存在一些可以被用來求解基于LI范數(shù)的正則化問題的數(shù)學方法，如迭代收縮方法，內(nèi)點法等等，但是，自發(fā)熒光斷層成像問題有病態(tài)性強，以及系統(tǒng)矩陣為實數(shù)稠密矩陣的特點，這使得以上的方法存在成像效率低，參數(shù)選擇困難，甚至只能得到局部最優(yōu)解而不能夠得到準確的重建結果。

發(fā)明內(nèi)容
本發(fā)明的目的是提供一種基于乘子法的自發(fā)熒光斷層成像重建方法。為實現(xiàn)上述目的，一種基于乘子法的自發(fā)熒光斷層成像重建方法，包括步驟SI使用有限元分析方法擴散方程進行離散化，基于LI范數(shù)的懲罰項建立無約束條件最優(yōu)化問題模型；S2得到所述無約束條件最優(yōu)化問題模型的對偶模型；S3建立所述對偶模型的增廣拉格朗日函數(shù)；
S4 :簡化增廣拉格朗日函數(shù)的最大值函數(shù)；S5使用截斷牛頓法求解增廣拉格朗日函數(shù)的最大值；S6將增廣拉格朗日函數(shù)的梯度作為目標向量的最快下降方向對目標向量進行更新;S7更新懲罰向量；S8計算目標函數(shù)值J(W)，如果I I (JWk-J(W)H) I /I |Φ I≥tol為真，計算k=k+Ι并跳至步驟S4，否則，結束計算，其中，tol為目標函數(shù)的收斂效率閾值。本發(fā)明能夠快速地在較大成像區(qū)域內(nèi)得到準確可靠的光源分布信息。并且，本發(fā)明除正則化參數(shù)的以外其他參數(shù)都可以實現(xiàn)自適應調整，從而提高了成像的魯棒性。


圖1，本發(fā)明方法的流程圖。圖2，小鼠軀干的非勻質模型示意圖以及對應的四面體網(wǎng)格。圖3，非勻質模型體內(nèi)光源位置示意圖以及表面測量值。圖4,兩種方法的成像結果示意圖(a) (b) (c)依次為Tikhonov方法(全域成像)、Tikhonov方法(局部成像)以及本發(fā)明所得到的成像結果。左側一列為使用iso-surface的三維視圖，右側一列為在最大值點處進行切片得到的熒光光強分布圖。圖5，兩個光源成像結果示意圖左側為使用iso-surface的三維視圖，右側為在最大值點處進行切片得到的熒光光強分布圖。圖6,三個光源成像結果示意圖左側為使用iso-surface的三維視圖,右側為在最大值點處進行切片得到的熒光光強分布圖。
具體實施例方式本發(fā)明所提出的方法充分考慮了自發(fā)熒光斷層成像病態(tài)性強，系統(tǒng)矩陣為稠密矩陣的特點，在構建和求解拉格朗日乘子函數(shù)的過程中運用了對偶技術、自適應的迭代收縮策略、以及預算子共軛梯度法等技術。因此，本發(fā)明不但能夠實現(xiàn)基于LI范數(shù)的光源重建，更重要的是降低了對參數(shù)選擇的敏感度和提高了重建的效率。下面將結合附圖詳細描述本發(fā)明的重建方法，應指出的是，所描述的實施例僅旨在便于對本發(fā)明的理解，而對其不起任何限定作用。圖1是本發(fā)明的重建方法的總體流程。下面對本發(fā)明的重建方法所涉及的關鍵步驟進行逐一詳細說明，具體形式如下面所述步驟101 :使用有限元分析方法將擴散方程進行離散化，結合成像區(qū)域的四面體網(wǎng)格剖分，得到光源分布和邊界測量數(shù)據(jù)之間的線性關系，然后使用基于LI范數(shù)的懲罰項建立起無約束條件最優(yōu)化問題模型min J(w) = -^Aw - Φ I +A||w|
22其中A是系統(tǒng)矩陣，Φ"1為通過CCD相機所得到的被觀測物體表面的光強值，λ是正則化參數(shù)，J(w)為最小二乘最優(yōu)化問題的目標函數(shù)，w表示光源的分布信息，其初始值設定為零；盡管輻射傳輸方程能夠對光子在生物組織中的傳輸進行精確地描述，但是該方程是一個包含多個變量的微分-積分方程，直接對該方程進行求解非常復雜。巨大的計算量限制了該方程在實際情況中的應用，所以在實際應用大都是使用該方程的簡化和近似模型。由于在近紅外和可見光波長范圍內(nèi)，大多數(shù)生物體組織都具有高散射體吸收的特性，即生物體對光子的散射作用遠遠大于吸收作用，所以可以使用擴散方程代替輻射傳輸方程來描述光子在生物體組織中的傳輸-V<D(x,2)VΦ(χ,λ)) + μα(χ,λ)Φ{χ,I) = S(x,2)(χ e Ω)其中Ω表示整個成像區(qū)域，Φ (r)表光流密度，S(r)表示光源強度，μ a(r)是組織吸收系數(shù)，μ ’s(r)表示約化組織反射系數(shù)，D(r) = 1/3(μ3ω + μ \(r))表示光學擴算系數(shù)。 在實際應用過程中，邊界上光學信號的采集工作都是在暗箱中進行的，所以假定沒有光子通過邊界0Ω進入成像區(qū)域Ω，這符合使用Robin邊界條件的要求Φ(χ, λ) + 2Α(χ; η, η,) (χ,λ)(\(χ}'^Φ(χ, 2))=0(χ e 5Ω)結合有限元分析方法，我們可以得到擴散方程及其邊界條件的弱解形式f (D(x, λ)(νφ(χ, 2))(V￥(x, λ)) + μα (χ, Λ)Φ(χ, λ)Ψ(χ, λ)) χ
J Ω+I ——-——Φ(χ λ)Ψ(χ,X)dx = f S(x，Ζ)Ψ(χ，Jl)dx其中，Ψ(χ，λ)為基函數(shù)。根據(jù)對成像區(qū)域進行四面體剖分得到的網(wǎng)格，結合對應的基函數(shù)，可以將上述的弱解形式離散化，得到下面的等式ΡΦ = FS其中，P是正定矩陣，Φ表示邊界結點上光強，F(xiàn)表示光源加權矩陣，S表示成像區(qū)域內(nèi)的光強分布。由于Φ中包含不可測量結點，所以去除掉F1F中與之對應的行向量，并得到AS = Φ111其中A e Rmxn表示系統(tǒng)矩陣，Φ111表示測量值。由于熒光斷層成像是一個有嚴重病態(tài)性的逆問題，現(xiàn)有技術都是需要融入一個懲罰項函數(shù)作為正則化項使得求解過程變得更加穩(wěn)定，本發(fā)明方法采用基于LI范數(shù)的懲罰項將光源的稀疏性特點作為先驗知識融入到光源的重建中，使得熒光斷層成像重建問題轉化為具有如下形式的優(yōu)化問題模型mill J(w) = —||,4w—Φ || +2||w|其中，w表示光源的分布信息，λ表示正則化參數(shù)。步驟102 :使用對偶變換將步驟I中所述的無約束條件最優(yōu)化問題模型進行轉化，得到相應的對偶模型max is (α, V) = - 1| - Φ If + ^ || ra If - C (^) subject to v_AT a =0該對偶模型是一個約束條件最優(yōu)化問題模型，目標解向量a的維數(shù)要遠小于w的維數(shù)，這將會顯著的提高成像的效率。由于步驟101中所提出的最優(yōu)化問題目標函數(shù)明顯是不可微分的，這使得現(xiàn)有很多基于L2范數(shù)的熒光光源重建方法不能直接用來求解該模型。針對上述問題，本發(fā)明使用對偶變換，將步驟101中的非約束條件最優(yōu)化問題模型進行轉化，并得到與之等價的約束條件最優(yōu)化問題模型。原模型中的目標變量w和對偶模型中的目標變量α的維數(shù)分別是M和N，M和N分別對應于系統(tǒng)矩陣的維數(shù)，還分別對應了表面熒光可測量結點的維數(shù)和成像區(qū)域所有結點的維數(shù)。顯然，運用對偶變換將明顯降低模型目標變量的維數(shù)。同時，轉化所得到的約束條件最優(yōu)化問題模型也為下一步使用增廣拉格朗日函數(shù)技術奠定了基礎。步驟103 :使用增廣拉格朗日函數(shù)函數(shù)技術，建立步驟2所述對偶模型的增廣拉格朗日函數(shù)L(a, v. w, μ) = E(a, r) - wT(ATa — v) — γ||^Γα - 其中μ為懲罰向量,w為拉格朗日乘子。在步驟102所得出的對偶最優(yōu)化目標函數(shù)maxEU，v)的基礎上添加增廣拉格朗日懲罰項-就可以得到模型的增廣拉格朗日函數(shù)；當μ =0時，得到的是模型的拉格朗日函數(shù)。此時，表示熒光光源分布信息的W成為了拉格朗日乘子。步驟104 :將V表示為α的削頂函數(shù)= CM；； + #-’增廣拉格朗日函數(shù)的最
H.
大值函數(shù)可以化簡為
2
I Il112 TfwmaxL(a,w,μ) = max--p-Φ η--ShrinkΛ(Α a-\——)
211丨丨2 2η 2從而降低了模型求解的復雜度。將步驟104中的增廣拉格朗日乘子函數(shù)看作是V的函數(shù)并展開為向量的形式，由于最大化函數(shù)中的每一項都能夠被展開為N維，并且都包含Vj，由此可以得到
JiTL(a, v, W，μ) = --||α — Φ 『+-|| m||2 —乞(令(vf — (-+^4 )2 + S^vj))
Z2,/=i LJU其中Vj表示向量的第j個元素。結合指示函數(shù)€的性質，vU)的最大值可以用包含α的削頂函數(shù)1= + 表示，并得到新的增廣拉格朗日函數(shù)
β
I2H2max L(a, w, μ) = max-—1| - Φ一.Shrink^ (Ara + —)從而降低了模型求解的復雜度。其中，ShirnkA (w)是軟閾值函數(shù)，其定義為
I IW,shirnkA(w) = (max(|%|-2,0)|^|), (, = 1，...，《)，削頂函數(shù)的定義為CLa O’) := (min |> I，λ)步驟105 :使用融合積極集策略的截斷牛頓法求解增廣拉格朗日函數(shù)的最大值maxLa ( a , w, u )。本發(fā)明的特性決定了重建熒光光源所需的計算量主要集中在求解步驟104所述 增廣拉格朗日函數(shù)的最大值，所以在本步驟中，通過綜合使用牛頓法和積極集策略提高最 大值的求解效率。下面對本步驟所涉及的子步驟進行詳細說明步驟201 :計算積極集A+，該矩陣是系統(tǒng)矩陣A的子矩陣，由
權利要求
1.一種基于乘子法的自發(fā)熒光斷層成像重建方法，包括步驟 SI使用有限元分析方法擴散方程進行離散化，基于LI范數(shù)的懲罰項建立無約束條件最優(yōu)化問題模型； S2得到所述無約束條件最優(yōu)化問題模型的對偶模型； S3建立所述對偶模型的增廣拉格朗日函數(shù)； S4 :簡化增廣拉格朗日函數(shù)的最大值函數(shù)； S5使用截斷牛頓法求解增廣拉格朗日函數(shù)的最大值； S6將增廣拉格朗日函數(shù)的梯度作為目標向量的最快下降方向對目標向量進行更新； S7更新懲罰向量； S8計算目標函數(shù)值J(W),如果
2.根據(jù)權利要求1所述的方法，其特征在于所述無約束條件最優(yōu)化問題模型由下式表達
3.根據(jù)權利要求2所述的方法，其特征在于所述對偶模型由下式表達
4.根據(jù)權利要求1所述的方法，其特征在于所述增廣拉格朗日函數(shù)由下式表達
5.根據(jù)權利要求4所述的方法，其特征在于通過削頂函數(shù)技術將ν表示成α的函數(shù)
6.根據(jù)權利要求1所述的方法，其特征在于所述使用截斷牛頓法求解增廣拉格朗日函數(shù)的最大值包括 計算積極集A+ ； 集合積極集計算出增廣拉格朗日函數(shù)的海森矩陣和雅克比矩陣 通過使用預共軛梯度法計算
7.根據(jù)權利要求6所述的方法，其特征在于所述雅克比矩陣和海森矩陣由下式計算得到Y2aLia) = Γ+μ, A+A ▼鄭、=α-Φπ! +ΑΕΤλ (wk + μ(Α α) 其中ShirnkA (w)為標準軟閾值函數(shù)，其表達式為Shirnki (w) = (max(|w 1-2,0) τ-^τ) (/ = 1，…，w) '1 Id 。
8.根據(jù)權利要求7所述的方法，其特征在于選取海森矩陣的對角元素得到預共軛梯度算子P: P = diagiV1!⑷)=diag(—Im —凡4^) O
9.根據(jù)權利要求6所述的方法，其特征在于積極集A+是系統(tǒng)矩陣A的子矩陣，包含了由 J+ = Ue {1,2, , η} I Xk+ μ kAT α > λ μ |J 所得到的積極集。
10.根據(jù)權利要求1所述的方法，其特征在于使用軟閾值函數(shù)更新目標向量。
全文摘要
一種基于乘子法的自發(fā)熒光斷層成像重建方法，使用有限元分析方法擴散方程進行離散化，基于L1范數(shù)的懲罰項建立無約束條件最優(yōu)化問題模型；得到所述無約束條件最優(yōu)化問題模型的對偶模型；建立所述對偶模型的增廣拉格朗日函數(shù)；簡化增廣拉格朗日函數(shù)的最大值函數(shù)；使用截斷牛頓法求解增廣拉格朗日函數(shù)的最大值；將增廣拉格朗日函數(shù)的梯度作為目標向量的最快下降方向對目標向量進行更新；更新懲罰向量；計算目標函數(shù)值J(w)，如果||(J(w)k-J(w)k-1)||/||Φm||≥tol為真，計算k＝k+1并跳至步驟S4，否則，結束計算，其中，tol為目標函數(shù)的收斂效率閾值。本發(fā)明能夠快速地在較大成像區(qū)域內(nèi)得到準確可靠的光源分布信息，除正則化參數(shù)的以外其他參數(shù)都可以實現(xiàn)自適應調整，從而提高了成像的魯棒性。
文檔編號A61B5/00GK102988026SQ20121052358
公開日2013年3月27日 申請日期2012年12月7日 優(yōu)先權日2012年12月7日
發(fā)明者田捷, 郭偉, 楊鑫 申請人:中國科學院自動化研究所