地表温度在基于光谱混合分析提取城市不透水层中的应用

地表温度在基于光谱混合分析提取城市不透水层中的应用颉继珍;王秀琴
【摘要】文中利用Landsat-7 ETM+影像的多光谱和热红外波段进行光谱分析和地表温度反演,旨在开发一套提取城市不透水层的完整解决方案,特别是利用地表温度,从低反照率盖度影像中去除非不透水层信息,以便提高城市不透水层的评估精度.【期刊名称】《矿山测量》
【年(卷),期】2018(046)004
【总页数】7页(P5-11)
【关键词】不透水层;V-I-S模型;光谱混合分析;地表温度
【作者】颉继珍;王秀琴
【作者单位】甘肃省基础地理信息中心,甘肃兰州 730000;甘肃省地图院,甘肃兰州 730000
【正文语种】中文
【中图分类】TD171;P237
不透水层是城市土地覆盖的主要特征，如道路、屋顶等。

大量的不透水层不仅对城市的水循环(如地表径流和地下水)，而且对城市热岛效应有很大的影响。

城市不透水层被认为是城市化和环境质量的重要指标[1]。

不透水层研究在城市规划、城市景观生态学、城市环境管理、海绵城市建设等方面具有重要意义。

一项重要的工作是获取准确的不透水层分布状况，虽然野外测量和航空摄影测量可
取得可靠的结果，但遥感技术由于成本低廉、速度快等优点而更具吸引力，尤其是利用中等空间分辨率卫星影像提取不透水层。

由于城市景观的复杂性和混合像素在中低空间分辨率卫星影像中的存在，传统的分类方法(如最大似然法)总是产生不理想的结果。

混合像素已成为遥感数据在提取城市环境和生态参数应用方面的瓶颈[2]。

一个有用的城市景观概念模型，即植被-不透水层-土壤(V-I-S)模型被提出后，取得了很大的进展[3]。

在V-I-S模型中，城市土地覆盖被表示为植被、不透水层、土壤的线性组合。

基于 V-I-S模型开发了许多提取城市不透水层的方法，如决策树分类[4]、多元回归[5] 、人工神经网络[6] 、CART算法[7]。

此外，线性或非线性光谱混合分析(SMA)能够准确地进行子像元分解，被认为是一种有效的方法[8]。

线性光谱混合分析(LSMA)根据V-I-S 模型原理处理混合像元，即任何混合像元是像元内所有分量的光谱线性组合[9]。

LSMA操作简单，被广泛
用于获取不透水层信息——被看作一个或多个端元的线性组合[10]。

然而，SMA或LSMA的缺点体现在卫星影像的时相或阴影产生的不确定性。

一般情况下，通过低反照率和高反照率盖度影像之和计算不透水层百分比[11]，由于低反照率盖度影像包含其他类型的信息，如建筑物阴影、树冠阴影、水等，以及高反照率盖度影像中包含一些高亮度覆盖，如沙土或高亮度土壤等，因此不透水层被过度估计。

从低和高反照率盖度影像中删除这些类型的组分非常重要。

不透水层区域的温度比自然特征区高10°～12°[12]。

因此，利用地表温度的差异，可较好地去
除低反照率盖度影像不透水层中的阴影区域。

地表温度可从遥感影像 (如TM/ETM +) 热红外波段中提取[12]。

本文探讨了利用LSMA并结合地表温度从ETM+影像中提取不透水层的方法。

1 研究区及所用的数据
1.1 研究区
选择香港最繁华的九龙半岛及港岛市区为研究区(见图1)，该研究区内有不同的土
地利用类型，包括工业区、住宅区、商业区、中央商务区(CBD)、公园、水域、草地、森林等。

在维多利亚港两侧，密集的高楼可以投射出大量的阴影分布在CBD，这有助于测试本文所用方法的可行性。

图1 研究区位置
1.2 数据集及预处理
研究区ETM +影像于2010年9月在晴朗的天气条件下获得，该影像包括8个波段，其中6个波段为多光谱波段，空间分辨率为30 m，1 个热红外波段空间分辨率为60 m，1 个全色波段空间分辨率为15 m。

使用1: 20 000数字线划图(DLG)和 10 m 格网数字高程模型(DEM)进行了正射校正，采用UTM投影。

正射影像均方根误差(RMSE)小于0.3个像素。

