基于Landsat8 OLI影像的山區(qū)植被地形落影校正方法研究

吳勇鋒，江 洪

（福州大學(xué)數(shù)字中國研究院，空間數(shù)據(jù)挖掘與信息共享教育部重點(diǎn)實(shí)驗(yàn)室，衛(wèi)星空間信息技術(shù)綜合應(yīng)用國家地方聯(lián)合工程研究中心，福建 福州 350108）

地形陰影是地形影響的一種，它的存在嚴(yán)重干擾山區(qū)植被信息提取精度，給山區(qū)土地覆被解譯以及生態(tài)參量的遙感反演帶來巨大挑戰(zhàn)［1］.目前，地形陰影校正方法主要可以概括為3類：1）基于山地輻射傳輸模型方法同時(shí)考慮太陽直射輻射、大氣散射輻射以及周圍地形反射輻射，效果較好，如Proy校正模型［2］、Sandmeier模型［3］和光學(xué)遙感反射率計(jì)算模型［4］，但模型參數(shù)較多且計(jì)算過程復(fù)雜，推廣難度大；2）基于DEM數(shù)據(jù)的經(jīng)驗(yàn)校正方法以反射率與cosi的統(tǒng)計(jì)關(guān)系為基礎(chǔ)進(jìn)行經(jīng)驗(yàn)校正，包括統(tǒng)計(jì)-經(jīng)驗(yàn)?zāi)Ｐ停ㄈ鏣eil?let-回歸校正［5］、b校正［6］、VECA模型［7］）、歸一化模型（如二階校正模型［8］、坡度匹配模型［9］）、朗伯體反射率模型（如余弦校正模型、C校正模型［3］、SCS模型［10］、SCS+C模型［11］等）以及非朗伯體反射率模型（如Minnaert模型［12］），但該類方法只對本影有效，對落影校正效果不佳；3）基于波段組合優(yōu)化計(jì)算模型的方法通過構(gòu)建特殊的植被指數(shù)來獲取消除地形影響后的植被信息，如波段比模型［13］、地形調(diào)節(jié)植被指數(shù)（Topography Adjusted Vegetation Index，TAVI）［14-15］、歸一化差值山地植被指數(shù)（Normalized Difference Mountain Vegeta?tion Index，NDMVI）［16-17］、陰影消除植被指數(shù)（Shadow Eliminated Vegetation Index，SEVI）［18-21］、植被區(qū)分陰影消除植被指數(shù)（Vegetation Distinguished and Shadow Eliminated Vegetation Index，VDSEVI）［22］等，但該類模型只能獲取落影的指數(shù)信息，缺乏多光譜信息.多光譜信息是遙感影像分類的基礎(chǔ)，對山區(qū)土地利用分類至關(guān)重要.

為了獲取消除落影影響的多波段反射率信息，提升山區(qū)土地利用分類精度，筆者以福建省連江縣內(nèi)某山地丘陵區(qū)域的Landsat8 OLI影像為例，提出一種地形落影校正方法.首先采用地表反射率計(jì)算的SE?VI作為陰影校正信息源，通過構(gòu)建光照區(qū)SCS+C校正后紅綠藍(lán)波段反射率與SEVI間的隨機(jī)森林回歸模型，進(jìn)而利用落影區(qū)SEVI校正地形落影，從而獲取消除地形落影影響的紅綠藍(lán)波段地表反射率的光譜信息.

1 研究區(qū)與數(shù)據(jù)源

1.1研究區(qū)概況研究區(qū)位于福建省福州市東部（圖1），地理坐標(biāo)范圍介于119°25′46″E～119°30′55″E，26°7′21″N～26°12′10″N之間，總面積為72 km2，最小坡度為0°，最大坡度大于60°，平均坡度21.7°，標(biāo)準(zhǔn)方差10.37°.該地區(qū)以山地為主，地形崎嶇復(fù)雜，區(qū)內(nèi)植被以常綠闊葉林和竹林為主.

圖1 研究區(qū)地理位置

1.2數(shù)據(jù)源及預(yù)處理本研究采用Landsat8 OLI衛(wèi)星影像、ASTERGDEM V2高程數(shù)據(jù)2種數(shù)據(jù)源，如表1所示.其中Landsat8 OLI衛(wèi)星影像過境日期為2019年12月11日，行列號(hào)為119/042，太陽高度角36.98°，太陽方位角155.56°.

表1 數(shù)據(jù)源

