CN117237565B - 基于高分辨率卫星立体影像的建筑物白模制作方法 - Google Patents

Abstract

本发明公开了一种基于高分辨率卫星立体影像的建筑物白模制作方法，包括如下步骤，S1、DSM滤波获取DTM；S2、DSM形态边缘提取；S3、DSM形态边缘与建筑物轮廓匹配；S4、建筑物轮廓线两两判断包含关系；S5、建筑物轮廓线内部高程范围统计；S6、三角网构建与保存。优点是：通过变换重叠窗口大小和基于RANSAC拟合算法的有效点确定方法解决DSM滤波中对过滤建筑物大小和基地坡度的自适应的问题。根据建筑物偏移量的分布特点和获取建筑物轮廓的数据基础不同，把建筑物偏移量的计算分为两大类，根据类别分别采用相应的搜索方法获取轮廓线节点的偏移量，提升偏移量估算的准确度。基于各个轮廓线包含关系，采用合理方式获取最佳高程。

Description

基于高分辨率卫星立体影像的建筑物白模制作方法

技术领域

本发明涉及摄影测量与遥感技术领域，尤其涉及一种基于高分辨率卫星立体影像的建筑物白模制作方法。

背景技术

大区域地表三维模型一般可以通过无人机、飞艇等中、低空平台搭载的激光雷达或多镜头相机获取的多视角影像提取，多镜头相机方案一般称为倾斜摄影。地表三维模型包括几何模型和纹理贴图两部分。倾斜摄影或激光雷达获取的点云是构建几何模型的基础。倾斜摄影对构建一些细部、镂空结构效果不佳，但成本较低，在城市级三维建模领域应用更加广泛。航空倾斜摄影获取的影像分辨率一般在1～20cm之间，2～5cm是较为常用的分辨率范围。

随着卫星平台、相机等相关技术的发展，通过卫星相机获取多视角高分辨率影像进行三维模型构建逐渐成为三维建模新的技术手段之一。多视角卫星影像在时间获取一致性、分辨率、观测角度等方面受轨道周期、轨道高度的制约较大，获取成本较高。高分七号、高分十四号等国产高分辨率卫星是两线阵立体成像卫星，前、后视夹角为31度，无法调整。即使多次侧摆成像，也无法克服前后视遮挡这一固有缺陷。

发明内容

本发明的目的在于提供一种基于高分辨率卫星立体影像的建筑物白模制作方法，从而解决现有技术中存在的前述问题。

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

一种基于高分辨率卫星立体影像的建筑物白模制作方法，包括如下步骤，

S1、DSM滤波获取DTM：

将原始DSM按照平面网格坐标划分为一定像素尺寸的正方形区域，对每个区域统计最小和最大高程，并将最小到最大高程之间的高程划分为多等段，统计位于每个高程段内符合共面性且面积超过一定门限的高程点数，获取一幅与原始DSM坐标基准和分辨率一致的DTM；

S2、DSM形态边缘提取：

利用原始DSM网格点邻域内的最小、最大高程分别构建最小高程DSM和最大高程DSM，分别计算最小高程DSM、最大高程DSM与原始DSM之间对应网格点的差值绝对值，并基于差值绝对值的大小生成最小高程二值图和最大高程二值图；对比两幅二值图的像素值，根据对比结果改动最小高程二值图中像素的像素值，获取合成二值图；

S3、DSM形态边缘与建筑物轮廓匹配：

基于加密前后所有建筑物轮廓点坐标在合成二值图上的位置确定统计值的大小，并根据建筑物轮廓线来源是否与生成DSM的立体影像同源，分别采用一维搜索和二维搜索获取每个轮廓线节点的偏移量，通过将原始DSM和DTM的左下角经度和左下角纬度减去相应的偏移量，实现DSM形态边缘与建筑物轮廓的匹配；

S4、建筑物轮廓线两两判断包含关系：

判断每个建筑物轮廓节点的坐标矩形范围与其余所有矩形范围之间的包含与被包含关系；当该矩形范围包含其他矩形范围，利用射线法判断被包含矩形范围对应轮廓的所有节点是否都位于当前矩形对应轮廓范围内，并基于判断结果记录包含与被包含轮廓的序号和标记；

S5、建筑物轮廓线内部高程范围统计：

统计每个建筑物轮廓所有顶点的坐标在原始DSM上的行列号和矩形范围，计算矩形范围左下角对应原始DSM行列和矩形范围右上角对应原始DSM行列；通过逐网格点判断矩形范围内的网格点是否位于建筑物轮廓范围内；

S6、三角网构建与保存：

基于当前轮廓线是否包含其他轮廓线构建顶面平面三角网，计算轮廓线节点高程中值，当顶部为倾斜面时，利用其平面坐标重新计算轮廓节点的高程，并基于相应高程构建筑物立面三角网。

优选的，步骤S1中，统计位于每个高程段内符合共面性且面积超过一定门限的高程点数，具体包括如下内容，

S11、将高程段内全部的有效高程点表示为一个三维坐标序列(X_i，Y_i，H_i)，给定另一个三维坐标序列(X_i，Y_i，H₀)，通过RANSAC算法计算两个三维坐标序列在一定误差门限m内的有效点数，若有效点数小于一定门限M，则该高程段的点标记加1；

其中，X_i＝L_i×100000；Y_i＝B_i×100000；i＝0，1，…，N-1；L_i、B_i、H_i分别为第i个点的经度、纬度、高度；

