CN105447832A - 一种基于探测器单元标定的ct图像伪影校正方法及应用 - Google Patents

Abstract

本发明涉及一种基于探测器单元标定的CT图像伪影校正方法，通过几个已知均匀材料的模体，按每个探测器单元标定多能投影到材料等效厚度的映射关系；在实际CT成像中可以按照此映射关系，将待测物体的多能投影数据转换成已知材料的等效厚度，并对转换后的数据作滤波处理；对处理后的数据按通常CT算法进行图像重建。本发明方法操作简单，所需模体易于制作，降低了模体的加工难度，且不要求模体的精确尺寸和精确放置，降低了实际操作难度，也不需要知道射线的能谱信息，降低了实现的复杂度，能同时消除岩心CT图像中的缓变环状伪影和高频环状伪影，能够减轻由于衰减器引起的图像CT值畸变，该伪影校正方法还可以用于校正其它CT成像中的图像伪影。

Description

一种基于探测器单元标定的CT图像伪影校正方法及应用

技术领域

本发明属于检测仪器设备技术领域，涉及一种X射线CT图像伪影校正方法及应用，尤其是可应用于岩心CT，以及工业、医学等一类对象的CT图像伪影校正中。

背景技术

近年来，数字岩心分析技术快速发展，在致密油气、页岩油气开发中的作用逐渐显现。CT技术可以直接重构岩心的三维空间结构信息，成为构建三维真实数字岩心的重要方法之一。为构建数字岩心模型，需要对样品进行CT成像，并对其成分分布进行三维图像分割。而CT图像伪影会导致分割错误，进而影响数字建模的精确性。

标准岩心通常为圆柱状样品。为得到数字岩心模型，需要对其进行CT成像。众所周知，通常用于CT成像的射线由多能光子组成，CT数据获取的数学模型是非线性的，而常用CT图像重建算法如FBP、ART等均是基于单能光子的线性数学模型。从而导致岩心CT成像会产生图像伪影如图像环状伪影、硬化伪影以及CT值畸变等。这些图像伪影涉及多种因素，包括射线源因素(如能谱、流强稳定性)、滤波片因素(材质、几何形态)、衰减器因素(材质、几何形态)、探测器因素(如探测器单元对不同能量光子的探测效率差异，各个探测器单元的探测效率差异)和图像重建算法。

有关CT图像环状伪影、硬化伪影、CT值畸变校正方法已有很多研究，但现有方法尚不能较好地消除缓变环状伪影、硬化伪影和CT值畸变。本发明将给出一种基于探测器单元标定的CT图像伪影校正方法，通过检索，尚未检索到与本发明专利申请相关的专利公开文献。

发明内容

本发明的目的在于克服现有技术的不足之处，提供一种基于探测器单元标定的CT图像伪影校正方法，该方法能较好地消除岩心CT图像中的缓变环状伪影和高频环状伪影，同时能减轻由于衰减器射束硬化而引起的图像CT值畸变，该伪影校正方法还可以用于校正其它CT成像中的图像伪影，如工业构件CT图像伪影、乳腺CT图像伪影等。

为了实现上述目的，本发明所采用的技术方案如下：

一种基于探测器单元标定的CT图像伪影校正方法，步骤如下：

⑴通过几个已知均匀材料的模体，按每个探测器单元获取一系列材料等效厚度下的多能投影数据的采样值；

⑵由该系列已知材料等效厚度下的多能投影数据的采样值，以及之间的固有性质拟合出每个探测器单元的多能投影与材料等效厚度的函数关系的反函数，即多能投影到材料等效厚度的映射关系；

⑶在实际CT成像中，获得待校正的样品的多能投影数据，按照上述拟合得到的反函数，将多能投影数据转换成已知材料的等效厚度，并对转换后的数据作Wavelet-FFT滤波处理；

⑷对处理后的数据按通常CT算法进行图像重建，即可得到伪影校正后的CT图像。

而且，对每个探测器单元标定多能投影到材料等效厚度的映射关系，不同探测器单元对应的映射关系不完全相同。

而且，所述步骤⑴中获取多能投影数据的采样值时，同时使用几个不同尺寸的模体进行扫描，或者，将模体放置在视野中的不同位置多次扫描，根据扫描数据的CT图像，利用分割的方法获取一系列材料等效厚度下的多能投影数据的采样值。

而且，所述步骤⑴中获取多能投影数据的采样值为利用已知厚度的均匀材料的板状模体的多能投影数据获取一系列材料等效厚度与多能投影数据的对应关系的采样值。

而且，所述步骤⑵的具体方法为：根据设备要求采集充分多的采样值，直接拟合每个探测器单元的多能投影数据到材料等效厚度的映射关系；或者，采集一部分等效厚度下的多能投影的采样值，利用材料等效厚度与多能投影的固有性质拟合出每个探测器单元的多能投影到材料等效厚度的映射关系。

而且，所述步骤⑴中几个已知均匀材料的模体是指已知模体的厚度尺寸、模体个数和模体的X射线吸收系数，其模体个数和厚度尺寸由成像设备和放置位置决定，但在给定成像设备和放置位置条件下模体个数和厚度尺寸并不唯一；或者，所述步骤⑴中模体的形状为柱状、锥状、球状、椭球状、圆台状、棱台状、球冠状或椭球冠状物体。

