基于超聲波斷層圖像的三維重構(gòu)方法

馬世博，劉宏偉，遲永波，王奕博，張 昭

(1.河北科技大學(xué) a.材料科學(xué)與工程學(xué)院，b.河北省材料近凈成形重點(diǎn)實(shí)驗(yàn)室，石家莊 050018；2.上海新孚美變速箱技術(shù)服務(wù)有限公司，上海 201100；3.河北京津冀再制造產(chǎn)業(yè)技術(shù)研究有限公司，河北 滄州 062455)

工業(yè)CT(computed tomography)技術(shù)又稱工業(yè)計(jì)算機(jī)斷層掃描成像技術(shù)，是一種重要的無(wú)損檢測(cè)技術(shù)，可以實(shí)現(xiàn)零件的斷層掃描及模型重構(gòu)，具有檢測(cè)直觀、數(shù)據(jù)精確和零件無(wú)損等特點(diǎn)[1].工業(yè)CT依靠放射性核素為發(fā)射源，因而會(huì)產(chǎn)生輻射，需在鉛制隔離房?jī)?nèi)檢測(cè)，進(jìn)而限制了使用環(huán)境，且軟、硬件成本較高[2-3].

隨著超聲波技術(shù)的不斷發(fā)展，在20世紀(jì)初出現(xiàn)了超聲波CT技術(shù)，受到眾多學(xué)者的廣泛關(guān)注，多位學(xué)者對(duì)超聲波CT技術(shù)及應(yīng)用開(kāi)展了相關(guān)研究.張蕾等[4-5]研究了小波包變換在超聲波CT特征提取中的應(yīng)用，實(shí)現(xiàn)了超聲波CT的特征提取，分析總結(jié)了超聲波CT技術(shù)的多種發(fā)展與應(yīng)用方向.劉帥鵬[6]使用超聲波CT技術(shù)對(duì)橋梁灌注樁進(jìn)行檢測(cè)，得到樁座的三維CT成像圖.鄒云等[7]通過(guò)激光超聲檢測(cè)方法測(cè)定了板材彈性模量.Mohd-Khairi等[8]將超聲波CT技術(shù)應(yīng)用于盒裝牛奶內(nèi)部異物檢測(cè)，構(gòu)建了斷層圖像重建算法.Abbaszadeh等[9]設(shè)計(jì)了鋼管輸送機(jī)的超聲波CT掃描系統(tǒng)，實(shí)現(xiàn)了管道超聲波CT斷層圖片獲取.目前已經(jīng)實(shí)現(xiàn)利用超聲波CT技術(shù)對(duì)典型零件斷層圖像進(jìn)行獲取，但超聲波CT圖像的特征分辨對(duì)檢測(cè)人員要求較高，需要具備一定理論知識(shí)和檢測(cè)經(jīng)驗(yàn).將超聲波CT圖像轉(zhuǎn)化為三維模型可以直觀表達(dá)零件形貌與內(nèi)部損傷情況.目前超聲波CT圖像三維重構(gòu)技術(shù)研究人員較少，研究方向存在空缺.基于超聲波CT斷層圖像三維重構(gòu)技術(shù)研究可以克服傳統(tǒng)CT技術(shù)的不足，降低零件可視化技術(shù)成本，拓展無(wú)損檢測(cè)人機(jī)交互渠道，對(duì)無(wú)損檢測(cè)技術(shù)與智能工業(yè)體系的發(fā)展具有一定意義和價(jià)值.

本文擬對(duì)超聲波CT技術(shù)檢測(cè)的人工預(yù)制缺陷零件斷層圖像進(jìn)行模型重構(gòu).通過(guò)噪音去除、異常信號(hào)調(diào)整、信號(hào)增強(qiáng)等方法對(duì)二維圖像進(jìn)行處理，建立輪廓提取算法處理圖像并優(yōu)化，得到可用于模型重構(gòu)的斷層圖像.建立圖像堆疊算法及模型重構(gòu)算法，對(duì)預(yù)制損傷零件的斷層圖像進(jìn)行模型重構(gòu).

