柳芳艷 孟靜 司廣濤

摘 要：針對(duì)限制視圖下光聲圖像的重建偽影問(wèn)題，提出一種改進(jìn)的固定點(diǎn)迭代量化光聲成像方法。首先，通過(guò)傳統(tǒng)的反投影重建算法重建由探測(cè)器探測(cè)到的原始光聲壓數(shù)據(jù)，得到原始的光聲壓圖像;接著，利用自適應(yīng)維納濾波算法對(duì)原始的光聲壓圖像進(jìn)行濾波去除重建圖像偽影;然后，通過(guò)光傳輸模型求解目標(biāo)成像區(qū)域的光通量;最后，進(jìn)行迭代計(jì)算，獲得目標(biāo)組織的光學(xué)吸收系數(shù)。此外，在求解光通量過(guò)程中引入Toast++軟件來(lái)實(shí)現(xiàn)光傳輸模型的前向求解，

提高量化成像的效率和精確性。仿體和活體實(shí)驗(yàn)結(jié)果表明，與傳統(tǒng)固定點(diǎn)迭代方法相比，所提方法能夠獲取更高質(zhì)量的光聲圖像，重建得到的深層量化光聲圖像中存在較少偽影;量化重建的深層目標(biāo)組織的光學(xué)吸收系數(shù)與淺層目標(biāo)組織的光學(xué)吸收系數(shù)的數(shù)值非常接近，前者約為后者的70%，能實(shí)現(xiàn)深層生物組織光學(xué)吸收系數(shù)的定量重建。

關(guān)鍵詞：量化光聲成像;光學(xué)吸收系數(shù);深層組織;維納濾波;Toast++;光子傳輸模型

中圖分類號(hào)：TP391.41

文獻(xiàn)標(biāo)志碼：A

Abstract： Focusing on the reconstruction artifact of photoacoustic images in restricted view， an improved fixed-point iterative quantitative photoacoustic imaging method was proposed. Firstly， the original photoacoustic pressure data detected by the detector were reconstructed by the traditional back projection reconstruction algorithm to obtain the original photoacoustic pressure image. Secondly， the original photoacoustic pressure image was filtered to remove the reconstruction artifact by adaptive Wiener filtering algorithm. Thirdly， the optical transmission model was used to solve the optical flux of the target imaging region. And finally， iterative calculation was performed to obtain the optical absorption coefficient of the target tissue. In addition， Toast++ software was introduced in the process of solving the optical flux to realize the forward solution of the optical transmission model， which improved the efficiency and accuracy of quantitative imaging. The phantom and in vivo experiments show that compared with the traditional fixed-point iterative method， the proposed method can obtain photoacoustic images with higher quality and there are fewer artifacts in the deep quantitative photoacoustic images reconstructed by the method. The optical absorption coefficient of the quantitatively reconstructed deep target tissue is very close to the optical absorption coefficient of the shallow target tissue， the former is about 70% of the latter. As a result， the quantitative reconstruction of the optical absorption coefficient of the deep biological tissue can be implemented by the proposed method.

Key words： quantitative photoacoustic imaging; optical absorption coefficient; deep tissue; Wiener filtering; Toast++; photon transmission model

0 引言

光聲成像（PhotoAcoustic Imaging， PAI）是近幾年發(fā)展起來(lái)的一種非入侵式和非電離式的新型生物醫(yī)學(xué)成像方式[1]，它融合了光學(xué)和超聲兩種成像模式的特點(diǎn)，兼具純光學(xué)成像的高對(duì)比度和純超聲成像的高分辨率優(yōu)點(diǎn)，已成為目前發(fā)展最快的生物醫(yī)學(xué)成像技術(shù)之一。光聲成像作為一種重要的生物醫(yī)學(xué)成像方式，為研究組織的結(jié)構(gòu)形態(tài)，生理功能和新陳代謝指標(biāo)提供了重要手段，在癌癥和心腦血管疾病的早期診斷中具有重要的應(yīng)用價(jià)值[2-3]。光聲成像的基本原理是：生物組織經(jīng)過(guò)合適的短脈沖激光照射后，吸收光能，組織內(nèi)部溫度瞬時(shí)升高，膨脹產(chǎn)生壓力，這種壓力以超聲波的方式在組織內(nèi)部傳播，最后被放置在生物組織表面的超聲換能器探測(cè)。根據(jù)采集到的超聲信號(hào)，經(jīng)過(guò)特定的重建方法可以得到生物組織光聲圖像[4]。

