生物組織聲學特性的空間變化通常是未知的，不同成分的組織往往具有不同的聲學特性。對于光聲層析成像（PAT），恒定聲速和聲學特性均勻分布的假設會導致重建圖像的細節模糊、目標錯位及圖像偽影等問題。本文對生物 PAT 圖像重建中組織聲學特性不均勻性問題（主要是聲速不均勻和聲衰減）的主要解決方法進行總結和歸納，討論各方法的優勢和不足，并展望未來可能的發展方向。
引用本文： 段爽, 孫正. 聲學特性不均勻組織的光聲層析圖像重建研究進展. 生物醫學工程學雜志, 2019, 36(3): 486-492. doi: 10.7507/1001-5515.201809014 復制
版權信息： ?四川大學華西醫院華西期刊社《生物醫學工程學雜志》版權所有，未經授權不得轉載、改編
引言
生物光聲層析成像（photoacoustic tomography，PAT）是一種多物理場耦合的新型復合功能成像方法，對疾病的早期檢測和準確診斷具有極高的研究價值[1]。其成像參數是組織的光吸收系數和散射系數，物理基礎是生物組織的光聲效應[2]。如圖 1 所示，用短脈沖激光照射組織，組織吸收激光能量產生熱彈性膨脹，并以超聲波的形式向表面傳播，即光聲信號[1]。超聲換能器檢測到光聲信號后，采用適當的圖像重建算法即可反演計算出組織內部的光吸收分布或初始聲壓分布[3]，而組織的光吸收特性與其血紅蛋白濃度和分子組成密切相關[4]。PAT 既具有光學成像的高對比度，又具有超聲成像的高穿透深度[5]，在血管成像[6]、血液流速檢測[7]、腫瘤檢測[8]及腦結構和功能成像[9]等領域具有重要意義。
圖1
PAT 成像原理示意圖
Figure1.
Schematics of PAT imaging
有效的光聲信號處理方法和高效的圖像重建算法是保證成像精度和速度的關鍵。目前較為成熟的 PAT 圖像重建算法有濾波反投影（filtered back-projection，FBP）、時間反演（time reversal，TR）、迭代重建法、相控聚焦法和基于傅里葉變換的重建算法等[10-11]。相對于純光學成像，PAT 的高成像深度和高分辨率源于超聲波在生物組織內傳播時的低散射。因此為了簡化問題，多數光聲圖像重建算法通常都假設組織的密度和超聲波在組織中的傳播速度是恒定的，不考慮聲散射，且成像組織是聲學特性均勻的非衰減聲學介質。
然而在實際應用中，考慮到人體組織的復雜性，不同成分的組織通常具有不同的聲學特性，而且其空間變化也是未知的。例如：軟組織內部的聲速差異可以達到 10%（皮下脂肪中的聲速是 1 400～1 450 m/s，而正常薄壁組織和基質組織中的聲速是 1 500～1 560 m/s[12]）。聲速差異會影響光聲信號到達超聲換能器的時間，進而影響根據傳播時間計算出的聲源到超聲換能器的距離。而恒定聲速的假設會導致重建圖像中目標組織的位置與其實際位置之間存在誤差，進而使圖像中存在細節模糊、目標錯位及嚴重偽影等問題[13]。此外，成像組織聲學特性的不均勻性會造成超聲波在組織中傳播時存在反射、散射、擴散和吸收等現象[14]，使部分光聲信號偏離原來的傳播方向，進而導致光聲信號的衰減，并且頻率越高衰減越嚴重，特別是在聲學特性不匹配的組織邊界處，如骨骼、肺、空腔等。此時，若假設成像組織的聲學特性均勻，則會導致重建圖像中包含嚴重的失真和偽影，降低圖像質量，同時可能降低成像深度。因此組織聲學特性不均勻問題是 PAT 圖像重建中需要解決的關鍵問題，對提高成像質量、減少圖像偽影和失真等具有重要意義。
本文對聲學特性不均勻（主要是聲速不均勻和聲衰減）組織的 PAT 圖像重建的主要方法進行綜述，討論各方法的優勢和不足，并展望未來可能的發展方向。
1 聲速不均勻的解決方法
在進行 PAT 成像前，采用超聲透射斷層掃描（ultrasonic transmission tomography，UTT）可獲取成像組織的聲速分布，但需要增加額外的成像設備和成像過程，限制了其應用范圍。而將有關聲速分布的先驗知識結合到圖像重建中，補償由于聲速變化引起的測量光聲場的像差，可以有效改善圖像的聚焦效果，提高圖像質量[13, 15-17]。常用的方法有：有限元法（finite element method，FEM）[18-19]、基于最優聚焦的平均聲速估計法[20]、全波迭代法[4, 21-23]、修正的時域重建法和修正的反卷積重建法[24]等，各方法的性能如表 1 所列。
表1
解決 PAT 成像中聲速不均勻方法的性能比較
Table1.
Performance of methods for PAT in tissues with variable speed of sound
1.1 有限元法
運用有限元方法將光聲波動方程所在的連續場分割成有限個子空間區域，然后應用插值函數表示每個子空間區域的解，將一個連續的無限自由度問題轉化成離散的有限自由度問題進行求解，可以只利用光聲信號同時重建出組織的光吸收分布和聲速分布的低頻成分。
該方法的重建精度高，求解域可以是任意形狀，可得到較好的目標位置信息[25]。缺點是計算量非常大，計算速度慢，且只能重建出待測組織圖像中的低頻成分[18]。
1.2 最優聚焦法
可將聚焦光聲信號絕對差的平均和作為聚焦品質因子 Q，評估光聲信號是否具有最優聚焦，其定義為：
 |
