基于豐度約束核非負(fù)矩陣分解的高光譜圖像非線性解混

智通祥，楊 斌，王 斌

(1.復(fù)旦大學(xué) 電磁波信息科學(xué)教育部重點(diǎn)實(shí)驗(yàn)室，上海 200433；2.復(fù)旦大學(xué) 信息科學(xué)與工程學(xué)院 智慧網(wǎng)絡(luò)與系統(tǒng)研究中心，上海 200433)

高光譜遙感圖像可以獲得目標(biāo)地區(qū)的空間分布和詳細(xì)的光譜信息，其光譜波段多達(dá)數(shù)百個(gè)，在環(huán)境監(jiān)測、地質(zhì)勘探、軍事偵察等方面有廣泛的應(yīng)用[1].但由于傳感器空間分辨率的限制，以及地物成分的復(fù)雜性，高光譜圖像中存在著大量的混合像元，嚴(yán)重影響了定量遙感的發(fā)展.為解決這一問題，就需要進(jìn)行光譜解混，從光譜圖像中提取出組成混合像元的純物質(zhì)光譜(即端元)和它們在混合像元中所占的比例(即豐度)[2].

在以往的高光譜解混研究中，常采用線性混合模型(Linear Mixing Model, LMM)來描述混合像元現(xiàn)象.這里，假設(shè)地物呈棋盤格網(wǎng)狀分布，并且光子只與場景中的一種物質(zhì)發(fā)生作用，這樣混合像元光譜就可以由端元光譜以不同的比例系數(shù)線性組合而成.由于LMM模型簡單，物理意義明確，近十年來，基于LMM的線性解混算法一直是國內(nèi)外學(xué)者研究的重點(diǎn)，并日趨成熟[3].但是在許多復(fù)雜的自然場景中存在著非線性混合效應(yīng)，如沙地、礦物地區(qū)的緊密混合現(xiàn)象，以及植被覆蓋、建筑林立地區(qū)的多層次混合現(xiàn)象.基于LMM的線性解混算法在這些非線性混合區(qū)域的解混結(jié)果會與實(shí)際情況存在較大誤差.因此，為了更準(zhǔn)確地進(jìn)行高光譜解混，非線性解混算法逐漸受到了學(xué)者的關(guān)注，多種非線性解混算法被提出，并用于非線性混合場景下的高光譜解混[4-5].

其中，針對沙地、礦物地區(qū)存在的緊密混合現(xiàn)象，Hapke模型[6]從輻射傳輸理論出發(fā)，推導(dǎo)出雙向反射率可以表示為與場景相關(guān)的粒子密度和大小以及單次散射反照率等物理參數(shù)的非線性函數(shù)，從而描述緊密混合情況下光線在微觀尺度上與不同粒子之間的多次散射情況.而植被覆蓋、建筑林立的地區(qū)則存在多層次混合現(xiàn)象.雙線性混合模型(Bilinear Mixing Model, BMM)在LMM的基礎(chǔ)上，考慮兩兩端元間存在的二次散射效應(yīng)并忽略更高次散射的影響，可以較為簡化地描述這種多層次混合現(xiàn)象.BMM主要由線性混合和非線性混合兩部分構(gòu)成，根據(jù)后者的不同可分為Nasimento模型(Nasimento Model, NM)[7]，F(xiàn)an模型(Fan Model, FM)[8]，廣義雙線性模型(Generalized Bilinear Model, GBM)[9]，多項(xiàng)式后驗(yàn)非線性模型(Polynomial Post-nonlinear Model, PPNM)[10]等.其中，由于GBM概括性強(qiáng)，且模型相對簡單，常作為BMM的代表用于非線性解混算法的研究.因此，本文重點(diǎn)關(guān)注Hapke模型和GBM這兩種不同的非線性混合模型的高光譜解混問題.