S12、第二次执行上述操作时各窗口沿影像列方向偏移半个窗口，第三次执行上述操作时沿影像行偏移半个窗口，第四次执行上述操作时行列两个方向均偏移半个窗口；每个窗口都进行有效点统计和标记，标记为4的点确定为需要过滤的点；

S13、将需要过滤的点的高程设为无效值，统计无效区域与周围有效区域邻接的有效高程点，按距离倒数加权法拟合无效值的高程，执行完毕后得到一幅与原始DSM坐标基准、分辨率一致的DTM。

优选的，步骤S2具体包括如下内容，

S21、统计每个原始DSM网格点5×5邻域内的最大、最小高程，利用所有最小高程形成一幅新的DSM，即最小高程DSM；利用所有最大高程形成一幅新的DSM，即最大高程DSM；最小高程DSM以及最大高程DSM的每个网格点与原始DSM的网格点坐标一一对应；

S22、将最小高程DSM与原始DSM对应网格点相减并取绝对值，若该绝对值大于门限dh，则保存标签为255，否则保存标签为0，形成一幅与原始DSM长宽一致的最小高程二值图；将最大高程DSM与原始DSM对应网格点相减并取绝对值，若该绝对值大于门限dh，则保存标签为255，否则保存标签为0，形成一幅与原始DSM长宽一致的最大高程二值图；

S23、逐像素对比最小高程二值图和最大高程二值图，若最小高程二值图中某像素值为0，而最大高程二值图对应的像素值为255，则将最小高程二值图中的该像素值设为255，遍历所有的像素值，获取一幅新的二值图，即合成二值图；完成原始DSM形态边缘提取。

优选的，步骤S3具体包括如下内容，

S31、将建筑物轮廓线顶点进行加密，即每两个相邻顶点之间距离若超过门限，则在这两个顶点连线中间加入一个加密点；设置统计值T为0，计算所有建筑物轮廓点坐标在合成二值图上的位置，若该点值为255，则T加1；

S33、若建筑物轮廓线来源与生成DSM的立体影像非同源，则采用二维搜索获取每个轮廓线的偏移量；若建筑物轮廓线与生成DSM的立体影像同源，则采用一维搜索获取每个轮廓线的偏移量；

S34、将原始DSM和DTM的左下角经度和左下角纬度分别减去相应的偏移量，从而实现DSM形态边缘与建筑物轮廓匹配。

优选的，采用二维搜索获取偏移量具体为，在一个R×R范围内逐像素移动所有建筑物轮廓线顶点坐标，即将所有顶点坐标加上一个偏移量(ΔL，ΔB)＝((i′-R/2)×dB，(j′-R/2)×dL)，dL和dB分别为原始DSM在列和行方向的分辨率，将得到R×R个统计值N_i′_j′，N_i′_j′中的最大值对应的偏移量为(ΔL_max，ΔB_max)；将原始DSM和DTM的左下角经度减去ΔL_max，左下角纬度减去ΔB_max，从而实现DSM形态边缘与建筑物轮廓匹配；i′和j′为行、列号；

若已知工作区域所有建筑物的轮廓相对DSM轮廓的偏移量是个系统量，则取所有计算偏移量的中位数作为所有轮廓线节点的偏移量。

优选的，利用生成DSM的立体影像中某一视角的影像采集建筑物轮廓线，令建筑物轮廓线与生成DSM的立体影像同源；具体地，

选择观测角度最小的影像作为基础原始影像，存在两种方式获取与DSM同一坐标系DTM轮廓线；一种是在基础原始影像上采集建筑物轮廓线，此时，轮廓线的坐标为像平面坐标，该轮廓线为像面轮廓线，将像面轮廓线的像面坐标经基础原始影像RPC投射到DTM上，即可获取轮廓线每个顶点的物方坐标，即得到DTM轮廓线；另一种是通过DTM和基础原始影像及其RPC获取正射影像，在正射影像上采集DTM轮廓线；获取DIM轮廓线后通过一维搜索获取每个轮廓线节点的偏移量；

采用一维搜索时，不仅通过每个搜索位置所有建筑物轮廓点坐标在合成二值图位置值为255的点数T的大小作为匹配度的唯一度量值，同时考虑每个搜索位置的轮廓线内部DSM高程值与搜索高度H_s的匹配度；通过计算轮廓内所有DSM网格点高程与搜索高度H_s的差值，将差值取绝对值排序后，计算平均值H_a，将T_a＝T-H_a作为匹配度的度量值；T_a最大时对应的步数为dT_a，则最佳匹配位置的偏移量为(ΔL*dT_a，ΔB*dT_a)，将轮廓线所有坐标加上该偏移量即为准确的偏移量。

优选的，一维搜索中，搜索方向用经度方向和纬度方向的搜索步长(ΔL′，ΔB′)来表示，通过RPC计算(ΔL′，ΔB′)，具体过程为，

A1、由一个物方坐标(L，B，H)通过基础影像RPC计算像面坐标(r，c)；

其中，L、B、H分别为为基础原始影像RPC的LONG_OFF、LAT_OFF、HEIGHT_OFF参数；r、c分别为行列坐标；

A2、由另一个物方坐标(L，B，H+1)通过基础影像RPC计算像面坐标(r₁，c₁)；

其中，+1代表搜索方向是向上搜索，若为-1则为向下搜索；

A3、由像面坐标(r₁，c₁)和H通过基础影像RPC计算物方坐标(L₁，B₁，H)；

A4、计算(r，c)和(r₁，c₁)之间的像面距离Δd；

