采用稀疏測量數(shù)據(jù)的有限角度光聲層析成像的研究進(jìn)展

孫正，閆向陽

(華北電力大學(xué)電子與通信工程系，河北保定071003)

0 引 言

光聲層析(Photoacoustic Tomography, PAT)成像是一種基于生物組織光聲效應(yīng)的非電離功能成像方法，其成像參數(shù)是組織的光吸收系數(shù)和散射系數(shù)，可以實(shí)現(xiàn)高分辨率和高對比度的軟組織深層成像[1]，PAT 成像可用于冠狀動脈粥樣硬化性疾病、乳腺癌和消化系統(tǒng)病變等的早期診斷。

PAT 成像的靈敏度和分辨率取決于組織的光吸收對比度、激光脈沖持續(xù)時間、超聲檢測系統(tǒng)的幾何結(jié)構(gòu)和帶寬以及圖像重建算法等。其中高質(zhì)量的圖像重建是提高PAT 成像分辨率的重要途徑，目前較為成熟的圖像重建算法包括濾波反投影法(Filtered Back-Projection, FBP) 、 時 間 反 演 法(Time-Reversal, TR)、基于傅里葉變換的算法和代數(shù)迭代法等[2]。這些算法一般都需要對光聲信號進(jìn)行全角度的掃描測量，以得到完備的光聲數(shù)據(jù)集。但是在實(shí)際應(yīng)用中，特別是內(nèi)窺式光聲成像，如血管內(nèi)光聲成像[3]，由于管腔內(nèi)封閉成像幾何的特殊性，并且受成像導(dǎo)管的機(jī)械結(jié)構(gòu)、空間位置及成像時間等的限制，超聲探測器往往只能在封閉的管腔內(nèi)進(jìn)行有限角度的掃描，導(dǎo)致采集到的光聲數(shù)據(jù)不完備。而根據(jù)稀疏的光聲測量數(shù)據(jù)精確計(jì)算未知光聲信號的過程是不穩(wěn)定的，難以求解其高頻成分。本文對在有限角度掃描條件下，采用在稀疏位置測量的不完備光聲信號重建PAT 圖像的主要方法進(jìn)行綜述，討論各方法的優(yōu)勢和不足，并展望未來可能的發(fā)展方向。

1 有限角度掃描PAT 成像

1.1 PAT 成像原理

PAT 成像原理示意圖如圖1 所示，采用短脈沖激光照射生物組織，組織吸收光能量后受熱膨脹產(chǎn)生瞬時壓力，并向外輻射寬帶(帶寬為10 kHz～100 MHz)超聲波，即光聲信號[1]。聲壓的幅值與脈沖激光的強(qiáng)度成正比，反映組織的光吸收特性。超聲換能器接收來自不同方向、不同位置的光聲信號，送入計(jì)算機(jī)后采用合適的算法即可反演得到組織內(nèi)部初始聲壓的空間分布，進(jìn)而得到空間光吸收分布圖，直觀顯示組織的內(nèi)部結(jié)構(gòu)。在此基礎(chǔ)上還可估算組織的光學(xué)特性參數(shù)(光吸收系數(shù)和散射系數(shù))的空間分布，反映組織的功能成分，即定量PAT(Quantitative PAT, qPAT)[4]。

圖1 PAT 成像原理示意圖Fig.1 Schematic diagram of PAT imaging

1.2 有限角度掃描問題

目前，多數(shù)PAT 圖像重建算法都需要對光聲信號進(jìn)行全角度的掃描測量(如二維圓周掃描)，以獲得完備的光聲測量數(shù)據(jù)集。但在實(shí)際應(yīng)用中，探測設(shè)備尺寸、機(jī)械結(jié)構(gòu)和空間位置等會使超聲換能器的掃描角度受限。有限角度PAT 成像的掃描方式主要包括兩種類型：有限角度圓周掃描和直線掃描。此外，一些不規(guī)則的掃描幾何，如未閉合的曲線掃描也屬于有限角度掃描。由于投影的角度范圍不能滿足數(shù)據(jù)完備性的條件，因此會導(dǎo)致重建圖像中出現(xiàn)嚴(yán)重的偽影和失真，降低圖像質(zhì)量。

