地表不透水面比例遥感反演研究综述

左家旗， 王泽根， 边金虎， 李爱农， 雷光斌， 张正健

(1.西南石油大学土木工程与建筑学院，成都 610500； 2.中国科学院水利部成都山地灾害与环境研究所，成都 610041)

0 引言

城市化是21世纪的一个显著特征，全球已有超过一半的人口居住于城市，预计到2050年将有2/3的人口居住于城市[1-2]。城市化是全球一体化和社会现代化的必由之路，在带来经济增长和生活便利的同时，也将带来住房紧张、交通拥挤和生态环境污染等诸多问题。监测城市化进程，是实现可持续发展和协调人类社会与自然环境之间关系的关键。城市化的一个重要特征就是大量的不透水面取代了耕地、植被和水体等自然景观。不透水面通常指能够隔离地表水渗透至土壤的表面，包括以裸岩为主的天然不透水面和人工不透水面[3]。城市化研究中的不透水面均指人工不透水面，包括道路、屋顶和停车场等。不透水面比例是指单位地表面积内不透水面所占的比例。不透水面的增加不仅能反映城市的发展与扩张，同时将引起城市表层结构变化，进而改变区域气候并影响城市生态系统服务功能。因此，准确估算不透水面比例对城市发展规划、生态环境评价及保护和资源管理等具有重要意义。

早期不透水面研究采用的数据主要来源于平面地图，Pang等[4]使用平面地图数据对香港2006年不透水面积进行估算； Wu等[5]基于时间序列平面地图数据监测美国爱荷华州中部地区1940—2011年间的不透水面变化。此类方法不适合地图数据缺失或不完整的地区，且制作地图数据的时间和人力成本高，不利于开展大范围的城市扩张监测和变化分析。遥感技术的发展以及卫星数据的开放获取，使利用遥感技术反演不透水面成为可能。遥感已成为当前获取大区域、长时期不透水面数据的重要技术手段。早期的不透水面遥感反演研究主要以不透水面制图为主[6-7]。Xian等[8]绘制了美国佛罗里达州坦帕湾1991—2002年间的不透水面空间分布图，并利用元胞自动机(cellular automata，CA)模型预测未来城市的范围； Patel等[9]使用多个光谱指数组合作为输入变量，构建了多种人工神经网络模型提取印度兰契地区的不透水面。此类不透水面制图将遥感像元“硬”分类为透水面和不透水面，忽略了遥感影像中普遍存在的混合像元问题，分类结果存在一定不确定性。为更准确地反映不透水面的实际状况，不透水面的研究逐渐由不透水面制图向像元内不透水面比例反演转变，从亚像元角度精细刻画城市不透水面的空间格局。程熙等[10]利用线性混合分解的方法将香港的土地覆盖类型分为高、低反射率不透水面，植被和裸土4类，将得到的高、低反射率不透水面的比例相加获得总不透水面比例； 高志宏等[11]用从QuickBird高空间分辨率影像获取的不透水面比例数据和Landsat TM影像的光谱数据构成训练样本集，基于回归树模型(classification and regression trees，CART)反演得到了2002—2006年间山东省泰安市30 m空间分辨率的不透水面比例。

线性混合分解[10]和机器学习[12-14]等不透水面反演方法的提出与应用，使遥感在不透水面反演精度和效率等方面都有了较大提高。不透水面比例的遥感反演经历了单时相到时间序列的发展，同时各种不透水面比例时间序列分析方法得以发展。本文将以此为脉络，总结单时相不透水面比例反演方法、时间序列不透水面比例反演方法和精度验证方法，并对各种方法的优缺点进行剖析； 最后对现阶段不透水面比例反演中存在的问题，提出针对性的建议，并展望未来的发展方向。

1 不透水面比例反演方法

现阶段不透水面比例的反演方法主要包括2类： 一类是以Ridd[15]提出的城市生态组分模型(vegetation-impervious-soil, V-I-S)为理论基础，将城市区域地表覆被类型表达为植被、不透水面与裸土等生态要素按不同比例的组合，采用光谱分解的方法反演不透水面比例； 第二类则是基于机器学习算法，构建不透水面比例与遥感影像各波段的光谱反射率或各类衍生指数的线性或非线性模型，从而反演不透水面比例。

