一种火箭橇图像测速方法

摘要

一种火箭橇图像测速方法，采用高速相机连续拍摄橇体的运行过程，在试验完成后通过对高速相机相对轨道的俯仰和方位信息修正、橇体特征点提取、运行速度计算，完成对橇体运行速度的测量。该方法具有测速精度高、思路简单明了，编程容易实现的优点，提高了火箭橇试验中关键位置的速度测量精度。

说明书

技术领域

本发明涉及一种火箭橇图像测速方法，尤其涉及一种采用高速照相机连续拍摄来测量橇体运行速度的方法，属于火箭橇测试领域。

背景技术

火箭橇试验具有产生大过载、高速度、强振动和冲击等综合条件的能力，可以最逼真地模拟导弹真实飞行环境。通过试验能够考核惯性测量系统在综合环境条件下的各项性能指标和精度，验证惯性测量系统误差模型在高动态条件下的正确性，特别是在大过载情况下，高次项放大作用，能够确定惯性测量系统高次误差项对导航性能的影响，是实现惯性测量系统动态性能验证的最佳途径。

在惯性测量系统火箭橇试验中，目前主要采用雷达测量设备、遮光板时空测量装置测量橇体的运行位置和速度。

雷达测量设备由雷达天线向火箭橇发射出连续等幅的电磁波并接收从橇体反射回来的电磁波信号，利用测量的多普勒频率，经处理获得橇体橇体相对于天线指向的径向速度以及橇体沿轨道的运行速度。但是，雷达测速存在以下两个缺点：

(1)适合于高速测量，低速测量精度低。在火箭橇试验的实际应用中在速度低于10m/s时不输出测量值；

(2)火箭橇固体发动机的尾焰效应使雷达测量的精度较低，在火箭橇一级发动机和二级发动机点火时会引起雷达测速的数据实效。

遮光板时空测量装置由橇载光电组件和沿火箭橇滑轨一侧每间隔一定距离安装的若干遮光板组成。光电组件安装在橇体上，由激光器和激光接收器组成，激光接收器连接有处理电路，实现把遮光板引起光通量的变化转变为电脉冲信号，并输出给数据采集系统。数据采集系统则记录该电脉冲的瞬时时刻，通过相邻遮光板沿轨道的距离可求出橇体的运行速度。但是，采用遮光板光电组件的缺点是：

(1)该光电组件的激光器对固体发动机点火后的光线较敏感，容易导致处理电路处于饱和状态，从而引起光电组件工作异常；

(2)相邻遮光板间距一般大于1米，采用相邻遮光板距离除以测量时间间隔求取橇体运行的速度的方法存在测量误差相对较大的缺点。

通过以上分析可以看出，随着火箭橇轨道后续由9km延长到16km，采用多级点火的需求越来越迫切，雷达测速和遮光板时空测量已不能适应多级点火段的要求。因此，需要研究一种新的测速方法以弥补上述两种方法的缺陷。

在火箭橇试验场，目前已有高速照像，但只用于观察发动机点火瞬时多级橇体分离是否正常等功能监测用途，本发明专利是研究一种在现有高速照像的基础上实现测速的新方法。

发明内容

本发明解决的技术问题是：找到了一种火箭橇图像测速方法，该方法能够精确的测量橇体运行的速度，解决多级发动机分离时雷达、遮光板受发动机点火影响而失效的问题，以提供高精度的速度信息。

本发明解决的技术方案是：

一种火箭橇图像测速方法，步骤如下：

(1)在火箭橇试验中，在多级发动机点火分离处的轨道旁架设高速相机，设置相机的拍摄频率f、原始照片中纵向边长的像素个数M₀和横向边长的像素个数N₀；测量相机相对轨道的地理信息和轨道上相邻扣件之间的距离L；M₀、N₀均为正整数，f、L为正实数；

(2)在火箭橇试验多级发动机点火分离的过程中，通过高速相机连续拍摄K幅照片，K为正整数；在每幅照片中由两维像素点附加颜色的方式构成，像素点在照片中的位置信息为原始照片像素点坐标系的坐标值；

(3)在每一幅照片中，将原始照片像素点坐标系中的坐标点(m,n)转换到原始照片中心点坐标系中的点(x,y)，其转换公式如下：

