薛澤鵬 徐家琦 劉炳坤 衛(wèi)子瑤 張本健 林輝

摘要：光聲成像利用生物組織的光聲效應(yīng)進(jìn)行成像，具有分辨率高、對(duì)比度強(qiáng)、穿透性深的優(yōu)勢(shì)，在腫瘤醫(yī)學(xué)檢測(cè)方面具有巨大發(fā)展?jié)摿Α１竟ぷ骰跓崧曓D(zhuǎn)換原理與軟件，對(duì)圓盤源、圖片源及具有高含氧量腫瘤區(qū)域的腦部CT圖片源的圖像重建進(jìn)行了仿真實(shí)現(xiàn)，對(duì)比了快速傅里葉變換和時(shí)間反轉(zhuǎn)重建算法。研究工作表明兩種圖像重建算法均可重建復(fù)雜圖像源。

關(guān)鍵詞：光聲成像;時(shí)間反轉(zhuǎn);快速傅里葉變換;含氧量

中圖分類號(hào)：TP18? ? ? ? 文獻(xiàn)標(biāo)識(shí)碼：A

文章編號(hào)：1009-3044（2019）23-0277-03

開放科學(xué)（資源服務(wù)）標(biāo)識(shí)碼（OSID）：

光聲層析成像（Photoacoustic Tomography）是基于光聲效應(yīng)發(fā)展而來的一種新型的無損成像技術(shù)。通過短波激光作為光源照射物體，在激勵(lì)激光脈沖的持續(xù)時(shí)間遠(yuǎn)小于物體對(duì)激光的吸收和傳導(dǎo)時(shí)間時(shí)，物體內(nèi)部會(huì)發(fā)生熱膨脹，進(jìn)而發(fā)出球面的超聲波[1]。醫(yī)學(xué)光聲成像根據(jù)生物組織對(duì)照射激光吸收的聲壓分布情況對(duì)組織結(jié)構(gòu)進(jìn)行成像。當(dāng)生物組織受到短脈沖電磁照射導(dǎo)致瞬間熱膨脹，從而發(fā)射相應(yīng)頻率和幅值的超聲波信號(hào)，然后用時(shí)間反轉(zhuǎn)算法或快速傅里葉變換算法重建圖像，將生物組織內(nèi)的微波能量吸收分布重建出來。

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

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

1 醫(yī)學(xué)光聲成像原理

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

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

假設(shè)在一個(gè)小范圍內(nèi)微波能量沉積時(shí)間很短，以至于可以忽略熱擴(kuò)散，則可以將溫度變化描述為

式中：

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

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

2 計(jì)算機(jī)仿真方法

2.1 k-Wave軟件

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

2.2 腫瘤組織的光吸收處理

k-Wave仿真函數(shù)主要解決的是超聲場(chǎng)的一階方程組，分別為動(dòng)量守恒、質(zhì)量守恒、聲壓與密度的關(guān)系函數(shù)。通過分析k-Wave基本方程，可以加深對(duì)k-Wave工具箱工作原理的理解，便于進(jìn)行編程實(shí)現(xiàn)相關(guān)的模擬仿真。

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

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

3 算例與結(jié)果

3.1圓盤源

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

（b） FFT算法重建圖

（c） TR算法重建圖 ]

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

3.2 圖片波源

為校驗(yàn)程序?qū)﹄x散的圖片波源的重建效果，本工作還在原程序中植入了圖片波源的讀入代碼，可讀入jpg/png/bmp 格式的圖片，并將圖片網(wǎng)格的灰度值作為各個(gè)離散波源的強(qiáng)度值， 重建了二維小熊貓圖像。

