基于超聲導波反演參數定征皮質骨骨質狀況研究

盛帆，宋小軍,2，樊天地，張璠，李義方

(1.上海電力大學電子與信息工程學院，上海 201300；2.復旦大學生物醫學工程中心，上海 200433；3.國網安徽省電力有限公司超高壓分公司，安徽合肥 230000；4.復旦大學工程與應用技術研究院，上海 200433)

0 引言

隨著人口老齡化的加劇，與衰老密切相關的骨質疏松越發引起人們的關注。當發生骨質疏松時，骨質將會減少，骨孔隙度增大，骨折風險增加[1]。目前骨質疏松的主要檢查手段有雙能X射線吸收法(DXA)[2]、計算機斷層掃描技術(CT)[3]。與之相比，定量超聲法具有無輻射、穿透力強、費用低、便攜等優點[4-5]，在骨質疏松檢測中具有良好的應用前景。

軸向透射(AT)技術已被用于評估皮質骨[6]。AT技術可以從皮質骨采集到兩種超聲信號。一種是第一到達波信號(FAS)，另一種是超聲導波(UWG)[7]。FAS對應的超聲波沿皮質骨表面傳播，而UWG則在整個皮質骨層傳播，對骨微觀結構更為敏感[8-9]。

Patrick 等[10]在Lamb 波中提取了對皮質骨厚度較敏感的A0 模式，并成功反演出了皮質骨厚度。Daniel Pereira等[11]在低頻超聲下(<60 kHz)對5具尸體橈骨(年齡53～88歲)進行研究，成功反演出密度、厚度、外徑和形狀因子(SF)。Tran等[12]開發了一種基于網格搜索算法的反演程序獲得了單層骨板的皮質骨厚度和體波速度。Li等[13]設計一種卷積神經網絡算法，搭建皮質骨參數和超聲導波之間的映射關系，成功反演了牛皮質骨厚度、縱波速度和橫波速度。

以上研究成功反演出了皮質骨厚度等參數，但并未對皮質骨參數和骨質疏松相關性展開研究。本文利用有限元方法構建不同孔隙度的皮質骨模型來模擬不同程度的骨質疏松，并利用Floquet-Bloch理論搭建理論頻散曲線數據庫，通過反演算法獲得皮質骨厚度，橫波速度和縱波速度參數，以研究骨質疏松與皮質骨參數之間的定量關系。最后，進行離體牛脛骨實驗驗證了反演算法的合理性。

1 理論與方法

1.1 皮質骨中的蘭姆波

蘭姆波(Lamb)是一種常見的在自由邊界固定板(或層)中傳播的超聲導波，具有頻散性和多模態性[14-15]。當超聲波在自由板中傳播時，波在板的邊界發生多次反射與轉換，橫波縱波之間相互干擾，最終形成穩定的波包，即不同模式的蘭姆波。對于厚度為2d的自由板，不同的導波組合成的多模式頻散曲線可由Rayleigh-Lamb[16]方程得出，對稱模式時：

反對稱模式時：

式中：k是波數，k=ω/Vp；Vp是Lamb波的相速度；ω是角頻率；系數p和q由式(3)給出；VL和VT分別是縱波速度和橫波速度。波數和頻率的關系曲線即頻散曲線，可通過求解方程式(1)和式(2)獲得。

1.2 時域信號處理及頻散曲線提取

為了利用超聲導波反演皮質骨厚度、縱波速度和橫波速度，需要從采集到的超聲信號g(x,t)中提取仿真實驗和離體實驗的頻散數據。使用二維傅里葉變換(2D-FT)將采集到的時間t和空間x域的超聲信號g(x,t)轉換為頻率ω和波數k域。模態能量的分布H(ω,k)即相應的頻散軌跡，H(ω,k)的表達式為

二維傅里葉變換方法容易受到噪聲和頻率混疊的影響，特別是在離體實驗中，由于受到探頭尺寸、數量以及測量范圍的影響，致使離體數據波數域的分辨率相對較低。為了提高波數域的分辨率，本文用Burg 算法獲得更高分辨率的ω-k域頻散數據[17]。

1.3 理論頻散曲線與數據庫

Floquet-Bloch理論已應用在不同結構的超聲頻散特性分析中[18]，相較于傳統理論頻散曲線計算的譜方法，具有建模簡單，計算簡單，頻散曲線模式較多的優點[19]。二維無限長板可以近似看作沿x軸方向上周期為L的周期性結構。Floquet-Bloch理論實際在計算單元兩側應用周期性邊界條件模擬無限長波導，將無限長薄板簡化為微單元進行計算，以達到簡化建模和計算的目的，表達式為

式中：usrc和udst是源邊界和終點邊界的位移；rsrc和rdst是源邊界和終點邊界的坐標。對于微小的單元模型，在給定的波數k下，計算出對應的特征頻率f，得到k-f數據對，作為理論頻散曲線的原始數據。理論頻散曲線數據庫在參數預設范圍內，按照設定的增量步長，從預設范圍內的最小值逐步增大到最大值。理論頻散曲線數據庫規格如下：皮質骨厚度h的范圍為2～5 mm，增量步長為0.5 mm，縱波速度VL的范圍為3 300～4 200 m·s-1，增量步長為50 m·s-1，橫波速度VT范圍為1 500～2 500 m·s-1，增量步長為25 m·s-1。