而且，具体步骤如下：

⑴数学模型：

岩心CT扫描***，由射线源、探测器、机械旋转***、衰减器以及控制和计算机构成，岩心CT的数学模型如下：

$I\left(u,\mathrm{\β}\right)={\∫}_{E_{\min }}^{E_{\max }}{I}_0\left(E,u\right)\mathrm{\γ}\left(E,u\right)\exp \left(-{\mathrm{\μ}}_f\left(E\right)r\left(u\right)\right)\exp \left(-\underset{l\∈L\left(u\right)}{\∫}\mathrm{\μ}\left(xR\left(\mathrm{\β}\right),E\right)dl\right)dE+\mathrm{\σ}\left(u\right),---\left(A\right)$

其中，x表示固定坐标系中的点，u为探测器坐标，L(u)表示从射线源到u的射线，β为岩心绕定轴旋转的角度，R(β)为旋转矩阵，μ(x,E)表示初始时刻岩心对能量为E的光子的线性衰减系数分布，μ_f(E)表示衰减器单位长度对能量E的光子的线性衰减系数，r(u)为射线到达探测器单元u所透过的衰减器的厚度，γ(E,u)表示探测器单元u的射线探测效率，I₀(E,u)代表能量为E的光子的入射强度，其中E_min和E_max分别表示光子能量的最小值和最大值，I(u,β)代表探测器单元u在角度β所采集的光子数，σ(u)表示其中包含的散射光子数；

当不放扫描物体时，探测器探测到的光子数可以表示为：

${\hat{I}}_0\left(u\right)={\∫}_{E_{\min }}^{E_{\max }}{I}_0(E,u)\mathrm{\γ}(E,u)\exp (-{\mathrm{\μ}}_f\left(E\right)r\left(u\right))dE+{\mathrm{\σ}}_0\left(u\right),---\left(B\right)$

于是，被测物体的多能投影数据表示为

$\begin{array}{c}p(u,\mathrm{\β})=-log\frac{I(u,\mathrm{\β})}{{\hat{I}}_0(E,u)}\\ =-\log \frac{{\∫}_{E_{\min }}^{E_{\max }}{I}_0\left(E,u\right)\mathrm{\γ}\left(E,u\right)\exp \left(-{\mathrm{\μ}}_f\left(E\right)r\left(u\right)\right)\exp \left(-\underset{l\∈L\left(u\right)}{\Σ}\mathrm{\μ}\left(xR\left(\mathrm{\β}\right),E\right)dl\right)dE+\mathrm{\σ}\left(u\right)}{{\∫}_{E_{\min }}^{E_{\max }}{I}_0\left(E,u\right)\mathrm{\γ}\left(E,u\right)\exp \left(-{\mathrm{\μ}}_f\left(E\right)r\left(u\right)\right)dE+{\mathrm{\σ}}_0\left(u\right)};\end{array}---\left(C\right)$

⑵多能投影数据与均匀材料厚度的函数关系：

当X射线由多能光子组成时，多能投影数据由公式(3)给出，当被测物体为单材质物体时，即μ(x,E)＝μ₀(E)ρ(x)，由公式(3)，可以得到

$p=p(t,u)\stackrel{def}{=}-log\frac{{\∫}_{E_{\min }}^{E_{\max }}\mathrm{\γ}(E,u)I_0(E,u)\exp (-{\mathrm{\μ}}_f(E)r\left(u\right))\exp (-{\mathrm{\μ}}_0(E\left)t\right)dE+\mathrm{\σ}\left(u\right)}{{\∫}_{E_{\min }}^{E_{\max }}\mathrm{\γ}(E,u)I_0(E,u)\exp (-{\mathrm{\μ}}_f(E)r\left(u\right))dE+{\mathrm{\σ}}_0\left(u\right)},---\left(D\right)$

其中

$t=\underset{l\∈L\left(u\right)}{\∫}\mathrm{\ρ}\left(xR\right(\mathrm{\β}\left)\right)dl,$

p＝p(t,u)反映了探测器单元u采集到的多能投影数据与该材料等效厚度的函数关系；

⑶由均匀材料圆柱状模体恢复t＝t(p,u)的方法：

①利用相同材料制作若干个均匀材质的模体；

②用CT扫描这些模体并由多能投影数据重建图像；

③通过对图像分割，按探测器单元确定每个多能投影对应的模体等效厚度获取一系列数据对，或者根据模体的已知厚度和其多能投影值获取每个探测器单元对应的一系列数据对；

④建立恢复t＝t(p,u)的优化数学模型，并从一系列模体等效厚度和多能投影数据对恢复t＝t(p,u)；

具体为：

记模体的密度函数为μ(x,E)＝μ₀(E)ρ(x)，通过CT扫描得到没有物体时的空扫描数据I₀(u_j)和加载物体后的扫描数据I(u_j,β_k),j∈J,k∈K，此处j和k分别是探测器单元序号和角度采样序号，于是可以得到一组多能投影数据j∈J,k∈K，由此数据直接重建一副CT图像其中可能含有噪声和CT值畸变，将图像分割，对每个p_k,j计算出t_k,j，于是得到U＝{(t_k,j,p_k,j),j∈J,k∈K}；对逐个探测器单元u_j，用多项式逼近t＝t(p,u_j)，即假设探测器u＝u_j对应的多项式为

