基于子空間表示和加權(quán)低秩張量正則化的高光譜圖像混合噪聲去除方法

周 航, 蘇延池, 李占山, 花昀嶠

(1.吉林大學(xué) 軟件學(xué)院, 長春 130012；2.吉林大學(xué) 人工智能學(xué)院, 長春 130012；3.吉林大學(xué) 計算機(jī)科學(xué)與技術(shù)學(xué)院, 長春 130012；4.吉林大學(xué) 資產(chǎn)管理處, 長春 130012)

高光譜圖像是一種高維圖像, 可反映真實(shí)場景的空間信息和光譜信息.但由于圖像獲取過程中的不可控因素, 例如傳感器的靈敏度、光子效應(yīng)和校準(zhǔn)誤差等, 圖像中會不可避免地存在噪聲.這些噪聲將嚴(yán)重影響圖像信息的傳達(dá)和后續(xù)處理, 限制了其在視覺領(lǐng)域的進(jìn)一步應(yīng)用.圖像去噪作為視覺領(lǐng)域的預(yù)處理步驟, 旨在去除圖像中的噪聲成分.非局部自相似性和全局光譜相關(guān)性是高光譜圖像去噪領(lǐng)域內(nèi)的兩個重要先驗(yàn)知識.Othman等[1]首次將二者相結(jié)合.Zhang等[2]先將高光譜圖像分解為多個互相重疊的全波段三維圖像塊, 然后沿光譜維度將三維圖像塊展開成矩陣, 在低秩矩陣恢復(fù)框架下進(jìn)行去噪.為進(jìn)一步提升去噪性能, 文獻(xiàn)[3]提出了低秩矩陣近似模型, 以適應(yīng)不同波段噪聲強(qiáng)度的差異, 顯著提高了低秩矩陣恢復(fù)的去噪性能.此外, 不同矩陣秩的近似被集成到去噪模型中, 如γ-范數(shù)[4]和加權(quán)Schatten p范數(shù)[5]等.文獻(xiàn)[6]基于光譜相關(guān)性提出了低秩矩陣分解方法用于圖像去噪.上述方法將低秩張量展開成矩陣進(jìn)行分析, 破壞了圖像的內(nèi)在三維結(jié)構(gòu), 不能充分利用高光譜圖像的先驗(yàn)信息.在圖像高維空間上的塊匹配操作也導(dǎo)致這些算法的復(fù)雜度較高.

針對上述問題, 本文提出一種基于子空間表示和加權(quán)低秩張量正則化的方法(subspace representation and weighted low-rank tensor regularization, SWLRTR).首先, 通過噪聲圖像學(xué)習(xí)低秩正交基, 將高光譜圖像投影到低維子空間中, 得到維度縮減的圖像；然后, 引入加權(quán)低秩張量正則化項(xiàng)和l1范數(shù)分別約束隨機(jī)噪聲和稀疏噪聲；最后, 針對提出的模型設(shè)計相應(yīng)的求解算法.實(shí)驗(yàn)結(jié)果表明, 該方法在模擬噪聲和真實(shí)數(shù)據(jù)集上效果均較好.

N階張量X∈I1×I2×…×IN的F范數(shù)定義為

l1范數(shù)定義為

N階張量X∈I1×I2×…×IN的模式n向量是一個以in為元素下標(biāo), 其他下標(biāo){i1,i2,…,iN}in全部固定不變的In維向量, 用符號記作Xi1…in-1: in+1…iN.張量的矩陣化是指將一個N階張量重新組成一個矩陣形式的變換, 分為橫向展開和縱向展開, 將N階張量X∈I1×I2×…×IN的模式n向量沿水平方向依次排列, 稱為張量的模式n橫向展開, 用符號X(n)表示,X(n)∈In×(I1…In-1In+1…IN).將張量的模式n向量沿縱向依次排列, 稱為張量的模式n縱向展開, 記為X(n),X(n)∈(I1…In-1In+1…IN)×In.一個N階張量X∈I1×I2×…×IN與一個Jn×In矩陣U(n)的n模式矩陣積用符號X×nU(n)表示, 是一個I1…In-1×Jn×In+1…IN維張量, 其元素定義為