1 模型重構(gòu)方法

超聲波CT模型重構(gòu)流程如圖1所示.通過(guò)超聲波CT圖像采集、二維圖像處理、三維模型建立實(shí)現(xiàn)零件模型重構(gòu).利用相共陣檢測(cè)方法進(jìn)行超聲波CT圖像采集；基于Matlab軟件構(gòu)建圖像二維處理算法，從灰度轉(zhuǎn)換、噪音處理、輪廓提取三個(gè)方面對(duì)圖像進(jìn)行處理；利用Matlab軟件構(gòu)建重構(gòu)算法來(lái)重構(gòu)損傷模型.

圖1 模型重構(gòu)流程圖

1.1 圖像采集原理

CT成像的實(shí)質(zhì)是對(duì)掃描得到的每個(gè)點(diǎn)的投影數(shù)據(jù)進(jìn)行計(jì)算得到各點(diǎn)的衰減系數(shù)，奧地利科學(xué)家拉東(Radon)提出的拉東變換與反變換實(shí)現(xiàn)了歷史上的首次CT圖像重建[10].這種方法同樣適用于超聲波CT成像，在性質(zhì)相同的材料中聲速變化不大，可以忽略超聲波的衍射和折射.基于拉東變換分析投影數(shù)據(jù)，對(duì)被測(cè)介質(zhì)參數(shù)(聲速、衰減系數(shù)等)與接收數(shù)據(jù)之間的線性關(guān)系進(jìn)行運(yùn)算分析[11].選擇二維線性超聲波作為發(fā)射源，通過(guò)拉東反變換對(duì)CT圖像進(jìn)行重建，如利用FBP(filtered back projection)算法、代數(shù)重建算法等重建被測(cè)介質(zhì)聲學(xué)參量(聲速、衰減系數(shù)等)的分布圖像[12].

圖2為超聲波CT成像系統(tǒng)，通過(guò)計(jì)算機(jī)發(fā)射指令控制超聲波發(fā)射器及支架的移動(dòng)，再采集發(fā)射器信號(hào)與接收器信號(hào)并轉(zhuǎn)入計(jì)算機(jī)進(jìn)行處理，得到零件超聲波CT斷層圖像.

圖2 超聲波CT成像系統(tǒng)

1.2 二維圖像處理方法

由于超聲波成像具有對(duì)比度差、噪聲污染高的固有特性，其圖像中組織特征模糊，噪聲干擾已成為影響超聲圖像質(zhì)量的關(guān)鍵因素之一[13].因此，基于Matlab軟件構(gòu)建自適應(yīng)濾波算法對(duì)圖像進(jìn)行降噪處理，并將圖像轉(zhuǎn)換為灰度圖像后提取出輪廓邊界.利用邊緣檢測(cè)算法對(duì)圖像信號(hào)輪廓進(jìn)行提取.目前常見(jiàn)的邊緣算子較多，主要有Robert、Sobel、Prewitt、LoG、Canny算子等，各種算子特點(diǎn)[14]如表1所示.

表1 算子特點(diǎn)比較

由于超聲波CT檢測(cè)得到的圖像邊界并不明確，特征位置灰度值存在梯度漸變且存在一定的低噪音.綜合考慮表1中各種算子的優(yōu)缺點(diǎn)，采用適用于低噪音與灰度漸變的Prewitt算子對(duì)圖像邊緣進(jìn)行提取.Prewitt算子原理是對(duì)圖像進(jìn)行閾值控制，將灰度值達(dá)到閾值的點(diǎn)均認(rèn)為是邊緣點(diǎn)[15]，即選取合適閾值K，若點(diǎn)G(i，j)≥K，則G(i，j)為邊緣圖像點(diǎn).

