基于貝葉斯準(zhǔn)則的雙側(cè)CT圖像重建研究

張 彧 程 華* 袁 野

(中國(guó)人民解放軍后勤工程學(xué)院，重慶 401331)

基于貝葉斯準(zhǔn)則的雙側(cè)CT圖像重建研究

張 彧 程 華* 袁 野

(中國(guó)人民解放軍后勤工程學(xué)院，重慶 401331)

對(duì)彈性波CT層析成像技術(shù)中增加檢測(cè)點(diǎn)陣列的數(shù)量的問(wèn)題進(jìn)行了研究，提出了基于貝葉斯準(zhǔn)則的缺陷識(shí)別方法，指出該方法達(dá)到了提高反演圖像分辨率和缺陷判別的準(zhǔn)確度的效果。

彈性波，反演算法，貝葉斯準(zhǔn)則，圖像重建

0 引言

彈性波CT層析成像技術(shù)[1-3]的基本原理是由傳感器結(jié)構(gòu)逐層提取射線信息，再通過(guò)算法重建二維、三維圖像。CT走時(shí)方程是高度病態(tài)性的大型稀疏矩陣，增加測(cè)量陣列，可以補(bǔ)充方程數(shù)量，提高求解精度?，F(xiàn)有算法未綜合考慮統(tǒng)計(jì)信息而認(rèn)為缺陷判別是確定性事件，而實(shí)際工程中，儀器精度、起跳點(diǎn)判讀、噪聲等因素都會(huì)導(dǎo)致結(jié)果不準(zhǔn)，利用概率統(tǒng)計(jì)將缺陷判別作為不確定性事件處理更切實(shí)際。

1 基于彈性波CT走時(shí)補(bǔ)充方程的反演算法

目前方程組的求解多采用間接算法[4]，常見(jiàn)的有反投影法(BPT)、奇異值分解法(SVD)、代數(shù)重建法(ART)和聯(lián)合迭代重建法(SIRT)等。其中，ART算法的計(jì)算速度和效率都比較高，但是其對(duì)初始值的依賴(lài)程度更高，且超過(guò)最佳迭代次數(shù)后計(jì)算結(jié)果不收斂；BPT法[7]計(jì)算簡(jiǎn)單，速度較快，但是結(jié)果不能重復(fù)，所得圖像的分辨率較低，常用于求解算法當(dāng)中所需的慢度初始值；SIRT算法[8]收斂性好、對(duì)初值的要求不高且對(duì)誤差不敏感，因而應(yīng)用廣泛。本文以SIRT算法作為基礎(chǔ)，通過(guò)補(bǔ)充走時(shí)方程提高圖像重建的分辨率。

1.1 彈性波CT反演算法

CT技術(shù)將構(gòu)件待測(cè)斷面劃分成若干個(gè)單元，發(fā)射、接收換能器分別在對(duì)應(yīng)的兩側(cè)傳遞信號(hào)，使每個(gè)網(wǎng)格都有測(cè)線通過(guò)。再將測(cè)得的走時(shí)數(shù)據(jù)反演重建得到圖像。彈性波的走時(shí)可用慢度s(r)(速度的倒數(shù))沿射線的積分表示如下：

(1)

將參數(shù)離散化可得：

(2)

其中，Si為第i個(gè)單元內(nèi)的平均慢度；lij為射線i在第j單元內(nèi)的長(zhǎng)度；m為離散單元數(shù)量。寫(xiě)成緊湊形式為：

Ax=B

(3)

其中，A為射線路徑矩陣；x為慢度向量；B為走時(shí)向量。圖像重構(gòu)的關(guān)鍵是求解慢度向量x。

SIRT法[5,6]是將測(cè)線走時(shí)誤差求和后依次平均到每個(gè)網(wǎng)格單元中，再對(duì)一個(gè)單元格中全部射線的修正值取平均，為該單元格內(nèi)的修正值：

(4)

SIRT算法的步驟為：

(5)

3)設(shè)定某一極小正數(shù)e，當(dāng)丨ti-Δti丨≤e或迭代一定次數(shù)后停止，否則進(jìn)入下一輪迭代。

