高聚物黏结炸药CT图像像素级裂纹识别方法研究

张 韬，宗和厚，陈 华，戴 斌，刘本德，骆吉洲

（1.中国工程物理研究院化工材料研究所，四川 绵阳 621999；2.哈尔滨工业大学计算机科学与技术学院，黑龙江 哈尔滨 150001）

1 引言

高聚物黏结炸药（Polymer bonded explosive，PBX）是由炸药晶体、黏合剂等构成的多组份复合含能材料，在应力影响、存放环境的物理条件变化等多因素的相互作用下，易形成不同程度、不同形态、不同扩展方式的裂纹［1-3］。这些裂纹一方面会使PBX 力学性能劣化，并最终导致材料破坏；另一方面作为能量聚集区域可能成为“热点”直接影响PBX 的力学性能、安全性能和爆轰性能［4-5］。PBX 裂纹的形成原因、发展和演化，及如何抑制裂纹的产生和发展等关系到PBX 的有效性和可用性，也关系到PBX 的力学性能和爆轰性能。为了评估PBX 的结构强度，建立强度准则和失效判据、分析起裂诱因以及裂纹对爆轰及环境适应性的影响，需要对裂纹进行精确的表征和定量分析。当前，对于PBX 内部裂纹的检测主要是采用X 射线层析成像（Computed Tomography，CT）技术，该技术对 PBX 裂纹的检测分析能力不足，主要表现在空间与密度分辨水平低，尺寸定量精度低，检测图像中的裂纹分割、识别能力欠缺，所以对PBX 在早期产生的细小裂纹无法及时诊断，无法对裂纹的宽度、走向、类型、位置等进行定量表征，也无法在结构出现开裂后从裂纹分布特征上研究开裂原因，解决上述问题的主要途径是建立高灵敏度、高精度、高适应性的裂纹识别方法。

当前HMX 基PBX 的CT 图像暗视野性、低对比度性、高噪声性等特点给裂纹识别带来了一定技术难度。现有的图像处理算法中可用于裂纹识别的算法包括阈值分割算法、基于灰度值梯度变化的检测算法和基于相位一致性的检测算法。这些算法从不同角度处理裂纹识别问题。阈值分割算法［6］是一类计算简单、效果最佳、应用最广的图像分割算法，它利用待提取目标与背景的灰度值差异，通过设置阈值将像素灰度值分类进而将目标和背景分离；这类算法及其改进主要适用于背景灰度一致、图像明暗度均匀和对比度较高的图像［7］，对于 HMX 基 PBX 裂纹提取效果不佳，细微裂纹无法提取，且提取的裂纹通常会出现较多的断点，连续性差。边缘检测方法将裂纹边缘视为图像中灰度值有阶跃性变化的像素集合，利用Sobel 算子、Canny 算子、Prewiit 算子、Robert 算子、Log 算子、Laplacian 算子、Kirsch算子［8］等各种边缘检测算子对裂纹进行提取；这类方法仅适用于边缘信息保存较为完整的图像，易将边缘信息较强的背景或噪声也当作裂纹；直接运用这些边缘算法来检测PBX 裂纹的效果也不理想；例如，比较研究发现［8］在这几种常用边缘算子中，Sobel算子检测到的边缘过宽，Log算子抗噪声能力过差，Canny算子在像素级边缘定位的精确性和抗噪声干扰性具有相对最好的效果。相位一致性方法通过定义相位一致性指标，将一致性指标超过阈值的像素视为边缘或裂纹像素［9-10］；这类方法能有效消除暗视野或低对比度对边缘检测的影响，从图像中提取到含角、线、纹理等信息的边缘，但其计算量大且对噪声异常敏感；由于HMX基PBX裂纹图像存在高噪声性，运用相位一致性方法检测裂纹效果很差，严重时甚至导致裂纹完全淹没在噪声中无法区分。针对上述问题，本研究结合统计学方法和连点成线的思想提出从PBX的CT图像提取裂纹的一种新方法——MDCE算法。该方法首先根据图像灰度值提取裂纹特征点，然后在裂纹骨架上利用马氏距离（Mahalanobis Distance）［11］判定裂纹扩展方向并删除不显著的特征点，最后从原始图像中拷贝像素灰度值来将骨架连接成完整裂纹。同时，将MDCE算法与Canny算法和相位一致性方法进行对比，MDCE算法可以提高裂纹识别精准度和改进裂纹表征方法。

2 裂纹提取算法

2.1 算法基本思想

本研究提出的MDCE 算法将CT 图像进行灰度转换、图像增强和噪声过滤等预处理操作后的灰度图像作为输入，其输出是图像中裂纹的识别结果。实验研究发现，在本研究的CT 图像上用常规方法无法识别到的裂纹表现出以下数字特征：

