建筑物轮廓方向计算和规则化的向量重组算法

刘昌振,马 威,马 红,魏世轩

1. 重庆市测绘科学技术研究院,重庆 401121; 2. 自然资源部智能城市时空信息与装备工程技术创新中心,重庆 401121; 3. 自然资源部国土空间规划监测评估预警重点实验室,重庆 401147; 4. 重庆交通大学智慧城市学院,重庆 400074

建筑物作为构筑城市的主体,是重要的地理空间要素。利用高分辨率遥感影像[1]、LiDAR点云[2]等各类数据提取建筑物轮廓是摄影测量与遥感领域的研究热点,也是地图制图[3]、智慧城市建设[4]等应用的重要数据源。然而,在原始数据获取和建筑轮廓提取过程中往往存在变形和误差,导致其几何特征无法满足相应的规则和应用要求。因此,在建筑物轮廓提取后,需要进一步对其几何进行规则化处理,构建其规则特征[5-6]。

目前,规则化方法主要分为主方向法和最小二乘法。主方向法以计算出的建筑物主方向为约束,将满足条件的轮廓边进行旋转,使其和主方向垂直或平行来实现直角化。常用的建筑物主方向计算方法有最小外接矩形法[7-10]、统计加权法[11-13]和墙均值法[14-15]。其中,最小外接矩形法是通过计算建筑物的最小外接矩形,并将其作为建筑物的方向。如文献[7]使用最小外接矩形估计概略主方向,然后利用Hough变换检测主方向,从而提取建筑物轮廓边。文献[10]以最小外接矩形为主方向,通过栅格填充的方式实现规则化。统计加权法是先确定一系列候选方向,然后计算建筑物各边对每个候选方向的贡献值,取贡献值最大的方向作为建筑物的方向。文献[2]通过加权平均值计算建筑物主方向,将轮廓边分为平行和垂直两类,以轮廓边中心点为轴心,旋转轮廓边实现规则化。墙均值法是以建筑物轮廓每条边方向与权重(即长度)乘积之和的平均值作为建筑物的方向。文献[16]使用墙均值法加强开放街道地图中建筑物的直角特征。另外,文献[17]结合Radon变换和主轴分析确定建筑物的主方向,然后将轮廓边分类并进行精确定位实现规则化建筑物轮廓。文献[18]使用主方向法优化遥感影像提取的建筑物轮廓。基于主方向法的建筑轮廓提取算法相对简单,且计算效率高,但现有方法在计算过程中通常只能得到一个方向结果,无法应对具有多个方向特征的建筑。

最小二乘法是以建筑物轮廓节点或轮廓边上的点作为观测值,将建筑物轮廓相邻边的垂直关系作为限制条件,使用总体最小二乘方法求解,获得直角化后满足限制条件的节点坐标,从而实现建筑物轮廓的直角化[6,19-20],其改进方法有混合最小二乘-总体最小二乘直接解法[21]、加权总体最小二乘法[22]和限制条件偏最小二乘法[23]。此外,也有研究将最小二乘方法应用在点云数据提取建筑物轮廓规则化[24-25]、室内地图边界校正[26]等方面。相比主方向法,基于最小二乘法的直角化结果误差更小,但也存在建模、推导和解算困难的问题[21]。

现有主方向法和最小二乘法解决了具有一个主方向的简单建筑物直角化问题,然而对于复杂建筑物及建筑群的规则化尚没有较好的应对方法。建筑物的规则特征主要表现为建筑物的方向及相邻轮廓边的夹角,当建筑物轮廓边的夹角全部为直角,则轮廓边的方向可由两个垂直的方向表示。然而,现实中许多建筑物与其相邻轮廓边夹角中存在非直角的情况,需要用多个方向来表达。目前,相关的研究尚缺乏对这种复杂建筑轮廓规则化的解决方法。此外,现有方法一般针对单个建筑物进行规则化,未考虑相邻建筑之间的群体关系,导致处理建筑群时整体规则化效果不佳。

