CN105136823A - 大口径管道壁外部ct局部扫描成像方法 - Google Patents

Abstract

本发明涉及一种大口径管道壁外部CT局部扫描成像方法，属于CT扫描成像技术领域。该方法将射线源和探测器设置在围绕待检测管道中心的圆形轨道上，并将探测器偏置放置；射线源和探测器沿着圆形轨道作圆周运动进行扫描，获得待检测管道外部环形区域的投影数据；将TVM-POCS重建算法与区域尺度拟合分割方法相结合，根据投影数据重建管道外部环形区域的图像。本发明提供的一种大口径管道壁外部CT局部扫描成像方法，扫描方式简单易行，扫描时间短，射线剂量低；该方法能够很好地处理投影数据截断问题，重建伪影大大减轻，并能很好地处理由于射线束硬化引起的重建图像灰度不均的问题，最后显示的管道外部环形区域的局部重建图像质量较好。

Description

大口径管道壁外部CT局部扫描成像方法

技术领域

本发明属于CT扫描成像技术领域，涉及一种大口径管道壁外部CT局部扫描成像方法。

背景技术

在实际的工业生产中，管道管壁等工业部件在制造和使用过程中不可避免地存在裂纹和缺陷，随着使用过程中承受交变载荷，裂纹和缺陷可能会在使用过程中突然引发构件的疲劳断裂，其危害性极大。因此，及时检测管壁内部的裂纹和缺陷对防止灾难事故的发生及减少经济损失具有重要的意义。

现有技术中，X射线CT成像检测能够无损、精确、快速地重建物体的内部缺陷。但在实际检测过程中，经常会出现下面的情况：待重建物体尺寸大于探测器尺寸而导致X射线无法完全覆盖物体；只对部件外层(如管壁等)的裂纹或者腐蚀感兴趣；物体直径太大或内部含有其他流动物质(如使用中的管道内部含有流动液体等其他物质)，X射线无法有效穿透物体的内部或受干扰。为了解决上述问题，通常将探测器在物体的两侧对称放置，只扫描感兴趣的物体外部环形区域，从而引起了大尺寸物体的外部重建问题。外部CT重建问题具有更短的扫描时间，更低的射线剂量，并且能够避免物体内部其他流动物质的影响，因此具有较高的应用价值。但受射线束扇角的限制，探测器对称放置扫描的物体尺寸有限，并且，传统的外部CT多采用物体转动的扫描方式，不适用于固定管道的扫描。管道的偏置扫描多见于DR(DigitalRadiography，数字化X射线摄影)中。

外部CT(ComputedTomography，计算机断层成像)重建问题是一种投影数据截断问题。由于投影数据的不完备性，传统的解析重建算法如FBP(Filteredback-projection，滤波反投影)算法的重建结果存在严重的条形伪影，无法达到实际应用的要求。FrankNatterer发展了一种正则化方法并指出在理想的条件下，外部重建问题的二维逆Radon变换有唯一解，但是当含有噪声时解是极其不稳定的。E.T.Quinto对外部CT重建问题的SVD(SingularValueDecomposition,奇异值分解)方法作了深入研究，但由于投影系数矩阵的奇异值可能很小，导致解的不适定性。因此上述方法很难用于实际含噪投影数据的图像重建。

发明内容

有鉴于此，本发明的目的在于提供一种大口径管道壁外部CT局部扫描成像方法，该成像方法扫描过程易于机械实现，扫描速度快，并且能得到高质量的重建图像。

为达到上述目的，本发明提供如下技术方案：

大口径管道壁外部CT局部扫描成像方法，该方法包括以下步骤：

步骤1)射线源和探测器设置在围绕待检测管道中心的圆形轨道上；探测器偏置放置，使射线束能够覆盖以管道中心为圆心，以r为半径的圆盘外的管道环形区域；射线源和探测器沿着圆形轨道作圆周运动进行扫描，获得待检测管道外部环形区域的投影数据；

步骤2)根据投影数据重建管道外部环形区域的图像；

步骤3)显示重建图像。

进一步，扫描开始时以射线源与原点连线所在的直线为y轴，并且以射线源指向原点的方向为正方向，x轴垂直于y轴，并且与y轴构成固定的笛卡尔坐标系O-xy；在射线源和探测器绕待检测管道作圆周运动的过程中，以坐标原点O建立旋转的笛卡尔坐标系O-ξη，η轴为扫描过程中射线源与原点连线所在的直线，并且以射线源指向原点的方向为正方向，ξ轴垂直于η轴，并且与η轴构成右手笛卡尔坐标系，x轴与ξ轴的夹角为旋转角度θ。

进一步，所述步骤2)根据投影数据重建管道外部环形区域的图像，具体包括以下步骤：

步骤2-1)凸集投影(POCS)；

步骤2-2)总变差最小化(TVM)；

步骤2-3)利用区域尺度拟合(RSF)模型对子区域进行平均化修正。

进一步，所述步骤2-1)凸集投影具体包括以下步骤：

步骤2-1-1)采用下面的加型代数迭代公式对待重建图像进行重建：

