空間和光譜信息保持的多光譜圖像融合算法

黃峰， 黃偉藍(lán)， 吳銜譽(yù)

(福州大學(xué)機(jī)械工程及自動(dòng)化學(xué)院， 福建 福州 350108)

0 引言

多光譜圖像是由多張光譜波段不同， 但空間分辨率相同的灰度圖像組成的數(shù)據(jù)立方體， 其在地物分類、 地質(zhì)調(diào)查以及環(huán)境質(zhì)量監(jiān)測(cè)等方面具有廣泛的應(yīng)用. 但是受多光譜成像器件硬件條件的限制， 多光譜成像傳感器往往只能獲取到低空間質(zhì)量的多光譜圖像， 而由全色圖像傳感器所獲取的全色圖像具有高空間分辨率. 為了獲得高空間質(zhì)量的多光譜圖像， 將多光譜和全色圖像進(jìn)行融合， 可實(shí)現(xiàn)優(yōu)勢(shì)互補(bǔ)， 在盡量保持多光譜信息的前提下， 提升融合圖像的空間質(zhì)量.

現(xiàn)有的多光譜圖像融合方法有成分替換法(component substitution， CS)、 多分辨率分析法(multi-resolution analysis， MRA)、 貝葉斯法和基于深度學(xué)習(xí)的方法. CS方法是指對(duì)多光譜圖像進(jìn)行投影， 以分離空間和光譜信息， 然后將包含空間信息的部分替換為全色圖像， 主要包括IHS[1]、 PCA[2]、 GS[3]、 AIHS[4]等， 該方法運(yùn)算耗時(shí)少， 但融合圖像易產(chǎn)生譜失真現(xiàn)象； MRA方法通過濾波將全色圖像的空間信息注入多光譜圖像， 主要包括Laplacian Pyramids[5]、 Wavelet[6]、 Discrete Wavelet Transform[7]等， 該方法運(yùn)算速度快且譜失真現(xiàn)象有所改善， 但可能會(huì)出現(xiàn)空間結(jié)構(gòu)失真的問題； 貝葉斯方法[8-10]用最大化后驗(yàn)概率的思路或變分的數(shù)學(xué)方法來(lái)建模， 能降低融合問題的病態(tài)程度， 但運(yùn)算耗時(shí)長(zhǎng)； 基于深度學(xué)習(xí)的方法[11-12]需要大量的數(shù)據(jù)訓(xùn)練模型， 限制了該方法在多光譜圖像融合方面的應(yīng)用.

上述融合方法存在如下兩個(gè)問題： 第一， 不能很好地描述多光譜與全色圖像之間的結(jié)構(gòu)相關(guān)性， 造成空間結(jié)構(gòu)失真； 第二， 不能很好地保持多光譜圖像的光譜信息， 造成光譜失真. 為了解決第一個(gè)問題， 本研究采用邊緣自適應(yīng)提取約束項(xiàng)和線性組合系數(shù)約束項(xiàng)來(lái)改善融合圖像的空間質(zhì)量； 為了解決第二個(gè)問題， 研究采用光譜信息保持約束項(xiàng)和波段比例關(guān)系保持約束項(xiàng)來(lái)保持多光譜圖像的光譜質(zhì)量. 聯(lián)立以上4個(gè)約束項(xiàng)， 即4個(gè)能量泛函， 建立總能量泛函， 使用梯度下降法進(jìn)行求解， 從而得到最優(yōu)值， 即最終的融合圖像.

1 約束項(xiàng)

本文方法主要基于以下4個(gè)約束項(xiàng)： ① 從全色圖像中自適應(yīng)提取邊緣細(xì)節(jié)信息并注入多光譜圖像中， 以建立邊緣自適應(yīng)提取約束項(xiàng)， 獲得第一個(gè)能量泛函； ② 根據(jù)全色圖像是融合圖像各波段的線性組合的假設(shè)， 估計(jì)多光譜圖像各波段線性組合系數(shù)， 以建立線性組合系數(shù)約束項(xiàng)， 獲得第二個(gè)能量泛函； ③ 根據(jù)低空間分辨率多光譜圖像是由融合圖像經(jīng)過模糊和下采樣處理后得到的假設(shè)， 建立光譜信息保持約束項(xiàng)， 獲得第三個(gè)能量泛函； ④ 根據(jù)融合前后多光譜圖像波段比例關(guān)系一致可減輕光譜失真的假設(shè)， 建立波段比例關(guān)系保持約束， 獲得第四個(gè)能量泛函. 將以上4個(gè)能量泛函合并到一個(gè)變分框架中， 從而獲得一個(gè)總的能量泛函.

