一種改進的遙感影像地形校正方法

臧　熹，楊　博，齊建偉，向夏蕓

(1. 武漢大學，湖北 武漢 430072； 2. 中國國土資源航空物探遙感中心，北京 100083)

An Improved Topographic Correction Method of Remote-sensing Images

ZANG Xi，YANG Bo，QI Jianwei，XIANG Xiayun

?

一種改進的遙感影像地形校正方法

臧熹1，楊博1，齊建偉2，向夏蕓1

(1. 武漢大學，湖北 武漢 430072； 2. 中國國土資源航空物探遙感中心，北京 100083)

An Improved Topographic Correction Method of Remote-sensing Images

ZANG Xi，YANG Bo，QI Jianwei，XIANG Xiayun

摘要：地形校正是遙感影像定量化應用環節之一，以往的地形校正研究多是針對一景影像中很小的局部影像塊來進行處理研究的，對整景大場景影像進行地形校正的研究尚不多?；诖?本文利用高分一號的寬視場相機拍攝的16 m分辨率的遙感影像,研究了大場景下地形校正方法,對C校正模型進行了改進，在C校正模型中加入了反射角的影響，并且驗證了改進模型的合理性；最后對改進的模型與余弦校正模型、傳統的C校正模型的處理結果進行了比較。通過分析，利用改進的模型,影像的標準差普遍變小,影像校正后的陰陽坡亮度值趨于一致的趨勢更明顯。試驗結果表明，對大場景、非星下點成像的遙感影像利用改進模型進行地形校正效果明顯增強。

關鍵詞：地形校正；寬幅；C校正模型；反射角

一、引言

衛星傳感器獲得的影像由于受到地形起伏即坡度和坡向變異的影響而導致陰陽坡影像亮度值的差異，尤其在丘陵地帶和山區，地形坡度、坡向和太陽光照幾何條件等對遙感影像的亮度值的影響是非常顯著的，朝向太陽的坡面會接收到更多的光照，在遙感影像上色彩自然要亮一些，背向太陽的陰面由于反射的是天空散射光和鄰近地物的散射光，在圖像上表現得要暗淡一些，會出現同物異譜或同譜異物的現象。復雜地形地區遙感影像的這種輻射畸變稱為地形效應，這種現象嚴重干擾了影像中地物的光譜信息，從而導致這些地區數字遙感影像自動分類的錯分、難分現象，嚴重影響著數字遙感影像的信息提取精度和使用[1]。因此必須對地形效應進行校正處理，即地形校正。目前的地形校正主要是基于太陽入射角的地形校正模型，其主要目的是消除影像亮度值與太陽入射角的關系，將影像斜面上的亮度值投影到水平面上來，從而達到地形校正的目的。

現有的地形校正模型很多，其中C校正模型因其簡單的使用條件和輸入參數及較好的普適性，成為應用最廣泛的校正模型之一。之前的研究對象大都是選取一景影像的一小部分影像塊，此時，影像覆蓋范圍比較小，像元數量比較少，像元的反射角對C校正模型經驗系數的擬合可能沒有太大的影響，而且可以不用考慮像元的反射角在影像塊內部的差異。但是當遙感影像的幅寬比較寬時，由于覆蓋區域廣，一景影像的像元數量非常多，在成像范圍內，像元的反射角及像元在景內位置的差異引起的反射角的變化對整景影像的校正方程的相關性及擬合經驗系數的計算就可能產生累積性和趨勢性的影響。尤其是當衛星不是星下點成像時，反射角的影響可能會更明顯。對于大場景遙感影像的地形C校正，考慮到反射角可能造成的影響，本文在C校正模型的基礎上加入了反射角，對C校正模型進行了改進；利用高分一號的寬視場相機拍攝的16 m分辨率的200 km幅寬的遙感影像對改進的方法進行了驗證，證明了反射角確實對校正結果有著比較大的影響，并且將校正結果與余弦模型和傳統的C模型進行了比較，發現校正后影像的方差有所降低，陰陽坡亮度值趨于一致的趨勢更加明顯。

二、改進的地形校正模型

1. C校正模型

C校正模型[2]是朗伯體模型，其假定地表具有朗伯體反射特性，即地表是一理想散射反射體，入射太陽輻射在各個方向和任意波長具有均等反射可能性。除了入射太陽輻射，地表像元接收到的輻射還包括鄰近地形的散射輻射和天空漫散射輻射，因此C校正中引入參數C用來描述大氣和鄰近地形對像元的散射輻射的貢獻，從而減弱余弦校正中的過校正問題。研究表明，對于任意波段影像的像素DN值和其對應的太陽入射角余弦值都遵循線性關系[3]。

