基于Sentinel-2A和Landsat8的城市不透水面的提取

赵怡，许剑辉，钟凯文，王云鹏，胡泓达，吴萍昊

(1.广东省科学院广州地理研究所/广东省遥感与地理信息应用重点实验室/广东省地理时空大数据工程实验室，广州 510070；2.南方海洋科学与工程广东省实验室(广州)，广州 511458；3.中国科学院广州地球化学研究所，广州 510640；4.中国科学院大学，北京 100049)

0 引言

城市不透水面是快速城市化的典型产物[1-3]。在城市土地利用中，植被和土壤等自然地表逐渐减少，取而代之的是不透水面。不透水面绝大部分是由人工材料制成，水体无法渗透，例如道路(水泥、沥青等材料)和建筑(水泥、岩石等材料)[4]。许多研究表明，不透水面是城市重要的土地覆盖类型，也是城市生态环境相关问题的重要参考指数，例如城市内涝和暴雨天气等[5]。因此，不透水面的遥感提取已经成为城市遥感的研究热点[6-13]。

在城市不透水面遥感提取中，线性光谱混合分解(linear spectral mixture analysis,LSMA)是常用的方法[14-15]，该方法基于亚像元尺度和植被-不透水面-土壤(vegetation-impervious-surface,V-I-S)模型[16]，假设遥感影像中混合像元的反射率是该像元中所有地物端元反射率的线性组合[17]。首先从混合像元内各类端元的纯净像元获取端元的光谱信息，然后通过光谱解混计算，可以得到混合像元中高反射率和低反射率不透水面的面积比例，得到不透水面的面积比例，即不透水面盖度。由于环境影响，以及在特定波段上不透水面的光谱特征与其他地物相似，导致不透水面的提取精度较低[18]。在LSMA解混过程中，端元选择是至关重要的一步[19]，原因在于LSMA是基于某一端元光谱特性，判断混像元内该类端元的面积比例，因此，所选取端元的光谱信息是进行解混混合像元的基础。

端元的纯净像元光谱信息来源有两个渠道：一是从标准地物光谱库中直接选取目标地物的光谱曲线，二是从实际目标影像中选择纯净像元，得到不同地物的光谱曲线。由于地物的光谱特征曲线会受到环境因素的影响，且不同地物之间也会产生相互作用，因此，标准地物光谱库中的光谱曲线并不能较好地表示目标地物，故而大部分研究选择第二种方式获取端元的光谱曲线[20]。传统的LSMA在端元中采用最小噪声分离变换(minimum noise fraction rotation，MNF)[21]将遥感影像的信息集中到前几个波段，达到去除噪声的目的；然后通过计算纯净像元指数(purity pixel index，PPI)[22]筛选出纯净度较高的像元；最后通过N维可视化选择纯净像元。此外，N-FINDAR[23]算法也常用于纯净端元的选择。

为了描述城市进程发展，长时间序列的遥感影像必不可少。在目前的遥感发展中，能满足长时间序列要求的影像大都是中等空间分辨率的。然而在中等空间分辨率的遥感影像(例如30 m空间分辨率的Landsat8 OLI影像)中通常存在大量的混合像元，纯净像元选择存在不确定性，其精度直接影响不透水面的提取精度。而在较高空间分辨率的遥感影像中(例如10 m空间分辨率Sentinel-2A MSI影像)，可以选择更多的纯净像元。在同一时期，目标地物在特定波段上的反射率是一定的，这就意味着，高空间分辨率中端元的光谱特性与中等空间空间分辨率中相似。为了获取更为可靠的纯净像元，本研究提出结合Sentinel-2A MSI影像改进从Landsat8 OLI 选择端元的光谱曲线。

