孫 正,鄭 蘭

定量光聲層析成像的研究進(jìn)展

孫 正*,鄭 蘭

(華北電力大學(xué) 電子與通信工程系，河北 保定 071003)

光聲層析(Photoacoustic tomography,PAT)成像結(jié)合了超聲成像的高分辨率和光學(xué)成像的高對(duì)比度的優(yōu)勢(shì)，是一種新型的生物醫(yī)學(xué)成像模式。PAT成像算法包含兩個(gè)逆問題，即根據(jù)組織產(chǎn)生的光聲信號(hào)構(gòu)建初始聲壓分布圖(即圖像重建)以及在此基礎(chǔ)上估算成像區(qū)域的光學(xué)特性參數(shù)。后者是一個(gè)非線性的不適定問題，通常稱為定量光聲層析(Quantitative photo-acoustic tomography，qPAT)成像。本文在介紹光聲成像原理的基礎(chǔ)上，對(duì)主要的qPAT算法進(jìn)行綜述，討論各自的優(yōu)勢(shì)和不足，并對(duì)未來可能的發(fā)展方向進(jìn)行展望。

定量光聲成像； 圖像重建； 逆問題； 光吸收系數(shù)； 散射系數(shù)； Gruneisen系數(shù)

1 引 言

光聲層析(Photoacoustic tomography,PAT)成像是近年來發(fā)展起來的一種非電離式的新型生物醫(yī)學(xué)成像方法。它采用脈沖激光對(duì)生物組織進(jìn)行激發(fā)，獲取組織的光吸收系數(shù)分布情況，因此繼承了光學(xué)成像的一系列優(yōu)點(diǎn)，如對(duì)組織損傷小、對(duì)比度高、可進(jìn)行功能成像和分子成像等[1]。PAT探測(cè)的是組織產(chǎn)生的超聲波信號(hào)，因此它繼承了超聲成像對(duì)深層組織成像的高分辨率的特點(diǎn)[2]。與X射線成像、超聲成像和磁共振成像等傳統(tǒng)的醫(yī)學(xué)成像技術(shù)相比，PAT技術(shù)的優(yōu)勢(shì)主要體現(xiàn)在：采用非電離波段，成像過程中不改變生物組織的屬性；較容易界定組織產(chǎn)生的光聲信號(hào)和組織的生理狀態(tài)之間的關(guān)系；可與超聲或光學(xué)成像技術(shù)相結(jié)合，以獲得更多的診療信息；可根據(jù)實(shí)際應(yīng)用的需要調(diào)整成像深度和成像分辨率等[3]。

目前，對(duì)生物PAT成像系統(tǒng)的改進(jìn)主要體現(xiàn)在三個(gè)方面：第一，探測(cè)組織產(chǎn)生的光聲信號(hào)的方法，除了采用傳統(tǒng)的基于壓電效應(yīng)的超聲換能器或者PVDF 膜的水聽器之外，還可利用光學(xué)方法實(shí)現(xiàn)光聲信號(hào)的探測(cè)。例如：探測(cè)光聲效應(yīng)引起的組織折射率的變化[4]、探測(cè)壓力傳感器的變形[5]或者采用附著在光纖束末端的FP(Fabry-Perot)聚合物薄膜作為超聲傳感器[6]；第二，提高成像速度，例如：通過提高脈沖激光光源的重復(fù)頻率，加快光聲信號(hào)的產(chǎn)生過程，縮短信號(hào)采集的等待時(shí)間[7]，或者提高光聲信號(hào)的采集速率從而節(jié)省數(shù)據(jù)的采集時(shí)間[8]等；第三，成像分辨率從超聲分辨水平發(fā)展到了光學(xué)分辨水平[9]。

PAT成像的逆問題包括聲學(xué)和光學(xué)兩方面：聲學(xué)逆問題是指根據(jù)探測(cè)器采集到的光聲信號(hào)(本質(zhì)是超聲波)重建組織內(nèi)部的初始聲壓分布或空間光吸收能量密度，即一般意義上的光聲圖像重建[10-12]；光學(xué)逆問題是指運(yùn)用合適的光傳輸模型與優(yōu)化算法，根據(jù)探測(cè)到的光聲信號(hào)或者光吸收能量密度，估算出組織的光學(xué)特性參數(shù)(包括光吸收系數(shù)和散射系數(shù))的空間分布[13]和熱膨脹特性參數(shù)(Gruneisen系數(shù))等[14-16]，得到對(duì)組織光學(xué)特性的定量評(píng)價(jià)，即定量光聲層析成像(Quantitative photoacoustic tomography,qPAT)。由于聲波在組織中傳播的時(shí)間比激光照射組織以及組織吸收光的時(shí)間長(zhǎng)約3個(gè)數(shù)量級(jí)，因此可以將聲學(xué)逆問題和光學(xué)逆問題分離開來分別求解[16]。聲學(xué)逆問題的重建結(jié)果，即光吸收能量密度是由局部的光吸收系數(shù)和光子數(shù)分布共同決定的，并不能反映組織的光吸收系數(shù)分布，而光吸收系數(shù)與組織的化學(xué)成分密切相關(guān)，正常和病變組織的光學(xué)特性參數(shù)之間通常有較明顯的差異，因此qPAT可為疾病的早期診斷提供更加準(zhǔn)確可靠的信息。

目前對(duì)qPAT的研究已成為光聲成像領(lǐng)域的熱點(diǎn)之一，尤其是同時(shí)重建光吸收系數(shù)和散射系數(shù)的空間分布。本文在簡(jiǎn)介光聲成像和qPAT原理的基礎(chǔ)上，對(duì)目前已提出的主要qPAT算法進(jìn)行總結(jié)和歸納，簡(jiǎn)介其中典型算法的原理，并討論各自的優(yōu)勢(shì)和不足。

2 光聲成像的原理

生物PAT成像是一種以超聲波作為媒介、以組織的光吸收系數(shù)和散射系數(shù)作為成像參數(shù)的生物光子成像方法，其物理基礎(chǔ)是生物組織的光聲效應(yīng)[3]。成像系統(tǒng)包括3個(gè)組成部分：光聲信號(hào)的產(chǎn)生單元、接收單元和后處理單元。一束短脈沖(～10 ns)激光經(jīng)過光學(xué)元件后照射到生物組織上，在組織內(nèi)部形成與組織光學(xué)特性參數(shù)相關(guān)的能量沉積分布。由于激光脈寬很窄，組織吸收的光子能量不能在短時(shí)間內(nèi)釋放，導(dǎo)致組織溫度瞬間的不均勻提升。周期性熱流使周圍的介質(zhì)產(chǎn)生熱膨脹從而激發(fā)出寬頻帶的超聲波，并迅速向組織邊界傳播。基于這種超聲波信號(hào)的特殊產(chǎn)生機(jī)理，為了區(qū)別于其他的超聲信號(hào)，通常稱其為光聲信號(hào)[17]。超聲換能器接收到光聲信號(hào)之后，經(jīng)計(jì)算機(jī)處理即可重建出組織內(nèi)部光能量沉積的分布圖，揭示病變組織的內(nèi)部信息。