光聲圖像本質(zhì)是一種光吸收累積圖像，是組織吸收系數(shù)和光通量共同作用的結(jié)果。在實(shí)際應(yīng)用中，獲取組織的光吸收系數(shù)分布（稱為量化光聲成像）具有重要的價(jià)值。首先，吸收系數(shù)分布與生物組織中吸收發(fā)色團(tuán)的空間分布直接相關(guān)，并且不依賴于光通量，這提供了更真實(shí)的圖像對(duì)比度分布，便于深層信號(hào)的高對(duì)比成像;其次，直接從多波長(zhǎng)光聲圖像獲得的光譜成像是有問(wèn)題的，因?yàn)樵谌魏挝恢玫牟ㄩL(zhǎng)依賴性實(shí)際上是用整個(gè)被照射的組織體積中的所有發(fā)色團(tuán)的光譜特征進(jìn)行編碼的[5]。這種空間光譜串?dāng)_的效果會(huì)嚴(yán)重破壞感興趣發(fā)色團(tuán)的光譜特征，損害識(shí)別其存在并量化其濃度的能力。然而，若獲取組織在多個(gè)波長(zhǎng)下的吸收系數(shù)分布圖像可避免上述情況?？梢?jiàn)，量化光聲成像對(duì)組織光譜成像具有重要的意義。結(jié)合光聲光譜技術(shù)，可以定量分析組織中成分的變化，精細(xì)地反映極其微小的組織病變及血紅蛋白濃度、血氧濃度、氧代謝率等一系列重要的生理參數(shù)[6-7]，這在癌癥和心血管疾病診療方面具有廣闊的應(yīng)用前景。針對(duì)量化光聲成像問(wèn)題，國(guó)內(nèi)外很多研究小組開(kāi)展了大量研究。早期，Rinat O. Esenaliev曾經(jīng)提出了Temporal Profile模型[8]，主要用于在光聲成像中獲取生物組織的吸收系數(shù)，并用于監(jiān)測(cè)大腦中血氧濃度的變化。隨后，Lihong. V. Wang研究小組提出了體外相似組織測(cè)量方法[9]，在體外測(cè)量經(jīng)過(guò)小鼠腦皮和頭后骨的光的分布情況，用來(lái)矯正光吸收分布，實(shí)現(xiàn)了在小鼠大腦中對(duì)氧合血紅蛋白濃度的非侵入性成像。2016年，北京大學(xué)李長(zhǎng)輝教授基于比爾定律矯正光通量，實(shí)現(xiàn)了活體人手關(guān)節(jié)量化成像[10]。以上方法雖然無(wú)需求解光傳輸方程便可得到高精度的活體量化光聲成像，但是只適用于淺表組織或者特定深度的斷層圖像的量化，在臨床應(yīng)用中具有一定的局限性。從2006年開(kāi)始，許多小組開(kāi)始開(kāi)展基于光傳輸模型的深層量化光聲成像方法研究。該類研究方法可分為三類：1）基于混合光聲和擴(kuò)散光學(xué)成像系統(tǒng)（Diffuse Optical Tomography， DOT）的量化光聲成像方式。2011年Adam Q. Baure與Lihong. V. Wang教授小組合作，應(yīng)用時(shí)域DOT系統(tǒng)和基于超聲陣列的光聲層析成像系統(tǒng)（Photoacoustic Tomography， PAT），實(shí)現(xiàn)了仿體組織的深層量化重建[11]。2012年，Huabei Jiang研究小組設(shè)計(jì)了在同一個(gè)平臺(tái)上的PAT/DOT雙模成像系統(tǒng)，基于這個(gè)系統(tǒng)，該小組已經(jīng)開(kāi)展了仿體和體外乳腺癌組織成像的研究[12]。這種集成的DOT/PAT的雙模成像量化方式，不需要知道組織背景先驗(yàn)信息，可以同時(shí)重建組織吸收和散射系數(shù)，但是需要兩套成像系統(tǒng)的數(shù)據(jù)采集和信號(hào)重建，系統(tǒng)的復(fù)雜度高，計(jì)算量大。2）單步重建方式。單光聲系統(tǒng)下，通過(guò)聯(lián)合求解光聲方程和光傳輸方程實(shí)現(xiàn)組織吸收系數(shù)重建。Huabei Jiang研究小組在假定背景散射系數(shù)和Gruneisien參數(shù)已知的情況下，應(yīng)用有限元方法和傳輸方程，量化求解人手關(guān)節(jié)組織的吸收系數(shù)[13]。2015年，Sarah Vallelian小組提出了一種用于定量重建光聲圖像的一步算法，該算法可用于同時(shí)從多個(gè)聲學(xué)數(shù)據(jù)集中恢復(fù)吸收系數(shù)和介質(zhì)的超聲速度，數(shù)值模擬實(shí)驗(yàn)說(shuō)明了所提出重建方法的可行性[14]。同年，Simon Rabanser小組提出了使用一步重建算法進(jìn)行定量光聲層析成像[15]，通過(guò)將輻射傳輸方程與聲波方程耦合，將定量光聲層析成像的圖像重建問(wèn)題制定為單個(gè)非線性逆問(wèn)題，數(shù)值仿真實(shí)驗(yàn)驗(yàn)證了算法的可行性。單步量化方法計(jì)算量大，實(shí)際應(yīng)用復(fù)雜性高，目前大多只用于數(shù)值仿真和仿體實(shí)驗(yàn)。3）兩步重建方式。先通過(guò)光聲成像系統(tǒng)求解出光吸收累積圖像;然后，根據(jù)光聲系統(tǒng)的光源分布，基于光傳輸模型求解光通量;最后，通過(guò)一定的方法在原始光聲圖像中去除光通量，得到組織光吸收系數(shù)分布。2008年，Banerjee等[16]描述了一種非迭代方法，從邊界壓力測(cè)量值重建的吸收能量圖中恢復(fù)光學(xué)吸收系數(shù)分布。該小組通過(guò)蒙特卡羅模擬獲得邊界聲壓測(cè)量值。2017年，劍橋大學(xué)的Brochu等[8]基于輻射傳輸方程，實(shí)施了校正光通量的迭代優(yōu)化，并用于活體小鼠的吸收系數(shù)量化成像。倫敦大學(xué)的Cox和Bread小組[17]提出了一種簡(jiǎn)單的固定點(diǎn)迭代模型，在背景組織散射系數(shù)均勻且已知的情況下，該模型能夠重構(gòu)組織的光學(xué)吸收系數(shù)。該方法的局限是需要知道組織背景的光學(xué)信息，優(yōu)點(diǎn)是簡(jiǎn)單、快速，在軟組織量化光聲方面具有良好的應(yīng)用前景。2009年， Cox小組[18]基于固定點(diǎn)迭代模型開(kāi)展多波長(zhǎng)光譜成像方法，通過(guò)數(shù)值仿真實(shí)驗(yàn)，量化求解具有不同吸收系數(shù)的載色體濃度。2017年，北京大學(xué)李長(zhǎng)輝教授小組基于手持式光聲成像系統(tǒng)，結(jié)合固定點(diǎn)迭代方法，實(shí)現(xiàn)了人體乳腺組織的深層活體量化光聲成像[19]。