此外，由于地物的光谱曲线受外界影响，且端元选择存在误差，因此，经过LSMA解混的不透水面盖度仍然存在错分现象，即在土壤分量中可能含有少量的不透水面，而在不透水面分量中也存在土壤或植被盖度。为了解决以上问题，本研究提出解混结果优化方案，即在像元尺度上，基于阈值分割法，利用归一化植被指数(normalized differential vegetation index，NDVI)和干旱裸土指数(dry bare-soil index，DBSI)识别不透水面、植被和土壤，从而分离出错分的像元，提高不透水面盖度提取精度。

1 研究区概况与数据源

1.1 研究区概况

研究区影像如图1所示。广州市位于E112°57′～114°30′，N22°26′～23°56′之间，是广东省的省会，属于国际中心城市，是国际商贸中心和综合交通枢纽，是粤港澳大湾区、泛珠江三角洲经济区的中心城市以及一带一路的枢纽城市。本研究选择广州市6个主城区，包括越秀区，海珠区，天河区，白云区，黄埔区，荔湾区(图1)。研究区域内土地利用类型丰富，包括大面积的植被区域(森林、农田)，裸土区域，不同密度的居民用地和商业用地。

图1 研究区Landsat8 B4(R)，B3(G)，B2(B)合成影像Fig.1 Study area image combined with Landsat8 B4(R)，B3(G)，B2(B)

1.2 数据源

Sentinel-2A 卫星发射于2015年6月23日，是“全球环境与安全监测”计划的第二颗卫星。该卫星的重访周期为十天，幅宽为290 km,携带的多光谱成像仪(multispectral instrument，MSI)有13个波段，从可见光和近红外到短波红外，具有不同的空间分辨率[24-25]。本研究选取2017年11月1日的Sentinel-2A MSI遥感影像的B2，B3，B4，B8，B11，B12作为研究目标，采用Sen2cor处理器进行大气校正。此外，通过哨兵系列卫星产品综合应用平台(Sentinel Application Platform，SNAP)，将B11，B12的遥感影像进行重采样，使其空间分辨率为10 m。

Landsat 8卫星发射于2013年2月11日，是美国陆地卫星计划的第8颗卫星。该卫星的重访周期为16 d，携带的陆地成像仪(operational land imager,OLI)有9个波段，其中包括7个常用波段(B1—B7，空间分辨率为30 m)，一个全色波段(B8，空间分辨率为15 m)和一个云检测波段(B9，空间分辨率为30 m)。本研究选取2017年10月23日的Landsat8 OLI遥感影像的B1—B7作为研究目标，通过辐射定标和大气校正，将影像的DN值转换为地表反射率。地表发射率产品是在美国地质调查局(United States Geological Survey，USGS)地球资源观测(Earth Resources Observation and Science，EROS)和科学中心科学处理机构(Center Science Processing Architecture，ESPA)(https://espa.cr.usgs.gov/)[26]的定制页面进行定制得到。

本研究采用的Sentinel-2A MSI 和Landsat8 OLI数据参数表1所示。

表1 Landsat8 OLI 和Sentinel-2A MSI 数据介绍Tab.1 Landsat8 OLI and Sentinel-2A MSI data

此外，采用2017年WorldView-2遥感影像(空间分辨率为1.8 m)，通过对样本区域内的不透水面进行矢量化，计算不透水面面积，作为精度验证数据。

2 研究方法

本研究所提出的方法分为3个步骤：首先，从Sentinel-2A和Landsat8 OLI 遥感影像中选取纯净像元，获取线性光谱解混所需的端元光谱曲线；然后利用传统线性光谱解混方法对两种数据进行光谱解混；最后，利用NDVI和DBSI两种光谱指数，对解混结果进行后处理，得到精度较高的城市不透水面比例。本文整体流程如图2所示。

图2 流程图Fig.2 Flow chart

2.1 端元选择优化方案

端元选择是LSMA过程中的重要步骤。理论上，作为光谱解混的基础，所选取的端元应是遥感影像上的纯净像元(只包括植被，裸土，高反射率不透水面与低反射率不透水面4类地物中的一类)，然后利用所选纯净像元的光谱信息，对混合像元进行光谱解混，计算得到混合像元内每一种地物的面积比例。因此，实验中所选端元的纯净度直接影响解混的精度。Sentinel-2A MSI影像的分辨率较高，像元大小为10 m×10 m，与Landsat8 OLI影像(像元尺寸为30 m×30 m)相比，其像元的纯净度较高，纯净像元的选取也比较容易。此外，在城市环境中，土壤分布比较破碎，很少有大面积的裸土呈现，因此，像元面积越小，裸土纯净像元就越多，也就是说在Sentinel-2A影像上裸土的纯净像元更多。

为了提高端元光谱的选去精度，本研究提出端元选取优化方案：结合Sentinel-2A MSI影像优化从Landsat8选取的端元光谱信息。具体过程如下：

1)采用传统方法对Landsat8 OLI影像进行端元选取，将所选纯净像元的光谱信息以文件的形式输出。