（1）图像整体明暗变化不显著，噪声像素的灰度值与裂纹像素灰度值接近，难以用现有算法区分；

（2）图像中待提取裂纹的灰度值梯度、相位一致性等数字特征与背景的数字特征相比不显著；

（3）人眼能够较准确地识别裂纹，主要是裂纹的部分像素的灰度值明显低于背景像素的平均灰度值并表现出很强的方向性和连续性；然而，噪声像素点的灰度值虽然也显著低于背景像素的平均灰度值，但其方向性和连续性差。

基于这一观察，本研究提出了一种两阶段自动裂纹检测算法，利用特征像素点的方向性和连续性来检测裂纹。第一个阶段利用统计方法从输入CT 图像中提取灰度值显著低于平均灰度值的特征像素点，提取结果中既包含裂纹骨架像素点，也包含了灰度值较小的噪声像素点。第二个阶段检测特征像素点的局部方向性，将局部方向性弱的特征像素点视为噪声剪除，而将局部方向性强的特征像素点用来确定裂纹的位置、方向和粗细，并通过从输入图像中拷贝特征像素点周围的像素灰度值来连接特征像素点得到完整的裂纹。

2.2 预处理

预处理是保障裂纹识别算法效果的重要步骤。主要包括灰度转换、图像增强和噪声过滤三个流程。灰度转换用于去除CT 图像中冗余的信息，图像增强旨在增强裂纹对比度，噪声过滤主要是降低噪声对裂纹检测的影响。在本文中，采用将任意像素（x，y）的三通道信息R（x，y）、G（x，y）、B（x，y）进行加权平均得到一个［0，255］内的灰度值f（x，y）的方法，转换过程为f（x，y）=0.299R（x，y）+0.587G（x，y）+0.114B（x，y）。采用先执行拉普拉斯锐化方法再执行Gamma 校正方法来进行图像增强［12］，最后选择非局部均值滤波［13］的方法过滤噪声平滑图像，较好保留裂纹检测所需数字特征。

2.3 特征像素点检测算法

特征像素点检测算法的输入是待检测的灰度图像，其输出是去除了非特征点像素影响并进行增强的图像，输入和输出图像具有相同大小。算法的基本思想是将输入CT 图像视为随机变量X∈［0，255］在图像像素空间上的分布，先统计得到X的均值E［X］和标准差σX，再利用Chebyshev 不等式（式（1））：

定义三个阈值Ti=E［X］-ti∙σX（i=1，2，3，t1＞t2＞t3＞1）将输入图像中所有像素根据灰度值分为四类，灰度值低于T1的像素被判定为特征像素点，灰度值介于T1和T2间的像素点被判定为强潜在特征像素点，灰度值介于T2和T3间的像素点被判定为弱潜在特征像素点，灰度值大于T3的像素点被判定为非特征像素点。最后，为了排除非特征像素点的影响并便于局部方向性的检测，将前三类像素点的灰度值拉伸到［0，255］。具体做法是，特征像素点的灰度值被线性拉伸到［0，80］，强潜在特征像素点的灰度值被线性拉伸到［80，120］，弱潜在特征像素点的灰度值被线性拉伸到［120，255］，非特征像素点被丢弃。

下面讨论如何确定阈值T1，T2，T3。一方面，由公式（1）和对称性，可以得出式（2）近似估计式：

另一方面，根据统计规律可知，随机变量X的取值大都在E［X］附近；并且，由于裂纹必然是稀疏、细密的，故CT 图像中裂纹的特征像素点占整幅CT 图像的比例通常不会太大（不超过10%）。据此，定义t1=81/2，t2=41/2，t3=21/2。这样，特征像素点占整幅图像像素点的比例不超过6.25%，特征像素点和强潜在特征像素点的总比例不超过12.5%，三类被保留的像素点的总比例不超过25%。

为改善图像加强的局部性，将直方图的概念运用到每个像素周围，建立层级的组合局域直方图来反映一个像素是在裂纹内部、边界上或远离裂纹。在对比度不显著的情况下，将这种信息融入算法，可以更显著地增强裂纹特征像素，弱化远离裂纹的像素，达到增强对比度的目的，提取数字特征不显著的干纹。

特征像素点检测算法具有很强的自适应性，尽管t1，t2，t3是根据统计规律和 PBX 的 CT 图像裂纹稀密性得出的。事实上，三类特征像素点的判定不仅依赖于t1，t2，t3的取值，还依赖于 CT 图像所有像素点的灰度均值和方差。因此，无论CT 图像整体偏暗或偏亮，特征像素点检测算法均能有效地检测出“等比例地”特征像素点，结果如图1 所示。

