烟支3D-μCT高质量成像及微观结构表征方法研究

王 亮，王天一，银董红，谢国勇

(1.湖南中烟工业有限责任公司技术中心，长沙 410007;2.北京航空航天大学，北京 100191)

0 引言

降低卷烟的焦油及其有害成分释放量是烟草行业高质量发展的必然要求，烟支结构中的烟丝分布差异对其燃烧过程中烟气的生成和传递具有重要的影响。因此，深入研究烟支微观结构信息，对改善和优化卷烟结构及加工工艺，进而实现降焦减害和提升卷烟品质具有指导作用。目前，烟草行业对卷烟及烟气的分析检测主要以化学分析和物理测量为主，包括卷烟燃烧后的烟气成分分析和原辅材料检测等[1-5]。然而，有关烟支结构的现代表征技术鲜有文献报道，因此，开展烟支微观结构表征方法研究对指导卷烟设计具有重要的应用价值。

近些年来，X射线3D-μCT(Three Dimensional Micro Computed Tomography)技术在无损检测、材料分析和生物医学领域的应用日益得到重视。3D-μCT技术利用微焦点X射线源和高投影放大比保证了高空间分辨率的信息重建，已被公认是替代传统破坏性层析方法的高端技术[6-8]。该技术不受被检物材质与内部复杂程度的限制，通过对被检测对象的一次扫描，获得样品不同视角下的高分辨率DR(Digital Radiography)图像，利用专用重建算法反演出样品内部结构的各向等分辨率三维图像。将3D-μCT技术应用于烟支内部微观结构数字化表征，提取微观物理特征量，可为研究烟支微观结构提供技术手段。

由于3D-μCT系统采用微焦点X射线源，其管电流仅为微安级，这就使得穿过物体到达探测器的射线光子数太少。以上原因最终导致重建图像的噪声较大，有用信息对比度低等问题。解决问题的办法就是抑制CT图像的噪声来提高图像的信噪比。目前，图像去噪的方法很多，较经典的如 Lu H等提出的利用惩罚似然估计法来平滑投影图[9]， Demirkaya K等提出各向异性扩散滤波方法[10]，Hsieh J提出根据噪声的局部特性自适应动态调整滤波器参数的滤波方法[11]。现有的图像降噪算法大多需要用到复杂的迭代运算、运行速度慢和迭代过程易受数据噪声和初值的影响而产生不收敛的情形，工程实用性不强。如何在降低卷烟图像噪声的同时也能很好地保留图像细节特征，提升降噪算法的鲁棒性及运行效率，是本文需要解决的问题。

研究表明，对于不同配方的卷烟，内部微观结构不尽相同，从而导致烟支物理性质乃至烟气质量的差异[12]。然而，目前对于烟支内部微观结构的研究相对较少，手段也较为单一。比如王亮等采用包埋切片法[13]，利用显微照相获得烟支断面的光学图像。但该方法为一次性破坏切片，存在较大局限性，而凭借3D-μCT技术获取的断层图像来对烟支进行特征无损提取则可有效避免上述问题。

综上所述，本文首先设计针对烟支类长物体的螺旋CT扫描重建方法来获得烟支的断层图像，根据断层图像的噪声特征和工程应用要求合理使用BM3D降噪算法得到清晰的断层图像。通过微观形态学方法和图像处理技术，基于烟丝断层图像提取烟丝的骨架结构和孔隙率等特征量，由此获得了较为精确的烟支微观结构，为改善烟草加工工艺提供了物理依据，对提高卷烟质量研究贡献了新的思路。

1 针对烟支长物体的螺旋CT扫描重建

