附录

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{array}{}\text{(1)}& \begin{array}{rl}\mu (E)& =\rho \times \frac{N_A}{A}\times (\frac{k\times Z^m}{E^3}+Z\times {\sigma }_{KN}(E))\\ & =\frac{\rho Z^m{N}_A}{A}\times \frac{k}{E^3}+\frac{\rho Z{N}_A}{A}\times {\sigma }_{KN}(E)\\ & ={\alpha }_{PE}\times \frac{k}{E^3}+{\alpha }_{CS}\times {\sigma }_{KN}(E)\end{array}\end{array}$$

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

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

$$\overrightarrow{\mu }=(\mu (20keV),\mu (21keV),\cdots ,\mu (140keV))$$

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

$$\overrightarrow{PE}=(\frac{k}{(20keV{)}^3},\frac{k}{(21keV{)}^3},\cdots ,\frac{k}{(140keV{)}^3})$$

康普顿散射对应向量为：

$$\overrightarrow{CS}=({\sigma }_{KN}(20keV),{\sigma }_{KN}(21keV),\cdots ,{\sigma }_{KN}(140keV))$$

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

$$\begin{array}{}\text{(2)}& \overrightarrow{\mu }={\alpha }_{PE}\times \overrightarrow{PE}+{\alpha }_{CS}\times \overrightarrow{CS}\end{array}$$

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

$$\begin{array}{}\text{(3)}& \overrightarrow{\mu }={\alpha }_{PMMA}\times {\overrightarrow{\mu }}_{PMMA}+{\alpha }_{Al}\times {\overrightarrow{\mu }}_{Al}\end{array}$$

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

3*. 材料分解

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

1. Image-domain 分解：

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

$$(\overrightarrow{{\mu }_1},\overrightarrow{{\mu }_2})\cdot ({\alpha }_1,{\alpha }_2)=\overrightarrow{\mu }$$

注意区别

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

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

$$({\alpha }_1,{\alpha }_2)=(\overrightarrow{{\mu }_1},\overrightarrow{{\mu }_2}{)}^{+}\cdot \overrightarrow{\mu }$$

注意

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

2. Projection-domain 分解：

提示

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

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

$$\begin{array}{}\text{(4)}& \mu (\overrightarrow{r},E)={\alpha }_{PMMA}(\overrightarrow{r})\times {\mu }_{PMMA}(E)+{\alpha }_{Al}(\overrightarrow{r})\times {\mu }_{Al}(E)\end{array}$$

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

$$\begin{array}{}\text{(5)}& \begin{array}{rl}N(l)& =N_0\int \mathrm{\Omega }(E)e^{-{\int }_l\mu (\overrightarrow{r},E)dl}dE\\ & =N_0\int \mathrm{\Omega }(E)e^{-{\int }_l{\alpha }_{PMMA}(\overrightarrow{r})\times {\mu }_{PMMA}(E)+{\alpha }_{Al}(\overrightarrow{r})\times {\mu }_{Al}(E)dl}dE\\ & =N_0\int \mathrm{\Omega }(E)e^{-{\mu }_{PMMA}(E)\times P_{l,PMMA}-{\mu }_{Al}(E)\times P_{l,Al}}dE\end{array}\end{array}$$

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

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

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

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

$$\begin{array}{}\text{(6)}& \begin{array}{r}{A}_1(l)=ln\frac{N_{1,0}}{N_1(l)}=-ln\int {\mathrm{\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 {\mathrm{\Omega }}_2(E)e^{-{\mu }_{PMMA}(E)\times P_{l,PMMA}-{\mu }_{Al}(E)\times P_{l,Al}}dE\end{array}\end{array}$$

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

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

$$\begin{array}{}\text{(7)}& \begin{array}{r}{P}_{l,PMMA}=f_{PMMA}(A_1(l),A_2(l);{\mathrm{\Omega }}_1(E),{\mathrm{\Omega }}_2(E))\\ P_{l,Al}=f_{Al}(A_1(l),A_2(l);{\mathrm{\Omega }}_1(E),{\mathrm{\Omega }}_2(E))\end{array}\end{array}$$

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

$$\begin{array}{}\text{(8)}& \begin{array}{r}{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{array}\end{array}$$

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

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

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

$$\begin{array}{}\text{(9)}& A_i(L_{PMMA},L_{Al})=-ln\frac{\int {\mathrm{\Omega }}_i(E)e^{-{\mu }_{PMMA}(E)\times P_{l,PMMA}-{\mu }_{Al}(E)\times P_{l,Al}}dE}{\int {\mathrm{\Omega }}_i(E)dE}\end{array}$$

其中 PMMA 长度取 0，20，40，…，200 mm，Al 长度取 0，10，20，30，40，50 mm。将 PMMA 和 Al 长度组合，一共有 $11\times 6=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(\overrightarrow{r}),{\alpha }_2(\overrightarrow{r})$。

注意

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

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

$$\begin{array}{}\text{(10)}& {\alpha }_i={\mathrm{\Re }}^{-1}{P}_{l,m}\end{array}$$

4. 投影值 post-log 计算

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

$$\begin{array}{}\text{(11)}& A_i(L_{PMMA},L_{Al})=-ln\frac{\int {\mathrm{\Omega }}_i(E)\cdot E\cdot e^{-{\mu }_{PMMA}(E)\times P_{l,PMMA}-{\mu }_{Al}(E)\times P_{l,Al}}dE}{\int {\mathrm{\Omega }}_i(E)\cdot EdE}\end{array}$$

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

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

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

5. 能谱估计

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

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

$$\begin{array}{}\text{(12)}& I=I_0\int \mathrm{\Omega }(E)e^{-{\mu }_{Al}(E)\times P_{l,Al}}dE\end{array}$$

$\mathrm{\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。则估计的归一化能谱 ${\mathrm{\Omega }}_0(E)$ 通过调整铝片和铜片等滤波片厚度，以尽可能满足式 (12)，即：

$$\begin{array}{}\text{(13)}& \begin{array}{r}\int {\mathrm{\Omega }}_0(E)e^{-{\mu }_{Al}(E)\times HVL}dE=\frac{1}{2}\\ \int {\mathrm{\Omega }}_0(E)e^{-{\mu }_{Al}(E)\times FVL}dE=\frac{1}{4}\end{array}\end{array}$$

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

提示

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