张菲菲 阮惠华 许剑辉 赵怡

（1.广东第二师范学院 计算机科学系，广东 广州 510303;2.广东省气象探测数据中心，广东 广州 510080;3.广州地理研究所//广东省遥感与地理信息应用重点实验室//广东省地理空间信息技术与应用公共实验室，广东 广州 510070;4.中国科学院 广州地球化学研究所，广东广州 510640;5.中国科学院大学，北京 100049）

摘要：以粤港澳大湾区为研究区域，以2017年11月1日Sentinel-2A高分辨率多光谱遥感影像为基础数据源，在BAEM指数基础上，结合归一化建筑物指数（NDBI）、归一化植被指数（NDVI）、归一化水体指数（MNDWI）和归一化裸土指数（RNDSI），提出一种针对城市群的高分辨率不透水面综合指数（CompositedISI），并利用Otsu算法实现粤港澳大湾区空间分辨率为10 m的城市不透水面识别与提取，并与传统的NDBI指数提取的不透水面结果进行精度验证与对比。结果表明，CompositedISI指数较好地避免单个指数提取不透水面时所面临的异物同谱等问题。NDBI指数提取不透水面的总体精度和Kappa系数分别为86.80%和0.68，CompositedISI指数提取不透水面的总体精度和Kappa系数达到了88.92%和0.74。相比传统的NDBI指数，CompositedISI指数提取的不透水面效果更优，这是因为CompositedISI指数综合考虑了不透水面与植被、水体、裸体之间差异。

关键词：粤港澳大湾区;Sentinel-2A;归一化建筑物指数;不透水面综合指数

中图分类号：TP75       文献标识码：A

文章编号：1009-3044（2019）27-0279-06

Abstract： To extract 10-m urban impervious surfaces in the Guangdong-Hong Kong-Macao Greater Bay Area （GBA）， the cloud-free Sentinel-2A imagery on 1 November 2017 was used in this study. The normalized difference built-up index （NDBI） was examined to the newer Sentinel-2A data for extracting high-resolution impervious surface during this study， and a new composited impervious surface index （CompositedISI） for impervious surfaces extraction has been proposed. Instead of using individual impervious surface index， CompositedISI may make full use of integrated advantages of NDBI， normalized difference vegetation index （NDVI）， modified normalized difference water index （MNDWI）， and ratio normalized difference soil index （RNDSI）. Then， an impervious surface image with the spatial resolution of 10 m was extracted from the CompositedISI image using the Otsus method. Through integration of NDBI， NDVI， MNDWI， and RNDSI， the CompositedISI was able to improve overall accuracy of impervious surface extraction compared to the NDBI method， with the overall accuracy of 88.92% and a Kappa value of 0.74. The overall accuracy and Kappa of the impervious surface extracted from the NDBI were 86.80% and 0.68， respectively. Results indicate that the CompositedISI was more accurate at mapping high-resolution impervious surfaces when applied to Sentinel-2A imagery as compared to the NDBI method.

Key words： Guangdong-Hong Kong-Macao Greater Bay Area;Sentinel-2A;Normalized difference built-up index;Composited impervious surface index

城市不透水面是指城市中各種不透水建筑材料所覆盖的表面[1]，现已变成了城市环境变化和人地相互作用的指示器与重要驱动力，是城市生态环境规划与保护的关键依据。不透水面的增加会影响城市的水循环，导致洪涝、内涝灾害风险的增加[2]，影响地表蒸腾以及感热潜热的存储和交换、形成更加强烈的城市局地气候，从而进一步加剧城市“热岛”的形成，导致城市生态系统的恶化[3-5]。因此，准确估算与提取高精度的精细不透水面信息对建设中国新型城镇化、生态城市、海绵城市有着重要的现实意义。

