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

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

0 引言

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

軸向透射(AT)技術(shù)已被用于評估皮質(zhì)骨[6]。AT技術(shù)可以從皮質(zhì)骨采集到兩種超聲信號。一種是第一到達(dá)波信號(FAS)，另一種是超聲導(dǎo)波(UWG)[7]。FAS對應(yīng)的超聲波沿皮質(zhì)骨表面?zhèn)鞑?，而UWG則在整個(gè)皮質(zhì)骨層傳播，對骨微觀結(jié)構(gòu)更為敏感[8-9]。

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

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

1 理論與方法

1.1 皮質(zhì)骨中的蘭姆波

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

反對稱模式時(shí)：

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

1.2 時(shí)域信號處理及頻散曲線提取

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

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

1.3 理論頻散曲線與數(shù)據(jù)庫

Floquet-Bloch理論已應(yīng)用在不同結(jié)構(gòu)的超聲頻散特性分析中[18]，相較于傳統(tǒng)理論頻散曲線計(jì)算的譜方法，具有建模簡單，計(jì)算簡單，頻散曲線模式較多的優(yōu)點(diǎn)[19]。二維無限長板可以近似看作沿x軸方向上周期為L的周期性結(jié)構(gòu)。Floquet-Bloch理論實(shí)際在計(jì)算單元兩側(cè)應(yīng)用周期性邊界條件模擬無限長波導(dǎo)，將無限長薄板簡化為微單元進(jìn)行計(jì)算，以達(dá)到簡化建模和計(jì)算的目的，表達(dá)式為

式中：usrc和udst是源邊界和終點(diǎn)邊界的位移；rsrc和rdst是源邊界和終點(diǎn)邊界的坐標(biāo)。對于微小的單元模型，在給定的波數(shù)k下，計(jì)算出對應(yīng)的特征頻率f，得到k-f數(shù)據(jù)對，作為理論頻散曲線的原始數(shù)據(jù)。理論頻散曲線數(shù)據(jù)庫在參數(shù)預(yù)設(shè)范圍內(nèi)，按照設(shè)定的增量步長，從預(yù)設(shè)范圍內(nèi)的最小值逐步增大到最大值。理論頻散曲線數(shù)據(jù)庫規(guī)格如下：皮質(zhì)骨厚度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 反演算法

經(jīng)過2D-FT 和Burg 算法處理后，實(shí)驗(yàn)頻譜能量將保存在二維矩陣A中，具體表述為

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

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

非線性目標(biāo)函數(shù)J的表達(dá)式為

圖1 參數(shù)反演流程圖Fig.1 Flow chart of parameter inversion

2 仿真與實(shí)驗(yàn)

2.1 仿真實(shí)驗(yàn)

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

表1 皮質(zhì)骨有限元仿真模型參數(shù)Table 1 Simulation model parameters for cortical bone

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

超聲波激勵(lì)由發(fā)射探頭發(fā)出，入射角45°，經(jīng)過楔塊傳入皮質(zhì)骨中。激勵(lì)信號中心頻率過高，會導(dǎo)致超聲信號嚴(yán)重衰減，激勵(lì)信號中心頻率過低，會使頻散曲線模態(tài)較少，不利于反演，故本文采用中心頻率1 MHz 的3 周期高斯包絡(luò)正弦波作為激勵(lì)，如圖3所示，表達(dá)式為

圖3 激勵(lì)信號波形Fig.3 Waveform of excitation signal

式中：H為信號幅值；fc為中心頻率；σ為標(biāo)準(zhǔn)差，設(shè)為6×10-7；t0為中心位置，設(shè)為1.5 μs。

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

2.2 離體實(shí)驗(yàn)

從市場購買6組牛脛骨，將表面軟組織刮除并進(jìn)行脫脂處理，用游標(biāo)卡尺測量牛皮質(zhì)厚度。分別為2.78、2.82、3.58、3.70、3.83 和3.92 mm?？v波速度與橫波速度無法直接測出，為確保反演結(jié)果的準(zhǔn)確性，后續(xù)會將離體實(shí)驗(yàn)提取的頻散曲線軌跡與反演結(jié)果對應(yīng)的理論頻散曲線進(jìn)行對比。

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

圖4 離體牛脛骨實(shí)驗(yàn)裝置圖Fig.4 Experimental setup diagram of the isolated bovine tibia

3 實(shí)驗(yàn)結(jié)果及分析

3.1 仿真實(shí)驗(yàn)結(jié)果分析

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

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

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

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

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

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

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

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

為了研究皮質(zhì)骨的橫波速度和縱波速度與骨質(zhì)疏松相關(guān)性的強(qiáng)弱。本文用QT、QL分別表示橫波速度和縱波速度的敏感度，計(jì)算公式為

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

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

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

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

圖8 離體牛脛骨實(shí)驗(yàn)得到的64組超聲導(dǎo)波信號Fig.8 The 64-groups of ultrasonic guided wave signals obtained from the bovine cortical bone ex-vivo experiment

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

表3 6組牛脛骨離體實(shí)驗(yàn)反演參數(shù)Table 3 Inversion parameters for the six sets of bovine cortical bone ex-vivo experiments

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

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

4 結(jié)論

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