結合Tikhonov正則化方法的近紅外漫射光血流成像技術

張曉娟，徐金榮，桂志國，白馨，左佳，劉祎，尚禹

1. 太原工業學院 電子工程系，山西 太原 030008；2. 中北大學 信息與通信工程學院，山西 太原 030051

引言

近紅外漫射光譜（NIRS）是一種無創檢測生物組織血氧含量的技術，為腦部缺血、骨骼肌痛和腫瘤等疾病提供診斷依據[1-3]。將NIRS技術擴展至成像領域，即漫射光學層析成像（DOT）技術，可以根據病變組織中氧含量與周圍組織的對比度，準確定位病灶位置[4-7]。測量的血氧含量，包括氧合血紅蛋白、脫氧血紅蛋白和總血紅蛋白濃度，可以評估組織氧供給量之間的平衡關系，是個靜態指標。血氧是通過血流對組織內血細胞進行氧氣供給的，血流變化可以動態反映組織微血管系統中血液動力變化趨勢。因此，血流的相對變化可以更早、更有效地預測疾病的產生與發展。

近年來，一種“動態”NIRS技術，即漫射相關光譜（DCS）技術快速發展，已經用于組織微血管血流變化的監測中。DCS技術利用光場的時間自相關函數（g1(τ)）計算紅細胞的運動狀態，進而無創地監測微血管中血流隨時間的變化值[8-9]。將DCS擴展至成像領域（DCT），就可以得到血流的空間圖像[10-11]。

與DOT技術類似，早期的DCT重建算法也是基于解析解的，同樣受限于組織的幾何形狀。最近幾年，有限元方法（FEM）也引入了DCT，并且已經應用在乳腺癌檢測中[12]。這兩種算法已經成為DCS和DCT的主流數值技術，但兩者共同的缺點是只能用自相關函數(g1(τ))曲線上的某一點作為原始數據進行重建，這將嚴重影響血流重建的精確性和穩定性。

在對當前方法深入分析的基礎上，本文引入了一種新型的N-階線性算法（NL算法）來實現血流成像。與傳統的解析表達式和有限元方法不同，NL算法不尋求偏微分方程的解，而是將g1(τ)的積分形式與其N階泰勒多項式相結合，并且通過光在組織內的傳輸信息來體現組織的非均勻性和不規則性[13]。此外，該算法還充分利用了多個時間點對應的g1(τ)值進行血流重建，這樣，可以有效增加圖像重建的準確性和穩定性。

將N-階線性算法應用到DCT血流成像中，由于測量數據遠遠少于待重建的圖像體素，其對應的線性方程組呈嚴重病態，因此，本文引入Tikhonov正則化方法進行圖像重建[14-15]。本文通過計算機仿真，在真實的磁共振人體頭部模型中進行了仿真驗證，重建過程中所需的光子路徑信息由蒙特卡羅仿真提供[16-17]。

1 方法

與DOT成像技術類似，DCT也需要在待測組織的表面放置若干光源—探測器（S-D）組合，用于收集逸出組織的光子，然后根據光子的自相關曲線進行血流的三維重建。

1.1 組織模型建立與S-D放置

本文采用的組織模型來自于開源的磁共振（MRI）數據庫（http://brainweb.bic.mni.mcgill.ca/brainweb/）中正常成人頭部影像。由MRI影像重建的頭部可視化模型如圖1a所示，光學探頭（包含S-D的光纖）放置于前額（藍色區域）。S-D的位置如圖1b所示，紅色為光源，共9個；黃色為探測器，共28個。每個光源發出的光子可以被任何一個探測器檢測到，因此共有252個S-D組合。光源和探測器的標號順序為從左往右，從下往上。

整個頭部模型的大小為18.1 cm×21.7 cm×18.1 cm，X、Y、Z的間隔均為0.05 cm，因此共有362×434×362個單元，這遠遠大于測量數據的個數（即252個S-D組合）。考慮光子對組織的入射深度和逸出范圍，取S-D覆蓋部分3 cm×3 cm×3 cm局部區域作為圖像重建的范圍 （即圖1中藍色陰影部分）。為進一步減少單元個數，將這部分區域重新抽樣，使體素大小為0.2 cm×0.2 cm×0.3 cm，則該局部區域單元個數縮減為15×15×10=2250個。為了檢驗算法的有效性，設置了兩種異常血流狀況（圖2～3）：一種是“十字形”異常血流組織，包含10個組織單元，圖2a為含“十字形”異常血流的頭部立體模型，圖2b為X-Y橫切面，圖2c為選定區域X-Y橫切面；另一種是“雙條形”異常血流，包含4個組織單元，圖3a為含“雙條形”異常血流的頭部立體模型，圖3b為X-Y橫切面，圖3c為選定區域X-Y橫切面。正常組織的血流值設為1×10-8cm2/s，異常組織單元的血流值設為5×10-8cm2/s，即異常組織的血流值為正常組織的5倍。

