CN114782354A - 爆炸动焦散斑特征识别提取及应力强度因子自适应算法 - Google Patents

Abstract

本发明公开了爆炸动焦散斑特征识别提取及应力强度因子自适应算法，具体包括以下步骤：S1:设计装药结构、制作药包并进行炮孔装药固定；S2：开展爆炸加载动态焦散线实验并拍摄记录整个过程；S3：批量读取测试得到的待分析爆炸加载实验焦散斑图像；S4：对批量待分析焦散斑图像进行Mask匀光处理；S5：将图像转换为灰度图并进行非局部均值滤波。本发明公开的爆炸动焦散斑特征识别提取及应力强度因子自适应算法具有实现了动态数字焦散实验焦散斑特征自动识别提取和应力强度因子等参数自适应计算、存储和输出过程，克服了传统分析精度差，分析工作量大，分析过程繁琐不具有自适应性且受分析人员主观意志影响严重的缺陷效果。

Description

爆炸动焦散斑特征识别提取及应力强度因子自适应算法

技术领域

本发明涉及动态断裂力学实验测量方法领域，尤其涉及爆炸动焦散斑特征识别提取及应力强度因子自适应算法。

背景技术

动态焦散线实验方法是利用纯几何光学的映射关系，将材料裂纹尖端奇异区域的复杂受力变形问题转化为简单且清晰明确的焦散斑图像，只需通过测量焦散斑的特征参数，便可得到裂纹尖端的应力强度因子等重要指标。因此，该方法已成为研究类岩石材料在爆炸荷载作用下裂纹尖端应力奇异场问题的重要手段。

焦散线实验过程中通过高速相机采用较小的像素场，匹配较高的拍摄速度来捕捉更加详细的爆炸裂纹的动态扩展过程。现阶段焦散斑相关特征及参数的确定通常是通过将高速相机拍摄到的待分析照片通过图像处理软件进行人工标定测量，数据测量工作量较大，测量误差受图像噪声及测量人员的主观意识影响较大，准确性、重复性及分析时效性均较差，直接影响到相关分析结论的科学性和可靠度。因此，亟需一种便捷高效的爆炸荷载下动态焦散线实验焦散斑识别和特征参数提取并实现批量图像数据分析处理的方法。

发明内容

本发明公开爆炸动焦散斑特征识别提取及应力强度因子自适应算法，旨在解决背景技术中提到的现阶段焦散斑相关特征及参数的确定通常是通过将高速相机拍摄到的待分析照片通过图像处理软件进行人工标定测量，数据测量工作量较大，测量误差受图像噪声及测量人员的主观意识影响较大，准确性、重复性及分析时效性均较差，直接影响到相关分析结论的科学性和可靠度的技术问题。

为了实现上述目的，本发明采用了如下技术方案：

爆炸动焦散斑特征识别提取及应力强度因子自适应算法，具体包括以下步骤：

S1:设计装药结构、制作药包并进行炮孔装药固定；

S2：开展爆炸加载动态焦散线实验并拍摄记录整个过程；

S3：批量读取测试得到的待分析爆炸加载实验焦散斑图像；

S4：对批量待分析焦散斑图像进行Mask匀光处理；

S5：将图像转换为灰度图并进行非局部均值滤波；

S6：对图像采用canny算法进行焦散斑边缘检测；

S7：采用Hough变换进行焦散斑识别并定位所识别圆形轮廓的中心坐标；

S8：将满足直径最大和最小圆方程等式的中心点坐标进行标识并存储；

S9：根据坐标定位在图像中绘制焦散斑最大、最小直径圆周轮廓；

S10：根据相关公式自动计算所提取焦散斑的直径、面积、周长等特征参数并存储输出；

S11：计算获取裂纹扩展速度、加速度、裂纹尖端应力强度因子数值随时间变化曲线并存储输出；

S12：计算确定并存储输出裂纹动态能量释放率随时间变化曲线。

在一个优选的方案中，所述S2步骤中，动态焦散线实验的仪器主要包括激光光源、扩束镜、场镜即凸透镜、炸药起爆装置、高速相机、数据采集计算机；所述激光光源采用绿色激光，且光源的波长为532mm，所述激光光源的激光功率在0～200mW之间可调，所述扩束镜将激光斑点发散，扩大光场，两个凸透镜作为场镜，直径为30cm，焦距为90cm；第一透镜将发射光场调整为平行光场，第二透镜将平行光场汇聚到相机镜头处；试件中炸药起爆采用MD-2000多通道脉冲点火器。高速相机采用Fastcam-SA5-16G相机。视场像素尺寸为320×192像素，实际拍摄区域约为6cm×10cm。相机参数设置要保证拍摄到的图片具有足够的视场记录爆生裂纹的扩展轨迹，同时图片亮度和分辨率满足焦散斑的测量需求。

在一个优选的方案中，所述S3步骤中，在所述S2中拍摄到的照片内的裂纹尖端位置按照一定的比例与实际的试件进行换算，通过程序自动测量每一时刻焦散斑的近似圆心，采用数值微分中心点公式可得每一时刻裂纹扩展速度，设测量的k组焦散斑近似圆心(xi,yi),每两幅焦散斑时间间隔为△t；

根据数值微分中心点公式计算：

其中：当i＝1时，

当i＝k时，

△t—相邻两张照片拍摄时间间隔。

在一个优选的方案中，所述Mask匀光处理的Mask匀光算法数学模型可表示为：

f(x,y)＝r(x,y)+δ(x,y)；