通過(guò)計(jì)算像素點(diǎn)上下左右臨點(diǎn)的灰度值插值判斷邊緣位置.在圖像空間上利用兩個(gè)不同方向模板(分別檢測(cè)水平邊緣與垂直邊緣)對(duì)圖像進(jìn)行鄰域卷積，其算子卷積核可以表達(dá)為

(1)

(2)

針對(duì)圖像中的像素點(diǎn)G(i，j)，Prewitt算子可以定義為

G(i)=|[f(i-1，j-1)+f(i-1，j)+f(i-1，j+1)]-[f(i+1，j-1)+f(i+1，j)+f(i+1，j+1)]|

(3)

G(j)=|[f(i-1，j+1)+f(i，j+1)+f(i+1，j+1)]-[f(i-1，j-1)+f(i，j-1)+f(i+1，j-1)]|

(4)

可見(jiàn)，Prewitt算子實(shí)際上是先進(jìn)行非歸一化的均值平滑，然后再進(jìn)行差分，這樣在進(jìn)行輪廓提取的同時(shí)也可以去除一定的噪音影響，實(shí)現(xiàn)對(duì)二維圖像輪廓的提取.

1.3 三維模型重構(gòu)方法

二維切片圖像的每張圖像均為獨(dú)立數(shù)據(jù)，要實(shí)現(xiàn)模型重構(gòu)首先需將二維數(shù)據(jù)進(jìn)行連通構(gòu)成三維數(shù)據(jù)組.因此，建立圖像堆疊算法對(duì)處理后的二維數(shù)據(jù)進(jìn)行排序并導(dǎo)入Matlab工作區(qū)，然后對(duì)數(shù)據(jù)進(jìn)行維度擴(kuò)張并確定圖像的Z軸位置，進(jìn)而實(shí)現(xiàn)圖像的三維堆疊[16].

移動(dòng)立方體算法的基本原理是在兩張圖像像素點(diǎn)的某個(gè)區(qū)域內(nèi)進(jìn)行采樣，若采樣點(diǎn)在x、y、z三個(gè)方向上的分布是均勻的，采樣間距分別為Δx、Δy、Δz，則數(shù)據(jù)可以用三維矩陣來(lái)表示[17].圖3為三維體素模型.

圖3 三維體素模型

每8個(gè)相鄰采樣點(diǎn)(體素的角點(diǎn))構(gòu)成一個(gè)立方體成為一個(gè)體素，由于采樣點(diǎn)是離散的，可以利用三線性插值法通過(guò)角點(diǎn)值計(jì)算體素內(nèi)的任意一點(diǎn)的數(shù)值[18].將區(qū)間點(diǎn)數(shù)值作為依據(jù)建立面片頂點(diǎn)與平面法線.移動(dòng)立方體算法構(gòu)建流程為：

1)將斷層圖像讀入內(nèi)存；

2)掃描兩張斷層圖像，分別利用兩張圖像的4個(gè)像素點(diǎn)創(chuàng)建一個(gè)立方體；

3)比較立方體頂點(diǎn)處的8個(gè)密度值與曲面常數(shù)并計(jì)算立方體的索引；

4)利用索引從預(yù)先計(jì)算的表格中查找邊列表；

5)通過(guò)線性插值找到曲面邊交點(diǎn)；

6)利用中心差計(jì)算每個(gè)立方體頂點(diǎn)的單位法線，對(duì)每個(gè)三角形頂點(diǎn)的法線進(jìn)行插值處理；

7)輸出三角形頂點(diǎn)和頂點(diǎn)法線.

2 超聲波CT圖像采集及處理

試驗(yàn)設(shè)備選用北極星辰挺桿C掃描成像系統(tǒng)，此系統(tǒng)具有較為完備的超聲波掃描功能，可以實(shí)現(xiàn)超聲波CT成像與3D層析掃查.

2.1 超聲波CT圖像采集

