基于光聲差分衰減頻譜的骨質(zhì)評(píng)估仿真研究

余東，解維婭，封婷，程茜

(1.同濟(jì)大學(xué)物理科學(xué)與工程學(xué)院聲學(xué)研究所，上海 200092；2.南京理工大學(xué)電子工程與光電技術(shù)學(xué)院，江蘇南京 210094)

0 引 言

生物醫(yī)學(xué)光聲是近年來(lái)發(fā)展起來(lái)的一種非侵入式和非電離式的新型生物組織檢測(cè)方法[1-2]。當(dāng)不同波長(zhǎng)的脈沖激光照射到生物組織時(shí)，組織中不同生物分子對(duì)光的吸收將產(chǎn)生超聲信號(hào)，這種由光激發(fā)產(chǎn)生的超聲信號(hào)被稱為光聲信號(hào)[3]。光聲技術(shù)在評(píng)估骨質(zhì)疏松癥方面具有顯著優(yōu)勢(shì)，它既能提供骨礦物質(zhì)、骨髓、血紅蛋白、膠原蛋白等分子信息又能評(píng)估骨骼微結(jié)構(gòu)信息和力學(xué)狀態(tài)，在疾病的診斷中顯示了巨大的潛力。

光聲技術(shù)在骨質(zhì)評(píng)估中的研究已有一些成果。密歇根大學(xué)的Wang課題組研究了使用熱光聲測(cè)量和光聲光譜分析兩種技術(shù)評(píng)估大鼠小梁骨的骨礦物質(zhì)密度和骨微結(jié)構(gòu)的可行性，并且利用光聲方法研究了不同波長(zhǎng)處骨頭的光聲吸收譜和各種化學(xué)成分的光聲譜[4-6]。Lashkari課題組通過(guò)雙背向散射超聲和光聲雷達(dá)系統(tǒng)評(píng)估皮質(zhì)骨和小梁骨的結(jié)構(gòu)和密度變化[7-8]。Steinberg課題組使用了雙模態(tài)多光譜光聲系統(tǒng)來(lái)量化骨髓中的血脂比，該血脂比與長(zhǎng)骨中的分子變化有關(guān)[9]。本課題組前期研究了使用光聲時(shí)頻譜分析方法評(píng)估骨礦物質(zhì)密度、骨微結(jié)構(gòu)和化學(xué)成分的可行性[10-11]，利用多波長(zhǎng)光聲分析技術(shù)在數(shù)值模擬和離體實(shí)驗(yàn)中評(píng)估了動(dòng)物模型和人體骨骼標(biāo)本中的膠原含量[12]，還提出了一種用于定位和分割人體跟骨光聲信號(hào)空間信息的超聲引導(dǎo)的光聲骨質(zhì)評(píng)估方法，從而可以從其他信號(hào)中分離出小梁骨的信號(hào)[13]。

在上述文獻(xiàn)報(bào)道中，通常采用激光擴(kuò)束后覆蓋整個(gè)骨樣本以激發(fā)光聲信號(hào)的方式。骨小梁、骨髓或者微血管等分布式的光吸收分子激發(fā)出的光聲信號(hào)在固液兩相多孔結(jié)構(gòu)上形成復(fù)雜的多傳播路徑，最終檢測(cè)到的信號(hào)是分布式骨微結(jié)構(gòu)與分布式光聲信號(hào)相互作用的結(jié)果，即在時(shí)間域與空間域上多維信息混疊的結(jié)果，這給定量分析骨微結(jié)構(gòu)增加了難度和復(fù)雜度。

為了更準(zhǔn)確地量化評(píng)估骨微結(jié)構(gòu)信息，本文對(duì)PAT系統(tǒng)進(jìn)行改進(jìn)：將擴(kuò)束式光源改為弱聚焦式光源以分離局部光聲信號(hào)，在輻照光斑中心與試樣中心連線上設(shè)置一對(duì)遠(yuǎn)、近檢測(cè)點(diǎn)以獲取差分衰減頻譜(Differential Attenuation Spectrum,DAS)，并通過(guò)軸對(duì)稱環(huán)形掃描方式考察骨微結(jié)構(gòu)各向異性對(duì)聲速和DAS 的影響。本文通過(guò)數(shù)值仿真計(jì)算和驗(yàn)證了松質(zhì)骨的孔隙率與PADAS 的相關(guān)性；提取了相應(yīng)的衰減頻譜特征參數(shù)，探討了其進(jìn)行骨質(zhì)評(píng)估和診斷的潛力。

