CN103530849A - 一种优化的K-edge成像方法 - Google Patents

Abstract

本发明公开了一种优化的K-edge成像方法，本方法在K-edge前后两个X射线能量段进行成像以提高已知材料成像对比度。K-edge前后两个X射线能量段的宽度决定了两个能谱CT图像的对比度和噪声水平，针对如何设置K-edge前后两个X射线能量段的宽度问题，本发明引入信号差异噪声比（Signal Difference to Noise Ratio,SDNR）作为最优化准则，在这个最优化准则的约束下，选取了最优的X射线能量段的宽度，然后进行K-edge成像。这种成像方法，能够在保证两个重建能谱CT图像感兴趣区域噪声最小化的同时，可以获得最大的对比度差异，从而达到提高已知材料成像对比度的目的。可以用于生物医学CT成像领域。

Description

一种优化的K-edge成像方法

技术领域

本发明属于X射线能谱CT成像领域，具体涉及一种优化的K-edge成像方法。

背景技术

K-edge是一个物理现象，也叫K吸收边缘。随着X射线能量的增加，X射线对物质的衰减系数逐渐减小，但某些物质对特定能量下的X射线光子吸收特别大，导致X射线光子衰减系数的突然增加，这种随着能量的增加，X射线衰减系数陡然增加的现象，称为K-edge。从量子物理学角度讲，物质原子内部离原子核最近电子层（K电子层）中的自由电子，容易与特定能量下X射线光子发生相互作用（光电吸收），导致这种物质对该能量下X射线光子吸收性特别大。不同元素其原子结构不同，因此相应的K-edge特性也不同。K-edge特性表现为在某个特定能量下，X射线衰减系数发生跳变。

在物质K-edge的前后，X射线衰减系数差异很大。因此，可以利用衰减系数跳变这一物理特性，选择不同X射线能量范围进行成像以提高已知材料成像对比度。这种成像方式，称为K-edge成像技术。随着X射线能谱CT技术的不断发展，K-edge成像技术已是X射线能谱CT的一个非常重要应用。传统或常规X-CT图像中，不同软组织的对比度较低，往往需要借助一些对比剂（造影剂）来提高图像中软组织对比度。纵然，在目前K-edge成像技术研究中，人们利用X射线能谱CT研究生物医学样本时，虽也需要借助一些对比剂，但已不再是简单的借助对比剂以提高软组织的对比度，而是充分利用了对比剂的K-edge特性进行K-edge成像，以便提供更多的认知信息。

发明内容

对于多色（多能量）X射线源，有限宽度的X射线能量段所包含的光子数是有限的，K-edge前后两个X射线能量段的宽度，决定了两个能谱CT图像的对比度和噪声水平。因此本发明需要解决的问题是如何设置K-edge前后两个X射线能量段的宽度，以保证两个重建能谱CT图像感兴趣区域噪声最小化的同时，获得最大的对比度差异，从而提高已知材料成像对比度。鉴于此，本发明的目的是提供一种优化的K-edge成像方法。

从核物理学角度分析，射线源产生的X射线光子数，是服从泊松随机分布的。假设射线源产生的X射线光子数为n，于是，X射线光子数n的概率密度函数，则可表示为

$P\left\{n\right\}=\frac{e^{-I_0}}{n!}{I}_0^n---\left(1\right)$

其中，I₀为X射线光子数n的泊松分布均值。

假设X射线源产生的X射线光子数为n，穿过检测物体后探测到的X射线光子数为m，未被物体吸收和散射而被探测器接受到的X射线光子数m概率，则可表示为

$P\left\{m\right|n\}=\left(\begin{array}{c}n\\ m\end{array}\right)P^m{(1-\mathrm{\ρ})}^{n-m},m=0,...,n---\left(2\right)$

其中,p=e^-g,

是X射线经过物体后的整体衰减系数，μ(l)表示衰减系数曲线，X射线光子数m概率函数，则可表示为

