CN102609908A - 基于基图像tv模型的ct射束硬化校正方法 - Google Patents

Abstract

本发明涉及一种CT透射成像射束硬化校正方法，特别是涉及一种基于基图像TV最小化模型的CT射束硬化校正方法。一种基于基图像TV模型的CT射束硬化校正方法，在建立带有可调参数的射束硬化校正模型的基础上，通过下述步骤完成：首先，在不同的可调参数条件下，原始投影数据经该模型预处理变换得到多组预处理投影序列，然后，分别对预处理投影序列进行重建得到一系列校正基图像，并以目标图像的全变分函数作为代价函数，通过迭代法求得加权系数最优解，最后，将得到的系列基图像加权求和，形成最终重建图像。与传统方法相比，不需要扫描件材质、射线源条件等先验知识，不受外在条件限制，通用性强；在确定加权系数时是对线性组合求解，因此具有计算复杂度小，运算时间短的优点。

Description

基于基图像TV模型的CT射束硬化校正方法

技术领域

本发明涉及一种CT透射成像射束硬化校正方法，特别是涉及一种基于基图像TV最小化模型的CT射束硬化校正方法。 

背景技术

CT(Computed Tomography，CT)技术是利用探测器对检测物体进行一系列不同角度的射线投影测量，通过三维图像重建算法直接计算出物体所有切片图像的一种成像技术，在工业无损检测和医学应用领域中有着广泛的应用。 

CT重建算法基于如下假设：所使用的X射线是单能的，射线强度随吸收物质厚度的增加按指数形式衰减，衰减系数为常数。而实际X射线源发出的射线是多能谱的，当多能谱射线与物质相互作用时，低能光子衰减量大于高能光子，因此射线平均能量随透射厚度的增加而升高，射束逐渐变“硬”，这种现象称为射束硬化。对于均匀铝柱而言，越靠近铝柱中心，射线穿过的厚度越厚，射束硬化现象越严重，射线的实际衰减量比理想衰减量低得多，这样经三维重建算法重建的铝柱图像就会出现中间暗边缘亮的现象，灰度曲线呈杯状。除了造成杯状伪影外，射束硬化还会导致CT图像中产生条状、带状等伪影， 

这些伪影的存在严重影响了图像质量，必须对其进行消除。 

目前已有许多硬化校正方法，但归纳起来主要可以分为单能法和双能法两大类。双能法中，物质总的衰减系数被分解两种能量下相应衰减系数的线性组合。这种方法的缺点是需要在不同的管电压下对被成像物采集两次投影数据，因此在实际应用中很少被采用。单能法易于实现，且校正效果良好，在实际中得到广泛应用。中国发明专利200710018779.X(公告号101126722，公告日2008年2月20日)提出了一种“基于配准模型仿真的锥束CT射束硬化校正方法”该方法首先从实际重建的零件序列切片图像中提取轮廓点集，然后采用初始-精确两步配准算法配准测量点集和CAD模型，并对配准后的CAD模型进行投影仿真，得到射线在各成像点上穿越零件的长度，最后采用新的指数函数拟合长度-灰度曲线并据此校正射束硬化伪影。该方法不足之处在于校正曲线对实验条件依赖性强，一条校正曲线只能针对一种材料在特定实验条件下使用，一旦实验条件或被测工件材料发生变化，则需要重新拟合校正曲线，当被扫描物结构或成分较为复杂时，该算法的复杂度和计算量会大大增加，影响该方法的适用范围和灵活性。 

发明内容

本发明针对现有技术不足，通过分析射束硬化对原始投影数据的影响，发现射束硬化造成原始投影中实际值低于理想值，且硬化程度越大偏离越严重的问题，提出一种基于基图像TV模型的CT射束硬化校正方法。 

本发明所采用的技术方案： 

一种基于基图像TV模型的CT射束硬化校正方法，在建立带有可调参数的射束硬化校正模型的基础上，通过下述步骤完成： 

首先，在不同的可调参数条件下，原始投影数据经该模型预处理变换得到多组预处理投影序列， 

然后，分别对预处理投影序列进行重建得到一系列校正基图像，并以目标图像的全变分函数作为代价函数，通过迭代法求得加权系数最优解， 

最后，将得到的系列基图像加权求和，形成最终重建图像。 