经转换后，用(x,y,r_x,y,g_x,y,b_x,y)描述原始照片中该像素点的位置和颜色信息；

所述原始照片中心点坐标系以照片中心为原点、ox轴水平向右的右手坐标系oxy；

(4)根据在多级发动机点火前的第一幅照片中轨道畸变信息修正照片，使得照片中每个像素所对应的空间距离相等；

(5)根据步骤(4)中得到的修正后无畸变的照片，提取每一幅照片中火箭橇橇体的第一特征点在修正后照片像素点坐标系om’n’中的位置(I_s,J_s)和该第一特征点在此位置的像素值s代表照片幅数；

(6)根据步骤(5)得到的第一特征点水平方向像素值J_s，以及步骤(1)中的参数，求取火箭橇橇体的速度序列；

(7)根据步骤(4)中得到的修正后无畸变的照片，提取每一幅照片中火箭撬橇体的其它特征点，重复步骤(5)和(6)，得到每个特征点对应的火箭橇橇体的速度序列，之后与第一特征点对应的速度序列进行融合，从而得到火箭橇橇体运行的全过程速度序列。

原始照片像素点坐标系omn的原点在照片左上角，om轴朝下，on轴朝右；像素点的颜色值由红色、绿色和蓝色不同比值构成，用五维数组(m,n,r_m,n,g_m,n,b_m,n)表示，其中，(m,n)描述了该像素点在照片中的位置，m的取值范围为从1至M₀的正整数，n的取值范围为从1至N₀的正整数；(r_m,n,g_m,n,b_m,n)描述了原始照片中该像素点的颜色信息，红色r_m,n、绿色g_m,n、蓝色b_m,n的取值范围从0至255。

所述步骤(4)根据在多级发动机点火前的第一幅照片中轨道畸变信息修正照片，具体为：

(3.1)在修正处理前的照片中，从左到右依次提取每个扣件的坐标值并进行编号，设有h≥4个扣件，以坐标值(x_p,y_p)表示，其中p＝1,2,…,h；

(3.2)对照片中的扣件序列(x_p,y_p)，p＝1,2,…,h；采用直线模型y_l＝k_lx+b_l，用最小二乘法计算系数k_l和b_l的值，具体计算公式如下：

其中，

(3.3)对照片中的扣件序列(x_p,y_p)，p＝1,2,…,h；计算相邻两个扣件的距离l_q和中心位置坐标(x_0,q,y_0,q)，q＝1,2,…,h-1；具体计算公式为

(3.4)取直线轨道的二次畸变模型采用最小二乘法计算系数α和β，具体公式如下：

其中

(3.5)根据系数α、k_l和b_l的值，设置M₁和N₁，并修正原始照片，具体公式如下：

(r_m',n',g_m',n',b_m',n')＝(r_m,n,g_m,n,b_m,n)

其中,数组(r_m,n,g_m,n,b_m,n)的下标m,n表示了原始照片像素点坐标系的坐标点位置；数组(r_m',n',g_m',n',b_m',n')的下标m’,n’表示了修正后照片像素点坐标系的坐标点位置，矫正后照片像素点坐标系om’n’是原始照片矫正后，以照片左上角为原点，om’轴垂直向下，on’轴水平向右的坐标系；M₁和N₁分别为矫正后照片的纵向和横向像素个数，令m′和n′取值大于0；(x’,y’)表示矫正后照片中心点坐标系ox’y’的坐标点，矫正后照片中心点坐标系ox’y’是原始照片矫正后，以照片中心为原点、ox’轴水平向右的右手坐标系；

(3.6)修正后的照片中，被拉伸的区域中存在没有被上色的点，在修正后照片像素点坐标系中对被拉伸区域进行修补，具体为：

其中，(r_m″,n″,g_m″,n″,b_m″,n″)表示修正后照片中没有上色点的颜色信息，m″和n″表示在修正后照片像素点坐标系没有上色点的坐标位置，λ表示该像素点周围八个像素点中，红色、绿色、蓝色的先像素值之和大于0的个数。

所述步骤(5)中提取每一幅照片中火箭橇橇体的第一特征点在修正后照片像素点坐标系om’n’中的位置(I_s,J_s)，和该第一特征点在此位置的像素值具体为：

