武汉市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）。