无人机影像DSM自动生成随机传播COLVLL算法

张春森，葛英伟，郭丙轩，张月莹

1. 西安科技大学测绘科学与技术学院，陕西 西安 710054； 2. 武汉大学测绘遥感信息工程国家重点实验室，湖北 武汉 430079

数字表面模型(digital surface model，DSM)是指包含了地表建筑物、树木和桥梁等高度的地面高程模型。它在数字高程模型(digital elevation model，DEM)的基础上，进一步涵盖了除地面以外的其他高程的信息。DSM数据广泛应用于各行各业，如在森林地区可以用于检测森林的生长情况，在城区可以用于检查城市的发展情况，在军事领域，可实现精确制导和自动化监控系统指挥作战等[1]。通常通过影像密集匹配算法求得视差或者深度图，再利用视差或深度图得到DSM。

无人机影像由于其影像几何变形大、待匹配数据量多等原因，在地物高差大、弱纹理区域影像匹配会产生一系列问题[2]。如在地物遮挡处或高差大的地物上匹配到的同名点过少，导致解算影像视差值或深度值不稳定[3]。加之，现有密集匹配算法在处理诸如弱纹理、大高差等对地观测影像方面表现不佳。虽然对于密集匹配国内外学者进行了大量相关研究与改进，文献[4]提出半全局匹配(semi-global matching,SGM)算法，利用双重惩罚系数实现能量方程的最优求解，该算法通过计算8个方向或者16个方向的代价聚合实现全局能量最小化。文献[5]在SGM算法中加入中值滤波对其进行改进，然而以上两种方法在弱纹理区域均表现不佳。SGM采用代价聚合的能量函数来确定匹配代价生成视差图后融合得到DSM，但存在内存占用率高，视差图中间结果流程复杂，计算量大，弱纹理、高差大时匹配效果差。铅垂线轨迹法(vertical line locus，VLL)是一种基于物方可直接得到地面点高程的影像匹配方法，该方法能够有效地避免辐射畸变对影像匹配准确性的影响[6]，然而VLL法在地物深度不连续或高差较大的区域存在误匹配[7-9]。Colmap密集匹配是基于随机物方patch面元扫描优化深度图的方法，它基于相邻像元的深度值法向量和当前像元的深度值法向量构建多个物方patch组合优化深度图，并使用光学一致性、几何一致性进行深度图优化，通过过滤、融合得到最终深度图[10-11]。基于物方patch面元的匹配算法由文献[12]提出，它通过引入随机采样及随机传播的机制，在影像间快速传播并搜索到相似影像块进行匹配。文献[13]提出了PatchMatch双目密集匹配算法。近几年来，随着深度学习的兴起，卷积神经网络被广泛地应用于密集匹配、同步定位与三维重建[14-17]。文献[15]提出一种端到端的监督学习方式，通过构建代价体来表示左右影像的对应关系，利用3D卷积神经网络对代价体进行学习得到视差图。然而由于深度学习算法运行所需的内(显)存较大，加之需耗费大量训练和人工标注时间，其适用性和泛化性较差。

针对传统密集匹配方法在弱纹理及高差较大区域表现不佳，以及密集匹配结果融合成DSM时导致信息丢失等问题，本文提出基于随机传播COLVLL的无人机影像DSM快速生成方法，通过分析patch面元影像匹配及随机传播机制的思想，给出结合VLL算法对随机生成的高程值进行迭代更新的无人机影像DSM快速生成算法。分别选取城市高差大、弱纹理区域与一般区域无人机影像及ISPRS WGII/4提供的Vaihingen数据集进行试验验证，结果证明了本文算法的有效性。

1 本文算法

本文算法包括：①利用空三解算后的影像确定DSM测区生成范围；②针对生成范围进行扫描迭代，确定测区内关联影像，利用候选影像筛选策略对候选影像进行筛选；③对候选影像进行高程、法向量等随机初始化；④对多个物方面元组合进行计算，利用VLL策略更新当前像素高程及法向量值；⑤在每次迭代结束后，利用相关系数匹配计算得出最优组合，并进行迭代更新；⑥通过多次迭代更新最终得到DSM等步骤。算法流程如图1所示。

图1 本文算法流程Fig.1 Flowchart of the proposed algorithm

