基于空间光谱信息协同的城市不透水层提取方法比较研究

范佳辉，张亚丽，李明诗,2*

(1. 南京林业大学林学院，江苏 南京 210037；2.南京林业大学，南方现代林业协同创新中心，江苏 南京 210037)

城市不透水层包括车道、人行道、停车场和屋顶等[1]，是城市化程度的指标之一[2]，与城市土地覆盖分类、城市利用规划、城市生态系统等有着密切的关系[3-4]；因此，对城市不透水层进行准确制图，有利于促进城市的可持续发展和管理。有学者使用Landsat、SPOT及ASTER等中空间分辨率影像提取城市不透水层信息[5]，但由于城市景观存在高度异质性和破碎化特征，利用这些中空间分辨率遥感数据进行城市不透水层制图时存在着混合像元问题。高空间分辨率遥感数据如Quickbird、IKONOS、Worldview等能在一定程度上较好地解决混合像元问题，但是在城市环境中应用这些数据时，城市高层建筑在图像上形成的本影和投落阴影等问题[6]影响图像解译的准确性。为此，基于中空间分辨率图像的多种亚像元分解算法被广泛使用。

光谱混合模型是进行混合像元分解的前提。按照光谱混合因素可以将光谱混合模型分为线性混合模型(linear mixing model，LMM)和非线性混合模型(nonlinear mixing model，NLMM)[7]。线性混合是指大尺度的光谱混合且入射光仅与一种地物发生作用[8]，Ridd[9]提出了植被-不透水层-土壤模型(vegetation-impervious surface-soil，V-I-S)，Wu等[10]提出使用光谱混合分析方法(spectral mixing analysis，SMA)来提取城市不透水层。目前的非线性光谱解混算法分为两类: 基于模型的方法和不基于模型的数据驱动的方法[11]。混合调制匹配滤波(mixture tuned matched filtering，MTMF)[12]是将LMM与匹配滤波(matched filtering，MF)相结合发展而来的。Zhuo等[13]利用两种遥感数据，采用MTMF算法提取得到长三角、珠三角的城市不透水层覆盖度；Altmann等[14]利用后验非线性模型(post-nonlinear model，PPNM)对内华达州的矿区进行提取; Fan等[15]比较了用于森林高光谱模拟数据的线性与非线性模型，并通过实验得到FM模型(Fan model，FM)；夏俊士等[16]利用人工神经网络对徐州的城市不透水层进行了提取。

进行光谱混合模型求解的前提是存在可靠的地物类型终端单元(endmember)。现有的端元提取方法主要分为两种：其一，假定影像内存在纯净像元，可直接从影像本身识别并提取端元的端元识别算法，包括像元纯度指数(pixel purity index，PPI)[17]、单形体增长算法(simple growing algorithm，SGA)[18]等；其二，假定影像内不存在纯净像元，需要根据影像生成虚拟端元的端元生成算法，包括最小体积变换(minimum volume transform，MVT)[19]、外包单行体收缩(simplex shrink-wrap algorithm，SSWA)[20]等。这些端元提取方法主要借助于高光谱图像精细的光谱信息来实现，但对于提取多光谱影像端元不容易达到预期效果，而中/高空间分辨率多光谱影像通常含有全色波段，具有比多光谱波段更高的空间细节，以往在识别终端单元时容易忽略它。为此，本研究利用2018年6月6日获取的覆盖南京市的Landsat 8 OLI多光谱影像数据，构建空间纯像元指数SPPI[21]和像元纯度指数PPI来提取端元。随后再采用线性混合光谱模型、混合调制匹配滤波、双线性混合光谱模型以及BP神经网络来提取南京市不透水层，并利用高空间分辨率数据的目视解译结果进行不透水层信息提取的精度验证，以期为南京及类似中纬度地区特大城市的不透水层信息提取找寻可靠的方法。

1 材料与方法

1.1 研究区概况

研究区为南京市(118°22′～119°14′E，31°14′～32°37′N)，南北长约150 km，东西宽约70 km，总面积6 587 km2，建成区面积1 916.81 km2[22]。南京属于亚热带季风气候，年均气温14.6～16.4 ℃，年均降雨量800～1 000 mm。

