CN104036550A - 基于形状语义的建筑立面激光雷达点云解译与重建的方法 - Google Patents

Abstract

建筑立面建模是城市建模的重要组成部分。本发明基于形状语义和增强学习优化策略，提出了一种建筑立面激光点云解译和重建的方法。首先基于数据驱动对点云进行建筑学语义上的预分类，接着设计了可以灵活应用于各种建筑立面风格的二分形状语义，最后将建筑立面的最佳形状语义求解的问题建模为马尔科夫决策过程，采用增强学习算法结合支持向量机预分类的结果优化求解，实现对建筑立面的解译和三维重建。本发明提出的方法不仅能够较为精确的对建筑立面进行精细建模，而且能够适应于各种风格建筑立面的解译，对于建筑立面点云缺失的情况也有很好的鲁棒性。

Description

基于形状语义的建筑立面激光雷达点云解译与重建的方法

一、技术领域

本发明涉及基于形状语义的建筑立面激光点云解译与重建的方法，属于空间信息技术领域。

二、背景技术

建筑立面的重建是数字城市建模的重要组成部分。传统的建筑立面建模方法主要是从二维的图像中来提取建筑的精细结构。近年来，随着激光雷达技术的发展，可以通过车载激光雷达或是机载激光雷达对地面进行快速的扫描从而获得城市模型的点云。与传统的图片相比，激光雷达点云是建筑模型的直接三维表示，包含了恢复建筑结构特征的几乎所有特征。一些研究显示了激光雷达数据可以高效地应用于城市的重建中。

本申请专利探索如何从激光雷达点云中提取建筑精细的结构特征，从而为建筑的自动化三维重建提供基础。数据驱动的方法利用建筑图像或激光雷达点云自身的信息，构造相应的特征属性自底向上对建筑进行分类和识别。数据驱动的方法虽然很精确，但是却无法很好显示建筑的语义特征(比如，建筑由门、窗、阳台、屋顶等建筑学意义上的物体组成)。为此人们提出了形状语义的概念用来自顶向下对建筑进行建筑学语义上的建模。采用形状语义的难点在于定义适合建筑的语义以及如何高效地求解最优语义。本申请专利将研究如何设计形状语义对建筑的点云立面进行解译并进行高效地求解，以实现对不规则且存在点云缺失的立面进行建模。

三、发明内容

1、目的：建筑立面的重建是数字城市建模的重要组成部分。传统的建模方法主要是基于人工或是半自动的重建，由于城市里的建筑数量庞大，且建筑形状复杂，人工建立城市模型需要建模者要有一定的经验，并耗费大量时间。目前存在一些半自动生成建筑模型的方法，但效果都不令人满意。与传统的基于图像的建模方法相比，激光雷达点云是建筑模型的直接三维表示。本申请专利实现了从激光雷达点云中自动提取建筑立面精细的结构特征，从而为城市建筑的自动化三维重建提供基础。

2、技术方案

基于形状语义的建筑立面激光点云解译与重建的方法，其特征在于，包括如下步骤(如图1)：

步骤一：点云噪声去除

点云中噪声的存在对于后面计算曲度会有较大的影响，因此，采用统计方法去除点云噪声。对于输入点云的每一点，计算它到邻近点的距离的均值μ和标准差σ。那些距离落在μ±α·σ之外的点则被视为噪声而移除，其中α是和选取的邻近点的数目相近的一个值。本申请专利中，α＝1.0，邻近点的数目为20。

步骤二：坐标系转换

为了方便后续的处理，把建筑立面所在的平面转换到X-Y平面中去，即，建筑的水平方向与X轴平行，竖直方向与Y轴平行。这样一来建筑立面点云的z坐标就代表了建筑结构的深度信息，并且建筑的水平和竖直方向与坐标轴平行，以便于下文对建筑立面进行网格剖分。

建筑立面的坐标系变换陈述为如下问题，给定建筑点云当前所处的基M₁＝{α₁，α₂，α₃}以及目标坐标系的基M₂＝{β₁，β₂，β₃}，求M₁到M₂的过渡矩阵。其中，β1＝{1，0，0}，β2＝{0，1，0}，β₃＝{0，0，1}。而要得到两组基之间的过渡矩阵P，只需要给出目标坐标系的基中的每个向量(β₁，β₂，β₃)在当前所处基下的坐标(p₁，p₂，p₃)，依次以这些坐标为列向量组成矩阵即是从M₁到M₂的过渡矩阵P。