為消除大氣和光照等因素對地物反射的影響，利用Landsat8 OLI衛(wèi)星影像數(shù)據(jù)元文件提供的輻射定標(biāo)參數(shù)將遙感影像的DN值轉(zhuǎn)換為輻射亮度值.采用ENVI5.3中的FLAASH大氣校正模塊對影像的輻射亮度進(jìn)行大氣校正，從而獲得更接近真實(shí)地物的地表反射率數(shù)據(jù)，并對DEM數(shù)據(jù)進(jìn)行投影變換，使其與Landsat8 OLI衛(wèi)星影像具有相同的投影坐標(biāo)系統(tǒng).

2 研究方法

為了實(shí)現(xiàn)山區(qū)遙感影像地形落影校正，首先利用隨機(jī)森林分類方法將研究區(qū)分為陰影區(qū)、光照區(qū)和水體，利用DEM數(shù)據(jù)和影像角度信息相結(jié)合提取本影，落影則由陰影區(qū)和本影的差運(yùn)算得到；其次，利用研究區(qū)地表反射率數(shù)據(jù)進(jìn)行SCS+C校正與SEVI計(jì)算，以光照區(qū)為掩膜對SCS+C校正結(jié)果和SEVI進(jìn)行隨機(jī)采樣，生成樣本集，進(jìn)而利用訓(xùn)練樣本集生成隨機(jī)森林回歸模型，并結(jié)合測試樣本集完成隨機(jī)森林回歸模型精度評(píng)價(jià).最后，結(jié)合落影區(qū)域的SEVI和隨機(jī)森林回歸模型，實(shí)現(xiàn)地形落影校正，從而獲取消除地形落影影響的地表反射率信息，之后對地形落影校正效果進(jìn)行評(píng)價(jià).研究的技術(shù)流程圖如圖2所示.

圖2 地形落影校正技術(shù)流程圖

2.1落影提取山區(qū)陰影依據(jù)其形成原因可分為本影和落影，如圖3所示.

圖3 地形陰影示意圖

本影可根據(jù)公式計(jì)算，落影由總體陰影和本影的差值求得.

其中，Sself表示陰影檢測結(jié)果，若為1則為本影，若為0則為非本影；σ表示地形坡度角，π表示180°，ω表示太陽方位角，β表示地形坡向角，γ表示太陽高度角.

2.2落影校正

2.2.1SCS+C校正SCS+C校正［11］計(jì)算公式

其中，LSCS+C表示SCS+C校正之后的像元值，LT表示SCS+C校正之前的像元值，θ表示太陽天頂角，σ表示地形坡度角，c表示地形校正參數(shù)，i表示太陽入射角.

cosi的計(jì)算公式

其中，i表示太陽入射角，σ表示地形坡度角，θ表示太陽天頂角，β表示地形坡向角，ω表示太陽方位角.

2.2.2SEVI計(jì)算SEVI［18］計(jì)算公式

其中，RVI表示比值植被指數(shù)，f(Δ)表示地形調(diào)節(jié)因子，SVI表示陰影植被指數(shù)，Bnir為近紅外波段反射率，Br為紅光波段反射率.

f(Δ)的計(jì)算步驟：

步驟1選擇具有明顯地形陰影效應(yīng)的樣區(qū)，確保陰陽坡陰影對等；

步驟2使f(Δ)從0開始，以0.001為間隔進(jìn)行循環(huán)疊加，同時(shí)計(jì)算SEVI與RVI的相關(guān)系數(shù)r1以及SEVI與SVI的相關(guān)系數(shù)r2；

步驟3當(dāng)r1和r2的差值無限接近時(shí)，得到最優(yōu)的f(Δ).

其中，r1為SEVI和RVI的相關(guān)系數(shù)，r2為SEVI和SVI的相關(guān)系數(shù)，n為參與f(Δ)計(jì)算的影像像元數(shù)，x為SE?VI，y1為RVI，y2為SVI.

2.2.3隨機(jī)森林回歸校正

1）在光照區(qū)隨機(jī)選取樣本點(diǎn)提取SEVI以及SCS+C校正后紅綠藍(lán)波段地表反射率數(shù)據(jù)，將SEVI數(shù)據(jù)作為自變量，將與其對應(yīng)的SCS+C校正后紅綠藍(lán)波段地表反射率數(shù)據(jù)作為因變量，制作原始樣本集.

2）使用Python3.7的Sklearn模塊將原始樣本集隨機(jī)分為訓(xùn)練樣本集和測試樣本集，其中訓(xùn)練樣本和測試樣本的比例為7∶3.在訓(xùn)練樣本集中通過Bootstrap重抽樣得到K個(gè)與訓(xùn)練樣本集相等的訓(xùn)練樣本并生成K棵回歸樹，組成隨機(jī)森林回歸模型［23］（見式9）.在Bootstrap重抽樣時(shí)，未被抽取的概率為當(dāng)N→∞時(shí)，表明每次有36.8%的數(shù)據(jù)未被抽取，這些數(shù)據(jù)被稱為袋外數(shù)據(jù)（Out of Bag，OOB）［24］，可用于估計(jì)預(yù)測誤差.

