基于DAS的光聲斷層圖像重建算法研究

張瑛，王浩全，李子圣，王小婷，李凱豐，侯清

（中北大學 信息與通信工程學院，山西太原，030051）

0 引言

近年來，醫學成像領域有了很大的進步。光學成像的優點是對比度很高，但是空間分辨率比較低，成像深度受到光傳播深度的限制；超聲成像雖然能獲得生物組織一定深度的高分辨率圖像，但是圖像的對比度比較低。光聲層析成像(Photoacoustic Tomography, PAT)能夠結合光學成像和超聲成像的優點[1～3]，逐漸成了醫學成像方面的研究熱點之一。

光聲斷層成像的重要步驟是圖像重建，常用的重建方案是利用組織邊界處記錄的壓力信息來獲得初始壓力分布。1995年，Krugger等人首次將光聲應用到生物組織成像領域[4]，在后來的時間中，國內外很多研究小組都開始將目光聚焦到光聲成像領域。在國外，研究光聲成像的小組有，(1)華盛頓大學的汪立宏教授小組，汪教授小組在光聲斷層成像方面做了很多研究[5]。(2)倫敦大學學院的Beard, Paul實驗組，該小組專注于研究生物醫學光譜成像，提出用一個基于高分子薄膜的Fabry-Perot干涉儀的光學檢測光聲信號的成像系統[6]。在國內，研究光聲成像的小組有復旦大學汪源源教授團隊和華南師范大學邢達教授團隊，他們在光聲成像系統開發以及光聲圖像重建等方面做了很多研究[7]。除此之外，廈門大學等許多所高校也在從事光聲成像技術的研究，并且取得了優秀的科研成果[8～10]。

在光聲圖像重建實際應用中，針對圖像重建需要快速處理等問題，本文采用延遲求和（DAS）算法，對圖像進行實時重建。

1 光聲斷層成像原理

PAT中的正向問題涉及收集給定初始壓力分布情況下，組織邊界上的壓力數據。聲波在生物組織中的傳播可以使用波動方程建模，如式(1)所示：

式中，c是介質中的聲速，β是熱膨脹系數，pC是比熱， (,)Hdt是單位時間內單位體積的吸收能量，空間位置由d表示。對于固定聲源，吸收能量可以寫成在應力約束條件下，光源的時間部分可用δ函數表示為式(1)中的相當于解決式(2)的初值問題[11]：

超聲換能器在位置0r處所探測的聲壓信號也即光聲信號作為公式(2)的解，表示為：

公式(3)表明，位置r時刻r的光聲信號是沿著半徑為ct的球面積分的初始聲壓分布的時間導數。不同的重建算法其實就是利用換能器探測到的通過不同方法去計算成像組織的初始聲壓分布

2 基于延遲求和的光聲圖像重建算法

當激發激光脈沖的持續時間短到足以滿足聲學和熱約束條件時，由于光激發而發射的壓力場隨空間r和時間t的變化可以表示為：

其中 ()′Hr是組織中每單位體積的吸收能量，Γ是無量綱Grüneisen參數，c是聲速。在下文中省略了常數項因為它不影響基于模型的重建。沿半徑為的球面S(r,t)進行積分。式(4)中正向模型的離散化是通過考慮笛卡爾網格上的N個體素來完成的，表示包圍所有光聲源的體積。每個圖像體素的位置由ri表示。空間H(r′)中任意位置處的吸收能量近似為插值函數K(r′)的加權疊加移動到不同的體素位置ri，即：

其中ih是體素i處的吸收。那么，式(4)可以表示為：

通過定義：

由式(3)，得：

Δ是體素大小，r=(x,y,z)。插值函數的非零值僅存在于體素i的鄰域，即當r′距離體素i小于一個體素時，hi僅影響H(r′)。因此，式(8)中的曲面積分僅為當曲面與曲面相交時的Δ鄰域，如圖1所示。讓測量位置r和體素位置ir之間的距離用s表示。因為體素大小Δ通常比s小幾個數量級，s中的球面積分，ir體素的鄰域可以通過平面積分來近似。在從ir到 ′r的方向上，積分值取決于距其表面d（用S′表示）到體素ir的距離，可以用參數化兩個角度α和β。因此，當ct=s-d時式(6)可以重新表述為：

圖1

圖1 笛卡爾網格上光聲正演模型的三維離散化。球面積分面S（r，t）與體素i的?-鄰域的交點。體素i和測量位置分別用ir和r表示。ir和r之間的距離用s表示，ir和球面之間的距離由d表示。方位角和仰角分別用α和β表示。

其中A是表示特定重建問題的模型矩陣。A的列表示與網格的每個體素相關的光聲信號。重建問題包括嵌入吸收矢量x，對于該矢量，理論模型與實驗測量的壓力信號b更好地匹配。

