基于弹丸图像的末敏弹阻力系数计算方法

胡志鹏，刘荣忠，郭 锐，朱光辉

（南京理工大学智能弹药技术国防重点学科实验室，江苏 南京 210094）

0 引言

末敏弹稳态扫描技术是末敏弹研制中的关键技术之一。双翼无伞末敏弹结合子弹本身的质量不对称及空气动力非对称形成稳态扫描运动，它具有体积小、落速高、受风影响小等优点，大大提高了末敏弹的作战效果，是目前世界各国竞相发展的新型灵巧弹药［1］。为分析双翼末敏弹的减速导旋特性，往往利用高塔投放的自由飞行试验方法，测试弹丸接近真实运动条件下的减速导旋性能。试验过程中，末敏弹的落速是其中一个重要指标参数，直接决定末敏弹的扫描效果，因此，作为实现末敏弹稳态扫描的基础性指标，获得无伞末敏弹的阻力系数［2］是试验研究过程中需要解决的重点问题之一。目前，弹丸阻力系数的求解方法主要包括试验研究［3］、数值仿真和理论计算［4］。由于末敏弹复杂的外形而使理论计算复杂化，仿真结果的准确性同样也需要试验结果的验证。进行风洞吹风试验可以得到精确的气动参数，但在初期研究中需要在多种模型设计中选取最优结构，选择风洞试验费用较高。传统弹道靶仪器准备时间长，事后处理底片复杂，需多次重复试验才能得到完整数据。将阻力系数的求解融进末敏弹的总体结构试验，节省人力物力成本的同时可以有效获得末敏弹飞行阻力系数。为此，本文提出以末敏弹弹道模型为基础，基于高速录像图像提取的末敏弹弹道参数辨识方法。

1 阻力系数辨识基本原理

1.1 末敏弹飞行运动学模型

试验求解末敏弹阻力系数过程中，将其简化为仅受空气阻力R和重力G作用下的单刚体运动。末敏弹在空气中的运动速度介于20～100 m／s之间，尾流阻力起主要作用的，此时阻力系数Cd可以视为一常数［5］。这时阻力与物体相对运动速度的二次方成正比，即

式（1）中，v为末敏弹飞行速度，S为末敏弹参考面积。末敏弹的运动方程可写为：

式（2）中，g是重力加速度，m为弹体质量。当阻力R大小与重力G相等时，末敏弹做匀速直线运动，此时末敏弹获得极限速度V：

［5］可得到极限速度表示的末敏弹弹道方程：

式（4）中，x为弹体下落实际高度，V为弹体极限速度，g为当地重力加速度，t为弹体飞行时间。

由式（4）可知，若已知末敏弹飞行时间序列值T［t1，t2，…，ti，…，tn］及 对 应 弹 道 高 度 向 量 X［x1，x2，…，xi，…，xn］，则可通过拟合的方法计算得到自由飞行试验中的末敏弹阻力系数值。

1.2 最小二乘拟合原理

假设有一组数据xi，yi，i＝1，2，…，N，且已知这组数据满足某一函数原型y^（x）＝f（a，x），其中a为待定系数向量，则最小二乘法［6］曲线拟合的目标就是求出这一组待定系数的值，使得目标函数

由末敏弹飞行时间序列值T［t1，t2，…，ti，…，tn］及对应弹道高度向量X［x1，x2，…，xi，…，xn］，通过合，可得到拟合的Cd初始估计值。

2 基于图像提取的试验数据处理

末敏弹自由飞行时间较短，天空背景几乎不变。通过对高速录像拍摄图像做差，可提取背景图像。通过对目标图像检测，得到末敏弹不同位置向量。由高速录像拍摄频率知末敏弹飞行时间序列，利用最小二乘拟合便得到末敏弹飞行时阻力系数。

2.1 运动图像提取方法

2.1.1 背景的提取

背景图像的提取主要有统计法、Surendra算法［7］等。Surendra算法能够自适应地获取背景图像，通过当前帧差分图像找到物体的运动区域，对运动区域内的背景保持不变，非运动区域的背景用当前帧进行替换更新，这样经过迭代可以提取出背景图像，Surendra算法如下：

将第1帧图像I0作为背景B0（x，y），迭代阈值τ，迭代次数m＝1，最大迭代次数Mmax。