其中，x為輸入向量，θk為獨(dú)立同分布隨機(jī)變量，K為回歸樹的個(gè)數(shù)，N為訓(xùn)練樣本集樣本數(shù).

3）回歸樹生長過程中，每個(gè)分裂節(jié)點(diǎn)隨機(jī)抽取所有變量中的M個(gè)變量作為當(dāng)前節(jié)點(diǎn)的特征子集（一般選取總變量的1/3，本文中M為1）進(jìn)行分裂，并使每棵回歸樹得到最大限度生長，分裂過程中不需要剪枝.

4）將所有回歸樹的預(yù)測平均值作為最終預(yù)測結(jié)果［25］.其最終分類決策如下

其中，H(x)表示K棵回歸樹得到的預(yù)測值的平均值，hk(x)表示在訓(xùn)練集上訓(xùn)練出的第k棵回歸樹得到的估計(jì)值.

5）使用Python3.7的bayes_opt模塊的貝葉斯優(yōu)化（Bayesian Optimization）調(diào)參方法尋找使隨機(jī)森林回歸模型表現(xiàn)最優(yōu)的超參數(shù)，并用最優(yōu)超參數(shù)進(jìn)行模型訓(xùn)練，分別得到紅綠藍(lán)波段SCS+C校正后地表反射率與SEVI的隨機(jī)森林回歸模型，結(jié)合落影區(qū)域的SEVI數(shù)據(jù)即可預(yù)測出落影區(qū)紅綠藍(lán)波段的地表反射率數(shù)據(jù).

2.3落影校正效果評(píng)價(jià)模型評(píng)價(jià)是建模過程中的核心工作之一，通過計(jì)算隨機(jī)森林回歸模型的決定系數(shù)（公式（11））來評(píng)價(jià)模型的預(yù)測能力.

其中，yi為實(shí)際觀測值為模型預(yù)估值為樣本平均數(shù)，n為樣本數(shù).

紅綠藍(lán)波段隨機(jī)森林回歸模型的預(yù)測精度（決定系數(shù)）如表2所示.

表2 Landsat8 OLI影像各波段隨機(jī)森林回歸模型的預(yù)測精度

采用目視分析、統(tǒng)計(jì)特征分析、光譜特征分析、相對誤差分析和反射率與cosi相關(guān)性分析等方法進(jìn)行地形校正效果評(píng)價(jià).相對誤差絕對值（Absolute Relative Error，ARE）計(jì)算公式

其中，ARE表示相對誤差絕對值，Bshadow表示陰影區(qū)域內(nèi)光譜波段的平均值，Bsunlit表示相鄰非陰影陽坡光譜波段的平均值.

3 結(jié)果分析

3.1目視分析原始影像、SCS+C校正后影像地以及本文方法校正后真彩色影像如圖4.對比影像校正效果可知，原始影像中地形陰影效應(yīng)十分明顯，陰影區(qū)地表反射率偏低，光譜信息缺失，如圖4a所示.經(jīng)SCS+C校正后，地形陰影效應(yīng)削弱，本影區(qū)域光譜信息得以恢復(fù)，但仍有少部分落影區(qū)域光譜信息難以恢復(fù)，如圖4b所示.經(jīng)本文方法校正后，地形陰影效應(yīng)基本消除，落影區(qū)光譜信息進(jìn)一步恢復(fù)，光照區(qū)和陰影區(qū)的光譜差異減小，如圖4c所示.

圖4 地形校正前后效果對比

3.2統(tǒng)計(jì)特征分析為了定量分析地形校正方法的校正效果，分別統(tǒng)計(jì)地形校正前后影像均值、標(biāo)準(zhǔn)差以及變異系數(shù)（Coefficient of Variation，CV）［26］，并進(jìn)行對比分析.影像標(biāo)準(zhǔn)差反映的是影像各波段的離散程度，地形校正后，各波段影像像元間的差異縮小.若校正后影像標(biāo)準(zhǔn)差小于原始影像，表明該校正方法具有一定的校正效果.影像變異系數(shù)是指其標(biāo)準(zhǔn)差和平均值之比，地形校正后變異系數(shù)越小，表明地形校正效果越好.由地形校正前后研究區(qū)影像各波段均值、標(biāo)準(zhǔn)差及變異系數(shù)可知，SCS+C校正后，影像各波段標(biāo)準(zhǔn)差及變異系數(shù)相對于原始影像均有所下降，說明SCS+C校正都有一定的校正效果.本文方法校正后，影像各波段標(biāo)準(zhǔn)差及變異系數(shù)均小于SCS+C校正，說明本文方法彌補(bǔ)了SCS+C校正的不足，使影像像元間的差異進(jìn)一步縮小，如表3所示.