$P\left\{m\right\}=\underset{n=m}{\overset{\∞}{\mathrm{\Σ}}}P\left\{m\right|n\}=\underset{n=m}{\overset{\∞}{\mathrm{\Σ}}}\frac{e^{-I_0}}{n!}{I}_0^n\left(\begin{array}{c}n\\ m\end{array}\right)P^m{(1-p)}^{n-m}=\frac{e^{-I_{0}p}}{m!}{\left(I_{0}p\right)}^m---\left(3\right)$

根据公式（3），可知X射线光子数m，仍然服从泊松分布，其泊松分布均值，则可表示为

E(m)=I₀e^-g      (4)

泊松分布方差为

Var(m)=I₀e^-g          (5)

当X射线束经过物体后，衰减后的X射线强度I，则可表示为

I=I₀e^-g+N_I         (6)其中，N_I为随机噪声，根据公式（6），计算出随机噪声N_I的均值和方差为

E(N_I)=0         (7)

Var(N_I)=I₀e^-g        (8)

在CT图像重建中，一个角度θ下的投影正弦图S_θ，通过测探到的X射线强度I₀和I，可以计算出S_θ为

$S_{\mathrm{\θ}}=-\ln \frac{I}{I_0}---\left(9\right)$

将公式（6）代入公式（9），可得

$\begin{array}{c}{S}_{\mathrm{\θ}}=-\ln \frac{I_0{e}^{-g}+N_1}{I_0}\\ =-\ln e^{-g}(1+\frac{N_1}{I_0{e}^{-g}})\\ =g-\ln (1+\frac{N_1}{I_0{e}^{-g}})\\ \≈g-\frac{N_I}{I_0{e}^{-g}}\end{array}---\left(10\right)$

根据公式（10），计算出投影正弦图S_θ的均值和方差为

E(S_θ)=g       (11)

$\mathrm{Var}\left(S_{\mathrm{\θ}}\right)=\frac{1}{I_0{e}^{-g}}---\left(12\right)$

与此同时，根据公式（10），可以得到投影正弦图S_θ的噪声分量N_S为

$N_S=\frac{N_I}{I_0{e}^{-g}}---\left(13\right)$

利用经典的解析重建算法——滤波反投影（Filtered Backprojection,FBP），可以得到重建图像f的数学表达式为

$\begin{array}{c}f=\underset{0}{\overset{\mathrm{\π}}{\∫}}\mathrm{d\θ}\underset{-\∞}{\overset{\∞}{\∫}}{S}_{\mathrm{\θ}}\left(t^{\′}\right)h(t-t^{\′})\mathrm{dt}\\ =\underset{0}{\overset{\mathrm{\π}}{\∫}}\mathrm{d\θ}\underset{-\∞}{\overset{\∞}{\∫}}(g\left(t^{\′}\right)+N_s\left(t^{\′}\right))h(t-t^{\′})d{t}^{\′}\end{array}---\left(14\right)$

其中，t=xcosθ+ysinθ，根据公式（14），计算出重建图像f的均值和方差为

$\begin{array}{c}E\left(f\right)=E\left(\underset{0}{\overset{\mathrm{\π}}{\∫}}\mathrm{d\θ}\underset{-\∞}{\overset{\∞}{\∫}}{S}_{\mathrm{\θ}}\left(t^{\′}\right)h(t-t^{\′})\mathrm{dt}\right)\\ =\underset{0}{\overset{\mathrm{\π}}{\∫}}\mathrm{d\θ}\underset{-\∞}{\overset{\∞}{\∫}}E\left(S_{\mathrm{\θ}}\left(t^{\′}\right)\right)h(t-t^{\′})\mathrm{dt}\\ =\underset{0}{\overset{\mathrm{\π}}{\∫}}\mathrm{d\θ}\underset{-\∞}{\overset{\∞}{\∫}}g\left(t^{\′}\right)h(t-t^{\′})d{t}^{\′}\end{array}---\left(15\right)$