此外，在某些應(yīng)用中，受成像目標(biāo)的幾何形狀、解剖位置、曝光時間和數(shù)據(jù)檢測時間等的限制，超聲探測器只能采集到不完備的稀疏光聲數(shù)據(jù)。例如連續(xù)光聲成像掃描，為了在不增加產(chǎn)品成本和保持高分辨率的同時加快數(shù)據(jù)采集速度，必須減少測量次數(shù)。稀疏測量會導(dǎo)致成像目標(biāo)的重建圖像沿著目標(biāo)中心連線的方向擴(kuò)散，并且擴(kuò)散的區(qū)域隨著掃描范圍的減小而增大，導(dǎo)致圖像模糊、細(xì)節(jié)缺失和偽影現(xiàn)象嚴(yán)重。若采用最小分辨間距來衡量重建圖像的分辨率，采用掃描角度衡量掃描范圍的大小，則最小分辨間距與掃描角度之間呈指數(shù)關(guān)系[5]：

其中：γ 是最小分辨間距；θ 是掃描角度；γ0是全方位掃描時的最小分辨間距；A 和B 是待定系數(shù)。

2 有限角度PAT圖像重建的主要方法

對于有限角度掃描采集的不完備光聲數(shù)據(jù)集，直接應(yīng)用傳統(tǒng)的圖像重建算法會導(dǎo)致圖像中存在盲區(qū)、模糊和偽影等問題。目前，用于有限角度稀疏測量數(shù)據(jù)的PAT 圖像重建方法主要有：插值法、基于圓弧Radon 變換的方法、基于Gerchberg-Papoulis外推的補(bǔ)償方法、迭代重建法、頻域反卷積法、基于壓縮感知的方法和基于深度學(xué)習(xí)的方法等。

2.1 插值法

針對由不完備數(shù)據(jù)造成的圖像偽影現(xiàn)象，最簡單的方法是插值法，即根據(jù)采集的光聲信號，利用特定的插值函數(shù)計(jì)算出未測位置上的信號值，進(jìn)而重建圖像[2]。插值法重建結(jié)果對比如圖2 所示，采用基于三維平面掃描的插值算法[6]可以提高合成的光聲信號的準(zhǔn)確性，降低計(jì)算成本。與常規(guī)的反投影法相比，它能夠減小旁瓣和偽影，而且由于考慮了光的衰減，因此可以實(shí)現(xiàn)連續(xù)重建。插值法的精度一般不高，而且還需要解決聲速不均勻問題和考慮超聲探測器尺寸的影響。

圖2 插值法重建結(jié)果對比[6]Fig.2 Comparison of reconstruction results with interpolation method[6]

2.2 基于圓弧Radon 變換的方法

基于圓弧Radon(Circular Arc Radon, CAR)變換的數(shù)值反演算法是一種高效的非迭代算法，它能夠解決全圓周Radon 變換反演算法中第一類Volterra積分方程得不到穩(wěn)定解的問題，其原理[7]是：首先，使用梯形積分(Trapezoidal Product Integration, TPI)法反演求解CAR 變換中出現(xiàn)的第一類Volterra 積分方程；然后采用截?cái)嗥娈愔捣纸?Truncated Singular Value Decomposition, TSVD)法對積分核的離散近似病態(tài)矩陣B 進(jìn)行分解，得到其逆矩陣B-1；最后，利用CAR 逆變換得到原始聲壓分布。在有限角度掃描的情況下，由于圓弧邊界混入大量偽影以及圓弧角度的限制，可能導(dǎo)致圖像中一些目標(biāo)的邊緣不可見。為此，可對矩陣B 進(jìn)行加權(quán)處理，以減少有限角度掃描造成的偽影。矩陣B 加權(quán)處理對重建結(jié)果的影響如圖3 所示。

圖3 矩陣B 加權(quán)處理對重建結(jié)果的影響[7]Fig.3 Effect of weighting B on the reconstruction results[7]

2.3 迭代重建算法

2.3.1 同步迭代重建算法