6个多光谱波段和1 个热红外波段通过邻近像
元法被重采样为30 m×30 m像元大小，用于生成不同端元的盖度影像，热红外
波段用于反演地表温度。

研究区彩红外航空正射影像与ETM +影像同年获取，空间分辨率为0.21 m，采用UTM投影，被重采样为1 m×1 m像元大小，用于端元选择的辅助数据和不透水
层结果验证及精度评估的参考数据。

2 地表温度反演
地表温度是城市热岛效应的重要指标，在覆盖层(UCL)热岛中，不透水层和自然特征(如水体、草地、森林等)的地表温度存在10°～12°的明显差异[12]，本文将根据这一特征从低反照率盖度影像中去除阴影，并从ETM+影像的热红外波段(第六波段)中提取地表温度。

从ETM+6中提取地表温度有三种方法：辐射传输方程、单
窗算法、单通道法。

实验结果表明，当过程变量、平均大气温度、地面辐射率、大气透射率的值为中等误差水平时，单窗算法获得的结果更好[13]。

本文采用单窗算法提取研究区地表温度，计算公式为：
TS={a6(1-C6-D6)+[b6(1-C6-D6)+C6+D6]T6-D6Ta}/C6
(1)
C6=ε6τ6
(2)
C6=(1-τ6)[1+(1-ε6)τ6]
(3)
式中，TS表示地表温度，单位为K；a6和b6是常数，a6=-67.355 351，
b6=0.458 606；Ta表示有效平均大气温度；ε6表示地面辐射率；τ6表示大气透射率；T6表示传感器端亮度温度。

根据香港的地理位置，平均大气温度(Ta)可通过T0的线性函数来确定，即：
Ta=17.976 9+0.917 15×T0
(4)
式中，T0表示卫星过境时的近地表气温，T0=300.2 K(来自香港天文台)，因此Ta=293.3 K。

地面辐射率(ε6)可通过NDVI来确定，在该研究区，根据参考数据将土地覆盖分为四种类型，表1列出了相应的 NDVI 值和地面辐射率。

表1 NDVI 值和相应的地面辐射率[14]NDVI地表覆盖类型地面辐射率<-0.185水体0.995-0.185NDVI<0.157城市用地0.9700.157NDVI≤0.727自然特征1.009 4+0.047×ln(NDVI)*>0.727森林0.986
大气透射率(τ6)可用方程(5)用于计算。

τ6=1.031 412- 0.115 36 w
(5)
式中，w表示水汽含量(单位：g/cm-2)，卫星过境时研究区的水汽含量为1.89 g/cm-2(来自香港天文台)。

传感器端亮度温度(T6)可按以下公式计算。

T6=K2/ln(1+K1/Lλ)
(6)
式中，K1和K2为常量，对ETM+影像而言，K1= 666.09 mW/(cm2·sr·μm)，K2=1 282.71 K；Lλ为光谱辐射亮度，可按照以下公式从ETM+影像的热红外波段DN值转换而得。

Lλ=0.037 058 8×DN+3.2
(7)
获得所有参数后，可用式(1)计算得到地表温度，图2为研究区地表温度图像。

显然，不同的土地覆盖类型具有不同的地表温度。

土地覆盖地表温度从低到高依次为水体、山顶、植被、郊区、住宅区、工业和商业区。

因此，根据地表温度的差异，可将不透水层与其他土地覆盖类型区分开来。

图2 从ETM+影像的热红外波段提取的地表温度
3 不透水层估算
3.1 线性光谱混合分析
LSMA被定义为：某一波段像元的光谱反射率是像元基本组件(端元)反射率的线性组合[9]。

LSMA用以下数学公式表示：
(8)
式中，i表示一个波段的像元数；λ表示波段数(= 1, 2,...,m)；Riλ表示第λ波段第i像元的光谱反射率；k表示端元数(k= 1, 2,..., n)；fki表示第i像元内端元k的占比；Ckλ表示第λ波段内端元k的光谱反射率；εiλ表示第λ波段的残差。

在式(8)中，通过将影像DN值转化为星上反射率可得到Riλ，定标参数以及传感器的增益值gains和偏置值bias，可直接从ETM+影像元数据中获得；当端元确定后Ckλ即为已知；fki为未知变量，通过求fki能得到盖度影像。

要解fki，端元
选择是关键一步。

选择端元应满足两个条件：(1)应具有代表性和独立性；(2)端元的数量不应大于波段数，即n≤m。

