基于非負(fù)矩陣分解的中紅外地表特性光譜重建方法

李殷娜, 李正強*, 鄭 楊, 侯偉真, 徐文斌, 3, 馬 ?,樊 程, 葛邦宇, 姚 前, 史 正

1. 中國科學(xué)院空天信息創(chuàng)新研究院, 國家環(huán)境保護(hù)衛(wèi)星遙感重點實驗室, 北京 100101 2. 中國科學(xué)院大學(xué), 北京 100049 3. 北京環(huán)境特性研究所光學(xué)輻射重點實驗室, 北京 100854

引 言

近年來, 高光譜技術(shù)取得了重大突破, 并成功支撐氣候、 生態(tài)、 環(huán)境等多領(lǐng)域應(yīng)用。 在軍事領(lǐng)域, 對高光譜數(shù)據(jù)的需求更是與日俱增, 如地表復(fù)雜物理參數(shù)估計及具有精細(xì)光譜特征的視覺相似材料的識別等[1]。 但高光譜數(shù)據(jù)相對多光譜數(shù)據(jù)而言較難獲取, 且數(shù)據(jù)量巨大, 故國內(nèi)外多數(shù)學(xué)者都嘗試通過光譜重建的方式來獲取高光譜數(shù)據(jù)。 部分學(xué)者利用混合像元分解或冠層輻射傳輸模型實現(xiàn)了對地物光譜的重建, 如趙永光等[2]提出一種基于冠層輻射傳輸物理機(jī)理并充分考慮像元異質(zhì)性的方法, 進(jìn)行地表反射率光譜重建, 但該類方法易受混合像元影響, 對地物的理想化假設(shè)會造成較大模擬誤差。 矩陣分解是實現(xiàn)大規(guī)模數(shù)據(jù)處理與分析的一種有效方式, Hou等[3-4]將主成分分析(principal components analysis, PCA)分解法應(yīng)用于光譜壓縮和重建的研究中并取得較好的結(jié)果, 但PCA得到的主成分向量不可避免地會有負(fù)數(shù)的存在, 喪失了直觀的物理意義, 且該方法針對病態(tài)問題的普適性不強。 相較于PCA分解法, 非負(fù)矩陣分解(nonnegative matrix factorization, NMF)法則采用了非負(fù)基向量的加權(quán)組合, 非負(fù)的數(shù)值結(jié)果更為直觀并具有可解釋的物理意義, 可以取代PCA分解等方法用來重建多光譜或高光譜地表反射率, 如官錚等[5]基于光譜重建約束的NMF法, 實現(xiàn)了高光譜圖像與全色圖像的融合。

現(xiàn)有研究主要是利用給定光譜范圍內(nèi)所有波段的信息進(jìn)行光譜重建, 對只利用其中少量典型波段信息進(jìn)行給定光譜范圍低秩病態(tài)的光譜重建方法還少有研究。 本工作基于約翰霍普金斯大學(xué)(Johns Hopkins University, JHU)地物波譜庫和多光譜衛(wèi)星數(shù)據(jù), 提出了一種利用非負(fù)矩陣分解法重建中紅外全波段地表反射率/發(fā)射率的方法, 并選取了4種典型地物類型(土壤、 植被、 人造材料和巖石)進(jìn)行方法驗證。 針對短波紅外和中紅外波段范圍內(nèi)的典型地物光譜重建方法進(jìn)行研究, 可為后續(xù)的全波段高光譜和多光譜遙感應(yīng)用提供關(guān)鍵支撐。

