汪月云,黃 微,王 睿

(上海大學(xué) 通信與信息工程學(xué)院，上海 200444)(*通信作者電子郵箱yueyunW@163.com)

0 引言

光學(xué)遙感影像往往受到薄云污染的影響，造成影像地物細(xì)節(jié)對比度下降，這給影像的應(yīng)用[1-3]帶來了很多困難。因此，如何有效去除降質(zhì)影像的薄云影響從而提高影像的利用率是遙感影像應(yīng)用的一項(xiàng)重要研究。

目前，光學(xué)遙感影像薄云去除主要有頻域法和空域法兩類。頻域法去云主要根據(jù)薄云信息處于影像低頻區(qū)域的原理，所以可以設(shè)計(jì)合適的高通濾波器以濾除低頻成分從而有效去云。比如:閻慶等[4]結(jié)合小波變換和同態(tài)濾波、馮春等[5]利用同態(tài)濾波都達(dá)到了有效去云效果。

空域法去云將薄云污染分為乘性污染和加性污染兩種模型。乘性污染經(jīng)典的去云方法有暗通道先驗(yàn)法[6-7]和薄云優(yōu)化變換(Haze Optimized Transformation, HOT)法[8-9]。暗通道先驗(yàn)法通過對清晰影像庫統(tǒng)計(jì)得到暗通道分布規(guī)律，即清晰影像至少在某一個(gè)顏色通道中有灰度值接近于0的像素，利用該先驗(yàn)規(guī)律得出云模型的大氣光和透射率從而進(jìn)行薄云去除。Pan等[10]研究發(fā)現(xiàn)無云遙感影像的暗通道數(shù)值很低，但不接近于0，通過改進(jìn)遙感影像暗通道分布從而有效去除薄云。HOT法是利用遙感影像紅光波段和藍(lán)光波段相關(guān)性高的特點(diǎn)尋找合適的晴空線，用像元到晴空線的距離表征該點(diǎn)云量進(jìn)而有效去云。He等[11]提出了AHOT(Advanced HOT)算法，先利用空間信息剔除易混淆地物,提高云檢測精度，再考慮大氣散射乘性效應(yīng)，提出基于AHOT的虛擬云點(diǎn)(Virtual Cloud Point, VCP)技術(shù)擴(kuò)大由薄云污染降低的數(shù)字量化值(Digital Number, DN)方差從而有效去云。Makarau等[12]利用加性薄云污染模型，在傳統(tǒng)的暗目標(biāo)減法(Dark Object Subtraction, DOS)基礎(chǔ)上[13]，通過局部搜索整個(gè)影像的暗目標(biāo)構(gòu)建影像的薄云厚度分布(Haze Thickness Map, HTM)，最后用原降質(zhì)影像減去薄云厚度分布得到無云影像。這些去云算法大多關(guān)注去云后影像細(xì)節(jié)的銳化、地物間對比度的增強(qiáng)，如何在去云的同時(shí)保證數(shù)據(jù)保真是一項(xiàng)具有挑戰(zhàn)性的問題。

近年來，已有學(xué)者提出了一些考慮保真的薄云去除算法[14-16]。Shen等[14]基于同態(tài)濾波法半自動化地確定不同波段的截止頻率，對薄云像素和無云像素分別處理從而獲得高保真的無云影像。Singh等[15]基于暗通道先驗(yàn)法提出改進(jìn)的無云影像恢復(fù)算法，通過改進(jìn)聯(lián)合三邊濾波器以提高大氣光的精細(xì)估計(jì)效果，并重新定義透射率從而降低無云影像的顏色失真程度。Liu等[16]基于加性薄云污染模型，將原始薄云厚度分布先經(jīng)過紋理去除，再進(jìn)行大面積地物反射抑制，然后用導(dǎo)向?yàn)V波器進(jìn)行平滑，最后用線性回歸求取不同波段的薄云厚度從而得到無云影像。文獻(xiàn)[16]算法的去云結(jié)果整體具有保真效果，但是依舊存在無云區(qū)域過度校正，以及影像中藍(lán)色地物色彩嚴(yán)重失真的問題。

