采用改進RPCA的遙感影像去云算法

石曉旭，夏克文，王寶珠，常 虹，武盼盼

(河北工業大學 電子與信息工程學院，天津 300401)

0 引 言

如何在遙感影像中去除云層的遮擋，重構遮擋區地貌細節信息，是遙感應用領域一個待處理的重要問題。目前，在遙感領域，不少學者根據云層的厚度提出了多種去云算法：針對薄云去除，一般將遙感影像轉換到頻域處理，常見有同態濾波法和小波分解等方法，但會丟失某些高頻細節信息，造成信息量損失，而且無法去除云區和無云區之間的突變云域和小片云。對于厚云覆蓋區域，幾乎難以獲取任何地貌信息，只能利用多時相、多傳感器遙感影像進行融合處理，但需要復雜的影像配準、校正技術，容易造成去云后影像出現畸變、灰度突變、紋理不連續等多種問題。由此可見，對成分復雜云區的聯合去云技術的研究具有十分重要的意義。

近幾年，采用壓縮感知[1]的稀疏模型發展突飛猛進，引起眾多學者的關注，在采樣過程中利用觀測矩陣稀疏化原有數據信號，去除冗余信息，使數據更加魯棒，去除了野點和稀疏噪聲對數據的影響。低秩矩陣恢復(low-rank matrix recovery，LRMR)[2]是壓縮感知稀疏模型的二維矩陣形式擴展，將觀測矩陣表述為低秩矩陣與稀疏矩陣之和，再通過將非凸優化問題轉化為凸優化問題求解低秩矩陣。目前，LRMR主要由魯棒主成分分析(RPCA)矩陣補全(MC)和低秩表示(LRR)等3類模型組成。

本文針對同一區域、不同時相遙感影中的稀疏云層遮擋，提出一種采用改進的RPCA去云算法，即對原有RPCA算法中l0范數采取一種新的分式函數逼近方式，同時引入加權核范數對奇異值閾值進行自適應的調節，滿足保留數據主要成分的實際需要，使改進的RPCA算法在處理遙感影像去云問題時更為有效，實現復雜的稀疏云區的聯合去云處理。

1 問題的提出

在魯棒主成分分析(robust PCA，RPCA)模型中，觀測矩陣D∈Rm×n往往具備低秩或近似低秩結構，但是同時存在較大的稀疏噪聲干擾，如果將矩陣D表示為一個低秩矩陣和一個稀疏矩陣之和，則從中恢復出原有矩陣的低秩結構的概率較大，能夠去除毀損矩陣結構的稀疏噪聲。恢復低秩矩陣A可描述為如下雙目標優化問題

(1)

通過引入折中因子λ(>0)，將低秩矩陣恢復的雙目標優化問題式(1)轉化為式(2)最小化問題

(2)

本文研究的多時相的遙感影像具有明顯的時間相關性，即圖像地物特征在一段時間內比較穩定，圖像序列之間具有極大的相似性，所以由多幅地物像素組成的觀測矩陣具有低秩特性；而處于動態變換狀態中的云這一要素是非常不穩定的目標，不具備低秩性，同時，在圖像上的分布范圍比較小，在由多幅遙感影像組成的觀測矩陣中分布更加稀疏，具有顯著的稀疏特性。所以，在去除遙感影像中稀疏云層時，可以利用觀測矩陣的低秩和稀疏兩種特性，構建去云模型，實現低秩矩陣恢復。

針對遙感影像去云的RPCA算法模型可描述為：所觀測的多時相遙感影像包含n個圖像序列，每個圖像像素值均可轉化為一個m維列向量，那么便可構造m×n維觀測矩陣D來表示輸入遙感影像序列，待恢復的低秩矩陣A表示具有極大相似性的地物信息，稀疏矩陣E表示分布范圍較小的云層部分，即D=A+E，于是形成式(2)所示的凸優化問題。通過有效的算法求解，恢復出A和E，就可以實現成分復雜的稀疏云區的聯合去云。

2 改進的RPCA算法研究

本節主要針對式(2)中的稀疏矩陣E和低秩矩陣A進行改進，其中，為增加稀疏矩陣E恢復成功率，構建基于分式函數的RPCA模型；為增加低秩矩陣A的低秩性，采用加權核范數最小化算法(WNNM)中自適應閾值對奇異值進行相應地收縮。

2.1 基于分式函數的稀疏矩陣優化算法

