一种基于多源遥感数据的建筑物智能化三维测图方法与流程
本发明属于建筑物三维测图领域，特别是涉及一种基于多源遥感数据的建筑物智能化三维测图方法。
背景技术：
传统立体测图方法的生产效率低下，依靠人工干预实现半自动的建筑物轮廓信息提取需要耗费较大的人工和时间成本。现有的大范围建筑物三维建模方法大多基于大场景点云信息进行构网重建，计算复杂度高，生成的模型不具有结构化的信息，且想要获得单体化的建筑物三维模型比较困难。
公开号为cn109685891a的中国专利申请了一种基于深度图像的建筑物三维建模与虚拟场景生成系统，该专利实现方法包括步骤：计算并重建精确的三维深度数据；对多种类型来源的三维数据进行优化和转换,得到对应的高速流式点云格式的模型,进行点云配准；图像的投影变形,建立各相邻画面像素的对应关系；将点云三角化形成的三角面片连接起来,生成建筑物体表面；使用二次曲面和复杂多边形工具建模,改变二次曲面的状态,针对不规则多边形进行处理,得到建筑物三维模型。该方法直接对点云数据进行构网、二次曲面规则化，适用于单建筑物或者较小范围建筑物的模型重建，此外计算量大且生成的建筑物模型不具备单体化、结构化信息，应用场景较小。
公开号为cn105354883b的中国专利申请了一种基于点云的3dsmax快速精细三维建模方法及系统，该专利实现方法包括步骤：点云预处理,包括先对机载或车载激光雷达所获取的点云数据进行点云去噪,然后进行快速聚类和构建kd树,再进行点云抽稀；基于3dsmax根据聚类后的每个类别的中心点提供全部点云的缩略图；生成局部点云,当用户选取建筑物角点时,通过所在类别的kd树获取邻域点云,采用角点提取算法从邻域点云中提取到其中所有的角点作为特征点,对邻域点云中点云数据进行平面拟合,得到平面的交线,作为特征线；根据特征点修正用户测点,当用户进行测线时,搜索绘制建筑的轮廓边线,修正用户测线；基于3dsmax根据建筑物轮廓完成建模。该专利涉及的方法使用三维点云数据为单一数据来源，且需要人工选取点云的特征角点，所需的人力成本和时间成本较大，结合现有的软件作为处理工具，限制条件较多。
文学东等人在《一种融合多源特征的建筑物三维模型重建方法》中提出利用结合机载点云,地面点云及倾斜多视纹理实现建筑物面元的拓扑重建，实现建筑物的三维建模。这种方法的缺点是需要人工的交互编辑实现精确的线特征约束，且对数据来源要求较多。
基于航空遥感影像数据具有边界清晰，分辨率高且细节信息丰富的特点；再结合激光点云具有的精度高，穿透力强，具备三维坐标和强度、颜色等信息，实现多源遥感数据的优势互补，为建筑物智能三维测图带来了新的实现方式。
技术实现要素：
本发明的目的是提供一种基于多源遥感数据的建筑物智能化三维测图方法。
本发明区别于传统的航空遥感影像数据密集匹配生成密集点云、三维构网、纹理映射的方式，本发明包含一种具有快速且智能化的建筑物屋顶二维矢量轮廓提取、结合oesm生成建筑物结构化线框模型等功能的方法及系统。采用本发明方法，可以大大节省人工成本的同时满足一定比例尺的三维测图的精度要求，同时保留了建筑物结构线之间的矢量拓扑关系，便于后续的统计、分析与实际应用，具有巨大的应用前景和拓展空间。
为达到上述目的，本发明提供基于多源遥感数据的建筑物智能化三维测图方法及系统，包括：
步骤1：利用目标测区的航空遥感影像数据和高精度pos数据进行影像高精度定位定姿、密集匹配生成高精度影像密集点云；
步骤2：将步骤1所述的高精度影像密集点云作为参考，对完成自配准和拼接后的激光点云进行配准，生成配准融合后的三维点云；
步骤3：利用配准融合后的三维点云，进行插值、栅格化得到高精度的dsm数据；对高精度的dsm数据进行滤波得到对应的dem数据；
步骤4：结合目标测区的航空遥感影像数据及其定向参数、步骤3中得到的高精度dsm数据、dem数据，生成传统dom影像和oesm模型；
步骤5：利用改进的8方向sobel边缘检测算子结合canny算子对步骤4得到的传统dom进行边缘提取，然后利用hough变换检测存在于边缘中的线特征，得到传统dom中的初始线特征集合；
步骤6：对步骤2中所述配准融合后的三维点云进行点云面片分割和建筑物屋顶点云顺序轮廓提取，得到建筑物屋顶点云的顺序轮廓；
步骤7：结合步骤6中得到的建筑物屋顶点云顺序轮廓，在步骤4中生成的dom上提取建筑物屋顶二维矢量轮廓；
步骤8：闭合步骤7得到的补全之后的顺序线特征集合得到目标测区的建筑物屋顶二维矢量轮廓，基于建筑物屋顶二维矢量轮廓、oesm、dem进行建筑物智能三维测图，生成结构化的建筑物三维线框模型；
进一步的，步骤4生成传统dom影像具体为：
对于dem数据上的每一点(xdem,ydem,zdem)，通过与对应影像的高精度pos数据建立共线条件方程解算其对应的影像坐标(xe,ye)，将该影像坐标对应的灰度值data(x,y)作为传统dom在该点的灰度值，当所有的dem单元完成共线方程解算后即可得到完整的传统dom；根据共线方程解算dem单元对应的影像坐标的公式如下：
式中，xe，ye表示dem某单元对应的影像平面坐标，xdem，ydem，zdem为dem上一点的坐标，其他字母表示影像的内外方位元素。
进一步的，步骤4生成oesm模型具体为：
对于dsm模型上的每一个点(xdsm,ydsm,zdsm)，通过与该栅格位置对应的航空遥感影像数据及高精度pos数据建立共线条件方程，反算得到该dsm点对应的影像坐标(xs,ys)；
然后根据该影像坐标(xs,ys)求取对应的dem点坐标(xdem,ydem,zdem)，具体如下：
步骤4.1,先假定该待求点的dem高程是整个dem模型的最大高程zmax，结合已知的影像坐标(xs,ys)，利用影像与dem之间建立的共线条件方程就可以解算该假设高程对应的dem平面坐标的计算值(x'dem,y'dem)；
步骤4.2,根据坐标计算值(x'dem,y'dem)查询dem上的真实高程值z'dem，计算前后两个高程值zmax与z'dem之间的差是否小于阈值；
步骤4.3,如果小于阈值则认为找到了真实的dem坐标，结束计算，否则进入步骤4.4；
步骤4.4,如果不满足阈值要求，将z'dem作为新的dem高程假设值，取代第一次设置的最大高程值zmax，再按照步骤4.1至步骤4.3所述步骤重新计算新的dem平面坐标(x”dem,y”dem)以及该坐标对应的dem真实高程z”dem；
直至高程假设值与真实值之间的差小于阈值则认为找到真实的dem坐标
(xdem,ydem,zdem)；
对于dsm数据上的此点(xdsm,ydsm,zdsm)对应的oesm坐标为(xdem,ydem,zdsm)；
对于落入同一个dem单元格内部的多个dsm高程仅保留高程值最大值作为最终的oesm坐标高程；
将所有的dsm单元完成计算后，对oesm进行插值生成完整的oesm；插值的方式采用反距离插值法，即对oesm上的每一个空洞查找其在x，y方向上最近的非空洞单元格，根据距离的倒数作为权重进行高程的加权求和，将求和的结果作为空洞单元格的最终高程计算值；
进一步的，改进的8方向sobel算子较传统的2方向、4方向对多方向边缘更加敏感；
进一步地，所述步骤6具体包括以下步骤：
步骤6.1，对配准融合后的三维点云进行滤波得到滤波后的三维点云，去除离散点，得到初步滤波后的三维点云；
设置滤波领域半径r、领域内最小点数n，删除领域内点数小于n的点；
步骤6.2，基于ransac原理对滤波后的点云进行去地面，得到去地面后的三维点云；
设置内点到模板平面的距离阈值为d、迭代次数为k、作为平面的最小点数为m，在多次迭代中，选取被保留点数最多的平面作为输出地面面片；
步骤6.2目的是通过ransac方法挑选出存在于三维点云中的地面面片；
所述地面面片的定义为：
ax+by+cz＝d
当平面方程的参数a，b，c，d确定，这个空间中的平面就被唯一确定，对三维点云中所有的点而言，落在这个平面上的点也就唯一确定即点到这个平面的距离小于设定的阈值就认为落在这个平面上，ransac原理就是不断地变更平面的a，b，c参数并计算每次落在平面上的点的数目，在达到设定的迭代次数后，选择点数最多的那一次平面参数作为结果即地面面片；
步骤6.3，对步骤6.2之后的去地面剩余点云进行基于区域生长方法的点云面片分割，得到包含非建筑物屋顶面片在内的所有点云面片；对于去地面后的点云，计算点云在各点的曲率和法向量并选取最小曲率点作为种子点，设置曲率阈值cur和法向量变化阈值nor进行种子点生长；cur、nor为经验值；
所述计算点云在各点的曲率和法向量并选取最小曲率点作为种子点为：
点云某一点的法向量计算：点云中某一点领域内k个点，对这k个点利用最小二乘法进行局部平面拟合得到局部的平面方程ax+by+cz＝d，其法向量即为(a,b,c)；
点云内某一点的曲率计算：以该点为中心取领域内k个点进行最小二乘法的二次曲面拟合得到局部拟合的二次曲面方程z(x,y)＝ax2+bxy+cy2，根据二次曲面方程即可求取对应的曲率，将该曲率作为点云在该点的曲率表达；
步骤6.4，对于步骤6.3面片分割得到的点云面片，进行基于面片特征的建筑物屋顶面片筛选和基于法线估计的建筑物屋顶面片轮廓提取，得到建筑物屋顶的点云边界；
进一步的，步骤6.4中从分割后的点云面片筛选建筑物屋顶面片的规则是：对于每个面片求取其对应的面片平均高程elei；即通过点云面片的平均高程等于点云内所有点的高程之和除以总的点数；
将此平均高程的计算结果与步骤6.2提取到的地面面片的平均高程ele0进行作差比对，差值大于阈值δele的面片保留；
计算每个点云面片的法向量nori，将此法向量的计算结果与步骤6.2提取到的地面面片的法向量nor0进行夹角大小的判断，小于夹角阈值的进行保留；
步骤6.5，对步骤6.4得到的无顺序的屋顶点云轮廓构建平面三角网，根据三角网提取点云轮廓的凹包，获得建筑物屋顶点云的顺序轮廓；
优选地，所述步骤7具体包括以下步骤：
步骤7.1，将步骤6中的建筑物屋顶点云轮廓贴合到步骤4中得到的传统dom上，按照缓冲半径buf_r生成影像缓冲区；
进一步的，生成缓冲区的方法为：步骤6点云处理之后得到了具体的单体建筑物屋顶边界点云，且具有正确的凹包顺序。将这些点映射到dom上(由于激光点云与影像已经完成了配准，参考步骤1中所述)，给定缓冲半径buf_r，一般根据影像分辨率设置为0.2m左右，根据该缓冲半径生成点云边界在dom影像上的缓冲区。缓冲区的生成规则如下：
min_dis(col,row)＝min(dis(p(col,row),pk)),k∈(0,bry_poi_sum)
0＜min_dis(col,row)＜buf_r,col∈(0,col_sum),row∈(0,row_sum)
其中，dis(p(col,row),pk)表示dom上一点p(col,row)与点云边界上第k点pk之间的距离，bry_poi_sum为点云边界总点数，min_dis(col,row)即为dom上第row行，第col列的点p(col,row)到点云边界的最小距离。col_sum，row_sum为dom总列，行数。
步骤7.2，按照步骤6.1得到的缓冲区对步骤4得到的线特征进行过滤，仅保留存在于缓冲区内部的线特征；
进一步的，步骤7.2中对处于缓冲区内部的线特征进行保留，判断在缓冲区内部的定义为：
判断该线特征上的每一个像素是否在缓冲区内，总计所有落在缓冲区内部的像素总数，最终计算这部分像素占整个线特征总像素数的比例是否大于阈值；
步骤7.3，对步骤7.2被保留的线特征按照方向向量和相对位置判断相似性，按照相似性对线特征进行分类，得到分类后的线特征集合；
进一步的，判断两线段特征相似的规则为：
方向——两线段之间的斜率变化或者线段的方向向量之间的夹角小于阈值(一般设置为5°)。
式中，line1,line2表示参与分类判断的两条线段。
位置——两线段如果满足方向上的相似条件，则判断这两条直线的位置关系。两线段中点连线在垂直于两线段方向向量方向上的投影长度小于阈值(设置为某固定值，根据精确性要求在图像上可设置为5-10个像素)，且两线段中点连线在平行于两线段方向向量方向上的投影长度小于阈值(设置为两线段长度和的1/2加上某固定值)。
m1m2·nt＜t_threshould
式中，m1m2表示参与分类的两线段中点组成的线段，nl，nt表示两线段单位方向向量和单位法向量，length_line1，length_line2表示线段line1，line2的长度，l_threshould，t_threshould表示两方向上的检测阈值。
步骤7.4，将步骤6得到的建筑物屋顶点云顺序轮廓按照距离阈值t_dis垂直投影到步骤7.3中分类后的线特征上，即对点云轮廓上的每一个点向所有线特征上作垂线，垂足落在线特征上且垂线的长度小于阈值视为投影成功；投影完成后，统计每一条线特征上的垂足数目，仅保留垂足数目不小于t_num的线特征；t_dis、t_num为经验值；这一步的目的是通过点云投影过滤大部分的非建筑物边界的线特征，保留置信度大的线特征集合。
步骤7.5，对步骤7.3、7.4后得到每类线特征进行最优线特征的挑选，即对每一类线特征仅保留作为建筑物边界置信度最大的线特征，挑选指标为：线特征长度、投影距离分布标准差、投影密度；指标权重为经验值；
进一步的，对步骤7.3、7.4后得到的每一类线段特征进行最优挑选的规则是：
线段长度，作为建筑物某一条真边界的一类线段中，长度越长其作为建筑物边界被保留的置信度越大；
各垂足与原始点云中的对应点之间的距离即投影距离的分布标准差，建筑物真边界理论上各点的投影距离一致，标准差为0；
垂足离散度即单位长度的线段所包含的垂足数量，越多则该边作为真边界的置信度越大。
按照以上三个规则对每一条候选线段进行打分，如下所示：
式中，scor1，scor2，scor3分别对应1-3规则所得分数，
步骤7.6，对步骤7.5得到的最优线特征集合按照步骤6.5中的建筑物屋顶点云顺序边界进行排序和具有独立方向的线特征检测，得到顺序的且不具有突变方向的dom线特征集合；
进一步的，检测具有独立方向的线特征的目的是去除那些具有突变方向且长度较短的线段特征，具体的检测步骤为：
按照斜率对线段进行分类。
对于无法归类的单一线段判断其长度。
如果线段长度小于2*line_length_threshould(每一条线特征应当具有的最小长度，为经验值)，直接舍弃。
步骤7.7，初次闭合步骤7.6得到的顺序线特征集合，检测是否存在缺边；检测存在缺边则进行缺边的补齐，得到顺序的且完成缺边补齐后的dom线特征集合；
进一步的，缺边判断规则如下：
第一缺边规则：前后线段间相邻且平行。
第二缺边规则：线段相邻且前一线段末尾点与后一线段起始点，即对应顺序点云集的点之间空缺点数超过阈值。
进一步的，对于缺边位置，需要进行补全或推理，缺边补全规则如下：
第一补全规则：在原始影像初筛线段中查找通过空缺边的前一线段尾点与后一线段起点的线段。
第二补全规则：第一补全规则查找到多条符合规则线段则按照如下规则筛选最优线段：
斜率属于最优线段集按照斜率分类后的某一类。
长度与空缺边的前一线段尾点与后一线段起点间连线长度接近。
第三补全规则：第一补全规则找不到符合规则的线段。对于属于第一缺边规则的缺边；作过前一线段结尾点与后一线段起点连线中点且垂直于前一线段和后一线段的垂线，两垂足连线作为推理缺失边。对于属于第二缺边规则的缺边；参考顺序点云中的空缺点，对这些空缺点进行dp规则化得到推理空缺边；
进一步的，步骤8中利用建筑物屋顶二维矢量轮廓和oesm结合获取建筑物结构化线框模型的具体做法为：
求二维轮廓的每个角点对应的建筑物高，将每个角点从平面坐标变换为三维坐标，z坐标实际上为建筑物高；
求建筑物高时，先查询角点在oesm上的高程zoesm、在dem上的高程zdem，则该角点对应的建筑物高为zoesm-zdem；
利用每个角点的建筑物高和平面坐标即可计算出整个建筑物的线框三维模型；由于已知建筑物屋顶角点的oesm坐标，该坐标包含dem上的x，y坐标以及与该dem单元保持共线关系的dsm单元的高程作为z值，可以根据如下所示公式计算出角点对应的真实三维坐标即dsm坐标。
其中，xs，ys，zs为摄影中心的坐标。在计算得到建筑物屋顶角点的dsm坐标后，结合建筑物高，可以得到屋顶角点对应的地面点的真实三维坐标，根据建筑物屋顶角点及其对应的地面点的真实三维坐标即可构建建筑物的线框三维模型。
和现有技术相比，本发明具有如下优点和有益效果：
有效的结合了航空遥感影像数据和激光点云这两种不同来源的遥感数据，实现了对目标测区的建筑物屋顶二维矢量轮廓的智能化提取，最终生成了结构化的建筑物线框三维模型。
实现了从轮廓提取到模型生成的智能化和高效化，大大节省了时间成本和人力成本，同时在数据质量得到保障的前提下也满足我国建筑物三维测图的精度需求；
生成的结构化三维模型保留了建筑物作为矢量模型具有的拓扑关系，区别于传统构建三角网丧失了该拓扑关系的缺点，具有十足的应用前景和拓展空间；
附图说明
图1为激光点云与影像配准的流程图。
图2为对激光点云进行屋顶面片分割的原理示意图。
图3为基于法线估计的建筑物屋顶点云面片边缘轮廓提取的原理示意图。
图4为基于deluanay平面三角网的凹包顺序边界提取原理示意图。
图5为改进的8方向sobel影像边缘提取算子示意图。
图6为根据数字表面模型(dsm)，数字高程模型(dem)和影像生成可量测影像高程同步模型(oesm)的流程示意图。
图7为点云边界指导影像缓冲区原理示意图。
图8为利用缓冲区对正射影像中线特征进行初步筛选的原理示意图。
图9为对正射影像线特征进行分类的判断准则示意图。
图10为点云边界按阈值投影到线特征进行再次筛选的示意图。
图11为点云边界按阈值投影到线特征的效果图。
图12为对同类线特征进行最优挑选的规则示意图。
图13为对同类线特征进行最优挑选的结果图。
图14为对最优挑选后线特征进行方向性检测从而去除异常短边的示意图。
图15为对线特征集合进行空缺边判断的判断规则和对应的补全规则示意图。
图16为线特征集合通过空缺边补全后的结果图。
图17为闭合顺序线特征集合得到建筑物闭合多边形边界的结果图。
图18为利用建筑物多边形边界结合oesm生成三维模型的结果图。
图19为本发明方法流程图。
具体实施方式
以下结合附图和实施例详细说明本发明的技术方案。
本发明是基于多源遥感数据的自动化三维测图方法及系统。首先对测区内激光点云与高精度航空遥感影像数据空三、密集匹配生成的影像密集点云的进行配准和融合，保证两者在三维坐标上的一致性。利用配准后的三维点云数据进行高精度插值生成高精度的dsm模型，将dsm进行滤波处理得到目标测区的dem模型，结合dem、dsm模型和原始航空遥感影像数据以及其对应的高精度定向参数生成保留了严格共线关系的传统dom影像和高精度oesm模型。再对配准融合后的三维点云进行点云分割与面片筛选，仅保留测区内的建筑物屋顶点云面片。对得到的建筑物屋顶面片进行凹包提取，获取建筑物点云顺序边界。对生成的传统dom进行边缘检测和线特征检测获取测区内所有线特征。利用建筑物点云顺序边界生成缓冲区对dom上提取的线特征进行筛选，保留包含所有可能是该建筑物边界的且含少量噪声的线特征集合。按照线段梯度与长度及其位置信息对前一步得到的线特征集合进行分类，去除孤立线段。随后按照阈值将建筑物点云顺序边界投影至上一步分类后的各线段，记录各线段与点云边界各点的对应关系，根据一定规则进行线段特征的最优筛选。对最优筛选后的线段特征进行排序和方向性检测。然后进行线段特征的空缺边检测和补齐，闭合线段集合得到闭合的建筑物多边形边界。最后根据得到的建筑物闭合边界结合oesm生成最终的建筑物三维模型。
下面结合图1至图19介绍本发明的具体实施方式。
本发明提供的实施例一种基于多源遥感影像的建筑物自动化三维测图方法，包括以下步骤：
步骤1：基于迭代最近点方法(icp，iterativeclosestpoint)对测区内的激光点云数据与航空遥感影像数据经过空三、密集匹配生成的影像点云进行配准，得到配准后的三维点云数据。
基于icp方法的点云配准目的是求取激光点云p转换到航空遥感影像数据生成的影像点云q的旋转平移矩阵，分别用r(旋转矩阵)，t(平移矩阵)表示；其变换关系可表示为：
pq＝r·pp+t
其中pq和pp表示点云q和p中的一对对应点，其具体实现包括以下子步骤：
步骤1.1，确定初始的最小化迭代最近点的目标函数。如下所示：
其中，
步骤1.2，确定初始的少量对应点集，并确定初始的旋转平移矩阵r0，t0。
步骤1.3，根据初始的旋转平移矩阵确定所有的对应点集，首先利用初始的r0，t0对激光点云p进行变换，得到变换后的点云p′，对比p′与影像点云q，如果存在距离小于一定阈值d的点对，则认为这两个点就是一组对应点。
步骤1.4，对r0，t0的优化。有了全部的对应点集后，利用最小二乘方法求解当前对应点下最优的旋转平移矩阵ri，ti。计算当前旋转平移矩阵下的目标函数值f(ri,ti)是否小于一定阈值。如果不满足条件则利用ri，ti继续参与步骤1.3，1.4进行迭代，直至收敛。
步骤2：ransac和区域生长方法结合的点云面片分割，得到建筑物屋顶的点云面片；
步骤2.1，对初始的三维点云，利用ransac方法去除地面面片g，得到去地面后的三维点云数据。需要设置三个参数，分别为：最大的迭代次数sac_max_iteration，点属于平面内点的最小距离阈值sac_dis_min，平面模型能作为地面面片应当具备的最小点数sac_min_poisnum。算法思路为：在迭代次数不超过最大阈值的情况下，不断变换模型参数计算每个模型下满足最小距离阈值的总内点数(内点数大于sac_min_poisnum保留否则舍弃该模型)，仅保留总内点数最多的模型作为最终输出。
步骤2.2，对去地面后的三维点云数据，利用区域生长方法进行点云面片分割和筛选，得到建筑物屋顶点云面片。需要设置的参数有法线差异阈值rg_nor_threshould，曲率差异阈值rg_curv_threshould，区域生长所得面片应当包含的点数最大与最小阈值rg_max_poisnum和rg_min_poisnum。算法思路为：选取初始种子点，判断种子点与周围一定范围内其余点的法线方向差异是否小于rg_nor_threshould，曲率变化是否小于rg_curv_threshould，均满足条件则加入种子点集。从初始种子点出发直至没有新的种子点产生，完成这一面片的聚类。判断该面片点数是否满足rg_max_poisnum和rg_min_poisnum范围要求，满足则保留，否则舍弃。设区域生长得到的面片集合为srg，判断该面片集合中任意面片srgi的平均高程hrgi与地面面片g的平均高程hg差值是否满足阈值(一般认为屋顶面片的平均高程应当高于地面至少3m)；判断该面片集合中任意面片srgi的法线nrgi与地面面片g的法线ng夹角是否满足阈值(一般认为屋顶面片与地面的夹角应当小于
δhmin＜hrgi-hg＜δhmax
步骤3：对于建筑物屋顶的点云面片，进行基于法线估计的点云边界提取，得到无顺序的建筑物屋顶点云边界。
基于法线估计的点云边界提取需要设置两个参数，分别是：法线估计的搜索半径es_radius，根据法线估计判断边界的角度阈值es_angle(一般设置pi/4)。算法思路为：按照搜索半径es_radius寻找某一点周围的半径内几个点，通过最小二乘结算这几个点组成的平面，以此平面的法向量作为点云在这一点上的法线估计结果。一般来说，设计搜索半径es_radius为点云空间分辨率的10倍效果较好；搜索半径太小噪声较大，搜索半径太大估计速度太慢。根据法线估计结果进行角度阈值判断，法方向变化超过阈值es_angle，说明该点为边界点。
步骤4：基于delaunay平面三角网的凹包算法提取点云顺序边界，得到建筑物屋顶点云顺序边界点集合poi_bry_ord_set。
步骤4.1，假设步骤3得到的点云边界点集为poi_bry_set，首先为点集poi_bry_set求取delaunay平面三角网bry_delau(不考虑各边界点的高程差异，仅考虑平面坐标)。
步骤4.2，为bry_delau初始化所有的边对象(edge)，并求取edge的长度以及邻接三角形集合。其中邻接2个三角形的边为内部边，1个三角形的为边界边，0个三角形的为计算过程中会退化的边。
步骤4.3，将所有长度大于阈值的边界加入队列，并当队列非空时循环下列过程：
(1)从队列中取出一条边，记为e1，得到e1的唯一邻接三角形t1。
(2)从t1中找出另外两条边，记为e2，e3。将他们的邻接三角形集合中删除t1。
(3)将e2，e3中新形成的长度大于阈值的边界加入队列。
(4)将e1置无效标记，若e2，e3有退化的，也置无效标记。
步骤4.4，收集所有的有效的边界边，形成边列表，输出，按照被保留的有效边所包含的有效点得到最终的建筑物屋顶点云的顺序轮廓点集。
步骤5，根据测区dem、影像数据及其高精度pos数据生成测区的传统dom；再利用传统canny算子结合改进的8方向sobel算子获取传统dom上的边缘轮廓，得到传统dom上的二值化边缘信息，融合两种算子的检测结果，增强边缘特征，减少噪声。
根据测区dem、影像数据及其高精度pos数据生成传统dom的具体步骤为，对于dem数据上的每一点(xdem,ydem,zdem)，通过与对应影像的高精度pos数据建立共线条件方程解算其对应的影像坐标(xe,ye)，将该影像坐标对应的灰度值data(x,y)作为传统dom在该点的灰度值，当所有的dem单元完成共线方程解算后即可得到完整的传统dom；根据共线方程解算dem单元对应的影像坐标的公式如下：
式中，xe，ye表示dem某单元对应的影像平面坐标，xdem，ydem，zdem为dem上一点的坐标，其他字母表示影像的内外方位元素。
步骤6，利用步骤5获取到的dom二值化后的边缘信息，利用hough变换检测存在于边缘中的线特征，获得影像线特征集合lines_set。
步骤7，根据测区dsm、dem结合影像数据生成可量测影像高程同步模型(oesm)。oesm保留了dem的平面坐标信息和dsm上与dem一一对应的坐标点对应的高程信息，oesm体现出严密的共线关系，将dem，dsm有机结合起来。
步骤7.1，利用dsm、影像及其高精度pos数据，根据影像内外方位元素建立dsm上各像素与影像上对应像素的共线条件方程，计算得到dsm上某一点对应的影像像素坐标。如下所示：
式中，x，y表示dsm上一点对应的影像平面坐标，xdsm，ydsm，zdsm为dsm上一点的坐标，其他字母表示影像的内外方位元素。
步骤7.2，根据步骤7.1得到的dsm上某一点对应的影像像素坐标，迭代求取dem上对应点的坐标。迭代求取的步骤为：
(1)选择dem上一随机高程dem_h1(建议取测区dem的最大高程)为待求dem点的高程，解算出该随机高程对应的dem点的平面坐标(dem_x1,dem_y1)。如下所示：
(2)根据上一步得到的dem点平面坐标(dem_x1,dem_y1)，查找dem上该点真实的高程dem_h2。
(3)比较高程差是否小于阈值，若差值大于阈值则利用新得到的高程继续求解新的dem点平面坐标，再查找该新的平面坐标对应的dem高程，比较高程差是否满足阈值要求，即迭代1，2，3步骤。
(4)直至满足阈值要求，就认为找到了真实的dem对应点。
对于落入同一个dem单元格内部的多个dsm高程仅保留高程值最大值作为最终的oesm坐标高程；
将所有的dsm单元完成计算后，对oesm进行插值生成完整的oesm；插值的方式采用反距离插值法，即对oesm上的每一个空洞查找其在x，y方向上最近的非空洞单元格，根据距离的倒数作为权重进行高程的加权求和，将求和的结果作为空洞单元格的最终高程计算值；
步骤8，利用步骤3，4得到的建筑物屋顶点云顺序边界指导生成dom影像缓冲区。点云处理之后得到了具体的单体建筑物屋顶边界点云，且具有正确的凹包顺序。将这些点映射到dom上(由于激光点云与影像已经完成了配准，参考步骤1中所述)，给定缓冲半径buffer_r，一般根据影像分辨率设置为0.2m左右，根据该缓冲半径生成点云边界在dom影像上的缓冲区。缓冲区的生成规则如下：
min_dis(col,row)＝min(dis(p(col,row),pk)),k∈(0,bry_poi_sum)
0＜min_dis(col,row)＜buffer_r,col∈(0,col_sum),row∈(0,row_sum)
其中，dis(p(col,row),pk)表示dom上一点p(col,row)与点云边界上第k点pk之间的距离，bry_poi_sum为点云边界总点数，min_dis(col,row)即为dom上第row行，第col列的点p(col,row)到点云边界的最小距离。col_sum，row_sum为dom总列，行数。
步骤9，利用dom影像缓冲区对步骤6中得到的所有线特征进行筛选，保留落在缓冲区内的线特征作为研究对象，得到筛选后的线特征集合lines_inbuf_set。筛选的规则如下：
式中，line表示lines_set中的某一具体线特征，within_line表示line落在缓冲区内的部分。
步骤10，对步骤9得到的线特征集合lines_inbuf_set进行分类，去除孤立的(所属类别包含线特征数过少)，长度小于阈值line_length_threshould(作为一条具有边界特征的线段应当具有的最小长度)的线段，得到分类后的传统dom线特征集合。分类规则为：
(1)方向——两线段之间的斜率变化或者线段的方向向量之间的夹角小于阈值(一般设置为5°)。
式中，line1,line2表示参与分类判断的两条线段。
(2)位置——两线段如果满足方向上的相似条件，则判断这两条直线的位置关系。两线段中点连线在垂直于两线段方向向量方向上的投影长度小于阈值(设置为某固定值，根据精确性要求在图像上可设置为5-10个像素)，且两线段中点连线在平行于两线段方向向量方向上的投影长度小于阈值(设置为两线段长度和的1/2加上某固定值)。
m1m2·nt＜t_threshould
式中，m1m2表示参与分类的两线段中点组成的线段，nl，nt表示两线段单位方向向量和单位法向量，length_line1，length_line2表示线段line1，line2的长度，l_threshould，t_threshould表示两方向上的检测阈值。
步骤11，将建筑物屋顶点云顺序边界点集合poi_bry_ord_set中的各点投影到步骤10分类后的各线特征上，保留投影距离小于阈值的投影垂足；完成投影后，记录每条线段上包含的垂足数量，保留垂足数不小于3的线段，其余舍弃。(一条线段作为有效的建筑物边界至少包含不少于3个边界点)
步骤12，完成点云边界到影像线段上的投影后，进行每一类线段(步骤10完成了线段的分类)的最优选择，得到每一类线特征中作为建筑物边界置信度最大的线特征。挑选规则如下：
(1)线段长度(作为建筑物某一条真边界的一类线段中，长度越长其作为建筑物边界被保留的置信度越大)。
(2)各垂足与原始点云中的对应点之间的距离(投影距离)的分布标准差(建筑物真边界理论上各点的投影距离一致，标准差为0)。
(3)垂足离散度即单位长度的线段所包含的垂足数量(越多则该边作为真边界的置信度越大)。
按照以上三个规则对每一条候选线段进行打分，如下所示：
式中，scor1，scor2，scor3分别对应1-3规则所得分数，
步骤13，对步骤12得到的最优筛选后的线特征进行排序和异常方向检测，得到具有顺序的且不具有突变方向的线特征集合。步骤12得到的线段集合不具备先后性，无法得知建筑物边界的拓扑关系，因此需要按照点云边界poi_bry_ord_set对各线段进行头尾判断和排序(参考步骤11的投影垂足顺序即可完成排序)。异常方向检测，目的是排除具有突变方向(与建筑物主方向不一致且不垂直于建筑物主方向)的较短线段。步骤如下：
(1)按照斜率对线段进行分类。
(2)对于无法归类的单一线段判断其长度。
(3)如果线段长度小于2*line_length_threshould，直接舍弃。
步骤14，对步骤13得到的线特征集合进行空缺边检测与补齐，得到补齐后的线特征集合。对于步骤13得到的具有顺序的最优边集合(完成了异常方向检测排除了错误短边)进行闭合性与完整性检测，判断是否存在缺边。判断规则：
(1)相邻线段平行(方向向量之间的夹角小于阈值)则认为存在缺边。
(2)前后线段之间存在6个点(参考步骤11的点云到线段的投影)以上的缺失。
对应的补全规则：
(1)记前一线段结尾点与后一线段起点连线的中点为点pmid，作直线过点pmid且垂直于前后线段，垂足为pf1，pf2，线段pf1pf2作为缺边的推理。
(2)对前后线段之间缺失的点使用道格拉斯-普朗克(dp)算法规则化得到缺边。
步骤15，按照顺序闭合步骤14得到的线段集合，得到最终的建筑物二维矢量轮廓。
步骤16，按照建筑物多边形的拐点坐标信息，查询对应坐标在oesm上的高程zoesm，在dem上的高程zdem，得到建筑物屋顶面在该拐点的建筑物高δz。且有：
δz＝zoesm-zdem
步骤17，根据建筑物屋顶面各个拐点的位置坐标和建筑物高，生成该建筑物的三维模型。由于已知建筑物屋顶角点的oesm坐标，该坐标包含dem上的x，y坐标以及与该dem单元保持共线关系的dsm单元的高程作为z值，可以根据如下所示公式计算出角点对应的真实三维坐标即dsm坐标。
其中，xs，ys，zs为摄影中心的坐标。在计算得到建筑物屋顶角点的dsm坐标后，结合建筑物高，可以得到屋顶角点对应的地面点的真实三维坐标，根据建筑物屋顶角点及其对应的地面点的真实三维坐标即可构建建筑物的线框三维模型。
以下提供本发明实施例所提供方法的应用实例具体过程来说明本发明的技术效果。
1.影像点云与激光点云的高精度配准融合
首先利用航空遥感影像数据和高精度的pos数据经过影像高精度空中三角测量、密集匹配生成高精度的影像密集点云数据，然后完成激光点云与影像点云的高精度配准。通过激光点云与影像密集点云的自动粗匹配，以及基于icp算法的精匹配方法(即利用icp迭代求取变换矩阵的方式)，实现激光点云与影像密集点云的高精度配准。如图1所示，最终的均方根误差控制在1.0像素以内，保证测区内的点云坐标信息准确性。
2.对配准融合后点云的处理
利用配准后的三维点云数据进行高精度插值生成高精度的dsm模型，将dsm进行滤波处理得到目标测区的dem模型，结合dem、dsm模型和原始影像以及其对应的高精度空中三角测量数据生成保留了严格共线关系的传统dom影像。
对配准融合后的点云进行面片分割。点云面片分割的目的是从初始点云中提取出建筑物屋顶点云面片。采用ransac结合区域生长的方式去除地面和墙体仅保留高程较高(局部区域内)、法向量变化不大且大概率平行或接近平行于地面的屋顶面片。如图2所示，非屋顶面片(例如墙面)不满足与地面法方向夹角要求和平均高程阈值要求(图中nor3，h3)。建筑物屋顶信息被完整保留，小的噪声信息分布离散，不会干扰单体建筑物的完整性。
基于法线估计的点云边界提取。点云面片分割完成后得到了初始的屋顶面片，首先从屋顶点云上计算出法线，再由法线结合数据估计出边界。如图3所示。首先进行法线求解(平面的法线是垂直于它的单位向量。在点云的表面的法线被定义为垂直于与点云表面相切的平面的向量)。对点云数据集的每个点的法线估计，可以看作是对表面法线的近似推断。因此该表面的判断就是你寻找的周围几个点或者半径内几个点组成的平面，该参数的设置一般设置为分辨率的10倍时，效果较好，主要是对于法线估计。邻域半径选择太小了，噪声较大，估计的法线就容易出错，而搜索邻域半径设置的太大估计速度就比较慢。
凹包算法求点云顺序边界。上一步得到的点云边界不是顺序点集，不具备多边形边界的拓扑关系，这里采用构建平面delaunay三角网的方式确定点云的顺序边界(这里不考虑点云的高程，只求取二维凹包顺序边界)。因为想要求取的边界就是delaunay三角网的子集，所以只需要从delaunay三角网的最外层边开始外不断的删去超过长度限制的边，当这个过程结束时，就能够得到一个大致符合预期的边界。如图4所示。
3.对生成的传统dom的处理
输入利用原始影像和dem生成的传统dom，采用改良sobel算子与canny算子结合获取影像的边缘轮廓。改良sobel算子8方向具有更高的检测灵敏度。算子结构如图5所示。两种算子结合对测区dom进行扫描。
dom线段特征提取。利用上一步生成的影像二值化边缘，利用hough变换检测存在于影像边缘中的线特征。
4.可量测影像同步高程模型(oesm)的生成
根据dsm、dem和影像数据生成oesm，结合了传统dom的坐标和dsm的高程，保留了严密的共线关系。如图6所示。
(1)根据共线条件方程确定dsm每个像素对应的影像像素坐标。
(2)根据步骤1中的得到的影像像素坐标迭代估计dem上对应的像素位置。迭代求取的步骤为：
a.选择一随机高程为待求dem点的高程，结算处dem点的xy
坐标。
b.查找该xy坐标对应的真实dem高程，计算两高程之差。
c.比较高程差是否小于阈值，若大于则利用新得到高程继续迭代。
d.直至满足阈值要求结束求解，对于落入同一dem单元格内的多个高程仅保留最大值作为最终的oesm高程。
(3)对于生成的oesm进行插值，弥补孔洞；插值采用反距离插值法，即对于oesm上的每一个空洞，查找其沿x，y方向上距离最近的非空洞单元格(一般情况下有4个)，根据到非空洞单元格的距离的倒数进行高程的加权求和，将求和的结果作为空洞单元的最终高程。
5.多源遥感数据融合推理建筑物边界
建筑物屋顶点云边界指导生成影像缓冲区。点云处理之后得到了具体的单体建筑物屋顶边界点云，给定缓冲半径，将点云映射到影像上进行缓冲区生成。如图7所示。
利用影像缓冲区进行线段特征的筛选。对于原始影像中所有的线段进行判断，其位于哪个缓冲区内部即归属于哪个建筑物。(一个线段70％均在缓冲区内部时认为其属于该缓冲区)。如图8所示。
对初筛后的线特征进行分类。初筛后线段存在诸多噪声，同一建筑物边界可能有多个线段组成，因此需要判断哪些线段属于同一类(在组成建筑物时大概率划分在同一边界)。如图9所示。分类判断准则：
(1)线段向量间夹角(斜率)小于阈值。
(2)线段间距小于阈值(线段间距定义为两线段中点连线在斜率方向和垂直斜率方向投影的加权和)。
按照距离阈值将点云顺序边界投影到分类后线段，如图10所示。考虑上一步得到的所有影像轮廓线段，将顺序点云边界的每个点依次投影到各条线段上，保留投影距离小于阈值的投影垂足；完成后，查看每一条线段包含的垂足数量，保留大于阈值个数的线段，完成线段的二次筛选(保留时只截取第一个垂足和最末尾垂足之间的线段部分)。结果如图11所示。
二次筛选后的线段按照类别进行最优选择，如图12所示。在以上步骤中分别完成了对影像线段的分类和二次筛选，去除了部分细小噪声。考虑真边界的唯一性，需要按照一定规则对各类别线段进行最优挑选，结果如图13所示。规则：
(1)线段长度(推理越长越有理由作为建筑物真边界)。
(2)各垂足与对应原始点云点之间距离(投影距离)的分布(建筑物真边界理论上各点投影距离一致)。
(3)垂足离散度即单位长度的线段所包含的垂足点数(越多则该边作为
真边界的置信度越大)。
对最优选择的线段进行方向性检测。对于最优线段集需要进行方向性检测，筛选出错误线段，如图14所示。步骤：
(1)计算各线段斜率，按照一定差异值进行分类。
(2)对于无法归类的单一线段判断其长度。
(3)长度小于阈值则舍弃。
空缺边检测与补齐，如图15所示。得到的顺序最优边集合(完成方向检测排除了错误短边)进行闭合性与完成行检测，判断是否存在缺边，结果如图16所示。缺边判断规则：
(1)前后线段间相邻且平行。
(2)线段相邻且前一线段末尾点与后一线段起始点(对应顺序点云集的点)之间空缺点数超过阈值。
对于缺边位置，需要进行补全或推理，缺边补全规则：
(1)在原始影像初筛线段中查找通过空缺边的前一线段尾点与后一线段起点的线段。
(2)补全规则(1)查找到多条符合规则线段则按照如下规则筛选最优线段：
a.斜率属于最优线段集按照斜率分类后的某一类。
b.长度与空缺边的前一线段尾点与后一线段起点间连线长度接近。
(3)补全规则(1)找不到符合规则的线段。对于属于缺边判断规则(1)的缺边；作过前一线段结尾点与后一线段起点连线中点且垂直于前一线段和后一线段的垂线，两垂足连线作为推理缺失边。对于属于缺边规则(2)的缺边；参考顺序点云中的空缺点，对这些空缺点进行dp规则化得到推理空缺边。
线段集合排序与闭合。最优筛选后得到的线段不具有先后性，无法得知多边形拓扑关系，因此需要按照点云顺序点集作为指导对各线段进行排序。排序之后即可闭合线段生成建筑物多边形边界。结果如图17所示。
6.建筑物三维模型生成
通过oesm和dem反算得到建筑物多边形各拐点的建筑物高。根据建筑物高和平面坐标生成建筑物三维模型。结果如图18所示。由于已知建筑物屋顶角点的oesm坐标，该坐标包含dem上的x，y坐标以及与该dem单元保持共线关系的dsm单元的高程作为z值，可以根据如下所示公式计算出角点对应的真实三维坐标即dsm坐标。
其中，xs，ys，zs为摄影中心的坐标。在计算得到建筑物屋顶角点的dsm坐标后，结合建筑物高，可以得到屋顶角点对应的地面点的真实三维坐标，根据建筑物屋顶角点及其对应的地面点的真实三维坐标即可构建建筑物的线框三维模型。
应当理解的是，本说明书未详细阐述的部分均属于现有技术。
应当理解的是，上述针对较佳实施例的描述较为详细，并不能因此而认为是对本发明专利保护范围的限制，本领域的普通技术人员在本发明的启示下，在不脱离本发明权利要求所保护的范围情况下，还可以做出替换或变形，均落入本发明的保护范围之内，本发明的请求保护范围应以所附权利要求为准。
- 还没有人留言评论。精彩留言会获得点赞！