CN109579784B - 基于数字表面模型的城市地区建筑物高度的自动获取方法 - Google Patents

Abstract

本发明公开了一种基于数字表面模型的城市地区建筑物高度的自动获取方法，其包括如下步骤：s1.通过空间层次分解和信息统计的方法确定非地面点；s2.非地面点剔除后，用周围高程拟合得到数字地形模型；s3.利用数字表面模型高程减去数字地形模型高程，得到城市地区建筑物的高度；s4.通过形态学变换方法去除城市地区建筑物高度中的噪声；通过超级像素分割将城市地区建筑物的高度分割为一定大小的区域；统计每个分割区域内的高度直方图；取高度直方图中峰值和最大值对应高度的平均值作为为该分割区域内建筑物的高度。本发明不仅解决了将建筑物的高程信息从包含地面高程的信息的数据中分离出来的技术问题；还解决了如何准确地获取建筑顶部的高度的技术问题。

Description

基于数字表面模型的城市地区建筑物高度的自动获取方法

技术领域

本发明属于数字表面模型(digital surface model，简称DSM)后处理及信息提取领域，特别涉及一种基于数字表面模型的城市地区建筑物高度的自动获取方法。

背景技术

传统的测绘成果主要为4D产品，包括数字高程模型(digital elevation model，简称DEM)、数字正射影像图(digital orthophoto map，简称DOM)、数字栅格地图(digitalraster graphic，简称DRG)以及数字线划地图(digital line graphic，简称DLG)。这些产品中主要包含地形高度信息和地物平面位置信息，不包含建筑物的高程信息，一些三维产品虽然含有地物的高程信息，但不包含建筑本身的高度信息。城市建筑的高度信息毫无疑问是非常有价值的信息，对城市规划、宏观经济分析等众多领域可以起到很大的帮助作用。

数字表面模型包含了地物、地形的高程信息，是获取建筑物高度的理想数据源。从数字表面模型中获取建筑物的高度信息难度很大，概括起来主要有两个难点：一是将建筑物的高程信息从包含地面高程的信息的数据中分离出来；二是准确地获取建筑顶部的高度。

发明内容

本发明的目的在于提出一种基于数字表面模型的城市地区建筑物高度的自动获取方法，以降低将建筑物的高程信息从包含地面高程的信息的数据中分离出来的难度。

本发明为了实现上述目的，采用如下技术方案：

基于数字表面模型的城市地区建筑物高度的自动获取方法，包括如下步骤：

s1.通过空间层次分解和信息统计的方法确定非地面点；

s2.非地面点剔除后，用周围高程拟合得到数字地形模型；

s3.利用数字表面模型高程减去数字地形模型高程，得到城市地区建筑物的高度。

优选地，步骤s1中，基于四度重叠空间层次分解和信息统计的方法确定非地面点，即：

利用重叠范围为每个窗口四分之一面积的四个重叠窗口内的高度进行空间层次分解；只有四个重叠窗口内都确定为非地面点的空间点，才作为最终的非地面点。

优选地，步骤s2具体为：

非地面点剔除后，对每个无效点向8邻域方向探测有效点，并记录对应的探测距离；

通过距离倒数加权法拟合非地面点的高程。

优选地，步骤s3之后还包括步骤：

s4.通过超级像素分割将城市地区建筑物的高度分割为一定大小的区域；

统计每个分割区域内的高度直方图；

取高度直方图中峰值和最大值对应高度的平均值作为为该分割区域内建筑物的高度。

优选地，在超级像素分割之前，通过形态学变换方法去除城市地区建筑物高度中的噪声。

优选地，步骤s4中，超级像素的数量与数字表面模型中像素的数量呈正比。

优选地，步骤s1具体为：

s1.1.设置初始参数，包括窗口参数初值cr0＝50，有效点数阈值初值minmx0＝800，高程层次划分参数hoff＝2.0米，迭代次数xhnum＝3；

s1.2.通过多次迭代，识别非地面点，对第mm次迭代：

设置窗口参数：cr＝cr0/(mm+1)；其中，mm＝1、2、3…；

有效点数阈值参数：minmx＝minmx0/((mm+1)*(mm+1))；

如果minmx小于25，则minmx设为25；

假设数字表面模型宽度为dcol，高度为drow；

列方向空间划分尺度cnum＝dcol/cr-1，

行方向空间划分尺度rnum＝drow/cr-1；

对四个重叠半个窗口尺度的窗口进行空间点数统计；

设起始行列坐标为(minr,minc)，对行方向第i个、列方向第j个起算窗口，对应的四个重叠窗口中第k个窗口的起始行列坐标为：

si＝minr+pyi_k+i*cr，sj＝minc+pyj_k+j*cr；