其中，M 和 N 分別是焦點和超聲探測器的數量，p 和 
分別是聚焦光聲信號的聲壓及其平均值。最優聚焦法的基本思想是：改變聲速，同時針對感興趣區域（region of interest，ROI）中的每條掃描線迭代計算 Q，使 Q 取得最小值的聲速即為 ROI 的最佳平均聲速，再采用標準重建法重建初始聲壓分布圖。
該方法可以提高重建圖像中目標的對比度和空間分辨率，計算簡便。局限性是需要有關目標粗略位置的先驗信息，而且如果 ROI 過小，可能增加迭代計算時間。
1.3 全波迭代重建方法
全波方法用于聲學非均勻有損介質的迭代 PAT 圖像重建，基本思想是建立基于精確光聲波動方程的離散成像模型[4]，利用 k-空間偽譜法計算出前向投影算子和后向投影算子，進而得到離散成像算子（即系統矩陣）及其伴隨或反投影算子，再結合正則化方法［如 TV（total variation）正則化］和迭代算法［如快速迭代收縮閾值算法（fast iterative shrinkage thresholding algorithm，FISTA）］求解偏最小二乘問題，重建初始聲壓分布圖。
該方法能夠有效地補償聲像差，減少重建圖像中的偽影和失真，與 TR 法相比有較高的分辨率。利用 k-空間偽譜法求解耦合方程，采用快速傅立葉變換法計算空間偏導數，比有限差分法具有更高的計算效率，且空間和時間采樣要求較少。局限性是圖像重建精度易受光聲信號中噪聲的影響，且計算時一般需要較多的迭代次數。
1.4 修正的時域重建法和修正的反卷積重建法
修正的時域重建法和修正的反卷積重建法都是在球面掃描條件下，利用在不同位置采集的光聲信號之間的相關性估計超聲波在空間兩點之間的傳播時間，并據此分別修正傳統的時域重建法和反卷積重建法[24]，補償聲速的不均勻性。兩種方法均無需預知聲速分布即可直接重建出目標圖像，且運算速度快。
修正的時域重建法的基本思想是：首先，依據在關于原點對稱的位置處探測到的光聲信號之間相關性最強的原理，估算超聲波在空間兩點之間的傳播時間；其次，利用比例函數關系，求得組織中任意位置處產生的光聲信號到達超聲換能器的時間；最后，根據傳播時間計算光吸收分布。與標準重建方法相比，采用該方法重建的圖像精度更高，邊緣更清晰，且具有良好的抗噪性能，對隨機噪聲不敏感。不足之處是存在理想化的假設，如忽略介質分界面對聲波的折射因素，認為光聲信號從組織內到超聲換能器是沿直線傳播的，導致估計的超聲波傳播時間存在一定誤差，而且聲速不均勻性越強或者待測組織越大，誤差也越大。同時，在計算超聲波傳播時間時，近似認為積分面是一個平面，可能降低重建圖像的質量。
修正的反卷積重建法的原理是：首先，構造一個函數 C（r）：
 |
其中，
，
是常數且
，r0 是超聲換能器的位置。采用修正的時域重建法估計組織中任意位置 r 處產生的光聲信號到達超聲換能器的時間 T（r, r0），并將其代入式（2）進行修正：
 |
其次，將組織的光吸收分布 A（r）和 C（r）分別看作一個線性時不變系統的輸入和輸出，即：
 |
其中 h（r）是系統的單位沖激響應。當輸入為 δ 函數時，即可得到 h（r）的表達式。則有：
 |
