基于梯度域的高輻射分辨率遙感圖像動態(tài)范圍壓縮算法

王振國，李少敏，曲廣兵

(中國人民解放軍96635部隊，天津 300350)

0 引言

現(xiàn)代遙感技術(shù)的日益發(fā)展，促使遙感信息獲取平臺和成像系統(tǒng)的水平不斷提高，特別是遙感圖像的空間分辨率、輻射分辨率、光譜分辨率和時間分辨率都在不斷提高，使得高質(zhì)量遙感數(shù)據(jù)的應用已經(jīng)從單純的定性應用向定量應用的方向發(fā)展［1］。目前，人們常用的航天航空數(shù)字傳感器(如IKONOS，ALOS，IRS，Quick Bird以及ADS40數(shù)字航空相機等)所獲取的數(shù)字遙感圖像的輻射分辨率都在11 bit或者14 bit以上［2］。利用高輻射分辨率遙感圖像(high radiation resolution remote sensing image，HRRRSI)，人們可以更加精細和準確地得到各類地物的細節(jié)結(jié)構(gòu)和光譜信息，從而可有效地降低混合像元對定量遙感分析結(jié)果的影響，進一步提高定量遙感分析的精度，為定量遙感分析的廣泛應用奠定基礎。然而，與一般的8 bit遙感圖像相比，HRRRSI圖像的利用在許多方面都存在特殊之處，其中最直觀的也是最常見的問題就是對HRRRSI的顯示問題?，F(xiàn)在常用的做法是將HRRRSI進行壓縮后顯示，常用的高動態(tài)域(high dynamic range，HDR)輻射圖像的顯示方法分為全局映射和空間變化2類。人們將前者作為色調(diào)再現(xiàn)曲線(tone reproduction curves，TRC)，將后者作為色調(diào)再現(xiàn)算子(tone reproduction operators，TRO)。對于多數(shù)的TRC方法，都是將HDR值進行線性分級，進而將圖像的像元值映射到設定的范圍之中(比如常見的圖像取值范圍［0，1］);現(xiàn)在常見的TRC方法主要有亮度校正和直方圖均衡化等，該方法可以較好地保留原始圖像像元之間的比例關(guān)系，但會降低所顯示圖像中比原始圖像動態(tài)范圍小的影像清晰度。針對上述問題，本文就HRRRSI顯示問題，提出了一種基于梯度域動態(tài)范圍壓縮的圖像增強和顯示算法，能夠明顯改善遙感圖像的顯示效果。通過與常用的線性拉伸處理顯示效果對比，進一步驗證了本文算法的可行性和有效性。

1 高動態(tài)范圍圖像的顯示

目前，HDR輻射分辨率的數(shù)字圖像越來越受到計算機圖像處理和應用領域的青睞，應用范圍越來越廣?？梢灶A測，隨著現(xiàn)代數(shù)字圖像技術(shù)的日益發(fā)展，未來的數(shù)字相機就能夠直接獲取HDR圖像和視頻［3－4］。

與常用的低動態(tài)域圖像相比，HDR圖像擁有諸多優(yōu)點［5］，且已有的應用都證明了HDR圖像具有廣闊的應用前景［6－7］。不過，HDR 圖像在應用中也遇到一定問題，就是在不同類型的普通顯示器上顯示的HDR圖像的動態(tài)域要比真實場景的動態(tài)域小很多;所以，針對現(xiàn)在普遍采用的低動態(tài)域顯示設備，就需要對其顯示HDR圖像的方法進行深入研究。

圖1示出輻射分辨率為14 bit的高輻射分辨率航空遙感圖像在低動態(tài)域顯示設備上的顯示結(jié)果。

圖1 輻射分辨率為14 bit的航空遙感圖像顯示結(jié)果Fig．1 Displaying results for airborne remote sensing image with 14 bit radiation resolution

從圖1中可以看出HDR圖像在普通顯示設備上進行顯示時出現(xiàn)的問題:在直接顯示的原始圖像(圖1(a))中，影像細節(jié)信息沒有得到很好的展現(xiàn);而經(jīng)過Photoshop軟件拉伸處理后的顯示結(jié)果(圖1(b))雖然比圖1(a)有了較大的改進，但所顯示的影像細節(jié)信息也不是特別完美。由此可見，在一般的顯示設備上對14 bit的HDR圖像進行顯示會存在一定的問題，其原因主要是一般顯示設備的動態(tài)范圍通常低于100∶1。

