CN106023312A

CN106023312A - Automatic 3D building model reconstruction method based on aviation LiDAR data

Info

Publication number: CN106023312A
Application number: CN201610319784.3A
Authority: CN
Inventors: 程亮; 许浩; 李满春; 王娅君; 谌颂; 魏晓燕; 袁; 袁一; 夏南; 孙越凡; 陈东
Original assignee: Yunnan Provincial Surveying And Mapping Information Archives (yunnan Basic Geographic Information Center); Nanjing University
Current assignee: Yunnan Provincial Surveying And Mapping Information Archives (yunnan Basic Geographic Information Center); Nanjing University
Priority date: 2016-05-13
Filing date: 2016-05-13
Publication date: 2016-10-12
Anticipated expiration: 2036-05-13
Also published as: CN106023312B

Abstract

The invention discloses an automatic three-dimensional building model reconstruction method based on aerial LiDAR data. The steps are: use reverse iterative mathematical morphology filtering and point cloud density-based methods to extract building roof point clouds from aerial LiDAR data; And optimize the roof patch; build a two-dimensional regular grid to resample different roof layers to obtain the interior points and edge points of the roof layer; optimize the edge points of different roof layers; connect the interior points and edge points of the roof layer to construct the roof surface and Finally, the 3D model reconstruction of the roof of the building is realized. Practice has proved that the invention can effectively reconstruct the three-dimensional model of the roof of the building, provides a new idea for the connection between different roof layers, and has high reconstruction accuracy of the three-dimensional model.

Description

Automatic reconstruction method of 3D building model based on aerial LiDAR data

技术领域technical field

本发明涉及一种基于航空LiDAR数据的三维建筑物模型自动重建方法，特别是涉及一种基于航空LiDAR数据的采取平滑策略和层间连接的三维建筑物模型自动重建方法。The invention relates to an automatic three-dimensional building model reconstruction method based on aerial LiDAR data, in particular to an automatic three-dimensional building model reconstruction method based on aerial LiDAR data using a smoothing strategy and interlayer connections.

背景技术Background technique

三维建筑物模型是建筑物三维结构信息的重要表达手段，在城市规划、灾害监测、通信设施建设和数字城市等领域具有非常广泛的应用。建筑物三维模型重建作为数字城市建设中最重要和最具挑战性的任务，其研究在过去几十年一直得到极大的关注。利用航空立体像对获取建筑物三维模型，是一种传统的摄影测量方法，依然在大量使用，但需要较多的人工干预，自动化程度不高。The 3D building model is an important means of expressing the 3D structure information of the building, and has a very wide range of applications in the fields of urban planning, disaster monitoring, communication facility construction and digital city. As the most important and challenging task in digital city construction, 3D model reconstruction of buildings has received great attention in the past decades. Using aerial stereo image pairs to obtain 3D models of buildings is a traditional photogrammetry method that is still widely used, but it requires more manual intervention and the degree of automation is not high.

近年来，航空LiDAR技术发展极快，且航空LiDAR数据已被广泛应用于地表探测(Vosselman G，2005)、特征检测(Tong L H，2013)、地物提取(Boyko A，2011)、三维模型重建(Cheng L，2012)等方面，显示出了巨大的应用前景。因此，航空LiDAR技术为建筑物的三维重建提供了一种可选方法。航空LiDAR设备能直接获取地面目标的三维信息，可提高自动化水平。但是，大量不规则点数据也给建筑物的模型重建带来了新的挑战。为此，如何有效利用航空LiDAR数据的优势，实现建筑物三维模型的自动化高质量重建，依然是一个值得深入研究的命题。In recent years, aerial LiDAR technology has developed rapidly, and aerial LiDAR data has been widely used in surface detection (Vosselman G, 2005), feature detection (Tong L H, 2013), ground object extraction (Boyko A, 2011), 3D model reconstruction (Cheng L, 2012) and other aspects, showing great application prospects. Therefore, aerial LiDAR technology provides an alternative method for 3D reconstruction of buildings. Aviation LiDAR equipment can directly obtain three-dimensional information of ground targets, which can improve the level of automation. However, a large amount of irregular point data also brings new challenges to the model reconstruction of buildings. For this reason, how to effectively use the advantages of aerial LiDAR data to realize the automatic high-quality reconstruction of 3D models of buildings is still a proposition worthy of in-depth study.

尽管利用航空LiDAR数据进行建筑物三维模型重建的研究越来越多，但是大多数的方法是通过提取建筑物屋顶轮廓线实现最终的模型重建。Zhou Q Y，Cheng L，Susaki J，程亮等人分别在《2.5D building modeling by discovering global regularities》、《Integration of LiDAR data and optical multi-view images for 3Dreconstruction of building roofs》、《Knowledge-based modeling of buildings indense urban areas by combining airborne LiDAR data and aerial images》、《集成多视航空影像与LiDAR数据重建三维建筑物模型》等文章中，提出使用轮廓线重建建筑物三维模型的方法。但通过轮廓线实现三维模型重建存在一些固有缺点：1)噪声影响。噪声的存在将使得提取得到的轮廓线不完整，即使通过拟合的方法形成完整轮廓线，与真实轮廓线相比则存在较大误差。2)点密度影响。点密度的大小直接影响轮廓线的提取，点密度较小则提取得到的轮廓线将不完整。3)层间连接关系。提取得到的不同屋顶层轮廓线仍然是离散状态，很难确定之间的连接关系。Although there are more and more researches on building 3D model reconstruction using aerial LiDAR data, most of the methods achieve the final model reconstruction by extracting building rooflines. Zhou Q Y, Cheng L, Susaki J, Cheng Liang and others respectively in "2.5D building modeling by discovering global regularities", "Integration of LiDAR data and optical multi-view images for 3D reconstruction of building roofs", "Knowledge-based modeling of In articles such as "buildings intense urban areas by combining airborne LiDAR data and aerial images" and "Integrating multi-view aerial images and LiDAR data to reconstruct 3D building models", a method of reconstructing 3D building models using contour lines is proposed. However, there are some inherent disadvantages in realizing 3D model reconstruction through contour lines: 1) Noise influence. The existence of noise will make the extracted contour line incomplete, even if the complete contour line is formed by the fitting method, there will be a large error compared with the real contour line. 2) The effect of point density. The size of the point density directly affects the extraction of the contour line. If the point density is small, the extracted contour line will be incomplete. 3) Inter-layer connection relationship. The extracted contour lines of different roof layers are still in a discrete state, and it is difficult to determine the connection relationship between them.