目前對光譜數(shù)據(jù)可視化的研究還處于傳統(tǒng)制圖階段, 往往只能進(jìn)行二維平面展示, 已無法滿足多源數(shù)據(jù)三維可視化表達(dá)的需求。 由于網(wǎng)絡(luò)的發(fā)展與信息量的激增, 網(wǎng)絡(luò)地理信息系統(tǒng)(web geographic information system, WebGIS)逐漸成為GIS的一個主流發(fā)展方向, 并且隨著對三維可視化的需求日益增長, 三維WebGIS應(yīng)運而生。 針對地表反射率重建結(jié)果可視化進(jìn)行研究, 基于WebGIS技術(shù)開發(fā)了一套以衛(wèi)星產(chǎn)品為中心, 將衛(wèi)星底圖、 地形數(shù)據(jù)與光譜重建結(jié)果等集成展示的二三維一體化可視化系統(tǒng), 使空間關(guān)系能得到更明確的展示, 表現(xiàn)力遠(yuǎn)超二維地理信息系統(tǒng), 擴(kuò)展了三維信息的應(yīng)用深度, 為衛(wèi)星產(chǎn)品的展示與研究提供新方案。

1 實驗部分

1.1 數(shù)據(jù)介紹

采用的光譜數(shù)據(jù)為約翰霍普金斯大學(xué)(JHU)光譜數(shù)據(jù)集, 該光譜庫提供了15個子庫, 包括巖石、 礦物、 月球土壤、 陸地土壤、 人造材料、 隕石、 植被、 冰雪等, 其中大部分測量結(jié)果來自JHU或美國地質(zhì)調(diào)查局[6]。

在中紅外范圍內(nèi), 地表特性的衛(wèi)星數(shù)據(jù)較少, 而中分辨率成像光譜儀(moderate-resolution imaging spectroradiometer, MODIS)能提供較全波段范圍的地表特性數(shù)據(jù), 因此選取MODIS數(shù)據(jù)進(jìn)行研究。 作為搭載于EOS(Earth Observing System)系列衛(wèi)星上的關(guān)鍵傳感器, MODIS因其較高的空間分辨率、 光譜分辨率、 時間分辨率優(yōu)勢, 在環(huán)境監(jiān)測、 土地利用、 作物估產(chǎn)等各個領(lǐng)域中得到廣泛應(yīng)用[7-9]。 本工作所采用的產(chǎn)品是MODIS的地表反射率產(chǎn)品MYD09A1和發(fā)射率產(chǎn)品MOD11C3。 其中, MYD09A1提供了根據(jù)大氣條件(例如氣體、 氣溶膠和瑞利散射)校正的AquaMODIS波段1～7的地表光譜反射率估計值; MOD11C3在0.05緯度/經(jīng)度氣候模擬網(wǎng)格(CMG)中提供了月度地表溫度和發(fā)射率(LST&E)值, 且LST&E值是通過對MOD11C1每日文件相應(yīng)月份的值進(jìn)行合成和平均得出的。

在2.0～5.0 μm范圍的短波紅外和中紅外波段, 由于MODIS只有2.105～2.155 μm的地表反射率產(chǎn)品, 以及3.66～3.84, 3.93～3.99和4.02～4.08 μm的地表發(fā)射率產(chǎn)品, 所以選取對應(yīng)的第7、 20、 22、 23波段作為模擬2.0～5.0 μm目標(biāo)通道的數(shù)據(jù)源通道(如表1所示), 其中, 中心波長2.13 μm為地表反射率產(chǎn)品的通道, 中心波長3.75、 3.96、 4.05 μm為地表發(fā)射率產(chǎn)品的通道。

表1 MODIS反射率/發(fā)射率數(shù)據(jù)獲取波段信息

1.2 非負(fù)矩陣分解方法

矩陣分解是一種有效降維方法, 可實現(xiàn)大規(guī)模數(shù)據(jù)處理與分析, 通過迭代分解的方式, 近似地將原始矩陣分解為2個較低維數(shù)矩陣的乘積, 進(jìn)而提取高光譜數(shù)據(jù)特征[10]。 非負(fù)矩陣分解(NMF)通過對基向量與對應(yīng)的系數(shù)增加非負(fù)約束, 以保證分解結(jié)果均為正值, 在高光譜圖像混合像元分解上取得應(yīng)用[11-13]; 基矩陣?yán)孟禂?shù)矩陣作為權(quán)重, 重構(gòu)原始矩陣。 對于一個n×d非負(fù)矩陣An×m(如:n表示光譜數(shù)據(jù)集的條數(shù),d表示每條光譜的波段個數(shù)), 存在一個非負(fù)矩陣Vn×k和Hn×k, 有

