基于COMSOL的超聲層析成像仿真及圖像重建

顧建飛，蘇明旭，蔡小舒

?

基于COMSOL的超聲層析成像仿真及圖像重建

顧建飛，蘇明旭，蔡小舒

(上海理工大學顆粒與兩相流測量研究所，上海市動力工程多相流動與傳熱重點實驗室，上海 200093)

為了驗證三角形和四邊形準則的二值邏輯反投影(Binary Logic Back Projection，BLBP)算法在超聲層析成像重建過程中的可行性，以超聲層析成像的重建區域中同時含有圓柱形和橢圓柱形障礙物為例，利用COMSOL對成像區域內存在障礙物的聲場進行仿真，得到聲場聲壓分布。然后，通過提取和處理所獲取的聲壓分布，獲得信號矩陣。最后，將信號矩陣代入二值邏輯反投影算法進行圖像重建。仿真結果表明：當入射波的波長λ遠小于障礙物的尺寸時，聲衍射的現象并不明顯。重建結果表明：兩種準則在幾何近似理論下可行，且四邊形準則的效果優于三角形準則。

超聲波；有限元；圖像重建；超聲過程層析成像

0 引言

在冶金、化工、電力等工業過程中經常出現兩相流或多相流介質，例如在冶金行業的選礦廠、冶煉廠和礦山充填工藝生產過程中，常常要求對漿狀的混合液中的固體(即金屬)量進行測量，以利于成本核算、質量控制和金屬回收。作為一種用以測量多相流流動的過程參數的技術，超聲過程層析成像技術(Ultrasonic Process Tomography，UPT)在20世紀80年代中后期逐漸形成和發展起來，并受到各界的廣泛關注[1-3]。層析成像的理論要求只有在投影角無限多，且每個投影角度下的投影線足夠多的情況下才能獲得對象的“精確重建”。目前，扇形束掃描方式由于有效投影數據量大，已經得到了廣泛的應用。其主要的實現方式是采用一種半柱面式的壓電晶片，換能器的結構如圖1所示，最大特點是其聲輻射/接收的角度范圍為180°，其聲場可以覆蓋檢測區域。但是，超聲層析成像技術仍有如下方面的制約：超聲波發射/接收同步和掃描切換技術不夠成熟，實現超聲波快速成像較為困難。此外，若超聲波波長與障礙物尺寸相當, 會產生聲散射、繞射等現象，且隨機發生，給測量分析帶來難度。

COMSOL Multiphysics是一款大型數值仿真軟件，通過有限元方法模擬，在科研和工程中能用偏微分方程(Partial Differential Equation，PDE)描述各種問題。它是集前處理、求解器和后處理于一體的有限元軟件，提供友好的圖形用戶界面和建模工具，簡化了建模工作。所以本文的工作將利用COMSOL軟件模擬柱面壓電換能器在超聲層析成像技術中的聲場。在層析成像的信號提取階段，利用COMSOL多物理場耦合軟件聲學模塊模擬聲場，得到理論上陣列傳感器的聲壓分布，這對UPT實驗系統的搭建具有一定的指導意義，為信號矩陣的精確度提供參照作用。同時也可以對改進的二值邏輯反投影算法做進一步說明和驗證。

1 模型建立

1.1 幾何模型

本文采用如圖2所示的二維幾何模型，建立與實驗系統相同的測量區域。其中，區域C0安裝了16個換能器的陣列裝置；區域C1是半徑為55 mm的成像區域。區域C2是一個半徑為10 mm的圓形區域，圓心位置在(-0.01，0.02) m，表示的是聚四氟乙烯圓柱障礙物。區域C3是一個逆時針旋轉45°的長半軸為18 mm、短半軸為12 mm、中點在(0，-0.02) m的聚四氟乙烯橢圓柱障礙物。障礙物聚四氟乙烯的聲速0=5200 m/s，密度0= 7865 kg/m3。障礙物外部為介質水，聲速s=1473 m/s，密度s= 997 kg/m3。假定某時刻，根據實際測量過程中采用的扇形束掃描模式，位于最左邊的換能器為發射探頭，剩下其余的換能器為接收探頭。由于換能器晶片的發射/接收面是半徑為5 mm的半圓柱面，因此設定成像區域外邊界5 mm位置，即(-0.06，0) m處有一柱面超聲波源，頻率=2 MHz，則聲波在障礙物區域的波長0=0/=2.6 mm，在介質水的區域波長s=s/= 0.736 mm。這里，計算出的波長對后續設定該區域的網格單元尺寸尤為重要。

1.2 超弱變分公式(Ultra Weak Variational Formulation，UWVF)模塊

