基于非负矩阵分解的中红外地表特性光谱重建方法

李殷娜, 李正强*, 郑 杨, 侯伟真, 徐文斌, 3, 马 䶮,樊 程, 葛邦宇, 姚 前, 史 正

1. 中国科学院空天信息创新研究院, 国家环境保护卫星遥感重点实验室, 北京 100101 2. 中国科学院大学, 北京 100049 3. 北京环境特性研究所光学辐射重点实验室, 北京 100854

引 言

近年来, 高光谱技术取得了重大突破, 并成功支撑气候、 生态、 环境等多领域应用。 在军事领域, 对高光谱数据的需求更是与日俱增, 如地表复杂物理参数估计及具有精细光谱特征的视觉相似材料的识别等[1]。 但高光谱数据相对多光谱数据而言较难获取, 且数据量巨大, 故国内外多数学者都尝试通过光谱重建的方式来获取高光谱数据。 部分学者利用混合像元分解或冠层辐射传输模型实现了对地物光谱的重建, 如赵永光等[2]提出一种基于冠层辐射传输物理机理并充分考虑像元异质性的方法, 进行地表反射率光谱重建, 但该类方法易受混合像元影响, 对地物的理想化假设会造成较大模拟误差。 矩阵分解是实现大规模数据处理与分析的一种有效方式, Hou等[3-4]将主成分分析(principal components analysis, PCA)分解法应用于光谱压缩和重建的研究中并取得较好的结果, 但PCA得到的主成分向量不可避免地会有负数的存在, 丧失了直观的物理意义, 且该方法针对病态问题的普适性不强。 相较于PCA分解法, 非负矩阵分解(nonnegative matrix factorization, NMF)法则采用了非负基向量的加权组合, 非负的数值结果更为直观并具有可解释的物理意义, 可以取代PCA分解等方法用来重建多光谱或高光谱地表反射率, 如官铮等[5]基于光谱重建约束的NMF法, 实现了高光谱图像与全色图像的融合。

现有研究主要是利用给定光谱范围内所有波段的信息进行光谱重建, 对只利用其中少量典型波段信息进行给定光谱范围低秩病态的光谱重建方法还少有研究。 本工作基于约翰霍普金斯大学(Johns Hopkins University, JHU)地物波谱库和多光谱卫星数据, 提出了一种利用非负矩阵分解法重建中红外全波段地表反射率/发射率的方法, 并选取了4种典型地物类型(土壤、 植被、 人造材料和岩石)进行方法验证。 针对短波红外和中红外波段范围内的典型地物光谱重建方法进行研究, 可为后续的全波段高光谱和多光谱遥感应用提供关键支撑。

目前对光谱数据可视化的研究还处于传统制图阶段, 往往只能进行二维平面展示, 已无法满足多源数据三维可视化表达的需求。 由于网络的发展与信息量的激增, 网络地理信息系统(web geographic information system, WebGIS)逐渐成为GIS的一个主流发展方向, 并且随着对三维可视化的需求日益增长, 三维WebGIS应运而生。 针对地表反射率重建结果可视化进行研究, 基于WebGIS技术开发了一套以卫星产品为中心, 将卫星底图、 地形数据与光谱重建结果等集成展示的二三维一体化可视化系统, 使空间关系能得到更明确的展示, 表现力远超二维地理信息系统, 扩展了三维信息的应用深度, 为卫星产品的展示与研究提供新方案。

1 实验部分

1.1 数据介绍

采用的光谱数据为约翰霍普金斯大学(JHU)光谱数据集, 该光谱库提供了15个子库, 包括岩石、 矿物、 月球土壤、 陆地土壤、 人造材料、 陨石、 植被、 冰雪等, 其中大部分测量结果来自JHU或美国地质调查局[6]。

在中红外范围内, 地表特性的卫星数据较少, 而中分辨率成像光谱仪(moderate-resolution imaging spectroradiometer, MODIS)能提供较全波段范围的地表特性数据, 因此选取MODIS数据进行研究。 作为搭载于EOS(Earth Observing System)系列卫星上的关键传感器, MODIS因其较高的空间分辨率、 光谱分辨率、 时间分辨率优势, 在环境监测、 土地利用、 作物估产等各个领域中得到广泛应用[7-9]。 本工作所采用的产品是MODIS的地表反射率产品MYD09A1和发射率产品MOD11C3。 其中, MYD09A1提供了根据大气条件(例如气体、 气溶胶和瑞利散射)校正的AquaMODIS波段1～7的地表光谱反射率估计值; MOD11C3在0.05纬度/经度气候模拟网格(CMG)中提供了月度地表温度和发射率(LST&E)值, 且LST&E值是通过对MOD11C1每日文件相应月份的值进行合成和平均得出的。

