林潔雯 陳 建 羅婷文 徐志搏

(1.自然資源部城市國土資源監(jiān)測與仿真重點實驗室, 深圳 518000;2.自然資源部超大城市自然資源時空大數(shù)據(jù)分析應(yīng)用重點實驗室, 上海 200063;3.中國農(nóng)業(yè)大學(xué)工學(xué)院, 北京 100083; 4.深圳市智能微小衛(wèi)星星座技術(shù)與應(yīng)用重點實驗室, 深圳 518107;5.浙江大學(xué)能源清潔利用國家重點實驗室, 杭州 310027;6.浙江大學(xué)浙江省農(nóng)業(yè)智能裝備與機(jī)器人重點實驗室, 杭州 310058)

0 引言

遙感技術(shù)已經(jīng)廣泛應(yīng)用于農(nóng)業(yè)遙感[1-4]、森林監(jiān)測、生態(tài)信息提取、海洋監(jiān)測等領(lǐng)域。遙感探測器記錄的能量與瞬時視場角(IFOV)對應(yīng)地面單元內(nèi)地表物質(zhì)反射或發(fā)射的能量有關(guān)。但是由于傳感器分辨率的限制,導(dǎo)致某一像元中包含多種地物,形成混合像元的現(xiàn)象[5]。通過光譜混合分析可以分解瞬時場角內(nèi)由一組固定端元混合而成的單個像元的光譜反射率觀測值,稱為解混。而能夠從中提取出來的每一類地物的光譜特征稱為端元,它決定了解混的精度。在純像元假設(shè)下,一組高光譜影像同時包含純像元和混合像元,純像元中只含有一個端元,混合像元中含有多個端元[6]。人工難以區(qū)分和辨認(rèn)混合像元和純像元,需通過算法自動提取實現(xiàn)解混。自動端元提取兼具便利性及準(zhǔn)確性:①便利性:人工無法確定表達(dá)某類地物信息的端元數(shù)量,而算法可以自動提取端元的數(shù)量。②準(zhǔn)確性:人工選擇沒有標(biāo)準(zhǔn)或指標(biāo)作為依據(jù),而自動端元提取是根據(jù)算法尋優(yōu)策略和指標(biāo)(如:均方根誤差(RMSE))選擇的,提取出的端元能夠更接近實測的地物光譜信息,解豐度再反混后更接近原圖信息,準(zhǔn)確性更強(qiáng)。

遙感影像上同物異譜或同譜異物的現(xiàn)象會造成同一混合像元端元在影像上光譜不唯一,產(chǎn)生光譜變異性[7]。當(dāng)影像中存在光譜變異現(xiàn)象時,如果只用一條端元光譜表示一類物質(zhì)的端元,容易造成誤差。因此,端元束比單條端元光譜更能反映地物信息。近幾年,國內(nèi)外專家和學(xué)者將端元提取研究轉(zhuǎn)向端元束的提取研究。端元束提取方法能夠獲得每一類物質(zhì)具有類內(nèi)光譜變化的若干端元。目前提取端元束的算法相對于端元提取算法較少,現(xiàn)有方法包括分塊提取和集群智能算法。分塊提取方式將截取的像素塊通過隨機(jī)(EBE)[8]、局部提取(LLG)[9]、聚類[10-13]等方法提取端元,通過篩掉奇異值再合并成端元束。左成歡[14]采用PPI[15]提取并經(jīng)過超像素分析保留的端元集,并根據(jù)植被指數(shù)將端元分類,通過聚類方法得到每類地物的端元束。陳立偉等[16]采用均值漂移算法完成聚類并提取端元束。上述方法依賴分塊操作或預(yù)分類的結(jié)果,對于沒有純像元的高光譜影像容易提取出錯誤的端元。在集群智能算法方面,文獻(xiàn)[17-18]提出了IQPSO和MODPSO。IQPSO通過極端重構(gòu)誤差判斷所提取的端元是否是無效解,并通過端元函數(shù)判斷是否結(jié)束迭代。該算法通過設(shè)定閾值控制算法迭代次數(shù),容易出現(xiàn)迭代冗余,算法效率較低。MODPSO將端元提取作為多目標(biāo)優(yōu)化分析,采用0/1圖像編碼方式和速度及位置更新方式。該方法獲取的端元束數(shù)量較少,無法準(zhǔn)確表達(dá)地物的光譜信息。綜上所述,采用集群智能算法效率低的主要原因在于:粒子編碼、更新方式及迭代機(jī)制耗時較多。而提取端元數(shù)量少的原因在于:粒子迭代更新忽略了決策空間的分布特性,陷入了局部最優(yōu)。因此,在端元束提取算法的設(shè)計上需要兼顧同類地物光譜的變化、提取精度、效率等,盡可能多地獲取有效端元束。

