基于遥感的城市典型地物信息的快速获取方法

李世伟，王召巴，杨建生

（1.中北大学 信息与通信工程学院，山西 太原030051；2.美国鲍尔州立大学 地理系，曼西 印第安那47036；3.中北大学 山西省光电信息与仪器工程技术研究中心，山西 太原030051）

0 引 言

近年来，我国城市进程加快，地物分布信息变化频繁，大量人工构建面成为了城市的主体，造成诸如城市热岛效应、城市内涝等问题时有发生。作为人类活动频繁的城市区域，其典型地物主要包括植被、水体、不透水面（包括道路、房屋等）和裸土区，对这些地物分布信息的及时掌握，是有效监测城市生态环境和制定规划决策的关键。随着科学技术的不断发展，基于遥感的城市地物信息获取方法得被广泛应用［1－4］。

目前，被广泛认可并能得到较好提取结果的方法是利用地物光谱特征进行组合运算的指数法。每一个指数都可以针对一类地物信息进行专题提取，通过波段运算和阈值划分，可以快速获得专题信息。这些指数中，比较知名的包括归一化植被指数（NDVI）［5］及其改进型指数［6］、归一化水体指数（NDWI）和改进型归一化水体指数（MNDWI）［7］、归一化不透水层指数（NDISI）［8，9］。但值得注意的是，这些指数运算结果只能对一种专题信息进行显示，无法使用户直接获得多类地物的分类信息，因此，本文就这一问题展开研究。

1 研究区简介和数据使用

为了信息收集和对实验结果进行验证的方便，本文选取了作者所在城市太原市。太原主城区位于太原河谷平原北部，属盆地结构，地处太行山与吕梁山之间，属暖温带大陆性季风气候，季节性明显。冬季通常寒冷干燥，到了春季气温迅速回升，雨量集中在夏季，昼夜温差较大；秋季天高气爽，气温变化起伏较小，全年光热资源充足，年均温度介于7.6－10.2℃之间，年均降水量为420－460毫米，年均降水量满足作物的基本生长需求，但降水量的季节差异较大，易造成干旱灾害和水土流失。主城区自北向南，汾河水从城区中部流过，在城区汾河中段的西面，分布着城区内唯一的湖泊－晋阳湖。

近年来，太原市经历了城市扩张后再绿化的过程，土地覆盖类型变化剧烈。为了对预想的方法进行验证，选取了一景LANDSAT－TM 影像，该影像的成像时间为2011年4月19 日（轨道：PATH－125，ROW－33），影像中心经度为112.311度，纬度36.046度。在这个月份，研究区的耕地作物还没有长出地面，因此对应区域全部为土壤覆盖，而除去农作物外的其它植被诸如树木和草地等已经开始生长，呈现出了较为明显的植被特征。此外，为了对最后的分类结果精度进行验证，选取了相同时相的一景SPOT 高分辨率影像作为参考。

为了利用指数运算的快速性优势，并同时获得多类地物信息，按照以下步骤进行实验：首先，使用指数法对专题信息进行提取；而后，针对每个指数的提取图像结果进行阈值划分和二值化处理；最后，在不破坏专题信息特征的基础上，利用伪彩色赋色的方法直接对专题信息的二值化图层进行图层堆叠和伪彩色赋色，以颜色作为分类依据，直接呈现多类地物分布信息。即：①使用指数分别对水体、植被和不透水层专题信息进行提取；②利用高分辨率影像，通过目视解译和对比对信息模糊地点进行辨别，得到每个专题信息图层的合理阈值，并使用二值化方法得到专题信息；③对得到的图层进行伪彩色合成实验，得到可以反映四类地物信息的最好色彩组合方式；④以高分辨率影像为参考，得到仅通过颜色信息获得的地物分类精度结果，以验证该方法的可行性。

2 数据预处理

2.1 数据预处理

所选TM 影像和SPOT 影像均为二级产品，已经过系统级辐射校正和几何校正。为了获得更为精确的结果，首先对影像进行了几何精校正，在采用多项式模型和双线性插值法进行了几何校正和图像重采样后，从而得到了精校正影像。其次，使用了1：25万山西省行政区边界矢量数据生成了太原市主城区的矢量边界，利用此边界生成研究区AOI数据层对影像进行裁剪。最终获得的太原市主城区影像如图1右图所示。

图1 太原市主城区遥感影像

2.2 反射率反演