此外，要使fki成为一个有意义的解，还应满足式(9)和式(10)：
(9)
0≤fki≤1
(10)
根据已知和未知的条件，采用约束条件下最小二乘法求解是可行的，并可用均方根误差 (RMSE) 评估其精度：
(11)
3.2 端元选择
根据研究区的实际土地利用或土地覆盖类型，通过参考高分辨率彩红外 DOM，最初确定了高反照率、低反照率、植被、土壤四种类型的端元。

利用最小噪声分离(MNF)变换进行数据降维，降低 ETM +影像6个波段之间的数据冗余和相关性。

在MNF变换前，对维多利亚港大片水域进行了掩膜处理。

按照信噪比对结果分量进行排序，MNF主分量用于生成散点图和计算像元纯度指数(PPI)。

影像端元可在每两个MNF分量(MNF1-MNF2、MNF1-MNF3、MNF2-MNF3)的特征空间散点图中选择，理论上，在两个MNF分量的特征空间中，所有像元都应位于三角形中，其顶点即为端元[15]。

然而在实际中，散点图通常是不规则的三角形，并且在顶点附近有很多点，为了提高端元选择的准确性，本文采用 PPI 对这些点进行滤波。

最终通过散点图和PPI确定了图像端元，他们与初始端元具有相似的光谱反射率特征(详见图3)。

图3 六个波段中四个所选端元的反射率
3.3 不透水层影像的生成
根据V-I-S模型，城市景观是植被、不透水层、土壤的线性组合，然而，作为人文景观，不透水层具有光谱可变性，因此不能通过一个端元来估算。

研究发现，不透水层和低反照率、高反照率这两个端元具有密切关系，不透水层可能位于或接近低反照率和高反照率端元在特征空间的连线[15]，这种关系可被表达为：
Rimp,b=flowRlow,b+fhighRhigh,b+eb
(12)
式中，Rimp,b、Rlow,b、Rhigh,b分别表示b波段中不透水层、低反照率、高反照率地物的反射光谱；flow和fhigh分别表示低反照率和高反照率盖度；eb表示残差。

双端元线性混合模型在哥伦布CBD的应用取得了令人满意的结果[15]。

本文也在
研究区的CBD区域对双端元线性混合模型进行了测试，并取得了一个可靠的结果，其平均RMSE小于0.02。

基于上述模型，可以通过低反照率和高反照率盖度相加
来估算不透水层。

然而，不透水层的估算精度受一些低反射率地物(如水、阴影)和高反射率地物(如云、沙土)的影响。

特别是低反射率地物，其组成往往比较复杂，可能包含不透水层、水体、建筑物阴影、森林中树冠阴影、潮湿的土壤等，直接将低反照率和高反照率盖度相加会导致过度估算不透水层。

因此，在根据约束线性混合模型估算不透水层前，从低反照率和高反照率盖度中分别剔除这些低反射率和高反射率地物显得尤为重要。

3.4 不透水层与其他类型地物的分离
为了提高不透水层估算的精度，本文采用一些有效的方法将不透水层与其它地物类型进行分离。

就沙土和干燥土壤而言，容易发现他们在土壤和高反照率盖度影像中具有较高的反射特性。

因此，通过在土壤盖度影像中确定的阈值可以从高反照率盖度影像中去除沙土、干燥土壤等。

在土壤中，沙土和干燥土壤的反射率高于其他类
型的土壤，阈值fS可利用彩红外航空影像确定。

根据两个条件对高反照率盖度影
像的像元重新赋值，即：如果土壤盖度影像的像元值大于fS，则赋值“0”，否则，保持不变。

对于低反照率盖度影像中的低反射率地物，能够利用DLG 数据对水体进行掩膜，但剔除阴影是一个复杂的过程，基于不透水层与自然特征之间存在温差的这一基础，本文提出了区分不透水层和高层建筑或树冠阴影下低反照率自然特征的方法，这可通过低反照率盖度影像和地表温度影像的特征空间来说明(见图4)。

图4 低反照率盖度影像和地表温度影像的特征空间
在散点图中，A区代表草地、森林等自然特征，在低反照率盖度影像和地表温度影像中像素值较小；B区代表低反照率盖度影像中高像素值的水体，地表温度影像中的低像素值；C区代表在地表温度影像中具有中像素值的住宅区，其像素值几乎覆盖整个低反照率盖度影像的值域；D区代表CBD、工业及交通运输区，在地表温
度影像中具有非常高的像素值，且几乎覆盖整个低反照率盖度影像的值域。

