基于交錯網(wǎng)格時域有限差分法的經(jīng)顱聚焦超聲模擬器實(shí)現(xiàn)

王祥達(dá) 滕世國 李建霖 呂根品

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

經(jīng)顱聚焦超聲是一項(xiàng)蓬勃發(fā)展的腦部疾病無創(chuàng)治療手段。其在顱內(nèi)腫瘤消融、神經(jīng)調(diào)控、血腦屏障打開、顱內(nèi)靶向給藥、顱內(nèi)血栓溶解等方面具有應(yīng)用前景。當(dāng)前，經(jīng)顱聚焦超聲的主要障礙是顱骨與其周圍介質(zhì)較大的聲阻抗差異、顱骨較強(qiáng)的非均勻性及聲吸收特性所帶來的顱骨及其周圍組織過熱和焦點(diǎn)偏移或散焦。目前主要通過數(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 方程，實(shí)現(xiàn)一套經(jīng)顱聚焦超聲場模擬器。

1 基本方程

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

對各項(xiàng)同性非均勻粘彈性介質(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ù)具有較強(qiáng)的非均勻性，其與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)驗(yàn)計(jì)算，其值可以設(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 模擬器的驗(yàn)證

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)格點(diǎn)間距Δx=Δy=Δz。在超過半個周期的計(jì)算之后，檢驗(yàn)歸一化聲壓誤差(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ù)值計(jì)算誤差隨著空間步長的增大而增大，聲壓的數(shù)值計(jì)算結(jié)果在空間上確實(shí)是四階精度的。

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

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

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

3.2 圓形活塞軸線聲壓分布

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

圓形活塞聲源輻射聲場的數(shù)值計(jì)算域如圖4所示。紅線包圍區(qū)域?yàn)檎麄€計(jì)算域，綠線和紅線之間區(qū)域?yàn)镻ML 域。數(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。為使計(jì)算結(jié)果收斂，時間步長設(shè)置滿足cmaxΔt/Δx=0.1，cmax是最大聲速。

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

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

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

4 結(jié)語

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