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

在中红外范围内, 地表特性的卫星数据较少, 而中分辨率成像光谱仪(moderate-resolution imaging spectroradiometer, MODIS)能提供较全波段范围的地表特性数据, 因此选取MODIS数据进行研究。 作为搭载于EOS(Earth Observing System)系列卫星上的关键传感器, MODIS因其较高的空间分辨率、 光谱分辨率、 时间分辨率优势, 在环境监测、 土地利用、 作物估产等各个领域中得到广泛应用^{[7, 8, 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

Table 1

表1(Table 1)

表1 MODIS反射率/发射率数据获取波段信息 Table 1 MODIS reflectivity/emissivity data acquisition band information 波段 波长范围/μm 中心波长 /μm 空间分辨 率/m 反射率/ 发射率 数据来源 Band7 2.105~2.155 2.13 500 反射率 MYD09A1 Band20 3.66~3.84 3.75 5 000 发射率 MOD11C3 Band22 3.93~3.99 3.96 5 000 发射率 MOD11C3 Band23 4.02~4.08 4.05 5 000 发射率 MOD11C3

表1 MODIS反射率/发射率数据获取波段信息 Table 1 MODIS reflectivity/emissivity data acquisition band information

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

An×d≈Vn×kHk×d(1)

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

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

minf(V, H)=‖V H‖F2s.t. V≥0, H≥0(2)

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

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

表2

Table 2

表2(Table 2)

表2 部分可用的非负矩阵分解Matlab软件包及约束条件和下载地址 Table 2 Some available non negative matrix decomposition MATLAB software packages, constraints and download addresses NMF algorithms Constraints Code links References Projected gradient methods for NMF (PGM-NMF) V≥ 0, H≥ 0 https://www.csie.ntu.edu.tw/~cjlin/nmf/ Lin, 2009 Block principal pivoting/Active set method (BPAS-NMF) V≥ 0, H≥ 0 https://www.cc.gatech.edu/~hpark/nmfsoftware.php Kim & Park, 2008, 2011 Robust collaborative NMF (Co-NMF) H≥ 0, ∑i=1ihi, l=1, ∀l https://www.lx.it.pt/~jun/demos.html Li et al., 2016 Group robust NMF for linear mixing model (G-rNMF) V≥ 0, H≥ 0, ∑i=1ihi, l=1, ∀l http://dobigeon.perso.enseeiht.fr/soft/rLMM.zip Fé votte & Dobigeon, 2015

表2 部分可用的非负矩阵分解Matlab软件包及约束条件和下载地址 Table 2 Some available non negative matrix decomposition MATLAB software packages, constraints and download addresses

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

r≈Vh(3)

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

rλ1︙rλd≈v1, 1…v1, k︙…︙vd, 1…vd, kh1︙hk(4)

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

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

r≈Vρh(5)

和

ρd1︙ρdm≈vd1, 1vd1, k︙…︙vdm, 1…vdm, kh1︙hk(6)

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

h≈Vρ+r(7)

式(7)中, Vρ+为矩阵V_ρ 的广义逆矩阵。

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

Fig.1

	Figure Option ViewDownloadNew Window
	图1 B/S架构Fig.1 B/S architecture

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

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

(1) JHU光谱数据预处理

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

图2

Fig.2

	Figure Option ViewDownloadNew Window
	图2 目标函数光谱响应函数(band7中心波长2.13 μm、 band20中心波长3.75 μm、 band22中心波长3.96 μm, band23中心波长4.05 μm)Fig.2 Objective function spectral response function (band7 center wavelength 2.13 μm, band20 center wavelength 3.75 μm, band22 center wavelength 3.96 μm, band23 center wavelength 4.05 μm)

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

射率数据。 基尔霍夫定律如式(8)所示

ρEmis+ρRefl=1(8)

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

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

εj=∑i=1nε(λi)f(λi)∑i=1nf(λi)(9)

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

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

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

图3

Fig.3

	Figure Option ViewDownloadNew Window
	图3 JHU波谱数据库提取的2~5 μm波段范围10 nm间隔的四种典型地表类型光谱数据集Fig.3 The spectral dataset of 4 typical surface types in range of 2~5 μm with the step of 10 nm extracted from JHU spectral data

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

图4

Fig.4

	Figure Option ViewDownloadNew Window
	图4 提取的2~5 μm波段范围10 nm间隔的用于NMF光谱重建的4个端元向量 (3) JHU 地物反射率光谱重建和结果验证Fig.4 The extracted spectra of 4 basis vectors for spectral reconstruction by NMF in the spectral range of 2 to 5 μm with the step of 10 nm

在获取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个波段数据可用的低秩病态的情况下, 可较好地满足光谱重建的精度要求。

图5

Fig.5

	Figure Option ViewDownloadNew Window
	图5 利用中心波长为2.13、 3.75、 3.96和4.05 μm的4个波段进行2~5 μm波段范围10 nm间隔的4个典型地表类型的反射率光谱重建结果Fig.5 The surface reflectance reconstructed results of 4 typical surface types in range of 2~5 μm with the step of 10 nm using the 4 multispectral reflectance bands centered at 2.13, 3.75, 3.96, 4.05 μm by the NMF method

图6

Fig.6

	Figure Option ViewDownloadNew Window
	图6 典型地表类型反射率光谱重建结果与真实结果的对比散点图Fig.6 The scatterplot of reconstructed reflectance versus the truth results for 4 typical surface types

在对光谱重建方法进行验证的基础上, 将其应用于广泛使用、 易于获取、 波段较全的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

Fig.7

	Figure Option ViewDownloadNew Window
	图7 MODIS全球月均地表反射率/发射率产品Fig.7 MODIS global monthly average surface reflectance/emissivity products

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

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

图8

Fig.8

	Figure Option ViewDownloadNew Window
	图8 全球范围内陆地5 km× 5 km分辨率的月均地表反射率Fig.8 Monthly average surface reflectances with 5 km× 5 km resolution on land in the world

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

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

图9

Fig.9

	Figure Option ViewDownloadNew Window
	图9 图层叠加展示Fig.9 Layer overlay display