发明名称 基于形状语义的建筑立面激光雷达点云解译与重建的方法
摘要 建筑立面建模是城市建模的重要组成部分。本发明基于形状语义和增强学优化策略,提出了一种建筑立面激光点云解译和重建的方法。首先基于数据驱动对点云进行建筑学语义上的预分类,接着设计了可以灵活应用于各种建筑立面风格的二分形状语义,最后将建筑立面的最佳形状语义求解的问题建模为马尔科夫决策过程,采用增强学算法结合支持向量机预分类的结果优化求解,实现对建筑立面的解译和三维重建。本发明提出的方法不仅能够较为精确的对建筑立面进行精细建模,而且能够适应于各种风格建筑立面的解译,对于建筑立面点云缺失的情况也有很好的鲁棒性。
申请公布号 CN104036550A 申请公布日期 2014.09.10
申请号 CN201410286813.1 申请日期 2014.06.25
申请人 北京师范大学 发明人 张立强;徐翔;张良
分类号 G06T17/05(2011.01)I 主分类号 G06T17/05(2011.01)I
代理机构 代理人
主权项 基于形状语义的建筑立面激光点云解译与重建的方法,其特征在于,包括如下步骤:步骤一:点云噪声去除对于输入点云的每一点,计算它到邻近点的距离的均值μ和标准差σ,那些距离落在μ±α·σ之外的点则被视为噪声而移除,其中α是和选取的邻近点的数目相近的一个值,本申请专利中,α=1.0,邻近点的数目为20;步骤二:坐标系转换把建筑立面所在的平面转换到X‑Y平面中去,即,建筑的水平方向与X轴平行,竖直方向与Y轴平行;建筑立面的坐标系变换陈述为如下问题,给定建筑点云当前所处的基M<sub>1</sub>={α<sub>1</sub>,α<sub>2</sub>,α<sub>3</sub>}以及目标坐标系的基M<sub>2</sub>={β<sub>1</sub>,β<sub>2</sub>,β<sub>3</sub>},求M<sub>1</sub>到M<sub>2</sub>的过渡矩阵,其中,β1={1,0,0},β2={0,1,0},β<sub>3</sub>={0,0,1},而要得到两组基之间的过渡矩阵P,只需要给出目标坐标系的基中的每个向量(β<sub>1</sub>,β<sub>2</sub>,β<sub>3</sub>)在当前所处基下的坐标(p<sub>1</sub>,p<sub>2</sub>,p<sub>3</sub>),依次以这些坐标为列向量组成矩阵即是从M<sub>1</sub>到M<sub>2</sub>的过渡矩阵P;基向量β<sub>1</sub>在当前基下对应的坐标p<sub>1</sub>是建筑立面水平方向的向量,基向量β<sub>2</sub>对应的坐标p<sub>2</sub>则是建筑立面竖直方向的向量,基向量p<sub>3</sub>对应的坐标则是建筑平面的法向量;采用随机抽样一致模型计算建筑立面的平面参数,继而得到其归一化的法向量,即为p<sub>3</sub>,对于p<sub>2</sub>,使用建筑墙面边缘线所处的向量表示,而p<sub>1</sub>则可以由p<sub>2</sub>和p<sub>3</sub>的外积求得,得到过渡矩阵P后,对原始点云中的点X应用变换P即可得到在新坐标系下的坐标X';步骤三:点云特征提取首先对建筑立面进行剖分,使得建筑立面由一系列三维的小网格组成,在本申请专利中,网格的大小设置为0.1米,设计了下面三种网格属性,用来通过支持向量机对建筑立面点云中的结构单元分类;(1)曲度定义网格内的曲度特征是网格内所有点的曲度的平均值,计算每点的曲度近似为计算该点和邻近的k个点近似组成的平面的曲度,二者可以通过分析邻近点的协相关矩阵的特征向量和特征矩阵得到;对于一点p<sub>i</sub>,它的协相关矩阵为:<maths num="0001" id="cmaths0001"><math><![CDATA[<mrow><mi>C</mi><mo>=</mo><mfrac><mn>1</mn><mi>k</mi></mfrac><munderover><mi>&Sigma;</mi><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>k</mi></munderover><mrow><mo>(</mo><msub><mi>p</mi><mi>i</mi></msub><mo>-</mo><mover><mi>p</mi><mo>&OverBar;</mo></mover><mo>)</mo></mrow><mo>&CenterDot;</mo><msup><mrow><mo>(</mo><msub><mi>p</mi><mi>i</mi></msub><mo>-</mo><mover><mi>p</mi><mo>&OverBar;</mo></mover><mo>)</mo></mrow><mi>T</mi></msup><mo>,</mo><mi>C</mi><mo>&CenterDot;</mo><mover><msub><mi>v</mi><mi>j</mi></msub><mo>&RightArrow;</mo></mover><mo>=</mo><msub><mi>&lambda;</mi><mi>j</mi></msub><mo>&CenterDot;</mo><mover><msub><mi>v</mi><mi>j</mi></msub><mo>&RightArrow;</mo></mover><mo>,</mo><mi>j</mi><mo>&Element;</mo><mo>{</mo><mn>0,1,2</mn><mo>}</mo><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>1</mn><mo>)</mo></mrow></mrow>]]></math><img file="FSA0000105507060000021.GIF" wi="1310" he="124" /></maths>其中,k是p<sub>i</sub>的邻近点的数目,<img file="FSA0000105507060000022.GIF" wi="47" he="74" />代表邻近点的三维中心坐标,λ<sub>j</sub>是协相关矩阵的第j个特征值,<img file="FSA0000105507060000023.GIF" wi="60" he="77" />是第j个特征向量;点p<sub>i</sub>的曲度为:γ<sub>pi</sub>=λ<sub>0</sub>/(λ<sub>0</sub>+λ<sub>1</sub>+λ<sub>2</sub>) (2)假设网格G中有n个点,那么该网格的曲度则表示为:<maths num="0002" id="cmaths0002"><math><![CDATA[<mrow><msub><mi>&gamma;</mi><mi>G</mi></msub><mo>=</mo><mfrac><mn>1</mn><mi>n</mi></mfrac><munderover><mi>&Sigma;</mi><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>n</mi></munderover><msub><mi>&gamma;</mi><mi>pi</mi></msub><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>3</mn><mo>)</mo></mrow></mrow>]]></math><img file="FSA0000105507060000024.GIF" wi="1259" he="119" /></maths>(2)密度设网格G中有n个点,点数最多网格内的点数为N,则该网格的密度d<sub>G</sub>表示为:d<sub>G</sub>=n/N (4)(3)深度假设网格G中有n个点,那么该网格的深度表示为:<maths num="0003" id="cmaths0003"><math><![CDATA[<mrow><msub><mi>s</mi><mi>G</mi></msub><mo>=</mo><mfrac><mn>1</mn><mi>n</mi></mfrac><munderover><mi>&Sigma;</mi><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>n</mi></munderover><msub><mi>z</mi><mi>i</mi></msub><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>5</mn><mo>)</mo></mrow></mrow>]]></math><img file="FSA0000105507060000025.GIF" wi="1182" he="120" /></maths>其中z<sub>i</sub>为点i的z坐标;以上三种特征将作为网格的特征向量用于步骤四的支持向量机分类中;步骤四:基于支持向量机的建筑立面点云分类首先选择训练集,训练集里的每一个样本都有一个类别标签及多个属性值,用数学的语言来描述就是:给定一个训练集(x<sub>i</sub>,y<sub>i</sub>),i=1,…,l,x<sub>i</sub>∈R<sup>n</sup>表示样本的特征向量,y∈{‑1,1}<sup>l</sup>表示样本的类别向量,支持向量机求解以下正定规划的最优解:<img file="FSA0000105507060000026.GIF" wi="1144" he="285" />训练向量的特征分别是曲度、密度和深度,本申请专利选择了半径基函数作为核函数,并设定γ=1,在步骤三中定义的三种网格属性分别归化到[‑1,1]之间,对样本集进行训练,最后进行预测;在分类完成之后,便可以得到了某个网格属于特定建筑结构的概率:<img file="FSA0000105507060000031.GIF" wi="1149" he="144" />步骤五:建筑立面的二分形状语义表达从给定的建筑立面点云中得到表达墙面结构的最佳剖分树,就可以重建建筑立面,用数学语言描述如下:给定建筑立面点云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),使得如下的目标函数最大:<maths num="0004" id="cmaths0004"><math><![CDATA[<mrow><munder><mi>arg</mi><mi>s</mi></munder><munder><mi>max</mi><mrow><mi>s</mi><mo>&Element;</mo><mi>L</mi><mrow><mo>(</mo><mi>G</mi><mo>)</mo></mrow></mrow></munder><munder><mi>&Sigma;</mi><mrow><mi>x</mi><mo>,</mo><mi>y</mi><mo>&Element;</mo><mi>facade</mi></mrow></munder><mi>p</mi><mrow><mo>(</mo><mi>s</mi><mrow><mo>(</mo><mi>x</mi><mo>,</mo><mi>y</mi><mo>)</mo></mrow><mo>|</mo><mi>x</mi><mo>,</mo><mi>y</mi><mo>)</mo></mrow><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>8</mn><mo>)</mo></mrow></mrow>]]></math><img file="FSA0000105507060000032.GIF" wi="1115" he="103" /></maths>其中,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),那么相应的回报值为,<maths num="0005" id="cmaths0005"><math><![CDATA[<mrow><msub><mi>R</mi><mi>t</mi></msub><mo>=</mo><munderover><mi>&Sigma;</mi><mrow><mi>i</mi><mo>=</mo><mi>x</mi></mrow><mrow><mi>x</mi><mo>+</mo><mi>w</mi></mrow></munderover><munderover><mi>&Sigma;</mi><mrow><mi>j</mi><mo>=</mo><mi>y</mi></mrow><mrow><mi>y</mi><mo>+</mo><mi>h</mi></mrow></munderover><mi>p</mi><mrow><mo>(</mo><mi>a</mi><mo>|</mo><mi>x</mi><mo>,</mo><mi>y</mi><mo>)</mo></mrow><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>9</mn><mo>)</mo></mrow></mrow>]]></math><img file="FSA0000105507060000041.GIF" wi="1158" he="141" /></maths>其中p(a|x,y)表示建筑立面(x,y)处是a类别的概率;设MDP决策过程的序列为T,那么长期的回报为:<maths num="0006" id="cmaths0006"><math><![CDATA[<mrow><mi>R</mi><mo>=</mo><munderover><mi>&Sigma;</mi><mrow><mi>t</mi><mo>=</mo><mn>0</mn></mrow><mi>T</mi></munderover><msub><mi>R</mi><mi>t</mi></msub><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>10</mn><mo>)</mo></mrow></mrow>]]></math><img file="FSA0000105507060000042.GIF" wi="1137" he="124" /></maths>由于MDP决策完毕后,把建筑立面剖分完毕,因此有,<maths num="0007" id="cmaths0007"><math><![CDATA[<mrow><mi>R</mi><mo>=</mo><munderover><mi>&Sigma;</mi><mrow><mi>t</mi><mo>=</mo><mn>0</mn></mrow><mi>T</mi></munderover><msub><mi>R</mi><mi>t</mi></msub><mo>=</mo><munder><mi>&Sigma;</mi><mrow><mi>x</mi><mo>,</mo><mi>y</mi><mo>&Element;</mo><mi>I</mi></mrow></munder><mi>p</mi><mrow><mo>(</mo><mi>x</mi><mo>,</mo><mi>y</mi><mo>,</mo><mi>s</mi><mrow><mo>(</mo><mi>x</mi><mo>,</mo><mi>y</mi><mo>)</mo></mrow><mo>)</mo></mrow><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>11</mn><mo>)</mo></mrow></mrow>]]></math><img file="FSA0000105507060000043.GIF" wi="1100" he="131" /></maths>MDP的最优解即是获取的最优形状语义;采用增强学习算法求解以上叙述的MDP问题,增强学习采用的是Q‑Learning算法;(2)建筑立面深度信息的恢复强化学习得到的最佳语义可以表示为剖分树的形式,其中剖分树的叶子节点是组成立面的所有建筑结构(终端形状),记为Leaf<sub>i</sub>(c,x,y,w,h),c表示叶子节点i所属的建筑类别,(x,y,w,h)定义了建筑结构的外包矩形R,通过这些叶子节点,计算各类建筑结构在建筑立面上的分布情况,并构造出建筑该建筑立面的二维平面,还需要得到各种建筑结构的深度以完成最后的三维建模,为此,查询点云中R内的所有点,利用随机抽样一致模型得到这些点所的平面,设该平面为:a<sub>i</sub>x+b<sub>i</sub>y+c<sub>i</sub>z+d<sub>i</sub>=0 (12)在上述坐标系转换中,已经求得建筑主平面所在的方程:Ax+By+Cx+D=0 (13)由于建筑结构的平面与主平面平行,则建筑类别c离主平面的深度d为:<maths num="0008" id="cmaths0008"><math><![CDATA[<mrow><mi>d</mi><mo>=</mo><mo>|</mo><mfrac><mrow><msub><mi>d</mi><mi>i</mi></msub><mo>-</mo><mi>D</mi></mrow><msqrt><msup><mi>A</mi><mn>2</mn></msup><mo>+</mo><msup><mi>B</mi><mn>2</mn></msup><mo>+</mo><msup><mi>C</mi><mn>2</mn></msup></msqrt></mfrac><mo>|</mo><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>14</mn><mo>)</mo></mrow></mrow>]]></math><img file="FSA0000105507060000044.GIF" wi="1247" he="149" /></maths>如此便恢复了建筑结构c的深度信息。
地址 100875 北京市海淀区新街口外大街19号科技处