$t=t_j(p;a_{0j},a_{1j},L,a_{nj})\stackrel{def}{=}{a}_{0j}+a_{1j}p+L+a_{nj}{p}^n$

由U＝{(t_k,j,p_k,j),j∈J,k∈K}恢复t＝t(p,u_j)的优化问题如下：

$\underset{\left\{a_0\left(u_j\right),a_1\left(u_j\right),L{a}_n\left(u_j\right)\right\}}{\arg \min }{\Σ}_{k=1}^K{\left(t_j\right(p_{k,j};a_{0j},a_{1j},L,a_{nj})-t_{k,j})}^2,$

s.t.t′_j(p；a_0j,a_1j,L,a_nj)≥0,t″_j(p；a_0j,a_1j,L,a_nj)≥0；(E)

⑷优化问题求解方法：

令α＝{a_0j,a_1j,L,a_nj}，则目标函数表示为多项式的一阶导数为 $g_k\left(\mathrm{\α}\right)=a_{1j}+2a_{2j}{p}_{k,j}+L+{\mathrm{na}}_{nj}{p}_{k,j}^{n-1},$ 多项式的二阶导数为 $h_k\left(\mathrm{\α}\right)=2a_{2j}+6a_{3j}{p}_{k,j}+L+n(n-1)a_{nj}{p}_{k,j}^{n-2},$ 则约束条件表示为：g_k(α)≥0，h_k(α)≥0，其中k＝1,2,L,K；问题归结为求解带不等式约束的优化问题：

$\underset{\mathrm{\α}}{\min }{E}_j\left(\mathrm{\α}\right)$

s.t.g_k,j(α)≥0，

h_k,j(α)≥0，k＝1,…K

求解上述问题得到每个探测器单元对应的模体等效厚度关于多能投影数据的函数映射关系。

如上所述的基于探测器单元标定的CT图像伪影校正方法在CT图像伪影校正中的应用。

而且，所述校正方法在多种不同的扫描方式中的应用；或者，所述校正方法在一类单材质占优物体的伪影校正中的应用；或者，所述校正方法在岩心三维CT成像图像伪影校正、柱状物体三维CT图像伪影校正、乳腺三维CT图像伪影校正或口腔CT图像伪影校正方面的应用。

本发明的优点和积极效果是：

本发明方法操作简单，所需模体易于制作，降低了模体的加工难度，且不要求模体的精确尺寸和精确放置，降低了实际操作难度，也不需要知道射线的能谱信息，降低了实现的复杂度，该方法能同时消除岩心CT图像中的缓变环状伪影和高频环状伪影，该方法能够减轻由于衰减器引起的图像CT值畸变，该伪影校正方法还可以用于校正其它CT成像中的图像伪影，如工业构件CT图像伪影、乳腺CT图像伪影等。

附图说明

图1为本发明中岩心的CT扫描示意图；其中，图1-1为衰减器靠近探测器一端的岩心的CT扫描示意图，图1-2为衰减器靠近射线源一端的岩心的CT扫描示意图，1为衰减器，2为射线源，3为被测件，4为探测器；

图2为本发明中衰减器和模体实体图；其中，图2-1为衰减器图，图2-2为模体图；

图3为本发明中模体CT图像及分割结果图；其中，图3-1为模体CT图像，图3-2为模体CT图像分割结果图；

图4为本发明中多能投影到材料等效厚度的映射关系图；其中，图4-1相应于u＝681，图4-2相应于u＝868；

图5为本发明中岩心CT图像伪影校正对比实验结果；其中，图5-1为由多能投影数据直接重建的图像，灰度窗[0,0.145]；图5-2为用Wavelet-FFT对多能投影数据滤波后重建的图像，灰度窗[0,0.145]；图5-3为用本发明方法对多能投影数据校正后重建的图像，灰度窗[0,0.062]；

图6为图5的局部放大图像；其中，图6-1为图5-1的局部放大图像，灰度窗[0.06,0.145]；图6-2为图5-2的局部放大图像，灰度窗[0.06,0.145]；图6-3为图5-3的局部放大图像，灰度窗[0.03,0.06]；

图7为本发明中铝校正模体和测试模体的CT图像；其中，图7-1为校正模体的CT图像，图7-2为测试模体的CT图像；

图8为本发明中铝模体CT图像伪影校正对比结果；其中，图8-1为由多能投影数据直接重建的图像，灰度窗[0,0.176]；图8-2用Wavelet-FFT对多能投影数据滤波后重建的图像，灰度窗[0,0.176]；图8-3为本发明方法对多能投影数据校正后重建的图像，灰度窗[0,0.087]；

图9为图8的局部放大图像；其中，图9-1为图8-1的局部放大图像，灰度窗[0,0.176]；图9-2图8-2的局部放大图像，灰度窗[0,0.176]，图9-3为图8-3的局部放大图像，灰度窗[0,0.087]。

具体实施方式

下面结合实施例，对本发明进一步说明；下述实施例是说明性的，不是限定性的，不能以下述实施例来限定本发明的保护范围。

本发明中所使用的设备，如无特殊规定，均为本领域内常用的设备；本发明中所使用的方法，如无特殊规定，均为本领域内常用的方法。

