基于粒子群優(yōu)化的非平滑非負矩陣分解算法

戴華平，王 旭，胡紅亮，王玉濤

(浙江大學工業(yè)控制研究所，杭州 310027)

1 概述

隨著成像光譜技術的不斷發(fā)展，通過成像光譜儀采集得到的遙感圖像包含了越來越豐富的空間、輻射和光譜信息。影響光譜技術進一步發(fā)展的關鍵問題是高光譜圖像中混合像元的廣泛存在，混合像元限制了遙感應用的精度，降低了物質分類的準確性。因此，有效分解混合像元對于高光譜技術的發(fā)展具有重要意義[1]。

非負矩陣分解(Nonnegative Matrix Factorization,NMF)[2]是一種新的矩陣分解方法，目前已經被用于高光譜的混合像元分解。然而非負矩陣分解只是提供了非負性約束，這遠不能很好地解決混合像元的分解問題。對于高光譜數(shù)據來說，除了非負性約束之外，稀疏性也是一個重要的約束條件，即任何端元的豐度分布一般都不會充滿整個圖像空間。文獻[3]提出的非平滑非負矩陣分解(nonsmooth Nonnegative Matrix Factorization, nsNMF)算法被用于高光譜圖像分解并顯示出良好的效果，但是NMF和nsMNF具有的共同缺陷是容易陷入局部最小值，限制了高光譜解混的精度。

粒子群優(yōu)化(Particle Swarm Optimization, PSO)[4]算法又稱微粒群算法，是由Kennedy J和Eberhart R C等人于 1995年開發(fā)的一種演化計算技術，來源于對一個簡化社會模型的模擬。粒子群算法可以解決多數(shù)的優(yōu)化問題，并且它具有簡單、高效和魯棒性強的優(yōu)點。

本文提出了一種基于PSO的nsNMF新算法，將nsNMF得到的結果作為PSO的初始值，每次迭代過后重新采用nsNMF進行更新，直到找到全局最優(yōu)值。

2 線性光譜混合模型

線性光譜混合模型(Liner Spectral Mixture Model,LSMM)[5]在高光譜圖像分析中有著廣泛的應用。在線性光譜混合模型中，像素的觀察值等于各端元的光譜特征按照它們的豐度進行線性組合，數(shù)學模型表述如下：

其中，矩陣 R ∈RL×N表示L個波段的遙感圖像，每個波段有N個像素。 M =[m1, m2,… ,mp]∈ RL×P，稱為端元光譜矩陣，mj對應著第 j個端元的光譜特征，P表示端元個數(shù)。 S =[s1, s2,… ,sp]T∈RP×N為豐度矩陣， sj=[sj1,sj2,… ,sjN]T ∈RN×1表示某一個端元上的所有像素的豐度。根據豐度的物理意義，豐度必須滿足和為一約束和非負約束：

3 非平滑約束的非負矩陣分解

對于高光譜數(shù)據來說，除了非負性約束和全加性約束之外，稀疏性也是一個重要的約束條件。一般情況下，任何端元的分布都不會充滿整個高光譜圖像的空間。混合像元中包含的只是幾種端元，而不是所有的端元。每個端元只是分布在高光譜圖像空間中的某些區(qū)域，具有一定的稀疏度。對于線性混合光譜模型(LSMM)來說，一個矩陣的稀疏性必然導致另一個矩陣的非稀疏性(平滑的)，使得最終的乘積可以更好的擬合觀測信號，這一點恰好加強了光譜信號的平滑性特征[3]。

給定一個L×N非負矩陣R，其中Rij≥0，非負矩陣分解算法找到2個非負矩陣M ∈RL×P和S ∈RP×N，所以得到：

為尋找一個近似的分解過程，使得R≈MS，通常定義R和MS的歐氏距離作為目標函數(shù)來保證逼近效果，即：

當且僅當R=MS時，上式取到最小值0。其中，Mij≥0,Sij≥0,? i, j 。

應用 nsNMF去分解混合像元，為達到全局稀疏性[3]，Montano 等人引入一個“平滑”矩陣C∈RP×P[4]到混合模型式(4)中：

其中，C是一個正對稱陣，其定義為：

其中，Ip是P×P維的單位矩陣；1p是元素全部為1的P×1維向量；稀疏因子θ滿足0 ≤ θ≤ 1。矩陣M的非稀疏性使得矩陣S為稀疏矩陣，同時矩陣S的非稀疏性也使得矩陣M為稀疏矩陣。在光譜解混問題中，只需要端元分布矩陣S稀疏的，因此，更新規(guī)則如下：

其中，α和β為規(guī)則化系數(shù)。

當且僅當M和S達到穩(wěn)定點時，目標函數(shù)式(5)不在變換。此外為了滿足端元分布的全加性約束，在每一步迭代之后還要對豐度矩陣進行歸一化運算：

4 非平滑非負矩陣分解的粒子群優(yōu)化

4.1 粒子群優(yōu)化算法