因此，本文針對HRRRSI顯示問題，提出了一種對HDR壓縮圖像進行顯示的增強算法，以便能夠在一般的顯示設備上對HDR圖像進行完美顯示。實驗結(jié)果表明，該算法理論簡單、運算效率高、效果較好。

2 梯度域動態(tài)壓縮

本文算法是建立在一個被人們所廣泛接受的假設［8］的基礎上，即人類的視覺系統(tǒng)對視網(wǎng)膜所接受到的亮度變化并不敏感，但對對比度的局部變化比較敏感。該算法還基于一個必要的前提，即對于HDR圖像而言，其亮度的顯著變化一定會在某些尺度上產(chǎn)生較大的亮度梯度，而細節(jié)信息則通常都對應比較小的亮度梯度。故本文算法的基本思路就是通過對不同尺度上較大亮度的識別，在保持其方向不變的情況下減小其量值;同時保證對較大尺度的衰減程度要大于對較小尺度的衰減程度，這樣，就可以在壓縮亮度變化范圍的同時，對圖像的細節(jié)信息給予保留。同時，對經(jīng)過壓縮后的HDR圖像還可以在衰減的梯度域中進行重建。

需要補充說明的是，本文算法的實現(xiàn)和計算過程都是針對圖像亮度值的對數(shù)來完成的，而非直接針對亮度值本身完成。這里采用圖像亮度值的對數(shù)來進行計算的原因主要為:①亮度的對數(shù)能夠較好地對感知到的光線亮度進行近似;②對數(shù)域中的梯度能夠?qū)诹炼扔虻木植苛炼取?/p>

為了對本文算法進行說明和解釋，首先對一維情況下的算法原理進行簡述。

給定一個高動態(tài)范圍的一維函數(shù)，用H(x)表示該函數(shù)的對數(shù)。如上所述，本文算法的目標是在保持較小局部變換的同時，壓縮H中較大的數(shù)值變化。此目標是通過把一個空變衰減映射函數(shù)Ф應用到導數(shù)H'(x)來完成的，即

式中:G為衰減后的導數(shù);H'為原始導數(shù);Ф為空變衰減映射函數(shù)。G的符號與H'相同，但H'的大小被一個由Ф所決定的因子改變;Ф設計用于更多地衰減較大的導數(shù)。事實上，如下面所解釋的那樣，Ф在不同尺度上代表不同的導數(shù)值。

現(xiàn)在可以通過對壓縮后的導數(shù)求積分來構(gòu)建一個簡化的動態(tài)域函數(shù)I(另有一個附加常數(shù)C)，即

最后，對上式的結(jié)果進行冪運算就可得到亮度函數(shù)。整個過程如圖2所示。

圖2 一維掃描線衰減前后導數(shù)變化對比Fig．2 Comparison between derivative changes before and after attenuation of one dimensional scanning line

圖2(a)為一條動態(tài)范圍為2 415∶1的高動態(tài)域一維掃描線(scanline);圖2(b)為圖2(a)的對數(shù)，即H(x)=lg(scanline);圖2(c)為圖2(b)的原始導數(shù);圖2(d)為衰減后的導數(shù)G(x);圖2(e)為重建的信號I(x)(如式(2)中所定義);圖2(f)為一條低動態(tài)域掃描線exp(I(x))，其新的動態(tài)范圍為7．5∶1。從圖2(a)和圖2(f)中2條掃描線的動態(tài)范圍對比中可以看出，后者的動態(tài)范圍遠遠小于前者，這樣更便于對掃描線的對比顯示;對應到二維圖像上，則動態(tài)范圍較低的圖像更容易在低動態(tài)域設備上得到較好的顯示效果。應該注意的是，每個示意圖中垂直軸的尺度并不一樣，只有圖2(c)和圖2(d)中的垂直軸的尺度一致，這是為了更好地對比衰減前后導數(shù)的變化。

為了將上述方法擴展到二維HDR函數(shù)中H(x，y)，本文采用梯度▽H代替導數(shù)。為了避免將空間失真引入到圖像中，只改變梯度的量值，而不改變其方向。因此，與一維情況相似，可以計算

與一維情況不同的是，不能簡單地通過對G求積分得到一個壓縮后的動態(tài)域圖像，因為該函數(shù)不一定是可積的;也就是說，可能并不存在使G=▽I的圖像I。事實上，勢函數(shù)(potential function)的梯度(以二維為例)必須是一個守恒場［9］。換言之，梯度▽I=(?I/?x，?I/?y)必須滿足

