王 淳,于長華,唐 昊,李海峰
(中國工程物理研究院計算機應用研究所,四川 綿陽 621999)
數字城市是指以計算機技術、多媒體技術和大規模存儲技術為基礎,以寬帶網絡為紐帶,運用遙感、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布爾運算,使得本文方法在建模效率、模型規模和可計算性等方面均能滿足高性能計算的要求。
本文采用OpenStreetMap[11]矢量數據作為構建數字城市模型的數據源。OpenStreetMap是一個開源地理信息項目,通過標簽和屬性值提供結構化的地圖描述,在路徑規劃、地圖導航、智能交通等領域都得到廣泛應用。目前,已有若干基于OpenStreetMap構建或瀏覽三維模型的軟件,如blender-osm[12]、OSM2World[13]和OSMBuildings[14]等。
基于OpenStreetMap建立三維模型的一般流程參見文獻[15]。本文的目標是可計算建模,需要處理原始數據中的非閉合、相接、相交、懸空等CAD模型缺陷,因此在一般流程的基礎上還要加入“可計算處理”步驟。因此,本文采用的基于OpenStreetMap矢量數據構建城市模型的流程包含4個步驟(如圖1所示):數據讀取、數據解析、可計算處理和模型輸出。每個步驟的具體內容如下所示:
(1)數據讀取。讀取原始數據,即從OpenStreetMap數據中讀取基本元素,從高程數據中讀取頂點坐標。
(2)數據解析。將原始數據轉化為可用于建模的城市單元體數據,包括解析建筑數據(含輪廓、高度、組合方式及其它屬性)和地形三角網。
(3)可計算處理。消除城市單元體中的CAD模型缺陷,使其能夠用于生成可計算城市模型。即消除建筑與建筑、建筑與地形之間的重疊、縫隙、曲率突變等幾何錯誤和細小特征,滿足數值模擬的要求。
(4)模型輸出。將各類城市單元體合并,形成城市模型,并通過圖形引擎輸出。
在上述各個步驟中,數據讀取、數據解析和模型輸出與一般的建模軟件相同。因此,本文重點研究的內容是可計算處理。
本文構建的是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節)。
CAD模型的表達采用的是拓撲和幾何分離的思想。幾何數據描述空間實體的位置,拓撲結構數據描述空間實體間的連接關系。在這種分離的設計方式中,其幾何數據與拓撲結構數據相對獨立實現,然后通過關聯操作建立它們之間的聯系。
根據CAD模型拓撲和幾何的特點,可以把CAD模型缺陷分為拓撲表達缺陷和幾何表達缺陷,如圖3所示。其中拓撲表達缺陷可以進一步分為非正則拓撲缺陷、拓撲異構缺陷和拓撲精度缺陷,幾何表達缺陷可進一步分為幾何異構缺陷、幾何精度缺陷和幾何潛在錯誤表達。詳細的CAD模型缺陷分類及處理方法參見文獻[8]。CAD模型缺陷嚴重破壞了幾何模型的完整性和精確性,會導致數值模擬計算失敗。

Figure 3 Categories of CAD model deficiencies圖3 CAD模型缺陷分類
由于LoD1城市模型的特殊結構,它所對應的CAD模型只會涉及CAD模型缺陷的部分類型。對各類CAD模型缺陷在LoD1城市模型中的表現分析如下:
(1)非正則拓撲缺陷。CAD系統中模型可以是描述物理模型的正則形體,也可以是在現實世界不存在的非正則形體。當 CAD 系統是非正則形體時,一部分拓撲缺陷也被包含在模型中。典型的非正則拓撲缺陷有懸面、懸邊和孤立的點。在LoD1城市模型中,如果建筑地面輪廓不是簡單閉合的多邊形,由其拉伸而成的建筑模型會產生非正則拓撲缺陷。另一種情況是建筑輪廓相交,此時頂點或邊的拓撲連接關系會發生錯誤。
(2)拓撲異構缺陷。拓撲異構是指CAD模型中的模型拓撲關系在計算機中的描述方法不同,因此在文件轉換時,會出現拓撲不相容等問題。本文的城市模型是基于圖形引擎直接生成的,不涉及多個CAD系統間的轉換,因此避免了拓撲異構缺陷。
(3)拓撲精度缺陷。拓撲精度缺陷是由系統精度問題造成的拓撲表達不一致,如曲面之間的裂縫或重疊、邊與邊的裂縫以及邊與頂點的裂縫。對LoD1城市模型來說,當多個建筑模型相接,或當建筑模型與地形模型之間存在重疊或裂縫時,會產生拓撲精度缺陷。
(4)幾何異構缺陷。幾何異構是指異構CAD系統中描述模型實體的數學定義不同,主要發生在復雜曲線曲面(如NURBS曲線曲面)上。現有的建筑模型和地形模型只包含平面,因此避免了二次曲面帶來的幾何異構缺陷。
(5)幾何精度缺陷。幾何精度缺陷是由系統精度問題造成的幾何表達不一致,主要是針對曲面和曲線而言。在LoD1城市模型中不存在這種情況。
(6)幾何潛在錯誤。幾何潛在錯誤主要由角度、尺寸、鄰近和曲率4類情況產生的。大部分的潛在錯誤都是拓撲正確,并能產生合法的實體模型,因此本文不對幾何潛在錯誤進行特別的處理。但是,為了減少數值模擬的計算量,需要消除模型中不影響計算結果的細小幾何特征,這也會消除絕大部分幾何潛在錯誤。
根據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所示的可計算處理流程中各個步驟的具體算法。


