附录

1. NIST 数据集简介

PHYSICAL MEASUREMENT LABORATORY

Summary 原文

Table 3：1keV ~ 20MeV 能量下，原子序数 Z=1 ~ 92 的元素，其物质衰减系数 \(\frac{\mu}{\rho}\) 和物质能量吸收系数\(\frac{\mu_{en}}{\rho}\)

Table 4：(常用) 1keV ~ 20MeV 能量下，放射领域化合物 / 混合物，其物质衰减系数 \(\frac{\mu}{\rho}\) 和物质能量吸收系数\(\frac{\mu_{en}}{\rho}\)

Table 1：\(\frac{Z}{A}\)，平均激发能量 \(I\)，密度

Table 2：计算 Table 4 的相关参数

2. 线性衰减系数的拆分

诊断能量范围内，物质的线性衰减系数主要由光电效应和康普顿散射组成：

\[ \begin{equation}\begin{split} \mu(E)&=\rho\times\frac{N_A}{A}\times\left(\frac{k\times Z^m} {E^3}+Z\times\sigma_{KN}(E)\right)\\ &=\frac{\rho Z^mN_A}{A}\times\frac{k}{E^3}+\frac{\rho ZN_A}{A}\times\sigma_{KN}(E)\\ &=\alpha_{PE}\times\frac{k}{E^3}+\alpha_{CS}\times\sigma_{KN}(E) \end{split}\end{equation}\tag{1} \]

其中，左侧为光电效应部分，右侧为康普顿散射部分，KN 为 Klein-Nishina。因此，在诊断能量范围内，\(\mu(E)\) 可由 \(\frac{k}{E^3}\) 和 \(\sigma_{KN}(E)\) 线性组合 (不考虑K-edge) 。

如果我们把 \(\mu(E)\) 看作高维向量：

\[ \vec\mu=(\mu(20keV),\mu(21keV),\cdots,\mu(140keV)) \]

同理，光电反应对应向量为：

\[ \overrightarrow{PE}=\left(\frac{k}{(20keV)^3}, \frac{k}{(21keV)^3},\cdots,\frac{k}{(140keV)^3}\right) \]

康普顿散射对应向量为：

\[ \overrightarrow{CS}=(\sigma_{KN}(20keV),\sigma_{KN}(21keV),\cdots,\sigma_{KN}(140keV)) \]

那么，不同材料的线性衰减系数向量可由 \(\overrightarrow{PE}\) 和 \(\overrightarrow{CS}\) 两个基向量经过线性组合获得：

\[ \vec\mu=\alpha_{PE}\times\overrightarrow{PE}+\alpha_{CS}\times\overrightarrow{CS}\tag{2} \]

即：诊断能量范围内，不同物质的 \(\vec\mu\) 组成的向量的秩为 2 (不考虑 K-edge)。因此，除了使用 \(\overrightarrow{PE}\) 和 \(\overrightarrow{CS}\) 作为基向量，理论上 \(\vec\mu\) 组成的向量内任意两个线性无关的向量也可以做基向量。例如我们选择有机玻璃和铝作为基材料，那么不同材料的线性衰减系数可写作：

\[ \vec\mu=\alpha_{PMMA}\times\vec\mu_{PMMA}+\alpha_{Al}\times\vec\mu_{Al}\tag{3} \]

类似的，其他常用的基材料还有水和骨头等等。

3*. 材料分解

材料分解可以发生在 图像域 和 投影域 两种尺度下：

1. Image-domain 分解：

由 2. 线性衰减系数的拆分 可知，任意物质的衰减系数 \(\vec\mu\) 可以由两种基材料线性表出：

\[ (\vec {\mu_1}, \vec {\mu_2})\cdot (\alpha_1,\alpha_2)=\vec\mu \]

注意区别

由 post-log 直接重建的图像中，每个像素点的值为 \(\mu\)。这并不是上述材料分解的高维向量，而是一个确切的值，理想情况下是某种材料在某个能谱下的加权平均：\(\mu_m=\frac{\int \Omega (E)\mu_{m}(E)dE}{\int \Omega (E) dE }\)。所以，一般对于掌握 100% 信息的 数字体模，才能进行较为精确的图像域材料分解。

那么，基材料的组合系数：

\[ (\alpha_1,\alpha_2)=(\vec {\mu_1}, \vec {\mu_2})^+\cdot\vec\mu \]

注意

图像域材料分解不依赖能谱信息，只需要 物质衰减系数向量 \(\frac{\vec\mu}{\rho}\)。

2. Projection-domain 分解：

提示

当拥有相同 kVp 下，不同能谱的投影数据、或拥有不同 kVp下投影数据，就可以进行投影域的材料分解。由于该方法需要同一扫描对象的不同能谱投影信息，故也称 双 (多) 能材料分解。