对于烟支这种长物体而言，3D-μCT扫描会面临两个难题：其一，由于采用了高投影几何放大比，锥束角一般在20°以上，重建图像的纵向模糊不可避免；其二，扫描系统的探测器成像高度有限，对烟支进行高放大比成像时会造成烟支高度方向(纵向)上的投影数据截断，不可能利用一次扫描获得的纵向截断投影数据重构出烟支的全局结构。这两个问题的解决要依赖于重建算法和扫描方式的改进。为了实现对烟支的高分辨率重建，需要采用较大的投影放大比进行成像，本文采用螺旋锥束CT，即采用螺旋扫描轨迹来扩展锥束CT的纵向视野。利用样品同步旋转+升降的方式实现螺旋轨迹扫描，完成对烟支高度的全部覆盖。在重建方法的选择上，比较经典的螺旋扫描重建算法有FDK近似重建算法和Katsevich精确重建算法等。作为近似重建方法，螺旋FDK算法在锥角较小时(锥角小于5°)的重建效果可以接受。但在锥角较大时，远离中心的断层图像伪影严重，CT图像的质量会随着锥角或螺距的增大而变差。相比之下，Katsevich精确重建算法[14]的图像质量不受锥角大小的影响。本论文中实验所用的X射线源型号为X-ray XWT-240-CT，射线源的锥角为30°，远大于5°，为大锥角；此外，本文在烟支螺旋扫描时选取了较大的螺距。所以本文在重建方法上选取了Katsevich算法。

本文采用基于平板探测器的CBCT(Cone-beam CT)系统，图2所示为扫描系统的空间坐标系。对于螺旋轨迹线内部的待重建区域，以螺旋轨迹的旋转轴为z轴，螺旋轨迹线的起始点所在水平面与z轴的交点作为原点建立笛卡尔坐标系o-xyz，待重建体素点的坐标为x=[x0,x1,x2]，射线源焦点的位置坐标为λ(s)，射线源焦点到探测器平面的距离为D，螺旋轨迹的旋转半径为R。螺旋轨迹线的表示方法为：

图2 螺旋扫描轨迹的CBCT系统的空间坐标系

(1)

2 高成像放大比下烟支CT图像降噪方法

由投影数据重建得到的CT图像信噪比偏低，重建出的断层图像及三维体数据容易出现失真，不利于后续提取烟支的结构信息。因此，应该对CT断层序列图像进行有效的降噪，通过图像后处理方法提高重建图像信噪比。

本文采用BM3D算法实现对卷烟图像的降噪，该算法有效地结合了非局部降噪和频域降噪两种思想，将相似的图像块组成三维图像组，利用三维图像组在频域具有较高稀疏性的特性减少噪声，不仅可以得到较高的峰值信噪比，并且具有较好的视觉效果。BM3D算法主要分为两个主要阶段:第一阶段产生降噪图像的基本估计,该阶段基于硬阈值处理，并且可以表现出相对良好的降噪性能，所以在要求不高的场景下可以单独使用;第二阶段是基于经验的维纳滤波，使用基本估计和噪声图像产生最终的降噪图像，进一步提高了降噪性能。两个阶段的处理流程基本相似，图3为BM3D算法的流程图[15]。

图3 BM3D算法流程图

BM3D算法通过块匹配的方法形成三维图像组，块匹配分组的意义在于分组之后可以对每个三维图像组使用更高维度的滤波处理，块匹配的过程类似于NL-Means算法[16]。将当前处理的图像块作为参考块，在图像中搜索与参考块相似的图像块，将相似的图像块堆叠起来形成三维图像组，块匹配分组的示意图如图4所示。根据BM3D算法的原理，相似图像块的搜索范围应该是整幅图像，但在实际应用中，考虑到整幅图像搜索效率低下的问题，对于给定的参考块，仅在一个固定大小的搜索窗口中进行搜索，可以大大降低搜索算法的复杂度，提高算法执行效率，搜索窗口的大小可根据实际的情况进行选取。

图4 三维块匹配示意图

协同滤波是BM3D算法中最为核心的一步，用来处理三维图像组。包括三个步骤：

(2)通过硬阈值操作，收缩变换频谱；

