基于Sentinel-2A/B的新疆典型城市不透水面提取及空间差异分析

段潘，张飞，3，刘长江

1.新疆大学 地理与遥感科学学院,乌鲁木齐 830046;2.新疆大学绿洲生态教育部重点实验室,乌鲁木齐 830046;3.新疆大学 智慧城市与环境建模自治区普通高校重点实验室,乌鲁木齐 830046

1 引 言

不透水面ISA（Impervious Surface Area）是指城市下垫面中由各种不透水建筑材料所覆盖的表面（徐涵秋，2008；匡文慧，2019）。由于全球快速城市化和气候变暖等严重影响了城市的人居环境、生态系统以及人类福祉等（Grimm 等，2008；Creutzig 等，2016）。与此同时，城市化的加速导致人工下垫面的连续扩张，城市不透水面不断取代以森林、草地、湿地等植被为主的自然地貌，由此导致城市内涝、城市热岛效应加剧以及城市生态系统服务减弱等（林云杉等，2007；Yuan 和Bauer，2007；徐涵秋，2009；周玄德等，2018）。因此，研究城市不透水面空间结构变化规律，对于城市生态环境保护和城市绿地空间规划都有着重要意义（Dou和Kuang，2020）。

国内外学者对不透水面的提取做了大量的探索和研究，早期有学者利用传统实地测量的方法获取不透水面的分布，但该方法费时费力，实时性差，难以实现城市不透水面的动态更新（Slonecker 等，2001）。近年来，随着遥感技术的发展，Wu 和Murray（2003）在Ridd（1995）提出的城市植被—不透水面—土壤的混合像元分解方法基础上，成功地将该算法应用到ETM+影像中，周存林和徐涵秋（2007）、王兆华等（2019）也利用光谱混合分析分别提取了兰州和福州城区的不透水面，但该方法提取精度很大程度上受限于端元提取的准确性。邵振峰等（2016）基于随机森林模型提取武汉城市不透水面，刘波等（2017）以高分二号为数据源采用支持向量机方法提取梧州城市不透水面，基于对象的机器学习方法，优化分割参数的难度较大，应用在大范围研究区域比较困难。与上述方法相比，光谱指数法无需设置参数，方便应用于实际生产中。Kawamura 等（1996）研究了不透水面的光谱特征，并提出了城市指数UI（Urban Index）来代表城市地区。之后，大量光谱指数被提出用于描述不透水面，其中代表性光谱指数有：Zha 等（2003）、徐涵秋和王美雅（2016）和穆亚超等（2018）分别提出的归一化建筑指数NDBI（Normalized Difference Built-up Index）、归一化差值不透水面指数NDISI （Normalized Difference Impervious Surface Index）和增强型归一化差值不透水面指数（ENDISI）。

不透水面是典型的土地覆盖类型，是衡量和分析城市发展和生态环境影响的重要指标（徐涵秋，2009）。目前不透水面在城市空间扩张变化方面有大量应用研究，Yang 和Liu（2005）提出用城市不透水面来体现城市扩张变化，Peng 等（2016）从不同功能区、环线结构、剖面线等多方面分析城市发展变化；由国内外研究可知，遥感技术的快速发展极大的促进了不透水面的研究（Deng 和Zhu，2020；Wang 等，2019）。但近年来不透水面提取多以30 m 的中等分辨率影像为基础，分析城市不透水面的面积和速率变化（包颖等，2020；张稼乐 等，2020；Gong 等，2020）。对城市不透水面空间差异变化规律研究较少。因此，本文以新疆典型城市（乌鲁木齐、喀什、哈密、克拉玛依）为研究区，借助10 m 分辨率的Sentinel-2A/B数据，通过L2范数归一化处理结合ENDISI以最大类间方差法自动确定最优阈值提取不透水面，利用剖面线法和盒维数法研究2017年—2019年城市不透水面空间差异变化。探讨干旱区绿洲城市不透水面的变化规律，有利于城市合理发展规划，改善城市生态环境，具有重要的科学意义和参考价值。

2 研究区概况

