焦陽 簡小華 向永嘉 崔崤峣

(中國科學院蘇州生物醫(yī)學工程技術(shù)研究所,蘇州 215163)

(2012年9月10日收到;2012年12月21日收到修改稿)

1 引言

早在十九世紀,Bell[1]就發(fā)現(xiàn)了光聲效應(photoacoustic effect),但直到近20年,隨著光纖、激光技術(shù)、壓電陶瓷技術(shù)以及計算機技術(shù)的發(fā)展,針對光聲效應的研究和應用才受到廣泛關(guān)注[2-5].光聲成像(photoacoustic imaging)就是利用光聲效應,以脈沖光作為激勵源、超聲信號作為載體的一種功能成像方法.該方法有機地結(jié)合了光學成像和聲學成像的優(yōu)點[6],具有無損、安全、成本低、且可以提供深層組織的高分辨率和高對比度的無損檢測[7].20世紀90年代,Kruger等[8,9]提出類似CT成像方法的逆三維Radon變換反投影重建算法,首次將光聲效應引入了生物組織成像領域.荷蘭Twente大學的Hoelen研究組[10]采用非聚焦的超聲傳感器以橫向線性掃描方式采集光聲信號,對模擬毛細血管樣品進行成像并引入延遲疊加的概念.2003年Wang等[11-14]利用光聲成像技術(shù)對老鼠腦部血管進行成像,提出給予精確解的重建算法,得到了與解剖效果一致的清晰腦部血管分布.中國華南師范大學邢達研究組[15-17]采用320陣元的醫(yī)用B超儀線陣探頭,以相控聚焦方式對光聲信號進行電子線掃采集,提高了數(shù)據(jù)采集與成像速度.上述這些光聲成像研究[18]多是利用光聲信號的光譜分析得到組織內(nèi)部的光吸收系數(shù),這種方法需要已知組織的光譜,并通過多光譜成像[19].少有研究者通過光聲信號頻譜特性反推出組織頻率特性,從而區(qū)分不同組織成分及結(jié)構(gòu)[20].為此,本文提出了光聲信號的雙譜分析方法,包括對組織的光聲光譜分析和光聲頻譜分析兩部分.兩者在區(qū)分組織的過程中互為參照,可進一步提高光聲成像的判別組分及結(jié)構(gòu)的準確性,提高不同組織間的對比度,以彌補單一光譜下光聲成像對組織識別的局限性[21].如在750 nm波長下,脂肪和血紅蛋白光吸收系數(shù)接近[22,23],難以通過該波長下光譜特性對兩者進行識別.對于光聲信號的頻譜分析還可以對信號的頻域信息作定量解釋,有助于研究或證明光聲效應的產(chǎn)生機理、量化光聲成像技術(shù),并以此為據(jù),還可進一步對不同組織的光聲頻域功能成像進行基礎理論研究與未來應用探索.

2 理論分析

2.1 光聲效應理論分析

光聲信號的實質(zhì)是超聲信號,對于超聲信號,聲壓是目前最為普遍采用的描述聲波性質(zhì)的物理量.在靜止、無源、各向同性的均勻理想介質(zhì)下,當脈沖激光照射在介質(zhì)上,由于介質(zhì)的吸收效應導致光能量的沉積和溫度的增加而產(chǎn)生光聲信號.在脈沖激光作用下,介質(zhì)溫度T變化最終導致壓力變化,此時,聲壓p可以表示為[24]

式中,t和r分別為時間與位置,β為熱膨脹系數(shù),ρ為介質(zhì)密度,vs代表聲波在介質(zhì)中傳播速度,T代表介質(zhì)溫度.等式右邊第一項為密度變化對聲壓的影響,第二項為溫度變化對聲壓的影響.

當光強為I0的脈沖激光照射組織時,其產(chǎn)生的聲壓為

式中φ(r)表示光能量吸收分布,η(t)是脈沖激光的時間函數(shù)[25].依據(jù)(2)式,最終通過超聲換能器測得實際聲壓數(shù)據(jù),反推出光能吸收分布,達到成像目的.

2.2 光聲信號振動理論分析