由于遥感图像在成像时，会受到成像角度、太阳高度角等因素的影响，所获取的图像数据均以灰度值进行显示，会造成 “同物异谱”和 “同谱异物”现象的发生，为了在一定程度上消除这些影响，最为有效的办法就是将灰度图像反演为反射率影像。为此，首先将研究区的TM 影像的第1－5波段和第7波段转换为反射率影像。其过程如下

（1）将6个波段的图像DN 值转换为辐射亮度值

式中：R——像元的辐射亮度值，DN 原始影像中的像元数字灰度，Gain、Bias——USGSLANDSAT－TM 用户手册中提供的对应波段的增益数据和偏置数据（见表1）。

表1 LANDSAT－TM 增益和偏置数据

（2）反射率计算

式中：ρ——像元的反射率，D——日地天文单位距离（通常取1），R——式（1）中的像元辐射亮度值，ESUNI——大气顶层太阳辐射平均值（见表2），——太阳天顶角。计算天顶角的计算公式为

式中：SUN＿ELEVATION——对应影像头文件中的太阳高度角，在所使用的影像头文件中，其值为57°。

表2 大气顶层太阳辐射平均值

根据以上的计算过程，利用ERDAS IMAGINE 遥感图像处理软件的MODLE MAKER 建立了相应的计算模型并进行计算，得到了反射率反演结果影像，结果如图2所示（RGB＝432）。

3 专题信息提取

为了在将来使用本方法时能够提高工作效率，利用ERDAS IMAGINE遥感图像处理软件的MODLE MAKER建立了相应的指数计算模型，从而在对其它区域或本区域其它时相使用时，快速生成对应的专题信息图层。

3.1 植被信息的获取

通过归一化植被指数NDVI进行植被信息的获取，其表达式为

式中：ρNIR——第4 波段的反射率，ρRed——第三波段的反射率，其结果如图3所示。其中灰度值相对较高的区域为植被区域，灰度值较低的为非植被区域。为了使结果更具有可读性，在结合高分辨率影像进行了多次植被信息提取的阈值选取后，对NDVI进行了二值化处理，处理阈值为0.374，大于阈值区域显示为白色，反之为黑色，二值化处理图像如图4所示。

图3 NDVI处理结果

3.2 水体信息的提取

水体信息的提取采用了由徐涵秋教授于2005年提出的改进型归一化水体指数（MNDWI），其表达式为

式中：Green——绿光波段，MIR——中红外波段，分别为LANDSAT－TM 的第2 波段和第5 波段，提取结果如图5所示，其中灰度值较高处为水体。同样为了使结果更具有可读性，在结合高分辨率影像进行了多次植被信息提取的阈值选取后，对MNDWI进行了二值化处理，处理阈值为0.197，大于阈值区域显示为白色，反之为黑色，二值化处理图像如图6所示。

图4 NDVI处理结果

图5 MNDWI阈值处理结果

图6 MNDWI阈值处理结果

3.3 不透水层信息的提取

不透水层信息的提取，目前已经有多种方法。这些方法中，大部分都是基于Ridd于1995 年提出的VIS 模型。但是这些模型通常不能体现水体信息，因此，在不透水层提取前需要先对水体信息进行掩膜操作，将水体信息提出后在进行下一步的操作，这种做法势必会将误差引入不透水层提取的工作中来［10］。2008年徐涵秋教授在总结和比较已有指数的基础上，重新通过对地物各自在影像上的光谱特征以及地物之间的光谱特征差异分析，找到了能够使地物差异最大化的波段进行组合计算，提出了归一化不透水层指数（NDISI）来对不透水层信息进行提取。其表达式如下

式中：TIR——热红外波段，对TM 数据就是第六波段。NIR、MIR1——TM第4和第5波段。本文首先使用了NDISI指数对影像进行处理，但经过对提取结果进行解译和实地验证，我们在NDISI提取结果中发现：太原武宿机场等新建不透水面的灰度值高于水体，而大量道路、楼宇等建筑的灰度出现了低于水体灰度的情况，甚至出现少量不透水层与水体灰度值相等的情况。

因此，即使不考虑少量不透水层与水体灰度值相等的情况，单纯使用NDISI提取不透水层也需要两个以上的阈值才可以将不透水层信息提取出来，并且除了每增加一个阈值会相应增加误差外，后续还需要对处理结果进行掩膜、合并等大量操作。基于这种考虑，我们采取了另外的方法进行了不透水面的提取。考虑到指数方法基于地物光谱特征，我们重新对所使用影像中的典型地物进行光谱特征采集（如图7所示），图中纵坐标为TM 影像波段编号，纵坐标为多个典型地物样点的灰度均值。通过对该特征曲线组进行分析后，得到了以下结果：