在斜面上，像素亮度值與光照系數(太陽入射角余弦值)滿足關系式

LT=a+b·cosi

(1)

式中，LT為校正前影像斜面的亮度值；i為斜面的太陽入射角；a和b是通過LT與cosi的關系擬合出的直線的截距和斜率。太陽入射角i與像元的地理經緯度、太陽高度角和方位角，以及像元的坡度和坡向有關系，計算公式為

cosi=cosθSZcosθS+sinθSZsinθScos(φS-φA)

(2)

式中，i為太陽入射角；θSZ為太陽天頂角；θS為像元坡度角；φS為太陽方位角；φA為像元坡向角。

對于水平面，太陽入射角就是太陽天頂角，其像素亮度值與光照系數的關系式為

LH=a+b·cosθSZ

(3)

式中，LH為水平面的像元亮度值。

將傾斜地面的直線投影到水平地面的直線，即

即可得到C校正模型的方程

(4)

式中，參數C是根據像元亮度值與光照系數之間的統計相關性，由線性回歸擬合得到的斜率b和截距a可計算得到C=a/b。

參數C并沒有明確的物理含義，但根據像元受到的散射輻射隨著擬合出來的截距a的增大而增大，通?？梢姽獠ǘ伪燃t外波段的C值要大一些，因為可見光波段的波段短，受到大氣和鄰近地形的散射輻射更嚴重[4]。C校正模型簡單易行，并且在一定程度上避免了低光照區域的過校正，可以對遙感影像有比較好的地形校正效果。

可以注意到，式(1)中表達的是地表像元在地面處接收到的輻射。但是地表像元在地面處接收到的輻射還需要經過反射才能到達傳感器，因此式(1)并不能很好地對傳感器接收到的輻射進行表示。當影像幅寬比較寬時，由于影像覆蓋范圍比較廣，景內像元的最左端和最右端的反射角的大小就會有比較大的差異。而當衛星不是星下點成像時，傳感器和地表像元連線與地表像元天頂線形成的反射角不再為0，這對景內像元的反射角的影響也是不可忽視的。因此有必要考慮反射角，從而更好地模擬到達傳感器處的輻射量。

2. 加入反射角的改進模型

(1) 加入反射角的模型原理

傳感器接收到的輻射包括以下幾部分：地表像元反射的太陽直射輻射、鄰近地形的散射輻射、天空漫散射輻射及程輻射[4]。由于是相對地形校正，校正的目的主要是消除影像亮度值與太陽入射角的關系，即校正太陽直射輻射的部分，而鄰近地形的散射輻射、天空漫散射輻射及程輻射這3個量對于傳感器接收到的總輻射來說是加性輻射量，因此統一將這3個輻射量作為一個變量a來表示。成像時刻太陽直射輻射在大氣下層的輻射量為b；考慮到成像時刻的太陽天頂角，則到達地表像元的太陽直接輻射量為b·cosi；考慮到像元的反射角，基于地表反射是朗伯體的假定，地表像元接收到的輻射在像元與傳感器連線這一反射方向上，實際被地表像元反射的太陽輻射則表示為b·cosi·cose。

因此，在斜面上，傳感器接收到的總輻射表示為

LT=a+b·cosi·coset

(5)

式中，LT為校正前影像斜面的亮度值；i為斜面像元的太陽入射角；et為斜面像元的反射角；a和b的求解與C校正模型中的類似，是通過LT與cosi·coset的關系，擬合出的直線的截距和斜率。

類似的，在水平面上，傳感器接收到的總輻射表示為

LH=a+b·cosθSZ·coseh

(6)

式中，LH表示水平面像元的亮度值；θSZ表示太陽天頂角；eh表示水平面的像元反射角。

將傾斜地面的直線投影到水平地面，即

即可得到以下C校正模型的方程

(7)

式中，參數C=a/b。

(2) 反射角的計算

① 平面像元反射角的計算

當衛星不是星下點成像時，衛星和地面目標的關系如圖1所示[5]。

圖1中，α是地面目標的天頂方向與地面目標和衛星連線方向所成的夾角，稱為衛星天頂角；β是衛星的星下點(即天底)至地面目標點的張角，稱為衛星星下點角。衛星天頂角與衛星的星下點角存在的幾何關系，可由正弦定理表示為

(8)

