基于无人机低空高分辨率影像的密集点云提取方法

摘要

本发明涉及一种基于无人机低空高分辨率影像的密集点云提取方法，包括以下步骤：将N-1张配对影像分别与基准影像进行配对，组成立体像对；以一个立体像对为单位，遍历基准图中的所有像素点，采用不等间隔深度采样算法，搜索相应配对影像中的同名核线空间；采用Daisy算法分别计算同名核线空间中每个点于基准图中对应像素点的Daisy描述量，获取同名点集；提取每个立体像对应的密集三维点云；遍历配对图像，提取出N-1套三维点云集；利用强制连续性检核算法获得最终的密集三维点云。与现有技术相比，本发明搜索速度快，具有更强的可靠性及执行效率，非常适用于密集点云之间的相关性量测，能有效提取基于无人机低空致密三维点云。

说明书

技术领域

本发明涉及一种点云提取方法，尤其是涉及一种基于无人机低空高分辨率影 像的密集点云提取方法。

背景技术

基于航空或卫星影像重建地表区域的数字高程模型(DEM)和数字表面模型 (DSM)一直都是航空摄影测量的主要任务和研究重点。随无人机技术的不断发 展及国内低空管制政策的逐渐放宽，将无人机作为传感器的搭载平台，于千米以下 低空飞行获取航摄数据的作业模式已被越来越多的国内企业及科研单位采用。采用 这种作业模式既可以充分发挥无人机的高机动性，又能弥补大飞机及卫星因飞行高 度过高而造成的地物细节信息缺失的缺点。

但由于目前国内并没有针对低空无人机高分辨率影像的密集点云提取方法， 现有的其它点云提取方法应用于此问题时效能不佳，均不适用。故目前一般的解决 方案是采用无人机搭载激光扫描仪在低空对地扫描获取密集点云，如中国专利 CN103426165A公开的一种地面激光点云与无人机影像重建点云的精配准方法，现 有技术根据激光点云数据提取点云道路矢量线，经过处理与遥感影像配准。但激光 扫描仪与单反相机相比成本过高，不适合推广。

发明内容

本发明的目的就是为了克服上述现有技术存在的缺陷而提供一种低成本、高 效率、高可靠的基于无人机低空高分辨率影像的密集点云提取方法。

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

一种基于无人机低空高分辨率影像的密集点云提取方法，包括以下步骤：

1)在N张无人机影像序列中，挑选一张影像作为基准图，其余作为配对影像；

2)将N-1张配对影像分别与基准影像进行配对，组成立体像对；

3)以一个立体像对为单位，遍历基准图中的所有像素点，对每个像素点采用 不等间隔深度采样算法，搜索相应配对影像中的同名核线空间；

4)采用Daisy算法分别计算同名核线空间中每个点与基准图中对应像素点的 Daisy描述量，并根据公式(1)计算两点间的相关性概率分布，获取同名点集：

$P\left(d\right)=\frac{1}{Z}\exp (-\frac{{\left|\right|D_{\mathrm{base}}-D_{\mathrm{match}}\left|\right|}^2}{\sigma })---\left(1\right)$

其中，Z为归一化系数，σ表示相关性概率分布的方差，D_base为基准图中像素 点的Daisy描述量，D_match为配对影像同名核线中点的Daisy描述量；

5)利用影像对应的内、外方位元素及同名点集，提取每个立体像对对应的密 集三维点云；

6)遍历配对图像，重复3)-5)，提取出N-1套三维点云集；

7)利用强制连续性检核算法剔除离群点，获得最终的密集三维点云。

所述步骤3)中，不等间隔深度采样算法为：根据基准图中目标点的像素坐标、 立体像对中各影像对应的内外方位元素、立体像对对应地面区域的深度范围和采样 间隔，获取相应配对影像中与所述目标点相对应的同名核线空间。

所述Daisy算法中，关键点位方向只沿着垂直于同名核线的方向。

所述Daisy算法中，参数设置如下：