同样，由 2. 线性衰减系数的拆分 可知，物质的线性衰减系数可以拆分为基材料的线性组合。使用有机玻璃和铝，\(\mu(\vec r,E)\) 可以拆分为空间和能量依赖两个部分：

\[ \mu(\vec r,E)=\alpha_{PMMA}(\vec r)\times\mu_{PMMA}(E)+\alpha_{Al}(\vec r)\times\mu_{Al}(E)\tag{4} \]

\(\alpha_{PMMA}(\vec r), \alpha_{Al}(\vec r)\)是位于位置 \(\vec r\) 的材料的线性衰减分解系数。因此，x 射线的衰减可写作:

\[ \begin{equation}\begin{split} N(l)&=N_0\int\Omega (E)e^{-\int_l \mu(\vec r,E)dl}dE\\ &=N_0\int\Omega (E)e^{-\int_l \alpha_{PMMA}(\vec r)\times\mu_{PMMA}(E)+\alpha_{Al}(\vec r)\times\mu_{Al}(E)dl}dE\\ &=N_0\int\Omega (E)e^{- \mu_{PMMA}(E)\times P_{l,PMMA}-\mu_{Al}(E)\times P_{l,Al}}dE \end{split}\end{equation}\tag{5} \]

其中，\(P_{l,PMMA}=\int_l \alpha_{PMMA}(\vec r)dl,\ P_{l,Al}=\int_l \alpha_{Al}(\vec r)dl\)，这两个值是对 x 射线穿过路径上各个物质系数 \(\alpha_{PMMA}(\vec r), \alpha_{Al}(\vec r)\) 的积分，与能量无关，满足 Radon 变换。

如果我们能得到每一个 x 射线穿过路径上的 \(P_{l,PMMA},\ P_{l,Al}\) 值，然后用这些值组成正弦图(哪里有正弦图？)来做重建，即可重建出物质内部的系数 \(\alpha_{PMMA}(\vec r), \alpha_{Al}(\vec r)\) 的分布，从而完成材料分解。

Q：怎么获得每一个 x 射线穿过路径上的 \(P_{l,PMMA},\ P_{l,Al}\) 值呢？

A：寻找 post-log 投影数据 \(ln\frac{N_0}{N(l)}\) 和 \(P_{l,m}\) 的函数关系。对于两个不同能谱 \(\Omega_1(E),\Omega_2(E)\) 下 x 射线衰减的数据，有：

\[ \begin{equation}\begin{split} A_1(l)=ln\frac{N_{1,0}}{N_1(l)}=-ln\int\Omega_1 (E)e^{- \mu_{PMMA}(E)\times P_{l,PMMA}-\mu_{Al}(E)\times P_{l,Al}}dE\\ A_2(l)=ln\frac{N_{2,0}}{N_2(l)}=-ln\int\Omega_2 (E)e^{- \mu_{PMMA}(E)\times P_{l,PMMA}-\mu_{Al}(E)\times P_{l,Al}}dE\\ \end{split}\end{equation} \tag{6} \]

其中，\(N_{1,0},N_{2,0}\) 是两个光谱衰减前的信号，\(N_1(l),N_2(l)\) 是衰减后信号。这样就产生了两个方程，原则上可解出 \(P_{l,PMMA},P_{l,Al}\) 两个未知数。

对于给定两个光谱 \(\Omega_1(E),\Omega_2(E)\)，则 \(P_{l,PMMA},P_{l,Al}\) 和 \(A_1(l),A_2(l)\) 存在函数关系：

\[ \begin{equation}\begin{split} P_{l,PMMA}=f_{PMMA}(A_1(l),A_2(l);\Omega_1(E),\Omega_2(E))\\ P_{l,Al}=f_{Al}(A_1(l),A_2(l);\Omega_1(E),\Omega_2(E)) \end{split}\end{equation} \tag{7} \]

那么，对这两个二元函数，分别用二阶泰勒展开来，就可以拟合它们之间的函数关系：

\[ \begin{equation}\begin{split} P_{l,PMMA}=\beta_{1,0}A_1(l)+\beta_{0,1}A_2(l)+\beta_{2,0}A_1^2(l)+\beta_{1,1}A_1(l)A_2(l)+\beta_{0,2}A_2^2(l)\cdots \\ P_{l,Al}=\gamma_{1,0}A_1(l)+\gamma_{0,1}A_2(l)+\gamma_{2,0}A_1^2(l)+\gamma_{1,1}A_1(l)A_2(l)+\gamma_{0,2}A_2^2(l)\cdots \end{split}\end{equation} \tag{8} \]

