随着无人机技术的迅速发展, 无人机的应用领域也越来越广. 目前, 高效的飞控算法, 结合陀螺仪、超声定位、GPS定位和计算机视觉算法, 使无人机能够应对较复杂的自然环境. 但是某些特殊情况下, 如室内环境, GPS失效, 无人机精确定位只能依靠陀螺仪和机载摄像头, 陀螺仪仅能保持无人机自身姿态, 无法解决无人机位置和朝向漂移的问题. 因此, 这种情况下, 无人机的精确定位就只能依靠计算机视觉算法. 但是, 计算机视觉定位严重依赖于图像质量, 在弱纹理情况下, 匹配误差和置信度较差, 使无人机的姿态参数较难估计, 对无人机的控制方面带来了一定的难度. 因此, 室内弱纹理环境下无人机参数的估计非常有必要, 准确的参数估计可以更好的控制无人机, 使应用场景更多. 图1是无人机参数估计流程图.

为了更好的提升无人机的稳定性, 国内外多名学者进行了研究. 王亭亭等^[1]采用基于原始ORB特征的LK光流法, 利用随机抽样一致滤波(RANSAC)来提高光流场的精度, 再融合里程计、惯导等信息, 实现了无人机室内导航, 拥有较高的实时性, 但其光流速度估计处理较为简单, 采用传统的前后单追踪策略, 滤除误匹配效果不理想; 齐咏生等^[2]针对无人机在停机坪上空飞行时, 采集地面视频, 视频中包含停机坪上的几何特征, 预处理后, 采用五步标志提取方法与距离三点法的角点检测算法相结合, 对视觉标志进行提取, 提取过后进行匹配, 得到单应性矩阵, 分解则可得到飞行的参数, 该方法针对场景是含有H标志的停机坪, 应用场合具有较明显的特征, 在室内面对弱纹理区域及重复纹理区域时, 没有明显的几何标志, 匹配过程中该方法错误率较高. 周克旻等^[3]使用目前流行的基于增强现实(augmented reality, AR)技术的ArUco标记算法和颜色空间域标记算法, 设计了一种多标记的无人机跟踪系统, 对无人机在室内飞行过程中进行参数估计.

针对以上部分文献存在的问题, 本文使用的算法在面对室内弱纹理或者重复纹理区域时, 可以准确的估计出无人机飞行时的平移参数及旋转参数. 本方法在面对室内弱纹理区域时, 使用特征提取较快的ORB光流法, 配合前后双向双追踪策略和RANSAC算法, 提高参数估计的准确性, 降低了误匹配概率; 在旋转参数估计时, 提取单应性矩阵, 分解进行估计, 得到准确的参数.

1 本文算法

本文主要针对的是室内弱纹理及重复纹理区域采集的视频, 所以算法主要分成以下3部分: 视频预处理、平移参数估计和旋转参数估计, 流程图如图2所示. 算法实现步骤如下.

步骤1. 预处理阶段: 在处理视频帧之前, 首先将彩色视频帧转换成灰度图, 进行后续处理, 节省时间. 拍摄过程中难免存在一定的噪声, 为了不对后续处理造成影响, 先对采集的视频帧进行去噪处理.

步骤2. 平移参数估计阶段: 根据室内弱纹理区域可能存在其他物体的情况, 对采集视频图像进行判断, 根据不同环境采用不同的算法进行估计.

步骤3. 旋转参数估计阶段: 通过ORB算法, 提取出4个特征点, 如地板的角点等, 结合OpenCV图像库可以很好的计算出一个拟合所有点的单应性矩阵, 分解得到旋转分量

2 预处理

无人机室内飞行过程中, 面对地板砖等弱纹理区域, 会采集到该区域的彩色视频, 首先转成灰度图像, 消除算法对彩色图像的依赖, 减少计算量, 降低程序运行时间.