聚合步骤是将协同滤波后的图像进行整合，得到最终的降噪图像。在图像块匹配分组的过程中，最终得到的三维图像组有可能是相交的，即同一个图像块可能会被匹配分组到多个图像组中，因此，图像中的一个像素可能会对应多个估计值，这就需要通过加权平均来确定最终的估计值，像素x的基础估计如公式(2)所示。

(2)

式中，X(P,x)为指示函数，若像素x属于图像块P，则值为1，否则值为0；函数uP,R(x)为像素x在滤波后图像块P中的值。第二步的聚合与第一步基本相同，只是加权的权重取决于维纳滤波的系数和噪声强度。

3 烟支特征量提取

3.1 烟丝骨架提取

表达平面区域结构形状的一种重要方法是把它简化成图形，这种简化可以通过骨架提取算法(也称细化算法)实现。烟支骨架的提取有利于深入了解烟丝的拓扑学结构，为优化烟支工艺提供理论基础。骨架提取细化算法输出的结果应满足如下要求[17]：(a)输出单像素细线，(b)细化结果为目标图像的中心，(c)不能破坏连通性，(d)不可删除端点，(e)具有较好的稳定性。

Zhang快速并行骨架提取方法对二值图像进行逐层删除边界点[18]，直到没有再可删除的点时，骨架提取结束。Zhang并行快速细化方法规定的骨架提取方法原理为：在如图5所示的3×3邻域窗口内，以点P为中心，P1到P8顺时针定义，并使用式(3)对该邻域的像素进行约束，N(P)指P1到P8中非零点的数目，Z0指P1到P8从0到1的变化次数，同时满足上述条件的点P就被删除。

P8P1P2P7PP3P6P5P4

图5Zhang快速并行细化算法模板

(3)

Zhang快速并行骨架提取方法分两步执行：

Step1：在图4所示的邻域中执行式(3)，将满足式(3)所有要求的点进行标记，当所有点检验完毕后，将所标记点删除。

Step2：将式(3)中的(c)、(d)换成P3×P5×P7=0、P1×P5×P7=0，同样对所有像素点进行检验，删除同时满足这四个条件的像素点。以上两个步骤完成一次迭代，直到没有像素点再满足式(3)的标记条件，迭代结束，骨架提取完毕。

3.2 骨架的彻底细化

利用Zhang并行快速细化方法提取的骨架，虽然初步满足了骨架提取的相关原则，但是这种方法的提取结果并不能保证满足单像素宽度的要求，这也说明这种方法在骨架提取的过程中存在死点，不能彻底细化。但是，通过对其算法分析可知，这些像素点又不满足上述骨架提取算法的删除条件，最终使不必要的像素点保留下来。为了使提取的骨架满足单像素宽度的要求，需要对提取的骨架进行彻底细化。本文采用纹线追踪方式[19]对非单像素宽度的骨架进行彻底细化，图6(a)为一段放大的骨架，可以看出烟丝的骨架具有一定的方向性，通过进一步判断目标像素3×3邻域的外围骨架点的属性，可以估计纹线的走向，因此采用如图6(b)所示的5×5的彻底细化模板，以图6(a)的骨架为例，骨架点P有P8和P1两点相连，而P4属于正确的骨架走向，对于P8点，纹线的连续可以向Q1～Q3及Q15、Q16方向延伸，而对于点P1，则可以向Q2～Q4方向延伸，而对于3×3邻域的外围，只有点P8存在另一个骨架点Q1，而点P1的外围不存在可能的骨架点。因此，可以判断骨架的走向应是沿着Q1、P8、P、P4、Q9，点P1为不彻底细化的结果，应该删除。

彻底细化算法的执行步骤可概括为：

(1)逐行扫描骨架图像，直至扫描到骨架点Pi；

(2)对点Pi统计其3×3邻域内的骨架点数目n，如果n≤2，则点Pi不存在不彻底细化情况，退出，执行步骤1；如果n≥3，执行步骤3；

