基于多層聲速模型的合成孔徑超聲皮質骨成像*

李云清 江晨 李穎 徐峰 許凱亮? 他得安? 黎仲勛

1) (復旦大學電子工程系,上海 200433)

2) (阿爾伯塔大學放射與診斷圖像系,埃德蒙頓 T6G2B7)

1 引 言

隨著人口老齡化加劇,骨質疏松作為一種常見的代謝性骨病在日益加劇公共衛生負擔,相關骨質檢測與成像方法得到了學術界的廣泛關注,成為近年來的研究熱點[1].目前常用的骨質診斷技術有雙能X線吸收測定法(dual-energy X-ray absorptiometry,DXA)、定量CT法及定量超聲法等[2].在骨質評價中,DXA仍是骨質狀況評價的金標準.由于超聲評價方法具有無輻射和便攜低廉等特點,更適于臨床普查與病患監測,近年來得到較快發展[3].

研究表明,長骨皮質骨能夠支持超聲導波的傳播,基于軸向傳播法所發展起來的超聲導波長骨檢測技術,能提供長骨皮質骨的材料特性和內部結構等信息,已被用于長骨皮質骨評價及骨質疏松診斷[4,5].超聲導波可用于長骨中頻散和衰減的測量[6,7],導波的模式轉換規律可用于長骨骨折評估[8],其相速度的變化可用于長骨疲勞評估[9].導波信號的Randon變換可用于長骨皮質骨厚度估計[10].臨床測試表明,點醫療超聲(point-of-care ultrasound,PoCU)檢測長骨骨折的靈敏度可達到64.7%—100%[11].劉洋等[12]利用超聲非線性和時間反轉技術對長骨骨裂成像進行了仿真研究,并分析了裂紋深度對成像結果的影響.然而,由于骨與軟組織間顯著的聲速差異,固定聲速的傳統超聲波束形成用于皮質骨成像仍存在困難.Li等[13]結合已知聲速模型對皮質骨進行了圖像重建.Renaud等[14]采用射線追蹤(ray-tracing)法尋找皮質骨最佳成像結果,并對皮質骨參數進行了估計.

合成孔徑技術最早應用于雷達系統,20世紀70年代起應用于超聲成像[15].相比傳統定點聚焦,合成孔徑采用動態聚焦,可提升整體探測區域的成像質量.已有大量研究將合成孔徑超聲成像技術用于血管內成像[16]、無創頸動脈彈性成像[17]、肝臟病變成像等軟組織檢測[18].然而骨成像相關研究仍有待完善,近年來Renaud等[14]的研究結果表明了合成孔徑超聲骨成像的可行性.此外,考慮皮質骨與軟組織間顯著的聲速差異,建立正確可靠的聲速模型是皮質骨超聲成像的關鍵.針對多層介質成像建立聲速模型最早應用于地球物理學[19],近年來廣泛應用于無損檢測領域[20].相位遷移法(phase shift migration,PSM)是一種可用于多層聲速模型的頻域重建算法,最早用于反射地震學中地球內部結構成像[21].Olofsson[22]將PSM用于無損檢測領域多層介質成像.

皮質骨可簡化為軟組織-皮質骨-骨髓三層介質模型,利用壓縮感知方法提取皮質骨上下表面回波,在無需皮質骨厚度先驗信息的情況下建立聲速模型.本文采用合成孔徑方法測量超聲信號,通過壓縮感知計算延時參數進而建立聲速模型,最后在合成孔徑數據和聲速模型的基礎上,利用PSM算法對皮質骨在頻域進行圖像重建.通過時域有限差分(finite-difference time-domain,FDTD)仿真和離體骨板實驗驗證上述方法成像的可行性和準確性.

2 基本原理

方法流程如圖1所示,主要包括合成孔徑測量超聲信號,壓縮感知計算延時參數進而估計皮質骨厚度建立聲速模型,及PSM算法重建圖像.下文中將一一展開介紹.

