基于OSP與NMF的光譜混合像元分解方法

徐 君，展愛云，劉志偉

（華東交通大學信息工程學院，江西南昌330013）

高光譜遙感技術可以同時提供空間域信息和光譜域信息，得到地物近似平滑的光譜曲線，因此近年來廣泛地應用在土地資源調(diào)查、環(huán)境監(jiān)測、軍事偵察等領域[1-3]。但大多數(shù)高光譜成像儀的空間分辨率并不高，這就導致在同一個像元中可能包含多種地物，所提取的光譜特征是多種地物光譜特征的混合，這種像元被稱為光譜混合像元。混合像元的存在使得傳統(tǒng)的基于像元級光譜統(tǒng)計特性的分類方法無法直接使用，這是因為直接將這類像元歸屬到哪一種典型地物都是不準確的，因此混合像元問題是遙感技術向定量化發(fā)展的重要障礙?；旌舷裨忻糠N地物的比例成為地物分布的豐度，反演混合像元豐度的過程稱為混合像元的解混。

一般光譜混合像元分為線性混合模型和非線性混合模型兩類[4]。非線性光譜混合模型比較復雜、物理意義并不明確。線性光譜混合模型建立與解釋相對容易，物理意義明確，精度能夠滿足大部分應用的要求。因此在一般情況下，可以用線性模型來描述混合像元的形成機理，大部分混合像元分解算法都是基于線性模型的。

非負矩陣分解[5-8]（nonnegative matrix factorization，NMF）是20世紀末出現(xiàn)的一種盲源分解算法。由于NMF與線性光譜混合模型的公式之間有著很高的相似性，這就為NMF在混合像元分解中的應用提供了可能。將NMF應用于高光譜混合像元的分解不需要假定純像元的存在，并且在提取端元的同時可以獲取每種端元對應的豐度圖。然而NMF的目標函數(shù)具有明顯的非凸性，因而存在大量局部極小。對于NMF問題V≈WH，（V，W，H均為非負矩陣），可以找到大量的非負可逆矩陣D及其逆矩陣D-1，有WH=(WD)(D-1H)成立，這樣可以得到很多對解WD和D-1H，這是NMF存在的最大的問題。

本文首先通過正交子空間投影（orthogonal subspace projection，OSP）技術估計端元的個數(shù)[9]，然后通過改進的外包單形體最小體積約束非負矩陣分解算法（minimum volume constrained nonnegative matrix factor?ization，MVC-NMF）對高光譜混合像元進行分解[10]，同時得到目標景物的端元光譜及豐度分布圖。

1 基于OSP的端元估計算法

線性光譜混合模型的數(shù)學公式：

式中：X是高光譜圖像中某個光譜混合像元的光譜矢量；A為m×n的端元矩陣；m為波段數(shù)目；n為端元數(shù)目；S為豐度矩陣；ε為噪聲或模型的誤差。

OSP算法的思想就是將式（1）轉(zhuǎn)換為式（2）的形式：

式（2）中的端元信號可以表示為目標信號D和背景信號U的組合，fD和fU分別為對應的豐度。高光譜圖像可視為由一組端元所張成的子空間中的高維數(shù)據(jù)點的集合，而OSP算法能實現(xiàn)目標和背景分離的關鍵就在于如何利用端矩陣A構造一個投影矩陣，從而能將光譜向量通過投影到空間的正交子空間中去。

本文所提算法需隨迭代次數(shù)的增加逐步提取端元，并將其逐個加入到端元矩陣A中，相應的也隨之增大。端元提取算法采用單形體生長算法（simplex growing algorithm，SGA），該算法能穩(wěn)定的逐個提取端元，保證提取結果的一致性。采用高光譜圖像中所有像元光譜矢量的平均值μ作為待投影的矢量來進行迭代，具體步驟如下：

1）對實驗所用的高光譜數(shù)據(jù)進行去噪及幾何校正等預處理；

3）設定迭代的終止條件閾值σ；

4）開始迭代，在每次迭代過程中都為端元矩陣Ai增加一個端元光譜ai，得到

2 改進的MVC-NMF混合像元分解算法

對于線性光譜混合模型而言，豐度矩陣S存在兩個限制條件：非負性約束（abundance nonnegative con?straint，ANC）與和為一約束（abundance sum-to-one constraint，ASC）。如果不考慮其他條件，僅在這兩個約束條件下NMF算法結果不具有唯一性。因此，在將NMF應用于光譜數(shù)據(jù)的線性解混時一般需要引入其他的約束條件。

根據(jù)凸面幾何學的理論，高光譜數(shù)據(jù)位于由端元光譜作為頂點的單形體內(nèi)，像元在單形體中的位置決定了混合像元各端元的豐度，而最優(yōu)的單形體應該外包特征空間中的所有點且具有最小的體積，因此在NMF光譜解混算法過程中可以加入最小單形體體積約束來優(yōu)化算法，實現(xiàn)混合像元的精確分解。

文獻10中的算法正是通過求最大單形體的體積作為約束求解NMF算法而得到各個端元的，其體積公式如下：

式中：an是第n個端元的列向量；V(A)是由這n個端元作為頂點構成的單形體的體積。由于用到了求行列式的運算，所以要求A必須為方陣，這樣向量an的維數(shù)必須是n-1，這就需要先對原始數(shù)據(jù)進行降維處理，因此計算可能導致偏差，這是此算法的不足之處。