圖2 給出了小熊貓?jiān)瓐D與重建圖，圖2（b）為FFT算法重建的圖像，圖2（c）為TR算法重建的圖像?？梢姡瑘D2（b）相比于圖2（c）出現(xiàn)了更多的偽影。因此對(duì)于圖片源，即在組織內(nèi)部為非均勻聲強(qiáng)的情況下，TR算法在圖像重建方面性能較優(yōu)。因?yàn)橄鄬?duì)其他重建算法，TR算法的假定初始條件更少，不需要設(shè)定聲介質(zhì)的初始速度。TR算法可以對(duì)任意形狀的封閉表面進(jìn)行測(cè)量，不受測(cè)量表面外的聲源影響，甚至在一些假設(shè)條件不成立如非均勻的聲速分布時(shí)，該算法仍可以得到較為理想的結(jié)果，因此TR算法被認(rèn)為是具備“最少約束條件”的光聲圖像重建算法[7]。

3.3 組織氧含量仿真

光聲信號(hào)強(qiáng)弱是由光能的吸收和分布決定的。血紅蛋白存在于生物組織中，氧的輸運(yùn)需要通過血紅蛋白。與氧結(jié)合的血紅蛋白稱為氧合血紅蛋白，表現(xiàn)為紅色;而脫氧的血紅蛋白呈藍(lán)紫色。血紅蛋白的含氧量與生物組織的代謝密切相關(guān)，在不同的組織中具有不同的氧飽和度。腫瘤區(qū)域含氧量高，吸收光能也相對(duì)較高。

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

在原始圖像中，由于背景（深色）吸收光的能力要比目標(biāo)區(qū)域高得多。而背景產(chǎn)生的光聲信號(hào)先于目標(biāo)區(qū)域的信號(hào)到達(dá)傳感器，導(dǎo)致目標(biāo)區(qū)域的光聲信號(hào)被背景區(qū)域的信號(hào)所覆蓋，導(dǎo)致重建圖像嚴(yán)重失真。為了避免這個(gè)問題，我們?cè)趯⒃紙D像輸入代碼之前，添加了一個(gè)灰度反轉(zhuǎn)程序。

灰度變換是圖像增強(qiáng)的一種經(jīng)典而有效的方法。其原理是將圖像中的每一個(gè)像素的灰度值，通過一個(gè)函數(shù)對(duì)應(yīng)到另一個(gè)灰度值上去，從而實(shí)現(xiàn)灰度的變換。常見的灰度變換有線性灰度變換和非線性灰度變換。對(duì)灰度圖像進(jìn)行線性灰度變換能將輸入圖像的灰度值的動(dòng)態(tài)范圍按線性關(guān)系公式拉伸擴(kuò)展至指定范圍或整個(gè)動(dòng)態(tài)范圍。

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

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

4 結(jié)論

光聲成像利用生物組織的光聲效應(yīng)進(jìn)行成像，在醫(yī)學(xué)癌癥檢測(cè)方面具有巨大發(fā)展?jié)摿?。本工作基于熱聲轉(zhuǎn)換原理與Matlab工具箱k-Wave軟件，對(duì)圓盤源、圖片源及具有高含氧量腫瘤區(qū)域的腦部CT圖片源的圖像重建進(jìn)行了仿真實(shí)現(xiàn)。對(duì)比了快速傅里葉變換和時(shí)間反轉(zhuǎn)重建算法的優(yōu)缺點(diǎn)。結(jié)果表明兩種圖像重建算法均可重建復(fù)雜圖像源，其中對(duì)于圖像尺寸較小，或圖像較為對(duì)稱（例如點(diǎn)波源、圓盤源），F(xiàn)FT算法較TR算法具有優(yōu)勢(shì)。而對(duì)于復(fù)雜的非均勻強(qiáng)度源TR算法的精度更高。

致謝：本工作得到合肥工業(yè)大學(xué)大學(xué)生創(chuàng)新訓(xùn)練計(jì)劃項(xiàng)目（編號(hào)：2018CXCYS209）的支持。

參考文獻(xiàn)：

[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].東北大學(xué).2011

[7] 孫明健. 光聲成像系統(tǒng)信號(hào)采集與圖像重建研究 [D]. 哈爾濱：哈爾濱工業(yè)大學(xué)，2011：3-5.

【通聯(lián)編輯：唐一東】