虛擬相控陣超聲散射CT成像中的散射面快速定位方法

宋文愛,郭禹姬,楊順民,陳以方

(1.中北大學軟件學院,太原 030051;2.清華大學機械工程學院機械系,北京 100084)

虛擬相控陣超聲散射CT成像方法不需要反投影重構,在一定程度上彌補了線性超聲CT方法由于線性假設或弱散射假設帶來的不足,提高了超聲檢測的空間分辨率。但是由于該方法需要綜合考慮幅值和相位信息,所以需要對虛擬相控陣中的每一個陣元發射、其余陣元(包括發射陣元本身)接收到的所有信號進行處理,信息量較大,從而影響了成像的速度。為了快速而且準確進行超聲CT成像,如何快速查找掃描區域(檢測區域)內存在的散射源,是提高成像速度和成像質量的關鍵[1-2]。

筆者提出了一種利用時空壓縮的方法,該方法可以快速地識別是否存在散射界面,并能確定散射面產生的信號在時間軸上的粗略位置,為提高CT成像速度奠定了一定的基礎。

1 虛擬相控陣超聲散射CT成像方法簡述

1.1 試驗系統方案設計

圖1所示的鈦合金CSK-IA試塊為檢測對象,圖中標出了試塊的尺寸以及陣列傳感器的布置區域以及可能引起散射的界面。陣列傳感器A具有49個陣元,陣元寬度為 1 mm。傳感器頻率為6.25 MH z,在被檢試塊中的波長近似為1 mm。界面1為半徑100 mm的圓弧界面,界面 2為深度91 mm的水平界面,界面3為高度6 mm的豎直界面,界面4為深度85 mm,寬度2 mm的水平界面,界面5為高度15 mm的豎直界面,界面6為深度100 mm的水平界面。

圖1 超聲散射CT的鈦合金CSK-IA試塊及探頭陣列布設

1.2 虛擬相控陣超聲散射CT成像方法

如圖2所示,假設目標點P為探測區域內任意一個虛擬掃描點(散射點),有N個陣元的矩形陣列探頭,當第i個陣元激勵,第j個陣元接收(j=1,2,…,N),可以采集到的N個波列為[3-4]：

圖2 超聲散射CT的相位反演成像方法原理示意圖

式中W i代表第i個陣元作為發射元時,所獲取的時間序列集(波列),W18用波形表示如圖3(a)。元素w i,j(n)中i,j=0,1,2,…,N-1;n=0,1,2,…,M-1,代表第i個陣元發射,第j個陣元接收的時間序列信號,M為時間序列的采樣點數。w18,15用波形表示如圖3(b)。當i=j時為自發自收時間序列。則信號矩陣中的單元wi,j(n)中所包含的空間點P(x,y)聲場的超聲波信號分量為：

式中w inij(n)為中心為lij,長度為2l的窗口函數。

一般l取為：

式中int()為取整函數;c為聲速;ri和rj如圖2所示,分別為陣元i和j到虛擬掃描點P(x,y)的距離;Δt為時間采樣間隔;λ為波長。

窗函數根據具體情況,一般作如下選擇：

發射陣列在空間一點產生的聲場,等于位于這一點的點聲源在相同延遲或相移的接收陣列產生的信號總和。則對空間任意虛擬掃描點P(x,y),其聲場近似可用下式計算：

根據需要對式(6)所表示的信號進行處理,取出成像參量,如最大幅值等,形成f(x,y),作為空間點P(x,y)的圖像灰度值。對探測區域的每一點(按一定的空間分辨率進行虛擬掃描得到的點)進行上述計算,即可實現對探測區域的相位反演初步成像。

上述方法由于多組信號求和,其統計平均效果消去了噪聲,如果目標點沒有散射信號到達探頭,則合成信號的相對幅值幾乎為0。

2 超聲波列時、空壓縮方法簡述

對于有N個陣元的探頭陣列,在探測區域內可獲得N個如圖3(a)所示的波列(時間序列集),每個波列由N個時間序列組成,全部點參與成像運算所需要的時間長,內存需求大,制約了成像速度。如果能實現散射源的快速定位,對于提高成像速度,減少偽像,提高成像質量至關重要[5]。

為了提高成像速度,進行散射源的快速定位[6],首先要將波列的冗余信息進行壓縮。

2.1 時間壓縮

當檢測或掃描區域內存在散射源時,接收信號的幅值增大,在相應的接收陣元上出現局部最大值點。最大值點出現以前的接收序列開始段或反方向最大值點出現以后的結尾段,被認為是冗余信號,可以進行壓縮。CT成像時被壓縮時間對應的空間點的灰度函數,按沒有散射源點接收幅值的均值成像。而掃描區域的大小由陣元大小以及陣元數大小決定。

如圖4所示,P為掃描區域內任意散射點,其與起始陣元間的距離為r0,與陣列正中陣元的距離為rN/2,與其距離最小的陣元之間的距離為ry,則與其距離最小陣元的陣元序列號為：

