基于子像元交換算法和地形數據的土地覆蓋亞像元制圖研究

虞舟魯， 王文超， 戎 奕， 沈掌泉

(1.浙江大學環境與資源學院,杭州 310058； 2.杭州學聯土地規劃設計咨詢有限公司,杭州 310030)

基于子像元交換算法和地形數據的土地覆蓋亞像元制圖研究

虞舟魯1， 王文超2， 戎 奕2， 沈掌泉1

(1.浙江大學環境與資源學院,杭州 310058； 2.杭州學聯土地規劃設計咨詢有限公司,杭州 310030)

在將遙感影像應用于土地覆蓋制圖的過程中，混合像元會產生負面影響，尤其是在空間分辨率較低的情況下。軟分類和超分辨率(亞像元/子像元)制圖技術可以解決上述問題。子像元交換算法是一種簡單而有效的亞像元制圖技術，但也存在計算效率不夠高和在超分辨率制圖因子(scale factor,S)較大時制圖精度較差的問題，其原因可能是亞像元制圖只是從軟分類結果中獲得信息。因此，將數字高程模型(digital elevation model,DEM)及其衍生的數據作為子像元交換算法的輔助信息，研究提高其制圖精度的有效性。研究結果表明： ①在DEM信息的輔助下，亞像元制圖的精度得到了改善，即使在S較大時也是有效的； ②同時使用高程和坡度信息時的效果要好于單獨應用高程或坡度信息； ③在S較大的情況下，制圖精度對鄰域范圍大小的敏感性要低于S較小時； ④在DEM信息的輔助下，計算效率得到提高。實驗結果表明，將DEM作為輔助信息進行亞像元制圖是有效和可行的。

亞像元制圖； 超分辨率制圖； 子像元交換算法； 遙感； 數字高程模型(DEM)； 土地覆蓋制圖

0 引言

土地覆蓋是生態系統中生物化學循環、水文變化及地表與大氣之間相互作用的體現[1]，其信息可用于理解和管理生態環境。因此，獲得準確和現勢的土地覆蓋信息在土地管理、規劃和監測中是非常重要的。由于遙感技術具有大范圍獲取信息的能力，因此遙感影像是獲得土地覆蓋信息的有效來源[2]。但是，衛星遙感影像普遍存在混合像元的問題，特別是在空間分辨率較低的情況下[3]。因此，在利用遙感數據提取土地覆蓋信息時，混合像元常成為一個不容忽視的問題。在處理混合像元時，與硬分類技術相比，軟分類技術可以避免信息的損失，因為其可以獲得各土地覆蓋類型在像元中所占的比例，并形成一系列對應各土地覆蓋類型的比例圖像[4-5]。可是，軟分類形成的比例圖像并沒有明確各類型在像元內的空間位置，而在很多應用中，確定各類型在像元內的空間位置并獲得詳細的土地覆蓋圖是非常必要和重要的。Atkinson[6]提出的亞像元制圖(亦稱子像元制圖、超分辨率制圖或銳化)的概念認為，根據軟分類所獲得的比例圖像形成銳化的土地覆蓋圖，可以綜合發揮軟分類和硬分類的優勢。目前，學者們已經提出了一系列的亞像元制圖技術，諸如子像元交換算法[7-8]、圖像銳化技術[9-10]、基于知識的分析技術[11]、Hopfield神經網絡[1,12-13]、反卷積濾波器[14]、線性優化技術[15]、遺傳算法[4]、前饋神經網絡[16]、基于馬爾可夫隨機場的技術[2]、基于子像元/像元空間吸引模型的技術[17]、基于子像元/像元空間吸引模型的改進子像元交換算法[18]、子像元/像元空間吸引模型與單親遺傳算法相結合的技術[19]等。子像元交換算法具有算法簡單和效率較高的特點，已成功地應用于馬來西亞的海岸線制圖[20]和英國Dorset的Christchurch地區的農村土地覆蓋制圖[21]； 但是當超分辨率因子(scale factor, S)較大時(即在空間分辨率較低的情況下)，其計算效率和制圖精度就會下降，其原因可能是只依靠軟分類結果中各類型的比例信息和子像元類型的隨機初始化所導致的大量子像元交換。另外，子像元類型的隨機初始化也會對最終的制圖結果產生影響。因此，通過輔助信息來改進該算法的初始化過程，有可能提高其計算效率和制圖精度。

