基于Landsat8数据和“珞珈一号”夜光数据的合肥建成区提取

2021-08-17吴庆双张伯琦蔺陆洲

安徽师范大学学报(自然科学版) 2021年4期

李强，吴庆双，周华，张伯琦，李俊，蔺陆洲

(1.安徽师范大学地理与旅游学院，安徽芜湖 241003；2.全图通位置网络有限公司，北京 100716；3.武汉吉威空间信息技术研究院，湖北武汉 430000)

中国目前正处于城市化高速发展时期，建成区范围作为城市化的重要指标之一，直接反映一个城市的经济、社会、文化、科技发展水平，同时也对局地气候、生态环境等一系列自然因素产生显著的影响。因此，获取准确的建成区对城市的区域规划、管理、经济发展具有指导意义。合肥市作为安徽省省会，正处于城市化快速发展阶段，从建设国家综合性科学中心到融入长三角一体化高质量发展过程中，合肥市的城市建设有了长足的发展。但是，如果要进一步加快合肥城市化进程，应对城市化过程中存在的问题，实现合肥市“大湖名城，创新高地”的宏伟目标，了解合肥市建成区的范围具有重大意义。

当前，运用遥感技术进行建成区提取的方法，总体上可以分为以下三种。

第一种是根据高分辨卫星影像(主要是Landsat系列影像)的光谱信息和纹理特征进行图像分类提取建成区：温奇通过分析高分辨率遥感影像中平原建成区的纹理特征和局部关键点特征,提出了基于多核学习、多尺度分割以及多假设投票的平原建成区提取方法[1]；胡德勇等采用监督分类、非监督分类和归一化裸露指数(NDBI)等方法提取城市用地信息[2]；马生清在Landsat8遥感影像融合的基础上对建成区进行提取[3]。虽然在提取算法和波段处理过程中不断提高图像分类的精度，但是由于裸地及部分地物的光谱特征与建成区极为相似，在进行图像分类过程中，难免会被分类为建成区。

第二种是运用美国DMSP卫星的OLS传感器获取的夜光数据作为主要数据源，通过一系列不同的方法获取精确阈值，根据阈值进行二值化提取建成区：Imhoff提出的突变检测法，即随着灯光阈值的增加，在达到某个阈值时，高亮区域的边长会发现突变，把此时阈值作为最佳阈值，进行二值化，提取建成区[4]；苏泳娴等基于DMSP/OLS数据采用邻域分析法进行了建成区的提取[5]；王海军等基于DMSP/OLS数据采用图像识别、分割等方法对成渝地区城市化进行检测[6]，相较于第一种方法，由于建成区和裸地在夜光数据上的亮度值存在明显的差异，可较为高效地排除裸地等地物对于建成区提取的影响，但是DMSP/OLS夜光数据的空间分辨率只有1km，同时夜光数据存在外溢效果，使得提取的建成区有较大误差。

第三种方法则是结合夜光数据和高分辨的多源遥感数据进行建成区提取：Cao等采用支持向量机的区域生长算法，结合DMSP/OLS夜间灯光数和SPOT NDVI产品提取了中国东部2000年25个城镇的建成区[7]；柴宝惠等利用DMSP/OLS数据采用阈值法确定建成区的大致范围，再在该范围内对Landsat影像进行分类，得到天津市的建成区分布[8]；王若曦等对Landsat数据进行监督分类完成建成区初步提取，对DMSP/OLS影像进行形态学运算，运用得到的阈值提取建成区[9]；刘建良等以福州市为例对阈值法、突变检测、VAUNI指数法、基于辅助指数法的空间比较法等提取建成区的方法进行对比分析[10]。利用夜光数据去除裸地，再用高分辨率卫星提高精度，这种结合多种数据源的方法大大提升了提取建成区的精度。

本文以合肥市作为研究区域，以“珞珈一号”夜光数据、Landsat8数据为主要数据源，通过多源数据融合的方式提取建成区，主要有以下几个创新点：

(1)采用国产的“珞珈一号”夜光数据代替国外的DMSP/OLS和NPP/VIIRS夜光数据进行建成区提取，“珞珈一号”夜光数据相比前者具有巨大优势，大大提高了提取的精度；

(2)在对Landsat8影像采用“三指数法”压缩数据维，减少数据冗余的过程中，根据实际建成区的特征，采用了不同以往的三个指数。

