基于弧段的遙感專題信息圖斑分段平滑方法

高 倩,章孝燦,胡 祺

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

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

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

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

針對上述問題，以矢量化后的圖斑為研究對象，考慮到相鄰圖斑的拓撲關系，以及圖斑平滑前后面積和位置保持的問題，提出基于弧段的遙感信息圖斑分段平滑方法。通過對矢量化后圖斑的弧段進行分段后再平滑，達到保持其拓撲關系，也保持原始圖斑面積和位置的結果。

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

基于弧段的遙感專題信息圖斑的分段平滑方法利用深度搜索的方式提取弧段，防止后續平滑出現拓撲不一致的結果，再根據鋸齒梯度對弧段分段，然后進行平滑，以降低矢量化后數據平滑造成的面積和位置的偏差。

1.1 弧段提取

柵格矢量化后的數據沒有拓撲結構，平滑時相鄰圖斑的公共邊平滑不一致，容易產生如圖1所示的“裂縫”，“重疊”等拓撲不一致的結果。但是在拓撲結構下，相鄰圖斑的公共邊共用一條弧段，僅平滑一次，則避免了平滑不一致產生的拓撲錯誤。為消除平滑時容易出現的拓撲錯誤，采取深度搜索的方式提取弧段，為后續平滑做準備。

圖1 有公共邊的平滑

實驗數據為常用ArcGIS的shp矢量文件，shp文件沒有拓撲結構，按數組順序存儲多邊形，多邊形外環的點順時針存儲，多邊形內環的點逆時針存儲。針對shp文件的特點，以及深度搜索方法[2-3, 16]的特點，提出以下方面修正，使其更適用于shp文件，流程如圖2所示。

圖2 弧段提取流程

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

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

圖3 圖斑中點的存儲

3)搜索到一段弧段時立即平滑替換掉原始弧段，并標記，再次搜索到該弧段時，不需要二次對比，直接跳過。

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

1.2 弧段分段平滑

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

1.2.1 平滑方式

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

圖4 鋸齒邊界分析

1.2.2 弧段分段

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

圖5 一般去鋸齒算法

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

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

步驟1：對弧段進行冗余點去除。遍歷弧段上所有點，連續3個點在一條直線上時，中間點為冗余點，需去除，然后進入步驟2。

步驟2：依次遍歷弧段上的點，計算相鄰4個點的關系，當條件：L(i-1,i)×L(i+1,i+2)>0,且|L(i,i+1)|的值為單位長時，分以下3種情況，不滿足條件時為線段標注0：

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

圖6 3種情況

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

圖7 弧段分段平滑

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

2 實驗與討論

將分辨率為30 m的Landsat8遙感影像經過非監督分類為5類，得到分類圖見圖8(a)，將其柵格矢量化后得到圖8(b)，以圖8(b)為實驗圖。

圖8 實驗數據

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

圖9 實驗結果

為了比較算法間的優缺點，實驗對比分析面積保持和位置保持效果以及平滑程度和壓縮率這4個方面的效果，統計結果見表1與圖10，分析方法如下：

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

2)面積保持效果分析。對平滑前后的面積作比較，計算各個圖斑的面積變化率，用平滑后的面積減去平滑前的面積得到面積差，面積差除以原始面積得到圖斑的面積變化率。實驗統計面積變化率的最大、最小值和均值，統計結果見圖10。

3)位置誤差分析。通過緩沖區限差[17]對位置誤差進行分析，緩沖區限差是對原始弧段建立緩沖區，統計平滑后弧段完全落入對應原始弧段緩沖區里的弧段數占總弧段數的比例，用此比例來評價位置的保持情況，比例值越大，位置保持越好。

緩沖區限差：

(1)

式中：Ln為落入相應原始弧段緩沖區的平滑后弧段數，L為總弧段數。

在統計落入緩沖區的弧段時，為保證相鄰弧段的緩沖區不會相交，對原始弧段建立半個像元距離的緩沖區。

4)壓縮程度。以平滑后的總點數與原始總點數之比的壓縮率來評價平滑前后的壓縮程度。

5)時間效率。統計4個方法計算同一數據的計算時長，以此來比較4個方法的時間效率。計算時長越短計算效率越高。

表1 4種方法平滑程度、位置誤差、壓縮率、時間效率的評價結果

圖10 面積變化率

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

3 結束語

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