Figure 5 Counter-clockwise polygon with self-intersection圖5 逆時針走向多邊形自相交
假設多邊形頂點數是n,判斷自相交將計算n條直線方程和n次直線相交,然后計算n個交點是否在多邊形內部。由于判斷點在多邊形內部的算法復雜度是O(n),該算法復雜度為O(n2)。
建筑輪廓簡化是指對給定的建筑輪廓,在減少建筑輪廓頂點數量的情況下,使簡化后的建筑輪廓與原建筑輪廓盡量保持形狀一致。因為建筑輪廓只有一種拓撲類型,所以本文從目標邊2條鄰邊的位置關系(平行、非平行)來考慮建筑輪廓簡化問題。2種位置關系的簡化如圖6a和圖6b所示,設簡化容差為Ts,當Sn的長度小于Ts時,將Sn設置為待處理邊。圖中加粗的邊為新增加的邊,虛線邊為待刪除的邊,文獻[17]詳細介紹了2種位置關系的簡化方法。

Figure 6 Contour simplification method圖6 輪廓簡化方法
基于以上方法可實現建筑輪廓簡化,在保證建筑模型在外觀特征不失真的情況下,減少三維建筑物模型的細節特征和面片數,以適應不同復雜度要求的數字城市可計算模型。
如圖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)。
在CAD系統中,實體間的拓撲相鄰關系經常出現誤差精度導致的問題。CAD系統判斷一個點是否位于一條直線上并不是由計算的絕對結果來決定,而是通過判斷該點到直線的距離是否小于系統偏差。如果在系統允許的偏差范圍內,則可以表示該點在曲線上,反之說明該點在曲線外。大多數商業CAD系統中使用的最小偏差大約為10-4甚至10-3。

Figure 8 Joint polygons caused by geometric precision圖8 幾何精度造成的相接
為了檢測相接的建筑輪廓,需要計算一個多邊形的所有頂點到另一個多邊形所有邊的最短距離。如果頂點不在邊上但最短距離小于偏差(如圖8所示),為避免拓撲精度缺陷需判定為相接。相接檢測的算法復雜度與相交檢測相同,都是O(n2)。
為避免城市建筑與相鄰建筑碰撞,可將建筑輪廓收縮指定的距離(即安全緩沖距離),即從建筑輪廓頂點沿角平分線方向向建筑輪廓內部平移一定的距離。
如圖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 內縮點計算
實現建筑與地形貼合的關鍵在于計算建筑所在的地表海拔高度,即從建筑頂點發出垂直射線向地形三角網投影,并計算投影點所在的三角面片。
通過射線與三角網求交計算建筑所在地表的海拔高度的方法如下所示。將射線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,則交點位于三角形內。通過以上算法完成所有三角網面片求交測試,得到射線與三角網的交點。
如果建筑通過了高度檢查,即建筑輪廓所有頂點的海拔高度一致且頂點在合法坐標范圍內,則將投影點所在的三角面片的高度作為建筑海拔高度。
對于大規模城市模型來說,可計算處理的計算復雜度成為影響算法性能的關鍵問題。假設建筑數量為M,建筑頂點數量為N,地形三角面片數量為D,則對所有建筑計算相接或相交關系的復雜度為O(M2N2),為所有建筑計算海拔高度的復雜度為O(MND)。大規模城市中通常建筑和地形三角面片的數量大約在105~106量級,如此規模的計算量實際上是無法執行的。為了解決建模效率的問題,本節提出大規模圖形快速檢索算法,以提高圖形檢索和處理的效率,降低可計算處理的復雜度。
在對建筑模型的可計算處理中,建筑地面輪廓的相交和相接檢測占用了大量時間。為此本文提出基于最大包圍盒的鄰近建筑輪廓快速檢測算法。通過快速提取與每棟建筑鄰近的建筑,減少相交檢測次數。
如圖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 鄰近建筑輪廓快速檢測
由于城市模型中包含大量的建筑頂點和地形三角面片,而普通的求交涉及大量冗余的射線和三角面片,計算量過高,為此本文提出基于層次化包圍盒的快速求交算法。
本文將整個三角網區域進行層次化分解,待計算出最小的包圍盒后,再對其內部的三角網面片進行求交測試。該算法剔除了包圍盒外部的面片,因此大大減少了求交測試的次數。具體算法步驟如下所示:
(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)。
本文在服務器單個結點上構建可計算模型。每個結點包括144核,主頻為2.3 GHz。對可計算處理的關鍵步驟(如建筑輪廓求交、簡化、地形貼合)采用OpenMP并行。
為了驗證模型的可計算性,本文選擇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 在不同簡化模型上的計算結果比較
以構建紐約城市模型為例,本文選取了710 km2的場景區域。該場景包含陸地、海面和山地,建筑數量約30萬棟。構建的可計算模型含217萬面片,建模時間約25 min。該實驗表明,經過快速可計算處理,基于原始矢量數據構建的城市模型在準確性、效率、模型規模等方面均能滿足大規模城市場景數值模擬的需求。
本文構建的模型不僅可以用于數值模擬計算,還可用于虛擬展示。圖13顯示了采用OpenSceneGraph圖形引擎對可計算處理之后的城市模型進行構建,并且在數字地球模塊OSGEarth中的渲染效果。

