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

王振國，李少敏，曲廣兵

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

0 引言

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

1 高動態范圍圖像的顯示

目前，HDR輻射分辨率的數字圖像越來越受到計算機圖像處理和應用領域的青睞，應用范圍越來越廣。可以預測，隨著現代數字圖像技術的日益發展，未來的數字相機就能夠直接獲取HDR圖像和視頻［3－4］。

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

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

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

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

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

2 梯度域動態壓縮

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

需要補充說明的是，本文算法的實現和計算過程都是針對圖像亮度值的對數來完成的，而非直接針對亮度值本身完成。這里采用圖像亮度值的對數來進行計算的原因主要為:①亮度的對數能夠較好地對感知到的光線亮度進行近似;②對數域中的梯度能夠對應于亮度域的局部亮度。

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

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

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

現在可以通過對壓縮后的導數求積分來構建一個簡化的動態域函數I(另有一個附加常數C)，即

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

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

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

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

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

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

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

最小化，式中

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

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

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

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

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

3 梯度衰減函數

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

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

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

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

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

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

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

4 算法應用

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

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

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

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

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

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

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

5 實驗結果與分析

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

表1 圖3中不同結果圖像顯示質量評價參數Tab．1 Evaluation parameter list of displaying quality for different mages in Fig．3

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

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

6 結論

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

2)經過有針對性的實驗表明，本文算法能夠實現對高動態范圍圖像的完美顯示，同時還能夠很好地保留圖像的細節信息和抑制圖像中存在的邊緣效應。

［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］ 朱述龍，朱寶山，王紅衛．遙感圖像處理與應用［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］ 關元秀，程曉陽．高分辨率衛星影像處理指南［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．