齊夢(mèng)雨， 趙麗麗， 劉 欣， 嚴(yán)壯志

1.上海大學(xué)通信與信息工程學(xué)院，上海200444

2.上海大學(xué)上海生物醫(yī)學(xué)工程研究所，上海200444

光聲成像技術(shù)是一種新興的生物醫(yī)學(xué)成像技術(shù)，它的基礎(chǔ)是光聲效應(yīng).當(dāng)生物組織被脈沖激光照射時(shí)，被照射區(qū)域因光能量的沉積而產(chǎn)生熱膨脹，從而產(chǎn)生微弱的光聲信號(hào).通過分布在生物組織周圍的超聲探測(cè)器陣列，或者按照指定軌道運(yùn)動(dòng)的單個(gè)超聲探測(cè)器[1]，即可探測(cè)并采集光聲信號(hào).利用采集到的光聲信號(hào)，通過一定的重建算法反演出待成像組織的光吸收分布，進(jìn)而得到該組織的光學(xué)特性.近年來，光聲成像技術(shù)在生物醫(yī)學(xué)的很多領(lǐng)域取得了顯著的成果，例如與納米傳感器結(jié)合的光聲成像技術(shù)[2-3]，生物分子的光聲成像和功能成像[4]以及多模態(tài)成像方式聯(lián)合的光聲成像[5]等.

高分辨率的圖像重建往往需要處于不同位置的超聲探測(cè)器采集足夠充分的光聲信號(hào).為了提高數(shù)據(jù)采集的速度并保證采集數(shù)據(jù)的完備性，不同位置的超聲探測(cè)器進(jìn)行數(shù)據(jù)采集時(shí)一般采用多角度并行采集的方式.在這樣的條件下研究出很多重建算法[6]，其中以迭代法最為常見.大多數(shù)迭代重建算法采用的是正則化L2范數(shù)的方法，但這種方法會(huì)引起圖像局部特征丟失，進(jìn)而導(dǎo)致過平滑現(xiàn)象[7]的出現(xiàn).圖像含有的不連續(xù)分量（例如輪廓和邊緣）越多，這種過平滑現(xiàn)象也就越明顯.因此，為了獲得高分辨率的重建結(jié)果，傳統(tǒng)的迭代法要求測(cè)量數(shù)據(jù)具有較高的完備性，這就大大增加了數(shù)據(jù)采集及成像時(shí)間.如何縮短數(shù)據(jù)采集時(shí)間或者只利用較少測(cè)量數(shù)據(jù)便可重建出理想的光聲圖像，已成為光聲成像技術(shù)面臨的重要挑戰(zhàn)之一.

近年來，壓縮感知理論在高分辨率光聲成像領(lǐng)域發(fā)展迅速.根據(jù)該理論，如果圖像具有可壓縮性或者稀疏性，則能夠在低于奈奎斯特采樣率的條件下，以較高的分辨率重建出該圖像[8].由于稀疏重建的算法大多基于L1范數(shù)正則化，因此越來越多的L1正則化方法被應(yīng)用于光聲成像領(lǐng)域[9].與基于L2范數(shù)正則化的重建方法相比，基于L1范數(shù)正則化的重建方法會(huì)保留圖像的更多細(xì)節(jié)，從而在測(cè)量數(shù)據(jù)不完備的情況下也可重建出分辨率較高的光聲圖像[10].其中，文獻(xiàn)[11]提出了一種基于L1范數(shù)正則化的交替方向算法（alternating direction method,ADM），可以利用很少的測(cè)量數(shù)據(jù)重建出較高分辨率的光聲圖像.由于這種算法可以節(jié)省大量的數(shù)據(jù)采集時(shí)間，因此成像速度得到了很大的提高.然而，對(duì)于更少的測(cè)量數(shù)據(jù)和更為復(fù)雜的組織結(jié)構(gòu)，這種算法依然存在著精度和分辨率不足的局限.