Figure 13 Displaying city model with OpenSceneGraph圖13 用OpenSceneGraph展示城市模型
為了統計輪廓簡化算法在建模中的作用,本文選擇了一組簡化邊長度參數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%。實驗表明,輪廓簡化消除了細小幾何特征,避免了絕大部分幾何潛在錯誤,使城市模型更適合于數值模擬。
原始數據中存在大量建筑相交或相接,而簡單地刪除這些建筑會造成模型失真。因此,本文對相交和相接的建筑采用輪廓收縮,在消除模型非正則缺陷的同時盡量減少刪除的建筑數量,以保持處理前后模型的一致性。
為了顯示輪廓簡化算法在建模中的作用,本文選擇平行線間距L=0.2,從6個城市中選擇了不同尺寸的區域,統計輪廓收縮前后的建筑相交次數,結果如表3所示。可以看到,在OpenStreetMap原始數據中,相交建筑的數量占建筑總數的比例很大(最多占71%),而采用輪廓收縮后相交次數大幅降低,最多下降95%。

Table 3 Effect of contour shrinkage
為了統計模型規模對建模時間的影響,本文從6個城市中選擇了不同尺寸的區域進行建模,對模型規模、建筑數量和建模時間進行統計,如表4所示。由于建模分為數據讀取、數據解析、可計算處理和模型構建4個階段,建模時間是這4個階段的總時間。可以看到,隨著模型面片數和建筑數量的增加,各個階段的執行時間都基本滿足線性增長,其中模型構建時間增長最快,而數據讀取、解析、可計算處理占用的時間相對較短。因此,本文的建模效率主要受限于圖形引擎的性能,不會隨著模型規模的增大而急劇增加。圖15顯示了隨著面片數和建筑數量的增加,建模時間增長的速度。

Figure 15 Relation between modeling time and model scale圖15 建模時間與模型規模關系
本文提出了大規模城市可計算模型的快速構建方法。通過分析LoD1城市模型中的CAD模型缺陷,本文確定了需要滿足的可計算條件。在可計算處理中,本文方法消除了建筑與建筑、建筑與地形之間的相交、相接和裂縫,避免了城市模型的CAD模型缺陷。本文還針對大規模場景提出了圖形快速檢索算法,避免了計算量的快速增長。實驗表明,構建710 km2、30余萬棟建筑、217萬面片的可計算模型約需要25 min,模型滿足可計算條件。本文方法計算效率高、可靠性強,滿足了大規模城市數值模擬的要求。

Table 4 Statistics on modeling time
未來的城市可計算建模可在以下方向開展深入研究:(1)增加城市單體的類型,如對道路、水域、植被等建模,使城市內容更加豐富;(2)提升模型復雜度,構建LoD2或LoD3的城市可計算模型;(3)改進模型的可計算性,如進一步消除模型中的潛在幾何錯誤,降低數值模擬計算的復雜度。