随着卫星对地观测技术的快速发展，越来越多的不透水面遥感提取方法被相继提出。当前，利用遥感技术提取不透水面信息的方法主要包括回归方法[6， 7]、光谱解混分析法[8-10]、指数法[11-13]、影像分类法[14，15]等。光谱解混分析方法主要致力于解决中低分辨率影像中存在的混合像元问题[16， 17]。由于其线性关系对地物的解释能力强，且结果可靠，线性光谱混合分析模型（Linear Spectral Mixture Analysis，LSMA）是目前最为常用的光谱混合分析模型[8， 9， 18]。然而，在LSMA 中，仅仅用1或2个固定端元无法有效地代表城市中各种不透水面物质，使得很多裸土会混在不透水面信息中;有限的多光谱影像的光谱分辨率不能有效地获取代表纯净像元的端元光谱特征[19]。分类回归树（classification and regress tree，CART）的方法首先通过高分辨率的数据建立训练样本，然后通过回归树的方法来对中低分辨率大范围的数据进行预测[20， 21]。但结果的准确性很大程度上依赖于高质量的训练样本[20]。指数法是通过遥感光谱指数来增强和提取不透水面信息，可以自动、快速地提取大范围的不透水面。常用的提取不透水面的光谱指数包括：归一化不透水面指数NDISI（Normalized Difference Impervious Surface Index）[1]、城市指数UI（Urban Index）[12]、归一化建筑物指数NDBI（Normalized Difference Built-up Index）[22]、IBI指数（Index-based Built-up Index）[23]、修正的归一化不透水面指数MNDISI（Normalized Difference Impervious Surface Index）[24]、NDII指数（Normalized Difference Impervious Index）[25]、EBBI指数（Enhanced Built-up， and Bareness Index）[26]、BCI指数（Biophysical Composition Index）[27]以及CBI指数（Combinational Built-up Index）[28]等。这些指数都可以用于遥感不透水面提取，并取得了一定的效果。

总的来说，用于不透水面提取的方法很多，但上述大部分方法主要集中于基于中、低分辨率影像的大区域不透水面制图。许多算法都应用到中红外、热红外这些有利于建筑不透水面识别的光谱波段。但是现有的高分辨率影像中除了少数具有中红外波段外，大部分只有可见光和近红外波段，从而限制了这些方法的使用。当前针对高分影像开发的算法还很少，这给高分辨率影像不透水面提取带来较大的困难[19]。

Sentinel-2卫星是全球环境与安全监视系统（GMES）中的多光谱遥感成像任务，该任务的实施由欧洲委员会和欧空局共同执行。该任务用于对全球陆地状况的监测，可用于农业估产，勘察土壤和水的覆盖，内陆河道和海岸地区。Sentinel-2包括两颗A和B星，可以每5天实现对全球一次观测，相较于Landsat 8卫星具有较高的空间分辨率和重访能力。Sentinel-2A/2B卫星已被广泛应用于水体提取、农作物监测和不透水面盖度提取，取得了不错效果。Sentinel-2 影像数据已应用于叶面积指数反演、监视植被健康、水体提取、不透水面盖度提取等[29-33]（郑阳等，2017;谢巧云，2017;方灿莹等，2017;Yang et al.，2017;Du et al.，2017）。就目前来看，针对高分辨率的Sentinel-2影像数据自动提取城市不透水面信息的遥感指数，目前还没有相关的成果发表。

本研究以粵港澳大湾区为研究区域，2017年11月1日Sentinel-2A无云遥感影像为数据源，对BAEM在算法上进行改进，提出一种针对城市群的高分辨率不透水面综合指数（CompositedISI），并利用Otsu算法提取粤港澳大湾区10 m分辨率的不透水面。

1 研究区域与数据

1.1 研究区

本研究以粤港澳大湾区为研究区域（东经112°～115°，北纬21.5°～24°）进行高分辨率不透水面提取研究，如图1所示。粤港澳大湾区海陆兼备，河海交互，涵盖了珠三角9市（广州、深圳、东莞、佛山、中山、珠海、江门、惠州和肇庆）和香港、澳门2个特别行政区。40多年的高速发展，粤港澳大湾区城镇化高度发达，土地扩张明显，建设用地面积增加显著，其中建设用地面积由2000年的4486 km2增长至2016年的6442 km2，人口约7000万[34]，是国家建设世界级城市群和参与全球竞争的重要空间载体，也是世界上最具发展潜力的湾区之一。

1.2 Sentinel-2A影像

