臧佳琦 許凱亮4)? 韓清見 陸起涌 梅永豐 他得安
1) (復旦大學信息科學與工程學院, 上海 200433)
2) (復旦大學腦科學研究院, 上海 200433)
3) (復旦大學材料科學系, 上海 200433)
4) (上海市醫學圖像處理與計算機輔助手術重點實驗室, 上海 200433)
微小血管及其血流實時成像對監測生物體血氧代謝等具有重要意義.在無微泡造影劑的情況下, 傳統超聲多普勒技術仍較難實現高信噪比的微小血管成像.本研究提出了一種無造影劑增強的超快超聲脊髓微血管成像方法.本研究從基于多角度復合平面波的高幀頻成像技術出發, 提出基于特征值分解的頻率-幅值雙閾值濾波法, 從而將脊髓組織信號和微血流信號分離, 可實現脊髓內微血流的動態成像.在體成像實驗結果表明, 無超聲造影劑時, 超快超聲多普勒成像技術仍可獲得較為清晰的大鼠脊髓內微血流的實時圖像, 并能夠清晰地呈現脊髓受損所致的微血流缺失狀況.定量分析結果表明, 增大復合平面波角度數可有效提高圖像的信噪比.綜上, 超快超聲多普勒成像技術有潛力被應用于脊髓內微血管成像及功能實時監測與動態評價, 相關結果可為脊髓功能成像方法的研究提供借鑒.
脊髓損傷往往會引起脊髓內血管網絡受損和局部血流灌注量缺失, 進而引發繼發性功能損傷[1-3], 如下肢功能障礙等, 這會給家庭和社會帶來巨大負擔[4].近年來, 針對脊髓損傷的血流成像與功能診斷方法已成為國內外學界的研究熱點與重點[5-7].
目前對在體微小血管成像及分析血流動力學變化的方法包括電子計算機斷層掃描血管造影(computed tomography angiography, CTA)[8,9]、磁共振血管造影(magnetic resonance angiography,MRA)[10]和超聲多普勒成像技術[11,12].CTA[13]和MRA[14]的血管成像分辨率可達到50和100 μm,但以上兩種技術存在掃描時間長、電離輻射和設備笨重等不足, 限制了其在實時血流成像領域的應用.
超聲成像具有無輻射、便攜與快速的特點, 當前超聲多普勒成像技術已被廣泛應用于人體血流成像以及包括血管內斑塊與血栓等疾病的臨床診斷與分析[15,16].近十年來, 基于平面波成像方法的超快超聲多普勒成像技術得到了較快發展[17,18].相較于傳統的波束聚焦成像方法, 平面波成像可大幅度提高成像幀頻[19](由每秒幾十幀提升至上萬幀),從而可用于觀測血流量、流速和血管直徑的瞬時微小變化并進行血流動力學分析.顧名思義, 平面波方法所發射的是非聚焦聲波; 因而, 與傳統聚焦波束的成像方法相比, 單幀平面波成像存在信噪比(signal to noise ratio, SNR)低、對比度差的缺點.就此, 2009年Montald等[20]提出了多角度平面波相干合成的方法, 實驗結果表明該方法能夠在保證高成像幀頻的同時, 提高成像結果的SNR、對比度和分辨率.2011年, Mace等[21]提出基于平面波發射的超聲功能成像方法, 分別在刺激胡須和癲癇發作情況下對大鼠腦血流進行了高時空分辨率成像,依據血流變化動態監測結果識別相關腦功能活動區.2013年, Mace等[22]深入介紹了超快超聲成像方法, 并與傳統多普勒成像方法相對比, 通過大鼠腦的在體實驗證明了該方法能有效提高微小血管血流變化的監測靈敏度.2017年, Demene等[23]將超快超聲多普勒成像技術應用于新生兒的腦血流成像, 在新生兒睡眠狀態下觀察了活躍態和安靜態腦區的血流變化, 為腦功能超聲成像技術的應用與發展開拓了前景.
相較于腦血流成像, 脊椎復雜的骨結構特征和呼吸運動的干擾給脊髓內微血流超聲成像帶來挑戰[24], 2018年, Khaing等[7]應用超聲微泡造影劑對椎板打開條件下的大鼠脊髓血流進行了對比度增強超聲成像(contrast-enhanced ultrasound,CEUS), 成功觀察到了脊髓中的微血管結構.近期研究表明, 超快超聲可在無造影劑條件下實現腦微血管成像[23], 但該技術在脊髓內的微血流成像仍待深入.另外, 相關多角度成像方法、雜波濾除技術以及多普勒頻偏分析等研究仍待探討.
本文提出了一種無需注射超聲造影劑的超快超聲多普勒成像方法, 該方法基于超快超聲平面波發射與接收序列, 利用基于爆炸反射器模型(exploding reflector model, ERM)的頻率-波數域遷移(fk-migration)算法進行圖像重建; 并利用基于特征值分解(eigenvalue decomposition, EVD)的頻率-幅值雙閾值濾波法進行動態血流信號與靜態組織信號的分離, 最終實現了動態血流的功率多普勒成像.本文進行了仿體實驗和大鼠脊髓在體實驗, 平面波成像幀頻大于萬幀每秒, 獲得了脊髓內微血流的高時空分辨率的成像結果, 并探討了超快多普勒血流成像中復合平面波角度數對成像質量的影響.
方法流程示意如圖1所示, 首先發射多角度平面波序列對感興趣區域進行高幀頻成像, 應用fkmigration算法進行波束合成并對多角度成像結果進行相干疊加.隨后, 對一段時間內得到的圖像進行EVD處理, 采用頻率-幅值雙閾值濾波法進行雜波濾除, 從而提取動態血流信號.最終, 對多幀血流圖像進行疊加平均獲得功率多普勒成像結果.

