基于薄云厚度分布評估的遙感影像高保真薄云去除方法

汪月云,黃 微,王 睿

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

0 引言

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

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

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

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

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

通過多幅影像的實驗發現，與文獻[14]算法去云結果、文獻[15]算法去云結果，以及文獻[16]算法去云結果相比，本文方法的結果影像去云效果良好，同時具有數據高保真的優點。

1 加性模型薄云去除

1.1 云成像模型

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

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

(1)

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

1.2 薄云估計

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

(2)

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

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

Himpr(x)=Hg(x)-minHg

(3)

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

(4)

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

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

圖1 文獻[16]算法的薄云估計和去云結果Fig.1 Haze estimation and haze removal results of algorithm in literature [16]

2 薄云去除改進

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

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

2.1 無云區域保真

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

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

(5)

其中：

(6)

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

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

t(x)∈(0,1]

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

(7)

2.2 藍色地物云量估計

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

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

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

(8)

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

(9)

Pf=P1∩P2

(10)

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

(11)

其中:Hblue是藍色地物薄云厚度;Hnew(x)是優化的薄云厚度分布;Dnear是在Hnew(x)中藍色地物附近人工選取出的薄云區域;n是Dnear中像素點個數。

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

(12)

為了去除Htotal(x)中的塊狀效應，本文應用導向濾波器[20]將其平滑得到最終的薄云厚度分布Hfinal(x)。

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

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

2.3 高保真無云影像

通過逐步優化薄云分布，單獨估計藍色地物云量，最終得到的無云影像J表示如下：

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

(13)

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

3 實驗結果及分析

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

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

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

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

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

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

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

圖像波段降質影像文獻[14]方法文獻[15]方法文獻[16]方法本文方法圖4藍波段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藍波段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藍波段15.1713.6313.7613.9514.35綠波段13.2913.7212.8113.0912.37紅波段10.029.459.709.9810.01平均值12.8312.2712.0912.3412.24

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

(14)

其中：Dclea是降質影像I(x)和恢復的無云影像J(x)的同一片無云區域；m是無云區內像素點個數。δ值越小保真效果越好，該評價指標MRE的統計值如表2所示。

表2 不同方法MRE統計值對比Tab. 2 Comparison of MRE statistical values of different methods

由表2可看出，四種去云方法都具有高度的保真特性，文獻[15]方法去云結果的數據保真性能相比其他方法略差一些。本文方法的無云區域保真性能最好。綜合兩個評價指標可以看出，本文方法能得到高保真的最優去云結果。

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

4 結語

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