一种无人机在压力管道中心轴线的速度估计方法

技术领域

本发明涉及无人机速度估计技术领域，尤其是涉及一种无人机在压力管道中心轴线的速度估计方法。

背景技术

压力管道是水电站系统的重要组成部分。压力管道用于将水库的水输送到水轮机中发电。压力管道是由钢和钢筋混凝土制成。水电站在运行一段时间后，压力管道内壁的状况就是未知的。因此需要对压力管道定期进行检查，看是否有损坏的情况，如开裂，腐蚀或者气蚀。人工的监测方法由于压力管道内由于可能存在渗水，混凝土掉块等安全隐患,对检查人员的人身安全会造成威胁。基于无人机的视觉检查方法是一个安全，高效的解决方式。它也解决了人员难以到达管道的斜直段，上弯段等检查部位的难题。

压力管道可以看成是一个简单的垂直截面的沿弯曲的轴线拉伸结构，如图2所示。由于在压力管道内，GPS信号缺失且地磁场受到干扰。因此，基于无人机的视觉检查方法必须要解决管道中的无人机定位导航问题。目前无人机在GPS-denied环境的定位导航方法主要有依赖于UWB、光流、视觉相机、3D激光雷达等传感器或者采用UWB-视觉-激光雷达-IMU多传感器融合定位的方式。

视觉SLAM随着以SIFT算子为代表的图像局部不变特征描述方法的兴起，图像特征点的提取与匹配能适应一定程度的图像形变及光照变化。SIFT算法提高了视觉SLAM在复杂环境中的适用性，但是它较大的计算量限制了提取效率，难以满足实时性要求。研究人员为了提高提取效率陆续开发了SURF、ORB等算子。这些算子尽管性能有不同程度降低，但是处理效率得到了数倍至数十倍的提升，使得实时视觉SLAM成为现实。目前基于视觉的定位研究许多仍然是在纹理非常良好的室内空间进行的。压力管道内暗光且重复的纹理给基于视觉的定位方式带来了很大程度的挑战。

压力管道内部结构较为单一，3D激光雷达沿中心轴线扫描的前后两帧点云信息均是管道壁面结构相似点云,缺乏几何特征。Filip等人在a featureless tunnelenvironment以F-LOAM、LeGO-LOAM、LIO-SAM等为代表的激光SLAM算法进行测试(state-ofthe-art lidar SLAM algorithms)。在隧道环境中，由于缺乏大量的边缘特征，上述SLAM算法性能显著下降。其中，LIO-SAM、LeGo-LOAM、F-LOAM算法的轨迹Root Mean SquareError(RMSE)分别为8.74m、9.12m、9.53m。在隧道中放置了大量具有很强的边缘特征的地标之后，LIO-SAM算法的RMSE显著降低到0.95m，但是LeGo-LOAM、F-LOAM算法的轨迹RMSE没有明显降低，它们分别为6.87m、4.34m。3D激光SLAM在压力管道内沿中心轴线的定位导航面临着缺乏大量显著的边缘结构特征的挑战。

综上所述，UWB基站是较难在斜直段或者上弯段等位置进行布设的，UWB定位方法还存在壁面反射干扰的问题，当无人机距离UWB基站过远或较为接近壁面时，UWB就不能良好的工作，甚至出现大幅度的信号突变或丢失。

发明内容

本发明的目的就是为了克服上述现有技术存在的缺陷而提供的一种无人机在压力管道中心轴线的速度估计方法，挖掘模糊图像中暗藏的速度信息，克服壁面反射干扰和速度信息丢失的问题。

本发明的目的可以通过以下技术方案来实现：

一种无人机在压力管道中心轴线的速度估计方法，包括以下步骤：

S1、基于无人机获取在压力管道的图像，提取图像中连续的运动模糊图像；

S2、提取运动模糊图像对应的频谱图，判断频谱图是否符合频谱图条件，若符合，则执行S3，反之则：

将图像中满足线性运动模型的n帧图像进行合成，合成得到的运动模糊图像为n帧图像根据所有线段的离散样本的加权平均值，所述线段与n帧图像的每个像素一一对应，然后执行S3；