對于GBM，傳統(tǒng)的線性端元提取算法如頂點(diǎn)成分分析(Vertex Component Analysis, VCA)[11]可以用來提取端元.但是在非線性混合場景下，線性端元提取算法提取端元的不準(zhǔn)確性會對非線性解混結(jié)果有比較大的影響.而對于Hapke模型，則無法利用線性端元提取算法提取端元.非負(fù)矩陣分解(Nonnegative Matrix Factorization, NMF)[12-13]在無監(jiān)督的高光譜線性解混中得到了廣泛應(yīng)用，它可以在端元未知的情況下，同時(shí)對端元和豐度進(jìn)行計(jì)算和求解.在非線性混合情況下，基于核函數(shù)的核非負(fù)矩陣分解(Kernel NMF, KNMF)[14]不僅可以無監(jiān)督地非線性解混高光譜數(shù)據(jù)，而且不需要知道具體的非線性表達(dá)形式，這樣就避免了Hapke模型和GBM中參數(shù)選擇的問題，可以適用于不同的非線性混合模型.目前基于核非負(fù)矩陣分解的高光譜非線性解混算法，主要包括： 不含純像元的核非負(fù)矩陣分解(non-pure-pixels KNMF, npKNMF)[15]，約束的核非負(fù)矩陣分解(Constrained KNMF, CKNMF)[16]，雙目標(biāo)非負(fù)矩陣分解(Bi-objective NMF)[17]等.npKNMF直接利用核非負(fù)矩陣分解進(jìn)行高光譜非線性解混，并沒有考慮到實(shí)際地物的空間分布特性，容易陷入到局部極小.CKNMF在npKNMF的基礎(chǔ)上添加了端元最小體積約束和豐度平滑約束，但是，端元最小體積的特性在高維特征空間中是否成立，是值得懷疑的.并且，它使用了豐度的2范數(shù)來實(shí)現(xiàn)豐度的平滑約束，這并不能反映豐度的局部平滑特性.Bi-objective NMF同時(shí)考慮原始空間和高維特征空間，但是二者所占的權(quán)重不確定，且沒有添加約束，也容易陷入局部極小.

為了解決上述非線性解混算法中存在的問題，本文提出一種基于豐度約束核非負(fù)矩陣分解的非線性光譜解混算法(Abundance Sparseness and Smoothness Constrained KNMF, ASSKNMF).該算法主要由核非負(fù)矩陣分解和豐度約束兩部分組成.首先，利用核非負(fù)矩陣分解將原始數(shù)據(jù)映射到高維空間進(jìn)行無監(jiān)督的高光譜解混.在無監(jiān)督解混的基礎(chǔ)上，依據(jù)地物分布的空間特性，對豐度添加稀疏約束和局部平滑約束，使陷入局部極小的概率變小.實(shí)驗(yàn)結(jié)果表明，所提出的算法可以提高Hapke模型和GBM這兩種非線性混合模型的解混精度.

1 問題描述

對于高光譜數(shù)據(jù)X=[x1,x2,…,xN]∈L×N，其每列xj∈L×1(j=1,2,…,N)都對應(yīng)一個(gè)具有L個(gè)波段的像元向量，共有N個(gè)這樣的像元.以A=[a1,a2,…,aP]∈L×P表示端元矩陣(P是端元數(shù))，sij是端元ai在像元xj中的豐度，ej是噪聲誤差.LMM假設(shè)下的像元向量是端元光譜的線性組合，考慮到實(shí)際的物理意義，豐度需滿足非負(fù)與“和為一”的約束條件[2]：

(1)

1.1 Hapke混合模型

在沙地、礦物緊密混合地區(qū)的遙感圖像數(shù)據(jù)中，依據(jù)輻射傳輸理論[6]，在粒子為球形且各向散射同性，表面反照比低的假設(shè)下，地物介質(zhì)表面雙向反射率r可以簡化地表示為：

(2)

其中：μ0,μ分別是光線入射角和出射角的余弦值；ω是單次散射反照率；H是多項(xiàng)散射函數(shù)，表達(dá)式是：

(3)

(4)

根據(jù)式(4)，可以在已知反射率x的情況下推導(dǎo)出單次散射反照率ω，R-1表示反函數(shù)：

(5)

根據(jù)Hapke模型，將像元的反射率x非線性映射到單次反照率ω后，可以用端元單次反照率的線性組合來表示：

(6)

其中：ωi是端元的單次散射反照率；sij是對應(yīng)的豐度系數(shù).端元在單次反照率按豐度系數(shù)線性組合后，經(jīng)過非線性映射變?yōu)楦吖庾V數(shù)據(jù)中的反射率.

1.2 GBM混合模型

針對植被覆蓋、建筑林立地區(qū)的非線性混和現(xiàn)象，GBM在LMM的基礎(chǔ)上考慮了兩兩端元間的二次散射項(xiàng)作為非線性效應(yīng)，它的表達(dá)式如下：