2.1.1 模型改進

在式(2)中，由于低秩矩陣恢復優化模型的秩函數和l0范數都是矩陣空間上的非連續函數，因此難以得到式(2)的最優解，所以需要對其進行凸松弛處理。

(3)

凸松弛為

(4)

而針對l0范數優化問題，現有的算法[3]主要有貪婪算法、迭代硬閾值算法、基追蹤算法、迭代重賦權最小二乘極小化等算法，但是這些算法誤差比較大，很難精確恢復稀疏解，并且在處理本文遙感影像序列這類大規模矩陣時，計算的復雜度會變大，計算速度也會相應的變慢。為了更好挖掘稀疏矩陣的稀疏度，增強恢復稀疏信號的能力，提高有噪音干擾的情況下算法的穩定性，本文考慮利用分式函數ρb(|t|)去替代不連續l0范數，分式函數圖像如圖1所示

(5)

式中：b∈(0,+∞)。

圖1 分式函數圖像

(6)

2.1.2 模型求解過程

為求解式(6)的優化問題，本文對式(6)構造相應地增廣拉格朗日函數，如式(7)所示

(7)

其中，Y∈Rm×n表示拉格朗日乘子，懲罰參數μ>0，·是標準內積，當Y=Yk，μ=μk，使用交替投影法[4]求解塊優化問題交替迭代更新矩陣A和E，直到滿足最優解收斂條件為止。

(8)

其中，奇異值算子Dε(Q)=USε(∑)VT，U∑VΤ為矩陣Q的奇異值分解，軟閾值Sε(Q)的第(i,j)個元素為max(|qij|-ε,0)sgn(qij)，ε為大于0的常數。

(9)

(10)

步驟4 更新參數μ

(11)

其中，ρ>1為常數；ε>0為比較小的正數。

不難發現式(9)仍是非凸優化問題，難以求解，為求解式(9)最優解，我們引入無約束DC規劃問題，利用DCA算法[5]，將非凸優化問題轉化為凸優化問題，快速求解大規模優化問題。

在DCA算法中，無約束DC規劃問題如式(12)所示

minf(x)=g(x)-h(x)s.t.x∈RN

(12)

其中，g和h是RN上的下半連續凸函數。

minf(x)=g(x)-h(x)?

min{g(x)-x,vk:x∈Rn,vk∈?h(xk)}

(13)

那么，對于式(9)的求解，可令

(14)

則

(15)

(16)

其實，迭代更新A和E時外層循環無需解出精確值，只需得到子問題的近似解就可以滿足最優解的收斂條件，即低秩矩陣A和稀疏矩陣E的非精確更新公式為

Ak+1=argmin Γ(A,Ek+1,Yk,μk)=D1/μk(D-Ek+1+Yk/μk)

(17)

Ak+1+Yk/μk)

(18)

其中

(19)

2.2 采用WNNM的低秩矩陣優化算法

上文已經介紹了基于分式函數的RPCA算法的實現過程，不難發現，借助奇異值算子更新低秩矩陣時，使用了固定常數ε對閾值進行收縮，但是在實際應用中，大的奇異值表征數據的主要成分。為了提高矩陣的低秩性，本文對奇異值根據大小進行了不同收縮，賦予相應的收縮閾值，更符合實際情況，實現更好的恢復效果。

Candes等提出了一種重加權范數最小化的思想，提升了在稀疏恢復方面的效果。之后，Gu等提出了加權核范數最小化算法(weighted nuclear norm minimization，WNNM)[6]，其優化模型表示如下

(20)

(21)

式中：c為大于0的常數。σi(L)的初始值設為

(22)

因此，在本文改進RPCA算法中更新低秩矩陣A時，修改式(8)中的奇異值算子如式(23)所示

SW(∑)=max(∑ii-wi,0)

(23)

2.3 改進RPCA的整體算法

正如上兩節所述，在本文改進的RPCA算法中，設計采用分式函數對l0范數逼近，結合DCA算法將非凸優化問題轉化為凸優化問題求解，進一步提高稀疏矩陣恢復成功率，更適用于大規模遙感影像去云；同時引入WNNM算法中的權值向量對奇異值算子實現自適應收縮，增加地物矩陣低秩性，使算法具有更好的低秩矩陣恢復效果。改進RPCA的整體算法流程如圖2所示。

圖2 改進RPCA的整體算法流程

改進的RPCA的整體算法如下所示：

