楊 垚 , 黃 聰, 王華軍

(成都理工大學(xué) a.地球物理學(xué)院,b.地球勘探與信息技術(shù)教育部重點(diǎn)實(shí)驗(yàn)室,成都 610059)

0 引言

高光譜遙感影像(HSI)是當(dāng)下遙感領(lǐng)域的研究熱點(diǎn)，其影像包含大量空間信息和光譜信息，眾多研究領(lǐng)域已廣泛使用HSI提取可用信息。然而，HSI在成像、傳輸?shù)冗^程中無可避免的受到噪聲的干擾，復(fù)雜的噪聲類型(高斯噪聲、脈沖噪聲、死線等),降低了圖像的分辨率并給科學(xué)研究增加了難度。應(yīng)對(duì)復(fù)雜噪聲的干擾，盡量提高信噪比且最大程度保留有用信息數(shù)據(jù)的去噪流程是必需的。各領(lǐng)域由于研究目的和方法的不同，影像的去噪速度和去噪效果等眾多因素各有所需。

為了更加適用于各領(lǐng)域，眾多HSI去噪算法相繼被提出，HSI數(shù)據(jù)有多個(gè)頻帶，每個(gè)頻帶可看作一張灰度圖像，2-D圖像去噪的方式是最直接的(如BM3D、WNNM)，此類方法利用了頻帶信息，忽略了頻帶之間的強(qiáng)相關(guān)性，實(shí)際中不僅不能有效地恢復(fù)圖像，而且還會(huì)有偽像和失真的現(xiàn)象[1-2]?？紤]HSI多頻帶相關(guān)特性，眾多基于譜間相關(guān)性的3-D去噪方法被提出。3-D去噪方法分為兩類，①基于變換域，該方法一般基于主成分分析，離散變換或小波變換(如基于主成分分析、基于離散余弦變換(DCT)、基于小波變換的濾波[3-5])；②基于空間域的方法，該方法更多地利用了空間和頻譜信息。稀疏先驗(yàn)?zāi)Ｐ褪褂米值鋵W(xué)習(xí)來有效利用HSI的RAC(the redundancy and correlation，冗余與相關(guān)性)特點(diǎn)，因HSI干凈數(shù)據(jù)可以由少量原子線性表示，而噪聲數(shù)據(jù)不能(如將高度相關(guān)的頻段聚集到一起充分利用全局和局部的RAC，再用頻譜和空間聯(lián)合稀疏表示和聯(lián)合字典學(xué)習(xí)的去噪辦法)[6]。全變分(TV)先驗(yàn)?zāi)Ｐ瓦\(yùn)用像素之間的平滑特征，方法在去噪和保護(hù)邊緣的統(tǒng)一表現(xiàn)突出，但是由于HSI噪聲類型復(fù)雜，一般與其他模型一起使用(如提出全變分正則化和低秩約束模型結(jié)合使用，低秩和稀疏矩陣的分解為TV正則化提供噪聲分量，TV正則化后獲得了增強(qiáng)的干凈圖像[7])。低秩先驗(yàn)?zāi)Ｐ蛷母呔S數(shù)據(jù)中發(fā)現(xiàn)和利用低維結(jié)構(gòu)來去噪，是HSI去噪最有力的約束模型?；诟蓛鬑SI低秩先驗(yàn)特性和非高斯噪聲稀疏性，建立了低秩矩陣恢復(fù)模型(low-rank matrix recovery，LRMR)，使用GoDec和增廣拉格朗日乘子法(ALM)求解提出的恢復(fù)模型，最終得到低秩的矩陣(干凈的HSI)[8-9]；LRMR模型沒有考慮不同頻帶噪聲強(qiáng)度不同的事實(shí)，進(jìn)而提出一種噪聲調(diào)整的迭代LRMA改進(jìn)方法(NAILRMA)，引入隨機(jī)奇異值分解來求解提出的迭代模型，通過每次迭代的反饋數(shù)據(jù)不斷更新噪聲強(qiáng)度來提升去噪能力[10]；構(gòu)建了l0梯度正則化低秩張量(LETL0)模型，使用ALM求解模型，用此方法降噪后HSI圖像邊緣大部分保留[11]。

由于頻譜信息冗余的充分使用，低秩先驗(yàn)?zāi)Ｐ驮贖SI去噪上大放光彩，然而在HSI空間上的相關(guān)性卻沒有充分利用。這里考慮了低秩模型的長(zhǎng)短處，提出一種基于子空間低秩和稀疏矩陣分解的迭代算法SLRSF(subspace low-rank matrix and sparse matrix factorization)，將HSI每個(gè)頻帶矢量化處理并將三維張量重塑為二維矩陣，再將二維矩陣映射到低維頻譜子空間，使用低秩和稀疏張量正則化的方法迭代處理特征圖像得到干凈特征圖像，再還原得到干凈圖像。該方法充分利用了RAC，算法效率高，去噪能力優(yōu)秀。

1 方法

