基于弧段的遙感專(zhuān)題信息圖斑分段平滑方法

高 倩,章孝燦,胡 祺

(浙江大學(xué) 地球科學(xué)學(xué)院，浙江 杭州 310027)

柵格數(shù)據(jù)和矢量數(shù)據(jù)是地理信息系統(tǒng)中常用的兩種數(shù)據(jù)。柵格數(shù)據(jù)是一種離散的點(diǎn)陣數(shù)據(jù)，每一個(gè)單元的屬性值代表實(shí)體屬性，具有數(shù)據(jù)結(jié)構(gòu)簡(jiǎn)單，易于進(jìn)行空間分析等優(yōu)點(diǎn)，但是數(shù)據(jù)存儲(chǔ)量太大，縮放尺度有限制。矢量數(shù)據(jù)由點(diǎn)的精確坐標(biāo)來(lái)表達(dá)點(diǎn)、線、面等地理要素，具有數(shù)據(jù)存儲(chǔ)量少，放大縮小不失真，精確定位等優(yōu)點(diǎn)。由于矢量數(shù)據(jù)有存儲(chǔ)空間小，便于觀察，具有空間對(duì)象，可以精確定位等較于柵格數(shù)據(jù)的優(yōu)點(diǎn)，在實(shí)際生產(chǎn)中常常需要將柵格數(shù)據(jù)矢量化。柵格矢量化通常分為兩種，一種是線狀數(shù)據(jù)的矢量化，比如等值線圖，另一種是圖斑的矢量化，比如遙感分類(lèi)圖，專(zhuān)題圖等都需要柵格矢量化。本文主要討論圖斑矢量化后的問(wèn)題。柵格矢量化后，由于柵格數(shù)據(jù)具有點(diǎn)陣特點(diǎn)，是離散性數(shù)據(jù)，矢量化后的數(shù)據(jù)呈現(xiàn)鋸齒狀，十分不美觀，與實(shí)際地物的連續(xù)邊界不符，且造成大量的數(shù)據(jù)冗余。為了矢量化后圖斑的美觀性、真實(shí)性，以及減少數(shù)據(jù)的冗余，節(jié)省數(shù)據(jù)的存儲(chǔ)空間，需要對(duì)矢量化后的數(shù)據(jù)進(jìn)行平滑。

目前多邊形平滑算法可以分為3類(lèi)。一類(lèi)是矢量抽稀算法[1-8]，將數(shù)據(jù)抽稀，以達(dá)到減少數(shù)據(jù)冗余且平滑的結(jié)果。主要方法有垂距法，光欄法，道格拉斯普克(Douglas-Peucker , D-P)算法，小波算法，遺傳算法等。其中D-P算法是最常用的一種，自提出以來(lái)不斷被改進(jìn)，如黃萬(wàn)里等提出的基于面積保持的D-P算法[9]，顧騰等提出的D-P算法與Li-Openshaw結(jié)合改進(jìn)的方法[7]。雖然矢量抽稀算法能有效平滑鋸齒，消除數(shù)據(jù)冗余，但是往往需要一個(gè)適合的閾值，閾值的大小難以把握。第二類(lèi)是曲線擬合算法[10-12]，如多項(xiàng)式擬合，B樣條擬合，貝塞爾曲線等，曲線擬合是平滑多邊形最常用方法。然而曲線擬合難以去除圖斑鋸齒，反而會(huì)出現(xiàn)“波浪”。第三類(lèi)是綜合法，將多類(lèi)方法結(jié)合使用，比如前面提到的D-P算法和Li-Openshaw結(jié)合的方法，D-P算法和B樣條結(jié)合的方法[13]等，此類(lèi)方法因結(jié)合多種方法計(jì)算緩慢。