针对上述问题,本文提出一种基于向量重组的建筑物轮廓方向计算和规则化方法。首先将建筑物的轮廓边转化为一组按照一定顺序组成的向量集合,根据向量的方向特征进行分组和变换,对向量组计算方向向量,突破现有方法只能描述建筑物单一主方向的限制。在此基础上,使用直角改正加强方向的直角特征,进一步修正建筑规则化的几何特征。最后,以方向向量为约束条件调整建筑物轮廓边,实现建筑物轮廓的规则化。

1 基于向量重组的规则化算法

1.1 方法概述

基于向量重组的建筑物轮廓规则化算法流程如图1所示,主要分为2个部分:方向向量计算和规则化。首先输入建筑物原始轮廓,将其转换为一组按照一定顺序组成的向量集合,并计算得到初始向量组。然后,根据向量的方向特征进行分组和变换,将每组向量描述为一个方向特征,进一步通过直角改正加强方向的直角特征,从而得到描述建筑物轮廓方向特征的方向向量。最后,以方向向量为约束条件,计算轮廓边旋转基点,对建筑物轮廓边进行调整,得到建筑物规则化结果。

图1 建筑物轮廓方向计算和规则化方法流程Fig.1 Flowchart of building footprints orientation calculation and regularization

1.2 基于向量重组的建筑物轮廓方向向量计算

基于向量重组的建筑物轮廓方向计算方法主要思想是将建筑物轮廓边视为向量,并设计一定规则对每个建筑物轮廓边进行重组。然后,依据方向特征对重组边进行分组,将同一组向量边变换为同一方向。进而计算每个分组方向,判断分组方向是否满足垂直关系,如果满足则进行直角改正,从而获得描述建筑物轮廓方向的一组方向向量。在向量重组和分组的过程中,相对平行的边分在同一组,计算得到的方向向量表示了同一方向的所有平行边,因此本文提出的向量重组方法得到的方向向量可以从整体上描述建筑物轮廓的方向特征。

1.2.1 向量重组

向量重组是将建筑物的轮廓边转化为向量,则建筑物轮廓可以看作为按照一定顺序组成的向量集合[27],如图2(a)所示,建筑多边形由P1,P2,…,P8共8个点组成,将多边形的边表示为向量,则多边形由向量P1P2,P2P3,…,P8P1组成的集合。然后,将所有向量的起点平移至同一点,因原始建筑物轮廓的规则特征,重组后的向量具有一定的方向特征。如图2(b)所示,将向量的起点平移到同一个点O,组成新的向量,如向量P1P2重组为向量OP2,向量P2P3重组为向量OP3。

图2 向量重组Fig.2 Vectors reconstruction

1.2.2 向量分组和变换

通过向量重组后,因原始建筑物的规则特征,重组边具有明显的方向特征。考虑到原始数据或提取误差的存在,方向接近向量边的方向值不完全相等。因此需设定条件阈值,将方位角小于条件阈值或旋转角度π后小于条件阈值的边分为同一组。根据建筑物的复杂程度,可得到多组向量,每一组向量表示建筑物轮廓的一个方向。如图3所示,重组后向量方位角取值为[0,2π),根据向量方位角进行分组,分组规则如下:设定一个角度阈值θ,向量的方位角α1、α2,若|α1-α2|<θ、|α1-α2-π|<θ或|α2-α1-π|<θ,则将向量分为一组。θ的取值和原始数据精度有较大的相关性,误差越大,则同一方向边的角度差值越大,相应θ的取值应大一些,向量重组后,同一方向的向量有一定的差值,角度阈值θ可依据该差值确定。根据分组规则及阈值设定,图3中的原始向量可分为3组:{OP2,OP5,OP8}、{OP1,OP4,OP6}、{OP3,OP7}。对分组后向量的进行变换,使同一组的向量方向相同,如第一组中,OP5、OP8的方向和OP2的方向相反,对OP5、OP8旋转角度π,分别得到OP5′、OP8′。对其余向量用同样的方法变换,最终向量分为3组:G1={OP2,OP5′,OP8′},G2={OP1,OP4′,OP6},G3={OP3′,OP7}。

图3 向量分组和变换Fig.3 Vector grouping and transformation

1.2.3 分组方向计算