可見(jiàn)，固定點(diǎn)迭代方法受到了研究者的關(guān)注。然而，基于該方法的深層活體量化光聲成像研究剛剛起步，還有很多問(wèn)題需要解決。其中，限制視圖數(shù)據(jù)采集模式下重建偽影對(duì)深層量化成像結(jié)果有較大影響。因此，針對(duì)限制視圖光聲圖像的重建偽影問(wèn)題，本文提出了一種改進(jìn)的固定點(diǎn)迭代量化光聲成像方法。該方法是在迭代工作中對(duì)重建光聲圖像進(jìn)行偽影濾波，然后進(jìn)行深層組織量化成像。此外，本文在光傳輸模型前項(xiàng)求解方面融入了Toast++軟件，進(jìn)一步提高了量化光聲成像的效率和準(zhǔn)確度。實(shí)驗(yàn)中將量化重建算法用于一個(gè)新的128陣元商用超聲陣列的仿體和活體成像，實(shí)驗(yàn)結(jié)果驗(yàn)證了本文方法的有效性和準(zhǔn)確性。

1 相關(guān)研究

1.1 光聲成像

根據(jù)光聲生成理論研究，當(dāng)激光脈沖的時(shí)間寬度足夠短時(shí)，光聲波在無(wú)損介質(zhì)中產(chǎn)生和傳播滿足以下波動(dòng)方程[20]：

其中：p0（r）表示由脈沖激光激發(fā)的初始聲壓;c是聲速;r是空間位置坐標(biāo);t表示時(shí)間;表示光聲壓在位置r上的梯度;δ表示一個(gè)與時(shí)間成正相關(guān)的熱激發(fā)函數(shù)。假設(shè)成像區(qū)域的組織在聲學(xué)上是均勻的，可以通過(guò)求解以上波動(dòng)方程來(lái)獲取生物組織的光聲壓，因此，p（r，t）可以被表示為：

1.2 量化光聲成像

根據(jù)光聲成像的原理，由于生物組織被短脈沖激光照射，會(huì)瞬時(shí)使得局部組織的溫度升高，在這種情況下，如果局部區(qū)域的生物組織通過(guò)吸收激光脈沖加熱所產(chǎn)生的能量，就會(huì)產(chǎn)生聲波。如果不考慮其他因素的影響，假設(shè)聲波在生物組織中先行傳播，則聲壓p的波動(dòng)方程[17]是：其中：β是體積熱膨脹系數(shù);Cp是恒壓比熱容;Vs代表聲速;H是單位時(shí)間內(nèi)在每體積組織中沉積的熱量;p和H通常取決于在組織中的位置r和時(shí)間t。

當(dāng)激光脈沖足夠短時(shí)，熱函數(shù)可以被建模為以下與時(shí)間成正相關(guān)的δ函數(shù)[17]：

其中：A（r）是在位置r處轉(zhuǎn)化為具體熱量的吸收能量密度（即在r處具體的光吸收）。如果在生物組織的局部區(qū)域中，在位置r處的吸收系數(shù)分布是μa（r），并且在該位置的光通量是φ（r， μa），則在該位置的吸收能量密度A（r）[19]可以表示為：

根據(jù)光聲效應(yīng)機(jī)理，在短脈沖激光照射條件下，可以將式（3）進(jìn)行化簡(jiǎn)，得到在組織激光脈沖照射后產(chǎn)生的聲壓，即初始的光聲壓力分布p0，其計(jì)算形式如式（6）所示：

其中：Γ是Grüeneisen系數(shù)，代表熱量轉(zhuǎn)換為壓力的效率，可以取0～1的任意常數(shù)[21]。在本文中不考慮Γ的影響，假定所有位置的Γ相同，且取值為1（即認(rèn)為所有的熱量全部轉(zhuǎn)換為探測(cè)器探測(cè)到的聲壓）。