1.2 数据及预处理

使用2018年6月6日获取的Landsat 8 OLI多光谱影像数据，来自美国地质调查局(USGS)地球资源观测和科学中心(EROS)(https://glovis.usgs.gov/)，轨道号p120r038。此外，还从Google Earth上获取了2018年的米级高空间分辨率影像数据，用来进行目视解译并验证各类提取方法的准确性。对于Landsat 8 OLI多光谱波段遥感影像，采用Landsat地表反射率编码(LaSRC)进行大气校正进而得到Landsat 8地表反射率数据[23]，对Landsat 8 OLI全色波段遥感影像进行辐射定标，消除传感器本身的误差。

水体容易和不透水层如旧建筑物、沥青等的低反照率类混淆，在提取城市不透水层信息前，首先进行水体掩膜处理。即利用改进的归一化差异水体指数(modified normalized difference water index，MNDWI)进行水体掩膜[24]。MNDWI(式中记为IMNDWI)计算公式为：

(1)

式中：ρg为绿光波段的地表反射率(在Landsat 8 OLI中为第3波段)，ρmir为中红外波段的地表反射率(在Landsat 8 OLI中为第6波段)。本研究中，经反复试验最终以0.032 882作为水体掩膜的阈值来剔除水体，剔除水体后的研究区假彩色图像如图1所示。

图1 剔除水体后的研究区假彩色图像Fig.1 The false colour composite image of the study area excluding water

1.3 城市不透水层提取

1.3.1 端元提取

假设某个多光谱波段像元为纯像元，则相对应的全色波段多个像元亮度值相同或者形似。相反，当该多光谱波段像元为混合像元时，相对应的全色波段像元亮度值会发生明显变化。因此，当多光谱波段像元为纯像元时，它对应的所有全色波段灰度值离散程度就低，反之就高。于是，引入一个空间变异系数(spatial coefficient of variation, SCV,式中记为ISCV)来表达某个多光谱波段像元内对应的全色波段像元灰度值的离散程度，其计算公式为：

(2)

SPPI的计算过程如下：①在整个图幅范围内对多光谱波段和全色波段选择300个控制点进行配准，使1个多光谱像元完全对应于全色波段的4个像元。②计算SCV并进行阈值分割，剔除SCV值较大的像元。在本研究中，根据经验尝试，选择0.05作为分割阈值，保留较接近端元的像元(空间变异阈值记为ISCVthreshold)。③对多光谱波段进行端元提取，为了保证所选端元的可靠性和有效性，首先采用最小噪声分离变换(minimum noise fraction，MNF)来降低原始多光谱影像数据噪声并消除原始波段之间的相关性；然后，将阈值分割后的SCV图像对MNF结果进行掩膜处理，剔除混合程度较大的像元。④再利用PPI获取纯像元，PPI是目前最常用的端元提取方法[17]，像元光谱纯度指数记为IPPI。经过掩膜后的SCV图像像元值作为权重系数，进一步扩大纯像元与混合像元的差异。对SCV图像进行归一化处理(SCVnormal，式中归一化处理后空间变异指数记为ISCVnormal)：

(3)

经过归一化处理后的SCVnormal图像与PPI相类似，其值越大则越接近纯像元。最后计算SPPI图像像元空间纯度指数值(ISPPI)：

ISPPI=IPPI·ISCVnormal·ISCVthreshold。

(4)

得到SPPI后，再利用ENVI软件结合研究区实际情况，将端元分为植被、高反照率不透水层、低反照率不透水层和土壤4类。

1.3.2 丰度估计

1) 线性光谱混合模型。线性光谱混合模型的原理是将像元在某一光谱波段的反射率假定为由构成像元的基本组分其反射率占像元面积比例为权重系数的线性组合。模型表达式为：

(5)

式中：xb为OLI影像第b波段的反射率，ρ为端元数，ai为端元i在像元内的面积权重，Si,b为端元i在波段b的反射率，еb为残差。可以利用最小二乘法进行丰度求解。