因此，本文在文獻(xiàn)[16]算法的基礎(chǔ)上提出了一種改進(jìn)的薄云去除方法，即基于薄云厚度分布(HTM)評估的高保真薄云去除方法。該方法首先利用基于暗通道先驗(yàn)薄云去除算法求得影像的透射率分布圖[10]，根據(jù)透射率分布圖中像素?cái)?shù)值接近于1來自動確定無云區(qū)域的范圍。然后,在HTM中，求取無云區(qū)域的平均值,將整個(gè)HTM減去無云區(qū)域平均值，使得HTM滿足無云區(qū)域的薄云厚度接近于零，以確保無云區(qū)域的保真度。針對藍(lán)色地物區(qū)域云量被過估計(jì)，造成去云結(jié)果中色彩失真的問題，本文利用藍(lán)色地物RGB波段的光譜特性，結(jié)合閾值法檢測出藍(lán)色地物，用藍(lán)色地物附近的HTM平均值作為藍(lán)色地物區(qū)域的云量估計(jì)。最后,用降質(zhì)影像減去不同波段優(yōu)化的HTM得到最終去云結(jié)果。

通過多幅影像的實(shí)驗(yàn)發(fā)現(xiàn)，與文獻(xiàn)[14]算法去云結(jié)果、文獻(xiàn)[15]算法去云結(jié)果，以及文獻(xiàn)[16]算法去云結(jié)果相比，本文方法的結(jié)果影像去云效果良好，同時(shí)具有數(shù)據(jù)高保真的優(yōu)點(diǎn)。

1 加性模型薄云去除

1.1 云成像模型

薄云污染可以認(rèn)為是加性污染[12，16],并且薄云污染對不同波段的影響不同[16-17],所以云成像模型可以描述為：

Ic(x)=Jc(x)+Hc(x);c∈{R,G,B}

(1)

其中：x代表像素的位置；I代表降質(zhì)影像；J代表無云影像；H代表薄云厚度分布；c是{R,G,B}中的任一波段。

1.2 薄云估計(jì)

在Makarau等[12]的研究中，薄云厚度分布(HTM)估計(jì)如下：

(2)

其中：Ib(x)、Ig(x)分別是降質(zhì)影像的藍(lán)光波段和綠光波段；W(x)是以x為中心、大小為w×w像素的窗口。本文實(shí)驗(yàn)中該窗口大小取3×3像素。

Liu等[16]發(fā)現(xiàn)由式(2)得到的Hraw(x)存在將地物反射光作為薄云組成成分，造成降質(zhì)影像被過度矯正的問題。為了減少地物反射光的影響，文獻(xiàn)[16]對Hraw(x)進(jìn)一步改進(jìn)：首先為了去除Hraw(x)中的地物紋理，文獻(xiàn)[16]采用了可以去除無用紋理和細(xì)節(jié)的全變差(Total Variation, TV)正則化圖像修復(fù)模型[18-19]得到平滑的薄云厚度分布Hs(x)；接著,為解決大面積高亮地物的薄云厚度過估計(jì)問題，文獻(xiàn)[16]用Cannny算子檢測出高亮地物邊緣，用局部窗口對Hs(x)中高亮地物區(qū)域作數(shù)值抑制得到Hd(x)；然后,應(yīng)用導(dǎo)向?yàn)V波器[20]去除Hd(x)中的塊狀效應(yīng)得到Hg(x)；最后,求取不同波段的薄云厚度分布。

Himpr(x)=Hg(x)-minHg

(3)

Hc(x)=kcHimpr(x)；c∈{R,G,B}

(4)

為了去除氣溶膠的影響，式(3)用Hg減去其最小值minHg得到最終的薄云厚度估計(jì)Himpr(x)。式(4)中：kc是不同波段的云估計(jì)強(qiáng)度，它是在Hg低于平均值的像素集與降質(zhì)影像I的不同波段中低于平均值的像素集的交集中，以Hg的像素為自變量，以影像的不同波段的像素為因變量，通過線性回歸得到；Hc(x)是不同波段的薄云厚度估計(jì)。得到Hc(x)后，通過式(1)即可得到無云影像J。