本文以新疆维吾尔自治区（简称新疆）典型城市（乌鲁木齐、喀什、哈密、克拉玛依）为研究范围。乌鲁木齐地处亚洲大陆的地理中心，是新疆的人口集中区域，支撑着新疆的经济、社会和文化建设。但同时也伴随水资源匮乏、地表干燥，裸地广泛分布，生态环境脆弱等问题，给乌鲁木齐的社会经济发展产生较大影响（Fang，2019）。其中乌鲁木齐选取范围为43°41′50.446″N—43°56′37.376″N，87°24′55.192″E—87°42′53.392″E；喀什位于中亚腹部，雨水匮乏，在南疆地区的政治、经济、文化处于中心地位，研究范围为39°25′28.739″N—39°32′07.093″N，75°57′37.859″E—76°03′02.598″E；哈密处于新疆的东部，也是新疆东部经济建设的重要节点城市，近年哈密的城市规划建设发展迅速，不透水面的变化比较明显，城市社会经济文化稳步发展。研究范围为42°47′34.235″N—42°51′39.926″N，93°28′07.370″E—93°33′53.954″E；克拉玛依在新疆北部地区，冬夏差异悬殊，在春夏季大风较多，近年来，克拉玛依注重城市内部优化发展，提升人居适宜性，合理有效进行城市建设。研究范围为45°32′48.856″N—45°37′57.822″N，84°49′57.144″E—84°53′22.862″E；由于新疆幅员辽阔人口分布稀疏，所以重点选取了首府乌鲁木齐，喀什，哈密，克拉玛依作为本文的研究区，该研究区为人类活动集中区域，不透水面是其主要构成部分，如图1所示。

3 数据来源及预处理

获取了2017年、2019年4 个城市共8 景Sentinel-2A/B 影像数据，属于日间成像影像数据，选取成像日期相近数据，该影像数据云量覆盖低于5%，数据成像质量较好，可很好的用于不透水面研究，具体数据标识，获取时间详见下表1。

表1 本研究所采用的Sentinel-2A/B影像数据Table 1 Sentinel-2A/B image data used in this research

数据来源于欧洲航天局（ESA） 数据平台（https://scihub.copernicus.eu/dhus/#/home［2020-06-16］），下载获取的数据是大气表观反射率产品，数据已经过正射校正和几何精校正。本文采用ESA 提供的Sen2Cor 插件进行数据预处理，对Sentinel-2A/B 数据进行大气校正，设置命令行为10 m 分辨率，处理完成后导入SNAP 中转换为ENVI 标准格式，由于哨兵数据处理为10 m 分辨率数据后，丢失Coastal aerosol、Water vapour 和SWIR-Cirrus 3 个60 m 分辨率波段，剩余10 个波段，分别是B1-Blue，B2-Green，B3-Red，B4—B6、B7 为植被红边波段，B8-NIR，B9-SWIR1，B10-SWIR2。然后根据研究区范围裁剪输出数据，最后得到覆盖整个研究区范围的数据。

本文验证数据采用PlanetScope（PS）影像数据（https://www.planet.com［2020-06-16］），PS 卫星产品有4 个多光谱波段，波段范围为455—860 nm，空间分辨率达3 m，满足本文不透水面验证需求。

4 研究方法

4.1 Sentinel-2A/B影像L2范数归一化处理

在干旱区城市不透水面研究中，有较多的针对裸土的消除方法，其中余红娇（2018）在干旱区城市中采用裸土指数BSI（Bare Soil Index）进行裸土信息剔除，但仍然无法消除干旱区城市下垫面的裸地、沙地等光谱特征易混淆地物。夏宇等（2019）采用Sentinel SAR 和光学影像数据为主，辅以夜间灯光数据、DEM 数据为阈值筛选的提取方法，虽然能够提高精度，但不同数据间尺度转换对结果的影响以及提取方法过于复杂。L2 范数（Nie等，2010；Ren等，2012）在特征选择等机器学习领域有广泛应用，本文采用L2 范数对Sentinel-2A/B 影像进行归一化，对地物的光谱信息进行处理。其定义如下：