光聲信號在組織中的傳播過程可以描述為無限自由度非保守系統(tǒng)在瞬時激勵下的非線性振動問題.其具體過程為：組織受激光照射發(fā)生光聲效應,產(chǎn)生的光聲信號在組織中傳播并衰減.衰減過程包括兩部分：其一是光聲信號在組織中傳播時,組織內(nèi)部產(chǎn)生交變的應變,通過組織內(nèi)摩擦而消耗能量,這種阻尼稱為結(jié)構(gòu)阻尼;其二是光聲信號在組織傳播過程中,組織與介質(zhì)(低黏度流體)發(fā)生摩擦,產(chǎn)生非線性阻尼.本文涉及的實驗中,超聲探測器采用緊貼組織的側(cè)向接收模式,暫不考慮光聲信號在介質(zhì)內(nèi)的衰減.

非保守系統(tǒng)由阻尼耗能導致的非保守系統(tǒng)振動微分方程為[26]

式中d1為組織與介質(zhì)摩擦產(chǎn)生的阻尼力,d2為組織結(jié)構(gòu)阻尼產(chǎn)生的阻尼力.其中阻尼力d1隨速度的變化為

式中,β1為低黏度阻尼系數(shù).當n=2時,-d(u˙)是低黏度流體非線性阻尼力,適用于描述物體在空氣或低黏度液體中做中高速運動所受到的阻力[26],即組織與介質(zhì)摩擦產(chǎn)生阻尼力.復雜的非線性阻尼的作用是耗能,從而可通過能量法將非線性阻尼等效為線性黏性阻尼,等效阻尼系數(shù)[26]：

其中,η1為系統(tǒng)損耗因子,k1為系統(tǒng)彈性特性,ω1是系統(tǒng)固有頻率,上式中損耗因子表示為

其中W是阻尼力在一個振動周期T內(nèi)消耗的能量,Vmax為系統(tǒng)最大彈性勢能,等于無阻時的系統(tǒng)最大彈性勢能.一個周期內(nèi),低黏度流體阻尼所消耗的能量為

將(7)式代入(5)式得到阻尼系數(shù)

另一方面,d2代表的結(jié)構(gòu)阻尼系統(tǒng)動響應的主要譜成分一般落在系統(tǒng)的共振頻段或激勵中主要譜成分所在頻段.因此阻尼系統(tǒng)的響應可以對每一頻段用一個黏性阻尼來等效該頻段上的結(jié)構(gòu)阻尼[27]：

其中ωr是該頻段的中心頻率.

試驗證明,在很寬的頻率范圍內(nèi),結(jié)構(gòu)阻尼所消耗的能量W2幾乎與振動頻率無關(guān),僅與振幅平方成正比[26]：

式中,由于結(jié)構(gòu)阻尼造成應力滯后于應變,β2亦可稱為遲滯阻尼系數(shù),aˉ是衰減過程中的等效平均振幅.把(6)和(10)式代入(9)式,得到

數(shù)值計算表明,在η≤0.4時這種近似對單位初速度所激發(fā)自由振動的固有頻率及振幅引起的相對誤差分別為2.02%和1.4%[27].

將(8)和(11)式帶入振動微分方程(3),得到最終近似后的光聲信號振動微分方程為

此系統(tǒng)所描述的范圍僅為組織樣品表面附近直接受激光照射的微小區(qū)域.(12)式為系統(tǒng)各頻段中心頻率 ωr(r=1,2,···)的隱式表達式,結(jié)合系統(tǒng)固有特性,可以得到系統(tǒng)的質(zhì)量、剛度矩陣、最大彈性勢能Vmax;根據(jù)組織及介質(zhì)的固有阻尼特性,可以得到低黏度阻尼系數(shù)β1和遲滯阻尼系數(shù)β2.最后求解振動方程可以得到組織樣品受光照區(qū)域的固有頻率理論解,具體求解過程可參照文獻[26]中無限自由度系統(tǒng)的振動方程求解方法,這里不作贅述.值得注意的是：根據(jù)光聲效應理論分析,不同激光波長對應組織不同的光吸收系數(shù),而最終導致光聲信號幅值不同,與頻率無關(guān).

3 實驗分析

3.1 聲效應實驗系統(tǒng)