本发明中，所述均匀材料的模体，通常是指不同直径圆柱状模体或者不同厚度板状模体，但不限于这些形状，具体个数视具体情况而定；所述材料等效厚度，可以用已知厚度的板状模体获得，也可以通过对圆柱状模体由多能投影数据重建的CT图像分割获得；所述模体材料选择与扫描物体的X射线吸收系数相同或者相近的材料；所述多能投影与材料等效厚度之间的固有性质是指其对应函数的一阶导数大于等于0；所述多能投影与材料等效厚度之间的固有性质是指其对应函数的二阶导数大于等于0；所述材料的等效厚度，是指该材料的自然厚度或者单能光子穿过自然厚度的衰减值。

实施例1

一种基于探测器单元标定的CT图像伪影校正方法，包括步骤：通过几个已知均匀材料的模体，按每个探测器单元获取一系列材料等效厚度下的多能投影的采样值；由该系列已知材料等效厚度下的多能投影采样值，以及之间的固有性质拟合出每个探测器单元的“多能投影与材料等效厚度的函数关系的反函数”，即“多能投影到材料等效厚度的映射”；在实际CT成像中可以按照此映射关系，将多能投影数据转换成已知材料的等效厚度，并对转换后的数据作Wavelet-FFT滤波处理；对处理后的数据按通常CT算法进行图像重建；

较优地，对每个探测器单元标定多能投影到材料等效厚度的映射关系，不同探测器单元对应的映射关系不完全相同。

较优地，可以同时用几个不同尺寸的模体进行扫描，也可以将模体放置在视野中的不同位置多次扫描。根据扫描数据的CT图像，利用分割的方法获取一系列材料等效厚度下的多能投影的采样值。

较优地，利用已知厚度的均匀材料的板状模体的多能投影数据获取一系列材料等效厚度与多能投影的对应关系的采样值。

较优地，可以根据设备要求采集充分多的采样值，直接用最小二乘方法拟合每个探测器单元的“多能投影到材料等效厚度的映射关系”；也可以采集一部分等效厚度下的多能投影的采样值，利用材料等效厚度与多能头一个的固有性质拟合出每个探测器单元的“多能投影到材料等效厚度的映射关系”。

较优地，分别采集几个已知厚度的铝板的多能投影数据获取一系列材料等效厚度下的多能投影的采样值。利用这些采样值拟合每个探测器单元的“多能投影到材料等效厚度的映射关系”。

较优地，铝板的尺寸由成像设备和放置位置决定。

较优地，模体形状可以是圆柱状、板状、椭圆状，但不限于这些形状。

本方法适用于多种不同的扫描方式；本方法适用于一类单材质占优物体的环状伪影校正，较优地，物体组成成分中某一类成分占总成分的50％以上。

本方法可以用于岩心三维CT成像图像伪影校正、柱状物体三维CT图像伪影校正、乳腺三维CT图像伪影校正、口腔CT图像伪影校正等方面的应用。

实施例2

一种基于探测器单元标定的岩心CT图像伪影校正方法，步骤如下：

一、数学模型

岩心CT扫描***，由如图1所示的射线源、探测器、机械旋转***、衰减器以及控制和计算机构成。岩心CT的数学模型如下：

$I\left(u,\mathrm{\β}\right)={\∫}_{E_{\min }}^{E_{\max }}{I}_0\left(E,u\right)\mathrm{\γ}\left(E,u\right)\exp \left(-{\mathrm{\μ}}_f\left(E\right)r\left(u\right)\right)\exp \left(-\underset{l\∈L\left(u\right)}{\∫}\mathrm{\μ}\left(xR\left(\mathrm{\β}\right),E\right)dl\right)dE+\mathrm{\σ}\left(u\right)---\left(A\right)$

其中，x表示固定坐标系中的点，u为探测器坐标，L(u)表示从射线源到u的射线，β为岩心绕定轴旋转的角度，R(β)为旋转矩阵，μ(x,E)表示初始时刻岩心对能量为E的光子的线性衰减系数分布，μ_f(E)表示衰减器单位长度对能量E的光子的线性衰减系数，r(u)为射线到达探测器单元u所透过的衰减器的厚度，γ(E,u)表示探测器单元u的射线探测效率，I₀(E,u)代表能量为E的光子的入射强度，其中E_min和E_max分别表示光子能量的最小值和最大值，I(u,β)代表探测器单元u在角度β所采集的光子数，σ(u)表示其中包含的散射光子数。如上所述，衰减器的作用是在保证岩心圆柱体中心部分相应的探测器单元有足够的计数，又保证射线经过岩心圆柱体边缘部分和未经过岩心所到达的探测器单元计数不过载。因此，衰减器设计为中间较薄边缘较厚的弧形。注意在模型(1)中，假定了探测器单元u的对入射光子的探测效率与入射光子强度呈线性关系。

当不放扫描物体时，探测器探测到的光子数可以表示为：

${\hat{I}}_0\left(u\right)={\∫}_{E_{\min }}^{E_{\max }}{I}_0(E,u)\mathrm{\γ}(E,u)\exp (-{\mathrm{\μ}}_f(E)r\left(u\right))dE+{\mathrm{\σ}}_0\left(u\right)---\left(B\right)$

