基于贝叶斯准则的双侧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