在提取建筑物轮廓线之前，主要是要实现建筑物屋顶面片的分割。当前，对建筑物屋顶面片的分割方法主要可以分为两大类：1)模型驱动方法。2)数据驱动方法。Mass(1999)，Oude(2009)，Hebel(2012)，Susaki(2013)，Huang(2013)，Henn(2013)等人分别使用第一种方法实现建筑物屋顶面片分割。这种方法的思路是首先确定待实验建筑物屋顶类型，然后在建筑物屋顶模型数据库中确定需要使用的屋顶模型。然而这种方法获得的面片分割结果准确性得不到保障，当待实验建筑物屋顶类型超出预先设定的屋顶类型时，分割结果将存在较大错误。第二种方法的思路是完全以数据驱动提取出屋顶面片，预先不对屋顶类型做任何的假设。2012年，Chen D等人利用渐进形态学滤波、区域生长算法和自适应的随机采样一致性算法完成建筑物屋顶面片的分割；2014年，Fan H C等人提出一种基于屋脊线层次分解的屋顶面片分割方法，利用连接性和共面性的特点，沿着屋脊线实现屋顶面片分割；2014年，Awrangjeb M等人将原始航空LiDAR点云分为地面点和非地面点，对于非地面点使用点的共面性和点的局部特征完成平面屋顶的分割。但这些方法受原始数据影响较大，得到的屋顶面片之间不存在连接关系，且屋顶面片的轮廓比较曲折。Before extracting the building outline, the main thing is to realize the segmentation of the building roof patch. At present, the segmentation methods for building roof patches can be mainly divided into two categories: 1) Model-driven methods. 2) Data-driven approach. Mass (1999), Oude (2009), Hebel (2012), Susaki (2013), Huang (2013), Henn (2013) and others used the first method to achieve building roof patch segmentation. The idea of this method is to first determine the roof type of the building to be tested, and then determine the roof model to be used in the building roof model database. However, the accuracy of the patch segmentation results obtained by this method cannot be guaranteed. When the roof type of the building to be tested exceeds the preset roof type, the segmentation results will have large errors. The idea of the second method is to extract roof patches completely data-driven, without making any assumptions about the roof type in advance. In 2012, Chen D et al. used progressive morphological filtering, region growing algorithm and adaptive random sampling consensus algorithm to complete the segmentation of building roof patches; in 2014, Fan H C et al. The roof patch segmentation method uses the characteristics of connectivity and coplanarity to achieve roof patch segmentation along the ridge line; in 2014, Awrangjeb M et al. divided the original aerial LiDAR point cloud into ground points and non-ground points. Ground points use the coplanarity of points and the local features of points to complete the segmentation of flat roofs. However, these methods are greatly affected by the original data, and there is no connection relationship between the obtained roof patches, and the outline of the roof patches is relatively tortuous.

通过航空LiDAR数据进行建筑物屋顶快速自动三维建模是自动建模发展方向。然而现阶段对这类方法的研究虽然较多，但是都存在一定的问题，如何充分利用航空LiDAR数据的优势，并实现准确度、精确度较高的自动三维重建仍然有待进一步的研究。Rapid automatic 3D modeling of building roofs through aerial LiDAR data is the development direction of automatic modeling. However, although there are many studies on such methods at this stage, there are still certain problems. How to make full use of the advantages of aerial LiDAR data and realize automatic 3D reconstruction with high accuracy and precision still needs further research.

发明内容Contents of the invention

本发明要解决的技术问题是：针对现有建筑物重建方法较多通过提取建筑物轮廓线实现，以及轮廓线提取存在的缺陷和轮廓线之间的拓扑关系较难确定的问题，本发明提供一种基于航空LiDAR数据的三维建筑物模型自动重建方法，该方法能够保证屋顶面片分割的准确性，并使用屋顶层重采样得到的内部点和边缘点快速、高效地实现了不同屋顶层之间的连接，高精度地重建了建筑物的三维模型，解决了不同屋顶层之间连接关系较难确定的问题。The technical problem to be solved by the present invention is: aiming at the problems that existing building reconstruction methods are mostly realized by extracting building contour lines, and the defects in contour line extraction and the topological relationship between contour lines are difficult to determine, the present invention provides An automatic 3D building model reconstruction method based on aerial LiDAR data, which can ensure the accuracy of roof patch segmentation, and use the internal points and edge points obtained by roof layer resampling to quickly and efficiently realize the relationship between different roof layers. The connection between them reconstructs the 3D model of the building with high precision, which solves the problem that it is difficult to determine the connection relationship between different roof layers.

本发明提供的一种基于航空LiDAR数据的三维建筑物模型自动重建方法，步骤如下：A method for automatic reconstruction of a three-dimensional building model based on aerial LiDAR data provided by the invention, the steps are as follows:

第一步、建筑物屋顶点云提取——提取建筑物点云，剔除建筑物墙面上的点云，得到建筑物屋顶点云；The first step, building roof point cloud extraction - extract the building point cloud, remove the point cloud on the building wall, and get the building roof point cloud;

第二步、分割屋顶面片——对建筑物屋顶点云进行屋顶面片分割，再对得到的屋顶面片进行平整；The second step is to segment the roof patch—to segment the roof patch of the building roof point cloud, and then level the obtained roof patch;

第三步、合并屋顶面片——当两个屋顶面片邻近并且所在平面存在交线，则将这两个屋顶面片合并形成屋顶层；The third step is to merge the roof patches - when two roof patches are adjacent and there is an intersection line in the plane, the two roof patches are merged to form a roof layer;

第四步、屋顶层重采样——对所有屋顶层进行重采样，获得屋顶层的内部点和边缘点；The fourth step, roof layer resampling - resample all roof layers to obtain the interior points and edge points of the roof layer;

第五步、建筑物模型重建——对属于同一屋顶层的内部点和边缘点构建三角网，形成屋顶面；对属于相邻屋顶且位于同一竖直平面上的边缘点构建三角网，形成屋顶面建筑物墙面，最终完成建筑物三维模型重建。The fifth step, building model reconstruction - construct a triangular network for the internal points and edge points belonging to the same roof layer to form a roof surface; construct a triangular network for edge points belonging to adjacent roofs and located on the same vertical plane to form a roof The building wall surface, and finally complete the reconstruction of the 3D model of the building.

本发明还具有一下进一步的特征：The present invention also has following further features:

1、所述第一步中，剔除建筑物墙面上的点云方法如下：以r为半径，对每一个LiDAR点搜寻邻域点，用搜寻得到的邻域点的数量除以邻域区域的体积得到每一个LiDAR点的点云密度，当LiDAR点的点云密度小于指定阈值时，该LiDAR点为墙面点，进行剔除。1. In the first step, the method of removing point clouds on the building wall is as follows: take r as the radius, search for neighborhood points for each LiDAR point, and divide the number of neighborhood points obtained by searching by the neighborhood area Get the point cloud density of each LiDAR point. When the point cloud density of the LiDAR point is less than the specified threshold, the LiDAR point is a wall point and is eliminated.

2、所述第二步中，选取种子区域，使用区域生长算法实现屋顶面片分割，种子区域的选取方法如下：首先估算每一个LiDAR点的法向量，以计算点p为例，寻找以点p为中心，半径r范围内的点集N_p，根据公式(1)和(2)计算得到三个特征值λ₁、λ₂、λ₃，三个特征值中最小的特征值λ_min所对应的特征向量即为点p的法向量。2. In the second step, the seed area is selected, and the roof patch is segmented using the region growing algorithm. The selection method of the seed area is as follows: first estimate the normal vector of each LiDAR point, take the calculation point p as an example, and find the point p is the center and the point set N _p within the radius r is calculated according to the formulas (1) and (2) to obtain three eigenvalues λ ₁ , λ ₂ , λ ₃ , and the smallest eigenvalue λ _min among the three eigenvalues is The corresponding eigenvector is the normal vector of point p.