本文提出了一種改進(jìn)的交替方向加權(quán)算法，可以進(jìn)一步提高在欠采樣情況下光聲圖像重建的質(zhì)量.交替方向加權(quán)算法基于傳統(tǒng)交替方向算法，在迭代過程中引入了結(jié)構(gòu)先驗(yàn)信息進(jìn)行約束，顯著抑制了偽影并降低了噪聲，引導(dǎo)迭代計(jì)算過程中的原始變量和對(duì)偶變量向正確解收斂.實(shí)驗(yàn)結(jié)果顯示，在采集數(shù)據(jù)較少的情況下，利用本文算法重建出的三維光聲圖像具有更好的精度和更高的空間分辨率.

1 三維光聲圖像的重建算法

1.1 光聲成像的前向模型

根據(jù)光聲效應(yīng)原理，目標(biāo)組織經(jīng)過脈沖激光照射產(chǎn)生光聲波并形成聲場(chǎng).在均勻和非粘滯的聲介質(zhì)中，光聲波的傳播公式為

式中，c為聲速，p(r,t)和H(r,t)分別描述的是位置r處、t時(shí)刻的聲壓和由于激光照射而產(chǎn)生的熱能量沉積，Cp為等壓比熱容，β為等壓膨脹系數(shù).

如果滿足如下熱學(xué)限制條件和壓力限制條件：1）脈沖激光持續(xù)時(shí)間小于組織的熱擴(kuò)散時(shí)間；2）脈沖激光的脈寬小于聲波在目標(biāo)組織中傳播所需的時(shí)間，則熱能量沉積函數(shù)H(r,t)就可以寫成光源的時(shí)間包絡(luò)函數(shù)與光能量沉積分布的乘積，即

式中，I(r,t)為光源能量的時(shí)間包絡(luò)函數(shù)，μa(r)為位置r處的光吸收系數(shù).如果將I(r,t)近似為沖激函數(shù)δ(t)，則根據(jù)式(1)求得位置r處、t時(shí)刻的瞬時(shí)聲壓為

式(3)為光聲信號(hào)的時(shí)域解析式.如果將介質(zhì)進(jìn)行離散，同時(shí)將超聲探測(cè)器檢測(cè)到的光聲信號(hào)表示為向量b，那么基于式(3)，可以將光聲成像的前向模型寫成離散表達(dá)式[12]為

式中，W稱為測(cè)量矩陣，可以理解為是光聲信號(hào)在介質(zhì)中的傳播，而p0為t=0 時(shí)刻聲場(chǎng)的初始值.

1.2 基于壓縮感知的圖像重建

一般而言，由于探測(cè)數(shù)據(jù)是不充分的，因此測(cè)量矩陣W是一個(gè)欠定矩陣.所以通過直接求解式(4)得到初始聲壓p0幾乎是不可能的.光聲圖像在通過一定的變換后具有可壓縮性或者稀疏性[7,13]，根據(jù)壓縮感知理論，光聲圖像可以通過較少的探測(cè)數(shù)據(jù)進(jìn)行重建.對(duì)于這樣的一種稀疏變換ψ:x=ψp0，求解初始聲壓p0可以轉(zhuǎn)化為最小化問題

在這種情況下，式(5)中的限制條件也發(fā)生了變化

求解式(7)的問題可以轉(zhuǎn)化為求解以下問題的最優(yōu)解：

式中，σ是正則化參數(shù).當(dāng)σ和ε都接近于0 時(shí)，式(7)和(8)將會(huì)趨近于式(5).

1.3 傳統(tǒng)交替方向算法

從數(shù)學(xué)的角度上，傳統(tǒng)的交替方向算法[14]在對(duì)L1范數(shù)問題進(jìn)行求解時(shí)，使用原始對(duì)偶算法（primal-dual algorithms）在每次迭代的過程中對(duì)原始變量和對(duì)偶變量同時(shí)進(jìn)行更新.式(8)的增廣拉格朗日子問題有如下形式：

式中，y為拉格朗日乘子，β(β ＞0)是懲罰項(xiàng).如果已知，那么通過式(9)就可以求得對(duì)于r來說，如果給定x=xk且y=yk，式(9)取得最小值時(shí)滿足

如果給定r=rk+1且y=yk，最小化式(10)相當(dāng)于求解以下問題

通過對(duì)式(11)的求解[15]就可以得到

式中，τ(τ ＞0)是近似參數(shù)，并且有

最后，在已知r=rk+1和x=xk+1的情況下，下一次迭代中的乘子進(jìn)行更新

