CN106017463A - 一种基于定位传感装置的飞行器定位方法 - Google Patents

Abstract

本发明涉及一种基于定位传感装置的飞行器定位方法，通过惯性传感单元获取飞行器姿态信息，利用超声波传感器收集相对高度信息，利用摄像头获取地面的实时图像信息；利用改进的光流法获取飞行器相对位移信息，利用网格信息获取飞行器全局位置与姿态信息；利用滤波算法实现多传感器信息融合获取飞行器准确鲁棒的定位方法。在无GPS信号等外界辅助导航情况下，能实现飞行器的长期准确定位，且鲁棒性好，无累积误差。

Description

一种基于定位传感装置的飞行器定位方法

技术领域

本发明涉及飞行器定位领域，具体涉及一种基于定位传感装置的飞行器定位方法。

背景技术

由于缺乏外界辅助导航，无人飞行器难以在未知环境下估计无人飞行器的位姿与运动，飞行器自主导航的过程中需要解决这个关键问题。而这个问题的解决方法与飞行器机载传感器的类型紧密联系，因此无人飞行器的硬件构架对相关算法的选择有直接影响。目前针对不同的飞行器硬件构建，为上述问题提出了不同解决方案。

对于地面移动机器人，姿态信息的估计值通常利用惯性测量单元IMU(InertialMeasurement Unit)解算获得。而位置信息的获取，则有赖于一类被称为里程计(Odometry)的导航方法。就其本质而言，里程计是一类采用了匹配与积分技术的算法。根据里程计匹配数据的类型，可将里程计划分为视觉里程计、激光里程计以及编码盘里程计等。里程计的解算精度不仅与匹配算法有关，对传感器的测量误差也十分敏感。无人飞行器由于负载有限，无法沿用许多在地面移动机器人上得到广泛应用的高精度IMU与里程计。为此，需要从重量、大小、实时性与精度等方面对机载传感器的选择进行权衡考虑。在近十年的相关研究中，比较有代表性的组合方案包括如下几种：

1)单/双目视觉系统与IMU的组合。视觉系统在机器人领域的应用较为成熟。单目视觉系统的运动恢复结构算法通过计算序列图像对应的基础矩阵，估算摄像机投影矩阵，可恢复场景特征点的投影结构，获得载体的姿态变化；在载体运动已知的情况下，能够恢复特征点在世界坐标系的真实位置，从而实现基于单目视觉的相对局部环境的位置估计。双目立体视觉系统通过匹配两个摄像头视野的图像特征，利用空间特征的投影关系，获得载体的相对位置信息。但该方法鲁棒性差，可能存在错匹配等情况，将直接影响位姿估计的准确性。

2)光流法与IMU的组合。光流法同样依赖于机载摄像头提供的图像灰度信息。但与1)所不同的是，光流法无需对图像进行特征提取与特征匹配，转而考虑图像灰度在时间上的变化率与摄像头和图像场景之间相对变化的关系。在灰度信息之外，光流法还要求提供图像像素对应的深度信息，从而完成对位姿的解算。但该方法存在累积误差，在飞行器长期飞行的情况下，其定位误差将逐渐累积，因此无法保证飞行器的长期稳定准确飞行定位。

3)激光测距仪与IMU的组合。激光测距仪可以感知环境的结构信息，得到环境的二维点云数据。利用迭代最近邻点ICP(Iterative Closest Point)算法，对相邻时间的两帧点云数据进行匹配，可求取载体在该时间段内的平移运动分量及旋转运动分量。对平移运动分量和旋转运动分量进行积分，进而可获得飞行器的位置和姿态估计值。可利用ICP算法对激光测距仪的二维点云数据进行匹配，得到飞行器在运动平面上的位置及飞行器朝向的估计值，并采用扩展卡尔曼滤波器EKF(Extended Kalman Filter)将ICP算法解算得到的飞行器位姿信息与IMU解算获得的飞行器姿态信息进行融合，实现了无人飞行器在局部场景下的位姿估计。但利用该方法需要在飞行器上安装激光测距仪，该装置体积较大，价格昂贵，并不适用于微飞行器平台的安装及推广使用。

4)深度摄像机(RGB-D摄像机)与IMU的组合。RGB-D摄像机是结合了摄像头与测距仪优点的一种新型传感器，其在获得视野图像灰度信息的同时，也可获得相应的深度信息。同双目立体视觉系统相比，RGB-D摄像机采样深度信息的效率更高，获得的深度信息更为精确。在无人飞行器上配备Microsoft Kinect，利用和双目立体视觉系统类似的原理，从场景图像中提取特征点及其深度信息并进行特征匹配，得到空间特征的投影关系，获得载体的相对运动信息。使用该传感器同样由于体积较大，并不适用于安装在微小飞行器平台中。

发明内容

本发明的目的在于针对现有技术的不足，提供一种基于定位传感装置的飞行器定位方法，本方法在无GPS信号等外界辅助导航情况下，能实现飞行器的长期准确定位，且鲁棒性好，无累积误差。

本发明所提供的技术方案为：一种基于定位传感装置的飞行器定位方法，包括如下步骤：

1)通过惯性传感单元中的加速度计和陀螺仪分别测得飞行器实时的加速度和角速度信息，通过磁力计确定方位，获得飞行器的实时姿态信息；所述的实时姿态信息包括偏航角、横滚角和俯仰角；

利用超声波传感器收集飞行器相对于地面的高度信息，并利用惯性传感单元获取的实时姿态信息进行姿态补偿，实时校正垂直高度信息；

利用摄像头，获取地面的实时图像信息；

2)对获取的实时图像信息作灰度化处理，获取灰度图，对每帧灰度图选择固定间隔固定数量的像素块，通过匹配相邻两帧灰度图，计算光流向量，并通过差分，获取光流速度；

3)同步惯性传感单元、超声波传感器和摄像头的时间戳，利用惯性传感单元获取的实时姿态信息，对光流速度进行姿态补偿，得到光流速度的水平面分量；利用超声波传感器获取的垂直高度信息对光流速度的水平面分量进行尺度还原，得到飞行器的水平速度；

4)对飞行器的水平速度积分得到相对位移信息，累加相对位移信息，即得水平的位置信息，结合超声波传感器得到的垂直高度信息，即可得到飞行器的位置信息。

所述的步骤2)中光流速度的获取方法为：

对每帧灰度图选择固定间隔固定数量的像素块，将前一帧灰度图中的像素块位置作为后一帧灰度图的像素块的初始位置，并建立一个搜索区域，阈值函数的公式如下：

$\underset{i=0}{\overset{N-1}{\Σ}}\underset{j=0}{\overset{N-1}{\Σ}}|G_1(i,j)-G_0(i+u,j+v)|---\left(1\right)$

式(1)中，u表示x轴方向的搜索像素个数，v表示y轴方向的搜索像素个数；G₀和G₁分别表示前一帧和后一帧像素块中某一个像素点的灰度值；i，j分别表示图像平面像素块中某一个像素点的x,y轴方向的坐标，N表示像素块的x轴或者y轴方向上的像素个数，G₁(i,j)表示前一帧像素块中每一像素点的灰度值，G₀(i+u,j+v)表示将前一帧图像中的像素块投射到后一帧图像中相同位置后并沿x轴、y轴方向分别移动u、v像素的像素块中每一像素点的灰度值；

通过令前一帧的像素块在后一阵的灰度图的搜索区域内搜索灰度差和最小，求得后一帧的像素块位置，两帧像素块的距离差，即为光流向量，通过差分，获取光流速度。

进一步，所述的步骤3)中获得飞行器的水平速度的方法为：

建立图像平面中的点与摄像机坐标系中的点的转换关系：

$u=f\frac{p_c}{z}---\left(2\right)$

式(2)中，u为图像平面中的点的坐标，p_c为摄像机坐标系的点的坐标，f为焦距，z为垂直高度；

对式(2)两边求导，即得光流速度：

$\frac{flow}{\mathrm{\Δ}time}=v=f\frac{{\mathrm{zv}}_c-v_z{p}_c}{z^2}---\left(3\right)$

式(3)中，v_c为摄像头坐标系下点的运动速度，v_z为垂直方向上的运动速度，v为图像平面的运动速度，即光流速度；

利用刚体运动学的公式：

v_c＝-v_w-ω×p_c (4)