(7)

其中：ai·ak=[i,1ak,1,ai,2ak,2,…,ai,Lak,L]T是端元之間的Hadamard乘積，用來描述光線在兩種端元物質(zhì)間發(fā)生的二次散射；參數(shù)γi,k,j用于控制非線性程度，使得GBM可以更加靈活的描述非線性混合效應(yīng).

1.3 非負(fù)矩陣分解

非負(fù)矩陣分解在模式識別和機(jī)器學(xué)習(xí)等方面取得了廣泛的應(yīng)用[12]，它可以從給定的非負(fù)矩陣中分解出兩個(gè)非負(fù)矩陣，使得兩個(gè)非負(fù)矩陣的乘積與給定矩陣盡量保持一致.沿用高光譜數(shù)據(jù)的表示符號，非負(fù)矩陣分解將矩陣X∈L×N分解為A∈L×P和S∈P×N，P是矩陣A和S的秩，它的目標(biāo)函數(shù)是基于歐式距離構(gòu)造的：

(8)

通過最小化該目標(biāo)函數(shù)，可以求得矩陣A和S.Lee和Seung提出了一種迭乘算法來對A和S進(jìn)行求解[18]：

A←A.*(XST)./(ASS),

(9)

S←S.*(ATX)./(ATAS).

(10)

其中： 矩陣A是端元矩陣，每列對應(yīng)于一個(gè)端元光譜；矩陣S是豐度矩陣，每列對應(yīng)于一個(gè)像元的豐度；.*和./分別對應(yīng)于矩陣的點(diǎn)乘和點(diǎn)除.

2 方法描述

本文利用豐度約束的核非負(fù)矩陣分解進(jìn)行無監(jiān)督的高光譜非線性解混.首先，利用核方法將原始高光譜數(shù)據(jù)映射到高維特征空間，在高維空間中用非負(fù)矩陣分解進(jìn)行線性解混，同時(shí)根據(jù)實(shí)際地物的分布特性，在豐度上添加稀疏和平滑約束，更好地改善解混效果.

2.1 核方法

考慮非線性映射：

φ:x∈L→φ(x)∈CM.

(11)

將原始的高光譜數(shù)據(jù)映射到高維空間中，使得數(shù)據(jù)在高維空間中變得線性可分.但通常情況下，非線性映射φ是很復(fù)雜和難以求得的，可以通過核函數(shù)的方法來避免這個(gè)問題[19].兩個(gè)向量在高維空間中的內(nèi)積可以用核函數(shù)的形式表示：

k(x,z)=〈φ(x),φ(z)〉x,z∈L.

(12)

根據(jù)式(12)可以定義核矩陣，對于一個(gè)向量集{x1,x2,…,xl}，核矩陣K定義為一個(gè)l×l的矩陣，它的矩陣元素為Kij=〈φ(xi),φ(xj)〉=k(xi,xj)，是一個(gè)對稱矩陣.一般情況下，滿足Mercer定理或正定性的函數(shù)可以作為核函數(shù)，本文采用RBF函數(shù)[20]：

(13)

其中σ是核參數(shù).

2.2 核非負(fù)矩陣分解

針對本文前面所提及的Hapke模型和GBM兩種不同的非線性混合模型，使用核非負(fù)矩陣分解將原始高光譜數(shù)據(jù)映射到高維空間，再利用非負(fù)矩陣分解進(jìn)行無監(jiān)督解混，而不需要考慮具體的非線性表達(dá)形式.核非負(fù)矩陣分解的表達(dá)形式是：

(14)

其中：Xφ=[(x1),φ(x2),…,φ(xN)]；Aφ=[(a1),φ(a2)，…，φ(aP)].通過非線性映射分別將像元和端元映射到特征空間，而豐度矩陣S反應(yīng)真實(shí)地物分布，保持不變.根據(jù)式(12)和核矩陣的定義，式(14)可以表示為：

(15)

2.3 豐度約束

