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

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

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

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

1 虛擬相控陣超聲散射CT成像方法簡(jiǎn)述

1.1 試驗(yàn)系統(tǒng)方案設(shè)計(jì)

圖1所示的鈦合金CSK-IA試塊為檢測(cè)對(duì)象,圖中標(biāo)出了試塊的尺寸以及陣列傳感器的布置區(qū)域以及可能引起散射的界面。陣列傳感器A具有49個(gè)陣元,陣元寬度為 1 mm。傳感器頻率為6.25 MH z,在被檢試塊中的波長(zhǎng)近似為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試塊及探頭陣列布設(shè)

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

如圖2所示,假設(shè)目標(biāo)點(diǎn)P為探測(cè)區(qū)域內(nèi)任意一個(gè)虛擬掃描點(diǎn)(散射點(diǎn)),有N個(gè)陣元的矩形陣列探頭,當(dāng)?shù)趇個(gè)陣元激勵(lì),第j個(gè)陣元接收(j=1,2,…,N),可以采集到的N個(gè)波列為[3-4]：

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

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

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

一般l取為：

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

窗函數(shù)根據(jù)具體情況,一般作如下選擇：

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

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

上述方法由于多組信號(hào)求和,其統(tǒng)計(jì)平均效果消去了噪聲,如果目標(biāo)點(diǎn)沒有散射信號(hào)到達(dá)探頭,則合成信號(hào)的相對(duì)幅值幾乎為0。

2 超聲波列時(shí)、空壓縮方法簡(jiǎn)述

對(duì)于有N個(gè)陣元的探頭陣列,在探測(cè)區(qū)域內(nèi)可獲得N個(gè)如圖3(a)所示的波列(時(shí)間序列集),每個(gè)波列由N個(gè)時(shí)間序列組成,全部點(diǎn)參與成像運(yùn)算所需要的時(shí)間長(zhǎng),內(nèi)存需求大,制約了成像速度。如果能實(shí)現(xiàn)散射源的快速定位,對(duì)于提高成像速度,減少偽像,提高成像質(zhì)量至關(guān)重要[5]。

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

2.1 時(shí)間壓縮

當(dāng)檢測(cè)或掃描區(qū)域內(nèi)存在散射源時(shí),接收信號(hào)的幅值增大,在相應(yīng)的接收陣元上出現(xiàn)局部最大值點(diǎn)。最大值點(diǎn)出現(xiàn)以前的接收序列開始段或反方向最大值點(diǎn)出現(xiàn)以后的結(jié)尾段,被認(rèn)為是冗余信號(hào),可以進(jìn)行壓縮。CT成像時(shí)被壓縮時(shí)間對(duì)應(yīng)的空間點(diǎn)的灰度函數(shù),按沒有散射源點(diǎn)接收幅值的均值成像。而掃描區(qū)域的大小由陣元大小以及陣元數(shù)大小決定。

如圖4所示,P為掃描區(qū)域內(nèi)任意散射點(diǎn),其與起始陣元間的距離為r0,與陣列正中陣元的距離為rN/2,與其距離最小的陣元之間的距離為ry,則與其距離最小陣元的陣元序列號(hào)為：

式中y0為起始陣元的坐標(biāo);y為與P點(diǎn)距離最小陣元的坐標(biāo);Δd為陣元之間的間距。

圖4 掃描區(qū)域內(nèi)任意散射點(diǎn)P與其距離最近接收陣元關(guān)系示意圖

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

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

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

對(duì)w iyiy(n)(n=0,1,2,…,M-1)采用搜索法尋找出現(xiàn)第一個(gè)最大值對(duì)應(yīng)的時(shí)間點(diǎn)liy,將得到的時(shí)間點(diǎn)左移2l,得：

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

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

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

這樣對(duì)于所獲得的所有波列的時(shí)間序列的時(shí)間值均壓縮到以n=nmin-2l為起點(diǎn),n=nmax+2l為終點(diǎn)的范圍內(nèi),并以該范圍內(nèi)的每個(gè)時(shí)間采樣點(diǎn)作為灰度圖像的坐標(biāo)y。

2.2 空間壓縮

對(duì)如圖3(a)所示的波列經(jīng)過時(shí)間壓縮后,再進(jìn)行空間壓縮,以采樣時(shí)間點(diǎn)為坐標(biāo)x,各陣元的相對(duì)起始陣元的位置(對(duì)于圖1陣列A,從左向右,以第一個(gè)陣元為起點(diǎn),每個(gè)陣元的相對(duì)位置作為這一陣元所接收信號(hào)的坐標(biāo))為坐標(biāo)y,對(duì)應(yīng)的采樣信號(hào)幅值為灰度值f(x,y)。進(jìn)行空間壓縮后,圖3(a)所示的波列轉(zhuǎn)換為二維灰度圖,如圖5所示。

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

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

對(duì)壓縮后的N幅灰度圖,分別進(jìn)一步進(jìn)行空間壓縮,以每一個(gè)發(fā)射陣元的接收信號(hào)為參考信號(hào),進(jìn)行互相關(guān)。即：

則mj為ρij(m)最大時(shí)對(duì)應(yīng)的時(shí)間差,將Wi=(w i1(n),w i2(n),…,wiN(n))中的所有w ij(n)(i≠j,n=nmin,…,nmax)平移變?yōu)閣ij(n+mj)。圖5平移后如圖6所示。

圖6 圖5進(jìn)行時(shí)間軸平移后的灰度圖

將上述移位后的二維信號(hào)進(jìn)一步壓縮為時(shí)間軸上的一維信號(hào),即：

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

圖7 圖6壓縮的一維信號(hào)圖

通過對(duì)壓縮信號(hào)的相鄰信號(hào)值比較求極大,極小值。即：

(1)W i(n)＞W(wǎng)i(n-1),且Wi(n)＞W(wǎng) i(n+1)時(shí),W i(n)為極大值。

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

依次類推,直到找到最后所需的峰谷值。根據(jù)極大值對(duì)應(yīng)的時(shí)間軸位置,即可快速準(zhǔn)確地找到散射界面的對(duì)應(yīng)位置。

3 試驗(yàn)驗(yàn)證

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

圖8 灰度成像

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

4 結(jié)論

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

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

[2]劉波,李朝榮.超聲CT成像方法及應(yīng)用[J].中國(guó)儀器儀表,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].北京：清華大學(xué).2003.

[5]彭虎.超聲成像算法導(dǎo)論[M].北京：中國(guó)科學(xué)技術(shù)大學(xué)出版,2008.

[6]劉超.超聲層析成像的理論與實(shí)現(xiàn)[D].杭州：浙江大學(xué),2003：6-7.