基于MODIS影像的亞象元積雪覆蓋率提取方法研究

積雪是地球表面最為活躍的自然要素之一,其特征(如積雪面積、雪深、積雪反照率等)是全球能量平衡模型中的重要輸入?yún)?shù)中分辨率成像光譜儀(MODIS)是EOS系統(tǒng)中主要的傳感器之一， 它具有較高的空間分辨率、光譜分辨率和時(shí)間分辨率、雖然MODIS全球每日500m空間分辨率雪產(chǎn)品較以往雪蓋制圖精度有了一定的提高 ,但是混合像元的問題仍比較突出.針對(duì)混合像元問題,一些學(xué)者提出了光譜線性分解模型通過選取恰當(dāng)?shù)牡乇砀采w類型組合,就能在一定程度上提高最終結(jié)果的精度[8].但是這種模型運(yùn)算量較大,不適合積雪蓋率的業(yè)務(wù)化制圖.Salomoson[9]等提出基于歸一化雪蓋指數(shù)(normalizeddifferencesnowindex,NDSI)的MODIS像元雪蓋率(FRA)反演簡(jiǎn)化方法、中國(guó)學(xué)者曹云剛[10]考慮地表覆蓋,建立了像元雪蓋率與雪蓋指數(shù)、植被指數(shù)之間的線性關(guān)系模型.然而精確反演雪蓋率依舊非常困難其中一個(gè)重要原因在于積雪本身的物理性質(zhì)也會(huì)對(duì)NDSI產(chǎn)生影響另一方面地表覆蓋的復(fù)雜性也會(huì)降低NDSI與積雪比例間的相關(guān)性.若對(duì)數(shù)據(jù)不加區(qū)分統(tǒng)一建模,則必然會(huì)對(duì)模型的精度產(chǎn)生影響.因此,本文嘗試將雪蓋率按NDSI分類,然后對(duì)不同類別分別進(jìn)行建模、選擇研究區(qū)域需要在一定程度上,研究區(qū)具有很好的代表性本文使用MODIS1B500m分辨率數(shù)據(jù)反演雪蓋率(1),1、2波段為MODIS250m原500m后得到.1同一天ETM+數(shù)據(jù)作為檢驗(yàn)數(shù)據(jù)(數(shù)據(jù)已經(jīng)過投影轉(zhuǎn)換UTM投影,30m).表1“l(fā)uhat/f：?+和MODIS影像數(shù)據(jù)信息校感器類型LandsaL.ETM+

成像時(shí)間2001/10?09

穫蓋范圍那，4197'38.4779N

便用波段沁in 2,4t5EOS'MDDK 2001/10/09數(shù)據(jù)預(yù)處理

100-4S13—103-3857E35-49715514N994468 4067E

500m 1,2,4?5.6數(shù)據(jù)預(yù)處理包括波段反射率計(jì)算以及投影轉(zhuǎn)換 TM由于ETM+數(shù)據(jù)已經(jīng)過投影轉(zhuǎn)換故只需計(jì)算其各波段的反射率值 .以ENVI413為例使用CalibratenUtilites中的LandsatTM功能計(jì)MODISMODIS1B數(shù)據(jù),MODIStoolsbow2tie工具對(duì)圖像去“蝴蝶結(jié)”效應(yīng)，接著使用MODIS1B影像自帶的控制點(diǎn)信息對(duì)圖像進(jìn)行幾何校正(1)計(jì)算各波段的反射率值:

.然后根據(jù)公式p=reflectanee-scalesx(DN2reflectanee-offsets),(1)式中,p射率;DN為像元值;refleetanee-scaleseflectanee-offsets為偏置系數(shù).得到波段反射率以后，從MODIS圖ETM+數(shù)據(jù)范圍的區(qū)域，并將其轉(zhuǎn)換為UTM投影.轉(zhuǎn)換過程中為保持波段反射率不變，故使用最鄰近插值法.為了與LandsatETM+數(shù)據(jù)進(jìn)行對(duì)比分析，還需要對(duì)裁剪后的MODIS影像進(jìn)行空間配準(zhǔn).MODIS1、4、3ETM4、3、2合成影像的對(duì)比分析選擇控制點(diǎn)MODIS影像進(jìn)行配準(zhǔn)得到相應(yīng)的MODIS圖像后，分別計(jì)算雪蓋指數(shù)(NDSI,(p4-p6)P(p4+p6))和植被指數(shù)(NDVI,(p2-p1)P(p2+p1)).雪蓋率計(jì)算ETM+數(shù)據(jù)，SNOWMAP方法［4］進(jìn)行分類.SNOWMAPNDSI,NDSI的閾值為0140NDSI0140,p4>0111時(shí)判定為雪像元,p4為L(zhǎng)andsatETM+4波段的反射率,1代表積雪像元,0代表非積雪像元，30m;然后3個(gè)研究區(qū)分MODISMODIS率MatlabMODISETM+積雪分類圖；MODIS像元對(duì)應(yīng)的區(qū)域中所包含的ETM+1的積雪像元個(gè)數(shù)MODIS像元雪蓋率(fraction).fraction=叢?"其中，nsMODISETM+二值積雪分類圖中積雪像元的個(gè)數(shù)；nMODISETM類圖的像元個(gè)數(shù).誤差分析(1)位于中間部分的NDSI值相對(duì)較稀疏，從而導(dǎo)致模型擬合結(jié)果受到影響；(2)從空間分布來(lái)看這些數(shù)據(jù)主要集中在高積雪區(qū)與非積雪區(qū)的過渡地帶而這些地區(qū)地形起伏較大所以山體的陰影也會(huì)對(duì)模型反演的精度產(chǎn)生影響;(3)在過渡地區(qū)地表覆蓋復(fù)雜僅僅使用NDVI并不能完全反映地表的影響;(4)MODIS數(shù)據(jù)和ETM+差異也會(huì)導(dǎo)致雪蓋率的計(jì)算產(chǎn)生誤差所以這些區(qū)域的雪蓋率值波動(dòng)較大。------------------------------MODIS?像元積雪覆蓋率提取方法1.pdf針對(duì)目前應(yīng)用廣泛的MODIS傳感器數(shù)據(jù),充分利用了雪蓋指數(shù)在積雪監(jiān)測(cè)中的重要性,并在考慮了地表覆蓋的情況下,建立了像元雪蓋率與雪蓋指數(shù)、植被指數(shù)之間的線性關(guān)系模型,并利用ETM+據(jù)對(duì)模型估算的雪蓋率進(jìn)行了驗(yàn)證.結(jié)果表明,該方法能有效地提取?像元尺度的積雪信息1 引言積雪是地表覆蓋的重要組成部分,積雪表面的高反射特性使其成為決定地球輻射平衡的一個(gè)關(guān)鍵因子,也使其成為區(qū)域與全球氣候變化研究的重要內(nèi)容之一.同時(shí),積雪的空間分布還是山區(qū)與季節(jié)性積雪區(qū)水文、氣象模型的重要輸入因子[1]因此,研究積雪分布及其特性,對(duì)于水資源利用、災(zāi)害分析、大氣環(huán)流分析以及氣候演變有著重要的意義[2近年來(lái),遙感資料已經(jīng)成為積雪監(jiān)測(cè)的主要信息來(lái)源,其中最主要的有AVHRR(AdvancedVeryHighResolutionRadiometer)、MODIS(ModerateResolutionImagingSpectroradiometer)、TM/ETM+(ThematicMapper/EnhancedThematicMapper)等.針對(duì)這些數(shù)據(jù)源,多種方法已經(jīng)用于積雪識(shí)別與雪蓋制圖,包括:目視判讀,多光譜圖像分類,反射率特性計(jì)算法,決策樹等[3-4].這些方法的一個(gè)共同特點(diǎn)就是,將遙感圖像分為兩個(gè)部分,即積雪區(qū)與非積雪區(qū).也就是說(shuō),在進(jìn)行積雪像元識(shí)別時(shí),都假設(shè)所有的像元是純像元,一個(gè)像元對(duì)應(yīng)一種地物類型.因此,當(dāng)遙感圖像的空間分辨率較低時(shí),混合像元的問題就比較突出,造成分類結(jié)果的精度不夠理想[5].針對(duì)混合像元問題,一些學(xué)者提出了光譜線性分解模型,通過選取恰當(dāng)?shù)牡乇砀采w類型組合,就能在一定程度上提高最終結(jié)果的精度[6]光譜線性分解模型在精度上有所提高,但是它通常使用最小二乘法來(lái)研究區(qū)域與數(shù)據(jù)準(zhǔn)2估算每個(gè)像元中的組分比例,當(dāng)遙感數(shù)據(jù)量較大時(shí),就會(huì)耗費(fèi)大量的計(jì)算時(shí)間,并且這種方法的自動(dòng)化程度也很低,很難滿足一些業(yè)務(wù)上的需要.針對(duì)這樣的問題,本文引入了兩個(gè)基本輸入?yún)?shù)（雪蓋指數(shù)與植被指數(shù)）,提出了一種簡(jiǎn)化的算法來(lái)計(jì)算單個(gè)像元中積雪所占的比例（即像元雪蓋率,snowfraction）.采用研究區(qū)內(nèi)的MODIS與ETM+數(shù)據(jù),對(duì)算法中的參數(shù)進(jìn)行了計(jì)算,最后對(duì)算法的結(jié)果進(jìn)行了比較和分析 研究區(qū)域與數(shù)據(jù)準(zhǔn)2域是高原主要的積雪區(qū)之一,且地形比較復(fù)雜,植被覆蓋類型多樣,4ETM+影像數(shù)據(jù)1（UTM投影30m）.同時(shí)選擇了包含以上區(qū)域兩個(gè)時(shí)期內(nèi)的MODIS1B500m（MODIS1、2、、61、2250m,4、6波段分辨率為500m,500m）.1相關(guān)ETM+數(shù)據(jù)信息Table1InformalionoftheETM+■dala黑域黑域成像日期軌道號(hào)太陽(yáng)方位角12<1)]1021pl36rO3&15],23M7太陽(yáng)高度角43.U44歹2(1)11021P136rO39150.225744.178^31023149.714(f44.70H2斗2(1)11U23pl34rO4O15],7295142.435ff2KTVI+各波段耳值Table2ValtHLofdifferentwavebantfeETM+波段號(hào)/<i】/(*1葉*匸]yfw*315514]OU5225.77K2,073 數(shù)據(jù)預(yù)處理數(shù)據(jù)預(yù)處理包括波段反射率計(jì)算以及投影轉(zhuǎn)換.對(duì)于本文所使用的ETM+數(shù)據(jù)，因其已經(jīng)過投影轉(zhuǎn)換，故只需要計(jì)算其各波段的反射率值，使用下面的公式將其像元灰度值轉(zhuǎn)換為反射率(I)￡(>xoosO式中:a、b為定標(biāo)常數(shù);L為輻射亮度；d為日地距離；E0為大氣頂層太陽(yáng)輻照度；！太陽(yáng)天頂角.MODIS1B數(shù)據(jù)，MODIStoolsbow-tie工具對(duì)圖像的邊緣重疊效應(yīng)進(jìn)行改正,接著使用MODIS1BGeoreferenee命令對(duì)圖像進(jìn)行了幾何校正（ENVI軟件為例）.然后，使用下面的公式將MODIS各波段的像元灰度值轉(zhuǎn)換為反射率值[7]:P=reflectan x（Z）A—refleclance_offsets）式中: 為反射率；DN為像元值各波段比例系數(shù)（refleetance_scales） 與偏置系數(shù)（refleetance_offsets）見表3（只列出了文中使用的幾個(gè)波段參數(shù)表3MODIS 據(jù)部分波段定標(biāo)參數(shù)值）值）.60+000()3238{)波段號(hào)比例系數(shù)傾直系數(shù)1(1(XW51370.00X)31970()4{)在得到波段反射率以后，MODISETM+數(shù)據(jù)范圍的區(qū)域UTM投影.MODIS圖像，使用下面的公式計(jì)算其雪蓋指數(shù)（NDSI）與植被指數(shù)（NDVIAOS/二NDVI=式中:PMODISi波段反射率.4 研究方法

