半城市化地区不透水面遥感提取优化方法研究

崔 雯 卓，林 晨，陈 曦，熊 俊 峰，侯 渲，冯 雪 娇，吴 紫 静,5

(1.中国科学院流域地理学重点实验室，中国科学院南京地理与湖泊研究所，江苏 南京 210008；2.中国科学院大学，北京 100049；3.中山大学地理科学与规划学院，广东 广州 510275；4.自然资源部海岸带开发与保护重点实验室，江苏 南京 210023；5.南京工业大学测绘科学与技术学院，江苏 南京 211816)

0 引言

中国已进入快速城市化阶段，大量具有城市化基础的半城市化地区逐渐转为城市地区[1-4]。不透水面(Impervious Surface Area,ISA)是典型的人工地表，目前已成为研究城市化进程及其发展强度最直观的指标[5，6]。与城市化地区相似，半城市化地区主要由不透水面覆盖[7-9]，但其建筑密度和功能与城市化地区存在明显差距。半城市化地区作为城市发展的后备用地，肩负着承接城市职能和支持城市发展的重任[10]，因此，快速获取半城市化地区的不透水面分布信息对城市规划和管理具有重要意义。

遥感技术因具有探测范围大、获取资料速度快、实时性强等特点，被广泛应用于不透水面提取[11]。目前，不透水面提取方法分为两大类：1)基于输入图像的分类方法[12]，包含基于像素或面向对象的分类方法[13-16]，如光谱混合分析法(Spectral Mixture Analysis,SMA[17])、影像分类法[18,19]、基于人工神经网络(Artificial Neural Network,ANN)[20]的分类方法等；2)通过阈值直接分割图像的方法(指数法)，即利用不透水面的光谱影像特征，通过构建遥感光谱指数提取不透水面信息[11,21]，常采用影像的原始波段构建不透水面指数，通过对不透水面与其他地类之间有明显区别的反射波段进行运算，以此抑制其他地类信息并增强不透水面信息，该方法无需预设训练样本，可操作性更强且自动化程度更高。目前，应用最广泛的不透水面提取指数是归一化建筑指数(Normalized Difference Build-up Index,NDBI)[22]，该指数初次应用于提取南京市不透水面且总体精度高达92.6%。由于城市地表具有复杂异质性，原始波段的组合无法精确地剔除不透水面以外的地物，一些学者采用其他指数组合或利用分量方法重构不透水面指数。例如：Sun等[23]提出组合建筑区指数(Combination Build-up Index,CBI)，通过组合主成分分析的第一分量和归一化差异水体指数(Normalized Difference Water Index,NDWI)[24]提取不透水面信息；Xu[25]根据改进型归一化差异水体指数(Modified Normalized Difference Water Index,MNDWI)、土壤调节植被指数(Soil-Adjusted Vegetation Index，SAVI)[26]和归一化建筑指数(NDBI)提出遥感建筑用地指数(Index-based Built-up Index,IBI)，并在福州地区得到有效应用；Bhatti等[12]结合归一化差异水体指数(NDWI)、归一化植被指数(Normalized Difference Vegetation Index,NDVI)[27]和改进型归一化差异水体指数(MNDWI)[28]构建了建筑指数(Built-up Area Extraction Method,BAEM)。

综上，现有研究针对不透水面提取构建了较为成熟的指数，并在京津冀[29-31]、长江三角洲[32,33]等高度城市化地区得到广泛应用，却极少应用于半城市化地区。半城市化地区建筑密度较低，不透水面周围存在大量耕地与裸土，不仅需要突出不透水面光谱信息，还需要抑制其他地物的光谱信息，而现有不透水面提取指数难以适用于半城市化地区，因此，本研究提出一种适用于半城市化地区的不透水面提取方法。

1 研究区概况与数据来源

1.1 研究区概况

目前我国城市正处于转型期，江苏省提出“1+3”功能区战略助推高质量发展，江淮生态经济区(图1a)是全省经济发展相对落后的区域，该区域城镇化水平差异较大，城镇体系结构不完整，城镇综合承载力低，亟须协同推进区域新型城镇化建设，是城市转型需要重视的区域。本文以区县级行政单元作为半城市化地区，选取江淮生态经济区的淮阴区和盱眙县为研究区。1)盱眙县位于江苏省中西部(图1b)，地处淮河中下游、洪泽湖南岸，作为丘陵区的代表，耕地多为水田，以坑塘养殖产业为主，区域经济发展状况位居江苏省中等水平，根据第七次全国人口普查数据，盱眙县城镇化率达61.66%，近5年提高近10.46%，处于逐渐向城市转型时期。2)淮阴区位于江苏省北部平原的中心(图1c)，耕地以旱地为主，处在江淮生态经济区、淮河生态经济带、大运河文化带战略交汇处，又位于长三角一体化发展、长江经济带和南京都市圈的影响范围内，目前经济发展位居江苏省中等偏高水平，根据第七次全国人口普查数据，淮阴区城镇化率达57.61%，近5年提高了12.95%。为更形象地体现半城市化地区在遥感影像上的特点，选取A、B、C、D 4个区域(图1d)进行展示，其中A和C地区除不透水面外还存在大量裸土，B和D地区建筑密度较低，周围被大量耕地包围。