1.2 彈性波CT走時(shí)方程補(bǔ)充

傳統(tǒng)的測(cè)點(diǎn)布置方法[7]是在測(cè)區(qū)兩側(cè)分別設(shè)置激發(fā)點(diǎn)和接收點(diǎn)，以8×8網(wǎng)格為例，陣列布置如圖1所示。現(xiàn)將測(cè)點(diǎn)布置補(bǔ)充為如圖2所示。

2 基于貝葉斯準(zhǔn)則的圖像重建算法

貝葉斯統(tǒng)計(jì)方法[8,9]將未知參數(shù)的先驗(yàn)信息與樣本信息綜合，推得未知參數(shù)的后驗(yàn)分布，再求得未知參數(shù)，實(shí)現(xiàn)對(duì)初始模型的修正，提高成像的分辨率。

考慮噪聲影響，式(3)寫(xiě)為：

B=AX+n

(6)

其中，A為m×n維的射線路徑矩陣;X為n×1維的慢度向量;B為m×1維的走時(shí)向量;n為m×1維的測(cè)量噪聲。假設(shè)隨機(jī)噪聲為標(biāo)準(zhǔn)正態(tài)分布，則其概率分布函數(shù)為：

(7)

從式(7)可得，統(tǒng)計(jì)逆問(wèn)題不僅是對(duì)未知參數(shù)的估計(jì)，更是概率分布。CT統(tǒng)計(jì)逆問(wèn)題就是獲得慢度向量x的后驗(yàn)分布，據(jù)貝葉斯公式[9]可得：

(8)

變換得到:

π(X|B)=C·π(x)·π(B|x)

(9)

其中，C為邊緣密度函數(shù)的積分常數(shù)。

由于CT層析成像對(duì)于噪聲很敏感，難以保證解空間的穩(wěn)定性，因此將先驗(yàn)分布π(S)用于降低噪聲的干擾，平滑圖像，通常為以下形式：

π(S)=e-αA(s)

(10)

其中,A(s)為對(duì)應(yīng)具體的正則化形式；α為正則化參數(shù)，表示正則化的程度。

據(jù)式(10)可推得后驗(yàn)分布函數(shù)為：

(11)

本文采用最大后驗(yàn)概率估計(jì)法對(duì)后驗(yàn)分布進(jìn)行處理，以此推斷單元網(wǎng)格中的慢度。

3 圖像重建

3.1 數(shù)值模型一

構(gòu)造數(shù)值模型一如圖3所示，將待測(cè)區(qū)域劃分為8×8的網(wǎng)格，將速度異常網(wǎng)格標(biāo)記為黑色，波速為3 000 m/s，其余正常區(qū)波速為4 500 m/s。對(duì)單方向陣列測(cè)量和本文測(cè)量方法進(jìn)行對(duì)比。反演成像結(jié)果如圖4所示。

不同陣列方向數(shù)值模擬結(jié)果對(duì)比(一)見(jiàn)表1。

表1 不同陣列方向數(shù)值模擬結(jié)果對(duì)比(一)

3.2 數(shù)值模型二

構(gòu)造數(shù)值模型二如圖5所示，將待測(cè)區(qū)域劃分為8×8的網(wǎng)格，其中白色正常區(qū)的波速值為4 500 m/s，黑色異常區(qū)波速為3 000 m/s，對(duì)原始方法和本文方法進(jìn)行對(duì)比。反演成像結(jié)果如圖6所示。

不同陣列方向數(shù)值模擬結(jié)果對(duì)比(二)見(jiàn)表2。

表2 不同陣列方向數(shù)值模擬結(jié)果對(duì)比(二)

4 構(gòu)件試驗(yàn)