圖1是我們進行光聲信號雙譜分析的實驗系統(tǒng)示意圖.其基本原理如下：OPO激光器(PhocusTM,OPOTEk Inc.)發(fā)出的脈沖激光經(jīng)過擴束鏡(兩個D=25.4 mm,F=25.4 mm雙凸透鏡,ThorLabs Inc.)擴束后,經(jīng)中性密度衰減片(mounted variable circular ND,D=50 mm,ThorLabs Inc.)衰減后耦合入光纖中,再由光纖引導照射到實驗樣品上,產(chǎn)生的光聲信號由側(cè)向緊貼實驗樣品的超聲換能器(V307-SU水浸探頭,Olympus Inc.)探測接收,探測到的光聲信號經(jīng)放大器(OLYMPUS-5678,帶寬50 kHz—40 MHz,增益40 dB)放大后接入脈沖函數(shù)發(fā)生接收器(DPR-500,JSR Ultrasonics),再由示波器(數(shù)字熒光示波器DP05034,Tek.Inc.最高采樣率為5 Gs/s,帶寬為500 MHz)顯示采集到的信號,并保存至計算機中.激光器波長變化范圍：680—950 nm,實驗從680 nm波長的激光束開始,每次實驗波長增加20 nm,其他各實驗參數(shù)、條件不變,以此研究不同波長的激光束對實驗樣品的光聲信號的影響.實驗中,分別取完整豬肝的四分之一和相似體積的豬皮下脂肪作為實驗樣品.

圖1 實驗系統(tǒng)示意圖

3.2 實驗結(jié)果光譜分析

對實驗獲得的不同波長激光激勵下產(chǎn)生的多組光聲信號取絕對值,并做平均,根據(jù)之前的理論分析,該值可以用來表征該波長激光產(chǎn)生的光聲信號的強度,實驗得到的光聲光譜與國外相關(guān)結(jié)果進行比較如圖2和圖3.

圖2 豬肝實驗樣品光譜與國外類似工作對比 (a)Prahl等實驗所得血紅蛋白的光譜[23];(b)豬肝實驗樣品光譜

由以上實驗結(jié)果可知,豬肝實驗樣品的光譜與去氧血紅蛋白光譜類似,豬肝實驗樣品血液含量很高,而血紅蛋白的光吸收系數(shù)較高.因此,我們推測豬肝實驗樣品的光聲信號主要是血紅蛋白受激光照射激發(fā)而出的.而豬皮下脂肪實驗樣品的光譜與國外類似結(jié)果大致相同.

圖3 豬皮下脂肪樣品光譜與國外類似工作對比 (a)Michels等實驗所得脂肪的光譜[22];(b)豬皮下脂肪實驗樣品光譜

圖4 豬肝對于不同波長的激光激勵下的光聲信號功率譜

3.3 實驗結(jié)果頻譜分析

3.3.1光聲信號的頻率估計

首先,需要對所得到的光聲信號其大致頻率進行估計,確定實驗仿體的固有頻率,以設計相應的濾波器以降低噪聲影響,提高信號信噪比上述實驗中豬肝在680—950 nm波長激光激勵下的光聲信號功率譜如圖4所示.

通過對其光聲信號功率譜的分析,不難發(fā)現(xiàn)其功率譜峰值對應頻率大致集中在30—90 kHz之間,以此設計相應的帶通濾波器,以提高光聲信號的信噪比,具體過程如下.

3.3.2光聲信號的信號處理

1)實驗采集到的原始數(shù)據(jù)(豬肝在950 nm波長激光照射下產(chǎn)生的光聲信號,此時激光器能量衰減最大,光聲信號信噪比最低)如圖5(a)所示;

圖5 光聲信號信號處理 (a)示波器采集的光聲信號;(b)濾波后的光聲信號

2)通過濾波器濾波得到光聲信號如圖5(b)所示;

3)對濾除背景噪聲的光聲信號圖5(b)用窗函數(shù)切趾后,通過快速傅里葉變換,即可得到光聲信號的功率譜,見圖6;

4)取出各信號功率譜最高峰值及次高峰值對應的頻率,統(tǒng)計后得到兩種實驗樣品的光聲信號頻率,見圖7.由圖可知,實驗豬肝樣品的低頻段固有頻率為47和75 kHz;豬皮下脂肪的低頻段固有頻率主要集中在30 kHz.兩者頻率有明顯差異,且與激光波長無關(guān),與理論結(jié)果一致.

圖6 濾波后光聲信號的功率譜

圖7 實驗樣品功率譜峰值對應頻率統(tǒng)計結(jié)果 (a)豬肝實驗樣品頻率分布統(tǒng)計;(b)豬皮下脂肪實驗樣品頻率分布統(tǒng)計

4 結(jié)論