地形是基本的地球物理形態，包含了一個地區地殼變化、地質構造和氣候變化等豐富的信息，在土地覆蓋分類中經常作為輔助的信息[22]。詳細和準確的數字高程模型(digital elevation model,DEM)信息可通過SPOT及ASTER等遙感數據獲取，也可通過機載激光掃描系統(light detecting and ranging,LiDAR)獲得，美國國家航空航天局(National Aeronautics and Space Administration,NASA)則通過航天飛機雷達地形測繪任務(shuttle Radar topography mission,SRTM)獲得了全球大部分陸地的地形信息。Nguyen等[22]提出將LiDAR獲得的DEM作為輔助信息通過Hopfield神經網絡(Hopfield neural network,HNN)進行超分辨率制圖，其研究結果表明，采用LiDAR獲得的DEM數據結合光學遙感數據，HNN可以在子像元尺度上準確地預測土地覆蓋類型，使制圖精度明顯提升，尤其是建筑物制圖。其應用的最大問題是在研究區需要有可用的LiDAR數據，否則無法將其應用于通過軟分類獲得的比例圖像結果中。基于上述情況，本文通過統計分析地形數據與土地覆蓋類型之間的關系，并將地形數據作為輔助信息用于子像元交換算法的初始化過程，來代替原算法中的隨機初始化過程，以研究將地形作為輔助信息進行基于子像元交換算法的超分辨率制圖的有效性和可行性。

1 研究區概況與數據源

研究區位于浙江省蘭溪市東北部，研究區土地覆蓋、高程和坡度如圖1所示。

(a) 土地覆蓋(b) 高程 (c) 坡度

圖1研究區土地覆蓋、高程和坡度

Fig.1Landcover,elevationandslopeofstudyarea

研究區面積約為420 hm2，覆蓋研究區的航空遙感影像包含512像元×512像元，空間分辨率為4 m。土地覆蓋類型圖來自于對航空遙感影像的目視解譯。土地覆蓋類型主要有耕地、林地、建設用地和水域。通過對1∶5萬比例尺地形圖數字化后建立不規則三角網(triangulated irregular network,TIN)獲得DEM數據，然后利用DEM計算得到坡度數據(以上運算均在ArcGIS9.3的spatial analysis擴展模塊中完成)。

在本文中，將由子像元制圖的參考圖像生成比例圖像的過程稱為退化，各像元的土地覆蓋類型的比例是以參考圖像為依據、按照S計算得到[4,17]。計算得到的比例圖像，作為子像元交換算法的輸入。以參考的土地覆蓋類型圖為基礎，通過退化獲得S=2,4,8,16，32這5組比例圖像。硬分類是傳統分類的結果，以比例圖像為基礎，將比例最大的土地覆蓋類型作為相應像元所有子像元的類型。圖2與圖3分別為研究區不同土地覆蓋類型的比例圖像和相應硬分類結果的示例，圖上所示為S=8時的結果。

(a) 耕地(b) 林地 (c) 建設用地 (d) 水域

圖2比例圖像示例

Fig.2Exampleoffractionimages

圖3 硬分類結果示例Fig.3 Example of hard classifcation result

2 實驗算法

2.1 子像元交換算法

子像元交換算法的基本思想是通過交換像元內子像元的空間位置使相鄰子像元之間的空間相關性最大。在算法的開始，按照各類型的比例，隨機賦予像元中各子像元的土地覆蓋類型。在初始化完成后，各像元中子像元的空間位置可以改變，而土地覆蓋類型的比例保持不變； 按照鄰近子像元之間及像元的空間相關性最大化的原則，改變像元中子像元的空間位置。該算法包含3個步驟：

1)對于每個像元，按照與距離成反比的規則計算每個子像元的空間吸引力。假設pi,j為像元Pa,b的一個子像元，pk為pi,j相鄰的子像元，pk屬于Pa,b或其相鄰的像元，則Pa,b的空間吸引力OPa,b為

(1)

式中：S為超分辨率制圖因子，每個像元包含S2個子像元；Opi,j為第i行、第j列子像元的空間吸引力，以其鄰域中的子像元為基礎，可得

