非參數(shù)Bayesian字典學(xué)習(xí)的遙感影像融合方法

李 麗

(成都大學(xué) 建筑與土木工程學(xué)院，四川 成都 610106)

0 引 言

近年來，高光譜(hyperspectral，HS)影像因其所具有的豐富光譜信息，而被廣泛應(yīng)用于目標(biāo)的追蹤、識(shí)別與分類，以及災(zāi)害的監(jiān)測與管理等領(lǐng)域[1-3].依據(jù)HS的遙感成像技術(shù)，其HS影像雖具有高光譜分辨率與圖譜合一的特點(diǎn)，但由于成像光譜儀設(shè)計(jì)和制造方面的技術(shù)瓶頸，HS影像數(shù)據(jù)在信噪比、空間分辨率及掃描幅寬等性能方面受到了一定程度的制約[4].相關(guān)研究發(fā)現(xiàn)，若能獲取同一場景的高空間分辨率影像，則通過影像融合方法便可獲得具有高空間分辨率的HS影像[5-6].對(duì)此，科研人員對(duì)HS影像和多光譜(multispectral，MS)影像或全色(panchromatic,PAN)影像融合做了大量研究.例如，Aiazzi等[7]利用廣義拉普拉斯金字塔(generalized laplacian pyramid，GLP)的全色銳化方法對(duì)HS-MS影像進(jìn)行了融合；Berné等[8]提出了一種適用于中紅外波段的非負(fù)矩陣分解(nonnegative matrix factorization，NMF)的HS-MS影像融合方法；Yokoya等[9]提出了一種利用耦合非負(fù)矩陣分解(coupled nonnegative matrix factorization，CNMF)進(jìn)行光譜解混的HS-MS影像融合方法；Akhtar等[10]利用Bayesian稀疏表示的方法對(duì)端元和豐度進(jìn)行學(xué)習(xí)來實(shí)現(xiàn)HS-MS影像的融合.Wei等[11]利用在線字典學(xué)習(xí)(online dictionary learning，ODL)進(jìn)行字典學(xué)習(xí)，提出一種基于稀疏表示的HS-MS影像融合算法.雖然上述研究成果取得了很好的HS-MS融合效果，但其融合結(jié)果對(duì)模型的參數(shù)設(shè)置存在較強(qiáng)敏感性，且無法通過引入更多先驗(yàn)知識(shí)來進(jìn)一步減少空間信息融入不足和光譜畸變等問題.基于此，考慮到非參數(shù)Bayesian方法能夠充分利用圖像的先驗(yàn)信息，以及對(duì)未知參數(shù)利用概率估計(jì)代替點(diǎn)估計(jì)，且在估計(jì)過程中引入了不確定性，進(jìn)而避免了不同參數(shù)間估計(jì)誤差的擴(kuò)大等特點(diǎn)[12-13]，本研究提出一種非參數(shù)Bayesian字典學(xué)習(xí)的HS-MS影像融合方法，其特點(diǎn)是可充分利用先驗(yàn)信息，更好地融入空間信息和保持光譜信息，且具有更小的空間結(jié)構(gòu)誤差及光譜畸變程度.

1 方 法

本研究的目的在于通過對(duì)高光譜分辨率的HS影像YH與高空間分辨率的MS影像YM進(jìn)行融合，以此來恢復(fù)具有高空間分辨率的HS影像X.由此，基于HS-MS融合方法的融合過程定義為算法1，具體為：

算法1:非參數(shù)Bayesian字典學(xué)習(xí)的高光譜與多光譜影像融合輸入:YH,YM,NH,NM,B,S,R, mλ,K;%利用高光譜影像YH確定子空間轉(zhuǎn)換矩陣H^1 H^←PCAYH, mλ();%計(jì)算U的粗略估計(jì)值 U2 利用文獻(xiàn)[16]中所述方法計(jì)算U的粗略估計(jì)值 U?μ^U|YM;3 fori=1 to mλdo%利用算法2學(xué)習(xí)字典元素4 Di←Beta-Bernoulli( Ui);%利用文獻(xiàn)[11]的OMP算法學(xué)習(xí)稀疏系數(shù)5 Ai←OMP( Di, Ui,K);6 end%交替優(yōu)化7 fort=1,2,...to stopping rule 1 do%利用SALSA[15]及LS[11]對(duì)U、A進(jìn)行交替優(yōu)化更新8 U^(t)∈U|LU,A^(t-1)()≤LU^(t-1),A^(t-1)(){};9 A^(t)∈A|LU^(t),A()≤LU^(t),A^(t-1)(){};10 end%重建出高空間分辨率的高光譜影像11 Recover X^=H^U^;輸出:X^

1.1 非參數(shù)Bayesian字典及稀疏系數(shù)學(xué)習(xí)

