CN103426200B

CN103426200B - 基于无人机航拍序列图像的树木三维重建方法

Info

Publication number: CN103426200B
Application number: CN201310375676.4A
Authority: CN
Inventors: 冯志勇; 陆泽萍; 徐超; 张东萍; 王丹丹
Original assignee: Tianjin University
Current assignee: Tianjin University
Priority date: 2013-08-26
Filing date: 2013-08-26
Publication date: 2017-05-03
Anticipated expiration: 2033-08-26
Also published as: CN103426200A

Abstract

本发明公开了一种基于无人机航拍序列图像的树木三维重建方法，首先通过分水岭分割方法提取树冠的特征区域，然后通过计算RGB颜色空间中的区域相关系数来匹配特征区域并提取匹配特征点对，再结合航拍摄像机的成像模型，通过双目立体视觉原理计算深度信息，最后通过L系统规则构建树木的三维模型。与现有技术相比，本发明考虑了树冠的纹理、颜色以及轮廓信息，充分提取出能反映树冠复杂结构的特征点集合。并且，考虑到无人机航拍过程中，摄像机存在的抖动、平移等，针对航拍摄像机复杂的成像模型来计算特征点集合的深度信息，并采取了合理的树木建模方法，构建出树木的合理的近似三维模型。

Description

基于无人机航拍序列图像的树木三维重建方法

技术领域

本专利涉及图像三维重建技术，特别是涉及基于无人机航拍序列图像的三维重建方法。

背景技术

无人机航拍是近年来空间数据获取的重要数据.来源。无人机航拍可以对同一个目标、场景进行多次成像，形成多视角、多位置成像。通过机载高分辨率光学摄像机，无人机航行一次能够获得大量图像序列，这些图像具有高分辨率、高重叠率，并且可以进行量测。应用序列图像三维重建的相关理论和技术，对无人机序列成像的处理，可以获得对地目标的三维模型。

图像三维重建方法的原理是通过摄像机从不同角度获取同一目标的两幅图像，然后利用立体视觉原理，由计算机重建目标的三维形状，恢复出其空间模型。而在航拍图像的三维重建中，由于无人机航拍的抖动、旋转、移动等都会引起航拍相机的震荡，从而使得航拍摄像机的成像模型更加复杂。其中，利用现有图像三维重建方法对树类物体的三维重建，多数只针对树木的骨干结构信息进行建模的，而很少考虑树冠轮廓和内部的纹理细节信息等，往往会导致重建出的树类物体的三维影像失真。此为本发明亟待解决的技术问题。

发明内容

基于现有技术中存在的问题，本发明提出了一种基于无人机航拍序列图像的树木三维重建方法，该完整的三维重建算法，首先通过分水岭分割方法提取树冠的特征区域，然后通过计算RGB颜色空间中的区域相关系数来匹配特征区域并提取匹配特征点对，再结合航拍摄像机的成像模型，通过双目立体视觉原理计算深度信息，最后通过L系统规则构建树木的三维模型。

本发明提出了一种基于无人机航拍序列图像的树木三维重建方法，该方法包括以下步骤：

步骤一、对无人机航拍获得的两张相邻的树木序列图像，从中分别截取大小相同同一树冠图像一、树冠图像二，用grabcut算法去除背景，仅保留树冠区域；

步骤二、对去除背景的树冠图像一，使用分水岭分割算法做分割，得到明暗区域划分，计算每个分割区域的形心作为在树冠图像一上提取的特征点；

步骤三、计算树冠图像一各特征点在RGB颜色空间下的特征区域与树冠图像二的特征区域的相关系数

c_ij＝Rc_ij+Gc_ij+Bc_ij

其中，Rc_ij，Gc_ij，Bc_ij表示以R、G、B三个分量分别计算区域相关系数；

其中，和分别为是树冠图像一、树冠图像二中第i和第j个特征点；wn、wm分别为滑动窗口的半长和半宽，n为窗口中的像素点个数，R₁(u_i，v_i)、G₁(u_i，v_i)、B₁(u_i，v_i)和R₂(u_i，v_i)、G₂(u_i，v_i)、B₂(u_i，v_i)分别为树冠图像一、树冠图像二上点(u_i，v_i)的R、G、B三个分量的值，σ_i、σ_j分别为两幅图像中第i和第j个特征点的图像滑动窗口内的标准差，其中，图像的各分量均值计算如下式：

R分量均值为：