2)利用选择感兴趣区域(region of interest，ROI)的方式对Sentinel-2A影像进行端元选取，选择植被、裸土、高反射率不透水面与低反射率不透水面4类地物的纯净像元，并统计其光谱信息，生成端元光谱曲线，并以文件的形式输出。

3)利用从Sentinel-2A MSI 影像中选取的的纯净像元在每个波段(B2,B3,B4,B8,B11,B12)上的反射率替换Landsat8 OLI影像中选取的纯净像元在相应波段(B2—B7)的反射率，生成一个新的端元光谱信息文件，作为线性光谱混合分解的基础(从表1可以看出，Sentinel-2A MSI影像的6个波段与Landsat8 OLI的前6个波段的波谱范围比较接近，可以视为相同的波段)。

2.2 线性光谱混合分解(LSMA)

诸多研究表明，LSMA具有较强的理论基础。该方法假设遥感影像中混合像元在某一波段上的反射率是组成混合像元的各类端元在该波段上的反射率的线性组合，像元反射率越接近某类端元的反射率，就说明这种端元在该像元中的面积比例越大[27]。LSMA公式为[28]：

(1)

(2)

式中：i=1,2…,M,M为波段数；n为端元数目；Ri为混合像元在i波段的反射率；fk为端元k的面积比例；ERi为i波段的计算残差；RMS为计算残差，以此评价解混结果的优劣[29]。

由于混合像元中多类端元的光谱之间会发生相互作用，此外，环境等因素也会影响端元的光谱信息。因此在解混过程中，会存在不透水面被错误估计的问题，例如在不透水面覆盖度低于20%的区域，采用LSMA解混混合像元时，不透水面比例会被严重高估[3]。

由于不透水面的低反射率部分与水体部分反射率都比较低，十分接近，为了减少解混误差，在数据预处理阶段，利用改进的归一化水体指数(modified normalized difference water index,MNDWI)[30]将水体去除，计算公式为：

(3)

2.3 解混结果优化方案

在LSMA结果中，解混误差普遍存在，具体表现为：①在土壤盖度中含有不透水面，即不透水面被错分为土壤；②在不透水面盖度中(高反射率与低反射率盖度之和)含有植被和土壤，即少量植被和土壤被错分为不透水面。为了降低以上误差，提高城市不透水面的提取精度，本研究结合像元尺度，利用DBSI和NDVI，通过二者的阈值，提取解混过程中被错分的像元，将其归还到相应的盖度中去。其中，DBSI和NDVI的计算公式分别为：

(4)

(5)

式中：BNIR为影像的近红外波段(Landsat8 B5)；BR为影像的红色波段(Landsat8 B4)。

优化过程分为3步，具体如下：

1)从土壤盖度分量中提取错分的不透水面盖度。利用DBSI阈值筛选土壤盖度中的不透水面像元，借助DBSI影像的直方图，通过实验确定0.1作为该步DBSI的阈值，即在土壤盖度中，DBSI值小于0.1的像元被认为是不透水面，将这部分像元的解混值认为是不透水面盖度，对应图2中的不透水面盖度(1)；

