基于扩充投影矩阵与盲去卷积的代数重建算法

王化祥，薛 倩

(天津大学电气与自动化工程学院，天津 300072)

工业过程中存在大量多相流测量问题，高精度、实时在线测量对生产过程的控制与管理有重要意义．图像重建算法是实现层析成像(computed tomography，CT)系统可视化测量的关键．随着 CT技术的广泛应用，一些新算法如小波变换、神经网络模型、遗传算法等被引入到图像重建中，取得了较好效果．然而对于工业多相流 CT系统，投影数据少导致的图像重构病态问题并未很好地解决．笔者针对设计的多相流 CT系统，改进了投影矩阵计算和迭代重建算法，并通过实验证明该方法的有效性．

1 多相流CT系统简介

多相流 CT系统由γ射线源、探测器阵列、数据采集处理及图像重建单元组成[1]．射线穿过多相流体后衰减，其衰减率与介质吸收系数有关，测量射线的衰减程度即射线投影，经数据处理及图像重构，可再现管内多相流体分布．

如图 1(a)所示，重建区域分成J nn= × 像素，依次编号为 1～J．在第 j像素中，射线吸收系数可视为常数，以像素值xj表示．j像素对i射线投影的贡献权值以加权因子rij表示，则射线i总的射线投影为

式中：I为射线数．为保证成像质量，I、J一般选取较大，所以投影矩阵计算复杂．但对每条射线，它仅与不超过2,n像素相交，因此R为大型稀疏矩阵．

图1 多相流CT系统Fig.1 Multiphase flow CT system

2 扩充投影矩阵

忽略射线宽度，并设其间距为τ，加权因子 rij以2种方法定义．

(2) rij＝i号射线在j号像素内的长度．

针对以上2种定义，提出投影数较少的CT系统扇束投影矩阵的2种算法，并比较二者的成像效果．

2.1 逐点插值算法

(1) 赋R初值为I×J的零矩阵；

(2) 将成像区域内的像素E的直角坐标(x, y)转化成极坐标(r,θ)；图 2中 S0代表射线源， C1C2代表探测器阵列，由S0发出的射线经过E，与C1′C2′交于A点，与C1C2交于B点；Os为C1C2的中点；C1′C2′//C1C2//EF．由△ A OS0～△ E FS0可得O与A间距离为

图2 等距射线扇形束参数关系Fig.2 Relations among parameters of fan beam

若i是 1～17间的整数，则对 rij赋值为 1，若i是1～17间的小数，则令 a = f loor(i)，(floor(i)取i的整数部分)，对 ra+1,j赋值为t = i - a ，而 ra,j＝1- t．

(4) 依次计算全部像素点，循环 J 次．

该方法可以图3解析．由于投影稀疏，在射线1、2间存在大量无射线穿过的像素，若令其加权因子为零，这些像素值将得不到修正，当该像素过多时将无法成像．现假设存在第1.4条射线穿过像素e，则e对射线1投影的贡献权值为 0.6，对射线 2投影的贡献权值为 0.4．这种赋值方法相当于卷积反投影重建算法中的线性插值，因此称其为插值算法．

(3) 由s算出经过像素jx的射线号

2.2 逐线循环算法

如图 3所示，将重建区域划分为 n×n方格，依次编号为1，2，…，n2．若射线斜率的绝对值不大于1，设距离参数 d表示射线与网格的交点到网格底边的距离，d1为当前的d，d2表示在x方向移动1个网格后的d．若射线斜率的绝对值大于 1，则 d为横向距离，每步在y方向移动1个网格，在图中相应表示为 d1′、．5个源与对应探测器阵列的组合依次编为 1～5组，第3组探测器1～17的纵坐标为-8×16～8×16，如图4所示．

图3 成像区域坐标Fig.3 Imaging region on xy-plane

图4 5源模型投影示意Fig.4 Skech map of projections in 5 sources model

步骤1在-128～128之间每隔τ＝2插入1条虚拟射线，总射线数由 5×17＝85增至 5×129＝645，实现了信息扩充，同时以射线穿过像素的长度为投影系数，较式(3)更为精确；

步骤2根据几何模型，计算每条射线方程：第3组的各条射线可由其起点与终点确定其方程，而 1、2、4、5组的各条射线可由第 3组的相应射线旋转相应角度获得，易得其点斜式方程．事实上，由于模型关于x轴对称，只需计算第1、2组以及第3组的一半射线穿过的网格及长度，其他射线穿过的网格号及长度可由其对称射线的计算结果直接获得．例如关于x轴对称的网格Ea、Eb，其关系式为