如圖1所示，從組織內(nèi)的吸收體吸收光子到探測(cè)器測(cè)得聲壓時(shí)間序列的整個(gè)過程稱為光聲成像的正問題，可分為光學(xué)正問題和聲學(xué)正問題兩部分：前者的結(jié)果是光吸收能量密度，后者的結(jié)果是探測(cè)器測(cè)得的聲壓信號(hào)。假設(shè)采用單一波長(zhǎng)λ的光源照射組織，組織的光吸收系數(shù)和散射系數(shù)分別為μa和μs，待測(cè)組織區(qū)域Ω內(nèi)一點(diǎn)r處的光吸收能量密度H(r,λ)為

H(r,λ)=μa(r,λ)Φ(r,λ;μa(r,λ),μs(r,λ)),

(1)

其中，H(r,λ)定義為單位體積單位時(shí)間內(nèi)的熱能轉(zhuǎn)換；r∈Ω，Ω?Rn(Rn為n維實(shí)數(shù)空間，n=2,3)為有界域，Φ是光能流率，則r處的初始聲壓為

(2)

圖1 光聲成像正問題框圖Fig.1 Photoacoustic imaging block diagram of the forward problem

3 定量光聲成像的原理

3.1 聲學(xué)逆問題

假設(shè)介質(zhì)的聲學(xué)特性均勻，在理想激光脈沖的均勻照射下，被照組織產(chǎn)生的三維光聲信號(hào)的幅值與脈沖激光的幅值成正比，光聲信號(hào)的特性由光能量的吸收分布決定。因此，可以根據(jù)探測(cè)器測(cè)量到的聲壓時(shí)間序列p(r,λ,t)重建初始聲壓的空間分布p0(r,λ)，進(jìn)而得到光吸收分布H(r,λ)，即PAT聲學(xué)逆問題，也就是通常所說的光聲圖像重建。

如果超聲波在密度均勻的無(wú)損介質(zhì)中的傳播速度是均勻的，并且光激發(fā)可被看作是瞬時(shí)的，那么超聲波傳播的物理模型可由以下3個(gè)方程表示[9]：

(3)

式中，p(r,t)為r∈Ω?Rn處、t∈R+時(shí)刻的聲壓；u(r,t)為介質(zhì)的振動(dòng)速度；c為聲速；ρ(r,t)為位置r處、時(shí)刻t時(shí)的聲學(xué)密度；ρ0是介質(zhì)密度。

從不同的應(yīng)用角度出發(fā)，借鑒其他成像技術(shù)，研究者已提出了多種圖像重建算法，例如濾波反投影法(Filtered back-projection,FBP)、時(shí)間反演法(Time-reversal,TR)、相控聚焦算法、基于傅里葉變換的重建算法、反卷積重建算法和迭代重建算法等[12]。

3.2 光學(xué)逆問題

PAT光學(xué)逆問題是指由初始聲壓分布p0(r,λ)估算組織的光吸收系數(shù)和散射系數(shù)的空間分布。在早期的研究中，通常假設(shè)組織的光散射系數(shù)是已知的，那么采用遞歸法[18]或者非遞歸的方法[19]即可重建出光吸收系數(shù)的分布。但是，該假設(shè)在多數(shù)情況下都是不成立的，更為通用的方法是基于誤差最小化的方法。由式(2)可知，若待測(cè)組織內(nèi)部的光吸收系數(shù)和散射系數(shù)是有界常數(shù)且滿足Lipschitz連續(xù)，且已知Gruneisen系數(shù)，光吸收能量密度的測(cè)量值為

(4)

那么，光學(xué)逆問題可表述為如下的非線性最小二乘問題[20-21]

C·H(r,λ;μa(r,λ),μs′(r,λ))‖2)， (5)

其中，μs′是組織的約化散射系數(shù)

μs′=μs(1-g)，

(6)

光能流率Φ是未知的，而且它與組織的光吸收系數(shù)和散射系數(shù)有關(guān)，同時(shí)PAT成像是三維高分辨率成像，所涉及的數(shù)據(jù)量極大，因此由光吸收能量密度的測(cè)量值重建組織的光學(xué)特性參數(shù)是一個(gè)大規(guī)模的非線性不適定問題，特別是當(dāng)同時(shí)重建光吸收系數(shù)和散射系數(shù)時(shí)[22]。

求解光學(xué)逆問題需要解決兩個(gè)關(guān)鍵問題[20]：第一，建立光在組織中傳輸?shù)那跋蚰Ｐ停Ｐ偷妮敵鼍褪枪馕漳芰棵芏鹊睦碚撝?；第二，選擇適當(dāng)?shù)膬?yōu)化策略，使光吸收能量密度的測(cè)量值與前向模型的輸出值之間的誤差最小，進(jìn)而得到光學(xué)參數(shù)的估計(jì)值。

3.3 光在組織中傳輸?shù)那跋蚰Ｐ?/p>

通常采用輻射傳輸方程(Radiative transfer equation，RTE)描述光子在混濁介質(zhì)中的遷移過程[22]。RTE是積分-微分方程，求解時(shí)常需要在空間域和角度域內(nèi)對(duì)方程進(jìn)行離散化，步驟較為繁瑣，因而通常對(duì)其進(jìn)行擴(kuò)散近似(Diffusion approximation，DA)。相比于DA，RTE能夠更準(zhǔn)確地描述光子在組織中的遷移過程，尤其是在非擴(kuò)散區(qū)域，但是其復(fù)雜的求解過程和較高的運(yùn)算成本限制了它的廣泛應(yīng)用。

3.3.1 輻射傳輸方程

RTE描述了在特定控制體(Control volume)內(nèi)的能量守恒：當(dāng)光在待測(cè)區(qū)域內(nèi)沿某一方向傳輸時(shí)，組織對(duì)光子的吸收、偏離光傳輸方向上光子的散射和超過待測(cè)區(qū)域的光子外流出都會(huì)導(dǎo)致能量的損失，而沿光傳輸方向的光子散射或介質(zhì)中的其他光源又會(huì)造成能量的增加。

(7)

(8)

假設(shè)除了光源Γs??Ω處以外，在待測(cè)區(qū)域邊界?Ω處光子不會(huì)向內(nèi)傳輸，那么方程(7)的邊界條件是

(9)

(10)

其中，Φ是光能流率，也稱為光輻射通量。

3.3.2 擴(kuò)散近似