1.1 光谱分解法反演不透水面比例

光谱混合分解模型中，将像元在某一波段的光谱反射率表示为占一定比例的各个基本组分(端元)反射率的组合。基于以下假设： 在瞬时视场下，各组分光谱按一定数学规则混合，其比例由相关组分光谱的丰度决定，即遥感器瞬时视场内的地面目标辐射亮度仅与各组分所占比例有关，而与其形状及分布等无关[16]。光谱分解的过程包括端元提取和定量化端元在影像像元中的百分比[17-20]2部分，其关键在于计算混合像元内各纯净端元的组成比例，并把混合像元的光谱分解成多种纯净端元的组合。由于不透水面类型多样，光谱分解法反演不透水面比例的研究重点在于如何提取各种不透水面类型的端元光谱值以及如何有效区分不透水面与裸土、裸岩等亮地表。

1.1.1 端元确定

根据端元的选择，可以将光谱分解模型分为固定端元分解[21-22]和可变端元分解[18，20]2类。固定端元分解是以V-I-S模型为理论依据，将城市地表视为植被、不透水面和裸土三者的组合。Wu等[23]基于V-I-S模型进一步将不透水面分为高反射率地物和低反射率地物，对Landsat ETM+影像进行光谱混合分解，最后将得到的高、低反射率地物分量影像相加，得到总不透水面比例。可变端元分解方法的提出是为了解决复杂城市景观中的同物异谱、同谱异物现象。其基本思想是建立具有一定数量的端元光谱库，对每一个像元运行所有端元组合，选择最能代表该像元的最优端元组合。Bian等[24]提出一种自适应端元选择线性光谱混合模型(adaptive endmember selection linear spectral mixture model，ASLSMM)，通过考虑每个像元的端元动态组合来提高估算精度，其估算结果的均方根误差(root-mean-squared error，RMSE)相比于固定端元分解降低了5.7%。

端元获取是混合分解的关键，端元的提取精度直接影响地物拟合优度和混合分解的精度。端元的获取途径有2类： 一是直接从遥感影像中提取，此方法优势在于获取相对容易，端元光谱和影像具有相同尺度，端元噪声可与影像噪声相互抵消，但此方法对空间分辨率要求较高，影像空间分辨率越低，提取效果越差； 二是从已知端元库中获取，其优势是可以获得较纯净的端元，但端元库中端元是在理想状态下得到的，与自然环境下的端元存在差距，且端元光谱缺少与实际影像的校正也会产生误差。常用的获取方法是通过最小噪声分离(minimum noise fraction，MNF)和纯净像元指数(pure pixel index，PPI)等辅助手段直接从影像中提取，程熙等[10]对影像进行MNF变换，基于变换后的前3个波段MNF1(0.45～0.52 μm)，MNF2(0.52～0.60 μm)和MNF3(0.63～0.69 μm)的特征分布(图1)，提取三角形散点图，三角形的顶点即为端元，运用线性混合分解方法得到天津市的不透水面比例数据。

(a) MNF1和MNF2反射率散点图 (b) MNF1和MNF3反射率散点图 (c) MNF2和MNF3反射率散点图

图1 MNF变换前3个波段的特征分布(改编自程熙等[10])

Fig.1FeaturedistributionofthefirstthreeMNFcomponents(adaptedfromChengetal.[10])

1.1.2 不透水面与亮地表的区分