于是，被测物体的多能投影数据表示为

$\begin{array}{c}p(u,\mathrm{\β})=-log\frac{I(u,\mathrm{\β})}{{\hat{I}}_0(E,u)}\\ =-\log \frac{{\∫}_{E_{\min }}^{E_{\max }}{I}_0\left(E,u\right)\mathrm{\γ}\left(E,u\right)\exp \left(-{\mathrm{\μ}}_f\left(E\right)r\left(u\right)\right)\exp \left(-\underset{l\∈L\left(u\right)}{\Σ}\mathrm{\μ}\left(xR\left(\mathrm{\β}\right),E\right)dl\right)dE+\mathrm{\σ}\left(u\right)}{{\∫}_{E_{\min }}^{E_{\max }}{I}_0\left(E,u\right)\mathrm{\γ}\left(E,u\right)\exp \left(-{\mathrm{\μ}}_f\left(E\right)r\left(u\right)\right)dE+{\mathrm{\σ}}_0\left(u\right)},\end{array}---\left(C\right)$

模体已知，表示模体的厚度就知道了，同时公式中的μ(x,E)和也知道了，从而对应的多色投影值I(u,β)就知道了。u取不同值对应不同的探测器单元，β取不同值对应不同的厚度值，从而一系列的β就得到一系列的厚度。

二、多能投影数据的校正方法

当X射线由多能光子组成时，多能投影数据由公式(C)给出。特别被测物体为单材质物体时，即μ(x,E)＝μ₀(E)ρ(x)，由公式(C)，本申请人可以得到

$p=p(t,u)\stackrel{def}{=}-log\frac{{\∫}_{E_{\min }}^{E_{\max }}\mathrm{\γ}(E,u)I_0(E,u)\exp (-{\mathrm{\μ}}_f(E)r\left(u\right))\exp (-{\mathrm{\μ}}_0(E\left)t\right)dE+\mathrm{\σ}\left(u\right)}{{\∫}_{E_{\min }}^{E_{\max }}\mathrm{\γ}(E,u)I_0(E,u)\exp (-{\mathrm{\μ}}_f(E)r\left(u\right))dE+{\mathrm{\σ}}_0\left(u\right)},---\left(D\right)$

其中

$t=\underset{l\∈L\left(u\right)}{\∫}\mathrm{\ρ}\left(xR\right(\mathrm{\β}\left)\right)dl.$

p＝p(t,u)反映了探测器单元u采集到的多能投影数据与该材料等效厚度的函数关系。

下面介绍本申请人的由均匀材料圆柱状模体恢复t＝t(p,u)的方法。基本思路是：(1)利用相同材料制作若干个均匀材质的模体；(2)用CT扫描这些模体并由多能投影数据重建图像；(3)通过对图像分割，按探测器单元确定每个多能投影对应的模体等效厚度；(4)建立恢复t＝t(p,u)的优化数学模型，并从一系列“模体等效厚度和多能投影”数据对恢复t＝t(p,u)。

记模体的密度函数为μ(x,E)＝μ₀(E)ρ(x)，通过CT扫描，本申请人可以得到没有物体时的空扫描数据I₀(u_j)和加载物体后的扫描数据I(u_j,β_k),j∈J,k∈K，此处j和k分别是探测器单元序号和角度采样序号。于是可以得到一组多能投影数据j∈J,k∈K.由此数据可以直接重建一副CT图像其中可能含有噪声和CT值畸变。将图像分割，对每个p_k,j计算出t_k,j，于是得到U＝{(t_k,j,p_k,j),j∈J,k∈K}。本发明中对逐个探测器单元u_j，用多项式逼近t＝t(p,u_j)，即假设探测器u＝u_j对应的多项式为

$t=t_j(p;a_{0j},a_{1j},L,a_{nj})\stackrel{def}{=}{a}_{0j}+a_{1j}p+L+a_{nj}{p}^n$

由U＝{(t_k,j,p_k,j),j∈J,k∈K}恢复t＝t(p,u_j)的优化问题如下：

$\underset{\left\{a_0\left(u_j\right),a_1\left(u_j\right),L{a}_n\left(u_j\right)\right\}}{\arg \min }{\Σ}_{k=1}^K{\left(t_j\right(p_{k,j};a_{0j},a_{1j},L,a_{nj})-t_{k,j})}^2,$

s.t.t′_j(p；a_0j,a_1j,L,a_nj)≥0,t″_j(p；a_0j,a_1j,L,a_nj)≥0.(E)

接下来，给出优化问题(E)的求解方法。令α＝{a_0j,a_1j,L,a_nj}，则目标函数表示为 $E_j={\Σ}_{k=1}^K{\left(t_j\right(p_{k,j};\mathrm{\α})-t_{k,j})}^2.$ 多项式的一阶导数为 $g_k\left(\mathrm{\α}\right)=a_{1j}+2a_{2j}{p}_{k,j}+L+{\mathrm{na}}_{nj}{p}_{k,j}^{n-1},$ 多项式的二阶导数为则约束条件表示为：g_k(α)≥0，h_k(α)≥0，其中问题归结为求解带不等式约束的优化问题：

$\underset{\mathrm{\α}}{\min }{E}_j\left(\mathrm{\α}\right)$

s.t.g_k,j(α)≥0