假設(shè)組織的約化散射系數(shù)遠(yuǎn)大于吸收系數(shù)，即光在組織中發(fā)生散射的幾率遠(yuǎn)高于吸收(對(duì)于人體的大多數(shù)組織，該假設(shè)是成立的)，并且光在組織中的傳輸和分布近似于各向同性，則可將RTE簡(jiǎn)化為擴(kuò)散方程(Diffusion equation,DE)，它是用球諧函數(shù)表示的RTE的一階相位近似。

在DA模型中，光子密度被簡(jiǎn)化為只與空間變量r有關(guān)，設(shè)有界域Ω的光滑邊界為Γ，則滿足Robin邊界條件的時(shí)間獨(dú)立的DE可寫作[14,20]

-·(D(r)Φ(r))+μa(r)Φ(r)=S(r),r∈Ω，

(11)

其中，S(r)是源項(xiàng)，D(r)是擴(kuò)散系數(shù)：

(12)

為了求式(11)的數(shù)值解，需要確定適當(dāng)?shù)倪吔鐥l件。例如可采用Dirichlet邊界條件[20]，即在外推邊界處(距離介質(zhì)邊界2D處)，令Φ=0。計(jì)算出Φ(r)后，則

H(r)=μa(r)Φ(r),

(13)

DE的求解方法相對(duì)較簡(jiǎn)單，只需在空間域內(nèi)對(duì)方程離散化即可，因而它是目前qPAT中最常用的前向模型[14]。在較為理想的情況下，例如各向同性的渾濁介質(zhì)中嵌入柱狀的不均勻物質(zhì)，求解DE可以得到光吸收能量密度的解析解[24]。對(duì)于實(shí)際的生物組織，則需要采用有限差分法或者有限元法得到DE的數(shù)值解。

DA僅在擴(kuò)散域內(nèi)有效，即在距離光源幾個(gè)遷移平均自由程(Transportmeanfreepaths)的區(qū)域內(nèi)μa=μs′，在該區(qū)域內(nèi)光由于發(fā)生多次散射，因而失去方向性而擴(kuò)散開。但是對(duì)于PAT成像，接近光源的區(qū)域和待測(cè)區(qū)域的邊界構(gòu)成了圖像的重要組成部分，通常包含診療疾病所需的重要信息，在這些區(qū)域內(nèi)，光的傳輸是高度各向異性的，因此DA是不成立的[23]。

3.4 主要的優(yōu)化方法

為了得到光在組織中傳輸?shù)那跋蚰Ｐ偷臄?shù)值解，一般需要先對(duì)其進(jìn)行有限元離散化，相比一般的單網(wǎng)格方法，采用雙網(wǎng)格方法可在保證重建精度的前提下，明顯縮短重建時(shí)間，提高計(jì)算效率[20,25]。在此基礎(chǔ)上求解式(5)的非線性最小二乘問題，使光吸收能量密度的測(cè)量值與其理論值之間的均方誤差最小，即可求得光吸收系數(shù)和散射系數(shù)的空間分布。最常用的方法是Levenberg-Marquardt(LM)方法[26]，即L2范數(shù)Gauss-Newton法[27]，通過計(jì)算Hessian矩陣的逆矩陣，可迭代地調(diào)整(μa,μs)的值。但是該計(jì)算過程非常耗時(shí)，因此出現(xiàn)了近似計(jì)算逆Hessian矩陣的算法，主要包括基于Jacobian矩陣的線性方法[21,28]、非線性梯度法[23,14,29-32]、Bregman迭代法[15]等，那么光學(xué)逆問題解的精度將很大程度上取決于所選數(shù)值模型的精度。

3.4.1 基于Jacobian矩陣的線性方法

考慮到Gauss-Newton法存在運(yùn)算成本過高的問題，可采用吸收系數(shù)和散射系數(shù)的Jacobian矩陣構(gòu)建Hessian矩陣的近似矩陣，用于在Gauss-Newton迭代中更新(μa,μs)的值。此類方法的優(yōu)點(diǎn)是原理簡(jiǎn)單，易于實(shí)現(xiàn)，而且它保留了Gauss-Newton法的二階收斂性。但是，PAT是三維成像，數(shù)據(jù)量巨大，因此得到的Jacobian矩陣是大型密集矩陣，對(duì)計(jì)算速度和存儲(chǔ)空間的要求都很高，而且計(jì)算該密集矩陣的逆矩陣也非常耗時(shí)，可能降低此種方法的實(shí)用性[23]。文獻(xiàn)[27]提出一種無(wú)矩陣(Matrix-free)求解近似Hessian矩陣的方法，即只計(jì)算矩陣向量積，可降低對(duì)存儲(chǔ)空間的要求，可是計(jì)算量仍然很大。

3.4.2 非線性梯度法

此類方法采用基于梯度的擬牛頓(Quasi-Newton)最小化法[33-35]求解式(5)的最優(yōu)化問題，避免計(jì)算Jacobian矩陣和Hessian矩陣，而只需計(jì)算式(5)中誤差函數(shù)的梯度，從而迭代更新光吸收系數(shù)和散射系數(shù)的值。通常需要引入正則項(xiàng)，降低問題的病態(tài)性，使之成為良態(tài)問題，如L2正則化。采用多位置激光照射的方法采集多個(gè)光聲信號(hào)數(shù)據(jù)集，在已知Gruneisen系數(shù)的前提下，采用此類方法可同時(shí)、唯一地確定光吸收系數(shù)和散射系數(shù)。與基于Jacobian矩陣的方法相比，該類算法具有超線性收斂性，且計(jì)算速度快，所需的存儲(chǔ)空間小，對(duì)三維圖像的重建更具實(shí)用性。但是其計(jì)算較復(fù)雜，一般需要較多的迭代次數(shù)。

3.4.3 正則化方法

對(duì)于病態(tài)問題，采用正則化方法(即引入正則項(xiàng))可降低問題的病態(tài)程度，把非線性問題轉(zhuǎn)化為某種條件下的線性問題。因此在解決qPAT這一非線性病態(tài)逆問題的過程中，選取合適的正則化方法可使算法更加簡(jiǎn)單，保證穩(wěn)定近似解的存在[36]。對(duì)正則化方法的選取依賴于光吸收系數(shù)和散射系數(shù)的先驗(yàn)知識(shí)。

(1) L1和L2正則化法

(2) Tikhonov正則化方法

基于變分原理的Tikhonov正則化方法[39]是研究不適定問題的最重要的正則化方法之一，即在求解過程中結(jié)合解的某些已知信息對(duì)解進(jìn)行限制。PAT是對(duì)不同的組織區(qū)域進(jìn)行成像，因此對(duì)組織結(jié)構(gòu)光滑性的約束也不同，這一信息可用于設(shè)定Tikhonov方法中解的光滑性假設(shè)，反映光學(xué)參數(shù)的光滑性[36]。對(duì)于光滑(連續(xù))變化的光學(xué)參數(shù)，可采用標(biāo)準(zhǔn)Tikhonov正則化方法[39]，設(shè)定解的光滑性，反映光學(xué)參數(shù)的光滑性。PAT非常適合于對(duì)血管這種非光滑性的結(jié)構(gòu)進(jìn)行成像，這種情況下，將組織的光吸收系數(shù)和散射系數(shù)看作是分段連續(xù)的函數(shù)更為適合[40-42]，文獻(xiàn)[36]設(shè)計(jì)了一種基于標(biāo)準(zhǔn)Level-set[43]的正則化方法，求解分段連續(xù)的光吸收系數(shù)和散射系數(shù)。此外，還可以考慮采用迭代的正則化方法[44]。

