CN105319228A - 一种用于多源x射线成像的pid控制图像重建方法 - Google Patents

Abstract

本发明涉及一种用于多源X射线成像的PID控制图像重建方法，其特性在于包括下列步骤：(1)在多源X射线层析成像***中，用X射线探测器阵列对X射线源发出的m束射线进行检测，构造测量矩阵M，并计算得到灵敏度矩阵S；(2)构造PID控制重建方法迭代公式，将松弛因子α设为0.1，控制器参数β设为-0.0035：(3)设定迭代次数和容忍误差开始迭代，容忍误差达到要求后停止迭代过程，计算得到重建图像的像素值矩阵Y，进而得到被测场域的物质分布图像。本发明将迭代式图像重建方法与控制理论中控制器的设计方法相结合，得到的图像重建结果比传统的图像重建方法产生的伪影少，视觉效果好，重建精度更高，并且可以应用于其他成像技术的逆问题求解过程中。

Description

一种用于多源X射线成像的PID控制图像重建方法

技术领域

本发明涉及一种用于多源X射线成像的PID控制图像重建方法，属于分布参数测量领域。

背景技术

X射线于1895年由德国科学家伦琴发现。X射线波长短，能量大，穿透能力强。当X射线穿透待测物场时，由于射线中的部分光子会被吸收，射线强度会发生衰减，且衰减程度与被测物场的性质有关。传统的X射线成像***采用扫描式的成像方式，普遍存在结构复杂、造价昂贵等缺点。2010年，Mudde在《PowderTechnology》第199卷55-59页发表的论文“Time-resolvedX-raytomographyofafluidizedbed”中提出一种时间分辨X射线层析成像***，采样频率为2500Hz，能够对直径大于25cm的流化床反应器进行测量，该***有三个X射线源和三个检测阵列组成，均匀分布在流化床周围。该***对流化床中气泡分布具有较好的检测效果。

在多源X射线层析成像中，X射线从射线源发出，穿过被测物场，由检测器探测，可以测得穿透路径上的衰减系数之和。经过多条通道结合，可对被测场域内物质分布图像进行重建。通常可认为测量数据和区域内物质分布的关系为线性。假设灵敏度矩阵为S，测量数据为M，则重建图像像素Y满足：

SY＝M(1)

这里，M是m×1维测量矩阵，矩阵元素为衰减系数在路径上的积分，Y是n×1维灰度数据矩阵，S是m×n维灵敏度矩阵。X射线成像的目的就是求得矩阵Y，求解过程存在以下难点：1)通常灵敏度矩阵S中m小于n，S的逆不存在；2)求解过程是病态的；3)由于X射线的固有特性，测量结果存在波动，会对重建图像产生影响。

常用的逆问题重建方法有LBP，SART，Landweber方法等。LBP方法是最早使用的一种简单成像方法，它将通过某点的所有投影射线进行累加，再反向估算出该点的密度值。从成像观点分析，它是不完全的雷登逆变换通常。SART方法由Andersen等人于1984年在《UltrasonicImaging》第6卷第1期81-94页的论文“Simultaneousalgebraicreconstructiontechnique(SART)：asuperiorimplementationoftheARTalgorithm”中提出，是对ART算法的一种改进。Landweber方法由Landweber于1951年在《AmericanJournalofMathematics》杂志第73卷第615-624页的论文“AnIterationFormulaforFredholmIntegralEquationsoftheFirstKind”中提出。

迭代式的图像重建算法中的迭代过程可以看做是一种反馈。PID控制方法是目前最主要的反馈控制形式，有超过90％的控制闭环采用PID控制。积分、比例、微分反馈过程就是基于对过去、现在、未来的控制误差。本发明提出一种PID控制图像重建方法，将图像重建方法与控制理论中控制器的设计方法相结合，得到的图像重建方法比传统的图像重建方法产生的伪影少，视觉效果好，重建精度更高，并且可以应用于其他成像技术的逆问题求解过程中。

发明内容