h_k,j(α)≥0，k＝1,…K

该问题可化为无约束优化问题进行求解。

求解该问题得到每个探测器单元对应的模体等效厚度关于多能投影数据的函数映射关系。

例如，利用罚函数方法将带约束的优化问题转化为无约束的优化问题，构造如下增广目标函数：

H_j(α,τ)＝E_j(α)+λP_j(α)^τ

其中λ＞0为惩罚因子， $P_j\left(\mathrm{\α}\right)=\frac{1}{g_{k,j}\left(\mathrm{\α}\right)}+L+\frac{1}{g_{K,j}\left(\mathrm{\α}\right)}+\frac{1}{h_{1,j}\left(\mathrm{\α}\right)}+L+\frac{1}{h_{K,j}\left(\mathrm{\α}\right)}$ 或者

P_j(α)＝ln(g_1,j(α))+L+ln(g_K,j(α))+ln(h_1,j(α))+L+ln(h_K,j(α))，该无约束优化问题可以使用最速下降方法进行求解。

本发明基于探测器单元标定的CT图像伪影校正方法的相关检测结果：

首先，本申请人测试了一个直径100mm的岩心样品。使用的NTB探测器的最大计数是4096。流强较大时，该探测器很容易饱和。为了避免这种情况的出现，在扫描岩芯时，将衰减器置于探测器前面。衰减器如图2-1所示，需要注意的是衰减器的材质与样品的材质没有相关关系。数据采集时，射线源的管电压为140kV，管电流为12mA。线探测器含有3710个探测器单元，每个探测器单元的尺寸为0.083mm。射线源焦点到转台中心的距离为730mm，射线源焦点到探测器的距离为916mm。旋转一周采集1800个角度的投影数据。

岩心样品往往由不用种类的化合物按不同的比例组成，并且其密度分布也不够均匀。不可能找到一种与其X射线吸收系数完全相同的材料进行标定。实验中用的岩心样品主要成分为碳酸盐，本申请人通过不同材质的模体对该岩心样品进行标定，实验表明铝的线性衰减系数与该岩心近似。以下岩心CT实验中，都用铝材质的模体来标定样品等效厚度与多能投影值之间的函数关系。

模体的尺寸和扫描时的放置位置会影响到每个探测器单元获取的数据范围。实际上，本申请人可以将不同尺寸的模体分别扫描，或者放置不同的位置分别扫描。同时使用多组数据，保证每个探测器单元都能够采集到足够多的不同厚度的数据。本发明利用带约束的优化模型求解该问题，本申请人只是选取了两个不同尺寸的铝圆柱同时扫描，如图2-2所示，其中较小铝块的直径为40mm，较大铝块的直径为80mm。图3-1是模体的CT图像，图3-2给出了其分割结果。

根据模体CT图像的分割结果获取等效厚度，并与模体的多能投影数据组成数据对。由于模体的放置位置和尺寸影响，不同探测器单元获取的数据对分布不同。本申请人给出了两个探测器单元对应的数据对拟合结果，其中u＝681靠近视野的边缘，u＝868靠近视野的中心，如图4所示。

直接使用带衰减器的多能投影数据重建结果如图5-1所示，其中带有明显的环状伪影。由于衰减器对岩心的厚度进行了有效补偿，从而使得岩心的灰度值畸变并不明显。图5-2给出了用Wavelet-FFT对多能投影数据滤波后重建的图像，其中滤波参数中小波分解层数为5，方差为2，小波基函数为‘db25’。从图像上看，高频环状伪影明显减轻，但是存在一些低频的环状伪影。图5-3给出了用本发明方法对多能投影数据校正后重建图像，其中先对多能投影数据进行标定，之后再用Wavelet-FFT进行滤波，图5-3中低频和高频环状伪影被有效去除。对图5的局部细节进行放大，结果如图6所示。

接下来，本申请人选择了两个直径相同的铝圆柱，直径为100mm，结构如图7所示。其中一个不带孔洞，用来标定不同探测器单元对应的函数关系，另外一个带有四个不同尺寸孔洞的圆柱作为测试模体。本申请人在测试模体的中心孔洞里面填充了一个装有水的塑料管，在边缘的孔洞里面填充了一段实心的塑料管。

本实验设备与上述岩心CT设备不同，使用的是YXLONY.LDA探测器，最大计数为65536，在扫描过程中没有出现探测器计数饱和的问题，因此没有使用衰减器。射线源(YXLONY.TU450D09tube)的管电压为185kV，管电流为3.5mA。该线探测器阵列(YXLONY.LDA)含有2604个探测器，每个探测器的尺寸为0.25mm。探测器的采样时间为50ms。射线源焦点到转台中心的距离为730mm，射线源焦点到探测器的距离为916mm。旋转一周采集1800个角度的投影数据。

图8-1给出了由测试模体的多能投影数据直接重建的CT图像，其中具有明显的环状伪影和杯状伪影。图8-2给出了用Wavelet-FFT对多能投影数据滤波后重建的CT图像，其中存在一些低频的环状伪影，并且还有明显的杯状伪影。图8-3给出了用本发明方法对多能投影数据校正后重建的CT图像，环状伪影和杯状伪影得到了明显的去除。

图9为图8的局部放大图像，从图像上可以清晰的分辨出水、塑料和铝之间的差异。

