剪切波對HIFU經顱聚焦形成溫度場影響的數值仿真研究

張艷秋 張 浩 孫天宇 菅喜岐

(天津醫科大學生物醫學工程與技術學院 天津 300070)

0 引言

高強度聚焦超聲(High-intensity focused ultrasound,HIFU)是將超聲波聚焦在靶區組織上使靶區組織溫度升高致死病變組織的治療技術，具有無創/微創、可重復治療等優勢，現已應用于子宮肌瘤[1]、前列腺癌[2]、乳腺癌[3]、胰腺癌[4]等軟組織腫瘤的臨床治療，HIFU經顱治療膠質瘤[5]、神經刺激[6]、運動障礙[7]腦部神經系統疾病也進入臨床試驗。經顱HIFU治療的挑戰是顱骨引起的超聲急劇衰減、顱骨與周圍軟組織阻抗差異大而引起的焦點失真、移位。2001年，Clement等[8]利用凹平面單陣元換能器對經顱超聲傳播路徑進行了數值仿真和離體人體顱骨實驗驗證，結果表明當超聲波入射角度大于20?時，HIFU經顱治療中剪切波的影響不可忽略。2004年，Concor等[9]利用開口直徑為300 mm的500陣元的半球狀相控換能器進行經顱聚焦，數值仿真結果表明，在未考慮剪切波條件下，大口徑半球形換能器可避免顱骨及其周邊組織損傷，并在顱內形成可治療焦域。2006年，White等[10]利用開口直徑為12.7 mm的單陣元換能器輻照人離體顱骨，通過垂直入射法測得縱波經顱傳播速度，測定結果表明，當入射角度大于40?時，超聲波傳播速度大于已測得的縱波傳播速度，剪切波的傳播可以提高能量傳輸效率。2011年，Pinton等[11]利用開口直徑為300 mm的512陣元相控換能器對HIFU經顱形成的溫度場進行了仿真研究，結果表明，在相同焦距條件下，顱內淺表組織內形成的可治療焦域體積小于深層組織。2015年，Jones等[12]利用開口直徑為300 mm的128陣元稀疏分布半球形相控換能器進行經顱聚焦，并基于CT進行像差校正，結果表明，在未考慮剪切波條件下，基于人顱骨CT的異質顱骨模型的像差校正優于均質顱骨模型。2015年，Ding等[13]在未考慮剪切波影響的條件下進行了HIFU經顱聲壓場和溫度場的數值仿真研究，研究結果表明，采用時間反轉法結合振幅調制相控，可以降低顱骨處的溫升。2018年，王祥達等[14]利用256陣元平面相控陣結合猴顱骨CT建立了經顱聚焦的數值仿真模型，結果表明，考慮剪切波的數值仿真模型形成焦域處聲強更大，當聚焦深度較大時可以不考慮剪切波。

本文建立由曲率半徑為150 mm的半球形相控換能器和人體頭顱CT數據構成的三維數值仿真模型，并利用時域有限差分法(Finite difference time domain,FDTD)結合Westervelt聲波非線性傳播方程、動量方程、運動守恒方程和Pennes生物熱傳導方程進行聲壓場和溫度場的數值仿真。分析研究不同聚焦角度、不同陣元數和相同輸入功率條件下剪切波對形成焦域面積的影響，為HIFU經顱治療方案制定提供理論依據。

1 基本方程式

1.1 聲壓場

1.1.1 腦組織和水內聲壓場

在水、人體軟組織和不考慮剪切波影響的骨組織內，采用Westervelt聲波非線性傳播方程[15]：

式(1)中：?為拉普拉斯算子；p為聲壓；cl和ρ分別為縱波聲速和密度；β=1+B/(2A)為非線性系數，B/A為介質的非線性參數；t為時間；擴散系數δ為

式(2)中：α為吸收系數，ω=2πf為角頻率，f為頻率。

1.1.2 骨組織內聲壓場

動量方程和質量守恒方程適用于含水70%的人體組織內的聲波傳播[16]。對于骨組織超聲波傳播可在質量守恒方程中加考慮剪切波的彈性項來實現[17]。未考慮骨黏度時骨組織內聲波傳播的動量方程和質量守恒方程[18]：

其中：υ為縱波傳播速度矢量；體積模量ζE=其中，剪切模量為剪切波傳播速度；彈性應力張量為