由于圖斑中鋸齒的存在，上述3種多邊形平滑方法無(wú)法較好地平滑圖斑。加之遙感圖斑通常具有實(shí)際意義，特別是土地利用類(lèi)型圖斑的面積屬性，在實(shí)際應(yīng)用中有重要意義，平滑后圖斑與原始圖斑面積和位置的保持十分重要。此外,若相鄰圖斑公共邊處理不一致，容易出現(xiàn)“裂縫”等拓?fù)溴e(cuò)誤的結(jié)果。因此，需要提出具有針對(duì)性的圖斑平滑方法。石軍南等人提出的矢量線條概率中值平滑方法[14]，宋正祥等人提出的顧及拓?fù)渑c尖角的矢量數(shù)據(jù)分組壓縮算法都是針對(duì)圖斑的平滑[15]，均達(dá)到了良好的平滑效果，但無(wú)法保持原始圖斑面積和位置。綜上圖斑的平滑存在3個(gè)問(wèn)題: ①圖斑呈現(xiàn)鋸齒狀；②平滑時(shí)難以保持圖斑面積及位置信息；③平滑時(shí)圖斑公共邊易出現(xiàn)拓?fù)溴e(cuò)誤。

針對(duì)上述問(wèn)題，以矢量化后的圖斑為研究對(duì)象，考慮到相鄰圖斑的拓?fù)潢P(guān)系，以及圖斑平滑前后面積和位置保持的問(wèn)題，提出基于弧段的遙感信息圖斑分段平滑方法。通過(guò)對(duì)矢量化后圖斑的弧段進(jìn)行分段后再平滑，達(dá)到保持其拓?fù)潢P(guān)系，也保持原始圖斑面積和位置的結(jié)果。

1 基于弧段的分段平滑方法

基于弧段的遙感專(zhuān)題信息圖斑的分段平滑方法利用深度搜索的方式提取弧段，防止后續(xù)平滑出現(xiàn)拓?fù)洳灰恢碌慕Y(jié)果，再根據(jù)鋸齒梯度對(duì)弧段分段，然后進(jìn)行平滑，以降低矢量化后數(shù)據(jù)平滑造成的面積和位置的偏差。

1.1 弧段提取

柵格矢量化后的數(shù)據(jù)沒(méi)有拓?fù)浣Y(jié)構(gòu)，平滑時(shí)相鄰圖斑的公共邊平滑不一致，容易產(chǎn)生如圖1所示的“裂縫”，“重疊”等拓?fù)洳灰恢碌慕Y(jié)果。但是在拓?fù)浣Y(jié)構(gòu)下，相鄰圖斑的公共邊共用一條弧段，僅平滑一次，則避免了平滑不一致產(chǎn)生的拓?fù)溴e(cuò)誤。為消除平滑時(shí)容易出現(xiàn)的拓?fù)溴e(cuò)誤，采取深度搜索的方式提取弧段，為后續(xù)平滑做準(zhǔn)備。

圖1 有公共邊的平滑

實(shí)驗(yàn)數(shù)據(jù)為常用ArcGIS的shp矢量文件，shp文件沒(méi)有拓?fù)浣Y(jié)構(gòu)，按數(shù)組順序存儲(chǔ)多邊形，多邊形外環(huán)的點(diǎn)順時(shí)針存儲(chǔ)，多邊形內(nèi)環(huán)的點(diǎn)逆時(shí)針存儲(chǔ)。針對(duì)shp文件的特點(diǎn)，以及深度搜索方法[2-3, 16]的特點(diǎn)，提出以下方面修正，使其更適用于shp文件，流程如圖2所示。

圖2 弧段提取流程

1)深度搜索采取“順序比較”和“逆序比較”的方法搜索公共邊。在shp文件中，公共邊在相鄰圖斑中必然是順序相反，因此只需要逆序比較。

2)搜索時(shí)起始點(diǎn)為第一個(gè)存儲(chǔ)點(diǎn)，不是弧段的端點(diǎn)。如圖3所示，圖斑AB在公共邊均有圖斑存儲(chǔ)的起始點(diǎn)與終止點(diǎn)，于是除了“逆序比較”找到end點(diǎn)，還要將start點(diǎn)向前移動(dòng)。即A多邊形的點(diǎn)1，2，3，…，比較B多邊形i,i-1,i-2，…，找到end點(diǎn)，還需要A多邊形的點(diǎn)n，n-1，n-2,…,比較B多邊形i,i+1,i+2,…,找到真正的start點(diǎn)。