${C C}_{p p} = = \frac{11}{| | {N N}_{P P} | |} {Σ Σ}_{i i = = 11}^{n no} (({q q}_{i i} - - \overset{&OverBar; &OverBar;}{p p})) {(({q q}_{i i} - - \overset{&OverBar; &OverBar;}{p p}))}^{T T} - - - - - - ((11))$

$\overset{&OverBar; &OverBar;}{p p} = = \frac{11}{n no} {Σ Σ}_{i i = = 11}^{n no} {q q}_{i i} - - - - - - ((22))$

式中q_i∈N_p，n是点集N_p中点的数量，C_p是点p的协方差矩阵，根据公式计算点p的曲率，当点p的曲率小于指定曲率阈值λ_T时，就可认为点p的邻域点N_p在一个平面上，选取曲率所对应的邻域点N_p作为满足要求的种子区域。where q _i ∈ N _p , n is the number of points in the point set N _p , C _p is the covariance matrix of point p, according to the formula Calculate the curvature of point p, when the curvature of point p When it is less than the specified curvature threshold λ _T , it can be considered that the neighborhood point N _p of point p is on a plane, and the curvature The corresponding neighborhood point N _p is used as the seed area that meets the requirements.

3、所述第二步中，区域生长算法需要确定两个标准：一个是局内点的个数，另一个是拟合平面的标准差；当LiDAR点到种子区域所在平面的距离小于指定距离阈值时，认为是局内点；所有点计算完毕后，再计算局内点的标准差，当标准差小于指定阈值时，认为满足屋顶面片分割要求。3. In the second step, the region growing algorithm needs to determine two criteria: one is the number of points in the local area, and the other is the standard deviation of the fitting plane; when the distance from the LiDAR point to the plane where the seed area is located is less than the specified distance threshold , it is considered as an intra-office point; after all the points are calculated, the standard deviation of the intra-office points is calculated, and when the standard deviation is less than the specified threshold, it is considered to meet the roof patch segmentation requirements.

4、所述第二步中，面片平整优化包括分割面片点云平整和局部突出点云平整：4. In the second step, the optimization of patch smoothing includes segmented patch point cloud smoothing and local prominent point cloud smoothing:

对于分割面片点云平整，使用最小二乘法拟合出分割面片点云方程，然后根据方程将点平整到拟合面片上，拟合平面方程根据公式(3)进行计算：For the leveling of the segmented patch point cloud, the least square method is used to fit the segmented patch point cloud equation, and then the points are leveled to the fitted patch according to the equation, and the fitting plane equation is calculated according to formula (3):

$|\begin{matrix} Σ Σ {x x}_{i i}^{22} & Σ Σ {x x}_{i i} {y the y}_{i i} & Σ Σ {x x}_{i i} \\ Σ Σ {x x}_{i i} {y the y}_{i i} & Σ Σ {y the y}_{i i}^{22} & Σ Σ {y the y}_{i i} \\ Σ Σ {x x}_{i i} & Σ Σ {y the y}_{i i} & n no \end{matrix}| (\begin{matrix} {a a}_{00} \\ {a a}_{11} \\ {a a}_{22} \end{matrix}) = = (\begin{matrix} Σ Σ {x x}_{i i} {z z}_{i i} \\ Σ Σ {y the y}_{i i} {z z}_{i i} \\ Σ Σ {z z}_{i i} \end{matrix}) - - - - - - ((33))$

式中x_i，y_i，z_i表示点的三维坐标，a₀，a₁，a₂表示拟合平面方程的系数；In the formula, x _i , y _i , _zi represent the three-dimensional coordinates of the point, a ₀ , a ₁ , a ₂ represent the coefficients of the fitting plane equation;

对于局部突出点云的平整，利用距离阈值对局部突出点云进行搜索，计算不在分割面片中的点到分割面片的距离，当距离小于阈值时，记为所需点，将所需点按照上述得到的拟合平面方程平整到分割面片上。For the smoothing of the local prominent point cloud, the distance threshold is used to search the local prominent point cloud, and the distance from the point not in the segmented patch to the segmented patch is calculated. When the distance is less than the threshold, it is recorded as the required point, and the required point is According to the fitting plane equation obtained above, it is flattened onto the split surface.

5、所述第四步中，屋顶层重采样方法如下：首先将点云投影至XY平面，并对XY平面构建二维格网，然后进行如下判断：5. In the fourth step, the roof layer resampling method is as follows: first, project the point cloud to the XY plane, and construct a two-dimensional grid on the XY plane, and then make the following judgments:

(a)、当四邻域格网单元中的LiDAR点与中心格网单元中的LiDAR点属于同一个屋顶层，此时中心格网单元的中心作为屋顶层的内部点，该内部点的高度为中心格网单元中LiDAR点的平均高度；(a) When the LiDAR points in the four-neighborhood grid unit and the LiDAR point in the central grid unit belong to the same roof layer, then the center of the central grid unit is used as the internal point of the roof layer, and the height of the internal point is Average height of the LiDAR points in the central grid cell;

(b)、当中心格网单元的四邻域格网单元中的某一个邻域格网单元中的LiDAR点属于两个或两个以上屋顶层时，设该邻域格网单元内点数最多的两个屋顶层分别为屋顶层j和屋顶层k，获取该邻域格网单元内屋顶层j的LiDAR点和屋顶层k的LiDAR点的分割线，然后计算该分割线与中心格网单元中线在该邻域格网单元内的交点坐标(x₀，y₀)，则(x₀，y₀，z_j)为屋顶层j的边缘点，z_j为屋顶层j的LiDAR点的平均高度，则(x₀，y₀，z_k)为屋顶层k的边缘点，z_k为屋顶层k的LiDAR点的平均高度；(b) When the LiDAR points in one of the four neighborhood grid units of the central grid unit belong to two or more roof layers, set the largest number of points in the neighborhood grid unit The two roof layers are roof layer j and roof layer k respectively, and the dividing line between the LiDAR point of roof layer j and the LiDAR point of roof layer k in the neighborhood grid unit is obtained, and then the dividing line and the center line of the central grid unit are calculated The intersection point coordinates (x ₀ , y ₀ ) in the neighborhood grid unit, then (x ₀ , y ₀ , z _j ) is the edge point of the roof layer j, and z _j is the average height of the LiDAR points of the roof layer j , then (x ₀ , y ₀ , z _k ) is the edge point of the roof layer k, and z _k is the average height of the LiDAR points of the roof layer k;

