CN101226642A - 基于ct数据一致性的投影射束硬化校正方法 - Google Patents

Abstract

本发明公开了一种基于数据相关性的CT投影数据射束硬化校正方法，该方法以医用诊断X线物理学成像模型为基础；根据CT数据一致性，构造CT投影数据之间极小化求解问题；结合上述极小化问题，得到中间参数的求解方法；为了实现更加精确的校正结果，根据初次校正重建结果和CT重投影，估计投影数据中的高密度物质的衰减比例。本发明适合于各类X线CT设备的射束硬化校正功能的实现。和现有的只采用水模校正的方法相比较，能够分别考虑不同密度组织的射束硬化效应的差异；和现有的骨校正的方法相比较，能够更为灵活地自动适应成像对象的不同，具有更好的校正精度。

Description

基于CT数据一致性的投影射束硬化校正方法

技术领域

本发明属于X线CT重建成像技术领域，涉及基于CT数据一致性的投影数据射束硬化校正方法。

背景技术

CT重建原理为Radon变换与反变换理论，其基础是投影数据来自于线性投影模型，比如对于二维断层成像来讲，有如下的投影公式：

p(θ，l)＝∫_L(θ，l)f(x，y)dxdy    (1)

其中θ为投影角度，l为投影数据的采样点，L(θ，l)为射线穿透路径，f(x，y)为图像函数。

虽然CT技术已经从第一代发展到当今的第五代，扫描轨迹从最早的直线，一直到现在的圆形和螺旋轨迹、乃至任意轨迹，投影方式有扇形和锥形束，其投影和重建理论一直基于线性投影模型。所进行的重建理论研究总是围绕着不同的扫描轨迹和投影方式进行的。

理论上，如果CT的X线源是单能的，那么上述的线性投影模型则是满足的，这是因为物质对能量为E₀的X线光子的衰减成指数衰减形式，取对数后即可满足所谓的线性投影模型，如下式：

I(θ，l)＝I₀exp{-∫_L(θ，l)μ(x，y，E₀)ρ(x，y)dxdy}    (2)

上式中，I₀为球管发出X线光子强度，μ(x，y，E₀)为成像对象的质量衰减系数，ρ(x，y)为成像对象体密度。

事实上，现代CT球管发出的X线光谱是具有一定能谱分布的宽谱函数。以头部CT常用的120kVp电压下球管发出的X线光谱为例，其X线光子能量从30keV到120keV按照一定的规律连续分布，而且其分布函数在58-68keV处还存在4根不连续的特征辐射光谱。如图(1)所示，图中上部图形是CT球管发出的X线光谱函数，下部图形是经过35mm铝衰减后的光谱函数(放大了5倍)。宽能谱条件下所形成的投影数据如下：

$I(\mathrm{\θ},l)=I_0{\∫}_0^{V}P(V,E)\exp \{-{\∫}_{L(\mathrm{\θ},l)}\mathrm{\μ}(x,y,E)\mathrm{\ρ}(x,y)\mathrm{dxdy}\}\mathrm{dE}---\left(3\right)$

其中，P(V，E)为球管电压取VkVp时，发出的X线光谱乘以探测器吸收谱的归一化函数，μ(x，y，E)仍为成像对象的质量衰减系数，但是考虑到该参数随X线光子能量E变化的关系。上式取对数后，不再满足上述的线性投影模型。以人体常见组织如骨皮质、肌肉和脂肪为例，不同厚度组织在120kVp电压下，对球管发出的宽谱X线光子的衰减数据取对数后函数关系如图(2)所示(假设X线探测器是CsI晶体，下同)。图(2)为骨组织、肌肉和脂肪在120kVp电压球管发出的光谱中衰减量和厚度的关系，其中衰减函数与对应线性关系的偏差已经放大了10倍。

线性投影模型要求上述函数是直线形式，但实际上该函数并不满足其要求。这样直接运用重建理论得到的结果是：对于同样的人体组织，不同空间区域的重建数值(即CT值)会出现不同，由于CT设备提供主要的信息反映在CT值上，因此这会给诊断带来至关重要的误差，这就是所谓的“射束硬化效应”，所形成的CT值误差被称为“射束硬化伪影”，是CT设备的主要重建伪影之一。

现有CT设备中，实现的所谓射束硬化校正方法多来自于工程实践。观察图(2)，在宽谱条件下，虽然投影数据与组织厚度不是线性关系，但与直线相差不大，因此可以考虑采用一个多项式来描述该函数，并且根据该多项式系数来实现线性化校正。为了减少病人辐射剂量，减轻射束硬化伪影的强度，CT设备在球管出线处加上一个不均匀厚度的Bowtie补偿滤波器，将光谱中一部份低能光子滤除，使穿透人体的X线投影数据的线性度尽可能地好些，再进行多项式校正，因此得到现有CT设备的射束硬化校正措施1。

校正措施1：将投影数据取对数后，利用一个已知的多项式进行线性化校正，而该多项式系数是通过事先的水模实验来获得，由于存在球管老化等问题需要定期校正。该校正措施也被称为水模校正或者水校正。

注意到图(2)所示，不同的人体组织具有不同的校正多项式系数，而在投影数据校正时，并不知道人体中不同组织成份的多少。为了避免人体不同结构带来的校正误差，现有的CT设备采用了校正措施2。

校正措施2：首先，采用水校正对取对数后的投影数据进行射束硬化预校正；然后，对投影数据进行CT重建得到一个预重建图像；对预重建图像进行简单的阈值分割，将整幅图像分割为软组织区域和骨组织区域；分别对分割后软组织区域和骨组织区域图像进行重投影，投影几何与采集几何相同；经验选择一个尺度因子，由重投影数据估计软组织面密度和骨组织面密度；根据已有能谱分布函数、软组织与骨组织的典型质量衰减系数、水校正多项式、软组织面密度和骨组织面密度，估计射束硬化预校正的偏差量；将该偏差量反馈到射束硬化预校正投影数据中，完成射束硬化校正。该校正措施也被称为骨校正。

事实上，现有医学CT的射束硬化校正方法的是上述两种措施的组合。显然，这样的方法存在以下问题。

现有CT设备射束硬化校正方法存在的缺点：

1.每一台CT需要进行定期的水模实验，以获取校正多项式系数，增加使用难度和成本。

2.为了适应人体各部分的差异，需要按照部位分别设计模体进行实验，来获得不同的校正多项式系数。

3.对于同一人体部位还存在个体差异，因此用统一系数的多项式校正可能会存在问题。这里最典型的是头部CT中的硬膜外伪影，目前尚无理想解决方案，在大部分的CT设备中，需要靠医生的临床经验来判断区分，而在其它的一些CT设备中，有采用后处理方法校正的方案，但存在重建数据不真实的风险。

4.为了克服以上问题而采用骨校正，虽然具有一定的自适应性，但是也存在一些重要的缺点，如：经验选择尺度因子具有很大难度等，而尺度因子的选择是否合适又对整体射束硬化校正性能有重要影响。

CT设备的核心虽然是重建算法，但为了实现准确的重建效果，众多的校正算法一直是各类CT设备关键核心技术。除了射束硬化校正，还有各种基准校正、探测器通道非线性校正等等。在这些校正中射束硬化校正与扫描的成像对象有关，其它则大多只与设备相关。与设备相关的校正无论复杂与否，均可以通过一个特定的程序精确地进行。而与成像对象有关的射束硬化校正，如果只采用模体校正方法，那么在理论上不可能做到精确。而骨校正方法虽然一定程度上考虑了成像对象，但是终归只是一种工程化的解决问题方案，缺乏严格的理论支持。