所述频谱图条件具体为：频谱图的中心为亮条纹，左右两侧为周期性等间距明暗条纹，且明暗条纹为平行的直线；所述频谱图条件通过建立退化图像模型，将频谱图二维离散化，结合退化图像模型的频域关系式和运动模糊函数的傅里叶变换，得到的运动模糊图像的频谱关系式分析得到；

S3、将运动模糊图像和其频谱图裁剪为正方形，提取裁剪后的运动模糊图像的频谱图的第一暗条纹倾角，得到第一运动模糊方向，所述第一运动模糊方向为无人机沿轴向的第一运动方向；

S4、从运动模糊图像的频谱图中提取中心亮条纹旁的两条暗条纹间距，得到模糊运动长度；

S5、基于模糊运动长度，根据相机成像原理计算在曝光时间范围内实际的运动距离，将在曝光时间范围内实际的运动距离除以相机的非曝光时间，得到无人机沿轴向的速度值，综合无人机沿轴向的第一运动方向和速度值得到估计的速度。

进一步地，提取中心亮条纹旁的两条暗条纹间距的具体步骤为：

将运动模糊图像的频谱图中的像素点进行动态阈值二值化；

基于拉东变换获取第二暗条纹倾角，确定第二运动方向；

将像素点沿第二运动方向进行积分，得到积分图；

对所获得的积分图进行滤波；

计算滤波后的积分图的两个波谷之间的距离，并均值化处理，得到两条暗条纹间距。

进一步地，模糊运动长度的表达式为：

其中，θ为运动模糊方向，M为运动模糊图像的长，N为运动模糊图像的宽，d为暗条纹间距。

进一步地，在曝光时间范围内实际的运动距离的表达式为：

其中，L表示曝光时间t

进一步地，无人机沿轴向的速度值的表达式为：

其中，v

进一步地，运动模糊图像的频谱图为正方形时，运动模糊方向与暗条纹倾角的关系式为：

其中，θ为运动模糊方向，β为暗条纹倾角，M为运动模糊图像的长，N为运动模糊图像的宽，M和N相等。

进一步地，暗条纹倾角满足以下关系式：

其中，M为运动模糊图像的长，N为运动模糊图像的宽，a和b为图像在x和y方向上的位移。

进一步地，建立的退化图像模型的表达式为：

g(x,y)＝f(x,y)*h(x,y)+n(x,y)

其中，g(x,y)是退化模糊图像的表征值，f(x,y)是原图像的表征值，h(x,y)是运动模糊函数，n(x,y)是高斯噪声。

进一步地，运动模糊函数的表达式为：

其中，h(x,y)是运动模糊函数，l表示图片的模糊长度，即运动模糊长度，θ为运动模糊方向。

进一步地，退化图像模型转换到频域时，设定n(x,y)＝0，n(x,y)是高斯噪声。

与现有技术相比，本发明具有以下有益效果：

(1)本发明考虑到无人机沿压力管道中心轴线方向上速度信息缺失的问题，相比于传统的视觉SLAM去模糊的方法，根据运动模糊图像的频谱图提取模糊运动长度和运动模糊方向，挖掘模糊图像中暗藏的速度信息，克服壁面反射干扰和速度信息丢失的问题。

(2)本发明针对单帧微小模糊图像无法估计速度的情况，本发明提出将满足线性运动模型的多帧图像进行合成，得到的运动模糊图像，可以较好的将单帧图像上的微小运动模糊合成明显的运动模糊。

附图说明

图1为本发明的流程图；

图2为压力管道模型图；

图3为相机成像关系图；

图4(a)为图片灰度图，图4(b)为图片灰度图对应的频率图；

图5为无人机拍摄的图片示例；

图6为选取的前后间隔30帧的两张测试图片，其中图6(a)为先拍摄的第一帧图片，图6(b)为后拍摄的第二帧图片；

图7为SIFT特征检测和匹配结果图，其中(a)为第一张提取的SIFT特征图，(b)第二张提取的SIFT特征图，(c)为SIFT特征原始匹配结果，(d)为RANSAC消除错误匹配后，剩余的匹配结果；

图8为ORB特征检测和匹配结果图，其中(a)第一张提取的ORB特征图，(b)第二张提取的ORB特征图，(c)ORB特征原始匹配结果，(d)RANSAC消除错误匹配后，剩余的匹配结果；