式中，f(x,y)代表分析图像；δ(x,y)代表背景影像；r(x,y)代表亮度分布均匀的影像；

所述Mask匀光处理的具体过程如下：

S41.背景影像的获取:背景影像的获取是MASK匀光处理最为关键的一步，其质量的好坏直接影响了匀光效果的好坏,考虑到影像亮度和反差分布的不均匀性，很难用一个数学模型来处理，通常使用低通滤波来获得可以反映影像颜色和亮度变化的低频信息，作为背景影像；

S42.影像相减运算：将原始影像与背景影像做相减运算，得到剔除了低频信息的影像，为了保持匀光前后影像整体亮度基本一致，相减运算的同时还要加上代表原始影像亮度均值的偏移量:；

S43.反差拉伸处理：由于相减运算会导致图像灰度值的动态范围缩小，影像反差减小，对上述结果进行拉伸处理，增加影像整体反差，得到最终的匀光结果，进行拉伸处理的数学模型选取如下:

式中，f(x,y)表示最终的结果图像像，f'(x,y)表示相减运算后的图像，m为拉伸前图像灰度值均值，β为拉伸因子，

为拉伸系数，具体公式如下：

式中，

表示(x，y)像素点的拉伸系数，α为调节参数，f_bk(x，y)表示背景影像中(x，y)点的灰度值，f_bkmax表示背景影像的最大灰度值，f_bkmin表示背景影像的最小灰度值；建立信息熵、清晰度、平均对比度和块标准差四个评价指标并设定相应的阈值，当同时满足上述四项指标设定阈值便认为匀光效果理想，否则调整相关参数进行重新计算直至达到效果要求为止；

所述S5的步骤中的非局部均值滤波算法(NLM)是一种较为优秀的降噪算法，该方法优势在于可以在整个图像像素空间内对相似像素进行加权平均，而不是局限于某个邻域，不仅可以有效降噪，而且可以保留细节，锐化边缘。其算法原理如下：

假设含噪图像为v＝{v(a)|a∈A}，去噪之后的图像为NL[v],每个像素a的加权值可用下式求出：

其中，w(a,b)表征像素a与像素b之间的相似度(高斯加权欧式距离)权值，满足0≤w(a,b)≤1且∑_bw(a,b)≤1，其表达式如下：

其中，

表示以像素a为中心的子块内像素和以像素b为中心的子块内像素的高斯加权距离的平方，用于衡量两个像素点之间的相似度，v(N(a))表示a周围的局部子块像素集合，h为滤波参数。

在一个优选的方案中，所述非局部均值滤波主要步骤如下：

S11.在含噪图像中建立一个边长为D的搜索块，再建立一个边长为d的相似块；

S12.在搜索块中选择一个像素点i，以i为中心构建一个相似块；在搜索块中任取一点j构建同样大小的相似块；依据高斯加权欧式距离计算权重系数w(i,j)；

S13.在搜索窗中遍历所有点j；

S14.将权重的最大值赋予点i，并将权重归一化，通过加权求和公式计算去噪后第i点的像素值：

其中，NL[v](i)为采用非局部均值算法得到的新的像素灰度值；w(i,j)为权重系数，v(j)为像素点j的灰度值,w(i,j)的具体表达式为：

式中，v(N_i)和v(N_j)表示以像素点i和像素点j为中心构建的相似块；e指数部分代表高斯加权欧式距离；a代表构建高斯核时的加权系数；h是构建权重的参数，控制衰减速度，h越小边缘特点保持得越好，但噪声严重，反之图像光滑但边缘保持的较差。

S15.最后遍历图像中所有的i点，完成去噪过程并通过PSNR检验图像去噪效果。

在一个优选的方案中，所述S6的步骤中，离散Canny算子边缘检测提取步骤如下：

S61.采用尺度为W的平滑算子s按行对图像做卷积运算，然后用相应的最优求导算子d对平滑后的算子按列做卷积运算，得到Gy；再用s按列对图像做卷积运算，然后d对平滑后的算子按行做卷积运算，得到Gx；

S62.计算图像水平和垂直的梯度幅值和方向。该方法中使用一阶偏导的有限差分来计算梯度幅值和方向。对于函数图像梯度大小和计算梯度方向计算如下：

方向梯度计算：

梯度幅值计算：

梯度方向计算：

S63.非极大值抑制：遍历G，把边缘的梯度方向大致分为水平、竖直、45°和135°四个方向，各个方向用不同的邻近像素进行比较，以决定局部极大值，若某个像素的灰度值与其梯度方向上前后两个像素的灰度值相比不是很大的，则将该像素值置为0，即非边缘，该过程称为非极大值抑制；

S64.双阈值操作：设置高低阈值，遍历G，凡是大于高阈值的设为边缘；凡是小于低阈值的设为非边缘，若检测结果介于高、低阈值之间，则要判断像素的邻接像素中是否存在超过高阈值的边缘像素，若有则设为边缘，否则设为非边缘；

S65：通过抑制孤立的弱边缘最终完成边缘检测。

在一个优选的方案中，所述S7的步骤中，Hough圆检测算法利用圆形本身的几何特性，将二维像素点映射到高维的Hough空间，以一定的半径步长对得到的圆心和半径做检测分析，同时，引入聚类分析的思想，合并相似度高的圆形，有效减少了提取结果中由于边缘模糊而造成的多个假圆的问题，提高了准确性并减少了计算工作量。

在一个优选的方案中，所述Hough圆检测算法的具体步骤为：