式(4)中，v_w表示世界坐标系下点的运动速度，ω表示角速度；

将公式(4)带入公式(3)，展开可得到x轴和y轴方向的光流速度与世界坐标系下的飞行器速度的关系：

$\begin{array}{c}{v}_u=f\frac{{\mathrm{zv}}_{cx}-v_z{x}_c}{z^2}=f\frac{z(-v_{wx}-w_{y}z+w_z{y}_c)-(-v_{wz}-w_x{y}_c+w_y{x}_c)x_c}{z^2}\\ =\frac{v_{wx}x-v_{wx}f}{z}-w_{y}f+w_{z}y+\frac{w_{x}uv-w_y{u}^2}{f}\end{array}---\left(5\right)$

$v_v=\frac{v_{wz}y-v_{wy}f}{z}+w_{x}f+w_{z}x+\frac{w_x{v}^2-w_{y}uv}{f}---\left(6\right)$

式(5)中，v_u表示图像平面x轴方向的光流速度，v_cx表示相机坐标系下x轴方向的速度，v_wx表示世界坐标系下x轴方向的速度，v_wz表示世界坐标系下z轴方向的速度，w_x，w_y，w_z分别表示x,y,z轴方向上的角速度，x_c表示相机坐标系下x轴方向的坐标，y_c表示相机坐标系下y轴方向的坐标；u，v分别表示图像平面x,y轴方向的坐标；

式(6)中，v_v表示图像平面y轴方向的光流速度，v_wy表示世界坐标系下y轴方向的速度；v_wz表示世界坐标系下z轴方向的速度；对垂直高度信息差分可以求得世界坐标系下z轴方向的速度，即

根据式(5)和(6)，以及已知的世界坐标系下z轴方向的速度，可以求得世界坐标系下x轴和y轴方向的飞行器的速度，即为飞行器的水平速度。

进一步，当飞行器定高飞行时，式(5)和(6)简化为：

$v_u=-\frac{v_{wx}f}{z}-w_{y}f+w_{z}v---\left(7\right)$

$v_v=-\frac{v_{wy}f}{z}+w_{x}f-w_{z}u---\left(8\right).$

作为改进，对所述的飞行器定位方法进行优化，所述的摄像头获取地面的实时图像信息为规则的网格信息，包括如下步骤：

i)利用摄像头，获取地面的实时图像信息，进一步处理实时图像信息中的规则的网格信息，使用霍夫变换算法，调整筛选直线的阈值，对其中的直线信息进行准确提取，从而获取实时图像信息中的直线交点坐标；

ii)同步惯性传感单元、超声波传感器和摄像头的时间戳，对飞行器的水平速度积分，计算飞行器的初始相对位移；给定飞行器的初始位置估计；利用惯性传感单元给定实时姿态信息，给定飞行器的初始姿态信息估计；使用摄像头的投影模型和实时图像信息中获取的直线交点坐标，获取地面网格交点的坐标；

iii)根据地面网格交点相邻距离恒定的条件，对求得的地面网格交点的坐标进行求精操作，获取准确的地面网格交点的坐标，再利用实时图像信息对应的坐标，建立匹配点对，使用直接线性变换算法或者迭代算法，求得飞行器准确的位置信息及姿态信息。

进一步，所述的步骤ii)中获取地面网格交点的实际坐标的方法为：

根据实时图像信息中的粗略位置估计，公式如下：

$p_k^w=p_{k-1}^w+v_w{d}_t---\left(9\right)$

式(9)中，v_w为水平光流速度，为前一时刻水平位置，为后一时刻水平位置，d_t为间隔时间，另外令超声波传感器在k-1时刻和k时刻给入的垂直高度信息分别为：和给定飞行器的初始位置估计，其中均为三维列向量，表示x,y,z轴方向的位置；

由惯性传感单元给定飞行器实时的姿态信息，包括偏航角、横滚角和俯仰角，如下：

roll＝γ,pitch＝θ,yaw＝ψ (10)

式(10)中，roll、pitch、yaw分别为该时刻的横滚角、俯仰角和偏航角；

对于步骤i)中获取的实时图像信息中的直线交点坐标(u,v)，都满足直线方程：

au+bv+c＝0 (11)；

摄像头的投影模型由以下公式表达：

式(12)中，λ表示相机坐标系下特征点的z轴坐标，(u,v,1)^T表示实时图像信息中特征点的齐次坐标，K表示摄像头的内参矩阵，f_x,f_y,c_x,c_y分别表示x,y轴方向的像素焦距以及主点的像素坐标；表示世界坐标系相对于相机坐标系的转换矩阵，即外参矩阵，前三行三列矩阵表示旋转矩阵，最后一列矩阵表示平移向量；(x_w,y_w,z_w,1)^T表示世界坐标系下特征点的齐次坐标；

将直线方程公式(11)代入摄像头的投影模型公式(12)，得到公式如下：

给定飞行器的初始姿态信息估计和初始位置估计，如下：

$\left(\begin{array}{ccc}{r}_{11}& r_{12}& r_{13}\\ r_{21}& r_{\mathrm{22}}& r_{\mathrm{23}}\\ r_{31}& r_{\mathrm{32}}& r_{\mathrm{33}}\end{array}\right)=\left(\begin{array}{ccc}\cos \mathrm{\γ}\cos \mathrm{\ψ}-\sin \mathrm{\γ}\sin \mathrm{\θ}\sin \mathrm{\ψ}& -\cos \mathrm{\θ}\sin \mathrm{\ψ}& \sin \mathrm{\γ}\cos \mathrm{\ψ}+\cos \mathrm{\γ}\sin \mathrm{\θ}\sin \mathrm{\ψ}\\ \cos \mathrm{\γ}\sin \mathrm{\ψ}+\sin \mathrm{\γ}\sin \mathrm{\θ}\cos \mathrm{\ψ}& \cos \mathrm{\θ}\cos \mathrm{\ψ}& \sin \mathrm{\γ}\sin \mathrm{\ψ}-\cos \mathrm{\γ}\sin \mathrm{\θ}\cos \mathrm{\ψ}\\ -\sin \mathrm{\γ}\cos \mathrm{\θ}& \sin \mathrm{\θ}& \cos \mathrm{\θ}\cos \mathrm{\γ}\end{array}\right)---\left(14\right)$

$\left(\begin{array}{c}{t}_x\\ t_y\\ t_z\end{array}\right)=\left(\begin{array}{c}p_k^w{}_x\\ p_k^w{}_y\\ p_k^w{}_z\end{array}\right)---\left(15\right)$

将给定的参数代入公式(13)，求解公式(13)的直线方程，设求出的世界坐标系的坐标点p_w为：