圖1 MRI影像重建的頭部可視化模型（a）及光源探測器放置位置（b）

1.2 MC仿真與自相關函數g1(τ)的計算

本文利用MC仿真（MCVM）來模擬光子在腦部的傳輸軌跡[15]，從而得到每個S-D組合檢測的光子數和光子傳輸路徑。光子MC仿真的參數設置如下：吸收系數μa=0.1 cm-1、散射系數 μs=80 cm-1、各項異性因子 g=0.9，折射系數n=1.4[15]。每個光源約107個光子出射，探測器采集到的光子信息包括光子坐標、傳輸方向、路徑長度和權重（即光子能量），這些信息用于計算光斑強度的自相關函數 g1(τ)。

其中P(s1,…sn)為光子傳輸路徑在n個單元中的歸一化分布，k0(i)為第i個單元的光波矢幅值，表示光子在第i個單元中的隨機步長，該值等于1/μs'(i)，τ為自相關函數的延時。＜Δri2(τ)表示第i個單元中散射粒子在時間段τ內的均方位移，其形式取決于具體的流動模型，根據以往的研究，布朗運動模型與實驗數據吻合較好，即 ＜Δri2(τ)＞=6αDB(i)τ。α為運動粒子與總粒子數的比值，αDB定義為生物組織中的血流值（BFI）[13-14]。

為了更有效地接近實際測量，在產生自相關函數g1(τ)的同時，加入了噪聲[8]，得到了含噪聲的自相關曲線，并利用含噪數據進行了血流重建。

1.3 N-階線性算法

將g1(τ)展開為N階泰勒級數g[13]：

并結合式(1)，可以得到一階(N=1，公式(3))和N階(N＞1，公式(4))的關于g1(τ)的展開式 ：

對 g1(τ)的 一 階 近 似 式 (3)，即為第j個S-D的相關曲線g1(τ,j)-1關于時間τ的斜率，表示為Sl(j)。

令A=CT，則血流值BFIs(αDB)可以直接由式(6)給出：

這 里，αDB=[αDB(1),…, αDB(n)]T，A=A(i,j)m×n，Sl=[ Sl(1),…, Sl(n)]T。

對g1(τ)的N階近似式(4)，由于方程兩邊均出現了αDB，因此需用迭代的方法進行求解：

公式(6)和(8)即為圖像重建的關鍵理論支撐[13-14]。

1.4 圖像重建方法

令X=αDB，B=Sl，則公式(5)和(7)可以由線性方程(8)表示：

這里，A=A(j,i)252×2250為所有S-D（共252組）針對所有組織單元（共2250個組織單元）得到的系數矩陣，X2250×1為 2250個組織單元，B252×1為得到的 252 條 g1(τ)曲線斜率，對這種嚴重不適定問題，一種較為有效的計算方法就是在求解過程中利用某些先驗信息進行約束，Tikhonov正則化方法即利用2-范數[16-17]，將問題轉化為：

X0為均勻組織的血流值BFIs（αDB）。λ是平衡線性方程組與懲罰項權重的正則化因子。Tikhonov迭代算法的公式如下所示：

這里，λ0, λ1, … ,λk一系列正則化因子，λ0為首次迭代的正則化因子，采用通常的方法，將L曲線的拐點作為λ0[18-19]，λ1, … ,λk則通過公式 (12)迭代得到 ：

通過使用這種非平穩迭代方法，vk更接近最優解，收斂速度更快。Tikhonov正則化的偽代碼如下：

給定：A，b，初始化：v0=ATB（反投影方法），λ0取值為L曲線的拐點，設定迭代的條件，計算：

非負約束：

更新：

另外，由于組織的厚度和S-D的距離問題，某些組織單元中只有很少的光子經過，使用這些數據將對圖像重建帶來不穩定性，所以，我們將組織單元中光子數少于15的單元剔除。這樣，2250個原始單元減少至1424個；得到的線性方程組為 A252×1424X1424×1=B252×1。重建完成后，剔除掉的單元以均勻組織的血流值（1×10-8cm2/s）補充。

2 結果分析

本節主要闡述了N-階線性算法結合Tikhonov正則化對選定組織區域的重建效果，在含噪聲和不含噪聲的條件下，分別對“十字形”和“雙條形”異常BFI進行了重建，并通過均方根誤差（RMSE）對重建效果進行了客觀評價。

2.1 不同形狀異常BFI重建效果

Tikhonov正則化算法重建的結果由圖2和圖3所示。

圖2 “十字形”異常血流值位置建模和重建效果圖

可以看出，不含噪聲的血流圖像重建效果精確程度高，噪聲顯著影響圖像重建的質量，尤其是深度較大的血流異質物。但是，無論是否有噪聲，都可以清楚地判斷出異常血流組織的形狀和位置。例如，“十字形”異常血流的形狀可以辨別出來（圖2）。對于異常血流的位置，Tikhono-DCT算法也區分出了距離僅為0.4 cm的兩個條形異常血流組織，即使兩個條形異常組織的寬度只有0.2 cm，也可以精確重建（圖3）。