S71.以图像左下角为坐标原点建立二维平面坐标系。对图像进行预处理，经过高斯去噪、二值化、边缘提取等步骤，将所建立坐标系下图像的边缘像素信息保存在集合S中；

S72.在集合S中随机地抽取一个像素点X，构造三个直角三角形，分别沿着水平方向和垂直方向构造之间三角形A1，沿着-45°和45°构造直角三角形A2，沿着-15°和75°构造直角三角形A3；

S73.A1,A2,A3的直径相交于一点B，以B为中心，得到X的对称点Y；

S74.如果A1,A2,A3的斜边和X与Y组成的直径满足圆形的几何性质，则将此圆的圆心坐标和半径r、直径D信息加入到集合C中；根据Hough变换识别提取到的圆心和半径，编制程序计算所识别标记的焦散斑的周长、面积并存储输出；具体公式为：

周长：L＝2πr＝πD

面积：

通过编制循环程序计算每一时刻焦散斑的直径D_max和D_min，可得裂纹尖端动态应力强度因子，动焦散试验方法中动态应力强度因子的值与焦散斑直径按照如下公式计算：

式中：

—I型裂纹尖端应力强度因子；

—Ⅱ型裂纹尖端应力强度因子；

D_max—焦散斑最大直径，由前述程序自动计算得到；

g—应力强度因子的数值系数，根据g～μ曲线确定；

F(v)—裂纹扩展速度的调节系数，通常在具有实际意义的裂纹速度下取值为1；

d_eff—试件有效厚度即材料实际厚度，根据试验材料实际厚度确定；

c—材料的应力光学常数，取为0.85×10-10m2/N；

z₀—参考平面到试件表面的距离，根据实验布置确定，通常取为100mm；

μ—Ⅰ、Ⅱ型动态应力强度因子之比，通过曲线得出。

在一个优选的方案中，所述动态能量释放率作为裂纹扩展过程中的能量量度指标，是控制扩展裂纹尖端的应力场、应变场的参量，是裂纹扩展过程中弹性能量释放，反映了裂纹扩展的驱动力；平面状态下能量释放率计算公式如下：

式中：

α_d＝1-(v/c₁)²

α_s＝1-(v/c₂)²

其中，c₁，c₂分别为膨胀波波速和剪切波波速，v为裂纹扩展速度，可利用裂纹尖端强度因子计算得到裂纹尖端的动态能量释放率，此外亦可根据分析需要将相关公式编制程序，得到扩展裂纹尖端动态能量释放率与裂纹长。

由上可知，爆炸动焦散斑特征识别提取及应力强度因子自适应算法，具体包括以下步骤：

S1:设计装药结构、制作药包并进行炮孔装药固定；

S2：开展爆炸加载动态焦散线实验并拍摄记录整个过程；

S3：批量读取测试得到的待分析爆炸加载实验焦散斑图像；

S4：对批量待分析焦散斑图像进行Mask匀光处理；

S5：将图像转换为灰度图并进行非局部均值滤波；

S6：对图像采用canny算法进行焦散斑边缘检测；

S7：采用Hough变换进行焦散斑识别并定位所识别圆形轮廓的中心坐标；

S8：将满足直径最大和最小圆方程等式的中心点坐标进行标识并存储；

S9：根据坐标定位在图像中绘制焦散斑最大、最小直径圆周轮廓；

S10：根据相关公式自动计算所提取焦散斑的直径、面积、周长等特征参数并存储输出；

S11：计算获取裂纹扩展速度、加速度、裂纹尖端应力强度因子数值随时间变化曲线并存储输出；

S12：计算确定并存储输出裂纹动态能量释放率随时间变化曲线。本发明提供的爆炸动焦散斑特征识别提取及应力强度因子自适应算法具体以下技术效果：

1、实现了动态数字焦散实验焦散斑特征自动识别提取和应力强度因子等参数自适应计算、存储和输出过程，克服了传统分析精度差，分析工作量大，分析过程繁琐不具有自适应性且受分析人员主观意志影响严重的缺陷；

2、通过组合算法解决了动态焦散线实验批量图像数据处理分析难题，改善了传统分析方法易受炮烟和图像背景噪声影响的缺陷，增强了数据处理和结果分析的时效性；

3、提供了一种适应未来自动化、智能化和精细化的焦散斑识别和特征提取算法，进一步提供了更为全面、科学的动态焦散线实验过程数字化技术方案，解决了传统人工测量方式测量精度差、可操作性差、受测量人员主观意识影响程度高的缺陷，为高效、安全、精细化爆破设计和参数优化选取提供了智能化、便捷高效的特征识别和提取算法。

附图说明

图1为本发明提出的爆炸动焦散斑特征识别提取及应力强度因子自适应算法的流程图。

图2为本发明提出的爆炸动焦散斑特征识别提取及应力强度因子自适应算法的Ⅰ型和Ⅰ-Ⅱ复合型裂纹焦散斑形状示意图。

图3为本发明提出的爆炸动焦散斑特征识别提取及应力强度因子自适应算法的数字激光动态焦散系统。

图4为本发明提出的爆炸动焦散斑特征识别提取及应力强度因子自适应算法的数字激光动态焦散线原理图。

图5为本发明提出的爆炸动焦散斑特征识别提取及应力强度因子自适应算法的爆炸实验试件及固定试件的加载固定装置结构示意图。

图6为本发明提出的爆炸动焦散斑特征识别提取及应力强度因子自适应算法的非局部均值滤波边长为D的搜索块示意。