An×d≈Vn×kHk×d

(1)

式(1)中,V為基矩陣;H為系數(shù)矩陣;k為分解得到矩陣的秩, 且k?min(n,d);A中的列向量可解釋為對基矩陣V中所有列向量的線性組合, 且權(quán)重系數(shù)為系數(shù)矩陣V中對應(yīng)列向量中的元素。

采用歐式(Euclidean)距離作為目標(biāo)函數(shù)衡量矩陣分解前后的逼近程度, 非負(fù)矩陣分解問題被表示成如式(2)最優(yōu)化問題

s.t.V≥0,H≥0

(2)

式(2)中, ‖‖F(xiàn)表示矩陣Frobenius范數(shù); s.t.為subjectto的縮寫, 對應(yīng)約束條件。

表2列出了4種典型的NMF算法和非負(fù)性約束, 經(jīng)過多次實驗, 最終采用了較為穩(wěn)定的BPAS-NMF方法。

表2 部分可用的非負(fù)矩陣分解Matlab軟件包及約束條件和下載地址

1.3 基于挑選波段的病態(tài)光譜重建方法

與基于PCA的光譜重建類似, 地表反射率的光譜向量可以表示為

r≈Vh

(3)

寫成具體的向量和矩陣形式為

(4)

式(4)中, 矩陣V的列向量即為提取的端元向量, 可通過對光譜數(shù)據(jù)集的非負(fù)矩陣分解處理得到;h為對應(yīng)的豐度權(quán)重系數(shù)。

對于2.0～5.0 μm光譜范圍的地表反射率/發(fā)射率數(shù)據(jù)集, 基于已提取的端元向量, 若利用衛(wèi)星遙感常用的幾個典型波段的地表反射率/發(fā)射率信息來重建連續(xù)波段的地表反射率/發(fā)射率, 就必須先獲取對應(yīng)的豐度權(quán)重系數(shù)向量h的值。 利用ρ來表示從r中挑選的m個波段地表反射率組成的向量(m?d), 已知的幾個波段對應(yīng)的下標(biāo)依次用d1,d2, …,dm來表示,Vρ表示V這幾個波段對應(yīng)的端元子矩陣, 目的就是利用ρ來確定h的近似值, 進(jìn)而進(jìn)行光譜重建。 這樣,ρ可以表示為

r≈Vρh

(5)

和

(6)

基于式(6), 豐度系數(shù)向量h可以利用一個最小二乘的形式進(jìn)行求解

(7)

1.4 二三維一體化可視化方法

WebGIS是網(wǎng)絡(luò)與地理信息系統(tǒng)的結(jié)合, 作為GIS發(fā)展的一個主流模式, 已廣泛應(yīng)用于災(zāi)害響應(yīng)、 室內(nèi)導(dǎo)航、 資源管理、 環(huán)境監(jiān)測等各個領(lǐng)域, 并逐漸與三維可視化技術(shù)相結(jié)合。

WebGL是HTML5標(biāo)準(zhǔn)的一部分, 是一種新的三維Web繪圖標(biāo)準(zhǔn), 遵從HTML5標(biāo)準(zhǔn)的瀏覽器自帶對WebGL的支持, 降低了用戶的使用成本, 目前火狐、 谷歌、 歐朋等瀏覽器均實現(xiàn)了對WebGL的支持, 用戶只需安裝瀏覽器, 即可在瀏覽器頁面中實現(xiàn)強大的三維功能, 為用戶使用三維技術(shù)創(chuàng)造了極大的便利。 WebGL應(yīng)用在包含HTML、 CSS與JavaScript文件的傳統(tǒng)Web應(yīng)用基礎(chǔ)上, 還包括了三維模型數(shù)據(jù)及著色器的源代碼, 用JavaScript調(diào)用WebGLAPI即可快速實現(xiàn)Web三維功能。