式中，H為衛星飛行的平均軌道高度；R為地球半徑。

由圖1可知，衛星天頂角α就是地面目標像元的反射角，因此有

(9)

圖1　衛星天頂角與地面目標星下點角的關系圖

對于線陣推掃的傳感器來說，傳感器是垂直于衛星軌道進行推掃成像，因此在計算景內像元的星下點角時，近似認為影像上每一掃描行內的像元的星下點角的大小不同，且與其在景內的位置有關；而影像上掃描行間的相同列的像元的星下點角相同。這樣只需計算通過中心像元這一掃描行上所有像元的星下點角即可得到一景內任意像元的星下點角，進一步可以得到一景內任意像元的反射角。

根據景內地面目標像元的位置差異計算每個像元的星下點角β，景內任意像元的星下點角與中心像元的星下點角的關系如圖2所示。

圖2中，O點表示景中心像元的位置；A表示與中心像元O同一掃描行內的任意像元A的位置；

M點是星下點；β0表示景中心像元的星下點角；β′表示像元A的星下點角。

圖2　任意像元星下點角與中心像元星下點角的關系圖

景中心像元的衛星天頂角可以從影像的元數據中讀出，記為α0，由式(7)可以計算出景中心像元的星下點角β0，進而可以近似計算出OM的距離

OM=H·tanβ0

(10)

又由像元A與中心像元O的位置關系，可以計算AO的距離，從而可以近似得到AM=AO+OM，因此像元A的星下點角β′為

(11)

再由式(8)可得任意像元A的衛星天頂角，即反射角α′。

② 斜面像元反射角的計算

圖3是太陽—地面—衛星的幾何關系圖[6]。

圖3　太陽—地面—衛星的幾何關系圖

根據圖3可以給出以下近似計算斜面像元反射角的關系：

圖3(a)表示當地面是平面時，太陽入射角為i，θsz是太陽天頂角，衛星天頂角為α，平面像元的反射角e=α。

設衛星的方位線與坡向線所成的不大于180°的夾角為δ，圖中θs表示斜面的坡度。

圖3(b)表示當地面是斜面，且90°<δ≤180°時，斜面的坡向背離衛星，斜面像元的反射角e=α+θs。

圖3(c)表示當地面是斜面，且0°≤δ<90°時，斜面的坡向朝向衛星，斜面像元的反射角e=α-θs。

當δ=90°時，斜面像元的反射角e=α。

綜合圖3(b)和圖3(c)可以看出，當衛星是星下點成像時，即α=0時，水平面像元的反射角為0，斜面像元的反射角等于像元的坡度。

③ 像元平面反射角的取值范圍

當遙感影像的幅寬比較寬時，像元在景內位置的差異引起的反射角的變化不容忽視。高分一號衛星的寬視場成像儀中的單相機拍攝的一景影像幅寬為200 km，平均飛行軌道高度為645 km。根據高分一號的寬視場數據的基本信息來估算一景影像內平面反射角的取值范圍。

a. 當衛星是星下點成像時，如圖4所示。

圖4　衛星是星下點成像時，像元反射角示意圖

通過計算可得，α的取值范圍為0～8.8°，可見雖然衛星是星下點成像，但當影像幅寬比較寬時，一景影像內平面反射角的變化還是很明顯的。

b. 當衛星不是星下點成像時，反射角的變化可能會更明顯。高分一號衛星的寬視場的成像儀由4臺相機并排組成，兩兩相機的主光軸之間有16°的夾角。因此在成像時并不是星下點成像的，其景中心像元的衛星天頂角α0最小可以達到0，最大可以達到30°以上。根據高分一號寬視場相機的特點，估算一景影像內平面反射角的取值范圍，根據式(9)—(11)有：當α0=10°時，一景影像的像元平面反射角α的取值范圍為0.3°～19.3°；當α0=20°時，一景影像的像元平面反射角α的取值范圍為10.8°～28.5°；當α0=30°時，一景影像的像元平面反射角α的取值范圍為21.6°～37.6°。

可以看出，當衛星不是星下點成像時，一景影像內的像元平面反射角的變化更加劇烈。因此，當對寬幅遙感影像進行地形校正，尤其是衛星不是星下點成像時，有必要將像元的反射角的影響考慮進去。

三、試驗驗證與分析

1. 試驗數據