高光譜影像噪聲類型眾多，常見且對(duì)數(shù)據(jù)干擾較大的噪聲有高斯噪聲、脈沖噪聲、壞點(diǎn)或條紋、死線等?？紤]這些噪聲數(shù)據(jù)特點(diǎn)(稀疏和低秩)可將HSI分解為：

Y=X+S+N

(1)

其中:Y為觀測(cè)數(shù)據(jù)；X為干凈的數(shù)據(jù)；S為稀疏噪聲(包含脈沖噪聲、條紋、死線)；N為高斯噪聲；Y、X、S、N∈Rh×p；m、n分別是HSI圖像空間的長(zhǎng)和寬；p是光譜波段數(shù)；h=m*n。

2.1 頻譜低秩分解

HSI光譜波段數(shù)較多且高度相關(guān)，為此使用一個(gè)k維子空間E來線性表示X。

X=EZ

s.t.ETE=Ik

(2)

其中:k?p；Ik∈Rk×k是單位矩陣；X∈Rh×p；E∈Rh×k，E由HySime[12]算法計(jì)算得出，Z被稱為RCI(representation coefficients image)，它的行稱為本征圖像，Z∈Rk×p。進(jìn)而HSI退化模型表示為式(3)。

Y=EZ+S+N

(3)

1.2 迭代模型

假設(shè)N是滿足高斯分布且均值為0，HSI的去噪問題表示為：

(4)

(5)

式(4)中第一項(xiàng)控制數(shù)據(jù)保真度，第二項(xiàng)控制約束稀疏噪聲，第三項(xiàng)控制干凈圖像的低秩特性。

式(4)的低秩正則化和稀疏正則化的拉格朗日模型為式(6)。

rank(X)≤r

(6)

由式(2)，RCI的去噪模型可以表示為式(7)。

rank(Z)≤r

(7)

1.3 優(yōu)化及實(shí)現(xiàn)

對(duì)于式(7)，我們使用增廣拉格朗日乘子法(augmented lagrange method)來解決它。

λ‖S2‖1+<Λ,Y2-Z-S2>+

(8)

Y2=ETY,S2=ETS

(9)

式中：Λ為拉格朗日乘數(shù)，式(8)的解法為對(duì)其中一個(gè)變量迭代時(shí)固定其余變量。

1.3.1Z的迭代

包含Z的子問題為：

(10)

參考文獻(xiàn)[13]中第二節(jié)給出了求解式(10)的具體方法。

(11)

1.3.2S2的迭代

包含S2的子問題為：

(12)

運(yùn)用下面的軟閾值運(yùn)算來解式(12)，其中T>0。

(13)

(14)

1.3.2Λ的迭代

Λ(k+1)=Λ(k)+μ(Y2-Z(k+1)-S(k+1))

(15)

基于SLRSF的低秩先驗(yàn)?zāi)Ｐ退惴鞒倘缦隆?/p>

輸入：加噪HSI數(shù)據(jù)Y∈Rm×n×p，參數(shù)lamda，迭代停止條件ε

輸出：干凈HSI數(shù)據(jù)X

1：Y=reshape(Y,m*n,p)’

2： 使用HySime從Y中計(jì)算子空間E和subdim

3：E=E[:,1:subdim]

4：Y2=E’ *Y;

5：Y2=reshape(Y2’,m,n,subdim)’

6：Y2_clean=RCI_LR(Y2,lamda)

6.1：Z=0，S2=0，∧=0，μ=0.01，μmax=106，ρ=1.5

6.2： 用式(11)、式(14)、和式(15)分別更新Z、S2、∧

6.3：μ=min(ρμ,μmax)

6.4： 判斷，若不滿足不等式則繼續(xù)從5.2開始執(zhí)行，若滿足則迭代完成，若超過一定迭代次數(shù)則自動(dòng)退出。

7：Y2_clean=reshape(Y2_clean’,m,n,subdim)’

8：X=E*Y2_clean

用于混合噪聲的高光譜圖像降噪。

2 算法測(cè)試與結(jié)果分析

為驗(yàn)證方法和代碼的可行性與效果，分別對(duì)Indian pines、KSC、Pavia、Pavia80、Salinas、Salinas A、Cuprite數(shù)據(jù)集加噪。加噪類型分為：①標(biāo)準(zhǔn)差相同的高斯噪聲；②標(biāo)準(zhǔn)差不同的高斯噪聲；③標(biāo)準(zhǔn)差不同的高斯噪聲，并隨機(jī)加入脈沖噪聲；④標(biāo)準(zhǔn)差不同的高斯噪聲，隨機(jī)加脈沖噪聲和死線。

這里選取了ALM RPCA、LRMR、BM4D方法對(duì)比本文提出的方法。量化指標(biāo)采用平均峰值信噪比(MPSNR)、平均結(jié)構(gòu)相似度(MSSIM)和光譜差異指標(biāo)(ERGAS)來測(cè)評(píng)這些方法。視覺上圖像細(xì)節(jié)保存尺度和噪聲消除能力大小也是一重要指標(biāo)。