Claims (9)

1.一种基于探测器单元标定的CT图像伪影校正方法，其特征在于：步骤如下：

⑴通过几个已知均匀材料的模体，按每个探测器单元获取一系列材料等效厚度下的多能投影数据的采样值；

⑵由该系列已知材料等效厚度下的多能投影数据的采样值，以及之间的固有性质拟合出每个探测器单元的多能投影与材料等效厚度的函数关系的反函数，即多能投影到材料等效厚度的映射关系；

⑶在实际CT成像中，获得待校正的样品的多能投影数据，按照上述拟合得到的反函数，将多能投影数据转换成已知材料的等效厚度，并对转换后的数据作Wavelet-FFT滤波处理；

⑷对处理后的数据按通常CT算法进行图像重建，即可得到伪影校正后的CT图像。

2.根据权利要求1所述的基于探测器单元标定的CT图像伪影校正方法，其特征在于：对每个探测器单元标定多能投影到材料等效厚度的映射关系，不同探测器单元对应的映射关系不完全相同。

3.根据权利要求1所述的基于探测器单元标定的CT图像伪影校正方法，其特征在于：所述步骤⑴中获取多能投影数据的采样值时，同时使用几个不同尺寸的模体进行扫描，或者，将模体放置在视野中的不同位置多次扫描，根据扫描数据的CT图像，利用分割的方法获取一系列材料等效厚度下的多能投影数据的采样值。

4.根据权利要求1所述的基于探测器单元标定的CT图像伪影校正方法，其特征在于：所述步骤⑴中获取多能投影数据的采样值为利用已知厚度的均匀材料的板状模体的多能投影数据获取一系列材料等效厚度与多能投影数据的对应关系的采样值。

5.根据权利要求1所述的基于探测器单元标定的CT图像伪影校正方法，其特征在于：所述步骤⑵的具体方法为：根据设备要求采集充分多的采样值，直接拟合每个探测器单元的多能投影数据到材料等效厚度的映射关系；或者，采集一部分等效厚度下的多能投影的采样值，利用材料等效厚度与多能投影的固有性质拟合出每个探测器单元的多能投影到材料等效厚度的映射关系。

6.根据权利要求1所述的基于探测器单元标定的CT图像伪影校正方法，其特征在于：所述步骤⑴中几个已知均匀材料的模体是指已知模体的厚度尺寸、模体个数和模体的X射线吸收系数，其模体个数和厚度尺寸由成像设备和放置位置决定，但在给定成像设备和放置位置条件下模体个数和厚度尺寸并不唯一；或者，所述步骤⑴中模体的形状为柱状、锥状、球状、椭球状、圆台状、棱台状、球冠状或椭球冠状物体。

7.根据权利要求1所述的基于探测器单元标定的CT图像伪影校正方法，其特征在于：具体步骤如下：

⑴数学模型：

岩心CT扫描***，由射线源、探测器、机械旋转***、衰减器以及控制和计算机构成，岩心CT成像的数学模型如下：

$I(u,\mathrm{\β})={\∫}_{E_{\min }}^{E_{max}}{I}_0(E,u)\mathrm{\γ}(E,u)\exp (-{\mathrm{\μ}}_f(E\left)r\right(u\left)\right)\exp (-\underset{l\∈L\left(u\right)}{\∫}\mathrm{\μ}(xR\left(\mathrm{\β}\right),E\left)dl\right)dE+\mathrm{\σ}\left(u\right)---\left(A\right)$

其中，x表示固定坐标系中的点，u为探测器坐标，L(u)表示从射线源到u的射线，β为岩心绕定轴旋转的角度，R(β)为旋转矩阵，μ(x,E)表示初始时刻岩心对能量为E的光子的线性衰减系数分布，μ_f(E)表示衰减器单位长度对能量E的光子的线性衰减系数，r(u)为射线到达探测器单元u所透过的衰减器的厚度，γ(E,u)表示探测器单元u的射线探测效率，I₀(E,u)代表能量为E的光子的入射强度，其中E_min和E_max分别表示光子能量的最小值和最大值，I(u,β)代表探测器单元u在角度β所采集的光子数，σ(u)表示其中包含的散射光子数；

当不放扫描物体时，探测器探测到的光子数可以表示为：

${\hat{I}}_0\left(u\right)={\∫}_{E_{\min }}^{E_{max}}{I}_0(E,u)\mathrm{\γ}(E,u)\exp (-{\mathrm{\μ}}_f(E\left)r\right(u\left)\right)dE+{\mathrm{\σ}}_0\left(u\right),---\left(B\right)$

于是，被测物体的多能投影数据表示为

$\begin{array}{c}p(u,\mathrm{\β})=-log\frac{I(u,\mathrm{\β})}{{\hat{I}}_0(E,u)}\\ =-\log \frac{{\∫}_{E_{\min }}^{E_{max}}{I}_0(E,u)\mathrm{\γ}(E,u)\exp (-{\mathrm{\μ}}_f(E\left)r\right(u\left)\right)\exp (-\underset{l\∈L\left(u\right)}{\∫}\mathrm{\μ}(xR\left(\mathrm{\β}\right),E\left)dl\right)dE+\mathrm{\σ}\left(u\right)}{{\∫}_{E_{\min }}^{E_{max}}{I}_0(E,u)\mathrm{\γ}(E,u)\exp (-{\mathrm{\μ}}_f(E\left)r\right(u\left)\right)dE+{\mathrm{\σ}}_0\left(u\right)}\end{array};---\left(C\right)$

