低頻超聲胸腔成像的有限元仿真模型

周英鋼 李運(yùn)杰

（沈陽工業(yè)大學(xué) 遼寧省沈陽市 110870）

1 引言

超聲成像由于可以進(jìn)行無創(chuàng)、重復(fù)診斷，并且其測量時(shí)間短，成為了體內(nèi)成像的主要方法之一[1]-[4]。傳統(tǒng)醫(yī)學(xué)超聲成像技術(shù)所使用的超聲在2-10MHz 范圍內(nèi)，通過超聲回波數(shù)據(jù)進(jìn)行成像。而在肺部醫(yī)學(xué)中，由于其內(nèi)部空氣的存在，導(dǎo)致2-10MHz 范圍內(nèi)的超聲無法滲透進(jìn)肺部。并且這種高頻超聲在在肺部發(fā)生強(qiáng)烈分散，無法提取到有用的回波信號(hào)，因此將低頻超聲用于透射式的胸腔成像成為一種可能[5]。

本文在COMSOL軟件中仿真了低頻超聲在胸腔中的傳播過程，在一發(fā)多收的激勵(lì)-接收模式下，得到了接收陣列的信號(hào)矩陣，為低頻超聲胸腔成像的應(yīng)用研究提供了一種仿真方法。

2 超聲波的胸腔傳播模型

2.1 模型簡化

與胸部和腹部的運(yùn)動(dòng)相比，本研究中所采用超聲波的波長較小。為了簡化問題，假定傳播介質(zhì)是靜止的，不考慮呼吸過程中肺部和心臟的運(yùn)動(dòng)。

使用comsol 軟件建立如圖1 所示的模型，忽略脊柱和心臟，為直觀清晰地顯示超聲波在胸腔中的聲壓分布，將雙肺合并為一個(gè)長軸長為0.15m，短軸長為0.1m 的橢圓，在肺部外設(shè)置環(huán)狀的骨骼肌，骨骼肌外部周圍均勻放置12 個(gè)超聲波換能器。

圖1：胸腔二維模型幾何結(jié)構(gòu)

2.2 超聲波的波動(dòng)方程

超聲波的波動(dòng)方程解釋了理想狀態(tài)下流體中超聲波的傳播過程，換言之，波動(dòng)方程描述了任意位置的質(zhì)點(diǎn)在任意時(shí)間下的運(yùn)動(dòng)情況。波動(dòng)方程的表達(dá)式如式（1）。

在式（1）中，介質(zhì)密度用符號(hào)ρ 表示，介質(zhì)運(yùn)動(dòng)速度用字母v 表示，聲壓用字母P 表示，E 為當(dāng)前穩(wěn)態(tài)條件下介質(zhì)的剛度常數(shù)。介質(zhì)密度、運(yùn)動(dòng)速度和聲壓受空間位置和時(shí)間t 共同作用。

3 基于COMSOL的仿真模型建立

3.1 幾何模型

通過在有限元軟件COMSOL Multiphysics 中添加“壓力聲學(xué)，瞬態(tài)”、“固體力學(xué)”和“靜電”等物理場，建立了如圖2 所示的二維胸腔模型。在胸腔外部均勻設(shè)置一圈12 個(gè)收發(fā)一體式的超聲波換能器，每兩個(gè)換能器間隔30°，通過水與胸腔耦合，仿真區(qū)域最外部設(shè)置為完美匹配層（PML）。

圖2：簡化后的二維胸腔模型

3.2 激勵(lì)模式

超聲波激勵(lì)源采用高斯調(diào)制的正弦波，表達(dá)式如式（2）。

其中，f=25k，gp1(t)為高斯脈沖，其標(biāo)準(zhǔn)差為2×10-4，峰值位置為5×10-5。

仿真時(shí)，采用一個(gè)換能器發(fā)射，正對(duì)的七個(gè)換能器接收，逐次激勵(lì)-接收的模式。發(fā)射換能器按順時(shí)針方向變換，依次進(jìn)行12次仿真計(jì)算來達(dá)到360°環(huán)形掃描。

3.3 網(wǎng)格剖分