但是，對于G而言，這種情況很少見。

解決上述問題的一個可能的辦法是將G垂直投影到一個有限的、能夠生成可積矢量場的規(guī)范正交基函數(shù)集上(如傅立葉基函數(shù))［10］。在本文算法中，使用了一個更直接、也更有效的方法，即尋找在最小二乘框架下梯度最接近G的函數(shù)I所有的二維勢函數(shù)空間。換言之，I應該能夠使得積分

最小化，式中

根據(jù)變分原理，使積分最小化的函數(shù)I滿足歐拉－拉格朗日方程，即

這是一個I中的偏微分方程。替換F可以得到

公式(8)兩邊除以2并化簡，可得到泊松方程

式中:Δ2為拉普拉斯算子，即

div G為向量場G的散度，定義為

3 梯度衰減函數(shù)

本文算法利用Φ(x，y)的一個因子在各個像元上對HDR的圖像梯度進行衰減，從而實現(xiàn)HDR壓縮;并希望衰減能夠依次進行，且較大梯度的收縮比小梯度的收縮更大。

真實圖像中包含所有尺度的邊緣。因此，為了可靠地探測所有重要的強度變化，必須采用多分辨率邊緣探測的方法。然而，不能簡單地在探測到的分辨率圖像上對各個梯度進行衰減，這可能會在銳利邊緣周圍產(chǎn)生光環(huán)效應(如第2節(jié)所提到的那樣)。解決的方法是將合適的衰減從被探測的那個尺度擴展到全分辨率圖像。因此，所有梯度處理都在單分辨率尺度上進行，這樣就不會產(chǎn)生光環(huán)效應。

首先，需要構(gòu)造一個高斯金字塔 H0，H1，…，Hd，其中:H0是全分辨率HDR圖像;Hd是金字塔中的最低一級圖像。d的選取標準為使Hd的寬度和高度不小于32。在第k層上利用中心差分公式來計算梯度，即

第k級上的一個尺度因數(shù)φk(x，y)由此尺度上的像元梯度值所決定，即

它是帶有2個參數(shù)的函數(shù)族。第一個參數(shù)α決定了哪個梯度值沒有改變。較大的梯度被衰減(假定β＜1)，而小于α的梯度值則被稍微放大。本文設置α為平均梯度值的0．1倍，β取值在0．8～0．9之間。

全分辨率梯度衰減函數(shù)Φ(x，y)采用自上而下的方式進行計算，計算方法是用線性內(nèi)插將比例因數(shù)φk(x，y)由某一級擴展到下一級，并利用逐點相乘法將其累積。這個過程可通過求解方程組E，即

來實現(xiàn)。式(14)中:d為最低一層;Фk為在第k級所累積的衰減函數(shù);L為線性內(nèi)插中的一個采樣算子。計算結(jié)果中，在最佳層次上每個像元的梯度衰減是由所有通過該像元的邊緣(來自不同級別)強度來決定的。

4 算法應用

為了求解式(9)中的微分方程，必須首先確定邊界條件。在本文的例子中，最直接的選擇是諾伊曼(Neumann)邊界條件，即▽I·n=0(垂直于邊界方向的導數(shù)為0)。

由于拉普拉斯算子Δ2和div都是線性算子，所以利用標準有限差使其近似服從一個線性方程系統(tǒng)。具體地說，可以利用近似方程

在全分辨率圖像中，將像元柵格間距定為1;而對于 梯度▽H，可使用前向差分來近似得到，即

對于div G，則可使用后向差分來近似得到，即

這樣，就確保了前向和后向差分的組合使用對div G的近似與使用拉普拉斯的中心差分方法是一致的。

在邊界處可采用同樣的定義，但要假設原始圖像柵格周圍的導數(shù)為0。例如，對于圖像左邊邊界的每個像元有方程成立。

有限差分方案可產(chǎn)生一個大的線性方程系統(tǒng)——每個方程對應于圖像中的一個像元，但相應的矩陣在每行只有5個非零元素，因為每個像元只和它周圍4個像元相鄰。本文通過全多格網(wǎng)算法(full multi－grid algorithm)和高斯－塞德爾平滑迭代法來解算該方程系統(tǒng)。這樣，要得到一個近似的解，需要大約O(n)步運算(n為圖像中的像元個數(shù))。另外一個選擇是使用“快速泊松解法”，它是用快速傅立葉變換來對拉普拉斯算子進行求逆;然而，這個算法的復雜度可能是O(nlgn)步運算。

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