本研究以Sentinel-2遥感影像为数据源，提取粤港澳大湾区空间分辨率为10 m的城市不透水面。Sentinel-2是高分辨率多光谱成像卫星，搭载多光谱成像仪MSI（multi spectral  instrument），可用于陆地监测，可提供植被、土壤和水覆盖、内陆水路及海岸区域等图像，还可用于紧急救援服务。Sentinel-2包括2A和2B两颗卫星，重访周期5天。Sentinel-2A包括从可见光和近红外和短波红外通道覆盖13个光谱波段，空间分辨率分为10 m、20 m和60 m，如表1所示。本研究选取2017年11月1日无云的Level-1C Sentinel-2A影像数据，数据来源于Copernicus Open Access Hub（https：//scihub.copernicus.eu）。利用Sen2cor软件包对Sentinel-2A的Level-1C产品进行大气校正、拼接等预处理，并将上述影像数据统一投影到UTM/WGS84中。采用最邻近插值方法对空间分辨率为20m的近红外波段（NIR，Band 8b）和中红外波段1（SWIR1，Band 11）进行空间插值处理，得到空间分辨率为10m的NIR和SWIR1波段反射率。

2 研究方法

2.1不透水面综合指数

BAEMOLI指数是由Bhatti和Tripathi针对Landsat 8 OLI影像数据提出的不透水面综合指数[35]。BAEMOLI指数结合了归一化建筑指数（NDBI）、归一化植被指数（NDVI）、改进的归一化水体指数（MNDWI），综合考虑了不透水面与植被和水体之间的差异，能够快速、自动地提取不透水面，并且在不透水面提取过程中不需要掩膜水体。BAEMOLI指数的计算公式为：

式中，[NDVIOLI]为利用Landsat 8近红外波段（NIR）和红波段（Red）计算的归一化植被指数;[MNDWIOLI]表示利用Landsat 8绿波段（Green）和中红外波段1（SWIR1）计算的改进归一化水体指数;[NDBIOLI]为归一化建筑物指数，利用Landsat 8的NIR、SWIR和热红外波段（TIR）联合计算得到：

[NDBIOLI=1stPCAof Bands6，7+1stPCAof Bands10，11-Band 51stPCAof Bands6，7+1stPCAof Bands10，11+Band 5]    （2）

式中，[1stPCAof Bands6，7]表示Landsat 8 的SWIR1和SWIR2的第一主成分，[1stPCAof Bands10，11]表示Landsat 8 的热红外波段10和11的第一主成分，Band 5表示Landsat 8的NIR波段。然而，由于Sentinel-2A影像数据并没有热红外波段，BAEMOLI指数并不适用于Sentinel-2A提取10 m分辨率的不透水面。参考BAEMOLI指数，本研究基于NDBI、NDVI、MNDWI和RNDSI（归一化裸土指数）构建出一个基于Sentinel-2A高分辨率遥感影像的城市不透水面综合指数：

式中，[NDBInor]、[MNDWInor]、[NDVInor]和[RNDSInor]为经过归一化处理的归一化建筑指数（NDBI）、归一化水体指数（MNDWI）、归一化植被指数（NDVI）和相对归一化裸土指数（RNDSI）。利用Sentinel-2A的SWIR1和NIR波段反射率计算NDBI[36]，并对NDBI指数进行归一化处理：

式中，[TC1nor]表示归一化的裸土亮度信息（Soil Brightness）。采用Landsat 8穗帽变换计算系数，结合Sentinel-2A影像的蓝波段、绿波段、红波段、近红外波段和中红外波段等6个波段信息计算Sentinel-2A影像的裸土亮度信息[TC1nor]：

结合Sentinel-2A遥感影像的6个波段信息，利用式（3）计算研究区不透水面综合指数，获取高分辨率的城市不透水面综合指数影像;利用Otsu算法对高分辨率的城市不透水面综合指数影像进行阈值提取[39]：

式中，m1和m0分别为影像中属于不透水面和透水面像元的平均值，m为整幅影像像元平均值，w1和w0分别为影像中属于不透水面和透水面的比例。

结合式（15）提取的阈值，对高分辨率的城市不透水面综合指数影像进行二值化处理，得到二值分割后的高分辨率城市不透水面影像IS（Impervious Surface）：

式中，（i，j）表示不透水面综合指数影像的像素行列号，threshold为Otsu阈值分割方法自适应选取的分割阈值。

2.2精度评价