申请人在对CT设备射束硬化校正的深入研究过程中，注意到这样的事实：在不同角度的投影成像中，成像对象中的每一个坐标点被多个射线穿过，也就是说各个不同角度的投影数据均包含了该点的信息，这意味着不同角度的投影数据是相关的。如果存在一个具有有限参数的校正方程，那么可能可以通过这种相关性解出校正方程的系数，从而实现校正的目标，而且这样的校正来自于成像对象的投影数据本身，具有自适应性，理论上会有更好的精度，有可能避免上述缺点。

在申请人己申请的专利中，详细描述了平行束CT数据一致性、二维扇形束CT数据一致性、以及在此基础上求解预先设定的非线性方程参数从而完成线性校正的方法。

另外现有研究结果表明，CT投影数据校正可以由多项式完成，理论上该方法可以用于宽谱XCT的投影数据校正。这些工作给基于自身一致性进行CT投影数据射束硬化校正研究工作，提供了可靠的理论基础。

发明内容

本发明的目的在于，提供一种基于CT数据一致性的投影数据射束硬化效应校正方法，该方法适用于各类CT设备重建的数据校正。

实现上述发明目的的技术解决方案是：

一种基于CT数据一致性的CT投影数据射束硬化效应校正方法，其特征在于，该方法根据CT数据一致性，对不同角度的投影数据进行射束硬化效应校正，具体包括以下步骤：

步骤1：水模校正

在确定线性化多项式阶段，直接将水模体作为成像对象，获取不同厚度模体所对应的宽能谱投影值，然后拟和计算多项式系数，得到如下表达式：

$g_i=a_0^i+a_1^i{\hat{g}}_i+a_2^i{\hat{g}}_i^2+L+a_N^i{\hat{g}}_i^N,---\left(1\right)$

其中，g_i和

分别为单能谱估计值和宽能谱投影值。

本发明文档中，上下标i用于区分不同X线投影路径。

在校正阶段，根据已经确定好的多项式系数α_n ⁱ(n＝0，1，2，L，N)和现场采集到的投影数据

计算单能谱投影值，即：

$g_i^f=a_0^i+a_1^i{\hat{g}}_i^f+a_2^i{\left({\hat{g}}_i^f\right)}^2+L+a_N^i{\left({\hat{g}}_i^f\right)}^N=P_i^w\left[{\hat{g}}_i^f\right].---\left(2\right)$

最后采用标准CT重建算法处理投影数据集{g_i ^f}。

该校正方法的有效性依赖于成像对象X线衰减特性的水等效假设。正常人体部位的组织可大致分为软组织和骨组织，软组织基本满足水等效假设，而骨组织与其有一定偏差。正因为如此造成水模校正一般仅对杯形伪影有效，而对条纹伪影无能为力。虽然如此仍可将水模校正作为预校正方法，提升整体校正效果。

步骤2：预分割、重投影与面密度比例计算

水模校正、预重建后得到校正图像f^w(x，y)。通过阈值分割得到软组织区域图像