算法2:基于非參數(shù)Bayesian的字典學(xué)習(xí)輸入: Ui,K,τ0,η0,c0,d0,e0,f0;1 從粗略估計(jì)值 Ui中選取N個(gè)觀測值;2 利用Gibbs抽樣對(duì)dk、sik、zik、πk、γ、γs等隱變量進(jìn)行初始化;3 fori=1,2,...to stopping rule 2 do 4 分別利用公式(6-8)對(duì)隱變量dk、sik、zik進(jìn)行優(yōu)化更新;5 利用公式(9)對(duì)πk進(jìn)行優(yōu)化更新;6 利用公式(10)對(duì)γ進(jìn)行優(yōu)化更新;7 利用公式(11)對(duì)γs進(jìn)行優(yōu)化更新;8 end9 計(jì)算dk的平均值;輸出: Di

算法2中，各隱變量對(duì)應(yīng)的抽樣公式分別為,

p(dk|-)∝

(6)

p(sik|-)∝

p(zik=1|-)∝

(8)

p(πk|-)∝

Beta(c0η0,c0(1-η0))

(9)

p(γ|-)∝

p(γS|-)∝

(11)

1.2 交替優(yōu)化的模型求解

在交替優(yōu)化過程中，算法1融合框架中的該優(yōu)化問題為一個(gè)標(biāo)準(zhǔn)的關(guān)于U和A的二次約束優(yōu)化問題.

由于矩陣間運(yùn)算維數(shù)較高，且相關(guān)算子難以對(duì)角化，因此，該類約束優(yōu)化問題通常很難進(jìn)行求解.為解決此問題，本研究選擇先固定稀疏系數(shù)A更新矩陣U，然后固定U更新A的交替優(yōu)化算法對(duì)U和A進(jìn)行聯(lián)合優(yōu)化更新.在該步驟中，本研究采用SALS算法和LS算法來分別解決目標(biāo)影像U及對(duì)應(yīng)稀疏系數(shù)A的約束優(yōu)化問題(見算法1中第9步).

2 結(jié)果與分析

為驗(yàn)證所提方法的有效性，本研究選取HS-MS融合方法中較為流行的GLP[7]、CNMF[9]和SRFM[11]3種算法與本文方法進(jìn)行仿真實(shí)驗(yàn)比較，并采用重構(gòu)信噪比(restored signal noise ratio，RSNR)、均方根誤差(root mean square error，RMSE)、光譜角填圖(spectral angle mapper，SAM)、通用圖像質(zhì)量評(píng)價(jià)指標(biāo)(universal imagequality index，UIQI)、相對(duì)全局融合誤差(erreur relative globale adimensionnelle de，ERGAS)和畸變程度(degree of distortion，DD)6個(gè)評(píng)價(jià)指標(biāo)[18]對(duì)融合結(jié)果進(jìn)行定量評(píng)價(jià).其中，RSNR、UIQI值越大，融合質(zhì)量越好；RMSE、SAM、DD及ERGAS值越小，融合質(zhì)量越高.

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

在實(shí)驗(yàn)中，參考影像采用由ROSIS(reflective optics system imaging spectrometer)傳感器于2013年獲取的位于意大利帕維亞地區(qū)的空間分辨率為1.3 m，像素大小為610×340的高空間分辨率高光譜影像.該原始影像共包含115個(gè)光譜波段,波長覆蓋范圍為0.430～0.838 μm[18].去除12個(gè)強(qiáng)噪聲波段及10個(gè)水汽吸收波段，本研究分別選取128×128×93(圖1中場景1)及256×256×93(圖1中場景2)大小的ROSIS影像數(shù)據(jù)進(jìn)行兩組仿真實(shí)驗(yàn).分別取紅、綠、藍(lán)3個(gè)波段，場景1及場景2所合成的假彩色影像如圖1所示(為便于顯示，以下場景1、場景2的相關(guān)影像均分別以原始像素的90%及50%比例顯示).

(a)場景1

在ROSIS影像數(shù)據(jù)的場景1及場景2的融合實(shí)驗(yàn)中，分別對(duì)ROSIS參考影像X的每個(gè)波段進(jìn)行5×5高斯空間濾波，然后在水平和垂直方向進(jìn)行4倍下采樣后得到HS影像YH.同時(shí)，使用IKONOS-like波譜響應(yīng)[18]對(duì)參考影像X進(jìn)行光譜降質(zhì)，得到含有4個(gè)光譜波段的MS影像YM.然后對(duì)YH和YM均添加加性零均值高斯噪聲，假定YH的前43個(gè)波段信噪比為35 dB，剩余50個(gè)波段信噪比為30dB，YM各波段影像信噪比均為30 dB.