利用向量重组方法获得每一组的向量边后,需要计算每一组的方向值代表该组的整体方向。合理的分组方向值要考虑每个向量边的长度和方向,向量边的长度不同,对分组方向的贡献应该也不同[14-15]。本文采用总体最小二乘方法获得每个分组的方向值:计算每个向量边终点到分组方向的距离,当距离的平方和最小时,则认为该方向为分组方向。

如图4所示,向量α1,α2,…,αn为一组向量,起点均为点O,假设其分组方向经过点p0,则其向量可表示为p0,向量v为分组方向p0的单位法向量。当分组方向满足总体最小二乘时,向量α1-p0,α2-p0,…,αn-p0在向量v投影长度的平方最小,即向量α1-p0,α2-p0,…,αn-p0和向量v点乘的平方和最小,可得目标函数

图4 分组方向计算Fig.4 Orientation calculation of group

(1)

表示为矩阵形式

(2)

对列向量p0求导可得

(3)

(4)

1.2.4 直角改正

原建筑物轮廓边本应垂直的方向向量夹角不一定为π/2,在获得每组方向向量后,仍需对分组方向进一步处理,使应垂直的方向向量严格垂直,本文称该过程为分组方向的直角改正。若两分组方向向量的方位角|α1-α2-π/2|<θ,则对其进行直角改正,对方位角较大的向量,顺时针旋转π/2,对旋转后的向量和方位角较小的向量求和,求和后向量的方向为直角改正后方位角较小向量的方向,逆时针旋转π/2后的方向为方位角较大的向量方向,分组向量长度取原向量在改正后方向向量的投影长度,最终得到直角改正后的分组向量。由式(1)—式(4)的推导可知,直角改正过程满足总体最小二乘。

对图5(a)所示的建筑物轮廓计算方向,计算分组向量和:OG1=OP2+OP5′+OP8′,OG2=OP1+OP4′+OP6,OG3=OP3′+OP7。分组向量OG1和OG2满足直角改正的条件,将向量OG2顺时针旋转π/2后得到向量OR2,对向量OG1和向量OR2求和,取求和后向量的方向和OG1在求和向量上的投影长度,得到OG1直角改正后的向量OGR1,将OGR1逆时针旋转π/2,并取OG2的在该旋转向量的投影长度,得到OG2直角改正后的向量OGR2,得到建筑物轮廓的最终方向向量:OGR1、OGR2和OG3。

图5 方向计算和直角改正Fig.5 Orientation calculation and right angle correction

1.3 建筑物轮廓规则化

在计算得到建筑物方向向量的基础上,构建建筑物规则化通用公式,基于方向向量的旋转基点进行节点坐标的调整,从而实现建筑物规则化。如图6所示,l为待调整边,l1和l2为其相邻边,l′为其通过向量重组计算出的方向,对建筑物轮廓规则化,就是求解一点P,将待调整边l绕点P旋转到和l′相同的方向,P为旋转基点。

图6 轮廓边旋转基点Fig.6 The rotation base point of footprints edge

S=(x1-x)2+[y1-(k1x+b1)]2+

(x2-x)2+[y2-(k2x+b2)]2

(5)

将y1=k1x1+b1和y2=k2x2+b2代入式(5),可得

(6)

对x求导可得

(7)

令导数等于0,可求得

(8)

设l′和l1的夹为α1,设l′和l2的夹为α2,l′的方位角为π/2,k1和α1的关系为k1=cotα1,k2和α2的关系为k2=cotα2,代入式(8)可得

(9)

式(9)表示出了求解点P横坐标的公式,以l′的方向为横坐标轴方向,同样的推导方法可得纵坐标公式

(10)

式(10)是假设l1和l2的斜率都存在情况下得到的,当斜率不存在时,夹角α=π/2,此时有cotα=0,代入式(10)即可得到斜率不存在时纵坐标的计算公式。

通过式(9)和式(10)即可求解出旋转基点P的坐标。式(9)和式(10)表明,当P1和P2坐标固定,则点P坐标只和相邻两条边的夹角有关。当坐标系发生旋转时,如图7所示,相邻边的夹角不会改变,点P相对于P1和P2的位置关系不会改变,式(9)和式(10)对旋转后的坐标系同样适用,因此式(9)和式(10)是解算旋转基点P的一般公式。

图7 坐标系旋转Fig.7 Coordinate system rotation