主流的三維WebGIS產(chǎn)品包括Cesium、 ArcGIS API for Javascript和Super Map iClient 3D for WebGL等, 考慮到開發(fā)成本和需求, 采用完全開源的、 基于WebGL的Cesium框架。 Cesium是一個成本低、 開發(fā)簡單、 支持多種數(shù)據(jù)可視化方法、 可實現(xiàn)二三維一體化, 面向三維地球的JavaScript庫, 可在不同平臺、 不同瀏覽器展示三維地球, 不需任何插件支持。 Cesium提供基于JavaScript語言的開發(fā)包, 可實現(xiàn)全球高精度的地形和影像服務(wù), 且支持多種場景模式(2D, 3.5D以及3D場景), 使用者可迅速創(chuàng)建零插件的虛擬地球網(wǎng)絡(luò)應(yīng)用, 實現(xiàn)二三維一體化[14]。

本研究采用的B/S架構(gòu), 本質(zhì)是將瀏覽器作為客戶端, 于服務(wù)器端進(jìn)行交互, 具體運作模式如圖1所示。

圖1 B/S架構(gòu)

瀏覽器端通過HTM5和CSS實現(xiàn)界面及交互設(shè)計, 包括圖層切換、 拖拽、 視角切換等, 并進(jìn)行場景渲染。 服務(wù)器端使用JS語言, 利用Cesium框架提供的接口, 調(diào)用已發(fā)布的瓦片地圖, 并進(jìn)行數(shù)據(jù)讀取、 分析及三維可視化。

2 結(jié)果與討論

2.1 基于JHU光譜數(shù)據(jù)集的NMF中紅外光譜重建

為論證本文所提出的光譜重建方法的可行性, 首先將該方法應(yīng)用于地物波譜庫數(shù)據(jù)集的重建上, 基于地物波譜庫中典型地物的波段關(guān)系, 建立2.5～5.0 μm中紅外波段的光譜重建模型, 并對重建結(jié)果進(jìn)行了驗證。

(1) JHU光譜數(shù)據(jù)預(yù)處理

MODIS傳感器數(shù)據(jù)源通道的光譜響應(yīng)函數(shù)如圖2所示。

圖2 目標(biāo)函數(shù)光譜響應(yīng)函數(shù)(band7中心波長2.13 μm、 band20中心波長3.75 μm、 band22中心波長3.96 μm, band23中心波長4.05 μm)

由于選取的各類地物樣本光譜分辨率不盡相同, 因此需要對地物光譜進(jìn)行光譜間隔插值, 采用線性插值將地物光譜間隔統(tǒng)一為10 nm, 進(jìn)而基于光譜響應(yīng)函數(shù)對地物光譜庫中地物光譜數(shù)據(jù)進(jìn)行重采樣。 光譜響應(yīng)函數(shù)是波長的函數(shù), 表示傳感器在每個波長處所接收的輻射能量與入射的輻射能量的比值。 光譜響應(yīng)函數(shù)表征的是不同波長處信號響應(yīng)的差異, 而傳感器光譜響應(yīng)函數(shù)在真實連續(xù)光譜上進(jìn)行加權(quán)平均采樣后得到光學(xué)遙感圖像。 在構(gòu)建轉(zhuǎn)換模型時, 會用到發(fā)射率數(shù)據(jù)。 而地物光譜庫中的數(shù)據(jù)均為地物反射率數(shù)據(jù), 因此可以通過基爾霍夫定律對發(fā)射率數(shù)據(jù)進(jìn)行轉(zhuǎn)換, 從而得到反射率數(shù)據(jù)。 基爾霍夫定律如式(8)所示

ρEmis+ρRefl=1

(8)

式(8)中,ρEmis為發(fā)射率,ρRefl為反射率。 因此, 利用該公式可實現(xiàn)反射率數(shù)據(jù)和發(fā)射率數(shù)據(jù)的轉(zhuǎn)換。

基于已建立好的地物光譜庫樣本信息, 利用MODIS傳感器的光譜響應(yīng)函數(shù), 對數(shù)據(jù)源通道的地物反射率數(shù)據(jù), 根據(jù)等效計算公式進(jìn)行重采樣, 等效計算公式如式(9)所示

(9)

