大规模城市可计算模型快速构建方法研究

王 淳,于长华,唐 昊,李海峰

(中国工程物理研究院计算机应用研究所，四川 绵阳 621999)

1 引言

数字城市是指以计算机技术、多媒体技术和大规模存储技术为基础，以宽带网络为纽带，运用遥感、GPS、GIS、遥测、仿真-虚拟等技术，对城市进行多分辨率、多尺度、多时空和多种类的三维描述[1]。随着智慧城市建设的开展，数字城市模型在城市规划与管理、环境模拟、应急响应、导航、虚拟旅游以及其他科学和大众应用方面的需求越来越高[2]，是提升城市竞争力、解决城市发展问题的重要支撑。

作为贯穿整个智慧城市建设的重要技术，数字城市建模是计算机视觉、计算机图形学、遥感测量等领域的热点问题，获得国内外研究者的大量关注[3]。从数据源来说，使用的数据包括从地面、航空和航天平台获取的影像和激光点云数据、众源地理数据、数字高程等。从建模方法来说，目前主流方法有多视图三维重建、基于激光点云的三维重建、基于矢量数据建模和多源数据融合建模等。从发展趋势来看，大规模数字城市建模正在朝自动化、精细化、语义化、集成化、标准化和开放共享等方向发展，在建模成本不断降低的同时增强建模效率和模型可重用性[4]。

本文关注的是数字城市可计算模型，即可用于进行数值模拟计算的数字城市模型。城市可计算模型扩展了数字城市的应用领域，使城市模型可用于自然灾害、战争破坏效应研究和环境模拟研究[5 - 7]，在解决重大领域问题、提高管理和决策的科学性和高效性等方面发挥着重大作用。

Figure 1 Diagram of city modeling using OpenStreetMap data图1 基于OpenStreetMap矢量数据的城市建模流程

为了给数值模拟软件提供大规模城市的可计算模型，不仅需要具有面向大空间尺度、海量数据、频繁更新的城市的快速建模能力，还需要使模型满足特定的可计算条件，如不允许存在实体相交、悬边、悬面、非闭合体、二次曲面等，且对不影响计算结果的细小特征进行化简。但是，由于数据采集过程中存在的误差，以及模型处理所做的必要的简化操作，导致直接生成的数字城市模型中存在大量非闭合、相接、相交、悬空等CAD模型缺陷[8,9]，使得城市模型不可用于数值模拟。由于城市模型的数据量巨大，不可能通过人工方式修复有缺陷的几何体。即使使用商用CAD软件的修复功能，图形的布尔运算仍需要耗费大量时间，无法满足快速构建大规模城市可计算模型的需要。

本文以矢量数据建模技术为基础，提出一种包含建筑和地形的LoD1城市模型[10]的快速构建方法，重点研究建模数据的可计算处理。通过分析LoD1城市模型中所包含的CAD模型缺陷，本文给出了几何图形需满足的可计算条件以及处理方法。这种在几何图形上进行可计算处理的方式避免了人工操作和CAD布尔运算，使得本文方法在建模效率、模型规模和可计算性等方面均能满足高性能计算的要求。

2 城市可计算建模流程

本文采用OpenStreetMap[11]矢量数据作为构建数字城市模型的数据源。OpenStreetMap是一个开源地理信息项目，通过标签和属性值提供结构化的地图描述，在路径规划、地图导航、智能交通等领域都得到广泛应用。目前，已有若干基于OpenStreetMap构建或浏览三维模型的软件，如blender-osm[12]、OSM2World[13]和OSMBuildings[14]等。

基于OpenStreetMap建立三维模型的一般流程参见文献[15]。本文的目标是可计算建模，需要处理原始数据中的非闭合、相接、相交、悬空等CAD模型缺陷，因此在一般流程的基础上还要加入“可计算处理”步骤。因此，本文采用的基于OpenStreetMap矢量数据构建城市模型的流程包含4个步骤(如图1所示)：数据读取、数据解析、可计算处理和模型输出。每个步骤的具体内容如下所示：

(1)数据读取。读取原始数据，即从OpenStreetMap数据中读取基本元素，从高程数据中读取顶点坐标。

