一种针对多元遥感数据的房屋矢量化方法
【专利摘要】本发明公开了一种针对多元遥感数据的房屋矢量化方法，包括以下步骤：获取需要处理的可以代表房屋位置、形状的信息；采用八邻域边缘追踪得到有序的边缘点，同时生成房屋边缘二值图像，边缘点为“1”，非边缘点为“0”；利用Radon变换结合主轴分析得到房屋的两个相互正交的主要方向；将房屋的边缘线分割成三类边缘线段：第一主方向类、第二主方向类、不定向类；对第一主方向类和第二主方向类的边缘线段的精确定位；计算相邻边缘线段的拐角点；生成矢量化的房屋边界线。本发明可避免房屋规则化过程中选择初始点和处理顺序的麻烦和不利影响，还可减弱房屋提取结果误差的影响，可获取与房屋外形较为吻合、有效、规则化的房屋轮廓信息。
【专利说明】
一种针对多元遥感数据的房屋矢量化方法
技术领域
[0001] 本发明涉及一种针对多元遥感数据的房屋矢量化方法，属于图像处理技术领域。
【背景技术】
[0002] 作为一种主要地理空间研究信息，房屋（也称建筑物)对二维地图的更新、智慧城 市的构建、违章建筑的监管等具有重要意义。遥感技术的不断发展使得自动矢量化提取遥 感数据中房屋信息成为可能。
[0003]多元遥感数据包括:影像数据、LiDAR数据、DEM、DSM等数据。由于多种因素的干扰， 目前从多元遥感数据中提取的房屋斑块(边缘线)往往与房屋实际形状不符，需要对其进行 后处理一一去噪、规则化等。
[0004] 中国专利文献CN 104156988 A公开的建筑物轮廓规则化方法，主要利用现代房屋 相邻轮廓边界与正交的特性，基于迭代最小外包矩形进行房屋轮廓信息规则化，主要针对 LiDAR数据中获得的房屋点集。
[0005] 中国专利文献CN 104200212 A公开的建筑物外边界线提取方法，主要利用道格拉 斯算法和最小二乘直线拟合法，结合房屋的两个主要方向，以拟合直线方向的长度、拟合直 线与两个主要方向之间的差值等阈值参数来完成房屋边界线的规则化，对初始房屋提取结 果的精度依赖性较高。
[0006] 中国专利文献CN 102938066 B公开的重建建筑物外轮廓多边形的方法，针对DSM 数据和影像数据，利用两种数据的互补特性对房屋轮廓进行重建。
[0007] 无论是LiDAR数据，还是DSM数据，均具有准确的高程信息。高程是房屋区别与其周 围地物(道路、车辆、植被)的重要特征，保障了房屋提取结果的精度不至于太低。然而，准确 的高程信息并不易获取，这也导致利用多元遥感数据(尤其是只采用影像数据的情况)提取 的房屋锯齿状变形较为严重，现有的房屋轮廓信息后处理方法适用性不强。因此，如何对不 同种类的多元遥感数据中提取的房屋(斑块、边缘线、边界线等)进行矢量化，减弱房屋提取 结果的影响显得尤为重要。
[0008] 现有技术尤其是规则化方法需要面临最佳起始点、最佳处理顺序的选择问题，不 同的起始点和处理顺序可能带来完全不同的规则化结果。现有技术需要设置太多的阈值， 自动化程度不高。此外，现有技术主要针对房屋提取结果精度较好的情况，直接对提取结果 进行后处理，而实际应用中房屋提取结果往往并不理想，目前仍然缺乏有效的针对不同种 类影像数据中房屋的后处理方法。
【发明内容】
[0009] 本发明正是针对现有技术存在的不足，提供一种针对多元遥感数据的房屋矢量化 方法，可避免房屋规则化过程中选择初始点和处理顺序的麻烦和不利影响，还可减弱房屋 提取结果误差的影响，可获取与房屋外形较为吻合、有效、规则化的房屋轮廓信息，可对多 元遥感数据(LiDAR数据、DEM、DSM、航空影像数据、卫星影像数据、无人机影像数据等）的提 取结果进行后处理。
[0010] 为解决上述问题，本发明所采取的技术方案如下： 一种针对多元遥感数据的房屋矢量化方法，包括以下步骤： 步骤1，获取需要处理的可以代表房屋位置、形状的信息； 步骤2,获取房屋边缘:采用八邻域边缘追踪得到有序的边缘点，同时生成房屋边缘二 值图像，边缘点为"Γ，非边缘点为"0" ； 步骤3,利用Radon变换结合主轴分析得到房屋的两个相互正交的主要方向； 步骤4,将房屋的边缘线分割成三类边缘线段:第一主方向类、第二主方向类、不定向 类； 步骤5,对第一主方向类和第二主方向类上的所有边缘点在其法方向上进行模板匹配， 寻找最佳匹配点来完成对第一主方向类和第二主方向类的边缘线段的精确定位； 步骤6,在准确获取第一主方向类和第二主方向类的边缘线段上边缘点的基础上，计算 相邻边缘线段的拐角点； 步骤7，生成矢量化的房屋边界线。
[0011] 本发明与现有技术相比较，本发明的实施效果如下： 本发明所述的一种针对多元遥感数据的房屋矢量化方法，以检测或提取的房屋斑块 (或房屋点集、边缘线)作为已知数据，在基于扩展理论对房屋边缘线段进行多类分割的 基础上，针对不同遥感数据的房屋提取结果分别进行相应的处理:对房屋提取精度较高的 情况，利用房屋形状这一先验知识，结合直线拟合求拐角点；对房屋提取精度不高的情况， 利用影像数据在边缘处有较高的采样率这一特征，构建自适应的边缘模板对每一条边缘线 段进行精确定位，在获取准确边缘线段的基础上，结合直线拟合求拐角点；可获取与房屋外 形较为吻合、有效、规则化的房屋轮廓信息。
[0012] 本发明所述的一种针对多元遥感数据的房屋矢量化方法主要具有以下特点和进 步之处： 1、利用房屋大多具有主要方向的特征，以Radon变换结合主轴分析取代传统的Hough变 换来检查房屋的两个垂直的主要方向，计算量少，检测结果快捷，有效。
[0013] 2、针对边缘线分段问题，引入多类分割思想将其转化为能量函数的最小化问题， 这样既考虑了每一个边缘点的方向信息，也利用了相邻边缘点趋于同一类的这一先验知 识，可以实现近似全局最优的分类结果，通过将原始拐角点及锯齿状变形区域分为不定方 向类，可避免了选择初始点和处理顺序的麻烦和不利影响。此外，整个边缘线段的多类分割 过程是全自动的，不需要对参数进行太多的调整。
[0014] 3、针对房屋边缘线的锯齿状变形导致房屋提取精度不高的情况，通过与影像数据 相结合，采用构建的自适应边缘模板对分类后的边缘线段分别进行精确定位，记录定位后 边缘线段的中心点坐标和其主要方向，可减弱了房屋提取结果误差的影响。
【附图说明】
[0015] 图1为本发明所述的一种针对多元遥感数据的房屋矢量化方法的技术路线图； 图2为本发明所述的一种针对多元遥感数据的房屋矢量化方法的流程示意图； 图3为本发明【具体实施方式】步骤1中获取房屋的示意图； 图4为本发明【具体实施方式】步骤2中获取房屋边缘的示意图； 图5为本发明【具体实施方式】步骤3中获取房屋主要方向的示意图； 图6为本发明【具体实施方式】步骤5中边缘模板构建的示意图； 图7为本发明【具体实施方式】步骤5中截取对应局部图像的示意图； 图8为本发明【具体实施方式】中房屋边缘线多类分割后的示意图； 图9为本发明【具体实施方式】中房屋矢量化结果的示意图； 图10为本发明【具体实施方式】中房屋矢量化结果与实拍影像叠加的示意图。
【具体实施方式】
[0016] 下面将结合具体的实施例来说明本发明的内容。
[0017] 如图1和图2所示，本实施例所述一种针对多元遥感数据的房屋矢量化方法具体流 程如下： 步骤1，获取需要处理的房屋。
[0018] 如图3所示，此处所述的房屋可以是采用各种手段得到的房屋候选点，也可以是房 屋斑块、房屋边界线等可以代表房屋位置、形状的信息:例如利用多尺度分割结合支持向量 机分割得到的房屋斑块、利用点云分割技术和α-shapes算法得到的房屋点集外轮廓线、利 用活动轮廓方法得到的房屋边缘线等。
[0019] 步骤2,获取房屋边缘。
[0020] 如图4所示，此处主要进行两步：一是采用八邻域边缘追踪得到有序的边缘点，主 要用于后续的步骤4、步骤5和步骤6;二是生成房屋边缘二值图像，边缘点为"1"，非边缘点 为 "0"。
[0021]步骤3,如图5所示，利用Radon变换结合主轴分析得到房屋的两个相互正交的主要 方向。
[0022] 具体步骤如下： 31)对初始轮廓线图像进行Radon变换得到一个累加数列，投影角度Θ为0°~179° ;累加 数列中，行表示投影角度Θ，列表示投影累计值。
[0023] 32)去除每一列中小于2的累计值(包含最大值），以减少计算量。
[0024] 33)将每一列的最大值和次大值相加作为这一列的的最终值，生成一个1x180的数 组。
[0025] 34)将数组分成两组:一组角度范围是0°~89°，另一组角度范围是90°~179° ;接 着将第二部分的投影累计值与对应的第一部分投影累计值相加，即0°+90°，1°+91°，……， 89°+179°，从而得到第一部分0°~89°新的角度累计值。
[0026] 35)最后依据角度累加值对角度θ(〇°~89°)进行重新排序，取最大值0max作为第一 主方向，此时第二主方向为9max+90°。
[0027] 步骤4,对房屋的边缘线进行边缘线段的多类分割。
[0028] 此处所述边缘线段多类分割是指将房屋的边缘线分割成三类边缘线段:第一主方 向类、第二主方向类、不定向类。
[0029] 具体步骤如下： 41)计算每个边缘点的局部方向：采用以当前点P为中心，取依次取与其相邻的半径R范 围内的轮廓点，即2R+1个点进行主成分分析，以第一主成分作为当前点P的法向量(a，b)，计 算当前点的局部方向。
[0030] local a-i?,e￡r{|j)- arct3ti{ls/a) (1) 42)计算连接每一个边缘点分别与两个顶点（θ，θ+90)的T链的权重。
其中，__和分别为点ρ属于第一主方向g类和第二主方向料雜类需要付出 的代价。系数科为断开T链的一个惩罚系数。参数为常数。
[0032] 43)计算连接点p与不确定类的T链的权重为：
其中，由于这一类的方向是不确定的，我们以一个固定值知作为断开点P与不确定类连 接需要付出的代价。
[0033] 44)计算相连边缘点之间N链的权重。
[0034] 根据相邻边缘点趋于同一类的这一先验知识，如果相邻边缘点被分为同一类，这 两点之间的N链并没有被破坏，那么N链的权重为0;如果相邻边缘点未被分为同一类，则必 须进行惩罚，惩罚值与这两个点的局部方向密切相关，相应的计算公式如下：
其中，系数4为相邻边缘点没有被分为同一类的惩罚系数，也是引入辅助节点的惩罚 系数。Wiyk)即为将相邻边缘点分成不同类别需要付出的代价。
[0035] 45)边缘线段分割的能量函数形式为：
式中，S表示一栋房屋所有边缘点的集合，N为边缘点的邻域系统，本文均采用八邻域。 ?是边缘线的一种分割结果。为分割结果中边缘点P所属类别，有三个标记值："Γ、"2"、 "3"，分别表示第一主方向类、第二主方向类、不定方向类。
[0036] 在构建好能量函数的基础上，利用α-扩展算法对能量函数进行解算以将房屋边缘 线分割为三类边缘线段。
[0037] 步骤5,利用模板匹配方法，对边缘线段进行精确定位。
[0038] 此处所述边缘线段是指在边缘线多类分割中被分为第一主方向类和第二主方向 类的边缘线段。因为不定方法类主要由锯齿状噪声、拐角点及其周围相邻点组成，因此只对 第一主方向类和第二主方向类进行精确定位。由于同一边缘线段中各边缘点出现偏差情况 的不同，采用对每一条边缘线段(第一主方向类和第二主方向类)上的所有边缘点在其法方 向上进行模板匹配，寻找最佳匹配点来完成对该边缘线段的精确定位。
[0039] 详细步骤： 51)边缘线段的选取：按照顺时针方向，依次选取每一条边缘线段％16廣。深为房屋边 缘线中第一主方向类和第二主方向类的边缘线段数目之和，假设边缘线段的方向为t 52)边缘模板的构建:如图6所示，首先构建初始边缘模板Mask;然后旋转模板至模板中 边缘方向为?，且II区域对应背景，I区域对应房屋。
[0040] 53)从第一个边缘点开始，依次选取每一个边缘点； 54 )如图7所示，将边缘模板Mask的中心点对准候选边缘点，截取对应的局部图像，计算 并记录模板与局部图像之间的相关系数P以及差值離。图6所示模板中，差值_可表示为：
式中，分别表示与模板I和II区域对应的局部图像区域的像素个数。/是图像窗 口。f和叙分别表示图像窗口中与模板I和II区域对应的像素集合。
[0041 ] 55)将模板对准候选边缘点，沿边缘线段的法方向（θ+90° )，在半径R(R=〖|j:，如SR= 〇.13，R=8)的范围内进行平移，同时计算模板与相应局部图像之间的相关系数:P以及差值 鮮。
[0042] 56)记录最佳匹配窗口中心的X，Y坐标，相关系数，差值，边缘线方向，边缘线所属 类别。最佳匹配窗口具有最大的差值_銳_，且其相关系数:P:大于指定阈值办（例如，办_ __=: 0.5)。如果P < Pr或者A/n <心/r，则认为未找到最佳匹配窗口，删除此时记录的信息。 [0043] 57)返回步骤53)，重复进行，直到处理完所有的候选边缘点为止。
[0044] 58)判断边缘线段是否应删除。在进行模板匹配的过程中，假设一条边缘线段上的 点全部被移除，也就是均没有找到符合条件的匹配点，如果该线段同时满足以下三个条件， 则保留其内原始点坐标： 该边缘线段的长度大于^，（即1.5米）； (D该边缘线段与相邻的边缘线段属于不同类； 急)其相邻的两个边缘线段属于同一类。
[0045] 59)返回步骤51)，选取下一条边缘线段，重复进行，直到处理完所有的边缘线段为 止。
[0046]步骤6,计算拐角点。
[0047] 在获取准确边缘线段上边缘点的基础上，计算相邻边缘线段的拐角点。
[0048] 61)针对所有定位后边缘线段，依次选取每一条边缘线； 62 )计算定位后边缘线段的重心（)。以边缘线段所属方向奋和重心（知％ )表示当 前边缘线段，记录直线参数(3???) i9); 63) 计算边缘线段的长度藏，!为图像空间中的线段长度，SR为影像的空间分 辨率； 64) 若￡签，则返回步骤61)。否则继续进行，计算每个点？Kx，y)到相应线段的距离色：
如果_方也(1=[0.5/SR]，则删除此时的边缘点(X，y)。
[0049] 65)返回步骤62)，直到所有边缘点满足Dis玄d ; 66)返回步骤61)，重复进行，直到处理完所有边缘线段为止。
[0050] 67)计算相邻边界线段之间的拐角点（威，辦），相邻两条边界线段的直线参数分别 为(??，_)，（褐鑛幽:，_:M):
其中，参数S是为保证分母不为零。
[0051] 步骤7，生成矢量化的房屋边界线。
[0052] 71)依次连接每个拐角点，获取所有边界点(《 3 )。点S1 (綠，議），S2(χ·3，_)之间 直线参数的计算公式如下：
其中，点(茶30是点S1和S2之间的房屋边缘点。
[0053] 72)坐标转换。由于实验过程是在数字图像坐标系中进行的，矢量化的房屋边界需 要重新赋予坐标。已知实验图的左上角（χ，γ)，所有房屋边界线转换后的坐标(氧彳υ为： X = X: + X - ?;爹=_v + F ~.:i。
[0054] 本实施例中房屋边缘线多类分割后的示意图、房屋矢量化结果的示意图、以及房 屋矢量化结果与实拍影像叠加的示意图可见附图8至10。
[0055] 以上内容是结合具体的实施例对本发明所作的详细说明，不能认定本发明具体实 施仅限于这些说明。对于本发明所属技术领域的技术人员来说，在不脱离本发明构思的前 提下，还可以做出若干简单推演或替换，都应当视为属于本发明保护的范围。
【主权项】
1. 一种针对多元遥感数据的房屋矢量化方法，其特征是，包括W下步骤： 步骤1，获取需要处理的可W代表房屋位置、形状的信息； 步骤2,获取房屋边缘:采用八邻域边缘追踪得到有序的边缘点，同时生成房屋边缘二 值图像，边缘点为"Γ，非边缘点为"0"; 步骤3,利用Radon变换结合主轴分析得到房屋的两个相互正交的主要方向； 步骤4,将房屋的边缘线分割成Ξ类边缘线段：第一主方向类、第二主方向类、不定向 类； 步骤5,对第一主方向类和第二主方向类上的所有边缘点在其法方向上进行模板匹配， 寻找最佳匹配点来完成对第一主方向类和第二主方向类的边缘线段的精确定位； 步骤6,在准确获取第一主方向类和第二主方向类的边缘线段上边缘点的基础上，计算 相邻边缘线段的拐角点； 步骤7，生成矢量化的房屋边界线。2. 如权利要求1所述的一种针对多元遥感数据的房屋矢量化方法，其特征是，步骤3包 括W下步骤： 31) 对初始轮廓线图像进行Radon变换得到一个累加数列，投影角度目为0°~179%累加 数列中，行表示投影角度Θ，列表示投影累计值； 32) 去除每一列中小于2的包含最大值的累计值，W减少计算量； 33) 将每一列的最大值和次大值相加作为运一列的的最终值，生成一个1x180的数组； 34) 将数组分成两组：一组角度范围是0°~89°，另一组角度范围是90°~179%接着将 第二部分的投影累计值与对应的第一部分投影累计值相加，即0° +90°，1° +91°，……，89° + 179°，从而得到第一部分0°~89°新的角度累计值； 35) 最后依据角度累加值对角度θ(〇°~89°)进行重新排序，取最大值作为第一主方 向，此时第二主方向为0max+9〇°。3. 如权利要求2所述的一种针对多元遥感数据的房屋矢量化方法，其特征是，步骤4包 括W下步骤： 41) 计算每个边缘点的局部方向：采用W当前点P为中屯、，取依次取与其相邻的半径R范 围内的轮廓点，即2R+1个点进行主成分分析，W第一主成分作为当前点P的法向量(a,b)，计 算当前点的局部方向： 扭C圧!為濤玲磅均释麵I ; 42) 计算连接每一个边缘点分别与两个顶点（θ，θ+90)的T链的权重：其中，兹灣郝務:梦争攤分别为点Ρ属于第一主方向i类和第二主方向挣辜得扭类需要付出 的代价，系数私为断开T链的一个惩罚系数，参数醉?为常数； 43) 计算连接点Ρ与不确定类的Τ链的权重为：其中，由于方向不确定因而W固定值.?作为断开点P与不确定类连接需要付出的代价； 44) 计算相连边缘点之间N链的权重：根据相邻边缘点趋于同一类的运一先验知识，如 果相邻边缘点被分为同一类，运两点之间的N链并没有被破坏，那么N链的权重为0;如果相 邻边缘点未被分为同一类，则必须进行惩罚，惩罚值与运两个点的局部方向密切相关，相应 的计算公式如下：其中，系数.Ss为相邻边缘点没有被分为同一类的惩罚系数，也是引入辅助节点的惩罚系 数，嗎動縣即为将相邻边缘点分成不同类别需要付出的代价； 45) 边缘线段分割的能量函数形式为：式中，S表示一栋房屋所有边缘点的集合，N为边缘点的邻域系统，本文均采用八邻域，/ 是边缘线的一种分割结果，疾为分割结果中边缘点P所属类别，有Ξ个标记值："Γ、"2"、 "沪，分别表示第一主方向类、第二主方向类、不定方向类； 在构建好能量函数的基础上，利用α-扩展算法对能量函数进行解算W将房屋边缘线分 割为Ξ类边缘线段。4.如权利要求3所述的一种针对多元遥感数据的房屋矢量化方法，其特征是，步骤5包 括W下步骤： 51) 边缘线段的选取:按照顺时针方向，依次选取每一条边缘线段i;-e <￥，旋为房屋边缘 线中第一主方向类和第二主方向类的边缘线段数目之和，假设边缘线段的方向为奋； 52) 边缘模板的构建:首先构建初始边缘模板Mask;然后旋转模板至模板中边缘方向为 @，且II区域对应背景，I区域对应房屋； 53) 从第一个边缘点开始，依次选取每一个边缘点； 54) 将边缘模板Mask的中屯、点对准候选边缘点，截取对应的局部图像，计算并记录模板 与局部图像之间的相关系数及差值幾f，差值褒/可表示为：式中，幾、驅分别表示与模板I和II区域对应的局部图像区域的像素个数，/是图像窗 口，《和猿分别表示图像窗口中与模板I和II区域对应的像素集合； 55) 将模板对准候选边缘点，沿边缘线段的法方向（θ+90° )，在半径R的范围内进行平 移，其中R=[fl，同时计算模板与相应局部图像之间的相关系数Ρ助及差值避f; 56) 记录最佳匹配窗口中屯、的Χ，Υ坐标，相关系数，差值，边缘线方向，边缘线所属类别， 最佳匹配窗口具有最大的差值且其相关系数Ρ大于指定阔值熱，如果肖《綻或者 <取，则认为未找到最佳匹配窗口，删除此时记录的信息； 57) 返回步骤53)，重复进行，直到处理完所有的候选边缘点为止； 58) 判断边缘线段是否应删除:在进行模板匹配的过程中，假设一条边缘线段上的点全 部被移除，也就是均没有找到符合条件的匹配点，如果该线段同时满足W下Ξ个条件，则保 留其内原始点坐标： 惠该边缘线段的长度大于i;;; 望)该边缘线段与相邻的边缘线段属于不同类； 堇其相邻的两个边缘线段属于同一类； 59)返回步骤51)，选取下一条边缘线段，重复进行，直到处理完所有的边缘线段为止。5. 如权利要求4所述的一种针对多元遥感数据的房屋矢量化方法，其特征是，步骤6包 括W下步骤： 61) 针对所有定位后边缘线段，依次选取每一条边缘线； 62) 计算定位后边缘线段的重屯、(?.鼓），W边缘线段所属方向愚和重屯、0:。>沉)表示当前 边缘线段，记录直线参数(％，:翁> S); 63) 计算边缘线段的长度￡，&. = !我繊，I为图像空间中的线段长度，SR为影像的空间分辨 率. 64) 若抗篮1濃，则返回步骤61)，否则继续进行，计算每个点g(x，y)到相应线段的距离凌：如果殺芬运，d=[0.5/SR]，则删除此时的边缘点(X，y); 65) 返回步骤62)，直到所有边缘点满足泣S立-d ; 66 )返回步骤61 )，重复进行，直到处理完所有边缘线段为止； 67)计算相邻边界线段之间的拐角点（X!':，於），相邻两条边界线段的直线参数分别为（ 斬躬.，馬），（灰編帶卒吉，瑟編）：其中，参数S是为保证分母不为零。6. 如权利要求5所述的一种针对多元遥感数据的房屋矢量化方法，其特征是，步骤7包 括W下步骤： 71) 依次连接每个拐角点，获取所有边界点（％方)，点S1(汹I，斜)，S2(紛，縛)之间直线 参数的计算公式如下：其中，点巧茶)是点S1和S2之间的房屋边缘点； 72) 坐标转换：由于实验过程是在数字图像坐标系中进行的，矢量化的房屋边界需要重 新赋予坐标，已知实验图的左上角α，Υ)，所有房屋边界线转换后的坐标（?,5Μ为： .?. = .?-斗方一主；X = V 冬 F - 1。
【文档编号】G06K9/00GK106096497SQ201610360575
【公开日】2016年11月9日
【申请日】2016年5月28日
【发明人】孙金彦, 王春林, 钱海明, 黄祚继, 周杰, 宋强
【申请人】安徽省（水利部淮河水利委员会）水利科学研究院