根據(jù)式（5）和（6），可以得出在位置r處的初始?jí)毫0（r）與吸收系數(shù)μa（r）和光通量φ（r）的關(guān)系式：通常，由于強(qiáng)烈的光散射和吸收，光通量φ（r）在組織中分布不均勻，重建的初始?jí)毫0（r）不能準(zhǔn)確地表示吸收系數(shù)μa（r）的分布。特別地，當(dāng)隨著組織深度的增大，φ（r）大幅降低時(shí)，這種不均勻的光通量導(dǎo)致深層組織中光聲信號(hào)的對(duì)比度減小，影響深層信號(hào)成像。

1.3 光通量求解

量化光聲成像中，最主要的困難是求解未知的光通量。雖然光在組織中的傳播滿足麥克斯韋方程，但是由于涉及空間尺度，在渾濁（高散射）介質(zhì)中傳播的光很快就達(dá)到了實(shí)際計(jì)算的極限。擴(kuò)散方程是光傳輸模型的一種近似形式，它在高散射介質(zhì)中可以高精度的表示光子在組織中的傳輸，擴(kuò)散方程的數(shù)學(xué)公式[22]如下：

其中：k（r）=[3（μa（r）+μs（r））]-1是生物組織在某一位置的擴(kuò)散系數(shù);φ（r，ω）是求解得到的光通量分布; μa（r）是生物組織在某一位置的光學(xué)吸收系數(shù)， μs（r）是生物組織在某一位置的散射系數(shù);c是光在介質(zhì)中的傳播速;ω是角頻率;ζ是在邊界處折射率不匹配項(xiàng);v是外部邊界標(biāo)準(zhǔn);m是在邊界處的某一位置;q是在Ω范圍內(nèi)邊界上的光源分布; i是一個(gè)虛數(shù)單位。

擴(kuò)散方程的求解可以通過(guò)有限元方法實(shí)現(xiàn)[23]。本文引入Toast++軟件，實(shí)現(xiàn)了擴(kuò)散方程的前向求解。首先，根據(jù)求解問(wèn)題得到該問(wèn)題的變分形式;其次，選擇適合成像區(qū)域的子空間區(qū)域，并根據(jù)子空間區(qū)域?qū)⒊上駞^(qū)域進(jìn)行離散化，在本文中，采用了三角剖分形式;再次，基于Toast++軟件中的dotSysmat（）函數(shù)，計(jì)算單元?jiǎng)偠炔⒔M合成總剛度矩陣K;最后，根據(jù)公式Kφ=q的變形公式φ=K-1q來(lái)求解光通量。

Toast++軟件是一個(gè)開(kāi)源軟件，是由倫敦大學(xué)學(xué)院的Martin Schweiger研究小組[22]開(kāi)發(fā)，主要用于解決DOT成像中的正向和逆向問(wèn)題。該軟件包括一組函數(shù)庫(kù)，用于模擬具有復(fù)雜邊界和不均勻內(nèi)部參數(shù)分布的高散射介質(zhì)中的近紅外光傳播，可以對(duì)穩(wěn)態(tài)、時(shí)域和頻域數(shù)據(jù)采集系統(tǒng)進(jìn)行建模。在使用過(guò)程中，可以直接鏈接到Toast ++庫(kù)以編譯DOT的應(yīng)用程序，并且可以快速地調(diào)用Toast++函數(shù)庫(kù)中的函數(shù)。

2 重建方法與實(shí)驗(yàn)設(shè)置

2.1 固定點(diǎn)迭代方法

傳統(tǒng)的固定點(diǎn)迭代重建算法是一個(gè)簡(jiǎn)單的迭代機(jī)制，適用于在重建過(guò)程中從生物組織的吸收能量中去除光通量對(duì)恢復(fù)生物組織的光學(xué)系數(shù)（吸收系數(shù)）的影響，從而根據(jù)假設(shè)的生物組織的光學(xué)系數(shù)來(lái)重建出生物組織真實(shí)的光學(xué)吸收系數(shù)的分布，得到生物組織真實(shí)的光學(xué)圖像[17]。

傳統(tǒng)的固定點(diǎn)迭代算法主要分為以下幾個(gè)步驟：

1）首先設(shè)置初始的光學(xué)吸收系數(shù)μa=0;

2）根據(jù)DOT模型求解光通量分布φ;

3）根據(jù)光通量和光學(xué)吸收系數(shù)計(jì)算光聲壓強(qiáng);

4）計(jì)算原始光聲壓強(qiáng)與重建光聲壓強(qiáng)之間的誤差;

5）迭代計(jì)算μa和φ，直到光聲壓強(qiáng)誤差小于設(shè)定閾值。

2.2 改進(jìn)的固定點(diǎn)迭代方法

傳統(tǒng)的固定點(diǎn)迭代方法重建光學(xué)吸收系數(shù)分布，計(jì)算過(guò)程簡(jiǎn)單并且能夠得到較為精確的重建結(jié)果;但在限制視圖光聲成像系統(tǒng)中，深層組織的重建光聲圖像存在較多的重建偽影，嚴(yán)重影響量化光聲圖像的精度和質(zhì)量。因此，針對(duì)這個(gè)問(wèn)題，本文提出了一種改進(jìn)的固定點(diǎn)迭代量化光聲成像方法，該方法融合了Toast++軟件，具體流程如圖1所示。

