基于交錯網格時域有限差分法的經顱聚焦超聲模擬器實現

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

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

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

1 基本方程

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

對各項同性非均勻粘彈性介質，應力和應變之間的能量守恒關系為：

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

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

為刻畫彈性波傳播，還必須結合應力-應變之間的動量守恒關系：

方程（1）和（3）構成考慮顱骨橫波的經顱聚焦超聲場所用的Kelvin-Voigt 控制方程。

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

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

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

2 數值方法

2.1 三維顱腦模型的建立

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

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

顱骨聲參數具有較強的非均勻性，其與CT 值H 之間具有如下關系：

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

2.2 Kelvin-Voigt方程的交錯網格時域有限差分格式

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

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

3 模擬器的驗證

3.1 周期性波的傳播

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

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

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

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

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

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

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

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

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

3.2 圓形活塞軸線聲壓分布

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

圓形活塞聲源輻射聲場的數值計算域如圖4所示。紅線包圍區域為整個計算域，綠線和紅線之間區域為PML 域。數值模擬設置如下：a=2λ，2 個數值域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。為使計算結果收斂，時間步長設置滿足cmaxΔt/Δx=0.1，cmax是最大聲速。

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

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

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

4 結語

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