目前的任务是计算(β₁，β₂，β₃)在当前所处基下的坐标(p₁，p₂，p₃)。由前述知，建筑的水平方向对其到X轴，建筑的竖直方向对齐到Y轴，因此，基向量β₁在当前基下对应的坐标p₁是建筑立面水平方向的向量，基向量β₂对应的坐标p₂则是建筑立面竖直方向的向量，基向量p₃对应的坐标则是建筑平面的法向量。

采用随机抽样一致模型计算建筑立面的平面参数，继而得到其归一化的法向量，即为p₃。对于p₂，使用建筑墙面边缘线所处的向量表示，而p₁则可以由p₂和p₃的外积求得。得到过渡矩阵P后，对原始点云中的点X应用变换P即可得到在新坐标系下的坐标X＇。

步骤三：点云特征提取

为了有效识别建筑立面的结构(如门、窗户、阳台)，设计了基于网格的点云特征提取方法。首先对建筑立面进行剖分，使得建筑立面由一系列三维的小网格组成。在本申请专利中，网格的大小设置为0.1米。这样对点云分类的任务就变成了对这些网格的分类。本申请专利设计了三种网格的属性，用来通过支持向量机对建筑立面点云中的结构单元分类。

(1)曲度

曲度反映了建筑局部结构的变化情况，对于提取建筑结构的变化有很大作用。定义网格内的曲度特征是网格内所有点的曲度的平均值。计算每点的曲度近似为计算该点和邻近的k个点近似组成的平面的曲度，二者可以通过分析邻近点的协相关矩阵的特征向量和特征矩阵得到。

对于一点p_i，它的协相关矩阵为：

$C=\frac{1}{k}\underset{i=1}{\overset{k}{\mathrm{\Σ}}}(p_i-\stackrel{\‾}{p})\·{(p_i-\stackrel{\‾}{p})}^T,C\·\stackrel{\→}{v_j}={\mathrm{\λ}}_j\·\stackrel{\→}{v_j},j\∈\left\{\mathrm{0,1,2}\right\}---\left(1\right)$

其中，k是p_i的邻近点的数目，代表邻近点的三维中心坐标，λ_j是协相关矩阵的第j个特征值，是第j个特征向量。

点p_i的曲度为：

γ_pi＝λ₀/(λ₀+λ₁+λ₂) (2)

假设网格G中有n个点，那么该网格的曲度则表示为：

${\mathrm{\γ}}_G=\frac{1}{n}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\γ}}_{\mathrm{pi}}---\left(3\right)$

(2)密度

除了曲度外，网格内点的密度变化同样能反映不同的建筑结构，另外，对于密度还能够明确地反映数据缺失或是背景网格。

设网格G中有n个点，点数最多网格内的点数为N，则该网格的密度d_G表示为：

d_G＝n/N (4)

(3)深度

由于建筑表面各个结构通常位于不同的平面，通过深度值可以较容易识别出不同的建筑立面结构。步骤二已经将建筑平面转换到了X-Y平面，因此点云的z坐标就表示了深度的信息。

假设网格G中有n个点，那么该网格的深度表示为：

$s_G=\frac{1}{n}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{z}_i---\left(5\right)$

其中z_i为点i的z坐标。

以上三种特征将作为网格的特征向量用于步骤四的支持向量机分类中。

步骤四：基于支持向量机的建筑立面点云分类

由于建筑立面的建筑学意义上的结构(如窗户、阳台、门、屋顶等)是从人类认知角度去定义的，采用非监督的方法往往很难得到满意的结果。为此，本申请专利采用了支持向量机的方法进行分类。

首先选择训练集，训练集里的每一个样本都有一个类别标签及多个属性值。用数学的语言来描述就是：给定一个训练集(x₁，y₁)，i＝1，…，l，x_i∈Rⁿ表示样本的特征向量，y∈{-1，1}^l表示样本的类别向量，支持向量机求解以下正定规划的最优解：

训练向量的特征分别是曲度、密度和深度。本申请专利选择了半径基函数作为核函数，并设定γ＝1。在步骤三中定义的三种网格属性分别归化到[-1，1]之间，对样本集进行训练，最后进行预测。

在分类完成之后，便可以得到了某个网格属于特定建筑结构的概率：

步骤五：建筑立面的二分形状语义表达

从给定的建筑立面点云中得到表达墙面结构的最佳剖分树，就可以重建建筑立面。用数学语言描述如下：

