基于Snakes算法的道路要素移位軟件模塊

劉遠(yuǎn)剛，郭慶勝，孫雅庚，鄭春燕

(1.武漢大學(xué) 資源與環(huán)境科學(xué)學(xué)院，湖北 武漢 430079； 2. 長(zhǎng)江大學(xué) 地球科學(xué)學(xué)院，湖北 荊州 434023；3. 武漢大學(xué) 測(cè)繪遙感信息工程國(guó)家重點(diǎn)實(shí)驗(yàn)室，湖北 武漢 430079； 4. 嘉應(yīng)學(xué)院 地理科學(xué)與旅游學(xué)院，廣東 梅州 514015)

一、引 言

地圖綜合中的移位是一種為了保持要素之間的空間關(guān)系，保證要素的清晰易辨性或適應(yīng)地圖上的其他要素，而進(jìn)行的解決地圖要素間的沖突的地圖綜合操作。國(guó)內(nèi)外關(guān)于地圖要素移位的研究已經(jīng)有了較長(zhǎng)的歷史，并提出了眾多有效的算法[1-9]。其中 Burghardt和Meier(1997年)提出的基于Snakes的地圖要素移位算法是一種非常適合線(xiàn)狀要素移位的全局最優(yōu)化算法，該算法把計(jì)算機(jī)視覺(jué)中已得到廣泛應(yīng)用的Snakes模型引入到地圖綜合領(lǐng)域，可有效地控制移位引起的后繼沖突問(wèn)題；Bader(2001年)、Galanda(2003年)和吳小芳(2005年)等多位學(xué)者對(duì)該算法進(jìn)行了改進(jìn)，并將其應(yīng)用于線(xiàn)狀要素移位和變形、多邊形要素的移位、放大和夸大等操作[10-13]。本文采用C#編程語(yǔ)言和ArcGIS Engine二次開(kāi)發(fā)組件，設(shè)計(jì)并實(shí)現(xiàn)了一種基于Snakes算法的道路要素移位軟件模塊。該模塊針對(duì)道路要素移位的應(yīng)用需求，采用面向?qū)ο蟮脑O(shè)計(jì)思想，實(shí)現(xiàn)了Snakes移位算法的數(shù)據(jù)模型、邏輯結(jié)構(gòu)和算法流程，從而為Snakes移位算法的研究和應(yīng)用提供軟件支持。

二、Snakes移位算法基本原理

在基于Snakes的地圖要素移位算法中，能量的描述均由最基本的移位量表示，將曲線(xiàn)要素的幾何特征變化產(chǎn)生的能量作為內(nèi)部能量，將空間沖突產(chǎn)生的排斥力作為外部能量，通過(guò)對(duì)內(nèi)外能的計(jì)算，獲得能量最小時(shí)的道路移位后的最優(yōu)形狀和位置。公式(1)是移位量的計(jì)算公式，其中l(wèi)表示線(xiàn)的長(zhǎng)度，(x0，y0)表示原始線(xiàn)的坐標(biāo)，(x,y)表示線(xiàn)移位后的坐標(biāo)，d(s)是移位量的參數(shù)表達(dá)；公式(2)是總能量計(jì)算公式，其中E(d)是總能量，Eint是內(nèi)部能量，Eext是外部能量；內(nèi)部能量由公式(3)求出，其中d′(s)和d″(s)分別是移位量d(s)關(guān)于s的一階導(dǎo)數(shù)和二階導(dǎo)數(shù)，反映了由于移位而產(chǎn)生的形狀變化的大小，α(s)和β(s)參數(shù)決定Snakes模型的彈性和剛性，反映模型的屬性，對(duì)移位效果有一定的控制作用；外部能量由產(chǎn)生沖突時(shí)地圖符號(hào)的重疊產(chǎn)生，其值與重疊的大小成正比，由它促使線(xiàn)移動(dòng)變形，從而解決沖突。

s→d(s)=(x(s)-x0(s),y(s)-y0(s))T,0≤s≤l

(1)

(2)

Snakes模型的原則是保持內(nèi)能和外能之和的總能量值最小，因此需求解公式(2)中的E(d)取最小值時(shí)曲線(xiàn)各處的移位量。利用歐拉方程、有限元方法，并經(jīng)過(guò)一系列變換，求得矩陣方程公式

Kd=f

(4)

式中，K是剛度矩陣；d是由曲線(xiàn)上各點(diǎn)處移位量及其一階導(dǎo)數(shù)組成的移位向量，是方程的未知數(shù)；f是由曲線(xiàn)上各點(diǎn)處所受外力計(jì)算得到的向量。根據(jù)有限元方法，曲線(xiàn)上每一線(xiàn)段上對(duì)應(yīng)的局部矩陣KL、dL、fL的計(jì)算方法如下(以x軸方向?yàn)槔?，y軸方向類(lèi)似)

