王加勝,楊 昆,朱彥輝,熊建紅

1. 云南師范大學(xué)信息學(xué)院，云南 昆明 650500； 2. 西部資源環(huán)境地理信息技術(shù)教育部工程研究中心，云南 昆明 650500

距離變換最早由Rosenfeld和Pfaltz提出，是指將每個柵格像元映射為到最近感興趣區(qū)域或目標(biāo)區(qū)域的距離，其輸入圖像一般為二值圖像[1]。其中的距離一般為歐氏距離。由于根據(jù)定義計算復(fù)雜度較高，文獻[2]提出了一種距離變換的快速實現(xiàn)算法——楔形距離變換。該方法僅考慮鄰域像元的影響，通過兩次掃描完成距離變換。該算法速度快，但無法得到精確的結(jié)果。此后，許多學(xué)者通過改變鄰域大小和優(yōu)化距離算子提高楔形距離變換精度[3-6]。

2004年，文獻[7]將距離變換引入到GIS中，討論了距離變換在GIS領(lǐng)域的應(yīng)用場景。目前距離變換已成為GIS中的一個基本工具[8-12]。在GIS鄰域距離變換的許多應(yīng)用場景中，常常需要考慮障礙區(qū)域的影響，計算非障礙柵格像元到最近目標(biāo)最短路徑長度(如航線設(shè)計、海上救助等)[13-16]。這種距離變換，本文稱之為繞障歐氏距離變換。目前，未見有專門針對繞障距離變換的研究。

元胞自動機(cellular automata，CA)是一種時間和空間都離散的動力學(xué)系統(tǒng)[17]，廣泛應(yīng)用于城市規(guī)劃[18-19]、土地利用[10]、傳染病傳播[20]、交通系統(tǒng)[21]、人群疏散等領(lǐng)域[22-24]。例如，文獻[25]構(gòu)建元胞自動機模型，分析了個體行為對人群疏散效果的影響；文獻[26]將元胞自動機與系統(tǒng)生物學(xué)結(jié)合研究癌癥病毒的演化過程；文獻[27]集成CA與隨機森林方法模擬城市擴張過程；文獻[28]構(gòu)建了城市和區(qū)域規(guī)劃元胞自動機模型。

元胞自動機的元胞與柵格數(shù)據(jù)的像元有相似之處，距離變換可以離散成為由目標(biāo)像元向外擴散和更新繞障距離值的過程，因此可用元胞狀態(tài)變化來反映距離變換過程。文獻[29]提出了一種基于元胞自動機的柵格空間復(fù)雜實體加權(quán)距離變換方法，但未考慮障礙的影響。本文擬提出一種基于元胞自動機的繞障距離變換方法，并以南海及鄰近海域為例，分析該方法計算結(jié)果的精度與適用性。

1 研究區(qū)與數(shù)據(jù)

1.1 研究區(qū)

本文選擇南海及其鄰域作為研究區(qū)(0°N—24°N，98°E—122°E)，以測試方法的精度與適用性。南海及鄰域位于亞洲東南部，海域包括南海和泰國灣的全部，蘇祿海和蘇拉威西海的一部分，以及馬六甲海峽、臺灣海峽、巴士海峽等眾多海峽，是太平洋與印度洋的連接地帶。陸地涉及中國南部沿海、南海諸島、中南半島、菲律賓群島和大巽他群島。該區(qū)域島礁眾多，大小不一，陸地海岸線狹長復(fù)雜，這些特性為分析繞障歐氏距離變換結(jié)果提供了良好的測試環(huán)境。

1.2 數(shù)據(jù)及預(yù)處理

本文使用的數(shù)據(jù)資料包括陸地島礁圖層和目標(biāo)點圖層。陸地島礁圖層使用中國及東南亞的行政區(qū)劃Shapefile矢量數(shù)據(jù)，通過國家地理信息公共服務(wù)平臺服務(wù)資源中全球矢量地圖服務(wù)中對南海及鄰域矢量化陸地范圍得到，比例尺為1∶1000萬，坐標(biāo)系統(tǒng)為2000國家大地坐標(biāo)系。目標(biāo)點圖層通過在研究區(qū)海域隨機選取7個點，用于分組測試不同數(shù)量目標(biāo)情況下繞障歐氏距離變換效果。