$f^s(x,y)=\left\{\begin{array}{cc}{f}^w(x,y),& f^w(x,y)<T\\ 0,& f^w(x,y)\&GreaterEqual;T\end{array}\right.,---\left(3\right)$

和骨组织区域图像

$f^b(x,y)=\left\{\begin{array}{cc}0,& f^w(x,y)<T\\ f^w(x,y),& f^w(x,y)\&GreaterEqual;T\end{array}\right.,---\left(4\right)$

其中阈值T的选择采用简单统计方式，该发明整体性能对该值的选取不敏感。在重投影之前需要为骨组织区域图像f^b(x，y)选择一个尺度因子λ₀，从而有

$f^b(x,y)\left\{\begin{array}{cc}0,& f^w(x,y)<T\\ f^w(x,y)/{\mathrm{\&lambda;}}_0,& f^w(x,y)\&GreaterEqual;T\end{array}\right..---\left(5\right)$

该尺度因子主要是用来调整骨组织区域重建值与体密度之间的映射关系，采用经验选择方式。

重投影采用与现场CT投影数据采集相同的几何模式。预分割和重投影占用较多运算时间和存储空间。重投影完成后得到第i条X线对应的软组织面密度ρ_i ^a，s和骨组织面密度ρ_i ^a，b，而面密度比例ri按如下公式计算

${\mathrm{\&gamma;}}_i=\left\{\begin{array}{cc}{\mathrm{\&rho;}}_i^{a,b}/({\mathrm{\&rho;}}_i^{a,s}+{\mathrm{\&rho;}}_i^{a,b}),& {\mathrm{\&rho;}}_i^{a,s}+{\mathrm{\&rho;}}_i^{a,b}\&NotEqual;0\\ 0,& {\mathrm{\&rho;}}_i^{a,s}+{\mathrm{\&rho;}}_i^{a,b}=0\end{array}\right..---\left(6\right)$

步骤3：H-L射束硬化校正方法

步骤3.1：General H-L校正方法

步骤3.1.1：训练数据采集与计算

本发明方法需要对一些训练数据进行采集和计算，以便得到单能谱投影与宽能谱投影之间的约束关系。在医学诊断X线谱段，由于X线光子与类似于人体组织的物质的相互作用主要表现为光电吸收和康普顿散射两种效应，所以宽能谱X线成像的投影数据可以进一步近似表示为

$I_i\&ap;I_0{\&Integral;}_0^V{P}_i(V,E)\exp \{-{\&Integral;}_{L_i}[\mathrm{bd}\left(l\right)\frac{12.4}{E^3}+\mathrm{bp}\left(l\right)f_{\mathrm{KN}}\left(E\right)\left]\mathrm{\&rho;}\left(l\right)\mathrm{dl}\right\}\mathrm{dE}.---\left(7\right)$

在采集训练数据过程中，仅考虑两种材料，它们的X线衰减特性分别与人体软组织和骨组织近似(可使用有机玻璃和纯铝)。因为不会造成混淆，在此仍以上标s与b区分。

则有：

$I_i$

$\&ap;I_0{\&Integral;}_0^V{P}_i(V,E)\exp \{-{\mathrm{\&rho;}}_i^s{L}_i^s[{\mathrm{bd}}_i^s\frac{12.4}{E^3}+{\mathrm{bp}}_i^s{f}_{\mathrm{KN}}\left(E\right)]-{\mathrm{\&rho;}}_i^b{L}_i^b[{\mathrm{bd}}_i^b\frac{12.4}{E^3}+{\mathrm{bp}}_i^b{f}_{\mathrm{KN}}\left(E\right)\left]\right\}\mathrm{dE}$

$=I_0{\&Integral;}_0^V{P}_i(V,E)\exp \{-{\mathrm{\&rho;}}_i^a(1-{\mathrm{\&gamma;}}_i)[{\mathrm{bd}}_i^s\frac{12.4}{E^3}+{\mathrm{bp}}_i^s{f}_{\mathrm{KN}}\left(E\right)]-{\mathrm{\&rho;}}_i^a{\mathrm{\&gamma;}}_i[{\mathrm{bd}}_i^b\frac{12.4}{E^3}+{\mathrm{bp}}_i^b{f}_{\mathrm{KN}}\left(E\right)\left]\right\}\mathrm{dE}$

$=I_0{\&Integral;}_0^V{P}_i(V,E)\exp \{-{\mathrm{\&rho;}}_i^a\{[(1-{\mathrm{\&gamma;}}_i){\mathrm{bd}}_i^s+{\mathrm{\&gamma;}}_i{\mathrm{bd}}_i^b]\frac{12.4}{E^3}+[(1-{\mathrm{\&gamma;}}_i){\mathrm{bp}}_i^s+{\mathrm{\&gamma;}}_i{\mathrm{bp}}_i^b]f_{\mathrm{KN}}\left(E\right)\left\}\right\}\mathrm{dE}$

$=I_0{\&Integral;}_0^V{P}_i(V,E)\exp \{-{\mathrm{\&rho;}}_i^a[{\mathrm{bd}}_i\frac{12.4}{E^3}+{\mathrm{bp}}_i{f}_{\mathrm{KN}}\left(E\right)\left]\right\}\mathrm{dE}$

其中：

${\mathrm{\&rho;}}_i^a={\mathrm{\&rho;}}_i^s{L}_i^s+{\mathrm{\&rho;}}_i^b{L}_i^b---\left(8\right)$

${\mathrm{\&gamma;}}_i={\mathrm{\&rho;}}_i^b{L}_i^b/({\mathrm{\&rho;}}_i^s{L}_i^s+{\mathrm{\&rho;}}_i^b{L}_i^b)---\left(9\right)$

这样

$\frac{{\mathrm{bp}}_i}{{\mathrm{bd}}_i}$

$=\frac{(1-{\mathrm{\&gamma;}}_i){\mathrm{bp}}_i^s+{\mathrm{\&gamma;}}_i{\mathrm{bp}}_i^b}{(1-{\mathrm{\&gamma;}}_i){\mathrm{bd}}_i^s+{\mathrm{\&gamma;}}_i{\mathrm{bd}}_i^b}---\left(10\right)$

$\&ap;\frac{(1-{\mathrm{\&gamma;}}_i)k{Z}_s^{3}b{d}_i^s+{\mathrm{\&gamma;}}_{i}k{Z}_b^3{\mathrm{bd}}_i^b}{(1-{\mathrm{\&gamma;}}_i){\mathrm{bd}}_i^s+{\mathrm{\&gamma;}}_i{\mathrm{bd}}_i^b}$

$\&ap;k[(1-{\mathrm{\&gamma;}}_i)Z_s^3+{\mathrm{\&gamma;}}_i{Z}_b^3]$

所以大致有如下比例关系

$\frac{{\mathrm{bp}}_i}{{\mathrm{bd}}_i}\&Proportional;Z_{\mathrm{eff},i}^3---\left(11\right)$

其中，Z_s、Z_b与Z_eff，i分别为软组织近似材料原子序数、骨组织近似材料原子序数与第i条X线经过所有材料的等效原子序数。这样用于训练的单能谱投影与宽能谱投影之间的约束关系可以选为：

$g_i=\underset{m=1}{\overset{\max m}{\mathrm{\&Sigma;}}}\underset{n=1}{\overset{\max n}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_{\mathrm{imn}}{\left(Z_{\mathrm{eff},i}^3\right)}^{n-1}{\left({\hat{g}}_i\right)}^m---\left(12\right)$

其中，g_i与Z_eff，i计算得到，而

从采集的训练数据I_i计算得到，从而可以确定出系数α_imn，将之存储下来用于后续数据处理。对于训练数据的采集，可以不强调几何模式。

步骤3.1.2：γ与Z³之间的映射函数逼近

对于X线成像现场数据，虽然缺乏较精确计算Z_eff，i的手段，但是可以令

$Z_{\mathrm{eff},i}^3=\underset{r=1}{\overset{\max r}{\mathrm{\&Sigma;}}}{\mathrm{\&beta;}}_r{\mathrm{\&gamma;}}_i^{r-1}---\left(13\right)$

则有

$g_i^f=\underset{m=1}{\overset{\max m}{\mathrm{\&Sigma;}}}\underset{n=1}{\overset{\max n}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_{\mathrm{imn}}{\left(\underset{r=1}{\overset{\max r}{\mathrm{\&Sigma;}}}{\mathrm{\&beta;}}_r{\mathrm{\&gamma;}}_i^{r-1}\right)}^{n-1}{\left({\hat{g}}_i^f\right)}^m---\left(14\right)$

其中，α_imn在训练阶段(步骤3.1.1)已确定，γ_i通过水模校正、预分割、重投影、面密度比例计算等一系列步骤得到，由X线成像现场数据计算得到，g_i ^f虽未知但其必满足发散束CT数据一致性条件。通过所有这些已知条件，可以反解出β_r，实现γ与Z³之间的映射函数逼近。

反解β_r的方法如下：

取maxn＝3，则有

$g_i^f=\underset{m=1}{\overset{\max m}{\mathrm{\&Sigma;}}}\underset{n=1}{\overset{3}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_{\mathrm{imn}}{\left(\underset{r=1}{\overset{\max r}{\mathrm{\&Sigma;}}}{\mathrm{\&beta;}}_r{\mathrm{\&gamma;}}_i^{r-1}\right)}^{n-1}{\left({\hat{g}}_i^f\right)}^m$

$=\underset{m=1}{\overset{\max m}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_{\mathrm{im}1}{\left({\hat{g}}_i^f\right)}^m+\underset{m=1}{\overset{\max m}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_{\mathrm{im}2}{\left({\hat{g}}_i^f\right)}^m\underset{r=1}{\overset{\max r}{\mathrm{\&Sigma;}}}{\mathrm{\&beta;}}_r{\mathrm{\&gamma;}}_i^{r-1}+\underset{m=1}{\overset{\max m}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_{\mathrm{im}3}{\left({\hat{g}}_i^f\right)}^m{\left(\underset{r=1}{\overset{\max r}{\mathrm{\&Sigma;}}}{\mathrm{\&beta;}}_r{\mathrm{\&gamma;}}_i^{r-1}\right)}^2---\left(15\right)$

$={\mathrm{\&Omega;}}_i+{\mathrm{\&Psi;}}_i\underset{r=1}{\overset{\max r}{\mathrm{\&Sigma;}}}{\mathrm{\&beta;}}_r{\mathrm{\&gamma;}}_i^{r-1}+{\mathrm{\&Gamma;}}_i{\left(\underset{r=1}{\overset{\max r}{\mathrm{\&Sigma;}}}{\mathrm{\&beta;}}_r{\mathrm{\&gamma;}}_i^{r-1}\right)}^2$

$={\mathrm{\&Omega;}}_i+\underset{r=1}{\overset{\max r}{\mathrm{\&Sigma;}}}{\mathrm{\&beta;}}_r\left({\mathrm{\&Psi;}}_i{\mathrm{\&gamma;}}_i^{r-1}\right)+\underset{r1=1}{\overset{\max r}{\mathrm{\&Sigma;}}}\underset{r2=1}{\overset{\max r}{\mathrm{\&Sigma;}}}{\mathrm{\&beta;}}_{r1}{\mathrm{\&beta;}}_{r2}\left({\mathrm{\&Gamma;}}_i{\mathrm{\&gamma;}}_i^{r1+r2-2}\right)$

与d＝0时的发散束CT数据一致性条件相结合有如下极小化问题：

${\mathrm{\&Phi;}}_1({\mathrm{\&beta;}}_r,m_{\mathrm{0,0}},{\mathrm{\&lambda;}}_0)=$

$\underset{\mathrm{\&beta;}}{\mathrm{\&Sigma;}}\{\underset{r1=1}{\overset{\max r}{\mathrm{\&Sigma;}}}\underset{r2=1}{\overset{\max r}{\mathrm{\&Sigma;}}}{\mathrm{\&beta;}}_{r1}{\mathrm{\&beta;}}_{r2}{\&Integral;}_{-\&infin;}^{+\&infin;}\mathrm{\&Gamma;}(t,\mathrm{\&beta;}-\mathrm{arctg}\frac{t}{D}){\mathrm{\&gamma;}}^{r1+r2-2}(t,\mathrm{\&beta;}-\mathrm{arctg}\frac{t}{D})\left(\frac{D^3}{{(D^2+t^2)}^{3/2}}\right)\mathrm{dt}$

$+\underset{r=1}{\overset{\max r}{\mathrm{\&Sigma;}}}{\mathrm{\&beta;}}_r{\&Integral;}_{-\&infin;}^{+\&infin;}\mathrm{\&Psi;}(t,\mathrm{\&beta;}-\mathrm{arctg}\frac{t}{D}){\mathrm{\&gamma;}}^{r-1}(t,\mathrm{\&beta;}-\mathrm{arctg}\frac{t}{D})\left(\frac{D^3}{{(D^2+t^2)}^{3/2}}\right)\mathrm{dt}---\left(16\right)$

$+{\&Integral;}_{-\&infin;}^{+\&infin;}\mathrm{\&Omega;}(t,\mathrm{\&beta;}-\mathrm{arctg}\frac{t}{D})\left(\frac{D^3}{{(D^2+t^2)}^{3/2}}\right)\mathrm{dt}-m_{\mathrm{0,0}}{\}}^2$

其中，m_0，0为不依赖于(t，β)的哑变量。

步骤3.2：H-Lλ₀校正方法

发散束宽能谱X线CT成像中，骨校正中λ₀可以通过解如下极小值问题确定：

${\mathrm{\&Phi;}}_2(m_{\mathrm{0,0}},{\mathrm{\&lambda;}}_0)=$

$\underset{\mathrm{\&beta;}}{\mathrm{\&Sigma;}}\{m_{\mathrm{0,0}}-m_{\mathrm{0,0}}^f+{\&Integral;}_{-\&infin;}^{+\&infin;}\{P_i^w[-\log \left({\&Integral;}_0^V{P}_i(V,E)e^{-{\mathrm{\&mu;}}_s{f}_i^s-{\mathrm{\&mu;}}_b{f}_i^b/{\mathrm{\&lambda;}}_0}\mathrm{dE}\right)]$

$-P_i^w\left[{\hat{g}}^f(t,\mathrm{\&beta;}-\mathrm{arctg}\frac{t}{D})\right]\}\left(\frac{D^3}{{(D^2+t^2)}^{3/2}}\right)\mathrm{dt}{\}}^2---\left(17\right)$

其中m_0，0为不依赖于(t，β)的哑变量，m_0，0 ^f可如下预先计算得到

$m_{\mathrm{0,0}}^f={\&Integral;}_{-\&infin;}^{+\&infin;}f(t,\mathrm{\&beta;}-\mathrm{arctg}\frac{t}{D})\left(\frac{D^3}{{(D^2+t^2)}^{3/2}}\right)\mathrm{dt}---\left(18\right)$

f(t，β)为水校正、预重建图像的重投影。

步骤3.3：H-L Bone校正方法：

发散束宽能谱X线CT成像中，骨校正中能谱函数可以通过解如下极小值确定：

${\mathrm{\&Phi;}}_3(m_{\mathrm{0,0}},c)=$

$\underset{\mathrm{\&beta;}}{\mathrm{\&Sigma;}}\{m_{\mathrm{0,0}}-m_{\mathrm{0,0}}^f+c({\&Integral;}_{-\&infin;}^{+\&infin;}{f}^s(t,\mathrm{\&beta;}-\mathrm{arctg}\frac{t}{D})\left(\frac{D^3}{{(D^2+t^2)}^{3/2}}\right)\mathrm{dt},{\&Integral;}_{-\&infin;}^{+\&infin;}{f}^b(t,\mathrm{\&beta;}-\mathrm{arctg}\frac{t}{D})\left(\frac{D^3}{{(D^2+t^2)}^{3/2}}\right)\mathrm{dt})$

$-{\&Integral;}_{-\&infin;}^{+\&infin;}{\hat{g}}^f(t,\mathrm{\&beta;}-\mathrm{arctg}\frac{t}{D})\left(\frac{D^3}{{(D^2+t^2)}^{3/2}}\right)\mathrm{dt}{\}}^2,---\left(19\right)$

其中m_0，0为不依赖于(t，β)的哑变量，m_0，0 ^f可如下预先计算得到

$m_{\mathrm{0,0}}^f={\&Integral;}_{-\&infin;}^{+\&infin;}f(t,\mathrm{\&beta;}-\mathrm{arctg}\frac{t}{D})\left(\frac{D^3}{{(D^2+t^2)}^{3/2}}\right)\mathrm{dt},---\left(20\right)$

c(f^s，f^b)＝c₁f^b+c₂f^s+c₃(f^b)²+c₄(f^s)²+c₅(f^bf^s)+c₆(f^b)^0.5。    (21)

步骤4：投影数据校正

根据已经确定出的中间参数(β_r，m_0，0，λ₀)、(m_0，0，λ₀)或(m_0，0，c)，计算出既满足发散束CT数据一致性又消除射束硬化效应影响的g_i ^f，用于CT图像重建中。

本发明的方法适合于各类X线CT设备的射束硬化校正功能的实现。射束硬化现象对于许多细小组织病变的诊断和工业CT等应用场合意义至关重要。由于射束硬化校正与成像对象有关，本发明方法与现有各类CT中，只采用模体校正的方法(如：水模校正)相比较，能够分别考虑不同密度组织的射束硬化效应的差异，自动适应成像对象的不同；与一定程度上考虑成像对象内容的校正算法(如：骨校正)相比较，可以自动选择尺度因子λ₀(如：General H-L与H-Lλ₀校正方法)，可以不需要能谱分布函数和尺度因子λ₀信息(如：H-L Bone校正方法)。整体而言，具有更好的校正性能。

附图说明

图1中上部图形是CT球管发出的X线光谱函数，下部图形是经过35mm铝衰减后的光谱函数(放大了5倍)；

图2是骨组织、肌肉和脂肪在120kVp电压球管发出的光谱中，衰减量和厚度的关系。粗实线代表肌肉组织，粗长虚线代表骨组织，粗短虚线代表脂肪；细实线代表肌肉组织衰减与线性关系的偏差，细长虚线代表骨组织衰减与线性关系的偏差，细短虚线代表脂肪衰减与线性关系的偏差，其中所有偏差曲线都已经放大了10倍。

图3是圆形扫描轨迹、扇形投影束的示意图，β为投影角度，S为射线源，g^f(t，β)为投影数据，O为旋转中心，D为S到O的距离；

图4～8中实验用的是FORBILD头部模体，用于说明射束硬化校正效果，其中亮的部分是骨皮组织，灰色部分是软组织，黑色部分是空气。重建方法采用滤波反投影(FBP)算法

图4是X线光子能量为66keV时的单能谱重建图像。

图5(a)是宽能谱重建图像。投影参数为：球管电压120kVp，附加一个铝材料的Bowtie补偿滤波器。图6～8中射束硬化校正前的投影参数与本图相同。

图5(b)是图4的单能谱重建图像与图5(a)的宽能谱重建图像之间的差图像，该差图像用于表示一个理想差图像应具备的校正性能与特点。

图6(a)是利用Geneal H-L校正方法获得的重建图像。

图6(b)是图4的单能谱重建图像与图6(a)的Geneal H-L校正重建图像之间的差图像。

图7(a)是利用H-Lλ₀校正方法获得的重建图像。

图7(b)是图4的单能谱重建图像与图7(a)的H-Lλ₀校正重建图像之间的差图像。

图8(a)是利用H-L Bone校正方法获得的重建图像。

图8(b)是图4的单能谱重建图像与图8(a)的H-L Bone校正重建图像之间的差图像。

以下结合附图对本发明作进一步的详细说明。

具体实施方式

1.圆形扫描轨迹、扇形束投影几何的CT数据一致性表达式

投影示意图见图(3)，如果

m_i，k＝∫∫_Cxⁱy^kf(x，y)dxdy  i≥0，k≥0，    (1)

$v_d\left(\mathrm{\&beta;}\right)={\&Integral;}_{-\&infin;}^{+\&infin;}{g}^f(t,\mathrm{\&beta;}-\mathrm{arctg}\frac{t}{D}){\left(\frac{\mathrm{tD}}{\sqrt{D^2+r^2}}\right)}^d\left(\frac{D^3}{{(D^2+t^2)}^{3/2}}\right)\mathrm{dt},---\left(2\right)$

$Q_d\left(\mathrm{\&beta;}\right)=\underset{r=1}{\overset{d}{\mathrm{\&Sigma;}}}\left(\begin{array}{c}d\\ r\end{array}\right)m_{r,d-r}{\cos }^r\left(\mathrm{\&beta;}\right){\sin }^{d-r}\left(\mathrm{\&beta;}\right),---\left(3\right)$

则有CT数据一致性表达式

v_d(β)＝Q_d(β)，    (4)

其中，0≤d≤N-1，而N为投影角度的个数。该CT数据一致性表达式将始终用于构造以下射束硬化校正方法对应得极小化求解问题。

2.General H-L投影射束硬化校正方法

2.1：训练数据采集与计算

由于需要对一些训练数据进行采集和计算，以便得到单能谱投影与宽能谱投影之间的约束关系。在医学诊断X线谱段，X线光子与类似于人体组织的物质的相互作用主要表现为光电吸收和康普顿散射两种效应。在采集训练数据过程中，采用X线衰减特性与人体软组织近似的有机玻璃和与骨组织近似的纯铝作为模体。则有：

$Z_{\mathrm{eff},i}^3=[(1-{\mathrm{\&gamma;}}_i)Z_s^3+{\mathrm{\&gamma;}}_i{Z}_b^3]---\left(5\right)$

其中，Z_s、Z_b与Z_eff，i分别为有机玻璃原子序数、纯铝原子序数与第i条X线经过所有材料的等效原子序数。

这样用于训练的单能谱投影与宽能谱投影之间的约束关系可以选为：

$g_i=\underset{m=1}{\overset{\max m}{\mathrm{\&Sigma;}}}\underset{n=1}{\overset{\max n}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_{\mathrm{imn}}{\left(Z_{\mathrm{eff},i}^3\right)}^{n-1}{\left({\hat{g}}_i\right)}^m,---\left(6\right)$

其中g_i与Z_eff，i计算得到，而

从采集的训练数据I_i计算得到，从而可以确定出系数α_imn，将之存储下来用于后续数据处理。对于训练数据的采集，不强调几何模式。

2.2：γ与Z³之间的映射函数逼近

X线成像现场数据的Z_eff，i无法精确确定，但是可以用γ_i对其进行逼近

$Z_{\mathrm{eff},i}^3=\underset{r=1}{\overset{\max r}{\mathrm{\&Sigma;}}}{\mathrm{\&beta;}}_r{\mathrm{\&gamma;}}_i^{r-1},---\left(7\right)$

则有

$g_i^f=\underset{m=1}{\overset{\max m}{\mathrm{\&Sigma;}}}\underset{n=1}{\overset{\max n}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_{\mathrm{imn}}{\left(\underset{r=1}{\overset{\max r}{\mathrm{\&Sigma;}}}{\mathrm{\&beta;}}_r{\mathrm{\&gamma;}}_i^{r-1}\right)}^{n-1}{\left({\hat{g}}_i^f\right)}^m,---\left(8\right)$

其中α_imn在训练阶段已确定，γ_i通过水模校正、预分割、重投影、面密度比例计算等一系列步骤得到，由X线成像现场数据计算得到，g_i ^f虽未知但其必满足发散束CT数据一致性条件。通过所有这些已知条件，可以反解出β_r，实现γ与Z³之间的映射函数逼近。

反解β_r过程对应于求解如下极小化问题

${\mathrm{\&Phi;}}_1({\mathrm{\&beta;}}_r,m_{\mathrm{0,0}},{\mathrm{\&lambda;}}_0)=$

$\underset{\mathrm{\&beta;}}{\mathrm{\&Sigma;}}\{\underset{r1=1}{\overset{\max r}{\mathrm{\&Sigma;}}}\underset{r2=1}{\overset{\max r}{\mathrm{\&Sigma;}}}{\mathrm{\&beta;}}_{r1}{\mathrm{\&beta;}}_{r2}{\&Integral;}_{-\&infin;}^{+\&infin;}\mathrm{\&Gamma;}(t,\mathrm{\&beta;}-\mathrm{arctg}\frac{t}{D}){\mathrm{\&gamma;}}^{r1+r2-2}(t,\mathrm{\&beta;}-\mathrm{arctg}\frac{t}{D})\left(\frac{D^3}{{(D^2+t^2)}^{3/2}}\right)\mathrm{dt}$

$+\underset{r=1}{\overset{\max r}{\mathrm{\&Sigma;}}}{\mathrm{\&beta;}}_r{\&Integral;}_{-\&infin;}^{+\&infin;}\mathrm{\&Psi;}(t,\mathrm{\&beta;}-\mathrm{arctg}\frac{t}{D}){\mathrm{\&gamma;}}^{r-1}(t,\mathrm{\&beta;}-\mathrm{arctg}\frac{t}{D})\left(\frac{D^3}{{(D^2+t^2)}^{3/2}}\right)\mathrm{dt},---\left(9\right)$

$+{\&Integral;}_{-\&infin;}^{+\&infin;}\mathrm{\&Omega;}(t,\mathrm{\&beta;}-\mathrm{arctg}\frac{t}{D})\left(\frac{D^3}{{(D^2+t^2)}^{3/2}}\right)\mathrm{dt}-m_{\mathrm{0,0}}{\}}^2$

其中m_0，0为不依赖于(t，β)的哑变量。

图6(a)是利用申请人提出的General H-L投影射束硬化校正方法校正后的重建图像，采用的重建方法是FBP，并采用R-L型RAMP滤波函数。为了说明校正效果，在图6(b)中显示图4单能谱重建图像与图6(a)Geneal H-L校正重建图像之间的差图像。比较图5(b)与6(b)可以看出，杯形伪影与条纹伪影都得到良好的校正，这说明该校正方法对硬膜外伪影这一CT设备难以克服的问题，将会有较好的解决方案途径。

3.H-Lλ₀投影射束硬化校正方法

发散束宽能谱X线CT成像中，骨校正中最优λ₀可以通过解如下极小值问题确定

${\mathrm{\&Phi;}}_2(m_{\mathrm{0,0}},{\mathrm{\&lambda;}}_0)=$

$\underset{\mathrm{\&beta;}}{\mathrm{\&Sigma;}}\{m_{\mathrm{0,0}}-m_{\mathrm{0,0}}^f+{\&Integral;}_{-\&infin;}^{+\&infin;}\{P_i^w[-\log \left({\&Integral;}_0^V{P}_i(V,E)e^{-{\mathrm{\&mu;}}_s{f}_i^s-{\mathrm{\&mu;}}_b{f}_i^b/{\mathrm{\&lambda;}}_0}\mathrm{dE}\right)]$

$-P_i^w\left[{\hat{g}}^f(t,\mathrm{\&beta;}-\mathrm{arctg}\frac{t}{D})\right]\}\left(\frac{D^3}{{(D^2+t^2)}^{3/2}}\right)\mathrm{dt}{\}}^2,---\left(10\right)$

其中m_0，0为不依赖于(t，β)的哑变量，m_0，0 ^f可如下预先计算得到

$m_{\mathrm{0,0}}^f={\&Integral;}_{-\&infin;}^{+\&infin;}f(t,\mathrm{\&beta;}-\mathrm{arctg}\frac{t}{D})\left(\frac{D^3}{{(D^2+t^2)}^{3/2}}\right)\mathrm{dt},---\left(11\right)$

f(t，β)为水校正、预重建图像的重投影。

当被成像物体为旋转圆对称或者近似旋转圆对称时，极小化问题Φ₂无法同时确定出m_0，0与λ₀。为此进一步将Φ₂转变为如下一对极小化问题Φ₂₁与Φ₂₂，极小化过程采用非线性最小二乘迭代交替进行

${\mathrm{\&Phi;}}_{21}^{\left(k\right)}\left({\mathrm{\&lambda;}}_0^{\left(k\right)}\right)=\underset{\mathrm{\&beta;}}{\mathrm{\&Sigma;}}\{m_{\mathrm{0,0}}^{(k-1)}-m_{\mathrm{0,0}}^f+{\&Integral;}_{-\&infin;}^{+\&infin;}\{P_i^w[-\log \left({\&Integral;}_0^V{P}_i(V,E)e^{-{\mathrm{\&mu;}}_s{f}_i^s-{\mathrm{\&mu;}}_b{f}_i^b/{\mathrm{\&lambda;}}_0^{\left(k\right)}}\mathrm{dE}\right)]$

$-P_i^w\left[{\hat{g}}^f(t,\mathrm{\&beta;}-\mathrm{arctg}\frac{t}{D})\right]\}\left(\frac{D^3}{{(D^2+t^2)}^{3/2}}\right)\mathrm{dt}{\}}^2,---\left(12\right)$

${\mathrm{\&Phi;}}_{22}^{\left(k\right)}\left(m_{\mathrm{0,0}}^{\left(k\right)}\right)=\underset{\mathrm{\&beta;}}{\mathrm{\&Sigma;}}\{m_{\mathrm{0,0}}^{\left(k\right)}-m_{\mathrm{0,0}}^f+{\&Integral;}_{-\&infin;}^{+\&infin;}\{P_i^w[-\log \left({\&Integral;}_0^V{P}_i(V,E)e^{-{\mathrm{\&mu;}}_s{f}_i^s-{\mathrm{\&mu;}}_b{f}_i^b/{\mathrm{\&lambda;}}_0^{\left(k\right)}}\mathrm{dE}\right)]$

$-P_i^w\left[{\hat{g}}^f(t,\mathrm{\&beta;}-\mathrm{arctg}\frac{t}{D})\right]\}\left(\frac{D^3}{{(D^2+t^2)}^{3/2}}\right)\mathrm{dt}{\}}^2,---\left(13\right)$

k＝1，2，3，L，max k，取 $m_{\mathrm{0,0}}^{\left(0\right)}=m_{\mathrm{0,0}}^f.$

最优λ₀确定后，就可以如同原始骨校正完成射束硬化校正和CT重建。

图7(a)是利用申请人提出的H-Lλ₀投影射束硬化校正方法校正后的重建图像，采用的重建方法是FBP，并采用R-L型RAMP滤波函数。为了说明校正效果，在图7(b)中显示图4单能谱重建图像与图7(a)H-Lλ₀校正重建图像之间的差图像。比较图5(b)与7(b)可以看出，杯形伪影与条纹伪影也都得到良好的校正。

4.H-L Bone投影射束硬化校正方法

发散束宽能谱X线CT成像中，如果能谱分布函数没有明确信息，骨校正中能谱函数可以通过解如下极小值问题逼近确定

${\mathrm{\&Phi;}}_3(m_{\mathrm{0,0}},c)=$

$\underset{\mathrm{\&beta;}}{\mathrm{\&Sigma;}}\{m_{\mathrm{0,0}}-m_{\mathrm{0,0}}^f+c({\&Integral;}_{-\&infin;}^{+\&infin;}{f}^s(t,\mathrm{\&beta;}-\mathrm{arctg}\frac{t}{D})\left(\frac{D^3}{{(D^2+t^2)}^{3/2}}\right)\mathrm{dt},{\&Integral;}_{-\&infin;}^{+\&infin;}{f}^b(t,\mathrm{\&beta;}-\mathrm{arctg}\frac{t}{D})\left(\frac{D^3}{{(D^2+t^2)}^{3/2}}\right)\mathrm{dt})$

$-{\&Integral;}_{-\&infin;}^{+\&infin;}{\hat{g}}^f(t,\mathrm{\&beta;}-\mathrm{arctg}\frac{t}{D})\left(\frac{D^3}{{(D^2+t^2)}^{3/2}}\right)\mathrm{dt}{\}}^2,---\left(14\right)$

其中m_0，0为不依赖于(t，β)的哑变量，m_0，0 ^f可如下预先计算得到

$m_{\mathrm{0,0}}^f={\&Integral;}_{-\&infin;}^{+\&infin;}f(t,\mathrm{\&beta;}-\mathrm{arctg}\frac{t}{D})\left(\frac{D^3}{{(D^2+t^2)}^{3/2}}\right)\mathrm{dt},---\left(15\right)$

c(f^s，f^b)＝c₁f^b+c₂f^s+c₃(f^b)²+c₄(f^s)²+c₅(f^bf^s)+c₆(f^b)^0.5。    (16)

图8(a)是利用申请人提出的H-L Bone投影射束硬化校正方法校正后的重建图像，采用的重建方法是FBP，并采用R-L型RAMP滤波函数。为了说明校正效果，在图8(b)中显示图4单能谱重建图像与图8(a)H-L Bone校正重建图像之间的差图像。比较图5(b)与8(b)可以看出，杯形伪影与条纹伪影也都得到良好的校正。

Claims (1)

1.一种基于CT数据一致性的投影数据射束硬化校正方法，该方法首先对图像进行水模校正，然后对校正后图像进行预分割、重投影与面密度比例计算，最后根据CT数据一致性对不同角度的投影数据进行射束硬化效应校正，其特征在于，所述的射束硬化效应校正具体包括以下步骤：

步骤1，H-L射束硬化校正方法

步骤1.1，General H-L校正方法

第一步，训练数据采集与计算：

对一些训练数据进行采集和计算，以便得到单能谱投影与宽能谱投影之间的约束关系，在医学诊断X线谱段，由于X线光子与类似于人体组织的物质的相互作用主要表现为光电吸收和康普顿散射两种效应，所以宽能谱X线成像的投影数据可以进一步近似表示为：

$I_i\&ap;I_0{\&Integral;}_0^V{P}_i(V,E)\exp \{-{\&Integral;}_{L_i}[\mathrm{bd}\left(l\right)\frac{12.4}{E^3}+\mathrm{bp}\left(l\right)f_{\mathrm{KN}}\left(E\right)\left]\mathrm{\&rho;}\left(l\right)\mathrm{dl}\right\}\mathrm{dE}$

在采集训练数据过程中，仅考虑两种材料，它们的X线衰减特性分别与人体软组织和骨组织近似；则有：

$I_i$

$\&ap;I_0{\&Integral;}_0^V{P}_i(V,E)\exp \{-{\mathrm{\&rho;}}_i^s{L}_i^s[{\mathrm{bd}}_i^s\frac{12.4}{E^3}+{\mathrm{bp}}_i^s{f}_{\mathrm{KN}}\left(E\right)]-{\mathrm{\&rho;}}_i^b{L}_i^b[{\mathrm{bd}}_i^b\frac{12.4}{E^3}+{\mathrm{bp}}_i^b{f}_{\mathrm{KN}}\left(E\right)\left]\right\}\mathrm{dE}$

$=I_0{\&Integral;}_0^V{P}_i(V,E)\exp \{-{\mathrm{\&rho;}}_i^a(1-{\mathrm{\&gamma;}}_i)[{\mathrm{bd}}_i^s\frac{12.4}{E^3}+{\mathrm{bp}}_i^s{f}_{\mathrm{KN}}\left(E\right)]-{\mathrm{\&rho;}}_i^a{\mathrm{\&gamma;}}_i[{\mathrm{bd}}_i^b\frac{12.4}{E^3}+{\mathrm{bp}}_i^b{f}_{\mathrm{KN}}\left(E\right)\left]\right\}\mathrm{dE}$

$=I_0{\&Integral;}_0^V{P}_i(V,E)\exp \{-{\mathrm{\&rho;}}_i^a\{[(1-{\mathrm{\&gamma;}}_i){\mathrm{bd}}_i^s+{\mathrm{\&gamma;}}_i{\mathrm{bd}}_i^b]\frac{12.4}{E^3}+[(1-{\mathrm{\&gamma;}}_i){\mathrm{bp}}_i^s+{\mathrm{\&gamma;}}_{i}b{p}_i^b]f_{\mathrm{KN}}\left(E\right)\left\}\right\}\mathrm{dE}$

$=I_0{\&Integral;}_0^V{P}_i(V,E)\exp \{-{\mathrm{\&rho;}}_i^a[{\mathrm{bd}}_i\frac{12.4}{E^3}+{\mathrm{bp}}_i{f}_{\mathrm{KN}}\left(E\right)\left]\right\}\mathrm{dE}$

其中：

${\mathrm{\&rho;}}_i^a={\mathrm{\&rho;}}_i^s{L}_i^s+{\mathrm{\&rho;}}_i^b{L}_i^b$

${\mathrm{\&gamma;}}_i={\mathrm{\&rho;}}_i^b{L}_i^b/({\mathrm{\&rho;}}_i^s{L}_i^s+{\mathrm{\&rho;}}_i^b{L}_i^b)$

这样

$\frac{{\mathrm{bp}}_i}{b{d}_i}$

$=\frac{(1-{\mathrm{\&gamma;}}_i)b{p}_i^s+{\mathrm{\&gamma;}}_{i}b{p}_i^b}{(1-{\mathrm{\&gamma;}}_i)b{d}_i^s+{\mathrm{\&gamma;}}_{i}b{d}_i^b}$

$\&ap;\frac{(1-{\mathrm{\&gamma;}}_i)k{Z}_s^{3}b{d}_i^s+{\mathrm{\&gamma;}}_{i}k{Z}_b^{3}b{d}_i^b}{(1-{\mathrm{\&gamma;}}_i)b{d}_i^s+{\mathrm{\&gamma;}}_{i}b{d}_i^b}$

$\&ap;k[(1-{\mathrm{\&gamma;}}_i)Z_s^3+{\mathrm{\&gamma;}}_i{Z}_b^3]$

所以有如下比例关系：

$\frac{{\mathrm{bp}}_i}{{\mathrm{bd}}_i}\&Proportional;Z_{\mathrm{eff},i}^3$

其中，Z_s、Z_b与Z_eff，i分别为软组织近似材料原子序数、骨组织近似材料原子序数与第i条X线经过所有材料的等效原子序数；用于训练的单能谱投影与宽能谱投影之间的约束关系为：

$g_i=\underset{m=1}{\overset{\max m}{\mathrm{\&Sigma;}}}\underset{n=1}{\overset{\max n}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_{\mathrm{imn}}{\left(Z_{\mathrm{eff},i}^3\right)}^{n-1}{\left({\hat{g}}_i\right)}^m$

其中g_i与Z_eff，i计算得到，而

从采集的训练数据I_i计算得到或者直接计算得到，从而确定出系数α_imn，将之存储下来用于后续数据处理；对于训练数掘的采集与计算不强调几何模式；该公式也可以采用适合于计算的其他逼近多项式形式；

第二步，γ与Z³之间的映射函数逼近：

对于X线成像现场数据，对Z_eff，i采用下式进行逼近表达

$Z_{\mathrm{eff},i}^3=\underset{r=1}{\overset{\max r}{\mathrm{\&Sigma;}}}{\mathrm{\&beta;}}_r{\mathrm{\&gamma;}}_i^{r-1}$

则有：

$g_i^f=\underset{m=1}{\overset{\max m}{\mathrm{\&Sigma;}}}\underset{n=1}{\overset{\max n}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_{\mathrm{imn}}{\left(\underset{r=1}{\overset{\max r}{\mathrm{\&Sigma;}}}{\mathrm{\&beta;}}_r{\mathrm{\&gamma;}}_i^{r-1}\right)}^{n-1}{\left({\hat{g}}_i^f\right)}^m$

其中，α_imn在训练阶段已确定，γ_i通过水模校正、预分割、重投影、面密度比例计算一系列步骤得到，

由X线成像现场数据计算得到，g_i ^f必满足发散束CT数据一致性条件，该公式也可以采用适合于计算的其他逼近多项式形式；

通过所有这些已知条件，反解出β_r，实现γ与Z³之间的映射函数逼近；反解β_r的方法如下：

取maxn＝3，有

$g_i^f=\underset{m=1}{\overset{\max m}{\mathrm{\&Sigma;}}}\underset{n=1}{\overset{3}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_{\mathrm{imn}}{\left(\underset{r=1}{\overset{\max r}{\mathrm{\&Sigma;}}}{\mathrm{\&beta;}}_r{\mathrm{\&gamma;}}_i^{r-1}\right)}^{n-1}{\left({\hat{g}}_i^f\right)}^m$

$=\underset{m=1}{\overset{\max m}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_{\mathrm{im}1}{\left({\hat{g}}_i^f\right)}^m+\underset{m=1}{\overset{\max m}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_{\mathrm{im}2}{\left({\hat{g}}_i^f\right)}^m\underset{r=1}{\overset{\max r}{\mathrm{\&Sigma;}}}{\mathrm{\&beta;}}_r{\mathrm{\&gamma;}}_i^{r-1}+\underset{m=1}{\overset{\max m}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_{\mathrm{im}3}{\left({\hat{g}}_i^f\right)}^m{\left(\underset{r=1}{\overset{\max r}{\mathrm{\&Sigma;}}}{\mathrm{\&beta;}}_r{\mathrm{\&gamma;}}_i^{r-1}\right)}^2$

$={\mathrm{\&Omega;}}_i+{\mathrm{\&Psi;}}_i\underset{r=1}{\overset{\max r}{\mathrm{\&Sigma;}}}{\mathrm{\&beta;}}_r{\mathrm{\&gamma;}}_i^{r-1}+{\mathrm{\&Gamma;}}_i{\left(\underset{r=1}{\overset{\max r}{\mathrm{\&Sigma;}}}{\mathrm{\&beta;}}_r{\mathrm{\&gamma;}}_i^{r-1}\right)}^2$

$={\mathrm{\&Omega;}}_i+\underset{r=1}{\overset{\max r}{\mathrm{\&Sigma;}}}{\mathrm{\&beta;}}_r\left({\mathrm{\&Psi;}}_i{\mathrm{\&gamma;}}_i^{r-1}\right)+\underset{r1=1}{\overset{\max r}{\mathrm{\&Sigma;}}}\underset{r2=1}{\overset{\max r}{\mathrm{\&Sigma;}}}{\mathrm{\&beta;}}_{r1}{\mathrm{\&beta;}}_{r2}\left({\mathrm{\&Gamma;}}_i{\mathrm{\&gamma;}}_i^{r1+r2-2}\right)$

与d＝0时的发散束CT数据一致性条件相结合，有如下极小化问题：

${\mathrm{\&Phi;}}_1({\mathrm{\&beta;}}_r,m_{\mathrm{0,0}},{\mathrm{\&lambda;}}_0)=$

$\underset{\mathrm{\&beta;}}{\mathrm{\&Sigma;}}\{\underset{r1=1}{\overset{\max r}{\mathrm{\&Sigma;}}}\underset{r2=1}{\overset{\max r}{\mathrm{\&Sigma;}}}{\mathrm{\&beta;}}_{r1}{\mathrm{\&beta;}}_{r2}{\&Integral;}_{-\&infin;}^{+\&infin;}\mathrm{\&Gamma;}(t,\mathrm{\&beta;}-\mathrm{arctg}\frac{t}{D}){\mathrm{\&gamma;}}^{r1+r2-2}(t,\mathrm{\&beta;}-\mathrm{arctg}\frac{t}{D})\left(\frac{D^3}{{(D^2+t^2)}^{3/2}}\right)\mathrm{dt}$

$+\underset{r=1}{\overset{\max r}{\mathrm{\&Sigma;}}}{\mathrm{\&beta;}}_r{\&Integral;}_{-\&infin;}^{+\&infin;}\mathrm{\&Psi;}(t,\mathrm{\&beta;}-\mathrm{arctg}\frac{t}{D}){\mathrm{\&gamma;}}^{r-1}(t,\mathrm{\&beta;}-\mathrm{arctg}\frac{t}{D})\left(\frac{D^3}{{(D^2+t^2)}^{3/2}}\right)\mathrm{dt}$

$+{\&Integral;}_{-\&infin;}^{+\&infin;}\mathrm{\&Omega;}(t,\mathrm{\&beta;}-\mathrm{arctg}\frac{t}{D})\left(\frac{D^3}{{(D^2+t^2)}^{3/2}}\right)\mathrm{dt}-m_{\mathrm{0,0}}{\}}^2$

其中，m_0，0为不依赖于(t，β)的哑变量；max n也可以取其他数值；

步骤1.2，H-Lλ₀校正方法：

发散束宽能谱X线CT成像中，骨校正中参数λ₀，通过解极小值问题确定：

${\mathrm{\&Phi;}}_2(m_{\mathrm{0,0}},{\mathrm{\&lambda;}}_0)=$

$\underset{\mathrm{\&beta;}}{\mathrm{\&Sigma;}}\{m_{\mathrm{0,0}}-m_{\mathrm{0,0}}^f+{\&Integral;}_{-\&infin;}^{+\&infin;}\{P_i^w[-\log \left({\&Integral;}_0^V{P}_i(V,E)e^{-{\mathrm{\&mu;}}_s{f}_i^s-{\mathrm{\&mu;}}_b{f}_i^b/{\mathrm{\&lambda;}}_0}\mathrm{dE}\right)]$

$-P_i^w\left[{\hat{g}}^f(t,\mathrm{\&beta;}-\mathrm{arctg}\frac{t}{D})\right]\}\left(\frac{D^3}{{(D^2+t^2)}^{3/2}}\right)\mathrm{dt}{\}}^2$

其中，m_0，0为不依赖于(t，β)的哑变量，m_0，0 ^f由下式预先计算得到

$m_{\mathrm{0,0}}^f={\&Integral;}_{-\&infin;}^{+\&infin;}f(t,\mathrm{\&beta;}-\mathrm{arctg}\frac{t}{D})\left(\frac{D^3}{{(D^2+t^2)}^{3/2}}\right)\mathrm{dt}$

f(t，β)为水校正、预重建图像的重投影；

步骤1.3，H-L Bone校正方法：

发散束宽能谱X线CT成像中，骨校正中能谱函数通过解如下极小值确定：

${\mathrm{\&Phi;}}_3(m_{\mathrm{0,0}},c)=$

$\underset{\mathrm{\&beta;}}{\mathrm{\&Sigma;}}\{m_{\mathrm{0,0}}-m_{\mathrm{0,0}}^f+c({\&Integral;}_{-\&infin;}^{+\&infin;}{f}^s(t,\mathrm{\&beta;}-\mathrm{arctg}\frac{t}{D})\left(\frac{D^3}{{(D^2+t^2)}^{3/2}}\right)\mathrm{dt},{\&Integral;}_{-\&infin;}^{+\&infin;}{f}^b(t,\mathrm{\&beta;}-\mathrm{arctg}\frac{t}{D})\left(\frac{D^3}{{(D^2+t^2)}^{3/2}}\right)\mathrm{dt})$

$-{\&Integral;}_{-\&infin;}^{+\&infin;}{\hat{g}}^f(t,\mathrm{\&beta;}-\mathrm{arctg}\frac{t}{D})\left(\frac{D^3}{{(D^2+t^2)}^{3/2}}\right)\mathrm{dt}{\}}^2$

其中，m_0，0为不依赖于(t，β)的哑变量，m_0，0 ^f按下式预先计算得到：

$m_{\mathrm{0,0}}^f={\&Integral;}_{-\&infin;}^{+\&infin;}f(t,\mathrm{\&beta;}-\mathrm{arctg}\frac{t}{D})\left(\frac{D^2}{{(D^2+t^2)}^{3/2}}\right)\mathrm{dt}$

式中，c(f^s，f^b)为任意容易求解的逼近多项式；

步骤2，投影数据校正：

根据确定出的中间参数(β_r，m_0，0，λ₀)、(m_0，0，λ₀)或(m_0，0，c)，计算出既满足发散束CT数据一致性又消除射束硬化效应影响的g_i ^f，用于CT图像重建中。

Images