光聲信號的雙譜分析方法研究*

焦陽 簡小華 向永嘉 崔崤峣

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

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

1 引言

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

2 理論分析

2.1 光聲效應(yīng)理論分析

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

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

當(dāng)光強(qiáng)為I0的脈沖激光照射組織時(shí),其產(chǎn)生的聲壓為

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

2.2 光聲信號振動(dòng)理論分析

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3 實(shí)驗(yàn)分析

3.1 聲效應(yīng)實(shí)驗(yàn)系統(tǒng)

圖1是我們進(jìn)行光聲信號雙譜分析的實(shí)驗(yàn)系統(tǒng)示意圖.其基本原理如下：OPO激光器(PhocusTM,OPOTEk Inc.)發(fā)出的脈沖激光經(jīng)過擴(kuò)束鏡(兩個(gè)D=25.4 mm,F=25.4 mm雙凸透鏡,ThorLabs Inc.)擴(kuò)束后,經(jīng)中性密度衰減片(mounted variable circular ND,D=50 mm,ThorLabs Inc.)衰減后耦合入光纖中,再由光纖引導(dǎo)照射到實(shí)驗(yàn)樣品上,產(chǎn)生的光聲信號由側(cè)向緊貼實(shí)驗(yàn)樣品的超聲換能器(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)顯示采集到的信號,并保存至計(jì)算機(jī)中.激光器波長變化范圍：680—950 nm,實(shí)驗(yàn)從680 nm波長的激光束開始,每次實(shí)驗(yàn)波長增加20 nm,其他各實(shí)驗(yàn)參數(shù)、條件不變,以此研究不同波長的激光束對實(shí)驗(yàn)樣品的光聲信號的影響.實(shí)驗(yàn)中,分別取完整豬肝的四分之一和相似體積的豬皮下脂肪作為實(shí)驗(yàn)樣品.

圖1 實(shí)驗(yàn)系統(tǒng)示意圖

3.2 實(shí)驗(yàn)結(jié)果光譜分析

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

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

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

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

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

3.3 實(shí)驗(yàn)結(jié)果頻譜分析

3.3.1光聲信號的頻率估計(jì)

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

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

3.3.2光聲信號的信號處理

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

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

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

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

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

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

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

4 結(jié)論

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

[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物理學(xué)報(bào)59 5850]

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

[5]Jian X H,Cui Y Y,Xiang Y J,Han Z L 2012 Acta Phys.Sin.61 217801(in Chinese)[簡小華,崔崤峣,向永嘉,韓志樂2012物理學(xué)報(bào)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物理學(xué)報(bào)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)[曾亞光,邢達(dá),付洪波,王毅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)[譚毅,邢達(dá),王毅,曾亞光,尹邦政2005光子學(xué)報(bào)34 1019]

[21]Qian SY,Xing D 2000 Acta Laser Biol.9 228(in Chinese)[錢盛友,邢達(dá)2000激光生物學(xué)報(bào)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博士學(xué)位論文(天津：天津大學(xué))]

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

[27]Hu H Y 1992 J.Vib.Eng.5 8(in Chinese)[胡海巖1992振動(dòng)工程學(xué)報(bào)5 8]

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