采用基于全變分(Total Variation, TV)的迭代重建算法可以大幅減少所需的測量數(shù)據(jù)，但是缺少某些角度的投影數(shù)據(jù)可能導(dǎo)致圖像邊緣模糊和明顯的偽影，特別是在數(shù)據(jù)嚴(yán)重不足時[8]。此外，還可采用代數(shù)重建算法(Algebraic Reconstruction Techniques, ART)，與FBP 相比，ART 可以顯著提高圖像的分辨率和對比度[9]，但是其計(jì)算復(fù)雜度較高，一定程度上會影響重建的速度和質(zhì)量。

為了解決上述問題，可采用同步迭代重建技術(shù)(Simultaneous Iterative Reconstruction Technique,SIRT)及其改進(jìn)算法(Modified SIRT, MSIRT)[10-12]。同步迭代就是以各像素對應(yīng)于所有迭代變化的平均值作為像素的修正值。SIRT 的原理是：首先，設(shè)定初始重建圖像；然后，計(jì)算通過圖像中各像素點(diǎn)的所有投影數(shù)據(jù)，以此修正圖像；最后，檢驗(yàn)是否滿足收斂準(zhǔn)則。ART 和SIRT 的重建時間如表1 和圖4 所示，相較于ART，SIRT 的迭代收斂速度更快，能夠以較少的迭代次數(shù)實(shí)現(xiàn)高精度的圖像重建，但重建時間較長。為此，需對SIRT 進(jìn)行改進(jìn)。

表1 ART 和SIRT 的重建時間對比[12]Table 1 Comparisonof reconstruction time of ART and SIRT[12]

圖4 ART 和SIRT 的收斂速度比較(λ 為松弛因子)[12]Fig.4 Comparison of convergence speed between ARTand SIRT (λ is the relaxation factor)[12]

MSIRT 方法中每次迭代的是各像素的修正值，即累加經(jīng)過該像素的所有投影數(shù)據(jù)的誤差值，而不是只與一條投影數(shù)據(jù)有關(guān)，可有效抑制測量數(shù)據(jù)中的噪聲。此外，通過線性搜索調(diào)整速度更新的步長，還可以提高迭代收斂的速度[13]，ART 和MSIRT 的性能和重建結(jié)果分別如表2 和圖5 所示。

表2 ART 和MSIRT 的性能對比[13]Table 2 Comparison of performance of ART and MSIRT[13]

圖5 ART 和MSIRT 的重建結(jié)果對比[13]Fig.5 Comparison of reconstruction results with ART and MSIRT[13]

2.3.2 與FBP 結(jié)合的迭代重建算法

(1) 結(jié)合有限域FBP 的實(shí)時迭代重建方法

FBP 的原理簡單、速度快且易實(shí)現(xiàn)。但對于有限角度掃描獲得的不完備數(shù)據(jù)，F(xiàn)BP 會導(dǎo)致圖像的空間分辨率和對比度降低，使目標(biāo)的形狀發(fā)生模糊和變形。與FBP 相比，SIRT 可以大大減少重建圖像中的偽影，并在數(shù)據(jù)不完整以及含有噪聲的情況下重建更高質(zhì)量的圖像，但其計(jì)算復(fù)雜度較高，所需的計(jì)算時間遠(yuǎn)大于FBP，可能會限制其在實(shí)時成像中的應(yīng)用。通過將SIRT 與有限域FBP 相結(jié)合，則可在保證重建精度的基礎(chǔ)上，顯著減少重建時間。其原理是：首先利用有限域FBP 重建初始圖像，然后使用SIRT 對初始圖像進(jìn)行優(yōu)化[14]，結(jié)果如圖6所示。另外，對于換能器陣列接收立體角有限的問題，則可將換能器陣列的方向性模式函數(shù)作為迭代過程中的加權(quán)因子，以進(jìn)一步改善圖像質(zhì)量。

圖6 與有限域FBP 相結(jié)合的SIRT 重建結(jié)果[14]Fig.6 Results of SIRT combined with FBP in limited-view[14]

(2) 組合迭代重建算法

圖7 3 種算法的重建結(jié)果對比Fig.7 Comparison of reconstruction results of three algorithms[15]

(3) 迭代加權(quán)FBP 算法