1 研究区及数据

1.1 研究区

安徽省省会合肥市，位于安徽省正中部，长江淮河之间、巢湖之滨，具有承东启西、贯通南北的重要区位优势，是国家级皖江城市带承接产业转移示范区核心城市、长三角城市经济协调会城市[11]，同时也是安徽的政治、经济、文化、信息、交通、金融和商贸中心。近些年，合肥经济迅速发展，也从不同程度上刺激着建成区的不断扩张。图1为合肥市地理区位图。

图1 合肥市地理区位图Fig.1 Geographical location map of Hefei

1.2 “珞珈一号”夜光数据

世界上第一颗夜光卫星是美国国防部发射的DMSP/OLS卫星，可获得1992—2013年空间分辨率1km的年度夜光数据；2011年，在NPP/VIIRS卫星发射成功后，可以获得2013年后的空间分辨率500米的月度夜光数据。但是随着国内外对夜光数据研究的不断深入，获取更高空间分辨率和时间分辨率的夜光数据显得尤为重要，同时为了摆脱我国对夜光数据源过于依赖国外的局面，中国于2018年6月成功发射了我国第一颗专业夜光遥感卫星——“珞珈一号”。其数据的空间分辨率相较于前面两款夜光数据有了较大的提升，可达到130米，表明“珞珈一号”可以捕捉到更为细微的夜光空间信息。在时间分辨率上，“珞珈一号”理论上15天就可以获得全球的一期夜光数据[12]。是目前研究地表人类活动、社会经济参量估算、城市监测等社会热点问题的重要数据源。“珞珈一号”卫星参数见表1。

表1 “珞珈一号”卫星参数Table 1 Parameters of luojia-1 satellite

1.3 Landsat8卫星数据

1972年以来，美国NASA陆地卫星计划已经发射了8颗卫星，Landsat8作为第八颗发射的卫星于2013年2月成功升空，卫星高度705千米，周期为16天，卫星上搭载的OLI陆地成像仪包含9个波段，在可见光波段的空间分辨率为30米，全色波段分辨率可达到15米。目前广泛运用于生态环境监测、土地利用变化、城市热岛效应等领域。表2是Landsat8波段介绍。

表2 Landsat8波段介绍Table 2 Introduction of Landsat8 bands

2 基本方法

本文的总体流程分为三个部分：第一部分是根据Landsat8的光谱特征和纹理特征获取建成区分类图，对Landsat8中第2、3、4、5、6、7(即Blue、Green、Red、NIR、SWIR1、SWIR2)波段与自身的第8(即Pan)波段进行Gram-Schmidt变换(以下简称GS变换)，获得6个空间分辨率为15米的全新波段。对新形成的波段进行波段计算，获取土壤调节植被指数(SAVI),改进的归一化差异水体指数(MNDWI)和改进归一化裸露指数(MNDBI)这三指数波段，分别代表合肥主要的土地覆盖类型：植被、水体和建成区，并对三个波段合并后的影像进行监督分类，获得建成区分类图。第二部分是根据统计年鉴数据对“珞珈一号”进行二分迭代，获取最佳阈值，并对夜光数据进行二值化处理获取建成区范围。第三部分则是对建成区范围和建成区分类图进行决策级上的遥感融合，运用夜光数据提取的建成区范围去裁剪建成区分类图获取最终建成区，具体流程图见图2。

图2 总体流程图

2.1 “三指数法”提取建成区

波段代数运算是遥感影像处理中常用的一种方法，根据地物在不同波段的灰度差异，通过不同波段的代数运算产生新的“波段”，用于突出特定的地物类型。植被、水体、建成区作为地表的主要地物，相关学者根据NDXI模型分别构建了归一化植被指数(Normalized Difference Vegetation Index,NDVI)、归一化水体指数(Normalized Difference Water Index,NDWI)、归一化建筑物指数(Normalized Difference Build-up Index,NDBI)代表植被、水体和建成区，其中裸地因为面积所占比重相对其他三种地物较少，且光谱特征与建成区相似，往往会被划分到建成区中。本文在结合前人理论和建成区内主要地物类型的光谱特征基础上，采用土壤调节植被指数(SAVI)、改进的归一化差异水体指数(MNDWI)和改进的归一化裸露指数(MNDBI)代表上述三类地物。