$\left\{\begin{array}{cc}{f}_j^{\left(0\right)}=0& \\ f_j^{\left(k\right)}=f_j^{(k-1)}+\mathrm{\λ}\frac{p_i-\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{w}_{ij}\·f_j^{(k-1)}}{\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{w}_{ij}}{w}_{ij}& \left(j=1,2,...,N;i=1,2,...,M\right)\end{array}\right.$

设待重建图像的像素点总数为N个，f表示待重建的数字图像，位于(s,t)处的像素灰度值表示为f_s,t，则f可以表示成H×W的图像矩阵f＝(f_s,t)，其中H为待重建图像的行数，W为待重建图像的列数，将图像的像素点逐点排列成向量其中N＝H×W，设经过待重建图像的扫描射线数为M条，将射线投影数据按射线逐条排列为向量W＝(w_ij)为M×N维的投影系数矩阵，其中w_ij表示第j个点对第i条射线投影数据的贡献率；

其中，表示第k次迭代时图像数据向量的第j个分量，N_count表示重建算法的最大迭代次数，表示图像数据向量的初始值的第j个分量，p_i表示第i条射线对应的投影数据，k为迭代次数，λ为松弛因子；

步骤2-1-2)引入非负性限制，得到图像数据的校正值：

$f_j^{\left(POCS\right)}=\left\{\begin{array}{cc}{f}_j^{\left(M_1\right)},& \begin{array}{cc}if& f_j^{\left(M_1\right)}>0\end{array}\\ 0,& else\end{array}\right.(j=1,2,...,N)$

其中，表示经过非负校正后得到的图像数据向量的第j个分量，表示按加型代数迭代公式经过M₁次迭代后的图像数据向量的第j个分量。

进一步，所述步骤2-2)总变差最小化TVM，具体包括以下步骤：

步骤2-2-1)将总变差最小化的梯度下降方向f^(TVM-GRAD)初始化为f^(TVM-GRAD)＝f^(POCS)，将下降程度d_POCS初始化为d_POCS＝||f⁽⁰⁾-f^(POCS)||；通过以下公式计算图像的总变差TV(f)及在图像(s,t)处的逼近形式的偏导数v_s,t：

$TV\left(f\right)=\underset{s,t}{\Σ}\sqrt{{(f_{s,t}-f_{s-1,t})}^2+{(f_{s,t}-f_{s,t-1})}^2+\mathrm{\τ}}$

$\begin{array}{c}{v}_{s,t}=\frac{\∂TV\left(f\right)}{\∂f_{s,t}}\≈\frac{\left(f_{s,t}-f_{s-1,t}\right)+\left(f_{s,t}-f_{s,t-1}\right)}{\sqrt{{\left(f_{s,t}-f_{s-1,t}\right)}^2+{\left(f_{s,t}-f_{s,t-1}\right)}^2+\mathrm{\τ}}}\\ -\frac{f_{s+1,t}-f_{s,t}}{\sqrt{{\left(f_{s+1,t}-f_{s,t}\right)}^2+{\left(f_{s+1,t}-f_{s+1,t-1}\right)}^2+\mathrm{\τ}}}\\ -\frac{f_{s,t+1}-f_{s,t}}{\sqrt{{\left(f_{s,t+1}-f_{s,t}\right)}^2+{\left(f_{s,t+1}-f_{s-1,t+1}\right)}^2+\mathrm{\τ}}}\end{array}$

步骤2-2-2)总变差梯度下降法迭代按下述公式进行：

$\left\{\begin{array}{c}{v}_{s,t}^{\left(i_1\right)}=\frac{\∂TV\left(f^{(TV-GRAD)}\right)}{\∂f_{s,t}^{(TV-GRAD)}}\\ f^{(TV-GRAD)}=f^{(TV-GRAD)}-{\mathrm{\αd}}_{POCS}\frac{v^{\left(i_1\right)}}{\left|\right|v^{\left(i_1\right)}\left|\right|}\end{array}\right.i_1=1,2,...,N_{gard}$

其中，N_gard表示总变差梯度下降法的迭代次数，TV(f)表示图像数据f的总变差，τ为很小的正常数，表示第i₁次迭代图像像素点(s,t)处的总变差梯度下降方向，表示第i₁次迭代后整幅图像每个像素点处的总变差梯度下降方向矩阵，即1≤s≤H,1≤t≤W，H为待重建图像的行数，W为待重建图像的列数，||·||表示向量的Frobenius范数，f_s,t表示位于(s,t)处的像素点的灰度值，f_s-1,t表示位于(s-1,t)处的像素点的灰度值，f_s,t-1表示位于(s,t-1)处的像素点的灰度值，f_s+1,t表示位于(s+1,t)处的像素点的灰度值，f_s+1,t-1表示位于(s+1,t-1)处的像素点的灰度值，f_s,t+1表示位于(s,t+1)处的像素点的灰度值，f_s-1,t+1表示位于(s+1,t+1)处的像素点的灰度值，α为权系数；令f⁽⁰⁾＝f^(TVM-GRAD)，判断是否达到总变差最小化中预设的迭代次数N_TVM，如果是，则跳转至下一步骤S2-3)，否则跳转至步骤S2-1)。