圖1 方法流程示意Fig.1.Flow chart of the proposed method.
文中所用超快超聲成像序列如圖2所示.該成像序列由若干個子序列構成, 每個子序列發射一組N角度的傾斜平面波并接收反射回波, 不斷重復子序列從而進行連續成像.為避免連續兩次接收到的回波信號發生混疊, 脈沖發射時間間隔要大于超聲波在成像區域往返一次所需的最長時間, 該時間所對應的最高極限幀頻K, 其計算公式如下:

式中 Δt表示超聲波在成像區域對角線長度上的往返走時,

其中d表示成像區域的深度,L表示超聲換能器陣列的總長度,c表示超聲波的傳播速度.

圖2 超快超聲成像序列示意圖Fig.2.Schematic diagram of ultrafast ultrasound imaging sequence.
2013年, Damien等[25]將Stolt提出的fkmigration算法應用于平面波成像, 并結合ERM模型進行波束合成, 從而實現圖像重建.相較于傳統的延時疊加算法, fk-migration算法能夠提高成像質量和計算效率, 計算過程中應用快速傅里葉變換算法可在保持圖像高SNR和橫向分辨率的前提下進一步提高計算速度[25].
2.2.1 ERM模型超聲波在傳播過程中遇到一系列散射子會發生散射, 反向散射回波會被超聲換能器陣列接收.超聲波發射與回波接收模型如圖3(a)所示,t=0時刻起, 某一陣元受到脈沖激勵而發射超聲波, 聲波傳播至(x1,z1)處的散射子, 產生反向散射回波并最終由(x0,0 )處的陣元接收, 全過程時長為[25]

其中c表示超聲波的傳播速度,θ表示傾斜平面波的角度, 該角度可由控制換能器陣列的脈沖發射延時加以實現.
另外, 建立一個ERM模型, 如圖3(b)所示,該模型中假設處的散射子是t=0 時刻產生反向散射回波的二次聲源[26], 超聲波傳播的總時長為[25]

其中c? 表示ERM模型中超聲波的傳播速度.
為使ERM模型與上述超聲波發射與回波接收模型等效, 令(4)式計算所得時長與(3)式及其一、二階導數式相等, 通過解方程組可得如下參數變換[25]:

2.2.2 頻率-波數域波束合成
首先, 對波動方程進行傅里葉變換可得亥姆霍茲方程[25]

其中p(x,z,t) 表示 (x,z) 處t時刻的聲場,表示z軸方向上的波數分量, 由下式計算所得[25]:

其中f表示超聲波的頻率,表示經(5)式參數變換后的聲波傳播速度,kx表示x軸方向上的波數分量.(8)式的解為[25]