整個(gè)流程主要包括三個(gè)階段：

1）重建原始光聲圖像。在這個(gè)過(guò)程中，首先，利用傳統(tǒng)的反投影重建算法重建經(jīng)探測(cè)器探測(cè)到的原始光聲壓數(shù)據(jù)得到原始的光聲壓圖像。然后，采用2-D自適應(yīng)維納濾波去除原始光聲壓圖像中的噪聲和偽影。自適應(yīng)維納濾波是一種使噪聲圖像和其濾波圖像之間的均方誤差最小的復(fù)原方法，并根據(jù)局部方差的大小調(diào)整濾波效果[24]。在該方法中，使用大小為m×n的鄰域（模板）估計(jì)局部圖像均值和標(biāo)準(zhǔn)偏差。如果鄰域選擇太小，則只能對(duì)噪聲濾除產(chǎn)生很小的作用;如果鄰域選擇太大，則會(huì)導(dǎo)致濾波后的圖像失真。因此，在本文中通過(guò)選擇幾組不同鄰域大小的模板進(jìn)行實(shí)驗(yàn)，最終選取12×12的模板作為維納濾波在光聲圖像去噪應(yīng)用中的鄰域大小。

2）光通量求解。在這個(gè)過(guò)程中，在組織光學(xué)散射系數(shù)已知的情況下，首先利用Toast++軟件中的toastMesh（）函數(shù)結(jié)合有限元算法對(duì)成像區(qū)域進(jìn)行網(wǎng)格剖分;其次，設(shè)置初始的光學(xué)吸收系數(shù)，在本文中，將目標(biāo)組織的光學(xué)系數(shù)假設(shè)為背景組織的光學(xué)吸收系數(shù);再次，根據(jù)toast++軟件中的dotSysmat（）函數(shù)，輸入剖分好的成像區(qū)域網(wǎng)格，設(shè)置的初始光學(xué)參數(shù)，系統(tǒng)頻率等參數(shù)，計(jì)算每個(gè)單元?jiǎng)偠炔⒔M合成總剛度矩陣K;最后，根據(jù)公式φ=K-1q得出第k次迭代中的光通量φ（k）。

3）更新光學(xué)參數(shù)分布。首先，根據(jù)1.2節(jié)中的式（7）求解第k次迭代中的光聲壓強(qiáng)p（k）c;其次，迭代更新，得到光學(xué)吸收系數(shù)μ（k+1）a，直到濾波后的光聲壓圖像與計(jì)算光聲壓圖像之間的誤差小于設(shè)定的閾值ε。

單陣元大小 300μm×900μm;3）128-通道數(shù)據(jù)采集卡。在實(shí)驗(yàn)中采用單側(cè)打光的照明方式，用一根1500μm的多模光纖打光，光經(jīng)過(guò)兩個(gè)透鏡在成像組織表面形成一個(gè)大小為38mm×5mm的矩形光斑[25]。為實(shí)現(xiàn)超聲信號(hào)的耦合探測(cè)，系統(tǒng)準(zhǔn)備了裝滿去離子水的水槽，水槽下部有一個(gè)由低密度聚乙烯薄膜封閉的窗口，用于激光照射和信號(hào)采集?；谠撓到y(tǒng)，每一幀二維圖像的數(shù)據(jù)采集需1次激光脈沖照射，三維數(shù)據(jù)可以通過(guò)機(jī)械掃描獲得。

2.3.2 仿體實(shí)驗(yàn)與活體實(shí)驗(yàn)

仿體模型為在瓊脂塊中插入針灸針。實(shí)驗(yàn)中，瓊脂塊中Agar的濃度為1.5g/100ml。在制作過(guò)程中，在長(zhǎng)方體容器中加入了800ml水，并加入12g Agar，加熱到100℃，然后進(jìn)行自然冷卻凝固，凝固后瓊脂塊的高度為72mm、寬度為80mm，將四根直徑大約為0.3mm的針灸針在不同深度斜著嵌入到長(zhǎng)方體瓊脂塊模型中，如圖2所示（x是換能器陣列的（橫向）方向，y是機(jī)械掃描方向，z是深度方向）。圖2中四個(gè)實(shí)心黑點(diǎn)表示要成像的四根針灸針。采用波長(zhǎng)為780nm的激光照射扦插鋼針的仿體表面，配合機(jī)械掃描采集光聲信號(hào)。

對(duì)于活體成像，將激光輸出調(diào)節(jié)至1064nm（一個(gè)等吸光點(diǎn)），在該點(diǎn)處氧合血紅蛋白和脫氧血紅蛋白均勻地吸收光。對(duì)于小動(dòng)物成像，使用Sprague Dawley大鼠（雌性）。