(3)此时需要从P1-P2，P2-P3，…P8-P1连续扫描点Pi的邻域，确定是否存在连接的骨架点，如果不存在，则说明此处不存在细化不彻底的现象，退出，执行步骤1；如果存在Pj-Pj+1或P8-P1连接，执行步骤4；

(4)扩展到点Pi的5×5邻域，判断点Pj的外围邻域是否存在骨架点，如果存在，则可判断骨架应该沿着这个方向延伸；否则判断该点为不彻底细化结果，应该删除；

(5)所有骨架点扫描结束，彻底细化算法执行完毕。

3.3 孔隙率的提取

对于烟支这种特殊的多孔介质而言，孔隙率是极为重要的物理参数之一。孔隙率是指多孔介质内微小空隙的总体积与多孔介质总体积的比值，是无量纲的物理量。孔隙率与多孔介质的形状、结构和排列有关。孔隙率是影响烟支内流体传输性能的重要参数，本文提出了基于图像分割的一系列处理方法来提取烟支的孔隙率。

在燃烧过程中，烟气从燃烧部位透过未燃烧烟丝段被人体吸入，所以烟支的孔隙率对研究烟气在烟支内部的传输具有重要意义。为了得到烟支断层图像的孔隙率，本文首先采用Otsu获得背景像素与烟叶像素的分割阈值，之后在图像边缘开始进行漫水填充，目的是分割出烟支区域和背景区域；获得烟支区域后进行两次图像腐蚀，目的是腐蚀掉最外面一层的包裹纸张，之后根据之前提取的Otsu分割阈值来判断烟支区域中对应的像素点是否为孔隙(小于该阈值的被认为是孔隙)，计算出被判定为孔隙的像素点数占烟支区域总像素点数的百分比，如公式(4)所示，即为该断层的孔隙率。

(4)

式中，p为孔隙率；N′为烟支区域孔隙部分像素点数；N为烟支区域所有像素点数。

为了得到某一段烟支的总体孔隙率，则需要统计烟支中所有断层的孔隙率，求和平均后即为该段烟支的总体孔隙率。孔隙率的计算流程如图7所示。

4 实验结果与分析

针对烟支的螺旋轨迹CT扫描几何参数如表1所示，其中探测器的像元尺寸为0.2mm。这里采用连续扫描模式，探测器帧频设置为2，整个扫描过程的时间为370s。图8(a)给出了Katsevich重建方法得到的重建断层图像。实验结果表明，本文针对烟支的螺旋CT扫描重建方法可以有效地扩展纵向成像视场，但也存在着噪声问题。

图8 降噪图像前后对比示意图(a)原始噪声图像;(b)BM3D降噪后的图像

如图8(a)所示，螺旋CT扫描重建后的原始图像受到的噪声影响较为严重。图8(b)为BM3D降噪后的结果，能够观察到肉眼可见的噪声已基本被滤除，同时烟丝的图像细节得到了很好地保留。为了更加客观地评价BM3D降噪方法对CT图像的降噪效果，本文选取了峰值信噪比(PSNR)、均方误差(MSE)、边缘保持系数(EPI)和结构相似度系数(SSIM)四种客观评价指标，定量分析结果如表2所示。

表2 BM3D降噪前后图像参数对比

由表2可知，本文采用的降噪算法在峰值信噪比、均方误差和结构相似度的表现较好，在有效滤除噪声的同时，能够很好地保留原图像中的边缘和纹理等细节信息，在视觉效果上也有着更好的表现。

为了进一步说明改进降噪模型的普适性，本文对烟支的600层CT图像进行了降噪处理，并将降噪结果进行了三维体数据重建，其重建效果如图9所示。可以看出，降噪前后的三维体数据模型清晰度对比明显。应用本文降噪方法能明显提升噪声烟支图像的质量。

图9 烟支CT降噪三维重建前后对比示意图 (a)原始CT图像三维重建;(b)降噪后CT图像三维重建

在获取了高质量CT图像的基础上，图10为使用Zhang快速并行细化方法提取得到的烟丝骨架，右侧为放大的骨架，可以观察到并不满足骨架提取单像素细线的要求，这说明了该方法提取骨架信息并不彻底。