本文為解決現(xiàn)有方法存在的問題,通過多目標(biāo)多模態(tài)任務(wù)實現(xiàn)端元束提取。首先,對光譜圖像進(jìn)行一維標(biāo)號編碼,采用改進(jìn)粒子群速度及位置更新方式更新粒子迭代信息。然后,引入環(huán)形拓?fù)浣Y(jié)構(gòu)改善粒子間的信息傳遞,提升信息搜索能力,并通過擁擠距離排序改善決策空間和目標(biāo)空間的擁擠距離,最后,將端元束提取問題轉(zhuǎn)換到多模態(tài)多目標(biāo)的問題上,提高算法效率的同時實現(xiàn)端元束提取的有效性。

1 線性光譜混合模型與優(yōu)化

光譜解混[19-20]依賴于場景混合模型,大尺度的光譜混合可認(rèn)為是線性混合問題[21]。本文主要圍繞線性混合模型進(jìn)行混合像元問題的探究,并在此基礎(chǔ)上提出基于多模態(tài)多目標(biāo)優(yōu)化的端元束提取方法。

1.1 線性光譜混合模型

(1)

式中ej——端元

aij——端元ej在像元ri中的權(quán)重,即豐度

εj——誤差項

q——端元數(shù)量

由于aij的物理意義是端元在像元中的面積比,需滿足“非負(fù)約束”以及“和為1約束”兩個條件,即

aij>0

(2)

(3)

本文將LSMM模型的求解方法建立在基于純像元假設(shè)的基礎(chǔ)上,所提取的端元均源于高光譜影像R中的像元,通過搜索像元提取端元束。

1.2 多模態(tài)多目標(biāo)優(yōu)化

多目標(biāo)優(yōu)化問題定義為

(4)

式中o——待優(yōu)化的目標(biāo)數(shù)

x——d維決策變量

gz(x)——不等式約束

hj(x)——等式約束

k——不等式約束數(shù)量

s——等式約束數(shù)量

決策變量需滿足k個不等式約束gz(x)≤0,z=1,2,…,k,以及s個等式約束hj(x)=0,j=1,2,…,s。同時滿足以上約束的兩個空間分別為決策空間R(d)和決策空間映射出的目標(biāo)空間R(o)。因此,在給定一個包含o個目標(biāo)的優(yōu)化問題時,當(dāng)且僅?fi(x1)≤fi(x2)∧?fj(x1)

多模態(tài)多目標(biāo)問題如圖1所示,以待優(yōu)化目標(biāo)數(shù)o=2為例,點p1、p2對應(yīng)的目標(biāo)空間值都為f(x1)。多模態(tài)特點就在于解的形態(tài)是多樣的,即在保證目標(biāo)值最優(yōu)的情況下可以找到不同的解。一方面可以更深入地了解該組合解的分布規(guī)律。另一方面,混合像元分解的目的在于求解出像元所含地物的光譜特征。而光譜變異性導(dǎo)致了同類地物在不同像元的光譜特征存在的差異。該情況下,沒有絕對的最優(yōu)解,決策空間提供多組解應(yīng)對光譜特征的變化。

圖1 多模態(tài)多目標(biāo)問題示意圖

通常,多模態(tài)多目標(biāo)問題滿足以下條件之一[23]:①至少有一個局部帕累托最優(yōu)解。②至少有兩個等效全局帕累托最優(yōu)解,它們對應(yīng)PF上同一點。對滿足多模態(tài)多目標(biāo)條件的高光譜影像進(jìn)行端元束提取。