不透水面综合指数提取的结果是不透水面和非不透水面的分类信息，可以采用传统的混淆矩阵和Kappa系数来检验分类精度。本研究通过Google高分影像进行交叉验证并随机选择均匀分布的5000个样本（如图2所示），采用总体精度、Kappa系数、生产者精度和使用者精度评价城市不透水面综合指数提取的精度。

3 结果与分析

结合Sentinel-2A遥感影像数据，本研究分别采用CompositedISI和NDBI指数提取空间分辨率为10 m的粤港澳大湾区不透水面，结果如图3所示。从图3可以看出，CompositedISI和NDBI指数能较好地提取粤港澳大湾区的不透水面，不透水面具有类似的倒U型空间分布，主要分布在广州、东莞、深圳、香港特别行政区、佛山、中山、珠海。然而，在肇庆的森林、裸土地区，NDBI会把部分属于植被、土壤的像元被划分为不透水面，而CompositedISI具有较好的结果，如图3左上角所示。

为了更客观地分析高分辨率不透水面提取结果，本研究采用总体精度OA（Overall Accuracy）、Kappa系数、生产者精度（Producers Accuracy）和使用者精度（Users Accuracy）比较CompositedISI和NDBI提取精度，精度结果如表2所示。由表2可以看出，NDBI的总体精度和Kappa系数分别为86.80%和0.68，生产者精度达到91.88%，而使用者精度仅有66.71%;而CompositedISI的总体精度和Kappa系数都比NDBI高，达到了88.92%和0.74。总的来说，CompositedISI用于提取空间分辨率为10 m的城市不透水面是可行的，精度也比较理想。

为了进一步验证本研究提出的方法提取粤港澳大湾区城市群不透水面提取情况，对提取结果进行分区可视化。选择3个典型样区，分别为建筑物密集区域、森林区域和农田区域，如图4所示，其中图4a为高分遥感影像，图4b为本研究提出的CompositedISI指数提取的不透水面，图4c为传统NDBI指数提取的不透水面。

从图4可以看出，CompositedISI指数和NDBI都能较好地提取不透水面。NDBI提取的结果存在一些不足，在树木冠层覆盖比较大的道路或河流沿岸，有部分植被覆盖的像元被分成了不透水面;在森林区域，部分裸地像元被划分为不透水面（如图4c所示），这是不合理的。这是因为NDBI不能很好地区分不透水面和裸土。CompositedISI指数能够从植被、水体、裸土中较好地提取不透水面，道路也能较好地提取。这是因为本研究提出的高分辨率CompositedISI指數考虑了不透水面与植被、水体、裸体之间的差异，避免了NDBI指数提取不透水面时所面临的异物同谱等问题。

[14] Heremans S， Orshoven J V. Machine learning methods for sub-pixel land-cover classification in the spatially heterogeneous region of Flanders （Belgium）： a multi-criteria comparison[J]. International Journal of Remote Sensing， 2015， 36（11）：2934-2962.

[15] Xuefei Hu， Qihao Weng. Estimating impervious surfaces from medium spatial resolution imagery： a comparison between fuzzy classification and LSMA[J]. International Journal of Remote Sensing， 2011， 32（20）：5645-5663.

[16] Lu D， Weng Q. Spectral Mixture Analysis of the Urban Landscape in Indianapolis with Landsat ETM+ Imagery[J]. Photogrammetric Engineering & Remote Sensing， 2004， 70（9）： 1053-1062.

[17] Somers， B， Verbesselt， J， Ampe， E. M， et al. Spectral mixture analysis to monitor defoliation in mixed-aged Eucalyptus globulus Labill plantations in southern Australia using Landsat 5-TM and EO-1 Hyperion data.[J]. International Journal of Applied Earth Observation and Geoinformation， 2010， 12（4）：270-277.

[18] Quarmby N A， Townshend J R G， Settle J J， et al. Linear mixture modelling applied to AVHRR data for crop area estimation[J]. International Journal of Remote Sensing， 1992， 13（3）：415-425.

[19] 徐涵秋，王美雅.地表不透水面信息遙感的主要方法分析[J].遥感学报，2016，20（05）：1270-1289.

[20] Liming Jiang， Mingsheng Liao， Hui Lin， et al. Synergistic use of optical and InSAR data for urban impervious surface mapping： a case study in Hong Kong[J]. International Journal of Remote Sensing， 2009， 30（11）：2781-2796.