(2)数据解析。将原始数据转化为可用于建模的城市单元体数据，包括解析建筑数据(含轮廓、高度、组合方式及其它属性)和地形三角网。

(3)可计算处理。消除城市单元体中的CAD模型缺陷，使其能够用于生成可计算城市模型。即消除建筑与建筑、建筑与地形之间的重叠、缝隙、曲率突变等几何错误和细小特征，满足数值模拟的要求。

(4)模型输出。将各类城市单元体合并，形成城市模型，并通过图形引擎输出。

在上述各个步骤中，数据读取、数据解析和模型输出与一般的建模软件相同。因此，本文重点研究的内容是可计算处理。

3 城市模型CAD模型缺陷分析

本文构建的是LoD1城市模型(如图2a所示)[10]，包含从OpenStreetMap数据中提取的建筑模型和从高程数据中提取的地形模型。建筑模型只包含地面轮廓和高度，高程地形由三角面片组成。LoD1城市模型是对真实城市的高度简化和抽象，既具有结构简单的优点，又能反映真实城市的几何特性，适合于大多数在大规模数字城市上计算的数值模拟程序。

Figure 2 Model structure图2 模型结构

由于数据采集过程中存在的误差，以及模型处理所做的必要的简化操作，生成的数字城市模型很可能存在CAD模型缺陷[8]。直接修复CAD模型缺陷的过程非常复杂和耗时，但考虑到LoD1城市模型(如图2a所示)的结构远比CAD模型(如图2b所示)的几何结构简单，如果在LoD1城市模型上进行可计算处理，将极大简化处理过程。这是由于:(1)LoD1城市模型只涉及一部分CAD模型缺陷；(2)对LoD1城市模型的可计算处理可转化为对二维图形的处理，将大大降低处理难度。基于以上思想，本文首先在LoD1城市模型所对应的二维图形上进行可计算处理，然后将其转化成CAD模型，就可保证大规模城市模型的可计算性和建模效率。

本节提出对LoD1城市模型进行可计算处理的解决方案。首先，给出CAD模型缺陷的分类和分析(3.1节)；然后，研究LoD1城市模型中存在的CAD模型缺陷及其应对措施(3.2节)；最后，给出使模型避免CAD模型缺陷、满足可计算条件的技术路线(3.3节)。

3.1 CAD模型缺陷分类

CAD模型的表达采用的是拓扑和几何分离的思想。几何数据描述空间实体的位置，拓扑结构数据描述空间实体间的连接关系。在这种分离的设计方式中，其几何数据与拓扑结构数据相对独立实现，然后通过关联操作建立它们之间的联系。

根据CAD模型拓扑和几何的特点，可以把CAD模型缺陷分为拓扑表达缺陷和几何表达缺陷，如图3所示。其中拓扑表达缺陷可以进一步分为非正则拓扑缺陷、拓扑异构缺陷和拓扑精度缺陷，几何表达缺陷可进一步分为几何异构缺陷、几何精度缺陷和几何潜在错误表达。详细的CAD模型缺陷分类及处理方法参见文献[8]。CAD模型缺陷严重破坏了几何模型的完整性和精确性，会导致数值模拟计算失败。

Figure 3 Categories of CAD model deficiencies图3 CAD模型缺陷分类

3.2 LoD1城市模型CAD模型缺陷分析

由于LoD1城市模型的特殊结构，它所对应的CAD模型只会涉及CAD模型缺陷的部分类型。对各类CAD模型缺陷在LoD1城市模型中的表现分析如下：

(1)非正则拓扑缺陷。CAD系统中模型可以是描述物理模型的正则形体，也可以是在现实世界不存在的非正则形体。当 CAD 系统是非正则形体时，一部分拓扑缺陷也被包含在模型中。典型的非正则拓扑缺陷有悬面、悬边和孤立的点。在LoD1城市模型中，如果建筑地面轮廓不是简单闭合的多边形，由其拉伸而成的建筑模型会产生非正则拓扑缺陷。另一种情况是建筑轮廓相交，此时顶点或边的拓扑连接关系会发生错误。

(2)拓扑异构缺陷。拓扑异构是指CAD模型中的模型拓扑关系在计算机中的描述方法不同，因此在文件转换时，会出现拓扑不相容等问题。本文的城市模型是基于图形引擎直接生成的，不涉及多个CAD系统间的转换，因此避免了拓扑异构缺陷。