给定建筑立面点云O，二分形状语义G＝(N，T，R，ω)，其中N是非终端形状语义的集合，如果一个带标号的矩形c(x，y，w，h)出现在一个规则的左边，那么它就是非终端形状语义；T是终端语义的集合，，(x，y，w，h)定义了矩形在坐标系中的位置、宽度w和高度h；如果一个带标号的矩形c(x，y，w，h)不能出现在一个规则的左边，那么它就是终端形状语义；ω是一个特殊的非终端形状语义，称为初始形状，剖分总是从初始形状开始进行。R是一个有限的剖分规则的集合。设L(G)为二分形状语义的所有可能的剖分树(剖分方式)。寻找到一个最佳的剖分语言s∈L(G)，使得如下的目标函数最大：

$\underset{s}{\arg }\underset{s\∈L\left(G\right)}{\max }\underset{x,y\∈\mathrm{facade}}{\mathrm{\Σ}}p\left(s(x,y)\right|x,y)---\left(8\right)$

其中，s(x，y)表示剖分s得到的网格(x，y)的建筑结构类别；p(s(x，y)|x，y)即为步骤四得到的网格(x，y)处建筑类别为s(x，y)的概率。从上式(8)可以看出，评价一个剖分的好坏在于它剖分得到的终端形状的类别与预分类的结果是否一致。本申请专利采用了逐网格进行累加的方式对建筑立面上的所有网格进行判断，看剖分的结果与预分类的结果是否一致，若一致则进行累加1，否则累加0，最终的累加值则反映了该剖分的好坏。最大化式(15)便可以得到建筑立面的最佳剖分。

步骤六：城市建筑立面三维建模

二分形状语义可以有效表达建筑立面结构的生成过程，只要构建建筑立面的二分形状语义剖分树后，就可以对该建筑立面的解译。从二分语义的剖分过程来看，可以把求解(8)式看成是一个马尔科夫决策过程(MDP)。

(1)最优形状语义的求解

二分形状语义产生建筑的过程用MDP来建模，下面是MDP里的各个要素在建筑剖分中代表的含义。

智能体：MDP中的智能体表示的是形状语义剖分建筑立面的过程。在剖分建筑立面的时候，每一步都要处理一个非终端形状，智能体则决定采用哪种规则及规则参数应用于该形状。

环境：马尔可夫过程中的环境与智能体交互，指的是随着剖分进行不断建立的剖分树。随着剖分的不断进行，剖分树不断完善，环境也就有了对模型更多的知识。

状态：为剖分树中的结点。该结点表示了当前剖分所对应的形状及其位置，用state(x，y，symbol)表示，其中(x，y)为当前状态所处的位置，symbol表示当前形状的名字，symbol只能是非终端形状。

行为：为智能体可以选择的语法规则，用action(rule，parameter)表示，其中rule表示语法，parameter表示剖分的宽度。

回报值：设马尔可夫决策在t时刻执行动作action(rule，parameter)得到的终端形状为A(x，y，w，h，a)，那么相应的回报值为，

$R_t=\underset{i=x}{\overset{x+w}{\mathrm{\Σ}}}\underset{j=y}{\overset{y+h}{\mathrm{\Σ}}}p\left(a\right|x,y)---\left(9\right)$

其中p(a|x，y)表示建筑立面(x，y)处是a类别的概率。

设MDP决策过程的序列为T，那么长期的回报为：

$R=\underset{t=0}{\overset{T}{\mathrm{\Σ}}}{R}_t---\left(10\right)$

由于MDP决策完毕后，把建筑立面剖分完毕，因此有，

$R=\underset{t=0}{\overset{T}{\mathrm{\Σ}}}{R}_t=\underset{x,y\∈I}{\mathrm{\Σ}}p(x,y,s(x,y))---\left(11\right)$

MDP的最优解即是获取的最优形状语义。

采用增强学习算法求解以上叙述的MDP问题。增强学习采用的是Q-Learning算法。

(2)建筑立面深度信息的恢复

强化学习得到的最佳语义可以表示为剖分树的形式，其中剖分树的叶子节点是组成立面的所有建筑结构(终端形状)，记为Leaf_i(c，x，y，w，h)，c表示叶子节点i所属的建筑类别，(x，y，w，h)定义了建筑结构的外包矩形R。通过这些叶子节点，计算各类建筑结构在建筑立面上的分布情况，并构造出建筑该建筑立面的二维平面。还需要得到各种建筑结构的深度以完成最后的三维建模，为此，查询点云中R内的所有点，利用随机抽样一致模型得到这些点所的平面，设该平面为：