图1 骨架提取算法结果图Fig.1 Result of Skeleton extraction algorithm

从图1 可以看出：（1）在裂纹与背景对比度显著增强的区域，增强后裂纹能够被识别出来，达到了预期设计效果；（2）通过层级的组合局域直方图，在对比度不显著的部分，对裂纹特征点进行显著加强，最终提取出数字特征不显著的干纹。

2.4 特征像素点剪除和连接算法

特征像素点剪除和连接算法的输入是两幅图像，一幅是PBX 的CT 图像经灰度转换、图像增强和噪声过滤后得到的图像Imgh，另一幅是特征点检测算法输出的图像Imgc；特征像素点剪除和连接算法的输出是裂纹图像Imgo。特征像素点剪除和连接算法的核心任务是检测Imgc中特征点的局部方向性，对于局部方向性强的特征点，定位裂纹位置、粗细和方向，然后相应地从Imgh中拷贝对应像素的灰度值到Imgo中，最终产生裂纹检测结果。

首先讨论如何刻画像素集合Set 的局部方向性。图像中的任意像素（x，y）可以视为平面上的一个点，因此像素集合可以等效于一个平面点集Set。本研究借助马氏距离（Mahalanobis Distance）［11］来评估平面点集的局部方向性。直观上，这相当于在Set 的统计学中心（x0，y0）放置一个面积最小的椭圆来覆盖Set；如果椭圆的长轴明显长于短轴，则表明Set 具有强方向性；否则，Set 具有弱方向性。像素点集Set 的局部方向性可以如下计算。首先，统计Set 中所有点的横、纵坐标均值x0和y0；再统计横、纵坐标的方差VarX，VarY和协方差VarXY；然后，构造协方差矩阵并计算其逆矩阵Σ-1（由Σ 的对称性可知Σ-1也是对称矩阵，进而Σ-1的两个特征值必然介于0 和1 之间）；然后，计算得到 Σ-1的两个特征值λ1，λ2（1≥λ1≥λ2＞0）和对应的单位特征向量v1，v2。此时，v1对应椭圆长轴的方向而v2对应椭圆短轴的方向。本文定义比值λ1/λ2为 Set 的局部方向性度量。当λ1/λ2≥1.5 时，判定 Set 的局部方向性较好；否则，判定Set 的局部方向性差。此外，Set 中任意一点（x，y）到统计中心（x0，y0）的马氏距离定义为：

利用（3）式计算得到Set 中各点的马氏距离之后，将所有点的马氏距离从小到大排序，取0.6 分位点对应的距离值称为Set 的拷贝距离Dset。

下面运用局部方向性剪除和连接特征点像素得到完整的裂纹。首先，将Imgc划分为大小适中的网格G，第i行第j列网格内的特征像素点记为集合G（i，j），然后运用马氏距离方法计算每个G（i，j）的局部方向性并将所有G（i，j）按局部方向性度量递减排序并依次进行后续步骤的处理。

（1）其中局部方向性度量小于1.5 的部分，表示覆盖该网格中所有点的椭圆近似一个标准圆，说明该网格内的点排列是散乱，可以认定网格内的点为噪声点，将这类特征点剪除掉；

（2）针对局部方向性度量大于1.5 的每个网格G（i，j），将其与周围 8 个方向上的网格，如图2 所示（G（i-1，j-1）、G（i-1，j）、G（i-1，j+1）、G（i，j-1）、G（i，j+1）、G（i+1，j-1）、G（i+1，j）、G（i+1，j+1））依次合并形成新的集合，然后计算新的集合的方向性：

a.若这8 个新的集合表现的方向性都很弱，表示该网格中可能有两种情况：一是存在一个极小的裂纹，二是由于噪声引起的计算错误。针对这种情况，由于裂纹识别任务中要求可以假阳性，不能出现假阴性的结果，所以对于这种情况本文不做处理，将其保留下来；

b.当8 个新集合中有一个或多个方向表现出较强的局部方向性，则可认为中心集合G（i，j）与表现较强方向性的对应的邻居集合G（p，q）（p∈［i-1，i+1］q∈［j-1，j+1］）具有相似方向的裂纹，对G（i，j）与G（p，q）中的每一点都计算出其与本身集合中心点的马氏距离m1与合并集合的中心点的马氏距离m2。若某点的m1小于本身集合的拷贝距离D1并且m2小于合并后集合的拷贝距离D2，表明该点是裂纹上的点，将Imgh中对应位置的点的像素值拷贝进Imgo，这样就能将特征点连接成完整裂纹；

图2 基于网格八方向的特征点剪除和连接方法示意图Fig.2 Illustration of Eight-orientation-grid-based Pruning and Connecting method of feature points