p_w＝(x_w',y_w',z_w') (16)

对世界坐标系中的直线交点(x_w,y_w,z_w)，同样建立直线方程：

a'x_w+b'y_w+c'z_w+d'＝0 (17)

利用规则的网格信息的特征，直线只分为两类，两类直线互相垂直，且认为地面为平面，则z_w＝0，将直线方程分为两类：

$\left\{\begin{array}{c}{y}_w=-\frac{d^{\′}}{b^{\′}},if{a}^{\′}=0\\ x_w=-\frac{d^{\′}}{a^{\′}},if{b}^{\′}=0\end{array}\right.---\left(18\right)$

由公式(18)对世界坐标系的坐标点进行求精操作，即水平和垂直的两类直线分布间隔均匀，地面的直线实际交点坐标为常值，以此可以对p_w＝(x_w',y_w',z_w')求出的坐标中的x_w'，y_w'取与之相差绝对值最小的常值x_w″,y_w″，z_w'认为是0；求精后的世界坐标系下的直线交点的坐标值为p_w＝(x_w″,y_w″,0)；

利用公式(13)代入求精后p_w的坐标，如下：

简化后得到：

$\begin{array}{c}\[{\mathrm{af}}_x{r}_{11}+{\mathrm{bf}}_y{r}_{21}+({\mathrm{ac}}_x+{\mathrm{bc}}_y+c)r_{31}\]{x_w}^{\′\′}+\[{\mathrm{af}}_x{r}_{12}+{\mathrm{bf}}_y{r}_{22}+({\mathrm{ac}}_x+{\mathrm{bc}}_y+c)r_{32}\]{y_w}^{\′\′}\\ +\[{\mathrm{af}}_x{t}_x+{\mathrm{bf}}_y{t}_y+({\mathrm{ac}}_x+{\mathrm{bc}}_y+c)t_z\]=0\end{array}---\left(20\right)$

由公式(18)知直线方程有两种，故将方程(19)分为两类：

$\begin{array}{c}\[{\mathrm{af}}_x{r}_{11}+{\mathrm{bf}}_y{r}_{21}+({\mathrm{ac}}_x+{\mathrm{bc}}_y+c)r_{31}\]=0\\ \[{\mathrm{af}}_x{r}_{12}+{\mathrm{bf}}_y{r}_{22}+({\mathrm{ac}}_x+{\mathrm{bc}}_y+c)r_{32}\]{y_w}^{\′\′}+\[{\mathrm{af}}_x{t}_x+{\mathrm{bf}}_y{t}_y+({\mathrm{ac}}_x+{\mathrm{bc}}_y+c)t_z\]=0\\ \[{\mathrm{af}}_x{r}_{12}+{\mathrm{bf}}_y{r}_{22}+({\mathrm{ac}}_x+{\mathrm{bc}}_y+c)r_{32}\]=0\\ \[{\mathrm{af}}_x{r}_{11}+{\mathrm{bf}}_y{r}_{21}+({\mathrm{ac}}_x+{\mathrm{bc}}_y+c)r_{31}\]{x_w}^{\′\′}+\[{\mathrm{af}}_x{t}_x+{\mathrm{bf}}_y{t}_y+({\mathrm{ac}}_x+{\mathrm{bc}}_y+c)t_z\]=0\end{array}---\left(21\right)$

提取未知数r₁₁,r₂₁,r₃₁,r₁₂,r₂₂,r₃₂,t_x,t_y,t_z，重写方程(13)，有：

$\left(\begin{array}{ccccccccc}{\mathrm{af}}_x& {\mathrm{bf}}_y& {\mathrm{ac}}_x+{\mathrm{bc}}_y+c& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& {\mathrm{af}}_x{y_w}^{\′\′}& {\mathrm{bf}}_y{y_w}^{\′\′}& ({\mathrm{ac}}_x+{\mathrm{bc}}_y+c){y_w}^{\′\′}& {\mathrm{af}}_x& {\mathrm{bf}}_y& {\mathrm{ac}}_x+{\mathrm{bc}}_y+c\\ 0& 0& 0& {\mathrm{af}}_x& {\mathrm{bf}}_y& {\mathrm{ac}}_x+{\mathrm{bc}}_y+c& 0& 0& 0\\ {\mathrm{af}}_x{x_w}^{\′\′}& {\mathrm{bf}}_y{x_w}^{\′\′}& ({\mathrm{ac}}_x+{\mathrm{bc}}_y+c){x_w}^{\′\′}& 0& 0& 0& {\mathrm{af}}_x& {\mathrm{bf}}_y& {\mathrm{ac}}_x+{\mathrm{bc}}_y+c\end{array}\right)\left(\begin{array}{c}{r}_{11}\\ r_{21}\\ r_{31}\\ r_{12}\\ r_{\mathrm{22}}\\ r_{\mathrm{32}}\\ t_x\\ t_y\\ t_z\end{array}\right)=0---\left(22\right);$

方程(22)是在取两类直线中的各一条建立的方程组，由于存在9个未知数，至少需要五条直线才能解该线性方程组；但直线数超过5条后，该方程组为超定方程组，对其进行奇异值分解即可解该方程组；

至此，解出旋转矩阵R和平移向量t，如下：

$R=\left(\begin{array}{ccc}{r}_{11},& r_{12}& r_{13}\\ r_{21}& r_{22}& r_{23}\\ r_{31}& r_{32}& r_{33}\end{array}\right),t=\left(\begin{array}{c}{t}_x\\ t_y\\ t_z\end{array}\right)---\left(23\right)$

在已求出r₁₁,r₂₁,r₃₁,r₁₂,r₂₂,r₃₂的情况下，可利用旋转矩阵是单位正交矩阵这一特性，求出r₁₃,r₂₃,r₃₃；

再将旋转矩阵转换为姿态角，得到姿态信息：

roll＝a tan2(r₃₂,r₃₃)

$pitch=atan2(-r_{31},\sqrt{{r_{32}}^2+{r_{33}}^2})---\left(24\right).$

yaw＝a tan2(r₂₁,r₁₁)

进一步改进，对所述的飞行器定位方法进一步优化，采用光流速度信息，网格定位的位置信息、惯性传感单元的加速度信息，建立状态方程和观测方程，对视觉算法和惯性传感单元获取的结果使用扩展卡尔曼滤波算法，进行融合处理，获取飞行器准确的位置及姿态信息。

所述的获取飞行器准确的位置及姿态信息具体的算法推导如下：

非线性连续状态方程如下：

x_k+1＝f(x_k,w_k)

对其进行离散化操作，状态量包括有网格定位得到的位置信息p，光流法提供的光流速度信息v和惯性传感单元提供的加速度a；Δt表示连续两次采样的间隔时间，A表示状态转移矩阵；

观测方程如下，分为视觉算法部分的观测方程，观测量位置p；惯性传感单元的观测方程，观测量取加速度a：

$\begin{array}{c}{z}_{v,k}=H_{v,k}{x}_k=\left[\begin{array}{ccc}I& 0& 0\end{array}\right]x_k\\ z_{i,k}=H_{i,k}{x}_k=\left[\begin{array}{ccc}0& 0& I\end{array}\right]x_k\end{array}---\left(26\right)$

其中公式(26)中，z_v,k，z_i,k分别是视觉算法和惯性传感单元的观测量，分别是世界坐标系下的水平方向的位置信息和加速度信息；

根据卡尔曼滤波算法的公式：

$\begin{array}{c}{P}_k^{-}={\mathrm{AP}}_{k-1}{A}^T+Q\\ K_{v,k}=P_k^{-}{H}_{v,k}^T{(H_{v,k}{P}_k^{-}{H}_{i,k}^T+R_v)}^{-1}\\ x_k=x_k^{-}+K_{v,k}(x_{vision}-H_{v,k}{x}_k^{-})\\ P_k=(I-K_{v,k}{H}_{v,k})P_k^{-}\end{array}---\left(27\right)$

$\begin{array}{c}{P}_k^{-}={\mathrm{AP}}_{k-1}{A}^T+Q\\ K_{i,k}=P_k^{-}{H}_{i,k}^T{(H_{i,k}{P}_k^{-}{H}_{i,k}^T+R_i)}^{-1}\\ x_k=x_k^{-}+K_{i,k}(x_{IMU}-H_{i,k}{x}_k^{-})\\ P_k=(I-K_{i,k}{H}_{i,k})P_k^{-}\end{array}---\left(28\right)$

可分别计算视觉算法部分和惯性传感单元部分的卡尔曼增益K_v,k，K_i,k，其中P_k-1表示对应于状态x_k-1的协方差，表示对应于状态的协方差，Q表示系统噪声的协方差矩阵，R_v，R_i分别表示视觉算法和惯性传感单元噪声的协方差矩阵；再结合预测值和测量值，分别求取视觉部分和惯性传感单元部分状态量的最优估值x_k，其中x_vision，x_IMU分别表示视觉算法的测量值和惯性传感单元的测量值；最后再次更新相对于状态量x_k的协方差矩阵。

同现有技术相比，本发明的有益效果体现在：

(1)本发明利用改进的光流法获取飞行器相对位移信息；利用网格信息获取飞行器全局位置与姿态信息；利用滤波算法实现多传感器信息融合获取飞行器准确鲁棒的定位方法。

(2)采用光流与网格信息相结合的组合定位方法，各自巧妙地将两者的优势结合起来，相比单纯使用光流信息，避免了累计误差，定位效果可保证长期较高的准确度。

(3)采用视觉算法与惯性传感单元信息融合的扩展卡尔曼滤波算法，不仅确保了提供飞行器飞行轨迹的平滑性，更进一步保证了定位算法的鲁棒性。

附图说明

图1为本发明飞行器的定位传感装置的结构示意图；

图2为本发明实现飞行器定位的场景图；

图3为本发明实现飞行器定位的传感数据处理流程图；

图4为本发明基于定位传感装置的飞行器定位方法的流程图。

具体实施方式

如图1所示为飞行器的定位传感装置，由微处理器(NUC)1，惯性传感单元(IMU)2，超声波传感器3，摄像头4组成。

微处理器采用Intel公司的NUC(NUC5i7RYH)系列，整重仅0.607Kg，该系列处理器具备体积纤小、多接口、处理速度快、功能强大、低功耗、散热快等特点，本发明由于摄像头的高帧率采集到的图像信息庞大，对数据处理的实时能力有很高的要求；且所有的数据处理要求均在linux(ubuntu)系统下进行；要求整套装置能安装于微型飞行器上，故基于以上要求，由于NUC5i7RYH的大容量数据处理能力、轻薄外观、跨平台的特性(全面支持Windows、Linux)，选择它非常合适。

摄像头与超声波模块，选用PX4FLOW，它是一款开源智能光学流动传感器，拥有752×480像素分辨率以及250Hz刷新率，计算光学流的过程中采用了4倍分级和裁剪算法，具备非常高的感光度，能够在室内或者室外暗光环境下工作，且硬件和软件结构均开源，适合对设计的算法进行重新编程，二次开发；

惯性传感单元采用LPMS系列，LP-RESEARCH公司该系列产品使用了先进的数据融合算法，为用户提供高精度高稳定性的姿态动作信息以及3轴加速度/角速度/地磁量/方位角等数据，为算法实现提供准确的姿态信息。

如图2所示：飞行器在大理石地面、地板等有规则和纹理地面条件下飞行示意图。下标w的坐标系意为世界坐标系，即在整个场景下建立的全局坐标系。飞行器上下标为c的坐标系意为机体坐标系，由于设计了多传感固连装置，惯性传感单元、摄像头、超声波传感器、机体均为刚性连接，因此各自坐标系只存在一个固定的坐标转换关系，该关系在设计固连装置时已给定，如此解决了各传感器之间的坐标转换问题。坐标轴为u,v的为图像坐标系，任意地面一点投影到图像平面将在该平面建立一个像素单位的坐标。相机为下视视野，能够向下观测到部分地面的网格和纹理信息。

如图3所示：定位传感装置中的传感器套件进行初始化操作。在微处理器下运行各节点，主要包括图像节点，用于接收地面图像信息；超声波节点，用于检测相对高度信息；惯性传感单元节点，用于接收整个模块的角速度，加速度信息；时间同步节点，用于对多传感器数据进行同步操作；处理运算节点，对同步后的数据进行定位算法运算处理。运算主节点将订阅其他节点的信息，若无数据通信，则保持通讯状态；若接收到传感器信息，将首先对各传感器的时间戳进行识别，以时间戳作为同步原则，对不同帧率的传感器数据首先进行同步处理。接下来，运算主节点将根据本实施例设计的算法，对处理后的数据进行运算，实时输出飞行器的位置与姿态，作为飞行器控制器的输入量。

如图4所示，基于定位传感装置的飞行器定位方法：包括如下步骤：

1)通过惯性传感单元中的加速度计和陀螺仪分别测得飞行器实时的加速度和角速度信息，通过磁力计确定方位，获得飞行器的实时姿态信息；所述的实时姿态信息包括偏航角、横滚角和俯仰角；

利用超声波传感器收集飞行器相对于地面的高度信息，并利用惯性传感单元获取的实时姿态信息进行姿态补偿，实时校正垂直高度信息；

利用摄像头，获取地面的实时图像信息；

2)对获取的实时图像信息作灰度化处理，获取灰度图，对每帧灰度图选择固定间隔固定数量的像素块，通过匹配相邻两帧灰度图，计算光流向量，并通过差分，获取光流速度；

3)同步惯性传感单元、超声波传感器和摄像头的时间戳，利用惯性传感单元获取的实时姿态信息，对光流速度进行姿态补偿，得到光流速度的水平面分量；利用超声波传感器获取的垂直高度信息对光流速度的水平面分量进行尺度还原，得到飞行器的水平速度；