a_ix+b_iy+c_iz+d_i＝0 (12)

在上述坐标系转换中，已经求得建筑主平面所在的方程：

Ax+By+Cx+D＝0 (13)

由于建筑结构的平面与主平面平行，则建筑类别c离主平面的深度d为：

$d=\left|\frac{d_i-D}{\sqrt{A^2+B^2+C^2}}\right|---\left(14\right)$

如此便恢复了建筑结构c的深度信息。

3、优点及功效：本申请专利基于形状语义和增强学习优化策略，提出了一种建筑立面激光雷达点云解译和重建的方法。首先基于数据驱动对点云进行建筑学语义上的预分类，接着设计了灵活应用于各种建筑立面风格的二分形状语义，最后将建筑立面的最佳形状语义求解的问题建模为马尔科夫决策过程，采用增强学习算法结合支持向量机预分类的结果优化求解，实现对建筑立面的解译和三维重建。本申请专利提出的方法不仅能够较为精确的对建筑立面进行精细建模，而且能够适应于各种风格建筑立面的解译，对于建筑立面点云缺失的情况也有很好的鲁棒性。

四、附图说明

图1本发明方法的流程示意图。

图2一个立面建模过程。(a)原始点云；(b)形状语义解译结果；(c)恢复深度建模结果。

图3建筑建模立面建模过程。(a)原始点云；(b)形状语义解译结果；(c)恢复深度建模结果。

图4北师大科技楼立面建模过程。(a)原始点云；(b)形状语义解译结果；(c)恢复深度建模结果。

图5点云严重缺失的立面建模过程。(a)原始点云；(b)形状语义解译结果；(c)恢复深度建模结果。

五、具体实施方式

基于形状语义的建筑立面激光点云解译与重建的方法，其特征在于，包括如下步骤(如图1)：

步骤一：点云噪声去除

点云中噪声的存在对于后面计算曲度会有较大的影响，因此，采用统计方法去除点云噪声。对于输入点云的每一点，计算它到邻近点的距离的均值μ和标准差σ。那些距离落在μ±α·σ之外的点则被视为噪声而移除，其中α是和选取的邻近点的数目相近的一个值。本申请专利中，α＝1.0，邻近点的数目为20。

步骤二：坐标系转换

为了方便后续的处理，把建筑立面所在的平面转换到X-Y平面中去，即，建筑的水平方向与X轴平行，竖直方向与Y轴平行。这样一来建筑立面点云的z坐标就代表了建筑结构的深度信息，并且建筑的水平和竖直方向与坐标轴平行，以便于下文对建筑立面进行网格剖分。

建筑立面的坐标系变换陈述为如下问题，给定建筑点云当前所处的基M₁＝{α₁，α₂，α₃}以及目标坐标系的基M₂＝{β₁，β₂，β₃}，求M₁到M₂的过渡矩阵。其中，β1＝{1，0，0}，β2＝{0，1，0}，β₃＝{0，0，1}。而要得到两组基之间的过渡矩阵P，只需要给出目标坐标系的基中的每个向量(β₁，β₂，β₃)在当前所处基下的坐标(p₁，p₂，p₃)，依次以这些坐标为列向量组成矩阵即是从M₁到M₂的过渡矩阵P。

目前的任务是计算(β₁，β₂，β₃)在当前所处基下的坐标(p₁，p₂，p₃)。由前述知，建筑的水平方向对其到X轴，建筑的竖直方向对齐到Y轴，因此，基向量β₁在当前基下对应的坐标p₁是建筑立面水平方向的向量，基向量β₂对应的坐标p₂则是建筑立面竖直方向的向量，基向量p₃对应的坐标则是建筑平面的法向量。

采用随机抽样一致模型计算建筑立面的平面参数，继而得到其归一化的法向量，即为p₃。对于p₂，使用建筑墙面边缘线所处的向量表示，而p₁则可以由p₂和p₃的外积求得。得到过渡矩阵P后，对原始点云中的点X应用变换P即可得到在新坐标系下的坐标X＇。

步骤三：点云特征提取

为了有效识别建筑立面的结构(如门、窗户、阳台)，设计了基于网格的点云特征提取方法。首先对建筑立面进行剖分，使得建筑立面由一系列三维的小网格组成。在本申请专利中，网格的大小设置为0.1米。这样对点云分类的任务就变成了对这些网格的分类。本申请专利设计了三种网格的属性，用来通过支持向量机对建筑立面点云中的结构单元分类。

(1)曲度