其中:j=1,2,…,Jn；ik=1,2,…,Ik；k=1,2,…,N.Tucker分解是一種高階奇異值分解, 是矩陣奇異值分解向高階的推廣, 每個I1×I2×…×IN維實(shí)張量X均可分解為n模式積X=G×1U(1)×2U(2)…×NU(N), 其中U(n)是一個In×Jn半正交矩陣, 與奇異值矩陣不同, 核心張量G不取對角結(jié)構(gòu), 一般是一個滿張量, 表示各模式矩陣U(n)(n=1,2,…,N)之間的相互作用.

1 模型構(gòu)建

1.1 問題描述

噪聲高光譜圖像Y∈n1×n2×n3(n1,n2分別為空間域上的長和寬,n3為頻譜數(shù)目)可被建模為干凈信號X與噪聲之和, 根據(jù)噪聲強(qiáng)度的分布特點(diǎn)可將噪聲分為隨機(jī)噪聲N和稀疏噪聲S.因此, 高光譜圖像混合噪聲去除模型可表示為

Y=X+N+S.

(1)

求解去噪問題的關(guān)鍵是基于高光譜圖像的先驗(yàn)信息建立合適的正則化項(xiàng).基于此, 高光譜圖像的去噪模型可表示為

(2)

(3)

1.2 子空間表示

高光譜圖像的每個光譜信號均可由少量光譜成分線性表示.假設(shè)高光譜圖像X存在于一個k維子空間Hk中, 干凈的高光譜圖像可表示為X=Z×3A, 其中A∈n3×k是包含簡化圖像Z∈n1×n2×k的基底的正交矩陣.因此, 高光譜圖像去噪模型(3)可表示為

He等[7]已經(jīng)證明, 矩陣A的正交性使得A中每列對應(yīng)的高光譜圖像特征之間具有較大的差異, 有效避免了Y中的噪聲信號在計算過程中出現(xiàn)在矩陣A中, 因此維度縮減后的圖像Z保留了與Y中相同的噪聲分布.此外, 將高光譜圖像分解為模-3張量矩陣積的形式保護(hù)了結(jié)構(gòu)信息的完整性, 并降低了算法的復(fù)雜度.

1.3 加權(quán)低秩張量正則化

本文用簡化圖像的先驗(yàn)信息通過組合相似的全光譜頻帶圖像塊得到一個具有明顯低秩性的張量.在實(shí)際應(yīng)用中, 為簡化計算一般在參考圖像塊的鄰域內(nèi)尋找.假設(shè)圖像塊的大小為p×p×k, 相似圖像塊的數(shù)目為q, 則構(gòu)造張量Zi的過程如下：

步驟1) 在參考圖像塊鄰域內(nèi)間隔特定步長選取圖像塊, 即搜索范圍;

步驟2) 計算參考圖像塊與鄰域內(nèi)圖像塊的歐氏距離, 衡量圖像塊之間的相似程度, 選取歐氏距離最小的前q個圖像塊組成相似圖像塊組;

步驟3) 將相似圖像塊組中的三維圖像塊沿第3模式縱向展開成p2×k大小的二維矩陣;

步驟4) 將步驟3)得到的q個二維矩陣堆疊為p2×q×k大小的三維張量Zi.

將上述過程描述為Zi=RiZ, 其中i是參考圖像塊左上角像素點(diǎn)坐標(biāo),Ri表示在簡化圖像Z上提取構(gòu)造Zi的一系列操作.不同的正則化項(xiàng)會影響去噪效果, 本文引入加權(quán)低秩張量正則化項(xiàng)約束Zi的低秩性.根據(jù)Zi估計其低秩近似, 即干凈張量Li可描述為

(5)

由于Zi的元素之間具有較強(qiáng)的相關(guān)性, 其核心張量的奇異值大多數(shù)接近于零.較大的奇異值對應(yīng)Zi主要的投影方向, 較小的奇異值一般大于干凈張量, 本質(zhì)上是噪聲導(dǎo)致的擾動.因此應(yīng)合理地對較大的奇異值施加較少的懲罰, 對較小的奇異值施加較多的懲罰.權(quán)重設(shè)置為

(6)

其中c>0是一個常數(shù),ε是一個保證除數(shù)不為零的較小常數(shù).基于以上分析, 本文提出的基于子空間迭代的加權(quán)低秩張量正則化高光譜圖像去噪模型形式如下：

