基于Sentinel-2数据的不透水面识别方法与应用研究

徐 瑞，屈 冉，刘 晓 龙，马 慧 榕，史 正 涛

(1.云南师范大学地理学部，云南 昆明 650500；2.生态环境部卫星环境应用中心，北京 100094；3.内蒙古工业大学信息工程学院，内蒙古 呼和浩特 010051)

0 引言

不透水面[1](如沥青道路、房屋、广场、停车场、塑料大棚等)在扩张过程中占用了大量耕地、林地及水域，不仅导致生物多样性锐减，影响人与自然和谐相处[2]，还影响城市水循环系统，改变地表径流，增加城市内涝灾害风险[3,4]，且加剧了城市“热岛”形成，导致城市生态系统恶化[5-7]。因此，快速、准确获取不透水面的空间分布信息对于制定生态环境保护政策、合理规划土地资源、构建新型生态城市以及海绵城市等具有重要意义。

遥感技术时效性高、覆盖范围广，被广泛应用于不透水面提取。面向对象法[8]、影像分类法[9]、光谱混合分析(Spectral Mixture Analysis，SMA)[10]、分类回归树(Classification and Regression Trees，CART)[11]、人工神经网络(Artificial Neural Network，ANN)[12,13]、深度学习[14]及指数法[15,16]常用于不透水面识别。面向对象法基于训练样本对影像进行分割分类，从而提取不透水面，主要适用于高分辨率遥感影像[8];影像分类法(如最大似然法、平行管道法、支持向量机等)依据不透水面光谱信息利用传统分类器进行自动分类，操作性强[17]，主要针对中低空间分辨率光学影像(如MODIS[18,19]、Landsat[20,21])，但存在大量混合像元，降低了提取精度;SMA将影像像元分为高、低反射率对象、土壤和植被4种端元，通过生成各端元的灰度图提取不透水面，能解决混合像元导致分类精度低的问题[10]，适用于大规模城市(如京津冀、长三角等城市群)的不透水面提取[22-24];ANN适用于识别范围较小的城市不透水面，且分类精度优于传统的统计分类方法[12,13]，但存在大量隐含层数和参数，求取最佳学习速率和动量系数耗时长[14];深度学习能辅助ANN处理参数，提取高层次的特征数据，常被应用于不透水面提取[25,26];CART算法通过训练样本形成的规则预测连续变量提取不透水面[11]。上述方法均建立在模型学习的基础上，对训练样本依赖性较强，一直面临模型训练成本与精度权衡这一主要问题[27];而指数法利用不透水面的光谱响应特征，通过构建遥感光谱指数提取不透水面信息[28]，无需训练样本，操作简便、计算效率高[15,16,28]。常用的遥感指数有归一化建筑指数(Normalized Difference Built-up Index，NDBI)[29]、遥感建筑用地指数(Index-based Built-up Index，IBI)[30]、归一化差值不透水面指数(Normalized Difference Impervious Surface Index，NDISI)[31]、生物物理综合指数(Biophysical Composition Index， BCI)[32]、比值居民地指数(Ratio Resident-area Index，RRI)[33]、不透水面综合指数(Combinational Build-up Index，CBI)[34]、改进型归一化差值不透水表面指数(Modified Normalized Difference Impervious Surface Index，MNDISI)[15]等。

Attarchi等[35-38]学者尝试应用合成孔径雷达(SAR)和夜间灯光等遥感影像提取不透水面。其中，SAR数据具有穿云透雾、不受天气影响等优点，但在复杂地形区提取不透水面仍面临挑战[39]；夜间灯光传感器(如DMSP-OLS、VIIRS)基于其独特的光学放大能力[40,41]，可以根据灯光亮度阈值提取不透水面，但在偏远地区和欠发达地区仍存在不足[42]。与常用的陆地卫星相比，Sentinel-2空间分辨率更高(Red、Green、Blue、NIR波段空间分辨率为10 m，SWIR波段空间分辨率为20 m)，重访周期短(单星重访周期为10 d，双星重访周期为5 d)。因此，本文基于Sentinel-2数据在不透水面提取方面的良好性能以及指数法的优点[43,44]，提出一种结合多种指数和主成分分析的不透水面识别方法(Impervious Surface Identification Method，ISIM-S2)，以满足快速、准确获取局部或大范围不透水面分布信息的需求。

1 研究区和数据

1.1 研究区概况