图1 研究区概况Fig.1 Overview of the study area

1.2 数据及预处理

研究区存在大量较窄道路，选择空间分辨率较高的遥感数据可以尽量减少混分情况。本文根据研究区实际，选用云量小于5%、空间分辨率为10 m的Sentinel-2A数据，数据来源于哥白尼数据开放访问中心(Copernicus Open Access Hub)，获取时间均为植被茂密、坑塘未排水的2021年7月30日。研究采用的波段为可见光波段(Band 2、Band 3、Band 4，10 m空间分辨率)、近红外波段(Band 8，10 m空间分辨率)和短波红外波段(Band 11、Band 12，20 m空间分辨率)。利用欧空局(ESA)发布的Sen2Cor工具箱对Sentinel-2A 1C级产品进行辐射定标、大气校正，最后将数据重采样至10 m。

2 半城市化地区不透水面提取方法

传统BAEM指数[12]内部的约束条件并不符合半城市化地区情况，应用于半城市化地区时，不透水面提取效果较差。因此，本文从以下两方面进行优化：1)考虑半城市化地区存在大量裸土，需要对现有的指数增加土壤调节植被指数(SAVI)进行优化，以区分裸土的影响；2)传统方法采用双窗口变步长阈值搜寻算法(DFPS)[26]提取阈值，考虑到半城市化地区存在混合像元的问题，对样本选取要求方面进行阈值法改进。

2.1 传统BAEM指数优化

BAEM指数将地物划分为水体、植被与不透水面，利用各地物最具代表性的指数(归一化建筑指数(NDBI)、归一化植被指数(NDVI)和改进的归一化差异水体指数(MNDWI))进行差值运算(式(1)-式(4))，以有效剔除不透水面以外的地物。然而，将传统BAEM指数方法应用于半城市化地区时，易导致不透水面提取结果存在大量裸土与不透水面混分。

BAEM=NDBI-NDVI-MNDWI

(1)

(2)

(3)

(4)

式中：ρNIR、ρRed、ρGreen、ρSWIR分别为近红外(Band 8)、红光(Band 4)、绿光(Band 3)、短波红外波段(Band 11)的地表反射率。

本研究将地表覆盖物分为不透水面、水体、植被和裸土4类典型地物，由其光谱曲线(图2)可以看出，不透水面与裸土的光谱特征相似，从绿光波段(Band 3)到短波红外波段(Band 11)呈递增趋势，到Band 12反射率降低，植被呈明显的红边效应，水体整体反射率较低。因半城市化地区存在大量裸土，需要在传统的BAEM指数中增加SAVI指数(式(5))进行优化，以区分裸土的影响，其中，土壤调节因子L取0.5，最终得到改良后的BAEM指数(MBAEM)(式(6))。

图2 典型地物反射率光谱曲线Fig.2 Reflective spectrum curves of typical surface features

(5)

MBAEM=NDBI-SAVI-NDVI-MNDWI

(6)

为验证加入SAVI能否区分裸土与其他典型地物，选取盱眙县和淮阴区具有代表性的典型像元进行分析。在每个区域中选择较密集的建成区a、被大量耕地包围的房屋b、道路c、水体d、植被e及裸土f，分别展示各个像元的光谱曲线以及4种指数(SAVI、NDVI、NDBI、MNDWI) (图3-图5)。由图3b和图4b可知，a、b、c 3类不透水面具有相似的光谱曲线特征，但其光谱曲线与水体d、植被e差别较大，与裸土f光谱曲线趋势一致，因此将其像元的4种指数进行对比。如图5所示，a、b、c 3类不透水面的4种指数变化趋势基本一致；水体d的MNDWI指数与其他指数区分明显；植被e的NDVI指数与其他指数区分明显；裸土f的光谱曲线与不透水面相似，但在指数方面存在一定差异，因此可以利用4种指数的组合抑制植被、水体与裸土的信息，突出不透水面信息。