(7)

2 優(yōu)化方法

本文引入一種基于迭代最小化策略的算法求解提出的模型.模型(7)的求解過程包括兩步： 求解加權(quán)低秩張量正則化項(xiàng)和簡化圖像.

2.1 求解加權(quán)低秩張量正則化項(xiàng)

本文從低秩張量Zi中估計張量Li, 忽略模型中與Li無關(guān)的變量, 得到如下子問題：

(8)

將Li替換為其Tucker分解的形式, 則式(8)可轉(zhuǎn)化為

(9)

其中°表示基于元素的乘積,j表示一個三階張量的模式索引.分別求解U1,U2,U3, 令W1=Zi(1)(U3?U2), 其中?表示Kronecker積, 則

U1=PQT,

(10)

(11)

得到4個變量后, 用模-3張量矩陣積恢復(fù)最優(yōu)Li.

2.2 求解簡化圖像Z

固定Li, 則式(7)可轉(zhuǎn)化為下列最小化問題：

(12)

式(12)的增廣Lagrange函數(shù)為

其中ξ{Ik}(·)是指示函數(shù).采用交替最小化方法迭代優(yōu)化S,Z,A, 求解式(13), 步驟如下.

1) 更新S： 固定Z和A, 則式(13)可被改寫為

(14)

其閉形式的解為Sλ2(Y-Z×3A), 其中Sτ(x)=sgn(x)max{|x|-τ,0}是軟閾值操作.

2) 更新Z： 固定S和A, 則式(13)可被改寫為

(15)

式(15)是一個二次優(yōu)化問題, 具有封閉解, 表示為

(16)

3) 更新A： 固定S和Z, 則式(13)可被改寫為

(17)

令M=(Y-S)(3)Z(3), 假設(shè)M奇異值分解的結(jié)果為M=UΛVT, 則

(18)

2.3 迭代正則化

迭代正則化為每次迭代的輸入加入一定比例的噪聲信號作為對結(jié)果的擾動.該策略已被廣泛應(yīng)用以提升去噪效果.在下一次迭代中更新噪聲圖像的公式為

Yn+1=αXn+(1-α)Y,

(19)

其中α是平衡下次迭代輸入的參數(shù),n是迭代次數(shù).子空間大小k的選擇與不同的數(shù)據(jù)集和噪聲強(qiáng)度有關(guān), 本文根據(jù)HySime[8]初始化k.隨著每次迭代噪聲強(qiáng)度越來越弱, 得到的去噪圖像越來越接近真實(shí)圖像,k應(yīng)逐漸增大使維度縮減后的圖像Z保留更多的干凈圖像信息.因此,k隨迭代次數(shù)變化的公式為

k=k+β×n,

(20)

其中β是一個決定步長的常數(shù).從而每次迭代后An+1可從輸入圖像中提取更多信息.去噪過程對應(yīng)的流程如圖1所示.

圖1 去噪過程Fig.1 Denoising process

2.4 復(fù)雜度分析

其中N1是估計S,Z,A時達(dá)到收斂所需的迭代次數(shù),N是去噪過程整體的迭代次數(shù).

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

下面用仿真實(shí)驗(yàn)和真實(shí)數(shù)據(jù)集上的實(shí)驗(yàn)說明本文方法的有效性.將本文方法與LRTDTV[9],NMoG[10],L1HyMixDe[11],SDeCNN[12],3DLogTNN[13],SNLRSF[14]高光譜圖像去噪方法進(jìn)行比較, 對比方法的參數(shù)設(shè)置參照源代碼和原文獻(xiàn)中的要求.所有測試高光譜圖像在去噪前都?xì)w一化到[0,1]波段, 所有算法在16 GB RAM的Intel(R) Core(TM) i7-8700 CPU 3.20 GHz環(huán)境下, 用MATLAB R2018a實(shí)現(xiàn).

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

