Landsat8不透水面遥感信息提取方法对比

刘畅， 杨康,2,3， 程亮,2,3， 李满春,2,3， 郭紫燕

(1.南京大学地理与海洋科学学院，南京 210023； 2.江苏省地理信息技术重点实验室，南京 210023； 3.中国南海研究协同创新中心，南京 210023)

0 引言

伴随着全球范围内的城市化，不透水面正快速取代自然地表，目前已成为一种十分关键的地表覆盖类型[1-2]。不透水面是指水体不能下渗至土壤的物质，包括自然不透水面和人工不透水面[3]。不透水面的时空分布，直接反映了城市扩张趋势[4-5]。不透水面阻止雨水下渗，导致地表径流量增大，引发城市内涝，威胁地下水资源补给[6-7]。同时，不透水面会作用于城市热岛效应，直接影响城市宜居性[8-10]。因此，准确提取不透水面，尤其是人工不透水面，对于及时监测城市动态具有重要意义。

不透水面遥感信息提取方法主要包括指数法[11-12]、分类回归树法[13]、支持向量机法[14-15]和光谱混合分析法[16-17]等。其中，指数法可操作性强，自动化程度较高，提取结果不受训练样本影响，能够快速提取不透水面遥感信息，现已被广泛应用[18-24]。各种不透水面指数设计的基本原理相似，均是通过波段组合增强不透水面与其他地表覆盖类型的差异，但具体使用的波段存在差异。例如，陈洁丽等[22]利用红光、近红外和短波红外波段设计了新建筑指数(new built-up index，NBI)，从Landsat5影像中提取了常州市不透水面信息，总体精度达到90%； 田玉刚等[23]利用蓝光和近红外波段设计了垂直不透水层指数(perpendicular impervious index，PII)，从Landsat8影像中提取了武汉市与北京市的不透水面遥感信息，总体精度高于96%。

由于使用的波段不同，不同指数提取不透水面遥感信息的结果存在显著差异，这对综合分析不透水面动态变化带来了较大不确定性[25]。同时，大部分不透水面指数只针对不透水面密集分布的城市区域[18,20]，在大区域应用的能力有待验证。对比分析不透水面指数的提取效果，是不透水面遥感信息提取研究的重要内容。为此，本研究归纳总结了现有的不透水面指数，选择具有不同地表覆盖类型特点的实验区，分析评价不透水面指数的提取精度，揭示现有不透水面指数的设计优势和存在问题。由于现有不透水面指数主要用于人工不透水面(建筑物和道路等)提取，进而监测城市动态变化，因此，本研究仅考虑人工不透水面，评价不透水面指数提取人工不透水面的精度，不考虑裸岩等自然不透水面。

1 不透水面指数对比分析

1.1 不透水面指数

本研究归纳总结了目前较常用的8种不透水面指数，按照其设计方法的不同，将其划分为3类。

1)利用卫星原始波段构建的不透水面指数。此类不透水面指数包括城市指数(urban index，UI)[19]、NBI[22]、归一化建筑指数(normalized difference built-up index，NDBI)[21]、PII[23]和比值居民地指数(ratio resident-area index，RRI)[24]。这类不透水面指数都是利用不透水面区别于其他地类的最强和最弱反射波段设计，以此增强不透水面遥感信息，抑制其他地类信息。图1展示了不同地类在Landsat8影像不同波段的大气表观反射率。

图1 Landsat8卫星遥感影像不同地类光谱特征Fig.1 Spectral characteristics of different land covertypes in Landsat8 satellite image

UI和NDBI通过选取近红外波段和短波红外波段设计，RRI和PII则利用了蓝光波段和近红外波段。虽然RRI和PII具有相同的光谱波段，但在波段组合形式上分别采用了比值组合型和线性组合型。此外，在波段权重的设置上，RRI中的蓝光波段和近红外波段具有相同权重，而PII能够根据实验区光谱自适应调整蓝光波段系数m、近红外波段系数n及常数C[23]。NBI选取了红光、近红外和短波红外波段，利用乘法扩大特征波段灰度值，进而增强不透水面遥感信息。