图7为本发明提出的爆炸动焦散斑特征识别提取及应力强度因子自适应算法的Hough变换圆识别算法流程。

图8为本发明提出的爆炸动焦散斑特征识别提取及应力强度因子自适应算法的I型裂纹原始测试焦散图像。

图9为本发明提出的爆炸动焦散斑特征识别提取及应力强度因子自适应算法的Ι型裂纹焦散图像Mask匀光处理提取的背景图像。

图10为本发明提出的爆炸动焦散斑特征识别提取及应力强度因子自适应算法的I型裂纹焦散图像Mask匀光后灰度处理图像。

图11为本发明提出的爆炸动焦散斑特征识别提取及应力强度因子自适应算法的I型裂纹焦散斑离散canny算法边缘检测结果。

图12为本发明提出的爆炸动焦散斑特征识别提取及应力强度因子自适应算法的Hough变换I型裂纹焦散斑圆心及直径标定及存储输出。

图13为本发明提出的爆炸动焦散斑特征识别提取及应力强度因子自适应算法的Ⅰ-Ⅱ复合型裂纹原始测试焦散图像。

图14为本发明提出的爆炸动焦散斑特征识别提取及应力强度因子自适应算法的Ⅰ-Ⅱ复合型裂纹Mask匀光处理提取的背景图像。

图15为本发明提出的爆炸动焦散斑特征识别提取及应力强度因子自适应算法的Ⅰ-Ⅱ复合型焦散图像Mask匀光后灰度处理图像。

图16为本发明提出的爆炸动焦散斑特征识别提取及应力强度因子自适应算法的Ⅰ-Ⅱ复合型裂纹焦散斑离散canny算法边缘检测结。

图17为本发明提出的爆炸动焦散斑特征识别提取及应力强度因子自适应算法的Hough变换Ⅰ-Ⅱ复合型裂纹焦散斑圆心及直径标定及存储输出。

图18为本发明提出的爆炸动焦散斑特征识别提取及应力强度因子自适应算法的裂纹尖端扩展速度随时间变换曲线。

图19为本发明提出的爆炸动焦散斑特征识别提取及应力强度因子自适应算法的裂纹尖端扩展速度随时间变换曲线。

图20为本发明提出的爆炸动焦散斑特征识别提取及应力强度因子自适应算法的裂纹尖端动态应力强度因子随时间变化曲线。

图21为本发明提出的爆炸动焦散斑特征识别提取及应力强度因子自适应算法的动态能量释放率随时间变换曲线。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。

动态焦散线方法实验原理为：基于几何光学原理，将试件放置于两场镜间的平行光场，在爆炸应力波作用下，试件产生裂纹，在裂纹尖端存在高应力奇异场，裂纹尖端附近的试件厚度d、密度ρ、折射率n等均发生急剧变化，入射光通过试件后将发生偏转，从而在裂纹尖端产生一个光线暗区，称之为焦散斑。在焦散斑暗区与周围亮区之间存在一条亮线，称为焦散线。根据弹性断裂力学理论，各向同性材料平面试件在爆炸荷载下，爆生裂纹主要为张开型(Ⅰ型)和复合型(Ⅰ-Ⅱ型)其中，Ⅰ型、复合型裂纹焦散斑如图1所示。

所述数字激光爆炸加载动态焦散线实验系统，如图2所示。主要包括激光光源、扩束镜、场镜(凸透镜)、炸药起爆装置、高速相机、数据采集计算机。光源采用绿色激光，波长为532mm，激光功率在0～200mW之间可调。扩束镜将激光斑点发散，扩大光场。两个凸透镜作为场镜，直径为30cm，焦距为90cm。透镜1将发射光场调整为平行光场，透镜2将平行光场汇聚到相机镜头处。试件中炸药起爆采用MD-2000多通道脉冲点火器。高速相机采用Fastcam-SA5-16G相机。视场像素尺寸为320×192像素，实际拍摄区域约为6cm×10cm。相机参数设置要保证拍摄到的图片具有足够的视场记录爆生裂纹的扩展轨迹，同时图片亮度和分辨率满足焦散斑的测量需求。

当试件受到炸药产生的爆炸荷载作用时，折射率随裂纹尖端奇异点附近厚度的变化而发生改变。当平行光场经过折射率改变的试件后，反射光和折射光都将偏离平行状态，折射光便会在相距Z0的平面上形成如图3所示的光线集中和阴影区，该阴影区即为焦散斑，偏离的光线再经过透镜2汇聚，在其焦距处的高速相机中成像。

参照图1-21，爆炸动焦散斑特征识别提取及应力强度因子自适应算法，爆炸动焦散斑特征识别提取及应力强度因子自适应算法，具体包括以下步骤：

S1:设计装药结构、制作药包并进行炮孔装药固定；

S2：开展爆炸加载动态焦散线实验并拍摄记录整个过程；

S3：批量读取测试得到的待分析爆炸加载实验焦散斑图像；

S4：对批量待分析焦散斑图像进行Mask匀光处理；

S5：将图像转换为灰度图并进行非局部均值滤波；

S6：对图像采用canny算法进行焦散斑边缘检测；

S7：采用Hough变换进行焦散斑识别并定位所识别圆形轮廓的中心坐标；

S8：将满足直径最大和最小圆方程等式的中心点坐标进行标识并存储；

S9：根据坐标定位在图像中绘制焦散斑最大、最小直径圆周轮廓；

S10：根据相关公式自动计算所提取焦散斑的直径、面积、周长等特征参数并存储输出；

