圓柱類構件超聲相控陣聚焦模型

孔惠元 王 鑒 李仰軍 韓 焱

(1 中北大學 山西省信息探測與處理重點實驗室 太原 030051)

(2 中北大學 山西省現代無損檢測工程技術研究中心 太原 030051)

0 引言

相控陣超聲檢測技術是近年來發展起來的新型超聲無損檢測技術，其基本原理是通過控制不同陣元的激發時間，使所有陣元在被測構件中產生的超聲場“相長干涉”和“相消干涉”，從而實現合成波束在預定位置聚焦[1]，在不移動換能器的前提下實現對被測構件大范圍檢測[2]。與常規單陣元換能器相比，相控陣換超聲檢測具有聲束可控、覆蓋率大、焦區能量強等優勢[3]。目前超聲相控陣在平面構件無損檢測中得到了廣泛的應用。然而，圓柱類構件在工程中大量存在，如固體火箭發動機、鋼管、傳動軸架等，在利用相控陣換能器對此類構件進行檢測時，受曲界面結構引起的入射波和回波時延的影響，掃描聲束的波陣面產生彎曲。目前在該類檢測應用中，往往采用迭代遍歷的方式進行延遲時間的計算，從而實現波陣面的控制。然而此類算法需大量迭代離散點，計算速度慢，且計算精度受離散點取值影響，從而導致實時性較弱、檢測效果差，影響了在該類構件檢測中的應用。為了提高相控陣儀器針對此類結構檢測時的效率，充分發揮相控陣儀器的優勢，本文建立了一種具有通用性的柱類構件相控陣波束掃描模型，該模型基于換能器、耦合介質、圓柱類構件的聲學、材料特性和幾何關系，通過輸入耦合介質和被檢測構件的聲速、曲率半徑及換能器陣列與曲面的距離，即可計算出換能器各陣元的收發時間，利用該時間進行換能器收發時間控制即可實現柱類構件的相控陣成像。

目前針對單一介質以及雙層平面結構的超聲相控陣聚焦延時法則已經非常成熟，對圓柱類構件的聚焦方案研究較少，主要難點是曲面使超聲波發生反射、折射等現象，已知聚焦點位置逆向求解聲束路徑時，聲束在界面的入射點確定困難，從而無法確定延遲時間。當前對超聲能量聚焦可以分為兩個方向，一是物理聚焦，通過相控陣換能器控制各激勵通道波形形狀、時間延遲來實現聲波能量在特定點的聚焦；二是后處理聚焦，通過后處理各換能器的檢測信號以實現能量在指定位置處的虛擬聚焦[4]。后處理方法意味著聚焦聲束不會物理地存在于被測構件中，而是通過用陣列數據的完整矩陣和相關成像算法合成其效果[5]。利用相控陣換能器控制波束在曲面構件內部聚焦已有部分學者進行了相關研究。徐娜等[6]將界面離散化，基于費馬定律對每一個離散點與陣元和聚焦點之間的聲程傳播時間進行計算，將最短傳播時間路徑作為聲線真實路徑，該方法需要遍歷界面離散點，計算量大；周正干等[7]基于折射定律對每一個入射角度進行遍歷，尋找折射聲線與聚焦點距離最近的聲線路徑進而確定聲線傳播時間，該方法同樣需要進行大量迭代計算，計算效率低；甘勇[8]指出，上述兩種求解方法的精度均取決于一個事先設定的增量?，?越小，解的精度越高，但求解速度越慢；高世凱等[9]將陣元垂直向下輻射聲波與界面的交點視為折射點，進而確定聲束路徑與傳播時間，該方法對折射點的確定可以理解為在曲界面上構建了多個虛擬源，只是一種假設，并不滿足折射定律；Zheng[10]研究了柔性相控陣聚焦方法，但柔性陣列換能器制造工藝復雜且價格昂貴，不具備普遍性[11]。利用回波信號對曲面構件進行虛擬聚焦也有部分學者進行了相關研究。關山月等[12]基于費馬定律對曲面結構全聚焦算法聲束傳播時間進行計算，該方法計算量同樣較大；Camacho等[13]將虛擬源與合成孔徑算法結合，通過陣元參數和近似幾何模型求取曲界面上虛擬源坐標，將聲源通過折射點到達聚焦點的路徑替換為聲源通過虛擬源到達聚焦點的路徑，進而實現聲傳播時間的計算與圖像重建，避免了復雜的折射點坐標計算，但要求工件曲率半徑遠大于陣元中心間距，否則求取虛擬源坐標誤差較大。

