不同圖像重建算法對定量化光聲重建效果的影響

宋寧寧，秦航，蔣紅兵，余晶

1.南京醫(yī)科大學附屬南京醫(yī)院（南京市第一醫(yī)院） 醫(yī)療設備處，江蘇 南京 210006；2.南京市衛(wèi)生信息中心，江蘇 南京 210003

引言

光聲成像結(jié)合了光學成像高對比度和超聲成像高分辨度的優(yōu)點。光聲成像的原理是用光照射生物組織，生物組織吸收一部分能量并轉(zhuǎn)化成熱能，使生物組織的溫度升高，引起熱膨脹并產(chǎn)生超聲波。生物組織吸收的能量正比于光子數(shù)密度和生物組織的吸收系數(shù)。光聲重建的第一步是用采集到的超聲波重建初始場強，初始場強是Grüneisen系數(shù)（反映生物組織的聲學特性）、光子數(shù)密度和生物組織的光吸收系數(shù)三者的乘積，只有生物組織的光學系數(shù)才能反映生物組織的本質(zhì)，所以光聲成像的第二步是分離出生物組織的光學系數(shù)（如光的吸收系數(shù)和散射系數(shù)）和聲學系數(shù)（Grüneisen系數(shù)等），光聲重建的第二步也稱定量化光聲重建。

對于定量化光聲重建，已經(jīng)有大量的研究成果，我們大體上歸為4類：

（1）重建光的吸收系數(shù)。假設光的散射系數(shù)是已知的，Ripol等[1]用一個點光源重建了光的吸收系數(shù)，考慮了光源的長度，光源脈沖的形狀，脈沖響應和超聲換能器的尺寸有限性。Banerjee等[2]從測量到的邊界超聲直接重建了光的吸收系數(shù)。

（2）重建光的吸收系數(shù)和光子數(shù)密度。Rosenthal等[3]提取了光的吸收系數(shù)和光子數(shù)密度，這個算法不是基于光的光擴散方程的具體解，所以對光源的形狀沒有具體要求。

（3）通過一個不同波長的光源同時重建光的吸收系數(shù)和散射系數(shù)，Guillaume等[4]已經(jīng)證明，用一個不同波長的光源可以同時重建光的吸收系數(shù)，散射系數(shù)和Grüneisen系數(shù)。

（4）通過多個光源同時重建光的吸收系數(shù)和散射系數(shù)。自上世紀以來，學者們做了大量的工作旨在提高光聲重建的定量化水平。通過從初始場強的分布中重建出生物組織的不同物理參數(shù)：光的吸收系數(shù)，光的散射系數(shù)和Grüneisen系數(shù)。Shao等[5]用多個光源從重建的初始場強中分離出光的吸收系數(shù)，光的散射系數(shù)和Grüneisen系數(shù)；在他們最新的工作中，直接從測量到的超聲波中分離出了光的吸收系數(shù)和散射系數(shù)。

本研究直接從測量到的超聲波中分離出了光的吸收系數(shù)和散射系數(shù)，同時考慮了不同圖像重建算法對定量化光聲重建效果的影響。

1　光聲成像物理理論

下面我們回顧一下光聲成像過程中的物理定義和光聲成像中的主要公式，包括光聲成像的正向過程和逆向過程（光聲重建過程）。

在光聲成像中，入射光的脈寬必須小于生物組織的熱弛豫時間和壓力擴散時間，人們把這兩個限制條件稱之為熱限制條件和壓力限制條件。只有滿足這兩個條件，在激光脈沖周期內(nèi)，生物組織吸收光能產(chǎn)生的熱量來不及外散只能用于引起自身的溫度升高，同時，溫度升高引起的熱彈性膨脹壓力也來不及外傳，只能形成壓力波向外傳輸[6-9]。滿足這兩個條件的光脈沖時間通常只有幾個納秒，所以，熱函數(shù)可以寫成一個正比于Dirac函數(shù)的瞬時函數(shù)。在光聲成像中，生物組織吸收的熱量和光照強度是成比例的。