式(9)中,λ為波長,ε(λi)和f(λi)分別是波長為λi時的輻射數(shù)據(jù)和響應(yīng)函數(shù),n為通道j范圍內(nèi)的輻射數(shù)據(jù)或響應(yīng)函數(shù)值的個數(shù)。

(2) JHU地物反射率光譜數(shù)據(jù)集端元光譜向量的提取

利用線性插值將JHU波譜數(shù)據(jù)庫2.0～5.0 μm波段范圍反射率結(jié)果重采樣到10 nm光譜間隔的301個波段, 具體如圖3所示。

圖3 JHU波譜數(shù)據(jù)庫提取的2～5 μm波段范圍10 nm間隔的四種典型地表類型光譜數(shù)據(jù)集

在此基礎(chǔ)上, 進(jìn)行非負(fù)矩陣分解處理, 提取得到對應(yīng)的端元向量, 如圖4所示, 這4條端元向量即為用于非負(fù)矩陣光譜重建可通用的基準(zhǔn)光譜向量。

圖4 提取的2～5 μm波段范圍10 nm間隔的用于NMF光譜重建的4個端元向量

(3) JHU 地物反射率光譜重建和結(jié)果驗證

在獲取4條端元向量光譜曲線的基礎(chǔ)上, 從JHU地表反射率光譜數(shù)據(jù)集中提取MODIS衛(wèi)星中心波長為2.13、 3.75、 3.96和4.05 μm這4個波段的地表反射率的子數(shù)據(jù)集, 并計算對應(yīng)的權(quán)重系數(shù)向量結(jié)果, 進(jìn)行2.0～5.0 μm光譜范圍內(nèi)的全波段的反射率光譜重建。 利用MODIS的4個波段反射率結(jié)果進(jìn)行全波段光譜重建的結(jié)果如圖5所示, 對應(yīng)的散點圖如圖6所示, 對應(yīng)的平均絕對誤差為0.01, 平均相對誤差為10%。 對應(yīng)在只有MODIS衛(wèi)星4個波段數(shù)據(jù)可用的低秩病態(tài)的情況下, 可較好地滿足光譜重建的精度要求。

2.2 基于衛(wèi)星觀測數(shù)據(jù)的全球中紅外任意波段光譜重建

在對光譜重建方法進(jìn)行驗證的基礎(chǔ)上, 將其應(yīng)用于廣泛使用、 易于獲取、 波段較全的MODIS反射率/發(fā)射率產(chǎn)品的光譜重建。 利用MOEIS的地表反射率和發(fā)射率產(chǎn)品, 推導(dǎo)計算到任意中紅外波段的地表反射率。

(1) MODIS衛(wèi)星觀測數(shù)據(jù)預(yù)處理

MODIS獲取的地表反射率數(shù)據(jù)和發(fā)射率數(shù)據(jù)存在數(shù)據(jù)無效或缺失的問題, 針對這一問題, 一方面基于往年同期、 同年前后時間的反射率補丁數(shù)據(jù)對原始反射率數(shù)據(jù)中的無效像元進(jìn)行修復(fù); 另一方面利用原始反射率數(shù)據(jù)中獲取到的海表反射率數(shù)據(jù), 計算其均值和標(biāo)準(zhǔn)差, 并基于正態(tài)分布模型對海洋范圍缺失數(shù)據(jù)進(jìn)行修復(fù)。 圖7展示了經(jīng)過預(yù)處理得到的MODIS全球月均地表反射率/發(fā)射率產(chǎn)品, 其中圖7(a)為中心波長2.13 μm的空間分辨率500 m的地表反射率產(chǎn)品的通道, 圖7(b—c)分別對應(yīng)中心波長3.75、 3.96、 4.05 μm的空間分辨率5 km的地表發(fā)射率產(chǎn)品。

圖7 MODIS全球月均地表反射率/發(fā)射率產(chǎn)品

(2) 全球中紅外波段光譜重建