4)对飞行器的水平速度积分得到相对位移信息，累加相对位移信息，即得水平的位置信息，结合超声波传感器得到的垂直高度信息，即可得到飞行器的位置信息。

所述的步骤2)中光流速度的获取方法为：

对每帧灰度图选择固定间隔固定数量的像素块，将前一帧灰度图中的像素块位置作为后一帧灰度图的像素块的初始位置，并建立一个搜索区域，阈值函数的公式如下：

$\underset{i=0}{\overset{N-1}{\Σ}}\underset{j=0}{\overset{N-1}{\Σ}}|G_1(i,j)-G_0(i+u,j+v)|---\left(1\right)$

式(1)中，u表示x轴方向的搜索像素个数，v表示y轴方向的搜索像素个数；G₀和G₁分别表示前一帧和后一帧像素块中某一个像素点的灰度值；i，j分别表示图像平面像素块中某一个像素点的x,y轴方向的坐标，N表示像素块的x轴或者y轴方向上的像素个数，G₁(i,j)表示前一帧像素块中每一像素点的灰度值，G₀(i+u,j+v)表示将前一帧图像中的像素块投射到后一帧图像中相同位置后并沿x轴、y轴方向分别移动u、v像素的像素块中每一像素点的灰度值；

通过令前一帧的像素块在后一阵的灰度图的搜索区域内搜索灰度差和最小，求得后一帧的像素块位置，两帧像素块的距离差，即为光流向量，通过差分，获取光流速度。

进一步，所述的步骤3)中获得飞行器的水平速度的方法为：

建立图像平面中的点与摄像机坐标系中的点的转换关系：

$u=f\frac{p_c}{z}---\left(2\right)$

式(2)中，u为图像平面中的点的坐标，p_c为摄像机坐标系的点的坐标，f为焦距，z为垂直高度；

对式(2)两边求导，即得光流速度：

$\frac{flow}{\mathrm{\Δ}time}=v=f\frac{{\mathrm{zv}}_c-v_z{p}_c}{z^2}---\left(3\right)$

式(3)中，v_c为摄像头坐标系下点的运动速度，v_z为垂直方向上的运动速度，v为图像平面的运动速度，即光流速度；

利用刚体运动学的公式：

v_c＝-v_w-ω×p_c (4)

式(4)中，v_w表示世界坐标系下点的运动速度，ω表示角速度；

将公式(4)带入公式(3)，展开可得到x轴和y轴方向的光流速度与世界坐标系下的飞行器速度的关系：

$\begin{array}{c}{v}_u=f\frac{{\mathrm{zv}}_{cx}-v_z{x}_c}{z^2}=f\frac{z(-v_{wx}-w_{y}z+w_z{y}_c)-(-v_{wz}-w_x{y}_c+w_y{x}_c)x_c}{z^2}\\ =\frac{v_{wx}x-v_{wx}f}{z}-w_{y}f+w_{z}y+\frac{w_{x}uv-w_y{u}^2}{f}\end{array}---\left(5\right)$

$v_v=\frac{v_{wz}y-v_{wy}f}{z}+w_{x}f+w_{z}x+\frac{w_x{v}^2-w_{y}uv}{f}---\left(6\right)$

式(5)中，v_u表示图像平面x轴方向的光流速度，v_cx表示相机坐标系下x轴方向的速度，v_wx表示世界坐标系下x轴方向的速度，v_wz表示世界坐标系下z轴方向的速度，w_x，w_y，w_z分别表示x,y,z轴方向上的角速度，x_c表示相机坐标系下x轴方向的坐标，y_c表示相机坐标系下y轴方向的坐标；u，v分别表示图像平面x,y轴方向的坐标；

式(6)中，v_v表示图像平面y轴方向的光流速度，v_wy表示世界坐标系下y轴方向的速度；v_wz表示世界坐标系下z轴方向的速度；对垂直高度信息差分可以求得世界坐标系下z轴方向的速度，即

根据式(5)和(6)，以及已知的世界坐标系下z轴方向的速度，可以求得世界坐标系下x轴和y轴方向的飞行器的速度，即为飞行器的水平速度。

进一步，当飞行器定高飞行时，式(5)和(6)简化为：

$v_u=-\frac{v_{wx}f}{z}-w_{y}f+w_{z}v---\left(7\right)$