其中，pyi₀＝0，pyi₁＝0，pyi₂＝cr/2，pyi₃＝cr/2；pyj₀＝0，pyj₁＝cr/2，pyj₂＝0，pyj₃＝cr/2；

四个窗口中每个窗口包含cr*cr条记录，每个记录有四个标记，分别为窗口行号m，窗口列号n，窗口内(m，n)元素对应的数字表面模高度fH，是否为非地面点标记iferr；

将是否为非地面点标记iferr初始化为false；

统计窗口内的最大高程hmax和最小高程hmin；

计算高程划分层数tsnum＝INT((hmax-hmin)/hoff+1)，INT表示取整；

如果tsnum大于50，将tsnum设为50，同时，hofftemp＝(hmax-hmin)/(tsnum-1)，

否则hoff为初始值；

将点数num＝cr*cr平分到每一层，每层点数为pn＝num/tsnum；

如果pn为0，将pn设为1；

对第m层，

设置当前层点数初值thisnum＝0；

当前层最小高程为thishmin＝hmin+m*hofftemp；

当前层最大高程为thishmax＝thishmin+hofftemp；

当前窗口的每个像素都有一个点数记录orgi；

对每个像素p，如果不是非地面点：

如果该点高度位于当前层次空间，记录：orgi_thisnum＝p；

点数thisnum增加1；

如果thisnum小于minmx且thisnum大于0，将位于该层的所有点设为非地面点，窗口记录jlnum增加1；

对数字表面模型中的每个点，记录其平面坐标和高程；

如果某个数字表面模型格网点的记录数为4，则为非地面点，否则为地面点。

本发明具有如下优点：

本发明通过空间层次分解和信息统计的方法确定非地面点，剔除非地面点后，用周围高程拟合得到数字地形模型，进而得到初步的城市地区建筑物的高度，解决了如何将建筑物的高程信息从包含地面高程的信息的数据中分离出来的技术问题。此外，本发明在得到城市地区建筑物的高度后，通过形态学变换去除建筑物高度的噪声，然后利用超级像素分割和高度直方图统计确定建筑物的高度，解决了建筑顶部的高度准确获取的技术问题。

附图说明

图1为本发明中基于数字表面模型的城市地区建筑物高度的自动获取方法的流程图；

图2为本发明实施例中非地面点确定方法的示意图；

图3为本发明实施例中数字表面模型的效果示意图；

图4为本发明实施例中数字地形模型的效果示意图；

图5为本发明实施例中建筑物高度示意图；

图6为本发明实施例中超级像素分割效果图；

图7为本发明实施例中非地面点的高程拟合示意图。

具体实施方式

下面结合附图以及具体实施方式对本发明作进一步详细说明：

如图1所示，基于数字表面模型的城市地区建筑物高度的自动获取方法，包括如下步骤：

s1.通过空间层次分解和信息统计的方法确定非地面点。

如图2所示，本实施例中的空间层次分解基于四度重叠空间层次分解，即：

利用重叠范围为每个窗口四分之一面积(即重叠范围的长、宽方向分别为每个重叠窗口长、宽方向的二分之一长度，如图2所示)的四个重叠窗口内的高度进行空间层次分解；只有四个重叠窗口内都确定为非地面点的空间点，才作为最终的非地面点。

在图2中示出的四个重叠窗口分别为重叠窗口1、重叠窗口2、重叠窗口3和重叠窗口4。

上述四个重叠窗口例如可以通过下述方式实现：

首先将图2中的重叠窗口1向右移动半个窗口得到重叠窗口2，将重叠窗口1向下移动半个窗口得到重叠窗口3，将重叠窗口3向右移动半个窗口得到重叠窗口4。

通过上述方式，得到的重叠范围的面积为每个重叠窗口面积的四分之一。

本实施例在每个重叠窗口内分别做一次独立的分析，只有四个重叠窗口内都确定为非地面点的空间点，才作为最终的非地面点。

上述设计避免了地形阶梯状变化区域被错分为非地面点，因而能更好的保留地形细节。

将当前窗口的高程范围分为多层，图2中的高程层次1是指从下到上的第一层。

基于上述四度重叠空间层次分解，下面详述步骤s1的具体过程：

s1.1.设置初始参数，包括窗口参数初值cr0＝50，有效点数阈值初值minmx0＝800，高程层次划分参数hoff＝2.0米，迭代次数xhnum＝3。

s1.2.通过多次迭代，识别非地面点，对第mm次迭代：

设置窗口参数：cr＝cr0/(mm+1)；其中，mm＝1、2、3…；

有效点数阈值参数：minmx＝minmx0/((mm+1)*(mm+1))；

如果minmx小于25，则minmx设为25。

假设数字表面模型宽度为dcol，高度为drow；

列方向空间划分尺度cnum＝dcol/cr-1，

