光聲成像（PAI）技術是一種發展迅速的新型復合生物醫學成像技術，可實現組織的結構成像和功能成像。在對活體的成像中，目標的呼吸、心跳或眼球轉動等可導致重建圖像中出現運動偽影，降低成像分辨率，增加獲得高質量光聲圖像的難度。本文分別對光聲顯微成像（PAM）和光聲層析成像（PAT）中校正運動偽影的方法進行歸納和總結，討論不同方法的優勢和不足，展望未來可能的發展方向。
引用本文： 杜潔潔, 孫正. 光聲顯微與光聲層析成像中運動偽影校正方法的研究進展. 生物醫學工程學雜志, 2021, 38(2): 369-378. doi: 10.7507/1001-5515.202009062 復制
引言
生物光聲成像（photoacoustic imaging，PAI）是一種多物理場耦合的新型成像技術，它利用納秒級短脈沖激光照射生物組織，組織內的發色團快速吸收光能量，使周圍組織的溫度上升繼而受熱膨脹，激發出寬帶超聲波（即光聲信號）并向組織表面傳播，位于目標周圍的超聲換能器采集光聲信號并送入計算機內，通過求解聲學逆問題重建目標的初始聲壓分布圖或光吸收能量分布圖[1]。在此基礎上，通過求解光學逆問題還可重建組織的光學特性參數和熱膨脹系數分布圖，實現定量光聲成像。PAI 技術結合了純光學成像的高對比度和超聲成像的高穿透深度的優點，可提供從細胞器到器官的多尺度的解剖、功能和分子對比信息，在生物醫學領域具有廣闊的應用前景[2]。按照成像模式的不同，PAI 可分為光聲顯微成像（photoacoustic microscopy，PAM）和光聲層析成像（photoacoustic tomography，PAT）。
在活體成像中，生物體的運動通常會導致圖像中產生偽影，降低圖像質量。生物體的運動包括自主運動和不自主運動，其中呼吸、心跳和眼球轉動等屬于不自主運動，身體移動屬于自主運動。針對自主運動，可通過使目標在成像階段保持靜默以避免產生運動偽影[3]。例如，在對小動物進行全身光聲成像時，通常采用不同的安裝方案，通過機械夾緊和固定減少運動偽影[4-6]。但是并非所有身體部位都適合夾持，例如眼睛，夾持會阻礙血流，可能降低圖像質量[3]。因而此類方法只適用于特定的應用場合和成像系統，不具有普適性，不能從根本上解決運動偽影問題。本文對 PAM 和 PAT 運動偽影校正方法的研究現狀進行綜述，分析各方法的優勢和不足，同時展望未來可能的發展方向。
1 PAM 運動校正方法
1.1 PAM 原理及目標運動對成像質量的影響
PAM 是一種聚焦型掃描成像技術，它采用聚焦型入射光或聚焦型超聲探測的方法進行成像，成像深度為毫米級[7-8]，現已應用于血管生物學、神經學、腫瘤學、皮膚學以及眼科學等領域的研究中[9-12]。如圖 1 所示，根據成像系統中光束焦斑與聲束焦斑的相對大小，PAM 分為光學分辨率光聲顯微成像（optical resolution PAM，OR-PAM）和聲學分辨率光聲顯微成像（acoustic resolution PAM，AR-PAM）[3, 7]，二者的區別在于光束和聲束聚焦的相對大小，圖像重建過程基本一致[10]。AR-PAM 的成像深度較大（1～10 mm）[13]，橫向分辨率取決于探頭焦點寬度，縱向分辨率取決于探頭帶寬。OR-PAM 的成像深度較小（1～2 mm），與光學成像接近，橫向分辨率為微米級甚至亞微米級，縱向分辨率與超聲換能器的頻率呈正比，特別適合觀察細胞和微血管水平的生理過程[10]。
圖1
PAM 聚焦方式對比
Figure1.
Comparison of focusing modes of PAM
PAM 的成像分辨率高，因而對運動偽影非常敏感。當對大面積區域成像時，由于樣品表面通常不平坦，因此需要針對不同分區調整成像系統的焦平面，獲得樣品的多幅圖像，再通過圖像拼接得到整個區域的圖像。獲取各分區的圖像時，對應區域的目標運動會導致不同分區圖像之間的不匹配，使拼接后的圖像存在扭曲和變形。此外，運動偽影還會影響圖像后處理（如信號增強、提取和特征識別等）的效果，導致關鍵信息丟失，進而降低定量成像的精度[3]。目前對 PAM 運動偽影的主要校正方法包括動態參考對象法、運動量化法、亞像素運動估計法、Demons 圖像配準法和基于深度學習的方法，如表 1 所列。
表1
PAM 運動校正的主要方法對比
Table1.
Comparison of main methods of motion correction in PAM
1.2 光聲介觀成像的運動校正方法
光柵掃描光聲介觀成像（raster-scan optoacoustic mesoscopy，RSOM）是一種聲學分辨率光聲顯微成像，其信號采集速度取決于激光脈沖的重復頻率（pulse repetition rate，PRR），當 PRR 為 1 kHz 時，掃描一塊 250 × 250 點的區域約需 60 s[14]。受信號采集速度的限制，成像目標的運動（如頻率為 30 次/分的呼吸運動）會改變聚焦超聲換能器與感興趣區域的相對位置，導致圖像質量下降。
1.2.1 動態參考對象法
Schwarz 等[14]提出一種基于動態參考對象的 RSOM 運動校正方法，假設皮膚表層是連續的，根據惠更斯-菲涅耳原理，連續皮膚表層發出的超聲波會產生一個連續次波前，然而呼吸運動造成皮膚表層相對于換能器的垂直運動，使次波前發生中斷（如圖 2 所示）。觀測三維波形圖中皮膚層產生的超聲波波前的垂直中斷，利用滑動平均濾波算法對受損表層進行平滑和濾波處理，生成一個光滑、連續的仿真皮膚表層。最后計算仿真表層和受損表層之間的偏移，并據此調整換能器的焦點。該方法可有效減少重建圖像中由皮膚表層相對于超聲換能器的垂直運動所致的偽影，提高信噪比和分辨率。但由于穿透深度和成像區域有限，較難找到相對于成像目標靜止的動態參考目標，因而其適用范圍有限。
圖2
采用動態參考對象法校正運動前后的皮膚表層[14]
a. 未校正的健康皮膚；b. 校正后的健康皮膚；c. 未校正的銀屑病斑塊；d. 校正后的銀屑病斑塊
Figure2. Skin surface detected in raw RSOM data before and after motion correction with the dynamic reference method[14]a. healthy skin uncorrected; b. healthy skin corrected; c. psoriasis plaque uncorrected; d. psoriasis plaque corrected
1.2.2 運動量化法
動態參考對象法假設目標運動是剛性的，當目標運動量較大時，非剛性運動會影響其校正效果。針對該問題，Aguirre 等[15]提出利用垂直位移對圖像質量的改善程度進行量化，校正由皮膚垂直運動所致的偽影。首先，將 RSOM 掃描頭固定在所需位置，采集 10 s 的光聲信號，激光的重復頻率與常規 RSOM 相同。采集與每個激光脈沖對應的一維光聲信號，將所有一維信號組合起來形成一個二維平面（即運動圖），用于觀察運動。計算第一個一維信號和其余一維信號之間互相關函數的最大值，得到位移函數：
 |