曲度反映了建筑局部结构的变化情况，对于提取建筑结构的变化有很大作用。定义网格内的曲度特征是网格内所有点的曲度的平均值。计算每点的曲度近似为计算该点和邻近的k个点近似组成的平面的曲度，二者可以通过分析邻近点的协相关矩阵的特征向量和特征矩阵得到。

对于一点p_i，它的协相关矩阵为：

$C=\frac{1}{k}\underset{i=1}{\overset{k}{\mathrm{\Σ}}}(p_i-\stackrel{\‾}{p})\·{(p_i-\stackrel{\‾}{p})}^T,C\·\stackrel{\→}{v_j}={\mathrm{\λ}}_j\·\stackrel{\→}{v_j},j\∈\left\{\mathrm{0,1,2}\right\}---\left(1\right)$

其中，k是p_i的邻近点的数目，代表邻近点的三维中心坐标，λ_j是协相关矩阵的第j个特征值，是第j个特征向量。

点p_i的曲度为：

γ_pi＝λ₀/(λ₀+λ₁+λ₂) (2)

假设网格G中有n个点，那么该网格的曲度则表示为：

${\mathrm{\γ}}_G=\frac{1}{n}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\γ}}_{\mathrm{pi}}---\left(3\right)$

(2)密度

除了曲度外，网格内点的密度变化同样能反映不同的建筑结构，另外，对于密度还能够明确地反映数据缺失或是背景网格。

设网格G中有n个点，点数最多网格内的点数为N，则该网格的密度d_G表示为：

d_G＝n/N (4)

(3)深度

由于建筑表面各个结构通常位于不同的平面，通过深度值可以较容易识别出不同的建筑立面结构。步骤二已经将建筑平面转换到了X-Y平面，因此点云的z坐标就表示了深度的信息。

假设网格G中有n个点，那么该网格的深度表示为：

$s_G=\frac{1}{n}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{z}_i---\left(5\right)$

其中z_i为点i的z坐标。

以上三种特征将作为网格的特征向量用于步骤四的支持向量机分类中。

步骤四：基于支持向量机的建筑立面点云分类

由于建筑立面的建筑学意义上的结构(如窗户、阳台、门、屋顶等)是从人类认知角度去定义的，采用非监督的方法往往很难得到满意的结果。为此，本申请专利采用了支持向量机的方法进行分类。

首先选择训练集，训练集里的每一个样本都有一个类别标签及多个属性值。用数学的语言来描述就是：给定一个训练集(x_i，y_i)，i＝1，…，l，x_i∈Rⁿ表示样本的特征向量，y∈{-1，1}^l表示样本的类别向量，支持向量机求解以下正定规划的最优解：

训练向量的特征分别是曲度、密度和深度。本申请专利选择了半径基函数作为核函数，并设定γ＝1。在步骤三中定义的三种网格属性分别归化到[-1，1]之间，对样本集进行训练，最后进行预测。

在分类完成之后，便可以得到了某个网格属于特定建筑结构的概率：

步骤五：建筑立面的二分形状语义表达

从给定的建筑立面点云中得到表达墙面结构的最佳剖分树，就可以重建建筑立面。用数学语言描述如下：

给定建筑立面点云O，二分形状语义G＝(N，T，R，ω)，其中N是非终端形状语义的集合，如果一个带标号的矩形c(x，y，w，h)出现在一个规则的左边，那么它就是非终端形状语义；T是终端语义的集合，，(x，y，w，h)定义了矩形在坐标系中的位置、宽度w和高度h；如果一个带标号的矩形c(x，y，w，h)不能出现在一个规则的左边，那么它就是终端形状语义；ω是一个特殊的非终端形状语义，称为初始形状，剖分总是从初始形状开始进行。R是一个有限的剖分规则的集合。设L(G)为二分形状语义的所有可能的剖分树(剖分方式)。寻找到一个最佳的剖分语言s∈L(G)，使得如下的目标函数最大：

$\underset{s}{\arg }\underset{s\∈L\left(G\right)}{\max }\underset{x,y\∈\mathrm{facade}}{\mathrm{\Σ}}p\left(s(x,y)\right|x,y)---\left(8\right)$

其中，s(x，y)表示剖分s得到的网格(x，y)的建筑结构类别；p(s(x，y)|x，y)即为步骤四得到的网格(x，y)处建筑类别为s(x，y)的概率。从上式(8)可以看出，评价一个剖分的好坏在于它剖分得到的终端形状的类别与预分类的结果是否一致。本申请专利采用了逐网格进行累加的方式对建筑立面上的所有网格进行判断，看剖分的结果与预分类的结果是否一致，若一致则进行累加1，否则累加0，最终的累加值则反映了该剖分的好坏。最大化式(15)便可以得到建筑立面的最佳剖分。

