algorithm
1 引言
固体火箭发动机因其结构简单、使用维护方便，广泛应用于军事领域。与液体火箭发动机以及吸气式发动机不同，固体发动机在生产过程中推进剂已直接浇筑于燃烧室内，无法使用拆解的方法对其进行全面检查或检测。因此，如何在不破坏发动机的前提下，对其内部可能存在的缺陷进行快速、准确检测，对于保持导弹武器作战状态、提高导使用安全性、降低日常检测和维护成本具有重要作用。
目前，用于固体火箭发动机无损检测的方法包括超声、射线、激光、红外、微波和ICT[1]。与其他无损检测技术相比，ICT检测技术具有直接得到物体内部结构的数字图像，检测结果的存储、传输、分析和处理方便等突出特点。目前，ICT技术已经广泛应用于各类工业领域，在对小型固体火箭发动机的检测过程中也展现出了良好的检测效果。但对于直径在2 m以上的大型固体火箭发动机，采用ICT检测方法,即使进行一个发动机断层面的重建就需要耗费大量时间，而完成发动机整体结构检测需要完成上千个断层面的重建。ICT检测技术应用于此类大型固体发动机检测,这存在时间及人力成本过大的问题。而在此类大型工件的ICT检测过程中，有时并不需要关心其内部结构只需对工件外部环形区域进行检测即可。史源源等[2]对固体火箭发动的切向CT重建技术进行了研究，采用基于待测物体先验图像的IRR迭代算法(iterative reconstruction reprojection,IRR)实现了切向投影数据的图像重建。熊亮等研究了基于最大似然期望最大化的ML-EM算法(maximum likelihood expectation maximization,ML-EM)。曾理等[4]针对石油管道的成像问题通过引入图像先验信息，提出了基于凸集投影全变分最小化的TVM-POCS算法(total variation minimization- projection onto convex sets,TVM-POCS)得到了较好的成像效果。
2 切向ICT检测方法
传统ICT检测采用全断面扫描成像，完成检测需要扫描并重建固体火箭发动机整个断面。与之相比，切向ICT检测只采集穿过物体外部的射线投影数据，在图像重建过程中也仅重建待测物体外部环形区域。在对大型固体火箭发动机进行检测时，采用切向ICT检测，可以在满足对脱粘等缺陷进行检测的前提下降低对设备硬件条件的要求，提高检测效率，因此具有较高的应用价值。
以窄扇束TR扫描ICT系统为例，其对固体火箭发动机进行检测过程中采用平移/旋转扫描方式。如图1(a)所示，在进行TR扫描过程中，射线源与探测器以某一固定步长从待检发动机一端扫描至另一端，完成一次扫描后射线源与探测器复位，检测系统根据设定参数将待检发动机旋转一固定角度，完成旋转后射线源与探测器再进行一次平移扫描，依此往复，直至得到待检发动机360°范围内的全部扫描数据。而采用切向ICT检测方法时，在射线源与探测器同步平移的过程中，不需要扫过整个发动机截面，只需要扫描发动机壳体部分。射线源与探测器同步平移一定步长后即复位，待检发动机旋转固定角度后射线源与探测器进行下一轮平移扫描，直至得到待检发动机360°范围内的切向扫描数据，如图1(b)所示。
如图1(a)阴影部分所示，对于传统的TR全局扫描方法在图像重建过程中得到的图像为待测发动机完整的断面图像。而切向扫描方法只对待测发动机外部环形区域断面图像进行重建，重建区域如图1(b)阴影部分所示。
图1 窄扇束TR全局扫描/切向扫描示意图
Fig.1 Fan-beam global scan & tangential scan
3 切向ICT图像重建算法原理
切向ICT检测方法存在的主要问题在于如何在采集数据基础上完成待测物体外部环形区域的图像重建。由于相较于全局扫描，切向ICT检测只扫描并采集待测物体的部分投影数据，使得在此基础上的图像重建成为一类投影数据不完备条件下的图像重建问题。其解决难点在于待测物体的投影数据是截断的，虽然理论上存在有唯一解，但是由于计算机数值计算误差以及噪声的影响，重建图像沿径向的边缘会出现严重的伪影。对于切向ICT图像的重建，现有重建算法中较典型的包括:ART代数迭代算法(algebraic reconstruction techniques，ART)，IRR迭代算法，ML-EM算法以及基于凸集投影全变分最小化的TVM-POCS算法。
3.1 ART迭代算法
图像重建过程中可以将原始图像视为一个n×n大小的灰度图像矩阵F，将图像矩阵F按行展开得到图像向量f，图像重建问题即求解如下方程组，即：
(1)
式(1)中：M为投影射线总数；pj为第j条射线穿过待测物体后的投影值；wi, j为投影系数反映第i个像素对第j条射线的影响，在数值上可以取为第j条射线与第i个像素块的相交长度。式(1)写为矩阵形式，有：
Rf=P
(2)
式(2)中： R为M×N系统投影系数矩阵； f为由灰度图像矩阵F展开后的图像向量；P为由对应射线投影值组成的投影数据向量。由于单条射线相交的像素有限，因此系统投影系数矩阵R为一大型稀疏矩阵。因此,方程(2)在求解利用常规的矩阵理论难以进行求解。
ART迭代算法其核心思想就是将系统矩阵R按行拆开，将每一行看作一个n维超平面。然后从初始值开始，依次计算迭代值在每一个超平面上的垂直投影，直到迭代收敛至所有超平面的共同交点。利用ART迭代算法求解方程(1),每次迭代过程可以表示为：
(3)
当完成M次投影得到f(M)为一次完整迭代，之后将f(M)作为初始值代回，重复上述过程直到结果满足终止条件。在ART算法的基础上，Andersen等[5]提出的联合代数重建算法SART(simultaneous algebraic reconstruction techniques，SART)。与ART算法不同，SART算法将初始值并行投影至所有超平面，然后计算所有投影落点的平均值，以此作为下一次投影的迭代值。相较于ART算法，SART算法增强了抗干扰能力,并拥有了进行并行计算的可能。
3.2 IRR迭代算法
IRR 是一种基于滤波反投影的迭代算法，由于外部问题的投影数据不完备，因此需要引入虚拟探测器并对采集得到的投影数据进行补充[6]。IRR 算法的本质在于通过迭代，不断改善对虚拟探测器投影数据的估计，使其更加接近真实情况。虚拟探测器布置示意图如图2。
图2 虚拟探测器布置示意图
Fig.2 Schematic diagram of virtual detector layout
设真实探测器得到的数据为P1，虚拟探测器的数据估计为P2，令P0={p|p∈P1∪P2}。将补完得到的投影数据代入式(2)中，即有：
Rf=P0
(4)
由于虚拟探测器数据的加入， 作为投影数据是完备的，可以利用FBP算法或代数迭代方法均可以得到待测物体的重建图像f1。进一步的将得到的重建图像f1代入式(2)中进行重新投影，得到新的投影数据
进一步，令：
P0=Rf1
(5)
(6)
完成式(4)、式(5)、式(6)，即完成一次迭代过程，由于每次迭代的投影数据均完全包含了原始数据，因此每次迭代后结果不会更加偏离真实情况。重复以上迭代过程直到满足终止条件得到重建图像。
3.3 ML-EM算法
对于求解方程组(1)，ML-EM算法的目标是找到一组 f1使得探测器实际测量得到的数据pi出现的概率最大。理想情况下探测器的测量数据为
二者的关系为:
(7)
式(7)中，ni为测量过程中产生的噪声信号。一般研究认为测量产生的误差噪声近似适用于泊松噪声模型，即：
(8)
式(8)取对数，得到极大似然代价函数，即fi由pi得到的概率为：
(9)
在式(9)基础上求得fi，使得ψ(f) 取得极值，求得的fi即为图像重建的最优估计值。上述方程使用解析方法难以求解，一般采用EM方法迭代求解。其迭代过程数学表示为:
(10)
由于ML-EM 算法在重建过程中既考虑了系统的物理模型，又考虑了观测数据的统计特性，重建出的图像要优于 FBP重建出的图像。虽然ML-EM算法在一定迭代中可以取得不错的效果，但是在实际使用中却存在一个问题，即当观测数据被噪声污染比较严重时，重建图像的效果并不是完全随着迭代的进行而越来越好，而是当迭代次数达到一定值之后，重建图像的质量会下降[6]。
3.4 TVM-POCS算法
TVM-POCS算法最早由Sidky等人应用到稀疏角CT图像重建问题中。Zeng L等将此方法应用于石油管道切向CT投影数据的处理，提出了基于子区域平滑的凸集投影总变差最小化迭代算法。对于灰度图像，其图像的全变分定义为：
(11)
TVM-POCS算法的目的在于求得解，使得其在满足式(2)基础上使得图像的全变分最小，即：
(12)
在实际求解过程中，由于测量误差等因素的影响不一定存在完全满足式(2)的求解结果。因此在实际求解过程中允许求解结果存在一定误差，即：
(13)
TVM-POCS算法的实现包含2个步骤：POCS步主要保证重建结果满足数据保真项；TVM步主要是使重建结果的TV减小。两步交替迭代进行，最终得到重建图像，求解过程如下。
步骤1 数据初始化，令f(0)=0；
步骤2 将初始值代入式(3)中，利用ART算法进行一次完整的迭代运算，得到fART；
步骤3 对fART进行非负约束，即
(14)
步骤4 计算增量因子d，即
(15)
步骤5 计算图像全变分梯度大小，图像fART的全变分梯度大小为：
(16)
步骤6 如式(17)所示，利用vi, j对图像进行全变分最小化修正，其中α为调节因子，即
(17)
步骤7 判断迭代结果是否满足迭代求解过程的终止条件
如满足，则终止迭代，如不满足，则令 f(0)=fTVM，代回步骤2,进行重新计算直到满足迭代条件输出结果。
4 仿真结果与分析
利用仿真模体的计算机仿真扫描数据，可以验证以上图像重建算法在固体火箭发动机切向ICT检测图像重建过程中的有效性并比较以上算法各自的优缺点。
4.1 仿真实验参数设置
本文仿真实验采用Matlab进行编程，程序运行环境采用I7-9750H处理器，RAM为16 GB，利用SART算法进行迭代计算过程中调用GPU加速运算过程。仿真实验主要参数如表1所示。
表1 切向ICT检测仿真实验参数设置
Table 1 Parameter setting of tangential ICT
detection simulation experiment
参数名称数值总体探测器数量750切向探测器数量200发动机单次旋转步长 1°重建图像大小512×512ICT扫描方式TR
仿真模体及其缺陷设置如图3所示，全局ICT检测仿真扫描结果如图4(a)所示，为验证切向ICT重建算法，在断面重建过程中只采用边缘部分探测器采集的数据,所采集数据如图4(b)所示。
图3 仿真模体及其缺陷设置示意图
Fig.3 Simulation model body and defect setting
图4 仿真扫描结果示意图
Fig.4 Simulation experiment data
4.2 图像重建结果
利用图4(b)所示的切向扫描数据，分别采用ART、IRR、ML-EM及TVM-POCS等4种算法，对仿真模体外部环形区域进行重建，重建结果如图5所示。
图5 仿真数据重建效果图
Fig.5 Image reconstruction results of simulation data
如图6所示，将以上重建图像中对应仿真模体的缺陷部分进行局部放大，可以较直观地比较不同成像算法的成像效果。
从图6中可以发现，相较于IRR算法及ML-EM算法，ART算法和TVM-POCS算法在对切向ICT扫描数据进行图像重建时效果更好。由于IRR算法引入了人为设定的虚拟探测器数据，在环形重建区域的切向上产生了严重的伪影。切向伪影的存在使得采用IRR算法重建图像难以对径向拓展的缺陷进行有效检测。ML-EM算法较好地还原了仿真模体内的缺陷特征，但重建图像存在严重的棋盘状噪点，这容易导致被测物体内可能的微小缺陷被噪声掩盖。
图6 仿真模体缺陷成像效果图
Fig.6 Simulation of the imaging effect of defects
in the mold
4.3 仿真结果分析
对于仿真实验得到的重建图像，可以利用图像均方误差(RMSE)、峰值信噪比(PSNR)和结构相似性(SSIM)等3种评价指标对其进行定量的比较分析。将原始仿真模体视为一个n×n大小的灰度图像矩阵F，其第i行第j列的像素灰度值记为fn·i-n+j。图像重建结果同样可以视为一个n×n大小的灰度矩阵F′，其第i行第j列的像素灰度值记为
二者之间的均方误差为：
(18)
重建结果与原始图像的均方误差越大，说明重建结果的误差越大。记图像内所有像素中的最大灰度值为fmax，重建结果与原始图像的峰值信噪比定义为：
(19)
SSIM是一种衡量两幅图像相似度的指标，其计算公式为：
(20)
式(20)中： μ、μ′分别为F、F′的均值；σ2、σ′2分别为F、F′的方差；σ″为F和F′的协方差。设L为图像的灰度范围，则 c1=(0.01L)2，c2=(0.03L)2。
当其他条件相同时，重建结果与原始图像的均方误差越小，峰值信噪比越大，结构相似性越大则说明重建效果越好。4种重建算法与原始图像的均方误差、峰值信噪比和结构相似性如表2所示。
表2 4种算法重建图像的RMSE、PSNR和SSIM
Table 2 RMSE,PSNR and SSIM of four reconstructed images
算法RMSEPSNRSSIMART0.066 354.268 10.738 6IRR0.064 054.991 10.891 2ML-EM0.140 839.207 30.666 1TVM-POCS0.058 556.787 00.936 6
通过比较4种算法所得到的重建图像与原始图像的RMSE、PSNR和SSIM数值可以发现，TVM-POCS算法在以上3种评价指标上均优于其余3种算法，这与图5、图6所示重建图像的直观比较效果是一致的。
对于切向ICT扫描得到的数据，ART算法、IRR算法、ML-EM算法及TVM-POCS算法均可以较好地重建出被测物体断面图像。从图6仿真模体内的模拟脱粘图像上看，切向ICT检测可以满足固体火箭发动机界面脱粘的检测要求。
从上述4种切向ICT重建算法的仿真重建结果上看，TVM-POCS算法的图像重建结果优于其他3种算法。ART算法及ML-EM算法的主要缺点在于图像噪点较多，影响对缺陷的识别判断，IRR算法的主要缺点在于重建图像在切向上伪影严重，由于脱粘在图像上是切向延伸的，因此切向伪影将导致无法根据重建图像准确计算脱粘缺陷的实际大小。
5 结论
将ICT检测技术应用于固体火箭发动机无损检测，对于掌握发动机质量情况和变化规律，保持和提高部队战斗力水平具有重大意义。对于大型固体火箭发动机，切向ICT检测对于发动机壳体各类复合材料之间脱粘问题的检测是可行的。本文主要针对切向ICT图像重建算法进行了梳理和研究，通过仿真实验分析了各个算法的优缺点，可为下一步切向ICT重建算法的研究提供参考。
[1] 陈庆贵,卢洪义,齐强,等.固体火箭发动机界面脱粘切向CT检测[J].固体火箭技术,2016,39(03):347-352.
Chen Q G,Lu H Y,Qi Q,et al.Tangential CT inspection of interface debonding of SRM[J].Journal of Solid Rocket Technology,2016,39(03):347-352.
[2] 史源源,王召巴,金永.大型固体火箭发动机脱粘CT检测技术[J].弹箭与制导学报,2006,26(01):928-933.
Shi Y Y,Wang Z B,Jin Y.CT technoligy for solid rocket motor detachment[J].Journal of Projectiles,Rockets,Missiles and Guidance,2006,26(01):928-933.
[3] 郑伟,邓安华,刘云峰,等.固体火箭发动机衬层与药柱脱粘高能X射线检测技术[J].海军航空工程学院学报,2014,29(04):355-358.
Deng W,Deng A H,Liu Y F,et al.High energy X-ray detection technology of solid rocket Motor lining and cartridge debonding[J].Journal of Naval Aeronautical and Astronautical University,2014,29(04):355-358.
[4] Liu L H,Zeng L.Defects segmentation in digital radiographyimages of castings on a railway freight car[J].Insight,2011,53(07):372-376.
[5] 杨富强,张定华,黄魁东,等.CT不完全投影数据重建算法综述[J].物理学报,2014,63(05):058701.
Yang F Q,Zhang D H,Huang K D,et al.Review of reconstruction algorithms for CT incomplete projection data[J].Journal of physics,2014,63(05):058701.
[6] 郭雨濛.CT外问题的重建模型及算法研究[D].重庆:重庆大学,2019.
Guo Y M.Image reconstruction models and algorithms for exterior problem of CT[D].Chong Qing:Chongqing University,2019.
[7] Sidky E Y,PanX.Image reconstruction in circular cone-beam computedtomography by constrained,total-variationminimization[J].Physic in Medicine and Biology,2008,53(17):4777.
[8] 张芳，张权,崔学英,等.基于片相似性和最大似然期望最大化的低剂量CT重建算法[J].计算机应用,2019,34(09):2716-2719.
Zhang F,Zhang Q,Cui X Y,et al.Low-dose computed tomography reconstruction algorithm based onpatch similarity and maximum likelihood expectation maximization[J].Journal of Computer Applications,2014,34(09):2716-2719.
[9] 郭雨濛,曾理.CT外问题的原理与算法综述[J].中国体视学与图像分析,2017,22(03):320-326.
Guo Y M,Zeng L.Review on the principle and algorithms of the exterior problem of CT[J].Chinese Journal of Stereology and Image Analysis,2017,22(03):320-326.
[10] 陈庆贵,卢洪义,于光辉，等.基于先验图像-压缩感知的CT局部重建算法[J].计算机辅助设计与图形学学报,2016,28(04):632-636.
Chen Q G,Lu H Y,Yu G H,et al.A CT local reconstruction algorithm based on prior image and compressed sensing[J].Journal of Computer-Aided Design & Computer Graphics,2016,28(04):632-636.
[11] 王正安,刘宾,华俊伟.基于标样比对测试的固体火箭发动机工业CT检测工艺[J].宇航材料工艺,2020(04):81-86.
Wang Z A,Liu B,Hua J W.Industrial CT detection technology of solid rocket motor based on standard sample comparison test[J].Aerospace Materials & Technology,2020(04):81-86.
[12] 陈庆贵,卢洪义,周 源.固体火箭发动机窄角扇束CT投影直接重建算法[J].固体火箭技术,2015,38(04):591-600.
Chen Q G,Lu H Y,Zhou Y.Direct reconstruction algorithm for narrow fan-beam CT projections of solid rocket motor[J].Journal of Solid Rocket Technology.2015,38(04):591-600.
[13] Sidky E Y,Kao C M,Pan X C.Accurate image reconstructionfrom few-views and limited-angle data in divergent-beam CT[J].Journal of X-Ray Science and Technology,2006,14(02):119-139.
[14] Shepp L A,Logan B F.The fourier reconstruction of a headsection[J].IEEE Transactions on Nuclear Science,1974,21(03):21-43.
[15] Palakkal S,Prabhu K M M.Poisson image denoising using fastdiscrete curvelet transform and wave form[J].Signal Processing,2012,92(09):2002-2017.
[16] 张顺利,张定华,赵歆波,等.基于最小区域的快速CT图像重建[J].计算机辅助设计与图形学学报,2009,21(02):160-164.
Zhang S L,Zhang D H,Zhao X B,et al.Fast CT imagereconstruction based on minimal reconstruction domain[J].Journal of Computer-Aided Design & Computer Graphics,2009,21(02):160-164.