CN102930540A - 城市建筑物轮廓检测的方法及系统 - Google Patents

Abstract

本发明公开了一种城市建筑物轮廓检测的方法及系统，其中方法包括以下步骤：S1、对LiDAR点云数据进行处理提取建筑物的初始轮廓；S2、将LiDAR点云数据与航空遥感影像数据进行数据配准；S3、根据数据配准的结果进行影像定向；S4、根据影像定向，检测影像边缘；S5、在初始轮廓周围建立缓冲区，并在该缓冲区内搜索影像边缘，利用几何约束条件识别像方轮廓；S6、对像方轮廓进行轮廓矢量化与规则化，获得最终建筑物轮廓。本发明能够精确检测大范围的城市建筑物轮廓。

Description

城市建筑物轮廓检测的方法及系统

技术领域

本发明涉及图像处理领域，尤其涉及一种城市建筑物轮廓检测的方法及系统。

背景技术

单纯利用遥感影像很难高质量高效率的实现城市建筑物轮廓检测，因为：1）由于遥感卫星的传感器成像过程非常复杂，致使建筑物轮廓在影像中不再是理想的线条，其宽度、灰度是不一致的，边界有时是模糊的、断裂的；2）传感器的摄影比例尺和影像分辨率决定了建筑物轮廓在影像上的几何特征；3）较高的建筑物容易遮挡其旁边的矮层，使与其邻近的建筑物信息完全损失。较高的树木容易使建筑物轮廓在影像上表现模糊，或是与其它灰度相似的地物混淆，尤其会在轮廓识别过程中造成道路断裂；4）屋顶面的其他目标，如空调、管线、围栏等会使得轮廓的灰度分布不均匀。

LiDAR（Light Detection And Ranging激光雷达）数据的空间分辨率依赖于几个因素，如飞行高度、速度和扫描频率.现阶段，空载LIDAR数据的点间距一般为米级，相对于摄影测量常用的厘米级分辨率航空影像，LIDAR数据空间分辨率低，单纯基于LIDAR获取的建筑物轮廓相对不够精细，且精度不高。利用规则化技术可以提高轮廓线的精度，但规则化方法也受制于LIDAR数据的空间分辨率、滤波与分类等处理的精度，还存在不确定性且容易产生轮廓整体偏移问题。LiDAR数据以其较高空间频率提供了精确的点坐标，然而从中无法获得断裂线。在城区，绝大多数的断裂线是建筑物的边界。因此，建筑物精确边界的确定是重点和难点。由LiDAR获得的不明确的建筑物轮廓会导致的正射影像和数字地形图建筑物边缘不精确，出现锯齿现象。

从影像中提取出的轮廓几何精度高、细节丰富，但处理的自动化程度不高，从LiDAR数据中提取的轮廓精确性与细节性不高，但LiDAR直接提供了3维坐标，有利于提高处理的自动化程度。

发明内容

本发明要解决的技术问题在于针对现有技术中无法精确确定建筑物边界的缺陷，提供一种能够精确检测大范围建筑物边界的城市建筑物轮廓检测的方法及系统。

本发明解决其技术问题所采用的技术方案是：

提供一种城市建筑物轮廓检测的方法，包括以下步骤：

S1、对LiDAR点云数据进行处理提取建筑物的初始轮廓；

S2、将所述LiDAR点云数据与航空遥感影像数据进行数据配准；

S3、根据数据配准的结果进行影像定向；

S4、根据所述影像定向，检测影像边缘；

S5、在所述初始轮廓周围建立缓冲区，并在该缓冲区内搜索所述影像边缘，利用几何约束条件识别像方轮廓；

S6、对所述像方轮廓进行轮廓矢量化与规则化，获得最终建筑物轮廓。

本发明所述的城市建筑物轮廓检测的方法中，步骤S1包括：

S11、采用线性预测法分离所述LiDAR点云数据；

S12、根据分离结果并利用扫描带法分割建筑物点云；

S13、根据分割好的建筑物点云，利用轮廓规则化分析方法分离树丛；

S14、根据分割好的建筑物点云以及分离的树丛，进行建筑物轮廓检测，得到建筑物的初始轮廓。

本发明所述的城市建筑物轮廓检测的方法中，在进行建筑物轮廓检测时，先确定检测轮廓的搜索方向，建筑物轮廓方向与所述搜索方向的夹角小于树丛轮廓方向与所述搜索方向的夹角。

本发明所述的城市建筑物轮廓检测的方法中，步骤S1还包括：

S15、对得到的建筑物初始轮廓进行直线段提取，得到优化后的建筑物初始轮廓。

本发明所述的城市建筑物轮廓检测的方法中，所述直线段提取具体为：

以预设的第一距离阈值搜索符合直线要求的直线片段；

对直线片段采用预设的角度阈值和第二距离阈值合并邻接直线片段为直线段。

本发明所述的城市建筑物轮廓检测的方法中，所述第一距离阈值为小于等于1.5个像素；所述第二距离阈值为4个像素；所述角度阈值为0.5弧度。

本发明解决其技术问题所采用的另一技术方案是：

提供一种城市建筑物轮廓检测的系统，包括：

建筑物初始轮廓提取模块，用于对LiDAR数据点云进行处理提取建筑物的初始轮廓；

配准模块，用于将所述LiDAR点云数据与航空遥感影像数据进行数据配准；

影像定向模块，用于根据数据配准的结果进行影像定向；

影像边缘检测模块，用于根据所述影像定向，检测影像边缘；

像方轮廓识别模块，用于在所述初始轮廓附近建立缓冲区，并在该缓冲区内搜索所述影像边缘，利用几何约束条件识别像方轮廓；