距離分析需在投影坐標(biāo)系統(tǒng)下進行。在導(dǎo)入模型前，將數(shù)據(jù)轉(zhuǎn)換為Lambert-等角圓錐投影坐標(biāo)系統(tǒng)，標(biāo)準(zhǔn)緯線為6°N與18°N，中央經(jīng)線為110°E，單位為千米。由于距離變換結(jié)果用柵格表示，需要將投影轉(zhuǎn)換后的行政區(qū)劃圖層轉(zhuǎn)換為柵格圖層(像元大小為5 km)，1表示陸地或島礁(即障礙區(qū))，0表示水域(正常區(qū)域)，轉(zhuǎn)換后的柵格大小為500×500。同時將目標(biāo)點位置處的柵格值設(shè)置為2。處理后的柵格在此稱為元胞的環(huán)境柵格。

2 研究方法

2.1 基本概念

在介紹模型之前，有必要對文中涉及的歐氏距離、歐氏距離變換、繞障歐氏距離變換、元胞自動機等基本概念予以明確。假設(shè)二值柵格圖像尺寸為m×n，Ci(xi,yi)、Cj(xj,yj)為圖像中的任意兩像元。O為目標(biāo)像元集合，B為障礙像元集合。

(1) 歐氏距離D(Ci,Cj)。用兩點坐標(biāo)差值平方和的平方根表示，即

(1)

(2) 歐氏距離變換。變換的結(jié)果使得每個像元值為該像元至最近目標(biāo)點的歐氏距離。設(shè)像元p的歐氏距離變換值用DT(p)表示。q為O中的任意一像元。則歐氏距離變換可由式(2)表示

DT(p)=min{D(p,q),q∈O}

(2)

(3) 繞障歐氏距離。表示任一非障礙集像元經(jīng)過非障礙區(qū)到達最近目標(biāo)像元最短路徑的長度。

(4) 繞障歐氏距離變換。變換結(jié)果使得任一非障礙集像元值為該像元到目標(biāo)像元集的繞障歐氏距離。設(shè)Lm(p…q)表示從p經(jīng)過非障礙區(qū)到q路徑的最小長度。像元p的繞障歐氏距離變換可表示為

EDT(p)=min{Lm(p…q),q∈O}

(3)

D(Ck,p)}

(4)

式中，C1,Ci,Ck?B。

(5) 元胞自動機。元胞自動機是一種時間、空間、狀態(tài)都離散，空間相互作用和時間因果關(guān)系都為局部的網(wǎng)格動力學(xué)模型，具有復(fù)雜系統(tǒng)時空演化過程的能力。每一元胞取有限的狀態(tài)、遵循同樣的規(guī)則、通過簡單的相互作用構(gòu)成動態(tài)系統(tǒng)的演化。元胞自動機可用一個四元組描述。A={Ld,S,N,f}，其中Ld標(biāo)識d維元胞空間；S為元胞有限狀態(tài)；N為鄰域元胞集合；f表示中心元胞的狀態(tài)轉(zhuǎn)換規(guī)則。

2.2 繞障歐氏距離變換CA模型

CA模型的關(guān)鍵在于元胞與狀態(tài)轉(zhuǎn)換規(guī)則的定義。在傳統(tǒng)的元胞自動機模型中，元胞的狀態(tài)是有限的整數(shù)。本文將元胞狀態(tài)定義為實數(shù)類型，以模擬距離變換計算過程。繞障歐氏距離變換元胞自動機模型架構(gòu)如圖1所示，包括模型初始化、元胞狀態(tài)轉(zhuǎn)換、模型表達與輸出3個部分。

圖1 距離分析元胞自動機模型架構(gòu)Fig.1 The technique route of distance analysis CA model

2.2.1 元胞定義

將研究區(qū)根據(jù)像元大小(柵格分辨率，用Sc表示)劃分為m×n正方形格網(wǎng)，每個單元即為一個元胞。本文將元胞的狀態(tài)定義為元胞中心位置至目標(biāo)點集的繞障歐氏距離(格數(shù))，用dt表示。此外，元胞還具有位置(即所在列x,所在行y)、狀態(tài)變化標(biāo)志(fs)，障礙區(qū)標(biāo)志(fl)兩個屬性。位置為(x，y)元胞表示為C(x，y)。元胞屬性的具體描述見表1。

表1 元胞屬性說明

2.2.2 模型初始化