(4.1)从第一幅照片开始，等间隔选取Z幅照片，找到这些照片在修正后照片像素点坐标系om’n’的橇体特征点坐标(I_s,J_s)；同时提取橇体的第一特征点在此位置的像素值其中代表坐标点(I_s,J_s)红色像素值，代表坐标点(I_s,J_s)绿色像素值，代表坐标点(I_s,J_s)蓝色像素值，s代表图片幅数；

找到这些等间隔照片的第一特征点，读取这些第一特征点的五维数组根据以上五维数组所在的范围确定橇体第一特征点颜色红色、绿色和蓝色的域值范围分别为(a₁,a₄),(a₂,a₅),(a₃,a₆)，以及撬体特征点坐标(I_s,J_s)所在的范围I_min≤I_s≤I_max,J_min≤J_s≤J_max；a₁、a₂、a₃、a₄、a₅、a₆均为正整数；

(4.2)从第一幅照片开始，读取每一副照片信息，依据I_min≤I_s≤I_max,J_min≤J_s≤J_max范围内坐标点(I_s,J_s)的颜色数据确定特征点的位置；若同时满足和并且J_s的值为最大时，则提取(I_s,J_s)，并且记录幅数s；若不能满足或中任意一个，则取J_s＝J_s-1。

所述步骤(6)根据水平方向像素值J_s，以及步骤(1)中的参数，求取火箭橇橇体的速度序列，具体步骤为：

(5.1)给定x₁(0)和x₂(0)的初始值0；

(5.2)根据照片中特征点位置的变化计算第s副照片中瞬时速度v_s，1≤s≤K–1，具体公式为：

v_s＝kfc(J_s+1-J_s)

式中c为单个像素点对应的长度，f为相机拍摄的频率，J_s+1和J_s分别表示第s+1张照片和第s张照片中撬体特征点的水平方向坐标，k为比例系数k表示单个像素点对应于实际空间中的距离与其成像元器件上实际长度的比值，v_s为瞬时速度；

(5.3)将所有照片的特征点的瞬时速度v_s按照时间顺序依次排列形成瞬时速度序列，对瞬时速度序列进行平滑处理后得到橇体的速度序列其中，速度的计算公式为

其中，b₀、A₁、A₂、B₁、B₂为离散化二阶低通滤波器的参数，x₁(0)和x₂(0)为离散化二阶低通滤波器的两个初始参数。

所述步骤(7)得到每个特征点对应的火箭橇橇体的速度序列，之后与第一特征点对应的速度序列进行融合，具体为：

(6.1)根据第一幅照片确定其他特征点和第一特征点位置关系，即其他特征点水平方向像素值J_s^τ和第一特征点的水平方向像素值J_s之间关系：

J_s^τ＝J_s+ω^τ

其中ω^τ是其他特征点和第一特征点水平方向像素值的差值，τ＝2，3，4…,J_s^τ表示第τ特征点水平方向像素值；

(6.2)重复步骤(5)得到其他特征点的速度序列，将所述其他特征点速度序列和第一特征点速度序列按时间先后顺序排列，若同一时刻存在多个速度，取这些同一时刻速度的算数平均值，由此完成融合，得到火箭橇橇体运行的全过程速度序列。

本发明与现有技术相比的优点如下：

(1)本发明能够精确得到火箭橇多级发动机点火分离时的橇体运行速度，克服雷达测速、遮光板光电组件测量不准确的问题，测量精度较高。

(2)这种方法可以半全自动运行，不需要加入过多的人工参与，能够实现照片俯仰和方位偏差信息的自动修正、发动机点火前后橇体特征点提取、橇体速度计算和噪声消除，理论清晰简单，编程容易。

(3)目前尚没有能够通过高速相机对橇体运行速度进行精确测量的方法，本方法是一种创新，为分析和评价火箭橇试验的设计合理性提供了一种处理办法。

(4)本发明将照片转换为数值矩阵，利用数值关系分析特征点位置变化，整个过程科学合理、严谨准确，同时能更快速精确的读取照片数据进行分析。同时提取照片特征点的方法，所需存取的数据量很小，即使处理大批量连续拍摄的照片也能在很短的时间分析出特征点位置变化

(5)目前使用的修正照片方法都存在不尽人意的地方，需要一定的工作环境、速度慢、效率低且价格昂贵。本发明对原始照片建立畸变模型直接求取运动方向的实际变化，有其明显的优势，大大提高修正效率。