在迭代過程中不斷更新這些原始變量和對(duì)偶變量，在一定迭代次數(shù)后就可以得到最終解.

2 改進(jìn)的交替方向加權(quán)算法

雖然基于L1范數(shù)正則化的傳統(tǒng)交替方向算法在圖像重建的過程中保留了圖像的更多細(xì)節(jié)，但是在探測(cè)器數(shù)目不足和測(cè)量數(shù)據(jù)較少的情況下，傳統(tǒng)交替方向法的重建圖像依然存在著比較嚴(yán)重的偽影.如何進(jìn)一步消除這些偽影，繼續(xù)提升重建圖像的分辨率和準(zhǔn)確度是本文的主要關(guān)注點(diǎn).由式(10)、(12)和(14)可知，在傳統(tǒng)交替方向法的迭代過程中，原始變量和對(duì)偶變量會(huì)不斷進(jìn)行更新并最終得到最優(yōu)解.如果某次迭代過程中的原始變量和對(duì)偶變量陷入了局部最優(yōu)解，就會(huì)使圖像在某些位置上產(chǎn)生噪點(diǎn).由于測(cè)量數(shù)據(jù)較少和噪聲干擾等原因，下一次迭代可能無法對(duì)這些噪點(diǎn)進(jìn)行有效糾正，從而一定程度上導(dǎo)致了重建結(jié)果中部分偽影的產(chǎn)生.如果能夠利用目標(biāo)物的結(jié)構(gòu)、形狀或大致位置等先驗(yàn)信息對(duì)傳統(tǒng)交替方向法的迭代過程進(jìn)行約束，對(duì)于那些出現(xiàn)在反常位置上且大概率為噪點(diǎn)的值，在下一次迭代過程中進(jìn)行適當(dāng)?shù)囊种?，將有助于抑制偽影從而提高重建圖像分辨率.

2.1 結(jié)構(gòu)先驗(yàn)信息

交替方向加權(quán)算法通過引入先驗(yàn)結(jié)構(gòu)信息來改善傳統(tǒng)交替方向算法的迭代過程，使迭代結(jié)果向正確的結(jié)果快速收斂.在光聲成像中，光聲信號(hào)或圖像由于自身的結(jié)構(gòu)或者組成成分，往往包含天然的分組或分塊先驗(yàn)信息.在交替方向加權(quán)法中，結(jié)構(gòu)先驗(yàn)信息包含了目標(biāo)組織的大致位置或形狀等信息，具體呈現(xiàn)形式為三維權(quán)重矩陣Wi

式中，Ti為權(quán)重矩陣中目標(biāo)組織大概率出現(xiàn)的位置，n為圖像離散后的網(wǎng)格數(shù)目，a和b為權(quán)重值，且a ＞1,0＜b≤1.某點(diǎn)處對(duì)應(yīng)的權(quán)重值大小代表了算法在迭代過程中對(duì)該點(diǎn)實(shí)際值的抑制程度，當(dāng)權(quán)重值大于1 表示該點(diǎn)的實(shí)際值在迭代過程中趨向于被抑制，值越大代表被抑制的程度越大；當(dāng)權(quán)重值小于1 表示該點(diǎn)的實(shí)際值在迭代過程中趨向于被增強(qiáng)，值越小代表被增強(qiáng)的程度越大；當(dāng)權(quán)重值等于1 時(shí)表示該點(diǎn)的實(shí)際值維持原狀態(tài)，既不被增強(qiáng)也不被抑制.由此可知，權(quán)重值a和b的大小實(shí)際上控制著迭代過程的速度和算法的收斂步長(zhǎng).在實(shí)際應(yīng)用中，往往會(huì)根據(jù)待成像組織的天然結(jié)構(gòu)、解剖特征，或者借助于其他成像技術(shù)來確定Ti的大致位置.例如文獻(xiàn)[16]在腦部成像時(shí)借助于腦部的解剖結(jié)構(gòu)以及不同部位的活躍程度對(duì)腦部進(jìn)行了分組并確定Ti的位置，而文獻(xiàn)[17]在熒光斷層成像的研究中利用結(jié)構(gòu)成像技術(shù)來獲得Ti的大致位置.a和b具體數(shù)值將作為算法的輸入?yún)?shù)加入算法的迭代過程中.本文在Ti位置處設(shè)置b=1，其他位置處設(shè)置a=10.