S11：计算获取裂纹扩展速度、加速度、裂纹尖端应力强度因子数值随时间变化曲线并存储输出；

S12：计算确定并存储输出裂纹动态能量释放率随时间变化曲线。

在一个优选的实施方式中，实验时将高速摄影仪的拍摄帧数设置为44000帧/秒，则拍摄到的相邻两张照片之间的间隔为22.7μs足够拍摄到裂纹的扩展轨迹和扩展过程。通过程度读取焦散斑图像并进行坐标原点标定，将拍摄到的照片内的裂纹尖端位置按照一定的比例与实际的试件进行换算，通过程序自动测量每一时刻焦散斑的近似圆心，采用数值微分中心点公式可得每一时刻裂纹扩展速度。设测量的k组焦散斑近似圆心(x_i,y_i),每两幅焦散斑时间间隔为△t。根据数值微分中心点公式计算：

其中：当i＝1时，

当i＝k时，

△t—相邻两张照片拍摄时间间隔。

在一个优选的实施方式中，。

参照图1和图3，在一个优选的实施方式中，所示PMMA试件的几何尺寸为200mm×200mm×5mm，在试件中心钻一个直径为10mm的圆孔作为炮孔，所述装药形式采用不耦合装药，不耦合系数为1.67，试验步骤为：

1.使用透明胶带粘在炮孔一侧并压入橡皮泥形成孔底，采用塑料吸管作为装药容器形成不耦合装药；

2.在孔内内装入100mg的Pb(N₃)₂之后插入自制的起爆探针，该起爆探针是由数码电子雷管元件端部焊接点火药头制成，使用橡皮泥对孔口进行封堵防止搬运过程中炸药洒落；

3.安装采用亚克力板制成的5mm厚度的防护板，防止在起爆时炸药、试件的碎片飞溅造成仪器、人员损伤；

4.在炮孔两侧采用夹具夹住绝缘橡胶皮将炮孔封堵，防止爆生气体逸散，夹具示意图如图4所示；

5.连接起爆探针上的起爆线到起爆器，通过监测线路、设置起爆时间、充电之后等待高速摄影仪实验员发出拍摄指令后松开起爆键，完成起爆，在起爆之前的瞬间高速摄影仪开始拍摄，精确捕捉裂纹尖端焦散斑形态特征，保证了拍摄起爆过程的完整性。由于装药结构具有对称性，为了满足视场需求，高速相机只记录一侧的焦散斑运动轨迹。

所述Mask匀光算法，其原理是选择一张称为Mask的密度或色彩随位置缓慢变化的模糊的透明薄膜作为正片，与原始图像(称为负片)按轮廓线叠加进行晒像获得一张色彩密度较为均匀的像片，随后进行晒印以增强像片总体反差。Mask匀光法的特点是消减影像中的大反差而增强影像中的小反差，从而使整张影像照度与色彩趋于一致。Mask匀光算法不仅可以提高影像的目视效果，还有利于影像的后续处理和应用。Mask匀光算法数学模型可表示为：

f(x,y)＝r(x,y)+δ(x,y)

式中，f(x,y)代表分析图像；δ(x,y)代表背景影像；r(x,y)代表亮度分布均匀的影像。

图像Mask匀光处理的具体过程如下：

1.背景影像的获取。背景影像的获取是MASK匀光处理最为关键的一步，其质量的好坏直接影响了匀光效果的好坏。考虑到影像亮度和反差分布的不均匀性，很难用一个数学模型来处理，通常使用低通滤波来获得可以反映影像颜色和亮度变化的低频信息，作为背景影像。

2.影像相减运算。将原始影像与背景影像做相减运算，得到剔除了低频信息的影像，为了保持匀光前后影像整体亮度基本一致，相减运算的同时还要加上代表原始影像亮度均值的偏移量。

3.反差拉伸处理。由于相减运算会导致图像灰度值的动态范围缩小，影像反差减小，对上述结果进行拉伸处理，增加影像整体反差，得到最终的匀光结果。进行拉伸处理的数学模型选取如下:

式中，f(x,y)表示最终的结果图像像，f'(x,y)表示相减运算后的图像，m为拉伸前图像灰度值均值，β为拉伸因子，

为拉伸系数，具体公式如下：

式中，

表示(x，y)像素点的拉伸系数，α为调节参数，f_bk(x，y)表示背景影像中(x，y)点的灰度值，f_bkmax表示背景影像的最大灰度值，f_bkmin表示背景影像的最小灰度值。

建立信息熵、清晰度、平均对比度和块标准差四个评价指标并设定相应的阈值，当同时满足上述四项指标设定阈值便认为匀光效果理想，否则调整相关参数进行重新计算直至达到效果要求为止。

非局部均值滤波算法(NLM)是一种较为优秀的降噪算法，该方法优势在于可以在整个图像像素空间内对相似像素进行加权平均，而不是局限于某个邻域，不仅可以有效降噪，而且可以保留细节，锐化边缘。其算法原理如下：

假设含噪图像为v＝{v(a)|a∈A}，去噪之后的图像为NL[v],每个像素a的加权值可用下式求出：

其中，w(a,b)表征像素a与像素b之间的相似度(高斯加权欧式距离)权值，满足0≤w(a,b)≤1且∑_bw(a,b)≤1，其表达式如下：

其中，

表示以像素a为中心的子块内像素和以像素b为中心的子块内像素的高斯加权距离的平方，用于衡量两个像素点之间的相似度，v(N(a))表示a周围的局部子块像素集合，h为滤波参数。