图3 盱眙县典型地物的像元光谱曲线Fig.3 Pixel spectrum curves of typical surface features in Xuyi County

图4 淮阴区典型地物的像元光谱曲线Fig.4 Pixel spectrum curves of typical surface features in Huaiyin District

图5 盱眙县和淮阴区典型地物像元的指数特征Fig.5 Pixel exponential characteristics of typical surface features in Xuyi County and Huaiyin District

2.2 阈值法改进

Chen等提出的双窗口变步长阈值搜寻算法(DFPS)是在土地利用/覆被变化监测中提取变化像元的半自动阈值搜寻算法[34]。首先选取典型训练样区的样本，其中选取的样本最好为“岛状”，即典型变化区被非典型变化像元包围，使典型区与非典型区存在明显的差异性；其次设定搜寻范围和步长，根据给定的步长在合理的搜索范围内寻找获得典型样区变化检测精度最高的阈值。搜索过程需要进行反复尝试，每次循环对应的范围和步长会依次缩小，保证搜索的变化阈值最优。该算法存在的问题是训练样本需要被其他地类训练样本包围，样本之间相互依赖，需要形成非独立的双窗口。但在半城市化地区，训练样本往往复杂交错，仍存在混合像元问题，即使10 m空间分辨率也不能保证样本完全独立，因此本研究采用大津法(OTSU算法)[35]提取阈值，使图像分割不受样本影响，计算相对简单。

OTSU算法[35]按图像的灰度特性分为背景和前景两部分，自动计算背景与前景的最大类间方差，从而将图像进行二值化，其最大优势在于不需要预设样本，自动获取最佳阈值，计算简单。本研究以MBAEM指数灰度图像为例对该算法进行展开介绍。理论上MBAEM图像的灰度值范围为-3～3，其中前景(不透水面)的灰度值大于背景(植被、水体以及裸土)的灰度值,利用OTSU算法分离不透水面和其他地类的计算公式见式(7)。

(7)

式中：σ2为前景和背景图像的方差；μ1和μ0分别为前景和背景的平均值；μ为整幅影像像元平均值；w1和w0分别为前景和背景在整幅影像中的占比。

根据式(7)提取的阈值threshold，对MBAEM指数灰度图像进行分割，得到半城市化地区不透水面影像(Impervious Surface,IS)(式(8))。本研究中MBAEM指数提取不透水面在盱眙县阈值为-0.48，淮阴区阈值为-0.45，BAEM指数提取不透水面在盱眙县和淮阴区的阈值分别为-0.28和-0.33。

(8)

式中：(i，j)表示MBAEM指数影像的像元行列号。

2.3 精度检验

为客观验证不透水面提取结果的准确性，本研究在Sentinel-2A原始影像上，采用目视解译方法和分层随机抽样方法选取验证样本[36]，并在Google Earth上选取与成像时间接近的影像对验证点进行比对和微调。最终，在淮阴区选取不透水面与非不透水面样本数量分别为3 964个和4 272个，在盱眙县选取不透水面与非不透水面数量分别为3 261个和3 796个。

3 不透水面提取结果与评价

3.1 优化方法与其他方法提取结果对比

3.1.1 优化方法与传统方法提取结果对比 基于Sentinel-2A数据，采用BAEM指数和MBAEM指数提取盱眙县和淮阴区的不透水面，分别用两种阈值法提取阈值(图6、图7)。从空间分布上看，盱眙县的不透水面主要集中在中部建成区，其他村镇向四周呈现放射状分布，村镇之间有道路连接；淮阴区不透水面也主要位于中部建成区，周围分布村镇，道路纵横交错。从提取结果整体上看，BAEM方法得到的不透水面提取结果分布更多，因此，本研究进一步采用总体精度(OA)、Kappa系数、制图精度 (PA)和用户精度 (UA)进行精度检验(表1)，以客观、具体展示4种方法的不透水面提取结果。由表1可知：盱眙县和淮阴区利用MBAEM指数和OTSU算法得到的结果均为最佳，其总体精度达89.50%和89.73%，Kappa系数为0.7907和0.7950；相比BAEM指数，MBAEM指数提取的总体精度提高约10%，Kappa系数提高约15%。综上，MBAEM指数提取结果比BAEM指数提取结果精度更高，OTSU算法的精度比DFPS算法的精度更高，因此本研究基于指数和阈值两方面进行优化的精度均得到提升。

图6 盱眙县不透水面提取结果对比Fig.6 Comparison of extraction results of impervious surface in Xuyi County