在有限元仿真中，需要對(duì)所建模型進(jìn)行網(wǎng)格剖分，所剖分的網(wǎng)格尺寸和幾何形狀以及設(shè)置的時(shí)間步進(jìn)對(duì)仿真中數(shù)值計(jì)算的精度均有一定影響。本模型下，流體域（即骨骼肌和肺）部分的網(wǎng)格尺寸設(shè)置為波長的五分之一。完美匹配層和超聲波換能器采用映射的方式進(jìn)行剖分，這種剖分方式在一定程度上可以加快后期運(yùn)算速度。

模型網(wǎng)格剖分圖如圖3 所示。

圖3：模型網(wǎng)格剖分圖

3.4 時(shí)間步進(jìn)

COMSOL 軟件中確定時(shí)間步進(jìn)的方法有向后差分公式、廣義α兩種方式。向后差分公式算法會(huì)產(chǎn)生散射，波形畸變隨計(jì)算時(shí)間成正比，所以在瞬態(tài)計(jì)算中，通常使用廣義α 的方法。該方法使用之前五個(gè)時(shí)間步進(jìn)的解，這樣一來有效地避免了波形畸變，同時(shí)還能夠預(yù)測下一個(gè)時(shí)間步進(jìn)的仿真結(jié)果。本文中采用廣義α 的方法計(jì)算，時(shí)間步進(jìn)控制為0.5us。

3.5 數(shù)據(jù)采集過程

圖4：接收信號(hào)波形圖

1 號(hào)換能器首先開始發(fā)射超聲波，正對(duì)的4-10 號(hào)換能器接收超聲波數(shù)據(jù)；2 號(hào)換能器發(fā)射超聲波，正對(duì)的5-11 號(hào)換能器接收超聲波數(shù)據(jù)，以此類推，直到12 號(hào)換能器發(fā)射，3-9 號(hào)換能器接收。

采樣步長設(shè)置為2us，采樣時(shí)間設(shè)置為1000us，這樣能夠保留信號(hào)的完整性。經(jīng)過12 次發(fā)射之后，得到每個(gè)接收換能器的透射原始數(shù)據(jù)。圖4 分別為第1 號(hào)、第7 號(hào)、第4 號(hào)和第10 號(hào)換能器發(fā)射時(shí)，正對(duì)的7 個(gè)換能器接收的信號(hào)波形。

超聲波的傳播速度與介質(zhì)有關(guān)，介質(zhì)密度越高超聲波在其中的傳播速度就越快，從圖4 中能夠知道，超聲波換能器在相同的位置發(fā)射時(shí)，不同位置的接收信號(hào)幅值不同，并且不同位置接收時(shí)超聲波穿過的障礙物不同，使得換能器接收陣列中各點(diǎn)接收到的聲波信號(hào)有一定不同。并且超聲波在對(duì)稱位置發(fā)射時(shí)，其對(duì)應(yīng)的接收換能器波形也存在著對(duì)稱性。說明建立的模型符合理想情況。

4 結(jié)論

本文使用COMSOL Multiphysics 平臺(tái)建立了一個(gè)簡化版的胸腔環(huán)境，它被 12 個(gè)超聲波換能器均勻包圍在四周。采用一個(gè)換能器發(fā)射、正對(duì)的七個(gè)換能器接收的方式仿真了超聲波在胸腔中的傳播過程，并提取出了12 次發(fā)射-接收過程中超聲換能器接收陣列的聲波數(shù)據(jù)。在超聲斷層成像的研究領(lǐng)域中，圖像重建主要有兩種實(shí)現(xiàn)方式，其中一種應(yīng)用渡越時(shí)間實(shí)現(xiàn)，即通過某種方式得到發(fā)射信號(hào)和接收信號(hào)最大幅值處的時(shí)間差；另一種圖像重建方式應(yīng)用幅值衰減實(shí)現(xiàn)，通過測量得到發(fā)射信號(hào)和接收信號(hào)的最大幅值之比。本文所述的模型能夠得到渡越時(shí)間和幅值衰減，可以為低頻超聲胸腔成像的應(yīng)用研究提供參考。