式中，α、β為α(s)、β(s)的常量形式；h為線(xiàn)段的長(zhǎng)度；x0、x1為線(xiàn)段起點(diǎn)和終點(diǎn)處的坐標(biāo)；d(x0)和d(x1)表示兩點(diǎn)處x方向移位量；d′(x0)和d′(x1)表示兩點(diǎn)處x方向移位量的一階導(dǎo)數(shù)；f(x0)和f(x1)表示兩點(diǎn)處所受外力在x方向上的分量。

對(duì)于整個(gè)道路網(wǎng)，需要對(duì)每一條線(xiàn)段分別計(jì)算局部矩陣，并依次將局部矩陣聚集到全局矩陣K、d、f中，最后分別對(duì)x和y方向上的全局矩陣方程求解，得到各點(diǎn)的移位向量。需要注意的是，K是奇異矩陣，無(wú)法求得其逆矩陣，導(dǎo)致無(wú)法直接對(duì)方程進(jìn)行求解。因此，需要在矩陣中增加邊界條件，將奇異矩陣K轉(zhuǎn)換為常規(guī)矩陣，使之可解。另外，如果沖突區(qū)域目標(biāo)不是很密集，沖突可能在一次操作中就得到解決；但實(shí)際上大多數(shù)的沖突區(qū)域情況復(fù)雜，參與沖突的目標(biāo)較多，一次操作并不能完全解決所有的沖突，此時(shí)就要多次操作，逐步解決沖突。采用迭代最優(yōu)化過(guò)程逐步解決所有沖突的計(jì)算方法如下

(E+γK)d(t)=d(t-1)+γf(t-1)

(8)

式中，E是單位矩陣；t為迭代次數(shù)；d(t)與d(t-1)分別是第t次的移位向量和第t-1次的移位向量；f(t-1)是第t-1次的各點(diǎn)處的受力向量；γ為迭代步長(zhǎng)。

三、道路要素移位軟件模塊的設(shè)計(jì)與實(shí)現(xiàn)

該軟件模塊在Visual Studio 2010開(kāi)發(fā)環(huán)境下，采用C#語(yǔ)言和ArcGIS Engine 10.0二次開(kāi)發(fā)組件進(jìn)行集成開(kāi)發(fā)。地圖顯示功能和數(shù)據(jù)存取功能借助ArcGIS Engine提供的相關(guān)組件實(shí)現(xiàn)，數(shù)據(jù)結(jié)構(gòu)和算法邏輯自主開(kāi)發(fā)完成，其總體上分為主程序和Snakes移位算法模塊兩部分(如圖1所示)。主程序是基于ArcGIS Engine組件的窗口應(yīng)用程序，實(shí)現(xiàn)了地圖顯示和瀏覽、圖層管理、算法全局參數(shù)交互設(shè)置和地圖要素集存取等功能；Snakes移位算法模塊用于實(shí)現(xiàn)整個(gè)算法的業(yè)務(wù)邏輯，主要包括道路網(wǎng)數(shù)據(jù)模型、算法全局參數(shù)列表和算法邏輯實(shí)現(xiàn)幾個(gè)部分。道路網(wǎng)數(shù)據(jù)模型定義了道路網(wǎng)中所包含的道路對(duì)象、道路頂點(diǎn)對(duì)象、道路端點(diǎn)對(duì)象，以及三者之間的關(guān)聯(lián)關(guān)系，且每種對(duì)象都以類(lèi)進(jìn)行封裝，并設(shè)計(jì)了各自的屬性和方法；算法全局參數(shù)列表包括整個(gè)算法的各項(xiàng)參數(shù)，如形狀參數(shù)、迭代步長(zhǎng)、迭代次數(shù)、道路符號(hào)設(shè)置、目標(biāo)比例尺等；算法邏輯實(shí)現(xiàn)是對(duì)Snakes移位算法的實(shí)現(xiàn)，包括計(jì)算剛度矩陣、計(jì)算外力向量，以及矩陣加、減、乘、除、求逆等基本運(yùn)算的實(shí)現(xiàn)。主程序與Snakes移位算法模塊之間通過(guò)可視化界面進(jìn)行數(shù)據(jù)交換。主程序調(diào)用算法模塊時(shí)，系統(tǒng)先向算法模塊中的全局參數(shù)列表寫(xiě)入用戶(hù)設(shè)置好的全局參數(shù)，并將ArcGIS Geodatabase中的原始要素集導(dǎo)入算法模塊，建立對(duì)應(yīng)的道路網(wǎng)數(shù)據(jù)集；然后將全局參數(shù)作為輸入條件運(yùn)行算法邏輯實(shí)現(xiàn)部分對(duì)道路網(wǎng)數(shù)據(jù)集進(jìn)行移位處理，處理后的結(jié)果再導(dǎo)出為新的ArcGIS Geodatabase要素集(結(jié)果數(shù)據(jù)集)顯示到地圖窗口中。

