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

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

距離變換最早由Rosenfeld和Pfaltz提出，是指將每個柵格像元映射為到最近感興趣區(qū)域或目標區(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ū)域的影響，計算非障礙柵格像元到最近目標最短路徑長度(如航線設(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ù)的像元有相似之處，距離變換可以離散成為由目標像元向外擴散和更新繞障距離值的過程，因此可用元胞狀態(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ù)資料包括陸地島礁圖層和目標點圖層。陸地島礁圖層使用中國及東南亞的行政區(qū)劃Shapefile矢量數(shù)據(jù)，通過國家地理信息公共服務(wù)平臺服務(wù)資源中全球矢量地圖服務(wù)中對南海及鄰域矢量化陸地范圍得到，比例尺為1∶1000萬，坐標系統(tǒng)為2000國家大地坐標系。目標點圖層通過在研究區(qū)海域隨機選取7個點，用于分組測試不同數(shù)量目標情況下繞障歐氏距離變換效果。

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

2 研究方法

2.1 基本概念

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

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

(1)

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

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

(2)

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

(4) 繞障歐氏距離變換。變換結(jié)果使得任一非障礙集像元值為該像元到目標像元集的繞障歐氏距離。設(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標識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)定義為元胞中心位置至目標點集的繞障歐氏距離(格數(shù))，用dt表示。此外，元胞還具有位置(即所在列x,所在行y)、狀態(tài)變化標志(fs)，障礙區(qū)標志(fl)兩個屬性。位置為(x，y)元胞表示為C(x，y)。元胞屬性的具體描述見表1。

表1 元胞屬性說明

2.2.2 模型初始化

模型初始化即模型運行前的元胞屬性初始化與參數(shù)設(shè)置。模型的參數(shù)包括環(huán)境柵格和像元大小Sc。讀入環(huán)境柵格數(shù)據(jù)后賦予元胞的障礙標志屬性fl。Sc可以根據(jù)實際要求的精度確定，影響著元胞的大小與數(shù)量；狀態(tài)dt的初始值分為兩種情況，目標點所在元胞狀態(tài)dt設(shè)為0，其余元胞為+∞(比最大距離大的數(shù)，如w×h)；狀態(tài)變化標志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)。距離計算過程可以看作由目標像元向周圍逐漸擴散的過程，目標與周圍的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)(至目標最短路徑長度)組成的3×3矩陣為Dt，鄰域至中心元胞距離的3×3矩陣為A(稱為距離算子，具有中心對稱性，中心為0)。元胞下一狀態(tài)值(至目標最短路徑長度)即為矩陣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)容包括海陸分布、目標位置和元胞狀態(tài)，其中最核心的是元胞的展示。海陸分布作為背景，根據(jù)柵格值將海洋和陸地分別設(shè)置為藍色與淡黃色；目標位置采用圓點表示；元胞狀態(tài)的顯示采用漸變顏色表示，涉及顏色的映射，通過運行一遍模型后獲取最大的狀態(tài)值，再進行顏色映射以獲取良好的可視效果。

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

3 結(jié)果與討論

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

將目標點數(shù)據(jù)與預(yù)處理后的環(huán)境柵格輸入模型，創(chuàng)建250 000個元胞組成的CA模型，采用表2最后一個算子，模擬繞障歐氏距離變換過程，結(jié)果如圖2所示?？梢钥闯?，模型動態(tài)地展現(xiàn)出距離變換的計算過程。計算的元胞數(shù)量圍繞目標點呈正方形向外逐漸增多，顏色也越來越深，說明元胞狀態(tài)隨著擴散更新。當(dāng)兩個目標點的擴散圈相遇時(如圖2中的B位置)，接縫處(T=150的B位置)過渡自然，沒有出現(xiàn)橫豎形狀的突變狀態(tài)。這是由于計算過的元胞依然會根據(jù)鄰域情況作狀態(tài)更新。當(dāng)擴散至障礙區(qū)域時，會繞過陸地，元胞狀態(tài)值(如圖2中的A點)反映了繞障歐氏距離。在T=240時，表示模型經(jīng)過240次迭代，從目標經(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é)果的準確度，將CA模擬結(jié)果與普通歐氏距離變換結(jié)果對比。為避免其他因素影響，只選擇1個目標點，分別對研究區(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)被障礙隔開，而無法繞障到達目標點，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，其余位置誤差則與其與目標連線與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

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

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