(3) 基于全變差(Total-variation,TV)的非光滑正則化方法

對(duì)于待求光學(xué)特性參數(shù)是分段常數(shù)的情形，帶約束的全變差(或全變分)正則化方法比L2正則化法更為適合[15,21]，全變分項(xiàng)可以良好地約束信號(hào)的平滑性。如果將TV與Bregman方法相結(jié)合，則解的精度和計(jì)算效率都會(huì)得到大幅提升。

3.4.4 Bregman迭代法

Bregman迭代法采用基于次梯度[45]的最小化法求解式(5)的最優(yōu)化問題，基本模型是ROF(Rudin-Osher-Fatemi)模型。該方法可以解決非線性的非凸優(yōu)化問題，把非凸問題轉(zhuǎn)換成凸問題去逼近，通常需結(jié)合其他優(yōu)化方法。

目前，Bregman迭代法在qPAT中的應(yīng)用還處于初步研究階段?；舅悸肥牵涸O(shè)定光吸收系數(shù)和散射系數(shù)的初值為0，在每次迭代過程中通過更新次梯度極小化Bregman距離[45]，達(dá)到目標(biāo)函數(shù)最小化，最終得到光吸收系數(shù)和散射系數(shù)的最優(yōu)解。該算法的計(jì)算效率高、收斂速度快。引入正則化后，正則化參數(shù)是不變的，且在迭代過程中系統(tǒng)穩(wěn)定。

在qPAT中常采用以下兩種Bregman迭代法求解最優(yōu)化問題[15]：

(1)基于Jacobian矩陣的線性Bregman迭代法

采用吸收系數(shù)和散射系數(shù)的Jacobian矩陣構(gòu)建Hessian矩陣的近似矩陣，用于在Bregman迭代法中迭代更新吸收系數(shù)和散射系數(shù)的值。在每次迭代過程中求解的是關(guān)于Jacobian矩陣的誤差函數(shù)最小值，比直接計(jì)算Jacobian矩陣簡(jiǎn)單。也可將TV正則化與Bregman迭代法結(jié)合，在目標(biāo)函數(shù)中增加由當(dāng)前迭代過程產(chǎn)生的誤差(其初始值為0)，用于修正下一次迭代，可提高結(jié)果的精度，但會(huì)增加迭代次數(shù)。

(2)基于梯度的Bregman迭代法

若式(5)可微，則可對(duì)式(5)采用L-BFGS(limited-memory BFGS)算法[15]得出光吸收系數(shù)和散射系數(shù)的初值，再引入正則化，進(jìn)而進(jìn)行Bregman迭代求解。L-BFGS算法通過對(duì)吸收系數(shù)和散射系數(shù)加權(quán)，克服同時(shí)重建吸收系數(shù)和散射系數(shù)對(duì)其梯度數(shù)據(jù)的敏感度問題，所需要的存儲(chǔ)空間小。

4 定量光聲成像的主要算法

4.1 按照PAT成像光源的不同分類

根據(jù)PAT成像時(shí)所用激光光源的不同可將qPAT方法分為單光源、多光源和光譜qPAT 3類。

(1)單光源qPAT

采用單一波長(zhǎng)的單個(gè)激光光源照射組織，在已知待測(cè)組織的光散射系數(shù)分布的情況下，可以重建出光吸收系數(shù)的空間分布。若未知組織的內(nèi)在散射特性，則不能唯一、準(zhǔn)確地同時(shí)重建出光吸收參數(shù)和散射系數(shù)的分布[15]。

(2)多光源qPAT(Multiple-source，MS-qPAT)[14,28,34,46-51]

采用相同波長(zhǎng)、不同位置的多個(gè)激光光源照射組織，采集多個(gè)光聲信號(hào)數(shù)據(jù)集，進(jìn)而得出多個(gè)初始聲壓數(shù)據(jù)集。在已知待測(cè)組織的Gruneisen系數(shù)的情況下，可唯一、準(zhǔn)確地同時(shí)重建出待測(cè)組織邊界和內(nèi)部的光吸收系數(shù)和散射系數(shù)的空間分布。

(3)光譜qPAT[16,19,30,52-54]

采用多波長(zhǎng)的激光光源在不同位置照射生物組織，獲得多個(gè)初始聲壓數(shù)據(jù)集，在已知待測(cè)組織的光散射特性與入射光波長(zhǎng)之間關(guān)系的前提下，可唯一地同時(shí)重建出光吸收系數(shù)、散射系數(shù)和Gruneisen系數(shù)的空間分布。其中求解光吸收系數(shù)和散射系數(shù)是非線性問題，求解Gruneisen系數(shù)是線性問題。此類方法的缺點(diǎn)是在每次不同波長(zhǎng)的激光照明時(shí)都需要重新測(cè)量初始聲壓分布圖，計(jì)算繁瑣。

4.2 按照所用數(shù)據(jù)源的不同分類

根據(jù)重建光學(xué)特性參數(shù)時(shí)所用數(shù)據(jù)源的不同，可將qPAT算法歸結(jié)為兩步算法和一步算法兩大類，其中對(duì)兩步算法的研究和優(yōu)化是目前的主流。

4.2.1 兩步算法

兩步算法指在得到聲學(xué)逆問題結(jié)果的基礎(chǔ)上，再求解光學(xué)逆問題。即假設(shè)從探測(cè)器采集到的聲壓時(shí)間序列中重建出的光吸收能量密度或初始聲壓分布是已知的，在此基礎(chǔ)上求解光學(xué)參數(shù)的空間分布。

目前已提出的多數(shù)qPAT算法都屬于兩步算法，其基本思想如§3.2～§3.4中所述，首先建立光在組織中傳輸?shù)那跋蚰Ｐ?如RTE或者DA[21,47,55])，并對(duì)其進(jìn)行有限元離散化，然后采用適當(dāng)?shù)膬?yōu)化算法使光吸收能量密度的測(cè)量值與其理論值之間的均方誤差最小，求得光吸收系數(shù)和散射系數(shù)的空間分布。通過對(duì)光吸收能量密度進(jìn)行尺度變換(如對(duì)數(shù)化處理)[21]，可以大大提高優(yōu)化算法的收斂速度。但是qPAT測(cè)得的光強(qiáng)度的動(dòng)態(tài)范圍可能非常大，需要優(yōu)化尺度變換的比例，以確保最優(yōu)化問題數(shù)值解的穩(wěn)定性。為了得到聲速不均勻情況下的解，可采用不同波長(zhǎng)的激光對(duì)組織進(jìn)行照射，即MS-qPAT，利用采集到的多個(gè)光聲信號(hào)數(shù)據(jù)集重建光吸收系數(shù)，在這個(gè)過程中超聲波的傳播速度是可以求解的。