最后，这些拟合的参数 \(\beta,\gamma\) 分别作用于真实投影数据 \(A_1(l),A_2(l)\)，即可得到 \(P_{l,PMMA},P_{l,Al}\) ，从而完成投影域材料分解。

重新总结一下，进行如下步骤：

• 根据式 (6)，利用多能谱下的 postlog 衰减公式，获取不同厚度 PMMA 和 Al 的衰减：

\[ A_i(L_{PMMA},L_{Al})=-ln\frac{\int \Omega_i (E)e^{- \mu_{PMMA}(E)\times P_{l,PMMA}-\mu_{Al}(E)\times P_{l,Al}} dE}{\int \Omega_i (E) dE}\tag{9} \]

其中 PMMA 长度取 0，20，40，…，200 mm，Al 长度取 0，10，20，30，40，50 mm。将 PMMA 和 Al 长度组合，一共有 \(11\times6=66\) 种组合。因此对于一个光谱，产生 66 个对应关系，两个光谱共 132 个对应关系。

• 根据式 (8) 使用二元多项式拟合 132 个对应关系，求解出 \(\beta,\gamma\)。
• 这些拟合的参数 \(\beta,\gamma\) 逐像素作用于真实投影数据 \(A_1(l),A_2(l)\)，即可得到由 \(P_{l,PMMA},P_{l,Al}\) 构成的投影图像。
• 重建 \(P_{l,PMMA},P_{l,Al}\) 所得图像即为各位置材料关于 PMMA 和 Al 的分解系数 \(\alpha_1(\vec r),\alpha_2(\vec r)\)。

注意

投影域材料分解依赖能谱信息以获得 \(A_i(l)\)，没有能谱的情况下进行准确的材料分解比较困难。不过，实验表明，根据 半值层 (HVL) 估计的能谱也能做出较准确的材料分解。

综上，材料分解分为 图像域 和 投影域 两种方案，分别适用于不同场景。材料分解以获得基材料分解系数 \(\alpha_i(\vec r)\) 或基材料投影长度 \(P_{l,m}\) 为最终目的。这两者与能量无关，满足 Radon 变换，我们认为它们是等价的：

\[ \alpha_i=\Re^{-1}P_{l,m}\tag{10} \]

4. 投影值 post-log 计算

在上述过程中，我们大量使用了带能谱的 postlog 值计算，如式 (9)。而目前广泛使用的 能量积分型探测器 EID，其 post-log 值计算应调整为：

\[ A_i(L_{PMMA},L_{Al})=-ln\frac{\int \Omega_i (E)\cdot E\cdot e^{- \mu_{PMMA}(E)\times P_{l,PMMA}-\mu_{Al}(E)\times P_{l,Al}} dE}{\int \Omega_i (E)\cdot E\ dE}\tag{11} \]

而对于 光子计数探测器 PCD，其 post-log 可仍按式 (9) 计算。

此外，对于仿真投影实验，如材料分解，也需要根据实际探测器种类选择 post-log 计算方式。

即，凡是牵扯到有关能谱的加权计算，都需要根据实际情况做出式 (9) 和 式 (11) 的选择。

5. 能谱估计

当仅有 CT 机器，但缺失能谱信息时，可以通过如下方法估计能谱：

• 根据比尔朗伯定律 (Beer-Lambert law)，假设投影空间仅有金属铝：

\[ I = I_0\int\Omega(E)e^{-\mu_{Al}(E)\times P_{l,Al}}dE\tag{12} \]

\(\Omega(E)\) 为归一化能谱。初始时，\(P_{l_0,Al}=0\)，通过不断增加铝片的厚度，直到 \(\frac{I}{I_0}=\frac{1}{2}\)，此时 \(P_{l,Al}\) 即为铝的 HVL 厚度。

什么是 HVL？

半值层 Half-Value Layer (HVL)：在探测器前不断增加某材料厚度，当 50% 的入射能量衰减时，其材料厚度被称为半值层。相应的还有 fourth-value layer (FVL) ，为 75% 入射能量衰减时的材料厚度等。

• 记录 CT 机器在某恒定 kVp 下的 HVL 和 FVL。则估计的归一化能谱 \(\Omega_0(E)\) 通过调整铝片和铜片等滤波片厚度，以尽可能满足式 (12)，即：

\[ \begin{equation}\begin{split} \int\Omega_0(E)e^{-\mu_{Al}(E)\times HVL}dE=\frac{1}{2}\\ \int\Omega_0(E)e^{-\mu_{Al}(E)\times FVL}dE=\frac{1}{4} \end{split}\end{equation} \tag{13} \]

此时估计的能谱就可以较准确的进行后续计算。

提示

一般模拟能谱的软件，如 SPEKTR 3.0 本身提供能谱的 HVL 等信息，可以略去式 (13) 的计算。