其中
、
和
分別是 A（r）、h（r）和 C（r）的傅里葉變換，對
進行逆傅里葉變換即可得到光吸收分布。當組織的聲速差異低于 10% 時，采用該方法重建的圖像精度較高且抗噪性能較好。通過使用快速傅里葉變換算法，可獲得極高的計算速度，并且適用于有限角度掃描的情況[24]。但它是一種近似算法，而且待測組織的尺寸越大，近似誤差也越大。
2 聲衰減的補償方法
超聲波在組織中傳播時會發生衰減，并且頻率越高衰減越嚴重。造成衰減的主要因素有超聲波的反射、散射、擴散和吸收等。其中，隨機散射對超聲波的傳播有很大影響，目前解決 PAT 中聲散射問題的主要方法有 Green 函數重構法[14]和基于聲散射體分布的統計加權方法[26-30]等。此外，補償聲衰減的主要方法有矯正不同頻率成分衰減的方法[31]、時變濾波法（time-variant filtering，TVF）[32]、迭代正則化法[33]、基于奇異值分解（singular-value decomposition，SVD）的方法[34-36]和基于 TR 的方法[37-41]等。上述方法優缺點的對比如表 2 所列。
表2
PAT 成像中聲衰減補償方法的優缺點對比
Table2.
Comparison of methods for compensating acoustic attenuation in PAT
2.1 Green 函數重構法
根據互易定理推導出的、適用于常見散射介質的 Green 函數有其局限性。對于采用超聲換能器陣列的散射介質 PAT 成像，可將成像區域劃分為若干網格，根據互易定理在待測點構建一個虛擬聲源，計算該點到每一個換能器的 Green 函數，再結合 TR 方法，得到組織中的初始聲壓分布，進而得到光吸收分布。
采用該方法重構組織內某點的 Green 函數時，無需在該點放置超聲換能器或實際聲源，而是虛構一個虛擬聲源，快速計算該點與散射域另一側測量點之間的 Green 函數，具有計算速度快的優勢，適用于實時 PAT 成像。但是當探測信號頻率降低時，重構的 Green 函數精度也明顯下降，且當有隨機噪聲擾動時，使用理想波形濾波器無法重構出 Green 函數。
2.2 基于聲散射體分布的統計加權方法
在成像目標和超聲探測器之間放置強吸收體，測量到的超聲信號首先用于估計成像目標中散射體的空間分布。然后，在此基礎上計算光聲信號不受反射或散射聲波影響的概率，并利用該概率函數對光聲信號測量值進行加權，從而利用強吸收體產生的、以更高概率從聲源直接傳播到探測器的光聲信號重建光吸收分布，校正由于聲波反射或散射引起的圖像偽影[26-27]。由于確定這種概率的精度取決于被測組織的結構信息，因此可以將有關高散射區近似位置的先驗知識引入重建過程，從而進一步提高層析圖像重建的質量[28-30]。
此類方法能夠去除重建圖像中由于聲散射引起的多數偽影，與假設聲散射體隨機分布的方法相比，能夠顯著提高強散射或反射條件下的光聲成像質量。但是重建圖像的精度受估計的聲散射體分布信息的直接影響，而且需要利用除成像目標之外的強吸收體，可能增加成像復雜度。
2.3 矯正不同頻率成分的聲衰減
根據超聲波衰減理論，由于介質的吸收和散射，光聲信號中的不同頻率成分隨傳播距離的衰減是不同的，高頻成分的衰減較多，低頻成分的衰減較少，并且不同組織對各種頻率超聲波的衰減也是不同的，因此，可對光聲信號不同頻率成分的衰減進行矯正。思路是：組織表面的光聲壓與光吸收分布滿足光聲波動方程，而實際測量的光聲信號是光聲壓與探測器脈沖響應的卷積，因此可以通過頻域的快速傅里葉變換得到光吸收分布與光聲信號測量值之間的關系。再結合超聲衰減公式和組織的衰減參數得到對光聲信號不同頻率成分進行矯正后的光聲圖像重建公式。此方法的計算速度快，可補償光聲信號中各種頻率成分隨距離的衰減，增強信號的高頻成分，突出吸收體的邊界變化和細微結構特征，提高成像分辨率。
2.4 時變濾波法
在假定組織的聲衰減參數分布均勻的前提下，對采集的每個聲壓時間序列進行頻域的時變濾波，修正聲衰減和色散。構造時變濾波器時，可利用 Tukey 窗口進行正則化，并根據聲壓測量值的時頻分布選擇濾波器的截止頻率，與設定固定截止頻率相比，可有效減少引入的高頻噪聲。
該方法適用于任何掃描幾何，且計算效率非常高。局限性在于聲衰減參數在介質內分布均勻的假設可能會引入一定的誤差，且用于正則化的時變窗口的性能也可再進一步優化，以改善濾波結果。
2.5 迭代正則化方法
迭代正則化方法是在基于最小二乘原理的廣義逆算法的基礎上改進而來的一種穩定的廣義逆算法[42]。用于補償聲衰減的迭代正則化 PAT 圖像重建方法是將聲衰減介質的 PAT 逆問題描述為 Hilbert 空間中的一個不適定算子方程[33]，并用迭代正則化法進行求解。其基本思路是：首先，根據衰減定理建立超聲波在衰減介質中傳播的一般模型，并在 Hilbert 空間中設計 PAT 成像的前向算子；然后，選擇適當的正則化方法計算前向算子的伴隨算子，同時對最小二乘問題施加正則化約束；最后，利用迭代法求解最小二乘問題，反演重建出初始聲壓分布圖。
選擇適當的正則化方法是此類方法的關鍵，例如 Landweber 正則化法主要用于計算伴隨算子，由于伴隨項的平滑效應，因而其性能非常穩定[33]。迭代 Tikhonov 正則化在一定條件下能夠達到最優收斂速度，但不適合大規模計算[43]。此外，它需要根據光滑性條件選取最小階數，而不必提高迭代次數。
迭代正則化法是一種通用的反演方法，可以應用于不同的衰減模型和測量幾何，以及稀疏測量數據的 PAT 圖像重建，對于弱衰減的情形，其重建結果甚至優于無衰減的情況。但是在某些應用中它的計算速度很慢，需要較多的迭代次數。
2.6 基于奇異值分解的方法
組織對超聲波的吸收與其頻率直接相關，在時域-頻域中，包含衰減項的光聲波動方程是具有復數波數的非齊次 Helmholtz 方程，直接求解的計算量非常大。而通過對組織中的聲衰減進行建模，可在圖像重建中對其進行校正。
例如，文獻[34]通過積分算子 L 將理想的、未衰減的聲壓信號 p 與在組織表面測量到的衰減的聲壓信號 patt 相關聯，即
 |
