武漢市1991~2013年城市擴張與地表覆蓋變化研究

孫安昌，陳 濤，張敘葭，牛瑞卿

（1.中國地質大學（武漢）地球物理與空間信息學院，430074）

武漢市1991~2013年城市擴張與地表覆蓋變化研究

孫安昌1，陳 濤1，張敘葭1，牛瑞卿1

（1.中國地質大學（武漢）地球物理與空間信息學院，430074）

利用1991～2013年間6個時相的武漢市Landsat遙感數據，使用多端元光譜混合分析方法進行混合像元分解，得到各時相的地表覆蓋豐度圖。根據豐度信息，計算土地變化動態度與土地覆蓋變化強度，得到地表覆蓋變化情況，分析武漢地區城市擴張與土地轉移特征。結果表明，武漢市主城區建設用地源自對城市外圍的綠地與內部湖泊地區的侵占，湖泊與綠地面積總體呈下降趨勢。

城市擴張；地表覆蓋變化；多端元光譜混合分析；武漢市；遙感

城市化過程將自然景觀轉化為建設用地，導致城市地表覆蓋物理特性變化迅速，造成城市內部環境發生改變[1]。例如，建設用地的快速增長會導致城市熱島效應加劇[2-3]。衛星遙感技術的發展為周期性大范圍城市變化監測提供了支持[4]。對于中、低空間分辨率的衛星影像，因城市地表復雜，場景中的部分地物小于影像像元的尺寸，不能被直接檢測[5]，傳感器獲取的不是單一地物的地表反射率，而是不同地物按照一定比例加權的總和，形成混合像元[6]。對于混合像元，使用傳統的圖像硬分類，直接將像元判為某一地物，必然會影響地物識別與圖像分類精度。

光譜混合分析（spectral mixture analysis，SMA）可將像元分解為不同的端元組分，求得其所占像元的百分含量，解決了傳統分類法對混合像元分析精度下降以及信息丟失的問題，得到更加準確的結果[4]，可用于解決混合像元在土地分類中遇到的問題。SMALL采用線性光譜混合模型研究城市植被空間與時間變化[7-8]；LU等人運用SMA技術提取城市不透水面，分類精度達到83．78%[9]。ROBERTS等人根據SMA提出了多端元光譜分解（multiple endmember spectral mixture analysis, MESMA），該方法通過調整端元數量，評估和選擇最優端元組合，以動態方式處理混合像元，被證明是一種更合理、精度更高的方法[6,10-13]。ROBERTS等人運用MESMA分別對美國Los Angeles、巴西Manaus的TM圖像進行城市土地覆蓋像元分解并測定城市形態物質組成，結果表明MESMA在亞像元水平定量獲取地表信息是有效的[10,13]。

MESMA在分析城市地表成分上應用廣泛，但少有使用該方法對基于時間序列的遙感影像進行定性定量的地表覆蓋變化研究。本文以武漢市為例，使用6個時相的Landsat TM遙感影像，運用MESMA進行混合像元分解，通過分解結果獲取地表信息，以此驗證該方法在多時相城市地表定量變化監測中的應用；同時分析總結武漢市近20 a來城市擴張的特點，為城市發展與規劃提供參考。

1 研究區概況

武漢市位于江漢平原東部，長江和漢水交匯處（113°41'E～115°05'E，29°50'N～31°22'N），市區總面積為8 494．41 km2，是我國中部地區第一大城市，也是重要的金融、商業、貿易、文化中心。武漢屬于亞熱帶季風氣候，雨量充沛、四季分明；在地形上北高南低，屬于殘丘性河湖沖積平原，北部為山地丘陵，其余均屬江漢平原，地勢平坦、河道縱橫、湖泊星羅棋布。

2 數據獲取與預處理

2.1 數據準備

本文使用的遙感數據為Landsat影像，共6景（表1）；另外還有全國行政省、市（縣）邊界圖及其他相關資料，所用圖像處理軟件為Envi5．0。

表1 遙感數據介紹

2.2 數據預處理

本文使用的Landsat數據為L1T級，已經過輻射校正與幾何校正。在6景影像內均勻選取同名像點進行配準，誤差控制在0．5個像素內，并用矢量數據進行裁剪。考慮到數據限制以及研究需求，僅選取武漢的主城區作為研究區域，即江岸區、江漢區、硚口區、漢陽區、武昌區、青山區、洪山區。

圖1 研究區遙感影像縮略圖（Landsat8 R、Band4 G、Band3 B、Band2）

3 算法及模型

3.1 線性光譜混合模型

光譜混合分析是將混合像元表現為地表純凈端元光譜組合形式[4]。其常用線性光譜模型定義為：像元在某一光譜波段的反射率是由構成像元基本組分的反射率及其所占像元面積比例為權重系數的線性組合。可用以下公式表達：