$v_v=-\frac{v_{wy}f}{z}+w_{x}f-w_{z}u---\left(8\right).$

对所述的飞行器定位方法进行优化，所述的摄像头获取地面的实时图像信息为规则的网格信息，包括如下步骤：

i)利用摄像头，获取地面的实时图像信息，进一步处理实时图像信息中的规则的网格信息，使用霍夫变换算法，调整筛选直线的阈值，对其中的直线信息进行准确提取，从而获取实时图像信息中的直线交点坐标；

ii)同步惯性传感单元、超声波传感器和摄像头的时间戳，对飞行器的水平速度积分，计算飞行器的初始相对位移；给定飞行器的初始位置估计；利用惯性传感单元给定实时姿态信息，给定飞行器的初始姿态信息估计；使用摄像头的投影模型和实时图像信息中获取的直线交点坐标，获取地面网格交点的坐标；

iii)根据地面网格交点相邻距离恒定的条件，对求得的地面网格交点的坐标进行求精操作，获取准确的地面网格交点的坐标，再利用实时图像信息对应的坐标，建立匹配点对，使用直接线性变换算法或者迭代算法，求得飞行器准确的位置信息及姿态信息。

进一步，所述的步骤ii)中获取地面网格交点的实际坐标的方法为：

根据实时图像信息中的粗略位置估计，公式如下：

$p_k^w=p_{k-1}^w+v_w{d}_t---\left(9\right)$

式(9)中，v_w为水平光流速度，为前一时刻水平位置，为后一时刻水平位置，d_t为间隔时间，另外令超声波传感器在k-1时刻和k时刻给入的垂直高度信息分别为：和给定飞行器的初始位置估计，其中和均为三维列向量，表示x,y,z轴方向的位置；

由惯性传感单元给定飞行器实时的姿态信息，包括偏航角、横滚角和俯仰角，如下：

roll＝γ,pitch＝θ,yaw＝ψ (10)

式(10)中，roll、pitch、yaw分别为该时刻的横滚角、俯仰角和偏航角；

对于步骤i)中获取的实时图像信息中的直线交点坐标(u,v)，都满足直线方程：

au+bv+c＝0 (11)；

摄像头的投影模型由以下公式表达：

式(12)中，λ表示相机坐标系下特征点的z轴坐标，(u,v,1)^T表示实时图像信息中特征点的齐次坐标，K表示摄像头的内参矩阵，f_x,f_y,c_x,c_y分别表示x,y轴方向的像素焦距以及主点的像素坐标；表示世界坐标系相对于相机坐标系的转换矩阵，即外参矩阵，前三行三列矩阵表示旋转矩阵，最后一列矩阵表示平移向量；(x_w,y_w,z_w,1)^T表示世界坐标系下特征点的齐次坐标；

将直线方程公式(11)代入摄像头的投影模型公式(12)，得到公式如下：

给定飞行器的初始姿态信息估计和初始位置估计，如下：

$\left(\begin{array}{ccc}{r}_{11}& r_{12}& r_{13}\\ r_{21}& r_{\mathrm{22}}& r_{\mathrm{23}}\\ r_{31}& r_{\mathrm{32}}& r_{\mathrm{33}}\end{array}\right)=\left(\begin{array}{ccc}\cos \mathrm{\γ}\cos \mathrm{\ψ}-\sin \mathrm{\γ}\sin \mathrm{\θ}\sin \mathrm{\ψ}& -\cos \mathrm{\θ}\sin \mathrm{\ψ}& \sin \mathrm{\γ}\cos \mathrm{\ψ}+\cos \mathrm{\γ}\sin \mathrm{\θ}\sin \mathrm{\ψ}\\ \cos \mathrm{\γ}\sin \mathrm{\ψ}+\sin \mathrm{\γ}\sin \mathrm{\θ}\cos \mathrm{\ψ}& \cos \mathrm{\θ}\cos \mathrm{\ψ}& \sin \mathrm{\γ}\sin \mathrm{\ψ}-\cos \mathrm{\γ}\sin \mathrm{\θ}\cos \mathrm{\ψ}\\ -\sin \mathrm{\γ}\cos \mathrm{\θ}& \sin \mathrm{\θ}& \cos \mathrm{\θ}\cos \mathrm{\γ}\end{array}\right)---\left(14\right)$

$\left(\begin{array}{c}{t}_x\\ t_y\\ t_z\end{array}\right)=\left(\begin{array}{c}p_k^w{}_x\\ p_k^w{}_y\\ p_k^w{}_z\end{array}\right)---\left(15\right)$

将给定的参数代入公式(13)，求解公式(13)的直线方程，设求出的世界坐标系的坐标点p_w为：

