基于l1范數與tv范數的復合正則化生物發光斷層成像重建方法
【專利摘要】本發明屬于醫學圖像處理領域，涉及一種基于L1范數與TV范數的復合正則化生物發光斷層成像重建方法。本發明聯合L1正則化和TV正則化的復合正則化方法，融合兩種正則化方法的優點，并突破單一正則化方法的局限性，進而提高了成像的質量，并通過基于回溯線搜索(Backtracking Line Search，BLS)的共軛梯度下降(Conjugate Gradient，CG)方法進行求解。實驗結果表明，本方法不但可以對熒光光源進行準確重建定位，同時有極高的計算效率。
【專利說明】
基于L1范數與TV范數的復合正則化生物發光斷層成像重建 方法
技術領域
[0001] 本發明屬于醫學圖像處理領域，涉及一種基于L1范數與TV范數的復合正則化生物 發光斷層成像重建方法。
【背景技術】
[0002] 隨著人類對生命奧秘的不斷探索以及現代計算機技術的飛速發展，各種醫學影像 技術和醫學成像設備也進入了快速發展期。從傳統的計算機斷層成像（Computed Tomography，CT)、超生成像（Ultrasound Imaging，USI)、核磁共振成像（Magnetic Resonance Imaging，MRI)、正電子成像(Positron Emission Tomography，PET)、單光子發 射斷層成像（Single-Photon Emission Computed Tomography，SPECT)發展到激發焚光斷 層成像（Fluorescent Molecular Tomography，FMT)和生物發光斷層成像 (Bioluminescence Tomography，BLT)等多種成像技術，這些成像技術的發展也極大地推動 了生命科學的進步。
[0003] 生物發光斷層成像（以下簡稱BLT)作為一種重要的成像技術，近年來已經成為光 學分子影像的一個重要分支。BLT憑借其高靈敏度、背景噪聲低等特性，通過在生物體表獲 得熒光信號從而重建生物體內熒光源的分布以及觀測細胞分子水平變化。BLT技術不需要 外在光源的激發，而是通過一種生物化學發光反應在體內發光。體內產生的熒光在生物組 織內部以某種規律傳播并不斷地與生物組織發生相互作用，并到達體表。最后，在生物組織 體表利用高靈敏度的探測器獲得的熒光圖像就可以重建出熒光光源在小動物體內的分布 情況，從而在本質上揭示在體分子的活動規律。
[0004] 由于生物組織具有強散射、弱吸收的特點，熒光光子無法在生物組織中沿直線傳 輸，而是經歷了大量的散射過程，導致BLT逆問題在數學上是一個高度病態的問題，外界微 小的測量擾動都會給重建結果帶來很大的變化。因此，降低BLT的病態性，如何唯一、準確地 重建熒光光源是BLT研究的重點及熱點。
[0005] 近年來，國內外研究學者為了提高BLT成像的準確性和可靠性做了各種嘗試。為了 降低BLT重建問題的病態性，在重建時可以使用正則化的求解方法將光源重建問題轉變成 一個非線性的最優化問題。正則化方法最早在上世紀60年代由蘇聯科學院院士 Tikhonov提 出，用來求解形如Fx = y方程的不適定問題。正則化的核心思想是在求解問題過程中引入先 驗信息，從而得到對原問題有意義的解。正則化函數構造主要針對正則化項的形式進行研 究，目前有L2范數、總變分(total variation，TV)范數和LP范數的正則化項形式。基于L2范 數的正則化項容易求解，但重建的圖像會過于平滑。基于稀疏特性的正則化方法(典型的的 是L1正則化）由于融入了光源的稀疏先驗信息，能夠提高重建圖像的質量。但是該方法會導 致重建的光源過于稀疏，進而降低成像的質量。TV正則化方法由于突出了光源的邊界信息 從而改善了成像質量，但該方法會消除重建圖像中一些小的特征信息和小的物體信息（如 稀疏光源）。因此，從上面的分析中可以看出，任何一種正則化方法都有自身固有的優勢與 缺陷。可以使用聯合L1正則化和TV正則化的復合正則化方法，融合兩種正則化方法的優點， 并突破單一正則化方法的局限性，進而提高了成像的質量。
[0006] 本發明提出一種基于L1范數與TV范數的復合正則化生物發光斷層成像重建方法， 并通過基于回溯線搜索（Backtracking Line Search，BLS)的共輒梯度下降（Conjugate Grad i ent，CG)方法進行求解。
【發明內容】
[0007] 在BLT的研究中，通常認為其工作在連續波的工作模式下。其在求解的過程中通常 不需要考慮時間的影響，認為生物發光光源長時間內能量穩定，因此可以選取穩態漫射方 程對光子在生物組織中的傳輸過程進行描述：
[0008] -▽.「￡)▽0(，)] + ""(6)(/.) = 5(/.) (1}
[0009] 其中r表示求解域中的一個點，Φ (r)表示r點處的光子通量流率，S(r)表示組織內 部光源的能量密度，D表示組織的散射系數，V與V，分別為梯度與散度的標識。考慮到生物 組織內部的非勻質性，也就是組織的光學特性參數在不同區域會有所區別，上述的漫射方 程可以改寫為如下形式：
[0010] -￥ · [ /) 0- ) νφ (r)] + βα (r ) Φ 0- ) = S (r ) (2)
[0011]其中D ( r )具體表示為r點處的擴散系數，其滿足
//" (r)分別表示生物組織的吸收系數與約化散射系 數。
[0012]在BLT實驗中，檢測到的是生物組織表面的光子流密度Q(r)，其可以表示為：
(3)
[0014] 其中，v(r)代表r位置處的外法線方向矢量，A(r)表示r位置處的邊界失配因子，在 后文會有詳細介紹。
[0015] 式(2)是一個經典的橢圓偏微分方程，根據微分方程的求解理論，對漫射方程的求 解還需要引入適當的邊界條件。常見的邊界條件包括外延邊界條件與Robin邊界條件兩種， 其中在光學成像中更常用到的是Robin邊界條件，其可以表示為如下形式：
[0016] 0(r) + 2,4(r)D(r)[VO(r)· v(r) | = 0 (4)
[0017] 在這里需要對邊界失配因子A(r)進行說明，其表示如下：
(5)
[0019] 其中R(r)代表的是擴散傳輸內反射系數，其值可由下式表示：
[0020] R(r) ^-1.4399n_2+0.7099n_1+0.6881+0.0636η (6)
[0021 ] η表示生物組織對空氣的相對折射率。
[0022]通常采用有限元的方法對漫射方程進行求解，其求解過程可以分為四個步驟:建 立積分方程、網格剖分建立有限元子空間、建立有限元方程和有限元方程的求解。
[0023]依據上述步驟，穩態漫射方程的求解可以如下表示：
[0024]在區域Ω內，假設光源S(r) = So(r)。利用分段連續的測試函數Ψ(Γ)乘以穩態漫 射方程的等式兩端，并在整個區域Ω內對等式兩端進行積分。Ψ(Γ)選自一個測試空間，它 及其所有一階導數都在區域Ω上可積，其邊界條件應與Φ相同。結合Robin邊界條件，使用 繼表示區域的邊界，采用Galerkin變分方法，可以得到如下的方程弱解形式：
[0026]對漫射方程的求解就可以轉化為尋找對任意Ψ(Γ)都可以滿足式(7)的Φ值。
[0027]將區域Ω離散為Ζ個頂點和Ν個單元，每個子區域可以表示為Ω(1)(1 = 1，2,…，Ν)。 此時，選擇私作為單元Ω(1)的節點基函數，Φ(Γ)即可以近似地表示為如下的分段多項式函 數：
L〇〇29J 在上式中，Φη(Γ)代表Φ(Γ)的有限元近似，在第m個節點'上對應的 值，仏代表區域Ω(1)的公共節點。將式(8)帶入到式(7)中，得到如下形式的弱解式：
[0031] 同理，光源項也可近似離散為如下形式：
[0032] 5(r) ? S" (r) = (r), re Ω (l〇) m=l
[0033] 在上式中，Sh(r)是光源S(r)的有限元近似，sm是第m個節點上光源的離散值， (r)為插值基函數。
[0034]對離散單元進行積分得到其有限元方程，合并后可得到如下矩陣形式的方程：
[0035] (K+C+B)?=M?=FS (11)
[0036] 其中，Φ為各節點處光子通量流率合并得來的矩陣，S為各節點處光源值合并得來 的矩陣，矩陣M=K+C+B，矩陣K、C、B、F的i行j列元素 k小叫、b小分別表示為如下形式的積 分：
[0038]如上文所述，釔(>·)、％(/-)分別為第i、j個有限元單元的節點基函數，γ Kr)為第j 個有限元單元插值基函數。
[0039]可以將式（11)進行整理得到如下的關系式：
[0040] (6=1^=^8 (13)
[0041] 在上式中，f =M4F，其說明了光子通量流率Φ與內部的光源能量密度之間滿足線 性關系。在BLT光源重建的過程中是根據得到的表面光子通量流率去尋找光源的分布。在實 際的實驗中，只有部分邊界光子能夠被捕獲，因此可以將Φ分為兩個部分，即：可測量的表 面節點上的光子通量流率Φ和不能測量的Φ'于是，在式（13)的矩陣方程中，去除系數矩 陣#中與Φ+相關的行得到系數矩陣A，可以得到如下的線性矩陣方程：
[0042] (14)
[0043]考慮到光的強散射特性與噪聲的影響會導致矩陣Α具有嚴重的病態性，同時由于 測量數據的不足，導致式（14)是一個病態的欠定方程。為了減弱方程的病態性，考慮到生物 發光斷層成像中發光光源的稀疏性，BLT重建問題可以使用正則化方法求解。
[0044] 基于正則化思想，可以將光源求解即BLT重建問題轉換為如下形式的最小二乘問 題：
[0045] minv f(S) = S", φ) +) (15)
[0046] 在上式中，Φ (S，Φ )代表數據擬合項，λ為正則化參數，以幻為正則化項，求解光源 S應使目標函數f (S)取最小值。
[0047] 除了常見的L1范數，本方法中還要用到TV范數，在有限元框架下，TV范數項| |S| |TV 可以表示成下面的形式： Μ
[_ II確= Σ，-產1 (16) p=l
[0049] 在上式中，N表示有限元節點數。在二維空間中，LP表示第p條邊的邊長，sp，lef^s p ，nght分別表示第P條邊的左、右端點;在三維空間中，LP表示第P個面的面積，SP， left和SP，nght 分別表示第P個面的左、右兩個體網格。
[0050] 針對生物發光斷層成像的有限元求解方式，本方法提出了一種基于L1范數與TV范 數的復合正則化方法對生物發光斷層成像重建問題求解，其表達式如下：
[0051 ] mi% f(S) = -Φ|[" + \ + Λ ||S||JF (Π )
[0052]在上式中，為求解光源分布矩陣S應使目標函數f(S)取最小值。在上式中， μ5-Φ|22為數據擬合項，用于表征求解值與測量值間的差異。hi |s| ^與~！ |s| |τν分別為L1 正則化與TV(總變分)正則化項，正則化項的引入可以降低方程的病態性。其中λ^λ2分別表 示L1范數正則化項與TV范數正則化項的正則化參數，一般情況下正則化參數不宜選取過 大，在本方法中λι與λ2的取值選為10- 3與nr4』為用于表不光源分布矩陣S與Φ間關系的系 數矩陣，Φ為測量得到的生物組織表面光子通量流率值。系數矩陣Α的計算，φ的測量及生 物組織的有限元劃分可以借助Matlab的nirfast工具箱完成。
[0053] 為了對式（17)進行求解，本方法采用了基于回溯線搜索（Backtracking Line Search，BLS)的共輒梯度下降(Conjugate Gradient，CG)迭代方法，其描述如下：
[0054] 1、初始化:設置初始迭代次數k = 0;設置回溯線搜索的初始步長to=l;設置光源 初始分布矩陣So,其矩陣元素均設置為0到1之間的值，光源初始分布矩陣初始值的設置會 影響迭代步數，但對最后的計算結果沒有影響;初始下降梯度& = ，下降梯度gk計算 公式為V/(5；) = 2, p+ 其中，V為求梯度的標識;gk代表第k 次迭代時的下降梯度;AT代表系數矩陣A的轉置;φ為測量得到的生物組織表面光子通量流 率值;Sk表示第k次迭代后解得的光源分布矩陣;初始改變值△ So設定為求得的初始下降梯 度求負，即Δ SQ = -gQ;通過輸入生物組織的結構模型及生物組織各區域的光學參數(吸收系 數、散射系數及折射率)到Matlab的nirfast工具箱計算得系數矩陣A，并通過nirfast工具 箱的前向仿真功能得到邊界測量值矩陣
[0055] 2、進行迭代的條件判斷：判斷下降梯度Iklg是否已經小于初始設定的閾值，迭代 次數k是否已經超過初始設定的最大迭代次數，滿足上述條件之一則結束迭代，并返回最終 的計算結果Sk;為了能夠減小誤差，初始設定的閾值一般應小于10-2();迭代次數不宜設置過 小，一般情況下應不小于1〇〇〇次，否則較為容易產生到達最大迭代次數但計算結果誤差較 大的情況。
[0056] 3、基于回溯線搜索的梯度下降迭代：如果滿足汽51{+丨1^51〇>汽5 1〇+€[^化1/八 Sk)，則改變搜索步長tk+i = Ptk，若不滿足則步長保持不變即tk+i = tk，其中tk與tk+i分別表示 第k與k+Ι次迭代計算出的步長;其中△ Sk為第k次迭代時計算出的光源分布矩陣改變值;α ￡(0，1)、0^(0，1)是線搜索的相關參數，<1的值一般取小于0.14的值一般取大于0.5為宜， 此時往往能夠獲得較好的迭代下降速度。
[0057] 4、迭代步驟:更新光源分布矩陣S的值Sk+1 = Sk+tk+1 Δ Sk，Sk+1為本次迭代計算得到 的光源分布矩陣，sk為上一次迭代計算得到的光源分布矩陣;更新下降梯度gi+1 =ν/ρ?+1); 計算
，其中gk+1為新計算得到的下降梯度，gk為上一次迭代時計算得到的下降梯 度;使用公式Δ Sk+i = -gk+i+ γ Δ Sk對改變值進行更新;更新迭代次數k = k+l并返回步驟2。
[0058] 使用傳統的梯度下降法往往速度較慢，通過使用上述的方法能夠大大提高方法的 計算效率。
【附圖說明】
[0059] 圖1為仿體初始的光源分布，圖中面積較大的黑色圓形區域為仿體區域，黑色區域 中面積較小的兩個白色區域為光源分布區域；
[0060] 圖2為通過本方法重建的光源分布結果。
【具體實施方式】
[0061] 下面根據具體實施示例與附圖對本發明進行說明。
[0062]首先，通過Matlab的nirfast建立仿體。本實驗使用圓形區域仿體，其包含兩個光 源區域，如圖1所示。
[0063] 為了避免逆行為（inverse crime)，前向仿真的有限元網格和重建的網格節點數 不同。一般情況下，為了保證重建時計算量不會過大，重建模型的有限元節點數往往小于前 向的有限元節點數。在本實驗中，前向網格的有限元節點數設置為3508個，共包含6807個面 元;重建網格的有限元節點數設置為1309個，共包含2491個面元。
[0064]由于BLT逆問題求解中基于單光譜的重建結果可能存在不唯一性，因此實驗中采 用了兩個不同的波段。不同光波段在仿體中擁有不同的光吸收系數與光散射系數，實驗選 擇的兩個譜段分別為600nm和630nm，此時仿體中的光學特性參數如表1所示。
[0065]表1仿體的光學特性參數
[0067] 分別對不同光譜下添加光源的仿體區域進行前向仿真，合并后可以得到邊界測量 值悉，并可以計算得到重建需要的系數矩陣A。
[0068] 使用式（17)所示的復合正則化方法進行光源重建，在實驗中設置個節點S的初始 值均為l〇_5，Ll范數正則化項與TV范數正則化項的正則化系數心與\ 2分別設置為10-3與ΚΓ4， 實驗證明在取該正則化系數值時能得到較好的重建效果。
[0069] 本發明使用基于回溯線搜索的共輒梯度下降迭代方法。為了保證重建結果的誤差 較小，實驗中設置最大迭代次數為20000,下降梯度的容忍閾值設為10- 3()。為了能夠得到較 好的梯度下降速度，線性搜素的相關參數設置如下：初始步長設為1，α與β的值分別設為 0.01與0.6。通過本方法計算，可以得出如圖2所示的光源重建結果。實驗結果表明，本方法 不但可以對熒光光源進行準確重建定位，同時有極高的計算效率。
【主權項】
1.基于LI范數與TV范數的復合正則化生物發光斷層成像重建方法，其特征在于:對生 物發光斷層成像重建問題求解，其表達式如下：(1) 在上式中，為求解光源分布矩陣S應使目標函數f(s)取最小值;在上式中，ιμ^-壺為數 據擬合項，用于表征求解值與測量值間的差異;λι II SII1與λ21 I SII TV分別為L1正則化與TV總 變分正則化項λι與λ2的取值選為10-3與ιοΛα為用于表示光源分布矩陣S與壺間關系的系數 矩陣，樂為測量得到的生物組織表面光子通量流率值； 采用W下步驟對式(1)求解； 1) 初始化:設置初始迭代次數k=0;設置回溯線捜索的初始步長t〇=l;設置光源初始分 布矩陣So,其矩陣元素均設置為0到1之間的值，初始下降梯度沉=^"5。)，下降梯度邑式十算 公式為巧(展寬)+AV|悶L+4V|悶Iw;其中，y為求梯度的標識;gk代表第k 次迭代時的下降梯度;AT代表系數矩陣A的轉置；壺為測量得到的生物組織表面光子通量流 率值;Sk表示第k次迭代后解得的光源分布矩陣;初始改變值Δ So設定為求得的初始下降梯 度求負，即ASo = -go;通過輸入生物組織的結構模型及生物組織各區域的光學參數到 Matlab的ni計ast工具箱計算得系數矩陣A，生物組織各區域的光學參數包括吸收系數、散 射系數及折射率，并通過nirfast工具箱的前向仿真功能得到邊界測量值矩陣壺； 2) 進行迭代的條件判斷:判斷下降梯度|&||^是否已經小于初始設定的闊值，迭代次數1^ 是否已經超過初始設定的最大迭代次數，滿足上述條件之一則結束迭代，并返回最終的計 算結果Sk;初始設定的闊值小于迭代次數不小于1000次，否則較為容易產生到達最大 迭代次數但計算結果誤差較大的情況； 3) 基于回溯線捜索的梯度下降迭代:如果滿足f(Sk+tkΔSk)>f(Sk)+αtk·(gkTASk)，則 改變捜索步長tk+l =化k，若不滿足則步長保持不變即tk+l = tk，其中tk與tk+汾別表示第k與k + 1次迭代計算出的步長;其中A Sk為第k次迭代時計算出的光源分布矩陣改變值;ae(〇， 1)、βε(〇，1)是線捜索的相關參數； 4) 迭代:更新光源分布矩陣S的值Sk+l = Skパk+lΔSk，Sk+l為本次迭代計算得到的光源分 布矩陣，Sk為上一次迭代計算得到的光源分布矩陣；更新下降梯度涼巧+1);計算，其中gk+i為新計算得到的下降梯度，阱為上一次迭代時計算得到的下降梯度;使 用公式Δ Sk+i = -gk+i+ 丫 Δ Sk對改變值進行更新;更新迭代次數k = k+l并返回步驟2)。
【文檔編號】A61B5/00GK106097441SQ201610475428
【公開日】2016年11月9日
【申請日】2016年6月25日
【發明人】馮金超, 李祎楠, 賈克斌
【申請人】北京工業大學