式中，xi是第i个像元的像元值，n是像元数目，X是图像矩阵。

如图2 所示，L2 范数归一化前，不透水面高值和低值光谱信息处于两个不同的区间，而裸土的光谱信息处于两者之间，且与不透水面信息混淆；L2 范数归一化后，不透水面和裸土的光谱信息空间被压缩，图2（b）中可见光波段（Blue、Red）裸土反射率与ISA 低值和高值分离明显；短波红外波段（SWIR1、SWIR2）裸土反射率与ISA低值分离明显，与ISA高值分离相对较弱。但四个波段处ISA 低值和ISA 高值都与裸土反射率差别明显，所以本文采用L2 范数归一化有助于减弱干旱区城市的裸土干扰信息，是下一步不透水面有效提取重要基础。

图2 Sentinel-2A/B L2范数归一化处理Fig.2 Sentinel-2A/B L2 norm normalization processing

4.2 增强型归一化差值不透水面指数（ENDISI）

本文在L2 范数归一化处理的影像基础上采用一种与其结合紧密的增强型归一化差值不透水面指数（ENDISI）进行不透水面提取。穆亚超等（2018）的增强型归一化差值不透水面指数（ENDISI） 选取了Blue、Red、NIR、SWIR1以 及SWIR2共5 个波段，该指数的不透水面增强因子根据反射率信息采用Blue和SWIR2波段，同时将Blue波段的影响因子扩大一倍；为了抑制非不透水面因子，选择Red、NIR、SWIR1波段来实现。最后，将以上波段综合进行差值运算，实现增强不透水面提取信息的结果。但单独用该指数进行干旱区城市不透水面提取并不能达到很好的效果。只有通过L2 范数归一化处理后，进一步分离不透水面与裸土光谱信息特征以凸显在Blue 和SWIR2波段二者的差异，使用ENDISI 指数进行干旱区城市不透水面提取才能够达到较好的提取结果。本文研究区范围内存在少量水体信息，通过McFeeters（1996）的归一化差异水体指数NDWI（Normalized Difference Water Index）对这一部分水体进行掩膜剔除。其中公式如下：

4.3 不透水面自适应阈值识别算法

为了检验上述方法的适用性，将其运用到新疆4个典型城市的不透水面提取，结果表明这一方法能够很好的提取不透水面信息，本文下一步将所提取的不透水面结果进行二值化处理，但由于新疆各城市间分布距离较远，针对不同区域的城市不透水面的分离阈值需调整，不能使用同一固定阈值，不同区域的城市不透水面阈值需根据情况率定，为了减少人为率定阈值造成的误差，本文采用最大类间方差法（OTSU）自动确定不透水面识别阈值（李琳琳，2012）。

最大类间方差法（Otsu，1979）是一种基于灰度图像的自适应阈值分割算法，该算法以图像灰度值分布特征将不透水面图像的灰度直方图分为前景和背景两部分，它的目的是使两部分之间的方差取得最大值，即分离性最大，使得不透水面主要信息展现在二值图上。设ENDISI 指数运算后的灰度图像等级为L，使得类间方差最大的灰度值即为所求的最佳阈值k。

式中，ω0和ω1为前景和背景出现的概率：

式中，ni为灰度值i出现的个数；N则是总的像元数目；k为划分前景和背景的值。

μ0和μ1为前景和背景的灰度均值：

式中，pi为灰度值i出现的概率。

μ为整幅图像的灰度均值：

4.4 精度评价方法

本文对不透水面提取结果的精度定量化评估中使用总体精度OA（Overall Accuracy）和Kappa系数（Kappa Coefficient）（赵英时，2003），其中公式如下：

总体精度OA 是具有概率意义的一个统计量，表述的是以本文方法提取的每一个不透水面随机像元pr与地面所对应区域的真实像元p相一致的概率。

Kappa 系数（Congalton，1991）也称为Khat统计，Khat是一种测定两幅图像之间吻合度或精度的指标。

式中，r是错误矩阵中总列数（即总的类别数）；xii是错误矩阵中第i行、第i列上像元数量（即正确提取的数目）；xi+和x+i分别是第i行和第i列的总像元数量；N是总的用于精度评估的像元数量。