高光譜影像特征復(fù)雜,雖然集群智能算法能夠?qū)崿F(xiàn)端元束提取[17-18],但在尋優(yōu)過程中未考慮到PS的分布,無法在廣泛空間中尋優(yōu),導(dǎo)致獲得有效端元數(shù)量較少。本文采用多模態(tài)多目標(biāo)解決端元束提取問題主要關(guān)注兩方面:①提高算法搜索能力,搜尋較多的帕累托最優(yōu)解。②盡量保留目標(biāo)空間距離較小而決策空間較大的解,同時保證決策空間和目標(biāo)空間的多樣性。

2 基于改進(jìn)離散粒子群的端元束提取

集群智能算法通過簡單信息處理單元之間的信息交流產(chǎn)生一種解決問題的能力。通常包括種群編碼及初始化、種群進(jìn)化與更新、篩選下一代種群、循環(huán)變異直至滿足終止條件4個主要步驟。粒子群優(yōu)化算法(Particle swarm optimization, PSO)是集群智能算法的主要代表,仿照鳥群覓食進(jìn)行算法設(shè)計,主要用于解決可行解空間連續(xù)的單目標(biāo)優(yōu)化問題。本文所探究的可行解空間是離散集合,且從多模態(tài)多目標(biāo)角度提取端元束。因此,在采用PSO算法的基礎(chǔ)上,均需要對粒子編碼、速度與位置更新、決策空間與目標(biāo)空間排序及目標(biāo)函數(shù)方面進(jìn)行調(diào)整及優(yōu)化,提出一種多模態(tài)多目標(biāo)優(yōu)化的端元束提取方法(Multi-modal and multi-objective particle swarm optimization by special crowding distance, MOPSOSCD)。

2.1 粒子編碼

對于大小為n×h×l的高光譜影像R,以圖像從左到右、從上至下對像元標(biāo)號,并對圖像中的像素位置編碼,像素點編號集合Z為

Z={z1,z2,…,zn×h}

(5)

通過該編碼方式,粒子的位置編碼還原至影像中形成多組端元解。其中,進(jìn)行種群進(jìn)化時某一組可行解為

X={x1,x2,…,xq} (xi∈Z)

(6)

式中X——粒子像元對應(yīng)的標(biāo)號組成的序列

每個X對應(yīng)一組端元E=[e1e2…eq]T。粒子初始化時,按照設(shè)定的粒子數(shù)隨機(jī)生成矩陣E進(jìn)行速度與位置更新及迭代。

2.2 改進(jìn)粒子速度與位置更新

PSO算法通過迭代尋優(yōu),在每次迭代中通過跟蹤個體最優(yōu)xpbest和全局最優(yōu)xgbest更新自身的速度和位置,粒子i第t+1代的速度和位置更新分別為

vi(t+1)=ωvi(t)+c1(xpbest,i(t)-xi(t))+

c2(xgbest,i(t)-xi(t))

(7)

xi(t+1)=xi(t)+vi(t+1)

(8)

式中vi(t+1)——粒子的更新速度

ω——慣性系數(shù)

c1、c2——學(xué)習(xí)因子

xi(t+1)——粒子的更新位置

本文的高光譜影像編碼方式導(dǎo)致粒子為離散形式更新,需在每個粒子位置更新完成后對該值取整。同時,為了保證所更新的粒子位置位于尋優(yōu)范圍內(nèi),將式(8)改進(jìn)為

(9)

式中 ceil()——取整函數(shù)

xr(t+1)——區(qū)間Z內(nèi)的隨機(jī)數(shù)

2.3 環(huán)形拓?fù)浣Y(jié)構(gòu)與擁擠距離

在多模態(tài)多目標(biāo)優(yōu)化中,粒子的選擇極為關(guān)鍵,本文采用MO_Ring_PSO_SCD[24]方法進(jìn)行粒子信息傳遞和特殊擁擠距離計算。針對信息傳遞,設(shè)置個體最優(yōu)xpbest的存檔PBA和鄰域最優(yōu)xnbest的存檔NBA,在鄰域最優(yōu)和子代最優(yōu)選擇中采用環(huán)形拓?fù)涓路绞?其更新規(guī)則為:對于初始化或者已經(jīng)進(jìn)行迭代更新的粒子群,當(dāng)粒子j為1～P(P表示粒子數(shù)量)時,拓?fù)涓卤憩F(xiàn)為粒子j僅與j-1和j+1粒子進(jìn)行信息交互,當(dāng)j=1時,拓?fù)涓卤憩F(xiàn)為粒子j與P和第2個粒子進(jìn)行信息交互,當(dāng)j=P時,粒子j與第1個和P-1粒子進(jìn)行信息交互,以此形成閉環(huán)拓?fù)浣Y(jié)構(gòu)信息傳遞。因此,將式(7)更新方法調(diào)整為

