基于Sentinel-2A与NPP-VIIRS夜间灯光数据的城市建成区提取

刘智丽，张启斌，岳德鹏，郝玉光，苏 凯

(1.北京林业大学林学院，北京 100083；2.中国林业科学研究院沙漠林业实验中心，巴彦淖尔 015200)

0 引言

近年来，中国的城市化水平大幅度提升。经国家统计局统计(http://www.stats.gov.cn/tjsj/)，2017年中国的城市化率已达58.52%，相比2010年增长了11.93%。然而快速城市化在提升人类生活水平的同时，也给城市带来了一系列环境及社会问题，如生态恶化、热岛效应、环境污染、交通拥挤等[1]。城市建成区是一片具有建成环境的连续区域，包括城市用地(道路、建筑及工业设施)及其周边的生态用地(公园绿地、行道树、水体等)，作为衡量城市发展的一个重要指标，其范围及空间分布不仅可以反映城市发展的规模及趋势，如城市扩张水平、人口密度等[2]，还对城市发展规划起着重要作用。因此，对城市建成区的提取得到越来越多的关注[3]。

随着遥感技术的不断发展，城市建成区提取的自动化程度与精度都在不断提升[4-5]。夜间灯光数据(nighttime light data，NTL)能够直接反映城市建成区的夜间亮度，分布连续[6]，但受城市边缘亮度影响较大，因而提取准确度不高。因此，许多研究将NTL与光学遥感影像结合起来提取城市建成区，例如Zhang等[7]将NTL与归一化植被指数(normalized difference vegetation index，NDVI)相结合，得出基于植被校正的城市夜间灯光指数(vegetation adjusted NTL urban index，VANUI)，有效丰富了城市边缘信息[8]，但VANUI会放大水体的反射率，混淆城市边缘地带的水体与城市建成区，且空间分辨率不够高，提取结果仍无法满足更高精度的研究需求。

欧洲航天局发射的Sentinel-2A光学遥感卫星覆盖了可见光、近红外以及短波红外范围内的13个波段，空间分辨率最高可达10 m[9]，能够准确识别城市建成区的边界，可为城市边界提取提供良好的数据支持。但目前将Sentinel-2A数据应用于城市建成区识别、提取的研究较少，因此，本文对VANUI指数加以改进，基于Sentinel-2A与NTL数据提出基于建筑校正的城市夜间灯光指数(building adjusted NTL urban index，BANUI)，以包头市南部为研究区，进行城市建成区的提取。

1 研究区概况与数据源

1.1 研究区概况

包头市位于内蒙古自治区中西部，总面积约为27 768 km2，平均海拔为2 000 m，中部地势较高，由中部向南北两侧地势逐渐变低，横贯其中部的阴山山脉将包头市分成中部山区、北部高原草场与南部平原地区3个部分，其中南部地区包括包头市中心城区与土默特右旗。近年来，包头市经济发展迅速，南部平原地区城市规模不断扩张，在2000—2016年间，包头市城区面积增长了71.10%。但当地气候干燥、地形起伏较大，生态环境脆弱，而城市的不断扩张又加剧了对生态环境的破坏，由此引发了一系列环境问题。包头市南部地区聚集了全市89.69%的人口，城市分布密集，具备进行城市建成区遥感监测的基础，而中部和北部地区大多为高原山地和草原，人口稀疏、城市数量屈指可数，缺乏城市遥感监测的条件，因此将包头市南部平原地区作为本次研究的研究区(图1)。

图1 研究区遥感影像Fig.1 Image of study area

1.2 数据源

本文所使用的数据有Sentinel-2A遥感影像数据、NPP-VIIRS NTL数据以及研究区2016年土地利用数据。遥感影像选用处理后的2016年6月份研究区Sentinel-2A Level-2A数据，来源于欧洲航天局网站(https://scihub.copernicus.eu/)，经Sen2Cor模型处理，包含12个波段。

目前使用较多的NTL数据有DMSP-OLS和NPP-VIIRS 2种，其中，NPP-VIIRS数据的空间、时间分辨率都比较高，并已经过辐射定标，更适宜应用于分析与研究[12]。因此，本文选用2016年6月份包头NPP-VIIRS数据，下载自美国国家海洋大气局网站(https://www.ngdc.noaa.gov/eog/viirs/download_dnb_composites.html)。数据预处理时将亮度值小于0的异常像元赋值为0。

研究区2016年土地利用调查数据来源于包头市环境保护局，但其中缺少城市建成区的地面实况数据，因此基于此数据通过手动标注来获取城市建成区的范围，作为本次研究的验证数据。

本文中使用的所有地理空间数据投影均为WGS1984 UTM投影。

2 研究方法

2.1 技术路线

本文提出BANUI指数用于研究区城市建成区的提取。首先，基于Sentinel-2A数据得到研究区的归一化建筑指数(normalized difference build-up index，NDBI)，将其与NPP-VIIRS数据结合得到BANUI数据；然后，基于BANUI利用分水岭分割算法提取包头市城市建成区；最后，对提取结果进行对比分析和精度验证。技术路线如图2所示。

图2 技术路线Fig.2 Technology roadmap

2.2 VANUI

为了将城市与农村更好地区分开来，增大两者亮度值的差异，Zhang等[14]根据DMSP/OLS数据与NDVI提出VANUI指数，计算公式为

VANUI=(1-NDVI)NTL，

(1)

式中NDVI为归一化植被指数。当VANUI接近0时，区域植被覆盖较好；VANUI的值较大时，植被相对稀疏。但是计算NDVI时使用的MODIS数据空间分辨率仅为500 m，NTL数据中的细节信息无法加以利用；且在VANUI中，非植被覆盖地区的亮度值会被过度放大，水体反射与其他异常亮度点也会对提取精度有所影响。

2.3 BANUI

2.3.1 NDBI的提取

由TM影像中提取的NDBI可以准确描述城市建筑用地信息，计算公式为

(2)

式中TM4和TM5分别为TM影像的近红外和短波红外波段的反射率。一般情况下，NDBI=0时，为耕地或林地；NDBI<0时，为水体；NDBI>0时，为建筑用地。

处理后的Sentinel-2A数据包含12个波段，不同波段的组合会导致城市建成区范围有所差异，因此分析不同地类(建筑、植被、水体、未利用地)上的光谱特征并确定最佳波段组合具有重要意义。本文利用感兴趣区(region of interest，ROI)工具选取研究区各地类的代表性像元，并利用ArcGIS10.2软件进行分区统计，结果如图3所示。

图3 Sentinel-2A数据各波段的平均亮度值Fig.3 Mean brightness of each band in Sentinel-2A

参照式(2)中所使用波段的波长，本文从Sentinel-2A数据的第8、第11与第12波段中选出合适的波段进行组合。由图3可知，在第11和12波段上，建筑用地的像元亮度值较高，水体和植被的亮度值远小于建筑用地；在第8波段上，植被的亮度值达到极值，水体与建筑的亮度值相对较小；在第11、第8波段或第12、第8波段的亮度差值中，建筑为正数，植被为负数，水体为接近于0的负值。利用这3个波段的组合构建了2个NDBI，分别记作NDBI11和NDBI12，公式为

(3)

(4)

式中：S8,S11,S12分别为Sentinel-2A数据中第8、第11与第12波段的反射率。

为了直观地对比NDBI11与NDBI12在不同地物上的表达差异，选取涵盖了建筑、植被、水体与未利用地4种类型的矩形区域，并截取对应的NDBI图像，结果如图4所示。图4显示，NDBI11和NDBI12都能够较好地反映建筑用地的范围，建筑的亮度值明显高于图中其他2种地类，水体的亮度值在图中并不高，植被的亮度值也明显低于建筑，且植被覆盖程度越高的地方亮度值越低。此外，2幅图植被亮度差异不明显，NDBI12中建筑的亮度略高，因此NDBI12中建筑与植被更容易区分；但NDBI12中水体的DN值略高于NDBI11的水体，这在一定程度上会干扰建筑用地的提取。为定量分析2种波段组合所获取的NDBI的优劣，本文计算出各种地类的NDBI概率分布并绘制其正态分布曲线，结果如图5所示。由图5可知，2幅图中建筑用地的NDBI值最高，水体次之，植被的NDBI值最低，其值远小于0。2幅图像中建筑用地与植被的差别都比较明显，但在图5(b)中，建筑与水体有较多重叠区域，此外建筑用地与未利用地在2幅图中都难以区分。

(a)NDBI11(b)NDBI12

图4 NDBI11和NDBI12的代表性区域

Fig.4RepresentativeareasofNDBI11andNDBI12

(a)NDBI11(b)NDBI12

图5 不同地类在NDBI中的概率分布及其正态分布曲线

Fig.5ProbabilitydistributionandnormaldistributioncurvesofdifferentlandtypesinNDBI

基于上图中正态分布曲线以及色谱分离分析理论可以计算建筑与其余各种地类之间的分离度R，目前该指数被广泛应用于色谱柱中2个相邻元素间分离程度的计算。一般当R<1时，2波峰有部分重叠；R=1时，其分离程度可达98%；R=1.5时，分离程度接近99.7%，此时2个相邻元素完全分离；R=2时，2波峰从基线起就完全分离。通常情况下，利用色谱峰的保留时间以及基线上的峰宽可计算分离度R，本文则借助各地类正态曲线峰值所对应的NDBI值计算分离度，公式为

(5)

式中：BR1和BR2分别为2条正态曲线的峰值所在NDBI的值；W1和W2分别为2条正态曲线在基线上的峰宽。计算结果如表1所示。

表1 NDBI11与NDBI12中建筑与各地类之间的分离度RTab.1 R between buildings and the other land types in NDBI11 and NDBI12

由表1计算结果可知，建筑与植被的分离度均大于2，表明2种地类的正态曲线完全分离，在NDBI中差异明显，但NDBI12中两者的分离度更高；建筑与水体在NDBI11中的分离度高于NDBI12中的分离度，因此NDBI11中建筑与水体差别明显；建筑与未利用地的分离度都比较低、不易区分。由此可以推断，除了建筑与未利用地，其他地类信息在NDBI11中都得到了不同程度的抑制，而NDBI12对水体信息的抑制作用并不明显，NDBI11与NTL数据的结合可以有效消除噪声的影响、并增强建筑信息的表达。因此本文选择Sentinel-2A数据中第8波段和第11波段来提取NDBI。

2.3.2 BANUI的提取

在VANUI中，水体的值相对偏离正常范围，且空间分辨率较低。为解决这一问题，本文利用NDBI与NPP-VIIRS数据构建了BANUI指数，公式为

BANUI=(1+NDBI)VNTL，

(6)

式中VNTL为NPP-VIIRS数据的亮度值。经过重采样后的NPP-VIIRS数据的空间分辨率为20 m，最终BANUI的分辨率同样为20 m。为检验BANUI的提取效果，同样利用Sentinel-2A与NPP-VIIRS数据计算得到VANUI，其空间分辨率为10 m。

选取E109.35°～110.85°范围内一横截面，该横截面可以最大程度地经过包头市南部城市建成区，信息丰富，基于此可检验BANUI中建筑信息的表达。横截面上VANUI与BANUI的值如图6所示。由图6可知，两者变化趋势比较相似；横截面经过包头市一公园，公园内部植被覆盖度较高，且因NDVI对植被更敏感，导致VANUI的值要比BANUI更低些。

图6 横截面上VANUI与BANUI的值Fig.6 Values of VANUI and BANUI on the cross section

表2数据为研究区内VANUI和BANUI的像元亮度值范围、像元亮度均值、像元亮度值标准差，以及水体、植被、建筑的像元亮度均值。由于未利用地距离城区较远，其面积只有153.47 km2(研究区总面积为4 965.78 km2)，对于城市建成区的提取影响甚微，因此未利用地的数据并未包含在其中。

表2 VANUI和BANUI的统计数据Tab.2 Statistical data of VANUI and BANUI

从表2可以看出，BANUI像元亮度值覆盖范围较广、标准差更高，其变化幅度更大，所涵盖的信息更多；水体在BANUI中的像元亮度均值较低，建筑在BANUI中的像元亮度均值较高，表明BANUI更好地凸显了建筑信息。

图7为包头市城市建成区边缘地带，包含了水体、植被以及连续分布建筑等信息，因此将其作为检验BANUI区分建筑用地效果的典型区域。

(a)典型区域位置 (b)典型区域影像

图7 典型区域

Fig.7Typicalarea

该区域水体、建筑、植被的概率分布以及正态分布情况如图8所示。

(a)VANUI (b)BANUI

图8 典型区域内不同地类在VANUI与BANUI中的概率分布及正态分布曲线

Fig.8ProbabilitydistributionandnormaldistributioncurvesofdifferentlandtypesofVANUIandBANUIinthetypicalarea

图8显示，BANUI中建筑与植被的像元亮度值要高于VANUI中建筑与植被的亮度值，水体的亮度值在两幅图像中相差不大，导致图8(b)中建筑与水体距离更远，二者更容易区分。

计算VANUI和BANUI中建筑与其他地类的分离度，结果显示，建筑与植被的分离度都超过了1.5(1.76和1.52)；而BANUI中建筑与水体的分离度(0.63)明显高于VANUI中二者的分离度(0.29)，表明BANUI相比VANUI能够更好地将建筑与水体区分开来。由此可见，BANUI可以丰富NPP-VIIRS的信息、更好地抑制水体反射造成的影响，并增强建筑信息的表达，因此将其视作提取城市建成区的有效数据。

2.4 分水岭分割算法

为检验BANUI在城市建成区提取应用中的优势，本文利用分水岭分割算法从NPP-VIIRS与VANUI数据中提取城市建成区。

分水岭分割算法是一种基于区域进行影像分割的有效方法，可用于城市建成区范围的提取[15]。BANUI数据的像元被视作地形表面，每个像元的亮度值为地形表面上相应像元的高程值，亮度值高的像元为高地，亮度值低的像元为低地，城市建成区的边界为地形中集水盆地的边界即分水岭[16]。

主要步骤如下：①首先，对影像进行平滑处理，去除BANUI影像中椒盐噪声的影响，自适应加权中值滤波方法可在去除噪声的同时保留更多信息，而高斯滤波、中值滤波等常用方法则无法满足这一需求；②利用形态重构提高影像中同类物体的同质性，并降低不同物体之间的异质性，由此得到形态重构影像；③计算影像的形态梯度，通常情况下像元亮度值变化程度较高的地区，影像梯度较大；④借助形态重构影像求得局部极大值并标记为前景对象，通过欧氏距离倒数法标记背景影像；⑤利用分水岭分割算法划分标有前景与背景对象的待处理影像，提取出城市建成区。

2.5 对比分析与精度验证

为了检验本文提取城市建成区所采取方法的有效性，选取包头市2016年的土地利用调查数据作为对照，通过人机交互方式在其基础上绘制出城市建成区的范围作为验证数据，此范围具有高精度、连续性高等特点。

通过实验结果与验证数据之间逐像元的对比建立混淆矩阵，进而计算总体精度、Kappa系数、生产者精度和用户精度等指标进行精度检验[17]。此外，需借助景观指数来描述结果空间结构上的特征，包括总面积(total area，TA)、景观形状指数(landscape shape index，LSI)和聚集指数(aggregation index，AI)等。

LSI可衡量整个景观的形状与面积相同的正方形或圆形之间的偏差程度，以此来表示该景观整体形态的复杂程度[18]，取值范围为[1，+∞)。LSI值越大，该景观形状与正方形的偏差越大。公式为

(7)

式中：0.25为正方形校正常数；E为景观中所有斑块的边缘总长度，100 m；A为景观总面积，hm2。

AI表示景观内部斑块之间相临的格网单元数目与在这些数中的最大值的比值，它反映了景观的聚集程度与连续性[19]。公式为

(8)

式中gi为景观内斑块i相邻的格网单元数。AI的取值范围为(0,100%]，当景观内部斑块无限分散时，AI的值趋向于0；随着斑块的逐渐集聚，AI的值也随之增加。

3 成果分析与检验

3.1 提取结果

利用分水岭分割算法由不同数据中提取的城市建成区(图9)。由图9可知，利从BANUI中提取的城市建成区与验证数据相似度较高，而从VANUI和NPP-VIIRS中提取的城市建成区相比于验证数据范围较大。由BANUI中提取的城市建成区包括城市用地，以及水体、公园等其他地类，空间分布连续、成分完整，与上文对城市建成区的定义相吻合。

(a)NTL (b)VANUI (c)BANUI

图9 由不同数据提取的城市建成区

Fig.9Extractedurbanbuilt-upareasfromdifferentdata

3.2 对比分析与精度验证

通过上文所述精度验证方法得到不同数据的提取精度，如图10所示。从图10可以看出，由BANUI提取出的城市建成区总体精度最高，为93.61%，Kappa系数为0.793 4，城市建成区的生产者精度与用户精度分别为85.34%和81.34%，总体上提取的精度较高；VANUI位列其次，NPP-VIIRS精度最低。

图10 不同数据源提取的城市建成区的精度对比Fig.10 Accuracy comparison of urban built-up areas extracted by different data sources

表3为不同数据源提取的城市建成区景观格局指数对比，其中，TA.Diff，LSI.Diff，AI.Diff分别表示TA，LSI，AI与验证数据的误差。由表4可知，BANUI提取结果的景观指数与验证数据十分接近，提取的城市建成区总面积比验证数据多了4.91%，LSI和AI的差距也较小，提取效果最佳；而其余2个数据源所提取结果与验证数据相差较大，面积分别多出了24.12%和16.55%，精度较低。

表3 不同数据源提取的城市建成区景观格局指数对比Tab.3 Comparison of landscape pattern indexes of urban built-up areas extracted by different data sources

4 结论

1)基于Sentinel-2A与NPP-VIIRS数据，对传统的VANUI进行改进，提出了用于提取城市建成区的BANUI指数。

2)通过对Sentinel-2A数据各个波段的特征进行分析，选择Sentinel-2A影像的第8和第11波段计算研究区NDBI，从而构建BANUI指数

3)实验结果表明，由BANUI所提取的研究区城市建成区的面积为971.25 km2，相比验证数据仅相差4.91%，总体精度高达93.61%，Kappa系数为0.793 4，明显优于由VANUI和NPP-VIIRS数据的提取结果。

4)BANUI在城市建成区的提取上具有结果完整、覆盖范围广、提取精度高等优点，应用前景广阔，可为进一步研究NTL数据对城市建成区的提取提供思路，也可用于城市发展的规划与监测。