1 光聲信號(hào)的數(shù)值模擬仿真與分析

1.1 松質(zhì)骨樣本模型

為了研究松質(zhì)骨的孔隙率φ對(duì)光聲信號(hào)衰減頻譜的影響，從人跟骨CT 圖中提取的直徑為13.11 mm的圓形松質(zhì)骨樣本作為本文的仿真對(duì)象，如圖1中1號(hào)骨組織所示，白色代表骨小梁，主要由骨礦物質(zhì)構(gòu)成，骨小梁間隙的黑色部分代表骨髓，圓形骨樣本外側(cè)黑色部分為耦合介質(zhì)(水)。為模擬骨質(zhì)疏松時(shí)孔隙率的變化，對(duì)上述圖像進(jìn)行二值化處理，提高二值化函數(shù)閾值，大于閾值的部分賦值1，即為白色，反之賦值0，為黑色，從而獲得骨小梁逐漸變細(xì)，孔隙率逐漸變大的樣本圖像，如圖1中2號(hào)～5號(hào)所示。

圖1 不同孔隙率的松質(zhì)骨圖片F(xiàn)ig.1 Pictures of cancellous bone with different porosities

圖1中5塊松質(zhì)骨樣本的孔隙率如表1所示。

表1 不同松質(zhì)骨樣本的孔隙率Table 1 Porosities of different cancellous bone samples

1.2 光聲激勵(lì)-檢測(cè)系統(tǒng)的改進(jìn)

本文將傳統(tǒng)的PAT 系統(tǒng)改進(jìn)為偏心激勵(lì)-差分檢測(cè)模式，如圖2所示。圖2中，黑色代表骨小梁，白色代表骨髓，紅色圓形區(qū)域?yàn)楣饴曉矗疑L(zhǎng)方形為檢測(cè)點(diǎn)A 和B 所在位置。如圖1 所示，采用光斑直徑為2 mm 的弱聚焦脈沖光源輻照在圖1中的骨頭樣本上激發(fā)局部的光聲信號(hào)，光斑與圓形樣本外緣內(nèi)切，輻照區(qū)域?yàn)楦信d趣區(qū)域(Region of Interest,ROI)。以光斑中心與骨樣本中心連線為檢測(cè)基準(zhǔn)線，在樣本外側(cè)沿基準(zhǔn)線設(shè)置一對(duì)近端和遠(yuǎn)端光聲信號(hào)探測(cè)點(diǎn)A和B以獲取差分光聲信號(hào)，A和B 距骨樣本中心距離均為8 mm；骨樣本可圍繞其中心旋轉(zhuǎn)，即對(duì)骨樣本進(jìn)行環(huán)形掃描檢測(cè)不同方向的ROI，以研究骨微結(jié)構(gòu)的各向異性對(duì)聲速和衰減頻譜的影響。激光波長(zhǎng)為骨礦物質(zhì)的特征吸收波長(zhǎng)，并假設(shè)僅由骨礦物質(zhì)吸收，骨髓等其他物質(zhì)不吸收該波長(zhǎng)激光的能量。

圖2 偏心激勵(lì)-差分檢測(cè)模式的光聲系統(tǒng)示意圖Fig.2 Photoacoustic system diagram of the eccentric excitation-differential detection mode

1.3 松質(zhì)骨樣本的光聲數(shù)值仿真參數(shù)

為研究松質(zhì)骨孔隙率對(duì)光聲信號(hào)的影響，本文利用Matlab軟件(R2018a)的k-wave工具包模擬松質(zhì)骨樣本在光源激勵(lì)下產(chǎn)生和傳播的光聲信號(hào)[14]。仿真參數(shù)設(shè)置如下：骨小梁中的聲速為3 200 m·s-1，骨髓中的聲速為1 500 m·s-1，骨小梁的密度為1 850 kg·m-3，骨髓的密度為1 000 kg·m-3[14]，骨松質(zhì)的聲衰減系數(shù)為9.94 dB·(MHz·cm)-1[15]，骨松質(zhì)四周耦合介質(zhì)的參數(shù)均采用水的各項(xiàng)聲學(xué)參數(shù)；采樣頻率為500 MHz，數(shù)值模擬仿真時(shí)長(zhǎng)為15 μs，確保遠(yuǎn)端B檢測(cè)點(diǎn)也可以采集到ROI產(chǎn)生的光聲信號(hào)的完整直達(dá)波。