本文針對工業中使用一維線性相控陣換能器檢測圓柱類構件聚焦法則展開研究，基于換能器、耦合介質、圓柱類構件材料特性和幾何關系以及聲線模型和折射定律，建立了耦合介質及被檢構件的聲速、曲面曲率半徑、陣列與曲面間的距離等關聯的延遲時間聚焦控制模型。通過計算構件內部聚焦時陣元發射聲束在界面入射點坐標，得到各陣元發射聲波到達聚焦點時間得到延遲時間。經有限元仿真實驗，驗證本文聚焦方案在效率上的優勢與可行性。

1 圓柱類構件掃描波束的影響因素分析

1.1 雙層結構平界面相控陣成像檢測波束形成原理

雙層結構平界面聚焦聲束路徑計算示意圖如圖1 所示，首層厚度h，聲速c1，第二層介質聲速c2，且滿足c1

圖1 雙層結構平界面聚焦聲束路徑計算示意圖Fig.1 Schematic diagram of path calculation of the two-layer flat interface focused sound beam

設換能器陣元數為N，中心間距為p，陣元i發出聲束以入射角α入射到界面點a(m,h)，然后以折射角β到達位于第二層結構內的聚焦點F(xf,yf)，則求解點a坐標過程如下。

設陣元i位置坐標(xi,0)，由已知條件得

入射角α滿足

折射角β滿足

入射角α、折射角β滿足折射定律

將式(2)、式(3)代入式(4)得

通過解算式(5)即可求得a(m,h)，則聲束從陣元i發射到達聚焦點F(xf,yf)的總傳播時間TiF為

重復上述步驟，計算所有陣元到達聚焦點F對應聲時構成集合{TF}，則第i個陣元的延遲時間可以表示為式(7)，通過調整各陣元的不同延遲時間實現聲場在點F波束聚焦。

1.2 圓柱類構件相控陣成像檢測的影響因素分析

目前工業中對圓柱類構件的檢測通常采用水浸耦合或定制楔塊的方式將聲場能量傳遞到被測構件中進行檢測[14]。采用相控陣成像檢測時，對某一聚焦點，計算不同陣元發射聲束經耦合液(或楔塊)-構件界面的入射點及時間延時至關重要[15]。不同于雙層結構平界面入射點的計算，聲束在耦合介質-圓柱類構件界面的入射點橫縱坐標都是變量，且在界面的每一個入射點法線方向與x軸方向不垂直，隨著入射點的變化，法線方向發生偏轉，基于折射定律求解入射點時，陣元、入射點、聚焦點坐標和入射角、折射角難以通過直角三角形建立幾何關系。當前許多學者求取折射點坐標方法可以分為基于費馬定律(或折射定律)的遍歷法和在界面構建虛擬源法，而基于折射定律和構件幾何關系直接求取折射點坐標，進而確定延遲時間暫無相關研究。

2 圓柱類構件掃描控制波束時間

2.1 聚焦聲束路徑及延遲時間計算

超聲波經耦合介質入射圓柱類構件時，界面形狀(凸面或凹面)、耦合介質和被測構件聲速相對大小都會影響透射聲束的傳播特性，使聲束在圓柱類構件內部自動發生擴散或者聚焦現象，本節以液體為耦合劑、超聲波入射凸面構件為例推導聚焦聲束路徑計算方法及延遲時間、聲束入射臨界條件。