圖1是根據(jù)文獻(xiàn)[16]算法得到的實(shí)驗(yàn)結(jié)果。其中:圖1(a)是大小為300×400像素的GF-2降質(zhì)影像;圖1(b)是最終的云厚度分布估計(jì)Himpr(x);圖1(c)是去云結(jié)果。由圖1可以看出，Himpr(x)中盡管云厚度分布范圍估計(jì)比較準(zhǔn)確，但是無云區(qū)域的云厚度值并不為0或接近于0，同時(shí)無云區(qū)域的藍(lán)色和白色物體估計(jì)的云厚度分布值也很高，這就造成了去云影像的無云區(qū)域存在過度矯正和藍(lán)色地物區(qū)域嚴(yán)重失真的問題。

圖1 文獻(xiàn)[16]算法的薄云估計(jì)和去云結(jié)果Fig.1 Haze estimation and haze removal results of algorithm in literature [16]

2 薄云去除改進(jìn)

為解決圖1去云結(jié)果中反映出來的無云區(qū)域過度矯正和藍(lán)色地物失真嚴(yán)重問題，基于加性薄云污染模型，本文在文獻(xiàn)[16]得到的最終薄云厚度分布Himpr(x)基礎(chǔ)上進(jìn)行改進(jìn)，在獲得良好去云結(jié)果的同時(shí)確保影像的高保真。本文方法的流程如圖2所示。

圖2 本文方法流程Fig. 2 Flow chart of the proposed method

2.1 無云區(qū)域保真

首先針對無云區(qū)域薄云厚度被過估計(jì)，造成無云區(qū)域數(shù)據(jù)失真的問題，本文對Himpr(x)作了如下改進(jìn)：

Hnew(x)=max(Himpr(x)-K,0)

(5)

其中：

(6)

式中：Dclean是Himpr(x)中的無云區(qū)域;M是Dclean中的像素個(gè)數(shù);K是無云區(qū)域的平均值。理論上，Hnew(x)的無云區(qū)域DN值接近或等于0，有云區(qū)域DN值大于0，為了避免由于Himpr(x)-K導(dǎo)致負(fù)數(shù)的出現(xiàn)，這里使用了最大值算子max。

無云區(qū)域Dclean的確定影響著最終無云影像無云區(qū)域的保真度。為了有效確定無云區(qū)域Dclean的范圍，本文利用基于暗通道先驗(yàn)薄云去除算法[6-7]中的透射率來估計(jì)該區(qū)域。其中透射率為：

t(x)∈(0,1]

球形工頻電場檢測傳感器采用球形電容檢測感應(yīng)電壓的方法來檢測電場,一維球形傳感器的結(jié)構(gòu)如圖1所示,將一個(gè)空心金屬球殼切割成兩部分,通過絕緣物質(zhì)將兩個(gè)半球連接在一起,并在兩個(gè)半球中間并聯(lián)一個(gè)測量電容,兩個(gè)半球分別構(gòu)成了傳感器的兩個(gè)電極[6-7]。

(7)

2.2 藍(lán)色地物云量估計(jì)

對于降質(zhì)影像中的藍(lán)色地物，由于像素在藍(lán)光波段的亮度值遠(yuǎn)大于綠光波段，因此根據(jù)式(1)得到的薄云厚度分布Hraw(x)中容易造成藍(lán)色地物薄云厚度被過估計(jì)，最終導(dǎo)致去云結(jié)果中藍(lán)色地物色彩失真。

通過大量實(shí)驗(yàn)數(shù)據(jù)發(fā)現(xiàn)，在優(yōu)化后的薄云厚度分布Hnew(x)中，白色地物和藍(lán)色地物區(qū)域的DN值都很高。因此，本文先利用閾值選定白色區(qū)域和藍(lán)色區(qū)域，再根據(jù)原始影像中藍(lán)色地物藍(lán)光波段像素值遠(yuǎn)大于紅光波段像素值的特點(diǎn)[21]，剔除白色地物區(qū)域最終得到藍(lán)色地物的像素集：

P1={x|Hnew(x)≥T}

(8)

P2={x|Ib(x)-Ir(x)>0}

(9)

Pf=P1∩P2

(10)

其中：x代表像素坐標(biāo)位置；P1是包含藍(lán)色地物和白色地物的像素集；P2是包含偏藍(lán)地物的像素集；Pf是藍(lán)色地物像素集。通過閾值T可以選定白色地物區(qū)域和藍(lán)色地物區(qū)域，閾值T的選定不能過大，否則影像中的藍(lán)色地物區(qū)域不能全部選中，閾值T也不能過小，否則容易包含非藍(lán)色地物。本文中閾值T取Hnew(x)中最大像素值的70%，即T=max(Hnew(x))*0.7。Ib(x)、Ir(x)分別是原始影像的藍(lán)光波段和紅光波段。由于薄云的分布一般是連續(xù)和平緩的，因此本文將藍(lán)色地物區(qū)域的薄云厚度用周圍區(qū)域薄云厚度的均值代替，其表達(dá)式為：