(6)本发明采用速度滤波平滑的方法，不仅大幅消除速度噪声的影响，而且使得的速度曲线更加合理、清晰，同时使用多个特征点分析火箭橇撬体速度变化，满足了在环境遮挡、光线变化、多级点火等苛刻外界条件也能精准测量的要求。

附图说明

图1为本发明实现流程图；

图2为火箭橇轨道和扣件在原始照片中位置。

图3为扣件序列直线模型；

图4为二次畸变模型示意图；

图5为矫正后的照片；

图6为橇体特征点I在照片中的位置；

图7为橇体特征点II在照片中的位置；

图8为橇体特征点I在不同照片中的像素点位置；

图9为橇体特征点II在不同照片中的像素点位置；

图10为橇体特征点I和特征点II经对齐后在不同照片中的像素点位置；

图11为经融合后橇体特征点I在不同照片中的像素点位置；

图12为对速度噪声的频谱分析结果；

图13为橇体运行速度的计算值和平滑值。

具体实施方式

本发明的实现原理是：采用高速相机连续拍摄橇体的运行过程，在试验完成后通过对高速相机相对轨道的俯仰和方位信息修正、橇体特征点提取、运行速度计算，完成对橇体运行速度的测量。该方法具有测速精度高、思路简单明了，编程容易实现的优点，提高了火箭橇试验中关键位置的速度测量精度。

下面给出一种火箭橇图像测速方法，其特征在于步骤如下：

(1)在火箭橇试验多级发动机点火分离处的轨道旁架设高速相机，设置相机的拍摄频率f、原始照片中纵向边长的像素个数M₀和横向边长的像素个数N₀等参数；测量相机相对轨道的地理信息和轨道上相邻扣件之间的距离L；M₀、N₀均为正整数，f、L为正实数；如图2所示，在照片底部火箭橇轨道上标记为紫色矩形方框的部分为扣件在火箭橇轨道的位置，白色方框内数值表示两个相邻扣件的像素点距离。

(2)在火箭橇试验多级发动机点火分离的过程中，通过高速相机连续拍摄K幅照片；在每幅照片中由两维像素点附加颜色的方式构成，像素点在照片中的位置信息为原始照片像素点坐标系的坐标值；原始照片像素点坐标系(omn)的原点在左上角，om轴朝下，on轴朝右；像素点的颜色值由红色、绿色和蓝色不同比值构成；用五维数组(m,n,r_m,n,g_m,n,b_m,n)来表示，其中，(m,n)描述了该像素点在照片中的位置，m的取值范围为从1至M₀的正整数，n的取值范围为从1至N₀的正整数；(r_m,n,g_m,n,b_m,n)描述了原始照片中该像素点的颜色信息，红色r_m,n、绿色g_m,n、蓝色b_m,n的取值范围从0至255。

(3)在每一幅照片中，将原始照片像素点坐标系中的坐标点(m,n)转换到以照片中心为原点、ox轴水平向右的右手坐标系oxy(原始照片中心点坐标系)中的点(x,y)，其转换公式如下：

经转换后，用(x,y,r_x,y,g_x,y,b_x,y)描述了原始照片中该像素点的位置和颜色信息；如图3所示将原始照片经过去色处理后，标记网格刻度线，原始照片中心点坐标系横轴和纵轴为图中交叉黑色粗线：

(4)根据在发动机点火前的第一幅照片中轨道畸变信息修正照片，使得照片中每个像素所对应的空间距离相等；修正照片的具体步骤为：

a.从左到右依次提取每个扣件的坐标值并进行编号，设有h(h≥4)个扣件，以坐标值(x_p,y_p)表示，其中p＝1,2,…,h；如图2所示，在照片底部火箭橇轨道上标记为紫色矩形方框的部分为扣件在火箭橇轨道的位置，对其逐个进行标记其坐标值(xp,yp)；

b.对照片中的扣件序列(x_p,y_p)，p＝1,2,…,h；采用直线模型y_l＝k_lx+b_l，

用最小二乘法计算系数k_l和b_l的值，具体计算公式如下：

其中，