描述子的作用半径R＝8，子区层数Q＝2，每层子区数T＝4，梯度方向数H＝4， 参与运算的子区数S＝9，描述子的最终维数Ds＝36。

所述步骤7)中，强制连续性检核算法具体步骤如下；

701)获取N-1套三维点云集，设定限差阈值error_threshold和连续性阈值 c_threshold；

702)选取其中一套数据作为基准点云数据集B；

703)选取除基准点云数据之外的另一套数据作为检核数据集C_i；

704)遍历基准点云数据集与检核数据集中的所有点云，根据公式(2)分别 计算基准点云数据集中每个点的偏差值e，当e<error_threshold时，将该点对应的 连续性计数加一；

$e=\sqrt{{(x_B-x_{C_i})}^2+{(y_B-y_{C_i})}^2+{(z_B-z_{C_i})}^2}---\left(2\right)$

705)根据步骤703)、704)遍历所有三维点云数据；

706)保留所有连续性计数小于c_threshold的三维点，将这些保留下来的点作 为最终三维点云输出。

与现有技术相比，本发明具有以下优点：

1、本发明使用不等间隔深度空间采样算法提高了核线搜索算法的速度。

2、本发明使用改进后的Daisy算法计算像素点间的相关性，具有更强的可靠 性及执行效率，非常适用于密集点云之间的相关性量测。

3、本发明使用单反相机就可以实现密集点云的提取，与无人机搭载激光扫描 仪在低空对地扫描获取密集点云的方法相比，成本大大降低。

附图说明

图1为本发明的流程示意图；

图2为SIFT描述子原理中已旋转至主方向的局部像素矩阵示意图；

图2中，圆圈的覆盖范围表示一个圆形的高斯滤波掩膜；每个小方块表示一 个像素区域；图中的箭头方向及长度表示网格中对应像素点的梯度方向和梯度大 小；

图3为SIFT描述子原理中各子区在八个不同方向的加权梯度累计直方图；

图3中，每个网格覆盖了4*4大小的区域，图中的箭头指向表示对应的梯度 方向，箭头长度表示该子区中所有点位在该方向上的梯度加权和；

图4为Daisy描述子的形状示意图；

图4中，圆圈表示关键点邻域内的子区范围；+表示当前子区的中心点位置； 子区的半径大小等于作用于该区域中高斯卷积核的标准差大小，以子区半径作为依 据可以将关键点周围的子区分成若干层，位于同一层内的子区用相同颜色表示，反 之则用不同的颜色表示；位于相同层内的子区中心点到关键点的距离相同，且子区 半径相同；离关键点越远的子区半径越大也即其对应的高斯卷积核标准差越大；图 中direction-j所标注的点位即是关键点所在位置。

具体实施方式

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

如图1所示，本实施例提供一种基于无人机低空高分辨率影像的密集点云提 取方法，包括以下步骤：

步骤S01，获取无人机影像序列内、外方位元素；

步骤S02，在N张无人机影像序列中，挑选一张影像作为基准图，其余作为 配对影像；

步骤S03，将N-1张配对影像分别与基准影像进行配对，组成立体像对；

步骤S04，选取一个立体像对，遍历基准图中的所有像素点，对每个像素点采 用不等间隔深度采样算法，搜索相应配对影像中的同名核线空间；

步骤S05，采用Daisy算法分别计算同名核线空间中每个点于基准图中对应像 素点的Daisy描述量，并根据公式(1)计算两点间的相关性概率分布，获取同名 点集：

$P\left(d\right)=\frac{1}{Z}\exp (-\frac{{\left|\right|D_x^i-D_{x^{\prime }\left(d\right)}\left|\right|}^2}{\sigma })---\left(1\right)$

其中，Z为归一化系数，σ表示相关性概率分布的方差，D_base为基准图中像素 点的Daisy描述量，D_match为配对影像同名核线中点的Daisy描述量；

步骤S06，利用影像对应的内、外方位元素及同名点集，提取每个立体像对对 应的密集三维点云；

步骤S07，所有立体像对是否遍历完毕，若是，则执行步骤S09，若否，则执 行步骤S08；