不透水面端元和裸土、裸岩等亮地表端元由于反射率值比较接近而难以区分。根据光谱不同，易混淆的裸土、裸岩等亮地表端元主要分为反射率相对较低的裸土、建筑工地等和反射率相对较高的沙地、裸露基岩等； 不透水面端元也分为低反射率和高反射率2个亚类，前者以沥青材料为主，后者以混凝土材料为主。仅利用PPI方法不能很好地区分这2类端元，且获得的纯像元数目较少。目前这2类地物区分的主要方法是由高空间分辨率影像进行目视判读。邓蕾等[25]利用PPI和手动选取相结合的方法提取纯像元，探讨了选取这2类端元时内部高、低反射率像元个数比例对分解结果的影响，结果表明： 裸土、裸岩等亮地表端元选取时兼顾低反射率的裸土和高反射率的裸露基岩，可以整体提高模型的分解适宜度和分解精度； 不透水面端元中高反射率像元比例提高有助于改善整体分解效果，却又不能全部提取高反射率纯像元。而且，这种手动选取纯像元的方法主观性较强，所选端元中亚类间的比例不能准确量化，而端元中亚类的比例又直接影响反演精度，从而使反演结果存在一定的不确定性。此外，在进行大区域如国家、省区尺度的不透水面反演时，手动选取纯像元的方法仍然存在很大的不确定性。

1.1.3 线性和非线性分解模型

线性光谱分解是最常用的一种光谱分解方法。该方法认为每一波段中单一像元的反射率，为各端元组分反射率与其各自比例的线性组合。在不透水面比例反演中多将水体进行掩模，用高、低反射率地物，植被和裸土4种光谱端元的线性组合来模拟影像的波谱特征，以高、低反射率地物比例之和作为总不透水面比例[23]。但由于大气多次散射和仪器校准等因素影响，光谱混合是高度非线性的[26]，因此利用非线性模型进行分解具有较高精度。如多项式非线性混合模型[17]、基于核函数的分解方法[27]，可涵盖广泛的非线性混合类型，且实现较为灵活。李慧等[28]提出基于支持向量机(support vector machine，SVM)策略的端元可变非线性混合分解方法，为每个像元构造优化端元集，提高分解精度，实验表明该方法分解结果与线性分解相比，平均绝对误差(mean absolute error，MAE)降低了约20%。虽然非线性分解模型因其广泛性和灵活性得到一定的应用，但与线性混合分解模型相比还很不成熟，仍需要进一步深化研究和应用。另外，线性与非线性分解模型在函数拟合的过程中，将同时计算拟合残差和RMSE，可根据误差值的大小及分布来评估模型的反演效果。利用高空间分辨率影像和光谱分解的不透水面比例结果进行对比分析，可对模型分解效果的真实性作进一步的检验。

1.2 机器学习模型反演不透水面比例

近年来，机器学习算法在不透水面比例反演研究中逐渐得到广泛应用。机器学习算法可以避免预先进行端元的提取和优选，相比于光谱混合分解来说，能够最小化人工时间成本、在一定程度上实现自动化，适用于大区域不透水面反演。机器学习算法在处理数据噪声、复杂数据源以及有限训练样本方面有较好的表现[29]。主要包括： 决策树、SVM、多层感知器(multi-layer perceptron，MLP)、随机森林(random forest，RF)等。机器学习算法反演不透水面比例精度很大程度上决定于训练样本的质量，训练样本主要包括由波段反射率和光谱指数构成的输入变量、由高空间分辨率影像提取得到的不透水面“真值”数据或现有土地覆盖数据构成的参考数据集。下面将从输入变量、参考数据和机器学习模型3个方面进行总结。

1.2.1 输入变量

波段反射率是反映不透水面在不同光谱波段吸收或反射能力的最基本输入变量。基于单季相或多季相的波段反射率特征建立机器学习模型是不透水面比例反演的常用方法。现有研究表明，选用生长季和非生长季的波段，辅以波段变换后的特征构建训练样本集，将有助于地物的区分识别，提高模型反演精度[11,13]。利用各类地物在各波段上的反射率差异，有针对性地构建光谱指数可以使不透水面比例反演精度进一步提高。各学者根据不透水面在热红外波段辐射率高、在近红外波段反射率低这一特性，提出了城市指数(urban index，UI)、归一化差值建筑指数(normalized difference built-up index，NDBI)、归一化差值不透水面指数(normalized difference impervious surface index，NDISI)等[30-32]。随着夜间灯光数据的应用，结合多光谱数据和夜间灯光数据的一系列指数被提出，如人类居住指数(human settlement index，HSI)、植被调节型夜间灯光城市指数(vegetation adjusted nighttime light urban index，VANUI)、归一化不透水面指数(normalized impervious surface index，NISI)等[33-35]。