(c)、当中心格网单元内的LiDAR点和四邻域格网单元中的某一个邻域格网单元中的LiDAR点分别属于两个不同屋顶层时，设两个不同屋顶层为屋顶层j和屋顶层k，计算中心格网单元和所述邻域格网单元之间的格网边线的中点坐标(x₁，y₁)，则(x₁，y₁，z_j)为屋顶层j的边缘点，z_j为屋顶层j的LiDAR点的平均高度，则(x₁，y₁，z_k)为屋顶层k的边缘点，z_k为屋顶层k的LiDAR点的平均高度。(c) When the LiDAR points in the central grid unit and the LiDAR points in one of the four neighborhood grid units belong to two different roof layers respectively, set the two different roof layers as the roof layer j and roof layer k, calculate the midpoint coordinates (x ₁ , y ₁ ) of the grid edge between the central grid unit and the neighborhood grid unit, then (x ₁ , y ₁ , z _j ) is the roof The edge point of layer j, z _j is the average height of LiDAR points of roof layer j, then (x ₁ , y ₁ , z _k ) is the edge point of roof layer k, z _k is the average height of LiDAR points of roof layer k .

6、所述二维格网的大小为LiDAR点云平均间距的2倍。6. The size of the two-dimensional grid is twice the average spacing of the LiDAR point cloud.

7、所述第四步中，获取邻域格网单元内屋顶层j的LiDAR点和屋顶层k的LiDAR点的分割线的方法是：将邻域格网单元内屋顶层j的LiDAR点和屋顶层k的LiDAR点视为两个类别，使用支持向量机算法获取分割线。7. In the fourth step, the method of obtaining the dividing line between the LiDAR point of the roof layer j in the neighborhood grid unit and the LiDAR point of the roof layer k is: the LiDAR point of the roof layer j in the neighborhood grid unit and The LiDAR points of the roof layer k are regarded as two categories, and the support vector machine algorithm is used to obtain the dividing line.

本发明的有益效果如下：The beneficial effects of the present invention are as follows:

(a)、本发明提出“种子区域选取—屋顶面片生长—面片平整优化”的策略，该策略在提取屋顶面片的同时，也对屋顶面片进行平整优化，有利于后续三维模型重建。(a), the present invention proposes a strategy of "seed area selection - roof surface growth - surface leveling optimization". This strategy not only extracts the roof surface, but also optimizes the leveling of the roof surface, which is beneficial to the subsequent 3D model reconstruction .

(b)、本发明利用屋顶层重采样计算得到内部点和边缘点，避免了传统方法提取屋顶轮廓线后，较难确定不同轮廓线之间的拓扑关系的问题。且内部点和边缘点的计算高效、便捷，避免了繁琐的工作，提高了自动化水平。(b), the present invention uses roof layer resampling calculation to obtain internal points and edge points, avoiding the problem that it is difficult to determine the topological relationship between different contour lines after the traditional method extracts the roof contour lines. Moreover, the calculation of internal points and edge points is efficient and convenient, avoiding tedious work and improving the automation level.

综上所述，本发明公开一种基于航空LiDAR数据的采取平滑策略和层间连接的三维建筑物模型自动重建方法，该方法通过建筑物屋顶点云提取、屋顶面片分割、屋顶面片合并、屋顶层重采样、建筑物模型重建五个步骤实现了建筑物屋顶的自动三维重建。试验证明，该方法具有较高的重建精度，能很好地反映建筑物顶部的几何结构特征。In summary, the present invention discloses a method for automatic reconstruction of a three-dimensional building model based on aerial LiDAR data using a smoothing strategy and interlayer connections. Automatic 3D reconstruction of building roofs is realized in five steps, namely, roof layer resampling and building model reconstruction. Experiments have proved that the method has high reconstruction accuracy and can well reflect the geometric structure characteristics of the top of the building.

附图说明Description of drawings

下面结合附图对本发明方法作进一步的说明。The method of the present invention will be further described below in conjunction with the accompanying drawings.

图1为本发明实施例的总体流程图。Fig. 1 is an overall flow chart of the embodiment of the present invention.

图2为本发明实施例区域航空LiDAR数据示意图。Fig. 2 is a schematic diagram of regional aviation LiDAR data according to an embodiment of the present invention.

图3为本发明实施例面片平整流程示意图。Fig. 3 is a schematic diagram of a smoothing process of a dough sheet according to an embodiment of the present invention.

图4为本发明实施例内部点和边缘点示意图。Fig. 4 is a schematic diagram of internal points and edge points according to an embodiment of the present invention.

图5a～图5e为本发明实施例内部点和边缘点计算示意图。5a to 5e are schematic diagrams of calculating internal points and edge points according to an embodiment of the present invention.

图6为本发明实施例内部点和边缘点最终结果示意图。Fig. 6 is a schematic diagram of the final result of internal points and edge points according to an embodiment of the present invention.

图7为本发明实施例整个区域模型重建结果图。Fig. 7 is a diagram showing the reconstruction result of the whole region model according to the embodiment of the present invention.

图8a～图8f为本发明实施例中6栋建筑物模型重建结果图。8a to 8f are reconstruction results of six building models in the embodiment of the present invention.

具体实施方式detailed description

下面根据附图详细阐述本发明，使本发明的技术路线和操作步骤更加清晰。The present invention will be described in detail below according to the accompanying drawings, so as to make the technical route and operation steps of the present invention clearer.

使用Optech ALTM Gemini激光扫描仪获取得到的航空LiDAR数据为实验数据，如图2所示。航空LiDAR数据平均点间距0.4m，高程精度15cm，平面精度30cm。The aerial LiDAR data obtained by using the Optech ALTM Gemini laser scanner is the experimental data, as shown in Figure 2. The average point spacing of aviation LiDAR data is 0.4m, the elevation accuracy is 15cm, and the plane accuracy is 30cm.

本实施例是一种基于航空LiDAR数据的三维建筑物模型自动重建方法(流程图见图1)，包括以下步骤：The present embodiment is a method for automatic reconstruction of a three-dimensional building model based on aerial LiDAR data (flow chart is shown in Fig. 1), comprising the following steps:

第一步、建筑物屋顶点云提取——原始航空LiDAR数据所包含的地物类型主要是建筑物、植被和地面三类。本发明基于专利CN201110432421.8中的反向迭代数学形态学滤波算法提取建筑物LiDAR点云。The first step is the point cloud extraction of building roofs - the types of ground objects contained in the original aerial LiDAR data are mainly buildings, vegetation and ground. The invention extracts building LiDAR point clouds based on the reverse iterative mathematical morphology filtering algorithm in the patent CN201110432421.8.

方法描述如下：对无规则分布的原始LiDAR点云进行重采样(重采样间距设定为1m)，对重采样后的等间距点进行反向迭代数学形态学滤波，以较小的步长逐渐减小滤波窗口的大小，对每个窗口都使用形态学“开”操作，同时不断对相邻两个窗口的滤波结果作差，当差值大于最小建筑物的高度时，对应的点就被标记为非地面点。提取得到的非地面点包括建筑物区域和密集树木区域，接着利用点云高程的粗糙度(高程的方差)，设定阈值将粗糙度较高的树木点云剔除。最后，提取出建筑物屋顶LiDAR点云。本实施例中，选择的最大窗口为106m，最小窗口为6m，窗口减小步长为10m，最小建筑物的高度为3m，高程方差阈值为0.4m。The method is described as follows: resample the original LiDAR point cloud with irregular distribution (the resampling interval is set to 1m), perform reverse iterative mathematical morphology filtering on the resampled equidistant points, and gradually Reduce the size of the filtering window, use the morphological "open" operation for each window, and continuously make a difference between the filtering results of two adjacent windows. When the difference is greater than the height of the minimum building, the corresponding point will be Marked as non-ground points. The extracted non-ground points include the building area and the dense tree area, and then use the roughness of the point cloud elevation (the variance of the elevation) to set a threshold to remove the tree point cloud with higher roughness. Finally, the building roof LiDAR point cloud is extracted. In this embodiment, the selected maximum window is 106 m, the minimum window is 6 m, the window reduction step is 10 m, the minimum building height is 3 m, and the elevation variance threshold is 0.4 m.

