ABSTRACT

To solve the cone beam X-ray scattering problem in the bone mineral density instrument, a Monte Carlo method-based scatter correction algorithm was designed in this paper. First, this article introduced the mathematical principle and formula derivation process of the algorithm, and used Monte Carlo software toolkit simulation to obtain simulation data. Then, the measured image deconvolution scattering kernel function is used to obtain the original image and achieve the effect of scattering correction. The algorithm does not need additional hardware devices, at the same time; it has the advantages of simple calculation and strong operability. Experimental results manifest that this system is accurate, useful and high efficient.

Keywords:Cone Beam X-Ray, Monte Carlo Simulation, Scatter Correction, Deconvolution

应用于骨密度仪的锥束X射线散射校正研究

时旭，赵兴群

东南大学，生物科学与医学工程学院，江苏 南京

收稿日期：2018年4月7日；录用日期：2018年4月20日；发布日期：2018年4月27日

摘 要

为解决骨密度仪中锥束X射线的散射问题，作者设计了一个基于蒙特卡洛方法的散射校正算法。本文首先介绍了该算法的数学原理与公式推导过程，并利用蒙特卡洛软件工具包模拟得到仿真数据。随后将得到的测量图像反卷积散射核函数得到原始图像，达到散射校正的效果。该算法不需要额外的硬件设备，同时计算简单，可操作性强。最后的实验表明该算法散射校正效果良好，具有较强的实用性。

关键词 :锥形束射线，蒙特卡洛模拟，散射校正，反卷积

Copyright © 2018 by authors and Hans Publishers Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY).

http://creativecommons.org/licenses/by/4.0/

1. 引言

双能X射线吸收法(DXA)经历了长期的发展，由于其高精度、低剂量和速度快的优点，是世界卫生组织(WTO)推荐的成熟骨质疏松测量方法。目前，双能X线骨密度仪广泛采用扇束扫描方式。为了提高扫描速度和精度，最新的DXA骨密度仪使用锥束X射线进行骨骼扫描，采用平板探测器进行成像，这种方式效率高、计量小、扫描过程中没有机械和扫描对象的移动，成像效果好 [1] 。然而新的技术总是面临一些瓶颈问题，我们必须解决锥束X射线散射和放大问题才能将此技术普及并推广到市场。

本文提出的算法旨在对X射线的散射现象进行校正。散射来源于X射线与物质原子的相互作用：光电效应、康普顿效应、瑞利散射。在光电效应中，X射线将全部能量传递给原子内的电子，并且导致电子挣脱原子核束缚。在康普顿效应中，X射线光子与自由电子相互作用，导致电子改变传播路径，变成散射光子。在瑞利散射中，入射光子与内层电子发生弹性散射，束缚电子吸收光子从而跃迁到高能级，同时释放一个散射光子，当射线能量低于200千伏时，瑞利散射可以忽略不计。本文所指散射主要由于康普顿效应产生，参考 [2] 。另外，平板探测器本身也会造成图像的不均匀，即中心接收光子较多四角较少，反应到图像上是中心较亮，周围较暗，从下文的图中不难发现。

散射的校正方法分为两大类：基于硬件方法的散射校正和基于软件方法的散射校正。其中基于硬件的散射校正又包括准直器法、过滤器法等，这些方法目的在于减少接触平板探测器的散射光子数量 [3] ，但同时产生了改变射线路径、以及增加硬件设备带来的成本上升和稳定性下降等问题。本文重点关注软件方法，提出并实现了一种新算法，该算法根据测量投影每个像素点信号强度，反卷积适当的核函数，从而降低散射改善图像质量 [4] 。核函数和模拟数据通过蒙特卡洛软件包EGSnrc仿真实验得到。

2. 理论基础

2.1. 公式推导

我们把平板探测器上的测量值记为(Im)，认为测量值是散射(Is)与原始图像(Ip)的和。

$\mathrm{Im}=\text{Ip}+\text{Is}$ (1)

散射部分可以用原始图像(Ip)和散射核函数(hs)的卷积表示。

$\text{Is}=\text{Ip}\ast \text{hs}$ (2)

因为原始图像(Ip)是未知的，我们用包含了原始图像和散射的测量图像近似的代替原始图像。

$\text{Is}\approx \text{Im}\ast \text{hs}$ (3)

所以，原始图像可以表示为测量图减去散射估计。

$\text{Ip}\approx \mathrm{Im}-\left(\mathrm{Im}\ast \text{hs}\right)$ (4)

如果我们将散射核函数(hs)重新定义为包含了散射部分和原始部分的新卷积核(Is)，可以写成：

$\text{Im}=\text{Ip}\ast \text{gs}$ (5)

通过傅里叶变换可以得到：

$F\left(\text{Im}\right)=F\left(\text{Ip}\ast \text{gs}\right)$