圖3 “雙條形”異常血流值位置建模和重建效果圖

通過定量分析，不含噪聲的重建圖像中，重建的最大異常值為4.65×10-8（接近真實的異常值（5×10-8））組織位于圖像正中間，邊緣部分逐漸減小，如圖2f所示。含噪聲重建的血流圖像，最大異常值為3.93×10-8，精確性有所降低，這符合常識。本次試驗為五次加噪之后的平均值，測試次數越多，均值越穩定，但是耗時越長。

2.2 均方誤差評價標準

為了對圖像重建效果進行定量分析，我們使用了均方根誤差方法，如式(16)所示：

均方根誤差反映了重建圖像與原圖像的相似度，該值越小，圖像重建效果越好。異常組織深度增加，重建圖像的RMSE柱狀圖，見圖4。圖4(a)為不含噪聲時，“十字形”的RMSE柱狀圖，從圖中可以看出，當異常血流組織位于頭皮下5 mm的全局（2250個組織單元）RMSE為13.07%，局部（10個組織單元）RMSE為23.59%；7 mm深度時，全局RMSE為13.81%，局部RMSE為35.83%。隨著深度的增加，全局誤差基本保持不變，但是局部誤差急劇增加，深度為11 mm時，局部誤差達到54.67%。全局誤差小，局部誤差大的原因在于Tikhonov算法之前的投影計算將整體均勻化了，被均勻化后，局部異常血流值明顯下降，導致局部誤差較大，而本次仿真中局部（10個組織單元）體積相比全部體積（2250個組織單元）很小，對整體影響不大，因此整體誤差較小。與投影計算相反，Tikhonov正則化可以有效地保持邊緣，最終使異常組織的血流得以恢復。

圖4 異常組織深度增加，重建圖像的RMSE柱狀圖

圖4(b)為含噪聲時，“十字形”的RMSE 柱狀圖，從圖中可以看出，當異常血流組織位于頭皮下5 mm的全局（2250個組織單元）RMSE為20.90%，局部（10個組織單元）RMSE為44.56%，遠高于不含噪聲時的誤差值。圖中全局誤差依然基本保持不變，分析原因如下：

（1）在Tikhonov正則化之前，采用反投影的方法作為初值進行正則化，相當于對異常血流值進行了平滑作用，使異常血流值降低，造成了局部誤差大。

（2）在本文的設置中，S-D數量（252）與體素個數（2250）的比率較小，甚至因為距離太遠，某些S-D之間探測到的光子數很少，甚至沒有光子被探測到，這個比率更小，導致光場自相關函數g1(τ)信號較弱，因此，重建圖像的RMSE值偏高。

（3）對g1(τ)信號添加了一定的噪聲，尤其是S-D距離較大時，自相關函數g1(τ)信號較弱，而噪聲數據大，對重建結果影響很大。

需要指出的是：44.56%是異常血流的誤差，而全局血流的誤差為20.90%；因為異常血流與周圍均勻組織的對比度較大，因此，這些誤差的具體數值不影響本文的結論，即我們提出的Tikhonov正則化圖像重建方法可以檢測出異常的血流。

另外，“雙條形”異常血流組織（深度5 mm）的全局（2250個組織單元）RMSE為9.19%，局部（4個組織單元）RMSE為37.87%。全局的RMSE在深度增加時依然基本維持不變，而且比“十字形”的全局RMSE小，原因是它的異常組織個數更少，只有4個。

3 結論

本文基于NL算法，并結合Tikhonov正則化算法提出了一種新穎的血流重建方法（Tikhonov-DCT）。根據這個方法，利用MC仿真得到的光子傳輸信息，被測組織的幾何形態和內部特征等信息得以體現，同時該方法也能夠充分利用測得的光學信號（即自相關曲線）。同時，Tikhonov正則化可以提高血流異質物的重建精確性和穩定性。

在采用Tikhonov正則化方法進行血流重建的過程中，發現正則化參數λ對重建效果影響很大，因此，在測試了一系列λ值以后，選定L-曲線的拐點作為λ0值，并采用迭代的方法，逐步減小λ值，以增加重建血流值的保真度。在今后的研究中將對正則化參數λ進行進一步優化，例如，使用自適應的方法來優化參數λ的數值，從而改善血流重建的質量

從重建效果可以看出，異常組織的定位較為準確，形狀邊緣保持較好，可以有效地區分異常組織與正常組織。與原始模型相比，Tikhonov-DCT重建圖像的局部RMSE比較大，隨著深度增加，急劇上升。全局RMSE較小，且隨著深度增加保持微量增加。今后的工作將引入壓縮感知思想優化S-D的位置和數量，增加圖像重建的準確性，降低RMSE值，這些算法也將在仿體中進行實驗驗證。