A5、计算搜索步长(ΔL′，ΔB′)：ΔL′＝(L₁-L)/Δd/2，ΔB′＝(B₁-B)/Δd/2；

A6、搜索步数为Hoff/Δd*2，每个搜索步长对应的高度变化为ΔH＝Δd/2；

其中，Hoff为搜索高程范围。

优选的，步骤S4具体为，统计每个建筑物轮廓节点的坐标矩形范围；对每个矩形范围，判断其与其余所有矩形范围的包含与被包含关系，若当前矩形范围包含其他矩形范围，则通过射线法进一步判断被包含矩形范围对应轮廓的所有节点是否全部位于当前矩形范围对应的轮廓范围内，若是，则为当前轮廓记录被包含轮廓的序号并标记0，同时为被包含轮廓记录当前轮廓的序号并标记1；反之，若当前轮廓被其他轮廓包含，则为当前轮廓记录包含轮廓的序号并标记1，同时为包含轮廓记录当前轮廓的序号并标记0。

优选的，步骤S6具体包括如下内容，

S61、如果当前轮廓不包含其它轮廓，利用高程均为H₀的建筑物轮廓节点三维坐标构建狄洛尼三角网TIN；若当前轮廓包含其它轮廓，则两个轮廓节点高程都设为H0，两组改变高程的轮廓节点都参与构建狄洛尼三角网TIN；

S62、利用冒泡法对轮廓线节点高程H_i排序，计算高程最大值H_max1和中值H_a1；

S63、利用仿射变换公式计算(L₀，B₀)、(L₀+0.00001，B₀)、(L₀，B₀+0.00001)三个平面坐标对应的高程H_0、H_1和H_2，计算H_1与H_0差值的绝对值H₀₁和H_2与H_0差值的绝对值H₀₂，计算H₀₁和H₀₂中的最大值H_max0，若H_max0大于0.5，则认为顶部为倾斜面，所有节点的高程根据其平面坐标通过仿射变换公式重新计算高程；否则，所有节点高程统一设为H_a；

S64、建筑物的立面三角网通过相邻两个节点对应的顶部高程和底部高程形成的四个共面点构建；若相邻两点的平面坐标为(L₀，B₀)和(L₁，B₁)，对应顶部高程分别为H₀和H₁，对应DTM高程分别为H₂和H₃，则四个顶点坐标为G₀(L₀，B₀，H₀)、G₁(L₁，B₁，H₁)、G₂(L₁，B₁，H₃)、G₃(L₀，B₀，H₂)，构成两个三角形分别为(G₀，G₁，G₃)和(G₀，G₂，G₃)；利用同样的方法生成其余邻接顶点对应的三角形；

S65、保存建筑物顶部三角网和立面三角网。

优选的，S61中，假设三角网包含N_t0个三角形，对第p个三角形三个顶点平面坐标取平均得到三角形重心平面坐标，通过射线法判断该重心平面坐标是否位于建筑物轮廓线内部，如果存在轮廓包含关系，则利用射线法判断该重心平面坐标是否位于被包含建筑物轮廓线内部；p＝0，1，…，N_t0-1；

在不存在轮廓包含关系的前提下，若重心点不在轮廓线内部，则该三角形不参与TIN构建；如果存在轮廓包含关系，若重心点不在外轮廓线内部，则该三角形不参与TIN构建，若重心点位于内轮廓线内部，则该三角形不参与TIN构建，但记录为内部三角形；

完成所有三角形的判断之后，假设剩余三角形个数为N_t，计算每个剩余三角形重心点平面坐标在原始DSM上的高程H_cg；通过冒泡法排序计算最小高程H_min，最大高程H_max和中间高程H_a；g＝0，1，…，N_t-1。

本发明的有益效果是：1、本发明通过变换重叠窗口大小和基于RANSAC拟合算法的有效点确定方法解决对过滤建筑物大小和基地坡度的自适应这一DSM滤波难点问题。2、本发明根据建筑物偏移量的分布特点和获取建筑物轮廓的数据基础不同，把建筑物偏移量的计算分为两大类，根据类别分别采用二维搜索和一维搜索的方式确定偏移量，一维搜索中除了建筑物轮廓和DSM形态边缘的匹配度之外，还有另外一个匹配度可以利用，那就是偏移量对应的高度变化和建筑物内部实际统计的高度的一致性，这两个额外的约束条件使得偏移量估计的准确度显著提升。3、本发明在分析各个轮廓线包含关系的基础上，通过合理的统计方法获取最佳高程。

附图说明

图1是本发明实施例中方法的流程图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施方式仅仅用以解释本发明，并不用于限定本发明。

如图1所示，本实施例中，提供了一种基于高分辨率卫星立体影像的建筑物白模制作方法，包括如下步骤，

一、DSM滤波获取DTM：

将原始DSM按照平面网格坐标划分为一定像素尺寸的正方形区域，假设该区域大小为m×m矩形窗口，对每个区域统计最小高程h_min和最大高程h_max，并将最小到最大高程之间的高程划分为n等段，统计位于每个高程段内符合共面性且面积超过一定门限的高程点数，获取一幅与原始DSM坐标基准和分辨率一致的DTM。

统计位于每个高程段内符合共面性且面积超过一定门限的高程点数，具体包括如下内容，

1、将高程段内全部的有效高程点表示为一个三维坐标序列(X_i，Yi，H_i)，给定另一个三维坐标序列(X_i，Y_i，H₀)，可以看出这是一个平面坐标与有效高程点平面坐标一一对应但高程均为H0的坐标序列。通过RANSAC算法计算两个三维坐标序列在一定误差门限m内的有效点数，若有效点数小于一定门限M，则该高程段的点标记加1；