行方向空间划分尺度rnum＝drow/cr-1；

对四个重叠半个窗口尺度的窗口进行空间点数统计。

设起始行列坐标为(minr,minc)，对行方向第i个、列方向第j个起算窗口，对应的四个重叠窗口中第k个窗口的起始行列坐标为：

si＝minr+pyi_k+i*cr，sj＝minc+pyj_k+j*cr；

其中，pyi₀＝0，pyi₁＝0，pyi₂＝cr/2，pyi₃＝cr/2；pyj₀＝0，pyj₁＝cr/2，pyj₂＝0，pyj₃＝cr/2。

四个窗口中每个窗口包含cr*cr条记录，每个记录有四个标记，分别为窗口行号m，窗口列号n，窗口内(m，n)元素对应的数字表面模高度fH，是否为非地面点标记iferr。

将是否为非地面点标记iferr初始化为false。

统计窗口内的最大高程hmax和最小高程hmin。

计算高程划分层数tsnum＝INT((hmax-hmin)/hoff+1)，INT表示取整。

如果tsnum大于50，将tsnum设为50，同时，hofftemp＝(hmax-hmin)/(tsnum-1)，

否则hoff为初始值。

将点数num＝cr*cr平分到每一层，每层点数为pn＝num/tsnum；

如果pn为0，将pn设为1。

对第m层，

设置当前层点数初值thisnum＝0；

当前层最小高程为thishmin＝hmin+m*hofftemp；

当前层最大高程为thishmax＝thishmin+hofftemp。

当前窗口的每个像素都有一个点数记录orgi。

对每个像素p，如果不是非地面点：

如果该点高度位于当前层次空间，记录：orgi_thisnum＝p；

点数thisnum增加1；

如果thisnum小于minmx且thisnum大于0，将位于该层的所有点设为非地面点，窗口记录jlnum增加1；

对数字表面模型中的每个点，记录其平面坐标和高程；

如果某个数字表面模型格网点的记录数为4，则为非地面点，否则为地面点。

通过上述方法可准确确定非地面点，并且能更好的保留地形细节。

s2.非地面点剔除后，用周围高程拟合得到数字地形模型(digital terrinemodel，简称DTM)。

该步骤s2具体为：

非地面点剔除后，对每个无效点向8邻域方向探测有效点，并记录对应的探测距离；

通过距离倒数加权法拟合非地面点的高程。

其中，无效点是非地面点剔除后留下的窟窿点，而有效点则是指地面点。

如图7所示，黑色区域为有效，白色区域为检测出的无效区域。

对于无效区域内的任意一点O，向8邻域方向探测有效点。

设八个方向有效点离O的距离为D1、D2、D3、D4、D5、D6、D7、D8。

探测到的有效点高程分别为h1、h2、h3、h4、h5、h6、h7、h8。

则O点的高程Oh为：

s3.利用数字表面模型高程减去数字地形模型高程，得到城市地区含噪声的建筑物高度。

此步骤s3得到的建筑物高度，为包含噪声的建筑物高度的初步信息。

s4.通过形态学学变换方法去除城市地区建筑物高度中的噪声。

本实施例例如可以通过形态学腐蚀膨胀处理，以去除建筑物高度中的毛刺噪声。

去噪后，通过超级像素分割将城市地区建筑物的高度分割为一定大小的区域；其中，超级像素的数量与数字表面模型中像素的数量呈正比。

例如，设定超级像素的数量为N′，数字表面模型中像素的数量为N，

则N′＝N/800。

统计每个分割区域内的高度直方图，将高度直方图分为多个高度等级。

取高度直方图中峰值和最大值对应高度的平均值作为为该分割区域内建筑物的高度。

本实施例采用超级像素分割的方法将非地面点进行直方图统计分析，获取的建筑高度具有区域代表性，有效削弱了噪声的影响。

此外，本实施例还对本发明方法进行了实验，分别如图3至图6所示。

其中，图3为原始的数字表面模型(即DSM)的效果示意图。

第一步，利用本实施例中的方法确定非地面点。

第二步，非地面点剔除后，利用周围地面点高度拟合非地面点高度，得到图4的DTM。

第三步，图3的数字表面模型高程减去图4的数字地形模型高程得到建筑高度值。

在第三步中得到的建筑高度值含有噪声，需要进行去噪处理。

第四步，通过形态学的腐蚀、膨胀操作去除噪声得到图5的建筑高度图。

然后将去噪后的建筑高度图进行超级像素分割，并统计每个分割区域内的峰值高度和最大高度，取二者的均值作为该分割区域内的建筑物高度，如图6所示。

通过上述实验可以看出，本发明不仅解决了如何将建筑物的高程信息从包含地面高程的信息的数据中分离出来的问题；还解决了如何准确地获取建筑顶部的高度的问题。

