一种无人机序列影像批处理三维重建方法与流程

本发明涉及一种无人机序列影像批处理三维重建方法，特别是融合低精度GPS/IMU信息的无人机序列影像批处理三维重建方法。

背景技术：

无人机能够连续获取重叠度大的高精度序列影像，但获取的影像会丢失深度信息。基于图像的三维重建，是指利用多幅数码相机图像全自动恢复出场景三维结构的方法与技术。近年来三维重建技术在视频、图像三维重建处理领域获得了巨大的成功，将其应用到无人机图像处理领域，对无人机图像进行全自动重建相关的应用，可以拓展无人机的应用范围，提高无人机的应用水平。但目前对于无人机序列影像三维重建的研究尚处于起步阶段，主要存在以下问题：(1)相对于地面影像，基于无人机序列影像的三维重建一般是大数据量大场景的三维重建；(2)大多直接将计算机视觉中成熟的算法应用于无人机序列影像三维重建中；(3)没有充分利用精度不高的辅助信息。

现在越来越多的成像系统都带有全球定位系统和惯性测量装置，可以获得包含三维世界地理坐标系(如WGS84)下地理坐标信息和相机姿态信息的序列图像。但是，这些系统赖高精度的地理定位设备，通过这些设备进行的标定以及获得的姿态和位置数据的精度一般都比图像的方式(例如，亚像素级的图像配准)要高。另一方面，目前的各种地理定位定向系统一般可以提供连续的但常常精度不高有时甚至不准确的位置和姿态信息，就像无人机所搭载的地理定位定姿系统。不幸的是，从这些设备得到的GPS/IMU数据达不到直接用于三维目标重建及导航等部分计算机视觉工作所要求的像素级的图像匹配精度要求。这些设备可以提供概略的相机姿态和位置信息，充分利用这些精度不高的辅助信息有望提高基于无人机序列影像三维重建的效率和精度。

技术实现要素：

本发明所要解决的技术问题是，克服现有技术的缺点，提供一种无人机序列影像批处理三维重建方法，本发明通过融合低精度GPS/IMU信息，利用无人机序列影像全自动地恢复出场景三维结构，得到重构点云的相对坐标。

本发明提供了一种无人机序列影像批处理三维重建方法，包括以下步骤：

步骤一、融合低精度GPS/INS信息的影像匹配：

借助于GPS/IMU系统提供的全局的位置信息，这些容易混淆的影响匹配效果的影像可以被过滤掉；假设一组影像I＝I_i，...，I_n以及对应的概略的位置姿态信息G＝G_i，...，G_n，以及一个可能匹配的视图对集的子集V_i；

㈠位置姿态信息

通过GPS和IMU获得的位置姿态信息G_i＝[R_i|t_i]，其中R_i是一个3×3的旋转矩阵，t_i是一个3维空间向量分别代表相机的位置和姿态角；标准GPS接收机得到的全局位置坐标是以经纬度和高程表示地球某一位置的WGS84坐标系；下一步，就是将GPS数据转换为地心地固坐标系，地心地固坐标系是一种可以在全球表示场景重建的笛卡尔坐标系统；相机的定向可以用偏航角、俯仰角和滚转角来表示，其中，偏航角是从地磁北起算的；这样，外部姿态G_i就包括了转换到了地心地固坐标系的GPS坐标和三个旋转角；再加上已知的相机的内方位参数，就可以得到每张影像I_i的完整的投影矩阵了：

该投影矩阵给出了相机的位置和姿态信息的概略值，这些概略值将用于后续的处理过程中；无人机所搭载的GPS/IMU系统由于系统跳变等原因有时会获取一些明显错误的位置姿态数据，所以这些概略值在用于后续的处理过程之前，必须首先进行预处理；预处理的主要内容是剔除明显错误的数据，用前后两张影像的位置姿态数据的平均值来近似代替；

需要注意的是，当错误数据出现航带的起始点或终点时，则用前后两条航带起始点或终点影像的位置姿态数据的平均值来代替；

㈡视图选择