其中，X_i＝L_i×100000；Y_i＝B_i×100000；i＝0，1，…，N-1；L_i、B_i、H_i分别为第i个点的经度、纬度、高度；

2、第二次执行上述操作时各窗口沿影像列方向偏移半个窗口，第三次执行上述操作时沿影像行偏移半个窗口，第四次执行上述操作时行列两个方向均偏移半个窗口；每个窗口都进行有效点统计和标记，标记为4的点确定为需要过滤的点；

3、将需要过滤的点的高程设为无效值，统计无效区域与周围有效区域邻接的有效高程点，按距离倒数加权法拟合无效值的高程，执行完毕后得到一幅与原始DSM坐标基准、分辨率一致的DTM。

本实施例中，采用多次迭代方法执行该过滤步骤，一共迭代执行三次，每次采用不同的参数，三次采用的m参数为100，70，30，M参数为1000，700，100。

二、DSM形态边缘提取：

利用原始DSM网格点邻域内的最小、最大高程分别构建最小高程DSM和最大高程DSM，分别计算最小高程DSM、最大高程DSM与原始DSM之间对应网格点的差值绝对值，并基于差值绝对值的大小生成最小高程二值图和最大高程二值图；对比两幅二值图的像素值，根据对比结果改动最小高程二值图中像素的像素值，获取合成二值图。

具体包括如下内容，

1、统计每个原始DSM网格点5×5邻域内的最大、最小高程，利用所有最小高程形成一幅新的DSM，即最小高程DSM；利用所有最大高程形成一幅新的DSM，即最大高程DSM；最小高程DSM以及最大高程DSM的每个网格点与原始DSM的网格点坐标一一对应；

2、将最小高程DSM与原始DSM对应网格点相减并取绝对值，若该绝对值大于门限dh，则保存标签为255，否则保存标签为0，形成一幅与原始DSM长宽一致的最小高程二值图；将最大高程DSM与原始DSM对应网格点相减并取绝对值，若该绝对值大于门限dh，则保存标签为255，否则保存标签为0，形成一幅与原始DSM长宽一致的最大高程二值图；

3、逐像素对比最小高程二值图和最大高程二值图，若最小高程二值图中某像素值为0，而最大高程二值图对应的像素值为255，则将最小高程二值图中的该像素值设为255，遍历所有的像素值，获取一幅新的二值图，即合成二值图；完成原始DSM形态边缘提取。

本实施例中，合成二值图体现了DSM的形态边缘信息，255值主要集中在建筑物边缘等高度变化较大的区域。

三、DSM形态边缘与建筑物轮廓匹配：

基于加密前后所有建筑物轮廓点坐标在合成二值图上的位置确定统计值的大小，并根据建筑物轮廓线来源是否与生成DSM的立体影像同源，分别采用一维搜索和二维搜索获取每个轮廓线节点的偏移量，通过将原始DSM和DTM的左下角经度和左下角纬度减去相应的偏移量，实现DSM形态边缘与建筑物轮廓的匹配。

具体包括如下内容，

1、将建筑物轮廓线顶点进行加密，即每两个相邻顶点之间距离若超过门限，则在这两个顶点连线中间加入一个加密点；设置统计值T为0，计算所有建筑物轮廓点(包括加密前后的)坐标在合成二值图上的位置(加0.5取整数)，若该点值为255，则T加1；

3、若建筑物轮廓线来源与生成DSM的立体影像非同源，则采用二维搜索获取每个轮廓线的偏移量；若建筑物轮廓线与生成DSM的立体影像同源，则采用一维搜索获取每个轮廓线的偏移量；

(1)、采用二维搜索获取偏移量具体为，在一个R×R范围内逐像素移动所有建筑物轮廓线顶点坐标，即将所有顶点坐标加上一个偏移量(ΔL，ΔB)＝((i′-R/2)×dB，(j′-R/2)×dL)，dL和dB分别为原始DSM在列和行方向的分辨率，将得到R×R个统计值N_i′_j′，N_i′_j′中的最大值对应的偏移量为(ΔL_max，ΔB_max)；将原始DSM和DTM的左下角经度减去ΔL_max，左下角纬度减去ΔB_max，从而实现DSM形态边缘与建筑物轮廓匹配，即消除了建筑物轮廓线与原始DSM中间的坐标偏移；i′和j′为行、列号，i′≥0，i′＜R，j′≥0，j′＜R。

上述做法的前提是每个轮廓线的偏移量都是不同的，如果已知工作区域所有建筑物的轮廓相对DSM轮廓的偏移是个系统量，则取所有计算偏移量的中位数作为所有轮廓线节点的偏移量。

(2)、利用生成DSM的立体影像中某一视角的影像采集建筑物轮廓线，令建筑物轮廓线与生成DSM的立体影像同源；具体地，选择观测角度最小的影像作为基础原始影像(基础原始影像可以是全色影像，也可以是与同视角多光谱影像融合后的RGB彩色影像)，存在两种方式获取与DSM同一坐标系DTM轮廓线；一种是在基础原始影像上采集建筑物轮廓线，此时，轮廓线的坐标为像平面坐标，该轮廓线为像面轮廓线，将像面轮廓线的像面坐标经基础原始影像RPC投射到DTM上，即可获取轮廓线每个顶点的物方坐标，即得到DTM轮廓线；另一种是通过DTM和基础原始影像及其RPC获取正射影像，在正射影像上采集DTM轮廓线；获取DTM轮廓线后通过一维搜索获取每个轮廓线节点的偏移量。