除了利用核非負(fù)矩陣分解處理非線性混合效應(yīng)外，還應(yīng)該考慮實(shí)際情況中地物的空間分布特性，如其所具有的稀疏特性和平滑特性，通過添加適當(dāng)?shù)募s束使陷入局部極小的概率變小.在實(shí)際的高光譜圖像中，單個(gè)像元xj一般情況下只由所有端元中的幾個(gè)端元混合而成，體現(xiàn)在豐度向量上就是sj只有少數(shù)元素不為0，這就具有稀疏的性質(zhì).錢沄濤等人提出的1/2范數(shù)約束[21]可以很好地表達(dá)這一性質(zhì)，它的表達(dá)式是：

(16)

表示豐度矩陣中每個(gè)元素的平方根的和.除了稀疏性質(zhì)之外，豐度還具有局部平滑的性質(zhì)，是因?yàn)樵谡鎸?shí)的場景中，局部區(qū)域內(nèi)地物的分布具有連續(xù)性，體現(xiàn)在豐度上就是具有一定的空間平滑性，在文獻(xiàn)[22]中，通過最小化像素和其鄰域像素豐度的差值來體現(xiàn)局部平滑性.其表達(dá)式如下：

(17)

其中：Spk是第p個(gè)端元在第k個(gè)像素的豐度值；Npk表示它的鄰域.將該端元的豐度向量轉(zhuǎn)換為二維圖像形式;I和J是圖像的尺寸，滿足N=I×J;Sijp是轉(zhuǎn)換后對應(yīng)的豐度值;Nijp是對應(yīng)的鄰域，通過函數(shù)g(x)來衡量豐度值的差異：

(18)

其中： e是自然常數(shù)；x的絕對值越小，g(x)的函數(shù)值越小.g(x)函數(shù)值的最小化等價(jià)于中心像素和鄰域像素豐度差值的最小化，可以起到局部平滑約束的目的.

2.4 ASSKNMF算法

綜合核非負(fù)矩陣分解和豐度約束，ASSKNMF算法的目標(biāo)函數(shù)可以表示為：

(19)

與式(8)中的問題類似，式(19)可以分為兩個(gè)子問題：

(21)

其中k是迭代次數(shù).由于核函數(shù)的存在，式(20)中端元矩陣A的更新不能使用式(9)中的迭乘方法，本文中采用投影梯度法(Projected Gradient, PG)[23]對端元矩陣A進(jìn)行更新，它的表達(dá)式是：

Ak+1=max(Ak-ηk)Af(Ak,Sk,0),

(22)

(23)

尺度因子ηk的計(jì)算采用非負(fù)矩陣分解中用到的回溯Armoji搜索法[24]，在每個(gè)迭代過程中，目標(biāo)函數(shù)值都滿足一定的下降：

f(Ak+1,Sk)-f(Ak,Sk)≤γηkvec(Af(Ak,Sk))Tvec(Ak+1-Ak),

(24)

(25)

其中： vec(g)是將矩陣轉(zhuǎn)為向量；γ是下降程度，通常設(shè)為0.01.

豐度矩陣的求解依然采用式(10)中使用的迭乘更新法，不同的是需要在高維特征空間中進(jìn)行，并且考慮稀疏約束和平滑約束：

(26)

其中： 矩陣Aφ和Xφ是端元矩陣和像元矩陣在高維空間的映射，矩陣S是豐度矩陣；.*和./分別對應(yīng)于矩陣的點(diǎn)乘和點(diǎn)除，h(S-SN)，g′(S-SN)的形式與g(S-SN)相似，具體如下：

(27)

(28)

綜合端元矩陣A和豐度矩陣S的求解過程，整個(gè)算法的流程圖如表1所示：

表1 ASSKNMF算法流程

關(guān)于端元矩陣A和豐度矩陣S的初始化問題，本文中A的初始化采用基于光譜信息散度(Spectral Information Dirergence, SID)[25]的方法.光譜信息散度SID的定義如下：

(29)

其中：xi和xj是任意兩個(gè)像元，它們歸一化的形式是：

(30)

其中L是波段數(shù).

端元初始化的具體步驟是：

a) 隨機(jī)從圖像中選取一個(gè)像元作第一個(gè)端元；

b) 計(jì)算第一個(gè)端元和其他像元的SID值，本文中SID的閾值設(shè)為0.05，第一個(gè)滿足閾值的像元作第二個(gè)端元；

c) 計(jì)算剩余像元和第一、第二個(gè)端元的SID，第一個(gè)滿足閾值的像元選作第三個(gè)端元，依次類推，找出所有端元.