在固體和流體中，線性時諧波場可以由Helmholtz和Navier方程來描述。但是，由于波的振蕩特性難以對波動方程作數值逼近。在能接受的精確度范圍內，對Helmholtz或Navier方程問題作近似求解，需要有一個相對密集的空間描述。對于標準的數值計算方法，比如低階有限差分法(Finite Difference，FD)和有限元法(Finite Element，FE)，根據經驗，每個波長取10個點是合適的。然而高頻超聲，由于聲波波長非常短，因此要求將網格劃分得更細，這將帶來計算復雜性問題。針對高頻聲波，文中利用COMSOL獨有的超弱變分公式可以解決上述問題，和低階有限元分析相比，UWVF能顯著減少計算時間，具有明顯的優勢。這種優勢來自UWVF的單元屬性，UWVF單元屬性包含了求解自由空間Helmholtz波方程的信息。對于二維、三維Helmholtz問題以及二維Navier問題，該方法被證明是有效的[4]。

1.3 網格劃分

網格的劃分對于模擬計算的收斂性和準確性有較大的影響。要給出可接受的準確性，聲場模塊的超弱變分公式需要一個相對均勻的網格，并通過最大單元尺寸參數max控制，而不是通過明顯的幾何細節。COMSOL軟件的超聲散射模型指出，當使用默認的UWVF設置時，定義的max等于或稍微大于2λ(為正整數)時，會在收斂性和準確性之間給出較好的平衡，可以獲得較好的仿真結果。由于本文建立的模型中介質水和障礙物區域中的波長并不一致，需要在COMSOL中對這兩個不同區域設置并定義不同大小的max。根據上述原則，分別對介質水和障礙物兩個求解域的max設置為0.0015 m和0.0052 m。網格劃分結果如圖3所示。

1.4 求解過程

模擬采用了連續波(continuous wave)形式，聲源發出的柱面波與第二類零階Hankel函數成正比：

其中：=[(-0)2+(-0)2]1/2指成像區域內任意一點坐標(，)離聲源位置(0,0)的距離；s=2/s，表示介質中的波數；函數J和Y分別指第一、二類Bessel函數，i為虛部單位。成像區域的外邊界施加輻射邊界條件，吸收柱面波。此外，要取得可接受的準確性，聲場模塊的超弱變分公式采用默認20階超弱Helmholtz單元。對求解域進行求解即可獲得聲場分布。

2 仿真結果分析

本文算例的模擬結果如圖4所示，結果以聲壓級云圖表示，顯示了超聲層析成像系統中液態介質和固態障礙物的聲壓級分布情況。從圖中可以看出，超聲波在水中的擴散致使聲壓沿傳播路徑逐漸減弱。另外，每一個柱體對超聲波存在明顯的遮擋作用，柱體后向聲壓非常弱，形成陰影區。考慮到扇形束掃描的幾何近似理論：只有當入射波的波長遠小于目標物的尺寸時，才能忽略由目標引起的衍射[5]，即>> 1(=2π /，表示聲波在介質中的波數；是障礙物的等效半徑)。所以，滿足該條件的2 MHz超聲波所測障礙物的等效半徑須遠大于0.12 mm。由于仿真幾何模型中設定的障礙物尺寸均遠大于該頻率下所測障礙物等效半徑的最小值，所以從該圖中還可以看到，聲波存在的散射現象，非常微弱，符合幾何近似理論。作者在實驗中也注意到了類似現象，其表現為障礙物后方對應的傳感器仍可以接收到極其微弱的超聲信號，可以忽略不計。然后，提取和處理所獲取的該聲壓分布結果，得到理論上陣列傳感器的聲壓分布。通過對代入二值邏輯反投影算法的信號矩陣進行閾值控制，得到一個信號矩陣(，)，其中為所使用的換能器個數。那么，第個發射位置對應的第個接收位置的信號為：

3 圖像重建算法

目前，采用二值邏輯反投影反演算法采用了上述的幾何聲學特性，例如基于三角形判斷的反投影反演算法[6]，其主要思想如圖5所示。當1號換能器發射時，4號換能器被其傳播路徑中物體遮擋，則該探頭的信號經過二值化處理后置0。這樣，由換能器間隔的中點a、b、c三點組成的區域表示目標物影響區間，圖像矩陣對應區域增加權重“1”。在此基礎上，本文提出了一種基于四邊形判斷的二值邏輯反投影反演算法，如圖6所示，當1號換能器發射時，4號換能器被其傳播路徑中物體遮擋，則該探頭的信號經過二值化處理后置0。這樣，由換能器間隔的中點a、b、c、d四個點組成的四邊形區域表示目標物影響區間，圖像矩陣對應區域增加權重“1”，以此類推，待所有掃描信號處理完畢后通過閾值對圖像進行目標和背景分割，獲得重建圖像。分別按三角形判斷準則和四邊形判斷準則對其進行反演，得到的結果如圖7和圖8所示。