(2)

式中：N為鄰域的子像元數；Z(pi,j,pk)為二值函數，若子像元pi,j與其鄰域中的子像元pk的土地覆蓋類型相同，則其值為1，否則為0；λk為權重系數，與距離有關，即

(3)

式中：a為基于對數的非線性系數；h(pi,j,pk)為子像元pi,j與其鄰域中的子像元pk之間的距離，即

(4)

2)以當前像元的子像元空間位置為基礎，逐像元進行評估。

3)在同一像元中，選擇2個具有不同土地覆蓋類型且空間吸引力最小的子像元，如果交換其位置后該像元的空間吸引力增加，則交換其位置； 否則不進行交換。

對上述3個步驟進行迭代，直到子像元的空間吸引力不再增加時，運算過程停止。

Atkinson[7-8]提出的子像元交換算法是針對2個土地覆蓋類型的情況，而Thornton等[21]則將其擴展到可用于多個土地覆蓋類型的超分辨率制圖。

2.2 基于DEM的初始化過程

將從DEM獲得的高程和坡度數據用于像元中子像元的初始化，以代替原算法中的隨機初始化過程。在本研究中，還分析了只用高程數據、只用坡度數據和同時使用高程與坡度數據3種情況。初始化過程包括4個步驟：

1)根據DEM數據和對應的土地覆蓋類型統計二者的關系。由于高程和坡度為連續值，因此需要將其轉換為類型值，在ArcMap的spatial analysis擴展模塊中應用natural break將高程和坡度分為15個子區間，各子區間的范圍如表1所示。

表1 各子區間高程和坡度數值范圍Tab.1 Data ranges of elevation and slope for partitioning

對分區后的高程和坡度數據與土地覆蓋類型之間的關系進行統計，假設土地覆蓋類型為i，高程或坡度的子區間為j，則它們之間的關系Cij為

(5)

式中nij為土地覆蓋類型i落在第j個高程或坡度子區間的子像元個數。因此，c是一個二維或三維矩陣，一個維度是高程和(或)坡度，另一個維度為土地覆蓋類型。圖4示出通過參考圖像數據和高程/坡度統計得到的cij。

(a) 土地覆蓋類型與高程的關系 (b) 土地覆蓋類型與坡度的關系

(c)各土地覆蓋類型與高程和坡度的關系圖4 土地覆蓋類型與DEM信息的相互關系Fig.4 Interrelationships between land cover types and DEM data

從圖4可以看出，不同的土地覆蓋類型與高程、坡度子區間之間存在明顯的差異。

（2）當時，為通項遞減的正項級數．因為，所以當時，；當時，．由定理2可知，當時，級數收斂；當時，級數發散．

2)以高程或坡度數據為基礎，獲得像元中子像元各土地覆蓋類型的置信值。假定pi,j是像元Pa,b中的一個子像元，根據pi,j的高程或坡度與土地覆蓋類型之間的關系cij，可以知道對于土地覆蓋類型l的置信值是cl,pi,j。

3)將置信值進行歸一化，即

(6)

以DEM數據為基礎，逐像元地進行初始化。

2.3 子像元交換算法參數

(a)N=1(b)N=2(c)N=4 (d)N=6(e)N=8

(S為4； 中間黑色點為子像元pi,j； 其他灰色點為其鄰域的子像元； 粗線條為像元范圍)

圖5鄰域大小定義的示例

Fig.5Illustrationofdefinitionofneighborhoodsize

2.4 算法比較與評價標準

分別采用調整的Kappa值、所需CPU時間、迭代次數和總交換次數4個指標來比較和評價算法的表現。

調整的Kappa值通過計算結果與參考圖像的比較獲得，用于評價超分辨率制圖的精度。與Kappa值不同的是只統計混合像元，因此對亞像元制圖精度的描述更敏感和更合理。所需CPU時間用于比較算法的計算效率，在相同計算環境下獲得。迭代次數和總交換次數也用于衡量算法的表現和效率，如在相同參數下，執行的迭代次數多，則所需的CPU時間就長； 同樣，如果總交換次數少而制圖精度高，則說明初始化過程的效果好、計算效率高。

在本研究中，各算法是由MATLAB的腳本m語言實現的，為了使計算結果更為可靠和更具可比性，對每組參數均重復5次，然后進行平均。