图9为SURF特征检测和匹配结果图，其中(a)第一张提取的SURF特征图，(b)第二张提取的SURF特征图，(c)SURF特征原始匹配结果，(d)RANSAC消除错误匹配后，剩余的匹配结果；

图10为不同帧数合成的运动模糊频谱图，其中(a)为2帧合成的运动模糊频谱图，(b)为5帧合成的运动模糊频谱图，(c)为10帧合成的运动模糊频谱图，(d)为20帧合成的运动模糊频谱图；

图11为5帧人工合成运动模糊图处理结果，其中(a)五帧合并效果图，(b)运动模糊灰度图，(c)运动模糊频谱图，(d)运动模糊频谱二值化图，(e)运动模糊频谱二值化积分图，(f)二值化积分差分阈值滤波图；

图12为速度对比图。

具体实施方式

下面结合附图和具体实施例对本发明进行详细说明。本实施例以本发明技术方案为前提进行实施，给出了详细的实施方式和具体的操作过程，但本发明的保护范围不限于下述的实施例。

本发明提出一种无人机在压力管道中心轴线的速度估计方法，方法的流程图如图1所示。包括以下步骤：

S1、基于无人机获取在压力管道的图像，提取图像中连续的运动模糊图像；

S2、提取运动模糊图像对应的频谱图，判断频谱图是否符合频谱图条件，若符合，则执行S3，反之则：

将图像中满足线性运动模型的n帧图像进行合成，合成得到的运动模糊图像为n帧图像根据所有线段的离散样本的加权平均值，线段与n帧图像的每个像素一一对应，然后执行S3；

频谱图条件具体为：频谱图的中心为亮条纹，左右两侧为周期性等间距明暗条纹，且明暗条纹为平行的直线；频谱图条件通过建立退化图像模型，将频谱图二维离散化，结合退化图像模型的频域关系式和运动模糊函数的傅里叶变换，得到的运动模糊图像的频谱关系式分析得到；

S3、将运动模糊图像和其频谱图裁剪为正方形，提取裁剪后的运动模糊图像的频谱图的第一暗条纹倾角，得到第一运动模糊方向，第一运动模糊方向为无人机沿轴向的第一运动方向；

S4、从运动模糊图像的频谱图中提取中心亮条纹旁的两条暗条纹间距，得到模糊运动长度；

S5、基于模糊运动长度，根据相机成像原理计算在曝光时间范围内实际的运动距离，将在曝光时间范围内实际的运动距离除以相机的非曝光时间，得到无人机沿轴向的速度值，综合无人机沿轴向的第一运动方向和速度值得到估计的速度。

S2中，频谱图条件的原理如下：

在有关模糊长度的研究中，模糊图和原图的关系通常被表示为公式(4)：

g(x,y)＝f(x,y)*h(x,y)+n(x,y)(4)

其中x，y是图像的长度和宽度的坐标。

g(x,y)是退化模糊图像的表征值，f(x,y)是原图像的表征值，h(x,y)是运动模糊函数，n(x,y)是高斯噪声，*表示两者的卷积。

在瞬时运动模型中,运动模糊函数h(x,y)可以用关系式(5)表示：

其中，θ为运动模糊方向，在曝光时间t

时域卷积是一个耗费时间的过程，因此将时域的卷积转换到频域的乘积来减少运算时间。

设定n(x,y)＝0,则退化图像模型转化为频域可以表示为关系式(6)：

G(u，v)＝F(u，v)·H(u，v)(6)

瞬时运动模糊函数的傅里叶变换可以表示为关系式(7)：

假定图像大小为M*N(长*宽)。将频谱图像二维离散化，联合(6)(7)可以推导出运动模糊图像的频谱关系式(8)：

其中u属于[0,N-1]，v属于[0,M-1]；令

当

S2中，若频谱图不符合频谱图条件，此时需要考虑从连续捕获的两张或者多张不模糊图像来合成运动模糊图。以两张图像的合成为例进行解释。曝光时间[s1,e1]时间段内拍摄的照片记作F

对于图像中每个像素F

其中Fi(p,q)是Fi在图像任意位置(p,q)的双线性插值结果,al,bl分别表示在上一曝光时间t1内,图像在x和y方向的位移。Wi(p,q,k)为每个像素与离散核的卷积，离散核由线预测模型中使用的加权双线性插值核的总和组成。