p_w＝(x_w',y_w',z_w') (16)

对世界坐标系中的直线交点(x_w,y_w,z_w)，同样建立直线方程：

a'x_w+b'y_w+c'z_w+d'＝0 (17)

利用规则的网格信息的特征，直线只分为两类，两类直线互相垂直，且认为地面为平面，则z_w＝0，将直线方程分为两类：

$\left\{\begin{array}{c}{y}_w=-\frac{d^{\′}}{b^{\′}},if{a}^{\′}=0\\ x_w=-\frac{d^{\′}}{a^{\′}},if{b}^{\′}=0\end{array}\right.---\left(18\right)$

由公式(18)对世界坐标系的坐标点进行求精操作，即水平和垂直的两类直线分布间隔均匀，地面的直线实际交点坐标为常值，以此可以对p_w＝(x_w',y_w',z_w')求出的坐标中的x_w'，y_w'取与之相差绝对值最小的常值x_w″,y_w″，z_w'认为是0；求精后的世界坐标系下的直线交点的坐标值为p_w＝(x_w″,y_w″,0)；

利用公式(13)代入求精后p_w的坐标，如下：

简化后得到：

$\begin{array}{c}\[{\mathrm{af}}_x{r}_{11}+{\mathrm{bf}}_y{r}_{21}+({\mathrm{ac}}_x+{\mathrm{bc}}_y+c)r_{31}\]{x_w}^{\′\′}+\[{\mathrm{af}}_x{r}_{12}+{\mathrm{bf}}_y{r}_{22}+({\mathrm{ac}}_x+{\mathrm{bc}}_y+c)r_{32}\]{y_w}^{\′\′}\\ +\[{\mathrm{af}}_x{t}_x+{\mathrm{bf}}_y{t}_y+({\mathrm{ac}}_x+{\mathrm{bc}}_y+c)t_z\]=0\end{array}---\left(20\right)$

由公式(18)知直线方程有两种，故将方程(19)分为两类：

$\begin{array}{c}\[{\mathrm{af}}_x{r}_{11}+{\mathrm{bf}}_y{r}_{21}+({\mathrm{ac}}_x+{\mathrm{bc}}_y+c)r_{31}\]=0\\ \[{\mathrm{af}}_x{r}_{12}+{\mathrm{bf}}_y{r}_{22}+({\mathrm{ac}}_x+{\mathrm{bc}}_y+c)r_{32}\]{y_w}^{\′\′}+\[{\mathrm{af}}_x{t}_x+{\mathrm{bf}}_y{t}_y+({\mathrm{ac}}_x+{\mathrm{bc}}_y+c)t_z\]=0\\ \[{\mathrm{af}}_x{r}_{12}+{\mathrm{bf}}_y{r}_{22}+({\mathrm{ac}}_x+{\mathrm{bc}}_y+c)r_{32}\]=0\\ \[{\mathrm{af}}_x{r}_{11}+{\mathrm{bf}}_y{r}_{21}+({\mathrm{ac}}_x+{\mathrm{bc}}_y+c)r_{31}\]{x_w}^{\′\′}+\[{\mathrm{af}}_x{t}_x+{\mathrm{bf}}_y{t}_y+({\mathrm{ac}}_x+{\mathrm{bc}}_y+c)t_z\]=0\end{array}---\left(21\right)$

提取未知数r₁₁,r₂₁,r₃₁,r₁₂,r₂₂,r₃₂,t_x,t_y,t_z，重写方程(13)，有：

$\left(\begin{array}{ccccccccc}{\mathrm{af}}_x& {\mathrm{bf}}_y& {\mathrm{ac}}_x+{\mathrm{bc}}_y+c& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& {\mathrm{af}}_x{y_w}^{\′\′}& {\mathrm{bf}}_y{y_w}^{\′\′}& ({\mathrm{ac}}_x+{\mathrm{bc}}_y+c){y_w}^{\′\′}& {\mathrm{af}}_x& {\mathrm{bf}}_y& {\mathrm{ac}}_x+{\mathrm{bc}}_y+c\\ 0& 0& 0& {\mathrm{af}}_x& {\mathrm{bf}}_y& {\mathrm{ac}}_x+{\mathrm{bc}}_y+c& 0& 0& 0\\ {\mathrm{af}}_x{x_w}^{\′\′}& {\mathrm{bf}}_y{x_w}^{\′\′}& ({\mathrm{ac}}_x+{\mathrm{bc}}_y+c){x_w}^{\′\′}& 0& 0& 0& {\mathrm{af}}_x& {\mathrm{bf}}_y& {\mathrm{ac}}_x+{\mathrm{bc}}_y+c\end{array}\right)\left(\begin{array}{c}{r}_{11}\\ r_{21}\\ r_{31}\\ r_{12}\\ r_{\mathrm{22}}\\ r_{\mathrm{32}}\\ t_x\\ t_y\\ t_z\end{array}\right)=0---\left(22\right);$

方程(22)是在取两类直线中的各一条建立的方程组，由于存在9个未知数，至少需要五条直线才能解该线性方程组；但直线数超过5条后，该方程组为超定方程组，对其进行奇异值分解即可解该方程组；

至此，解出旋转矩阵R和平移向量t，如下：

$R=\left(\begin{array}{ccc}{r}_{11},& r_{12}& r_{13}\\ r_{21}& r_{22}& r_{23}\\ r_{31}& r_{32}& r_{33}\end{array}\right),t=\left(\begin{array}{c}{t}_x\\ t_y\\ t_z\end{array}\right)---\left(23\right)$

在已求出r₁₁,r₂₁,r₃₁,r₁₂,r₂₂,r₃₂的情况下，可利用旋转矩阵是单位正交矩阵这一特性，求出r₁₃,r₂₃,r₃₃；

再将旋转矩阵转换为姿态角，得到姿态信息：

roll＝a tan2(r₃₂,r₃₃)

$pitch=atan2(-r_{31},\sqrt{{r_{32}}^2+{r_{33}}^2})---\left(24\right).$

yaw＝a tan2(r₂₁,r₁₁)

对所述的飞行器定位方法进一步优化，采用光流速度信息，网格定位的位置信息、惯性传感单元的加速度信息，建立状态方程和观测方程，对视觉算法和惯性传感单元获取的结果使用扩展卡尔曼滤波算法，进行融合处理，获取飞行器准确的位置及姿态信息。

所述的获取飞行器准确的位置及姿态信息具体的算法推导如下：

非线性连续状态方程如下：

x_k+1＝f(x_k,w_k)

对其进行离散化操作，状态量包括有网格定位得到的位置信息p，光流法提供的光流速度信息v和惯性传感单元提供的加速度a；Δt表示连续两次采样的间隔时间，A表示状态转移矩阵；

观测方程如下，分为视觉算法部分的观测方程，观测量位置p；惯性传感单元的观测方程，观测量取加速度a：

$\begin{array}{c}{z}_{v,k}=H_{v,k}{x}_k=\left[\begin{array}{ccc}I& 0& 0\end{array}\right]x_k\\ z_{i,k}=H_{i,k}{x}_k=\left[\begin{array}{ccc}0& 0& I\end{array}\right]x_k\end{array}---\left(26\right)$

其中公式(26)中，z_v,k，z_i,k分别是视觉算法和惯性传感单元的观测量，分别是世界坐标系下的水平方向的位置信息和加速度信息；

根据卡尔曼滤波算法的公式：

$\begin{array}{c}{P}_k^{-}={\mathrm{AP}}_{k-1}{A}^T+Q\\ K_{v,k}=P_k^{-}{H}_{v,k}^T{(H_{v,k}{P}_k^{-}{H}_{i,k}^T+R_v)}^{-1}\\ x_k=x_k^{-}+K_{v,k}(x_{vision}-H_{v,k}{x}_k^{-})\\ P_k=(I-K_{v,k}{H}_{v,k})P_k^{-}\end{array}---\left(27\right)$

$\begin{array}{c}{P}_k^{-}={\mathrm{AP}}_{k-1}{A}^T+Q\\ K_{i,k}=P_k^{-}{H}_{i,k}^T{(H_{i,k}{P}_k^{-}{H}_{i,k}^T+R_i)}^{-1}\\ x_k=x_k^{-}+K_{i,k}(x_{IMU}-H_{i,k}{x}_k^{-})\\ P_k=(I-K_{i,k}{H}_{i,k})P_k^{-}\end{array}---\left(28\right)$

可分别计算视觉算法部分和惯性传感单元部分的卡尔曼增益K_v,k，K_i,k，其中P_k-1表示对应于状态x_k-1的协方差，表示对应于状态的协方差，Q表示系统噪声的协方差矩阵；R_v，R_i分别表示视觉算法和惯性传感单元噪声的协方差矩阵；再结合预测值和测量值，分别求取视觉部分和惯性传感单元部分状态量的最优估值x_k，其中x_vision，x_IMU分别表示视觉算法的测量值和惯性传感单元的测量值；最后再次更新相对于状态量x_k的协方差矩阵。

Claims (8)

1.一种基于定位传感装置的飞行器定位方法，其特征在于，包括如下步骤：

1)通过惯性传感单元中的加速度计和陀螺仪分别测得飞行器实时的加速度和角速度信息，通过磁力计确定方位，获得飞行器的实时姿态信息；所述的实时姿态信息包括偏航角、横滚角和俯仰角；

利用超声波传感器收集飞行器相对于地面的高度信息，并利用惯性传感单元获取的实时姿态信息进行姿态补偿，实时校正垂直高度信息；

利用摄像头，获取地面的实时图像信息；

2)对获取的实时图像信息作灰度化处理，获取灰度图，对每帧灰度图选择固定间隔固定数量的像素块，通过匹配相邻两帧灰度图，计算光流向量，并通过差分，获取光流速度；

3)同步惯性传感单元、超声波传感器和摄像头的时间戳，利用惯性传感单元获取的实时姿态信息，对光流速度进行姿态补偿，得到光流速度的水平面分量；利用超声波传感器获取的垂直高度信息对光流速度的水平面分量进行尺度还原，得到飞行器的水平速度；

4)对飞行器的水平速度积分得到相对位移信息，累加相对位移信息，即得水平的位置信息，结合超声波传感器得到的垂直高度信息，即可得到飞行器的位置信息。

2.根据权利要求1所述的基于定位传感装置的飞行器定位方法，其特征在于，所述的步骤2)中光流速度的获取方法为：

对每帧灰度图选择固定间隔固定数量的像素块，将前一帧灰度图中的像素块位置作为后一帧灰度图的像素块的初始位置，并建立一个搜索区域，阈值函数的公式如下：

$\underset{i=0}{\overset{N-1}{\Σ}}\underset{j=0}{\overset{N-1}{\Σ}}|G_1(i,j)-G_0(i+u,j+v)|---\left(1\right)$

式(1)中，u表示x轴方向的搜索像素个数，v表示y轴方向的搜索像素个数；G₀和G₁分别表示前一帧和后一帧像素块中某一个像素点的灰度值；i，j分别表示图像平面像素块中某一个像素点的x,y轴方向的坐标，N表示像素块的x轴或者y轴方向上的像素个数，G₁(i,j)表示前一帧像素块中每一像素点的灰度值，G₀(i+u,j+v)表示将前一帧图像中的像素块投射到后一帧图像中相同位置后并沿x轴、y轴方向分别移动u、v像素的像素块中每一像素点的灰度值；

通过令前一帧的像素块在后一阵的灰度图的搜索区域内搜索灰度差和最小，求得后一帧的像素块位置，两帧像素块的距离差，即为光流向量，通过差分，获取光流速度。

3.根据权利要求2所述的基于定位传感装置的飞行器定位方法，其特征在于，所述的步骤3)中获得飞行器的水平速度的方法为：

建立图像平面中的点与摄像机坐标系中的点的转换关系：

$u=f\frac{p_c}{z}---\left(2\right)$