2) 混合调制匹配滤波。混合调制匹配滤波通过对目标波谱的最优线性检测与评估，进而反演得到每个像元中不同端元的丰度。MTMF利用MF进行丰度计算及背景压制后，利用MT识别假阳性并进行拒绝。模型简易表达式为：

IMFI=IMFV×IMNF。

(6)

式中：IMFI是丰度结果，IMFV是MF形成的投影向量，IMNF为上述最小噪声分离结果。

3) 双线性光谱混合模型。双线性混合模型考虑不同地物类别间的二次散射作用，并通过端元间的Hadamard积和非线性系数来表示非线性散射效应。模型表达式为：

(7)

式中：aj为端元j在像元内的面积权重，Sj,b为端元j在波段b上的反射率。该模型中各参数的意义和线性光谱混合模型中的参数意义相同。由于该模型为非线性的，一般无法直接求解，可采用高斯-牛顿法这一非线性最小二乘迭代法进行计算。

4)BP神经网络模型。人工神经网络具有自组织、自学习和高容错性等特点，故而能够解决复杂的非线性问题。最常见的是前向多层网络反向传播神经网络算法(BP神经网络算法)，该方法一般可以分为输入层、隐含层、输出层3层。输入层为波段数，在本研究中为7。隐含层神经元个数(M)可以采用式(8)进行计算：

(8)

式中：A为分类数，B为波段数。在本研究中，选用Sigmoid函数作为激活函数。

1.4 城市不透水层覆盖度计算

城市中的不透水层不是单一的土地覆盖类型，而是光谱变化较大的多种地物类型的总和。于是，采用Wu等[10]提出的模型，在高反照率、低反照率和不透水表层之间建立一种关系来刻画真实的不透水层。

IIS,b=Ilow·Ilow,b+Ihigh·Ihigh,b+eb。

(9)

式中：IIS,b、Ilow,b、Ihigh,b分别代表第b波段不透水层、低反照度和高反照度端元的反射率，Ilow、Ihigh分别为像元内部低反照度、高反照度地物的比率，еb为残差。

1.5 精度验证

利用研究区2018年Google Earth影像目视解译建筑物及道路结果来验证不同模型提取的不透水层丰度图的精度。在不透水层丰度图上，随机生成50个90 m×90 m的样本，计算每个样本的不透水层比率。同时，在Google Earth影像上选择对应区域进行目视解译，提取对应区域的不透水层面积比，再通过式(10)计算不透水层的窗口平均精度(IKp)。

(10)

式中：PiG为第i个低空间分辨率的影像像元中目标端元的百分比，即测量值；PiZ为第i个高空间分辨率的影像区域中目标端元的百分比，即参考值；n为验证样本数。

2 结果与分析

2.1 基于SCV的提取结果

对全色波段计算其空间变异系数SCV，如图2a所示，在空间内相对均质的区域，SCV值较低，这些区域较为纯净(如图中的紫金山、老山以及扬子江)。在多种地物交错的地带，SCV值较高，像元混合程度较大(如图中部的南京主城区)。随后利用阈值0.05对SCV图像进行掩膜，得到较为均质的区域，如图2b所示，同时经过掩膜后，计算量进一步减小。最后计算得到归一化SCV图像值，如图2c所示，当归一化SCV图像值越接近于1时，表明越接近于纯像元。

图2 SCV图像、SCV阈值掩膜图像和归一化SCV图像Fig.2 SCV image,threshold mask of SCV and normalized SCV image

2.2 端元提取结果

将SPPI指数的计算结果与PPI指数的计算结果相比较可以发现，纯净端元个数从PPI总共的48 727个点降到了SPPI总共的22 952个点。利用ENVI中的二维散点图提取植被、高反照率不透水层、低反照率不透水层和土壤4类地物端元，提取的二维散点图和各个地物类型的光谱曲线结果如图3、图4所示。

图3 MNF特征空间中使用所选端元绘制的散点图(显示目标端元的纯度)Fig.3 Scattergrams plotted in the MNF feature spaces by using the selected endmembers(showing the purity of the target endmembers)

