吳 凱 馬藝文
(1. 濟南市勘察測繪研究院, 山東 濟南 250013; 2. 武漢市國土資源和規(guī)劃信息中心, 湖北 武漢 430014)
在現(xiàn)實世界中,兩個目標物之間可以直接通行的情況較少,一般都需要繞過若干障礙物才可到達。地理空間問題中有很多典型的障礙空間問題,例如:軍事問題中如何選擇合適的進軍路線以避開難以逾越的地理障礙及人為障礙[1];日常生活中交通工具繞開擁堵路段或施工路段的最短路徑規(guī)劃[2];災害救援中“安全島”問題[3];無人駕駛飛機的防撞避障和導航設計問題[4];科學研究中機器人避障路徑規(guī)劃問題[5];迷宮求解問題[6-7];障礙空間最近鄰查詢問題[8]。障礙空間,指空間中具有各種障礙物,此處的障礙物并非單指點或線或多邊形,而是指全形態(tài)圖形,它是點、線、面以及它們的組合,障礙本身不計算距離,并且不傳遞距離。
距離變換是計算并標識空間點(對參照體)距離的變換或過程。障礙距離變換,即在空間中存在障礙物時,距離的傳播需要繞過障礙物,它屬于條件距離變換[9-10]。這類問題比較抽象卻來源于實際地理世界的需要,障礙距離變換是解決此類地理空間問題的有效途徑。本文主要對比了兩種障礙距離變換柵格生成方法,并進行了精度分析和對比。
距離變換是將包含實體特征和空間背景兩種像元的二值圖像轉變?yōu)榫嚯x圖像的變換。在距離圖像中,每一個像素值表示該像素到其最近的一個實體像素的距離,具體體現(xiàn)為每個實體的距離波不斷地往外空間進行擴張,直到與鄰近實體的距離波相遇。距離變換的結果往往是空間等距離線,任一點的趨源梯度線(即該線上任一點的切線方向都垂直于過該點的等距離線),即該點趨源的最短路徑。
當考慮障礙空間問題時,上述距離變換則應擴展為DTO(Distance Transformation with Obstacles),即在障礙空間中進行距離變換。空間中有生成元還有若干障礙,生成元傳播的距離波需要繞過障礙進行傳播。由以上定義可知,一般意義上的距離變換可以看作是障礙物個數(shù)為零的障礙距離變換。
地圖代數(shù)的障礙距離變換(MA-DTO)在二維障礙空間中是通過柵格路徑距離變換實現(xiàn)的。可根據(jù)精度要求確定開窗大小,采用3×3,5×5等奇數(shù)(2n+1)×(2n+1)模板。模板柵格路徑定義為模板的中心至其他各個柵格中心間的有限個方向線段的組合,模板柵格路徑所對應的長度稱為模板柵格路徑長度。
本文主要利用ArcGIS的成本距離(Cost Distance)工具,來計算每個像元到成本面上最近源的最小累計成本距離,即ArcGIS-DTO。在使用Cost Distance工具進行障礙距離變換時,生成元(Source Data)可以為矢量數(shù)據(jù)也可以為柵格數(shù)據(jù),如果是柵格數(shù)據(jù)則具有 NoData 值的像元不包括在源集內,成本柵格數(shù)據(jù)(Cost Raster)中的NoData部分充當成本面中的障礙。
成本距離的算法如下:像元八鄰域從像元0移動到四個與其邊連接(1、3、5、7四像元)的近鄰之一,取像元1時,跨越連接線移動到相鄰結點的成本為a1=(cost0+cost1)/2,其中:cost0代表像元0的成本,cost1代表像元1的成本,a1為從像元0到像元1連接線的總成本。如果沿對角線移動,例如從像元0到像元8,則連接線上的行程成本為a2=1.414 214(cost0+cost8)/2,其中cost0代表像元0的成本,cost8代表像元8的成本,a2為從像元0到像元8連接線的總成本,如圖1所示。
圖1 像元八鄰域
MA-DTO算法是一個統(tǒng)一的、通用的算法,具有優(yōu)良的數(shù)據(jù)初始化、可視化性能,并且障礙物和發(fā)生元為全形態(tài)圖。柵格空間下兩定點之間的距離由其所屬的兩柵格中心間的距離定義。在障礙空間下由于連通概念的改變,障礙空間下的柵格路徑定義為指定兩定點a、b所屬柵格中心間并由相通柵格中心間線段組成的有序線段集。柵格路徑所對應的累計長度稱為柵格路徑長度。一般情況下,兩定點之間的柵格路徑并不唯一,往往有多條,所有路徑中最短累計長度的柵格路徑最有意義。實際上基于地圖代數(shù)的距離變換給出了整個障礙空間所有點的趨源距離,并且給出了多生成元的整個障礙空間所有點的趨最近元距離。
ArcGIS-DTO的實現(xiàn)工具,成本加權距離通過將距離等同為成本因子來修改“歐氏”距離,該成本即為經(jīng)過任何指定像元所需的成本。成本分配工具可以基于累積行程成本來識別最近(或成本最小的)源的像元。成本回溯鏈接工具將提供一個道路地圖,用于標識從任何像元開始,沿著最小成本路徑,返回到最近源的路徑。
下面以一個點的距離變換為例,對比其距離波的形態(tài),如圖2所示,圖2(a)為一個點的歐氏距離變換圖,其產(chǎn)生的距離波為同心圓;圖2(b)為ArcGIS中Cost Distance所使用的3×3模板距離變換,可以看出由于3×3模板的八鄰域結構,生成的距離波是以等邊八邊形的方式逼近歐氏距離變換的圓形;圖2(c)為地圖代數(shù)方法一個點的距離變換結果,使用的是17×17模板,從距離波的形態(tài)上來看,已經(jīng)十分精確地逼近歐氏距離變換結果。
圖2 一個點的距離變換對比
用ArcGIS的距離變換結果疊加歐氏距離變換的同心圓,來觀察ArcGIS實現(xiàn)結果的誤差分布情況。如圖3所示,為ArcGIS的3×3模板距離變換結果。
圖3 誤差分布示意圖
疊加歐氏距離變換距離波同心圓可以看出,在八方向上,如線段OP所在位置,ArcGIS的結果和歐氏距離相等;在兩個相鄰八方向的中點,如線段OR所在位置,誤差最大;誤差隨著向八方向的靠近而減小。設線段OP、OQ、OR所處位置處的誤差分別為ΔOP、ΔOQ、ΔOR,則ΔOP<ΔOQ<ΔOR。
以線狀生成元在障礙空間中的距離變換試驗來進行兩種DTO實現(xiàn)方法的精度對比,如圖4所示,圖4(a)中直線段生成元,螺旋區(qū)域內部為距離變換空間,外部為障礙空間;圖4(b)為ArcGIS實現(xiàn)的DTO圖;圖4(d)為地圖代數(shù)方法實現(xiàn)的DTO圖。圖4(c)和圖4(e)分別為兩種方法實現(xiàn)結果的局部放大圖,放大部分為螺旋彎曲單元的轉彎區(qū)域。對比圖4(c)和4(e)可以看出:地圖代數(shù)方法的距離波更平滑、更細膩,而ArcGIS的距離波在轉彎區(qū)域容易形成折線,距離波紋不夠平滑。
圖4 線狀生成元障礙距離變換精度對比
如圖5所示,圖5(a)中點為生成元,共7個,線和面為障礙物;圖5(b)為ArcGIS實現(xiàn)的多點生成元DTO圖;圖5(c)為地圖代數(shù)方法實現(xiàn)的DTO圖。兩種方法所得距離波的形態(tài)不同之處前面已經(jīng)進行了討論。在此,布設距離變換圖上n個隨機點作為檢查點,通過對比它們所在位置處兩種距離變換結果的距離值,分析距離差值的變化規(guī)律。
圖5 多點生成元DTO
在圖幅范圍內隨機選取40個檢查點,如圖6所示,三角形點在此用P1,P2,…,P40來表示。圖6(a)為檢查點與ArcGIS-DTO結果的疊加圖,圖6(b)為檢查點與MA-DTO結果的疊加圖。
圖6 隨機點與DTO結果的疊加
采集40個檢查點所在位置,ArcGIS-DTO和MA-DTO所得到的兩個距離變換結果圖中的距離值。如表1所示,D1為ArcGIS所得結果的距離值,D2為地圖代數(shù)方法所得結果的距離值。再求出兩種距離值的差值(ΔD),ΔD=D1-D2。對40個檢查點按照距離值從小到大的順序進行排序,排序后的結果如表2所示。
表1 隨機點處距離值及距離差值 單位:m
表2 按距離值從小到大排序 單位:m
由于MA-DTO方法已十分接近于歐氏距離變換,所以本文將MA-DTO方法設為標準,考察ArcGIS-DTO的誤差。則表2中距離差值ΔD可表示為ArcGIS-DTO的誤差值。利用誤差值(ΔD)和ArcGIS-DTO的距離值(D1)繪制成折線圖,即誤差隨距離值變大的變化趨勢圖。如圖7所示,圖7(a)中的0值點即為上述的八方向上的點,在統(tǒng)計變化趨勢前我們先將0值點去除,統(tǒng)計其余37個點的趨勢圖,如圖7(b)所示。
圖7 ArcGIS-DTO結果誤差值隨距離值的變化趨勢圖
對圖7(b)所示的誤差值隨距離值的變化趨勢圖進行線性回歸,如圖8所示,線性回歸函數(shù)為:
圖8 誤差隨距離值變化趨勢線性回歸
y=1.926 54+0.062 38x(x為距離值,y為誤差值),相關系數(shù)R=0.811 48??梢钥闯?雖然距離差值有波動,但整體上來講,隨著距離值的增大,也就是距離生成元越遠,ArcGIS-DTO的誤差越大。
通過可視化和量化兩種方式的精度對比,本文得出以下結論:ArcGIS-DTO和MA-DTO在實現(xiàn)障礙空間的距離變換時都是模擬無障礙空間中歐氏距離變換的思想。地圖代數(shù)算法使用的模板大小可以為3×3、5×5、奇數(shù)(2n+1)×(2n+1),隨著n的增大誤差會越來越小,越接近歐氏距離,精度也越高,但算法的復雜性也會增加;ArcGIS中使用的算法是固定的3×3模板,所以靈活性比較差,精度相對也較低些。