很明显，人工地物(如建筑物、道路、停车场)比自然地物(如草地、森林、水体)的温度更高。

在阴影区域，自然地物的温度一定低于人工地物。

因此，温度是区分阴影区域不透水层和自然地物的良好指标。

本文通过低反照率盖度影像、地表温度影像以及航空影像的特征空间，确定了从低反照率盖度影像中去除阴影的温度阈值。

具体做法是根据两个条件给低反照率盖度影像重新赋值，即：如果地表温度影像的像元值小于fT，则赋值“0”，否则，保持不变。

从高反照率盖度影像和低反照率盖度影像中分别去除高反射率地物(如沙土、干燥
土壤)和低反射率地物(如水体和阴影)后，就可通过将改正后的高反照率和低反照
率盖度影像相加来估算不透水层。

4 结果和分析
4.1 LSMA结果
利用约束条件下最小二乘法获得高反照率、低反照率、土壤和植被四个端元盖度影像(见图5)。

LSMA精度良好，平均RMSE为 0.004 7，最大误差为 0.088，大多
数RMSE小于0.018，较大的误差基本上集中在一些高反射区，特别是一些具有
特殊材料屋顶的建筑物(如体育馆)，这表明一些特殊材料不能用所选的四个端元的线性组合很好地表达。

用盖度影像的像素值表示像元中每个端元所占的比例，其取值范围为0到 1，像素越亮，端元占比越高。

高反照率盖度影像在商业区及一些特殊建筑物(如位于湾仔
的香港会展中心)具有较高的像素值。

土壤主要分布在居住区、郊区和植被稀少的
山顶。

在高反照率和土壤盖度影像中，有些像元具有较高的像素值，通过与参考数据比较，证明他们是砂土或干燥土壤。

植被盖度影像在公园、绿地和有森林及浓密植被的郊区具有较高的像素值。

低反照率盖度影像的空间分布比较复杂且难以描述，因为低反照率组成比较复杂，包括道路、跑道、沥青表面材料、海港、水体、湿润土壤、建筑物或树冠阴影等。

图5 四个端元盖度影像
4.2 不透水层提取结果
图6显示了两种不同方法获得的结果，在采用改进方法前的不透水层影像中(图
6a)，发现部分局部区域的不透水层分布不合理且比例过高，如维多利亚港两侧的
海岸线和港口、高层建筑之间的绿地、郊区的林地、沙土、干燥裸露地表等，这是因为一些低反射率地物和一些高反射率地物被分别包含在低反照率和高反照率盖度影像中。

采用改进的方法后，估算的不透水层最终结果(图6b)是令人满意的，这
说明与不透水层不相关的其他覆盖类型被成功分离。

最终的不透水层影像空间分布合理，不透水层占比由高到低依次为：CBD、商业区、工业区、高密度住宅区、
低密度住宅区、郊区。

图6 用不同方法提取的不透水层影像
4.3 精度评定
利用彩红外航空影像进行精度评估，采用随机抽样法获得了不透水层精度评定的参考数据，随机选取了72个样本，大小为5像素×5像素(150 m×150 m)，以减少航空影像和不透水层影像之间几何误差的影响，用式(13)计算RMSE，评价不透水层估算精度。

(13)
式中，表示第i个样本不透水层估算盖度；Vi表示从航空影像获得的实际不透水层盖度；N是样本总数，N=72。

评估结果表明，不透水层估算精度很高，所有样本RSME为10.26%，详细的精度分析如图7所示。

使用地表温度前，精度情况如下：(1)RSME为 18.04%；(2)在72个样本中，有
11个样本残差超过±0.3，最大残差为0.465，有56个样本残差在±0.2之间，占
样本总数的约78%；(3)当实际盖度小于0.6时，几乎所有的不透水层盖度被高估，当实际盖度介于0.6 和0.7之间时，估算基本准确，当实际盖度大于0.8时，又被低估。

使用地表温度后，精度情况如下：(1)RSME为10.26%，整体精度提高了约43%；(2)在所有样本中，只有1个样本残差值超过了，最大的残差为-0.349，67个样本残差在之间，占样本总数的 93%；(3)实际盖度小于0.3 的不透水层略有高估，0.3 到0.7之间时，估算准确，大于0.7 时，略被低估，这一趋势与之前的类似研究是一致的。