vi(t+1)=ωvi(t)+c1(xpbest,i(t)-xi(t))+

c2(xnbest,i(t)-xi(t))

(10)

對于多模態(tài)多目標(biāo)優(yōu)化算法,需同時考慮決策空間和目標(biāo)空間中個體的優(yōu)劣問題。特殊擁擠距離Si可表示為

(11)

式中Ci,x——第i個粒子的決策距離

Cavg,x——決策距離平均值

Ci,f(x)——第i個粒子的目標(biāo)空間距離

Cavg,f(x)——目標(biāo)空間距離平均值

決策空間在空間分布上為離散問題,其數(shù)值范圍為區(qū)間Z內(nèi)的正整數(shù)。中間粒子的決策距離表示為

(12)

其中

d=q

式中d——決策空間維度

δ——決策空間加權(quán)項

通常為保證決策空間計算的有效性,令δ=1。如圖2a所示,以二維決策空間計算為例,如獲取粒子②在決策空間中的距離,則通過左右粒子①、③以及空間中最大最小值計算得

圖2 決策空間與目標(biāo)空間擁擠距離計算方法示意圖

(13)

如果計算邊界粒子的決策距離,邊界粒子i=P時計算方法為

(14)

粒子④在決策空間中的距離可通過自身及臨界粒子計算,計算式為

(15)

同理,邊界粒子j=1時計算方法為

(16)

目標(biāo)空間在空間分布上可認(rèn)為是連續(xù)問題,由于目標(biāo)函數(shù)值的區(qū)間不同,通常將目標(biāo)空間值歸一化后進(jìn)行距離計算。計算中間粒子的目標(biāo)空間距離方法與決策空間類似,即

(17)

如圖2b所示,粒子②在決策空間中的距離同樣通過左右粒子①、③以及空間中最大、最小值計算得到,計算方法為

(18)

計算目標(biāo)空間的邊界點與決策空間不同,本文解混問題為最小優(yōu)化問題,當(dāng)粒子某目標(biāo)值最小時,將該方向的值設(shè)置為1,反之設(shè)置為0。例如:粒子④在目標(biāo)空間中的距離為C4,f(x)=0+1=1。

2.4 目標(biāo)函數(shù)

(19)

2.5 算法流程

基于多模態(tài)多目標(biāo)優(yōu)化的端元束提取方法流程如圖3所示。結(jié)合上文所描述的MOPSOSCD方法,具體步驟為:①輸入大小為n×h×l的高光譜數(shù)據(jù)R,按照數(shù)據(jù)從上至下,從左至右標(biāo)號編碼,進(jìn)行粒子初始化并設(shè)置相關(guān)參數(shù):定向選擇概率pm、粒子數(shù)量P及迭代次數(shù)M。②采用目標(biāo)函數(shù)評估初始化粒子的兩個目標(biāo)函數(shù)值。③通過環(huán)形拓?fù)浣Y(jié)構(gòu)更新NBA,使粒子與其鄰近交換信息。④根據(jù)式(11)計算NBA的距離信息S:包含式(12)、(14)、(16)所計算的決策空間距離及式(17)所計算的目標(biāo)空間距離,再按照距離信息S對NBA排序。⑤按照定向選擇概率pm隨機(jī)產(chǎn)生更新索引并判斷PBA的首行元素是否在更新索引中:若是,按照式(9)生成隨機(jī)數(shù)xr,執(zhí)行步驟⑥;若否,xi保持不變,執(zhí)行步驟⑥。⑥選擇NBA及PBA首行作為xnbest和xpbest進(jìn)行速度與位置更新。⑦評估位置更新粒子的兩個目標(biāo)函數(shù)值。⑧更新PBA:使粒子與其鄰近交換信息。⑨根據(jù)式(11)計算PBA的距離信息SCD:包含決策空間距離、目標(biāo)空間距離,再按照距離信息SCD對PBA排序。⑩判斷是否達(dá)到收斂極限或達(dá)到最大迭代次數(shù)M;若是,執(zhí)行步驟;若否,返回執(zhí)行步驟④。輸出NBA中的非支配解集,提取解集中的候選端元束,即為非支配解組成帕累托最優(yōu)解集PS,并得到相應(yīng)帕累托前沿PF。

