大口径管道壁外部ct局部扫描成像方法
【技术领域】
[0001] 本发明属于CT扫描成像技术领域，涉及一种大口径管道壁外部CT局部扫描成像 方法。
【背景技术】
[0002] 在实际的工业生产中，管道管壁等工业部件在制造和使用过程中不可避免地存在 裂纹和缺陷，随着使用过程中承受交变载荷，裂纹和缺陷可能会在使用过程中突然引发构 件的疲劳断裂，其危害性极大。因此，及时检测管壁内部的裂纹和缺陷对防止灾难事故的发 生及减少经济损失具有重要的意义。
[0003] 现有技术中，X射线CT成像检测能够无损、精确、快速地重建物体的内部缺陷。但 在实际检测过程中，经常会出现下面的情况：待重建物体尺寸大于探测器尺寸而导致X射 线无法完全覆盖物体；只对部件外层（如管壁等）的裂纹或者腐蚀感兴趣；物体直径太大 或内部含有其他流动物质（如使用中的管道内部含有流动液体等其他物质），X射线无法有 效穿透物体的内部或受干扰。为了解决上述问题，通常将探测器在物体的两侧对称放置，只 扫描感兴趣的物体外部环形区域，从而引起了大尺寸物体的外部重建问题。外部CT重建问 题具有更短的扫描时间，更低的射线剂量，并且能够避免物体内部其他流动物质的影响，因 此具有较高的应用价值。但受射线束扇角的限制，探测器对称放置扫描的物体尺寸有限，并 且，传统的外部CT多采用物体转动的扫描方式，不适用于固定管道的扫描。管道的偏置扫 描多见于DR(Digital Radiography，数字化X射线摄影）中。
[0004] 外部CT(Computed Tomography，计算机断层成像）重建问题是一种投影数据截断 问题。由于投影数据的不完备性，传统的解析重建算法如FBP(Filtered back-projection， 滤波反投影）算法的重建结果存在严重的条形伪影，无法达到实际应用的要求。Frank Natterer发展了一种正则化方法并指出在理想的条件下，外部重建问题的二维逆Radon 变换有唯一解，但是当含有噪声时解是极其不稳定的。E. T. Quinto对外部CT重建问题的 SVD(Singular Value Decomposition,奇异值分解）方法作了深入研究，但由于投影系数矩 阵的奇异值可能很小，导致解的不适定性。因此上述方法很难用于实际含噪投影数据的图 像重建。
【发明内容】
[0005] 有鉴于此，本发明的目的在于提供一种大口径管道壁外部CT局部扫描成像方法， 该成像方法扫描过程易于机械实现，扫描速度快，并且能得到高质量的重建图像。
[0006] 为达到上述目的，本发明提供如下技术方案：
[0007] 大口径管道壁外部CT局部扫描成像方法，该方法包括以下步骤：
[0008] 步骤1)射线源和探测器设置在围绕待检测管道中心的圆形轨道上；探测器偏置 放置，使射线束能够覆盖以管道中心为圆心，以r为半径的圆盘外的管道环形区域；射线源 和探测器沿着圆形轨道作圆周运动进行扫描，获得待检测管道外部环形区域的投影数据；
[0009] 步骤2)根据投影数据重建管道外部环形区域的图像；
[0010] 步骤3)显示重建图像。
[0011] 进一步，扫描开始时以射线源与原点连线所在的直线为y轴，并且以射线源指向 原点的方向为正方向，X轴垂直于y轴，并且与y轴构成固定的笛卡尔坐标系〇~xy ;在射 线源和探测器绕待检测管道作圆周运动的过程中，以坐标原点0建立旋转的笛卡尔坐标系 〇-ξ η，η轴为扫描过程中射线源与原点连线所在的直线，并且以射线源指向原点的方向 为正方向，ξ轴垂直于η轴，并且与η轴构成右手笛卡尔坐标系，X轴与ξ轴的夹角为旋 转角度θ。
[0012] 进一步，所述步骤2)根据投影数据重建管道外部环形区域的图像，具体包括以下 步骤：
[0013] 步骤2-1)凸集投影（POCS);
[0014] 步骤2-2)总变差最小化（TVM);
[0015] 步骤2-3)利用区域尺度拟合（RSF)模型对子区域进行平均化修正。
[0016] 进一步，所述步骤2-1)凸集投影具体包括以下步骤：
[0017] 步骤2-1-1)采用下面的加型代数迭代公式对待重建图像进行重建：
[0018]
[0019] 设待重建图像的像素点总数为N个，f表示待重建的数字图像，位于（s，t)处的像 素灰度值表示为f Sit，则f可以表示成HXW的图像矩阵f = (fSit)，其中H为待重建图像的 行数，W为待重建图像的列数，将图像的像素点逐点排列成向量？ = 其中N =HXW，设经过待重建图像的扫描射线数为M条，将射线投影数据按射线逐条排列为向量 75 = [/W7w,,,：r_ , W = (Wl])为MXN维的投影系数矩阵，其中Wl]表示第j个点对第i条 射线投影数据的贡献率；
[0020] 其中，, _ = 1,2，一具_;)表示第k次迭代时图像数据向量歹的第j个分量， Nraunt表示重建算法的最大迭代次数，/f表示图像数据向量I的初始值的第j个分量，P1 表示第i条射线对应的投影数据，k为迭代次数，λ为松弛因子；
[0021] 步骤2-1-2)引入非负性限制，得到图像数据的校正值：
[0022]
[0023] 其中，表示经过非负校正后得到的图像数据向量I啲第j个分量，/7(~表 示按加型代数迭代公式经过M1次迭代后的图像数据向量/的第j个分量。
[0024] 进一步，所述步骤2-2)总变差最小化TVM，具体包括以下步骤：
[0025] 步骤2-2-1)将总变差最小化的梯度下降方向f(TVM _初始化为f (TVM _ = f (_， 将下降程度(1_初始化为Cirocs= I If I ;通过以下公式计算图像的总变差TV(f)及 在图像（S，t)处的逼近形式的偏导数W
[0028] 步骤2-2-2)总变差梯度下降法迭代按下述公式进行：
[0026]
[0027]
[0029]
[0030] 其中，Ngy表示总变差梯度下降法的迭代次数，TV(f)表示图像数据f的总变差， τ为很小的正常数，表示第h次迭代图像像素点（s，t)处的总变差梯度下降方向， 表示第I1次迭代后整幅图像每个像素点处的总变差梯度下降方向矩阵，即v(;|i 1彡s彡H，1彡t彡W，H为待重建图像的行数，W为待重建图像的列数，M · I I表示向量的 Frobenius范数，fs,t表示位于（s，t)处的像素点的灰度值，:^^表示位于（s-l，t)处的像 素点的灰度值，4^表示位于（s，t-l)处的像素点的灰度值，f s+lit表示位于（s+l，t)处的 像素点的灰度值，fs+litl表示位于（s+l，t-l)处的像素点的灰度值，f Sit+1表示位于（s，t+l) 处的像素点的灰度值，fs lit+1表示位于（s+1，t+1)处的像素点的灰度值，α为权系数；令f w =，判断是否达到总变差最小化中预设的迭代次数Ntvm，如果是，则跳转至下一步骤 S2-3)，否则跳转至步骤S2-1)。
[0031] 进一步，所述步骤利用区域尺度拟合模型对子区域进行平均化修正，具体包括以 下步骤：
[0032] 2-3-1)使用区域尺度拟合RSF活动轮廓模型提取待重建图像的边缘，通过求解以 下梯度流演化方程得到水平集函数： UiN 丄 丄 Ο?(5 乙 O Λ J ^ 吁/ O X
[0033]
[0034] 其中，X，y为表示图像中像素点位置的二维坐标向量，f (X)表示图像在X处的灰 度值，&00和fb(x)分别为图像X处的局部区域内，轮廓线的内部和外部的像素加权强度 均值，Uy)和f b(y)分别为图像y处的局部区域内，轮廓线的内部和外部的像素加权强度 均值，
F为Gauss核函数，〇为尺度参数，*表示卷积运算， Φ (X)为水平集函数，为水平集函数Φ (X)的梯度，He (ζ)为Heaviside函数的正则化 形式，S E (z)为维Dirac测度的正则化形式，ε为正常数，div(·)表示散度算子，F2表示 Laplace算子，λ λ 2>〇, μ，v>〇是各项的权值系数，t为引入的时间辅助变量；
[0035] 步骤2-3-2)在得到RSF活动轮廓模型的水平集函数后，利用水平集函数将图像划 分成不同的子区域，分别用每个子区域内的像素点灰度的平均值代替该子区域内各像素点 的灰度值；判定是否达到设定的迭代次数凡，如果是，则结束迭代，否则跳转至步骤S2-1)。
[0036]