当前帧的差分图像

由二值图像Dk更新背景图像Bk，即

式（7）中，Bk（x，y），Dk为背景图像和差分二值图像在（x，y）的灰度值；Ik为第k帧图像；α为迭代速度系数。迭代完成后的Bk（x，y）可视为背景图像。

2.1.2 运动目标的检测

在高速录像场景和摄像头之间的位置保持相对不变时，在摄取的序列图像里，背景图像的大小和位置近似不变。由于末敏弹的运动，相邻两帧的图像之间发生了变化。将相邻两帧图像的对应像素的灰度值相减得到的图像称为差分图像。在差分图像中各个点的像素值代表的是该位置上的像素在前后两帧中灰度值的变化。为了描述图像在前后两帧中的变化情况，用下述模型来描述具有运动物体的图像［8－9］：

式（8）中，fk（x，y）为第k帧图像，Bk（x，y）表示静止的背景，mk（x，y）表示其中的运动物体，nk（x，y）为背景噪声。将相邻两帧图像fk（x，y）、fk＋1（x，y）相减得到差分图像Dk，k＋1（x，y）。

式（9）中，s（x，y）是第k帧及k＋1帧中未变化的背景区域相减得到的图像（包括噪声），其像素特点是灰度值的绝对值很小，mdk（x，y）是第k帧图像中运动物体所占的区域。bdk（x，y）是物体经过运动在第k＋1中所占的区域，两者合并在一起就是由物体运动所引起的运动变化区域，其特点是灰度值的绝对值相对大一些，即该区域中像素的灰度值在前后两帧中的变化相对较为明显。由于噪声的影响，会使一些属于背景上的点被错误地检测为前景的运动目标，也会使前景目标点被错误地检测为背景点。同时由于背景上物体的轻微扰动也会使这些背景点被错误地识别为运动目标点。为了消除这些影响，需要对获得前景和背景的差分图像作一些处理，因此要选取一个阈值τ来减小噪声产生的影响并且加强运动变化产生的影响，由此得到二值的差值图像：

2.2 末敏弹运动过程离散数据

末敏弹自由飞行过程中，可认为天空背景不变。通过2.1.1节方法得到天空背景Bk（x，y）。进一步提取运动目标，序列选取末敏弹运动图像，与背景Bk（x，y）差分并去噪后得到末敏弹运动图像bk（x，y），进而得到运动目标的上下左右像素点极值。在本文中关心的是末敏弹在铅垂方向的运动，故可以由上下像素点确定的目标轮廓来确定其运动中心点Xm。

为得到末敏弹飞行的弹道实际高度，如图1所示，可在高速录像拍摄的背景中设置两个与地面高度确定的标尺。距离地面较近的标志物为A，距离地面实际距离为Ha，另一个记为B，其距离地面实际距离记为Hb。A的铅垂像素坐标为Xa，B的像素坐标为Xb。末敏弹的实际高度记为Hm，其在图像中的像素坐标为Xm。

由此可以计算末敏弹实际高度值

当高速录像以频率f拍摄处于自由飞行运动中的末敏弹，可计算在拍摄第n张图片至第n＋i张图片的末敏弹飞行时间为（（n＋i）－n）／f，设第一幅图像拍摄于末敏弹投出后的t1s，根据拍摄频率f可以确定第n张图像中末敏弹飞行时间为t1＋n／f。

图1 试场及标志物Fig.1 Field of view and makers

记末敏弹初始弹道高度为H0，故末敏弹在ti时间对应的飞行距离xi均可求。在若干图像中间隔抽取一定数量n张图像，依次求得末敏弹飞行时间向量T［t1，t2，…，ti，…，tn］及对应距离向量 X［x1，x2，…，xi，…，xn］，即为用于拟合的离散点。

2.3 最小二乘拟合法提取阻力系数

即取末敏弹飞行到标志A与B的中点速度v－及坐标，代入式（4）后得到V，即求得Cd。

初始时间t1为未知，因此可以通过逐渐增大试取t1的值，由拟合误差来判断是否达到要求。当误差小于设定值时，认为所取t1与试验中末敏弹飞行到图像拍摄时刻吻合。得到V值后将其做为初值继续重复前面过程，直到V值变化比规定值小为止，通常取V的1%。