无人机飞行过程中, 由于飞行会使得采集的视频模糊, 具有部分噪声, 本算法采用双边滤波的方法去除噪声. 均值滤波^[4]直接对邻域内像素值取平均, 对图像进行模糊, 简单快速, 飞行过程本身会产生一部分的模糊, 若使用均值滤波会使得视频帧更加模糊; 高斯滤波^[4]是一种适用性较强的去噪算法, 相比均值滤波, 高斯滤波考虑到邻域内像素点和距离的关系, 即距离中心像素越近证明和该像素点越相似, 权值对应高一点, 这样效果较好. 双边滤波^[5], 结果实际上由他邻域内的像素点决定. 类似于均值滤波后的结果是根据相邻像素值的均值取得的, 只不过各像素值权重为1. 双边滤波的权重则是根据相邻像素的距离和相似度得到的. 所以, 双边滤波相比高斯滤波和均值滤波, 既能对像素平滑也能保持一定的细节信息. 室内采集视频时, 主要面对地板, 水泥地面等弱纹理区域, 会包含部分缝隙, 这种区域在我们后期参数估计时是非常好的特征, 双边滤波可以在去噪的基础上保存下边缘信息, 但是由于双边滤波计算时间较长, 不能保证实时性, 因此选用快速双边滤波^[5], 对视频帧进行去噪, 经过多次实验,x方向参数取10,y方向参数取30, 得到的滤波效果是最好的, 即保留了边缘, 又去除了噪声.

3 参数估计

无人机飞行过程中, 参数估计主要包含两方面: 平移参数估计和旋转参数估计. 针对平移参数估计, 首先判断视频类型, 是均为弱纹理的室内地面, 还是存在非弱纹理区域, 根据区域不同的特征选择不同的算法进行估计.

小型无人机室内飞行时, 摄像头采集地面的视频, 面对室内的环境时, 存在部分地板和大理石地面, 这些区域属于重复纹理区域, 本文将该区域划分为弱纹理区域. 室内地面大都是纹理均匀的地板, 但是偶尔会存在其他区域, 如图3所示. 针对不同种类的图像, 运动估计的方法也有所差别, 所以在估计之前, 应先判断视频种类, 根据视频选择对应的方法, 节约时间, 加快参数估计速度. 因为含有非弱纹理区域的图像与均为弱纹理区域的图像, 所以需要对所属区域进行判断, 本方法使用灰度共生矩阵(GLCM)^[6]中的特征值进行判断. 由于灰度共生矩阵计算量较大, 所以应在不影响图像纹理特征的前提下, 对灰度进行压缩, 压缩成16个灰度级, 以1为步距、8×8的滑动窗口, 采集0°, 45°, 90°, 135°四个方向的灰度共生矩阵, 计算每个矩阵的熵和反差分矩阵特征值, 对所属区域进行判断, 再选择后续算法进行参数估计.

3.1 平移参数估计

通过以上对视频种类的判断, 存在两种类型视频, 针对这两种类型的视频, 我们采用两种方法进行运动矢量估计. 当视频均为弱纹理区域时, 采用ORB光流特征算法; 当图像存在其他特征时, 采用改进Snake模型, 将对应区域圈出来, 计算重心, 进而估计出图像的运动矢量.

3.1.1 非弱纹理区域参数估计

非弱纹理区域通俗来说就是弱纹理区域包含一些其他东西, 如室内地面上会存在书本、凳子等物体, 将弱纹理的地板进行了遮掩, 此时选用Snake算法, 将该区域分割出来, 进而计算重心偏移得出此时的无人机飞行参数.

(1) Snake模型

实际上, 图像分割有大量算法, 例如阈值法、区域分割法、边缘分割法等^[7]. 最简单的阈值法将图像根据阈值分割为二值图像, 忽略了图像中的纹理信息, 只关注图像本身像素的值, 因此在很多时候效果差强人意; 区域分割方法需要给定一个种子点, 进行区域分割, 但是会造成区域不连续的效果; 边缘分割法面对的图像需要是特征丰富的, 边缘明显的, 这样分割出的效果才是理想的. 上述几种图像分割方法对噪声敏感, 无人机采集视频会存在部分噪声, 因此几种分割方法都会产生不理想的分割效果.