本文中所用的3景試驗數據橫跨河北、河南、山西和陜西4省，景序列號分別為153580、153595和153610，是高分一號衛星寬視場成像儀于2013年11月3日11:26(北京時間)拍攝的，其星下點地面像元分辨率優于16 m。試驗所用數據是L1A級別，經過系統的輻射校正，有4個波段(紅、綠、藍、近紅外)，每景影像尺寸為12 000像素×13 400像素，其基本信息見表1。試驗中所用的全球DEM的空間分辨率為1″(約30 m)，垂直精度20 m，水平精度30 m。

2. 試驗過程

1) 數據準備。根據每景影像的經緯度范圍，選取同區域的全球DEM數據切片，對DEM數據進行無縫鑲嵌，從而可以完全覆蓋影像范圍。根據影像的經緯度范圍內插出每個像素的經緯度，選用的內插方法是雙線性內插。通過影像覆蓋區域的DEM數據，根據像素的經緯度，計算像素在DEM中的格網的行列號，得到像素的坡度和坡向。從影像元數據中可以獲得成像時的太陽天頂角和太陽方位角，進一步根據式(2)計算出太陽入射角。最后根據第二章中介紹的斜面像元反射角的計算方法，計算每個像素的斜面反射角。

2) 模型驗證。用景號為153580、153595和153610的3景影像，驗證本文加入反射角的改進地形校正模型，直線擬合的相關系數有所提高，證明改進模型的合理性，從而可以進一步應用改進的模型對影像進行校正。

3) 校正模型。針對景號為153610的影像分別采用余弦校正[7]、C校正和本文的改進方法進行校正，比較試驗結果并進行評價。

表1　試驗數據基本信息

3. 試驗驗證與結果分析

(1) 改進的地形校正模型的合理性驗證

式(1)描述的是斜面像元的亮度值與太陽入射角的線性關系，LT為因變量，cosi為自變量。對景號為153580、153595和153610的3景影像的4個波段分別應用式(1)進行最小二乘擬合，擬合出來的直線的相關系數見表2。

表2　根據式(1)擬合的三景影像各波段的相關系數

式(5)描述的是考慮反射角后斜面像元的亮度值與太陽入射角及像元的反射角的關系，LT為因變量，cosi·coset為自變量。對景號為153580、153595和153610的3景影像分別應用式(5)進行最小二乘擬合，擬合出來的直線的相關系數見表3。

表3　根據式(5)擬合的三景影像各波段的相關系數

對比表2和表3可以看出，在模型中加入反射角的影響后，3景影像擬合出的直線相關系數有大幅度的提升。這說明影像像元的亮度值不僅與太陽入射角有關，還與像元的反射角有著密切的聯系，反射角的大小對最終到達傳感器處的輻射量有著很大的影響。進一步表明加入反射角的改進的地形校正模型的合理性。

(2) 改進的地形校正模型的試驗結果評價

① 視覺評價

由于高分一號寬視場相機拍攝的影像尺寸比較大，因此在這里截取整景地形校正后影像中的一塊影像，對這一區域影像塊進行放大顯示。圖5(a)、(b)、(c)、(d)分別為校正前影像、余弦校正的影像、C校正的影像、改進方法校正的影像。

圖5　原影像及3種方法校正后的影像圖

從圖5可以看出，校正前的影像有明顯的地勢起伏，立體感很強烈，陰坡陽坡的亮度值差異很大；校正后，3種地形校正方法基本上都達到了地形校正的效果，地形起伏對像元亮度值的影響基本上都消除了；但從余弦校正后的影像上可以看到，在太陽入射角很低的地方出現了大量的亮區，即出現過校正的現象；而在C校正和改進的校正模型的校正后影像上則沒有過校正的情況。仔細對比圖5(c)和圖5(d)可以發現，本文的改進方法校正后影像的陰影地區的地表信息恢復得更好，陰陽坡像元的亮度值趨于一致的趨勢更加明顯。

② 影像統計信息評價

為了進一步對影像的校正質量進行定量的客觀的評價，本文列出了影像校正后最大值、最小值、均值及標準差等統計量(見表4)。地形校正的效果越好，陰陽坡的像元的亮度值越接近，陰陽坡趨于一致的趨勢更加明顯，影像的標準差就更小。

表4　原影像與3種方法校正后的影像統計量

從表4中可以看出，余弦校正的最大值和均值都明顯大于未校正前的影像，校正后的影像的標準差不降反升，這是由于出現了大量的過校正現象所導致的。C校正和本文方法校正后，均值與原影像相比都沒有太大的變化，而影像的標準差比原影像有所降低，說明這兩種方法都使陰陽坡的像元亮度值更加接近，達到了地形校正的目的。從標準差減小的幅度來看，本文方法校正后影像的標準差減小的幅度更大，說明本文方法對地形陰影的消除、陰陽坡亮度值的均一化效果更好。

