基于壓縮采樣匹配追蹤的生物發光斷層成像

蒲 鑫，侯榆青，朱 力，郭紅波

(西北大學信息科學與技術學院，陜西西安 710127)

隨著醫學影像技術的發展，功能成像已為大家所熟知，分子影像可以從細胞水平對疾病做出提前預警，并且能夠進行在體的醫學實驗，實時反映生物體的病理過程，這些優點使它成為了醫學界研究的熱點。生物發光斷層成像BLT(bioluminescence tomography)作為分子影像的一種成像方法，具有成本低、無輻射、操作簡單等優點［1］。在BLT的研究中，光源重建是個典型的逆向問題，即通過動物體表含噪聲的表面測量光子密度，結合生物體解剖結構信息和光學參數信息等，采用一定重建算法來確定其內部的光源分布，定位病灶位置［2］。光在生物組織中傳播時的多重散射，導致了BLT逆問題的病態性［3－4］，并且需要重建的光源信息量要遠大于邊界測量的光信息量，使得BLT成為典型的不適定問題。重建算法的核心就是求解這個不適定問題。為了解決這個問題，在文獻［5］中提出了一種光源可行區域方法，將其作為先驗信息，這種方法在很大程度上降低了其病態性。然而光源可行區域是根據光源表面分布和動物表面的異質結構來推斷的，當更深層的多光源存在的時候，這種推斷并不一定正確。為此文獻［6］中提出了一種多光譜的方法，將熒光素酶的輻射分為不同波段，增加了外部測量信息，提高了測量效率，在一定程度上降低了其不適定性。然而多光譜的數據量大，增加了垃圾信息，給BLT重建的計算效率帶來了一定壓力。基于此，學者們融合了可行區域、解剖信息以及多光譜等先驗信息以及測量噪聲，將重建問題轉換為一個最小二乘問題，利用正則化方法來求解。TIkhonov作為L2正則化的代表［7］首先提供一個相對平滑的解，但其重建效果卻存在偽影，降低了定位的準確性。后來，L1正則化引起了研究者的注意，由于其可以從較少的邊界測量得到較好的重建圖像，克服了L2范數求解的過平滑性［8］。L1范數是一個凸優化問題，可以采用現有的優化手段求解，比如一些追蹤算法［9］。另一方面研究表明，L1正則化算法即使在測量數據非常有限的情況下也能有很好的光源重建［10］。關于L1正則化的問題在壓縮感知領域已經有了廣泛的研究。

本文依據傳統壓縮感知理論，提出一種新穎的壓縮采樣匹配追蹤(CoSaMP)［11］算法，將光源重建問題看作一個基追蹤問題，并結合自適應有限元思想。應用算法本身的回溯特性，在重建過程中，對離散網格上各節點光信息向量不斷的比對，尋找包含最佳光源信息的點。與一般的L1正則化方法相比，該算法不僅提高了光源重建的質量，同時也提高了雙光源的重建速度。

1 基于壓縮采樣匹配追蹤的光源重建

1．1 重建模型

光在組織中的傳播可以由輻射傳輸模型來描述

對其進行有限元離散，拉格朗日二次插值，最后就可得到表面測量和動物體內光源分布的線性矩陣方程

加上測量噪聲，我們將其轉化為一個約束性最小二乘問題

A∈Rm×n是權值矩陣，b∈Rm是邊界測量的光密度，x∈Rn是生物體內光源分布，m?n，e是測量噪聲。

1.2 光源重建

當我們進行BLT的重建時，受采集設備采集方式的限制和光源分布的特點，不會對仿體全部表面進行探測，只會對感興趣的部位進行探測，所采集到的光能量分布相對于需要重建的全部光源能量分布表現出很大的稀疏性，這使得BLT的重建表現出稀疏重建的特性［7］。而壓縮采樣也是以信號的稀疏性和可壓縮性為前提，即信號當中只有s個非零系數，或者小系數經過變換之后都歸零，這是用CoSaMP方法來重建BLT的基礎。設s為信號的稀疏度，其數學定義如下