1.1 邊緣自適應(yīng)提取約束

傳統(tǒng)的多光譜圖像融合算法通過將全色圖像的邊緣信息注入多光譜圖像, 從而獲得具有高頻細(xì)節(jié)信息的融合圖像. 但當(dāng)注入的邊緣信息較少時(shí)， 會(huì)造成融合圖像細(xì)節(jié)信息模糊， 當(dāng)注入的邊緣信息過多時(shí)， 則會(huì)造成融合圖像光譜失真. 本研究采用自適應(yīng)提取邊緣方法[4]來(lái)控制注入多光譜圖像的邊緣信息量， 即

Fi(x)=Mi(x)+h(x)(P(x)-I(x))

(1)

(2)

(3)

其中：x為像素位置，x∈Ω，Ω∈Rn為開放有界區(qū)間域；Mi為第i個(gè)波段的多光譜圖像；ai=1/m；I表示多光譜圖像的平均灰度；P為全色圖像； ▽P為全色圖像的梯度；ε為控制分母非零的小值；v為權(quán)衡參數(shù)， 控制全色圖像邊緣信息的注入量；h(x)為邊緣檢測(cè)函數(shù)， 用來(lái)控制注入的邊緣信息量(P-I)；Fi為第i個(gè)波段的融合圖像， 1≤i≤m.

由式(1)得到的能量泛函為

(4)

1.2 線性組合系數(shù)約束

全色圖像可以看作融合圖像各波段的線性組合[8]， 但線性組合系數(shù)如果相同， 則不能很好地表達(dá)多光譜圖像與全色圖像之間的結(jié)構(gòu)相關(guān)性. 為避免空間結(jié)構(gòu)失真， 本研究采用自適應(yīng)估計(jì)線性組合系數(shù)方法[4]計(jì)算線性組合系數(shù)αi， 即

(5)

(6)

其中：αi為第i個(gè)波段的線性組合系數(shù)， 1≤i≤m；P▽為降質(zhì)后的全色圖像.

式(5)的第一項(xiàng)保證產(chǎn)生的系數(shù)組合{α1, …,αi, …,αm}近似融合圖像各波段的線性組合； 第二項(xiàng)保證求得的系數(shù)αi是非負(fù)的，γ為拉格朗日乘子.使用梯度下降法求解式(5)， 即

(7)

(8)

(9)

由式(6)得到的能量泛函為

(10)

其中：Fi為第i個(gè)波段的融合圖像；P為全色圖像.

1.3 光譜信息保持約束

假設(shè)低空間分辨率的多光譜圖像是由高空間分辨率的多光譜圖像經(jīng)過模糊和下采樣處理后得到的[8]， 每個(gè)波段的模糊核存在一定的差異， 需要分波段進(jìn)行估計(jì)， 以保持光譜信息.

本研究采用Krishnan[13]提出的算法進(jìn)行模糊核的估計(jì)， 即

(11)

其中：ki為第i個(gè)波段的模糊核， 1≤i≤m；λ和φ控制著對(duì)應(yīng)項(xiàng)的權(quán)重大小； *表示卷積運(yùn)算.

表1 模糊核估計(jì)

在此假設(shè)下， 得到的能量泛函為

(12)

其中：Ds為下采樣矩陣.

1.4 波段比例關(guān)系保持約束

多光譜圖像波段之間具有一定的比例關(guān)系， 融合前后圖像波段比例關(guān)系保持恒定， 能減輕圖像的光譜失真程度.

采用波段比例關(guān)系保持約束來(lái)表征這種恒定的比例關(guān)系[14]， 即

(13)

由式(13)得到的能量泛函為

(14)

2 算法框架

算法框架如圖1所示.

圖1 算法框架Fig.1 Framework of the proposed method

2.1 總能量泛函

綜合4個(gè)能量泛函， 本研究使用總能量泛函對(duì)圖像進(jìn)行融合， 即

其中：γ,λ,μ,η≥0為正則化系數(shù)，γ和λ可調(diào)節(jié)空間信息約束的程度，μ和η可調(diào)節(jié)光譜信息約束的程度.

2.2 求解過程

