利用深度學習擴展非接觸擴散光學層析成像系統中探測器的動態范圍

趙 星，劉桐君，陳 平，王京凡，劉偉偉，2

（1. 南開大學 現代光學研究所，天津 300350；2. 天津市微尺度光學信息技術科學重點實驗室，天津300350）

1 引 言

擴散光學層析成像（Diffuse Optical Tomog?raphy，DOT）［1］是一種能夠快速無損探測組織吸收散射系數分布的技術，被廣泛應用于腦功能成像［2-3］，乳腺癌診斷［4-5］，關節成像［6-7］，農產品質量檢測［8］等領域。其本質是求解一個包含光子傳輸模型的逆問題，具有內在的非線性及病態性［9］。而獲取更多的信息用于重建能夠有效地提高DOT 重建圖像的空間分辨率和對比度，并消除偽影。高密度DOT 系統是最直接的方案，如Eggebrecht 等人［2］在 2014 年搭建的接觸式高密度DOT 系統，該系統對人腦部進行功能成像并和功能核磁共振進行對比，展示了高密度DOT 系統的有效性與準確性。雖然該系統具有很高的準確率，但系統組件非常多且復雜。非接觸式DOT 系統則能有效克服這個問題，其光源和探測器不再以光纖探頭的方式接觸式地分布于組織表面，光源是直接在空間中傳輸的激光束，探測器則使用EM-CCD［10-11］等面成像裝置，特點是能夠一次性大面積地采集組織表面的能流分布。Takamasa［12］等人使用納秒級脈沖激光和時間提取的圖像傳感器制作了非接觸式時域高密度DOT 系統，該系統無需使用密集的光纖探頭就可以輕松獲取高密度測量點。使用這套系統分別對模擬頭部吸收系數的仿體和真人的腦部進行了測量，研究了人腦在完成簡單計算任務時的活動情況。當然該系統也存在著限制，探測到的部分光信號會被具有較高吸收系數和散射系數的組織或介質衰減，這意味著隨著探測范圍的加大，光的衰減將進一步加強，數據的動態范圍也進一步變大。因此，在已有的非接觸式DOT 系統中，為了獲得足夠精確的探測數據，一般采用調整硬件參數的方法，如選擇較小尺寸或較小吸收系數的探測目標以減小對探測器動態范圍的需求，或選擇使用較為昂貴的高動態范圍CCD來滿足探測需要。雖然這些方法保證了測量精度，但卻限制了非接觸DOT 系統的應用空間及前景。使用算法修復有誤差的數據則能夠在保證DOT 重建精度的同時，降低非接觸DOT 系統對于探測器動態范圍的需求。

深度學習是機器學習的一種，適用于擁有大量數據的非線性問題的建模，已經應用在圖像分辨率增強［13］、模糊圖像恢復［14］、光纖傳輸后散斑恢復［15］、高光譜分類［16］和遙感圖像檢索［17］等領域，并展示了非常強的數據恢復能力和建模能力。本文采用深度學習的方法，使用可進行DOT 計算和重建的NIRFAST 開源程序對不同二維吸收散射特性的介質進行計算，通過數據截斷方法模擬低動態范圍探測器的信號獲取過程，生成大量訓練樣本。在此基礎上搭建了一維全連接網絡并完成了模型訓練，最后，通過仿真和實驗方法對模型的性能進行了分析及驗證。該方法能夠有效地修復受探測器動態范圍限制而獲得的誤差較大的探測數據，并成功重建出介質中吸收系數的分布特性。

2 基本原理

2.1 DOT 簡介及重建算法

DOT 是一種利用光通過吸收散射介質表面能流密度分布和光子傳輸模型計算重建介質內部吸收系數和散射系數分布的技術。其基本的測量過程示意圖如圖1 所示。圖中的介質具有一定的吸收散射特性，介質上表面的點表示光源的入射位置，下表面的點表示探測位置。當光依次從介質上表面的光源位置入射時，進入介質的光發生吸收及多次散射，最后從介質下表面出射。利用點探測器或面探測器探測從介質下表面出射的擴散光的能流密度，在得到各探測位置的一組數據后，利用DOT 求解算法即可計算重建介質內部的吸收系數分布。

