李云清 江晨 李穎 徐峰 許凱亮? 他得安? 黎仲勛
1) (復(fù)旦大學(xué)電子工程系,上海 200433)
2) (阿爾伯塔大學(xué)放射與診斷圖像系,埃德蒙頓 T6G2B7)
隨著人口老齡化加劇,骨質(zhì)疏松作為一種常見的代謝性骨病在日益加劇公共衛(wèi)生負(fù)擔(dān),相關(guān)骨質(zhì)檢測與成像方法得到了學(xué)術(shù)界的廣泛關(guān)注,成為近年來的研究熱點(diǎn)[1].目前常用的骨質(zhì)診斷技術(shù)有雙能X線吸收測定法(dual-energy X-ray absorptiometry,DXA)、定量CT法及定量超聲法等[2].在骨質(zhì)評價(jià)中,DXA仍是骨質(zhì)狀況評價(jià)的金標(biāo)準(zhǔn).由于超聲評價(jià)方法具有無輻射和便攜低廉等特點(diǎn),更適于臨床普查與病患監(jiān)測,近年來得到較快發(fā)展[3].
研究表明,長骨皮質(zhì)骨能夠支持超聲導(dǎo)波的傳播,基于軸向傳播法所發(fā)展起來的超聲導(dǎo)波長骨檢測技術(shù),能提供長骨皮質(zhì)骨的材料特性和內(nèi)部結(jié)構(gòu)等信息,已被用于長骨皮質(zhì)骨評價(jià)及骨質(zhì)疏松診斷[4,5].超聲導(dǎo)波可用于長骨中頻散和衰減的測量[6,7],導(dǎo)波的模式轉(zhuǎn)換規(guī)律可用于長骨骨折評估[8],其相速度的變化可用于長骨疲勞評估[9].導(dǎo)波信號的Randon變換可用于長骨皮質(zhì)骨厚度估計(jì)[10].臨床測試表明,點(diǎn)醫(yī)療超聲(point-of-care ultrasound,PoCU)檢測長骨骨折的靈敏度可達(dá)到64.7%—100%[11].劉洋等[12]利用超聲非線性和時(shí)間反轉(zhuǎn)技術(shù)對長骨骨裂成像進(jìn)行了仿真研究,并分析了裂紋深度對成像結(jié)果的影響.然而,由于骨與軟組織間顯著的聲速差異,固定聲速的傳統(tǒng)超聲波束形成用于皮質(zhì)骨成像仍存在困難.Li等[13]結(jié)合已知聲速模型對皮質(zhì)骨進(jìn)行了圖像重建.Renaud等[14]采用射線追蹤(ray-tracing)法尋找皮質(zhì)骨最佳成像結(jié)果,并對皮質(zhì)骨參數(shù)進(jìn)行了估計(jì).
合成孔徑技術(shù)最早應(yīng)用于雷達(dá)系統(tǒng),20世紀(jì)70年代起應(yīng)用于超聲成像[15].相比傳統(tǒng)定點(diǎn)聚焦,合成孔徑采用動(dòng)態(tài)聚焦,可提升整體探測區(qū)域的成像質(zhì)量.已有大量研究將合成孔徑超聲成像技術(shù)用于血管內(nèi)成像[16]、無創(chuàng)頸動(dòng)脈彈性成像[17]、肝臟病變成像等軟組織檢測[18].然而骨成像相關(guān)研究仍有待完善,近年來Renaud等[14]的研究結(jié)果表明了合成孔徑超聲骨成像的可行性.此外,考慮皮質(zhì)骨與軟組織間顯著的聲速差異,建立正確可靠的聲速模型是皮質(zhì)骨超聲成像的關(guān)鍵.針對多層介質(zhì)成像建立聲速模型最早應(yīng)用于地球物理學(xué)[19],近年來廣泛應(yīng)用于無損檢測領(lǐng)域[20].相位遷移法(phase shift migration,PSM)是一種可用于多層聲速模型的頻域重建算法,最早用于反射地震學(xué)中地球內(nèi)部結(jié)構(gòu)成像[21].Olofsson[22]將PSM用于無損檢測領(lǐng)域多層介質(zhì)成像.
皮質(zhì)骨可簡化為軟組織-皮質(zhì)骨-骨髓三層介質(zhì)模型,利用壓縮感知方法提取皮質(zhì)骨上下表面回波,在無需皮質(zhì)骨厚度先驗(yàn)信息的情況下建立聲速模型.本文采用合成孔徑方法測量超聲信號,通過壓縮感知計(jì)算延時(shí)參數(shù)進(jìn)而建立聲速模型,最后在合成孔徑數(shù)據(jù)和聲速模型的基礎(chǔ)上,利用PSM算法對皮質(zhì)骨在頻域進(jìn)行圖像重建.通過時(shí)域有限差分(finite-difference time-domain,FDTD)仿真和離體骨板實(shí)驗(yàn)驗(yàn)證上述方法成像的可行性和準(zhǔn)確性.
方法流程如圖1所示,主要包括合成孔徑測量超聲信號,壓縮感知計(jì)算延時(shí)參數(shù)進(jìn)而估計(jì)皮質(zhì)骨厚度建立聲速模型,及PSM算法重建圖像.下文中將一一展開介紹.
圖1 方法流程圖Fig.1.Flow chart of the proposed method.
傳統(tǒng)波束形成多采用定點(diǎn)聚焦,即發(fā)射和接收時(shí)只對一點(diǎn)進(jìn)行聚焦.定點(diǎn)聚焦僅在焦點(diǎn)處有較高分辨率,而焦點(diǎn)以外區(qū)域分辨率較低.此外,定點(diǎn)聚焦導(dǎo)致成像的方位分辨率受限于換能器參數(shù)[23].通常采用大孔徑換能器或提高換能器工作頻率以得到更高的方位分辨率.然而過大孔徑換能器缺少實(shí)用性,且換能器工作頻率與最大可探測深度成反比,無法同時(shí)滿足方位分辨率高和探測范圍廣的要求.不同于定點(diǎn)聚焦,圖2(a)所示的合成孔徑超聲通過小孔徑換能器的水平移動(dòng)依次發(fā)射接收數(shù)據(jù),研究表明合成孔徑超聲可有效解決定點(diǎn)聚焦存在的分辨率問題[23].
有關(guān)成像方位分辨率的計(jì)算如下[23]:
定點(diǎn)聚焦
合成孔徑
(1)和(2)式中,βf為定點(diǎn)聚焦的半功率波束角,λ為工作波長,d為換能器孔徑直徑,ρa(bǔ)為方位分辨率,R為成像目標(biāo)距離超聲換能器的斜距.(3)—(5)式中,βs為合成孔徑的半功率波束角,Lsm為最大綜合長度,ρs為合成孔徑方位分辨率,z0為成像目標(biāo)距離超聲換能器的垂直距離.對比(2)和(5)式可知,合成孔徑換能器尺寸越小方位分辨率越高,且方位分辨率與距離無關(guān),保持任意探測區(qū)域分辨率的一致性.
仿真與實(shí)驗(yàn)均采用圖2(b)所示發(fā)射合成孔徑超聲對信號進(jìn)行測量.線性陣列共有N個(gè)陣元,每個(gè)陣元依次發(fā)射,由全部陣元接收.發(fā)射合成孔徑相比于合成孔徑聚焦,可以得到更加完整的數(shù)據(jù)信息,實(shí)現(xiàn)發(fā)射和接收的全動(dòng)態(tài)聚焦,提高成像質(zhì)量和幀率[24].發(fā)射合成孔徑每次發(fā)射接收得到1組數(shù)據(jù),通過N次循環(huán)可以得到完備集數(shù)據(jù),用于后續(xù)圖像重建.
圖2 合成孔徑超聲原理 (a)合成孔徑聚焦; (b)發(fā)射合成孔徑Fig.2.Principle of synthetic aperture ultrasound: (a) Synthetic aperture focusing technique; (b) synthetic transmit aperture.
壓縮感知廣泛應(yīng)用于具有稀疏變換域的信號[25].假設(shè)原信號x∈RN為N維實(shí)信號,構(gòu)造正交基ψ=[ψ1,ψ2...ψN]∈RN×N,則x在正交基ψ下的表示為[25]
(6)式的矩陣形式如下:
其中θ=[θ1,θ2,···,θN]T為x在ψ上的系數(shù)向量.假設(shè)系數(shù)向量θ中只有K個(gè)非零量,即x可由K個(gè)基的線性組合進(jìn)行重構(gòu),則x是K稀疏的.
考慮皮質(zhì)骨表面回波信號波包的時(shí)域稀疏性,采用壓縮感知計(jì)算延時(shí)參數(shù),并建立聲速模型.回波信號可視作發(fā)射脈沖經(jīng)過延時(shí)和幅度變化后的線性疊加,因此設(shè)計(jì)不同延時(shí)發(fā)射脈沖作為壓縮感知基底,延時(shí)間隔設(shè)置為采樣間隔.PSM算法要求數(shù)據(jù)是0延時(shí)形式,首先采用壓縮感知提取各個(gè)接收通道間的延時(shí),將完備集數(shù)據(jù)調(diào)整為0延時(shí)形式,以便直接應(yīng)用PSM算法.進(jìn)而提取骨與軟組織交界處的兩次反射回波所對應(yīng)基底,根據(jù)基底間延時(shí)和已知骨聲速估計(jì)骨板厚度,從而建立聲速模型.
PSM是基于爆炸反射模型(exploding reflector model)的一種頻域算法.在爆炸反射模型中,假設(shè)聲速為v,當(dāng)發(fā)射陣元和接收陣元位置固定時(shí),聲波從發(fā)射陣元至探測點(diǎn)的路徑及探測點(diǎn)至接收陣元的路徑一致,因此反射波可視為由反射表面的一系列點(diǎn)源在t=0時(shí)刻產(chǎn)生,此時(shí)聲速為
首先建立圖3(a)所示x-z坐標(biāo)下的單層介質(zhì)模型,根據(jù)波動(dòng)方程有[22]
其中p(t,x,z) 表示 (x,z) 處t時(shí)刻的聲場.(8)式的解為[22]
其中,復(fù)數(shù)形式P表示幅度,?表示角頻率,kx和kz分別表示x方向和z方向的波數(shù).根據(jù)(8)和(9)式可以得到[22]
考慮回波傳播方向?yàn)楱Cz,解得kz為[22]
觀察可知(12)式為傅里葉反變換形式,相應(yīng)傅里葉變換形式為[22]
圖3 PSM模型 (a)單層介質(zhì)模型; (b)多層介質(zhì)模型Fig.3.PSM model: (a) Single layer model; (b) multi-layer model.
根據(jù)爆炸反射模型,爆炸瞬間(t=0)聲場最為集中,聚焦程度最大,因此要求數(shù)據(jù)為0延時(shí)形式.
在單層介質(zhì)模型基礎(chǔ)上設(shè)計(jì)多層介質(zhì)模型.如圖3(b)所示,Z0為線陣放置位置,dl,vl,Zl分別表示第l(l=1,2,3)層介質(zhì)的厚度、聲速和交界處位置,在多層介質(zhì)模型中(11)式改寫為[22]
對于各層介質(zhì)內(nèi) (?z
對于交界面處的聲場,
因此PSM可對多層介質(zhì)模型任意深度z成像.
目前FDTD方法廣泛應(yīng)用于電磁學(xué)、光學(xué)、聲學(xué)等領(lǐng)域[26].仿真中建立二維皮質(zhì)骨模型,采用FDTD數(shù)值仿真方式模擬超聲波動(dòng).FDTD方法將計(jì)算區(qū)間離散化為眾多微分單元,在時(shí)間和空間域上進(jìn)行差分化,進(jìn)而得到超聲波動(dòng)控制方程的數(shù)值解[27].超聲波動(dòng)控制方程為
其中ρ為材料密度;w為某點(diǎn)處位移矢量;λ,μ分別為第一、第二拉梅常數(shù).FDTD仿真中,皮質(zhì)骨模型設(shè)置為三層介質(zhì)模型,上層為3 mm厚軟組織,中間為3.4 mm厚皮質(zhì)骨,下層為3.6 mm厚骨髓.軟組織及骨髓的密度設(shè)置為1000 kg/m3,第一拉梅常數(shù)設(shè)置為2.241 GPa,第二拉梅常數(shù)設(shè)置為0 GPa,聲速設(shè)置為1500 m/s[28].皮質(zhì)骨密度設(shè)置為1850 kg/m3,第一拉梅常數(shù)設(shè)置為9.306 GPa,第二拉梅常數(shù)設(shè)置為3.127 GPa,聲速設(shè)置為2900 m/s[29].仿真采用128陣元線性陣列,線陣放置于軟組織表面,陣元間隔0.3 mm,中心頻率6.25 MHz,網(wǎng)格尺寸設(shè)置為0.02 mm×0.02 mm,迭代頻率250 MHz.為減少多次回波的干擾,在模型四周設(shè)置4 mm厚完全匹配層(perfectly matched layer).發(fā)射脈沖為高斯包絡(luò)調(diào)制的正弦波.采用發(fā)射合成孔徑方法測量信號.
實(shí)驗(yàn)中選用3.4 mm厚牛脛骨骨板作為實(shí)驗(yàn)材料.測量聲速為2900—3100 m/s,與仿真參數(shù)基本一致.為更好地模擬在體實(shí)驗(yàn),用3%瓊脂制成軟組織,覆蓋于骨板四周.實(shí)驗(yàn)裝置見圖4,采用Verasonics系統(tǒng)(Vantage 128 or 256,Verasonics Inc,WA,USA)進(jìn)行信號測量,線陣緊貼軟組織表面,探頭型號為L11-4v,包含128個(gè)陣元,陣元間距0.3 mm,中心頻率6.25 MHz,采樣頻率25 MHz.發(fā)射脈沖為高斯包絡(luò)調(diào)制的正弦波.為減少成像深度影響,實(shí)驗(yàn)中采用線性時(shí)間增益補(bǔ)償.Verasonics系統(tǒng)控制陣元依次發(fā)射脈沖,每次發(fā)射所有128個(gè)陣元進(jìn)行接收,數(shù)據(jù)經(jīng)總線傳輸?shù)接?jì)算機(jī).
圖4 實(shí)驗(yàn)裝置示意圖Fig.4.Experiment setup.
圖5為皮質(zhì)骨FDTD仿真結(jié)果,其中AN為歸一化幅度.圖5(a)為仿真采用的發(fā)射脈沖經(jīng)過不同延時(shí)組成的壓縮感知基底.圖5(b)為FDTD仿真的單通道接收信號,在4.6和7.0 μs附近可觀察到兩次明顯的反射回波,分別對應(yīng)軟組織和骨板上表面的分界面以及骨板下表面和軟組織的分界面.因超聲存在衰減,第二回波幅度明顯小于第一回波.圖5(c)為第一個(gè)陣元發(fā)射時(shí),前30個(gè)通道的接收信號,因回波信號到達(dá)各個(gè)陣元的路徑不同,通道間波形存在延時(shí).壓縮感知計(jì)算各通道間接收延時(shí)參數(shù),并對信號進(jìn)行延時(shí)調(diào)整,結(jié)果如圖5(d)所示.與圖5(c)對比可得,調(diào)整后信號的第一回波延時(shí)已對齊,可進(jìn)行后續(xù)圖像重建.
圖6為采用PSM算法對FDTD仿真數(shù)據(jù)進(jìn)行圖像重建的結(jié)果.未建立聲速模型的重建結(jié)果如圖6(a)所示,雖可獲得清晰的上下骨板界面,但因聲速未經(jīng)矯正,骨板厚度存在誤差,平均骨板厚度約為1.8 mm,相對誤差47.1%.經(jīng)壓縮感知建立的多層聲速模型調(diào)整后的重建結(jié)果如圖6(b)所示,平均骨板厚度約為3.4 mm,相對誤差為0%,使用多層聲速模型可實(shí)現(xiàn)骨板厚度的正確成像.
為驗(yàn)證壓縮感知建立聲速模型的準(zhǔn)確性,采用FDTD方法對3—5 mm不同厚度的皮質(zhì)骨進(jìn)行多次仿真.壓縮感知估計(jì)的皮質(zhì)骨厚度及誤差如表1所列.
壓縮感知估計(jì)皮質(zhì)骨厚度的平均相對誤差為4.9%,方差為13.5%.
圖7為牛脛骨骨板的實(shí)驗(yàn)結(jié)果.圖7(a)為實(shí)驗(yàn)采用的發(fā)射脈沖經(jīng)過不同延時(shí)組成的壓縮感知基底.圖7(b)為實(shí)驗(yàn)的單通道接收信號,在16.3和18.7 μs附近可觀察到兩次明顯的反射回波,分別對應(yīng)軟組織和骨板上表面的分界面以及骨板下表面和軟組織的分界面.因超聲存在衰減,第二回波的幅度明顯小于第一回波.圖7(c)為第一個(gè)陣元發(fā)射時(shí),前30個(gè)通道的接收信號,因回波信號到達(dá)各個(gè)陣元的路徑不同,通道間波形存在延時(shí).由壓縮感知計(jì)算各通道間接收延時(shí)參數(shù),并對相應(yīng)信號進(jìn)行延時(shí)調(diào)整,結(jié)果如圖7(d)所示.與圖7(c)對比可得,調(diào)整后信號的第一回波延時(shí)已對齊,從而可進(jìn)行后續(xù)圖像重建.
圖5 仿真結(jié)果 (a)壓縮感知前30個(gè)基底; (b)單通道接收信號; (c)單次發(fā)射全部通道接收信號; (d)壓縮感知調(diào)整的單次發(fā)射全部通道接收信號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 仿真重建結(jié)果 (a)未建立聲速模型的仿真重建結(jié)果; (b)建立聲速模型的仿真重建結(jié)果Fig.6.Simulated reconstructed results: (a) Simulated reconstructed result without velocity model; (b) simulated reconstructed result with velocity model.
表1 皮質(zhì)骨厚度估計(jì)及誤差Table 1.Estimation and relative error of cortical bone thickness.
圖7 實(shí)驗(yàn)結(jié)果 (a)壓縮感知前30個(gè)基底; (b)單通道接收信號; (c)單次發(fā)射全部通道接收信號; (d)壓縮感知調(diào)整的單次發(fā)射全部通道接收信號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算法對實(shí)驗(yàn)數(shù)據(jù)進(jìn)行圖像重建的結(jié)果,圖8(a)與圖8(b)分別對應(yīng)未建立聲速模型的重建圖像以及基于多層聲速模型的重建圖像.如圖8(a)所示,雖可看到清晰的上下骨板界面,但由于沒有進(jìn)行聲速矯正,骨板厚度存在誤差,平均骨板厚度約為1.8 mm,相對誤差47.1%,形態(tài)重建不準(zhǔn)確.而多層聲速模型的應(yīng)用實(shí)現(xiàn)了準(zhǔn)確的板厚成像.對同一骨板進(jìn)行重復(fù)實(shí)驗(yàn),估計(jì)得到平均骨板厚度約為3.5 mm,平均相對誤差為3.6%,方差為5.4%.可見聲速模型可實(shí)現(xiàn)較為準(zhǔn)確的骨板成像.
圖8 實(shí)驗(yàn)重建結(jié)果 (a)未建立聲速模型的實(shí)驗(yàn)重建結(jié)果; (b)建立聲速模型的實(shí)驗(yàn)重建結(jié)果Fig.8.Experiment reconstructed results: (a) Experiment reconstructed result without velocity model; (b) experiment reconstructed result with velocity model.
由于皮質(zhì)骨和軟組織間顯著的聲速差異,傳統(tǒng)超聲波束形成無法應(yīng)用于骨成像,建立正確聲速模型是解決聲速差異問題的關(guān)鍵.目前仿真及實(shí)驗(yàn)中聲速模型的建立需要事先測量真實(shí)骨板尺寸,聲速模型的實(shí)際應(yīng)用受到限制.因此本文提出了基于壓縮感知建立聲速模型的方法,通過壓縮感知計(jì)算延時(shí)參數(shù)從而在未測量骨板尺寸的情況下建立較為準(zhǔn)確的聲速模型.
圖9 壓縮感知與Hilbert變換比較 (a)原始信號及帶噪信號; (b)針對原始信號和帶噪信號的壓縮感知結(jié)果; (c)針對原始信號和帶噪信號的Hilbert變換結(jié)果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變換提取延時(shí)的結(jié)果.圖9(a)中黑線表示原始信號,原始信號由兩個(gè)不同延時(shí)的發(fā)射脈沖疊加組成,紅線為信噪比20 dB的帶噪信號.圖9(b)為壓縮感知提取結(jié)果,其中峰值對應(yīng)發(fā)射脈沖起始點(diǎn).圖9(c)為Hilbert變換提取結(jié)果,其中峰值對應(yīng)發(fā)射脈沖包絡(luò)頂點(diǎn).比較圖9(b)和圖9(c)可知,壓縮感知的提取結(jié)果受噪聲影響較小,延時(shí)提取準(zhǔn)確,而Hilbert變換在有噪聲條件下提取的信號峰值相比無噪聲條件下發(fā)生了顯著的偏移.壓縮感知是一個(gè)求解最優(yōu)解的過程,在計(jì)算中允許誤差項(xiàng)存在,因此壓縮感知法比Hilbert變換更適用于有噪聲條件下的延時(shí)估計(jì).且壓縮感知提取結(jié)果對應(yīng)脈沖起始點(diǎn),可直接用于波形到達(dá)時(shí)間的計(jì)算.皮質(zhì)骨仿真和實(shí)驗(yàn)結(jié)果均驗(yàn)證了壓縮感知提取延時(shí)的準(zhǔn)確性,厚度估計(jì)的相對誤差分別為4.9%和3.6%.如圖6(b)和圖8(b)所示,建立聲速模型的仿真重建結(jié)果和實(shí)驗(yàn)重建結(jié)果均可以看到清晰的皮質(zhì)骨上下界面,且相比圖6(a)和圖8(a)未建立聲速模型的結(jié)果,皮質(zhì)骨重建厚度準(zhǔn)確.與仿真重建結(jié)果比較,實(shí)驗(yàn)重建結(jié)果存在明顯偽影,這是實(shí)驗(yàn)過程中存在噪聲干擾導(dǎo)致的.
在未知實(shí)際皮質(zhì)骨尺寸的情況下,采用壓縮感知方法估計(jì)模型各層延時(shí),結(jié)合已知聲速可計(jì)算模型各層厚度,正確建立聲速模型.目前本文方法只適用于多層規(guī)則介質(zhì)模型.射線追蹤法可用于多層非規(guī)則介質(zhì)的時(shí)域重建.Qin等[30]研究了多層非規(guī)則介質(zhì)的頻域重建.后續(xù)研究將圍繞非規(guī)則骨板成像展開,在現(xiàn)有方法上進(jìn)行改進(jìn),解決水平方向聲速變化問題.
本文基于合成孔徑法采集超聲信號,通過壓縮感知計(jì)算延時(shí)參數(shù)建立聲速模型,采用PSM算法進(jìn)行圖像重建,實(shí)現(xiàn)皮質(zhì)骨超聲成像.仿真和實(shí)驗(yàn)結(jié)果表明,本文采用的壓縮感知方法可準(zhǔn)確估計(jì)皮質(zhì)骨厚度,從而建立多層聲速模型.重建圖像可觀察到清晰的皮質(zhì)骨上下界面,且形態(tài)正確.本文方法相比傳統(tǒng)超聲波束形成,解決了固定聲速重建圖像的挑戰(zhàn),對于皮質(zhì)骨超聲成像的發(fā)展有一定借鑒意義.相比于Hilbert變換,壓縮感知方法更適于計(jì)算帶噪實(shí)驗(yàn)信號的延時(shí)參數(shù).本文方法目前適于規(guī)則骨板成像,考慮到骨的多種形態(tài)及復(fù)雜結(jié)構(gòu),之后的工作應(yīng)當(dāng)建立更為實(shí)際的骨模型,尋找更加通用的成像方法.同時(shí)探究在體實(shí)驗(yàn)的可行性.