模型初始化即模型運行前的元胞屬性初始化與參數(shù)設(shè)置。模型的參數(shù)包括環(huán)境柵格和像元大小Sc。讀入環(huán)境柵格數(shù)據(jù)后賦予元胞的障礙標(biāo)志屬性fl。Sc可以根據(jù)實際要求的精度確定，影響著元胞的大小與數(shù)量；狀態(tài)dt的初始值分為兩種情況，目標(biāo)點所在元胞狀態(tài)dt設(shè)為0，其余元胞為+∞(比最大距離大的數(shù)，如w×h)；狀態(tài)變化標(biāo)志ft設(shè)置為0；x、y則根據(jù)元胞位置確定。元胞的屬性中fl、x、y初始化后并不會改變。

2.2.3 元胞狀態(tài)轉(zhuǎn)換

元胞狀態(tài)轉(zhuǎn)換是整個模型的核心，根據(jù)相鄰元胞的狀態(tài)，計算當(dāng)前元胞的狀態(tài)。距離計算過程可以看作由目標(biāo)像元向周圍逐漸擴散的過程，目標(biāo)與周圍的8個像元的距離可根據(jù)歐氏距離計算得出，再加上周圍元胞的最小狀態(tài)值，即可得到元胞的狀態(tài)。本文每個時刻的元胞狀態(tài)轉(zhuǎn)換就是更新元胞距離值。狀態(tài)轉(zhuǎn)換涉及鄰域定義和轉(zhuǎn)換規(guī)則。

(1) 鄰域。本文采用鄰近的8個元胞表示當(dāng)前元胞的鄰域，即元胞C(x,y)的鄰域元胞有C(x-1,y-1)、C(1,y-1)、C(x+1,y-1)、C(x-1,y)、C(x+1,y)、C(x-1,y+1)、C(x+1,y)、C(x+1,y+1)。

(2) 狀態(tài)轉(zhuǎn)換規(guī)則：

陸域元胞：不進行狀態(tài)轉(zhuǎn)換。

鄰域皆未計算的元胞：不進行狀態(tài)轉(zhuǎn)換。

其他元胞：設(shè)其與鄰域狀態(tài)(至目標(biāo)最短路徑長度)組成的3×3矩陣為Dt，鄰域至中心元胞距離的3×3矩陣為A(稱為距離算子，具有中心對稱性，中心為0)。元胞下一狀態(tài)值(至目標(biāo)最短路徑長度)即為矩陣Dt+A的最小元素值。這一規(guī)則轉(zhuǎn)化過程可用式(5)—式(7)表示

(5)

(6)

dt+1=min(Dt+A)

(7)

距離算子的選擇會影響到變換的精度和速度。在楔形距離變換中，常見的距離算子及其最大絕對錯誤率見表2。其中最小的距離算子為Butt-Maragos優(yōu)化算子。本文使用該算子距離變換。

表2 常用3×3距離算子的最大絕對錯誤率[7]