建筑物轮廓获得模块，用于对所述像方轮廓进行轮廓矢量化与规则化，获得最终建筑物轮廓。

本发明所述的城市建筑物轮廓检测的系统中，所述建筑物初始轮廓提取模块还包括：

点云数据分离模块，用于采用线性预测法分离所述LiDAR点云数据；

建筑物点云分割模块，用于根据分离结果并利用扫描带法分割建筑物点云；

树丛分离模块，用于根据分割好的建筑物点云，利用轮廓规则化分析方法分离树丛；

建筑物轮廓检测模块，用于根据分割好的建筑物点云以及分离的树丛，进行建筑物轮廓检测，得到建筑物的初始轮廓。

本发明所述的城市建筑物轮廓检测的系统中，所述建筑物初始轮廓提取模块还包括：

初始轮廓优化模块，用于对所述建筑物轮廓检测模块得到的建筑物初始轮廓进行直线段提取，得到优化后的建筑物初始轮廓。

本发明所述的城市建筑物轮廓检测的系统中，所述建筑物轮廓检测模块还用于在进行建筑物轮廓检测时，先确定检测轮廓的搜索方向，建筑物轮廓方向与所述搜索方向的夹角小于树丛轮廓方向与所述搜索方向的夹角。

本发明产生的有益效果是：本发明通过利用LiDAR数据易于自动化检测目标的特点，获得建筑物初始轮廓，再在初始轮廓缓四周建立缓冲区内，并在缓冲区内搜索符合实际的影像边缘，利用几何约束条件识别像方轮廓；对所述像方轮廓进行规则化最终获得精确、完整的建筑物轮廓。本发明可提高建筑物轮廓检测的自动化程度，可以实现大范围快速检测；且提取的建筑物轮廓具有最好的质量，避免了噪声和影像质量导致的算法局部极值和迭代死循环问题，可有效的保证算法的稳健性；另外，本发明对于LiDAR数据和遥感影像的来源、格式无特殊要求，具有很大的普适性。

附图说明

下面将结合附图及实施例对本发明作进一步说明，附图中：

图1是本发明实施例的城市建筑物轮廓检测的方法流程图；

图2是图1中步骤S1的具体流程图一；

图3是图1中步骤S1的具体流程图二；

图4a是某地区密度为1点/m²激光扫描数据激光扫描数据图；

图4b是和图4a中相同地区的遥感影像图；

图5是对图4a进行了建筑物点云分割后的示意图；

图6a是图4a中建筑物的初始轮廓图；

图6b是图6a中虚线方框内的局部放大图；

图6c是图6b中局部填充图；

图7a是初步像素聚类产生的边缘片段

图7b是直线片段聚类之后的轮廓线段

图7c是消除过聚类和若聚类之后的结果

图8a是图6a中提取直线段像素后的局部示意图；

图8b是图8a中直线段经过二步法优化后的示意图；

图9a是本发明实施例在初始轮廓周围建立缓冲区的示意图；

图9b是与图9a对应区域的影像边缘示意图；

图9c是经几何约束条件识别的像方轮廓；

图10是本发明实施例中像素p的3种可能搜索方向的示意图；

图11a是图4a与图4b进行数据配准后的LiDAR点云图的配准点示意图；

图11b是图4a与图4b进行数据配准后的航空遥感影像图的配准点示意图；

图12是本发明实施例城市建筑物轮廓检测的系统结构示意图；

图13是图12中建筑物初始轮廓提取模块10的具体结构示意图一；

图14是图12中建筑物初始轮廓提取模块10的具体结构示意图二。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅用以解释本发明，并不用于限定本发明。

如图1所示，本发明实施例城市建筑物轮廓检测的方法包括以下步骤：

S1、对LiDAR点云数据进行处理提取建筑物的初始轮廓；

S2、将LiDAR点云数据与航空遥感影像数据进行数据配准；

S3、根据数据配准的结果进行影像定向；

S4、根据影像定向，检测影像边缘；

S5、在初始轮廓周围建立缓冲区，并在该缓冲区内搜索影像边缘，利用几何约束条件识别像方轮廓；

S6、规则化像方轮廓进行获得建筑物轮廓。

在本发明的一个实施例中，如图2所示，步骤S1具体包括：

S11、采用线性预测法分离LiDAR点云数据；

将测区分割为等间隔的小方块，方格大小依据点的密度和分布特征可预先设定，在处理每一个目标方格单元时，同时考虑其8领域小方格的数据点，那么就能根据9个相邻单元格的数据点利用最小二乘法计算该区域的趋势面系数。拟合的趋势面表达为：

Z_i＝a₀+a₁X_i+b₂Y_i                                  (1)

每个数据点高程的最小二乘平差与数据点高程的拟合残差为：

l_i＝Z_i-obs_i                                       (2)

假设l_i呈正态分布，且标准偏差为σ_z，容忍因子为k，则设定阈值为：

T_z＝kσ_z                                          (3)

如果单元方格内数据点的残差大于阈值，则该点为粗差，在进行下一次计算时被过滤，重复上述操作，直到剔除所有粗差和趋势面稳定。线性最小二乘的协方差函数由下式定义：

$C(P_i,P_k)=C\left(0\right)\·e^{-1.30103{\left(\frac{P_i{P}_k}{B}\right)}^2}---\left(4\right)$

式中，两点P_i、P_k间的协方差取决于两点间的间距，间距越小，协方差值越大，反之越小；参数C(0)为信号协方差函数的峰值，它表示高程偏差中系统误差和随机误差成分的比例关系，如果C(0)＝1.0，就认为没有随机误差，一般的C(0)＝0.99；参数B为坡度参数因子。协方差函数的峰值越大，方差越小，滤波效果就越若。