为了识别可能存在共同对应特征点的影像，本发明为每张影像I_i选择相应的具有足够相似度的候选匹配影像集T_i＝T₁...T_k。接下来，影像集将根据影像对应的GPS/IMU信息得到的概略的重叠区域标准剔除影像；如果场景的精细三维模型是可见的，影像间的重叠区域可以很容易地通过视图I_i和I_j之间相互投影的视锥体来得到；如果场景的精细三维模型是不可显示的，只能通过估计最大景深S_i来限制影像I_i可见的区域；例如，给定了数字高程模型(DEM)，估算的地面高程可以限制相机拍摄地面的最大景深范围，地面高程可以由公开的全球30米分辨率的DEM数据得到；而且，最大景深值S_ij可以通过影像对＜I_i，I_j＞的基线恢复出来。定义S_ij：

S_ij＝t.d(G_i，G_j)，

其中，d(.，.)表示欧氏距离，t是一个决定重建所需精度的参数。给定了这些约束条件，可以通过影像对＜I_i，I_j＞的重建计算最大景深值S：

S＝min(S_ij，S_i，S_j)，

而且，这些影像必须要有重叠，为了计算一个粗略的重叠标准定义一个平行于影像I_i且到相机中心点G_i的距离为S的平面πⁱ；R_i和R_j表示视图I_i和I_j投影在平面πⁱ上的图像范围；影像的重叠度可以通过下式计算：

其中，a(.)表示投影矩形的面积；

因为SIFT等特征描述子只能适用于旋转角度小于30度的情况，这就需要视图的投影矩阵和标准平面πⁱ之间的夹角小于最大旋转角度α，否则就设置为0；对于每组影像对＜I_i，I_j＞计算其重叠区域其中I_j∈T_i；如果重叠区域值大于设定的阈值，I_j就被加入到子集V_i以用于后续的精细匹配；也就是说每一幅图像I_i只与同时满足以下两个条件的图像I_j进行匹配：

其中，表示图像I_i和I_j拍摄时的方向角，threshold为设定的阈值；

步骤二、建立极线图：

首先，对每幅图像提取尺度不变特征点；本发明的方法是采用高效的SIFT提取算子和描述算子，这种方法对于宽基线的图像匹配具有很好的适应性，尤其本发明用的是公开的SiftGPU软件；图像对之间的特征点匹配应用基于CUBLAS矩阵运算库的GPU加速的图像匹配方法；

在对每个候选的视图I_i匹配上相关的影像集V_i之后，利用五点算法进行几何验证；因为根据特征点描述子的匹配经常会出现错误的外点，采用RANSAC算法进行剔除；匹配的输出结果是用极线图表示的结构图，极线图由对应影像的顶点集v＝{I₁...I_N}和边界集ε＝{e_ij|i，j∈v}组成并是成对重建的，也就是由视图i和j之间的相对定向e_ij＝＜P_i，P_j＞和各自的影像三角化后的点集组成；其中，

P_i＝K_i[I|O]，，P_j＝K_j[R|t]；

步骤三、计算全局一致性的旋转集：

给定了极线图接下来就是确定相机的初始位置和定向信息；根据两幅图像之间的约束，两个相机的绝对位置姿态(R_i，t_i)和(R_j，t_j)需要满足旋转一致性和平移方向一致性

首先，视图对i和j之间的相对旋转集{R_ij}通过解超定方程组可以升级为全局一致性的旋转集{R_i}，

R_ijR_i＝R_j，

上述是以R_i必须为标准正交的为限制条件的；然后利用SVD(奇异值分解方法)分解使R_i满足正交约束，得到最终解可以通过解出系统初始的近似旋转矩阵(不受R_i必须为标准正交的条件限制)并用Frobenius范数将近似旋转矩阵投影到最接近的旋转矩阵的方法解决这个限制问题；

并不是所有的对极几何都是一样重要的，也就是通过外极几何得到的R_ij不是同等重要的；本文通过下面的公式来计算权重：

其中，N＝|F_ij|是视图i和j之间的内点数量，c_i，c_j是特征覆盖范围值，

其中，α是所需内点的最小数量，表示整个影像的面积，A(F_ij，r)是特征点F_ij覆盖范围经过以为圆半径进行扩张操作后的面积。除了原始的内点数量可以决定相对定向的可靠性外，覆盖范围标准更多地考虑对应特征点的空间分布；结果，具有恰当分布对应点的聚合的视图对会比那些有同样多的对应点但随机分布的视图对所占权重要高；因此，全局一致性的旋转集可以拓展为加了权重的形式，