$\begin{array}{c}\mathrm{Var}\left(f\right)=\mathrm{Var}\left(\underset{0}{\overset{\mathrm{\π}}{\∫}}\mathrm{d\θ}\underset{-\∞}{\overset{\∞}{\∫}}{S}_{\mathrm{\θ}}\left(t^{\′}\right)h(t-t^{\′})\mathrm{dt}\right)\\ =\underset{0}{\overset{\mathrm{\π}}{\∫}}\mathrm{d\θ}\underset{-\∞}{\overset{\∞}{\∫}}\mathrm{Var}\left(S_{\mathrm{\θ}}\left(t^{\′}\right)\right)h^2(t-t^{\′})\mathrm{dt}\\ =\underset{0}{\overset{\mathrm{\π}}{\∫}}\mathrm{d\θ}\underset{-\∞}{\overset{\∞}{\∫}}\frac{1}{I_0{e}^{-g\left(t^{\′}\right)}}{h}^2(t-t^{\′})d{t}^{\′}\end{array}---\left(16\right)$

由此，重建图像f，可以分解成期望图像E(f)和图像噪声N_f，即

f=E(f)+N_f(17)

其中，图像噪声N_f的均值和方差可表示为

E(N_f)=0(18)

$\mathrm{Var}\left(N_f\right)=\underset{0}{\overset{\mathrm{\π}}{\∫}}\mathrm{d\θ}\underset{-\∞}{\overset{\∞}{\∫}}\frac{1}{I_0{e}^{-g\left(t^{\′}\right)}}{h}^2(t-t^{\′})d{t}^{\′}---\left(19\right)$

上述对CT成像的基本理论进行了分析，这为本发明基于X射线能谱CT的K-edge成像技术而提出的优化的K-edge成像方法提供了理论基础。

有鉴于此，为实现上述目的，本发明的方法引入信号差异噪声比（SDNR）作为寻找K-edge前后两个最佳X射线能量段宽度的准则，在这个准则的约束下，选取了最佳的X射线能量段的宽度w，以此进行K-edge成像。

假设多色（多能量）X射线源产生的不同能量的X射线光子，它是服从函数C(E)分布的，函数C(E)的分布曲线。在K-edge两侧，取两个X射线能量段进行成像。这里，假设两个X射线能量段宽度相等，均为w，于是，两个X射线能量段光子数，则分别为

$I_R=\underset{K^{+}}{\overset{K+w}{\∫}}C\left(E\right)\mathrm{dE}---\left(20\right)$

$I_L=\underset{K-w}{\overset{K^{-}}{\∫}}C\left(E\right)\mathrm{dE}---\left(21\right)$

已知X射线对已知材料（对比剂）的理论衰减曲线为a(E)，如图1所示。于是，两个能量段的衰减系数μ_R(E)和μ_L(E)，则可表示为

${\mathrm{\μ}}_R\left(E\right)=\frac{1}{w}\underset{K^{+}}{\overset{K+w}{\∫}}a\left(E\right)\mathrm{dE}---\left(22\right)$

${\mathrm{\μ}}_L\left(E\right)=\frac{1}{w}\underset{K-w}{\overset{K^{-}}{\∫}}a\left(E\right)\mathrm{dE}---\left(23\right)$

下面，将利用信号差异噪声比（SDNR）作为在K-edge成像中寻找最佳X射线能量段宽度的准则。这里，信号差异噪声比（SDNR）可定义为

$\mathrm{SDNR}=\frac{{\stackrel{\‾}{\mathrm{\μ}}}_R-{\stackrel{\‾}{\mathrm{\μ}}}_L}{\sqrt{{\mathrm{\σ}}_R^2+{\mathrm{\σ}}_L^2}}=\frac{E\left(f_R\right)-E\left(f_L\right)}{\sqrt{\mathrm{Var}\left(N_{f_R}\right)+\mathrm{Var}\left({N_f}_L\right)}}=\frac{E\left(f_R\right)-E\left(f_L\right)}{\sqrt{\mathrm{Var}\left(f_R\right)+\mathrm{Var}\left(f_L\right)}}---\left(24\right)$