3 實驗結果與分析

表2示出本文提出的初始化算法與原算法中隨機初始化算法的超分辨率制圖精度與計算所需的CPU時間。

表2 各初始化算法的調整Kappa值和計算所需CPU時間Tab.2 Adjusted Kappa values derived by different initializing algorithms and their consuming CPU time

①HC為硬分類； ②R為隨機初始化； ③H為基于高程的初始化； ④S為基于坡度的初始化； ⑤H&S為基于高程和坡度的初始化。

從表2可知，在S≤8的情況下，本文算法制圖精度要低于硬分類； 而在S>8時，其精度高于硬分類。除S=2外，本文算法的制圖精度均好于隨機初始化算法； 在S>8時，優勢更為明顯。在3種DEM輔助初始化算法中，S>8時高程數據輔助初始化算法的精度最高。與隨機初始化算法相比，DEM輔助初始化所需的處理時間要長，而且隨著S的加大，所需時間也明顯增加。

圖6是不同初始化算法的結果，圖中線條表示參考的土地覆蓋類型圖中各土地覆蓋類型邊界。

(a) 隨機初始化

(b) 基于高程輔助的初始化

圖6-1 不同初始化算法5種S的制圖結果Fig.6-1 Sub-pixel mapping results with 5 degraded S generated by different initializing algorithms

(c) 基于坡度輔助的初始化

(d) 基于高程和坡度輔助的初始化

圖6-2 不同初始化算法5種S的制圖結果Fig.6-2 Sub-pixel mapping results with 5 degraded S generated by different initializing algorithms

通過目視比較可以看出，基于DEM輔助初始化算法的制圖結果要好于隨機初始化算法。因此，應用基于DEM輔助的初始化算法可以為后續的子像元交換提供更好的基礎。圖7示出基于子像元交換算法的制圖精度與硬分類和原算法制圖精度的對比。在S>4時，本文算法能取得更好的制圖精度。在DEM輔助的3種算法中，應用高程輔助算法的精度最好，而僅基于坡度輔助算法的精度最低。基于DEM信息輔助算法的不足在于，當S≤8時，對鄰域大小這個參數的變化比較敏感。

圖7 不同算法的子像元制圖精度Fig.7 Sub-pixel mapping accuracy for different algorithms

圖8示出不同參數情況下各算法子像元制圖所需的CPU時間。

圖8不同算法子像元制圖所需的CPU時間

Fig.8CPUtimeinsub-pixelmappingcomputationwithdifferentalgorithms

盡管初始化過程所需的計算時間較多，但基于DEM輔助初始化的算法過程所需CPU時間要遠少于原算法； 在基于DEM輔助的3種算法中，基于高程輔助算法所需的CPU時間最少，而基于坡度輔助算法所需的CPU時間最多。其原因是本文提出的算法無論是迭代次數還是總交換次數均明顯少于原算法(圖9和圖10)。說明本文算法可以提高制圖精度，同時提高運算效率。

圖9不同算法子像元制圖運算過程中的迭代次數

Fig.9Numberofiterationinsub-pixelmappingcomputationwithdifferentalgorithms

圖10 不同算法子像元制圖運算過程中的總交換次數Fig.10 Total swapping numbers in sub-pixel mapping computation with different algorithms

圖11示出硬分類和4種子像元交換算法的超分辨率制圖的最終結果。

(a) 硬分類

(b) 基于隨機初始化的子像元交換算法

(c) 基于高程輔助初始化的子像元交換算法

(d) 基于坡度輔助初始化的子像元交換算法

圖11-1不同算法子像元制圖最終結果

Fig.11-1Sub-pixelmappingresultswith5degradedscalefactorsderivedbydifferentalgorithms

(e) 基于高程和坡度輔助初始化的子像元交換算法

圖11-2不同算法子像元制圖最終結果

Fig.11-2Sub-pixelmappingresultswith5degradedscalefactorsderivedbydifferentalgorithms

從圖11可以看出，在S=32時，所有子像元交換算法的結果均好于硬分類，但本文算法獲得的結果更加合理； 在大部分S值下，特別是當S較大時，本文提出的算法均好于原算法。

4 結論

