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

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

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

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

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

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

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

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

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

1.1 圖像采集原理

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

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

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

1.2 二維圖像處理方法

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

表1 算子特點比較

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

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

(1)

(2)

針對圖像中的像素點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)

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

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

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

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

圖3 三維體素模型

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

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

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

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

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

5)通過線性插值找到曲面邊交點；

6)利用中心差計算每個立方體頂點的單位法線，對每個三角形頂點的法線進行插值處理；

7)輸出三角形頂點和頂點法線.

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

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

2.1 超聲波CT圖像采集

圖4為零件模型與采集結(jié)果.預(yù)制損傷零件尺寸為50 mm×50 mm×125 mm，其中人工預(yù)制圓孔缺陷直徑為1 mm、深度為5 mm，選擇試件中均勻分布的4個圓孔作為檢測對象，其位置分別為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)化

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

算法1：輸入Sxy，Smax；

窗口像素值按大小排列；

計算Zmin，Zmax，Zmed，Zxy；

計算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ù)超聲波信號特性可知，圖像轉(zhuǎn)為灰度圖像后缺陷位置灰度值不均勻，利用Prewitt算子運算得到的圖像邊緣存在丟失，計算結(jié)果中存在雜點且不連貫，圖像斷裂較多，不能達到預(yù)期結(jié)果.因此，再利用腐蝕膨脹開運算方法對輪廓進行進一步處理.

圖5 圖像降噪效果

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

算法2：讀取圖像；

定義圖像灰度值范圍；

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

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

對圖像像素點G(i，j)進行運算；

ifG(i，j)≥T

(i，j)為邊緣點

else

(i，j)為非邊緣點

輸出邊緣點圖像；%Prewitt算子運算結(jié)束

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

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

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

腐蝕運算I；

膨脹運算J；

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

輸出處理圖像；%圖像開運算

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

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

2.3 模型重構(gòu)

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

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

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

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

a=文件名編號；

IMAGES=定義一個全為零的空矩陣；

fork=1：length(names)

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

image=寫入(nm)

IMAGES(：；：；k)=image

end

保存IMAGES.mat；

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

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

for(對每一個物體)

{

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

for(對每一個體素單元)

{

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

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

}

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

for(對INCIDENCE中的每一個平面)

{

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

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

{

查表找到三角形或矩形的邊對應(yīng)于Merger中的編號，寫入Merger；

}

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

}

將PolyList中涉及到的頂點加入PointList，同時建立頂點的逆向索引；

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

{

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

找到“T”型點，加入該多邊形，同時進行標(biāo)記；

進行三角劃分并將三角面片加入FaceList；

}

清除PolyList；

清除INCIDENCE；

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

清除FaceList；

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

清除PointList；

}

對所有頂點計算其法矢量；

輸出三維模型；

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

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

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

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

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

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

3 結(jié) 論

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

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

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

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