其中，为R_i的列(k＝1，2，3)；上式可以通过稀疏的最小二乘算子，例如应用ARPACK库，解算出来；

步骤四、初始化相机中心点：

要实现相机的中心点在地心地固坐标系下的初始化，需要进行一个转换，就是旋转矩阵R_i必须转换为适应GPS的方式，这可以通过将相对定向v_ij调整为相应的GPS定向

其中，v_ij是全局坐标系中图像I_i和图像I_j之间的相对平移，是GPS坐标系中图像I_i和图像I_j之间的相对平移；这是一个典型的正交的Procrustes问题，可以通过奇异值分解法解算R；

步骤五、生成对应特征点轨迹：

极线图存储了一个相对定向集和视图对＜I_i，I_j＞之间的对应特征点；每张影像I_i都和邻近一定数量的影像匹配，匹配的信息被存储在本地的节点；需要注意的是，是一个单向图，I_i→I_j的匹配并不一定包含I_j←I_i的匹配；然后，对于极线图中的每张图节点I_i，节点被聚合成轨迹(track)其中f＝＜x^k，y^k＞表示特征点在影像I_k中的坐标位置；也就是说，根据图像匹配关系，寻找每幅图像中的每个特征点在其他匹配图像中对应的特征点，所有这些特征点，构成一个点轨迹，对应现实世界中的一个3D点；因为点轨迹是为每张影像而建，并存储在本地，起初，点轨迹集是冗余的，例如，影像I_k的一个特征点f可以属于几个不同的轨迹；然后，点轨迹应用光束法平差进行整体优化；从实用化的角度来看，由于在进行整体优化时会涉及更多的参数，冗余的估算并不可取，因此，本发明尽可能使用最小的轨迹集来表达；为了这个目的，本发明确定的一个轨迹的子集，它只包括极线图上每一个匹配的对应特征点一次；这是一个集合覆盖问题，是众所周知的最早的非完全多项式问题之一；本发明运用贪婪搜索算法来确定轨迹集的最小子集，这个最小子集随后被用来初始化稀疏的3D结构；

步骤六、初始化3D结构：

通过前面的处理步骤，可以得到相机的方向信息集合(即标定和姿态)和点轨迹还需要根据每个点轨迹确定3D点的坐标；假如相机的方向信息集总体上达不到像素级的精度且中还有外点，基于的线性三角化将导致随机的大的重建错误，即3D结构初始化错误；事实上，可以发现直接的三角化方法不能保证足够的结构初始化的精度，甚至经常连cheirality约束都不能满足；但是，极线图能提供视图对之间像素级或亚像素级精度的相机方向信息，可以满足两视图三角化的精度；因此，在每个点轨迹中基于相对定向对具有最长基线的视图对(长基线可以保证相对较低的GPS坐标误差)进行两视图三角化，也就是选取每个点轨迹中特征点所在图像的GPS坐标相差最大的两个特征点做三角化，得到初始3D点；

步骤七、光束法平差：

给定一个量测值集合，光束法平差可以通过最小化重投影误差来优化相机的定向和结构恢复；

其中，x_ij是在未知相机P_i中未知3D点X_j的观测值，v_ij是一个二进制的变量(当点X_j在影像P_i中可见时为1，否则为0)；事实上，光束法平差涉及到通过最小化非线性实值函数的平方和来调整每个3D点和相机中心点之间的射线束；光束法平差是一个庞大的，但稀疏的几何参数值估算问题，它适用于缺失数据的情况(即并不要求每一个3D点必须在每个相机中可见)；在存在高斯噪音的情况下，非线性的最小二乘法可以实现最大似然估计(必要的限制条件是初始值必须足够接近于整体最小值)；但是图像量测值高斯噪音的假设只是设定的理想状况，对于真实世界中基于自然特征的匹配技术的运动恢复结构问题是不确切的；在本发明中，特征点轨迹是基于极线约束几何调整过的，但是误匹配仍有可能出现；因此，最小平方和并不是一个合适的标准，需要一个鲁棒的目标函数来处理外点。在光束法平差的基本实现形式是最小化矢量的平方和Σ_i||∈||²，其中因此，鲁棒的目标函数可以以重设误差矢量∈′_i＝w_i∈_i的权重值来实现，如下式所示，