其中，Riλ為第λ波段第i像元的反射率（已知）；N為端元數目；fi為對應第i像元某一端元的豐度；ρiλ為各端元的反射率；ελ為殘余誤差值。

評價模型用殘差ελ或均方根誤差RMSE表示：

其中，M是傳感器波段數。

線性模型從混合像元中分離和提取各組分的平均光譜響應，通過求解線性方程來反解端元在像元中所占的比例（豐度），從而將所有像元分解成這些基本組分的分量。模型計算的結果表現為各端元的分量值和以均方根誤差表示的殘余誤差圖像。

3.2 端元提取

端元必須準確全面地代表地面材料（地表覆蓋類型）[15]。建立應用于混合像元分解的端元光譜庫的挑戰是雙面的。首先，每類地物的光譜庫應該包含足夠數量的光譜，這樣才能充分代表地面材料的光譜變化。然而，隨著端元總數的增加，又會導致計算效率下降[14]。為達到以上目的，更好地選擇出能夠代表土地覆蓋類別的端元，本文參考RIDD提出的概念[17]，結合研究區域的實際情況，定義端元類別為水體、植被、建設用地、裸地。具體說明如下：

1）水體包括江河湖泊、水田、漁場。

2）植被包括林地、草場、農田。

3）建設用地包括所有建成的人工建筑、道路等。4）裸地包括土壤、未建成的建設用地。

3.3 最優端元光譜篩選

將端元集制作成光譜庫，并通過參數EAR（endmember average RMSE）從庫中篩選出最優端元光譜。EAR由DENNISON等在2003年提出[15]，該參數能夠評判出地表覆蓋類型中最具有代表性的端元。計算公式如下：

式中，A代表端元類別；n為類別A的光譜數量；Ai為具體的光譜。

EAR是用來評價每個端元模擬分解其地表類型能力的參數，用于衡量端元在該類別分解過程中的表現能力，不能用來衡量端元的純凈程度。EAR值較低的端元比EAR值較高的端元在分解過程中有更好的表現能力，即EAR值較低的端元是最能代表這一地物類別的端元。因此，篩選的最優端元應具有較低的EAR值。考慮到數據空間分辨率和光譜分辨率的限制，以及研究中對研究區域影像的判讀和了解，不同時期數據中同一地表覆蓋類別所篩選的光譜數不相同。各端元光譜數量見表2。

表2 不同時相地表覆蓋模型數與像元成功分解比例

3.4 混合像元分解

本文根據線性光譜混合模型，采用雙端元模型，將像元分解為2個組分。增加一個端元個數雖然可以減少RMSE，但對結果的提升不大[15]。過多端元會造成端元之間光譜的混淆，隨著端元個數增加，模型數量也隨之增加，計算效率就會下降。雙端元即代表實際地物的端元光譜與考慮到傳感器照明的變化而加上的陰影（shade，SHD）端元。本次研究中，設置分解的RMSE小于 0．025，若RMSE大于該閾值，則劃為未識別區域。

3.5 地表信息參數

引入兩個地表信息參數：土地變化動態度和土地覆蓋變化強度。

1）土地變化動態度。土地變化動態度公式如下：

式中，AV表示某土地利用類型的變化率；Aε和As分別表示前后兩個時相的面積，該面積為像元的豐度乘以該像元所代表的面積，T表示兩個時相的時間差（單位：年）。該系數能定量描述土地利用變化狀況，土地變化動態度越大說明該地表覆蓋類別在相應時間段內的變化越劇烈。

2）土地覆蓋變化強度。參考MICHISHITA等人[16]的做法，引入土地覆蓋變化強度（LCCI）。土地覆蓋變化以像素為單位，反映了兩景影像時間內每個土地類別每天的變化情況。公式如下：

式中，LCFt1和LCFt2為兩個時相的LCF圖；?t為兩個時相間隔的天數；A為像元所代表的面積，即900 m2，LCCI的單位為m2/day。

4 結果與討論

4.1 精度評價

分解成功像元所占影像的比例如表2所示。1991～2009年期間5個時相的數據分解成功的像元所占比例大于99%，2013年成功率為94．14%。圖2為1991年分解結果。

由于缺乏地表覆蓋豐度數據，無法從豐度上給出精度評價。因此，本文將每個地表覆蓋類別在影像上隨機選擇500個像元，計算得到1991～2013年混淆矩陣。結果表明，總體精度最高達到90%，最低為75%。這4個地表覆蓋類型中，裸地的混淆情況比較嚴重，在該光譜分辨率下裸地與某些建設用地的光譜相近。但總體而言，結果滿足進一步研究需求。

