基于高速摄影系统的冲击波超压场测量方法与流程

1.本发明涉及一种基于高速摄影系统的非接触测量冲击波超压场的测量方法。属于冲击波超压测试及数据分析领域。


背景技术：

2.冲击波超压测试是各类殉爆等安全性试验考核中常见的测试需求，常规测量方法是以爆心为圆心，在被试件周围由近及远布设超压传感器，爆炸后冲击波传播过程中会通过超压传感器感压面，从而记录该传感器布设位置处，该时刻点的冲击波超压数值，试验后将数据回收，根据采集到的各个点出的冲击波超压数值进行数据拟合得到冲击波超压数据场。
3.传感器测量方法是现在测量冲击波超压的主要方法，该方法可以获得殉爆点近场的有效数据，可进行冲击波超压数据场分析。但是传感器测量方法存在以下问题：1、传感器布设的方向、高度及角度对测试数据准确性存在较大影响；2、需要为传感器设计专用工装，安装难度大；3、因传感器均布设在近场位置处，被测点反应产生的试验碎片、石块等物品对测量可能产生不可控影响；4、试验环境恶劣，传感器存在较大损坏的可能性，会导致试验数据获得率达不到100％且试验成本增加。
4.在殉爆等需要进行冲击波超压采集的试验中，高速摄像系统作为辅助测量手段进行被测产品观测和视频数据采集，采集到的视频及图片数据中可以观测到冲击波波阵面的传播，因此通过对高速摄影系统数据进行分析提取来进行冲击波超压数据测试，具有良好的数据条件，同时高速摄影系统架设距离爆心较远，不受试验时环境影响，且采集数据完整，对冲击波超压数据分析提供了可靠基础。


技术实现要素：