圖1 方法流程圖Fig.1.Flow chart of the proposed method.

2.1 基于合成孔徑方法的超聲信號測量

傳統波束形成多采用定點聚焦,即發射和接收時只對一點進行聚焦.定點聚焦僅在焦點處有較高分辨率,而焦點以外區域分辨率較低.此外,定點聚焦導致成像的方位分辨率受限于換能器參數[23].通常采用大孔徑換能器或提高換能器工作頻率以得到更高的方位分辨率.然而過大孔徑換能器缺少實用性,且換能器工作頻率與最大可探測深度成反比,無法同時滿足方位分辨率高和探測范圍廣的要求.不同于定點聚焦,圖2(a)所示的合成孔徑超聲通過小孔徑換能器的水平移動依次發射接收數據,研究表明合成孔徑超聲可有效解決定點聚焦存在的分辨率問題[23].

有關成像方位分辨率的計算如下[23]:

定點聚焦

合成孔徑

(1)和(2)式中,βf為定點聚焦的半功率波束角,λ為工作波長,d為換能器孔徑直徑,ρa為方位分辨率,R為成像目標距離超聲換能器的斜距.(3)—(5)式中,βs為合成孔徑的半功率波束角,Lsm為最大綜合長度,ρs為合成孔徑方位分辨率,z0為成像目標距離超聲換能器的垂直距離.對比(2)和(5)式可知,合成孔徑換能器尺寸越小方位分辨率越高,且方位分辨率與距離無關,保持任意探測區域分辨率的一致性.

仿真與實驗均采用圖2(b)所示發射合成孔徑超聲對信號進行測量.線性陣列共有N個陣元,每個陣元依次發射,由全部陣元接收.發射合成孔徑相比于合成孔徑聚焦,可以得到更加完整的數據信息,實現發射和接收的全動態聚焦,提高成像質量和幀率[24].發射合成孔徑每次發射接收得到1組數據,通過N次循環可以得到完備集數據,用于后續圖像重建.

圖2 合成孔徑超聲原理 (a)合成孔徑聚焦; (b)發射合成孔徑Fig.2.Principle of synthetic aperture ultrasound: (a) Synthetic aperture focusing technique; (b) synthetic transmit aperture.

2.2 壓縮感知延時參數計算與聲速模型建立

壓縮感知廣泛應用于具有稀疏變換域的信號[25].假設原信號x∈RN為N維實信號,構造正交基ψ=[ψ1,ψ2...ψN]∈RN×N,則x在正交基ψ下的表示為[25]

(6)式的矩陣形式如下:

其中θ=[θ1,θ2,···,θN]T為x在ψ上的系數向量.假設系數向量θ中只有K個非零量,即x可由K個基的線性組合進行重構,則x是K稀疏的.

考慮皮質骨表面回波信號波包的時域稀疏性,采用壓縮感知計算延時參數,并建立聲速模型.回波信號可視作發射脈沖經過延時和幅度變化后的線性疊加,因此設計不同延時發射脈沖作為壓縮感知基底,延時間隔設置為采樣間隔.PSM算法要求數據是0延時形式,首先采用壓縮感知提取各個接收通道間的延時,將完備集數據調整為0延時形式,以便直接應用PSM算法.進而提取骨與軟組織交界處的兩次反射回波所對應基底,根據基底間延時和已知骨聲速估計骨板厚度,從而建立聲速模型.

2.3 多層介質PSM原理

PSM是基于爆炸反射模型(exploding reflector model)的一種頻域算法.在爆炸反射模型中,假設聲速為v,當發射陣元和接收陣元位置固定時,聲波從發射陣元至探測點的路徑及探測點至接收陣元的路徑一致,因此反射波可視為由反射表面的一系列點源在t=0時刻產生,此時聲速為

首先建立圖3(a)所示x-z坐標下的單層介質模型,根據波動方程有[22]