2)利用复合波段构建的不透水面指数。此类不透水面指数的代表是归一化差值不透水面指数(normalized difference impervious surface index，NDISI)[12]。这是第一个直接针对不透水面特性建立的不透水面指数，它组合了可见光、近红外、短波红外波段及热红外波段数据提取不透水面。其中，可见光波段包括蓝光、绿光、红光波段。NDISI还能够根据实验区的实际情况采用改进的归一化差值水体指数(modified normalized difference water index，MNDWI)[26]代替可见光波段。在使用NDISI提取不透水面遥感信息时，需要分别测试可见光波段和MNDWI指数，以获得更好的不透水面提取结果。

3)利用指数波段或分量构建的不透水面指数。此类不透水面指数包括基于指数的建筑指数(index-based built-up index，IBI)[20]、生物物理组分指数(biophysical composition index，BCI)[18]。Xu[20]基于城市中的地表覆盖类型[21]，以NDBI指数、土壤调节植被指数(soil adjustment vegetation index，SAVI)[27]和MNDWI指数分别代表不透水面、植被和水体，利用这些指数代替传统的卫星原始波段设计了IBI。Deng等[18]则利用缨帽变换的亮度(TC1)、绿度(TC2)、湿度(TC3)分量分别代表不透水面、植被和水体，设计了BCI提取不透水面遥感信息。

不透水面指数统计情况如表1所示。

表1 不透水面指数汇总Tab.1 Summary of impervious surface indexes

①公式中BLUE，RED，NIR，SWIR1，SWIR2和TIR分别为蓝光、红光、近红外、短波红外1、短波红外2和热红外波段的像素值；VISIBLE为可见光波段的像素值。

1.2 实验数据