输入变量的选择和组合是决定模型反演精度和计算效率的关键。选择变量时既需要保证所选波段和光谱指数使不透水面与其他地物的区分度最大，也需兼顾考虑模型计算效率，进行必要数据筛选避免冗余。对输入变量进行组合建立多个机器学习模型，再对比各模型精度和计算时间，选取最优模型是当前的一种普遍方式[12,14]。Hu等[36]利用Landsat TM数据、国防气象卫星线性扫描业务系统(the defence meteorological satellite program/operational linescan system，DMSP/OLS)夜间灯光数据、数字高程模型(digital elevation model，DEM)及其相应的坡度数据，以北京市为研究区，对输入变量进行组合建立了多个CART模型，并进行了模型精度对比分析。各组合模型及其精度分别如表1所示，其结果表明并非输入变量越多反演精度越高，选择对反演贡献较大的变量进行建模可以在保证精度的同时节约计算资源。同样，在类似研究中有学者对不同输入变量和不同机器学习模型开展地物识别的效果进行了对比分析，结果表明在保证输入变量相同的情况下，不同机器学习模型的选择对地物提取精度的差异较小，而合理的选取输入变量则可以使提取精度显著提高[37]。因此，在不透水面比例反演过程中需要结合先验知识和地理特情，选取对地物识别贡献较大的变量建立定量反演模型。

表1 各CART模型和相应输入变量[36]Tab.1 CART models and corresponding input variables

注： On b1-7和Off b1-7分别包括生长季和落叶季影像的第1—7波段； On TC和Off TC分别包括生长季和落叶季影像经过缨帽变换后的亮度、绿度、湿度； Light为夜间灯光数据； Slope为坡度；AE为绝对误差；RE为相对误差；R为相关系数。

1.2.2 参考数据

参考数据即作为不透水面比例反演模型训练和反演精度验证的数据。获取参考数据主要有2条途径： 一是利用现有土地覆盖数据集或地方政府测绘的大比例尺数字化城市地图，如Sexton等[38]将马里兰州数字化的道路、建筑物、停车场和人行道等地图数据作为不透水面数据； 二是以高空间分辨率影像为数据源，对城市地物进行分类，绘制不透水面空间分布图[39-40]，如Yang等[13]将IKONOS高空间分辨率影像进行监督分类获得空间分辨率为1 m的不透水面分类结果图。将以上2种途径获得的不透水面分布图栅格化，然后进行聚合使之与反演所用影像的空间分辨率一致，得到不透水面比例参考数据集。

1.2.3 机器学习模型

机器学习模型能够较好地拟合线性或非线性函数，对于输入变量的变化具有较强的鲁棒性，不透水面比例的反演精度也较好。但现有研究多基于某一种机器学习模型，选择不同的模型参数或输入变量组合进行建模，分析不同条件下的模型表现，选择最优的模型进行反演[36，41]。对同一研究区基于多种机器学习模型分别进行不透水面比例反演的研究较为缺乏，难以定论何种模型最优，仅能从各研究反演结果大致分析，几种模型反演精度基本处在相近的范围[14,42-43]。表2总结了各模型的机理、必要参数及当前反演不透水面比例的精度水平。机器学习模型也存在共同的不足： ①训练样本集的大小、所选样本是否具有代表性对输出结果的精度有很大影响[29]； ②各模型必要参数的设定需经过多次尝试、对比才能选择出最优模型； ③由于年间光谱差异较大，仅利用单年训练数据得到的模型难以直接应用于其他年份的不透水面比例反演。

表2 不透水面机器学习反演模型对比Tab.2 Comparison of machine learning retrieval models of impervious surface

2 时间序列不透水面比例反演与分析

2.1 基于单时相影像逐年反演不透水面比例