当然，以上说明仅仅为本发明的较佳实施例，本发明并不限于列举上述实施例，应当说明的是，任何熟悉本领域的技术人员在本说明书的教导下，所做出的所有等同替代、明显变形形式，均落在本说明书的实质范围之内，理应受到本发明的保护。

Claims (6)

1.基于数字表面模型的城市地区建筑物高度的自动获取方法，其特征在于，

包括如下步骤：

s1.通过空间层次分解和信息统计的方法确定非地面点；

所述步骤s1中，基于四度重叠空间层次分解和信息统计的方法确定非地面点，即：

利用重叠范围为每个窗口四分之一面积的四个重叠窗口内的高度进行空间层次分解；只有四个重叠窗口内都确定为非地面点的空间点，才作为最终的非地面点；

s2.非地面点剔除后，用周围高程拟合得到数字地形模型；

s3.利用数字表面模型高程减去数字地形模型高程，得到城市地区建筑物的高度。

2.根据权利要求1所述的城市地区建筑物高度的自动获取方法，其特征在于，

所述步骤s2具体为：

非地面点剔除后，对每个无效点向8邻域方向探测有效点，并记录对应的探测距离；

通过距离倒数加权法拟合非地面点的高程。

3.根据权利要求1所述的城市地区建筑物高度的自动获取方法，其特征在于，

所述步骤s3之后还包括步骤：

s4.通过超级像素分割将城市地区建筑物的高度分割为一定大小的区域；

统计每个分割区域内的高度直方图；

取高度直方图中峰值和最大值对应高度的平均值作为该分割区域内建筑物的高度。

4.根据权利要求3所述的城市地区建筑物高度的自动获取方法，其特征在于，

在超级像素分割之前，通过形态学变换方法去除城市地区建筑物高度中的噪声。

5.根据权利要求3所述的城市地区建筑物高度的自动获取方法，其特征在于，

所述步骤s4中，超级像素的数量与数字表面模型中像素的数量呈正比。

6.根据权利要求1所述的城市地区建筑物高度的自动获取方法，其特征在于，

所述步骤s1具体为：

s1.1.设置初始参数，包括窗口参数初值cr0＝50，有效点数阈值初值minmx0＝800，高程层次划分参数hoff＝2.0米，迭代次数xhnum＝3；

s1.2.通过多次迭代，识别非地面点，对第mm次迭代：

设置窗口参数：cr＝cr0/(mm+1)；其中，mm＝1、2、3…；

有效点数阈值参数：minmx＝minmx0/((mm+1)*(mm+1))；

如果minmx小于25，则minmx设为25；

假设数字表面模型宽度为dcol，高度为drow；

列方向空间划分尺度cnum＝dcol/cr-1，

行方向空间划分尺度rnum＝drow/cr-1；

对四个重叠半个窗口尺度的窗口进行空间点数统计；

设起始行列坐标为(minr,minc)，对行方向第i个、列方向第j个起算窗口，对应的四个重叠窗口中第k个窗口的起始行列坐标为：

si＝minr+pyi_k+i*cr，sj＝minc+pyj_k+j*cr；

其中，pyi₀＝0，pyi₁＝0，pyi₂＝cr/2，pyi₃＝cr/2；pyj₀＝0，pyj₁＝cr/2，pyj₂＝0，pyj₃＝cr/2；

四个窗口中每个窗口包含cr*cr条记录，每个记录有四个标记，分别为窗口行号m，窗口列号n，窗口内(m，n)元素对应的数字表面模高度fH，是否为非地面点标记iferr；

将是否为非地面点标记iferr初始化为false；

统计窗口内的最大高程hmax和最小高程hmin；

计算高程划分层数tsnum＝INT((hmax-hmin)/hoff+1)，INT表示取整；

如果tsnum大于50，将tsnum设为50，同时，hofftemp＝(hmax-hmin)/(tsnum-1)，

否则hoff为初始值；

将点数num＝cr*cr平分到每一层，每层点数为pn＝num/tsnum；

如果pn为0，将pn设为1；

对第m层，

设置当前层点数初值thisnum＝0；

当前层最小高程为thishmin＝hmin+m*hofftemp；

当前层最大高程为thishmax＝thishmin+hofftemp；

当前窗口的每个像素都有一个点数记录orgi；

对每个像素p，如果不是非地面点：

如果该点高度位于当前层次空间，记录：orgi_thisnum＝p；

点数thisnum增加1；

如果thisnum小于minmx且thisnum大于0，将位于该层的所有点设为非地面点，窗口记录jlnum增加1；

对数字表面模型中的每个点，记录其平面坐标和高程；

如果某个数字表面模型格网点的记录数为4，则为非地面点，否则为地面点。

Images