其中，f_R和f_R为K-edge两侧两个X射线能量段重建图像。根据公式(15)、(16)、(22)和(23)，可以利用X射线能量段宽度w，借以表示重建图像f和图像噪声N_f的方差。进而，信号差异噪声比（SDNR），可以用X射线能量段宽度w来表示，即

$\mathrm{SDNR}=\frac{\underset{0}{\overset{\mathrm{\π}}{\∫}}\mathrm{d\θ}{g}_R(w,t^{\′})h(t-t^{\′})d{t}^{\′}-\underset{0}{\overset{\mathrm{\π}}{\∫}}\mathrm{d\θ}\underset{-\∞}{\overset{\∞}{\∫}}{g}_L(w,t^{\′})h(t-t^{\′})d{t}^{\′}}{\sqrt{\underset{0}{\overset{\mathrm{\π}}{\∫}}\mathrm{d\θ}\underset{-\∞}{\overset{\∞}{\∫}}\frac{1}{e^{-g_R(w,t^{\′})\frac{1}{w}\underset{K^{+}}{\overset{K+w}{\∫}}C\left(E\right)\mathrm{dE}}}{h}^2(t-t^{\′})d{t}^{\′}+\underset{0}{\overset{\mathrm{\π}}{\∫}}\mathrm{d\θ}\underset{-\∞}{\overset{\∞}{\∫}}\frac{1}{e^{-g_L(w,t^{\′})}\frac{1}{w}\underset{K-w}{\overset{K^{-}}{\∫}}C\left(E\right)\mathrm{dE}}{h}^2(t-t^{\′})d{t}^{\′}}}---\left(25\right)$

其中，g_R(w,t)和g_L(w,t)分别为K-edge两侧两个X射线能量段重建图像的正弦图，可表示为

$g_R(w,t)=\∫({\mathrm{\μ}}_R(w,t)+{\mathrm{\μ}}_{\mathrm{Rother}}(w,t))\mathrm{dl}---\left(26\right)$

$g_L(w,t)=\∫({\mathrm{\μ}}_L(w,t)+{\mathrm{\μ}}_{\mathrm{Lother}}(w,t))\mathrm{dl}---\left(27\right)$

在公式（26）和（27）中，μ_Rother(w,t)和μ_Lother(w,t)分别表示在K-edge两侧两个X射线能量段的非已知材料（对比剂）区域投影，μ_R(w,t)和μ_L(w,t)分别表示在K-edge两侧两个X射线能量段的已知材料（对比剂）区域投影。

公式（25）还可以看作数学研究中的一个最优化问题，即，如何选取X射线能量段宽度w使信号差异噪声比（SDNR）最大，从而获得最佳X射线能量段宽度。在计算信号差异噪声比（SDNR）之前，需要决定X射线能量段宽度w的取值范围。理论上，在对比剂K-edge处，跳变前的X射线衰减系数μ_L，应该小于跳变后的X射线衰减系数μ_R。由此，可以得到以下的不等式：

${\mathrm{\μ}}_L=\frac{1}{w}\underset{K-w}{\overset{K^{-}}{\∫}}a\left(E\right)\mathrm{dE}\≤{\mathrm{\μ}}_R=\frac{1}{w}\underset{K^{+}}{\overset{K+w}{\∫}}a\left(E\right)\mathrm{dE}---\left(28\right)$

理论上，X射线能量段宽度w最小取值为0，即，在对比剂K-edge处进行成像。X射线能量段宽度w的最大取值，可以根据公式（28）的边界条件计算出。公式（28）的边界条件，可以转化为以下的方程：

$y\left(w\right)={\mathrm{\μ}}_R-{\mathrm{\μ}}_L=\frac{1}{w}\underset{K^{+}}{\overset{K+w}{\∫}}a\left(E\right)\mathrm{dE}-\frac{1}{w}\underset{K-w}{\overset{K^{-}}{\∫}}a\left(E\right)\mathrm{dE}=0---\left(29\right)$