进一步，所述步骤利用区域尺度拟合模型对子区域进行平均化修正，具体包括以下步骤：

2-3-1)使用区域尺度拟合RSF活动轮廓模型提取待重建图像的边缘，通过求解以下梯度流演化方程得到水平集函数：

$\left\{\begin{array}{c}{f}_o\left(x\right)=\frac{K_{\mathrm{\σ}}\left(x\right)*\[H_{\mathrm{\ϵ}}\left(\mathrm{\φ}\right(x\left)\right)f\left(x\right)\]}{K_{\mathrm{\σ}}\left(x\right)*H_{\mathrm{\ϵ}}\left(\mathrm{\φ}\right(x\left)\right)};\\ f_b\left(x\right)=\frac{K_{\mathrm{\σ}}\left(x\right)*\[(1-H_{\mathrm{\ϵ}}(\mathrm{\φ}\left(x\right))\left)f\right(x)\]}{K_{\mathrm{\σ}}\left(x\right)*(1-H_{\mathrm{\ϵ}}(\mathrm{\φ}\left(x\right)))};\\ \frac{\∂\mathrm{\φ}}{\∂t}=-{\mathrm{\δ}}_{\mathrm{\ϵ}}\left(\mathrm{\φ}\right)\left({\mathrm{\λ}}_1{e}_1\right(x)-{\mathrm{\λ}}_2{e}_2(x\left)\right)+{\mathrm{\μ\δ}}_{\mathrm{\ϵ}}\left(\mathrm{\φ}\right)div\left(\frac{\▿\mathrm{\φ}}{\left|\▿\mathrm{\φ}\right|}\right)+v({\▿}^2\mathrm{\φ}-div(\frac{\▿\mathrm{\φ}}{\left|\▿\mathrm{\φ}\right|}));\\ e_1\left(x\right)=\∫K_{\mathrm{\σ}}(y-x){|f\left(x\right)-f_o\left(y\right)|}^{2}dy;e_2\left(x\right)=\∫K_{\mathrm{\σ}}(y-x){|f\left(x\right)-f_b\left(y\right)|}^{2}dy;\\ {\mathrm{\δ}}_{\mathrm{\ϵ}}\left(z\right)=\frac{1}{\mathrm{\π}}\·\frac{\mathrm{\ϵ}}{{\mathrm{\ϵ}}^2+z^2};\\ H_{\mathrm{\ϵ}}\left(z\right)=\frac{1}{2}\[1+\frac{2}{\mathrm{\π}}\arctan \left(\frac{z}{\mathrm{\ϵ}}\right)\];\end{array}\right.$

其中，x,y为表示图像中像素点位置的二维坐标向量，f(x)表示图像在x处的灰度值，f_o(x)和f_b(x)分别为图像x处的局部区域内，轮廓线的内部和外部的像素加权强度均值，f_o(y)和f_b(y)分别为图像y处的局部区域内，轮廓线的内部和外部的像素加权强度均值，为Gauss核函数，σ为尺度参数，*表示卷积运算，φ(x)为水平集函数，为水平集函数φ(x)的梯度，H_ε(z)为Heaviside函数的正则化形式，δ_ε(z)为维Dirac测度的正则化形式，ε为正常数，div(·)表示散度算子，表示Laplace算子，λ₁,λ₂>0,μ,v>0是各项的权值系数，t为引入的时间辅助变量；

步骤2-3-2)在得到RSF活动轮廓模型的水平集函数后，利用水平集函数将图像划分成不同的子区域，分别用每个子区域内的像素点灰度的平均值代替该子区域内各像素点的灰度值；判定是否达到设定的迭代次数N_e，如果是,则结束迭代，否则跳转至步骤S2-1)。

本发明的有益效果在于：本发明提供的一种大口径管道壁外部CT局部扫描成像方法，利用小尺寸的探测器和现有的CT机的机构实现大口径管道壁外部环形区域内裂纹和缺陷的检测，扫描前只需将射线源和探测器设置在围绕管道的圆形轨道上，并使探测器偏置放置，射线源和探测器绕待检测管道作圆周旋转，得到管道外部环形区域的投影数据，扫描过程易于机械实现。然后通过将TVM-POCS(TotalVariationMinimization-ProjectionontoConvexSets,总变差最小化-投影到凸集)重建算法与RSF(Region-ScalableFitting,区域尺度拟合)分割方法相结合，得到待检测管道横截面外部环形区域的重建图像；该方法能够很好地处理投影数据截断问题，重建伪影大大减轻，能很好地处理由于射线束硬化引起的重建图像灰度不均的问题，并且重建速度较快，重建图像质量好，分辨率高；可用较低能量的射线及较小尺寸的探测器扫描重建较大直径的管道外部环形区域的图像；可用于对石油化工等行业的大口径管道(含固定管道)进行CT扫描成像并重建管壁的外部环形区域图像。

附图说明

