基于光聲差分衰減頻譜的骨質評估仿真研究

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

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

0 引 言

生物醫學光聲是近年來發展起來的一種非侵入式和非電離式的新型生物組織檢測方法[1-2]。當不同波長的脈沖激光照射到生物組織時，組織中不同生物分子對光的吸收將產生超聲信號，這種由光激發產生的超聲信號被稱為光聲信號[3]。光聲技術在評估骨質疏松癥方面具有顯著優勢，它既能提供骨礦物質、骨髓、血紅蛋白、膠原蛋白等分子信息又能評估骨骼微結構信息和力學狀態，在疾病的診斷中顯示了巨大的潛力。

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

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

為了更準確地量化評估骨微結構信息，本文對PAT系統進行改進：將擴束式光源改為弱聚焦式光源以分離局部光聲信號，在輻照光斑中心與試樣中心連線上設置一對遠、近檢測點以獲取差分衰減頻譜(Differential Attenuation Spectrum,DAS)，并通過軸對稱環形掃描方式考察骨微結構各向異性對聲速和DAS 的影響。本文通過數值仿真計算和驗證了松質骨的孔隙率與PADAS 的相關性；提取了相應的衰減頻譜特征參數，探討了其進行骨質評估和診斷的潛力。

1 光聲信號的數值模擬仿真與分析

1.1 松質骨樣本模型

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

圖1 不同孔隙率的松質骨圖片Fig.1 Pictures of cancellous bone with different porosities

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

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

1.2 光聲激勵-檢測系統的改進

本文將傳統的PAT 系統改進為偏心激勵-差分檢測模式，如圖2所示。圖2中，黑色代表骨小梁，白色代表骨髓，紅色圓形區域為光聲源，灰色長方形為檢測點A 和B 所在位置。如圖1 所示，采用光斑直徑為2 mm 的弱聚焦脈沖光源輻照在圖1中的骨頭樣本上激發局部的光聲信號，光斑與圓形樣本外緣內切，輻照區域為感興趣區域(Region of Interest,ROI)。以光斑中心與骨樣本中心連線為檢測基準線，在樣本外側沿基準線設置一對近端和遠端光聲信號探測點A和B以獲取差分光聲信號，A和B 距骨樣本中心距離均為8 mm；骨樣本可圍繞其中心旋轉，即對骨樣本進行環形掃描檢測不同方向的ROI，以研究骨微結構的各向異性對聲速和衰減頻譜的影響。激光波長為骨礦物質的特征吸收波長，并假設僅由骨礦物質吸收，骨髓等其他物質不吸收該波長激光的能量。

圖2 偏心激勵-差分檢測模式的光聲系統示意圖Fig.2 Photoacoustic system diagram of the eccentric excitation-differential detection mode

1.3 松質骨樣本的光聲數值仿真參數

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

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

1.4 快波聲速的估算及有效信號范圍的確定

圖3 給出了1 號骨組織樣本在某一個角度下，A、B檢測點接收到的時域信號。

圖3 在A、B檢測點的有效信號時長示意圖Fig.3 Schematic diagram of the effective signal durations detected at the points A and B

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

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

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

1.5 光聲差分衰減頻譜

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

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

2 仿真結果與分析

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

圖4 A、B檢測點有效信號功率譜Fig.4 Effective signal power spectra measured at the points of A and B of different bone samples

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

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

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

為了對上述定性結論進行量化分析，進一步提取了每塊骨頭樣本的線性衰減頻段的擬合直線斜率和飽和衰減頻率作為光聲差分衰減頻譜的特征參數，分析其與孔隙率的關系。分析結果如表3和圖6所示，擬合直線斜率和飽和衰減頻率與骨樣本的孔隙率均呈強線性相關。因此，有望根據這兩個特征參數進行骨質孔隙率的定量評估。

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

圖6 差分衰減頻譜擬合曲線中線性衰減段和飽和衰減段與骨樣本孔隙率關系的線性回歸圖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 結 論

本文針對PAT中松質骨多孔結構造成的骨質信息提取難點，設計了PAT偏心激勵-差分檢測系統，以提取差分衰減頻譜，從而獲取骨微結構信息。利用Matlab k-wave工具包進行數值仿真計算和驗證。不同孔隙率骨樣本的仿真結果表明：光聲差分衰減頻譜的線性衰減頻段擬合直線斜率和衰減飽和頻率這兩個特征參數均與孔隙率有強的相關性。因此，基于光聲差分衰減頻譜的分析方法有望用于骨質定量評估和診斷。