圓柱類構件聚焦聲束路徑計算示意圖如圖2 所示，被測構件表面為圓弧，曲率為k，曲率半徑為r=1/k，圓心為O，聲速c2，相控陣換能器位于被測構件正上方，以陣列換能器中心為原點、水平向右為x軸且為正方向、垂直向下為y軸且為正方向建立直角坐標系，換能器中心位置與被測構件中心軸線處界面垂直距離為h，中間位置充斥著液體耦合劑，耦合劑聲速c1，且滿足c1

圖2 曲面構件聚焦聲束路徑計算示意圖Fig.2 Schematic diagram of focused beam path calculation for curved component

設換能器陣元數為N、中心間距為p，則聚焦點F(xf,yf)的聚焦聲束路徑及延遲時間可通過如圖2所示的幾何關系求得。

(1) 基于被檢構件內聚焦點位置的曲面入射點求解

設陣元i發出聲束以入射角α入射到界面點a(m,n)，然后以折射角β到達位于構件內的聚焦點F(xf,yf)。盡管從陣元i到a再到F有無數條路徑，但是聲波只沿遵循折射定律的直線傳播，此時聲線從陣元i到聚焦點F的路徑是唯一的，即入射角α是唯一的，這決定了點a的坐標[16]。

設陣元i位置坐標(xi,0)，xi滿足式(1)，由圖2中幾何關系可知入射角α為向量(-m,h+r-n)和向量(xi-m,-n)夾角的補角，即

入射角α、折射角β滿足折射定律

將sinα、sinβ向量坐標表達式代入式(12)，并與點a所在曲率圓上方程聯立得方程組

通過解算式(13)即可求得a(m,n)。

(2) 解得入射點a的坐標后，聲束在耦合劑中聲程l1與傳播時間T1為

(3) 被測構件中聲程l2與傳播時間T2為

(4) 陣元i發出聲波到達聚集點F傳播時間TiF為

重復上述步驟，計算所有陣元到達聚焦點F對應聲時構成集合{TF}，則第i個陣元的延遲時間可以表示為

通過式(19)即可解出各陣元發射聲束在點F聚焦的延遲時間，同樣對于任一位置聚焦點坐標代入上述步驟都可求出各陣元的延時，通過調整各陣元的不同延遲時間實現聲場在預定點能量增強。本節以液體耦合劑為例進行聲束路徑計算及延遲時間的推導，對于楔塊耦合下延遲時間的計算同樣適用，此外，本文以晶片陣列中心與圓柱類構件圓心處于同一垂直方向為例進行推導，而在推導過程中并不嚴格要求此兩點處于同一垂直方向，只需獲取各陣元坐標即可實現聲延遲時間計算、能量聚焦，因此陣元中心與構件圓心存在偏置時也適用。

2.2 聲束入射臨界條件

陣元發射聲線與構件左右兩側會有一點相切，此時入射角為90?，兩切線范圍外的入射聲束不會進入構件，此外如圖2所示，入射角大于第一臨界角θ時折射聲束無法進入構件，設此時入射聲束與界面兩側交點為p、q，則有效聲束入射點位于p與q之間(這里不考慮波束角大小和波型轉換現象，假設每個陣元發射聲束可以覆蓋整個構件，以縱波聲線路徑進行研究，不考慮橫波影響)。設第i個陣元位置坐標(xi,0)，xi滿足式(1)，設臨界點p位置坐標(xp,yp)，點q坐標(xq,yq)，由折射定律有

同時入射角θ為向量(-xp,h+r-yp)和向量(xi-xp,-yp)夾角的補角，即入射角θ滿足

將sinθ向量坐標表達式代入式(20)，并與點p所在曲率圓上方程聯立得方程組

