一种基于体元的机载LIDAR建筑物检测方法与流程

本发明涉及遥感数据处理领域，尤其涉及一种基于体元的机载LIDAR建筑物检测方法。

背景技术：
：

建筑物自动检测技术一直是三维(3Dimension，3D)城市建模领域的研究热点。机载激光雷达(LightDetectionAndRanging，LIDAR)点云数据包含了大量关于建筑物的三维信息，十分适用于建筑物检测。经典的建筑物检测方法可分为以下几类：基于拟合的方法、数学形态学方法、数字图像处理方法、模式识别方法及融合LIDAR数据和其它类型的航空影像或GIS数据的方法。上述方法采用的数据结构形式主要有：离散点云数据或栅格数据。点云数据为真3D数据结构，但其难以利用空间邻域信息；栅格数据是将真3D的点云规则化为2.5D数据，其无法表达多次回波数据。可见，经典建筑物检测方法所采用的数据结构均不利于发挥机载LIDAR真3D的技术优势。

技术实现要素：
：

针对现有技术的缺陷，本发明提供一种基于体元的机载LIDAR建筑物检测方法。该方法包括：

步骤1：读取原始机载LIDAR点云数据；

步骤2：将除原始机载LIDAR点云数据规则化为二值3D体元数据集；

步骤2.1：从原始机载LIDAR点云数据中剔除异常数据，得到去除异常数据集；

步骤2.2：将去除异常数据集规则化为二值3D体元数据集；

步骤3：从去除异常数据集中分离出非地面数据集，并映射到3D体元格网，得到非地面体元；

步骤4：用建筑物边缘点的阶跃特性从非地面体元中搜寻建筑物边缘体元作为种子体元，标记与其3D连通的目标体元作为建筑物体元数据集，完成基于体元的机载LIDAR建筑物检测。

进一步地，步骤2.1具体包括：

步骤2.1.1：将原始机载LIDAR点云数据所在空间划分为M×N×U三维格网，并将各原始机载LIDAR点云数据映射到各个格网单元，包含原始机载LIDAR点云数据的格网称为黑格网，不包含原始机载LIDAR点云数据的格网称为白格网；

步骤2.1.2：定位M×N×U三维格网内M×N个立柱内的高程最高与高程最低的黑格网作为候选异常格网单元得到候选异常数据集；

步骤2.1.3：对候选异常数据集中各个候选异常格网单元，比较其和周围给定邻域内黑格网的平均高程的高程差，若高程差大于给定阈值T_ed，则该候选异常格网单元内包含的原始机载LIDAR点云数据为异常数据，进行剔除，否则保留该候选异常格网单元内包含的原始机载LIDAR点云数据，最终获得去除异常数据集。

进一步地，步骤2.2具体包括：

步骤2.2.1：用去除异常数据集的轴向包围盒表示三维空间范围；

步骤2.2.2：计算体元分辨率即体元大小，x、y方向上的体元分辨率(Δx,Δy)依据去除异常数据集中数据点的标称点间距确定，z方向的体元分辨率Δz依据去除异常数据集中数据点的平均点间距确定；

步骤2.2.3：依据x、y、z方向上的体元分辨率对轴向包围盒进行划分得到3D体元格网，每一个3D体元格网单元称为体元；

步骤2.2.4：将去除异常数据集映射到3D体元格网，进而根据3D体元格网中是否包含有去除异常数据分别为3D体元格网中各体元赋值1和0，得到二值3D体元数据集，1和0分别代表目标体元和背景体元，完成对去除异常数据集的规则化。

进一步地，步骤3包括：

步骤3.1：采用不规则三角网的加密滤波算法对去除异常数据集进行滤波得到地面数据集和非地面数据集；

步骤3.2：将地面数据集和非地面数据集分别映射到3D体元格网中，并分别标记为地面体元和非地面体元。

进一步地，步骤4包括：

步骤4.1：依据建筑物边缘点的阶跃特性从非地面体元中搜索建筑物边缘体元作为种子体元V_k；

步骤4.2：依次从V_k的给定邻域尺度内的目标体元出发，深度优先遍历所有邻接目标体元，直至3D体元格网中和V_k有路径连通的所有目标体元均被标记，并将与V_k 3D连通的目标体元作为建筑物体元数据集。

进一步地，步骤4.2包括：

4.2.1：初始化，设置存储种子体元的初始栈，并标记这些种子体元为建筑物体元；

4.2.2：从初始栈栈顶弹出一个元素，获取其给定邻域内未标记的目标体元，标记为建筑物体元并存入初始栈中；