其中，A（τ,vt）是在時刻 τ 采集的一維信號，v 是超聲波在水中的傳播速度，t 是光聲信號從聲源到探測器的傳播時間，vt 是反距離（reciprocal distance），“*”是互相關運算。然后，利用正常光柵掃描方式采集 RSOM 正弦圖，并計算累積互相關表面（cross-correlation surface，CS）函數：
 |
其中，xi 和 yj 分別是快掃描軸和慢掃描軸上采集網格的位置，（xi,yj）是正弦圖中一條 A 掃描線的位置，tk 是每條 A 掃描線的時間，b（tk, xli, yN）表示時刻 tk、坐標（xi ? l/2, yj）到（xi+l/2, yj）范圍內的 A 掃描線對應的 B 掃描面。最后，計算一個合成平滑 CS 函數，將正弦圖每一行的垂直位移量化為合成 CS 和原始 CS 之間的反距離，繼而根據探測器到皮膚的實際距離重建圖像。如圖 3 所示，校正后圖像的清晰度明顯提高。該算法是全自動的，無需進行圖像分割，對中低運動區域的校正效果優于高運動區域。當對黑色素瘤成像時，病灶表面對激發光的強烈吸收可能導致陰影效應，此時可采用較長波長的激發光減少陰影效應。
1.2.3 亞像素運動估計法
運動量化法可實現像素級的運動校正，但其忽略了亞像素級的錯位對圖像質量的影響。Cheng 等[16]提出了一種亞像素級運動估計方法，首先進行亞像素級全局運動校正，即將前一幀 B 掃描圖像作為參考圖像，估計后一幀的運動向量，然后利用樣條插值法校正運動，最后進行亞像素級垂直運動校正，即通過累積互相關差分（cross-correlation differential，CD）函數計算 B 掃描面上 A 掃描線的總垂直位移：
 |
其中，xi 和 yj 是采集網格在 x 軸和 y 軸方向上的位置，yN 是當前 B 掃描圖像在 y 方向上的位置，cd（xi,yN）是在當前 B 掃描圖像中位于 xi 處的 A 掃描線的總位移，▽x 是 x 方向上的差分算子，dzik 表示（xi, yN）處的 A 線與（xi,yN ? 1）處的 A 線之間亞像素級相對位移的估計值。如圖 4 所示，進行全局運動校正后，明顯改善了水平和垂直方向的偽影，但仍存在條紋偽影，而進行垂直運動校正后，皮膚表面的光滑程度明顯提高。該算法具有較高的校正精度，校正后的圖像可以清晰顯示血管邊界，甚至是單根毛細血管。但是其所需的校正時間較長，且只能校正線性運動偽影，對旋轉所致偽影的校正能力有限。
1.3 基于 Demons 圖像配準的運動補償
Demons 算法是一種基于光流理論的圖像配準方法，通常用于非剛性圖像配準，其基本思想是假設運動目標的灰度不隨時間變化，那么圖像配準可以看作是浮動圖像中的各像素向參考圖像逐步擴散的過程，擴散速度取決于參考圖像的灰度梯度信息[17]。針對 OR-PAM 圖像中運動偽影的校正問題，Zhao 等[3]提出一種改進的 Demons 圖像配準算法，其原理是：首先通過比較相鄰層之間的互信息（mutual information，MI）值[18]，從原始三維光聲信號中識別和提取被運動損壞的 B 掃描圖像 K，對 K 的非運動損壞相鄰層，即 K-1 和 K+1 層進行樣條插值得到一個新層 K1；然后，將 K 層作為浮動圖像，K1 層作為參考圖像，利用 Demons 算法對 K 層和 K1 層進行非線性配準，使 K 層向 K1 層變形，得到一個新層 K2，根據光流方程得到位移 μ：
 |
