王祥達 滕世國 李建霖 呂根品

（1.乳源瑤族自治縣東陽光實業(yè)發(fā)展有限公司 廣東省韶關(guān)市 512721）（2.乳源東陽光機械有限公司 廣東省韶關(guān)市 512721 3.韶關(guān)東陽光自動化設(shè)備有限公司 廣東省韶關(guān)市 512721）

經(jīng)顱聚焦超聲是一項蓬勃發(fā)展的腦部疾病無創(chuàng)治療手段。其在顱內(nèi)腫瘤消融、神經(jīng)調(diào)控、血腦屏障打開、顱內(nèi)靶向給藥、顱內(nèi)血栓溶解等方面具有應(yīng)用前景。當(dāng)前，經(jīng)顱聚焦超聲的主要障礙是顱骨與其周圍介質(zhì)較大的聲阻抗差異、顱骨較強的非均勻性及聲吸收特性所帶來的顱骨及其周圍組織過熱和焦點偏移或散焦。目前主要通過數(shù)值模擬來預(yù)測顱骨給經(jīng)顱聚焦超聲療法帶來的困難，并給出相應(yīng)的解決辦法。傳統(tǒng)的經(jīng)顱聚焦超聲模擬采用的控制方程主要是HIFU 中常用的粘性流體Westervelt 方程。該方程只考慮縱波傳播，并未考慮顱骨中存在的橫波。近年來，經(jīng)顱聚焦超聲模擬對顱骨橫波的考慮逐漸增多，包括利用交錯網(wǎng)格時域有限差分法數(shù)值求解Biot 方程、利用時域偽譜法數(shù)值求解Kelvin-Voigt 方程所建立的考慮橫波的經(jīng)顱聚焦超聲模擬器。本文則利用交錯網(wǎng)格時域有限差分法數(shù)值求解Kelvin-Voigt 方程，實現(xiàn)一套經(jīng)顱聚焦超聲場模擬器。

1 基本方程

1.1 各向同性非均勻粘彈性介質(zhì)中超聲波的傳播

對各項同性非均勻粘彈性介質(zhì)，應(yīng)力和應(yīng)變之間的能量守恒關(guān)系為：

T、ε、V 是應(yīng)力、應(yīng)變、振速。λ、μ 是拉密常數(shù)，χ、η 是縱、橫波粘性系數(shù)，滿足：

ρ0是密度，cp、cs是縱、橫波波速，αp、αs是縱、橫波吸收系數(shù)，ω=2πf0、f0是角頻率、頻率。

為刻畫彈性波傳播，還必須結(jié)合應(yīng)力-應(yīng)變之間的動量守恒關(guān)系：

方程（1）和（3）構(gòu)成考慮顱骨橫波的經(jīng)顱聚焦超聲場所用的Kelvin-Voigt 控制方程。

1.2 完美匹配層（PML）吸收邊界

本文使用不分裂卷積完全匹配層(CPML)作為模擬吸收邊界。對空間微分算子操作如下：

參數(shù)κx≥1、αp≥，dx是衰減因子，i 是虛數(shù)單位。滿足迭代式：

2 數(shù)值方法

2.1 三維顱腦模型的建立

三維顱骨模型的建立是基于顱骨CT 掃描文件，如圖1所示。

圖1：基于顱腦CT 掃描的三維顱腦模型建立

顱骨聲參數(shù)具有較強的非均勻性，其與CT 值H 之間具有如下關(guān)系：

φ是顱骨孔隙率，ρskull,min、cp,skull,min、αp,skull,min是顱骨中最小的密度、縱波聲速、縱波吸收系數(shù)；ρskull,max、cp,skull,max、αp,skull,max是顱骨中最大的密度、縱波聲速、縱波吸收系數(shù)。對于顱骨內(nèi)橫波聲速和吸收系數(shù)，基于經(jīng)驗計算，其值可以設(shè)置為：cs=4/7cp，αs=90/85αp。

2.2 Kelvin-Voigt方程的交錯網(wǎng)格時域有限差分格式

本文使用空間四階精度、時間二階精度的交錯網(wǎng)格時域有限差分法數(shù)值求解Kelvin-Voigt 方程。不同時刻應(yīng)力和振速分別定義為Tij|n、Vi|n+1/2，下標(biāo)i、j 代表笛卡爾坐標(biāo)，上標(biāo)n 表示時間坐標(biāo)(t0+nΔt)，Δt 為時間間隔。方程（1）、（3）結(jié)合CPML 的時間離散格式為：

Tφij、Vφij是CPML 吸收邊界引入的吸收應(yīng)力、振速的輔助變量。

3 模擬器的驗證

3.1 周期性波的傳播