(3)拓扑精度缺陷。拓扑精度缺陷是由系统精度问题造成的拓扑表达不一致，如曲面之间的裂缝或重叠、边与边的裂缝以及边与顶点的裂缝。对LoD1城市模型来说，当多个建筑模型相接，或当建筑模型与地形模型之间存在重叠或裂缝时，会产生拓扑精度缺陷。

(4)几何异构缺陷。几何异构是指异构CAD系统中描述模型实体的数学定义不同，主要发生在复杂曲线曲面(如NURBS曲线曲面)上。现有的建筑模型和地形模型只包含平面，因此避免了二次曲面带来的几何异构缺陷。

(5)几何精度缺陷。几何精度缺陷是由系统精度问题造成的几何表达不一致，主要是针对曲面和曲线而言。在LoD1城市模型中不存在这种情况。

(6)几何潜在错误。几何潜在错误主要由角度、尺寸、邻近和曲率4类情况产生的。大部分的潜在错误都是拓扑正确，并能产生合法的实体模型，因此本文不对几何潜在错误进行特别的处理。但是，为了减少数值模拟的计算量，需要消除模型中不影响计算结果的细小几何特征，这也会消除绝大部分几何潜在错误。

3.3 可计算处理技术路线

根据2.2节对CAD模型缺陷的分类和分析，可归纳出LoD1城市模型需满足的可计算条件：

(1)建筑底面轮廓为简单多边形，即顶点数大于或等于3且不存在自相交。

(2)建筑轮廓之间不存在相交。

(3)建筑轮廓之间不存在几何精度造成的相接。

(4)建筑模型底部与地形贴合，即既无裂缝也不重叠。

本文采用以下4类方法使城市模型满足可计算条件(如表1所示)：

Table 1 Analysis on CAD model deficiencies in LoD1 city model

(1)轮廓简化。由于原始建筑轮廓含有大量冗余点，对轮廓进行简化不仅可以减少冗余，还可以消除绝大部分几何潜在错误。

(2)轮廓收缩。建筑轮廓相交或相接是由数据采集误差造成的，将相交建筑的轮廓收缩一小段距离，就可以避免绝大多数相交或相接。

(3)地形贴合。将建筑向地形三角网投影，计算海拔高度，保证建筑各顶点的海拔高度一致，可避免建筑与地形之间的重叠或裂缝。

(4)删除。如果用以上方法仍不能满足可计算条件，直接将建筑删除，不对其建模。

基于以上分析，完整的城市模型可计算处理包括轮廓简化、自相交检测、相交检测、相接检测、地形贴合、高度检查、轮廓收缩等步骤。可计算处理流程如图4所示。

Figure 4 Computable processing pipeline in this paper图4 本文采用的可计算处理流程

4 城市模型可计算处理

本节给出图4所示的可计算处理流程中各个步骤的具体算法。

4.1 建筑轮廓自相交检测

Figure 5 Counter-clockwise polygon with self-intersection图5 逆时针走向多边形自相交

假设多边形顶点数是n，判断自相交将计算n条直线方程和n次直线相交，然后计算n个交点是否在多边形内部。由于判断点在多边形内部的算法复杂度是O(n)，该算法复杂度为O(n2)。

4.2 建筑轮廓简化

建筑轮廓简化是指对给定的建筑轮廓，在减少建筑轮廓顶点数量的情况下，使简化后的建筑轮廓与原建筑轮廓尽量保持形状一致。因为建筑轮廓只有一种拓扑类型，所以本文从目标边2条邻边的位置关系(平行、非平行)来考虑建筑轮廓简化问题。2种位置关系的简化如图6a和图6b所示，设简化容差为Ts,当Sn的长度小于Ts时，将Sn设置为待处理边。图中加粗的边为新增加的边，虚线边为待删除的边，文献[17]详细介绍了2种位置关系的简化方法。

Figure 6 Contour simplification method图6 轮廓简化方法

基于以上方法可实现建筑轮廓简化，在保证建筑模型在外观特征不失真的情况下，减少三维建筑物模型的细节特征和面片数，以适应不同复杂度要求的数字城市可计算模型。