为了使本发明的目的、技术方案和优点更加清楚，下面将结合附图对本发明作进一步的详细描述，其中：

图1为本发明所述方法的流程图；

图2为大口径管道壁外部检测结构示意图；

图3为CT扫描的坐标系统示意图；

图4为外部CT扫描范围示意图。

具体实施方式

下面将结合附图，对本发明的优选实施例进行详细的描述。

本发明提供的大口径管道壁外部CT局部扫描成像方法，如图1所示，具体包括以下步骤：

步骤1)射线源1和探测器2设置在围绕待检测管道4中心的圆形轨道3上，扫描过程中，将探测器偏置放置，使射线束能够覆盖以管道中心为圆心，以r为半径的圆盘外的管道环形区域，如图4所示，使射线源和探测器沿着围绕待检测管道的圆形轨道作圆周运动，形成以旋转中心为圆心的圆周轨迹，通过扫描和数据采集获得管道外部环形区域的扫描数据，用于重建管道外部环形区域的图像。

射线源发出的射线束不能完全覆盖管道，偏置后的探测器只能得到管道外部环形区域的投影数据，探测器得到的投影数据是射线束完全覆盖管道时得到的投影数据的一部分。

扫描开始时以射线源与原点连线所在的直线为y轴，并且以射线源指向原点的方向为正方向，x轴垂直于y轴，并且与y轴构成固定的笛卡尔坐标系O-xy(x轴由y轴绕原点顺时针旋转90度得到，如图3所示)，在射线源和探测器绕待检测管道作圆周运动的过程中，以坐标原点O建立旋转的笛卡尔坐标系O-ξη，η轴为扫描过程中射线源与原点连线所在的直线，并且以射线源指向原点的方向为正方向，ξ轴垂直于η轴，并且与η轴构成右手笛卡尔坐标系(ξ轴由η轴绕原点顺时针旋转90度得到)，x轴与ξ轴的夹角为旋转角度θ(如图3所示)。

步骤2)根据投影数据重建管道外部环形区域的图像；

假设管道感兴趣区域(RegionofInterest，ROI)为：

$ROI=\left\{(x,y)\right|r^2<x^2+y^2<r_{max}^2\}$

其中r,r_max分别为管道外部环形区域的内径和外径。在感兴趣区域外，重建图像的值为零。根据公式计算可以得到真实的投影数据P(θ,ξ)，其中θ为旋转角度，ξ为旋转坐标系O-ξη下的横坐标值，P(θ,ξ)为旋转角度为θ时，在旋转坐标O-ξη下横坐标为ξ处的投影数据值，为旋转角度为θ，在旋转坐标O-ξη下横坐标为ξ处的入射射线未衰减的辐射强度，I^θ,ξ为旋转角度为θ，在旋转坐标O-ξη下横坐标为ξ处的入射射线穿过管道衰减后的辐射强度，和I^θ,ξ均可以通过测量得到。将f(x,y)和P(θ,ξ)离散化得到向量和其中f(x,y)表示在固定坐标系O-xy下(x,y)处的图像灰度值，N表示待重建图像中像素点的总个数，M表示扫描射线的总条数。CT系统可以表示成离散-离散的形式：其中W＝(w_ij)为M×N维的投影系数矩阵，w_ij表示第j个点对第i条射线投影数据的贡献率。重建管道外部环形区域的图像，具体包括以下步骤：

步骤2-1)凸集投影POCS(ProjectionontoConvexSets,POCS)；

步骤2-1-1)采用下面的加型代数迭代公式对待重建图像进行重建：

$\left\{\begin{array}{cc}{f}_j^{\left(0\right)}=0& \\ f_j^{\left(k\right)}=f_j^{(k-1)}+\mathrm{\&lambda;}\frac{p_i-\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}{w}_{ij}\&CenterDot;f_j^{(k-1)}}{\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}{w}_{ij}}{w}_{ij}& \left(j=1,2,...,N;i=1,2,...,M\right)\end{array}\right.$

设待重建图像的像素点总数为N个，f表示待重建的数字图像，位于(s,t)处的像素灰度值表示为f_s,t，则f可以表示成H×W的图像矩阵f＝(f_s,t)，其中H为待重建图像的行数，W为待重建图像的列数，将图像的像素点逐点排列成向量其中N＝H×W，设经过待重建图像的扫描射线数为M条，将射线投影数据按射线逐条排列为向量W＝(w_ij)为M×N维的投影系数矩阵，其中w_ij表示第j个点对第i条射线投影数据的贡献率；

其中，表示第k次迭代时图像数据向量的第j个分量，N_count表示重建算法的最大迭代次数，例如可取N_count＝300。表示图像数据向量的初始值的第j个分量，p_i表示第i条射线对应的投影数据，k为迭代次数，λ为松弛因子，例如可取λ＝1；

步骤2-1-2)引入非负性限制，得到图像数据的校正值：