算子 L 可根據遵循冪律吸收（power law absorption）的有損波動方程在時域-頻域中關聯衰減與未衰減超聲信號得到。對矩陣方程的特征值進行閾值化處理，采用 SVD 計算L的偽逆算子，進而從光聲信號測量值中估計理想的未衰減信號。最后，使用標準重建方法從未衰減信號中重建出組織的光吸收分布圖，實現對聲衰減的補償。但是由于衰減與未衰減信號之間的積分方程一般沒有很好的約束條件，因此這種反演通常是病態的。
為了解決上述問題，可根據光聲信號測量值的傅里葉變換與無衰減理想信號的傅里葉變換得到關聯方程，在對方程離散化得到矩陣方程之后，采用截斷奇異值分解（truncated singular value decomposition，TSVD）計算逆矩陣，進而從衰減的測量信號中獲得對理想信號的近似，實現一步補償[34]，消除聲衰減對光聲信號測量值的影響，再應用標準重建方法（如 TR 等）重建光聲圖像。
此外，文獻[36]提出一種將聲壓信號測量值 patt 直接與重建的初始聲壓信號 p0（包括衰減補償）相關聯的方法，關聯矩陣是時間、衰減系數和二維空間頻率的函數。采用 TSVD 計算矩陣方程的正則化偽逆，進而重建初始聲壓分布圖和光吸收分布圖。該方法適用于平面探測器幾何結構，能夠表征重構衰減介質中不同深度處光吸收分布的不適定性。
2.7 基于 TR 的方法
大孔徑超聲探測器的成像分辨率受到超聲波帶寬和衰減的限制，而且頻率越高衰減越嚴重，如果忽略這種影響，會導致重建圖像中小目標的分辨率嚴重下降。針對此問題，文獻[37-39]中設計了一種簡單的裝置，利用一個旋轉軸將平面探測器分隔成垂直于旋轉軸的線探測器。然后，利用在固定旋轉角度掃描的一組光聲信號，使用二維 TR 方法重建初始聲壓的線性投影。最后，計算與旋轉軸正交的平面中的二維逆 Radon 變換，精確重建成像目標的三維圖像。該方法可以精確地解決任意封閉檢測面的光聲逆問題，且能部分補償聲衰減對圖像重建的影響。其局限性在于有關聲衰減的兩個波動方程的時間反演是不穩定的，不能完全補償聲衰減的影響，重建圖像中存在一定的模糊。
此外，還可將 TR 重建方法與有損波動方程相結合，即反轉斯托克斯有損波動方程中吸收項的符號，那么聲波在組織中的傳播就是增強而不是衰減的。對高于預設截止頻率的信號進行濾波，使重建過程正則化，避免含噪聲的高頻信號增長過快。最后對聲學正問題的數值模型進行時間反演，在每個時間步長上，將測量的時變聲壓信號以時間反轉的順序作為 Dirichlet 邊界條件。當時間 t 達到 0 時，即可重建出初始聲壓分布。
將 TR 圖像重建和吸收狀態方程相結合，使用基于分數拉普拉斯算子的吸收狀態方程建立成像的正向模型，分別考慮組織中滿足頻率冪律的聲吸收和色散，反轉吸收比例系數的符號，同時色散項保持不變，可以補償聲吸收。通過在時間上反復迭代求解測量數據與聲吸收和色散比例系數之間的耦合聲學方程，重建光吸收分布圖。為了對重構過程進行正則化，可利用 Tukey 窗口對吸收和色散項進行頻域濾波，提高迭代的收斂速度。
此類方法適用于具有一般冪律吸收行為的有損聲學介質（具有一維、二維、三維以及任意冪律吸收參數），補償滿足冪律模型的聲衰減，提高重建圖像的分辨率，對更深層的組織特征成像。其局限性在于必須使用 TR 方法重建圖像，與其他快速重建算法相比效率較低，限制了其在實時成像中的應用。而且組織的聲吸收系數和冪律指數多數是在有限的頻率范圍內（通常低于 10 MHz），不能保證這些值在適用于 PAT 的更寬頻率范圍內也有效。
3 結語
圖像重建方法的選擇直接影響 PAT 成像質量。不同成分的組織通常具有不同的聲學特性，且其空間變化是未知的，對組織中聲速恒定和聲學特性分布均勻的理想化假設會導致重建圖像中存在細節模糊、目標錯位及嚴重偽影等問題。針對具有不均勻聲學特性組織的 PAT 成像，目前已提出的圖像重建算法能夠有效解決由組織的聲速不均勻、散射和吸收等造成的圖像質量下降問題，減少偽影，改善圖像的聚焦效果。在后續的研究中，可將修正的時域重建法與標準圖像重建算法（如 TR 等）相結合，補償組織的聲速不均勻性，以獲得更好的重建效果。在解決聲散射問題時，Green 函數重構法具有較明顯的優勢，為了更好地解決實時成像中的聲散射問題，可將 Green 函數重構法與基于散射體分布的統計加權方法有效地結合，更好地發揮二者的優勢。對于聲衰減介質的成像，可考慮將奇異值分解法與基于 TR 的方法相結合，從而以更快的速度和更高的精度來重建圖像，補償聲衰減。同時，在補償聲衰減時，均應考慮組織中聲速分布的不均勻性，在估算聲速分布的基礎上，再進行聲衰減的補償。
近年來，深度學習算法，特別是卷積神經網絡（convolutional neural network，CNN）已迅速成為醫學圖像分析的一種新方法[44-47]，在 PAT 圖像重建方面也顯示出了巨大潛力，可以快速重建出準確、高質量的圖像[48-49]。對于 PAT 成像中的聲衰減補償問題，可以將 CNN 與標準重建方法相結合，即首先利用大量數據集對 CNN 進行訓練，得到合適的網絡權值；然后利用標準重建方法，根據采集到的光聲測量數據重建出包含失真與偽影的低質量圖像；最后用訓練好的 CNN 進行圖像后處理，去除重建圖像中的偽影，提高圖像質量[50-51]。與其他方法相比，基于 CNN 的方法不僅計算速度快，且重建圖像的質量更高，是未來研究的主要方向。
除了聲速不均勻、聲散射和聲吸收以外，還有很多造成 PAT 圖像重建誤差的因素，如聲波在傳播過程中的折射、激光在組織表面的不均勻照射、有限角度掃描和稀疏測量等，都是提高 PAT 成像質量需要解決的關鍵問題。
引言
生物光聲層析成像（photoacoustic tomography，PAT）是一種多物理場耦合的新型復合功能成像方法，對疾病的早期檢測和準確診斷具有極高的研究價值[1]。其成像參數是組織的光吸收系數和散射系數，物理基礎是生物組織的光聲效應[2]。如圖 1 所示，用短脈沖激光照射組織，組織吸收激光能量產生熱彈性膨脹，并以超聲波的形式向表面傳播，即光聲信號[1]。超聲換能器檢測到光聲信號后，采用適當的圖像重建算法即可反演計算出組織內部的光吸收分布或初始聲壓分布[3]，而組織的光吸收特性與其血紅蛋白濃度和分子組成密切相關[4]。PAT 既具有光學成像的高對比度，又具有超聲成像的高穿透深度[5]，在血管成像[6]、血液流速檢測[7]、腫瘤檢測[8]及腦結構和功能成像[9]等領域具有重要意義。
圖1
PAT 成像原理示意圖
Figure1.
Schematics of PAT imaging
有效的光聲信號處理方法和高效的圖像重建算法是保證成像精度和速度的關鍵。目前較為成熟的 PAT 圖像重建算法有濾波反投影（filtered back-projection，FBP）、時間反演（time reversal，TR）、迭代重建法、相控聚焦法和基于傅里葉變換的重建算法等[10-11]。相對于純光學成像，PAT 的高成像深度和高分辨率源于超聲波在生物組織內傳播時的低散射。因此為了簡化問題，多數光聲圖像重建算法通常都假設組織的密度和超聲波在組織中的傳播速度是恒定的，不考慮聲散射，且成像組織是聲學特性均勻的非衰減聲學介質。
然而在實際應用中，考慮到人體組織的復雜性，不同成分的組織通常具有不同的聲學特性，而且其空間變化也是未知的。例如：軟組織內部的聲速差異可以達到 10%（皮下脂肪中的聲速是 1 400～1 450 m/s，而正常薄壁組織和基質組織中的聲速是 1 500～1 560 m/s[12]）。聲速差異會影響光聲信號到達超聲換能器的時間，進而影響根據傳播時間計算出的聲源到超聲換能器的距離。而恒定聲速的假設會導致重建圖像中目標組織的位置與其實際位置之間存在誤差，進而使圖像中存在細節模糊、目標錯位及嚴重偽影等問題[13]。此外，成像組織聲學特性的不均勻性會造成超聲波在組織中傳播時存在反射、散射、擴散和吸收等現象[14]，使部分光聲信號偏離原來的傳播方向，進而導致光聲信號的衰減，并且頻率越高衰減越嚴重，特別是在聲學特性不匹配的組織邊界處，如骨骼、肺、空腔等。此時，若假設成像組織的聲學特性均勻，則會導致重建圖像中包含嚴重的失真和偽影，降低圖像質量，同時可能降低成像深度。因此組織聲學特性不均勻問題是 PAT 圖像重建中需要解決的關鍵問題，對提高成像質量、減少圖像偽影和失真等具有重要意義。
本文對聲學特性不均勻（主要是聲速不均勻和聲衰減）組織的 PAT 圖像重建的主要方法進行綜述，討論各方法的優勢和不足，并展望未來可能的發展方向。
1 聲速不均勻的解決方法
在進行 PAT 成像前，采用超聲透射斷層掃描（ultrasonic transmission tomography，UTT）可獲取成像組織的聲速分布，但需要增加額外的成像設備和成像過程，限制了其應用范圍。而將有關聲速分布的先驗知識結合到圖像重建中，補償由于聲速變化引起的測量光聲場的像差，可以有效改善圖像的聚焦效果，提高圖像質量[13, 15-17]。常用的方法有：有限元法（finite element method，FEM）[18-19]、基于最優聚焦的平均聲速估計法[20]、全波迭代法[4, 21-23]、修正的時域重建法和修正的反卷積重建法[24]等，各方法的性能如表 1 所列。
表1
解決 PAT 成像中聲速不均勻方法的性能比較
Table1.
Performance of methods for PAT in tissues with variable speed of sound
1.1 有限元法
運用有限元方法將光聲波動方程所在的連續場分割成有限個子空間區域，然后應用插值函數表示每個子空間區域的解，將一個連續的無限自由度問題轉化成離散的有限自由度問題進行求解，可以只利用光聲信號同時重建出組織的光吸收分布和聲速分布的低頻成分。
該方法的重建精度高，求解域可以是任意形狀，可得到較好的目標位置信息[25]。缺點是計算量非常大，計算速度慢，且只能重建出待測組織圖像中的低頻成分[18]。
1.2 最優聚焦法
可將聚焦光聲信號絕對差的平均和作為聚焦品質因子 Q，評估光聲信號是否具有最優聚焦，其定義為：
|
其中，M 和 N 分別是焦點和超聲探測器的數量，p 和 分別是聚焦光聲信號的聲壓及其平均值。最優聚焦法的基本思想是：改變聲速，同時針對感興趣區域（region of interest，ROI）中的每條掃描線迭代計算 Q，使 Q 取得最小值的聲速即為 ROI 的最佳平均聲速，再采用標準重建法重建初始聲壓分布圖。
該方法可以提高重建圖像中目標的對比度和空間分辨率，計算簡便。局限性是需要有關目標粗略位置的先驗信息，而且如果 ROI 過小，可能增加迭代計算時間。
1.3 全波迭代重建方法
全波方法用于聲學非均勻有損介質的迭代 PAT 圖像重建，基本思想是建立基于精確光聲波動方程的離散成像模型[4]，利用 k-空間偽譜法計算出前向投影算子和后向投影算子，進而得到離散成像算子（即系統矩陣）及其伴隨或反投影算子，再結合正則化方法［如 TV（total variation）正則化］和迭代算法［如快速迭代收縮閾值算法（fast iterative shrinkage thresholding algorithm，FISTA）］求解偏最小二乘問題，重建初始聲壓分布圖。
該方法能夠有效地補償聲像差，減少重建圖像中的偽影和失真，與 TR 法相比有較高的分辨率。利用 k-空間偽譜法求解耦合方程，采用快速傅立葉變換法計算空間偏導數，比有限差分法具有更高的計算效率，且空間和時間采樣要求較少。局限性是圖像重建精度易受光聲信號中噪聲的影響，且計算時一般需要較多的迭代次數。
1.4 修正的時域重建法和修正的反卷積重建法
修正的時域重建法和修正的反卷積重建法都是在球面掃描條件下，利用在不同位置采集的光聲信號之間的相關性估計超聲波在空間兩點之間的傳播時間，并據此分別修正傳統的時域重建法和反卷積重建法[24]，補償聲速的不均勻性。兩種方法均無需預知聲速分布即可直接重建出目標圖像，且運算速度快。
修正的時域重建法的基本思想是：首先，依據在關于原點對稱的位置處探測到的光聲信號之間相關性最強的原理，估算超聲波在空間兩點之間的傳播時間；其次，利用比例函數關系，求得組織中任意位置處產生的光聲信號到達超聲換能器的時間；最后，根據傳播時間計算光吸收分布。與標準重建方法相比，采用該方法重建的圖像精度更高，邊緣更清晰，且具有良好的抗噪性能，對隨機噪聲不敏感。不足之處是存在理想化的假設，如忽略介質分界面對聲波的折射因素，認為光聲信號從組織內到超聲換能器是沿直線傳播的，導致估計的超聲波傳播時間存在一定誤差，而且聲速不均勻性越強或者待測組織越大，誤差也越大。同時，在計算超聲波傳播時間時，近似認為積分面是一個平面，可能降低重建圖像的質量。
修正的反卷積重建法的原理是：首先，構造一個函數 C（r）：
|
其中，， 是常數且，r0 是超聲換能器的位置。采用修正的時域重建法估計組織中任意位置 r 處產生的光聲信號到達超聲換能器的時間 T（r, r0），并將其代入式（2）進行修正：
|
其次，將組織的光吸收分布 A（r）和 C（r）分別看作一個線性時不變系統的輸入和輸出，即：
|
其中 h（r）是系統的單位沖激響應。當輸入為 δ 函數時，即可得到 h（r）的表達式。則有：
|
其中、和分別是 A（r）、h（r）和 C（r）的傅里葉變換，對進行逆傅里葉變換即可得到光吸收分布。當組織的聲速差異低于 10% 時，采用該方法重建的圖像精度較高且抗噪性能較好。通過使用快速傅里葉變換算法，可獲得極高的計算速度，并且適用于有限角度掃描的情況[24]。但它是一種近似算法，而且待測組織的尺寸越大，近似誤差也越大。
2 聲衰減的補償方法
超聲波在組織中傳播時會發生衰減，并且頻率越高衰減越嚴重。造成衰減的主要因素有超聲波的反射、散射、擴散和吸收等。其中，隨機散射對超聲波的傳播有很大影響，目前解決 PAT 中聲散射問題的主要方法有 Green 函數重構法[14]和基于聲散射體分布的統計加權方法[26-30]等。此外，補償聲衰減的主要方法有矯正不同頻率成分衰減的方法[31]、時變濾波法（time-variant filtering，TVF）[32]、迭代正則化法[33]、基于奇異值分解（singular-value decomposition，SVD）的方法[34-36]和基于 TR 的方法[37-41]等。上述方法優缺點的對比如表 2 所列。
表2
PAT 成像中聲衰減補償方法的優缺點對比
Table2.
Comparison of methods for compensating acoustic attenuation in PAT
2.1 Green 函數重構法
根據互易定理推導出的、適用于常見散射介質的 Green 函數有其局限性。對于采用超聲換能器陣列的散射介質 PAT 成像，可將成像區域劃分為若干網格，根據互易定理在待測點構建一個虛擬聲源，計算該點到每一個換能器的 Green 函數，再結合 TR 方法，得到組織中的初始聲壓分布，進而得到光吸收分布。
采用該方法重構組織內某點的 Green 函數時，無需在該點放置超聲換能器或實際聲源，而是虛構一個虛擬聲源，快速計算該點與散射域另一側測量點之間的 Green 函數，具有計算速度快的優勢，適用于實時 PAT 成像。但是當探測信號頻率降低時，重構的 Green 函數精度也明顯下降，且當有隨機噪聲擾動時，使用理想波形濾波器無法重構出 Green 函數。
2.2 基于聲散射體分布的統計加權方法
在成像目標和超聲探測器之間放置強吸收體，測量到的超聲信號首先用于估計成像目標中散射體的空間分布。然后，在此基礎上計算光聲信號不受反射或散射聲波影響的概率，并利用該概率函數對光聲信號測量值進行加權，從而利用強吸收體產生的、以更高概率從聲源直接傳播到探測器的光聲信號重建光吸收分布，校正由于聲波反射或散射引起的圖像偽影[26-27]。由于確定這種概率的精度取決于被測組織的結構信息，因此可以將有關高散射區近似位置的先驗知識引入重建過程，從而進一步提高層析圖像重建的質量[28-30]。
此類方法能夠去除重建圖像中由于聲散射引起的多數偽影，與假設聲散射體隨機分布的方法相比，能夠顯著提高強散射或反射條件下的光聲成像質量。但是重建圖像的精度受估計的聲散射體分布信息的直接影響，而且需要利用除成像目標之外的強吸收體，可能增加成像復雜度。
2.3 矯正不同頻率成分的聲衰減
根據超聲波衰減理論，由于介質的吸收和散射，光聲信號中的不同頻率成分隨傳播距離的衰減是不同的，高頻成分的衰減較多，低頻成分的衰減較少，并且不同組織對各種頻率超聲波的衰減也是不同的，因此，可對光聲信號不同頻率成分的衰減進行矯正。思路是：組織表面的光聲壓與光吸收分布滿足光聲波動方程，而實際測量的光聲信號是光聲壓與探測器脈沖響應的卷積，因此可以通過頻域的快速傅里葉變換得到光吸收分布與光聲信號測量值之間的關系。再結合超聲衰減公式和組織的衰減參數得到對光聲信號不同頻率成分進行矯正后的光聲圖像重建公式。此方法的計算速度快，可補償光聲信號中各種頻率成分隨距離的衰減，增強信號的高頻成分，突出吸收體的邊界變化和細微結構特征，提高成像分辨率。
2.4 時變濾波法
在假定組織的聲衰減參數分布均勻的前提下，對采集的每個聲壓時間序列進行頻域的時變濾波，修正聲衰減和色散。構造時變濾波器時，可利用 Tukey 窗口進行正則化，并根據聲壓測量值的時頻分布選擇濾波器的截止頻率，與設定固定截止頻率相比，可有效減少引入的高頻噪聲。
該方法適用于任何掃描幾何，且計算效率非常高。局限性在于聲衰減參數在介質內分布均勻的假設可能會引入一定的誤差，且用于正則化的時變窗口的性能也可再進一步優化，以改善濾波結果。
2.5 迭代正則化方法
迭代正則化方法是在基于最小二乘原理的廣義逆算法的基礎上改進而來的一種穩定的廣義逆算法[42]。用于補償聲衰減的迭代正則化 PAT 圖像重建方法是將聲衰減介質的 PAT 逆問題描述為 Hilbert 空間中的一個不適定算子方程[33]，并用迭代正則化法進行求解。其基本思路是：首先，根據衰減定理建立超聲波在衰減介質中傳播的一般模型，并在 Hilbert 空間中設計 PAT 成像的前向算子；然后，選擇適當的正則化方法計算前向算子的伴隨算子，同時對最小二乘問題施加正則化約束；最后，利用迭代法求解最小二乘問題，反演重建出初始聲壓分布圖。
選擇適當的正則化方法是此類方法的關鍵，例如 Landweber 正則化法主要用于計算伴隨算子，由于伴隨項的平滑效應，因而其性能非常穩定[33]。迭代 Tikhonov 正則化在一定條件下能夠達到最優收斂速度，但不適合大規模計算[43]。此外，它需要根據光滑性條件選取最小階數，而不必提高迭代次數。
迭代正則化法是一種通用的反演方法，可以應用于不同的衰減模型和測量幾何，以及稀疏測量數據的 PAT 圖像重建，對于弱衰減的情形，其重建結果甚至優于無衰減的情況。但是在某些應用中它的計算速度很慢，需要較多的迭代次數。
2.6 基于奇異值分解的方法
組織對超聲波的吸收與其頻率直接相關，在時域-頻域中，包含衰減項的光聲波動方程是具有復數波數的非齊次 Helmholtz 方程，直接求解的計算量非常大。而通過對組織中的聲衰減進行建模，可在圖像重建中對其進行校正。
例如，文獻[34]通過積分算子 L 將理想的、未衰減的聲壓信號 p 與在組織表面測量到的衰減的聲壓信號 patt 相關聯，即
|
算子 L 可根據遵循冪律吸收（power law absorption）的有損波動方程在時域-頻域中關聯衰減與未衰減超聲信號得到。對矩陣方程的特征值進行閾值化處理，采用 SVD 計算L的偽逆算子，進而從光聲信號測量值中估計理想的未衰減信號。最后，使用標準重建方法從未衰減信號中重建出組織的光吸收分布圖，實現對聲衰減的補償。但是由于衰減與未衰減信號之間的積分方程一般沒有很好的約束條件，因此這種反演通常是病態的。
為了解決上述問題，可根據光聲信號測量值的傅里葉變換與無衰減理想信號的傅里葉變換得到關聯方程，在對方程離散化得到矩陣方程之后，采用截斷奇異值分解（truncated singular value decomposition，TSVD）計算逆矩陣，進而從衰減的測量信號中獲得對理想信號的近似，實現一步補償[34]，消除聲衰減對光聲信號測量值的影響，再應用標準重建方法（如 TR 等）重建光聲圖像。
此外，文獻[36]提出一種將聲壓信號測量值 patt 直接與重建的初始聲壓信號 p0（包括衰減補償）相關聯的方法，關聯矩陣是時間、衰減系數和二維空間頻率的函數。采用 TSVD 計算矩陣方程的正則化偽逆，進而重建初始聲壓分布圖和光吸收分布圖。該方法適用于平面探測器幾何結構，能夠表征重構衰減介質中不同深度處光吸收分布的不適定性。
2.7 基于 TR 的方法
大孔徑超聲探測器的成像分辨率受到超聲波帶寬和衰減的限制，而且頻率越高衰減越嚴重，如果忽略這種影響，會導致重建圖像中小目標的分辨率嚴重下降。針對此問題，文獻[37-39]中設計了一種簡單的裝置，利用一個旋轉軸將平面探測器分隔成垂直于旋轉軸的線探測器。然后，利用在固定旋轉角度掃描的一組光聲信號，使用二維 TR 方法重建初始聲壓的線性投影。最后，計算與旋轉軸正交的平面中的二維逆 Radon 變換，精確重建成像目標的三維圖像。該方法可以精確地解決任意封閉檢測面的光聲逆問題，且能部分補償聲衰減對圖像重建的影響。其局限性在于有關聲衰減的兩個波動方程的時間反演是不穩定的，不能完全補償聲衰減的影響，重建圖像中存在一定的模糊。
此外，還可將 TR 重建方法與有損波動方程相結合，即反轉斯托克斯有損波動方程中吸收項的符號，那么聲波在組織中的傳播就是增強而不是衰減的。對高于預設截止頻率的信號進行濾波，使重建過程正則化，避免含噪聲的高頻信號增長過快。最后對聲學正問題的數值模型進行時間反演，在每個時間步長上，將測量的時變聲壓信號以時間反轉的順序作為 Dirichlet 邊界條件。當時間 t 達到 0 時，即可重建出初始聲壓分布。
將 TR 圖像重建和吸收狀態方程相結合，使用基于分數拉普拉斯算子的吸收狀態方程建立成像的正向模型，分別考慮組織中滿足頻率冪律的聲吸收和色散，反轉吸收比例系數的符號，同時色散項保持不變，可以補償聲吸收。通過在時間上反復迭代求解測量數據與聲吸收和色散比例系數之間的耦合聲學方程，重建光吸收分布圖。為了對重構過程進行正則化，可利用 Tukey 窗口對吸收和色散項進行頻域濾波，提高迭代的收斂速度。
此類方法適用于具有一般冪律吸收行為的有損聲學介質（具有一維、二維、三維以及任意冪律吸收參數），補償滿足冪律模型的聲衰減，提高重建圖像的分辨率，對更深層的組織特征成像。其局限性在于必須使用 TR 方法重建圖像，與其他快速重建算法相比效率較低，限制了其在實時成像中的應用。而且組織的聲吸收系數和冪律指數多數是在有限的頻率范圍內（通常低于 10 MHz），不能保證這些值在適用于 PAT 的更寬頻率范圍內也有效。
3 結語
圖像重建方法的選擇直接影響 PAT 成像質量。不同成分的組織通常具有不同的聲學特性，且其空間變化是未知的，對組織中聲速恒定和聲學特性分布均勻的理想化假設會導致重建圖像中存在細節模糊、目標錯位及嚴重偽影等問題。針對具有不均勻聲學特性組織的 PAT 成像，目前已提出的圖像重建算法能夠有效解決由組織的聲速不均勻、散射和吸收等造成的圖像質量下降問題，減少偽影，改善圖像的聚焦效果。在后續的研究中，可將修正的時域重建法與標準圖像重建算法（如 TR 等）相結合，補償組織的聲速不均勻性，以獲得更好的重建效果。在解決聲散射問題時，Green 函數重構法具有較明顯的優勢，為了更好地解決實時成像中的聲散射問題，可將 Green 函數重構法與基于散射體分布的統計加權方法有效地結合，更好地發揮二者的優勢。對于聲衰減介質的成像，可考慮將奇異值分解法與基于 TR 的方法相結合，從而以更快的速度和更高的精度來重建圖像，補償聲衰減。同時，在補償聲衰減時，均應考慮組織中聲速分布的不均勻性，在估算聲速分布的基礎上，再進行聲衰減的補償。
近年來，深度學習算法，特別是卷積神經網絡（convolutional neural network，CNN）已迅速成為醫學圖像分析的一種新方法[44-47]，在 PAT 圖像重建方面也顯示出了巨大潛力，可以快速重建出準確、高質量的圖像[48-49]。對于 PAT 成像中的聲衰減補償問題，可以將 CNN 與標準重建方法相結合，即首先利用大量數據集對 CNN 進行訓練，得到合適的網絡權值；然后利用標準重建方法，根據采集到的光聲測量數據重建出包含失真與偽影的低質量圖像；最后用訓練好的 CNN 進行圖像后處理，去除重建圖像中的偽影，提高圖像質量[50-51]。與其他方法相比，基于 CNN 的方法不僅計算速度快，且重建圖像的質量更高，是未來研究的主要方向。
除了聲速不均勻、聲散射和聲吸收以外，還有很多造成 PAT 圖像重建誤差的因素，如聲波在傳播過程中的折射、激光在組織表面的不均勻照射、有限角度掃描和稀疏測量等，都是提高 PAT 成像質量需要解決的關鍵問題。
