醫學光聲成像的物理原理與仿真實現

薛澤鵬 徐家琦 劉炳坤 衛子瑤 張本健 林輝

摘要：光聲成像利用生物組織的光聲效應進行成像，具有分辨率高、對比度強、穿透性深的優勢，在腫瘤醫學檢測方面具有巨大發展潛力。本工作基于熱聲轉換原理與軟件，對圓盤源、圖片源及具有高含氧量腫瘤區域的腦部CT圖片源的圖像重建進行了仿真實現，對比了快速傅里葉變換和時間反轉重建算法。研究工作表明兩種圖像重建算法均可重建復雜圖像源。

關鍵詞：光聲成像;時間反轉;快速傅里葉變換;含氧量

中圖分類號：TP18? ? ? ? 文獻標識碼：A

文章編號：1009-3044（2019）23-0277-03

開放科學（資源服務）標識碼（OSID）：

光聲層析成像（Photoacoustic Tomography）是基于光聲效應發展而來的一種新型的無損成像技術。通過短波激光作為光源照射物體，在激勵激光脈沖的持續時間遠小于物體對激光的吸收和傳導時間時，物體內部會發生熱膨脹，進而發出球面的超聲波[1]。醫學光聲成像根據生物組織對照射激光吸收的聲壓分布情況對組織結構進行成像。當生物組織受到短脈沖電磁照射導致瞬間熱膨脹，從而發射相應頻率和幅值的超聲波信號，然后用時間反轉算法或快速傅里葉變換算法重建圖像，將生物組織內的微波能量吸收分布重建出來。

通過對醫學光聲圖像的分析可用于對患者腫瘤等疾病的診斷[2]。由于癌組織內不同的化學成分具有不同的光學吸收特性，光聲成像利用光致超聲的成像機制，選擇特定的激光波長對癌組織內化學成分進行檢測，從而對癌組織進行高對比度的結構成像。腫瘤細胞的一個重要標志是過量的氧燃燒（即高代謝），高代謝導致組織內血紅蛋白濃度、血氧濃度和氧代謝率等發生改變[3]。光聲成像可對癌組織的這些參數進行定量分析，從而實現癌組織功能成像[4]。因此光聲成像技術在軟組織腫瘤的早期無損診斷、定位及分析上是行之有效的檢測方法，在未來腫瘤的臨床應用上具有潛在的應用價值。

本文首先介紹了光聲成像的基本原理和重建算法，并采用改編的k-Wave 軟件對圓盤源、圖片源及具有高含氧量腫瘤區域的腦部CT圖片源進行了圖像重建，對比了快速傅里葉變換（Fast-Fourier-Transform， FFT）和時間反轉（Time-Reversal， TR）重建算法各自的優勢。

1 醫學光聲成像原理

生物組織受到短脈沖電磁照射會瞬時吸熱膨脹，從而發射相應頻率和幅值的超聲波信號，通過超聲傳感器掃描收集樣品發出的熱聲信號，用圖像重建算法就可以把生物組織的微波能量吸收分布重建出來。熱聲波動方程（1）是熱聲成像的基本依據：

式中：p（r，t） 是時刻t位置r處的聲壓;為體膨脹系數;T（r， t） 是媒質由于吸收微波能量而引起的溫升;和分別為媒質的密度和聲速。

假設在一個小范圍內微波能量沉積時間很短，以至于可以忽略熱擴散，則可以將溫度變化描述為

式中：

將式（2）～（4）代入方程（1）中，得

熱聲或光聲成像本質上就是利用探測到的熱聲信號p（r， t） 來反演出生物組織中的微波能量吸收分布（SAR）。

2 計算機仿真方法

2.1 k-Wave軟件

k-Wave是用于Matlab和C++的開源聲學工具箱，由Bradley Treeby和Ben Cox（倫敦大學學院University College London）和Jiri Jaros（布爾諾科技大學Brno University of Technology）開發。該軟件設計用于復雜和組織真實介質中的時域聲學和超聲模擬。仿真功能基于k空間偽譜方法，具有快速、易用的特點。該工具箱包括能夠模擬聲波傳播壓力和速度源進行建模的高級時域模型，優化c++代碼，最大限度地提高大型仿真計算的性能[5]。