S3中，设暗条纹倾角为β，倾角满足如下关系式(9)：

运动模糊方向与该倾角存在关系式(10)：

当原图像经过裁剪为正方形时，即M＝N，频谱图像上暗条纹和运动模糊方向相互垂直.提取出频谱图像暗条纹的倾角即可根据关系式(10)计算出无人机沿轴向的运动方向。

S4中，根据傅里叶变换的平移特性，常常将u＝0,v＝0点移到频谱中心位置，直线

此时只要从运动模糊频谱图中提取到中心亮条纹旁的两条暗条纹间距2d,即可由式(12)求出模糊运动长度，再根据模糊运动长度可以计算出实际的运动距离从而推算无人机在该时间的瞬时速度。

按照如下步骤提取中心亮条纹旁的两条暗条纹间距：

1.将像素点进行动态阈值二值化；

2.基于Radon变换计算暗条纹倾角

3.将像素点沿运动方向进行积分；

4.对所获得的积分图进行滤波；

5.计算两个波谷之间的距离，并均值化处理。

S5中，根据得到的模糊长度，根据图3的相机成像原理可以计算得到在曝光时间t

式中：L表示曝光时间t

相机拍摄两帧图片的时间间隔包括曝光时间和非曝光时间两部分。相机的曝光时间t

则在两帧图片之间，实际的位移距离Lr可以表示为：Lr＝L/r。

轴线方向上(默认为x方向)，速度估计可以表示为：

下面进行实验验证：

在水电站的压力管道内进行了飞行测试，由于该压力管道内部开阔，它的直径为12.4m。即使在无人机上安装了六个补光灯，相机拍摄的图像仍然较暗，且部分区域为全黑。所拍摄的图片如图5所示。

从图5可以看出，虽然压力管道壁面存在一些焊缝，但是这些焊缝几乎是相同的。在后面的三种特征描述子提取实验中可以证实这些焊缝的特征很难提取到，即使提取到一定数量的特征，也会产生大量的错误匹配。拍摄的图片中存在无人机脚架的遮挡以及脚架投影在壁面的暗纹，这些都对视觉定位算法的稳定性和精度构成了较大挑战。接下来采用三个state-of-the-art视觉特征描述子(SIFT,SURF,ORB)进行测试。本实验中，选取了相比图5光照更均匀,存在四处微小特征的两张测试图片。这四处微小特征分别标号为1，2，3，4，它们之间存在如图6所示的对应关系。为减少无人机的脚架遮挡区域对算法的影响，预处理时将这部分都裁剪了。图6(a)为先拍摄的第一帧图片，图6(b)为后拍摄的第二帧图片

先采用上述三种视觉特征点检测算法SIFT,SURF,ORB时，提取管道壁面的特征，再对提取到的特征进行匹配，最后根据匹配阈值筛选错误匹配，经过RANSAC消除错误匹配后显示剩余匹配结果。

图7为SIFT特征检测和匹配结果图，其中(a)为第一张提取的SIFT特征图，(b)第二张提取的SIFT特征图，(c)为SIFT特征原始匹配结果，(d)为RANSAC消除错误匹配后，剩余的匹配结果。从图7可以观察到，在管道壁面只能提取到非常少量的SIFT特征点，焊缝处几乎没有提取到特征.对照图6的四处特征对应关系，(c)图中的匹配大部分是错误的，其中1号特征没有匹配，左边的3号特征不仅能匹配到右边的3号特征，还能匹配到4号以及其它处特征。经过筛选错误匹配后，仅剩四处匹配，且仅有2号的特征匹配正确。

图8为ORB特征检测和匹配结果图，其中(a)第一张提取的ORB特征图，(b)第二张提取的ORB特征图，(c)ORB特征原始匹配结果，(d)RANSAC消除错误匹配后，剩余的匹配结果。从图8可以观察到，在管道壁面提取到较多的ORB特征点.对照图6的四处特征对应关系，(c)中出现了大量的错误匹配，其中1，2，4号特征几乎没有对应匹配，3号特征有较多错误匹配。经过筛选错误匹配后，仅有3号的特征匹配正确。在(d)图左边的3号特征点分布较散，而右边的3号特征点聚在一点。