在提取得到的建筑物LiDAR点云中，仍然存在墙面点云，本发明使用基于密度的方法剔除墙面点云。In the extracted building LiDAR point cloud, there are still wall point clouds, and the present invention uses a density-based method to remove the wall point cloud.

方法描述如下：对于一个点p，可在r半径范围内得到点p的邻近点N_p，对于点p的密度可以定义为N_p中点的数量除以r半径范围的体积。当点p的密度小于密度阈值时，就认为点p是墙面点。反之，则为建筑物屋顶点。本实施例中，r半径大小为1m，密度阈值为2个点/m³。The method is described as follows: for a point p, the neighboring points N _p of point p can be obtained within the radius r, and the density of point p can be defined as the number of points in N _p divided by the volume of the radius r. When the density of point p is less than the density threshold, point p is considered to be a wall point. Otherwise, it is the roof point of the building. In this embodiment, the r radius is 1 m, and the density threshold is 2 points/m ³ .

第二步、分割屋顶面片——提取出建筑物屋顶LiDAR点云之后，需要对建筑物的不通屋顶面片进行提取。屋顶面片提取采用“种子区域选取—屋顶面片生长—面片平整优化”的策略，首先通过点曲率的估算选取种子区域，然后使用区域生长算法实现面片生长，最后对得到的屋顶面片进行平整优化，如图3所示。The second step is to segment the roof patch - after extracting the LiDAR point cloud of the building roof, it is necessary to extract the unreachable roof patch of the building. The roof patch extraction adopts the strategy of "seed area selection - roof patch growth - surface patch leveling optimization". First, the seed area is selected by point curvature estimation, and then the region growing algorithm is used to realize the patch growth. Finally, the obtained roof patch Perform leveling optimization, as shown in Figure 3.

建筑物屋顶面片平整的具体步骤如下：The specific steps for building roof surface leveling are as follows:

a1)、利用曲率估算选取种子区域——对于点的曲率估算，首先需要估算点的法向量。例如估算点p的法向量，需要确定点p的邻域点N_p＝{q|q∈P,d(p,q)<r}，其中r是搜索邻域点的搜索半径，它决定了邻域的大小。可以通过计算领域点的协方差矩阵来获得点的法向量，协方差矩阵定义如下：a1), using curvature estimation to select the seed area - for the curvature estimation of a point, the normal vector of the point needs to be estimated first. For example, to estimate the normal vector of point p, it is necessary to determine the neighborhood point N _p ={q|q∈P,d(p,q)<r} of point p, where r is the search radius of the search neighborhood point, which determines The size of the neighborhood. The normal vector of a point can be obtained by calculating the covariance matrix of the domain point, and the covariance matrix is defined as follows:

${C C}_{p p} = = \frac{11}{| | {N N}_{P P} | |} {Σ Σ}_{i i = = 11}^{n no} (({q q}_{i i} - - \overset{&OverBar; &OverBar;}{p p})) {(({q q}_{i i} - - \overset{&OverBar; &OverBar;}{p p}))}^{T T};; \overset{&OverBar; &OverBar;}{p p} = = \frac{11}{n no} {Σ Σ}_{i i = = 11}^{n no} {q q}_{i i}$

式中q_i∈N_p，n为点集N_p中点的数量，C_p为点p的协方差矩阵。通过上面的协方差矩阵可以计算得到三个特征值，分别为λ₁、λ₂、λ₃。假设λ₁<λ₂<λ₃，最小的特征值所对应的特征向量就为点p的法向量。然后再通过计算点的曲率来选取种子区域，曲率定义如下：In the formula, q _i ∈ N _p , n is the number of points in the point set N _p , and C _p is the covariance matrix of point p. Through the above covariance matrix, three eigenvalues can be calculated, which are λ ₁ , λ ₂ , and λ ₃ . Assuming λ ₁ <λ ₂ <λ ₃ , the eigenvector corresponding to the smallest eigenvalue is the normal vector of point p. Then select the seed area by calculating the curvature of the point. The curvature is defined as follows:

$\overset{&OverBar; &OverBar;}{λ λ} = = {λ λ}_{11} / / (({λ λ}_{11} + + {λ λ}_{22} + + {λ λ}_{33}))$

当小于曲率阈值λ_T时，就可以认为点p的邻域点在一个平面上，选取较小曲率所对应的邻域点作为满足要求的种子区域。本实施例中，曲率阈值为0.005。when When it is less than the curvature threshold λ _T , it can be considered that the neighborhood points of point p are on the same plane, and a smaller curvature is selected The corresponding neighborhood points are used as the seed area that meets the requirements. In this embodiment, the curvature threshold is 0.005.

a2)、面片生长——种子区域获得后，即获得了初始平面参数。这里定义两个满足面片生长的标准：一个是局内点的个数，另一个是分割平面的标准差。当点p_i到初始平面的距离小于距离阈值时，认为点p_i属于初始平面所在的平面，即为局内点，当没有点再满足到初始平面的距离小于距离阈值时，面片生长完毕，统计局内点的个数。并计算生长获得的面片的标准差，当标准差小于设定阈值时，认为屋顶面片提取完成，并符合提取要求。本实施例中，距离阈值为0.5m，标准差为0.95m。a2), patch growth—after the seed area is obtained, the initial plane parameters are obtained. Two criteria are defined here to meet the patch growth: one is the number of interior points, and the other is the standard deviation of the split plane. When the distance from point p _i to the initial plane is less than the distance threshold, it is considered that point p _i belongs to the plane where the initial plane is located, that is, it is an intra-local point. When no point satisfies the distance to the initial plane less than the distance threshold, the patch growth is complete. The number of points in the statistics bureau. And calculate the standard deviation of the surface patches obtained by growth. When the standard deviation is less than the set threshold, it is considered that the roof surface patch extraction is completed and meets the extraction requirements. In this embodiment, the distance threshold is 0.5m, and the standard deviation is 0.95m.

a3)、面片平整优化——优化主要针对两点：一是对分割面片的三维点集的平整；二是对细部干扰信息的平整。a3) Surface patch smoothing optimization—optimization is mainly aimed at two points: one is the smoothing of the three-dimensional point set of the divided patch; the other is the smoothing of detailed interference information.

本实施例在上述步骤a3)中的面片平整优化的具体方法如下：In the present embodiment, the specific method for smoothing and optimizing the surface in the above-mentioned step a3) is as follows:

b1)、三维点集的平整——在实际情况中，即使是平面屋顶，点云也不完全处于同一个二维平面上，即在提取得到的屋顶面片上，点存在波动。可以使用平面方程z＝a₀x+a₁y+a₂进行拟合，对于参数a₀，a₁，a₂的计算可以使用如下线性方程组进行计算：b1) Flattening of the 3D point set——In actual situations, even if it is a flat roof, the point cloud is not completely on the same two-dimensional plane, that is, there are fluctuations in the points on the extracted roof patch. The plane equation z=a ₀ x+a ₁ y+a ₂ can be used for fitting, and the calculation of parameters a ₀ , a ₁ and a ₂ can be performed using the following linear equations:

$|\begin{matrix} Σ Σ {x x}_{i i}^{22} & Σ Σ {x x}_{i i} {y the y}_{i i} & Σ Σ {x x}_{i i} \\ Σ Σ {x x}_{i i} {y the y}_{i i} & Σ Σ {y the y}_{i i}^{22} & Σ Σ {y the y}_{i i} \\ Σ Σ {x x}_{i i} & Σ Σ {y the y}_{i i} & n no \end{matrix}| (\begin{matrix} {a a}_{00} \\ {a a}_{11} \\ {a a}_{22} \end{matrix}) = = (\begin{matrix} Σ Σ {x x}_{i i} {z z}_{i i} \\ Σ Σ {y the y}_{i i} {z z}_{i i} \\ Σ Σ {z z}_{i i} \end{matrix})$

得到平面方程之后，可以把分割面片上的点按照得到的平面方程平整到同一个平面上。After obtaining the plane equation, the points on the split surface can be flattened to the same plane according to the obtained plane equation.

b2)、细部干扰信息的平整——对于提取得到的屋顶面片，存在点云缺失，这一部分点云需要被平整到对应屋顶面片上，对于每一个已知的分割面片向上寻找在屋顶面片提取结果中不存在的LiDAR点，找到后根据对应屋顶面片的平面方程将其平整到屋顶面片上。本实施例中向上寻找点的距离为2m。b2) Smoothing of detailed interference information - for the extracted roof patch, there is a lack of point cloud, this part of the point cloud needs to be leveled to the corresponding roof patch, and for each known segmented patch, search upwards on the roof surface The LiDAR points that do not exist in the patch extraction results are found and flattened to the roof patch according to the plane equation of the corresponding roof patch. In this embodiment, the distance to search for the point upwards is 2m.

第三步、合并屋顶面片——在屋顶层重采样之前，需要对所有屋顶面片按照一定规则进行合并以形成屋顶层。当两个屋顶面片满足两个条件：(1)两个屋顶面片所在平面存在交线；2)两个屋顶面片相互邻近时，两个屋顶面片可合并成一个屋顶层。The third step is to merge roof patches - before the roof layer is resampled, all roof patches need to be merged according to certain rules to form a roof layer. When two roof patches meet two conditions: (1) there is an intersection line in the plane where the two roof patches are located; 2) when the two roof patches are adjacent to each other, the two roof patches can be merged into one roof layer.

第四步、屋顶层重采样——计算屋顶层的内部点和边缘点，对k个屋顶层计算k个内部点或者边缘点，这k个点具有相同的x，y坐标，不同的高度值，如图4所示为内部点和边缘点示意图。屋顶层重采样基于二维规则格网进行。The fourth step, roof layer resampling-calculate the interior points and edge points of the roof layer, and calculate k interior points or edge points for k roof layers. These k points have the same x, y coordinates, and different height values , as shown in Figure 4 for the schematic diagram of internal points and edge points. Roof layer resampling is based on a 2D regular grid.

屋顶层重采样的具体步骤如下：The specific steps of roof layer resampling are as follows:

c1)、二维规则格网构建——二维规则格网的大小由LiDAR点云的在x和y方向上的四至坐标来确定，本实例中格网单元大小为1m。计算每一个点对应所属的格网单元，最终建立格网索引。c1), two-dimensional regular grid construction——the size of the two-dimensional regular grid is determined by the four-dimensional coordinates of the LiDAR point cloud in the x and y directions, and the grid unit size in this example is 1m. Calculate the corresponding grid unit of each point, and finally establish the grid index.

c2)、内部点和边缘点计算——每一个格网单元中存在一个及以上数量的内部点或者边缘点，在计算的过程中，需要考虑该格网单元的四个邻域格网单元。内部点和边缘点的计算有五种情况，①5个格网单元中的点都属于同一个屋顶层，②存在不同屋顶层的邻域格网单元在中心格网单元的左边或右边，③存在不同屋顶层的邻域格网单元在中心格网单元的上部或下部，④边缘点的水平位置在中心格网单元的边界上，⑤邻域格网单元中存在属于两个以上的不同屋顶层。c2) Calculation of internal points and edge points - there are one or more internal points or edge points in each grid unit, and the four neighboring grid units of the grid unit need to be considered during the calculation process. There are five situations for the calculation of internal points and edge points, ① the points in the five grid cells belong to the same roof layer, ② there are adjacent grid cells with different roof layers on the left or right of the central grid cell, ③ there are Neighborhood grid units of different roof layers are above or below the central grid unit, ④ the horizontal position of the edge points is on the boundary of the central grid unit, ⑤ there are more than two different roof layers in the neighborhood grid units .

c3)、不同屋顶层边缘点优化——边缘点的优化采用主成分分析方法。首先，通过原始建筑物屋顶点云边界点计算得到建筑物的主方向，然后通过迭代的方法将边缘点拟合到主方向对应的直线上。最终得到的建筑物屋顶层的内部点和边缘点如图6所示。c3) Optimization of edge points of different roof layers - the optimization of edge points adopts the method of principal component analysis. First, the main direction of the building is obtained by calculating the boundary points of the original building roof point cloud, and then the edge points are fitted to the straight line corresponding to the main direction by an iterative method. The interior points and edge points of the final building roof layer are shown in Fig. 6.

针对c2)中五种不同情况的具体步骤如下：The specific steps for the five different situations in c2) are as follows:

d1)、针对情况①(图5a)——5个格网单元中的点都属于同一个屋顶层，说明不存在墙面。将中心格网单元的中心确定为中心格网单元内部点的x，y坐标，该内部点的高度为中心格网单元中LiDAR点的平均高度。d1), for the case ① (Figure 5a)—the points in the 5 grid cells all belong to the same roof layer, indicating that there is no wall. Determine the center of the central grid cell as the x, y coordinates of a point inside the central grid cell whose height is the average height of the LiDAR points in the central grid cell.