图7 基于不同方法估算的不透水层精度评定结果
5 结论
本文基于LSMA模型，成功地从ETM+影像中提取了城市不透水层信息。

结果表
明，作为城市生态环境的重要指标，不透水层可利用遥感技术准确地估算，相比野外测量或航空摄影测量更廉价、更快捷、更高效。

本文提出了利用中空间分辨率卫星影像的多光谱和热红外波段提取城市不透水层的完整解决方案，利用高分辨率航空影像进行了精度评估，总体RMSE为10.26%，结果可靠。

通过分析估算过程和结果，可得出以下结论：(1)在LSMA中，端元选择尤为关键，直接影响不透水层的估算精度。

在本研究中，PPI的引入能够克服从具有不规则三角形形状的散点图中选择端元的不确定性，从而提高影像端元的精度和可靠性；(2)利用地表温度可以有效地从低反照率盖度影像中去除低反射性地物；(3)土壤盖度
影像可作为阈值，将具有高反射率的非不透水层从高反射率盖度影像中分离出来。

结果表明，这些方法是有效的，可大大提高不透水层的估算精度。

参考文献：
【相关文献】
[1] Arnild, C.L.Jr.,Gibbons,et al. Impervious surface coverage: the emergence of a key environmental indicator[J].Journal of the American Planning Association, 1996(62): 243-258.
[2] Cracknell, A.P.. Synergy in remote sensing- what’s in a pixel[J].Interna tional Journal of Remote Sensing, 1998(19): 2025-2047.
[3] Ridd, M. K.. Exploring a V-I-S (Vegetation-Impervious Surface-Soil) model for urban ecosystem analysis through remote sensing: Comparative anatomy for cities[J]. International Journal of Remote Sensing, 1995(16): 2165-2185.
[4] Dougherty, M., Dymond, R.L.,et al.. Evaluation of impervious surface estimates in a rapidly urbanizing watershed[J]. Photogrammetric Engineering and Remote Sensing,
2004(70):1275-1284.
[5] Bauer, M.E., Heiner, N.J., Doyle, J.K.,et al. Impervious surface mapping and change monitoring using Landsat remote sensing[C]. ASPRS Annual Conference Proceedings, Denver, Colorado, 2004.
[6] Flanagan, M.,Civco, D.L.. Subpixel impervious surface mapping[C]. ASPRS Annual Conference Proceedings, St. Louis, Missouri,2001.
[7] Yang, L., Huang, C., Homer, C.G., et al..An approach for mapping large-scale impervious surface: synergistic use of Landsat-7 ETM+ and high spatial resolution imagery[J]. Canadian Journal of Remote Sensing, 2003(29): 230-240.
[8] Roberts, D.A., Batista, et al.. Change identification using multitemporal spectral mixture analysis: applications in eastern Amazonia[M]. In: Remote Sensing Change Detection: Environmental Monitoring Methods and Applications, Ann Arbor Press, Ann Arbor, MI, 1998.
[9] Adams, J.B., Sabol,et al.. Classification of multispectral images based on fractions of endmembers: application to land cover change in the Brazilian Amazon[J]. Remote Sensing of Environment, 1995(52): 137-154.
[10] Rashed, T., Weeks, J., Roberts, et al.. Measuring the Physical Compositions of Urban Morphology using Multiple Endmember Spectral Mixture Models[J]. Photogrammetric Engineering & Remote Sensing,2003(69): 1011-1020.
[11] 岳文泽,吴次芳.基于混合光谱分解的城市不透水面分布估算[J].遥感学报, 2007(6): 914-922.
[12] Weng, Q., Lu, D.,J. Schubring. Estimation of Land Surface Temperature- Vegetation Abundance Relationship for Urban Heat Island Studies[J]. Remote Sensing of Environment, 2004, 89(4): 467-483.
[13] Sobrino, J.A., Jime′nez-Mun～oza, J.C., et al.. Land surface temperature retrieval from LANDSAT TM 5[J].Remote Sensing of Environment, 2004(90): 434-440.
[14] Van de Griend AA, Owe M.. On the relationship between thermal emissivity and the normalized difference vegetation index for natural surfaces[J]. Int J Remote Sens,
1993,14(6):1119-1131.
[15] Wu, C. ，Murray, A.T.. Estimating impervious surface distribution by spectral mixture analysis[J]. Remote Sensing of Environment, 2003(84): 493-505.。