其中，f 和 r 分別是浮動圖像 F 和參考圖像 R 在像素點 p 處的灰度值，▽r 是參考圖像在 p 點處的灰度梯度。最后，計算 K1、K2 層與 K? 1、K+1 層之間的 MI 值，用 MI 最大值所在層（K1 或 K2）替換原始數據中的 K 層。如圖 5 所示，運動校正前血管圖像中存在鋸齒狀的運動偽影，校正后圖像的清晰度明顯提高。該方法不需要動態參考對象，并且只需一組三維光聲信號數據即可，但是數據采集和處理速度較慢。利用浸水式微電機掃描系統可提高數據采集速度[19]，此外利用高性能圖形處理單元（graphic processing unit，GPU）進行并行計算，也可以縮短數據處理時間，實現實時成像。
1.4 基于深度學習的運動校正
近年來，深度學習（deep learning，DL）特別是卷積神經網絡（convolutional neural network，CNN）已經迅速成為醫學圖像分析與處理的首選方法，在快速重建高質量的圖像方面也表現出了巨大潛力[20]。例如 Chen 等[21]采用 DL 方法實現了 OR-PAM 圖像中運動偽影的校正，其原理是：首先利用希爾伯特變換提取每個深度方向上光聲信號的包絡，并取其軸向投影的最大幅值，獲得最大幅值投影（maximum amplitude projection，MAP）圖像；然后搭建 CNN 模型，根據先驗知識設計損失函數，使用反向傳播算法更新網絡參數（例如卷積核和神經元偏移向量），通過最小化損失函數對模型進行訓練，提升模型的學習能力；最后利用該 CNN 模型進行圖像后處理，校正運動偽影，結果如圖 6 所示。增大卷積核可提高 CNN 的性能，但會降低處理效率，因此卷積核大小的選擇應是性能與效率之間權衡的結果。雖然該方法是為 OR-PAM 設計的，但是可以通過改變訓練集實現其他成像模式的運動偽影校正。然而 CNN 模型的搭建及訓練需要大量成像數據，到目前為止還沒有針對光聲圖像重建的公開數據集，可用的臨床病例數據也十分有限。
2 PAT 運動校正方法
2.1 PAT 原理及目標運動的影響
如圖 7 所示，PAT 是用非聚焦短脈沖激光作為激勵源照射待成像組織，利用具有掃描功能的超聲換能器或換能器陣列采集組織產生的光聲信號，組織內不同深度的光吸收體產生的光聲信號到達超聲換能器的時間不同，利用時間分辨技術可獲得多個不同層析面的光聲信號，進而重建組織的二維或三維初始聲壓分布圖或光吸收能量分布圖[22]。PAT 的成像深度可達數厘米，空間分辨率可達數十微米，在冠狀動脈粥樣硬化性疾病、乳腺癌和消化系統病變等疾病的早期診斷方面具有良好的應用前景[23]。
圖7
PAT 成像原理示意圖
Figure7.
Schematic diagram of PAT principle
由于 PAT 具有在單個納秒級激光脈沖作用下從整個成像區域采集圖像數據的能力，因此它能夠有效消除單幀圖像中的運動偽影。但是，在多光譜光聲層析成像（multi-spectral optoacoustic tomography，MSOT）或基于掃描的成像系統中，多幀圖像分析通常是必不可少的，因此需要采用有效的方法校正由運動引起的圖像失真。例如，可以減少連續激光脈沖之間的延遲抑制運動[24-25]，但需要專用的激光器系統。除此之外，主要的運動校正方法包括呼吸門控技術、運動聚類法、基于混合成像模型的方法和幀平均-幀間運動補償法，如表 2 所列。
表2
PAT 運動校正方法的對比
Table2.
Comparison of existing methods of motion correction for PAT
2.2 呼吸門控技術
門控技術是醫學成像中減少由心跳和呼吸所致運動偽影的常用方法，分為前瞻性和回顧性兩種方式。前瞻性門控是指采用專門的心電或呼吸門控采集裝置，只在每個心動或呼吸周期中的相同相位處采集圖像或信號，由于每個心動或呼吸周期只采集一幀圖像，因而與連續采集的方式相比，這種方式延長了成像時間。回顧性門控是指在整個心動或呼吸周期內連續采集圖像或信號，同時記錄心電信號或呼吸波形。待掃描結束后，將采集的圖像或信號與心電或呼吸信號進行對照和分析，選擇在相同相位處采集的圖像或信號，該技術在許多需要精確定位器官的成像中是必不可少的[26-27]。
小動物全身 PAT 系統的幀頻通常小于 1 Hz，而小鼠的呼吸頻率可高達 3 Hz，所以在對胸部以下區域進行 PAT 掃描時，呼吸是產生運動偽影的主要原因[28]，導致圖像模糊和分辨率下降。回顧性呼吸門控技術是消除 PAT 呼吸運動偽影的常用方法，包括硬件門控和軟件門控兩種方式。回顧性硬件呼吸門控是在成像目標自由呼吸時采集光聲信號，同時利用外部設備監測呼吸波形，通過分析光聲信號數據，選擇出具有相同呼吸相位的信號[27]。例如，Xia 等[29]利用回顧性呼吸門控技術對小鼠的肝臟和腎臟 PAT 圖像進行運動校正，將小鼠完全或部分浸入水中，利用壓力傳感器持續監測水位波動，并據此監測呼吸波形。然后，根據呼吸相位對光聲信號進行排序和聚類，采用在同一呼吸相位采集的光聲數據重建圖像，結果如圖 8 所示。該方法對呼吸波形的監測不受圖像噪聲的影響，而且數據處理的計算量較低，可直接觀察呼吸周期，提高瞬時分辨率。但是若目標呼吸比較急促或者不均勻，往往會出現漏感應或錯誤感應，導致觸發精度降低，因此需要對目標進行呼吸訓練。
回顧性軟件呼吸門控不需要借助外部硬件設備監測呼吸波形，而是根據采集的光聲信號計算呼吸門控信號[30]。例如，Ron 等[28]提出一種自門控技術，首先將超聲換能器陣列中各陣元采集的光聲信號排列成二維矩陣，其中行表示各信號的 m 個時間樣本，列對應換能器陣列的 n 個陣元，在每個換能器位置共采集 k 幀，并將該幀序列排列成一個 m × n 行和 k 列的二維矩陣，表示在該換能器位置的整個幀序列；然后利用主成分分析法計算相鄰幀之間的自相關矩陣，并對其采用二階 k-means 算法進行聚類[31-32]，將 k 幀分為運動幀和靜態幀兩類；最后利用濾波反投影算法重建三維圖像，結果如圖 9 所示。該方法可在重建圖像之前進行運動校正，提高圖像的清晰度和信噪比，避免對圖像進行后處理。但是有的生物體兩次呼吸之間的停頓超過單次呼吸長度，需根據呼吸特征的先驗知識建立運動幀和靜態幀之間的區分標準。此外，呼吸門控信號提取的準確性會影響運動校正的效果。
圖9
采用自呼吸門控技術對小鼠 PAT 圖像的運動偽影抑制結果[28]
a. 選擇出的幀；b. 丟棄的幀（箭頭指示產生變形的結構）；c. 圖 a 中虛線標記位置的聲壓幅值曲線，其中虛線曲線對應采用所有幀進行圖像重建，實線曲線對應采用選擇幀進行圖像重建
Figure9. Motion rejection results for spiral volumetric optoacoustic tomography[28]a. selected frame; b. rejected frame (distorted structures are marked with arrows); c. amplitude profiles marked in Fig. a (dashed lines) for images reconstructed from all the frames (dashed lines) versus selected frames (solid lines)
2.3 基于運動聚類的校正方法
MSOT 根據不同組織的光譜吸收特性，利用不同波長的脈沖激光照射組織，產生具有組織特異性的光聲信號，可更好地識別和標定組織成分[33]。對于心臟 MSOT，信號采集過程中的心跳運動會導致圖像產生模糊，分辨率降低[34-35]。Taruttis 等[34]提出采用運動聚類技術對心臟 MSOT 圖像進行運動校正，其原理是：首先利用反投影算法重建每個波長對應的 100 幀圖像；然后利用 k-means 算法對 100 幀圖像進行聚類，再對聚類后的圖像求平均，考慮到心室運動存在兩個轉折點，因此聚類的集群數量取值為 2；最后對平均后的圖像進行線性光譜解混。如圖 10 所示，執行聚類的圖像中由運動造成的模糊程度明顯降低，且增加集群數量可進一步降低運動模糊程度，但也會降低信噪比。隨著集群數量的增加，集群之間的距離會減小，也提高了對 k-means 算法聚類能力的要求，因此需要結合圖像模糊程度、信噪比和 k-means 算法的性能設定集群數量。
圖10
由光譜解混獲得的脫氧血紅蛋白 MSOT 圖像[34]
a. 直接對各波長的 100 幀圖像進行平均得到的光譜解混圖像；b. 只對各波長的較大集群圖像進行平均得到的光譜解混圖像；c. 只對各波長的較小集群圖像進行平均得到的光譜解混圖像
Figure10. Motion correction for unmixed MSOT images of deoxyhemoglobin contribution[34]a. unmixed image after all 100 images were averaged for each wavelength; b. unmixed image after averaging only the larger clusters for each wavelength; c. unmixed image after averaging only the smaller clusters for each wavelength
2.4 基于混合成像模型的運動校正
Chung 等[36]提出一種基于混合成像模型的圖像重建方法，建立考慮呼吸運動的成像模型，并求解使模型輸出值與光聲測量數據之間誤差最小的非線性最小二乘問題，同時實現運動校正和圖像重建。其原理是：假設聲速恒定，成像目標在連續兩次照明和測量之間的變化非常緩慢，用聲壓的圓 Radon 變換（circular radon transform，CRT）[37]表示聲壓測量數據：
 |