c.对照片中的扣件序列(x_p,y_p)，p＝1,2,…,h；计算相邻两个扣件的距离l_q和中心位置坐标(x_0,q,y_0,q)，q＝1,2,…,h-1；具体计算公式为

d.取直线轨道的二次畸变模型采用最小二乘法计算系数α和β，具体公式如下：

其中

e.根据系数α、k_l和b_l的值，设置M₁和N₁，和矫正原始照片，

具体公式如下：

(r_m',n',g_m',n',b_m',n')＝(r_m,n,g_m,n,b_m,n)

其中数组(r_m,n,g_m,n,b_m,n)，下标m,n表示了原始照片像素点坐标系的坐标点位置；数组(r_m',n',g_m',n',b_m',n')，下标m’,n’表示了矫正后照片像素点坐标系的坐标点位置，矫正后照片像素点坐标系(om’n’)是原始照片矫正后，以照片左上角为原点，om’轴垂直向下，on’轴水平向右的坐标系；M₁和N₁分别为矫正后照片的纵向和横向像素个数，取值应当保证m′和n′大于0；(x’,y’)表示了矫正后照片中心点坐标系(ox’y’)的坐标点，矫正后照片中心点坐标系(ox’y’)是原始照片矫正后，以照片中心为原点、ox’轴水平向右的右手坐标系。

f.被拉伸的区域中存在没有被上色的点，在矫正后照片像素点坐标系中矫正后的照片修补方法为：

其中，(r_m″,n″,g_m″,n″,b_m″,n″)表示修正后照片中没有上色点的颜色信息，m″和n″表示在修正后照片像素点坐标系没有上色点的坐标位置，λ表示该像素点周围八个像素点中，红色、绿色、蓝色的先像素值之和大于0的个数。如图5所示经过修正后的照片，其中每个像素点大小所代表的空间实际距离相等，照片的修补部分为照片的外周围黑色部分；

(5)根据步骤(4)中得到的修正后无畸变的照片，提取每一幅照片中火箭橇橇体的第一特征点在修正后照片像素点坐标系om’n’中的位置(I_s,J_s)和该第一特征点在此位置的像素值s代表照片幅数；如图6所示，第一特征点在火箭橇的黑色中线部分最右端；

(6)根据步骤(5)得到的水平方向像素值，以及步骤(1)中的参数，求取橇体的瞬时速度，具体步骤为：

a.给定x₁(0)和x₂(0)的初始值0；

b.根据照片中特征点位置的变化计算第s副照片中瞬时速度v_s，1≤s≤K–1，具体公式为：

v_s＝kfc(J_s+1-J_s)

式中c为单个像素点对应的长度，f为相机拍摄的频率，J_s+1和J_s分别表示第s+1张照片和第s张照片中撬体特征点的水平方向坐标，k为比例系数k表示单个像素点对应于实际空间中的距离与其成像元器件上实际长度的比值，v_s为瞬时速度；

c.将所有照片的特征点的瞬时速度v_s按照时间顺序依次排列形成瞬时速度序列，对瞬时速度序列进行平滑处理后得到橇体的速度序列其中，速度的计算公式为

其中，b₀、A₁、A₂、B₁、B₂为离散化二阶低通滤波器的参数，x₁(0)和x₂(0)为离散化二阶低通滤波器的两个初始参数。如图12为对速度噪声的频谱分析结果，可以看出，噪声主要频段为大于10Hz的高频段；因此，采用二阶低通滤波器对速度进行滤波以消除高频噪声

(7)根据步骤(4)中得到的矫正后无畸变的照片，提取每一幅照片中橇体的其它特征点，重复步骤(5)和(6)，得到速度后与利用特征点I得到的速度进行融合，从而得到橇体运行的全过程速度序列。如图13为橇体运行速度的计算值和平滑值，同时平滑值清晰地表明火箭橇的全过程运动速度变化。

所述步骤(5)中提取特征点的实现方法：

(1)从第一幅照片开始，等间隔选取Z幅照片，找到这些照片在矫正后照片像素点坐标系(om’n’)的橇体特征点坐标(I_s,J_s)；同时提取橇体的特征点在此位置的像素值其中代表坐标点(I_s,J_s)红色像素值，代表坐标点(I_s,J_s)绿色像素值，代表坐标点(I_s,J_s)蓝色像素值，s代表图片幅数。根据这些等间隔照片确定撬体特征点颜色域值范围分别为(a₁,a₄),(a₂,a₅),(a₃,a₆)，以及撬体特征点坐标(I_s,J_s)所在的范围I_min≤I_s≤I_max,J_min≤J_s≤J_max。