4.5 不透水面空间差异分析

本文采用剖面线法和盒维数法对新疆典型城市2017年和2019年不透水面提取结果进行空间差异分析。

选取南北和东西方向剖面线，以城市不透水面中心区域为交点，沿剖面线对不透水面ENDISI指数值不同时间、不同区位的空间差异进行分析。剖面线法（傅滨桢，2017；李涵等，2019）能更细致的研究不透水面ENDISI 指数在空间形态的局地差异以及热点区域的扩展幅度。

盒维数（The Box-counting Dimension）也称计盒维数，是一种被学者广泛应用的分形维数（Orzechowski，1998；Miyata 和Watanabe，2003）。在分形理论应用（俞晓牮等，2019；林松等，2020）研究中被学者提出的大多数关于维数的定义都属于盒维数的一种变式，而且盒维数是基于相同形状集的覆盖而确定的，因此可以采用计算机技术来实现（杨书申和邵龙义，2006）。对于复杂的城市空间结构，本文拟采用计盒法来计算盒维数：首先，以边长为r的盒子将整个图像空间覆盖，由于图像信息分布不均，盒子可能为空，也可能仅仅涵盖了图像的部分信息。其次，统计涵盖图像信息的盒子数量N，然后改变盒子边长r的大小，来求得r与N相应的对数关系，以此来得到一系列数据对（-logr，logN）。最后，根据所得到的数据对，采用最小二乘法拟合，则该图像空间的盒维数以形成的直线的斜率来衡量（Sarkar 和Chaudhuri，1994；Li等，2009）。具体算法如下：

设A是Rn（n维欧氏空间）的任意非空有界子集，对于任意一个r＞0，N表示覆盖A所需要的边长为r的n维盒子的最小数目。如果存在一个数d，使得当r→0时，有

则

式中，k为常数，0＜r＜1；d为盒维数。

5 结果与分析

5.1 不透水面提取结果及精度评价

利用OTSU算法对ENDISI指数的灰度图像自动计算分类阈值。OTSU根据ENDISI指数的直方图自动确定使两类之间方差最大的阈值。根据计算结果，乌鲁木齐2019年和2017年阈值分别为0.0667，0.0784；喀什2019年和2017年阈值分别为0.0314，0.0353；哈密2019年和2017年年阈值分别为0.0353，0.0510；克拉玛依2019年和2017年阈值分别为0.0210，0.0363。图3 中黑色部分表示非不透水面，白色部分表示不透水面。最终获取了2019年和2017年的乌鲁木齐，喀什，哈密，克拉玛依主城区不透水面提取结果，为了对不透水面提取有一个准确的评估，以2019年10月乌鲁木齐不透水面提取结果为例，进行精度评价，利用ArcGIS随机点生成500 个验证点，参考同一时相的Planetscope卫星影像人工判别验证点是否为不透水面。然后在定量化的精度评估中使用总体精度和Kappa系数（Jensen 和Lulla，1987；Cheng 等，2020），总体精度为86.60%，Kappa系数为0.73。结果见表2。

表2 2019年10月乌鲁木齐不透水面提取结果精度评价Table 2 Accuracy evaluation of extraction results of impervious surface in Urumqi in October 2019

图3 研究区不透水面提取结果Fig.3 Results of impervious surface extraction in the study area

为了进一步展示本文研究方法在不同区域的提取精度，本文选取了2019年乌鲁木齐研究区的5 个典型局部区域（图4），图4（a）—（e）分别表示公园区、商业区、道路区、工厂区和裸土区，（f）—（j）表示对应局部区域不透水面提取结果。由图4（f）可看出，对于公园的内部道路提取效果较差，主要是因为公园内部道路宽度较小导致混合像元的存在；图4（g）在提取结果中，红色屋顶的光谱信息在红波段与近红外波段反射率较高，在蓝光波段与短波红外波段反射率相对较低，导致其与植被反射信息相似，采用ENDISI 指数法未能成功提取红色屋顶的不透水面信息；图4（i）中对于工厂区不透水面能够准确提取；图4（h）和4（j）中对植被和裸土有较好的剔除效果。