P_i点的线性预测值由下式估计：

u_i＝c^T·C^-1·z                                (5)

式中：

c^T＝(C(P_i,P₁),C(P_i,P₂),…,C(P_i,P_n))

z＝[l₁,l₂,…,l_n]^T

协方差矩阵的对角线元素为残差，假定所有观测值等权，即V_zz＝1.0。该方法要进行两次阈值处理：①比较l_i和T_z；②比较z_i和l_pre。

可以理解的是也可以采用其他方式分离LiDAR点云数据，如平面拟合法，坡度法、聚类方法等。

S12、根据分离结果并利用扫描带法分割建筑物点云；

通过上述分离LiDAR点云数据，可以将LiDAR点云数据分为地面点云和非地面点云。对于离散的非地面点云，在扫描线上会出现分散的片段，可由此区分并标示各个建筑物；在相邻扫描线上的属于同一建筑物的片段必然是彼此交叉的，由此可获得同一建筑物的所有点云。

利用扫描带法分割建筑物点云具体又包括以下步骤：

1）分割相邻扫描线

设定分割阈值T，对某一行C_i，扫描获得非地面点序列U_i＝[u_i1,u_i2,…,u_in]，若存在：

u_ij-u_i(j-1)＞T                                  (6)

则获得分割片段s，以u_i(j-1)为前一片段的终点，以u_ij为后一片段的起点。要特别注意片段在行边界的情况。最终，扫描获得片段集合S_i＝[s_i1,s_i2，…,s_in]，同样的，行C_i+1扫描获得片段集合S_i+1＝[s_(i+1)1，s_(i+1)2,…,s_(i+1)n]。C_i若为初始行，则S_i中片段以其自然顺序编号，否则C_i中片段通过交叉判断获得全局编号，而C_i+1目前按自然顺序编号，需要通过与C_i的交叉判断获得全局编号。此处，全局编号为建筑物编号，即最终分割结果的表示方法。

2)片段分析

为了获得一个建筑物的所有点云，需要判断相邻扫描线C_i，C_i+1中片段是否交叉，并依此使C_i+1中片段获得全局编号。设当前全局编号为Sum，C_i中片段数为N_i，C_i+1中片段数为N_i+1，其中Sum包含对C_i中N_i个片段交叉判断的结果。首先分析扫描行内片段的存在情况对交叉判断的影响，分为以下几种情况：

①Sum＝0，则自然N_i＝0，

若N_i+1＝0；即从开始到C_i+1的扫描都为空，不需要进行交叉判断，也不需要编号；

若N_i+1＞0；从开始到C_i的扫描都为空，不需要进行交叉判断，C_i+1中片段编号以其自然序号即可；

②Sum＞0，

若N_i+1＝0，此时不需要进行交叉判断，也不需要编号。

若N_i+1＞0，且N_i＝0；即C_i中片段数为0，此时不需要进行交叉判断，C_i+1中片段编号在Sum上加以其自然序号即可；

若N_i+1＞0，且N_i＞0；此时需要交叉判断，C_i+1中片段编号以交叉结果排序；

根据以上分析，从理论上，仅当N_i+1＞0，且N_i＞0时，需要交叉判断，然而这是实际中最多的情况。

3)交叉判断

设C_i中某片段S_ij长度为l₁，C_i+1中片段S_(i+1)l长度为l₂，计算S_ij起点到S_(i+1)l起点距离d₁，以及S_ij终点到S_(i+1)l终点距离d₂，则可得两片段交叉的判C_i+1决式：

|d₁|+|d₂|≤l₁+l₂                              (7)

对C_i中每个片段，逐个判断与C_i+1中的交叉情况，若满足式（7），则认为得到其交叉片段。设当前片段为S_ij，全局编号为N₁，其交叉片段S_(i+1)l的编号为N₂，则需要分析以下几种情况以确定S_(i+1)l，S_ij以及S_ij向前交叉的片段的全局编号。

①S_(i+1)l未被标记，即N₂为S_(i+1)l在C_i+1中的自然序号，则将S_ij的全局编号赋予S_(i+1)l，同时标记S_(i+1)l已交叉；

②S_(i+1)l已被标记，即N₂为S_(i+1)l全局编号，此时需要比较N₁、N₂大小，根据比较结果分三种情况处理：

若N₁＞N₂，将N₂赋予S_ij扫描C_i之前所有行，对于所有编号大于N₁的像素编号值减1，对于所有编号等于N₁的像素，编号修改为N₂，总编号Sum减1；

若N₁＜N₂，将N₁赋予S_(i+1)l，扫描C_i之前所有行，对于所有编号大于N₂的像素编号值减1，对于所有编号等于N₂的像素，编号修改为N₁，Sum减1；

若N₁＝N₂，则不需要特别处理。

若C_i+1中存在未被标记片段，即存在新的点云区域的起点，则其编号应在Sum基础上排序。

4)定位分割区域

从第一行开始，循环运行步骤1）-3），直到最大行，则实现对图像所有像素的分割和标示。

设分割区域的初始参数为P_i＝[0,rows，cols,0,0]，其中rows，cols分别表示图像的最大行、列；P_i(1)表示分割区域的像素数，即该区域的面积，这是去噪和考量区域重要性的主要依据；P_i(2),P_i(3)为区域在图像坐标系左上角坐标，P_i(4),P_i(5)为区域在图像坐标系右下角坐标；P_i(2)-P(5)标示了区域的最小外接矩形（MBR），是定位分割区域的重要参数。

对完成循环的图像Seg逐个像素扫描，若Seg(i,j)＝n，即像素Seg(i,j)属于第n个分割区域，则有：

P_n(1)＝P_n(1)+1;P_n(2)＝min(P_n(2),i);P_n(3)＝min(P_n(3),j);