基于BLT的稀疏性，可以把重建過程看作一個信號分解問題，對A的列向量vi(1≤i≤N)進行歸一化，其集合就構成了一個完備字典D，然后就可以將測量數據看作是D中s個向量的線性組合。這s個向量的索引構成一個集合Γ，因此，b就可以表示為b=AΓxΓ，AΓ是Γ中A的索引的列向量組成的子陣。xΓ表示Γ中x的索引組成的向量。由于y在字典中的稀疏性，可以尋找最少的一組向量{vi|i∈ˉΓ}，使得殘差r=b－AˉΓxˉΓ最小。

在CoSaMP進行迭代計算時，每次迭代過程當中同時選取多個向量，而后又根據一定的條件刪除之前選取的一些向量，保留至少一個用作信號重建，它保證在每次選擇時候選集當中的向量數目都少于3 s個，迭代支撐集有s個，每次迭代的時候刪掉不多于2s個向量。

第k次迭代的結果，都是一個s－稀疏的近似解ak，它和待重建信號之間滿足下面的關系

v是由于噪聲原因而產生的能量損耗

1.3 算法描述

下面是用CoSaMP方法進行光源重建的迭代過程，對其描述如下:

輸入:剛度矩陣A，含噪聲的測量數據向量y，稀疏度s;

輸出:待估信號是x，一個s－稀疏逼近a。

1)初始解a0=0，初始殘差r0=y，計數器t=1;

2)計算信號代理值 c=ATrt－1;

3)鑒定成分大小，將支撐集合并 Ω=sup p(c2s)，T= Ω ∪ sup p(at－1);

4)利用最小二乘法估算信號 b|T=A+Ty，b|Tc=0來獲得逼近解at=bm;

5)更新殘差 rt=y－Aat;

6)停止條件滿足，輸出a=at，若否，則令t=t+1，返回第2)步。

2 實驗

2.1 實驗設置

為了驗證本文提出的方法對復雜動物模型的重建性能，在此我們使用美國南加州大學提供的通過冷凍切片、CT和PET數據構造的三維數字鼠模型［12］進行仿真重建實驗。截取其中高32mm的腹部軀干部分作為研究對象，如圖1(a)所示，此數字鼠圖譜(mouse atlas)的具體光學特征信息如表1所示。

表1 各器官光學參數Tab．1 Optical properties for the organs

對此小鼠軀干模型進行有限元離散、網格剖分，其中四面體單元個數為112 795，網格點數為21 277，邊界單元數為17 942，三角形17 942，邊緣單元數4 992，邊界元素點2 203，相對于總數，邊界點很稀疏，如圖1(b)所示。

圖1 三維數字鼠模型Fig．1 3D digitalmousemodel

2．2 實驗對比

一個好的算法，不僅能在單光源的時候精確的重建出光源，也能在雙光源的時候表現出良好的重建特性，為此我們設置了單光源和雙光源兩組實驗。在單光源實驗中，一個半徑為0.5 mm，高為1mm的圓柱狀光源被放入肝部，其中心位置為(18.1 mm，6.3 mm，15.4 mm)。初始的光源總能量和密度分別為0.785 nW和1nW/mm3。

我們用L1正則化方法中的截斷牛頓法(L1Ls)和本文所提出的CoSaMP算法對單光源進行重建，在不同的高斯噪聲水平下，對重建位置誤差、重建能量密度大小以及算法速度進行對比，如圖2所示，(a)(b)為2種算法初次網格上的重建效果，(c)(d)為2種算法自適應細分后的網格上的重建效果，圖中黑圈為真實光源位置。重建數據如表格2所示。

圖2 單光源的重建Fig．2 The reconstruction of single source

在雙光源實驗中，我們在數字鼠體內設置兩個光源，光源的各項參數和單光源相同。光源位置分別為(11.6，10.8，16.4)和 (11.6，6.3，16.4)，2種算法重建效果如圖3所示，(a)(b)為初次網格上的重建效果，(c)(d)為自適應細分后的網格上的重建效果，圖中最內側黑斑為真實光源位置，重建數據如表2表3所示。