3 误差分析与试验验证

3.1 测试环境及条件

试验末敏弹从高塔顶部100 m位置点投放，初始速度为零。试验用Red lake dk100高速录像机，距弹体投放点在地面投影中心为100 m距离处以1 000幅／s的频率拍摄，生成图像像素为640×480。标志物为铅垂放置且有固定尺寸刻度的绳索，如图1所示，通过测量底部标志物A点的实际距离确定各标志物的高度。为准确记录弹体下落到地面时间，在末敏弹上装有记录弹体运动特征的黑匣子。

3.2 误差分析

试验结果误差的主要来源有弹体下落点散布引起的随机误差和读取图像中弹体中心位置与实际弹体质心位置不重合而引起的误差。

由弹体下落点散布引起的误差如图2所示，由图可见单次测量弹体下落位置散布引起的随机误差极限值，根据三角形相似原理，表达为下式：

式（12）表明，弹体下落位置散布引起的测量误差与视场高度H、高速录像与标志物距离x1和弹体在高速录像与标志物连线上的距离分布x2有关，且绝对误差与x2和H成正比，与x1与反比。

图2 弹体落点散布引起误差Fig.2 The err or of TSP i mpact dispersion

由图像提取弹体质心位置是误差的另一主要来源，理论上提取点像素值应为弹体质心位置，实际由于弹体姿态不断变化使图像提取过程中得到的坐标与弹体实际质心位置不重合。设弹体在图像中铅直方向像素数为P，则这一误差最大值为1／2P。

由式（12）的关系及误差传递理论可知测量中的随机误差E的极限值为［10］

3.3 试验验证

对结果的验证将从三个方面进行分析，首先利用计算流体力学方法对所用末敏弹的气动外形进行数值计算，得到末敏弹阻力系数，与飞行运动过程所测数据的拟合得到阻力系数进行对比，结果见表1。

由表1可见，试验计算的阻力系数与仿真结果有较大误差，大部分仿真计算结果要较拟合结果大，源于数值计算也有一定的不确定度。为进一步验证结果的准确性，同时直观地进行对比，将拟合结果得到的阻力系数代入末敏弹三自由度运动方程，得到这一阻力系数下的末敏弹弹道轨迹，如图3所示。

表1 试验结果与仿真计算阻力系数Tab.1 Test C d and si mulation C d result

图3 弹丸飞行轨迹与理论弹道曲线图Fig.3 Trajectory of the TSP and theory ballistic

由图3可见末敏弹试验飞行轨迹与理论弹道曲线吻合较好。

弹道程序所计算的末敏弹由初始高度飞行至地面所用时间，称为理论飞行时间。而弹体所载黑匣子记录的时间为标准，两者进行对照，结果列于表2。

表2 理论计算与实测飞行时间误差Tab.2 The ti me error bet ween co mpute and test

由表2中测试结果与计算结果的对比可以发现误差最大值在2%左右，证明计算结果比较准确地反应了弹体阻力系数。

分别计算末敏弹经过两相邻标志物间的平均速度，并由仿真计算的阻力系数计算理论速度曲线，如图4所示。

图4 实测平均速度与理论速度曲线Fig.4 Test mean velocity and theory velocity

由图4可见，末敏弹经过各标志物点平均速度分布与理论计算结果较为接近，速度变化趋势亦相符。图4为理论速度曲线与试验速度值对比，由平均速度法测得的速度和理论速度曲线同样很接近。

取标志物A与B中点，通过AB距离及高速录像图片间隔求取末敏弹通过这一点时的速度，作为试验中末敏弹的平均速度，而由拟合结果得到的阻力系数计算得到的通过这点的速度为理论速度，结果如表3所示。

表3 理论计算速度与平均速度结果Tab.3 The t heoretical velocity and mean velocity of test

结合图4曲线与表3的数据可知，平均速度和根据计算结果而推算的理论速度值符合良好，进一步说明此种方法的可行性与可靠性。

试验中由于弹体的落点随机分布引起测量误差以及图像提取弹体质心位置时也会产生误差，这些差误在处理过程中不会产生较大影响而且可以通过一定的措施来减少误差的影响，通过对结果的比较可以看到，在处理试验末敏弹阻力系数过程中这种方法还是很可靠的，与实际测量结果吻合。