由于m維單形體包含于m維平行多面體中，因此它們的體積對應關系為為m維單形體體積，為m維平行多面體體積。

假設a1，a2，...，an是光譜圖像圖像中的n個像元的光譜矢量，如令則在特征空間中，以這n個點為頂點的平行多面體體積為，因此以這n個點為頂點的n-1維單形體體積J(Vn)可通過下式計算：

式（5）表示的單形體體積計算公式對于任何維數(shù)的高維空間都是成立的。據(jù)此可以構建改進后的最小體積約束MVC-NMF目標函數(shù)：

式（6）中λ=00.5，取上述算法的初始化主要包括端元數(shù)目估計和對矩陣A進行初始化兩部分。本文采用OSP算法估計端元數(shù)目p，端元矩陣A可以采用隨機從高光譜圖像數(shù)據(jù)中選擇p個像元光譜進行初始化。

求解過程采用投影梯度算法以加速式（6）的迭代，具體方法是將式（6）分解為兩個子問題來求解，每個子問題固定A（或S）對S（或A）進行更新：

迭代過程可以采用最大迭代次數(shù)、誤差閾值、用戶交互等作為終止條件，本文采用相鄰兩次迭代目標函數(shù)值的比值小于某個閾值0.99的方法來實現(xiàn)算法的終止。

3 實驗

實驗采用Cuprite地區(qū)的AVIRIS數(shù)據(jù)進行光譜解混實驗（如圖2），AVIRIS數(shù)據(jù)包含224個波段（波長范圍為400～2 500 nm），是一種高質(zhì)量、低噪聲的高光譜數(shù)據(jù)。這里我們選取其中50個連續(xù)的波段（1 978～2 478 nm）的圖像進行算法的驗證。

首先利用OSP方法估算出端元的個數(shù)，實驗中設置迭代終止條件C＜0.01，最終估算出有9個端元。

利用本文提出的改進的最小體積約束MVC-NMF算法對整幅圖像進行光譜解混，圖2顯示了解混后得到的地表礦物的豐富分布圖。與Swayze和Clark等人已經(jīng)給出的該地區(qū)真實地物的分布報告[11]相比較，可確定這些端元各自對應的礦物。這說明本文所提出的算法基本反映了此區(qū)域的地物分布狀況。

圖1 Cuprite地區(qū)的AVIRIS數(shù)據(jù)假彩色合成圖Fig.1 False-color image ofAVIRIS data in Cuprite

為了進一步比較，我們將美國地質(zhì)調(diào)查局（United States Geological Survey，USGS）庫中的對應礦物光譜作為參考光譜，將參考光譜在相應的波段范圍內(nèi)能量求積分，轉(zhuǎn)換為光譜向量，并求取解混所得的端元與它們之間的光譜角，如表1所示。

從表1可以看出，用本文方法解混后所得的端元與USGS光譜庫中的標準光譜相似度較高，總體上優(yōu)于頂點成分分析（vertex component analysis，VCA）方法所提取的端元光譜。

表1 所提取的端元與USGS光譜庫中標準光譜之間的光譜角Tab.1 The SAD between endmembers extracted and the standard spectrum

圖2 不同礦物的豐度分布圖Fig.2 Abundance maps of different minerals

4 結論

本文提出了一種基于OSP和NMF方法的高光譜混合像元分解方法。首先用OSP方法估計圖像中端元的個數(shù)，從而更合理的對NMF算法進行初始化；針對MVC-NMF算法中單形體體積計算公式較為復雜而且需要降維的情況，提出了一種簡單的單形體體積的計算方法作為約束條件構建NMF的目標函數(shù)。實驗結果表明該方法解混效果跟實地調(diào)查報告分布情況基本吻合。

[1]賈森.非監(jiān)督的高光譜圖像解混技術研究[D].杭州：浙江大學，2007.

[2]張良培，張立福.高光譜遙感[M].武漢：武漢大學出版社，2005：15-18.

[3]錢樂祥，泮學芹，趙芊.中國高光譜成像遙感應用研究進展[J].國土資源遙感，2004，7（2）：1-6.

[4]KESHAVAN，MUSTRAD J F.Spectral unmixing[J].IEEE Signal Processing Magazine，2002，19（1）：44-57.

[5]JIA S，QIAN Y.Constrained nonnegative matrix factorization for hyperspectral unmixing[J].IEEE Trans.Geocsi.Remote sens.，2009，47（1）：161-173.

[6]LIU X，XIA W，WANG B，ZHANG L.An approach based on constrained nonnegative matrix factorization to unmix hyperspectral data J].IEEE Trans.Geocsi.Remote sens.，2011，49（2）：757-772.

[7]鄭軼，蔡體健.稀疏表示的人臉識別及其優(yōu)化算法[J].華東交通大學學報，2012，29（1）：11-14.

[8]PAUCA V，PIPER J，PLEMMONS R.Nonnegative matrix factorization for spectral data analysis[J].Linear Algebra and Applications，2006，416（1）：29-47.

[9]陳偉，余旭初，王鶴.基于OSP的端元個數(shù)估計方法[J].計算機工程，2009，35（24）：280-281.

[10]MIAO L，QI H.Endmember extraction from highly mixed data using minimum volume constrained nonnegative matrix factorization[J].IEEE Trans Geosci Remote Sens，2007，45（3）：765-777.

[11]SWAYZE G A.The hydrothermal and structural history of the cuprite mining district，southwestern Nevada：an integrated geological and geophysical approach[D].Boulder：University of Colorado，1997：399.