表2 單光源重建數據Tab．2 The data of single source reconstruction

圖3 雙光源重建Fig．3 The reconstruction of single source

我們發現，在不同的噪聲水平下，在初次和細分后的網格上，CoSaMP方法的重建能量密度均要比L1Ls方法好，重建位置偏差更小，且在計算速度方面我們發現CoSaMP方法明顯要優于L1Ls方法。

從圖和數據中可以看出，在初次網格和自適應細分一次后，CoSaMP方法對2個光源位置的定位比較準確，重建的雙光源能量密度和重建位置都要優于L1Ls方法，這是因為CoSaMP算法受約束等距特性的啟發，在每次迭代中都要計算信號代理中前2s個元所對應的支撐集以及前次計算的逼近解，將它們合并，作為本次候選支撐集，最后用最小二乘法求出其系數，這種不斷“刪除”和“添加”的回溯思想提高了計算的精確性。

3 結語

在BLT問題中，由于生物組織的影響，使得光子的傳播過程變得異常復雜，不可避免地添加了很多噪聲，因此，選取一種既能穩定重建光源位置，又能很好抑制噪聲的方法是重建光源算法的核心。本文為了解決BLT重建中存在的困難，基于傳統的壓縮感知理論，提出了一種新穎的壓縮采樣匹配追蹤算法。并通過數字鼠數值仿真實驗得以驗證，結果顯示:在不同噪聲水平下，不論是單光源還是雙光源，這種計算方法都能夠快速而準確地實現光源重建，重建精度和速度都較L1正則化方法有所提高，并且實驗證明CoSaMP方法具有很好的魯棒性。

表3 雙光源重建數據Tab．3 The data of double source reconstruction

［1］ WILLMANN JK，van BRUGGEN N，DINKELBORG L M，et al．Molecular imaging in drug development［J］．Nature Reviews Drug Discovery，2008，7(7):591-607．

［2］ GU X，ZHANG Q，LARCOM L，et al．Three-dimensional bioluminescence tomography with model-based reconstruction［J］．Optics Express，2004，12(17):3996-4000．

［3］ CONGW，WANG G，KUMAR D，et al．Practical reconstruction method for bioluminescence tomography［J］．Opt Express，2005，13(18):6756-6771．

［4］ WANG G，Li Y，JIANG M．Uniqueness theorems in bioluminescence tomography［J］．Medical Physics，2004，31:2289．

［5］ SHEN H，GOLDSTEIN A S，WANG G．Biomedical imaging and image processing in tissue engineering［M］//Tissue Engineering，Berlin Heidelberg:Springer，2011:155-178．

［6］ DEHGHANIH，DAVISSC，JIANG S，et al．Spectrally resolved bioluminescence optical tomography［J］．Optics Letters，2006，31(3):365-367．

［7］ LU Y，ZHANG X，DOURAGHY A，et al．Source reconstruction for spectrally-resolved bioluminescence tomography with sparse a priori information［J］．Optics Express，2009，17(10):8062-8080．

［8］ GAO H，ZHAO H．Multilevel bioluminescence tomography based on radiative transfer equation Part 1:l1 regularization［J］．Opt Express，2010，18(3):1854-1871．

［9］ DONOHO D L，HUO X．Uncertainty principles and ideal atomic decomposition［J］．Information Theory，IEEE Transactions on，2001，47(7):2845-2862．

［10］ CANDESE J，ROMBERG JK，TAO T．Stable signal recovery from incomplete and inaccuratemeasurements［J］．Communications on pure and applied mathematics，2006，59(8):1207-1223．

［11］ NEEDELL D，TROPP JA．CoSaMP:Iterative signal recovery from incomplete and inaccurate samples［J］．Applied and Computational Harmonic Analysis，2009，26(3):301-321．

［12］ DOGDAS B，STOUT D，CHATZIIOANNOU A F，et al．Digimouse:a 3D whole body mouse atlas from CT and cryosection data［J］．Physics in Medicine and Biology，2007，52:577-587．