4.2.3：如果初始栈中为空，则3D体元格网中和V_k有路径连通的所有目标体元均被标记；否则，进入步骤4.2.2。

由上述技术方案可知，本发明提出的基于体元的机载LIDAR建筑物检测方法，以3D连通性构建理论为基础，使得点云数据中的目标信息检测从点云聚类等传统方式转换成基于体元空间邻域关系的搜索标记方式，很好的利用了3D体元数据中各体元间隐含的邻域关系，有助于基于体元理论的机载LIDAR点云数据处理及应用的发展。实验基于ISPRS提供城区机载LIDAR点云数据定量评价计算精度，对大型、密集、不规则形状及其它屋顶类型比较特殊的复杂建筑物的检测结果其完整度可达到90％以上，正确率可达85％以上，可有效实现对建筑物的检测。

附图说明：

图1为本发明实施例提供的基于体元的机载LIDAR建筑物检测方法流程图；

图2为本发明具体实施例中步骤2的具体流程图；

图3为本发明具体实施例中步骤2中的格网边长计算原理图；

图4为本发明具体实施例中步骤2中的平均点间距计算原理图；

图5为本发明具体实施例中步骤3的具体流程图；

图6为本发明具体实施例中步骤4的具体流程图；

图7为本发明具体实施例中步骤4中的建筑物边缘点的特性示意图；

图8为本发明具体实施例中步骤4中的邻域尺度示意图，其中，(a)为6邻域，(b)为18邻域，(c)为26邻域，(d)为56邻域。

具体实施方式：

下面结合附图和实施例，对本发明的具体实施方式作进一步详细描述。以下实施例用于说明本发明，但不用来限制本发明的范围。

图1示出了一种基于体元的机载LIDAR建筑物检测方法流程图，包括：

步骤1：读取原始机载LIDAR点云数据；

本实施例中，定义原始机载LIDAR点云数据集P＝{p_i(x_i,y_i,z_i),i＝1,…,n},其中，i是原始机载LIDAR点云数据的索引，n是原始机载LIDAR点云数据的个数，p_i是第i个原始机载LIDAR点云数据，其坐标为(x_i,y_i,z_i)。

步骤2：将原始机载LIDAR点云数据中规则化为二值3D体元数据集；

步骤3：从去除异常数据集中分离出非地面数据集，并映射到3D体元格网，得到非地面体元；

步骤4：用建筑物边缘点的阶跃特性从非地面体元中搜寻建筑物边缘体元作为种子体元，标记与其3D连通的目标体元作为建筑物体元数据集，完成基于体元的机载LIDAR建筑物检测。

图2示出了基于体元的机载LIDAR建筑物检测方法步骤2的具体流程图，包括：

步骤2.1：从原始机载LIDAR点云数据中剔除异常数据，得到去除异常数据集；

步骤2.2：将去除异常数据集规则化为二值3D体元数据集。

其中，步骤2.1包括：

步骤2.1.1：将原始机载LIDAR点云数据划分为M×N×U三维格网，并将各原始机载LIDAR点云数据映射到各个格网单元，包含原始机载LIDAR点云数据的格网称为黑格网，不包含原始机载LIDAR点云数据的格网称为白格网；

(1)计算原始机载LIDAR点云数据集P的轴向包围盒。轴向包围盒由其左下角坐标(x_min,y_min,z_min)和右上角坐标(x_max,y_max,z_max)确定，其中，(x_max,y_max,z_max)＝max{(x_i,y_i,z_i),i＝1,2,…,n}，(x_min,y_min,z_min)＝min{(x_i,y_i,z_i),i＝1,2,…,n}分别代表数据集P中x、y和z坐标的最大值和最小值，其中，i是原始机载LIDAR点云数据的索引，n是原始机载LIDAR点云数据的个数。

(2)将轴向包围盒划分为三维格网，其中格网长度取原始机载LIDAR点云数据平均点间距d：

$d=\sqrt{\frac{A\left(C\right(S_{xy}))}{n}}---\left(1\right)$

如图3所示，点集S_xy＝{(x_i,y_i),i＝1,…,n}为数据集P＝{p_i(x_i,y_i,z_i),i＝1,…,n}在XOY平面上的投影，其中，C(S_xy)是点集S_xy的凸壳，A(C(S_xy))是凸壳C(S_xy)的面积，由此可得M×N×U个三维格网单元，M、N、U由公式(2)确定。

其中，为向下取整操作符。

