田淑愛 丁 婷 田志鑫 楊 錄
聚焦超聲治療[1?3](Focused ultrasound surgery, FUS)是利用超聲束良好的方向性、可穿透性和可匯聚性,將體外超聲聚焦到體內靶病變區域,通過熱效應、空化效應和機械效應等達到治療效果。其中空化效應是腫瘤熱消融、體外碎石.、超聲溶栓以及超聲藥物控制釋放等超聲治療的關鍵機制。為了控制和利用空化效應以實現FUS的高效精準治療,必須對FUS過程中的聲空化進行監控成像。
目前,現有的聲空化檢測和成像的方法主要是光學法和聲學法。光學法主要是通過高速/超高速攝影[4?5](High-speed photography, HSP)、聲致發光[6?7](Sonoluminescence, SL)以及聲致化學發光[8](Sonochemiluminescence, SCL)和光衰減法等。光學方法可以直觀記錄空化泡的行為和時空動態分布,但對于某一斷層面的空化分布無法給出具體信息,且受限于透明介質。聲學法包括主動空化成像[9](Active cavitation imaging, ACI)和被動空化成像[10](Passive cavitation imaging, PCI)。Farny 等[11]和Gyongy 等[12]使用線陣換能器,線陣陣元不發射只接收的方式進行數據采集空化信號,然后通過波束合成算法計算時間能量進行二維重建得到被動空化成像,然而對于提供準確的空間位置信息的一級重建算法有待提升。主動空化成像包括普通B 模式空化成像和超快速空化成像方法。由于聲空化具有瞬態特性,普通B 模式超聲成像采用逐線掃描模式,成像速度較慢,無法捕獲空化微泡瞬態特性。Gateau 等[13]提出了一種超快主動空化成像(Ultrafast active cavitation imaging, UACI)技術,該技術采用高幀速率的平面波發射模式,但是基于平面波的空化成像質量較差,需要依賴于后端波束合成算法以改進圖像橫向分辨率和信噪比。目前常用的后端波束合成算法包括延遲疊加(Delay-and-sum, DAS)算法、最小方差(Minimum variance, MV)算法和相干系數最小方差(Coherence factor based minimum variance,MVCF)算法。
空化微泡母小波(Cavitation bubble wavelet transform, CBWT)技術是近年來提出的一種新型平面波超快速超聲空化主動成像方法,它可在保留平面波高幀頻低能量的同時,提高空化噪聲比(Cavitation-to-noise ratio, CNR)。Liu 等[14]利用RPNNP 模型構造母小波,并且用波束合成后的信號進行小波變換得到空化圖像。Bai 等[15]基于特征空間的脈沖反轉次諧波和超快速超聲空化成像相結合,保證了波束合成圖像的質量,提高了圖像的分辨率。黃玉蓉等[16]在超聲診斷的造影成像中利用Doinikov 模式構建造影微泡母小波,先進行連續小波變換,再由小波系數代替原始回波信號,最后對其進行后續波束合成得到造影圖像。本文提出一種結合CBWT、波束合成和平方差值求和減影(Sum-of-squared differences,SSD)的超快速主動空化成像方法,以獲得高質量的空化圖像。首先,采用基于平面波的超快速超聲空化主動成像方法,獲取空化微泡原始回波信號;其次,基于RPNNP模型構建空化微泡母小波,并對空化微泡原始回波信號進行連續小波變換;再次,對獲取的小波系數進行多種波束合成算法,并結合SSD 減影算法,提高空化圖像質量。為驗證方法有效性,本文討論了CBWT、3種波束合成算法(DAS、MV、MVCF)以及SSD數字減影算法對空化圖像CNR 值的影響,實驗結果表明CBWT-MVCF-SSD空化圖像最優,其CNR值比DAS空化圖像高31.73 dB。

圖1 基于超快速主動空化實驗系統示意圖Fig.1 Schematic diagram of the ultra-fast active cavitation experiment system
實驗系統包括聚焦超聲系統和Verasonics 平面波超聲采集系統兩部分,兩者之間通過一個任意波形發生器進行時序控制。如圖1所示,聚焦超聲換能器固定在透明有機玻璃水槽(50 cm×40 cm×30 cm)的側壁,水槽中充滿自來水,并保持在室溫(20?C±2?C)。水槽壁和底部放置有吸聲材料,用以減少聚焦超聲聲束的多次反射干擾。聚焦超聲系統是由聚焦超聲換能器和功率放大器構成。聚焦超聲換能器是單陣元凹面探頭,其中心頻率為1.2 MHz,孔徑150 mm,聚焦超聲聲功率為72 W。聚焦超聲作用時間是1 ms,作用停止后觸發信號輸入到Verasonics 系統時,選擇的線陣探頭的中心頻率為5 MHz,帶寬為80%,陣元數為128,采樣率為40 MHz,發射并接收平面波信號,采集得到空化微泡原始回波信號,該回波信號在數值仿真軟件中進行處理。
本文提出結合CBWT、波束合成和SSD 減影的超快速主動空化成像方法,如圖2所示,給出了本文3 部分算法:CBWT 算法、波束合成算法(DAS、MV和MVCF)以及SSD數字減影。