由于函数a(E)是拟合高阶函数，难于从方程（29）中直接求解w的值。本发明利用牛顿迭代法求解，求解过程如下：

根据泰勒公式展开，函数y(w)可有m阶多项式，表示为

$y\left(w\right)=y\left(w_n\right)+y^{\′}\left(w_n\right)(w-w_n)+\frac{y^{\′\′}\left(\mathrm{\ξ}\right)}{2!}{(w-w_n)}^2+...+\frac{y^m\left(\mathrm{\ξ}\right)}{m!}{(w-w_n)}^m---\left(30\right)$

通过对函数y(w)线性近似表达，可以得到牛顿迭代公式

$w_{n+1}=w_n-\frac{y\left(w_n\right)}{y^{\′}\left(w_n\right)},n=\mathrm{0,1},...---\left(31\right)$

将公式（29）代入公式（31），可得到新的迭代形式

$w_{n+1}=w_n-\frac{\underset{K^{+}}{\overset{K+w_n}{\∫}}a\left(x\right)\mathrm{dx}-\underset{K-w_n}{\overset{K^{-}}{\∫}}a\left(x\right)\mathrm{dx}}{a(K+w_n)-a(K-w_n)+\frac{1}{w_n}\underset{K-w_n}{\overset{K^{-}}{\∫}}a\left(x\right)\mathrm{dx}-\frac{1}{w_n}\underset{K^{+}}{\overset{K+w_n}{\∫}}a\left(x\right)\mathrm{dx}}---\left(32\right)$

由此，可以把公式（32）看作一个无约束的优化问题，求解出X射线能量段宽度w最大取值w_max，进而，可以得到X射线能量段宽度w的取值范围（0，w_max）。在得到X射线能量段宽度w的取值范围后，可将公式（25）看作为一个有约束优化问题，在X射线能量段宽度w的取值范围，找到一个最佳的X射线能量段宽度w，使信号差异噪声比（SDNR）最大。

本发明的有益效果是这种K-edge成像方法，能够在保证两个重建能谱CT图像感兴趣区域噪声最小化的同时，可以获得最大的对比度差异，从而达到提高已知材料成像对比度的目的。可以用于生物医学CT成像，使所成图像的病灶区域更清晰，便于医生进行观察诊断。

附图说明

图1为本发明实施例某种已知材料（对比剂）的理论衰减曲线图；

图2为本发明实施例人类胸腔模型图；

图3为本发明实施例胸腔模型组成材料的线性衰减系数曲线图；

图4为本发明实施例能量段宽度(w)与信号差异噪声比（SDNR）关系曲线图；

图5为利用本发明方法得到的人类胸腔模型K-edge成像结果，在对比剂测试区域信号差异噪声比最大时在钆溶液(浓度为0.5%)K-edge两侧重建出的图像，A为K-edge左侧的成像结果，B为K-edge右侧的成像结果。

具体实施方式

下面通过实施例进一步说明本发明，并不因此将本发明限制在所述的实施例范围之中。本领域普通技术人员根据本发明构思，而做出的简单变换，应当在本发明所要求保护的范围内。

以下结合附图，具体说明本发明的构思，以及在此构思下的工作过程。

为了验证基于信号差异噪声比（SDNR）最大化的K-edge成像方法，本发明进行了数值模拟仿真分析。在仿真中，分析了一个数值仿真模型，即，注射有对比剂（钆溶液）的胸腔模型，如图2所示。

构建的该模型，具体描述如下：

该胸腔模型，它是针对临床前研究而设计的。其是参照Forbild胸腔模型（http://www.imp.uni-erlangen.de/forbild）而定义的。该模型图像为400×400的数字矩阵（每个像素为0.05cm×0.05cm），主要包含心脏区域、软组织区域、肺部区域、脊椎骨区域和对比度增强区域（感兴趣区域（Region ofInterest,ROI））。在数值模拟仿真实验中，使用钆溶液（浓度为0.5％）作为对比剂，在心脏内部一个感兴趣区域内分析K-edge成像特点。X射线对该胸腔模型中不同组成材料的线性衰减系数曲线，如图3所示。这些衰减曲线是根据美国国家标准局（National Institute of Standards and Technology,NIST,http://www.nist.gov/pml/data/xraycoef/index.cfm）提供的X射线衰减数据拟合而成的。