$f_j^{\left(POCS\right)}=\left\{\begin{array}{cc}{f}_j^{\left(M_1\right)},& \begin{array}{cc}if& f_j^{\left(M_1\right)}>0\end{array}\\ 0,& else\end{array}\right.(j=1,2,...,N)$

其中，表示经过非负校正后得到的图像数据向量的第j个分量，即经过凸集投影步骤后得到的图像数据向量的第j个分量，表示按加型代数迭代公式经过M₁次迭代后的图像数据向量的第j个分量。

步骤2-2)总变差最小化TVM(TotalVariationMinimization,TVM)；

步骤2-2-1)将总变差最小化的梯度下降方向f^(TVM-GRAD)初始化为f^(TVM-GRAD)＝f^(POCS)，将下降程度d_POCS初始化为d_POCS＝||f⁽⁰⁾-f^(POCS)||；然后计算图像的总变差TV(f)及在图像(s,t)处的逼近形式的偏导数v_s,t如下：

$TV\left(f\right)=\underset{s,t}{\&Sigma;}\sqrt{{(f_{s,t}-f_{s-1,t})}^2+{(f_{s,t}-f_{s,t-1})}^2+\mathrm{\&tau;}}$

$\begin{array}{c}{v}_{s,t}=\frac{\&part;TV\left(f\right)}{\&part;f_{s,t}}\&ap;\frac{\left(f_{s,t}-f_{s-1,t}\right)+\left(f_{s,t}-f_{s,t-1}\right)}{\sqrt{{\left(f_{s,t}-f_{s-1,t}\right)}^2+{\left(f_{s,t}-f_{s,t-1}\right)}^2+\mathrm{\&tau;}}}\\ -\frac{f_{s+1,t}-f_{s,t}}{\sqrt{{\left(f_{s+1,t}-f_{s,t}\right)}^2+{\left(f_{s+1,t}-f_{s+1,t-1}\right)}^2+\mathrm{\&tau;}}}\\ -\frac{f_{s,t+1}-f_{s,t}}{\sqrt{{\left(f_{s,t+1}-f_{s,t}\right)}^2+{\left(f_{s,t+1}-f_{s-1,t+1}\right)}^2+\mathrm{\&tau;}}}\end{array}$

步骤2-2-2)总变差梯度下降法迭代按下述公式进行：

$\left\{\begin{array}{c}{v}_{s,t}^{\left(i_1\right)}=\frac{\&part;TV\left(f^{(TV-GRAD)}\right)}{\&part;f_{s,t}^{(TV-GRAD)}}\\ f^{(TV-GRAD)}=f^{(TV-GRAD)}-{\mathrm{\&alpha;d}}_{POCS}\frac{v^{\left(i_1\right)}}{\left|\right|v^{\left(i_1\right)}\left|\right|}\end{array}\right.i_1=1,2,...,N_{gard}$

其中，N_gard表示总变差梯度下降法的迭代次数，例如可取N_gard＝20，TV(f)表示图像数据f的总变差，τ为很小的正常数，例如可取τ＝0.00000001，表示第i₁次迭代图像像素点(s,t)处的总变差梯度下降方向，表示第i₁次迭代后整幅图像每个像素点处的总变差梯度下降方向矩阵，即H为待重建图像的行数，W为待重建图像的列数，||·||表示向量的Frobenius范数，f_s,t表示位于(s,t)处的像素点的灰度值，f_s-1,t表示位于(s-1,t)处的像素点的灰度值，f_s,t-1表示位于(s,t-1)处的像素点的灰度值，f_s+1,t表示位于(s+1,t)处的像素点的灰度值，f_s+1,t-1表示位于(s+1,t-1)处的像素点的灰度值，f_s,t+1表示位于(s,t+1)处的像素点的灰度值，f_s-1,t+1表示位于(s+1,t+1)处的像素点的灰度值，α为权系数，例如可取α＝0.2；令f⁽⁰⁾＝f^(TVM-GRAD)，判断是否达到总变差最小化中预设的迭代次数N_TVM，例如可取N_TVM＝50，如果是，则跳转至下一步骤S2-3)，否则跳转至步骤S2-1)。

步骤2-3)利用区域尺度拟合(Region-ScalableFitting，RSF)模型对子区域进行平均化修正

2-3-1)使用RSF活动轮廓模型提取待重建图像的边缘，通过求解以下梯度流演化方程得到水平集函数：