圖3 圖斑中點(diǎn)的存儲(chǔ)

3)搜索到一段弧段時(shí)立即平滑替換掉原始弧段，并標(biāo)記，再次搜索到該弧段時(shí)，不需要二次對(duì)比，直接跳過(guò)。

4)弧段僅需平滑一次，若是公共邊的弧段，在另一圖斑中只需將平滑后弧段逆序存儲(chǔ)替換掉原始弧段。

1.2 弧段分段平滑

柵格矢量化后，由于遙感圖像的點(diǎn)陣特點(diǎn)，柵格數(shù)據(jù)對(duì)邊界的表達(dá)有限，矢量化后的圖斑存在大量鋸齒。圖斑不美觀和數(shù)據(jù)冗余都由鋸齒造成，于是平滑圖斑只需對(duì)鋸齒平滑，即可達(dá)到平滑效果。采取將弧段分段為平滑鋸齒單元的方式找到需要平滑的鋸齒，采用鋸齒中線代替鋸齒的方式進(jìn)行平滑，達(dá)到既平滑圖斑，又保持圖斑面積和位置的結(jié)果。

1.2.1 平滑方式

如圖4所示AB兩圖斑，公共邊為鋸齒。提到柵格數(shù)據(jù)轉(zhuǎn)換為矢量數(shù)據(jù)時(shí)，線條形狀變化，其特點(diǎn)分布概率不確定，柵格數(shù)據(jù)矢量化的曲線是一條柵格線條的中線時(shí)方差最小。在圖斑中同樣，如圖4所示，鋸齒為AB兩圖斑的交接處，線1和線3包圍區(qū)域?yàn)殇忼X所屬范圍，同時(shí)為既有圖斑A也有圖斑B的混淆地帶。在線1和線3包圍區(qū)域內(nèi)的線都能將此區(qū)域等分為兩份來(lái)替代鋸齒，但是中間線線2為鋸齒邊最佳擬合的直線，其殘差最小。線2可使該交叉地帶一分為二，使斑AB位置不偏移，面積也能保持一致。不難發(fā)現(xiàn)，若需要平滑后AB兩圖斑面積保持一致，中線需平分鋸齒，為保證左右分到的面積相等，中線的起始點(diǎn)在水平線(垂直線)上，終止點(diǎn)也應(yīng)在水平線(垂直線)上。如圖4所示，可以保證AB圖斑面積不變，這樣中線左右兩邊的點(diǎn)數(shù)相同，面積也相同。

圖4 鋸齒邊界分析

1.2.2 弧段分段

1)平滑鋸齒單元。單純地使用中線進(jìn)行鋸齒平滑，保持原狀的邊也進(jìn)行處理，且未考慮到平滑后圖斑面積保持的問(wèn)題。如圖5(b)中弧段的2～3處，此處數(shù)據(jù)不進(jìn)行處理，也能去除鋸齒且更能保持原始數(shù)據(jù)。為了減小平滑后數(shù)據(jù)的偏差，在平滑時(shí)此類(lèi)數(shù)據(jù)不應(yīng)平滑。由此，提出將鋸齒弧段分段為平滑鋸齒單元，平滑鋸齒單元為連續(xù)的相同梯度的鋸齒，僅對(duì)平滑鋸齒單元進(jìn)行平滑處理，這樣就能在去除鋸齒的同時(shí)盡量減少平滑對(duì)數(shù)據(jù)造成的偏差。

圖5 一般去鋸齒算法

相關(guān)概念定義如下，平滑鋸齒單元：相同梯度的鋸齒邊分割為一個(gè)平滑鋸齒單元，如圖5(b)中弧段的1～2處為一個(gè)平滑鋸齒單元。鋸齒邊：相鄰三邊，中間邊為單位像元長(zhǎng)，前后兩邊走向一致，則為鋸齒邊，如圖6(a)所示。鋸齒梯度：相鄰三邊，中間邊為單位像元長(zhǎng)，前一條邊與中間邊長(zhǎng)度的比值為鋸齒梯度。