其中P表示聲場強度幅值.
由Stolt所提出的fk-migration算法中, 頻域映射方式為[25]

其中sgn函數為符號函數.頻率域傅里葉變換之后, 需對傾斜平面波的發射延時進行相位補償, 其二維傅里葉變換的公式為[25]

圖3 超聲波傳播模型 (a) 超聲波發射與回波接收模型; (b) ERM模型Fig.3.Ultrasonic propagation model: (a) Ultrasonic transmitting and echo receiving model; (b) exploding reflector model (ERM).

其中φ(x,z,f) 是對回波信號進行了時頻域傅里葉變換以及相位補償的結果,φ(kx,z,f) 是進行了頻率-波數域二維傅里葉變換的結果, 二維傅里葉逆變換的公式為[25]

超聲探頭接收到的全部回波信號由組織、血流回波信號和噪聲信號組成.應用基于EVD的頻率-幅值雙閾值濾波法[27]進行血流信號的提取時, 首先對多幀圖像進行特征值分解, 計算公式為[28]

其中A為 (nx×nz,nt) 的二維矩陣, 由(14)式大小為 (nx,nz,nt) 三維時空矩陣p?(x,z,t) 整合得到;AT表示A的轉置矩陣,λk表示第k個特征值,ek表示第k個特征向量.每一個特征向量的多普勒頻移fD計算公式為[28]

其中 P RF 表示脈沖重復頻率,R(1) 表示一個單位的延時自相關, 其計算公式為[28]

因為靜態組織信號對應特征值大、多普勒頻移低的特征向量成分, 動態血流信號對應特征值小、多普勒頻移高的特征向量成分[28], 所以可設定一個頻率閾值和一個幅度閾值, 將多普勒頻移小于頻率閾值且特征值大于幅度閾值的特征向量加以濾除.最后進行矩陣重構, 得到動態血流信號為p?flow(x,z,t).
當超聲波通過血管時, 一部分能量被紅細胞散射并由超聲換能器陣列接收反向散射回波.多普勒超聲成像通過重復發射超聲脈沖來監測血流中紅細胞的運動從而得到血流變化的信息.功率多普勒模式下平均強度的計算公式為[22]

多通道超聲數據發射與采集在可編程的Vantage 256相控超聲實驗平臺(Verasonics Inc,WA, USA)上進行.探頭型號為L22-14vX, 包含128個陣元, 中心頻率為15.625 MHz, 陣元間距為0.1 mm.信號采樣頻率為62.5 MHz, 數據經由總線傳輸至計算機.
仿體采用長6 cm、寬5 cm、高5 cm的瓊脂塊, 在6—8 mm深處制備直徑為2 mm的管路以模擬血管.用注射泵以0.1 mL/s的速度將混有散射子的懸濁液注入中管道中模擬血流.仿體測量中采用每秒3500次的超快超聲脈沖發射序列, 該序列由500個子序列構成, 各子序列由—10°—10°之間等間隔分布的7個傾斜平面波組成, 復合成像幀頻為每秒500幀.成像區域為15.0 mm × 12.7 mm的矩形成像區域.
在體實驗對象為Sprague-Dawley大鼠(雄性,7周齡, 體重約350 g).通過手術將T13-L2段椎弓板和棘突切除, 克服超聲經過椎骨所致能量衰減.用異氟烷麻醉大鼠(5%用于誘導麻醉, 2.5%—3.0%用于術中維持麻醉).隨后, 制作刺入型脊髓損傷模型.實驗結束時, 大鼠脫位處死.實驗中采用每秒14040次的超快超聲平面波發射序列, 該序列由520個子序列構成, 每個子序列發射—10°—10°之間等間隔分布的27個傾斜角度平面波, 因此, 復合成像幀頻為每秒520幀.實際成像區域為14.0 mm × 12.7 mm的矩形成像區域(實驗設定的成像深度為5.9 mm, 實際成像深度是設定成像區域的對角線長度).實驗裝置示意如圖4所示.