利用周期性波的傳播上述模擬器。忽略，方程（1）、（3）退化為無黏流體聲波方程：

p=-Txx=-Tyy=-Tzz是聲壓，Txy=Tyz=Tzz=0。

在均勻流體中，方程（8）一組周期性波的解析解為：

i、j、k 是x、y、z 方向單位向量，k 是波數(shù)，p0和V0滿足：p0+ρ0cpV0=0。

方程（9）顯示聲壓和振速具有周期性空間分布，可將數(shù)值域設(shè)置為棱長為一個波長的正方體，同時用周期性邊界條件取代CPML 吸收邊界。f0=1MHz，ρ0=1000kg/m3，cp=1500m/s，V0=1m/s，網(wǎng)格點間距Δx=Δy=Δz。在超過半個周期的計算之后，檢驗歸一化聲壓誤差(pn-pa)/ ρ0cp，pn是聲壓數(shù)值解，pa是方程（9）中的聲壓解析解。以整個數(shù)值域中所有時刻[1,nt]最大2-范數(shù)和∞-范數(shù)表征該誤差，其定義分別為：

設(shè)置時間步長滿足Δt/Tp=0.0001 并保持不變，Tp=1/f0是周期。歸一化聲壓誤差的2-范數(shù)和∞-范數(shù)隨空間步長的變化及其擬合曲線如圖2所示。聲壓的數(shù)值計算誤差隨著空間步長的增大而增大，聲壓的數(shù)值計算結(jié)果在空間上確實是四階精度的。

圖2：歸一化聲壓的誤差的2-范數(shù)和∞-范數(shù)隨著空間步長的變化

同理，設(shè)置空間步長滿足Δx/λ=0.01并保持不變，λ= cp/f0是波長。歸一化聲壓誤差的2-范數(shù)和∞-范數(shù)隨時間步長的變化及其擬合曲線如圖3所示。聲壓的數(shù)值計算誤差同樣也是隨著時間步長的增大而增大，聲壓的數(shù)值計算結(jié)果在時間上也確實是二階精度的。

圖3：歸一化聲壓的誤差的2-范數(shù)和∞-范數(shù)隨著時間步長的變化

3.2 圓形活塞軸線聲壓分布

無限大平面障板上半徑為a 的圓形活塞聲源，圓心位于坐標(biāo)原點，活塞表面以Vz|z=0,ρ≤a=Vacos(ωt)均勻振動并向+ 軸輻射聲場。(ρ, φ,z)為柱坐標(biāo)系坐標(biāo)，Va是活塞表面振速幅值。則該圓形活塞聲源在軸線上(x=y=0，z ≥0)的歸一化（歸一化因子為2ρ0cpVa）聲壓分布為：

圓形活塞聲源輻射聲場的數(shù)值計算域如圖4所示。紅線包圍區(qū)域為整個計算域，綠線和紅線之間區(qū)域為PML 域。數(shù)值模擬設(shè)置如下：a=2λ，2 個數(shù)值域Lx*Ly*Lz=12λ*12λ*27λ & 12λ*12λ*31λ}，空間步長Δx=Δy=Δz=λ/10，PML厚度3λ，f0=1MHz，ρ0=1000kg/m3，cp=1500m/s，Va=1m/s。為使計算結(jié)果收斂，時間步長設(shè)置滿足cmaxΔt/Δx=0.1，cmax是最大聲速。

圖4：圓形活塞聲源輻射超聲場的數(shù)值計算域示意圖

Lz=27λ 和Lz=31λ 的計算域數(shù)值模擬得到的圓形活塞聲源軸線上的歸一化聲壓幅值分布與方程（11）的理論解析分布結(jié)果對比如圖5所示。數(shù)值模擬結(jié)果和解析結(jié)果具有很好的一致性，PML 吸收邊界對邊界內(nèi)聲場的吸收效果非常好，幾乎沒有數(shù)值反射波。

圖5：圓形活塞聲源的軸線上的歸一化聲壓幅值

4 結(jié)語

本文基于大腦的CT 掃描文件重建三維顱腦模型，分別利用空間四階精度、時間二階精度的交錯網(wǎng)格時域有限差分算法求解Kelvin-Voigt 方程，從而建立了經(jīng)顱聚焦超聲場模擬器。利用周期性聲波的傳播問題和圓形活塞輻射聲場問題對經(jīng)顱聚焦超聲場模擬器進行了核實。結(jié)果表明，模擬器確實對應(yīng)空間四階精度、時間二階精度的時域有限差分法；圓形活塞聲源軸線上聲壓分布的模擬值與解析解的一致性很好；CPML 吸收邊界很好的消除了聲波在計算區(qū)域邊界的反射。從而證實了經(jīng)顱聚焦超聲場準(zhǔn)確性和有效性。