不透水面的分布与扩张受地形地貌、气候(降水)的影响，与地理环境特征密切相关。为验证ISIM-S2方法的普适性，选取中国内陆不同气候区、地貌类型、城市环境的4个典型城市(北京、长春、广州和昆明)作为研究区。北京位于华北平原北部(115°42′～117°24′E，39°24′～41°36′N)，地势西北高、东南低，属于暖温带半湿润半干旱季风气候，冬季严寒干燥，夏季高温多雨；长春地处东北平原地区(124°18′～127°5′E，43°5′～45°15′N)，属于温带大陆性半湿润季风气候，雨热同期，四季分明，地势起伏小(地表相对高差不超过50 m)，有利于城市基础设施建设；广州地处南部沿海地区(112°57′～114°3′E，22°26′～23°56′N)，中部以中低山和丘陵为主，东部以山地为主，南部以平原为主，属于海洋性亚热带季风气候，冬暖夏凉，雷雨较多；昆明地处云贵高原中部(102°10′～103°40′E，24°23′～26°22′N)，总体地势北高南低，属于亚热带高原山地季风气候，全年温差较小，干、湿季分明。

1.2 数据及预处理

1.2.1 实验数据 1)Sentinel-2影像源自哥白尼开放访问中心(https://scihub.copernicus.eu/)，选用无云或云量较少的影像，北京、长春、广州影像获取卫星为Sentinel-2A，获取时间分别为2017年7月7日、2017年6月28日、2017年12月31日，昆明影像获取卫星为Sentinel-2B，获取时间为2018年1月7日。本文使用欧空局发布的Sen2Cor处理模块(http://step.esa.int/main/third-party-plugins-2/sen2cor/)对Sentinel-2 1C级产品进行辐射定标和大气校正，并采用最近邻法将20 m空间分辨率的短波红外(SWIR)波段数据重采样至10 m。2)高程数据(空间分辨率为90 m)和气候分区数据分别来源于https://srtm.csi.cgiar.org/srtmdata/和KÖPPEN-GEIGER气候分类世界地图(http://koeppen-geiger.vu-wien.ac.at/)。

1.2.2 验证数据 为评价不透水面提取结果的准确性，本文按随机分层抽样原则选择验证样本。尽管当前有不透水面分类产品(如GAIA[45]、GLASS-GLC[46])和高分辨率影像(如Google Earth中的遥感影像)可作为验证样本，但裸土和城市不透水面的快速变化使上述产品的时效性不足，不能直接作为本文提取结果的验证参考。为此，本文通过对Sentinel-2影像进行目视解译选取验证样本，然后根据成像日期接近的Google Earth影像对验证样本进行调整，以确保不透水面验证样本的准确性。最终获取北京、长春、广州和昆明4个研究区的验证样本数分别为217个、198个、213个和214个。

2 不透水面提取方法

ISIM-S2方法主要包括不透水面靶区(含裸土和不透水面)提取和靶区内不透水面识别两部分，具体流程如图1所示。

图1 ISIM-S2方法流程Fig.1 Flow chart of ISIM-S2 method

2.1 不透水面靶区提取

遥感指数可量化地表生物的物理特征，增强地类间的差异性，但由于不同地表覆盖类型(如裸土和不透水面)存在混淆，导致单一遥感指数的有效性受到限制[47]。为最大限度地提高遥感指数性能，本文采用多个遥感指数逐步分离不同土地覆盖类型。1)水体识别。改进的归一化水体指数(Modified Normalized Difference Water Index，MNDWI)[48]不仅能有效提取城市区域内的水体，还能区分水体中的阴影以及水体中的悬浮沉积物等细微特征，因此，本文选择MNDWI指数提取水体(式(1))，即将MNDWI值大于0的像元识别为水体。2)云/云阴影识别。采用Fmask算法[49]识别Sentinel-2 1C级影像中的云/云阴影，该算法下载自GitHub官网(https://github.com/GERSL/Fmask)，算法参数均为默认值。将上述已识别的水体、云/云阴影进行掩膜后，研究区内剩余植被、裸土和不透水面3种地类。3)不透水面靶区提取。NDBI指数主要根据不同地类的光谱特征差异增强不透水面信息[29]，虽然低覆盖植被、裸土及不透水面之间存在异物同谱现象，导致NDBI指数将部分裸土与不透水面和低覆盖植被混淆，但NDBI指数具有易操性[47]；而SAVI[50]能抑制裸土和背景噪声，与NDVI相比，更适用于分离植被和裸土。因此，本文采用NDBI-SAVI[47]提取裸土和不透水面(式(2)、式(3))。

(1)

(2)

(3)

式中：ρGreen、ρSWIR、ρNIR、ρRed分别为绿光(Green)、短波红外(SWIR)、近红外(NIR)和红光(Red)波段的地表反射率；L为土壤调节因子，本文取0.5。

NDBI-SAVI图像的灰度值范围为-2.5～2.5，其中不透水面靶区像元的灰度值大于植被像元的灰度值,为分离不透水面靶区和植被，应设置最佳阈值。OTSU[51]算法通过阈值将待处理图像分为前景与背景，取最佳阈值时，分类效果较好。因此，本研究采用OTSU算法确定NDBI-SAVI灰度图像的最佳阈值。

2.2 靶区内不透水面识别

主成分分析(PCA)是一种有效的特征提取工具，可在特征数量显著减少的同时，保证每个主成分丢失的原始数据信息较少，且不同主成分包含的信息唯一。因此，本文利用PCA和OTSU阈值分割方法将裸土和不透水面分离。实验表明：北京、长春采用PCA2，广州、昆明采用PCA3，OTSU阈值均为0.01。PCA系数可能随影像的变化而变化(如场景范围、土地覆盖类型、波段数、图像采集日期等)，表明PCA衍生的特征不稳定，主要是因为地表覆盖物(如植被、水体)在某些特定波段或所有波段具有高反射率(或低反射率)特征，而ISIM-S2方法基于分层分类思想去除裸土和不透水面以外的地物，可降低PCA衍生特征的不稳定性。

3 不透水面提取结果评价

本文采用BCI、CBI、RRI指数基于OTSU阈值分割法提取研究区不透水面并与ISIM-S2方法进行对比，采用总体精度(OA)、用户精度(UA)、生产者精度(PA)和Kappa系数进行精度评估(图2)，并对各研究区全局和局部区域不透水面空间分布进行比较分析(图3)。

图2 ISIM-S2与3种指数在不同研究区的分类精度对比Fig.2 Classification accuracies comparison of ISIM-S2 with three indexes in different study areas

由图2可知，ISIM-S2方法在4个研究区的总体分类精度(OA)均大于90%，Kappa系数均大于0.80，且该方法的OA和Kappa系数均高于其他3种指数，提取不透水面的精度最高，显著减少了裸土与不透水面的混淆。但ISIM-S2方法也存在如下不足：会将部分暗色不透水面(水泥、沥青等)误分为透水层(阴影、水体)；面积较小的不透水面边界常与其他地物形成混合像元，从而易被误分为透水层(图3长春研究区中)；部分有色屋顶(以蓝色为主)与植被在红光波段光谱值差异较小[52]，易被误分为透水层。

图3 ISIM-S2与3种指数在不同研究区的不透水面提取结果对比Fig.3 Result comparison of impervious surface extraction using ISIM-S2 with three indexes in different study areas

由4个研究区内提取结果对比(图3)可知：1)由于云与不透水面等部分高反射率地物的反射率相似，导致3种指数均会将云识别为不透水面，BCI和RRI指数的误识别情况更严重，而ISIM-S2在提取靶区时先用Fmask算法将云/云阴影去除，降低了云对不透水面的影响，识别错误较少。2)3种指数均存在将水体误识别为不透水面的情况，其中RRI指数的误识别最严重，而ISIM-S2在靶区提取阶段将水体去除，降低了误分率。3)4种方法均存在将不透水面与裸土混淆的情况，其中，CBI将裸土误识别为不透水面的概率最高，ISIM-S2针对靶区范围的裸土和不透水面进行主成分分析，混淆程度最低。