一维搜索中，搜索方向用经度方向和纬度方向的搜索步长(ΔL′，ΔB′)来表示，通过RPC计算(ΔL′，ΔB′)，具体过程为，

A1、由一个物方坐标(L，B，H)通过基础影像RPC计算像面坐标(r，c)；

其中，L、B、H分别为为基础原始影像RPC的LONG_OFF、LAT_OFF、HEIGHT_OFF参数；r、c分别为行列坐标；

A2、由另一个物方坐标(L，B，H+1)通过基础影像RPC计算像面坐标(r₁，c₁)；

其中，+1代表搜索方向是向上搜索，若为-1则为向下搜索；

A3、由像面坐标(r₁，c₁)和H通过基础影像RPC计算物方坐标(L₁，B₁，H)；

A4、计算(r，c)和(r₁，c₁)之间的像面距离Δd；

A5、计算搜索步长(ΔL′，ΔB′)：ΔL′＝(L₁-L)/Δd/2，ΔB′＝(B₁-B)/Δd/2；

A6、搜索步数为Hoff/Δd*2，每个搜索步长对应的高度变化为ΔH＝Δd/2；

其中，Hoff为搜索高程范围。

本实施例中，采用一维搜索时，不仅通过每个搜索位置所有建筑物轮廓点(包括加密前后的)坐标在合成二值图位置(加0.5取整数)值为255的点数T的大小作为匹配度的唯一度量值，同时考虑每个搜索位置的轮廓线内部DSM高程值与搜索高度H_s的匹配度；通过计算轮廓内所有DSM网格点高程与搜索高度H_s的差值，将差值取绝对值排序后，计算平均值H_a，将T_a＝T-H_a作为匹配度的度量值；T_a最大时对应的步数为dT_a，则最佳匹配位置的偏移量为(ΔL*dT_a，ΔB*dT_a)，将轮廓线所有坐标加上该偏移量即为准确的偏移量。

4、将原始DSM和DTM的左下角经度和左下角纬度分别减去相应的偏移量，从而实现DSM形态边缘与建筑物轮廓匹配。

四、建筑物轮廓线两两判断包含关系：

判断每个建筑物轮廓节点的坐标矩形范围与其余所有矩形范围之间的包含与被包含关系；当该矩形范围包含其他矩形范围，利用射线法判断被包含矩形范围对应轮廓的所有节点是否都位于当前矩形对应轮廓范围内，并基于判断结果记录包含与被包含轮廓的序号和标记。

具体为，统计每个建筑物轮廓节点的坐标矩形范围；对每个矩形范围，判断其与其余所有矩形范围的包含与被包含关系，若当前矩形范围包含其他矩形范围，则通过射线法进一步判断被包含矩形范围对应轮廓的所有节点是否全部位于当前矩形范围对应的轮廓范围内，若是，则为当前轮廓记录被包含轮廓的序号并标记0，同时为被包含轮廓记录当前轮廓的序号并标记1；反之，若当前轮廓被其他轮廓包含，则为当前轮廓记录包含轮廓的序号并标记1，同时为包含轮廓记录当前轮廓的序号并标记0。

五、建筑物轮廓线内部高程范围统计：

统计每个建筑物轮廓所有顶点的坐标在原始DSM上的行列号和矩形范围，计算矩形范围左下角对应原始DSM行列(r_min，c_min)和矩形范围右上角对应原始DSM行列(r_max，c_max)；通过逐网格点判断矩形范围内的网格点是否位于建筑物轮廓范围内。

具体地：假设第i个轮廓有M2个顶点，第j(j＝0，1，...，M2-1)个顶点的坐标为(L_i，B_i)，则该顶点在DSM上的行列号为(r_j，c_j)。用射线法判断矩形范围内的格网点是否位于(r_j，c_j)点序列形成的多边形内。假设一共有Q个点位于该多边形内，称这Q个点为内点。通过RANSAC算法判断这些内点的共面性。首先将Q个内点标表示为三维坐标序列的形式，假设第i个内点的坐标为(L_i，B_i，H_i)。将L_i和B_i都乘以100000转换为近似为以米为单位的数字，转换后的坐标表示为(L_1i，B_1i，H_i)。RANSAC算法用于计算两组坐标点间的仿射变换参数及一定误差范围内符合该仿射变换的点位。(L_1i，B_1i，H_i)是真实点位，与之对应的仿射变换点位通过将Hi统一设为H₀获得，平面坐标不变。也就是通过RANSAC算法获取(L_1i，B_1i，H_i)三维点序列和(L_1i，B_1i，H₀)三维点序列间的仿射变换参数和符合仿射变换的点。仿射变换形式为：H＝a×L+b×B+c，L和B表示经度、纬度坐标，H表示高程，a、b、c为仿射变换参数。

六、三角网构建与保存：

基于当前轮廓线是否包含其他轮廓线构建顶面平面三角网，计算轮廓线节点高程中值，当顶部为倾斜面时，利用其平面坐标重新计算轮廓节点的高程，并基于相应高程构建筑物立面三角网。

具体包括如下内容，

1、构建顶部平面三角网：如果当前轮廓不包含其它轮廓，利用高程均为H₀的建筑物轮廓节点三维坐标构建狄洛尼三角网TIN；若当前轮廓包含其它轮廓，则两个轮廓节点高程都设为H0，两组改变高程的轮廓节点都参与构建狄洛尼三角网TIN；