5.要解决的技术问题
6.冲击波超压测试目前均采取冲击传感器近点测试方法，该方法存在传感器损坏、传感器支架摆设位置限制、试验品碎片等多重不可控因素干扰，测试成本高，数据获得率不能达到100％等问题。
7.高速摄影系统在冲击波超压测试中的应用，可以获取到冲击波波阵面传播的视频数据，目前仅将高速摄影系统用于提供数据。为实现从高速摄影系统获取的冲击波视频数据中计算冲击波超压场计算，本发明需要解决的技术问题包括：
8.1、采用高速摄影系统拍摄的冲击波传播视频，如何将可视化视频数据转化为可进行数值计算的数值数据。
9.2、冲击波波阵面数据识别问题：如何在视频中，将冲击波传播的波阵面从视频中分离出来，转换为时间——距离曲线。
10.3、如何通过冲击波波阵面的时间——距离曲线获得冲击波不同时刻点距离——速度曲线。
11.4、如何通过冲击波距离——速度曲线计算各个时刻冲击波超压场数据。
12.5、如何将计算得到的冲击波超压数据转化为冲击波超压场特性。
13.本发明的技术方案为：
14.所述一种基于高速摄影系统的冲击波超压场测量方法，包括以下步骤：
15.步骤1：拍摄冲击波超压测试视频，获取完整的冲击波传播画面，获取的图像中，图像坐标系的横坐标平行于地平线，纵坐标为高度方向；
16.步骤2：提取冲击波波阵面：
17.步骤2.1：在高速摄像系统采集到的冲击波传播图像中，选择冲击波波阵面与爆炸所产生的火球烟雾完全分离的图像作为处理图像；
18.步骤2.2：将所述处理图像转化为灰度图像；
19.步骤2.3：采用边缘检测方法对灰度图像进行边缘检测提取，得到边缘图像；
20.步骤2.4：通过采用后一帧边缘图像减去前一帧边缘图像的方式消除背景边缘信息；
21.步骤2.5：根据冲击波波阵面边缘在外侧，火球及烟雾边缘在内侧的原理，对经过步骤2.4处理后的边缘图像进行有效边缘提取，得到冲击波有效波阵面边缘曲线；有效波阵面边缘曲线在该图像的图像矩阵中数据值为1，图像矩阵中的其他像素点数据值置0；
22.步骤3：确定冲击波波阵面时间
‑
距离序列：
23.选取冲击波传播方向线：利用图像中已知的爆心坐标(x0，y0)，并根据边缘检测结果中冲击波波阵面传播趋势，选定远端像素坐标(x
n
，y
n
)，该点坐标选定满足条件：连接爆心与该点所确定的直线，与冲击波传播过程中的波阵面均相交；
24.该直线方程表达式为：
[0025][0026]
表示为：ax+by+c＝0
[0027]
其中：
[0028][0029]
b＝
‑1[0030][0031]
对于该图像所对应的经过步骤2处理得到的图像矩阵，提取该矩阵中所有数值为1的像素坐标，形成冲击波波阵面序列(x
′
i
，y
′
i
)，并计算冲击波波阵面序列中各像素点到冲击波传播方向线的距离：
[0032][0033]
得到该图像中，冲击波波阵面与冲击波传播方向线相交像素点坐标为min(d
i
)对应的坐标(x
′
i
，y
′
i
)；
[0034]
进一步得到该图像中，冲击波波阵面与冲击波传播方向线相交像素点同爆心的实际距离为
[0035][0036]
对含有冲击波波阵面的图像均进行上述处理，其中选择的远端像素坐标(x
n
，y
n
)是相同的，得到时间
‑
距离序列[(t
j
，l
j
)|j＝1，2
…
]；
[0037]
对时间
‑
距离序列逐个进行判断，当满足条件l
j
＜l
j+1
，则认为第j个点是冲击波波阵面上的有效点，直至不满足条件l
j
＜l
j+1
，则得到之前的有效的时间
‑
距离序列为：[(t
j
，l
j
)|j＝1..k]，其中k表示有效图片数量，t
j
为对应图片j的相对时间，l
j
为图片j中冲击波波阵面与爆心的实际距离；
[0038]
步骤4：确定冲击波波阵面距离
‑
速度序列：
[0039]
对步骤3得到的有效的时间
‑
距离序列进行差分计算可得到冲击波的速度序列，进而得到对应的冲击波波阵面距离
‑
速度序列；
[0040]
步骤5：确定冲击波波阵面距离
‑
超压值序列：
[0041]
根据波阵面压力变化公式
[0042][0043]
得到冲击波波阵面速度
‑
超压值序列，进而得到冲击波波阵面距离
‑
超压值序列，其中m
s
表示冲击波相对于波前气流运动的马赫数，根据公式确定，p为峰前超压值，k为波前和波后空气比热比，d为波前介质本身速度，v为冲击波波阵面运动速度，c为测试环境温度下的空气中的声速，p0为测试环境大气压；
[0044]
步骤6、确定冲击波超压场特性：
[0045]
将得到的冲击波波阵面距离
‑
超压值序列代入萨道夫斯基经验公式，解算萨道夫斯基经验公式系数，从而得到冲击波超压场特性公式。
[0046]
进一步的，在进行冲击波超压测试前，根据产品装药量，通过萨道夫斯基经验公式和火球经验公式，计算冲击波传播范围以及火球直径；依据计算结果，将高速摄影系统架设在距离产品一定距离处，并调整高速摄影系统水平及视野，确保高速摄影系统可以获取完整的冲击波传播画面。
[0047]
进一步的，距离转换系数其中l为在产品同一竖直平面内设立的两个标杆之间的实际距离，l1为图像中两个标杆间的距离l1；所述竖直平面垂直于高速摄影系统的光轴。
[0048]
进一步的，高速摄影系统的拍摄帧率不小于5000帧，像素不小于1280*800。
[0049]
进一步的，步骤2.3中，边缘检测方法采用canny边缘检测方法。
[0050]
进一步的，步骤2.5中，得到外侧冲击波波阵面边缘的数据点，形成有效波阵面边缘曲线的过程为：
[0051]
确定图像中冲击波波阵面纵坐标的最小值y
min
和最大值y
max
；在纵坐标[y
min
，y
max
]范围里，对于每个纵坐标点，均在水平方向从图像一侧向另一侧找到该水平线上第一个像素值为1的数据点并保存，从而得到冲击波一侧的有效波阵面边缘曲线。
[0052]
进一步的，步骤6中，萨道夫斯基经验公式的秦勒展开表达式为：
[0053][0054]
式中w为冲击波的tnt当量，将得到的冲击波波阵面距离
‑
超压值序列代入公式，通过最小二乘拟合，得到系数a1、a2、a3，从而得到冲击波超压场特性公式。
[0055]
有益效果
[0056]
本发明提出的算法优化了仅使用边缘检测带来的无关信息干扰，从图像中提取出了有效的冲击波波阵面曲线，并给出了切实可行的边缘提取算法。
[0057]
通过本发明算法，实现了将高速摄影系统采集到的冲击波波阵面转换成了可进行数值计算的数值问题，并通过边缘检测、提取和数据筛选算法提取到了有效波阵面时间——距离曲线，计算得到冲击波传播速度，最后通过拟合得到冲击波场的建立。本发明提出的算法优化了仅使用边缘检测带来的无关信息干扰，从图像中提取处了有效的冲击波波阵面曲线，并给出了切实可行的边缘提取算法。
[0058]
本发明的附加方面和优点将在下面的描述中部分给出，部分将从下面的描述中变得明显，或通过本发明的实践了解到。
附图说明
[0059]
本发明的上述和/或附加的方面和优点从结合下面附图对实施例的描述中将变得明显和容易理解，其中：
[0060]
图1：本发明的算法流程图；
[0061]
图2：高速摄影系统采集到的冲击波图片；
[0062]
图3：canny算子进行边缘检测后图片；
[0063]
图4：去除背景边缘后图片；
[0064]
图5：有效冲击波左侧边缘；
[0065]
图6：时间
‑
距离趋势图；
[0066]
图7：时间
‑
速度趋势图；
[0067]
图8：距离
‑
速度趋势图；
[0068]
图9：速度
‑
超压趋势图；
[0069]
图10：距离
‑
超压趋势图；
[0070]
图11：数据拟合结果。
具体实施方式
[0071]
下面详细描述本发明的实施例，所述实施例是示例性的，旨在用于解释本发明，而不能理解为对本发明的限制。
[0072]
本实施例提出的基于高速摄影系统的冲击波超压场测量方法包括以下内容：
[0073]
1、冲击波超压视频拍摄：
[0074]
冲击波超压测试试验前，根据产品装药量，通过萨道夫斯基经验公式和火球经验公式，计算冲击波传播范围以及火球直径。
[0075]
根据计算结果，将高速摄影系统架设在距离产品一定距离处，并调整高速摄影系统水平及视野，确保高速摄影系统可以获取完整的冲击波传播画面。获取的图像中，图像坐
标系的横坐标平行于地平线，纵坐标为高度方向。
[0076]
同时还要获取图像中的像素距离，可以通过在产品同一竖直平面内设立两个标杆，该竖直平面垂直于高速摄影系统的光轴；测量两个标杆间的距离l，并得到图像中两个标杆间的距离l1，从而得到距离转换系数为
[0077]
而为了能够准确反映冲击波传播过程，高速摄影系统的拍摄帧率不小于5000帧，像素不小于1280*800。
[0078]
2、冲击波波阵面提取：
[0079]
在高速摄像系统采集到的冲击波传播画面中的一帧图像中，冲击波边缘为一条肉眼可识别的浅黑色细线，因此可以采用图像分析手段提取冲击波边缘特征并进行位置跟踪计算。而冲击波传播的过程，波阵面处空气会发生压缩导致两侧密度发生显著变化，从而引起视觉上可观测的冲击波传播过程，该变化在图像上的显示即为图像矩阵数值的变换。
[0080]
2.1：图像筛选：选择冲击波波阵面与爆炸所产生的火球烟雾完全分离的图像作为处理图像，即从图像中可以看到冲击波波阵面的图像。
[0081]
2.2：图像转换：采集图像为rgb彩色图像，图像由r、g、b三个矩阵组成，在同一点处，三个数值变化不具有统一性，为得到冲击波波阵面的变化趋势，将图像转化为灰度图像，即完成了彩色图像中各点的归一化处理，将rgb数值变化转换为灰度值变化。
[0082]
从灰度图像中可观测到：背景山脉、爆炸火球、燃烧烟雾及冲击波波阵面均为灰度突变点，即：在上述位置处，左右两侧灰度值与该位置处存在明显差异，所以冲击波波阵面提取转化为图像中灰度突变点位置识别与提取。
[0083]
2.3：图像边缘识别：可以采用canny边缘检测方法对灰度图像进行边缘检测提取，得到边缘图像，在边缘检测图像中，冲击波波阵面、火球烟雾边缘及山脉等边缘均被检测出来。且该边缘图像矩阵组成内容为：0和1，即白色部分为检测到的边缘信息，数值为1；黑色部分为非边缘部分，数值为0。
[0084]
2.4：背景边缘去除：在得到的边缘图像中，冲击波波阵面和火球烟雾边缘为变化部分，山脉等环境边缘为固定且不关注部分。其中山脉等环境边缘信息属于背景边缘，在不同图像中稳定存在，所以通过后一帧边缘检测图像减去前一帧边缘检测图像，即可消除大部分固定的背景边缘信息。剩余边缘主要为变化的冲击波波阵面和火球烟雾边缘。
[0085]
2.5：有效边缘提取：在爆炸产生后，冲击波传播速度大于火球及烟雾传播速度，因此，上述边缘图像中，冲击波波阵面边缘在外侧，火球及烟雾边缘在内侧。
[0086]
可以采用以下过程得到外侧冲击波波阵面边缘的数据点，形成有效波阵面边缘曲线：
[0087]
通过观察确定图像中冲击波波阵面纵坐标的最小值y
min
和最大值y
max
；由于图像坐标系横坐标与地平线平行，所以在纵坐标[y
min
，y
max
]范围里，对于每个纵坐标点，均在水平方向从一侧向另一侧找到该水平线上第一个像素值为1的数据点并保存，例如对于纵坐标点y
min
，在水平方向从左向右找到该纵坐标点对应水平线上第一个像素值为1的数据点并保存，这样得到了冲击波左侧部分的有效波阵面边缘曲线。有效波阵面边缘曲线在该图像的图像矩阵中数据值为1，图像矩阵中的其他像素点数据值置0当然由于噪点存在，得到的有效波阵面边缘曲线中也可能存在噪点，但不影响后续处理过程；此外，本发明所针对的冲击
波超压测试是在户外场景下的冲击波超压测试，冲击波各向均匀，因此只取冲击波左侧部分的有效波阵面边缘曲线即可满足要求。
[0088]
3、确定冲击波波阵面时间
‑
距离序列：
[0089]
选取冲击波传播方向线：
[0090]
利用图像中已知的爆心坐标(x0，y0)，并根据边缘检测结果中冲击波波阵面传播趋势，选定远端像素坐标(x
n
，y
n
)，该点坐标选定满足条件：连接爆心与该点所确定的直线，与冲击波传播过程中的波阵面均相交。
[0091]
该直线方程表达式为：
[0092][0093]
可表示为：ax+by+c＝0
[0094]
其中：
[0095][0096]
b＝
‑1[0097][0098]
对于该图像中经过步骤2处理得到的图像矩阵，提取该矩阵中所有数值为1的像素坐标，形成冲击波波阵面序列(x
′
i
，y
′
i
)，并计算冲击波波阵面序列中各像素点到冲击波传播方向线的距离：
[0099][0100]
可以认为该图像中，冲击波波阵面与冲击波传播方向线相交像素点坐标为min(d
i
)对应的坐标(x
′
i
，y
′
i
)
[0101]
实际上，只有当d
i
＝0时，得到的交点才是实际交点，但由于像素点均是整数，所以得到的(x
′
i
,y
′
i
)是最接近实际交点的像素点，而且因为图像的像素值较大，所以偏移角度实际上是极小的，可以忽略误差影响。这样得到该图像中，冲击波波阵面与冲击波传播方向线相交像素点同爆心的实际距离为
[0102][0103]
对含有冲击波波阵面的图像均进行上述处理，其中选择的远端像素坐标(x
n
,y
n
)是相同的，这样就得到时间
‑
距离序列[(t
j
,l
j
)|j＝1,2
…
n]。
[0104]
因为沿冲击波传播方向，波阵面与爆心距离逐渐增大，由于序列中，时间是依次增大的，因此当满足条件l
j
＜l
j+1
，则认为第j个点是冲击波波阵面上的有效点，根据时间顺序依次判断，直至不满足条件l
j
＜l
j+1
，则得到之前的有效的时间
‑
距离序列为：[(t
j
,l
j
)|j＝1
…
k]，其中k表示有效图片数量，t
j
为对应图片j的相对时间，l
j
为图片j中冲击波波阵面与爆心的实际距离。由于拍摄帧率很高，因此，即使因为出现噪点而导致不满足条件l
j
＜l
j+1
，但之前已经得到的有效距离序列数量足够进行后续运算。
[0105]
4、确定冲击波波阵面距离
‑
速度序列：
[0106]
对步骤3得到的有效的时间
‑
距离序列进行差分计算可得到冲击波的速度序列，进而得到对应的冲击波波阵面距离
‑
速度序列。其中速度计算公式为：
[0107][0108]
5、确定冲击波波阵面距离
‑
超压值序列：
[0109]
由等熵条件下一维运动的连续方程、动量方程和能量方程可导出冲击波的状态参数变化基本关系式，包括
[0110]
波阵面温度变化：
[0111][0112]
波阵面密度变化：
[0113][0114]
波阵面压力变化：
[0115][0116]
其中m
s
表示冲击波相对于波前气流运动的马赫数，根据公式确定，p为峰前超压值，k为波前和波后空气比热比，d为波前介质本身速度，v为冲击波波阵面运动速度，c为测试环境温度下的空气中的声速，p0为测试环境大气压。
[0117]
这样利用波阵面压力变化公式得到冲击波波阵面速度
‑
超压值序列，进而得到冲击波波阵面距离
‑
超压值序列。
[0118]
6、确定冲击波超压场特性：
[0119]
冲击波超压传播符合萨道夫斯基经验公式，该公式秦勒展开表达式为：
[0120][0121]
式中w为冲击波的tnt当量，将得到的冲击波波阵面距离
‑
超压值序列代入公式，通过最小二乘拟合，得到系数a1、a2、a3，从而得到冲击波超压场特性公式。
[0122]
本实施例中，通过数据拟合可以得到，测试冲击波超压场分布规律为：
[0123][0124]
此计算结果与超压传感器近场测试数据吻合度良好。
[0125]
尽管上面已经示出和描述了本发明的实施例，可以理解的是，上述实施例是示例性的，不能理解为对本发明的限制，本领域的普通技术人员在不脱离本发明的原理和宗旨的情况下在本发明的范围内可以对上述实施例进行变化、修改、替换和变型。