1.1 最佳候选影像筛选

本文最佳候选影像筛选包括数据获取、影像组像对打分、影像组像对筛选等步骤。

(1) 数据获取：准确的空三可提高匹配初值的精度，并且保证匹配的可靠性及降低迭代计算次数。针对空三加密后的数据确定DSM测区范围，以及测区内最大、最小高程值。在最大最小高程范围内，随机确定所有影像对应的高程值图。

(2) 影像组像对打分：以计算参考影像与每张匹配影像的交会角得分、公共可视点个数作为其打分标准。

设两条光线在物方空间中的表示分别为U(a1,b1,c1)、V(a2,b2,c2)，则交会角θ为

(1)

(3) 影像组像对筛选：基于交会角阈值、公共可视点个数、匹配影像阈值，通过式(2)对匹配影像组进行评价打分，确定最终筛选得到的最佳匹配影像组

gR(I)=∑f∈FR∩FIWa(f)Ws(f)pc(f)

(2)

式中，f为某一物方点；FR为基准影像可视物方点的集合；FI为其匹配影像可视物方点的集合；Wa(f)为在物方点f处，匹配影像I的交会角得分；Ws(f)为在物方点f处，匹配影像I的分辨率相似性得分；Pc(f)为在物方点f处，匹配影像I的公共可视点个数。

1.2 随机PatchMatch扫描迭代

PatchMatch算法最早被用于提高两个图像区域之间的近似最近邻域的计算效率，通过大量随机采样找到良好的patch匹配。如图2所示，算法流程包括核线校正、深度信息随机初始化、深度信息全局传播和深度信息随机搜索。

图2 PatchMatch算法流程Fig.2 Flowchart of PatchMatch algorithm

PatchMatch算法使用匹配代价估计的深度图作为其融合密集点云的中间结果，其深度图的生成、优化、融合过程，流程复杂，计算量大，无法直接匹配得到DSM成果。

1.3 基于VLL物方高程自适应更新

本文算法改进随机传播深度值的计算方式，通过引入VLL算法，将深度值传播改进为高程值传播。利用随机PatchMatch对影像进行扫描迭代，同时结合VLL算法对随机生成的高程值进行迭代更新得到DSM。算法的核心是扫描迭代更新高程变化，自适应计算物方patch面元。

如图3所示，物方地面三维点P(Xc,Yc,Zc)投影到参考影像I0上的像点为p0(x0,y0)，即以像点p0(x0,y0)为中心的影像窗口w0为物方patch面元的投影。过物点P(Xc,Yc,Zc)的物方patch面元由过物点P的法向量N所决定，而法向量N由法线方向角(α,β)确定。

法向量N与法线方向角(α,β)之间的关系为

(3)

图3 像方与物方patch面元投影关系Fig.3 The projection relationship between the image side and the object side patch

若物方patch面元法向量为N(a,b,c)，则patch面元方程为

a(X-Xc)+b(Y-Yc)+c(Z-Zc)=0

(4)

过投影中心S0(XS0,YS0,ZS0)，以及参考像点p0(x0,y0)的投影光线方程为

(5)

式中，(u0,v0,w0)为参考像点p0(x0,y0)在像空辅助坐标系下的坐标；λ为投影系数；(X,Y,Z)为投影光线与物方patch面元的交点坐标。像空间辅助坐标(u0,v0,w0)与参考像点p0(x0,y0)之间的关系为

(6)

式中，R为由影像外方位角元素φ、ω、κ所构成的旋转矩阵。

物方patch面元坐标系与像方坐标系间的单应矩阵H为

(7)

式中，K为影像内参数矩阵；(x0,y0)为物方patch面元中心坐标(自定义)；Rpatch是将patch方纳入像方间坐标系下的旋转矩阵；N为法向量；S为主距；C为影像外方位线元素；而(XspYspZsp)为平移向量。

基于单应矩阵根据匹配点位的变化自适应计算物方patch面元的高程及随机法向量的初值。通过迭代计算最优解，确定匹配的高程值。

归纳本文物方patch面元高程自适应计算步骤如下：

(1) 通过物方patch面元坐标系与像方坐标系间的单应矩阵H实现像点p0与物方patch面元(法向量方向角为α、β)的投影。