(8)

P_n(4)＝max(P_n(4),i);P_n(5)＝max(P_n(5),j)

最后可得到每个分割区域的像素数，MBR参数，完成对分割区域的统计和定位。

本发明实施例以图4a所示的某地区密度点位1点/m²激光扫描数据为例，该地区的建筑物数量较大，其结构和拓扑关系较复杂，具有一定的典型性。对其利用扫描带法分割建筑物点云后得到如图5所示的点云分割图。

可以理解的是，也可以采用其他方式分割建筑物点云，如聚类法、坡度法，平面拟合法等。

S13、根据分割好的建筑物点云，利用轮廓规则化分析方法分离树丛；

本发明实施例通过轮廓规则化分析的方法，利用建筑物轮廓与树丛轮廓的极大差异性分离。建筑物轮廓基本由垂直、平行直线段构成，有较强的规律性可循，而树丛轮廓缺少连续规则直线段。轮廓规则化参数为：

$\mathrm{\τ}=\frac{\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{s}{\mathrm{\Σ}}}{l}_{\mathrm{ij}}}{\underset{k=1}{\overset{n}{\mathrm{\Σ}}}{l}_k}---\left(9\right)$

其中，n为区域的所有线段数，l_ij为第i组平行线段中第j条线段长度，m为平行线段组数，s为该组线段数。因为建筑物极具规则性，τ的阈值T_τ可设置较大，当τ＞T_τ，判定区域为建筑物，反之为树丛，本发明T_τ取0.8。为了使式（9）免受分割噪声影响，规定若边界线段集合L＝(l₁,l₂，，，l_n)满足max(L)＜20，区域判定为树丛，不需采用式（9）判定；考虑到正多边形建筑物（塔）存在，若std(L)≈0，区域也判定为建筑。

式（9）仅适合于分离与建筑物离散的树丛，若树丛与建筑物相连，则在区域轮廓上会表现为连续不规则短小线段，此时从区域轮廓分割出这些线段，以上述方法判定这些线段包围区域属性，必要时可利用该区域的高度统计信息判断。

分离树丛也可以采用其他的方法，如：多次回波差值法。

S14、根据分割好的建筑物点云以及分离的树丛，进行建筑物轮廓检测，得到建筑物的初始轮廓。

各种边缘算子，着重在灰度急剧变化，边缘点在较高目标边缘和在较低目标边缘的可能性是一样的，因此上述算子中最好的检测结果均在50%附近。由于在线性预测分离地面点和非地面点时，已采用最小二乘法区分了目标间的高差，因此保留的区域符合当前目标的高度，其轮廓即为区域最外围的点。对区域内任意像素p，若编号为PN，即f(p)＝PN，其4邻域为A，则其是否轮廓的判定式为：

$P\left(\mathrm{\τ}\right)=\frac{{\mathrm{\Σ}}_{q\∈A}f\left(q\right)}{{\mathrm{\Σ}}_{q\∈A}\mathrm{PN}}---\left(10\right)$

若p(τ)＜1，像素p即可判定为轮廓像素。实验证明，式（10）满足canny最优边缘3准则，边缘检测率100%，定位准确，保证单像素宽度连续。通过轮廓的判定式可以区分像素p是否为轮廓像素，从而得到建筑物的初始轮廓。

对如图5所示的点云分割图利用轮廓规则化分析方法分离树丛后，再进行上述建筑物轮廓的检测，得到如图6a、图6b和图6c所示的建筑物的初始轮廓。其中图6b为图6a中虚线框内的放大图，图6c是获得初始轮廓后的部分填充图。图6b中所标识的圆圈部分，可以看到采用式（10）可以对相邻建筑物的轮廓进行很好的区分，而现有技术中的其他方法在有可能无法对相邻建筑物的轮廓进行较好的区分。

本发明实施例中，在进行建筑物轮廓检测时，先确定检测轮廓的搜索方向是顺时针或逆时针，建筑物轮廓方向与搜索方向的夹角小于树丛轮廓方向与搜索方向的夹角。

为了确保轮廓围绕方向的唯一性，方便后续的操作，定义轮廓起点p(i,j)的搜索邻域为9邻域中的第i+1行，这样保证在图像坐标系逆时针方向，笛卡尔坐标系顺时针方向的搜索。因轮廓检测保证了单向素响应和连续，因此理论上轮廓检测可顺利从起点出发按逆时针方向围绕区域，最后回到起点。但是由于与建筑物连接的树丛的影响，导致可能出现在某一点出现多个搜索方向的情况。考虑到连接树丛出现在建筑物外部，可得建筑物轮廓方向与搜索方向夹角小于树丛轮廓方向与搜索方向夹角。