其中，H(r,t)是熱函數(shù)，定義為在空間r，時間點t轉(zhuǎn)化成的熱能。(r)是光子數(shù)密度。μa(r)是吸收系數(shù)。光子數(shù)密度可通過求解光擴散方程得到。

D(r)是擴散系數(shù)，Q(r)是由脈沖激光的照射所產(chǎn)生的光源分布。生物組織吸收的熱量轉(zhuǎn)化成熱能，組織的溫度上升，產(chǎn)生熱膨脹，初始場強產(chǎn)生[10-16]：p0(r,u)=Γ(r)μa(r)(r,u)，Γ(r)是Grüneisen系數(shù)，表示了生物組織吸收的能量轉(zhuǎn)化為聲壓的效率[10-16]。

式(3)是外傳的超聲波的表達式。為了避免計算對超聲波對時間的偏導，我們計算的pmodel(r,u,t)是超聲波對時間的積分乘以聲速的平方。

其中u是生物組織的光學參數(shù)。u=[μaD]T。對于有異常的生物組織來說，生物組織的光學參數(shù)可以寫成u=u0+δu，其中u0是正常組織的光學參數(shù)，δu是在正常組織光學參數(shù)的基礎上疊加的光學參數(shù)的異常。如果δu＜＜u0，我們可以把測量到的超聲波表達式展開成一階泰勒級數(shù)的形式。

這里Δ=p-pmodel，H是Hessian矩陣。

2　仿真實驗及結(jié)果

我們以人體乳腺組織為例，進行了下面的仿真實驗（圖1），選取了圖1所示6 cm×6 cm的一個方形測試區(qū)域，其中位于中心的2 cm×2 cm的區(qū)域用來做光聲重建。重建區(qū)域正常組織的吸收系數(shù)和散射系數(shù)分別是μa0=0.1 cm-1和D0=0.03 cm。在正常組織的基礎上加入了一個擴散系數(shù)的異常δD=0.003 cm（位于中心，尺寸是0.5 cm×0.5 cm），兩個吸收系數(shù)的異常δμa=0.01 cm-1（尺寸分別是0.3 cm×1.1 cm和0.3 cm×0.3 cm，偏離中心0.6 cm）。選擇的異常相對較?。ㄕ＝M織的10%以內(nèi)），這樣能夠滿足Born approximation。如果異常相對正常組織較大，需要用非線性重建算法。光源位于重建區(qū)域0.3 cm處，超聲傳感器位于重建區(qū)2 cm處，即測試區(qū)域的邊緣。

我們測試了兩種不同的圖像重建算法對重建效果的影響：奇異值分解算法（Singular Value Decomposition，SVD）和代數(shù)迭代算法（Algebraic Reconstruction Technique，ART）。

在奇異值分解算法中，Hessian矩陣分解成：H=U×S×VT，S是和H同維的對角矩陣，其中非負奇異值iσ是以降序排序的。U和V是單位矩陣。為了避免稀疏矩陣H引起的不穩(wěn)定性，采用Tikhonov正則化，用代替了S-1的對角值。λ是Tikhonov正則化算子。

代數(shù)迭代重建算法用迭代方法尋求δu的解：

式(7)中，λn是一個加速算子，在下面所有的仿真中，λn取常數(shù)值1。

圖1　測試區(qū)域

在下面的仿真實驗中，我們用了4個點光源和60個超聲傳感器，超聲傳感器均勻的分布在測試區(qū)域的邊緣（一側(cè)15個傳感器）。測試了兩種不同的圖像重建算法：SVD和ART在不同的噪聲水平下對定量化圖像重建效果的影響。在不同的噪聲水平下重建的吸收系數(shù)和散射系數(shù)的異常，見圖2。

圖2　不同算法下重建的光的吸收系數(shù)和散射系數(shù)