迭代自適應(yīng)加權(quán) FBP(Iterative Adaptive Weighted FBP, IAWFBP)[16]的原理是：首先，利用基于有限角度的FBP 重建初始圖像；然后，根據(jù)初始圖像和光聲信號測量值計(jì)算投影信號的加權(quán)系數(shù)；其次，利用初始圖像和加權(quán)系數(shù)得到模擬的光聲信號，并計(jì)算其與測量值之間的殘差；最后，對殘差信號做反投影運(yùn)算得到殘差圖像，并將其與上一次迭代生成的光吸收分布圖相加獲得修正圖像。重復(fù)上述過程，并設(shè)置最大迭代次數(shù)和迭代誤差作為迭代結(jié)束條件。FBP、ART-FBP 和IAWFBP 算法的性能對比如圖8 所示。IAWFBP 方法可以有效減少圖像中的偽影，且在精度、對比度和對噪聲的魯棒性等方面都有很好的表現(xiàn)。但是，該算法存在成像系統(tǒng)掃描時間過長的缺陷，可以使用多陣元探測器陣列和多光譜高頻激光源減少掃描時間。

采用正則化迭代加權(quán)FBP(Regularized Iterative Weighted FBP, RIWFBP)[17]還可進(jìn)一步提高重建精度。首先，使用FBP 從有限角度測量數(shù)據(jù)中重建初始光吸收分布圖；然后，在每一次迭代中，計(jì)算光聲信號的測量值和模擬值之間的差值，以此修正圖像，并采用正則化方法提高迭代的收斂速度。實(shí)驗(yàn)證明，采用該算法在前五次迭代中即可有效抑制IAWFBP 中存在的有限角度掃描所致的偽影。而且該方法在精度和對噪聲的魯棒性方面均優(yōu)于IAWFBP 和傳統(tǒng)的迭代方法，能夠顯著減少所需的超聲換能器數(shù)量和掃描時間，IFBP、IAWFBP、RIWFBP 算法的性能對比如圖9 所示。

圖8 FBP、ART-FBP 和IAWFBP 算法的性能對比[16]Fig.8 Comparison of FBP, ART-FBP and IAWFBP[16]

(4) 濾波平均反投影迭代算法

FBP 是將去除探測器脈沖響應(yīng)的光聲信號測量值直接反投影到對采集信號有貢獻(xiàn)的測量點(diǎn)上，會造成不同位置測量點(diǎn)的投影不均勻現(xiàn)象。利用濾波平均反投影算法可解決該問題[18]：首先，對光聲信號測量值進(jìn)行反卷積去除探測器脈沖響應(yīng)的影響，統(tǒng)計(jì)對采集信號有貢獻(xiàn)的投影點(diǎn)個數(shù)，求得實(shí)際光聲信號均值，并將平均后的信號對每一個測量點(diǎn)進(jìn)行反投影重建；然后，對成像區(qū)域中每一點(diǎn)的像素值用該點(diǎn)的投影累加次數(shù)進(jìn)行平均。采用該方法重建的初始聲壓分布圖的結(jié)構(gòu)與實(shí)際的聲場分布結(jié)構(gòu)大致相同，各點(diǎn)的聲壓幅值和實(shí)際值較接近。

圖9 IFBP、IAWFBP 和RIWFBP 性能對比[17]Fig.9 Comparison of performance of IFBP, IAWFBP and RIWFBP[17]

在傳統(tǒng)的迭代重建中，聲壓分布的初始值通常采用空矩陣或者直接反投影的結(jié)果，可能增加迭代次數(shù)，延長計(jì)算時間。濾波平均反投影迭代(Filtered Mean-Back-Projection Iterative, FMBPI)算法[19]能夠彌補(bǔ)傳統(tǒng)迭代算法在計(jì)算時間上的不足，其原理是：首先，利用濾波平均反投影算法從不完備的光聲數(shù)據(jù)中重建初始聲壓分布，作為迭代的初始值；然后，采用迭代重建算法修正初始聲壓分布，減小光聲信號的計(jì)算值和測量值之間的誤差，去除由反投影引入的偽影，F(xiàn)BP 和FMBPI 結(jié)果對比如圖10所示。FMBPI 算法的適應(yīng)性較強(qiáng)，由于采用濾波均值反投影算法對光聲信號進(jìn)行預(yù)處理，因而可以消除反投影的不均勻現(xiàn)象。另外，相對于傳統(tǒng)的迭代重建算法，F(xiàn)MBPI 利用初始聲壓分布與實(shí)際分布之間的逼近關(guān)系以及權(quán)重因子來減少迭代次數(shù)，可以在相對較短的時間內(nèi)重建出精度更高的光聲圖像。