2)平滑鋸齒單元分割方法。由上述思想，提出一種平滑鋸齒單元分割方法。設(shè)弧段上點(diǎn)的坐標(biāo)為(Xi,Yi)，相鄰點(diǎn)橫坐標(biāo)差值|Xi=Xi+1|X，縱坐標(biāo)差值|Yi=Yi+1|Y。因?yàn)殇忼X弧段都是由水平或者豎直線段組成，所以相鄰兩點(diǎn)有且僅有一個(gè)為0，因此相鄰兩點(diǎn)組成的線段長(zhǎng)度L(i,i+1)=|ΔXi|+Yi,，且L(i,i+1)的正負(fù)代表其走向。平滑鋸齒單元分割方法如下：

步驟1：對(duì)弧段進(jìn)行冗余點(diǎn)去除。遍歷弧段上所有點(diǎn)，連續(xù)3個(gè)點(diǎn)在一條直線上時(shí)，中間點(diǎn)為冗余點(diǎn)，需去除，然后進(jìn)入步驟2。

步驟2：依次遍歷弧段上的點(diǎn)，計(jì)算相鄰4個(gè)點(diǎn)的關(guān)系，當(dāng)條件：L(i-1,i)×L(i+1,i+2)>0,且|L(i,i+1)|的值為單位長(zhǎng)時(shí)，分以下3種情況，不滿(mǎn)足條件時(shí)為線段標(biāo)注0：

②如圖6(b)所示，若|L(i+1,i)

圖6 3種情況

步驟3：標(biāo)注完成后，當(dāng)標(biāo)注連續(xù)數(shù)字相同且大于0的邊分為一個(gè)平滑鋸齒單元。如圖7(a)弧段可以分為三段，平滑鋸齒單元1～2，平滑鋸齒單元3～4，非平滑鋸齒單元2～3。

圖7 弧段分段平滑

步驟4：平滑。平滑鋸齒單元進(jìn)行平滑處理。非平滑鋸齒單元保持原始數(shù)據(jù)。如圖7(b)所示弧段分段平滑后。

2 實(shí)驗(yàn)與討論

將分辨率為30 m的Landsat8遙感影像經(jīng)過(guò)非監(jiān)督分類(lèi)為5類(lèi)，得到分類(lèi)圖見(jiàn)圖8(a)，將其柵格矢量化后得到圖8(b)，以圖8(b)為實(shí)驗(yàn)圖。

圖8 實(shí)驗(yàn)數(shù)據(jù)

為了表明方法的有效性，實(shí)驗(yàn)對(duì)比3種用在圖斑平滑中的算法，分別是概率中值算法[14]、D-P算法[2]( D-P算法使用的壓縮閾值為30)、去尖角算法[15]。3種算法及文中提出的方法均在Visual Studio 2010下用C++語(yǔ)言利用GDAL庫(kù)編程實(shí)現(xiàn)，對(duì) shp 文件格式的矢量實(shí)驗(yàn)圖進(jìn)行處理。平滑結(jié)果如圖9所示。

圖9 實(shí)驗(yàn)結(jié)果

為了比較算法間的優(yōu)缺點(diǎn)，實(shí)驗(yàn)對(duì)比分析面積保持和位置保持效果以及平滑程度和壓縮率這4個(gè)方面的效果，統(tǒng)計(jì)結(jié)果見(jiàn)表1與圖10，分析方法如下：

1)平滑程度分析。平滑程度除了可以從目視效果看出，還可以從尖角(0°～90°或者270°～360°的角)數(shù)量得出，尖角數(shù)越少越平滑。

2)面積保持效果分析。對(duì)平滑前后的面積作比較，計(jì)算各個(gè)圖斑的面積變化率，用平滑后的面積減去平滑前的面積得到面積差，面積差除以原始面積得到圖斑的面積變化率。實(shí)驗(yàn)統(tǒng)計(jì)面積變化率的最大、最小值和均值，統(tǒng)計(jì)結(jié)果見(jiàn)圖10。