2)将LSMA解混结果中的高反照率和低反照率不透水面盖度相加，得到LSMA解混遥感影像得到的不透水面盖度，然后将上一操作中的不透水面盖度(1)与之相加，得到不透水面盖度(2)；

3)去除上一操作中的不透水面盖度(2)中含有的少量植被和土壤盖度，得到最终的不透水面盖度。这部分误差来源于LSMA本身的解混误差。理论上，NVDI值越接近1，说明目标像元越接近植被，负值代表水体，0代表土壤或者沙漠化土地[32]。为了保证识别的像元为纯净的植被和土壤，借助DBSI和NDVI影像的直方图，确定DBSI和NDVI阈值分别为0.2和0.3，即DBSI<0.2且NDVI>0.3的像元被认为是植被，DBSI>0.2且NDVI<0.3的像元被认为是土壤。

在以往研究中，DBSI>0.26的像元被认为是裸土，NDVI值接近于1被认为是植被。在步骤1)中，从土壤盖度中提取出不透水面，为了尽量保留LSMA结果的土壤盖度中土壤的真实面积比例，DBSI值应稍小于理论值，因此DBSI值选择为0.1。而在步骤3)中，从不透水面中去除植被和土壤，为了尽量保留真实不透水面的面积比例，DBSI和NDVI阈值应稍大于理论值，因此NDVI阈值被确定为0.3；然而在实际实验过程中，经过反复试验得知，DBSI阈值为0.2能较好地区分不透水面与其他地物。

2.4 精度验证

为了验证本文所提出的方法对提取城市不透水面的有效性，随机选取了171个样本区域，样本大小为480 m×480 m，并且采用WorldView-2 影像对样本区域进行矢量化，得到不透水面的样本真实面积比例。通过系统误差(system error，SE)、平均绝对误差(mean absolute error，MAE)、均方根误差(root mean square error，RMSE)和决定系数R2来衡量实验结果的精度。前3个指标的计算公式分别为：

(6)

(7)

(8)

“喝西凤、吃泡馍、吼秦腔”是被概括成了“关中人的形象”，当然，这也是游人必不可少的原生态体验。性烈醇劲的西凤酒，劲道味重大碗泡馍，街市上都随处可见。烈酒好面，在北方并不罕见，唯有这秦腔却是西安人和关中汉所独有特征。“八百里秦川尘土飞扬，三千万楞娃齐吼秦腔。”听秦腔、吼秦腔至今仍是西安人日常生活中的一大乐事，和兵马俑一样，秦腔也可被视为陕西的象征。

3 结果与分析

3.1 端元选取

图3展示了不同方法选择的端元光谱曲线。在端元选取过程中，由于影像自身的原因，以及环境因素，采用传统端元选择方法选取的土壤和高反射率不透水面这两类端元容易混淆，且精度不高，如图3(a)所示。结合Sentinel-2A影像选择较纯净的像元，优化传统方法选择的端元光谱信息，如图3(b)所示，降低了各类地物端元光谱的混淆问题，且与实际地物光谱信息更接近，精度更高。

(a)传统端元选择方法 (b)端元选取优化方案

3.2 结果比较分析

通过采集到的端元光谱曲线，利用LSMA对Landsat8 OLI遥感影像进行像元解混，得到如图4所示的不透水面盖度分布图。通过对比实验分析，结果显示在城市密度较高的区域，图4(b)中不透水面盖度比图4(a)明显增加。图4(c)图采用了解混结果优化方案，结果显示在森林等透水区域，不透水面盖度大幅降低，接近于0，与实际不透水面分布十分接近。

(a)传统LSMA (b)结合端元优化方案的LSMA (c)结合端元优化方案和解混结果优化方案的LSMA