G分量均值为：

B分量均值为：

图像的各分量的图像滑动窗口内的标准差计算如下式：

R分量的图像滑动窗口内的标准差：

G分量的图像滑动窗口内的标准差：

B分量的图像滑动窗口内的标准差：

在树冠图像二中找到树冠图像一特征区域的匹配区域，以匹配区域对的形心作为匹配特征点，通过最相关和次相关比例法去除误匹配点，具体步骤如下：

以树冠图像一的分水岭分割的每个区域的外接矩形作为该区域匹配时的窗口大小；

以同样大小的窗口在树冠图像二上扫描，按上述计算方法计算两幅图像上窗口区域内的RGB区域相关系数，记录得到的最大相关系数和次大相关系数；

计算最大相关系数和次大相关系数的比例k，设定一个k的阈值K，如果k>K，则把该匹配认为是误匹配舍弃，否则，认为该匹配正确并接受；

通过上步得到的两幅树冠图像上的匹配区域对，计算他们的形心作为匹配得到的匹配特征点对。

步骤四、通过计算机立体视觉原理，对图像的特征点集合的一系列计算获得图像特征点的三维信息，包括直接获得图像的特征点的二维参数和计算得到的深度参数，深度参数的计算包括三步，即确定航拍摄像机的成像模型，即摄像机的内外参数；通过摄像机成像模型，确定从二维点到三维点的映射关系；通过映射关系获取图像的三维点集合；

步骤五、从三维点集合获取树木的骨架信息，包括以下处理：

1、找到Z坐标最高点，并且找到最高点和最低点的dz，dz作为树高，把最高点所在的垂直轴作为树的中心主干，其深度为1，加入树干集合TrunkList，并把最高点加入树干点集IncludePoints；

2、扫描不在IncludePoints集合中的特征点，找到当前点集中离树干集合距离最小的点p_i以及相应的树干t_i，限制t_i的深度不超过6，并且，分支数不超过3；

3、从该树干t_i向点p_i引出一个新的树干t_new与ti的夹角为30°，t_new的深度为t_i的深度加1（限制树的最大深度不超过4），加入TrunkList时插入t_i后面，并把点p_i加入已有树干包含的点集合IncludePoints；

4、回到2，直到所有三维散点都在IncludePoints集合中为止；

将上述获得的骨架信息，利用L系统规则，首先根据得到的骨架信息生成L系统语言的编码串，并读入编码串构建树木的骨架结构，并在枝干上加上叶子形成完整的三维树木模型。

与现有技术相比，本发明考虑到树冠的纹理、颜色以及轮廓信息，充分提取出能反映树冠复杂结构的特征点集合。并且，考虑到无人机航拍过程中，摄像机存在的抖动、平移等，针对航拍摄像机复杂的成像模型来计算特征点集合的深度信息。并采取了合理的树木建模方法，构建出树木的近似三维模型。因此，使用本方法，可以从航拍序列图像中，获得树木的三维信息并构建合理的树木三维模型。

附图说明

图1为本发明的基于无人机航拍序列图像的树木三维重建方法的整体流程图；

图2为特征点提取和匹配的流程图；

图3为基于RGB颜色空间下的特征区域匹配方法原理图；

图4为三维L系统的坐标系；

图5为由三维散点集合构建树木三维模型示意图；

图6为从相邻两张无人机航拍图像中截取出的同一个树冠的图像；

图7为对树冠图像应用grabcut方法去除背景的图像；

图8为按基于区域相关系数的区域匹配方法匹配后得到的两幅图像的特征点对集合图；

图9为最相关和次相关比例取不同阈值时获得的特征点对集合以及相应的三维树木骨架结构图，其中，图(a)为K的阈值为0.93时获得的特征点集构建的树木模型，图(b)为K的阈值为0.92时获得的特征点集构建的树木模型，图(c)为K的阈值为0.91时获得的特征点集构建的树木模型，图(d)为K的阈值为0.85时获得的特征点集构建的树木模型；

图10为对经重建的三维骨架结构的枝干上加上树叶的三维树木模型，并与原树冠图像对应的图像示意图，其中，图(a)为树木模型的侧面视图，图(b)为树木模型的俯视图，图(c)为树木的原图像。

具体实施方式

下面结合附图及实施例，进一步详细说明本发明的具体实现方式。