图1 烟支螺旋CT扫描方案

图10 Zhang快速并行细化方法提取得到的烟丝骨架

使用前文所述细化算法对提取的烟丝骨架进行了彻底细化，彻底细化结果如图11(b)所示。经过与图11(a)提取的原始骨架对比不难发现，采用基于纹线追踪的彻底细化算法有效地删除了多余的骨架点，使最终提取的骨架满足了骨架提取的原则。

图11 彻底细化前后对比图(a)原始骨架局部放大图;(b)彻底细化后局部放大图

图12为截取的某段烟支的孔隙率变化趋势图，显示了该段烟支的孔隙率分布情况。其中该段烟支的平均孔隙率为53.89%，最大孔隙率为61.84%，对应断层位于烟支中部，最小孔隙率为41.46%，对应断层位于烟支前端。结合传统测量方法得到的结果，证明了本文方法的正确性。通过观察与统计，结果对应的烟丝断层图像，孔隙率的变化规律与断层图像中烟丝密度的变化规律是相符的，这也验证了本文统计结果的合理性。

图12 某段烟支孔隙率变化趋势图

5 结论

本文以X射线CT检测技术为依托，开展烟支微观数字化表征的研究工作，设计了卷烟类长物体的螺旋CT扫描重建方法，进一步对CT断层图像进行BM3D降噪处理，在降低噪声的同时较好地保留了烟丝的图像细节，获得了较为清晰的断层图像。以高质量烟支CT图像为基础，采用Zhang快速并行骨架提取方法对烟支骨架进行提取和细化，基于骨架提取的不彻底现象，进一步采用了基于纹线追踪的彻底细化方法，得到了脉络清晰的烟丝拓扑学结构。以图像处理技术手段为基础，通过阈值分割、漫水填充和图像腐蚀等方法，计算得到烟支的孔隙率特征量。应用本文方法提取的烟支微观结构信息，可以为优化烟支结构和控制产品质量等提供理论依据与技术支持。

参考文献(References)

[1]FowlesJ,BatesM.Thechemicalconstituentsincigarettesandcigarettesmoke:Prioritiesforharmreduction[M]//AReporttotheNewZealandMinistryofHealth.KenepuruScienceCentre,Prirua, 2000: 1-67.

[2]ThomasAdam,JohnMcAughey,ChristophMocker,etal.Influenceoffilterventilationonthechemicalcompositionofcigarettemainstreamsmoke[J].AnalyticaChimicaActa, 2010,(657) :36-44.

[3]WangJ,PengX,XieY,etal.Fastanalysisofselectedcompoundsininhaledandexhaledvaporphaseofcigarettesmoketoevaluatecomponentsretainedintheupperrespiratorytract[J].RapidCommunMassSpectrom, 2021, 35:e8996.

[4] 申晓锋，李华杰，李善莲.烟丝结构表征方法研究[J]. 中国烟草学报，2010，4(16)：20-25.

ShenXiaofeng,LiHuajie,LiShanlian.Studyoncharacterizationofcuttobaccoparticlesizedistribution[J].ActaTabacariaSinica, 2010, 4(16):20-25.(inChinese)

[5] 朱令宇，李永杰.滤棒重量、压降稳定性对卷烟重量、吸阻稳定性的影响[J]. 福建分析测试，2009，18(2)：83-85.

ZhuLingyu,LiYongjie.Effectsofthestabilityofweightandresistanceoffilterrodonthestabilityofcigaretteweightanddrawresistance[J].FujianAnalysis&Testing, 2009, 18(2):83-85.(inChinese)

[6]LetitiaSchoeman,PaulWilliams,AntonduPlessis,etal.X-raymicro-computedtomography(μCT)fornon-destructivecharacterisationoffoodmicrostructure[J].TrendsinFoodScience&Technology, 2016, 47.