Snake模型属于可变模型, 使用闭合的参数和曲线对目标区域进行描绘^[7]. 简单来讲, Snake模型就是一条可变形的参数曲线及相应的能量函数, 以最小化能量目标函数为目标, 控制参数曲线变形, 具有最小能量的闭合曲线就是目标轮廓^[7]. Snake模型由一组节点组成, 如式(1)所示, x(s),y(s)代表图像中作为控制点的坐标, s是描述边界的自变量:

$ v(s) = [x(s) - y(s)] $

(1)

在Snakes的控制点上定义能量函数如式(2)所示(反映能量与轮廓之间的关系):

$ {E_{{\text{total}}}} = \int {\left(\alpha {{\left| {\frac{{\overline \partial }}{{{\partial _s}}}v} \right|}^2} + \beta {{\left| {\frac{{\overline \partial }}{{{\partial _s}}}v} \right|}^2} + {E_{\rm {ext}}}(\overline {v(s)} )\right)ds} $

(2)

其中, α,β参数理和口操纵着模型的物理行为和局部连续性. 根据Snake模型的函数, 可以知道, 想要得到目标的轮廓, 需要通过能量函数动态化, 进行多次迭代, 逐渐逼近, 最终得到理想的轮廓. Snake模型的优点就是不论画面质量如何, 都可以得到封闭、光滑的目标轮廓, 相对于其他图像分割算法, 更适合处理差异性较大、结构复杂的图像.

但是相对较费时, 本文使用改进的Snake算法——分段迭代B-Snake算法^[6]. 改进算法的原理是通过一次迭代找到目标的曲线单元, 并进行分离, 在下一次迭代中, 只对逼近曲线单元进行迭代; 当所有曲线迭代结果接近目标时, 对各部分曲线组合生成目标区域的轮廓. 本次采用文献[8]的改进分段迭代B-Snake算法, 进行5次分段迭代, 可以得到理想的轮廓, 如图4所示, 由于现实生活中会存在许多复杂轮廓的物体, 因此本次选用440×350的轮廓较不规则的海星图像进行实验, 提取出图像中海星的轮廓特征. 改进算法快速收敛到了目标轮廓, 节约时间, 减少了迭代次数, 提升速度.

(2)重心计算

根据上述Snake模型算法, 计算两帧图像的区域重心, 通过重心的比较可以得出此时无人机的平移分量. 如图5所示, 为前后两帧图像中重心所在区域截图, 数字代表该位置坐标, 阴影部分即为经过计算的重心点, 假设上一帧重心点坐标为(34,36), 下一帧为(34,37), 采用式(3)对坐标相减, 即可得到此时无人机位移量, 即1 cm.

$ H = H_1({{x_1,y_1}}) - H_2({{x_2,y_2}}) $

(3)

其中, $H_1({{x_1,y_1}})$ 和 $H_2({{x_2,y_2}})$ 为前后两帧视频图像Snake算法圈出的轮廓区域的重心, H为两帧重心的偏移量.

3.1.2 弱纹理区域参数估计

室内地板多为重复纹理和弱纹理区域, 特征不明显, 使用起初的特征匹配产生过多的误匹配, 导致匹配后的参数估计产生大量错误, 因此不可使用多数文献中使用的特征匹配估计参数的方法^[9]. 本文在ORB光流算法的基础上, 结合无人机特定时刻的飞行速度, 使用一种适用于无人机面对室内弱纹理地板区域时的细分匹配规则算法—改进ORB-LK算法, 对平移参数进行估计.

(1) ORB特征提取

ORB特征描述具有较高的实时性, 是由描述子和关键点两部分组成. 改进的FAST角点作为关键点, BRIEF^[9]作为描述子. 首先将采集到的视频帧, 建立图像金字塔, 采集每层金字塔中的FAST角点特征, 由于FAST特征不具备旋转不变性, 本文对FAST特征添加方向因子; 然后, 利用BRIEF算子对特征点进行描述.