$\left\{\begin{array}{c}{f}_o\left(x\right)=\frac{K_{\mathrm{\&sigma;}}\left(x\right)*\&lsqb;H_{\mathrm{\&epsiv;}}\left(\mathrm{\&phi;}\right(x\left)\right)f\left(x\right)\&rsqb;}{K_{\mathrm{\&sigma;}}\left(x\right)*H_{\mathrm{\&epsiv;}}\left(\mathrm{\&phi;}\right(x\left)\right)};\\ f_b\left(x\right)=\frac{K_{\mathrm{\&sigma;}}\left(x\right)*\&lsqb;(1-H_{\mathrm{\&epsiv;}}(\mathrm{\&phi;}\left(x\right))\left)f\right(x)\&rsqb;}{K_{\mathrm{\&sigma;}}\left(x\right)*(1-H_{\mathrm{\&epsiv;}}(\mathrm{\&phi;}\left(x\right)))};\\ \frac{\&part;\mathrm{\&phi;}}{\&part;t}=-{\mathrm{\&delta;}}_{\mathrm{\&epsiv;}}\left(\mathrm{\&phi;}\right)\left({\mathrm{\&lambda;}}_1{e}_1\right(x)-{\mathrm{\&lambda;}}_2{e}_2\left(x\right))+{\mathrm{\&mu;\&delta;}}_{\mathrm{\&epsiv;}}\left(\mathrm{\&phi;}\right)div\left(\frac{\&dtri;\mathrm{\&phi;}}{|\&dtri;\mathrm{\&phi;}|}\right)+v({\&dtri;}^2\mathrm{\&phi;}-div(\frac{\&dtri;\mathrm{\&phi;}}{|\&dtri;\mathrm{\&phi;}|}));\\ e_1\left(x\right)=\&Integral;K_{\mathrm{\&sigma;}}(y-x)|f\left(x\right)-f_o\left(y\right){|}^{2}dy;e_2\left(x\right)=\&Integral;K_{\mathrm{\&sigma;}}(y-x)|f\left(x\right)-f_b\left(y\right){|}^{2}dy;\\ {\mathrm{\&delta;}}_{\mathrm{\&epsiv;}}\left(z\right)=\frac{1}{\mathrm{\&pi;}}\&CenterDot;\frac{\mathrm{\&epsiv;}}{{\mathrm{\&epsiv;}}^2+z^2};\\ H_{\mathrm{\&epsiv;}}\left(z\right)=\frac{1}{2}\&lsqb;1+\frac{2}{\mathrm{\&pi;}}\arctan \left(\frac{z}{\mathrm{\&epsiv;}}\right)\&rsqb;;\end{array}\right.$

其中，x,y为表示图像中像素点位置的二维坐标向量，例如可取x＝(s,t),1≤s≤H,1≤t≤W，f(x)表示图像在x处的灰度值，f_o(x)和f_b(x)分别为图像x处的局部区域内，轮廓线的内部和外部的像素加权强度均值，f_o(y)和f_b(y)分别为图像y处的局部区域内，轮廓线的内部和外部的像素加权强度均值，为Gauss核函数，σ为尺度参数，例如可取σ＝3.0，*表示卷积运算，φ(x)为水平集函数，为水平集函数φ(x)的梯度，H_ε(z)为Heaviside函数的正则化形式，δ_ε(z)为Dirac测度的正则化形式，ε为正常数，例如可取ε＝1，div(·)表示散度算子，表示Laplace算子，λ₁,λ₂>0,μ,v>0是各项的权值系数，例如可取λ₁＝λ₂＝1，μ＝0.003×255×255,ν＝1.0，t为引入的时间辅助变量；

步骤2-3-2)在得到RSF活动轮廓模型的水平集函数后，利用这个水平集函数将图像划分成不同的子区域，分别用每个子区域内的像素点灰度的平均值代替该子区域内各像素点的灰度值；判定是否达到设定的迭代次数N_e，如果是,则结束迭代，否则跳转至步骤S2-1)。

步骤3)显示重建图像。

本发明将TVM-POCS(TotalVariationMinimization-ProjectionontoConvexSets,总变差最小化-投影到凸集)重建算法与RSF(Region-ScalableFitting,区域尺度拟合)分割模型结合，实现大口径管道外部环形区域的重建。在TVM-POCS重建算法的基础上引入RSF模型对重建的中间结果进行分割，能够大大地减轻图像边缘处的重建伪影，并且能够很好地处理由于射线束硬化引起的重建图像的灰度不均匀性，得到的重建图像质量高，分辨率高。

最后说明的是，以上优选实施例仅用以说明本发明的技术方案而非限制，尽管通过上述优选实施例已经对本发明进行了详细的描述，但本领域技术人员应当理解，可以在形式上和细节上对其作出各种各样的改变，而不偏离本发明权利要求书所限定的范围。

Claims (6)

1.大口径管道壁外部CT局部扫描成像方法，其特征在于：该方法包括以下步骤：

步骤1)射线源和探测器设置在围绕待检测管道中心的圆形轨道上；探测器偏置放置，使射线束能够覆盖以管道中心为圆心，以r为半径的圆盘外的管道环形区域；射线源和探测器沿着圆形轨道作圆周运动进行扫描，获得待检测管道外部环形区域的投影数据；

步骤2)根据投影数据重建管道外部环形区域的图像；

步骤3)显示重建图像。

2.根据权利要求1所述的大口径管道壁外部CT局部扫描成像方法，其特征在于：扫描开始时以射线源与原点连线所在的直线为y轴，并且以射线源指向原点的方向为正方向，x轴垂直于y轴，x轴与y轴构成固定的笛卡尔坐标系O-xy；在射线源和探测器围绕待检测管道作圆周运动的过程中，以坐标原点O建立旋转的笛卡尔坐标系O-ξη，η轴为扫描过程中射线源与原点连线所在的直线，并且以射线源指向原点的方向为正方向，ξ轴垂直于η轴，ξ轴与η轴构成右手笛卡尔坐标系；x轴与ξ轴的夹角为旋转角度θ。