步骤六：城市建筑立面三维建模

二分形状语义可以有效表达建筑立面结构的生成过程，只要构建建筑立面的二分形状语义剖分树后，就可以对该建筑立面的解译。从二分语义的剖分过程来看，可以把求解(8)式看成是一个马尔科夫决策过程(MDP)。

(1)最优形状语义的求解

二分形状语义产生建筑的过程用MDP来建模，下面是MDP里的各个要素在建筑剖分中代表的含义。

智能体：MDP中的智能体表示的是形状语义剖分建筑立面的过程。在剖分建筑立面的时候，每一步都要处理一个非终端形状，智能体则决定采用哪种规则及规则参数应用于该形状。

环境：马尔可夫过程中的环境与智能体交互，指的是随着剖分进行不断建立的剖分树。随着剖分的不断进行，剖分树不断完善，环境也就有了对模型更多的知识。

状态：为剖分树中的结点。该结点表示了当前剖分所对应的形状及其位置，用state(x，y，symbol)表示，其中(x，y)为当前状态所处的位置，symbol表示当前形状的名字，symbol只能是非终端形状。

行为：为智能体可以选择的语法规则，用action(rule，parameter)表示，其中rule表示语法，parameter表示剖分的宽度。

回报值：设马尔可夫决策在t时刻执行动作action(rule，parameter)得到的终端形状为A(x，y，w，h，a)，那么相应的回报值为，

$R_t=\underset{i=x}{\overset{x+w}{\mathrm{\Σ}}}\underset{j=y}{\overset{y+h}{\mathrm{\Σ}}}p\left(a\right|x,y)---\left(9\right)$

其中p(a|x，y)表示建筑立面(x，y)处是a类别的概率。

设MDP决策过程的序列为T，那么长期的回报为：

$R=\underset{t=0}{\overset{T}{\mathrm{\Σ}}}{R}_t---\left(10\right)$

由于MDP决策完毕后，把建筑立面剖分完毕，因此有，

$R=\underset{t=0}{\overset{T}{\mathrm{\Σ}}}{R}_t=\underset{x,y\∈I}{\mathrm{\Σ}}p(x,y,s(x,y))---\left(11\right)$

MDP的最优解即是获取的最优形状语义。

采用增强学习算法求解以上叙述的MDP问题。增强学习采用的是Q-Learning算法。

(2)建筑立面深度信息的恢复

强化学习得到的最佳语义可以表示为剖分树的形式，其中剖分树的叶子节点是组成立面的所有建筑结构(终端形状)，记为Leaf_i(c，x，y，w，h)，c表示叶子节点i所属的建筑类别，(x，y，w，h)定义了建筑结构的外包矩形R。通过这些叶子节点，计算各类建筑结构在建筑立面上的分布情况，并构造出建筑该建筑立面的二维平面。还需要得到各种建筑结构的深度以完成最后的三维建模，为此，查询点云中R内的所有点，利用随机抽样一致模型得到这些点所的平面，设该平面为：

a_ix+b_iy+c_iz+d_i＝0 (12)

在上述坐标系转换中，已经求得建筑主平面所在的方程：

Ax+By+Cx+D＝0 (13)

由于建筑结构的平面与主平面平行，则建筑类别c离主平面的深度d为：

$d=\left|\frac{d_i-D}{\sqrt{A^2+B^2+C^2}}\right|---\left(14\right)$

如此便恢复了建筑结构c的深度信息。

实施例1：

图2显示了一个建筑结构数目较少，但是层间是不对称的立面，下面一层的窗户分布明显要不同于其它的楼层。本申请专利方法建模得到的结果很好地体现了原始建筑的结构特征。图3-4显示了对更为复杂的墙面进行建模的结果，可以看出，对于墙面结构众多的情况，算法也具有很好的表现。图5显示了在点云严重缺失的情况下，本申请专利提出的算法也能很好地对其建模。

Claims (1)

1.基于形状语义的建筑立面激光点云解译与重建的方法，其特征在于，包括如下步骤：

步骤一：点云噪声去除

对于输入点云的每一点，计算它到邻近点的距离的均值μ和标准差σ，那些距离落在μ±α·σ之外的点则被视为噪声而移除，其中α是和选取的邻近点的数目相近的一个值，本申请专利中，α＝1.0，邻近点的数目为20；