3)位置誤差分析。通過(guò)緩沖區(qū)限差[17]對(duì)位置誤差進(jìn)行分析，緩沖區(qū)限差是對(duì)原始弧段建立緩沖區(qū)，統(tǒng)計(jì)平滑后弧段完全落入對(duì)應(yīng)原始弧段緩沖區(qū)里的弧段數(shù)占總弧段數(shù)的比例，用此比例來(lái)評(píng)價(jià)位置的保持情況，比例值越大，位置保持越好。

緩沖區(qū)限差：

(1)

式中：Ln為落入相應(yīng)原始弧段緩沖區(qū)的平滑后弧段數(shù)，L為總弧段數(shù)。

在統(tǒng)計(jì)落入緩沖區(qū)的弧段時(shí)，為保證相鄰弧段的緩沖區(qū)不會(huì)相交，對(duì)原始弧段建立半個(gè)像元距離的緩沖區(qū)。

4)壓縮程度。以平滑后的總點(diǎn)數(shù)與原始總點(diǎn)數(shù)之比的壓縮率來(lái)評(píng)價(jià)平滑前后的壓縮程度。

5)時(shí)間效率。統(tǒng)計(jì)4個(gè)方法計(jì)算同一數(shù)據(jù)的計(jì)算時(shí)長(zhǎng)，以此來(lái)比較4個(gè)方法的時(shí)間效率。計(jì)算時(shí)長(zhǎng)越短計(jì)算效率越高。

表1 4種方法平滑程度、位置誤差、壓縮率、時(shí)間效率的評(píng)價(jià)結(jié)果

圖10 面積變化率

觀察上述實(shí)驗(yàn)結(jié)果，圖9可以清楚觀察到4個(gè)算法平滑結(jié)果，僅目視效果而言，D-P算法呈現(xiàn)大量尖角，并不是很美觀。去尖角算法雖然去除了D-P算法的大部分尖角，但是由不同公共邊共同形成的尖角卻未能去除。通過(guò)表1中的對(duì)尖角的統(tǒng)計(jì)，可以看出分段平滑方法平滑程度最高，D-P算法平滑程度最低。由圖10 和表1中面積和位置的統(tǒng)計(jì)結(jié)果可以看到，分段平滑方法平滑前后面積保持了一致，而且弧段均在緩沖區(qū)內(nèi)，分段平滑方法面積與位置保持是最優(yōu)的。D-P算法是位置和面積保持最差的，但是在壓縮程度上D-P算法有著絕對(duì)的優(yōu)勢(shì)?？梢?jiàn)D-P算法對(duì)圖斑的壓縮平滑，在面積和位置信息上有所犧牲。另外兩個(gè)算法的壓縮率和面積位置保持都介于D-P算法和分段平滑方法之間。時(shí)間效率上去尖角算法是在D-P算法基礎(chǔ)上進(jìn)行了尖角去除，時(shí)間效率上比D-P低，本文算法較另外3種方法稍稍高一點(diǎn)。綜上所述,幾種方法雖然都能有效平滑，但是在平滑程度，面積保持，位置保持以及計(jì)算效率方面文中提出的方法最優(yōu)。

3 結(jié)束語(yǔ)

本文針對(duì)遙感信息圖斑邊界呈現(xiàn)鋸齒狀，提出一種基于弧段的遙感專(zhuān)題信息圖斑分段平滑方法，用深度搜索的方式提取圖斑弧段，根據(jù)鋸齒弧段的梯度分段為平滑鋸齒單元，對(duì)平滑鋸齒單元進(jìn)行平滑處理，以達(dá)到平滑弧段、圖斑面積與位置保持之間平衡，為遙感分類(lèi)圖，專(zhuān)題圖等矢量化后提供一種平滑方法。該方法減輕了以往方法平滑后對(duì)原始數(shù)據(jù)造成的偏差，但平滑后數(shù)據(jù)冗余消除程度遠(yuǎn)不及矢量抽稀算法。