其中p(t,x,z) 表示 (x,z) 處t時刻的聲場.(8)式的解為[22]

其中,復數形式P表示幅度,?表示角頻率,kx和kz分別表示x方向和z方向的波數.根據(8)和(9)式可以得到[22]

考慮回波傳播方向為–z,解得kz為[22]

觀察可知(12)式為傅里葉反變換形式,相應傅里葉變換形式為[22]

圖3 PSM模型 (a)單層介質模型; (b)多層介質模型Fig.3.PSM model: (a) Single layer model; (b) multi-layer model.

根據爆炸反射模型,爆炸瞬間(t=0)聲場最為集中,聚焦程度最大,因此要求數據為0延時形式.

在單層介質模型基礎上設計多層介質模型.如圖3(b)所示,Z0為線陣放置位置,dl,vl,Zl分別表示第l(l=1,2,3)層介質的厚度、聲速和交界處位置,在多層介質模型中(11)式改寫為[22]

對于各層介質內 (?z

對于交界面處的聲場,

因此PSM可對多層介質模型任意深度z成像.

2.4 FDTD仿真

目前FDTD方法廣泛應用于電磁學、光學、聲學等領域[26].仿真中建立二維皮質骨模型,采用FDTD數值仿真方式模擬超聲波動.FDTD方法將計算區間離散化為眾多微分單元,在時間和空間域上進行差分化,進而得到超聲波動控制方程的數值解[27].超聲波動控制方程為

其中ρ為材料密度;w為某點處位移矢量;λ,μ分別為第一、第二拉梅常數.FDTD仿真中,皮質骨模型設置為三層介質模型,上層為3 mm厚軟組織,中間為3.4 mm厚皮質骨,下層為3.6 mm厚骨髓.軟組織及骨髓的密度設置為1000 kg/m3,第一拉梅常數設置為2.241 GPa,第二拉梅常數設置為0 GPa,聲速設置為1500 m/s[28].皮質骨密度設置為1850 kg/m3,第一拉梅常數設置為9.306 GPa,第二拉梅常數設置為3.127 GPa,聲速設置為2900 m/s[29].仿真采用128陣元線性陣列,線陣放置于軟組織表面,陣元間隔0.3 mm,中心頻率6.25 MHz,網格尺寸設置為0.02 mm×0.02 mm,迭代頻率250 MHz.為減少多次回波的干擾,在模型四周設置4 mm厚完全匹配層(perfectly matched layer).發射脈沖為高斯包絡調制的正弦波.采用發射合成孔徑方法測量信號.

2.5 實驗設置

實驗中選用3.4 mm厚牛脛骨骨板作為實驗材料.測量聲速為2900—3100 m/s,與仿真參數基本一致.為更好地模擬在體實驗,用3%瓊脂制成軟組織,覆蓋于骨板四周.實驗裝置見圖4,采用Verasonics系統(Vantage 128 or 256,Verasonics Inc,WA,USA)進行信號測量,線陣緊貼軟組織表面,探頭型號為L11-4v,包含128個陣元,陣元間距0.3 mm,中心頻率6.25 MHz,采樣頻率25 MHz.發射脈沖為高斯包絡調制的正弦波.為減少成像深度影響,實驗中采用線性時間增益補償.Verasonics系統控制陣元依次發射脈沖,每次發射所有128個陣元進行接收,數據經總線傳輸到計算機.

圖4 實驗裝置示意圖Fig.4.Experiment setup.

3 仿真及實驗結果

3.1 仿真結果

圖5為皮質骨FDTD仿真結果,其中AN為歸一化幅度.圖5(a)為仿真采用的發射脈沖經過不同延時組成的壓縮感知基底.圖5(b)為FDTD仿真的單通道接收信號,在4.6和7.0 μs附近可觀察到兩次明顯的反射回波,分別對應軟組織和骨板上表面的分界面以及骨板下表面和軟組織的分界面.因超聲存在衰減,第二回波幅度明顯小于第一回波.圖5(c)為第一個陣元發射時,前30個通道的接收信號,因回波信號到達各個陣元的路徑不同,通道間波形存在延時.壓縮感知計算各通道間接收延時參數,并對信號進行延時調整,結果如圖5(d)所示.與圖5(c)對比可得,調整后信號的第一回波延時已對齊,可進行后續圖像重建.

圖6為采用PSM算法對FDTD仿真數據進行圖像重建的結果.未建立聲速模型的重建結果如圖6(a)所示,雖可獲得清晰的上下骨板界面,但因聲速未經矯正,骨板厚度存在誤差,平均骨板厚度約為1.8 mm,相對誤差47.1%.經壓縮感知建立的多層聲速模型調整后的重建結果如圖6(b)所示,平均骨板厚度約為3.4 mm,相對誤差為0%,使用多層聲速模型可實現骨板厚度的正確成像.

為驗證壓縮感知建立聲速模型的準確性,采用FDTD方法對3—5 mm不同厚度的皮質骨進行多次仿真.壓縮感知估計的皮質骨厚度及誤差如表1所列.

壓縮感知估計皮質骨厚度的平均相對誤差為4.9%,方差為13.5%.

3.2 實驗結果

圖7為牛脛骨骨板的實驗結果.圖7(a)為實驗采用的發射脈沖經過不同延時組成的壓縮感知基底.圖7(b)為實驗的單通道接收信號,在16.3和18.7 μs附近可觀察到兩次明顯的反射回波,分別對應軟組織和骨板上表面的分界面以及骨板下表面和軟組織的分界面.因超聲存在衰減,第二回波的幅度明顯小于第一回波.圖7(c)為第一個陣元發射時,前30個通道的接收信號,因回波信號到達各個陣元的路徑不同,通道間波形存在延時.由壓縮感知計算各通道間接收延時參數,并對相應信號進行延時調整,結果如圖7(d)所示.與圖7(c)對比可得,調整后信號的第一回波延時已對齊,從而可進行后續圖像重建.

圖5 仿真結果 (a)壓縮感知前30個基底; (b)單通道接收信號; (c)單次發射全部通道接收信號; (d)壓縮感知調整的單次發射全部通道接收信號Fig.5.Simulated results: (a) The first 30 bases of compressed sensing; (b) received signal of single element; (c) received signals of all elements; (d) compressed sensing based temporally adjusted received signals of all elements.

圖6 仿真重建結果 (a)未建立聲速模型的仿真重建結果; (b)建立聲速模型的仿真重建結果Fig.6.Simulated reconstructed results: (a) Simulated reconstructed result without velocity model; (b) simulated reconstructed result with velocity model.

表1 皮質骨厚度估計及誤差Table 1.Estimation and relative error of cortical bone thickness.

圖7 實驗結果 (a)壓縮感知前30個基底; (b)單通道接收信號; (c)單次發射全部通道接收信號; (d)壓縮感知調整的單次發射全部通道接收信號Fig.7.Experiment results: (a) The first 30 bases of compressed sensing; (b) received signal of single element; (c) received signals of all elements; (d) compressed sensing based temporally adjusted received signals of all elements.

圖8為采用PSM算法對實驗數據進行圖像重建的結果,圖8(a)與圖8(b)分別對應未建立聲速模型的重建圖像以及基于多層聲速模型的重建圖像.如圖8(a)所示,雖可看到清晰的上下骨板界面,但由于沒有進行聲速矯正,骨板厚度存在誤差,平均骨板厚度約為1.8 mm,相對誤差47.1%,形態重建不準確.而多層聲速模型的應用實現了準確的板厚成像.對同一骨板進行重復實驗,估計得到平均骨板厚度約為3.5 mm,平均相對誤差為3.6%,方差為5.4%.可見聲速模型可實現較為準確的骨板成像.

圖8 實驗重建結果 (a)未建立聲速模型的實驗重建結果; (b)建立聲速模型的實驗重建結果Fig.8.Experiment reconstructed results: (a) Experiment reconstructed result without velocity model; (b) experiment reconstructed result with velocity model.

4 討論部分

由于皮質骨和軟組織間顯著的聲速差異,傳統超聲波束形成無法應用于骨成像,建立正確聲速模型是解決聲速差異問題的關鍵.目前仿真及實驗中聲速模型的建立需要事先測量真實骨板尺寸,聲速模型的實際應用受到限制.因此本文提出了基于壓縮感知建立聲速模型的方法,通過壓縮感知計算延時參數從而在未測量骨板尺寸的情況下建立較為準確的聲速模型.

圖9 壓縮感知與Hilbert變換比較 (a)原始信號及帶噪信號; (b)針對原始信號和帶噪信號的壓縮感知結果; (c)針對原始信號和帶噪信號的Hilbert變換結果Fig.9.Comparison of compressed sensing and Hilbert transform: (a) Origin signal and noisy signal; (b) result of compressed sensing; (c) result of Hilbert transform.

圖9比較了壓縮感知和Hilbert變換提取延時的結果.圖9(a)中黑線表示原始信號,原始信號由兩個不同延時的發射脈沖疊加組成,紅線為信噪比20 dB的帶噪信號.圖9(b)為壓縮感知提取結果,其中峰值對應發射脈沖起始點.圖9(c)為Hilbert變換提取結果,其中峰值對應發射脈沖包絡頂點.比較圖9(b)和圖9(c)可知,壓縮感知的提取結果受噪聲影響較小,延時提取準確,而Hilbert變換在有噪聲條件下提取的信號峰值相比無噪聲條件下發生了顯著的偏移.壓縮感知是一個求解最優解的過程,在計算中允許誤差項存在,因此壓縮感知法比Hilbert變換更適用于有噪聲條件下的延時估計.且壓縮感知提取結果對應脈沖起始點,可直接用于波形到達時間的計算.皮質骨仿真和實驗結果均驗證了壓縮感知提取延時的準確性,厚度估計的相對誤差分別為4.9%和3.6%.如圖6(b)和圖8(b)所示,建立聲速模型的仿真重建結果和實驗重建結果均可以看到清晰的皮質骨上下界面,且相比圖6(a)和圖8(a)未建立聲速模型的結果,皮質骨重建厚度準確.與仿真重建結果比較,實驗重建結果存在明顯偽影,這是實驗過程中存在噪聲干擾導致的.

在未知實際皮質骨尺寸的情況下,采用壓縮感知方法估計模型各層延時,結合已知聲速可計算模型各層厚度,正確建立聲速模型.目前本文方法只適用于多層規則介質模型.射線追蹤法可用于多層非規則介質的時域重建.Qin等[30]研究了多層非規則介質的頻域重建.后續研究將圍繞非規則骨板成像展開,在現有方法上進行改進,解決水平方向聲速變化問題.

5 結 論

本文基于合成孔徑法采集超聲信號,通過壓縮感知計算延時參數建立聲速模型,采用PSM算法進行圖像重建,實現皮質骨超聲成像.仿真和實驗結果表明,本文采用的壓縮感知方法可準確估計皮質骨厚度,從而建立多層聲速模型.重建圖像可觀察到清晰的皮質骨上下界面,且形態正確.本文方法相比傳統超聲波束形成,解決了固定聲速重建圖像的挑戰,對于皮質骨超聲成像的發展有一定借鑒意義.相比于Hilbert變換,壓縮感知方法更適于計算帶噪實驗信號的延時參數.本文方法目前適于規則骨板成像,考慮到骨的多種形態及復雜結構,之后的工作應當建立更為實際的骨模型,尋找更加通用的成像方法.同時探究在體實驗的可行性.