圖2 分解結果示例（1991）

4.2 地表覆蓋變化分析

根據土地變化動態度的計算結果（表3），從1991～2013年，研究區域內建設用地面積增長迅速，年均增長率為4．36%，植被、水體以及裸地呈負增長，變化率分別為-2．33%、-1．1%、-3．45%。

建設用地在1996～2000年變化最快，達到6．82%，隨后增長速度放緩，2005～2009年在研究區域出現了負增長，減少了1．69%。建設用地的整體變化率是所有土地類型中最大的，其快速增長也反映了城市20多年來的迅猛發展。

在綠地面積減少的大趨勢下，2000～2005年變化率減小，僅為-0．74%；在2005～2009年，則有小幅度的上升。對比綠地與建設用地，變化相反且強度相對一致。

水體面積也呈下降趨勢，但在1996～2000年有明顯的恢復。經查閱資料，武漢市為了經濟發展，采取了“退耕還漁”的相關策略，養殖水面增加，導致湖泊也有較大增長，這與曹麗琴[2]的研究結果一致。至2009年，水體面積與1996年基本持平，2009～2013年則降到更低。

裸地變化雜亂，其整體呈減少趨勢，但在1991～1996年和2005～2009年間又有明顯增長。

表3 1991～2013年武漢主城區土地變化動態度/%

4.3 城區擴張與土地轉移分析

為了探究城市擴張以及土地轉移的情況，本文提取研究區域2013年建設用地部分，并對其他影像進行掩膜處理，制成時間序列的土地覆蓋豐度變化圖（圖3）。其中，紅色表示建設用地，綠色、藍色和灰色分別表示綠地、水體和裸地在2013年前轉換為建設用地的部分；地表覆蓋的豐度越高，對應的顏色越深。

圖3 武漢主城區土地覆蓋豐度變化圖

該時間序列圖主要闡述了武漢地區建設用地的增長模式以及土地來源。從空間上看，建設用地的擴張是以武漢的地理中心逐層向外，以圈層式的結構發展。漢口、漢陽和武昌的擴張特點各不相同（圖3a）。漢口地區依托老城區，呈弧線向西北地區逐層發展；漢陽地區的變化緩慢，主要沿著長江與漢江方向擴張；武昌地區變化較大，在2009年以前，一方面向東擴張，另一方面青山工業地區向南和武昌建成的老城區相連，2009～2013年向外圍郊區擴張明顯。

建設用地擴張侵占了植被以及水體面積。對于圖中藍色被侵占的水體，主要為城市內部湖泊，且存在于2005年以前，尤其是2000～2005年（圖3c、3d），大量藍色區域消失，填湖情況嚴重；但在2005年以后，建設用地對水體的侵占不明顯。

4.4 城市擴張模式分析

根據土地覆蓋變化強度，將建設用地的發展分為兩種形式：一種是新建設用地的擴張；另一種是城區內的改造或再發展。新建設用地的開發是其他土地變化為建設用地的覆蓋類型，而再發展是將已有的建設用地進行改造，造成其豐度變化。

由表4可以看出，建設用地變化最快的時間段為1996～2000年，植被、水體、裸地快速減少，其中植被和水體的變化強度在歷年中最大；1991～1996年和2009～2013年建設用地增長速度適中，相應的其他土地覆蓋變化強度為負，這3個時間段主要為新城區的擴張。2000～2005年建設用地的變化強度僅次于1996～2000年，但植被和水體變化強度弱，即該段時間武漢城市的變化主要為城區內的改造和再發展。在2005～2009年，建設用地的面積出現了負增長，本文認為這也與城市內部的改造與再發展有關系。

表4 土地覆蓋變化強度表/m2/day

5 結 語

本文使用6個時相的武漢Landsat遙感影像為數據源，使用多端元光譜混合分析方法進行雙端元模型的混合像元分解，得到結論如下：

1）通過目視解譯直接選取端元，使用EAR參數進行端元光譜篩選，并采用雙端元模型進行多端元光譜混合分析，成功提取了地表覆蓋類別與豐度信息，精度能夠滿足研究需求。根據結果計算土地利用動態度與LCCI，繪制時間序列LCF圖與LCCI圖，將土地覆蓋類別、土地覆蓋豐度與空間信息相結合，較好地監測并分析了1991～2013年武漢城市擴張情況與特點。

2）1991～2013年這23 a中，武漢市主城區地表變化劇烈，建設用地城市區域擴張明顯，植被和水體總體減少，裸地則隨年份不同變化較大。其中，建設用地和植被區域在2005～2009年分別減少和增加，水體在1996～2000年間也有增長，這3個變化與整體趨勢相背。