假设三角网包含N_t0个三角形，对第p个三角形三个顶点平面坐标取平均得到三角形重心平面坐标，通过射线法判断该重心平面坐标是否位于建筑物轮廓线内部，如果存在轮廓包含关系，则利用射线法判断该重心平面坐标是否位于被包含建筑物轮廓线内部；p＝0，1，…，N_t0-1；

在不存在轮廓包含关系的前提下，若重心点不在轮廓线内部，则该三角形不参与TIN构建；如果存在轮廓包含关系，若重心点不在外轮廓线内部，则该三角形不参与TIN构建，若重心点位于内轮廓线内部，则该三角形不参与TIN构建，但记录为内部三角形；

完成所有三角形的判断之后，假设剩余三角形个数为N_t，计算每个剩余三角形重心点平面坐标在原始DSM上的高程H_cg；通过冒泡法排序计算最小高程H_min，最大高程H_max和中间高程H_a；g＝0，1，…，N_t-1。

2、计算轮廓线节点高程中值：利用冒泡法对轮廓线节点高程H_i排序，计算高程最大值H_max1和中值H_a1；

3、倾斜度计算：利用仿射变换公式计算(L₀，B₀)、(L₀+0.00001，B₀)、(L₀，B₀+0.00001)三个平面坐标对应的高程H_0、H_1和H_2，计算H_1与H_0差值的绝对值H₀₁和H_2与H_0差值的绝对值H₀₂，计算H₀₁和H₀₂中的最大值H_max0，若H_max0大于0.5，则认为顶部为倾斜面，所有节点的高程根据其平面坐标通过仿射变换公式重新计算高程；否则，所有节点高程统一设为H_a；

4、构建建筑物立面三角网：建筑物的立面三角网通过相邻两个节点对应的顶部高程和底部高程形成的四个共面点构建；若相邻两点的平面坐标为(L₀，B₀)和(L₁，B₁)，对应顶部高程分别为H₀和H₁，对应DTM高程分别为H₂和H₃，则四个顶点坐标为G₀(L₀，B₀，H₀)、G₁(L₁，B₁，H₁)、G₂(L₁，B₁，H₃)、G₃(L₀，B₀，H₂)，构成两个三角形分别为(G₀，G₁，G₃)和(G₀，G₂，G₃)；利用同样的方法生成其余邻接顶点对应的三角形；

5、白模数据保存：保存建筑物顶部三角网和立面三角网。

通过采用本发明公开的上述技术方案，得到了如下有益的效果：

本发明提供了一种基于高分辨率卫星立体影像的建筑物白模制作方法，本发明通过变换重叠窗口大小和基于RANSAC拟合算法的有效点确定方法解决对过滤建筑物大小和基地坡度的自适应这一DSM滤波难点问题。本发明根据建筑物偏移量的分布特点和获取建筑物轮廓的数据基础不同，把建筑物偏移量的计算分为两大类，根据类别分别采用二维搜索和一维搜索的方式确定偏移量，一维搜索中除了建筑物轮廓和DSM形态边缘的匹配度之外，还有另外一个匹配度可以利用，那就是偏移量对应的高度变化和建筑物内部实际统计的高度的一致性，这两个额外的约束条件使得偏移量估计的准确度显著提升。本发明在分析各个轮廓线包含关系的基础上，通过合理的统计方法获取最佳高程。

以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视本发明的保护范围。

Claims (9)

1.一种基于高分辨率卫星立体影像的建筑物白模制作方法，其特征在于：包括如下步骤，

S1、DSM滤波获取DTM：

将原始DSM按照平面网格坐标划分为一定像素尺寸的正方形区域，对每个区域统计最小和最大高程，并将最小到最大高程之间的高程划分为多等段，统计位于每个高程段内符合共面性且面积超过一定门限的高程点数，获取一幅与原始DSM坐标基准和分辨率一致的DTM；步骤S1中，统计位于每个高程段内符合共面性且面积超过一定门限的高程点数，具体包括如下内容，

S11、将高程段内全部的有效高程点表示为一个三维坐标序列(X_i,Y_i,H_i)，给定另一个三维坐标序列(X_i,Y_i,H₀)，通过RANSAC算法计算两个三维坐标序列在一定误差门限m内的有效点数，若有效点数小于一定门限M，则该高程段的点标记加1；

其中，X_i＝L_i×100000；Y_i＝B_i×100000；i＝0,1,…,N-1；L_i、B_i、H_i分别为第i个点的经度、纬度、高度；

S12、第二次执行S11操作时各窗口沿影像列方向偏移半个窗口，第三次执行S11操作时沿影像行偏移半个窗口，第四次执行S11操作时行列两个方向均偏移半个窗口；每个窗口都进行有效点统计和标记，标记为4的点确定为需要过滤的点；

S13、将需要过滤的点的高程设为无效值，统计无效区域与周围有效区域邻接的有效高程点，按距离倒数加权法拟合无效值的高程，执行完毕后得到一幅与原始DSM坐标基准、分辨率一致的DTM；

S2、DSM形态边缘提取：

利用原始DSM网格点邻域内的最小、最大高程分别构建最小高程DSM和最大高程DSM，分别计算最小高程DSM、最大高程DSM与原始DSM之间对应网格点的差值绝对值，并基于差值绝对值的大小生成最小高程二值图和最大高程二值图；对比两幅二值图的像素值，根据对比结果改动最小高程二值图中像素的像素值，获取合成二值图；

S3、DSM形态边缘与建筑物轮廓匹配：