4.3 建筑轮廓相交检测

如图7所示，2个简单多边形A和B的关系可分为6种情况：不相交、相接、叠加、包含、内部和相等[18]。对于可计算模型的建筑轮廓来说，只允许图7a和图7d所示的2种情况：

Figure 7 Spatial relations between polygons图7 多边形空间关系示意

(1)2栋建筑的底面轮廓不相交。此时只需判断A的所有边都在B的外部。当且仅当A的所有顶点都在B的外部，且A的每条边与B都不相交时，此情况成立。

(2)建筑底面轮廓包含空洞，即建筑的底面轮廓A包含空洞B。此时只需判断B的所有边都在A的内部。当且仅当B的每个顶点都在A的内部，且B的每条边都与A不相交时，此情况成立。

判断多边形关系的算法复杂度是O(n2)。

4.4 建筑轮廓相接检测

在CAD系统中，实体间的拓扑相邻关系经常出现误差精度导致的问题。CAD系统判断一个点是否位于一条直线上并不是由计算的绝对结果来决定，而是通过判断该点到直线的距离是否小于系统偏差。如果在系统允许的偏差范围内，则可以表示该点在曲线上，反之说明该点在曲线外。大多数商业CAD系统中使用的最小偏差大约为10-4甚至10-3。

Figure 8 Joint polygons caused by geometric precision图8 几何精度造成的相接

为了检测相接的建筑轮廓，需要计算一个多边形的所有顶点到另一个多边形所有边的最短距离。如果顶点不在边上但最短距离小于偏差(如图8所示)，为避免拓扑精度缺陷需判定为相接。相接检测的算法复杂度与相交检测相同，都是O(n2)。

4.5 建筑轮廓收缩

为避免城市建筑与相邻建筑碰撞，可将建筑轮廓收缩指定的距离(即安全缓冲距离)，即从建筑轮廓顶点沿角平分线方向向建筑轮廓内部平移一定的距离。

如图9所示，多边形的相邻2条边L1和L2，v1和v2分别是L1和L2的单位向量，2条边在交点Pi处作平行于L1和L2、平行线间距为L且位于多边形内部的2条边，交于Qi。Qi的计算公式如下所示：

其中,sinθ=v1×v2。对建筑轮廓的n个顶点进行收缩，算法的复杂度为O(n)。

Figure 9 Computation of compressed vertex图9 内缩点计算

4.6 建筑与地形贴合

实现建筑与地形贴合的关键在于计算建筑所在的地表海拔高度，即从建筑顶点发出垂直射线向地形三角网投影，并计算投影点所在的三角面片。

通过射线与三角网求交计算建筑所在地表的海拔高度的方法如下所示。将射线R定义为一个基点P和一个方向d，其参数方程是：

R(t)=P+td

假设三角形Q的顶点序列为{V0,V1,V2}，对于三角形内的任意一点，它的位置可通过三角形顶点定义，其参数方程如下：

Q(a,b,c)=aV0+bV1+cV2

其中，a+b+c=1。计算射线R与三角面片Q的交点只需将上述2个方程联立，经整理为：

这是一个由3个未知量组成的线性方程。求解t,b,c的值，如果0≤b≤1，0≤c≤1且b+c≤1，则交点位于三角形内。通过以上算法完成所有三角网面片求交测试，得到射线与三角网的交点。

如果建筑通过了高度检查，即建筑轮廓所有顶点的海拔高度一致且顶点在合法坐标范围内，则将投影点所在的三角面片的高度作为建筑海拔高度。

5 大规模图形快速检索

对于大规模城市模型来说，可计算处理的计算复杂度成为影响算法性能的关键问题。假设建筑数量为M，建筑顶点数量为N，地形三角面片数量为D，则对所有建筑计算相接或相交关系的复杂度为O(M2N2)，为所有建筑计算海拔高度的复杂度为O(MND)。大规模城市中通常建筑和地形三角面片的数量大约在105～106量级，如此规模的计算量实际上是无法执行的。为了解决建模效率的问题，本节提出大规模图形快速检索算法，以提高图形检索和处理的效率，降低可计算处理的复杂度。

5.1 邻近建筑轮廓快速检测