(2)从第一幅照片开始，读取每一副照片信息，依据I_min≤I_s≤I_max,J_min≤J_s≤J_max范围内坐标点(I_s,J_s)的颜色数据确定特征点的位置；若同时满足且J_s的值为最大时，则提取(I_s,J_s)，并且记录幅数s；若不能满足任意一个，则取J_s＝J_s-1。如图8为特征点I在各幅照片中的像素点位置，其中像素点位置是以矫正后照片像素点坐标系中的坐标取出，曲线中有不平滑的部分主要由于环境噪声和光线影响，取J_s＝J_s-1，但是并不影响曲线总体趋势和精度；如图11为经融合后，橇体特征点在各幅照片中的像素点位置，其中像素点位置是以矫正后照片像素点坐标系中的坐标取出，特征点II通过平移后，将所述特征点II像素点位置和特征点I像素点位置按时间先后顺序排列，若同一时刻存在多个像素点位置，取这些同一时刻像素点位置的算数平均值。

通过上述方法即可完成火箭橇橇体运行速度的测量。

实施例：

在一次火箭橇试验中，设拍摄频率f＝2000幅/秒、单个像素点对应的长度c＝22μm、比例系数k为＝853.65。图2为火箭橇轨道和扣件在原始照片中的位置，并且通过图像处理将照片处理为无背景噪声，在照片底部火箭橇轨道上标记为紫色矩形方框的部分为扣件在火箭橇轨道的位置，白色方框内数值表示两个相邻扣件的像素点距离。；图3为扣件序列在原始照片的点位置建立的直线模型，坐标系为以照片中心为原点、ox轴水平向右的右手坐标系oxy(原始照片中心点坐标系)，坐标轴用照片中间交叉黑色粗线表示，建立的直线为y＝-0.01057x-101.2；图4为二次畸变模型示意图表示像素点位置和相邻扣件距离导数的关系，图中空心圆圈代表的是已知的相邻扣件距离导数和像素点位置的关系，黑色曲线为依照空心圆圈位置形成的拟合曲线，建立的畸变模型为图5为矫正后的照片，其中每个像素点所对应的空间实际距离相等，并且修补了被拉伸的区域中存在没有被上色的点的部分；图6为特征点I在橇体上的位置，此特征点为火箭橇橇体上中间黑色线状部分的最右端，与周围像素点颜色值有很大的差值，易于确定范围取出其坐标，满足读取特征点的操作要求，但是，特征点I会在图7照片中被遮挡，因此，选取特征点II，此特征点为火箭橇橇体上黑色三角锥状顶部外壳白色中心点，其在橇体上的位置见图7；图8为特征点I在各幅照片中的像素点位置，其中像素点位置是以矫正后照片像素点坐标系中的坐标取出，曲线中有不平滑的部分主要由于环境噪声和光线影响，但是并不影响曲线总体趋势和精度，图9为特征点II在各幅照片中的像素点位置，其中像素点位置是以矫正后照片像素点坐标系中的坐标取出，并且从第70幅到第140幅与特征点I的位置通过平移后可以进行比较验证；图10为对特征点I和特征点II像素点位置经对齐后在各幅照片中的像素点位置；如图11为经融合后，橇体特征点在各幅照片中的像素点位置，其中像素点位置是以矫正后照片像素点坐标系中的坐标取出，特征点II通过平移后，将所述特征点II像素点位置和特征点I像素点位置按时间先后顺序排列，若同一时刻存在多个像素点位置，取这些同一时刻像素点位置的算数平均值；图12为对速度噪声的频谱分析结果，可以看出，噪声主要频段为大于10Hz的高频段；因此，采用二阶低通滤波器对速度进行滤波以消除高频噪声，图13为橇体运行速度的计算值和平滑值，可以看出计算值有等幅的波动，保持稳定趋势，同时平滑值清晰地表明火箭橇的运动速度变化趋势，至此完成火箭橇图像测速全过程。

本发明未详细描述内容为本领域技术人员公知技术。