在數(shù)據(jù)集Washington DC(WDC),Pavia University(PaviaU)和Pavia Center(PaviaC)上進(jìn)行仿真實(shí)驗(yàn).真實(shí)場景下的高光譜圖像通常被不同類型的噪聲污染, 包括高斯噪聲、脈沖噪聲和條紋等.為模擬真實(shí)高光譜圖像的噪聲分布, 將不同類型的噪聲以不同的形式組合加入到高光譜圖像中： 第一種是具有相同分布的高斯噪聲加入到3個模擬數(shù)據(jù)集的每個頻帶中, 其均值為0, 方差為[0.1,0.2]；第二種是具有不同強(qiáng)度的高斯噪聲加入到3個模擬數(shù)據(jù)集的每個頻帶中, 其均值為0, 噪聲方差σ為[0.1,0.2]；第三種將與第二種情況相同的高斯噪聲加入到3個高光譜數(shù)據(jù)集的每個頻帶中, 此外, 隨機(jī)選取20個頻帶加入20%的脈沖噪聲；第四種是在原始圖像中加入與第三種情況相同的高斯噪聲和脈沖噪聲.本文從含有脈沖噪聲的頻帶和其他頻帶中分別選取10個頻帶加入寬度為1～3的死線.

本文用4個評價指標(biāo)度量去噪效果： 平均峰值信噪比(mean peak signal-to-noise ratio, MPSNR)、平均結(jié)構(gòu)相似度(mean structural similarity, MSSIM)、相對無量綱全局誤差(relative dimensionless global error synthesis, ERGAS)和平均光譜角(mean spectral angle, MSA).表1列出了不同算法在上述3個數(shù)據(jù)集上的去噪效果, 包括在4種模擬噪聲情況下不同評價指標(biāo)的值.由表1可見, 本文算法的去噪結(jié)果在大多數(shù)情況下都最好.

表1 不同算法在3個數(shù)據(jù)集上的定量實(shí)驗(yàn)結(jié)果

續(xù)表1

圖2為不同算法在數(shù)據(jù)集PaviaC頻帶64上的去噪結(jié)果.由圖2可見, 本文算法更好地保留了圖像的細(xì)節(jié)信息并去除了噪聲.表2列出了不同去噪算法的運(yùn)行時間.由表2可見, 本文算法在計算時間方面具有優(yōu)勢, 包括不同算法在3個數(shù)據(jù)集上的平均運(yùn)行時間.

圖2 不同算法對數(shù)據(jù)集PaviaC在加入第4種噪聲后的去噪結(jié)果Fig.2 Denoising results of different algorithms on PaviaC dataset after adding the fourth noise

表2 不同去噪算法的運(yùn)行時間

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

由于真實(shí)數(shù)據(jù)集KSC沒有干凈的圖像作為參考, 因此本文采用帶有徑向基函數(shù)(radial basis function, RBF)核的支持向量機(jī)得到的分類結(jié)果定量評估本文算法在真實(shí)數(shù)據(jù)集上的去噪性能.表3列出了不同算法每個類別的分類準(zhǔn)確率, 并給出了3個評價指標(biāo), 分別為總體準(zhǔn)確率(overall accuracy, OA)、平均準(zhǔn)確率(average accuracy, AA)、Kappa統(tǒng)計量.由表3可見, 本文算法在3個評價指標(biāo)上都取得了優(yōu)于其他算法的分類結(jié)果, 表明本文算法在真實(shí)數(shù)據(jù)集上同樣具有較好的去噪效果.圖3為不同算法在數(shù)據(jù)集KSC頻帶89上的去噪效果.由圖3可見, 本文算法去噪后的圖像具有最好的視覺效果.

表3 基于不同去噪算法數(shù)據(jù)集KSC的分類結(jié)果

圖3 不同算法在數(shù)據(jù)集KSC上的去噪結(jié)果 Table 3 Denoising results of different algorithms on KSC dataset

綜上所述, 本文提出了一種基于子空間表示的加權(quán)低秩張量正則化的高光譜圖像去噪方法, 針對高光譜圖像中的高斯噪聲和稀疏噪聲進(jìn)行建模.相比于其他去噪算法, 本文模型較好地保留了張量的固有結(jié)構(gòu), 提高了高光譜圖像先驗(yàn)信息的利用效率, 并且計算復(fù)雜度低于多數(shù)去噪算法.實(shí)驗(yàn)結(jié)果表明, 本文方法在模擬噪聲圖像和真實(shí)高光譜圖像上都取得了較好的實(shí)驗(yàn)結(jié)果.