在2.0～5.0 μm范围的短波红外和中红外波段, 由于MODIS只有2.105～2.155 μm的地表反射率产品, 以及3.66～3.84, 3.93～3.99和4.02～4.08 μm的地表发射率产品, 所以选取对应的第7、 20、 22、 23波段作为模拟2.0～5.0 μm目标通道的数据源通道(如表1所示), 其中, 中心波长2.13 μm为地表反射率产品的通道, 中心波长3.75、 3.96、 4.05 μm为地表发射率产品的通道。

表1 MODIS反射率/发射率数据获取波段信息

1.2 非负矩阵分解方法

矩阵分解是一种有效降维方法, 可实现大规模数据处理与分析, 通过迭代分解的方式, 近似地将原始矩阵分解为2个较低维数矩阵的乘积, 进而提取高光谱数据特征[10]。 非负矩阵分解(NMF)通过对基向量与对应的系数增加非负约束, 以保证分解结果均为正值, 在高光谱图像混合像元分解上取得应用[11-13]; 基矩阵利用系数矩阵作为权重, 重构原始矩阵。 对于一个n×d非负矩阵An×m(如:n表示光谱数据集的条数,d表示每条光谱的波段个数), 存在一个非负矩阵Vn×k和Hn×k, 有

An×d≈Vn×kHk×d

(1)

式(1)中,V为基矩阵;H为系数矩阵;k为分解得到矩阵的秩, 且k≪min(n,d);A中的列向量可解释为对基矩阵V中所有列向量的线性组合, 且权重系数为系数矩阵V中对应列向量中的元素。

采用欧式(Euclidean)距离作为目标函数衡量矩阵分解前后的逼近程度, 非负矩阵分解问题被表示成如式(2)最优化问题

s.t.V≥0,H≥0

(2)

式(2)中, ‖‖F表示矩阵Frobenius范数; s.t.为subjectto的缩写, 对应约束条件。

表2列出了4种典型的NMF算法和非负性约束, 经过多次实验, 最终采用了较为稳定的BPAS-NMF方法。

表2 部分可用的非负矩阵分解Matlab软件包及约束条件和下载地址

1.3 基于挑选波段的病态光谱重建方法

与基于PCA的光谱重建类似, 地表反射率的光谱向量可以表示为

r≈Vh

(3)

写成具体的向量和矩阵形式为

(4)

式(4)中, 矩阵V的列向量即为提取的端元向量, 可通过对光谱数据集的非负矩阵分解处理得到;h为对应的丰度权重系数。

对于2.0～5.0 μm光谱范围的地表反射率/发射率数据集, 基于已提取的端元向量, 若利用卫星遥感常用的几个典型波段的地表反射率/发射率信息来重建连续波段的地表反射率/发射率, 就必须先获取对应的丰度权重系数向量h的值。 利用ρ来表示从r中挑选的m个波段地表反射率组成的向量(m≪d), 已知的几个波段对应的下标依次用d1,d2, …,dm来表示,Vρ表示V这几个波段对应的端元子矩阵, 目的就是利用ρ来确定h的近似值, 进而进行光谱重建。 这样,ρ可以表示为

r≈Vρh

(5)

和

(6)

基于式(6), 丰度系数向量h可以利用一个最小二乘的形式进行求解

(7)

1.4 二三维一体化可视化方法

WebGIS是网络与地理信息系统的结合, 作为GIS发展的一个主流模式, 已广泛应用于灾害响应、 室内导航、 资源管理、 环境监测等各个领域, 并逐渐与三维可视化技术相结合。

WebGL是HTML5标准的一部分, 是一种新的三维Web绘图标准, 遵从HTML5标准的浏览器自带对WebGL的支持, 降低了用户的使用成本, 目前火狐、 谷歌、 欧朋等浏览器均实现了对WebGL的支持, 用户只需安装浏览器, 即可在浏览器页面中实现强大的三维功能, 为用户使用三维技术创造了极大的便利。 WebGL应用在包含HTML、 CSS与JavaScript文件的传统Web应用基础上, 还包括了三维模型数据及着色器的源代码, 用JavaScript调用WebGLAPI即可快速实现Web三维功能。