本发明进行了数值模拟分析，以此论证本发明提出的优化的K-edge成像方法。模拟分析，可分为以下几个步骤：

①利用公式（22）和（23），给定一个X射线能量段宽度w，计算出对比剂区域的平均衰减系数。

②利用平行束几何关系，对仿真模型图像进行360°的等角度（1°间隔扫描）模拟投影，得到投影正弦图g(x)。

③根据公式（11）和（12），对获取的投影正弦图g(x)，添加高斯噪声。

④通过FBP方法，重建出对比剂K-edge两侧两个X射线能量段的断层图像，并计算出重建图像中对比剂区域的均值和方差。最后，利用公式（25），计算出信号差异噪声比（SDNR）。

⑤在X射线能量段宽度w的取值范围（0，w_max）内，反复执行以上四个步骤，直至计算出一个最佳的X射线能量段宽度w，使信号差异噪声比（SDNR）最大。

利用以上5个步骤，对注射有对比剂（钆溶液）的胸腔模型进行分析。依据胸腔模型中对比剂材料（钆溶液）的X射线吸收特性进行K-edge成像，得到钆溶液K-edge成像中的X射线能量段宽度w与信号差异噪声比（SDNR）对应曲线，如图4所示。最后，总结出了K-edge成像中对比剂区域信号差异噪声比（SDNR）最大时所对应最佳的X射线能量段宽度w，如表1所示。根据表1中列出的最佳的X射线能量段宽度，通过FBP方法，在钆溶液K-edge两侧两个X射线能量段进行成像，其成像结果如图5所示。

从图5可以看出，两个重建能谱CT图像的感兴趣区域（对比剂区域）对比差异较大，且噪声水平较小。

表1对比剂区域信号差异噪声比（SDNR）最大情况下所对应的最佳能量段宽度取值（在K-edge成像中）

Claims (3)

1.一种优化的K-edge成像方法，其特征在于:在K-edge两侧取两个X射线能量段进行K-edge成像，其中所述两个X射线能量段宽度相等，且利用信号差异噪声比作为最优化约束准则，获取使信号差异噪声比最大的最佳的X射线能量段宽度w，以此进行CT成像。

2.根据权利要求1所述一种优化的K-edge成像方法，其特征在于，本方法具体过程如下：

假设多色X射线源产生的不同能量的X射线光子，它是服从函数C(E)分布的，在K-edge两侧，取两个X射线能量段进行成像，且假设两个X射线能量段宽度相等，均为w，于是，两个X射线能量段光子数分别为

$I_R=\underset{K^{+}}{\overset{K+w}{\∫}}C\left(E\right)\mathrm{dE}$ ①

$I_L=\underset{K-w}{\overset{K^{-}}{\∫}}C\left(E\right)\mathrm{dE}$ ②

假设X射线对已知材料的理论衰减曲线为a(E)，于是，在K-edge两侧的两个能量段的衰减系数μ_R(E)和μ_L(E)，则可表示为

${\mathrm{\μ}}_R\left(E\right)=\frac{1}{w}\underset{K^{+}}{\overset{K+w}{\∫}}a\left(E\right)\mathrm{dE}$ ③

${\mathrm{\μ}}_L\left(E\right)=\frac{1}{w}\underset{K-w}{\overset{K^{-}}{\∫}}a\left(E\right)\mathrm{dE}$ ④

利用经典的解析重建算法——滤波反投影（Filtered Backprojection,FBP），可以得到重建图像f，其均值和方差为的数学表达式为