本发明的目的在于提出一种用于多源X射线成像的PID控制图像重建方法，包括下列步骤：

步骤一、在多源X射线层析成像***中，用X射线探测器阵列对X射线源发出的m束射线进行检测，得到的m个测量数据组成m×1维测量矩阵M，设Y是n×1维灰度数据矩阵，由下式得到m×n维灵敏度矩阵S

S＝{S_ij},1≤i≤m,1≤j≤n(2)

其中S_ij是第i束X射线通过第j个像素点的长度；

则重建图像像素矩阵Y满足

SY＝M(3)

设Y的第k个估计为Y_k，则Landweber方法通常写为：

Y_k＝Y_k-1+αS^T(M-SY_k-1)(4)

这里，α是控制收敛速度的松弛因子，通常满足：

||αS^TS||＜2(5)

T表示共轭转置，||·||是对矩阵范数。α的值越大，收敛速度越快。

将方程(4)写为z变换形式，有：

Y(z)＝z^-1Y(z)+αS^T(M-S_z ^-1Y(z))(6)

则：

Y(_z)＝[I-(I-αS^TS)z^-1]^-1αS^TM(7)

步骤二、在信号处理中，典型的低通滤波器有：

$F\left(z\right)=\frac{1-\mathrm{\β}}{1-{\mathrm{\βz}}^{-1}}---\left(8\right)$

这里，β>0时，F(z)是典型的低通滤波器；β<0时，F(z)可以看做是高通滤波器。

若将F(z)加入方程(7)的右边，则Landweber算法可以修正为：

$Y\left(z\right)=\frac{1-\mathrm{\&beta;}}{1-{\mathrm{\&beta;z}}^{-1}}{\&lsqb;I-(I-{\mathrm{\&alpha;S}}^{T}S)z^{-1}\&rsqb;}^{-1}{\mathrm{\&alpha;S}}^{T}M---\left(9\right)$

这里，扮演着PID控制器的角色。将方程(9)写为迭代方式，得到PID控制重建算法迭代公式：

Y_k＝[(1-β)I-αS^TS]Y_k-1-β(I-αS^TS)Y_k-2+(1-β)αS^TM

(10)

＝(1-β)Y_k-1-βY_k-2+αS^T[S(βY_k-2-Y_k-1)+(1-β)M]

I是与S^TS同等大小的单位矩阵，α和β两个参数需要进行调整以得到重建算法最佳性能。通过一系列测试，将α设为0.1，β设为-0.0035，此时F(z)体现为微分滤波器的性能。

步骤三、设定最大迭代次数和容忍误差，将测量数据M代入方程(10)进行迭代，当容忍误差达到要求后，停止迭代，得到重建图像的像素值矩阵Y，进而得到被测场域的物质分布图像。

本发明将图像重建方法与控制理论中控制器的设计方法相结合，得到的图像重建方法比传统的图像重建方法产生的伪影少，视觉效果好，重建精度更高，并且可以应用于其他成像技术的逆问题求解过程中。

附图说明

图1为三源X射线成像***结构图；

图2为原始图像和用三种图像重建方法重建得到的图像；

图3用三种图像重建方法重建得到的图像灰度值直方图。

具体实施方式

本发明，即一种用于多源X射线成像的PID控制图像重建方法，包括下列步骤：

步骤一、在多源X射线层析成像***中，用X射线探测器阵列对X射线源发出的m束射线进行检测，得到的m个测量数据组成m×1维测量矩阵M，设Y是n×1维灰度数据矩阵，由下式得到m×n维灵敏度矩阵S

S＝{S_ij},1≤i≤m,1≤j≤n(11)

其中S_ij是第i束X射线通过第j个像素点的长度；

步骤二、由步骤一得到的灵敏度矩阵S，构造PID控制重建方法迭代公式：

Y_k＝[(1-β)I-αS^TS]Y_k-1-β(I-αS^TS)Y_k-2+(1-β)αS^TM

(12)

＝(1-β)Y_k-1-βY_k-2+αS^T[S(βY_k-2-Y_k-1)+(1-β)M]