如图10所示，当前像素q出现3个可能搜索方向，此时需要计算pq与3个方向顺时针最小的夹角。以q原点建立坐标系(x',y′)，则各个方向数据偏移量(Δx,Δy)在新坐标系为：

$\left[\begin{array}{c}\mathrm{\Δ}{x}^{\′}\\ \mathrm{\Δ}{y}^{\′}\end{array}\right]=\left[\begin{array}{cc}\cos \mathrm{\θ}& \sin \mathrm{\θ}\\ -\sin \mathrm{\θ}& \cos \mathrm{\θ}\end{array}\right]\left[\begin{array}{c}\mathrm{\Δx}\\ \mathrm{\Δy}\end{array}\right]---\left(11\right)$

从而得到各个方向相对与pq的夹角：

${\mathrm{\Δ\θ}}^{\′}=\arctan \left(\frac{\mathrm{\Δ}{y}^{\′}}{\mathrm{\Δ}{x}^{\′}}\right)---\left(12\right)$

针对搜索方向在笛卡尔坐标系顺时针的特点，定义Δθ′取值范围为[-π,π]。由此，Δθ′最小的方向即为最可能建筑物轮廓方向，图10中θ₃为Δθ′最小的方向。这个选择建筑物轮廓方向的过程，也是消除部分连接树丛的过程。根据上述搜索方向的判定依据和方法可以有效避免树丛和噪声的影像，使得检测到的建筑物初始轮廓更精确。

如图3所示，本发明的一个实施例中，步骤S1还包括：

S15、对得到的建筑物初始轮廓进行直线段提取，得到优化后的建筑物初始轮廓。

LIDAR建筑物点云无法提供精确的轮廓线，因此检测得到的区域轮廓大多呈锯齿状，如图8a所示。那么在提取直线的时候，必然要求在角度或距离上有一定的容忍阈值。由此，可能导致直线段提取过度问题：一方面各种直线提取方法都通过加权考虑了像素偏离直线距离的影响，但若阈值范围内某个方向的偏移累积占了绝对优势，将使得直线逐渐偏离，那么原本阈值之外的像素会逐渐满足阈值要求，直线提取了过多像素；另一方面轮廓出现剧烈转折之后的几个像素也满足距离阈值或角度阈值，将会被加入当前直线段，而实际上它们是新直线段起点。这种错误对于采用最小二乘方法平差获得直线段参数的方法影响巨大，直线准确性和精度无法保证。

本发明实施例中，为了避免直线段提取过度问题，采用二步法对直线段进行提取，具体为：

a以预设的第一距离阈值搜索符合直线要求的直线片段；

b对直线片段采用预设的角度阈值和第二距离阈值合并邻接直线片段为直线段。

对图8a采用二步法进行优化后得到如图8b所示的建筑物初始轮廓。

在本发明的一个较佳实施例中，可用严格的第一距离阈值搜索最合乎直线要求的直线片段，第一距离阈值取像素偏移直线距离不大于1.5个像素；因阈值很小，可以避免提取剧烈转折之后的像素，也可避免逐渐偏离的积累，最终获得符合严格直线要求的直线片段；

对直线片段采用较大的角度阈值和距离阈值合并邻接直线片段为直线段，可取第二距离阈值为4个像素，角度阈值0.5弧度。因第一步a中集合了具有相同性质的像素，剧烈转折前后的像素也各自被集合到其所在直线片段，避免了单个像素无法逃离阈值的问题，使得错误提取的可能性大大降低了，最终提取的直线段具有更高的准确性。

步骤S2中，将LiDAR点云数据与航空遥感影像数据进行数据配准具体为：

首先提取建筑物角点特征，主要包括3个步骤：①检测LiDAR数据和航空影像边缘；②提取两条垂直边缘的交点为建筑物角点；③将角点分配到各个水平面。

利用同水平面3个LiDAR角点p_i、p_h和p_u及其航空影像中的匹配点q_j、q_k和q_v构成一个六元组(p_i,q_j,p_h,q_k,p_u,q_v)。令δ_ijhk(u,v)为当p_i与q_j配对且p_h与q_k配对时，由p_i、p_h和p_u组成的三角形与由q_j、q_k和q_v构成的三角形的相似度。定义Δp_ih为p_i和p_h间的距离，类似的定义其他距离，令：

${\mathrm{\α}}_{\mathrm{ijhk}}=\frac{\left|\right|\mathrm{\Δ}{p}_{\mathrm{ih}}\left|\right|}{\left|\right|\mathrm{\Δ}{p}_{\mathrm{jk}}\left|\right|},{\mathrm{\α}}_{\mathrm{ijuv}}=\frac{\left|\right|\mathrm{\Δ}{p}_{\mathrm{iu}}\left|\right|}{\left|\right|\mathrm{\Δ}{p}_{\mathrm{jv}}\left|\right|},{\mathrm{\α}}_{\mathrm{hkuv}}=\frac{\left|\right|\mathrm{\Δ}{p}_{\mathrm{hu}}\left|\right|}{\left|\right|\mathrm{\Δ}{p}_{\mathrm{kv}}\left|\right|}---\left(13\right)$

$\stackrel{\‾}{\mathrm{\α}}=({\mathrm{\α}}_{\mathrm{ijhk}}+{\mathrm{\α}}_{\mathrm{ijuv}}+{\mathrm{\α}}_{\mathrm{hkuv}})/3---\left(14\right)$

则：

${\mathrm{\δ}}_{\mathrm{ijhk}}(u,v)=\sqrt{{({\mathrm{\α}}_{\mathrm{ijhk}}-\stackrel{\‾}{\mathrm{\α}})}^2+{({\mathrm{\α}}_{\mathrm{ijuv}}-\stackrel{\‾}{\mathrm{\α}})}^2+{({\mathrm{\α}}_{\mathrm{hkuv}}-\stackrel{\‾}{\mathrm{\α}})}^2}---\left(15\right)$

若δ_ijhk(u,v)为0，则表示两个三角形完全相似，则点对(p_u,q_v)给予(p_i,q_j,p_h,q_k)以最大的支持度。随着δ_ijhk(u,v)的增大，其支持度减小，则可由相似度得到支持度的定义：

$\mathrm{\φ}\left({\mathrm{\δ}}_{\mathrm{ijhk}}(u,v)\right)=\frac{1}{1+{\left|{\mathrm{\δ}}_{\mathrm{ijhk}}(u,v)\right|}^2}---\left(16\right)$

当p_i与q_j配对且p_h与q_k配对时，p_u仅与和其联系的，对(p_i,q_j,p_h,q_k)的支持度最大的q_v相配对，即：

$m(u,v)=\underset{v\≠j,v\≠k}{\max }\mathrm{\φ}\left({\mathrm{\δ}}_{\mathrm{ijhk}}(u,v)\right)---\left(17\right)$

取所有p_u的平均值，得到(p_i,q_j,p_h,q_k)的初始匹配度量值为：

$S^0(p_i,q_j,p_h,q_k)=\frac{1}{m-2}\underset{u\≠i,u\≠h}{\mathrm{\Σ}}\underset{v\≠j,v\≠k}{\max }\mathrm{\φ}\left({\mathrm{\δ}}_{\mathrm{ijhk}}(u,v)\right)---\left(18\right)$

则(p_i,q_j)的初始匹配度量值为：

$S^0(p_i,q_j)=\underset{h,k}{\max }{S}^0(p_i,q_j,p_h,q_k)---\left(19\right)$

在第r次迭代时，(p_u，q_v)对p_i,q_j,p_h,q_k的匹配度量值同时依赖于p_u与q_v间的位置差别及他们的S^r-1(p_u,q_v)值，为了体现出这两个因素的相互作用，取他们的最小值，有：

$S^r(p_i,q_j,p_h,q_k)=\frac{1}{m-2}\underset{u\≠i,u\≠h}{\mathrm{\Σ}}\underset{v\≠j,v\≠k}{\max }\min \left[S^{r-1}(p_u,q_v)\mathrm{\φ}\left({\mathrm{\δ}}_{\mathrm{ijhk}}(u,v)\right)\right]---\left(20\right)$

以及：

$S^r(p_i,q_j)=\underset{h,k}{\max }{S}^r(p_i,q_j,p_h,q_k)---\left(21\right)$

为了避免由于阈值选取不当而带来的问题，可采用下面的决定迭代条件。

取 $d_r=\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}\underset{h\≠i}{\overset{m}{\underset{h=1}{\mathrm{\Σ}}}}\underset{k\≠j}{\overset{n}{\underset{k=1}{\mathrm{\Σ}}}}|S^r(p_i,q_j,p_h,q_k)-S^{r-1}(p_i,q_j,p_h,q_k)|,$ 当d_r＜ε，迭代终止。ε为预设门限。