表3 地形校正前后影像統(tǒng)計(jì)參數(shù)對比

圖5 局部地形校正前后效果對比

3.3光譜特征分析依據(jù)陰影檢測結(jié)果，選取30組本影、落影及相鄰非陰影陽坡樣本，并分別統(tǒng)計(jì)30組樣本在原始影像、SCS+C校正以及本文方法校正后紅綠藍(lán)三波段的光譜均值.統(tǒng)計(jì)結(jié)果顯示，原始影像的紅綠藍(lán)波段光譜值在本影和落影處均低于相鄰非陰影陽坡.SCS+C校正后本影處的紅綠藍(lán)波段光譜值恢復(fù)到陽坡水平，而落影依舊低于陽坡，證明SCS+C校正對本影的校正效果良好，但對落影的校正效果欠佳.經(jīng)本文方法校正后，落影處的紅綠藍(lán)波段光譜值校正到陽坡水平，證明本文方法能夠?qū)β溆斑M(jìn)行校正.

3.4本影、落影相對誤差分析地形校正前后，紅綠藍(lán)波段對應(yīng)的本影、落影相對誤差絕對值統(tǒng)計(jì)結(jié)果表明，如表4所示：

表4 紅綠藍(lán)波段本影落影區(qū)域相對誤差絕對值

1）未經(jīng)地形校正的各光譜波段本影、落影與其相鄰陽坡的相對誤差較大，相對誤差絕對值在30%～50%左右.

2）SCS+C校正后，紅綠藍(lán)波段本影與其相鄰陽坡的相對誤差大幅降低，分別從40.40%、43.43%、29.28%降至8.75%、13.35%、6.28%，而落影與其相鄰陽坡的相對誤差降幅相對較小，分別從48.76%、51.30%、38.50%降至38.25%、41.23%、29.84%，表明SCS+C校正對本影校正效果較好，而對落影校正效果欠佳.

3）本文方法校正后，紅綠藍(lán)波段落影與其相鄰陽坡的相對誤差大幅降低，分別從48.76%、51.30%、38.50%降至0.43%、1.82%、1.91%，表明本文方法在落影區(qū)域能夠取得較好的校正效果.

圖6 30組樣本紅綠藍(lán)波段折線圖

3.5反射率與cosi相關(guān)性分析由于地形效應(yīng)影響，原始遙感影像地表反射率與cosi之間存在一定相關(guān)性.地形校正后，此相關(guān)性將被削弱甚至消除，故地形校正前后兩者決定系數(shù)的變化被廣泛用于檢驗(yàn)地形校正的效果［21］.地形校正前后各波段反射率與cosi的斜率及決定系數(shù)如表5所示.影像校正前，各波段像元值與其cosi都存在一定相關(guān)性.SCS+C校正后，研究區(qū)影像各波段斜率和決定系數(shù)均顯著降低，表明SCS+C校正方法能夠削弱地形陰影效應(yīng).本文方法校正后，斜率和決定系數(shù)的下降幅度都高于SCS+C校正，表明本文方法能夠彌補(bǔ)SCS+C校正的不足，使整體校正效果更優(yōu).

表5 地表反射率與光照系數(shù)線性相關(guān)性分析

4 結(jié)束語

結(jié)合SCS+C模型和SEVI指數(shù)完成山區(qū)植被地形落影可見光波段校正，并對本文方法及傳統(tǒng)地形校正方法進(jìn)行比較.未經(jīng)地形校正的原始影像地形陰影效應(yīng)明顯，本影和落影與相鄰陽坡的相對誤差在30%～50%左右.因此，利用遙感影像研究崎嶇山區(qū)植被時(shí)，要考慮地形陰影效應(yīng).SCS+C等傳統(tǒng)地形校正對本影的校正效果較好，但對落影校正效果欠佳.本文方法對山區(qū)落影區(qū)域植被可見光多光譜信息校正效果明顯，彌補(bǔ)了傳統(tǒng)地形校正方法對落影校正失效的缺陷和SEVI只能獲取消除地形陰影干擾的單一指數(shù)信息的不足.本文地形落影校正方法目前對可見光波段效果較好，對其他波段的適用性還有待進(jìn)一步研究.