d2)、针对情况②(图5b)——存在某一个邻域格网单元中的LiDAR点属于两个不同的屋顶层，且邻域格网单元所处的位置在中心格网单元的左边或右边。设两个不同屋顶层为屋顶层j和屋顶层k。首先，需要在二维平面上使用支持向量机算法计算得到两个屋顶层的分割线。然后，以计算得到的分割线与中心格网单元的水平中线的交点作为两个屋顶层边缘点的x，y坐标(x₀，y₀)，则(x₀，y₀，z_j)为屋顶层j的边缘点，z_j为屋顶层j的LiDAR点的平均高度，则(x₀，y₀，z_k)为屋顶层k的边缘点，z_k为屋顶层k的LiDAR点的平均高度。d2), for the case ② (Figure 5b) - there is a LiDAR point in a certain neighborhood grid unit belonging to two different roof layers, and the location of the neighborhood grid unit is on the left side of the central grid unit or right. Let two different roof layers be roof layer j and roof layer k. First, it is necessary to use the support vector machine algorithm to calculate the dividing line of the two roof layers on the two-dimensional plane. Then, take the intersection of the calculated dividing line and the horizontal midline of the central grid unit as the x, y coordinates (x ₀ , y ₀ ) of the edge points of the two roof layers, then (x ₀ , y ₀ , z _j ) is The edge point of roof layer j, z _j is the average height of LiDAR points of roof layer j, then (x ₀ , y ₀ , z _k ) is the edge point of roof layer k, z _k is the average height of LiDAR points of roof layer k high.

d3)、针对情况③(图5c)——存在某一个邻域格网单元中的LiDAR点属于两个不同的屋顶层，且邻域格网单元所处的位置在中心格网单元的上部或下部。设两个不同屋顶层为屋顶层j和屋顶层k。同理，首先需要在二维平面上使用支持向量机算法计算得到两个屋顶层的分割线。然后，以计算得到的分割线与中心格网单元的竖直中线的交点作为两个屋顶层边缘点的x，y坐标。各屋顶层边缘点的高度为该格网单元中对应屋顶层LiDAR点的平均高度。d3), for the case ③ (Figure 5c) - there is a LiDAR point in a certain neighborhood grid unit that belongs to two different roof layers, and the location of the neighborhood grid unit is in the upper part of the central grid unit or lower part. Let two different roof layers be roof layer j and roof layer k. In the same way, it is first necessary to use the support vector machine algorithm to calculate the dividing line of the two roof layers on the two-dimensional plane. Then, the intersection of the calculated dividing line and the vertical midline of the central grid unit is used as the x, y coordinates of the edge points of the two roof layers. The height of each roof layer edge point is the average height of the corresponding roof layer LiDAR points in the grid unit.

d4)、针对情况④(图5d)——中心格网单元中的LiDAR点与邻域格网单元中的LiDAR点恰好处于两个不同的屋顶层。设两个不同屋顶层为屋顶层j和屋顶层k，计算中心格网单元和所述邻域格网单元之间的格网边线的中点坐标(x₁，y₁)，则(x₁，y₁，z_j)为屋顶层j的边缘点，z_j为屋顶层j的LiDAR点的平均高度，则(x₁，y₁，z_k)为屋顶层k的边缘点，z_k为屋顶层k的LiDAR点的平均高度。d4), for case ④ (Fig. 5d) - the LiDAR points in the center grid cell and the LiDAR points in the neighborhood grid cells are exactly in two different roof layers. Let two different roof layers be roof layer j and roof layer k, calculate the midpoint coordinates (x ₁ , y ₁ ) of the grid edge between the central grid unit and the neighborhood grid unit, then (x ₁ , y ₁ , z _j ) is the edge point of roof layer j, z _j is the average height of LiDAR points of roof layer j, then (x ₁ , y ₁ , z _k ) is the edge point of roof layer k, z _k is Average height of LiDAR points for roof layer k.

d5)、针对情况⑤(图5e)——当某一个邻域格网单元中的LiDAR点属于两个以上的屋顶层，以LiDAR点数最多的两个屋顶层确定不同屋顶层的边缘点的x，y坐标，具体情况可参考d2)或d3)。各屋顶层边缘点的高度为该格网单元中对应屋顶层LiDAR点的平均高度。d5), for the situation ⑤ (Figure 5e) - when the LiDAR points in a certain neighborhood grid cell belong to more than two roof layers, the x of the edge points of different roof layers are determined by the two roof layers with the largest number of LiDAR points , y-coordinate, please refer to d2) or d3) for details. The height of each roof layer edge point is the average height of the corresponding roof layer LiDAR points in the grid unit.

第四步、建筑物模型重建——模型的重建包括建筑物屋顶面的重建和墙面的重建。The fourth step, building model reconstruction - the reconstruction of the model includes the reconstruction of the roof surface of the building and the reconstruction of the wall surface.

针对建筑物模型重建，具体步骤如下：For building model reconstruction, the specific steps are as follows:

e1)、屋顶面重建——对属于同一屋顶层的内部点和边缘点构建三角网，形成屋顶面。e1) Roof surface reconstruction—construct a triangular network for the interior points and edge points belonging to the same roof layer to form a roof surface.

e2)、墙面重建——对属于相邻屋顶层且位于同一竖直平面上的边缘点构建三角网，形成屋顶面建筑物墙面，最终完成建筑物三维模型重建。e2), wall reconstruction - construct a triangular network for the edge points belonging to adjacent roof layers and located on the same vertical plane to form the roof surface of the building wall, and finally complete the reconstruction of the three-dimensional model of the building.

本实施例以建筑物屋顶航空LiDAR数据为真实值。以整个区域和6栋建筑物来计算模型与真实点之间的偏移距离来对重建模型的精度进行评估。整个区域模型如图7所示，6栋建筑物模型如图8所示。In this embodiment, the aerial LiDAR data on the roof of the building is used as the real value. The accuracy of the reconstructed model was evaluated by calculating the offset distance between the model and the real point with the whole area and 6 buildings. The model of the whole area is shown in Figure 7, and the model of the six buildings is shown in Figure 8.

表1 建筑物模型评价结果Table 1 Evaluation results of building model

从表1的统计结果来看，从单独的6栋建筑物模型和整个实验区的建筑物模型与真实值的偏移平均距离来看，大多处于0.04m左右，且小于0.3m的百分比超过96％，进一步说明本发明对建筑物屋顶的三维模型重建具有较高精度。From the statistical results in Table 1, from the perspective of the average offset distance between the individual 6 building models and the building models of the entire experimental area and the real value, most of them are around 0.04m, and the percentage of less than 0.3m exceeds 96 %, which further illustrates that the present invention has relatively high precision for reconstruction of the three-dimensional model of the roof of the building.

除上述实施例外，本发明还可以有其他实施方式。凡采用等同替换或等效变换形成的技术方案，均落在本发明要求的保护范围。In addition to the above-mentioned embodiments, the present invention can also have other implementations. All technical solutions formed by equivalent replacement or equivalent transformation fall within the scope of protection required by the present invention.

Claims

1. A method for automatic reconstruction of a three-dimensional building model based on aerial LiDAR data, the steps are as follows:

The first step, building roof point cloud extraction - extract the building point cloud, remove the point cloud on the building wall, and get the building roof point cloud;

The second step is to segment the roof patch—to segment the roof patch of the building roof point cloud, and then level the obtained roof patch;

The third step is to merge the roof patches - when two roof patches are adjacent and there is an intersection line on the plane, the two roof patches are merged to form a roof layer;

