CT与DR双能X射线物质识别算法实现与应用（工业选煤、稀土分拣、毒爆检测、垃圾分类等）

CT与DR双能X射线物质识别算法实现与应用（工业选煤、稀土分拣、毒爆检测、垃圾分类等）

       算法涉及到的理论和模型，来自CNKI相关技术论文（不包括本人改进部分），如涉及版权问题，请联系本人。

        物质识别的技术很多，比如质谱仪，色谱仪（利用物质原子/离子改变火焰颜色，通过光谱分辨物质），分析化学等各种方法，但这些物质识别方法都要一定的专业知识，使用方法也不方便，大部分需要采集部分样品进行分析，远不如双能射线非入侵式简便。笔者长期进行双能射线物质识别算法的应用研究，开发多种应用设备，识别的效果还算可以，但依然有改进之处，本人写下博文希望抛砖引玉，共同从事这方面工作的高手不吝赐教，提出宝贵的意见。下面介绍我在双能物质识别算法中使用过的方法，并展现工业选煤的效果。

依据的理论 X射线通过物体后，由于物理效应（康普顿效应、电子对效应、光电效应、相干散射等）会发生衰减，衰减规律服从比尔律：

---①，---②

其中，μ衰减系数（与Z，ρ和E有关），l为物体的长度（厚度），I和I0为射线穿透物体前和穿透物体后衰减后的能量，e是自然底数（为什么是e，我也不知道，世界简单又复杂）。

双能能进行物质识别的根本原因是：相同厚度不同物质对同一射线的衰减量不一样；等效厚度一致的不同材料的物体，对另一种能量的射线衰减量也不一样。衰减是能量E和物质Z，厚度ρ的函数。Z用于物质判别（色调H），ρ用于物质形态（厚度和形态，比如液态水和水蒸气）判别（饱和度V）。

CT和DR在物质识别算法上的不同  CT和DR成像对于物质识别算法最大的不同是：DR成像是探测器得到物体射线方向总厚度上的衰减系数的线积分，基于线积分做物质分解（等效原子系数Z和厚度ρ）；CT成像是探测器得到物体射线方向总厚度上的衰减系数的线积分，但是把线积分通过重建算法，重新投射到物体对应的空间位置上（等效于积分逆过程），基于重新投射的衰减系数做物质分解（等效原子系数Z和电子密度ρ），这样可以避免重叠问题。

在单向投影上的物质识别方法  在KV能级条件下，衰减规律主要是康普顿效应和光电效应（查阅相应文献），MV能级条件下，衰减规律主要是康普顿效应和电子对效应（查阅相应文献）。所以，在KV能级下和MV级能级下，我们的DR物质识别算法不一样。KV能级主要是射线管产生，比如我们的乘用车检查系统，安检机检查系统，MV能级主要是直线加速器或者回旋加速器产生，比如我们的固定式成像检查系统，组合移动式检查系统和车载式成像检查系统。

下面，我以KV能级为例（CT机的射线源跟乘用车的射线源一样，都是基于射线管产生，能级在100KV~350KV之间）说下我们的物质识别方法。

在KV能级下，主要产生的效应是康普顿效应和光电效应。衰减系数可以表示为

-----------------③

其中，fK（E）为只与能量有关与材料无关的康普顿效应作用分量，知晓能量E就可以直接算出，具体的公式可以通过相关文献查询（项目中我们用10.0）。fD（E）为只与能量有关，与材料无关的光电效应作用分量，知晓能量E就可以直接算出，具体的公式可以通过相关文献查询（项目中我们用0.05）。为啥用加号把这两个效应连起来，请参考比尔定律。ak和ad又叫康普顿和光电效应系数。

根据④⑤两公式，我们可以算出有效原子系数：--⑥

等效电子密度怎么表述？很简单，随便取一个分量，我们取简单一点的ak（DR成像中，我们直接取高能图当做电子密度，这么牛X，把我自己都震惊了，为什么这样取，后面会说到，DR直接用线积分），则可以得到电子密度：---⑦，其中W1,W2是系数（现在不用关心，标定的时候会对所有的系数自动调整）。

到了这一步，我们万事俱备，只欠东风，东风是啥，就是ak和ad。我们怎么通过高低能探测器数据求出这两个值呢，这又回到牛逼的比尔定律②！

先说一下探测器，我们用的探测器是闪烁体探测器，闪烁体探测器是一种射线强度探测器，而非时间积分探测器。理解这一点很重要，在计算的时候，我们不考虑时间，探测器得到的值就是衰减系数的线积分量，根据比尔定律②，我们得到如下公式：