由此得出，上式满足Σ_iC(||∈_i||)＝Σ_i||∈_i||²的要求，其中，

权值w_i常被称为衰减因子，因为它用于降低外点的影响；

步骤八、稠密点云重建：

采用PMVS算法进行稠密点云的重建；

步骤9、纹理映射：

对重建的稠密点云基于Poisson算法构建场景或物体的网格表面，并通过遮挡分析，实现纹理自动映射。

本发明的应用基础是带有低精度GPS/IMU信息的无人机序列影像，这里低精度是指定位精度一般在10米左右，姿态角精度一般在5度左右，这样的定姿定位数据精度不能满足三维重建中图像匹配和相机定向的精度要求，无法直接依赖这些低精度的定姿定位数据进行三维重建；但是，这些低精度的GPS/IMU数据可以提供序列影像的概略地理位置和相互关系，可以提供相机的概略定向信息，这些信息都可以作为辅助信息提高基于无人机序列影像三维重建的效率和精度。

附图说明

图1为本发明的方法流程图；

图2为实例中无人机影像数据的待匹配集合图；

图3为实例中生成的极线图；

图4为实例中平差后的相机位置和姿态信息图；

图5为实例中稠密点云重建的结果图；

图6为实例中三维重建结果图。

具体实施方式

实施例1

本实施例提供一种无人机序列影像批处理三维重建方法，以下将结合附图1-6对本发明的技术方案进行详细描述。

如图1所示，首先要融合低精度GPS/INS信息对序列影像进行图像特征匹配，然后通过建立极线图、计算全局一致性的旋转集、初始化相机中心点、生成对应特征点轨迹、初始化3D结构和光束法平差等步骤完成运动恢复结构重建过程，最后通过稠密点云重建和自动纹理映射，得到了三维重建模型。

以下将通过具体的算例对本发明的技术方案的实现过程予以说明。

以下结合一组由搭载低精度的GPS/IMU组合传感器的Quickeye II型无人机获取的某校园的图像数据集，对本发明的技术方案进行详细描述。

拍摄使用的相机为Canon EOS 5D Mark II，无人机照片影像的分辨率为5616×3744。无人机飞行高度为700米，平台上挂载了定位精度在10米左右的动态单点定位的GPS，和精度在5度左右的陀螺仪，共获得153幅无人机照片影像。

本发明的技术方案通过如下步骤实现：

(1)融合低精度GPS/INS信息的影像匹配

利用低精度的GPS信息和IMU信息，获取每幅图像的待匹配的图像集合，将每幅图像的对应连接绘制成矩阵的形式，得到如图2所示的待匹配集合图。根据飞行实验前的无人机航线规划数据，图像航向重叠度最大约为70％，旁向重叠度最大约为30％，利用本文提出的视图选择方法，并将影像的重叠度阈值设为50％，可以确定最大邻近视图数d＝18，即每幅图像最多取18幅图像进行图像匹配，就可以将图像匹配时间复杂度从O(153×153)降到了O(153×18)。

(2)建立极线图

匹配的输出结果是用极线图表示的结构图，如图3所示，极线图由对应影像的顶点集v＝{I₁...I_N}和边界集ε＝{e_ij|i，j∈v}组成并是成对重建的，也就是由视图i和j之间的相对定向e_ij＝＜P_i，P_j＞和各自的影像三角化后的点集组成。

(3)计算全局一致性的旋转集

视图对i和j之间的相对旋转集{R_ij}通过解超定方程组可以升级为全局一致性的旋转集{R_i}

R_ijR_i＝R_j

上述是以R_i必须为标准正交的为限制条件的，利用SVD(奇异值分解方法)分解使R_i满足正交约束，得到最终解

根据计算，这153张影像全局一致性的旋转集为：

……

(4)初始化相机中心点

要实现相机的中心点在地心地固坐标系下的初始化，需要进行一个转换，就是旋转矩阵R_i必须转换为适应GPS的方式，这可以通过将相对定向v_ij调整为相应的GPS定向