为更好地监测城市化进程及其对城市环境的影响，学者在探索单一时期不透水面比例反演方法的同时，也关注了不透水面比例的动态变化趋势。常用的研究方法是基于单时相影像的反演方法，即对每一年份的不透水面比例单独进行反演，构建时间序列不透水面比例数据集，再结合人口、经济和政府规划等数据进行不透水面变化分析[44]。此类研究虽然使用了多个时间或者时间序列数据，但均是基于单一时相的影像来反演不透水面比例，对年际物候变化产生的影响考虑不足，没有充分利用时间特征信息，难以捕捉不透水面在连续时间上的细微变化。

2.2 基于多时相影像逐年反演不透水面比例

受气候条件和降雨空间分布的影响，不同季节、年份的遥感影像在光谱特征上存在差异，最显著的表现就是植被的物候变化。比如在旱季，休耕地等高反射率地物常与不透水面混淆，往往难以通过单一时相的数据区分。利用多季相、连续年份的影像数据开展不透水面比例反演是解决这一问题的有效途径。如叠加生长季和非生长季影像，基于叠加影像构建训练样本集，利用植被季节变化特征加强机器学习模型对季节性休耕农田和不透水面的识别[41]； 再如使用时间窗口计算连续年份的反射率中值和NDVI最大值，利用其带有的生长季信息加强地物识别，同时起到了滤波的效果，减少残余环境噪声和大气物候的影响[38]。

总的来说，利用多时相影像进行不透水面比例反演具有以下优势： ①长时间序列的观测能够捕捉到微弱的不透水面变化，更精细地表征不透水面变化趋势； ②连续年份或同年多季相影像提供的植被信息能够加强不透水面和裸地的区分度； ③利用时间滤波窗口进行滤波有助于消除部分残余噪声。然而，基于多时相影像的反演也存在明显的缺陷： 云层覆盖、传感器老化或损坏等原因导致影像不可用，难以获取日期接近、光谱一致性较好的密集时间序列数据集。

2.3 不透水面比例变化趋势分析方法

2.3.1 差值函数模型

差值函数模型即利用任意2 a的不透水面比例进行求差，是获得时间序列不透水面变化趋势最直接、简单的方法。Zhang等[45]对舟山岛1986—2011年间的不透水面变化进行监测，利用1986年和2011年这2期的不透水面比例求差，最终得到25 a间不透水面平均每年增长2.70 km2； Li等[46]利用线性光谱分解得到哈尔滨市1984—2010年间的不透水面比例数据，将每年的不透水面比例分为0～1的10个区间，同时对各区间赋予不同的权重，以求得当年哈尔滨市的城市重心，用重心移动的方向和距离来表现不透水面的空间变化，用年际差值作为衡量不透水面变化速率的指标。差值函数模型能直观地展现不透水面的总体增长趋势，但难以反映具体的变化过程。同时，求差可能由于某年的数据存在误差，而导致所得到的不透水面变化数据不够准确。

2.3.2 拟合函数模型

由于差值分析方法存在着一些问题： ①特定年份的误差将传播到任何使用该数据产生的不透水面变化数据中[47]； ②误差传播将导致变化结果的内在不一致； ③不透水面变化完成需要一段较长时间(如1～3 a)。因此，使用符合不透水面变化特征的函数模型能更细致、准确地描述不透水面在连续时间和空间上的动态变化。Song等[48]基于不透水面变化是不可逆的这一假设，将不透水面的变化分为3个阶段： 变化前稳定阶段、急剧增加或逐渐增加的变化阶段、变化后稳定阶段。采用Logistic函数来模拟这3个阶段，不透水面比例的公式为

(1)

式中：ISPt为t时间的不透水面比例，是随时间t变化的函数；a为不透水面变化量，b为变化的速率；c为变化的时间；d为变化前或后的ISP值。

基于上述模型得到4个参数后，分别用变化量、变化中间年份和变化持续时长来描述不透水面的变化(图2)。其中变化量用a表示； 变化中间年份用c/b来表示； 计算变化持续时长需确定变化开始和结束的年份，可从变化中间年份开始向前和向后搜索，以年变化率小于5%时对应的年份为变化开始和结束年。拟合函数的方法可以完整地描绘不透水面变化过程，但使用何种函数需要根据研究区域不透水面的变化特征来具体确定。