步骤S08，获取下一立体像对，返回步骤S04；

步骤S09，提取出N-1套三维点云集，利用强制连续性检核算法剔除离群点；

步骤S10，获得最终的密集三维点云。

上述方法的关键步骤具体描述如下：

1、不等间隔深度采样算法

本方法采用不等间隔深度空间采样算法搜索目标点对应的同名核线空间，与 传统的核线搜索算法相比速度有较大提升。不等间隔深度空间采样算法简述如下：

功能：在一个立体像对中，已知左片某点的像素坐标X及两张影像的内、外 方位元素及目标地物区域的深度范围，在右片中生成与左片X点对应的同名核线 空间。

概述：

设一张影像对应相机的内方位元素为K，外方位元素为R，该相机在世界坐标 系下的坐标为C。那么就可建立世界坐标系下一点与该点在影像上像点之间的对应 关系，如下所示：

λx＝KR(X-C)   (1)

其中，表示物方点X在影像上的相片坐标；λ表示该点相对于相机平面的深 度值。

利用公式(1)可以反求出在任意深度λ上与像点x对应的物方点三维坐标， 如公式(2)所示。

X(λ)＝λR^TK^-1x+C   (2)

假设现在用两个架设在不同位置上的相机对物点X成像。两个相机的参数分 别为(K₀,R₀,C₀)和(K₁,R₁,C₁)，其中K表示相机内方位元素，R表示相机的旋转矩 阵，C表示相机在世界坐标系下的坐标。物点X在左片上的像点是x_l，为了在右片 中获得一定范围内与xl对应的等间隔同名核线点，需要左片摄影光束上进行不等 间隔采样，假设某次采样的左片深度间隔为dλ，核线间隔为r，右片深度间隔为 dω。则结合公式(2)可得公式(3)：

$(\omega +\mathrm{d\omega })\left(\begin{array}{c}u+\mathrm{rdu}\\ v+\mathrm{rdv}\\ 1\end{array}\right)=(\lambda +\mathrm{d\lambda })a+b---\left(3\right)$

其中ω表示物方采样点距右相机的深度范围；[u v 1]^T表示物方采样点在右片上 的像点坐标，$a={\left(\begin{array}{ccc}{a}_0& a_1& a_2\end{array}\right)}^T=K_1{R}_1{R}_0^T{K}_0^{-1}x,b={\left(\begin{array}{ccc}{b}_0& b_1& b_2\end{array}\right)}^T=K_1{R}_1(C_0-C_1).$

设物方点的深度范围为[λ_min,λ_max]，其对应到右片核线方向的起始、终止坐标 为(u_min,v_min)、(u_max,v_max)。核线采样间隔对应的坐标轴间隔为du,dv。则可获得如 下关系(4)：

$\mathrm{dl}=\sqrt{{(u_{\max }-u_{\min })}^2+{(v_{\max }-v_{\min })}^2}$

du＝(u_max-u_min)/dl            (4)

dv＝(v_max-v_min)/dl

结合(3)、(4)式可得(5)：

$\mathrm{ad\lambda }=\omega \left(\begin{array}{c}\mathrm{rdu}\\ \mathrm{rdv}\\ 0\end{array}\right)+\mathrm{d\omega }\left(\begin{array}{c}u+\mathrm{rdu}\\ v+\mathrm{rdv}\\ 1\end{array}\right)---\left(5\right)$

整理(5)式可得(6)：

$\mathrm{d\lambda }=\frac{\mathrm{\omega rdu}}{a_0-a_2(u+\mathrm{rdu})}=\frac{\mathrm{\omega rdv}}{a_1-a_2(v+\mathrm{rdv})}---\left(6\right)$

根据公式(6)就能对右片上的同名核线空间进行不等间隔采样。

输入：

①立体像对中两种影像分别对应的内、外方位元素(K₀,R₀,C₀)和(K₁,R₁,C₁)；

②左片中目标点的像素坐标L(x,y)；

③立体像对对应地面区域的深度范围(λ_min,λ_max)；

④采样间隔r；

输出：右片中与左片对应的同名核线空间。