延遲求和算法有兩種遍歷方式：基于傳感器的遍歷和基于重建像素的遍歷。基于傳感器的遍歷方法是基于傳感器來計算傳感器在重建圖像上的每個時刻的信號區域，根據傳感器的空間響應范圍，將強度疊加在重建區域上。基于重構像素點的遍歷是根據像素點計算與每個傳感器的相位角，看是否在傳感器的空間接收范圍內。如果在范圍內，計算與傳感器的距離、計算延時以及相應的強度。最后，對每個范圍內的傳感器強度進行累加，得到像素的最終光吸收強度。

延遲求和法是通過式(11)來反演組織內的光吸收分區情況：

本文采用的算法根據探測器與生物組織實驗點之間的距離，計算超聲傳輸到探測器的延遲，最終確定生物組織發射光波時聲源的強度，并將每個探測點處的波束形成信號求和。在延遲求和這一算法中，主要使用超聲換能器檢測實驗樣本和生物組織的聲壓信號，根據傳感器位置與需要重建生物組織位置之間的距離差，使生物組織某一點吸收的光量被反轉。

本文實驗采用線性陣列探測器，如圖2所示。假設光聲信號由線性陣列探測器檢測，波束形成信號和DAS的接收信號之間的關系可以定義為：

式中，SDAS(x,y)是探測點(x,y)處的波束形成信號，S(i,τ)是時間t處第i個探測器元件的接收信號。延遲時間t(x,y,i)表示在探測點(x,y)處產生的聲波傳播到第i個探測器元件的時間，它可以表示為：

式中，l(x,y,i)是探測點(x,y)和第i個探測器元件之間的距離，c是聲速。該算法通過選取適當的數據點的平均值形成圖像像素后，可以使用式(11)獲得重建圖像。除像素點(x,y)外，與 PA 源（圖 2 中的S(x,y)）具有相同距離的所有像素都將受到強信號的影響。由于有效信號與不必要的數據相加，因此，在重建圖像中不可避免地會出現旁瓣和偽影。

圖2 信號波束形成圖

由圖2表示，位置（x，y）處的波束形成信號與第i個探測器元件接收信號之間的關系。旁瓣標記在PA源S的兩側。

本文的延遲求和算法采用的基于像素的遍歷方法，其偽代碼如下：

輸入：傳感器信號pn(t)，其中t代表時刻，n代表傳感器序號；N為傳感器個數；初始延遲t；聲速c：傳感器位置;權重ω= 1;nx,ny輸出圖像橫縱坐標范圍；傳感器空間響應范圍thea。

輸出：組織分布圖像A(r)，其中r代表像素位置。

計算：

延遲求和流程圖如圖3 所示。

圖3 延遲求和流程圖

3 實驗結果與分析

本文對光聲成像的成像環境進行了空間特征仿真，設置相關參數以便進行光聲數據采集：設置成像網格為256×256pixels，分辨率為150μm；設置設聲速為1540 m/s，聲衰減為0.5 dB/MHz/cm，Grüneisen為0.15；設置的線性陣列是由128個間距為200微米的探測器組成，將這一陣列放置在固定平面，探測器中心頻率為5MHz。

將光聲數據采集信號應用于真值圖像計算。在光聲數據采集信號中加入高斯白噪聲，來模擬實際環境。之后采用80%帶寬脈沖響應函數（IRF）對每個采集到的光聲信號進行濾波。

以上仿真采用的仿真軟件為Matlab R2019a，在處理器 為Inter(R)Core(TM)i5-4210H、CPU@2.90GHz、內 存8GB的個人計算機上進行的。

圖4(a)為原始模型圖A（像素191×191），內有一個小男孩在圖片右側站著，手拿燈籠。圖4(b)為所有探測器探測到的模型圖A的聲壓信號。

圖4 原始模型圖A

圖5是模型圖A基于DAS的重建效果，可以看出，圖像銳化了重建邊緣，輪廓清晰，能夠清楚地看到“小男孩手拿燈籠”這一景象，重建效果良好。

圖5 基于DAS的重建結果

圖6 (a)為原始模型圖B（像素256×256），內為散放的大米。圖6(b)為所有探測器探測到的模型圖B的聲壓信號。

圖6 原始模型圖B

圖7是模型圖B基于DAS的重建效果，可以看出，重建圖像效果沒有較上一組（模型圖A）重建效果好，存在相鄰微觀結構之間的混疊問題。這是由于來自一個目標的信號可能會受到來自相鄰目標信號的影響，從而導致圖像失真和混疊，并惡化圖像質量。這是由于在圖像信息過多時，使用DAS算法重建圖像通常會出現高旁瓣和強烈偽影[12～13]。

圖7 基于DAS的重建結果

4 結論

通過上述實驗結果可以說明，DAS基于對波束形成信號的簡單延遲和求和計算，能夠快速響應，非常適合實時光聲斷層成像。然而，當圖像信息過多時，在重建過程中有可能使用了一些不必要的數據進行求和，會使得基于DAS的光聲斷層圖像重建結果有可能出現高旁瓣和強烈偽影。解決高旁瓣和強烈偽影這一問題，是DAS算法接下來的一個研究方向。