(2) LK光流算法

LK (Lucas-Kanade)光流算法^[10]是利用两帧之差来追踪特征点瞬时速度的估计方法, 通过迭代计算可以得到一种稀疏光流场. 但是LK光流法必须满足3条前提假设^[10]: (1)相邻帧之间的亮度、颜色不发生变化; (2)相邻帧之间的像素点运动是足够小且连续的; (3)相邻视频帧之间, 像素的运动一致, 拥有相同的光流. 因此, 可得亮度守恒约束方程^[10]如式(4):

$ I_xv_x + I_yv_y + I_t = 0 $

(4)

其中, $ v_x $ 和 $ v_y $ 分别表示图像在x和y方向上对应的光流值; $ I_x $ 和 $ I_y $ 表示像素点在x和y方向上对应的偏导数; $I_{{t}}$ 表示图像像素点亮度关于时间的导数.

(3)改进ORB-LK算法

无人机飞行过程中, 光流传感器开始启用^[11], 首先进行特征提取, 在接下来计算光流值. 当面向室内弱纹理区域以及重复纹理区域时, 存在部分区域含有相似的纹理或者特征较少, 在参数估计过程中在匹配过程中若使用暴力匹配法这种匹配算法, 会造成大量的错误匹配, 使估计的参数不准确, 因此本次选用FLANN匹配算法.

在原始ORB光流算法^[12]基础上, 结合无人机飞行的速度 $ {v_u} $ , 可以计算出光流的阈值如式(5)所示:

$ \vartriangle = v = v_u \cdot f/h $

(5)

其中, h为无人机的飞行高度, 根据具体情况而定, 假设为一个定值, f表示摄像头焦距, 根据经验, 光流阈值为5. 当光流阈值大于 $ \Delta $ 时, 使用ORB提取特征, 使用双向单追踪策略对特征点进行筛选, 称为初步匹配; 当光流阈值小于 $ \Delta $ 时, 使用FLANN完成特征匹配, 匹配过程中可能会存在误匹配情况, 使用RANSAC进行剔除, 进行二次匹配. 初步匹配的目的在于找到视频帧中尽可能多的特征点, 为二次匹配做准备, 提高二次匹配的速度; 二次匹配目的是优化特征匹配, 可以得到准确的参数估计, 从而对无人机进行精确的定位. (1)初步匹配: 首先对图像使用ORB算法进行特征提取, 得到特征点集M, 再使用金字塔LK算法找到下一帧的图像特征点集N(前向跟踪), 然后利用金字塔将特征点集N找到前一帧对应的点集O(后向跟踪), 由于噪声原因, 使得M和O的点集不能完全重合, 存在一定的误差距离, 遍历M和O点集, 计算距离的误差p, 选取距离小于σ的特征点集, 再结合点集N, 得出筛选后的点集 $M_iN_{i + 1}$ . (2)二次匹配: 当初步匹配过程中光流阈值小于 $ \Delta $ 时, 进行二次匹配. 根据初步匹配点集 $M_iN_{i + 1}$ 可以得到第i帧的特征点集 $A_i$ , ORB特征继续提取第i+1帧对应的特征点集 $A_i$ . 然后使用FLANN进行匹配, 匹配会存在误匹配, 此时使用RANSAC算法剔除误匹配点. 匹配效果如图6所示.

根据上述的结果, 再采用前后双向追踪策略, 如图7所示. 使用第i帧的特征点集 $A_i$ 进行前向追踪, 得点集 $F_{i + 1}$ , 再接下来后向追踪, 获得点集 $B_i$ , 同样由于噪声的存在, 导致 $A_i$ 和 $B_i$ 不能完全重合, 存在一定的误差, 距离误差记为 $d_i$ ; 同理, 可以得到 $A_{i + 1}$ 的前向跟踪点集 $ {F_i} $ 与后向跟踪点集 $ {B_{i + 1}} $ , 同样可以得到距离误差 $ {d_{i + 1}} $ . 光流法存在3个假设, 所以可以得第i帧点集和 $ {F_i} $ 点集距离误差为 $ {d_j} $ ; 第i+1帧, $A_{i + 1}$ 点集和 $ {F_{i + 1}} $ 点集距离误差为 $ {d_{_{j + 1}}} $ . 经过实验, 当误差阈值σ取10时, 效果最好, 进而得出匹配过滤函数W, 如式(6)所示:

$ W = \left\{ \begin{array}{l} {\rm {false}}\,\,\,{\rm {if}}({d_i}{\& }{{{d}}_{{{i}} + {\text{1}}}}{\& }{{{d}}_{{j}}}{\& }{{{d}}_{{{j}} + {\text{1}}}}) > \sigma \\ {\rm {true}}\,\,\,\,{\rm {if}}({d_i}{\& }{{{d}}_{{{i}} + {\text{1}}}}{\& }{{{d}}_{{j}}}{\& }{{{d}}_{{{j}} + {\text{1}}}}) \le \sigma \\ \end{array} \right. $

(6)

3.2 平移参数估计算法描述

算法1. 平移参数估计

输入: 拍摄视频S, 前后相邻帧像素 $\scriptstyle{H}_{1}\left({x}_{1},{y}_{1}\right) $ , $ \scriptstyle{H}_{2}\left({x}_{2},{y}_{2}\right) $

输出: 平移参数H

1 for each frame in S

2 { 计算每帧GLCM的熵和反差特征值m,n

3 if(m>3,n<15)

4 { 由式(1)和式(2)通过B-Snake, 得到的非纹理区域轮廓

5 由式(3)计算前后两帧轮廓重心偏移量H}

6 else

7 { 提取ORB特征值

8 由式(5)计算光流值Δ}

9 if(Δ>5)

10 {通过光流与惯导系统得到参数}

11 else

12 {FLANN特征匹配

13 RANSAC剔除错误匹配点

14 由式(5)计算光流值Δ

15 通过光流与惯导系统得到参数}}}}

3.3 旋转参数估计

无人机在飞行过程中, 不仅会平移产生平移参数, 同时也会存在一定的旋转运动产生旋转参数, 所以参数估计过程中, 对旋转参数也要进行估计^[13]. 单目摄像机拍摄视频对前后两帧特征点进行匹配, 到处旋转矩阵, 即可得到旋转参数. 在飞行过程中, 主要旋转角度分为绕Z轴的航向角(ψ), 绕X轴的俯仰角(Ф)和绕Y轴的翻滚角(θ)^[14], 如图8所示.

根据特征匹配, 可以得到前后两帧图像间的单应性矩阵H, 对矩阵H进行分解即可得到旋转参数, 如式(7)所示:

$ {M_{n + 1}} = H{M_n} $

(7)

其中, ${M_{{{n}} + 1}}$ 和 ${M_{{n}}}$ 分别为第n+1帧与第n帧图像的像素值, 矩阵H中包含旋转矩阵 $ R({r_x},{r_y},{r_{\textit{z}}}) $ . 根据实验数据处理结果, 对旋转矩阵进行解算, 由于数据误差, 旋转矩阵不能精准满足 $R = R{R^{\rm T}} = {R^{\rm T}}R = I$ . 因此对旋转矩阵进行SVD分解^[14], 令 $R = UD{V^{\rm T}}$ , D =I, 使旋转矩阵R强制转为满足条件的旋转矩阵 ${R_1} = UI{V^{\rm T}}$ . 特征点在机体坐标系和世界坐标系下对应关系如式(8):