圖1 擴散光學層析成像測量過程示意圖Fig.1 Measurement process of diffuse optical tomogra?phy

DOT 求解算法包含正問題及反問題的求解。正問題求解指的是根據已知的光源探測器分布以及介質參數分布，求解光子傳輸模型以獲取介質表面的能流密度分布。反問題求解指的是通過已知光源探測器分布、介質表面的出射光能流密度分布以及光子傳輸模型求解介質中的參數分布。由此可見，反問題是一個包含正問題的優化求解問題，求解算法流程如圖2 所示。

常見的正問題求解模型有麥克斯韋方程、輻射傳輸方程和擴散方程等。擴散方程是光源各向同性情況下從輻射傳輸方程簡化而來的，相對輻射傳輸方程來說具有計算難度低、計算時間快等優點。NIRFAST 軟件包［18］是一款被廣泛使用的DOT 算法開源軟件包，所使用的正問題求解模型即為擴散方程。這里以NIRFAST 軟件包使用的正問題求解模型及反問題算法為例介紹圖2 中的流程。首先，輸入介質的尺寸參數和出射表面測得的能流密度ФM，然后設定初始的吸收散射系數，代入軟件包中的擴散方程進行正向問題求解，得到計算的能流密度分布ФC，擴散方程一般為：

圖2 擴散光學層析成像重建算法流程Fig.2 Flow chart of diffuse optical tomography recon?struction algorithm

式中：k是約化散射系數，Ф是強度I對角度的積分，也就是能流密度，μa是吸收系數，ω是光的調制頻率，對于穩態DOT 來說，調制頻率為0，q0是源項。求解該方程一般使用數值方法，包括有限元、蒙特卡洛等。

如果計算所得和測量所得的能流密度分布的差值大于設定的閾值，則根據差值使用式（2）所示的列文伯格-馬夸特法確定迭代的參數更新值，進行反向問題求解。

式中：是正則系數，I表示單位矩陣，δФ代表能流密度測量值和計算值之間的差值，δμ代表吸收系數及散射系數的更新值，J 是雅閣比矩陣，是正向計算值對于參數的導數矩陣可以表示為：

式中：μa和k分別代表吸收系數和約化散射系數，ds為光源位置數量與探測點位置數量的乘積，N表示網格節點數。根據式（2）所得的參數更新值重新進行正向求解，直到能流密度計算值和測量值之間的誤差小于閾值。此時的吸收和散射系數分布即為計算重建結果。從式（3）中可以看到，在實際求解過程中，采用了能流密度數值的對數，這要求DOT 系統探測器必須具有較大的動態范圍，才能利用軟件算法實現高精度的計算重建。

2.2 深度學習網絡模型及損失函數

深度學習是機器學習的一種，它使用大量樣本進行模型訓練，建立起輸入樣本集和目標樣本集之間的模型。在本文中參與訓練的數據是一維數據，且輸入數據和目標數據長度相同。為了實現數據的訓練，選擇了全連接網絡作為深度學習網絡模型。理論上全連接網絡對于復雜問題的建模能力略好于卷積神經網絡［19］，且能夠靈活地調整網絡中每一層的大小，能夠更容易地實現下采樣和上采樣。使用的網絡模型示意圖如圖3所示。

圖3 全連接網絡示意圖Fig.3 Schematic diagram of fully connected network

網絡的輸入為一個1×40 的向量，經過3 次下采樣后變為1×10 的向量，之后再經過3 次上采樣變回1×40 的向量。訓練的目標是修復由于探測器動態范圍受限而測量誤差較大的數據，所以在實際訓練中，使用NIRFAST 開源程序對不同二維吸收散射特性的介質進行透射能流密度的計算，通過數據截斷方法模擬低動態范圍探測器的信號獲取過程，生成大量訓練樣本。考慮到這些數據動態范圍大的特點，為了使所有數據在訓練過程中具有相近的權重，同時為了使訓練結果適用于DOT 重建算法，采用式（4）所示的對數處理方法，對由NIRFAST 計算得到的透射能流密度和數據截斷后的透射能流密度進行預處理，分別得到網絡的目標數據和輸入數據，即：