根据计算，这153张影像的初始化相机中心点坐标为：

ImageIndex，CameraPos_x，CameraPos_y，CameraPos_z

0，-3.33200979232788，2.05530166625977，-2.15661430358887

1，-1.72906887531281，2.11266350746155，-1.87004625797272

2，-0.306412011384964，2.10093069076538，-1.73969519138336

3，-0.497344762086868，1.97058546543121，-2.0764467716217

4，1.33266019821167，1.89519202709198，-2.15731501579285

5，1.38418865203857，1.90700125694275，-2.0815737247467

6，3.32213163375854，2.05892443656921，-1.49860370159149

7，3.29912209510803，2.09036254882813，-1.41904950141907

8，3.39662861824036，2.08454847335815，-1.40254628658295

9，3.45938682556152，2.06369590759277，-1.46558439731598

10，6.04636240005493，2.01808667182922，-1.26724004745483

11，8.23372077941895，2.00806713104248，-1.07459223270416

12，11.5055103302002，1.93985688686371，-0.930897831916809

13，13.8602714538574，1.91319286823273，-0.736517250537872

……

146，16.170618057251，1.88167405128479，-0.515337228775024

147，17.8826656341553，1.84787964820862，-0.456986546516418

148，20.083517074585，1.80206346511841，-0.362621396780014

149，21.4840564727783，1.80437672138214，-0.227693825960159

150，22.4221382141113，1.85853540897369，0.0575104840099812

151，23.7377300262451，1.75472116470337，-0.0763002559542656

152，25.1186275482178，1.78453361988068，0.136511072516441

(5)生成对应特征点轨迹

根据图像匹配关系，寻找每幅图像中的每个特征点在其他匹配图像中对应的特征点，所有这些特征点，构成一个点轨迹，对应现实世界中的一个3D点，并运用贪婪搜索算法来确定轨迹集的最小子集。

根据计算，轨迹集的最小子集为：

ImageIndex，Point_id，Pos_x，Pos_y，

C000000，0，-0.885，2.106

C000000，1，-0.339，2.076

C000000，2，2.472，1.961

……

C000001，0，-0.883，0.503

C000001，1，-0.352，0.577

C000001，2，2.449，0.810

……

C000152，14876，-3.258，-2.092

C000152，14976，-2.754，-1.976

C000152，15629，-2.643，-1.972

(6)初始化3D结构

在每个点轨迹中基于相对定向对具有最长基线的视图对(长基线可以保证相对较低的GPS坐标误差)进行两视图三角化，也就是选取每个点轨迹中特征点所在图像的GPS坐标相差最大的两个特征点做三角化，得到初始3D点。

根据计算，得到7862个3D点初始化坐标为：

(7)光束法平差

对于上述步骤得到解算结果，光束法平差可以通过最小化重投影误差来优化相机的定向和结构恢复，平差后的相机位置和姿态解算结果如图4所示，具体数值如下：

PhotoID，X，Y，Z，Omega，Phi，Kappa

0，0.0254479340093484，0.3415414621610514，0.0250366517649538，175.7103213867737800，2.0291928596802138，2.0673573899635556

1，0.0410902629089479，0.1477144927140597，0.0294442811860536，178.3745342696717300，4.0489057627394587，-3.1093510454025033

2，0.0684955405233049，-0.0479907859590106，0.0364048464871272，-176.2261276627340600，2.9858672357529628，-12.5320937664505990

……

152，2.1318594371907205，0.2284795097077054，0.3511295664437739，-179.6044286840305200，13.2343419289542870，158.7878915203033300

(8)稠密点云重建

采用PMVS算法进行稠密点云的重建，重建结果如图5所示。

根据计算，共得到1640119个稠密点云，坐标为：

(9)纹理映射

对重建的稠密点云基于Poisson算法构建场景(或物体)的网格表面，并通过遮挡分析，实现纹理自动映射，结果如图6所示。

以上实施例仅为说明本发明的技术思想，不能以此限定本发明的保护范围，凡是按照本发明提出的技术思想，在技术方案基础上所做的任何改动，均落入本发明保护范围之内。