[7]ZhaoBudi,WangJianfeng.3Dquantitativeshapeanalysisonform,roundnessandcompactnesswithμCT[J].PowderTechnology, 2016, 291: 262-275.

[8] 熊瑛，刘海强，杜本莉，等.微焦点CT在陶瓷基复合材料上的检测应用[J]. 航空制造技术，2018，61(19)：58-63.

XiongYing,LiuHaiqiang,DuBenli,etal.ApplicationofMicro-FocusCToninspectionofceramicmatrixcomposites[J].AeronauticalManufacturingTechnology, 2018, 61(19):58-63.(inChinese)

[9]LuH,LiX,HsiaoIT,etal.Analyticaltreatmentofnoiseforlow-doseCTprojectiondatabypenalizedweightedleast-squaresmoothingintheK-Ldomain[J].SPIEMedImag, 2002, 4682:146-152.

[10]DemirkayaK.Reductionofnoiseandimageartifactsincomputedtomographybynonlinearfiltrationoftheprojectionimages[J].SPIEMedImag, 2001, 4322:917-923.

[11]HsiehJ.AdaptivestreakartifactreductionincomputedtomographyresultingfromexcessiveX-rayphotonnoise[J].MedPhys, 1998, 25:2139-2147.

[12] 于川芳，王兵，申玉军，等.国内外卷烟烟支结构的分析对比[C]//全面建设小康社会:中国科技工作者的历史责任，中国科协2003年学术年会论文集(上)，2003.

YuChuanfang,WangBing,ShenYujun,etal.Analysisandcomparisonofcigarettestructureathomeandabroad[C]//BuildingaModeratelyProsperousSocietyinanAll-roundWay:TheHistoricalResponsibilityofChineseScienceandTechnologyWorkers:Proceedingsofthe2003AnnualConferenceoftheChinaAssociationforScienceandTechnology(Part1), 2003.(inChinese)

[13] 王亮，银董红，谢国勇，等.卷烟燃烧过程中烟丝形态结构变化的表征方法[J]. 烟草科技, 2018，51(9)：67-72.

WangLiang,YinDonghong,XieGuoyong,etal.Amethodforcharacterizingvariationsofmorphologicalstructureofcutfillerduringcigaretteburning[J].TobaccoScience&Technology, 2018, 51(9):67-72.(inChinese)

[14]YuH,WangG.StudiesonimplementationoftheKatsevichalgorithmforspiralcone-beamCT[J].JournalofX-rayScienceTechnology, 2004, 12(2): 97-116.

[15]DabovK,FoiA,KatkovnikV,etal.Imagedenoisingbysparse3-Dtransform-domaincollaborativefiltering[J].IEEETransactionsonImageProcessing, 2007, 16(8): 2080-2095.

[16]BuadesA,CollB,MorelJ-M.Anon-localalgorithmforimagedenoising[C]//Proceedingsofthe2005IEEEComputerSocietyConferenceonComputerVisionandPatternRecognition(CVPR′05).

[17] 李迎，段汕.基于数学形态学的二值图像细化算法研究[J].中南民族大学学报(自然科学版)，2005，24(4)：96-99.

LiYing,DuanShan.Studyofbinaryimagethinningbasedonmathematicalmorphology[J].JournalofSouth-CentralMinzuUniversity(NaturalScienceEdition), 2005, 24(4):96-99.(inChinese)

[18] 吴选忠.Zhang快速并行细化算法扩展[J].福建工程学院学报，2006，4(1)：89-92.

WuXuanzhong.ExtendingofZhang’sfastparallelthinningalgorithm[J].JournalofFujianUniversityofTechnology, 2006, 4(1):89-92.(inChinese)

[19] 王晶，李景萃.指纹图像快速细化的改进算法及应用[J]. 仪器仪表学报，2008，29(7)：1535-1539.

WangJing,LiJingcui.Improvedfastthinningalgorithmforfingerprintimageanditsapplication[J].ChineseJournalofScientificInstrument, 2008, 29(7):1535-1539.(inChinese)