圖2 CBWT 結合波束合成以及SSD 數字減影研究流程圖Fig.2 CBWT combined beam synthesis and SSD digital subtraction research flowchart
2.1.1 構建空化微泡母小波
本文所研究的空化微泡屬于無包膜的自由氣泡,并且假設初始微泡是靜止不動,微泡振動過程中一直保持球形而沒有形狀改變的理想狀態。RPNNP 模型[17]假設符合這一理論,其模型表達式為

利用四階Runge-Kutta 方法對與時間有關的半徑進行數值求解,可以得到微泡振動半徑隨時間的變化曲線R(t)。振動微泡輻射出的聲壓曲線P(t)計算公式如下:

其中,r是驅動聲壓到微泡中心的距離。
在常溫20?下,在液體水中的RPNNP 模型中各符號的含義以及仿真需要的參數如表1所示。根據空化泡尺寸分布[18]相關研究,這里將初始半徑設置為1.0 μm。
通過對公式(1)和公式(2)求解,可預測出空化微泡回波聲壓曲線,如圖3所示。將預測的聲壓曲線P(t)進行歸一化,作為空化微泡母小波。聲學參數是平面波傳輸的聲場分布參數,驅動聲壓是實際測量值,如圖3(a)所示。

表1 RPNNP 模型中符號的含義與仿真計算時使用的參數Table 1 Meaning of symbols in RPNNP model and parameters used in simulation calculation

圖3 空化微泡RPNNP 模型仿真結果圖Fig.3 Simulation results of cavitation microbubble RPNNP model
2.1.2 空化微泡母小波技術
本文采用基于平面波超快速超聲空化主動成像,平面波高幀頻采集到空化回波信號E(t),將空化微泡母小波φ(t)作為小波基,對回波信號E(t)進行連續小波變換。關于回波信號E(t)的連續小波變換可以定義為

基于空化回波信號與構建的空化微泡母小波具有相關性,本文將兩者進行連續小波變換,其實質是將尺度變化下的空化回波信號與空化微泡母小波信號進行卷積運算,得到一系列關于尺度的小波系數,小波系數將代替原始回波信號。本文參考Liu等[14]的相關研究,將小波尺度設置為23。
平面波發射的不聚焦性,導致空化圖像質量降低,為了提高CBWT 后的圖像質量,需要后端進行波束合成。本文采用DAS、MV 和MVCF 三種波束合成算法,DAS算法是通過對不同信道接收的回波信號進行特定的延遲再疊加求和,得到目標點的聚焦信號;MV 算法是在回波信號沒有失真的情況下,在特定方向上輸出能量最低,得到最優的加權值。
MVCF[19]波束合成是在最小方差的基礎上加入相干系數,進一步提高了圖像質量。選取目標點并根據目標點位置計算波束合成的有效陣元數大小,記為M。根據目標點位置對M個陣元通道信號計算延時,得到延時后的信號,記為X(t):

把陣元數為M的陣列分為長度為L的子陣,子陣的總數為M -L+1,計算X(t)空間平均相關矩陣R(t):

用R(t)+γI來代替R(t),實現對角線加載,其中I為單位矩陣,γ=?·trace{R(t)},trace{R(t)}為相關矩陣的跡,?為算法加入的空間噪聲與信號功率比。取a為單位向量,并用式(4)的R(t)計算最優加權系數w(t):

計算X(t)的相干系數(Coherence factor, CF),CF定義為相干方向的能量與陣元信號總能量的比值:


數字減影是一種基于B 超視頻數據的處理方法,通過減影法可以消除兩幅圖像之間的差異,從而消除了背景噪聲帶來的高回聲影響。本文采用平方差求和數字減影算法,即SSD 數字減影算法,計算公式如下:

其中,In(z,x)、Ib(z,x)分別表示聚焦照射不同時間的B 超圖像,z為深度,x為掃描方向,n為計算的窗寬。
圖4是尺度為23 時CBWT 前后3 種波束合成算法空化圖像對比結果。圖4(a)是空化原始數據直接成像效果圖。圖4(c)、圖4(e)和圖4(g)是CBWT前,3 種不同的波束合成算法對圖像的質量的影響,DAS 和MV 波束合成算法周圍噪聲影響較大;MV 算法提高了空間分辨率,但是沒有提高對比度;MVCF 算法是在MV 算法基礎上加入了相關系數,圖像空間分辨率和對比度均有所提高,空化圖像質量最佳。圖4(d)、圖4(f)和圖4(h)是CBWT后,CBWT 技術結合3 種波束合成算法后的空化圖像。將圖4(d)、圖4(f)和圖4(h)與圖4(c)、圖4(e)、圖4(g)進行比較,圖像的噪聲有所減少,空化圖像質量有所提高,說明CBWT 對周圍組織和噪聲有抑制作用。CBWT 是應用小波變換的解相關算法,將預測的回波信號與空化信號進行小波變換構建空化微泡母小波,得到的小波系數與原始的回波信號具有很高的相似性。波束合成算法效果中,在CBWT 前后,MVCF 算法最佳,整體最佳效果為CBWT-MVCF。
空化噪聲比(CNR)是評價空化圖像質量的重要指標,其表達式為