2.2 基于結(jié)構(gòu)先驗(yàn)的交替方向加權(quán)算法

本算法中的結(jié)構(gòu)先驗(yàn)信息以三維權(quán)重矩陣W′i的形式并作為式(8)的一個(gè)約束條件，加入到每一次L1范數(shù)正則化的過程中.這個(gè)約束條件可以描述為

具體來說，加入結(jié)構(gòu)先驗(yàn)信息后式(8)將具有以下形式：

對(duì)式(17)的求解過程也就是利用傳統(tǒng)交替方向法求解的過程.該式將作為算法外循環(huán)過程中的一次迭代，第l次外循環(huán)的結(jié)果將和結(jié)構(gòu)先驗(yàn)信息一起作為第l+1 次外循環(huán)的迭代約束條件：

通過合理設(shè)置循環(huán)迭代容差來確定交替方向加權(quán)算法的終止條件，其中，內(nèi)循環(huán)迭代容差為tol1，外循環(huán)迭代容差為tol2，且tol1與tol2決定算法的總迭代次數(shù).

2.3 本文算法主要步驟

本文算法的主要步驟如下：

輸入測(cè)量矩陣W和測(cè)量數(shù)據(jù)b，內(nèi)循環(huán)和外循環(huán)終止條件tol1=e-4tol2=e-3.

輸出初始聲壓分布的ψ稀疏逼近x.

初始化x=0，設(shè)置權(quán)重矩陣W′i=Wi，內(nèi)循環(huán)迭代次數(shù)k=1，外循環(huán)迭代次數(shù)l=1.

步驟1使用原始對(duì)偶算法計(jì)算式(17)，并通過式(10)、(12)和(14)更新原始變量和對(duì)偶變量r,x,y.

步驟2如果滿足返回步驟1 繼續(xù)進(jìn)行內(nèi)循環(huán)迭代；否則，令xl=xk，并繼續(xù)執(zhí)行步驟3.

步驟3令

步驟4對(duì)式(17)加入結(jié)構(gòu)先驗(yàn)約束，更新權(quán)重矩陣，即令

步驟5如果滿足返回步驟1；否則外循環(huán)結(jié)束，輸出xl.

3 實(shí)驗(yàn)結(jié)果與討論

為了說明交替方向加權(quán)算法在欠采樣條件下進(jìn)行光聲圖像重建的優(yōu)越性，本文進(jìn)行了仿真實(shí)驗(yàn)，將分別使用最小平方QR（least squares QR，LSQR）算法、傳統(tǒng)交替方向算法和本文提出的交替方向加權(quán)算法進(jìn)行光聲圖像的重建.其中，LSQR 算法是基于L2范數(shù)正則化進(jìn)行迭代重建，是一種求解大型線性方程組的經(jīng)典方法.為了全面驗(yàn)證算法性能，將通過兩個(gè)不同的方面進(jìn)行仿真實(shí)驗(yàn)：1）在欠采樣重建實(shí)驗(yàn)中，通過減少探測(cè)器數(shù)目并利用更少的測(cè)量數(shù)據(jù)進(jìn)行重建，以驗(yàn)證本文算法在欠采樣條件下進(jìn)行三維光聲圖像重建的優(yōu)勢(shì)；2）在噪聲魯棒性實(shí)驗(yàn)中，對(duì)測(cè)量數(shù)據(jù)加入不同強(qiáng)度的噪聲后進(jìn)行重建，通過重建結(jié)果的比較來考察本文算法的噪聲魯棒性.仿真實(shí)驗(yàn)中的三維權(quán)重矩陣根據(jù)仿體的解剖學(xué)結(jié)構(gòu)進(jìn)行構(gòu)建.探測(cè)器都按照理想的超聲探測(cè)器進(jìn)行處理.所有的仿真實(shí)驗(yàn)在MATLAB 平臺(tái)上進(jìn)行，光聲信號(hào)的傳播和聲場(chǎng)的仿真使用了開源工具箱k-Wave[18].

3.1 欠采樣條件下的光聲圖像重建