同時(shí)，為了研究骨松質(zhì)的光聲衰減頻譜是否存在各向異性，在仿真模型中，保持激光和檢測(cè)點(diǎn)的位置不變，使每塊圓形骨組織樣本繞自身中心旋轉(zhuǎn)，以30°為間隔獲取骨組織樣本的12組光聲信號(hào)。

1.4 快波聲速的估算及有效信號(hào)范圍的確定

圖3 給出了1 號(hào)骨組織樣本在某一個(gè)角度下，A、B檢測(cè)點(diǎn)接收到的時(shí)域信號(hào)。

圖3 在A、B檢測(cè)點(diǎn)的有效信號(hào)時(shí)長(zhǎng)示意圖Fig.3 Schematic diagram of the effective signal durations detected at the points A and B

由于光聲源在空間上有一定分布，光斑不同位置處的信號(hào)到達(dá)A檢測(cè)點(diǎn)有時(shí)間差，為了確定A檢測(cè)點(diǎn)有效信號(hào)的時(shí)間范圍，首先根據(jù)B檢測(cè)點(diǎn)接收到的直達(dá)波初始時(shí)刻及其到光聲源的距離，計(jì)算出5 塊骨組織樣本的不同方向上各12 組信號(hào)的聲速，每塊骨組織樣本的最大、最小和平均聲速如表2所示。由表2 可知，在不同方向上的聲速差別較小，波動(dòng)均小于1%。根據(jù)每塊樣本的聲速和光斑直徑可估算光聲直達(dá)信號(hào)到達(dá)A檢測(cè)點(diǎn)的時(shí)間范圍。

表2 不同骨組織樣本中的聲速Table 2 Sound speeds in different bone tissue samples

由表2中的數(shù)值模擬仿真結(jié)果可以看出，骨組織樣本的平均聲速為1 600～1 800 m·s-1，結(jié)合模擬仿真中所用光斑直徑為2 mm，由此可計(jì)算出在無(wú)衰減情況下，ROI 產(chǎn)生的光聲信號(hào)傳播時(shí)間小于1.25 μs。因此，選取A檢測(cè)點(diǎn)接收到的前1.25 μs內(nèi)的有效信號(hào)作為ROI產(chǎn)生的光聲直達(dá)信號(hào)，如圖3(a)中虛線框中選中的信號(hào)所示。同時(shí)，選取B檢測(cè)點(diǎn)接收信號(hào)的第一個(gè)波包作為ROI對(duì)應(yīng)的有效直達(dá)信號(hào)，如圖3(b)所示。

1.5 光聲差分衰減頻譜

為了求得每塊骨組織樣本的光聲衰減頻譜，本文采用差分法，分別計(jì)算A和B檢測(cè)點(diǎn)接收到的直達(dá)波的光聲功率譜密度，再利用差分法計(jì)算骨組織樣本的衰減頻譜。具體方法如下：(1) 選取A 檢測(cè)點(diǎn)接收到的前1.25 μs 內(nèi)的有效信號(hào)作為無(wú)衰減直達(dá)信號(hào)；(2) 選取B 檢測(cè)點(diǎn)接收到的第一個(gè)波包作為有衰減直達(dá)信號(hào)，相比A點(diǎn)的接收信號(hào)，在黏滯衰減系數(shù)不變的前提下，衰減主要來(lái)自于光聲傳播路徑上骨樣本微結(jié)構(gòu)引起的散射衰減；(3) 分別計(jì)算兩個(gè)信號(hào)的功率譜密度，并計(jì)算差分衰減頻譜。(4) 根據(jù)式(1)計(jì)算每塊骨頭樣本在12個(gè)角度下衰減頻譜衰減系數(shù)的平均值及其誤差：

其中：αDAS為樣本的衰減系數(shù)(dB·(Hz·mm)-1)；PB(f)為B 檢測(cè)點(diǎn)接收的光聲功率譜密度，PA(f)為A檢測(cè)點(diǎn)接收的光聲功率譜密度；D為B檢測(cè)點(diǎn)接收到的信號(hào)傳播距離與A檢測(cè)點(diǎn)接收到的信號(hào)傳播距離的差值，數(shù)值上等于骨頭的直徑減去光斑的直徑(如圖2中所示)。

2 仿真結(jié)果與分析