利用松弛匹配结果可以得到LiDAR角点与航片粗略的变换关系，可以估算投影中心的坐标。设有粗略的投影中心O(X_s,Y_s,Z_s)，LiDAR角点A(x,y,z)，角点处的地面高为z₀，则角点A处的投影变形可估算为：

$\mathrm{\Δw}=\frac{z-z_0}{Z_s-z}\sqrt{{(X_s-x)}^2+{(Y_s-y)}^2}---\left(22\right)$

这里如果Zs的数值超出一般规范，则利用规范确定大概变化范围。

角点A相对于投影中心的方位估算为：

$\mathrm{\θ}=\arctan \frac{Y_s-y}{X_s-x}---\left(23\right)$

分别给定Δw和θ较大的变化范围，可得到每个LiDAR角点寻找匹配影像角点的局部搜索空间，此时不需要基于6元组的松弛法在图像全局寻找相似三角形。定义端点皆为角点的边缘为完整边缘，此时采取的匹配条件如下：

(1)LiDAR角点连接垂直边缘与影像角点连接垂直边缘的方位差近似等于松弛匹配解出的LiDAR数据与航空影像的方位差；

(2)LiDAR角点与影像角点各自至少连接一条完整边缘；

(3)完整边缘的比值与同水平面或邻近水平面的相似三角形比值近似相等；

利用搜索空间和匹配条件的双重约束，出现错误匹配的可能性很低，可通过平差过程检测和剔除之。

如有必要，还可根据匹配优化结果计算出的精确投影模型进一步缩小搜索空间，利用更宽松的匹配条件，搜索更多的匹配点对。

LiDAR数据与航空影像之间的严密的变换模型为描述中心投影的共线方程,利用单像空间后方交会求解外方位元素需要预知影像的内方位元素。本发明着力解决不具备内、外方位元素的航空影像与LiDAR数据配准，因此采用直接变换模型(DLT)表达投影变换关系，即：

$x+\frac{l_{1}X+l_{2}Y+l_{3}Z+l_4}{l_{9}X+l_{10}Y+l_{11}Z+1}=0$ $\left(24\right)$

$y+\frac{l_{5}X+l_{6}Y+l_{7}Z+l_8}{l_{9}X+l_{10}Y+l_{11}Z+1}=0$

其中，(x,y)为像方坐标值，(X,Y,Z)为物方坐标值，(l₁,l₂,…,l₁₁)为变换系数。解算变换系数(l₁,l₂,…,l₁₁)的关系式为：

Xl₁+Yl₂+Zl₃+l₄+0+0+0+0+xXl₉+xYl₁₀+xZl₁₁+x＝0

(25)

0+0+0+0+Xl₅+Yl₆+Zl₇+l₈+yXl₉+yYl₁₀+yZl₁₁+y＝0

当观测值大于6个控制点时，可以利用式(25)列出线性方程组，并略去一个方程式，可列出11个方程，平差求解变换系数(l₁,l₂,…,l₁₁)的近似值。

利用直接线性变换参数，可以求解遥感影像的内外方位元素，包括(X_s,Y_s,Z_s)摄影中心大地坐标，(x,y,-f)已知的像点的像空间坐标，(X,Y,Z)像点对应的待求大地坐标， $R=\left[\begin{array}{ccc}{a}_1& a_2& a_3\\ b_1& b_2& b_3\\ c_1& c_2& c_3\end{array}\right]$ 像片的旋转矩阵，是三个方向的旋转角的函数。数字影像与像平面尚需内定向配准，一般采用仿射变换模型即可，目的是确定像素行列(i,j)在像空间的坐标(x_i,y_j,-f)。影像定向完成之后，可以利用下式计算任何像素的空间坐标。