所述的基于基图像TV最小化模型的CT射束硬化校正方法，根据射束硬化对投影数据的影响进行分析，并在此基础上建立射束硬化校正模型： 

单能X射线穿过均匀物质的强度衰减公式，用比尔定律表示，如式(1)所示， 

I₁＝I₀e^-μx    (1) 

其中I₀表示单能射线初始强度；I₁表示透射射线强度；μ为物质衰减系数，与物质的密度和射线的能量有关；x表示射线所经过物质的长度； 

那么，由式(1)可得： 

$\mathrm{\μx}=-\ln \frac{I_1}{I_0}---\left(2\right)$

即为我们通常所说的投影数据P，(2)式变为： 

P＝μx    (3) 

由式(3)可以看出：对于单能射线，投影数据与射线穿越厚度之间呈线性关系； 

多能谱X射线情况下，式(1)变为： 

$I_1=I_0\underset{E}{\∫}S\left(E\right)\exp (-\underset{x}{\∫}\mathrm{\μ}(x^{\′},E){\mathrm{dx}}^{\′})\mathrm{dE}---\left(4\right)$

其中S(E)为入射光谱的光子密度函数，μ(x′，E)为被测截面上x点处的物质关于能量E的线性衰减系数；因此，多能谱X射线产生的投影值为 

$Q=\ln (I_1/I_0)=\ln \underset{E}{\∫}S\left(E\right)\exp (-\underset{x}{\∫}\mathrm{\μ}(x^{\′},E){\mathrm{dx}}^{\′})\mathrm{dE}---\left(5\right)$

由式(5)可知，当射线具有多能谱分布时，投影数据和射线经过物质的厚度表现为非线性关系， 

由上述分析可知，由于射束硬化的影响，透射厚度越厚，多能谱射线投影值越偏离单能射线投影值，即实际投影值越大则其偏离理想投影值越远； 

因此，对射束硬化进行校正可以定性为：将较大投影值经过校正变的更大，保持较小投影值基本不变。 

所述的基于基图像TV模型的CT射束硬化校正方法，通过下述步骤建立射束硬化校正模型： 

设P为多能谱射线条件下得到的原始投影，P′为P校正之后相应的投影数据，定义P′与P之间存在如下函数关系： 

P′＝S(P)    (6) 

其中，S为预校正函数，S是线性函数、指数函数或幂函数，结合幂函数与指数函数，设计函数S的形式如下： 

S(P)＝P    (7) 

加入不定参数n(n∈N⁺)，即S_n(P)＝Pⁿ(n∈N⁺)，利用S_n(P)拟合校正曲线，对n取不同的值则得到不同的S_n(P)，表示以不同的程度对投影值进行校正，这样会出现欠校正和过校正的情况。 

因此，赋以S_n(P)不同的权重系数c_n，来平衡欠校正和过校正的现象。S_n(P)加权求和后作为最终的预处理投影数据，即， 

$S\left(P\right)=\underset{n=1}{\mathrm{\Σ}}{c}_n{S}_n\left(P\right)=\underset{n=1}{\overset{N}{\mathrm{\Σ}}}{c}_n{P}^n---\left(8\right)$

式中系数c_n通过代价函数来确定， 

对预处理后的投影数据进行重建便得到相应的基图像，表示为： 

f_n(r)＝R^-1S_n(P)    (9) 

式中R^-1为Radon逆变换；f_n(r)是三维空间函数；r∈R³且r≡(x，y，z)，那么，最终重建图像为 

$f\left(r\right)=R^{-1}{P}^{\′}=R^{-1}S\left(P\right)=\underset{n=1}{\overset{N}{\mathrm{\Σ}}}{c}_n{f}_n\left(r\right)---\left(10\right)$

即，最终重建图像可由基图像经线性组合得到。 

所述的基于基图像TV模型的CT射束硬化校正方法，使用TV作为代价函数来评价式中最终重建图像f(r)，求解线性加权系数最优解，通过求解式： 

c＝argmin(∫G(c，r)d³r-a·c)    (11) 

来确定c，式中加入约束条件a·c＝const以避免c＝0的情况， 

令 $M_i(c,r)={\∂}_x{f}_i\left(r\right)\left({\∂}_{x}f\left(r\right)\right)+{\∂}_y{f}_i\left(r\right)\left({\∂}_{y}f\left(r\right)\right)+{\∂}_z{f}_i\left(r\right)\left({\∂}_{z}f\left(r\right)\right),$