$\begin{array}{c}E\left(f\right)=E\left(\underset{0}{\overset{\mathrm{\π}}{\∫}}\mathrm{d\θ}\underset{-\∞}{\overset{\∞}{\∫}}{S}_{\mathrm{\θ}}\left(t^{\′}\right)h(t-t^{\′})\mathrm{dt}\right)\\ =\underset{0}{\overset{\mathrm{\π}}{\∫}}\mathrm{d\θ}\underset{-\∞}{\overset{\∞}{\∫}}E\left(S_{\mathrm{\θ}}\left(t^{\′}\right)\right)h(t-t^{\′})\mathrm{dt}\\ =\underset{0}{\overset{\mathrm{\π}}{\∫}}\mathrm{d\θ}\underset{-\∞}{\overset{\∞}{\∫}}g\left(t^{\′}\right)h(t-t^{\′})d{t}^{\′}\end{array}$ ⑤

$\begin{array}{c}\mathrm{Var}\left(f\right)=\mathrm{Var}\left(\underset{0}{\overset{\mathrm{\π}}{\∫}}\mathrm{d\θ}\underset{-\∞}{\overset{\∞}{\∫}}{S}_{\mathrm{\θ}}\left(t^{\′}\right)h(t-t^{\′})\mathrm{dt}\right)\\ =\underset{0}{\overset{\mathrm{\π}}{\∫}}\mathrm{d\θ}\underset{-\∞}{\overset{\∞}{\∫}}\mathrm{Var}\left(S_{\mathrm{\θ}}\left(t^{\′}\right)\right)h^2(t-t^{\′})\mathrm{dt}\\ =\underset{0}{\overset{\mathrm{\π}}{\∫}}\mathrm{d\θ}\underset{-\∞}{\overset{\∞}{\∫}}\frac{1}{I_0{e}^{-g\left(t^{\′}\right)}}{h}^2(t-t^{\′})d{t}^{\′}\end{array}$ ⑥

式⑤、⑥中S_θ表示投影正弦图，g表示X射线经过物体后的整体衰减系数，f表示重建图像，E(f)为重建图像的均值，Var(f)为重建图像的方差；

下面，将利用信号差异噪声比作为最优化约束准则，基于此，获取使信号差异噪声比最大的最佳的X射线能量段宽度w，以此进行CT成像；这里，信号差异噪声比定义为

$\mathrm{SDNR}=\frac{{\stackrel{\‾}{\mathrm{\μ}}}_R-{\stackrel{\‾}{\mathrm{\μ}}}_L}{\sqrt{{\mathrm{\σ}}_R^2+{\mathrm{\σ}}_L^2}}=\frac{E\left(f_R\right)-E\left(f_L\right)}{\sqrt{\mathrm{Var}\left(N_{f_R}\right)+\mathrm{Var}\left({N_f}_L\right)}}=\frac{E\left(f_R\right)-E\left(f_L\right)}{\sqrt{\mathrm{Var}\left(f_R\right)+\mathrm{Var}\left(f_L\right)}}$ ⑦

其中，f_R和f_R为K-edge两侧两个X射线能量段重建图像；根据公式③、④、⑤和⑥；利用X射线能量段宽度w，借以表示重建图像f和图像噪声N_f的方差；进而，信号差异噪声比SDNR，可以用X射线能量段宽度w来表示，即

$\mathrm{SDNR}=\frac{\underset{0}{\overset{\mathrm{\π}}{\∫}}\mathrm{d\θ}{g}_R(w,t^{\′})h(t-t^{\′})d{t}^{\′}-\underset{0}{\overset{\mathrm{\π}}{\∫}}\mathrm{d\θ}\underset{-\∞}{\overset{\∞}{\∫}}{g}_L(w,t^{\′})h(t-t^{\′})d{t}^{\′}}{\sqrt{\underset{0}{\overset{\mathrm{\π}}{\∫}}\mathrm{d\θ}\underset{-\∞}{\overset{\∞}{\∫}}\frac{1}{e^{-g_R(w,t^{\′})\frac{1}{w}\underset{K^{+}}{\overset{K+w}{\∫}}C\left(E\right)\mathrm{dE}}}{h}^2(t-t^{\′})d{t}^{\′}+\underset{0}{\overset{\mathrm{\π}}{\∫}}\mathrm{d\θ}\underset{-\∞}{\overset{\∞}{\∫}}\frac{1}{e^{-g_L(w,t^{\′})}\frac{1}{w}\underset{K-w}{\overset{K^{-}}{\∫}}C\left(E\right)\mathrm{dE}}{h}^2(t-t^{\′})d{t}^{\′}}}$ ⑧