在一个优选的实施方式中，非局部均值滤波主要步骤如下：

1.在含噪图像中建立一个边长为D的搜索块，再建立一个边长为d的相似块；

2.在搜索块中选择一个像素点i，以i为中心构建一个相似块；在搜索块中任取一点j构建同样大小的相似块；依据高斯加权欧式距离计算权重系数w(i,j)；

3.在搜索窗中遍历所有点j；

4.将权重的最大值赋予点i，并将权重归一化，通过加权求和公式计算去噪后第i点的像素值：

其中，NL[v](i)为采用非局部均值算法得到的新的像素灰度值；w(i,j)为权重系数，v(j)为像素点j的灰度值。w(i,j)的具体表达式为：

式中，v(N_i)和v(N_j)表示以像素点i和像素点j为中心构建的相似块；e指数部分代表高斯加权欧式距离；a代表构建高斯核时的加权系数；h是构建权重的参数，控制衰减速度。h越小边缘特点保持得越好，但噪声严重，反之图像光滑但边缘保持的较差。

5.最后遍历图像中所有的i点，完成去噪过程并通过PSNR检验图像去噪效果。

在目前常用的边缘检测方法中，离散Canny算子边缘检测算法是具有严格定义的，可以提供良好可靠检测的方法之一。由于它具有满足边缘检测的三个标准和实现过程简单的优势，成为边缘检测最流行的算法之一。离散Canny算子边缘检测提取步骤如下：

1.采用尺度为W的平滑算子s按行对图像做卷积运算，然后用相应的最优求导算子d对平滑后的算子按列做卷积运算，得到G_y；再用s按列对图像做卷积运算，然后d对平滑后的算子按行做卷积运算，得到G_x；

2.计算图像水平和垂直的梯度幅值和方向。该方法中使用一阶偏导的有限差分来计算梯度幅值和方向。对于函数图像梯度大小和计算梯度方向计算如下：

方向梯度计算：

梯度幅值计算：

梯度方向计算：

3.非极大值抑制：遍历G，把边缘的梯度方向大致分为水平、竖直、45°和135°四个方向。各个方向用不同的邻近像素进行比较，以决定局部极大值。若某个像素的灰度值与其梯度方向上前后两个像素的灰度值相比不是很大的，则将该像素值置为0，即非边缘，该过程称为非极大值抑制；

4.双阈值操作：设置高低阈值，遍历G，凡是大于高阈值的设为边缘；凡是小于低阈值的设为非边缘。若检测结果介于高、低阈值之间，则要判断像素的邻接像素中是否存在超过高阈值的边缘像素，若有则设为边缘，否则设为非边缘；

5.通过抑制孤立的弱边缘最终完成边缘检测。

Hough圆检测算法利用圆形本身的几何特性，将二维像素点映射到高维的Hough空间，以一定的半径步长对得到的圆心和半径做检测分析。同时，引入聚类分析的思想，合并相似度高的圆形，有效减少了提取结果中由于边缘模糊而造成的多个假圆的问题，提高了准确性并减少了计算工作量。Hough圆检测算法的具体步骤为：

1.以图像左下角为坐标原点建立二维平面坐标系。对图像进行预处理，经过高斯去噪、二值化、边缘提取等步骤，将所建立坐标系下图像的边缘像素信息保存在集合S中；

2.在集合S中随机地抽取一个像素点X，构造三个直角三角形，分别沿着水平方向和垂直方向构造之间三角形A₁，沿着-45°和45°构造直角三角形A₂，沿着-15°和75°构造直角三角形A₃；

3.A₁,A₂,A₃的直径相交于一点B，以B为中心，得到X的对称点Y；

4.如果A₁,A₂,A₃的斜边和X与Y组成的直径满足圆形的几何性质，则将此圆的圆心坐标和半径r、直径D信息加入到集合C中。

根据Hough变换识别提取到的圆心和半径，编制程序计算所识别标记的焦散斑的周长、面积并存储输出。具体公式为：

周长：L＝2πr＝πD

面积：

通过编制循环程序计算每一时刻焦散斑的直径D_max和D_min，可得裂纹尖端动态应力强度因子。动焦散试验方法中动态应力强度因子的值与焦散斑直径按照如下公式计算：

式中：

—I型裂纹尖端应力强度因子；

—Ⅱ型裂纹尖端应力强度因子；

D_max—焦散斑最大直径，由前述程序自动计算得到；

g—应力强度因子的数值系数，根据g～μ曲线确定；

F(v)—裂纹扩展速度的调节系数，通常在具有实际意义的裂纹速度下取值为1；

d_eff—试件有效厚度即材料实际厚度，根据试验材料实际厚度确定；

c—材料的应力光学常数，取为0.85×10^-10m²/N；

z₀—参考平面到试件表面的距离，根据实验布置确定，通常取为100mm。

μ—Ⅰ、Ⅱ型动态应力强度因子之比，通过(D_max-D_min)/D_max～μ曲线得出。

动态能量释放率作为裂纹扩展过程中的能量量度指标，是控制扩展裂纹尖端的应力场、应变场的参量，是裂纹扩展过程中弹性能量释放，反映了裂纹扩展的驱动力。平面状态下能量释放率计算公式如下：

式中：

α_d＝1-(v/c₁)²

α_s＝1-(v/c₂)²

其中，c₁，c₂分别为膨胀波波速和剪切波波速，v为裂纹扩展速度。因此，可利用裂纹尖端强度因子计算得到裂纹尖端的动态能量释放率。