首先建立如圖1(a)所示的仿體幾何模型，仿體分為3 個(gè)部分，其中外部的大圓柱體為組織介質(zhì)，內(nèi)部大小相同且一上一下放置的兩個(gè)小圓柱體為1 號(hào)和2 號(hào)目標(biāo)物.組織介質(zhì)的底面直徑和高均為30 mm，總體積為15×15×30π mm3；內(nèi)部的兩個(gè)目標(biāo)物體的底面直徑和高分別為6 mm 和4 mm，體積均為3×3×4π mm3.介質(zhì)和目標(biāo)物的光吸收系數(shù)大小分別設(shè)置為Φ1= 1.0 mJ/mm3和Φ2= 5.0 mJ/mm3.48 個(gè)高斯光源均勻分布在以大圓柱體中心為球心、直徑為45 mm 的球體表面.光源點(diǎn)和光擴(kuò)散的過程通過有限元方法進(jìn)行仿真.有限元模型總共有1 834 個(gè)格點(diǎn)和8 249 個(gè)正四面體網(wǎng)格.在聲場(chǎng)的仿真和圖像的重建過程中，仿體則被離散總數(shù)為27 000 的均勻網(wǎng)格.整個(gè)仿體聲速均勻，聲速大小設(shè)置為1 500 m/s.

圖1 仿體模型與探測(cè)器陣列(單位:mm)Figure1 Geometry of phantom and position of the transducers (unit:mm)

首先利用24 個(gè)探測(cè)器的測(cè)量數(shù)據(jù)進(jìn)行3 種重建算法的重建，將這24 個(gè)探測(cè)器每6 個(gè)為一組分為4 組，每組均勻分布在仿體z=3 mm，z=11 mm，z=19 mm，z=27 mm 4 個(gè)平面的圓柱體表面，如圖1(b)所示.為了便于比較，圖2中選取了重建結(jié)果中縱向通過1 號(hào)和2號(hào)兩個(gè)目標(biāo)物的y= 15 mm 平面.從圖2(b)～(d)的重建結(jié)果中可以看到，LSQR、傳統(tǒng)交替方向法以及本文算法都較清晰地重建出了兩個(gè)目標(biāo)物.這說明3 種算法在測(cè)量數(shù)據(jù)較完備的情況下都能取得較為理想的重建結(jié)果.但相比之下，LSQR 算法的偽影更為嚴(yán)重.此外作為一種L2正則化的重建算法，LSQR算法重建出的兩個(gè)目標(biāo)物由于算法產(chǎn)生的過平滑現(xiàn)象而導(dǎo)致邊緣較模糊.圖2(c)中基于L1正則化的交替方向法則保留了圖像更多的細(xì)節(jié)，兩個(gè)目標(biāo)物的邊緣更加清晰且呈現(xiàn)為兩個(gè)大小相同的亮團(tuán).在圖2(d)中，本文算法除了重建出兩個(gè)大小一致的亮團(tuán)外，背景噪聲和偽影被大幅地抑制也使重建結(jié)果具有更好的空間分辨率.這體現(xiàn)出了結(jié)構(gòu)先驗(yàn)信息在迭代過程中對(duì)于背景部分的約束作用，兩個(gè)目標(biāo)之間的背景部分根據(jù)結(jié)構(gòu)先驗(yàn)信息而設(shè)置了較小的權(quán)重，所以在迭代過程中這些部分產(chǎn)生的偽影和噪聲會(huì)被不斷削弱.