式中 re m(Ea- 1 ,n)表示取 ( Ea- 1 )/n 的余数．

步骤3根据射线斜率与源位置判断射线i与重建区域相交边界，算出射线i与边界的交点坐标，由此得到射线i穿过的第 1个网格(记为0E)及最后 1个网格(记为1E)；

步骤4射线斜率k分为k＞1，0≤k≤1，-1 ≤ k ＜0，k ＜ -1，分别通过4种不同的循环，由 d1、k、δ求出 d2，根据 d2求出射线穿过的网格号及长度，网格号存入矩阵 G，长度存入矩阵 L．以第 2组k＞ 1 的情况为例，计算流程如图 5所示，其中s =，j初值为1．其他情况可参照编程．

图5 射线i穿过网格号及长度计算流程Fig.5 Flow chart of loop for ray i

步骤 5若当前网格号 E = E1，射线i计算完成，否则计算新的 d2，求下一个网格号和长度；

步骤6依次计算各条射线，循环 I 次；

步骤 7以 G中元素为列标，L 中对应元素为值，将 G重排成 85行的矩阵 R，或用插值法对测量矢量P进行扩充，将射线投影值加权后赋予相邻虚拟射线．

对每组虚拟112条射线，参考平行射束的投影矩阵算法[2-3]计算射线穿过的网格号及被网格所截长度，较之于逐点循环，由于I远少于J，并利用实际系统的对称性简化计算(以 5源系统为例，计算量减少1/2；若系统结构关于x和y轴均对称，则可减少3/4)，极大提高了计算速度，且像素越多，速度提高越明显．

由于平均每条射线穿过n个像素，而上述 2种算法使投影矩阵具有远大于85n个非零值，故可视为投影矩阵的扩充式算法．

3 盲去卷积复原

设 f为原始图像，g为退化图像，h为点扩展函数(point spread function，PSF)，“*”表示卷积，则图像退化模型为

图像复原即由g、h对f进行最佳估计．将 CT系统视为一个图像退化系统，重建图像为退化图像，则可将其代入图像复原算法．由于CT系统的点扩展函数(PSF)未知，故选择盲去卷积算法．盲去卷积仅利用退化图像，即可同时估算PSF和清晰图像[4]．

设f、g均为离散概率频率函数，那么f、g在像素i上的值if，ig可以认为是事件(假定收集到的单位光子为一个事件)在该像素点上发生的频率数．由贝叶斯条件概率公式得

4 改进重建算法

4.1 代数重建算法

Herman提出的代数重建算法(algebrai reconstruction technique，ART)迭代公式为式中：ik= k ( modI )+1，k为迭代次数，ik为射线号；pi为测得的射线i的投影，以探测的粒子数表示[5-6]．

基于式(11)重建，经3I ～ 8I次迭代，结果较好[5].为进一步提高成像质量与实时性，将盲去卷积算法融入代数重建算法的迭代过程(记为 BD-ART)，流程如图6所示．其中，ε为根据精度要求指定的正数，KI为总迭代次数．实验证明该方法经1I ～ 2I次迭代即可取得优于ART的效果．

图6 加入盲去卷积的ART算法流程Fig.6 Flow chart of BD-ART

4.2 同乘代数重建算法

Byrne提出的同乘代数重建算法(SMART)迭代公式为

最小交叉熵算法适用于投影数据较少的场合，重建质量优，但需要非均匀分布的先验信息．在缺乏先验信息情况下将其应用于多相流 CT系统，需进行适当改进．将盲去卷积算法融入 SMART的迭代过程(记为 BD-SMART)，以复原图像为先验信息，进行迭代计算，流程如图7所示．其中，K为迭代次数．

图7 加入盲去卷积的SMART算法流程Fig.7 Flow chart of BD-SMART

4.3 图像处理

图像处理步骤如下．

步骤 1将管道外部的像素赋值成任意背景值，以消除物场外的伪影．

步骤2利用盲去卷积估算的 PSF，对图像进行Richardson-Lucy 复原[9]．

步骤 3根据成像结果，选择适当阈值分割图像[10]．

5 实验结果与分析

计算机配置为Pentium(R)4 2.93GHz CPU、504 MB内存，模拟 4种油水两相流流型(见图 8)并测得投影数据，以MATLAB7.1为平台验证本文方法．

图8 油水两相流流型Fig.8 Flow regimes of oil-water two phase flow