式中：y為網絡目標數據和輸入數據，Ф為經過NIRFAST 計算出的能流密度及其截斷處理后的值。網絡的損失函數被設置為有權重的均方差（Mean Square Error，MSE），即有：

式中：L即為損失值，yture即訓練的目標值，M為數據長度，ytrain為經過網絡模型輸出的數據。該損失函數可以使訓練過程中絕對值差異較大的數據具備一致的模型訓練收斂性。

3 仿真及實驗結果

3.1 深度學習網絡模型訓練

在模型訓練中，吸收散射介質為65 mm×15 mm 的矩形仿體，參考了實際應用中的DOT 系統的光源密度和探測器密度［2-3，6］，把光源位置數設置為4 個，探測點數量設置為10 個，分別等間距地分布于介質上下表面，總探測次數等于光源與探測器位置數量的乘積，即40 次。訓練樣本中仿體背景的吸收系數分布于0.3~0.5 mm-1之間，異質體的吸收系數分布于0.4~1 mm-1之間。通過設置不同的背景吸收系數、異質體的大小、位置以及吸收系數等參數，計算出的能流密度分布在10-5~10-25內，經過式（4）預處理后，共產生了含有4 000個樣本的訓練集、100個樣本的驗證集以及7個樣本的測試集。

為了模擬低動態范圍探測器（動態范圍為103）的信號獲取過程，所有能流密度探測數據小于最大值千分之一的數據都會被人為截斷，作為一維網絡的輸入值，利用上述網絡模型和損失函數開展訓練。在實際訓練中，為了展示不同樣本量的訓練集對于網絡預測能力的影響，分別使用了1 000 個樣本和4 000 個樣本構成訓練集進行網絡訓練，并采用同樣的驗證集進行檢測驗證，訓練過程中訓練集和驗證集的損失值下降如圖4所示。

圖4 不同樣本數量訓練集的網絡訓練誤差Fig.4 Errors of training network with different training sets

從圖4 中可以發現，網絡訓練過程中，不論是1 000 樣本量還是4 000 樣本量，訓練集和驗證集的誤差均呈現下降趨勢，這反映了模型訓練的有效性。

3.2 仿真驗證結果

為對模型性能進行驗證和分析，在網絡模型訓練完成后，把沒有參與訓練的7 個測試樣本輸入到網絡中，可獲得如圖5 所示的其中3 個樣本的測試結果。圖5 展示了按照式（4）預處理前樣本能流密度的完整數據、截斷數據和修復數據。從圖中可以看到，網絡模型修復后的數據與截斷前的完整數據非常接近，這說明訓練好的網絡修復了被截斷的數據。為了進一步驗證網絡模型對低動態范圍探測器測量數據的修復效果，使用圖2 所示流程進行了DOT 重建，重建效果如圖6和圖7 所示。

圖5 不同樣本截斷數據的修復結果Fig.5 Recovery results of truncated data of different samples

圖7 對不同樣本修復后數據重建的吸收系數分布效果比較Fig.7 Comparison of reconstructed absorption coefficient distributions of different samples

圖6 為利用圖5 中所示的被截斷處理的數據重建得到的吸收系數分布結果。可以發現，不論仿體和異質體的參數如何設置，模擬低動態范圍探測器采集的數據重建結果具有相同的特點，即吸收系數分布中會出現一些高吸收系數（10 mm-1左右）的異常區域，這些高吸收區域的分布呈現空間對稱性。該結果說明使用低動態范圍探測器采集吸收散射光信號所產生的數據誤差嚴重影響了DOT 的重建效果，對DOT 計算重建來說是不容忽視的。

