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

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

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

0 引言

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

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

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

1 光聲斷層成像原理

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

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

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

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

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

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

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

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

通過定義：

由式(3)，得：

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

圖1

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

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

延遲求和算法有兩種遍歷方式：基于傳感器的遍歷和基于重建像素的遍歷?；趥鞲衅鞯谋闅v方法是基于傳感器來計算傳感器在重建圖像上的每個時刻的信號區(qū)域，根據(jù)傳感器的空間響應(yīng)范圍，將強度疊加在重建區(qū)域上?；谥貥?gòu)像素點的遍歷是根據(jù)像素點計算與每個傳感器的相位角，看是否在傳感器的空間接收范圍內(nèi)。如果在范圍內(nèi)，計算與傳感器的距離、計算延時以及相應(yīng)的強度。最后，對每個范圍內(nèi)的傳感器強度進行累加，得到像素的最終光吸收強度。

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

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

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

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

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

圖2 信號波束形成圖

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

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

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

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

計算：

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

圖3 延遲求和流程圖

3 實驗結(jié)果與分析

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

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

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

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

圖4 原始模型圖A

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

圖5 基于DAS的重建結(jié)果

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

圖6 原始模型圖B

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

圖7 基于DAS的重建結(jié)果

4 結(jié)論

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