[21] Limin Yang， Chengquan Huang， Collin G Homer， et al. An approach for mapping large-area impervious surfaces： synergistic use of Landsat-7 ETM+ and high spatial resolution imagery[J]. Canadian Journal of Remote Sensing， 2003， 29（2）：230-240.

[22] Zha Y ， Gao J ， Ni S . Use of normalized difference built-up index in automatically mapping urban areas from TM imagery[J]. International Journal of Remote Sensing， 2003， 24（3）：583-594.

[23] Xu， H. A new index for delineating built‐up land features in satellite imagery[J]. International Journal of Remote Sensing， 2008， 29（14）：4269-4276.

[24] Liu C ， Shao Z ， Chen M ， et al. MNDISI： a multi-source composition index for impervious surface area estimation at the individual city scale[J]. Remote Sensing Letters， 2013， 4（8）：803-812.

[25] Wang Z ， Gang C ， Li X ， et al. Application of a normalized difference impervious index （NDII） to extract urban impervious surface features based on Landsat TM images[J]. International Journal of Remote Sensing， 2015， 36（4）：1-15.

[26] As-Syakur A R ， Adnyana I W S ， Arthana I W ， et al. Enhanced Built-Up and Bareness Index （EBBI） for Mapping Built-Up and Bare Land in an Urban Area[J]. Remote Sensing， 2012， 4（12）：2957-2970.

[27] Deng C ， Wu C . BCI： A biophysical composition index for remote sensing of urban environments[J]. Remote Sensing of Environment， 2012， 127（none）：247-259.

[28] Sun G ， Chen X ， Jia X ， et al. Combinational Build-Up Index （CBI） for Effective Impervious Surface Mapping in Urban Areas[J]. IEEE Journal of Selected Topics in Applied Earth Observations & Remote Sensing， 2017， 9（5）：2081-2092.

[29] 郑阳， 吴炳方， 张淼. Sentinel-2数据的冬小麦地上干生物量估算及评价[J]. 遥感学报， 2017， 21（2）：318-328.

[30] 谢巧云. 考虑红边特性的多平台遥感数据叶面积指数反演方法研究[D].中国科学院大学（中国科学院遥感与数字地球研究所），2017.

[31] Yang X， Zhao S， Qin X， et al. Mapping of Urban Surface Water Bodies from Sentinel-2 MSI Imagery at 10 m Resolution via NDWI-Based Image Sharpening[J]. Remote Sensing， 2017， 9（6）：596.

[32] Du Y， Zhang Y， Ling F， et al. Water Bodies' Mapping from Sentinel-2 Imagery with Modified Normalized Difference Water Index at 10-m Spatial Resolution Produced by Sharpening the SWIR Band[J]. Remote Sensing， 2016， 354（8）：1-19.

[33] Xu R ， Liu J ， Xu J . Extraction of High-Precision Urban Impervious Surfaces from Sentinel-2 Multispectral Imagery via Modified Linear Spectral Mixture Analysis[J]. Sensors， 2018， 18（9）.

[34] 杨智威， 陈颖彪， 吴志峰， et al. 粤港澳大湾区建设用地扩张与城市热岛扩张耦合态势研究[J]. 地球信息科学学报， 2018， 20（11）：56-67.

[35] Bhatti S S ， Tripathi N K . Built-up area extraction using Landsat 8 OLI imagery[J]. GIScience & Remote Sensing， 2014， 51（4）：445-467.

[36] Zha Y ， Gao J ， Ni S . Use of normalized difference built-up index in automatically mapping urban areas from TM imagery[J]. International Journal of Remote Sensing， 2003， 24（3）：583-594.

[37] 徐涵秋. 利用改進的归一化差异水体指数（MNDWI）提取水体信息的研究[J]. 遥感学报， 2005， 9（5）：79-85.

[38] RNDSI： A ratio normalized difference soil index for remote sensing of urban/suburban environments[J]. International Journal of Applied Earth Observation and Geoinformation， 2015， 39：40-48.

[39] Ostu N . A threshold selection method from gray-histogram[J]. IEEE Transactions on Systems， Man， and Cybernetics， 2007， 9（1）：62-66.

【通联编辑：唐一东】