(3)将原始机载LIDAR点云数据映射到各个格网单元，获取各原始机载LIDAR点云数据的格网索引，包含原始机载LIDAR点云数据的格网被称为黑格网，不包含原始机载LIDAR点云数据的格网被称为白格网。

将各数据映射到各个格网单元：

其中，(m_i,n_i,u_i)代表原始机载LIDAR点云数据p_i所在格网单元的索引。

步骤2.1.2：定位M×N×U三维格网内M×N个立柱内的高程最高与高程最低的黑格网作为候选异常格网单元得到候选异常数据集；

步骤2.1.3：对候选异常数据集中各个候选异常格网单元，比较其和周围给定邻域内黑格网单元的平均高程的高程差，若高程差大于给定阈值T_ed，则候选异常格网单元内包含的原始机载LIDAR点云数据为异常数据，进行剔除，否则保留候选异常格网单元内包含的原始机载LIDAR点云数据，最终获得去除异常数据集。

在本实施方式中，T_ed为常数(如1米)，其取值需根据原始机载LIDAR点云数据的空间分布情况确定。

其中，步骤2.2包含以步骤：

步骤2.2.1：用去除异常数据集的轴向包围盒表示三维空间范围；

去除异常数据集记做Q＝{q_i'(x_i',y_i',z_i'),i'＝1,…,t}，其中，i'是去除异常数据集中数据的索引，t是去除异常数据集中数据的个数，q_i'是去除异常数据集中第i'个数据，其坐标为(x_i',y_i',z_i')。用Q的轴向包围盒表示三维空间范围，轴向包围盒的确定参见步骤2.1.1；

步骤2.2.2：计算体元分辨率即体元大小，体元x、y方向上的分辨率(Δx,Δy)依据去除异常数据集中数据点的标称点间距确定，体元z方向的分辨率Δz依据去除异常数据集中数据点的平均点间距确定；

体元x、y方向的分辨率(Δx,Δy)依据去除异常数据集中数据点的标称点间距确定：

S'_xy＝{(x_i',y_i'),i'＝1,…,t}

Δx_i'＝min{|x_i'-x_i”|；i”＝1,…,t,i”≠1}

Δy_i'＝min{|y_i'-y′_i”|；i”＝1,…,t,i”≠1}

P_sx＝arg{#{Δx_i'<P_x；i'＝1,…,t}＝0.95t}

P_sy＝arg{#{Δy_i'<P_y；i'＝1,…,t}＝0.95t}

Δx＝Δy＝max{P_sx,P_sy} (4)

体元z方向的分辨率Δz依据去除异常数据集中数据点的平均点间距确定；

$\mathrm{\&Delta;}z=min(\sqrt{\frac{A\left(C\right(S_{xz}^{\&prime;}))}{t}},\sqrt{\frac{A\left(C\right(S_{yz}^{\&prime;}))}{t}})---\left(5\right)$

其中，S'_xy＝{(x_i',y_i'),i'＝1,…,t}，S'_xz＝{(x_i',z_i'),i'＝1,…,t}，S'_yz＝{(y_i',z_i'),i'＝1,…,t}。