注：a為算子上下左右的值；b為4個`對角的值。

2.2.4 結(jié)果表達與輸出

模型初始化后，對每一時刻t，元胞根據(jù)狀態(tài)轉(zhuǎn)換規(guī)則進行繞障歐氏距離變換和數(shù)據(jù)更新，直到所有的元胞狀態(tài)都穩(wěn)定為止(對任意元胞，fs=0)，就得到模型運行結(jié)果。將結(jié)果乘以像元大小Sc即可得到該像元大小對應(yīng)的繞障歐氏距離變換結(jié)果。

假設(shè)所有元胞集合為C，則模型停止的條件可表述為：對?c∈C，有dt+1=dt，則fs=0，模型停止。

元胞自動機的優(yōu)勢在于能夠展示元胞演化的過程。本文應(yīng)用可視化工具Processing3.0展示模型演化與距離動態(tài)計算過程。模型中可視化的內(nèi)容包括海陸分布、目標(biāo)位置和元胞狀態(tài)，其中最核心的是元胞的展示。海陸分布作為背景，根據(jù)柵格值將海洋和陸地分別設(shè)置為藍色與淡黃色；目標(biāo)位置采用圓點表示；元胞狀態(tài)的顯示采用漸變顏色表示，涉及顏色的映射，通過運行一遍模型后獲取最大的狀態(tài)值，再進行顏色映射以獲取良好的可視效果。

模型的目的除了觀察計算過程，最重要的是結(jié)果分析。將元胞的最終狀態(tài)值按照行列號組合形成一個二維數(shù)組，根據(jù)像元大小Sc、左上角投影坐標(biāo)，行列數(shù)，輸出為地理信息系統(tǒng)(geographical information system,GIS)軟件支持的ASCII碼格式文件，以便后期分析處理。

3 結(jié)果與討論

3.1 繞障歐氏距離變換模擬結(jié)果

將目標(biāo)點數(shù)據(jù)與預(yù)處理后的環(huán)境柵格輸入模型，創(chuàng)建250 000個元胞組成的CA模型，采用表2最后一個算子，模擬繞障歐氏距離變換過程，結(jié)果如圖2所示?？梢钥闯?，模型動態(tài)地展現(xiàn)出距離變換的計算過程。計算的元胞數(shù)量圍繞目標(biāo)點呈正方形向外逐漸增多，顏色也越來越深，說明元胞狀態(tài)隨著擴散更新。當(dāng)兩個目標(biāo)點的擴散圈相遇時(如圖2中的B位置)，接縫處(T=150的B位置)過渡自然，沒有出現(xiàn)橫豎形狀的突變狀態(tài)。這是由于計算過的元胞依然會根據(jù)鄰域情況作狀態(tài)更新。當(dāng)擴散至障礙區(qū)域時，會繞過陸地，元胞狀態(tài)值(如圖2中的A點)反映了繞障歐氏距離。在T=240時，表示模型經(jīng)過240次迭代，從目標(biāo)經(jīng)海域能到達的區(qū)域元胞完成距離計算，共運行了23.363 s。以上結(jié)果表明，基于CA的海上距離分析模型能夠動態(tài)展示計算過程，距離計算自動考慮到了繞障距離花費，說明模型可行。

圖2 CA模擬過程結(jié)果Fig.2 Results of CA simulation process

3.2 模擬結(jié)果與普通歐氏距離變換結(jié)果對比

為了分析模擬結(jié)果的準(zhǔn)確度，將CA模擬結(jié)果與普通歐氏距離變換結(jié)果對比。為避免其他因素影響，只選擇1個目標(biāo)點，分別對研究區(qū)作CA繞障距離變換(圖3(a))與歐氏距離變換(圖3(b))，再對兩個結(jié)果進行差值計算(圖3(c))。為直觀區(qū)分，圖3用等值區(qū)域法表示?？梢钥闯?，模擬結(jié)果與歐氏距離分析結(jié)果具有相似的空間分布特征，大部分區(qū)域數(shù)值差距在50 km以內(nèi)，CA模擬結(jié)果比普通歐氏距離普遍偏高(圖3(c)中輻條之外的區(qū)域)。但在局部區(qū)域存在較明顯差異，如圖3中的A、B、C3個位置，普通歐氏距離變換由于不考慮障礙導(dǎo)致誤差偏大，而本文方法顧及到了這一點。這些區(qū)域二者差距在300 km以上，差距最大的出現(xiàn)在馬六甲海峽(圖3中的C位置)。

另外，比較結(jié)果中還存在一些低于-50 km的負值，這是由于這些位置在研究區(qū)內(nèi)被障礙隔開，而無法繞障到達目標(biāo)點，CA模型中未參與計算(值為-1)而歐氏距離變換結(jié)果存在值所造成的。綜上，基于CA的距離分析方法在海上距離計算上存在很大的優(yōu)勢，合理考慮了繞陸地距離和無法到達的區(qū)域。

3.3 基于CA的繞障距離變換結(jié)果誤差分析

本文中的模型如果不考慮障礙分布，運算結(jié)果便可與歐氏距離變換結(jié)果對比，分析本文方法的誤差。在不考慮障礙分布情況下，普通歐氏距離變換結(jié)果即為真值。這時，設(shè)CA模擬結(jié)果為Ra，真值為Rb，則誤差E可由式(4)表示

E=(Ra-Rb)/Ra

(8)

據(jù)此計算得到的模擬結(jié)果、真值和誤差分布如圖4所示。從圖4(a)與圖4(b)對比可以看出，本文模擬結(jié)果等值線呈現(xiàn)出正八邊形的特征，這是由于鄰域考慮的是8鄰域所導(dǎo)致的。而真值是呈圓型分布，這成為了誤差的主要來源。從圖4(c)可看到誤差的分布具有傘狀特征。正八邊形邊的中點(即米字形表示的方向)處誤差為0，其余位置誤差則與其與目標(biāo)連線與8個方向最小夾角有關(guān)，夾角越大則誤差也越大。從誤差統(tǒng)計結(jié)果可以看出CA距離分析方法的誤差范圍為0～3.96%。

圖3 模擬結(jié)果與歐氏距離分析結(jié)果對比Fig.3 Comparison results of CA simulation and Euclidean distance transformation

圖4 距離分析CA模型誤差分析Fig.4 Error analysis of CA based distance analysis

在有障礙的情況下，繞障距離變換值為到目標(biāo)最短路徑的長度Lp，可以轉(zhuǎn)化為多條線段的長度和。即

(9)

(10)

因此，本文提出的繞障歐氏距離變換的最大誤差率也為3.96%。

3.4 元胞大小和轉(zhuǎn)換規(guī)則對結(jié)果的影響

3.4.1 元胞大小對繞障距離變換結(jié)果的影響

本文元胞大小與地理信息系統(tǒng)中柵格的像元大小相對應(yīng)。由研究方法描述可知，本文方法并不受元胞大小的限制，任意元胞大小都可以使用。但元胞大小對結(jié)果運算速度和計算精度有影響。

從運算速度來看，對特定的研究區(qū)域，元胞越大，則該區(qū)域元胞數(shù)量越少，運算速度越快；元胞越小，則該區(qū)域劃分后的元胞數(shù)量越多，運算速度越慢。

從計算精度來看，最終的繞障歐氏距離變換結(jié)果等于元胞自動機模型模擬結(jié)果乘以元胞大小(分辨率)。因此，計算精度可由式(11)計算

ε=min(a,b)×Sc

(11)

式中，ε表示計算精度；a為距離算子上下左右位置的值；b為距離算子4角的值；Sc為元胞大小。

在式(11)中，距離算子的值相對固定，因此計算精度與元胞大小成正比。元胞越大，精度值越大，精度越低。元胞越小，精度值越小，精度越高。例如本文試驗中：a=0.961 94，b=1.360 39，Sc=5 km，因此計算精度為4.809 7 km。

運算速度和計算精度相互制約，在具體應(yīng)用中，像元大小的設(shè)置需在運算速度和計算精度中權(quán)衡選擇。

3.4.2 轉(zhuǎn)換規(guī)則對繞障距離變換結(jié)果的影響

圖5 距離算子對最大絕對錯誤率影響示意圖Fig.5 The influence of distance to maximum absolute error ratio

(1) 當(dāng)x

(2) 當(dāng)x=cos(π/8)時，圓O為八邊形的外接圓，最大絕對錯誤率體現(xiàn)在O與邊垂直方向。

(3) 當(dāng)cos(π/8)

(4) 當(dāng)x=1時，圓O為八邊形的內(nèi)切圓，最大絕對錯誤率體現(xiàn)在O與頂點連線方向。

(5) 當(dāng)x>1時，圓O在八邊形的內(nèi)部，最大絕對錯誤率體現(xiàn)在O與頂點連線方向。

通過幾何分析得出最大絕對錯誤率y與x的關(guān)系如下

(12)

進一步簡化為

(13)

從以上公式可看出，當(dāng)0[1+cos(π/8)]/2時，y=f(x)為遞增函數(shù)。當(dāng)x=[1+cos(π/8)]/2≈0.961 94時，得到y(tǒng)的最小值[1-cos(π/8)]/[1+cos(π/8)]≈0.396=3.96%。該取值即為本文使用的Butt-Maragos優(yōu)化算子。為此，本方法轉(zhuǎn)換規(guī)則的最優(yōu)算子為Butt-Maragos優(yōu)化算子，結(jié)果誤差≤3.96%。

4 結(jié)論與展望

本文以南海為例，基于海陸分布數(shù)據(jù)(障礙分布)和目標(biāo)點數(shù)據(jù)，提出了一種基于CA的繞障歐氏距離變換模型，模擬后得到變換結(jié)果，并對其繞障效果和精度進行分析。本文得到以下結(jié)論：①CA為繞障距離變換提供了一種解決方案，直觀動態(tài)地展示繞障歐氏距離變換的計算過程，是一種計算可視化的實現(xiàn)方式；②CA繞障距離變換模型在運行中能夠自動更新，兩目標(biāo)的距離擴散接縫處過渡自然；③受格網(wǎng)和鄰域的影響，基于CA的繞障歐氏距離變換結(jié)果為繞障至目標(biāo)最短路徑長度的近似值，采用Butt-Maragos優(yōu)化算子，誤差比例低于3.96%。計算結(jié)果可應(yīng)用于航線設(shè)計、海上救助等領(lǐng)域。

由于誤差的存在，如何通過對變換結(jié)果后續(xù)處理以改進基于CA的繞障歐氏距離變換的精度，有待進一步研究。