通過解算式(23)即可求得p(xp,yp)，在yp≤h+r范圍內臨界點p、q坐標同時滿足式(23)方程組，因此可以解出兩組實根，即為臨界點p、點q坐標。

2.3 入射點方程組求解

在上述延遲時間的計算中，重點在于對非線性方程組式(13)的求解，目前求解非線性方程組的方法有很多，如直接降維法、最小二乘法、牛頓迭代法等，本文使用數據處理軟件中相關工具包對方程組進行編程求解，算法思路及核心函數如圖3所示。

圖3 方程組求解算法流程圖Fig.3 Flow of equation system solving algorithm

fsolve 函數是通過選擇不同的參數來選取不同的算法，如牛頓迭代法或Levenberg-Marquardt 算法，再根據初始值和方程組的導數信息逐步接近方程組的解，本文通過選擇Levenberg-Marquardt 算法尋找未知變量的值。Levenberg-Marquardt 算法是使用最廣泛的非線性最小二乘法，它的關鍵是用模型函數對待估參數在其領域內做線性近似，忽略掉二階以上的導數項，從而轉化為線性最小二乘問題，它具有收斂速度快的優點。此外使用fsolve函數時提供一個合適的初始值非常重要，初始值選擇不合理可能導致算法無法收斂。初始值的確定與陣元坐標、聚焦點坐標、檢測對象參數、坐標系建立方式等有關，結合實際相控陣換能器各陣元中心間距較小(本文選用中心間距0.6 mm相控陣換能器)，則不同陣元對于初始值的選取影響較小，且本文以換能器中心建立坐標系，對于左右兩側的陣元，入射點分別位于構件中心左右兩側，因此本文選擇待檢測圓柱面上超聲換能器陣列中心對應位置(0,5)作為初始值。

基于費馬定律的迭代算法計算延遲時間時，入射點的坐標精度取決于界面的離散精度，離散精度越小，計算精度越高，但計算時間越長。分析本文的推導過程與方程組的求解過程可知，本文方法避免了費馬定律遍歷算法中界面離散精度的影響。然而在方程組求解過程中，所選算法不同、初始值選取不同也會影響方程組求解的精度與效率。實際檢測過程中，受陣元坐標測量誤差、換能器頻率存在微小偏差導致近場區計算誤差、構件聲速與理論值存在偏差等因素影響，也會導致聲束聚焦點坐標實際值與理論值存在微小偏差。

3 仿真驗證與結果分析

3.1 仿真條件

為了驗證本文方法在計算效率上的優勢與有效性，本文利用數據處理軟件進行數據計算，利用有限元仿真軟件對如圖4 所示的圓柱類構件進行仿真。模型結構的長30 mm、高20 mm，上半部分為耦合介質，下半部分為曲率半徑r=50 mm的半圓，材料為鋼，為節省計算時間，只截取36?圓心角對應弧長作為耦合介質與鋼的交界面。分別對耦合介質為水和有機玻璃楔塊進行驗證，材料參數如表1 所示，根據各向同性固體材料中縱波聲速計算公式

表1 仿真使用耦合介質與鋼的參數Table 1 Simulate parameters using coupling media with steel

圖4 仿真模型結構圖Fig.4 Simulation model structure diagram

可得縱波聲速，式(24)中E為楊氏模量，ρ為密度，σ為泊松比。相控陣陣列為工程中使用比較普遍的16陣元相控陣換能器，中心頻率5 MHz，陣元編號從左向右依次遞增，陣列參數為陣元寬度d=0.5 mm，中心間距p=0.6 mm，陣列與鋼的距離h=5 mm(指陣列中心與構件中心軸線處界面垂直距離)。

仿真中網格劃分大小會直接影響到計算的精度，網格劃分越小，計算結果越精確，但計算時間較長，對計算機性能要求較高，較大的網格則會導致計算結果不精確。結合本文研究內容，網格劃分滿足縱波波長包含12個網格即可，采用網格類型為自由三角形網格。由于超聲波在兩種介質中聲速不同，網格劃分需要分開計算，最大網格大小為