在處理好MODIS短波紅外和中紅外4個波段的全球月均地表反射率/發(fā)射率產(chǎn)品的基礎(chǔ)上, 利用基于NMF和光譜數(shù)據(jù)集所提取的2.0～5.0 μm用于反射率光譜重建的端元向量, 可以計算得到每個像元對應(yīng)的權(quán)重系數(shù)向量, 進(jìn)而可進(jìn)行任意波段的光譜重建, 得到的全球范圍內(nèi)陸地5 km×5 km分辨率的月均地表反射率如圖8所示, 通過全球陸地地表反射率的光譜重建, 可為后續(xù)中紅外相關(guān)波段的背景輻射研究提供重要支撐。

圖8 全球范圍內(nèi)陸地5 km×5 km分辨率的月均地表反射率

2.3 基于WebGIS的二三維一體化可視化

基于Cesium框架, 采用B/S架構(gòu), 搭建了二三維一體化可視化系統(tǒng), 將衛(wèi)星底圖、 地形數(shù)據(jù)與光譜重建結(jié)果等集成展示。 實現(xiàn)了圖層切換、 拖拽、 視角切換等地圖控制功能及瓦片地圖調(diào)用, 數(shù)據(jù)讀取、 分析及三維可視化等數(shù)據(jù)處理功能。

將全球中紅外光譜重建結(jié)果進(jìn)行投影轉(zhuǎn)換, 建立符號系統(tǒng)并進(jìn)行影像切片處理發(fā)布服務(wù), 即可在系統(tǒng)中調(diào)用, 將地形服務(wù)進(jìn)行夸張?zhí)幚砗笈c之疊加, 可直觀展現(xiàn)出地表反射率與地形之間的關(guān)系如圖9所示。 圖9(a)展示了圖層疊加后地表反射率重建結(jié)果在三維地球上的展示效果, 圖9(b)近距離展示了地表反射率圖層覆蓋于地形上的效果, 可以直觀看出, 該區(qū)域內(nèi)高程與地表反射率總體呈負(fù)相關(guān)。

圖9 圖層疊加展示

3 結(jié) 論

提出了一種基于地物波譜庫和多光譜衛(wèi)星數(shù)據(jù), 利用非負(fù)矩陣分解方法重建高光譜地表反射率的方法, 實現(xiàn)了2.5～5.0 μm光譜范圍內(nèi)10 nm間隔的地表反射率光譜重建, 同時為綜合評價該光譜重建方法, 選取了4種典型地物類型(土壤、 植被、 人造材料和巖石)進(jìn)行方法驗證。 研究結(jié)果表明:

(1)本文提出的NMF光譜重建方法, 相較于PCA等方法, 增加了非負(fù)約束, 具有可解釋的物理意義, 得到的重建結(jié)果較為穩(wěn)定和精確。 利用本方法得到的典型地表類型反射率光譜重建結(jié)果與真實結(jié)果的平均絕對誤差為0.01, 平均相對誤差為10%, 誤差較小, 有效論證了該方法的可行性, 對應(yīng)在只有MODIS衛(wèi)星4個波段數(shù)據(jù)可用的低秩病態(tài)的情況下, 可較好地滿足光譜重建的精度要求。

(2)基于本文提出的方法, 利用MODIS的地表反射率和發(fā)射率產(chǎn)品, 可獲得全球范圍內(nèi)地表反射率或發(fā)射率高光譜數(shù)據(jù), 能夠在一定程度上彌補衛(wèi)星手段難以獲取紅外高光譜地表反射率數(shù)據(jù)的問題, 并且該方法可進(jìn)一步擴(kuò)展到可見光光譜范圍內(nèi)的地表反射率重建, 為衛(wèi)星遙感的相應(yīng)應(yīng)用提供支撐。

(3)為將地表反射率與對應(yīng)地表覆蓋和地形更直觀地進(jìn)行可視化分析, 基于Cesium開發(fā)庫進(jìn)行設(shè)計, 采用B/S架構(gòu), 搭建了二三維一體化可視化系統(tǒng), 將衛(wèi)星底圖、 地形數(shù)據(jù)與光譜重建結(jié)果等集成展示, 能直觀進(jìn)行多因素分析, 擴(kuò)展了三維信息應(yīng)用的深度, 為衛(wèi)星產(chǎn)品的展示與驗證提供新方案。