式(2)中，u为图像平面中的点的坐标，p_c为摄像机坐标系的点的坐标，f为焦距，z为垂直高度；

对式(2)两边求导，即得光流速度：

$\frac{flow}{\mathrm{\Δ}time}=v=f\frac{{\mathrm{zv}}_c-v_z{p}_c}{z^2}---\left(3\right)$

式(3)中，v_c为摄像头坐标系下点的运动速度，v_z为垂直方向上的运动速度，v为图像平面的运动速度，即光流速度；

利用刚体运动学的公式：

v_c＝-v_w-ω×p_c (4)

式(4)中，v_w表示世界坐标系下点的运动速度，ω表示角速度；

将公式(4)带入公式(3)，展开可得到x轴和y轴方向的光流速度与世界坐标系下的飞行器速度的关系：

$\begin{array}{c}{v}_u=f\frac{{\mathrm{zv}}_{cx}-v_z{x}_c}{z^2}=f\frac{z(-v_{wx}-w_{y}z+w_z{y}_c)-(-v_{wz}-w_x{y}_c+w_y{x}_c)x_c}{z^2}---\left(5\right)\\ =\frac{v_{wx}x-v_{wx}f}{z}-w_{y}f+w_{z}y+\frac{w_{x}uv-w_y{u}^2}{f}\end{array}$

$v_v=\frac{v_{wz}y-v_{wy}f}{z}+w_{x}f+w_{z}x+\frac{w_x{v}^2-w_{y}uv}{f}---\left(6\right)$

式(5)中，v_u表示图像平面x轴方向的光流速度，v_cx表示相机坐标系下x轴方向的速度，v_wx表示世界坐标系下x轴方向的速度，v_wz表示世界坐标系下z轴方向的速度，w_x，w_y，w_z分别表示x,y,z轴方向上的角速度，x_c表示相机坐标系下x轴方向的坐标，y_c表示相机坐标系下y轴方向的坐标；u，v分别表示图像平面x,y轴方向的坐标；

式(6)中，v_v表示图像平面y轴方向的光流速度，v_wy表示世界坐标系下y轴方向的速度；v_wz表示世界坐标系下z轴方向的速度；对垂直高度信息差分可以求得世界坐标系下z轴方向的速度，即

根据式(5)和(6)以及已知的世界坐标系下z轴方向的速度，可以求得世界坐标系下x轴和y轴方向的飞行器的速度，即为飞行器的水平速度。

4.根据权利要求3所述的基于定位传感装置的飞行器定位方法，其特征在于，当飞行器定高飞行时，式(5)和(6)简化为：

$v_u=-\frac{v_{wx}f}{z}-w_{y}f+w_{z}v---\left(7\right)$

$v_v=-\frac{v_{wy}f}{z}+w_{x}f-w_{z}u---\left(8\right).$

5.根据权利要求1所述的基于定位传感装置的飞行器定位方法，其特征在于，对所述的飞行器定位方法进行优化，所述的摄像头获取地面的实时图像信息为规则的网格信息，包括如下步骤：

i)利用摄像头，获取地面的实时图像信息，进一步处理实时图像信息中的规则的网格信息，使用霍夫变换算法，调整筛选直线的阈值，对其中的直线信息进行准确提取，从而获取实时图像信息中的直线交点坐标；

ii)同步惯性传感单元、超声波传感器和摄像头的时间戳，对飞行器的水平速度积分，计算飞行器的初始相对位移；给定飞行器的初始位置估计；利用惯性传感单元给定实时姿态信息，给定飞行器的初始姿态信息估计；使用摄像头的投影模型和实时图像信息中获取的直线交点坐标，获取地面网格交点的坐标；

iii)根据地面网格交点相邻距离恒定的条件，对求得的地面网格交点的坐标进行求精操作，获取准确的地面网格交点的坐标，再利用实时图像信息对应的坐标，建立匹配点对，使用直接线性变换算法或者迭代算法，求得飞行器准确的位置信息及姿态信息。

6.根据权利要求5所述的基于定位传感装置的飞行器定位方法，其特征在于，所述的步骤ii)中获取地面网格交点的实际坐标的方法为：

根据实时图像信息中的粗略位置估计，公式如下：

$p_k^w=p_{k-1}^w+v_w{d}_t---\left(9\right)$

式(9)中，v_w为水平光流速度，为前一时刻水平位置，为后一时刻水平位置，d_t为间隔时间，另外令超声波传感器在k-1时刻和k时刻给入的垂直高度信息分别为：和给定飞行器的初始位置估计，其中和均为三维列向量，表示x,y,z轴方向的位置；

由惯性传感单元给定飞行器实时的姿态信息，包括偏航角、横滚角和俯仰角，如下：

roll＝γ,pitch＝θ,yaw＝ψ (10)

式(10)中，roll、pitch、yaw分别为该时刻的横滚角、俯仰角和偏航角；

对于步骤i)中获取的实时图像信息中的直线交点坐标(u,v)，都满足直线方程：

au+bv+c＝0 (11)；

摄像头的投影模型由以下公式表达：

式(12)中，λ表示相机坐标系下特征点的z轴坐标，(u,v,1)^T表示实时图像信息中特征点的齐次坐标，K表示摄像头的内参矩阵，f_x,f_y,c_x,c_y分别表示x,y轴方向的像素焦距以及主点的像素坐标；表示世界坐标系相对于相机坐标系的转换矩阵，即外参矩阵，前三行三列矩阵表示旋转矩阵，最后一列矩阵表示平移向量；(x_w,y_w,z_w,1)^T表示世界坐标系下特征点的齐次坐标；

将直线方程公式(11)代入摄像头的投影模型公式(12)，得到公式如下：

给定飞行器的初始姿态信息估计和初始位置估计，如下：

$\left(\begin{array}{ccc}{r}_{11}& r_{12}& r_{13}\\ r_{21}& r_{\mathrm{22}}& r_{\mathrm{23}}\\ r_{31}& r_{\mathrm{32}}& r_{\mathrm{33}}\end{array}\right)=\left(\begin{array}{ccc}\cos \mathrm{\γ}\cos \mathrm{\ψ}-\sin \mathrm{\γ}\sin \mathrm{\θ}\sin \mathrm{\ψ}& -\cos \mathrm{\θ}\sin \mathrm{\ψ}& \sin \mathrm{\γ}\cos \mathrm{\ψ}+\cos \mathrm{\γ}\sin \mathrm{\θ}\sin \mathrm{\ψ}\\ \cos \mathrm{\γ}\sin \mathrm{\ψ}+\sin \mathrm{\γ}\sin \mathrm{\θ}\cos \mathrm{\ψ}& \cos \mathrm{\θ}\cos \mathrm{\ψ}& \sin \mathrm{\γ}\sin \mathrm{\ψ}-\cos \mathrm{\γ}\sin \mathrm{\θ}\cos \mathrm{\ψ}\\ -\sin \mathrm{\γ}\cos \mathrm{\θ}& \sin \mathrm{\θ}& \cos \mathrm{\θ}\cos \mathrm{\γ}\end{array}\right)---\left(14\right)$

$\left(\begin{array}{c}{t}_x\\ t_y\\ t_z\end{array}\right)=\left(\begin{array}{c}p_k^w{}_x\\ p_k^w{}_y\\ p_k^w{}_z\end{array}\right)---\left(15\right)$

将给定的参数代入公式(13)，求解公式(13)的直线方程，设求出的世界坐标系的坐标点p_w为：

