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

劉遠剛，郭慶勝，孫雅庚，鄭春燕

(1.武漢大學 資源與環境科學學院，湖北 武漢 430079； 2. 長江大學 地球科學學院，湖北 荊州 434023；3. 武漢大學 測繪遙感信息工程國家重點實驗室，湖北 武漢 430079； 4. 嘉應學院 地理科學與旅游學院，廣東 梅州 514015)

一、引 言

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

二、Snakes移位算法基本原理

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

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

(1)

(2)

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

Kd=f

(4)

式中，K是剛度矩陣；d是由曲線上各點處移位量及其一階導數組成的移位向量，是方程的未知數；f是由曲線上各點處所受外力計算得到的向量。根據有限元方法，曲線上每一線段上對應的局部矩陣KL、dL、fL的計算方法如下(以x軸方向為例，y軸方向類似)

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

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

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

(8)

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

三、道路要素移位軟件模塊的設計與實現

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

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

圖1 總體結構圖

圖2 UML類圖

表1 AlgSnakes屬性說明

表2 AlgSnakes方法說明

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

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

圖3 算法流程圖

圖4 Snakes移位算法參數設置

四、結束語

基于Snakes模型的地圖要素移位算法是一種適合線狀要移位的全局最優化算法，該算法具有成熟的理論基礎。本文采用C#編程語言和ArcGIS Engine二次開發組件，設計并實現了一種基于Snakes算法的道路要素移位軟件模塊。使用該模塊對道路網進行移位，在算法中的形狀參數設置合適時，可同時解決道路網中多條道路之間的沖突問題，且能較好地保持道路的形狀特征。通過本軟件模塊的開發和試驗數據的測試，進一步驗證了基于Snakes的移位算法全局性最優化的優點，同時也發現了該算法在參數設置方面自動化程度低的缺點，為算法的改進指明了方向；而且模塊采用了面向對象的程序設計方法和獨立于ArcGIS平臺的數據模型，具有較好的可移植性和可擴展性，有利于今后進一步的改進和應用。

參考文獻：

[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] 艾廷華.基于場論分析的建筑群的移位[J].測繪學報,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] 吳小芳,杜清運, 胡月明，等. 基于改進Snake模型的道路網空間沖突處理[J].測繪學報,2008,37(2):224-229.