3 热力加载试验与CT检测

3.1 试验样品

本试验的检测样品为HMX 基PBX，并将PBX 制备成Φ10 mm×10 mm 的圆柱形试件。

3.2 热力加载压缩试验

热力加载压缩试验在Instro 材料试验机上完成，将检测样品放置于材料试验机上进行压缩试验，首先对样品施加5 MPa 的载荷，然后逐渐增大直到样品产生裂纹，试验温度为50 ℃。力学试验结束后将样品从材料试验机取下，采用CT 进行扫描。

3.3 CT检测

CT 扫描成像实验在nanoVoxel 设备（Sanying Precision Engineering Research Center，Tianjin，China）上完成，试验所采用的参数为电压100 kV，电流100 μA，单幅投影图采样曝光时间0.3 ms，样品旋转360°，投影图采样数为1080，每幅投影图为同一角度的4 幅投影图叠加平均得到。本次试验中所得到的成像的空间分辨率为25 μm/pixel。

图3 CT 原理示意图Fig.3 CT principle diagram

4 结果与分析

从CT 检测结果可知经过热力耦合加载的PBX 损伤破坏模式是以微裂纹为主的脆性损伤，裂纹形态主要特征表现为微细网状、低对比度、干纹分支等。

图4 中展示了HMX 基高聚物粘结炸药的CT 图像，图中1 处的是主干裂纹，主干裂纹图像数字特征明显，有较高的对比度，通常算法即能检测出；2 处的是网状微细裂纹，该裂纹宽度通常低于三个像素，且对比度较低，从数字特征来看，该裂纹很多地方的数字特征和背景噪声一致，难以区分；3 处是低对比度裂纹，这类裂纹数字特征与背景噪声相近，对比度很低，通常像素灰度值相差10 以内，很能提取出来；4 处标识的是主干裂纹分支，由于突然从特征明显的主干裂纹过渡到特征模糊的分支裂纹，所以图像处理算法在该处提取的裂纹容易出现断裂。

分别采用Canny 边缘检测方法、基于相位一致性的方法与本研究MDCE 算法对上述三种裂纹进行识别处理，结果如图5～图7 所示。

图4 CT 检测结果Fig.4 CT results

图5 三种算法对微细网状裂纹的识别结果Fig.5 Network fine crack recognition results from three algorithm

图6 三种算法对低对比度裂纹的识别结果Fig.6 Low-contrast crack recognition results from three algorithms

从图5 可以看出，在微细网状裂纹的检测结果中，Canny 算法（图5b）只检测出部分干纹，对于干纹周围的微细裂纹根本检测不出，基于相位一致性的算法（图5c）能检测出干纹和微细裂纹，但是背景噪声很多，检测结果被淹没在噪声之中，相对而言，MDCE 算法（图5d）检测结果的精度最高，能体现出原图中微细网状裂纹的空间位置及形态，且对干纹部分的检测结果也是最连续的。

从图6 可以看出，在低对比度裂纹的检测结果中，基于相位一致性的算法能识别主要干纹，但是低对比度的裂纹完全湮没在噪声中，难以分辨。本研究MDCE 算法和Canny 边缘算子都能较好的识别出低对比度的裂纹。

图7 显示，三种算法都能检测出斜向分支的裂纹，但Canny 算法在裂纹分支处不连续，基于相位一致性的算法的背景噪声难以去除，将裂纹特征掩盖。本文算法检测出的斜向分支裂纹准确、清晰、连续。

图7 斜向分支裂纹识别结果Fig.7 Oblique branch crack recognition results

5 结论

（1）提出了PBX 炸药CT 裂纹图像上像素级裂纹提取的新方法——MDCE 算法。该方法针对CT 图像特征，先根据灰度值直接提取裂纹特征点，然后在裂纹特征骨架上利用马氏距离判定裂纹扩展方向，避免噪声对裂纹提取的影响，并从原始图像中拷贝像素灰度值来将骨架连接成完整裂纹，最终识别出清晰、完整、连续的裂纹。

（2）在采用三种方法对微细网状、低对比度、干纹分支三种不同形态的HMX 基PBX 裂纹进行识别，结果表明：对于边缘检测和相位一致性方法等常规方法无法识别的一些细微裂纹，本方法能较准确和有效地识别。

（3）在MDCE 算法的基础上进一步研究，有望改善裂纹描述方法，提高裂纹表征能力，分析裂纹空间分布差异对PBX 起爆性能的影响，探索不同裂纹形态对PBX 起爆性能的影响规律，为外界力热作用下PBX 起裂诱因以及裂纹对环境适应性的影响评估提供新的研究基础。