在麻醉的過(guò)程中，實(shí)施腹腔內(nèi)注射氯胺酮（85mg/kg）和甲苯噻嗪（15mg/kg）混合物來(lái)深度麻醉大鼠，等待一段時(shí)間確認(rèn)大鼠已處于麻醉狀態(tài)（呼吸、心跳平穩(wěn)，按壓腳趾無(wú)疼痛反射）。在實(shí)驗(yàn)過(guò)程中，首先利用溫和脫毛膏除去大鼠腹部皮膚表面上的毛發(fā)（為了防止表皮毛發(fā)對(duì)實(shí)驗(yàn)結(jié)果的干擾），然后在大鼠腹部表皮上面涂抹一層薄薄的超聲耦合膠（避免探測(cè)聲壓直接在空氣中傳播，提高耦合效率），用動(dòng)物支架將大鼠固定到成像系統(tǒng)中，并給大鼠通上氧氣，在大鼠身下配備加熱墊，使其體溫始終保持在37±0.5℃。將大鼠腹部皮膚表面的光通量設(shè)定為每脈沖約～19mJ/cm2。實(shí)驗(yàn)持續(xù)1h，在整個(gè)實(shí)驗(yàn)過(guò)程中，本研究確保應(yīng)用此大小的脈沖能量對(duì)組織沒(méi)有熱損傷。實(shí)驗(yàn)結(jié)束后，動(dòng)物自然恢復(fù)體能，繼續(xù)放回籠中飼養(yǎng)。本文描述的動(dòng)物均在中國(guó)科學(xué)院深圳先進(jìn)技術(shù)研究院和曲阜師范大學(xué)動(dòng)物實(shí)驗(yàn)倫理動(dòng)物委員會(huì)批準(zhǔn)的順應(yīng)性方案中進(jìn)行。

3 實(shí)驗(yàn)結(jié)果與分析

3.1 仿體實(shí)驗(yàn)結(jié)果

圖3展示了濾波后的仿體光聲圖像以及重建后的光學(xué)吸收系數(shù)分布圖像，在本實(shí)驗(yàn)中只展示了一幀B-scan仿體實(shí)驗(yàn)數(shù)據(jù)。由于780nm波長(zhǎng)下，水是主要的光學(xué)吸收物質(zhì)，瓊脂本身既沒(méi)有散射也沒(méi)有吸收。因此，在本文實(shí)驗(yàn)中，假設(shè)瓊脂塊的吸收系數(shù)與散射系數(shù)在此波長(zhǎng)下與水的吸收和散射系數(shù)相同，分別設(shè)為μa=0.02346cm-1， μs=1.0cm-1。通過(guò)將濾波后的光聲圖像與光通量圖分開(kāi)來(lái)實(shí)現(xiàn)光通量補(bǔ)償從而得到目標(biāo)物體的光學(xué)吸收系數(shù)分布圖像。圖3（a）是經(jīng)過(guò)濾波后的光聲壓圖像，從中可以看出：由于光通量分布不均勻，導(dǎo)致在相同光學(xué)系數(shù)下獲得的位于不同深度的四根針灸針的光學(xué)對(duì)比度不同;處于較深位置的兩根針灸針的對(duì)比度較差（圖中白色圓圈所示）。通過(guò)傳統(tǒng)的固定點(diǎn)迭代算法與本文算法，量化重建得到了具有光通量補(bǔ)償?shù)墓鈱W(xué)吸收系數(shù)圖像如圖3（b）、（c）所示。圖3（b）、（c）是利用傳統(tǒng)的固定點(diǎn)迭代算法和本文算法量化得到的光學(xué)吸收系數(shù)分布圖像。從圖3（c）標(biāo)記的兩個(gè)位置可以看出，在較深位置的兩個(gè)針灸針與較淺位置的對(duì)比度基本相同，也就是說(shuō)，經(jīng)過(guò)本文算法也可以定量地恢復(fù)深層仿體目標(biāo)的光學(xué)吸收系數(shù)，獲得目標(biāo)物體相對(duì)真實(shí)的光學(xué)吸收分布圖像，實(shí)現(xiàn)對(duì)深層組織的量化;并且在重建光學(xué)吸收系數(shù)分布的過(guò)程中，通過(guò)對(duì)原始的光聲壓信號(hào)進(jìn)行濾波，重建得到的光學(xué)吸收系數(shù)分布中存在較少的噪聲和偽影（圖3（b）和（c）中白色方框所示）。

為了定量評(píng)估本文算法的效果，用圖3（d）、（e）分別展示了位于四個(gè)不同位置的濾波后的光聲壓數(shù)值分布與重建得到的吸收系數(shù)數(shù)值分布。這四個(gè)位置分別在圖3（a）中用數(shù)字標(biāo)記。從圖3（d）中可以清楚地看出，在位置2、3、4這三個(gè)位置點(diǎn)上，隨著組織深度加深，光聲壓的數(shù)值逐漸減小，并且數(shù)值之間差別較大。但是經(jīng)過(guò)本文算法所重建得到的吸收系數(shù)之間的數(shù)值差別較?。▓D3（e）所示）?？梢?jiàn)，本文算法能有效補(bǔ)償光通量，能夠定量恢復(fù)組織的光學(xué)吸收系數(shù)分布。理論上，相同的物質(zhì)具有相同的光學(xué)吸收系數(shù)，針灸針屬于鋼材料，鋼是鐵碳合金，主要由鐵元素和碳元素構(gòu)成，且鐵元素含量較高，因此在本文中把在780nm波長(zhǎng)下鐵的吸收系數(shù)0.6cm-1當(dāng)成量化標(biāo)準(zhǔn)。通過(guò)圖3（e）可以看出雖然本文算法能夠在這四個(gè)位置上得到較為接近的光學(xué)吸收系數(shù)，但是每個(gè)位置上的吸收系數(shù)值仍然存在誤差。這是因?yàn)樵谕垦a(bǔ)償?shù)倪^(guò)程中，目前還沒(méi)有一種有效的算法能夠進(jìn)行精確的通量補(bǔ)償，本文研究的算法在進(jìn)行通量補(bǔ)償時(shí)也存在誤差，因此本文算法也只是一種相對(duì)量化，重建得到的是目標(biāo)物體光學(xué)吸收系數(shù)的相對(duì)值。