圖6 被截斷數據的重建吸收系數分布Fig.6 Reconstructed absorption coefficient distribution of truncated data

圖7 為分別利用經過圖5 中所示完整的數據和網絡模型修復后的數據重建得到的吸收系數分布結果。雖然圖7（b）和圖7（c）中吸收系數重建結果在對比度和空間分布上和圖7（a）中的目標值重建結果有偏差，但與圖6 中使用截斷數據重建得到的吸收系數分布結果相比，使用經過網絡模型修復后的預測數據能夠成功進行重建。為了進一步分析重建結果，使用均方根誤差（Root Mean Square Error，RMSE）和皮爾森相關系數（Pierson Correlation Coefficient，PCC）對圖7 中的吸收系數分布進行量化比較，均方根誤差越小，皮爾森相關系數越接近于1，這表示由網絡模型修復數據重建的吸收系數分布越接近于截斷前完整數據的重建結果。數據如表1所示。

表1 不同訓練集樣本量訓練的網絡修復后數據重建的吸收系數分布量化評價結果Tab.1 Quantitative evaluation results of absorption coef?ficient of recovery data of network trained by dif?ferent training sets

從表1 中可以發現，使用更多訓練樣本訓練出的模型具有更好的數據修復性能，而且修復的數據重建出的吸收系數和截斷前完整數據重建出的吸收系數是強相關的。

由此可見，使用深度學習方法對模擬動態范圍為103的探測器采集的有誤差的DOT 能流密度數據進行修復，能夠將探測器的動態范圍擴展至1020，經由DOT 重建算法對數據計算后，可以獲得較好的吸收系數分布重建結果。模擬結果說明，該方法可以有效擴展探測器的動態范圍，減小低動態范圍探測器采集信號對DOT 重建結果精度的影響。

3.3 實驗結果

在仿真分析的基礎上，進一步對利用深度學習方法擴展非接觸DOT 系統中探測器的動態范圍進行實驗驗證。首先，參考文獻報道［10］搭建了非接觸式穩態DOT 系統，系統示意圖和實物圖分別如圖8（a）和圖8（b）所示。整個實驗系統由三部分組成，分別是在一維方向移動的光源，具有一定吸收散射系數的仿體，以及對仿體后表面進行成像的探測器。

光源部分由光纖輸出的激光器（B&W TEK INC. BRM-785-0.55-100-0.22-SMA）、準直透鏡（Ocean Insight 74-VIS）、夾持器和電動平移臺（Thorlabs LTS300/M）組成。探測器部分由鏡頭和CCD（DAHENG MV-EM510M/C）組成，其中CCD 的灰階為8 位，動態范圍為256。實驗中仿體是采用環氧樹脂作為基質，摻雜印度墨水和氧化鋅粉末，并加入固化劑進行固化得到的長方體固體，仿體上鉆孔填充一定濃度的印度墨水溶液作為異質體，其橫截面結構尺寸及吸收系數分布如圖8（c）所示。仿體溶液配置過程中，使用吸收光譜法可以確定摻雜印度墨水后的環氧樹脂溶液吸收系數，通過調整印度墨水的摻雜濃度，確定溶液吸收系數為0.1 mm-1。根據氧化鋅粉末的參數以及Mie 散射理論［20］可計算獲得溶液散射系數的理論值，通過調整氧化鋅粉末的摻雜濃度，近似確定溶液散射系數為1 mm-1。異質體溶液配置過程中，采用與仿體相同的方法確定其吸收系數為0.3 mm-1。

圖8 實驗光路及仿體參數分布Fig.8 Experimental optical path and distribution of phan?tom parameter

在實驗測量過程中，光源及探測器的數量和分布情況同模擬仿真中一致。激光束依次從仿體前表面預設的位置入射，經仿體吸收散射后從后表面出射并被CCD 相機探測。為獲得仿體后表面出射光的能流密度分布，通過實驗標定法建立了仿體后表面能流密度和CCD 像素灰度值的關系［21］。按照實驗中仿體設定的吸收系數及散射系數進行計算，可知能流密度的理論動態范圍大約為1010，這遠大于實驗中CCD 的動態范圍。根據3.2 節中仿真模擬得到的結論，在這種情況下想要獲得無重建誤差的結果，需要對采集到的數據進行修復。