2.2 參數(shù)設(shè)置

2.2.1 低維子空間學(xué)習(xí)

2.2.2 字典學(xué)習(xí)

2.3 結(jié)果分析

ROSIS數(shù)據(jù)的場景1與場景2在GLP[7]、CNMF[9]和SRFM[11]及本文方法下的融合結(jié)果如圖2、圖3所示.其中，3種對(duì)比算法的實(shí)驗(yàn)參數(shù)均參考原文獻(xiàn)取值.

(a)參考影像

由圖2、圖3可看出，本文方法的融合效果明顯優(yōu)于GLP[7]、CNMF[9]算法，略優(yōu)于SRFM[11]算法，尤其表現(xiàn)在圖2中白色地物區(qū)域處，以及圖3中建筑物邊緣紋理處，本文方法的細(xì)節(jié)恢復(fù)能力更強(qiáng).其主要原因在于利用非參數(shù)Bayesian方法進(jìn)行字典學(xué)習(xí)，能夠通過引入更多的先驗(yàn)信息來提高字典元素對(duì)目標(biāo)影像結(jié)構(gòu)和紋理信息的表達(dá)能力，從而進(jìn)一步提升融合結(jié)果的質(zhì)量.

(a)參考影像

為更加直觀地觀察不同融合算法下的融合結(jié)果在細(xì)節(jié)處的差異，本研究引入差值圖[10]，即參考影像與融合影像之間的絕對(duì)誤差圖來進(jìn)行主觀評(píng)價(jià)，點(diǎn)越亮代表差值越大，暗點(diǎn)則表示差異性小.場景1與場景2在不同算法下融合結(jié)果的絕對(duì)誤差如圖4、圖5所示.

(a)參考影像

由圖4、圖5可看出，本文方法的融合結(jié)果的絕對(duì)誤差均遠(yuǎn)小于GLP、CNMF算法，而略小于SRFM算法，且誤差分布的趨勢符合圖2、圖3的目視對(duì)比結(jié)果.此也進(jìn)一步證明了本文方法的有效性.

(a)參考影像

此外，除了上述主觀評(píng)價(jià)，表1、表2分別給出場景1與場景2在各算法下融合結(jié)果的客觀評(píng)價(jià)指標(biāo)數(shù)據(jù).

表1 ROSIS場景1不同算法下的定量評(píng)價(jià)指標(biāo)數(shù)據(jù)

表2 ROSIS場景2不同算法下的定量評(píng)價(jià)指標(biāo)數(shù)據(jù)

由表1、表2中PSNR、RMSE、SAM、UIQI、ERGAS和DD值可看出，本方法的融合效果仍為最優(yōu).此也進(jìn)一步證明了客觀評(píng)價(jià)結(jié)果與上述主觀評(píng)價(jià)結(jié)果的一致性.

3 結(jié) 論

針對(duì)高光譜影像難以兼具高空間分辨率及高光譜分辨率進(jìn)而在一定程度上限制了其應(yīng)用領(lǐng)域的問題，本研究在文獻(xiàn)[11]所提融合算法的基礎(chǔ)上提出一種基于非參數(shù)Bayesian字典學(xué)習(xí)的高光譜與多光譜影像融合方法.該方法在子空間內(nèi)進(jìn)行影像融合：首先，利用Beta Bernoulli過程對(duì)投影至低維子空間的高光譜與多光譜影像進(jìn)行字典學(xué)習(xí)，獲得字典元素的后驗(yàn)分布；然后，利用OMP算法進(jìn)行稀疏系數(shù)學(xué)習(xí)；最后，利用ADMM算法結(jié)合SALS算法與LS算法對(duì)目標(biāo)影像和對(duì)應(yīng)的稀疏系數(shù)進(jìn)行交替優(yōu)化更新，從而獲得具有更高空間分辨率的高光譜影像.與當(dāng)前流行的幾種HS-MS算法相比，本研究方法融合結(jié)果的空間結(jié)構(gòu)誤差和光譜畸變程度相對(duì)較小，兩組仿真實(shí)驗(yàn)中其融合結(jié)果均具有最優(yōu)的目視效果及定量評(píng)價(jià)指標(biāo)結(jié)果.其主要原因在于本研究所采用的非參數(shù)Bayesian方法能夠在字典學(xué)習(xí)的過程中引入更多先驗(yàn)信息，從而提高字典元素的表達(dá)精度以及目標(biāo)影像的估計(jì)精度.值得注意的是，本研究方法采用Gibbs進(jìn)行迭代抽樣，字典元素的后驗(yàn)分布依賴于大量的迭代抽樣，因此其計(jì)算效率相對(duì)較低，下一步的工作將探討引入變分貝葉斯等快速近似解求解方法以提高算法的計(jì)算效率.