3.2 活體實(shí)驗(yàn)結(jié)果

由于皮膚下的主要光聲信號(hào)來(lái)源是血管，所以在本文研究中假設(shè)除了血管以外的其他背景組織的光學(xué)系數(shù)都為軟組織的光學(xué)系數(shù)。由于1064nm波長(zhǎng)下，軟組織中水是主要的吸收物質(zhì)，因此，組織背景的光學(xué)吸收系數(shù)設(shè)置為1064nm波長(zhǎng)下水的光學(xué)吸收系數(shù)，即μa=0.12cm-1[19]。該波長(zhǎng)下軟組織的散射系數(shù)為μs=7.1563cm-1[26]。

利用傳統(tǒng)的固定點(diǎn)迭代算法與本文算法量化重建的實(shí)驗(yàn)結(jié)果比較如圖4所示。圖4中展示了四幀不同的B-scan活體實(shí)驗(yàn)數(shù)據(jù)。其中圖4（a）是濾波后的光聲壓圖像。在實(shí)驗(yàn)結(jié)果展示圖中，本文標(biāo)記了兩個(gè)不同深度的血管（圖4中白色虛線圓圈所示）。如前所述，隨著深度的增加，由于光通量的影響，重建的光聲信號(hào)對(duì)比度變得越來(lái)越弱，當(dāng)達(dá)到一定深度，血管基本上被淹沒(méi)在背景中。盡管通過(guò)傳統(tǒng)的固定點(diǎn)迭代算法能夠恢復(fù)目標(biāo)組織的光學(xué)吸收系數(shù)（如圖4（b）所示），但是量化重建的結(jié)果中存在較多偽影。通過(guò)圖4（c）可以看出，本文算法能夠有效抑制光通量的影響，得到生物組織相對(duì)真實(shí)的光學(xué)吸收分布圖像。進(jìn)一步說(shuō)明，本文方法可以實(shí)現(xiàn)對(duì)深層組織的光學(xué)系數(shù)的量化。通過(guò)比較圖4（b）和（c）白色虛線方框中的內(nèi)容，可以看出，當(dāng)使用改進(jìn)的固定點(diǎn)迭代方法時(shí)，背景中的大部分偽影被移除，提高了圖像的保真度，獲得較高質(zhì)量的光學(xué)吸收系數(shù)圖像，實(shí)現(xiàn)更加精確的量化重建。

同樣為了評(píng)估本文方法的優(yōu)勢(shì)，在圖5中分別展示了經(jīng)過(guò)濾波后的這四幀不同B-Scan的光聲壓數(shù)值分布與重建得到的吸收系數(shù)值分布。其中，主要展示了位于四個(gè)不同深度的信號(hào)點(diǎn)的數(shù)值分布情況，這四個(gè)位置分別在圖4（a）中用數(shù)字標(biāo)記。圖5（a）是經(jīng)過(guò)濾波后的光聲壓數(shù)值分布，圖5（b）是利用本文算法量化重建后獲得的光學(xué)吸收系數(shù)數(shù)值分布。從數(shù)值分布情況可以看出，隨著組織深度加深，位于不同位置的光聲壓的數(shù)值之間差別較大（如圖5（a）所示），尤其是當(dāng)目標(biāo)血管處于較深位置時(shí)，光聲壓值幾乎變得非常小。但是通過(guò)改進(jìn)的固定算法重建得到的目標(biāo)血管的光學(xué)吸收系數(shù)之間的數(shù)值差別較?。ㄈ鐖D5（b）所示）。當(dāng)目標(biāo)血管的位置較深時(shí)，重建獲得的光學(xué)吸收系數(shù)基本上與較淺位置血管的吸收系數(shù)相同?？梢?jiàn)，通過(guò)活體實(shí)驗(yàn)，說(shuō)明本文算法能有效補(bǔ)償光通量并能夠有效地恢復(fù)深層生物組織的光學(xué)吸收系數(shù)分布，在量化深層組織光學(xué)吸收系數(shù)具有一定的優(yōu)勢(shì)。

4 收斂性分析

在本文實(shí)驗(yàn)中還計(jì)算了當(dāng)前重建的光聲壓值與濾波后的光聲壓值之間的誤差。將從當(dāng)前重建得到的光聲壓與濾波后的光聲壓圖之間的差異與迭代次數(shù)用作迭代停止標(biāo)準(zhǔn)，如2.2節(jié)中圖1所示。在本文實(shí)驗(yàn)中，當(dāng)重建得到的光聲壓與濾波后的光聲壓圖之間的差異變得非常?。葱∮诒疚膶?shí)驗(yàn)所設(shè)計(jì)的閾值），并且當(dāng)?shù)螖?shù)達(dá)到最大迭代次數(shù)時(shí)，迭代過(guò)程將停止。

