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

徐 君，展愛云，劉志偉

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

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

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

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

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

1 基于OSP的端元估計算法

線性光譜混合模型的數學公式：

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

OSP算法的思想就是將式（1）轉換為式（2）的形式：

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3 實驗

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

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

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

圖1 Cuprite地區的AVIRIS數據假彩色合成圖Fig.1 False-color image ofAVIRIS data in Cuprite

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

從表1可以看出，用本文方法解混后所得的端元與USGS光譜庫中的標準光譜相似度較高，總體上優于頂點成分分析（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方法估計圖像中端元的個數，從而更合理的對NMF算法進行初始化；針對MVC-NMF算法中單形體體積計算公式較為復雜而且需要降維的情況，提出了一種簡單的單形體體積的計算方法作為約束條件構建NMF的目標函數。實驗結果表明該方法解混效果跟實地調查報告分布情況基本吻合。

[1]賈森.非監督的高光譜圖像解混技術研究[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]鄭軼，蔡體健.稀疏表示的人臉識別及其優化算法[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的端元個數估計方法[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.