步驟1 初始化Y0，A0=0，μ0，k=0，b=1.3

步驟2 while 不收斂 do

步驟4 更新權值向量W

步驟5SW(∑)=max(∑ii-wi,0)

基于此，本文選取東三省及西南地區部分省市A股上市公司2014～2015年數據為樣本，實證檢驗影響企業研發支出資本化的相關因素。與以往研究相比，本文的主要貢獻在于創新性地提出了系統風險這一因素對研發支出的影響，同時為了檢驗不同行業不同因素對資本化的影響程度，將企業行業類別分為高新技術企業和非高新技術企業，對比分析不同企業環境下研發支出資本化的不同結果，為客觀觀察企業會計行為提供了一定的理論及實踐依據。

步驟6Ak+1=USW(∑)VT

步驟7

步驟9Yk+1=Yk+μk(D-Ak+1-Ek+1)

步驟10 更新μk

步驟11k=k+1

步驟12 end while

3 實驗結果與分析

3.1 遙感影像去云算法設計

針對遙感影像云層遮擋區域的復雜情況，本文設計一種采用改進RPCA算法的遙感影像去云算法。

由于技術限制，獲取同一區域的多幅遙感影像時，地理空間上會存在一定角度偏移。此外，由于光照、局部地理環境、拍攝時間等多種因素的干擾，會導致影像亮度不均。首先，基于精度和色彩一致性的要求，本文的改進RPCA去云算法對遙感圖像序列做預處理，即采用SIFT算法的圖像配準[7]以及Wallis勻色[8]。其次，針對云區，通過構造低秩觀測矩陣，利用改進的RPCA算法去除；最后重構出無云遙感影像序列，實現了遙感影像復雜稀疏云區的聯合去云處理。

具體去云算法流程如圖3所示。

圖3 采用改進RPCA的遙感影像去云算法流程

3.2 實驗結果與分析

為了驗證本文采用改進RPCA的去云算法在復雜稀疏云區的去云效果，本文采用LANDSAT 8衛星采集的一組不同時相的10張遙感影像數據作為實驗原始影像進行去云處理實驗，與傳統去云算法(同態濾波法[9]、小波分解法[10])以及常見RPCA算法如APG[11]和IALM[12]進行對比實驗，原始影像數據如圖4所示。

本文實驗環境為：Windows 8系統，Intel(R)Core(TM)i7-4720HQ CPU@2.60 GHz,8 GB內存的PC機器，算法采用MATLAB R2014b編程實現。

為比較復雜稀疏云區的去云效果，選取圖4原始遙感影像中圖4(c)和圖4(j)兩幅云特征鮮明的影像進行分析，去云結果如圖5和圖6所示。

圖4 原始影像數據

圖5 影像c去云結果對比

圖6 影像j去云結果對比

從圖5和圖6中可以明顯看出，同態濾波法和小波分解法對部分薄云區域有一定去除能力，但對小塊云層無法有效去除；而RPCA算法(APG、IALM)及本文改進的RPCA算法對于稀疏云區(小片薄云、小片厚云、突變云區以及成分復雜的混合云區)，有非常明顯的去除效果。其中，本文改進的RPCA算法處理結果圖在云層去除、景物信息增強及清晰區域的保護方面的表現均優于APG算法和IALM算法。

為了定性評價遙感圖像去云效果，我們采用圖像灰度均值、標準差和熵這3項指標參數來分析。

(1)均值

(24)

圖像均值表征圖像平均灰度值，稀疏云區灰度值比較高，因此去云之后均值會下降。式(24)中f(x,y)表示大小為M×N的圖像矩陣各像素點灰度值。

(2)標準差

(25)

標準差表示各像素點灰度值偏離均值的程度，去云之后遙感圖像各像素點之間灰度值更加接近，標準差相應變小。

(3)熵

(26)

熵用來度量圖像所含信息量的多少，在去云處理中，必須保證主要信息不丟失，所以熵也是一個很重要的的評價指標。式(26)中Pi表示第i級灰度值的出現概率。

根據上述3個指標，我們針對圖5和圖6，計算出具體指標參數值，分別見表1和表2。

表1 影像c去云結果比較

表2 影像j去云結果比較