其中，I是与S^TS同等大小的单位矩阵，Y_k为Y的第k个估计矩阵，T算符表示共轭转置，α是松弛因子，β为控制器参数，将α设为0.1，β设为-0.0035，此时F(z)体现为微分滤波器的性能；

步骤三、设定最大迭代次数和容忍误差，将测量数据M代入方程(12)进行迭代，当容忍误差达到要求后，停止迭代，得到重建图像的像素值矩阵Y，进而得到被测场域的物质分布图像；

为了验证所提出的方法，进行如下实验验证。将三源X射线层析成像***布置在充满Al₂O₃颗粒的流化床上，如图1所示，颗粒直径为0.25mm。而后，将两个内径分别为52mm和23.5mm，壁厚为2mm的中空有机玻璃棒放置于被测场域中，如图2所示。对测量***进行空场和满场标定，而后分别用SART、Landweber和PID控制方法分别进行图像重建。

图像重建结果如图2(a)到(i)所示，可以看到，对比重建图像和原摄影图像，三种图像重建方法均可以重建出两个圆柱物体的位置和大小。但传统的SART和Landweber方法重建结果存在较多伪影，而PID控制重建方法重建物质区域外的伪影明显减少，图像质量更好。重建结果的灰度值直方图如图3所示。这里，尽管原始图像是二值化的，但重建结果如果未加阈值的话，结果中会存在[0,1]之间分布的灰度值。可以看到，PID控制方法可以得到更多的0和1灰度值，体现出更好的重建效果。

以上对本发明及其实施方式的描述，并不局限于此，附图中所示仅是本发明的实施方式之一。在不脱离本发明创造宗旨的情况下，不经创造性地设计出与该技术方案类似的结构或实施例，均属本发明保护范围。

Claims (1)

1.一种用于多源X射线成像的PID控制图像重建方法，其特性在于包括下列步骤：

步骤一、在多源X射线层析成像***中，用X射线探测器阵列对X射线源发出的m束射线进行检测，得到的m个测量数据组成m×1维测量矩阵M，设Y是n×1维灰度数据矩阵，由下式得到m×n维灵敏度矩阵S

S＝{S_ij},1≤i≤m,1≤j≤n(1)

其中S_ij是第i束X射线通过第j个像素点的长度；

步骤二、由步骤一得到的灵敏度矩阵S，构造PID控制重建方法迭代公式：

$\begin{array}{c}{Y}_k=\&lsqb;\left(1-\mathrm{\&beta;}\right)I-{\mathrm{\&alpha;S}}^{T}S\&rsqb;Y_{k-1}-\mathrm{\&beta;}\left(I-{\mathrm{\&alpha;S}}^{T}S\right)Y_{k-2}+\left(1-\mathrm{\&beta;}\right){\mathrm{\&alpha;S}}^{T}M\\ \left(1-\mathrm{\&beta;}\right)Y_{k-1}-{\mathrm{\&beta;Y}}_{k-2}+{\mathrm{\&alpha;S}}^T\&lsqb;S\left({\mathrm{\&beta;Y}}_{k-2}-Y_{k-1}\right)+\left(1-\mathrm{\&beta;}\right)M\&rsqb;\end{array}---\left(2\right)$

其中，I是与S^TS同等大小的单位矩阵，Y_k为Y的第k个估计矩阵，T算符表示共轭转置，α是松弛因子，β为控制器参数，将α设为0.1，β设为-0.0035，此时F(z)体现为微分滤波器的性能；

步骤三、设定最大迭代次数和容忍误差，将测量数据M代入方程(2)进行迭代，当容忍误差达到要求后，停止迭代，得到重建图像的像素值矩阵Y，进而得到被测场域的物质分布图像；

本发明将图像重建方法与控制理论中控制器的设计方法相结合，得到的图像重建方法比传统的图像重建方法产生的伪影少，视觉效果好，重建精度更高，并且可以应用于其他成像技术的逆问题求解过程中。