$M_{\mathrm{ij}}(c,r)=\left({\∂}_x{f}_i\left(r\right)\right)({\∂}_x{f}_j\left(r\right)+\left({\∂}_y{f}_i\left(r\right)\right)({\∂}_y{f}_j\left(r\right)+\left({\∂}_z{f}_i\left(r\right)\right)({\∂}_z{f}_j\left(r\right),$ 则： 

$G^2(c,r)=\underset{i}{\overset{N}{\mathrm{\Σ}}}{c}_i{M}_i(c,r)=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{c}_i{c}_j{M}_{\mathrm{ij}}(c,r)---\left(12\right)$

设式最小值为0，得： 

∫G(c，r)d³r＝a·c    (13) 

将式代入得， 

$\∫G(c,r)d^{3}r=\∫\frac{\underset{i}{\overset{N}{\mathrm{\Σ}}}{c}_i{M}_i(c,r)}{G(c,r)}{d}^{3}r=\underset{i}{\overset{N}{\mathrm{\Σ}}}{c}_i\∫\frac{M_i(c,r)}{G(c,r)}{d}^{3}r=\underset{i}{\overset{N}{\mathrm{\Σ}}}{a}_i{c}_i---\left(14\right)$

$\∫G(c,r)d^{3}r=\∫\frac{\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{c}_i{c}_j{M}_{\mathrm{ij}}(c,r)}{G(c,r)}{d}^{3}r=a\·c---\left(15\right)$

由式知， 

$a_i=\∫\frac{M_i(c,r)}{G(c,r)}{d}^{3}r$

令 $B_{\mathrm{ij}}=\∫\frac{\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{c}_i{c}_j{M}_{\mathrm{ij}}(c,r)}{G(c,r)}{d}^{3}r,$ 则由式知：B＝a·c，因此，c＝B^-1a， 

使用迭代法求解线性系统c＝B^-1a，操作步骤如下： 

1)给c赋估计值 

求出c； 

2)利用约束条件a·c＝const对求出的c进行修正； 

3)将修正之后的c代入(11)式继续下一次迭代； 

4)重复步骤2、3，不断修正估计值，直到c值达到收敛条件。 

本发明的有益积极效果： 

1、本发明基于基图像TV模型的CT射束硬化校正方法，建立了射束硬化校正模型，直接对原始投影数据进行运算，与传统方法相比，不需要扫描件材质、射线源条件等先验知识，不受外在条件限制，通用性强；在确定加权系数时是对线性组合求解，因此具有计算复杂度小，运算时间短的优点；以代价函数最小化的方法求取加权系数，可以使算法在抑制伪影的同时尽可能降低噪声。 

2、本发明基于基图像TV模型的CT射束硬化校正方法，能够有效抑制射束硬化造成的杯状和条纹伪影，改善图像质量，从而提高信息识别和提取的准确率。实验结果表明，该算法对真实的铝柱体模和工业检测件CT图像中的杯状和条状伪影均有明显的抑制作用。 

附图说明

图1：投影数据与透射厚度关系的影响分析示意图； 

图2：采用本发明校正方法校正过程的中间结果：铝柱校正基图像； 

图3(a)～图3(c)分别为采用本发明校正方法得到的校正基图像：铝柱硬化校正效果图； 

图4：采用本发明校正方法校正前后的CT图像：光缆硬化校正效果对比示意图。 

具体实施方式

实施例一：CT成像系统中，X射束硬化导致重建图像中出现杯状和条状伪影，严重影响了图像质量。本文通过分析射束硬化对原始投影数据的影响，发现射束硬化造成原始投影中实际值低于理想值，且硬化程度越大偏离越严重。本发明基于基图像TV模型的CT射束硬化校正方法，在建立带有可调参数的射束硬化校正模型的基础上，通过下述步骤完成：首先，在不同的可调参数条件下，原始投影数据经该模型预处理变换得到多组预处理投影序列，然后，分别对预处理投影序列进行重建得到一系列校正基图像，并以目标图像的全变分函数作为代价函数，通过迭代法求得加权系数最优解，最后，将得到的系列基图像加权求和，形成最终重建图像。 