主流的三维WebGIS产品包括Cesium、 ArcGIS API for Javascript和Super Map iClient 3D for WebGL等, 考虑到开发成本和需求, 采用完全开源的、 基于WebGL的Cesium框架。 Cesium是一个成本低、 开发简单、 支持多种数据可视化方法、 可实现二三维一体化, 面向三维地球的JavaScript库, 可在不同平台、 不同浏览器展示三维地球, 不需任何插件支持。 Cesium提供基于JavaScript语言的开发包, 可实现全球高精度的地形和影像服务, 且支持多种场景模式(2D, 3.5D以及3D场景), 使用者可迅速创建零插件的虚拟地球网络应用, 实现二三维一体化[14]。

本研究采用的B/S架构, 本质是将浏览器作为客户端, 于服务器端进行交互, 具体运作模式如图1所示。

图1 B/S架构

浏览器端通过HTM5和CSS实现界面及交互设计, 包括图层切换、 拖拽、 视角切换等, 并进行场景渲染。 服务器端使用JS语言, 利用Cesium框架提供的接口, 调用已发布的瓦片地图, 并进行数据读取、 分析及三维可视化。

2 结果与讨论

2.1 基于JHU光谱数据集的NMF中红外光谱重建

为论证本文所提出的光谱重建方法的可行性, 首先将该方法应用于地物波谱库数据集的重建上, 基于地物波谱库中典型地物的波段关系, 建立2.5～5.0 μm中红外波段的光谱重建模型, 并对重建结果进行了验证。

(1) JHU光谱数据预处理

MODIS传感器数据源通道的光谱响应函数如图2所示。

图2 目标函数光谱响应函数(band7中心波长2.13 μm、 band20中心波长3.75 μm、 band22中心波长3.96 μm, band23中心波长4.05 μm)

由于选取的各类地物样本光谱分辨率不尽相同, 因此需要对地物光谱进行光谱间隔插值, 采用线性插值将地物光谱间隔统一为10 nm, 进而基于光谱响应函数对地物光谱库中地物光谱数据进行重采样。 光谱响应函数是波长的函数, 表示传感器在每个波长处所接收的辐射能量与入射的辐射能量的比值。 光谱响应函数表征的是不同波长处信号响应的差异, 而传感器光谱响应函数在真实连续光谱上进行加权平均采样后得到光学遥感图像。 在构建转换模型时, 会用到发射率数据。 而地物光谱库中的数据均为地物反射率数据, 因此可以通过基尔霍夫定律对发射率数据进行转换, 从而得到反射率数据。 基尔霍夫定律如式(8)所示

ρEmis+ρRefl=1

(8)

式(8)中,ρEmis为发射率,ρRefl为反射率。 因此, 利用该公式可实现反射率数据和发射率数据的转换。

基于已建立好的地物光谱库样本信息, 利用MODIS传感器的光谱响应函数, 对数据源通道的地物反射率数据, 根据等效计算公式进行重采样, 等效计算公式如式(9)所示

(9)

式(9)中,λ为波长,ε(λi)和f(λi)分别是波长为λi时的辐射数据和响应函数,n为通道j范围内的辐射数据或响应函数值的个数。

(2) JHU地物反射率光谱数据集端元光谱向量的提取

利用线性插值将JHU波谱数据库2.0～5.0 μm波段范围反射率结果重采样到10 nm光谱间隔的301个波段, 具体如图3所示。

图3 JHU波谱数据库提取的2～5 μm波段范围10 nm间隔的四种典型地表类型光谱数据集

在此基础上, 进行非负矩阵分解处理, 提取得到对应的端元向量, 如图4所示, 这4条端元向量即为用于非负矩阵光谱重建可通用的基准光谱向量。

图4 提取的2～5 μm波段范围10 nm间隔的用于NMF光谱重建的4个端元向量

(3) JHU 地物反射率光谱重建和结果验证

在获取4条端元向量光谱曲线的基础上, 从JHU地表反射率光谱数据集中提取MODIS卫星中心波长为2.13、 3.75、 3.96和4.05 μm这4个波段的地表反射率的子数据集, 并计算对应的权重系数向量结果, 进行2.0～5.0 μm光谱范围内的全波段的反射率光谱重建。 利用MODIS的4个波段反射率结果进行全波段光谱重建的结果如图5所示, 对应的散点图如图6所示, 对应的平均绝对误差为0.01, 平均相对误差为10%。 对应在只有MODIS卫星4个波段数据可用的低秩病态的情况下, 可较好地满足光谱重建的精度要求。