其中,Icavitation指的是ROIs 區域空化的平均強度,ROIs 區域選擇為空化微泡區域;Inoise指的是同等面積大小下周圍噪聲的平均強度。
表2是CBWT 前后空化圖像CNR 值的對比。通過表中的數據發現,使用CBWT 后,進行波束合成對于圖像CNR有一定提升。DAS算法的CNR值提高了0.27 dB,MV算法的CNR值提高了0.36 dB,MVCF算法的CNR值提高了1.17 dB。CBWT中構造母小波與空化微泡信號的相關性較高,CBWT 技術抑制了周圍噪聲和組織信號。通過研究分析,在CBWT后,CBWT-MVCF效果最佳,空化圖像質量較好。

表2 CWBT 前后空化圖像的CNR 值Table 2 The CNRs of cavitation images before and after CWBT
圖5是CBWT 對于3 種波束合成算法結合SSD 減影圖像的對比,圖5(a)、圖5(c)和圖5(e)是CBWT 前波束合成進行了SSD 數字減影,圖5(b)、圖5(d)和圖5(f)是CBWT后波束合成進行SSD 數字減影。SSD減影算法很大程度上抑制了周圍組織和噪聲,圖像分辨率得到提高,進一步提高了空化圖像的質量。

圖4 CBWT 前后不同波束合成算法空化圖像對比Fig.4 Comparison of cavitation images with different beam synthesis algorithms before and after CBWT
將SSD 數字減影后的圖像與DAS、 MV、MVCF 空化圖像進行對比,SSD 數字減影算法在消除了背景噪聲后的圖像明顯優于只單獨進行波束合成算法時的空化圖像。CBWT 結合波束合成算法與數字減影算法處理后的空化圖像,去除背景噪聲的干擾,空化圖像的質量同樣優于只進行波束合成算法的空化圖像。但與波束合成后進行數字減影的圖像效果差異較小,所以本文圖像質量評價指標參考表3中CNR值。CNR值顯示CBWTMVCF-SSD效果最佳。
表3是結合CBWT 和SSD前后空化圖像的CNR值。通過與表2進行比較,DAS-SSD、MVSSD 以及MVCF-SSD 三種算法,CNR 值分別提高了15.99 dB、13.51 dB、10.19 dB;CBWT-DASSSD、CBWT-MV-SSD、CBWT-MVCF-SSD 三種算法,CNR 值分別提高了16.34 dB、15.07 dB、17.71 dB。從數據可以看出,通過SSD 減影算法,提高了空化噪聲比,可以驗證SSD 算法是有效的。CBWT-MVCF-SSD 的CNR 值最高,比DAS 空化圖像高31.73 dB,三者相結合的效果最佳。

圖5 CBWT 對于不同波束合成算法結合SSD 減影圖像的影響Fig.5 The effect of CBWT on different beam synthesis algorithms combined with SSD subtraction images

表3 結合CBWT 和SSD 前后空化圖像的CNR 值Table 3 The CNRs of cavitation images before and after CBWT and SSD
圖6是基于CBWT-MVCF-SSD的空化動態監控序列圖,分別是聚焦超聲作用50 μs、1 ms、2 ms、5 ms 時的空化圖像,圖中高亮區域為泡群,顯示了不同聚焦超聲作用時刻空化分布情況,從而達到對空化活動的動態監控。

圖6 CBWT-MVCF-SSD 最佳效果不同聚焦時刻下空化泡群分布序列圖Fig.6 CBWT-MVCF-SSD distribution sequence diagram of cavitation bubble group at different focusing moments for best results
本文提出一種結合CBWT、波束合成和SSD數字減影的超快速主動空化成像方法,主要討論了CBWT、3 種波束合成算法(DAS、MV、MVCF)以及SSD數字減影算法對空化圖像CNR值的影響,實驗結果表明,CBWT結合波束合成算法可提高空化圖像的質量和CNR值,且CBWT-MVCF的效果最佳。進一步結合SSD 數字減影,CBWT-MVCFSSD 空化圖像最優,其CNR 值可比DAS 空化圖像高31.73 dB,可實現對空化活動的動態監控。