无人机大场景在线密集点云和DSM生成算法

杨佳琪, 范大昭, 杨佳宾,, 杨幸彬, 纪 松

(1. 信息工程大学地理空间信息学院,河南 郑州 450001; 2. 北京字跳网络技术有限公司,北京 100086)

无人机遥感是近年来兴起的一种新型遥感系统,具有快速、灵活、低成本、高分辨率等特点[1]。传统雷达(light detection and ranging, LiDAR)点云由于其成本较高,难以得到广泛应用。随着航摄影像质量的不断提高,利用无人机影像处理软件[2]获取三维信息已成为一种新型测绘手段。然而由于无人机影像的数据量大、软件处理算法复杂,数据处理常依赖于离线工作流[3],难以在线获取所需三维信息,因此在自然灾害的快速应急响应中难以发挥良好的作用。

近年来,利用同步定位与地图构建(simultaneous localization and mapping,SLAM)技术在线获取无人机影像位姿普及[4]。相较于LiDAR点云,影像密集匹配点云具有丰富的纹理特征信息,更易于高效处理[5],且成本相对较低。为了实现高质量的影像密集匹配,文献[6]提出一种改进代价的半全局匹配算法,以获取视差图,通过限制视差范围减少了计算资源。文献[7]采用双目视觉方法实时获取三维点云,实现了对远距离目标三维信息的实时获取。文献[8]提出一种基于纹理优化的动态约束视差半全局匹配(semi-global matching,SGM)算法,通过对半全局立体匹配代价积聚过程的优化,大大减少了深度估计过程的耗时。文献[9]提出一种基于多视影像间的重叠关系与方向确定高分辨率航空多视影像密集匹配算法,完成了密集点云的生成与融合。物体边缘往往具有较大的视差变化,文献[10]提出了一种顾及影像纹理信息的视差图膨胀腐蚀算法,结合金字塔影像进一步约束视差搜索范围,减少了算法所需的内存开销和计算时间。上述几种方法尽管在一定程度上优化了算法耗时和效果,然而仅能考虑两视图信息,在应急响应中难以满足需求。文献[11]针对无人机短时间应对自然灾害设计了一个在线三维重建框架,使用开源的块匹配算法逐步增量生成密集的地理参考三维点云和数字表面模型(digital surface model, DSM),但在影像分辨率过高时无法实时运行。文献[12]设计了一种实时测图框架,基于视差平面的快速立体匹配算法实现了影像密集匹配过程,然而获取的深度和点云信息精度有限。为了最大限度减少计算时间、提高计算精度,文献[13]提出了一种基于物方的半全局匹配方法直接获取DSM,但在计算效率和精度方面优化较少。文献[14]提出了一种基于半全局铅垂线轨迹匹配的直接DSM生成方法,可以考虑多视信息在物方空间中直接生成DSM,但无法实现在线计算。

为了解决上述问题,本文提出一种无人机在线密集点云和DSM生成算法,结合增量式的多视图半全局在线密集匹配算法获取实时点云与DSM。首先,利用SLAM和实时差分定位(real time kinematic,RTK)松耦合方式在线估计影像位姿,在此基础上选择当前帧的多视图参考帧。然后,在深度空间计算离散采样深度,将采样深度投影到多视图参考帧上计算联合匹配代价,并考虑多视图遮挡关系对代价进行优化。最后,基于增量式半全局积聚策略得到积聚代价体,利用“胜者为王”算法结合抛物线拟合得到最优深度图,利用深度图融合得到密集点云和DSM。利用3组典型地区的无人机影像对本文方法进行在线验证,基于定性和定量评价验证本文方法生成密集点云和DSM的有效性。

1 无人机在线密集点云和DSM生成流程

基于无人机影像在线生成密集点云和DSM的流程主要包括:影像在线位姿估计、影像在线密集匹配获取深度图、深度图融合和DSM生成4部分。为了在线估计影像位姿,在ORB SLAM3[15]系统框架基础上进行了改进。无人机系统上搭载的RTK可实时获取无人机位置和轨迹,通过对视觉SLAM和RTK获取的位置信息进行松耦合绑定,从而在线获取地理坐标框架下的无人机影像位姿信息。

1.1 影像在线密集匹配

针对仅利用两视图信息可能会因视角变化存在遮挡变形和直接将场景深度作为深度空间搜索范围进行求解时存在较大困难的问题,本文提出一种增量式的在线影像密集匹配方法,利用多视图信息直接在深度空间计算联合代价。另外,为节省在线深度估计计算时间,增量式地将前一帧深度映射至待求解帧,增量式约束代价体计算和积聚过程。具体算法如下:

(2)多视图联合代价体计算。为了计算当前帧深度图,直接在深度空间进行解算,基本思路为利用图像位姿将不同深度投影至多帧图像计算匹配代价。首先根据无人机飞行时相对航高H和相机大致倾角φ,利用z=Harcsinφ确定深度空间搜索范围[zmin,zmax];然后利用式(1)将最小和最大深度投影至多个参考帧图像,计算两个投影点之间像素距离的最大值。

(1)

式中,zl为给定深度值;p为当前图像像素点齐次坐标,p′l为像点p在参考帧图像的投影点坐标;K、R、T分别为相机内参和相机到世界坐标系的旋转和平移矩阵。为在深度空间求解,本文在深度空间内进行非等间隔采样,离散采样深度值zl(l=0, 1,…,L)计算公式如下

(2)

按照式(2)将像点p按照不同zl投影至参考帧图像得到极线上像点坐标p′l,后续过程中在D轴采样方向上将每个采样深度位置对应一个匹配代价值,计算过程如图1所示。

图1 多视图联合代价体计算过程

为了计算像点p在不同深度下的匹配代价,本文对文献[16]的census代价进行改进,采用中心对称思想避免中心像素点噪声对代价的影响,同时节省了一半的计算耗时。考虑视图间像素遮挡概率,联合多视图信息计算代价权重,通过权重联合匹配代价计算代价矩阵C。中心对称census比特链的计算函数式为

P(u+i,v+j)]

(3)

(4)

式中,WCTfj(u′,v′)为P(u,v)所对应投影点的中心像素灰度P′(u′,v′)的比特值。

(3) 增量式代价积聚与深度计算(考虑深度图映射)。考虑到Census局部代价本身的约束,为了获得更高的深度图精度与完整度,结合半全局匹配方法中扫描线优化策略,对代价矩阵C进行全局能量的聚集,使代价矩阵更加精确可靠。鉴于本文方法在深度空间上直接计算多视图联合代价,因此将传统全局能量函数E(l)[17]修改为

(5)

式中,lp为像点p的深度采样位置;lq为像点p临近像点q的深度采样位置;c(p,lp)为像点p在lp处的多视图联合匹配代价;T()为判断函数,其目的是对像点p临近视差差值等于1和大于1的像点q代价分别增加P1、P2惩罚项。

(6)

在代价积聚过程中,对fi帧图像上有先验深度的区域按照一定的容差范围对采样位置进行代价积聚;对于无先验深度的区域,按照正常积聚策略完成代价积聚,具体计算过程为

(7)

以上计算仅考虑了一个方向r的代价聚合值Lr(p,lp),聚合后的信息量较少,因此考虑从8个方向对代价矩阵进行代价积聚(如图2所示),最终生成积聚后的代价矩阵。

(8)

图3显示了利用两视图半全局立体匹配方法[17](简称为两视图法)和本文方法生成的深度图和点云。可以看出,本文方法相较于两视图法更好地保护了深度图的边缘与棱角信息(如图3(a)所示);此外,在弱纹理处(如图3(b)所示)能够较大程度上改善代价计算的准确度,从而使深度图更加平滑精确。

1.2 增量式点云与DSM生成

为了在线获取测区密集点云,首先将前一帧fi-1深度投影至当前帧fi得到先验深度;然后将fi估计的深度与先验深度作比较,剔除相对误差大于5%阈值范围的噪声,对剔除后的深度进行等权融合得到新的深度图;最后按照式(8)将新的深度图增量式地投影至三维空间,即可得到稠密点云结果。

(9)

2 试验与分析

试验选取乡村、山地、城区3组典型无人机影像测试本文算法生成稠密点云和DSM的效果。乡村区域包含弱纹理建筑物和边缘轮廓清晰的伐木场,用于评定本文方法在弱纹理和物体边缘的匹配效果;山地区域地形起伏较大,不同视角拍摄的图像存在一定变形和遮挡区域,用于评价本文方法在高差大的山地区域点云生成的稳健性;城区数据用于评定本文方法在城区的效果。3个测区拍摄的影像位置(SLAM在线位姿估计得到)和本文方法生成的稠密点云如图4所示。

图4 3个测区拍摄时相机位置和稠密点云结果