式中，r 是換能器與變形區域之間的距離，? 是換能器陣列旋轉的角度，z（?）=（cos?, sin?）是測量位置，f? 是初始聲壓 f 在 t? 時刻的變形，Φ（?, x）是點 x 到點 y 的變形映射，Ψ（?,·）是 Φ（?,·）的逆映射[38]，Jx［Φ（?, x）］是變量 x→y = Φ（?, x）變形映射的雅可比行列式的絕對值。那么，PAT 圖像重建問題可簡化為逆 CRT 問題。假設由呼吸產生的運動是垂直拉伸，建立光聲成像的離散模型：
 |
其中，g =［g1 ··· gn］T∈Rn × m 是測量值矩陣，n 是換能器的測量位置數，f 是初始聲壓矩陣，γ =［γ1 ··· γn］T 是運動參數矩陣，γi 是 zi 處的運動參數，e =［e1 ··· en］T∈Rn × m 是獨立同分布的高斯噪聲矩陣，A（γ）=［A1K（γ1） ··· AnK（γn）］T∈Rnm × N，Ai 是每個換能器位置對應的投影矩陣，Aif 表示 f 在以 zi 為中心的圓上的 CRT，fi = K（γi）f 是第 i 個變形后的聲壓矩陣，K（γi）∈RN × N 表示一維拉伸/膨脹變形。當運動參數是變量時，f 的估計問題定義為如下非線性最小二乘問題：
 |
式中，λ 是正則化參數。利用變量投影算法[39]分離線性變量和非線性變量，進而采用高斯牛頓迭代法[40]對非線性參數γ進行優化，其中利用復合最小二乘 QR 分解（hybrid least square QR factorization，HyBR）算法自動選擇正則化參數，并求解關于 f 的線性最小二乘問題，結果如圖 11 所示。該方法不僅可應用于垂直拉伸模型，還可以擴展到其他變形模型，例如旋轉或者平移。
圖11
基于混合成像模型的 PAT 圖像重建結果[36]
a. 忽略運動；b. 采用 LSQR 算法求解線性問題；c. 采用 HyBR 算法求解線性問題
Figure11. PAT images reconstructed based on a hybrid model[36]a. initial regularized reconstruction ignoring motion; b. reconstruction obtained by using LSQR to solve the linear problem; c. reconstruction obtained by using HyBR to solve the linear problem
2.5 幀平均-幀間運動補償
對于深層組織的 PAT，由于產生的聲壓信號較弱，導致重建圖像質量較差。為了提高信噪比，可以利用高功率脈沖激光來提高聲壓，但是這種方法存在安全風險（皮膚輻射暴露 < 20 mJ/cm2）。對于采用安全激光脈沖能量的 PAT，可利用幀平均（frame averaging，FA）技術，對相同幾何結構的圖像進行平均以提高圖像的信噪比。但是因為激光脈沖重復頻率的限制以及成像目標運動的影響，采用 FA 技術很難實現高分辨率實時成像。針對該問題，Kim 等[41]將 FA 技術與幀間運動補償（inter-frame motion compensation，IFMC）相結合，即 FA-IFMC 技術，用于補償深層組織 PAT 圖像中的運動偽影。其原理是：首先，利用塊匹配算法[42]和三步搜索法進行運動估計，從最近獲得的圖像中選取局部參考區，在先驗圖像中選取對比區，二者大小相同，用平均絕對誤差（mean absolute difference，MAD）作為匹配準則：
 |