圖4為零件模型與采集結(jié)果.預(yù)制損傷零件尺寸為50 mm×50 mm×125 mm，其中人工預(yù)制圓孔缺陷直徑為1 mm、深度為5 mm，選擇試件中均勻分布的4個(gè)圓孔作為檢測(cè)對(duì)象，其位置分別為10 mm×10 mm、20 mm×20 mm、30 mm×30 mm、40 mm×40 mm處.圖4b為不同位置的掃描斷層圖像，可以準(zhǔn)確反映不同斷層的損傷形貌.圖4c、d為在頂部與正面C掃描圖像損傷在模型中的具體位置.

圖4 零件模型與采集結(jié)果

2.2 二維圖像處理及算法優(yōu)化

對(duì)采集到的超聲波CT圖像進(jìn)行預(yù)處理，先將圖像轉(zhuǎn)變?yōu)榛叶葓D像并使用自適應(yīng)中值濾波方法對(duì)圖像進(jìn)行降噪處理.去除圖像中由于外界噪音所引起的異常信號(hào)，借助Matlab軟件編寫(xiě)自適應(yīng)中值濾波算法.自適應(yīng)中值濾波器的基本原理[19]是比較一定領(lǐng)域內(nèi)像素值的大小，取出其中值作為該領(lǐng)域中心像素的新值，同時(shí)根據(jù)預(yù)設(shè)好的條件，動(dòng)態(tài)改變中值濾波器的窗口尺寸，以同時(shí)兼顧去噪聲作用和保護(hù)細(xì)節(jié)效果，相應(yīng)算法偽代碼可以表示為

算法1：輸入Sxy，Smax；

窗口像素值按大小排列；

計(jì)算Zmin，Zmax，Zmed，Zxy；

計(jì)算A1=Zmed-Zmin，A2=Zmed-Zmax，

B1=Zxy-Zmin，B2=Zxy-Zmax；

ifA1>0且A2<0

ifB1>0且B2<0

printfZxy

else

printfZxy-Zmed

else

增大Sxy

ifSxy>Smax

printfZxy

else

重復(fù)算法

其中：Zmin=min(Sxy)為窗口Sxy中的最小灰度值；Zmax=max(Sxy)為窗口Sxy中的最大灰度值；Zmed=med(Sxy)為窗口Sxy中的灰度中值；Zxy為坐標(biāo)(x，y)處的灰度值；Smax為Sxy允許的最大尺寸.

圖像降噪處理效果如圖5所示.完成圖像噪聲去除后，利用Prewitt算子進(jìn)行輪廓提取.根據(jù)超聲波信號(hào)特性可知，圖像轉(zhuǎn)為灰度圖像后缺陷位置灰度值不均勻，利用Prewitt算子運(yùn)算得到的圖像邊緣存在丟失，計(jì)算結(jié)果中存在雜點(diǎn)且不連貫，圖像斷裂較多，不能達(dá)到預(yù)期結(jié)果.因此，再利用腐蝕膨脹開(kāi)運(yùn)算方法對(duì)輪廓進(jìn)行進(jìn)一步處理.

圖5 圖像降噪效果

膨脹算法與腐蝕算法原理[20]相似，通俗來(lái)講腐蝕即是刪除對(duì)象邊界某些像素；而膨脹則是為圖像中的對(duì)象邊界添加像素.對(duì)腐蝕與膨脹圖像進(jìn)行交集運(yùn)算就可以得到輪廓信息.腐蝕膨脹開(kāi)運(yùn)算方法是先腐蝕后膨脹的運(yùn)算方法，可以實(shí)現(xiàn)圖像中的小像素去除，從而達(dá)到一定去噪效果.利用膨脹后的圖像與腐蝕圖像相減即可得到膨脹腐蝕處理后的邊界輪廓模型.輪廓提取算法邏輯可以表示為

算法2：讀取圖像；

定義圖像灰度值范圍；