2.1 Indian pines

Indian pines數(shù)據(jù)集是由AVIRIS傳感器在印第安納州西北部采集的數(shù)據(jù)，波長(zhǎng)范圍為0.4*10-6m到2.5*10-6m，224個(gè)光譜帶，每個(gè)光譜帶145*145像素。經(jīng)過去除吸水條帶的數(shù)據(jù)處理后，還剩220個(gè)光譜。

對(duì)于混合噪聲的去噪后圖像的所有波段，BM4D的去噪效果明顯不佳，高斯噪聲與死線均無法去除，LRMR去死線的效果不太好，其余方法均能去除死線、高斯噪聲和椒鹽噪聲的影響。圖1(b)、圖1 (d)、圖1 (e)和圖2(b)、圖2 (e)的右下圖區(qū)域恢復(fù)不是很清晰，提出的SLRSF基本能保留原圖像的細(xì)節(jié)信息且去噪效果良好。

圖1 Indian pines數(shù)據(jù)集第129波段死線與高斯噪聲情況實(shí)驗(yàn)結(jié)果對(duì)比

圖2 Indian pines數(shù)據(jù)集第29波段椒鹽噪聲與高斯噪聲情況實(shí)驗(yàn)結(jié)果對(duì)比

2.2 KSC

Kennedy Space Center(KSC)是NASA AVIRIS儀器在弗羅里達(dá)州肯尼迪航天中心采集的224個(gè)10 nm帶寬的高光譜數(shù)據(jù)，中心波長(zhǎng)范圍在400 nm到2 500 nm之間。去除吸水條帶和低SNR頻段后，剩下176個(gè)頻段可用。

對(duì)于Indian pines數(shù)據(jù)集(表1)，只有高斯噪聲的情況下SLRSF和BM4D的量化指標(biāo)都表現(xiàn)不錯(cuò)，其中SLRSF綜合表現(xiàn)更為優(yōu)秀。當(dāng)加入脈沖噪聲或死線混合噪聲時(shí)SLRSF各項(xiàng)量化指標(biāo)都比其余三種方法好，且花費(fèi)的時(shí)間遠(yuǎn)遠(yuǎn)低于它們。

表1 基于Indian pines數(shù)據(jù)集的算法比較

對(duì)于KSC數(shù)據(jù)集(表2)(數(shù)據(jù)量是Indian pines的12倍)，當(dāng)數(shù)據(jù)集增大時(shí)，SLRSF算法運(yùn)行所需內(nèi)存約為其余方法的0.1倍，且運(yùn)行時(shí)間遠(yuǎn)小于其余方法。經(jīng)過對(duì)各頻段的圖像對(duì)比后發(fā)現(xiàn)，BM4D在純高斯噪聲情況下表現(xiàn)優(yōu)秀，混合噪聲的去噪效果不好；ALM RPCA在各種噪聲下平均結(jié)構(gòu)相似度很高；SLRSF在4種噪聲情況下量化指標(biāo)均在中等偏上水平，且在噪聲類型越復(fù)雜的數(shù)據(jù)集下表現(xiàn)優(yōu)于其余方法。

表2 基于KSC數(shù)據(jù)集的算法比較

當(dāng)將算法運(yùn)用于數(shù)據(jù)集Salinas(AVIRIS成像光譜儀對(duì)美國(guó)加利福尼亞州Salinas山谷所成的像)和Pavia University(德國(guó)機(jī)載反射光學(xué)光譜成像儀對(duì)意大利的帕維亞城所成的像)時(shí)，發(fā)現(xiàn)若一個(gè)原始波段被噪聲干擾特別嚴(yán)重(基本無圖像顯示)而相鄰波段有數(shù)據(jù)時(shí)，這個(gè)波段在去噪后會(huì)有不同程度的假圖像生成，可能是過度利用了光譜的相關(guān)性?？偟膩碚f，對(duì)于不同數(shù)據(jù)集SLRSF均能快速有效地去除人工合成噪聲，各光譜去噪視覺效果表現(xiàn)不錯(cuò)，各量化指標(biāo)均表現(xiàn)突出。

3 結(jié)論

筆者基于低秩和稀疏矩陣分解方法，針對(duì)該類方法在空間相關(guān)性信息上利用的不足，提出一種基于子空間低維投影特征圖像低秩和稀疏矩陣分解的迭代方法。通過噪聲迭代的方式逐步消除高斯噪聲和非高斯噪聲的影響，得到干凈圖像。測(cè)試模擬數(shù)據(jù)集比較結(jié)果后得出，該方法不僅能快速有效去除高斯噪聲、椒鹽噪聲和死線的影響，恢復(fù)出高質(zhì)量的干凈圖像，并且在定量評(píng)價(jià)指標(biāo)的表現(xiàn)上也非常優(yōu)秀，最突出的是算法空間復(fù)雜度遠(yuǎn)低于其余方法，用于HSI圖像即時(shí)處理和批量處理等的選擇上有很大的優(yōu)先性。