矢量場(chǎng)法[48,53-54]也屬于兩步算法，其基本思想是：用相同波長(zhǎng)的激光光源從兩個(gè)不同位置照射組織，收集兩個(gè)光聲信號(hào)數(shù)據(jù)。假定聲速均勻，待測(cè)區(qū)域內(nèi)組織的Gruneisen系數(shù)以及邊界處的光吸收和散射系數(shù)是已知的，待測(cè)組織內(nèi)部的光吸收系數(shù)和散射系數(shù)是有界常數(shù)且滿足Lipschitz連續(xù)，組織吸收的光子能量全部轉(zhuǎn)化為初始聲壓。通過建立基于Dirichlet邊界條件的DA模型，將根據(jù)兩個(gè)初始聲壓數(shù)據(jù)集構(gòu)建的向量場(chǎng)代入光擴(kuò)散方程，通過等量變換，求解出光吸收系數(shù)和散射系數(shù)。該算法的優(yōu)點(diǎn)是：把一個(gè)非線性問題分解成兩個(gè)線性問題，非迭代地求解光學(xué)參數(shù)，計(jì)算速度快且高效。其局限之處在于：需在待測(cè)區(qū)域的邊界光散射系數(shù)已知的情況下才能使用；只能解決滿足Dirichlet邊界條件的光散射問題；不能同時(shí)確定地重建光吸收系數(shù)、散射系數(shù)和Gruneisen系數(shù)，但假定待測(cè)區(qū)域的某一系數(shù)值已知時(shí)，可以唯一地確定另外兩個(gè)系數(shù)。

混合數(shù)值重建法[53-55]是在矢量場(chǎng)法的基礎(chǔ)上優(yōu)化而來的，原理是：假定組織的Gruneisen系數(shù)是一個(gè)不確定的常數(shù)，聲速均勻，被測(cè)組織吸收的光能量全部轉(zhuǎn)化為初始聲壓。通過擴(kuò)散光學(xué)層析(Diffuse optical tomography,DOT)成像模型獲取待測(cè)區(qū)域的邊界光吸收系數(shù)和散射系數(shù)。建立基于Dirichlet邊界條件的DA模型，采用非線性最優(yōu)化方法[36,49,55]更新邊界光吸收系數(shù)和散射系數(shù)。采用矢量場(chǎng)法求解出區(qū)域內(nèi)部的光吸收系數(shù)和散射系數(shù)，并迭代更新，直至通過RTE和DA模型生成的數(shù)據(jù)與實(shí)測(cè)數(shù)據(jù)一一匹配。該算法除了具有矢量場(chǎng)法的優(yōu)點(diǎn)以外，由于用非線性優(yōu)化問題來處理PAT逆問題，因而減小了未知空間，而利用矢量場(chǎng)法可以減小優(yōu)化范圍；結(jié)合DOT成像，求解出待測(cè)區(qū)域的邊界系數(shù)值，克服了矢量場(chǎng)法的缺陷。

此外，可以利用DOT成像測(cè)得光能流率[56]，或者利用已知光吸收譜的光學(xué)對(duì)比劑來定量測(cè)量待測(cè)組織上的局部光能流率[57]。將PAT與聲光調(diào)制結(jié)合起來，亦可去除初始聲壓測(cè)量中光能流率的影響[58]。

兩步算法的主要不足在于：根據(jù)探測(cè)器接收的光聲信號(hào)重建光吸收能量分布是存在誤差的，所以在此基礎(chǔ)上重建出的光學(xué)參數(shù)也是不準(zhǔn)確的；超聲波在組織內(nèi)的傳播速度是變化的，并不能準(zhǔn)確獲得聲速，而兩步算法僅在聲速已知的情況下重建光學(xué)參數(shù)的誤差較??；光學(xué)噪聲和光散射問題也是不可忽略的。

4.2.2 一步算法

一步算法指根據(jù)超聲探測(cè)器接收到的組織產(chǎn)生的光聲信號(hào)直接重建光學(xué)特性參數(shù)，包括一步反演算法[59]、單級(jí)法[60]和免校正法[61-62]等。

一步反演算法的基本思想是采用相同波長(zhǎng)、不同位置的光源照射待測(cè)組織，采集多個(gè)聲壓數(shù)據(jù)集 。假定已知組織的光散射系數(shù)和Gruneisen系數(shù)，將求解光吸收系數(shù)分布的逆問題分為線性和非線性兩種情況分別進(jìn)行分析。對(duì)于線性逆問題，假定已知背景組織的光吸收系數(shù)，采用Born近似法[63]和Landweber迭代法[64]，求解待測(cè)組織和背景組織的光吸收系數(shù)之間的差值，進(jìn)而重建光吸收系數(shù)分布。這種方法一般適用于聲速變化且未知、光吸收系數(shù)的變化相對(duì)較小的情況。對(duì)于非線性逆問題，采用最優(yōu)控制法，設(shè)定一個(gè)目標(biāo)函數(shù)，用Levenberg-Marquardt方法使目標(biāo)函數(shù)最小，求解光吸收系數(shù)的最優(yōu)解。在不同組織的邊界處，超聲波傳播速度和組織光學(xué)系數(shù)的相對(duì)變化一般較小，因而可以看作是線性逆問題。該算法的優(yōu)點(diǎn)是：可在超聲波速未知的情況下求解光吸收系數(shù)，同時(shí)還可以求解出超聲波的傳播速度。但是不能同時(shí)重建光吸收系數(shù)、散射系數(shù)和Gruneisen系數(shù)[52]。

單級(jí)法的基本思想是在單一波長(zhǎng)、不同位置的激光照射條件下[53]，采集超聲探測(cè)器產(chǎn)生的聲壓數(shù)據(jù)集。假設(shè)已知Gruneisen系數(shù)，且聲速均勻，使用RTE前向模型，采用有限元離散化；在包含成像目標(biāo)的封閉區(qū)域內(nèi)，使用廣義Tikhonov正則化方法[39]求解出RTE算子參數(shù)；使用近端梯度算法使Tikhonov函數(shù)最小化重建光吸收系數(shù)。在已知光散射系數(shù)的情況下，利用該方法可以唯一地重建出光吸收系數(shù)。也可在不同波長(zhǎng)的激光照明條件下，利用該算法重建光學(xué)參數(shù)。