特别地,当α1=α2时,即待调整边和相邻两个边的夹角相同时,有

(11)

(12)

式(11)和式(12)表示,当待调整边和相邻两个边夹角相同时,待调整边的中点即为旋转基点。

图8 建筑物规则化方法Fig.8 Building regularization method

2 试验与分析

为验证本文方法的有效性和效率,本文首先对复杂建筑和建筑群进行规则化,然后和主方向法进行对比分析,最后对算法效率分析。

2.1 复杂建筑物规则化

现有方法难以实现具有多方向特征复杂建筑物的规则化,选取形状较复杂的建筑物轮廓进行方向计算和规则化试验。图9(a)为原始建筑轮廓经过向量重组变换后的重组边和方向向量对比(为便于展示,方向向量的长度进行了缩放),原始重组边的方向有一定的差别,经过方向计算和直角改正后,最终得到3组相互垂直的方向向量。以方向向量为控制条件,对建筑轮廓进行规则化处理,结果如图9(b)所示,试验结果证明方法对于复杂建筑同样适用,在加强建筑物规则特征的同时保持了原始形状特征。

图9 复杂建筑物方向向量和规则化结果Fig.9 Complex building orientation vectors and regularization results

2.2 建筑群规则化

邻近的建筑物方向通常具有相似性,尤其是多个建筑组成的建筑群,其中部分建筑物往往具有相同的方向。如图10(a)所示为多个建筑组成的建筑群和单个建筑物的方向向量,可将建筑群作为一个整体,得到建筑群的方向向量,结果如图10(b)所示,最终得到3组垂直的方向向量,建筑方向向量的角度值见表1。

表1 建筑群方向向量的角度值

图10 建筑群方向向量Fig.10 Orientation vectors of building group

从建筑方向向量中可得到建筑物群的方向关系,如建筑1、2、3具有近似的方向特征,然而,由于原始采集误差的存在,逐一计算单个建筑物方向会导致相同方向建筑物的方向向量存在一定的差别。将建筑群作为整体考虑,计算得到的方向向量更能代表建筑群的方向特征,以方向3为例,有方向3的建筑包括建筑4—建筑10共7个建筑,每个建筑方向3的具体方向值不相等,因此,考虑更多建筑的整体方向值更加准确。

以得到的整体方向向量为控制条件,对每个建筑轮廓进行规则化,规则化结果如图11所示。相比于对单个建筑规则化,以建筑群为整体得到的规则化结果,具有接近方向特征的建筑轮廓方向完全相同,因此规则化结果更加合理。

图11 建筑群规则化结果Fig.11 Regularization results of building group

2.3 对比分析

为进一步验证本文方法,将本文方法与主方向法进行对比分析。基于最小外接矩形法计算得到的方向较为稳健且应用较多[7-10,18],本文使用最小外接矩形获得建筑物主方向。主方向法只能实现建筑物轮廓的直角化,以建筑物轮廓的直角化来对比两种方法,在直角化的过程中,应尽量减少节点的移动距离,因此通过节点移动距离的均方根误差评价两种方法。试验建筑物轮廓有6个节点且相邻轮廓边相互垂直,节点坐标见表2,以该建筑物的节点坐标作为真实值,通过加入高斯误差的方式获得模拟数据。为充分对比两种方法,共进行1000次模拟试验,加入均值为0、标准差为0.5 m的高斯误差。加入高斯误差后,相邻轮廓边不再满足垂直条件,对模拟结果进行直角化处理,角度阈值取值为π/12,计算直角化前后节点移动距离的均方根误差。原始模拟结果节点均方根误差分布如图12(a)所示,为直观对比两种方法误差的分布规律,以主方向法均方根误差排序结果如图12(b)所示,以本文方法均方根误差排序结果如图12(c)所示。

表2 试验建筑物节点坐标

图12 主方向法和本文方法1000次模拟节点均方根误差Fig.12 Coordinate RMSE of main orientation method and the proposed method for 1000 simulations

在1000次模拟的直角化结果中,主方向法节点移动最大均方根误差为0.972 3 m,平均值为0.324 9 m,本文方法最大均方根误差为0.736 8 m,平均值为0.305 1 m。由图12(b)可知,以主方向法均方根误差从小到大排序,本文方法的误差曲线主要在主方向法误差曲线的下方波动;由图12(c)可知,以本文方法均方根误差从小到大排序,主方向的误差曲线在本文方法误差曲线的上方波动。