2.2 基于卫星观测数据的全球中红外任意波段光谱重建

在对光谱重建方法进行验证的基础上, 将其应用于广泛使用、 易于获取、 波段较全的MODIS反射率/发射率产品的光谱重建。 利用MOEIS的地表反射率和发射率产品, 推导计算到任意中红外波段的地表反射率。

(1) MODIS卫星观测数据预处理

MODIS获取的地表反射率数据和发射率数据存在数据无效或缺失的问题, 针对这一问题, 一方面基于往年同期、 同年前后时间的反射率补丁数据对原始反射率数据中的无效像元进行修复; 另一方面利用原始反射率数据中获取到的海表反射率数据, 计算其均值和标准差, 并基于正态分布模型对海洋范围缺失数据进行修复。 图7展示了经过预处理得到的MODIS全球月均地表反射率/发射率产品, 其中图7(a)为中心波长2.13 μm的空间分辨率500 m的地表反射率产品的通道, 图7(b—c)分别对应中心波长3.75、 3.96、 4.05 μm的空间分辨率5 km的地表发射率产品。

图7 MODIS全球月均地表反射率/发射率产品

(2) 全球中红外波段光谱重建

在处理好MODIS短波红外和中红外4个波段的全球月均地表反射率/发射率产品的基础上, 利用基于NMF和光谱数据集所提取的2.0～5.0 μm用于反射率光谱重建的端元向量, 可以计算得到每个像元对应的权重系数向量, 进而可进行任意波段的光谱重建, 得到的全球范围内陆地5 km×5 km分辨率的月均地表反射率如图8所示, 通过全球陆地地表反射率的光谱重建, 可为后续中红外相关波段的背景辐射研究提供重要支撑。

图8 全球范围内陆地5 km×5 km分辨率的月均地表反射率

2.3 基于WebGIS的二三维一体化可视化

基于Cesium框架, 采用B/S架构, 搭建了二三维一体化可视化系统, 将卫星底图、 地形数据与光谱重建结果等集成展示。 实现了图层切换、 拖拽、 视角切换等地图控制功能及瓦片地图调用, 数据读取、 分析及三维可视化等数据处理功能。

将全球中红外光谱重建结果进行投影转换, 建立符号系统并进行影像切片处理发布服务, 即可在系统中调用, 将地形服务进行夸张处理后与之叠加, 可直观展现出地表反射率与地形之间的关系如图9所示。 图9(a)展示了图层叠加后地表反射率重建结果在三维地球上的展示效果, 图9(b)近距离展示了地表反射率图层覆盖于地形上的效果, 可以直观看出, 该区域内高程与地表反射率总体呈负相关。

图9 图层叠加展示

3 结 论

提出了一种基于地物波谱库和多光谱卫星数据, 利用非负矩阵分解方法重建高光谱地表反射率的方法, 实现了2.5～5.0 μm光谱范围内10 nm间隔的地表反射率光谱重建, 同时为综合评价该光谱重建方法, 选取了4种典型地物类型(土壤、 植被、 人造材料和岩石)进行方法验证。 研究结果表明:

(1)本文提出的NMF光谱重建方法, 相较于PCA等方法, 增加了非负约束, 具有可解释的物理意义, 得到的重建结果较为稳定和精确。 利用本方法得到的典型地表类型反射率光谱重建结果与真实结果的平均绝对误差为0.01, 平均相对误差为10%, 误差较小, 有效论证了该方法的可行性, 对应在只有MODIS卫星4个波段数据可用的低秩病态的情况下, 可较好地满足光谱重建的精度要求。

(2)基于本文提出的方法, 利用MODIS的地表反射率和发射率产品, 可获得全球范围内地表反射率或发射率高光谱数据, 能够在一定程度上弥补卫星手段难以获取红外高光谱地表反射率数据的问题, 并且该方法可进一步扩展到可见光光谱范围内的地表反射率重建, 为卫星遥感的相应应用提供支撑。

(3)为将地表反射率与对应地表覆盖和地形更直观地进行可视化分析, 基于Cesium开发库进行设计, 采用B/S架构, 搭建了二三维一体化可视化系统, 将卫星底图、 地形数据与光谱重建结果等集成展示, 能直观进行多因素分析, 扩展了三维信息应用的深度, 为卫星产品的展示与验证提供新方案。