考慮到初始值對非負(fù)矩陣分解的結(jié)果有較大影響這個(gè)因素，SID的閾值取為0.05，使得初始化的端元之間具有一定的不相似性，豐度矩陣S的初始化采用FCLS算法[26].

3 實(shí)驗(yàn)結(jié)果與分析

本節(jié)將提出的ASSKNMF算法與基于高斯核的核非負(fù)矩陣分解算法npKNMF[15]，端元和豐度約束的核非負(fù)矩陣分解算法CKNMF[16]，雙目標(biāo)非負(fù)矩陣分解算法Bi-objective NMF[17]進(jìn)行性能比較，同時(shí)將本文所用的稀疏約束和平滑約束與非負(fù)矩陣分解結(jié)合作為一種線性方法，記為CNMF. Bi-objective NMF中涉及到一個(gè)權(quán)重系數(shù)θ，這里選取0.2,0.5,0.8三種情況作為比較.評價(jià)指標(biāo)主要有端元的光譜角距離(Spectral Angle Distance)DSA和豐度的均方根誤差(Root Mean Square Error)ERMS：

(31)

(32)

3.1 模擬數(shù)據(jù)實(shí)驗(yàn)

從美國地質(zhì)調(diào)查局(United States Geological Survey, USGS)光譜庫中選取了6種不同地物的光譜作為端元生成模擬數(shù)據(jù)，該光譜庫中地物光譜的波長在0.38～2.5μm之間，分辨率10nm，共有224個(gè)波段，圖1顯示了這6種端元的光譜曲線.針對Hapke和GBM這兩種非線性混合模型生成數(shù)據(jù)的方式如下： 1) 首先生成豐度，豐度生成方式和文獻(xiàn)[21]中類似，將36×36大小的圖像分割為互不重疊的大小為6×6的正方形塊，并將這些塊隨機(jī)分配給不同的端元；用7×7的低通濾波器對地物的豐度進(jìn)行濾波，產(chǎn)生混合像元；豐度值超過純像元程度時(shí)，該像元的豐度值由所有端元平均.圖2顯示了6個(gè)端元對應(yīng)的豐度圖.2) 生成Hapke模型數(shù)據(jù)時(shí)，先將端元光譜反射率按照式(5)轉(zhuǎn)為單次散射反照率，然后與1)中產(chǎn)生的豐度相乘得到混合的單次散射反照率，再按式(4)變?yōu)榛旌瞎庾V反射率數(shù)據(jù).3) 生成GBM數(shù)據(jù)時(shí)，將1)中產(chǎn)生的豐度值帶入式(7)產(chǎn)生雙線性數(shù)據(jù)，參數(shù)γ在[0,1]中隨機(jī)取值.生成模擬數(shù)據(jù)后，再添加不同信噪比的高斯白噪聲.

圖1 USGS光譜庫中的6條光譜Fig.1 Six endmember spectra from USGS spectral library

圖2 6個(gè)端元的豐度圖Fig.2 Abundance maps of six endmembers

實(shí)驗(yàn)1參數(shù)確定實(shí)驗(yàn)： 本文提出的算法中主要有稀疏約束正則參數(shù)α，平滑約束正則參數(shù)β，核參數(shù)σ以及平滑程度參數(shù)γm.其中平滑程度參數(shù)γm按照文獻(xiàn)[22]的取法取0.5，剩余3個(gè)參數(shù)的選取通過實(shí)驗(yàn)來決定.第一步固定核參數(shù)σ=6，變化α和β來觀察結(jié)果.圖像大小是36×36，6個(gè)端元，純像元程度為0.8，加入40dB的高斯白噪聲，分別針對Hapke模型和GBM進(jìn)行實(shí)驗(yàn)，取10次試驗(yàn)的平均結(jié)果.圖3顯示了α和β對Hapke模型和GBM結(jié)果的影響.從圖上可以看出，在兩種不同的非線性模型下，α和β分別取0.001和0.1，可以使DSA和ERMS達(dá)到比較好的結(jié)果.第二步固定α和β的值，變化σ來觀察結(jié)果.實(shí)驗(yàn)設(shè)置不變，取10次試驗(yàn)的平均結(jié)果.

圖3 α和β對結(jié)果的影響

圖4顯示了σ對兩個(gè)模型結(jié)果的影響，可以看出σ在區(qū)間[3,10]上時(shí)結(jié)果比較好，經(jīng)過權(quán)衡后，最終σ選定為4.