未扩充的投影矩阵记为 R0，逐点插值法计算的扩充投影矩阵记为 R1，逐线循环法计算的扩充投影矩阵记为R2．

图9为采用R0基于ART重建的层流图像，几乎无法由该图像识别流型．对比表 1可见，扩充投影矩阵后成像效果显著改善．

图9 采用R0重建层流图像Fig.9 Laminar flow image reconstructed with R0

像素足够多时，逐线循环法才比逐点插值法的速度明显快，故表 1中，重建像素为 200×200的图像．对比发现，采用 R2明显提高了实时性；加入盲去卷积运算后，成像质量明显提高；虽然盲去卷积运算使迭代耗时增加，但由于所需迭代次数显著减少，成像速度得到提高．与采用R1的ART相比，采用R2的BD-ART在实时性与成像质量方面均有明显优势．

表1 基于ART与BD-ART的成像结果Tab.1 Results obtained using ART and BD-ART

考虑本文中成像区域尺寸与分辨率要求，取 n＝58即可，对于58×58像素的图像，R1与R2耗时相当，而 R1计算简单，故表 2中采用 R1．对比表 2第 2、4行，可见 BD-SMART明显提高了成像质量．第 3行为使用文献[1]提出的3倍信息扩充SMART算法(记为IE-SMART)重建的结果，对比BD-SMART的重建图像，后者重建的伪影更少，速度更快，而前者使图像边缘趋于平滑．将2种方法结合(记为IE-BDSMART)，对较复杂流型重建的图像质量更佳，但实时性有所下降．

BD-SMART的重建结果经图像处理后，提取二值图像面积，计算截面的分相含率，如表3所示．

表2 基于SMART与改进SMART的成像结果Tab.2 Results obtained using SMART and improved SMART

表3 由重建图像计算分相含率Tab.3 Phase volume fraction obtained with reconstructed images%

6 结 语

本文提出针对5源17探测器的工业多相流CT系统的扩充投影矩阵计算，并将盲去卷积算法融入代数重建算法与同乘代数重建算法的迭代过程．利用改进的算法对多相流 CT系统进行图像重建，将采用相同投影矩阵不同迭代过程、不同投影矩阵相同迭代过程的重建图像精度和收敛速度进行比较．实验结果表明，本文方法利用少量投影数据即可取得满意的成像效果，能够准确识别流型，重建图像的分相含率误差大致在±1%以内．

［1］王化祥，王 琦，郝魁红. 最小交叉熵图像重建算法[J]. 仪器仪表学报，2009，30(12)：2574-2579.

Wang Huaxiang，Wang Qi，Hao Kuihong. Reconstruction images based on minimum cross-entropy[J].Chinese Journal of Scientific Instrument，2009，30(12)：2574-2579(in Chinese).

［2］Herman G，Meyer L. Algebraic reconstruction can be made computationally efficient[J]. IEEE Transactions onMedical Imaging，1993，12(3)：600-609.

［3］张顺利，张定华，李 山，等. ART算法快速图像重建研究[J]. 计算机工程与应用，2006，42(24)：1-3.

Zhang Shunli，Zhang Dinghua，Li Shan，et al. Research of fast image reconstruction on ART algorithm[J].Computer Engineering and Applications，2006，42(24)：1-3(in Chinese).

［4］Levente Kovács，Tamás Szirányi. Relative focus map estimation using blind deconvolution[J].Optics Letters，2005，30(22)：3021-3023.

［5］庄天戈. CT原理与算法[M]. 上海：上海交通大学出版社，1992.

Zhuang Tiange.Principles and Algorithms of CT[M].Shanghai：Shanghai Jiao Tong University Press，1992(in Chinese).

［6］Herman G T，Lent A. Iterative reconstruction algorithm[J].Computers in Biology and Medicine，1976，6(4)：273-294.

［7］Byrne C L. Iterative image reconstruction algorithms based on cross-entropy minimization[J].IEEE Transactions on Image Processing，1993，2(1)：96-103.

［8］Boer P T，Kroese D P，Mannor S，et al. A tutorial on the cross-entropy method[J].Annals of Operations Research，2005，134(1)：19-67.

［9］Masschaele B，Dierick M，Van Hoorebeke L，et al. Neutron CT enhancement by iterative de-blurring of neutron transmission images[J].Nuclear Instruments and Methods in Physics Research A，2005，542(1/2/3)：361-366.

［10］Golosio B，Masala G L，Piccioli A，et al. A novel multithreshold method for nodule detection in lung CT[J].Medical Physics，2009，36(8)：3607-3618.