实验结果表明，该算法对真实的铝柱体模和工业检测件CT图像中的杯状和条状伪影均有明显的抑制作用。 

实施例二：本实施例的基于基图像TV模型的CT射束硬化校正方法，首先根据射束硬化对投影数据的影响进行分析，然后在此基础上建立射束硬化校正模型，采用下述方法进行射束硬化对投影数据的影响分析： 

单能X射线穿过均匀物质的强度衰减公式可用比尔定律表示，如式(1)所示， 

I₁＝I₀e^-μx    (1) 

其中I₀表示单能射线初始强度；I₁表示透射射线强度；μ为物质衰减系数，与物质的密度和 射线的能量有关；x表示射线所经过物质的长度。那么，由式(1)可得： 

$\mathrm{\μx}=-\ln \frac{I_1}{I_0}---\left(2\right)$

即为我们通常所说的投影数据P，(2)式变为： 

P＝μx    (3) 

由式(3)可以看出：对于单能射线，投影数据与射线穿越厚度之间呈线性关系。 

多能谱X射线情况下，式(1)变为： 

$I_1=I_0\underset{E}{\∫}S\left(E\right)\exp (-\underset{x}{\∫}\mathrm{\μ}(x^{\′},E){\mathrm{dx}}^{\′})\mathrm{dE}---\left(4\right)$

其中S(E)为入射光谱的光子密度函数，μ(x′，E)为被测截面上x点处的物质关于能量E的线性衰减系数。因此，多能谱X射线产生的投影值为 

$Q=\ln (I_1/I_0)=\ln \underset{E}{\∫}S\left(E\right)\exp (-\underset{x}{\∫}\mathrm{\μ}(x^{\′},E){\mathrm{dx}}^{\′})\mathrm{dE}---\left(5\right)$

由式(5)可知，当射线具有多能谱分布时，投影数据和射线经过物质的厚度表现为非线性关系，如附图1投影数据与透射厚度关系所示。由前述分析过程及附图1可知，由于射束硬化的影响，透射厚度越厚，多能谱射线投影值越偏离单能射线投影值。即实际投影值越大则其偏离理想投影值越远。因此，对射束硬化进行校正可以定性为：将较大投影值经过校正变的更大，保持较小投影值基本不变。射束硬化对投影值影响的定量研究是复杂的。蒙卡模拟可以做到定量估，但需要对扫描件中所包含的所有材质进行厚度尺寸的测量，而这种测量是难以达到精确的，而且蒙卡模拟耗时较长，实用难度大。 

实施例三：本实施例的基于基图像TV模型的CT射束硬化校正方法，在实施例二的射束硬化对投影数据的影响进行分析的基础上，通过下述步骤建立硬化校正模型： 

设P为多能谱射线条件下得到的原始投影，P′为P校正之后相应的投影数据，定义P′与P之间存在如下函数关系： 

P′＝S(P)    (6) 

其中，S为预校正函数，S是线性函数、指数函数或幂函数，结合幂函数与指数函数，设计函数S的形式如下： 

S(P)＝P    (7) 

加入不定参数n(n ∈N⁺)，即S_n(P)＝Pⁿ(n∈N⁺)，利用S_n(P)拟合校正曲线，对n取不同 的值则得到不同的S_n(P)，表示以不同的程度对投影值进行校正，这样会出现欠校正和过校正的情况。 

因此，赋以S_n(P)不同的权重系数c_n，来平衡欠校正和过校正的现象。S_n(P)加权求和后作为最终的预处理投影数据，即， 

$S\left(P\right)=\underset{n=1}{\mathrm{\Σ}}{c}_n{S}_n\left(P\right)=\underset{n=1}{\overset{N}{\mathrm{\Σ}}}{c}_n{P}^n---\left(8\right)$

式中系数c_n通过代价函数来确定， 

对预处理后的投影数据进行重建便得到相应的基图像，表示为： 

f_n(r)＝R^-1S_n(P)    (9) 

式中R^-1为Radon逆变换；f_n(r)是三维空间函数；r∈R³且r≡(x，y，z)。那么，最终重建图像为 

$f\left(r\right)=R^{-1}{P}^{\′}=R^{-1}S\left(P\right)=\underset{n=1}{\overset{N}{\mathrm{\Σ}}}{c}_n{f}_n\left(r\right)---\left(10\right)$