图2 Logistic函数模拟不透水面变化(改编自Song等[48])Fig.2 Simulating the ISP by using Logisticfunction(adapted from Song et al.[48])

3 精度验证方法

常用的不透水面比例反演精度验证方法包括： ①以训练数据中的不透水面比例参考数据构成验证集，用n折交叉验证的方法求取平均偏差误差(mean bias error，MBE)，MAE和RMSE[38]； ②在所得不透水面比例影像上随机选取验证点，以高空间分辨率遥感影像为数据源对验证点处的不透水面比例进行解译，将解译得到的不透水面数据和预测所得数据求差得到误差值[41,49]。这2种方法都存在各自的优势和不足，交叉验证利用已建立的数据集进行验证，较为方便，但验证数据和模型训练数据属同源数据，会影响验证结果的客观性； 采用其他来源的高空间分辨率影像进行验证，避免了使用同源数据，结果更加客观，但流程较为复杂且某些年份的高空间分辨率影像不易获取。

4 结论与展望

近20 a来，不透水面研究进展迅速，从分类制图到定量反演，从获取单一时期数据到时间序列变化监测，从简单的城市扩张定量化到人与环境的综合评价分析，在各方面都取得了重大突破。基于上述相关研究进展，综合国内外不透水面比例反演方法及发展趋势，笔者对今后遥感反演不透水面比例的方法及其应用有如下思考和认识。

1)端元获取是光谱混合分解的关键，是决定反演精度的主要因素。探究影像变换的新方法从中提取端元或寻求端元库与实际影像间的校正方法以保证纯净光谱和实际影像契合，是提高混合像元分解精度的2条途径。

2)建立充足且具有代表性的训练样本集是机器学习模型反演不透水面比例的必要条件。当前的研究大多以大城市或特大城市作为研究区，未将模型推广至具有不同环境条件的更大尺度区域。针对不同气候、地形条件下的区域进行多输入变量组合建模，通过模型对比确立各类型区域的输入变量选择方案，是实现大区域不透水面比例反演的有效方法。另外，从高空间分辨率影像提取不透水面并聚合获得不透水面比例是获取参考数据的重要手段，其提取精度决定了训练样本的质量。因此，解决高空间分辨率影像中阴影造成的高估现象和植被遮盖不透水面造成的低估问题，是进一步提高机器学习模型反演精度的突破方向。

3)机器学习模型对训练样本的选择要求较高，需要一定的人工参与，未来需加强样本自动选取方法的研究，实现高效、自动的不透水面反演。深度学习算法的发展降低了训练样本选取难度，算法对原始输入数据进行变换产生新的特征，并自动寻找特征间的关系建立函数模型，未来可针对新算法进行试验并对各方法对比分析。此外，发展不透水面等地物识别模型应该结合地物识别的机理，围绕“在识别中定量，在定量中识别”的基本思路。

4)裸土和不透水面的混淆仍是不透水面比例反演的一大难题。可基于多季相影像进行反演，如综合植被生长旺盛的夏季影像和冬季影像所获得的反演结果，将有利于裸土和不透水面的区分。深入挖掘多源辅助数据的潜力，如利用激光雷达(light detection and ranging，LiDAR)的高程信息区分建筑物和裸土； 利用Google Map抓取的路网信息提高道路提取精度； 利用人类感兴趣点(point of interest，POI)分布图获取城市范围剔除部分裸土等。

5)全球城市化发展对城市定量遥感提出了新的挑战，基于遥感的城市生态环境参数定量描述成为研究热点。结合国家发展规划和区域合作战略构想，不透水面的定量分析需要从城市扩张与土地覆被变化，拓展到生态、人居和热环境等多个主题，为城市环境健康、发展规划和资源管理等提供支持。在研究的空间尺度上，形成城市、城市群、国家、区域经济共同体以及全球等多尺度研究格局，分析不同尺度上城市化的影响； 在时间尺度上，充分利用长时间序列遥感影像的优势，实现对城市化定性、定时的综合研究。