步骤二：坐标系转换

把建筑立面所在的平面转换到X-Y平面中去，即，建筑的水平方向与X轴平行，竖直方向与Y轴平行；

建筑立面的坐标系变换陈述为如下问题，给定建筑点云当前所处的基M₁＝{α₁，α₂，α₃}以及目标坐标系的基M₂＝{β₁，β₂，β₃}，求M₁到M₂的过渡矩阵，其中，β1＝{1，0，0}，β2＝{0，1，0}，β₃＝{0，0，1}，而要得到两组基之间的过渡矩阵P，只需要给出目标坐标系的基中的每个向量(β₁，β₂，β₃)在当前所处基下的坐标(p₁，p₂，p₃)，依次以这些坐标为列向量组成矩阵即是从M₁到M₂的过渡矩阵P；

基向量β₁在当前基下对应的坐标p₁是建筑立面水平方向的向量，基向量β₂对应的坐标p₂则是建筑立面竖直方向的向量，基向量p₃对应的坐标则是建筑平面的法向量；

采用随机抽样一致模型计算建筑立面的平面参数，继而得到其归一化的法向量，即为p₃，对于p₂，使用建筑墙面边缘线所处的向量表示，而p₁则可以由p₂和p₃的外积求得，得到过渡矩阵P后，对原始点云中的点X应用变换P即可得到在新坐标系下的坐标X＇；

步骤三：点云特征提取

首先对建筑立面进行剖分，使得建筑立面由一系列三维的小网格组成，在本申请专利中，网格的大小设置为0.1米，设计了下面三种网格属性，用来通过支持向量机对建筑立面点云中的结构单元分类；

(1)曲度

定义网格内的曲度特征是网格内所有点的曲度的平均值，计算每点的曲度近似为计算该点和邻近的k个点近似组成的平面的曲度，二者可以通过分析邻近点的协相关矩阵的特征向量和特征矩阵得到；

对于一点p_i，它的协相关矩阵为：

$C=\frac{1}{k}\underset{i=1}{\overset{k}{\mathrm{\Σ}}}(p_i-\stackrel{\‾}{p})\·{(p_i-\stackrel{\‾}{p})}^T,C\·\stackrel{\→}{v_j}={\mathrm{\λ}}_j\·\stackrel{\→}{v_j},j\∈\left\{\mathrm{0,1,2}\right\}---\left(1\right)$

其中，k是p_i的邻近点的数目，代表邻近点的三维中心坐标，λ_j是协相关矩阵的第j个特征值，是第j个特征向量；

点p_i的曲度为：

γ_pi＝λ₀/(λ₀+λ₁+λ₂) (2)

假设网格G中有n个点，那么该网格的曲度则表示为：

${\mathrm{\γ}}_G=\frac{1}{n}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\γ}}_{\mathrm{pi}}---\left(3\right)$

(2)密度

设网格G中有n个点，点数最多网格内的点数为N，则该网格的密度d_G表示为：

d_G＝n/N (4)

(3)深度

假设网格G中有n个点，那么该网格的深度表示为：

$s_G=\frac{1}{n}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{z}_i---\left(5\right)$

其中z_i为点i的z坐标；

以上三种特征将作为网格的特征向量用于步骤四的支持向量机分类中；

步骤四：基于支持向量机的建筑立面点云分类

首先选择训练集，训练集里的每一个样本都有一个类别标签及多个属性值，用数学的语言来描述就是：给定一个训练集(x_i，y_i)，i＝1，…，l，x_i∈Rⁿ表示样本的特征向量，y∈{-1，1}^l表示样本的类别向量，支持向量机求解以下正定规划的最优解：

训练向量的特征分别是曲度、密度和深度，本申请专利选择了半径基函数作为核函数，并设定γ＝1，在步骤三中定义的三种网格属性分别归化到[-1，1]之间，对样本集进行训练，最后进行预测；

在分类完成之后，便可以得到了某个网格属于特定建筑结构的概率：

步骤五：建筑立面的二分形状语义表达

从给定的建筑立面点云中得到表达墙面结构的最佳剖分树，就可以重建建筑立面，用数学语言描述如下：