③ 像元亮度值與光照系數的相關性評價

地形校正的主要目的是消除影像亮度值與太陽入射角的關系。相關系數的平方即為決定系數，決定系數的大小決定了相關的密切程度，決定系數說明模型中的自變量對因變量的影響程度。因此如果像元的亮度值不受太陽入射角的影響，則線性擬合方程的決定系數應該接近0，且越接近0說明校正的效果越好。表5列出了影像校正前和用3種方法校正后影像的各個波段的像元亮度值與光照系數的判定系數。

從表5中可以看出，余弦校正的決定系數比未校正之前要高出許多。這是因為余弦校正中出現了嚴重的過校正現象，在很多低光照區域出現了比原有像素亮度值高出許多的非正常像元導致的，已經出現了高度的負相關現象。C校正和本文方法校正后的決定系數比未校正之前都要小，說明兩種方法降低了相關的密切程度，模型中光照系數對影像亮度值的影響程度降低。比較C校正模型和本文方法，除了第一波段，本文方法在其他幾個波段的決定系數均小于C校正模型，而且第一波段的決定系數也非常接近0，說明整體上本文方法是優于C校正模型的。

表5原影像與3種方法校正后的影像的各波段擬合決定系數

band校正方法原影像余弦校正C校正改進方法校正band10.0154060.5543990.0002470.007099band20.0236270.3849950.0009340.000797band30.0243930.2279580.0014190.000017band40.0454170.2408060.0023950.000495

四、結論

1) 本文針對線陣推掃傳感器的成像特點，給出了像元平面反射角及斜面反射角的計算方法，對衛星星下點成像時和非星下點成像時的平面反射角的范圍進行了估算，利用多景影像驗證了反射角對到達傳感器處的輻射量有影響這一普遍性規律。

2) 本文模型是在C校正的模型上改進的，既保留了C校正的基本思想，又添加了反射角的影響，對非星下點成像的大場景寬幅遙感影像進行了地形校正，并從多個角度對校正結果進行了分析和評價。結果表明，本文的改進方法是合理、有效的，是優于傳統的C校正的。

參考文獻：

[1]劉三超, 張萬昌, 蔣建軍, 等. 用 TM 影像和 DEM 獲取黑河流域地表反射率和反照率[J]. 地理科學, 2003, 23(5): 585-591.

[2]TEILLET P M, GUINDON B, GOODENOUGH D G. On the Slope-aspect Correction of Multispectral Scanner Data[J]. Canadian Journal of Remote Sensing, 1982, 8(2): 84-106.

[3]MCDONALD E R, WU X, CACCETTA P A, et al.Illumination Correction of Landsat TM Data in South East NSW[M].[S.l.]: Environment Australia, 2002.

[4]JENSEN J R. 遙感數字影像處理導論[M]. 陳曉玲, 龔威,李平湘，等，譯.3版.北京：機械工業出版社, 2007：177-195.

[5]樊鵬山, 熊偉, 李智. 載荷側擺情況下衛星覆蓋區域計算方法研究[C]∥系統仿真技術及其應用學術會議論文集.宜昌：中國自動化學會系統仿真專業委員會，2009.

[6]GAO M L, ZHAO W J, GONG Z N et al.Topographic Correction of ZY-3 Satellite Images and Its Effects on Estimation of Shrub Leaf Biomass in Mountainous Areas[J]. Remote Sensing, 2014, 6(4):2745-2764.

[7]TEILLET P M,GUINDON B,GOODENONGH D G. On the Slope-aspect Correction of Multispectral Scanner Data[J].Photogrammetric Engineering and Remote Sensing, 1980, 9(46): 1183-1189.

引文格式： 臧熹，楊博，齊建偉,等. 一種改進的遙感影像地形校正方法[J].測繪通報，2015(1)：75-80.DOI:10.13474/j.cnki.11-2246.2015.0015

作者簡介：臧熹(1990—)，女，碩士生，主要從事遙感影像處理方面的研究。E-mail： zx901002@sina.com

基金項目：國家973計劃(2014CB744201)；國家863計劃(2011AA120203)；國家自然科學基金(41371430；91438203)

收稿日期：2014-07-30

中圖分類號：P237

文獻標識碼：B

文章編號：0494-0911(2015)01-0075-06