圖4 大鼠脊髓血流超聲成像實驗裝置示意Fig.4.Schematic diagram of experimental ultrasonic imaging set-up for blood flow of rat spinal cord.
對1 s內獲得的回波數據進行波束合成與多角度平面波相干疊加, 共得到500幀復合B模式圖像.為提取仿體血流信號, 對500幀圖像數據進行EVD處理和特征向量的多普勒頻移分析.特征值閾值選取會影響超快超聲功率多譜勒成像結果的對比度和SNR, 相關選取標準仍有賴于經驗性參數[28].通常, 動態血流信號對應多普勒頻移高且特征值小的信號成分, 靜態組織信號對應多普勒頻移低且特征值大的信號成分.圖5(a)為歸一化多普勒頻移對應的特征向量個數的直方圖, 縱坐標為特征向量個數.由圖可得, 低頻分量集中分布在歸一化多普勒頻移小于0.1的區間, 因此將截止頻率設為0.1以濾除低頻組織信號.在體實驗中也采用了該閾值參數.圖5(b)為數據集的特征向量與特征值的對應關系, 圖5(c)為特征向量與歸一化多普勒頻移的對應關系, 圖5(d)為特征值與歸一化多普勒頻移的對應關系, 圖中虛線為歸一化多普勒頻移為0.1的閾值線, 實線為幅值為—80 dB的閾值線.將多普勒頻移小于0.1且特征值大于—80 dB的特征向量濾除, 重建后的成像結果(第400幀)如圖6所示.圖6(a)—(d)分別是濾波前的仿體成像結果、濾波后的仿體血流圖、濾波后的仿體軟組織圖和仿體血流的功率多普勒成像結果.如圖6(d)所示, 通過上述雜波濾除方法可以得到6—8 mm深度處仿體血流的功率多普勒成像結果.

圖5 特征值分解與多普勒頻移分析結果 (a) 歸一化多普勒頻移對應特征向量個數的直方圖; (b)特征向量的特征值; (c) 特征向量的歸一化多普勒頻移; (d) 特征值對應的歸一化多普勒頻移Fig.5.Eigenvalue decomposition and Doppler shift analysis results: (a) Histogram of the number of eigenvectors corresponding to normalized Doppler shifts; (b) eigenvalues of eigenvectors; (c) normalized Doppler shifts of eigenvectors; (d) eigenvalues versus normalized Doppler shifts.

圖6 仿體血流成像結果(第400幀) (a) 濾波前的成像結果; (b) 濾波后的血流成像結果; (c) 濾波后的軟組織成像結果; (d) 功率多普勒成像結果Fig.6.Imaging results of the 400th frame of the phantom blood flow: (a) Original image before clutter filtering; (b) blood flow image after clutter filtering; (c) soft tissue image after clutter filtering; (d) power Doppler imaging result.
對1 s內得到的回波數據進行波束合成與多角度平面波相干疊加, 共得到520幀復合B模式圖像.對520幀圖像數據進行EVD處理和特征向量的多普勒頻移分析, 結果如圖7所示.圖7(a)為歸一化多普勒頻移對應特征向量個數的直方圖, 縱坐標為特征向量的個數, 相較于仿體血流, 實際血流有不同的流速, 因此其多普勒頻移的分布區間更寬.若頻率閾值選取過高, 會損失低速血流信號成分, 導致有效信息缺失; 本實驗采用與仿體實驗一致的經驗頻率閾值, 即取歸一化多普勒頻移為0.1.圖7(b)為數據集的特征向量與特征值的對應關系,圖7(c)為特征向量與歸一化多普勒頻移的對應關系, 圖7(d)為特征值與歸一化多普勒頻移的對應關系, 圖中虛線為歸一化多普勒頻移為0.1的閾值線, 實線為幅值為—130 dB的閾值線.將多普勒頻移小于0.1且特征值幅值大于—130 dB的特征向量濾除.圖8為多角度平面復合和濾波前后的成像結果.圖8(a), (b)分別為發射單角度傾斜平面波得到的波束合成后的結果和發射27個傾斜角度平面波相干復合成像的結果, 圖8(c)為濾波后的多幀脊髓血流圖像, 從中可見微血流的變化情況.
為討論復合平面波成像中傾斜平面波角度數在超快超聲多普勒血流成像中的作用, 分別設定每個子序列中發射傾斜平面波的個數為N= 3[—1°—1°],N= 9 [—3°—3°],N= 17 [—7°—7°],N= 27 [—10°—10°](傾斜平面波的角度在設定區間內均勻分布), 復合成像幀頻均為520幀/s, 成像結果如圖9所示.圖9(a)—9(d)的1—3圖分別為上述4種情況濾波之前、濾波之后和功率多普勒成像結果.當N= 3 [—1°—1°]時, 成像質量差, 無法辨別脊髓內部微血管分布(如圖9(a)).當N= 9[—3°—3°]時, 多普勒成像方法能夠對脊髓內部微血管結構進行成像, 但難以呈現完整的微血管結構(如圖9(b)).隨著成像角度數的增大, 當成像角度N= 17 [—7°—7°]和N= 27[—10°—10°]時, 在多普勒成像結果中可以看到較為清晰的微血管結構(如圖9(c),(d)), 其中微血管分布主要由靠近背側和靠近腹腔的兩部分血供構成, 該成像結果與組織病理學中脊髓血管的解剖結構形態相似[29].隨著角度數的增大, 圖像的對比度和SNR得到改善.并能觀察到水平方向—2.5 mm, 豎直方向1.5—2.0 mm處由細針刺入脊髓所致的局部微血流缺失現象.