牛肚菌。牛肚菌生長于海拔900－2200米的松樹混交林中或砍伐不久的邊緣地帶，多產(chǎn)于云南省。牛肚菌香甜可口、營養(yǎng)豐富，中醫(yī)認(rèn)為其具有養(yǎng)血和中、祛風(fēng)散寒、舒筋活血等功效，是婦科良藥，同時還有抗流感、治感冒的作用。

2.3.3 基于梯度投影的迭代算法

基于稀疏探測器陣列和迭代重建算法的PAT成像系統(tǒng)為高幀率三維PAT 成像的開發(fā)提供了可能性，應(yīng)用基于梯度投影的迭代圖像重建算法[20]可以提高成像效率，獲取以單束激光作為光源、幀率為10 Hz 的三維PAT 圖像(圖11)。該算法的原理是：首先，設(shè)定解的初始值；然后，由約束條件確定初始值的凸約束集邊界上梯度的投影，求出可行的下降方向和步長；接著，線性搜索最優(yōu)步長，并以此為條件求得新的迭代點(diǎn)；最后，進(jìn)行收斂測試，檢驗(yàn)是否可以終止迭代。

圖10 FBP 和FMBPI 結(jié)果對比[19]Fig.10 Comparison of reconstruction results with FBP and FMBPI[19]

圖11 SIRT 和基于梯度投影的迭代算法的重建結(jié)果對比[20]Fig.11 Comparison of reconstruction results with SIRT and gradient projection-based iterative algorithm[20]

2.4 基于Gerchberg-Papoulis 外推的補(bǔ)償方法

由信號f(t)的部分已知數(shù)據(jù)g(t)重建f(t)在整個定義域上的值稱為帶限函數(shù)外推問題。Gerchberg-Papoulis(GP)外推算法是解決帶限函數(shù)外推問題的一種迭代算法[21]，其理論依據(jù)是長球面波函數(shù)展開及其相關(guān)的譜估計(jì)理論[22]。用GP 外推算法可解決有限角度PAT 成像中由不完備采樣數(shù)據(jù)造成的盲區(qū)、模糊和偽影等問題[8]。首先，利用快速傅里葉變換得到光聲信號測量值的頻譜，在頻域中做反投影運(yùn)算，再利用傅里葉逆變換得到初始聲壓分布，進(jìn)而求出空間光吸收分布，以此作為迭代的初始值；然后，利用GP 外推算法對盲區(qū)信號進(jìn)行補(bǔ)償，以修正圖像，同時結(jié)合全變分(Total Variation, TV)正則化提高迭代的收斂速度。與迭代重建再投影(Iterative Reconstruction Reprojection, IRR)算法相比，該算法可以加快補(bǔ)償速度，減少計(jì)算時間，解決稀疏探測器陣列采集信號所導(dǎo)致的偽影、失真和模糊問題，增強(qiáng)圖像邊緣，更好地保留細(xì)節(jié)信息，3 種算法的性能對比如圖12 所示。

2.5 基于壓縮感知的方法

2.5.1 壓縮感知重建算法

壓縮感知(Compressed Sensing, CS)理論打破了奈奎斯特抽樣定理的限制，認(rèn)為只要信號是稀疏的或可壓縮的，就能夠通過數(shù)據(jù)壓縮節(jié)省對大量無用信息的采樣，然后求解特定的優(yōu)化問題，即可利用少量的壓縮數(shù)據(jù)恢復(fù)原信號。

實(shí)際人體組織中，多數(shù)光吸收體的時空相關(guān)性較低，所以常規(guī)的采樣數(shù)據(jù)通常是高度冗余的。由于大部分光聲圖像具有稀疏性，或者針對某個變換基函數(shù)滿足稀疏性條件，因此利用少量隨機(jī)采樣信號能夠?qū)崿F(xiàn)圖像的重建。即把原始光聲信號在特定測量矩陣上的投影作為測量數(shù)據(jù)，通過稀疏約束重構(gòu)算法恢復(fù)原始信號，最終利用少量數(shù)據(jù)和低成本硬件實(shí)現(xiàn)高分辨率光聲成像[23]。