如图4所示，C(S'_xz)是点集S'_xz的凸壳，C(S'_yz)是点集S'_yz的凸壳，其中，S'_xy、S'_xz、S'_yz分别为去除异常数据集Q在XOY、XOZ、YOZ平面上的投影，A(C(S'_xz))、A(C(S'_yz))是是凸壳C(S'_xz)、C(S'_yz)的面积。式(5)中取最小值是因为其代表所建立的3D体元格网和原始点云间存在更少的精度损失。

步骤2.2.3：依据x、y、z方向上的体元分辨率对轴向包围盒进行划分得到3D体元格网，每一个3D体元格网单元称为体元；

基于体元分辨率(Δx,Δy,Δz)就可以将轴向包围盒划分为3D体元格网，用3D体元阵列表示。设V是3D体元阵列中的体元集合，

V＝{v_j(r_j,c_j,l_j),j＝1,…,m}, (6)

其中，j是体元索引；m是体元数；v_j是第j个体元的体元值；(r_j,c_j,l_j)是第j个体元在体元阵列中的坐标(行、列和层号)。X方向上的体元个数为R，Y方向上的体元个数为C，Z方向上的体元个数为L。其中，R、C、L由公式7得出；

其中，为向上取整操作符。

由此可以得出，体元数m为：

m＝R*C*L (8)

步骤2.2.4：将异常数据剔除后的去除异常数据集映射到3D体元格网，进而根据3D体元格网中是否包含有去除异常数据分别为3D体元格网中各体元值赋值为1和0，得到二值3D体元数据集，1和0分别代表目标体元和背景体元，得到二值3D体元数据，完成对去除异常数据集的规则化。

体元赋值由公式9得出：

其中，为向下取整操作符。

图5示出了基于体元的机载LIDAR建筑物检测方法步骤3的具体流程图，包含以下步骤：

步骤3.1：采用不规则三角网(Triangulated IrregularNetwork，TIN)的加密滤波算法对去除异常数据集进行滤波得到地面数据集和非地面数据集；

首先选取去除异常数据集中的一些高程最低点为初始点构建初始不规则三角形，若去除异常数据集中的点到最近三角形的距离以及去除异常数据集中的点与最近三角形顶点的连线与该三角形的夹角均小于给定的阈值，则将该去除异常数据集中的点加密到该三角网中，并标记为地面数据，通过迭代加密过程，直到没有新点加入三角网时运算终止。完成从去除异常数据集中分离出地面数据集和非地面数据集。

步骤3.2：将地面数据集和非地面数据集分别映射到二值3D体元格网中并相应标记为地面体元(如标记为“G”)和非地面体元(如标记为“NG”)。

图6示出了基于体元的机载LIDAR建筑物检测方法步骤4的具体流程图，包括：

步骤4.1：依据建筑物边缘点的阶跃特性从非地面体元中搜索建筑物边缘体元作为种子体元V_k；

步骤4.2：依次从V_k的给定邻域尺度内的目标体元出发，深度优先遍历所有邻接目标体元，直至二值3D体元格网中和V_k有路径连通的所有体元均被标记，并将与V_k 3D连通的目标体元作为建筑物体元数据集。

其中，步骤4.1的具体步骤如下：

建筑物边缘点的特性为：在其水平邻域内与局部最低点的高程值有明显的阶跃变化，且其空间最近邻接点仅覆盖了其周围一半的区域。根据以上特性建立从非地面体元中搜寻建筑物边缘体元的判别方案：

4.1.1：对在V中搜寻k个与v_j欧氏距离最近的非地面体元，记作N_j＝{n_s(r_s,c_s,l_s),s＝1,…,k}，如图7所示，将v_j及N_j投影至XOY平面得二维离散点集，对应记做V_pj(r_j,c_j)及N_pj＝{n_ps(r_s,c_s),s＝1,…,k}。若N_pj中相邻的两点与V_pj的连线的夹角大于给定的阈值，如120°，则判定v_j可能位于建筑物边缘上；否则，判定v_j必然不位于建筑物边缘上。

其中，n_s为V中与v_j欧氏距离最近的非地面体元，s为n_s的体元索引，k为n_s的个数，(r_s,c_s,l_s)为n_s的体元坐标，N_j为n_s的集合。

4.1.2：对可能位于建筑物边缘上的体元v_j，在V中搜寻c个与其水平邻域最近的目标体元N_j′＝{n_s′(r_s,c_s,l_s),s＝1,…,c}，如图7所示，若v_j与N′_j中高程最低体元的高差大于给定的阈值，如3米，则初步判定v_j位于建筑物边缘上。初步判定该体元v_j为建筑物边缘体元，并作为种子体元V_k。

步骤4.1.3：剔除种子体元中的非建筑物边缘体元。

具体方法为：确定的建筑物边缘体元会存在部分错误判定，错误判断受原始机载LIDAR点云数据点的密度和分布的影响，例如部分植被被误判为建筑物边缘。剔除错误判断的建筑物边缘体元方案如下：首先，统计建筑物边缘体元的一定邻域(如5×5×5)内目标体元个数，若个数小于等于1个，则剔除该建筑物边缘体元，否则，统计该建筑物边缘体元与其周围目标体元的高程差，若大于某一阈值(如2×Δz)则剔除。

图8(a)示出了本实施例中步骤4.2所述建筑物体元集合的构建规则，对任意种子体元V_k，依次从V_k的6邻域内(或18、26、56等其它邻域尺度，其中，图8(b)、(c)、(d)依次为18邻域、26邻域、56邻域)的目标体元出发，深度优先遍历所有邻接目标体元，直至二值3D体元阵列中和V_k有路径相通的所有目标体元均被标记。详细步骤如下：

4.2.1：初始化，设置存储种子体元的初始栈，并标记这些种子体元为建筑物体元(如标记为“B”)；

4.2.2：从初始栈栈顶弹出一个元素，获取其给定邻域内未标记的目标体元，标记为建筑物体元并存入初始栈中；

4.2.3：如果初始栈中为空，则3D体元格网中和V_k有路径相通的所有目标体元均被标记；否则，进入步骤4.2.2。

在标记过程中应用不同的邻域尺度会得到不同的建筑物提取结果。邻域尺度过小会导致建筑物检测不完整；邻域尺度过大则影响效率及准确性。最佳邻域尺度在实验中确定。

本发明提出的建筑物检测方法是基于体元表示的，而计算机数据库中提供的数据则是离散的LIDAR激光点云数据，为和计算机数据库中提供的数据做对比以评价本发明提出的方法精度，首先统计本方法所检测的建筑物体元内包含的原始机载LIDAR点云数据的个数，然后和计算机数据库中提供的数据进行比对进而用I类误差(将建筑物脚点错分为非建筑物脚点比例)、II类误差(将非建筑物脚点错分为建筑物脚点比例)、总误差(错分的建筑物脚点的比例)、正确率(正确检测的建筑物脚点数占检测结果中建筑物脚点总数的比例)和完整度(正确检测的建筑物脚点数占标准数据中建筑物脚点总数的比例)来定量评价本发明所提出的建筑物检测方法的有效性。

本发明可以在CPU Core(TM)i5-24003.10GHz、内存4GB、Windows 7旗舰版系统上使用MATLAB 7.11.0平台编程实现该方法，并进一步通过对该方法的精度评定验证方法的有效性。

本实施例中采用国际摄影测量与遥感协会(International Society for Photogrammetry and Remote Sensing，ISPRS)第三工作组提供的(http://www.itc.nl/isprswgIII-3/filtertest/)专门用于滤波算法测试的城区样本数据作为实验数据，以检验方法的有效性和可行性。样本数据由OptechALTM机载激光扫描仪获取，覆盖Stuttgart和Vangen/Enz市中心区域，包括了不同的建筑物类型，样本数据见表1。

本实施例中，采用手工分类的方式从参考数据(被准确分类为地面点集和非地面点集的样本数据，由网站提供)的非地面点集中分离出建筑物点集，将该数据作为建筑物检测的标准数据，以评价本发明方法的计算精度。

表1样本数据的特性及基本参数

表2为本实施例中邻域尺度为6、8、26、56以及80时，对7个样本数据进行检测，对应的建筑物检测结果的总误差。该表中的数据旨在考查不同领域尺度对建筑物检测结果的影响。

表2不同邻域尺度的建筑物检测结果的总误差对比(％)

由表2可知，从总误差指标来看，6、18、26、56及80邻域的平均总误差分别为32.60％、17.43％、14.80％、4.82％及6.67％。同时，由于实验数据中包含了大型、不规则形状及其它屋顶类型比较特殊的复杂建筑物，因而可得出结论：56邻域是该建筑物检测方法的最佳邻域尺度。

表3为本实施例中，以手工分类结果为标准对7个样本数据的56邻域尺度下的建筑物检测精度进行的定量评价。可以看出：(1)I类误差、II类误差和总误差等指标的误差范围分别为1-9％,0-12％,0-8％，这表明本发明表现为最小化总误差；(2)从总误差指标看，样本11总误差最大(7.24％)，该地区误差大的原因有二：由于浓密的植被分布，有一处植被边缘被错判为建筑物边缘(即种子体元)，从而在检测结果中引入了II类误差，该项误差可通过对建筑物检测结果进行后处理消除掉，如依据“分割后单个建筑物的点云密度或面积”等条件；其二是建筑物和植被相邻构成了3D连通集合，因而植被被误判为建筑物，从而在检测结果中引入了II类误差。(3)对大型、密集、不规则形状及其它屋顶类型比较特殊的复杂建筑物完整度可达到90％以上，正确率可达85％以上。从而验证了本发明提出的方法的有效性。

表356邻域尺度下的建筑物检测精度

本发明提供的基于体元的机载LIDAR建筑物检测方法，以3D连通性构建理论为基础，使得点云数据中的目标信息检测从点云聚类等传统方式转换成基于体元空间邻域关系的搜索标记方式，很好的利用了3D体元数据中各体元间隐含的邻域关系，有助于基于体元理论的机载LIDAR点云数据处理及应用的发展。本方法对大型、密集、不规则形状及其它屋顶类型比较特殊的复杂建筑物的检测结果其完整度可达到90％以上，正确率可达85％以上，可有效实现对建筑物的检测。

最后应说明的是：以上各实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述各实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分或者全部技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明权利要求所限定的范围。