基于多模態多目標優化的端元束提取方法研究

林潔雯 陳 建 羅婷文 徐志搏

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

0 引言

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

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

本文為解決現有方法存在的問題,通過多目標多模態任務實現端元束提取。首先,對光譜圖像進行一維標號編碼,采用改進粒子群速度及位置更新方式更新粒子迭代信息。然后,引入環形拓撲結構改善粒子間的信息傳遞,提升信息搜索能力,并通過擁擠距離排序改善決策空間和目標空間的擁擠距離,最后,將端元束提取問題轉換到多模態多目標的問題上,提高算法效率的同時實現端元束提取的有效性。

1 線性光譜混合模型與優化

光譜解混[19-20]依賴于場景混合模型,大尺度的光譜混合可認為是線性混合問題[21]。本文主要圍繞線性混合模型進行混合像元問題的探究,并在此基礎上提出基于多模態多目標優化的端元束提取方法。

1.1 線性光譜混合模型

(1)

式中ej——端元

aij——端元ej在像元ri中的權重,即豐度

εj——誤差項

q——端元數量

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

aij>0

(2)

(3)

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

1.2 多模態多目標優化

多目標優化問題定義為

(4)

式中o——待優化的目標數

x——d維決策變量

gz(x)——不等式約束

hj(x)——等式約束

k——不等式約束數量

s——等式約束數量

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

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

圖1 多模態多目標問題示意圖

通常,多模態多目標問題滿足以下條件之一[23]:①至少有一個局部帕累托最優解。②至少有兩個等效全局帕累托最優解,它們對應PF上同一點。對滿足多模態多目標條件的高光譜影像進行端元束提取。

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

2 基于改進離散粒子群的端元束提取

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

2.1 粒子編碼

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

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

(5)

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

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

(6)

式中X——粒子像元對應的標號組成的序列

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

2.2 改進粒子速度與位置更新

PSO算法通過迭代尋優,在每次迭代中通過跟蹤個體最優xpbest和全局最優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)——粒子的更新速度

ω——慣性系數

c1、c2——學習因子

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

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

(9)

式中 ceil()——取整函數

xr(t+1)——區間Z內的隨機數

2.3 環形拓撲結構與擁擠距離

在多模態多目標優化中,粒子的選擇極為關鍵,本文采用MO_Ring_PSO_SCD[24]方法進行粒子信息傳遞和特殊擁擠距離計算。針對信息傳遞,設置個體最優xpbest的存檔PBA和鄰域最優xnbest的存檔NBA,在鄰域最優和子代最優選擇中采用環形拓撲更新方式,其更新規則為:對于初始化或者已經進行迭代更新的粒子群,當粒子j為1～P(P表示粒子數量)時,拓撲更新表現為粒子j僅與j-1和j+1粒子進行信息交互,當j=1時,拓撲更新表現為粒子j與P和第2個粒子進行信息交互,當j=P時,粒子j與第1個和P-1粒子進行信息交互,以此形成閉環拓撲結構信息傳遞。因此,將式(7)更新方法調整為

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

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

(10)

對于多模態多目標優化算法,需同時考慮決策空間和目標空間中個體的優劣問題。特殊擁擠距離Si可表示為

(11)

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

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

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

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

決策空間在空間分布上為離散問題,其數值范圍為區間Z內的正整數。中間粒子的決策距離表示為

(12)

其中

d=q

式中d——決策空間維度

δ——決策空間加權項

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

圖2 決策空間與目標空間擁擠距離計算方法示意圖

(13)

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

(14)

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

(15)

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

(16)

目標空間在空間分布上可認為是連續問題,由于目標函數值的區間不同,通常將目標空間值歸一化后進行距離計算。計算中間粒子的目標空間距離方法與決策空間類似,即

(17)

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

(18)

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

2.4 目標函數

(19)

2.5 算法流程

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

圖3 基于多模態多目標優化的端元束提取方法流程圖

3 試驗結果與分析

3.1 試驗數據

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

3.2 評價指標

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

(20)

式中b1、b2——光譜向量

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

3.3 試驗及討論

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

本文方法基于純像元假設進行端元束提取,像素標注結果如圖6所示。算法考慮到光譜變異性,從多模態多目標的角度實現更大范圍尋優。因地物重疊、分散造成的混合像元較多,MUUFL數據的環境相對復雜。通過算法尋優結果可發現,對于1#屋頂、2#草、3#樹、4#陰影端元采用MOPSOSCD算法能夠實現大區域內搜索,盡管圖6中5#瀝青分布較分散,MOPSOSCD算法的提取更傾向于尋優到5#瀝青的純像元信息。由于MOPSOSCD屬于多目標算法,個別端元存在錯分現象。總體來說算法在保證擴大搜索空間,避免陷入局部最優的效果較好。

圖6 MOPSOSCD算法端元提取結果

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

表1 MUUFL數據精度對比結果

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

圖7 MUUFL數據均方根誤差

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

表2 各端元的SAD值

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

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

4 結束語

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