1.本发明属于图像处理技术领域，具体涉及一种基于图像处理技术的破片飞散特性测试方法及系统。
背景技术：2.在爆炸场中，武器战斗部破片速度、形状、分布等参数是评价其毁伤效能的重要依据。由于战斗部破片的类型、干扰因素、测试环境、测试参数等要求不尽相同，对破片的测试原理和测试手段都提出更高的要求。
3.目前靶场的破片测速方法主要有区截测速法、雷达测速法和高速摄影法。区截测速法是在国内靶场测试中较多采用的一种方法，主要采用定距测时的方法进行测量。由于这种方法的布设和数据采集费时费力，目前国内破片测试正向光电磁等非接触测量方法发展。雷达测速法由于雷达回波分析较为复杂，只适应于较为分散的低速飞行破片。利用高速摄影相机作为光电采集器件，能够根据破片的运动图像实现对破片目标空间位置的测量，具有测量系统布局方便、实景记录、非接触式测量、测量精度高、分析速度快等特点。因此本技术提出了一种基于图像处理的破片飞散特性测试方法。
技术实现要素：4.本发明提供了一种基于图像处理技术的破片飞散特性的测试方法及系统。
5.实现本发明的技术解决方案为：一种基于图像处理技术的破片飞散特性测试方法，包括以下步骤：
6.步骤1：根据战斗部当量计算破片分布半径，按照破片分布半径安全区域布置双目相机，对相机进行标定与参数设置，采集破片飞行的轨迹序列图像；
7.步骤2：对采集到的轨迹序列图像进行破片的自动识别与坐标提取；
8.步骤3：对提取到的破片在两个相机视场中互相匹配，利用提取的破片在图像中的二维坐标计算破片在空间的三维坐标；
9.步骤4：在三维空间中对破片运动方程进行拟合，根据拟合方程和已知的破片三维空间坐标求解破片运动参数。
10.进一步的，按照破片分布半径安全区域布置双目相机，对相机进行标定与参数设置，包括以下步骤：
11.步骤s11，采用两个高速摄像机，且两个相机型号相同，都满足测试要求的高帧频，按照破片分布半径安全区域布置双目相机，且使两个摄像机的视场交汇于探测区域；
12.步骤s12，在试验中使用棋盘格图像进行标定，两台摄像机对同一棋盘格进行同时拍摄，求解两个摄像机坐标系之间的转换关系，其中棋盘格上任意点的三维空间位置参数通过gps、全站仪获得，该任意点在二维图像中的位置通过特征点提取算法获得；
13.步骤s13，为2台摄像机分别连接gps授时设备，通过irig-b码将2个设备同步为gps时间，同时，将相机设定b码同步模式，采集完成后，只需将2 个图像序列第1张图像的时间
对齐，保证时间的同步，若存在丢帧现象，其中一个相机采集的视频信息找不到另一个对应的第一幅图像，则按照时间序列向后取帧，取相同时刻对应图像都被存储的时刻作为开始，保证时间同步；
14.步骤s14，获取破片飞行的轨迹序列图像。
15.进一步的，对采集到的序列图像进行破片的自动识别与坐标提取，包括以下步骤：
16.步骤s21，对双目相机采集到的破片飞行轨迹序列图像进行预处理去噪，针对小波系数的分布情况，对不同的频带分别进行处理，设定阈值为软阈值，当小波系数的绝对值小于给定阈值时，令其为0，大于阈值时，令其都减去阈值，基本原理如下：
[0017][0018]
式(1)中，sign(w)为符号函数,w是原始图像小波系数，w
λ
是小波变换后重新生成的小波系数，λ是阈值；
[0019]
步骤s22，设fk、f
k+1
、f
k+2
为三帧相邻的图像，对三帧图像间分别进行差分，得到：
[0020][0021]
式(2)中，s1、s2、s3为经差分得到的目标前景灰度矩阵；
[0022]
步骤s23，对差分后的前景进行叠加，就可以增大目标所在位置的灰度值，提高目标与噪声的灰度差。对差分后前景灰度矩阵进行叠加，得到：
[0023]dk
＝s1+s2+s3ꢀꢀ
(3)
[0024]
式(3)中，dk为差分后前景灰度矩阵的叠加矩阵；
[0025]
步骤s24，为了进一步增大目标前景同噪声的对比度，对dk每个像素点进行灰度值平方处理，得到：
[0026][0027]
式(4)中，gk为叠加矩阵的灰度值平方处理矩阵；
[0028]
步骤s25，采用改进的局部阈值法，将目标前景从噪声中分割出来，先将灰度值平方处理得到的图像划分为子图像，再对每一个子图像进行如下操作：
[0029]
获取子图像f(i,j)的直方图；
[0030]
按灰度级从小到大遍历直方图，找到第一个峰和最后一个峰，相应的灰度值分别记为t0和t1，记tmin、tmax分别为子图像的最小、最大灰度；
[0031]
设f(i,j)为子图像中像素点的灰度值，若f(i,j)》t0且f(i,j)《t1，则用
[0032]
计算新的灰度值代替子图像灰度值，式中t为针对该子图像的最大类间方差法确定的阈值；
[0033]
将所有按照上述处理子图像拼接成原始大小，得到背景与目标分离的图像；
[0034]
步骤s26，进行形态学处理，去除差分和分割处理过程中出现的孤立噪点和修补空洞；
[0035]
步骤s27，对图像中的破片的数量作统计，并依据破片在每一条轨迹中飞行时出现
在图像中的形状对破片形状作推测描述。
[0036]
进一步的，对提取到的破片在两个相机视场中互相匹配，利用提取的二维坐标计算破片的三维坐标，包括以下步骤：
[0037]
步骤s31，对得到的背景与目标分离的图像利用harris角点提取算法提取角点并用zncc策略对角点进行初始匹配；
[0038]
步骤s32，进行sift特征点检测与匹配，利用由sift特征点计算得到的基础矩阵对sift匹配点以及初始匹配角点对步骤s31处理的结果进行误匹配的剔除；
[0039]
步骤s33，对于误匹配剔除后的匹配图像，由匹配对利用八点法计算基础矩阵，由基础矩阵得到本质矩阵；
[0040]
步骤s34，对本质矩阵进行分解，求得相机外参数；
[0041]
步骤s35，对优化后的sift特征点以及角点进行三维空间点三角化；
[0042]
步骤s36，通过标定得到左右相机的内参a
l
与ar，两相机间旋转矩阵r，根据单目成像原理和坐标变换原理，有如下方程：
[0043][0044]
式(5)中，xw，yw，zw为破片在世界坐标系坐标，u1,v1和u2,v2分别为破片在左右相机的像素坐标系坐标，s1,s2分别为左右相机的工作距离参量，t
l
，tr分别为左右相机从世界坐标系转化到相机坐标系的偏移向量，r
l
,rr分别为左右相机从世界坐标系转化到相机坐标系的旋转矩阵；
[0045]
假设左相机坐标系是测量坐标系，即世界坐标系的位置，此时左边相机坐标系到世界坐标系的旋转矩阵为r
l
＝e，平移向量为t
l
＝0，即左相机坐标系和世界坐标系重合，平移向量为0，将相机内参矩阵和外参矩阵进行化简得到下式：
[0046][0047]
用最小二乘法计算三维点坐标，基于以下假设：空间任意点p在两个摄像机 c1与c2上的图像点平p1与p2已经从两个图像中分别检测出来，即已知p1与p2为空间同一点p的对应点；摄像机c1与c2已标定，对应的投影矩阵分别为m1与m2，于是有：
[0048][0049]
其中，(u1,v1,1)与(u2,v2,1)分别为p1与p2点在各自图像中的图像齐次坐标； (x,y,z,1)为p点在世界坐标系下的齐次坐标；m
ijk
为mk的第i行第j列元素，消去s1与s2得到关于xw，yw，zw的四个线性方程，即可求出三维点p的坐标。
[0050]
进一步的，根据拟合方程和已知的破片三维空间坐标求解破片运动参数，包括以下步骤：
[0051]
步骤s41，构造破片在空气中的运动微分方程：
[0052][0053]
式(8)中，m为破片质量，c
x
为破片迎面阻力系数，ρ为当地空气密度，s为破片垂直于飞行方向上的迎风面积，t为破片飞行时间，v为破片瞬时速度；
[0054]
步骤s42，获取弹体静爆前坐标和爆炸时刻，结合获得的破片运动轨迹，采取间接解算的方法解算；
[0055]
首先，将轨迹的第一个时刻点作为相对时刻，其位置作为相对起始位置，则有：
[0056][0057]
式中：t1为破片轨迹第一个数据点对应的时刻，x
′
为t1时刻该破片与爆炸点的距离，v1对应于t1时刻该破片的速度，v0对应于t0时刻该破片的速度， x为破片在某一时刻距离爆炸点的距离，α为速度衰减系数，则有：
[0058][0059]
式(10)中，x(t)为破片在任意t时刻飞行时对应飞行距离，等效为：
[0060][0061]
式(11)中：s
tr
为破片相对t1时刻坐标的距离，tr为相对于t1时刻的时间差，利用获取得到的破片轨迹数据，即可拟合得到速度v1和速度衰减系数α；
[0062]
计算战斗部质心至上述破片飞行的拟合直线的垂足坐标，由于破片拟合直线不一定经过战斗部质心，因此以垂足坐标作为破片运动的起始点，在此基础上，计算轨迹第一个数据点与垂足间距离即为x
′
；
[0063]
将速度v1、速度衰减系数α、距离x
′
代入式(9)，解算得到初速v0。
[0064]
一种基于图像处理技术的破片飞散特性测试系统，通过所述的图像处理技术的破片飞散特性测试方法，实现基于图像处理技术的破片飞散特性测试。
[0065]
一种计算机设备，包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序，所述处理器执行所述计算机程序时，通过所述的图像处理技术的破片飞散特
性测试方法，实现基于图像处理技术的破片飞散特性测试。
[0066]
一种计算机可读存储介质，其上存储有计算机程序，所述计算机程序被处理器执行时，通过所述的图像处理技术的破片飞散特性测试方法，实现基于图像处理技术的破片飞散特性测试。
[0067]
本发明与现有技术相比，其显著优点在于：1)在三维空间内解算破片运动飞散特性，更加精确；2)基于双目视觉原理的非接触式测量系统，较传统靶板法布设简易，数据采集更为迅速；3)破片目标提取时采用局部最大类间方差法，能有效避免工作现场光照不均匀导致的目标提取丢失情况。
附图说明
[0068]
图1为本发明的系统流程图。
[0069]
图2为本发明的双目相机布局图。
[0070]
图3为平方增强的三帧差分法流程图。
具体实施方式
[0071]
为了使本技术的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本技术进行进一步详细说明。应当理解，此处描述的具体实施过程仅仅用以解释本技术，并不用于限定本技术。
[0072]
如图1所示，一种基于图像处理技术的破片飞散参数确定方法，步骤如下：
[0073]
步骤1，根据战斗部当量计算破片分布半径，按照破片分布半径安全区域布置双目相机，如图2所示，对相机进行标定与参数设置，采集破片飞行时经过两个相机交汇视场的图像。具体操作如下：
[0074]
根据战斗部当量计算破片分布半径，按照破片分布半径安全区域布置双目相机，考虑到目标破片初速在800m/s-1500m/s，选用高速摄影机拍摄频率为1000 s-1，镜头焦距约为165mm，图像分辨率为1280
×
800像素，像元尺寸为20μm，双相机横向公共视场约为10m，纵向公共视场约为5m，单像元空间分辨率约为 7mm。两相机距离目标均为70m，交会角为60
°
。
[0075]
为2台高速相机分别连接gps授时设备，通过irig-b码将2个设备同步为 gps时间。同时，将相机设定b码同步模式。采集完成后，只需将2个图像序列第1张图像的时间对齐。若存在丢帧现象，其中一个相机采集的视频信息找不到另一个对应的第一幅图像，则按照时间序列向后取帧，取相同时刻对应图像都被存储的时刻作为开始，保证时间同步。
[0076]
双目相机的触发可以选择手动触发或声音/振动触发。手动触发需要通过网线将高速相机与控制计算机连接，建立局域网，通过计算机手动控制触发。
[0077]
步骤2，对收集到的序列图像进行破片的自动识别与坐标提取；
[0078]
对双目相机采集到的破片图像进行预处理去噪，可针对小波系数的分布情况，设计滤波器对不同的频带分别进行处理。设定阈值为软阈值，当小波系数的绝对值小于给定阈值时，令其为0，大于阈值时，令其都减去阈值，基本原理如下：
[0079]
[0080]
式(1)中，sign(w)为符号函数，w是原始图像小波系数，w
λ
是小波变换后重新生成的小波系数，λ是阈值。
[0081]
设fk、f
k+1
、f
k+2
为三帧相邻的图像，对三帧图像间分别进行差分可得：
[0082][0083]
式(2)中,s1、s2、s3为经差分得到的目标前景灰度矩阵。
[0084]
由于噪声是随机的，不同位置的噪声大小不同，同一位置在不同时刻的噪声大小也不尽相同，而破片目标在三张邻帧图像中的位置基本是确定的，对差分后的前景进行叠加，就可以增大目标所在位置的灰度值，提高目标与噪声的灰度差。对差分后前景灰度矩阵进行叠加可得：
[0085]dk
＝s1+s2+s3ꢀꢀ
(3)
[0086]
式(3)中，dk为差分后前景灰度矩阵的叠加矩阵；
[0087]
步骤s24，为了进一步增大目标前景同噪声的对比度，对dk每个像素点进行灰度值平方处理，可得：
[0088][0089]
式(4)中，gk为叠加矩阵的灰度值平方处理矩阵；
[0090]
目标前景同噪声的对比度得到了增强，需要设定一个阈值，将目标前景从噪声中分割出来。由于图像光照不均匀、有突发噪声或背景灰度变化较大时，简单的二值化阈值并不能满足破片与背景分离的要求，因此采用改进的局部阈值法。局部阈值过程如下：
[0091]
两个ccd相机拍到的图像尺寸是一样的，灰度值平方处理得到的图像和原始图像均为1280*800，根据尺寸将灰度值平方处理得到的图像划分为128*80 的子图像，每一幅灰度平方图像划分为100张子图像；对每一个子图像进行如下操作：
[0092]
获取子图像f(i,j)的直方图；
[0093]
按灰度级从小到大遍历直方图，找到第一个峰和最后一个峰，相应的灰度值分别记为t0和t1。记tmin、tmax分别为子图像的最小、最大灰度。
[0094]
设f(i,j)为子图像中像素点的灰度值，若f(i,j)》t0且f(i,j)《t1,则用
[0095]
计算新的灰度值代替子图像灰度值，式中t为针对该子图像的最大类间方差法确定的阈值。
[0096]
将所有按照上述处理子图像拼接成原始大小，得到背景与目标分离的图像。
[0097]
进一步地，类间方差定义为：
[0098]
σ2＝θ(μ-μ1)2+(1-θ)(μ
2-μ)2＝θ(1-θ)(μ
1-μ2)2[0099]
式中，c1类的灰度值i∈[0，m]，c2类的灰度值i∈[m+1,l-1]，μ表示图像的总均值，μ1表示c1的类均值，μ2表示c2的类均值,θ为灰度值在c1中的出现概率。使σ2最大的阈值m即为最佳阈值。
[0100]
形态学处理去除差分和分割处理过程中出现的孤立噪点和修补空洞。
[0101]
通过上述步骤可以将破片从采集到的图像背景中提取，进而得到破片在每一帧图像中的二维坐标信息，对图像中的破片的数量作统计，并依据破片在每一条轨迹中飞行时出现在图像中的形状对破片形状推测描述。
[0102]
步骤3：对提取到的破片在两个相机视场中互相匹配，利用提取的二维坐标计算破片的三维坐标；
[0103]
利用harris角点提取算法提取角点并用zncc策略对角点进行初始匹配，利用由sift特征点计算得到的基础矩阵对sift匹配点以及初始匹配角点进行误匹配的剔除；
[0104]
进行sift特征点检测与匹配；
[0105]
由匹配对利用八点法计算基础矩阵，由基础矩阵可以得到本质矩阵；
[0106]
对本质矩阵进行分解即可求得相机外参数；
[0107]
对优化后的sift特征点以及角点进行三维空间点三角化；
[0108]
通过标定得到内参a
l
与ar,两相机间旋转矩阵r，平移向量t。根据单目成像原理和坐标变换原理，有如下方程：
[0109][0110]
式(5)中，xw，yw，zw为破片在世界坐标系坐标，u1,v1和u2,v2分别为破片在左右相机的像素坐标系坐标，s1,s2分别为左右相机的工作距离参量，t
l
，tr分别为左右相机从世界坐标系转化到相机坐标系的偏移向量，r
l
,rr分别为左右相机从世界坐标系转化到相机坐标系的旋转矩阵；
[0111]
假设左相机坐标系是测量坐标系，即世界坐标系的位置，此时左边相机坐标系到世界坐标系的旋转矩阵为r
l
＝e，平移向量为t
l
＝0，即左相机坐标系和世界坐标系重合，平移向量为0，将相机内参矩阵和外参矩阵进行化简得到下式：
[0112][0113]
用最小二乘法计算三维点坐标，我们要明确一下几点基本假设：空间任意点 p在两个摄像机c1与c2上的图像点平p1与p2已经从两个图像中分别检测出来。即已知p1与p2为空间同一点p的对应点。至于对应点是如何找到的，这属于图像配准问题。还有一个假设就是，摄像机c1与c2已标定，他们的投影矩阵分别为m1与m2，于是有：
[0114][0115]
其中，(u1,v1,1)与(u2,v2,1)分别为p1与p2点在各自图像中的图像齐次坐标； (x,y,z,1)为p点在世界坐标系下的齐次坐标；m
ijk
为mk的第i行第j列元素。式 7消去s1与s2得到关于xw，yw，zw的四个线性方程，此时即可求出三维点p的坐标。
[0116]
步骤四，根据图像序列与帧数模拟出破片飞行的轨迹，并对破片的数量与位置统计。
[0117]
在做出一定假设情况下构造破片在空气中的运动微分方程，假设如下：
[0118]
(1)忽略受到的重力作用，破片只在空气阻力的作用下做减速运动；
[0119]
(2)忽略破片旋转带来的影响，认为空气阻力系数不发生变化，空气阻力作用于破片质心；
[0120]
(3)忽略静爆瞬间空气气温变化对空气密度的影响。
[0121]
破片在空气阻力的作用下做减速运动，运动微分方程为：
[0122][0123]
式(8)中，m为破片质量，c
x
为破片迎面阻力系数，ρ为当地空气密度，s为破片垂直于飞行方向上的迎风面积，t为破片飞行时间，v为破片瞬时速度；
[0124]
获取弹体静爆前坐标和爆炸时刻，结合获得的破片运动轨迹，采取间接解算的方法解算；
[0125]
首先，将轨迹的第一个时刻点作为相对时刻，其位置作为相对起始位置，则有：
[0126][0127]
式中：t1为破片轨迹第一个数据点对应的时刻，x
′
为t1时刻该破片与爆炸点的距离，v1对应于t1时刻该破片的速度，v0对应于t0时刻该破片的速度， x为破片在某一时刻距离爆炸点的距离，α为速度衰减系数，则有：
[0128][0129]
式(10)中，x(t)为破片在任意t时刻飞行时对应飞行距离，等效为：
[0130][0131]
式(11)中：s
tr
为破片相对t1时刻坐标的距离，tr为相对于t1时刻的时间差，利用获取得到的破片轨迹数据，即可拟合得到速度v1和速度衰减系数α；
[0132]
计算战斗部质心至上述破片飞行的拟合直线的垂足坐标，由于破片拟合直线不一定经过战斗部质心，因此以垂足坐标作为破片运动的起始点，在此基础上，计算轨迹第一个数据点与垂足间距离即为x
′
；
[0133]
将速度v1、速度衰减系数α、距离x
′
代入式(9)，解算得到初速v0。
[0134]
综上所述，本发明可以高精度的测量破片飞散特性，具有布局方便精简，测量精度高，分析速度快等特点，为新型目标毁伤计算与分析提供有效的理论依据。
[0135]
以上实施例的各技术特征可以进行任意的组合，为使描述简洁，未对上述实施例中的各个技术特征所有可能的组合都进行描述，然而，只要这些技术特征的组合不存在矛盾，都应当认为是本说明书记载的范围。
[0136]
以上所述实施例仅表达了本技术的几种实施方式，其描述较为具体和详细，但并不能因此而理解为对发明专利范围的限制。应当指出的是，对于本领域的普通技术人员来说，在不脱离本技术构思的前提下，还可以做出若干变形和改进，这些都属于本技术的保护范围。因此，本技术专利的保护范围应以所附权利要求为准。
技术特征：1.一种基于图像处理技术的破片飞散特性测试方法，其特征在于，包括以下步骤：步骤1：根据战斗部当量计算破片分布半径，按照破片分布半径安全区域布置双目相机，对相机进行标定与参数设置，采集破片飞行的轨迹序列图像；步骤2：对采集到的轨迹序列图像进行破片的自动识别与坐标提取；步骤3：对提取到的破片在两个相机视场中互相匹配，利用提取的破片在图像中的二维坐标计算破片在空间的三维坐标；步骤4：在三维空间中对破片运动方程进行拟合，根据拟合方程和已知的破片三维空间坐标求解破片运动参数。2.根据权利要求1所述的基于图像处理技术的破片飞散特性测试方法，其特征在于，按照破片分布半径安全区域布置双目相机，对相机进行标定与参数设置，包括以下步骤：步骤s11，采用两个高速摄像机，且两个相机型号相同，都满足测试要求的高帧频，按照破片分布半径安全区域布置双目相机，且使两个摄像机的视场交汇于探测区域；步骤s12，在试验中使用棋盘格图像进行标定，两台摄像机对同一棋盘格进行同时拍摄，求解两个摄像机坐标系之间的转换关系，其中棋盘格上任意点的三维空间位置参数通过gps、全站仪获得，该任意点在二维图像中的位置通过特征点提取算法获得；步骤s13，为2台摄像机分别连接gps授时设备，通过irig-b码将2个设备同步为gps时间，同时，将相机设定b码同步模式，采集完成后，只需将2个图像序列第1张图像的时间对齐，保证时间的同步，若存在丢帧现象，其中一个相机采集的视频信息找不到另一个对应的第一幅图像，则按照时间序列向后取帧，取相同时刻对应图像都被存储的时刻作为开始，保证时间同步；步骤s14，获取破片飞行的轨迹序列图像。3.根据权利要求1所述的基于图像处理技术的破片飞散特性测试方法，其特征在于，对采集到的序列图像进行破片的自动识别与坐标提取，包括以下步骤：步骤s21，对双目相机采集到的破片飞行轨迹序列图像进行预处理去噪，针对小波系数的分布情况，对不同的频带分别进行处理，设定阈值为软阈值，当小波系数的绝对值小于给定阈值时，令其为0，大于阈值时，令其都减去阈值，基本原理如下：式(1)中，sign(w)为符号函数,w是原始图像小波系数，w
λ
是小波变换后重新生成的小波系数，λ是阈值；步骤s22，设f
k
、f
k+1
、f
k+2
为三帧相邻的图像，对三帧图像间分别进行差分，得到：式(2)中，s1、s2、s3为经差分得到的目标前景灰度矩阵；步骤s23，对差分后的前景进行叠加，就可以增大目标所在位置的灰度值，提高目标与噪声的灰度差。对差分后前景灰度矩阵进行叠加，得到：d
k
＝s1+s2+s3ꢀꢀꢀꢀ
(3)式(3)中，d
k
为差分后前景灰度矩阵的叠加矩阵；
步骤s24，为了进一步增大目标前景同噪声的对比度，对d
k
每个像素点进行灰度值平方处理，得到：式(4)中，g
k
为叠加矩阵的灰度值平方处理矩阵；步骤s25，采用改进的局部阈值法，将目标前景从噪声中分割出来，先将灰度值平方处理得到的图像划分为子图像，再对每一个子图像进行如下操作：获取子图像f(i,j)的直方图；按灰度级从小到大遍历直方图，找到第一个峰和最后一个峰，相应的灰度值分别记为t0和t1，记tmin、tmax分别为子图像的最小、最大灰度；设f(i,j)为子图像中像素点的灰度值，若f(i,j)＞t0且f(i,j)＜t1，则用计算新的灰度值代替子图像灰度值，式中t为针对该子图像的最大类间方差法确定的阈值；将所有按照上述处理子图像拼接成原始大小，得到背景与目标分离的图像；步骤s26，进行形态学处理，去除差分和分割处理过程中出现的孤立噪点和修补空洞；步骤s27，对图像中的破片的数量作统计，并依据破片在每一条轨迹中飞行时出现在图像中的形状对破片形状作推测描述。4.根据权利要求1所述的基于图像处理技术的破片飞散特性测试方法，其特征在于，对提取到的破片在两个相机视场中互相匹配，利用提取的二维坐标计算破片的三维坐标，包括以下步骤：步骤s31，对得到的背景与目标分离的图像利用harris角点提取算法提取角点并用zncc策略对角点进行初始匹配；步骤s32，进行sift特征点检测与匹配，利用由sift特征点计算得到的基础矩阵对sift匹配点以及初始匹配角点对步骤s31处理的结果进行误匹配的剔除；步骤s33，对于误匹配剔除后的匹配图像，由匹配对利用八点法计算基础矩阵，由基础矩阵得到本质矩阵；步骤s34，对本质矩阵进行分解，求得相机外参数；步骤s35，对优化后的sift特征点以及角点进行三维空间点三角化；步骤s36，通过标定得到左右相机的内参a
l
与a
r
，两相机间旋转矩阵r，根据单目成像原理和坐标变换原理，有如下方程：式(5)中，x
w
，y
w
，z
w
为破片在世界坐标系坐标，u1,v1和u2,v2分别为破片在左右相机的像
素坐标系坐标，s1,s2分别为左右相机的工作距离参量，t
l
，t
r
分别为左右相机从世界坐标系转化到相机坐标系的偏移向量，r
l
,r
r
分别为左右相机从世界坐标系转化到相机坐标系的旋转矩阵；假设左相机坐标系是测量坐标系，即世界坐标系的位置，此时左边相机坐标系到世界坐标系的旋转矩阵为r
l
＝e，平移向量为t
l
＝0，即左相机坐标系和世界坐标系重合，平移向量为0，将相机内参矩阵和外参矩阵进行化简得到下式：用最小二乘法计算三维点坐标，基于以下假设：空间任意点p在两个摄像机c1与c2上的图像点平p1与p2已经从两个图像中分别检测出来，即已知p1与p2为空间同一点p的对应点；摄像机c1与c2已标定，对应的投影矩阵分别为m1与m2，于是有：其中，(u1,v1,1)与(u2,v2,1)分别为p1与p2点在各自图像中的图像齐次坐标；(x,y,z,1)为p点在世界坐标系下的齐次坐标；m
ijk
为m
k
的第i行第j列元素，消去s1与s2得到关于x
w
，y
w
，z
w
的四个线性方程，即可求出三维点p的坐标。5.根据权利要求1所述的基于图像处理技术的破片飞散特性测试方法，其特征在于，根据拟合方程和已知的破片三维空间坐标求解破片运动参数，包括以下步骤：步骤s41，构造破片在空气中的运动微分方程：式(8)中，m为破片质量，c
x
为破片迎面阻力系数，ρ为当地空气密度，s为破片垂直于飞行方向上的迎风面积，t为破片飞行时间，v为破片瞬时速度；步骤s42，获取弹体静爆前坐标和爆炸时刻，结合获得的破片运动轨迹，采取间接解算的方法解算；首先，将轨迹的第一个时刻点作为相对时刻，其位置作为相对起始位置，则有：式中：t1为破片轨迹第一个数据点对应的时刻，x
′
为t1时刻该破片与爆炸点的距离，v1对应于t1时刻该破片的速度，v0对应于t0时刻该破片的速度，x为破片在某一时刻距离爆炸
点的距离，α为速度衰减系数，则有：式(10)中，x(t)为破片在任意t时刻飞行时对应飞行距离，等效为：式(11)中：s
tr
为破片相对t1时刻坐标的距离，t
r
为相对于t1时刻的时间差，利用获取得到的破片轨迹数据，即可拟合得到速度v1和速度衰减系数α；计算战斗部质心至上述破片飞行的拟合直线的垂足坐标，由于破片拟合直线不一定经过战斗部质心，因此以垂足坐标作为破片运动的起始点，在此基础上，计算轨迹第一个数据点与垂足间距离即为x
′
；将速度v1、速度衰减系数α、距离x
′
代入式(9)，解算得到初速v0。6.一种基于图像处理技术的破片飞散特性测试系统，其特征在于，通过权利要求1-5任一项所述的图像处理技术的破片飞散特性测试方法，实现基于图像处理技术的破片飞散特性测试。7.一种计算机设备，包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序，所述处理器执行所述计算机程序时，通过权利要求1-5任一项所述的图像处理技术的破片飞散特性测试方法，实现基于图像处理技术的破片飞散特性测试。8.一种计算机可读存储介质，其上存储有计算机程序，所述计算机程序被处理器执行时，通过权利要求1-5任一项所述的图像处理技术的破片飞散特性测试方法，实现基于图像处理技术的破片飞散特性测试。
技术总结本发明提出一种基于图像处理技术的破片飞散特性测试方法及系统，根据战斗部当量计算破片分布半径，按照破片分布半径安全区域布置双目相机，对相机进行标定与参数设置，采集破片飞行的轨迹序列图像；对采集到的轨迹序列图像进行破片的自动识别与坐标提取；对提取到的破片在两个相机视场中互相匹配，利用提取的破片在图像中的二维坐标计算破片在空间的三维坐标；在三维空间中对破片运动方程进行拟合，根据拟合方程和已知的破片三维空间坐标求解破片运动参数。本发明具有测量系统布局方便、实景记录、非接触式测量、测量精度高、分析速度快等特点，为新型目标毁伤计算与分析提供有效的理论依据。的理论依据。的理论依据。
技术研发人员：罗红娥 贺臻 孔筱芳 顾金良 夏言
受保护的技术使用者：南京理工大学
技术研发日：2022.06.10
技术公布日：2022/11/1