通过精度评价，如图5和表2所示，结果表明，通过端元选择优化方案，将高空间分辨率较高的Sentinel-2A MSI影像上选取的纯净像元的光谱信息代替相应波段上从Landsat8 OLI影像中选取的端元的光谱信息，经过LSMA解混后，不透水面盖度的精度得到提高，从线性拟合结果来看，图5(b)中的观测值(利用不同LSMA提取的不透水面盖度)与真值比图5(a)更接近，拟合线周围的样本点更为密集。如表2，结合优化方案的LSMA所提取的不透水面盖度的R2比传统LSMA高16%。

(a)传统LSMA (b)结合端元优化方案的LSMA (c)结合端元优化方案和解混 结果优化方案的LSMA

表2 不同方法提取的城市不透水面盖度精度分析Tab.2 accuracy assessment of urban IS fraction extracted by different methods

综上所述，本研究所提出的端元优化方案可以使LSMA更加准确地识别不透水面，降低了不透水面盖度被低估的问题(在不透水面覆盖度较高的区域)，提高城市化程度较高的区域中不透水面盖度的提取精度。

但是图5(a)和(b)同时存在一个问题，即在不透水面实际盖度较低的区域(横轴接近于0的位置)，观测值明显比真值大。也就是说，在不透水面覆盖度较低的区域，传统LSMA和结合端元光谱优化方案的LSMA会高估不透水面盖度。图5(c)显示了结合端元优化方案和解混结果优化方案的LSMA提取的水面盖度与真值的线性拟合结果，该图表明，在不透水面盖度较低的区域，观测值明显降低。如表2所示，结合端元优化方案和解混结果优化方案的LSMA所提取的不透水面盖度的R2比传统LSMA高20%，比仅结合端元优化方案的LSMA高3%，且SE为-0.066，说明整体系统误差较小，该方法所得到的不透水面盖度与实际情况接近。以上结果说明通过解混结果优化方案，大大改善了不透水面被高估的问题(在不透水面盖度较低的区域)。

图6展示了通过不同方式获取的不透水面盖度细节图，通过对比图6(b)和(c)，图6(f)和(g)来看，在城市化较快的区域，结合端元优化方案可以提高不透水面盖度，从而提高提取精度。从图6(d)和(h)来看，结合解混结果优化方案，可以明显去除森林和土壤等透水区域内不透水面盖度，同时保持原有的不透水面盖度，有效解决不透水面被高估的问题。

(a)Sentinel-2AMSI影像 (b)传统LSMA提取的不透水面盖度 (c)结合端元选择优化方案的LSMA提取的不透水面盖度 (d)结合端元选择优化方案和解混结果优化方案的LSMA提取的不透水面盖度

4 结论

1)本研究针对Landsat8 OLI遥感影像的端元选择问题提出了端元选择优化方案；针对LSMA后处理问题，结合解混结果优化方案，以此提高不透水面的提取精度。

2)与Landsat系列影像相比，Sentinel-2A MSI影像的空间分辨率较高，可以快速选择相对比较纯净的像元作为LSMA的端元，然后将选取像元的光谱信息代替对应在波段上从Landsat8 OLI影像上选择的端元的光谱信息，优化了光谱解混所必需的端元光谱信息，从而达到了提高LSMA对混合像元解混精度的目的。

3)对于LSMA后的各地物盖度分量，结合解混结果优化方案，利用DBSI和NDVI两个指数作为识别植被和土壤的重要参数，可以比较准确地从土壤盖度分量中提取被错分的不透水面盖度分量，且剔除了不透水面盖度分量中被错分的植被和土壤盖度分量，整体提高了不透水面盖度的提取精度。

以上两种优化方案解决了传统LSMA在不透水面盖度较高的区域(例如繁华的商业区，聚集的居民地等区域)被低估的问题，以及在不透水面盖度较低的区域(例如大面积的森林区域等)被高估的问题。然而，本研究在DBSI和NDVI阈值选择过程中，大部分操作依赖统计实验，未来可以从两个指数的阈值入手，为不透水面盖度提取提供可靠的理论方法，从而达到提高不透水面盖度提取精度的目的。