此外亦可根据分析需要将相关公式编制程序，得到扩展裂纹尖端动态能量释放率与裂纹长度、扩展速度、扩展加速度和动态应力强度因子之间的对应关系曲线。

以上所述，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，根据本发明的技术方案及其发明构思加以等同替换或改变，都应涵盖在本发明的保护范围之内。

Claims (10)

1.爆炸动焦散斑特征识别提取及应力强度因子自适应算法，其特征在于，具体包括以下步骤：

S1:设计装药结构、制作药包并进行炮孔装药固定；

S2：开展爆炸加载动态焦散线实验并拍摄记录整个过程；

S3：批量读取测试得到的待分析爆炸加载实验焦散斑图像；

S4：对批量待分析焦散斑图像进行Mask匀光处理；

S5：将图像转换为灰度图并进行非局部均值滤波；

S6：对图像采用canny算法进行焦散斑边缘检测；

S7：采用Hough变换进行焦散斑识别并定位所识别圆形轮廓的中心坐标；

S8：将满足直径最大和最小圆方程等式的中心点坐标进行标识并存储；

S9：根据坐标定位在图像中绘制焦散斑最大、最小直径圆周轮廓；

S10：根据相关公式自动计算所提取焦散斑的直径、面积、周长等特征参数并存储输出；

S11：计算获取裂纹扩展速度、加速度、裂纹尖端应力强度因子数值随时间变化曲线并存储输出；

S12：计算确定并存储输出裂纹动态能量释放率随时间变化曲线。

2.根据权利要求1所述的爆炸动焦散斑特征识别提取及应力强度因子自适应算法，其特征在于，所述S2步骤中，动态焦散线实验的仪器主要包括激光光源、扩束镜、场镜即凸透镜、炸药起爆装置、高速相机、数据采集计算机；所述激光光源采用绿色激光，且光源的波长为532mm，所述激光光源的激光功率在0～200mW之间可调，所述扩束镜将激光斑点发散，扩大光场，两个凸透镜作为场镜，直径为30cm，焦距为90cm；第一透镜将发射光场调整为平行光场，第二透镜将平行光场汇聚到相机镜头处。

3.根据权利要求1所述的爆炸动焦散斑特征识别提取及应力强度因子自适应算法，其特征在于，所述S3步骤中，在所述S2中拍摄到的照片内的裂纹尖端位置按照一定的比例与实际的试件进行换算，通过程序自动测量每一时刻焦散斑的近似圆心，采用数值微分中心点公式可得每一时刻裂纹扩展速度，设测量的k组焦散斑近似圆心(xi,yi),每两幅焦散斑时间间隔为△t；

根据数值微分中心点公式计算：

其中：当i＝1时，

当i＝k时，

△t—相邻两张照片拍摄时间间隔。

4.根据权利要求1所述的爆炸动焦散斑特征识别提取及应力强度因子自适应算法，其特征在于，所述Mask匀光处理的Mask匀光算法数学模型可表示为：

f(x,y)＝r(x,y)+δ(x,y)；

式中，f(x,y)代表分析图像；δ(x,y)代表背景影像；r(x,y)代表亮度分布均匀的影像。

5.根据权利要求4所述的爆炸动焦散斑特征识别提取及应力强度因子自适应算法，其特征在于，所述Mask匀光处理的具体过程如下：

S41.背景影像的获取:背景影像的获取使用低通滤波来获得可以反映影像颜色和亮度变化的低频信息，作为背景影像；

S42.影像相减运算：将原始影像与背景影像做相减运算，得到剔除了低频信息的影像，为了保持匀光前后影像整体亮度基本一致，相减运算的同时还要加上代表原始影像亮度均值的偏移量:；

S43.反差拉伸处理：由于相减运算会导致图像灰度值的动态范围缩小，影像反差减小，对上述结果进行拉伸处理，增加影像整体反差，得到最终的匀光结果，进行拉伸处理的数学模型选取如下:

式中，f(x,y)表示最终的结果图像像，f'(x,y)表示相减运算后的图像，m为拉伸前图像灰度值均值，β为拉伸因子，

为拉伸系数，具体公式如下：

式中，

表示(x，y)像素点的拉伸系数，α为调节参数，f_bk(x，y)表示背景影像中(x，y)点的灰度值，f_bkmax表示背景影像的最大灰度值，f_bkmin表示背景影像的最小灰度值；建立信息熵、清晰度、平均对比度和块标准差四个评价指标并设定相应的阈值，当同时满足上述四项指标设定阈值便认为匀光效果理想，否则调整相关参数进行重新计算直至达到效果要求为止；

所述S5的步骤中的非局部均值滤波算的算法原理如下：

假设含噪图像为v＝{v(a)|a∈A}，去噪之后的图像为NL[v],每个像素a的加权值可用下式求出：

其中，w(a,b)表征像素a与像素b之间的相似度(高斯加权欧式距离)权值，满足0≤w(a,b)≤1且∑_bw(a,b)≤1，其表达式如下：

其中，

表示以像素a为中心的子块内像素和以像素b为中心的子块内像素的高斯加权距离的平方，用于衡量两个像素点之间的相似度，v(N(a))表示a周围的局部子块像素集合，h为滤波参数。

6.根据权利要求5所述的爆炸动焦散斑特征识别提取及应力强度因子自适应算法，其特征在于，所述非局部均值滤波主要步骤如下：

S11.在含噪图像中建立一个边长为D的搜索块，再建立一个边长为d的相似块；