图7 四类地物光谱特征

不透水层在第3波段和第4波段的差异均大于其它3种地物包括水体、植被和土壤，若将4波段和3波段相除，可达到增大地物灰度差异的目的；不透水层在5 波段和6波段的光谱特征与其它3种地物呈现出相反的情况，即不透水层第5波段灰度大于其第6波段，而其它3种地物变现为第6波段灰度大于其地物波段；在第4波段和第5波段，水体也呈现出于其它3种地物的相反特征。

基于以上结果，并考虑到第6波段与其它6个波段分辨率不一致，且MNDWI的提取结果中水体呈现了较好的提取效果，因此采取了以下方案对不透水面进行提取：

（1）首先对影像的第4和第3波段进行比值操作（RI，式（7）），对原图像的不透水层和水体信息同时进行增强，结果如图8所示。

通过式（7）计算后，不透水层信息和水体信息同时得到增强，且水体灰度略高于高于不透水层，并且所用不透水层的灰度值大小相当，没有出现NDISI中出现的不透水层灰度既有大于水体的部分又有小于水体灰度的情况的发生。

图8 RI结果影像

（2）由于在RI图层中，水体和不透水层均具有较高的灰度，在MNDWI中的水体信息得到增强而不透水层信息被抑制，为此，首先将RI图层在进行阈值选取后的二值化处理，处理阈值为1.159，而后，将该二值化图层与MNDWI二值化图层进行差值运算后，最终得到了不透水层的信息图层（如图9所示）。

4 伪彩色合成分类和精度评价

为了实现分类，对所获取的专题图层进行分析，考虑到对NDVI、MNDWI、RISI减MNDWI的3个提取结果已经进行了二值化处理，又考虑到ERDAS软件可以快速实现对3个图层进行堆叠和伪彩色赋色的操作，因此直接对3个二值化图层进行了图层堆叠，并分别赋予不同的伪彩色，这样做的原因是：

在植被专题信息、水体专题信息和不透水层专题信息层的二值化图像中，其各自的专题信息区域灰度值为1，而其它信息区域灰度为0，因此，若将3个专题信息图层赋予不同的伪彩色并合成后，各专题信息在合成结果中的颜色应具有较大差异。而且，由于土壤在3个专题图层上的灰度值都为零，无论对哪一图层赋予除黑色外的伪彩色后，土壤都将呈现为黑色。因此，在使用堆叠三图层并伪彩色赋色后，四类地物应该可以从颜色进行区分。

图9 不透水层提取结果

基于以上考虑，在将3 个图层进行图层堆叠操作后，对NDVI二值化图层赋为绿色，不透水面二值化图层赋为红色，水体二值化图层赋为蓝色（如图10所示）。

为了对本文得到的结果进行验证，得到该方法的可行程度，在最后得到的伪彩色合成结果中随机选取了200个点，在结合高分辨率影像进行对比后，获得的具体分类精度见表3。

表3 分类精度评价

可以看出，该方法可得到83%的总体分类精度，为了使评价结果更具有客观性，通过式 （8）对分类结果进行了KAPPA 系数的计算。

式中：N——参与精度统计的像元数量，Xii——分类正确的像元数量Xi＋是表格中第i行像元的数量，X＋i——表格中第i列像元的数量，所得KAPPA 系数结果为0.803。

图10 伪彩色合成结果

通过参考同时期的高分辨率影像，本文使用的方法所造成的误分情况多发生在汾河南端的河道和沿岸发生，通过实地考察和询问，并结合相关资料分析发现，在影像获取时段内，汾河太原段中断和南段水流较小，河道裸露，小型死水面较多，其面积小于影像中一个像元所对应的面积，造成了遥感影像中混合像元的出现，因此，在一定程度上影响了分类精度。

5 结束语