p_w＝(x_w',y_w',z_w') (16)

对世界坐标系中的直线交点(x_w,y_w,z_w)，同样建立直线方程：

a'x_w+b'y_w+c'z_w+d'＝0 (17)

利用规则的网格信息的特征，直线只分为两类，两类直线互相垂直，且认为地面为平面，则z_w＝0，将直线方程分为两类：

$\left\{\begin{array}{c}{y}_w=-\frac{d^{\′}}{b^{\′}},if{a}^{\′}=0\\ x_w=-\frac{d^{\′}}{a^{\′}},if{b}^{\′}=0\end{array}\right.---\left(18\right)$

由公式(18)对世界坐标系的坐标点进行求精操作，即水平和垂直的两类直线分布间隔均匀，地面的直线实际交点坐标为常值，以此可以对p_w＝(x_w',y_w',z_w')求出的坐标中的x_w'，y_w'取与之相差绝对值最小的常值x_w”,y_w”，z_w'认为是0；求精后的世界坐标系下的直线交点的坐标值为p_w＝(x_w”,y_w”,0)；

利用公式(13)代入求精后p_w的坐标，如下：

简化后得到：

$\begin{array}{c}\[{\mathrm{af}}_x{r}_{11}+{\mathrm{bf}}_y{r}_{21}+({\mathrm{ac}}_x+{\mathrm{bc}}_y+c)r_{31}\]{x_w}^{\′\′}+\[{\mathrm{af}}_x{r}_{12}+{\mathrm{bf}}_y{r}_{22}+({\mathrm{ac}}_x+{\mathrm{bc}}_y+c)r_{32}\]{y_w}^{\′\′}\\ +\[{\mathrm{af}}_x{t}_x+{\mathrm{bf}}_y{t}_y+({\mathrm{ac}}_x+{\mathrm{bc}}_y+c)t_z\]=0\end{array}---\left(20\right)$

由公式(18)知直线方程有两种，故将方程(19)分为两类：

$\begin{array}{c}\[{\mathrm{af}}_x{r}_{11}+{\mathrm{bf}}_y{r}_{21}+({\mathrm{ac}}_x+{\mathrm{bc}}_y+c)r_{31}\]=0\\ \[{\mathrm{af}}_x{r}_{12}+{\mathrm{bf}}_y{r}_{22}+({\mathrm{ac}}_x+{\mathrm{bc}}_y+c)r_{32}\]{y_w}^{\′\′}+\[{\mathrm{af}}_x{t}_x+{\mathrm{bf}}_y{t}_y+({\mathrm{ac}}_x+{\mathrm{bc}}_y+c)t_z\]=0\\ \[{\mathrm{af}}_x{r}_{12}+{\mathrm{bf}}_y{r}_{22}+({\mathrm{ac}}_x+{\mathrm{bc}}_y+c)r_{32}\]=0\\ \[{\mathrm{af}}_x{r}_{11}+{\mathrm{bf}}_y{r}_{21}+({\mathrm{ac}}_x+{\mathrm{bc}}_y+c)r_{31}\]{x_w}^{\′\′}+\[{\mathrm{af}}_x{t}_x+{\mathrm{bf}}_y{t}_y+({\mathrm{ac}}_x+{\mathrm{bc}}_y+c)t_z\]=0\end{array}---\left(21\right)$

提取未知数r₁₁,r₂₁,r₃₁,r₁₂,r₂₂,r₃₂,t_x,t_y,t_z，重写方程(13)，有：

$\left(\begin{array}{ccccccccc}{\mathrm{af}}_x& {\mathrm{bf}}_y& {\mathrm{ac}}_x+{\mathrm{bc}}_y+c& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& {\mathrm{af}}_x{x_w}^{\′\′}& {\mathrm{bf}}_y{y_w}^{\′\′}& ({\mathrm{ac}}_x+{\mathrm{bc}}_y+c){x_w}^{\′\′}& {\mathrm{af}}_x& {\mathrm{bf}}_y& {\mathrm{ac}}_x+{\mathrm{bc}}_y+c\\ 0& 0& 0& {\mathrm{af}}_x& {\mathrm{bf}}_y& {\mathrm{ac}}_x+{\mathrm{bc}}_y+c& 0& 0& 0\\ {\mathrm{af}}_x{x_w}^{\′\′}& {\mathrm{bf}}_y{x_w}^{\′\′}& ({\mathrm{ac}}_x+{\mathrm{bc}}_y+c){x_w}^{\′\′}& 0& 0& 0& {\mathrm{af}}_x& {\mathrm{bf}}_y& {\mathrm{ac}}_x+{\mathrm{bc}}_y+c\end{array}\right)\left(\begin{array}{c}{r}_{11}\\ r_{21}\\ r_{31}\\ r_{12}\\ r_{22}\\ r_{32}\\ t_x\\ t_y\\ t_z\end{array}\right)=0---\left(22\right);$

方程(22)是在取两类直线中的各一条建立的方程组，由于存在9个未知数，至少需要五条直线才能解该线性方程组；但直线数超过5条后，该方程组为超定方程组，对其进行奇异值分解即可解该方程组；

至此，解出旋转矩阵R和平移向量t，如下：

$R=\left(\begin{array}{ccc}{r}_{11},& r_{12}& r_{13}\\ r_{21}& r_{22}& r_{23}\\ r_{31}& r_{32}& r_{33}\end{array}\right),t=\left(\begin{array}{c}{t}_x\\ t_y\\ t_z\end{array}\right)---\left(23\right)$

在已求出r₁₁,r₂₁,r₃₁,r₁₂,r₂₂,r₃₂的情况下，可利用旋转矩阵是单位正交矩阵这一特性，求出r₁₃,r₂₃,r₃₃；

再将旋转矩阵转换为姿态角，得到姿态信息：

$\begin{array}{c}roll=a\tan 2(r_{32},r_{33})\\ pitch=a\tan 2(-r_{31},\sqrt{{r_{32}}^2+{r_{33}}^2})\\ yaw=a\tan 2(r_{21},r_{11})\end{array}---\left(24\right).$

7.根据权利要求6所述的基于定位传感装置的飞行器定位方法，其特征在于，对所述的飞行器定位方法进一步优化，采用光流速度信息，网格定位的位置信息、惯性传感单元的加速度信息，建立状态方程和观测方程，对视觉算法和惯性传感单元获取的结果使用扩展卡尔曼滤波算法，进行融合处理，获取飞行器准确的位置及姿态信息。

8.跟据权利要求7所述的基于定位传感装置的飞行器定位方法，其特征在于，所述的获取飞行器准确的位置及姿态信息具体的算法推导如下：

非线性连续状态方程如下：

对其进行离散化操作，状态量包括有网格定位得到的位置信息p，光流法提供的光流速度信息v和惯性传感单元提供的加速度a；Δt表示连续两次采样的间隔时间，A表示状态转移矩阵；

观测方程如下，分为视觉算法部分的观测方程，观测量位置p；惯性传感单元的观测方程，观测量取加速度a：

$\begin{array}{c}{z}_{v,k}=H_{v,k}{x}_k=\left[\begin{array}{ccc}I& 0& 0\end{array}\right]x_k\\ z_{i,k}=H_{i,k}{x}_k=\left[\begin{array}{ccc}0& 0& I\end{array}\right]x_k\end{array}---\left(26\right)$

其中公式(26)中，z_v,k，z_i,k分别是视觉算法和惯性传感单元的观测量，分别是世界坐标系下的水平方向的位置信息和加速度信息；

根据卡尔曼滤波算法的公式：

$\begin{array}{c}{P}_k^{-}={\mathrm{AP}}_{k-1}{A}^T+Q\\ K_{v,k}=P_k^{-}{H}_{v,k}^T{(H_{v,k}{P}_k^{-}{H}_{i,k}^T+R_v)}^{-1}\\ x_k=x_k^{-}+K_{v,k}(x_{vision}-H_{v,k}{x}_k^{-})\\ P_k=(1-K_{v,k}{H}_{v,k})P_k^{-}\end{array}---\left(27\right)$

$\begin{array}{c}{P}_k^{-}={\mathrm{AP}}_{k-1}{A}^T+Q\\ K_{i,k}=P_k^{-}{H}_{i,k}^T{(H_{i,k}{P}_k^{-}{H}_{i,k}^T+R_i)}^{-1}\\ x_k=x_k^{-}+K_{i,k}(x_{IMU}-H_{i,k}{x}_k^{-})\\ P_k=(1-K_{i,k}{H}_{i,k})P_k^{-}\end{array}---\left(28\right)$

可分别计算视觉算法部分和惯性传感单元部分的卡尔曼增益K_v,k，K_i,k，其中P_k-1表示对应于状态x_k-1的协方差，表示对应于状态的协方差，Q表示系统噪声的协方差矩阵，R_v，R_i分别表示视觉算法和惯性传感单元噪声的协方差矩阵；再结合预测值和测量值，分别求取视觉部分和惯性传感单元部分状态量的最优估值x_k，其中x_vision，x_IMU分别表示视觉算法的测量值和惯性传感单元的测量值；最后再次更新相对于状态量x_k的协方差矩阵。