2.2 腫瘤組織的光吸收處理

k-Wave仿真函數主要解決的是超聲場的一階方程組，分別為動量守恒、質量守恒、聲壓與密度的關系函數。通過分析k-Wave基本方程，可以加深對k-Wave工具箱工作原理的理解，便于進行編程實現相關的模擬仿真。

在運行過程中，由于網格劃分問題容易使運行繁瑣復雜，同時為了解決數值分散問題，需要設置小時間步長。由于偽譜法采用了有限差分計算方法，并吸收了傅里葉級數的思想，經分析K-space算法能夠解決這一問題[6]。采用傅里葉級數具有兩大優勢，即通過傅里葉變換可以精確計算傅里葉分量的幅值，且由于采用正弦曲線作為基礎函數而可以使用更少的節點來保證計算精度。

結合上述相關可視化算法和理論，基于k-Wave仿真工具箱，以工具箱提供的2D-Time- Reversal-Reconstruction-For-A-Linear-Sensor例子為基礎，對其進行了代碼改編和模型植入。研究了快速傅里葉變換和時間反轉兩種圖像重建算法對圓盤源、圖片源，以及具有高含氧量腫瘤區域的腦部CT圖片源的重建效果。

3 算例與結果

3.1圓盤源

設水中有兩個圓盤聲源同時傳播，在空間相遇并相互穿越，最后聲信號被置于水模體邊緣的線性傳感器探測到。兩個圓盤波源半徑分別為2.5mm 和 4mm， 在2D平面上相距42mm。

（b） FFT算法重建圖

（c） TR算法重建圖 ]

圖1給出了圓盤源的原圖與重建圖，可見TR算法重建圖（圖1（c）） 相對于FFT算法重建圖（圖1（b）） 出現了更多的偽影，需要后期加工去除。在重建速度方面，FFT算法是使用傅里葉變化的公式進行運算，而TR算法相當于把傳播到水模體邊緣的聲信號作為新的波源，重新進行一次傳播仿真，因此使得TR算法的時間幾乎是FFT時間的兩倍。對于圖像尺寸較小，或圖像較為對稱（例如點波源、圓盤源），FFT算法較TR算法具有優勢。

3.2 圖片波源

為校驗程序對離散的圖片波源的重建效果，本工作還在原程序中植入了圖片波源的讀入代碼，可讀入jpg/png/bmp 格式的圖片，并將圖片網格的灰度值作為各個離散波源的強度值， 重建了二維小熊貓圖像。

圖2 給出了小熊貓原圖與重建圖，圖2（b）為FFT算法重建的圖像，圖2（c）為TR算法重建的圖像。可見，圖2（b）相比于圖2（c）出現了更多的偽影。因此對于圖片源，即在組織內部為非均勻聲強的情況下，TR算法在圖像重建方面性能較優。因為相對其他重建算法，TR算法的假定初始條件更少，不需要設定聲介質的初始速度。TR算法可以對任意形狀的封閉表面進行測量，不受測量表面外的聲源影響，甚至在一些假設條件不成立如非均勻的聲速分布時，該算法仍可以得到較為理想的結果，因此TR算法被認為是具備“最少約束條件”的光聲圖像重建算法[7]。

3.3 組織氧含量仿真

光聲信號強弱是由光能的吸收和分布決定的。血紅蛋白存在于生物組織中，氧的輸運需要通過血紅蛋白。與氧結合的血紅蛋白稱為氧合血紅蛋白，表現為紅色;而脫氧的血紅蛋白呈藍紫色。血紅蛋白的含氧量與生物組織的代謝密切相關，在不同的組織中具有不同的氧飽和度。腫瘤區域含氧量高，吸收光能也相對較高。

我們通過改變圖片區域的顏色來模擬含氧量的差別，即使用Matlab通過更改圖片的RGB256位數值將圖片部分區域的顏色更改為藍色和紅色來模擬對光吸收能力的高低。