采用建模方式建立了多个指数模型，并通过运算得到了相应的专题信息，采用阈值划分和二值化处理的方式后，在没有使用任何分类算法的情况下，利用伪彩色合成方法同时得到了包括水体、植被、不透水面和裸土共四类地物的分布信息。其中，利用NDVI和MNDWI分别对植被和水体进行了提取，结合不透水层和水体两者在MNDWI图层以及经过光谱分析计算建立的RI图层中的特点，对二值化后的RI图层和MNDWI图层进行差值运算后得到了不透水层信息。在将得到的分类结果图像与高分辨率参考图像对比后，得到了分类结果的总体分类精度为83%，Kappa系数0.803。结果表明，该方法通过色彩区分就可直接快速得到研究区四类地物的分类结果和分布影像，省去了常规分类方法中必需的类别合并、分类赋色等过程，模型实现容易，阈值改变仅需在模型中改变数值即可，实现了四类城市典型地物分布信息的快速获取，可为地物分布信息的快速监测提供参考技术手段。

［1］ZHAO Ming.Analysis on the evolvement of Wenzhou City’s morphological elements ［J］.Urban Studies，2009，16 （7）：28－32 （in Chinese）.［赵明.温州城市形态构成要素演化研究［J］.城市发展研究，2009，16 （7）：28－32.］

［2］SONG Yang，WANG Shijun，YANG Yanru，et al.Study on the urban internal spatial structure in northeast China ［J］.Areal Research and Development，2009，28 （4）：18－23 （in Chinese）.［宋飏，王士君，杨艳茹，等.东北地区城市内部空间结构研究 ［J］.地域研究与开发，2009，28 （4）：18－23.］

［3］MU Fengyun，ZHANG Zengxiang，TAN Wenbin.Analysis on the spatial－temporal characteristics of Guangzhou City’s spatial morphologic evolution ［J］.Geo－Information Science，2007，9 （5）：94－98 （in Chinese）. ［牟凤云，张增祥，谭文彬.广州城市空间形态特征与时空演化分析 ［J］.地球信息科学，2007，9 （5）：94－98.］

［4］LIU Yaxuan，ZHANG Xiaolei，LEI Jun，et al.Spatial expansion and driving forces of oasis cities in Xinjiang，China ［J］.Journal of Desert Research，2011，31 （4）：1015－1021 （in Chinese）.［刘雅轩，张小雷，雷军，等.新疆绿洲城市空间扩展特征及其驱动力分析 ［J］.中国沙漠，2011，31 （4）：1015－1021.］

［5］SONG Dongmei，ZHANG Qian，YANG Xiuchun，et al.Spatial and temporal characteristics of MODIS vegetation index in the source region of three rivers on Qinghai－Tibet plateau in China［J］.Geographical Research，2011，30 （11）：2067－2075（in Chinese）. ［宋冬梅，张茜，杨秀春，等.三江源区MODIS植被指数时空分布特征 ［J］.地理研究，2011，30 （11）：2067－2075.］

［6］ZHU Gaolong，LIU Yibo，JU Weimin，et al.Evaluation of topographic effects on four commonly used vegetation indices 4［J］.Journal of Remote Sensing，2013，17 （1）：210－234 （in Chinese）.［朱高龙，柳艺博，居为民，等.4种常用植被指数的地形效应评估 ［J］.遥感学报，2013，17 （1）：210－234.］

［7］XU Hanqiu.Comment on the enhanced water index（EW I）：A D iscussion on the creation of awater index ［J］.Geo－Information Science，2008，10 （6）：776－780 （in Chinese）.［徐涵秋.从增强型水体指数分析遥感水体指数的创建 ［J］.地球信息科学，2008，10 （6）：776－780.］

［8］XU Hanqiu.A new remote sensing index for fastly extracting impervious surface information ［J］.Geomatics and Information Science of Wuhan University，2008，33 （11）：1150－1153 （in Chinese）.［徐涵秋.一种快速提取不透水面的新型遥感指数［J］.武汉大学学报·信息科学版，2008，33 （11）：1150－1153.］

［9］XU Hanqiu.Analysis of impervious surface and its impact on urban heat environment using the normalized difference impervious surface index（NDISI） ［J］.Photogrammetric Engineering ＆Remote Sensing，2010，76 （5）：557－565.

［10］MU Fengyun，ZHANG Zengxiang，TAN Wenbin.Analysis on the spatial－temporal characteristics of Guangzhou City’s spatial morphologic evolution ［J］.Geo－Information Science，2007，9 （5）：94－98 （in Chinese）.［牟凤云，张增祥，谭文彬.广州城市空间形态特征与时空演化分析 ［J］.地球信息科学，2007，9 （5）：94－98.］