式(5)中：I為單位張量；形變張量e為

1.2 溫度場

Pennes生物體組織內熱傳導方程[19]為

式(7)中：Cr為比熱；T為溫度；rec為熱傳導率；q為單位體積發熱量；WB為血流灌注率；CB為血流比熱。

2 數值仿真模型及仿真參數

圖1為利用46歲男性健康志愿者人體頭顱CT數據和開口直徑為300 mm的256陣元半球形隨機分布相控換能器建立的經顱輻照數值仿真模型。其中，換能器的隨機分布陣列在相鄰陣元間距不小于1.0 mm的條件下依次選定各陣元中心位置，陣元半徑為8.5 mm。聚焦深度為45 mm，數值仿真區域大小為200 mm×300 mm×300 mm，聲軸為z軸。數值仿真的空間步長dx=dy=dz=0.3 mm，數值仿真的時間步長dt=10 ns，模型邊界采用Mur一階邊界吸收條件進行處理。水體、顱骨和腦組織各仿真參數[18,20]如表1所示。根據不同的聚焦角度θ在換能器基底選取部分陣元進行激勵，總功率為150 W時，陣元激勵參數如表2所示。

圖1 256陣元半球形隨機分布相控換能器經顱輻照數值仿真模型(單位：mm)Fig.1 Numerical simulation model of transcranial irradiation of 256 elements hemispherical randomly distributed phase controlled transducer(Unit:mm)

表1 水體、顱骨和腦組織仿真參數Table 1 Water,skull and brain tissue simulation parameters

表2 陣元激勵參數Table 2 Array element excitation parameters

3 結果

以圖1所示數值仿真模型，利用時間反轉法獲取陣元激勵信號，結合相位調控和幅值補償進行經顱輻照，輸入總功率為150 W，陣元激勵頻率選用0.7 MHz。

3.1 剪切波的影響

圖2 不同聚焦角度換能器在幾何焦點處聚焦的焦平面溫度場Fig.2 Focal plane temperature field with different focus angle transducers focused at geometric focus

不同聚焦角度換能器在形成焦域的最大溫度達到65?C時的溫度場如圖2所示。由圖2可知，隨換能器聚焦角度θ減小，在幾何焦點處形成的焦域面積逐漸增大，達到65?C所需時間逐漸延長；當聚焦角度60?6θ6 150?考慮剪切波的影響時，相同的聚焦角度條件下，溫度場達到65?C所需時間tmax更短，旁瓣更少，顱骨處溫升不明顯；當聚焦角度θ=30?時，考慮剪切波的溫度場與未考慮剪切波的溫度場達到相同溫度所需時間相同；當聚焦角度θ<90?時，焦點前移且在顱骨處產生較高溫度；當聚焦角度θ=30?時，顱骨處溫升較大。圖3為與圖2對應的聲軸溫度曲線。由圖3可知，在相同聚焦角度下，考慮剪切波形成的溫度場在顱骨處的溫度更高。圖4為與圖2對應的不同聚焦角度下焦平面面積和輻照時間。由圖4可知，隨著聚焦角度的增大，考慮剪切波形成的焦域面積小于未考慮剪切波形成的焦域面積。圖5為與圖2對應的不同聚焦角度換能器形成的焦域最高溫度與預設焦點在聲軸上的偏離距離。由圖5可知，聚焦角度減小，考慮剪切波與未考慮剪切波形成最高溫度在聲軸上的偏離增大，當聚焦角度θ>90?時，兩者間幾乎無差異。

圖3 不同聚焦角度換能器在幾何焦點處聚焦的聲軸溫度曲線Fig.3 Sound axis temperature curve of different focus angle transducers focused at geometric focus

圖4 不同聚焦角度換能器在幾何焦點處聚焦的焦平面面積柱狀圖和輻照時間折線圖Fig.4 Focal plane area histogram and irradiation time line graph of different focus angle transducers focused at geometric focus

圖5 不同聚焦角度換能器形成的焦域與預設焦域在z軸上的偏離距離Fig.5 Deviation distance between the focal region formed by the transducers with different focus angles and the preset focal region on the z-axis

3.2 聚焦角度與焦域面積的相關性

基于圖4的數值仿真數據，利用冪指數函數形式，在Matlab平臺上得到的擬合式為