圖12 TV-GD, TV-VB, TV-GPEF 算法的性能對比[8]Fig.12 Comparison of performance for TV-GD, TV-VB,TV-GPEF[8]

例如：在文獻(xiàn)[24]中提出了一種基于CS 的光吸收分布圖重建方法，利用全采樣光聲數(shù)據(jù)的1/4即可重建光吸收分布圖。其原理是：初始光吸收分布通常由組織邊界的光滑部分和跳變點(diǎn)構(gòu)成，所以其拉普拉斯變換是稀疏的(或者至少是被壓縮的)。根據(jù)光聲信號測量值的二階時間導(dǎo)數(shù)與初始光吸收分布的拉普拉斯變換之間的對應(yīng)關(guān)系，重建出光吸收分布圖[25]。針對脈沖激光隨機(jī)照射條件下的PAT 成像，文獻(xiàn)[25]提出利用CS 算法減少隨機(jī)照射次數(shù)來提高數(shù)據(jù)采集速度，僅利用兩個夾角為90°的探測器采集的信號測量值就能精確地重建光聲圖像。

此外，激光脈沖的重復(fù)頻率和超聲換能器的數(shù)量會限制光聲信號的采集速度。通過采用線性集成探測器進(jìn)行多角度掃描，設(shè)置合理的探頭分布角度以及測量矩陣等，可以提高信號的采集和處理速度，有效彌補(bǔ)遠(yuǎn)場成像分辨率，消除圖像偽影[26]。而CS 重建算法可以壓縮采樣數(shù)據(jù)，減少超聲換能器的數(shù)目，加速數(shù)據(jù)采集，降低硬件成本[27]，其圖像重建效果優(yōu)于反投影(Back-Projection, BP)和迭代重建[27-28]，如圖13 所示。另外，在提高對比噪聲比(contrast-to-noise ratio, CNR)和重建速度等方面，頻域CS 比時域CS 更具優(yōu)勢[29]。

圖13 反投影、迭代重建和CS 重建算法的性能對比[27]Fig.13 Comparison of performance of BP,IR and CS[27]

2.5.2 改進(jìn)的壓縮感知重建算法

(1) 結(jié)合快速交替方向法的壓縮感知算法

在使用CS 算法重建圖像的過程中，利用快速交替方向法(Fast Alternating Direction Method,FADM)可以很好地求解L0-正則化問題。首先，利用變量分裂技術(shù)將稀疏信號的L0-正則化問題轉(zhuǎn)化為約束優(yōu)化問題；然后，采用一步Gauss-Seidel 的思想，利用乘子函數(shù)對優(yōu)化問題中的變量做極小化處理；其次，二次更新變量和乘子；最后，利用反正交變換重建原始信號。與FBP 和其它CS 算法相比，F(xiàn)ADM 可以更快地減小相對誤差，在計(jì)算效率和數(shù)據(jù)保真方面效果理想。此外，該算法對噪聲有很好的魯棒性[31-32]，F(xiàn)BP 算法與FADM 重建結(jié)果的對比如圖14 所示。

圖14 FBP 算法與FADM 重建結(jié)果的對比[30]Fig.14 Comparison of reconstruction results with FBP and FADM[30]

(2) 具有部分已知支集的壓縮感知算法

具有部分已知支集的 CS(CS with Partially Known Support, CS-PKS)算法[32-35]是將部分已知信號的空間位置作為CS 重建中的先驗(yàn)信息，降低待重建信號的復(fù)雜度。與傳統(tǒng)CS 算法相比，CS-PKS可以用稀疏欠采樣數(shù)據(jù)的1/3 重建具有更高CNR 的圖像，減少高保真重建所需的測量次數(shù)，并能夠以相對較少的迭代次數(shù)快速收斂，傳統(tǒng)CS 方法與CS-PSK 方法對人手血管的成像結(jié)果如圖15 所示。

(3) 合標(biāo)準(zhǔn)重建方法的壓縮感知算法