$F\left(\text{Im}\right)=F\left(\text{Ip}\right)\times F(gs)$

$F\left(\text{Ip}\right)=\frac{F\left(\text{Im}\right)}{F(gs)}$ (6)

这种方法的优点是不需要对原始图像进行估计。早期的散射校正技术大多采用平方、指数 [5] 或高斯散射 [6] 核形状。经过大量实验后本文改进成了下面的模型：

$\text{gs}\left(r\right)=\{\begin{array}{cc}\frac{s}{\left(1-s\right)2kr}{\text{e}}^{-r/k}& r\ne 0\\ k& r=0\end{array}$ (7)

其中，s和k分别被定义为散射比例和散射半径(X射线平均散射距离)。这两个参数控制了散射核函数的形状 [7] 。当然，这两个参数收到厚度的影响，厚度越厚散射比例越大，散射半径也越远。参数r是在极坐标下的距离，也就是说本文假设散射核是旋转对称的。在实际应用中，需要把极坐标替转换为直角坐标系。在非原点处，得到下面的公式：

$\text{gs}\left(x,y\right)=\frac{s}{\left(1-s\right)2k\sqrt{x^2+y^2}}{\text{e}}^{-\sqrt{x^2+y^2}/k}$ (8)

2.2. 算法流程

算法流程如图1所示，首先读取平板探测器响应得到测量投影图。测量投影图像噪点较多，不利于算法的应用，所以运用维纳滤波对源数据进行预处理。接着对每个像素点的信号强度进行信号估计，正如上文提到的平板探测器造成的图像亮度不均匀的情况，这种情况存在于所有使用平板探测器的图像中。所以，本文模拟了空模体情况下探测器的响应图，利用测量投影图与空模体响应拟合函数在每个点的比值来对信号进行估计。随后根据信号比值的大小把像素分成若干类，在骨密度测量中分成骨和软组织两大类即可，在其他应用中可根据具体情况设计。给分类之后图像匹配对应散射核函数，不同的核函数s和k参数由实验获得。然后把成对的散射核与图像按照公式(6)进行反卷积，最后把所有分类反卷积结果叠加得到散射校正后图像，即原始图。

3. 实验过程

3.1. 射线源模拟

蒙特卡罗方法也称为随机抽样技巧或统计试验方法，是本文获得投影数据主要方法和实验基础。EGSnrc是斯坦福直线加速器中心开发的一款精确计算粒子运输的蒙特卡洛模拟软件。该程序包中包含了多种实现特定功能的子程序，比如BEAMnrc用来模拟医用加速器或球管的几何结构，Dosxyznrc用来计算直角坐标系下三维体模剂量分布，Dosrznrc用来计算极坐标下三维体模剂量分布，PEGS4用来生成计算所需的各种材料密度和吸收率等数据，以上三个子程序也是本文主要使用的部分。

双能X射线吸收法利用高低两种能量的射线来消除软组织对骨密度计算的影响，但使用同样的散射校正方法，这里以80 kev的高能射线为例。仿真过程使用的锥束X射线源如图2所示。除了使用锥束射线源获得被照射模体的测量投影数据之外，还需要获取铅盘实验数据，从而计算出核函数(公式(8))中的s和k两个关键参数。铅盘实验是指将不同大小的铅盘放在射线源与模体之间，铅盘厚度为3 mm以确保

图1. 算法流程图

图2. 锥形束射线源示意图

能吸收非散射射线，测量并找到铅盘半径与圆心处探测器的响应值的线性关系，在拟合得到的直线中，k为直线斜率的倒数，s是以自然对数e为底截距负数为指数的值，该数学模型参考(Seibert 1984)。

3.2. 模体和探测器模拟

本文的实验包括两组，第一组实验由不同厚度的有机玻璃(PMMA)组成，有机玻璃密度与软组织接近，用来模拟软组织。第一组实验也称铅盘实验，目的是获取探测器响应，从而推算出散射核函数。铅盘半径从0.1 cm到0.7 cm共七组实验，PMMA厚度从2 cm、4 cm、6 cm和8 cm共四组。探测器分辨率为0.1 cm，所以探测器原点处的响应实际上就是0.1 cm以内的剂量数据。这里计算出的参数s和k作为基准，可根据实际情况小幅调整。第二组实验则设计为不同密度的材料模拟人体脊柱正位的组织结构，称为仿脊柱模体，使用锥形束照射可得测量投影图，使用平行束射线照射得到模体示意图。当前市场上主流的平板探测器像素大小在100~200微米之间，为了模拟真实的探测器，仿真探测器的分辨率被设置为512 × 512，尺寸为10.24 × 10.24厘米，通过计算得出像素大小为200微米，两组实验的探测器材料均为非晶硒，硒层上有200 um的铜作为过滤层，其他参数见表1。

两组实验的模拟分别在EGSnrc的子程序Dosrznrc和Dosxyznrc中完成，运行之后得到egsdose和3ddose文件包含了所有体素的剂量数据，具体的数据格式见EGSnrc官方手册。值得注意的是，本文开头提到了锥束射线的放大现象，这会导致被测对象的畸变，对测量投影图像产生干扰。幸运的是，Dosxyznrc与Dosrznrc软件包都可以直接得到体模吸收剂量分布数据。为了防止这种放大现象对实验结果的干扰，下文所使用数据为模体吸收的剂量分布。

4. 结果分析

4.1. 散射核计算

本文通过编写Matlab程序对仿真数据进行计算。读取剂量分布的源数据，并进行拟合之后得到图3，其中横坐标为铅盘实验中铅盘半径，纵坐标为对应半径铅盘下原点处探测器响应与没有铅盘时探测器响应比值的对数。图中的点为仿真数据，使用最小二乘法进行拟合，得到结果误差平方和(sse)值为0.0058，确定系数(R-square)值为0.96基本满足线性关系。由拟合得出的直线推算出在80 kev能量4 cm厚有机玻璃模体下，s约为0.35，k约为1.42。因此散射核函数中距离与幅值的关系如图4所示，当r为0时认为幅值为k。

表1. 成像系统参数表

图3. 铅盘半径与信号比关系

图4. 散射核函数距离与幅值关系

4.2. 测量投影处理

图5展示了上文提到的来自平板探测器的亮度不均匀现象，可以观察到图的中心位置亮度较高，四角较暗。图6是图5使用二次多项式曲面拟合的结果，将拟合之后的光滑曲面与测量投影图作比较，从而减少这种现象对图像质量的影响。横坐标是像素点位置，纵坐标代表归一化之后的灰度值。在MATLAB中白色灰度值为1，黑色为0。

图7是仿脊柱模体的示意图，此图由平行束射线照射模体得到，然而即使使用平行束射线，在模体内部仍有散射产生，所以这张图旨在说明模体的形状。中间三块亮度较高的规则矩形区域由仿骨材料组成，其余部分是有机玻璃材料，也就是说最理想图像的中心截面应该是三个等高等宽的脉冲，且曲线上下波动越小越好，最理想的情况为直线，但这是实际应用中很难实现。图8是锥形束射线照射模体得到的测量投影图，本文的算法正是对此图进行校正。图9为图7卷积散射核之后的结果，此结果被认为是散射估计，即本文要校正的散射满足此图的分布。通过此图不难发现图像亮度不均和散射现象。图10是

图5. 空模体探测器响应

图6. 空模体探测器响应拟合图

图7. 模体示意图

图8. 测量投影图

图9. 散射估计图

图10. 校正后结果

最终的处理结果，为了更直观的展示结果，下文采用提取图像中心切面的方式对结果进行具体分析，纵坐标为点的灰度值，横坐标代表点的位置。如图11所示，图中的理想直线通过分别计算图7中软组织和骨部分灰度均值，将均值作为该部分的理想值，此图旨在帮助读者理解模体结构。

图12所展示的是图8的列中心截面，图像中的大量噪声对算法的处理结果不利，所以这里有必要先用维纳滤波进行预处理，如图13所示的去噪效果令人满意，但是，我们仍然可以看到图中仍然存在大量散射噪声，这种散射噪声体现在曲线的剧烈上下波动，并且图中有轻微的中间亮度高两边亮度低的现象存在。图14是图10的中心截面，在对比度下降不大的情况下，对噪声的抑制效果明显，并且对平板探测器的亮度不均匀现象也有一定改善。接下来的工作中需要继续改进优化算法，降低不同材料的交界处的散射。

图11. 理想模体中心切面图

图12. 维纳滤波前图像中心切面

图13. 维纳滤波后图像中心切面

图14. 校正后图像中心切面

4.3. 结束语

本文通过蒙特卡洛仿真实验，在前人的基础上改进了反卷积算法，利用仿真实验数据实现了对散射噪声的抑制，对于使用锥束X射线进行医学成像的领域具有一定的贡献，创新性首先体现在对算法的改进，新加入的维纳滤波对于噪声大的图像预处理效果明显，另外针对平板探测器亮度不均的现象也有一定改善，这归功于设计了恰当的实验，最后对于散射核数学模型的改进真正降低了散射噪声。不过，本文所提算法具有一定局限性，这种局限性体现在对于各项参数要求符合实验具体情况，特别是散射核的两个关键参数。

文章引用

时 旭,赵兴群. 应用于骨密度仪的锥束X射线散射校正研究
Study on Cone Beam X-Ray Scattering Correction Applied to Bone Densitometer[J]. 生物医学, 2018, 08(02): 34-43. https://doi.org/10.12677/HJBM.2018.82005

参考文献