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

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

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

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

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

l1范數(shù)定義為

N階張量X∈I1×I2×…×IN的模式n向量是一個以in為元素下標, 其他下標{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維實張量X均可分解為n模式積X=G×1U(1)×2U(2)…×NU(N), 其中U(n)是一個In×Jn半正交矩陣, 與奇異值矩陣不同, 核心張量G不取對角結構, 一般是一個滿張量, 表示各模式矩陣U(n)(n=1,2,…,N)之間的相互作用.

1 模型構建

1.1 問題描述

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

Y=X+N+S.

(1)

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

(2)

(3)

1.2 子空間表示

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

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

1.3 加權低秩張量正則化

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

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

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

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

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

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

(5)

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

(6)

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

(7)

2 優(yōu)化方法

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

2.1 求解加權低秩張量正則化項

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

(8)

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

(9)

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

U1=PQT,

(10)

(11)

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

2.2 求解簡化圖像Z

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

(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), 假設M奇異值分解的結果為M=UΛVT, 則

(18)

2.3 迭代正則化

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

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

(19)

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

k=k+β×n,

(20)

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

圖1 去噪過程Fig.1 Denoising process

2.4 復雜度分析

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

3 實驗結果

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

3.1 模擬數(shù)據(jù)實驗

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

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

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

續(xù)表1

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

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

表2 不同去噪算法的運行時間

3.2 真實數(shù)據(jù)實驗

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

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

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

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