圖7 特征值分解與多普勒頻移分析結果 (a) 歸一化多普勒頻移對應特征向量個數的直方圖; (b)特征向量的特征值; (c) 特征向量的歸一化多普勒頻移; (d) 特征值對應的歸一化多普勒頻移Fig.7.Eigenvalue decomposition and Doppler shift analysis results: (a) Histogram of the number of eigenvectors corresponding to normalized Doppler shifts; (b) eigenvalues of eigenvectors; (c) normalized Doppler shifts of eigenvectors; (d) eigenvalues versus normalized Doppler shifts.

圖8 基于多角度平面波復合成像的大鼠脊髓血流成像結果 (a) 單角度平面波發射成像; (b) 多角度平面波復合成像; (c) 雜波濾除結果.每一次發射超聲平面波的時間間隔為71.225 μs, 27次發射傾斜平面波與接收反射回波的總時長為1.923 ms, 對多角度信號經相干復合可獲得單幀超聲圖像, 其所對應的成像幀率為每秒520幀Fig.8.Blood flow imaging results of rat spinal cord based on multi-angle compounding method: (a) Beamforming results after a single emission; (b) multi-angle compounding images; (c) images after clutter filtering.The time interval between each emission is 71.225 μs.Each compounded frame is obtained using 27 steering-angle plane-waves within a period of 1.923 ms.Consequently the frame rate is 520 frames per second.

圖9 不同角度復合平面波成像結果對比圖(復合幀頻均為每秒520幀) (a) 3個角度[—1°—1°]傾斜平面波復合成像結果;(b) 9個角度[—3°—3°]傾斜平面波復合成像結果; (c) 17個角度[—7°—7°]傾斜平面波復合成像結果; (d) 27個角度[—10°—10°]傾斜平面波復合成像結果.1為單幀原始B模式圖像, 2為雜波濾除之后的成像結果, 其中可見微血流變化, 3為1 s內采集數據得到的功率多普勒血流圖(1, 2色標單位為dB, 3為歸一化數值的多普勒成像結果)Fig.9.Comparison of compounded images with different numbers of steering angles (composite frame rate is 520 frames per second): (a) Images compounded of data from emitting 3 [—1°—1°] steering plane-waves; (b) images compounded of data from emitting 9 [—3°—3°] steering plane-waves; (c) images compounded of data from emitting 17 [—7°—7°] steering plane-waves; (d) images compounded of data from emitting 27 [—10°—10°] steering plane-waves.Images labeled 1 are original B-mode images; images labeled 2 are imaging results after clutter filtering in which changes of blood flow can be observed; images labeled 3 are power Doppler images of micro-vessels (data was obtained within 1 s).
在圖9的功率多普勒血流圖中選定一個感興趣區域(如圖9(b)—(d)編號3的圖中矩形虛線框所示), 進而分析結果的對比分辨率.如圖10所示,分別將角度數N= 9 [—3°—3°],N= 17 [—7°—7°],N= 27 [—10°—10°]的成像結果中該感興趣區域放大.圖10(b)給出了圖10(a)中虛線深度處的幅度曲線, 進而比較微血管成像分辨率.由圖可得, 當N= 27時, 可由箭頭所指處的幅度峰值分辨兩相鄰微血管; 而當N= 9時, 因圖像分辨率降低, 故無法分辨此處的兩相鄰微血管.
由圖10可得, 隨復合平面波角度數增大, 血流信號相較于背景噪聲得到增強, 以下對圖像的SNR加以量化分析.當N= 27時最大的角度范圍為—10°—10°,N= 9時的角度范圍約為—3°—3°.從圖10(a)中選擇3個尺寸為0.2 mm × 0.2 mm的感興趣區域(編號為1, 2, 3), 白色方框為血流信號區域, 最臨近的灰色方框為參考背景區域.SNR的計算公式如下[30]:

其中Iflow(i) ,Ibkgd(i) 分別為血流信號區域和參考背景信號區域中第i個像素點的幅值,nflow和nbkgd分別為血流信號區域和背景信號區域像素點的個數, SNR是所選區域血流信號幅度與背景噪聲信號幅度均方根的比值.不同角度傾斜平面波的情況下, 圖10(a)中3個感興趣區域的SNR變化曲線如圖11所示.

圖10 對比分辨率隨復合平面波角度數增加的變化情況 (a) 當角度數N = 9, 17, 27時, 圖9(b)-(d)編號3的圖中矩形虛線框中圖塊的放大結果; (b) 當角度數N = 9, 17, 27時, 圖10(a)中虛線深度處的幅值曲線圖Fig.10.Change of contrast resolution with the increase of the number of steering angles: (a) Enlarged image blocks in the rectangular dashed box in No.3 figure in Fig.9 (b)-(d) (angle numbers are 9, 17, 27, respectively); (b) amplitude curve of the dotted line position in Fig.10(a) (angle numbers are 9, 17, 27, respectively).
如圖11所示, 3個感興趣區域的SNR都會隨著平面波角度數的增加而增加.對于1號感興趣區域, 發射9個角度[—3°—3°]傾斜平面波時, SNR為15.15 dB; 發射27個角度[—10°—10°]的傾斜平面波時, SNR為18.99 dB.類似地, 對于2號和3號感興趣區域, 當發射9個角度[—3°—3°]傾斜平面波時, SNR較小, 分別為13.66 dB和12.96 dB;當發射27個角度[—10°—10°]傾斜平面波時, SNR明顯提高, 分別為19.95 dB和18.49 dB.結果表明, 復合平面波角度數和角度范圍的增大可顯著改善多普勒成像結果的SNR.

圖11 SNR隨復合平面波角度數增加的變化結果Fig.11.SNR versus number of steering angles.
在實際應用中, 骨衰減和相位畸變對經椎骨的超聲成像影響不可忽略, 需要進行畸變校正與衰減補償.目前代表性的方法包括基于已知速度模型的相位畸變校正方法[31]和基于未知模型的全波反演技術[32]等, 相關技術已被用于經顱骨的腦組織超聲成像.后續研究可引入相關技術實現經椎骨的脊髓超聲功能成像.
本文將超快多普勒技術應用于大鼠脊髓微血流實時成像, 通過發射高幀頻、多角度平面波進行相干復合, 并采用fk-migration算法進行圖像重建.對高幀頻復合圖像進行EVD處理和多普勒頻移分析, 并利用頻率-幅值雙閾值濾波法進行雜波濾除和動態血流信號的提取, 最終得到了可反映大鼠脊髓微血流的高質量功率多普勒圖像.實驗中采用14040次/s的傾斜平面波發射頻率, 通過基于EVD的頻率-幅值雙閾值濾波法, 得到了脊髓微血管的高質量圖像.結果表明, 增大復合平面波角度數可有效改善圖像SNR, 與9個角度的傾斜平面波成像相比, 27個角度結果使得血流成像SNR值增強5 dB左右.在大鼠脊髓微血管的多普勒成像結果中, 能夠清晰地看到脊髓內部的微血管分布,進而可由局部微血流缺失觀察脊髓內微損傷區域.本文采用的超快超聲多普勒方法可在無需注射超聲造影劑的情況下, 實現脊髓內微血流成像.未來研究將圍繞著脊髓超聲功能成像展開, 特別是進行脊髓血流動力學分析和神經-血流耦合作用的研究.