為了驗(yàn)證本文方法的收斂性，針對(duì)迭代過(guò)程繪制了兩個(gè)不同B-scan光聲圖像的濾波后的光聲壓與重建光聲壓隨迭代次數(shù)增加的誤差曲線。如圖6所示：圖（a）是仿體目標(biāo)濾波后的光聲壓與重建光聲壓隨迭代次數(shù)變化的誤差收斂曲線;

圖（b）是大鼠腹部血管誤差收斂曲線。誤差值在前兩次迭代中迅速減小，大約在第5次迭代后收斂。

原始聲壓與重建光聲壓之間的差異接近于零，表明本文方法能夠快速收斂。

5 結(jié)語(yǔ)

針對(duì)限制視圖下光聲圖像的重建偽影問(wèn)題，在本文中主要提出了一種改進(jìn)的固定點(diǎn)迭代量化光聲成像方法。對(duì)于復(fù)雜的深層生物組織，本文方法能夠有效抑制量化重建組織光學(xué)吸收系數(shù)圖像中的噪聲和偽影，獲取高質(zhì)量的組織光學(xué)吸收系數(shù)圖像。相對(duì)于傳統(tǒng)的固定點(diǎn)迭代方法，本文方法進(jìn)一步提高了量化光聲成像的效率和準(zhǔn)確度。

但是，本文的研究工作只對(duì)量化光聲成像方法進(jìn)行了初步的研究，仍然還存在未解決的問(wèn)題。在今后的工作中，我們將對(duì)3-D光聲圖像進(jìn)行量化重建，并且也會(huì)在本文算法中引入更高級(jí)的重建偽影去除方法，進(jìn)一步提高深層光聲圖像的量化重建質(zhì)量。

參考文獻(xiàn)（References）

[1] LIN L， HU P， SHI J， et al. Single-breath-hold photoacoustic computed tomography of the breast [J]. Nature Communications， 2018， 9（1）：2352.

[2] FERNANDES D A， KOLIOS M C. Intrinsically absorbing photoacoustic and ultrasound contrast Agents for cancer therapy and imaging [J]. Nanotechnology， 2018， 29（50）： 1361-6528.

[3] MATSUMOTO Y， ASAO Y， SEKIGUCHI H， et al. Visualising peripheral arterioles and venules through high-resolution and large-area photoacoustic imaging， Scientific Reports 8[R]. Nature， 2018： Article number： 14930.

[4] BARDSLEY P， REN K， ZHANG R. Quantitative photoacoustic imaging of two-photon absorption [J]. Journal of Biomedical Optics， 2018， 23（1）：1-11.

[5] LEE C H， FOLZ J， ZHANG W， et al. Ion-selective nanosensor for photoacoustic and fluorescence imaging of potassium [J]. Analytical Chemistry， 2017，89（24）：7943-7949.

[6] DAHAL S， CULLUM B M. Characterization of multiphoton photoacoustic spectroscopy for subsurface brain tissue diagnosis and imaging [J]. Journal of Biomedical Optics， 2016， 21（4）：No.47001.

[7] LIM H T， MURUKESHAN V M. Hyperspectral photoacoustic spectroscopy of highly-absorbing samples for diagnostic ocular imaging applications [J]. International Journal of Optomechatronics， 2017， 11（1）：36-46.

[8] BROCHU F M， BRUNKER J， JOSEPH J， et al. Towards quantitative evaluation of tissue absorption coefficients using light fluence correction in optoacoustic tomography[J]. IEEE Transactions on Medical Imaging， 2017， 36（1）：322-331.

[21] ZHANG J， ZHANG J， LI W. Degradable hollow mesoporous silicon/carbon nanoparticles for photoacoustic imaging-guided highly effective chemo-thermal tumor therapy in vitro and in vivo [J]. Theranostics， 2017， 7（12）： 3007-3020.

[22] SCHWEIGER M， ARRIDGE S. The Toast++ software suite for forward and inverse modeling in optical tomography [J]. Journal of Biomedical Optics， 2014， 19（4）： 040801.

[23] SHAN T， QI J， JIANG M， et al. GPU-based acceleration and mesh optimization of finite-element-method-based quantitative photoacoustic tomography： a step towards clinical applications [J]. Applied Optics， 2017， 56（15）： 4426-4432.

[24] 易三莉，賀建峰.基于BEMD與自適應(yīng)維納濾波的圖像降噪[J].計(jì)算機(jī)工程與應(yīng)用，2013，49（10）：156-158.

（YI S L， HE J F. Image denoising method based on BEMD and adaptive Wiener filter [J]. Computer Engineering and Applications， 2013， 49（10）： 156-158.）

[25] LI M， LIU C， GONG X， et al. Linear array-based real-time photoacoustic imaging system with a compact coaxial excitation handheld probe for noninvasive sentinel lymph node mapping [J]. Biomedical Optics Express， 2018， 9（4）： 1408-1422.

[26] JACQUES S L. Optical properties of biological tissues： a review [J]. Physics in Medicine and Biology， 2013， 58（11）： R37-61.