给定建筑立面点云O，二分形状语义G＝(N，T，R，ω)，其中N是非终端形状语义的集合，如果一个带标号的矩形c(x，y，w，h)出现在一个规则的左边，那么它就是非终端形状语义；T是终端语义的集合，(x，y，w，h)定义了矩形在坐标系中的位置、宽度w和高度h；如果一个带标号的矩形c(x，y，w，h)不能出现在一个规则的左边，那么它就是终端形状语义；ω是一个特殊的非终端形状语义，称为初始形状，剖分总是从初始形状开始进行，R是一个有限的剖分规则的集合，设L(G)为二分形状语义的所有可能的剖分树(剖分方式)，寻找到一个最佳的剖分语言s∈L(G)，使得如下的目标函数最大：

$\underset{s}{\arg }\underset{s\∈L\left(G\right)}{\max }\underset{x,y\∈\mathrm{facade}}{\mathrm{\Σ}}p\left(s(x,y)\right|x,y)---\left(8\right)$

其中，s(x，y)表示剖分s得到的网格(x，y)的建筑结构类别；p(s(x，y)|x，y)即为步骤四得到的网格(x，y)处建筑类别为s(x，y)的概率，本申请专利采用了逐网格进行累加的方式对建筑立面上的所有网格进行判断，看剖分的结果与预分类的结果是否一致，若一致则进行累加1，否则累加0，最终的累加值则反映了该剖分的好坏，最大化式(15)得到建筑立面的最佳剖分；

步骤六：城市建筑立面三维建模

从二分语义的剖分过程来看，把求解(8)式看成是一个马尔科夫决策过程(MDP)；

(1)最优形状语义的求解

二分形状语义产生建筑的过程用MDP来建模，下面是MDP里的各个要素在建筑剖分中代表的含义；

智能体：MDP中的智能体表示的是形状语义剖分建筑立面的过程，在剖分建筑立面的时候，每一步都要处理一个非终端形状，智能体则决定采用哪种规则及规则参数应用于该形状；

环境：马尔可夫过程中的环境与智能体交互，指的是随着剖分进行不断建立的剖分树，随着剖分的不断进行，剖分树不断完善，环境也就有了对模型更多的知识；

状态：为剖分树中的结点，该结点表示了当前剖分所对应的形状及其位置，用state(x，y，symbol)表示，其中(x，y)为当前状态所处的位置，symbol表示当前形状的名字，symbol只能是非终端形状；

行为：为智能体可以选择的语法规则，用action(rule，parameter)表示，其中rule表示语法，parameter表示剖分的宽度；

回报值：设马尔可夫决策在t时刻执行动作action(rule，parameter)得到的终端形状为A(x，y，w，h，a)，那么相应的回报值为，

$R_t=\underset{i=x}{\overset{x+w}{\mathrm{\Σ}}}\underset{j=y}{\overset{y+h}{\mathrm{\Σ}}}p\left(a\right|x,y)---\left(9\right)$

其中p(a|x，y)表示建筑立面(x，y)处是a类别的概率；

设MDP决策过程的序列为T，那么长期的回报为：

$R=\underset{t=0}{\overset{T}{\mathrm{\Σ}}}{R}_t---\left(10\right)$

由于MDP决策完毕后，把建筑立面剖分完毕，因此有，

$R=\underset{t=0}{\overset{T}{\mathrm{\Σ}}}{R}_t=\underset{x,y\∈I}{\mathrm{\Σ}}p(x,y,s(x,y))---\left(11\right)$

MDP的最优解即是获取的最优形状语义；

采用增强学习算法求解以上叙述的MDP问题，增强学习采用的是Q-Learning算法；

(2)建筑立面深度信息的恢复

强化学习得到的最佳语义可以表示为剖分树的形式，其中剖分树的叶子节点是组成立面的所有建筑结构(终端形状)，记为Leaf_i(c，x，y，w，h)，c表示叶子节点i所属的建筑类别，(x，y，w，h)定义了建筑结构的外包矩形R，通过这些叶子节点，计算各类建筑结构在建筑立面上的分布情况，并构造出建筑该建筑立面的二维平面，还需要得到各种建筑结构的深度以完成最后的三维建模，为此，查询点云中R内的所有点，利用随机抽样一致模型得到这些点所的平面，设该平面为：

a_ix+b_iy+c_iz+d_i＝0 (12)

在上述坐标系转换中，已经求得建筑主平面所在的方程：

Ax+By+Cx+D＝0 (13)

由于建筑结构的平面与主平面平行，则建筑类别c离主平面的深度d为：

$d=\left|\frac{d_i-D}{\sqrt{A^2+B^2+C^2}}\right|---\left(14\right)$

如此便恢复了建筑结构c的深度信息。