(11)

其中:Hblue是藍(lán)色地物薄云厚度;Hnew(x)是優(yōu)化的薄云厚度分布;Dnear是在Hnew(x)中藍(lán)色地物附近人工選取出的薄云區(qū)域;n是Dnear中像素點(diǎn)個(gè)數(shù)。

經(jīng)過以上的優(yōu)化，可得到完整的薄云厚度分布Htotal(x)：

(12)

為了去除Htotal(x)中的塊狀效應(yīng)，本文應(yīng)用導(dǎo)向?yàn)V波器[20]將其平滑得到最終的薄云厚度分布Hfinal(x)。

本文方法的每一個(gè)優(yōu)化步驟結(jié)果如圖3所示。其中：圖3(b)是根據(jù)透射率得到的無云區(qū)域二值圖像，白色表示無云區(qū)域。圖3(c)是優(yōu)化后的薄云厚度分布Hnew(x)?？梢钥闯?，Hnew(x)在無云區(qū)域的DN值接近或等于0，這能保證恢復(fù)的無云影像在無云區(qū)域的高保真度。圖3(d)是藍(lán)色地物的檢測結(jié)果, 可以發(fā)現(xiàn)，圖中的藍(lán)色地物位置檢測比較準(zhǔn)確。圖3(e)是最終的薄云厚度分布Hfinal(x)。與圖3(c)相比，圖3(e)的Hfinal(x)在藍(lán)色地物區(qū)域的DN值相比Hnew(x)小很多，可以避免藍(lán)色地物云量被過估計(jì)。

圖3 本文方法各步驟結(jié)果Fig. 3 Results of each step of the proposed method

2.3 高保真無云影像

通過逐步優(yōu)化薄云分布，單獨(dú)估計(jì)藍(lán)色地物云量，最終得到的無云影像J表示如下：

Jc(x)=Ic(x)-kcHfinal(x)

(13)

其中：I是降質(zhì)影像；kc是根據(jù)文獻(xiàn)[16]算法中通過線性回歸得到的不同波段云估計(jì)強(qiáng)度；Hfinal(x)是最終優(yōu)化的薄云厚度分布。圖3(f)是本文改進(jìn)方法得到的無云影像，對比圖3(a)的降質(zhì)影像可以發(fā)現(xiàn)，圖3(f)的無云影像不但具有良好的去云效果，而且無云區(qū)域的地物色彩基本不變，藍(lán)色地物色彩自然，即本文方法得到的無云云影像具有高保真的特點(diǎn)。

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

本文采用三幅大小均是512×512遙感影像作為實(shí)驗(yàn)對象。圖4是分辨率為30 m的Landsat7遙感影像；圖5是分辨率為5.8 m的ZY3遙感影像；圖6是分辨率為2.0 m的WV-2遙感影像。數(shù)據(jù)包括了城市、農(nóng)田、草地等多種地物類型。本文采用文獻(xiàn)[14]方法、文獻(xiàn)[15]方法以及文獻(xiàn)[16]方法與本文方法進(jìn)行比較。

圖4 不同方法去云結(jié)果比較(Landsat7遙感影像)Fig. 4 Haze removal result comparison of different methods (Landsat7 remote sensing image)

圖5 不同方法去云結(jié)果比較(ZY3遙感影像)Fig. 5 Haze removal result comparison of different methods (ZY3 remote sensing image)

從圖4～6可看出，4種算法都具有良好的去云效果。文獻(xiàn)[14]去云方法會使降質(zhì)影像的云區(qū)變得有些模糊，使得地物細(xì)節(jié)不夠清晰，且去云影像的無云區(qū)邊界色彩過渡不自然。文獻(xiàn)[15]去云方法對影像地物數(shù)據(jù)的保真度較差，偏綠地物太過鮮艷，色彩不自然。文獻(xiàn)[16]去云方法對非藍(lán)色地物整體有一定的保真效果，然而在藍(lán)色地物區(qū)域(圖5(d)和圖6(d)中方塊)去云后色彩與原圖中的藍(lán)色地物(圖5(a)和圖6(a)中方塊)相比過于鮮艷，失真嚴(yán)重。本文方法得到的去云結(jié)果不但在無云區(qū)域有較高的保真度，而且影像中的藍(lán)色地物(圖5(e)和圖6(e)中方塊)去云后色彩自然，達(dá)到高保真去云目的。

