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

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

1 概述

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

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

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

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

2 線性光譜混合模型

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

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

3 非平滑約束的非負(fù)矩陣分解

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

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

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

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

應(yīng)用 nsNMF去分解混合像元，為達(dá)到全局稀疏性[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ù)。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

5 實驗與結(jié)果分析

5.1 合成數(shù)據(jù)實驗

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

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

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

表1 3個算法的準(zhǔn)確性比較

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

5.2 真實數(shù)據(jù)實驗

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

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

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

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

圖2 白云母的豐度結(jié)果

圖3 榍石的豐度結(jié)果

圖4 明礬石的豐度結(jié)果

6 結(jié)束語

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

[1]賈 森.非監(jiān)督的高光譜圖像解混技術(shù)研究[D].杭州:浙江大學(xué), 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]基于非負(fù)矩陣分解的高光譜遙感圖像混合像元分解[J].紅外與毫米波學(xué)報, 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.