為了比較3 種算法在欠采樣條件下的重建結(jié)果，將每層的6 個(gè)探測(cè)器減少至2 個(gè)，同樣分布在z= 3 mm，z= 11 mm，z= 19 mm，z= 27 mm 4 個(gè)平面上，如圖1(c)所示.相比于使用24 個(gè)探測(cè)器，8 個(gè)探測(cè)器不僅在采集數(shù)據(jù)的數(shù)據(jù)量上較小，而且由于只在8 個(gè)位置處進(jìn)行采集，所包含的位置信息也更少.這種條件下使用3 種重建算法得到的結(jié)果如圖3(b)～(d)所示.從圖中可知，3 種算法的重建結(jié)果都出現(xiàn)了不同程度的惡化.對(duì)比圖2(a)和圖3(a)可知，從LSQR 算法的重建結(jié)果中只能分辨出下方的1 號(hào)目標(biāo)物，上方的2 號(hào)目標(biāo)物則被嚴(yán)重的偽影干擾而無法分辨.LSQR 算法的過平滑現(xiàn)象由于欠采樣而加劇，導(dǎo)致了很多圖像細(xì)節(jié)的丟失.這一點(diǎn)在基于L1正則化的傳統(tǒng)交替方向法中得到了一定程度的優(yōu)化.從圖3(c)可知，傳統(tǒng)交替方向法得到的重建結(jié)果雖然也存在著比較嚴(yán)重的偽影，但是上方2 號(hào)目標(biāo)物的部分細(xì)節(jié)仍然被保存了下來.相比之下，本文算法的重建結(jié)果雖然由于欠采樣仍然產(chǎn)生了一定的偽影和噪點(diǎn)，但是兩個(gè)目標(biāo)物的位置準(zhǔn)確，大小也基本一致.由于結(jié)構(gòu)先驗(yàn)信息的約束，三維權(quán)重矩陣中權(quán)重值較大位置處的偽影在迭代過程中持續(xù)被抑制，在重建結(jié)果上體現(xiàn)為背景噪聲很小，分辨率更優(yōu).

圖2 使用24 探測(cè)器進(jìn)行不同算法的重建結(jié)果對(duì)比(y =15 mm)Figure2 Comparison of reconstructed image with 24 sensors from different method (y =15 mm)

圖3 使用8 探測(cè)器進(jìn)行不同算法的重建結(jié)果對(duì)比(y =15 mm)Figure3 Comparison of reconstructed image with 8 sensors from different method (y =15 mm)

對(duì)于使用8 個(gè)探測(cè)器的欠采樣條件下的重建結(jié)果，除了選取圖3中的y=15 mm 截面外，在圖4和5 的重建結(jié)果中還選取了通過目標(biāo)物中心的z= 6 mm 和z= 24 mm 兩個(gè)平面進(jìn)行展示.對(duì)比圖4(b)～(d)可知，在z= 6 mm 平面上，使用LSQR 算法得到的重建結(jié)果水平分辨率最差且偽影嚴(yán)重，圖4(c)中使用傳統(tǒng)交替方向算法得到的重建結(jié)果則對(duì)偽影有了一定程度的削弱.而在圖4(d)中，使用本文算法得到的重建結(jié)果在z=6 mm 平面上有效地抑制了偽影，水平分辨率有了顯著提升.z= 24 mm 平面的重建結(jié)果與z= 6 mm類似，本文算法的重建結(jié)果有效抑制了偽影，因此水平分辨率相較于LSQR 和傳統(tǒng)交替方向法有明顯提升，如圖5(d)所示.對(duì)比圖4(b)、圖5(b)以及圖4(c)、圖5(c)，無論是使用LSQR 算法還是使用傳統(tǒng)交替方向法，其重建結(jié)果在兩個(gè)平面上的水平分辨率都存在著差別，且在z=6 mm 平面上的水平分辨率都明顯優(yōu)于各自在z=24 mm 平面的水平分辨率.此外，LSQR 算法和傳統(tǒng)交替方向法在這兩個(gè)平面上重建出的目標(biāo)物也存在著大小和形狀上的差異.相比之下，本文算法的重建結(jié)果在兩個(gè)平面上不僅擁有相似的水平分辨率，而且能夠重建出大小基本一致的目標(biāo)物，這也體現(xiàn)了本文算法在欠采樣條件下的穩(wěn)定性.

為了進(jìn)一步衡量不同算法的重建性能，通過將不同算法的重建圖像歸一化后，與兩個(gè)目標(biāo)物所在的真實(shí)位置相比，得到了不同算法的均方誤差（mean square error,MSE）和峰值信噪比（peak signal to noise radio,PSNR），如表1所示.

圖4 不同算法重建結(jié)果的水平分辨率對(duì)比(z =6 mm)Figure4 Horizontal resolution comparison of reconstructed image by different method(z =6 mm)