設(shè)置灰度閾值T；

定義梯度算子卷積模板Prewittx，y；

對(duì)圖像像素點(diǎn)G(i，j)進(jìn)行運(yùn)算；

ifG(i，j)≥T

(i，j)為邊緣點(diǎn)

else

(i，j)為非邊緣點(diǎn)

輸出邊緣點(diǎn)圖像；%Prewitt算子運(yùn)算結(jié)束

設(shè)置灰度閾值P；

根據(jù)P將灰度圖像轉(zhuǎn)為二值圖像；

設(shè)置膨脹腐蝕結(jié)構(gòu)元素SE；

腐蝕運(yùn)算I；

膨脹運(yùn)算J；

膨脹與腐蝕圖像求差K=J-I；

輸出處理圖像；%圖像開(kāi)運(yùn)算

利用Matlab軟件編寫(xiě)Prewitt算子算法的處理結(jié)果如圖6a所示，通過(guò)圖像開(kāi)運(yùn)算方法優(yōu)化輪廓算法處理后的輪廓圖像如圖6b所示.可見(jiàn)，通過(guò)優(yōu)化算法處理的圖像無(wú)雜點(diǎn)、斷點(diǎn)等現(xiàn)象，圖像輪廓連續(xù)邊界清晰，可以用于三維模型重構(gòu)處理.

圖6 輪廓提取結(jié)果

2.3 模型重構(gòu)

對(duì)處理后的二維圖像進(jìn)行數(shù)據(jù)對(duì)齊與堆疊，添加z軸方向數(shù)據(jù)對(duì)圖像進(jìn)行排列，將二維數(shù)據(jù)轉(zhuǎn)化為三維數(shù)據(jù)并實(shí)現(xiàn)二維圖像之間的數(shù)據(jù)信息鏈接，其基本原理是對(duì)圖像進(jìn)行編號(hào)與讀取，采用圖像名作為標(biāo)記，以此標(biāo)記作為圖層信息并通過(guò)循環(huán)語(yǔ)句將多個(gè)數(shù)據(jù)堆疊在一起，具體數(shù)據(jù)堆疊算法可以表示為

算法3：讀取圖片位置與格式；

標(biāo)記文件名樣式；

length(names)=文件標(biāo)號(hào)最大值；

a=文件名編號(hào)；

IMAGES=定義一個(gè)全為零的空矩陣；

fork=1：length(names)

nm=[當(dāng)前編號(hào)的文件信息]

image=寫(xiě)入(nm)

IMAGES(：；：；k)=image

end

保存IMAGES.mat；

隨后對(duì)堆疊數(shù)據(jù)進(jìn)行體素分割，計(jì)算出各個(gè)體素中的像素點(diǎn)值，將各頂點(diǎn)與法線的相對(duì)位置連接為等值面，所有體素中等值面與法線信息構(gòu)成整個(gè)數(shù)據(jù)的三維模型.利用Matlab軟件構(gòu)建的模型重構(gòu)算法可以表示為

算法4：讀取數(shù)據(jù)文件(IMAGES.mat)；

for(對(duì)每一個(gè)物體)