圖15 人手血管的PAT 圖像[33]Fig.15 PAT images of the subcutaneous vasculature of human hand[33]

目前的三維PAT 成像系統(tǒng)能夠?qū)崿F(xiàn)高幀率的成像，但不能同時提供高空間和時間分辨率的圖像，這限制了其動態(tài)成像(即四維PAT)的能力。結(jié)合標(biāo)準(zhǔn)重建法的CS 算法為解決這一問題提供了可能[25,37]。原理是：采用CS 算法從空間欠采樣的光聲信號中恢復(fù)完整的光聲信號，然后利用標(biāo)準(zhǔn)重建算法重建圖像。此外，還可通過結(jié)合TV 正則化方法，加速迭代收斂。這種復(fù)合算法能夠減少空間測量次數(shù)，降低算法的復(fù)雜度，重建出具有較好空間分辨率和對比度的圖像，反投影法和結(jié)合反投影的CS 重建結(jié)果如圖16 所示[37]。

圖16 反投影法和結(jié)合反投影的CS 重建結(jié)果[36]Fig.16 Reconstruction results with BP and CS combined with BP[36]

2.6 基于主成分分析的方法

主成分分析(Principal Component Analysis,PCA)是一種常用的統(tǒng)計(jì)方法，可以將高維數(shù)據(jù)空間轉(zhuǎn)化為低維特征空間，在去除冗余、特征提取以及數(shù)據(jù)壓縮等方面效果顯著[37]，已被廣泛應(yīng)用于圖像處理、基因數(shù)據(jù)挖掘和醫(yī)學(xué)成像系統(tǒng)中[38-41]。

基于PCA 的三維光聲層析(PCA-PAT)成像[42]使用部分全采樣數(shù)據(jù)來計(jì)算全部數(shù)據(jù)的主要成分。將一些均勻分布的全采樣圖像作為PCA 的訓(xùn)練樣本，通過導(dǎo)出的PCA 數(shù)據(jù)恢復(fù)其它稀疏采樣的圖像，無需迭代過程便可根據(jù)稀疏采樣數(shù)據(jù)快速重建高質(zhì)量的三維光聲圖像。

PCA-PAT 成像是一種低成本的快速三維成像方法，可以有效地降低數(shù)據(jù)采集和圖像重建的時間，與BP 相比，測量次數(shù)可減少約50%。在重建時間縮短約40%的情況下，可重建出與基于全采樣的BP 質(zhì)量相同的圖像。另外，當(dāng)超聲換能器數(shù)量相同時，PCA 的圖像重建速度約是CS 的8 倍，且成像精度更高，分別如表3 和圖17 所示。

表3 BP，CS 和PCA 算法重建時間的比較[42]Table 3 Comparison of reconstruction time of BP, CS and PCA methods [42]

圖17 應(yīng)用BP、CS 和PCA 的重建圖像[42]Fig.17 PAT images reconstructed with BP, CS and PCA methods[42]

2.7 基于深度學(xué)習(xí)的方法

深度學(xué)習(xí)(deep learning)是機(jī)器學(xué)習(xí)中對數(shù)據(jù)進(jìn)行表征學(xué)習(xí)的一種方法，在模式識別和機(jī)器學(xué)習(xí)領(lǐng)域具有廣泛的應(yīng)用。近年來，深度學(xué)習(xí)算法特別是卷積神經(jīng)網(wǎng)絡(luò)(Convolutional Neural Network,CNN)已經(jīng)迅速成為醫(yī)學(xué)圖像分析的首選方法[43]，它在快速重建高質(zhì)量的圖像方面也具有巨大潛力。針對稀疏采樣的PAT 成像，通過將CNN 與標(biāo)準(zhǔn)重建方法相結(jié)合，可以得到高分辨率的光聲圖像[44]。與經(jīng)典的迭代重建算法相比，基于CNN 的重建算法在速度和質(zhì)量等方面也具有較大優(yōu)勢。