免校正法的原理是：假設(shè)待測(cè)區(qū)域內(nèi)的介質(zhì)均勻且已知光散射系數(shù)，采用有限元離散化的方法，使用DA前向模型，根據(jù)待測(cè)邊界的聲壓測(cè)量值和計(jì)算結(jié)果得出光吸收系數(shù)的測(cè)量值；然后，使用牛頓法和Marquardt-Tikhonov正則法[42]迭代更新光吸收系數(shù)，尋找目標(biāo)函數(shù)最小化的最優(yōu)值，重建出光吸收系數(shù)分布圖。該算法的優(yōu)點(diǎn)是它僅利用聲壓的歸一化邊界測(cè)量值和入射光的強(qiáng)度，不依賴于任何校正程序；采用標(biāo)準(zhǔn)光聲數(shù)據(jù)，消除了超聲探測(cè)器靈敏度的影響。此外，還可以將該算法應(yīng)用于光譜光聲測(cè)量中。

4.2.3 小結(jié)

兩步法需要首先根據(jù)光聲信號(hào)重建組織內(nèi)部的光吸收能量分布或光聲壓分布圖，再據(jù)此重建組織的光學(xué)參數(shù)。在早期的研究中，多采用一種波長(zhǎng)的激光照射組織進(jìn)行光聲成像。在后續(xù)的研究中發(fā)現(xiàn)，不同波長(zhǎng)的激光照射組織形成的光吸收能量分布和光聲壓分布圖不同，這樣會(huì)增加重建光學(xué)參數(shù)時(shí)的變量。對(duì)于光譜qPAT問題，也可以采用兩步算法解決，分別使用多種波長(zhǎng)的激光照射組織，再重復(fù)多次重建過程，但這樣會(huì)使步驟過于繁瑣，延長(zhǎng)重建過程。目前對(duì)兩步算法的優(yōu)化主要是針對(duì)光譜qPAT或者M(jìn)S-qPAT。由于使用不同光源照射時(shí)，生物組織的光學(xué)參數(shù)是不會(huì)改變的，因而對(duì)于MS-qPAT，采用一步算法重建光學(xué)參數(shù)可解決上述兩步算法的重建過程中變量增加的問題，相對(duì)較準(zhǔn)確地重構(gòu)組織的光學(xué)參數(shù)。此外，探測(cè)器采集的光聲信號(hào)是通過測(cè)量區(qū)域的邊界部分獲得的，據(jù)此重建出的初始聲壓分布圖是不準(zhǔn)確的，進(jìn)一步重建出的光學(xué)參數(shù)也存在較大誤差，采用一步算法則可以避免重建初始聲壓分布圖時(shí)產(chǎn)生的誤差，重建出相對(duì)較準(zhǔn)確的光學(xué)參數(shù)。

5 結(jié) 論

PAT成像為研究生物組織的形態(tài)結(jié)構(gòu)、生理和病理特征以及代謝功能等提供了重要的手段。生物組織的光吸收特性與其結(jié)構(gòu)功能和病理特征等緊密相關(guān)，qPAT根據(jù)生物組織產(chǎn)生的光聲信號(hào)恢復(fù)其光學(xué)特性參數(shù)，是一個(gè)大規(guī)模、非線性、病態(tài)的逆問題。

目前，qPAT面臨的挑戰(zhàn)主要包括：(1)聲學(xué)方面，進(jìn)一步提高重建圖像的精度。(2)光學(xué)方面，在多波長(zhǎng)激光照射下實(shí)現(xiàn)參數(shù)重建，考慮未知波長(zhǎng)的依賴性影響，以及全面求解光學(xué)和聲學(xué)逆問題的規(guī)模。(3)目前的研究工作多是側(cè)重在相對(duì)理想的條件下重建組織的光學(xué)參數(shù)，需進(jìn)一步考慮實(shí)際情況的復(fù)雜性，例如組織Gruneisen系數(shù)的變化、超聲波在組織中的非勻速傳播、光在組織表面的非均勻分布以及噪聲等，使重構(gòu)出的光學(xué)參數(shù)分布圖更接近真實(shí)情況。(4)將二維圖像重建算法擴(kuò)展到三維，恢復(fù)出光學(xué)參數(shù)的三維分布圖。(5)從數(shù)學(xué)計(jì)算的角度來講，要實(shí)現(xiàn)對(duì)光吸收系數(shù)和散射系數(shù)空間分布的高分辨率重建，所需的計(jì)算成本非常高，尤其是重建光學(xué)參數(shù)的三維空間分布時(shí)，對(duì)計(jì)算機(jī)的存儲(chǔ)空間和運(yùn)行速度都有很高的要求。

除了光學(xué)參數(shù)之外，光聲成像還能利用其他成像參量進(jìn)行多尺度成像，例如：化學(xué)參量成像、黏彈參量成像、溫度參量成像、速度參量成像、聲學(xué)功率譜參量成像和物化譜參量成像等[65-66]，為疾病的臨床診治提供更為豐富的組織結(jié)構(gòu)和功能信息。

[1] YAO J,WANG L V.Breakthroughs in photonics photoacoustic tomography [J].IEEEPhoton.J.,2014,6(2):0701006.

[2] NAM S Y,CHUNG E,SUGGSL J,etal..Combined ultrasound and photoacoustic imaging to noninvasively assess burn injury and selectively monitor a regenerative tissue-engineered construct [J].TissueEng.Part C:Methods,2015,21(6):557-566.

[3] WANG L V,GAO L.Photoacoustic microscopy and computed tomography: from bench to bedside [J].Annu.Rev.Biomed.Eng.,2014,11(16):155-185.

[4] NUSTER R,GRATT S,PASSLER K,etal..Photoacoustic section imaging using an elliptical acoustic mirror and optical detection [J].J.Biomed.Opt.,2012,17(3):99-106.

[5] LAUFER J,JOHNSON P,ZHANG E,etal..Invivopreclinical photoacoustic imaging of tumor vasculature development and therapy [J].J.Biomed.Opt.,2012,17(5):449-450.

[6] ANSARI R,ZHANG E,MATHEWS S,etal..Photoacoustic endoscopy probe using a coherent fiber-optic bundle [J]SPIEBios.,2016,142 (10):97080L

[7] BILLEH Y N,LIU M,BUMA T.Spectroscopic photoacousti microscopy using a photonic crystal fiber supereontinuum source [J].Opt.Express,2010,18(18):18519-18524.

[8] MALONE E,POWELL S,COX B T,etal..Reconstruction-classification method for quantitative photoacoustic tomography [J].J.Biomed.Opt.,2015,20(12):126004.

[9] SONG H,WANG L V.Optical-resolution photoacoustic microscopy auscultation of biological systems at the cellular level [J].Biophys.J.,2013,105(4):841-847.

[10] LUTZWEILER C,DEN-BEN X L,RAZANSKY D.Expediting model-based optoacoustic reconstructions with tomographic symmetries [J].Med.Phys.2014,41(1):013302.