粒子群優(yōu)化算法首先初始化一群隨機粒子，然后通過迭代找到該種群的最優(yōu)解。在每一次迭代中，每個粒子通過跟蹤2個極值來更新自己。一個是粒子本身的最優(yōu)解，稱個體極值Pbest，另一個是整個種群的最優(yōu)解，稱全局極值Gbest，在本文中Gbest是一個矩陣且每個元素值相等。

種群中每個粒子在找到上述2個極值以后，根據以下公式來更新速度和位置[6]：

其中，V粒子速度；Pnow為粒子當前位置；Pbest和Gbest見前述定義；w加權系數(shù)，取值在0.1到0.9之間；rand()是(0, 1)之間的隨機數(shù)；c1和c2稱為學習因子，通常取2。

粒子通過不斷的學習和更新，最終找到解空間中最優(yōu)解所在的位置，整個搜索過程結束。最后輸出的Gbest就是全局最優(yōu)解。

粒子群優(yōu)化算法結構簡單、運行速度快，沒有交叉和變異運算。但粒子群在搜索時有時會在全局最優(yōu)解附近出現(xiàn)“振蕩”，為了避免“振蕩”的出現(xiàn)，做出如下改進[7]：隨著迭代的進行，速度更新公式中的加權系數(shù)w由wini線性減小到wend，即：

其中，iter是當前迭代的代數(shù)；itermax是總的迭代次數(shù)。

4.2 非平滑非負矩陣分解的粒子群優(yōu)化算法

基于 PSO 的約束的非負矩陣(Constrained Nonnegative Matrix Factorization, CNMF)算法[8]將頂點成分分析算法(Vertex Component Analysis, VCA)產生的結果作為該算法的初始值，但是VCA算法要求高光譜圖像中至少存在每種端元的一個純像元，實際中，這樣的要求往往得不到滿足，這樣就降低了光譜解混的精度。PSO-CNMF算法利用 PSO同時對端元光譜矩陣和豐度矩陣進行搜索，然后通過 CNMF來更新，解空間維數(shù)相對較大，增加了計算復雜度，并且容易陷入局部最優(yōu)值。真實高光譜圖像中，一般圖像的像元中不可能包括所有端元，只會包含其中的幾種端元，導致豐度矩陣的稀疏性，但 PSO-CNMF算法不能很好搜索稀疏性矩陣解空間。

本文提出一種基于 PSO 的 nsNMF算法(PSO-nsNMF)，該算法采用 nsNMF輸出的結果作為初始值，nsNMF算法沒有VCA算法中純像元的要求，更加接近真實情況。利用 nsNMF產生的初始值，可以有效避免 PSO的盲目搜索，提高找到全局最優(yōu)值的概率。PSO-nsNMF算法利用PSO搜索端元矩陣，通過 nsNMF更新端元特征矩陣和豐度矩陣，相對于PSO-CNMF算法有效的縮小了解空間，降低了計算復雜度，提高了獲得全局最優(yōu)值的概率。PSO- nsNMF算法中 nsNMF引入“平滑”矩陣 C，可保證了更新過程中豐度矩陣的稀疏性，與PSO-CNMF算法相比，更符合真實情況。

在 PSO-nsNMF算法中，首先利用傳統(tǒng) nsNMF算法迭代得到的結果去初始化 PSO-nsNMF算法的M和S。然后通過粒子群優(yōu)化產生的M和S讓nsNMF算法進行迭代，產生的結果M送給PSO進行優(yōu)化，重復整個過程直到找到全局最優(yōu)值。

在PSO-nsNMF中的PSO部分，本文選擇讓M粒子群中每個粒子的位置，通過與元素在 0.002～0.005之間的隨機矩陣相加增加粒子的多樣性。粒子群中每個粒子的速度隨機初始化。粒子群中每個粒子的最優(yōu)值記為 PMi_best，粒子群的全局最優(yōu)值記為 PMg_best，式(4)為適應值函數(shù)。

PSO-nsNMF算法的具體步驟如下：

(1)利用nsNMF算法輸出的結果去初始化算法中的M和S。

(2)利用算法中的nsNMF部分去更新M和S，并輸出結果。

(3)利用算法中的PSO部分去更新M，返回M的最優(yōu)值 PMg_best和S的最優(yōu)值 PSg_best。

(4)如果已經滿足停止條件，則整個迭代過程結束，返回 PMg_best和 PSg_best。否則跳轉到步驟(2)繼續(xù)迭代。

5 實驗與結果分析

5.1 合成數(shù)據實驗

實驗所用的模擬數(shù)據為 224波段 3端元模擬圖像，是從USGS數(shù)字光譜數(shù)據庫中選擇的3個純物質光譜，分別為光鹵石、銨明礬石、黑云母。每種光譜都為224× 1維數(shù)據，其中224為光譜波段數(shù)。然后利用 Dirichlet分布[9]產生要合成的混合像元的豐度矩陣，大小為3× 1 296，其中1 296為要合成的像元的總數(shù)，并作歸一化處理。為使得實驗更具實效性，在合成數(shù)據中添加SNR為20 dB的隨機噪聲。實驗分別比較標準 NMF、nsNMF和 PSO-nsNMF算法在合成數(shù)據中的結果。

