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

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

0 引 言

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

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

本文針對同一區(qū)域、不同時相遙感影中的稀疏云層遮擋，提出一種采用改進的RPCA去云算法，即對原有RPCA算法中l(wèi)0范數(shù)采取一種新的分式函數(shù)逼近方式，同時引入加權(quán)核范數(shù)對奇異值閾值進行自適應的調(diào)節(jié)，滿足保留數(shù)據(jù)主要成分的實際需要，使改進的RPCA算法在處理遙感影像去云問題時更為有效，實現(xiàn)復雜的稀疏云區(qū)的聯(lián)合去云處理。

1 問題的提出

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

(1)

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

(2)

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

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

2 改進的RPCA算法研究

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

2.1 基于分式函數(shù)的稀疏矩陣優(yōu)化算法

2.1.1 模型改進

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

(3)

凸松弛為

(4)

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

(5)

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

圖1 分式函數(shù)圖像

(6)

2.1.2 模型求解過程

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

(7)

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

(8)

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

(9)

(10)

步驟4 更新參數(shù)μ

(11)

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

不難發(fā)現(xiàn)式(9)仍是非凸優(yōu)化問題，難以求解，為求解式(9)最優(yōu)解，我們引入無約束DC規(guī)劃問題，利用DCA算法[5]，將非凸優(yōu)化問題轉(zhuǎn)化為凸優(yōu)化問題，快速求解大規(guī)模優(yōu)化問題。

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

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

(12)

其中，g和h是RN上的下半連續(xù)凸函數(shù)。

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

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

(13)

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

(14)

則

(15)

(16)

其實，迭代更新A和E時外層循環(huán)無需解出精確值，只需得到子問題的近似解就可以滿足最優(yōu)解的收斂條件，即低秩矩陣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的低秩矩陣優(yōu)化算法

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

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

(20)

(21)

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

(22)

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

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

(23)

2.3 改進RPCA的整體算法

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

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

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

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

步驟2 while 不收斂 do

步驟4 更新權(quán)值向量W

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

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

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

步驟7

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

步驟10 更新μk

步驟11k=k+1

步驟12 end while

3 實驗結(jié)果與分析

3.1 遙感影像去云算法設(shè)計

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

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

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

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

3.2 實驗結(jié)果與分析

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

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

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

圖4 原始影像數(shù)據(jù)

圖5 影像c去云結(jié)果對比

圖6 影像j去云結(jié)果對比

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

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

(1)均值

(24)

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

(2)標準差

(25)

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

(3)熵

(26)

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

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

表1 影像c去云結(jié)果比較

表2 影像j去云結(jié)果比較

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

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

表3 RPCA算法性能比較

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

4 結(jié)束語

本文利用RPCA算法模型來解決遙感影像復雜的稀疏云區(qū)的去云問題，將傳統(tǒng)去云處理中需要解決的云層厚度問題轉(zhuǎn)化為恢復遙感影像的低秩和稀疏成分；同時，對原有算法中l(wèi)0范數(shù)采取一種基于分式函數(shù)的新的逼近方式，引入了加權(quán)核范數(shù)對奇異值閾值進行自適應的調(diào)節(jié)，提高稀疏成分恢復的成功率，增大低秩矩陣的低秩性，使改進算法更適用于大規(guī)模實際問題。

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

參考文獻：

[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).[雷莉霞,張躍進,黃德昌.壓縮感知在醫(yī)學圖像壓縮中的應用[J].計算機工程與設(shè)計,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.