在基于Snakes的道路要素的移位算法中，必須先計(jì)算出整個(gè)道路網(wǎng)的剛度矩陣和受力向量，然后解矩陣方程得到移位向量。對(duì)于前者，不僅需要用到道路網(wǎng)中所有道路頂點(diǎn)的坐標(biāo)數(shù)據(jù)，而且需要知道道路之間的拓?fù)潢P(guān)聯(lián)信息，否則無(wú)法完成不同道路的剛度矩陣的聚合。根據(jù)上述剛度矩陣和受力向量的數(shù)值計(jì)算方法的需要，筆者設(shè)計(jì)了一種簡(jiǎn)單的道路網(wǎng)數(shù)據(jù)模型。圖2中的道路網(wǎng)模型主要定義了RoadNetWork、Road、PointCoord和ConnNode幾個(gè)類(lèi)，分別用于描述與道路網(wǎng)、道路、道路網(wǎng)中的頂點(diǎn)和端點(diǎn)幾類(lèi)對(duì)象的屬性和方法。這幾個(gè)類(lèi)之間相互關(guān)聯(lián)形成一個(gè)包含了基本拓?fù)湫畔⒑妥鴺?biāo)信息的道路網(wǎng)模型。在該模型中，一個(gè)道路網(wǎng)由道路列表、頂點(diǎn)列表和端點(diǎn)列表構(gòu)成；一條道路包含一個(gè)起點(diǎn)和一個(gè)終點(diǎn)(統(tǒng)稱(chēng)為端點(diǎn))，包含多個(gè)(兩個(gè)以上)頂點(diǎn)；每一個(gè)端點(diǎn)同時(shí)也是一個(gè)頂點(diǎn)，它可以作為道路的起點(diǎn)或終點(diǎn)，與多條(一條以上)道路相關(guān)聯(lián)。另外，RoadGrade類(lèi)用于描述道路的等級(jí)屬性特征，在算法中它是設(shè)置單條道路的形狀參數(shù)的基本依據(jù)。圖2中的AlgSnakes類(lèi)是Snakes移位算法邏輯的實(shí)現(xiàn)，該類(lèi)中定義了Snakes算法的全局參數(shù)列表與所有的計(jì)算方法和計(jì)算流程，具體見(jiàn)表1和表2。其中Run_Alg_Snakes方法通過(guò)對(duì)其他計(jì)算方法的調(diào)用實(shí)現(xiàn)了Snakes移位算法的基本流程。

圖1 總體結(jié)構(gòu)圖

圖2 UML類(lèi)圖

表1 AlgSnakes屬性說(shuō)明

表2 AlgSnakes方法說(shuō)明

圖3是Snakes算法用于道路要素移位的基本流程，具體步驟包括：①初始化，設(shè)置各項(xiàng)全局參數(shù)的初始值和創(chuàng)建道路網(wǎng)數(shù)據(jù)模型。全局參數(shù)包括形狀參數(shù)、最大迭代次數(shù)、迭代步長(zhǎng)、目標(biāo)比例尺、道路符號(hào)寬度、地圖目標(biāo)間最小間隔等；道路網(wǎng)數(shù)據(jù)模型通過(guò)導(dǎo)入的道路要素集數(shù)據(jù)得到。②計(jì)算道路網(wǎng)剛度矩陣，先根據(jù)式(5)計(jì)算出組成道路要素的所有線(xiàn)段的局部剛度矩陣，然后根據(jù)道路頂點(diǎn)下標(biāo)將其依次累加到全局剛度矩陣中對(duì)應(yīng)的元素上，最終得到道路網(wǎng)的全局剛度矩陣。③計(jì)算各點(diǎn)的受力，如果所有頂點(diǎn)的受力均為0，則直接結(jié)束，否則執(zhí)行下一步。④計(jì)算道路網(wǎng)受力向量，先根據(jù)式(7) 計(jì)算出組成道路要素的所有線(xiàn)段的局部受力向量，然后根據(jù)道路頂點(diǎn)下標(biāo)將它們依次累加到全局受力向量中對(duì)應(yīng)的元素上，最終得到道路網(wǎng)的全局受力向量。⑤跟據(jù)式(8)解矩陣方程，得到新的移位向量。⑥用新的移位向量更新道路網(wǎng)中各坐標(biāo)點(diǎn)的值，判斷是否達(dá)到迭代最大次數(shù)，若達(dá)到則結(jié)束，否則迭代次數(shù)增加1次，然后轉(zhuǎn)到步驟③繼續(xù)執(zhí)行。