S12.在搜索块中选择一个像素点i，以i为中心构建一个相似块；在搜索块中任取一点j构建同样大小的相似块；依据高斯加权欧式距离计算权重系数w(i,j)；

S13.在搜索窗中遍历所有点j；

S14.将权重的最大值赋予点i，并将权重归一化，通过加权求和公式计算去噪后第i点的像素值：

其中，NL[v](i)为采用非局部均值算法得到的新的像素灰度值；w(i,j)为权重系数，v(j)为像素点j的灰度值,w(i,j)的具体表达式为：

式中，v(N_i)和v(N_j)表示以像素点i和像素点j为中心构建的相似块；e指数部分代表高斯加权欧式距离；a代表构建高斯核时的加权系数；h是构建权重的参数，控制衰减速度，h越小边缘特点保持得越好，但噪声严重，反之图像光滑但边缘保持的较差；

S15.最后遍历图像中所有的i点，完成去噪过程并通过PSNR检验图像去噪效果。

7.根据权利要求1所述的爆炸动焦散斑特征识别提取及应力强度因子自适应算法，其特征在于，所述S6的步骤中，离散Canny算子边缘检测提取步骤如下：

S61.采用尺度为W的平滑算子s按行对图像做卷积运算，然后用相应的最优求导算子d对平滑后的算子按列做卷积运算，得到Gy；再用s按列对图像做卷积运算，然后d对平滑后的算子按行做卷积运算，得到Gx；

S62.计算图像水平和垂直的梯度幅值和方向，该方法中使用一阶偏导的有限差分来计算梯度幅值和方向，对于函数图像梯度大小和计算梯度方向计算如下：

方向梯度计算：

梯度幅值计算：

梯度方向计算：

S63.非极大值抑制：遍历G，把边缘的梯度方向大致分为水平、竖直、45°和135°四个方向，各个方向用不同的邻近像素进行比较，以决定局部极大值，若某个像素的灰度值与其梯度方向上前后两个像素的灰度值相比不是很大的，则将该像素值置为0，即非边缘，该过程称为非极大值抑制；

S64.双阈值操作：设置高低阈值，遍历G，凡是大于高阈值的设为边缘；凡是小于低阈值的设为非边缘，若检测结果介于高、低阈值之间，则要判断像素的邻接像素中是否存在超过高阈值的边缘像素，若有则设为边缘，否则设为非边缘；

S65：通过抑制孤立的弱边缘最终完成边缘检测。

8.根据权利要求1所述的爆炸动焦散斑特征识别提取及应力强度因子自适应算法，其特征在于，所述S7的步骤中，Hough圆检测算法利用圆形本身的几何特性，将二维像素点映射到高维的Hough空间，以一定的半径步长对得到的圆心和半径做检测分析，同时，引入聚类分析的思想，合并相似度高的圆形。

9.根据权利要求8所述的爆炸动焦散斑特征识别提取及应力强度因子自适应算法，其特征在于，所述Hough圆检测算法的具体步骤为：

S71.以图像左下角为坐标原点建立二维平面坐标系，对图像进行预处理，经过高斯去噪、二值化、边缘提取等步骤，将所建立坐标系下图像的边缘像素信息保存在集合S中；

S72.在集合S中随机地抽取一个像素点X，构造三个直角三角形，分别沿着水平方向和垂直方向构造之间三角形A1，沿着-45°和45°构造直角三角形A2，沿着-15°和75°构造直角三角形A3；

S73.A1,A2,A3的直径相交于一点B，以B为中心，得到X的对称点Y；

S74.如果A1,A2,A3的斜边和X与Y组成的直径满足圆形的几何性质，则将此圆的圆心坐标和半径r、直径D信息加入到集合C中；根据Hough变换识别提取到的圆心和半径，编制程序计算所识别标记的焦散斑的周长、面积并存储输出；具体公式为：

周长：L＝2πr＝πD

面积：

通过编制循环程序计算每一时刻焦散斑的直径D_max和D_min，可得裂纹尖端动态应力强度因子，动焦散试验方法中动态应力强度因子的值与焦散斑直径按照如下公式计算：

式中：

—I型裂纹尖端应力强度因子；

—Ⅱ型裂纹尖端应力强度因子；

D_max—焦散斑最大直径，由前述程序自动计算得到；

g—应力强度因子的数值系数，根据g～μ曲线确定；

F(v)—裂纹扩展速度的调节系数，通常在具有实际意义的裂纹速度下取值为1；

d_eff—试件有效厚度即材料实际厚度，根据试验材料实际厚度确定；

c—材料的应力光学常数，取为0.85×10-10m2/N；

z₀—参考平面到试件表面的距离，根据实验布置确定，通常取为100mm；

μ—Ⅰ、Ⅱ型动态应力强度因子之比，通过曲线得出。

10.根据权利要求9所述的爆炸动焦散斑特征识别提取及应力强度因子自适应算法，其特征在于，所述动态能量释放率作为裂纹扩展过程中的能量量度指标，是控制扩展裂纹尖端的应力场、应变场的参量，是裂纹扩展过程中弹性能量释放，反映了裂纹扩展的驱动力；平面状态下能量释放率计算公式如下：

式中：

α_d＝1-(v/c₁)²

α_s＝1-(v/c₂)²

其中，c₁，c₂分别为膨胀波波速和剪切波波速，v为裂纹扩展速度，可利用裂纹尖端强度因子计算得到裂纹尖端的动态能量释放率。

Images