盡管子像元交換算法是一種簡單、有效的亞像元制圖算法，但在超分辨率制圖因子(S)較大時，仍然存在計算量大和制圖精度較低的問題。本研究基于地形信息和土地覆蓋類型之間的內在聯系，提出了一種新的算法，將原算法中的隨機初始化過程替換為基于DEM輔助的初始化過程。研究結果表明，本文提出的算法可以獲得更好的制圖結果，尤其是在S較大的情況下，優勢更加明顯； 而且，因為減少了迭代次數和總交換次數，其計算效率更高，所需的計算時間更少。在以地形作為輔助信息時，無論是制圖精度還是計算效率，基于高程輔助初始化的效果都最好。

本文所提出算法的不足是，在應用本算法前，需要先獲得地形信息與土地覆蓋類型之間的關系，因此初始化過程和制圖精度可能受其影響，對此需要進一步研究和評估。

志謝： 感謝美國密歇根州立大學全球變化與對地觀測中心Qi Jiaguo教授在研究過程中給予的指導和幫助。

[1] Tatem A J,Lewis H G,Atkinson P M,et al.Super-resolution land cover pattern prediction using a Hopfield neural network[J].Remote Sensing of Environment,2002,79(1):1-14.

[2] Kasetkasem T,Arora M K,Varshney P K.Super-resolution land cover mapping using a Markov random field based approach[J].Remote Sensing of Environment,2005,96(3/4):302-314.

[3] Foody G M.Remote Sensing Image Analysis:Including the Spatial Domain[M].Norwell,MA:Kluwer,2004:37-49.

[4] Mertens K C,Verbeke L P C,Ducheyne E I,et al.Using genetic algorithms in sub-pixel mapping[J].International Journal of Remote Sensing,2003,24(21):4241-4247.

[5] 任 武,葛 詠.遙感影像亞像元制圖方法研究進展綜述[J].遙感技術與應用,2011,26(1):33-44.

Ren W,Ge Y.Progress on sub-pixel mapping methods for remotely sensed images[J].Remote Sensing Technology and Application,2011,26(1):33-44.

[6] Atkinson P M.Mapping sub-pixel boundaries from remotely sensed images[M]//Kemp Z.Innovations in GIS IV.London:Taylor and Francis,1997:166-180.

[7] Atkinson P M.Super-resolution target mapping from soft-classified remotely sensed imagery[C]//Proceedings of the 5th International Conference on GeoComputation.Leeds,UK:University of Leeds,2001.

[8] Atkinson P M.Sub-pixel target mapping from soft-classified,remotely sensed imagery[J].Photogrammetric Engineering and Remote Sensing,2005,71(7):839-846.

[9] Foody G M.Sharpening fuzzy classification output to refine the representation of sub-pixel land cover distribution[J].International Journal of Remote Sensing,1998,19(13):2593-2599.

[10] Grossa H N,Schotta J R.Application of spectral mixture analysis and image fusion techniques for image sharpening[J].Remote Sensing of Environment,1998,63(2):85-94.

[11] Schneider W.Land use mapping with subpixel accuracy from Landsat TM image data[C]//Proceedings of 25th International Symposium,Remote Sensing and Global Environmental Change.Austria,Graz:[s.n.],1993,2:155-161.

[12] Tatem A J,Lewis H G,Atkinson P M,et al.Super-resolution target identification from remotely sensed images using a Hopfield neural network[J].IEEE Transaction on Geoscience and Remote Sensing,2001,39(4):781-796.

[13] Nguyen M Q,Atkinson P M,Lewis H G.Superresolution mapping using a Hopfield neural network with fused images[J].IEEE Transactions on Geoscience and Remote Sensing,2006,44(3):736-749.

[14] Ruiz C P,López F J A.Restoring SPOT images using PSF-derived deconvolution filters[J].International Journal of Remote Sensing,2002,23(12):2379-2391.

[15] Verhoeye J.De Wulf R.Land cover mapping at sub-pixel scales using linear optimization techniques[J].Remote Sensing of Environment,2002,79(1):96-104.

[16] Mertens K C,Verbeke L P C,Westra T,et al.Sub-pixel mapping and sub-pixel sharpening using neural network predicted wavelet coefficients[J].Remote Sensing of Environment,2004,91(2):225-236.