其中：S1為考慮剪切波的焦域面積，S2為未考慮剪切波時的焦域面積。

圖6為基于圖4的焦域面積的仿真值擬合曲線。由圖6可知，考慮剪切波的擬合值與仿真值之間的和方差為3.215，均方根為0.8965，可決系數為0.9999；未考慮剪切波的擬合值與仿真值之間的和方差為64.01，均方根為4，可決系數為0.9941。考慮剪切波情況下，焦域面積的擬合值與仿真值之間的和方差更小，均方根更接近于0，確定系數更接近于1。

圖6 不同聚焦角度換能器在幾何焦點處聚焦形成的焦域面積擬合曲線Fig.6 Focal domain area fitting curve formed by different focus angle transducers at geometric focus

4 結論

本文基于半球形相控換能器經顱輻照模型，建立了未考慮剪切波傳播和考慮剪切波傳播的兩種經顱輻照的數值仿真模型，并對激勵不同聚焦角度陣元的情況下形成的溫度場進行研究，得到以下結果：

(1)隨換能器聚焦角度θ減小，在幾何焦點處形成的焦域面積逐漸增大，焦點前移程度越大，考慮剪切波形成的溫度場達到65?C所需時間逐漸延長。

(2)在相同聚焦角度條件下，考慮剪切波的溫度場達65?C所需時間更短，旁瓣更少。

(3)隨換能器聚焦角度θ減小，考慮剪切波的模型形成的焦域面積變化范圍更大。

(4)經顱考慮剪切波情況下，焦域面積的仿真值與擬合值之間的和方差更小，均方根更接近于0，確定系數更接近于1。

由上述結果得到如下結論：

(1)在輸入功率和焦域最高溫度相同的條件下，聚焦角度越小的換能器形成焦域越大，焦點前移距離越大。

(2)考慮剪切波的傳播可以提高能量傳輸效率，但對焦點前移幾乎沒有影響。

(3)剪切波在顱骨處傳播會導致顱骨內熱沉積相對增多。

(4)冪指數函數形式擬合優度高，可預測不同聚焦角度換能器形成的焦域面積。

綜上所述，考慮剪切波的溫度場的溫升速度更快，形成焦域面積更小且旁瓣數量更少，但顱骨處熱沉積更多。對于不同患者以及同一患者不同位置上的腫瘤治療時，應根據腫瘤位置與相控換能器之間的相對位置以及腫瘤大小等，考慮剪切波對形成焦域的影響，以保障HIFU腦腫瘤治療的有效性和安全性。

5 討論

將圖1中的半球形換能器和頭顱結構細化為圖7所示x-z與y-z平面圖，設定頂骨部分為同心圓弧。聚焦角度θ在30?～150?之間時，最大入射角度λmax為8?～30?，隨聚焦角度增大，最大入射角度增大，顱骨處相對溫升減小且旁瓣減少。這與2013年Narumi等[21]通過改變超聲波從水中入射丙烯酸板的入射角度，提出的超聲波入射角度在0～34.3?時，隨超聲入射角度增大，剪切波透射率增加的結論一致。

(1)圖1所示的9～193陣元相控換能器數值仿真模型中，當聚焦角度為60?～90?時，超聲最大入射角度為8?～20?之間，剪切波對形成溫度場的影響不可忽略，這與文獻[8]的單陣元換能器進行經顱輻照所得結果不一致。

(2)本文基于人體真實顱骨CT和256陣元半球形相控換能器，在最大入射角度為8?～30?范圍內進行了研究，其研究結果與文獻[14]結果相一致。

(3)本文利用最小二乘法對不同聚焦角度換能器經顱聚焦形成焦域面積進行擬合時，出現與實際不符的負數；選用反比例函數和冪指數函數法進行擬合時，冪指數函數形式擬合精度更高。

(4)本研究采用一名志愿者頭部CT數據建立三維數值仿真模型，研究了剪切波對HIFU經顱聚焦形成溫度場影響。而不同患者的顱骨厚度、密度、曲率的差異很大，為了研究不同患者之間的差異性，下一步將導入多位志愿者的頭顱CT數據進行相關數值仿真模型，通過統計的方法研究顱骨厚度、密度、曲率等對剪切波傳播的影響。

圖7 超聲最大入射角度計算示意圖Fig.7 Schematic diagram of the maximum incident angle of ultrasound