1.4 反演算法

經過2D-FT 和Burg 算法處理后，實驗頻譜能量將保存在二維矩陣A中，具體表述為

式中：m和n分別是頻率和波數的范圍；aij是(i,j)點對應的能量值；理論頻散曲線在波數域和頻域可以表示為矩陣B：

如果理論頻散曲線經過點(i,j)，則bij=1,如果理論頻散曲線不經過點(i,j)，則bij=0。

非線性目標函數J的表達式為

圖1 參數反演流程圖Fig.1 Flow chart of parameter inversion

2 仿真與實驗

2.1 仿真實驗

本文以有限元法為基礎，通過仿真軟件COMSOL 求解偏微分方程來實現真實物理場景的仿真。選擇固體力學模塊和壓力聲學模塊中的時域求解器，建立二維仿真模型，如圖2所示，中間為帶孔皮質骨，長度為50 mm，為簡單起見，將孔隙設為均勻分布，孔隙度范圍設置為0～27%，間隔3%，共10組，用以模擬不同程度的骨質疏松。完美匹配層(Perfectly Matched Layer,PML)是一種人工吸收介質，可以實現對入射波的準完美吸收，以減少邊界反射，將它設置于模型的兩側。皮質骨參數設置如表1所示，與人體參數類似。

表1 皮質骨有限元仿真模型參數Table 1 Simulation model parameters for cortical bone

圖2 皮質骨有限元仿真模型Fig.2 Simulation model of cortical bone

超聲波激勵由發射探頭發出，入射角45°，經過楔塊傳入皮質骨中。激勵信號中心頻率過高，會導致超聲信號嚴重衰減，激勵信號中心頻率過低，會使頻散曲線模態較少，不利于反演，故本文采用中心頻率1 MHz 的3 周期高斯包絡正弦波作為激勵，如圖3所示，表達式為

圖3 激勵信號波形Fig.3 Waveform of excitation signal

式中：H為信號幅值；fc為中心頻率；σ為標準差，設為6×10-7；t0為中心位置，設為1.5 μs。

接收探頭共60 組，相鄰兩組間距為0.5 mm，與皮質骨表面垂直，為確保計算的穩定性、精度和計算效率，將采樣頻率設為10 MHz，共采集60組軸向時域信號。網格預定義設置為極細化，最大單元為524 μm，最小單元格為3 μm。

2.2 離體實驗

從市場購買6組牛脛骨，將表面軟組織刮除并進行脫脂處理，用游標卡尺測量牛皮質厚度。分別為2.78、2.82、3.58、3.70、3.83 和3.92 mm。縱波速度與橫波速度無法直接測出，為確保反演結果的準確性，后續會將離體實驗提取的頻散曲線軌跡與反演結果對應的理論頻散曲線進行對比。

實驗裝置如圖4所示，單個發射探頭和多組接收探頭陣列在試樣表面的軸向上對齊。使用中心頻率為1 MHz、直徑為15 mm的外部發射探頭。接收陣列是一個具有128組探頭(本實驗僅使用前64組)的多元件陣列。由任意波形發生器(美國)產生3 周期高斯包絡正弦信號作為激勵，經由功率放大器后輸入到發射換能器，接收陣列連接到可編程多通道超聲波實驗平臺，采樣頻率為10 MHz。發射換能器保持固定，接收陣列使用前64 個接收換能器，接收探頭間距為0.5 mm，共采集64組超聲數據。

圖4 離體牛脛骨實驗裝置圖Fig.4 Experimental setup diagram of the isolated bovine tibia

3 實驗結果及分析

3.1 仿真實驗結果分析

對于孔隙度為0的皮質骨模型，將接收到的60組超聲導波信號作歸一化處理，結果如圖5(a)所示。其中橫坐標為超聲信號在發射探頭和接收探頭之間傳導時間，縱軸依次為60 組接收探頭通道編號。對60 組時域信號先后分別通過2D-FT 和Burg算法處理，得到頻散數據，結果如圖5(b)所示。圖5中紅色區間能量高，能量沿著頻散軌跡向兩側遞減。將頻散數據與數據庫中的理論頻散曲線進行匹配反演，非線性目標函數的最小值Jmin所對應的皮質骨厚度、縱波速度、橫波速度分別是4 mm、3 650 m·s-1、1 850 m·s-1，其中，縱波速度的相對誤差為1.35%，厚度和橫波速度與仿真模型的設置參數完全相同，驗證了反演算法的可行性。

圖5 皮質骨仿真模型中的時域信號和頻域信號Fig.5 Time and frequency domain signals in a simulated cortical bone model