⑵多能投影数据与均匀材料厚度的函数关系

当X射线由多能光子组成时，多能投影数据由公式(C)给出，当被测物体为单材质物体时，即μ(x,E)＝μ₀(E)ρ(x)，由公式(C)，可以得到

$p=p(t,u)\stackrel{def}{=}-log\frac{{\∫}_{E_{\min }}^{E_{\max }}\mathrm{\γ}(E,u)I_0(E,u)\exp (-{\mathrm{\μ}}_f(E\left)r\right(u\left)\right)\exp (-{\mathrm{\μ}}_0(E\left)t\right)dE+\mathrm{\σ}\left(u\right)}{{\∫}_{E_{\min }}^{E_{max}}\mathrm{\γ}(E,u)I_0(E,u)\exp (-{\mathrm{\μ}}_f(E\left)r\right(u\left)\right)dE+{\mathrm{\σ}}_0\left(u\right)},---\left(D\right)$

其中

$t=\underset{l\∈L\left(u\right)}{\∫}\mathrm{\ρ}\left(xR\right(\mathrm{\β}\left)\right)dl;$

p＝p(t,u)反映了探测器单元u采集到的多能投影数据与该材料等效厚度的函数关系；

⑶由均匀材料圆柱状模体恢复t＝t(p,u)的方法：

①利用相同材料制作若干个均匀材质的模体；

②用CT扫描这些模体并由多能投影数据重建图像；

③通过对图像分割，按探测器单元确定每个多能投影对应的模体等效厚度获取一系列数据对，或者根据模体的已知厚度和其多能投影值获取每个探测器单元对应的一系列数据对；

④建立恢复t＝t(p,u)的优化数学模型，并从一系列模体等效厚度和多能投影数据对恢复t＝t(p,u)；

具体为：

记模体的密度函数为μ(x,E)＝μ₀(E)ρ(x)，通过CT扫描得到没有物体时的空扫描数据I₀(u_j)和加载物体后的扫描数据I(u_j,β_k),j∈J,k∈K，此处j和k分别是探测器单元序号和角度采样序号，于是可以得到一组多能投影数据j∈J,k∈K，由此数据直接重建一副CT图像其中可能含有噪声和CT值畸变，将图像分割，对每个p_k,j计算出t_k,j，于是得到U＝{(t_k,j,p_k,j),j∈J,k∈K}；对逐个探测器单元u_j，用多项式逼近t＝t(p,u_j)，即假设探测器u＝u_j对应的多项式为

$t=t_j(p;a_{0j},a_{1j},L,a_{nj})\stackrel{def}{=}{a}_{0j}+a_{1j}p+L+a_{nj}{p}^n$

由U＝{(t_k,j,p_k,j),j∈J,k∈K}恢复t＝t(p,u_j)的优化问题如下：

$\underset{\{a_0\left(u_j\right),a_1\left(u_j\right),L{a}_n\left(u_j\right)\}}{\arg \min }{\mathrm{\Σ}}_{k=1}^K{\left(t_j\right(p_{k,j};a_{0j},a_{1j},L,a_{nj})-t_{k,j})}^2,$

$\begin{array}{cc}s.t.& t_j^{\′}(p;a_{0j},a_{1j},L,a_{nj})\≥0,t_j^{\′\′}(p;a_{0j},a_{1j},L,a_{nj})\≥0\end{array};---\left(E\right)$

⑷优化问题求解方法：

令α＝{a_0j,a_1j,L,a_nj}，则目标函数表示为多项式的一阶导数为 $g_k\left(\mathrm{\α}\right)=a_{1j}+2a_{2j}{p}_{k,j}+L+{\mathrm{na}}_{nj}{p}_{k,j}^{n-1},$ 多项式的二阶导数为 $h_k\left(\mathrm{\α}\right)=2a_{2j}+6a_{3j}{p}_{k,j}+L+n(n-1)a_{nj}{p}_{k,j}^{n-2},$ 则约束条件表示为：g_k(α)≥0，h_k(α)≥0，其中k＝1,2,L,K；问题归结为求解带不等式约束的优化问题：

$\underset{\mathrm{\α}}{min}{E}_j\left(\mathrm{\α}\right)$

s.t.g_k,j(α)≥0，

h_k,j(α)≥0，k＝1,…K

求解上述问题得到每个探测器单元对应的模体等效厚度关于多能投影数据的函数映射关系。

8.如权利要求1至7所述的基于探测器单元标定的CT图像伪影校正方法在CT图像伪影校正中的应用。

9.根据权利要求8所述的基于探测器单元标定的CT图像伪影校正方法在CT图像伪影校正中的应用，其特征在于：所述校正方法在多种不同的扫描方式中的应用；或者，所述校正方法在一类单材质占优物体的伪影校正中的应用；或者，所述校正方法在岩心三维CT成像图像伪影校正、柱状物体三维CT图像伪影校正、乳腺三维CT图像伪影校正或口腔CT图像伪影校正方面的应用。