在对建筑模型的可计算处理中，建筑地面轮廓的相交和相接检测占用了大量时间。为此本文提出基于最大包围盒的邻近建筑轮廓快速检测算法。通过快速提取与每栋建筑邻近的建筑，减少相交检测次数。

如图10所示，黑点代表城市区域内建筑位置。首先计算所有建筑的包围盒，记所有包围盒的最大长度和宽度为(lmax,wmax)。假设当前计算的建筑(五角星所示)的包围盒长和宽为(l,w)，则x方向上只有与该建筑包围盒中心距离小于(l+lmax)/2的建筑可能与其相交(垂直带状灰色区域所示)。同理，y方向上与该建筑包围盒中心距离小于(w+wmax)/2的建筑可能与其相交(水平带状灰色区域所示)。因此，相交的建筑一定在x和y方向上同时满足包围盒中心距离要求，即落在灰色阴影区域。该区域内的建筑数量与城市建筑总量无关，取决于城市建筑的最大密度和最大包围盒范围，可以认为是一个常数。因此，一栋建筑只需要与常数个其它建筑作相交检测，就可得到与其相交的所有建筑。这样可将建筑相接和检测算法的复杂度降至O(MN2)。

Figure 10 Fast detection of neighboring building contours图10 邻近建筑轮廓快速检测

5.2 射线与三角网快速求交

由于城市模型中包含大量的建筑顶点和地形三角面片，而普通的求交涉及大量冗余的射线和三角面片，计算量过高，为此本文提出基于层次化包围盒的快速求交算法。

本文将整个三角网区域进行层次化分解，待计算出最小的包围盒后，再对其内部的三角网面片进行求交测试。该算法剔除了包围盒外部的面片，因此大大减少了求交测试的次数。具体算法步骤如下所示：

(1) 将三角网投影到xy平面上，构建所有三角网的二维矩形包围盒。

(2) 在矩形包围盒的基础上，按x和y2个方向将其分割成4个子矩形。

(3) 若某一子区域的三角形数量大于给定的阈值m，则返回(2)，将该子矩形进一步剖分成4个更小的矩形。直到每个矩形区域里的三角面片数都小于m(如图11所示)。

Figure 11 Recursive splitting to quadtree图11 四叉树递归分割

(4) 对于划分的每个子矩形，计算其包含的三角网的三维长方体包围盒，并组成包围盒树。

(5) 当计算射线与三角网交点时，首先从上至下计算射线与包围盒树的每个节点的交点。如果射线与包围盒不相交，则剔除包围盒内的所有三角面片，由此得到最小的包围盒区域，再将射线与区域块内的三角面片求交。

这样构造的层次化包围盒树的层数为O(logD)。此时计算一条射线与三角网的交点需要从顶向下计算射线与包围盒树的相交次数是O(logD)，计算所有建筑顶点海拔高度的复杂度降为O(MNlogD)。

对绝大多数建筑来说，其所有顶点都位于相同的地形三角面片内。因此，可采用启发式算法进一步降低复杂度。在获取与第1个建筑顶点相交的三角面片后，其它建筑顶点首先直接与该三角面片求交，如果不相交再与地形包围盒树求交。这将使绝大多数建筑顶点求交的复杂度降为O(1)，理想情况下所有建筑顶点海拔高度的复杂度为O(MlogD)。

6 实验

本文在服务器单个结点上构建可计算模型。每个结点包括144核，主频为2.3 GHz。对可计算处理的关键步骤(如建筑轮廓求交、简化、地形贴合)采用OpenMP并行。

6.1 模型可计算性验证

为了验证模型的可计算性，本文选择2 000 m×2 000 m范围的某城市区域，在同一区域构建了4个经过可计算处理(选择不同的简化边长度参数，Ts=0 m,4 m,8 m,12 m)的城市模型。本文设计了运输炸药车辆在城市街道发生意外突然爆炸的计算场景，利用自主研制程序，模拟计算了2 000 kg炸药在5 m高度爆炸后冲击波传输过程及其对城市建筑物的毁伤效应。在后处理中，选取同一时刻，利用超压峰值探查计算处理结果，如图12所示。可以看出，在相同的计算条件下，4个城市模型的计算结果差别非常小，说明在不同参数下的可计算处理基本不影响计算结果。经过检验，本文构建的城市模型不存在CAD模型缺陷，模型可计算性满足数值模拟的要求。

