乳腺超聲斷層成像系統(tǒng)的有限元仿真分析?

張 玉 張國軍 裴 毓 張文棟

(中北大學動態(tài)測試省部共建實驗室 太原 030051)

0 引言

乳腺癌是全世界女性最常見的惡性腫瘤之一。近年來，乳腺癌的發(fā)病率呈上升、年輕化趨勢，已成為當前社會的重大公共衛(wèi)生問題，因此，乳腺癌的精確診斷和及時治療尤為重要。目前，進行乳腺癌篩查的方法主要有X 線鉬靶、X 射線計算機斷層掃描成像(X-ray computed tomography, X-CT)、磁共振成像(Magnetic resonance imaging, MRI)、多普勒超聲[1?2]。然而，這些方法大多具有輻射性或者給患者帶來疼痛。近年來，超聲計算機斷層掃描(Ultrasound computed tomography, USCT)作為一種新型的成像技術(shù)[3]，在乳腺癌早期精確診斷方面具有良好的應用前景。超聲CT 以超聲波在人體不同組織內(nèi)的聲波傳播速度差異和人體內(nèi)的線性衰減系數(shù)為基礎(chǔ)[4]，利用超聲波對物體進行透射得到投影數(shù)據(jù)，對組織內(nèi)部的參數(shù)進行重建，從而得到乳腺圖像。

COMSOL Multiphysics 軟件以有限元法為基礎(chǔ)，通過求解偏微分方程來實現(xiàn)真實物理現(xiàn)象的仿真。本文利用COMSOL 軟件進行乳腺超聲斷層成像仿真，采用一發(fā)多收的掃描方式，得到256個環(huán)形陣列換能器的聲壓分布，提取了接收信號，并利用等角扇束濾波反投影重建算法進行驗證。

1 基于COMSOL的仿真分析

1.1 幾何建模

通 過COMSOL Multiphysics 中“AC/DC 模塊”、“聲學模塊”和“結(jié)構(gòu)力學模塊”接口創(chuàng)建，建立了如圖1所示的二維幾何模型。仿真測試在半徑75 mm 的圓形區(qū)域內(nèi)進行，內(nèi)置為水，256個壓電超聲換能器以環(huán)形排列的方式均勻分布，每兩個換能器間隔1.4?，置于水中。在(?20 mm, 10 mm)處設置一個半徑為12 mm 的圓形區(qū)域來代表軟組織，在(23 mm,?30 mm)處設置一個半徑10 mm 的圓形鐵塊障礙物來代表腫瘤。由于超聲波在軟組織和腫瘤中的傳播速度接近，本次實驗為了直觀、清晰地顯示超聲波在乳腺組織中的聲壓分布，因此，用鐵塊來代表腫瘤。復合壓電超聲換能器采用鋁層-粘接層-壓電陶瓷(PZT-5H)-粘接層-鋁層的結(jié)構(gòu)，粘接層設置的目的是阻抗匹配，如圖2所示。

圖1 幾何模型Fig.1 Geometric model

圖2 復合式超聲換能器結(jié)構(gòu)Fig.2 Structure of composite ultrasonic transducer

仿真時，采用一個換能器發(fā)射、其余255個換能器接收，即一發(fā)多收的掃描方式進行，使得接收換能器覆蓋面積達到最大，獲得足夠多的投影數(shù)據(jù)，以進行高質(zhì)量的圖像重建。發(fā)射換能器按順時針方向變換，依次進行256 次仿真來達到360?環(huán)形掃描。環(huán)形分布的換能器相較于線陣分布，數(shù)據(jù)采集速度更快，避免了因機械旋轉(zhuǎn)而產(chǎn)生的干擾，實驗數(shù)據(jù)更精確。發(fā)射換能器發(fā)射一個頻率為3 MHz的正弦脈沖波，周期為5，幅值為40 V。材料屬性見表1。

表1 主要材料屬性Table 1 Main material property

1.2 邊界條件

波動方程是用來描述聲學、光學、電磁學的波動情況的一類重要的偏微分方程，用COMSOL 中“壓力聲學，瞬態(tài)”接口求解瞬態(tài)波動方程[5]：

其中，t表示時間，ρ是流體密度，c是聲速，pt是總聲壓，qd是偶極源，Qm是單極源，pb是背景壓力波。輻射邊界選擇平面波輻射，完美匹配層中的典型波速設為水中的聲速即1500 m/s，聲壓級參考壓力使用水的參考壓力，即Pref,SPL=1 μPa。

固體域和水域的雙向自動耦合在水域和換能器邊界處實現(xiàn)，邊界條件為

其中，n代表表面法線方向，utt是結(jié)構(gòu)加速度，F(xiàn)A是作用在結(jié)構(gòu)上的載荷(每單位面積的力)。

換能器的壓電效應在PZT-5域中實現(xiàn)，通過求解線性方程，使得固體力學方程和靜電方程相耦合，方程通過將應力和應變與電場和電位移相耦合來模擬壓電效應。

1.3 網(wǎng)格劃分

在有限元仿真中，需要對所建模型進行網(wǎng)格劃分，COMSOL Multiphysics中提供交互式網(wǎng)格劃分環(huán)境，每一個網(wǎng)格劃分操作都將添加到網(wǎng)格劃分序列中，當網(wǎng)格劃分序列中所有操作完成以后，得到最終的網(wǎng)格。時間步長以及剖分網(wǎng)格的大小和形狀很大程度上限制了數(shù)值計算的精度。在流體域中，每個波長使用5～6 個網(wǎng)格單元，通常設置最大單元尺寸為波長/5，即網(wǎng)格尺寸=波長/5。因為此模型中水、軟組織和鐵塊的波長不同，最大單元網(wǎng)格大小分別為0.1 mm、0.101 mm、0.345 mm。完成網(wǎng)格剖分之后網(wǎng)格頂點數(shù)為26738，三角形單元數(shù)為48102，邊單元數(shù)為10242，頂點單元數(shù)為1796，整個模型網(wǎng)格劃分和局部劃分如圖3和圖4所示。