式中y0為起始陣元的坐標;y為與P點距離最小陣元的坐標;Δd為陣元之間的間距。

圖4 掃描區域內任意散射點P與其距離最近接收陣元關系示意圖

設陣列元共有N個,根據式(1),取wi,i(n),n=0,1,2,…,M-1;如果N力偶數,取i=0,N/2;如果N為奇數,取i=0,[(N-1)/2]+1。

對w i,i(n)(n=0,1,2,…,M-1)采用搜索法,尋找出現第一個最大值的時間點l0,lN/2,再根據式(3)可求得r0,r(N/2)。則根據圖4,可解得：

將式(8)代入式(7)可得iy。

對w iyiy(n)(n=0,1,2,…,M-1)采用搜索法尋找出現第一個最大值對應的時間點liy,將得到的時間點左移2l,得：

同理,可以沿時間軸的反方向搜索第一個最大值對應的時間點l01,l(N/2)1,再根據式(7)和(8),求得iy,對wiyiy(n)(n=0,1,2,…,M-1)沿時間軸的反方向,采用搜索法尋找出現第一個最大值對應的時間點l1iy,將得到的時間點右移2l,得：

l的計算方法見式(4)。

由于接收到的散射信號存在一定的相位差,為減小計算誤差、時間壓縮后不丟失有用的信息,所以nmin和nmax分別向左、向右移動2l。

這樣對于所獲得的所有波列的時間序列的時間值均壓縮到以n=nmin-2l為起點,n=nmax+2l為終點的范圍內,并以該范圍內的每個時間采樣點作為灰度圖像的坐標y。

2.2 空間壓縮

對如圖3(a)所示的波列經過時間壓縮后,再進行空間壓縮,以采樣時間點為坐標x,各陣元的相對起始陣元的位置(對于圖1陣列A,從左向右,以第一個陣元為起點,每個陣元的相對位置作為這一陣元所接收信號的坐標)為坐標y,對應的采樣信號幅值為灰度值f(x,y)。進行空間壓縮后,圖3(a)所示的波列轉換為二維灰度圖,如圖5所示。

圖5 圖3(a)所示的波列時、空壓縮后的灰度圖

由圖5可以看出,由散射源形成的散射波前一目了然。

對壓縮后的N幅灰度圖,分別進一步進行空間壓縮,以每一個發射陣元的接收信號為參考信號,進行互相關。即：

則mj為ρij(m)最大時對應的時間差,將Wi=(w i1(n),w i2(n),…,wiN(n))中的所有w ij(n)(i≠j,n=nmin,…,nmax)平移變為wij(n+mj)。圖5平移后如圖6所示。

圖6 圖5進行時間軸平移后的灰度圖

將上述移位后的二維信號進一步壓縮為時間軸上的一維信號,即：

圖6所示灰度圖移位壓縮后的波形見圖7。

圖7 圖6壓縮的一維信號圖

通過對壓縮信號的相鄰信號值比較求極大,極小值。即：

(1)W i(n)＞Wi(n-1),且Wi(n)＞W i(n+1)時,W i(n)為極大值。

(2)W i(n)＜Wi(n-1),且Wi(n)＜W i(n+1)時,W i(n)為極小值。

依次類推,直到找到最后所需的峰谷值。根據極大值對應的時間軸位置,即可快速準確地找到散射界面的對應位置。

3 試驗驗證

根據上述快速定位的方法,對圖1中界面2、界面4和界面6進行了初步的快速定位,其局部成像如圖8所示。

圖8 灰度成像

由圖8可以看出,初步定位基本符合實際散射面的位置。經過實際計算表明,所確定的散射面為圖1所示的界面2、界面4和界面6,誤差最大為2mm。

4 結論

提出的利用時空壓縮的方法,將虛擬相控陣所得到的超聲陣列信號集(波列),變換為位置、時間的二維灰度圖,再將其投影到時間軸上。通過研究一維灰度投影的特征,達到快速地識別是否存在散射界面,并確定散射面產生的信號在時間軸上粗略位置的方法可行。通過這種方法可以快速重建材料內部的缺陷斷層圖像。

[1]倪文磊.超聲CT理論與方法綜述[J].CT理論與應用研究,2004,13(1)：50-55.

[2]劉波,李朝榮.超聲CT成像方法及應用[J].中國儀器儀表,2007(2)：28-31.

[3]Rohde A H,Veidt M,Rose L R F,etal.A computer simulation study of imaging flexural inhomogeneities using p late-w ave diffraction tomography[J].U ltrasonics,2008,48(1)：6-15.

[4]楊奕.超聲相控陣檢查方法研究[D].北京：清華大學.2003.

[5]彭虎.超聲成像算法導論[M].北京：中國科學技術大學出版,2008.

[6]劉超.超聲層析成像的理論與實現[D].杭州：浙江大學,2003：6-7.