從圖2我們可以看出，重建結(jié)果對噪聲是非常敏感的，不同的噪聲水平下重建效果差別很大，在一定噪聲水平范圍內(nèi)（本文研究的噪聲范圍均≤1%），隨著噪聲的增大，ART的重建效果似乎更好一些，不同噪聲水平下重建圖像和原圖像之間的平方差，見表1。在一定噪聲水平范圍內(nèi)，ART的重建效果要比SVD好。在以后的重建中，我們選擇ART作為重建算法。

表1　不同算法和噪聲水平下重建圖像和原圖像之間的平方差

3　結(jié)論

我們用了多個光源來做定量化光聲重建，此方法可以同時重建光的吸收系數(shù)和散射系數(shù)。在定量化光聲重建過程中，我們測試了兩種不同的線性圖像重建算法：SVD和ART。測試結(jié)果表明，將一定噪聲水平融進去，ART比SVD效果好。同時，從結(jié)果中我們可以看出，這兩種線性算法對噪聲都是非常敏感的，在噪聲進一步加大的情況下，我們應考慮其他非線性算法來提高重建圖像質(zhì)量。

[參考文獻]

[1] Ripoll J,Ntziachristos V.Quantitative point source photoacoustic inversion formulas for scattering and absorbing media[J].Phys Rev E,2005,71.

[2] Banerjee B,Bagchi S,Vasu RM.Quantitative photoacoustic tomography from boundary pressure measurements:noniterative recovery of optical absorption coef fi cient from the reconstructed absorbed energy map[J].J Opt Soc Am,2008,25:2347-2356.

[3] Rosenthal A,Razansky D,Ntziachristos V.Quantitative optoacoustic signal extraction using sparse signal represen-tation[J].Med Imaging IEEE Trans,2009,28:1997-2006.

[4] Guillaume B,Ren K.On multi-spectral quantitative photoacoustic tomography in diffusive regime[J].Inverse Probl Imag,2012.

[5] Shao P,Cox B,Zemp RJ.Estimating optical absorption,scattering, and Grueneisen distributions with multipleillumination photoacoustic tomography[J].Appl Opt,2011,50:3145-3154.

[6] 肖嘉瑩.定量光聲成像技術(shù)及在骨關節(jié)炎診斷的研究[D].長沙:中南大學,2011.

[7] 向良中,邢達,谷懷民,等.改進的同步迭代算法在光聲血管成像中的應用[J].物理學報,2007,56(7):3911-3916.

[8] 吳丹,陶超,劉曉峻.光聲成像中延遲求和方法和反重影重構(gòu)方法的比較[J].無損檢測,2011,33(9):36-39.

[9] 楊迪武,黃仲,曾呂明,等.基于環(huán)形陣列探測器的快速光聲成像[J].激光生物學報,2015,24(5):423-427.

[10] 簡曉華,崔崤峣,向永嘉,等.自適應多光譜光聲成像技術(shù)研究[J].物理學報,2012,61(21):456-461.

[11] 黃盛松,劉博,吳登龍.光聲成像在泌尿系疾病的研究進展[J].外科研究與新技術(shù),2016,5(1):49-52.

[12] 關天培,方馳華.光聲成像技術(shù)及其在原發(fā)性肝癌邊界界定中的應用[J].中華肝臟外科手術(shù)學電子雜志,2016,5(2):65-67.

[13] 吳寧,任秋實,李長輝.光聲層析成像研究進展[J].中國醫(yī)療設備,2015,3(2):16-20.

[14] 孫正,苑園,王健健.血管內(nèi)光聲成像的研究進展[J].中國生物醫(yī)學工程學報,2015,(2):221-228.

[15] 翁國星,陳遠翔,唐嘉銘,等.光聲成像技術(shù)評估慢性心肌缺血的活體研究.中國激光醫(yī)學雜志,2016,25(5):256-259.

[16] Herzog E,Taruttis A,Beziere N,et al.多譜線光聲斷層攝影術(shù)的癌癥異質(zhì)性光學成像[J].中國醫(yī)療設備,2012,27(9):21-26.