[17] Mertens K C,De Baets B,Verbeke L P C,et al.A sub-pixel mapping algorithm based on sub-pixel/pixel spatial attraction models[J].International Journal of Remote Sensing,2006,27(15):3293-3310.

[18] Shen Z Q,Qi J G,Wang K.Modification of pixel-swapping algorithm with initialization from a sub-pixel/pixel spatial attraction model[J].Photogrammetric Engineering and Remote Sensing,2009,75(5):557-567.

[19] 沈掌泉,虞舟魯.基于單親遺傳算法和子像元/像元空間吸引模型的亞像元制圖研究[J].遙感技術與應用,2015,30(1):129-134.

Shen Z Q,Yu Z L.Sub-pixel mapping with partheno-genetic algorithm and sub-pixel/pixel spatial attraction model[J].Remote Sensing Technology and Application,2015,30(1):129-134.

[20] Muslim A M,Foody G M,Atkinson P M.Localized soft classification for super-resolution mapping of the shoreline[J].International Journal of Remote Sensing,2006,27(11):2271-2285.

[21] Thornton M W,Atkinson P M,Holland D A.Sub-pixel mapping of rural land cover objects from fine spatial resolution satellite sensor imagery using super-resolution pixel-swapping[J].International Journal of Remote Sensing,2006,27(3):473-491.

[22] Nguyen M Q,Atkinson P M,Lewis H G.Superresolution mapping using a Hopfield neural network with LiDAR data[J].IEEE Geoscience and Remote Sensing Letters,2005,2(3):366-370.

Sub-pixelmappingoflandcoverusingsub-pixelswappingalgorithmandtopographicdata

YU Zhoulu1， WANG Wenchao2， RONG Yi2， SHEN Zhangquan1

(1.CollegeofEnvironmentalandResourceSciences，ZhejiangUniversity，Hangzhou310058，China；2.HangzhouFederationofLandPlanningandDesignConsultingCo.Ltd.，Hangzhou310030，China)

When remote sensing images are used to provide information for land cover mapping, it is negatively affected by the occurrence of mixed pixels in the remote sensing images, particularly in the case of coarse spatial resolutions. Soft classification and super-resolution(sub-pixel) mapping techniques can solve this kind of problems. The pixel-swapping(PS)algorithm is a simple and efficient technique for sub-pixel mapping. However, its computation is inefficient and yields poor mapping accuracy when the super-resolution scale factor (S)is large. A possible reason for this is that it only relies on the information from the fraction images. In this study, the digital elevation model(DEM) and their derivative data are employed as supplementary information for the PS algorithm so as to improve super-resolution mapping accuracy. Some conclusions have been reached: ① The sub-pixel mapping accuracy could be improved with the assistance of the DEM even if the scale factor is large； ② The mapping accuracy by incorporating both elevation and slope information is better than that of using elevation or slope data alone; ③ Mapping accuracy is less sensitive to the number of neighbors when scale factor is large；④ The computing efficiency is improved when incorporating DEM in pixel-swapping. Thus, it is feasible to use DEM as supplemental information for sub-pixel mapping.

sub-pixel mapping； super-resolution mapping； sub- pixel swapping algorithm； remote sensing； digital elevation model(DEM)； land cover mapping

10.6046/gtzyyg.2017.04.14

虞舟魯,王文超,戎奕,等.基于子像元交換算法和地形數據的土地覆蓋亞像元制圖研究[J].國土資源遙感,2017,29(4)：88-97.(Yu Z L,Wang W C,Rong Y,et al.Sub-pixel mapping of land cover using sub-pixel swapping algorithm and topographic data[J].Remote Sensing for Land and Resources,2017,29(4)：88-97.)

TP 751.1

A

1001-070X(2017)04-0088-10

2016-05-05；

2016-06-28

國家科技支撐計劃項目“旱區多遙感平臺農田信息精準獲取技術集成與服務”(編號： 2012BAH29B04)資助。

虞舟魯(1986-),男,研究實習員,主要從事農業遙感與信息技術應用、土地利用等方面的研究。Email： yuzl@zju.edu.cn。

沈掌泉(1969-),男,副教授,主要從事農業遙感與信息技術、計算機應用等方面的研究。Email： zhqshen@zju.edu.cn。

(責任編輯：張仙)