本文首次結(jié)合光聲效應和振動理論分析,對不同實驗樣品在短脈沖激光照射下產(chǎn)生的光聲信號進行雙譜分析,發(fā)現(xiàn)光聲光譜信號的差異主要源于組織內(nèi)不同組分光吸收特性的不同,因此其可以用來進行組織組分、邊界識別等;而在光聲頻譜方面,實驗表明豬肝的固有頻率集中在47和75 kHz;豬皮下脂肪的低頻段固有頻率主要集中在30 kHz.兩者頻率有明顯差異,且與激光波長無關(guān),說明通過光聲頻譜分析區(qū)別不同物質(zhì)具有可行性,與理論分析一致.實驗結(jié)果表明光聲信號雙譜分析法可進一步增強光聲信號對組織的識別能力.光聲信號的光譜與頻譜分析方法還可以用于探索光聲頻率探傷等技術(shù).[28],在臨床診斷中具有重要的理論研究意義和廣闊的應用前景.

[1]Bell A G 1980 Am.J.Sci.Arts Ser.3 305

[2]Li J 2010 Modern Instrum.16 31(in Chinese)[李軍2010現(xiàn)代儀器16 31]

[3]Wu D,Tao C,Liu X J2010 Acta Phys.Sin.59 5850(in Chinese)[吳丹,陶超,劉曉峻2010物理學報59 5850]

[4]Xiang L Z,Xing D,Gu H M,Yang SH,Zeng L M 2007 Acta Phys.Sin.56 3916(in Chinese)[向良忠,邢達,谷懷民,楊思華,曾呂明2007物理學報56 3916]

[5]Jian X H,Cui Y Y,Xiang Y J,Han Z L 2012 Acta Phys.Sin.61 217801(in Chinese)[簡小華,崔崤峣,向永嘉,韓志樂2012物理學報61 217801]

[6]Li C,Wang L V 2009 Phys.Med.Biol.54 R59

[7]Xu X H,Li H 2008 Acta Phys.Sin.57 4628(in Chinese)[徐曉輝,李暉2008物理學報57 4628]

[8]Kruger RA,Liu P,Appledorn CR 1995 Med.Phys.22 1605

[9]Kruger RA,Stantz K,Kiser Jr WL 2002 Proc.SPIE 521

[10]Hoelen C,De Mul F,Pongers R,Dekker A 1998 Opt.Lett.23 648

[11]Ku G,Wang L V 2005 Opt.Lett.30 507

[12]Ku G,Wang X,Stoica G,Wang L V 2004 Phys.Med.Biol.49 1329

[13]Ku G,Wang X,Xie X,Stoica G,Wang L V 2005 Appl.Opt.44 770

[14]Xu M,Wang L V 2006 Rev.Sci.Instrum.77 041101

[15]Yang D,Xing D,Gu H,Tan Y,Zeng L 2005 Appl.Phys.Lett.87 194101

[16]Yin B,Xing D,Wang Y,Zeng Y,Tan Y,Chen Q 2004 Phys.Med.Biol.49 1339

[17]Zeng Y G,Xing D,Fu H B,Wang Y 2005 Chin.J.Lasers 32 97(in Chinese)[曾亞光,邢達,付洪波,王毅2005中國激光32 97]

[18]Yang D W,Xing D,Zhao X H 2010 Chin.Phys.Lett.27 144

[19]Zhang J,Yang SH 2010 Chin.J.Lasers 38 104001(in Chinese)[張建,楊思華2010中國激光38 104001]

[20]Tan Y,Xing D,Wang Y,Zeng Y G,Yi B Z 2005 Acta Photon.Sin.34 1019(in Chinese)[譚毅,邢達,王毅,曾亞光,尹邦政2005光子學報34 1019]

[21]Qian SY,Xing D 2000 Acta Laser Biol.9 228(in Chinese)[錢盛友,邢達2000激光生物學報9 228]

[22]Michels R,Foschum F,Kienle A 2008 Opt.Express16 5907

[23]Prahl S 1999 Oregon Medical Laser Center 15 277

[24]Hoelen C,de Mul FFM 1999 J.Acoust.Soc.Am.106 695

[25]Song Z Y 2009 Ph.D.Dissertation(Tianjin：Tianjin University)(in Chinese)[宋智源2009博士學位論文(天津：天津大學)]

[26]Hu HY,Sun JH,Chen HH 1998 Mechanical Vibration and Shock(1st Ed.)(Beijing：Aviation Industry Press)(in Chinese)[胡海巖,孫久厚,陳懷海1998機械振動與沖擊(第1版)(北京：航空工業(yè)出版社)]

[27]Hu H Y 1992 J.Vib.Eng.5 8(in Chinese)[胡海巖1992振動工程學報5 8]

[28]Fan H D,Fu JM,Zhou X M,Zhao H 2004 Nondestruc.Test.26 551(in Chinese)[范海東,傅金明,周小明,趙晗2004無損檢測26 551]