其中，g_R(w,t)和g_L(w,t)分别为K-edge两侧两个X射线能量段重建图像的正弦图，可表示为

$g_R(w,t)=\∫({\mathrm{\μ}}_R(w,t)+{\mathrm{\μ}}_{\mathrm{Rother}}(w,t))\mathrm{dl}$ ⑨

$g_L(w,t)=\∫({\mathrm{\μ}}_L(w,t)+{\mathrm{\μ}}_{\mathrm{Lother}}(w,t))\mathrm{dl}$ ⑩

在公式⑨和⑩中，μ_Rother(w,t)和μ_Lother(w,t)分别表示在K-edge两侧两个X射线能量段的非已知材料区域衰减系数。

3.根据权利要求2所述一种优化的K-edge成像方法，其特征在于：在计算信号差异噪声比之前，需要确定X射线能量段宽度的取值范围：

理论上，在已知材料K-edge处，跳变前的X射线衰减系数μ_L小于跳变后的X射线衰减系数μ_R；由此，可以得到以下的不等式：

${\mathrm{\μ}}_L=\frac{1}{w}\underset{K-w}{\overset{K^{-}}{\∫}}a\left(E\right)\mathrm{dE}\≤{\mathrm{\μ}}_R=\frac{1}{w}\underset{K^{+}}{\overset{K+w}{\∫}}a\left(E\right)\mathrm{dE}$

理论上，X射线能量段宽度w最小取值为0，即，在对比剂K-edge处进行成像；X射线能量段宽度w的最大取值，根据公式

的边界条件计算出；公式

的边界条件，可以转化为以下的方程：

$y\left(w\right)={\mathrm{\μ}}_R-{\mathrm{\μ}}_L=\frac{1}{w}\underset{K^{+}}{\overset{K+w}{\∫}}a\left(E\right)\mathrm{dE}-\frac{1}{w}\underset{K-w}{\overset{K^{-}}{\∫}}a\left(E\right)\mathrm{dE}=0$

由于函数a(E)是拟合高阶函数，难于从方程

中直接求解w的值；利用牛顿迭代法求解，求解过程如下：

根据泰勒公式展开，函数y(w)可有m阶多项式，表示为

$y\left(w\right)=y\left(w_n\right)+y^{\′}\left(w_n\right)(w-w_n)+\frac{y^{\′\′}\left(\mathrm{\ξ}\right)}{2!}{(w-w_n)}^2+...+\frac{y^m\left(\mathrm{\ξ}\right)}{m!}{(w-w_n)}^m$

通过对函数y(w)线性近似表达，可以得到牛顿迭代公式

$w_{n+1}=w_n-\frac{y\left(w_n\right)}{y^{\′}\left(w_n\right)},n=\mathrm{0,1},...$

将公式

代入公式

可得到新的迭代形式

$w_{n+1}=w_n-\frac{\underset{K^{+}}{\overset{K+w_n}{\∫}}a\left(x\right)\mathrm{dx}-\underset{K-w_n}{\overset{K^{-}}{\∫}}a\left(x\right)\mathrm{dx}}{a(K+w_n)-a(K-w_n)+\frac{1}{w_n}\underset{K-w_n}{\overset{K^{-}}{\∫}}a\left(x\right)\mathrm{dx}-\frac{1}{w_n}\underset{K^{+}}{\overset{K+w_n}{\∫}}a\left(x\right)\mathrm{dx}}$

由此求解出X射线能量段宽度w最大取值w_max，进而，可以得到X射线能量段宽度w的取值范围（0，w_max）。

Images