图9为SURF特征检测和匹配结果图，其中(a)第一张提取的SURF特征图，(b)第二张提取的SURF特征图，(c)SURF特征原始匹配结果，(d)RANSAC消除错误匹配后，剩余的匹配结果。从图9可以观察到，四处微小特征均被检测到。(c)图中可以看到特征匹配中存在部分错误匹配，2号特征部分错误匹配到1号和3号特征，3号特征错误匹配到4号特征。对于其它close-to-texture-less区域壁面，SURF较难提取到特征。经过筛选错误匹配后，1号和2号的特征匹配正确。

从上述三种特征描述子提取和匹配实验结果可以看出，压力管道内暗光环境且壁面大范围弱纹理给视觉特征提取和匹配带来了极大挑战。在实际测试中，上述三种视觉特征定位算法都很难达到期望的结果。

下面采用本发明的方法进行速度的估计：

运动模糊频谱图应该是一个中心亮条纹，两侧是明暗等间距条纹。对图4(a)做傅里叶变换转化为运动模糊频谱图，再对其二值化处理。

运动模糊频谱图中只有中心亮条纹，但没有出现预期的等间距明暗条纹，经过二值化处理后显示为团云。这是由于前后帧拍摄时间间隔较短且无人机的速度较慢，运动模糊不明显就导致了团云状的二值化频谱图。

既然无人机拍摄的单张图片的运动模糊不明显，考虑到相机拍摄时的帧率为30fps，可以利用多帧图片的信息来合成为一张运动模糊图片。此时就涉及到选取多少帧图片合成是较为合适的问题。叠加帧数的数量选取，受限于物理设备自身的帧数、最低所需的更新频率以及合并之后的效果呈现。对视频流进行分解，并选择2帧、5帧、10帧、20帧等不同的帧数进行合成运动模糊图。再对该合成运动模糊图做傅里叶变换，实验结果如图10所示。图10为不同帧数合成的运动模糊频谱图，其中(a)为2帧合成的运动模糊频谱图，(b)为5帧合成的运动模糊频谱图，(c)为10帧合成的运动模糊频谱图，(d)为20帧合成的运动模糊频谱图。

从图10可以看出，实际合成帧数的选取，极大的影响到频谱的结果。当帧数过多时，中心亮条纹还比较清晰，但是两侧的明暗条纹已经不明显了。起始帧和终止帧的状态差异可能会比较明显，从而影响到最终频谱图。相反，当帧数过少时，并不能得到明显的明暗交替条纹，条纹范围较宽。所以选择合适的合成帧数，是获得良好结果的保障。本实验采用5帧合成一个人工模糊图像帧。

选取另一张图片，进行5帧人工合成，得到的运动模糊图处理结果如图11所示。其中(a)五帧合并效果图，(b)运动模糊灰度图，(c)运动模糊频谱图，(d)运动模糊频谱二值化图，(e)运动模糊频谱二值化积分图，(f)二值化积分差分阈值滤波图。从合成图可以看出合成的人工模糊是运动模糊，但模糊的运动方向和模糊长度并不是一个确切的结果。对合成多帧运动模糊频谱图分析来看，频谱图呈现出与理论模型情况相近的明暗交替的条纹，符合理想预期，可以进行结果分析和计算。

在压力管道的下水平段内搭建四个UWB基站，在飞机上搭载UWB信标，这样五个UWB可以组成一个四棱锥。即使一个UWB基站信号丢失,它们仍然可以组成一个三棱锥。我们只需要解算出三棱锥或者四棱锥的高，即可求解出无人机在中心轴线的位移与速度。本文将此解算的速度作为真值。基于运动模糊的速度估计与上述真值进行比对，评估本研究提出算法的精度与稳定性，如图12所示。从图12可以看出，两条曲线基本重合，最大误差在5.4％，RMSE为0.025m/s，速度变化趋势近似，本文所提出的算法具有较好的准确性与稳定性。

以上详细描述了本发明的较佳具体实施例。应当理解，本领域的普通技术人员无需创造性劳动就可以根据本发明的构思作出诸多修改和变化。因此，凡本技术领域中技术人员依本发明的构思在现有技术的基础上通过逻辑分析、推理或者有限的实验可以得到的技术方案，皆应在由权利要求书所确定的保护范围内。