综上，相比BCI、CBI、RRI指数，ISIM-S2方法降低了云、水体对不透水面的影响，有效改善了裸土和不透水面提取混淆情况，且在4个研究区内总体识别精度均属最高。

4 结论

不透水面与裸土的混淆一直是基于遥感数据识别不透水面的主要难题，为降低裸土对不透水面提取的影响，本文提出一种基于Sentinel-2数据的不透水面识别方法(ISIM-S2)。首先利用Fmask、MNDWI、NDBI和SAVI 指数剔除云/云阴影、水体和植被，得到不透水面靶区；然后采用主成分分析法对靶区内的裸土和不透水面进行分离，得到不透水面提取结果。将该方法应用于4个不同地理环境下的研究区(北京、长春、广州和昆明)，并与BCI、CBI、RRI指数的提取结果进行对比，得出以下结论：1)ISIM-S2方法在4个研究区内的不透水面总体分类精度(OA)均大于90%，Kappa系数均大于0.80，不透水面提取精度较高。2)相比BCI、CBI、RRI指数，ISIM-S2预先去除云/云阴影和水体，有效抑制了高(低)反射率地物像元值变化导致PCA衍生特征的不稳定性，不仅降低了云、水体对不透水面的影响，还能有效改善裸土与不透水面提取混淆问题，在4个研究区提取的总体精度(OA)和Kappa系数均高于BCI、CBI、RRI指数。

本文提出的ISIM-S2方法无需训练样本，操作简单、高效且精度较高，可为不透水面识别提供一种新途径。但由于不同地理环境下城市分布、不透水面像元占比和裸土反射特性存在差异，如何实现靶区内PCA裸土和不透水面区分阈值的自适应设置，后期需要深入研究。