图4 局部区域不透水面提取结果Fig.4 Extraction result of impervious surface in local area

5.2 不透水面空间差异分析

5.2.1 不透水面剖面线ENDISI指数值分析

采用剖面线法对研究区不透水面各方向ENDISI 指数值空间变化分析。即从研究区中选取东西方向、南北方向上ENDISI 指数值进行分析，首先是确定东西方向、南北方向的交点，然后分别提取经过该点的南北、东西垂直剖面线上ENDISI 指数值，分析不同剖面线上ENDISI 指数值变化。如图5、6、7、8，其中图5（a）中横坐标Number of pixels 表示像元个数，即0—3000表示以城市北部边缘为起点向南部边缘沿剖面线扩展3000 个像元，纵坐标Data value 表示ENDISI 指数值，其中图5（b）中横坐标Number of pixels 表示像元个数，即0—2500表示以城市西部边缘为起点向东部边缘沿剖面线扩展2500 个像元，纵坐标Data value 表示ENDISI指数值。然后再以每个城市整体区域作为研究对象，分析剖面线上的ENDISI指数值的空间变化规律。

结合图5、6、7、8，对比新疆典型城市2017年、2019年ENDISI 指数值在剖面线上的变化特征可知：乌鲁木齐ENDISI指数值整体增加明显，其中北部和中西部增加较为显著；喀什ENDISI 指数值整体呈增长趋势，中部及北部区域增加较为明显；哈密北部和中东部ENDISI指数值增加显著，西部变化较小；克拉玛依ENDISI 指数值总体变化较小，北部和中西部区域均稍有增长。综上，乌鲁木齐由于建成区面积最大，城市构筑物结构复杂，所以ENDISI 指数值最高；克拉玛依建成区最小，城市结构相对简单，ENDISI 指数值最低；各城市ENDISI 指数值中心区域最高，向城市边缘逐渐减弱。

图5 乌鲁木齐2017年和2019年不透水面剖面线ENDISI指数值变化Fig.5 Impervious surface profile line ENDISI index value change in Urumqi 2017 and 2019

图6 喀什2017年和2019年不透水面剖面线ENDISI指数值变化Fig.6 Impervious surface profile line ENDISI index value change in Kashgar 2017 and 2019

5.2.2 不透水面盒维数分析

本文考虑到不透水面形状的复杂性，简单的统计学难以描述。因此，采用盒维数法分析不同时相不透水面的空间结构特征。采用Matlab计算乌鲁木齐、喀什、哈密、克拉玛依4 个城市2017年和2019年不透水面的盒维数值，结果如表3。由于新疆城市的空间结构复杂程度的不同，所得到的盒维数结果也不一致，通过计算得出新疆典型城市的盒维数值主要分布在1.731—1.8847。以不透水面二值化图（图3）来计算分形维数，并分析城市不透水面的空间结构变化规律。结合图3 和表3，可以明显的看出受人类活动干扰，随着城市的基础设施建设不断增加，2019年新疆典型城市的盒维数值较2017年都呈现增加趋势，城市空间结构复杂度不断增加。其中哈密盒维数变化幅度最大为3.05%，城市空间结构变化较为剧烈；克拉玛依盒维数变化幅度最小为1.41%，城市空间结构发展较为稳定。可以看出新疆典型城市不透水面的分形特征较为显著，由此可知，盒维数作为表征新疆典型城市不透水面随时间变化的特征参数是可行的。

表3 盒维数及变化幅度Table 3 Box dimension and change range

图7 哈密2017年和2019年不透水面剖面线ENDISI指数值变化Fig.7 Impervious surface profile line ENDISI index value change in Hami 2017 and 2019

图8 乌鲁木齐2017年和2019年不透水面剖面线ENDISI指数值变化Fig.8 Impervious surface profile line ENDISI index value change in Karamay 2017 and 2019

6 讨 论