基于CNN 的PAT 圖像重建算法可以分為兩類：(1) 首先采用快速、簡單的標(biāo)準(zhǔn)重建算法從稀疏采樣的光聲數(shù)據(jù)中重建出含有偽影的低質(zhì)量圖像，然后用訓(xùn)練后的CNN 進(jìn)行圖像后處理，去除偽影[45-46]，重建結(jié)果如圖18 所示；(2) 基于模型的深度學(xué)習(xí)，即將PAT 成像的前向物理模型包含到CNN 的訓(xùn)練和圖像重建中。例如文獻(xiàn)[44]提出的方法將稀疏光聲數(shù)據(jù)的梯度信息和迭代算法相結(jié)合，首先計(jì)算光聲信號測量值的梯度，通過訓(xùn)練使CNN 學(xué)習(xí)目標(biāo)圖像的先驗(yàn)知識(如結(jié)構(gòu)特點(diǎn)等)；然后，將訓(xùn)練后的CNN 應(yīng)用于光聲測量數(shù)據(jù)，并使用迭代算法重建光聲圖像，重建結(jié)果如圖19 所示。

相對于點(diǎn)探測器，采用集成線性探測器可實(shí)現(xiàn)對三維結(jié)構(gòu)的高分辨PAT 成像。最近，使用64 個并行線性探測器陣列的三維PAT 成像系統(tǒng)也已實(shí)現(xiàn)[47]。針對此類系統(tǒng)，將動態(tài)孔徑長度(Dynamic Aperture Length, DAL)校正算法與CNN 方法相結(jié)合，能夠?qū)崟r生成高分辨率的三維投影圖像[48]。

圖18 CNN 后處理的重建結(jié)果[45]Fig.18 Reconstruction results of CNN post-processing[45]

圖19 人手掌的光聲圖像[44]Fig.19 The photoacoustic tomography of a human palm[44]

3 有限角度定量光聲層析成像

對于有限角度掃描的多光源 qPAT (Multiple-Source, MS-qPAT)，不完備的測量數(shù)據(jù)會導(dǎo)致無法穩(wěn)定重建初始聲壓分布，故傳統(tǒng)的兩步算法[4]不再適用。此時可采用一步算法[49]，即利用多個光源照射得到的多組聲壓數(shù)據(jù)直接重建光吸收系數(shù)和散射系數(shù)。例如，文獻(xiàn)[50-51]提出一種實(shí)用的MS-qPAT 圖像重建方法，首先確定超聲換能器的空間沖激響應(yīng)(Spatial Impulse Response, SIR)和聲電沖激響應(yīng)(Acousto-Electric Impulse Response, EIR)，模擬聲學(xué)測量；然后，采集光源附近區(qū)域的聲壓信號，在光聲耦合前向傳輸模型中結(jié)合空間SIR、EIR以及有限角度掃描，利用基于交替方向法(Alternating Direction Method, ADM)的擬牛頓法獲得光吸收系數(shù)和散射系數(shù)的精確解。

4 結(jié) 語

PAT 成像的效率和圖像質(zhì)量與圖像重建算法有直接的關(guān)系。目前，多數(shù)成熟的標(biāo)準(zhǔn)重建算法都是基于完備數(shù)據(jù)的，但在實(shí)際應(yīng)用中受探測器陣元尺寸、機(jī)械結(jié)構(gòu)、接收角度以及帶寬等的限制，往往不能滿足全角度掃描和完備數(shù)據(jù)采集的條件。若仍采用標(biāo)準(zhǔn)重建算法，則可能導(dǎo)致圖像中的偽影嚴(yán)重，空間分辨率降低。近些年，針對有限角度掃描和稀疏測量數(shù)據(jù)的問題，各種新的圖像重建算法層出不窮。這些算法能夠利用高度不完備的測量數(shù)據(jù)恢復(fù)高質(zhì)量的圖像，它們對掃描角度和數(shù)據(jù)稀疏度的要求如表4 和表5 中所示。

減少測量次數(shù)、降低系統(tǒng)成本、提高數(shù)據(jù)采集速度和成像速度是今后PAT 研究的主要發(fā)展方向，這對于實(shí)現(xiàn)三維光聲圖像重建以及高性能、高分辨率和高性價比的PAT 成像系統(tǒng)開發(fā)具有重要意義。

表4 主要算法的有效掃描角度范圍Table 4 Effective angle range of the algorithms mentioned in the paper

圖5 主要算法的數(shù)據(jù)稀疏情況Table 5 The sparse data of the algorithms mentioned in this paper