圖4 σ對結(jié)果的影響Fig.4 Results with respect to σ

實(shí)驗(yàn)2端元數(shù)目和信噪比的影響分析： 在該實(shí)驗(yàn)中比較了各算法在不同端元數(shù)目和信噪比下解混的精度.圖像大小是36×36，純像元程度是0.8，端元數(shù)目取3、4、5、6.在端元數(shù)目固定的情況下，分別對Hapke模型和GBM數(shù)據(jù)進(jìn)行解混，信噪比取40dB和20dB.表2～表5(見第436～437頁)分別顯示了端元數(shù)目為3、4、5、6情況下各算法的解混結(jié)果.從表中的數(shù)據(jù)可以看出，在端元提取方面： 隨著端元數(shù)目的增加，非線性程度越復(fù)雜，各算法端元提取的精度會逐漸下降，本文提出的ASSKNMF算法和其他算法相比一直保持最優(yōu)的結(jié)果.僅利用核非負(fù)矩陣分解的npKNMF的解混精度隨著端元數(shù)目的增多下降很明顯，線性的CNMF的解混結(jié)果在端元數(shù)目為3時(shí)還尚可，端元數(shù)目增多后也不理想.CKNMF端元提取的結(jié)果并不理想，這也間接證明了其添加的端元最小體積約束并不能提高解混精度.

表2 端元數(shù)目為3時(shí)各算法結(jié)果比較(×10-2)

注： 每種評價(jià)指標(biāo)的最優(yōu)結(jié)果由粗體標(biāo)出.

表3 端元數(shù)目為4時(shí)各算法結(jié)果比較(×10-2)

注： 每種評價(jià)指標(biāo)的最優(yōu)結(jié)果由粗體標(biāo)出.

表4 端元數(shù)目為5時(shí)各算法結(jié)果比較(×10-2)

注： 每種評價(jià)指標(biāo)的最優(yōu)結(jié)果由粗體標(biāo)出.

表5 端元數(shù)目為6時(shí)各算法結(jié)果比較(×10-2)

注： 每種評價(jià)指標(biāo)的最優(yōu)結(jié)果由粗體標(biāo)出.

Bi-objective NMF通過結(jié)合線性和非線性，結(jié)果總體上比CNMF，npKNMF和CKNMF要好，但是和權(quán)重因子的取值有較大的關(guān)系.在豐度估計(jì)方面： 從數(shù)據(jù)上可以看出沒有添加豐度約束的npKNMF和Bi-objective NMF的豐度估計(jì)精度要低于添加豐度約束的CNMF，CKNMF和ASSKNMF.CKNMF豐度估計(jì)的結(jié)果雖比CNMF要好，但是其使用的2范數(shù)約束是從整體上平滑豐度，在空間分辨率低時(shí)會模糊化豐度圖，這與實(shí)際地物分布時(shí)存在的局部平滑特性是矛盾的，而本文提出的ASSKNMF既考慮了豐度分布的稀疏特性，又考慮了豐度分布的局部平滑特性，使得最后估計(jì)的豐度更接近于真實(shí)值.

3.2 實(shí)際數(shù)據(jù)實(shí)驗(yàn)

本文采用兩個(gè)實(shí)際地區(qū)的高光譜數(shù)據(jù)來分析和評價(jià)各算法的解混結(jié)果.第一個(gè)數(shù)據(jù)是機(jī)載可見光/紅外成像光譜儀(Airborne Visible/Infrared Imaging Spectrometer, AVIRIS)獲取的美國內(nèi)華達(dá)州Cuprite礦物地區(qū)的高光譜圖像，圖像大小是250×191(如圖5(a)所示)，有224個(gè)波段，解混前需剔除掉信噪比低和水吸收波段(1～2,104～113,148～167,221～224)，還剩余188個(gè)波段.第二個(gè)數(shù)據(jù)是HYDICE華盛頓國家廣場Washington DC Mall圖像(如圖5(b)所示)，圖像大小是1280×307，有210個(gè)波段，解混前同樣剔除掉信噪比低和水吸收波段(1～4,76,87,101～111,136～153,198～210)，還剩余162個(gè)波段.Cuprite地區(qū)富含多種礦物質(zhì)，會存在緊密混合現(xiàn)象，HYDICE地區(qū)有建筑、草地、樹木，會存在多層次混合現(xiàn)象，且這兩個(gè)場景有詳細(xì)的地物真實(shí)調(diào)查參考[27]同時(shí)以往也經(jīng)常被用作解混研究的數(shù)據(jù).