基于加密前后所有建筑物轮廓点坐标在合成二值图上的位置确定统计值的大小，并根据建筑物轮廓线来源是否与生成DSM的立体影像同源，分别采用一维搜索和二维搜索获取每个轮廓线节点的偏移量，通过将原始DSM和DTM的左下角经度和左下角纬度减去相应的偏移量，实现DSM形态边缘与建筑物轮廓的匹配；

S4、建筑物轮廓线两两判断包含关系：

判断每个建筑物轮廓节点的坐标矩形范围与其余所有矩形范围之间的包含与被包含关系；当该矩形范围包含其他矩形范围，利用射线法判断被包含矩形范围对应轮廓的所有节点是否都位于当前矩形对应轮廓范围内，并基于判断结果记录包含与被包含轮廓的序号和标记；

S5、建筑物轮廓线内部高程范围统计：

统计每个建筑物轮廓所有顶点的坐标在原始DSM上的行列号和矩形范围，计算矩形范围左下角对应原始DSM行列和矩形范围右上角对应原始DSM行列；通过逐网格点判断矩形范围内的网格点是否位于建筑物轮廓范围内；

S6、三角网构建与保存：

基于当前轮廓线是否包含其他轮廓线构建顶面平面三角网，计算轮廓线节点高程中值，当顶部为倾斜面时，利用其平面坐标重新计算轮廓节点的高程，并基于相应高程构建建筑物立面三角网；保存顶面平面三角网和立面三角网信息，获取建筑物白模数据。

2.根据权利要求1所述的基于高分辨率卫星立体影像的建筑物白模制作方法，其特征在于：步骤S2具体包括如下内容，

S21、统计每个原始DSM网格点5×5邻域内的最大、最小高程，利用所有最小高程形成一幅新的DSM，即最小高程DSM；利用所有最大高程形成一幅新的DSM，即最大高程DSM；最小高程DSM以及最大高程DSM的每个网格点与原始DSM的网格点坐标一一对应；

S22、将最小高程DSM与原始DSM对应网格点相减并取绝对值，若该绝对值大于门限dh，则保存标签为255，否则保存标签为0，形成一幅与原始DSM长宽一致的最小高程二值图；将最大高程DSM与原始DSM对应网格点相减并取绝对值，若该绝对值大于门限dh，则保存标签为255，否则保存标签为0，形成一幅与原始DSM长宽一致的最大高程二值图；

S23、逐像素对比最小高程二值图和最大高程二值图，若最小高程二值图中某像素值为0，而最大高程二值图对应的像素值为255，则将最小高程二值图中的该像素值设为255，遍历所有的像素值，获取一幅新的二值图，即合成二值图；完成原始DSM形态边缘提取。

3.根据权利要求1所述的基于高分辨率卫星立体影像的建筑物白模制作方法，其特征在于：步骤S3具体包括如下内容，

S31、将建筑物轮廓线顶点进行加密，即每两个相邻顶点之间距离若超过门限，则在这两个顶点连线中间加入一个加密点；设置统计值T为0，计算所有建筑物轮廓点坐标在合成二值图上的位置，若该点值为255，则T加1；

S33、若建筑物轮廓线来源与生成DSM的立体影像非同源，则采用二维搜索获取每个轮廓线的偏移量；若建筑物轮廓线与生成DSM的立体影像同源，则采用一维搜索获取每个轮廓线的偏移量；

S34、将原始DSM和DTM的左下角经度和左下角纬度分别减去相应的偏移量，从而实现DSM形态边缘与建筑物轮廓匹配。

4.根据权利要求3所述的基于高分辨率卫星立体影像的建筑物白模制作方法，其特征在于：采用二维搜索获取偏移量具体为，在一个R×R范围内逐像素移动所有建筑物轮廓线顶点坐标，即将所有顶点坐标加上一个偏移量(ΔL,ΔB)＝((i′-R/2)×dB,(j′-R/2)×dL)，dL和dB分别为原始DSM在列和行方向的分辨率，将得到R×R个统计值N_i′_j′，N_i′_j′中的最大值对应的偏移量为(ΔL_max,ΔB_max)；将原始DSM和DTM的左下角经度减去ΔL_max，左下角纬度减去ΔB_max，从而实现DSM形态边缘与建筑物轮廓匹配；i′和j′为行、列号；

若已知工作区域所有建筑物的轮廓相对DSM轮廓的偏移量是个系统量，则取所有计算偏移量的中位数作为所有轮廓线节点的偏移量。

5.根据权利要求3所述的基于高分辨率卫星立体影像的建筑物白模制作方法，其特征在于：利用生成DSM的立体影像中某一视角的影像采集建筑物轮廓线，令建筑物轮廓线与生成DSM的立体影像同源；具体地，

选择观测角度最小的影像作为基础原始影像，存在两种方式获取与DSM同一坐标系DTM轮廓线；一种是在基础原始影像上采集建筑物轮廓线，此时，轮廓线的坐标为像平面坐标，该轮廓线为像面轮廓线，将像面轮廓线的像面坐标经基础原始影像RPC投射到DTM上，即可获取轮廓线每个顶点的物方坐标，即得到DTM轮廓线；另一种是通过DTM和基础原始影像及其RPC获取正射影像，在正射影像上采集DTM轮廓线；获取DTM轮廓线后通过一维搜索获取每个轮廓线节点的偏移量；