2.1.1 土壤调节植被指数(SAVI) 植被在可见光波段中对光的吸收明显，尤其是在红波段，同时在近红外波段又具有强烈的反射效果，利用两个波段上的巨大反差，学者提出NDVI指数代表植被，但是NDVI指数对于土壤背景变化较为敏感，只适用于高植被覆盖度地区，因此Huete根据归一化植被指数(NDVI)的原理提出了土壤调节植被指数(SAVI)[13]，其探测植被覆盖率的下限可低至15%。合肥市正处于大发展时期，城市周边大片拆迁和待开发区，建成区边缘地区的植被覆盖率较低，采用SAVI代表植被更具合理性,其计算公式为

SAVI=[(NIR-Red)*(1+L)]/(NIR+Red+L)

(1)

其中L介于0到1之间,通常取0.5，因而本文采用下面的公式来代表植被:

SAVI=[(NIR-Red)*1.5]/(NIR+Red+0.5)

(2)

其中Red和NIR分别为OLI的第4波段(红光波段)和第5波段(近红外波段)。

2.1.2 改进的归一化差异水体指数(MNDWI) 水体的反射率从可见光波段到短波红外2波段的过程中会不断地减弱，Mcfeeters根据NDVI指数的原理提出了归一化水体指数NDWI[14]，但是用NDWI提取有较多建筑物为背景的水体过程中，往往会把建筑物混为水体，使得水体的范围有所扩大，因此徐涵秋提出改进的归一化差异水体指数MNDWI，并且通过大量实验表明MNDWI在提取建成区内部水体过程中有其独特优势[15]，合肥市主城区濒临巢湖，城市内部也有多处水体，采用MNDWI代表水体可以有效地减少建成区被划分为水体的情况，其计算公式为

MNDWI=(Green-SWIR1)/(Green+SWIR1)

(3)

其中Green和SWIR1分别为OLI的第3波段(绿光波段)和第6波段(短波红外1波段)。

2.1.3 改进归一化裸露指数(MNDBI) 采用改进归一化裸露指数MNDBI代表建成区，通过蓝光波段和短波红外2波段进行归一化计算后得到的波段，仅裸露地表(建成区和裸地)结果为正值，其余地物覆盖类型均为负值[16]，采用该指数可大大增加建成区与其余地物的区分度，虽然裸地会被突出，但整体面积较少，且后期采用夜光数据去除，其计算公式为

MNDBI=(SWIR2-Blue)/(SWIR2+Blue)

(4)

其中Blue和SWIR2分别为OLI的第2波段(蓝光波段)和第7波段(短波红外2波段)。

上述的三指数波段，分别表示研究区主要的土地覆盖类型:植被、水体、建成区，对上述获得的三个波段进行合成，对合成后的图像进行监督分类获取建成区分类图。

2.2 二分迭代法提取粗略建成区如何获取夜光数据的最佳阈值，一直是运用夜光数据提取建成区的难点和重点，前人提出了突变检测法、经验阈值法、较高分辨率数据比较法等多种方法，这些方法都各自有其局限性，突变检测法没有考虑城镇发展的区域差异性，而后两种方法则主观性很强，不同的人获得的结果往往是不同的。本文采用了二分迭代的方法获取最佳阈值[17]，具体步骤如下：

第一步，获得提取建成区域夜光数据的最大DN值DNmax和最小DN值DNmin，并对其求均值：

DNt=int [(DNmax+DNmin)/2]

(5)

如果DNt=DNmax，则认为DNt为最佳阈值，如果不满足条件，把此时的DNt作为参数传入下一步。

第二步，根据提供的DNt对夜光数据进行二值化，提取此时与统计年鉴上规定区域相对应的高亮区域作为建成区，计算出其面积area，并与统计年鉴上的建成区面积AREA进行对比，如果

|area-AREA|≤0.01*AREA

(6)

则认为此时的DNt为最佳阈值，如果不满足条件进行下一步。

第三步，判断area和AREA值的大小。如果area

获得最佳阈值后，对影像进行二值化，获取夜光数据高亮区域边界。此时存在两个主要问题：一是周边类似巢湖市(合肥市下属地级市)等地也会被提取出，其与主城区并不相连，应该去除；二是由于夜光数据对于微弱灯光的捕捉效果较好，使得延伸出建成区外部的道路由于路灯的原因被提取出，也需要去除。经过相应处理后，得到粗略的建成区边界。