目前不透水面提取的指数法较多，其中徐涵秋（2008）提出的NDISI指数最具有代表性，但干旱区城市周边高反射率裸土信息对这一指数干扰较大；由于采用的Sentinel-2A/B数据不具有热红外波段，也是无法选用NDISI 指数的一个重要原因。本文采用支持Sentinel-2A/B 计算的ENDISI 指数进行不透水面提取，但干旱区城市裸土对提取结果影响较大。因此，本文根据L2 范数归一化处理后的光谱信息，选择穆亚超等（2018）的增强型归一化差值不透水面指数（ENDISI）提取不透水面。以城市与裸土交界区为研究对象，对比ENDISI、NDBI 和UI 3 种方法的不透水面分类效果（图9），并采用目视解译法进行精度验证。结果表明ENDISI 提取效果较好，城市轮廓和裸土区道路提取精度最高（83.18%），NDBI 提取效果相对较差（78.38%），UI提取精度最低（60.18%）。

图9 3种方法不透水面提取结果对比Fig.9 Comparison of extraction results of impervious surface with three methods

以往对于不透水面空间差异的研究，普遍采用长时间序列30 m 的中等分辨率遥感数据对不透水面进行简单的分析，包括不透水面的面积和速率变化（包颖等，2020；张稼乐等，2020；Gong等，2020）。本文对不透水面进行剖面线分析，从城市的东西南北4 个方向分析ENDISI 指数值的变化特征，进而分析其在空间形态的局地差异及热点区域扩展幅度，以更好反映城市不透水面的发展状况。

如图10 所示，10 m 空间分辨率的Sentinel-2A/B 数据能够较好反映城市中的地物形状特征，而30 m 空间分辨率的Landsat 8 数据不能很好地刻画城市居住小区、商业建筑、城市道路等形状特征，用其进行盒维数分析意义不大，因此，本文采用10 m 空间分辨率的Sentinel-2A/B 数据进行不透水面空间结构差异分析，能够达到较好的研究效果。盒维数易于程序化计算，广泛应用于各类图形图像研究，受到了国内外学者的持续关注（Sarkar和Chaudhuri，1994；Li 等，2009；俞晓牮 等，2019；林松等，2020）。本文运用盒维数法反映城市空间结构的复杂程度，揭示不同时期城市空间结构的差异性，展示城市内涵式发展动态，利于决策者及时调整城市发展规划，促进城市内涵式与扩张式协调发展。

图10 乌鲁木齐高铁站区域Fig.10 Urumqi high-speed railway station area

受数据源所限，本文未能进行更长时间的不透水面空间差异规律研究。下一步应考虑融合多源数据进行长时间序列不透水面的差异性研究，结合社会人文因子、经济发展因子和自然环境因子等，多方位揭示城市不透水面变化规律，剖析城市不透水面动态变化的驱动机制，为干旱区城市生态环境保护提供理论依据。

7 结 论

本文通过L2 范数归一化处理结合ENDISI 提取了2017年和2019年新疆典型城市不透水面，并对提取结果进行精度评价和不透水面空间差异分析，得到如下结论：

（1）对Sentinel-2A/B 影像L2 范数归一化处理结合ENDISI 在减弱干旱区城市裸土干扰信息的同时能够有效实现不透水面的提取。

（2）利用OTSU 对不透水面提取阈值进行自动率定能够有效识别不透水面。以2019年乌鲁木齐为例，不透水面提取结果总体精度为86.60%，Kappa 系数为0.73，表明该算法能够很好的实现不透水面阈值率定。

（3）通过对不透水面空间差异分析可知：从剖面线角度分析得出乌鲁木齐北部区域ENDISI 指数值明显增加，喀什在中部和北部区域ENDISI 指数值增加明显，哈密则在北部和中东部ENDISI 指数值增加显著，克拉玛依在北部和中西部区域ENDISI 指数值增加较少；从不透水面盒维数分析中得出，新疆典型城市的盒维数值都呈现增加趋势，城市空间结构复杂度不断增强，其中，哈密盒维数值最大，乌鲁木齐盒维数值最低，而哈密的盒维数值变化幅度最大，克拉玛依的盒维数值变化最小。