為了驗證本文提出的基于梯度域的高動態(tài)范圍壓縮算法的效果，本文通過對一景由ADS40數(shù)字航空相機攝取的高輻射分辨率遙感圖像進行試驗和分析，圖像大小為512像元×512像元，圖像的輻射分辨率為14 bit。試驗結(jié)果如表1、圖1和圖3所示。

表1 圖3中不同結(jié)果圖像顯示質(zhì)量評價參數(shù)Tab．1 Evaluation parameter list of displaying quality for different mages in Fig．3

圖3 輻射分辨率為14 bit的航空遙感圖像本文方法顯示結(jié)果Fig．3 Displaying result using the method in this paper for airborne remote sensing image with 14 bit radiation resolution

比較圖1和圖3可以看出，在直接顯示的原始圖像(圖1(a))中幾乎看不到影像細節(jié)信息;在利用Photoshop軟件中的線性拉伸功能處理過的圖像(圖1(b))中雖然可以看出部分影像細節(jié)信息，但陰影部分的信息卻顯示不出來;而在利用本文算法處理后的圖像(圖3)中則可以清楚地看到影像信息和陰影部分的細節(jié)結(jié)構(gòu)信息(如油罐和灌木林陰影范圍內(nèi)的細節(jié)結(jié)構(gòu))。表1中列出的不同結(jié)果圖像顯示質(zhì)量評價參數(shù)也驗證了本文算法的有效性。

6 結(jié)論

1)本文針對高輻射分辨率遙感圖像在普通顯示設備上的顯示問題，通過對高輻射分辨率的特性分析，提出了對其圖像顯示結(jié)果進行增強的有效算法。該算法通過對高輻射分辨率圖像中大梯度域的梯度衰減來實現(xiàn)，能夠在常用的顯示設備上對高輻射分辨率遙感圖像的細節(jié)信息進行精確顯示。

2)經(jīng)過有針對性的實驗表明，本文算法能夠?qū)崿F(xiàn)對高動態(tài)范圍圖像的完美顯示，同時還能夠很好地保留圖像的細節(jié)信息和抑制圖像中存在的邊緣效應。

［1］ Mugnier L M，Robert C，Conan JM，et al．Myopic deconvolution from wave－front sensing［J］．Journal of Optics Engineering，2001，18(4):862－872．

［2］ 朱述龍，朱寶山，王紅衛(wèi)．遙感圖像處理與應用［M］．北京:科學出版社，2006．Zhu SL，Zhu B S，Wang H W．The processing and application of remote sensing images［M］．Beijing:Science Press，2006．

［3］ Schechner Y Y，Nayar S K．Generalized mosaicing［C］//Proceedings of the 8th IEEE International Conferences on Computer Vision．Vancouver，BC:IEEE，2001，I:17－24．

［4］ Nayar SK，Mitsunaga T．High dynamic range imaging:Spatially varying pixel exposures［C］//Proceedings of IEEE Conference on Computer Vision and Pattern Recognition．Hilton Head Island，SC:IEEE，2000:472－479．

［5］ 關(guān)元秀，程曉陽．高分辨率衛(wèi)星影像處理指南［M］．北京:科學出版社，2008．Guan Y S，Cheng X Y．The processing guide for high resolution satellite images［M］．Beijing:Science Press，2008．

［6］ Debevec P．Rendering synthetic objects into real scenes:Bridging traditional and image－based graphics with global illumination and high dynamic range photography［C］//In Proc．ACM SIGGRAPH 98，MCohen，Ed，1998:189－198．．

［7］ Cohen J，Tchou C，Hawkins T，et al．Real－time high dynamic range texture mapping［C］//Gortler S J，Myszkowski K．Eds Springer－Verlag，in Rendering Techniques，2001:313－320．

［8］ Dicarlo JM，Wan dell B A．Rendering high dynamic range images［C］//In Proceedings of the SPIE:Image Sensors，2001，3965:392－401．

［9］ Harris JW，Stocker H．Handbook of mathematics and computational science［M］．New York:Springer－Verlag，1998．

［10］ Frankot R T，Chellappa R．A method for enforcing integrability in shape from shading algorithms［C］//IEEE Transactions on Pattern Analysis and Machine Intelligence，1988，10(4):439－451．