根据本发明的技术方案，在图像获取上，本实施例采用无人机在某山区航拍获得图像，根据树冠大小选择合适的树冠图像大小能包含树冠的全部区域，保证不丢失树冠信息，所以开发者选择128×128的树冠图像，如图6所示。为了尽量消除非树冠区域对重建过程的影响又能不丢失树冠的信息，采用grabcut算法去除树冠图像的背景，如图7所示。

下面对特征点提取和匹配、深度计算以及三维建模三个过程的实施方式进行详细描述：

一、特征点提取和匹配

特征点提取和匹配是三维重建过程中非常重要的一个过程，因为所提取的特征点的优劣和对特征点匹配的准确性直接影响了最终重建出来的三维模型的好坏。

根据图2所示的流程图，本发明的实施分为以下几个步骤：

1、对树冠图像1采用基于梯度的分水岭分割，首先，对树冠图像使用sobel算子进行纵方向的梯度运算；然后，做开闭运算对梯度图像进行平滑处理；最后对平滑后的图像使用分水岭分割方法分割获得图像的特征区域。

2、记录分水岭分割的特征区域信息，以特征区域的外接矩形大小作为滑动窗口大小，在树冠图像2中扫描，计算两幅图像窗口区域的区域相关系数大小，以区域相关系数最大的区域对为最佳匹配区域对，匹配方法如图3所示。

3、提取匹配特征区域对的形心作为两幅树冠图像的匹配特征点，从而获得二维特征点对集合；

4、采用最相关和次相关比例法去除误匹配点，试验不同比例阈值情况下删除的误匹配点的情况，通过对结果的比较选择最合适的阈值来删除误匹配点，如图10所示的示例。

区域相关系数计算如下：

c_ij＝Rc_ij+Gc_ij+Bc_ij

其中，Rc_ij，Gc_ij，Bc_ij表示以R、G、B三个分量分别计算区域相关系数，以R分量为例，计算如下式：

其中，和是两幅图像中第i和第j个特征点；wn、wm分别表示滑动窗口的半长和半宽，n为窗口中的像素点个数，R₁(u_i，v_i)、G₁(u_i，v_i)、B₁(u_i，v_i)和R₂(u_i，v_i)、G₂(u_i，v_i)、B₂(u_i，v_i)分别为图像1、2上点(u_i，v_i)的R、G、B三个分量的值。以R分量为例，图像的分量均值计算如下式：

图像滑动窗口内的标准差：

二、深度计算

例如，采用小孔相机模型，设在该投影模型下三维空间点X＝(X，Y，Z)^T在图像平面上的像点为x^T＝(x，y)^T，则有

其中，和分别表示点x和X的齐次坐标，λ为方程两边齐次向量之间的比例因子，也称为射影深度；K为相机的内参数矩阵，而R和t为外参数，分别为摄像机在世界坐标系下的旋转矩阵和平移向量。

根据摄像机的成像模型，利用双目立体视觉原理来计算特征点的三维坐标，也就是模拟人双眼视觉功能。对两幅图像分别存在一个映射关系，联立两个映射关系得到关于世界坐标(X_w，Y_w，Z_w)的方程组如下：

其中，摄像机的内参数矩阵忽略摄像机存在的微小旋转，所以有对于摄像机的平移向量，取第一张图像拍摄时的位置为参考位置，所以T₁＝(000)^T，T₂＝(T_xT_y0)^T。所以，展开上式并整理得世界坐标系坐标值公式如下：

三、三维建模

首先是从三维散点集合获取树木的骨架枝干信息的方法。设定一个树干的集合TrunkList和这些树干包含的三维散点集合IncludePoints集合。获取骨架信息的方法如下：

2、扫描不在IncludePoints集合中的特征点，找到当前点集中离树干集合距离最小的点p_i以及相应的树干t_i（限制t_i的深度不超过6，并且，分支数不超过3）；

4、回到2，直到所有三维散点都在IncludePoints集合中为止。

获得树木的骨架枝干信息后，利用L系统的规则，采用OpenGL三维图形技术构建出树木的三维模型。

Claims

1.一种基于无人机航拍序列图像的树木三维重建方法，其特征在于，该方法包括以下步骤：

c_ij＝Rc_ij+Gc_ij+Bc_ij