本研究使用Landsat8卫星遥感影像作为实验数据。该卫星于2013年发射，重访周期为16 d，包含陆地成像仪(operational land imager,OLI)和热红外传感器(thermal infrared sensor,TIRS)传感器，能够提供30 m空间分辨率影像，目前已被广泛应用于地表覆盖动态变化分析[28]。由于大气校正对指数计算的结果影响有限，目前已被广泛使用的MNDWI水体指数、BCI和NDISI等不透水面指数都没有进行大气校正[18,26,29]，因此本研究使用Landsat8卫星原始影像，没有对其进行大气校正处理。为了对比不同指数采用Landsat8影像提取不透水面的效果，研究使用1 m空间分辨率的不透水面产品(http: //www.chesapeake. org)作为不透水面真值数据，验证不同指数不透水面提取精度。该数据产品为2013年美国东北部切萨皮克湾地区的地表覆盖数据，融合美国1 m空间分辨率的NAIP影像和LiDAR点云数据制作而成，是全球空间分辨率最高的土地覆盖数据集之一。本研究从该地表覆盖数据产品中提取出不透水面类型，包含建筑物和道路等2类人工不透水面，获取了高空间分辨率不透水面数据。从官方提供的数据产品精度评价和Google Earth影像目视解译分析可知，该不透水面数据产品分类精度较高，因此，本研究将其作为真值数据用于验证不透水面指数提取精度。为了与Landsat8影像的空间分辨率保持一致，利用众数法将1 m不透水面真值数据重采样至30 m。

1.3 指数对比分析方法

不透水面指数能够增强Landsat8影像不透水面遥感信息，但为了进一步区分不透水面与非不透水面，仍需要设定合理阈值生成不透水面二值图像，获取“不透水面”与“非不透水面”类别信息。通过比较二值图像与不透水面真值数据，一个像元将会出现4种可能结果： 真正(true positive，TP)、真负(true negative，TN)、假正(false positive，FP)、假负(false negative，FN)。TP或TN分别代表二值图像中被正确分为“不透水面”或“非不透水面”的像元， FP代表被错误判断成“不透水面”的“非不透水面”像元，FN代表被错误判断成“非不透水面”的“不透水面”像元。在此基础上，计算真正率(true positive rate，TPR)和假正率(false positive rate，FPR)能够评价不同不透水面指数的提取精度。TPR表示被正确分为“不透水面”像元(TP)占真值为“不透水面”像元(TP+FN)的比例，反映指数提取不透水面的正确率。FPR表示被错误分为“不透水面”中的“非不透水面”像元(FP)占真值为“非不透水面”像元(FP+TN)的比例，反映指数提取不透水面的错误率。

TPR的计算公式为

(1)

FPR的计算公式为

(2)

此外，总体精度(overall accuracy，OA)表示被正确分类的像元(TP+TN)占全部像元的比例，这一指标能够反映不透水面提取结果的总体情况，其计算公式为

(3)

本研究绘制接收者操作特征曲线(receiver operating characteristic curve，ROC)以比较不透水面指数在不同全局阈值下的提取精度。由于在生成不透水面二值图像时，阈值选择会对提取精度产生重要影响[21,30-31]。如果设定的阈值过低，会高估地表不透水面真实分布，在提取结果中混入非不透水面信息。如果阈值过高，部分不透水面不能被有效提取，将会低估地表不透水面真实分布。ROC曲线能够定量评估分类方法在选择不同阈值下的表现优劣，现已广泛应用于机器学习和数据挖掘领域[32]。ROC曲线以FPR为横轴，TPR为纵轴，是不同阈值下对应TPR和FPR点的连线。如果ROC曲线越接近坐标左上角((0，1)点)，即TPR越高、FPR越低，表示不透水面指数的提取精度越高。如果某一不透水面指数的ROC曲线能够完全包括其他指数的ROC曲线，可判断该指数的不透水面提取表现优于其他指数。若出现不同指数的ROC曲线互相交叉的情况，还需计算ROC曲线下面积(area under curve，AUC)进一步判断不透水面指数的优劣。AUC是ROC曲线下的区域面积，通常取值在0.5～1之间。AUC越大，不透水面指数的提取精度越高。

研究从ROC曲线中可以获取不透水面提取的最优全局阈值[33]。绘制ROC曲线时，最初设定的阈值较低，此时TPR与FPR值均较大。随着阈值不断提高，TPR与FPR值会同时减小。因此，最优全局阈值的确定需要权衡TPR与FPR。本研究使用Youden Index(J)方法[34]确定不透水面提取的最优全局阈值。该方法中P值能够评价二值分类的表现，其计算公式为

P=TPR-FPR。

(4)

P值越高，表明二值分类效果越好。当P值最大时，分类效果最好，以该点的阈值划分不透水面遥感信息，能够在降低错误分类的同时准确提取不透水面，进而获得不透水面指数的最优提取效果。

2 实验结果

2.1 方法实现

在Visual Studio平台下，利用Python2.7开发环境实现前文所述的8种不透水面指数提取方法。在计算不透水面指数时，大多研究选择了原始影像的灰度值，如NDBI，RRI，NDISI，IBI和NBI指数，仅有PII指数采用表观反射率计算。因此本研究统一采用原始影像的灰度值计算不透水面指数。根据不同指数要求，对RRI，BCI和PII覆盖MNDWI水体掩模。设置SAVI指数的经验常数为0.5[22,27]，BCI中所使用的Landsat8缨帽变换系数由Baig等[35]提供，设置PII中常数C=0.02。由于PII在不同实验区的蓝光波段系数m和近红外波段系数n非常接近，研究分别取田玉刚等[23]所选实验区m和n的平均值(m=0.905，n=0.435)。通过分别测试NDISI的可见光波段和MNDWI指数，当以MNDWI指数作为输入参数时，实验区的不透水面提取精度最高，因此本研究只选用MNDWI计算NDISI。此外，对不同不透水面指数的提取结果进行归一化处理，设ROC曲线的阈值间隔为0.01。

2.2 实验区

本研究选择美国华盛顿(实验区1)和宾夕法尼亚州东部(实验区2)作为实验区，分别获取了2个实验区的Landsat8影像作为输入数据(图2)。实验区1影像获取时间为2013年4月21日，实验区2影像获取时间为2014年5月26日，图2(b)和(d)均采用B7(R)，B6(G)，B4(B)假彩色波段组合，图2(c)和(e)为高空间分辨率不透水面真值数据，其中红色代表不透水面。2个实验区的面积较大，适合测试不同不透水面指数在大区域的应用表现。同时，2个实验区具有不同的地表覆盖特点，能够反映不透水面指数在不同地表覆盖条件下的提取效果，具有较强的代表性。实验区1的地表覆盖类型以植被、人工不透水面和水体为主，裸地面积较小。研究将裸岩与未种植农作物地表裸露的农田均视为裸地。该区的不透水面分布特点多样，中心城区的不透水面密集、面积大，郊区的不透水面稀疏、面积小。实验区2的主要地表覆盖类型是植被，裸地面积较大，而人工不透水面分布稀疏、面积小。实验区1和实验区2的面积分别为2 464像素×2 225像素和1 366像素×1 341像素。同时期获取并经重采样后的高空间分辨率不透水面真值数据作为精度验证数据。

(a) 实验区区位图 (b) 实验区1遥感影像 (c) 实验区1真值数据

(d) 实验区2遥感影像 (e) 实验区2真值数据

图2 实验区分布及Landsat8影像

Fig.2StudyareasdistributionandLandsat8satelliteimages

2.3 不透水面提取精度评价

实验获得了8种不透水面指数在2个实验区的ROC曲线(图3)。结果表明，PII在2个实验区所得的AUC值均最高，分别为0.895和0.839，说明PII提取不透水面遥感信息的精度最高； 其次是BCI和RRI，对应AUC值略低于PII； 而NDBI，IBI与NDISI这3个指数提取不透水面的精度较低，对应AUC值均低于0.8(表2)。此外，对比实验区1和实验区2，实验区1的不透水面提取精度高于实验区2，PII的AUC值从0.895降低至0.839，RRI的AUC值降低了0.050。这说明，相比裸地分布较广的区域，现有不透水面指数更适用于不透水面分布较为密集的区域。

(a) 实验区1 (b) 实验区2

图3 不透水面指数的ROC曲线 Fig.3 ROC curves of impervious surface indexes表2 不透水面指数在实验区1和实验区2的AUCTab.2 AUC of impervious surface indexes in study areas 1 and 2

为了进一步对比不同指数在各自最优阈值条件下的最佳不透水面提取结果，综合权衡TPR和FPR，利用Youden Index(J)方法确定了不同不透水面指数的最优全局阈值(图3)，并以该阈值划分不透水面和非不透水面，获取了2个实验区的不透水面二值图像(图4和图5)，其中白色代表不透水面。

(a) 真值数据 (b) PII(c) BCI

(d) RRI (e) NDBI(f) UI

(g) IBI (h) NDISI (i) NBI

图4 实验区1中最优不透水面提取效果

Fig.4Optimalimpervioussurfaceextractionresultsofstudyarea1

(a) 真值数据 (b) PII(c) BCI

(d) RRI (e) NDBI(f) UI

(g) IBI (h) NDISI (i) NBI

图5 实验区2中最优不透水面提取效果

Fig.5Optimalimpervioussurfaceextractionresultsofstudyarea2

在实验区1，PII，BCI和RRI的最优不透水提取结果与不透水面真值数据的吻合度较高，所对应的不透水面提取OA分别是89.6%，87.4%和87.5%，高于其他指数的OA(表3)。在这3种指数中，PII获得了最高的不透水面提取总体精度，提取不透水面的效果最好，原因在于相比BCI和RRI，PII的FPR值更低(9.5%)，表明PII提取不透水面的错误率更低。这与PII的设计原理有关，尽管PII和RRI的输入波段相同，但PII采用了波段线性组合形式，相比于RRI的波段比值组合形式，PII能够根据实验区的土壤线调整波段权重，减弱了土壤噪声影响[23]。而NDBI，IBI和NDISI的最优不透水面提取效果较差，对应的OA均低于75%。它们表现较差的主要原因是FPR值较高。IBI和NDBI的FPR值分别为28.2%和26.7%，远高于PII的FPR值。在实验区2，尽管PII，BCI和RRI的不透水面提取OA仍然高于其他5种不透水面指数，分别是77.3%，77.4%，78.3%，但相比于实验区1的实验结果，它们的不透水面提取精度均明显降低。PII指数的TPR值降低到77.1%，而FPR值升高到22.7%。除了PII，BCI和RRI，其他5种不透水面指数的OA也明显下降，平均降幅约9.6%。

表3 最优阈值不透水面指数的精度评价结果Tab.3 Accuracy evaluation of impervious surface indexes using the optimal threshold (%)

实验测试的8种指数提取出的不透水面范围，均大于地表真实不透水面分布，但也漏提了部分不透水面。在实验区1，UI，NDBI，IBI和NDISI将部分水体误提取为不透水面，RRI和IBI则漏提了郊区面积较小的不透水面。而在实验区2，所有指数均未能准确提取人工不透水面信息，提取结果中混入了大面积裸岩等自然不透水面，该结果表明不透水面指数难以正确区分人工和自然不透水面，这将影响不透水面指数监测大区域城市动态变化的能力。造成不透水面指数误提的主要原因为： 一方面，部分指数(UI，NDBI，IBI和NDISI)在设计时，低估了水体对不透水面的干扰作用，未作掩模处理[36-37]； 另一方面，由于不透水面建设材料的主要成分多为石和砂，这与裸地的组成成分基本相同[38]，造成了不透水面和裸地的光谱特征十分相似[39-40](图1)，尽管有的指数在此方面作了改进尝试，PII在组合波段的基础上设置土壤线，NDISI在设计时引入了热红外波段，但是实验结果表明，现有指数未能在光谱特征上有效区分裸地和不透水面，进而影响不透水面指数在大区域的推广应用。此外，部分不透水面指数漏提了郊区的不透水面，这是由于郊区不透水面面积较小，容易与周边地类产生混合像元效应，不透水面光谱特征不显著。

为了进一步展示不透水面指数在不同区域的提取精度，研究选取了2个实验区的6个典型局部区域(表4)。表4中的6行分别表示从2个实验区内选取的6个局部区域(图4—5)，第2列为局部区域的Landsat8卫星遥感影像，第3列为所对应的不透水面真值数据，第4—11列分别为各个不透水面指数的最优提取结果(白色代表不透水面)。从表4中可以看出，区域1的不透水面分布密集，现有指数对该区域中心城区的不透水面提取效果较好，其中PII，BCI和RRI与不透水面真值数据吻合度最高，NDBI，UI，IBI和NDISI将部分水体误提取为不透水面，NDISI漏提了少量不透水面。区域1较好的提取结果符合大部分不透水面指数只考虑城市地区的设计思路。区域2内不透水面的面积较小，除NDISI外，大部分指数均能准确提取不透水面。在区域3，现有指数都将位于区域中心的裸地误提取为不透水面。区域4和区域5内分布着裸地，不透水面指数的表现均受到了裸地的干扰。而对于裸地大面积分布的区域6，现有指数均未有效增强不透水面信息，人工不透水面提取结果受到裸岩等自然不透水面的干扰，造成大面积裸地被误提取为人工不透水面，不透水面提取精度较低，这将影响不透水面指数在地表覆盖类型较为复杂的大区域中的应用。

表4 不透水面指数在局部区域的提取结果Tab.4 Example zoomed images showing the impervious surface extraction results

3 结论

本文选择具有不同地表覆盖特点的实验区，利用Landsat8卫星遥感影像，较为全面地测试了现有8种主要不透水面指数的提取精度，得到如下结论：

1)在不透水面分布密集、裸地面积小的区域，8种不透水面指数中PII指数的不透水面提取精度最高(约90%)，其次是RRI和BCI指数(约87%)，而NDISI，NDBI和IBI等指数的提取精度较低(约73%)。

2)在不透水面分布稀疏、裸地面积较大的区域，现有不透水面指数未能有效剔除裸地干扰，不透水面提取精度普遍较低(约71%)。

3)在裸地广泛分布的区域，不透水面指数提取精度偏低，不利于监测大区域的城市不透水面动态变化，反映出目前不透水面指数设计的局限性。

对于大区域的不透水面制图，指数法操作简单、不需要训练样本，但因其无法有效区分不透水面和裸地，限制了指数法在大区域不透水面制图中的应用。今后对不透水面指数的改进，可从指数设计原理出发，在利用地类光谱特征的基础上融入多源数据(例如夜间灯光数据)，以有效识别不透水面。