进一步和真实值比较,计算原始模拟数据、主方向法和本文方法直角化结果节点坐标相对于真实值节点坐标的均方根误差,原始模拟结果节点均方根误差分布如图13(a)所示,为展示原始数据和直角化前后误差的变化分布规律,对原始模拟数据误差排序结果如图13(b)所示。

图13 相对于建筑物节点真实值均方根误差Fig.13 RMSE relative to real building coordinates

在本文1000次模拟试验中,相对于真实数据,原始模拟数据节点均方根误差最大值为1.028 2 m,平均值为0.481 2 m;主方向法直角化后节点均方根误差最大值为0.892 1 m,平均值为0.385 1 m;本文方法直角化后节点均方根误差最大值为0.822 1 m,平均值为0.366 9 m。从1000次模拟结果节点的均方根误差平均值来看,主方向法将原始模拟数据的均方根误差减小了19.97%,本文方法将原始模拟数据的均方根误差减小了23.75%。进一步对误差分布图分析可以看出,主方向法直角化后大部分建筑的均方根误差得到了减小,但出现部分结果均方根误差增大的情况,在1000次模拟试验中,共出现80次,增大的最大值为0.159 6 m。本文方法直角化的均方根误差分布曲线主要在原始模拟数据均方根误差曲线的下方波动,只有极个别出现直角化后均方根误差增大的情况,在1000次模拟试验中,共出现2次,增大的最大值为0.001 2 m。因此,本文方法通过直角化方法对建筑轮廓直角化不仅可以从视觉上得到更好的结果,还可以进一步提高原始数据的整体精度。在精度提高和整体稳定性方面,本文方法优于主方向法。

2.4 算法效率分析

本文方法的时间复杂度主要与建筑物轮廓的方向个数有关,具体每个步骤的时间复杂度及说明见表3。

表3 算法的时间复杂度

试验对约16 km2的重庆市中心城区建筑物进行规则化,建筑物轮廓如图14所示。试验区建筑物轮廓数据取自1∶2000比例尺地形图,采用立体航测法采集,立体航片分辨率为0.2 m,建筑物采集精度为0.2 m,其中包含建筑物11 214个,共有139 658条边,平均每个建筑物12.5条边。测试计算机配置为Intel(R) Core(TM) i5-10400F CPU @2.90 GHz,内存8 GB,角度阈值取值为π/20。

图14 试验区建筑物轮廓Fig.14 Building footprints of experiment area

对每个建筑物计算方向向量和规则化,处理总时间为2 min 38 s,平均每个建筑物耗时为0.014 s,算法效率较高。建筑物方向向量个数统计见表4,方向向量总个数为31 290个,平均每个建筑物2.8个方向向量,方向向量为2个的建筑物共8111个,占建筑物总数的72.33%,方向向量总个数为4个及以下的占建筑物总数的88.67%。

表4 试验数据方向向量个数统计

3 总 结

本文针对复杂建筑物和建筑群轮廓的方向定量描述和规则化等问题,从建筑物轮廓的几何规则特征分析出发,提出了一种基于向量重组的建筑物轮廓方向计算和规则化方法。在方向计算方面,通过向量重组获得描述建筑物轮廓方向特征的重组边,经过向量变换和方向计算,实现了复杂建筑物方向的定量描述。在此基础上,以方向向量为控制条件,通过计算旋转基点,进一步调整节点坐标,实现了复杂建筑物和建筑群的规则化。

通过推导验证,本文方法在方向计算和规则化过程中均满足总体最小二乘,有效降低了直接构建总体最小二乘模型的复杂解算困难。通过与主方向法对比试验可知,本文提出的向量重组方法在节点整体移动距离和稳定性方面更具有优势,通过规则化算法可进一步提高数据的整体精度。通过复杂建筑物和建筑群的试验,进一步验证了本文方法在复杂建筑物和建筑群规则化的适用性,为建筑群空间分布特征识别、智慧城市构建、城市建筑变化检测等应用提供一种解决思路。