3）武漢的城市擴張以長江與漢江的交匯處為中心，逐層向外以圈層式的結構發展。在1991～2000年和2009～2013年主要為城市向外擴張，在2000～2009年更多表現為內部重建。漢口、漢陽和武昌三鎮擴張特點不同，漢口依托建成的老城區，呈弧線向西北地區逐層發展；漢陽地區的變化緩慢，主要沿著長江與漢江擴張；武昌地區和青山地區逐漸發展相連，隨后沿東軸呈跳躍式增長。參考文獻

[1] HUANG Wei, ZENG Yongnian,LI Songnian． An Ananlysis of Urban Expansion and Its Associated Thermal Characteristics Using Landsat Imagery [J]． Geocarto International, 2015(30):93-103

[2] 曹麗琴．城市熱環境及其影像因子的定量遙感研究[D]．武漢:武漢大學,2009

[3] 陳濤,孫安昌,王鑫,等．基于遙感的武漢地區LUCC對熱島效應的影響研究[J]．地理空間信息,2015,13(1):62-65

[4] 趙英時．遙感應用分析原理與方法[M]．北京:科學出版社,2003

[5] ALAN H S, CURTIS E W, JAMES A S． On the Nature of Models in Remote Sensing [J]．Remote Sensing of Environment, 1986, 20(2): 121-139

[6] TAREK R, JOHN R W, DAR R, et al． Measuring the Physical Composition of Urban． Morphology Using Multiple Endmember Spectral Mixture Models[J]．Photogrammetric Engineering and Remote Sensing, 2003(69):1 011-1 020

[7] SMALL C． Multitemporal Analysis of Urban Reflectance [J]．Remote Sensing of Environment, 2002(81):427-442

[8] SMALL C． Estimation of Urban Vegetation Abundance by Spectral Mixture Analysis [J]． International Journal of Remote Sensing, 2001, 22(7):1 305-1 334

[9] LU Dengsheng, WENG Qihao． Urban Classification Using Full Spectral Information of Landsat ETM+Imagery in Marion County, Indiana [J]． Photogrammetric Engineering & Remote Sensing, 2005,71(11):1 275-1 284

[10] POWDLL R L,ROBERTS D A,DENNISON P E, et al． Subpixel Mapping of Urban Land Cover Using Multiple Endmember Spectral Mixture Analysis: Manaus, Brazil [J]． Remote Sensing of Environment, 2007(106):253-267

[11] WU C． Normalized Spectral Mixture Analysis for Monitoring Urban Composition Using ETM+Imagery [J] ．Remote Sensing of Environment, 2004(93): 480-492

[12] FRANKE J, ROBERTS D A, HALLIGAN K, et al．Hierarchical Multiple Endmembers Spectral Mixture Analysis (MESMA) of Hyperspectral Imagery for Urban Environments [J]．Remote Sensing of Environment,2009(113):1 712-1 723

[13] ROBERTS D A, GARDNER M, CHURCH R, et al．Mapping Chaparral in the Santa Monica Mountains Using Multiple Endmember Spectral Mixture Models [J]． Remote Sensing of Environment, 1998, 65(3): 267-279

[14] ROBERTS D A, NUMATA I, HOLMES K, et al． Large Area Mapping of Land-cover Change in Rondonia Using Multitemporal Spectral Mixture Analysis and Decision Tree Classifiers [J]． Journal of Geophysical Research Atmospheres, 2002, 107(D20): LBA 40-1-LBA 40-18

[15] DENNISON P E, ROBERTS D A． Endmember Selection for Multiple Endmember Spectral Mixture Analysis Using Endmember Average RMSE [J]． Remote Sensing of Environment, 2003(87):123-135

[16] RYO M, JIANG Zhiben, XU Bing． Monitoring Two Decades of Urbanization in the Poyang Lake Area, China Through Spectral Unmixing [J]．Remote Sensing of Environment, 2012(117):3-18

[17] RIDD M K． Exploring a V-I-S (vegetation-Impervious Surface-soil) Model for Urban Ecosystem Analysis Through Remote Sensing: Comparative Anatomy for Cities [J]．International Journal of Remote Sensing, 1995(16): 2 165-2 185

P273

B

1672-4623（2016）12-0033-04

10.3969/j.issn.1672-4623.2016.12.012

孫安昌，碩士研究生，研究方向為環境遙感、遙感地質、高光譜遙感。

2015-08-10。

項目來源：國家高技術研究發展計劃資助項目（2012AA121303）；數字制圖與國土信息應用工程國家測繪地理信息局重點實驗室開放研究基金資助項目（GCWD201202）；湖北省自然科學基金資助項目（2012FFB06501）；長江科學院開放研究基金資助項目 （CKWV2013221/KY）。