設(shè)計(jì)并制作兩個(gè)800 mm×800 mm×150 mm的混凝土試件，試件一的缺陷尺寸為200 mm×200 mm×150 mm；試件二的缺陷尺寸為200 mm×200 mm×150 mm,100 mm×100 mm×150 mm；分別如圖7，圖8所示。試驗(yàn)采用DH5960超高速動(dòng)態(tài)信號(hào)采集儀、KDL力錘以及帶磁座的加速度傳感器作為試驗(yàn)設(shè)備。將待測(cè)構(gòu)件劃分成8×8的網(wǎng)格，設(shè)置兩個(gè)方向，共計(jì)16對(duì)測(cè)點(diǎn)。限于篇幅，只列出試件一的實(shí)測(cè)走時(shí)如表3所示。

對(duì)圖7，圖8采用單方向陣列測(cè)量方法和本文的方法進(jìn)行計(jì)算，重建結(jié)果如圖9，圖10所示。

圖9，圖10中的白色方框?yàn)轭A(yù)設(shè)缺陷，可以看出，本文所述的基于貝葉斯準(zhǔn)則的雙側(cè)CT圖像重建方法相較于傳統(tǒng)方法而言成像效果更好，更吻合實(shí)際缺陷，且更好的抑制了偽像。

表3 試件一實(shí)測(cè)走時(shí)

5 結(jié)語(yǔ)

本文補(bǔ)充了彈性波CT的測(cè)點(diǎn)陣列，利用雙方向的射線增加走時(shí)方程的數(shù)量，改善了系數(shù)矩陣的病態(tài)性，最后利用貝葉斯準(zhǔn)則對(duì)模型進(jìn)行修正，加強(qiáng)圖像，提高判別準(zhǔn)確度。數(shù)值模擬和試驗(yàn)結(jié)果均顯示，本文方法在精度、抑制偽像方面均有一定效果，反演算法重建圖像的效果得到提高。

[1] 繆 侖.CT技術(shù)在混凝土超聲探傷中的應(yīng)用[D].長(zhǎng)沙：湖南大學(xué),2001.

[2] 顏 勇,林德宏.超聲波CT技術(shù)在某大橋樁基檢測(cè)中的應(yīng)用[J].建筑,2012(16):89-90.

[3] 蘇林王,李平杰,肖永順,等.CT掃描技術(shù)在混凝土結(jié)構(gòu)檢測(cè)中的應(yīng)用[J].水運(yùn)工程,2015(12):28-31.

[4] 黃 靚,黃政宇,汪 優(yōu).結(jié)構(gòu)混凝土超聲波層析成像的反演算法研究[J].湖南大學(xué)學(xué)報(bào)(自然科學(xué)版),2006(5):26-30.

[5] 牛法富,許令周.BPT算法SIRT算法的一種加權(quán)研究[J].青島大學(xué)學(xué)報(bào)(自然科學(xué)版),2011,24(2):28-32.

[6] 黃曉寒,王仲剛,高遠(yuǎn)富,等.基于主元加權(quán)和自控步長(zhǎng)的反演成像聯(lián)合迭代法[J].后勤工程學(xué)院學(xué)報(bào),2014,30(4):85-89.

[7] 王浩全,韓 焱,殷 黎.基于超聲CT的混凝土質(zhì)量陣列檢測(cè)方法研究[J].計(jì)算機(jī)工程與應(yīng)用,2010(15):239-241.

[8] Iman Elyasi, Mohammad Ali Pourmina. Reduction of speckle noise ultrasound images based on TV regularization and modified Bayes shrink techniques[J].Optik-International Journal for Light and Electron Optics,2016(4)：33-34.

[9] 張建新.基于貝葉斯方法的有限元模型修正研究[D].重慶：重慶大學(xué),2014.

ThemethodofbilateralCTimagereconstructionbasedonBayesianprinciple

ZhangYuChengHua*YuanYe

(LogisticalEngineeringUniversityofPLA，Chongqing401331，China)

Increasing the number of testing point array in the elastic wave CT technology was studied. The defects recognition method based on Bayesian principle was proposed. The accuracy of the inversion image resolution and defect discrimination was enhanced.

elastic wave, inversion algorithm, Bayesian principle, image reconstruction

1009-6825(2017)23-0040-03

2017-06-06

張 彧(1993- )，女，在讀碩士; 袁 野(1992- )，男，在讀碩士

程 華(1958- )，男，教授

O242

：A