圖5 真實(shí)高光譜遙感圖像Fig.5 Real hyperspectral remote sensing images

從第一幅Cuprite圖像中截取大小為50×40的子區(qū)域作為測試數(shù)據(jù)，如圖5(a)所示. 根據(jù)研究和實(shí)際調(diào)查結(jié)果[27]可知，該子區(qū)域中主要包括鈣鐵榴石，高嶺石和蒙脫石3種物質(zhì).表6中列出了各算法對該區(qū)域解混的DSA值，從數(shù)據(jù)可以看出，CKNMF的解混結(jié)果劣于其他所有方法，證明其使用的端元約束的正確性有待考究.npKNMF和Bi-objective NMF兩種非線性方法的結(jié)果要優(yōu)于線性的CNMF，而本文提出的算法兼顧考慮了非線性效應(yīng)和地物分布特性，得出的結(jié)果要優(yōu)于其余的方法.圖6中列出了本文算法求解的3種端元對應(yīng)的豐度圖，能夠比較真實(shí)的反映地物分布.

表6 不同算法在Cuprite數(shù)據(jù)上的端元提取結(jié)果與參考光譜的光譜角距離比較

注： 評價(jià)指標(biāo)的最優(yōu)結(jié)果由粗體標(biāo)出.

對于第二幅HYDICE圖像，從中截取大小為100×100的子區(qū)域用來測試算法的性能，如圖5(b)所示.該地區(qū)主要有5種地物： 樹木，屋頂，水體，道路和草地.表7中列出了各算法對該區(qū)域5種地物解混的DSA值.從結(jié)果可以看出，本文提出的方法ASSKNMF在屋頂，道路和草地上結(jié)果最好，樹木是Bi-objective方法最好，水體是CNMF最好，平均結(jié)果上本文方法是最好的.同時(shí)可以發(fā)現(xiàn)端元數(shù)為5時(shí)，CKNMF結(jié)果最不好且下降明顯，Bi-objective NMF要優(yōu)于CNMF，ASSKNMF算法端元提取的結(jié)果最好.根據(jù)以往用作高光譜解混的研究結(jié)果[22]可知，圖7中列出的5種端元對應(yīng)的豐度圖更接近實(shí)際地物分布.

圖6 Cuprite數(shù)據(jù)ASSKNMF算法提取的3種端元的豐度圖Fig.6 Abundance maps of three endmembers extracted by ASSKNMF of Cuprite data set

算法DSA樹木屋頂水體道路草地平均CNMF0.07780.12120.08020.15560.30820.1486CKNMF0.10060.33560.71090.42750.16060.3470Bi-objectiveNMFθ=0.20.03690.16080.27780.10880.05740.1283θ=0.50.02920.14070.18050.15660.11500.1244θ=0.80.04660.13790.10610.17020.16180.1245ASSKNMF0.05960.11820.14970.09240.05530.0950

注： 評價(jià)指標(biāo)的最優(yōu)結(jié)果由粗體標(biāo)出.

圖7 HYDICE數(shù)據(jù)ASSKNMF算法提取的5種端元的豐度圖Fig.7 Abundance maps of five endmembers extracted by ASSKNMF of HYDICE data set

4 結(jié) 論

提出了一種基于豐度約束核非負(fù)矩陣分解的高光譜非線性解混算法.與現(xiàn)有的基于核非負(fù)矩陣分解的算法相比，本算法先利用核方法將原始高光譜數(shù)據(jù)投影到高維特征空間，在特征空間中利用非負(fù)矩陣分解進(jìn)行解混，同時(shí)考慮到實(shí)際地物的分布特性，在豐度上添加稀疏約束和平滑約束，使陷入局部極小值的概率變小.由于非負(fù)矩陣分解的使用，無需考慮具體的非線性表達(dá)形式，對緊密混合情況下的Hapke模型和多層次混合下的GBM都可以適用.模擬數(shù)據(jù)和真實(shí)遙感圖像的實(shí)驗(yàn)結(jié)果表明，本文提出的算法與其他相關(guān)方法相比，可以提高非線性混合情況下光譜解混的精度.