(51本文以ETM+數(shù)據(jù)獲取的雪蓋率作為#真值？，(51MODISNDSI、NDVI之間的線性關(guān)系.，雪蓋率真值？計(jì)算需要經(jīng)過以下幾個(gè)步驟首先通過監(jiān)督分類的ETM+數(shù)據(jù)進(jìn)行積雪信息的提取，MODISETM+1的像元數(shù)（積雪像元）；使用式（6）來(lái)計(jì)算MODIS像元雪蓋率raction）:fraction= 積？;像元數(shù)像元總數(shù)

(6/2的散點(diǎn)圖可以看出雪蓋率存在于兩個(gè)主要的區(qū)間，即高值區(qū)（0.8）和低值區(qū)（0.1）,其中高值區(qū)表明了積雪是一種大面積比較連續(xù)的地表覆蓋,而低值區(qū)的存在除了研究區(qū)內(nèi)積雪較少的原因以外，還在一定程度上受計(jì)算誤差的影響（其中包括MODIS與ETM+數(shù)據(jù)之間的坐標(biāo)配準(zhǔn)誤差）.針對(duì)NDSI與雪蓋率之間的關(guān)系，Bartonetal.[8]Salomonsonetal.[9]（fraction=0.18+0.37 ?NDSI+0.26 ?（NDSI）2）和直線（fraction=0.06+1.21 ?NDSI）其中fraction=a+bDSI（A）NDSI=a+b?fraction（B）,2的3種算法，AB.2其在低值區(qū)的同樣A,其在高值B,則能得到造成以上誤差的影響NDSI其中積雪顆（NDSI值）;對(duì)于積雪的低值區(qū)，由于地表覆蓋的復(fù)雜性NDSI與積雪比例之間的相關(guān)性降低.通過以上的分析本文認(rèn)為只使用單因子來(lái)進(jìn)行雪蓋率的反演還存在一定的不足因此NDVI值來(lái)進(jìn)的方法來(lái)計(jì)算雪蓋率:fraction=</+/>xNDS1+rxADVI(7)12得到的積雪信息作為雪蓋率的#真值？，通過線性回歸的方法，本文對(duì)上式中的3個(gè)參數(shù)進(jìn)行了擬合計(jì)算為了不增加計(jì)算的復(fù)雜性最后的參數(shù)只是對(duì)上面的各區(qū)域擬合參數(shù)進(jìn)行了平均 (對(duì)Salomonson提出的算法，B作變形，NDSI示雪蓋率A的參數(shù)進(jìn)行平均X5 使用各區(qū)域算法平均后得到的最終結(jié)果Table5Th<^(ruitioiiKforiliffc^rrntalgoritLms最后擬合結(jié)果Baiton fradinti-0 24+it63xADS/+0.31+SalcmonAon

fraetinn=0.1K+0T

97xAI)Si帛文算法f(X+(K1x\-Q?A5 算法結(jié)果驗(yàn)證及其精度分析4MODIS日積雪圖與本文算法計(jì)算的雪蓋率圖其中，MODIS積雪圖使用像元值200來(lái)表示積雪像元，它只是在像元層次上給出了區(qū)域的積雪分布情況而使用雪蓋率的方法可以得到像元積雪所占的比例，它能清楚地表現(xiàn)積雪分布的強(qiáng)弱情況，尤其是在積雪邊緣區(qū)域雪蓋率更能反應(yīng)積雪覆蓋漸變的細(xì)節(jié)信息. 針對(duì)圖4所示的區(qū)域，使用本文算法計(jì)算了雪蓋率并對(duì)結(jié)果進(jìn)行了統(tǒng)計(jì)和分析，5從圖可以看出，研究區(qū)域內(nèi)雪蓋率估算誤差分布的主要區(qū)間集中在[-0.25,0.1]. 研究區(qū)積雪覆蓋率為0.598,由算法計(jì)算的像元雪蓋率統(tǒng)計(jì)得到的積雪覆蓋率為 0.509,14.8%.NDSI（6）NDSI0.6率.造成這種結(jié)果的原因可能(h)無(wú)駕積雪 0 0.5 1.04MODISH積雪圖(a)與雪蓋率圖(b)Pig.4