乳腺超聲斷層成像系統的有限元仿真分析?

張 玉 張國軍 裴 毓 張文棟

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

0 引言

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

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

1 基于COMSOL的仿真分析

1.1 幾何建模

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

圖1 幾何模型Fig.1 Geometric model

圖2 復合式超聲換能器結構Fig.2 Structure of composite ultrasonic transducer

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

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

1.2 邊界條件

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

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

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

其中，n代表表面法線方向，utt是結構加速度，FA是作用在結構上的載荷(每單位面積的力)。

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

1.3 網格劃分

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

圖3 模型網格剖分圖Fig.3 Meshing of model

圖4 局部網格剖分圖Fig.4 Local meshing

2 計算結果與分析

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

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

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

由聲衍射理論[8]：

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

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

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

圖6 接收信號Fig.6 Received signal

3 圖像重建

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

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

圖7 等角扇束系統示意圖Fig.7 Schematic diagram of isometric fan beam system

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

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

4 結論

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