其中，fc為信號中心頻率，本文設為5 MHz。時間步進的選取同樣會影響計算精度，仿真中最大時間步進要求小于單個網格中超聲波傳播時間，本文設為T0/20，T0為激勵信號周期，滿足仿真要求。為減小模型中其他邊界反射波對聚焦聲場的影響，將耦合介質左右邊界、鋼中下邊界設為吸收層。

由于聲場只有在近場才能有效聚焦，相控陣線性陣列換能器近場大小L為

式(26)中，α為陣元寬高比系數，D為陣列孔徑大小。將陣列換能器參數代入式(26)可得近場區為16 mm，因此將聚焦點設于近場區附近。

3.2 水浸耦合仿真驗證

首先將耦合介質設為水，聚焦點設于被測構件中間軸線處，焦點深度距換能器下表面10 mm。根據本文方程組求解思路進行編程，同時為驗證本文求解模型計算效率的優勢，對本文所提方法與費馬定律分別進行編程。費馬定律的原理為設定某聚焦點，針對某一陣元，將該陣元與聚焦點橫坐標范圍內的界面離散化，依次計算局部界面各離散點與陣元和聚焦點之間距離之和并計算對應聲線傳播時間，所有聲線路徑中傳播時間最小值即為聲線真實傳播時間，對應路徑為真實傳播路徑，對比結果如表2所示。

表2 兩種算法計算結果對比Table 2 Comparison of the calculation results of the two algorithms

由表2 可知，本文方法在計算次數上遠遠小于費馬定律算法，且本文方法在計算效率上提升顯著，對所有陣元進行計算所得相關參數如表3所示。

表3 各陣元編號、水中聲程、被測構件中聲程和延遲時間Table 3 The number of each array,the underwater sound path,the sound path in the measured component and the delay time

作為對比，繪制了沒有施加延時的瞬態聲場，仿真超聲波傳播過程中的聲場分布，仿真結果如圖5 所示。圖5(a)表示所有陣元同時激發超聲波，圖5(b)表示隨著深度增加波陣面在鋼內逐漸擴散。

圖5 陣元無延時瞬態聲場分布Fig.5 The transient sound field distribution of the array without delay

將表3 中的延時依次加在相控陣換能器各陣元，使聲波同時到達聚焦點位置，仿真結果如圖6 所示。圖6(a)表示換能器各陣元按照表3 的延時依次激發超聲波，圖6(b)表示在5.2 μs 左右在預設聚焦點位置能量增強實現聚焦。繪制聚焦聲場鋼內中心軸線處聲壓分布如圖7所示(數據經過歸一化處理)，由圖7 可見，沿深度方向，聲壓逐漸增大，在10 mm 深度左右所有陣元發射聲波相位一致能量疊加實現聚焦，聲壓達到最大值，而后隨著深度增大陣元輻射聲場再次發散，聲壓幅值逐漸減小，仿真和理論基本吻合。

圖6 中心軸線處聚焦瞬態聲場分布Fig.6 Focus transient sound field distribution at the central axis

圖7 鋼內中心軸線聲壓分布Fig.7 Sound pressure distribution in the central axis of steel

將聚焦點設于換能器中心位置右側5 mm、深度10 mm (偏轉加聚焦)，根據本文原理計算得各陣元激發延遲時間，同時對本文方法與費馬定律算法分別進行編程，對比結果如表4所示。

表4 兩種算法計算結果對比Table 4 Comparison of the calculation results of the two algorithms

由表4 可知本文方法迭代次數較少，且計算效率有效提升，對各陣元進行計算所得參數如表5所示。

表5 各陣元編號、水中聲程、被測構件中聲程和延遲時間Table 5 The number of each array,the underwater sound path,the sound path in the measured component and the delay time