對于標準NMF，M和S為隨機初始值，但滿足非負性和全加性約束。對于nsNMF算法，α和β分別為0和1，保證S的平滑性。θ取值為0.5保證一定的稀疏性。對于PSO-nsNMF、nsNMF算法產生的初始值作為粒子群的初始粒子。最大迭代次數(shù)設置為100，wini和wend分別設置為 0.9和0.4。學習因子c1和c2分別設置為 1.565和 1.565。

為評估 3種算法的有效性，利用光譜信息散度(Spectral Information Divergence, SID)和光譜角距離(Spectral Angle Distance SAD)來評估真實端元和估計端元之間的相似性，均方根誤差(Root Mean Square Error RMSE)[1]來評估真實豐度值和估計豐度值之間的相似性。最后的量化結果通過 10次實驗的平均值得到的，如表1所示。

表1 3個算法的準確性比較

從表 1可以看出，PSO-nsNMF所有指標均優(yōu)于標準NMF和nsNMF，因此，通過PSO-nsNMF分解的混合像元的端元值和豐度值都更接近真實值。

5.2 真實數(shù)據實驗

實驗所用數(shù)據是 AVIRIS在美國 Cuprite地區(qū)[9]采集得到的真實高光譜數(shù)據集，本文選用大小為250×190像素的子圖作為實驗對象，如圖1所示。在所有的224個波段中去除低噪聲比和水蒸氣吸收波段后(1～2, 104～113, 148～167和221～224)，剩余 188個波段[10]，即 L=188。

圖1 美國Cuprite地區(qū)高光譜數(shù)據集(波段80)

利用虛擬維度[11]的方法確定這塊圖像中包含的端元個數(shù)為10，即P=10。通過與光譜庫中的光譜對比得知這 10中端元分別為明礬石、榍石、白云母、綠脫石、鎂鋁石、膠嶺石、藍線、鐵鈣榴石、高嶺石、水銨長石。在實驗中選取白云母、榍石、明礬石3個典型的礦石作為例子來驗證算法的有效性。

由標準 NMF、nsNMF和 PSO-nsNMF分解得到端元豐度結果見圖2～圖4，顯示了3種礦石的端元分布圖，可以看出由PSO-nsNMF得到的端元分布圖比標準NMF和nsNMF得到的結果更接近真實值。

圖2 白云母的豐度結果

圖3 榍石的豐度結果

圖4 明礬石的豐度結果

6 結束語

本文提出一種基于粒子群優(yōu)化的非平滑非負矩陣分解的新算法 PSO-nsNMF，通過合成數(shù)據和真實的實驗結果得知，PSO-nsNMF可以獲得更好的全局最優(yōu)解，其端元光譜和豐度值更接近真實值。下一步將研究如何更好地選取系數(shù)w、c1、c2、wini和wend，并改進算法結構，加強算法的全局搜索能力。

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

[2]Lee D D, Seung H S.Algorithms for Non-negative Matrix Factorization[C]//Advances in Neural Information Processing Systems.Denver, USA: [s.n.], 2000: 556-562.

[3]Montano A P, Carazo J M, Kochi K, et al.Nonsmooth Nonnegative Matrix Factorization(nsNMF)[J].IEEE Transactions on Pattern Analysis and Machine Intelligence,2006, 28(3): 403-415.

[4]Kennedy J, Ebethart R C.Swarm Intelligence[M].San Francisco, USA: Morgan Kaufmann Publishers, 2001.

[5]Keshava N, Mustard J F.Spectral Un-mixing[J].IEEE Signal Processing Magazine, 2002, 19(3): 44-57.

[6]Shi Yuhui, Eberhart R.Parameter Selection in Particle Swarm Optimization[C]//Proceedings of EP’98.London,UK: Springer-Verlag, 1998: 1958-1962.

[7]Shi Yuhui, Eberhart R.A Modified Particle S-Warm Optimizer[C]//Proceedings of Conference on Evolutionary Computation.Piscataway, USA: IEEE Press, 1998: 69-73.

[8]Cui Jiantao, Li Xiaorun.Unsupervised Hyperspectral Unmixing Based on Constrainted Nonnegative Matrix Factorization and Partcle Swarm Optimization[C]//Proceedings of GCIS’10.[S.l.]: IEEE Computer Society,2010: 376-380.

[9]Nascimento J M P, Dias J M B.Vertex Component Analysis: A Fast Algorithm to Unmix Hyperspectral Data[J].IEEE Transactions on Geoscience and Romote Sensing, 2005, 43(4): 898-910.

[10]基于非負矩陣分解的高光譜遙感圖像混合像元分解[J].紅外與毫米波學報, 2011, 30(1): 27-32.

[11]Chang C I, Du Qian.Estimation of Number of Spectrally Distinct Signal Source in Hyperspectral Imagery[J].IEEE Transactions on Geoscience and Remote Sensing, 2004,42(3): 608-619.