Figure 12 Comparison of computation results on different simplified models图12 在不同简化模型上的计算结果比较

6.2 大规模场景建模

以构建纽约城市模型为例，本文选取了710 km2的场景区域。该场景包含陆地、海面和山地，建筑数量约30万栋。构建的可计算模型含217万面片，建模时间约25 min。该实验表明，经过快速可计算处理，基于原始矢量数据构建的城市模型在准确性、效率、模型规模等方面均能满足大规模城市场景数值模拟的需求。

本文构建的模型不仅可以用于数值模拟计算，还可用于虚拟展示。图13显示了采用OpenSceneGraph图形引擎对可计算处理之后的城市模型进行构建，并且在数字地球模块OSGEarth中的渲染效果。

Figure 13 Displaying city model with OpenSceneGraph图13 用OpenSceneGraph展示城市模型

6.3 轮廓简化对模型的影响

为了统计轮廓简化算法在建模中的作用，本文选择了一组简化边长度参数Ts=1 m,2 m,4 m,6 m，其效果如图14所示。

Figure 14 Effect of building contour simplification图14 建筑轮廓简化效果

本文对纽约城市模型开展轮廓简化实验，实验结果如表2所示。

Table 2 Statistics on the effect of contour simplification

随着简化边长度参数Ts的增大，建筑顶点数最多降低33.5%。实验表明，轮廓简化消除了细小几何特征，避免了绝大部分几何潜在错误，使城市模型更适合于数值模拟。

6.4 轮廓收缩对模型的影响

原始数据中存在大量建筑相交或相接，而简单地删除这些建筑会造成模型失真。因此,本文对相交和相接的建筑采用轮廓收缩，在消除模型非正则缺陷的同时尽量减少删除的建筑数量，以保持处理前后模型的一致性。

为了显示轮廓简化算法在建模中的作用，本文选择平行线间距L=0.2，从6个城市中选择了不同尺寸的区域，统计轮廓收缩前后的建筑相交次数，结果如表3所示。可以看到，在OpenStreetMap原始数据中，相交建筑的数量占建筑总数的比例很大(最多占71%)，而采用轮廓收缩后相交次数大幅降低，最多下降95%。

Table 3 Effect of contour shrinkage

6.5 模型规模对建模时间的影响

为了统计模型规模对建模时间的影响，本文从6个城市中选择了不同尺寸的区域进行建模，对模型规模、建筑数量和建模时间进行统计，如表4所示。由于建模分为数据读取、数据解析、可计算处理和模型构建4个阶段，建模时间是这4个阶段的总时间。可以看到，随着模型面片数和建筑数量的增加，各个阶段的执行时间都基本满足线性增长，其中模型构建时间增长最快，而数据读取、解析、可计算处理占用的时间相对较短。因此，本文的建模效率主要受限于图形引擎的性能，不会随着模型规模的增大而急剧增加。图15显示了随着面片数和建筑数量的增加，建模时间增长的速度。

Figure 15 Relation between modeling time and model scale图15 建模时间与模型规模关系

7 结束语

本文提出了大规模城市可计算模型的快速构建方法。通过分析LoD1城市模型中的CAD模型缺陷，本文确定了需要满足的可计算条件。在可计算处理中，本文方法消除了建筑与建筑、建筑与地形之间的相交、相接和裂缝，避免了城市模型的CAD模型缺陷。本文还针对大规模场景提出了图形快速检索算法，避免了计算量的快速增长。实验表明，构建710 km2、30余万栋建筑、217万面片的可计算模型约需要25 min，模型满足可计算条件。本文方法计算效率高、可靠性强，满足了大规模城市数值模拟的要求。

Table 4 Statistics on modeling time

未来的城市可计算建模可在以下方向开展深入研究：(1)增加城市单体的类型，如对道路、水域、植被等建模，使城市内容更加丰富；(2)提升模型复杂度，构建LoD2或LoD3的城市可计算模型；(3)改进模型的可计算性，如进一步消除模型中的潜在几何错误，降低数值模拟计算的复杂度。