即，最终重建图像可由基图像经线性组合得到。 

通过调整加权系数c_n可以降低最终重建图像中的噪声水平。 

实施例四：本实施例的基于基图像TV模型的CT射束硬化校正方法，与实施例三不同的是，使用TV作为代价函数来评价式中最终重建图像f(r)，求解线性加权系数最优解。 

TV是测量图像平滑度的一项重要指标。使图像的TV最小则可以抑制噪声。射束硬化引起的伪影给CT图像引入了大量噪声，降低了图像平滑度。因此，本文使用TV作为目标函数来评价最终重建图像f(r)，能使f(r)的TV最小的c作为最优解。即通过求解式： 

c＝argmin(∫G(c，r)d³r-a·c)    (11) 

来确定c。式中加入约束条件a·c＝const以避免c＝0的情况。 

令 $M_i(c,r)={\∂}_x{f}_i\left(r\right)\left({\∂}_{x}f\left(r\right)\right)+{\∂}_y{f}_i\left(r\right)\left({\∂}_{y}f\left(r\right)\right)+{\∂}_z{f}_i\left(r\right)\left({\∂}_{z}f\left(r\right)\right),$

$M_{\mathrm{ij}}(c,r)=\left({\∂}_x{f}_i\left(r\right)\right)({\∂}_x{f}_j\left(r\right)+\left({\∂}_y{f}_i\left(r\right)\right)({\∂}_y{f}_j\left(r\right)+\left({\∂}_z{f}_i\left(r\right)\right)({\∂}_z{f}_j\left(r\right),$ 则： 

$G^2(c,r)=\underset{i}{\overset{N}{\mathrm{\Σ}}}{c}_i{M}_i(c,r)=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{c}_i{c}_j{M}_{\mathrm{ij}}(c,r)---\left(12\right)$

设式最小值为0，得： 

∫G(c，r)d³r＝a·c    (13) 

将式代入得， 

$\∫G(c,r)d^{3}r=\∫\frac{\underset{i}{\overset{N}{\mathrm{\Σ}}}{c}_i{M}_i(c,r)}{G(c,r)}{d}^{3}r=\underset{i}{\overset{N}{\mathrm{\Σ}}}{c}_i\∫\frac{M_i(c,r)}{G(c,r)}{d}^{3}r=\underset{i}{\overset{N}{\mathrm{\Σ}}}{a}_i{c}_i---\left(14\right)$

$\∫G(c,r)d^{3}r=\∫\frac{\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{c}_i{c}_j{M}_{\mathrm{ij}}(c,r)}{G(c,r)}{d}^{3}r=a\·c---\left(15\right)$

由式知， 

$a_i=\∫\frac{M_i(c,r)}{G(c,r)}{d}^{3}r$

令 $B_{\mathrm{ij}}=\∫\frac{\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{c}_i{c}_j{M}_{\mathrm{ij}}(c,r)}{G(c,r)}{d}^{3}r,$

则由式知：B＝a·c。因此，c＝B^-1a。使用迭代法求解线性系统c＝B^-1a。首先给c赋一估计值 

求出c后，利用约束条件a·c＝const对求出的c进行修正，然后修正之后的c代入式继续下一次迭代。重复上述步骤，不断修正估计值，直到c值达到收敛条件。 

我们根据本文提出的方法，开展了如下实验。实验结果表明，该算法对真实的铝柱体模和工业检测件CT图像中的杯状和条状伪影均有明显的抑制作用。 

实验平台为：Thales Hawkeye130微焦斑X射线源和Thales Pixium RF 4343平板探测器所构建的工业CT平台。其中，射线源小焦斑模式可达到7μm，中焦斑模式20μm；探测器尺寸为43cm*43cm，分辨率为0.148mm。 

首先，本文选用直径2.5cm的材质相对均匀的铝柱作为实验对象以验证该方法对典型的杯状伪影的校正效果。 

实验条件如下：射线源选用中焦斑模式，管电压120kV、管电流200μa，探测器采集帧率3fps，360度圆周扫描，采集360帧投影图像。 

设N＝5，附图2是本算法校正过程的中间结果，取重建三维数据的正中间切片，且都在相同的窗宽窗位水平下显示。图3(a)～图3(e)是校正基图像；其中，图3(a)是n＝1时的基图像，也是原始数据直接重建图像，图中有明显的杯状伪影。利用这些基图像，根据上文中所述的求解加权系数方法可以求得校正系数：c＝(0.02420.82280.06940.02170.0329)’。 

附图3(a)是最终校正图像，图中几乎没有杯状伪影。为了对比校正前后的效果，附图3(b)展示了原始图像，并在附图3(c)中利用MATLAB工具画出附图3(a)和附图3(b)中间行的图像灰度值曲线。从校正前后的对比图中我们可以看出，杯状伪影得到很好的抑制。 

另外，本文选用了工业检测中一实物——光缆作为扫描对象来验证本文算法抑制条状伪影的能力。在对光缆进行CT检测时，检测目标——光纤的信号本身就很微弱，加之硬化伪影对其影响，会增大信号辨识难度，影响检测效率。图4是校正前后CT图像。图4(a)是未进行射束硬化校正直接重建得到的CT图像，从图中我们可以看到，射束硬化造成明显的明暗条纹，尤其是图中的暗条纹，直接覆盖了光纤信息(箭头所指)；图4(b)是使用本文方法进行射束硬化校正后重建得到的CT图像，从图中可以看到，校正之后的光缆图像中，条纹伪影得到一定程度的抑制，较好地恢复出了光纤信息；图4(c)是图4(a)、图4(b)中间列的图像灰度值曲线，由校正前后曲线对比可以看出，校正之后杯状伪影也得到较好的校正，图像边缘对比度有所提高。 

Claims (4)

1.一种基于基图像TV模型的CT射束硬化校正方法，其特征是：在建立带有可调参数的射束硬化校正模型的基础上，通过下述步骤完成：

首先，在不同的可调参数条件下，原始投影数据经该模型预处理变换得到多组预处理投影序列；

然后，分别对预处理投影序列进行重建得到一系列校正基图像，并以目标图像的全变分函数作为代价函数，通过迭代法求得加权系数最优解；

最后，将得到的系列基图像加权求和，形成最终重建图像。

2.根据权利要求1所述的基于基图像TV模型的CT射束硬化校正方法，其特征是：根据射束硬化对投影数据的影响进行分析，并在此基础上建立射束硬化校正模型：

单能X射线穿过均匀物质的强度衰减公式，用比尔定律表示，如式(1)所示，

I₁＝I₀e^-μx    (1)

其中I₀表示单能射线初始强度；I₁表示透射射线强度；μ为物质衰减系数，与物质的密度和射线的能量有关；x表示射线所经过物质的长度；

那么，由式(1)可得：

$\mathrm{\μx}=-\ln \frac{I_1}{I_0}---\left(2\right)$

即为我们通常所说的投影数据P，(2)式变为：

P＝μx    (3)

由式(3)可以看出：对于单能射线，投影数据与射线穿越厚度之间呈线性关系；

多能谱X射线情况下，式(1)变为：

$I_1=I_0\underset{E}{\∫}S\left(E\right)\exp (-\underset{x}{\∫}\mathrm{\μ}(x^{\′},E){\mathrm{dx}}^{\′})\mathrm{dE}---\left(4\right)$

其中S(E)为入射光谱的光子密度函数，μ(x′，E)为被测截面上x点处的物质关于能量E的线性衰减系数；因此，多能谱X射线产生的投影值为

$Q=\ln (I_1/I_0)=\ln \underset{E}{\∫}S\left(E\right)\exp (-\underset{x}{\∫}\mathrm{\μ}(x^{\′},E){\mathrm{dx}}^{\′})\mathrm{dE}---\left(5\right)$

由式(5)可知，当射线具有多能谱分布时，投影数据和射线经过物质的厚度表现为非线性关系，

由上述分析可知，由于射束硬化的影响，透射厚度越厚，多能谱射线投影值越偏离单能射线投影值，即实际投影值越大则其偏离理想投影值越远；因此，对射束硬化进行校正定性：将较大投影值经过校正变的更大，保持较小投影值基本不变。

3.根据权利要求1或2所述的基于基图像TV模型的CT射束硬化校正方法，其特征是：通过下述步骤建立射束硬化校正模型：

设P为多能谱射线条件下得到的原始投影，P′为P校正之后相应的投影数据，定义P′与P之间存在如下函数关系：

P′＝S(P)    (6)

其中，S为预校正函数，S是线性函数、指数函数或幂函数，结合幂函数与指数函数，设计函数S的形式如下：

S(P)＝P    (7)

加入不定参数n(n∈N⁺)，即S_n(P)＝Pⁿ(n∈N⁺)，利用S_n(P)拟合校正曲线，对n取不同的值则得到不同的S_n(P)，表示以不同的程度对投影值进行校正，

针对出现欠校正和过校正的情况，赋以S_n(P)不同的权重系数c_n，来平衡欠校正和过校正的现象；S_n(P)加权求和后作为最终的预处理投影数据，即，

$S\left(P\right)=\underset{n=1}{\mathrm{\Σ}}{c}_n{S}_n\left(P\right)=\underset{n=1}{\overset{N}{\mathrm{\Σ}}}{c}_n{P}^n---\left(8\right)$

式中系数c_n通过代价函数来确定，

对预处理后的投影数据进行重建便得到相应的基图像，表示为：

f_n(r)＝R^-1S_n(P)    (9)

式中R^-1为Radon逆变换；f_n(r)是三维空间函数；r∈R³且r≡(x，y，z)，那么，最终重建图像为

$f\left(r\right)=R^{-1}{P}^{\′}=R^{-1}S\left(P\right)=\underset{n=1}{\overset{N}{\mathrm{\Σ}}}{c}_n{f}_n\left(r\right)---\left(10\right)$

即，最终重建图像可由基图像经线性组合得到。

4.根据权利要求3所述的基于基图像TV模型的CT射束硬化校正方法，其特征是：使用TV作为代价函数来评价式中最终重建图像f(r)，求解线性加权系数最优解，通过求解式：

c＝argmin(∫G(c，r)d³r-a·c)    (11)

来确定c，式中加入约束条件a·c＝const以避免c＝0的情况，

令 $M_i(c,r)={\∂}_x{f}_i\left(r\right)\left({\∂}_{x}f\left(r\right)\right)+{\∂}_y{f}_i\left(r\right)\left({\∂}_{y}f\left(r\right)\right)+{\∂}_z{f}_i\left(r\right)\left({\∂}_{z}f\left(r\right)\right),$

$M_{\mathrm{ij}}(c,r)=\left({\∂}_x{f}_i\left(r\right)\right)({\∂}_x{f}_j\left(r\right)+\left({\∂}_y{f}_i\left(r\right)\right)({\∂}_y{f}_j\left(r\right)+\left({\∂}_z{f}_i\left(r\right)\right)({\∂}_z{f}_j\left(r\right),$ 则：

$G^2(c,r)=\underset{i}{\overset{N}{\mathrm{\Σ}}}{c}_i{M}_i(c,r)=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{c}_i{c}_j{M}_{\mathrm{ij}}(c,r)---\left(12\right)$

设式最小值为0，得：

∫G(c，r)d³r＝a·c    (13)

将式代入得，

$\∫G(c,r)d^{3}r=\∫\frac{\underset{i}{\overset{N}{\mathrm{\Σ}}}{c}_i{M}_i(c,r)}{G(c,r)}{d}^{3}r=\underset{i}{\overset{N}{\mathrm{\Σ}}}{c}_i\∫\frac{M_i(c,r)}{G(c,r)}{d}^{3}r=\underset{i}{\overset{N}{\mathrm{\Σ}}}{a}_i{c}_i---\left(14\right)$

$\∫G(c,r)d^{3}r=\∫\frac{\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{c}_i{c}_j{M}_{\mathrm{ij}}(c,r)}{G(c,r)}{d}^{3}r=a\·c---\left(15\right)$

由式知，

$a_i=\∫\frac{M_i(c,r)}{G(c,r)}{d}^{3}r$

令 $B_{\mathrm{ij}}=\∫\frac{\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{c}_i{c}_j{M}_{\mathrm{ij}}(c,r)}{G(c,r)}{d}^{3}r,$ 则由式知：B＝a·c，因此，c＝B^-1a，

使用迭代法求解线性系统c＝B^-1a，操作步骤如下：

1)给c赋估计值求出c；

2)利用约束条件a·c＝const对求出的c进行修正；

3)将修正之后的c代入(11)式继续下一次迭代；

4)重复步骤2)、3)，不断修正估计值，直到c值达到收敛条件。

Images