图4 PPI与SPPI提取的4种端元反射波谱曲线Fig.4 The reflectance of four endmembers derived from PPI and SPPI

2.3 解混精度验证及不透水层丰度制图

通过对两种端元提取方法利用4种解混模型提取得到的不透水层丰度图与50个随机样本窗口进行验证，各个模型精度验证结果如表1所示。从表1可以看出:SPPI_BPNN的验证精度结果最高，为90.45%;SPPI_BMM的验证精度结果次之，为89.50%；而PPI_LSMM的验证精度结果最低，为80.62%。同时，基于SPPI的4种解混模型提取精度要优于各自对应的基于PPI的4种解混模型。此外，与仵宏基[3]利用各种不透水层指数提取结果相比较，SPPI提取结果精度也较优。

本研究选择精度更高的SPPI端元提取方法提取4类地物端元，再利用BP神经网络来对南京城市不透水层进行制图，结果如图5所示。从图5可以看出，在空间上，不透水层的范围由主城区向外辐射。在南京主城区，城市不透水层覆盖度比较高，均在0.8以上；而在农村及边远地区，城市不透水层覆盖度较低，均在0.4以下。同时，通过计算得到南京市总体的不透水层覆盖度为28.09%，与南京市统计局发布的2018年建筑用地占总面积的29.1%相比[22]，降低1.01%，这主要是由于在基于遥感光谱信号定义的不透水层与地面统计结果中的不透水层在概念理解上不完全一致所致。

表1 不同方法提取的不透水层覆盖度验证结果

图5 2018年南京市不透水层覆盖度图Fig.5 The abundance map of the extracted impervious surfaces in Nanjing in 2018

3 讨 论

1)以V-I-S模型为理论基础，将城市划分为植被、高反照率不透水层、低反照率不透水层和土壤4类，通过引入空间纯像元指数SPPI来提取端元。结果表明，空间纯像元指数SPPI与纯像元指数PPI相比，其提取结果更高。传统的仅基于多光谱波段提取端元的方式，由于多光谱图像离散的宽波段使得每个像元包含的光谱信息较少，因此很难反映像元内部不同地物间的精细差别[25-26]。而全色波段由于其较高的空间分辨率和更丰富的空间信息，可辅助多光谱波段提取端元，可将光谱信息和空间信息相结合，能更好地准确提取端元。

2)采用线性混合光谱模型、混合调制匹配滤波、双线性混合光谱模型和BP神经网络进行不透水层丰度提取，均可以有效提取城市不透水层。其中，BP神经网络提取精度优于另外3种算法。线性混合光谱模型和混合调制匹配滤波虽然能有效利用其物理模型，但是对于复杂的地物混合模型，其精度结果较低[10]。采用双线性混合光谱模型通过考虑二次反射以更好地模拟混合像元，但是通常会基于虚拟端元并将丰度、非线性参数等同时作为未知参数进行双线性混合光谱的拟合求解。在计算处理过程中加入参数进行约束会使得运算过于复杂而陷入局部极小，因此结果的精度以及算法的效率均难以保证[27]。BP神经网络由于其具有自组织、自学习和高容错性等特点，可以很好地进行混合像元分解。

3)在应用方面，本研究选择精度更高的SPPI端元提取方法提取4类地物端元后再利用BP神经网络来对南京城市不透水层进行制图，可为南京市土地利用变化、热岛效应、生态系统等方面的研究提供有力的支撑。除此之外，通过提取和分析城市不透水层基本信息与变化规律，能够准确地反映城市的分布与规模、判断城市化进程与发展方向、客观地评价城市生态环境，从而为政府有关部门客观评价城市状况、合理规划城市发展、制定相关政策法规、可持续地进行城市管理提供建议。

随着技术的不断进步和发展，更高空间分辨率和光谱分辨率的数据将是后期研究的重点。同时，各种精度更高的非线性混合分解模型也是未来研究的重点，以获得更高精度的不透水层覆盖度结果，为城市生态、规划、土地的合理利用等提供决策支持。