本文算法在Windows 10 64位系统下基于VC++2015开发,其中影像位姿估计过程结合了ORB-SLAM3在线获取,影像密集匹配、增量式点云融合与DSM生成过程为本文核心算法,密集匹配和点云融合过程使用OpenCL进行了并行加速。试验使用的硬件平台为华硕笔记本,处理器为Intel® CoreTMi7-8700 CPU、最高主频为3.2 GHz,CPU内存为15.6 GB,图形显卡为NVIDIA GTX1060。

2.1 点云融合和DSM生成效果

为验证本文方法生成的稠密点云效果,在3组数据的基础上选取局部典型地物稠密点云进行放大显示,如图5—图7所示。由图5可知,无论是伐木场还是弱纹理房屋屋顶的三维点云都具有较好的精度和完整度;由图6可知,即使是在地形突变的区域(图6中②③)点云形状和走势也符合实际地表形态,说明本文方法能够在一定程度上抗拒视角和遮挡影响,稳健性较好;由图7可知,点云完整度较高,由于影像是从视频中抽帧得到,因此点云精度相比于乡村和山地略低,但从整体上看本文方法生成的点云在建筑物和树木上具有较好的完整度和精度,符合实际地物形态。

图5 乡村区域稠密点云结果

图6 山地区域稠密点云结果

图7 城区稠密点云结果

图8为本文增量式DSM生成的效果图,可以看出无论是在弱纹理房屋的乡村区域还是地面起伏较大的山地区域,本文方法生成的DSM均能保持完整地表形态,且平滑性较好;从建筑物的边缘细节和地形断面处可以看出,本文方法生成的DSM具有清晰的轮廓。

图8 3个测区DSM

2.2 点云精度定量评价与分析

为进一步定量评价本文方法生成的稠密点云绝对精度,导入3组数据的位姿数据到SURE软件,离线生成稠密点云且与本文生成结果进行叠加,通过绝对值误差评定本文方法精度。图9为3组数据的误差分布。

图9 3组数据误差分布

由图9可以看出,本文方法生成的稠密点云在大部分区域内误差较小,乡村地区的误差最小,相比于SURE软件生成的稠密点云,误差小于3 cm的像素数占比超过90%。

对图9误差进行统计分析,剔除不符合正态分布3σ的误差后,计算一定范围内最大绝对误差、绝对值均值误差和均方根误差来定量评定本文方法生成的点云精度,结果见表1。由均方根误差可知,本文方法生成的点云在乡村地区较小,精度最高,在山区和密集建筑物的城区误差相对较大;由最大绝对误差和绝对值均值误差可知,受视频影像质量影响城区最大绝对误差和绝对值均值误差最大,山地次之。

表1 3个测区的点云精度评定结果 m

2.3 计算效率对比分析

为定量评价本文算法的性能,将本文方法计算的位姿结果导入商业软件SURE和开源软件COLMAP后离线生成稠密点云和DSM,然后分别统计两种方法与本文方法的计算耗时,结果见表2,需要说明的是,3种方法均在相同的硬件环境下测试,使用COLMAP[19]时关闭了几何一致性深度计算。

表2 算法的计算性能对比结果

由表2统计结果可知,本文方法平均一帧计算耗时仅约为0.4 s,具有较高的计算效率。由于本文方法和SURE软件(密集匹配采用类似半全局立体匹配的思想)的输入输出数据相同,因此具有较好的可比性。由3组数据的对比结果可知,本文方法在精度没有损失太多的情况下计算效率要远高于SURE软件(计算效率是SURE软件的约25倍),平均计算效率基本可以达到2～3帧/s,足以满足无人机在线生成稠密点云和DSM的要求。

3 结 论

传统无人机影像获取密集点云和DSM的方法耗时过多、点云质量和处理速度难以有较好的平衡,无法实现对自然灾害的应急响应。本文针对此问题提出了一种无人机大场景在线密集点云和DSM生成方法。通过直接在深度空间内引入多视图中心对称census代价,并考虑视图遮挡关系计算联合代价体,减少一半计算耗时的同时提升了代价准确性,在此基础上提出增量式代价积聚策略,联系当前帧与历史帧深度进一步优化耗时且提高深度图精度,将相邻帧深度进行检查并剔除噪声后,快速在线生成密集点云与DSM。选择典型的乡村、山地、城区3组测试数据进行试验。结果表明,本文方法能够以2～3帧/s的平均处理速度在线生成测区稠密点云和DSM信息,同时生成的点云精度与SURE软件相比,可控制在为5 cm误差范围内,DSM上地物边缘清晰、平滑度较好,可以为自然灾害的快速应急响应提供精确实时的三维信息。