The fourth step, roof layer resampling - resample all roof layers to obtain the interior points and edge points of the roof layer;

The fifth step, building model reconstruction—construct a triangular network for the internal points and edge points belonging to the same roof layer to form a roof surface; construct a triangular network for edge points belonging to adjacent roof layers and located on the same vertical plane to form a The roof surface and the wall surface of the building are finally completed to reconstruct the 3D model of the building.

2. a kind of automatic reconstruction method of three-dimensional building model based on aerial LiDAR data according to claim 1, it is characterized in that: in the first step, the method of rejecting the point cloud on the building wall is as follows: take r as Radius, search for neighbor points for each LiDAR point, divide the number of searched neighbor points by the volume of the neighborhood area to get the point cloud density of each LiDAR point, when the point cloud density of LiDAR points is less than the specified threshold, The LiDAR point is a wall point, which is removed.

3. a kind of automatic reconstruction method of three-dimensional building model based on aerial LiDAR data according to claim 1, it is characterized in that: in described second step, select seed region, use region growing algorithm to realize roof patch segmentation, seed The selection method of the area is as follows: first estimate the normal vector of each LiDAR point, take the calculation of point p as an example, find the point set N _p within the radius r with the point p as the center, and calculate according to formulas (1) and (2) Three eigenvalues λ ₁ , λ ₂ , λ ₃ are obtained, and the eigenvector corresponding to the smallest eigenvalue λ _min among the three eigenvalues is the normal vector of point p.

{C C}_{p p} = = \frac{11}{| | {N N}_{P P} | |} {Σ Σ}_{i i = = 11}^{n no} (({q q}_{i i} - - \overset{&OverBar; &OverBar;}{p p})) {(({q q}_{i i} - - \overset{&OverBar; &OverBar;}{p p}))}^{T T} - - - - - - ((11))

\overset{&OverBar; &OverBar;}{p p} = = \frac{11}{n no} {Σ Σ}_{i i = = 11}^{n no} {q q}_{i i} - - - - - - ((22))

where q _i ∈ N _p , n is the number of points in the point set N _p , C _p is the covariance matrix of point p, according to the formula Calculate the curvature of point p, when the curvature of point p When it is less than the specified curvature threshold λ _T , it can be considered that the neighborhood point N _p of point p is on a plane, and the curvature The corresponding neighborhood point N _p is used as the seed area that meets the requirements.

4. a kind of method for automatically reconstructing a three-dimensional building model based on aerial LiDAR data according to claim 3, is characterized in that: in the described second step, region growing algorithm needs to determine two standards: one is the number of internal points number, and the other is the standard deviation of the fitting plane; when the distance from the LiDAR point to the plane where the seed area is located is less than the specified distance threshold, it is considered as an intra-office point; after all points are calculated, the standard deviation of the intra-office point is calculated again, when the standard deviation When it is less than the specified threshold, it is considered to meet the roof patch segmentation requirements.

5. A kind of automatic reconstruction method of three-dimensional building model based on aerial LiDAR data according to claim 4, it is characterized in that: in the second step, the optimization of patch smoothing comprises segmentation patch point cloud smoothing and local prominent points Cloud leveling:

For the leveling of the segmented patch point cloud, the least square method is used to fit the segmented patch point cloud equation, and then the points are leveled to the fitted patch according to the equation, and the fitting plane equation is calculated according to formula (3):

|\begin{matrix} Σ Σ {x x}_{i i}^{22} & Σ Σ {x x}_{i i} {y the y}_{i i} & Σ Σ {x x}_{i i} \\ Σ Σ {x x}_{i i} {y the y}_{i i} & Σ Σ {y the y}_{i i}^{22} & Σ Σ {y the y}_{i i} \\ Σ Σ {x x}_{i i} & Σ Σ {y the y}_{i i} & n no \end{matrix}| (\begin{matrix} {a a}_{00} \\ {a a}_{11} \\ {a a}_{22} \end{matrix}) = = (\begin{matrix} Σ Σ {x x}_{i i} {z z}_{i i} \\ Σ Σ {y the y}_{i i} {z z}_{i i} \\ Σ Σ {z z}_{i i} \end{matrix}) - - - - - - ((33))

In the formula, x _i , y _i , _zi represent the three-dimensional coordinates of the point, a ₀ , a ₁ , a ₂ represent the coefficients of the fitting plane equation;

For the smoothing of the local prominent point cloud, the distance threshold is used to search the local prominent point cloud, and the distance from the point not in the segmented patch to the segmented patch is calculated. When the distance is less than the threshold, it is recorded as the required point, and the required point is According to the fitting plane equation obtained above, it is flattened onto the split surface.

6. A kind of automatic reconstruction method of three-dimensional building model based on aerial LiDAR data according to claim 1, it is characterized in that: in the described 4th step, roof layer resampling method is as follows: first point cloud is projected to XY plane , and build a two-dimensional grid on the XY plane, and then make the following judgments:

(a) When the LiDAR points in the four-neighborhood grid unit and the LiDAR point in the central grid unit belong to the same roof layer, the center of the central grid unit is used as the internal point of the roof layer, and the height of the internal point is Average height of the LiDAR points in the central grid cell;

(b) When the LiDAR points in one of the four neighborhood grid units of the central grid unit belong to two or more roof layers, set the largest number of points in the neighborhood grid unit The two roof layers are roof layer j and roof layer k respectively, and the dividing line between the LiDAR point of roof layer j and the LiDAR point of roof layer k in the neighborhood grid unit is obtained, and then the dividing line and the center line of the central grid unit are calculated The intersection point coordinates (x ₀ , y ₀ ) in the neighborhood grid unit, then (x ₀ , y ₀ , z _j ) is the edge point of the roof layer j, and z _j is the average height of the LiDAR points of the roof layer j , then (x ₀ , y ₀ , z _k ) is the edge point of the roof layer k, and z _k is the average height of the LiDAR points of the roof layer k;

(c) When the LiDAR points in the central grid unit and the LiDAR points in one of the four neighborhood grid units belong to two different roof layers respectively, set the two different roof layers as the roof layer j and roof layer k, calculate the midpoint coordinates (x ₁ , y ₁ ) of the grid edge between the central grid unit and the neighborhood grid unit, then (x ₁ , y ₁ , z _j ) is the roof The edge point of layer j, z _j is the average height of LiDAR points of roof layer j, then (x ₁ , y ₁ , z _k ) is the edge point of roof layer k, z _k is the average height of LiDAR points of roof layer k .

7. A method for automatically reconstructing a three-dimensional building model based on aerial LiDAR data according to claim 6, wherein the size of the two-dimensional grid is twice the average spacing of the LiDAR point cloud.

8. A kind of automatic reconstruction method of three-dimensional building model based on aerial LiDAR data according to claim 6, it is characterized in that: in the described 4th step, obtain the LiDAR point and the roof of roof layer j in the neighborhood grid unit The method of dividing the LiDAR points of layer k is as follows: consider the LiDAR points of roof layer j and the LiDAR points of roof layer k in the neighborhood grid unit as two categories, and use the support vector machine algorithm to obtain the dividing line.