$ \left\{ \begin{array}{l} M_n^{'} = C_n^b({M_n} + {t_0}) \\ C_n^b = C_c^b{R_1}C_{n0}^c \\ \end{array} \right. $

(8)

其中, $M_{{n}}^{'}$ 为特征点在机体坐标系下的三维坐标, ${M_{{n}}}$ 和 $ {t_0} $ 为世界坐标系下的三维坐标和平移分量, $ C_n^b $ 是世界坐标系到体轴系的转换矩阵, $ C_c^b $ 是相机坐标系到体轴系的转换矩阵, $ C_{n0}^c $ 是源图像中相机坐标系到体轴系的转换矩阵, 是一个常值矩阵. 所以拍摄目标源图像时, 图像平面与目标平面平行, 通过分解 $ C_n^b $ 即可得到3个轴对应的姿态角, 如式(9)所示:

$ C_n^b = \left(\begin{array}{*{20}{c}} {{c_{11}}}&{{c_{12}}}&{{c_{13}}} \\ {{c_{21}}}&{{c_{22}}}&{{c_{23}}} \\ {{c_{31}}}&{{c_{32}}}&{{c_{33}}} \end{array}\right) = \left(\begin{array}{*{20}{c}} {\sin \varphi \sin \theta \sin \psi + \cos \theta \cos \psi }&{\sin \psi \cos \varphi }&{\cos \theta \sin \varphi \sin \psi - \sin \theta \cos \psi } \\ {\sin \varphi \sin \theta \cos \psi - \cos \theta \sin \psi }&{\cos \psi \cos \varphi }&{\cos \theta \sin \varphi \cos \psi + \sin \theta \sin \psi } \\ {\sin \theta \cos \varphi }&{ - \sin \varphi }&{\cos \psi \cos \theta } \end{array}\right) $

(9)

3.4 旋转参数估计算法描述

算法2. 旋转参数估计

输入: 拍摄视频S

输出: 旋转参数Ф,θ,ψ

1 for each frame in S

2 {提取ORB特征值

3 由式(7)计算单应性矩阵H

4 对单应性矩阵进行分解, 得到旋转矩阵R

5 对R进行SVD分解

6 由式(8)进行世界坐标系与相机坐标系的转换

7 由式(9)分解姿态角Ф,θ,ψ}

4 实验结果

根据上述算法, 本实验使用分辨率为440×350, 焦距为50 mm的摄像机作为图像采集设备, 使用电脑( CPU为Intel i7 3520M, 主频达到2.9 GHz, 内存为4 GB)作为图像处理计算机, 结合OpenCV 3.3.1和Ubuntu 16.04环境下进行仿真.

针对室内弱纹理条件下, 选择了多个场景进行实验. 为了验证算法的准确性, 在飞行器平稳飞行过程中, 对无人机制造一定的移动, 通过算法进行估计, 在平移和旋转参数已知情况下, 进行参数估计, 将估计量与偏移量进行比较. 针对不同的偏移量进行估计, 得到如表1–表3的结果.

由表1–表3可以证明, 本方法在无人机飞行过程中, 平均处理时间达到了25 ms, 准确率为95%. 文献[3]提出一种多跟踪的无人机姿态参数估计系统, 在参数估计方向属于较新的一种方法, 本文相较于该方法, 运行时间平均缩短了10 ms, 准确率有大幅度提升; 文献[14]采用改进的ORB光流算法进行参数估计, 本文在该基础上进行进一步改进, 增加了基本的去噪方法, 使得准确率有大幅度的提升.

5 结语

通过图像的预处理, 对采集的视频图像进行去噪, 结合图像的灰度共生矩阵对视频种类进行判断, 根据视频纹理重复类型选择使用改进ORB光流算法进行平移参数估计, 或者是使用Snake模型进行估计, 联合重心偏移计算得出运动平移分量. 本文提出的方法结合灰度共生矩阵进行判断, 两种方法相结合, 避免了使用单一方法, 全面匹配, 估计平移分量的问题. 结合单应性矩阵, 分解得到旋转参数. 本方法对多种室内弱纹理或重复纹理区域进行参数估计, 扩大了处理范围, 降低了处理时间.

本文算法在精确估计出无人机移动参数情况下, 可结合飞控系统, 将估计出的参数传入系统, 在室内GPS信号较弱情况下, 对无人机进行稳定控制, 保证室内采集到稳定清晰的图像材料.