將其余9組不同孔隙度皮質骨模型的超聲時域數據分別經2D-FT和Burg算法得到實驗頻散曲線，再利用反演算法獲得皮質骨的厚度、縱波速度和橫波速度。不同孔隙度皮質骨模型的反演結果如表2所示。

表2 不同孔隙度皮質骨模型的反演結果Table 2 Inversion results of cortical bone models with different porosities

10組不同孔隙度皮質骨模型的實驗頻散能量軌跡(彩色)與其反演參數理論頻散曲線(黑色)的對比，前者與后者基本保持一致，如圖6所示。圖6(a)顯示了實驗頻散能量軌跡(彩色)與皮質骨反演參數理論頻散曲線(黑色)的對比，前者與后者相一致，進一步驗證了反演算法的準確性。

圖6 不同孔隙度仿真實驗頻散曲線與對應的理論頻散曲線Fig.6 Comparison of the experimental dispersion trajectories with the corresponding theoretical dispersion curves under different porosities

當老年人患有骨質疏松時，與健康骨骼相比，其骨孔隙度將變大，這意味著其皮質層將變得多孔(充滿組織液)，密度降低。縱波可以在固體、液體和氣體中傳播，但固體中的縱波速度高于其他兩種介質[20]。相比之下，橫波不能在液體和氣體中傳播，只能在固體介質中傳播。不同孔隙度下反演得到的縱波速度和橫波速度如圖7所示。可以看出，隨著皮質骨中孔隙度的增加，反演得到的橫波速度和縱波速度都隨之下降，說明皮質骨的橫波速度和縱波速度與骨質疏松具有較強的相關性。

圖7 孔隙度與反演的縱波和橫波速度的關系Fig.7 Relationships of the inversed longitudinal and transverse velocities with porosities

為了研究皮質骨的橫波速度和縱波速度與骨質疏松相關性的強弱。本文用QT、QL分別表示橫波速度和縱波速度的敏感度，計算公式為

式中：VT0和VL0分別是孔隙度為0時反演得到的橫波速度和縱波速度；VT27和VL27分別是孔隙度為27%時反演得到的橫波速度和縱波速度。

計算得到橫波速度敏感度為19.0%，縱波速度敏感度為5.5%，表明在增加相同孔隙度的情況下，橫速度衰減更為明顯，即橫波速度對于孔隙度的變化比縱波速度更敏感。其原因可能是骨質疏松癥發生時，皮質骨中被組織液填充的孔隙變大變多，縱波在孔隙中傳播速度減慢，而橫波不能在孔隙中傳播，這也表明橫波速度在骨質疏松臨床診斷中的應用潛力更大。

3.2 離體實驗結果分析

圖8 為離體牛脛骨的64 組接收超聲時域信號(皮質骨厚度為2.78 mm)，每組時域信號都進行歸一化處理。將此64 組時域信號經過2D-FT 和Burg算法處理后提取出頻散能量。將頻散能量數據與數據庫中的理論頻散曲線進行反演，得到Jmin所對應的皮質骨參數分別是皮質骨厚度為3 mm、縱波速度為3 950 m·s-1和橫波速度為1 800 m·s-1。

圖8 離體牛脛骨實驗得到的64組超聲導波信號Fig.8 The 64-groups of ultrasonic guided wave signals obtained from the bovine cortical bone ex-vivo experiment

對3 組離體牛脛骨實驗數據的處理結果如表3所示。由表3可知，反演得到的皮質骨厚度與其實際測量值的誤差分別為7.9%，6.0%，2.2%，5.7%，4.3%，2.0%，平均相對誤差4.0%。

表3 6組牛脛骨離體實驗反演參數Table 3 Inversion parameters for the six sets of bovine cortical bone ex-vivo experiments

實驗頻散曲線(彩色)與反演參數對應的理論頻散曲線(黑色)如圖9所示。前者與后者吻合較好。以上結果說明本文的反演算法可以較準確地獲得與骨質疏松相關的厚度參數，驗證了本文算法在真實皮質骨參數反演上的可行性和準確性。

圖9 牛脛骨離體實驗得到的頻散能量軌跡與反演參數對應的理論頻散曲線對比Fig.9 Comparison of the dispersion energy trajectories obtained from the bovine cortical bone ex-vivo experiment with the corresponding theoretical dispersion curves

4 結論

針對皮質骨的骨質狀況，本文提出了一種基于超聲導波的反演算法，構建了不同孔隙度的皮質骨模型，利用Floquet-Bloch理論建立了單層皮質骨理論頻散曲線數據庫。仿真實驗結果顯示，隨著皮質骨中孔隙度的增加，反演得到的橫波速度和縱波速度都隨之下降，說明皮質骨的橫波速度和縱波速度與骨質疏松具有較強的相關性，且橫波速度相較縱波速度更敏感，臨床檢測潛力更大。離體牛脛骨實驗驗證了本文算法在真實皮質骨參數反演上的可行性和準確性。在今后的研究中，我們將與醫院合作，采集骨質疏松患者的超聲導波信號，從而進一步分析皮質骨參數與人體骨質疏松的關系。