不等间隔深度空间采样算法的伪代码如表1所示：

表1

2、改进的Daisy算法

本方法中采用了改进后的Daisy算法计算像素点间的相关性。与其他匹配测度 相比，具有更强的可靠性及执行效率，非常适用于密集点云之间的相关性量测。改 进后的Daisy算子简介如下：

●Daisy算法概述

Daisy描述子是由SIFT和GLOH描述子演变而来的，而GLOH又和SIFT关 系紧密，所以归根结底Daisy是经由SIFT改进后得到的。Daisy与SIFT相比最大 的优势是计算速度快，这主要得益于Daisy对影像局部区域描述策略的改进以及省 略了尺度空间构造和关键点主方向评估这两个关键步骤。而正是由于这两个步骤， 才使得SIFT描述子具有旋转，尺度以及对光照的不变性。所以，Daisy算子在理 论上既不具备尺度不变性也不具备旋转不变性。由于Daisy是经由Sift改进而来， 这里首先介绍Sift描述子的实现原理，进而推导出Daisy描述子的相应改进策略。

●SIFT描述子原理概述

sift描述子的最大特点是具备尺度和旋转的不变性。原作者David G.Lowe通 过构造影像尺度空间实现尺度不变性；通过计算关键点一定领域范围内的梯向直方 图实现旋转不变性；通过金字塔结构中的影像间差分(DOG)及最终描述向量的 归一化实现对光照变化的鲁棒性。因此只有在特定尺度影像中的特定方向下，计算 出某关键点的128维向量才能真正对该点实现上文中所提及的两个不变性描述。而 由于Daisy算子本身并不包含影像尺度和点位主方向的计算，故本节主要介绍sift 描述子128维向量的生成过程。

以图像中某个关键点为例。在该点位对应尺度的影像上定位到该点位置，同 时以当前点为中心，该点主方向为正北方向开辟一个局部影像窗口，此窗口中的所 有像素值都需要经原始影像双线性内插后获得，如图1所示。图1中网格的中心点 即为待描述的关键点。

将关键点邻域分成2*2的4个子区域，每个子区域中包含4*4个像素点。分 别计算每个像素点在八个不同方向上的梯度值，这样每个子区域就拥有八张对应的 梯度图。利用标准差为当前影像模糊度1.5倍的高斯滤波器平滑每个子区中对应的 八张梯度图(等效于对子区中不同像素点的梯度值进行加权，每个点位的权重与该 点位置离子区中心点位置的距离成反比)，获得八张卷积后梯度图。将每张卷积后 梯度图中的所有带权梯度值相加，就获得了每个子区在八个不同方向的梯度描述 量，如图2所示。在图1中，一共划分了4个子区，每个子区拥有8个方向的梯度 描述量。即在此例中，一个关键点有2*2*8＝32个描述量。实际上，sift算子是将 关键点的周围邻域划分成4*4的子区，也即最终会获得4*4*8＝128维的描述量。 为了使描述子对光照变化具有一定的鲁棒性，还需将128维向量进行归一化处理， 最终生产的归一化描述向量也即使sift闻名于世的128维特征向量。

●Daisy描述子原理概述

Daisy描述子与SIFT相比，同样是将关键点周围的领域分成若干子区，并将 每个子区N个不同方向上的加权梯度和作为描述量，不同之处在于子区的划分方 法、描述量的计算方式以及最终特征向量的维数(Daisy描述子的最终维数完全依 赖于初始参数的设定)。图3展示了Daisy描述子的大致形状，表2展示了Daisy 描述子的各项参数及含义。

表2 Daisy算子参数一览

●Daisy算法流程

将当前待描述的关键点记为P，以P点为中心，R个像素为半径开辟一个圆形 窗口。以为角度间隔，计算窗口内所有像素在H个方向上的梯度值，并将梯 度值赋予当前像素值，生成在该方向上的梯度图G₀，其计算公式为：

$G_0={\left(\frac{\partial I}{\partial H}\right)}^{+}---\left(2\right)$