图7 淮阴区不透水面提取结果对比Fig.7 Comparison of extraction results of impervious surface in Huaiyin District

表1 不透水面精度验证结果Table 1 Accuracy validation of impervious surface extraction results

3.1.2 MBAEM指数与其他常用指数提取结果对比 在研究区选取A、B、C、D 4个区域，对NDBI、IBI和CBI指数提取结果与MBAEM指数提取结果进行对比(图8)，所有指数阈值均采用OTSU法确定。可以看出，MBAEM指数的提取结果更接近真实情况；IBI指数的提取结果分布相对较少且稀疏，存在大量漏提；CBI指数的提取结果更密集，但也有许多非不透水面被错分为不透水面；在细节上，NDBI指数则漏分或错分了一些半城市化地区中的典型地物(如裸地)。由盱眙县和淮阴区的提取精度(表2)可知：MBAEM指数的总体精度、Kappa系数、制图精度和用户精度均最高。

图8 MBAEM与3种指数不透水面提取结果对比Fig.8 Comparison of extraction results of impervious surfaces using MBAEM and three indexes

表2 MBAEM指数与其他3个指数提取结果精度指标对比Table 2 Accuracy comparison of impervious surface extraction results using MBAEM and three indexes

3.2 优化方法与传统方法提取结果细节对比

通过对比发现，阈值上OTSU算法得到的结果优于DFPS算法的结果，因此选取OTSU算法的结果进行细节对比。将原始BAEM指数提取结果与优化后的MBAEM指数提取结果局部进行对比(图9、图10)，发现优化后的指数提取结果在水体和裸土方面效果更好。1)水体(图9)。水体在真彩色影像呈蓝色系线条，而沥青混凝土等道路在真彩色影像上呈灰色条状特征，因此二者呈现的特征相似。图9A中，道路呈浅色线状分布，河流呈深蓝色，MBAEM指数可以准确提取道路，但BAEM指数将水体误分为道路；图9B与图9C中，MBAEM指数均将水体与不透水面更好地区分开，而BAEM指数却将水体误分为不透水面。因此，MBAEM指数较好地解决了水体被误分为不透水面的问题且能清晰提取出道路。2)裸土(图10)。裸土在真彩色合成影像中呈现棕色系，大部分呈矩形块状分布，少量呈条状零星分布，且颜色偏亮。传统提取方法未考虑裸土与不透水面各波段均具有较高反射率的特点，因此未将二者区分开。如图10A所示，MBAEM指数较好地区分开建设用地周围未开发利用的裸土区域，而BAEM指数却将裸土与不透水面混分；在图10B和图10C中，MBAEM指数能有效将耕地周围未耕种的裸土与不透水面区分开。由此可见，通过增加土壤调节植被指数(SAVI)，能较准确地提取半城市化地区不透水面。

图9 水体的局部对比Fig.9 Local contrast for water

图10 裸土的局部对比Fig.10 Local contrast for bare soil

4 结论与展望

本文以盱眙县和淮阴区作为半城市化地区的典型区，利用Sentinel-2A遥感数据，在原有BAEM指数的基础上优化指数选择和阈值方法，提出一种新的不透水面综合指数MBAEM，该指数综合考虑了不透水面与水体、裸土和植被之间的差异，并采用OTSU阈值法提取半城市化地区的不透水面。结论如下：1)利用MBAEM指数和OTSU阈值法提取效果最好，盱眙县和淮阴区的总体精度分别为89.50%和89.73%，Kappa系数分别为0.7907和0.7950； 2) 相较于传统BAEM指数提取方法，MBAEM指数较好解决了部分不透水面与水体、裸土混分的问题；相比其他提取不透水面指数，总体精度和Kappa系数均有不同程度的提升，且优化方法的提取结果在空间分布上更接近真实结果。

本研究提出的MBAEM指数可以有效识别遥感影像下具有复杂环境的半城市化地区的不透水面，部分解决了传统指数裸土与不透水面混分的问题，且能挖掘不透水面的特征信息，进而较好地获取半城市化地区不透水面。但本研究仍存在局限性：1)当选用不同的影像时，利用OTSU算法得到的每期影像阈值不同，存在不确定性，未来需综合考虑多区域多时序，以得到一个相对具有代表性的阈值；2)该方法目前仅在小区域得到实证，在更大区域尺度的应用有待进一步研究。

感谢国家重点研发计划(2021YFD1700600)项目为本研究提供资金支持；感谢国家科技资源共享服务平台-国家地球系统科学数据中心湖泊-流域分中心(http://lake.geodata.cn)提供数据支持！