{Rc}_{i j} = \frac{Σ_{k = - w n}^{w n} Σ_{l = - w m}^{w m} {[R_{1} (u_{i}^{(I)} + k, v_{i}^{(I)} + l) - {\overset{&OverBar;}{R}}_{1} (u_{i}^{(I)}, v_{i}^{(I)})] \times [R_{2} (u_{j}^{(I I)} + k, v_{j}^{(I I)} + l) - {\overset{&OverBar;}{R}}_{2} (u_{j}^{(I I)}, v_{j}^{(I I)})]}}{n \sqrt{σ_{i}^{2} (R_{1}) \times σ_{j}^{2} (R_{2})}}

{Gc}_{i j} = \frac{Σ_{k = - w n}^{w n} Σ_{l = - w m}^{w m} {[G_{1} (u_{i}^{(I)} + k, v_{i}^{(I)} + l) - {\overset{&OverBar;}{G}}_{1} (u_{i}^{(I)}, v_{i}^{(I)})] \times [G_{2} (u_{j}^{(I I)} + k, v_{j}^{(I I)} + l) - {\overset{&OverBar;}{G}}_{2} (u_{j}^{(I I)}, v_{j}^{(I I)})]}}{n \sqrt{σ_{i}^{2} (G_{1}) \times σ_{j}^{2} (G_{2})}}

{Bc}_{i j} = \frac{Σ_{k = - w n}^{w n} Σ_{l = - w m}^{w m} {[B_{1} (u_{i}^{(I)} + k, v_{i}^{(I)} + l) - {\overset{&OverBar;}{B}}_{1} (u_{i}^{(I)}, v_{i}^{(I)})] \times [B_{2} (u_{j}^{(I I)} + k, v_{j}^{(I I)} + l) - {\overset{&OverBar;}{B}}_{2} (u_{j}^{(I I)}, v_{j}^{(I I)})]}}{n \sqrt{σ_{i}^{2} (B_{1}) \times σ_{j}^{2} (B_{2})}}

R分量均值为：

\overset{&OverBar;}{R} (u, v) = \frac{1}{n} Σ_{k = - w n}^{w n} Σ_{l = - w m}^{w m} R_{1} (u + k, v + l)

G分量均值为：

\overset{&OverBar;}{G} (u, v) = \frac{1}{n} Σ_{k = - w n}^{w n} Σ_{l = - w m}^{w m} G_{1} (u + k, v + l)

B分量均值为：

\overset{&OverBar;}{B} (u, v) = \frac{1}{n} Σ_{k = - w n}^{w n} Σ_{l = - w m}^{w m} B_{1} (u + k, v + l)

图像的各分量的图像滑动窗口内的标准差计算如下式：

R分量的图像滑动窗口内的标准差：

σ (R) = \sqrt{\frac{1}{n} Σ_{k = - w n}^{w n} Σ_{l = - w m}^{w m} R^{2} (u + k, v + l) - {\overset{&OverBar;}{R}}^{2} (u, v)}

G分量的图像滑动窗口内的标准差：

σ (G) = \sqrt{\frac{1}{n} Σ_{k = - w n}^{w n} Σ_{l = - w m}^{w m} G^{2} (u + k, v + l) - {\overset{&OverBar;}{G}}^{2} (u, v)}

B分量的图像滑动窗口内的标准差：

σ (B) = \sqrt{\frac{1}{n} Σ_{k = - w n}^{w n} Σ_{l = - w m}^{w m} B^{2} (u + k, v + l) - {\overset{&OverBar;}{B}}^{2} (u, v)}

以同样大小的窗口在树冠图像二上扫描，按上述R、G、B分量计算方法计算两幅图像上窗口区域内的RGB区域相关系数，记录得到的最大相关系数和次大相关系数；

通过得到的两幅树冠图像上的匹配区域对，计算他们的形心作为匹配得到的匹配特征点对；

(1)找到Z坐标最高点，并且找到最高点和最低点的dz，dz作为树高，把最高点所在的垂直轴作为树的中心主干，其深度为1，加入树干集合TrunkList，并把最高点加入树干点集IncludePoints；

(2)扫描不在IncludePoints集合中的特征点，找到当前点集中离树干集合距离最小的点p_i以及相应的树干t_i，限制t_i的深度不超过6，并且，分支数不超过3；

(3)从该树干t_i向点p_i引出一个新的树干t_new与ti的夹角为30°，t_new的深度为t_i的深度加1，加入TrunkList时插入t_i后面，并把点p_i加入已有树干包含的点集合IncludePoints；

(4)回到(2)，直到所有三维散点都在IncludePoints集合中为止；