為了定量評估本文方法的去云效果，采用廣義奇異值分解(Generalized Singular Value Decomposition, GSVD)作為評價(jià)因子[22]。GSVD是用于評估圖像受到模糊或噪聲影響的程度,薄云覆蓋是光學(xué)遙感影像被模糊的主要原因。GSVD值越小，則影像的去云效果越好。

不同方法的GSVD統(tǒng)計(jì)值如表1所示。從表1可看出，四種去云方法中，本文方法對應(yīng)的三波段GSVD平均值總體較小，表明本文方法具有理想的去云效果。

圖6 不同方法去云結(jié)果比較(WV-2遙感影像)Fig. 6 Haze removal result comparison of different methods (WV-2 remote sensing image)表1 不同方法的GSVD統(tǒng)計(jì)值對比Tab. 1 Comparison of GSVD statistical values of different methods

圖像波段降質(zhì)影像文獻(xiàn)[14]方法文獻(xiàn)[15]方法文獻(xiàn)[16]方法本文方法圖4藍(lán)波段25.8524.6225.1824.7925.14綠波段24.5824.2323.9824.0424.08紅波段23.3523.5623.3123.2123.20平均值24.6024.1424.1624.0124.14圖5藍(lán)波段22.5120.8118.6719.1812.15綠波段17.3514.8615.4716.9317.30紅波段19.6819.1518.5019.6919.44平均值19.8518.2717.5518.6016.30圖6藍(lán)波段15.1713.6313.7613.9514.35綠波段13.2913.7212.8113.0912.37紅波段10.029.459.709.9810.01平均值12.8312.2712.0912.3412.24

為了定量評價(jià)降質(zhì)影像無云區(qū)域的保真度，本文采用平均相對誤差(Mean Relative Error, MRE)作為評價(jià)指標(biāo)，結(jié)果記作δ:

(14)

其中：Dclea是降質(zhì)影像I(x)和恢復(fù)的無云影像J(x)的同一片無云區(qū)域；m是無云區(qū)內(nèi)像素點(diǎn)個(gè)數(shù)。δ值越小保真效果越好，該評價(jià)指標(biāo)MRE的統(tǒng)計(jì)值如表2所示。

表2 不同方法MRE統(tǒng)計(jì)值對比Tab. 2 Comparison of MRE statistical values of different methods

由表2可看出，四種去云方法都具有高度的保真特性，文獻(xiàn)[15]方法去云結(jié)果的數(shù)據(jù)保真性能相比其他方法略差一些。本文方法的無云區(qū)域保真性能最好。綜合兩個(gè)評價(jià)指標(biāo)可以看出，本文方法能得到高保真的最優(yōu)去云結(jié)果。

通過對不同數(shù)據(jù)范圍、不同分辨率大小、多種地物類型的薄云污染影像數(shù)據(jù)實(shí)驗(yàn)得出，本文方法得到的無云影像都具有數(shù)據(jù)高保真效果，表明本文方法對不同薄云污染影像去云具有高度適用性。

4 結(jié)語

本文提出了一種基于薄云厚度評估的高保真去薄云方法，不僅能夠滿足無云影像中無云區(qū)域數(shù)據(jù)的高保真度，還能夠解決傳統(tǒng)算法藍(lán)色地物色彩失真問題。針對影響算法的兩大因素：無云區(qū)域Dclean的確定和閾值T的選取，本文通過透射率分布圖中像素?cái)?shù)值接近于1實(shí)現(xiàn)了自動化確定無云區(qū)域的范圍，然而在確定藍(lán)色區(qū)域的過程中，本文是以經(jīng)驗(yàn)閾值T選取藍(lán)色地物和白色地物區(qū)域，自適應(yīng)求取最優(yōu)閾值T將是我們后續(xù)的改進(jìn)工作。