圖5 不同算法重建結(jié)果的水平分辨率對(duì)比(z =24 mm)Figure5 Horizontal resolution comparison of reconstructed image by different method(z =24 mm)

表1 3 種重建算法的性能比較Table1 Performance comparison by three reconstruction algorithms

從表1中可知，在y= 15 mm 平面上，與LSQR 算法相比，雖然傳統(tǒng)交替方向法的結(jié)果MSE 較低而PSNR 較高，不過這種差別并不明顯.反觀本文算法的重建結(jié)果，MSE 值明顯下降（約為90）而PSNR 值出現(xiàn)了較大程度的上升（約為4 dB），說明本文算法的重建結(jié)果有著最優(yōu)的垂直分辨率.在z=6 mm 平面上，傳統(tǒng)交替方向法的MSE 值與LSQR 算法相比有了一定程度的下降（均為32 左右），同時(shí)PSNR 提升了約1 dB，說明傳統(tǒng)交替方向法的重建圖像質(zhì)量方面要略優(yōu)于LSQR 算法.本文算法由于在傳統(tǒng)交替方向法基礎(chǔ)上通過結(jié)構(gòu)先驗(yàn)信息進(jìn)行進(jìn)一步優(yōu)化，所以MSE 值出現(xiàn)了明顯下降而PSNR 值出現(xiàn)了明顯的上升（約為6 dB），這充分證明了本文算法對(duì)于2 號(hào)目標(biāo)物的重建結(jié)果，在水平分辨率上比前兩種方法有了大幅提升.此外，傳統(tǒng)交替方向法在z=24 mm 平面上的PSNR 值反而略低于LSQR 算法，說明傳統(tǒng)交替方向法的重建結(jié)果在水平分辨率上并不穩(wěn)定，而本文算法的重建結(jié)果在z= 24 mm平面上依舊有著最低的MSE 和最高的PSNR，充分說明了本文算法的穩(wěn)定性.

綜上，在欠采樣的條件下，LSQR 方法和傳統(tǒng)交替方向法重建出的目標(biāo)物在位置和大小上都存在著或多或少的偏差，而本文算法則在正確的位置上重建出了兩個(gè)大小相同的目標(biāo)物，和預(yù)期結(jié)果相符.此外在這3 種算法中，使用本文算法重建出的結(jié)果在多個(gè)平面上都具有最小的MSE值和最大的PSNR 值，同時(shí)具有最好的垂直和水平分辨率.圖2～4 的(b)和(c)中目標(biāo)物周圍的大面積偽影在圖2～4 的(d)中基本消失，亮團(tuán)準(zhǔn)確地定位在兩個(gè)目標(biāo)物的真實(shí)位置處，這也造成了視覺上的較大差異，充分表明本文算法中引入的結(jié)構(gòu)先驗(yàn)信息起到的糾正目標(biāo)物位置，抑制偽影和提高重建結(jié)果分辨率的作用.

3.2 噪聲的魯棒性實(shí)驗(yàn)

算法良好的噪聲魯棒性在其實(shí)際應(yīng)用中具有很大意義.為了研究噪聲對(duì)不同算法重建結(jié)果的影響，分別在測(cè)量數(shù)據(jù)中加入不同強(qiáng)度的高斯噪聲，其作用是模擬真實(shí)實(shí)驗(yàn)中的系統(tǒng)噪聲.本實(shí)驗(yàn)的參數(shù)和設(shè)置與前述相同，測(cè)量數(shù)據(jù)是使用24 個(gè)探測(cè)器采集得到.三維光聲圖像的重建過程分別在信噪比（signal to noise radio,SNR）為80 dB 和120 dB 的條件下進(jìn)行，其中信噪比越小，說明加入噪聲的強(qiáng)度越大；反之信噪比越大，則說明加入的噪聲強(qiáng)度越小.3種算法在不同噪聲條件下的重建結(jié)果如圖6～8 所示.

圖6 不同噪聲條件下LSQR 算法的重建結(jié)果Figure6 Reconstructed image using LSQR method under different noise conditions

圖7 不同噪聲條件下ADM 算法的重建結(jié)果Figure7 Reconstructed image using ADM under different noise conditions