可以看出，采用本文算法對遙感影像處理后，圖像灰度均值和標準差均有很大程度下降，說明去云效果明顯，而且熵值略有增加，表明圖像信息量得到了保護，圖像更加清晰，符合預估。同時，在均值、標準差和熵3項指標上，本文改進的RPCA算法均占優勢，處理后的圖像具有最低的均值和標準差，同時熵值最大，彌補了傳統算法易造成的信息量的丟失。

下面分析APG、IALM及本文改進的RPCA算法在運行時間、迭代次數的數據比較，見表3。

表3 RPCA算法性能比較

改進的RPCA算法運算速度更快，需要更少的迭代次數就可以實現收斂，具有更好的建模效果。綜合實驗結果表明，本文采用改進的RPCA遙感影像去云算法，在復雜的稀疏云區去云方面比傳統去云算法有更好的表現；同時，在去云效果和運行時間方面，均優于常見RPCA算法，具有更好的實際應用價值。

4 結束語

本文利用RPCA算法模型來解決遙感影像復雜的稀疏云區的去云問題，將傳統去云處理中需要解決的云層厚度問題轉化為恢復遙感影像的低秩和稀疏成分；同時，對原有算法中l0范數采取一種基于分式函數的新的逼近方式，引入了加權核范數對奇異值閾值進行自適應的調節，提高稀疏成分恢復的成功率，增大低秩矩陣的低秩性，使改進算法更適用于大規模實際問題。

實驗結果及其分析表明，本文提出的采用改進RPCA的遙感影像去云算法，能夠有效解決稀疏云區(小片薄云、小片厚云、突變云區以及成分復雜的混合云區)的去云問題，同時還能保護圖像的信息量，算法穩定，在主觀視覺效果和客觀評價標準方面優于其它方法，有效實現了遙感影像復雜稀疏云區的去云處理。

參考文獻：

[1]LEI Lixia,ZHANG Yuejin,HUANG Dechang.Compressed sensing applications in medical image compression[J].Journal of Computer Engineering and Design,2014,35(11):3898-3902(in Chinese).[雷莉霞,張躍進,黃德昌.壓縮感知在醫學圖像壓縮中的應用[J].計算機工程與設計,2014,35(11):3898-3902.]

[2]Chen Y,Jalali A,Sanghavi S,et al.Low-rank matrix reco-very from errors and erasures[J].IEEE Transactions on Information Theory,2013,59(7):4324-4337.

[3]Qaisar S,Bilal R M,Iqbal W,et al.Compressive sensing:From theory to applications,a survey[J].Journal of Communications & Networks,2013,15(5):443-456.

[4]Deng W,Yin W.On the global and linear convergence of the generalized alternating direction method of multipliers[J].Journal of Scientific Computing,2016,66(3):889-916.

[5]Thi HAL,Le H M,Tao P D.Feature selection in machine lear-ning: An exact penalty approach using a difference of convex function algorithm[J].Machine Learning,2015,101(1):163-186.

[6]Gu S,Zhang L,Zuo W,et al.Weighted nuclear norm minimization with application to image denoising[C]//IEEE Conference on Computer Vision and Pattern Recognition.Columbus:IEEE,2014:2862-2869.

[7]Kupfer B,Netanyahu N S,Shimshoni I.An efficient sift-based mode-seeking algorithm for sub-pixel registration of remotely sensed images[J].IEEE Geoscience & Remote Sensing Letters,2015,12(2):379-383.

[8]ZHU Qiaoyun,DA Xing.Dodging method for multi-source remote sensing images based on wallis filter[J].Geomatics & Spatial Information Technology,2012,35(10):130-132(in Chinese).[朱巧云,答星.基于Wallis濾波器的異源遙感影像勻光方法[J].測繪與空間地理信息,2012,35(10):130-132.]

[9]Shen H,Li H,Qian Y,et al.An effective thin cloud removal procedure for visible remote sensing images[J].ISPRS Journal of Photogrammetry & Remote Sensing,2014,96(11):224-235.

[10]Wang X X,Jiang L S,Chen Y P,et al.Thin cloud removal from remote sensing images with adaptive thresholds of wavelet transforms[J].Journal of the University of Electronic Science & Technology of China,2013,42(3):390-393.

[11]Xue Q,Wang H,Cui Z,et al.Electrical capacitance tomography using an accelerated proximal gradient algorithm[J].Review of Scientific Instruments,2012,83(4):043704.

[12]Bhardwaj A,Raman S.Robust PCA-based solution to image composition using augmented Lagrange multiplier(ALM)[J].The Visual Computer,2016,32(5):591-600.