為了使總能量泛函最小， 需要依次解決以下3個(gè)問題.

① 總能量泛函離散化. 數(shù)字圖像相當(dāng)于一個(gè)矩陣， 數(shù)字圖像的數(shù)據(jù)是離散數(shù)據(jù)， 設(shè)全色圖像P、 上采樣的多光譜圖像M和融合圖像F的分辨率為S×T， 定義一個(gè)S×T的矩陣H， 則矩陣上的每個(gè)點(diǎn)表示圖像上的一個(gè)像素x∈H.式(15)總能量泛函的離散形式為

(16)

② 總能量泛函極小化. 總能量泛函E(F)對(duì)應(yīng)每個(gè)波段Fi的歐拉-拉格朗日方程為

(17)

總能量泛函的極小化與相應(yīng)的歐拉-拉格朗日方程是等價(jià)的， 因此， 總能量泛函E(F)對(duì)應(yīng)的每個(gè)波段Fi(1≤i≤m)的歐拉-拉格朗日方程為

③ 求解融合圖像F. 通過梯度下降法求解式(18)， 使得

則有

表2 空間和光譜信息保持的多光譜圖像融合算法

續(xù)表2

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

使用QuickBird和Pavia University的圖像數(shù)據(jù)驗(yàn)證本文算法， 將本文算法與SFIM[15]、 MTF_GLP[16]、 MTF_GLP_HPM[17]、 PCA[2]、 GS[3]、 GSA[18]、 AIHS[4]和GFPCA[19]等多光譜圖像融合算法進(jìn)行比較, 并從定性和定量?jī)蓚€(gè)方面對(duì)各算法進(jìn)行綜合分析， 如圖2、 3所示.

(a) 參考圖像 (b) 全色圖像 (c) 經(jīng)過插值處理的多光譜圖像

(d) SFIM (e) MTF_GLP (f) MTF_GLP_HPM

(g) PCA (h) GS (i) GSA

(j) AIHS (k) GFPCA (l) 本文方法圖2 QuickBird融合結(jié)果Fig.2 Fusion results of QuickBird

分別將QuickBird和Pavia University中的多光譜圖像作為參考圖像， 大小為256 px×256 px， 如圖2(a)和圖3(a)所示； 對(duì)QuickBird和Pavia University中的多光譜圖像進(jìn)行光譜退化處理， 即對(duì)多光譜圖像中具有相同坐標(biāo)位置的像素值求平均， 從而得到全色圖像， 大小為256 px×256 px， 如圖2(b)和圖3(b)所示； 對(duì)QuickBird和Pavia University中的多光譜圖像分別進(jìn)行模糊和下采樣處理， 得到低空間分辨率的多光譜圖像， 大小為64 px×64 px， 將低空間分辨率的多光譜圖像插值放大至全色圖像大小， 大小為256 px×256 px， 如圖2(c)和圖3(c)所示. 本節(jié)利用得到的全色圖像和低空間分辨率的多光譜圖像進(jìn)行融合， 將融合圖像與參考圖像進(jìn)行對(duì)比. 經(jīng)實(shí)驗(yàn)驗(yàn)證， 當(dāng)γ=1、λ=0.5、μ=2.5、η=100、 Δt=10-8、ε′=5×10-3時(shí)， 融合圖像的空間分辨率可以得到提升、 光譜失真問題可以得到有效的減輕.

(a) 參考圖像 (b) 全色圖像 (c) 經(jīng)過插值處理的多光譜圖像

(d) SFIM (e) MTF_GLP (f) MTF_GLP_HPM

(g) PCA (h) GS (i) GSA

(j) AIHS (k) GFPCA (l) 本文方法

3.1 主觀定性評(píng)價(jià)