3.根据权利要求1所述的大口径管道壁外部CT局部扫描成像方法，其特征在于：所述步骤2)根据投影数据重建管道外部环形区域的图像，具体包括以下步骤：

步骤2-1)凸集投影；

步骤2-2)总变差最小化；

步骤2-3)利用区域尺度拟合模型对子区域进行平均化修正。

4.根据权利要求3所述的大口径管道壁外部CT局部扫描成像方法，其特征在于：所述步骤2-1)凸集投影具体包括以下步骤：

步骤2-1-1)采用下面的加型代数迭代公式对待重建图像进行重建：

$\left\{\begin{array}{cc}{f}_j^{\left(0\right)}=0& \\ f_j^{\left(k\right)}=f_j^{(k-1)}+\mathrm{\&lambda;}\frac{p_i-\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}{w}_{ij}\&CenterDot;f_j^{(k-1)}}{\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}{w}_{ij}}{w}_{ij}& (j=1,2,...,N;i=1,2,...,M)\end{array}\right.$

设待重建图像的像素点总数为N个，f表示待重建的数字图像，位于(s,t)处的像素灰度值表示为f_s,t，则f可以表示成H×W的图像矩阵f＝(f_s,t)，其中H为待重建图像的行数，W为待重建图像的列数，将图像的像素点逐点排列成向量其中N＝H×W，设经过待重建图像的扫描射线数为M条，将射线投影数据按射线逐条排列为向量W＝(w_ij)为M×N维的投影系数矩阵，其中w_ij表示第j个点对第i条射线投影数据的贡献率；

其中，(k＝1,2,…,N_count)表示第k次迭代时图像数据向量的第j个分量，N_count表示重建算法的最大迭代次数，表示图像数据向量的初始值的第j个分量，p_i表示第i条射线对应的投影数据，k为迭代次数，λ为松弛因子；

步骤2-1-2)引入非负性限制，得到图像数据的校正值：

$f_j^{\left(POCS\right)}=\left\{\begin{array}{cc}{f}_j^{\left(M_1\right)},& \begin{array}{cc}if& f_j^{\left(M_1\right)}>0\end{array}\\ 0,& else\end{array}\right.,(j=1,2,...,N)$

其中，表示经过非负校正后得到的图像数据向量的第j个分量，表示按加型代数迭代公式经过M₁次迭代后的图像数据向量的第j个分量。

5.根据权利要求3所述的大口径管道壁外部CT局部扫描成像方法，其特征在于：所述步骤2-2)总变差最小化TVM，具体包括以下步骤：

步骤2-2-1)将总变差最小化的梯度下降方向f^(TVM-GRAD)初始化为f^(TVM-GRAD)＝f^(POCS)，将下降程度d_POCS初始化为d_POCS＝||f⁽⁰⁾-f^(POCS)||；通过以下公式计算图像的总变差TV(f)及在图像(s,t)处的逼近形式的偏导数v_s,t：

$TV\left(f\right)=\underset{s,t}{\&Sigma;}\sqrt{{(f_{s,t}-f_{s-1,t})}^2+{(f_{s,t}-f_{s,t-1})}^2+\mathrm{\&tau;}}$

$\begin{array}{c}{v}_{s,t}=\frac{\&part;TV\left(f\right)}{\&part;f_{s,t}}\&ap;\frac{(f_{s,t}-f_{s-1,t})+(f_{s,t}-f_{s,t-1})}{\sqrt{{(f_{s,t}-f_{s-1,t})}^2+{(f_{s,t}-f_{s,t-1})}^2+\mathrm{\&tau;}}}\\ -\frac{f_{s+1,t}-f_{s,t}}{\sqrt{{(f_{s+1,t}-f_{s,t})}^2+{(f_{s+1,t}-f_{s+1,t-1})}^2+\mathrm{\&tau;}}}\\ -\frac{f_{s,t+1}-f_{s,t}}{\sqrt{{(f_{s,t+1}-f_{s,t})}^2+{(f_{s,t+1}-f_{s-1,t+1})}^2+\mathrm{\&tau;}}}\end{array}$

步骤2-2-2)总变差梯度下降法迭代按下述公式进行：