從重建結(jié)果中可以看到，隨著SNR 逐漸變小，3 種算法在圖像的分辨率上都出現(xiàn)了不同程度的下降.根據(jù)圖6(a)～(c)可知，隨著噪聲的增加LSQR 算法重建結(jié)果中的偽影迅速增多，在SNR= 120 dB 時(shí)噪點(diǎn)和偽影已經(jīng)嚴(yán)重影響了目標(biāo)物的分辨率，而在SNR= 80 dB 的條件下，兩個(gè)目標(biāo)物都由于偽影和噪點(diǎn)的影響已經(jīng)不再呈現(xiàn)團(tuán)狀.而在圖7(b)和(c)中，雖然隨著SNR 的不斷減小，傳統(tǒng)交替方向法重建出的目標(biāo)物附近出現(xiàn)了較嚴(yán)重的偽影以及拖尾，但依舊保持了目標(biāo)物的大致形狀和位置信息，這里再次體現(xiàn)出了基于L1正則化方法在保留圖像細(xì)節(jié)和邊緣方面的優(yōu)勢(shì).在使用了本文重建算法的圖8(b)和(c)中，由于加入噪聲的緣故，兩個(gè)目標(biāo)物在相對(duì)的方向上均出現(xiàn)了一定的偽影，尤其當(dāng)信噪比下降至80 dB 時(shí)這種現(xiàn)象更加明顯，但是相比圖6(c)和圖7(c)而言，仍能夠明顯分辨出兩個(gè)目標(biāo)物的位置和大小.在本文算法的迭代過程中，包含著結(jié)構(gòu)先驗(yàn)信息的權(quán)重矩陣在遠(yuǎn)離目標(biāo)物的位置處具有較大的權(quán)值，所以這些部位產(chǎn)生的噪聲在多次迭代后逐漸被抑制甚至消除，這也是圖8(c)中偽影主要集中在兩個(gè)目標(biāo)物附近以及兩個(gè)目標(biāo)物的中間部位、而在遠(yuǎn)離目標(biāo)物的位置十分微弱的原因.從圖6(c)、圖7(c)以及圖8(c)的對(duì)比中看到，這些產(chǎn)生在目標(biāo)物及其附近位置處的偽影一般不會(huì)導(dǎo)致分辨率急劇惡化.綜上，相比LSQR 和傳統(tǒng)交替方向法，本文算法在不同的噪聲條件下都可以得到分辨率較高的重建結(jié)果，因此對(duì)于噪聲具有較好的魯棒性.

圖8 本文算法在不同噪聲條件下的重建結(jié)果Figure8 Reconstructed image using proposed method under different noise conditions

4 結(jié) 語

光聲成像是一種融合了超聲成像和純光學(xué)成像的多模態(tài)成像技術(shù)，兼具高成像深度和高對(duì)比度的特點(diǎn)，應(yīng)用場(chǎng)景多樣，因此具有光明的發(fā)展前景.然而，在成為一項(xiàng)成熟的且可普遍用于臨床的成像技術(shù)之前，光聲成像仍面臨諸多挑戰(zhàn).其中最主要的就是因數(shù)據(jù)采集量較大導(dǎo)致的數(shù)據(jù)采集時(shí)間過長(zhǎng)以及較苛刻的硬件要求.為此，本文提出了一種L1范數(shù)正則化的交替方向加權(quán)算法.該算法對(duì)傳統(tǒng)交替方向算法進(jìn)行了改進(jìn)，通過設(shè)計(jì)仿真實(shí)驗(yàn)對(duì)比了不同算法在欠采樣條件和噪聲條件下的重建結(jié)果，驗(yàn)證了本文算法的優(yōu)勢(shì).實(shí)驗(yàn)結(jié)果證明，交替方向加權(quán)算法可以利用較少的測(cè)量數(shù)據(jù)重建出分辨率和精度都較高的光聲圖像，同時(shí)也具有較好的噪聲魯棒性.這對(duì)于減少光聲成像在實(shí)際應(yīng)用中的數(shù)據(jù)采集時(shí)間，提高成像速度具有積極意義.下一步工作將考慮本文算法在真實(shí)實(shí)驗(yàn)中的應(yīng)用，并對(duì)本文算法的性能和應(yīng)用場(chǎng)景等進(jìn)行更深入的研究.