圖6 四邊形判斷法

將兩種方式的重建結果與仿真對象進行對比發現，聲波的散射和邊緣效應沒有對兩種算法的圖像重建帶來較大影響；障礙物的位置信息、形狀信息和尺寸信息，都較為準確。說明了二值邏輯反投影算法所用的幾何聲學理論在入射波的波長遠小于目標物的尺寸時具備可行性。圖7和圖8雖然可以看出兩種判斷準則的重建效果大致一樣，但仔細對比上述兩種算法，可以發現細微的差別。主要表現在：兩種準則都能較好地得到圓柱障礙物的截面形狀，但是四邊形準則得到的圓柱障礙物的尺寸更接近真實模型，而三角形準則測得的偏大；此外，四邊形準則得到的橢圓柱障礙物的截面信息也稍微接近真實模型。

4 結論

在超聲層析成像領域，采用COMSOL軟件建立不同形狀、尺寸的柱體障礙物的有限元模型，并利用超弱變分公式模擬仿真成像區域中聲場的各點聲壓級，提取各個換能器位置的聲壓級數值大小可獲得超聲信號矩陣。

結合三角形判斷和四邊形判斷準則的二值邏輯反投影算法的反演結果表明，兩種判斷準則的重建效果大致一樣，同時也驗證了當入射波的波長遠小于目標物時該反演算法的可行性。

[1] Ibrahim S, Yunus M A M, Khairi M T M, et al. A review on ultrasonic process tomography system[J]. Jurnal Teknologi, 2014, 70(3): 1-5.

[2] Abbaszadeh J, Rahim H A, Rahim R A, et al. Ultrasonic tomography system: optimizing the frequency in a metal pipe conveyor[C]//Applied Mechanics & Materials, 2013, 284: 572-576.

[3] Hafiz M, Rahiman F, Ruzairi, et al. Ultrasonic tomography - image reconstruction algorithms[J]. International Journal of Innovative Computing Information & Control Ijicic, 2012, 8(1): 527-538.

[4] Huttunen T. The ultra weak variational formulation for ultrasound transmission problems[M]. Kuopion yliopisto, 2004.

[5] Xu L J, Xu L A. Ultrasound tomography system used for monitoring bubbly gas/liquid two-phase flow[J]. Ultrasonics Ferroelectrics & Frequency Control IEEE Transactions on, 1997, 44(1): 67-76.

[6] 陳孝震, 蘇明旭, 蔡小舒. 液固兩相系中的超聲波層析成像實驗研究[C]//中國工程熱物理學會. 2012.

CHEN Xiaozheng, SU Mingxu, CAI Xiaosu. Liquid-solid two-phase system of ultrasonic tomographic imaging experiments[C]//Chinese Society of Engineering Thermophysics. 2012.

TheCOMSOL basedsound filed simulation and image reconstruction of ultrasonic tomography

GU Jian-fei, SU Ming-xu, CAI Xiao-shu

(Institute of Particle and Two-phase Flow Measurement, University of Shanghai for Science and Technology,Shanghai Key Laboratory of Multiphase Flow and Heat Transfer in Power Engineering, Shanghai 200093, China)

In order to verify the feasibility of the Binary Logic Back Projection (BLBP) algorithm of triangular and quadrilateral criteria in ultrasonic tomography reconstruction, a typical UPT reconstruction area, which simultaneously contains a cylindrical obstacle and an elliptic cylindrical obstacle,is taken as an example. COMSOL is used to make simulation of sound field when the obstacles exist, and an acoustic pressure distribution is obtained. Then, by extracting and processing the obtained sound pressure distribution, a signal matrix is formed. Finally, the signal matrix is substituted into the binary logic back projection algorithm for image reconstruction. It could be clearly observed that the effect of diffraction is insignificant from the computation results whenthe wavelength of the incident wave is smaller than the size of the obstacle. The reconstruction results show that the two criteria are feasible under the geometric acoustic, and the reconstructed quality of the quadrilateral criterion is better than that of the triangular one.

ultrasound; finite element; image reconstruction; Ultrasonic Process Tomography(UPT)

TB551

A

1000-3630(2016)-03-0231-04

10.16300/j.cnki.1000-3630.2016.03.009

2015-08-31;

2015-11-10

國家自然科學基金資助項目(51176128, 51206113)

顧建飛(1990－), 男, 江蘇海門人, 碩士研究生, 研究方向為超聲過程層析成像。

蘇明旭, E-mail: sumx@usst.edu.cn