4 结论

本文提出了基于弹丸图像的末敏弹阻力系数计算方法，通过Surendra方法目标背景图像，并对高速录像拍摄图像进行运动目标提取，获得运动目标的时间、飞行弹道向量，结合末敏弹运动模型，利用最小二乘拟合方法求解自由飞行状态下末敏弹阻力系数。试验结果表明：在不影响末敏弹高塔自由飞行试验中其它气动特性测量的情况下，能够快速处理高速录像获得的大量数据。基于高速录像的图像处理方法计算阻力系数结果准确，与实际测得结果偏差较小，最大测量偏差在2%以内。结合末敏弹姿态测量结果和空气动力系数数值模拟结果可以进一步提高这一方法的精度，从而为无伞末敏弹设计提供帮助。

参考文献：

5）了解不同国家地区低硫油加注监管要求。不同国家对低硫燃油的监管有不同要求，船东要提前了解船舶挂靠港口海事方面的监管要求，掌握港口检查内容和要求，规范船舶燃油使用记录，在航海日志等法定记录上做好低硫燃油采购的时间、地点、供应商等记录，保存好相应的单据和检测化验报告，便于主管机关随时查阅，也便于事后追溯。

［1］史金光，韩子鹏，舒敬荣，等.双翼无伞末敏弹稳态扫描运动数学模型［J］.弹道学报，2010，22（2）：24-27.SHI Jinguang，HAN Zipeng，SHU Jingrong，et al.Stable scanning model f or non-parachute ter minal sensitive projectile with t wo fins［J］.Jour nal of Ballistics，2010，22（2）：24-27.

［2］刘世平，韩子鹏，余劲天，等.末敏弹伞系统弹道模拟试验数据的处理［J］.弹道学报，1997，9（3）：70-74.LIU Shiping，HAN Zipeng，YU Jintian，et al.Parameters estimation and trajectory calculation of parachute-sub munition system fro m the data of free flight test［J］.Jour nal of Ballistics，1997，9（3）：70-74.

［3］First Austin Howard，Donald Elger.Rotary wing decelerators f or miniature at mospheric entr y probes［C］／／20th AIAA Aerodynamic Decelerator Systems Technology Conference and Seminar.Seattle，Washington，2009：2996.

［4］Peter Crimi.Finite element analysis of a samara-wing decelerator［J］.Jour nal of Aircraft，1996，33（4）：793-802

［5］张延教.高等动力学［M］.南京：南京理工大学出版社，2003.

参考文献：

［6］万建伟，熊栋梁，皇甫堪.弹丸阻力系数CD的Robust估计［J］.兵工学报，1991（1）：21-28.WAN Jian wei，XIONG Dongliang，HUANG Fukan.Robust Estimation for the Coefficient of Resistance CD［J］.Acta Ar mamentarii，1991（1）：21-28.

［7］王正勤，刘富强.自适应背景提取算法的比较［J］.计算机工程，2008，34（23）：220-223.WANG Zhengqin，LIU Fuqiang.Comparison of adaptive extract algorith m of object scene［J］.Computer Engineering，2008，34（23）：220-223.

［8］范莹，郭成安.一种运动图像的检测与识别技术［J］.大连理工大学学报，2004，44（1）：122-126.FAN Ying，GUO Chengan.A method f or detection and recognition of moving objects in sequential i mages［J］.Jour nal of Dalian University of Technology，2004，44（1）：122-126.

［9］刘贵喜，邵明礼，刘先红，等.真实场景下视频运动目标自动提取方法［J］.光学学报，2006，26（8）：1 150-1 155.LIU Guixi，SHAO Mingli，LIU Xianhong，et al.Video moving object auto-extraction in real scene［J］.Acta Optica Sinica，2006，26（8）：1 150-1 155.

［10］蒋恩臣，蒋亦元，刘道顺.高速摄影拍摄频率与判读的误差分析［J］.东北农业大学学报，2000，31（4）：381-384.JIANG Enchen，JIANG Yiyuan，LIU Daoshun.Studies on surveying error and photo frequency of ultra-speed cinematography［J］.Journal of Northeast Agricultural University，2000，31（4）：381-384.