[11] KUCHMENT P,KUNYANSKY L.Mathematics of photoacoustic and thermoacousic tomography [J].Eur.J.Appl.Math.,2009,2(19):817-866.

[12] 孫正,韓朵朵,王健健.血管內(nèi)光聲成像圖像重建的研究現(xiàn)狀 [J].光電工程,2015,42(3):20-27.SUN Z,HAN D D,WANG J J.Review of image reconstruction for intravascular photoacoustic imaging [J].Optoelect.Eng.,2015,42(3):20-27.(in Chinese)

[13] SONG N,DEUMIC C,DAS A.Considering sources and detectors distributions for quantitative photoacoustic tomography [J].Biomed.Opt.Express,2014,5(11):3960-3974.

[14] GAO H,OSHER S,ZHAO H.MathematicalModelinginBiomedicalImagingⅡ [M] Berlin/Heidelberg: Springer,2012 131-158.

[15] GAO H,ZHAO H,OSHER S.Bregman methods in quantitative photoacoustic tomography [J].Cam.Rep.,2010,30(6):3043-3054.

[16] COX B,LAUFER J G,ARRIDGES R,etal..Quantitative spectroscopic photoacoustic imaging: a review [J].J.Biomed.Opt.,2012,17(6):ID 061202.

[17] 孫正,苑園，王健健.血管內(nèi)光聲成像的研究進(jìn)展 [J].中國(guó)生物醫(yī)學(xué)工程學(xué)報(bào),2015,34(2):221-228.SUN Z,YUAN Y,WANG J J.Progress of intravascular photoacoustic imaging [J]Chin.J.Biomed.Eng.,2015,34(2):221-228.(in Chinese).

[18] COX B T,ARRIDGE S R,KOSTLI K P,etal..Two-dimensional quantitative photoacoustic image reconstruction of absorption distributions in scattering media by use of a simple iterative method [J].Appl.Opt.,2006,45(8):1866-1875

[19] BANERJEE B,BAGCHI S,VASU R M,etal..Quantitative photoacoustic tomography from boundary pressure measurements: non-iterative recovery of optical absorption coefficient from the reconstructed absorbed energy map [J].J.Opt.Soc.Am.A.2008,25(9):2347-2356.

[20] LI S,MONTCEL B,YUAN Z,etal..Multigrid-based reconstruction algorithm for quantitative photoacoustic tomography [J].Biomed.Opt.Express,2015,6(7):2424-2434.

[21] TARVAINEN T,COX B T,KAIPIO J P,etal..Reconstructing absorption and scattering distributions in quantitative photoacoustic tomography [J].InverseProbl.,2012,28(8):1067-1079.

[22] TARVAINEN T.ComputationalMethodsforLightTransportinOpticalTomography[D].Kuopio: University of Kuopio,2006.

[23] SARATOON T,TARVAINEN T,COX B T,etal..A gradient-based method for quantitative photoacoustic tomography using the radiative transfer equation [J].InverseProbl.,2013,29(7):75006-75024.

[24] LI S,MONTCEL B,LIU W,etal..Analytical model of optical fluence inside multiple cylindrical inhomogeneities embedded in an otherwise homogeneous turbid medium for quantitative photoacoustic imaging [J].Opt.Express,2014,22(17):20500-20514.

[25] 肖嘉瑩.定量光聲成像技術(shù)及在骨關(guān)節(jié)炎診斷的研究 [D].長(zhǎng)沙:中南大學(xué),2011.XIAO J Y.QuantitativePhotoaeoustieTomographyandItsApplicationtoTheDiagnosisofOsteoarthritis[D] Changsha: Central South University,2011.(in Chinese)

[26] KUCHMENT P.TheMathematicalLegacyofLeonEhrenpreis[M].Milan: Springer,2012:183.

[27] SCHWEIGER M,ARRIDGE S R,NISSILA I.Gauss-Newton method for image reconstruction in diffuse optical tomography [J].Phys.Med.Biol.,2005,50(10):2365-2386.

[28] COX B T,TARVAINEN T,ARRIDGE S.Multiple illumination quantitative photoacoustic tomography using transport and diffusion models [C]ProceedingsofInternationalConferenceonTomographyandInverseTransportTheory,USA,MathematicalSoc.,2011.

[29] COX B T,ARRIDGE S R,BEARD P C.Gradient-based quantitative photoacoustic image reconstruction for molecular imaging [C].ProceedingsofSPIEInternationalConferenceonPhotonsPlusUltrasound:ImagingandSensing.SPIE-IntSoc..OpticalEngineering，SanJose,California,UnitedStates,2007.

[30] BAL G,REN K.On multi-spectral quantitative photoacoustic tomography in diffusive regime [J].InverseProbl.,2012,28(2):025010-25022.

[31] MAMONOV AV,REN K.Quantitative photoacoustic imaging in radiative transport regime [J].Commun.Math.Sci.,2012,12(2):201-234

[32] SOONTHORNSARATOON T.Gradient-basedMethodsforQuantitativePhotoacousticTomography[D].London: University College London,2014.

[33] NOCEDAL J.NumericalOptimization(SpringerSeriesinOperationsResearchandFinancialEngineering)[M].New York: Springer,2006:134.

[34] BAL G,Uhlmann G.Inverse diffusion theory of photoacoustics [J].InverseProbl.,2009,6(8):1407-1426.

[35] KLOSE A,HIELSCHERA H.Optical tomographic image reconstruction with quasi-Newton methods [J].InverseProbl.,2003,19(2):387-409.

[36] CEZARO A D,CEZARO F T D.Regularization approaches for quantitative photoacoustic tomography using the radiative transfer equation [J].J.Math.Anal.Appl.,2015,429(1):415-438.

[37] GAO H,ZHAO H.Multilevel bioluminescence tomography based on radiative transfer equation Part 1: l1 regularization [J].Opt.Express，2010,18(3):1854-1871.

[38] GAO H,ZHAO H.Multilevel bioluminescence tomography based on radiative transfer equation Part 2: total variation and l1 data fidelity [J].Opt.Express,2010,18(3):2894-2912.

[39] GRASMAIR M,GROSSAUER H,Haltmeier M,etal..VariationalMethodsinImaging[M].New York: Springer,2009：294.

[40] NEATER W,SCHERZER O.Quantitative photoacoustic tomography with piecewise constant material parameters [J].SIAM.J.ImagingSci.,2014,7(3):1755-1774.

[41] BERETTA E,MUSZKIETA M,NAETAR W,etal..A variational method for quantitative photoacoustic tomography with piecewise constant coefficients [J].PreprintArXiv,2015:1509.04982.

[42] CEZARO A D,LEITO A,TAI X C.On piecewise constant level-set (pcls) methods for the identification of discontinuous parameters in ill-posed problems [J].InverseProbl.,2013,29(1):015003.