圖4和圖5為本文所實(shí)現(xiàn)的軟件系統(tǒng)運(yùn)行截圖。圖中所用試驗(yàn)數(shù)據(jù)是通過(guò)ArcGIS矢量化得到的來(lái)源于Google Map上某山區(qū)的部分道路網(wǎng)。在進(jìn)行移位操作之前，已經(jīng)完成線(xiàn)目標(biāo)的化簡(jiǎn)、彎曲合并等處理。圖4是對(duì)算法全局參數(shù)進(jìn)行設(shè)置的對(duì)話(huà)框，本例中目標(biāo)比例尺為1∶50萬(wàn)，形狀參數(shù)α和β分別設(shè)置為10 000 000和1 000 000，圖上最小間隔為0.2 mm，最大迭代次數(shù)為2，迭代步長(zhǎng)為0.1，各類(lèi)要素符號(hào)寬度設(shè)置分別為高速公路2.0 mm、國(guó)道1.8 mm、省道1.2 mm、河流中心線(xiàn)0.7 mm。圖5是采用該模塊對(duì)試驗(yàn)數(shù)據(jù)實(shí)施移位操作得到的最終效果。為了便于對(duì)比，將移位前后的道路網(wǎng)同時(shí)顯示于地圖窗口中，從中可以看出，原本處于沖突區(qū)域的道路路段已移出沖突范圍，且道路的整體形狀保持較好，道路間的拓?fù)潢P(guān)系也未被破壞，整體效果良好。

圖3 算法流程圖

圖4 Snakes移位算法參數(shù)設(shè)置

四、結(jié)束語(yǔ)

基于Snakes模型的地圖要素移位算法是一種適合線(xiàn)狀要移位的全局最優(yōu)化算法，該算法具有成熟的理論基礎(chǔ)。本文采用C#編程語(yǔ)言和ArcGIS Engine二次開(kāi)發(fā)組件，設(shè)計(jì)并實(shí)現(xiàn)了一種基于Snakes算法的道路要素移位軟件模塊。使用該模塊對(duì)道路網(wǎng)進(jìn)行移位，在算法中的形狀參數(shù)設(shè)置合適時(shí)，可同時(shí)解決道路網(wǎng)中多條道路之間的沖突問(wèn)題，且能較好地保持道路的形狀特征。通過(guò)本軟件模塊的開(kāi)發(fā)和試驗(yàn)數(shù)據(jù)的測(cè)試，進(jìn)一步驗(yàn)證了基于Snakes的移位算法全局性最優(yōu)化的優(yōu)點(diǎn)，同時(shí)也發(fā)現(xiàn)了該算法在參數(shù)設(shè)置方面自動(dòng)化程度低的缺點(diǎn)，為算法的改進(jìn)指明了方向；而且模塊采用了面向?qū)ο蟮某绦蛟O(shè)計(jì)方法和獨(dú)立于ArcGIS平臺(tái)的數(shù)據(jù)模型，具有較好的可移植性和可擴(kuò)展性，有利于今后進(jìn)一步的改進(jìn)和應(yīng)用。

參考文獻(xiàn)：

[1] LICHTNER W.Computer-assisted Processes of Cartographic Generalization in Topographic Maps[J]. Geo-Processing, 1979,1(1):183-199.

[3] RUAS A. A Method for Building Displacement in Automated Map Generalisation[J]. International Journal of Geographic Information Science, 1998,12(7):789-803.

[4] H?JHOLT P. Solving Local and Global Space Conflicts in Map Generalization: Using a Finite Element Method[J].Cartography and Geographic Information Science, 2000, 27(1):65-73.

[5] WATE J M, JONES C B. Conflict Reduction in Map Generalization Using Iterative Improvement[J]. GeoInformatica,1998, 2 (4):383-407.

[6] HARRIE L E. the Constraint Method for Solving Spatial Conflicts in Cartographic Generalization[J].Cartography and GIS, 1999, 26(1):55-69.

[7] BADER M. Energy Minimization Methods for Feature Displacement in Map Generalization[D]. Zurich: University of Zurich, 2001.

[8] WILSON I D, WARE J M, WARE J A.A Genetic Algorithm Approach to Cartographicmap Generalization[J]. Computers in Industry,2003,52 (3):291-304.

[9] 艾廷華.基于場(chǎng)論分析的建筑群的移位[J].測(cè)繪學(xué)報(bào),2004,35(1):89-94.

[10] BADER M. Energy Minimization Methods for Feature Displacement in Map Generalization[D]. Zurich: University of Zurich, 2001.

[11] GALANDA M.Automated Polygon Generalization in a Multi Agent System[D]. Zurich:University of Zurich,2003.

[12] 吳小芳,杜清運(yùn), 胡月明，等. 基于改進(jìn)Snake模型的道路網(wǎng)空間沖突處理[J].測(cè)繪學(xué)報(bào),2008,37(2):224-229.