$\left\{\begin{array}{c}X=(Z-Z_s)\frac{a_{1}x+a_{2}y-a_{3}f}{c_{1}x+c_{2}y-c_{2}f}+X_s\\ Y=(Z-Z_s)\frac{b_{1}x+b_{2}y-b_{3}f}{c_{1}x+c_{2}y-c_{2}f}+Y_s\end{array}\right.---\left(26\right)$

本发明以如图4b所示的某地区的遥感影像图与如图4a所示的激光扫描数据进行数据配准，得到配准后的图片如图11a和图11b所示，其中两幅图中的十字标识为匹配点。图11a和图11b进行数据配准后，采用上述计算，可以得到如下表1的解算影像定向参数。

表1.解算影像定向参数

根据上表1中的影像定向参数和式（26）可计算航空遥感影像中任何像素的空间坐标，进而可以检测出航空遥感影像的影像边缘。

步骤S5在初始轮廓周围建立缓冲区，并在该缓冲区内搜索影像边缘，利用几何约束条件识别像方轮廓；该步骤具体包括：

（1）在LIDAR的初始轮廓附近建立缓冲区，搜索影像边缘；

（2）以几何约束识别像方轮廓。

第（2）步以几何约束识别像方轮廓具体为：设某点云轮廓直线段为L，在其缓冲区B(L)（如图9a所示）存在一个影像边缘集合E＝(e₁,e₂,…,e_n)（如图9b所示），E近似平行L，且min(length(E))＞T，此处T为减弱偶然误差影响设置的边缘长度门槛。

若n＝0；无法识别得到像方轮廓，以L替代之；

若n＝1；则认为识别得到唯一的像方轮廓；

若n＞1；因为较长的平行于点云轮廓的边缘，一般只可能因建筑物内部护墙等结构及轮廓本身产生，建筑物外部附近产生规则结构的可能性极低。检测判定像方轮廓的算法如下：

若p₁(i₁,j₁),p₂(i₂,j₂)为E中边缘e_i的端点，则按行方向扫描p₁,p₂与O的边界

的交点，得到各自的交点集合P_1c和P_2c：

P_1c＝{(I₁,J₁)|I₁＝i₁,J₁＝(j₁₁,j₁₂,…,j_1k)}

(27)

P_2c＝{(I₂,J₂)|I₂＝i₂,J₂＝(j₂₁,j₂₂,…,j_2k)}

因E近似平行L，且影像边缘大多断裂，短于LIDAR的初始轮廓，可得：若存在j_1i＞j₁，且mod((k-i+1)/2)＝0，或存在j_2i＞j₂，且mod((m-i+1)/2)＝0；则

即边缘e_i在建筑物的外侧。

对E进行上述判断，得到全部外侧边缘的集合Eo＝(eo₁,eo₂,…,eo_l),其中l≤n。若p_i1(i₁,j₁),p_i2(i₂,j₂)为边缘eo_i的端点，分别计算端点到L的距离d_i1，d_i2。考虑到E近似平行L，则eo_i到L的距离可由d_i＝(d_i1+d_i2)/2近似得到。对Eo得到集合D＝(d₁,d₂,…,d_i)，若d_j＝max(D)，则eo_j则判定为L对应的像方轮廓。此处，若l＝0，则无法搜索到像方轮廓，以L代替之。最终获得如图9c所示的像方轮廓。

接着是矢量化的工作，利用K均值聚类轮廓线中像素，计算线段描述参数和边界点；然后利用建筑物大多相邻边垂直，且存在大量平行轮廓的特点，建立单个建筑物整体的几何约束关系，通过整体平差纠正某些线段可能存在偏差，实现建筑物整体的规则化。

但因点云投影缺乏精确断裂线，建筑物轮廓与实际位置可能存在较多偏差，需要整体几何约束下的规则化过程，而像方轮廓也很难保证平行、垂直的规则，其本身也需要规则化的过程，因此规则化误差缺乏有效的定量评估手段。但可确定的是，规则轮廓与真实轮廓的偏差小于影像采样密度，由此可以估算规则化过程的精度损失范围。从点云投影提取初始轮廓的精度验证可利用基于点云的初始轮廓和轮廓投影相对于像方轮廓的偏移计算。基于此，可结合像方轮廓和基于点云的初始轮廓的偏移评估LiDAR提取轮廓误差。考虑到此处影像采样密度远低于LiDAR采样密度，本发明不考虑像方轮廓规则化过程的精度损失。

本发明实施例城市建筑物轮廓检测的系统主要用于实现上述城市建筑物轮廓检测的方法，如图12所示，该系统具体包括：

建筑物初始轮廓提取模块10，用于对LiDAR数据点云进行处理提取建筑物的初始轮廓；

配准模块20，用于将LiDAR点云数据与航空遥感影像数据进行数据配准；

影像定向模块30，用于根据数据配准的结果进行影像定向；

影像边缘检测模块40，用于根据影像定向，检测影像边缘；

像方轮廓识别模块50，用于在初始轮廓附近建立缓冲区，并在该缓冲区内搜索影像边缘，利用几何约束条件识别像方轮廓；

建筑物轮廓获得模块60，用于对像方轮廓进行轮廓矢量化与规则化，获得最终建筑物轮廓。

进一步地，本发明实施例中，如图13所示，建筑物初始轮廓提取模块10还包括：

点云数据分离模块102，用于采用线性预测法分离LiDAR点云数据；

建筑物点云分割模块104，用于根据分离结果并利用扫描带法分割建筑物点云；

树丛分离模块106，用于根据分割好的建筑物点云，利用轮廓规则化分析方法分离树丛；

建筑物轮廓检测模块108，用于根据分割好的建筑物点云以及分离的树丛，进行建筑物轮廓检测，得到建筑物的初始轮廓。

进一步地，本发明实施例中，建筑物轮廓检测模块108还用于在进行建筑物轮廓检测时，先确定检测轮廓的搜索方向，建筑物轮廓方向与搜索方向的夹角小于树丛轮廓方向与搜索方向的夹角。

进一步地，本发明实施例中，如图14所示，建筑物初始轮廓提取模块10上文实施例的基础上还包括：

初始轮廓优化模块110，用于对建筑物轮廓检测模块得到的建筑物初始轮廓进行直线段提取，得到优化后的建筑物初始轮廓。

本发明通过利用LiDAR数据易于自动化检测目标的特点，获得建筑物初始轮廓，再在初始轮廓缓四周建立缓冲区内，并在缓冲区内搜索符合实际的影像边缘，利用几何约束条件识别像方轮廓；对像方轮廓进行规则化最终获得精确、完整的建筑物轮廓。本发明可提高建筑物轮廓检测的自动化程度，可以实现大范围快速检测；且提取的建筑物轮廓具有最好的质量，避免了噪声和影像质量导致的算法局部极值和迭代死循环问题，可有效的保证算法的稳健性；另外，本发明对于LiDAR数据和遥感影像的来源、格式无特殊要求，具有很大的普适性。

应当理解的是，对本领域普通技术人员来说，可以根据上述说明加以改进或变换，而所有这些改进和变换都应属于本发明所附权利要求的保护范围。

Claims (10)

1.一种城市建筑物轮廓检测的方法，其特征在于，包括以下步骤：

S1、对LiDAR点云数据进行处理提取建筑物的初始轮廓；

S2、将所述LiDAR点云数据与航空遥感影像数据进行数据配准；

S3、根据数据配准的结果进行影像定向；

S4、根据所述影像定向，检测影像边缘；

S5、在所述初始轮廓周围建立缓冲区，并在该缓冲区内搜索所述影像边缘，利用几何约束条件识别像方轮廓；

S6、对所述像方轮廓进行轮廓矢量化与规则化，获得最终建筑物轮廓。

2.根据权利要求1所述的城市建筑物轮廓检测的方法，其特征在于，步骤S1包括：

S11、采用线性预测法分离所述LiDAR点云数据；

S12、根据分离结果并利用扫描带法分割建筑物点云；

S13、根据分割好的建筑物点云，利用轮廓规则化分析方法分离树丛；

S14、根据分割好的建筑物点云以及分离的树丛，进行建筑物轮廓检测，得到建筑物的初始轮廓。

3.根据权利要求2所述的城市建筑物轮廓检测的方法，其特征在于，在进行建筑物轮廓检测时，先确定检测轮廓的搜索方向，建筑物轮廓方向与所述搜索方向的夹角小于树丛轮廓方向与所述搜索方向的夹角。

4.根据权利要求3所述的城市建筑物轮廓检测的方法，其特征在于，步骤S1还包括：

S15、对得到的建筑物初始轮廓进行直线段提取，得到优化后的建筑物初始轮廓。

5.根据权利要求4所述的城市建筑物轮廓检测的方法，其特征在于，所述直线段提取具体为：

以预设的第一距离阈值搜索符合直线要求的直线片段；

对直线片段采用预设的角度阈值和第二距离阈值合并邻接直线片段为直线段。

6.根据权利要求1所述的城市建筑物轮廓检测的方法，其特征在于，所述第一距离阈值为小于等于1.5个像素；所述第二距离阈值为4个像素；所述角度阈值为0.5弧度。

7.一种城市建筑物轮廓检测的系统，其特征在于，包括：

建筑物初始轮廓提取模块，用于对LiDAR数据点云进行处理提取建筑物的初始轮廓；

配准模块，用于将所述LiDAR点云数据与航空遥感影像数据进行数据配准；

影像定向模块，用于根据数据配准的结果进行影像定向；

影像边缘检测模块，用于根据所述影像定向，检测影像边缘；

像方轮廓识别模块，用于在所述初始轮廓附近建立缓冲区，并在该缓冲区内搜索所述影像边缘，利用几何约束条件识别像方轮廓；

建筑物轮廓获得模块，用于对所述像方轮廓进行轮廓矢量化与规则化，获得最终建筑物轮廓。

8.根据权利要求7所述的城市建筑物轮廓检测的系统，其特征在于，所述建筑物初始轮廓提取模块还包括：

点云数据分离模块，用于采用线性预测法分离所述LiDAR点云数据；

建筑物点云分割模块，用于根据分离结果并利用扫描带法分割建筑物点云；

树丛分离模块，用于根据分割好的建筑物点云，利用轮廓规则化分析方法分离树丛；

建筑物轮廓检测模块，用于根据分割好的建筑物点云以及分离的树丛，进行建筑物轮廓检测，得到建筑物的初始轮廓。

9.根据权利要求8所述的城市建筑物轮廓检测的系统，其特征在于,所述建筑物初始轮廓提取模块还包括：

初始轮廓优化模块，用于对所述建筑物轮廓检测模块得到的建筑物初始轮廓进行直线段提取，得到优化后的建筑物初始轮廓。

10.根据权利要求8所述的城市建筑物轮廓检测的系统，其特征在于, 所述建筑物轮廓检测模块还用于在进行建筑物轮廓检测时，先确定检测轮廓的搜索方向，建筑物轮廓方向与所述搜索方向的夹角小于树丛轮廓方向与所述搜索方向的夹角。

Images