圖2(d)～(l)和圖3(d)～(l)分別表示SFIM、 MTF_ GLP、 MTF_ GLP_ HPM、 PCA、 GS、 GSA、 AIHS、 GFPCA及本文方法在QuickBird和Pavia University圖像數(shù)據(jù)中的融合結(jié)果. 如圖2(d)、 (f)所示， SFIM和MTF_ GLP_ HPM融合圖像的池塘區(qū)域會(huì)引入一些原圖中不存在的信息； 如圖2(e)所示， MTF_GLP融合圖像的池塘區(qū)域顏色發(fā)生輕微變化， 光譜信息略微失真； 如圖2(g)、 (h)所示， PCA和GS融合圖像的顏色和亮度發(fā)生了較大的改變， 這意味著PCA和GS的融合結(jié)果光譜失真嚴(yán)重； 如圖2(i)、 (j)所示， GSA和AIHS雖然引入了自適應(yīng)的概念， 光譜失真的程度有所緩解， 但GSA的融合結(jié)果產(chǎn)生了重影現(xiàn)象， AIHS的融合結(jié)果中融合圖像的池塘區(qū)域顏色發(fā)生了明顯的改變； 如圖2(k)所示， GFPCA的融合結(jié)果模糊程度嚴(yán)重， 空間分辨率極低； 如圖2(l)所示， 本文方法的融合結(jié)果在空間信息和光譜信息上與參考圖像相近， 產(chǎn)生了較好的融合效果. 如圖3(d)所示， SFIM融合圖像的亮度發(fā)生了變化； 如圖3(g)所示， PCA的融合結(jié)果變模糊； 如圖3(h)所示， GS的融合結(jié)果輕微模糊； 如圖3(k)所示， GFPCA的融合結(jié)果模糊程度嚴(yán)重， 空間分辨率極低； 如圖3(e)、 (f)、 (i)、 (j)、 (l)所示， 對(duì)比MTF_GLP、 MTF_GLP_HPM、 GSA、 AIHS和本文方法的融合結(jié)果， 本文方法的融合圖像顏色與參考圖像更為接近， 這說(shuō)明本文方法對(duì)光譜信息保存得更好， 同時(shí)， 本文方法的融合圖像與圖3(c)相比， 空間分辨率有了較大的提升.

3.2 客觀定量評(píng)價(jià)

為了對(duì)融合圖像的質(zhì)量進(jìn)行客觀評(píng)價(jià)， 采用如下4個(gè)客觀評(píng)價(jià)指標(biāo)來(lái)評(píng)價(jià)融合圖像的空間和光譜質(zhì)量.

1) 相關(guān)系數(shù)(correlation coefficients, CC)， 它表示圖像幾何失真的程度， 即

(21)

(22)

CC的值越接近于1， 融合圖像的幾何失真程度越小.

2) 光譜角映射(spectral angle mapper, SAM)， 它表示圖像光譜失真的程度， 即

(23)

(24)

3) 均方根誤差(root mean squared error, RMSE)， 它表示參考圖像與融合圖像對(duì)應(yīng)像素之間的誤差， 即

(25)

4) 全局相對(duì)光譜損失(relative dimensionless global error in synthesis, ERGAS)， 它表示融合圖像光譜失真程度， 即

(26)

以上4個(gè)指標(biāo)的理想最優(yōu)值分別寫在對(duì)應(yīng)指標(biāo)后的括號(hào)中， QuickBird和Pavia University融合結(jié)果的客觀定量對(duì)比如表3、 4所示.

表3 QuickBird融合結(jié)果

表4 Pavia University融合結(jié)果

表3、 4中最優(yōu)結(jié)果用黑體表示， 次優(yōu)結(jié)果加下劃線. 從表中可以看出， 本文方法在QuickBird圖像數(shù)據(jù)中各項(xiàng)評(píng)價(jià)指標(biāo)均是最優(yōu)值， 在Pavia University圖像數(shù)據(jù)中SAM指標(biāo)處于次優(yōu)的位置， 但其他指標(biāo)均是最優(yōu)值， 這表明本文方法能在盡量減少光譜扭曲程度的前提下， 使融合圖像的空間分辨率得到提升.

4 結(jié)語(yǔ)

本研究提出4個(gè)約束項(xiàng)， 即邊緣自適應(yīng)提取約束項(xiàng)、 線性組合系數(shù)約束項(xiàng)、 光譜信息保持約束項(xiàng)和波段比例關(guān)系保持約束項(xiàng)， 聯(lián)立4個(gè)約束項(xiàng)， 建立總能量泛函并求解最優(yōu)值， 以改善融合圖像的空間和光譜質(zhì)量. 相比其他方法， 本文方法的融合圖像視覺效果更為清晰. 綜上所述， 本文方法能在盡量保持多光譜圖像光譜信息的前提下， 提升融合圖像的空間質(zhì)量， 總體融合質(zhì)量?jī)?yōu)于其他方法. 下一步工作將側(cè)重于正則化系數(shù)的自適應(yīng)選擇， 以改善多光譜圖像融合的質(zhì)量.