{

掃描兩層鏈表結(jié)構(gòu)INCIDENCE；

for(對(duì)每一個(gè)體素單元)

{

查表找到三角面片分布情況；

根據(jù)平面方向和所處位置將三角面片加入INCIDENCE

}

初始化三角面片鏈表FaceList、頂點(diǎn)鏈表PointList和多邊形鏈表PolyList；

for(對(duì)INCIDENCE中的每一個(gè)平面)

{

清空用于合并的二維數(shù)組Merger；

for(對(duì)于該平面上的每個(gè)三角形或矩形)

{

查表找到三角形或矩形的邊對(duì)應(yīng)于Merger中的編號(hào)，寫(xiě)入Merger；

}

掃描Merger，將圖形劃分為凸多邊形，加入PolyList；

}

將PolyList中涉及到的頂點(diǎn)加入PointList，同時(shí)建立頂點(diǎn)的逆向索引；

for(PolyList中的每個(gè)凸多邊形)

{

檢查其邊界上(不含端點(diǎn))是否有點(diǎn)存在；

找到“T”型點(diǎn)，加入該多邊形，同時(shí)進(jìn)行標(biāo)記；

進(jìn)行三角劃分并將三角面片加入FaceList；

}

清除PolyList；

清除INCIDENCE；

將FaceList中的數(shù)據(jù)轉(zhuǎn)移到數(shù)組Face-Array中；

清除FaceList；

將PointList中的數(shù)據(jù)轉(zhuǎn)移到數(shù)組VertexArray中，同時(shí)進(jìn)行插值；

清除PointList；

}

對(duì)所有頂點(diǎn)計(jì)算其法矢量；

輸出三維模型；

圖7a為采用算法3進(jìn)行數(shù)據(jù)堆疊的結(jié)果，經(jīng)過(guò)數(shù)據(jù)堆疊后，二維圖像間的通道已經(jīng)建立，可以借助側(cè)視圖堆疊數(shù)據(jù)觀察到正視方向上的圖像特征.圖7b為采用算法4處理后的模型重構(gòu)結(jié)果，模型外部輪廓與內(nèi)部損傷情況清晰可見(jiàn).

圖7 三維重構(gòu)模型

2.4 模型數(shù)據(jù)對(duì)比

借助建模軟件進(jìn)行數(shù)據(jù)測(cè)量，利用二維視圖工具截取模型的中間斷面部位，結(jié)果如圖8所示(單位：mm).由圖8可見(jiàn)，斷面圖中的缺陷形貌、外形特征及缺陷位置與實(shí)際相符，模型外形尺寸與損傷部位的形位偏差均在1 mm之內(nèi)，即模型外部輪廓與損傷位置還原精度較高，但重構(gòu)模型在損傷形貌還原方面還不夠精確.

圖8 重構(gòu)模型斷面圖

由分析結(jié)果可知，超聲波CT技術(shù)檢測(cè)斷層圖像方法檢測(cè)成本較低、安全性高、檢測(cè)迅速，可以作為一種傳統(tǒng)CT圖像的替代方法實(shí)現(xiàn)模型重構(gòu).雖然重構(gòu)技術(shù)仍存在一定不足，但超聲波CT斷層圖像為模型重構(gòu)技術(shù)帶來(lái)新的發(fā)展方向和更大的應(yīng)用空間，對(duì)模型重構(gòu)算法的開(kāi)發(fā)與研究具有一定的學(xué)術(shù)價(jià)值與廣闊的應(yīng)用前景.

3 結(jié) 論

通過(guò)以上分析可以得到如下結(jié)論：

1)通過(guò)Matlab軟件構(gòu)建輪廓提取與三維重構(gòu)等算法，實(shí)現(xiàn)了超聲波CT斷層圖像的模型重構(gòu)，驗(yàn)證了本文提出的模型重構(gòu)方法的可行性.

2)Prewitt和Robert算子輪廓提取算法對(duì)超聲波CT圖像處理效果不理想，綜合分析后利用膨脹腐蝕算法對(duì)Prewitt算子進(jìn)行優(yōu)化，改進(jìn)后的輪廓提取算法可得到清晰完整的輪廓信息，從而實(shí)現(xiàn)超聲波CT圖像的特征信息提取.

3)利用構(gòu)建的模型重構(gòu)算法對(duì)超聲波CT圖像進(jìn)行模型重構(gòu)的結(jié)果顯示，零件外部和內(nèi)部形貌完整且尺寸準(zhǔn)確，外形尺寸與損傷位置偏差約為1 mm，但零件的損傷形貌重構(gòu)不夠精細(xì)，有待改進(jìn)與優(yōu)化.