李 力, 唐曉霏, 劉 雯, 陳 城
(1. 湖北省氣象信息與技術(shù)保障中心, 武漢 430074; 2. 湖北省測繪工程院, 武漢 430074)
無人機遙感技術(shù)目前已廣泛應用于氣象環(huán)境監(jiān)測、 國土測繪勘察、 水文地質(zhì)救災等領域, 但由于無人機主要為中低空飛行, 而水汽云霧和霧霾天氣都會影響可見光和紅外線在大氣中的傳播, 因此使其獲取的圖像模糊有霧[1]. Oakley等[2-3]首次提出了有霧圖像的退化處理, 目前對單幅圖像的去霧算法研究已有很多結(jié)果[4-5]. He等[6]采用暗通道先驗(dark channel prior, DCP)原理量化大氣透射率, 獲取大氣透射率圖層(transmission map, TM), 再用軟摳圖方法得到相對準確的大氣透射率, 建立物理模型恢復清晰影像. 但該方法對場景中存在大量白色區(qū)域的情況去霧效果不佳, 相關算法復雜程度較高[7-8]、 計算速度較慢, 難以適應航拍圖像數(shù)據(jù)量大、 及時性強等特點.
通過分析基于暗通道先驗的圖像去霧算法, 發(fā)現(xiàn)運行中耗時較多的兩個環(huán)節(jié)是計算暗通道圖像和使用引導濾波精細化透射率圖層. 本文針對引導濾波精細化大氣透射率圖層時間較長的問題, 考慮采用與濾波窗口無關、 算法時間復雜度為O(1)的中值濾波精細化透射率圖層. 改進算法通過計算機編碼, 利用無參考的客觀圖像質(zhì)量評價方法進行量化評價, 將信息熵、 邊緣強度和方差作為評價圖像質(zhì)量的3個指標, 以在去霧圖像質(zhì)量和后期處理應用影響較小的條件下節(jié)省計算時間, 提高航拍圖像暗通道先驗去霧算法的處理速度.
在計算機視覺圖形學領域, 有霧圖像的退化模型[9]表示為
I(x)=J(x)t(x)+A(1-t(x)),
(1)
其中:I(x)為傳感器接收到的場景信號, 即輸入的有霧影像; 場景輻射J(x)為信號處理后的清晰圖像;t(x)為大氣透射率,J(x)t(x)稱為直接衰減項, 用于量化場景輻射和傳播中的信號損失;A為環(huán)境光照強度,A(1-t(x))稱為環(huán)境光照, 表示圖像接收的大氣散射光強[10], 是引起色彩偏移和云霧效果的直接原因.
基于統(tǒng)計大量清晰圖像得到的暗通道先驗理論, 是指大部分不含天空優(yōu)質(zhì)圖像的所有像素在R,G,B 3個通道中最少存在一個顏色通道灰度值相當?shù)蜕踔邻吔?, 即在一定的微小區(qū)域內(nèi)最小輻射強度值極低. 對于一副圖像J, 可定義
(2)
其中Ω(x)是中心為像素x的任一微小區(qū)域,y是該區(qū)域內(nèi)任一像素,Jc(y)為像素y在c通道的灰度值. 在不含天空區(qū)域的其他范圍內(nèi), 清晰圖像的Jdark通道灰度值均極小甚至趨近于0, 圖像J的暗通道值即為Jdark, 該統(tǒng)計理論稱為暗通道先驗.
去霧目標是將傳感器接收的有霧圖像利用獲得的有用信息通過去霧還原出清晰圖像. 暗通道先驗圖像去霧方法先根據(jù)暗通道原理獲取先驗知識, 再用有霧圖像退化模型實現(xiàn)圖像去霧.
1.3.1 估算環(huán)境光照強度
暗通道圖像中灰度值越高的區(qū)域云霧越厚, 在輸入圖像中位于該部分區(qū)域的像素灰度值越接近于環(huán)境光照強度. 估算環(huán)境光照強度首先需找到暗通道圖像中灰度值最高并占圖像總像素數(shù)量0.1%的像素點, 記錄其對應的坐標索引, 然后根據(jù)坐標索引在輸入的有霧圖像中找到對應像素點, 計算有霧圖像中對應像素點的灰度平均值作為環(huán)境光照強度A.
1.3.2 估算大氣透射率
(3)
按下式再次運算3個顏色通道的最小值, 可得以像素x為中心的濾波窗口內(nèi)的灰度最小值:
(4)
根據(jù)暗通道先驗原理和式(2)可得
(5)
于是, 由式(4)和式(5)可計算出大氣透射率為
(6)
1.3.3 圖像去霧處理
暗通道先驗條件用于量化云霧厚度和全部像素的輻射還原量, 恢復出優(yōu)質(zhì)清晰的圖像. 通過式(1)的有霧圖像退化模型和環(huán)境光照強度A、 大氣透射率t(x), 可進行單幅圖像去霧處理, 表示為
(7)
通過對暗通道先驗算法的剖析和代碼實驗發(fā)現(xiàn), 計算暗通道圖像和使用引導濾波精細化透射率圖層時運行耗時較多, 下面通過改進大氣透射率圖層精細化處理方法提高去霧運算的速度.
根據(jù)上述暗通道圖像得到的大氣透射率仍不能滿足去霧精度要求, 此時暗通道圖像只能大致展現(xiàn)云霧的厚薄情況. 由于在求取暗通道過程中進行了濾波窗口操作, 若用式(7)直接去霧會產(chǎn)生光暈現(xiàn)象. 在暗通道圖層中, 任意像素暗通道值都是以該像素為中心的n×n窗口區(qū)域內(nèi)所有像素在R,G,B通道的最低灰度值, 若圖層中有一個像素最低灰度值很小, 則以該像素為中心的n×n區(qū)域內(nèi)所有像素暗通道值都將很小, 導致區(qū)域灰度偏低, 形成擴張現(xiàn)象, 因此去霧前需精細化處理大氣透射率圖層.
引導濾波(guided filter, GF)[11-12]的處理效果與雙邊濾波相似, 但算法時間復雜度為O(n). 使用引導濾波的輸出圖像與引導圖像是局部線性關系, 處理圖像的平滑效果和邊緣保護效果更顯著. 引導濾波使輸出圖像與引導圖像具有相似的紋理結(jié)構(gòu)和場景信息[13-14], 在圖像去霧算法中可用于精細化大氣透射率圖層, 處理過程如圖1所示.
圖1 引導濾波示意圖Fig.1 Schematic diagram of guided filtering
線性平移變量濾波(linear translation-variant filtering, LTVF)中包括輸入圖像p、 輸出圖像q和引導圖像I. 輸入圖像與引導圖像是預先設定的兩幅相同或相異圖像, 特定像素i的濾波過程為
(8)
其中i,j表示像素坐標索引,Wij為引導圖像相關的濾波核, 其不依賴于圖像p, 引導濾波器與輸入圖像p線性相關[15-16].
設引導濾波是由輸出圖像q與引導圖像I構(gòu)成的線性相關模型, 表示為
qi=akIi+bk, ?i∈wk,
(9)
其中線性系數(shù)(ak,bk)在以像素k為中心的濾波窗口區(qū)域wk內(nèi)是定值. 由式(9)可知, 輸出圖像q的梯度隨引導圖像I的梯度變化而改變, 引導濾波可將梯度和紋理信息傳遞給輸出圖像, 使兩種圖像的結(jié)構(gòu)和場景信息非常接近, 因此引導濾波可用于精細化大氣透射率圖層、 圖像自動分割等領域. 為計算線性系數(shù)(ak,bk), 假設輸出圖像q是輸入圖像p減去噪聲部分, 表示為
qi=pi-ni,
(10)
式中ni表示噪聲.
在線性模型式(9)的基礎上, 要使輸入圖像p中噪聲占比最小, 需要輸入與輸出圖像間的差異最小. 用期望值E(ak,bk)量化以像素k為中心的濾波窗口區(qū)域內(nèi)輸入圖像p與輸出圖像q之間的差異, 表示為
(11)
其中ε為正則化參數(shù), 用于確定ak在式中所占比率. 如需濾波窗口區(qū)域wk內(nèi)噪聲占比最小, 則求解使式(11)取最小值的線性參數(shù)對(ak,bk)為
(12)
(13)
(14)
將線性系數(shù)(ak,bk)通過式(9)計算得到輸出圖像q, 再用式(8)可得出引導濾波內(nèi)核為
(15)
在圖像去霧過程中, 使用引導濾波對計算出的大氣透射率圖層精細化處理前后對比效果如圖2所示.
圖2 大氣透射率圖層經(jīng)引導濾波處理前后對比Fig.2 Comparison of atmospheric transmittance layers before and after guided filtering
暗通道先驗是在一個微小的局部區(qū)域內(nèi), 需要很低的最小輻射強度[17-18], 暗通道濾波用濾波方式找到其窗口區(qū)域內(nèi)輻射強度的最小值[19-21]. 將暗通道濾波視為二維黑白圖像的最小值濾波向三維空間彩色圖像的擴展, 三維空間指R,G,B顏色空間. 暗通道濾波要找到濾波窗口區(qū)域內(nèi)3個顏色空間中輻射強度最小的像素, 并將該像素的輻射強度作為暗通道值. 為避免光暈現(xiàn)象, 需選取一個適合的濾波器, 在實現(xiàn)圖像平滑處理的同時保護其邊緣信息. 中值濾波器相比引導濾波更易實現(xiàn), 且算法時間復雜度較低, 能節(jié)約圖像去霧時間. 由于兩個直方圖累加是一個O(1)操作, 使得對每增加一個像素, 直方圖只需累加一個固定的數(shù), 從而獲得O(1)的復雜度. 實現(xiàn)算法時間復雜度為O(1)的最小值濾波過程如下.
1) 初始化核直方圖. 若圖像有n列, 每個直方圖對應于圖像中的一列, 則需處理n個直方圖, 這些直方圖包含了以某個像素為中心的(2r+1)個像素信息(r為濾波窗口的半徑). 以每列直方圖第一個像素為中心的(2r+1)個像素信息為圖像初始狀態(tài), 更新處理時, 這些直方圖隨滑塊移動并完成更新.
2) 更新核直方圖. 隨著滑塊移動, 將每列直方圖減去最左側(cè)像素的直方圖信息, 如圖3(A)所示; 同時增加最右側(cè)像素的直方圖信息, 如圖3(B)所示, 每移動一列都是O(1)操作.
圖3 更新核直方圖Fig.3 Update kernel histogram
3) 獲取最小值. 根據(jù)更新的核直方圖取最小值.
通過上述計算, 所有單像素操作都與窗口大小無關, 即算法空間復雜度為O(1). 該算法偽代碼如下.
輸入: 圖像X, 大小為m×n, 濾波窗口半徑為r;
輸出: 與圖像X大小相同的圖像Y;
初始化核直方圖H及每列的直方圖h1,h2,…,hn
Fori=1 tomdo
Forj=1 tondo
RemoveXi-r-1,j+rfromhj+r
AddXi+r,j+rtohj+r
H←H+hj+r-hj+r-1
Yi,j← min(H)
End For
End For.
要將上述灰度圖像最小值濾波后擴展適用于彩色圖像, 可分別對R,G,B 3個波段圖層最小值濾波, 然后比較3個波段圖層最小值濾波后的像素點, 取3個像素值中的最小值作為暗通道值, 即實現(xiàn)了與濾波窗口無關, 算法時間復雜度為O(1)的暗通道濾波.
使用該算法時間復雜度為O(1)的中值濾波器雖然對圖像邊緣保護效果并不突出, 但其是極簡單且具有邊緣保護效果的濾波器. 實驗表明, 在圖中云霧分布較均勻時, 使用中值濾波器代替引導濾波效果較好; 當云霧分布不均勻時, 在暗通道先驗去霧中使用中值濾波器的效果欠佳, 采用引導濾波. 基于暗通道先驗的圖像去霧改進算法流程如圖4所示.
圖4 暗通道先驗去霧算法速度改進流程Fig.4 Flow chart of speed improvement of dark channel prior fog removal algorithm
當圖像中云霧均勻分布時, 用中值濾波代替引導濾波精細化大氣透射率圖層, 以提高精細化處理速度. 實驗在湖北省恩施市航拍圖像中選取68幅云霧較均勻的圖像作為原始數(shù)據(jù)進行驗證. 圖片均為JPEG格式, 包含了R,G,B 3個顏色通道. 實驗使用的計算機配置為: CoreTMi7-5500U @2.4 GHz處理器, 16 GB內(nèi)存, Windows10 64位操作系統(tǒng). 為使算法改進前后的去霧效果偏差較小, 在不影響輸出圖像后期處理的基礎上盡量提高計算速度, 實驗設定了統(tǒng)一參數(shù)批量處理全部原始有霧圖像.
實驗使用無參考的客觀圖像評價方法, 將圖像信息熵、 邊緣強度、 方差作為評價指標. 以輸入如圖5所示的云霧均勻圖像為例, 文件大小7.3 MB, 使用中值濾波代替引導濾波在暗通道先驗去霧算法中精細化處理大氣透射率圖層.
圖5 實驗選取的有霧圖像Fig.5 Foggy image selected for experiment
圖6為分別使用中值濾波和引導濾波精細化的大氣透射率圖層, 其中光暈現(xiàn)象主要由計算暗通道時的塊狀處理結(jié)果導致. 由圖6可見, 使用中值濾波精細化大氣透射率圖層效果不如引導濾波, 但中值濾波仍很好地模糊了因計算暗通道圖層而產(chǎn)生的塊狀處理結(jié)果, 并在一定程度保護了圖像邊緣信息.
圖6 使用中值濾波和引導濾波精細化后的透射率圖層Fig.6 Transmittance layer refined by median filtering and guided filtering
用圖6的兩幅大氣透射率圖層分別進行暗通道先驗去霧處理, 還原出的清晰圖像如圖7所示, 雖然使用中值濾波與引導濾波去霧后的兩幅清晰圖像色調(diào)存在一定差異, 但仍還原出了較清晰的圖像.
圖7 使用中值濾波和引導濾波去霧后的清晰圖像對比Fig.7 Comparison of clear images after defogging using median filtering and guided filtering
實驗選取68幅云霧較均勻的原始航拍有霧圖像(大小為5 180×5 400像素)進行暗通道實驗, 其中部分原始有霧圖像如圖8所示. 將圖8原始圖像對應的大氣透射率圖層使用中值濾波精細化, 得到如圖9所示優(yōu)化后的大氣透射率圖層.
圖8 原始航拍有霧圖像Fig.8 Original aerial foggy images
圖9 中值濾波處理后對應的透射率圖層Fig.9 Corresponding transmittance layer after median filtering processing
用中值濾波將選取的68幅云霧均勻圖像經(jīng)過暗通道先驗去霧后, 對實驗輸出的全部清晰圖像進行客觀圖像質(zhì)量評價, 結(jié)果列于表1. 由表1可見, 當輸入圖像云霧均勻時, 用中值濾波代替引導濾波對大氣透射率圖層精細化后的去霧效果有一定影響, 但仍能還原出清晰圖像. 當去霧圖像應用不受影響時, 用中值濾波代替引導濾波精細化處理大氣透射率圖層顯著提高了暗通道先驗去霧算法的運算速度.
表1 精細化大氣透射率模塊改進速度評價指標Table 1 Evaluation indexes of refined atmospheric transmittance module for improving speed
綜上所述, 本文針對無人機航拍圖像中云霧情況多樣、 數(shù)據(jù)量大等特點, 通過計算機編碼實驗, 提出了一種簡單且具有邊緣保護效果、 與濾波窗口無關、 時間復雜度為O(1)的中值濾波器, 采用中值濾波代替引導濾波精細化處理大氣透射率圖層, 利用無參考的客觀圖像質(zhì)量評價方法量化評價算法改進前后的去霧圖像效果. 結(jié)果表明, 當輸入圖像中云霧較均勻時, 中值濾波更易實現(xiàn), 且算法時間復雜度較低, 能節(jié)約去霧時間, 實現(xiàn)圖像平滑處理. 用中值濾波相比引導濾波去霧后的圖像效果有影響, 但仍可還原出清晰圖像, 在不影響應用時, 改進算法的計算速度明顯提高. 由于中值濾波僅很好地模糊了計算暗通道圖層產(chǎn)生的塊狀處理結(jié)果, 并未真正加入圖像的梯度和紋理信息, 所以中值濾波對云霧分布不均的圖像處理效果不佳, 精細化處理仍需使用引導濾波.