---⑧

---⑨

和分别是高能和低能探测器得到的值，尽管不是严格的等号关系，但一定是常系数的线性关系！

公式⑧⑨可以简化成如下公式：

---⑩⑪

通过解线性方程组⑩⑪，得到Ak和Ad。这两个值分别是ak和ad的线积分，在射线与物体相交的长度的积分：，。在DR成像时，由于直接把Ak当电子密度（我直接取高能数据），Ak是ak的线积分，不用进行重新分布计算，所以DR的先天缺陷就是电子密度（厚度）重叠，这点跟CT有本质的不一样。Ak和AD的比值就是等效原子系数Zeff，Ak就是等效原子密度ρe。得到Zeff和ρe。

标定 为什么要标定，有这么精确的数学模型，不是直接可以计算出Z和ρ吗？太天真了！理论上是这样的，但真实世界非常复杂。比如对于上面出现的K1、K2、W1、W2等各种系数，我们都无从得知，更有很多我们未知的原因（公式的变异、能谱的不稳定等等），在这种情况下，我们计算得出的原子系数和电子密度，早就放飞自我，距离真实的原子系数和电子密度相差很大，所以需要借助于测试体，重新刻度得到的Zeff和ρ值。从而根据创建的映射表，查询得到物质识别信息。只有重新刻度后，才能进行三分类物质识别。

   标定刻度的方法 ：略

实测进行物质识别  ：

矿井中挖出的煤炭（含有优质煤炭和矸石）              物质识别算法标记的优质煤（橙色）和矸石（绿色）

                        

                                                                 双能效果图

CT上应用的物质识别方法

1. 1双效应分解模型

CT机的射线源跟乘用车的射线源一样，都是基于射线管产生，能级在100KV~350KV之间。在这个能级范围内，双效应主要是康普顿效应和光电效应。如③式所示。CT的某个断面的一个方向（平行束180*N个投影图，扇形束（180+2α）*N个投影图）上的投影图，经过⑩⑪进行计算，得到Ak和Ad，重复对某个断面所有方向上的投影进行计算，分别得到序列Ak和Ad矩阵，然后根据某个断面的所有角度的Ak和Ad进行反投影重建，也就是Ak和Ad进行反积分得到ak和ad的分布。致此，得到了某个断面的ak和ad分布图，根据⑥⑦式得到Zeff和ρ值图。此时的Zeff和ρ值图都是相对的量，需要进行刻度，如何刻度，请看下一节。

标定刻度流程：略

实测进行物质识别：略

1. 2基物质分解模型

基物质分解原理上类似于双效应分解，基物质分解模型描述为物质的线性衰减系数可以表示成两个已知线性衰减系数的材料的线性组合：

---⑪

为什么选取碳和铁，后面会知道，如果不选取最大和最小原子系数的分解基，得到分解系数大概率出现负数，给反投影带来不利！

Kc和Kfe分别对应基物质的分解系数，后面需要计算的就是这两个参数，μc(E)和μAl(E)分别表示基物质在能量E下的线性衰减系数。μc(E)和μAl(E)分别是物质密度与质量衰减系数的乘积，质量衰减系数可以通过不同厚度测试体分别测出。因此这两个量一旦确定基物质，这μc(E)和μAl(E)值就确定了。

跟双效应分解方法类似，根据比尔定律②，我们得到如下公式：

---⑫

---⑬

和分别是高能和低能探测器得到的值，尽管不是严格的等号关系，但也一定是常系数的线性关系。

公式⑫⑬可以简化成如下公式：

---⑭⑮

通过解线性方程组⑭⑮，得到Kc和KFe。这两个值分别是kc和kfe的线积分.

在CT某个断面图像的各个角度的投影（CT断面重建前的投影数据），通过重建算法，进行反积分，把Kc和KFe重新投影到断面图像的空间上，得到kc和kfe断面图像数据。电子密度ρe和等效原子系数公式如下：

---⑯

---⑰

其中,ρ为碳的物理密度，Z为碳原子系数，M为碳原子质量。n取4.

得到kc和kfe断面图像数据通过⑯⑰计算得到电子密度和有效原子系数图。下一步进行标定刻度，方法和基于效应分解的一模一样！

标定刻度流程：略

实测进行物质识别：

CT三维双能图像

 欢迎技术交流，共同改进。