將表5 中的延時數據依次加到相控陣換能器各陣元，使聲波在預定點聚焦，仿真結果如圖8 所示。圖8(a)表示所有換能器各陣元按照表5的延時依次激發超聲波，圖8(b)表示在5.5 μs 左右在預設聚焦點位置能量增強實現聚焦。

圖8 偏轉聚焦瞬態聲場聲布Fig.8 Deflection focusing transient sound field sound distribution

由焦點坐標可得聚焦點位于換能器中心位置右側向深度方向偏轉63?，繪制該方向聲壓分布，如圖9 所示(數據經過歸一化處理)，由圖9 可知隨著深度增大，聲壓逐漸增大，在10 mm 深度左右所有陣元發射的超聲波相位一致能量疊加實現聚焦，聲壓達到最大值，之后隨著深度增加聲場再次發散聲壓逐漸減小，仿真結果與理論基本吻合。

圖9 偏轉63°方向聚焦聲場聲壓分布Fig.9 Deflection 63°direction focused sound field sound pressure distribution

此外計算了其他位置點聚焦的延時時間并進行聲場仿真驗證，與上述聚焦結果類似，聲場可以在指定位置實現聚焦。

3.3 楔塊耦合仿真驗證

將耦合劑換為有機玻璃制楔塊，由于偏轉聚焦情況更為普遍，對于只聚焦的情況這里不再闡述。將聚焦點設于換能器中心位置右側4 mm、深度10 mm，根據上述原理得各陣元延遲時間，同時對本文方法與費馬定律算法分別進行編程，對比結果如表6所示。

表6 兩種算法計算結果對比Table 6 Comparison of the calculation results of the two algorithms

由表6 可知本文方法迭代次數較少，且計算效率有效提高，對各陣元進行計算所得參數如表7所示。

表7 各陣元編號、有機玻璃中聲程、被測構件中聲程和延遲時間Table 7 The number of each array,the plexiglass middle sound path,the measured component middle sound path and the delay time

將表7中的延時數據加到換能器各陣元仿真聲場分布，如圖10 所示。圖10(a)所有換能器各陣元按照表7 的延時依次激發超聲波，圖10(b)表示在3.9 μs 左右所有陣元輻射超聲波在預設焦點實現聚焦。

由焦點坐標可得聚焦點位于換能器中心位置右側向深度方向偏轉68?，繪制該方向聲壓分布，如圖11所示(數據經過歸一化處理)，在10 mm深度左右聲壓值達到最大，之后隨著深度增加聲壓減小，與預期結果基本吻合。此外計算了其他聚焦點的延時法則并進行聲場仿真驗證，與上述聚焦結果類似，聲束都可以在指定位置實現聚焦。

圖11 偏轉68°方向聚焦聲場聲壓分布Fig.11 Deflection direction 68° focused sound field sound pressure distribution

4 結論

針對工業中相控陣超聲檢測圓柱類結構時延遲時間計算量大的問題，利用解析法對聲束聚焦延遲控制時間進行了理論推導并進行仿真驗證，建立了基于圓柱類構件特征和耦合介質特性的超聲相控陣掃描成像的聚焦模型?；趽Q能器、耦合介質、曲面的聲學、材料特性和幾何關系，利用聲線模型和折射定律建立了耦合介質及被檢構件的聲速、構件表面曲率半徑、陣列與曲面間的距離等關聯的延遲時間聚焦模型。利用該模型計算出延遲時間，對各陣元激發超聲波時間進行控制，進而實現掃描聲束的聚焦。本文以水和有機玻璃制楔塊作為耦合介質，對鋼曲面的相控陣聲束聚焦進行了仿真實驗，結果表明本文方案在計算效率上提升顯著，同時聲束可以較理想地在指定位置實現聚焦，驗證了本模型在計算效率上的優勢與有效性。