3 实验与评价

3.1 运用GS变换后的“三指数法”提取建成区结果评价

分别对Landsat8原始影像、GS变换影像、GS变换后的“三指数法”影像进行监督分类，得到图3的建成区分类图，三幅影像图主城区周边都存在许多零碎的斑点，在与Google高清影像进行对比后，可以确认大部分为土壤虚警。与第一幅原始数据的监督分类影像图相比，后两幅进行GS变换后的分类影像中建成区集中于主城区部分，建成区的总体轮廓也更为清晰，土壤虚警也较少。后两幅影像从整体上并没有多大差别，但是仔细对比能发现运用“三指数法”的分类影像图中，由于采用MNDBI指数代表建成区，增大了裸地和建成区的区分度，主城区周边的土壤虚警率大大降低。

图3 Landsat8不同处理方法建成区分类图

在对三幅影像进行分类精度评价过程中，原始影像分类的精度最低(总体精度75.95%，Kappa系数0.5125)，GS变化后精度明显提高(总体精度87.28%，Kappa系数0.7450)，但还是低于GS变化的“三指数法”精度(总体精度93.39%，Kappa系数0.8670)，详细的分类精度情况见表2。

3.2 运用“珞珈一号”夜光数据提取粗略建成区边界结果评价

为了更好地展现“珞珈一号”夜光数据提取建成区的优势，我们分别对NPP/VIIRS夜光数据和“珞珈一号”夜光数据采用二分迭代法获取合肥粗略建成区。由于缺少2018年的DMSP/OLS夜光数据，且前人经常使用的DMSP/OLS夜光数据完全没有优势，因此并未采用该数据进行建成区边界提取效果进行对比。通过查阅《中国建设统计年鉴》，得到合肥市2018年建成区面积为460km2(统计范围仅为合肥市瑶海区、庐阳区、蜀山区(含小庙镇)、包河区、高新区、经开区(含高刘镇)、新站区、政务区)，在此数据基础上分别对合肥市2018年NPP/VIIRS和“珞珈一号”夜光数据进行二分迭代，NPP/VIIRS夜光数据在8次迭代后确认最佳阈值是19，“珞珈一号”夜光数据通过10次迭代得到最佳阈值是14359，根据阈值提取粗略建成区边界。

图4是两种夜光数据提取边界效果比对图，“珞珈一号”夜光数据提取的建成区边界与建成区分类图匹配效果更好，两者的高度匹配性，也从一定程度达到了相互论证的效果，证明建成区分类图和建成区边界提取的正确性。

图4 不同夜光数据提取边界效果比对图

3.3 多源数据融合的效果评价

从监督分类的影像中可以发现许多非主城区的地区也被划分为建成区，而建成区是指城市行政区内实际已成片开发建设、市政公用设施和公共设施基本具备的地区，因此合肥市下属四县一市的建成区应去除。同时使用“三指数法”提取建成区时，虽然选取指数时尽可能地加大了裸地及相关地物与建成区的区分度，但是从图5中b图可发现仍存在把城市外部的部分裸地、农田(可能是大棚种植区)等区域划分为建成区的现象。用建成区边界去裁剪建成区分类图获得建成区，一方面，可去除四县一市的建成区，同时也去除了部分错分区域。这种结合多源数据获得的建成区(d图)相较于单纯使用“珞珈一号”夜光数据获取的建成区(c图)，分辨率提升了近8倍，在建成区内部，夜光数据由于分辨率和夜光外溢效果无法表现的细节信息能够得以表达。表3是详细精度分类表，建成区分类精度(总体精度95.10%，Kappa系数0.9011)有所提高。

表3 建成区分类精度评价Table 3 Classification accuracy evaluation of built-up area

图5 多源遥感数据的效果对比图

表4 结合多源数据提取建成区分类精度评价

4 结论

基于Landsat8和“珞珈一号”夜光数据的建成区提取过程中对Landsat8数据进行Gram-Schmidt变换和“三指数法”等操作获取建成区分类图，并使用夜光数据去除其中裸地等地物，经过相关处理得到图6(合肥市最终建成区)，这种结合多源数据的方式，大大提高了建成区提取的精度。相较于以前的方式，主要价值有以下几点。