采用一维搜索时，不仅通过每个搜索位置所有建筑物轮廓点坐标在合成二值图位置值为255的点数T的大小作为匹配度的唯一度量值，同时考虑每个搜索位置的轮廓线内部DSM高程值与搜索高度H_s的匹配度；通过计算轮廓内所有DSM网格点高程与搜索高度H_s的差值，将差值取绝对值排序后，计算平均值H_a，将T_a＝T-H_a作为匹配度的度量值；T_a最大时对应的步数为dT_a，则最佳匹配位置的偏移量为(ΔL*dT_a,ΔB*dT_a)，将轮廓线所有坐标加上该偏移量即为准确的偏移量。

6.根据权利要求5所述的基于高分辨率卫星立体影像的建筑物白模制作方法，其特征在于：一维搜索中，搜索方向用经度方向和纬度方向的搜索步长(ΔL′,ΔB′)来表示，通过RPC计算(ΔL′,ΔB′)，具体过程为，

A1、由一个物方坐标(L,B,H)通过基础影像RPC计算像面坐标(r,c)；

其中，L、B、H分别为为基础原始影像RPC的LONG_OFF、LAT_OFF、HEIGHT_OFF参数；r、c分别为行列坐标；

A2、由另一个物方坐标(L,B,H+1)通过基础影像RPC计算像面坐标(r₁,c₁)；

其中，+1代表搜索方向是向上搜索，若为-1则为向下搜索；

A3、由像面坐标(r₁,c₁)和H通过基础影像RPC计算物方坐标(L₁,B₁,H)；

A4、计算(r,c)和(r₁,c₁)之间的像面距离Δd；

A5、计算搜索步长(ΔL′,ΔB′)：ΔL^′＝(L₁-L)/Δd/2，ΔB^′＝(B₁-B)/Δd/2；

A6、搜索步数为Hoff/Δd*2,每个搜索步长对应的高度变化为ΔH＝Δd/2；

其中，Hoff为搜索高程范围。

7.根据权利要求1所述的基于高分辨率卫星立体影像的建筑物白模制作方法，其特征在于：步骤S4具体为，统计每个建筑物轮廓节点的坐标矩形范围；对每个矩形范围，判断其与其余所有矩形范围的包含与被包含关系，若当前矩形范围包含其他矩形范围，则通过射线法进一步判断被包含矩形范围对应轮廓的所有节点是否全部位于当前矩形范围对应的轮廓范围内，若是，则为当前轮廓记录被包含轮廓的序号并标记0，同时为被包含轮廓记录当前轮廓的序号并标记1；反之，若当前轮廓被其他轮廓包含，则为当前轮廓记录包含轮廓的序号并标记1，同时为包含轮廓记录当前轮廓的序号并标记0。

8.根据权利要求5所述的基于高分辨率卫星立体影像的建筑物白模制作方法，其特征在于：步骤S6具体包括如下内容，

S61、如果当前轮廓不包含其它轮廓，利用高程均为H₀的建筑物轮廓节点三维坐标构建狄洛尼三角网TIN；若当前轮廓包含其它轮廓，则两个轮廓节点高程都设为H0，两组改变高程的轮廓节点都参与构建狄洛尼三角网TIN；

S62、利用冒泡法对轮廓线节点高程H_i排序，计算高程最大值H_max1和中值H_a1；

S63、利用仿射变换公式计算(L₀,B₀)、(L₀+0.00001,B₀)、(L₀,B₀+0.00001)三个平面坐标对应的高程H_0、H_1和H_2，计算H_1与H_0差值的绝对值H₀₁和H_2与H_0差值的绝对值H₀₂，计算H₀₁和H₀₂中的最大值H_max0，若H_max0大于0.5，则认为顶部为倾斜面，所有节点的高程根据其平面坐标通过仿射变换公式重新计算高程；否则，所有节点高程统一设为H_a；

S64、建筑物的立面三角网通过相邻两个节点对应的顶部高程和底部高程形成的四个共面点构建；若相邻两点的平面坐标为(L₀,B₀)和(L₁,B₁)，对应顶部高程分别为H₀和H₁，对应DTM高程分别为H₂和H₃，则四个顶点坐标为G₀(L₀,B₀,H₀)、G₁(L₁,B₁,H₁)、G₂(L₁,B₁,H₃)、G₃(L₀,B₀,H₂)，构成两个三角形分别为(G₀,G₁,G₃)和(G₀,G₂,G₃)；利用同样的方法生成其余邻接顶点对应的三角形；

S65、保存建筑物顶面平面三角网和立面三角网。

9.根据权利要求8所述的基于高分辨率卫星立体影像的建筑物白模制作方法，其特征在于：S61中，假设三角网包含N_t0个三角形，对第p个三角形三个顶点平面坐标取平均得到三角形重心平面坐标，通过射线法判断该重心平面坐标是否位于建筑物轮廓线内部，如果存在轮廓包含关系，则利用射线法判断该重心平面坐标是否位于被包含建筑物轮廓线内部；p＝0,1,…,N_t0-1；

在不存在轮廓包含关系的前提下，若重心点不在轮廓线内部，则该三角形不参与TIN构建；如果存在轮廓包含关系，若重心点不在外轮廓线内部，则该三角形不参与TIN构建，若重心点位于内轮廓线内部，则该三角形不参与TIN构建，但记录为内部三角形；

完成所有三角形的判断之后，假设剩余三角形个数为N_t，计算每个剩余三角形重心点平面坐标在原始DSM上的高程H_cg；通过冒泡法排序计算最小高程H_min，最大高程H_max和中间高程H_a；g＝0,1,…,N_t-1。