圖3 模型網(wǎng)格剖分圖Fig.3 Meshing of model

圖4 局部網(wǎng)格剖分圖Fig.4 Local meshing

2 計算結(jié)果與分析

在網(wǎng)格剖分完成之后，選擇并行稀疏直接求解器(Multifrontal massive parallel sparse direct solver, MUMPS)進行求解計算，該求解器具有完整的線性系統(tǒng)矩陣，能夠?qū)崿F(xiàn)實數(shù)和復數(shù)的運算，將求解器中解的最大頻率設置為fmax,sol=3 MHz。

確定時間步進的方法有向后差分公式、廣義α。在瞬態(tài)計算中，通常使用廣義α的方法。因為向后差分公式算法會產(chǎn)生散射，波形畸變隨計算時間成正比；廣義α可以有效規(guī)避這種干擾，它用的是前5個時間步長的解，并且可以對下個時間步長的解進行預測[6]。本次計算采用廣義α的方法，手動控制時間步長，采用默認的時間步長即周期/60，這樣的計算結(jié)果精確高，能達到最佳性能。

仿真完成后，對仿真結(jié)果進行后處理，得到在乳腺超聲斷層成像中，超聲波在乳腺中的聲壓分布情況。圖5是第一個換能器發(fā)射后，時間為27.8 μs、41.7 μs、69.4 μs、97.2 μs時的分布云圖。軟組織、鐵塊的聲速和密度不同，其聲阻抗系數(shù)也不同，聲速、密度越大，聲阻抗越大[7]，超聲換能器接收到的信號越小。從圖5中可以看出，超聲波在水中逐漸衰減，當超聲波遇到鐵塊時，能量有明顯的衰減。

由聲衍射理論[8]：

式(5)中，k、λ分別代表聲波在介質(zhì)中的波數(shù)和波長，a代表障礙物半徑。聲衍射的強弱與聲波波長和障礙物尺寸相關(guān)，ka ?1 時，聲衍射現(xiàn)象弱；ka=1時，聲衍射現(xiàn)象強。超聲波頻率為3 MHz 時，為了避免聲衍射現(xiàn)象造成的干擾，障礙物尺寸應遠大于0.08 mm，所以本模型的設定符合理論，所接收的信號也更加精確。

圖5 聲壓變化云圖Fig.5 Cloud chart of sound pressure change

在超聲斷層成像中，主要是利用渡越時間和幅值衰減進行圖像重建，發(fā)射信號和接收信號最大幅值處的時間差為渡越時間，最大幅值之比為衰減的程度[9]。仿真結(jié)束后，得到每個接收換能器的終端電壓和透射原始數(shù)據(jù)。本次實驗利用幅值衰減進行圖像重建，一次發(fā)射可以得到255組接收信號，記錄下每個接收信號前5 個周期中的最大幅值，將幅值衰減取對數(shù)得到投影數(shù)據(jù)，經(jīng)過256次發(fā)射之后，得到一組完整的接收信號矩陣。圖6分別為第1 個、第64 個、第128 個和第192 個換能器發(fā)射時，正對面5個換能器接收的信號。超聲波在不同介質(zhì)中傳播速度不同，在密度高的介質(zhì)中傳播速度快，分析看出，同一個換能器發(fā)射，不同位置接收的信號幅值不同，而由于障礙物對超聲波的遮擋，使得換能器接收到的信號差別較大。

圖6 接收信號Fig.6 Received signal

3 圖像重建

等角扇束系統(tǒng)示意圖如圖7所示，射線源位置S0，探測器分布在圓弧上，各探測器單元性能一致，發(fā)射源中心線左右扇面張開的角度都相同，扇形束呈等角度分布，射線S0E的位置由(β,γ)決定，等角扇束重建公式[10]為

其中，重建圖像點M(r,?)，S0M長度為L，S0E的投影為pf(γ,β)，D表示S0到旋轉(zhuǎn)中心的距離。

圖7 等角扇束系統(tǒng)示意圖Fig.7 Schematic diagram of isometric fan beam system

本次實驗中256個超聲換能器等角度環(huán)形分布于乳腺四周，計算得出接收換能器旋轉(zhuǎn)角度增量為1.4?，發(fā)射換能器與任一接收換能器射線間隔λ為0.7?，符合等角扇束濾波反投影算法的重建規(guī)律，可以利用此算法進行圖像重建。將256 次發(fā)射后得到的接收信號矩陣代入算法中進行重建，結(jié)果見圖8?？梢钥闯觯亟▓D像精度較高，分析發(fā)現(xiàn)與待重建圖像即實驗設計的模型一致，能分辨出鐵塊的位置。

圖8 等角扇束濾波反投影重建圖像Fig.8 Reconstructed image by equal angle fan beam

4 結(jié)論

本文利用COMSOL Multiphysics 多物理場仿真軟件，模擬了一個256 個換能器呈環(huán)形分布的乳腺癌檢測環(huán)境。采用一個換能器發(fā)射、其余255 個換能器接收的方式進行仿真模擬，提取了256 組超聲換能器接收信號組成信號矩陣。通過等角扇束濾波反投影算法進行圖像重建，根據(jù)重建圖像，可以較清晰地分辨出軟組織和鐵塊障礙物的位置、形狀和大小，經(jīng)過比對，重建結(jié)果與仿真模型一致，且重建圖像精度較高。驗證了換能器呈360?排列、進行環(huán)形掃描的方式可行。