利用3.1 中訓練得到的網絡模型，對實驗采集到的能流密度數據進行修復，結果如圖9 所示。從圖9 可以發現，實際測量值的動態范圍大約為200，修復后數據的動態范圍被擴展到109左右。這說明深度學習方法修復了低動態范圍探測器采集的能流密度數據，并擴展了探測器的動態范圍。使用圖9 中的數據進行了吸收系數重建，得到的吸收系數分布如圖10 所示。

圖9 實驗采集值及經過網絡修復后的值Fig.9 Measured data and data after network recovery

從圖10（c）可以發現，使用實驗中探測得到的能流密度進行DOT 重建得到的吸收系數分布結果與圖8（c）中的預設參數相比有較大的誤差，而且吸收系數的空間分布與圖6 中仿真截斷后數據的重建結果相似。這從實驗上驗證了探測器動態范圍過小會導致探測到的能流密度的誤差較大，從而嚴重影響DOT 的重建精度。圖10（a）和圖10（b）中的結果則說明了盡管使用不同樣本量的訓練集訓練的網絡模型對于數據的修復能力有差異，但是都可以顯著提升DOT 重建吸收系數分布結果的精度。在此基礎上計算了圖10（a）和圖 10（b）中的結果與圖 8（c）中設定值的RMSE，分別為0.545 和0.542。雖然重建的吸收系數分布和圖8（c）的理論預設值之間存在一定的誤差，但考慮到實驗中運動精度、器件噪聲、仿體固化等不確定性因素的影響，該誤差是在可接受范圍內的。

圖10 實驗數據及其網絡修復后的值的DOT 重建結果對比Fig.10 Comparison of DOT reconstruction results of ex?perimental data and its predicted values

實驗結果表明，深度學習方法對實驗中動態范圍為256 的探測器采集的有誤差的DOT 能流密度數據進行修復后，能夠將探測器的動態范圍擴展至109，經由DOT 重建算法對數據計算后，吸收系數分布的重建結果準確。這從實驗上證明深度學習方法可以有效擴展探測器的動態范圍，減小低動態范圍探測器采集信號對DOT 重建結果精度的影響。

4 結 論

本文提出了一種利用深度學習擴展非接觸擴散光學層析成像系統中探測器動態范圍的方法。首先，使用NIRFAST 軟件包生成了一系列模擬低動態范圍采集的DOT 數據。其次，使用這些數據完成了一維全連接網絡模型的訓練。之后分別通過模擬仿真和實驗對該模型擴展非接觸DOT 系統中探測器動態范圍的能力進行了驗證。在模擬中，該方法修復了模擬動態范圍為103的探測器采集到的含有誤差的測試集數據，顯著減小了數據誤差對DOT 重建結果精度的影響，理論上將探測器的動態范圍擴展到1020；在實驗中，制備了實驗仿體并搭建了一套非接觸DOT 系統，使用動態范圍為256 的CCD 相機完成了DOT 系統的實驗數據采集。利用本文提出的方法對采集到的數據修復后，可重建得到與預設仿體吸收系數分布基本一致的重建結果，將DOT 系統中探測器的動態范圍從256 擴展到109。雖然模擬和實驗中通過網絡模型修復的數據重建出的吸收系數結果的精度還有待提高，但該方法能夠有效修復低動態范圍探測器所采集的有誤差的DOT 數據，減小它對DOT 重建結果精度的影響，為利用非接觸DOT 系統探測高吸收散射介質或大面積探測提供了一個可行的解決方案，使得非接觸式DOT 系統有了更廣闊的應用空間。

目前，該方法只對單層的仿體進行了驗證。如果進一步修改DOT 求解模型，使之能適用于多層不同介質組成的仿體，該方法還可以應用到更復雜的多層仿體或生物組織上。