(2) 根据物方patch面元法向量随机值，求方向角α、β。

(4) 根据得到的匹配最优组合信息更新当前点高程信息。

(5) 根据筛选出的最优组合确定当前像素的法向量、高程等信息。

(6) 继续传播迭代，最后得到当前匹配点高程值。

2 试验分析

2.1 试验数据

为验证基于无人机影像数据采用本文算法在高差较大及弱纹理区域DSM生成效果，试验数据(图4)分别选取安徽蚌埠某地18幅无人机影像，测区内高楼林立高差较大；四川绵阳某机场18幅无人机影像，测区内建筑物影像呈弱纹理；江苏扬中某地77幅无人机影像，测区为城市一般区域。3组影像数据相关信息见表1。算法实现基于Win7 64位操作系统，开发环境为VS 2017，采用C++语言程序实现。试验硬件配置CPU为Intel Core 3.3 GHz，内存大小为16 GB。

表1 试验数据相关信息

图4 试验用无人机影像数据(部分)Fig.4 UAV image data for test (part)

2.2 大高差区域试验

针对第1组无人机影像数据，分别采用基于半全局SGM密集匹配和本文算法生成DSM，以空三结束后的特征点作为控制点，统计两组生成DSM后的平均误差及均方根误差(表2)。图5为两种算法得到的DSM结果对比，图6为两种算法得到的阴影结果对比。其中，均方根误差计算公式为

(8)

表2 匹配精度及效率比较

改变面元窗口的尺寸大小，比较相关窗口尺寸对影像匹配结果的影响(表3)。

表3 本文算法窗口尺寸对匹配结果的影响

图5 大高差建筑物区域DSM生成对比Fig.5 Comparison of DSM generation in building area with large elevation difference

图6 大高差建筑物阴影对比Fig.6 Shadow contrast of buildings with large elevation difference

2.3 弱纹理区域试验

分别以SGM密集匹配、视差图融合得到的深度图为基础生成DSM，以及本文基于随机高程值传播直接生成DSM方法对第2组弱纹理区域数据进行试验(图7)。图8为两种算法得到的DSM结果对比。图9为两种算法得到的阴影结果对比。弱纹理区域匹配精度及效率对比见表4。

图7 弱纹理影像(局部)Fig.7 Weak texture image (local)

图8 弱纹理区域DSM生成对比Fig.8 Weak texture area DSM generation comparison

图9 弱纹理区域阴影对比Fig.9 Shadow contrast in weak texture areas

表4 弱纹理区域匹配精度及效率对比

2.4 不同软件DSM生成对比试验

基于第3组数据，采用本文算法与当前国内外主流商业软件(SURE，PhotoScan，Pix4D)生成DSM进行试验对比，如图10和图11所示。

图10 不同算法DSM对比(整体效果)Fig.10 Comparison of DSM effects of different methods(overall effect)

图11 不同算法DSM对比(局部)Fig.11 Comparison of DSM effects of different methods(local)

采用真正射影像(图12)在测区范围选取一定数量特征点(10个)作为DSM定量分析检查点(表5)，其精度统计见表6。

图12 测区真正射影像Fig.12 True orthophoto (TOP) of test area

表5 定量分析检查点统计(部分)

表6 不同软件DSM生成精度统计

2.5 基于ISPRS数据集的试验

为进一步验证本文算法的有效性，使用ISPRS WGII/4提供的ISPRS Test Project on Urban Classification and 3D Building Reconstruc-tion Vaihingen数据集，对本文算法进行试验验证。该数据集包括由Leica ALS50系统获取的机载激光扫描仪数据(平均点密度约为4个/m2)、数字表面模型(DSM)数据及真正射影像TOP数据(TOP和DSM的地面采样距离均为9 cm)。研究区域包含典型的欧洲建筑类型，其结构形态各异，包括山墙、坡屋顶及其混合结构建筑。试验以Vaihingen数据集中激光扫描数据为真值，分别就本文算法生成的DSM与数据集中DSM，以及各自生成的真正射影像进行可视化对比与精度分析。

(1) DSM生成对比。图13(a)、(b)、(c)是数据集中选取的3个研究区域1、2、3的DSM，图14(a)、(b)、(c)是使用本文算法生成对应区域DSM，图15(a)、(b)、(c)对应激光扫描数据生成DSM。对比图13、图14、图15可以看出，由本文算法生成的DSM边界清晰，更接近于激光扫描数据生成DSM数据。