其中，N 和 M 分別是待匹配圖像塊在軸向和橫向上的像素數，Cij 和 Rij 分別是對比區和參考區中像素（i, j）的灰度值。將參考區及其鄰近區域定義為搜索區域，在參考區及其八鄰域方向上以步長逐次減半的方式進行搜索，在搜索區域中找到與對比區最相似的子區域，參考區和相似子區域之間的坐標偏移就是估計出的運動矢量。然后，根據運動矢量，利用 IFMC 比較圖像序列中的相鄰幀實現運動偽影的校正。
如圖 12 所示，與傳統 FA 相比，FA-IFMC 提高了成像的峰值信噪比和軸向分辨率。與幀運動補償（frame motion compensation，FMC）[41]相比，由于 FMC 選取的參考圖像和對比圖像不是相鄰幀，兩幀之間的運動矢量隨著時間的增加而積累，而搜索步長又限制了搜索區域的大小，導致運動矢量很容易超出搜索范圍，而 FA-IFMC 的參考圖像和對比圖像取自相鄰幀，避免了幀平均-幀運動補償（frame averaging-frame motion compensation，FA-FMC）中的問題。但是步長決定了運動估計的精度和搜索范圍，需謹慎選擇，后續研究中可利用變形塊匹配算法進行運動估計，并利用浮動步長提高運動補償的能力。
3 結論
在光聲成像中，目標的呼吸、心跳、脈搏、眼球轉動和身體移動等會導致圖像中產生運動偽影，降低成像分辨率，對運動偽影進行抑制和補償對于提高成像質量具有重要意義。本文分別對 PAM 和 PAT 中運動偽影的主要校正方法進行了歸納和總結。目前，該領域仍面臨諸多問題，例如：為了提高運動聚類法的準確性，需要不斷提高 k-means 算法的聚類性能；多數方法不具有一般性，例如基于 Demons 圖像配準的方法只有在成像區域中血管數量較多的情況下才能取得理想的校正效果，呼吸門控法只適用于校正由呼吸運動產生的偽影；多數方法只能校正垂直方向上由一維運動產生的偽影，針對其他方向上或者由二維運動產生的偽影的校正是未來主要的研究方向；雖然采用運動校正方法可明顯改善圖像質量，但是為了實現實時成像，仍需提高成像系統的數據采集和處理速度。
此外，現有的運動校正方法均是針對體外掃描成像系統的，尚未檢索到針對內窺式光聲成像[43]的運動校正方法。例如，冠狀動脈血管內光聲（intravascular photoacoustic，IVPA）圖像序列[44]中存在由周期性的心臟運動、呼吸運動和動脈內搏動的血流造成的運動偽影，在橫截面圖像中表現為相鄰幀之間血管橫截面的錯位，在縱向視圖中表現為血管壁呈現鋸齒狀邊緣，嚴重影響了 IVPA 圖像序列的視覺效果，降低了血管形態結構和組織成分定量分析以及血管三維重建的精度。因此，針對內窺式光聲成像中運動偽影的校正將是未來該領域的重要研究方向之一。
利益沖突聲明：本文全體作者均聲明不存在利益沖突。
引言
生物光聲成像（photoacoustic imaging，PAI）是一種多物理場耦合的新型成像技術，它利用納秒級短脈沖激光照射生物組織，組織內的發色團快速吸收光能量，使周圍組織的溫度上升繼而受熱膨脹，激發出寬帶超聲波（即光聲信號）并向組織表面傳播，位于目標周圍的超聲換能器采集光聲信號并送入計算機內，通過求解聲學逆問題重建目標的初始聲壓分布圖或光吸收能量分布圖[1]。在此基礎上，通過求解光學逆問題還可重建組織的光學特性參數和熱膨脹系數分布圖，實現定量光聲成像。PAI 技術結合了純光學成像的高對比度和超聲成像的高穿透深度的優點，可提供從細胞器到器官的多尺度的解剖、功能和分子對比信息，在生物醫學領域具有廣闊的應用前景[2]。按照成像模式的不同，PAI 可分為光聲顯微成像（photoacoustic microscopy，PAM）和光聲層析成像（photoacoustic tomography，PAT）。
在活體成像中，生物體的運動通常會導致圖像中產生偽影，降低圖像質量。生物體的運動包括自主運動和不自主運動，其中呼吸、心跳和眼球轉動等屬于不自主運動，身體移動屬于自主運動。針對自主運動，可通過使目標在成像階段保持靜默以避免產生運動偽影[3]。例如，在對小動物進行全身光聲成像時，通常采用不同的安裝方案，通過機械夾緊和固定減少運動偽影[4-6]。但是并非所有身體部位都適合夾持，例如眼睛，夾持會阻礙血流，可能降低圖像質量[3]。因而此類方法只適用于特定的應用場合和成像系統，不具有普適性，不能從根本上解決運動偽影問題。本文對 PAM 和 PAT 運動偽影校正方法的研究現狀進行綜述，分析各方法的優勢和不足，同時展望未來可能的發展方向。
1 PAM 運動校正方法
1.1 PAM 原理及目標運動對成像質量的影響
PAM 是一種聚焦型掃描成像技術，它采用聚焦型入射光或聚焦型超聲探測的方法進行成像，成像深度為毫米級[7-8]，現已應用于血管生物學、神經學、腫瘤學、皮膚學以及眼科學等領域的研究中[9-12]。如圖 1 所示，根據成像系統中光束焦斑與聲束焦斑的相對大小，PAM 分為光學分辨率光聲顯微成像（optical resolution PAM，OR-PAM）和聲學分辨率光聲顯微成像（acoustic resolution PAM，AR-PAM）[3, 7]，二者的區別在于光束和聲束聚焦的相對大小，圖像重建過程基本一致[10]。AR-PAM 的成像深度較大（1～10 mm）[13]，橫向分辨率取決于探頭焦點寬度，縱向分辨率取決于探頭帶寬。OR-PAM 的成像深度較小（1～2 mm），與光學成像接近，橫向分辨率為微米級甚至亞微米級，縱向分辨率與超聲換能器的頻率呈正比，特別適合觀察細胞和微血管水平的生理過程[10]。
圖1
PAM 聚焦方式對比
Figure1.
Comparison of focusing modes of PAM
PAM 的成像分辨率高，因而對運動偽影非常敏感。當對大面積區域成像時，由于樣品表面通常不平坦，因此需要針對不同分區調整成像系統的焦平面，獲得樣品的多幅圖像，再通過圖像拼接得到整個區域的圖像。獲取各分區的圖像時，對應區域的目標運動會導致不同分區圖像之間的不匹配，使拼接后的圖像存在扭曲和變形。此外，運動偽影還會影響圖像后處理（如信號增強、提取和特征識別等）的效果，導致關鍵信息丟失，進而降低定量成像的精度[3]。目前對 PAM 運動偽影的主要校正方法包括動態參考對象法、運動量化法、亞像素運動估計法、Demons 圖像配準法和基于深度學習的方法，如表 1 所列。
表1
PAM 運動校正的主要方法對比
Table1.
Comparison of main methods of motion correction in PAM
1.2 光聲介觀成像的運動校正方法
光柵掃描光聲介觀成像（raster-scan optoacoustic mesoscopy，RSOM）是一種聲學分辨率光聲顯微成像，其信號采集速度取決于激光脈沖的重復頻率（pulse repetition rate，PRR），當 PRR 為 1 kHz 時，掃描一塊 250 × 250 點的區域約需 60 s[14]。受信號采集速度的限制，成像目標的運動（如頻率為 30 次/分的呼吸運動）會改變聚焦超聲換能器與感興趣區域的相對位置，導致圖像質量下降。
1.2.1 動態參考對象法
Schwarz 等[14]提出一種基于動態參考對象的 RSOM 運動校正方法，假設皮膚表層是連續的，根據惠更斯-菲涅耳原理，連續皮膚表層發出的超聲波會產生一個連續次波前，然而呼吸運動造成皮膚表層相對于換能器的垂直運動，使次波前發生中斷（如圖 2 所示）。觀測三維波形圖中皮膚層產生的超聲波波前的垂直中斷，利用滑動平均濾波算法對受損表層進行平滑和濾波處理，生成一個光滑、連續的仿真皮膚表層。最后計算仿真表層和受損表層之間的偏移，并據此調整換能器的焦點。該方法可有效減少重建圖像中由皮膚表層相對于超聲換能器的垂直運動所致的偽影，提高信噪比和分辨率。但由于穿透深度和成像區域有限，較難找到相對于成像目標靜止的動態參考目標，因而其適用范圍有限。
圖2
采用動態參考對象法校正運動前后的皮膚表層[14]
a. 未校正的健康皮膚；b. 校正后的健康皮膚；c. 未校正的銀屑病斑塊；d. 校正后的銀屑病斑塊
Figure2. Skin surface detected in raw RSOM data before and after motion correction with the dynamic reference method[14]a. healthy skin uncorrected; b. healthy skin corrected; c. psoriasis plaque uncorrected; d. psoriasis plaque corrected
1.2.2 運動量化法
動態參考對象法假設目標運動是剛性的，當目標運動量較大時，非剛性運動會影響其校正效果。針對該問題，Aguirre 等[15]提出利用垂直位移對圖像質量的改善程度進行量化，校正由皮膚垂直運動所致的偽影。首先，將 RSOM 掃描頭固定在所需位置，采集 10 s 的光聲信號，激光的重復頻率與常規 RSOM 相同。采集與每個激光脈沖對應的一維光聲信號，將所有一維信號組合起來形成一個二維平面（即運動圖），用于觀察運動。計算第一個一維信號和其余一維信號之間互相關函數的最大值，得到位移函數：
|
其中，A（τ,vt）是在時刻 τ 采集的一維信號，v 是超聲波在水中的傳播速度，t 是光聲信號從聲源到探測器的傳播時間，vt 是反距離（reciprocal distance），“*”是互相關運算。然后，利用正常光柵掃描方式采集 RSOM 正弦圖，并計算累積互相關表面（cross-correlation surface，CS）函數：
|
其中，xi 和 yj 分別是快掃描軸和慢掃描軸上采集網格的位置，（xi,yj）是正弦圖中一條 A 掃描線的位置，tk 是每條 A 掃描線的時間，b（tk, xli, yN）表示時刻 tk、坐標（xi ? l/2, yj）到（xi+l/2, yj）范圍內的 A 掃描線對應的 B 掃描面。最后，計算一個合成平滑 CS 函數，將正弦圖每一行的垂直位移量化為合成 CS 和原始 CS 之間的反距離，繼而根據探測器到皮膚的實際距離重建圖像。如圖 3 所示，校正后圖像的清晰度明顯提高。該算法是全自動的，無需進行圖像分割，對中低運動區域的校正效果優于高運動區域。當對黑色素瘤成像時，病灶表面對激發光的強烈吸收可能導致陰影效應，此時可采用較長波長的激發光減少陰影效應。
1.2.3 亞像素運動估計法
運動量化法可實現像素級的運動校正，但其忽略了亞像素級的錯位對圖像質量的影響。Cheng 等[16]提出了一種亞像素級運動估計方法，首先進行亞像素級全局運動校正，即將前一幀 B 掃描圖像作為參考圖像，估計后一幀的運動向量，然后利用樣條插值法校正運動，最后進行亞像素級垂直運動校正，即通過累積互相關差分（cross-correlation differential，CD）函數計算 B 掃描面上 A 掃描線的總垂直位移：
|
其中，xi 和 yj 是采集網格在 x 軸和 y 軸方向上的位置，yN 是當前 B 掃描圖像在 y 方向上的位置，cd（xi,yN）是在當前 B 掃描圖像中位于 xi 處的 A 掃描線的總位移，▽x 是 x 方向上的差分算子，dzik 表示（xi, yN）處的 A 線與（xi,yN ? 1）處的 A 線之間亞像素級相對位移的估計值。如圖 4 所示，進行全局運動校正后，明顯改善了水平和垂直方向的偽影，但仍存在條紋偽影，而進行垂直運動校正后，皮膚表面的光滑程度明顯提高。該算法具有較高的校正精度，校正后的圖像可以清晰顯示血管邊界，甚至是單根毛細血管。但是其所需的校正時間較長，且只能校正線性運動偽影，對旋轉所致偽影的校正能力有限。
1.3 基于 Demons 圖像配準的運動補償
Demons 算法是一種基于光流理論的圖像配準方法，通常用于非剛性圖像配準，其基本思想是假設運動目標的灰度不隨時間變化，那么圖像配準可以看作是浮動圖像中的各像素向參考圖像逐步擴散的過程，擴散速度取決于參考圖像的灰度梯度信息[17]。針對 OR-PAM 圖像中運動偽影的校正問題，Zhao 等[3]提出一種改進的 Demons 圖像配準算法，其原理是：首先通過比較相鄰層之間的互信息（mutual information，MI）值[18]，從原始三維光聲信號中識別和提取被運動損壞的 B 掃描圖像 K，對 K 的非運動損壞相鄰層，即 K-1 和 K+1 層進行樣條插值得到一個新層 K1；然后，將 K 層作為浮動圖像，K1 層作為參考圖像，利用 Demons 算法對 K 層和 K1 層進行非線性配準，使 K 層向 K1 層變形，得到一個新層 K2，根據光流方程得到位移 μ：
|
其中，f 和 r 分別是浮動圖像 F 和參考圖像 R 在像素點 p 處的灰度值，▽r 是參考圖像在 p 點處的灰度梯度。最后，計算 K1、K2 層與 K? 1、K+1 層之間的 MI 值，用 MI 最大值所在層（K1 或 K2）替換原始數據中的 K 層。如圖 5 所示，運動校正前血管圖像中存在鋸齒狀的運動偽影，校正后圖像的清晰度明顯提高。該方法不需要動態參考對象，并且只需一組三維光聲信號數據即可，但是數據采集和處理速度較慢。利用浸水式微電機掃描系統可提高數據采集速度[19]，此外利用高性能圖形處理單元（graphic processing unit，GPU）進行并行計算，也可以縮短數據處理時間，實現實時成像。
1.4 基于深度學習的運動校正
近年來，深度學習（deep learning，DL）特別是卷積神經網絡（convolutional neural network，CNN）已經迅速成為醫學圖像分析與處理的首選方法，在快速重建高質量的圖像方面也表現出了巨大潛力[20]。例如 Chen 等[21]采用 DL 方法實現了 OR-PAM 圖像中運動偽影的校正，其原理是：首先利用希爾伯特變換提取每個深度方向上光聲信號的包絡，并取其軸向投影的最大幅值，獲得最大幅值投影（maximum amplitude projection，MAP）圖像；然后搭建 CNN 模型，根據先驗知識設計損失函數，使用反向傳播算法更新網絡參數（例如卷積核和神經元偏移向量），通過最小化損失函數對模型進行訓練，提升模型的學習能力；最后利用該 CNN 模型進行圖像后處理，校正運動偽影，結果如圖 6 所示。增大卷積核可提高 CNN 的性能，但會降低處理效率，因此卷積核大小的選擇應是性能與效率之間權衡的結果。雖然該方法是為 OR-PAM 設計的，但是可以通過改變訓練集實現其他成像模式的運動偽影校正。然而 CNN 模型的搭建及訓練需要大量成像數據，到目前為止還沒有針對光聲圖像重建的公開數據集，可用的臨床病例數據也十分有限。
2 PAT 運動校正方法
2.1 PAT 原理及目標運動的影響
如圖 7 所示，PAT 是用非聚焦短脈沖激光作為激勵源照射待成像組織，利用具有掃描功能的超聲換能器或換能器陣列采集組織產生的光聲信號，組織內不同深度的光吸收體產生的光聲信號到達超聲換能器的時間不同，利用時間分辨技術可獲得多個不同層析面的光聲信號，進而重建組織的二維或三維初始聲壓分布圖或光吸收能量分布圖[22]。PAT 的成像深度可達數厘米，空間分辨率可達數十微米，在冠狀動脈粥樣硬化性疾病、乳腺癌和消化系統病變等疾病的早期診斷方面具有良好的應用前景[23]。
圖7
PAT 成像原理示意圖
Figure7.
Schematic diagram of PAT principle
由于 PAT 具有在單個納秒級激光脈沖作用下從整個成像區域采集圖像數據的能力，因此它能夠有效消除單幀圖像中的運動偽影。但是，在多光譜光聲層析成像（multi-spectral optoacoustic tomography，MSOT）或基于掃描的成像系統中，多幀圖像分析通常是必不可少的，因此需要采用有效的方法校正由運動引起的圖像失真。例如，可以減少連續激光脈沖之間的延遲抑制運動[24-25]，但需要專用的激光器系統。除此之外，主要的運動校正方法包括呼吸門控技術、運動聚類法、基于混合成像模型的方法和幀平均-幀間運動補償法，如表 2 所列。
表2
PAT 運動校正方法的對比
Table2.
Comparison of existing methods of motion correction for PAT
2.2 呼吸門控技術
門控技術是醫學成像中減少由心跳和呼吸所致運動偽影的常用方法，分為前瞻性和回顧性兩種方式。前瞻性門控是指采用專門的心電或呼吸門控采集裝置，只在每個心動或呼吸周期中的相同相位處采集圖像或信號，由于每個心動或呼吸周期只采集一幀圖像，因而與連續采集的方式相比，這種方式延長了成像時間。回顧性門控是指在整個心動或呼吸周期內連續采集圖像或信號，同時記錄心電信號或呼吸波形。待掃描結束后，將采集的圖像或信號與心電或呼吸信號進行對照和分析，選擇在相同相位處采集的圖像或信號，該技術在許多需要精確定位器官的成像中是必不可少的[26-27]。
小動物全身 PAT 系統的幀頻通常小于 1 Hz，而小鼠的呼吸頻率可高達 3 Hz，所以在對胸部以下區域進行 PAT 掃描時，呼吸是產生運動偽影的主要原因[28]，導致圖像模糊和分辨率下降。回顧性呼吸門控技術是消除 PAT 呼吸運動偽影的常用方法，包括硬件門控和軟件門控兩種方式。回顧性硬件呼吸門控是在成像目標自由呼吸時采集光聲信號，同時利用外部設備監測呼吸波形，通過分析光聲信號數據，選擇出具有相同呼吸相位的信號[27]。例如，Xia 等[29]利用回顧性呼吸門控技術對小鼠的肝臟和腎臟 PAT 圖像進行運動校正，將小鼠完全或部分浸入水中，利用壓力傳感器持續監測水位波動，并據此監測呼吸波形。然后，根據呼吸相位對光聲信號進行排序和聚類，采用在同一呼吸相位采集的光聲數據重建圖像，結果如圖 8 所示。該方法對呼吸波形的監測不受圖像噪聲的影響，而且數據處理的計算量較低，可直接觀察呼吸周期，提高瞬時分辨率。但是若目標呼吸比較急促或者不均勻，往往會出現漏感應或錯誤感應，導致觸發精度降低，因此需要對目標進行呼吸訓練。
回顧性軟件呼吸門控不需要借助外部硬件設備監測呼吸波形，而是根據采集的光聲信號計算呼吸門控信號[30]。例如，Ron 等[28]提出一種自門控技術，首先將超聲換能器陣列中各陣元采集的光聲信號排列成二維矩陣，其中行表示各信號的 m 個時間樣本，列對應換能器陣列的 n 個陣元，在每個換能器位置共采集 k 幀，并將該幀序列排列成一個 m × n 行和 k 列的二維矩陣，表示在該換能器位置的整個幀序列；然后利用主成分分析法計算相鄰幀之間的自相關矩陣，并對其采用二階 k-means 算法進行聚類[31-32]，將 k 幀分為運動幀和靜態幀兩類；最后利用濾波反投影算法重建三維圖像，結果如圖 9 所示。該方法可在重建圖像之前進行運動校正，提高圖像的清晰度和信噪比，避免對圖像進行后處理。但是有的生物體兩次呼吸之間的停頓超過單次呼吸長度，需根據呼吸特征的先驗知識建立運動幀和靜態幀之間的區分標準。此外，呼吸門控信號提取的準確性會影響運動校正的效果。
圖9
采用自呼吸門控技術對小鼠 PAT 圖像的運動偽影抑制結果[28]
a. 選擇出的幀；b. 丟棄的幀（箭頭指示產生變形的結構）；c. 圖 a 中虛線標記位置的聲壓幅值曲線，其中虛線曲線對應采用所有幀進行圖像重建，實線曲線對應采用選擇幀進行圖像重建
Figure9. Motion rejection results for spiral volumetric optoacoustic tomography[28]a. selected frame; b. rejected frame (distorted structures are marked with arrows); c. amplitude profiles marked in Fig. a (dashed lines) for images reconstructed from all the frames (dashed lines) versus selected frames (solid lines)
2.3 基于運動聚類的校正方法
MSOT 根據不同組織的光譜吸收特性，利用不同波長的脈沖激光照射組織，產生具有組織特異性的光聲信號，可更好地識別和標定組織成分[33]。對于心臟 MSOT，信號采集過程中的心跳運動會導致圖像產生模糊，分辨率降低[34-35]。Taruttis 等[34]提出采用運動聚類技術對心臟 MSOT 圖像進行運動校正，其原理是：首先利用反投影算法重建每個波長對應的 100 幀圖像；然后利用 k-means 算法對 100 幀圖像進行聚類，再對聚類后的圖像求平均，考慮到心室運動存在兩個轉折點，因此聚類的集群數量取值為 2；最后對平均后的圖像進行線性光譜解混。如圖 10 所示，執行聚類的圖像中由運動造成的模糊程度明顯降低，且增加集群數量可進一步降低運動模糊程度，但也會降低信噪比。隨著集群數量的增加，集群之間的距離會減小，也提高了對 k-means 算法聚類能力的要求，因此需要結合圖像模糊程度、信噪比和 k-means 算法的性能設定集群數量。
圖10
由光譜解混獲得的脫氧血紅蛋白 MSOT 圖像[34]
a. 直接對各波長的 100 幀圖像進行平均得到的光譜解混圖像；b. 只對各波長的較大集群圖像進行平均得到的光譜解混圖像；c. 只對各波長的較小集群圖像進行平均得到的光譜解混圖像
Figure10. Motion correction for unmixed MSOT images of deoxyhemoglobin contribution[34]a. unmixed image after all 100 images were averaged for each wavelength; b. unmixed image after averaging only the larger clusters for each wavelength; c. unmixed image after averaging only the smaller clusters for each wavelength
2.4 基于混合成像模型的運動校正
Chung 等[36]提出一種基于混合成像模型的圖像重建方法，建立考慮呼吸運動的成像模型，并求解使模型輸出值與光聲測量數據之間誤差最小的非線性最小二乘問題，同時實現運動校正和圖像重建。其原理是：假設聲速恒定，成像目標在連續兩次照明和測量之間的變化非常緩慢，用聲壓的圓 Radon 變換（circular radon transform，CRT）[37]表示聲壓測量數據：
|
式中，r 是換能器與變形區域之間的距離，? 是換能器陣列旋轉的角度，z（?）=（cos?, sin?）是測量位置，f? 是初始聲壓 f 在 t? 時刻的變形，Φ（?, x）是點 x 到點 y 的變形映射，Ψ（?,·）是 Φ（?,·）的逆映射[38]，Jx［Φ（?, x）］是變量 x→y = Φ（?, x）變形映射的雅可比行列式的絕對值。那么，PAT 圖像重建問題可簡化為逆 CRT 問題。假設由呼吸產生的運動是垂直拉伸，建立光聲成像的離散模型：
|
其中，g =［g1 ··· gn］T∈Rn × m 是測量值矩陣，n 是換能器的測量位置數，f 是初始聲壓矩陣，γ =［γ1 ··· γn］T 是運動參數矩陣，γi 是 zi 處的運動參數，e =［e1 ··· en］T∈Rn × m 是獨立同分布的高斯噪聲矩陣，A（γ）=［A1K（γ1） ··· AnK（γn）］T∈Rnm × N，Ai 是每個換能器位置對應的投影矩陣，Aif 表示 f 在以 zi 為中心的圓上的 CRT，fi = K（γi）f 是第 i 個變形后的聲壓矩陣，K（γi）∈RN × N 表示一維拉伸/膨脹變形。當運動參數是變量時，f 的估計問題定義為如下非線性最小二乘問題：
|
式中，λ 是正則化參數。利用變量投影算法[39]分離線性變量和非線性變量，進而采用高斯牛頓迭代法[40]對非線性參數γ進行優化，其中利用復合最小二乘 QR 分解（hybrid least square QR factorization，HyBR）算法自動選擇正則化參數，并求解關于 f 的線性最小二乘問題，結果如圖 11 所示。該方法不僅可應用于垂直拉伸模型，還可以擴展到其他變形模型，例如旋轉或者平移。
圖11
基于混合成像模型的 PAT 圖像重建結果[36]
a. 忽略運動；b. 采用 LSQR 算法求解線性問題；c. 采用 HyBR 算法求解線性問題
Figure11. PAT images reconstructed based on a hybrid model[36]a. initial regularized reconstruction ignoring motion; b. reconstruction obtained by using LSQR to solve the linear problem; c. reconstruction obtained by using HyBR to solve the linear problem
2.5 幀平均-幀間運動補償
對于深層組織的 PAT，由于產生的聲壓信號較弱，導致重建圖像質量較差。為了提高信噪比，可以利用高功率脈沖激光來提高聲壓，但是這種方法存在安全風險（皮膚輻射暴露 < 20 mJ/cm2）。對于采用安全激光脈沖能量的 PAT，可利用幀平均（frame averaging，FA）技術，對相同幾何結構的圖像進行平均以提高圖像的信噪比。但是因為激光脈沖重復頻率的限制以及成像目標運動的影響，采用 FA 技術很難實現高分辨率實時成像。針對該問題，Kim 等[41]將 FA 技術與幀間運動補償（inter-frame motion compensation，IFMC）相結合，即 FA-IFMC 技術，用于補償深層組織 PAT 圖像中的運動偽影。其原理是：首先，利用塊匹配算法[42]和三步搜索法進行運動估計，從最近獲得的圖像中選取局部參考區，在先驗圖像中選取對比區，二者大小相同，用平均絕對誤差（mean absolute difference，MAD）作為匹配準則：
|
其中，N 和 M 分別是待匹配圖像塊在軸向和橫向上的像素數，Cij 和 Rij 分別是對比區和參考區中像素（i, j）的灰度值。將參考區及其鄰近區域定義為搜索區域，在參考區及其八鄰域方向上以步長逐次減半的方式進行搜索，在搜索區域中找到與對比區最相似的子區域，參考區和相似子區域之間的坐標偏移就是估計出的運動矢量。然后，根據運動矢量，利用 IFMC 比較圖像序列中的相鄰幀實現運動偽影的校正。
如圖 12 所示，與傳統 FA 相比，FA-IFMC 提高了成像的峰值信噪比和軸向分辨率。與幀運動補償（frame motion compensation，FMC）[41]相比，由于 FMC 選取的參考圖像和對比圖像不是相鄰幀，兩幀之間的運動矢量隨著時間的增加而積累，而搜索步長又限制了搜索區域的大小，導致運動矢量很容易超出搜索范圍，而 FA-IFMC 的參考圖像和對比圖像取自相鄰幀，避免了幀平均-幀運動補償（frame averaging-frame motion compensation，FA-FMC）中的問題。但是步長決定了運動估計的精度和搜索范圍，需謹慎選擇，后續研究中可利用變形塊匹配算法進行運動估計，并利用浮動步長提高運動補償的能力。
3 結論
在光聲成像中，目標的呼吸、心跳、脈搏、眼球轉動和身體移動等會導致圖像中產生運動偽影，降低成像分辨率，對運動偽影進行抑制和補償對于提高成像質量具有重要意義。本文分別對 PAM 和 PAT 中運動偽影的主要校正方法進行了歸納和總結。目前，該領域仍面臨諸多問題，例如：為了提高運動聚類法的準確性，需要不斷提高 k-means 算法的聚類性能；多數方法不具有一般性，例如基于 Demons 圖像配準的方法只有在成像區域中血管數量較多的情況下才能取得理想的校正效果，呼吸門控法只適用于校正由呼吸運動產生的偽影；多數方法只能校正垂直方向上由一維運動產生的偽影，針對其他方向上或者由二維運動產生的偽影的校正是未來主要的研究方向；雖然采用運動校正方法可明顯改善圖像質量，但是為了實現實時成像，仍需提高成像系統的數據采集和處理速度。
此外，現有的運動校正方法均是針對體外掃描成像系統的，尚未檢索到針對內窺式光聲成像[43]的運動校正方法。例如，冠狀動脈血管內光聲（intravascular photoacoustic，IVPA）圖像序列[44]中存在由周期性的心臟運動、呼吸運動和動脈內搏動的血流造成的運動偽影，在橫截面圖像中表現為相鄰幀之間血管橫截面的錯位，在縱向視圖中表現為血管壁呈現鋸齒狀邊緣，嚴重影響了 IVPA 圖像序列的視覺效果，降低了血管形態結構和組織成分定量分析以及血管三維重建的精度。因此，針對內窺式光聲成像中運動偽影的校正將是未來該領域的重要研究方向之一。
利益沖突聲明：本文全體作者均聲明不存在利益沖突。