在原始圖像中，由于背景（深色）吸收光的能力要比目標區域高得多。而背景產生的光聲信號先于目標區域的信號到達傳感器，導致目標區域的光聲信號被背景區域的信號所覆蓋，導致重建圖像嚴重失真。為了避免這個問題，我們在將原始圖像輸入代碼之前，添加了一個灰度反轉程序。

灰度變換是圖像增強的一種經典而有效的方法。其原理是將圖像中的每一個像素的灰度值，通過一個函數對應到另一個灰度值上去，從而實現灰度的變換。常見的灰度變換有線性灰度變換和非線性灰度變換。對灰度圖像進行線性灰度變換能將輸入圖像的灰度值的動態范圍按線性關系公式拉伸擴展至指定范圍或整個動態范圍。

我們對一腦部腫瘤CT圖像進行了圖像重建。在模擬重建過程中，利用Matlab 選中原圖黑色區域（rgb（0，0，0））替換為藍色（rgb（0，0，255））（圖中為深灰色）（圖3（a））和紅色（rgb（255，0，0））（圖中為淺灰色）（圖3（b）），利用灰度值的差異模擬藍色、紅色兩種顏色對光能的吸收強弱區別，來模擬低氧濃度和高氧濃度兩種情況，然后進行TR算法重建得到圖3（c）和（d），觀察氧濃度對重建的影響。

可見，淺色組的成像內容多在黑區的一部分，這也證明在光聲成像中，光吸收能較強的區域可以更明顯地顯現出來。因此由于腫瘤細胞含氧量較高，腫瘤區域也可以通過模型很好地重建成像。淺色組重建圖像中腫瘤區域部分灰度值要高于深灰色組。重建圖像的灰度值與組織氧濃度水平之間呈非線性關系。

4 結論

光聲成像利用生物組織的光聲效應進行成像，在醫學癌癥檢測方面具有巨大發展潛力。本工作基于熱聲轉換原理與Matlab工具箱k-Wave軟件，對圓盤源、圖片源及具有高含氧量腫瘤區域的腦部CT圖片源的圖像重建進行了仿真實現。對比了快速傅里葉變換和時間反轉重建算法的優缺點。結果表明兩種圖像重建算法均可重建復雜圖像源，其中對于圖像尺寸較小，或圖像較為對稱（例如點波源、圓盤源），FFT算法較TR算法具有優勢。而對于復雜的非均勻強度源TR算法的精度更高。

致謝：本工作得到合肥工業大學大學生創新訓練計劃項目（編號：2018CXCYS209）的支持。

參考文獻：

[1] Bell A G. On the production and reproduction on of sound by light[J]. American Journal of Science， 1880（118）： 305-324.

[2] Fakhrejahani E， Torii M， Kitai T， et al. Clinical Report on the First Prototype of a Photoacoustic Tomography System with Dual Illumination for Breast Cancer Imaging[J]. Plos? One， 2015， 10（10）：e0139113.

[3] Gottschalk S， Fehm T F， Deanben X L， et al. Noninvasive real-time visualization of multiple cerebral hemodynamic parameters in whole mouse brains using five-dimensional optoacoustic tomography[J]. Journal of Cerebral Blood Flow & Metabolism Official Journal of the International Sociely of Cerebral Blood Flow & Metabolism， 2015， 35（4）： 531-535.

[4] Wang L， Maslov K， Xing W， et al. Video-rate functional photoacoustic microscopy at depths[J]. Journal of Biomedical Optics， 2012， 17（10）： 1060071-1060075.

[5] Treeby BE， Cox BT. k-Wave： MATLAB toolbox for the simulation and reconstruction of photoacoustic wave fields. J. Biomed. Opt. 2010; 15（2）：021314.

[6] 胡文娟.基于k空間的并行磁共振成像算法的研究[D].東北大學.2011

[7] 孫明健. 光聲成像系統信號采集與圖像重建研究 [D]. 哈爾濱：哈爾濱工業大學，2011：3-5.

【通聯編輯：唐一東】