$\left\{\begin{array}{c}{v}_{s,t}^{\left(i_1\right)}=\frac{\&part;TV\left(f^{(TV-GRAD)}\right)}{\&part;f_{s,t}^{(TV-GRAD)}}\\ f^{(TV-GRAD)}=f^{(TV-GRAD)}-{\mathrm{\&alpha;d}}_{POCS}\frac{v^{\left(i_1\right)}}{\left|\right|v^{\left(i_1\right)}\left|\right|}\end{array}\right.,i_1=1,2,...,N_{gard}$

其中，N_gard表示总变差梯度下降法的迭代次数，TV(f)表示图像数据f的总变差，τ为极小的正常数，表示第i₁次迭代图像像素点(s,t)处的总变差梯度下降方向，表示第i₁次迭代后整幅图像每个像素点处的总变差梯度下降方向矩阵，即1≤s≤H,1≤t≤W，H为待重建图像的行数，W为待重建图像的列数，||·||表示向量的Frobenius范数，f_s,t表示位于(s,t)处的像素点的灰度值，f_s-1,t表示位于(s-1,t)处的像素点的灰度值，f_s,t-1表示位于(s,t-1)处的像素点的灰度值，f_s+1,t表示位于(s+1,t)处的像素点的灰度值，f_s+1,t-1表示位于(s+1,t-1)处的像素点的灰度值，f_s,t+1表示位于(s,t+1)处的像素点的灰度值，f_s-1,t+1表示位于(s+1,t+1)处的像素点的灰度值，α为权系数；令f⁽⁰⁾＝f^(TVM-GRAD)，判断是否达到总变差最小化中预设的迭代次数N_TVM，如果是，则跳转至下一步骤S2-3)，否则跳转至步骤S2-1)。

6.根据权利要求3所述的大口径管道壁外部CT局部扫描成像方法，其特征在于：所述步骤利用区域尺度拟合模型对子区域进行平均化修正，具体包括以下步骤：

2-3-1)使用区域尺度拟合RSF活动轮廓模型提取待重建图像的边缘，通过求解以下梯度流演化方程得到水平集函数：

$\left\{\begin{array}{c}{f}_o\left(x\right)=\frac{K_{\mathrm{\&sigma;}}\left(x\right)*\&lsqb;H_{\mathrm{\&epsiv;}}\left(\mathrm{\&phi;}\right(x\left)\right)f\left(x\right)\&rsqb;}{K_{\mathrm{\&sigma;}}\left(x\right)*H_{\mathrm{\&epsiv;}}\left(\mathrm{\&phi;}\right(x\left)\right)};\\ f_b\left(x\right)=\frac{K_{\mathrm{\&sigma;}}\left(x\right)*\&lsqb;(1-H_{\mathrm{\&epsiv;}}(\mathrm{\&phi;}\left(x\right)))f\left(x\right)\&rsqb;}{K_{\mathrm{\&sigma;}}\left(x\right)*(1-H_{\mathrm{\&epsiv;}}(\mathrm{\&phi;}\left(x\right)))};\\ \frac{\&part;\mathrm{\&phi;}}{\&part;t}=-{\mathrm{\&delta;}}_{\mathrm{\&epsiv;}}\left(\mathrm{\&phi;}\right)({\mathrm{\&lambda;}}_1{e}_1\left(x\right)-{\mathrm{\&lambda;}}_2{e}_2\left(x\right))+{\mathrm{\&mu;\&delta;}}_{\mathrm{\&epsiv;}}\left(\mathrm{\&phi;}\right)div\left(\frac{\&dtri;\mathrm{\&phi;}}{|\&dtri;\mathrm{\&phi;}|}\right)+v({\&dtri;}^2\mathrm{\&phi;}-div(\frac{\&dtri;\mathrm{\&phi;}}{|\&dtri;\mathrm{\&phi;}|}));\\ e_1\left(x\right)=\&Integral;K_{\mathrm{\&sigma;}}(y-x)|f\left(x\right)-f_o\left(y\right){|}^{2}dy;e_2\left(x\right)=\&Integral;K_{\mathrm{\&sigma;}}(y-x)|f\left(x\right)-f_b\left(y\right){|}^{2}dy;\\ {\mathrm{\&delta;}}_{\mathrm{\&epsiv;}}\left(z\right)=\frac{1}{\mathrm{\&pi;}}\&CenterDot;\frac{\mathrm{\&epsiv;}}{{\mathrm{\&epsiv;}}^2+z^2};\\ H_{\mathrm{\&epsiv;}}\left(z\right)=\frac{1}{2}\&lsqb;1+\frac{2}{\mathrm{\&pi;}}\arctan \left(\frac{z}{\mathrm{\&epsiv;}}\right)\&rsqb;;\end{array}\right.$

其中，x,y为表示图像中像素点位置的二维坐标向量，f(x)表示图像在x处的灰度值，f_o(x)和f_b(x)分别为图像x处的局部区域内，轮廓线的内部和外部的像素加权强度均值，f_o(y)和f_b(y)分别为图像y处的局部区域内，轮廓线的内部和外部的像素加权强度均值，x∈R²为Gauss核函数，σ为尺度参数，*表示卷积运算，φ(x)为水平集函数，为水平集函数φ(x)的梯度，H_ε(z)为Heaviside函数的正则化形式，δ_ε(z)为维Dirac测度的正则化形式，ε为正常数，div(·)表示散度算子，表示Laplace算子，λ₁,λ₂>0,μ,v>0是各项的权值系数，t为引入的时间辅助变量；

步骤2-3-2)在得到RSF活动轮廓模型的水平集函数后，利用水平集函数将图像划分成不同的子区域，分别用每个子区域内的像素点灰度的平均值代替该子区域内各像素点的灰度值；判定是否达到设定的迭代次数N_e，如果是,则结束迭代，否则跳转至步骤S2-1)。