圖3 基于多模態(tài)多目標(biāo)優(yōu)化的端元束提取方法流程圖

3 試驗結(jié)果與分析

3.1 試驗數(shù)據(jù)

如圖4a所示,采用MUUFL高光譜數(shù)據(jù)(https://github.com/GatorSense/MUUFLGulfport)進(jìn)行本文算法驗證[25-26]。MUUFL數(shù)據(jù)搭載高光譜儀的有人機(jī)飛行高度為1 066.8 m,空間分辨率為1 m×1 m,該影像共包括64個波段,覆蓋波長443.9～967.2 nm,光譜分辨率為19.1 nm,本文選取其中圖像尺寸為90像素×130像素的ROI區(qū)域進(jìn)行試驗驗證[27]。如圖4b、4c所示,ROI區(qū)域內(nèi)含1#屋頂、2#草、3#樹、4#陰影、5#瀝青5種端元。

3.2 評價指標(biāo)

為驗證本文算法的有效性,除了采用基于無約束最小二乘法(ucls)解豐度并反混后與原圖的RMSE評價指標(biāo)外,還采用光譜角距離(SAD)進(jìn)行評價。光譜角表示兩個光譜向量之間的光譜角信息,計算式為

(20)

式中b1、b2——光譜向量

將提取端元與實測光譜代入計算,兩者相似程度越高,SAD越小,則說明端元提取效果越好。

3.3 試驗及討論

為驗證本文方法有效性,設(shè)置MOPSOSCD算法定向選擇概率pm=0.2、粒子數(shù)量P=30、迭代次數(shù)M=400。不同算法在MUUFL數(shù)據(jù)提取端元效果如圖5所示,對于 1#屋頂、2#草、3#樹、4#陰影、5#瀝青獲取的端元數(shù)量分別為5、7、7、7、5。

本文方法基于純像元假設(shè)進(jìn)行端元束提取,像素標(biāo)注結(jié)果如圖6所示。算法考慮到光譜變異性,從多模態(tài)多目標(biāo)的角度實現(xiàn)更大范圍尋優(yōu)。因地物重疊、分散造成的混合像元較多,MUUFL數(shù)據(jù)的環(huán)境相對復(fù)雜。通過算法尋優(yōu)結(jié)果可發(fā)現(xiàn),對于1#屋頂、2#草、3#樹、4#陰影端元采用MOPSOSCD算法能夠?qū)崿F(xiàn)大區(qū)域內(nèi)搜索,盡管圖6中5#瀝青分布較分散,MOPSOSCD算法的提取更傾向于尋優(yōu)到5#瀝青的純像元信息。由于MOPSOSCD屬于多目標(biāo)算法,個別端元存在錯分現(xiàn)象?？傮w來說算法在保證擴(kuò)大搜索空間,避免陷入局部最優(yōu)的效果較好。

圖6 MOPSOSCD算法端元提取結(jié)果

多模態(tài)多目標(biāo)問題面臨確定最終解[23],對于MOPSOSCD算法提取的端元束,采用基于無約束最小豐度的ISMA[28]算法逐個像素提取最優(yōu)端元組合,計算評價指標(biāo)。同時,MOPSOSCD算法所提取端元束中每一組解都可以作為一組端元提取結(jié)果,將MOPSOSCD算法所得端元束中最小mSAD(min mSAD)的端元組合,用于計算評價指標(biāo)。以VCA(頂點成分分析法)[29]和同樣基于迭代方式進(jìn)行端元提取的DPSO[30]作為對比方法,分別在MUUFL數(shù)據(jù)上進(jìn)行對比試驗。VCA算法以線性光譜混合模型的幾何學(xué)描述為基礎(chǔ),通過反復(fù)尋找正交子空間并計算圖像矩陣在正交子空間中的投影向量二范數(shù)并逐次提取端元。不同算法在MUUFL數(shù)據(jù)上的評價結(jié)果如表1所示。

表1 MUUFL數(shù)據(jù)精度對比結(jié)果

由表1可知,VCA、DPSO及MOPSOSCD 3種算法均能提取5種地物的端元信息。MOPSOSCD算法采用ISMA獲取的最終解為每個像元分配了最優(yōu)的端元組合,具有最優(yōu)RMSE。由于MOPSOSCD算法提取的31條端元均參與了每個像元最優(yōu)端元組合的求解,無法具體求解mSAD及5個端元對應(yīng)的SAD,因此,31條端元的mSAD值為0.111 2。端元束中最小mSAD(min mSAD)所構(gòu)成的端元組合的RMSE及mSAD分別為0.069 0和0.065 3,均優(yōu)于VCA及DPSO。其中,RMSE指標(biāo)如圖7所示。

圖7 MUUFL數(shù)據(jù)均方根誤差

由圖7a可知,MOPSOSCD算法所提取的端元束被充分分配至每個像素中,說明每個像素的RMSE達(dá)到了相對最優(yōu)值,同時驗證了該數(shù)據(jù)光譜變異性的存在及端元束提取的必要性。通過圖7b可知,在采用最小mSAD(min mSAD)所構(gòu)成的端元組合時,對1#屋頂、2#草、3#樹、4#陰影4種端元提取的效果較好,而RMSE較高的區(qū)域主要分布在5#瀝青端元上。從VCA和DPSO的均方根誤差(圖7c、7d)中可發(fā)現(xiàn),VCA算法均方根誤差較大,無法有效區(qū)分各端元信息。DPSO算法提取的屋頂邊緣以及具有植被分布的像元具有較高的均方根誤差,尤其對2#草與3#樹混合部分端元提取效果較差。因此,本文算法相對其他方法在地物邊緣或地物混合區(qū)域的端元提取效果較好。mSAD指標(biāo)對應(yīng)各端元的SAD值如表2所示。

表2 各端元的SAD值

通過表2可發(fā)現(xiàn),MOPSOSCD算法所提取的端元束最小mSAD達(dá)到最優(yōu)效果,尤其對2#草和4#陰影端元提取效果優(yōu)于VCA及DPSO算法。在1#屋頂和3#樹的SAD與DPSO算法接近。采用VCA提取的端元與真實光譜差異性較大,造成各端元SAD較大。

進(jìn)一步的,本文對比了3種算法的運(yùn)行時間。當(dāng)采用參數(shù)組合pm=0.2、P=30、M=400進(jìn)行迭代時,MOPSOSCD算法耗時32 741 s。相同參數(shù)組合的DPSO算法運(yùn)行耗時80 123 s,VCA算法運(yùn)行耗時3.47 s。對于經(jīng)典的VCA算法來說,其優(yōu)勢在于,相對DPSO和MOPSOSCD 2種算法耗時較短。DPSO和MOPSOSCD 2種算法更加傾向于尋優(yōu)效果和精度的提升。相較于DPSO,MOPSOSCD可以快速地提取更多可供決策的候選端元,構(gòu)成端元束。

4 結(jié)束語

提出了一種基于多模態(tài)多目標(biāo)優(yōu)化的端元束提取方法。首先,將高光譜圖像進(jìn)行標(biāo)號編碼,并通過設(shè)定鄰域最優(yōu)存檔NBA及個體最優(yōu)存檔PBA對粒子群速度及位置更新方式進(jìn)行改進(jìn)優(yōu)化。同時,采用基于索引的環(huán)形拓?fù)浣Y(jié)構(gòu)實現(xiàn)不同鄰域的個體進(jìn)行信息交互,防止粒子隨機(jī)更新,并改進(jìn)決策空間擁擠距離,提高粒子在決策空間的多樣性。最后,在MUUFL數(shù)據(jù)上驗證算法有效性。當(dāng)設(shè)置參數(shù)組合為pm=0.2、P=30、M=400時,MOPSOSCD算法的RMSE為0.008 8,提取的端元束的mSAD為0.111 2。該算法相較于其他算法在保證效率的同時獲得更多有效端元束,且具有更高的精度。