其中，I表示原始影像，则表示在H方向求取影像中各点的梯度值，计算 式(a)⁺＝max(a,0)，这样每个窗口都会获得H张对应的梯度图。

子区位置的确定实际上就是子区中心点位置的确定。公式(3)以关键点位置 为中心，采用极坐标系计算不同子区中心点的位置坐标。

$\left(\begin{array}{c}{r}_i=\frac{R(i+1)}{Q}\\ {\theta }_j=\frac{2\mathrm{\pi j}}{T}\end{array}\right)---\left(3\right)$

其中，i表示子区所在的层数；r_i表示当前子区中心到关键点之间的距离；θ_j表 示该点与关键点之间连线与水平方向的夹角，0≤j＜T，j∈N；R、Q、T为Daisy 算子参数，详情请参见表2。

每个子区的半径大小等同于作用于该子区上的高斯滤波器的标准差大小，其值 由式(4)计算得到。

${\sigma }_i=\frac{R(i+1)}{2Q}---\left(4\right)$

其中，i表示当前子区所在的层数；R、Q皆为Daisy算子参数，σ_i表示当前 子区的半径。

采用标准差与子区半径大小相等的高斯滤波器对每个子区进行平滑，也即对子 区内所有梯度加权。将同一个子区内加权后的梯度相加获得一个统计量M。由于 一共有H个方向的梯度图像，每个梯度图像中有S个子区，故最终形成的特征向 量共有H*S维。

原始Daisy算法中采用的归一化方式是将每个子区所对应8个方向的统计量作 为一组进行归一化，然后再将各组归一化后的描述量组成特征向量。这样做的主要 目的是为了尽量保证影像中遮蔽区域附近点位的正确描述。而在实际应用中应当根 据具体需求来调整归一化策略。

●Daisy算法在本方法应用中的改进

Daisy算子在本方法中的改进主要体现在以下四个方面：

①尺度影像的选择

因为Daisy算子并不是基于尺度空间来对点位进行描述的。故Daisy算子本身 不具有尺度不变性。为了尽量减少影像尺度变化对最终描述量的影响，需要事先评 估影像间的尺度差异，将尺度差异较小的影像对作为目标计算Daisy描述量。

②关键点位方向的确定

由前文可知，Daisy描述子的计算过程中并不涉及点位主方向的评估，故纯净 的Daisy算子也不具备旋转不变性。为了使Daisy描述子具备该性能，需要获得关 键点的主方向，并将局部窗口旋转至主方向后再计算特征向量。在本方法中只需沿 着垂直于核线方向计算特征向量即可。

③描述子归一化策略的调整

在本方法中，为了增加描述子的稳定性，描述子将进行全局归一化，而不是像 原有算法中采用以子区为单位局部归一化策略。

④参数的调整

在本方法中，各参数设置为R＝8，Q＝2，T＝4，H＝4，S＝9，Ds＝36，也即最终 的描述子为36维。

3、强制连续性检核算法

本方法利用强制连续性检核算法剔除出离群点，获得最终的密集三维点云。 强制连续性检核算法具体为：

功能：利用同一地物区域对应的多套三维点云数据，将其中一组数据作为基 准，剔除离群点。

输入：同一地区对应的多组三维点云数据。

输出：剔除粗差后的密集三维点云数据。

步骤：

1)输入点云数据，设定限差阈值error_threshold,设定连续性阈值c_threshold。

2)选取其中一套数据作为基准点云数据B。

3)选取除基准点云数据之外的另一套数据作为检核数据_Ci。

4)遍历基准数据集与检核数据集中的所有点云，根据公式(5)分别计算每 个点的偏差值e，当e＜error_threshold时，就将该点对应的连续性计数加一。

$e=\sqrt{{(x_B-x_{C_i})}^2+{(y_B-y_{C_i})}^2+{(z_B-z_{C_i})}^2}---\left(5\right)$

5)遍历所有三维点云数据，按3)-4)步重复处理。

6)最后将所有连续性计数小于c_threshold的三维点保留，将这些保留下来的 点作为最终三维点云输出。