基于貝葉斯準則的雙側CT圖像重建研究

張 彧 程 華* 袁 野

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

基于貝葉斯準則的雙側CT圖像重建研究

張 彧 程 華* 袁 野

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

對彈性波CT層析成像技術中增加檢測點陣列的數量的問題進行了研究，提出了基于貝葉斯準則的缺陷識別方法，指出該方法達到了提高反演圖像分辨率和缺陷判別的準確度的效果。

彈性波，反演算法，貝葉斯準則，圖像重建

0 引言

彈性波CT層析成像技術[1-3]的基本原理是由傳感器結構逐層提取射線信息，再通過算法重建二維、三維圖像。CT走時方程是高度病態性的大型稀疏矩陣，增加測量陣列，可以補充方程數量，提高求解精度。現有算法未綜合考慮統計信息而認為缺陷判別是確定性事件，而實際工程中，儀器精度、起跳點判讀、噪聲等因素都會導致結果不準，利用概率統計將缺陷判別作為不確定性事件處理更切實際。

1 基于彈性波CT走時補充方程的反演算法

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

1.1 彈性波CT反演算法

CT技術將構件待測斷面劃分成若干個單元，發射、接收換能器分別在對應的兩側傳遞信號，使每個網格都有測線通過。再將測得的走時數據反演重建得到圖像。彈性波的走時可用慢度s(r)(速度的倒數)沿射線的積分表示如下：

(1)

將參數離散化可得：

(2)

其中，Si為第i個單元內的平均慢度；lij為射線i在第j單元內的長度；m為離散單元數量。寫成緊湊形式為：

Ax=B

(3)

其中，A為射線路徑矩陣；x為慢度向量；B為走時向量。圖像重構的關鍵是求解慢度向量x。

SIRT法[5,6]是將測線走時誤差求和后依次平均到每個網格單元中，再對一個單元格中全部射線的修正值取平均，為該單元格內的修正值：

(4)

SIRT算法的步驟為：

(5)

3)設定某一極小正數e，當丨ti-Δti丨≤e或迭代一定次數后停止，否則進入下一輪迭代。

1.2 彈性波CT走時方程補充

傳統的測點布置方法[7]是在測區兩側分別設置激發點和接收點，以8×8網格為例，陣列布置如圖1所示。現將測點布置補充為如圖2所示。

2 基于貝葉斯準則的圖像重建算法

貝葉斯統計方法[8,9]將未知參數的先驗信息與樣本信息綜合，推得未知參數的后驗分布，再求得未知參數，實現對初始模型的修正，提高成像的分辨率。

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

B=AX+n

(6)

其中，A為m×n維的射線路徑矩陣;X為n×1維的慢度向量;B為m×1維的走時向量;n為m×1維的測量噪聲。假設隨機噪聲為標準正態分布，則其概率分布函數為：

(7)

從式(7)可得，統計逆問題不僅是對未知參數的估計，更是概率分布。CT統計逆問題就是獲得慢度向量x的后驗分布，據貝葉斯公式[9]可得：

(8)

變換得到:

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

(9)

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

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

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

(10)

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

據式(10)可推得后驗分布函數為：

(11)

本文采用最大后驗概率估計法對后驗分布進行處理，以此推斷單元網格中的慢度。

3 圖像重建

3.1 數值模型一

構造數值模型一如圖3所示，將待測區域劃分為8×8的網格，將速度異常網格標記為黑色，波速為3 000 m/s，其余正常區波速為4 500 m/s。對單方向陣列測量和本文測量方法進行對比。反演成像結果如圖4所示。

不同陣列方向數值模擬結果對比(一)見表1。

表1 不同陣列方向數值模擬結果對比(一)

3.2 數值模型二

構造數值模型二如圖5所示，將待測區域劃分為8×8的網格，其中白色正常區的波速值為4 500 m/s，黑色異常區波速為3 000 m/s，對原始方法和本文方法進行對比。反演成像結果如圖6所示。

不同陣列方向數值模擬結果對比(二)見表2。

表2 不同陣列方向數值模擬結果對比(二)

4 構件試驗

設計并制作兩個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所示。試驗采用DH5960超高速動態信號采集儀、KDL力錘以及帶磁座的加速度傳感器作為試驗設備。將待測構件劃分成8×8的網格，設置兩個方向，共計16對測點。限于篇幅，只列出試件一的實測走時如表3所示。

對圖7，圖8采用單方向陣列測量方法和本文的方法進行計算，重建結果如圖9，圖10所示。

圖9，圖10中的白色方框為預設缺陷，可以看出，本文所述的基于貝葉斯準則的雙側CT圖像重建方法相較于傳統方法而言成像效果更好，更吻合實際缺陷，且更好的抑制了偽像。

表3 試件一實測走時

5 結語

本文補充了彈性波CT的測點陣列，利用雙方向的射線增加走時方程的數量，改善了系數矩陣的病態性，最后利用貝葉斯準則對模型進行修正，加強圖像，提高判別準確度。數值模擬和試驗結果均顯示，本文方法在精度、抑制偽像方面均有一定效果，反演算法重建圖像的效果得到提高。

[1] 繆 侖.CT技術在混凝土超聲探傷中的應用[D].長沙：湖南大學,2001.

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

[3] 蘇林王,李平杰,肖永順,等.CT掃描技術在混凝土結構檢測中的應用[J].水運工程,2015(12):28-31.

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

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

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

[7] 王浩全,韓 焱,殷 黎.基于超聲CT的混凝土質量陣列檢測方法研究[J].計算機工程與應用,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].重慶：重慶大學,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