[43] JIANG H,YUAN Z,YIN L,etal..Quantitative photoacoustic tomography: recovery of optical absorption coefficient maps of heterogeneous media [C].SPIE8thConferenceonBiomedicalThermoacoustics,Optoacoustics,andAcousto-optics.PhotonsPlusUltrasound:ImagingandSensing,SPIE-IntSoc.OpticalEngineering,SanJose,California,UnitedStates,2007.

[44] KALTENBACHER B,NEUBAUER A,SCHERZER O.Iterative regularization methods for nonlinear ill-posed problems [C].RadonSeriesonComputationalandAppliedMathematics，Berlin,2008.

[45] ELVETUN O L,NIELSEN B F.The split Bregman algorithm applied to PDE-constrained optimization problems with total variation regularization [J].Comput.Optim.Appl.,2016,64(3):699-724.

[46] AMMARI H,BOSSY E,JUGNON V,etal..Reconstruction of the optical absorption coefficient of a small absorber from the absorbed energy density [J].SIAMJ.Appl.Math.,2011,71(3):676-693.

[47] BERGOUNIOUX M,SCHWINDT E L.On the uniqueness and stability of an inverse problem in photo-acoustic tomography [J].J.Math.Anal.Appl.,2015,431(2):1138-1152.

[48] ZEMP R J.Quantitative photoacoustic tomography with multiple optical sources [J].Appl.Opt.,2010,49(18):3566-3572.

[49] BAL G,REN K.Multi-source quantitative photoacoustic tomography in a diffusive regime [J].Inverseprobl.,2011,27(7):75003-75022.

[50] BAL G,UHLMANN G.Inverse diffusion theory for photoacoustics [J].InverseProbl.,2010，26(8):25021-25032.

[51] SHAO P,HARRISON T,ZEMP R J.Iterative algorithm for multiple illumination photoacoustic tomography (MIPAT) using ultrasound channel data [J].Biomed.Opt.Express,2012,3(12):3240-3249

[52] COX B T,ARRIDGE S R,BEARD P C.Estimating chromophore distributions from multiwavelength photoacoustic images [J].J.Opt.Soc.Am.A,2009,26(2):443-455.

[53] BAL G,JOLLIVET A,JUGNON V.Inverse transport theory of photoacoustics [J].InverseProbl.,2010,26(2):25011-25045.

[54] BAL G,REN K.On multi-spectral quantitative photoacoustic tomography in diffusive regime [J].InverseProbl.,2012,28(2):25010-25022.

[55] REN K,GAO H,ZHAO H.A hybrid reconstruction method for quantitative PAT [J].SIAMJ.ImagingSci.2013,6(1):32-55.

[56] BAUER A Q,NOTHDURFT R E,ERPELDING TN,etal..Quantitative photoacoustic imaging: correcting for heterogeneous light fluence distributions using diffuse optical tomography [J].J.Biomed.Opt.,2011,16(9):096016.

[57] RAJIAN J R,CARSON P L,WANG X.Quantitative photoacoustic measurement of tissue optical absorption spectrum aided by an optical contrast agent [J].Opt.Express,2009,17(6):4879-4889.

[58] DaOUDI K,HUSSAIN A,HONDEBRINK E,etal..Correcting photoacoustic signals for fluence variations using acousto-optic modulation [J].Opt.Express,2012,20(13):14117-14129.

[59] DING T,REN K,VALLELIAN S.A one-step reconstruction algorithm for quantitative photoacoustic imaging [J].InverseProbl.,2015,31(9):65003-65023.

[60] HALTMEIER M,NEUMANN L,RABANSER S.Single-stage reconstruction algorithm for quantitative photoacoustic tomography [J].InverseProbl.,2015,31(6):65005-65028.

[61] YUAN Z,JIANG H B.A calibration-free,one-step method for quantitative photoacoustic tomography [J].Med.Phys.,2012,39(11):6895-6899.

[62] JIANG H B,ZHANG Q Z,YUAN Z,etal..Simultaneous reconstruction of acoustic and optical properties of heterogeneous media by quantitative photoacoustic tomography [J].Opt.Express,2006,14(15):6749-6754.

[63] RONDI L,SANTOSA F.Enhanced electrical impedance tomographyviathe mumford-shah functional [J].EsaimContr.Optim.,2000,6(6):517538.

[64] KIRSCH A.AnIntroductiontoTheMathematicalTheoryofInverseProblems[M].New York: Springer,1996:234

[65] 殷杰，陶超，劉曉峻.多參量光聲成像及其在生物醫(yī)學(xué)領(lǐng)域的應(yīng)用 [J].物理學(xué)報(bào),2015,64(9):098102.YIN J,TAO C,LIU X J.Multi-parameter photoacoustic imaging and its application in biomedicine [J].ActaPhys.Sinica,2015,64(9):098102.(in Chinese)

[66] 苗少峰,楊虹,黃遠(yuǎn)輝,等.光聲成像研究進(jìn)展 [J].中國(guó)光學(xué),2015,8(5):699-713.MIAO S F,YANG H,HUANG Y H,etal..Research progresses of photoacoustic imaging [J].Chin.Opt.,2015,8(5):699-713.(in Chinese)

孫正(1977-)，女，河北保定人，博士，教授，碩士生導(dǎo)師，2004年于天津大學(xué)獲得博士學(xué)位，主要從事生物醫(yī)學(xué)信號(hào)處理、醫(yī)學(xué)成像和圖像處理等方面的研究。

E-mail: sunzheng_tju@163.com

Review on Progress of Quantitative Photoacoustic Tomography

SUN Zheng*,ZHENG Lan

(DepartmentofElectronicandCommunicationEngineering,NorthChinaElectricPowerUniversity,Baoding071003,China)

Photoacoustic tomography (PAT),an emerging medical imaging modality,combines the high resolution of ultrasonic imaging and high contrast of optical imaging.Current research on PAT includes two inverse problems,i.e.,constructing the distribution of initial acoustic pressures according to the photo-acoustic signals generated by the tissues and estimating the optical absorption and scattering coefficients of the tissues within the imaging region based on the results of the first inversion.The latter,known as quantitative photoacoustic tomography (qPAT),is in general a nonlinear ill-posed problem.This paper summarizes current algorithms for solving the qPAT inversion.Related advantages and limits as well as perspective studies are discussed.

quantitative photoacoustic tomography (qPAT); image reconstruction; inverse problem; optical absorption coefficient; diffusion coefficient; Gruneisen coefficient

1000-7032(2017)09-1222-11

2017-03-24；

2017-04-27

國(guó)家自然科學(xué)基金(61372042)； 中央高?；究蒲袠I(yè)務(wù)費(fèi)專項(xiàng)資金(2014ZD31)資助項(xiàng)目

O439

A

10.3788/fgxb20173809.1222

*CorrespondingAuthor,E-mail:sunzheng_tju@163.com

Supported by National Natural Science Foundation of China(61372042); Special Fund for Basic Scientific Research Business in Central Universities(2014ZD31)