图13 数据集DSMFig.13 DSM of dataset

图14 本文方法生成DSMFig.14 DSM generated by the proposed algorithm

图15 LAS数据生成DSMFig.15 DSM generated by LAS

(2) 生成真正射影像对比。图16(a)、(b)是数据集中选取的两个研究区域4和5匹配结果生成真正射影像，图17(a)、(b)是采用本文算法生成对应区域的真正射影像，对比发现本文算法生成真正射影像效果较优。

图16 数据集真正射影像Fig.16 True orthophoto of dataset

图17 本文算法生成真正射影像 Fig.17 True orthophoto generated by the paper algorithm

以Vaihingen数据集中激光扫描数据为真值，分别与本文算法生成的DSM及数据集中匹配生成(自带)DSM进行对比，在DSM上随机获取采样数据点，计算不同地物处采样点的均方根误差与最大值，精度统计见表7。由表7可以看出，本文算法生成DSM的精度在山墙、坡屋顶及混合建筑物等处均优于数据集匹配生成DSM的精度。

表7 与ISPRS数据对比精度分析

2.6 试验分析

(1) 本文随机传播COLVLL在patch方采用ZNCC双边滤波加权相关系数法匹配得到较为准确的同名像点。由表3可以看出，窗口大小在一定范围内时(小于31×31像素)，窗口越大，匹配到的正确值的概率也就越高。随着物方patch窗口的变大，计算的平均相关系数随之变大，即匹配到正确值的精度在逐渐提高。试验发现，当窗口尺寸取21×21像素时，可以取得较好的匹配的精度与匹配效率，说明patch面元窗口大小选取较为重要，如果选择过小会导致高程异常值多，过大会导致匹配效率过低。当patch窗口过大时，窗口内的纹理信息增加，空间信息复杂，导致窗口相关系数求得的数值变小，匹配像点的准确度降低，故窗口尺寸的选择对匹配结果的影响较大。

(2) 由图5可以看出，本文算法能正确匹配出高差较大的建筑物，并且建筑物边缘完整。从建筑物阴影图上看(图6)，传统方法匹配的视差出现大量无效值，在进行深度图融合及高程值融合时，前者导致楼顶高程出现大量无效值。因此相较于其他方法，本文算法可以较好地处理高差较大区域。

(3) 由于本文算法能够较好地匹配得到正确高程值，由图8弱纹理(原弱纹理影像如图7所示)影像DSM生成对比可以看出生成的DSM边缘平滑，纹理清楚。同时由表4可以看出，本文算法生成DSM精度及匹配效率均高于传统算法，说明本文算法可以较好地处理弱纹理区域影像。

(4) 由表6可以看出，本文算法DSM生成精度高于商业软件，由图11可以看出，本文算法在建筑物顶面处较其他3款商业软件生成DSM更平滑，边缘更清晰。

(5) 由图13、图14、图15可以看出，本文算法生成的DSM边界清晰，更接近于激光扫描数据生成DSM数据。由表7可以看出，本文算法在山墙、坡屋顶及混合建筑物等处生成DSM精度均优于Vaihingen数据集匹配生成DSM精度。由图16、图17可以看出，采用本文算法生成真正射影像优于Vaihingen数据集中匹配结果生成真正射影像。

3 结 论

借鉴Colmap基于随机物方patch面元扫描优化深度图进而生成DSM的算法，本文提出了基于随机传播COLVLL无人机影像DSM生成算法，即通过引入VLL算法，将深度值传播改进为高程值传播。利用随机PatchMatch对影像进行扫描迭代，同时结合VLL算法对随机生成的高程值进行多次迭代更新得到DSM。由于省略了深度图等中间过程，因此克服了在弱纹理区域、高差范围过大区域可能产生的误匹配，以及流程复杂，计算量大，内存占用率高等问题，在最终生成的DSM精度及效率方面较当前国内外主流商业软件(SURE，PhotoScan，Pix4D)有了一定的改善。对于大场景超大数据量无人机影像DSM自动生成效率与精度的提升，将是笔者后续的研究工作。