圖4 給出了5 塊骨組織樣本在12 個(gè)角度下A、B 檢測(cè)點(diǎn)接收到有效信號(hào)的平均光聲功率譜密度。從圖4中可以看出：(1) 骨組織樣本產(chǎn)生的PA源是寬頻帶信號(hào)；(2) 經(jīng)骨組織衰減后的信號(hào)能量主要集中在低頻段，高頻時(shí)信號(hào)衰減嚴(yán)重；(3) 骨質(zhì)越疏松，經(jīng)骨組織衰減后的信號(hào)保留的高頻成分越多。

圖4 A、B檢測(cè)點(diǎn)有效信號(hào)功率譜Fig.4 Effective signal power spectra measured at the points of A and B of different bone samples

根據(jù)式(1)計(jì)算的5塊骨頭樣本的差分衰減頻譜系數(shù)如圖5所示。藍(lán)實(shí)線為每個(gè)樣本在12個(gè)方向上差分衰減頻譜的平均值，并給出了各方向差分衰減頻譜的波動(dòng)范圍。由圖5可見，差分衰減頻譜呈現(xiàn)顯著的線性衰減段(linear attenuation band,LAB)和飽和衰減段(saturation attenuation band,SAB)。由于松質(zhì)骨樣本中的聲速大約為1 700 m·s-1，光斑直徑約為2 mm，骨小梁厚度方向激發(fā)的光聲信號(hào)的最低頻率約為0.85 MHz，因此本文將0.85 MHz 作為線性衰減段的初始頻率。同時(shí)，本文將每塊骨組織樣本B 檢測(cè)點(diǎn)的直達(dá)信號(hào)功率譜由最大值衰減40 dB對(duì)應(yīng)的頻率記為該樣本的飽和衰減頻率(Satu‐ration Attenuation Frequency,SAF)，即線性衰減段和飽和衰減段的臨界頻率。對(duì)線性衰減段進(jìn)行線性擬合，擬合結(jié)果如圖5(b)所示。

圖5 不同骨頭樣本差分衰減頻譜及其擬合圖Fig.5 Differential attenuation spectra and their fitting diagrams of different bone samples

通過(guò)對(duì)比圖5中不同骨頭樣品的光聲差分衰減頻譜可以分析出：(1) 當(dāng)孔隙率一定時(shí)，光聲信號(hào)的衰減系數(shù)隨頻率增大，直至達(dá)到飽和衰減點(diǎn)，最大衰減系數(shù)接近4.5 dB·(Hz·mm)-1；(2) 當(dāng)孔隙率φ增大時(shí)，線性衰減頻段的擬合直線斜率隨之減小，即線性衰減頻段變寬、飽和衰減頻率增大。

為了對(duì)上述定性結(jié)論進(jìn)行量化分析，進(jìn)一步提取了每塊骨頭樣本的線性衰減頻段的擬合直線斜率和飽和衰減頻率作為光聲差分衰減頻譜的特征參數(shù)，分析其與孔隙率的關(guān)系。分析結(jié)果如表3和圖6所示，擬合直線斜率和飽和衰減頻率與骨樣本的孔隙率均呈強(qiáng)線性相關(guān)。因此，有望根據(jù)這兩個(gè)特征參數(shù)進(jìn)行骨質(zhì)孔隙率的定量評(píng)估。

表3 不同骨組織樣本衰減特性Table 3 Attenuation characteristics of bone tissue samples

圖6 差分衰減頻譜擬合曲線中線性衰減段和飽和衰減段與骨樣本孔隙率關(guān)系的線性回歸圖Fig.6 The linear regression plots of the relationships of LAB and SAP in the fitting line of the differential attenuation spectra with bone porosity

3 結(jié) 論

本文針對(duì)PAT中松質(zhì)骨多孔結(jié)構(gòu)造成的骨質(zhì)信息提取難點(diǎn)，設(shè)計(jì)了PAT偏心激勵(lì)-差分檢測(cè)系統(tǒng)，以提取差分衰減頻譜，從而獲取骨微結(jié)構(gòu)信息。利用Matlab k-wave工具包進(jìn)行數(shù)值仿真計(jì)算和驗(yàn)證。不同孔隙率骨樣本的仿真結(jié)果表明：光聲差分衰減頻譜的線性衰減頻段擬合直線斜率和衰減飽和頻率這兩個(gè)特征參數(shù)均與孔隙率有強(qiáng)的相關(guān)性。因此，基于光聲差分衰減頻譜的分析方法有望用于骨質(zhì)定量評(píng)估和診斷。